# 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 DUST Test
#
#-----------------------------------------------------------------------------
#
# This simulation performs dust emissions and radiative effects.
#
#=============================================================================

INDIR=${basedir}/../../TESTSUITE
OUTDIR=${basedir}/experiments/horeka_art_dust_rad_base
ICONFOLDER=${basedir}
ARTFOLDER=${basedir}/externals/art
EXP=DUST_RAD
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


#-----------------------------------------------------------------------------
# Input files, grid and external parameter files, initial/forcing GRIB2 data
#-----------------------------------------------------------------------------
#
ln -sf ${ICONFOLDER}/data/ECHAM6_CldOptProps.nc ECHAM6_CldOptProps.nc
ln -sf ${ICONFOLDER}/data/rrtmg_lw.nc rrtmg_lw.nc

ln -sf $INDIR/$EXP/icon_grid_0024_R02B06_G.nc iconR2B06_DOM01.nc
ln -sf $INDIR/$EXP/icon_grid_0024_R02B06_G-grfinfo.nc iconR2B06_DOM01-grfinfo.nc
ln -sf $INDIR/$EXP/icon_extpar_0024_R02B06_G_20180625_tiles.nc extpar_iconR2B06_DOM01.nc

ln -sf $INDIR/$EXP/domain2_DOM01.nc iconR2B07_DOM02.nc
ln -sf $INDIR/$EXP/external_parameter_icon_domain2_DOM01_tiles.nc extpar_iconR2B07_DOM02.nc

ln -sf $INDIR/$EXP/igfff00013000.grb dwdFG_R2B06_DOM01.grb
ln -sf $INDIR/$EXP/igaf20190622000000.grb dwdana_R2B06_DOM01.grb
ln -sf ${ICONFOLDER}/run/ana_varnames_map_file.txt map_file.ana

ln -sf $INDIR/$EXP/igfff00000000_dust.nc ART_IAE_iconR2B06-grid_TCNR.nc
ln -sf $INDIR/$EXP/ART_STY_iconR2B06-grid_TCNR.nc ART_STY_iconR2B06-grid_TCNR.nc
ln -sf $INDIR/$EXP/ART_STY_iconR2B07-grid_TCNR.nc ART_STY_iconR2B07-grid_TCNR.nc

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=1
dt_restart=`expr ${ndays_restart} \* 86400`
#
#-----------------------------------------------------------------------------
# model timing
start_date=${start_date:="2019-06-22T00:00:00Z"}
end_date=${end_date:="2019-06-22T03:00:00Z"}
output_start=${start_date:="2019-06-22T00:00:00Z"}
output_end=${end_date:="2019-06-22T03:00:00Z"}
output_interval="PT1H"
modelTimeStep="PT6M"
leadtime="PT3H"
checkpoint_interval="P10D"
#
#-----------------------------------------------------------------------------
# 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
#-----------------------------------------------------------------------------

iconv="1,1"
iturb="1,1"
lfeedbacks=".false.,.false."
num_levs="90,60"
dynamics_parent_grid_ids="0,1"
dyn_grid_files="'iconR2B06_DOM01.nc','iconR2B07_DOM02.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
#
# ------------------------

cat > ${atmo_namelist} << EOF
&parallel_nml
 nproma         = 8  ! 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  ! best value for CRAY (slightly faster than option 1)
/
&grid_nml
 dynamics_grid_filename  = ${dyn_grid_files}
 dynamics_parent_grid_id = ${dynamics_parent_grid_ids}
 lfeedback               = ${lfeedbacks}
 ifeedback_type          = 2                   !
!lredgrid_phys           = .false.
!radiation_grid_filename = ${radiation_grid_filename}
 lredgrid_phys           = .false.
 start_time              = 0.,10800.,14400.,16200.,18000.
/
&initicon_nml
 init_mode      		= 5           				! operation mode 2: IFS
 dwdfg_filename  		= "dwdFG_R2B06_DOM01.grb" 	! initial data filename
 dwdana_filename        = "dwdana_R2B06_DOM01.grb"
 ana_varnames_map_file  = "map_file.ana" 			! dictionary mapping internal names onto GRIB2 shortNames
 lread_ana      	 	= .TRUE.     				! no analysis data will be read
 ltile_coldstart 		= .FALSE.      				! coldstart for surface tiles
 ltile_init      		= .FALSE.     				! set it to .TRUE. if FG data originate from run without tiles
 dt_iau      			= 10800.       				! time interval (in s) during which the IAU procedure is performed
 dt_shift    			= -5400.       				! time (in s) by which the model start is shifted ahead, given by ini_datetime_string.
 iterate_iau 			= .true.
 zpbl1       			= 500.
 zpbl2       			= 1000.
 !use_lakeiceana 		= .true.
 lp2cintp_incr=.true.,
 lp2cintp_sfcana=.true.,
 check_ana(1)%list      = 'T_SO','W_SO','P','QV','T','U','V','FRESHSNW','H_SNOW'
 check_ana(2)%list      = 'T_SO','W_SO','FRESHSNW','H_SNOW'
/
&run_nml
 num_lev        = ${num_levs}
 lvert_nest     = .true.       ! use vertical nesting if a nest is active
 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   = 1            ! 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  = .false., .true.  ! the second entry refers to the nested domain (if present)
 efdt_min_raylfric = 7200.
 itype_z0         = 2
 icpl_aero_conv   = 0
 icpl_aero_gscp   = 0
 ! resolution-dependent settings - please choose the appropriate one
 dt_rad    		= 1440.
 dt_conv   		= 360.
 dt_sso    		= 720.
 dt_gwd    		= 720.
/
&nwp_tuning_nml
 tune_zceff_min = 0.05 ! ** default value to be used for R3B7; use 0.05 for R2B6 in order to get similar temperature biases in upper troposphere **
 itune_albedo   = 1     ! somewhat reduced albedo (w.r.t. MODIS data) over Sahara in order to reduce cold bias
/
&turbdiff_nml
 a_hshr        = 2.0
 frcsmot       = 0.2   ! these 2 switches together apply vertical smoothing of the TKE source terms
 icldm_turb    = 2
 imode_frcsmot = 2     ! in the tropics (only), which reduces the moist bias in the tropical lower troposphere
 imode_tkesso  = 2
! use horizontal shear production terms with 1/SQRT(Ri) scaling to prevent unwanted side effects:
 itype_sher    = 2
 ltkeshs       = .true.
 ltkesso       = .true.
 pat_len       = 750.
 q_crit        = 2.0
 rat_sea       = 0.8 ! :cp-cv-bugfix, old value: 7.0
 tkhmin        = 0.5
 tkmmin        = 0.75
 tur_len       = 300.
! cp-cv-bugfix:
 rlam_heat     = 10.0
 alpha1        = 0.125
/
&lnd_nml
 ntiles         = 3      !!! 1 for assimilation cycle and forecast
 nlev_snow      = 1      !!! 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.
 frlake_thrhld	= 0.05
 itype_lndtbl   = 3  ! minimizes moist/cold bias in lower tropical troposphere
 itype_root     = 2
/
&radiation_nml
 irad_o3       = 7
 irad_aero     = 9
 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 = "${ICONFOLDER}/externals/ecrad/data"
/
&nonhydrostatic_nml
 iadv_rhotheta  = 2
 ivctype        = 2
 itime_scheme   = 4
 exner_expol    = 0.333
 vwind_offctr   = 0.2
 damp_height    = 50000.
 rayleigh_coeff = 0.10
 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.004
 igradp_method  = 3
 l_zdiffu_t     = .true.
 thslp_zdiffu   = 0.02
 thhgtd_zdiffu  = 125.
 htop_moist_proc= 22500.
 hbot_qvsubstep = 22500. ! use 19000. with R3B7
/
&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      = 75000.
 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.
/
&transport_nml
 ctracer_list  = '12345'
 ivadv_tracer  = 3,3,3,3,3,3,3,3,3,3,3
 itype_hlimit  = 3,4,4,4,4,3,3,3,3,3,3
 ihadv_tracer  = 52,2,2,2,2,22,22,22,22,22,22
 iadv_tke      = 0
/
&diffusion_nml
 hdiff_order      = 5
 itype_vn_diffu   = 1
 itype_t_diffu    = 2
 hdiff_efdt_ratio = 24.0
 hdiff_smag_fac   = 0.025
 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
 generatingCenter 		= 78
 generatingSubcenter 	= 255
/
&extpar_nml
 itopo          		= 1
 n_iter_smooth_topo 	= 1
 heightdiff_threshold 	= 3000.
/
&io_nml
 itype_pres_msl = 4  ! IFS method with bug fix for self-consistency between SLP and geopotential
 itype_rh       = 1  ! RH w.r.t. water
/
&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}"
 steps_per_file               =  1
 include_last                 =  .TRUE.
 output_filename              = 'icon-art-${EXP}-aero-reg'            ! file name base
 ml_varlist                   = 'temp','sob_s','sodifd_s','z_mc','group:ART_AEROSOL'
 remap                        =  1
 reg_lon_def                  = -180.,0.5,179.5
 reg_lat_def                  = 90.,-0.5, -90.
/
&output_nml
 filetype                     =  ${FILETYPE}        ! output format: 2=GRIB2, 4=NETCDFv2
 dom                          =  2                 ! write all domains
 output_start    			  = "${output_start}"
 output_end                   = "${output_end}"
 output_interval              = "${output_interval}"
 steps_per_file               =  1
 include_last                 =  .TRUE.
 output_filename              = 'icon-art-${EXP}-aero-reg'            ! file name base
 ml_varlist                   = 'temp','sob_s','sodifd_s','z_mc','group:ART_AEROSOL'
 remap                        =  1
 reg_lon_def                  = -30.,0.25,45
 reg_lat_def                  = 60.,-0.25, 0.
/
&art_nml
 lart_diag_out   = .TRUE.
 lart_aerosol    = .TRUE.
 iart_ari        = 1
 iart_init_aero  = 5
 iart_seasalt    = 0
 iart_dust       = 1
 iart_nonsph	 = 1
 cart_io_suffix  = 'TCNR','TCNR'
 cart_aerosol_xml = '${ARTFOLDER}/runctrl_examples/xml_ctrl/${EXP}/tracers_dustrad.xml'
 cart_input_folder = '${OUTDIR}'
 cart_modes_xml  = '${ARTFOLDER}/runctrl_examples/xml_ctrl/${EXP}/modes_dustrad.xml'
 cart_coag_xml   = '${ARTFOLDER}/runctrl_examples/xml_ctrl/${EXP}/coagulate_dustrad.xml'
 cart_aero_emiss_xml = '${ARTFOLDER}/runctrl_examples/xml_ctrl/${EXP}/aero_emiss_dustrad.xml'
 cart_diagnostics_xml = '${ARTFOLDER}/runctrl_examples/xml_ctrl/${EXP}/diagnostics_aerosol_with_grib2.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
#-----------------------------------------------------------------------------
