CCPP Scientific Documentation
v5.0.0
subroutine satmedmfvdifq::satmedmfvdifq_run ( 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(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)  garea,
integer, dimension(im), intent(in)  islimsk,
real(kind=kind_phys), dimension(im), intent(in)  snwdph_lnd,
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,
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,
integer, intent(in)  ntoz,
real(kind=kind_phys), dimension(:,:), intent(inout)  du3dt,
real(kind=kind_phys), dimension(:,:), intent(inout)  dv3dt,
real(kind=kind_phys), dimension(:,:), intent(inout)  dt3dt,
real(kind=kind_phys), dimension(:,:), intent(inout)  dq3dt,
real(kind=kind_phys), dimension(:,:), intent(inout)  do3dt,
logical, intent(in)  gen_tend,
logical, intent(in)  ldiag3d,
logical, intent(in)  qdiag3d,
character(len=*), intent(out)  errmsg,
integer, intent(out)  errflg 
)

Argument Table

satmedmfvdifq_run argument table
local_namestandard_namelong_nameunitstypedimensionskindintentoptional
imhorizontal_loop_extenthorizontal loop extentcountinteger()inFalse
kmvertical_dimensionvertical layer dimensioncountinteger()inFalse
ntracnumber_of_vertical_diffusion_tracersnumber of tracers to diffuse verticallycountinteger()inFalse
ntcwindex_for_liquid_cloud_condensatetracer index for cloud condensate (or liquid water)indexinteger()inFalse
ntiwindex_for_ice_cloud_condensate_vertical_diffusion_tracertracer index for ice water in the vertically diffused tracer arrayindexinteger()inFalse
ntkeindex_for_turbulent_kinetic_energy_vertical_diffusion_tracerindex for turbulent kinetic energy in the vertically diffused tracer arrayindexinteger()inFalse
gravgravitational_accelerationgravitational accelerationm s-2real()kind_physinFalse
rdgas_constant_dry_airideal gas constant for dry airJ kg-1 K-1real()kind_physinFalse
cpspecific_heat_of_dry_air_at_constant_pressurespecific heat of dry air at constant pressureJ kg-1 K-1real()kind_physinFalse
rvgas_constant_water_vaporideal gas constant for water vaporJ kg-1 K-1real()kind_physinFalse
hvaplatent_heat_of_vaporization_of_water_at_0Clatent heat of evaporation/sublimationJ kg-1real()kind_physinFalse
hfuslatent_heat_of_fusion_of_water_at_0Clatent heat of fusionJ kg-1real()kind_physinFalse
fvratio_of_vapor_to_dry_air_gas_constants_minus_one(rv/rd) - 1 (rv = ideal gas constant for water vapor)nonereal()kind_physinFalse
epsratio_of_dry_air_to_water_vapor_gas_constantsrd/rvnonereal()kind_physinFalse
epsm1ratio_of_dry_air_to_water_vapor_gas_constants_minus_one(rd/rv) - 1nonereal()kind_physinFalse
dvtendency_of_y_wind_due_to_model_physicsupdated tendency of the y windm s-2real(horizontal_dimension, vertical_dimension)kind_physinoutFalse
dutendency_of_x_wind_due_to_model_physicsupdated tendency of the x windm s-2real(horizontal_dimension, vertical_dimension)kind_physinoutFalse
tdttendency_of_air_temperature_due_to_model_physicsupdated tendency of the temperatureK s-1real(horizontal_dimension, vertical_dimension)kind_physinoutFalse
rtgtendency_of_vertically_diffused_tracer_concentrationupdated tendency of the tracers due to vertical diffusion in PBL schemekg kg-1 s-1real(horizontal_dimension, vertical_dimension, number_of_vertical_diffusion_tracers)kind_physinoutFalse
u1x_windx component of layer windm s-1real(horizontal_dimension, vertical_dimension)kind_physinFalse
v1y_windy component of layer windm s-1real(horizontal_dimension, vertical_dimension)kind_physinFalse
t1air_temperaturelayer mean air temperatureKreal(horizontal_dimension, vertical_dimension)kind_physinFalse
q1vertically_diffused_tracer_concentrationtracer concentration diffused by PBL schemekg kg-1real(horizontal_dimension, vertical_dimension, number_of_vertical_diffusion_tracers)kind_physinFalse
swhtendency_of_air_temperature_due_to_shortwave_heating_on_radiation_time_steptotal sky shortwave heating rateK s-1real(horizontal_dimension, vertical_dimension)kind_physinFalse
hlwtendency_of_air_temperature_due_to_longwave_heating_on_radiation_time_steptotal sky longwave heating rateK s-1real(horizontal_dimension, vertical_dimension)kind_physinFalse
xmuzenith_angle_temporal_adjustment_factor_for_shortwave_fluxeszenith angle temporal adjustment factor for shortwavenonereal(horizontal_dimension)kind_physinFalse
gareacell_areaarea of the grid cellm2real(horizontal_dimension)kind_physinFalse
islimsksea_land_ice_masksea/land/ice mask (=0/1/2)flaginteger(horizontal_dimension)inFalse
snwdph_lndsurface_snow_thickness_water_equivalent_over_landwater equivalent snow depth over landmmreal(horizontal_dimension)kind_physinFalse
pskdimensionless_exner_function_at_lowest_model_interfacedimensionless Exner function at the surface interfacenonereal(horizontal_dimension)kind_physinFalse
rbsoilbulk_richardson_number_at_lowest_model_levelbulk Richardson number at the surfacenonereal(horizontal_dimension)kind_physinFalse
zorlsurface_roughness_lengthsurface roughness length in cmcmreal(horizontal_dimension)kind_physinFalse
u10mx_wind_at_10mx component of wind at 10 mm s-1real(horizontal_dimension)kind_physinFalse
v10my_wind_at_10my component of wind at 10 mm s-1real(horizontal_dimension)kind_physinFalse
fmMonin_Obukhov_similarity_function_for_momentumMonin-Obukhov similarity function for momentumnonereal(horizontal_dimension)kind_physinFalse
fhMonin_Obukhov_similarity_function_for_heatMonin-Obukhov similarity function for heatnonereal(horizontal_dimension)kind_physinFalse
tseasurface_skin_temperaturesurface skin temperatureKreal(horizontal_dimension)kind_physinFalse
heatkinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughnesskinematic surface upward sensible heat fluxK m s-1real(horizontal_dimension)kind_physinFalse
evapkinematic_surface_upward_latent_heat_flux_reduced_by_surface_roughnesskinematic surface upward latent heat fluxkg kg-1 m s-1real(horizontal_dimension)kind_physinFalse
stresssurface_wind_stresssurface wind stressm2 s-2real(horizontal_dimension)kind_physinFalse
spd1wind_speed_at_lowest_model_layerwind speed at lowest model levelm s-1real(horizontal_dimension)kind_physinFalse
kpblvertical_index_at_top_of_atmosphere_boundary_layerPBL top model level indexindexinteger(horizontal_dimension)outFalse
prsiair_pressure_at_interfaceair pressure at model layer interfacesPareal(horizontal_dimension, vertical_dimension_plus_one)kind_physinFalse
delair_pressure_difference_between_midlayerspres(k) - pres(k+1)Pareal(horizontal_dimension, vertical_dimension)kind_physinFalse
prslair_pressuremean layer pressurePareal(horizontal_dimension, vertical_dimension)kind_physinFalse
prslkdimensionless_exner_function_at_model_layersExner function at layersnonereal(horizontal_dimension, vertical_dimension)kind_physinFalse
phiigeopotential_at_interfacegeopotential at model layer interfacesm2 s-2real(horizontal_dimension, vertical_dimension_plus_one)kind_physinFalse
philgeopotentialgeopotential at model layer centersm2 s-2real(horizontal_dimension, vertical_dimension)kind_physinFalse
delttime_step_for_physicstime step for physicssreal()kind_physinFalse
dspheatflag_TKE_dissipation_heatingflag for using TKE dissipation heatingflaglogical()inFalse
dusfcinstantaneous_surface_x_momentum_fluxx momentum fluxPareal(horizontal_dimension)kind_physoutFalse
dvsfcinstantaneous_surface_y_momentum_fluxy momentum fluxPareal(horizontal_dimension)kind_physoutFalse
dtsfcinstantaneous_surface_upward_sensible_heat_fluxsurface upward sensible heat fluxW m-2real(horizontal_dimension)kind_physoutFalse
dqsfcinstantaneous_surface_upward_latent_heat_fluxsurface upward latent heat fluxW m-2real(horizontal_dimension)kind_physoutFalse
hpblatmosphere_boundary_layer_thicknessPBL thicknessmreal(horizontal_dimension)kind_physoutFalse
kinverindex_of_highest_temperature_inversionindex of highest temperature inversionindexinteger(horizontal_dimension)inFalse
xkzm_matmosphere_momentum_diffusivity_backgroundbackground value of momentum diffusivitym2 s-1real()kind_physinFalse
xkzm_hatmosphere_heat_diffusivity_backgroundbackground value of heat diffusivitym2 s-1real()kind_physinFalse
xkzm_sdiffusivity_background_sigma_levelsigma level threshold for background diffusivitynonereal()kind_physinFalse
dspfactke_dissipative_heating_factortke dissipative heating factornonereal()kind_physinFalse
bl_upfrupdraft_fraction_in_boundary_layer_mass_flux_schemeupdraft fraction in boundary layer mass flux schemenonereal()kind_physinFalse
bl_dnfrdowndraft_fraction_in_boundary_layer_mass_flux_schemedowndraft fraction in boundary layer mass flux schemenonereal()kind_physinFalse
ntozindex_for_ozonetracer index for ozone mixing ratioindexinteger()inFalse
du3dtcumulative_change_in_x_wind_due_to_PBLcumulative change in x wind due to PBLm s-1real(horizontal_loop_extent, vertical_dimension)kind_physinoutFalse
dv3dtcumulative_change_in_y_wind_due_to_PBLcumulative change in y wind due to PBLm s-1real(horizontal_loop_extent, vertical_dimension)kind_physinoutFalse
dt3dtcumulative_change_in_temperature_due_to_PBLcumulative change in temperature due to PBLKreal(horizontal_loop_extent, vertical_dimension)kind_physinoutFalse
dq3dtcumulative_change_in_water_vapor_specific_humidity_due_to_PBLcumulative change in water vapor specific humidity due to PBLkg kg-1real(horizontal_loop_extent, vertical_dimension)kind_physinoutFalse
do3dtcumulative_change_in_ozone_mixing_ratio_due_to_PBLcumulative change in ozone mixing ratio due to PBLkg kg-1real(horizontal_loop_extent, vertical_dimension)kind_physinoutFalse
gen_tendflag_for_generic_planetary_boundary_layer_tendencytrue if GFS_PBL_generic should calculate tendenciesflaglogical()inFalse
ldiag3dflag_diagnostics_3Dflag for 3d diagnostic fieldsflaglogical()inoutFalse
qdiag3dflag_tracer_diagnostics_3Dflag for 3d tracer diagnostic fieldsflaglogical()inoutFalse
errmsgccpp_error_messageerror message for error handling in CCPPnonecharacter()len=*outFalse
errflgccpp_error_flagerror flag for error handling in CCPPflaginteger()outFalse

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

  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) [163] ).
  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) [181]
  • 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) [174] 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) [176]

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

to be less than or equal to xkzinv

Compute an asymtotic mixing length

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

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: