Atmospheric Chemistry

In this article it is described how to perform different kinds of atmospheric chemistry simulations. This includes the description of simulations with a simplified chemistry and MECCA-based (full) chemistry, their nameless settings, possible modules to make use of and information about initialization data. Further, there are given some examples of typical simulation you can do with ICON-ART including atmospheric chemistry.

Simplified Chemistry

When we talk about simplified calculated chemistry in ICON-ART, we mean that the concentration of the gases we want to simulate is calculated with a parametrization. Here production and depletion rates are calculated to solve the differential equation

${\frac {\mathrm {dc_{i}} }{/}}{\mathrm {dt} }=P_{i}-{\frac {c_{i}}{/}}{\tau _{i}}$

numerically. Here, $c_{i}$ describes the number concentration of a certain tracer, $c_{i}$ describes the chemical production and $\tau _{i}$is the belonging life time of tracer $i$. For the namelist settings you are able to use for atmospheric chemistry, check out the ART-namelist parameters (see ART namelists). The procedure of creating an ICON-ART simulation in Atmospheric Chemistry always comes back to switching on a namelist parameter and providing the path of the respective XML-file. How to create these for several cases, please check the examples below in the Configurations part.

To learn more about technical details of simplified chemistry, see also Weimer et. al. (2017)

Note: When enabling simplified chemistry with the switch lart_chemtracer = .TRUE., you can improve your runtime but the simulated concentration values are less exact compared to MECCA-based chemistry.

MECCA-based Chemistry

General Information

The MECCA(=Module Efficiently Calculating the Chemistry of the Atmosphere) based chemistry describes a full gas phase chemistry that can be applied as an extension to the parametrized Simplified Chemistry (see above). MECCA based chemistry is generally more exact in the concentration values but the overall runtime is longer compared to purely simplified chemistry simulations. MECCA itself is originally a submodule of the CAABA box model where an air parcel is described as a box and outgoing from this model all exchange processes in- and outward of the box are calculated. As MECCA is part of this model, it contains a wide collection of the most important reactions, including Ozone-, Methane-, HOx-, NOx-, Carbonhydrogen-, Halogene- and Sulfur chemistry. MECCA is available in a supplement, available to download for free and containing all auxiliaries to perform MECCA-simulations.

Including MECCA-based Chemistry in a ICON-ART Simulation

(Note: It is recommended to perform all the following steps in the shell environment.)

The above mentioned collection of the gase phase chemistry reactions can be found in the supplement in the gas.eqn (path: caaba3.0/Mecca/gas.eqn). Additionally it is also possible to edit existing reactions as well as creating new reactions with the help of "Replacement-files" (see an example in the Configurations part). Inside the gas.eqn every reaction is marked with a certain code. To select the specific reactions for the machanism labels can be set to your belonging reactions or, more easily, a new Gas-Equation-file gas_Mechanism1.eqn can be created, containing only the wanted reactions. (Note: Never edit the original gas.eqn! Better copy it in the first place and then rename and edit it, depending on the respective scientific goal.) After that the following steps have to be fulfilled to create the code of your specific mechanism and to be able to execute an ICON-ART simulation with MECCA-based chemistry:

• set up a batch file: all previously set information about the mechanism can be selected and stated here (an example can be found below or also inside the supplement in /caaba3.0/mecca/batch/example.bat).
• execute ./mecca inside the folder /caaba3.0/mecca. Here the previously created batch file has to be selected and the Fortran files with the mechanism are created.
• since the created Fortran code is only located inside Mecca and not in ICON-ART so far, a transfer has to be carried out. A script that performs this transfer can be obtained via git clone https://gitlab.dkrz.de/art/mecca preproc.git.
• in a new directory Mecca_preproc has been generated and the script create_icon_code4.sh can be found inside of it. By executing /.create_icon_code4.sh -h paths to the Mecca- and ICON home directories can be provided as well as a name for the XML-file that is going to be linked in the unscript later.
• the Mecca-XML-file is now generated and can be found in ICON in /icon home>/runctrl examples/xml ctrl.

Now, in the respective runscript the namelist parameter lart_mecca has be set to .TRUE and for cart_mecca_xml the path to the Mecca file can be provided. Important: As a final step, the ICON code has to be recompiled with the command ./config/dkrz/levante.intel --enable-art --enable-ecrad and after executed make -j 8.

Configurations

Soon some examples for Atmospheric Chemistry simulations will be shown here.

- work in progress -

Example: How to create a simulation with simplified chemistry?

In this first example it is shown how to perform a simulation with simplified chemistry in ICON-ART. Emission data will be applied on the simulation as well. 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:

${\ce {O3 + hv -> O2 + O(^1D)}}$ ${\ce {O(^1D) + H2O -> 2OH}}$

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

${\ce {O3 + hv -> O2 + O(^1D)}}$ ${\ce {O(^1D) + H2O -> 2OH}}$

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

${\ce {OH + CH4 -> H2O + CH3 -> ... -> CO + HO2}}$ ${\ce {OH + CO -> H + CO2}}$ ${\ce {OH + CO -> HOCO}}$

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 R2B04-grid the following emission data are the most suitable ones for the respective trace gases:

• ${\ce {CH4}}$: anthropogenic (EDGAR-432 monthly), biomass-burning (GFED3), biogenic (MEGAN-MACC)
• ${\ce {CO}}$: 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.

Let's start with the runscript that has to be prepared. Here it is named ohsim_simple_icon.run but you can call it differently as well. Inside of that, 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.

#!/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.

! 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 ICON-namelist parameters have to be set. For a regular ICON-ART-Simulation the following settings are recommended - if not stated differently.

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

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 wee 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 have to be fulfilled.

&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

Example: How to create a simulation with MECCA-based chemistry?

-Work in progress-