Simplified Chemistry

From icon-art guide
Revision as of 12:10, 14 July 2023 by Tim R (talk | contribs) (Update on Runscript settings)
Jump to navigation Jump to search

- work in progress -

In this first example it is shown how to perform a simulation of with simplified chemistry in ICON-ART. This tutorial teaches you...

  • the basics of setting up an ICON runscript with ICON-ART settings
  • the use of the most simple ICON-ART namelist parameter
  • the implementation of the desired chemical species in a simulation by setting up a chemtracer xml-data for simplified chemistry simulations
  • the implementation of emission data in a simulation

Emission data will be applied on the simulation as well.

Configuration case

The depicted case is about simulating the tropospheric hydroxyl radical (OH), one of the most important oxidants of the atmosphere. It's main source in the lower troposphere is the photolysis of ozone and its consequent reaction of an excited oxygen atom with the surrounding water vapor:

Additionally the excited Oxygen atom reacts further with Nitrogen and Oxygen:

The main sink of OH in the Troposphere is methane and carbon monooxide:

Now the OH concentrations are calculated with the respective kinetic and photolysis constants, based on chemical kinetic laws: Failed to parse (syntax error): {\displaystyle [\ce{OH}]=\frac{\mathrm{2[\ce{O(^1D)}]}k_{\ce{H2O}}[\ce{H2O}]}{k_{\ce{CH4}}[\ce{CH4}]+(k_{\ce{CO,1}}+k_{\ce{CO,2}})[\ce{CO}]}} with Failed to parse (syntax error): {\displaystyle [\ce{O(^1D)}]=\frac{J_{\ce{O3}}[\ce{O3}]}{k_{\ce{O2}}[\ce{O2}]+k_{\ce{N2}}[\ce{N2}]+k_{\ce{H2O}}[\ce{H2O}]}} .

Additionally emission data of the main sinks of OH are implemented. Since the simulation is performed on a R2B05-grid the following emission data are the most suitable ones for the respective trace gases:

  • : anthropogenic (EDGAR-432 monthly), biomass-burning (GFED3), biogenic (MEGAN-MACC)
  • : anthropogenic (EDGAR-432 monthly)

For more information on recommended emission data see the abstract, dealing with Emission Data. Since emission data are relatively large, they can also just be left out.

Setting up the Runscript

Let's start with the runscript that has to be prepared. Please note that the in the following explained parts have to be printed in one runscript-file with the naming designation "xyz.run". Here it is named ohsim_simple_icon.run but of course you can call it differently as well.

Inside of that, first check that all your directories are correct, probably they have to be adjusted. Abbreviations used here are the following:

  • CENTER: Your organization
  • EXPNAME: name of your ICON-Simulation
  • OUTDIR: Directory where the simulation output will be stored
  • ARTFOLDER: Directory where the ICON-ART code is stored
  • INDIR: Directory where the necessary Input data are stored
  • EXP:
  • lart: For ICON-ART Simulation that has to be switched to Oheim_simple_icon.run.

Part 1: Runscript Directory Settings (Example configuration) 

#!/bin/bash
CENTER=IMK
workspace=/hkfs/work/workspace/scratch/hp8526-ws_icon_oh
basedir=${workspace}/icon-kit-testsuite
icon_data_poolFolder=/hkfs/work/workspace/scratch/fb4738-dwd_ozone/icon/INPUT/AMIP/amip_input
EXPNAME=ohsim_icon_simple_atom1
OUTDIR=${workspace}/output/${EXPNAME}
ICONFOLDER=/home/hk-project-iconart/hp8526/icon-kit
ARTFOLDER=${ICONFOLDER}/externals/art
INDIR=/hkfs/work/workspace/scratch/fb4738-dwd_ozone/icon/INPUT
EXP=ECHAM_AMIP_LIFETIME
lart=.True.

FILETYPE=4
COMPILER=intel
restart=.False.
read_restart_namelists=.False.

# Remove folder  from OUTDIR for postprocessing output
OUTDIR_PREFIX=${workspace}

# Create output directory and go to this directory

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

cd $OUTDIR

#input for global domain
ln -sf ${INDIR}/../INPUT/GRID/icon_grid_0014_R02B05_G.nc iconR2B05_DOM01.nc
ln -sf ${INDIR}/../INPUT/EXTPAR/icon_extpar_0014_R02B05_G.nc extpar_iconR2B05_DOM01.nc
ln -sf /hkfs/work/workspace/scratch/fb4738-dwd_ozone/icon/output/0014_R02B05/uc1_ifs_t1279_grb2_remap_rev832_0014_R02B05_2018010100.nc ifs2icon_R2B05_DOM01.nc

ln -sf $ICONFOLDER/data/rrtmg_lw.nc rrtmg_lw.nc
ln -sf $ICONFOLDER/data/ECHAM6_CldOptProps.nc ECHAM6_CldOptProps.nc

ln -sf ${ARTFOLDER}/runctrl_examples/init_ctrl/mozart_coord.nc ${OUTDIR}/mozart_coord.nc
ln -sf ${ARTFOLDER}/runctrl_examples/init_ctrl/Linoz2004Br.dat ${OUTDIR}/Linoz2004Br.dat
ln -sf ${ARTFOLDER}/runctrl_examples/init_ctrl/Simnoy2002.dat ${OUTDIR}/Simnoy2002.dat

ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_scat-aer.dat      ${OUTDIR}/FJX_scat-aer.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_j2j.dat           ${OUTDIR}/FJX_j2j.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_scat-cld.dat      ${OUTDIR}/FJX_scat-cld.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_scat-ssa.dat      ${OUTDIR}/FJX_scat-ssa.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_scat-UMa.dat      ${OUTDIR}/FJX_scat-UMa.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_spec_extended.dat ${OUTDIR}/FJX_spec_extended.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_spec.dat          ${OUTDIR}/FJX_spec.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/atmos_std.dat         ${OUTDIR}/atmos_std.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/atmos_h2och4.dat      ${OUTDIR}/atmos_h2och4.dat

ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_j2j_extended.dat          ${OUTDIR}/FJX_j2j_extended.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_spec_extended.dat ${OUTDIR}/FJX_spec_extended.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_spec_extended_lyman.dat ${OUTDIR}/FJX_spec_extended_lyman.dat

# this if condition is necessary because otherwise
# a new link in ${INDIR}/${EXP}/emiss_minimal is generated
# linking to itself
if [ ! -L ${OUTDIR}/emissions ]; then
ln -sd /home/hk-project-iconart/hp8526/emissions       ${OUTDIR}/emissions
fi

# this if condition is necessary because otherwise
# a new link in ${INDIR}/${EXP}/emiss_minimal is generated
# linking to itself
if [ ! -L ${OUTDIR}/emissions ]; then
ln -sd /home/hk-project-iconart/hp8526/emissions       ${OUTDIR}/emissions
fi

# the namelist filename
atmo_namelist=NAMELIST_${EXPNAME}
#
#-----------------------------------------------------------------------------

Additionally in the next lines of code you set the timing. In this simulation we only simulate a few days. Because OH is dependent from the solar radiation, the output interval is set to 10 hours to calculate OH to every time of the day.

Part 2: Runscript Timing Settings (Example configuration) 

! run_nml: general switches ----------
&# model timing

start_date=${start_date:="2016-07-29T00:00:00Z"}
end_date=${end_date:="2016-08-23T00:00:00Z"}
output_start=${start_date:="2016-07-29T00:00:00Z"}
output_end=${end_date:="2016-08-23T00:00:00Z"}
output_interval="PT10H"
modelTimeStep="PT6M"
leadtime="P10D"
checkpoint_interval="P30D"

Further, all the namelist parameters (from the regular ICON model without ART-extension) have to be set. For a regular ICON-ART-Simulation the following settings are recommended - if not stated differently.

Part 3: Runscript ICON-Parameter and -Namelist Settings (Example configuration) 

# 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
declare -a atmo_dyn_grids=("iconR2B04_DOM01.nc" "iconR2B05_DOM02.nc" "iconR2B06_DOM03.nc")
# "iconR2B08_DOM03.nc"
#atmo_dyn_grids="iconR2B07_DOM02.nc","iconR2B08_DOM03.nc"
#atmo_rad_grids="iconR2B06_DOM01.nc"
#-----------------------------------------------------------------------------

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

no_of_models=1

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
/
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   = 0   ! up to one PE per output stream is possible
 itype_comm     = 1
 iorder_sendrecv = 3  ! best value for CRAY (slightly faster than option 1)
/
&grid_nml
 dynamics_grid_filename  = 'iconR2B05_DOM01.nc'
 dynamics_parent_grid_id = 0
!radiation_grid_filename = ${radiation_grid_filename}
 lredgrid_phys           = .false.
 lfeedback               = .true.
 ifeedback_type          = 2
/
&initicon_nml
 lconsistency_checks      = .false.
 init_mode   = 2           ! operation mode 2: IFS
 zpbl1       = 500. 
 zpbl2       = 1000. 
! l_sst_in    = .true.
/
&run_nml
 num_lev        = 90
 lvert_nest     = .true.       ! use vertical nesting if a nest is active
! nsteps         = ${nsteps}    ! 50 ! 1200 ! 7200 !
! dtime          = ${dtime}     ! timestep in seconds
 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   = 10            ! can be increased up to 10 for detailed timer output
 output         = "nml"
 lart           = ${lart}
/
&nwp_phy_nml
 inwp_gscp       = 1
 inwp_convection = 1
 inwp_radiation  = 1
 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   = 1
 icpl_aero_gscp   = 1
 ! resolution-dependent settings - please choose the appropriate one
 dt_rad    = 2160.
 dt_conv   = 720.
 dt_sso    = 1440.
 dt_gwd    = 1440.
/
&nwp_tuning_nml
 tune_zceff_min = 0.075 ! ** 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
 rat_sea = 8.5  ! ** 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      = 3      !!! 1 for assimilation cycle and forecast
 lmulti_snow    = .true. !!! .false. for assimilation cycle and forecast
 itype_heatcond = 2
 idiag_snowfrac = 2
 lsnowtile      = .false.  !! later on .true. if GRIB encoding issues are solved
 lseaice        = .true.
 llake          = .false.
 itype_lndtbl   = 3  ! minimizes moist/cold bias in lower tropical troposphere
 itype_root     = 2
/
&radiation_nml
 irad_o3       = 10
 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
/
&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
 lhdiff_rcf     = .true.
 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
 l_open_ubc     = .false.
 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
 idiv_method    = 1
 divavg_cntrwgt = 0.50
 lcoriolis      = .TRUE.
/
&transport_nml
!                qv, qc, qi, qr, qs
 itype_vlimit  = 1,1,1,1,1
 ivadv_tracer = 3, 3, 3, 3, 3
 itype_hlimit = 3, 4, 4, 4 , 4
 ihadv_tracer = 52, 2,2,2,2
 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.
/
&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-${EXPNAME}-chem'                ! file name base
 ml_varlist                   = 'temp','pres','group:ART_CHEMISTRY','u','v','OH_Nconc'
 output_grid                  =  .TRUE.
 remap                        = 1
 reg_lon_def                  = -180.,1,180.
 reg_lat_def                  = -90.,1,90.
/

Please note in the last namelist section "output_nml" that you can set all output variables that you need to postprocess your data later. All assigned variables here will be written in the output netCDF-files as well. To learn more about post processing your data, check out a later chapter of this article or the Postprocessing article.

Now, we're getting to the ICON-ART settings. To enable chemistry in an ICON-ART Simulation inn general, the switch lart_chem has to be set to .TRUE.. With lart_diag_out output of the diagnostic fields can be enabled. Due to setting lart_chem=.TRUE. either lart_chemtracer or lart_mecca have to be set to .TRUE.. Because we want to perform a simulation with simplified chemistry, we have to switch on lart_chemtracer. If this namelist parameter is set to .TRUE., also cart_chemtracer_xml has to be fulfilled. Here you enter the path of your xml-file which describes the tracers occurring and their properties in the simulation. How to create this xml-file is explained in the next chapter. Because this xml-file will contain information about included emission data of certain chemical species, in cart_emiss_xml_file the path of a second emission-xml-file has to be set. An example configuration for this part is shown in the following:

Part 4: Runscript ICON-ART Settings (Example configuration) 

&art_nml
 lart_chem       = .TRUE.
 lart_diag_out   = .TRUE.
 lart_aerosol    = .FALSE.
 lart_mecca      = .FALSE.
 lart_chemtracer = .TRUE.
 cart_emiss_xml_file   = '${ARTFOLDER}/runctrl_examples/emiss_ctrl/emissions_R2B05_0014_cs.xml'
 cart_chemtracer_xml   = '/home/hk-project-iconart/hp8526/icon-kit/externals/art/runctrl_examples/xml_ctrl/chemtracer_reimus.xml'
 cart_input_folder     = '${OUTDIR}'
 cart_io_suffix        = '0014'
 iart_init_gas         =  0
/
EOF

Please note that there are also several other namelist parameter you can select from (see Namelist article) but to perform our we're done for this part.

Depending on the used HPC-System, some parameter concerning the running job like maximum running time and used nodes can be set. For this case study the following settings can be copied. Note that this is valid for the HoreKa HPC system and that it can differ to other systems.

Part 5: Runscript job Settings (Example configuration) 

cp ${ICONFOLDER}/bin/icon ./icon.exe

cat > job_ICON << ENDFILE
#!/bin/bash -x
#SBATCH --nodes=4
#SBATCH --time=12:00:00
#SBATCH --ntasks-per-node=76
#SBATCH --partition=cpuonly
#SBATCH -A hk-project-iconart
###SBATCH --constraint=LSDF

module load compiler/intel/2022.0.2 mpi/openmpi/4.0 lib/netcdf/4.9_serial lib/hdf5/1.12_serial lib/netcdf-fortran/4.5_serial lib/eccodes/2.25.0 numlib/mkl/2022.0.2

mpirun --bind-to core --map-by core --report-bindings ./icon.exe

ENDFILE

chmod +x job_ICON
sbatch job_ICON

To conclude and to double check, in the following box the complete runscript is shown once again.

Complete example configuration of the runscript 

#!/bin/bash
CENTER=IMK
workspace=/hkfs/work/workspace/scratch/hp8526-ws_icon_oh
basedir=${workspace}/icon-kit-testsuite
icon_data_poolFolder=/hkfs/work/workspace/scratch/fb4738-dwd_ozone/icon/INPUT/AMIP/amip_input
EXPNAME=ohsim_icon_mecca_atom1
OUTDIR=${workspace}/output/${EXPNAME}
ICONFOLDER=/home/hk-project-iconart/hp8526/icon-kit
ARTFOLDER=${ICONFOLDER}/externals/art
INDIR=/hkfs/work/workspace/scratch/fb4738-dwd_ozone/icon/INPUT
EXP=ECHAM_AMIP_LIFETIME
lart=.True.

FILETYPE=4
COMPILER=intel
restart=.False.
read_restart_namelists=.False.


# Remove folder  from OUTDIR for postprocessing output
OUTDIR_PREFIX=${workspace}

# Create output directory and go to this directory

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

cd $OUTDIR

#input for global domain
ln -sf ${INDIR}/../INPUT/GRID/icon_grid_0014_R02B05_G.nc iconR2B05_DOM01.nc
ln -sf ${INDIR}/../INPUT/EXTPAR/icon_extpar_0014_R02B05_G.nc extpar_iconR2B05_DOM01.nc
ln -sf /hkfs/work/workspace/scratch/fb4738-dwd_ozone/icon/output/0014_R02B05/uc1_ifs_t1279_grb2_remap_rev832_0014_R02B05_2018010100.nc ifs2icon_R2B05_DOM01.nc

ln -sf $ICONFOLDER/data/rrtmg_lw.nc rrtmg_lw.nc
ln -sf $ICONFOLDER/data/ECHAM6_CldOptProps.nc ECHAM6_CldOptProps.nc

ln -sf ${ARTFOLDER}/runctrl_examples/init_ctrl/mozart_coord.nc ${OUTDIR}/mozart_coord.nc
ln -sf ${ARTFOLDER}/runctrl_examples/init_ctrl/Linoz2004Br.dat ${OUTDIR}/Linoz2004Br.dat
ln -sf ${ARTFOLDER}/runctrl_examples/init_ctrl/Simnoy2002.dat ${OUTDIR}/Simnoy2002.dat

ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_scat-aer.dat      ${OUTDIR}/FJX_scat-aer.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_j2j.dat           ${OUTDIR}/FJX_j2j.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_scat-cld.dat      ${OUTDIR}/FJX_scat-cld.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_scat-ssa.dat      ${OUTDIR}/FJX_scat-ssa.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_scat-UMa.dat      ${OUTDIR}/FJX_scat-UMa.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_spec_extended.dat ${OUTDIR}/FJX_spec_extended.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_spec.dat          ${OUTDIR}/FJX_spec.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/atmos_std.dat         ${OUTDIR}/atmos_std.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/atmos_h2och4.dat      ${OUTDIR}/atmos_h2och4.dat

ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_j2j_extended.dat          ${OUTDIR}/FJX_j2j_extended.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_spec_extended.dat ${OUTDIR}/FJX_spec_extended.dat
ln -sf ${ARTFOLDER}/runctrl_examples/photo_ctrl/FJX_spec_extended_lyman.dat ${OUTDIR}/FJX_spec_extended_lyman.dat


# this if condition is necessary because otherwise
# a new link in ${INDIR}/${EXP}/emiss_minimal is generated
# linking to itself
if [ ! -L ${OUTDIR}/emissions ]; then
ln -sd /home/hk-project-iconart/hp8526/emissions       ${OUTDIR}/emissions
fi


# the namelist filename
atmo_namelist=NAMELIST_${EXPNAME}
#
#-----------------------------------------------------------------------------
# global timing
ndays_restart=366
dt_restart=`expr ${ndays_restart} \* 86400`
#
#-----------------------------------------------------------------------------
# model timing

start_date=${start_date:="2016-07-29T00:00:00Z"}
end_date=${end_date:="2016-08-23T00:00:00Z"}
output_start=${start_date:="2016-07-29T00:00:00Z"}
output_end=${end_date:="2016-08-23T00:00:00Z"}
output_interval="PT10H"
modelTimeStep="PT6M"
leadtime="P10D"
checkpoint_interval="P30D"

#
#-----------------------------------------------------------------------------
# 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
declare -a atmo_dyn_grids=("iconR2B04_DOM01.nc" "iconR2B05_DOM02.nc" "iconR2B06_DOM03.nc")
# "iconR2B08_DOM03.nc"
#atmo_dyn_grids="iconR2B07_DOM02.nc","iconR2B08_DOM03.nc"
#atmo_rad_grids="iconR2B06_DOM01.nc"
#-----------------------------------------------------------------------------

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

no_of_models=1

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
/
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   = 0   ! up to one PE per output stream is possible
 itype_comm     = 1
 iorder_sendrecv = 3  ! best value for CRAY (slightly faster than option 1)
/
&grid_nml
 dynamics_grid_filename  = 'iconR2B05_DOM01.nc'
 dynamics_parent_grid_id = 0
!radiation_grid_filename = ${radiation_grid_filename}
 lredgrid_phys           = .false.
 lfeedback               = .true.
 ifeedback_type          = 2
/
&initicon_nml
 lconsistency_checks      = .false.
 init_mode   = 2           ! operation mode 2: IFS
 zpbl1       = 500. 
 zpbl2       = 1000. 
! l_sst_in    = .true.
/
&run_nml
 num_lev        = 90
 lvert_nest     = .true.       ! use vertical nesting if a nest is active
! nsteps         = ${nsteps}    ! 50 ! 1200 ! 7200 !
! dtime          = ${dtime}     ! timestep in seconds
 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   = 10            ! can be increased up to 10 for detailed timer output
 output         = "nml"
 lart           = ${lart}
/
&nwp_phy_nml
 inwp_gscp       = 1
 inwp_convection = 1
 inwp_radiation  = 1
 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   = 1
 icpl_aero_gscp   = 1
 ! resolution-dependent settings - please choose the appropriate one
 dt_rad    = 2160.
 dt_conv   = 720.
 dt_sso    = 1440.
 dt_gwd    = 1440.
/
&nwp_tuning_nml
 tune_zceff_min = 0.075 ! ** 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
 rat_sea = 8.5  ! ** 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      = 3      !!! 1 for assimilation cycle and forecast
 lmulti_snow    = .true. !!! .false. for assimilation cycle and forecast
 itype_heatcond = 2
 idiag_snowfrac = 2
 lsnowtile      = .false.  !! later on .true. if GRIB encoding issues are solved
 lseaice        = .true.
 llake          = .false.
 itype_lndtbl   = 3  ! minimizes moist/cold bias in lower tropical troposphere
 itype_root     = 2
/
&radiation_nml
 irad_o3       = 10
 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
/
&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
 lhdiff_rcf     = .true.
 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
 l_open_ubc     = .false.
 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
 idiv_method    = 1
 divavg_cntrwgt = 0.50
 lcoriolis      = .TRUE.
/
&transport_nml
!                qv, qc, qi, qr, qs
 itype_vlimit  = 1,1,1,1,1
 ivadv_tracer = 3, 3, 3, 3, 3
 itype_hlimit = 3, 4, 4, 4 , 4
 ihadv_tracer = 52, 2,2,2,2
 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.
/
&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-${EXPNAME}-chem'                ! file name base
 ml_varlist                   = 'z_ifc','temp','pres','group:ART_CHEMISTRY','u','v','OH_Nconc'
 output_grid                  =  .TRUE.
 remap                        = 1
 reg_lon_def                  = -180.,1,180.
 reg_lat_def                  = -90.,1,90.
/
&art_nml
 lart_diag_out   = .TRUE.
 lart_aerosol    = .FALSE.
 lart_chem       = .TRUE.
 lart_mecca      = .TRUE.
 lart_chemtracer = .TRUE.
 cart_emiss_xml_file   = '${ARTFOLDER}/runctrl_examples/emiss_ctrl/emissions_R2B05_0014_cs.xml'
 cart_mecca_xml        = '${ARTFOLDER}/runctrl_examples/xml_ctrl/tracers_oh_reimus.xml'
 cart_chemtracer_xml   = '/home/hk-project-iconart/hp8526/icon-kit/externals/art/runctrl_examples/xml_ctrl/chemtracer_reimus.xml'
 cart_input_folder     = '${OUTDIR}'
 cart_io_suffix        = '0014'
 iart_init_gas         =  0
/
EOF


	
cp ${ICONFOLDER}/bin/icon ./icon.exe
	   
cat > job_ICON << ENDFILE
#!/bin/bash -x
#SBATCH --nodes=4
#SBATCH --time=12:00:00
#SBATCH --ntasks-per-node=76
#SBATCH --partition=cpuonly
#SBATCH -A hk-project-iconart
###SBATCH --constraint=LSDF

module load compiler/intel/2022.0.2 mpi/openmpi/4.0 lib/netcdf/4.9_serial lib/hdf5/1.12_serial lib/netcdf-fortran/4.5_serial lib/eccodes/2.25.0 numlib/mkl/2022.0.2

mpirun --bind-to core --map-by core --report-bindings ./icon.exe

ENDFILE

chmod +x job_ICON
sbatch job_ICON

Setting up the xml-files

An xml-file describes the chemical components of the simulation which means that all trace gases or aerosols and their properties that are relevant for the simulation are listed here. Since we perform a simulation with simplified ICON-ART chemistry we need the matching chemtracer-xml-file. Additionally we need to create an emission-xml-file since emission data of certain chemical species is included in our simulation.

Chemtracer-xml-file

This file contains all the necessary information to describe the chemical mechanism with the respective important chemical species. The following information are given per chemical species:

  • properties like mol weight, units or lifetime of the tracer
  • main sink
  • main reaction product after reacting with the main sink after a given lifetime
  • partly information about included emission data (ANT=anthropogenic emission, BIO=biogenic emission, BBE=biomass burning emission, more information about available emission data, see Emission Data)

Part 4: Chemtracer-xml-file for OH chemistry (Example configuration) 

<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE tracers SYSTEM "tracers.dtd">

<tracers>
  <chemtracer id="TRCH4" full="FALSE" chemtr="TRUE">
    <tag001 type="char">chemtr</tag001>
    <mol_weight type="real">1.604E-2</mol_weight>
    <?source_lifetime Hayman et al., ACP, 2017 ?>
    <lifetime type="real">286977600</lifetime>
    <transport type="char"> stdchem </transport>
    <unit type="char">mol mol-1</unit>
    <c_solve type="char">OH</c_solve>
    <products type="char">TRCO</products>
    <emiss_ANT type="char" inum_levs="1">CH4_ANT_EDGAR432-monthly</emiss_ANT>
    <emiss_BBE type="char" inum_levs="1">CH4_BBE_GFED3</emiss_BBE>
    <emiss_BIO type="char" inum_levs="1">CH4_BIO_MEGAN-MACC</emiss_BIO>
    <init_mode type="int">0</init_mode>
    <init_name type="char">CH4</init_name>
  </chemtracer>
  <chemtracer id="TRC2H6" full="FALSE" chemtr="TRUE">
    <tag001 type="char">chemtr</tag001>
    <mol_weight type="real">3.006E-2</mol_weight>
    <?source_lifetime Hodnebrog et al, Atmos. Sci. Lett., 2018 ?>
    <lifetime type="real">5011200</lifetime>
    <transport type="char"> stdchem </transport>
    <unit type="char">mol mol-1</unit>
    <c_solve type="char">OH</c_solve>
    <init_mode type="int">0</init_mode>
    <init_name type="char">C2H6</init_name>
  </chemtracer>
  <chemtracer id="TRC3H8" full="FALSE" chemtr="TRUE">
    <tag001 type="char">chemtr</tag001>
    <mol_weight type="real">4.40956E-2</mol_weight>
    <?source_lifetime Hodnebrog et al, Atmos. Sci. Lett., 2018 ?>
    <lifetime type="real">1123200</lifetime>
    <transport type="char"> stdchem </transport>
    <unit type="char">mol mol-1</unit>
    <c_solve type="char">OH</c_solve>
    <products type="char">0.736*TRCH3COCH3</products>
    <init_mode type="int">0</init_mode>
    <init_name type="char">C3H8</init_name>
  </chemtracer>
  <chemtracer id="TRC5H8"  full="FALSE" chemtr="TRUE">
    <tag001 type="char">chemtr</tag001>
    <mol_weight type="real">6.812E-2</mol_weight>
    <?source_lifetime Weimer (2015), p. 16 ?>
    <lifetime type="real">8640</lifetime>
    <transport type="char"> stdchem </transport>
    <unit type="char">mol mol-1</unit>
    <init_mode type="int">0</init_mode>
    <c_solve type="char">OH</c_solve>
  </chemtracer>
  <chemtracer id="TRCH3COCH3" full="FALSE" chemtr="TRUE">
    <tag001 type="char">chemtr</tag001>
    <mol_weight type="real">5.808E-2</mol_weight>
    <?source_lifetime Weimer (2015), p. 9 ?>
    <lifetime type="real">1728000</lifetime>
    <transport type="char"> stdchem </transport>
    <unit type="char">mol mol-1</unit>
    <init_mode type="int">0</init_mode>
    <c_solve type="char">OH</c_solve>
  </chemtracer>
  <chemtracer id="TRCO" full="FALSE" chemtr="TRUE">
    <tag001 type="char">chemtr</tag001>
    <mol_weight type="real">2.801E-2</mol_weight>
    <?source_lifetime Ehhalt et al., IPCC, 2001, Chapter 4 ?>
    <lifetime type="real">5184000</lifetime>
    <transport type="char"> stdchem </transport>
    <unit type="char">mol mol-1</unit>
    <c_solve type="char">OH</c_solve>
    <products type="char">TRCO2</products>
    <emiss_ANT type="char" inum_levs="1">CO_ANT_EDGAR432-monthly</emiss_ANT>
    <init_mode type="int">0</init_mode>
    <init_name type="char">CO</init_name>
  </chemtracer>
  <chemtracer id="TRCO2" full="FALSE" chemtr="TRUE">
    <tag001 type="char">chemtr</tag001>
    <mol_weight type="real">4.401E-2</mol_weight>
    <?source_lifetime Houghton et al., IPCC, Cambridge University Press, 2001 ?>
    <lifetime type="real">3153600000</lifetime>
    <transport type="char"> stdchem </transport>
    <unit type="char">mol mol-1</unit>
    <init_mode type="int">0</init_mode>
    <c_solve type="char">lt</c_solve>
    <products type="char">TRCO</products>
  </chemtracer>
  <chemtracer id="TRNH3" full="FALSE" chemtr="TRUE">
    <tag001 type="char">chemtr</tag001>
    <mol_weight type="real">1.70E-2</mol_weight>
    <?source_lifetime Pinder et al., GRL, 2008?>
    <lifetime type="real">86400</lifetime>
    <transport type="char"> stdchem </transport>
    <unit type="char">mol mol-1</unit>
    <init_mode type="int">0</init_mode>
    <c_solve type="char">OH</c_solve>
    <products type="char">TRNO2</products>
  </chemtracer>
  <chemtracer id="TRNO2" full="FALSE" chemtr="TRUE">
    <tag001 type="char">chemtr</tag001>
    <mol_weight type="real">4.601E-2</mol_weight>
    <?source_lifetime Wypych, 2017, Atlas of Material Damage?>
    <lifetime type="real">4730400000</lifetime>
    <transport type="char"> stdchem </transport>
    <unit type="char">mol mol-1</unit>
    <init_mode type="int">0</init_mode>
    <c_solve type="char">OH</c_solve>
    <products type="char">TRHNO3</products>
  </chemtracer>
  <chemtracer id="TRSO2" full="FALSE" chemtr="TRUE">
    <tag001 type="char">chemtr</tag001>
    <mol_weight type="real">6.40E-2</mol_weight>
    <?source_lifetime Von Glasow, Chemical Geology, 2009 ?>
    <lifetime type="real">1209600</lifetime>
    <transport type="char"> stdchem </transport>
    <unit type="char">mol mol-1</unit>
    <init_mode type="int">0</init_mode>
    <c_solve type="char">OH</c_solve>
    <products type="char">TRH2SO4</products>
  </chemtracer>
  <chemtracer id="TROCS" full="FALSE" chemtr="TRUE">
    <tag001 type="char">chemtr</tag001>
    <mol_weight type="real">6.01E-2</mol_weight>
    <?source_lifetime Ullwer (2017) ?>
    <lifetime type="real">504576000</lifetime>
    <transport type="char"> stdchem </transport>
    <unit type="char">mol mol-1</unit>
    <init_mode type="int">0</init_mode>
    <c_solve type="char">OH</c_solve>
    <products type="char">TRSO2</products>
  </chemtracer>
  <chemtracer id="TRDMS" full="FALSE" chemtr="TRUE">
    <tag001 type="char">chemtr</tag001>
    <mol_weight type="real">6.21E-2</mol_weight>
    <?source_lifetime Ullwer (2017) ?>
    <lifetime type="real">216000</lifetime>
    <transport type="char"> stdchem </transport>
    <unit type="char">mol mol-1</unit>
    <init_mode type="int">0</init_mode>
    <c_solve type="char">OH</c_solve>
    <products type="char">0.993*TRSO2;0.007*TROCS</products>
  </chemtracer>
  <chemtracer id="TRHNO3" full="FALSE" chemtr="TRUE">
    <tag001 type="char">chemtr</tag001>
    <mol_weight type="real">6.30E-2</mol_weight>
    <?source_lifetime Day et al., ACP, 2008?>
    <lifetime type="real">21600</lifetime>
    <transport type="char"> stdchem </transport>
    <init_name type="char"> HNO3 </init_name>
    <unit type="char">mol mol-1</unit>
    <init_mode type="int">0</init_mode>
    <c_solve type="char">lt</c_solve>
  </chemtracer>
  <chemtracer id="TRH2SO4" full="FALSE" chemtr="TRUE">
    <tag001 type="char">chemtr</tag001>
    <mol_weight type="real">9.80E-2</mol_weight>
    <?source_lifetime Fiedler et al., ACP, 2005?>
    <lifetime type="real">1800</lifetime>
    <transport type="char"> stdchem </transport>
    <unit type="char">mol mol-1</unit>
    <init_mode type="int">0</init_mode>
    <c_solve type="char">lt</c_solve>
  </chemtracer>
  <chemtracer id="TRO3" full="FALSE" chemtr="TRUE">
    <tag001 type="char">chemtr</tag001>
    <mol_weight type="real">4.800E-2</mol_weight>
    <?source_lifetime Ehhalt et al., IPCC, 2001, Chapter 4 ?>
    <lifetime type="real">1555200</lifetime>
    <transport type="char"> stdchem </transport>
    <init_mode type="int"> 0 </init_mode>
    <init_name type="char">O3</init_name>
    <unit type="char">mol mol-1</unit>
    <c_solve type="char">linoz</c_solve>
  </chemtracer>
</tracers>

Emission-xml-file

To link the information from the chemtracer-xml-file that we want to include emission data of a specific chemical species with the specific information about the respective emission, ICON-ART is working with, we need a second xml-file: The emission-xml-file. It is normally a standard grid-dependant xml-file you can use where nothing has to be changed. Since this Simulation is performed on a R2B05 grid, we use the emissions_R2B05_0014_cs.xml. Have a look at the article about Emission Data to check out all the available emission-xml-files for different resolutions.

Running the simulation

Double check all filled in paths and namelist - especially the ART-namelists. If every namelist parameter in the runscript is filled in correctly, the runscript has to be saved. Afterwards by typing 

/.ohsim_simple_icon.run

a job can be submitted to the respective HPC-System. Type the terminal command 

squeue

to view a list of your submitted and currently running and jobs. By changing in the output directory (which is according to our runscript /hkfs/work/workspace/scratch/hp8526-ws_icon_oh/output/ohsim_icon_simple_atom1 you can check the slurm file for possible errors and run times after your job has been run through.

In the output directory you can also find all output data for postprocessing in netCDF format.

Retrieved from "https://www.icon-art.kit.edu/userguide/index.php?title=Simplified_Chemistry&oldid=238"