Example: convert DHARMA LES 2D 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 function_extract_2d.py 

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_hm_largedomain'
my_outfile = 'DHARMA_Lx25_dx100_FixN_2D_large.nc'

# FixN with no ice test
my_readdir = 'case0313_diag_ice0_nomiz_dx100_specZ0'
my_outfile = 'DHARMA_Lx25_dx100_FixN_noice_2D.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')]

print(NCFILES_STR)

input_filenames = []
count = 0
for fn in NCFILES: 
    if ('alt_0' in NCFILES_STR[count]):
        input_filenames.append(fn)
    count += 1
    
input_filenames = sorted(input_filenames)
['/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_021984.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_010091.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_007202.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_001440.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_018909.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_028317.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_011548.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_017385.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_3d_029815.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_026734.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_PLT_029815.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_000000.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_014470.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_020442.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_004320.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_008645.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_015926.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_029815.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_005760.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_002880.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_023591.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_013007.tgz', '/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma_alt_025104.tgz']

Read 2D files#

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

nx = dharma_params['geometry'].nx
Lx = dharma_params['geometry'].L_x
ERROR 1: PROJ: proj_create_from_database: Open of /opt/conda/share/proj failed
/data/home/floriantornow/dharma/sandbox/case0313_diag_ice0_nomiz_dx100_specZ0/0-20h/dharma.cdf
xy_res = Lx/nx
xy_res
100.0
%run function_extract_2d.py 
time_vec = []

counter = 0
for input_filename in input_filenames:
    
    ## obtain timestamp
    time_vec.append(obtain_time(path,input_filename))

    ## 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)
        
    ## unpack tgz and extract relevant fields 
    opd_cloud = read_2d(FILE=input_filename,vbase = 'opd_drops',retain=True)
    opd_ice   = read_2d(FILE=input_filename,vbase = 'opd_ice',retain=True,there=True)
    alb       = (opd_cloud + opd_ice)/(13 + opd_cloud + opd_ice)
    opd       = opd_cloud + opd_ice
    
    lwp       = read_2d(FILE=input_filename,vbase = 'lwp',retain=True,there=True)/1000
    iwp       = read_2d(FILE=input_filename,vbase = 'iwp',retain=True,there=True)/1000
    pflux     = read_2d(FILE=input_filename,vbase = 'pflux',retain=True,there=True)
    sfcTflx   = read_2d(FILE=input_filename,vbase = 'sfc_T_flx',retain=True,there=True)
    sfcQflx   = read_2d(FILE=input_filename,vbase = 'sfc_qv_flx',retain=True,there=True)
    u10m      = read_2d(FILE=input_filename,vbase = 'u_layer1',retain=True,there=True) + dharma_params['translate'].u
    v10m      = read_2d(FILE=input_filename,vbase = 'v_layer1',retain=True,there=True) + dharma_params['translate'].v
    LWup      = read_2d(FILE=input_filename,vbase = 'Fup_toa_10um',there=True)
    
    ## accumulate
    if counter == 0:
        sfcTflx_col = sfcTflx
        sfcQflx_col = sfcQflx
        u10m_col    = u10m
        v10m_col    = v10m
        pflux_col   = pflux
        lwp_col     = lwp
        iwp_col     = iwp
        alb_col     = alb
        opd_col     = opd
        LWup_col    = LWup
    else:
        sfcTflx_col = np.append(sfcTflx_col,sfcTflx,axis=0)
        sfcQflx_col = np.append(sfcQflx_col,sfcQflx,axis=0)
        u10m_col    = np.append(u10m_col,u10m,axis=0)
        v10m_col    = np.append(v10m_col,v10m,axis=0)
        pflux_col   = np.append(pflux_col,pflux,axis=0)
        lwp_col     = np.append(lwp_col,lwp,axis=0)
        iwp_col     = np.append(iwp_col,iwp,axis=0)
        alb_col     = np.append(alb_col,alb,axis=0)
        opd_col     = np.append(opd_col,opd,axis=0)
        LWup_col    = np.append(LWup_col,LWup,axis=0)
    
    counter += 1
 
    
|      0 |       0.0 |      0.0 | 0.0000E+00 | 0.0000E+00 | 0.3255 | 0.0000 |

|   1440 |    3600.0 |      0.4 | 0.2395E-07 | 0.4198E-08 | 0.3247 | 0.6000 |

|   2880 |    7200.0 |      0.4 | 0.2652E-07 | 0.4070E-08 | 0.3141 | 0.6000 |

|   4320 |   10800.0 |      0.4 | 0.2791E-07 | 0.3844E-08 | 0.4053 | 0.6000 |

|   5760 |   14400.0 |      0.5 | 0.3063E-07 | 0.3442E-08 | 0.5153 | 0.6000 |

|   7202 |   18000.0 |      0.5 | 0.3884E-07 | 0.2439E-08 | 0.1469 | 0.6000 |

|   8645 |   21600.0 |      0.5 | 0.3184E-07 | 0.2072E-08 | 0.2324 | 0.6000 |

|  10091 |   25200.0 |      0.5 | 0.3594E-07 | 0.1698E-08 | 0.3456 | 0.6000 |

|  11548 |   28800.0 |      0.5 | 0.2809E-07 | 0.1481E-08 | 0.0009 | 0.0562 |

|  13007 |   32400.0 |      0.5 | 0.3115E-07 | 0.1557E-08 | 0.2232 | 0.6000 |

|  14470 |   36000.0 |      0.5 | 0.3068E-07 | 0.1722E-08 | 0.0551 | 0.6000 |

|  15926 |   39600.0 |      0.5 | 0.2986E-07 | 0.1897E-08 | 0.1307 | 0.6000 |

|  17385 |   43200.0 |      0.5 | 0.3778E-07 | 0.2163E-08 | 0.4093 | 0.6000 |

|  18909 |   46800.0 |      0.5 | 0.3518E-07 | 0.2405E-08 | 0.0846 | 0.6000 |

|  20442 |   50400.0 |      0.5 | 0.3853E-07 | 0.2736E-08 | 0.5808 | 0.6000 |

|  21984 |   54000.0 |      0.5 | 0.4572E-07 | 0.2869E-08 | 0.4125 | 0.6000 |

|  23591 |   57600.0 |      0.5 | 0.3958E-07 | 0.4079E-08 | 0.5121 | 0.6000 |

|  25104 |   61200.0 |      0.5 | 0.4455E-07 | 0.4341E-08 | 0.3821 | 0.6000 |

|  26734 |   64800.0 |      0.5 | 0.4428E-07 | 0.4690E-08 | 0.0968 | 0.6000 |

|  28317 |   68400.0 |      0.5 | 0.3679E-07 | 0.5194E-08 | 0.3057 | 0.6000 |

|  29815 |   72000.0 |      0.5 | 0.3417E-07 | 0.5125E-08 | 0.1688 | 0.6000 |
## convert narrowband Flux (spectrally integrated) into Brightness Temperatureh = 6.62607015*(10**-34) ## J Hz-1
h = 6.62607015*(10**-34) ## J Hz-1
k = 1.380649*(10**-23)   ## J⋅K−1
c = 2.99*(10**8)         ## m s-1
lamb = 10.5*(10**-6) 


BT_col = h*c/(k*lamb*np.log(1 + (2*h*(c**2)/((LWup_col/2.12/np.pi)*(10**6)*(lamb**5)))))
    

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

pd.set_option('display.max_rows', None)
vars_mean_list
standard_name variable_id units dimensions comment (10-min instantaneous)
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 surface_upward_sensible_heat_flux hfss W m-2 time, x, y
4 surface_upward_latent_heat_flux hfls W m-2 time, x, y
5 surface_friction_velocity ustar m s-1 time, x, y
6 surface_eastward_wind us m s-1 time, x, y at 10-m for comparison with SAR satellite meas...
7 surface_northward_wind vs m s-1 time, x, y at 10-m for comparison with SAR satellite meas...
8 precipitation_flux_at_surface pr kg m-2 s-1 time, x, y all hydrometeors
9 atmosphere_mass_content_of_liquid_water lwp kg m-2 time, x, y all liquid hydrometeors
10 atmosphere_mass_content_of_ice iwp kg m-2 time, x, y all ice hydrometeors
11 atmosphere_optical_thickness opt 1 time, x, y all hydrometeors for comparison with VIIRS; mi...
12 pseudo-albedo alb 1 time, x, y opt/(opt+13), where opt = all hydrometeor opti...
13 toa_outgoing_longwave_flux_atmospheric_window olr11 K time, x, y flux at top-of-atmosphere in window containing...
shape(alb_col)
shape(lwp_col)
(21, 256, 256)
# 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_2D.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

nt = 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[3:]:
    std_name = vars_mean_list.standard_name.iat[index]
    var_name = vars_mean_list['variable_id'].iat[index]
    #print(var_name)
    if (var_name == 'hfss') | (var_name == 'hfls'):
        continue
    if var_name == 'hfss':
        mod_field = sfcTflx_col
    if var_name == 'hfls':
        mod_field = sfcQflx_col
    if var_name == 'us':
        mod_field = u10m_col
    if var_name == 'vs':
        mod_field = v10m_col
    if var_name == 'pr':
        mod_field = pflux_col
    if var_name == 'lwp':
        mod_field = lwp_col
    if var_name == 'iwp':
        mod_field = iwp_col
    if var_name == 'opt':
        mod_field = opd_col
    if var_name == 'alb':
        mod_field = alb_col
    if var_name == 'olr11':
        mod_field = BT_col#LWup_col
    if var_name == 'ustar':
        continue
    
    if vars_mean_list.dimensions.iat[index]=='time, x, y':
        new_snd = dephy_file.createVariable(var_name, np.float64, ('time','x','y'))
        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[:] = mod_field

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: Florian Tornow (florian.tornow@nasa.gov) and Ann Fridlind (ann.fridlind@nasa.gov)
    version: 2024-12-18 21:39:36
    format_version: DEPHY-derivative
    script: convert_DHARMA_LES_output_to_dephy_format_2D.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), time(21)
    variables(dimensions): float64 x(x), float64 y(y), float64 time(time), float64 us(time, x, y), float64 vs(time, x, y), float64 pr(time, x, y), float64 lwp(time, x, y), float64 iwp(time, x, y), float64 opt(time, x, y), float64 alb(time, x, y), float64 olr11(time, x, y)
    groups: 
dephy_check = xr.open_dataset(dephy_filename)
dephy_check
<xarray.Dataset> Size: 88MB
Dimensions:  (x: 256, y: 256, time: 21)
Coordinates:
  * x        (x) float64 2kB -1.28e+04 -1.27e+04 -1.26e+04 ... 1.26e+04 1.27e+04
  * y        (y) float64 2kB -1.28e+04 -1.27e+04 -1.26e+04 ... 1.26e+04 1.27e+04
  * time     (time) datetime64[ns] 168B 2020-03-12T22:00:00 ... 2020-03-13T18...
Data variables:
    us       (time, x, y) float64 11MB ...
    vs       (time, x, y) float64 11MB ...
    pr       (time, x, y) float64 11MB ...
    lwp      (time, x, y) float64 11MB ...
    iwp      (time, x, y) float64 11MB ...
    opt      (time, x, y) float64 11MB ...
    alb      (time, x, y) float64 11MB ...
    olr11    (time, x, y) float64 11MB ...
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-12-18 21:39:36
    format_version:  DEPHY-derivative
    script:          convert_DHARMA_LES_output_to_dephy_format_2D.ipynb
    startDate:       2020-03-12 22:00:00.0
    force_geo:       1
    surfaceType:     ocean (after spin-up)
    surfaceForcing:  ts (after spin-up)
dephy_filename
'./../../output_les/dharma/sandbox/DHARMA_Lx25_dx100_FixN_noice_2D.nc'
#dephy_check['olr11'].plot(row='time')
dephy_check['alb'].plot(row='time',col_wrap=7)
<xarray.plot.facetgrid.FacetGrid at 0x7f65f39ec220>
../../_images/b51cbb3e62a34668419089892be3c5b3771b03a53c1c3aaf48a8d8ff0f4a5ca6.png