diff --git a/src/data/air_composition.F90 b/src/data/air_composition.F90 index d69685a9..0da0a5b2 100644 --- a/src/data/air_composition.F90 +++ b/src/data/air_composition.F90 @@ -1,4 +1,5 @@ -! air_composition module defines major species of the atmosphere and manages the physical properties that are dependent on the composition of air +! air_composition module defines major species of the atmosphere and manages +! the physical properties that are dependent on the composition of air module air_composition use ccpp_kinds, only: kind_phys @@ -13,7 +14,8 @@ module air_composition save public :: air_composition_init - public :: air_composition_update + public :: dry_air_composition_update + public :: water_composition_update ! get_cp_dry: (generalized) heat capacity for dry air public :: get_cp_dry ! get_cp: (generalized) heat capacity @@ -511,54 +513,66 @@ end subroutine air_composition_init !=========================================================================== !----------------------------------------------------------------------- - ! air_composition_update: Update the physics "constants" that vary + ! dry_air_composition_update: Update the physics "constants" that vary !------------------------------------------------------------------------- !=========================================================================== - subroutine air_composition_update(mmr, ncol, vcoord, to_moist_factor) - use dyn_tests_utils, only: vc_height, vc_moist_pressure, vc_dry_pressure + subroutine dry_air_composition_update(mmr, ncol, to_dry_factor) - real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array + !(mmr = dry mixing ratio, if not, use to_dry_factor to convert!) + real(kind_phys), intent(in) :: mmr(:,:,:) ! mixing ratios for species dependent dry air integer, intent(in) :: ncol ! number of columns - integer, intent(in) :: vcoord ! vertical coordinate specifier for dycore - real(kind_phys), optional, intent(in) :: to_moist_factor(:,:) - - character(len=*), parameter :: subname = 'air_composition_update' + real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) call get_R_dry(mmr(:ncol, :, :), thermodynamic_active_species_idx, & - rairv(:ncol, :), fact=to_moist_factor) + rairv(:ncol, :), fact=to_dry_factor) call get_cp_dry(mmr(:ncol,:,:), thermodynamic_active_species_idx, & - cpairv(:ncol,:), fact=to_moist_factor) + cpairv(:ncol,:), fact=to_dry_factor) call get_mbarv(mmr(:ncol,:,:), thermodynamic_active_species_idx, & - mbarv(:ncol,:), fact=to_moist_factor) + mbarv(:ncol,:), fact=to_dry_factor) cappav(:ncol,:) = rairv(:ncol,:) / cpairv(:ncol,:) + end subroutine air_composition_update + + !=========================================================================== + !--------------------------------------------------------------------------- + ! water_composition_update: Update generalized cp or cv depending on dycore + !--------------------------------------------------------------------------- + !=========================================================================== + subroutine water_composition_update(mmr, ncol, vcoord, to_dry_factor) + use dyn_tests_utils, only: vc_height, vc_moist_pressure, vc_dry_pressure + + real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: vcoord ! vertical coordinate specifier for dycore + real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) + + character(len=*), parameter :: subname = 'water_composition_update' + ! update enthalpy or internal energy scaling factor for energy consistency with CAM physics - ! FIXME hplin 9/5/24: why is air_composition_update in CAM-SIMA to_moist_factor - ! whereas water_composition_update in CAM to_dry_factor. - ! get_cp_dry implementation appears to be the same?, but - ! get_cp implementation appears to take an 1/factor. - ! Nothing in CAM-SIMA passes to_moist_factor to cam_thermo_update for now, so it is not possible - ! to know which is correct. In CAM, physpkg.F90 calls cam_thermo_water_update with to_dry_factor=pdel/pdeldry. + ! cp_or_cv_dycore is now a registry variable (physics_types) in CAM-SIMA instead of in-module + if (vcoord == vc_moist_pressure) then ! FV: moist pressure vertical coordinate does not need update. else if (vcoord == vc_dry_pressure) then ! SE call get_cp(mmr(:ncol,:,:), .false., cp_or_cv_dycore(:ncol,:), & - dp_dry=to_moist_factor, active_species_idx_dycore=thermodynamic_active_species_idx) + factor=to_dry_factor, active_species_idx_dycore=thermodynamic_active_species_idx, & + cpdry=cpairv(:ncol,:)) else if (vcoord == vc_height) then ! MPAS - ! Note hplin 9/5/24: the below code is "unused and untested" and reproduced verbatim (together with this warning) - ! from CAM air_composition.F90 - call get_R(mmr(:ncol,:,:), thermodynamic_active_species_idx, cp_or_cv_dycore(:ncol,:), & - fact=to_moist_factor, R_dry=rairv(:ncol,:)) + call get_R(mmr(:ncol,:,:), thermodynamic_active_species_idx, & + cp_or_cv_dycore(:ncol,:), fact=to_dry_factor, Rdry=rairv(:ncol,:)) + + ! internal energy coefficient for MPAS + ! (equation 92 in Eldred et al. 2023; https://rmets.onlinelibrary.wiley.com/doi/epdf/10.1002/qj.4353) cp_or_cv_dycore(:ncol,:) = cp_or_cv_dycore(:ncol,:) * (cpairv(:ncol,:) - rairv(:ncol,:)) / rairv(:ncol,:) else call endrun(subname//': dycore vcoord not supported') end if - end subroutine air_composition_update + end subroutine water_composition_update !=========================================================================== !*************************************************************************** @@ -668,27 +682,33 @@ end subroutine get_cp_dry_2hd ! !*************************************************************************** ! - subroutine get_cp_1hd(tracer, inv_cp, cp, dp_dry, active_species_idx_dycore) + subroutine get_cp_1hd(tracer, inv_cp, cp, factor, active_species_idx_dycore, cpdry) use cam_abortutils, only: endrun use string_utils, only: to_str ! Dummy arguments ! tracer: Tracer array + ! + ! factor not present then tracer must be dry mixing ratio + ! if factor present tracer*factor must be dry mixing ratio + ! real(kind_phys), intent(in) :: tracer(:,:,:) - real(kind_phys), optional, intent(in) :: dp_dry(:,:) ! inv_cp: output inverse cp instead of cp logical, intent(in) :: inv_cp real(kind_phys), intent(out) :: cp(:,:) + ! dp: if provided then tracer is mass not mixing ratio + real(kind_phys), optional, intent(in) :: factor(:,:) ! active_species_idx_dycore: array of indices for index of ! thermodynamic active species in dycore tracer array ! (if different from physics index) integer, optional, intent(in) :: active_species_idx_dycore(:) + real(kind_phys), optional, intent(in) :: cpdry(:,:) ! Local variables integer :: qdx, itrac real(kind_phys) :: sum_species(SIZE(cp, 1), SIZE(cp, 2)) real(kind_phys) :: sum_cp(SIZE(cp, 1), SIZE(cp, 2)) - real(kind_phys) :: factor(SIZE(cp, 1), SIZE(cp, 2)) + real(kind_phys) :: factor_local(SIZE(cp, 1), SIZE(cp, 2)) integer :: idx_local(thermodynamic_active_species_num) character(len=*), parameter :: subname = 'get_cp_1hd: ' @@ -704,28 +724,38 @@ subroutine get_cp_1hd(tracer, inv_cp, cp, dp_dry, active_species_idx_dycore) idx_local = thermodynamic_active_species_idx end if - if (present(dp_dry)) then - factor = 1.0_kind_phys / dp_dry + if (present(factor)) then + factor_local = factor else - factor = 1.0_kind_phys + factor_local = 1.0_kind_phys end if + sum_species = 1.0_kind_phys ! all dry air species sum to 1 do qdx = dry_air_species_num + 1, thermodynamic_active_species_num itrac = idx_local(qdx) - sum_species(:,:) = sum_species(:,:) + & - (tracer(:,:,itrac) * factor(:,:)) + sum_species(:,:) = sum_species(:,:) + (tracer(:,:,itrac) * factor_local(:,:)) end do - ! Get heat capacity at constant pressure (Cp) for dry air: - call get_cp_dry(tracer, idx_local, sum_cp, fact=factor) + if (dry_air_species_num == 0) then + ! FIXME hplin 9/13/24: this if-clause is always in CAM but was not in CAM-SIMA + sum_cp = thermodynamic_active_species_cp(0) + else if (present(cpdry)) then + ! + ! if cpdry is known don't recompute + ! + sum_cp = cpdry + else + ! Get heat capacity at constant pressure (Cp) for dry air: + call get_cp_dry(tracer, idx_local, sum_cp, fact=factor_local) + end if ! Add water species to Cp: do qdx = dry_air_species_num + 1, thermodynamic_active_species_num itrac = idx_local(qdx) sum_cp(:,:) = sum_cp(:,:) + & - (thermodynamic_active_species_cp(qdx) * tracer(:,:,itrac) * & - factor(:,:)) + (thermodynamic_active_species_cp(qdx) * tracer(:,:,itrac) * factor_local(:,:)) end do + if (inv_cp) then cp = sum_species / sum_cp else @@ -736,7 +766,7 @@ end subroutine get_cp_1hd !=========================================================================== - subroutine get_cp_2hd(tracer, inv_cp, cp, dp_dry, active_species_idx_dycore) + subroutine get_cp_2hd(tracer, inv_cp, cp, factor, active_species_idx_dycore, cpdry) ! Version of get_cp for arrays that have a second horizontal index use cam_abortutils, only: endrun use string_utils, only: to_str @@ -744,14 +774,15 @@ subroutine get_cp_2hd(tracer, inv_cp, cp, dp_dry, active_species_idx_dycore) ! Dummy arguments ! tracer: Tracer array real(kind_phys), intent(in) :: tracer(:,:,:,:) - real(kind_phys), optional, intent(in) :: dp_dry(:,:,:) ! inv_cp: output inverse cp instead of cp logical, intent(in) :: inv_cp real(kind_phys), intent(out) :: cp(:,:,:) + real(kind_phys), optional, intent(in) :: factor(:,:,:) ! active_species_idx_dycore: array of indicies for index of ! thermodynamic active species in dycore tracer array ! (if different from physics index) integer, optional, intent(in) :: active_species_idx_dycore(:) + real(kind_phys), optional, intent(in) :: cpdry(:,:,:) ! Local variables integer :: jdx @@ -759,11 +790,17 @@ subroutine get_cp_2hd(tracer, inv_cp, cp, dp_dry, active_species_idx_dycore) character(len=*), parameter :: subname = 'get_cp_2hd: ' do jdx = 1, SIZE(cp, 2) - if (present(dp_dry)) then - call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :), & - dp_dry=dp_dry(:, jdx, :), active_species_idx_dycore=active_species_idx_dycore) + if (present(factor).and.present(cpdry)) then + call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),& + factor=factor(:, jdx, :), active_species_idx_dycore=active_species_idx_dycore, cpdry=cpdry(:,jdx,:)) + else if (present(factor)) then + call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),& + factor=factor(:, jdx, :), active_species_idx_dycore=active_species_idx_dycore) + else if (present(cpdry)) then + call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),& + active_species_idx_dycore=active_species_idx_dycore, cpdry=cpdry(:,jdx,:)) else - call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :), & + call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),& active_species_idx_dycore=active_species_idx_dycore) end if end do @@ -872,9 +909,10 @@ end subroutine get_R_dry_2hd ! !*************************************************************************** ! - subroutine get_R_1hd(tracer, active_species_idx, R, fact) + subroutine get_R_1hd(tracer, active_species_idx, R, fact, Rdry) use cam_abortutils, only: endrun use string_utils, only: to_str + use physconst, only: rair ! Dummy arguments ! tracer: !tracer array @@ -885,6 +923,7 @@ subroutine get_R_1hd(tracer, active_species_idx, R, fact) real(kind_phys), intent(out) :: R(:, :) ! fact: optional factor for converting tracer to dry mixing ratio real(kind_phys), optional, intent(in) :: fact(:, :) + real(kind_phys), optional, intent(in) :: Rdry(:, :) ! Local variables integer :: qdx, itrac @@ -903,12 +942,19 @@ subroutine get_R_1hd(tracer, active_species_idx, R, fact) call endrun(subname//"SIZE mismatch in dimension 2 "// & to_str(SIZE(fact, 2))//' /= '//to_str(SIZE(factor, 2))) end if - call get_R_dry(tracer, active_species_idx, R, fact=fact) factor = fact(:,:) else - call get_R_dry(tracer, active_species_idx, R) factor = 1.0_kind_phys end if + + if (dry_air_species_num == 0) then + R = rair + else if (present(Rdry)) then + R = Rdry + else + call get_R_dry(tracer, active_species_idx, R, fact=factor) + end if + idx_local = active_species_idx sum_species = 1.0_kind_phys ! all dry air species sum to 1 do qdx = dry_air_species_num + 1, thermodynamic_active_species_num @@ -963,7 +1009,7 @@ end subroutine get_R_2hd !************************************************************************************************************************* ! subroutine get_mbarv_1hd(tracer, active_species_idx, mbarv_in, fact) - use physconst, only: mwdry, rair, cpair + use physconst, only: mwdry real(kind_phys), intent(in) :: tracer(:,:,:) !tracer array integer, intent(in) :: active_species_idx(:) !index of active species in tracer real(kind_phys), intent(out) :: mbarv_in(:,:) !molecular weight of dry air diff --git a/src/data/cam_thermo.F90 b/src/data/cam_thermo.F90 index 233f92c7..2b0af969 100644 --- a/src/data/cam_thermo.F90 +++ b/src/data/cam_thermo.F90 @@ -33,8 +33,10 @@ module cam_thermo ! cam_thermo_init: Initialize constituent dependent properties public :: cam_thermo_init - ! cam_thermo_update: Update constituent dependent properties - public :: cam_thermo_update + ! cam_thermo_dry_air_update: Update dry air composition dependent properties + public :: cam_thermo_dry_air_update + ! cam_thermo_water_update: Update water dependent properties + public :: cam_thermo_water_update ! get_enthalpy: enthalpy quantity = dp*cp*T public :: get_enthalpy ! get_virtual_temp: virtual temperature @@ -208,51 +210,67 @@ end subroutine cam_thermo_init !=========================================================================== + ! !*************************************************************************** ! - ! cam_thermo_update: update species dependent constants for physics + ! cam_thermo_dry_air_update: update dry air species dependent constants for physics ! !*************************************************************************** ! - subroutine cam_thermo_update(mmr, T, ncol, update_thermo_variables, to_moist_factor) - use air_composition, only: air_composition_update, update_zvirv - use string_utils, only: to_str - !----------------------------------------------------------------------- - ! Update the physics "constants" that vary - !------------------------------------------------------------------------- + subroutine cam_thermo_dry_air_update(mmr, T, ncol, update_thermo_variables, to_dry_factor) + use air_composition, only: dry_air_composition_update + use air_composition, only: update_zvirv + use string_utils, only: int2str - !------------------------------Arguments---------------------------------- + real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array (mmr = dry mixing ratio, if not use to_dry_factor to convert) + real(kind_phys), intent(in) :: T(:,:) ! temperature + integer, intent(in) :: ncol ! number of columns + logical, intent(in) :: update_thermo_variables ! true: calculate composition-dependent thermo variables + ! false: do not calculate composition-dependent thermo variables + real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) ! if mmr moist convert - real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array - real(kind_phys), intent(in) :: T(:,:) ! temperature - integer, intent(in) :: ncol ! number of columns - logical, intent(in) :: update_thermo_variables ! true: calculate composition-dependent thermo variables - ! false: do not calculate composition-dependent thermo variables - - real(kind_phys), optional, intent(in) :: to_moist_factor(:,:) - ! - !---------------------------Local storage------------------------------- - real(kind_phys):: sponge_factor(SIZE(mmr, 2)) - character(len=*), parameter :: subname = 'cam_thermo_update: ' + ! Local vars + real(kind_phys) :: sponge_factor(SIZE(mmr, 2)) + character(len=*), parameter :: subname = 'cam_thermo_dry_air_update: ' - if (.not. update_thermo_variables) then - return + if (present(to_dry_factor)) then + if (SIZE(to_dry_factor, 1) /= ncol) then + call endrun(subname//'DIM 1 of to_dry_factor is'//to_str(SIZE(to_dry_factor,1))//'but should be'//to_str(ncol)) + end if end if - if (present(to_moist_factor)) then - if (SIZE(to_moist_factor, 1) /= ncol) then - call endrun(subname//'DIM 1 of to_moist_factor is'//to_str(SIZE(to_moist_factor,1))//'but should be'//to_str(ncol)) - end if - end if sponge_factor = 1.0_kind_phys - call air_composition_update(mmr, ncol, to_moist_factor=to_moist_factor) + call dry_air_composition_update(mmr, ncol, to_dry_factor=to_dry_factor) call get_molecular_diff_coef(T(:ncol,:), .true., sponge_factor, kmvis(:ncol,:), & - kmcnd(:ncol,:), tracer=mmr(:ncol,:,:), fact=to_moist_factor, & + kmcnd(:ncol,:), tracer=mmr(:ncol,:,:), fact=to_dry_factor, & active_species_idx_dycore=thermodynamic_active_species_idx) + + ! Calculate zvirv for WACCM-X. call update_zvirv() - end subroutine cam_thermo_update + end subroutine cam_thermo_dry_air_update + + ! + !*************************************************************************** + ! + ! cam_thermo_water_update: update water species dependent constants for physics + ! + !*************************************************************************** + ! + subroutine cam_thermo_water_update(mmr, ncol, vcoord, to_dry_factor) + use air_composition, only: water_composition_update + !----------------------------------------------------------------------- + ! Update the physics "constants" that vary + !------------------------------------------------------------------------- + + real(kind_phys), intent(in) :: mmr(:,:,:) ! constituents array + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: vcoord + real(kind_phys), optional, intent(in) :: to_dry_factor(:,:) + + call water_composition_update(mmr, ncol, vcoord, to_dry_factor=to_dry_factor) + end subroutine cam_thermo_water_update !=========================================================================== diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index 8b56e9d9..64efb9b2 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -583,7 +583,10 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) use runtime_obj, only: wv_stdname use physics_types, only: lagrangian_vertical use physconst, only: cpair, gravit, zvir, cappa - use cam_thermo, only: cam_thermo_update + use cam_thermo, only: cam_thermo_dry_air_update, cam_thermo_water_update + use air_composition, only: thermodynamic_active_species_num + use air_composition, only: thermodynamic_active_species_idx + use air_composition, only: dry_air_species_num use physics_types, only: cpairv, rairv, zvirv use physics_grid, only: columns_on_task use geopotential_temp, only: geopotential_temp_run @@ -607,7 +610,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) !constituent properties pointer type(ccpp_constituent_prop_ptr_t), pointer :: const_prop_ptr(:) - integer :: m, i, k + integer :: m, i, k, m_cnst integer :: ix_q !Needed for "geopotential_temp" CCPP scheme @@ -622,6 +625,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) real(r8) :: mmrSum_O_O2_H ! Sum of mass mixing ratios for O, O2, and H real(r8), parameter :: mmrMin=1.e-20_r8 ! lower limit of o2, o, and h mixing ratios real(r8), parameter :: N2mmrMin=1.e-6_r8 ! lower limit of o2, o, and h mixing ratios + real(r8), parameter :: H2lim=6.e-5_r8 ! H2 limiter: 10x global H2 MMR (Roble, 1995) !---------------------------------------------------------------------------- ! Nullify pointers @@ -683,14 +687,23 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) end do ! wet pressure variables (should be removed from physics!) + factor_array(:,:) = 1.0_kind_phys + !$omp parallel do num_threads(horz_num_threads) private (k, i, m_cnst) + do m_cnst = dry_air_species_num + 1, thermodynamic_active_species_num + m = thermodynamic_active_species_idx(m_cnst) + ! for now m is just == ix_q (one entry) + do k = 1, nlev + do i = 1, pcols + ! at this point all q's are dry + factor_array(i,k) = factor_array(i,k) + const_data_ptr(i,k,m) + end do + end do + end do !$omp parallel do num_threads(horz_num_threads) private (k, i) - do k=1,nlev - do i=1, pcols - ! to be consistent with total energy formula in physic's check_energy module only - ! include water vapor in moist dp - factor_array(i,k) = 1._kind_phys+const_data_ptr(i,k,ix_q) - phys_state%pdel(i,k) = phys_state%pdeldry(i,k)*factor_array(i,k) + do k = 1, nlev + do i = 1, pcols + phys_state%pdel(i,k) = phys_state%pdeldry(i,k) * factor_array(i,k) end do end do @@ -721,31 +734,12 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) do k = 1, nlev do i = 1, pcols phys_state%rpdel(i,k) = 1._kind_phys/phys_state%pdel(i,k) - phys_state%exner(i,k) = (phys_state%pint(i,pver+1)/phys_state%pmid(i,k))**cappa - end do - end do - - ! all tracers (including moisture) are in dry mixing ratio units - ! physics expect water variables moist - factor_array(:,1:nlev) = 1._kind_phys/factor_array(:,1:nlev) - - !$omp parallel do num_threads(horz_num_threads) private (m, k, i) - do m=1, num_advected - do k = 1, nlev - do i=1, pcols - !This should ideally check if a constituent is a wet - !mixing ratio or not, but until that is working properly - !in the CCPP framework just check for the water species status - !instead, which is all that CAM physics requires: - if (const_is_water_species(m)) then - const_data_ptr(i,k,m) = factor_array(i,k)*const_data_ptr(i,k,m) - end if - end do end do end do !------------------------------------------------------------ - ! Ensure O2 + O + H (N2) mmr greater than one. + ! Apply limiters to mixing ratios of major species (waccmx): + ! Ensure N2 = 1 - (O2 + O + H) mmr is greater than 0 ! Check for unusually large H2 values and set to lower value. !------------------------------------------------------------ if (cam_runtime_opts%waccmx_option() == 'ionosphere' .or. & @@ -769,8 +763,8 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) endif - if(const_data_ptr(i,k,ixh2) .gt. 6.e-5_r8) then - const_data_ptr(i,k,ixh2) = 6.e-5_r8 + if(const_data_ptr(i,k,ixh2) > H2lim) then + const_data_ptr(i,k,ixh2) = H2lim endif end do @@ -789,9 +783,60 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend) ! Compute molecular viscosity(kmvis) and conductivity(kmcnd). ! Update zvirv registry variable; calculated for WACCM-X. !----------------------------------------------------------------------------- + if (dry_air_species_num > 0) then + call cam_thermo_dry_air_update( & + mmr = const_data_ptr, & ! dry MMR + T = phys_state%t, & + ncol = pcols, & + update_thermo_variables = cam_runtime_opts%update_thermodynamic_variables() & + ) + else + zvirv(:,:) = zvir + end if - call cam_thermo_update(const_data_ptr, phys_state%t, pcols, & - cam_runtime_opts%update_thermodynamic_variables()) + ! + ! update cp_or_cv_dycore in module air_composition. + ! (note: at this point q is dry) + ! + call cam_thermo_water_update( & + mmr = const_data_ptr, & ! dry MMR + ncol = pcols, & + vcoord = vc_dry_pressure & + ) + + !$omp parallel do num_threads(horz_num_threads) private (k, i) + do k = 1, nlev + do i = 1, pcols + ! CAM version: + ! phys_state%exner(i,k) = (phys_state%pint(i,pver+1)/phys_state%pmid(i,k))**cappav(i,k) + + ! CAM-SIMA version (uses constant cappa): + phys_state%exner(i,k) = (phys_state%pint(i,pver+1)/phys_state%pmid(i,k))**cappa + end do + end do + + ! ========= Q is dry ^^^ ---- Q is moist vvv ========= ! + + ! + ! CAM physics expects that: water tracers (including moisture) are moist; the rest dry mixing ratio + ! at this point Q is converted to moist. + ! + factor_array(:,1:nlev) = 1._kind_phys/factor_array(:,1:nlev) + + !$omp parallel do num_threads(horz_num_threads) private (m, k, i) + do m = 1, num_advected + do k = 1, nlev + do i = 1, pcols + ! This should ideally check if a constituent is a wet + ! mixing ratio or not, but until that is working properly + ! in the CCPP framework just check for the water species status + ! instead, which is all that CAM physics requires: + if (const_is_water_species(m)) then + const_data_ptr(i,k,m) = factor_array(i,k)*const_data_ptr(i,k,m) + end if + end do + end do + end do !Call geopotential_temp CCPP scheme: call geopotential_temp_run(pver, lagrangian_vertical, pver, 1, &