Hi,
I would need a way to get vertical mass fluxes (in kg/m2/s) of air through w-levels of WRF, consistent with mass conservation.
I try to figure out how well WRF output fields satisfy the continuity equation for air-mass.
I have took a WRF output from AMPS (any other WRF output would do):
And run the following code that integrates the continuity equation along vertical to get air flux divergence across half of the domain, and across full vertical of the domain.
The half-domain flux divergence ( or mass flux across the w level in the middle of the domain) looks like
For the whole domain, I would expect it to vanish, but it stays well beyond the numerical errors.
I wonder if it is some bug in my script to calculate it, my misunderstanding of continuity in WRF, my overoptimism about numerical errors, or some issue with conserving air mass in the model?
Any help would be appreciated.
UPD 2024-05-31: Corrected the magic to download the file
I would need a way to get vertical mass fluxes (in kg/m2/s) of air through w-levels of WRF, consistent with mass conservation.
I try to figure out how well WRF output fields satisfy the continuity equation for air-mass.
I have took a WRF output from AMPS (any other WRF output would do):
Code:
ncks -v U,V,XLAT,XLONG,PSFC,P_TOP,ZNW,P_STRAT http://amrdcdata.ssec.wisc.edu/thredds/dodsC/yoppshwrfnetcdf/202206/2022060400/wrfout_d01_2022060400_f066.nc wrfout_d01_2022060400_f066.ncAnd run the following code that integrates the continuity equation along vertical to get air flux divergence across half of the domain, and across full vertical of the domain.
Python:
#!/usr/bin/env python3
import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as plt
import netCDF4 as nc4
import datetime as dt
import numpy as np
import os
import sys
def cmap_and_norm_clevs(cmapname, clevs):
# Input -- name of mpl colormap
# np.array of cleca
#Generates colormap and norm with given clevs
# As in grads
basecm=mpl.cm.get_cmap(cmapname)
nlevs=len(clevs)
cols=basecm(np.linspace(0, 1, nlevs+1))
norm = mpl.colors.BoundaryNorm(clevs, nlevs-1)
#cm_out = pylab.cm.from_list("%s_discr"%cmapname, cols[1:-1], nlevs-1)
cm_out = mpl.colors.ListedColormap(cols[1:-1], name="%s_discr"%cmapname)
cm_out.set_under(cols[0])
cm_out.set_over(cols[-1])
cm_out.set_bad([0,0,0,0]) #Transparent
return cm_out, norm
infname="wrfout_d01_2022060400_f066.nc"
with nc4.Dataset(infname) as inf:
inf.set_auto_mask(False) ## No masked arrays, please!
nz = inf.dimensions['bottom_top'].size
it = 0
dx = inf.DX
dy = inf.DY
u = inf.variables['U'][it,:,:,:]
v = inf.variables['V'][it,:,:,:]
eta = inf.variables['ZNW'][it,:]
ptop = inf.variables['P_TOP'][it] #P_TOP:description = "PRESSURE TOP OF THE MODEL" ;
psfc = inf.variables['PSFC'][it,:,:] # PSFC:description = "SFC PRESSURE" ;
# zeroth-level skipped
psfc0 = (psfc[np.newaxis,:,:] - ptop)
psfcu = np.zeros(u.shape[1:])
psfcv = np.zeros(v.shape[1:])
# remap surface pressure to U and V
psfcu[:,0] = psfc[:,0]
psfcu[:,-1] = psfc[:,-1]
psfcu[:,1:-1] = 0.5*(psfc[:,:-1] + psfc[:,1:])
psfcv[0,:] = psfc[0,:]
psfcv[-1,:] = psfc[-1,:]
psfcv[1:-1,:] = 0.5*(psfc[:-1,:] + psfc[1:,:])
u *= psfcu / 9.81 ## kg/m2/s
v *= psfcv / 9.81 ## kg/m2/s
div = (u[:,:,:-1] - u[:,:,1:] +
v[:,:-1,:] - v[:,1:,:]) / dx
deta = eta[1:] - eta[:-1]
div *= deta[:,np.newaxis,np.newaxis] # rho*dz
#grid%pb(i,k,j) = grid%znu(k)*(p_surf - p_top) + p_top
fig = plt.figure()
clevs = [-2, -1, -.5, -.2, -.1, -.05, .05, .1, .2, .5, 1, 2]
cmap,norm = cmap_and_norm_clevs("jet", clevs)
for iLev in [ nz//2, nz-1 ]:
ax = fig.add_axes((.1,.1,.8,.8)) #(left, bottom, width, height)
ax.imshow(np.sum(div[0:(iLev+1),:,:], axis=0), cmap=cmap, norm=norm)
ax.set_title("vertical mass flux at w level %d, kg/m2/s"%(iLev))
axcb = fig.add_axes((.85,.1,.03,.8))
cb = mpl.colorbar.ColorbarBase(axcb, norm=norm, cmap=cmap,
orientation='vertical', extend='both')
axcb.set_yticks(clevs)
figname = "div_lev_0-%02d.png"%(iLev)
fig.savefig(figname, dpi=600)
fig.clf()For the whole domain, I would expect it to vanish, but it stays well beyond the numerical errors.
I wonder if it is some bug in my script to calculate it, my misunderstanding of continuity in WRF, my overoptimism about numerical errors, or some issue with conserving air mass in the model?
Any help would be appreciated.
UPD 2024-05-31: Corrected the magic to download the file
Last edited:
