cam6_3_134: Update atmospheric_physics external

Externals_CAM.cfg

@@ -55,7 +55,7 @@ tag = ALI_ARMS_v1.0.1
 required = True
 
 [atmos_phys]
-tag = atmos_phys0_00_011
+tag = atmos_phys0_01_000
 protocol = git
 repo_url = https://github.com/NCAR/atmospheric_physics
 required = True

bld/configure

@@ -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";
 
 

src/dynamics/se/dp_coupling.F90

@@ -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
@@ -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
@@ -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    , &

src/physics/cam/geopotential.F90

@@ -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
 !---------------------------------------------------------------------------------
 
@@ -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 (CAM version)
+use geopotential_temp,         only: geopotential_temp_run !CCPP version
 !------------------------------Arguments--------------------------------
 !
 ! Input arguments
@@ -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
@@ -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)
@@ -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)
           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

src/physics/cam/physpkg.F90

@@ -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(:)
@@ -978,6 +980,10 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out )
 
     end if
 
+    ! Initialize CAM CCPP constituent properties array
+    ! for use in CCPP-ized physics schemes:
+    call ccpp_const_props_init()
+
     ! Initialize qneg3 and qneg4
     call qneg_init()
 
@@ -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')
@@ -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)
@@ -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)

src/physics/cam/ref_pres.F90

@@ -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
@@ -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

src/physics/cam_dev/physpkg.F90

@@ -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(:)
@@ -947,6 +949,10 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out )
 
     end if
 
+    ! Initialize CAM CCPP constituent properties array
+    ! for use in CCPP-ized physics schemes:
+    call ccpp_const_props_init()
+
     ! Initialize qneg3 and qneg4
     call qneg_init()