Example: convert DHARMA LES output to DEPHY format#

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

Contributed by Ann Fridlind from NASA/GISS

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_les/YOUR_MODEL_NAME/sandbox/YOUR_RUN_NAME. These can be committed and removed at any time.

It is requested to name your baseline small-domain run directory as ‘Lx25km_dx100m’ (in place of YOUR_RUN_NAME), so that it can readily be automatically compared with other runs using the same test specification.

# specify start time of simulation and simulation name
start_dtime = '2020-03-12 22:00:00.0'

# specify input data directory name, traceable to source machine, and
# specify output file name (see file naming convention in TOC)

# FixN with no ice test
my_readdir = 'case0313_diag_ice0_nomiz_dx100_specZ0'
my_outfile = 'DHARMA_Lx25_dx100_FixN_noice.nc'

# FixN with ice test
my_readdir = 'case0313_diag_ice25_nomiz_dx100_specZ0_sbhomo'
my_outfile = 'DHARMA_Lx25_dx100_FixN.nc'

# ProgNa with ice test
my_readdir = 'case0313_diag_ice25_nomiz_dx100_specZ0_hm_largedomain'
my_outfile = 'DHARMA_Lx125_dx100_FixN.nc'

# FixN with ice test and default z0
#my_readdir = 'case0313_diag_ice25_nomiz_dx100' 
#my_outfile = 'DHARMA_Lx25_dx100_FixN_def_z0.nc'

# ProgNa with ice test
#my_readdir = 'case0313_prog_ice25_nomiz_dx100_specZ0'
#my_outfile = 'DHARMA_Lx25_dx100_ProgNa.nc'

# specify local source directories (additional subdirectories if restart was required)

my_rundir = '/data/home/floriantornow/dharma/sandbox/' + my_readdir + '/'

my_outdirs = sorted([f for f in os.listdir(my_rundir) if not f.startswith('.')], key=str.lower)
print(my_outdirs)

# specify Github scratch directory where processed model output will be committed
my_gitdir = '../../output_les/dharma/sandbox/'
['0-6h', '6-20h']

Read DHARMA data#

Read set-up parameters#

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

# read in DHARMA parameter settings
input_filename = my_rundir + my_outdirs[0] + '/dharma.cdf'
dharma_params = xr.open_dataset(input_filename)
print(input_filename)

# check if the run contains ice variables
do_ice = bool(dharma_params['Cond'].do_ice)
print('do_ice = ',do_ice)

# check for prognostic aerosol
do_progNa = bool(dharma_params['Cond'].do_prog_na)
print('do_progNa = ',do_progNa)

# full parameter list
dharma_params
/data/home/floriantornow/dharma/sandbox/case0313_diag_ice25_nomiz_dx100_specZ0_hm_largedomain/0-6h/dharma.cdf
do_ice =  True
do_progNa =  False
<xarray.Dataset> Size: 104B
Dimensions:      ()
Data variables: (12/26)
    geometry     int32 4B ...
    timing       int32 4B ...
    options      int32 4B ...
    zstretch     int32 4B ...
    assimilate   int32 4B ...
    restart      int32 4B ...
    ...           ...
    Radiation    int32 4B ...
    Surface      int32 4B ...
    Subsidence   int32 4B ...
    Turb         int32 4B ...
    InvAnalysis  int32 4B ...
    carma        int32 4B ...

Read domain-mean profiles#

# concatenate DHARMA domain-mean instantaneous profiles and take 10-min average:
# resample-average before concatenating and removing duplicates
for index, elem in enumerate(my_outdirs):
    input_filename = my_rundir + elem + '/dharma.soundings.cdf'
    print(input_filename)
    if index==0:
        dharma_snds = xr.open_dataset(input_filename)
        dharma_snds['time'] = pd.to_datetime(dharma_snds['time'].values, unit='s', origin=pd.Timestamp(start_dtime))
        dharma_snds = dharma_snds.resample(time="600s",closed="right",label="right").mean()
    else:
        dharma_snds2 = xr.open_dataset(input_filename)
        dharma_snds2['time'] = pd.to_datetime(dharma_snds2['time'].values, unit='s', origin=pd.Timestamp(start_dtime))
        dharma_snds2 = dharma_snds2.resample(time="600s",closed="right",label="right").mean()
        dharma_snds = xr.concat([dharma_snds,dharma_snds2],dim='time')
dharma_snds = dharma_snds.drop_duplicates('time',keep='first')
dharma_snds
/data/home/floriantornow/dharma/sandbox/case0313_diag_ice25_nomiz_dx100_specZ0_hm_largedomain/0-6h/dharma.soundings.cdf
/data/home/floriantornow/dharma/sandbox/case0313_diag_ice25_nomiz_dx100_specZ0_hm_largedomain/6-20h/dharma.soundings.cdf
<xarray.Dataset> Size: 50MB
Dimensions:       (zt: 159, zw: 160, time: 121)
Coordinates:
  * zt            (zt) float64 1kB 10.0 32.5 60.0 ... 6.92e+03 6.975e+03
  * zw            (zw) float64 1kB 0.0 20.0 45.0 ... 6.89e+03 6.95e+03 7e+03
  * time          (time) datetime64[ns] 968B 2020-03-12T22:00:00 ... 2020-03-...
Data variables: (12/325)
    jact          (time, zt) float64 154kB 0.4543 0.5679 0.6814 ... 1.363 1.136
    jacw          (time, zw) float64 155kB 0.3975 0.5111 0.6246 ... 1.249 1.704
    rhobar        (time, zt) float64 154kB 1.399 1.396 1.392 ... 0.592 0.5882
    u             (time, zt) float64 154kB 2.775 2.799 2.844 ... 4.846 4.877
    u2            (time, zt) float64 154kB 3.852 3.917 4.043 ... 11.74 11.89
    Su_rk         (time, zt) float64 154kB 0.0 0.0 0.0 ... 0.0003637 0.0003731
    ...            ...
    hw            (time, zw) float64 155kB 0.0 0.0 0.0 ... -6.454e-09 0.0
    qw            (time, zw) float64 155kB 0.0 0.0 0.0 ... 1.259e-11 0.0
    txz_tot       (time, zw) float64 155kB -0.1378 -0.1055 ... -3.064e-07 0.0
    tyz_tot       (time, zw) float64 155kB 1.335 1.019 0.4581 ... -1.823e-06 0.0
    qhz_tot       (time, zw) float64 155kB -8.216e-05 -6.272e-05 ... 0.0
    qqz_tot       (time, zw) float64 155kB 5.931e-06 4.528e-06 ... 1.259e-11 0.0
Attributes:
    nx:        1280
    ny:        1280
    nz:        159
    L_x:       128000.0
    L_y:       128000.0
    H:         7000.0
    theta_00:  247.957

Profile unit conversions and sundry#

# create a dummy sounding and initialize some variables needed
dummy_snd = dharma_snds['qc']*0.
nz = dharma_params['geometry'].nz
dz = dharma_snds['zw'].data[1:nz+1]-dharma_snds['zw'].data[0:nz]
cp = 1004.
Lhe = 2.50e6
Lhs = Lhe + 3.34e5

# compute some intermediate quantities for use below
Fql_turb = dharma_snds['Fqc_turb'].data+dharma_snds['Fqr_turb'].data
if do_ice:
    Fqi_turb = dharma_snds['Fqic_turb'].data+dharma_snds['Fqif_turb'].data+dharma_snds['Fqid_turb'].data
wql_tot = 0.5*(Fql_turb[:,0:nz]+Fql_turb[:,1:nz+1])/dharma_snds['rhobar'].data+dharma_snds['WQL'].data

if do_ice:
    wqi_tot = 0.5*(Fqi_turb[:,0:nz]+Fqi_turb[:,1:nz+1])/dharma_snds['rhobar'].data+dharma_snds['WQI'].data
else:
    wqi_tot = np.nan*wql_tot
PFql = dharma_snds['PFqc'].data+dharma_snds['PFqr'].data
if do_ice:
    PFqi = dharma_snds['PFqic'].data+dharma_snds['PFqif'].data+dharma_snds['PFqid'].data
else:
    PFqi = np.nan*PFql
wpfl = 0.5*(PFql[:,0:nz]+PFql[:,1:nz+1])*-1./3600./dharma_snds['rhobar'].data
wpfi = 0.5*(PFqi[:,0:nz]+PFqi[:,1:nz+1])*-1./3600./dharma_snds['rhobar'].data
PFqc = dharma_snds['PFqc'].data 
PFqr = dharma_snds['PFqr'].data 
if do_ice:
    PFqic = dharma_snds['PFqic'].data 
    PFqif = dharma_snds['PFqif'].data 
    PFqid = dharma_snds['PFqid'].data 
if do_progNa:
    ssa_sfc = (dharma_snds['Sna_1_sfc'].data[:,0]+dharma_snds['Sna_2_sfc'].data[:,0]+dharma_snds['Sna_3_sfc'].data[:,0])*dharma_snds['zw'].data[1]
Flwd = dharma_snds['Flw_dn'].data
Flwu = dharma_snds['Flw_up'].data
Fnlw = Flwu - Flwd
Suvar_adv = dharma_snds['Su2_adv'].data - dharma_snds['Subar2_adv'].data
Svvar_adv = dharma_snds['Sv2_adv'].data - dharma_snds['Svbar2_adv'].data
Swvar_adv = dharma_snds['Sw2_adv'].data - dharma_snds['Swbar2_adv'].data
Stke_a = Suvar_adv + Svvar_adv + 0.5*(Swvar_adv[:,0:nz]+Swvar_adv[:,1:nz+1])
Stke_adv_dis = dharma_snds['Stke_adv'].data + dharma_snds['Sprod'].data
Smke = (dharma_snds['u'].data-dharma_params.translate.u)*dharma_snds['Suavg_SGS'].data + (dharma_snds['v'].data-dharma_params.translate.v)*dharma_snds['Svavg_SGS'].data
Ske = dharma_snds['Su2avg_SGS'].data+dharma_snds['Sv2avg_SGS'].data+dharma_snds['Sw2avg_SGS'].data
Stke_dis = Smke - Ske

# append new variables to the data structure
dharma_snds = dharma_snds.assign(theta = dummy_snd + (dharma_snds['th'].data+1)*dharma_snds.theta_00)
dharma_snds = dharma_snds.assign(pi = dummy_snd + dharma_snds['T'].data/dharma_snds['theta'].data)
dharma_snds = dharma_snds.assign(pressure = dummy_snd + np.power(dharma_snds['pi'].data,7./2)*np.power(10.,5))
dharma_snds = dharma_snds.assign(PF = dummy_snd + 0.5*(PFqc[:,0:nz]+PFqc[:,1:nz+1]) + 0.5*(PFqr[:,0:nz]+PFqr[:,1:nz+1]))
dharma_snds = dharma_snds.assign(nlcic = dummy_snd + dharma_snds['nc_cld'].data*1.e6/dharma_snds['rhobar'].data)
if do_ice:
    dharma_snds = dharma_snds.assign(niic = dummy_snd + dharma_snds['ni_cld'].data*1.e6/dharma_snds['rhobar'].data)
if do_ice:
    dharma_snds = dharma_snds.assign(PFi = dummy_snd + 0.5*(PFqic[:,0:nz]+PFqic[:,1:nz+1]) + 0.5*(PFqif[:,0:nz]+PFqif[:,1:nz+1]) 
                                     + 0.5*(PFqid[:,0:nz]+PFqid[:,1:nz+1]) )
    dharma_snds['PF'] += dharma_snds['PFi']
else:
    dharma_snds['RHI'] = np.nan*dharma_snds['RH']
    dharma_snds['PFi'] = np.nan*dharma_snds['PF']
dharma_snds = dharma_snds.assign(uw_zt = dummy_snd + 0.5*(dharma_snds['txz_tot'].data[:,0:nz]+dharma_snds['txz_tot'].data[:,1:nz+1]))
dharma_snds = dharma_snds.assign(vw_zt = dummy_snd + 0.5*(dharma_snds['tyz_tot'].data[:,0:nz]+dharma_snds['tyz_tot'].data[:,1:nz+1]))
dharma_snds = dharma_snds.assign(w2_zt = dummy_snd + (dharma_snds['w2'].data[:,0:nz]+dharma_snds['w2'].data[:,1:nz+1])) # w2 = 0.5*w'2
dharma_snds = dharma_snds.assign(wth_zt = dummy_snd + 0.5*(dharma_snds['qhz_tot'].data[:,0:nz] + 
                                        dharma_snds['qhz_tot'].data[:,1:nz+1])*dharma_snds.theta_00)
dharma_snds = dharma_snds.assign(wthli_zt = dummy_snd + dharma_snds['wth_zt'].data
                    - wql_tot*Lhe/(dharma_snds['pi'].data*cp) - wqi_tot*Lhs/(dharma_snds['pi'].data*cp)
                    - wpfl*Lhe/(dharma_snds['pi'].data*cp) - wpfi*Lhs/(dharma_snds['pi'].data*cp))                                 
dharma_snds = dharma_snds.assign(wqv_zt = dummy_snd + 0.5*(dharma_snds['qqz_tot'].data[:,0:nz]+dharma_snds['qqz_tot'].data[:,1:nz+1]))
dharma_snds = dharma_snds.assign(wqt_zt = dummy_snd + dharma_snds['wqv_zt'].data + wql_tot + wqi_tot + wpfl + wpfi)
dharma_snds = dharma_snds.assign(eps = dummy_snd + Stke_dis + Stke_adv_dis)
dharma_snds = dharma_snds.assign(LWdn = dummy_snd + 0.5*(Flwd[:,0:nz]+Flwd[:,1:nz+1]))
dharma_snds = dharma_snds.assign(LWup = dummy_snd + 0.5*(Flwu[:,0:nz]+Flwu[:,1:nz+1]))
dharma_snds = dharma_snds.assign(HRlw = dummy_snd + 0.5*(Fnlw[:,0:nz]+Fnlw[:,1:nz+1])/dz/dharma_snds['rhobar'].data)
dharma_snds = dharma_snds.assign(dth_micro = dummy_snd + dharma_snds['Sth_micro'].data + dharma_snds['Sth_cond'].data)
dharma_snds = dharma_snds.assign(dq_micro = dummy_snd + dharma_snds['Sqv_micro'].data + dharma_snds['Sqv_cond'].data)
#dharma_snds = dharma_snds.assign(dth_turb = dummy_snd + (dharma_snds['qhz_tot'].data[:,0:nz] - 
#                                        dharma_snds['qhz_tot'].data[:,1:nz+1])*dharma_snds.theta_00)
#dharma_snds = dharma_snds.assign(dq_turb = dummy_snd + dharma_snds['qqz_tot'].data[:,0:nz] - dharma_snds['qqz_tot'].data[:,1:nz+1])
dharma_snds = dharma_snds.assign(dth_turb = dummy_snd + (dharma_snds['qhz_tot'].data[:,0:nz] - 
                                                         dharma_snds['qhz_tot'].data[:,1:nz+1])/dz)
dharma_snds = dharma_snds.assign(dq_turb = dummy_snd + (dharma_snds['qqz_tot'].data[:,0:nz] - dharma_snds['qqz_tot'].data[:,1:nz+1])/dz)

 
if do_progNa:
    dharma_snds = dharma_snds.assign(na_loss_liq = dummy_snd + dharma_snds['na_loss_prof'].data - dharma_snds['na_loss_ice'].data)
    dharma_snds = dharma_snds.assign(dna_mixing = dummy_snd + dharma_snds['Sna_1_adv'].data + dharma_snds['Sna_2_adv'].data + dharma_snds['Sna_3_adv'].data + 
                                dharma_snds['Sna_1_sfc'].data + dharma_snds['Sna_2_sfc'].data + dharma_snds['Sna_3_sfc'].data)

dharma_snds
<xarray.Dataset> Size: 53MB
Dimensions:       (zt: 159, zw: 160, time: 121)
Coordinates:
  * zt            (zt) float64 1kB 10.0 32.5 60.0 ... 6.92e+03 6.975e+03
  * zw            (zw) float64 1kB 0.0 20.0 45.0 ... 6.89e+03 6.95e+03 7e+03
  * time          (time) datetime64[ns] 968B 2020-03-12T22:00:00 ... 2020-03-...
Data variables: (12/347)
    jact          (time, zt) float64 154kB 0.4543 0.5679 0.6814 ... 1.363 1.136
    jacw          (time, zw) float64 155kB 0.3975 0.5111 0.6246 ... 1.249 1.704
    rhobar        (time, zt) float64 154kB 1.399 1.396 1.392 ... 0.592 0.5882
    u             (time, zt) float64 154kB 2.775 2.799 2.844 ... 4.846 4.877
    u2            (time, zt) float64 154kB 3.852 3.917 4.043 ... 11.74 11.89
    Su_rk         (time, zt) float64 154kB 0.0 0.0 0.0 ... 0.0003637 0.0003731
    ...            ...
    LWup          (time, zt) float64 154kB 211.1 211.2 211.2 ... 196.7 196.4
    HRlw          (time, zt) float64 154kB 3.237 2.605 2.186 ... 4.354 5.262
    dth_micro     (time, zt) float64 154kB 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    dq_micro      (time, zt) float64 154kB 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    dth_turb      (time, zt) float64 154kB -9.717e-07 -1.381e-06 ... -1.291e-10
    dq_turb       (time, zt) float64 154kB 7.016e-08 9.97e-08 ... 2.518e-13
Attributes:
    nx:        1280
    ny:        1280
    nz:        159
    L_x:       128000.0
    L_y:       128000.0
    H:         7000.0
    theta_00:  247.957

Read domain-mean scalars#

for index, elem in enumerate(my_outdirs):
    input_filename = my_rundir + elem + '/dharma.scalars.cdf'
    print(input_filename)
    if index==0:
        dharma_scas = xr.open_dataset(input_filename)
        dharma_scas['time'] = pd.to_datetime(dharma_scas['time'].values, unit='s', origin=pd.Timestamp(start_dtime))
        dharma_scas = dharma_scas.resample(time="600s",closed="right",label="right").mean()
    else:
        dharma_scas2 = xr.open_dataset(input_filename)
        dharma_scas2['time'] = pd.to_datetime(dharma_scas2['time'].values, unit='s', origin=pd.Timestamp(start_dtime))
        dharma_scas2 = dharma_scas2.resample(time="600s",closed="right",label="right").mean()
        dharma_scas = xr.concat([dharma_scas,dharma_scas2],dim='time')
dharma_scas = dharma_scas.drop_duplicates('time',keep='first')
dharma_scas
/data/home/floriantornow/dharma/sandbox/case0313_diag_ice25_nomiz_dx100_specZ0_hm_largedomain/0-6h/dharma.scalars.cdf
/data/home/floriantornow/dharma/sandbox/case0313_diag_ice25_nomiz_dx100_specZ0_hm_largedomain/6-20h/dharma.scalars.cdf
<xarray.Dataset> Size: 110kB
Dimensions:                     (time: 121)
Coordinates:
  * time                        (time) datetime64[ns] 968B 2020-03-12T22:00:0...
Data variables: (12/113)
    iwp                         (time) float64 968B 0.0 0.0 0.0 ... 380.4 368.9
    rms_iwp                     (time) float64 968B 0.0 0.0 0.0 ... 689.2 678.3
    iwpf                        (time) float64 968B 0.0 0.0 0.0 ... 287.1 278.7
    iwpd                        (time) float64 968B 0.0 0.0 0.0 ... 50.22 46.96
    iwpc                        (time) float64 968B 0.0 0.0 0.0 ... 43.03 43.2
    max_qi                      (time) float64 968B 0.0 0.0 ... 0.002166
    ...                          ...
    predictorstep_call          (time) float64 968B 0.849 134.4 ... 2.457e+04
    setfirstphi_call            (time) float64 968B 0.8492 0.8492 ... 0.0 0.0
    forwardeuler_call           (time) float64 968B 0.01256 1.697 ... 280.2
    compsadvect_call            (time) float64 968B 0.0 49.12 ... 8.254e+03
    processsources_call         (time) float64 968B 9.094 165.0 ... 3.026e+04
    computepressureupdate_call  (time) float64 968B 0.2887 30.88 ... 4.952e+03
# calculate some additional variables requested
dummy_sca = dharma_scas['lwp']*0.
dharma_scas = dharma_scas.assign(Psurf = dummy_sca + dharma_params['sounding'].Psurf*100.)
if do_ice:
    dharma_scas = dharma_scas.assign(opd_tot = dummy_sca + dharma_scas['opd_drops'].data + dharma_scas['opd_ice'].data)
else:
    dharma_scas = dharma_scas.assign(opd_tot = dummy_sca + dharma_scas['opd_drops'].data)
    #dharma_scas = dharma_scas.assign(RHI = np.nan*dharma_scas['opd_tot'])
dharma_scas = dharma_scas.assign(LWdnSFC = dummy_sca + dharma_snds['Flw_dn'].data[:,0])
dharma_scas = dharma_scas.assign(LWupSFC = dummy_sca + dharma_snds['Flw_up'].data[:,0])
dharma_scas = dharma_scas.assign(avg_precip_ice = dummy_sca + dharma_scas['avg_precip'].data 
                                 - dharma_snds['PFqc'].data[:,0] - dharma_snds['PFqr'].data[:,0])
if do_progNa:
    dharma_scas = dharma_scas.assign(ssaf = dummy_sca + ssa_sfc)
dharma_scas
<xarray.Dataset> Size: 115kB
Dimensions:                     (time: 121)
Coordinates:
  * time                        (time) datetime64[ns] 968B 2020-03-12T22:00:0...
Data variables: (12/118)
    iwp                         (time) float64 968B 0.0 0.0 0.0 ... 380.4 368.9
    rms_iwp                     (time) float64 968B 0.0 0.0 0.0 ... 689.2 678.3
    iwpf                        (time) float64 968B 0.0 0.0 0.0 ... 287.1 278.7
    iwpd                        (time) float64 968B 0.0 0.0 0.0 ... 50.22 46.96
    iwpc                        (time) float64 968B 0.0 0.0 0.0 ... 43.03 43.2
    max_qi                      (time) float64 968B 0.0 0.0 ... 0.002166
    ...                          ...
    computepressureupdate_call  (time) float64 968B 0.2887 30.88 ... 4.952e+03
    Psurf                       (time) float64 968B 9.954e+04 ... 9.954e+04
    opd_tot                     (time) float64 968B 0.0 0.0 0.0 ... 5.99 5.835
    LWdnSFC                     (time) float64 968B 120.6 120.6 ... 252.0 251.4
    LWupSFC                     (time) float64 968B 211.0 211.0 ... 339.7 338.4
    avg_precip_ice              (time) float64 968B 0.0 0.0 ... 0.3631 0.3323

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=0&format=xlsx',
                              sheet_name='Mean')

pd.set_option('display.max_rows', None)
vars_mean_list
standard_name variable_id units dimensions comment (10-min average reported at endpoints, green=minimum)
0 time time s dimension, seconds since 2020-03-12 18:00:00
1 height zf m dimension, altitude of mid-level points above ...
2 layer_top_height ze m dimension, altitude of layer top points above ...
3 surface_pressure ps Pa time
4 surface_temperature ts K time
5 surface_friction_velocity ustar m s-1 time
6 surface_roughness_length_for_momentum_in_air z0 m time
7 surface_roughness_length_for_heat_in_air z0h m time
8 surface_roughness_length_for_humidity_in_air z0q m time
9 surface_upward_sensible_heat_flux hfss W m-2 time
10 surface_upward_latent_heat_flux hfls W m-2 time
11 obukhov_length ol m time
12 atmosphere_mass_content_of_liquid_cloud_water lwpc kg m-2 time default breakdown to cloud water and rain in a...
13 atmosphere_mass_content_of_rain_water lwpr kg m-2 time
14 atmosphere_mass_content_of_ice_water iwp kg m-2 time all ice-phase hydrometeors
15 cloud_area_fraction clt 1 time fraction of atmospheric columns with total hyd...
16 optical_depth od 1 time scene (all sky); mid-visible, all hydrometeors...
17 optical_depth_of_liquid_cloud odlc 1 time as for optical_depth but cloud droplets only
18 precipitation_flux_at_surface pr kg m-2 s-1 time liquid and ice phase, all hydrometeors
19 precipitation_flux_at_surface_in_ice_phase pri kg m-2 s-1 time ice phase hydrometeors only
20 toa_outgoing_longwave_flux rlut W m-2 time at top of atmosphere
21 surface_downwelling_longwave_flux rlds W m-2 time
22 surface_upwelling_longwave_flux rlus W m-2 time
23 surface_sea_spray_number_flux ssaf m-2 s-1 time when using prognostic aerosol; emission only (...
24 air_pressure pa Pa time, height
25 eastward_wind ua m s-1 time, height
26 northward_wind va m s-1 time, height
27 air_dry_density rhoa kg m-3 time, height
28 air_temperature ta K time, height
29 water_vapor_mixing_ratio qv kg kg-1 time, height per kg dry air
30 relative_humidity hur 1 time, height relative to liquid
31 relative_humidity_over_ice huri 1 time, height relative to ice
32 air_potential_temperature theta K time, height
33 specific_turbulent_kinetic_energy_resolved tke_res m2 s-2 time, height resolved only; compute mean values horizontall...
34 specific_turbulent_kinetic_energy_sgs tke_sgs m2 s-2 time, height SGS only
35 mass_mixing_ratio_of_cloud_liquid_water_in_air qlc kg kg-1 time, height scene (all sky) per kg dry air; default breakd...
36 mass_mixing_ratio_of_rain_water_in_air qlr kg kg-1 time, height
37 mass_mixing_ratio_of_cloud_ice_in_air qic kg kg-1 time, height default breakdown to cloud_ice, snow and graup...
38 mass_mixing_ratio_of_snow_in_air qis kg kg-1 time, height
39 mass_mixing_ratio_of_graupel_in_air qig kg kg-1 time, height
40 number_of_liquid_cloud_droplets_in_air nlc kg-1 time, height scene (all sky) per kg dry air; default breakd...
41 number_of_rain_drops_in_air nlr kg-1 time, height
42 number_of_cloud_ice_crystals_in_air nic kg-1 time, height default breakdown to cloud_ice, snow and graup...
43 number_of_snow_crystals_in_air nis kg-1 time, height
44 number_of_graupel_crystals_in_air nig kg-1 time, height
45 number_of_total_aerosol_mode1 na1 kg-1 time, height when using prognostic aerosol; scene (all sky)...
46 number_of_total_aerosol_mode2 na2 kg-1 time, height accumulation mode
47 number_of_total_aerosol_mode3 na3 kg-1 time, height sea spray mode
48 number_of_liquid_cloud_droplets_in_cloud nlcic kg-1 time, height per kg dry air where cloud water > 0.01 g/kg; ...
49 number_of_ice_crystals_in_cloud niic kg-1 time, height total ice hydrometeors per kg dry air where cl...
50 dissipation_rate_of_turbulent_kinetic_energy eps m2 s-3 time, height total (including numerical diffusion contribut...
51 zonal_momentum_flux uw m2 s-2 time, height total (resolved plus SGS); for resolved compon...
52 meridional_momentum_flux vw m2 s-2 time, height total (resolved plus SGS); for resolved compon...
53 variance_of_upward_air_velocity w2 m2 s-2 time, height total (resolved plus SGS); for resolved compon...
54 vertical_flux_potential_temperature wth K m s-1 time, height total (resolved plus SGS); for resolved compon...
55 vertical_flux_liquid_ice_water_potential_tempe... vf_thli K m s-1 time, height total (resolved plus SGS); include sedimentati...
56 vertical_flux_water_vapor wqv kg kg-1 m s-1 time, height total (resolved plus SGS); for resolved compon...
57 vertical_flux_total_water vf_qt kg kg-1 m s-1 time, height total (resolved plus SGS); vapor+all liquid+al...
58 area_fraction_of_liquid_cloud flc 1 time, height fraction of cells with cloud water threshold o...
59 precipitation_flux_in_air prf kg m-2 s-1 time, height liquid and ice phase, all hydrometeors
60 precipitation_flux_in_air_in_ice_phase prfi kg m-2 s-1 time, height
61 downwelling_longwave_flux_in_air rld W m-2 time, height
62 upwelling_longwave_flux_in_air rlu W m-2 time, height
63 tendency_of_air_potential_temperature_due_to_r... dth_rad K s-1 time, height
64 tendency_of_air_potential_temperature_due_to_m... dth_micro K s-1 time, height including net condensation and precipitation
65 tendency_of_air_potential_temperature_due_to_m... dth_turb K s-1 time, height including surface fluxes
66 tendency_of_water_vapor_mixing_ratio_due_to_mi... dq_micro kg kg-1 s-1 time, height including net condensation and precipitation
67 tendency_of_water_vapor_mixing_ratio_due_to_mi... dq_turb kg kg-1 s-1 time, height including surface fluxes
68 tendency_of_aerosol_number_due_to_warm_microph... dna_micro_warm kg-1 s-1 time, height activated and unactivated aerosol (sum over al...
69 tendency_of_aerosol_number_due_to_cold_microph... dna_micro_cold kg-1 s-1 time, height activated and unactivated aerosol (sum over al...
70 tendency_of_aerosol_number_due_to_mixing dna_turb kg-1 s-1 time, height activated and unactivated aerosol (sum over al...
71 tendency_of_ice_number_due_to_heterogeneous_fr... dni_het kg-1 s-1 time, height sum of primary ice nucleation due to activatio...
72 tendency_of_ice_number_due_to_secondary_ice_pr... dni_sip kg-1 s-1 time, height sum of secondary ice formation processes (e.g....
73 tendency_of_ice_number_due_to_homogeneous_free... dni_hom kg-1 s-1 time, height ice nucleation source due to homogoeneous free...

Match DHARMA variables to requested outputs#

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

# drop comments
vars_mean_list = vars_mean_list#.drop(columns='comment (10-min average reported at endpoints, green=minimum)')

# 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)
# identify requested variables with only time dimension
vars_mean_scas = vars_mean_list[vars_mean_list['dimensions']=='time']

# match to DHARMA variable names and specify conversion factors
for index in vars_mean_scas.index:
    standard_name = vars_mean_list.standard_name.iat[index]
    if standard_name=='surface_pressure': 
        vars_mean_list.model_name.iat[index] = 'Psurf'
    if standard_name=='surface_temperature': 
        vars_mean_list.model_name.iat[index] = 'avg_T_sfc'
    if standard_name=='surface_friction_velocity': 
        vars_mean_list.model_name.iat[index] = 'avg_ustar'
    if standard_name=='surface_roughness_length_for_momentum_in_air':
        vars_mean_list.model_name.iat[index] = 'avg_z0'
    if standard_name=='surface_roughness_length_for_heat_in_air':
        vars_mean_list.model_name.iat[index] = 'avg_z0h'
    if standard_name=='surface_roughness_length_for_humidity_in_air':
        # same as roughness length for heat in DHARMA
        vars_mean_list.model_name.iat[index] = 'avg_z0h'
    if standard_name=='surface_upward_sensible_heat_flux': 
        vars_mean_list.model_name.iat[index] = 'avg_T_flx'
    if standard_name=='surface_upward_latent_heat_flux': 
        vars_mean_list.model_name.iat[index] = 'avg_qv_flx'
    if standard_name=='obukhov_length': 
        vars_mean_list.model_name.iat[index] = 'avg_obk'
    if standard_name=='atmosphere_mass_content_of_liquid_cloud_water': 
        vars_mean_list.model_name.iat[index] = 'cwp'
        vars_mean_list.conv_factor.iat[index] = 1/1000.
    if standard_name=='atmosphere_mass_content_of_rain_water': 
        vars_mean_list.model_name.iat[index] = 'rwp'
        vars_mean_list.conv_factor.iat[index] = 1/1000.
    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] = 1/1000.
    if standard_name=='cloud_area_fraction': 
        vars_mean_list.model_name.iat[index] = 'colf_opd'
    if standard_name=='optical_depth': 
        vars_mean_list.model_name.iat[index] = 'opd_tot'
    if standard_name=='optical_depth_of_liquid_cloud': 
        vars_mean_list.model_name.iat[index] = 'opd_cloud'
    if standard_name=='precipitation_flux_at_surface': 
        vars_mean_list.model_name.iat[index] = 'avg_precip'
        vars_mean_list.conv_factor.iat[index] = 1/3600.
    if do_ice:
        if standard_name=='precipitation_flux_at_surface_in_ice_phase': 
            vars_mean_list.model_name.iat[index] = 'avg_precip_ice'
            vars_mean_list.conv_factor.iat[index] = 1/3600.
    if standard_name=='optical_depth_of_cloud_droplets': 
        vars_mean_list.model_name.iat[index] = 'opd_cloud'
    if standard_name=='toa_outgoing_longwave_flux': 
        vars_mean_list.model_name.iat[index] = 'LWupTOA'
    if standard_name=='surface_downwelling_longwave_flux': 
        vars_mean_list.model_name.iat[index] = 'LWdnSFC'  
    if standard_name=='surface_upwelling_longwave_flux': 
        vars_mean_list.model_name.iat[index] = 'LWupSFC'
    if do_progNa:
        if standard_name=='surface_sea_spray_number_flux': 
            vars_mean_list.model_name.iat[index] = 'ssaf'
# identify requested variables with time and vertical dimensions
vars_mean_snds = vars_mean_list[vars_mean_list['dimensions']=='time, height']

# match to DHARMA variable names and specify conversion factors
for index in vars_mean_snds.index:
    standard_name = vars_mean_list.standard_name.iat[index]
    if standard_name=='air_pressure': 
        vars_mean_list.model_name.iat[index] = 'pressure'
    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] = 'rhobar'
    if standard_name=='air_temperature': 
        vars_mean_list.model_name.iat[index] = 'T'
    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] = 1/100.
    if standard_name=='relative_humidity_over_ice': 
        vars_mean_list.model_name.iat[index] = 'RHI'
        vars_mean_list.conv_factor.iat[index] = 1/100.
    if standard_name=='air_potential_temperature': 
        vars_mean_list.model_name.iat[index] = 'theta'
    if standard_name=='specific_turbulent_kinetic_energy_resolved': 
        vars_mean_list.model_name.iat[index] = 'tkeavg'
    if standard_name=='specific_turbulent_kinetic_energy_sgs': 
        vars_mean_list.model_name.iat[index] = 'tke_smag'
    if standard_name=='mass_mixing_ratio_of_cloud_liquid_water_in_air': 
        vars_mean_list.model_name.iat[index] = 'qc'
    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_cloud_ice_in_air': 
            vars_mean_list.model_name.iat[index] = 'qic'
        if standard_name=='mass_mixing_ratio_of_snow_in_air': 
            vars_mean_list.model_name.iat[index] = 'qif'
        if standard_name=='mass_mixing_ratio_of_graupel_in_air': 
            vars_mean_list.model_name.iat[index] = 'qid'
    if standard_name=='number_of_liquid_cloud_droplets_in_air': 
        vars_mean_list.model_name.iat[index] = 'nc'
    if standard_name=='number_of_rain_drops_in_air': 
        vars_mean_list.model_name.iat[index] = 'nr'
    if do_ice:
        if standard_name=='number_of_cloud_ice_crystals_in_air': 
            vars_mean_list.model_name.iat[index] = 'nic'
        if standard_name=='number_of_snow_crystals_in_air': 
            vars_mean_list.model_name.iat[index] = 'nif'
        if standard_name=='number_of_graupel_crystals_in_air': 
            vars_mean_list.model_name.iat[index] = 'nid'    
    if do_progNa:
        if standard_name=='number_of_total_aerosol_mode1': 
            vars_mean_list.model_name.iat[index] = 'na_1'
        if standard_name=='number_of_total_aerosol_mode2': 
            vars_mean_list.model_name.iat[index] = 'na_2'
        if standard_name=='number_of_total_aerosol_mode3': 
            vars_mean_list.model_name.iat[index] = 'na_3'
    if standard_name=='number_of_liquid_cloud_droplets_in_cloud': 
        vars_mean_list.model_name.iat[index] = 'nlcic'
    if do_ice:
        if standard_name=='number_of_ice_crystals_in_cloud': 
            vars_mean_list.model_name.iat[index] = 'niic'
    if standard_name=='dissipation_rate_of_turbulent_kinetic_energy': 
        vars_mean_list.model_name.iat[index] = 'eps'
    if standard_name=='zonal_momentum_flux': 
        vars_mean_list.model_name.iat[index] = 'uw_zt'
    if standard_name=='meridional_momentum_flux': 
        vars_mean_list.model_name.iat[index] = 'vw_zt'
    if standard_name=='variance_of_upward_air_velocity': 
        vars_mean_list.model_name.iat[index] = 'w2_zt'
    if standard_name=='vertical_flux_potential_temperature': 
        vars_mean_list.model_name.iat[index] = 'wth_zt'
    if standard_name=='vertical_flux_liquid_ice_water_potential_temperature': 
        vars_mean_list.model_name.iat[index] = 'wthli_zt'
    if standard_name=='vertical_flux_water_vapor': 
        vars_mean_list.model_name.iat[index] = 'wqv_zt'
    if standard_name=='vertical_flux_total_water': 
        vars_mean_list.model_name.iat[index] = 'wqt_zt'
    if standard_name=='area_fraction_of_liquid_cloud': 
        vars_mean_list.model_name.iat[index] = 'cloud_f'
    if standard_name=='precipitation_flux_in_air': 
        vars_mean_list.model_name.iat[index] = 'PF'
        vars_mean_list.conv_factor.iat[index] = 1/(3600.)
    if standard_name=='precipitation_flux_in_air_in_ice_phase': 
        vars_mean_list.model_name.iat[index] = 'PFi'
        vars_mean_list.conv_factor.iat[index] = 1/(3600.)
    if standard_name=='downwelling_longwave_flux_in_air': 
        vars_mean_list.model_name.iat[index] = 'LWdn'
    if standard_name=='upwelling_longwave_flux_in_air': 
        vars_mean_list.model_name.iat[index] = 'LWup'
    if standard_name=='tendency_of_air_potential_temperature_due_to_radiative_heating': 
        vars_mean_list.model_name.iat[index] = 'Srad'
        vars_mean_list.conv_factor.iat[index] = 1/3600.
    if standard_name=='tendency_of_air_potential_temperature_due_to_microphysics': 
        vars_mean_list.model_name.iat[index] = 'dth_micro'
        vars_mean_list.conv_factor.iat[index] = 1/3600.
    if standard_name=='tendency_of_air_potential_temperature_due_to_mixing': 
        vars_mean_list.model_name.iat[index] = 'dth_turb'
    if standard_name=='tendency_of_water_vapor_mixing_ratio_due_to_microphysics': 
        vars_mean_list.model_name.iat[index] = 'dq_micro'
        vars_mean_list.conv_factor.iat[index] = 1/3.6e6
    if standard_name=='tendency_of_water_vapor_mixing_ratio_due_to_mixing': 
        vars_mean_list.model_name.iat[index] = 'dq_turb'
    if do_progNa:
        if standard_name=='tendency_of_aerosol_number_due_to_warm_microphysics': 
            vars_mean_list.model_name.iat[index] = 'na_loss_liq'
            vars_mean_list.conv_factor.iat[index] = -1.
        if standard_name=='tendency_of_aerosol_number_due_to_mixing': 
            vars_mean_list.model_name.iat[index] = 'dna_mixing'
        if do_ice:
            if standard_name=='tendency_of_aerosol_number_due_to_cold_microphysics': 
                vars_mean_list.model_name.iat[index] = 'na_loss_ice'
                vars_mean_list.conv_factor.iat[index] = -1.
    if do_ice:
        if standard_name=='tendency_of_ice_number_due_to_heterogeneous_freezing': 
            vars_mean_list.model_name.iat[index] = 'Sice_het'
        if standard_name=='tendency_of_ice_number_due_to_secondary_ice_production': 
            vars_mean_list.model_name.iat[index] = 'Sice_sec'
        if standard_name=='tendency_of_ice_number_due_to_homogeneous_freezing': 
            vars_mean_list.model_name.iat[index] = 'Sice_hom'

vars_mean_list[3:] # echo variables (first rows are dimensions)
standard_name variable_id units dimensions comment (10-min average reported at endpoints, green=minimum) model_name conv_factor
3 surface_pressure ps Pa time Psurf 1.000000e+00
4 surface_temperature ts K time avg_T_sfc 1.000000e+00
5 surface_friction_velocity ustar m s-1 time avg_ustar 1.000000e+00
6 surface_roughness_length_for_momentum_in_air z0 m time avg_z0 1.000000e+00
7 surface_roughness_length_for_heat_in_air z0h m time avg_z0h 1.000000e+00
8 surface_roughness_length_for_humidity_in_air z0q m time avg_z0h 1.000000e+00
9 surface_upward_sensible_heat_flux hfss W m-2 time avg_T_flx 1.000000e+00
10 surface_upward_latent_heat_flux hfls W m-2 time avg_qv_flx 1.000000e+00
11 obukhov_length ol m time avg_obk 1.000000e+00
12 atmosphere_mass_content_of_liquid_cloud_water lwpc kg m-2 time default breakdown to cloud water and rain in a... cwp 1.000000e-03
13 atmosphere_mass_content_of_rain_water lwpr kg m-2 time rwp 1.000000e-03
14 atmosphere_mass_content_of_ice_water iwp kg m-2 time all ice-phase hydrometeors iwp 1.000000e-03
15 cloud_area_fraction clt 1 time fraction of atmospheric columns with total hyd... colf_opd 1.000000e+00
16 optical_depth od 1 time scene (all sky); mid-visible, all hydrometeors... opd_tot 1.000000e+00
17 optical_depth_of_liquid_cloud odlc 1 time as for optical_depth but cloud droplets only opd_cloud 1.000000e+00
18 precipitation_flux_at_surface pr kg m-2 s-1 time liquid and ice phase, all hydrometeors avg_precip 2.777778e-04
19 precipitation_flux_at_surface_in_ice_phase pri kg m-2 s-1 time ice phase hydrometeors only avg_precip_ice 2.777778e-04
20 toa_outgoing_longwave_flux rlut W m-2 time at top of atmosphere LWupTOA 1.000000e+00
21 surface_downwelling_longwave_flux rlds W m-2 time LWdnSFC 1.000000e+00
22 surface_upwelling_longwave_flux rlus W m-2 time LWupSFC 1.000000e+00
23 surface_sea_spray_number_flux ssaf m-2 s-1 time when using prognostic aerosol; emission only (... missing data 1.000000e+00
24 air_pressure pa Pa time, height pressure 1.000000e+00
25 eastward_wind ua m s-1 time, height u 1.000000e+00
26 northward_wind va m s-1 time, height v 1.000000e+00
27 air_dry_density rhoa kg m-3 time, height rhobar 1.000000e+00
28 air_temperature ta K time, height T 1.000000e+00
29 water_vapor_mixing_ratio qv kg kg-1 time, height per kg dry air qv 1.000000e+00
30 relative_humidity hur 1 time, height relative to liquid RH 1.000000e-02
31 relative_humidity_over_ice huri 1 time, height relative to ice RHI 1.000000e-02
32 air_potential_temperature theta K time, height theta 1.000000e+00
33 specific_turbulent_kinetic_energy_resolved tke_res m2 s-2 time, height resolved only; compute mean values horizontall... tkeavg 1.000000e+00
34 specific_turbulent_kinetic_energy_sgs tke_sgs m2 s-2 time, height SGS only tke_smag 1.000000e+00
35 mass_mixing_ratio_of_cloud_liquid_water_in_air qlc kg kg-1 time, height scene (all sky) per kg dry air; default breakd... qc 1.000000e+00
36 mass_mixing_ratio_of_rain_water_in_air qlr kg kg-1 time, height qr 1.000000e+00
37 mass_mixing_ratio_of_cloud_ice_in_air qic kg kg-1 time, height default breakdown to cloud_ice, snow and graup... qic 1.000000e+00
38 mass_mixing_ratio_of_snow_in_air qis kg kg-1 time, height qif 1.000000e+00
39 mass_mixing_ratio_of_graupel_in_air qig kg kg-1 time, height qid 1.000000e+00
40 number_of_liquid_cloud_droplets_in_air nlc kg-1 time, height scene (all sky) per kg dry air; default breakd... nc 1.000000e+00
41 number_of_rain_drops_in_air nlr kg-1 time, height nr 1.000000e+00
42 number_of_cloud_ice_crystals_in_air nic kg-1 time, height default breakdown to cloud_ice, snow and graup... nic 1.000000e+00
43 number_of_snow_crystals_in_air nis kg-1 time, height nif 1.000000e+00
44 number_of_graupel_crystals_in_air nig kg-1 time, height nid 1.000000e+00
45 number_of_total_aerosol_mode1 na1 kg-1 time, height when using prognostic aerosol; scene (all sky)... missing data 1.000000e+00
46 number_of_total_aerosol_mode2 na2 kg-1 time, height accumulation mode missing data 1.000000e+00
47 number_of_total_aerosol_mode3 na3 kg-1 time, height sea spray mode missing data 1.000000e+00
48 number_of_liquid_cloud_droplets_in_cloud nlcic kg-1 time, height per kg dry air where cloud water > 0.01 g/kg; ... nlcic 1.000000e+00
49 number_of_ice_crystals_in_cloud niic kg-1 time, height total ice hydrometeors per kg dry air where cl... niic 1.000000e+00
50 dissipation_rate_of_turbulent_kinetic_energy eps m2 s-3 time, height total (including numerical diffusion contribut... eps 1.000000e+00
51 zonal_momentum_flux uw m2 s-2 time, height total (resolved plus SGS); for resolved compon... uw_zt 1.000000e+00
52 meridional_momentum_flux vw m2 s-2 time, height total (resolved plus SGS); for resolved compon... vw_zt 1.000000e+00
53 variance_of_upward_air_velocity w2 m2 s-2 time, height total (resolved plus SGS); for resolved compon... w2_zt 1.000000e+00
54 vertical_flux_potential_temperature wth K m s-1 time, height total (resolved plus SGS); for resolved compon... wth_zt 1.000000e+00
55 vertical_flux_liquid_ice_water_potential_tempe... vf_thli K m s-1 time, height total (resolved plus SGS); include sedimentati... wthli_zt 1.000000e+00
56 vertical_flux_water_vapor wqv kg kg-1 m s-1 time, height total (resolved plus SGS); for resolved compon... wqv_zt 1.000000e+00
57 vertical_flux_total_water vf_qt kg kg-1 m s-1 time, height total (resolved plus SGS); vapor+all liquid+al... wqt_zt 1.000000e+00
58 area_fraction_of_liquid_cloud flc 1 time, height fraction of cells with cloud water threshold o... cloud_f 1.000000e+00
59 precipitation_flux_in_air prf kg m-2 s-1 time, height liquid and ice phase, all hydrometeors PF 2.777778e-04
60 precipitation_flux_in_air_in_ice_phase prfi kg m-2 s-1 time, height PFi 2.777778e-04
61 downwelling_longwave_flux_in_air rld W m-2 time, height LWdn 1.000000e+00
62 upwelling_longwave_flux_in_air rlu W m-2 time, height LWup 1.000000e+00
63 tendency_of_air_potential_temperature_due_to_r... dth_rad K s-1 time, height Srad 2.777778e-04
64 tendency_of_air_potential_temperature_due_to_m... dth_micro K s-1 time, height including net condensation and precipitation dth_micro 2.777778e-04
65 tendency_of_air_potential_temperature_due_to_m... dth_turb K s-1 time, height including surface fluxes dth_turb 1.000000e+00
66 tendency_of_water_vapor_mixing_ratio_due_to_mi... dq_micro kg kg-1 s-1 time, height including net condensation and precipitation dq_micro 2.777778e-07
67 tendency_of_water_vapor_mixing_ratio_due_to_mi... dq_turb kg kg-1 s-1 time, height including surface fluxes dq_turb 1.000000e+00
68 tendency_of_aerosol_number_due_to_warm_microph... dna_micro_warm kg-1 s-1 time, height activated and unactivated aerosol (sum over al... missing data 1.000000e+00
69 tendency_of_aerosol_number_due_to_cold_microph... dna_micro_cold kg-1 s-1 time, height activated and unactivated aerosol (sum over al... missing data 1.000000e+00
70 tendency_of_aerosol_number_due_to_mixing dna_turb kg-1 s-1 time, height activated and unactivated aerosol (sum over al... missing data 1.000000e+00
71 tendency_of_ice_number_due_to_heterogeneous_fr... dni_het kg-1 s-1 time, height sum of primary ice nucleation due to activatio... Sice_het 1.000000e+00
72 tendency_of_ice_number_due_to_secondary_ice_pr... dni_sip kg-1 s-1 time, height sum of secondary ice formation processes (e.g.... Sice_sec 1.000000e+00
73 tendency_of_ice_number_due_to_homogeneous_free... dni_hom kg-1 s-1 time, height ice nucleation source due to homogoeneous free... Sice_hom 1.000000e+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”.

# create DEPHY output file
dephy_filename = './' + my_gitdir + my_outfile
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')

# create global attributes
dephy_file.title='DHARMA LES results for COMBLE-MIP case: fixed Nd and Ni'
dephy_file.reference='https://github.com/ARM-Development/comble-mip'
dephy_file.authors='Ann Fridlind (ann.fridlind@nasa.gov) and Florian Tornow (florian.tornow@nasa.gov)'
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_DHARMA_LES_output_to_dephy_format.ipynb'
dephy_file.startDate=start_dtime
dephy_file.force_geo=1
dephy_file.surfaceType='ocean (after spin-up)'
dephy_file.surfaceForcing='ts (after spin-up)'
dephy_file.lat=str(dharma_params['Coriolis'].lat) + ' deg N'
dephy_file.dx=str(dharma_params['geometry'].L_x/dharma_params['geometry'].nx) + ' m'
dephy_file.dy=str(dharma_params['geometry'].L_y/dharma_params['geometry'].ny) + ' m'
dephy_file.dz='see zf variable'
dephy_file.nx=dharma_params['geometry'].nx
dephy_file.ny=dharma_params['geometry'].ny
dephy_file.nz=dharma_params['geometry'].nz

# create dimensions
nz = dharma_snds.sizes['zt']
zf = dephy_file.createDimension('zf', nz)
zf = dephy_file.createVariable('zf', np.float64, ('zf',))
zf.units = 'm'
zf.long_name = 'height'
zf[:] = dharma_snds['zt'].data

ze = dephy_file.createDimension('ze', nz)
ze = dephy_file.createVariable('ze', np.float64, ('ze',))
ze.units = 'm'
ze.long_name = 'layer_top_height'
ze[:] = dharma_snds['zw'].data[1:]

nt = dharma_snds.sizes['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
delta_t = (dharma_snds['time'].data[1]-dharma_snds['time'].data[0])/np.timedelta64(1, "s")
time[:] = np.arange(nt)*delta_t

# 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
        if vars_mean_list.model_name.iat[index]!='missing data' and mod_name in dharma_scas:
            if (var_name == 'z0'):
                new_sca[:] = np.round(dharma_scas[mod_name].data*c_factor,4)
            elif (var_name == 'z0h')| (var_name == 'z0q'):
                new_sca[:] = np.round(dharma_scas[mod_name].data*c_factor,7)
            else:
                new_sca[:] = dharma_scas[mod_name].data*c_factor
    if vars_mean_list.dimensions.iat[index]=='time, height':
        new_snd = dephy_file.createVariable(var_name, np.float64, ('time','zf'))
        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' and mod_name in dharma_snds: 
            new_snd[:] = dharma_snds[mod_name].data*c_factor

print(dephy_file)
dephy_file.close()
<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF3_CLASSIC data model, file format NETCDF3):
    title: DHARMA LES results for COMBLE-MIP case: fixed Nd and Ni
    reference: https://github.com/ARM-Development/comble-mip
    authors: Ann Fridlind (ann.fridlind@nasa.gov) and Florian Tornow (florian.tornow@nasa.gov)
    source: /data/home/floriantornow/dharma/sandbox/case0313_diag_ice25_nomiz_dx100_specZ0_hm_largedomain/6-20h/dharma.scalars.cdf
    version: 2024-12-19 18:43:26
    format_version: DEPHY SCM format version 1.6
    script: convert_DHARMA_LES_output_to_dephy_format.ipynb
    startDate: 2020-03-12 22:00:00.0
    force_geo: 1
    surfaceType: ocean (after spin-up)
    surfaceForcing: ts (after spin-up)
    lat: 73.0 deg N
    dx: 100.0 m
    dy: 100.0 m
    dz: see zf variable
    nx: 1280
    ny: 1280
    nz: 159
    dimensions(sizes): zf(159), ze(159), time(121)
    variables(dimensions): float64 zf(zf), float64 ze(ze), float64 time(time), 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 lwpc(time), float64 lwpr(time), float64 iwp(time), float64 clt(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 pa(time, zf), float64 ua(time, zf), float64 va(time, zf), float64 rhoa(time, zf), float64 ta(time, zf), float64 qv(time, zf), float64 hur(time, zf), float64 huri(time, zf), float64 theta(time, zf), float64 tke_res(time, zf), float64 tke_sgs(time, zf), float64 qlc(time, zf), float64 qlr(time, zf), float64 qic(time, zf), float64 qis(time, zf), float64 qig(time, zf), float64 nlc(time, zf), float64 nlr(time, zf), float64 nic(time, zf), float64 nis(time, zf), float64 nig(time, zf), float64 na1(time, zf), float64 na2(time, zf), float64 na3(time, zf), float64 nlcic(time, zf), float64 niic(time, zf), float64 eps(time, zf), float64 uw(time, zf), float64 vw(time, zf), float64 w2(time, zf), float64 wth(time, zf), float64 vf_thli(time, zf), float64 wqv(time, zf), float64 vf_qt(time, zf), float64 flc(time, zf), float64 prf(time, zf), float64 prfi(time, zf), float64 rld(time, zf), float64 rlu(time, zf), float64 dth_rad(time, zf), float64 dth_micro(time, zf), float64 dth_turb(time, zf), float64 dq_micro(time, zf), float64 dq_turb(time, zf), float64 dna_micro_warm(time, zf), float64 dna_micro_cold(time, zf), float64 dna_turb(time, zf), float64 dni_het(time, zf), float64 dni_sip(time, zf), float64 dni_hom(time, zf)
    groups: 

Check output file#

dephy_check = xr.open_dataset(dephy_filename)
dephy_check
<xarray.Dataset> Size: 8MB
Dimensions:         (zf: 159, ze: 159, time: 121)
Coordinates:
  * zf              (zf) float64 1kB 10.0 32.5 60.0 ... 6.92e+03 6.975e+03
  * ze              (ze) float64 1kB 20.0 45.0 75.0 ... 6.89e+03 6.95e+03 7e+03
  * time            (time) datetime64[ns] 968B 2020-03-12T22:00:00 ... 2020-0...
Data variables: (12/71)
    ps              (time) float64 968B ...
    ts              (time) float64 968B ...
    ustar           (time) float64 968B ...
    z0              (time) float64 968B ...
    z0h             (time) float64 968B ...
    z0q             (time) float64 968B ...
    ...              ...
    dna_micro_warm  (time, zf) float64 154kB ...
    dna_micro_cold  (time, zf) float64 154kB ...
    dna_turb        (time, zf) float64 154kB ...
    dni_het         (time, zf) float64 154kB ...
    dni_sip         (time, zf) float64 154kB ...
    dni_hom         (time, zf) float64 154kB ...
Attributes: (12/18)
    title:           DHARMA LES results for COMBLE-MIP case: fixed Nd and Ni
    reference:       https://github.com/ARM-Development/comble-mip
    authors:         Ann Fridlind (ann.fridlind@nasa.gov) and Florian Tornow ...
    source:          /data/home/floriantornow/dharma/sandbox/case0313_diag_ic...
    version:         2024-12-19 18:43:26
    format_version:  DEPHY SCM format version 1.6
    ...              ...
    dx:              100.0 m
    dy:              100.0 m
    dz:              see zf variable
    nx:              1280
    ny:              1280
    nz:              159
dephy_check['z0h'].values
array([5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06,
       5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06,
       5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06,
       5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06,
       5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06,
       5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06,
       5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06,
       5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06,
       5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06,
       5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06,
       5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06,
       5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06,
       5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06,
       5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06,
       5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06,
       5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06,
       5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06, 5.5e-06,
       5.5e-06, 5.5e-06])
import math
math.isnan(float('nan'))
float('NaN')
nan