Convert DEPHY forcing to UCLALES-SALSA forcing

Contents

Convert DEPHY forcing to UCLALES-SALSA forcing#

Based on the example 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
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 = '../../../inputs/'

Create input files for UCLALES-SALSA in ASCII format#

Write input sounding file#

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

# number of lines to be read
n_z = dephy.dims['lev']
n_z_les = n_z # Values up to 7000 m

# arrays to be reported in columns
z_arr = dephy.coords['lev'].values
p_arr = dephy['pressure'].values.squeeze()/100. # hPa = mbar
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 = mbar
#thetas = dephy['thetas'].values.squeeze() # K

# The first points: surface pressure and then theta, qt, u, and v for the first level 
# *below* surface (surface is interpolated). thetas looks a bit low, so set the surface 
# to th_arr[0]
#vars_arr = np.array([ps,thetas*2-th_arr[1],qt_arr[0],u_arr[0],v_arr[0]]) # theta is low!
vars_arr = np.array([ps,th_arr[0],qt_arr[0],u_arr[0],v_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)
    if z_arr[kk]>7000.0:
        n_z_les=kk+1
        break

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

Grid coordinates#

# generate vertical grid file (zm_grid_in)
zw_arr = dephy['zw_grid'].values.squeeze() # m
n_zp = dephy.dims['zw_grid']
#
zw_file=open(my_write_dir+'zm_grid_in',mode='wt',encoding='utf-8')
for kk in dephy['zw_grid'].values:
    zw_file.write( str(kk) ) 
    zw_file.write('\n')
zw_file.close()

# also zt grid
zw_file=open(my_write_dir+'zt_grid_in',mode='wt',encoding='utf-8')
zw_file.write( str(zw_arr[0]-0.5*(zw_arr[1]-zw_arr[0])) )
for kk in range(n_zp-1):
    zw_file.write('\n')
    zw_file.write( str(0.5*(zw_arr[kk+1]+zw_arr[kk])) ) 
zw_file.close()

Write background profiles for radiative transfer#

# 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

temp_arr = dephy['temp'].values.squeeze() # K
temp_arr = temp_arr[::-1]

qt_arr = dephy['qv'].values.squeeze() # kg/kg
qt_arr = qt_arr[::-1]

o3_arr = dephy['o3'].values.squeeze() # kg kg-1
o3_arr = o3_arr[::-1]

# add a surface point
str_firstline = ("%.3f %s") % (temp_arr[n_z-1],n_z)
lines_pressure.append(str_firstline)

# do profile
for kk in range(n_z):
    str_vars = ("%11.6f %8.3f  %e  %e 0") % (p_arr[kk],temp_arr[kk],qt_arr[kk],o3_arr[kk])
    lines_pressure.append(str_vars)

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

Wind forcings#

# Read data
u_force = dephy['ug'][:,:].values # m/s
v_force = dephy['vg'][:,:].values # m/s
n_t = dephy.dims['time']

# Open files
forc_file_u=open(my_write_dir+'input_u_forcing.dat',mode='wt',encoding='utf-8')
forc_file_v=open(my_write_dir+'input_v_forcing.dat',mode='wt',encoding='utf-8')

# The first line: dimensions
forc_file_u.write(("%i %i ")%(n_z_les,n_t))
forc_file_v.write(("%i %i ")%(n_z_les,n_t))

# Note: time is obtained from the sst file

# Save the data values
for j in range(n_z_les):
    forc_file_u.write("\n%10.3f "%z_arr[j])
    forc_file_v.write("\n%10.3f "%z_arr[j])
    for i in range(n_t):
        forc_file_u.write("%10f "%u_force[i,j])
        forc_file_v.write("%10f "%v_force[i,j])
forc_file_u.close()
forc_file_v.close()

Surface temperature forcing#

# Read data
ts_arr = dephy['ts'][:].values # K
n_t = dephy.dims['time']

# Open file
sst_file=open(my_write_dir+'input_sst.dat',mode='wt',encoding='utf-8')

# The first line: dimensions
sst_file.write("%i"%n_t)

# Save the data values: elapsed time (s), sst (K)
for kk in range(n_t):
    delta = ((dephy['time'][kk]-dephy['time'][0]).dt.seconds).values # time between profiles (s)
    sst_file.write(("\n%i %10.3f")%(delta,ts_arr[kk]))
#
sst_file.close()

# All done
dephy.close()