I am trying to recreate the values of geopotential height in ERA5 that have been interpolated below model terrain, but am thus far unsuccessful. I am following the documentation in Chapter 3 “Vertical interpolations and extrapolations” in the latest documentation version (46t1r1), specifically equations 1, 3, and 6 for geopotential height.
As the attached figure shows, my recreated values (black) do not match the ERA5 values (red) below the surface (dashed horizontal line). My python code used to perform these calculations is shown below.
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib as mpl
# Read in ERA5 datasets -------------------------------------------------------------
era5_T = xr.open_dataset('./data/e5.oper.an.pl.128_130_t.ll025sc.2022010100_2022010123.nc').isel(time=12)
era5_Z = xr.open_dataset('./data/e5.oper.an.pl.128_129_z.ll025sc.2022010100_2022010123.nc').isel(time=12)
era5_PHIS = xr.open_dataset('./data/e5.oper.invariant.128_129_z.ll025sc.1979010100_1979010100.nc')
era5_PS = xr.open_dataset('./data/e5.oper.an.sfc.128_134_sp.ll025sc.2022010100_2022013123.nc').isel(time=12)
levels = era5_T.level
# Static Parameters -----------------------------------------------------------------
# Unit reminder: J = kg m**2 s**-2
dTdz_st = 0.0065 # units = K/m
Rd = 287.058 # units = J/(K*kg)
g = 9.80665 # units = m/(s**2)
cp = 1004 # units = J/(K*kg)
p0 = 100000 # units = Pa
# Define Exner function -------------------------------------------------------------
def exner(p_Pa):
return (p_Pa/p0)**(Rd/cp)
# Grid point-specific Parameters ----------------------------------------------------
Rlat, Rlon = 207, 1010 # 38.25°N, 107.5°W
Psrf = era5_PS['SP'][Rlat,Rlon].values # Sfc pres at grid point, units = Pa
PLe = (levels[levels < (Psrf/100)][-1].values)*100. # Lowest pressure level with real (non-interpolated) values, units = Pa
PIs = exner(Psrf) # Exner function value at surface pressure
phis = era5_PHIS['Z'][0,Rlat,Rlon].values # Sfc geopotential, units = m**2 s**-2
# Calculating Tsurf ------------------------------------------------------------------------
TLe = era5_T['T'][:,Rlat,Rlon].sel(level=(PLe/100)).values # Temperature at lowest pressure level with real (non-interpolated) values, units = K
Tsurf = TLe + (dTdz_st * (Rd / g) * ((PIs / exner(PLe)) - 1) * TLe) # Surface temperature, units = K (calculated with equation 1 on pg. 19)
# Calculate Zextrap for each pressure level below terrain -----------------------
plevs = (levels[levels > (Psrf/100)].values)*100 # all pres lvls below terrain, units = Pa
Zextrap = np.zeros(len(plevs)) # store calculated values in this numpy array
for i, p in enumerate(plevs):
# Calculate y
y = (dTdz_st) * (Rd / g) * np.log10((PIs / exner(p))) # unitless
# Calculate Zextrap, units = m**2 s**-2
Zextrap[i] = (phis) - ( Rd * Tsurf * np.log10((PIs / exner(p))) * (1 + (y/2) + ((y**2)/6)) )
