import numpy as np
[docs]
def vertical_turbulent_diffusivity(transf, inout_datastore, ddi, mapper):
"""Compute the vertical turbulent diffusivity Kzz at CHIMERE layer interfaces.
Inside the PBL, uses an O'Brien-type parabolic profile scaled by a
convective velocity that depends on stability, with a minimum diffusivity
enhanced in cloudy layers (using ``rhmaxx`` from
:func:`~pycif.plugins.transforms.complex.diagmet.utils.low_cloud_top.low_cloud_top`).
Above the PBL, uses a mixing-length closure with a moist-corrected gradient
Richardson number and the Louis (1982) stability correction. See
:doc:`/documentation/doc-models/chimere/diagmet` (section 11) for the full
derivation.
Args:
transf (Plugin): diagmet transform instance.
inout_datastore (dict): mutable datastore.
ddi (datetime): current sub-simulation date.
mapper (dict): transform mapper.
"""
# Parameters
vkmindry = 0.1 # Minimum Kz in the dry boundary layer (m2/s)
vkminwet = 5.0 # Minimum Kz in cloudy boundary layer (m2/s)
vkminup = 0.1 # Minimum Kz above PBL (m2/s)
vkmax = 500. # Maximum Kz
crhx = 0.90 # Min RH for cloud BLH enhancement
vkarm = 0.4 # von Karman constant
g = 9.81 # Gravity constant
R = 287.04 # R constant
Cp = 1005. # Cp
Lv = 2.45e6 # Lv
Rv = 461.5 # Rv
rlam = 150. # Upper air mixing length
# Inputs
hght = inout_datastore["inputs"][("meteo", "hght")][ddi]["spec"].values
ustar = inout_datastore["outputs"][("meteo", "usta")][ddi]["spec"].values
wstar = inout_datastore["outputs"][("meteo", "wsta")][ddi]["spec"].values
obuklen = inout_datastore["outputs"][("meteo", "obuk")][ddi]["spec"].values
rhmaxx = transf.diag_misc["rhmaxx"]
alti = transf.diag_misc["alti"]
temp = transf.diag_misc["temp"]
sphu = transf.diag_misc["sphu"]
po = transf.diag_misc["po"]
th = transf.diag_misc["th"]
rh = transf.diag_misc["rh"]
winm = transf.diag_misc["winm"]
winz = transf.diag_misc["winz"]
# Compute
ntime, nlev, nlat, nlon = alti[:, :-1].shape
ref_t, ref_lev, ref_i, ref_j = np.meshgrid(np.arange(ntime), np.arange(nlev),
np.arange(nlat), np.arange(nlon),
indexing="ij")
vkminbl = vkmindry + (vkminwet - vkmindry) * (rhmaxx - crhx) / (1. - crhx)
zn = alti[:, 1:] / hght
ep = np.minimum(0.1, zn)
zsl = zn * hght / obuklen
wc = np.where(obuklen > 0,
ustar / (1. + 4.7 * zsl),
(ustar ** 3 + ep * 2.8 * wstar ** 3) ** 0.3333)
# kzzz
kzzz = 0. * zn
# zn < 1
mask = zn <= 1
t, l, i, j = np.where(mask)
kzzz[mask] = np.minimum(
vkmax, np.maximum(vkminbl[t, i, j],
vkarm * wc[mask] * hght[t, 0, i, j]
* zn[mask] * (1 - zn[mask]) ** 2))
# Kz at interface l (top of layer l+1, alti index l+1) is computed from
# the levels l+1 and l+2 bracketing it; the top interface reuses the last
# available pair
iup = np.minimum(np.arange(2, nlev + 2), nlev)
ilo = iup - 1
# Blend towards vkminup if the next interface is above the PBL top
alti_up = alti[:, iup]
alti_lo = alti[:, ilo]
mask = mask & (alti_up >= hght) & (ref_lev < nlev - 1)
t, l, i, j = np.where(mask)
fk = (hght[t, 0, i, j] - alti_lo[mask]) / (alti_up[mask] - alti_lo[mask])
kzzz[mask] = kzzz[mask] * fk + (1 - fk) * vkminup
# zn > 1: Richardson closure everywhere, moist correction only in humid layers
dzz = alti_up - alti_lo
ss = 1e-6 + ((winz[:, iup] - winz[:, ilo]) ** 2
+ (winm[:, iup] - winm[:, ilo]) ** 2) / (dzz * dzz)
rig = g * (po[:, iup] - po[:, ilo]) / (dzz * ss) / th[:, [0]]
moist = (rh[:, iup] > crhx) | (rh[:, ilo] > crhx)
qbar = 0.5 * (sphu[:, iup][moist] + sphu[:, ilo][moist])
tbar = 0.5 * (temp[:, iup][moist] + temp[:, ilo][moist])
alph = Lv * qbar / (R * tbar)
chi = Lv * Lv * qbar / (Cp * Rv * tbar * tbar)
rig[moist] = (1 + alph) \
* (rig[moist] - (g * g * (chi - alph)) / (1 + chi) / ss[moist] / Cp / tbar)
dk = np.sqrt(ss) / (1. / (vkarm * alti[:, 1:]) + 1. / rlam) ** 2
upkz = np.where(
rig < 0,
dk * (1. - 8. * rig / (1. + 1.286 * np.sqrt(np.maximum(-rig, 0.)))),
dk / (1. + 5 * np.maximum(rig, 0.)) ** 2
)
above = zn > 1
kzzz[above] = np.minimum(vkmax, np.maximum(vkminup, upkz[above]))
# Save output
ref_dataarray = 0. * inout_datastore["inputs"][("meteo", "winm")][ddi]["spec"]
inout_datastore["outputs"][("meteo", "kzzz")][ddi]["spec"] = ref_dataarray + kzzz
