nwpexp.run_ICON_27_R2B6_mvstream_qbudget

# ICON
#
# ------------------------------------------
# Copyright (C) 2004-2024, 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
# ------------------------------------------

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

make_and_change_to_experiment_dir # function in ../add_run_routines

# Combine START and MODEL if START_MODEL is not already set.
# START_MODEL is used to ease the execution of a machine that uses a complex
# mpirun command with multiple binaries
START_MODEL="${START_MODEL:=$START $MODEL}"

# set icon_data_poolFolder
icon_data_poolFolder="${icon_data_rootFolder:-/pool/data/ICON}/buildbot_data/nwp"

# base directory for ICON sources and binary:
ICONDIR="${basedir}" # basedir can be set with create_target_header

EXP="NWP"    # experiment identifier
EDIR="exp16"
# directory for grid and extpar files
GRIDDIR="${icon_data_poolFolder}/grids/public/edzw"
# GRIDDIR links to /hpc/rwork0/routfor/routfox/icon/grids/public/edzw on RCL
EXTPDIR="$GRIDDIR" # external parameter directory

# absolute path to input directory - for the time being
INDIR=${icon_data_poolFolder}/Checksuite_data/EPS_R2B6N7

# ----------------------------------------------------------------------
# copy input data: grids, external parameters
# ----------------------------------------------------------------------

# the directory for the experiment will be created, if not already there
if [ ! -d $EXPDIR ]; then
    mkdir -p $EXPDIR
fi
cd $EXPDIR

ln -sf $GRIDDIR/icon_grid_0023_R02B05_R.nc iconR2B05_DOM00.nc
ln -sf $GRIDDIR/icon_grid_0023_R02B05_R-grfinfo.nc iconR2B05_DOM00-grfinfo.nc
ln -sf $GRIDDIR/icon_grid_0024_R02B06_G.nc iconR2B06_DOM01.nc
ln -sf $GRIDDIR/icon_grid_0024_R02B06_G-grfinfo.nc iconR2B06_DOM01-grfinfo.nc
ln -sf $GRIDDIR/icon_grid_0028_R02B07_N02.nc iconR2B07_DOM02.nc
ln -sf $GRIDDIR/icon_grid_0028_R02B07_N02-grfinfo.nc iconR2B07_DOM02-grfinfo.nc

ln -sf $EXTPDIR/icon_extpar_0024_R02B06_G_20180625_tiles.nc extpar_iconR2B06_DOM01.nc
ln -sf $EXTPDIR/icon_extpar_0028_R02B07_N02_20180625_tiles.nc extpar_iconR2B07_DOM02.nc


cymdg=2017041500
# `datconv` is only available on the dwd machine. Thus we give cymdgm3 directly.
# cymdgm3=$(datconv -Cy -d -3 ${cymdg})
cymdgm3=2017041421

echo dates $cymdg $cymdgm3


# files needed for radiation
ln -sf ${ICONDIR}/externals/ecrad/data ecrad_data


# ----------------------------------------------------------------------
# global namelist settings
# ----------------------------------------------------------------------

# the namelist filename
atmo_namelist=NAMELIST_${EXP}

# global timing
start_date="2017-04-14T22:30:00Z"
ndays_restart=60
dt_restart=`expr ${ndays_restart} \* 86400`

dtime=360
nsteps=30
nsteps_per_hour=10

# the grid parameters
atmo_dyn_grids="iconR2B06_DOM01.nc iconR2B07_DOM02.nc"
atmo_rad_grids="iconR2B05_DOM00.nc"

# output configuration
output_bounds=(
  "0., $((dtime * nsteps))., ${dtime}"
  "0., $((dtime * nsteps))., 3600."
  "0., $((dtime * nsteps))., ${dtime}"
  )
output_operation=(
  ""
  "mean"
  ""
  )
output_filename=(
  "raw"
  "hr_mean"
  "raw_gscp3"
  )

inwp_gscp=(
  1
  1
  3
  )
extra_vars_flux=(
  ""
  ""
  ""
  )
extra_vars_storage=(
  ""
  ""
  ""
  )

cp -p $MODEL icon


 ln -sf $INDIR/igfff${cymdgm3}-0130.grb fg_R02B06_DOM01.grb
 ln -sf $INDIR/iefff${cymdgm3}-0130.grb fg_R02B07_DOM02.grb



# ----------------------------------------------------------------------
# create ICON master namelist
# ----------------------------------------------------------------------

# For a complete list see Namelist_overview and Namelist_overview.pdf

cat > icon_master.namelist << EOF
&master_nml
 lrestart               = .false.
/
&time_nml
 ini_datetime_string = "$start_date"
 dt_restart          = $dt_restart
/
&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
/
EOF


# ----------------------------------------------------------------------
# model namelists
# ----------------------------------------------------------------------

# 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

# BEGIN RUN LOOP
for irun in 0 1 2; do

# Variables to configure the mpirun command on the RCL
num_io_procs=2
num_restart_procs=0
num_prefetch_proc=0
num_io_procs_radar=0

cat > ${atmo_namelist} << EOF
&parallel_nml
 nproma          =  ${nproma}
 nproma_sub      =  ${nproma_sub} ! loop chunk length for radiation
 p_test_run      = .false.
 l_test_openmp   = .true.
 l_log_checks    = .true.
 num_io_procs       = ${num_io_procs}
 num_restart_procs  = ${num_restart_procs}
 num_prefetch_proc  = ${num_prefetch_proc}
 num_io_procs_radar = ${num_io_procs_radar}
 iorder_sendrecv =  3
 proc0_shift     = ${proc0_shift}
 use_omp_input   = .true.
/
&grid_nml
 dynamics_grid_filename  = ${dynamics_grid_filename}
 radiation_grid_filename = ${radiation_grid_filename}
 dynamics_parent_grid_id = 0,1
 lredgrid_phys           = .true.,.true.,
 lfeedback               = .false.        !!!.true.
 ifeedback_type          = 2
 start_time              = 0., -5400.,
 end_time                = 0., 432000.    ! ** bewirkt, dass das Nest nur bis vv=120h laueft **
/
&run_nml
 num_lev        = 90, 60
 lvert_nest     = .true.
 nsteps         = ${nsteps}
 dtime          = ${dtime}     ! timestep in seconds
 ldynamics      = .TRUE.       ! dynamics
 ltransport     = .true.
 ntracer        = 5            ! default: 0
 iforcing       = 3
 ltestcase      = .false.
 msg_level      = 7 ! detailed report during integration
 ltimer         = .true.
 timers_level   = 10
 check_uuid_gracefully = .true.
 output         = "nml" ! "none" to turn off output
/
&initicon_nml
 init_mode                    = 5
 lread_ana                    = .false.
 dt_iau                       = 0
 iterate_iau                  = .false.
 zpbl1                        = 500.
 zpbl2                        = 1000.
 dwdfg_filename               = "<path>fg_R<nroot0>B<jlev>_DOM<idom>.grb"
 ana_varnames_map_file        = "map_file.ana"
 ltile_coldstart              = .false.
 lvert_remap_fg               = .false.
 use_lakeiceana               = .true.
 lp2cintp_incr                = .TRUE.
 lp2cintp_sfcana              = .TRUE.
/
&io_nml
 itype_pres_msl               = 5        ! (1) 3: IFS-type extrapolation
 itype_rh                     = 1         ! (1) 2: mixed phase (water and ice)
 output_nml_dict              = 'map_file.fc'
 dt_checkpoint           = 2592000.  ! 30 days
/
&nwp_phy_nml
inwp_gscp       = ${inwp_gscp[irun]}
inwp_convection = 1
inwp_radiation  = 4
inwp_cldcover   = 1 ! 3 geaendert ab Test2
inwp_turb       = 1
inwp_satad      = 1
inwp_sso        = 1
inwp_gwd        = 1
inwp_surface    = 1
latm_above_top  = .false.,.true.
ldetrain_conv_prec = .false.,.false.
efdt_min_raylfric = 7200.
itype_z0         = 2
icapdcycl        = 3
icpl_aero_conv   = 1
icpl_aero_gscp   = 1
icpl_o3_tp       = 1
dt_rad           = 2160.
dt_conv          = 720.,360.
dt_sso           = 1440.,720.
dt_gwd           = 1440.,
/
&nwp_tuning_nml
itune_albedo  = 1
tune_gkwake   = 1.5
tune_gfrcrit  = 0.425
tune_grcrit   = 0.5
tune_gkdrag   = 0.16
tune_zvz0i    = 0.85
tune_box_liq_asy = 3.25
tune_minsnowfrac = 0.2
tune_gfluxlaun  = 3.75e-3
tune_rcucov = 0.075
tune_rhebc_land = 0.825
/
&turbdiff_nml
tkhmin       = 0.6
tkhmin_strat = 1.0
tkmmin       = 0.75
tkmmin_strat = 1.5
pat_len      =  750.
c_diff       =  0.2
rat_sea      =  0.8
rlam_heat    = 10.
ltkesso      = .true.
frcsmot      = 0.2
imode_frcsmot = 2
icldm_turb   = 1
itype_sher   = 1
ltkeshs      = .true.
a_hshr       = 2.0
/
&lnd_nml
  ntiles         = 3
  nlev_snow      = 1
  lmulti_snow    = .false.
  itype_heatcond = 3
  idiag_snowfrac = 20
  itype_snowevap = 2
  itype_canopy   = 2
  lsnowtile      = .true.
  lseaice        = .true.
  llake          = .true.
  itype_lndtbl   = 4
  itype_evsl     = 4
  itype_root     = 2
  cwimax_ml      = 5.e-4
  c_soil         = 1.25
  c_soil_urb     = 0.5
  sstice_mode    = 2
  lprog_albsi    = .true.
/
&radiation_nml
 irad_o3     = 79
 irad_aero   = 6
 albedo_type = 2 ! Modis albedo
 direct_albedo = 4
 direct_albedo_water = 3
 albedo_whitecap = 1
 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 = './ecrad_data'
 isolrad    = 1
 ecrad_isolver = ${radiation_ecrad_isolver}
 ! Solver version. 2: McICA with OpenACC, 0: McICA
/
&nonhydrostatic_nml
 itime_scheme    = 4
 vwind_offctr    = 0.2
 damp_height     = 44000.
 rayleigh_coeff  = 0.5
 divdamp_order   = 24
 divdamp_fac     = 0.004
 divdamp_type    = 32
 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.
 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      = .TRUE.
/
&transport_nml
  ctracer_list  = '12345'
 ivadv_tracer   = 3,3,3,3,3
 itype_hlimit   = 3,4,4,4,4,
 ihadv_tracer   = 52,2,2,2,2,
 llsq_svd       = .false.
 beta_fct       = 1.005
 iadv_tke       = 0
/
&diffusion_nml
 hdiff_order      = 5
 itype_vn_diffu   = 1
 itype_t_diffu    = 2
 hdiff_efdt_ratio = 24.
 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.
rbf_scale_mode_ll = 2
support_baryctr_intp=.true.,
/
&gridref_nml
 grf_intmethod_e  = 6
 grf_scalfbk      = 2
 denom_diffu_v    = 150.
/
&extpar_nml
 itopo                = 1
 n_iter_smooth_topo   = 1,1
 hgtdiff_max_smooth_topo = 750.,750.,
 heightdiff_threshold = 3000.
 itype_vegetation_cycle = 2
/
&ensemble_pert_nml
use_ensemble_pert=.false.,
/
&gribout_nml
 preset                          = "ensemble"
 perturbationNumber              = 1
 localTypeOfEnsembleForecast     = 101,
 numberOfForecastsInEnsemble     = 40,
 ldate_grib_act                  = .true.
 lgribout_24bit                  = .false.
 tablesVersion                   = 14
 backgroundProcess               = 0
 localNumberOfExperiment         = 52
 productionStatusOfProcessedData = 2
 generatingProcessIdentifier     = 1, 2                              ! 2 .. nest
/
&synsat_nml
lsynsat=.false. ! .false.,.true., ! requires RTTOV library - left out for this test
/
&output_nml
 filetype  = 4
 dom       = -1
 remap     = 0
 operation = '${output_operation[irun]}'
 mode      = 1
 steps_per_file = ${nsteps}
 output_bounds = ${output_bounds[irun]}
 output_filename = '${output_filename[irun]}'
 ml_varlist = 'group:dwd_fg_atm_vars', 'rh', 'pres_msl', 'tqv'
/
&output_nml
 filetype  = 4
 dom       = -1
 remap     = 0
 mode      = 1
 steps_per_file = 1
 output_bounds = 0., 0., 9999.
 output_filename = 'initial'
 ml_varlist = 'group:dwd_fg_atm_vars', 'rh', 'pres_msl', 'tqv'
/
! native output for water conservation testing  - after IAU - operation_mean fluxes
&output_nml
 output_filename         = 'water_budget_flux_${irun}'       ! file name base
 steps_per_file          = ${nsteps}
 output_bounds           = 5400., $((dtime * nsteps))., ${dtime}
 operation               = "mean"                    ! only working with remap=0, not accumulated!
 filetype                = 4                         ! output format: 2=GRIB2, 4=NETCDFv2
 dom                     = 1
 mode                    = 1                         ! 1: forecast mode (relative t-axis)
 remap                   = 0                         ! 1: latlon,  0: native grid
 output_grid             = .TRUE.
 ml_varlist              = 'rain_gsp_rate','snow_gsp_rate','ice_gsp_rate',
                           'rain_con_rate_corr','snow_con_rate_corr','tot_prec_rate',
                           'qhfl_s','qcfl_s','qifl_s'${extra_vars_flux[irun]}
/
! native output for water conservation testing - after IAU - normal variables
&output_nml
 output_filename         = 'water_budget_state_${irun}'      ! file name base
 steps_per_file          = ${nsteps}
 output_bounds           = 5400., $((dtime * nsteps))., ${dtime}
 filetype                = 4                         ! output format: 2=GRIB2, 4=NETCDFv2
 dom                     = 1
 mode                    = 1                         ! 1: forecast mode (relative t-axis)
 remap                   = 0                         ! 1: latlon,  0: native grid
 output_grid             = .TRUE.
 ml_varlist              = 'tqv','tqc','tqi','tqs','tqr'${extra_vars_storage[irun]}
                         ! 'tqg'
/
EOF

cat > map_file.ana << EOF2
# internal name     GRIB2 shortName
theta_v             THETA_V
rho                 DEN
vn                  VN
u                   U
v                   V
w                   W
tke                 TKE
temp                T
pres                P
qv                  QV
qc                  QC
qi                  QI
qr                  QR
qs                  QS
t_g                 T_G
qv_s                QV_S
fr_seaice           FR_ICE
alb_si              ALB_SEAICE
t_ice               T_ICE
h_ice               H_ICE
t_snow              T_SNOW
freshsnow           FRESHSNW
snowfrac_lc         SNOWC
w_snow              W_SNOW
rho_snow            RHO_SNOW
h_snow              H_SNOW
w_i                 W_I
w_so                W_SO
w_so_ice            W_SO_ICE
t_so                T_SO
smi                 SMI
gz0                 Z0
pres_sfc            PS
z_ifc               HHL

t_mnw_lk            T_MNW_LK
t_wml_lk            T_WML_LK
h_ml_lk             H_ML_LK
t_bot_lk            T_BOT_LK
c_t_lk              C_T_LK
t_b1_lk             T_B1_LK
h_b1_lk             H_B1_LK
EOF2


cat > map_file.fc << EOF3
# GRIB2 shortName   internal name
  alb_rad           ALBDIF
  alb_seaice        ALB_SI
  apab_s            ASWFLX_PAR_SFC
  aswdifd_s         ASODIFD_S
  aswdifu_s         ASODIFU_S
  aswdir_s          ASODIRD_S
  cape_con          CAPE
  den               RHO
  dtke_con          DDT_TKE_PCONV
  dtke_hsh          DDT_TKE_HSH
  fi                GEOPOT
  fr_ice            FR_SEAICE
  freshsnw          FRESHSNOW
  hhl               Z_IFC
  hsurf             TOPOGRAPHY_C
  p                 PRES
  pmsl              PRES_MSL
  ps                PRES_SFC
  relhum            RH
  relhum_2m         RH_2M
  snowc             SNOWFRAC_LC
  sobs_rad          SOB_S
  t                 TEMP
  t_sea             T_SEASFC
  thbs_rad          THB_S
  vmax_10m          GUST10
  z0                GZ0
EOF3



#
# configure START_MODEL_function
#
ICON_COMPONENT1_VH_procs=$((num_restart_procs + num_io_procs + num_prefetch_proc + num_io_procs_radar))

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

$START_MODEL
EXIT_STATUS=$?
if [ "$EXIT_STATUS" -ne "0" ]; then
  echo "ICON EXIT_STATUS: $EXIT_STATUS"
  exit $EXIT_STATUS
fi

done
# END RUN LOOP

# ----------------------------------------------------------------------
# EPILOGUE: Compare the output to a given reference
#     "... in every code we have some trouble
#          But when you worry you make it double ..."
# ----------------------------------------------------------------------

set +x

# ----------------------------------------------------------------------
# Test water budget
# ----------------------------------------------------------------------

echo 'Test water conservation ...'

$ICONDIR/run/checksuite.nwp/nwpexp.run_ICON_27_R2B6_mvstream_qbudget.check.sh

EXIT_STATUS=$?
if [ "$EXIT_STATUS" -ne "0" ]; then
  echo "Water conservation test not passed: $EXIT_STATUS"
  exit $EXIT_STATUS
fi


# ----------------------------------------------------------------------
# Test mvstream output
# ----------------------------------------------------------------------

echo 'Comparing means from raw output and mvstream...'

result=()
diffs=()
for dom in 01 02; do

  echo "DOM$dom: preparing time means from raw outputs..."
  $cdo -s timmean -seltimestep,2/$((nsteps+1)) raw_DOM${dom}_ML_0001.nc cdo_tim_mean_raw_DOM${dom}.nc > /dev/null
  $cdo -s timmean -seltimestep,2/$((nsteps/nsteps_per_hour+1)) hr_mean_DOM${dom}_ML_0001.nc cdo_tim_mean_hr_mean_DOM${dom}.nc > /dev/null

  echo "DOM$dom: checking pressure and temperature..."
  set +e
  diff="$($cdo diffn,names=left,abslim=1,rellim=1e-5 -selvar,pres,pres_msl,pres_sfc,theta_v,temp,t_2m,td_2m \
    cdo_tim_mean_raw_DOM${dom}.nc cdo_tim_mean_hr_mean_DOM${dom}.nc)"
  (( $? == 0 )) || {
    result+=( "DOM${dom}: pressure & temp: mean(raw) != mean(mvstream hourly mean) -- relative tolerance 1e-5" )
    diffs+=( "$diff" )
  }
  echo "DOM$dom: checking remaining variables..."
  diff="$($cdo diffn,names=left,abslim=1e-5,rellim=1 -delvar,pres,pres_msl,pres_sfc,theta_v,temp,t_2m,td_2m \
    cdo_tim_mean_raw_DOM${dom}.nc cdo_tim_mean_hr_mean_DOM${dom}.nc)"
  (( $? == 0 )) || {
    result+=( "DOM${dom}: not pressure or temp: mean(raw) != mean(mvstream hourly mean) -- absolute tolerance 1e-5" )
    diffs+=( "$diff" )
  }
  set -e
done

if (( ${#result[@]} == 0 )); then
  echo "OK!"
else
  echo "Differences detected:"
  for ((i=0; i < ${#result[@]}; i++)); do
    echo "- ${result[i]}"
    echo "${diffs[i]}"
    echo ''
  done
fi

#-- exit with exit code from mvstream test

EXIT_STATUS=${#result[@]}
if [ "$EXIT_STATUS" -ne "0" ]; then
  echo "mvstream test not passed: $EXIT_STATUS"
  exit $EXIT_STATUS
fi