# 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 Aerosols Test
#
#-----------------------------------------------------------------------------
#
# This simulation performs tests for various types of aerosols (dust, sea salt, wildfire).
#
#=============================================================================

INDIR=${basedir}/../../TESTSUITE
OUTDIR=${basedir}/experiments/horeka_art_allaero_norad_base
ICONFOLDER=${basedir}
ARTFOLDER=${basedir}/externals/art
EXP=ALLAERO_NORAD
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/rrtmg_lw.nc rrtmg_lw.nc
ln -sf ${ICONFOLDER}/data/ECHAM6_CldOptProps.nc ECHAM6_CldOptProps.nc

# Grid and external parameter files
ln -sf ${INDIR}/${EXP}/iconR2B06_DOM01.nc iconR2B06_DOM01.nc
ln -sf ${INDIR}/${EXP}/iconR2B06_DOM01-grfinfo.nc iconR2B06_DOM01-grfinfo.nc
ln -sf ${INDIR}/${EXP}/extpar_iconR2B06_DOM01.nc extpar_iconR2B06_DOM01.nc
ln -sf ${INDIR}/${EXP}/ART_IAE_iconR2B06-grid_TCNR.nc ART_IAE_iconR2B06-grid_TCNR.nc
ln -sf ${INDIR}/${EXP}/ART_STY_iconR2B06-grid_TCNR.nc  ART_STY_iconR2B06-grid_TCNR.nc

# Input GRIB2 analysis and forecast data from DWD
ln -sf ${INDIR}/${EXP}/dwdFG_R2B06_DOM01.grb dwdFG_R2B06_DOM01.grb
ln -sf ${INDIR}/${EXP}/dwdana_R2B06_DOM01.grb dwdana_R2B06_DOM01.grb
ln -sf ${INDIR}/${EXP}/ART_BCF_iconR2B06-grid_TCNR.nc ART_BCF_iconR2B06-grid-2019-06-21_TCNR.nc
ln -sf ${ICONFOLDER}/run/ana_varnames_map_file.txt map_file.ana

# Prescribed aerosol optical properties
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:="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="P1D"
#
#-----------------------------------------------------------------------------
# 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="'iconR2B06-DOM01.nc"
#-----------------------------------------------------------------------------

iconv="1"
iturb="1"
lfeedbacks=".false."
num_levs="90"
dynamics_parent_grid_ids="0"
dyn_grid_files="'iconR2B06_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                    = 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.
/
&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            = 0            ! 0 = no radiation
 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.  ! 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
 tkhmin                    = 0.75  ! new default since rev. 16527
 tkmmin                    = 0.75  !
 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
/
&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                 = 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           = "${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}"
 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               = -180.,0.5,179.5
 reg_lat_def               = 90.,-0.5, -90.
/
&art_nml
 lart_diag_out             = .TRUE.
 lart_aerosol              = .TRUE.
 iart_ari                  = 0
 iart_init_aero            = 5
 iart_seasalt              = 0
 iart_dust                 = 1
 iart_fire                 = 1
 iart_volcano              = 0
 iart_aero_washout         = 1
 iart_nonsph	           = 0
 iart_isorropia            = 0
 cart_io_suffix            = 'TCNR'
 cart_aerosol_xml          = '${ARTFOLDER}/runctrl_examples/xml_ctrl/${EXP}/tracers_allaero_norad.xml'
 !cart_coag_xml             = '${ARTFOLDER}/runctrl_examples/xml_ctrl/${EXP}/coagulate_allaero_norad.xml'
 cart_input_folder         = '.'
 cart_modes_xml            = '${ARTFOLDER}/runctrl_examples/xml_ctrl/${EXP}/modes_allaero_norad.xml'
 cart_aero_emiss_xml       = '${ARTFOLDER}/runctrl_examples/xml_ctrl/${EXP}/aero_emiss_allaero_norad.xml'
 cart_diagnostics_xml      = '${ARTFOLDER}/runctrl_examples/xml_ctrl/${EXP}/diagnostics_allaero_norad.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
#-----------------------------------------------------------------------------
