diff --git a/schemes/hack_shallow/hack_convect_shallow.meta b/schemes/hack_shallow/hack_convect_shallow.meta
index 3b90a151..8d07760a 100644
--- a/schemes/hack_shallow/hack_convect_shallow.meta
+++ b/schemes/hack_shallow/hack_convect_shallow.meta
@@ -185,7 +185,7 @@
dimensions = (horizontal_loop_extent, vertical_layer_dimension)
intent = in
[ qpert_in ]
- standard_name = convective_water_vapor_wrt_moist_air_and_condensed_water_perturbation_due_to_pbl_eddies
+ standard_name = convective_water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water_perturbation_due_to_pbl_eddies
units = kg kg-1
type = real | kind = kind_phys
dimensions = (horizontal_loop_extent)
diff --git a/schemes/holtslag_boville/holtslag_boville_diff.F90 b/schemes/holtslag_boville/holtslag_boville_diff.F90
new file mode 100644
index 00000000..a7a1fef4
--- /dev/null
+++ b/schemes/holtslag_boville/holtslag_boville_diff.F90
@@ -0,0 +1,718 @@
+! Module to compute mixing coefficients associated with turbulence in the
+! planetary boundary layer and elsewhere. PBL coefficients are based on Holtslag
+! and Boville, 1991.
+!
+! Original authors:
+! Standardized: J. Rosinski, June 1992
+! Reviewed: P. Rasch, B. Boville, August 1992
+! Reviewed: P. Rasch, April 1996
+! Reviewed: B. Boville, April 1996
+! Rewritten: B. Boville, May 2000
+! Rewritten: B. Stevens, August 2000
+! Modularized: J. McCaa, September 2004
+! Functional Updates from M. Waxmonsky, February 2025
+! CCPP-ized: Haipeng Lin, May 2025
+module holtslag_boville_diff
+ use ccpp_kinds, only: kind_phys
+
+ implicit none
+ private
+ save
+
+ ! public CCPP-compliant subroutines
+ public :: holtslag_boville_diff_init
+ public :: hb_pbl_independent_coefficients_run ! formerly trbintd
+ public :: hb_pbl_dependent_coefficients_run ! formerly pblintd (also used by CLUBB)
+ public :: hb_diff_exchange_coefficients_run ! formerly austausch_pbl and compute_hb_diff logic
+ public :: hb_diff_free_atm_exchange_coefficients_run ! formerly compute_hb_free_atm_diff logic
+ public :: holtslag_boville_diff_finalize
+
+ ! tuning parameters used by HB and related subroutines
+ ! PBL limits
+ real(kind_phys), parameter :: pblmaxp = 4.e4_kind_phys ! PBL max depth [Pa]
+ real(kind_phys), parameter :: zkmin = 0.01_kind_phys ! Minimum kneutral*f(ri)
+
+ ! PBL parameters
+ real(kind_phys), parameter :: onet = 1._kind_phys/3._kind_phys ! 1/3 power in wind gradient expression
+ real(kind_phys), parameter :: betam = 15.0_kind_phys ! Constant in wind gradient expression
+ real(kind_phys), parameter :: betas = 5.0_kind_phys ! Constant in surface layer gradient expression
+ real(kind_phys), parameter :: betah = 15.0_kind_phys ! Constant in temperature gradient expression
+ real(kind_phys), parameter :: fakn = 7.2_kind_phys ! Constant in turbulent prandtl number
+ real(kind_phys), parameter :: fak = 8.5_kind_phys ! Constant in surface temperature excess
+ real(kind_phys), parameter :: ricr = 0.3_kind_phys ! Critical richardson number
+ real(kind_phys), parameter :: sffrac = 0.1_kind_phys ! Surface layer fraction of boundary layer
+ real(kind_phys), parameter :: binm = betam*sffrac ! betam * sffrac
+ real(kind_phys), parameter :: binh = betah*sffrac ! betah * sffrac
+
+ ! derived parameters from initialization
+ real(kind_phys), allocatable :: ml2(:) ! mixing lengths squared [m^2] at interfaces
+ real(kind_phys) :: ccon ! fak * sffrac * karman
+ integer :: npbl ! maximum # of levels in PBL from surface
+ integer :: ntop_turb ! top level to which turbulent vertical diffusion is applied
+ integer :: nbot_turb ! bottom level to which turbulent vertical diffusion is applied
+
+contains
+
+ ! Initialization subroutine to set module parameters, including mixing range
+!> \section arg_table_holtslag_boville_diff_init Argument Table
+!! \htmlinclude arg_table_holtslag_boville_diff_init.html
+ subroutine holtslag_boville_diff_init( &
+ amIRoot, iulog, &
+ pver, pverp, &
+ karman, &
+ pref_mid, &
+ is_hbr_pbl_scheme, &
+ ntop_turb_in, &
+ errmsg, errflg)
+
+ ! Input arguments
+ logical, intent(in) :: amIRoot ! are we on the MPI root task?
+ integer, intent(in) :: iulog ! log output unit
+ integer, intent(in) :: pver
+ integer, intent(in) :: pverp
+ real(kind_phys), intent(in) :: karman ! von_karman_constant [1]
+ real(kind_phys), intent(in) :: pref_mid(:) ! reference_pressure_in_atmosphere_layer [Pa]
+ logical, intent(in) :: is_hbr_pbl_scheme ! is HBR = true; is HB = false [flag]
+ integer, intent(in) :: ntop_turb_in ! top turbulence level [index]
+
+ ! Output arguments
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ ! Local variables
+ integer :: k
+
+ errmsg = ''
+ errflg = 0
+
+ ccon = fak * sffrac * karman
+
+ ! ntop_turb must be 1 or <= nbot_molec (lowest vertical level with molecular diffusion)
+ ! In WACCM-X, when 'waccmx_mode' is 'ionosphere' or 'neutral', this should be set to
+ ! press_lim_idx(ntop_eddy_pres, top=.true.) where ntop_eddy_pres = 1.e-7_kind_phys [Pa]
+ !
+ ! FIXME: this is not as clean as I would like it to be
+ ! but I want to avoid depending on host model press_lim_idx; so the index is passed
+ ! into this initialization subroutine. This also avoids having to query something like
+ ! 'waccmx_is' or 'waccmx_mode' that is host model specific inside this scheme
+ ! to make it portable, but this means the flow on the host model side has to
+ ! specify this top turbulence index. (hplin, 5/8/25)
+ ntop_turb = ntop_turb_in
+
+ nbot_turb = pver
+
+ ! allocate and populate mixing lengths squared to grid vertical interfaces
+ allocate(ml2(pverp), stat=errflg, errmsg=errmsg)
+ if(errflg /= 0) return
+
+ ml2(ntop_turb) = 0._kind_phys
+ do k = ntop_turb+1, nbot_turb
+ ml2(k) = 30.0_kind_phys**2 ! HB scheme: length scale = 30m
+ if(is_hbr_pbl_scheme) then
+ ml2(k) = 1.0_kind_phys**2 ! HBR scheme: length scale = 1m
+ end if
+ end do
+ ml2(nbot_turb+1) = 0._kind_phys
+
+ ! Limit pbl height to regions below 400 mb
+ ! npbl = max # of levels (from bottom) in pbl
+ npbl = 0
+ do k = nbot_turb,ntop_turb,-1
+ if (pref_mid(k) >= pblmaxp) then
+ npbl = npbl + 1
+ end if
+ end do
+ npbl = max(npbl, 1)
+
+ if(amIRoot) then
+ write(iulog,*) 'Holtslag-Boville PBL: PBL height will be limited to bottom ', npbl, &
+ ' model levels. Top is ', pref_mid(pver+1-npbl), ' pascals'
+ write(iulog,*) 'Holtslag-Boville PBL: top level of turbulence is ', ntop_turb, &
+ ' and bottom is ', nbot_turb
+ end if
+
+ end subroutine holtslag_boville_diff_init
+
+ ! Computation of time-dependent, PBL height-independent variables for HB mixing.
+ ! Original author: B. Stevens, rewrite August 2000
+!> \section arg_table_hb_pbl_independent_coefficients_run Argument Table
+!! \htmlinclude arg_table_hb_pbl_independent_coefficients_run.html
+ pure subroutine hb_pbl_independent_coefficients_run( &
+ ncol, pver, &
+ zvir, rair, cpair, gravit, karman, &
+ exner, t, &
+ q_wv, &
+ z, &
+ pmid, &
+ u, v, &
+ taux, tauy, &
+ shflx, q_wv_flx, &
+ ! below output
+ thv, ustar, &
+ khfs, kqfs, kbfs, &
+ obklen, s2, ri, &
+ errmsg, errflg)
+ use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_friction_velocity, calc_obukhov_length, &
+ calc_ideal_gas_rrho, &
+ calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, &
+ calc_kinematic_buoyancy_flux
+
+ ! Input arguments
+ integer, intent(in) :: ncol ! # of atmospheric columns
+ integer, intent(in) :: pver ! # of vertical levels
+ real(kind_phys), intent(in) :: zvir
+ real(kind_phys), intent(in) :: rair
+ real(kind_phys), intent(in) :: cpair
+ real(kind_phys), intent(in) :: gravit
+ real(kind_phys), intent(in) :: karman
+ real(kind_phys), intent(in) :: exner(:,:) ! exner [1]
+ real(kind_phys), intent(in) :: t (:,:) ! temperature [K]
+ real(kind_phys), intent(in) :: q_wv(:,:) ! specific humidity [kg kg-1]
+ real(kind_phys), intent(in) :: z (:,:) ! height above surface [m]
+ real(kind_phys), intent(in) :: u (:,:) ! zonal velocity [m s-1]
+ real(kind_phys), intent(in) :: v (:,:) ! meridional velocity [m s-1]
+ real(kind_phys), intent(in) :: taux (:) ! zonal stress [N m-2]
+ real(kind_phys), intent(in) :: tauy (:) ! meridional stress [N m-2]
+ real(kind_phys), intent(in) :: shflx(:) ! sensible heat flux [W m-2]
+ real(kind_phys), intent(in) :: q_wv_flx(:) ! upward water vapor flux at surface [kg m-2 s-1]
+ real(kind_phys), intent(in) :: pmid(:,:) ! midpoint pressures
+
+ ! Output arguments
+ real(kind_phys), intent(out) :: thv(:,:) ! virtual potential temperature [K]
+ real(kind_phys), intent(out) :: ustar(:) ! surface friction velocity [m s-1]
+ real(kind_phys), intent(out) :: khfs(:) ! kinematic surface heat flux [K m s-1]
+ real(kind_phys), intent(out) :: kqfs(:) ! kinematic surface water vapor flux [kg kg-1 m s-1]
+ real(kind_phys), intent(out) :: kbfs(:) ! surface kinematic buoyancy flux [m^2 s-3]
+ real(kind_phys), intent(out) :: obklen(:) ! Obukhov length [m]
+ real(kind_phys), intent(out) :: s2(:,:) ! shear squared [s-2]
+ real(kind_phys), intent(out) :: ri(:,:) ! richardson number: n2/s2 [1]
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ integer :: i, k
+ real(kind_phys) :: dvdz2 ! velocity shear squared [m^2 s-2]
+ real(kind_phys) :: dz ! delta z between midpoints [m]
+ real(kind_phys) :: rrho(ncol) ! 1 / bottom level density [m^3 kg-1]
+ real(kind_phys) :: n2(ncol, pver) ! brunt vaisaila frequency [s-2]
+ real(kind_phys) :: th(ncol, pver) ! potential temperature [K]
+
+ errmsg = ''
+ errflg = 0
+
+ ! calculate potential temperature [K]
+ th(:ncol,:) = t(:ncol,:) * exner(:ncol,:)
+
+ ! virtual temperature
+ thv(:,:) = 0._kind_phys
+ thv(:ncol,ntop_turb:) = calc_virtual_temperature(th(:ncol,ntop_turb:), q_wv(:ncol,ntop_turb:), zvir)
+
+ ! Compute ustar, Obukhov length, and kinematic surface fluxes.
+ rrho(:ncol) = calc_ideal_gas_rrho(rair, t(:ncol,pver), pmid(:ncol,pver))
+ ustar(:ncol) = calc_friction_velocity(taux(:ncol),tauy(:ncol), rrho(:ncol))
+ khfs(:ncol) = calc_kinematic_heat_flux(shflx(:ncol), rrho(:ncol), cpair)
+ kqfs(:ncol) = calc_kinematic_water_vapor_flux(q_wv_flx(:ncol), rrho(:ncol))
+ kbfs(:ncol) = calc_kinematic_buoyancy_flux(khfs(:ncol), zvir, th(:ncol,pver), kqfs(:ncol))
+ obklen(:ncol) = calc_obukhov_length(thv(:ncol,pver), ustar(:ncol), gravit, karman, kbfs(:ncol))
+
+ ! Compute s^2 (shear squared), n^2 (brunt vaisaila frequency), and ri (Richardson number = n^2/s^2)
+ ! using virtual temperature
+ ! (formerly trbintd)
+ do k = ntop_turb,nbot_turb-1
+ do i = 1,ncol
+ dvdz2 = (u(i,k)-u(i,k+1))**2 + (v(i,k)-v(i,k+1))**2
+ dvdz2 = max(dvdz2,1.e-36_kind_phys)
+ dz = z(i,k) - z(i,k+1)
+ s2(i,k) = dvdz2/(dz**2)
+ n2(i,k) = gravit*2.0_kind_phys*(thv(i,k) - thv(i,k+1))/((thv(i,k) + thv(i,k+1))*dz)
+ ri(i,k) = n2(i,k)/s2(i,k)
+ end do
+ end do
+ end subroutine hb_pbl_independent_coefficients_run
+
+ ! Compute time-dependent, PBL height-dependent variables
+ ! (formerly pblintd)
+ ! The PBL depth follows:
+ ! Holtslag, A.A.M., and B.A. Boville, 1993:
+ ! Local versus Nonlocal Boundary-Layer Diffusion in a Global Climate Model.
+ ! J. Clim., vol. 6., p. 1825--1842. https://doi.org/10.1175/1520-0442(1993)006<1825:LVNBLD>2.0.CO;2
+ !
+ ! Updated by Holtslag and Hack to exclude the surface layer from the
+ ! definition of the boundary layer Richardson number. Ri is now defined
+ ! across the outer layer of the pbl (between the top of the surface
+ ! layer and the pbl top) instead of the full pbl (between the surface and
+ ! the pbl top). For simiplicity, the surface layer is assumed to be the
+ ! region below the first model level (otherwise the boundary layer depth
+ ! determination would require iteration).
+ !
+ ! Original author: B. Stevens, August 2000, extracted from pbldiff
+!> \section arg_table_hb_pbl_dependent_coefficients_run Argument Table
+!! \htmlinclude arg_table_hb_pbl_dependent_coefficients_run.html
+ pure subroutine hb_pbl_dependent_coefficients_run( &
+ ncol, pver, pverp, &
+ gravit, &
+ z, zi, &
+ u, v, &
+ cldn, &
+ ! below inputs from pbl_independent_coefficients
+ thv, ustar, kbfs, obklen, &
+ ! below output
+ pblh, &
+ wstar, bge, &
+ errmsg, errflg)
+
+ ! Input arguments
+ integer, intent(in) :: ncol ! # of atmospheric columns
+ integer, intent(in) :: pver ! # of vertical levels
+ integer, intent(in) :: pverp ! # of vertical level interfaces
+ real(kind_phys), intent(in) :: gravit
+ real(kind_phys), intent(in) :: z (:,:) ! height above surface [m]
+ real(kind_phys), intent(in) :: zi (:,:) ! height above surface [m], interfaces
+ real(kind_phys), intent(in) :: u (:,:) ! zonal velocity [m s-1]
+ real(kind_phys), intent(in) :: v (:,:) ! meridional velocity [m s-1]
+ real(kind_phys), intent(in) :: cldn(:,:) ! stratiform cloud fraction [fraction]
+
+ ! Input arguments (output from hb_pbl_independent_coefficients)
+ real(kind_phys), intent(in) :: thv (:,:) ! virtual potential temperature [K]
+ real(kind_phys), intent(in) :: ustar(:) ! surface friction velocity [m s-1]
+ real(kind_phys), intent(in) :: kbfs(:) ! surface kinematic buoyancy flux [m^2 s-3]
+ real(kind_phys), intent(in) :: obklen(:) ! Obukhov length [m]
+
+ ! Output arguments
+ real(kind_phys), intent(out) :: pblh(:) ! boundary-layer height [m]
+ real(kind_phys), intent(out) :: wstar(:) ! convective scale velocity [m s-1]
+ real(kind_phys), intent(out) :: bge(:) ! buoyancy gradient enhancement [m s-2]
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ ! Local variables
+ integer :: i, k
+ real(kind_phys) :: phiminv(ncol) ! inverse phi function for momentum
+ real(kind_phys) :: rino(ncol,pver) ! bulk Richardson no. from level to ref lev
+ real(kind_phys) :: tlv(ncol) ! ref. level potential tmp + tmp excess
+ real(kind_phys) :: vvk ! velocity magnitude squared
+
+ logical :: check(ncol) ! false if Richardson number > critical
+ logical :: ocncldcheck(ncol) ! true if ocean surface (not implemented) and cloud in lowest layer
+
+ ! Local parameters
+ real(kind_phys), parameter :: tiny = 1.e-36_kind_phys ! lower bound for wind magnitude
+ real(kind_phys), parameter :: fac = 100._kind_phys ! ustar parameter in height diagnosis
+
+ errmsg = ''
+ errflg = 0
+
+ ! Compute Obukhov length virtual temperature flux and various arrays for use later:
+ do i=1,ncol
+ check(i) = .true.
+ rino(i,pver) = 0.0_kind_phys
+ pblh(i) = z(i,pver)
+ end do
+
+ ! PBL height calculation: Scan upward until the Richardson number between
+ ! the first level and the current level exceeds the "critical" value.
+ do k=pver-1,pver-npbl+1,-1
+ do i=1,ncol
+ if (check(i)) then
+ vvk = (u(i,k) - u(i,pver))**2 + (v(i,k) - v(i,pver))**2 + fac*ustar(i)**2
+ vvk = max(vvk,tiny)
+ rino(i,k) = gravit*(thv(i,k) - thv(i,pver))*(z(i,k)-z(i,pver))/(thv(i,pver)*vvk)
+ ! Modified for boundary layer height diagnosis: Bert Holtslag, June 1994
+ ! >>>>>>>>> (Use ricr = 0.3 in this formulation)
+ if (rino(i,k) >= ricr) then
+ pblh(i) = z(i,k+1) + (ricr - rino(i,k+1))/(rino(i,k) - rino(i,k+1)) * &
+ (z(i,k) - z(i,k+1))
+ check(i) = .false.
+ end if
+ end if
+ end do
+ end do
+
+ ! Estimate an effective surface temperature to account for surface fluctuations
+ do i=1,ncol
+ if (check(i)) pblh(i) = z(i,pverp-npbl)
+ check(i) = (kbfs(i) > 0._kind_phys)
+ if (check(i)) then
+ phiminv(i) = (1._kind_phys - binm*pblh(i)/obklen(i))**onet
+ rino(i,pver) = 0.0_kind_phys
+ tlv(i) = thv(i,pver) + kbfs(i)*fak/(ustar(i)*phiminv(i))
+ end if
+ end do
+
+ ! Improve pblh estimate for unstable conditions using the convective temperature excess:
+ do i = 1,ncol
+ bge(i) = 1.e-8_kind_phys
+ end do
+
+ do k=pver-1,pver-npbl+1,-1
+ do i=1,ncol
+ if (check(i)) then
+ vvk = (u(i,k) - u(i,pver))**2 + (v(i,k) - v(i,pver))**2 + fac*ustar(i)**2
+ vvk = max(vvk,tiny)
+ rino(i,k) = gravit*(thv(i,k) - tlv(i))*(z(i,k)-z(i,pver))/(thv(i,pver)*vvk)
+ if (rino(i,k) >= ricr) then
+ pblh(i) = z(i,k+1) + (ricr - rino(i,k+1))/(rino(i,k) - rino(i,k+1)) * &
+ (z(i,k) - z(i,k+1))
+ bge(i) = 2._kind_phys*gravit/(thv(i,k)+thv(i,k+1))*(thv(i,k)-thv(i,k+1))/(z(i,k)-z(i,k+1))*pblh(i)
+ if (bge(i).lt.0._kind_phys) then
+ bge(i) = 1.e-8_kind_phys
+ endif
+ check(i) = .false.
+ end if
+ end if
+ end do
+ end do
+
+ ! PBL height must be greater than some minimum mechanical mixing depth
+ ! Several investigators have proposed minimum mechanical mixing depth
+ ! relationships as a function of the local friction velocity, u*. We
+ ! make use of a linear relationship of the form h = c u* where c=700.
+ ! The scaling arguments that give rise to this relationship most often
+ ! represent the coefficient c as some constant over the local coriolis
+ ! parameter. Here we make use of the experimental results of Koracin
+ ! and Berkowicz (1988) [BLM, Vol 43]: https://doi.org/10.1007/BF00153969
+ ! for wich they recommend 0.07/f
+ ! where f was evaluated at 39.5 N and 52 N. Thus we use a typical mid
+ ! latitude value for f so that c = 0.07/f = 700. Also, do not allow
+ ! PBL to exceed some maximum (npbl) number of allowable points
+ do i=1,ncol
+ if (check(i)) pblh(i) = z(i,pverp-npbl)
+ pblh(i) = max(pblh(i),700.0_kind_phys*ustar(i))
+ wstar(i) = (max(0._kind_phys,kbfs(i))*gravit*pblh(i)/thv(i,pver))**onet
+ end do
+
+ ! Final requirement on PBL height is that it must be greater than the depth
+ ! of the lowest model level over ocean if there is any cloud diagnosed in
+ ! the lowest model level. This is to deal with the inadequacies of the
+ ! current "dry" formulation of the boundary layer, where this test is
+ ! used to identify circumstances where there is marine stratus in the
+ ! lowest level, and to provide a weak ventilation of the layer to avoid
+ ! a pathology in the cloud scheme (locking in low-level stratiform cloud)
+ ! If over an ocean surface, and any cloud is diagnosed in the
+ ! lowest level, set pblh to 50 meters higher than top interface of lowest level
+ !
+ ! jrm This is being applied everywhere (not just ocean)!
+ do i=1,ncol
+ ocncldcheck(i) = .false.
+ if (cldn(i,pver).ge.0.0_kind_phys) ocncldcheck(i) = .true.
+ if (ocncldcheck(i)) pblh(i) = max(pblh(i),zi(i,pver) + 50._kind_phys)
+ end do
+ end subroutine hb_pbl_dependent_coefficients_run
+
+ ! Atmosphere boundary layer computation.
+ !
+ ! Nonlocal scheme that determines eddy diffusivities based on a
+ ! specified boundary layer height and a turbulent velocity scale;
+ ! also, countergradient effects for heat and moisture, and constituents
+ ! are included, along with temperature and humidity perturbations which
+ ! measure the strength of convective thermals in the lower part of the
+ ! atmospheric boundary layer.
+ !
+ ! For more information, see Holtslag, A.A.M., and B.A. Boville, 1993:
+ ! Local versus Nonlocal Boundary-Layer Diffusion in a Global Climate Model.
+ ! J. Clim., vol. 6., p. 1825--1842. https://doi.org/10.1175/1520-0442(1993)006<1825:LVNBLD>2.0.CO;2
+ !
+ ! Updated by Holtslag and Hack to exclude the surface layer from the
+ ! definition of the boundary layer Richardson number. Ri is now defined
+ ! across the outer layer of the pbl (between the top of the surface
+ ! layer and the pbl top) instead of the full pbl (between the surface and
+ ! the pbl top). For simiplicity, the surface layer is assumed to be the
+ ! region below the first model level (otherwise the boundary layer depth
+ ! determination would require iteration).
+ !
+ ! Original authors: B. Boville, B. Stevens, rewrite August 2000
+!> \section arg_table_hb_diff_exchange_coefficients_run Argument Table
+!! \htmlinclude arg_table_hb_diff_exchange_coefficients_run.html
+ pure subroutine hb_diff_exchange_coefficients_run( &
+ ncol, pver, pverp, &
+ karman, cpair, &
+ z, &
+ is_hbr_pbl_scheme, &
+ ! input from hb_pbl_independent_coefficients
+ kqfs, khfs, kbfs, &
+ ustar, obklen, &
+ s2, ri, &
+ ! input from hb_pbl_dependent_coefficients
+ pblh, wstar, bge, &
+ ! below output
+ kvm, kvh, kvq, &
+ cgh, cgs, &
+ tpert, qpert, &
+ tke, &
+ errmsg, errflg)
+
+ use atmos_phys_pbl_utils, only: calc_eddy_flux_coefficient
+
+ ! Input arguments
+ integer, intent(in) :: ncol ! # of atmospheric columns
+ integer, intent(in) :: pver ! # of vertical levels
+ integer, intent(in) :: pverp ! # of vertical level interfaces
+ real(kind_phys), intent(in) :: karman
+ real(kind_phys), intent(in) :: cpair
+ real(kind_phys), intent(in) :: z (:,:) ! height above surface [m]
+ logical, intent(in) :: is_hbr_pbl_scheme
+
+ real(kind_phys), intent(in) :: khfs(:) ! kinematic surface heat flux [K m s-1]
+ real(kind_phys), intent(in) :: kqfs(:) ! kinematic surface constituent flux [kg kg-1 m s-1]
+ real(kind_phys), intent(in) :: kbfs(:) ! surface kinematic buoyancy flux [m^2 s-3]
+ real(kind_phys), intent(in) :: ustar(:) ! surface friction velocity [m s-1]
+ real(kind_phys), intent(in) :: obklen(:) ! Obukhov length [m]
+ real(kind_phys), intent(in) :: s2(:,:) ! shear squared [s-2]
+ real(kind_phys), intent(in) :: ri(:,:) ! richardson number: n2/s2 [1]
+
+ real(kind_phys), intent(in) :: pblh(:) ! boundary-layer height [m]
+ real(kind_phys), intent(in) :: wstar(:) ! convective scale velocity [m s-1]
+ real(kind_phys), intent(in) :: bge(:) ! buoyancy gradient enhancement [m s-2]
+
+ ! Output variables
+ real(kind_phys), intent(out) :: kvm(:,:) ! eddy diffusivity for momentum [m^2 s-1], interfaces
+ real(kind_phys), intent(out) :: kvh(:,:) ! eddy diffusivity for heat [m^2 s-1], interfaces
+ real(kind_phys), intent(out) :: kvq(:,:) ! eddy diffusivity for constituents [m^2 s-1], interfaces
+ real(kind_phys), intent(out) :: cgh(:,:) ! counter-gradient term for heat [J kg-1 m-1], interfaces
+ real(kind_phys), intent(out) :: cgs(:,:) ! counter-gradient star (cg/flux) [s m-2], interfaces
+ real(kind_phys), intent(out) :: tpert(:) ! convective temperature excess [K]
+ real(kind_phys), intent(out) :: qpert(:) ! convective humidity excess [kg kg-1]
+ real(kind_phys), intent(out) :: tke(:,:) ! (estimated) turbulent kinetic energy [m^2 s-2], interfaces
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ ! Local variables
+ integer :: i, k
+ real(kind_phys) :: kvf(ncol,pverp) ! free atmospheric eddy diffusivity [m^2 s-1]
+ integer :: ktopbl(ncol) ! index of first midpoint inside PBL (diagnostic) [index]
+
+ real(kind_phys) :: phiminv(ncol) ! inverse phi function for momentum
+ real(kind_phys) :: phihinv(ncol) ! inverse phi function for heat
+ real(kind_phys) :: wm(ncol) ! turbulent velocity scale for momentum
+ real(kind_phys) :: zp(ncol) ! current level height + one level up
+ real(kind_phys) :: fak1(ncol) ! k*ustar*pblh
+ real(kind_phys) :: fak2(ncol) ! k*wm*pblh
+ real(kind_phys) :: fak3(ncol) ! fakn*wstar/wm
+ real(kind_phys) :: pblk(ncol) ! level eddy diffusivity for momentum
+ real(kind_phys) :: pr(ncol) ! Prandtl number for eddy diffusivities
+ real(kind_phys) :: zl(ncol) ! zmzp / Obukhov length
+ real(kind_phys) :: zh(ncol) ! zmzp / pblh
+ real(kind_phys) :: zzh(ncol) ! (1-(zmzp/pblh))**2
+ real(kind_phys) :: zmzp ! level height halfway between zm and zp
+ real(kind_phys) :: term ! intermediate calculation
+ real(kind_phys) :: kve ! diffusivity at entrainment layer in unstable cases [m s-2]
+
+ logical :: unstable(ncol) ! points with unstable pbl (positive virtual heat flux) [count]
+ logical :: pblpt(ncol) ! points within pbl
+
+ errmsg = ''
+ errflg = 0
+
+ ! Get atmosphere exchange coefficients
+ kvf(:ncol,:) = 0.0_kind_phys
+ do k = ntop_turb, nbot_turb-1
+ do i = 1, ncol
+ kvf(i,k+1) = calc_eddy_flux_coefficient(ml2(k), ri(i, k), s2(i, k))
+ end do
+ end do
+
+ ! Compute PBL exchange coefficients (formerly austausch_pbl)
+ kvm = 0._kind_phys
+ kvh = 0._kind_phys
+ kve = 0._kind_phys
+ cgh = 0._kind_phys
+ cgs = 0._kind_phys
+ tpert = 0._kind_phys
+ qpert = 0._kind_phys
+ ktopbl = 0._kind_phys
+ tke = 0._kind_phys
+
+ ! Initialize height independent arrays
+ do i=1,ncol
+ unstable(i) = (kbfs(i) > 0._kind_phys)
+ pblk(i) = 0.0_kind_phys
+ fak1(i) = ustar(i)*pblh(i)*karman
+ if (unstable(i)) then
+ phiminv(i) = (1._kind_phys - binm*pblh(i)/obklen(i))**onet
+ phihinv(i) = sqrt(1._kind_phys - binh*pblh(i)/obklen(i))
+ wm(i) = ustar(i)*phiminv(i)
+ fak2(i) = wm(i)*pblh(i)*karman
+ fak3(i) = fakn*wstar(i)/wm(i)
+ tpert(i) = max(khfs(i)*fak/wm(i),0._kind_phys)
+ qpert(i) = max(kqfs(i)*fak/wm(i),0._kind_phys)
+ else
+ tpert(i) = max(khfs(i)*fak/ustar(i),0._kind_phys)
+ qpert(i) = max(kqfs(i)*fak/ustar(i),0._kind_phys)
+ end if
+ end do
+
+ ! Initialize output arrays with free atmosphere values
+ do k=1,pverp
+ do i=1,ncol
+ kvm(i,k) = kvf(i,k)
+ kvh(i,k) = kvf(i,k)
+ cgh(i,k) = 0.0_kind_phys
+ cgs(i,k) = 0.0_kind_phys
+ end do
+ end do
+
+ ! Main level loop to compute the diffusivities and counter-gradient terms. These terms are
+ ! only calculated at points determined to be in the interior of the pbl (pblpt(i)==.true.),
+ ! and then calculations are directed toward regime: stable vs unstable, surface vs outer
+ ! layer.
+ do k=pver,pver-npbl+2,-1
+ do i=1,ncol
+ pblpt(i) = (z(i,k) < pblh(i))
+ if (pblpt(i)) then
+ ktopbl(i) = k
+ zp(i) = z(i,k-1)
+ if (zkmin == 0.0_kind_phys .and. zp(i) > pblh(i)) then
+ zp(i) = pblh(i)
+ endif
+ zmzp = 0.5_kind_phys*(z(i,k) + zp(i)) ! we think this is an approximation to the interface height (where KVs are calculated)
+ zh(i) = zmzp/pblh(i)
+ zl(i) = zmzp/obklen(i)
+ zzh(i) = zh(i)*max(0._kind_phys,(1._kind_phys - zh(i)))**2
+ if (unstable(i)) then
+ if (zh(i) < sffrac) then
+ term = (1._kind_phys - betam*zl(i))**onet
+ pblk(i) = fak1(i)*zzh(i)*term
+ pr(i) = term/sqrt(1._kind_phys - betah*zl(i))
+ else
+ pblk(i) = fak2(i)*zzh(i)
+ pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
+ cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
+ cgh(i,k) = khfs(i)*cgs(i,k)*cpair
+ end if
+ else
+ if (zl(i) <= 1._kind_phys) then
+ pblk(i) = fak1(i)*zzh(i)/(1._kind_phys + betas*zl(i))
+ else
+ pblk(i) = fak1(i)*zzh(i)/(betas + zl(i))
+ end if
+ pr(i) = 1._kind_phys
+ end if
+ kvm(i,k) = max(pblk(i),kvf(i,k))
+ kvh(i,k) = max(pblk(i)/pr(i),kvf(i,k))
+ end if
+ end do
+ end do
+
+ ! Check whether last allowed midpoint is within PBL
+ ! HBR scheme only
+ if(is_hbr_pbl_scheme) then
+ ! apply new diffusivity at entrainment zone
+ do i = 1,ncol
+ if (bge(i) > 1.e-7_kind_phys) then
+ k = ktopbl(i)
+ kve = 0.2_kind_phys*(wstar(i)**3+5._kind_phys*ustar(i)**3)/bge(i)
+ kvm(i,k) = kve
+ kvh(i,k) = kve
+ end if
+ end do
+ end if
+
+ ! Crude estimate of tke (tke=0 above boundary layer)
+ do k = max(pverp-npbl,2),pverp
+ do i = 1, ncol
+ if (z(i,k-1) < pblh(i)) then
+ tke(i,k) = (kvm(i,k) / pblh(i))**2
+ endif
+ end do
+ end do
+
+ ! Finalize
+ kvq(:ncol,:) = kvh(:ncol,:)
+ end subroutine hb_diff_exchange_coefficients_run
+
+ ! A version of hb_diff_exchange_coefficients
+ ! that only computes free atmosphere exchanges (no PBL computations)
+ !
+ ! Original authors: B. Stevens, rewrite August 2000
+ ! Thomas Toniazzo, Peter H. Lauritzen, June 2023
+!> \section arg_table_hb_diff_free_atm_exchange_coefficients_run Argument Table
+!! \htmlinclude arg_table_hb_diff_free_atm_exchange_coefficients_run.html
+ pure subroutine hb_diff_free_atm_exchange_coefficients_run( &
+ ncol, pver, pverp, &
+ ! input from hb_pbl_independent_coefficients
+ s2, ri, &
+ ! top of CLUBB (or layer above where HB applies)
+ bottom_boundary, &
+ ! below output
+ kvm, kvh, kvq, &
+ cgh, cgs, &
+ errmsg, errflg)
+
+ use atmos_phys_pbl_utils, only: calc_free_atm_eddy_flux_coefficient
+
+ ! Input arguments
+ integer, intent(in) :: ncol ! # of atmospheric columns
+ integer, intent(in) :: pver ! # of vertical levels
+ integer, intent(in) :: pverp ! # of vertical level interfaces
+ real(kind_phys), intent(in) :: s2(:,:) ! shear squared [s-2]
+ real(kind_phys), intent(in) :: ri(:,:) ! richardson number: n2/s2 [1]
+
+ !REMOVECAM: change to integer once CAM snapshots are no longer used.
+ real(kind_phys), intent(in) :: bottom_boundary(:) ! level above which HB runs [index]
+
+ ! Output variables
+ real(kind_phys), intent(out) :: kvm(:,:) ! eddy diffusivity for momentum [m^2 s-1], interfaces
+ real(kind_phys), intent(out) :: kvh(:,:) ! eddy diffusivity for heat [m^2 s-1], interfaces
+ real(kind_phys), intent(out) :: kvq(:,:) ! eddy diffusivity for constituents [m^2 s-1], interfaces
+ real(kind_phys), intent(out) :: cgh(:,:) ! counter-gradient term for heat [J kg-1 m-1], interfaces
+ real(kind_phys), intent(out) :: cgs(:,:) ! counter-gradient star (cg/flux) [s m-2], interfaces
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ ! Local variables
+ integer :: i, k
+ real(kind_phys) :: kvf(ncol,pverp) ! free atmospheric eddy diffusivity [m^2 s-1]
+ integer :: bottom_boundary_int(ncol)
+
+ errmsg = ''
+ errflg = 0
+
+ ! Get free atmosphere exchange coefficients
+ kvf(:ncol,:) = 0.0_kind_phys
+ do k = ntop_turb, nbot_turb-1
+ do i = 1, ncol
+ kvf(i,k+1) = calc_free_atm_eddy_flux_coefficient(ml2(k), ri(i, k), s2(i, k))
+ end do
+ end do
+
+ kvq(:ncol,:) = kvf(:ncol,:)
+ kvm(:ncol,:) = kvf(:ncol,:)
+ kvh(:ncol,:) = kvf(:ncol,:)
+ cgh(:ncol,:) = 0._kind_phys
+ cgs(:ncol,:) = 0._kind_phys
+
+ ! HB coefficients will be zeroed out in the layers below
+ ! the bottom_boundary, which is the topmost layer where CLUBB is active.
+ ! 1 --- TOA ---
+ ! 2 .. HB (free atm) ..
+ ! 3
+ ! ...
+ ! --- bottom_boundary --- | levels here
+ ! ... .. CLUBB .. | have HB coefficients
+ ! ... | zeroed out as CLUBB is active.
+ ! pverp --- surface ---
+ bottom_boundary_int(:ncol) = int(bottom_boundary(:ncol))
+ do i = 1, ncol
+ do k = bottom_boundary_int(i), pverp
+ kvm(i,k) = 0._kind_phys
+ kvh(i,k) = 0._kind_phys
+ kvq(i,k) = 0._kind_phys
+ end do
+ end do
+ end subroutine hb_diff_free_atm_exchange_coefficients_run
+
+!> \section arg_table_holtslag_boville_diff_finalize Argument Table
+!! \htmlinclude arg_table_holtslag_boville_diff_finalize.html
+ subroutine holtslag_boville_diff_finalize(errmsg, errflg)
+ ! Output arguments
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ errmsg = ''
+ errflg = 0
+
+ if(allocated(ml2)) then
+ deallocate(ml2)
+ endif
+ end subroutine holtslag_boville_diff_finalize
+
+end module holtslag_boville_diff
diff --git a/schemes/holtslag_boville/holtslag_boville_diff.meta b/schemes/holtslag_boville/holtslag_boville_diff.meta
new file mode 100644
index 00000000..75df7fd4
--- /dev/null
+++ b/schemes/holtslag_boville/holtslag_boville_diff.meta
@@ -0,0 +1,638 @@
+[ccpp-table-properties]
+ name = holtslag_boville_diff
+ type = scheme
+
+[ccpp-arg-table]
+ name = holtslag_boville_diff_init
+ type = scheme
+[ amIRoot ]
+ standard_name = flag_for_mpi_root
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = in
+[ iulog ]
+ standard_name = log_output_unit
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pverp ]
+ standard_name = vertical_interface_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ karman ]
+ standard_name = von_karman_constant
+ units = 1
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ pref_mid ]
+ standard_name = reference_pressure_in_atmosphere_layer
+ units = Pa
+ type = real | kind = kind_phys
+ dimensions = (vertical_layer_dimension)
+ intent = in
+[ is_hbr_pbl_scheme ]
+ standard_name = flag_for_hbr_configuration_of_holtslag_boville_mixing_scheme
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = in
+[ ntop_turb_in ]
+ standard_name = vertical_layer_index_of_eddy_vertical_diffusion_top
+ units = index
+ type = integer
+ dimensions = ()
+ intent = in
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-arg-table]
+ name = holtslag_boville_diff_finalize
+ type = scheme
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+
+[ccpp-table-properties]
+ name = hb_pbl_independent_coefficients
+ type = scheme
+
+[ccpp-arg-table]
+ name = hb_pbl_independent_coefficients_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ zvir ]
+ standard_name = ratio_of_water_vapor_to_dry_air_gas_constants_minus_one
+ units = 1
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ rair ]
+ standard_name = gas_constant_of_dry_air
+ units = J kg-1 K-1
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ cpair ]
+ standard_name = specific_heat_of_dry_air_at_constant_pressure
+ units = J kg-1 K-1
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ gravit ]
+ standard_name = standard_gravitational_acceleration
+ units = m s-2
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ karman ]
+ standard_name = von_karman_constant
+ units = 1
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ exner ]
+ standard_name = reciprocal_of_dimensionless_exner_function_wrt_surface_air_pressure
+ units = 1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ t ]
+ standard_name = air_temperature
+ units = K
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ q_wv ]
+ standard_name = water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water
+ units = kg kg-1
+ advected = True
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ z ]
+ standard_name = geopotential_height_wrt_surface
+ units = m
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ pmid ]
+ standard_name = air_pressure
+ units = Pa
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ u ]
+ standard_name = eastward_wind
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ v ]
+ standard_name = northward_wind
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ taux ]
+ standard_name = eastward_total_stress_at_surface_for_holtslag_boville_boundary_layer_scheme
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ tauy ]
+ standard_name = northward_total_stress_at_surface_for_holtslag_boville_boundary_layer_scheme
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ shflx ]
+ standard_name = surface_upward_sensible_heat_flux_from_coupler
+ units = W m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ q_wv_flx ]
+ standard_name = surface_upward_water_vapor_flux
+ # surface_upward_ccpp_constituent_fluxes but only for water vapor.
+ units = kg m-2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ thv ]
+ standard_name = virtual_potential_temperature_wrt_surface_pressure_at_bottom_layer
+ units = K
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = out
+[ ustar ]
+ standard_name = surface_friction_velocity
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ khfs ]
+ standard_name = kinematic_sensible_heat_flux_at_surface
+ units = K m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ kqfs ]
+ standard_name = kinematic_water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water_flux_at_surface
+ units = kg kg-1 m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ kbfs ]
+ standard_name = kinematic_buoyancy_flux_at_surface
+ # m s-2 m s-1
+ units = m2 s-3
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ obklen ]
+ standard_name = obukhov_length
+ units = m
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ s2 ]
+ standard_name = square_of_wind_shear
+ units = s-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = out
+[ ri ]
+ standard_name = gradient_richardson_number
+ units = 1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-table-properties]
+ name = hb_pbl_dependent_coefficients
+ type = scheme
+
+[ccpp-arg-table]
+ name = hb_pbl_dependent_coefficients_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pverp ]
+ standard_name = vertical_interface_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ gravit ]
+ standard_name = standard_gravitational_acceleration
+ units = m s-2
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ z ]
+ standard_name = geopotential_height_wrt_surface
+ units = m
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ zi ]
+ standard_name = geopotential_height_wrt_surface_at_interface
+ units = m
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ u ]
+ standard_name = eastward_wind
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ v ]
+ standard_name = northward_wind
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ cldn ]
+ standard_name = stratiform_cloud_area_fraction
+ units = fraction
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ thv ]
+ standard_name = virtual_potential_temperature_wrt_surface_pressure_at_bottom_layer
+ units = K
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ ustar ]
+ standard_name = surface_friction_velocity
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ kbfs ]
+ standard_name = kinematic_buoyancy_flux_at_surface
+ # m s-2 m s-1
+ units = m2 s-3
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ obklen ]
+ standard_name = obukhov_length
+ units = m
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ pblh ]
+ standard_name = atmosphere_boundary_layer_thickness
+ units = m
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ wstar ]
+ standard_name = atmosphere_boundary_layer_convective_velocity_scale
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ bge ]
+ standard_name = atmosphere_boundary_layer_buoyancy_gradient_enhancement_for_holtslag_boville_mixing_scheme
+ units = m s-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-table-properties]
+ name = hb_diff_exchange_coefficients
+ type = scheme
+
+[ccpp-arg-table]
+ name = hb_diff_exchange_coefficients_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pverp ]
+ standard_name = vertical_interface_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ karman ]
+ standard_name = von_karman_constant
+ units = 1
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ cpair ]
+ standard_name = specific_heat_of_dry_air_at_constant_pressure
+ units = J kg-1 K-1
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ z ]
+ standard_name = geopotential_height_wrt_surface
+ units = m
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ is_hbr_pbl_scheme ]
+ standard_name = flag_for_hbr_configuration_of_holtslag_boville_mixing_scheme
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = in
+[ kqfs ]
+ standard_name = kinematic_water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water_flux_at_surface
+ units = kg kg-1 m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ khfs ]
+ standard_name = kinematic_sensible_heat_flux_at_surface
+ units = K m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ kbfs ]
+ standard_name = kinematic_buoyancy_flux_at_surface
+ # m s-2 m s-1
+ units = m2 s-3
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ ustar ]
+ standard_name = surface_friction_velocity
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ obklen ]
+ standard_name = obukhov_length
+ units = m
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ s2 ]
+ standard_name = square_of_wind_shear
+ units = s-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ ri ]
+ standard_name = gradient_richardson_number
+ units = 1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ pblh ]
+ standard_name = atmosphere_boundary_layer_thickness
+ units = m
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ wstar ]
+ standard_name = atmosphere_boundary_layer_convective_velocity_scale
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ bge ]
+ standard_name = atmosphere_boundary_layer_buoyancy_gradient_enhancement_for_holtslag_boville_mixing_scheme
+ units = m s-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ kvm ]
+ standard_name = eddy_momentum_diffusivity_at_interfaces
+ units = m2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = out
+[ kvh ]
+ standard_name = eddy_heat_diffusivity_at_interfaces
+ units = m2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = out
+[ kvq ]
+ standard_name = eddy_constituent_diffusivity_at_interfaces
+ units = m2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = out
+[ cgh ]
+ standard_name = counter_gradient_coefficient_for_vertical_diffusion_of_heat_at_interfaces
+ units = J kg-1 m-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = out
+[ cgs ]
+ standard_name = counter_gradient_coefficient_for_vertical_diffusion_of_constituents_at_interfaces
+ units = s m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = out
+[ tpert ]
+ standard_name = convective_temperature_perturbation_due_to_pbl_eddies
+ units = K
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ qpert ]
+ standard_name = convective_water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water_perturbation_due_to_pbl_eddies
+ units = kg kg-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ tke ]
+ standard_name = specific_turbulent_kinetic_energy_at_interfaces
+ units = m2 s-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-table-properties]
+ name = hb_diff_free_atm_exchange_coefficients
+ type = scheme
+
+[ccpp-arg-table]
+ name = hb_diff_free_atm_exchange_coefficients_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pverp ]
+ standard_name = vertical_interface_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ s2 ]
+ standard_name = square_of_wind_shear
+ units = s-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ ri ]
+ standard_name = gradient_richardson_number
+ units = 1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ bottom_boundary ]
+ standard_name = vertical_layer_index_of_clubb_top_tbd
+ units = index
+ #REMOVECAM: change to integer after CAM snapshots are no longer used
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ kvm ]
+ standard_name = eddy_momentum_diffusivity_at_interfaces
+ units = m2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = out
+[ kvh ]
+ standard_name = eddy_heat_diffusivity_at_interfaces
+ units = m2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = out
+[ kvq ]
+ standard_name = eddy_constituent_diffusivity_at_interfaces
+ units = m2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = out
+[ cgh ]
+ standard_name = counter_gradient_coefficient_for_vertical_diffusion_of_heat_at_interfaces
+ units = J kg-1 m-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = out
+[ cgs ]
+ standard_name = counter_gradient_coefficient_for_vertical_diffusion_of_constituents_at_interfaces
+ units = s m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
diff --git a/schemes/holtslag_boville/holtslag_boville_diff_interstitials.F90 b/schemes/holtslag_boville/holtslag_boville_diff_interstitials.F90
new file mode 100644
index 00000000..dcf12bb2
--- /dev/null
+++ b/schemes/holtslag_boville/holtslag_boville_diff_interstitials.F90
@@ -0,0 +1,294 @@
+! Interstitial scheme to prepare inputs for vertical diffusion solver
+! from surface coupler fluxes and stress data
+module holtslag_boville_diff_interstitials
+ use ccpp_kinds, only: kind_phys
+
+ implicit none
+ private
+ save
+
+ ! CCPP-compliant public interfaces
+ public :: hb_diff_set_vertical_diffusion_top_init
+ public :: hb_diff_set_total_surface_stress_run
+ public :: hb_diff_prepare_vertical_diffusion_inputs_run
+ public :: hb_diff_prepare_vertical_diffusion_inputs_timestep_final
+ public :: hb_free_atm_diff_prepare_vertical_diffusion_inputs_run
+
+contains
+
+ ! Interstitial to set top of vertical diffusion.
+!> \section arg_table_hb_diff_set_vertical_diffusion_top_init Argument Table
+!! \htmlinclude hb_diff_set_vertical_diffusion_top_init.html
+ subroutine hb_diff_set_vertical_diffusion_top_init( &
+ ntop_eddy, &
+ errmsg, errflg)
+
+ ! Host model dependency
+ ! See note below. If removal of this dependency is desired in a low-top
+ ! model configuration, this call can simply be replaced with ntop_eddy = 1.
+ use ref_pres, only: press_lim_idx
+
+ ! Output arguments
+ integer, intent(out) :: ntop_eddy ! Vertical layer index of vertical diffusion top [index]
+ character(len=512), intent(out) :: errmsg ! Error message
+ integer, intent(out) :: errflg ! Error flag
+
+ ! Local parameters
+ real(kind_phys), parameter :: ntop_eddy_pres = 1.e-7_kind_phys ! Pressure below which eddy diffusion is not done in WACCM-X [Pa]
+
+ errmsg = ''
+ errflg = 0
+
+ ! Set top level for vertical diffusion.
+ ! When WACCM-X is active, the top of the model is sufficiently high that
+ ! ntop_eddy is not at level 1 (top-of-atmosphere).
+ ! In other configurations, the top of the model pressure is larger than
+ ! ntop_eddy_pres, thus ntop_eddy computes to 1.
+ ntop_eddy = press_lim_idx(ntop_eddy_pres, top=.true.)
+
+ end subroutine hb_diff_set_vertical_diffusion_top_init
+
+ ! Set total surface stresses for input into the HB PBL scheme.
+ !
+ ! NOTE: Temporarily, TMS and Beljaars are "hard-coupled" into this subroutine because
+ ! the current focus is CAM4 and TMS/Beljaars drag have not been CCPPized.
+ ! In the future, after TMS/Beljaars is implemented, this subroutine should only initialize
+ ! surface stresses from the coupler, then TMS and Beljaars can work on adding the respective
+ ! surface stresses to the tautotx/tautoty used by the PBL scheme.
+!> \section arg_table_hb_diff_set_total_surface_stress_run Argument Table
+!! \htmlinclude hb_diff_set_total_surface_stress_run.html
+ subroutine hb_diff_set_total_surface_stress_run( &
+ ncol, &
+ wsx_from_coupler, wsy_from_coupler, &
+ tautmsx, tautmsy, &
+ taubljx, taubljy, &
+ tautotx, tautoty, &
+ errmsg, errflg)
+
+ ! Input arguments
+ integer, intent(in) :: ncol
+ real(kind_phys), intent(in) :: wsx_from_coupler(:) ! Surface eastward wind stress from coupler [Pa]
+ real(kind_phys), intent(in) :: wsy_from_coupler(:) ! Surface northward wind stress from coupler [Pa]
+ real(kind_phys), intent(in) :: tautmsx(:) ! Eastward turbulent mountain surface stress [Pa]
+ real(kind_phys), intent(in) :: tautmsy(:) ! Northward turbulent mountain surface stress [Pa]
+ real(kind_phys), intent(in) :: taubljx(:) ! Eastward Beljaars surface stress [Pa]
+ real(kind_phys), intent(in) :: taubljy(:) ! Northward Beljaars surface stress [Pa]
+
+ ! Output arguments
+ real(kind_phys), intent(out) :: tautotx(:) ! Eastward total stress at surface for boundary layer scheme [Pa]
+ real(kind_phys), intent(out) :: tautoty(:) ! Northward total stress at surface for boundary layer scheme [Pa]
+ character(len=512), intent(out) :: errmsg
+ integer, intent(out) :: errflg
+
+ errmsg = ''
+ errflg = 0
+
+ ! Initialize total surface stresses from coupler
+ ! These are used for HB diffusion scheme and later PBL diagnostics but
+ ! not for the vertical diffusion solver, which uses surface stresses from the coupler
+ ! or just zero (in the case of CLUBB)
+ tautotx(:ncol) = wsx_from_coupler(:ncol)
+ tautoty(:ncol) = wsy_from_coupler(:ncol)
+
+ ! Add turbulent mountain stress to total surface stress
+ tautotx(:ncol) = tautotx(:ncol) + tautmsx(:ncol)
+ tautoty(:ncol) = tautoty(:ncol) + tautmsy(:ncol)
+
+ ! Add Beljaars integrated drag to total surface stress
+ tautotx(:ncol) = tautotx(:ncol) + taubljx(:ncol)
+ tautoty(:ncol) = tautoty(:ncol) + taubljy(:ncol)
+
+ end subroutine hb_diff_set_total_surface_stress_run
+
+ ! Interstitial for full HB (CAM4) to handle inputs from coupler and pass them
+ ! to vertical diffusion solver, as well as vertical coordinate set up.
+!> \section arg_table_hb_diff_prepare_vertical_diffusion_inputs_run Argument Table
+!! \htmlinclude hb_diff_prepare_vertical_diffusion_inputs_run.html
+ subroutine hb_diff_prepare_vertical_diffusion_inputs_run( &
+ ncol, pverp, pcnst, &
+ const_props, &
+ wsx_from_coupler, wsy_from_coupler, &
+ shf_from_coupler, &
+ cflx_from_coupler, &
+ pint, &
+ ! below output
+ taux, tauy, &
+ shflux, &
+ cflux, &
+ itaures, &
+ p, &
+ q_wv_cflx, &
+ errmsg, errflg)
+
+ use coords_1d, only: Coords1D
+
+ ! framework dependency for const_props
+ use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
+
+ ! dependency to get constituent index
+ use ccpp_const_utils, only: ccpp_const_get_idx
+
+ ! Input arguments
+ integer, intent(in) :: ncol ! Number of atmospheric columns [count]
+ integer, intent(in) :: pverp ! Number of vertical interfaces [count]
+ integer, intent(in) :: pcnst ! Number of CCPP constituents [count]
+ type(ccpp_constituent_prop_ptr_t), &
+ intent(in) :: const_props(:) ! CCPP constituent properties pointer
+ real(kind_phys), intent(in) :: wsx_from_coupler(:) ! Surface eastward wind stress from coupler [Pa]
+ real(kind_phys), intent(in) :: wsy_from_coupler(:) ! Surface northward wind stress from coupler [Pa]
+ real(kind_phys), intent(in) :: shf_from_coupler(:) ! Surface upward sensible heat flux from coupler [W m-2]
+ real(kind_phys), intent(in) :: cflx_from_coupler(:,:) ! Surface upward constituent fluxes from coupler [kg m-2 s-1]
+ real(kind_phys), intent(in) :: pint(:,:) ! Air pressure at interfaces [Pa]
+
+ ! Output arguments
+ real(kind_phys), intent(out) :: taux(:) ! Eastward stress at surface for vertical diffusion [Pa]
+ real(kind_phys), intent(out) :: tauy(:) ! Northward stress at surface for vertical diffusion [Pa]
+ real(kind_phys), intent(out) :: shflux(:) ! Surface upward sensible heat flux for vertical diffusion [W m-2]
+ real(kind_phys), intent(out) :: cflux(:,:) ! Surface upward constituent fluxes for vertical diffusion [kg m-2 s-1]
+ logical, intent(out) :: itaures ! Flag for updating residual stress at surface in vertical diffusion [flag]
+ type(coords1d), intent(out) :: p ! Vertical moist pressure coordinates for vertical diffusion [Pa]
+ real(kind_phys), intent(out) :: q_wv_cflx(:) ! Surface upward water vapor flux [kg m-2 s-1]
+ character(len=512), intent(out) :: errmsg ! Error message
+ integer, intent(out) :: errflg ! Error flag
+
+ ! Local variables
+ integer :: const_wv_idx ! Water vapor constituent index
+
+ errmsg = ''
+ errflg = 0
+
+ ! Copy surface fluxes from coupler to vertical diffusion input arrays
+ taux(:ncol) = wsx_from_coupler(:ncol)
+ tauy(:ncol) = wsy_from_coupler(:ncol)
+ shflux(:ncol) = shf_from_coupler(:ncol)
+ cflux(:ncol,:pcnst) = cflx_from_coupler(:ncol,:pcnst)
+
+ ! Set flag for updating residual stress to true
+ itaures = .true.
+
+ ! Initialize pressure coordinate object for vertical diffusion solver
+ p = Coords1D(pint(:ncol,:pverp))
+
+ ! Get water vapor constituent index
+ call ccpp_const_get_idx(const_props, &
+ 'water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water', &
+ const_wv_idx, errmsg, errflg)
+ if (errflg /= 0) return
+
+ ! Extract water vapor flux for use in HB.
+ q_wv_cflx(:ncol) = cflx_from_coupler(:ncol, const_wv_idx)
+
+ end subroutine hb_diff_prepare_vertical_diffusion_inputs_run
+
+ ! Interstitial to clean up vertical coordinate after use.
+!> \section arg_table_hb_diff_prepare_vertical_diffusion_inputs_timestep_final Argument Table
+!! \htmlinclude hb_diff_prepare_vertical_diffusion_inputs_timestep_final.html
+ subroutine hb_diff_prepare_vertical_diffusion_inputs_timestep_final(p, errmsg, errflg)
+ use coords_1d, only: Coords1D
+
+ type(coords1d), intent(inout) :: p ! Vertical moist pressure coordinates for vertical diffusion [Pa]
+ character(len=512), intent(out) :: errmsg ! Error message
+ integer, intent(out) :: errflg ! Error flag
+
+ errmsg = ''
+ errflg = 0
+
+ call p%finalize()
+ end subroutine hb_diff_prepare_vertical_diffusion_inputs_timestep_final
+
+ ! Interstitial for free atmosphere version of HB used above CLUBB which will allow
+ ! the diffusion solver to handle non-water vapor surface fluxes (CAM6)
+ ! or no surface fluxes (CAM7)
+!> \section arg_table_hb_free_atm_diff_prepare_vertical_diffusion_inputs_run Argument Table
+!! \htmlinclude hb_free_atm_diff_prepare_vertical_diffusion_inputs_run.html
+ subroutine hb_free_atm_diff_prepare_vertical_diffusion_inputs_run( &
+ ncol, pverp, pcnst, &
+ const_props, &
+ apply_nonwv_cflx, &
+ cflx_from_coupler, &
+ pint, &
+ ! below output
+ taux, tauy, &
+ shflux, &
+ cflux, &
+ itaures, &
+ p, &
+ q_wv_cflx, &
+ errmsg, errflg)
+
+ use coords_1d, only: Coords1D
+
+ ! framework dependency for const_props
+ use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
+
+ ! dependency to get constituent index
+ use ccpp_const_utils, only: ccpp_const_get_idx
+
+ ! Input arguments
+ integer, intent(in) :: ncol ! Number of atmospheric columns [count]
+ integer, intent(in) :: pverp ! Number of vertical interfaces [count]
+ integer, intent(in) :: pcnst ! Number of CCPP constituents [count]
+ type(ccpp_constituent_prop_ptr_t), &
+ intent(in) :: const_props(:) ! CCPP constituent properties pointer
+ logical, intent(in) :: apply_nonwv_cflx ! Flag for applying constituent fluxes excluding water vapor [flag]
+ real(kind_phys), intent(in) :: cflx_from_coupler(:,:) ! Surface upward constituent fluxes from coupler [kg m-2 s-1]
+ real(kind_phys), intent(in) :: pint(:,:) ! Air pressure at interfaces [Pa]
+
+ ! Output arguments
+ real(kind_phys), intent(out) :: taux(:) ! Eastward stress at surface for vertical diffusion [Pa]
+ real(kind_phys), intent(out) :: tauy(:) ! Northward stress at surface for vertical diffusion [Pa]
+ real(kind_phys), intent(out) :: shflux(:) ! Surface upward sensible heat flux for vertical diffusion [W m-2]
+ real(kind_phys), intent(out) :: cflux(:,:) ! Surface upward constituent fluxes for vertical diffusion [kg m-2 s-1]
+ logical, intent(out) :: itaures ! Flag for updating residual stress at surface in vertical diffusion [flag]
+ type(coords1d), intent(out) :: p ! Vertical moist pressure coordinates for vertical diffusion [Pa]
+ real(kind_phys), intent(out) :: q_wv_cflx(:) ! Surface upward water vapor flux (for PBL scheme) [kg m-2 s-1]
+ character(len=512), intent(out) :: errmsg ! Error message
+ integer, intent(out) :: errflg ! Error flag
+
+ ! Local variables
+ integer :: const_wv_idx ! Water vapor constituent index
+
+ errmsg = ''
+ errflg = 0
+
+ ! Check constituents list and locate water vapor index
+ ! (not assumed to be 1)
+ call ccpp_const_get_idx(const_props, &
+ 'water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water', &
+ const_wv_idx, errmsg, errflg)
+ if (errflg /= 0) return
+
+ ! CLUBB applies some fluxes itself, but we still want constituent
+ ! fluxes applied here (except water vapor).
+ taux(:ncol) = 0._kind_phys
+ tauy(:ncol) = 0._kind_phys
+ shflux(:ncol) = 0._kind_phys
+
+ ! If apply_nonwv_cflx is on, then non-water vapor constituent fluxes from the coupler
+ ! are applied using vertical diffusion.
+ if (apply_nonwv_cflx) then
+ ! Copy non-water vapor constituent fluxes from coupler
+ cflux(:ncol, :) = cflx_from_coupler(:ncol, :)
+ ! But still zero out water vapor flux
+ cflux(:ncol, const_wv_idx) = 0._kind_phys
+ else
+ ! Surface fluxes applied in CLUBB emissions module (CAM7)
+ ! so vertical diffusion applies no fluxes.
+ cflux(:ncol, :) = 0._kind_phys
+ end if
+
+ ! Set flag for updating residual stress to true
+ itaures = .true.
+
+ ! Initialize pressure coordinate object for vertical diffusion solver
+ p = Coords1D(pint(:ncol,:pverp))
+
+ ! Extract water vapor flux for use in the HB scheme to calculate
+ ! kinematic water vapor fluxes.
+ ! This is separate from the cflux above, which is provided to the diffusion
+ ! solver for flux application.
+ q_wv_cflx(:ncol) = cflx_from_coupler(:ncol, const_wv_idx)
+
+ end subroutine hb_free_atm_diff_prepare_vertical_diffusion_inputs_run
+
+end module holtslag_boville_diff_interstitials
diff --git a/schemes/holtslag_boville/holtslag_boville_diff_interstitials.meta b/schemes/holtslag_boville/holtslag_boville_diff_interstitials.meta
new file mode 100644
index 00000000..d9eeb452
--- /dev/null
+++ b/schemes/holtslag_boville/holtslag_boville_diff_interstitials.meta
@@ -0,0 +1,344 @@
+[ccpp-table-properties]
+ name = hb_diff_set_vertical_diffusion_top
+ type = scheme
+ dependencies = ../../../../data/ref_pres.F90
+
+[ccpp-arg-table]
+ name = hb_diff_set_vertical_diffusion_top_init
+ type = scheme
+[ ntop_eddy ]
+ standard_name = vertical_layer_index_of_eddy_vertical_diffusion_top
+ units = index
+ dimensions = ()
+ type = integer
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ dimensions = ()
+ type = character | kind = len=512
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ dimensions = ()
+ type = integer
+ intent = out
+
+[ccpp-table-properties]
+ name = hb_diff_set_total_surface_stress
+ type = scheme
+
+[ccpp-arg-table]
+ name = hb_diff_set_total_surface_stress_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+[ wsx_from_coupler ]
+ standard_name = surface_eastward_wind_stress_from_coupler
+ units = N m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = in
+[ wsy_from_coupler ]
+ standard_name = surface_northward_wind_stress_from_coupler
+ units = N m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = in
+[ tautmsx ]
+ standard_name = eastward_turbulent_mountain_surface_stress_tbd
+ units = N m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = in
+[ tautmsy ]
+ standard_name = northward_turbulent_mountain_surface_stress_tbd
+ units = N m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = in
+[ taubljx ]
+ standard_name = eastward_beljaars_surface_stress_tbd
+ units = N m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = in
+[ taubljy ]
+ standard_name = northward_beljaars_surface_stress_tbd
+ units = N m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = in
+[ tautotx ]
+ standard_name = eastward_total_stress_at_surface_for_holtslag_boville_boundary_layer_scheme
+ units = N m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = out
+[ tautoty ]
+ standard_name = northward_total_stress_at_surface_for_holtslag_boville_boundary_layer_scheme
+ units = N m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ dimensions = ()
+ type = character | kind = len=512
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ dimensions = ()
+ type = integer
+ intent = out
+
+[ccpp-table-properties]
+ name = hb_diff_prepare_vertical_diffusion_inputs
+ type = scheme
+ dependencies = ../../to_be_ccppized/coords_1d.F90
+
+[ccpp-arg-table]
+ name = hb_diff_prepare_vertical_diffusion_inputs_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+[ pverp ]
+ standard_name = vertical_interface_dimension
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+[ pcnst ]
+ standard_name = number_of_ccpp_constituents
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+[ const_props ]
+ standard_name = ccpp_constituent_properties
+ units = none
+ type = ccpp_constituent_prop_ptr_t
+ dimensions = (number_of_ccpp_constituents)
+ intent = in
+[ wsx_from_coupler ]
+ standard_name = surface_eastward_wind_stress_from_coupler
+ units = N m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = in
+[ wsy_from_coupler ]
+ standard_name = surface_northward_wind_stress_from_coupler
+ units = N m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = in
+[ shf_from_coupler ]
+ standard_name = surface_upward_sensible_heat_flux_from_coupler
+ units = W m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = in
+[ cflx_from_coupler ]
+ standard_name = surface_upward_ccpp_constituent_fluxes_from_coupler
+ units = kg m-2 s-1
+ dimensions = (horizontal_loop_extent, number_of_ccpp_constituents)
+ type = real | kind = kind_phys
+ intent = in
+[ pint ]
+ standard_name = air_pressure_at_interface
+ units = Pa
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ type = real | kind = kind_phys
+ intent = in
+[ taux ]
+ standard_name = eastward_stress_at_surface_for_vertical_diffusion
+ units = N m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = out
+[ tauy ]
+ standard_name = northward_stress_at_surface_for_vertical_diffusion
+ units = N m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = out
+[ shflux ]
+ standard_name = surface_upward_sensible_heat_flux_for_vertical_diffusion
+ units = W m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = out
+[ cflux ]
+ standard_name = surface_upward_ccpp_constituent_fluxes_for_vertical_diffusion
+ units = kg m-2 s-1
+ dimensions = (horizontal_loop_extent, number_of_ccpp_constituents)
+ type = real | kind = kind_phys
+ intent = out
+[ itaures ]
+ standard_name = update_residual_stress_at_surface_in_vertical_diffusion
+ units = flag
+ dimensions = ()
+ type = logical
+ intent = out
+[ p ]
+ standard_name = vertical_moist_pressure_coordinates_for_vertical_diffusion
+ units = none
+ dimensions = ()
+ type = coords1d
+ intent = out
+[ q_wv_cflx ]
+ standard_name = surface_upward_water_vapor_flux
+ units = kg m-2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ dimensions = ()
+ type = character | kind = len=512
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ dimensions = ()
+ type = integer
+ intent = out
+
+[ccpp-arg-table]
+ name = hb_diff_prepare_vertical_diffusion_inputs_timestep_final
+ type = scheme
+[ p ]
+ standard_name = vertical_moist_pressure_coordinates_for_vertical_diffusion
+ units = none
+ dimensions = ()
+ type = coords1d
+ intent = inout
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ dimensions = ()
+ type = character | kind = len=512
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ dimensions = ()
+ type = integer
+ intent = out
+
+[ccpp-table-properties]
+ name = hb_free_atm_diff_prepare_vertical_diffusion_inputs
+ type = scheme
+ dependencies = ../../to_be_ccppized/coords_1d.F90
+
+[ccpp-arg-table]
+ name = hb_free_atm_diff_prepare_vertical_diffusion_inputs_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+[ pverp ]
+ standard_name = vertical_interface_dimension
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+[ pcnst ]
+ standard_name = number_of_ccpp_constituents
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+[ const_props ]
+ standard_name = ccpp_constituent_properties
+ units = none
+ type = ccpp_constituent_prop_ptr_t
+ dimensions = (number_of_ccpp_constituents)
+ intent = in
+[ apply_nonwv_cflx ]
+ standard_name = do_apply_ccpp_constituent_fluxes_excluding_water_vapor_in_vertical_diffusion
+ units = flag
+ dimensions = ()
+ type = logical
+ intent = in
+[ cflx_from_coupler ]
+ standard_name = surface_upward_ccpp_constituent_fluxes_from_coupler
+ units = kg m-2 s-1
+ dimensions = (horizontal_loop_extent, number_of_ccpp_constituents)
+ type = real | kind = kind_phys
+ intent = in
+[ pint ]
+ standard_name = air_pressure_at_interface
+ units = Pa
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ type = real | kind = kind_phys
+ intent = in
+[ taux ]
+ standard_name = eastward_stress_at_surface_for_vertical_diffusion
+ units = N m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = out
+[ tauy ]
+ standard_name = northward_stress_at_surface_for_vertical_diffusion
+ units = N m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = out
+[ shflux ]
+ standard_name = surface_upward_sensible_heat_flux_for_vertical_diffusion
+ units = W m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = out
+[ cflux ]
+ standard_name = surface_upward_ccpp_constituent_fluxes_for_vertical_diffusion
+ units = kg m-2 s-1
+ dimensions = (horizontal_loop_extent, number_of_ccpp_constituents)
+ type = real | kind = kind_phys
+ intent = out
+[ itaures ]
+ standard_name = update_residual_stress_at_surface_in_vertical_diffusion
+ units = flag
+ dimensions = ()
+ type = logical
+ intent = out
+[ p ]
+ standard_name = vertical_moist_pressure_coordinates_for_vertical_diffusion
+ units = none
+ dimensions = ()
+ type = coords1d
+ intent = out
+[ q_wv_cflx ]
+ standard_name = surface_upward_water_vapor_flux
+ units = kg m-2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ dimensions = ()
+ type = character | kind = len=512
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ dimensions = ()
+ type = integer
+ intent = out
diff --git a/schemes/holtslag_boville/holtslag_boville_diff_options.F90 b/schemes/holtslag_boville/holtslag_boville_diff_options.F90
new file mode 100644
index 00000000..38f7ecff
--- /dev/null
+++ b/schemes/holtslag_boville/holtslag_boville_diff_options.F90
@@ -0,0 +1,43 @@
+! Module to load namelist options for the Holstlag-Boville boundary layer scheme
+module holtslag_boville_diff_options
+ use ccpp_kinds, only: kind_phys
+
+ implicit none
+ private
+ save
+
+ ! public CCPP-compliant subroutines
+ public :: holtslag_boville_diff_options_init
+
+contains
+
+!> \section arg_table_holtslag_boville_diff_options_init Argument Table
+!! \htmlinclude arg_table_holtslag_boville_diff_options_init.html
+ subroutine holtslag_boville_diff_options_init( &
+ amIRoot, iulog, &
+ is_hbr_pbl_scheme, &
+ errmsg, errflg)
+
+ ! Input arguments
+ logical, intent(in) :: amIRoot ! are we on the MPI root task?
+ integer, intent(in) :: iulog ! log output unit
+ logical, intent(in) :: is_hbr_pbl_scheme ! is HBR = true; is HB = false [flag]
+
+ ! Output arguments
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ errmsg = ''
+ errflg = 0
+
+ if(amIRoot) then
+ if(is_hbr_pbl_scheme) then
+ write(iulog,*) 'Holtslag-Boville PBL: initializing as HBR PBL scheme.'
+ else
+ write(iulog,*) 'Holtslag-Boville PBL: initializing as HB PBL scheme.'
+ endif
+ endif
+
+ end subroutine holtslag_boville_diff_options_init
+
+end module holtslag_boville_diff_options
diff --git a/schemes/holtslag_boville/holtslag_boville_diff_options.meta b/schemes/holtslag_boville/holtslag_boville_diff_options.meta
new file mode 100644
index 00000000..dc31cb7f
--- /dev/null
+++ b/schemes/holtslag_boville/holtslag_boville_diff_options.meta
@@ -0,0 +1,37 @@
+[ccpp-table-properties]
+ name = holtslag_boville_diff_options
+ type = scheme
+
+[ccpp-arg-table]
+ name = holtslag_boville_diff_options_init
+ type = scheme
+[ amIRoot ]
+ standard_name = flag_for_mpi_root
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = in
+[ iulog ]
+ standard_name = log_output_unit
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = in
+[ is_hbr_pbl_scheme ]
+ standard_name = flag_for_hbr_configuration_of_holtslag_boville_mixing_scheme
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = in
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
diff --git a/schemes/holtslag_boville/holtslag_boville_diff_options_namelist.xml b/schemes/holtslag_boville/holtslag_boville_diff_options_namelist.xml
new file mode 100644
index 00000000..4db4a597
--- /dev/null
+++ b/schemes/holtslag_boville/holtslag_boville_diff_options_namelist.xml
@@ -0,0 +1,92 @@
+
+
+
+
+
+
+
+
+
+ logical
+ pbl
+ hb_nl
+ flag_for_hbr_configuration_of_holtslag_boville_mixing_scheme
+ flag
+
+ Flag to use the Rasch modified version of the Holtslag-Boville PBL scheme.
+
+
+ .false.
+
+
+
diff --git a/schemes/sima_diagnostics/diffusion_solver_diagnostics.F90 b/schemes/sima_diagnostics/diffusion_solver_diagnostics.F90
new file mode 100644
index 00000000..54c3f375
--- /dev/null
+++ b/schemes/sima_diagnostics/diffusion_solver_diagnostics.F90
@@ -0,0 +1,445 @@
+! Diagnostic scheme for vertical diffusion solver
+module diffusion_solver_diagnostics
+ use ccpp_kinds, only: kind_phys
+
+ implicit none
+ private
+ save
+
+ ! public CCPP-compliant subroutines
+ public :: vertical_diffusion_tendencies_diagnostics_init
+ public :: vertical_diffusion_tendencies_diagnostics_run
+
+contains
+
+!> \section arg_table_vertical_diffusion_tendencies_diagnostics_init Argument Table
+!! \htmlinclude vertical_diffusion_tendencies_diagnostics_init.html
+ subroutine vertical_diffusion_tendencies_diagnostics_init(const_props, errmsg, errflg)
+
+ ! framework dependency for const_props
+ use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
+
+ ! dependency to get constituent index
+ use ccpp_const_utils, only: ccpp_const_get_idx
+
+ use cam_history, only: history_add_field
+ use cam_history_support, only: horiz_only
+
+ type(ccpp_constituent_prop_ptr_t), intent(in) :: const_props(:)
+ character(len=512), intent(out) :: errmsg
+ integer, intent(out) :: errflg
+
+ errmsg = ''
+ errflg = 0
+
+ ! History add field calls
+ call history_add_field('QT', 'Total water mixing ratio', 'lev', 'inst', 'kg/kg')
+ call history_add_field('SL', 'Liquid water static energy', 'lev', 'inst', 'J/kg')
+ call history_add_field('SLV', 'Liquid water virtual static energy', 'lev', 'inst', 'J/kg')
+ call history_add_field('SLFLX', 'Liquid static energy flux', 'ilev', 'inst', 'W/m2')
+ call history_add_field('QTFLX', 'Total water flux', 'ilev', 'inst', 'kg/m2/s')
+ call history_add_field('UFLX', 'Zonal momentum flux', 'ilev', 'inst', 'N/m2')
+ call history_add_field('VFLX', 'Meridional momentum flux', 'ilev', 'inst', 'N/m2')
+
+ ! Pre-PBL profiles
+ call history_add_field('qt_pre_PBL', 'Total water mixing ratio before PBL', 'lev', 'inst', 'kg/kg')
+ call history_add_field('sl_pre_PBL', 'Liquid water static energy before PBL', 'lev', 'inst', 'J/kg')
+ call history_add_field('slv_pre_PBL', 'Liquid water virtual static energy before PBL', 'lev', 'inst', 'J/kg')
+ call history_add_field('u_pre_PBL', 'Zonal wind before PBL', 'lev', 'inst', 'm/s')
+ call history_add_field('v_pre_PBL', 'Meridional wind before PBL', 'lev', 'inst', 'm/s')
+ call history_add_field('qv_pre_PBL', 'Water vapor mixing ratio before PBL', 'lev', 'inst', 'kg/kg')
+ call history_add_field('ql_pre_PBL', 'Cloud liquid water mixing ratio before PBL', 'lev', 'inst', 'kg/kg')
+ call history_add_field('qi_pre_PBL', 'Cloud ice water mixing ratio before PBL', 'lev', 'inst', 'kg/kg')
+ call history_add_field('t_pre_PBL', 'Temperature before PBL', 'lev', 'inst', 'K')
+ call history_add_field('rh_pre_PBL', 'Relative humidity before PBL', 'lev', 'inst', '%')
+
+ ! Post-PBL profiles
+ call history_add_field('qt_aft_PBL', 'Total water mixing ratio after PBL', 'lev', 'inst', 'kg/kg')
+ call history_add_field('sl_aft_PBL', 'Liquid water static energy after PBL', 'lev', 'inst', 'J/kg')
+ call history_add_field('slv_aft_PBL', 'Liquid water virtual static energy after PBL', 'lev', 'inst', 'J/kg')
+ call history_add_field('u_aft_PBL', 'Zonal wind after PBL', 'lev', 'inst', 'm/s')
+ call history_add_field('v_aft_PBL', 'Meridional wind after PBL', 'lev', 'inst', 'm/s')
+ call history_add_field('qv_aft_PBL', 'Water vapor mixing ratio after PBL', 'lev', 'inst', 'kg/kg')
+ call history_add_field('ql_aft_PBL', 'Cloud liquid water mixing ratio after PBL', 'lev', 'inst', 'kg/kg')
+ call history_add_field('qi_aft_PBL', 'Cloud ice water mixing ratio after PBL', 'lev', 'inst', 'kg/kg')
+ call history_add_field('t_aft_PBL', 'Temperature after PBL', 'lev', 'inst', 'K')
+ call history_add_field('rh_aft_PBL', 'Relative humidity after PBL', 'lev', 'inst', '%')
+
+ ! PBL fluxes
+ call history_add_field('slflx_PBL', 'Liquid static energy flux by PBL', 'ilev', 'inst', 'W/m2')
+ call history_add_field('qtflx_PBL', 'Total water flux by PBL', 'ilev', 'inst', 'kg/m2/s')
+ call history_add_field('uflx_PBL', 'Zonal momentum flux by PBL', 'ilev', 'inst', 'N/m2')
+ call history_add_field('vflx_PBL', 'Meridional momentum flux by PBL', 'ilev', 'inst', 'N/m2')
+
+ ! Counter-gradient fluxes
+ call history_add_field('slflx_cg_PBL','Counter-gradient liquid static energy flux by PBL', 'ilev', 'inst', 'W/m2')
+ call history_add_field('qtflx_cg_PBL','Counter-gradient total water flux by PBL', 'ilev', 'inst', 'kg/m2/s')
+ call history_add_field('uflx_cg_PBL', 'Counter-gradient zonal momentum flux by PBL', 'ilev', 'inst', 'N/m2')
+ call history_add_field('vflx_cg_PBL', 'Counter-gradient meridional momentum flux by PBL', 'ilev', 'inst', 'N/m2')
+
+ ! PBL tendencies
+ call history_add_field('qtten_PBL', 'Total water mixing ratio tendency by PBL', 'lev', 'inst', 'kg/kg/s')
+ call history_add_field('slten_PBL', 'Liquid static energy tendency by PBL', 'lev', 'inst', 'J/kg/s')
+ call history_add_field('uten_PBL', 'Zonal wind tendency by PBL', 'lev', 'inst', 'm/s2')
+ call history_add_field('vten_PBL', 'Meridional wind tendency by PBL', 'lev', 'inst', 'm/s2')
+ call history_add_field('qvten_PBL', 'Water vapor tendency by PBL', 'lev', 'inst', 'kg/kg/s')
+ call history_add_field('qlten_PBL', 'Cloud liquid water tendency by PBL', 'lev', 'inst', 'kg/kg/s')
+ call history_add_field('qiten_PBL', 'Cloud ice water tendency by PBL', 'lev', 'inst', 'kg/kg/s')
+ call history_add_field('tten_PBL', 'Temperature tendency by PBL', 'lev', 'inst', 'K/s')
+ call history_add_field('rhten_PBL', 'Relative humidity tendency by PBL', 'lev', 'inst', '%/s')
+
+ ! Vertical diffusion tendencies
+ call history_add_field('DTVKE', 'dT/dt vertical diffusion KE dissipation', 'lev', 'inst', 'K s-1')
+ call history_add_field('DTV', 'Temperature vertical diffusion', 'lev', 'inst', 'K s-1')
+ call history_add_field('DUV', 'Zonal wind vertical diffusion', 'lev', 'inst', 'm s-2')
+ call history_add_field('DVV', 'Meridional wind vertical diffusion', 'lev', 'inst', 'm s-2')
+
+ ! These fields are constituent dependent.
+ ! TODO: use the constituent properties pointer to initialize these
+ ! diagnostics for all constituents. The "short name" of the constituents
+ ! is not currently available, so we initialize some manually here.
+ call history_add_field('VERT_DIFF_Q', 'Vertical diffusion of water vapor', 'lev', 'inst', 'kg/kg/s')
+ call history_add_field('VERT_DIFF_CLDLIQ', 'Vertical diffusion of cloud liquid water', 'lev', 'inst', 'kg/kg/s')
+ call history_add_field('VERT_DIFF_CLDICE', 'Vertical diffusion of cloud ice water', 'lev', 'inst', 'kg/kg/s')
+
+ end subroutine vertical_diffusion_tendencies_diagnostics_init
+
+!> \section arg_table_vertical_diffusion_tendencies_diagnostics_run Argument Table
+!! \htmlinclude vertical_diffusion_tendencies_diagnostics_run.html
+ subroutine vertical_diffusion_tendencies_diagnostics_run( &
+ ncol, pver, pverp, dt, &
+ const_props, &
+ latvap, latice, zvir, cpair, gravit, rair, &
+ pmid, pint, zi, zm, &
+ kvh, kvm, cgh, cgs, &
+ cam_in_shf, cam_in_cflx, &
+ tautotx, tautoty, &
+ t0, &
+ q0, s0, u0, v0, &
+ q1, s1, u1, v1, &
+ dtk, &
+ tend_s, tend_u, tend_v, tend_q, &
+ errmsg, errflg)
+
+ use cam_history, only: history_out_field
+ use runtime_obj, only: wv_stdname
+
+ ! framework dependency for const_props
+ use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
+
+ ! dependency to get constituent index
+ use ccpp_const_utils, only: ccpp_const_get_idx
+
+ ! utility subroutine for calculating diagnostics based on outputs
+ ! from the vertical diffusion solver.
+ use vertical_diffusion_diagnostic_utils, only: vertical_diffusion_diagnostic_profiles
+
+ ! Input arguments
+ integer, intent(in) :: ncol
+ integer, intent(in) :: pver
+ integer, intent(in) :: pverp
+ real(kind_phys), intent(in) :: dt
+ type(ccpp_constituent_prop_ptr_t), intent(in) :: const_props(:)
+ real(kind_phys), intent(in) :: latvap ! Latent heat of vaporization [J kg-1]
+ real(kind_phys), intent(in) :: latice ! Latent heat of fusion [J kg-1]
+ real(kind_phys), intent(in) :: zvir ! rh2o/rair - 1 [1]
+ real(kind_phys), intent(in) :: cpair ! Specific heat of dry air at constant pressure [J kg-1 K-1]
+ real(kind_phys), intent(in) :: gravit ! Acceleration due to gravity [m s-2]
+ real(kind_phys), intent(in) :: rair ! Gas constant for dry air [J kg-1 K-1]
+ real(kind_phys), intent(in) :: pmid(:, :) ! Midpoint pressure [Pa]
+ real(kind_phys), intent(in) :: pint(:, :) ! Interface pressure [Pa]
+ real(kind_phys), intent(in) :: zi(:, :) ! Geopotential height at interfaces [m]
+ real(kind_phys), intent(in) :: zm(:, :) ! Geopotential height at midpoints [m]
+
+ ! Diffusion coefficients and PBL diagnostics
+ real(kind_phys), intent(in) :: kvh(:, :) ! Eddy diffusivity for heat at interfaces [m2 s-1]
+ real(kind_phys), intent(in) :: kvm(:, :) ! Eddy diffusivity for momentum at interfaces [m2 s-1]
+ real(kind_phys), intent(in) :: cgh(:, :) ! Counter-gradient term for heat [K m s-1]
+ real(kind_phys), intent(in) :: cgs(:, :) ! Counter-gradient star [s m-2]
+
+ ! Surface fluxes from coupler
+ real(kind_phys), intent(in) :: cam_in_shf(:) ! Surface sensible heat flux [W m-2]
+ real(kind_phys), intent(in) :: cam_in_cflx(:,:) ! Surface constituent fluxes [kg m-2 s-1]
+
+ ! Surface stresses
+ real(kind_phys), intent(in) :: tautotx(:) ! Total zonal surface stress [N m-2]
+ real(kind_phys), intent(in) :: tautoty(:) ! Total meridional surface stress [N m-2]
+
+ ! Initial state (before vertical diffusion)
+ real(kind_phys), intent(in) :: t0(:, :) ! Temperature before diffusion [K]
+ real(kind_phys), intent(in) :: q0(:, :, :) ! Constituent mixing ratios before diffusion [kg kg-1]
+ real(kind_phys), intent(in) :: s0(:, :) ! Dry static energy before diffusion [J kg-1]
+ real(kind_phys), intent(in) :: u0(:, :) ! Zonal wind before diffusion [m s-1]
+ real(kind_phys), intent(in) :: v0(:, :) ! Meridional wind before diffusion [m s-1]
+
+ ! "After state" (provisionally updated) quantities from diffusion solver
+ real(kind_phys), intent(in) :: q1(:, :, :) ! Constituent mixing ratios after diffusion [J kg-1]
+ real(kind_phys), intent(in) :: s1(:, :) ! Dry static energy after diffusion [J kg-1]
+ real(kind_phys), intent(in) :: u1(:, :) ! Zonal wind after diffusion [m s-1]
+ real(kind_phys), intent(in) :: v1(:, :) ! Meridional wind after diffusion [m s-1]
+
+ ! Output from vertical diffusion kinetic energy (KE) dissipation
+ ! (in vertical_diffusion_diffuse_horizontal_momentum_run)
+ real(kind_phys), intent(in) :: dtk(:, :) ! T tendency from KE dissipation [J kg-1]
+
+ ! Tendencies from vertical diffusion
+ real(kind_phys), intent(in) :: tend_s(:, :) ! Dry static energy tendency [J kg-1 s-1]
+ real(kind_phys), intent(in) :: tend_u(:, :) ! Zonal wind tendency [m s-2]
+ real(kind_phys), intent(in) :: tend_v(:, :) ! Meridional wind tendency [m s-2]
+ real(kind_phys), intent(in) :: tend_q(:, :, :) ! Constituent tendency from vertical diffusion [kg kg-1 s-1]
+
+ ! Error handling
+ character(len=512), intent(out) :: errmsg
+ integer, intent(out) :: errflg
+
+ ! Local variables for constituent indices
+ integer :: const_wv_idx, const_cldliq_idx, const_cldice_idx
+
+ ! Local variables for diagnostic profiles
+ real(kind_phys) :: qt_pre_PBL(ncol, pver) ! Total water mixing ratio before PBL [kg kg-1]
+ real(kind_phys) :: sl_pre_PBL(ncol, pver) ! Liquid water static energy before PBL [J kg-1]
+ real(kind_phys) :: slv_pre_PBL(ncol, pver) ! Virtual liquid water static energy before PBL [J kg-1]
+ real(kind_phys) :: rh_pre_PBL(ncol, pver) ! Relative humidity before PBL [%]
+ real(kind_phys) :: qt_aft_PBL(ncol, pver) ! Total water mixing ratio after PBL [kg kg-1]
+ real(kind_phys) :: sl_aft_PBL(ncol, pver) ! Liquid water static energy after PBL [J kg-1]
+ real(kind_phys) :: slv_aft_PBL(ncol, pver) ! Virtual liquid water static energy after PBL [J kg-1]
+ real(kind_phys) :: qv_aft_PBL(ncol, pver) ! Water vapor mixing ratio after PBL [kg kg-1]
+ real(kind_phys) :: ql_aft_PBL(ncol, pver) ! Cloud liquid water mixing ratio after PBL [kg kg-1]
+ real(kind_phys) :: qi_aft_PBL(ncol, pver) ! Cloud ice water mixing ratio after PBL [kg kg-1]
+ real(kind_phys) :: t_aft_PBL(ncol, pver) ! Temperature after PBL [K]
+ real(kind_phys) :: rh_aft_PBL(ncol, pver) ! Relative humidity after PBL [%]
+ real(kind_phys) :: u_aft_PBL(ncol, pver) ! Zonal wind after PBL [m s-1]
+ real(kind_phys) :: v_aft_PBL(ncol, pver) ! Meridional wind after PBL [m s-1]
+ real(kind_phys) :: slflx(ncol, pverp) ! Liquid static energy flux at interfaces [W m-2]
+ real(kind_phys) :: qtflx(ncol, pverp) ! Total water flux at interfaces [kg m-2 s-1]
+ real(kind_phys) :: uflx(ncol, pverp) ! Zonal momentum flux at interfaces [N m-2]
+ real(kind_phys) :: vflx(ncol, pverp) ! Meridional momentum flux at interfaces [N m-2]
+ real(kind_phys) :: slflx_cg(ncol, pverp) ! Counter-gradient liquid static energy flux [W m-2]
+ real(kind_phys) :: qtflx_cg(ncol, pverp) ! Counter-gradient total water flux [kg m-2 s-1]
+ real(kind_phys) :: uflx_cg(ncol, pverp) ! Counter-gradient zonal momentum flux [N m-2]
+ real(kind_phys) :: vflx_cg(ncol, pverp) ! Counter-gradient meridional momentum flux [N m-2]
+ real(kind_phys) :: slten(ncol, pver) ! Liquid water static energy tendency [J kg-1 s-1]
+ real(kind_phys) :: qtten(ncol, pver) ! Total water mixing ratio tendency [kg kg-1 s-1]
+ real(kind_phys) :: tten(ncol, pver) ! Temperature tendency [K s-1]
+ real(kind_phys) :: rhten(ncol, pver) ! Relative humidity tendency [% s-1]
+ real(kind_phys) :: dtvke(ncol, pver) ! Normalized KE dissipation heating for history [K s-1]
+ real(kind_phys) :: dtv(ncol, pver) ! Normalized temperature tendency for history [K s-1]
+
+ ! Subset constituent upward ccpp constituent fluxes from coupler
+ real(kind_phys) :: cam_in_cflx_wv(ncol)
+ real(kind_phys) :: cam_in_cflx_cldliq(ncol)
+ real(kind_phys) :: cam_in_cflx_cldice(ncol)
+
+ ! Subset q0, q1, tendencies for constituents from coupler
+ real(kind_phys) :: q0_wv(ncol, pver)
+ real(kind_phys) :: q0_cldliq(ncol, pver)
+ real(kind_phys) :: q0_cldice(ncol, pver)
+ real(kind_phys) :: q1_cldliq(ncol, pver)
+ real(kind_phys) :: q1_cldice(ncol, pver)
+ real(kind_phys) :: tend_q_wv(ncol, pver)
+ real(kind_phys) :: tend_q_cldliq(ncol, pver)
+ real(kind_phys) :: tend_q_cldice(ncol, pver)
+
+ errmsg = ''
+ errflg = 0
+
+ ! Get constituent indices for wv, cldliq, cldice
+ call ccpp_const_get_idx(const_props, &
+ wv_stdname, &
+ const_wv_idx, errmsg, errflg)
+ if (errflg /= 0) return
+
+ call ccpp_const_get_idx(const_props, &
+ 'cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water', &
+ const_cldliq_idx, errmsg, errflg)
+ if (errflg /= 0) return
+
+ call ccpp_const_get_idx(const_props, &
+ 'cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water', &
+ const_cldice_idx, errmsg, errflg)
+ if (errflg /= 0) return
+
+ ! Extract constituent fluxes from coupler upward constituent fluxes
+ if(const_wv_idx > 0) then
+ cam_in_cflx_wv(:ncol) = cam_in_cflx(:ncol, const_wv_idx)
+ q0_wv(:ncol, :) = q0(:ncol, :, const_wv_idx)
+ tend_q_wv(:ncol, :) = tend_q(:ncol, :, const_wv_idx)
+ else
+ cam_in_cflx_wv(:ncol) = 0._kind_phys
+ q0_wv(:ncol, :) = 0._kind_phys
+ tend_q_wv(:ncol, :) = 0._kind_phys
+ endif
+
+ if(const_cldliq_idx > 0) then
+ cam_in_cflx_cldliq(:ncol) = cam_in_cflx(:ncol, const_cldliq_idx)
+ q0_cldliq(:ncol, :) = q0(:ncol, :, const_cldliq_idx)
+ q1_cldliq(:ncol, :) = q1(:ncol, :, const_cldliq_idx)
+ tend_q_cldliq(:ncol, :) = tend_q(:ncol, :, const_cldliq_idx)
+ else
+ cam_in_cflx_cldliq(:ncol) = 0._kind_phys
+ q0_cldliq(:ncol, :) = 0._kind_phys
+ q1_cldliq(:ncol, :) = 0._kind_phys
+ tend_q_cldliq(:ncol, :) = 0._kind_phys
+ endif
+
+ if(const_cldice_idx > 0) then
+ cam_in_cflx_cldice(:ncol) = cam_in_cflx(:ncol, const_cldice_idx)
+ q0_cldice(:ncol, :) = q0(:ncol, :, const_cldice_idx)
+ q1_cldice(:ncol, :) = q1(:ncol, :, const_cldice_idx)
+ tend_q_cldice(:ncol, :) = tend_q(:ncol, :, const_cldice_idx)
+ else
+ cam_in_cflx_cldice(:ncol) = 0._kind_phys
+ q0_cldice(:ncol, :) = 0._kind_phys
+ q1_cldice(:ncol, :) = 0._kind_phys
+ tend_q_cldice(:ncol, :) = 0._kind_phys
+ endif
+
+ ! Call the diagnostic profile calculation subroutine
+ call vertical_diffusion_diagnostic_profiles( &
+ ncol = ncol, &
+ pver = pver, &
+ pverp = pverp, &
+ dt = dt, &
+ latvap = latvap, &
+ latice = latice, &
+ zvir = zvir, &
+ cpair = cpair, &
+ gravit = gravit, &
+ rair = rair, &
+ pint = pint(:ncol, :pverp), &
+ pmid = pmid(:ncol, :pver), &
+ zi = zi(:ncol, :pverp), &
+ zm = zm(:ncol, :pver), &
+ kvh = kvh(:ncol, :pverp), &
+ kvm = kvm(:ncol, :pverp), &
+ cgh = cgh(:ncol, :pverp), &
+ cgs = cgs(:ncol, :pverp), &
+ cam_in_shf = cam_in_shf(:ncol), &
+ cam_in_cflx_wv = cam_in_cflx_wv(:ncol), &
+ cam_in_cflx_cldliq = cam_in_cflx_cldliq(:ncol), &
+ cam_in_cflx_cldice = cam_in_cflx_cldice(:ncol), &
+ tautotx = tautotx(:ncol), & ! total surface stresses to PBL coefficient scheme, not from provisional
+ tautoty = tautoty(:ncol), & ! total surface stresses to PBL coefficient scheme, not from provisional
+ t0 = t0(:ncol, :pver), &
+ q0_wv = q0_wv(:ncol, :pver), &
+ q0_cldliq = q0_cldliq(:ncol, :pver), &
+ q0_cldice = q0_cldice(:ncol, :pver), &
+ s0 = s0(:ncol, :pver), &
+ u0 = u0(:ncol, :pver), &
+ v0 = v0(:ncol, :pver), &
+ s1 = s1(:ncol, :pver), &
+ u1 = u1(:ncol, :pver), &
+ v1 = v1(:ncol, :pver), &
+ tend_s = tend_s(:ncol, :pver), &
+ tend_u = tend_u(:ncol, :pver), &
+ tend_v = tend_v(:ncol, :pver), &
+ tend_q_wv = tend_q_wv(:ncol, :pver), &
+ tend_q_cldliq = tend_q_cldliq(:ncol, :pver), &
+ tend_q_cldice = tend_q_cldice(:ncol, :pver), &
+ ! below output:
+ qt_pre_PBL = qt_pre_PBL(:ncol, :pver), &
+ sl_pre_PBL = sl_pre_PBL(:ncol, :pver), &
+ slv_pre_PBL = slv_pre_PBL(:ncol, :pver), &
+ rh_pre_PBL = rh_pre_PBL(:ncol, :pver), &
+ qt_aft_PBL = qt_aft_PBL(:ncol, :pver), &
+ sl_aft_PBL = sl_aft_PBL(:ncol, :pver), &
+ slv_aft_PBL = slv_aft_PBL(:ncol, :pver), &
+ qv_aft_PBL = qv_aft_PBL(:ncol, :pver), &
+ ql_aft_PBL = ql_aft_PBL(:ncol, :pver), &
+ qi_aft_PBL = qi_aft_PBL(:ncol, :pver), &
+ t_aft_PBL = t_aft_PBL(:ncol, :pver), &
+ rh_aft_PBL = rh_aft_PBL(:ncol, :pver), &
+ u_aft_PBL = u_aft_PBL(:ncol, :pver), &
+ v_aft_PBL = v_aft_PBL(:ncol, :pver), &
+ slflx = slflx(:ncol, :pverp), &
+ qtflx = qtflx(:ncol, :pverp), &
+ uflx = uflx(:ncol, :pverp), &
+ vflx = vflx(:ncol, :pverp), &
+ slflx_cg = slflx_cg(:ncol, :pverp), &
+ qtflx_cg = qtflx_cg(:ncol, :pverp), &
+ uflx_cg = uflx_cg(:ncol, :pverp), &
+ vflx_cg = vflx_cg(:ncol, :pverp), &
+ slten = slten(:ncol, :pver), &
+ qtten = qtten(:ncol, :pver), &
+ tten = tten(:ncol, :pver), &
+ rhten = rhten(:ncol, :pver))
+
+ ! Standard output variables
+ call history_out_field('QT', qt_aft_PBL)
+ call history_out_field('SL', sl_aft_PBL)
+ call history_out_field('SLV', slv_aft_PBL)
+ call history_out_field('SLFLX', slflx)
+ call history_out_field('QTFLX', qtflx)
+ call history_out_field('UFLX', uflx)
+ call history_out_field('VFLX', vflx)
+
+ ! Post-PBL profiles
+ call history_out_field('sl_aft_PBL', sl_aft_PBL)
+ call history_out_field('qt_aft_PBL', qt_aft_PBL)
+ call history_out_field('slv_aft_PBL', slv_aft_PBL)
+ call history_out_field('u_aft_PBL', u_aft_PBL)
+ call history_out_field('v_aft_PBL', v_aft_PBL)
+ call history_out_field('qv_aft_PBL', qv_aft_PBL)
+ call history_out_field('ql_aft_PBL', ql_aft_PBL)
+ call history_out_field('qi_aft_PBL', qi_aft_PBL)
+ call history_out_field('t_aft_PBL', t_aft_PBL)
+ call history_out_field('rh_aft_PBL', rh_aft_PBL)
+
+ ! PBL fluxes
+ call history_out_field('slflx_PBL', slflx)
+ call history_out_field('qtflx_PBL', qtflx)
+ call history_out_field('uflx_PBL', uflx)
+ call history_out_field('vflx_PBL', vflx)
+
+ ! Counter-gradient fluxes
+ call history_out_field('slflx_cg_PBL', slflx_cg)
+ call history_out_field('qtflx_cg_PBL', qtflx_cg)
+ call history_out_field('uflx_cg_PBL', uflx_cg)
+ call history_out_field('vflx_cg_PBL', vflx_cg)
+
+ ! PBL tendencies
+ call history_out_field('slten_PBL', slten)
+ call history_out_field('qtten_PBL', qtten)
+ call history_out_field('uten_PBL', tend_u)
+ call history_out_field('vten_PBL', tend_v)
+ call history_out_field('qvten_PBL', tend_q_wv)
+ call history_out_field('qlten_PBL', tend_q_cldliq)
+ call history_out_field('qiten_PBL', tend_q_cldice)
+ call history_out_field('tten_PBL', tten)
+ call history_out_field('rhten_PBL', rhten)
+
+ ! Pre-PBL profiles
+ call history_out_field('qt_pre_PBL', qt_pre_PBL)
+ call history_out_field('sl_pre_PBL', sl_pre_PBL)
+ call history_out_field('slv_pre_PBL', slv_pre_PBL)
+ call history_out_field('u_pre_PBL', u0)
+ call history_out_field('v_pre_PBL', v0)
+ call history_out_field('qv_pre_PBL', q0_wv)
+ call history_out_field('ql_pre_PBL', q0_cldliq)
+ call history_out_field('qi_pre_PBL', q0_cldice)
+ call history_out_field('t_pre_PBL', t0)
+ call history_out_field('rh_pre_PBL', rh_pre_PBL)
+
+ ! Vertical diffusion tendencies
+ call history_out_field('DUV', tend_u)
+ call history_out_field('DVV', tend_v)
+ call history_out_field('VERT_DIFF_Q', tend_q_wv)
+ call history_out_field('VERT_DIFF_CLDLIQ', tend_q_cldliq)
+ call history_out_field('VERT_DIFF_CLDICE', tend_q_cldice)
+
+ ! Normalize kinetic energy dissipation heating for history output
+ ! Convert from [J kg-1] to [K s-1]
+ dtvke(:ncol, :pver) = dtk(:ncol, :pver) / cpair / dt
+
+ ! Normalize dry static energy tendency for history output
+ ! Convert from [J kg-1 s-1] to [K s-1]
+ dtv(:ncol, :pver) = tend_s(:ncol, :pver) / cpair
+
+ ! Output diagnostics
+ call history_out_field('DTVKE', dtvke)
+ call history_out_field('DTV', dtv)
+
+ end subroutine vertical_diffusion_tendencies_diagnostics_run
+
+
+end module diffusion_solver_diagnostics
diff --git a/schemes/sima_diagnostics/diffusion_solver_diagnostics.meta b/schemes/sima_diagnostics/diffusion_solver_diagnostics.meta
new file mode 100644
index 00000000..fa52a731
--- /dev/null
+++ b/schemes/sima_diagnostics/diffusion_solver_diagnostics.meta
@@ -0,0 +1,268 @@
+[ccpp-table-properties]
+ name = vertical_diffusion_tendencies_diagnostics
+ type = scheme
+ dependencies = ../vertical_diffusion/vertical_diffusion_diagnostic_utils.F90
+
+[ccpp-arg-table]
+ name = vertical_diffusion_tendencies_diagnostics_init
+ type = scheme
+[ const_props ]
+ standard_name = ccpp_constituent_properties
+ units = none
+ type = ccpp_constituent_prop_ptr_t
+ dimensions = (number_of_ccpp_constituents)
+ intent = in
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-arg-table]
+ name = vertical_diffusion_tendencies_diagnostics_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pverp ]
+ standard_name = vertical_interface_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ dt ]
+ standard_name = timestep_for_physics
+ units = s
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ const_props ]
+ standard_name = ccpp_constituent_properties
+ units = none
+ type = ccpp_constituent_prop_ptr_t
+ dimensions = (number_of_ccpp_constituents)
+ intent = in
+[ latvap ]
+ standard_name = latent_heat_of_vaporization_of_water_at_0c
+ units = J kg-1
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ latice ]
+ standard_name = latent_heat_of_fusion_of_water_at_0c
+ units = J kg-1
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ zvir ]
+ standard_name = ratio_of_water_vapor_to_dry_air_gas_constants_minus_one
+ units = 1
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ cpair ]
+ standard_name = specific_heat_of_dry_air_at_constant_pressure
+ units = J kg-1 K-1
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ gravit ]
+ standard_name = standard_gravitational_acceleration
+ units = m s-2
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ rair ]
+ standard_name = gas_constant_of_dry_air
+ units = J kg-1 K-1
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ pmid ]
+ standard_name = air_pressure
+ units = Pa
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ pint ]
+ standard_name = air_pressure_at_interface
+ units = Pa
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ zi ]
+ standard_name = geopotential_height_wrt_surface_at_interface
+ units = m
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ zm ]
+ standard_name = geopotential_height_wrt_surface
+ units = m
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ kvh ]
+ standard_name = eddy_heat_diffusivity_at_interfaces
+ units = m2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ kvm ]
+ standard_name = eddy_momentum_diffusivity_at_interfaces
+ units = m2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ cgh ]
+ standard_name = counter_gradient_coefficient_for_vertical_diffusion_of_heat_at_interfaces
+ units = J kg-1 m-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ cgs ]
+ standard_name = counter_gradient_coefficient_for_vertical_diffusion_of_constituents_at_interfaces
+ units = s m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ cam_in_shf ]
+ standard_name = surface_upward_sensible_heat_flux_from_coupler
+ units = W m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ cam_in_cflx ]
+ standard_name = surface_upward_ccpp_constituent_fluxes_from_coupler
+ units = kg m-2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, number_of_ccpp_constituents)
+ intent = in
+[ tautotx ]
+ # these are the taux input to the underlying boundary layer scheme providing the inputs
+ # to the diffusion solver.
+ standard_name = eastward_total_stress_at_surface_for_holtslag_boville_boundary_layer_scheme
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ tautoty ]
+ # these are the tauy input to the underlying boundary layer scheme providing the inputs
+ # to the diffusion solver.
+ standard_name = northward_total_stress_at_surface_for_holtslag_boville_boundary_layer_scheme
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ t0 ]
+ standard_name = air_temperature
+ units = K
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ q0 ]
+ standard_name = ccpp_constituents
+ units = none
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccpp_constituents)
+ intent = in
+[ s0 ]
+ standard_name = dry_static_energy
+ units = J kg-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ u0 ]
+ standard_name = eastward_wind
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ v0 ]
+ standard_name = northward_wind
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ q1 ]
+ standard_name = ccpp_constituents_after_vertical_diffusion
+ units = none
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccpp_constituents)
+ intent = in
+[ s1 ]
+ standard_name = dry_static_energy_after_vertical_diffusion
+ units = J kg-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ u1 ]
+ standard_name = eastward_wind_after_vertical_diffusion
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ v1 ]
+ standard_name = northward_wind_after_vertical_diffusion
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ dtk ]
+ standard_name = kinetic_energy_dissipation_in_atmosphere_boundary_layer
+ units = J kg-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ tend_s ]
+ standard_name = tendency_of_dry_air_enthalpy_at_constant_pressure
+ units = J kg-1 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ tend_u ]
+ standard_name = tendency_of_eastward_wind
+ units = m s-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ tend_v ]
+ standard_name = tendency_of_northward_wind
+ units = m s-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ tend_q ]
+ standard_name = ccpp_constituent_tendencies
+ units = none
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccpp_constituents)
+ intent = in
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
diff --git a/schemes/sima_diagnostics/holtslag_boville_diff_diagnostics.F90 b/schemes/sima_diagnostics/holtslag_boville_diff_diagnostics.F90
new file mode 100644
index 00000000..5343c1ce
--- /dev/null
+++ b/schemes/sima_diagnostics/holtslag_boville_diff_diagnostics.F90
@@ -0,0 +1,93 @@
+! Diagnostics for cloud fraction
+module holtslag_boville_diff_diagnostics
+ use ccpp_kinds, only: kind_phys
+
+ implicit none
+ private
+ save
+
+ public :: holtslag_boville_diff_diagnostics_init
+ public :: holtslag_boville_diff_diagnostics_run
+
+contains
+
+ !> \section arg_table_holtslag_boville_diff_diagnostics_init Argument Table
+ !! \htmlinclude holtslag_boville_diff_diagnostics_init.html
+ subroutine holtslag_boville_diff_diagnostics_init(errmsg, errflg)
+ use cam_history, only: history_add_field
+ use cam_history_support, only: horiz_only
+
+ character(len=512), intent(out) :: errmsg
+ integer, intent(out) :: errflg
+
+ ! Local variables:
+
+ errmsg = ''
+ errflg = 0
+
+ ! History add field calls
+ call history_add_field('HB_ri', 'Richardson Number (HB scheme)', 'lev', 'inst', '1')
+
+ call history_add_field('USTAR', 'Surface friction velocity', horiz_only, 'inst', 'm s-1')
+ call history_add_field('obklen', 'Obukhov length', horiz_only, 'avg', 'm')
+ call history_add_field('PBLH', 'Planetary boundary layer height', horiz_only, 'avg', 'm')
+
+ call history_add_field('KVH', 'Vertical diffusion diffusivity coefficient for heat/moisture', 'ilev', 'avg', 'm2 s-1')
+ call history_add_field('KVM', 'Vertical diffusion diffusivity coefficient for momentum', 'ilev', 'avg', 'm2 s-1')
+ call history_add_field('CGS', 'Counter-gradient coefficient on surface kinematic fluxes', 'ilev', 'avg', 's m-2')
+
+ call history_add_field('TKE', 'Specific turbulent kinetic energy at interfaces', 'ilev', 'avg', 'm2 s-2')
+ call history_add_field('TPERT', 'Perturbation temperature (eddies in PBL)', horiz_only, 'avg', 'K')
+ call history_add_field('QPERT', 'Perturbation specific humidity (eddies in PBL)', horiz_only, 'avg', 'kg kg-1')
+
+ end subroutine holtslag_boville_diff_diagnostics_init
+
+ !> \section arg_table_holtslag_boville_diff_diagnostics_run Argument Table
+ !! \htmlinclude holtslag_boville_diff_diagnostics_run.html
+ subroutine holtslag_boville_diff_diagnostics_run( &
+ ri, &
+ ustar, obklen, &
+ pblh, &
+ kvh, kvm, cgs, &
+ tke, &
+ tpert, qpert, &
+ errmsg, errflg)
+
+ use cam_history, only: history_out_field
+
+ ! Input parameters
+ real(kind_phys), intent(in) :: ri(:,:)
+ real(kind_phys), intent(in) :: ustar(:)
+ real(kind_phys), intent(in) :: obklen(:)
+ real(kind_phys), intent(in) :: pblh(:)
+
+ real(kind_phys), intent(in) :: kvh(:,:)
+ real(kind_phys), intent(in) :: kvm(:,:)
+ real(kind_phys), intent(in) :: cgs(:,:)
+ real(kind_phys), intent(in) :: tke(:,:)
+
+ real(kind_phys), intent(in) :: tpert(:)
+ real(kind_phys), intent(in) :: qpert(:)
+
+ character(len=512), intent(out) :: errmsg
+ integer, intent(out) :: errflg
+
+ errmsg = ''
+ errflg = 0
+
+ call history_out_field('HB_ri', ri)
+ call history_out_field('USTAR', ustar)
+ call history_out_field('obklen', obklen)
+ call history_out_field('PBLH', pblh)
+
+ call history_out_field('KVH', kvh)
+ call history_out_field('KVM', kvm)
+ call history_out_field('CGS', cgs)
+ call history_out_field('TKE', tke)
+
+ call history_out_field('TPERT', tpert)
+ call history_out_field('QPERT', qpert)
+
+ end subroutine holtslag_boville_diff_diagnostics_run
+
+end module holtslag_boville_diff_diagnostics
diff --git a/schemes/sima_diagnostics/holtslag_boville_diff_diagnostics.meta b/schemes/sima_diagnostics/holtslag_boville_diff_diagnostics.meta
new file mode 100644
index 00000000..e1fb9a7c
--- /dev/null
+++ b/schemes/sima_diagnostics/holtslag_boville_diff_diagnostics.meta
@@ -0,0 +1,95 @@
+[ccpp-table-properties]
+ name = holtslag_boville_diff_diagnostics
+ type = scheme
+
+[ccpp-arg-table]
+ name = holtslag_boville_diff_diagnostics_init
+ type = scheme
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-arg-table]
+ name = holtslag_boville_diff_diagnostics_run
+ type = scheme
+[ ri ]
+ standard_name = gradient_richardson_number
+ units = 1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ ustar ]
+ standard_name = surface_friction_velocity
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ obklen ]
+ standard_name = obukhov_length
+ units = m
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ pblh ]
+ standard_name = atmosphere_boundary_layer_thickness
+ units = m
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ kvh ]
+ standard_name = eddy_heat_diffusivity_at_interfaces
+ units = m2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ kvm ]
+ standard_name = eddy_momentum_diffusivity_at_interfaces
+ units = m2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ cgs ]
+ standard_name = counter_gradient_coefficient_for_vertical_diffusion_of_constituents_at_interfaces
+ units = s m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ tke ]
+ standard_name = specific_turbulent_kinetic_energy_at_interfaces
+ units = m2 s-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ tpert ]
+ standard_name = convective_temperature_perturbation_due_to_pbl_eddies
+ units = K
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ qpert ]
+ standard_name = convective_water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water_perturbation_due_to_pbl_eddies
+ units = kg kg-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
diff --git a/schemes/vertical_diffusion/diffusion_solver.F90 b/schemes/vertical_diffusion/diffusion_solver.F90
new file mode 100644
index 00000000..5bc6cadf
--- /dev/null
+++ b/schemes/vertical_diffusion/diffusion_solver.F90
@@ -0,0 +1,1028 @@
+! Vertical diffusion module.
+! Solves vertical diffusion equations using a tri-diagonal solver
+! The module will also apply countergradient fluxes
+! Includes interstitials for applying vertical diffusion tendencies and
+! other intermediate calculations formerly in vertical_diffusion_tend.
+!
+! This is a reduced functionality version;
+! molecular diffusion is unsupported until WACCM-X is ported.
+!
+! Original authors: B. Boville and others, 1991-2004
+! Modularized: J. McCaa, September 2004.
+! Updates: Sungsu Park, August 2006 - January 2010
+! Reduced/CCPPized: Haipeng Lin, May 2025
+module diffusion_solver
+ use ccpp_kinds, only: kind_phys
+
+ implicit none
+ private
+ save
+
+ ! CCPP-compliant public interfaces
+
+ public :: vertical_diffusion_set_temperature_at_toa_default_run
+ public :: vertical_diffusion_interpolate_to_interfaces_run ! Interpolate t, rho, rair to interfaces
+ public :: implicit_surface_stress_add_drag_coefficient_run ! Add implicit turbulent surface stress (if do_iss)
+ public :: vertical_diffusion_wind_damping_rate_run ! Compute wind damping rate from drag coeff
+ ! Beljaars will add wind damping rates at this point
+
+ public :: vertical_diffusion_diffuse_horizontal_momentum_timestep_init
+ public :: vertical_diffusion_diffuse_horizontal_momentum_run ! Vertical diffusion solver (horizontal momentum)
+
+ public :: vertical_diffusion_set_dry_static_energy_at_toa_zero_run
+ public :: vertical_diffusion_set_dry_static_energy_at_toa_molecdiff_run
+ public :: vertical_diffusion_diffuse_dry_static_energy_run ! Vertical diffusion solver (dry static energy)
+
+ public :: vertical_diffusion_diffuse_tracers_init
+ public :: vertical_diffusion_diffuse_tracers_run ! Vertical diffusion solver (tracers)
+
+ public :: vertical_diffusion_tendencies_run ! Using output from solver compute phys. tend
+
+ ! Parameters for implicit surface stress treatment
+ real(kind_phys), parameter :: wsmin = 1._kind_phys ! Minimum sfc wind speed for estimating frictional
+ ! transfer velocity ksrf. [m s-1]
+ real(kind_phys), parameter :: ksrfmin = 1.e-4_kind_phys ! Minimum surface drag coefficient [kg s-1 m-2]
+ real(kind_phys), parameter :: timeres = 7200._kind_phys ! Relaxation time scale of residual stress (>= dt) [s]
+
+contains
+
+ ! Set (default) temperature used at top of atmosphere for interpolation.
+ ! In the future, this can use the upper boundary condition provided for temperature,
+ ! or a different extrapolation in WACCM-X mode.
+!> \section arg_table_vertical_diffusion_set_temperature_at_toa_default_run Argument Table
+!! \htmlinclude arg_table_vertical_diffusion_set_temperature_at_toa_default_run.html
+ pure subroutine vertical_diffusion_set_temperature_at_toa_default_run( &
+ ncol, pver, &
+ t, &
+ t_toai, &
+ errmsg, errflg)
+
+ ! Input arguments
+ integer, intent(in) :: ncol
+ integer, intent(in) :: pver
+ real(kind_phys), intent(in) :: t(:, :) ! Temperature [K]
+
+ ! Output arguments
+ real(kind_phys), intent(out) :: t_toai(:) ! Temperature at interface above TOA for vertical diffusion [K]
+ character(len=512), intent(out) :: errmsg
+ integer, intent(out) :: errflg
+
+ errmsg = ''
+ errflg = 0
+
+ ! Default
+ t_toai(:ncol) = t(:ncol, 1) ! Use temperature at top model level [K]
+
+ ! Eventually, if ubc_fixed_temp is available, use ubc_t.
+
+ ! Eventually, in WACCM-X, use extrapolated version:
+ ! if (waccmx_mode) then
+ ! ! For WACCM-X, set ubc temperature to extrapolate from next two lower interface level temperatures
+ ! ! the original formulation is:
+ ! ! t_toai(:ncol) = 1.5_r8*tint(:ncol,2)-.5_r8*tint(:ncol,3)
+ ! ! this appears to be:
+ ! ! = tint(:ncol,2) + 0.5_r8*(tint(:ncol,2) - tint(:ncol,3))
+ ! ! assuming that the extrapolated gradient is 1/2 of the temperature gradient between the lower interfaces
+ ! ! because the interpolation will be done later in the CCPP-ized scheme, formulate this in terms of
+ ! ! the temperature (at midpoints):
+ ! t_toai(:ncol) = 1.5_r8*(state%t(:ncol,2)+state%t(:ncol,1))/2._r8-.5_r8*(state%t(:ncol,3)+state%t(:ncol,2))/2._r8
+
+ end subroutine vertical_diffusion_set_temperature_at_toa_default_run
+
+ ! Interpolates temperature, air density (moist and dry),
+ ! and sets gas constant (not constituent dependent in non-WACCM-X mode)
+ ! at interfaces for vertical diffusion
+!> \section arg_table_vertical_diffusion_interpolate_to_interfaces_run Argument Table
+!! \htmlinclude arg_table_vertical_diffusion_interpolate_to_interfaces_run.html
+ pure subroutine vertical_diffusion_interpolate_to_interfaces_run( &
+ ncol, pver, pverp, &
+ gravit, &
+ rair, rairv, &
+ flag_for_constituent_dependent_gas_constant, &
+ t, &
+ t_toai, &
+ pint, pintdry, &
+ ! below output
+ ti, rairi, rhoi, rhoi_dry, &
+ dpidz_sq, &
+ errmsg, errflg)
+
+ ! Input arguments
+ integer, intent(in) :: ncol ! Number of columns
+ integer, intent(in) :: pver ! Number of vertical layers
+ integer, intent(in) :: pverp ! Number of vertical interfaces (pver+1)
+ real(kind_phys), intent(in) :: gravit
+ real(kind_phys), intent(in) :: rair ! Gas constant for dry air [J kg-1 K-1]
+ real(kind_phys), intent(in) :: rairv(:,:) ! composition_dependent_gas_constant_of_dry_air [J kg-1 K-1]
+ logical, intent(in) :: flag_for_constituent_dependent_gas_constant
+ real(kind_phys), intent(in) :: t(:,:) ! Air temperature at midpoints [K]
+ real(kind_phys), intent(in) :: t_toai(:) ! Air temperature to use at interface above TOA [K]
+ real(kind_phys), intent(in) :: pint(:,:) ! Air pressure at interfaces [Pa]
+ real(kind_phys), intent(in) :: pintdry(:,:) ! Dry air pressure at interfaces [Pa]
+
+ ! Output variables
+ real(kind_phys), intent(out) :: ti(:,:) ! Air temperature at interfaces [K]
+ real(kind_phys), intent(out) :: rairi(:,:) ! Gas constant for dry air at interfaces [J kg-1 K-1]
+ real(kind_phys), intent(out) :: rhoi(:,:) ! Air density at interfaces [kg m-3]
+ real(kind_phys), intent(out) :: rhoi_dry(:,:) ! Dry air density at interfaces [kg m-3]
+ real(kind_phys), intent(out) :: dpidz_sq(:,:) ! Square of derivative of pressure with height (moist) [kg2 m-4 s-4]
+ character(len=512), intent(out) :: errmsg ! Error message
+ integer, intent(out) :: errflg ! Error flag
+
+ ! Local variables
+ integer :: i, k
+
+ errmsg = ''
+ errflg = 0
+
+ ! Interpolate temperature to interfaces
+ ti(:ncol,1) = t_toai(:ncol)
+ do k = 2, pver
+ do i = 1, ncol
+ ti(i,k) = 0.5_kind_phys * (t(i,k) + t(i,k-1))
+ end do
+ end do
+ ! Assume temperature at surface interface is the same as the midpoint value in the lowest layer:
+ ti(:ncol,pverp) = t(:ncol,pver)
+
+ ! Set gas constant at interfaces
+ if(flag_for_constituent_dependent_gas_constant) then
+ rairi(:ncol,1) = rairv(:ncol,1)
+ do k = 2, pver
+ do i = 1, ncol
+ rairi(i,k) = 0.5_kind_phys * (rairv(i,k)+rairv(i,k-1))
+ end do
+ end do
+
+ ! Assume gas constant at surface interface is the same as the midpoint value in the lowest layer:
+ rairi(:ncol,pverp) = rairv(:ncol,pver)
+ else
+ rairi(:ncol,:pverp) = rair
+ endif
+
+ ! Compute rho at interfaces.
+ do k = 1, pverp
+ do i = 1, ncol
+ rhoi(i,k) = pint(i,k) / (rairi(i,k)*ti(i,k))
+ end do
+ end do
+
+ ! Compute rho_dry at interfaces.
+ do k = 1, pverp
+ do i = 1, ncol
+ rhoi_dry(i,k) = pintdry(i,k) / (rairi(i,k)*ti(i,k))
+ end do
+ end do
+
+ ! Compute square of derivative of pressure with height at interfaces.
+ ! Note that this equation assumes hydrostatic balance, i.e. the *derivative* dp/dz is g*rho
+ dpidz_sq = gravit*rhoi(:ncol, :)
+ dpidz_sq = dpidz_sq*dpidz_sq
+
+ ! FIXME: The following two lines are kept in only to preserve answers;
+ ! they really should be taken out completely.
+ dpidz_sq(:, 1) = gravit*(pint(:, 1)/(rairv(:ncol, 1)*t(:ncol, 1)))
+ dpidz_sq(:, 1) = dpidz_sq(:, 1)*dpidz_sq(:, 1)
+
+ end subroutine vertical_diffusion_interpolate_to_interfaces_run
+
+ ! Add turbulent surface stress to total surface drag coefficient, if implicit
+!> \section arg_table_implicit_surface_stress_add_drag_coefficient_run Argument Table
+!! \htmlinclude arg_table_implicit_surface_stress_add_drag_coefficient_run.html
+ pure subroutine implicit_surface_stress_add_drag_coefficient_run( &
+ ncol, pver, &
+ do_iss, &
+ taux, tauy, &
+ u0, v0, &
+ ! below output
+ ksrf, &
+ errmsg, errflg)
+
+ ! Input arguments
+ integer, intent(in) :: ncol
+ integer, intent(in) :: pver
+ logical, intent(in) :: do_iss ! Use implicit turbulent surface stress computation
+ real(kind_phys), intent(in) :: taux(:) ! Surface zonal stress [N m-2]
+ ! Input u-momentum per unit time per unit area into the atmosphere.
+ real(kind_phys), intent(in) :: tauy(:) ! Surface meridional stress [N m-2]
+ ! Input v-momentum per unit time per unit area into the atmosphere.
+
+ real(kind_phys), intent(in) :: u0(:,:) ! Input u-wind [m s-1]
+ real(kind_phys), intent(in) :: v0(:,:) ! Input v-wind [m s-1]
+
+ ! Input/output arguments
+ real(kind_phys), intent(inout) :: ksrf(:) ! total surface drag coefficient [kg m-2 s-1]
+
+ ! Output arguments
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ ! Local variables
+ integer :: i
+ real(kind_phys) :: ws(ncol) ! Lowest-level wind speed [m s-1]
+ real(kind_phys) :: ksrfturb(ncol) ! Surface drag coefficient of 'normal' stress. > 0.
+ ! Virtual mass input per unit time per unit area [kg m-2 s-1]
+ real(kind_phys) :: tau(ncol) ! Turbulent surface stress (not including mountain stress)
+
+ errmsg = ''
+ errflg = 0
+
+ ! Do 'Implicit Surface Stress' treatment for numerical stability
+ ! in the lowest model layer.
+ if (do_iss) then
+ ! Add surface drag coefficient for implicit diffusion
+ do i = 1, ncol
+ ws(i) = max(sqrt(u0(i, pver)**2 + v0(i, pver)**2), wsmin)
+ tau(i) = sqrt(taux(i)**2 + tauy(i)**2)
+ ksrfturb(i) = max(tau(i)/ws(i), ksrfmin)
+ end do
+ ksrf(:ncol) = ksrf(:ncol) + ksrfturb(:ncol) ! Do all surface stress implicitly
+ endif
+
+ end subroutine implicit_surface_stress_add_drag_coefficient_run
+
+ ! Computes rate at which wind is exponentially damped by surface stress
+ ! based on total surface stress (after all surface stress is added)
+!> \section arg_table_vertical_diffusion_wind_damping_rate_run Argument Table
+!! \htmlinclude arg_table_vertical_diffusion_wind_damping_rate_run.html
+ pure subroutine vertical_diffusion_wind_damping_rate_run( &
+ ncol, pver, &
+ gravit, &
+ p, &
+ ksrf, &
+ tau_damp_rate, &
+ errmsg, errflg)
+
+ use coords_1d, only: Coords1D
+
+ ! Input arguments
+ integer, intent(in) :: ncol
+ integer, intent(in) :: pver
+ real(kind_phys), intent(in) :: gravit
+ type(coords1d), intent(in) :: p ! Pressure coordinates [Pa]
+ real(kind_phys), intent(in) :: ksrf(:) ! total surface drag coefficient [kg m-2 s-1]
+
+ ! Output arguments
+ real(kind_phys), intent(out) :: tau_damp_rate(:,:) ! Rate at which external (surface) stress damps wind speeds [s-1]
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ errmsg = ''
+ errflg = 0
+
+ ! In most layers, no damping at all.
+ tau_damp_rate(:ncol,:pver) = 0._kind_phys
+
+ ! Compute damp rate based on total surface drag coeff
+ !
+ ! Physical interpretation:
+ ! ksrf is stress per unit wind speed.
+ ! p%del / gravit is approximately the mass in the layer per unit of
+ ! surface area.
+ ! Therefore, gravit*ksrf/p%del is the acceleration of wind per unit
+ ! wind speed, i.e. the rate at which wind is exponentially damped by
+ ! surface stress.
+ tau_damp_rate(:ncol, pver) = -gravit*ksrf(:ncol)*p%rdel(:ncol, pver)
+
+ end subroutine vertical_diffusion_wind_damping_rate_run
+
+ ! Zero surface stresses at timestep initialization
+ ! Any surface stresses, e.g., Turbulent Mountain Stress
+ ! will be added to this later as they are computed.
+!> \section arg_table_vertical_diffusion_diffuse_horizontal_momentum_timestep_init Argument Table
+!! \htmlinclude arg_table_vertical_diffusion_diffuse_horizontal_momentum_timestep_init.html
+ pure subroutine vertical_diffusion_diffuse_horizontal_momentum_timestep_init( &
+ ncol, pver, &
+ ksrf, &
+ tau_damp_rate, &
+ errmsg, errflg)
+
+ integer, intent(in) :: ncol
+ integer, intent(in) :: pver
+
+ ! Output (initial zeroed) arguments
+ real(kind_phys), intent(out) :: ksrf(:) ! total surface drag coefficient [kg m-2 s-1]
+ real(kind_phys), intent(out) :: tau_damp_rate(:,:) ! Rate at which external (surface) stress damps wind speeds [s-1]
+
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ errmsg = ''
+ errflg = 0
+
+ ksrf(:ncol) = 0._kind_phys
+
+ ! In most layers, no damping at all.
+ tau_damp_rate(:ncol,:pver) = 0._kind_phys
+
+ end subroutine vertical_diffusion_diffuse_horizontal_momentum_timestep_init
+
+ ! Driver routine to compute vertical diffusion of horizontal momentum,
+ ! including the loss of momentum due to frictional (KE to heat) dissipation.
+ ! Turbulent diffusivities and boundary layer nonlocal transport terms are
+ ! obtained from the turbulence module.
+!> \section arg_table_vertical_diffusion_diffuse_horizontal_momentum_run Argument Table
+!! \htmlinclude arg_table_vertical_diffusion_diffuse_horizontal_momentum_run.html
+ subroutine vertical_diffusion_diffuse_horizontal_momentum_run( &
+ ncol, pver, pverp, &
+ dt, &
+ rair, gravit, &
+ do_iss, &
+ am_correction, &
+ itaures, &
+ do_beljaars, &
+ p, t, rhoi, &
+ taux, tauy, &
+ tau_damp_rate, &
+ ksrftms, &
+ dragblj, &
+ kvm, &
+ dpidz_sq, &
+ u0, v0, dse0, &
+ tauresx, tauresy, &
+ dtk, &
+ tautmsx, tautmsy, &
+ u1, v1, dse1, &
+ errmsg, errflg)
+
+ use coords_1d, only: Coords1D
+ use vertical_diffusion_solver, only: fin_vol_solve
+
+ ! Input Arguments
+ integer, intent(in) :: ncol
+ integer, intent(in) :: pver
+ integer, intent(in) :: pverp
+ real(kind_phys), intent(in) :: dt ! physics timestep [s]
+ real(kind_phys), intent(in) :: rair
+ real(kind_phys), intent(in) :: gravit
+ logical, intent(in) :: do_iss ! Use implicit turbulent surface stress computation
+ logical, intent(in) :: am_correction ! Do angular momentum conservation correction
+ logical, intent(in) :: itaures ! Flag for updating tauresx tauresy in this subroutine.
+ logical, intent(in) :: do_beljaars ! Flag indicating Beljaars drag
+ type(coords1d), intent(in) :: p ! Pressure coordinates [Pa]
+ real(kind_phys), intent(in) :: t(:, :) ! Temperature [K]
+ real(kind_phys), intent(in) :: rhoi(:, :) ! Density of air at interfaces [kg m-3]
+ real(kind_phys), intent(in) :: taux(:) ! Surface zonal stress [N m-2]
+ ! Input u-momentum per unit time per unit area into the atmosphere.
+ real(kind_phys), intent(in) :: tauy(:) ! Surface meridional stress [N m-2]
+ ! Input v-momentum per unit time per unit area into the atmosphere.
+ real(kind_phys), intent(in) :: tau_damp_rate(:,:) ! Rate at which external (surface) stress damps wind speeds [s-1]
+ real(kind_phys), intent(in) :: ksrftms(:) ! Surface drag coefficient for turbulent mountain stress. > 0. [kg m-2 s-1]
+ real(kind_phys), intent(in) :: dragblj(:, :) ! Drag profile from Beljaars SGO form drag > 0. [s-1]
+ real(kind_phys), intent(in) :: kvm(:, :) ! Eddy viscosity (Eddy diffusivity for momentum) [m^2 s-1], interfaces
+ real(kind_phys), intent(in) :: dpidz_sq(:,:) ! Square of derivative of pressure with height (moist) [kg2 m-4 s-4], interfaces
+
+ real(kind_phys), intent(in) :: u0(:,:) ! Input u-wind [m s-1]
+ real(kind_phys), intent(in) :: v0(:,:) ! Input v-wind [m s-1]
+ real(kind_phys), intent(in) :: dse0(:,:) ! Input Dry static energy [J kg-1]
+
+ ! Input: Reserved surface stress at previous time step
+ ! Output: Reserved surface stress at current time step
+ real(kind_phys), intent(inout) :: tauresx(:)
+ real(kind_phys), intent(inout) :: tauresy(:)
+
+ ! Output Arguments
+ real(kind_phys), intent(out) :: dtk(:, :) ! T tendency from KE dissipation [J kg-1]
+
+ ! Recomputed TMS stresses using provisionally updated winds:
+ real(kind_phys), intent(out) :: tautmsx(:) ! Implicit zonal turbulent mountain surface stress [N m-2]
+ real(kind_phys), intent(out) :: tautmsy(:) ! Implicit meridional turbulent mountain surface stress [N m-2]
+
+ ! Outputs from vertical diffusion will be converted to physics tendencies in separate scheme
+ ! These provisional outputs are used to compute diagnostics, and thus need to be kept in non-tend form
+ real(kind_phys), intent(out) :: u1(:,:) ! After vertical diffusion u-wind [m s-1]
+ real(kind_phys), intent(out) :: v1(:,:) ! After vertical diffusion v-wind [m s-1]
+ real(kind_phys), intent(out) :: dse1(:,:) ! After vertical diffusion Dry static energy [J kg-1]
+
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ ! Local variables
+ integer :: i, k
+
+ real(kind_phys) :: flux_to_state_conversion_factor(ncol)
+ real(kind_phys) :: tmpi1(ncol, pverp) ! Interface KE dissipation
+ real(kind_phys) :: tmpi2(ncol, pverp) ! dt*(g*rho)**2/dp at interfaces
+ real(kind_phys) :: ramda ! dt/timeres [unitless]
+
+ real(kind_phys) :: keg_in(ncol, pver) ! KE on entry to subroutine
+ real(kind_phys) :: keg_out(ncol, pver) ! KE after U and V dissipation/diffusion
+
+ real(kind_phys) :: tautotx(ncol) ! Total surface stress (zonal)
+ real(kind_phys) :: tautoty(ncol) ! Total surface stress (meridional)
+
+ real(kind_phys) :: dinp_u(ncol, pverp) ! Vertical difference at interfaces, input u
+ real(kind_phys) :: dinp_v(ncol, pverp) ! Vertical difference at interfaces, input v
+ real(kind_phys) :: dout_u ! Vertical difference at interfaces, output u
+ real(kind_phys) :: dout_v ! Vertical difference at interfaces, output v
+
+ real(kind_phys) :: usum_in(ncol) ! Vertical integral of input u-momentum. Total zonal
+ ! momentum per unit area in column (sum of u*dp/g) [kg m s-1 m-2]
+ real(kind_phys) :: vsum_in(ncol) ! Vertical integral of input v-momentum. Total meridional
+ ! momentum per unit area in column (sum of v*dp/g) [kg m s-1 m-2]
+ real(kind_phys) :: usum_out(ncol) ! Vertical integral of u-momentum after doing implicit diffusion
+ real(kind_phys) :: vsum_out(ncol) ! Vertical integral of v-momentum after doing implicit diffusion
+
+ real(kind_phys) :: tauimpx(ncol) ! Actual net stress added at the current step other than mountain stress
+ real(kind_phys) :: tauimpy(ncol) ! Actual net stress added at the current step other than mountain stress
+ real(kind_phys) :: taubljx(ncol) ! recomputed explicit/residual beljaars stress
+ real(kind_phys) :: taubljy(ncol) ! recomputed explicit/residual beljaars stress
+
+ errmsg = ''
+ errflg = 0
+
+ ! Set initial iterative outputs to input values
+ u1(:ncol,:pver) = u0(:ncol,:pver)
+ v1(:ncol,:pver) = v0(:ncol,:pver)
+ dse1(:ncol,:pver) = dse0(:ncol,:pver)
+
+ ! necessary temporaries used in computation
+ flux_to_state_conversion_factor(:ncol) = dt*gravit*p%rdel(:, pver)
+
+ ! second term here uses dpidz_sq(:,1) as defined because
+ ! the actual dpidz_sq term is "kept in only to preserve answers"
+ ! but it is actually not used here ...
+ tmpi2(:ncol, 1) = dt*(gravit*rhoi(:ncol,1)*gravit*rhoi(:ncol,1))/(p%mid(:, 1) - p%ifc(:, 1))
+ tmpi2(:ncol, 2:pver) = dt*dpidz_sq(:, 2:pver)*p%rdst
+
+ !---------------------------- !
+ ! Diffuse Horizontal Momentum !
+ !---------------------------- !
+ ! Compute the vertical upward differences of the input u,v for KE dissipation
+ ! at each interface.
+ ! Velocity = 0 at surface, so difference at the bottom interface is -u,v(pver)
+ ! These 'dinp_u, dinp_v' are computed using the non-diffused input wind.
+
+ do i = 1, ncol
+ dinp_u(i, 1) = 0._kind_phys
+ dinp_v(i, 1) = 0._kind_phys
+ dinp_u(i, pver + 1) = -u0(i, pver)
+ dinp_v(i, pver + 1) = -v0(i, pver)
+ end do
+ do k = 2, pver
+ do i = 1, ncol
+ dinp_u(i, k) = u0(i, k) - u0(i, k - 1)
+ dinp_v(i, k) = v0(i, k) - v0(i, k - 1)
+ end do
+ end do
+
+ ! -------------------------------------------------------------- !
+ ! Do 'Implicit Surface Stress' treatment for numerical stability !
+ ! in the lowest model layer. !
+ ! -------------------------------------------------------------- !
+ if (do_iss) then
+ ! Vertical integration of input momentum.
+ ! This is total horizontal momentum per unit area [ kg*m/s/m2 ] in each column.
+ ! Note (u,v) are the raw input to the PBL scheme, not the
+ ! provisionally-marched ones within the iteration loop of the PBL scheme.
+
+ do i = 1, ncol
+ usum_in(i) = 0._kind_phys
+ vsum_in(i) = 0._kind_phys
+ do k = 1, pver
+ usum_in(i) = usum_in(i) + (1._kind_phys/gravit)*u0(i, k)*p%del(i, k)
+ vsum_in(i) = vsum_in(i) + (1._kind_phys/gravit)*v0(i, k)*p%del(i, k)
+ end do
+ end do
+
+ ! Add residual stress of previous time step explicitly into the lowest
+ ! model layer with a relaxation time scale of 'timeres'.
+
+ if (am_correction) then
+ ! preserve time-mean torque
+ ramda = 1._kind_phys
+ else
+ ramda = dt/timeres
+ end if
+
+ u1(:ncol, pver) = u1(:ncol, pver) + flux_to_state_conversion_factor(:ncol)*tauresx(:ncol)*ramda
+ v1(:ncol, pver) = v1(:ncol, pver) + flux_to_state_conversion_factor(:ncol)*tauresy(:ncol)*ramda
+
+ else ! .not. do_iss
+ ! In this case, do 'turbulent mountain stress' implicitly,
+ ! but do 'normal turbulent stress' explicitly.
+ ! In this case, there is no 'residual stress' as long as 'tms' is
+ ! treated in a fully implicit way, which is true.
+
+ ! Do 'normal stress' explicitly
+ u1(:ncol, pver) = u1(:ncol, pver) + flux_to_state_conversion_factor(:ncol)*taux(:ncol)
+ v1(:ncol, pver) = v1(:ncol, pver) + flux_to_state_conversion_factor(:ncol)*tauy(:ncol)
+ end if ! End of 'do iss' (implicit surface stress)
+
+ ! --------------------------------------------------------------------------------------- !
+ ! Diffuse horizontal momentum implicitly using tri-diagnonal matrix. !
+ ! The 'u,v' are input-output: the output 'u,v' are implicitly diffused winds. !
+ ! For implicit 'normal' stress : ksrf = ksrftms + ksrfturb, !
+ ! u(pver) : explicitly include 'residual normal' stress !
+ ! For explicit 'normal' stress : ksrf = ksrftms !
+ ! u(pver) : explicitly include 'normal' stress !
+ ! Note that in all the two cases above, 'tms' is fully implicitly treated. !
+ ! --------------------------------------------------------------------------------------- !
+
+ v1(:ncol, :) = fin_vol_solve(dt, p, v1(:ncol, :), ncol, pver, &
+ coef_q=tau_damp_rate(:ncol,:pver), &
+ coef_q_diff=kvm(:ncol, :)*dpidz_sq(:ncol, :))
+
+ u1(:ncol, :) = fin_vol_solve(dt, p, u1(:ncol, :), ncol, pver, &
+ coef_q=tau_damp_rate(:ncol,:pver), &
+ coef_q_diff=kvm(:ncol, :)*dpidz_sq(:ncol, :))
+
+ ! ---------------------------------------------------------------------- !
+ ! Calculate 'total' ( tautotx ) and 'tms' ( tautmsx ) stresses that !
+ ! have been actually added into the atmosphere at the current time step. !
+ ! Also, update residual stress, if required. !
+ ! ---------------------------------------------------------------------- !
+
+ do i = 1, ncol
+ ! Compute the implicit 'tms' using the updated winds.
+ ! Below 'tautmsx(i),tautmsy(i)' are pure implicit mountain stresses
+ ! that has been actually added into the atmosphere both for explicit
+ ! and implicit approach.
+
+ tautmsx(i) = -ksrftms(i)*u1(i, pver)
+ tautmsy(i) = -ksrftms(i)*v1(i, pver)
+
+ ! We want to add vertically-integrated Beljaars drag to residual stress.
+ ! So this has to be calculated locally.
+ ! We may want to rethink the residual drag calculation performed here on. (jtb)
+ taubljx(i) = 0._kind_phys
+ taubljy(i) = 0._kind_phys
+ do k = 1, pver
+ taubljx(i) = taubljx(i) + (1._kind_phys/gravit)*dragblj(i, k)*u1(i, k)*p%del(i, k)
+ taubljy(i) = taubljy(i) + (1._kind_phys/gravit)*dragblj(i, k)*v1(i, k)*p%del(i, k)
+ end do
+
+ if (do_iss) then
+ ! Compute vertical integration of final horizontal momentum
+ usum_out(i) = 0._kind_phys
+ vsum_out(i) = 0._kind_phys
+ do k = 1, pver
+ usum_out(i) = usum_out(i) + (1._kind_phys/gravit)*u1(i, k)*p%del(i, k)
+ vsum_out(i) = vsum_out(i) + (1._kind_phys/gravit)*v1(i, k)*p%del(i, k)
+ end do
+
+ ! Compute net stress added into the atmosphere at the current time step.
+ ! Note that the difference between 'usum_in' and 'usum_out' are induced
+ ! by 'explicit residual stress + implicit total stress' for implicit case, while
+ ! by 'explicit normal stress + implicit tms stress' for explicit case.
+ ! Here, 'tautotx(i)' is net stress added into the air at the current time step.
+ tauimpx(i) = (usum_out(i) - usum_in(i))/dt
+ tauimpy(i) = (vsum_out(i) - vsum_in(i))/dt
+
+ tautotx(i) = tauimpx(i)
+ tautoty(i) = tauimpy(i)
+
+ ! Compute residual stress and update if required.
+ ! Note that the total stress we should have added at the current step is
+ ! the sum of 'taux(i) - ksrftms(i)*u(i,pver) + tauresx(i)'.
+ if (itaures) then
+ tauresx(i) = taux(i) + tauresx(i) - tauimpx(i) + tautmsx(i) + taubljx(i)
+ tauresy(i) = tauy(i) + tauresy(i) - tauimpy(i) + tautmsy(i) + taubljy(i)
+ end if
+ else
+ ! .not. do_iss
+ tautotx(i) = taux(i) + tautmsx(i)
+ tautoty(i) = tauy(i) + tautmsy(i)
+ tauresx(i) = 0._kind_phys
+ tauresy(i) = 0._kind_phys
+ end if ! End of 'do_iss' if
+ end do ! End of 'do i = 1, ncol' loop
+
+ ! ------------------------------------ !
+ ! Calculate kinetic energy dissipation !
+ ! ------------------------------------ !
+
+ ! Modification : In future, this should be set exactly same as
+ ! the ones in the convection schemes
+
+ if (do_beljaars) then
+ ! 2. Add Kinetic Energy change across dissipation to Static Energy
+ do k = 1, pver
+ do i = 1, ncol
+ keg_in(i, k) = 0.5_kind_phys*(u0(i, k)*u0(i, k) + v0(i, k)*v0(i, k))
+ keg_out(i, k) = 0.5_kind_phys*(u1(i, k)*u1(i, k) + v1(i, k)*v1(i, k))
+ end do
+ end do
+
+ do k = 1, pver
+ do i = 1, ncol
+ dtk(i, k) = keg_in(i, k) - keg_out(i, k)
+ dse1(i, k) = dse1(i, k) + dtk(i, k) ! + dkeblj(i,k)
+ end do
+ end do
+ else
+ ! .not. do_beljaars
+
+ ! 1. Compute dissipation term at interfaces
+ ! Note that 'u,v' are already diffused wind, and 'tautotx,tautoty' are
+ ! implicit stress that has been actually added. On the other hand,
+ ! 'dinp_u, dinp_v' were computed using non-diffused input wind.
+
+ ! Modification : I should check whether non-consistency between 'u' and 'dinp_u'
+ ! is correctly intended approach. I think so.
+ k = pver + 1
+ do i = 1, ncol
+ tmpi1(i, 1) = 0._kind_phys
+ tmpi1(i, k) = 0.5_kind_phys*dt*gravit* &
+ ((-u1(i, k - 1) + dinp_u(i, k))*tautotx(i) + (-v1(i, k - 1) + dinp_v(i, k))*tautoty(i))
+ end do
+
+ do k = 2, pver
+ do i = 1, ncol
+ dout_u = u1(i, k) - u1(i, k - 1)
+ dout_v = v1(i, k) - v1(i, k - 1)
+ tmpi1(i, k) = 0.25_kind_phys*tmpi2(i, k)*kvm(i, k)* &
+ (dout_u**2 + dout_v**2 + dout_u*dinp_u(i, k) + dout_v*dinp_v(i, k))
+ end do
+ end do
+
+ ! 2. Compute dissipation term at midpoints, add to dry static energy
+ do k = 1, pver
+ do i = 1, ncol
+ dtk(i, k) = (tmpi1(i, k + 1) + tmpi1(i, k))*p%rdel(i, k)
+ dse1(i, k) = dse1(i, k) + dtk(i, k)
+ end do
+ end do
+ end if
+
+ ! End of kinetic energy dissipation due to horizontal momentum diffusion
+ end subroutine vertical_diffusion_diffuse_horizontal_momentum_run
+
+ ! Set dry static energy top boundary condition - zero
+ ! used when (1) no molecular diffusion is active, or (2) molecular diffusion is active and not WACCM-X.
+ !
+ ! Note: this is an artifact to the way the vertical diffusion calculations are done;
+ ! it does not represent physical reality.
+!> \section arg_table_vertical_diffusion_set_dry_static_energy_at_toa_zero_run Argument Table
+!! \htmlinclude arg_table_vertical_diffusion_set_dry_static_energy_at_toa_zero_run.html
+ pure subroutine vertical_diffusion_set_dry_static_energy_at_toa_zero_run( &
+ ncol, &
+ dse_top, &
+ errmsg, errflg)
+
+ ! Input arguments
+ integer, intent(in) :: ncol
+
+ ! Output arguments
+ real(kind_phys), intent(out) :: dse_top(:) ! Dry static energy top boundary condition.
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ errmsg = ''
+ errflg = 0
+
+ dse_top(:ncol) = 0._kind_phys
+
+ end subroutine vertical_diffusion_set_dry_static_energy_at_toa_zero_run
+
+ ! Set dry static energy top boundary condition - molecular diffusion non-WACCM-X version
+!> \section arg_table_vertical_diffusion_set_dry_static_energy_at_toa_molecdiff_run Argument Table
+!! \htmlinclude arg_table_vertical_diffusion_set_dry_static_energy_at_toa_molecdiff_run.html
+ pure subroutine vertical_diffusion_set_dry_static_energy_at_toa_molecdiff_run( &
+ ncol, &
+ gravit, &
+ cpairv, &
+ zi, &
+ tint, &
+ dse_top, &
+ errmsg, errflg)
+
+ ! Input arguments
+ integer, intent(in) :: ncol
+ real(kind_phys), intent(in) :: gravit
+ real(kind_phys), intent(in) :: cpairv(:,:) ! specific heat of dry air at constant pressure [J kg-1 K-1]
+ real(kind_phys), intent(in) :: zi(:,:) ! geopotential_height_wrt_surface_at_interface [m]
+ real(kind_phys), intent(in) :: tint(:,:) ! Air temperature at interfaces [K]
+
+ ! Output arguments
+ real(kind_phys), intent(out) :: dse_top(:) ! Dry static energy top boundary condition.
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ errmsg = ''
+ errflg = 0
+
+ dse_top(:ncol) = cpairv(:ncol,1) * tint(:ncol,1) + gravit * zi(:ncol,1)
+
+ end subroutine vertical_diffusion_set_dry_static_energy_at_toa_molecdiff_run
+
+ ! Driver routine to compute vertical diffusion of dry static energy.
+ ! The new temperature is computed from the diffused dry static energy
+ ! after tendencies are applied.
+ ! Turbulent diffusivities and boundary layer nonlocal transport terms are
+ ! obtained from the turbulence module.
+!> \section arg_table_vertical_diffusion_diffuse_dry_static_energy_run Argument Table
+!! \htmlinclude arg_table_vertical_diffusion_diffuse_dry_static_energy_run.html
+ subroutine vertical_diffusion_diffuse_dry_static_energy_run( &
+ ncol, pver, &
+ dt, &
+ gravit, &
+ p, rhoi, &
+ shflx, &
+ dse_top, &
+ kvh, cgh, &
+ dpidz_sq, &
+ ! input/output (may have been updated by horizontal momentum diffusion)
+ dse, &
+ errmsg, errflg)
+
+ use coords_1d, only: Coords1D
+ use linear_1d_operators, only: BoundaryType, BoundaryFixedLayer, &
+ BoundaryData
+ use vertical_diffusion_solver, only: fin_vol_solve
+
+ ! Input Arguments
+ integer, intent(in) :: ncol
+ integer, intent(in) :: pver
+ real(kind_phys), intent(in) :: dt ! physics timestep [s]
+
+ real(kind_phys), intent(in) :: gravit
+ type(coords1d), intent(in) :: p ! Pressure coordinates [Pa]
+ real(kind_phys), intent(in) :: rhoi(:, :) ! Density of air at interfaces [kg m-3]
+ real(kind_phys), intent(in) :: shflx(:) ! Surface sensible heat flux [W m-2]
+ real(kind_phys), intent(in) :: dse_top(:) ! Dry static energy top boundary condition.
+ real(kind_phys), intent(in) :: kvh(:, :) ! Eddy diffusivity for heat [m^2 s-1], interfaces
+ real(kind_phys), intent(in) :: cgh(:, :) ! Counter-gradient term for heat [J kg-1 m-1], interfaces
+ real(kind_phys), intent(in) :: dpidz_sq(:,:) ! Square of derivative of pressure with height (moist) [kg2 m-4 s-4]
+
+ ! Input/output arguments
+ real(kind_phys), intent(inout) :: dse(:,:) ! After vertical diffusion Dry static energy [J kg-1]
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ ! Local variables
+ integer :: k
+ type(BoundaryType) :: interface_boundary ! Boundary layer objects
+
+ real(kind_phys) :: flux_to_state_conversion_factor(ncol)
+
+ errmsg = ''
+ errflg = 0
+
+ ! Boundary condition for a fixed concentration directly on a boundary
+ ! interface (i.e. a boundary layer of size 0).
+ interface_boundary = BoundaryFixedLayer(spread(0._kind_phys, 1, ncol))
+
+ flux_to_state_conversion_factor(:ncol) = dt*gravit*p%rdel(:, pver)
+
+ ! Modification : In future, we should diffuse the fully conservative
+ ! moist static energy, not the dry static energy.
+
+ ! Add counter-gradient to input static energy profiles
+ do k = 1, pver
+ dse(:ncol, k) = dse(:ncol, k) + dt*p%rdel(:, k)*gravit* &
+ (rhoi(:ncol, k + 1)*kvh(:ncol, k + 1)*cgh(:ncol, k + 1) &
+ - rhoi(:ncol, k)*kvh(:ncol, k)*cgh(:ncol, k))
+ end do
+
+ ! Add the explicit surface fluxes to the lowest layer
+ dse(:ncol, pver) = dse(:ncol, pver) + flux_to_state_conversion_factor(:ncol)*shflx(:ncol)
+
+ !---------------------------------------------------
+ ! Solve for temperature using thermal conductivity
+ !---------------------------------------------------
+
+ ! Boundary layer thickness of "0._kind_phys" signifies that the boundary
+ ! condition is defined directly on the top interface.
+ dse(:ncol, :) = fin_vol_solve(dt, p, dse(:ncol, :), ncol, pver, &
+ coef_q_diff=kvh(:ncol, :)*dpidz_sq(:ncol, :), &
+ upper_bndry=interface_boundary, &
+ l_cond=BoundaryData(dse_top(:ncol)))
+
+ end subroutine vertical_diffusion_diffuse_dry_static_energy_run
+
+ ! Set tracers to diffuse using vertical diffusion.
+!> \section arg_table_vertical_diffusion_diffuse_tracers_init Argument Table
+!! \htmlinclude arg_table_vertical_diffusion_diffuse_tracers_init.html
+ pure subroutine vertical_diffusion_diffuse_tracers_init( &
+ ncnst, &
+ do_diffusion_const, &
+ errmsg, errflg)
+
+ ! Input arguments
+ integer, intent(in) :: ncnst ! # of constituents to diffuse. In eddy_diff, only wv. Others, all advected CCPP constituents.
+
+ ! Output arguments
+ logical, intent(out) :: do_diffusion_const(:) ! diffuse constituents (size ncnst) [flag]
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ errmsg = ''
+ errflg = 0
+
+ ! Check if logical array is of the expected size
+ if (size(do_diffusion_const) .ne. ncnst) then
+ errflg = 1
+ write(errmsg, *) 'compute_vdiff: do_diffusion_const size ', &
+ size(do_diffusion_const), ' is not equal to ncnst ', ncnst
+ return
+ end if
+
+ ! Diffuse all constituents using moist coordinates.
+ do_diffusion_const(:ncnst) = .true.
+
+ ! TODO: in the future (when dropmixnuc is implemented),
+ ! handle that if constituent is treated in dropmixnuc
+ ! (vertically mixed by ndrop activation process) then
+ ! do not handle vertical diffusion in this module.
+
+ end subroutine vertical_diffusion_diffuse_tracers_init
+
+ ! Diffuse tracers (no molecular diffusion)
+!> \section arg_table_vertical_diffusion_diffuse_tracers_run Argument Table
+!! \htmlinclude arg_table_vertical_diffusion_diffuse_tracers_run.html
+ subroutine vertical_diffusion_diffuse_tracers_run( &
+ ncol, pver, &
+ ncnst, &
+ dt, &
+ rair, &
+ gravit, &
+ do_diffusion_const, &
+ p, &
+ t, &
+ rhoi, &
+ cflx, &
+ kvh, kvq, cgs, &
+ qmincg, &
+ dpidz_sq, &
+ ubc_mmr, &
+ cnst_fixed_ubc, &
+ q0, &
+ ! below output
+ q1, &
+ errmsg, errflg)
+
+ use coords_1d, only: Coords1D
+ use linear_1d_operators, only: TriDiagDecomp
+ use vdiff_lu_solver, only: fin_vol_lu_decomp
+
+ ! Input arguments
+ integer, intent(in) :: ncol ! Number of atmospheric columns
+ integer, intent(in) :: pver
+ integer, intent(in) :: ncnst ! # of constituents to diffuse. In eddy_diff, only wv. Others, all advected CCPP constituents.
+ real(kind_phys), intent(in) :: dt ! physics timestep [s]
+ real(kind_phys), intent(in) :: rair
+ real(kind_phys), intent(in) :: gravit
+ logical, intent(in) :: do_diffusion_const(:) ! diffuse constituents (size ncnst) [flag]
+
+ type(coords1d), intent(in) :: p ! Pressure coordinates [Pa]
+ real(kind_phys), intent(in) :: t(:, :) ! Temperature [K]
+ real(kind_phys), intent(in) :: rhoi(:, :) ! Density of air at interfaces [kg m-3]
+
+ real(kind_phys), intent(in) :: cflx(:, :) ! Surface constituent flux [kg m-2 s-1] (ncol,ncnst)
+
+ real(kind_phys), intent(in) :: kvh(:, :) ! Eddy diffusivity for heat [m^2 s-1], interfaces
+ real(kind_phys), intent(in) :: kvq(:, :) ! Eddy diffusivity for constituents [m^2 s-1], interfaces
+ real(kind_phys), intent(in) :: cgs(:, :) ! Counter-gradient star [cg/flux] [s m-2], interfaces
+
+ real(kind_phys), intent(in) :: qmincg(:) ! Minimum constituent mixing ratios from cg fluxes (ncnst)
+ real(kind_phys), intent(in) :: dpidz_sq(:,:) ! Square of derivative of pressure with height (moist) [kg2 m-4 s-4]
+
+ ! Upper boundary properties (for when molecular diffusion is off)
+ real(kind_phys), intent(in) :: ubc_mmr(:, :) ! Upper boundary condition of constituents [none]
+ logical, intent(in) :: cnst_fixed_ubc(:)! Whether upper boundary condition is fixed
+
+ real(kind_phys), intent(in) :: q0(:,:,:) ! Input Constituents [kg kg-1]
+
+ ! Output arguments
+ real(kind_phys), intent(out) :: q1(:,:,:) ! After vertical diffusion Constituents [kg kg-1]
+
+
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ ! Local variables
+ integer :: m, k
+ type(TriDiagDecomp) :: no_molec_decomp ! LU decomposition information.
+ real(kind_phys) :: flux_to_state_conversion_factor(ncol)
+ logical :: lqtst(ncol) ! Adjust vertical profiles
+ real(kind_phys) :: qtm(ncol, pver) ! Temporary copy of q
+ real(kind_phys) :: rrho(ncol) ! 1./bottom level density
+
+ errmsg = ''
+ errflg = 0
+
+ ! Set initial iterative outputs to input values
+ q1(:ncol,:pver,:ncnst) = q0(:ncol,:pver,:ncnst)
+
+ ! necessary temporaries used in computation
+ rrho(:ncol) = rair*t(:ncol, pver)/p%mid(:, pver)
+ flux_to_state_conversion_factor(:ncol) = dt*gravit*p%rdel(:, pver)
+
+ ! Loop through constituents
+ no_molec_decomp = fin_vol_lu_decomp(dt, p, &
+ coef_q_diff=kvq(:ncol, :)*dpidz_sq(:ncol, :))
+
+ do m = 1, ncnst
+ if (do_diffusion_const(m)) then
+ ! Add the nonlocal transport terms to constituents in the PBL.
+ ! Check for neg q's in each constituent and put the original vertical
+ ! profile back if a neg value is found. A neg value implies that the
+ ! quasi-equilibrium conditions assumed for the countergradient term are
+ ! strongly violated.
+ qtm(:ncol, :pver) = q1(:ncol, :pver, m)
+
+ do k = 1, pver
+ q1(:ncol, k, m) = q1(:ncol, k, m) + &
+ dt*p%rdel(:, k)*gravit*(cflx(:ncol, m)*rrho(:ncol))* &
+ (rhoi(:ncol, k + 1)*kvh(:ncol, k + 1)*cgs(:ncol, k + 1) &
+ - rhoi(:ncol, k)*kvh(:ncol, k)*cgs(:ncol, k))
+ end do
+ lqtst(:ncol) = all(q1(:ncol, 1:pver, m) >= qmincg(m), 2)
+ do k = 1, pver
+ q1(:ncol, k, m) = merge(q1(:ncol, k, m), qtm(:ncol, k), lqtst(:ncol))
+ end do
+
+ ! Add the explicit surface fluxes to the lowest layer
+ q1(:ncol, pver, m) = q1(:ncol, pver, m) + flux_to_state_conversion_factor(:ncol)*cflx(:ncol, m)
+
+ ! not doing molecular diffusion
+ ! explicitly set mmr in top layer for cases where molecular diffusion is not active
+ if (cnst_fixed_ubc(m)) then
+ ! assumes TOA is at vertical level 1 here.
+ q1(:ncol, 1, m) = ubc_mmr(:ncol, m)
+ end if
+ call no_molec_decomp%left_div(q1(:ncol, :, m))
+ end if
+ end do
+
+ call no_molec_decomp%finalize()
+
+ end subroutine vertical_diffusion_diffuse_tracers_run
+
+ ! Based on provisional output from vertical diffusion,
+ ! calculate final physics tendencies to be applied by the tendency updater.
+!> \section arg_table_vertical_diffusion_tendencies_run Argument Table
+!! \htmlinclude arg_table_vertical_diffusion_tendencies_run.html
+ subroutine vertical_diffusion_tendencies_run( &
+ ncol, pver, pcnst, &
+ const_props, &
+ dt, &
+ pdel, pdeldry, &
+ u0, v0, s0, q0, & ! actual values at beginning of vdiff
+ u1, v1, s1, q1, & ! provisional values after vdiff, not actual
+ ! below output
+ tend_s, tend_u, tend_v, tend_q, &
+ scheme_name, &
+ errmsg, errflg)
+
+ ! framework dependency for const_props
+ use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
+
+ ! Input arguments
+ integer, intent(in) :: ncol ! Number of atmospheric columns [count]
+ integer, intent(in) :: pver ! Number of vertical levels [count]
+ integer, intent(in) :: pcnst ! Number of ccpp constituents [count]
+ type(ccpp_constituent_prop_ptr_t), &
+ intent(in) :: const_props(:) ! CCPP constituent properties pointer
+ real(kind_phys), intent(in) :: dt ! Timestep [s]
+ real(kind_phys), intent(in) :: pdel(:,:) ! Pressure thickness of layers [Pa]
+ real(kind_phys), intent(in) :: pdeldry(:,:) ! Dry pressure thickness of layers [Pa]
+ real(kind_phys), intent(in) :: u0(:,:) ! Initial eastward wind [m s-1]
+ real(kind_phys), intent(in) :: v0(:,:) ! Initial northward wind [m s-1]
+ real(kind_phys), intent(in) :: s0(:,:) ! Initial dry static energy [J kg-1]
+ real(kind_phys), intent(in) :: q0(:,:,:) ! Initial constituent mixing ratios [kg kg-1]
+ real(kind_phys), intent(in) :: u1(:,:) ! Provisional eastward wind after diffusion [m s-1]
+ real(kind_phys), intent(in) :: v1(:,:) ! Provisional northward wind after diffusion [m s-1]
+ real(kind_phys), intent(in) :: s1(:,:) ! Provisional dry static energy after diffusion [J kg-1]
+ real(kind_phys), intent(in) :: q1(:,:,:) ! Provisional constituent mixing ratios after diffusion [kg kg-1]
+
+ ! Output arguments
+ real(kind_phys), intent(out) :: tend_s(:,:) ! Tendency of dry air enthalpy at constant pressure [J kg-1 s-1]
+ real(kind_phys), intent(out) :: tend_u(:,:) ! Eastward wind tendency [m s-2]
+ real(kind_phys), intent(out) :: tend_v(:,:) ! Northward wind tendency [m s-2]
+ real(kind_phys), intent(out) :: tend_q(:,:,:) ! Constituent mixing ratio tendencies [kg kg-1 s-1]
+ character(len=64), intent(out) :: scheme_name ! scheme name
+ character(len=512), intent(out) :: errmsg ! Error message
+ integer, intent(out) :: errflg ! Error flag
+
+ integer :: m
+ logical :: const_is_dry
+
+ ! for bit-to-bitness with CAM.
+ real(kind_phys) :: rdt
+
+ scheme_name = "vertical_diffusion_tendencies"
+
+ rdt = 1._kind_phys/dt
+
+ errmsg = ''
+ errflg = 0
+
+ ! calculate physics tendencies
+ tend_s(:ncol,:) = (s1(:ncol,:) - s0(:ncol,:)) * rdt
+ tend_u(:ncol,:) = (u1(:ncol,:) - u0(:ncol,:)) * rdt
+ tend_v(:ncol,:) = (v1(:ncol,:) - v0(:ncol,:)) * rdt
+ tend_q(:ncol,:pver,:) = (q1(:ncol,:pver,:) - q0(:ncol,:pver,:)) * rdt
+
+ ! convert tendencies of dry constituents to dry basis
+ do m = 1, pcnst
+ call const_props(m)%is_dry(const_is_dry, errflg, errmsg)
+ if(const_is_dry) then
+ tend_q(:ncol,:pver,m) = tend_q(:ncol,:pver,m) * pdel(:ncol,:pver) / pdeldry(:ncol,:pver)
+ endif
+ enddo
+
+ end subroutine vertical_diffusion_tendencies_run
+
+end module diffusion_solver
diff --git a/schemes/vertical_diffusion/diffusion_solver.meta b/schemes/vertical_diffusion/diffusion_solver.meta
new file mode 100644
index 00000000..1063cf92
--- /dev/null
+++ b/schemes/vertical_diffusion/diffusion_solver.meta
@@ -0,0 +1,1045 @@
+[ccpp-table-properties]
+ name = vertical_diffusion_set_temperature_at_toa_default
+ type = scheme
+
+[ccpp-arg-table]
+ name = vertical_diffusion_set_temperature_at_toa_default_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ t ]
+ standard_name = air_temperature
+ units = K
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ t_toai ]
+ standard_name = air_temperature_at_interface_above_top_of_atmosphere_for_vertical_diffusion
+ units = K
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-table-properties]
+ name = vertical_diffusion_interpolate_to_interfaces
+ type = scheme
+
+[ccpp-arg-table]
+ name = vertical_diffusion_interpolate_to_interfaces_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pverp ]
+ standard_name = vertical_interface_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ gravit ]
+ standard_name = standard_gravitational_acceleration
+ units = m s-2
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ rair ]
+ standard_name = gas_constant_of_dry_air
+ units = J kg-1 K-1
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ rairv ]
+ standard_name = composition_dependent_gas_constant_of_dry_air
+ units = J kg-1 K-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ flag_for_constituent_dependent_gas_constant ]
+ standard_name = use_composition_dependent_gas_constant_in_vertical_diffusion
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = in
+[ t ]
+ standard_name = air_temperature
+ units = K
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ t_toai ]
+ standard_name = air_temperature_at_interface_above_top_of_atmosphere_for_vertical_diffusion
+ units = K
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ pint ]
+ standard_name = air_pressure_at_interface
+ units = Pa
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ pintdry ]
+ standard_name = air_pressure_of_dry_air_at_interface
+ units = Pa
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ ti ]
+ standard_name = air_temperature_at_interfaces_for_vertical_diffusion
+ units = K
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = out
+[ rairi ]
+ standard_name = gas_constant_of_dry_air_at_interfaces_for_vertical_diffusion
+ units = J K-1 kg-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = out
+[ rhoi ]
+ standard_name = air_density_at_interfaces_for_vertical_diffusion
+ units = kg m-3
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = out
+[ rhoi_dry ]
+ standard_name = dry_air_density_at_interfaces_for_vertical_diffusion
+ units = kg m-3
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = out
+[ dpidz_sq ]
+ standard_name = square_of_derivative_of_air_pressure_with_respect_to_height_at_interfaces_for_vertical_diffusion
+ units = kg2 m-4 s-4
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-table-properties]
+ name = implicit_surface_stress_add_drag_coefficient
+ type = scheme
+
+[ccpp-arg-table]
+ name = implicit_surface_stress_add_drag_coefficient_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ do_iss ]
+ standard_name = do_implicit_total_surface_stress_in_vertical_diffusion
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = in
+[ taux ]
+ standard_name = eastward_stress_at_surface_for_vertical_diffusion
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ tauy ]
+ standard_name = northward_stress_at_surface_for_vertical_diffusion
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ u0 ]
+ standard_name = eastward_wind
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ v0 ]
+ standard_name = northward_wind
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ ksrf ]
+ standard_name = total_drag_coefficient_at_surface
+ units = kg m-2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = inout
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-table-properties]
+ name = vertical_diffusion_wind_damping_rate
+ type = scheme
+
+[ccpp-arg-table]
+ name = vertical_diffusion_wind_damping_rate_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ gravit ]
+ standard_name = standard_gravitational_acceleration
+ units = m s-2
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ p ]
+ standard_name = vertical_moist_pressure_coordinates_for_vertical_diffusion
+ units = none
+ type = coords1d
+ dimensions = ()
+ intent = in
+[ ksrf ]
+ standard_name = total_drag_coefficient_at_surface
+ units = kg m-2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ tau_damp_rate ]
+ standard_name = horizontal_wind_damping_rate_due_to_total_surface_drag
+ units = s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-table-properties]
+ name = vertical_diffusion_diffuse_horizontal_momentum
+ type = scheme
+
+[ccpp-arg-table]
+ name = vertical_diffusion_diffuse_horizontal_momentum_timestep_init
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ ksrf ]
+ standard_name = total_drag_coefficient_at_surface
+ units = kg m-2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_dimension)
+ intent = out
+[ tau_damp_rate ]
+ standard_name = horizontal_wind_damping_rate_due_to_total_surface_drag
+ units = s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_dimension, vertical_layer_dimension)
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-arg-table]
+ name = vertical_diffusion_diffuse_horizontal_momentum_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pverp ]
+ standard_name = vertical_interface_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ dt ]
+ standard_name = timestep_for_physics
+ units = s
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ rair ]
+ standard_name = gas_constant_of_dry_air
+ units = J kg-1 K-1
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ gravit ]
+ standard_name = standard_gravitational_acceleration
+ units = m s-2
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ do_iss ]
+ standard_name = do_implicit_total_surface_stress_in_vertical_diffusion
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = in
+[ am_correction ]
+ standard_name = do_angular_momentum_correction_in_vertical_diffusion
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = in
+[ itaures ]
+ standard_name = update_residual_stress_at_surface_in_vertical_diffusion
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = in
+[ do_beljaars ]
+ standard_name = do_beljaars_drag_in_vertical_diffusion_tbd
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = in
+[ p ]
+ standard_name = vertical_moist_pressure_coordinates_for_vertical_diffusion
+ units = none
+ type = coords1d
+ dimensions = ()
+ intent = in
+[ t ]
+ standard_name = air_temperature
+ units = K
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ rhoi ]
+ standard_name = air_density_at_interfaces_for_vertical_diffusion
+ units = kg m-3
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ taux ]
+ standard_name = eastward_stress_at_surface_for_vertical_diffusion
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ tauy ]
+ standard_name = northward_stress_at_surface_for_vertical_diffusion
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ tau_damp_rate ]
+ standard_name = horizontal_wind_damping_rate_due_to_total_surface_drag
+ units = s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ ksrftms ]
+ standard_name = turbulent_orographic_form_drag_coefficent_at_surface
+ units = kg s-1 m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ dragblj ]
+ standard_name = turbulent_orographic_form_drag_coefficent
+ units = s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ kvm ]
+ standard_name = eddy_momentum_diffusivity_at_interfaces
+ units = m2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ dpidz_sq ]
+ standard_name = square_of_derivative_of_air_pressure_with_respect_to_height_at_interfaces_for_vertical_diffusion
+ units = kg2 m-4 s-4
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ u0 ]
+ standard_name = eastward_wind
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ v0 ]
+ standard_name = northward_wind
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ dse0 ]
+ standard_name = dry_static_energy
+ units = J kg-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ tauresx ]
+ standard_name = eastward_reserved_stress_at_surface_on_previous_timestep
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = inout
+[ tauresy ]
+ standard_name = northward_reserved_stress_at_surface_on_previous_timestep
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = inout
+[ dtk ]
+ standard_name = kinetic_energy_dissipation_in_atmosphere_boundary_layer
+ units = J kg-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = out
+[ tautmsx ]
+ standard_name = eastward_turbulent_mountain_surface_stress_tbd
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ tautmsy ]
+ standard_name = northward_turbulent_mountain_surface_stress_tbd
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ u1 ]
+ standard_name = eastward_wind_after_vertical_diffusion
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = out
+[ v1 ]
+ standard_name = northward_wind_after_vertical_diffusion
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = out
+[ dse1 ]
+ standard_name = dry_static_energy_after_vertical_diffusion
+ units = J kg-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-table-properties]
+ name = vertical_diffusion_set_dry_static_energy_at_toa_zero
+ type = scheme
+
+[ccpp-arg-table]
+ name = vertical_diffusion_set_dry_static_energy_at_toa_zero_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ dse_top ]
+ standard_name = dry_static_energy_at_top_of_atmosphere_tbd
+ units = J kg-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-table-properties]
+ name = vertical_diffusion_set_dry_static_energy_at_toa_molecdiff
+ type = scheme
+
+[ccpp-arg-table]
+ name = vertical_diffusion_set_dry_static_energy_at_toa_molecdiff_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ gravit ]
+ standard_name = standard_gravitational_acceleration
+ units = m s-2
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ cpairv ]
+ standard_name = composition_dependent_specific_heat_of_dry_air_at_constant_pressure
+ units = J kg-1 K-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ zi ]
+ standard_name = geopotential_height_wrt_surface_at_interface
+ units = m
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ tint ]
+ standard_name = air_temperature_at_interfaces_for_vertical_diffusion
+ units = K
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ dse_top ]
+ standard_name = dry_static_energy_at_top_of_atmosphere_tbd
+ units = J kg-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-table-properties]
+ name = vertical_diffusion_diffuse_dry_static_energy
+ type = scheme
+
+[ccpp-arg-table]
+ name = vertical_diffusion_diffuse_dry_static_energy_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ dt ]
+ standard_name = timestep_for_physics
+ units = s
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ gravit ]
+ standard_name = standard_gravitational_acceleration
+ units = m s-2
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ p ]
+ standard_name = vertical_moist_pressure_coordinates_for_vertical_diffusion
+ units = none
+ type = coords1d
+ dimensions = ()
+ intent = in
+[ rhoi ]
+ standard_name = air_density_at_interfaces_for_vertical_diffusion
+ units = kg m-3
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ shflx ]
+ standard_name = surface_upward_sensible_heat_flux_for_vertical_diffusion
+ units = W m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ dse_top ]
+ standard_name = dry_static_energy_at_top_of_atmosphere_tbd
+ units = J kg-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ kvh ]
+ standard_name = eddy_heat_diffusivity_at_interfaces
+ units = m2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ cgh ]
+ standard_name = counter_gradient_coefficient_for_vertical_diffusion_of_heat_at_interfaces
+ units = J kg-1 m-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ dpidz_sq ]
+ standard_name = square_of_derivative_of_air_pressure_with_respect_to_height_at_interfaces_for_vertical_diffusion
+ units = kg2 m-4 s-4
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ dse ]
+ # this is inout as previously horizontal momentum diffusion already
+ # modified the provisional after-vdiff dse.
+ standard_name = dry_static_energy_after_vertical_diffusion
+ units = J kg-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = inout
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-table-properties]
+ name = vertical_diffusion_diffuse_tracers
+ type = scheme
+
+
+[ccpp-arg-table]
+ name = vertical_diffusion_diffuse_tracers_init
+ type = scheme
+[ ncnst ]
+ # note: this is set to ccpp_constituents for now because
+ # q is sized number_of_ccpp_constituents as well, and cannot
+ # yet be a custom dimension. when CAM5 is ported, ncnst = 1,
+ # and q would have to have a custom dimension of, e.g.,
+ # number_of_ccpp_constituents_for_vertical_diffusion_solver.
+ # but this is not yet supported, so ncnst is all constituents for now.
+ standard_name = number_of_ccpp_constituents
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ do_diffusion_const ]
+ # only dimension should be ncnst, not pcnst.
+ standard_name = flag_for_ccpp_constituent_vertical_diffusion
+ units = flag
+ type = logical
+ dimensions = (number_of_ccpp_constituents)
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-arg-table]
+ name = vertical_diffusion_diffuse_tracers_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ ncnst ]
+ # note: this is set to ccpp_constituents for now because
+ # q is sized number_of_ccpp_constituents as well, and cannot
+ # yet be a custom dimension. when CAM5 is ported, ncnst = 1,
+ # and q would have to have a custom dimension of, e.g.,
+ # number_of_ccpp_constituents_for_vertical_diffusion_solver.
+ # but this is not yet supported, so ncnst is all constituents for now.
+ standard_name = number_of_ccpp_constituents
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ dt ]
+ standard_name = timestep_for_physics
+ units = s
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ rair ]
+ standard_name = gas_constant_of_dry_air
+ units = J kg-1 K-1
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ gravit ]
+ standard_name = standard_gravitational_acceleration
+ units = m s-2
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ do_diffusion_const ]
+ # only dimension should be ncnst, not pcnst.
+ standard_name = flag_for_ccpp_constituent_vertical_diffusion
+ units = flag
+ type = logical
+ dimensions = (number_of_ccpp_constituents)
+ intent = in
+[ p ]
+ standard_name = vertical_moist_pressure_coordinates_for_vertical_diffusion
+ units = none
+ type = coords1d
+ dimensions = ()
+ intent = in
+[ t ]
+ standard_name = air_temperature
+ units = K
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ rhoi ]
+ standard_name = air_density_at_interfaces_for_vertical_diffusion
+ units = kg m-3
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ cflx ]
+ # same note applies for the ncnst argument above that the second dimension
+ # is actually ncnst, not pcnst.
+ #
+ # when this is made ncnst, the standard name could be changed to
+ # surface_upward_ccpp_constituent_fluxes_for_vertical_diffusion
+ standard_name = surface_upward_ccpp_constituent_fluxes_from_coupler
+ units = kg m-2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, number_of_ccpp_constituents)
+ intent = in
+[ kvh ]
+ standard_name = eddy_heat_diffusivity_at_interfaces
+ units = m2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ kvq ]
+ standard_name = eddy_constituent_diffusivity_at_interfaces
+ units = m2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ cgs ]
+ standard_name = counter_gradient_coefficient_for_vertical_diffusion_of_constituents_at_interfaces
+ units = s m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ qmincg ]
+ # same note applies for the ncnst argument above that the only dimension
+ # is actually ncnst, not pcnst.
+ standard_name = ccpp_constituent_minimum_values
+ units = kg kg-1
+ type = real | kind = kind_phys
+ dimensions = (number_of_ccpp_constituents)
+ intent = in
+[ dpidz_sq ]
+ standard_name = square_of_derivative_of_air_pressure_with_respect_to_height_at_interfaces_for_vertical_diffusion
+ units = kg2 m-4 s-4
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ intent = in
+[ ubc_mmr ]
+ # same note applies for the ncnst argument above that the second dimension
+ # is actually ncnst, not pcnst.
+ standard_name = upper_boundary_condition_of_ccpp_constituents
+ units = none
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, number_of_ccpp_constituents)
+ intent = in
+[ cnst_fixed_ubc ]
+ # same note applies for the ncnst argument above that the only dimension
+ # is actually ncnst, not pcnst.
+ standard_name = use_fixed_upper_boundary_condition_of_ccpp_constituents_tbd
+ units = flag
+ type = logical
+ dimensions = (number_of_ccpp_constituents)
+ intent = in
+[ q0 ]
+ standard_name = ccpp_constituents
+ units = none
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccpp_constituents)
+ intent = in
+[ q1 ]
+ standard_name = ccpp_constituents_after_vertical_diffusion
+ units = none
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccpp_constituents)
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-table-properties]
+ name = vertical_diffusion_tendencies
+ type = scheme
+
+[ccpp-arg-table]
+ name = vertical_diffusion_tendencies_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pcnst ]
+ standard_name = number_of_ccpp_constituents
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ const_props ]
+ standard_name = ccpp_constituent_properties
+ units = none
+ type = ccpp_constituent_prop_ptr_t
+ dimensions = (number_of_ccpp_constituents)
+ intent = in
+[ dt ]
+ standard_name = timestep_for_physics
+ units = s
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ pdel ]
+ standard_name = air_pressure_thickness
+ units = Pa
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ pdeldry ]
+ standard_name = air_pressure_thickness_of_dry_air
+ units = Pa
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ u0 ]
+ standard_name = eastward_wind
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ v0 ]
+ standard_name = northward_wind
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ s0 ]
+ standard_name = dry_static_energy
+ units = J kg-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ q0 ]
+ standard_name = ccpp_constituents
+ units = none
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccpp_constituents)
+ intent = in
+[ u1 ]
+ standard_name = eastward_wind_after_vertical_diffusion
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ v1 ]
+ standard_name = northward_wind_after_vertical_diffusion
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ s1 ]
+ standard_name = dry_static_energy_after_vertical_diffusion
+ units = J kg-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ q1 ]
+ standard_name = ccpp_constituents_after_vertical_diffusion
+ units = none
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccpp_constituents)
+ intent = in
+[ tend_s ]
+ standard_name = tendency_of_dry_air_enthalpy_at_constant_pressure
+ units = J kg-1 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = out
+[ tend_u ]
+ standard_name = tendency_of_eastward_wind
+ units = m s-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = out
+[ tend_v ]
+ standard_name = tendency_of_northward_wind
+ units = m s-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = out
+[ tend_q ]
+ standard_name = ccpp_constituent_tendencies
+ units = none
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccpp_constituents)
+ intent = out
+[ scheme_name ]
+ standard_name = scheme_name
+ units = none
+ type = character | kind = len=64
+ dimensions = ()
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
diff --git a/schemes/vertical_diffusion/diffusion_stubs.F90 b/schemes/vertical_diffusion/diffusion_stubs.F90
new file mode 100644
index 00000000..9f5357d3
--- /dev/null
+++ b/schemes/vertical_diffusion/diffusion_stubs.F90
@@ -0,0 +1,332 @@
+! Stubs for schemes that have not been CCPPized but their interfaces
+! to vertical diffusion have been.
+! These schemes will eventually be moved to the corresponding scheme
+! directory after their actual computations have been CCPPized.
+!
+! CCPPized: Haipeng Lin, May 2025
+module diffusion_stubs
+ use ccpp_kinds, only: kind_phys
+
+ implicit none
+ private
+ save
+
+ ! CCPP-compliant subroutines
+ public :: zero_upper_boundary_condition_init
+
+ public :: tms_beljaars_zero_stub_run
+ public :: turbulent_mountain_stress_add_drag_coefficient_run
+ public :: beljaars_add_wind_damping_rate_run
+
+ public :: beljaars_add_updated_residual_stress_run
+ public :: turbulent_mountain_stress_add_updated_surface_stress_run
+
+ public :: vertical_diffusion_not_use_rairv_init
+
+contains
+
+ ! Stub for zero upper boundary conditions before UBC module is CCPPized
+ !
+ ! NOTE: when UBC is CCPPized, the units for ubc_mmr will depend on the units of the
+ ! ccpp constituents (host model dependent) and not necessarily kg kg-1 mass mixing ratio;
+ ! the upper boundary condition is applied directly by overwriting the post-diffusion value
+ ! of the constituents array at the upper boundary of the model, q1(:ncol,1,m), with the value
+ ! in ubc_mmr. Thus care needs to be taken to make sure the units match; at the CCPP level,
+ ! the units specified here for ubc_mmr are "none", same as the units in the ccpp_constituents
+ ! array.
+!> \section arg_table_zero_upper_boundary_condition_init Argument Table
+!! \htmlinclude zero_upper_boundary_condition_init.html
+ subroutine zero_upper_boundary_condition_init( &
+ ncol, pcnst, &
+ ! below output
+ ubc_mmr, &
+ cnst_fixed_ubc, &
+ errmsg, errflg)
+
+ ! Input arguments
+ integer, intent(in) :: ncol
+ integer, intent(in) :: pcnst
+
+ ! Output arguments
+ real(kind_phys), intent(out) :: ubc_mmr(:,:) ! Upper boundary condition mass mixing ratios [none]
+ logical, intent(out) :: cnst_fixed_ubc(:) ! Flag for fixed upper boundary condition of constituents [flag]
+ character(len=512), intent(out) :: errmsg ! Error message
+ integer, intent(out) :: errflg ! Error flag
+
+ errmsg = ''
+ errflg = 0
+
+ ! Set all upper boundary condition mixing ratios to zero
+ ubc_mmr(:ncol, :pcnst) = 0._kind_phys
+
+ ! Set all fixed upper boundary condition flags to false
+ cnst_fixed_ubc(:pcnst) = .false.
+
+ end subroutine zero_upper_boundary_condition_init
+
+ ! Stub for TMS/Beljaars to be set to zero while they are not implemented.
+!> \section arg_table_tms_beljaars_zero_stub_run Argument Table
+!! \htmlinclude tms_beljaars_zero_stub_run.html
+ subroutine tms_beljaars_zero_stub_run( &
+ ncol, pver, &
+ ksrftms, &
+ tautmsx, tautmsy, &
+ do_beljaars, &
+ dragblj, &
+ taubljx, taubljy, &
+ errmsg, errflg)
+
+ ! Input arguments
+ integer, intent(in) :: ncol
+ integer, intent(in) :: pver
+
+ ! Output arguments
+ real(kind_phys), intent(out) :: ksrftms(:) ! Surface drag coefficient for turbulent mountain stress. > 0. [kg m-2 s-1]
+ real(kind_phys), intent(out) :: tautmsx(:) ! Eastward turbulent mountain surface stress [N m-2]
+ real(kind_phys), intent(out) :: tautmsy(:) ! Northward turbulent mountain surface stress [N m-2]
+ logical, intent(out) :: do_beljaars
+ real(kind_phys), intent(out) :: dragblj(:,:)! Drag profile from Beljaars SGO form drag > 0. [s-1]
+ real(kind_phys), intent(out) :: taubljx(:) ! Eastward Beljaars surface stress [N m-2]
+ real(kind_phys), intent(out) :: taubljy(:) ! Northward Beljaars surface stress [N m-2]
+ character(len=512), intent(out) :: errmsg ! Error message
+ integer, intent(out) :: errflg ! Error flag
+
+ errmsg = ''
+ errflg = 0
+
+ ! Set TMS drag coefficient to zero (stub implementation)
+ ksrftms(:ncol) = 0._kind_phys
+
+ ! Set all TMS and Beljaars stresses to zero (stub implementation)
+ tautmsx(:ncol) = 0._kind_phys
+ tautmsy(:ncol) = 0._kind_phys
+ taubljx(:ncol) = 0._kind_phys
+ taubljy(:ncol) = 0._kind_phys
+
+ ! Set Beljaars 3-D drag profile to zero (stub implementation)
+ dragblj(:ncol,:pver) = 0._kind_phys
+
+ ! Set do_beljaars flag to false
+ do_beljaars = .false.
+
+ end subroutine tms_beljaars_zero_stub_run
+
+ ! Add turbulent mountain stress to the total surface drag coefficient
+!> \section arg_table_turbulent_mountain_stress_add_drag_coefficient_run Argument Table
+!! \htmlinclude arg_table_turbulent_mountain_stress_add_drag_coefficient_run.html
+ subroutine turbulent_mountain_stress_add_drag_coefficient_run( &
+ ncol, pver, &
+ ksrftms, &
+ ksrf, &
+ errmsg, errflg)
+
+ ! Input arguments
+ integer, intent(in) :: ncol
+ integer, intent(in) :: pver
+ real(kind_phys), intent(in) :: ksrftms(:) ! Surface drag coefficient for turbulent mountain stress. > 0. [kg m-2 s-1]
+
+ ! Input/output arguments
+ real(kind_phys), intent(inout) :: ksrf(:) ! total surface drag coefficient [kg m-2 s-1]
+
+ ! Output arguments
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ errmsg = ''
+ errflg = 0
+
+ ksrf(:ncol) = ksrf(:ncol) + ksrftms(:ncol)
+
+ end subroutine turbulent_mountain_stress_add_drag_coefficient_run
+
+ ! Add Beljaars drag to the wind damping rate for vertical diffusion.
+ ! Has to run after vertical_diffusion_wind_damping_rate
+!> \section arg_table_beljaars_add_wind_damping_rate_run Argument Table
+!! \htmlinclude arg_table_beljaars_add_wind_damping_rate_run.html
+ subroutine beljaars_add_wind_damping_rate_run( &
+ ncol, pver, &
+ dragblj, &
+ tau_damp_rate, &
+ errmsg, errflg)
+
+ ! Input arguments
+ integer, intent(in) :: ncol
+ integer, intent(in) :: pver
+ real(kind_phys), intent(in) :: dragblj(:, :) ! Drag profile from Beljaars SGO form drag > 0. [s-1]
+
+ ! Input/output arguments
+ real(kind_phys), intent(inout) :: tau_damp_rate(:,:) ! Rate at which external (surface) stress damps wind speeds [s-1]
+
+ ! Output arguments
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ errmsg = ''
+ errflg = 0
+
+ ! Beljaars et al SGO scheme incorporated here. It
+ ! appears as a "3D" tau_damp_rate specification.
+ tau_damp_rate(:ncol,:pver) = tau_damp_rate(:ncol,:pver) + dragblj(:ncol,:pver)
+
+ end subroutine beljaars_add_wind_damping_rate_run
+
+ ! Add Beljaars stress using updated provisional winds to surface stresses used
+ ! for horizontal momentum diffusion - kinetic energy dissipation.
+!> \section arg_table_beljaars_add_updated_residual_stress_run Argument Table
+!! \htmlinclude arg_table_beljaars_add_updated_residual_stress_run.html
+ subroutine beljaars_add_updated_residual_stress_run( &
+ ncol, pver, &
+ gravit, &
+ p, &
+ do_iss, &
+ itaures, &
+ dragblj, &
+ u1, v1, & ! provisionally updated winds
+ ! input/output
+ tauresx, tauresy, &
+ errmsg, errflg)
+
+ use coords_1d, only: Coords1D
+
+ ! Input arguments
+ integer, intent(in) :: ncol
+ integer, intent(in) :: pver
+ real(kind_phys), intent(in) :: gravit
+ type(coords1d), intent(in) :: p ! Pressure coordinates [Pa]
+ logical, intent(in) :: do_iss ! Flag for implicit surface stress (namelist from vdiff)
+ logical, intent(in) :: itaures ! Flag for updating tauresx tauresy in this subroutine.
+ real(kind_phys), intent(in) :: dragblj(:, :) ! Drag profile from Beljaars SGO form drag > 0. [s-1]
+ real(kind_phys), intent(in) :: u1(:,:) ! After vertical diffusion u-wind [m s-1]
+ real(kind_phys), intent(in) :: v1(:,:) ! After vertical diffusion v-wind [m s-1]
+
+ ! Input/output arguments
+ real(kind_phys), intent(inout) :: tauresx(:) ! Partially updated residual surface stress
+ real(kind_phys), intent(inout) :: tauresy(:)
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ integer :: i, k
+ real(kind_phys) :: taubljx(ncol) ! recomputed explicit/residual beljaars stress
+ real(kind_phys) :: taubljy(ncol) ! recomputed explicit/residual beljaars stress
+
+ errmsg = ''
+ errflg = 0
+
+ if(.not. do_iss) then
+ ! Not added to residual stress if implicit surface stress is off
+ return
+ endif
+
+ if(.not. itaures) then
+ ! Not added to residual stress if residual stress is not requested to be updated.
+ return
+ endif
+
+ do i = 1, ncol
+ ! Add vertically-integrated Beljaars drag to residual stress.
+ ! note: in vertical_diffusion_tend, tautotx has taubljx added to it
+ ! but this is computed separately using UPDATED u, v (and is not written back to pbuf)
+ ! there is a FIXME in beljaars_drag that notes maybe updated u, v could be used there
+ ! but that is too early (before vdiff). Keeping this for now. hplin 5/22/25
+ taubljx(i) = 0._kind_phys
+ taubljy(i) = 0._kind_phys
+ do k = 1, pver
+ taubljx(i) = taubljx(i) + (1._kind_phys/gravit)*dragblj(i, k)*u1(i, k)*p%del(i, k)
+ taubljy(i) = taubljy(i) + (1._kind_phys/gravit)*dragblj(i, k)*v1(i, k)*p%del(i, k)
+ end do
+ enddo
+
+ tauresx(:ncol) = tauresx(:ncol) + taubljx(:ncol)
+ tauresy(:ncol) = tauresy(:ncol) + taubljy(:ncol)
+
+ end subroutine beljaars_add_updated_residual_stress_run
+
+ ! Add TMS using updated provisional winds to total surface stress in
+ ! horizontal momentum diffusion - kinetic energy dissipation.
+!> \section arg_table_turbulent_mountain_stress_add_updated_surface_stress_run Argument Table
+!! \htmlinclude arg_table_turbulent_mountain_stress_add_updated_surface_stress_run.html
+ subroutine turbulent_mountain_stress_add_updated_surface_stress_run( &
+ ncol, pver, &
+ do_iss, itaures, &
+ ksrftms, &
+ u1, v1, & ! provisionally updated winds
+ ! input/output
+ tauresx, tauresy, &
+ tautotx, tautoty, &
+ ! output
+ tautmsx, tautmsy, &
+ errmsg, errflg)
+
+ ! Input Arguments
+ integer, intent(in) :: ncol ! Number of atmospheric columns
+ integer, intent(in) :: pver
+ logical, intent(in) :: do_iss ! Flag for implicit surface stress (namelist from vdiff)
+ logical, intent(in) :: itaures ! Flag for updating tauresx tauresy in this subroutine.
+ real(kind_phys), intent(in) :: ksrftms(:) ! Surface drag coefficient for turbulent mountain stress. > 0. [kg m-2 s-1]
+ real(kind_phys), intent(in) :: u1(:,:) ! After vertical diffusion u-wind [m s-1]
+ real(kind_phys), intent(in) :: v1(:,:) ! After vertical diffusion v-wind [m s-1]
+
+ ! Input/output arguments
+ real(kind_phys), intent(inout) :: tauresx(:) ! Partially updated residual surface stress
+ real(kind_phys), intent(inout) :: tauresy(:)
+ real(kind_phys), intent(inout) :: tautotx(:) ! Total surface stress (zonal)
+ real(kind_phys), intent(inout) :: tautoty(:) ! Total surface stress (meridional)
+
+ ! Output arguments
+ real(kind_phys), intent(out) :: tautmsx(:) ! Implicit zonal turbulent mountain surface stress [N m-2]
+ real(kind_phys), intent(out) :: tautmsy(:) ! Implicit meridional turbulent mountain surface stress [N m-2]
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ errmsg = ''
+ errflg = 0
+
+ ! Compute the implicit 'tms' using the updated winds.
+ ! Below 'tautmsx(i),tautmsy(i)' are pure implicit mountain stresses
+ ! that has been actually added into the atmosphere both for explicit
+ ! and implicit approach.
+ tautmsx(:ncol) = -ksrftms(:ncol)*u1(:ncol, pver)
+ tautmsy(:ncol) = -ksrftms(:ncol)*v1(:ncol, pver)
+
+ if(do_iss) then
+ ! Implicit surface stress:
+ if(.not. itaures) then
+ ! Not added to residual stress if residual stress is not requested to be updated.
+ return
+ endif
+
+ tauresx(:ncol) = tauresx(:ncol) + tautmsx(:ncol)
+ tauresy(:ncol) = tauresy(:ncol) + tautmsy(:ncol)
+
+ ! leave tautotx, tautoty unchanged (they are not updated by TMS if do_iss.)
+ else ! .not. do_iss
+ ! Not implicit surface stress. Add to total surface stress
+ tautotx(:ncol) = tautotx(:ncol) + tautmsx(:ncol)
+ tautoty(:ncol) = tautoty(:ncol) + tautmsy(:ncol)
+ endif
+ end subroutine turbulent_mountain_stress_add_updated_surface_stress_run
+
+ ! Set rairv use flag for vertical_diffusion_interpolate_to_interfaces
+ ! to be false since we are not in WACCM-X mode.
+!> \section arg_table_vertical_diffusion_not_use_rairv_init Argument Table
+!! \htmlinclude vertical_diffusion_not_use_rairv_init.html
+ subroutine vertical_diffusion_not_use_rairv_init( &
+ use_rairv, &
+ errmsg, errflg)
+
+ ! Output arguments
+ logical, intent(out) :: use_rairv ! Flag for constituent-dependent gas constant [flag]
+ character(len=512), intent(out) :: errmsg
+ integer, intent(out) :: errflg
+
+ errmsg = ''
+ errflg = 0
+
+ ! Set to .false. since we are not in WACCM-X mode
+ ! Only use constituent-dependent gas constant when in WACCM-X mode
+ use_rairv = .false.
+
+ end subroutine vertical_diffusion_not_use_rairv_init
+
+
+end module diffusion_stubs
diff --git a/schemes/vertical_diffusion/diffusion_stubs.meta b/schemes/vertical_diffusion/diffusion_stubs.meta
new file mode 100644
index 00000000..ff2836e0
--- /dev/null
+++ b/schemes/vertical_diffusion/diffusion_stubs.meta
@@ -0,0 +1,415 @@
+[ccpp-table-properties]
+ name = zero_upper_boundary_condition
+ type = scheme
+
+[ccpp-arg-table]
+ name = zero_upper_boundary_condition_init
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_dimension
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+[ pcnst ]
+ standard_name = number_of_ccpp_constituents
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+[ ubc_mmr ]
+ standard_name = upper_boundary_condition_of_ccpp_constituents
+ units = none
+ dimensions = (horizontal_dimension, number_of_ccpp_constituents)
+ type = real | kind = kind_phys
+ intent = out
+[ cnst_fixed_ubc ]
+ standard_name = use_fixed_upper_boundary_condition_of_ccpp_constituents_tbd
+ units = flag
+ dimensions = (number_of_ccpp_constituents)
+ type = logical
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ dimensions = ()
+ type = character | kind = len=512
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ dimensions = ()
+ type = integer
+ intent = out
+
+[ccpp-table-properties]
+ name = tms_beljaars_zero_stub
+ type = scheme
+
+[ccpp-arg-table]
+ name = tms_beljaars_zero_stub_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ ksrftms ]
+ standard_name = turbulent_orographic_form_drag_coefficent_at_surface
+ units = kg s-1 m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ tautmsx ]
+ standard_name = eastward_turbulent_mountain_surface_stress_tbd
+ units = N m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = out
+[ tautmsy ]
+ standard_name = northward_turbulent_mountain_surface_stress_tbd
+ units = N m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = out
+[ do_beljaars ]
+ standard_name = do_beljaars_drag_in_vertical_diffusion_tbd
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = out
+[ dragblj ]
+ standard_name = turbulent_orographic_form_drag_coefficent
+ units = s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = out
+[ taubljx ]
+ standard_name = eastward_beljaars_surface_stress_tbd
+ units = N m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = out
+[ taubljy ]
+ standard_name = northward_beljaars_surface_stress_tbd
+ units = N m-2
+ dimensions = (horizontal_loop_extent)
+ type = real | kind = kind_phys
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ dimensions = ()
+ type = character | kind = len=512
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ dimensions = ()
+ type = integer
+ intent = out
+
+[ccpp-table-properties]
+ name = turbulent_mountain_stress_add_drag_coefficient
+ type = scheme
+
+[ccpp-arg-table]
+ name = turbulent_mountain_stress_add_drag_coefficient_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ ksrftms ]
+ standard_name = turbulent_orographic_form_drag_coefficent_at_surface
+ units = kg s-1 m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ ksrf ]
+ standard_name = total_drag_coefficient_at_surface
+ units = kg m-2 s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = inout
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ dimensions = ()
+ type = character | kind = len=512
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ dimensions = ()
+ type = integer
+ intent = out
+
+[ccpp-table-properties]
+ name = beljaars_add_wind_damping_rate
+ type = scheme
+
+[ccpp-arg-table]
+ name = beljaars_add_wind_damping_rate_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ dragblj ]
+ standard_name = turbulent_orographic_form_drag_coefficent
+ units = s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ tau_damp_rate ]
+ standard_name = horizontal_wind_damping_rate_due_to_total_surface_drag
+ units = s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = inout
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-table-properties]
+ name = beljaars_add_updated_residual_stress
+ type = scheme
+
+[ccpp-arg-table]
+ name = beljaars_add_updated_residual_stress_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ gravit ]
+ standard_name = standard_gravitational_acceleration
+ units = m s-2
+ type = real | kind = kind_phys
+ dimensions = ()
+ intent = in
+[ p ]
+ standard_name = vertical_moist_pressure_coordinates_for_vertical_diffusion
+ units = none
+ type = coords1d
+ dimensions = ()
+ intent = in
+[ do_iss ]
+ standard_name = do_implicit_total_surface_stress_in_vertical_diffusion
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = in
+[ itaures ]
+ standard_name = update_residual_stress_at_surface_in_vertical_diffusion
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = in
+[ dragblj ]
+ standard_name = turbulent_orographic_form_drag_coefficent
+ units = s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ u1 ]
+ standard_name = eastward_wind_after_vertical_diffusion
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ v1 ]
+ standard_name = northward_wind_after_vertical_diffusion
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ tauresx ]
+ standard_name = eastward_reserved_stress_at_surface_on_previous_timestep
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = inout
+[ tauresy ]
+ standard_name = northward_reserved_stress_at_surface_on_previous_timestep
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = inout
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-table-properties]
+ name = turbulent_mountain_stress_add_updated_surface_stress
+ type = scheme
+
+[ccpp-arg-table]
+ name = turbulent_mountain_stress_add_updated_surface_stress_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ pver ]
+ standard_name = vertical_layer_dimension
+ units = count
+ type = integer
+ dimensions = ()
+ intent = in
+[ do_iss ]
+ standard_name = do_implicit_total_surface_stress_in_vertical_diffusion
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = in
+[ itaures ]
+ standard_name = update_residual_stress_at_surface_in_vertical_diffusion
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = in
+[ ksrftms ]
+ standard_name = turbulent_orographic_form_drag_coefficent_at_surface
+ units = kg s-1 m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = in
+[ u1 ]
+ standard_name = eastward_wind_after_vertical_diffusion
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ v1 ]
+ standard_name = northward_wind_after_vertical_diffusion
+ units = m s-1
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+ intent = in
+[ tauresx ]
+ standard_name = eastward_reserved_stress_at_surface_on_previous_timestep
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = inout
+[ tauresy ]
+ standard_name = northward_reserved_stress_at_surface_on_previous_timestep
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = inout
+[ tautotx ]
+ standard_name = eastward_stress_at_surface_for_vertical_diffusion
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = inout
+[ tautoty ]
+ standard_name = northward_stress_at_surface_for_vertical_diffusion
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = inout
+[ tautmsx ]
+ standard_name = eastward_turbulent_mountain_surface_stress_tbd
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ tautmsy ]
+ standard_name = northward_turbulent_mountain_surface_stress_tbd
+ units = N m-2
+ type = real | kind = kind_phys
+ dimensions = (horizontal_loop_extent)
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
+
+[ccpp-table-properties]
+ name = vertical_diffusion_not_use_rairv
+ type = scheme
+
+[ccpp-arg-table]
+ name = vertical_diffusion_not_use_rairv_init
+ type = scheme
+[ use_rairv ]
+ standard_name = use_composition_dependent_gas_constant_in_vertical_diffusion
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = out
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
diff --git a/schemes/vertical_diffusion/vertical_diffusion_diagnostic_utils.F90 b/schemes/vertical_diffusion/vertical_diffusion_diagnostic_utils.F90
new file mode 100644
index 00000000..c38d4557
--- /dev/null
+++ b/schemes/vertical_diffusion/vertical_diffusion_diagnostic_utils.F90
@@ -0,0 +1,284 @@
+! Diagnostic helper functions for vertical diffusion
+module vertical_diffusion_diagnostic_utils
+ use ccpp_kinds, only: kind_phys
+
+ implicit none
+ private
+
+ public :: vertical_diffusion_diagnostic_profiles
+
+contains
+
+ ! Computes profile output before and after vertical diffusion
+ ! for PBL analysis
+ !
+ ! This subroutine was designed to be reusable in current CAM and SIMA
+ subroutine vertical_diffusion_diagnostic_profiles( &
+ ncol, pver, pverp, &
+ dt, &
+ latvap, latice, zvir, cpair, gravit, rair, &
+ pint, pmid, zi, zm, &
+ kvh, kvm, cgh, cgs, &
+ cam_in_shf, &
+ cam_in_cflx_wv, cam_in_cflx_cldliq, cam_in_cflx_cldice, &
+ tautotx, tautoty, &
+ t0, q0_wv, q0_cldliq, q0_cldice, s0, u0, v0, &
+ s1, u1, v1, &
+ tend_s, tend_u, tend_v, tend_q_wv, tend_q_cldliq, tend_q_cldice, &
+ qt_pre_PBL, sl_pre_PBL, slv_pre_PBL, rh_pre_PBL, &
+ qt_aft_PBL, sl_aft_PBL, slv_aft_PBL, &
+ qv_aft_PBL, ql_aft_PBL, qi_aft_PBL, &
+ t_aft_PBL, rh_aft_PBL, &
+ u_aft_PBL, v_aft_PBL, &
+ slten, qtten, tten, rhten, &
+ slflx, qtflx, uflx, vflx, &
+ slflx_cg, qtflx_cg, uflx_cg, vflx_cg)
+
+ ! To-be-ccppized dependency.
+ use wv_saturation, only: qsat
+
+ ! Input arguments
+ integer, intent(in) :: ncol
+ integer, intent(in) :: pver
+ integer, intent(in) :: pverp
+ real(kind_phys), intent(in) :: dt ! Physics timestep [s]
+ real(kind_phys), intent(in) :: latvap ! Latent heat of vaporization [J kg-1]
+ real(kind_phys), intent(in) :: latice ! Latent heat of fusion [J kg-1]
+ real(kind_phys), intent(in) :: zvir ! rh2o/rair - 1 [1]
+ real(kind_phys), intent(in) :: cpair ! Specific heat of dry air at constant pressure [J kg-1 K-1]
+ real(kind_phys), intent(in) :: gravit ! Acceleration due to gravity [m s-2]
+ real(kind_phys), intent(in) :: rair ! Gas constant for dry air [J kg-1 K-1]
+ real(kind_phys), intent(in) :: pint(:, :) ! Interface pressure [Pa]
+ real(kind_phys), intent(in) :: pmid(:, :) ! Midpoint pressure [Pa]
+ real(kind_phys), intent(in) :: zi(:, :) ! Geopotential height at interfaces [m]
+ real(kind_phys), intent(in) :: zm(:, :) ! Geopotential height at midpoints [m]
+
+ ! Diffusion coefficients and counter-gradient terms
+ real(kind_phys), intent(in) :: kvh(:, :) ! Eddy diffusivity for heat at interfaces [m2 s-1]
+ real(kind_phys), intent(in) :: kvm(:, :) ! Eddy diffusivity for momentum at interfaces [m2 s-1]
+ real(kind_phys), intent(in) :: cgh(:, :) ! Counter-gradient term for heat [K m s-1]
+ real(kind_phys), intent(in) :: cgs(:, :) ! Counter-gradient star [s m-2]
+
+ ! Surface fluxes from coupler
+ real(kind_phys), intent(in) :: cam_in_shf(:) ! Surface sensible heat flux [W m-2]
+
+ ! NOTE: These surface fluxes need to be disassembled from cam_in_cflx into wv, cldliq, cldice
+ real(kind_phys), intent(in) :: cam_in_cflx_wv(:) ! Surface water vapor flux [kg m-2 s-1]
+ real(kind_phys), intent(in) :: cam_in_cflx_cldliq(:) ! Surface cloud liquid flux [kg m-2 s-1]
+ real(kind_phys), intent(in) :: cam_in_cflx_cldice(:) ! Surface cloud ice flux [kg m-2 s-1]
+
+ ! Total surface stresses
+ ! NOTE: pre-update (not using provisional winds for TMS/Beljaars orographic surface drag)
+ real(kind_phys), intent(in) :: tautotx(:) ! Total zonal surface stress [N m-2]
+ real(kind_phys), intent(in) :: tautoty(:) ! Total meridional surface stress [N m-2]
+
+ ! Initial state (before vertical diffusion)
+ real(kind_phys), intent(in) :: t0(:, :) ! Temperature before diffusion [K]
+ real(kind_phys), intent(in) :: q0_wv(:, :) ! Water vapor mixing ratio before diffusion [kg kg-1]
+ real(kind_phys), intent(in) :: q0_cldliq(:, :) ! Cloud liquid water mixing ratio before diffusion [kg kg-1]
+ real(kind_phys), intent(in) :: q0_cldice(:, :) ! Cloud ice water mixing ratio before diffusion [kg kg-1]
+ real(kind_phys), intent(in) :: s0(:, :) ! Dry static energy before diffusion [J kg-1]
+ real(kind_phys), intent(in) :: u0(:, :) ! Zonal wind before diffusion [m s-1]
+ real(kind_phys), intent(in) :: v0(:, :) ! Meridional wind before diffusion [m s-1]
+
+ ! Final state (after vertical diffusion)
+ real(kind_phys), intent(in) :: s1(:, :) ! Dry static energy after diffusion [J kg-1]
+ real(kind_phys), intent(in) :: u1(:, :) ! Zonal wind after diffusion [m s-1]
+ real(kind_phys), intent(in) :: v1(:, :) ! Meridional wind after diffusion [m s-1]
+
+ ! Tendencies attributable to vertical diffusion
+ real(kind_phys), intent(in) :: tend_s(:, :) ! Dry static energy tendency [J kg-1 s-1]
+ real(kind_phys), intent(in) :: tend_u(:, :) ! Zonal wind tendency [m s-2]
+ real(kind_phys), intent(in) :: tend_v(:, :) ! Meridional wind tendency [m s-2]
+
+ ! Constituent tendencies from vertical diffusion
+ real(kind_phys), intent(in) :: tend_q_wv(:, :) ! Water vapor tendency [kg kg-1 s-1]
+ real(kind_phys), intent(in) :: tend_q_cldliq(:, :) ! Cloud liquid water tendency [kg kg-1 s-1]
+ real(kind_phys), intent(in) :: tend_q_cldice(:, :) ! Cloud ice water tendency [kg kg-1 s-1]
+
+ ! Output diagnostic profiles - before vertical diffusion
+ real(kind_phys), intent(out) :: qt_pre_PBL(:, :) ! Total water mixing ratio before PBL [kg kg-1]
+ real(kind_phys), intent(out) :: sl_pre_PBL(:, :) ! Liquid water static energy before PBL [J kg-1]
+ real(kind_phys), intent(out) :: slv_pre_PBL(:, :) ! Virtual liquid water static energy before PBL [J kg-1]
+ real(kind_phys), intent(out) :: rh_pre_PBL(:, :) ! Relative humidity before PBL [percent]
+
+ ! Output diagnostic profiles - after vertical diffusion
+ real(kind_phys), intent(out) :: qt_aft_PBL(:, :) ! Total water mixing ratio after PBL [kg kg-1]
+ real(kind_phys), intent(out) :: sl_aft_PBL(:, :) ! Liquid water static energy after PBL [J kg-1]
+ real(kind_phys), intent(out) :: slv_aft_PBL(:, :) ! Virtual liquid water static energy after PBL [J kg-1]
+ real(kind_phys), intent(out) :: qv_aft_PBL(:, :) ! Water vapor mixing ratio after PBL [kg kg-1]
+ real(kind_phys), intent(out) :: ql_aft_PBL(:, :) ! Cloud liquid water mixing ratio after PBL [kg kg-1]
+ real(kind_phys), intent(out) :: qi_aft_PBL(:, :) ! Cloud ice water mixing ratio after PBL [kg kg-1]
+ real(kind_phys), intent(out) :: t_aft_PBL(:, :) ! Temperature after PBL [K]
+ real(kind_phys), intent(out) :: rh_aft_PBL(:, :) ! Relative humidity after PBL [percent]
+ real(kind_phys), intent(out) :: u_aft_PBL(:, :) ! Zonal wind after PBL [m s-1]
+ real(kind_phys), intent(out) :: v_aft_PBL(:, :) ! Meridional wind after PBL [m s-1]
+
+ ! Output tendency diagnostics
+ real(kind_phys), intent(out) :: slten(:, :) ! Liquid water static energy tendency [J kg-1 s-1]
+ real(kind_phys), intent(out) :: qtten(:, :) ! Total water mixing ratio tendency [kg kg-1 s-1]
+ real(kind_phys), intent(out) :: tten(:, :) ! Temperature tendency [K s-1]
+ real(kind_phys), intent(out) :: rhten(:, :) ! Relative humidity tendency [% s-1]
+
+ ! Output flux diagnostics
+ real(kind_phys), intent(out) :: slflx(:, :) ! Liquid static energy flux at interfaces [W m-2]
+ real(kind_phys), intent(out) :: qtflx(:, :) ! Total water flux at interfaces [kg m-2 s-1]
+ real(kind_phys), intent(out) :: uflx(:, :) ! Zonal momentum flux at interfaces [N m-2]
+ real(kind_phys), intent(out) :: vflx(:, :) ! Meridional momentum flux at interfaces [N m-2]
+ real(kind_phys), intent(out) :: slflx_cg(:, :) ! Counter-gradient liquid static energy flux at interfaces [W m-2]
+ real(kind_phys), intent(out) :: qtflx_cg(:, :) ! Counter-gradient total water flux at interfaces [kg m-2 s-1]
+ real(kind_phys), intent(out) :: uflx_cg(:, :) ! Counter-gradient zonal momentum flux at interfaces [N m-2]
+ real(kind_phys), intent(out) :: vflx_cg(:, :) ! Counter-gradient meridional momentum flux at interfaces [N m-2]
+
+ ! Local variables
+ integer :: i, k
+ real(kind_phys) :: rdt ! Reciprocal of timestep [s-1]
+ real(kind_phys) :: rhoair ! Air density [kg m-3]
+ real(kind_phys) :: sat_pre(ncol, pver) ! Saturation vapor pressure before PBL [Pa]
+ real(kind_phys) :: sat_aft(ncol, pver) ! Saturation vapor pressure after PBL [Pa]
+ real(kind_phys) :: tem2_pre(ncol, pver) ! Saturation specific humidity before PBL [kg kg-1]
+ real(kind_phys) :: tem2_aft(ncol, pver) ! Saturation specific humidity after PBL [kg kg-1]
+ real(kind_phys) :: s_aft_PBL(ncol, pver) ! Dry static energy after PBL [J kg-1]
+
+ rdt = 1.0_kind_phys / dt
+
+ ! ====================================================
+ ! Compute profiles before vertical diffusion (pre-PBL)
+ ! ====================================================
+
+ ! Liquid water static energy before PBL
+ sl_pre_PBL(:ncol, :pver) = s0(:ncol, :pver) &
+ - latvap * q0_cldliq(:ncol, :pver) &
+ - (latvap + latice) * q0_cldice(:ncol, :pver)
+
+ ! Total water mixing ratio before PBL
+ qt_pre_PBL(:ncol, :pver) = q0_wv(:ncol, :pver) &
+ + q0_cldliq(:ncol, :pver) &
+ + q0_cldice(:ncol, :pver)
+
+ ! Virtual liquid water static energy before PBL
+ slv_pre_PBL(:ncol, :pver) = sl_pre_PBL(:ncol, :pver) &
+ * (1.0_kind_phys + zvir * qt_pre_PBL(:ncol, :pver))
+
+ ! Compute saturation vapor pressure and relative humidity before PBL
+ do k = 1, pver
+ call qsat(t0(:ncol, k), pmid(:ncol, k), tem2_pre(:ncol, k), sat_pre(:ncol, k), ncol)
+ end do
+ rh_pre_PBL(:ncol, :pver) = q0_wv(:ncol, :pver) / sat_pre(:ncol, :pver) * 100.0_kind_phys
+
+ ! ====================================================
+ ! Compute profiles after vertical diffusion (post-PBL)
+ ! ====================================================
+
+ ! Apply tendencies to get final state variables
+ !
+ ! Note: this is indeed computed twice because at the point this utility subroutine runs,
+ ! the tendencies have not yet been applied.
+ qv_aft_PBL(:ncol, :pver) = q0_wv(:ncol, :pver) + tend_q_wv(:ncol, :pver) * dt
+ ql_aft_PBL(:ncol, :pver) = q0_cldliq(:ncol, :pver) + tend_q_cldliq(:ncol, :pver) * dt
+ qi_aft_PBL(:ncol, :pver) = q0_cldice(:ncol, :pver) + tend_q_cldice(:ncol, :pver) * dt
+ u_aft_PBL(:ncol, :pver) = u0(:ncol, :pver) + tend_u(:ncol, :pver) * dt
+ v_aft_PBL(:ncol, :pver) = v0(:ncol, :pver) + tend_v(:ncol, :pver) * dt
+
+ ! Dry static energy after PBL
+ s_aft_PBL(:ncol, :pver) = s0(:ncol, :pver) + tend_s(:ncol, :pver) * dt
+
+ ! Liquid water static energy after PBL
+ sl_aft_PBL(:ncol, :pver) = s1(:ncol, :pver) &
+ - latvap * ql_aft_PBL(:ncol, :pver) &
+ - (latvap + latice) * qi_aft_PBL(:ncol, :pver)
+
+ ! Total water mixing ratio after PBL
+ qt_aft_PBL(:ncol, :pver) = qv_aft_PBL(:ncol, :pver) &
+ + ql_aft_PBL(:ncol, :pver) &
+ + qi_aft_PBL(:ncol, :pver)
+
+ ! Virtual liquid water static energy after PBL
+ slv_aft_PBL(:ncol, :pver) = sl_aft_PBL(:ncol, :pver) &
+ * (1.0_kind_phys + zvir * qt_aft_PBL(:ncol, :pver))
+
+ ! Temperature after PBL (derived from dry static energy)
+ t_aft_PBL(:ncol, :pver) = (s_aft_PBL(:ncol, :pver) - gravit * zm(:ncol, :pver)) / cpair
+
+ ! Compute saturation vapor pressure and relative humidity after PBL
+ do k = 1, pver
+ call qsat(t_aft_PBL(:ncol, k), pmid(:ncol, k), tem2_aft(:ncol, k), sat_aft(:ncol, k), ncol)
+ end do
+ rh_aft_PBL(:ncol, :pver) = qv_aft_PBL(:ncol, :pver) / sat_aft(:ncol, :pver) * 100.0_kind_phys
+
+ ! ====================================================
+ ! Compute tendency diagnostics
+ ! ====================================================
+
+ ! Liquid water static energy tendency
+ slten(:ncol, :pver) = (sl_aft_PBL(:ncol, :pver) - sl_pre_PBL(:ncol, :pver)) * rdt
+
+ ! Total water mixing ratio tendency
+ qtten(:ncol, :pver) = (qt_aft_PBL(:ncol, :pver) - qt_pre_PBL(:ncol, :pver)) * rdt
+
+ ! Temperature tendency
+ tten(:ncol, :pver) = (t_aft_PBL(:ncol, :pver) - t0(:ncol, :pver)) * rdt
+
+ ! Relative humidity tendency
+ rhten(:ncol, :pver) = (rh_aft_PBL(:ncol, :pver) - rh_pre_PBL(:ncol, :pver)) * rdt
+
+ ! ===================================================================
+ ! Compute flux diagnostics at interfaces
+ ! ===================================================================
+
+ ! Initialize top interface (no flux at TOA)
+ slflx(:ncol, 1) = 0.0_kind_phys
+ qtflx(:ncol, 1) = 0.0_kind_phys
+ uflx(:ncol, 1) = 0.0_kind_phys
+ vflx(:ncol, 1) = 0.0_kind_phys
+ slflx_cg(:ncol, 1) = 0.0_kind_phys
+ qtflx_cg(:ncol, 1) = 0.0_kind_phys
+ uflx_cg(:ncol, 1) = 0.0_kind_phys
+ vflx_cg(:ncol, 1) = 0.0_kind_phys
+
+ ! Compute fluxes at model interfaces
+ do k = 2, pver
+ do i = 1, ncol
+ ! Air density at interface k [kg m-3]
+ rhoair = pint(i, k) / (rair * ((0.5_kind_phys * (slv_aft_PBL(i, k) + slv_aft_PBL(i, k-1)) &
+ - gravit * zi(i, k)) / cpair))
+
+ ! Liquid static energy flux [W m-2]
+ slflx(i, k) = kvh(i, k) * &
+ (cgh(i, k) - rhoair * (sl_aft_PBL(i, k-1) - sl_aft_PBL(i, k)) / (zm(i, k-1) - zm(i, k)))
+
+ ! Total water flux [kg m-2 s-1]
+ qtflx(i, k) = kvh(i, k) * &
+ ( - rhoair * (qt_aft_PBL(i, k-1) - qt_aft_PBL(i, k)) / (zm(i, k-1) - zm(i, k)) &
+ + rhoair * (cam_in_cflx_wv(i) + cam_in_cflx_cldliq(i) + cam_in_cflx_cldice(i)) * cgs(i, k) )
+
+ ! Zonal momentum flux [N m-2]
+ uflx(i, k) = kvm(i, k) * &
+ ( - rhoair * (u_aft_PBL(i, k-1) - u_aft_PBL(i, k)) / (zm(i, k-1) - zm(i, k)) )
+
+ ! Meridional momentum flux [N m-2]
+ vflx(i, k) = kvm(i, k) * &
+ ( - rhoair * (v_aft_PBL(i, k-1) - v_aft_PBL(i, k)) / (zm(i, k-1) - zm(i, k)) )
+
+ ! Counter-gradient fluxes
+ slflx_cg(i, k) = kvh(i, k) * cgh(i, k)
+ qtflx_cg(i, k) = kvh(i, k) * rhoair * (cam_in_cflx_wv(i) + cam_in_cflx_cldliq(i) &
+ + cam_in_cflx_cldice(i)) * cgs(i, k)
+ uflx_cg(i, k) = 0.0_kind_phys
+ vflx_cg(i, k) = 0.0_kind_phys
+ end do
+ end do
+
+ ! Set fluxes at bottom interface using surface fluxes
+ slflx(:ncol, pverp) = cam_in_shf(:ncol)
+ qtflx(:ncol, pverp) = cam_in_cflx_wv(:ncol)
+ uflx(:ncol, pverp) = tautotx(:ncol)
+ vflx(:ncol, pverp) = tautoty(:ncol)
+ slflx_cg(:ncol, pverp) = 0.0_kind_phys
+ qtflx_cg(:ncol, pverp) = 0.0_kind_phys
+ uflx_cg(:ncol, pverp) = 0.0_kind_phys
+ vflx_cg(:ncol, pverp) = 0.0_kind_phys
+
+ end subroutine vertical_diffusion_diagnostic_profiles
+
+end module vertical_diffusion_diagnostic_utils
diff --git a/schemes/vertical_diffusion/vertical_diffusion_options.F90 b/schemes/vertical_diffusion/vertical_diffusion_options.F90
new file mode 100644
index 00000000..7648a88f
--- /dev/null
+++ b/schemes/vertical_diffusion/vertical_diffusion_options.F90
@@ -0,0 +1,42 @@
+! Module to load namelist options for vertical diffusion solver scheme.
+module vertical_diffusion_options
+
+ use ccpp_kinds, only: kind_phys
+
+ implicit none
+ private
+ save
+
+ ! CCPP-compliant public interfaces
+ public :: vertical_diffusion_options_init
+
+contains
+
+!> \section arg_table_vertical_diffusion_options_init Argument Table
+!! \htmlinclude arg_table_vertical_diffusion_options_init.html
+ subroutine vertical_diffusion_options_init( &
+ amIRoot, iulog, &
+ do_iss_in, &
+ am_correction_in, &
+ errmsg, errflg)
+
+ ! Input arguments
+ logical, intent(in) :: amIRoot ! are we on the MPI root task?
+ integer, intent(in) :: iulog ! log output unit
+ logical, intent(in) :: do_iss_in ! Use implicit turbulent surface stress computation
+ logical, intent(in) :: am_correction_in ! Do angular momentum conservation correction
+
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ errmsg = ''
+ errflg = 0
+
+ if(amIRoot) then
+ write(iulog,*) "vertical diffusion solver: do_iss ", do_iss_in
+ write(iulog,*) "vertical diffusion solver: am_correction ", am_correction_in
+ endif
+
+ end subroutine vertical_diffusion_options_init
+
+end module vertical_diffusion_options
diff --git a/schemes/vertical_diffusion/vertical_diffusion_options.meta b/schemes/vertical_diffusion/vertical_diffusion_options.meta
new file mode 100644
index 00000000..80197e0e
--- /dev/null
+++ b/schemes/vertical_diffusion/vertical_diffusion_options.meta
@@ -0,0 +1,43 @@
+[ccpp-table-properties]
+ name = vertical_diffusion_options
+ type = scheme
+
+[ccpp-arg-table]
+ name = vertical_diffusion_options_init
+ type = scheme
+[ amIRoot ]
+ standard_name = flag_for_mpi_root
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = in
+[ iulog ]
+ standard_name = log_output_unit
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = in
+[ do_iss_in ]
+ standard_name = do_implicit_total_surface_stress_in_vertical_diffusion
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = in
+[ am_correction_in ]
+ standard_name = do_angular_momentum_correction_in_vertical_diffusion
+ units = flag
+ type = logical
+ dimensions = ()
+ intent = in
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ type = character | kind = len=512
+ dimensions = ()
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ type = integer
+ dimensions = ()
+ intent = out
diff --git a/schemes/vertical_diffusion/vertical_diffusion_options_namelist.xml b/schemes/vertical_diffusion/vertical_diffusion_options_namelist.xml
new file mode 100644
index 00000000..cb3176f9
--- /dev/null
+++ b/schemes/vertical_diffusion/vertical_diffusion_options_namelist.xml
@@ -0,0 +1,108 @@
+
+
+
+
+
+
+
+
+ logical
+ pbl
+ vertical_diffusion_nl
+ do_implicit_total_surface_stress_in_vertical_diffusion
+ flag
+
+ Logical switch to turn on implicit turbulent surface stress calculation in
+ diffusion solver routine. Default: FALSE.
+
+
+ false
+ true
+ true
+
+
+
+ logical
+ pbl
+ vertical_diffusion_nl
+ do_angular_momentum_correction_in_vertical_diffusion
+ flag
+
+ Flag to turn on corrections that improve angular momentum conservation.
+ Default: FALSE.
+
+
+ false
+
+
+
diff --git a/schemes/vertical_diffusion/vertical_diffusion_sponge_layer.F90 b/schemes/vertical_diffusion/vertical_diffusion_sponge_layer.F90
new file mode 100644
index 00000000..50a0da23
--- /dev/null
+++ b/schemes/vertical_diffusion/vertical_diffusion_sponge_layer.F90
@@ -0,0 +1,138 @@
+! Vertical diffusion sponge layer scheme
+! Adds artificial sponge layer vertical diffusion at the top of the atmosphere
+!
+! CCPP-ized: Haipeng Lin, June 2025
+module vertical_diffusion_sponge_layer
+ use ccpp_kinds, only: kind_phys
+
+ implicit none
+ private
+ save
+
+ ! CCPP-compliant public interfaces
+ public :: vertical_diffusion_sponge_layer_init
+ public :: vertical_diffusion_sponge_layer_run
+ public :: vertical_diffusion_sponge_layer_final
+
+ ! Module variables for sponge layer parameters
+ real(kind_phys), allocatable :: kvm_sponge(:) ! sponge layer diffusion coefficients [m^2 s-1]
+
+contains
+
+!> \section arg_table_vertical_diffusion_sponge_layer_init Argument Table
+!! \htmlinclude vertical_diffusion_sponge_layer_init.html
+ subroutine vertical_diffusion_sponge_layer_init( &
+ amIRoot, iulog, &
+ ptop_ref, &
+ errmsg, errflg)
+
+ ! Input arguments
+ logical, intent(in) :: amIRoot ! are we on the MPI root task?
+ integer, intent(in) :: iulog ! log output unit
+ real(kind_phys), intent(in) :: ptop_ref ! reference top pressure [Pa]
+
+ ! Output arguments
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ ! Local variables
+ integer :: k
+
+ errmsg = ''
+ errflg = 0
+
+ ! Add sponge layer vertical diffusion based on model top pressure
+ if (ptop_ref > 1.e-1_kind_phys .and. ptop_ref < 100.0_kind_phys) then
+ !
+ ! CAM7 FMT (but not CAM6 top (~225 Pa) or CAM7 low top or lower)
+ !
+ allocate(kvm_sponge(4), stat=errflg, errmsg=errmsg)
+ if (errflg /= 0) then
+ return
+ end if
+ kvm_sponge(1) = 2.e6_kind_phys
+ kvm_sponge(2) = 2.e6_kind_phys
+ kvm_sponge(3) = 0.5e6_kind_phys
+ kvm_sponge(4) = 0.1e6_kind_phys
+ else if (ptop_ref > 1.e-4_kind_phys) then
+ !
+ ! WACCM and WACCM-X
+ !
+ allocate(kvm_sponge(6), stat=errflg, errmsg=errmsg)
+ if (errflg /= 0) then
+ return
+ end if
+ kvm_sponge(1) = 2.e6_kind_phys
+ kvm_sponge(2) = 2.e6_kind_phys
+ kvm_sponge(3) = 1.5e6_kind_phys
+ kvm_sponge(4) = 1.0e6_kind_phys
+ kvm_sponge(5) = 0.5e6_kind_phys
+ kvm_sponge(6) = 0.1e6_kind_phys
+ end if
+
+ if (amIRoot) then
+ if (allocated(kvm_sponge)) then
+ write(iulog, *) 'Artificial sponge layer vertical diffusion added:'
+ do k = 1, size(kvm_sponge(:))
+ write(iulog, '(a44,i2,a17,e7.2,a8)') 'vertical diffusion coefficient at interface ', k, &
+ ' is increased by ', kvm_sponge(k), ' m2 s-1'
+ end do
+ else
+ write(iulog, *) 'No sponge layer vertical diffusion applied (ptop_ref = ', ptop_ref, ' Pa)'
+ end if
+ end if
+
+ end subroutine vertical_diffusion_sponge_layer_init
+
+!> \section arg_table_vertical_diffusion_sponge_layer_run Argument Table
+!! \htmlinclude vertical_diffusion_sponge_layer_run.html
+ subroutine vertical_diffusion_sponge_layer_run( &
+ ncol, pverp, &
+ kvm, &
+ errmsg, errflg)
+
+ ! Input arguments
+ integer, intent(in) :: ncol
+ integer, intent(in) :: pverp
+
+ ! Input/output arguments
+ real(kind_phys), intent(inout) :: kvm(:,:) ! Eddy diffusivity for momentum [m^2 s-1], interfaces
+
+ ! Output arguments
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ ! Local variables
+ integer :: k
+
+ errmsg = ''
+ errflg = 0
+
+ ! Add sponge layer vertical diffusion
+ if (allocated(kvm_sponge)) then
+ do k = 1, size(kvm_sponge(:))
+ kvm(:ncol, 1) = kvm(:ncol, 1) + kvm_sponge(k)
+ end do
+ end if
+
+ end subroutine vertical_diffusion_sponge_layer_run
+
+!> \section arg_table_vertical_diffusion_sponge_layer_final Argument Table
+!! \htmlinclude vertical_diffusion_sponge_layer_final.html
+ subroutine vertical_diffusion_sponge_layer_final( &
+ errmsg, errflg)
+
+ ! Output arguments
+ character(len=512), intent(out) :: errmsg ! error message
+ integer, intent(out) :: errflg ! error flag
+
+ errmsg = ''
+ errflg = 0
+
+ if (allocated(kvm_sponge)) then
+ deallocate(kvm_sponge)
+ end if
+
+ end subroutine vertical_diffusion_sponge_layer_final
+
+end module vertical_diffusion_sponge_layer
diff --git a/schemes/vertical_diffusion/vertical_diffusion_sponge_layer.meta b/schemes/vertical_diffusion/vertical_diffusion_sponge_layer.meta
new file mode 100644
index 00000000..d091f95c
--- /dev/null
+++ b/schemes/vertical_diffusion/vertical_diffusion_sponge_layer.meta
@@ -0,0 +1,87 @@
+[ccpp-table-properties]
+ name = vertical_diffusion_sponge_layer
+ type = scheme
+
+[ccpp-arg-table]
+ name = vertical_diffusion_sponge_layer_init
+ type = scheme
+[ amIRoot ]
+ standard_name = flag_for_mpi_root
+ units = flag
+ dimensions = ()
+ type = logical
+ intent = in
+[ iulog ]
+ standard_name = log_output_unit
+ units = 1
+ dimensions = ()
+ type = integer
+ intent = in
+[ ptop_ref ]
+ standard_name = air_pressure_at_top_of_atmosphere_model
+ units = Pa
+ dimensions = ()
+ type = real | kind = kind_phys
+ intent = in
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ dimensions = ()
+ type = character | kind = len=512
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ dimensions = ()
+ type = integer
+ intent = out
+
+[ccpp-arg-table]
+ name = vertical_diffusion_sponge_layer_run
+ type = scheme
+[ ncol ]
+ standard_name = horizontal_loop_extent
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+[ pverp ]
+ standard_name = vertical_interface_dimension
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+[ kvm ]
+ standard_name = eddy_momentum_diffusivity_at_interfaces
+ units = m2 s-1
+ dimensions = (horizontal_loop_extent, vertical_interface_dimension)
+ type = real | kind = kind_phys
+ intent = inout
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ dimensions = ()
+ type = character | kind = len=512
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ dimensions = ()
+ type = integer
+ intent = out
+
+[ccpp-arg-table]
+ name = vertical_diffusion_sponge_layer_final
+ type = scheme
+[ errmsg ]
+ standard_name = ccpp_error_message
+ units = none
+ dimensions = ()
+ type = character | kind = len=512
+ intent = out
+[ errflg ]
+ standard_name = ccpp_error_code
+ units = 1
+ dimensions = ()
+ type = integer
+ intent = out
diff --git a/suites/suite_cam4.xml b/suites/suite_cam4.xml
index 7c3fc7f8..a92bf774 100644
--- a/suites/suite_cam4.xml
+++ b/suites/suite_cam4.xml
@@ -7,9 +7,10 @@
Shallow convection Hack
Macrophysics RK
Microphysics RK
- Radiation RRTMGP (not implemented)
- Chemistry None (not implemented)
- Vertical Diffusion HB (not implemented)
+ Radiation RRTMGP (not implemented)
+ Chemistry None (not implemented)
+ Vertical Diffusion HB
+ Gravity Wave Drag Orographic (not implemented)
-->
@@ -227,10 +228,96 @@
+
+
+
+ holtslag_boville_diff_options
+ vertical_diffusion_options
+
+
+ zero_upper_boundary_condition
+
+
+ hb_diff_set_vertical_diffusion_top
+
+
+ tms_beljaars_zero_stub
+ hb_diff_set_total_surface_stress
+
+
+ hb_diff_prepare_vertical_diffusion_inputs
+
+
+ holtslag_boville_diff
+
+
+ hb_pbl_independent_coefficients
+ hb_pbl_dependent_coefficients
+ hb_diff_exchange_coefficients
+
+
+ vertical_diffusion_sponge_layer
+
+
+ holtslag_boville_diff_diagnostics
+
+
+
+ vertical_diffusion_not_use_rairv
+ vertical_diffusion_set_temperature_at_toa_default
+ vertical_diffusion_interpolate_to_interfaces
+
+
+ implicit_surface_stress_add_drag_coefficient
+
+
+
+ vertical_diffusion_wind_damping_rate
+
+
+
+
+ vertical_diffusion_diffuse_horizontal_momentum
+
+
+ vertical_diffusion_set_dry_static_energy_at_toa_zero
+
+
+ vertical_diffusion_diffuse_dry_static_energy
+
+
+ vertical_diffusion_diffuse_tracers
+
+
+ vertical_diffusion_tendencies
+
+
+ vertical_diffusion_tendencies_diagnostics
+
+
+ apply_tendency_of_northward_wind
+ apply_tendency_of_eastward_wind
+ apply_constituent_tendencies
+ apply_heating_rate
+ qneg
+ geopotential_temp
+ update_dry_static_energy
+
+
+
check_energy_save_teout
diff --git a/test/test_suites/suite_vdiff_holtslag_boville.xml b/test/test_suites/suite_vdiff_holtslag_boville.xml
new file mode 100644
index 00000000..31b581c6
--- /dev/null
+++ b/test/test_suites/suite_vdiff_holtslag_boville.xml
@@ -0,0 +1,90 @@
+
+
+
+
+ initialize_constituents
+ to_be_ccppized_temporary
+
+
+ holtslag_boville_diff_options
+ vertical_diffusion_options
+
+
+ zero_upper_boundary_condition
+
+
+ hb_diff_set_vertical_diffusion_top
+
+
+ tms_beljaars_zero_stub
+ hb_diff_set_total_surface_stress
+
+
+ hb_diff_prepare_vertical_diffusion_inputs
+
+
+ holtslag_boville_diff
+
+
+ hb_pbl_independent_coefficients
+ hb_pbl_dependent_coefficients
+ hb_diff_exchange_coefficients
+
+
+ vertical_diffusion_sponge_layer
+
+
+ holtslag_boville_diff_diagnostics
+
+
+
+ vertical_diffusion_not_use_rairv
+ vertical_diffusion_set_temperature_at_toa_default
+ vertical_diffusion_interpolate_to_interfaces
+
+
+ implicit_surface_stress_add_drag_coefficient
+
+
+
+ vertical_diffusion_wind_damping_rate
+
+
+
+
+ vertical_diffusion_diffuse_horizontal_momentum
+
+
+ vertical_diffusion_set_dry_static_energy_at_toa_zero
+
+
+ vertical_diffusion_diffuse_dry_static_energy
+
+
+ vertical_diffusion_diffuse_tracers
+
+
+ vertical_diffusion_tendencies
+
+
+ vertical_diffusion_tendencies_diagnostics
+
+
+ apply_tendency_of_northward_wind
+ apply_tendency_of_eastward_wind
+ apply_constituent_tendencies
+ apply_heating_rate
+ qneg
+ geopotential_temp
+ update_dry_static_energy
+
+
diff --git a/to_be_ccppized/coords_1d.F90 b/to_be_ccppized/coords_1d.F90
index 6ed75d9f..40a8c0de 100644
--- a/to_be_ccppized/coords_1d.F90
+++ b/to_be_ccppized/coords_1d.F90
@@ -10,9 +10,11 @@ module coords_1d
private
save
- public :: Coords1D
-
- type :: Coords1D
+ public :: coords1d
+
+!! \section arg_table_coords1d
+!! \htmlinclude coords1d.html
+ type :: coords1d
! Number of sets of coordinates in the object.
integer :: n = 0
! Number of coordinates in each set.
@@ -36,7 +38,7 @@ module coords_1d
contains
procedure :: section
procedure :: finalize
- end type Coords1D
+ end type coords1d
interface Coords1D
module procedure new_Coords1D_from_fields
@@ -149,4 +151,4 @@ end subroutine guarded_deallocate
end subroutine finalize
end module coords_1d
-
\ No newline at end of file
+
diff --git a/to_be_ccppized/coords_1d.meta b/to_be_ccppized/coords_1d.meta
new file mode 100644
index 00000000..74371e7c
--- /dev/null
+++ b/to_be_ccppized/coords_1d.meta
@@ -0,0 +1,7 @@
+[ccpp-table-properties]
+ name = coords1d
+ type = ddt
+
+[ccpp-arg-table]
+ name = coords1d
+ type = ddt
diff --git a/to_be_ccppized/vdiff_lu_solver.F90 b/to_be_ccppized/vdiff_lu_solver.F90
new file mode 100644
index 00000000..b4826060
--- /dev/null
+++ b/to_be_ccppized/vdiff_lu_solver.F90
@@ -0,0 +1,207 @@
+module vdiff_lu_solver
+
+! This module provides a function returning the matrix decomposition for
+! an implicit finite volume solver for vertical diffusion. It accepts
+! diffusion coefficients, time/grid spacing, and boundary condition
+! objects, and returns a TriDiagDecomp object that can be used to diffuse
+! an array for one time step with the "left_div" method.
+
+use coords_1d, only: Coords1D
+use linear_1d_operators, only: TriDiagOp, operator(+), TriDiagDecomp
+
+implicit none
+private
+save
+
+! Public interfaces
+public :: vd_lu_decomp
+public :: fin_vol_lu_decomp
+
+! 8-byte real.
+integer, parameter :: r8 = selected_real_kind(12)
+
+contains
+
+! ========================================================================!
+
+! Designed to solve the equation:
+! dq/dt = c1 q'' + c2 q' + c q
+
+function vd_lu_decomp(dt, dp, coef_q, coef_q_d, coef_q_d2, upper_bndry, &
+ lower_bndry) result(decomp)
+
+ use linear_1d_operators, only: &
+ identity_operator, &
+ diagonal_operator, &
+ first_derivative, &
+ second_derivative, &
+ BoundaryType
+
+ ! ---------------------- !
+ ! Input-Output Arguments !
+ ! ---------------------- !
+
+ ! Time step.
+ real(r8), intent(in) :: dt
+ ! Grid spacing (deltas).
+ real(r8), USE_CONTIGUOUS intent(in) :: dp(:,:)
+
+ ! Coefficients for q, q', and q''.
+ real(r8), USE_CONTIGUOUS intent(in), optional :: coef_q(:,:), &
+ coef_q_d(:,:), coef_q_d2(:,:)
+
+ ! Boundary conditions (optional, default to 0 flux through boundary).
+ class(BoundaryType), target, intent(in), optional :: &
+ upper_bndry, lower_bndry
+
+ ! Output decomposition.
+ type(TriDiagDecomp) :: decomp
+
+ ! --------------- !
+ ! Local Variables !
+ ! --------------- !
+
+ ! Operator objects.
+ type(TriDiagOp) :: add_term
+ type(TriDiagOp) :: net_operator
+
+ ! ----------------------- !
+ ! Main Computation Begins !
+ ! ----------------------- !
+
+ if (present(coef_q)) then
+ net_operator = diagonal_operator(1._r8 - dt*coef_q)
+ else
+ net_operator = identity_operator(size(dp, 1), size(dp, 2) + 1)
+ end if
+
+ if (present(coef_q_d)) then
+ add_term = first_derivative(dp, upper_bndry, lower_bndry)
+ call add_term%lmult_as_diag(-dt*coef_q_d)
+ call net_operator%add(add_term)
+ end if
+
+ if (present(coef_q_d2)) then
+ add_term = second_derivative(dp, upper_bndry, lower_bndry)
+ call add_term%lmult_as_diag(-dt*coef_q_d2)
+ call net_operator%add(add_term)
+ end if
+
+ decomp = TriDiagDecomp(net_operator)
+
+ call net_operator%finalize()
+ call add_term%finalize()
+
+end function vd_lu_decomp
+
+! ========================================================================!
+
+! Designed to solve the equation:
+!
+! w * dq/dt = d/dp (D q' - v q) + c q
+!
+! where q is a grid-cell average, and p is the vertical coordinate
+! (presumably pressure).
+!
+! In this function, coef_q_weight == w, coef_q_diff == D,
+! coef_q_adv == v, and coef_q == c. All these are optional; omitting a
+! coefficient is equivalent to setting the entire array to 0.
+!
+! coef_q_diff and coef_q_adv are defined at the level interfaces, while
+! coef_q and coef_q_weight are grid-cell averages.
+
+function fin_vol_lu_decomp(dt, p, coef_q, coef_q_diff, coef_q_adv, &
+ coef_q_weight, upper_bndry, lower_bndry, graft_decomp) result(decomp)
+
+ use linear_1d_operators, only: &
+ zero_operator, &
+ diagonal_operator, &
+ diffusion_operator, &
+ advection_operator, &
+ BoundaryType
+
+ ! ---------------------- !
+ ! Input-Output Arguments !
+ ! ---------------------- !
+
+ ! Time step.
+ real(r8), intent(in) :: dt
+ ! Grid spacings.
+ type(Coords1D), intent(in) :: p
+
+ ! Coefficients for diffusion and advection.
+ !
+ ! The sizes must be consistent among all the coefficients that are
+ ! actually present, i.e. coef_q_diff and coef_q_adv should be one level
+ ! bigger than coef_q and coef_q_weight, and have the same column number.
+ real(r8), USE_CONTIGUOUS intent(in), optional :: coef_q(:,:), &
+ coef_q_diff(:,:), coef_q_adv(:,:), coef_q_weight(:,:)
+
+ ! Boundary conditions (optional, default to 0 flux through boundary).
+ class(BoundaryType), target, intent(in), optional :: &
+ upper_bndry, lower_bndry
+
+ ! Decomposition to graft onto. If this is provided, you can pass in
+ ! smaller coefficients.
+ type(TriDiagDecomp), intent(in), optional :: graft_decomp
+
+ ! Output decomposition.
+ type(TriDiagDecomp) :: decomp
+
+ ! --------------- !
+ ! Local Variables !
+ ! --------------- !
+
+ ! Operator objects.
+ type(TriDiagOp) :: add_term
+ type(TriDiagOp) :: net_operator
+
+ ! ----------------------- !
+ ! Main Computation Begins !
+ ! ----------------------- !
+
+ ! A diffusion term is probably present, so start with that. Otherwise
+ ! start with an operator of all 0s.
+
+ if (present(coef_q_diff)) then
+ net_operator = diffusion_operator(p, coef_q_diff, &
+ upper_bndry, lower_bndry)
+ else
+ net_operator = zero_operator(p%n, p%d)
+ end if
+
+ ! Constant term (damping).
+ if (present(coef_q)) then
+ add_term = diagonal_operator(coef_q)
+ call net_operator%add(add_term)
+ end if
+
+ ! Effective advection.
+ if (present(coef_q_adv)) then
+ add_term = advection_operator(p, coef_q_adv, &
+ upper_bndry, lower_bndry)
+ call net_operator%add(add_term)
+ end if
+
+ ! We want I-dt*(w^-1)*A for a single time step, implicit method, where
+ ! A is the right-hand-side operator (i.e. what net_operator is now).
+ if (present(coef_q_weight)) then
+ call net_operator%lmult_as_diag(-dt/coef_q_weight)
+ else
+ call net_operator%lmult_as_diag(-dt)
+ end if
+ call net_operator%add_to_diag(1._r8)
+
+ ! Decompose, grafting on an optional input decomp. The graft is a way to
+ ! avoid re-calculating the ending (bottom) levels when the coefficients
+ ! have only changed at the beginning (top), e.g. for different
+ ! constituents in the molecular diffusion.
+ decomp = TriDiagDecomp(net_operator, graft_decomp=graft_decomp)
+
+ ! Ensure local objects are deallocated.
+ call net_operator%finalize()
+ call add_term%finalize()
+
+end function fin_vol_lu_decomp
+
+end module vdiff_lu_solver