subroutine satmedmfvdifq::satmedmfvdifq_run	(	integer, intent(in)	ix,
		integer, intent(in)	im,
		integer, intent(in)	km,
		integer, intent(in)	ntrac,
		integer, intent(in)	ntcw,
		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(im,km), intent(inout)	dv,
		real(kind=kind_phys), dimension(im,km), intent(inout)	du,
		real(kind=kind_phys), dimension(im,km), intent(inout)	tdt,
		real(kind=kind_phys), dimension(im,km,ntrac), intent(inout)	rtg,
		real(kind=kind_phys), dimension(ix,km), intent(in)	u1,
		real(kind=kind_phys), dimension(ix,km), intent(in)	v1,
		real(kind=kind_phys), dimension(ix,km), intent(in)	t1,
		real(kind=kind_phys), dimension(ix,km,ntrac), intent(in)	q1,
		real(kind=kind_phys), dimension(ix,km), intent(in)	swh,
		real(kind=kind_phys), dimension(ix,km), intent(in)	hlw,
		real(kind=kind_phys), dimension(im), intent(in)	xmu,
		real(kind=kind_phys), dimension(im), intent(in)	garea,
		real(kind=kind_phys), dimension(ix), intent(in)	psk,
		real(kind=kind_phys), dimension(im), intent(in)	rbsoil,
		real(kind=kind_phys), dimension(im), intent(in)	zorl,
		real(kind=kind_phys), dimension(im), intent(in)	u10m,
		real(kind=kind_phys), dimension(im), intent(in)	v10m,
		real(kind=kind_phys), dimension(im), intent(in)	fm,
		real(kind=kind_phys), dimension(im), intent(in)	fh,
		real(kind=kind_phys), dimension(im), intent(in)	tsea,
		real(kind=kind_phys), dimension(im), intent(in)	heat,
		real(kind=kind_phys), dimension(im), intent(in)	evap,
		real(kind=kind_phys), dimension(im), intent(in)	stress,
		real(kind=kind_phys), dimension(im), intent(in)	spd1,
		integer, dimension(im), intent(out)	kpbl,
		real(kind=kind_phys), dimension(ix,km+1), intent(in)	prsi,
		real(kind=kind_phys), dimension(ix,km), intent(in)	del,
		real(kind=kind_phys), dimension(ix,km), intent(in)	prsl,
		real(kind=kind_phys), dimension(ix,km), intent(in)	prslk,
		real(kind=kind_phys), dimension(ix,km+1), intent(in)	phii,
		real(kind=kind_phys), dimension(ix,km), intent(in)	phil,
		real(kind=kind_phys), intent(in)	delt,
		logical, intent(in)	dspheat,
		real(kind=kind_phys), dimension(im), intent(out)	dusfc,
		real(kind=kind_phys), dimension(im), intent(out)	dvsfc,
		real(kind=kind_phys), dimension(im), intent(out)	dtsfc,
		real(kind=kind_phys), dimension(im), intent(out)	dqsfc,
		real(kind=kind_phys), dimension(im), intent(out)	hpbl,
		integer, dimension(im), 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,
		character(len=*), intent(out)	errmsg,
		integer, intent(out)	errflg
	)

Argument Table

local_name	standard_name	long_name	units	type	dimensions	kind	intent	optional
ix	horizontal_dimension	horizontal dimension	count	integer	()		in	False
im	horizontal_loop_extent	horizontal loop extent	count	integer	()		in	False
km	vertical_dimension	vertical layer dimension	count	integer	()		in	False
ntrac	number_of_vertical_diffusion_tracers	number of tracers to diffuse vertically	count	integer	()		in	False
ntcw	index_for_liquid_cloud_condensate	tracer index for cloud condensate (or liquid water)	index	integer	()		in	False
ntiw	index_for_ice_cloud_condensate_vertical_diffusion_tracer	tracer index for ice water in the vertically diffused tracer array	index	integer	()		in	False
ntke	index_for_turbulent_kinetic_energy_vertical_diffusion_tracer	index for turbulent kinetic energy in the vertically diffused tracer array	index	integer	()		in	False
grav	gravitational_acceleration	gravitational acceleration	m s-2	real	()	kind_phys	in	False
rd	gas_constant_dry_air	ideal gas constant for dry air	J kg-1 K-1	real	()	kind_phys	in	False
cp	specific_heat_of_dry_air_at_constant_pressure	specific heat of dry air at constant pressure	J kg-1 K-1	real	()	kind_phys	in	False
rv	gas_constant_water_vapor	ideal gas constant for water vapor	J kg-1 K-1	real	()	kind_phys	in	False
hvap	latent_heat_of_vaporization_of_water_at_0C	latent heat of evaporation/sublimation	J kg-1	real	()	kind_phys	in	False
hfus	latent_heat_of_fusion_of_water_at_0C	latent heat of fusion	J kg-1	real	()	kind_phys	in	False
fv	ratio_of_vapor_to_dry_air_gas_constants_minus_one	(rv/rd) - 1 (rv = ideal gas constant for water vapor)	none	real	()	kind_phys	in	False
eps	ratio_of_dry_air_to_water_vapor_gas_constants	rd/rv	none	real	()	kind_phys	in	False
epsm1	ratio_of_dry_air_to_water_vapor_gas_constants_minus_one	(rd/rv) - 1	none	real	()	kind_phys	in	False
dv	tendency_of_y_wind_due_to_model_physics	updated tendency of the y wind	m s-2	real	(horizontal_dimension, vertical_dimension)	kind_phys	inout	False
du	tendency_of_x_wind_due_to_model_physics	updated tendency of the x wind	m s-2	real	(horizontal_dimension, vertical_dimension)	kind_phys	inout	False
tdt	tendency_of_air_temperature_due_to_model_physics	updated tendency of the temperature	K s-1	real	(horizontal_dimension, vertical_dimension)	kind_phys	inout	False
rtg	tendency_of_vertically_diffused_tracer_concentration	updated tendency of the tracers due to vertical diffusion in PBL scheme	kg kg-1 s-1	real	(horizontal_dimension, vertical_dimension, number_of_vertical_diffusion_tracers)	kind_phys	inout	False
u1	x_wind	x component of layer wind	m s-1	real	(horizontal_dimension, vertical_dimension)	kind_phys	in	False
v1	y_wind	y component of layer wind	m s-1	real	(horizontal_dimension, vertical_dimension)	kind_phys	in	False
t1	air_temperature	layer mean air temperature	K	real	(horizontal_dimension, vertical_dimension)	kind_phys	in	False
q1	vertically_diffused_tracer_concentration	tracer concentration diffused by PBL scheme	kg kg-1	real	(horizontal_dimension, vertical_dimension, number_of_vertical_diffusion_tracers)	kind_phys	in	False
swh	tendency_of_air_temperature_due_to_shortwave_heating_on_radiation_time_step	total sky shortwave heating rate	K s-1	real	(horizontal_dimension, adjusted_vertical_layer_dimension_for_radiation)	kind_phys	in	False
hlw	tendency_of_air_temperature_due_to_longwave_heating_on_radiation_time_step	total sky longwave heating rate	K s-1	real	(horizontal_dimension, adjusted_vertical_layer_dimension_for_radiation)	kind_phys	in	False
xmu	zenith_angle_temporal_adjustment_factor_for_shortwave_fluxes	zenith angle temporal adjustment factor for shortwave	none	real	(horizontal_dimension)	kind_phys	in	False
garea	cell_area	area of the grid cell	m2	real	(horizontal_dimension)	kind_phys	in	False
psk	dimensionless_exner_function_at_lowest_model_interface	dimensionless Exner function at the surface interface	none	real	(horizontal_dimension)	kind_phys	in	False
rbsoil	bulk_richardson_number_at_lowest_model_level	bulk Richardson number at the surface	none	real	(horizontal_dimension)	kind_phys	in	False
zorl	surface_roughness_length	surface roughness length in cm	cm	real	(horizontal_dimension)	kind_phys	in	False
u10m	x_wind_at_10m	x component of wind at 10 m	m s-1	real	(horizontal_dimension)	kind_phys	in	False
v10m	y_wind_at_10m	y component of wind at 10 m	m s-1	real	(horizontal_dimension)	kind_phys	in	False
fm	Monin_Obukhov_similarity_function_for_momentum	Monin-Obukhov similarity function for momentum	none	real	(horizontal_dimension)	kind_phys	in	False
fh	Monin_Obukhov_similarity_function_for_heat	Monin-Obukhov similarity function for heat	none	real	(horizontal_dimension)	kind_phys	in	False
tsea	surface_skin_temperature	surface skin temperature	K	real	(horizontal_dimension)	kind_phys	in	False
heat	kinematic_surface_upward_sensible_heat_flux	kinematic surface upward sensible heat flux	K m s-1	real	(horizontal_dimension)	kind_phys	in	False
evap	kinematic_surface_upward_latent_heat_flux	kinematic surface upward latent heat flux	kg kg-1 m s-1	real	(horizontal_dimension)	kind_phys	in	False
stress	surface_wind_stress	surface wind stress	m2 s-2	real	(horizontal_dimension)	kind_phys	in	False
spd1	wind_speed_at_lowest_model_layer	wind speed at lowest model level	m s-1	real	(horizontal_dimension)	kind_phys	in	False
kpbl	vertical_index_at_top_of_atmosphere_boundary_layer	PBL top model level index	index	integer	(horizontal_dimension)		out	False
prsi	air_pressure_at_interface	air pressure at model layer interfaces	Pa	real	(horizontal_dimension, vertical_dimension_plus_one)	kind_phys	in	False
del	air_pressure_difference_between_midlayers	pres(k) - pres(k+1)	Pa	real	(horizontal_dimension, vertical_dimension)	kind_phys	in	False
prsl	air_pressure	mean layer pressure	Pa	real	(horizontal_dimension, vertical_dimension)	kind_phys	in	False
prslk	dimensionless_exner_function_at_model_layers	Exner function at layers	none	real	(horizontal_dimension, vertical_dimension)	kind_phys	in	False
phii	geopotential_at_interface	geopotential at model layer interfaces	m2 s-2	real	(horizontal_dimension, vertical_dimension_plus_one)	kind_phys	in	False
phil	geopotential	geopotential at model layer centers	m2 s-2	real	(horizontal_dimension, vertical_dimension)	kind_phys	in	False
delt	time_step_for_physics	time step for physics	s	real	()	kind_phys	in	False
dspheat	flag_TKE_dissipation_heating	flag for using TKE dissipation heating	flag	logical	()		in	False
dusfc	instantaneous_surface_x_momentum_flux	x momentum flux	Pa	real	(horizontal_dimension)	kind_phys	out	False
dvsfc	instantaneous_surface_y_momentum_flux	y momentum flux	Pa	real	(horizontal_dimension)	kind_phys	out	False
dtsfc	instantaneous_surface_upward_sensible_heat_flux	surface upward sensible heat flux	W m-2	real	(horizontal_dimension)	kind_phys	out	False
dqsfc	instantaneous_surface_upward_latent_heat_flux	surface upward latent heat flux	W m-2	real	(horizontal_dimension)	kind_phys	out	False
hpbl	atmosphere_boundary_layer_thickness	PBL thickness	m	real	(horizontal_dimension)	kind_phys	out	False
kinver	index_of_highest_temperature_inversion	index of highest temperature inversion	index	integer	(horizontal_dimension)		in	False
xkzm_m	atmosphere_momentum_diffusivity_background	background value of momentum diffusivity	m2 s-1	real	()	kind_phys	in	False
xkzm_h	atmosphere_heat_diffusivity_background	background value of heat diffusivity	m2 s-1	real	()	kind_phys	in	False
xkzm_s	diffusivity_background_sigma_level	sigma level threshold for background diffusivity	none	real	()	kind_phys	in	False
dspfac	tke_dissipative_heating_factor	tke dissipative heating factor	none	real	()	kind_phys	in	False
bl_upfr	updraft_fraction_in_boundary_layer_mass_flux_scheme	updraft fraction in boundary layer mass flux scheme	none	real	()	kind_phys	in	False
bl_dnfr	downdraft_fraction_in_boundary_layer_mass_flux_scheme	downdraft fraction in boundary layer mass flux scheme	none	real	()	kind_phys	in	False
errmsg	ccpp_error_message	error message for error handling in CCPP	none	character	()	len=*	out	False
errflg	ccpp_error_flag	error flag for error handling in CCPP	flag	integer	()		out	False

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) [73] .

1. The local turbulent mixing is represented by an eddy-diffusivity scheme which is a function of a prognostic TKE.
2. For the convective boundary layer, nonlocal transport by large eddies (mfpbltq.f), is represented using a mass flux approach (Siebesma et al.(2007) [157] ).
3. 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 as a function of horizontal grid size with xkzm_h & xkzm_m for gdx >= 25km and 0.01 for gdx=5m, i.e.,
  xkzm_hx = 0.01 + (xkzm_h - 0.01)/(xkgdx-5.) * (gdx-5.)
  xkzm_mx = 0.01 + (xkzm_h - 0.01)/(xkgdx-5.) * (gdx-5.)
• Some output variables and logical flags are initialized
• 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) [177]
• 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) [169] 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) [171]

  \[ 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) [169] (also equation 8 from Hong and Pan (1996) [85]):
• 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 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) [28] (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) [85] 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) [85] and are used to calculate the "scaled virtual temperature excess near the surface" (equation 9 in Hong and Pan (1996) [85]) 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) [157]) 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

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

The background vertical diffusivities in the inversion layers are limited

Compute an asymtotic mixing length

• Following Bougeault and Lacarrere(1989), the characteristic length scale ( \(l_2\)) (eqn 10 in Han et al.(2019) [73]) 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) [129] (eqn 9 of Han et al.(2019) [73])

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 funcphys::fpvs(), mfpbltq(), mfscuq(), tridi2(), tridin(), and tridit().

Here is the call graph for this function: