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

# ----------------------------------------------------------------------------
#
# ICON-ART local Pollen SPP test
#
#-----------------------------------------------------------------------------
#
# This simulates multiple pollen types and Cs-137 calculation.
#
#=============================================================================

INDIR=${basedir}/../../TESTSUITE
OUTDIR=${basedir}/experiments/horeka_art_pollen_spp_base
ICONFOLDER=${basedir}
ARTFOLDER=${basedir}/externals/art
EXP=POLLEN_SPP
lart=.True.

FILETYPE=5
restart=.False.
read_restart_namelists=.False.

# Create output directory and go to this directory

if [ ! -d $OUTDIR ]; then
    mkdir -p $OUTDIR
fi

cd $OUTDIR

ln -sf ${INDIR}/${EXP}/rrtmg_lw.nc rrtmg_lw.nc
ln -sf ${INDIR}/${EXP}/ECHAM6_CldOptProps.nc ECHAM6_CldOptProps.nc

ln -sf ${INDIR}/${EXP}/iconR3B8_DOM01_lbc.nc iconR3B8_DOM01_lbc.nc
ln -sf ${INDIR}/${EXP}/iconR3B8_DOM01.nc iconR3B8_DOM01.nc
ln -sf ${INDIR}/${EXP}/extpar_iconR3B8_DOM01.nc extpar_iconR3B8_DOM01.nc
ln -sf ${INDIR}/${EXP}/ART_POC_iconR3B08-grid_0049.grb ART_POC_iconR3B08-grid_0049.grb
ln -sf ${INDIR}/${EXP}/ART_POV_iconR3B08-grid_0049.grb ART_POV_iconR3B08-grid_0049.grb
ln -sf ${INDIR}/${EXP}/R03B08_DOM01_ilaf20210322000000.grb R03B08_DOM01_ilaf20210322000000.grb
ln -sf ${INDIR}/${EXP}/RADRRTM RADRRTM
ln -sf ${INDIR}/${EXP}/RADSRTM RADSRTM
ln -sf ${INDIR}/${EXP}/socrates_droplet_scattering_rrtm.nc socrates_droplet_scattering_rrtm.nc
ln -sf ${INDIR}/${EXP}/fu_ice_scattering_rrtm.nc fu_ice_scattering_rrtm.nc
ln -sf ${INDIR}/${EXP}/mcica_gamma.nc mcica_gamma.nc
ln -sf ${INDIR}/${EXP}/map_file.ana map_file.ana
ln -sf ${INDIR}/${EXP}/map_file.lbc map_file.lbc
ln -sf ${INDIR}/${EXP}/map_file.fc map_file.fc

ln -sf ${ARTFOLDER}/runctrl_examples/init_ctrl/PrescribedOpticalProperties_CMD_USE.nc ${OUTDIR}/prescribed_opt_props.nc


# the namelist filename
atmo_namelist=NAMELIST_${EXP}
#
#-----------------------------------------------------------------------------
# global timing
ndays_restart=60
dt_restart=`expr ${ndays_restart} \* 86400`
#
#-----------------------------------------------------------------------------
# model timing
start_date=${start_date:="2021-03-22T00:00:00Z"}
end_date=${end_date:="2021-03-22T02:30:00Z"}
output_start=${start_date:="2021-03-22T00:00:00Z"}
output_end=${end_date:="2021-03-22T02:30:00Z"}
output_interval="PT15M"
modelTimeStep="PT1M"
leadtime="PT2H30M"
checkpoint_interval="P2D"
#
#-----------------------------------------------------------------------------
# model parameters
model_equations=3             # equation system
#                     1=hydrost. atm. T
#                     1=hydrost. atm. theta dp
#                     3=non-hydrost. atm.,
#                     0=shallow water model
#                    -1=hydrost. ocean
#-----------------------------------------------------------------------------
# the grid parameters
atmo_dyn_grids="'iconR3B8_DOM01.nc'"
#-----------------------------------------------------------------------------

iconv="1"
iturb="1"
lfeedbacks=".FALSE."
num_levs="60"
dynamics_parent_grid_ids="0"
dyn_grid_files="'iconR3B8_DOM01.nc'"

# create ICON master namelist
# ------------------------
# For a complete list see Namelist_overview and Namelist_overview.pdf

cat > icon_master.namelist << EOF
&master_nml
 lRestart = .FALSE.
/
&master_time_control_nml
experimentStartDate = "$start_date"
experimentStopDate  = "$end_date"
forecastLeadTime = "$leadtime"
checkpointTimeIntval = "$checkpoint_interval"
/
&master_model_nml
  model_type		=	1
  model_name		=	"ATMO"
  model_namelist_filename="$atmo_namelist"
  model_min_rank	=	1
  model_max_rank	=	65536
  model_inc_rank	=	1
/
&time_nml
 ini_datetime_string = "$start_date"
 dt_restart          = $dt_restart
/
EOF

#-----------------------------------------------------------------------------
#

#-----------------------------------------------------------------------------
#
# write ICON namelist parameters
# ------------------------
# For a complete list see Namelist_overview and Namelist_overview.pdf
#
# ------------------------
# reconstrcuct 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


# ------------------------

cat > ${atmo_namelist} << EOF
&parallel_nml
 nproma                    = 16  ! optimal setting 8 for CRAY; use 16 or 24 for IBM
 p_test_run                = .FALSE.
 l_test_openmp             = .FALSE.
 l_log_checks              = .FALSE.
 num_io_procs              = 1   ! up to one PE per output stream is possible
 iorder_sendrecv           = 3
/
&grid_nml
 dynamics_grid_filename    = ${dyn_grid_files}
 dynamics_parent_grid_id   = ${dynamics_parent_grid_ids}
 lfeedback                 = ${lfeedbacks}
 ifeedback_type            = 2                   !
 lredgrid_phys             = .FALSE.
 l_limited_area            = .TRUE.
 start_time                = 0.
/
&initicon_nml
 init_mode                 = 7
 dwdfg_filename            = "./R<nroot0>B<jlev>_DOM<idom>_ilaf20210322000000.grb"
 ana_varnames_map_file     = "map_file.ana" 			! dictionary mapping internal names onto GRIB2 shortNames
 lread_ana                 = .FALSE.    				  ! no analysis data will be read
 ltile_coldstart           = .TRUE.       				! coldstart for surface tiles
 dt_iau                    = 600          				! time interval (in s) during which the IAU procedure is performed
 dt_shift                  = -0           				! time (in s) by which the model start is shifted ahead, given by ini_datetime_string.
 dt_ana                    = 3600
 iterate_iau               = .FALSE.
 zpbl1                     = 500.
 zpbl2                     = 1000.
 use_lakeiceana            = .TRUE.
 lp2cintp_incr             = .TRUE.
 lp2cintp_sfcana           = .TRUE.
 icpl_da_sfcevap           = 2
/
&limarea_nml
 init_latbc_from_fg        = .TRUE.
 itype_latbc               = 0
 latbc_boundary_grid       = 'iconR3B8_DOM01_lbc.nc'
 latbc_varnames_map_file   = 'map_file.lbc'
/
&run_nml
 num_lev                   = ${num_levs}
 modelTimeStep             = "${modelTimeStep}"
 ldynamics                 = .TRUE.       ! dynamics
 ltransport                = .TRUE.
 iforcing                  = 3            ! NWP forcing
 ltestcase                 = .FALSE.      !.FALSE. run with real data
 msg_level                 = 7            ! print maximum wind speeds every 5 time steps
 ltimer                    = .TRUE.       ! set .TRUE. for timer output
 timers_level              = 5            ! can be increased up to 10 for detailed timer output
 output                    = "nml"
 lart                      = .TRUE.
 check_uuid_gracefully     = .TRUE.
/
&nwp_phy_nml
 inwp_gscp                 = 1
 inwp_convection           = ${iconv}
 inwp_radiation            = 4
 inwp_cldcover             = 1
 inwp_turb                 = 1
 inwp_satad                = 1
 inwp_sso                  = 1
 inwp_gwd                  = 1
 inwp_surface              = 1
 icapdcycl                 = 3 				! apply CAPE modification to improve diurnalcycle over tropical land (optimizes NWP scores)
 latm_above_top            = .TRUE.
 efdt_min_raylfric         = 3600.
 itype_z0                  = 2
 icpl_aero_conv            = 1
 icpl_aero_gscp            = 1
 ! resolution-dependent settings - please choose the appropriate one
 dt_rad    	           = 720.
 dt_conv   	           = 180.
 dt_sso    	           = 360.
 dt_gwd    	           = 360.
 ldetrain_conv_prec    = .TRUE.
 mu_rain               = 0.5
 rain_n0_factor        = 0.1
/
&nwp_tuning_nml
 tune_zceff_min = 0.05
 itune_albedo   = 0
 tune_gkdrag      = 0.09
 tune_gkwake      = 1.8
 tune_gfrcrit     = 0.333
 tune_dust_abs    = 1.
 tune_zvz0i       = 0.85
 tune_box_liq_asy = 3.25
 tune_rhebc_land  = 0.825
 tune_rcucov      = 0.075
 tune_gfluxlaun   = 3.75e-3
 icpl_turb_clc    = 1
/
&turbdiff_nml
 tkhmin                    = 0.75  ! new default since rev. 16527
 tkhmin_strat              = 1.0
 tkmmin                    = 0.75  !
 tkmmin_strat              = 4.0
 pat_len                   = 100.
 c_diff                    = 0.2
 rlam_heat                 = 10.
 rat_sea                   = 0.875 ! ** new since r20191: 8.5 for R3B7, 8.0 for R2B6 in order to get similar temperature biases in the tropics **
 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
 alpha0                    = 0.0123
 alpha0_max                = 0.0335
 alpha1                    = 0.125
 icldm_turb                = 1
/
&lnd_nml
 ntiles                    = 3      !!! 1 for assimilation cycle and forecast
 nlev_snow                 = 3      !!! 1 for assimilation cycle and forecast
 lmulti_snow               = .FALSE. !!! .FALSE. for assimilation cycle and forecast
 itype_heatcond            = 2
 idiag_snowfrac            = 20
 lsnowtile                 = .TRUE.  !! later on .TRUE. if GRIB encoding issues are solved
 lseaice                   = .TRUE.
 llake                     = .TRUE.
 itype_lndtbl              = 3  ! minimizes moist/cold bias in lower tropical troposphere
 itype_root                = 2
 itype_canopy              = 2
 itype_evsl                = 4
 itype_trvg                = 3
 cwimax_ml                 = 5.e-4
 c_soil                    = 1.25
 c_soil_urb                = 0.5
 sstice_mode               = 2
 itype_snowevap            = 3
/
&radiation_nml
 irad_o3                   = 79
 irad_aero                 = 6
 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
 ecrad_data_path           = "."
/
&nonhydrostatic_nml
 iadv_rhotheta             = 2
 ivctype                   = 2
 itime_scheme              = 4
 exner_expol               = 0.333
 vwind_offctr              = 0.2
 damp_height               = 14000.
 rayleigh_coeff            = 1.00
 divdamp_order             = 24    ! for data assimilation runs, '2' provides extra-strong filtering of gravity waves
 divdamp_type              = 32    !!! optional: 2 for assimilation cycle if very strong gravity-wave filtering is desired
 divdamp_fac               = 0.003
 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.
 max_lay_thckn             = 400.   ! maximum layer thickness below htop_thcknlimit
 htop_thcknlimit           = 14000. ! this implies that the upcoming COSMO-EU nest will have 60 levels
 top_height                = 23000.
 stretch_fac               = 0.9
 decay_scale_1             = 4000.
 decay_scale_2             = 2500.
 decay_exp                 = 1.2
 flat_height               = 16000.
/
&dynamics_nml
 iequations                = 3
 divavg_cntrwgt            = 0.50
 lcoriolis                 = .TRUE.
/
&diffusion_nml
 hdiff_order               = 5
 itype_vn_diffu            = 1
 itype_t_diffu             = 2
 hdiff_efdt_ratio          = 30.0
 hdiff_smag_fac            = 0.02
 lhdiff_vn                 = .TRUE.
 lhdiff_temp               = .TRUE.
/
&interpol_nml
 nudge_zone_width  	   = 8
 lsq_high_ord      	   = 3
 l_intp_c2l        	   = .TRUE.
 l_mono_c2l        	   = .TRUE.
 support_baryctr_intp = .false.
/
&gribout_nml
 preset                          = "deterministic"
 ldate_grib_act                  = .TRUE.
 lgribout_24bit                  = .FALSE.
 backgroundProcess               = 0
 localNumberOfExperiment         = 0
 productionStatusOfProcessedData = 1
 localTypeOfEnsembleForecast     = -1
 numberOfForecastsInEnsemble     = -1
 perturbationNumber              = -1
 generatingProcessIdentifier     = 11
 tablesVersion                   = 19
/
&extpar_nml
 itopo          	         = 1
 n_iter_smooth_topo 	     = 1
 heightdiff_threshold      = 3000.
/
&io_nml
 itype_pres_msl            = 5  ! IFS method with bug fix for self-consistency between SLP and geopotential
 itype_rh                  = 1  ! RH w.r.t. water
 output_nml_dict           = 'map_file.fc'
 lmask_boundary            = .FALSE.
/
&output_nml
 filetype                  =  ${FILETYPE}        ! output format: 2=GRIB2, 4=NETCDFv2
 dom                       =  1                 ! write all domains
 output_start    	         = "${output_start}"
 output_end                = "${output_end}"
 output_interval           = "${output_interval}"
 mode                      =  1
 steps_per_file            =  1
 include_last              =  .TRUE.
 output_filename           = 'icon-art-${EXP}-aero-reg'            ! file name base
 ml_varlist                = 'z_ifc','rho','temp', 'group:ART_AEROSOL'
 remap                     =  1
 reg_lon_def               = -15.,0.25,35.
 reg_lat_def               =  30.,0.25,65.
/
&art_nml
 lart_diag_out             = .TRUE.
 lart_aerosol              = .TRUE.
 lart_chem                 = .FALSE.
 lart_chemtracer           = .FALSE.
 lart_pntSrc               = .TRUE.
 iart_ari                  = 0
 iart_aci_warm             = 0
 iart_aci_cold             = 0
 iart_init_aero            = 0
 iart_seasalt              = 0
 iart_dust                 = 0
 iart_fire                 = 0
 iart_volcano              = 0
 iart_pollen               = 1
 iart_radioact             = 1
 iart_aero_washout         = 1
 iart_nonsph	             = 0
 iart_isorropia            = 0
 cart_aerosol_xml          = '${ARTFOLDER}/runctrl_examples/xml_ctrl/${EXP}/tracers_poll_spp.xml'
 cart_input_folder         = '${OUTDIR}'
 cart_modes_xml            = '${ARTFOLDER}/runctrl_examples/xml_ctrl/${EXP}/modes_poll_spp.xml'
 cart_diagnostics_xml      = '${ARTFOLDER}/runctrl_examples/xml_ctrl/${EXP}/diagnostics_poll_spp.xml'
 cart_pntSrc_xml           = '${ARTFOLDER}/runctrl_examples/xml_ctrl/${EXP}/pntSrc_poll_spp.xml'
 cart_opt_props_nc         = '${OUTDIR}/prescribed_opt_props.nc'
/

EOF

rm -f finish.status

#-----------------
# start the model
#-----------------
ls -l ${MODEL}
check_error $? "${MODEL} does not exist?"
#
ldd ${MODEL}

date
$START $MODEL || exit 1
date

#-----------------------------------------------------------------------------
# Check run status
#-----------------------------------------------------------------------------
#
if [[ -f finish.status ]] ; then
  finish_status=$(echo $(cat finish.status))
else
  echo "Run failed or did not complete properly"
  exit 1
fi

echo "Finish status: ${finish_status}"

#-----------------------------------------------------------------------------
# exit 0
#-----------------------------------------------------------------------------
