From c5fde0a6f0118c3ce874f80ea3f2e11e935a1eed Mon Sep 17 00:00:00 2001 From: Jesse Nusbaumer Date: Wed, 24 Sep 2025 21:07:39 -0600 Subject: [PATCH 1/2] Add Github action to use git-fleximod for externals checking (#306) Originator(s): nusbaume Description (include issue title and the keyword ['closes', 'fixes', 'resolves'] and issue number): Adds a new Github Action that downloads and runs git-fleximod to ensure that the `.gitmodules` file is correct and matches the actual git submodules that are present in the repository. List all namelist files that were added or changed: N/A List all files eliminated and why: N/A List all files added and what they do: A .github/workflows/fleximod_test.yaml - Add new git-fleximod externals-checking Github Action workflow. List all existing files that have been modified, and describe the changes: (Helpful git command: `git diff --name-status development...`) List all automated tests that failed, as well as an explanation for why they weren't fixed: ALL PASS Is this an answer-changing PR? If so, is it a new physics package, algorithm change, tuning change, etc? No, no scientific source code was modified. If yes to the above question, describe how this code was validated with the new/modified features: --- .github/workflows/fleximod_test.yaml | 90 ++++++++++++++++++++++++++++ 1 file changed, 90 insertions(+) create mode 100644 .github/workflows/fleximod_test.yaml diff --git a/.github/workflows/fleximod_test.yaml b/.github/workflows/fleximod_test.yaml new file mode 100644 index 00000000..574984c7 --- /dev/null +++ b/.github/workflows/fleximod_test.yaml @@ -0,0 +1,90 @@ +# Github Action to ensure that +# the "fxtag" and "url" values +# in the ".gitmodules" file +# match the hashes or tags that +# are committed for the git submodules +# in this repository. + +# If this test is failing, verify that +# the commit hash or tag in the failing +# submodule folder is an exact match of +# the values specified for the corresponding +# submodule in the ".gitmodules" file. +# Both the submodule itself (referenced by Git), +# and the fxtag and url entires in the +# ".gitmodules" file, must be updated if +# a submodule is updated in this repository. + +# Finally, please note that 'git-fleximod' +# is not actually needed for this +# repo, it is simply being used +# here to run the test described +# above. + +on: pull_request + +jobs: + fleximod-test: + runs-on: ubuntu-latest + steps: + - name: checkout atmospheric_physics + uses: actions/checkout@v4 + + - name: checkout latest python version + uses: actions/setup-python@v4 + with: + python-version: "3.x" + + #Note: We are using the head of main for git-fleximod + # below for two reasons: + # 1. We will automatically get any git-fleximod bug fixes or updates + # 2. We will automatically test with those updates, which could + # allow us to catch any issues before updating the host models. + - name: checkout git-fleximod + uses: actions/checkout@v4 + with: + repository: ESMCI/git-fleximod + fetch-depth: 1 + ref: main + path: .lib/git-fleximod #Must use relative path here. + + - name: create git-fleximod script + run: | + echo + echo "Creating git-fleximod python script..." + + #Create new directory: + mkdir $GITHUB_WORKSPACE/bin + + #Write python script to expected location: + cat < $GITHUB_WORKSPACE/bin/git-fleximod + #!/usr/bin/env python3 + import sys + import os + sys.path.insert(0,os.path.abspath(os.path.join(os.path.dirname(__file__),"..",".lib","git-fleximod"))) + from git_fleximod.git_fleximod import main + + if __name__ == '__main__': + sys.exit(main()) + EOF + + #Make sure python script is executable: + chmod a+x $GITHUB_WORKSPACE/bin/git-fleximod + + echo "...git-fleximod python script generated." + + - name: run git-fleximod + run: | + $GITHUB_WORKSPACE/bin/git-fleximod update -o + echo + echo "Update complete, checking status" + echo + $GITHUB_WORKSPACE/bin/git-fleximod test + - name: check cleanliness + run: | + echo + echo "Checking if git fleximod matches expected externals" + echo + git add . + git diff --exit-code -- . ':(exclude).lib' ':(exclude)bin' + git diff --cached --exit-code -- . ':(exclude).lib' ':(exclude)bin' From 49716256e23e277d4095cee9efacbcbbeff0c112 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 26 Sep 2025 13:38:34 -0400 Subject: [PATCH 2/2] Complete CCPPization of Holtslag-Boville PBL scheme and vertical diffusion solver (#270) Companion PRs: https://github.com/ESCOMP/CAM/pull/1361 and https://github.com/ESCOMP/CAM-SIMA/pull/410 Originator(s): @jimmielin Description (include issue title and the keyword ['closes', 'fixes', 'resolves'] and issue number): * Partial cleanup and refactoring of HB PBL scheme for CCPPization. * Completes #239. Moves HB (`hb_diff.F90`) both `compute_hb_diff` and `compute_hb_free_atm_diff` to atmospheric_physics `holtslag_boville.F90`. * Interstitials necessary for inputs into HB in `holtslag_boville_interstitials.F90`. * (is there an issue for diffusion solver?) Moves vertical diffusion solver (`compute_vdiff`) without molecular diffusion (i.e., not WACCM-X) into CCPPized schemes in `diffusion_solver.F90`. * Creates stub (zeroed outputs) interstitial schemes for Turbulent Mountain Stress (TMS) and Beljaars orographic surface drag and other code for WACCM (e.g., UBC; use of composition-dependent rair) in `diffusion_stubs.F90` -- schemes in this file will eventually be moved to their corresponding schemes, e.g., tms, beljaars, ubc, ... when they are fully CCPPized. List all namelist files that were added or changed: ``` A schemes/holtslag_boville/holtslag_boville_diff_options_namelist.xml A schemes/vertical_diffusion/vertical_diffusion_options_namelist.xml - separate "namelist" scheme for HB. - separate "namelist" scheme for diffusion solver. ``` List all files eliminated and why: N/A List all files added and what they do: ``` A to_be_ccppized/coords_1d.meta - metadata file for coords1d DDT. A to_be_ccppized/vdiff_lu_solver.F90 - LU solver for diffusion_solver.F90, moved from CAM A schemes/holtslag_boville/holtslag_boville_diff.F90 A schemes/holtslag_boville/holtslag_boville_diff.meta A schemes/holtslag_boville/holtslag_boville_diff_options.F90 A schemes/holtslag_boville/holtslag_boville_diff_options.meta - CCPPized Holtslag-Boville boundary layer scheme (originally hb_diff.F90) A schemes/holtslag_boville/holtslag_boville_diff_interstitials.F90 A schemes/holtslag_boville/holtslag_boville_diff_interstitials.meta - CCPPized interstitial schemes supporting inputs into HB. A schemes/vertical_diffusion/diffusion_solver.F90 A schemes/vertical_diffusion/diffusion_solver.meta A schemes/vertical_diffusion/vertical_diffusion_options.F90 A schemes/vertical_diffusion/vertical_diffusion_options.meta - CCPPized vertical diffusion solver without support for molecular diffusion (originally diffusion_solver.F90 in CAM) A schemes/vertical_diffusion/diffusion_stubs.F90 A schemes/vertical_diffusion/diffusion_stubs.meta - CCPPized stub for UBC (upper boundary conditions) - CCPPized stub for TMS (Turbulent Mountain Stress) surface orographic drag. - CCPPized stub for Beljaars surface orographic drag. - CCPPized stub for WACCM to use composition-dependent rair in vertical diffusion. A schemes/vertical_diffusion/vertical_diffusion_sponge_layer.F90 A schemes/vertical_diffusion/vertical_diffusion_sponge_layer.meta - CCPPized interstitial for sponge layer kvm (originally in vertical_diffusion_tend) A test/test_suites/suite_vdiff_holtslag_boville.xml - draft SDF for vertical diffusion using HB PBL scheme. ``` List all existing files that have been modified, and describe the changes: (Helpful git command: `git diff --name-status development...`) ``` M to_be_ccppized/coords_1d.F90 - make coords1d lowercase to workaround upstream framework limitation ``` List all automated tests that failed, as well as an explanation for why they weren't fixed: - CAM-SIMA will not build without removing `to_be_ccppized/nlte_fomichev.F90` due to dependency on CAM-only module `chem_surfvals` Is this an answer-changing PR? If so, is it a new physics package, algorithm change, tuning change, etc? B4B with CAM - new physics package If yes to the above question, describe how this code was validated with the new/modified features: * CAM regression tests: bit-for-bit unchanged as of cam6_4_107; compared against Derecho Intel only, Izumi tests will be ran before merging. * CAM-SIMA regression tests using new tests defined in https://github.com/ESCOMP/CAM-SIMA/pull/410 bit-for-bit unchanged, except for Derecho GNU physics tests: - new baselines in hb_vdiff derecho gnu test (new test) - updated baselines in rk_stratiform test due to change in snapshots necessary due to CAM pbuf field change. and updated baselines in CAM4 tests due to adding HB PBL and vertical diffusion solver. --------- Co-authored-by: Jesse Nusbaumer --- .../hack_shallow/hack_convect_shallow.meta | 2 +- .../holtslag_boville_diff.F90 | 718 +++++++++++ .../holtslag_boville_diff.meta | 638 ++++++++++ .../holtslag_boville_diff_interstitials.F90 | 294 +++++ .../holtslag_boville_diff_interstitials.meta | 344 ++++++ .../holtslag_boville_diff_options.F90 | 43 + .../holtslag_boville_diff_options.meta | 37 + ...holtslag_boville_diff_options_namelist.xml | 92 ++ .../diffusion_solver_diagnostics.F90 | 445 +++++++ .../diffusion_solver_diagnostics.meta | 268 +++++ .../holtslag_boville_diff_diagnostics.F90 | 93 ++ .../holtslag_boville_diff_diagnostics.meta | 95 ++ .../vertical_diffusion/diffusion_solver.F90 | 1028 ++++++++++++++++ .../vertical_diffusion/diffusion_solver.meta | 1045 +++++++++++++++++ .../vertical_diffusion/diffusion_stubs.F90 | 332 ++++++ .../vertical_diffusion/diffusion_stubs.meta | 415 +++++++ .../vertical_diffusion_diagnostic_utils.F90 | 284 +++++ .../vertical_diffusion_options.F90 | 42 + .../vertical_diffusion_options.meta | 43 + .../vertical_diffusion_options_namelist.xml | 108 ++ .../vertical_diffusion_sponge_layer.F90 | 138 +++ .../vertical_diffusion_sponge_layer.meta | 87 ++ suites/suite_cam4.xml | 93 +- .../suite_vdiff_holtslag_boville.xml | 90 ++ to_be_ccppized/coords_1d.F90 | 12 +- to_be_ccppized/coords_1d.meta | 7 + to_be_ccppized/vdiff_lu_solver.F90 | 207 ++++ 27 files changed, 6991 insertions(+), 9 deletions(-) create mode 100644 schemes/holtslag_boville/holtslag_boville_diff.F90 create mode 100644 schemes/holtslag_boville/holtslag_boville_diff.meta create mode 100644 schemes/holtslag_boville/holtslag_boville_diff_interstitials.F90 create mode 100644 schemes/holtslag_boville/holtslag_boville_diff_interstitials.meta create mode 100644 schemes/holtslag_boville/holtslag_boville_diff_options.F90 create mode 100644 schemes/holtslag_boville/holtslag_boville_diff_options.meta create mode 100644 schemes/holtslag_boville/holtslag_boville_diff_options_namelist.xml create mode 100644 schemes/sima_diagnostics/diffusion_solver_diagnostics.F90 create mode 100644 schemes/sima_diagnostics/diffusion_solver_diagnostics.meta create mode 100644 schemes/sima_diagnostics/holtslag_boville_diff_diagnostics.F90 create mode 100644 schemes/sima_diagnostics/holtslag_boville_diff_diagnostics.meta create mode 100644 schemes/vertical_diffusion/diffusion_solver.F90 create mode 100644 schemes/vertical_diffusion/diffusion_solver.meta create mode 100644 schemes/vertical_diffusion/diffusion_stubs.F90 create mode 100644 schemes/vertical_diffusion/diffusion_stubs.meta create mode 100644 schemes/vertical_diffusion/vertical_diffusion_diagnostic_utils.F90 create mode 100644 schemes/vertical_diffusion/vertical_diffusion_options.F90 create mode 100644 schemes/vertical_diffusion/vertical_diffusion_options.meta create mode 100644 schemes/vertical_diffusion/vertical_diffusion_options_namelist.xml create mode 100644 schemes/vertical_diffusion/vertical_diffusion_sponge_layer.F90 create mode 100644 schemes/vertical_diffusion/vertical_diffusion_sponge_layer.meta create mode 100644 test/test_suites/suite_vdiff_holtslag_boville.xml create mode 100644 to_be_ccppized/coords_1d.meta create mode 100644 to_be_ccppized/vdiff_lu_solver.F90 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