(using icon-2.3.00, following the RCEMIP protocol)

The GridGenerator (git@git.mpimet.mpg.de:GridGenerator.git) provides a script template /run/grid_creators/grid.create_torus_grid.

cp /run/grid_creators/grid.create_torus_grid into /run

make_runscript grid.create_torus_grid

Set appropriate values in the last line for number of rows and columns as well the edge length in meter, e.g.:

y_noOfRows = 22
x_noOfColumns = 20
edge_length = 5000 m

a) Download “profile_analytical_o3.nc”:https://code.mpimet.mpg.de/attachments/19451/profile_analytical_o3.nc.
This file contains the 1-dimensional vertical profile of ozone concentration (in ppmv), as prescribed by the RCEMIP protocol.

b) expand to the full grid

cdo enlarge,Torus_Triangles_20x22_5000m.nc \
	profile_analytical_o3.nc \
	rcemip_analytical_o3_20x22_5000m.nc

When processed with newer versions of cdo, it may be required to rename cell to ncells for the NWP Setup:

 
 ncrename -O -d cell,ncells rcemip_analytical_o3_20x22_5000m.nc \
        rcemip_analytical_o3_20x22_5000m_nwp_ncells.nc
 mv rcemip_analytical_o3_20x22_5000m_ncells.nc rcemip_analytical_o3_20x22_5000m.nc
        

c) adapt the run script: set irad_o3 = 4 in radiation_nml and link the file

 add_link_file ${OZONEDIR}/rcemip_analytical_o3_20x22_5000m.nc o3_icon_DOM01.nc

To have the same initial thermodynamic profile as defined by the RCEMIP protocol, two new testcases need to be defined in the ICON source code, 'RCEMIP_analytical' and 'RCEMIP_ncdf'. 'RCEMIP_analytical' initializes the base state using the analytical profiles of T,q defined by the RCEMIP protocol, and 'RCEMIP_ncdf' reads in mean profiles of T,q from a netcdf file of a previous simulation, assuming an identical grid.

The new testcases are defined in src/testcases/mo_nh_testcases.f90. There the new subroutines 'init_torus_rcemip_analytical_sounding' and 'init_torus_rcemip_ncdf_sounding' are called, contained in src/testcases/mo_nh_torus_exp.f90.

In total, the following src files have been adapted:
testcases/mo_nh_testcases.f90
testcases/mo_nh_torus_exp.f90

The solar constant, the zenith angle and the surface albedo need to be adapted. In the source code, the solar constant needs to be factorized into the individual bands for our newly defined testcases (in mo_nwp_phy_init.f90). Besides, the albedo needs to be defined as a namelist parameter (in mo_nh_testcases_nml.f90), and needs to be transferred to an ICON variable for the two testcases, replacing the albedos calculated as a function of soiltype and zenith angle (in mo_nwp_rad_interface.f90).

In total, the following src files have been adapted:
atm_phy_nwp/mo_nwp_phy_init.f90
namelists/mo_nh_testcases_nml.f90
atm_phy_nwp/mo_nwp_rad_interface.f90

In the run script, the nh_testcase_nml namelist then can be set like this (following the RCEMIP protocol):

 &nh_testcase_nml
  nh_test_name   = 'RCEMIP_analytical'     ! test case identifier
  ape_sst_case   = 'sst_const'
  w_perturb      = 0.1        ! m/s, default=0.05
  th_perturb     = 0.3        ! K, default=0.2
  zenithang      = 42.05 ! degrees
  sol_const      = 551.58
  albedo_set     = 0.07 ! fixed surface albedo; default is 0.07
 /

In the radiation_nml, izenith = 5 needs to be set to have no diurnal cycle and the aerosol/gas mixing ratios need to be set according to the RCEMIP protocol.

 &radiation_nml
  irad_o3    = 4 ! read in external ozone file
  izenith    = 5 ! 5: prescribed zenith angle (no diurnal cycle); 6: equinox diurnal cycle
  irad_aero  = 0 ! switch off aerosols
  irad_cfc11 = 0 
  irad_cfc12 = 0
  irad_co2   = 2
  irad_ch4   = 2
  irad_n2o   = 2
  vmr_co2    =  348.e-06 ! RCEMIP values; modern-day values
  vmr_ch4    = 1650.e-09
  vmr_n2o    =  306.e-09
 /   

To have a vertical grid similar to the RCEMIP protocol, num_lev = 75 needs to be set in run_nml, and the sleve_nml needs to be adapted:

 &sleve_nml
  top_height = 48000.0 ! these settings yield an actual model top of about 34 km
  min_lay_thckn = 75. ! thickness of lowest model layer
  max_lay_thckn   = 500.    ! max thickness
  htop_thcknlimit = 40000.0 ! dont exceed max_lay_thckn below this height
  stretch_fac     = 1.2     ! values > 1 increase dz near model top
 /

In the dynamics_nml namelist, lcoriolis = .FALSE. needs to be set.

In the les_nml namelist, isrfc_type = 5 (fixed SST) and, e.g., sst = 300 need to be set. In this namelist, a threshold value for the minimum surface wind speed can also be set, which is asked for by the RCEMIP protocol (min_sfc_wind = 1.0).

A complete run script as it has been used for one of the RCEMIP simulations (ICON LEM, small domain, 300K) can be downloaded here.

To switch from ICON LEM to ICON NWP, three modifications to the run script are necessary: In nwp_phy_nml, the Raschendorfer turbulence parameterization (inwp_turb = 1) and the diagnostic cloud cover parameterization by Martin Koehler (inwp_cldcover = 1) should be used, and in nh_testcase_nml, the fixed value for the SST has to be set as ape_sst_val. Note that ape_sst_val has to be set in degrees Celsius!

WARNING: ape_sst_val is only used by the model if the new testcases are included in atm_phy_nwp/mo_nwp_phy_init.f90:

   IF (ltestcase .AND. (nh_test_name == 'APE_nwp' .OR. nh_test_name == 'dcmip_tc_52' .OR. nh_test_name == 'RCEMIP_analytical'  .OR. nh_test_name == 'RCEMIP_ncdf') ) THEN  !Tobias Becker added new testcases for RCEMIP
     ! t_g = ape_sst1
     DO jc = i_startidx, i_endidx
       zlat = p_patch%cells%center(jc,jb)%lat
       p_prog_lnd_now%t_g  (jc,jb)   = ape_sst(ape_sst_case,zlat) ! set SST
       p_prog_lnd_new%t_g  (jc,jb)   = ape_sst(ape_sst_case,zlat)
       p_prog_lnd_now%t_g_t(jc,jb,1) = ape_sst(ape_sst_case,zlat)
       p_prog_lnd_new%t_g_t(jc,jb,1) = ape_sst(ape_sst_case,zlat)
       ! Humidity at water surface = humidity at saturation
       p_diag_lnd%qv_s(jc,jb)     = &
         &  spec_humi(sat_pres_water(p_prog_lnd_now%t_g(jc,jb)),p_diag%pres_sfc(jc,jb))
       p_diag_lnd%qv_s_t(jc,jb,1) = p_diag_lnd%qv_s(jc,jb)
     END DO

The longitudes and latitudes on the Torus grid must not be interpreted. Instead we have to remap the data onto a regular grid ranging from 180°W to 180°E and 10°S to 10°N. Remapped data can then be displayed in the index space, not in the geographical domain. The grid description of the regular grid is stored in grid_icon_torus_nx20ny22_5km:

gridtype = lonlat
xsize    =   20
ysize    =   20
xfirst   = -180.0  ! start longitude (fixed)
xinc     =   18    ! 360 / xsize
yfirst   =  -10.0  ! start latitude (fixed)
yinc     =    1    !  20 / ysize

Generate a weight file once:

 cdo gendis,grid_icon_torus_nx20ny22_5km <infile> remapweights_nx20ny22_5km.nc

Remap the data using the weight file:

 cdo remap,grid_icon_torus_nx20ny22_5km,remapweights_nx20ny22_5km.nc <infile> <outfile>