Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

cam6_3_134: Update atmospheric_physics external #891

Merged
merged 18 commits into from
Oct 31, 2023
Merged
Show file tree
Hide file tree
Changes from 7 commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions Externals_CAM.cfg
Original file line number Diff line number Diff line change
Expand Up @@ -55,9 +55,9 @@ tag = ALI_ARMS_v1.0.1
required = True

[atmos_phys]
tag = atmos_phys0_00_011
branch = updated_standard_names
cacraigucar marked this conversation as resolved.
Show resolved Hide resolved
protocol = git
repo_url = https://github.com/NCAR/atmospheric_physics
repo_url = https://github.com/nusbaume/atmospheric_physics
required = True
local_path = src/atmos_phys

Expand Down
1 change: 1 addition & 0 deletions bld/configure
Original file line number Diff line number Diff line change
Expand Up @@ -2269,6 +2269,7 @@ sub write_filepath
print $fh "$camsrcdir/src/cpl/$cpl\n";
print $fh "$camsrcdir/src/control\n";
print $fh "$camsrcdir/src/utils\n";
print $fh "$camsrcdir/src/utils/cam_ccpp\n";
print $fh "$camsrcdir/src/atmos_phys/utilities\n";


Expand Down
27 changes: 16 additions & 11 deletions src/dynamics/se/dp_coupling.F90
Original file line number Diff line number Diff line change
Expand Up @@ -546,15 +546,17 @@ subroutine derived_phys_dry(phys_state, phys_tend, pbuf2d)
use shr_const_mod, only: shr_const_rwv
use phys_control, only: waccmx_is
use geopotential, only: geopotential_t
use static_energy, only: update_dry_static_energy_run
use check_energy, only: check_energy_timestep_init
use hycoef, only: hyai, ps0
use shr_vmath_mod, only: shr_vmath_log
use qneg_module, only: qneg3
use dyn_tests_utils, only: vc_dry_pressure
use shr_kind_mod, only: shr_kind_cx
! arguments
type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state
type(physics_tend ), intent(inout), dimension(begchunk:endchunk) :: phys_tend
type(physics_buffer_desc), pointer :: pbuf2d(:,:)
type(physics_buffer_desc), pointer :: pbuf2d(:,:)

! local variables
integer :: lchnk
Expand All @@ -564,6 +566,10 @@ subroutine derived_phys_dry(phys_state, phys_tend, pbuf2d)

integer :: m, i, k, ncol, m_cnst
type(physics_buffer_desc), pointer :: pbuf_chnk(:)

!Needed for "update_dry_static_energy" CCPP scheme
integer :: errflg
character(len=shr_kind_cx) :: errmsg
!----------------------------------------------------------------------------

! Evaluate derived quantities
Expand Down Expand Up @@ -686,18 +692,17 @@ subroutine derived_phys_dry(phys_state, phys_tend, pbuf2d)
end do

! Compute initial geopotential heights - based on full pressure
call geopotential_t (phys_state(lchnk)%lnpint, phys_state(lchnk)%lnpmid , phys_state(lchnk)%pint , &
phys_state(lchnk)%pmid , phys_state(lchnk)%pdel , phys_state(lchnk)%rpdel , &
phys_state(lchnk)%t , phys_state(lchnk)%q(:,:,:), rairv(:,:,lchnk), gravit, zvirv , &
phys_state(lchnk)%zi , phys_state(lchnk)%zm , ncol )
call geopotential_t(phys_state(lchnk)%lnpint, phys_state(lchnk)%lnpmid , phys_state(lchnk)%pint, &
phys_state(lchnk)%pmid , phys_state(lchnk)%pdel , phys_state(lchnk)%rpdel , &
phys_state(lchnk)%t , phys_state(lchnk)%q(:,:,:), rairv(:,:,lchnk), gravit, zvirv , &
phys_state(lchnk)%zi , phys_state(lchnk)%zm , ncol)

! Compute initial dry static energy, include surface geopotential
do k = 1, pver
do i = 1, ncol
phys_state(lchnk)%s(i,k) = cpairv(i,k,lchnk)*phys_state(lchnk)%t(i,k) &
+ gravit*phys_state(lchnk)%zm(i,k) + phys_state(lchnk)%phis(i)
end do
end do
call update_dry_static_energy_run(pver, gravit, phys_state(lchnk)%t(1:ncol,:), &
phys_state(lchnk)%zm(1:ncol,:), &
phys_state(lchnk)%phis(1:ncol), &
phys_state(lchnk)%s(1:ncol,:), &
cpairv(1:ncol,:,lchnk), errflg, errmsg)

! Ensure tracers are all positive
call qneg3('D_P_COUPLING',lchnk ,ncol ,pcols ,pver , &
Expand Down
144 changes: 54 additions & 90 deletions src/physics/cam/geopotential.F90
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ module geopotential
! The hydrostatic matrix elements must be consistent with the dynamics algorithm.
! The diagonal element is the itegration weight from interface k+1 to midpoint k.
! The offdiagonal element is the weight between interfaces.
!
!
! Author: B.Boville, Feb 2001 from earlier code by Boville and S.J. Lin
!---------------------------------------------------------------------------------

Expand All @@ -29,17 +29,18 @@ subroutine geopotential_t( &
t , q , rair , gravit , zvir , &
zi , zm , ncol )

!-----------------------------------------------------------------------
!
! Purpose:
! Compute the geopotential height (above the surface) at the midpoints and
!-----------------------------------------------------------------------
!
! Purpose:
! Compute the geopotential height (above the surface) at the midpoints and
! interfaces using the input temperatures and pressures.
!
!-----------------------------------------------------------------------

use ppgrid, only : pcols
use air_composition, only: thermodynamic_active_species_num
use air_composition, only: thermodynamic_active_species_idx, dry_air_species_num
use ppgrid, only: pcols
use constituents, only: pcnst, cnst_get_ind
use ccpp_constituent_prop_mod, only: ccpp_const_props !CCPP constituent properties array (stub version)
cacraigucar marked this conversation as resolved.
Show resolved Hide resolved
use geopotential_temp, only: geopotential_temp_run !CCPP version
!------------------------------Arguments--------------------------------
!
! Input arguments
Expand All @@ -65,6 +66,8 @@ subroutine geopotential_t( &
!
!---------------------------Local variables-----------------------------
!
logical :: lagrang ! Lagrangian vertical coordinate flag
integer :: ixq ! state constituent array index for water vapor
integer :: i,k,idx ! Lon, level indices, water species index
real(r8) :: hkk(ncol) ! diagonal element of hydrostatic matrix
real(r8) :: hkl(ncol) ! off-diagonal element
Expand All @@ -74,28 +77,35 @@ subroutine geopotential_t( &
real(r8) :: qfac(ncol,pver) ! factor to convert from wet to dry mixing ratio
real(r8) :: sum_dry_mixing_ratio(ncol,pver)! sum of dry water mixing ratios

!CCPP-required variables (not used):
integer :: errflg
character(len=512) :: errmsg

!
!-----------------------------------------------------------------------
!
rog(:ncol,:) = rair(:ncol,:) / gravit

! The surface height is zero by definition.

do i = 1,ncol
zi(i,pverp) = 0.0_r8
end do

! Compute zi, zm from bottom up.
! Note, zi(i,k) is the interface above zm(i,k)
!Determine index for water vapor mass mixing ratio
call cnst_get_ind('Q', ixq)

!
! original code for backwards compatability with FV and EUL
!
if (.not.(dycore_is('MPAS') .or. dycore_is('SE'))) then

!dry air gas constant over gravity
rog(:ncol,:) = rair(:ncol,:) / gravit

! The surface height is zero by definition.
do i = 1,ncol
zi(i,pverp) = 0.0_r8
end do

! Compute zi, zm from bottom up.
! Note, zi(i,k) is the interface above zm(i,k)
do k = pver, 1, -1

! First set hydrostatic elements consistent with dynamics

if ((dycore_is('LR') .or. dycore_is('FV3'))) then
do i = 1,ncol
hkl(i) = piln(i,k+1) - piln(i,k)
Expand All @@ -107,83 +117,37 @@ subroutine geopotential_t( &
hkk(i) = 0.5_r8 * hkl(i)
end do
end if

! Now compute tv, zm, zi

do i = 1,ncol
tvfac = 1._r8 + zvir(i,k) * q(i,k,1)
tv = t(i,k) * tvfac

zm(i,k) = zi(i,k+1) + rog(i,k) * tv * hkk(i)
zi(i,k) = zi(i,k+1) + rog(i,k) * tv * hkl(i)
end do
end do
else
!
! For the computation of generalized virtual temperature (equation 16
! in Lauritzen et al. (2018); https://doi.org/10.1029/2017MS001257)
!
! Compute factor for converting wet to dry mixing ratio (eq.7)
!
qfac = 1.0_r8
do idx = dry_air_species_num + 1,thermodynamic_active_species_num
do k = 1,pver
do i = 1,ncol
qfac(i,k) = qfac(i,k)-q(i,k,thermodynamic_active_species_idx(idx))
end do
end do
end do
qfac = 1.0_r8/qfac

! Compute sum of dry water mixing ratios
sum_dry_mixing_ratio = 1.0_r8
do idx = dry_air_species_num + 1,thermodynamic_active_species_num
do k = 1,pver
do i = 1,ncol
sum_dry_mixing_ratio(i,k) = sum_dry_mixing_ratio(i,k)&
+q(i,k,thermodynamic_active_species_idx(idx))*qfac(i,k)
end do
end do
end do
sum_dry_mixing_ratio(:,:) = 1.0_r8/sum_dry_mixing_ratio(:,:)
! Compute zi, zm from bottom up.
! Note, zi(i,k) is the interface above zm(i,k)
do k = pver, 1, -1

! First set hydrostatic elements consistent with dynamics

!
! the outcommented code is left for when/if we will support
! FV3 and/or FV with condensate loading
!

! if ((dycore_is('LR') .or. dycore_is('FV3'))) then
! do i = 1,ncol
! hkl(i) = piln(i,k+1) - piln(i,k)
! hkk(i) = 1._r8 - pint(i,k) * hkl(i) * rpdel(i,k)
! end do
! else!MPAS, SE or EUL
!
! For SE : pmid = 0.5*(pint(k+1)+pint(k))
! For MPAS : pmid is computed from theta_m, rhodry, etc.
!
do i = 1,ncol
hkl(i) = pdel(i,k) / pmid(i,k)
hkk(i) = 0.5_r8 * hkl(i)
end do
! end if


! Now compute tv, zm, zi

do i = 1,ncol
tvfac = (1._r8 + (zvir(i,k)+1.0_r8) * q(i,k,1)*qfac(i,k))*sum_dry_mixing_ratio(i,k)
tvfac = 1._r8 + zvir(i,k) * q(i,k,ixq)
cacraigucar marked this conversation as resolved.
Show resolved Hide resolved
tv = t(i,k) * tvfac

zm(i,k) = zi(i,k+1) + rog(i,k) * tv * hkk(i)
zi(i,k) = zi(i,k+1) + rog(i,k) * tv * hkl(i)
end do
end do
else !Using MPAS or SE dycore

!Determine vertical coordinate type,
!NOTE: Currently the FV (LR) or FV3 dycores
! do not allow for condensate loading,
! so for now 'lagrang' will always be FALSE.
if ((dycore_is('LR') .or. dycore_is('FV3'))) then
lagrang = .true.
else
lagrang = .false.
end if

!Use CCPP version of geopotential_t:
call geopotential_temp_run(pver, lagrang, pver, 1, pverp, 1, &
pcnst, piln(1:ncol,:), pint(1:ncol,:), pmid(1:ncol,:), &
pdel(1:ncol,:), rpdel(1:ncol,:), t(1:ncol,:), &
q(1:ncol,:,ixq), q(1:ncol,:,:), ccpp_const_props, &
rair(1:ncol,:), gravit, zvir(1:ncol,:), zi(1:ncol,:), &
zm(1:ncol,:), ncol, errflg, errmsg)

end if
return
end subroutine geopotential_t
end module geopotential
12 changes: 9 additions & 3 deletions src/physics/cam/physpkg.F90
Original file line number Diff line number Diff line change
Expand Up @@ -786,6 +786,8 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out )
use phys_grid_ctem, only: phys_grid_ctem_init
use cam_budget, only: cam_budget_init

use ccpp_constituent_prop_mod, only: ccpp_const_props_init

! Input/output arguments
type(physics_state), pointer :: phys_state(:)
type(physics_tend ), pointer :: phys_tend(:)
Expand Down Expand Up @@ -978,6 +980,10 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out )

end if

! Initialize stub CCPP constituent properties array
cacraigucar marked this conversation as resolved.
Show resolved Hide resolved
! for use in CCPP-ized physics schemes:
call ccpp_const_props_init()

! Initialize qneg3 and qneg4
call qneg_init()

Expand All @@ -989,7 +995,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out )

! Initialize the budget capability
call cam_budget_init()

! addfld calls for U, V tendency budget variables that are output in
! tphysac, tphysbc
call addfld ( 'UTEND_DCONV', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by deep convection')
Expand Down Expand Up @@ -1938,7 +1944,7 @@ subroutine tphysac (ztodt, cam_in, &
else
!
! for moist-mixing ratio based dycores
!
!
! Note: this operation will NOT be reverted with set_wet_to_dry after set_dry_to_wet call
!
call set_dry_to_wet(state)
Expand All @@ -1960,7 +1966,7 @@ subroutine tphysac (ztodt, cam_in, &
if (vc_dycore == vc_height.or.vc_dycore == vc_dry_pressure) then
!
! MPAS and SE specific scaling of temperature for enforcing energy consistency
! (and to make sure that temperature dependent diagnostic tendencies
! (and to make sure that temperature dependent diagnostic tendencies
! are computed correctly; e.g. dtcore)
!
scaling(1:ncol,:) = cpairv(:ncol,:,lchnk)/cp_or_cv_dycore(:ncol,:,lchnk)
Expand Down
8 changes: 4 additions & 4 deletions src/physics/cam/ref_pres.F90
Original file line number Diff line number Diff line change
@@ -1,14 +1,14 @@
module ref_pres
!--------------------------------------------------------------------------
!
!
! Provides access to reference pressures for use by the physics
! parameterizations. The pressures are provided by the dynamical core
! since it determines the grid used by the physics.
!
!
! Note that the init method for this module is called before the init
! method in physpkg; therefore, most physics modules can use these
! reference pressures during their init phases.
!
!
!--------------------------------------------------------------------------

use shr_kind_mod, only: r8=>shr_kind_r8
Expand All @@ -25,7 +25,7 @@ module ref_pres
! surface pressure ('eta' coordinate)

real(r8), protected :: ptop_ref ! Top of model
real(r8), protected :: psurf_ref ! Surface pressure
real(r8), protected :: psurf_ref ! reference pressure

! Number of top levels using pure pressure representation
integer, protected :: num_pr_lev
Expand Down
6 changes: 6 additions & 0 deletions src/physics/cam_dev/physpkg.F90
Original file line number Diff line number Diff line change
Expand Up @@ -769,6 +769,8 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out )
use cam_budget, only: cam_budget_init
use phys_grid_ctem, only: phys_grid_ctem_init

use ccpp_constituent_prop_mod, only: ccpp_const_props_init

! Input/output arguments
type(physics_state), pointer :: phys_state(:)
type(physics_tend ), pointer :: phys_tend(:)
Expand Down Expand Up @@ -947,6 +949,10 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out )

end if

! Initialize stub CCPP constituent properties array
cacraigucar marked this conversation as resolved.
Show resolved Hide resolved
! for use in CCPP-ized physics schemes:
call ccpp_const_props_init()

! Initialize qneg3 and qneg4
call qneg_init()

Expand Down
Loading