#!/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
# ----------------------------------------------------------------------------

# ----------------------------------------------------------------------------
#
# LES run script for real case:
#   input from ICON to rerun with SCM closely point from global ICON
#
# Martin Koehler, Jul. 2021
# ----------------------------------------------------------------------------
#PBS -S /bin/bash
#PBS -q sx_norm
#PBS -v NE=16,CPE=8
#PBS -l elapstim_req=24:00:00
#PBS --venode=${NE}          # NE: number VE nodes, CPE: number cores on each VE
#PBS --venum-lhost=4         # Number of VE per logical host
#PBS -l cpunum_job=3         # number x86 cores per logical host (1ini+#IO)
#PBS -l coresz_prc=0
#PBS -T necmpi_hydra         # parallel job mit NEC-MPIs
#PBS --use-hca=2             # number of HCA per logical host (2)
##PBS -N SCM_ICONini
#PBS -o /hpc/uwork/${USER}/wq/iconLES.%s
#PBS -j o
# ----------------------------------------------------------------------------

date
module load mpi/3.5.0
set -x
ulimit -s unlimited
. /etc/profile

# OpenMP settings
export OMP_SCHEDULE="static"
export OMP_DYNAMIC="false"
export OMP_NUM_THREADS=1                   # check total number of requested cores!

# Run information
let PPN=${CPE}*${NE}/${OMP_NUM_THREADS}    # no. of MPI procs. per NQS job
let NE1=${NE}-1
mpi_atm_io_procs=2                         # >0 for atmosphere plus land (not working for monitoring)
echo "PPN etc. $NN $NE $CPE $PPN $mpi_atm_io_procs"

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

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

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

# base directory for ICON sources and binary:
ICONDIR=${PBS_O_WORKDIR%/*/*/*}    # local icon directory (up from scripts/scm/run)
#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

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

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


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

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

inidate='2020071200'               # initial date
start_date="2020-07-12T00:00:00Z"

location='lat52.2_lon14.1'         # Lindenberg: lat52.2_lon14.1
corio_lat=52.2

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

# the grid parameters
#atmo_dyn_grids="Torus_Triangles_100x100_100m"
atmo_dyn_grids="Torus_Triangles_50x50_200m"
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}.nc',"
done
radiation_grid_filename=""
for gridfile in ${atmo_rad_grids}; do
  radiation_grid_filename="${radiation_grid_filename} '${gridfile}',"
done

# ----------------------------------------------------------------------------
# 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}

# grid
ln -sf ${GRIDDIR}/${atmo_dyn_grids}.nc .

# initial condition and forcing
ln -sf ${SCMDATA}/init_data/init_SCM_data_ICON_${inidate}_${location}.nc init_SCM.nc

# SCM extpar file
ln -sf ${SCMDATA}/extpar/extpar_${atmo_dyn_grids}_${location}.nc extpar_${atmo_dyn_grids}.nc
# 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'

# ----------------------------------------------------------------------------
# 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          =  752
 p_test_run      = .false.
 l_test_openmp   = .false.
 l_log_checks    = .true.
 num_io_procs    =  ${mpi_atm_io_procs}    ! number of I/O processors
 num_restart_procs = 0
 num_prefetch_proc = 0
 iorder_sendrecv =  3
 proc0_shift     = 1
 use_omp_input   = .true.
/

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

&scm_nml
 i_scm_netcdf  = 1                    ! read initial profiles and forcings from netcdf
 lscm_read_tke = .FALSE.              ! read initial tke from netcdf
 lscm_read_z0  = .TRUE.               ! read initial z0 from netcdf
 lscm_icon_ini = .TRUE.               ! read initial conditions from ICON output
 scm_sfc_mom   = 0                    ! 0: TURBTRANS: no prescribed u*
 scm_sfc_temp  = 0                    ! 0: TERRA:     no prescribed sensible heat flux at surface
 scm_sfc_qv    = 0                    ! 0: TERRA:     no prescribed latent heat flux at surface
 lscm_random_noise = .TRUE.           ! add random noise at init
/

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

&run_nml
 num_lev        = 90           ! number of full levels of vertical grid
 dtime          = ${dtime}     ! timestep in seconds
 nsteps         = ${nsteps}
 ldynamics      = .TRUE.       ! compute adiabatic dynamic tendencies
 ltransport     = .TRUE.
 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'
 check_uuid_gracefully = .TRUE.
/

&nwp_phy_nml
 inwp_gscp       = 1
 inwp_convection = 0 ! 1:Tiedtke/Bechtold, 0:off
 inwp_radiation  = 4 ! 1:RRTM radiation, 4:ecRad
 inwp_cldcover   = 5 ! 3: clouds from COSMO SGS cloud scheme 0: no cloud 5: grid-scale clouds
 inwp_turb       = 5 ! 1: TKE diffusion and transfer, 5: Smagorinsky
 inwp_satad      = 1
 inwp_sso        = 0
 inwp_gwd        = 0
 inwp_surface    = 1 ! 0: none; 1: TERRA   (0: simple ocean, sea-ice albedo!)
 icapdcycl       = 3 ! apply CAPE modification to improve diurnalcycle over tropical land (optimizes NWP scores)
 latm_above_top  = .TRUE.  ! needed for radiation routine
 itype_z0        = 2
 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
/

&turbdiff_nml
 tkhmin                       = 0.6
 tkhmin_strat                 = 1.0
 tkmmin                       = 0.75
 tkmmin_strat                 = 4
 alpha0                       = 0.0123
 alpha0_max                   = 0.0335
 alpha1                       = 0.125
 pat_len                      = 750.
 c_diff                       = 0.2
 rlam_heat                    = 10.0
 rat_sea                      = 0.8
 ltkesso                      = .true.  ! SSO dissipation energy used in TKE equation
 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
 itype_sher                   = 3       ! use horizontal shear production terms with 1/SQRT(Ri) scaling to prevent unwanted side effects
 ltkeshs                      = .true.
 a_hshr                       = 2.0
 icldm_turb                   = 1       ! 2: Gauss clouds for turbulence    1: grid scale clouds
 icldm_tran                   = 2       ! 2: Gauss clouds for surface layer 1: grid scale clouds
 lconst_z0                    = .TRUE.
 const_z0                     = 0.035
/

&lnd_nml
  ntiles         = 1
  nlev_snow      = 3
  lmulti_snow    = .false.
  itype_heatcond = 3
  idiag_snowfrac = 20
  lsnowtile      = .true.
  lseaice        = .true.
  llake          = .true.
  itype_lndtbl   = 4
  itype_evsl     = 4
  itype_trvg     = 3
  itype_root     = 2
  cwimax_ml      = 5.e-4
  c_soil         = 1.25
  c_soil_urb     = 0.5
  sstice_mode    = 2
  lprog_albsi    = .true.
  itype_snowevap = 2
/

&radiation_nml
 irad_o3                      = 9
 irad_aero                    = 6
 izenith                      = 4           ! 4: NWP default, 3: no annual cycle
 albedo_type                  = 2 ! Modis albedo
 vmr_co2                      = 390.e-06 ! values representative for 2012
 vmr_ch4                      = 1800.e-09
 vmr_n2o                      = 322.0e-09
 vmr_o2                       = 0.20946
 vmr_cfc11                    = 240.e-12
 vmr_cfc12                    = 532.e-12
 direct_albedo                = 4
 direct_albedo_water          = 3
 albedo_whitecap              = 1
 ecrad_llw_cloud_scat         = .true.
 ecRad_data_path              = '${ecRad_data_path}'
/

&ls_forcing_nml
 is_subsidence_moment         = .FALSE.
 is_subsidence_heat           = .FALSE.
 is_advection                 = .FALSE.
 is_advection_uv              = .FALSE.
 is_advection_tq              = .FALSE.
 is_geowind                   = .FALSE.
 is_rad_forcing               = .FALSE.
 is_nudging                   = .TRUE.
 is_nudging_uv                = .TRUE.
 is_nudging_t                 = .TRUE.         ! TRUE best T profiles, FALSE best T2m ???
 is_nudging_q                 = .TRUE.
 dt_relax_uv                  = 7200.0         ! 10800.0
 dt_relax_t                   = 7200.0         ! 10800.0
 dt_relax_q                   = 7200.0         ! 10800.0
/

&nonhydrostatic_nml
 iadv_rhotheta                = 2
 ivctype                      = 2
 itime_scheme                 = 4
 exner_expol                  = 0.333
 vwind_offctr                 = 0.2         ! 0.2 for R2B6 and higher resolution, 0.3 for lower resolution
 damp_height                  = 20000.
 rayleigh_coeff               = 0.0
 divdamp_order                = 24          ! 2 ass, 24 fc
 divdamp_type                 = 32          ! optional: 2 assimilation cycle, 32 forecast
 divdamp_fac                  = 0.004       ! 0.004 for R2B6; recommendation for R3B7: 0.003
 divdamp_trans_start          = 12500
 divdamp_trans_end            = 17500
 igradp_method                = 3
 l_zdiffu_t                   = .true.
 thslp_zdiffu                 = 0.02
 thhgtd_zdiffu                = 125.
 htop_moist_proc              = 22500.
 hbot_qvsubstep               = 16000.
/

&sleve_nml
 min_lay_thckn                = 20.         ! lowest level thickness (between half-levels)
 max_lay_thckn                = 400.        ! maximum layer thickness below htop_thcknlimit
 htop_thcknlimit              = 14000.
 top_height                   = 75000.
 stretch_fac                  = 0.9
 decay_scale_1                = 4000.
 decay_scale_2                = 2500.
 decay_exp                    = 1.2
 flat_height                  = 16000.
/

&dynamics_nml
 divavg_cntrwgt               = 0.50
 lcoriolis                    = .false.   ! SCM attention: .TRUE. distroys U-profile!!
/

&transport_nml
 ivadv_tracer                 = 3,3,3,3,3
 itype_hlimit                 = 3,4,4,4,4,0
 ihadv_tracer                 = 52,2,2,2,2,0
/

&diffusion_nml
 hdiff_order                  = 5
 itype_vn_diffu               = 1
 itype_t_diffu                = 2
 hdiff_efdt_ratio             = 24.0   ! for R2B6; recommendation for R3B7: 30.0
 hdiff_smag_fac               = 0.025  ! for R2B6; recommendation for R3B7: 0.02
 lhdiff_vn                    = .true.
 lhdiff_temp                  = .true.
/

&les_nml
 isrfc_type                   = 1      ! TERRA coupling
/

&extpar_nml
 itopo                        = 1   ! 0: analytical topo; 1: topography/ext. data read from file - if TERRA is on!!!
!n_iter_smooth_topo           = 1,1
!hgtdiff_max_smooth_topo      = 750.,750.,
!heightdiff_threshold         = 3000.
/

&output_nml
 output_time_unit =  1                        ! 1: seconds
 output_bounds    =  0., 10000000., 10800.    ! start, end, increment
 mode             =  1                        ! 1: forecast
 steps_per_file   = 10000
 include_last     = .TRUE.
 output_filename  = 'les_out_ML'
 filename_format  = "<output_filename>_<levtype>_<datetime2>"
 ml_varlist       = 'u','v','theta_v','tot_qv_dia','tot_qc_dia','tot_qi_dia','clc'
 output_grid      = .TRUE.
/

&output_nml
 output_time_unit =  1                        ! 1: seconds
 output_bounds    =  0., 10000000., 900.      ! start, end, increment
 mode             =  1                        ! 1: forecast
 steps_per_file   = 10000
 include_last     = .TRUE.
 output_filename  = 'les_out_SFC'
 filename_format  = "<output_filename>_<levtype>_<datetime2>"
 ml_varlist       = 'shfl_s', 'lhfl_s','u_10m', 'v_10m', 'qv_2m','t_2m', 'clct',
                    'tqv_dia', 'tqc_dia', 'tqi_dia','tot_prec','sp_10m'
 output_grid      = .TRUE.
/

EOF

  # 'group:all'
  # 'z_ifc','z_mc','u','v','w','temp','pres','rho','theta_v','pres_sfc','div',
  # 'qv','qc','qi','qs','qr','rh',
  # 'ashfl_s', 'alhfl_s', 'athb_s', 'athb_t', 'asob_s', 'asob_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','z_mc','lhfl_s','shfl_s','umfl_s','vmfl_s','tcm','tch','clc','tke','rcld','qhfl_s',
  # 'sob_s', 'thb_s','sob_t'


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

export VE_ERRCTL_ALLOCATE=MSG
export NMPI_PROGINF=YES
export VE_TRACEBACK=VERBOSE
export NMPI_SEPSELECT=3
export GMON_OUT_PREFIX=scal_prof
export VE_FPE_ENABLE=DIV,FOF,INV
export GFORTRAN_UNBUFFERED_PRECONNECTED=y
export NMPI_EXPORT="GFORTRAN_UNBUFFERED_PRECONNECTED"

date1=`date`

# good: NE=1, 4, 8 or 16, CPE=8, venum-lhost=8, cpunum_job=2, use-hca=2

/opt/nec/ve/mpi/3.5.0/bin/runtime/mpirun                                    \
  -v -vh     -node 0        -np 1      -env OMP_NUM_THREADS 1 ${MODEL_VH} : \
  -x -venode -node 0-${NE1} -np ${PPN} -env OMP_NUM_THREADS 1 ${MODEL_VE} : \
     -vh     -node 0        -np 2      -env OMP_NUM_THREADS 1 ${MODEL_VH}

# --- other mpirun options:

# debug output: each processor with separate output in EXPDIR
# add /opt/nec/ve/bin/mpisep.sh before ${MODEL_V*}

# BAD: NE=4, CPE=8, venum-lhost=2, cpunum_job=2, use-hca=2
#/opt/nec/ve/bin/mpirun -v    -vh     -node 0        -np 1      -env OMP_NUM_THREADS 1 ${MODEL_VH} : \
#                          -x -venode -node 0-${NE1} -np ${PPN} -env OMP_NUM_THREADS 1 ${MODEL_VE} : \
#                             -vh     -node 0        -np 4      -env OMP_NUM_THREADS 1 ${MODEL_VH}

# NE=8
#mpirun -v -vh     -node 0        -np 1                   -env OMP_NUM_THREADS 1 /opt/nec/ve/bin/mpisep.sh ${MODEL_VH} : \
#       -x -venode -node 0-${NE1} -np ${PPN}              -env OMP_NUM_THREADS 1 /opt/nec/ve/bin/mpisep.sh ${MODEL_VE} : \
#          -vh     -node 0        -np ${mpi_atm_io_procs} -env OMP_NUM_THREADS 1 /opt/nec/ve/bin/mpisep.sh ${MODEL_VH}
# : \
#          -vh     -node 1        -np 3      ${MODEL_VH}

# NE=16, CPE=8, venum-lhost=8, cpunum_job=2, use-hca=2

#/opt/nec/ve/bin/mpirun -v    -vh     -node 0        -np 1      -env OMP_NUM_THREADS 1  ${MODEL_VH} : \
#                         -x -venode  -node 0-${NE1} -np ${PPN} -env OMP_NUM_THREADS 1  ${MODEL_VE} : \
#                             -vh     -node 0        -np 2      -env OMP_NUM_THREADS 1  ${MODEL_VH}

# NE=32
#mpiexec -v    -vh     -node 0        -np 1         ${MODEL_SCAL} : \
#        -x    -venode -node 0-${NE1} -np ${PPN}    ${MODEL}      : \
#        -v    -vh     -nn 4          -nnp 3        ${MODEL_SCAL}

#mpiexec -v    -node 0   -np 1        -vh          ${MODEL_SCAL} : \
#        -x    -venode   -node 0-${NE1} -np ${PPN} ${MODEL}      : \
#        -v    -nn 4     -nnp 3       -vh          ${MODEL_SCAL}


echo '---- time before/after: '$data1'   '`date`


# SCM mean over columns:

outfile=les_out_ML_`echo ${inidate} |cut -c 1-8`*Z.nc
meanfile=`echo ${outfile}| sed 's/.nc/_mean.nc/'`
cdodir=/hpc/sw/cdo/2.5.3/x86/gnu/bin/
${cdodir}cdo fldavg ${outfile} ${meanfile}

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