convert ICON SCM output to DEPHY format

convert ICON SCM output to DEPHY format#

Code to read ICON SCM output files and write to DEPHY format (NetCDF)

Contributed by Martin Köhler from DWD (based on Ann. Fridlind’s example)

Import libraries#

import xarray as xr
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import csv
import os
import netCDF4
import datetime as dt
from netCDF4 import Dataset

### Functions for optical depth calculation

def intialise_sizeDistri():
    global a_CD,b_CD,nuCD,muCD,a_RD,b_RD,nuRD,muRD,a_ID,b_ID,nuID,muID,a_SD,b_SD,nuSD,muSD
    global a_GD,b_GD,nuGD,muGD,xxCD,xxRD,xxID,xxSD,xxGD,ddCD,ddRD,ddID,ddSD,ddGD

    #collection of parameters of size distributions
    #cloud= cloud_nue1mue1:
    a_CD = .124 #m/kg^beta      #D=ax^b
    b_CD = 1./3.
    nuCD = 1.                    #nu of f(x) cloud
    muCD = 1.                    #mu  "  " cloud

    #rain= rainSBB:
    a_RD = 0.124
    b_RD = 1./3.
    nuRD = 1.0                 
    muRD = 1./3.                  

    #ice= ice_cosmo5:
    a_ID = 0.835
    b_ID = 0.39
    nuID = 0.0              
    muID = 1./3. 

    #snow= snowSBB:
    a_SD = 5.13
    b_SD = 0.5
    nuSD = 0.0               
    muSD = 0.5

    #graupel= graupelhail_cosmo5:
    a_GD = 0.142
    b_GD = 0.314
    nuGD = 1.0                 
    muGD = 1./3. 

    #size ranges [x_min,x_max]
    xxCD = [4.2e-15, 2.6e-10] # kg
    xxRD = [2.6e-10, 6.5e-5] # kg
    xxID = [1.0e-12, 1.0e-5] # kg
    xxSD = [1.0e-10, 2.0e-5] # kg
    xxGD = [4.19e-09, 5.3e-4] # kg

    ddCD = a_CD*np.array(xxCD)**b_CD
    ddRD = a_RD*np.array(xxRD)**b_RD
    ddID = a_ID*np.array(xxID)**b_ID
    ddSD = a_SD*np.array(xxSD)**b_SD
    ddGD = a_GD*np.array(xxGD)**b_GD

def gamma_dist(x, A, nu, ll, mu):
    return A[None,:,:]*(x**nu)*np.exp(-ll[None,:,:]*(x**mu))

def transform_x2r(a_CD,b_CD,nuCD,muCD,qc,nc):
    xmCD = np.float64(qc)/np.float64(nc)
    dmCD = a_CD*xmCD**b_CD
    #lambda and A of size distribution (Seifert & Beheng 2006)
    lCD = (xmCD*gamma((nuCD+1)/muCD)/gamma((nuCD+2)/muCD))**(-muCD)
    Acd = muCD*nc*(lCD**((nuCD+1)/muCD))/gamma((nuCD+1)/muCD)
    #conversion constants to get to f(r)=K*r**psi*exp(-delta*r**xi)
    Kcd = (Acd/b_CD)*(2/a_CD)**((nuCD+1)/b_CD)
    dltCD = lCD*(2/a_CD)**(muCD/b_CD)
    psiCD = (nuCD+1)/b_CD-1
    xiCD = muCD/b_CD
    return Kcd,dltCD,psiCD,xiCD

def Ap_r(Dlim, K, psi, dlt, xi, n=1000):
    r_bounds = np.linspace(*Dlim/2,n)
    r_diff = np.diff(r_bounds)
    r_bin = (r_bounds[:-1] + r_bounds[1:]) / 2
    p = gamma_dist(r_bin[:, None, None], K, psi, dlt, xi)
    
    return np.sum(math.pi*r_bin[:, None, None]**2*p*r_diff[:, None, None], axis=0)

### Other functions for variable computations

#calculate vapour pressure
def calc_e(p,q):
    return q*p/(0.622 +(1-0.622)*q)

def calc_e_is(T):
     return np.exp(9.550426 - (5723.265 / T) + 3.53068 * np.log(T) - 0.00728332 * T)

def potential_temperature(temperature, pressure, reference_pressure=100000):
    """
    Calculate potential temperature from temperature and pressure.
    
    Args:
    - temperature: Temperature in Kelvin
    - pressure: Pressure in Pascals
    - reference_pressure: Reference pressure in Pascals (default: 100000 Pa)
    
    Returns:
    - Potential temperature in Kelvin
    """
    R = 287.05  # Specific gas constant for dry air (J/(kg·K))
    Cp = 1005  # Specific heat capacity at constant pressure for dry air (J/(kg·K))
    
    potential_temp = temperature * (reference_pressure / pressure) ** (R / Cp)
    
    return potential_temp

Specify directory locations#

If on the ARM JupyterHub, it is recommended to create and specify a local directory that is outside of the COMBLE-MIP repository to upload raw model output files in your model’s format.

Processed domain-mean outputs are invited for commit to the GitHub repository on a user-specified branch under /comble-mip/output_scm/YOUR_MODEL_NAME/sandbox/YOUR_OUTPUT_FILE_NAME. These can be committed and removed at any time.

If you are able to make a run without ice, it is requested to append ‘noice’ to YOUR_OUTPUT_FILE_NAME, so that it can readily be automatically compared with the baseline and other liquid-only runs.

# specify input and output file names: versions of ICON-SCM with and without ice

# Phys: fixed CCN, default z0
#my_input_suffix  = 'FixN_def_z0.nc'
#my_output_suffix = 'FixN_def_z0.nc'

# Phys: fixed CCN, COMBLE z0
#my_input_suffix  = 'FixN.nc'
#my_output_suffix = 'FixN.nc'

# Phys: fixed CCN, COMBLE z0, no ice
my_input_suffix  = 'FixN_noice.nc'
my_output_suffix = 'FixN_noice.nc'

# Phys without ice
# my_input_suffix = 'entppe_noice.nc'
# my_output_suffix = 'Phys_FixN_noice.nc'

# Tun1
# my_input_suffix = 't698ml.nc'
# my_output_suffix = 'Tun1_FixN.nc'

# Tun1 without ice
# my_input_suffix = 't698ml_noice.nc'
# my_output_suffix = 'Tun1_FixN_noice.nc'

# Tun2
# my_input_suffix = 't705ml.nc'
# my_output_suffix = 'Tun2_FixN.nc'

# Tun2 without ice
# my_input_suffix = 't705ml_noice.nc'
# my_output_suffix = 'Tun2_FixN_noice.nc'

# specify local source directories (with subdirectories for spin up over ice and restart over water)
my_input_dir = '/user-data-home/comble-mip/output_scm/icon-scm/input/'

# specify Github scratch directory where processed model output will be committed (automate later)
my_output_filename = 'ICON-SCM_' + my_output_suffix
my_gitdir = '/user-data-home/comble-mip/output_scm/icon-scm/sandbox/'

Read ICON-SCM data#

Read single file containing all output data#

input_filename = my_input_dir + 'ICON-SCM_' + my_input_suffix
model_data = xr.open_dataset(input_filename)

# check if the run contains ice variables
do_ice = bool(max(model_data['tqi'].values)>0.)
print('do_ice = ',do_ice)

# data variable equality
icon_snds = model_data

# full parameter list
model_data

do_ice =  True

<xarray.Dataset> Size: 1MB
Dimensions:           (time: 121, height: 90, bnds: 2, height_2: 91,
                       height_3: 91)
Coordinates:
  * time              (time) datetime64[ns] 968B 2020-03-12T22:00:00 ... 2020...
    clon              float64 8B ...
    clat              float64 8B ...
  * height            (height) float64 720B 1.0 2.0 3.0 4.0 ... 88.0 89.0 90.0
  * height_2          (height_2) float64 728B 1.0 2.0 3.0 4.0 ... 89.0 90.0 91.0
  * height_3          (height_3) float64 728B 1.0 2.0 3.0 4.0 ... 89.0 90.0 91.0
Dimensions without coordinates: bnds
Data variables: (12/58)
    height_bnds       (height, bnds) float64 1kB ...
    height_2_bnds     (height_2, bnds) float64 1kB ...
    gz0               (time) float32 484B ...
    pres_sfc          (time) float32 484B ...
    t_g               (time) float32 484B ...
    umfl_s            (time) float32 484B ...
    ...                ...
    ddt_qv_gscp       (time, height) float32 44kB ...
    ddt_qv_turb       (time, height) float32 44kB ...
    ddt_qv_conv       (time, height) float32 44kB ...
    lwflx_up          (time, height_3) float32 44kB ...
    lwflx_dn          (time, height_3) float32 44kB ...
    acdnc             (height) float32 360B ...
Attributes:
    CDI:          Climate Data Interface version 2.3.0 (https://mpimet.mpg.de...
    Conventions:  CF-1.6
    uuidOfVGrid:  e0670a27-30fa-0503-fcfe-c68dbca12220
    source:       git@gitlab.dkrz.de:icon/icon-nwp.git@d18361cc8f8c43adc5d328...
    institution:  Max Planck Institute for Meteorology/Deutscher Wetterdienst
    title:        ICON simulation
    history:      Fri Jun 28 15:09:05 2024: ncwa -a ncells out_SCM_COMBLE_Fix...
    references:   see MPIM/DWD publications
    comment:      unknown user name on vh1l105 (Linux 4.18.0-477.15.1.el8_8.x...
    CDO:          Climate Data Operators version 2.3.0 (https://mpimet.mpg.de...
    NCO:          netCDF Operators version 4.8.1 (Homepage = http://nco.sf.ne...
### Calculate and append additional variables

#get array dimension of [time,height] array
dims=icon_snds["rho"].data.shape
dims

#surface friction velocity
#calculate using rho near surface
ustar=np.sqrt(np.sqrt(icon_snds["umfl_s"].data**2 + icon_snds["vmfl_s"].data**2)/icon_snds["rho"].data[:,dims[1]-1])
ustar[0]=np.nan      # fix 0 uv fluxes at first time-stepfor division in obukhov length calculation
model_data=model_data.assign(variables={'ustar': (('time'),ustar)})

#calculate surface roughness length
z0 = 9.0e-4 + 0.0*ustar
if 'def_z0' in my_output_filename:
    z0h = np.copy(z0)
else:
    z0h = 5.5e-6+0.0*np.copy(z0)

model_data=model_data.assign(variables={'z0':  (('time'),z0)})
model_data=model_data.assign(variables={'z0h': (('time'),z0h)})
model_data=model_data.assign(variables={'z0q': (('time'),z0h)})

#obukhov length
k = 0.4 #van Karman constant
g = 9.80665  #gravitational acceleration
rd= 287.04   #specific gas constant of dry air
cpd = 1004.64 # heat capcity of dry air at const. presssure
p0ref = 100000.0 # atm. reference pressure [Pa]
exner = np.exp(rd/cpd*np.log(icon_snds["pres_sfc"].data/p0ref))
theta_sfc = icon_snds["t_g"].data/exner

obk_length = - theta_sfc/ustar**3/(k*g*icon_snds["shfl_s"].data)
model_data=model_data.assign(variables={'ol': (('time'),obk_length)})

#atmosphere mass content of rain water
dz = (icon_snds['z_ifc'][0:-1].data-icon_snds['z_ifc'][1:].data)
rwp=np.sum(icon_snds['rho']*icon_snds['qr']*dz,axis=1)
model_data=model_data.assign(variables={'lwpr': (('time'),rwp.data)})

#atmosphere mass content of total ice water
iwp_tot=np.sum(icon_snds['rho']*(icon_snds['tot_qi_dia']+icon_snds['qs']+icon_snds['qg'])*dz,axis=1)
model_data=model_data.assign(variables={'iwp': (('time'),iwp_tot.data)})

#atmosphere mass content of total ice water
qi_tot=icon_snds['tot_qi_dia'] + icon_snds['qs'] + icon_snds['qg']
model_data=model_data.assign(variables={'qi_tot': (('time','height'),qi_tot.data)})

#calculate ice-phase precipitation flux at surface
pri=icon_snds['ice_gsp_rate'].data + icon_snds['snow_gsp_rate'].data + icon_snds['graupel_gsp_rate'].data
model_data=model_data.assign(variables={'pri': (('time'),pri)})

#calculate relative humidity over ice
rhi=calc_e(icon_snds['pres'].data,icon_snds['qv'].data)/calc_e_is(icon_snds['temp'].data)
model_data=model_data.assign(variables={'rhi': (('time','height'),rhi)})

#calculate potential temperature
theta=potential_temperature(icon_snds['temp'].data, icon_snds['pres'].data)
model_data=model_data.assign(variables={'theta': (('time','height'),theta)})

#calculate resolved TKE
#uprime=icon_snds["u"].data - icon_snds["u"].data.mean(axis=2)[:,:,None]
#vprime=icon_snds["v"].data - icon_snds["v"].data.mean(axis=2)[:,:,None]
#wprime=icon_snds["w"].data - icon_snds["w"].data.mean(axis=2)[:,:,None]
#tke_res=(0.5*np.sqrt(uprime**2 + vprime**2 + wprime**2)).mean(axis=2)
#wind_snds=wind_snds.assign(variables={'tke_res': (('time','height'),tke_res)})

#incloud droplet number concentration
#calculate nd fldmean only where qc>0.01g/kg
#nl_incloud=np.nanmean(np.where(icon_3D["qc"].data>1.0e-05,icon_3D["qnc"].data,np.nan),axis=2)
#set np.nan to "0"
#nl_incloud=np.where(np.isnan(nl_incloud),0,nl_incloud)
#icon_snds=icon_snds.assign(variables={'nlcic': (('time','height'),nl_incloud)})

#incloud total ice crystal number concentration
#computed only where qc>0.01g/kg & qi>1.0e-10
#ni_incloud=np.nanmean(np.where(np.logical_and(icon_3D["qc"].data>1.0e-05, 
#                               icon_3D["qi"].data+icon_3D["qs"].data+icon_3D["qg"].data > 1.0e-10),
#                               icon_3D["qni"].data+icon_3D["qns"].data+icon_3D["qng"].data,np.nan),axis=2)
#icon_snds=icon_snds.assign(variables={'niic': (('time','height'),ni_incloud)})

#calculate upward longwave flux
model_data=model_data.assign(variables={'thd_s': (('time'), icon_snds["thu_s"].data + icon_snds["thb_s"].data)})

model_data

<xarray.Dataset> Size: 1MB
Dimensions:           (time: 121, height: 90, bnds: 2, height_2: 91,
                       height_3: 91)
Coordinates:
  * time              (time) datetime64[ns] 968B 2020-03-12T22:00:00 ... 2020...
    clon              float64 8B ...
    clat              float64 8B ...
  * height            (height) float64 720B 1.0 2.0 3.0 4.0 ... 88.0 89.0 90.0
  * height_2          (height_2) float64 728B 1.0 2.0 3.0 4.0 ... 89.0 90.0 91.0
  * height_3          (height_3) float64 728B 1.0 2.0 3.0 4.0 ... 89.0 90.0 91.0
Dimensions without coordinates: bnds
Data variables: (12/70)
    height_bnds       (height, bnds) float64 1kB ...
    height_2_bnds     (height_2, bnds) float64 1kB ...
    gz0               (time) float32 484B ...
    pres_sfc          (time) float32 484B ...
    t_g               (time) float32 484B ...
    umfl_s            (time) float32 484B 0.0 0.1609 0.1061 ... 0.02222 0.022
    ...                ...
    iwp               (time) float32 484B 0.0 0.0 0.0 ... 0.03866 0.04319
    qi_tot            (time, height) float32 44kB 0.0 0.0 0.0 ... 0.0 0.0 0.0
    pri               (time) float32 484B 0.0 0.0 0.0 ... 1.022e-05 1.028e-05
    rhi               (time, height) float32 44kB 1.768e-06 2.588e-06 ... 0.5328
    theta             (time, height) float64 87kB 2.15e+03 2.015e+03 ... 266.0
    thd_s             (time) float32 484B 131.1 131.1 131.1 ... 283.0 282.8
Attributes:
    CDI:          Climate Data Interface version 2.3.0 (https://mpimet.mpg.de...
    Conventions:  CF-1.6
    uuidOfVGrid:  e0670a27-30fa-0503-fcfe-c68dbca12220
    source:       git@gitlab.dkrz.de:icon/icon-nwp.git@d18361cc8f8c43adc5d328...
    institution:  Max Planck Institute for Meteorology/Deutscher Wetterdienst
    title:        ICON simulation
    history:      Fri Jun 28 15:09:05 2024: ncwa -a ncells out_SCM_COMBLE_Fix...
    references:   see MPIM/DWD publications
    comment:      unknown user name on vh1l105 (Linux 4.18.0-477.15.1.el8_8.x...
    CDO:          Climate Data Operators version 2.3.0 (https://mpimet.mpg.de...
    NCO:          netCDF Operators version 4.8.1 (Homepage = http://nco.sf.ne...
#### Prepare output file in DEPHY format

Read requested variables list#

Variable description, naming, units, and dimensions.

# read list of requested variables
vars_mean_list = pd.read_excel('https://docs.google.com/spreadsheets/d/1Vl8jYGviet7EtXZuQiitrx4NSkV1x27aJAhxxjBb9zI/export?gid=1026157027&format=xlsx',
                              sheet_name='SCM')
pd.set_option('display.max_rows', None)
vars_mean_list
standard_name variable_id units dimensions comment (reported at end of each model physics time step, green=minimum, red=granularity enabling EMC2)
0 time time s – dimension, seconds since 2020-03-12 18:00:00
1 pressure_layer layer 1 – dimension, pressure layer number from 1 at sur...
2 air_pressure pa Pa time, layer pressure at mid-level points (native model lev...
3 layer_top_pressure pe Pa time, layer dimension, pressure at layer top points (used ...
4 surface_pressure ps Pa time –
5 surface_temperature ts K time –
6 surface_friction_velocity ustar m s-1 time –
7 surface_roughness_length_for_momentum_in_air z0 m time –
8 surface_roughness_length_for_heat_in_air z0h m time –
9 surface_roughness_length_for_humidity_in_air z0q m time –
10 surface_upward_sensible_heat_flux hfss W m-2 time –
11 surface_upward_latent_heat_flux hfls W m-2 time –
12 obukhov_length ol m time –
13 pbl_height pblh m time PBL scheme layer thickness (if available)
14 inversion_height zi m time sharpest vertical gradient in air_potential_te...
15 atmosphere_mass_content_of_liquid_cloud_water lwpc kg m-2 time scene (all sky); cloud water path in all class...
16 atmosphere_mass_content_of_rain_water lwpr kg m-2 time scene (all sky); rain water path in all classe...
17 atmosphere_mass_content_of_ice_water iwp kg m-2 time scene (all sky); all ice-phase hydrometeors in...
18 area_fraction_cover_of_hydrometeors cf 1 time all hydrometeors and cloud types (e.g., all ph...
19 area_fraction_cover_of_liquid_cloud cflc 1 time liquid cloud cover without precipitation, incl...
20 area_fraction_cover_of_convective_hydrometeors cfc 1 time all hydrometeors, default breakdown into conve...
21 optical_depth od 1 time scene (all sky); mid-visible, all hydrometeors...
22 optical_depth_of_liquid_cloud odlc 1 time scene (all sky); mid-visible, cloud liquid onl...
23 precipitation_flux_at_surface pr kg m-2 s-1 time scene (all sky); all hydrometeors
24 precipitation_flux_of_ice_at_surface pri kg m-2 s-1 time scene (all sky); all ice phase hydrometeors
25 toa_outgoing_longwave_flux rlut W m-2 time –
26 surface_downwelling_longwave_flux rlds W m-2 time –
27 surface_upwelling_longwave_flux rlus W m-2 time –
28 surface_sea_spray_number_flux ssaf m-2 s-1 time when using prognostic aerosol; emission only (...
29 height zf m time, layer altitude of layer mid-level points above sea s...
30 eastward_wind ua m s-1 time, layer –
31 northward_wind va m s-1 time, layer –
32 air_dry_density rhoa kg m-3 time, layer per kg dry air
33 air_temperature ta K time, layer –
34 water_vapor_mixing_ratio qv kg kg-1 time, layer –
35 relative_humidity hur 1 time, layer relative to liquid
36 relative_humidity_over_ice huri 1 time, layer relative to ice
37 air_potential_temperature theta K time, layer –
38 mass_mixing_ratio_of_cloud_liquid_water_in_air qlc kg kg-1 time, layer scene (all sky) per kg dry air; cloud water pa...
39 mass_mixing_ratio_of_rain_water_in_air qlr kg kg-1 time, layer rain water path only in all classes (e.g., con...
40 mass_mixing_ratio_of_ice_water_in_air qi kg kg-1 time, layer all ice water path in all classes (e.g., conve...
41 area_fraction_of_hydrometeors fh 1 time, layer all hydrometeors and cloud types (e.g., all ph...
42 area_fraction_of_liquid_cloud flc 1 time, layer liquid cloud cover without precipitation, incl...
43 area_fraction_of_convective_hydrometeors fc 1 time, layer all hydrometeors, default breakdown into conve...
44 precipitation_flux_in_air prf kg m-2 s-1 time, layer scene (all sky); all hydrometeors
45 precipitation_flux_in_air_in_ice_phase prfi kg m-2 s-1 time, layer scene (all sky); all ice phase hydrometeors
46 specific_turbulent_kinetic_energy tke m2 s-2 time, layer –
47 dissipation_rate_of_turbulent_kinetic_energy eps m2 s-3 time, layer report as negative
48 zonal_momentum_flux uw m2 s-2 time, layer parameterized turbulent flux
49 meridional_momentum_flux vw m2 s-2 time, layer parameterized turbulent flux
50 variance_of_upward_air_velocity w2 m2 s-2 time, layer parameterized turbulent flux
51 vertical_flux_potential_temperature wth K m s-1 time, layer parameterized turbulent flux
52 vertical_flux_liquid_ice_water_potential_tempe... vf_thli K m s-1 time, layer parameterized turbulent flux; include sediment...
53 vertical_flux_water_vapor wqv kg kg-1 m s-1 time, layer parameterized turbulent flux
54 vertical_flux_total_water vf_qt kg kg-1 m s-1 time, layer parameterized turbulent flux; vapor+all liquid...
55 convection_updraft_mass_flux cmfu kg m-2 s-1 time, layer –
56 convection_downdraft_mass_flux cmfd kg m-2 s-1 time, layer –
57 downwelling_longwave_flux_in_air rld W m-2 time, layer –
58 upwelling_longwave_flux_in_air rlu W m-2 time, layer –
59 tendency_of_air_potential_temperature_due_to_r... dth_rad K s-1 time, layer scene (all sky)
60 tendency_of_air_potential_temperature_due_to_m... dth_micro K s-1 time, layer including net condensation and precipitation i...
61 tendency_of_air_potential_temperature_due_to_m... dth_turb K s-1 time, layer including surface fluxes
62 tendency_of_water_vapor_mixing_ratio_due_to_mi... dq_micro s-1 time, layer including net condensation and precipitation i...
63 tendency_of_water_vapor_mixing_ratio_due_to_mi... dq_turb s-1 time, layer including surface fluxes
64 number_of_total_aerosol_mode1 na1 kg-1 time, layer when using prognostic aerosol; scene (all sky)...
65 number_of_total_aerosol_mode2 na2 kg-1 time, layer accumulation mode
66 number_of_total_aerosol_mode3 na3 kg-1 time, layer sea spray mode
67 tendency_of_aerosol_number_due_to_warm_microph... dna_micro_warm kg-1 s-1 time, layer activated and unactivated aerosol (sum over al...
68 tendency_of_aerosol_number_due_to_cold_microph... dna_micro_cold kg-1 s-1 time, layer activated and unactivated aerosol (sum over al...
69 tendency_of_aerosol_number_due_to_mixing dna_turb kg-1 s-1 time, layer activated and unactivated aerosol (sum over al...
70 tendency_of_ice_number_due_to_heterogeneous_fr... dni_het kg-1 s-1 time, layer sum of primary ice nucleation due to activatio...
71 tendency_of_ice_number_due_to_secondary_ice_pr... dni_sip kg-1 s-1 time, layer sum of secondary ice formation processes (e.g....
72 tendency_of_ice_number_due_to_homogeneous_free... dni_hom kg-1 s-1 time, layer ice nucleation source due to homogoeneous free...
73 mass_mixing_ratio_of_liquid_cloud_water_in_air... qlcs kg kg-1 time, layer scene (all sky) per kg dry air; default breakd...
74 mass_mixing_ratio_of_rain_water_in_air_stratiform qlrs kg kg-1 time, layer –
75 mass_mixing_ratio_of_ice_cloud_in_air_stratiform qics kg kg-1 time, layer default breakdown as for liquid; if other ice-...
76 mass_mixing_ratio_of_ice_precipitation_in_air_... qips kg kg-1 time, layer –
77 mass_mixing_ratio_of_liquid_cloud_water_in_air... qlcc kg kg-1 time, layer –
78 mass_mixing_ratio_of_rain_water_in_air_convective qlrc kg kg-1 time, layer –
79 mass_mixing_ratio_of_ice_cloud_in_air_convective qicc kg kg-1 time, layer –
80 mass_mixing_ratio_of_ice_precipitation_in_air_... qipc kg kg-1 time, layer –
81 number_of_liquid_cloud_droplets_in_air_stratiform nlcs kg-1 time, layer scene (all sky) per kg dry air; if other categ...
82 number_of_rain_drops_in_air_stratiform nlrs kg-1 time, layer –
83 number_of_ice_cloud_crystals_in_air_stratiform nics kg-1 time, layer if other ice-phase categories are used, provid...
84 number_of_ice_precipitation_crystals_in_air_st... nips kg-1 time, layer –
85 effective_radius_of_liquid_cloud_droplets_conv... relcc m time, layer EMC2 uses effective radius for any hydrometeor...
86 effective_radius_of_rain_convective relrc m time, layer –
87 effective_radius_of_ice_cloud_convective reicc m time, layer –
88 effective_radius_of_ice_precipitation_convective reipc m time, layer –
89 area_fraction_of_liquid_cloud_stratiform flcs 1 time, layer EMC2 uses area fraction profiles for all hydro...
90 area_fraction_of_rain_stratiform flrs 1 time, layer –
91 area_fraction_of_ice_cloud_stratiform fics 1 time, layer if other ice categories are used, provide addi...
92 area_fraction_of_ice_precipitation_stratiform fips 1 time, layer –
93 area_fraction_of_liquid_cloud_convective flcc 1 time, layer –
94 area_fraction_of_rain_convective flrc 1 time, layer –
95 area_fraction_of_ice_cloud_convective ficc 1 time, layer –
96 area_fraction_of_ice_precipitation_convective fipc 1 time, layer –
97 mass_weighted_fall_speed_of_liquid_cloud_water... vmlcs m s-1 time, layer EMC2 uses mass-weighted fall-speed profiles fo...
98 mass_weighted_fall_speed_of_rain_stratiform vmlrs m s-1 time, layer –
99 mass_weighted_fall_speed_of_ice_cloud_stratiform vmics m s-1 time, layer if other ice-phase categories are used, provid...
100 mass_weighted_fall_speed_of_ice_precipitation_... vmips m s-1 time, layer –
101 mass_weighted_fall_speed_of_liquid_cloud_water... vmlcc m s-1 time, layer –
102 mass_weighted_fall_speed_of_rain_convective vmlrc m s-1 time, layer –
103 mass_weighted_fall_speed_of_cloud_ice_crystals... vmicc m s-1 time, layer –
104 mass_weighted_fall_speed_of_ice_precipitation_... vmipc m s-1 time, layer –

Match ICON variables to requested outputs#

Expand the table to include columns that indicate ICON model variable names and any conversion factor.

# drop comments
vars_mean_list = vars_mean_list.drop(columns='comment (reported at end of each model physics time step, green=minimum, red=granularity enabling EMC2)')

# add columns to contain model output name and units conversion factors
vars_mean_list = vars_mean_list.assign(model_name='missing data',conv_factor=1.0)

# match to ICON variable names and specify conversion factors
for index in vars_mean_list.index:
    standard_name = vars_mean_list.standard_name.iat[index]
    if standard_name=='air_pressure': 
        vars_mean_list.model_name.iat[index] = 'pres'
    if standard_name=='layer_top_pressure': 
        vars_mean_list.model_name.iat[index] = 'pres_ifc'
    if standard_name=='surface_pressure': 
        vars_mean_list.model_name.iat[index] = 'pres_sfc'
#        vars_mean_list.conv_factor.iat[index] = 100.
    if standard_name=='surface_temperature': 
        vars_mean_list.model_name.iat[index] = 't_g'
    if standard_name=='surface_friction_velocity': 
        vars_mean_list.model_name.iat[index] = 'ustar'
    if standard_name=='surface_roughness_length_for_momentum_in_air': 
        vars_mean_list.model_name.iat[index] = 'z0'
    if standard_name=='surface_roughness_length_for_heat_in_air': 
        vars_mean_list.model_name.iat[index] = 'z0h'
    if standard_name=='surface_roughness_length_for_humidity_in_air': 
        vars_mean_list.model_name.iat[index] = 'z0q'
    if standard_name=='surface_upward_sensible_heat_flux': 
        vars_mean_list.model_name.iat[index] = 'shfl_s'
        vars_mean_list.conv_factor.iat[index] = -1.
    if standard_name=='surface_upward_latent_heat_flux': 
        vars_mean_list.model_name.iat[index] = 'lhfl_s'
        vars_mean_list.conv_factor.iat[index] = -1.
    if standard_name=='obukhov_length': 
        vars_mean_list.model_name.iat[index] = 'ol'
 #   if standard_name=='pbl_height': 
 #       vars_mean_list.model_name.iat[index] = 'pblht_bp'
 #   if standard_name=='inversion_height': 
 #       vars_mean_list.model_name.iat[index] = 'pblht_th'
    if standard_name=='atmosphere_mass_content_of_liquid_cloud_water': 
        vars_mean_list.model_name.iat[index] = 'tqc_dia'
#        vars_mean_list.conv_factor.iat[index] = 0.001
    if standard_name=='atmosphere_mass_content_of_rain_water': 
        vars_mean_list.model_name.iat[index] = 'tqr'
#        vars_mean_list.conv_factor.iat[index] = 0.001
    if do_ice:
        if standard_name=='atmosphere_mass_content_of_ice_water': 
            vars_mean_list.model_name.iat[index] = 'iwp'
#            vars_mean_list.conv_factor.iat[index] = 0.001
    if standard_name=='area_fraction_cover_of_hydrometeors': 
        vars_mean_list.model_name.iat[index] = 'clct'
#    if standard_name=='area_fraction_cover_of_liquid_cloud': 
#        vars_mean_list.model_name.iat[index] = ''
    if standard_name=='area_fraction_cover_of_convective_hydrometeors': 
        vars_mean_list.model_name.iat[index] = 'clct'
#   if standard_name=='optical_depth': 
#       vars_mean_list.model_name.iat[index] = 'tau'
#    if standard_name=='optical_depth_of_liquid_water': 
#        vars_mean_list.model_name.iat[index] = ''
    if standard_name=='precipitation_flux_at_surface': 
        vars_mean_list.model_name.iat[index] = 'tot_prec_rate'
#        vars_mean_list.conv_factor.iat[index] = 1./86400
    if standard_name=='precipitation_flux_of_ice_at_surface': 
        vars_mean_list.model_name.iat[index] = 'pri'
#        vars_mean_list.conv_factor.iat[index] = 1./86400
    if standard_name=='toa_outgoing_longwave_flux': 
        vars_mean_list.model_name.iat[index] = 'thb_t'
        vars_mean_list.conv_factor.iat[index] = -1.0
    if standard_name=='surface_downwelling_longwave_flux': 
        vars_mean_list.model_name.iat[index] = 'thd_s'  
    if standard_name=='surface_upwelling_longwave_flux': 
        vars_mean_list.model_name.iat[index] = 'thu_s'  
    if standard_name=='height': 
        vars_mean_list.model_name.iat[index] = 'z_mc'
    if standard_name=='eastward_wind': 
        vars_mean_list.model_name.iat[index] = 'u'
    if standard_name=='northward_wind': 
        vars_mean_list.model_name.iat[index] = 'v'
    if standard_name=='air_dry_density': 
        vars_mean_list.model_name.iat[index] = 'rho'
    if standard_name=='air_temperature': 
        vars_mean_list.model_name.iat[index] = 'temp'
    if standard_name=='water_vapor_mixing_ratio': 
        vars_mean_list.model_name.iat[index] = 'qv'
    if standard_name=='relative_humidity': 
        vars_mean_list.model_name.iat[index] = 'rh'
        vars_mean_list.conv_factor.iat[index] = 0.01
    if standard_name=='relative_humidity_over_ice': 
        vars_mean_list.model_name.iat[index] = 'rhi'
#        vars_mean_list.conv_factor.iat[index] = 0.01
    if standard_name=='air_potential_temperature': 
        vars_mean_list.model_name.iat[index] = 'theta'
    if standard_name=='mass_mixing_ratio_of_cloud_liquid_water_in_air': 
        vars_mean_list.model_name.iat[index] = 'tot_qc_dia'
    if standard_name=='mass_mixing_ratio_of_rain_water_in_air': 
        vars_mean_list.model_name.iat[index] = 'qr'
    if do_ice: 
        if standard_name=='mass_mixing_ratio_of_ice_water_in_air': 
            vars_mean_list.model_name.iat[index] = 'qi_tot'
    if standard_name=='area_fraction_of_hydrometeors': 
        vars_mean_list.model_name.iat[index] = 'clc'
#    if standard_name=='area_fraction_of_liquid_cloud': 
#        vars_mean_list.model_name.iat[index] = 'lcf'
#    if standard_name=='area_fraction_of_convective_hydrometeors': 
#        vars_mean_list.model_name.iat[index] = 'clc'
#    if standard_name=='precipitation_flux_in_air': 
#        vars_mean_list.model_name.iat[index] = 'tot_prec_rate'
#        vars_mean_list.conv_factor.iat[index] = 1./86400
#    if do_ice:
#        if standard_name=='precipitation_flux_in_air_in_ice_phase': 
#            vars_mean_list.model_name.iat[index] = 'pri'
#            vars_mean_list.conv_factor.iat[index] = 1./86400
    if standard_name=='specific_turbulent_kinetic_energy': 
        vars_mean_list.model_name.iat[index] = 'tke'
#    if standard_name=='disspation_rate_of_turbulent_kinetic_energy': 
#        vars_mean_list.model_name.iat[index] = 'dissip_tke_turb'
#        vars_mean_list.conv_factor.iat[index] = -1.
#    if standard_name=='zonal_momentum_flux': 
#        vars_mean_list.model_name.iat[index] = 'uw_turb'
#    if standard_name=='meridional_momentum_flux': 
#        vars_mean_list.model_name.iat[index] = 'vw_turb'
#    if standard_name=='variance_of_upward_air_velocity': 
#        vars_mean_list.model_name.iat[index] = 'w2_turb'
#    if standard_name=='vertical_flux_potential_temperature': 
#        vars_mean_list.model_name.iat[index] = 'wth_turb'
#    if standard_name=='vertical_flux_liquid_water_potential_temperature': 
#        vars_mean_list.model_name.iat[index] = ''
#    if standard_name=='vertical_flux_water_vapor': 
#        vars_mean_list.model_name.iat[index] = 'wq_turb'
#    if standard_name=='vertical_flux_total_water': 
#        vars_mean_list.model_name.iat[index] = 'wqt_turb'
#    if standard_name=='convection_updraft_mass_flux': 
#        vars_mean_list.model_name.iat[index] = 'lwdp'
#    if standard_name=='convection_downdraft_mass_flux': 
#        vars_mean_list.model_name.iat[index] = 'lwdp'
    if standard_name=='downwelling_longwave_flux_in_air': 
        vars_mean_list.model_name.iat[index] = 'lwflx_dn'
    if standard_name=='upwelling_longwave_flux_in_air': 
        vars_mean_list.model_name.iat[index] = 'lwflx_up'
    if standard_name=='tendency_of_air_potential_temperature_due_to_radiative_heating': 
        vars_mean_list.model_name.iat[index] = 'ddt_temp_radlw'
#        vars_mean_list.conv_factor.iat[index] = 1./86400
    if standard_name=='tendency_of_air_potential_temperature_due_to_microphysics': 
        vars_mean_list.model_name.iat[index] = 'ddt_temp_gscp'
#        vars_mean_list.conv_factor.iat[index] = 1./86400
    if standard_name=='tendency_of_air_potential_temperature_due_to_mixing': 
        vars_mean_list.model_name.iat[index] = 'ddt_temp_turb'
#        vars_mean_list.conv_factor.iat[index] = 1./86400
    if standard_name=='tendency_of_water_vapor_mixing_ratio_due_to_microphysics': 
        vars_mean_list.model_name.iat[index] = 'ddt_qv_gscp'
#        vars_mean_list.conv_factor.iat[index] = 1./86400
    if standard_name=='tendency_of_water_vapor_mixing_ratio_due_to_mixing': 
        vars_mean_list.model_name.iat[index] = 'ddt_qv_turb'
#        vars_mean_list.conv_factor.iat[index] = 1./86400
#    if standard_name=='mass_mixing_ratio_of_liquid_cloud_water_in_air_stratiform': 
#        vars_mean_list.model_name.iat[index] = 'qcl'
#    if standard_name=='mass_mixing_ratio_of_rain_water_in_air_stratiform': 
#        vars_mean_list.model_name.iat[index] = 'qpl'
#    if do_ice:
#        if standard_name=='mass_mixing_ratio_of_ice_cloud_in_air_stratiform': 
#            vars_mean_list.model_name.iat[index] = 'qci'
#        if standard_name=='mass_mixing_ratio_of_ice_precipitation_in_air_stratiform': 
#            vars_mean_list.model_name.iat[index] = 'qpi'
#    if standard_name=='mass_mixing_ratio_of_liquid_cloud_water_in_air_convective': 
#        vars_mean_list.model_name.iat[index] = 'QCLmc'
#    if standard_name=='mass_mixing_ratio_of_rain_water_in_air_convective': 
#        vars_mean_list.model_name.iat[index] = 'QPLmc'
#    if do_ice:
#        if standard_name=='mass_mixing_ratio_of_ice_cloud_in_air_convective': 
#            vars_mean_list.model_name.iat[index] = 'QCImc'
#        if standard_name=='mass_mixing_ratio_of_ice_precipitation_in_air_convective': 
#            vars_mean_list.model_name.iat[index] = 'QPImc'
    if standard_name=='number_of_liquid_cloud_droplets_in_air_stratiform': 
        vars_mean_list.model_name.iat[index] = 'acdnc'
#    if standard_name=='number_of_rain_drops_in_air_stratiform': 
#        vars_mean_list.model_name.iat[index] = 'npl'
#    if do_ice:
#        if standard_name=='number_of_ice_cloud_crystals_in_air_stratiform': 
#            vars_mean_list.model_name.iat[index] = 'nci'
#        if standard_name=='number_of_ice_precipitation_crystals_in_air_stratiform': 
#            vars_mean_list.model_name.iat[index] = 'npi'
#   if standard_name=='effective_radius_of_liquid_cloud_droplets_convective': 
#        vars_mean_list.model_name.iat[index] = 're_mccl'
#    if standard_name=='effective_radius_of_rain_convective': 
#        vars_mean_list.model_name.iat[index] = 're_mcpl'
#    if do_ice:
#        if standard_name=='effective_radius_of_ice_cloud_convective': 
#            vars_mean_list.model_name.iat[index] = 're_mcci'
#        if standard_name=='effective_radius_of_ice_precipitation_convective': 
#            vars_mean_list.model_name.iat[index] = 're_mcpi'
#    if standard_name=='area_fraction_of_liquid_cloud_stratiform': 
#        vars_mean_list.model_name.iat[index] = 'cldsscl'
#    if standard_name=='area_fraction_of_rain_stratiform': 
#        vars_mean_list.model_name.iat[index] = 'cldsspl'
#    if do_ice:
#        if standard_name=='area_fraction_of_ice_cloud_stratiform': 
#            vars_mean_list.model_name.iat[index] = 'cldssci'
#        if standard_name=='area_fraction_of_ice_precipitation_stratiform': 
#            vars_mean_list.model_name.iat[index] = 'cldsspi'
#    if standard_name=='area_fraction_of_liquid_cloud_convective': 
#        vars_mean_list.model_name.iat[index] = 'cldmccl'
#    if standard_name=='area_fraction_of_rain_convective': 
#        vars_mean_list.model_name.iat[index] = 'cldmcpl'
#    if do_ice:
#        if standard_name=='area_fraction_of_ice_cloud_convective': 
#            vars_mean_list.model_name.iat[index] = 'cldmcci'
#        if standard_name=='area_fraction_of_ice_precipitation_convective': 
#            vars_mean_list.model_name.iat[index] = 'cldmcpi'
#    if standard_name=='mass_weighted_fall_speed_of_liquid_cloud_water_stratiform': 
#        vars_mean_list.model_name.iat[index] = 'vm_sscl'
#    if standard_name=='mass_weighted_fall_speed_of_rain_stratiform': 
#        vars_mean_list.model_name.iat[index] = 'vm_sspl'
#    if do_ice:
#        if standard_name=='mass_weighted_fall_speed_of_ice_cloud_stratiform': 
#            vars_mean_list.model_name.iat[index] = 'vm_ssci'
#        if standard_name=='mass_weighted_fall_speed_of_ice_precipitation_stratiform': 
#            vars_mean_list.model_name.iat[index] = 'vm_sspi'
#    if standard_name=='mass_weighted_fall_speed_of_liquid_cloud_water_convective': 
#        vars_mean_list.model_name.iat[index] = 'vm_mccl'
#    if standard_name=='mass_weighted_fall_speed_of_rain_convective': 
#        vars_mean_list.model_name.iat[index] = 'vm_mcpl'
#    if do_ice:
#        if standard_name=='mass_weighted_fall_speed_of_cloud_ice_crystals_convective': 
#            vars_mean_list.model_name.iat[index] = 'vm_mcci'
#        if standard_name=='mass_weighted_fall_speed_of_ice_precipitation_convective': 
#            vars_mean_list.model_name.iat[index] = 'vm_mcpi'

vars_mean_list[2:] # echo variables (first two rows are dimensions)
standard_name variable_id units dimensions model_name conv_factor
2 air_pressure pa Pa time, layer pres 1.00
3 layer_top_pressure pe Pa time, layer pres_ifc 1.00
4 surface_pressure ps Pa time pres_sfc 1.00
5 surface_temperature ts K time t_g 1.00
6 surface_friction_velocity ustar m s-1 time ustar 1.00
7 surface_roughness_length_for_momentum_in_air z0 m time z0 1.00
8 surface_roughness_length_for_heat_in_air z0h m time z0h 1.00
9 surface_roughness_length_for_humidity_in_air z0q m time z0q 1.00
10 surface_upward_sensible_heat_flux hfss W m-2 time shfl_s -1.00
11 surface_upward_latent_heat_flux hfls W m-2 time lhfl_s -1.00
12 obukhov_length ol m time ol 1.00
13 pbl_height pblh m time missing data 1.00
14 inversion_height zi m time missing data 1.00
15 atmosphere_mass_content_of_liquid_cloud_water lwpc kg m-2 time tqc_dia 1.00
16 atmosphere_mass_content_of_rain_water lwpr kg m-2 time tqr 1.00
17 atmosphere_mass_content_of_ice_water iwp kg m-2 time iwp 1.00
18 area_fraction_cover_of_hydrometeors cf 1 time clct 1.00
19 area_fraction_cover_of_liquid_cloud cflc 1 time missing data 1.00
20 area_fraction_cover_of_convective_hydrometeors cfc 1 time clct 1.00
21 optical_depth od 1 time missing data 1.00
22 optical_depth_of_liquid_cloud odlc 1 time missing data 1.00
23 precipitation_flux_at_surface pr kg m-2 s-1 time tot_prec_rate 1.00
24 precipitation_flux_of_ice_at_surface pri kg m-2 s-1 time pri 1.00
25 toa_outgoing_longwave_flux rlut W m-2 time thb_t -1.00
26 surface_downwelling_longwave_flux rlds W m-2 time thd_s 1.00
27 surface_upwelling_longwave_flux rlus W m-2 time thu_s 1.00
28 surface_sea_spray_number_flux ssaf m-2 s-1 time missing data 1.00
29 height zf m time, layer z_mc 1.00
30 eastward_wind ua m s-1 time, layer u 1.00
31 northward_wind va m s-1 time, layer v 1.00
32 air_dry_density rhoa kg m-3 time, layer rho 1.00
33 air_temperature ta K time, layer temp 1.00
34 water_vapor_mixing_ratio qv kg kg-1 time, layer qv 1.00
35 relative_humidity hur 1 time, layer rh 0.01
36 relative_humidity_over_ice huri 1 time, layer rhi 1.00
37 air_potential_temperature theta K time, layer theta 1.00
38 mass_mixing_ratio_of_cloud_liquid_water_in_air qlc kg kg-1 time, layer tot_qc_dia 1.00
39 mass_mixing_ratio_of_rain_water_in_air qlr kg kg-1 time, layer qr 1.00
40 mass_mixing_ratio_of_ice_water_in_air qi kg kg-1 time, layer qi_tot 1.00
41 area_fraction_of_hydrometeors fh 1 time, layer clc 1.00
42 area_fraction_of_liquid_cloud flc 1 time, layer missing data 1.00
43 area_fraction_of_convective_hydrometeors fc 1 time, layer missing data 1.00
44 precipitation_flux_in_air prf kg m-2 s-1 time, layer missing data 1.00
45 precipitation_flux_in_air_in_ice_phase prfi kg m-2 s-1 time, layer missing data 1.00
46 specific_turbulent_kinetic_energy tke m2 s-2 time, layer tke 1.00
47 dissipation_rate_of_turbulent_kinetic_energy eps m2 s-3 time, layer missing data 1.00
48 zonal_momentum_flux uw m2 s-2 time, layer missing data 1.00
49 meridional_momentum_flux vw m2 s-2 time, layer missing data 1.00
50 variance_of_upward_air_velocity w2 m2 s-2 time, layer missing data 1.00
51 vertical_flux_potential_temperature wth K m s-1 time, layer missing data 1.00
52 vertical_flux_liquid_ice_water_potential_tempe... vf_thli K m s-1 time, layer missing data 1.00
53 vertical_flux_water_vapor wqv kg kg-1 m s-1 time, layer missing data 1.00
54 vertical_flux_total_water vf_qt kg kg-1 m s-1 time, layer missing data 1.00
55 convection_updraft_mass_flux cmfu kg m-2 s-1 time, layer missing data 1.00
56 convection_downdraft_mass_flux cmfd kg m-2 s-1 time, layer missing data 1.00
57 downwelling_longwave_flux_in_air rld W m-2 time, layer lwflx_dn 1.00
58 upwelling_longwave_flux_in_air rlu W m-2 time, layer lwflx_up 1.00
59 tendency_of_air_potential_temperature_due_to_r... dth_rad K s-1 time, layer ddt_temp_radlw 1.00
60 tendency_of_air_potential_temperature_due_to_m... dth_micro K s-1 time, layer ddt_temp_gscp 1.00
61 tendency_of_air_potential_temperature_due_to_m... dth_turb K s-1 time, layer ddt_temp_turb 1.00
62 tendency_of_water_vapor_mixing_ratio_due_to_mi... dq_micro s-1 time, layer ddt_qv_gscp 1.00
63 tendency_of_water_vapor_mixing_ratio_due_to_mi... dq_turb s-1 time, layer ddt_qv_turb 1.00
64 number_of_total_aerosol_mode1 na1 kg-1 time, layer missing data 1.00
65 number_of_total_aerosol_mode2 na2 kg-1 time, layer missing data 1.00
66 number_of_total_aerosol_mode3 na3 kg-1 time, layer missing data 1.00
67 tendency_of_aerosol_number_due_to_warm_microph... dna_micro_warm kg-1 s-1 time, layer missing data 1.00
68 tendency_of_aerosol_number_due_to_cold_microph... dna_micro_cold kg-1 s-1 time, layer missing data 1.00
69 tendency_of_aerosol_number_due_to_mixing dna_turb kg-1 s-1 time, layer missing data 1.00
70 tendency_of_ice_number_due_to_heterogeneous_fr... dni_het kg-1 s-1 time, layer missing data 1.00
71 tendency_of_ice_number_due_to_secondary_ice_pr... dni_sip kg-1 s-1 time, layer missing data 1.00
72 tendency_of_ice_number_due_to_homogeneous_free... dni_hom kg-1 s-1 time, layer missing data 1.00
73 mass_mixing_ratio_of_liquid_cloud_water_in_air... qlcs kg kg-1 time, layer missing data 1.00
74 mass_mixing_ratio_of_rain_water_in_air_stratiform qlrs kg kg-1 time, layer missing data 1.00
75 mass_mixing_ratio_of_ice_cloud_in_air_stratiform qics kg kg-1 time, layer missing data 1.00
76 mass_mixing_ratio_of_ice_precipitation_in_air_... qips kg kg-1 time, layer missing data 1.00
77 mass_mixing_ratio_of_liquid_cloud_water_in_air... qlcc kg kg-1 time, layer missing data 1.00
78 mass_mixing_ratio_of_rain_water_in_air_convective qlrc kg kg-1 time, layer missing data 1.00
79 mass_mixing_ratio_of_ice_cloud_in_air_convective qicc kg kg-1 time, layer missing data 1.00
80 mass_mixing_ratio_of_ice_precipitation_in_air_... qipc kg kg-1 time, layer missing data 1.00
81 number_of_liquid_cloud_droplets_in_air_stratiform nlcs kg-1 time, layer acdnc 1.00
82 number_of_rain_drops_in_air_stratiform nlrs kg-1 time, layer missing data 1.00
83 number_of_ice_cloud_crystals_in_air_stratiform nics kg-1 time, layer missing data 1.00
84 number_of_ice_precipitation_crystals_in_air_st... nips kg-1 time, layer missing data 1.00
85 effective_radius_of_liquid_cloud_droplets_conv... relcc m time, layer missing data 1.00
86 effective_radius_of_rain_convective relrc m time, layer missing data 1.00
87 effective_radius_of_ice_cloud_convective reicc m time, layer missing data 1.00
88 effective_radius_of_ice_precipitation_convective reipc m time, layer missing data 1.00
89 area_fraction_of_liquid_cloud_stratiform flcs 1 time, layer missing data 1.00
90 area_fraction_of_rain_stratiform flrs 1 time, layer missing data 1.00
91 area_fraction_of_ice_cloud_stratiform fics 1 time, layer missing data 1.00
92 area_fraction_of_ice_precipitation_stratiform fips 1 time, layer missing data 1.00
93 area_fraction_of_liquid_cloud_convective flcc 1 time, layer missing data 1.00
94 area_fraction_of_rain_convective flrc 1 time, layer missing data 1.00
95 area_fraction_of_ice_cloud_convective ficc 1 time, layer missing data 1.00
96 area_fraction_of_ice_precipitation_convective fipc 1 time, layer missing data 1.00
97 mass_weighted_fall_speed_of_liquid_cloud_water... vmlcs m s-1 time, layer missing data 1.00
98 mass_weighted_fall_speed_of_rain_stratiform vmlrs m s-1 time, layer missing data 1.00
99 mass_weighted_fall_speed_of_ice_cloud_stratiform vmics m s-1 time, layer missing data 1.00
100 mass_weighted_fall_speed_of_ice_precipitation_... vmips m s-1 time, layer missing data 1.00
101 mass_weighted_fall_speed_of_liquid_cloud_water... vmlcc m s-1 time, layer missing data 1.00
102 mass_weighted_fall_speed_of_rain_convective vmlrc m s-1 time, layer missing data 1.00
103 mass_weighted_fall_speed_of_cloud_ice_crystals... vmicc m s-1 time, layer missing data 1.00
104 mass_weighted_fall_speed_of_ice_precipitation_... vmipc m s-1 time, layer missing data 1.00

Create DEPHY output file#

Write a single file to contain all domain-mean scalar and profile outputs. This code expects the write directory to be pre-existing (already created by the user). In the case that this output will be committed to the comble-mip GitHub repository, see above “Specify directory locations”.

#print(model_data)
variable_names = list(model_data.data_vars.keys())
print(variable_names)

['height_bnds', 'height_2_bnds', 'gz0', 'pres_sfc', 't_g', 'umfl_s', 'vmfl_s', 'shfl_s', 'lhfl_s', 'tcm', 'tch', 'tqv', 'tqc', 'tqi', 'tqr', 'tqs', 'tqg', 'tot_prec', 'tot_prec_rate', 'ice_gsp_rate', 'snow_gsp_rate', 'graupel_gsp_rate', 'thu_s', 'thb_s', 'thb_t', 'sob_s', 'tqv_dia', 'tqc_dia', 'tqi_dia', 'clct', 'z_mc', 'z_ifc', 'pres', 'pres_ifc', 'u', 'v', 'rho', 'temp', 'tke', 'exner', 'qv', 'tot_qc_dia', 'tot_qi_dia', 'qr', 'qs', 'qg', 'clc', 'rh', 'ddt_temp_radlw', 'ddt_temp_gscp', 'ddt_temp_turb', 'ddt_temp_pconv', 'ddt_qv_gscp', 'ddt_qv_turb', 'ddt_qv_conv', 'lwflx_up', 'lwflx_dn', 'acdnc', 'ustar', 'z0', 'z0h', 'z0q', 'ol', 'lwpr', 'iwp', 'qi_tot', 'pri', 'rhi', 'theta', 'thd_s']

model_data

<xarray.Dataset> Size: 1MB
Dimensions:           (time: 121, height: 90, bnds: 2, height_2: 91,
                       height_3: 91)
Coordinates:
  * time              (time) datetime64[ns] 968B 2020-03-12T22:00:00 ... 2020...
    clon              float64 8B ...
    clat              float64 8B ...
  * height            (height) float64 720B 1.0 2.0 3.0 4.0 ... 88.0 89.0 90.0
  * height_2          (height_2) float64 728B 1.0 2.0 3.0 4.0 ... 89.0 90.0 91.0
  * height_3          (height_3) float64 728B 1.0 2.0 3.0 4.0 ... 89.0 90.0 91.0
Dimensions without coordinates: bnds
Data variables: (12/70)
    height_bnds       (height, bnds) float64 1kB ...
    height_2_bnds     (height_2, bnds) float64 1kB ...
    gz0               (time) float32 484B ...
    pres_sfc          (time) float32 484B ...
    t_g               (time) float32 484B ...
    umfl_s            (time) float32 484B 0.0 0.1609 0.1061 ... 0.02222 0.022
    ...                ...
    iwp               (time) float32 484B 0.0 0.0 0.0 ... 0.03866 0.04319
    qi_tot            (time, height) float32 44kB 0.0 0.0 0.0 ... 0.0 0.0 0.0
    pri               (time) float32 484B 0.0 0.0 0.0 ... 1.022e-05 1.028e-05
    rhi               (time, height) float32 44kB 1.768e-06 2.588e-06 ... 0.5328
    theta             (time, height) float64 87kB 2.15e+03 2.015e+03 ... 266.0
    thd_s             (time) float32 484B 131.1 131.1 131.1 ... 283.0 282.8
Attributes:
    CDI:          Climate Data Interface version 2.3.0 (https://mpimet.mpg.de...
    Conventions:  CF-1.6
    uuidOfVGrid:  e0670a27-30fa-0503-fcfe-c68dbca12220
    source:       git@gitlab.dkrz.de:icon/icon-nwp.git@d18361cc8f8c43adc5d328...
    institution:  Max Planck Institute for Meteorology/Deutscher Wetterdienst
    title:        ICON simulation
    history:      Fri Jun 28 15:09:05 2024: ncwa -a ncells out_SCM_COMBLE_Fix...
    references:   see MPIM/DWD publications
    comment:      unknown user name on vh1l105 (Linux 4.18.0-477.15.1.el8_8.x...
    CDO:          Climate Data Operators version 2.3.0 (https://mpimet.mpg.de...
    NCO:          netCDF Operators version 4.8.1 (Homepage = http://nco.sf.ne...
# create DEPHY output file
dephy_filename = my_gitdir + my_output_filename
if os.path.exists(dephy_filename):
    os.remove(dephy_filename)
    print('The file ' + dephy_filename + ' has been deleted successfully')    
dephy_file = Dataset(dephy_filename,mode='w',format='NETCDF3_CLASSIC')
start_date = '2020-03-12T22:00:00Z'

# create global attributes
dephy_file.title='ICON SCM results for COMBLE-MIP case: fixed stratiform Nd and Ni'
dephy_file.reference='https://github.com/ARM-Development/comble-mip'
dephy_file.authors='Martin Köhler (martin.koehler@dwd.de), Anna Possner (apossner@iau.uni-frankfurt.de)'
dephy_file.source=input_filename
dephy_file.version=dt.datetime.now().strftime('%Y-%m-%d %H:%M:%S')
dephy_file.format_version='DEPHY SCM format version 1.6'
dephy_file.script='convert_ModelE3_SCM_output_to_dephy_format.ipynb'
dephy_file.startDate=start_date
dephy_file.force_geo=1
dephy_file.surfaceType='ocean'
dephy_file.surfaceForcing='ts'
dephy_file.lat='74.5 deg N'
dephy_file.dp='see pressure variable'
dims=icon_snds["pres"].data.shape
dephy_file.np=dims[1]

# create dimensions
nt = model_data.dims['time']
time = dephy_file.createDimension('time', nt)
time = dephy_file.createVariable('time', np.float64, ('time',))
time.units = 'seconds since ' + dephy_file.startDate
time.long_name = 'time'
# find time step and build time in seconds
time_str1 = str(model_data['time'].data[0]).split('.')[0]  # Removes fractional seconds
time_str2 = str(model_data['time'].data[1]).split('.')[0]  # Removes fractional seconds
time1 = dt.datetime.strptime(time_str1,'%Y-%m-%dT%H:%M:%S')
time2 = dt.datetime.strptime(time_str2,'%Y-%m-%dT%H:%M:%S')
delta_t = (time2-time1).total_seconds()
time[:] = (np.arange(nt)+1.)*delta_t

# vertical dimension: model layers called 'pressure_layer'
nl = dims[1]      # Number of levels
layer           = dephy_file.createDimension('layer', nl)
layer           = dephy_file.createVariable ('layer', np.int32, ('layer',))
layer.units     = 'unitless'
layer.long_name = 'model layer'
layer[:]        = np.arange(1, nl + 1)   # values from 1 to 90, 1 at the surface and 90 at the top


# Old fixed pressure layers
#nl = model_data.dims['pres']
#pa = dephy_file.createDimension('pa', nl)
#pa = dephy_file.createVariable('pa', np.float64, ('pa',))
#pa.units = 'Pa'
#pa.long_name = 'pressure'
#pa[:] = model_data['pres'].data*100.

# create and fill variables
for index in vars_mean_list.index[2:]:
    std_name = vars_mean_list.standard_name.iat[index]
#   print(std_name) # debug
    var_name = vars_mean_list.variable_id.iat[index]
    mod_name = vars_mean_list.model_name.iat[index]
    c_factor = vars_mean_list.conv_factor.iat[index]
    if vars_mean_list.dimensions.iat[index]=='time':
        new_sca           = dephy_file.createVariable(var_name, np.float64, ('time'))
        new_sca.units     = vars_mean_list.units.iat[index]
        new_sca.long_name = std_name
        #print('variable:', mod_name, model_data[mod_name].data)
        if vars_mean_list.model_name.iat[index]!='missing data':
            new_sca[:]    = model_data[mod_name].data*c_factor
    if vars_mean_list.dimensions.iat[index]=='time, layer':
        new_snd           = dephy_file.createVariable(var_name, np.float64, ('time','layer'))
        new_snd.units     = vars_mean_list.units.iat[index]
        new_snd.long_name = std_name
        if vars_mean_list.model_name.iat[index]!='missing data':
            shape = model_data[mod_name].data.shape
            if shape[-1] == nl+1:
              new_snd[:]    = model_data[mod_name].data[:,0:nl]*c_factor
            else:
              new_snd[:]    = model_data[mod_name].data*c_factor

print(dephy_file)
dephy_file.close()

/tmp/ipykernel_105/2056566154.py:27: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  nt = model_data.dims['time']

<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF3_CLASSIC data model, file format NETCDF3):
    title: ICON SCM results for COMBLE-MIP case: fixed stratiform Nd and Ni
    reference: https://github.com/ARM-Development/comble-mip
    authors: Martin Köhler (martin.koehler@dwd.de), Anna Possner (apossner@iau.uni-frankfurt.de)
    source: /user-data-home/comble-mip/output_scm/icon-scm/input/ICON-SCM_FixN_noice.nc
    version: 2024-06-28 22:37:39
    format_version: DEPHY SCM format version 1.6
    script: convert_ModelE3_SCM_output_to_dephy_format.ipynb
    startDate: 2020-03-12T22:00:00Z
    force_geo: 1
    surfaceType: ocean
    surfaceForcing: ts
    lat: 74.5 deg N
    dp: see pressure variable
    np: 90
    dimensions(sizes): time(121), layer(90)
    variables(dimensions): float64 time(time), int32 layer(layer), float64 pa(time, layer), float64 pe(time, layer), float64 ps(time), float64 ts(time), float64 ustar(time), float64 z0(time), float64 z0h(time), float64 z0q(time), float64 hfss(time), float64 hfls(time), float64 ol(time), float64 pblh(time), float64 zi(time), float64 lwpc(time), float64 lwpr(time), float64 iwp(time), float64 cf(time), float64 cflc(time), float64 cfc(time), float64 od(time), float64 odlc(time), float64 pr(time), float64 pri(time), float64 rlut(time), float64 rlds(time), float64 rlus(time), float64 ssaf(time), float64 zf(time, layer), float64 ua(time, layer), float64 va(time, layer), float64 rhoa(time, layer), float64 ta(time, layer), float64 qv(time, layer), float64 hur(time, layer), float64 huri(time, layer), float64 theta(time, layer), float64 qlc(time, layer), float64 qlr(time, layer), float64 qi(time, layer), float64 fh(time, layer), float64 flc(time, layer), float64 fc(time, layer), float64 prf(time, layer), float64 prfi(time, layer), float64 tke(time, layer), float64 eps(time, layer), float64 uw(time, layer), float64 vw(time, layer), float64 w2(time, layer), float64 wth(time, layer), float64 vf_thli(time, layer), float64 wqv(time, layer), float64 vf_qt(time, layer), float64 cmfu(time, layer), float64 cmfd(time, layer), float64 rld(time, layer), float64 rlu(time, layer), float64 dth_rad(time, layer), float64 dth_micro(time, layer), float64 dth_turb(time, layer), float64 dq_micro(time, layer), float64 dq_turb(time, layer), float64 na1(time, layer), float64 na2(time, layer), float64 na3(time, layer), float64 dna_micro_warm(time, layer), float64 dna_micro_cold(time, layer), float64 dna_turb(time, layer), float64 dni_het(time, layer), float64 dni_sip(time, layer), float64 dni_hom(time, layer), float64 qlcs(time, layer), float64 qlrs(time, layer), float64 qics(time, layer), float64 qips(time, layer), float64 qlcc(time, layer), float64 qlrc(time, layer), float64 qicc(time, layer), float64 qipc(time, layer), float64 nlcs(time, layer), float64 nlrs(time, layer), float64 nics(time, layer), float64 nips(time, layer), float64 relcc(time, layer), float64 relrc(time, layer), float64 reicc(time, layer), float64 reipc(time, layer), float64 flcs(time, layer), float64 flrs(time, layer), float64 fics(time, layer), float64 fips(time, layer), float64 flcc(time, layer), float64 flrc(time, layer), float64 ficc(time, layer), float64 fipc(time, layer), float64 vmlcs(time, layer), float64 vmlrs(time, layer), float64 vmics(time, layer), float64 vmips(time, layer), float64 vmlcc(time, layer), float64 vmlrc(time, layer), float64 vmicc(time, layer), float64 vmipc(time, layer)
    groups:

Check output file#

dephy_check = xr.open_dataset(dephy_filename)
dephy_check

<xarray.Dataset> Size: 7MB
Dimensions:         (time: 121, layer: 90)
Coordinates:
  * time            (time) datetime64[ns] 968B 2020-03-12T22:10:00 ... 2020-0...
  * layer           (layer) int32 360B 1 2 3 4 5 6 7 8 ... 84 85 86 87 88 89 90
Data variables: (12/103)
    pa              (time, layer) float64 87kB ...
    pe              (time, layer) float64 87kB ...
    ps              (time) float64 968B ...
    ts              (time) float64 968B ...
    ustar           (time) float64 968B ...
    z0              (time) float64 968B ...
    ...              ...
    vmics           (time, layer) float64 87kB ...
    vmips           (time, layer) float64 87kB ...
    vmlcc           (time, layer) float64 87kB ...
    vmlrc           (time, layer) float64 87kB ...
    vmicc           (time, layer) float64 87kB ...
    vmipc           (time, layer) float64 87kB ...
Attributes: (12/14)
    title:           ICON SCM results for COMBLE-MIP case: fixed stratiform N...
    reference:       https://github.com/ARM-Development/comble-mip
    authors:         Martin Köhler (martin.koehler@dwd.de), Anna Possner (apo...
    source:          /user-data-home/comble-mip/output_scm/icon-scm/input/ICO...
    version:         2024-06-28 22:37:39
    format_version:  DEPHY SCM format version 1.6
    ...              ...
    force_geo:       1
    surfaceType:     ocean
    surfaceForcing:  ts
    lat:             74.5 deg N
    dp:              see pressure variable
    np:              90