 |
CCPP Scientific Documentation
v5.0.0
|
|
|
CCPP Scientific Documentation
v5.0.0
|
|
| subroutine hedmf::hedmf_run | ( | integer, intent(in) | im, |
| integer, intent(in) | km, | ||
| integer, intent(in) | ntrac, | ||
| integer, intent(in) | ntcw, | ||
| 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) | tau, | ||
| real(kind=kind_phys), dimension(im,km,ntrac), intent(inout) | rtg, | ||
| real(kind=kind_phys), dimension(im,km), intent(in) | u1, | ||
| real(kind=kind_phys), dimension(im,km), intent(in) | v1, | ||
| real(kind=kind_phys), dimension(im,km), intent(in) | t1, | ||
| real(kind=kind_phys), dimension(im,km,ntrac), intent(in) | q1, | ||
| real(kind=kind_phys), dimension(im,km), intent(in) | swh, | ||
| real(kind=kind_phys), dimension(im,km), intent(in) | hlw, | ||
| real(kind=kind_phys), dimension(im), intent(in) | xmu, | ||
| real(kind=kind_phys), dimension(im), 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(im,km+1), intent(in) | prsi, | ||
| real(kind=kind_phys), dimension(im,km), intent(in) | del, | ||
| real(kind=kind_phys), dimension(im,km), intent(in) | prsl, | ||
| real(kind=kind_phys), dimension(im,km), intent(in) | prslk, | ||
| real(kind=kind_phys), dimension(im,km+1), intent(in) | phii, | ||
| real(kind=kind_phys), dimension(im,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, | ||
| real(kind=kind_phys), dimension(im), intent(inout) | hgamt, | ||
| real(kind=kind_phys), dimension(im), intent(inout) | hgamq, | ||
| real(kind=kind_phys), dimension(im,km-1), intent(out) | dkt, | ||
| 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, | ||
| logical, intent(in) | lprnt, | ||
| integer, intent(in) | ipr, | ||
| real(kind=kind_phys), intent(in) | xkzminv, | ||
| real(kind=kind_phys), intent(in) | moninq_fac, | ||
| logical, intent(in) | lssav, | ||
| logical, intent(in) | ldiag3d, | ||
| logical, intent(in) | qdiag3d, | ||
| integer, intent(in) | ntoz, | ||
| real(kind=kind_phys), dimension(:,:), intent(inout) | du3dt_PBL, | ||
| real(kind=kind_phys), dimension(:,:), intent(inout) | dv3dt_PBL, | ||
| real(kind=kind_phys), dimension(:,:), intent(inout) | dt3dt_PBL, | ||
| real(kind=kind_phys), dimension(:,:), intent(inout) | dq3dt_PBL, | ||
| real(kind=kind_phys), dimension(:,:), intent(inout) | do3dt_PBL, | ||
| logical, intent(in) | flag_for_pbl_generic_tend, | ||
| character(len=*), intent(out) | errmsg, | ||
| integer, intent(out) | errflg | ||
| ) |
The calculation of the boundary layer height follows Troen and Mahrt (1986) [174] 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) [174] (also equation 8 from Hong and Pan (1996) [87]):
\[ 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)
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) [28] (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) [87] 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) [87]. 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) [174]
\[ w_s = (u_*^3 + 7\epsilon k w_*^3)^{1/3} \]
Next, the counter-gradient terms for temperature and humidity are calculated using equation 4 of Hong and Pan (1996) [87] and are used to calculate the "scaled virtual temperature excess near the surface" (equation 9 in Hong and Pan (1996) [87]) 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.
For an unstable PBL, the Prandtl number is calculated according to Hong and Pan (1996) [87], equation 10, whereas for a stable boundary layer, the Prandtl number is simply \(Pr = \frac{\phi_h}{\phi_m}\).
Below the PBL top, the diffusion coefficients ( \(K_m\) and \(K_h\)) are calculated according to equation 2 in Hong and Pan (1996) [87] 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.
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) [117] :
\[ 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) [87] as the previous documentation suggests. They follow equation 14 of Louis (1979) [117] 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) [87]
\[ \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, 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
If a stratocumulus layer has been identified in the PBL, the diffusion coefficients in the PBL are modified in the following way.
\[ 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.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) [111] ].
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.
Following Han et al. (2016) [78] , turbulence dissipation contributes to the tendency of temperature in the following way. First, turbulence dissipation is calculated by equation 17 of Han et al. (2016) [78] for the PBL and equation 16 for the surface layer.
Next, the temperature tendency is updated following equation 14.
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.
References physcons::con_cp, physcons::con_fvirt, physcons::con_g, physcons::con_hvap, physcons::con_rd, funcphys::fpvs(), mfpbl(), tridi2(), and tridin().
1.8.11