#!/bin/bash

# ICON
#
# ------------------------------------------
# Copyright (C) 2004-2026, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
# Contact information: icon-model.org
# See AUTHORS.TXT for a list of authors
# See LICENSES/ for license information
# SPDX-License-Identifier: BSD-3-Clause
# ------------------------------------------

# ----------------------------------------------------------------------------
#
# SCM run script for ARM DEPHY case:
#
# Details: - Input data is only defined to 20km.  Therefore use ICON-D2 vertical grid
#
# ----------------------------------------------------------------------------
#PBS -S /bin/bash
#PBS -q sx_norm
#PBS -v NE=1,CPE=1
#PBS -l elapstim_req=01:00:00
#PBS --venode=${NE}          # NE: number VE nodes, CPE: number cores on each VE
#PBS --venum-lhost=1         # Number of VE per logical host
#PBS --cpunum-lhost=2        # Number of VH-x86-Cores per log. host
#PBS -T necmpi               # necessary for parallel jobs  with NEC-MPI
#PBS -o /hpc/uwork/${USER}/wq/iconSCM.%s
#PBS -j o
# ----------------------------------------------------------------------------

date
module load mpi
set -x

# OpenMP settings
export OMP_SCHEDULE="static"
export OMP_DYNAMIC="false"
export OMP_NUM_THREADS=1                   # check total number of requested cores!
export VE_TRACEBACK=VERBOSE
export VE_FPE_ENABLE=DIV,FOF,INV

# Run information
let PPN=${CPE}*${NE}/${OMP_NUM_THREADS}    # no. of MPI procs. per NQS job

ID=`echo $PBS_JOBID | cut -d: -f2 | cut -d. -f1`
ulimit -s unlimited

# ----------------------------------------------------------------------------
# specifiy experiment (idealized simulation)
# ----------------------------------------------------------------------------
EXPNAME=SCM_ARM_dephy_v1

# ----------------------------------------------------------------------------
# path definitions
# ----------------------------------------------------------------------------

# base directory for ICON sources and binary:
#ICONDIR=${PBS_O_WORKDIR%/*}       # local icon directory (up from SCM-ideal)
ICONDIR="/hpc/uhome/mkoehler/icon/icon-nwp-test2"

# SCM data directory (grids, init data, extpar)
SCMDATA=/hpc/uwork/mkoehler/scm/data       # at DWD on NEC

# directory with input grids:
GRIDDIR=${SCMDATA}/grid

# directory with init files:
INITDIR=${SCMDATA}/init_data/dephy_v1

# absolute path to output directory for results:
EXPDIR=${WORK}/run-icon/scm/${EXPNAME}

# path to model binary, including the executable:
MODEL_VE=${ICONDIR}/build/VE/bin/icon

# ----------------------------------------------------------------------------
# copy input data: grids, external parameters, model
# ----------------------------------------------------------------------------

# the directory for the experiment will be created, if not already there
if [ ! -d $EXPDIR ]; then
    mkdir -p $EXPDIR
fi
cd ${EXPDIR}

# files needed for radiation
ln -sf ${ICONDIR}/data/ECHAM6_CldOptProps.nc .
ln -sf ${ICONDIR}/data/rrtmg_lw.nc .
ln -sf ${ICONDIR}/data/rrtmg_sw.nc .

ecRad_data_path=${ICONDIR}'/externals/ecrad/data'


# ----------------------------------------------------------------------------
# model timing
# ----------------------------------------------------------------------------

# 0.4 s for 70 m res., 0.5 for 100 m; 30 s for 2.5 km?
dtime=60
dt_checkpoint=`expr 36 \* 86400`  # write restart file every hours (when lrestart = TRUE)
nhours=14.5
#nhours=1
#nsteps=`expr ${nhours} \* 3600 / ${dtime}`
nsteps=`echo ${nhours} \* 3600 / ${dtime} | bc`

start_date="1997-06-21T00:00:00Z"
end_date="1997-06-21T14:30:00Z"

# ----------------------------------------------------------------------------
# output
# ----------------------------------------------------------------------------
DT_DATA=`expr 1 \* 3600`      # output each n hours
DT_DATA=${dtime}              # output every time step
#n_in_ofile=60                # number of time steps per output file
n_in_ofile=10000              # number of time steps per output file

# ----------------------------------------------------------------------------
# grid namelist settings
# ----------------------------------------------------------------------------

# the grid parameters
atmo_dyn_grids="Torus_Triangles_4x4_2500m.nc"
atmo_rad_grids=""

# reconstruct the grid parameters in namelist form
dynamics_grid_filename=""
for gridfile in ${atmo_dyn_grids}; do
  dynamics_grid_filename="${dynamics_grid_filename} '${gridfile}',"
done
radiation_grid_filename=""
for gridfile in ${atmo_rad_grids}; do
  radiation_grid_filename="${radiation_grid_filename} '${gridfile}',"
done

ln -sf ${GRIDDIR}/${atmo_dyn_grids} .

# initial condition and forcing
ln -sf ${INITDIR}/ARMCU_REF_SCM_driver.nc init_SCM.nc

# ----------------------------------------------------------------------------
# create ICON master namelist
# ----------------------------------------------------------------------------

cat > icon_master.namelist << EOF

&master_nml
 lrestart                    =                     .FALSE.        ! .TRUE.=current experiment is resumed
/
&master_model_nml
 model_type                  =                          1         ! identifies which component to run (atmosphere,ocean,...)
 model_name                  =                      "ATMO"        ! character string for naming this component.
 model_namelist_filename     =       "NAMELIST_${EXPNAME}"        ! file name containing the model namelists
 model_min_rank              =                          1         ! start MPI rank for this model
 model_max_rank              =                      65536         ! end MPI rank for this model
 model_inc_rank              =                          1         ! stride of MPI ranks
/

! time_nml: specification of date and time------------------------------------
&time_nml
 ini_datetime_string         =               "$start_date"        ! initial date and time of the simulation
 end_datetime_string         =                 "$end_date"        ! initial date and time of the simulation
/

EOF

# ----------------------------------------------------------------------------
# model namelists
# ----------------------------------------------------------------------------
# For a complete list see Namelist_overview.pdf

cat > NAMELIST_${EXPNAME} << EOF

&parallel_nml
 nproma         =  32
 p_test_run     = .false.
 num_io_procs   =  1         ! number of I/O processors
 l_test_openmp  = .false.
 l_log_checks   = .false.
/

&grid_nml
 dynamics_grid_filename = "${atmo_dyn_grids}",
 !corio_lat             = 36.00
 corio_lat              = 35.00
 is_plane_torus         = .TRUE.
 l_scm_mode             = .TRUE.                    ! main logical to turn on SCM mode
/                                                   ! note: if ldynamics off this transitions to LES or CRM

&SCM_nml
 i_scm_netcdf = 2           ! read initial profiles and forcings from netcdf - universal format
 lscm_read_tke = .FALSE.    ! read initial tke from netcdf
 lscm_read_z0  = .FALSE.    ! read initial z0 from netcdf
 scm_sfc_mom   = 0          ! 0: TURBTRANS: no prescribed u*
 scm_sfc_temp  = 1          ! 0: TERRA:     no prescribed sensible heat flux at surface
 scm_sfc_qv    = 3          ! 0: TERRA:     no prescribed latent heat flux at surface
/                           ! 3: qv,sat(Tg)

&io_nml
 dt_checkpoint  = ${dt_checkpoint}
 lkeep_in_sync  = .true.
 lflux_avg      = .FALSE.      ! false: accumulated fluxes
/

&run_nml
 num_lev        = 65           ! number of full levels of vertical grid
 dtime          = ${dtime}     ! timestep in seconds
 nsteps         = ${nsteps}
 ldynamics      = .FALSE.      ! compute adiabatic dynamic tendencies
 ltransport     = .FALSE.
 ntracer        = 5            ! default: 0
 iforcing       = 3            ! 3: NWP forcing; 6:inhecham forcing
 ltestcase      = .TRUE.       ! run testcase
 ltimer         = .FALSE.      !
 msg_level      = 12           ! detailed report during integration
 output         = 'nml','totint'
/

&nwp_phy_nml
 inwp_gscp       = 1
 inwp_convection = 1 ! 1:Tiedtke/Bechtold, 0:off
 inwp_radiation  = 0 ! 1:RRTM radiation
 inwp_cldcover   = 1 ! 3: clouds from COSMO SGS cloud scheme 0:off
 inwp_turb       = 1 ! 1: TKE diffusion and transfer
 inwp_satad      = 1
 inwp_sso        = 0
 inwp_gwd        = 0
 inwp_surface    = 0 ! 0: none; 1: TERRA
 icapdcycl	 = 3 ! apply CAPE modification to improve diurnalcycle over tropical land (optimizes NWP scores)
 latm_above_top  = .TRUE.  ! needed for radiation routine
 itype_z0        = 1
 dt_rad	         = 1800.        ! Default: 1800   ! M. Koehler: 1440
 dt_conv         = 600.         ! Default: 600    ! M. Koehler: 360
 dt_sso	         = 600.         ! Default: 1200   ! M. Koehler: 720
 dt_gwd	         = 600.         ! Default: 1200   ! M. Koehler: 720
/

&radiation_nml
 irad_o3         = 9       ! ozone climatology (ecrad: 0, 7 GEMS, 9 MACC, 79, 97)
 ecRad_data_path = '${ICONDIR}/externals/ecrad/data'
/

&ls_forcing_nml
 is_subsidence_moment = .FALSE.
 is_subsidence_heat   = .FALSE.
 is_advection         = .TRUE.
 is_advection_uv      = .FALSE.
 is_advection_tq      = .TRUE.
 is_geowind           = .TRUE.
 is_rad_forcing       = .FALSE.
 is_nudging           = .FALSE.
 is_nudging_uv        = .FALSE.
 is_nudging_t         = .FALSE.
 is_nudging_q         = .FALSE.
 is_sim_rad           = .FALSE.
/

&turbdiff_nml
 tkhmin        = 0.6
 tkmmin        = 0.75
 pat_len       = 750.
 c_diff        = 0.2
 rat_sea       = 0.8
 ltkesso       = .TRUE.
!frcsmot       = 0.2    ! these 2 switches together apply vertical smoothing of the TKE source terms
!imode_frcsmot = 2      ! in the tropics (only), which reduces the moist bias in the tropical lower troposphere
 ! use horizontal shear production terms with 1/SQRT(Ri) scaling to prevent unwanted side effects:
 itype_sher    = 0
 ltkeshs       = .false.
 a_hshr        = 2.0
 !lconst_z0    = .TRUE.
 !const_z0     = 0.035
 icldm_tran    = 2
 icldm_turb    = 2
 rlam_heat     = 10.0
 rat_sea       = 0.8
 alpha1        = 0.125
/

&nwp_tuning_nml
!tune_entrorg    = 7.4E-3 ! is 1.85e-3 normally
lcalib_clcov     = .FALSE. !no tuning of cloud cover!
tune_rdepths     = 60000.
!tune_texc       = 0.2 !0.125
!tune_qexc       = 1.e-4 !1.25e-2
tune_rcucov      = 0.075
tune_rhebc_land  = 0.825
tune_zvz0i       = 0.85
tune_gust_factor = 7.0
icpl_turb_clc    = 2
tune_box_liq     = 0.04
tune_box_liq_asy = 4.0
/

&diffusion_nml
 lhdiff_temp  = .TRUE.
 lhdiff_vn    = .TRUE.
 lhdiff_w     = .TRUE.
 hdiff_order  = 4
/

&nonhydrostatic_nml
  damp_height     = 12250.
  divdamp_fac     = 0.004   ! 0.004 for R19B07
 !divdamp_order   =  ! 2 ass, 24 fc - may become unnecessary with incremental analysis update
  divdamp_order   = 24
 !divdamp_type    =   ! optional: 2 assimilation cycle, 3 forecast
  divdamp_type    = 32
  exner_expol     = 0.333
  hbot_qvsubstep  = 22500 ! r3b7: 19000.  else: 22500.
  htop_moist_proc = 22500.
  iadv_rhotheta   = 2
  igradp_method   = 3
  itime_scheme    = 4
  ivctype         = 2
  l_open_ubc      = .false.
  l_zdiffu_t      = .true.
  rayleigh_coeff  = 5.0
  thhgtd_zdiffu   = 125.
  thslp_zdiffu    = 0.02
  vwind_offctr    = 0.2
/

&sleve_nml
 decay_exp       = 1.2
 decay_scale_1   = 4000.
 decay_scale_2   = 2500.
 flat_height     = 16000.
 itype_laydistr  = 1
 min_lay_thckn   = 20.
 stretch_fac     = 0.65
 top_height      = 22000.
/

&extpar_nml
 itopo          = 0         ! 0: analytical topo; 1: topography/ext. data read from file
/

&dynamics_nml
 iequations     = 3         ! equation system
 idiv_method    = 1
 divavg_cntrwgt = 0.50
 lcoriolis      = .TRUE.
/

&output_nml
 output_time_unit =  1                        ! 1: seconds
 output_bounds    =  0., 10000000., 3600.     ! start, end, increment
 mode             =  1                        ! 1: forecast
 steps_per_file   = ${n_in_ofile}
 include_last     = .TRUE.
 output_filename  = 'out_${EXPNAME}'
 filename_format  = "<output_filename>_<levtype>_<datetime2>"
 ml_varlist       = 'z_ifc','z_mc','u','v','w','temp','pres','rho','theta_v','pres_sfc','div',
                    'qv','qc','qi','qs','qr','rh','exner','tot_qv_dia','tot_qc_dia','tot_qi_dia',
                    'tqv_dia','tqc_dia','tqi_dia','clc',
                    'accshfl_s', 'acclhfl_s', 'accthb_s', 'accthb_t', 'accsob_s', 'accsob_t',
                    'ddt_temp_radsw', 'ddt_temp_radlw', 'ddt_temp_turb', 'ddt_temp_drag', 'ddt_temp_pconv',
                    'ddt_qv_turb','ddt_qc_turb','ddt_qv_conv','ddt_qc_conv','u_10m', 'v_10m', 't_2m', 't_g',
                    'qv_s','lhfl_s','shfl_s','umfl_s','vmfl_s','tcm','tch','tke','rcld','qhfl_s'
 output_grid      = .TRUE.
/

EOF


# --- 80m resolution to 7000m
# &run_nml
#  num_lev        = 87        ! number of full levels of vertical grid
# &nonhydrostatic_nml
#  ivctype        = 2         ! sleve vertical coordinate
#  damp_height    = 5000.     ! top_height-damp_height should be about 15km
#  l_open_ubc     = .false.   ! top open upper boundary condition. might help to go higher
#  rayleigh_coeff = 0.0
#  htop_moist_proc= 22500.
#  hbot_qvsubstep	= 22500.    ! Default: 22500  ! M. Koehler: 19000; at least as large as htop_moist_proc
# /
#
# &sleve_nml
#  top_height    = 7000.
#  min_lay_thckn = -40.       ! thickness of lowest model layer
# /


# ----------------------------------------------------------------------------
# run the model!
# ----------------------------------------------------------------------------

# VE with openMP only (PPN = 1!)
date1=`date`

mpiexec -v -x -venode -np ${PPN} ${MODEL_VE}

echo $date1 ; date


# SCM mean over columns:

outfile=out_${EXPNAME}_ML_`echo ${inidate} |cut -c 1-8`*Z.nc
meanfile=`echo ${outfile}| sed 's/.nc/_mean.nc/'`
cdo fldavg ${outfile} ${meanfile}
