====== ICON Holocene Simulations ======

This documentation helps you to:
   - run a spin up simulation forced with Holocene conditions from ECHAM runs
   - run a high-resolution nesting simulation starting from restart files from the spin up simulation

==== Basic Preparation ====

To start a ICON simulation with Holocene conditions you have to download the ICON binaries for ICON-NWP, especially //icon-2.5.0-rc-orbit//. A detailed description on how to start with ICON can be found in the [[models:pot-pourri:how_to:icon_quick_start_guide|ICON quick start guide]].

<code bash>
git clone --recursive git@git.mpimet.mpimet.de:icon-les icon-2.5.0-rc-orbit
cd icon-2.5.0-rc-orbit

patch -p1 < patch-async_latbc_finalize.txt
cd externals/jsbach
patch -p1 ../../patch-jsbach_orbit.txt
cd ../..

./configure --with-fortran=intel17 --with-openmp
./build_command
</code>

The last step is to create the runscripts. To do so, you need the name of an experimental descriptor file in ./run. In case of the holocene simulations these are:

<code bash>
exp.holocene_spinup
exp.holocene_nested
</code>

In your icon-directory execute:

<code bash>
./make_runscripts holocene_spinup
./make_runscripts holocene_nested
</code>





==== Running the Holocene spin up simulation ====

In the following it will be explained how to start the 30-year spin up simulation or in general - how to force ICON with ECHAM initial and boundary data.

To run a ICON simulation you need:

   * Initial file
   * grid file
   * External parameters
   * boundary data if your simulation is limited in its area



=== Grids and external parameters ===

The grid file(s) can be created either with the DWD online grid generator. The grids from this source were unfortunately connected to problems preforming nesting simulations. The better way is to create the grids with the grid generator you can compile on mistral. 
How to download the GridGenerator can be found in the [[models:pot-pourri:how_to:icon_quick_start_guide|ICON quick start guide]].
The setup to create the grids for the Holocene Simulations:

<hidden grid.create_nesting_grids.run>
<code bash>
set_spr_optimization_grids
base_file_name="Global_Icos_<resolution>.nc"
max_resolution=40000
create_icosahedron_grids

#creat first patch
# cut_lonlat_rectangle_grid Global_Icos_0039km.nc nest_patch.nc -14.5 4 75.5 44
cut_lonlat_rectangle_grid Global_Icos_0039km.nc nest_patch.nc 14.5 4  85   54

#=============================================================================
#  create nests
  rm holocene_*_nested.nc
  cp nest_patch.nc holocene_40km_nested.nc

cut=$lonlat_rectangle_condition
inner=$inner_cells_condition
inner_cells_depth=13

  no_of_patches=5
  patchNameList=(     holocene_40km_nested.nc     holocene_20km_nested.nc      holocene_10km_nested.nc     holocene_5km_nested.nc      holocene_2_5km_nested.nc   )
  patchIDList=(           1               2                3               4                 5    )
  patchParentIDList=(     0               1                2               3                 4    )
  patchShapeList=(      none            $inner           $inner          $cut             $inner  )
  patchCenterXList=(    none            14.5             14.5            8                 8    )
  patchCenterYList=(    none              4                 4             20                20    )
  rectangleXradiousList=(none            70.5              68             45                42    )
  rectangleYradiousList=(none             39               36.5           20                18    )
   patchOptimization=(   .false.        .false.           .false.         .false.          .false. )
  boundary_indexing_depth=12

  set_hdcp2_optimization_grids
  create_hierarchical_grids

</code>
</hidden>



The external parameters need to be created seperately. For this you can ask Daniel Klocke from DWD.

=== Prepare the boundary data ===
To prepare the boundary data, you need to install the DWD icontools. These are needed to remap the boundary data to the //boundary frame grid//. Check also the [[models:pot-pourri:how_to:icon_quick_start_guide|ICON quick start guide]] on this under "How to?".

<code bash>
git clone git@git.mpimet.mpg.de:dwd_icon_tools.git
</code>

There is a detailed documentation on the DWD ICON tools which you can find in your icontools directory:
<code bash>
dwd_icon_tools/doc/icontools doc.pdf.
</code>


The boundary data are necessary to run limited area simulations. In our case they are red 6-hourly from ICON.
We want to use **ECHAM data** to create the boundary files.

In the [[https://code.mpimet.mpg.de/projects/holocene/wiki/List_of_simulations|List of Simulations]] you find the paths to all conducted holocene simulations and their output. We use the **slo0021a** run. 
ICON needs the boundary data files for each time step seperately.
The following steps help you to prepare the boundary data:

1. __convert grib to nc and convert echam code numbers to variable names__
<code bash>
cdo -f nc -t echam copy infile file1
</code>

2. __convert spectral to grid-point grid: needed because ECHAM is a spectral grid model__
<code bash>
cdo sp2gp file1 file2
</code>

3. __use the afterburner to get the vertical velocity from u and v__
<code bash>
cdo after file2 file3 <<EON
CODE=131,132,135 TYPE=20
EON
</code>

4. __extract only the needed variables for the boundary data__
<code bash>
ncks -v hyai,hyam,hybi,hybm,lev,t,tsw,aps,q,xi,xl,tsurf,geosp,var131,var132,var135 file3 file4
</code>

5. __step in between (not necessary) and rename the variables: you can do this step directly in the script for mapping the boundary data to the boundary frame grid__
<code bash>
cdo chname,var131,U -chname,var132,V -chname,var135,OMEGA -chname,tsw,SST -chname,aps,PS -chname,q,QV -chname,t,T -chname,xi,QI -chname,xl,QC -chname,geosp,GEOSP file4 file5
</code>

6. __because there were problems with the time axis in ICON do:__
<code bash>
cdo -a setcalendar,standard file5 file6
ncatted -a units,time,m,c,day as %Y%m%d.%f file6 file7
</code>

With the following script you can split the yearly files into one-time step nc-files and seperate the needed variables:

<hidden boundary_data.py>
<code bash>
import numpy as np
import netCDF4 as nc4
from   netCDF4 import Dataset
import datetime
import os


MainFolder="/work/mh0731/m300674/ICON/boundary_data/original_data/only_one/"
for (path, dirs, files) in os.walk(MainFolder):
    for ncfile in files:
        if ncfile[-3:]=='.nc':
            ncfile=os.path.join(path,ncfile)
            ncfile=Dataset(ncfile, 'r', 'NETCDF4')
            lon_in = ncfile.variables['lon'][:]
            lat_in = ncfile.variables['lat'][:]
            time_in = ncfile.variables['time'][:]
            lev_in = ncfile.variables['lev'][:]
            T_in  = ncfile.variables['T'][:]     # Air Temperature
            QV_in  = ncfile.variables['QV'][:]     # specific humidity
            GEOSP_in = ncfile.variables['GEOSP'][:] # surface geopotential
            QI_in = ncfile.variables['QI'][:]   # cloud ice
            QC_in = ncfile.variables['QC'][:]   # cloud water
            U_in = ncfile.variables['U'][:]
            V_in = ncfile.variables['V'][:]
            OMEGA_in = ncfile.variables['OMEGA'][:]
            PS_in = ncfile.variables['PS'][:]
            SST_in = ncfile.variables['SST'][:]
            ncfile.close()

   # ----- CREATE 6-HOURLY FILES ---------------------------

            for t in range(0,len(time_in)):

               number_dec = time_in[t]-int(time_in[t])
               number_int = int(time_in[t])
               print(number_dec)

               if number_dec < 0.5:
                   res = str(number_int) + "_0"+ str(int(number_dec * 24))
               else:
                   res = str(number_int) + "_" + str(int(number_dec * 24))

               print('t',t,time_in[t])

               f_out = nc4.Dataset('/work/mh0731/m300674/ICON/boundary_data/splitted_boundary_data/boundary_conditions_%s.nc' % res,'w', format='NETCDF4_CLASSIC')
               print('boundary_conditions_%s.nc' % res)

               f_out.createDimension('lon', len(lon_in))
               f_out.createDimension('lat', len(lat_in))
               f_out.createDimension('time', None)
               f_out.createDimension('lev',len(lev_in))
               
                #------ from first file -----------------------------------
               lat = f_out.createVariable('lat', np.float,'lat')
               lat[:] = lat_in

               lon = f_out.createVariable('lon', np.float,'lon')
               lon[:] = lon_in

               time = f_out.createVariable('time',np.double,'time')
               time[0] = time_in[t]

               lev = f_out.createVariable('lev',np.double,'lev')
               lev[:] = lev_in      

               T = f_out.createVariable('T',np.float,('time','lev','lat','lon'))
               T[0,:,:,:] = T_in[t,:,:,:]

               PS = f_out.createVariable('PS',np.float,('time','lat','lon'))
               PS[0,:,:] = PS_in[t,:,:]

               QV = f_out.createVariable('QV',np.float,('time','lev','lat','lon'))
               QV[0,:,:,:] = QV_in[t,:,:,:]

               U = f_out.createVariable('U',np.float,('time','lev','lat','lon'))
               U[0,:,:,:] = U_in[t,:,:,:]

               V = f_out.createVariable('V',np.float,('time','lev','lat','lon'))
               V[0,:,:,:] = V_in[t,:,:,:]

               OMEGA = f_out.createVariable('OMEGA',np.float,('time','lev','lat','lon'))
               OMEGA[0,:,:,:] = OMEGA_in[t,:,:,:]

               QI = f_out.createVariable('QI',np.float,('time','lev','lat','lon'))
               QI[0,:,:,:] = QI_in[t,:,:,:]

               QC = f_out.createVariable('QC',np.float,('time','lev','lat','lon'))
               QC[0,:,:,:] = QC_in[t,:,:,:]
               
               GEOSP = f_out.createVariable('GEOSP',np.float,('time','lat','lon'))
               GEOSP[0,:,:] = GEOSP_in[t,:,:]

               SST = f_out.createVariable('SST',np.float,('time','lat','lon'))
               SST[0,:,:] = SST_in[t,:,:]

               #----local attributes ---------
               lon.units = 'degrees_east'
               lat.units = 'degrees_north'
               T.long_name = 'Air Temperature'
               T.units = 'K'
               PS.long_name = 'Surface pressure'
               PS.units = 'Pa'
               QV.long_name = 'Specific humidity'
               QV.units = 'kg/kg'
               U.long_name = 'horizontal wind component u'
               U.units = 'm/s'
               V.long_name = 'horizontal wind component v'
               V.units = 'm/s' 
               OMEGA.long_name = 'vertical wind velocity'
               OMEGA.units = 'Pa/s'
               QI.long_name = 'cloud ice'
               QI.units = 'kg/kg'
               QC.long_name = 'cloud water'
               QC.units = 'kg/kg'
               GEOSP.longname = 'Surface geopotential'
               GEOSP.units = 'm2 s-2'
               SST.longname = 'Sea Surface Temperature'
               SST.units = 'K'
 
               #--------------------------------------------------------------------

               f_out.close()
               
</code>
</hidden>


The variables which are needed are : **U, V, W, PS, T, SST, GEOSP, QC, QI, QV**
Because of the high time resolution the boundary data are needed, it is useful to minimize the file size. Therefore we map the boundary data on a boundary frame grid. To do so, we use the DWD ICON tools. Here is a setup of the script:
The script can be found under the following path: **''/work/mh0731/m300674/ICON/boundary_data/xce_remap_lbc_ECHAM''**

In all of the scripts you have to adjust the paths to your own directories and you can adjust variable names, attributes etc.


=== Prepare the initial file ===

The initial file for the Holocene simulation consists of data from the ECHAM holocene Data, IFS data and some variables which are set to zero.
You can use your first boundary file, **before** you remapped it onto the boundary frame. 
Because there are no variables for soil moisture in the ECHAM data (as they are needed in ICON) we use IFS data for the soil initialization. Variables you need to extract from an IFS file: **SMIL1, SMIL2, SMIL3, SMIL4, STL1, STL2, STL3, STL4, SKT, CI, LSM**. 
Variables you can set to zero are mainly ice and snow variables, which can be neglected in our domain: **W_SNOW, T_SNOW, W_I. RHO_SNOW**

<WRAP important>
Do not forget to adjust the time axis to be identical in all files
</WRAP>


==== Running the high resolution nesting simulation ====

You already prepared the grids and external parameters, as well as the boundary data. The latter are needed for the coarsest resolution only. The higher resolution domains are forced  by the coarser ones respectively.
First, you should check when the soil moisture equilibrates and chose a date to start the nesting simulation after that point. 
The first-guess file need to be remapped to all of your domains. 

The nested simulation will be initialized from first guess variables. You can write out these variables using **"group:dwd_fg_atm_vars","group:dwd_fg_sfc_vars"** in your output_nml. 
<WRAP important>
You have to write out only the one instant time step from which you want to start your nesting simulation
</WRAP>

Create the runscript for the nesting simulation mentioned above. With this you can start the nested Holocene simulations easily.