Example: convert DHARMA LES 2D output to DEPHY format

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_sbhomo'
my_outfile = 'DHARMA_Lx25_dx100_FixN_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')]

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)

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

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

xy_res = Lx/nx
xy_res

100.0

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')
    opd_ice   = read_2d(FILE=input_filename,vbase = 'opd_ice')
    alb       = (opd_cloud + opd_ice)/(13 + opd_cloud + opd_ice)
    opd       = opd_cloud + opd_ice
    
    lwp       = read_2d(FILE=input_filename,vbase = 'lwp')/1000
    iwp       = read_2d(FILE=input_filename,vbase = 'iwp')/1000
    pflux     = read_2d(FILE=input_filename,vbase = 'pflux')
    sfcTflx   = read_2d(FILE=input_filename,vbase = 'sfc_T_flx')
    sfcQflx   = read_2d(FILE=input_filename,vbase = 'sfc_qv_flx')
    u10m      = read_2d(FILE=input_filename,vbase = 'u_layer1') + dharma_params['translate'].u
    v10m      = read_2d(FILE=input_filename,vbase = 'v_layer1') + dharma_params['translate'].v
    LWup      = read_2d(FILE=input_filename,vbase = 'Fup_toa_10um')
    
    ## 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.6 | 0.2410E-07 | 0.4196E-08 | 0.3253 | 0.6000 |

|   2880 |    7200.0 |      0.6 | 0.2695E-07 | 0.4070E-08 | 0.3137 | 0.6000 |

|   4320 |   10800.0 |      0.7 | 0.2916E-07 | 0.3844E-08 | 0.3881 | 0.6000 |

|   5760 |   14400.0 |      0.7 | 0.3864E-07 | 0.3393E-08 | 0.5640 | 0.6000 |

|   7200 |   18000.0 |      0.7 | 0.2956E-07 | 0.2480E-08 | 0.4532 | 0.6000 |

|   8640 |   21600.0 |      0.7 | 0.3239E-07 | 0.2085E-08 | 0.5125 | 0.6000 |

|  10080 |   25200.0 |      0.7 | 0.2713E-07 | 0.1618E-08 | 0.5182 | 0.6000 |

|  11520 |   28800.0 |      0.7 | 0.2600E-07 | 0.1432E-08 | 0.5423 | 0.6000 |

|  12960 |   32400.0 |      0.7 | 0.2968E-07 | 0.1463E-08 | 0.4731 | 0.6000 |

|  14400 |   36000.0 |      0.7 | 0.2777E-07 | 0.1661E-08 | 0.4961 | 0.6000 |

|  15841 |   39600.0 |      0.7 | 0.2618E-07 | 0.1788E-08 | 0.2711 | 0.6000 |

|  17282 |   43200.0 |      0.7 | 0.2922E-07 | 0.1958E-08 | 0.1285 | 0.6000 |

|  18723 |   46800.0 |      0.7 | 0.2857E-07 | 0.2273E-08 | 0.1562 | 0.6000 |

|  20199 |   50400.0 |      0.7 | 0.3986E-07 | 0.2530E-08 | 0.1862 | 0.6000 |

|  21679 |   54000.0 |      0.8 | 0.3477E-07 | 0.2724E-08 | 0.2045 | 0.6000 |

|  23194 |   57600.0 |      0.7 | 0.3028E-07 | 0.4011E-08 | 0.4434 | 0.6000 |

|  24678 |   61200.0 |      0.8 | 0.3973E-07 | 0.4202E-08 | 0.2893 | 0.6000 |

|  26199 |   64800.0 |      0.7 | 0.3157E-07 | 0.4565E-08 | 0.2191 | 0.6000 |

|  27650 |   68400.0 |      0.8 | 0.4095E-07 | 0.5086E-08 | 0.3644 | 0.6000 |

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

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

# 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]
    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()

The file ./../../output_les/dharma/sandbox/DHARMA_Lx25_dx100_FixN_2D.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-23 22:05:59
    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 hfss(time, x, y), float64 hfls(time, x, y), 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_filename

'./../../output_les/dharma/sandbox/DHARMA_Lx25_dx100_FixN_2D.nc'

#dephy_check['olr11'].plot(row='time')
dephy_check['olr11'].plot(row='time',col_wrap=7)

<xarray.plot.facetgrid.FacetGrid at 0x7fe5dac745b0>

../../_images/b752bda65258f83332dde79668792470d57435adad6929900e52912bf63351cf.png