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))
mask = mask & (alti[:, 1:] >= hght) & (ref_lev < nlev - 1)
t, l, i, j = np.where(mask)
fk = (hght[t, 0, i, j] - alti[:, :-1][mask]) \
/ (alti[:, 1:][mask] - alti[:, :-1][mask])
kzzz[mask] = kzzz[mask] * fk + (1 - fk) * vkminup
# else
dzz = np.diff(alti, axis=1)
ss = 1e-6 + ((winz[:, 1:] - winz[:, :-1]) ** 2
+ (winm[:, 1:] - winm[:, :-1]) ** 2) / (dzz * dzz)
rig = g * (po[:, 1:] - po[:, :-1]) / (dzz * ss) / th[:, [0]]
alph = 0.
chi = 0.
mask = (rh[:, 1:] > crhx) | (rh[:, :-1] > crhx)
qbar = 0.5 * (sphu[:, 1:][mask] + sphu[:, :-1][mask])
tbar = 0.5 * (temp[:, 1:][mask] + temp[:, :-1][mask])
alph = Lv * qbar / (R * tbar)
chi = Lv * Lv * qbar / (Cp * Rv * tbar * tbar)
rig[mask] = (1 + alph) \
* (rig[mask] - (g * g * (chi - alph)) / (1 + chi) / ss[mask] / Cp / tbar)
dk = np.sqrt(ss[mask]) / (1. / (vkarm * alti[:, 1:][mask]) + 1. / rlam) ** 2
upkz = np.where(
rig[mask] < 0,
dk * (1. - 8. * rig[mask] / (1. + 1.286 * np.sqrt(-rig[mask]))),
dk / (1. + 5 * rig[mask]) ** 2
)
kzzz[mask] = np.minimum(vkmax, np.maximum(vkminup, upkz))
# Save output
ref_dataarray = 0. * inout_datastore["inputs"][("meteo", "winm")][ddi]["spec"]
inout_datastore["outputs"][("meteo", "kzzz")][ddi]["spec"] = ref_dataarray + kzzz
