◆ satmedmfvdifq_run()

subroutine satmedmfvdifq::satmedmfvdifq_run	(	integer, intent(in)	im,
		integer, intent(in)	km,
		integer, intent(in)	ntrac,
		integer, intent(in)	ntcw,
		integer, intent(in)	ntrw,
		integer, intent(in)	ntiw,
		integer, intent(in)	ntke,
		real(kind=kind_phys), intent(in)	grav,
		real(kind=kind_phys), intent(in)	rd,
		real(kind=kind_phys), intent(in)	cp,
		real(kind=kind_phys), intent(in)	rv,
		real(kind=kind_phys), intent(in)	hvap,
		real(kind=kind_phys), intent(in)	hfus,
		real(kind=kind_phys), intent(in)	fv,
		real(kind=kind_phys), intent(in)	eps,
		real(kind=kind_phys), intent(in)	epsm1,
		real(kind=kind_phys), dimension(:,:), intent(inout)	dv,
		real(kind=kind_phys), dimension(:,:), intent(inout)	du,
		real(kind=kind_phys), dimension(:,:), intent(inout)	tdt,
		real(kind=kind_phys), dimension(:,:,:), intent(inout)	rtg,
		real(kind=kind_phys), dimension(:,:), intent(in)	u1,
		real(kind=kind_phys), dimension(:,:), intent(in)	v1,
		real(kind=kind_phys), dimension(:,:), intent(in)	t1,
		real(kind=kind_phys), dimension(:,:,:), intent(in)	q1,
		real(kind=kind_phys), dimension(:,:), intent(in)	swh,
		real(kind=kind_phys), dimension(:,:), intent(in)	hlw,
		real(kind=kind_phys), dimension(:), intent(in)	xmu,
		real(kind=kind_phys), dimension(:), intent(in)	garea,
		real(kind=kind_phys), dimension(:), intent(in)	zvfun,
		real(kind=kind_phys), dimension(:), intent(in)	sigmaf,
		real(kind=kind_phys), dimension(:), intent(in)	psk,
		real(kind=kind_phys), dimension(:), intent(in)	rbsoil,
		real(kind=kind_phys), dimension(:), intent(in)	zorl,
		real(kind=kind_phys), dimension(:), intent(in)	u10m,
		real(kind=kind_phys), dimension(:), intent(in)	v10m,
		real(kind=kind_phys), dimension(:), intent(in)	fm,
		real(kind=kind_phys), dimension(:), intent(in)	fh,
		real(kind=kind_phys), dimension(:), intent(in)	tsea,
		real(kind=kind_phys), dimension(:), intent(in)	heat,
		real(kind=kind_phys), dimension(:), intent(in)	evap,
		real(kind=kind_phys), dimension(:), intent(in)	stress,
		real(kind=kind_phys), dimension(:), intent(in)	spd1,
		integer, dimension(:), intent(out)	kpbl,
		real(kind=kind_phys), dimension(:,:), intent(in)	prsi,
		real(kind=kind_phys), dimension(:,:), intent(in)	del,
		real(kind=kind_phys), dimension(:,:), intent(in)	prsl,
		real(kind=kind_phys), dimension(:,:), intent(in)	prslk,
		real(kind=kind_phys), dimension(:,:), intent(in)	phii,
		real(kind=kind_phys), dimension(:,:), intent(in)	phil,
		real(kind=kind_phys), intent(in)	delt,
		logical, intent(in)	dspheat,
		real(kind=kind_phys), dimension(:), intent(out)	dusfc,
		real(kind=kind_phys), dimension(:), intent(out)	dvsfc,
		real(kind=kind_phys), dimension(:), intent(out)	dtsfc,
		real(kind=kind_phys), dimension(:), intent(out)	dqsfc,
		real(kind=kind_phys), dimension(:), intent(out)	hpbl,
		real(kind=kind_phys), dimension(:,:), intent(out)	dkt,
		real(kind=kind_phys), dimension(:,:), intent(out)	dku,
		integer, dimension(:), intent(in)	kinver,
		real(kind=kind_phys), intent(in)	xkzm_m,
		real(kind=kind_phys), intent(in)	xkzm_h,
		real(kind=kind_phys), intent(in)	xkzm_s,
		real(kind=kind_phys), intent(in)	dspfac,
		real(kind=kind_phys), intent(in)	bl_upfr,
		real(kind=kind_phys), intent(in)	bl_dnfr,
		real(kind=kind_phys), intent(in)	rlmx,
		real(kind=kind_phys), intent(in)	elmx,
		integer, intent(in)	sfc_rlm,
		integer, intent(in)	tc_pbl,
		integer, intent(in)	ntqv,
		real(kind=kind_phys), dimension(:,:,:), intent(inout)	dtend,
		integer, dimension(:,:), intent(in)	dtidx,
		integer, intent(in)	index_of_temperature,
		integer, intent(in)	index_of_x_wind,
		integer, intent(in)	index_of_y_wind,
		integer, intent(in)	index_of_process_pbl,
		logical, intent(in)	gen_tend,
		logical, intent(in)	ldiag3d,
		character(len=*), intent(out)	errmsg,
		integer, intent(out)	errflg
	)

Argument Table

GFS satmedmfvdifq General Algorithm

satmedmfvdifq_run() computes subgrid vertical turbulence mixing using the scale-aware TKE-based moist eddy-diffusion mass-flux (EDMF) parameterization of Han and Bretherton (2019) [80] .

The local turbulent mixing is represented by an eddy-diffusivity scheme which is a function of a prognostic TKE.
For the convective boundary layer, nonlocal transport by large eddies (mfpbltq.f), is represented using a mass flux approach (Siebesma et al.(2007) [178] ).
A mass-flux approach is also used to represent the stratocumulus-top-induced turbulence (mfscuq.f).

GFS satmedmfvdifq Detailed Algorithm

Compute preliminary variables from input arguments

Compute physical height of the layer centers and interfaces from the geopotential height (zi and zl)
Compute horizontal grid size (gdx)
Initialize tke value at vertical layer centers and interfaces from tracer (tke and tkeh)
Compute reciprocal of $ \Delta z $ (rdzt)
Compute reciprocal of pressure (tx1, tx2)
Compute minimum turbulent mixing length (rlmnz)
Compute background vertical diffusivities for scalars and momentum (xkzo and xkzmo)
set background diffusivities with xkzm_h & xkzm_m for gdx >= xkgdx and as a function of horizontal grid size for gdx < xkgdx
xkzm_hx = xkzm_h * (gdx / xkgdx)
xkzm_mx = xkzm_m * (gdx / xkgdx)
Some output variables and logical flags are initialized
Compute a function for green vegetation fraction and surface roughness. Entrainment rate in updraft is a function of vegetation fraction and surface roughness length
Compute $\theta$(theta), and $q_l$(qlx), $\theta_e$(thetae), $\theta_v$(thvx), $\theta_{l,v}$ (thlvx) including ice water
Compute an empirical cloud fraction based on Xu and Randall (1996) [198]
Compute buoyancy modified by clouds
Initialize diffusion coefficients to 0 and calculate the total radiative heating rate (dku, dkt, radx)
Compute stable/unstable PBL flag (pblflg) based on the total surface energy flux (false if the total surface energy flux is into the surface)

Calculate the PBL height

The calculation of the boundary layer height follows Troen and Mahrt (1986) [190] section 3. The approach is to find the level in the column where a modified bulk Richardson number exceeds a critical value.

Compute critical bulk Richardson number ( $Rb_{cr}$) (crb)
For the unstable PBL, crb is a constant (0.25)
For the stable boundary layer (SBL), $Rb_{cr}$ varies with the surface Rossby number, $R_{0}$, as given by Vickers and Mahrt (2004) [193]
\[ Rb_{cr}=0.16(10^{-7}R_{0})^{-0.18} \]

\[ R_{0}=\frac{U_{10}}{f_{0}z_{0}} \]
where $U_{10}$ is the wind speed at 10m above the ground surface, $f_0$ the Coriolis parameter, and $z_{0}$ the surface roughness length. To avoid too much variation, we restrict $Rb_{cr}$ to vary within the range of 0.15~0.35

Compute $\frac{\Delta t}{\Delta z}$ , $u_*$

Compute buoyancy $\frac{\partial \theta_v}{\partial z}$ (bf) and the wind shear squared (shr2)
Given the thermal's properties and the critical Richardson number, a loop is executed to find the first level above the surface (kpblx) where the modified Richardson number is greater than the critical Richardson number, using equation 10a from Troen and Mahrt (1996) [190] (also equation 8 from Hong and Pan (1996) [95]):
Once the level is found, some linear interpolation is performed to find the exact height of the boundary layer top (where $R_{i} > Rb_{cr}$) and the PBL height (hpbl and kpbl) and the PBL top index are saved.
Compute mean tke within pbl
Compute wind shear term as a sink term for updraft and downdraft velocity

Compute Monin-Obukhov similarity parameters

Calculate the Monin-Obukhov nondimensional stability paramter, commonly referred to as $\zeta$ using the following equation from Businger et al.(1971) [33] (eqn 28):
\[ \zeta = Ri_{sfc}\frac{F_m^2}{F_h} = \frac{z}{L} \]
where $F_m$ and $F_h$ are surface Monin-Obukhov stability functions calculated in sfc_diff.f and $L$ is the Obukhov length.

Calculate the nondimensional gradients of momentum and temperature ( $\phi_m$ (phim) and $\phi_h$(phih)) are calculated using eqns 5 and 6 from Hong and Pan (1996) [95] depending on the surface layer stability:

For the unstable and neutral conditions:
\[ \phi_m=(1-16\frac{0.1h}{L})^{-1/4} \phi_h=(1-16\frac{0.1h}{L})^{-1/2} \]
For the stable regime
\[ \phi_m=\phi_t=(1+5\frac{0.1h}{L}) \]

The $z/L$ (zol) is used as the stability criterion for the PBL.Currently, strong unstable (convective) PBL for $z/L < -0.02$ and weakly and moderately unstable PBL for $0>z/L>-0.02$

Compute the velocity scale $w_s$ (wscale) (eqn 22 of Han et al. 2019). It is represented by the value scaled at the top of the surface layer:
\[ w_s=(u_*^3+7\alpha\kappa w_*^3)^{1/3} \]
where $u_*$ (ustar) is the surface friction velocity, $\alpha$ is the ratio of the surface layer height to the PBL height (specified as sfcfrac =0.1), $\kappa =0.4$ is the von Karman constant, and $w_*$ is the convective velocity scale defined as eqn23 of Han et al.(2019):
\[ w_{*}=[(g/T)\overline{(w'\theta_v^{'})}_0h]^{1/3} \]

The counter-gradient terms for temperature and humidity are calculated.

Equation 4 of Hong and Pan (1996) [95] and are used to calculate the "scaled virtual temperature excess near the surface" (equation 9 in Hong and Pan (1996) [95]) for use in the mass-flux algorithm.

Determine whether stratocumulus layers exist and compute quantities

Starting at the PBL top and going downward, if the level is less than 2.5 km and $q_l\geq q_{lcr}$ then set kcld = k (find the cloud top index in the PBL. If no cloud water above the threshold is hound, scuflg is set to F.

Starting at the PBL top and going downward, if the level is less than the cloud top, find the level of the minimum radiative heating rate wihin the cloud. If the level of the minimum is the lowest model level or the minimum radiative heating rate is positive, then set scuflg to F.

Compute components for mass flux mixing by large thermals

If the PBL is convective, the updraft properties are initialized to be the same as the state variables.

Call mfpbltq(), which is an EDMF parameterization (Siebesma et al.(2007) [178]) to take into account nonlocal transport by large eddies. For details of the mfpbltq subroutine, step into its documentation ::mfpbltq

Call mfscuq(), which is a new mass-flux parameterization for stratocumulus-top-induced turbulence mixing. For details of the mfscuq subroutine, step into its documentation ::mfscuq
unify mass fluxes with Cu
taper off MF in high-wind conditions

Compute Prandtl number \f$P_r\f$ (prn) and exchange coefficient varying with height

Compute an asymtotic mixing length

Following Bougeault and Lacarrere(1989), the characteristic length scale ( $l_2$) (eqn 10 in Han et al.(2019) [80]) is given by:
\[ l_2=min(l_{up},l_{down}) \]
and dissipation length scale $l_d$ is given by:
\[ l_d=(l_{up}l_{down})^{1/2} \]
where $l_{up}$ and $l_{down}$ are the distances that a parcel having an initial TKE can travel upward and downward before being stopped by buoyancy effects.
Compute the surface layer length scale ( $l_1$) following Nakanishi (2001) [143] (eqn 9 of Han et al.(2019) [80])
If sfc_rlm=1, use zk for elm within surface layer

Compute eddy diffusivities

Compute TKE.

Compute a minimum TKE deduced from background diffusivity for momentum.

Compute buoyancy and shear productions of TKE

First predict tke due to tke production & dissipation(diss)
Compute updraft & downdraft properties for TKE
Compute tridiagonal matrix elements for turbulent kinetic energy
Call tridit() to solve tridiagonal problem for TKE
Recover the tendency of tke

Compute tridiagonal matrix elements for heat and moisture

Call tridin() to solve tridiagonal problem for heat and moisture
Recover the tendencies of heat and moisture

Add TKE dissipative heating to temperature tendency

Compute tridiagonal matrix elements for momentum

Call tridi2() to solve tridiagonal problem for momentum
Recover the tendencies of momentum

Save PBL height for diagnostic purpose

References mfpbltq_mod::mfpbltq(), mfscuq_mod::mfscuq(), tridi_mod::tridi2(), tridi_mod::tridin(), and tridi_mod::tridit().

Here is the call graph for this function: