Example: convert DEPHY forcing to WRF-LES forcing

Example: convert DEPHY forcing to WRF-LES forcing#

This code reads the COMBLE-MIP LES/SCM forcing file that is supplied in the DEPHY format and writes input files that are formatted for WRF-LES (ASCII and NetCDF)

Contributed by Tim Juliano from NCAR, last updated on 12/1/23

Import libraries#

import numpy as np
import pandas as pd
import xarray as xr
from netCDF4 import Dataset, stringtochar
import os
import datetime as dt

Specify local write directory#

If on the ARM JupyterHub, it is recommended to create and specify a local directory that is outside of the COMBLE-MIP repository to write input files in your model’s format

If on the ARM JupyterHub, these files can be conveniently downloaded to your local machine and will remain when you log out

my_write_dir = '../../../my_input_files/'

Create input files for WRF-LES in ASCII or NetCDF format#

Write initial sounding, ozone profile, forcing and vertical grid files

Write input sounding file for WRF LES in ASCII format#

# use a list to accumulate output ASCII lines
lines_sounding = []

# number of lines to be read
n_z = dephy.dims['lev']

# arrays to be reported in columns
z_arr = dephy.coords['lev'].values
th_arr = dephy['theta'].values.squeeze() # K
qt_arr = dephy['qv'].values.squeeze()*1000. # g/kg
u_arr = dephy['u'].values.squeeze() # m/s
v_arr = dephy['v'].values.squeeze() # m/s
ps = dephy['ps'].values.squeeze()/100. # hPa
thetas = dephy['thetas'].values.squeeze() # K

# add a surface point (required by WRF code)
vars_arr = np.array([ps,thetas,qt_arr[0]])
str_firstline = np.array2string(vars_arr,formatter={'float_kind':lambda vars_arr:"%11.5f" % vars_arr})[1:-1]
lines_sounding.append(str_firstline)

# do profile
for kk in range(n_z):
    vars_arr = np.array([z_arr[kk],th_arr[kk],qt_arr[kk],u_arr[kk],v_arr[kk]])
    str_vars = np.array2string(vars_arr,formatter={'float_kind':lambda vars_arr:"%11.5f" % vars_arr})[1:-1]
    lines_sounding.append(str_vars)

# write list contents to ASCII file
filename_sounding_LES = 'input_sounding'
with open(filename_sounding_LES,mode='wt',encoding='utf-8') as sounding_file:
    sounding_file.write('\n'.join(lines_sounding))

Write ozone file for RRTMG#

Pressure file#

# use a list to accumulate output ASCII lines
lines_pressure = []

# arrays to be reported in columns
p_arr = dephy['pressure'].values.squeeze() # Pa
p_arr = p_arr[::-1]/100. # hPa

# do profile
for kk in range(n_z):
    vars_arr = np.array([p_arr[kk]])
    str_vars = np.array2string(vars_arr,formatter={'float_kind':lambda vars_arr:"%11.4f" % vars_arr})[1:-1]
    lines_pressure.append(str_vars)

# write list contents to ASCII file
filename_pressure_LES = 'ozone_plev.formatted'
with open(filename_pressure_LES,mode='wt',encoding='utf-8') as pressure_file:
    pressure_file.write('\n'.join(lines_pressure))

Ozone file#

# use a list to accumulate output ASCII lines
lines_ozone = []

n_months = 12
n_lat = 64

# arrays to be reported in columns
o3_arr = dephy['o3'].values.squeeze() # kg kg-1
o3_arr = o3_arr[::-1]

# do profile
for ii in range(n_months):
    for jj in range(n_lat):
        for kk in range(n_z):
            vars_arr = np.array([o3_arr[kk]])
            str_vars = np.array2string(vars_arr,formatter={'float_kind':lambda vars_arr:"%.4E" % vars_arr})[1:-1]
            lines_ozone.append(str_vars)

# write list contents to ASCII file
filename_ozone_LES = 'ozone.formatted'
with open(filename_ozone_LES,mode='wt',encoding='utf-8') as ozone_file:
    ozone_file.write('\n'.join(lines_ozone))

Write forcing file#

Read forcing winds from forcing file#

u_force = dephy['ug'][:,:].values # m/s
v_force = dephy['vg'][:,:].values # m/s

Repeat z_ls 1D array to be 2D array for writing#

nhrs = dephy.dims['time']
z_ls_write = np.tile(z_arr,(np.shape(u_force)[0],1))[0:nhrs,:]

Compute wind forcing tendency#

delt = ((dephy['time'][1]-dephy['time'][0]).dt.seconds).values # time between profiles (s)

u_tend = np.zeros([np.shape(u_force)[0],np.shape(u_force)[1]])
v_tend = np.zeros([np.shape(u_force)[0],np.shape(u_force)[1]])
z_tend = np.zeros([np.shape(u_force)[0],np.shape(u_force)[1]])

u_tend[0:-1,:] = (u_force[1:,:] - u_force[0:-1,:]) / delt
u_tend[-1,:] = u_tend[-2,:]
v_tend[0:-1,:] = (v_force[1:,:] - v_force[0:-1,:]) / delt
v_tend[-1,:] = v_tend[-2,:]
z_tend[0:-1,:] = (z_ls_write[1:,:] - z_ls_write[0:-1,:]) / delt
z_tend[-1,:] = z_tend[-2,:]

Forcing NetCDF name and delete file if already exists#

savename = 'input_ls_forcing.nc'

if os.path.exists(savename):
    os.remove(savename)
    print('The file ' + savename + ' has been deleted successfully')

The file input_ls_forcing.nc has been deleted successfully

Create new netcdf file#

try: ncfile.close()  # just to be safe, make sure dataset is not already open.
except: pass
ncfile = Dataset('./' + savename,mode='w',format='NETCDF3_CLASSIC') 
#print(ncfile)

Create dimensions#

levs = np.shape(u_force)[1]
lev_dim = ncfile.createDimension('force_layers', levs)      # level axis
time_dim = ncfile.createDimension('Time', None)    # unlimited axis (can be appended to)
datestrlen = ncfile.createDimension('DateStrLen', 19)

Create global attributes#

ncfile.title='AUXILIARY FORCING FOR CRM'

Create variables#

Dimensions#

time = ncfile.createVariable('Times', 'S1', ('Time','DateStrLen',))

Add times#

for i in np.arange(len(dephy['time'].values)):
    yyyy = pd.to_datetime(dephy['time'].values[i]).year
    m = pd.to_datetime(dephy['time'].values[i]).month
    dd = pd.to_datetime(dephy['time'].values[i]).day
    hh = pd.to_datetime(dephy['time'].values[i]).hour
    hold = dt.datetime(yyyy,m,dd,hh)
    hold2 = hold.strftime('%Y-%m-%d_%H:%M:%S')
    
    time[i,:] = stringtochar(np.array(hold2, 'S19'))

Time-varying wind forcing profiles#

u_ls = ncfile.createVariable('U_LS', np.float32, ('Time','force_layers',),fill_value=-999.)
u_ls.FieldType = 104
u_ls.MemoryOrder = 'Z  '
u_ls.description = 'large scale zonal wind velocity'
u_ls.units = ''
u_ls.stagger = ''

v_ls = ncfile.createVariable('V_LS', np.float32, ('Time','force_layers',),fill_value=-999.)
v_ls.FieldType = 104
v_ls.MemoryOrder = 'Z  '
v_ls.description = 'large scale meridional wind velocity'
v_ls.units = ''
v_ls.stagger = ''

z_ls = ncfile.createVariable('Z_LS', np.float32, ('Time','force_layers',),fill_value=-999.)
z_ls.FieldType = 104
z_ls.MemoryOrder = 'Z  '
z_ls.description = 'height of forcing time series'
z_ls.units = ''
z_ls.stagger = ''

u_ls_tend = ncfile.createVariable('U_LS_TEND', np.float32, ('Time','force_layers',),fill_value=-999.)
u_ls_tend.FieldType = 104
u_ls_tend.MemoryOrder = 'Z  '
u_ls_tend.description = 'tendency of zonal wind'
u_ls_tend.units = ''
u_ls_tend.stagger = ''

v_ls_tend = ncfile.createVariable('V_LS_TEND', np.float32, ('Time','force_layers',),fill_value=-999.)
v_ls_tend.FieldType = 104
v_ls_tend.MemoryOrder = 'Z  '
v_ls_tend.description = 'tendency of meridional wind'
v_ls_tend.units = ''
v_ls_tend.stagger = ''

z_ls_tend = ncfile.createVariable('Z_LS_TEND', np.float32, ('Time','force_layers',),fill_value=-999.)
z_ls_tend.FieldType = 104
z_ls_tend.MemoryOrder = 'Z  '
z_ls_tend.description = 'height of forcing time series'
z_ls_tend.units = ''
z_ls_tend.stagger = ''

Write data#

u_ls[:] = u_force
v_ls[:] = v_force
z_ls[:] = z_ls_write
u_ls_tend[:] = u_tend
v_ls_tend[:] = v_tend
z_ls_tend[:] = z_tend

Close the file#

# first print the Dataset object to see what we've got
print(ncfile)
# close the Dataset.
ncfile.close(); print('Dataset is closed!')

<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF3_CLASSIC data model, file format NETCDF3):
    title: AUXILIARY FORCING FOR CRM
    dimensions(sizes): force_layers(136), Time(21), DateStrLen(19)
    variables(dimensions): |S1 Times(Time, DateStrLen), float32 U_LS(Time, force_layers), float32 V_LS(Time, force_layers), float32 Z_LS(Time, force_layers), float32 U_LS_TEND(Time, force_layers), float32 V_LS_TEND(Time, force_layers), float32 Z_LS_TEND(Time, force_layers)
    groups: 
Dataset is closed!

Write surface forcing file for WRF LES in netCDF format#

Read skin temperature from forcing file#

ts_arr = dephy['ts'][:].values # K

Compute skin temperature tendency#

ts_tend = np.zeros(nhrs)

ts_tend[0:-1] = (ts_arr[1:] - ts_arr[0:-1]) / delt
ts_tend[-1] = ts_tend[-2]

Forcing NetCDF name and delete file if already exists#

savename = 'input_sfc_forcing.nc'

if os.path.exists(savename):
    os.remove(savename)
    print('The file ' + savename + ' has been deleted successfully')

The file input_sfc_forcing.nc has been deleted successfully

Create new netcdf file#

try: ncfile.close()  # just to be safe, make sure dataset is not already open.
except: pass
ncfile = Dataset('./' + savename,mode='w',format='NETCDF3_CLASSIC') 
#print(ncfile)

Create dimensions#

time_dim = ncfile.createDimension('Time', None)    # unlimited axis (can be appended to)
datestrlen = ncfile.createDimension('DateStrLen', 19)

Create global attributes#

ncfile.title='AUXILIARY FORCING FOR CRM'

Create variables#

Dimensions#

time = ncfile.createVariable('Times', 'S1', ('Time','DateStrLen',))

Add times#

for i in np.arange(len(dephy['time'].values)):
    yyyy = pd.to_datetime(dephy['time'].values[i]).year
    m = pd.to_datetime(dephy['time'].values[i]).month
    dd = pd.to_datetime(dephy['time'].values[i]).day
    hh = pd.to_datetime(dephy['time'].values[i]).hour
    hold = dt.datetime(yyyy,m,dd,hh)
    hold2 = hold.strftime('%Y-%m-%d_%H:%M:%S')
    
    time[i,:] = stringtochar(np.array(hold2, 'S19'))

Time-varying skin temperature#

pre_tsk = ncfile.createVariable('PRE_TSK', np.float32, ('Time',),fill_value=-999.)
pre_tsk.FieldType = 104
pre_tsk.MemoryOrder = 'Z  '
pre_tsk.description = 'skin temperature'
pre_tsk.units = ''
pre_tsk.stagger = ''

pre_tsk_tend = ncfile.createVariable('PRE_TSK_TEND', np.float32, ('Time',),fill_value=-999.)
pre_tsk_tend.FieldType = 104
pre_tsk_tend.MemoryOrder = 'Z  '
pre_tsk_tend.description = 'skin temperature'
pre_tsk_tend.units = ''
pre_tsk_tend.stagger = ''

Write data#

pre_tsk[:] = ts_arr
pre_tsk_tend[:] = ts_tend

Close the file#

# first print the Dataset object to see what we've got
print(ncfile)
# close the Dataset.
ncfile.close(); print('Dataset is closed!')

<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF3_CLASSIC data model, file format NETCDF3):
    title: AUXILIARY FORCING FOR CRM
    dimensions(sizes): Time(21), DateStrLen(19)
    variables(dimensions): |S1 Times(Time, DateStrLen), float32 PRE_TSK(Time), float32 PRE_TSK_TEND(Time)
    groups: 
Dataset is closed!

Write vertical grid file for WRF LES in ASCII format#

# create vertical grid file contents as list
lines_zw = []

n_zw = n_z = dephy.dims['zw_grid'] # number of heights

for kk in range(n_zw):
    lines_zw.append(str(dephy['zw_grid'][kk].values))

# write vertical grid to ASCII file
filename_zw_LES = 'input_zw_grid'
with open(filename_zw_LES,mode='wt',encoding='utf-8') as zw_file:
    zw_file.write('\n'.join(lines_zw))

Example: convert DEPHY forcing to WRF-LES forcing

Contents

Example: convert DEPHY forcing to WRF-LES forcing#

Import libraries#

Specify local write directory#

Read the COMBLE-MIP forcing file#

Create input files for WRF-LES in ASCII or NetCDF format#

Write input sounding file for WRF LES in ASCII format#

Write ozone file for RRTMG#

Pressure file#

Ozone file#

Write forcing file#

Read forcing winds from forcing file#

Repeat z_ls 1D array to be 2D array for writing#

Compute wind forcing tendency#

Forcing NetCDF name and delete file if already exists#

Create new netcdf file#

Create dimensions#

Create global attributes#

Create variables#

Dimensions#

Add times#

Time-varying wind forcing profiles#

Write data#

Close the file#

Write surface forcing file for WRF LES in netCDF format#

Read skin temperature from forcing file#

Compute skin temperature tendency#

Forcing NetCDF name and delete file if already exists#

Create new netcdf file#

Create dimensions#

Create global attributes#

Create variables#

Dimensions#

Add times#

Time-varying skin temperature#

Write data#

Close the file#

Write vertical grid file for WRF LES in ASCII format#