GFS Operational Physics Documentation  Revision: 81451
Hybrid Eddy-diffusivity Mass-flux Scheme

The Hybrid EDMF scheme is a first-order turbulent transport scheme used for subgrid-scale vertical turbulent mixing in the PBL and above. It blends the traditional first-order approach that has been used and improved over the last several years with a more recent scheme that uses a mass-flux approach to calculate the countergradient diffusion terms. More...

Detailed Description

The Hybrid EDMF scheme is a first-order turbulent transport scheme used for subgrid-scale vertical turbulent mixing in the PBL and above. It blends the traditional first-order approach that has been used and improved over the last several years with a more recent scheme that uses a mass-flux approach to calculate the countergradient diffusion terms.

The PBL scheme's main task is to calculate tendencies of temperature, moisture, and momentum due to vertical diffusion throughout the column (not just the PBL). The scheme is an amalgamation of decades of work, starting from the initial first-order PBL scheme of Troen and Mahrt (1986) [53], implemented according to Hong and Pan (1996) [26] and modified by Han and Pan (2011) [21] and Han et al. (2015) [22] to include top-down mixing due to stratocumulus layers from Lock et al. (2000) [38] and replacement of counter-gradient terms with a mass flux scheme according to Siebesma et al. (2007) [49] and Soares et al. (2004) [50]. Recently, heating due to TKE dissipation was also added according to Han et al. (2015) [22].

Calling Hierarchy Diagram

Hybrid_EDMF_Flowchart.png
Diagram depicting how the Hybrid EDMF PBL scheme is called from the GSM physics time loop

Intraphysics Communication

This space is reserved for a description of how this scheme uses information from other scheme types and/or how information calculated in this scheme is used in other scheme types.

Functions/Subroutines

subroutine tridi2 (l, n, cl, cm, cu, r1, r2, au, a1, a2)
 Routine to solve the tridiagonal system to calculate temperature and moisture at \( t + \Delta t \); part of two-part process to calculate time tendencies due to vertical diffusion. More...
 
subroutine tridin (l, n, nt, cl, cm, cu, r1, r2, au, a1, a2)
 Routine to solve the tridiagonal system to calculate u- and v-momentum at \( t + \Delta t \); part of two-part process to calculate time tendencies due to vertical diffusion. More...
 
subroutine mfpbl (im, ix, km, ntrac, delt, cnvflg, zl, zm, thvx, q1, t1, u1, v1, hpbl, kpbl, sflx, ustar, wstar, xmf, tcko, qcko, ucko, vcko)
 This subroutine is used for calculating the mass flux and updraft properties. More...
 
subroutine moninedmf (ix, im, km, ntrac, ntcw, dv, du, tau, rtg, u1, v1, t1, q1, swh, hlw, xmu, psk, rbsoil, zorl, u10m, v10m, fm, fh, tsea, qss, heat, evap, stress, spd1, kpbl, prsi, del, prsl, prslk, phii, phil, delt, dspheat, dusfc, dvsfc, dtsfc, dqsfc, hpbl, hgamt, hgamq, dkt, kinver, xkzm_m, xkzm_h, xkzm_s, lprnt, ipr)
 This subroutine contains all of logic for the Hybrid EDMF PBL scheme except for the calculation of the updraft properties and mass flux. More...
 

Function/Subroutine Documentation

subroutine mfpbl ( integer  im,
integer  ix,
integer  km,
integer  ntrac,
real(kind=kind_phys)  delt,
logical, dimension(im)  cnvflg,
real(kind=kind_phys), dimension(im,km)  zl,
real(kind=kind_phys), dimension(im,km+1)  zm,
real(kind=kind_phys), dimension(im,km)  thvx,
real(kind=kind_phys), dimension(ix,km,ntrac)  q1,
real(kind=kind_phys), dimension(ix,km)  t1,
real(kind=kind_phys), dimension(ix,km)  u1,
real(kind=kind_phys), dimension(ix,km)  v1,
real(kind=kind_phys), dimension(im)  hpbl,
integer, dimension(im)  kpbl,
real(kind=kind_phys), dimension(im)  sflx,
real(kind=kind_phys), dimension(im)  ustar,
real(kind=kind_phys), dimension(im)  wstar,
real(kind=kind_phys), dimension(im,km)  xmf,
real(kind=kind_phys), dimension(im,km)  tcko,
real(kind=kind_phys), dimension(im,km,ntrac)  qcko,
real(kind=kind_phys), dimension(im,km)  ucko,
real(kind=kind_phys), dimension(im,km)  vcko 
)

This subroutine is used for calculating the mass flux and updraft properties.

The mfpbl routines works as follows: if the PBL is convective, first, the ascending parcel entrainment rate is calculated as a function of height. Next, a surface parcel is initiated according to surface layer properties and the updraft buoyancy is calculated as a function of height. Next, using the buoyancy and entrainment values, the parcel vertical velocity is calculated using a well known steady-state budget equation. With the profile of updraft vertical velocity, the PBL height is recalculated as the height where the updraft vertical velocity returns to 0, and the entrainment profile is updated with the new PBL height. Finally, the mass flux profile is calculated using the updraft vertical velocity and assumed updraft fraction and the updraft properties are calculated using the updated entrainment profile, surface values, and environmental profiles.

Parameters
[in]imnumber of used points
[in]ixhorizontal dimension
[in]kmvertical layer dimension
[in]ntracnumber of tracers
[in]deltphysics time step
[in]cnvflgflag to denote a strongly unstable (convective) PBL
[in]zlheight of grid centers
[in]zmheight of grid interfaces
[in]thvxvirtual potential temperature at grid centers ( \( K \))
[in]q1layer mean tracer concentration (units?)
[in]t1layer mean temperature ( \( K \))
[in]u1u component of layer wind ( \( m s^{-1} \))
[in]v1v component of layer wind ( \( m s^{-1} \))
[in,out]hpblPBL top height (m)
[in,out]kpblPBL top index
[in]sflxtotal surface heat flux (units?)
[in]ustarsurface friction velocity
[in]wstarconvective velocity scale
[out]xmfupdraft mass flux
[in,out]tckoupdraft temperature ( \( K \))
[in,out]qckoupdraft tracer concentration (units?)
[in,out]uckoupdraft u component of horizontal momentum ( \( m s^{-1} \))
[in,out]vckoupdraft v component of horizontal momentum ( \( m s^{-1} \))

General Algorithm

  1. Determine an updraft parcel's entrainment rate, buoyancy, and vertical velocity.
  2. Recalculate the PBL height (previously calculated in moninedmf) and the parcel's entrainment rate.
  3. Calculate the mass flux profile and updraft properties.

Detailed Algorithm

Since the mfpbl subroutine is called regardless of whether the PBL is convective, a check of the convective PBL flag is performed and the subroutine returns back to moninedmf (with the output variables set to the initialized values) if the PBL is not convective.

Determine an updraft parcel's entrainment rate, buoyancy, and vertical velocity.

Calculate the entrainment rate according to equation 16 in Siebesma et al. (2007) [49] for all levels (xlamue) and a default entrainment rate (xlamax) for use above the PBL top.

Using equations 17 and 7 from Siebesma et al (2007) [49] along with \(u_*\), \(w_*\), and the previously diagnosed PBL height, the initial \(\theta_v\) of the updraft (and its surface buoyancy) is calculated.

From the second level to the middle of the vertical domain, the updraft virtual potential temperature is calculated using the entraining updraft equation as in equation 10 of Siebesma et al (2007) [49], discretized as

\[ \frac{\theta_{v,u}^k - \theta_{v,u}^{k-1}}{\Delta z}=-\epsilon^{k-1}\left[\frac{1}{2}\left(\theta_{v,u}^k + \theta_{v,u}^{k-1}\right)-\frac{1}{2}\left(\overline{\theta_{v}}^k + \overline{\theta_v}^{k-1}\right)\right] \]

where the superscript \(k\) denotes model level, and subscript \(u\) denotes an updraft property, and the overbar denotes the grid-scale mean value.

Rather than use the vertical velocity equation given as equation 15 in Siebesma et al (2007) [49] (which parameterizes the pressure term in terms of the updraft vertical velocity itself), this scheme uses the more widely used form of the steady state vertical velocity equation given as equation 6 in Soares et al. (2004) [50] discretized as

\[ \frac{w_{u,k}^2 - w_{u,k-1}^2}{\Delta z} = -2b_1\frac{1}{2}\left(\epsilon_k + \epsilon_{k-1}\right)\frac{1}{2}\left(w_{u,k}^2 + w_{u,k-1}^2\right) + 2b_2B \]

The constants used in the scheme are labeled \(bb1 = 2b_1\) and \(bb2 = 2b_2\) and are tuned to be equal to 1.8 and 3.5, respectively, close to the values proposed by Soares et al. (2004) [50] .

Recalculate the PBL height and the parcel's entrainment rate.

Find the level where the updraft vertical velocity is less than zero and linearly interpolate to find the height where it would be exactly zero. Set the PBL height to this determined height.

Recalculate the entrainment rate as before except use the updated value of the PBL height.

Calculate the mass flux profile and updraft properties.

Calculate the mass flux:

\[ M = a_uw_u \]

where \(a_u\) is the tunable parameter that represents the fractional area of updrafts (currently set to 0.08). Limit the computed mass flux to be less than \(\frac{\Delta z}{\Delta t}\). This is different than what is done in Siebesma et al. (2007) [49] where the mass flux is the product of a tunable constant and the diagnosed standard deviation of \(w\).

The updraft properties are calculated according to the entraining updraft equation

\[ \frac{\partial \phi}{\partial z}=-\epsilon\left(\phi_u - \overline{\phi}\right) \]

where \(\phi\) is \(T\) or \(q\). The equation is discretized according to

\[ \frac{\phi_{u,k} - \phi_{u,k-1}}{\Delta z}=-\epsilon_{k-1}\left[\frac{1}{2}\left(\phi_{u,k} + \phi_{u,k-1}\right)-\frac{1}{2}\left(\overline{\phi}_k + \overline{\phi}_{k-1}\right)\right] \]

The exception is for the horizontal momentum components, which have been modified to account for the updraft-induced pressure gradient force, and use the following equation, following Han and Pan (2006) [20]

\[ \frac{\partial v}{\partial z} = -\epsilon\left(v_u - \overline{v}\right)+d_1\frac{\partial \overline{v}}{\partial z} \]

where \(d_1=0.55\) is a tunable constant.

Definition at line 41 of file mfpbl.f.

References physcons::con_cp, and physcons::con_g.

Referenced by moninedmf().

Here is the caller graph for this function:

subroutine moninedmf ( integer  ix,
integer  im,
integer  km,
integer  ntrac,
integer  ntcw,
real(kind=kind_phys), dimension(im,km)  dv,
real(kind=kind_phys), dimension(im,km)  du,
real(kind=kind_phys), dimension(im,km)  tau,
real(kind=kind_phys), dimension(im,km,ntrac)  rtg,
real(kind=kind_phys), dimension(ix,km)  u1,
real(kind=kind_phys), dimension(ix,km)  v1,
real(kind=kind_phys), dimension(ix,km)  t1,
real(kind=kind_phys), dimension(ix,km,ntrac)  q1,
real(kind=kind_phys), dimension(ix,km)  swh,
real(kind=kind_phys), dimension(ix,km)  hlw,
real(kind=kind_phys), dimension(im)  xmu,
real(kind=kind_phys), dimension(im)  psk,
real(kind=kind_phys), dimension(im)  rbsoil,
real(kind=kind_phys), dimension(im)  zorl,
real(kind=kind_phys), dimension(im)  u10m,
real(kind=kind_phys), dimension(im)  v10m,
real(kind=kind_phys), dimension(im)  fm,
real(kind=kind_phys), dimension(im)  fh,
real(kind=kind_phys), dimension(im)  tsea,
real(kind=kind_phys), dimension(im)  qss,
real(kind=kind_phys), dimension(im)  heat,
real(kind=kind_phys), dimension(im)  evap,
real(kind=kind_phys), dimension(im)  stress,
real(kind=kind_phys), dimension(im)  spd1,
integer, dimension(im)  kpbl,
real(kind=kind_phys), dimension(ix,km+1)  prsi,
real(kind=kind_phys), dimension(ix,km)  del,
real(kind=kind_phys), dimension(ix,km)  prsl,
real(kind=kind_phys), dimension(ix,km)  prslk,
real(kind=kind_phys), dimension(ix,km+1)  phii,
real(kind=kind_phys), dimension(ix,km)  phil,
real(kind=kind_phys)  delt,
logical  dspheat,
real(kind=kind_phys), dimension(im)  dusfc,
real(kind=kind_phys), dimension(im)  dvsfc,
real(kind=kind_phys), dimension(im)  dtsfc,
real(kind=kind_phys), dimension(im)  dqsfc,
real(kind=kind_phys), dimension(im)  hpbl,
real(kind=kind_phys), dimension(im)  hgamt,
real(kind=kind_phys), dimension(im)  hgamq,
real(kind=kind_phys), dimension(im,km-1)  dkt,
integer, dimension(im)  kinver,
real(kind=kind_phys)  xkzm_m,
real(kind=kind_phys)  xkzm_h,
real(kind=kind_phys)  xkzm_s,
logical  lprnt,
integer  ipr 
)

This subroutine contains all of logic for the Hybrid EDMF PBL scheme except for the calculation of the updraft properties and mass flux.

The scheme works on a basic level by calculating background diffusion coefficients and updating them according to which processes are occurring in the column. The most important difference in diffusion coefficients occurs between those levels in the PBL and those above the PBL, so the PBL height calculation is of utmost importance. An initial estimate is calculated in a "predictor" step in order to calculate Monin-Obukhov similarity values and a corrector step recalculates the PBL height based on updated surface thermal characteristics. Using the PBL height and the similarity parameters, the diffusion coefficients are updated below the PBL top based on Hong and Pan (1996) [26] (including counter-gradient terms). Diffusion coefficients in the free troposphere (above the PBL top) are calculated according to Louis (1979) [40] with updated Richardson number-dependent functions. If it is diagnosed that PBL top-down mixing is occurring according to Lock et al. (2000) [38] , then then diffusion coefficients are updated accordingly. Finally, for convective boundary layers (defined as when the Obukhov length exceeds a threshold), the counter-gradient terms are replaced using the mass flux scheme of Siebesma et al. (2007) [49] . In order to return time tendencies, a fully implicit solution is found using tridiagonal matrices, and time tendencies are "backed out." Before returning, the time tendency of temperature is updated to reflect heating due to TKE dissipation following Han et al. (2015) [22] .

Parameters
[in]ixhorizontal dimension
[in]imnumber of used points
[in]kmvertical layer dimension
[in]ntracnumber of tracers
[in]ntcwcloud condensate index in the tracer array
[in,out]dvv-momentum tendency ( \( m s^{-2} \))
[in,out]duu-momentum tendency ( \( m s^{-2} \))
[in,out]tautemperature tendency ( \( K s^{-1} \))
[in,out]rtgmoisture tendency ( \( kg kg^{-1} s^{-1} \))
[in]u1u component of layer wind ( \( m s^{-1} \))
[in]v1v component of layer wind ( \( m s^{-1} \))
[in]t1layer mean temperature ( \( K \))
[in]q1layer mean tracer concentration (units?)
[in]swhtotal sky shortwave heating rate ( \( K s^-1 \))
[in]hlwtotal sky longwave heating rate ( \( K s^-1 \))
[in]xmutime step zenith angle adjust factor for shortwave
[in]pskExner function at surface interface?
[in]rbsoilsurface bulk Richardson number
[in]zorlsurface roughness (units?)
[in]u10m10-m u wind ( \( m s^{-1} \))
[in]v10m10-m v wind ( \( m s^{-1} \))
[in]fmfm parameter from PBL scheme
[in]fhfh parameter from PBL scheme
[in]tseaground surface temperature (K)
[in]qsssurface saturation humidity (units?)
[in]heatsurface sensible heat flux (units?)
[in]evapevaporation from latent heat flux (units?)
[in]stresssurface wind stress? ( \( cm*v^2\) in sfc_diff subroutine) (units?)
[in]spd1surface wind speed? (units?)
[out]kpblPBL top index
[in]prsipressure at layer interfaces (units?)
[in]delpressure difference between level k and k+1 (units?)
[in]prslmean layer pressure (units?)
[in]prslkExner function at layer
[in]phiiinterface geopotential height (units?)
[in]phillayer geopotential height (units?)
[in]deltphysics time step (s)
[in]dspheatflag for TKE dissipative heating
[out]dusfcsurface u-momentum tendency (units?)
[out]dvsfcsurface v-momentum tendency (units?)
[out]dtsfcsurface temperature tendency (units?)
[out]dqsfcsurface moisture tendency (units?)
[out]hpblPBL top height (m)
[out]hgamtcounter gradient mixing term for temperature (units?)
[out]hgamqcounter gradient mixing term for moisture (units?)
[out]dktdiffusion coefficient for temperature (units?)
[in]kinverindex location of temperature inversion
[in]xkzm_mbackground vertical diffusion coefficient for momentum (units?)
[in]xkzm_hbackground vertical diffusion coefficeint for heat, moisture (units?)
[in]xkzm_ssigma threshold for background momentum diffusion (units?)
[in]lprntflag to print some output
[in]iprindex of point to print

General Algorithm

  1. Compute preliminary variables from input arguments.
  2. Calculate the first estimate of the PBL height ("Predictor step").
  3. Calculate Monin-Obukhov similarity parameters.
  4. Update thermal properties of surface parcel and recompute PBL height ("Corrector step").
  5. Determine whether stratocumulus layers exist and compute quantities needed for enhanced diffusion.
  6. Calculate the inverse Prandtl number.
  7. Compute diffusion coefficients below the PBL top.
  8. Compute diffusion coefficients above the PBL top.
  9. If the PBL is convective, call the mass flux scheme to replace the countergradient terms.
  10. Compute enhanced diffusion coefficients related to stratocumulus-topped PBLs.
  11. Solve for the temperature and moisture tendencies due to vertical mixing.
  12. Calculate heating due to TKE dissipation and add to the tendency for temperature.
  13. Solve for the horizontal momentum tendencies and add them to output tendency terms.

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 reciprocal of \( \Delta z \) (rdzt)
  • Compute reciprocal of pressure (tx1, tx2)
  • Compute background vertical diffusivities for scalars and momentum (xkzo and xkzmo)
  • The background scalar vertical diffusivity is limited to be less than or equal to xkzminv
  • Some output variables and logical flags are initialized
  • Compute \(\theta\) (theta), \(q_l\) (qlx), \(q_t\) (qtx), \(\theta_e\) (thetae), \(\theta_v\) (thvx), \(\theta_{l,v}\) (thlvx)
  • Initialize diffusion coefficients to 0 and calculate the total radiative heating rate (dku, dkt, radx)
  • Set lcld to first index above 2.5km
  • Compute \(\frac{\partial \theta_v}{\partial z}\) (bf) and the wind shear squared (shr2)
  • Calculate \(\frac{g}{\theta}\) (govrth), \(\beta = \frac{\Delta t}{\Delta z}\) (beta), \(u_*\) (ustar), total surface flux (sflux), and set pblflag to false if the total surface energy flux is into the surface

Calculate the first estimate of the PBL height (``Predictor step")

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

The temperature of the thermal is of primary importance. For the initial estimate of the PBL height, the thermal is assumed to have one of two temperatures. If the boundary layer is stable, the thermal is assumed to have a temperature equal to the surface virtual temperature. Otherwise, the thermal is assumed to have the same virtual potential temperature as the lowest model level. For the stable case, the critical bulk Richardson number becomes a function of the wind speed and roughness length, otherwise it is set to a tunable constant.

Given the thermal's properties and the critical Richardson number, a loop is executed to find the first level above the surface where the modified Richardson number is greater than the critical Richardson number, using equation 10a from Troen and Mahrt (1986) [53] (also equation 8 from Hong and Pan (1996) [26]):

\[ h = Ri\frac{T_0\left|\vec{v}(h)\right|^2}{g\left(\theta_v(h) - \theta_s\right)} \]

where \(h\) is the PBL height, \(Ri\) is the Richardson number, \(T_0\) is the virtual potential temperature near the surface, \(\left|\vec{v}\right|\) is the wind speed, and \(\theta_s\) is for the thermal. Rearranging this equation to calculate the modified Richardson number at each level, k, for comparison with the critical value yields:

\[ Ri_k = gz(k)\frac{\left(\theta_v(k) - \theta_s\right)}{\theta_v(1)*\vec{v}(k)} \]

Once the level is found, some linear interpolation is performed to find the exact height of the boundary layer top (where \(Ri = Ri_{cr}\)) and the PBL height and the PBL top index are saved (hpblx and kpblx, respectively)

Calculate Monin-Obukhov similarity parameters

Using the initial guess for the PBL height, Monin-Obukhov similarity parameters are calculated. They are needed to refine the PBL height calculation and for calculating diffusion coefficients.

First, calculate the Monin-Obukhov nondimensional stability parameter, commonly referred to as \(\zeta\) using the following equation from Businger et al. (1971) [6] (equation 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. Then, the nondimensional gradients of momentum and temperature (phim and phih) are calculated using equations 5 and 6 from Hong and Pan (1996) [26] depending on the surface layer stability. Then, the velocity scale valid for the surface layer ( \(w_s\), wscale) is calculated using equation 3 from Hong and Pan (1996) [26]. For the neutral and unstable PBL above the surface layer, the convective velocity scale, \(w_*\), is calculated according to:

\[ w_* = \left(\frac{g}{\theta_0}h\overline{w'\theta_0'}\right)^{1/3} \]

and the mixed layer velocity scale is then calculated with equation 6 from Troen and Mahrt (1986) [53]

\[ w_s = (u_*^3 + 7\epsilon k w_*^3)^{1/3} \]

Update thermal properties of surface parcel and recompute PBL height ("Corrector step").

Next, the counter-gradient terms for temperature and humidity are calculated using equation 4 of Hong and Pan (1996) [26] and are used to calculate the "scaled virtual temperature excess near the surface" (equation 9 in Hong and Pan (1996) [26]) so that the properties of the thermal are updated to recalculate the PBL height.

The PBL height calculation follows the same procedure as the predictor step, except that it uses an updated virtual potential temperature for the thermal.

Determine whether stratocumulus layers exist and compute quantities needed for enhanced diffusion

  • Starting at the PBL top and going downward, if the level is less than 2.5 km and \(q_l>q_{l,cr}\) then set kcld = k (find the cloud top index in the PBL). If no cloud water above the threshold is found, 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 within 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.
  • Starting at the PBL top and going downward, count the number of levels below the minimum radiative heating rate level that have cloud water above the threshold. If there are none, then set the scuflg to F.
  • Find the height of the interface where the minimum in radiative heating rate is located. If this height is less than the second model interface height, then set the scuflg to F.
  • Calculate the hypothetical \(\theta_v\) at the minimum radiative heating level that a parcel would reach due to radiative cooling after a typical cloud turnover time spent at that level.
  • Determine the distance that a parcel would sink downwards starting from the level of minimum radiative heating rate by comparing the hypothetical minimum \(\theta_v\) calculated above with the environmental \(\theta_v\).
  • Calculate the cloud thickness, where the cloud top is the in-cloud minimum radiative heating level and the bottom is determined previously.
  • Find the largest between the cloud thickness and the distance of a sinking parcel, then determine the smallest of that number and the height of the minimum in radiative heating rate. Set this number to \(zd\). Using \(zd\), calculate the characteristic velocity scale of cloud-top radiative cooling-driven turbulence.

Calculate the inverse Prandtl number

For an unstable PBL, the Prandtl number is calculated according to Hong and Pan (1996) [26], equation 10, whereas for a stable boundary layer, the Prandtl number is simply \(Pr = \frac{\phi_h}{\phi_m}\).

Compute diffusion coefficients below the PBL top

Below the PBL top, the diffusion coefficients ( \(K_m\) and \(K_h\)) are calculated according to equation 2 in Hong and Pan (1996) [26] where a different value for \(w_s\) (PBL vertical velocity scale) is used depending on the PBL stability. \(K_h\) is calculated from \(K_m\) using the Prandtl number. The calculated diffusion coefficients are checked so that they are bounded by maximum values and the local background diffusion coefficients.

Compute diffusion coefficients above the PBL top

Diffusion coefficients above the PBL top are computed as a function of local stability (gradient Richardson number), shear, and a length scale from Louis (1979) [40] :

\[ K_{m,h}=l^2f_{m,h}(Ri_g)\left|\frac{\partial U}{\partial z}\right| \]

The functions used ( \(f_{m,h}\)) depend on the local stability. First, the gradient Richardson number is calculated as

\[ Ri_g=\frac{\frac{g}{T}\frac{\partial \theta_v}{\partial z}}{\frac{\partial U}{\partial z}^2} \]

where \(U\) is the horizontal wind. For the unstable case ( \(Ri_g < 0\)), the Richardson number-dependent functions are given by

\[ f_h(Ri_g) = 1 + \frac{8\left|Ri_g\right|}{1 + 1.286\sqrt{\left|Ri_g\right|}}\\ \]

\[ f_m(Ri_g) = 1 + \frac{8\left|Ri_g\right|}{1 + 1.746\sqrt{\left|Ri_g\right|}}\\ \]

For the stable case, the following formulas are used

\[ f_h(Ri_g) = \frac{1}{\left(1 + 5Ri_g\right)^2}\\ \]

\[ Pr = \frac{K_h}{K_m} = 1 + 2.1Ri_g \]

The source for the formulas used for the Richardson number-dependent functions is unclear. They are different than those used in Hong and Pan (1996) [26] as the previous documentation suggests. They follow equation 14 of Louis (1979) [40] for the unstable case, but it is unclear where the values of the coefficients \(b\) and \(c\) from that equation used in this scheme originate. Finally, the length scale, \(l\) is calculated according to the following formula from Hong and Pan (1996) [26]

\[ \frac{1}{l} = \frac{1}{kz} + \frac{1}{l_0}\\ \]

\[ or\\ \]

\[ l=\frac{l_0kz}{l_0+kz} \]

where \(l_0\) is currently 30 m for stable conditions and 150 m for unstable. Finally, the diffusion coefficients are kept in a range bounded by the background diffusion and the maximum allowable values.

If the PBL is convective, call the mass flux scheme to replace the countergradient terms.

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

For details of the mfpbl subroutine, step into its documentation mfpbl

Compute enhanced diffusion coefficients related to stratocumulus-topped PBLs

If a stratocumulus layer has been identified in the PBL, the diffusion coefficients in the PBL are modified in the following way.

  1. First, the criteria for CTEI is checked, using the threshold from equation 13 of Macvean and Mason (1990) [41]. If the criteria is met, the cloud top diffusion is increased:

    \[ K_h^{Sc} = -c\frac{\Delta F_R}{\rho c_p}\frac{1}{\frac{\partial \theta_v}{\partial z}} \]

    where the constant \(c\) is set to 0.2 if the CTEI criterion is not met and 1.0 if it is.
  2. Calculate the diffusion coefficients due to stratocumulus mixing according to equation 5 in Lock et al. (2000) [38] for every level below the stratocumulus top using the characteristic stratocumulus velocity scale previously calculated. The diffusion coefficient for momentum is calculated assuming a constant inverse Prandtl number of 0.75.

After \(K_h^{Sc}\) has been determined from the surface to the top of the stratocumulus layer, it is added to the value for the diffusion coefficient calculated previously using surface-based mixing [see equation 6 of Lock et al. (2000) [38] ].

Solve for the temperature and moisture tendencies due to vertical mixing.

The tendencies of heat, moisture, and momentum due to vertical diffusion are calculated using a two-part process. First, a solution is obtained using an implicit time-stepping scheme, then the time tendency terms are "backed out". The tridiagonal matrix elements for the implicit solution for temperature and moisture are prepared in this section, with differing algorithms depending on whether the PBL was convective (substituting the mass flux term for counter-gradient term), unstable but not convective (using the computed counter-gradient terms), or stable (no counter-gradient terms).

The tridiagonal system is solved by calling the internal tridin subroutine.

After returning with the solution, the tendencies for temperature and moisture are recovered.

Calculate heating due to TKE dissipation and add to the tendency for temperature

Following Han et al. (2015) [22] , turbulence dissipation contributes to the tendency of temperature in the following way. First, turbulence dissipation is calculated by equation 17 of Han et al. (2015) [22] for the PBL and equation 16 for the surface layer.

Next, the temperature tendency is updated following equation 14.

Solve for the horizontal momentum tendencies and add them to the output tendency terms

As with the temperature and moisture tendencies, the horizontal momentum tendencies are calculated by solving tridiagonal matrices after the matrices are prepared in this section.

Finally, the tendencies are recovered from the tridiagonal solutions.

Definition at line 96 of file moninedmf.f.

References physcons::con_cp, physcons::con_g, physcons::con_hvap, physcons::con_rd, mfpbl(), tridi2(), and tridin().

Referenced by gbphys().

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine tridi2 ( integer  l,
integer  n,
real(kind=kind_phys), dimension(l,2:n)  cl,
real(kind=kind_phys), dimension(l,n)  cm,
real(kind=kind_phys), dimension(l,n-1)  cu,
real(kind=kind_phys), dimension(l,n)  r1,
real(kind=kind_phys), dimension(l,n)  r2,
real(kind=kind_phys), dimension(l,n-1)  au,
real(kind=kind_phys), dimension(l,n)  a1,
real(kind=kind_phys), dimension(l,n)  a2 
)

Routine to solve the tridiagonal system to calculate temperature and moisture at \( t + \Delta t \); part of two-part process to calculate time tendencies due to vertical diffusion.

Origin of subroutine unknown.

Definition at line 1196 of file moninedmf.f.

Referenced by moninedmf().

Here is the caller graph for this function:

subroutine tridin ( integer  l,
integer  n,
integer  nt,
real(kind=kind_phys), dimension(l,2:n)  cl,
real(kind=kind_phys), dimension(l,n)  cm,
real(kind=kind_phys), dimension(l,n-1)  cu,
real(kind=kind_phys), dimension(l,n)  r1,
real(kind=kind_phys), dimension(l,n*nt)  r2,
real(kind=kind_phys), dimension(l,n-1)  au,
real(kind=kind_phys), dimension(l,n)  a1,
real(kind=kind_phys), dimension(l,n*nt)  a2 
)

Routine to solve the tridiagonal system to calculate u- and v-momentum at \( t + \Delta t \); part of two-part process to calculate time tendencies due to vertical diffusion.

Origin of subroutine unknown.

Definition at line 1239 of file moninedmf.f.

Referenced by moninedmf().

Here is the caller graph for this function: