Example: convert AOSCM output to DEPHY format

Example: convert AOSCM output to DEPHY format#

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

Contributed by Ann Fridlind from NASA/GISS, Michail Karalis from Stockholm University

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

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 ModelE3 with and without ice

# Phys
my_input_suffix = 'COMBLE_scm_in.nc'
my_prog_suffix = 'progvar.nc'
my_diag_suffix = 'diagvar.nc'
my_diag2_suffix = 'diagvar2.nc'
my_merged_suffix = 'aoscm.nc'

# Phys without ice
# my_input_suffix = 'entppe_noice.nc'
exper = 'FixN_def_z0'
my_output_suffix = 'AOSCM_'+exper+'_alt.nc'



# specify local source directories (with subdirectories for spin up over ice and restart over water)
# my_input_dir = '/data/home/fridlind/modelE3/'
my_input_dir = '~/AOSCM_Karalis/'+exper+'/'
dsprog=xr.open_dataset(my_input_dir + my_prog_suffix,decode_times=False)
dsdiag=xr.open_dataset(my_input_dir + my_diag_suffix,decode_times=False)
dsdiag2=xr.open_dataset(my_input_dir + my_diag2_suffix,decode_times=False)


output = xr.merge([dsprog, dsdiag, dsdiag2])
output.to_netcdf(my_input_dir + my_merged_suffix)
# specify Github scratch directory where processed model output will be committed (automate later)
my_output_filename =  'aoscm.nc'
my_gitdir = '../../output_scm/AOSCM/sandbox/'
# my_gitdir = 'dephy/'

ERROR 1: PROJ: proj_create_from_database: Open of /opt/conda/share/proj failed

Read ModelE3 data#

Read single file containing all output data#

Note: ERROR 1: PROJ… message can be ignored here.

input_filename = my_input_dir + my_merged_suffix


modele_data = xr.open_dataset(input_filename,decode_times=False)

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

# full parameter list

modele_data
# scm_in

do_ice =  True

<xarray.Dataset>
Dimensions:              (nlev: 137, nlevp1: 138, nlevs: 4, time: 161,
                          ntiles: 9, norg: 4, ncextr: 137)
Coordinates:
  * nlev                 (nlev) int32 1 2 3 4 5 6 7 ... 132 133 134 135 136 137
  * nlevp1               (nlevp1) int32 1 2 3 4 5 6 ... 133 134 135 136 137 138
  * nlevs                (nlevs) int32 1 2 3 4
  * time                 (time) float32 0.0 450.0 900.0 ... 7.155e+04 7.2e+04
  * ntiles               (ntiles) int32 -2147483647 -2147483647 ... -2147483647
  * norg                 (norg) int32 -2147483647 -2147483647 ... -2147483647
  * ncextr               (ncextr) int32 -2147483647 -2147483647 ... -2147483647
Data variables: (12/223)
    pressure_f           (time, nlev) float32 ...
    pressure_h           (time, nlevp1) float32 ...
    height_f             (time, nlev) float32 ...
    height_h             (time, nlevp1) float32 ...
    relative_humidity    (time, nlev) float32 ...
    t                    (time, nlev) float32 ...
    ...                   ...
    extra_col_27         (time, ncextr) float32 ...
    extra_col_28         (time, ncextr) float32 ...
    extra_col_29         (time, ncextr) float32 ...
    extra_col_30         (time, ncextr) float32 ...
    extra_col_31         (time, ncextr) float32 ...
    extra_col_32         (time, ncextr) float32 ...
Attributes:
    title:         SCM: 36r1  Sim: trref_wind
    modelID:       36r1
    simulationID:  trref_wind
    dataID:        SCM_OUTPUT
    start_day:     20200412
    start_hour:    79200

Interpolate to fixed pressure vertical grid#

for var in modele_data.variables:
    if len(modele_data[var].data.shape) >1:
        if modele_data[var].data.shape[1] == 137:
            for tt in range(modele_data[var].shape[0]):
                modele_data[var].data[tt,:] = np.interp(modele_data['pressure_f'].data[0,:], modele_data['pressure_f'].data[tt,:], modele_data[var].data[tt,:])
        elif modele_data[var].data.shape[1] == 138:
            for tt in range(modele_data[var].shape[0]):
                modele_data[var].data[tt,:] = np.interp(modele_data['pressure_h'].data[0,:], modele_data['pressure_h'].data[tt,:], modele_data[var].data[tt,:])

Calculate and append additional variables#

additional_var_list = ['surface_pressure', 'surface_friction_velocity',
                        'obukhov_length', 'inversion_height',
                        'mass_mixing_ratio_of_cloud_liquid_water_in_air',
                        'mass_mixing_ratio_of_ice_water_in_air',
                        'mass_mixing_ratio_of_rain_water_in_air',
                        'tendency_of_air_potential_temperature_due_to_radiative_heating',
                        'tendency_of_air_potential_temperature_due_to_microphysics',
                        'tendency_of_air_potential_temperature_due_to_mixing',
                        'atmosphere_mass_content_of_rain_water',
                        'air_dry_density', 
                        'surface_upwelling_longwave_flux',
                        'surface_sea_spray_number_flux',
                        'precipitation_flux_at_surface',
                        'precipitation_flux_of_ice_at_surface',
                        'vertical_flux_total_water',
                        'vertical_flux_water_vapor',
                        'zonal_momentum_flux',
                        'meridional_momentum_flux',
                        'water_vapor_mixing_ratio']


def additional_variables(var_name):

    scm_in = xr.open_dataset(my_input_dir + my_input_suffix)
    modele_data = xr.open_dataset(input_filename)
    
    Rair = 287.1
    Cp = 1003       
    kappa = 0.4     
    g = 9.8         #ms-2
    T_ice = 250.16  #K
    T_0 = 273.16    #K
    N_tot = 50      #cm-3 over ocean
    rho_liq = 1     #kgm-3
    P_0 = 100000    #Pa
    k = 2/7         
    
    ## surface pressure
    ps = scm_in.ps.data[0]*np.ones(len(modele_data.time.data))



    inv_height = np.zeros(len(modele_data.time.data))
    rain_wat_path = np.zeros(len(modele_data.time.data))

    
    ## dry air density
    rho = modele_data.pressure_f.data/(Rair*modele_data.t.data)

    ## rain water path
    for i in range(rain_wat_path.shape[0]):
        layer = modele_data.height_h.data[i,:-1] - modele_data.height_h.data[i,1:]
        rain_wat_path[i] = np.sum(modele_data.qr.data[i,:]*rho[i,:]*layer)

    ## liquid cloud fraction
    liq_cl_frct = np.zeros(modele_data.cloud_fraction.data.shape)
    for i in range(modele_data.cloud_fraction.data.shape[0]):
        ice_cld = (modele_data.t.data[i,:] <= T_ice) & (modele_data.cloud_fraction.data[i,:] > 0)
        mixed_cld = (modele_data.t.data[i,:] > T_ice) & (modele_data.t.data[i,:] < T_0) & (modele_data.cloud_fraction.data[i,:] > 0)
        warm_cld = (modele_data.t.data[i,:] > T_0) & (modele_data.cloud_fraction.data[i,:] > 0)
        liq_cl_frct[i,ice_cld] = 0
        liq_cl_frct[i,mixed_cld] = ((modele_data.t.data[i,mixed_cld] - T_ice)/(T_0 - T_ice))**2
        liq_cl_frct[i,warm_cld] = 1

    
    # upwelling longwave
    lwu = -(modele_data.sfc_lwrad.data - modele_data.sfc_lwrad_down.data)
    
    # surf. sea-spray
    sspray = np.zeros(modele_data.time.data.shape[0])
    
    # prec. fluxes
    prec = (modele_data.conv_rain.data + modele_data.stra_rain.data + 
            modele_data.conv_snow.data + modele_data.stra_snow.data )
    
    snow = modele_data.conv_snow.data + modele_data.stra_snow.data
    
    # turb. flux profiles
    turb_tot_water = (modele_data.turb_flx_wv.data + modele_data.turb_flx_liq.data + modele_data.turb_flx_ice.data)
    turb_tot_water_full = np.zeros(modele_data.t.data.shape)
    turb_tot_wv_full = np.zeros(modele_data.t.data.shape)
    for i in range(modele_data.time.data.shape[0]):
#         print(np.flipud(modele_data.height_f.data[i,:]))
        turb_tot_water_full[i,:] = -np.flipud(np.interp(np.flipud(modele_data.height_f.data[i,:]),
                                                       np.flipud(modele_data.height_h.data[i,:]),
                                                       np.flipud(turb_tot_water[i,:])))
        turb_tot_wv_full[i,:] = -np.flipud(np.interp(np.flipud(modele_data.height_f.data[i,:]),
                                                       np.flipud(modele_data.height_h.data[i,:]),
                                                       np.flipud(modele_data.turb_flx_wv.data[i,:])))

    
    
    #mixing ratios
    w = modele_data.q.data/(1-modele_data.q.data)
    w_l = modele_data.ql.data/(1-modele_data.ql.data)
    w_i = modele_data.qi.data/(1-modele_data.qi.data)
    w_r = modele_data.qr.data/(1-modele_data.qr.data)
    w_sn = modele_data.qsn.data/(1-modele_data.qsn.data)
    
    # momentum fluxes
    wu_flux = np.zeros(modele_data.t.data.shape)
    wv_flux = np.zeros(modele_data.t.data.shape)
    for i in range(modele_data.time.data.shape[0]):
        wu_flux[i,:] = -np.flipud(np.interp(np.flipud(modele_data.height_f.data[i,:]),
                                           np.flipud(modele_data.height_h.data[i,:]),
                                           np.flipud(modele_data.turb_flx_u.data[i,:])))
        wv_flux[i,:] = -np.flipud(np.interp(np.flipud(modele_data.height_f.data[i,:]),
                                           np.flipud(modele_data.height_h.data[i,:]),
                                           np.flipud(modele_data.turb_flx_v.data[i,:])))

    ## optical depths - According to IFS c43 parameterizations
    ## not computed for this experiment

#     r_e_liq = (3*modele_data.ql.data/(4*np.pi*rho_liq*0.8*N_tot))**(1/3)
#     [c_0, c_1, c_2, c_3] = [326.3, 12.42, 0.197, 0.0012]
#     r_e_ice = c_0 + c_1*T_c + c_2* T_c**2 + c_3*T_c**3
    
#     [a_0, a_1] = [0.02672, 1.32] #μmm2g-1
#     b_ext_liq = modele_data.ql.data*(a_0_liq + a_1_liq/r_e_liq)  # LWC in gm-3
#     [a_0, a_1] = [-9.45458*10**(-5), 2.52061] #μmm2g-1
#     b_ext_ice = modele_data.qi.data*(a_0_ice + a_1_ice/r_e_ice)
#     # same for rain and snow
#     b_all = b_ext_liq + b_ext_ice

#     delta = np.trapz(b_all,dp)

    ## temperature tendencies to theta tendencies
    pres_tend = np.zeros(modele_data.pressure_f.data.shape)
    
    for j in range(pres_tend.shape[1]):
        pres_tend[:,j] = np.gradient(modele_data.pressure_f.data[:,j], modele_data.time.data)
        
    theta_rad_dot = modele_data.extra_col_05.data*(P_0/modele_data.pressure_f.data)**(k) - k*modele_data.t.data*(P_0**k)*(modele_data.pressure_f.data**(-k-1))*pres_tend
    theta_turb_dot = modele_data.extra_col_08.data*(P_0/modele_data.pressure_f.data)**(k) - k*modele_data.t.data*(P_0**k)*(modele_data.pressure_f.data**(-k-1))*pres_tend
    theta_micro_dot = modele_data.extra_col_19.data*(P_0/modele_data.pressure_f.data)**(k) - k*modele_data.t.data*(P_0**k)*(modele_data.pressure_f.data**(-k-1))*pres_tend
    
    ## q tendencies
    q_turb_dot = modele_data.extra_col_09.data
    q_micro_dot = modele_data.extra_col_20.data
    
    
    ## inversion height
    t = modele_data.t.data
    t[modele_data.pressure_f.data/100<400] = np.nan
    for i in range(inv_height.shape[0]):
        t_grad = np.gradient(t[i,:],modele_data.height_f.data[i,:])
        i_grad_max = np.nonzero(t_grad==np.nanmax(t_grad))
        inv_height[i] = modele_data.height_f.data[i, i_grad_max[0]]


    ## u_star
    stress = np.sqrt(modele_data.u_sfc_strss.data**2 + modele_data.v_sfc_strss.data**2)
    u_star = np.sqrt(stress/rho[0,-1])

    ## L_star
    L_star = (u_star**3)/((kappa*g/modele_data.pot_temperature.data[:,-1])*(modele_data.sfc_sen_flx.data/(Cp*rho[:,-1])))


    if var_name == 'surface_pressure':
        return ps
    elif var_name == 'surface_friction_velocity':
        return u_star
    elif var_name == 'obukhov_length':
        return L_star
    elif var_name == 'inversion_height':
        return inv_height
    elif var_name == 'water_vapor_mixing_ratio':
        return w
    elif var_name == 'mass_mixing_ratio_of_cloud_liquid_water_in_air':
        return w_l
    elif var_name == 'mass_mixing_ratio_of_ice_water_in_air':
        return w_i + w_sn
    elif var_name == 'mass_mixing_ratio_of_rain_water_in_air':
        return w_r
    elif var_name == 'tendency_of_air_potential_temperature_due_to_radiative_heating':
        return theta_rad_dot
    elif var_name == 'tendency_of_air_potential_temperature_due_to_microphysics':
        return theta_micro_dot
    elif var_name == 'tendency_of_air_potential_temperature_due_to_mixing':
        return theta_turb_dot
    elif var_name == 'atmosphere_mass_content_of_rain_water':
        return rain_wat_path
    elif var_name == 'air_dry_density':
        return rho
    elif var_name == 'surface_upwelling_longwave_flux':
        return lwu
    elif var_name == 'surface_sea_spray_number_flux':
        return sspray
    elif var_name == 'precipitation_flux_at_surface':
        return prec
    elif var_name == 'vertical_flux_total_water':
        return turb_tot_water_full
    elif var_name == 'vertical_flux_water_vapor':
        return turb_tot_wv_full
    elif var_name == 'zonal_momentum_flux':
        return wu_flux
    elif var_name == 'meridional_momentum_flux':
        return wv_flux

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 ModelE3 variables to requested outputs#

Expand the table to include columns that indicate ModelE3 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 ModelE3 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=='surface_temperature':                                     #done
        vars_mean_list.model_name.iat[index] = 't_skin'


    if standard_name=='surface_roughness_length_for_momentum_in_air':             #done
        vars_mean_list.model_name.iat[index] = 'rough_len_mom'
    if standard_name=='surface_roughness_length_for_heat_in_air':                 #done
        vars_mean_list.model_name.iat[index] = 'rough_len_heat'
        
    if standard_name=='surface_roughness_length_for_humidity_in_air':           #to_do 
        vars_mean_list.model_name.iat[index] = 'rough_len_heat'

    if standard_name=='surface_upward_sensible_heat_flux':                        #done
        vars_mean_list.model_name.iat[index] = 'sfc_sen_flx'
        vars_mean_list.conv_factor.iat[index] = -1.
    if standard_name=='surface_upward_latent_heat_flux':                          #done
        vars_mean_list.model_name.iat[index] = 'sfc_lat_flx' 
        vars_mean_list.conv_factor.iat[index] = -1.
        
        
    if standard_name=='pbl_height':                                               #done
        vars_mean_list.model_name.iat[index] = 'pbl_height'     


    if standard_name=='atmosphere_mass_content_of_liquid_cloud_water':            #done
        vars_mean_list.model_name.iat[index] = 'liq_wat_path'                     
        vars_mean_list.conv_factor.iat[index] = 1
        

    if do_ice:
        if standard_name=='atmosphere_mass_content_of_ice_water':                  #done
            vars_mean_list.model_name.iat[index] = 'ice_wat_path'
            vars_mean_list.conv_factor.iat[index] = 1
    if standard_name=='area_fraction_cover_of_hydrometeors':                        ##done
        vars_mean_list.model_name.iat[index] = 'total_cloud'
        
        
    if standard_name=='toa_outgoing_longwave_flux':                                #done
        vars_mean_list.model_name.iat[index] = 'top_lwrad'
        vars_mean_list.conv_factor.iat[index] = -1
    if standard_name=='surface_downwelling_longwave_flux':                         #done
        vars_mean_list.model_name.iat[index] = 'sfc_lwrad_down'  
        vars_mean_list.conv_factor.iat[index] = 1

    if standard_name=='height':                                                     #done
        vars_mean_list.model_name.iat[index] = 'height_f'          

    if standard_name=='eastward_wind':                                              #done
        vars_mean_list.model_name.iat[index] = 'u'
    if standard_name=='northward_wind':                                             #done
        vars_mean_list.model_name.iat[index] = 'v'

    if standard_name=='air_temperature':                                             #done
        vars_mean_list.model_name.iat[index] = 't'

    if standard_name=='relative_humidity':                                            #done
        vars_mean_list.model_name.iat[index] = 'relative_humidity'
        vars_mean_list.conv_factor.iat[index] = 1


    if standard_name=='air_potential_temperature':                                   #done
        vars_mean_list.model_name.iat[index] = 'pot_temperature'
        
        


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

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”.

# create DEPHY output file
dephy_filename = my_gitdir + my_output_suffix
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='AOSCM results for COMBLE-MIP case: fixed stratiform Nd'
dephy_file.reference='https://github.com/ARM-Development/comble-mip'
dephy_file.authors='Michail Karalis (michail.karalis@misu.su.se), Gunilla Svensson (gunilla@misu.su.se)'
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_AOSCM_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'
dephy_file.np=modele_data.sizes['nlev']
dephy_file.alt='Nd = 50 cm-3 over the ocean and diagnostic Ni'

# create dimensions
nt = modele_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
# time1 = dt.datetime.strptime(str(modele_data['time'].data[0]),'%Y-%m-%d %H:%M:%S')
# time2 = dt.datetime.strptime(str(modele_data['time'].data[1]),'%Y-%m-%d %H:%M:%S')
# delta_t = (time2-time1).total_seconds()


delta_t = modele_data['time'].data[1]
# time[:] = (np.arange(nt)+1.)*delta_t
time[:] = modele_data['time'].data

nl = modele_data.dims['nlev']
# pa = dephy_file.createDimension('pressure', nl)
# pa = dephy_file.createVariable('pressure', np.float64, ('pressure',))
# pa.units = 'Pa'
# pa.long_name = 'pressure'
# pa[:] = modele_data['pressure_f'].data[0,:]*100.

lvl = dephy_file.createDimension('pa', nl)
pa = dephy_file.createVariable('pa', np.float64, ('pa',))
pa.units = 'Pa'
pa.long_name = 'pressure'
pa[:] = modele_data['pressure_f'].data[0,:]


# create and fill variables
for index in vars_mean_list.index[3:]:
    std_name = vars_mean_list.standard_name.iat[index]
    print(std_name)
#   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
        if vars_mean_list.model_name.iat[index]!='missing data':
            new_sca[:] = modele_data[mod_name].data*c_factor
        elif np.any(np.array(additional_var_list) == std_name) :
            new_sca[:] = additional_variables(std_name)
    if vars_mean_list.dimensions.iat[index]=='time, layer':
        new_snd = dephy_file.createVariable(var_name, np.float64, ('time','pa'))
        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': 
            new_snd[:] = modele_data[mod_name].data*c_factor
        elif np.any(np.array(additional_var_list) == std_name) :
            new_snd[:] = additional_variables(std_name)
    


dephy_file.close()

The file ../../output_scm/AOSCM/sandbox/AOSCM_FixN_def_z0_alt.nc has been deleted successfully
layer_top_pressure
surface_pressure
surface_temperature
surface_friction_velocity
surface_roughness_length_for_momentum_in_air
surface_roughness_length_for_heat_in_air
surface_roughness_length_for_humidity_in_air
surface_upward_sensible_heat_flux
surface_upward_latent_heat_flux
obukhov_length
pbl_height
inversion_height
atmosphere_mass_content_of_liquid_cloud_water
atmosphere_mass_content_of_rain_water
atmosphere_mass_content_of_ice_water
area_fraction_cover_of_hydrometeors
area_fraction_cover_of_liquid_cloud
area_fraction_cover_of_convective_hydrometeors
optical_depth
optical_depth_of_liquid_cloud
precipitation_flux_at_surface
precipitation_flux_of_ice_at_surface
toa_outgoing_longwave_flux
surface_downwelling_longwave_flux
surface_upwelling_longwave_flux
surface_sea_spray_number_flux
height
eastward_wind
northward_wind
air_dry_density
air_temperature
water_vapor_mixing_ratio
relative_humidity
relative_humidity_over_ice
air_potential_temperature
mass_mixing_ratio_of_cloud_liquid_water_in_air
mass_mixing_ratio_of_rain_water_in_air
mass_mixing_ratio_of_ice_water_in_air
area_fraction_of_hydrometeors
area_fraction_of_liquid_cloud
area_fraction_of_convective_hydrometeors
precipitation_flux_in_air
precipitation_flux_in_air_in_ice_phase
specific_turbulent_kinetic_energy
dissipation_rate_of_turbulent_kinetic_energy
zonal_momentum_flux
meridional_momentum_flux
variance_of_upward_air_velocity
vertical_flux_potential_temperature
vertical_flux_liquid_ice_water_potential_temperature
vertical_flux_water_vapor
vertical_flux_total_water
convection_updraft_mass_flux
convection_downdraft_mass_flux
downwelling_longwave_flux_in_air
upwelling_longwave_flux_in_air
tendency_of_air_potential_temperature_due_to_radiative_heating
tendency_of_air_potential_temperature_due_to_microphysics
tendency_of_air_potential_temperature_due_to_mixing
tendency_of_water_vapor_mixing_ratio_due_to_microphysics
tendency_of_water_vapor_mixing_ratio_due_to_mixing
number_of_total_aerosol_mode1
number_of_total_aerosol_mode2
number_of_total_aerosol_mode3
tendency_of_aerosol_number_due_to_warm_microphysics
tendency_of_aerosol_number_due_to_cold_microphysics
tendency_of_aerosol_number_due_to_mixing
tendency_of_ice_number_due_to_heterogeneous_freezing
tendency_of_ice_number_due_to_secondary_ice_production
tendency_of_ice_number_due_to_homogeneous_freezing
mass_mixing_ratio_of_liquid_cloud_water_in_air_stratiform
mass_mixing_ratio_of_rain_water_in_air_stratiform
mass_mixing_ratio_of_ice_cloud_in_air_stratiform
mass_mixing_ratio_of_ice_precipitation_in_air_stratiform
mass_mixing_ratio_of_liquid_cloud_water_in_air_convective
mass_mixing_ratio_of_rain_water_in_air_convective
mass_mixing_ratio_of_ice_cloud_in_air_convective
mass_mixing_ratio_of_ice_precipitation_in_air_convective
number_of_liquid_cloud_droplets_in_air_stratiform
number_of_rain_drops_in_air_stratiform
number_of_ice_cloud_crystals_in_air_stratiform
number_of_ice_precipitation_crystals_in_air_stratiform
effective_radius_of_liquid_cloud_droplets_convective
effective_radius_of_rain_convective
effective_radius_of_ice_cloud_convective
effective_radius_of_ice_precipitation_convective
area_fraction_of_liquid_cloud_stratiform
area_fraction_of_rain_stratiform
area_fraction_of_ice_cloud_stratiform
area_fraction_of_ice_precipitation_stratiform
area_fraction_of_liquid_cloud_convective
area_fraction_of_rain_convective
area_fraction_of_ice_cloud_convective
area_fraction_of_ice_precipitation_convective
mass_weighted_fall_speed_of_liquid_cloud_water_stratiform
mass_weighted_fall_speed_of_rain_stratiform
mass_weighted_fall_speed_of_ice_cloud_stratiform
mass_weighted_fall_speed_of_ice_precipitation_stratiform
mass_weighted_fall_speed_of_liquid_cloud_water_convective
mass_weighted_fall_speed_of_rain_convective
mass_weighted_fall_speed_of_cloud_ice_crystals_convective
mass_weighted_fall_speed_of_ice_precipitation_convective

Check output file#

dephy_check = xr.open_dataset(dephy_filename)

dephy_check

<xarray.Dataset>
Dimensions:         (time: 161, pa: 137)
Coordinates:
  * time            (time) datetime64[ns] 2020-03-12T22:00:00 ... 2020-03-13T...
  * pa              (pa) float64 1.0 2.551 3.884 ... 9.918e+04 9.943e+04
Data variables: (12/102)
    pe              (time, pa) float64 ...
    ps              (time) float64 ...
    ts              (time) float64 ...
    ustar           (time) float64 ...
    z0              (time) float64 ...
    z0h             (time) float64 ...
    ...              ...
    vmics           (time, pa) float64 ...
    vmips           (time, pa) float64 ...
    vmlcc           (time, pa) float64 ...
    vmlrc           (time, pa) float64 ...
    vmicc           (time, pa) float64 ...
    vmipc           (time, pa) float64 ...
Attributes: (12/15)
    title:           AOSCM results for COMBLE-MIP case: fixed stratiform Nd
    reference:       https://github.com/ARM-Development/comble-mip
    authors:         Michail Karalis (michail.karalis@misu.su.se), Gunilla Sv...
    source:          ~/AOSCM_Karalis/FixN_def_z0/aoscm.nc
    version:         2024-04-15 10:33:58
    format_version:  DEPHY SCM format version 1.6
    ...              ...
    surfaceType:     ocean
    surfaceForcing:  ts
    lat:             74.5 deg N
    dp:              see pressure variable
    np:              137
    alt:             Nd = 50 cm-3 over the ocean and diagnostic Ni