Example: convert DHARMA LES 3D output to DEPHY format

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

Contributed by Florian Tornow Columbia University & NASA/GISS and 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
import re
import sys

import pathlib

%run ../plotting/functions_plotting.py
%run function_extract_2d.py 

R_const = 287 #J kg^-1 K^-1

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)

# ProgNa with ice test
my_readdir = 'case0313_diag_ice25_nomiz_dx100_specZ0_sbhomo'
my_outfile = 'DHARMA_Lx25_dx100_FixN_3D.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-20h']

Read DHARMA data

Locate 2D files

## locate all files that have the expected filename (*.tgz)
path = my_rundir + my_outdirs[0]
direc = pathlib.Path(path)
NCFILES = list(direc.rglob("*tgz"))
NCFILES_STR = [str(p) for p in pathlib.Path(path).rglob('*.tgz')]

NCFILES_STR

input_filenames = []
count = 0
for fn in NCFILES: 
    if ('PLT' in NCFILES_STR[count]):
        input_filenames.append(fn)
    count += 1

input_filenames

[PosixPath('/data/home/floriantornow/dharma/sandbox/case0313_diag_ice25_nomiz_dx100_specZ0_sbhomo/0-20h/dharma_PLT_029100.tgz')]

Read 2D files

## obtain model horizontal resolution
cdf_filename = my_rundir + my_outdirs[0] + '/dharma.cdf'
dharma_params = xr.open_dataset(cdf_filename)
print(cdf_filename)

nx = dharma_params['geometry'].nx
Lx = dharma_params['geometry'].L_x

## obtain vertical resolution
nz = dharma_params['geometry'].nz

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

/data/home/floriantornow/dharma/sandbox/case0313_diag_ice25_nomiz_dx100_specZ0_sbhomo/0-20h/dharma.cdf

## obtain auxiliary info from sounding file
sound_filename = my_rundir + my_outdirs[0] + '/dharma.soundings.cdf'
dharma_sound = xr.open_dataset(sound_filename)
theta_00 = dharma_sound.theta_00

xy_res = Lx/nx
xy_res

100.0

## read profiles to complete pressure profile
var_vec_1d = ['lwpc'] # variables with ERA5 (longer time axis)
var_vec_2d = ['rhoa','ta','theta']
df_col_1d,df_col_2d = load_sims('../../output_les/',var_vec_1d,var_vec_2d,t_shift=-2,keyword='DHARMA_Lx25_dx100_FixN.')

init_prof = df_col_2d[abs(round(df_col_2d.time)-df_col_2d['time'].min())<2.]

Loading variables: f(time)
../../output_les/dharma/sandbox/DHARMA_Lx25_dx100_FixN.nc
Loading variables: f(time,height)
../../output_les/dharma/sandbox/DHARMA_Lx25_dx100_FixN.nc

## prepare pressure field translation
rhobar = dharma_sound['rhobar'].isel(time=0).values
th_snd = theta_00*(1. + dharma_sound['th'].isel(time=0).values)
T_snd = dharma_sound['T'].isel(time=0).values

pii = T_snd/th_snd
penv_mks = 1e6*pii**(1/0.286)*0.1

time_vec = []

counter = 0
for input_filename in input_filenames[0:1]:
        
    ## also check if partner file ("alt_3d") is available
    input_filename_2 = Path(str(input_filename).replace('PLT', 'alt_3d'))
    
    ## obtain timestamp
    time_vec.append(obtain_time(path,input_filename,type='PLT_'))

    ## ascertain horizontal resolution and prepare horizontal grid
    xy_res = Lx/nx
    x_vec = np.arange(-Lx/2, Lx/2, xy_res)
    y_vec = np.arange(-Lx/2, Lx/2, xy_res)
    z_vec = dharma_sound['zt'].values
        
    ## unpack tgz and extract relevant fields     
    pa_3d  = threed_loader(input_filename,vbase = 'press')
    pa_1d  = (pa_3d.mean(axis=2).mean(axis=1)*rhobar + penv_mks) # Pa
    pa_3d_act = np.tile(penv_mks.reshape(len(penv_mks),1,1),[nx,nx]) + np.tile(rhobar.reshape(len(rhobar),1,1),[nx,nx])*pa_3d # Pa
    
    ua_3d  = threed_loader(input_filename,vbase = 'u') #+ dharma_params['translate'].u
    va_3d  = threed_loader(input_filename,vbase = 'v') #+ dharma_params['translate'].v
    wa_3d  = threed_loader(input_filename,vbase = 'w')
    th_3d  = theta_00*(1. + threed_loader(input_filename,vbase = 'th'))
    ta_3d  = np.tile(pii.reshape(len(pii),1,1),[nx,nx]) * th_3d 
    rhoa_3d = pa_3d_act/ta_3d/R_const
    
    qv_3d   = threed_loader(input_filename,vbase = 'qv')
    qlc_3d  = threed_loader(input_filename,vbase = 'qc')
    qlr_3d  = threed_loader(input_filename,vbase = 'qr')
    qic_3d  = threed_loader(input_filename,vbase = 'qic')
    qis_3d  = threed_loader(input_filename,vbase = 'qif')
    qig_3d  = threed_loader(input_filename,vbase = 'qid')
    nlc_3d  = threed_loader(input_filename,vbase = 'nc')
    nlr_3d  = threed_loader(input_filename,vbase = 'nr')
    nic_3d  = threed_loader(input_filename,vbase = 'nic')
    nis_3d  = threed_loader(input_filename,vbase = 'nif')
    nig_3d  = threed_loader(input_filename,vbase = 'nid')
    eps_3d     = threed_loader(input_filename_2,vbase = 'kmd2')
    tend_cond_3d = threed_loader(input_filename_2,vbase = 'net_cond')
    tend_dep_3d  = threed_loader(input_filename_2,vbase = 'net_dep')
    tend_rim_3d  = threed_loader(input_filename_2,vbase = 'rime_prod')

    ## accumulate
    if counter == 0:
        pa_col = pa_1d
        ta_col = ta_3d
        rhoa_col = rhoa_3d
        ua_col = ua_3d
        va_col = va_3d
        wa_col = wa_3d
        qv_col = qv_3d
        qlc_col = qlc_3d
        qlr_col = qlr_3d
        qic_col = qic_3d
        qis_col = qis_3d
        qig_col = qig_3d
        nlc_col = nlc_3d
        nlr_col = nlr_3d
        nic_col = nic_3d
        nis_col = nis_3d
        nig_col = nig_3d
        eps_col = eps_3d
        tend_cond_col = tend_cond_3d
        tend_dep_col  = tend_dep_3d
        tend_rim_col  = tend_rim_3d
    else:
        pa_col = np.stack([pa_col,pa_1d],axis=0)
        ta_col = np.stack([ta_col,ta_3d],axis=0)
        rhoa_col = np.stack([rhoa_col,rhoa_3d],axis=0)
        ua_col = np.stack([ua_col,ua_3d],axis=0)
        va_col = np.stack([va_col,va_3d],axis=0)
        wa_col = np.stack([wa_col,wa_3d],axis=0)
        qv_col = np.stack([qv_col,qv_3d],axis=0)
        qlc_col = np.stack([qlc_col,qlc_3d],axis=0)
        qlr_col = np.stack([qlr_col,qlr_3d],axis=0)
        qic_col = np.stack([qic_col,qic_3d],axis=0)
        qis_col = np.stack([qis_col,qis_3d],axis=0)
        qig_col = np.stack([qig_col,qig_3d],axis=0)
        nlc_col = np.stack([nlc_col,nlc_3d],axis=0)
        nlr_col = np.stack([nlr_col,nlr_3d],axis=0)
        nic_col = np.stack([nic_col,nic_3d],axis=0)
        nis_col = np.stack([nis_col,nis_3d],axis=0)
        nig_col = np.stack([nig_col,nig_3d],axis=0)
        eps_col = np.stack([eps_col,eps_3d],axis=0)
        tend_cond_col = np.stack([tend_cond_col,tend_cond_3d],axis=0)
        tend_dep_col  = np.stack([tend_dep_col,tend_dep_3d],axis=0)
        tend_rim_col  = np.stack([tend_rim_col,tend_rim_3d],axis=0)
    counter += 1
 

|  29100 |   72000.0 |      0.8 | 0.4401E-07 | 0.4925E-08 | 0.0455 | 0.5321 |

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=1233994833&format=xlsx',
                              sheet_name='3D')

pd.set_option('display.max_rows', None)
vars_mean_list

	standard_name	variable_id	units	dimensions	comment (60-min instantaneous, red=required for EMC2)
0	time	time	s	–	dimension, seconds since 2020-03-12 18:00:00
1	x	x	m	–	dimension, grid cell center
2	y	y	m	–	dimension, grid cell center
3	height	zf	m	–	dimension, altitude of mid-level points above ...
4	air_pressure	pa	Pa	time, height	–
5	eastward_wind	ua	m s-1	time, height, x, y	–
6	northward_wind	va	m s-1	time, height, x, y	–
7	upward_air_velocity	wa	m s-1	time, height, x, y	–
8	air_temperature	ta	K	time, height, x, y	–
9	air_dry_density	rhoa	kg m-3	time, height, x, y	–
10	water_vapor_mixing_ratio	qv	kg kg-1	time, height, x, y	per kg dry air
11	mass_mixing_ratio_of_cloud_liquid_water_in_air	qlc	kg kg-1	time, height, x, y	per kg dry air; default breakdown to cloud wat...
12	mass_mixing_ratio_of_rain_water_in_air	qlr	kg kg-1	time, height, x, y	–
13	mass_mixing_ratio_of_cloud_ice_in_air	qic	kg kg-1	time, height, x, y	default breakdown to cloud_ice, snow and graup...
14	mass_mixing_ratio_of_snow_in_air	qis	kg kg-1	time, height, x, y	–
15	mass_mixing_ratio_of_graupel_in_air	qig	kg kg-1	time, height, x, y	–
16	number_of_dry_aerosol_in_air	na	kg-1	time, height, x, y	per kg dry air; when using prognostic aerosol
17	number_of_cloud_droplets_in_air	nlc	kg-1	time, height, x, y	when using prognostic aerosol
18	number_of_rain_droplets_in_air	nlr	kg-1	time, height, x, y	–
19	number_of_cloud_ice_crystals_in_air	nic	kg-1	time, height, x, y	default breakdown to cloud_ice, snow and graup...
20	number_of_snow_crystals_in_air	nis	kg-1	time, height, x, y	–
21	number_of_graupel_crystals_in_air	nig	kg-1	time, height, x, y	–
22	dissipation_rate_of_turbulent_kinetic_energy	eps	m2 s-3	time, height, x, y	including numerical diffusion contributions if...
23	mass_weighted_fall_speed_of_liquid_cloud_water	vmlc	m s-1	time, height, x, y	EMC2 uses mass-weighted fall-speed profiles fo...
24	mass_weighted_fall_speed_of_rain	vmlr	m s-1	time, height, x, y	–
25	mass_weighted_fall_speed_of_ice_cloud	vmic	m s-1	time, height, x, y	if other ice-phase categories are used, provid...
26	mass_weighted_fall_speed_of_snow	vmis	m s-1	time, height, x, y	–
27	mass_weighted_fall_speed_of_graupel	vmig	m s-1	time, height, x, y	–
28	tendency_of_condensation	tend_cond	s-1	time, height, x, y	combined tendencies of cloud liquid and rain w...
29	tendency_of_deposition	tend_dep	s-1	time, height, x, y	tendency of ice mixing ratio due to deposition...
30	tendency_of_riming	tend_rim	s-1	time, height, x, y	tendency of ice mixing ratio due to riming; su...

# 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='Florian Tornow (florian.tornow@nasa.gov) and Ann Fridlind (ann.fridlind@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-derivative'
dephy_file.script='convert_DHARMA_LES_output_to_dephy_format_3D.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
xvec = dephy_file.createDimension('x', nx)
xvec = dephy_file.createVariable('x', np.float64, ('x',))
xvec.units = 'm'
xvec.long_name = 'x'
xvec[:] = x_vec

yvec = dephy_file.createDimension('y', nx)
yvec = dephy_file.createVariable('y', np.float64, ('y',))
yvec.units = 'm'
yvec.long_name = 'y'
yvec[:] = y_vec

zvec = dephy_file.createDimension('height', nz)
zvec = dephy_file.createVariable('height', np.float64, ('height',))
zvec.units = 'm'
zvec.long_name = 'y'
zvec[:] = z_vec

nt = len(time_vec)#len(input_filenames)
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'
time[:] = time_vec[:]

# create and fill variables
for index in vars_mean_list.index[4:]:
    std_name = vars_mean_list.standard_name.iat[index]
    var_name = vars_mean_list['variable_id'].iat[index]
    if var_name == 'ta':
        mod_field = ta_col
    if var_name == 'pa':
        mod_field = pa_col
    if var_name == 'rhoa':
        mod_field = rhoa_col
    if var_name == 'ua':
        mod_field = ua_col
    if var_name == 'va':
        mod_field = va_col
    if var_name == 'wa':
        mod_field = wa_col
    if var_name == 'qv':
        mod_field = qv_col
    if var_name == 'qlc':
        mod_field = qlc_col
    if var_name == 'qlr':
        mod_field = qlr_col
    if var_name == 'qic':
        mod_field = qic_col
    if var_name == 'qis':
        mod_field = qis_col
    if var_name == 'qig':
        mod_field = qig_col
    if var_name == 'nlc':
        mod_field = nlc_col
    if var_name == 'nlr':
        mod_field = nlr_col
    if var_name == 'nic':
        mod_field = nic_col
    if var_name == 'nis':
        mod_field = nis_col
    if var_name == 'nig':
        mod_field = nig_col
    if var_name == 'eps':
        mod_field = eps_col
    if var_name == 'tend_cond':
        mod_field = tend_cond_col
    if var_name == 'tend_dep':
        mod_field = tend_dep_col
    if var_name == 'tend_rim':
        mod_field = tend_rim_col
    if (var_name in ['na','vmlc','vmlr','vmic','vmis','vmig']):
        continue
    
    if vars_mean_list.dimensions.iat[index]=='time, height':
        new_snd = dephy_file.createVariable(var_name, np.float64, ('time','height'))
        new_snd.units = vars_mean_list.units.iat[index]
        new_snd.long_name = std_name
        new_snd[:] = mod_field
        
    if vars_mean_list.dimensions.iat[index]=='time, height, x, y':
        new_snd = dephy_file.createVariable(var_name, np.float64, ('time','height','x','y'))
        new_snd.units = vars_mean_list.units.iat[index]
        new_snd.long_name = std_name
        new_snd[:] = mod_field

print(dephy_file)
dephy_file.close()

The file ./../../output_les/dharma/sandbox/DHARMA_Lx25_dx100_FixN_3D.nc has been deleted successfully
<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: Florian Tornow (florian.tornow@nasa.gov) and Ann Fridlind (ann.fridlind@nasa.gov)
    version: 2024-09-25 16:28:30
    format_version: DEPHY-derivative
    script: convert_DHARMA_LES_output_to_dephy_format_3D.ipynb
    startDate: 2020-03-12 22:00:00.0
    force_geo: 1
    surfaceType: ocean (after spin-up)
    surfaceForcing: ts (after spin-up)
    dimensions(sizes): x(256), y(256), height(159), time(1)
    variables(dimensions): float64 x(x), float64 y(y), float64 height(height), float64 time(time), float64 pa(time, height), float64 ua(time, height, x, y), float64 va(time, height, x, y), float64 wa(time, height, x, y), float64 ta(time, height, x, y), float64 rhoa(time, height, x, y), float64 qv(time, height, x, y), float64 qlc(time, height, x, y), float64 qlr(time, height, x, y), float64 qic(time, height, x, y), float64 qis(time, height, x, y), float64 qig(time, height, x, y), float64 nlc(time, height, x, y), float64 nlr(time, height, x, y), float64 nic(time, height, x, y), float64 nis(time, height, x, y), float64 nig(time, height, x, y), float64 eps(time, height, x, y), float64 tend_cond(time, height, x, y), float64 tend_dep(time, height, x, y), float64 tend_rim(time, height, x, y)
    groups: 

dephy_check = xr.open_dataset(dephy_filename)
dephy_check

<xarray.Dataset> Size: 2GB
Dimensions:    (x: 256, y: 256, height: 159, time: 1)
Coordinates:
  * x          (x) float64 2kB -1.28e+04 -1.27e+04 ... 1.26e+04 1.27e+04
  * y          (y) float64 2kB -1.28e+04 -1.27e+04 ... 1.26e+04 1.27e+04
  * height     (height) float64 1kB 10.0 32.5 60.0 ... 6.92e+03 6.975e+03
  * time       (time) datetime64[ns] 8B 2020-03-13T18:00:00
Data variables: (12/21)
    pa         (time, height) float64 1kB ...
    ua         (time, height, x, y) float64 83MB ...
    va         (time, height, x, y) float64 83MB ...
    wa         (time, height, x, y) float64 83MB ...
    ta         (time, height, x, y) float64 83MB ...
    rhoa       (time, height, x, y) float64 83MB ...
    ...         ...
    nis        (time, height, x, y) float64 83MB ...
    nig        (time, height, x, y) float64 83MB ...
    eps        (time, height, x, y) float64 83MB ...
    tend_cond  (time, height, x, y) float64 83MB ...
    tend_dep   (time, height, x, y) float64 83MB ...
    tend_rim   (time, height, x, y) float64 83MB ...
Attributes:
    title:           DHARMA LES results for COMBLE-MIP case: fixed Nd and Ni
    reference:       https://github.com/ARM-Development/comble-mip
    authors:         Florian Tornow (florian.tornow@nasa.gov) and Ann Fridlin...
    version:         2024-09-25 16:28:30
    format_version:  DEPHY-derivative
    script:          convert_DHARMA_LES_output_to_dephy_format_3D.ipynb
    startDate:       2020-03-12 22:00:00.0
    force_geo:       1
    surfaceType:     ocean (after spin-up)
    surfaceForcing:  ts (after spin-up)

xarray.Dataset

• Dimensions:
  • x: 256
  • y: 256
  • height: 159
  • time: 1
• Coordinates: (4)
  • x
    (x)
    float64
    -1.28e+04 -1.27e+04 ... 1.27e+04

    units :

    m

    long_name :

    x

    array([-12800., -12700., -12600., ...,  12500.,  12600.,  12700.])

  • y
    (y)
    float64
    -1.28e+04 -1.27e+04 ... 1.27e+04

    units :

    m

    long_name :

    y

    array([-12800., -12700., -12600., ...,  12500.,  12600.,  12700.])

  • height
    (height)
    float64
    10.0 32.5 ... 6.92e+03 6.975e+03

    units :

    m

    long_name :

    y

    array([  10. ,   32.5,   60. ,   92.5,  130. ,  170. ,  210. ,  250. ,  290. ,
            330. ,  370. ,  410. ,  450. ,  490. ,  530. ,  570. ,  610. ,  650. ,
            690. ,  730. ,  770. ,  810. ,  850. ,  890. ,  930. ,  970. , 1010. ,
           1050. , 1090. , 1130. , 1170. , 1210. , 1250. , 1290. , 1330. , 1370. ,
           1410. , 1450. , 1490. , 1530. , 1570. , 1610. , 1650. , 1690. , 1730. ,
           1770. , 1810. , 1850. , 1890. , 1930. , 1970. , 2010. , 2050. , 2090. ,
           2130. , 2170. , 2210. , 2250. , 2290. , 2330. , 2370. , 2410. , 2450. ,
           2490. , 2530. , 2570. , 2610. , 2650. , 2690. , 2730. , 2770. , 2810. ,
           2850. , 2890. , 2930. , 2970. , 3010. , 3050. , 3090. , 3130. , 3170. ,
           3210. , 3250. , 3290. , 3330. , 3370. , 3410. , 3450. , 3490. , 3530. ,
           3570. , 3610. , 3650. , 3690. , 3730. , 3770. , 3810. , 3850. , 3890. ,
           3930. , 3970. , 4010. , 4050. , 4090. , 4130. , 4170. , 4210. , 4250. ,
           4290. , 4330. , 4370. , 4410. , 4450. , 4490. , 4530. , 4570. , 4610. ,
           4650. , 4690. , 4730. , 4770. , 4810. , 4850. , 4890. , 4930. , 4970. ,
           5010. , 5050. , 5090. , 5130. , 5170. , 5210. , 5250. , 5290. , 5330. ,
           5370. , 5410. , 5450. , 5490. , 5530. , 5580. , 5650. , 5730. , 5810. ,
           5890. , 5970. , 6050. , 6130. , 6210. , 6290. , 6370. , 6450. , 6530. ,
           6610. , 6690. , 6770. , 6850. , 6920. , 6975. ])

  • time
    (time)
    datetime64[ns]
    2020-03-13T18:00:00

    long_name :

    time

    array(['2020-03-13T18:00:00.000000000'], dtype='datetime64[ns]')

• Data variables: (21)
  • pa
    (time, height)
    float64
    ...

    units :

    Pa

    long_name :

    air_pressure

    [159 values with dtype=float64]

  • ua
    (time, height, x, y)
    float64
    ...

    units :

    m s-1

    long_name :

    eastward_wind

    [10420224 values with dtype=float64]

  • va
    (time, height, x, y)
    float64
    ...

    units :

    m s-1

    long_name :

    northward_wind

    [10420224 values with dtype=float64]

  • wa
    (time, height, x, y)
    float64
    ...

    units :

    m s-1

    long_name :

    upward_air_velocity

    [10420224 values with dtype=float64]

  • ta
    (time, height, x, y)
    float64
    ...

    units :

    K

    long_name :

    air_temperature

    [10420224 values with dtype=float64]

  • rhoa
    (time, height, x, y)
    float64
    ...

    units :

    kg m-3

    long_name :

    air_dry_density

    [10420224 values with dtype=float64]

  • qv
    (time, height, x, y)
    float64
    ...

    units :

    kg kg-1

    long_name :

    water_vapor_mixing_ratio

    [10420224 values with dtype=float64]

  • qlc
    (time, height, x, y)
    float64
    ...

    units :

    kg kg-1

    long_name :

    mass_mixing_ratio_of_cloud_liquid_water_in_air

    [10420224 values with dtype=float64]

  • qlr
    (time, height, x, y)
    float64
    ...

    units :

    kg kg-1

    long_name :

    mass_mixing_ratio_of_rain_water_in_air

    [10420224 values with dtype=float64]

  • qic
    (time, height, x, y)
    float64
    ...

    units :

    kg kg-1

    long_name :

    mass_mixing_ratio_of_cloud_ice_in_air

    [10420224 values with dtype=float64]

  • qis
    (time, height, x, y)
    float64
    ...

    units :

    kg kg-1

    long_name :

    mass_mixing_ratio_of_snow_in_air

    [10420224 values with dtype=float64]

  • qig
    (time, height, x, y)
    float64
    ...

    units :

    kg kg-1

    long_name :

    mass_mixing_ratio_of_graupel_in_air

    [10420224 values with dtype=float64]

  • nlc
    (time, height, x, y)
    float64
    ...

    units :

    kg-1

    long_name :

    number_of_cloud_droplets_in_air

    [10420224 values with dtype=float64]

  • nlr
    (time, height, x, y)
    float64
    ...

    units :

    kg-1

    long_name :

    number_of_rain_droplets_in_air

    [10420224 values with dtype=float64]

  • nic
    (time, height, x, y)
    float64
    ...

    units :

    kg-1

    long_name :

    number_of_cloud_ice_crystals_in_air

    [10420224 values with dtype=float64]

  • nis
    (time, height, x, y)
    float64
    ...

    units :

    kg-1

    long_name :

    number_of_snow_crystals_in_air

    [10420224 values with dtype=float64]

  • nig
    (time, height, x, y)
    float64
    ...

    units :

    kg-1

    long_name :

    number_of_graupel_crystals_in_air

    [10420224 values with dtype=float64]

  • eps
    (time, height, x, y)
    float64
    ...

    units :

    m2 s-3

    long_name :

    dissipation_rate_of_turbulent_kinetic_energy

    [10420224 values with dtype=float64]

  • tend_cond
    (time, height, x, y)
    float64
    ...

    units :

    s-1

    long_name :

    tendency_of_condensation

    [10420224 values with dtype=float64]

  • tend_dep
    (time, height, x, y)
    float64
    ...

    units :

    s-1

    long_name :

    tendency_of_deposition

    [10420224 values with dtype=float64]

  • tend_rim
    (time, height, x, y)
    float64
    ...

    units :

    s-1

    long_name :

    tendency_of_riming

    [10420224 values with dtype=float64]

• Indexes: (4)
  • x
    PandasIndex

    PandasIndex(Float64Index([-12800.0, -12700.0, -12600.0, -12500.0, -12400.0, -12300.0,
                  -12200.0, -12100.0, -12000.0, -11900.0,
                  ...
                   11800.0,  11900.0,  12000.0,  12100.0,  12200.0,  12300.0,
                   12400.0,  12500.0,  12600.0,  12700.0],
                 dtype='float64', name='x', length=256))

  • y
    PandasIndex

    PandasIndex(Float64Index([-12800.0, -12700.0, -12600.0, -12500.0, -12400.0, -12300.0,
                  -12200.0, -12100.0, -12000.0, -11900.0,
                  ...
                   11800.0,  11900.0,  12000.0,  12100.0,  12200.0,  12300.0,
                   12400.0,  12500.0,  12600.0,  12700.0],
                 dtype='float64', name='y', length=256))

  • height
    PandasIndex

    PandasIndex(Float64Index([  10.0,   32.5,   60.0,   92.5,  130.0,  170.0,  210.0,  250.0,
                   290.0,  330.0,
                  ...
                  6290.0, 6370.0, 6450.0, 6530.0, 6610.0, 6690.0, 6770.0, 6850.0,
                  6920.0, 6975.0],
                 dtype='float64', name='height', length=159))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2020-03-13 18:00:00'], dtype='datetime64[ns]', name='time', freq=None))

• Attributes: (10)

  title :

  DHARMA LES results for COMBLE-MIP case: fixed Nd and Ni

  reference :

  https://github.com/ARM-Development/comble-mip

  authors :

  Florian Tornow (florian.tornow@nasa.gov) and Ann Fridlind (ann.fridlind@nasa.gov)

  version :

  2024-09-25 16:28:30

  format_version :

  DEPHY-derivative

  script :

  convert_DHARMA_LES_output_to_dephy_format_3D.ipynb

  startDate :

  2020-03-12 22:00:00.0

  force_geo :

  1

  surfaceType :

  ocean (after spin-up)

  surfaceForcing :

  ts (after spin-up)

(dephy_check['ua'].mean(dim=['x','y'])-dharma_params['translate'].u).plot()

dharma_params['translate'].v

qlr_3d
nlr_3d

array([[[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]],

       [[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]],

       [[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]],

       ...,

       [[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]],

       [[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]],

       [[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]]])

rho_i = 9.e2                       # general density of ice used here, kg/m3
alpha = max([2.2e4,-4.99e3-4.94e4*np.log10(iwc_s*1e3)])
mass_s = np.pi*4.*rho_i/alpha**3   # kg (based on distribution, half of value used in code as 'typical')
num_s = iwc_s/mass_s               # m-3 total and m-4 leading constant
N_icef(ix,iy,iz) = num_s*1e-3      # m-3*1e-3 --> L-1
num_lrg_s = num_s*(1 - igamma(5,alpha*d_thresh)/gamma(5))  # m-3
N_100_ice(ix,iy,iz) = N_100_ice(ix,iy,iz) + num_lrg_s*1e-3 # L-1
area_s = np.pi/4*num_s*720/(alpha**2*24.)                  # m2/m3
R2_ice(ix,iy,iz) = R2_ice(ix,iy,iz) + area_s*1e9/np.pi     # m2/m3 --> um2/L

z_s = num_s*gamma(8)/alpha**6*f_refr*1e18 # mm6/m3
vz_s = z_s*vel_s*(0.3/rhobar(iz))**0.5    # mm6/m3 * m/s estimated

  Cell In[2], line 5
    N_icef(ix,iy,iz) = num_s*1e-3      # m-3*1e-3 --> L-1
    ^
SyntaxError: cannot assign to function call here. Maybe you meant '==' instead of '='?

ng_ice = 3
do_rain = True
do_morr = True


nd_thresh = 1e-9
ni_thresh = 1e-12

## rewrite to gather from dharma.cdf
use_s08_gamr = 0
use_s08_limr = 0
disp_type = 0
droplet_con = 20.0
droplet_sig = 1.2
is_qr = 1


varnames = ['w','LWC','RWC','Nr','mu_r','Nwr','Dmr','D0r','Dmtr','lamr',$
            'dBZ','dBZr','DopV', 'DopVr','PR','PRr']
varnames = ['Nr','Dmtr','lamr']
    
vname      = [ 'w', moisture, 'qc' ]
i_w        = 0
i_moisture = 1
i_qc       = 2
i_last     = i_qc

vname  = [ vname, 'qic', 'qif', 'qid' ]
i_qic  = i_last + 1
i_qif  = i_last + 2
i_qid  = i_last + 3
i_last = i_last + 3

vname  = [ vname, 'qr' ]
i_qr   = i_last + 1
i_last = i_last + 1

vname = [ vname, 'nc', 'nr', 'nic', 'nif', 'nid' ]
i_nc  = i_last + 1
i_nr  = i_last + 2
i_nic = i_last + 3
i_nif = i_last + 4
i_nid = i_last + 5


vname = [ 'qr', 'nr' ]
i_qr = 0
i_nr = 1

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[1], line 1
----> 1 gamma

NameError: name 'gamma' is not defined

        self.vel_param_a = {'cl': 3e7, 'ci': 700., 'pl': 841.997, 'pi': 11.72}
        self.vel_param_b = {'cl': 2. * ureg.dimensionless,
                            'ci': 1. * ureg.dimensionless,
                            'pl': 0.8 * ureg.dimensionless,
                            'pi': 0.41 * ureg.dimensionless}
        
        v_tmp = model.vel_param_a[hyd_type] * p_diam ** model.vel_param_b[hyd_type]