#!/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 GOAmazon case on NEC
#
# ----------------------------------------------------------------------------
#PBS -q sx_norm
#PBS -v NE=8,CPE=8,NMPI_ROOT=/opt/nec/ve/mpi/2.3.1
#PBS -v NMPI_MALLOC_HEAP_EXPANSION_SIZE=512
#PBS -v NMPI_MALLOC_MMAP_THRESHOLD=512
#PBS -v MPI_IB_VBUF_TOTAL_SIZE=131072
#PBS -v NMPI_DAEMON_PATH=${NMPI_ROOT}/libexec/mpid
#PBS -l elapstim_req=01:00:00
#PBS --venode=${NE}
#PBS --venum-lhost=4         # Number of VE per logical host
#PBS -l cpunum_job=2
#PBS -l coresz_prc=0
#PBS -T necmpi
#PBS --use-hca=2             # Number of HCA per logical host
##PBS -N icon-scm-real
#PBS -o /hpc/uwork/mkoehler/wq/iconCRM.%s
#PBS -j o
# ----------------------------------------------------------------------------

source ${NMPI_ROOT}/bin/necmpivars.sh

set -x

# 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
echo "PPN etc. $NN $NE $CPE $PPN"

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

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

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

# base directory for ICON sources and binary:
ICONDIR=/hpc/uhome/mkoehler/icon/icon-nwp-test1

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

# directory with input grids:
GRIDDIR=/hpc/uhome/dklocke

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

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

# debug traceback
export ATP_ENABLED=1

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


# link binary
ln -sf $MODEL icon


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

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

#-----------------------------------------------------------------------------
# loop 5day forecasts from every day 00UTC
#-----------------------------------------------------------------------------
#for gridfile in ${atmo_dyn_grids}; do

start_date=`cat start_date`
end_date=`cat end_date`

start_date="2014-02-20T00:00:00Z"    # 5day case study #DK
end_date="2014-02-20T06:00:00Z"      # ...             #DK

#start_date="2014-02-15T00:00:00Z"    # 5day from start
#end_date="2014-02-20T00:00:00Z"      # ...

#start_date="2014-02-15T00:00:00Z"   # official start time
#end_date="2014-03-23T21:00:00Z"     # official end time
#end_date="2014-03-01T00:00:00Z"
#end_date="2014-02-15T00:02:00Z"     # one time step

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

#forcings
#ln -sf ${SCMDATA}/init_data/init_SCM_GASS_DCP_GoAmazon.nc init_SCM.nc
#ln -sf ${SCMDATA}/init_data/init_SCM_GASS_DCP_GoAmazon_IOP1.nc init_SCM.nc

# cut period from input (on lce or rcl)
#  cdodir=/hpc/sw/cdo/1.9.8/x86/gnu/bin
#  ncodir=/usr/bin
#  $cdodir/cdo seldate,`echo ${start_date} | cut -c 1-19`,`echo ${end_date} | cut -c 1-19` \
#    ${SCMDATA}/init_data/init_SCM_GASS_DCP_GoAmazon_IOP1.nc init_SCM.nc
#  $ncodir/ncks -v longitude,latitude,FR_LAND,PLCOV_MX,LAI_MX,ROOTDP,RSMIN,SOILTYP,Z0,EMIS_RAD,TOPO \
#    ${SCMDATA}/init_data/init_SCM_GASS_DCP_GoAmazon_IOP1.nc init_SCM_2.nc
#  $ncodir/ncks -A init_SCM_2.nc init_SCM.nc
#  $ncodir/ncks -d time,0,360 init_SCM.nc init_SCM.nc
#  $ncodir/ncks -3 init_SCM.nc out.nc
#  $ncodir/ncrename -d time,nt out.nc
#  $ncodir/ncks -4 out.nc init_SCM.nc
#  $ncodir/ncatted -O -a units,time,a,c,'seconds since 2014-2-15 0:00:00 0:00' 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              = -3.11         ! Manaus 3.11S, 60.02W
 is_plane_torus         = .TRUE.
 l_scm_mode = .TRUE.       ! main logical to turn on SCM mode
/

&nh_testcase_nml
 nh_test_name   = 'SCM'     ! test case identifier
/

&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_random_noise  = .TRUE. !add random noise at init
 scm_sfc_mom   = 2          ! 2: prescribed u*
 scm_sfc_temp  = 2          ! 2: prescribed sensible heat flux at surface
 scm_sfc_qv    = 2          ! 2: prescribed latent heat flux at surface
!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
/

&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'
/

&nwp_phy_nml
 inwp_gscp       = 1
 inwp_convection = 0 ! 1:Tiedtke/Bechtold, 0:off
 inwp_radiation  = 1 ! 1:RRTM radiation, 4:ecRad
 inwp_cldcover   = 1 ! 3: clouds from COSMO SGS cloud scheme 0: no cloud 5: grid-scale clouds
 inwp_turb       = 1 ! 1: TKE diffusion and transfer
 inwp_satad      = 1
 inwp_sso        = 0
 inwp_gwd        = 0
 inwp_surface    = 0 ! 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        = 1
!dt_rad	         = 1800.        ! Default: 1800   ! M. Koehler: 1440
 dt_rad	         = 3600.        ! 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
/

&lnd_nml
 ntiles	         = 1       ! Default: 3
 lmulti_snow     = .false.
 itype_heatcond	 = 2
 idiag_snowfrac	 = 20      ! Default: 1 ! M. Koehler: 2
 lsnowtile       = .true.  ! later on .true. if GRIB encoding issues are solved
 lseaice         = .true.
 llake	         = .true.
 itype_lndtbl    = 4       ! minimizes moist/cold bias in lower tropical troposphere
 itype_root      = 2
 sstice_mode     = 2       ! requires extpar >= 20170202
/

&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   = .TRUE.
!is_subsidence_heat   = .FALSE.
 is_advection         = .TRUE.
 is_advection_uv      = .FALSE.
 is_advection_tq      = .TRUE.
!is_advection_tq      = .FALSE.
 is_geowind           = .FALSE.
 is_rad_forcing       = .FALSE.
 is_nudging           = .TRUE.
 is_nudging_uv        = .TRUE.
 is_nudging_t         = .FALSE.
 is_nudging_q         = .FALSE.
 dt_relax             = 3600.0         ! 10800.0
/

&turbdiff_nml
 tkhmin        = 0.75   ! new default since rev. 16527
 tkmmin        = 0.75   !           "
 pat_len       = 750.
 c_diff        = 0.2
 rat_sea       = 8.5    ! new value since rev. 25646 (reduced saturation vapor pressure over salt water)
 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    = 3
 ltkeshs       = .true.
 a_hshr        = 2.0
 !lconst_z0    = .TRUE.
 !const_z0     = 0.035
/

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

&nonhydrostatic_nml
 ivctype        = 2         ! sleve vertical coordinate
 damp_height    = 20000.    ! top_height-damp_height should be about 15km
 ndyn_substeps  = ${ndyn_substeps}
 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
 min_lay_thckn	= 20.       ! Default: 50     ! M. Koehler: 20
 max_lay_thckn	= 400.      ! maximum layer thickness below htop_thcknlimit
 htop_thcknlimit= 14000.    !
!top_height     = 75000.
 top_height     = 30000.
 stretch_fac	= 0.9       ! Default: 1      ! M. Koehler: 0.9
 decay_scale_1	= 4000.
 decay_scale_2	= 2500.
 decay_exp	= 1.2
 flat_height	= 16000.
/

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

&dynamics_nml
 divavg_cntrwgt = 0.50
!lcoriolis      = .TRUE.
 lcoriolis      = .FALSE.
/

&output_nml
 output_time_unit =  1                        ! 1: seconds
 output_bounds    =  0., 10000000., 21600.    ! start, end, increment
 mode             =  1                        ! 1: forecast
 filetype         =  4                        ! test nc4-output with 5
 steps_per_file   = ${n_in_ofile}
 include_last     = .TRUE.
 output_filename  = 'scm_out'
 filename_format  = "<output_filename>_<levtype>_<datetime2>"
 ml_varlist       = 'group:all'
 output_grid      = .TRUE.
/
&output_nml
 output_start     = "${start_date}"
 output_end       = "${end_date}"
 output_interval  = "PT01H"
!output_interval  = "PT60S"
!file_interval    = "P01H"
 mode             =  1                         ! 1: forecast
 output_time_unit =  1                         ! 1: seconds
 steps_per_file   = ${n_in_ofile}
 include_last     = .TRUE.
 output_filename  = 'scm_GASS_out'
 filename_format  = "<output_filename>_<levtype>_<datetime2>"
!ml_varlist       = 'group:all'
 ml_varlist       = 'pres_sfc','accshfl_s','acclhfl_s','accthb_s','accthb_t','accsob_s','accsob_t',
                    'accthd_s', 'accthu_s', 'accsod_s', 'accsod_t', 'accsou_t', 'accsou_t',
                    'accumfl_s','accvmfl_s','cape','rain_gsp','rain_con','snow_gsp','snow_con',
                    'u_10m','v_10m','t_2m','td_2m','t_g','qv_s','tot_prec',
                    'clct','clch','clcm','clcl','tqv_dia','tqc_dia','tqi_dia','tqr','tqs'
 pl_varlist       = 'u','v','w','temp','rho','geopot','rh','clc','tot_qv_dia','tot_qc_dia','tot_qi_dia',
                    'ddt_temp_dyn', 'ddt_temp_radsw', 'ddt_temp_radlw', 'ddt_temp_turb', 'ddt_temp_drag',
                    'ddt_temp_pconv','ddt_qv_turb', 'ddt_qc_turb', 'ddt_qi_turb', 'ddt_qv_conv',
                    'ddt_temp_gscp','ddt_qv_gscp','ddt_qc_gscp','ddt_qi_gscp'
 p_levels         =  2500, 5000, 7500,10000,12500,15000,17500,20000,22500,25000,27500,30000,
                    32500,35000,37500,40000,42500,45000,47500,50000,52500,55000,57500,60000,
                    62500,65000,67500,70000,72500,75000,77500,80000,82500,85000,87500,90000,
                    92500,95000,97500,100000
 output_grid      = .TRUE.
/
&output_nml
 output_start     = "${start_date}"
 output_end       = "${end_date}"
 output_interval  = "PT30M"
 mode             =  1                         ! 1: forecast
 output_time_unit =  1                         ! 1: seconds
 steps_per_file   = ${n_in_ofile}
 include_last     = .TRUE.
 output_filename  = 'scm_GASS_flx'
 filename_format  = "<output_filename>_<levtype>_<datetime2>"
 ml_varlist       = 'accshfl_s','acclhfl_s','accthb_s','accthb_t','accsob_s','accsob_t',
                    'accthd_s', 'accthu_s', 'accsod_t', 'accsou_t', 'accsou_t','accsod_s'
                    'accumfl_s','accvmfl_s','rain_gsp','rain_con','snow_gsp','snow_con','tot_prec'
 output_grid      = .TRUE.
/

EOF

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

ls -l $MODEL

export VE_ERRCTL_ALLOCATE=MSG
export NMPI_PROGINF=YES
export VE_TRACEBACK=ALL
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`

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

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