Address review comments

cime_config/config_component.xml

@@ -162,7 +162,6 @@
 
       <!-- Simple models -->
       <value compset="_CAM%ADIAB">--physics-suites adiabatic</value>
-<!--      <value compset="_CAM%DABIP04">-phys adiabatic</value> -->
       <value compset="_CAM%TJ16">--physics-suites tj2016 --analytic_ic</value>
 <!--      <value compset="_CAM%KESSLER">-phys kessler -chem terminator -analytic_ic</value> -->
       <value compset="_CAM%KESSLER">--physics-suites kessler --analytic_ic</value>

src/data/air_composition.F90

@@ -544,6 +544,7 @@ end subroutine dry_air_composition_update
    !===========================================================================
    subroutine water_composition_update(mmr, ncol, energy_formula, to_dry_factor)
       use cam_thermo_formula, only: ENERGY_FORMULA_DYCORE_FV, ENERGY_FORMULA_DYCORE_SE, ENERGY_FORMULA_DYCORE_MPAS
+      use string_utils,       only: stringify
 
       real(kind_phys),           intent(in) :: mmr(:,:,:)     ! constituents array
       integer,                   intent(in) :: ncol           ! number of columns
@@ -559,24 +560,21 @@ subroutine water_composition_update(mmr, ncol, energy_formula, to_dry_factor)
          ! FV: moist pressure vertical coordinate does not need update.
       else if (energy_formula == ENERGY_FORMULA_DYCORE_SE) then
          ! SE
-
-         ! **TEMP** TODO hplin 9/17/24:
-         ! for compatibility with CAM-SIMA that allocates thermodynamic_active_species_idx(0:num_advected)
-         ! (whereas CAM only allocates 1-indexed) I subset it here. This should be verified during code
-         ! review.
+         ! Note: species index subset to 1: because SIMA currently uses index 0. See #334.
          call get_cp(mmr(:ncol,:,:), .false., cp_or_cv_dycore(:ncol,:), &
                      factor=to_dry_factor, active_species_idx_dycore=thermodynamic_active_species_idx(1:), &
                      cpdry=cpairv(:ncol,:))
       else if (energy_formula == ENERGY_FORMULA_DYCORE_MPAS) then
          ! MPAS
+         ! Note: species index subset to 1: because SIMA currently uses index 0. See #334.
          call get_R(mmr(:ncol,:,:), thermodynamic_active_species_idx(1:), &
                     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)
+         ! (equation 92 in Eldred et al. 2023; doi: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 energy formula not supported')
+         call endrun(subname//': dycore energy formula (value = '//stringify((/energy_formula/))//') not supported')
       end if
 
    end subroutine water_composition_update
@@ -696,14 +694,14 @@ subroutine get_cp_1hd(tracer, inv_cp, cp, factor, active_species_idx_dycore, cpd
       ! 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
+      ! if factor not present then tracer must be a dry mixing ratio
+      ! if factor present tracer*factor must be a dry mixing ratio
       !
       real(kind_phys),           intent(in)  :: tracer(:,:,:)
       ! 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
+      ! if provided then tracer is not a mass 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

src/data/cam_thermo.F90

@@ -221,14 +221,14 @@ end subroutine cam_thermo_init
    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: to_str
+      use string_utils,    only: stringify
 
       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), optional, intent(in) :: to_dry_factor(:,:) ! if mmr wet or moist convert
 
       ! Local vars
       real(kind_phys) :: sponge_factor(SIZE(mmr, 2))
@@ -240,7 +240,7 @@ subroutine cam_thermo_dry_air_update(mmr, T, ncol, update_thermo_variables, to_d
 
       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))
+          call endrun(subname//'DIM 1 of to_dry_factor is'//stringify((/SIZE(to_dry_factor,1)/))//'but should be'//stringify((/ncol/)))
         end if
       end if
 
@@ -263,17 +263,29 @@ end subroutine cam_thermo_dry_air_update
     !
     !***************************************************************************
     !
-    subroutine cam_thermo_water_update(mmr, ncol, energy_formula, to_dry_factor)
+    subroutine cam_thermo_water_update(mmr, ncol, pver, energy_formula, 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) :: pver       ! number of vertical levels
       integer,                   intent(in) :: energy_formula
       real(kind_phys), optional, intent(in) :: to_dry_factor(:,:)
 
+      character(len=*), parameter :: subname = 'cam_thermo_water_update: '
+
+      if (present(to_dry_factor)) then
+        if (SIZE(to_dry_factor, 1) /= ncol) then
+          call endrun(subname//'DIM 1 of to_dry_factor is'//stringify((/SIZE(to_dry_factor,1)/))//'but should be'//stringify((/ncol/)))
+        end if
+        if (SIZE(to_dry_factor, 2) /= pver) then
+          call endrun(subname//'DIM 2 of to_dry_factor is'//stringify((/SIZE(to_dry_factor,2)/))//'but should be'//stringify((/pver/)))
+        end if
+      end if
+
       call water_composition_update(mmr, ncol, energy_formula, to_dry_factor=to_dry_factor)
     end subroutine cam_thermo_water_update
 

src/data/cam_thermo_formula.F90

@@ -1,5 +1,7 @@
 module cam_thermo_formula
 
+   use runtime_obj, only: unset_int
+
    implicit none
    private
    save
@@ -17,9 +19,9 @@ module cam_thermo_formula
    !! \htmlinclude cam_thermo_formula.html
    ! energy_formula_dycore: energy formula used for dynamical core
    ! written by the dynamical core
-   integer, public :: energy_formula_dycore
+   integer, public :: energy_formula_dycore  = unset_int
    ! energy_formula_physics: energy formula used for physics
-   integer, public :: energy_formula_physics = ENERGY_FORMULA_DYCORE_FV
+   integer, public :: energy_formula_physics = unset_int
 
    ! Public subroutines
    public :: cam_thermo_formula_init

src/data/registry.xml

@@ -264,15 +264,15 @@
     <variable local_name="dycore_energy_consistency_adjust"
               standard_name="flag_for_dycore_energy_consistency_adjustment"
               units="flag" type="logical">
-      <long_name>flag indicating if dynamical core energy not consistent with CAM physics and to perform adjustment of temperature and temperature tendency</long_name>
+      <long_name>flag indicating if dynamical core energy is not consistent with CAM physics and to perform adjustment of temperature and temperature tendency</long_name>
       <initial_value>.false.</initial_value>
       <initial_value dycore="SE,MPAS">.true.</initial_value>
     </variable>
     <!-- Timestep properties -->
     <variable local_name="is_first_timestep"
               standard_name="is_first_timestep"
               units="flag" type="logical">
-      <long_name>flag indicating if is first timestep of initial run</long_name>
+      <long_name>flag indicating if it is the first timestep of an initial run</long_name>
     </variable>
     <!-- State tendencies -->
     <variable local_name="dTdt_total"

src/dynamics/none/dyn_grid.F90

@@ -62,7 +62,6 @@ subroutine model_grid_init()
       use pio,              only: PIO_BCAST_ERROR, pio_seterrorhandling
       use pio,              only: pio_get_var, pio_freedecomp
       use pio,              only: pio_read_darray
-      use pio,              only: pio_inq_att
       use spmd_utils,       only: npes, iam
       use cam_pio_utils,    only: cam_pio_handle_error, cam_pio_find_var
       use cam_pio_utils,    only: cam_pio_var_info, pio_subsystem

src/dynamics/se/dp_coupling.F90

@@ -806,6 +806,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
    call cam_thermo_water_update( &
         mmr            = const_data_ptr,          & ! dry MMR
         ncol           = pcols,                   &
+        pver           = pver,                    &
         energy_formula = ENERGY_FORMULA_DYCORE_SE &
    )
 

src/physics/utils/physics_grid.F90

@@ -378,8 +378,6 @@ end subroutine phys_grid_init
    !========================================================================
 
    real(r8) function get_dlat_p(index)
-      use cam_logfile,    only: iulog
-      use cam_abortutils, only: endrun
       ! latitude of a physics column in degrees
 
       ! Dummy argument
@@ -398,8 +396,6 @@ end function get_dlat_p
    !========================================================================
 
    real(r8) function get_dlon_p(index)
-      use cam_logfile,    only: iulog
-      use cam_abortutils, only: endrun
       ! longitude of a physics column in degrees
 
       ! Dummy argument
@@ -418,8 +414,6 @@ end function get_dlon_p
    !========================================================================
 
    real(r8) function get_rlat_p(index)
-      use cam_logfile,    only: iulog
-      use cam_abortutils, only: endrun
       ! latitude of a physics column in radians
 
       ! Dummy argument
@@ -438,8 +432,6 @@ end function get_rlat_p
    !========================================================================
 
    real(r8) function get_rlon_p(index)
-      use cam_logfile,    only: iulog
-      use cam_abortutils, only: endrun
       ! longitude of a physics column in radians
 
       ! Dummy argument
@@ -458,8 +450,6 @@ end function get_rlon_p
    !========================================================================
 
    real(r8) function get_area_p(index)
-      use cam_logfile,    only: iulog
-      use cam_abortutils, only: endrun
       ! area of a physics column in radians squared
 
       ! Dummy argument
@@ -478,8 +468,6 @@ end function get_area_p
    !========================================================================
 
    real(r8) function get_wght_p(index)
-      use cam_logfile,    only: iulog
-      use cam_abortutils, only: endrun
       ! weight of a physics column in radians squared
 
       ! Dummy argument
@@ -498,8 +486,6 @@ end function get_wght_p
    !========================================================================
 
    subroutine get_rlat_all_p(rlatdim, rlats)
-      use cam_logfile,    only: iulog
-      use cam_abortutils, only: endrun
       !-----------------------------------------------------------------------
       !
       ! get_rlat_all_p: Return all latitudes (in radians) on task.
@@ -528,8 +514,6 @@ end subroutine get_rlat_all_p
    !========================================================================
 
    subroutine get_rlon_all_p(rlondim, rlons)
-      use cam_logfile,    only: iulog
-      use cam_abortutils, only: endrun
       !-----------------------------------------------------------------------
       !
       ! get_rlon_all_p: Return all longitudes (in radians) on task.
@@ -558,8 +542,6 @@ end subroutine get_rlon_all_p
    !========================================================================
 
    subroutine get_wght_all_p(wghtdim, wghts)
-      use cam_logfile,    only: iulog
-      use cam_abortutils, only: endrun
       !-----------------------------------------------------------------------
       !
       ! get_wght_all_p: Return all weights on task.
@@ -588,7 +570,6 @@ end subroutine get_wght_all_p
    !========================================================================
 
    subroutine get_dyn_col_p(index, blk_num, blk_ind)
-      use cam_logfile,    only: iulog
       use cam_abortutils, only: endrun
       ! Return the dynamics local block number and block offset(s) for
       ! the physics column indicated by <index>.
@@ -620,8 +601,6 @@ end subroutine get_dyn_col_p
    !========================================================================
 
    integer function global_index_p(index)
-      use cam_logfile,    only: iulog
-      use cam_abortutils, only: endrun
       ! global column index of a physics column
 
       ! Dummy argument
@@ -638,8 +617,6 @@ integer function global_index_p(index)
    end function global_index_p
 
    integer function local_index_p(index)
-      use cam_logfile,    only: iulog
-      use cam_abortutils, only: endrun
       ! global column index of a physics column
 
       ! Dummy argument

src/utils/gmean_mod.F90

@@ -115,16 +115,16 @@ subroutine gmean_scl (arr, gmean)
       !
       ! Arguments
       !
-      real(r8), intent(in) :: arr(pcols)
-      ! Input array, chunked
-      real(r8), intent(out):: gmean      ! global means
+      real(r8), intent(in)  :: arr(pcols)
+      ! Input array
+      real(r8), intent(out) :: gmean      ! global means
       !
       ! Local workspace
       !
-      integer, parameter :: nflds = 1
-      real(r8)           :: gmean_array(nflds)
-      real(r8)           :: array(pcols, nflds)
-      integer            :: ncols, lchnk
+      integer, parameter    :: nflds = 1
+      real(r8)              :: gmean_array(nflds)
+      real(r8)              :: array(pcols, nflds)
+      integer               :: ncols, lchnk
 
       array(:ncols, 1) = arr(:ncols)
       call gmean_arr(array, gmean_array, nflds)
@@ -192,7 +192,7 @@ end subroutine gmean_float_norepro
    !
    !========================================================================
    !
-   subroutine gmean_fixed_repro (arr, arr_gmean, rel_diff, nflds)
+   subroutine gmean_fixed_repro(arr, arr_gmean, rel_diff, nflds)
       !-----------------------------------------------------------------------
       !
       ! Purpose:
@@ -206,6 +206,7 @@ subroutine gmean_fixed_repro (arr, arr_gmean, rel_diff, nflds)
       use physics_grid,     only: ngcols_p => num_global_phys_cols
       use physconst,        only: pi
       use shr_reprosum_mod, only: shr_reprosum_calc
+      use cam_abortutils,   only: check_allocate
       !
       ! Arguments
       !
@@ -218,14 +219,20 @@ subroutine gmean_fixed_repro (arr, arr_gmean, rel_diff, nflds)
       !
       ! Local workspace
       !
-      integer               :: icol, ifld ! column, field indices
+      real(r8),         parameter :: pi4 = 4.0_r8 * pi
+      character(len=*), parameter :: subname = 'gmean_fixed_repro: '
+
+      integer               :: icol, ifld        ! column, field indices
+      integer               :: errflg
 
       real(r8)              :: wght(pcols)       ! integration weights
       real(r8), allocatable :: xfld(:,:)         ! weighted summands
-      !
-      !-----------------------------------------------------------------------
-      !
-      allocate(xfld(pcols, nflds))
+
+      errflg = 0
+
+      allocate(xfld(pcols, nflds), stat=errflg)
+      call check_allocate(errflg, subname, 'xfld(pcols, nflds)', &
+                          file=__FILE__, line=__LINE__)
 
       ! pre-weight summands
       call get_wght_all_p(pcols, wght)
@@ -249,7 +256,7 @@ subroutine gmean_fixed_repro (arr, arr_gmean, rel_diff, nflds)
 
       deallocate(xfld)
       ! final normalization
-      arr_gmean(:) = arr_gmean(:) / (4.0_r8 * pi)
+      arr_gmean(:) = arr_gmean(:) / pi4
 
    end subroutine gmean_fixed_repro