Merge branch 'development' into vdiff-update

.gitmodules

@@ -1,7 +1,7 @@
 [submodule "ccpp-framework"]
         path          = ccpp_framework
         url           = https://github.com/NCAR/ccpp-framework
-        fxtag         = 2024-07-19-dev
+        fxtag         = 2024-10-31-dev
         fxrequired    = AlwaysRequired
 	fxDONOTUSEurl = https://github.com/NCAR/ccpp-framework
 [submodule "history"]

src/control/cam_comp.F90

@@ -86,6 +86,7 @@ subroutine cam_init(caseid, ctitle, model_doi_url,                         &
       use cam_initfiles,        only: cam_initfiles_open
       use dyn_grid,             only: model_grid_init
       use phys_comp,            only: phys_init
+      use phys_comp,            only: phys_register
       use dyn_comp,             only: dyn_init
 !      use cam_restart,          only: cam_read_restart
       use camsrfexch,           only: hub2atm_alloc, atm2hub_alloc
@@ -186,6 +187,10 @@ subroutine cam_init(caseid, ctitle, model_doi_url,                         &
       ! Open initial or restart file, and topo file if specified.
       call cam_initfiles_open()
 
+      ! Call CCPP physics register phase (must happen before
+      !   cam_register_constituents)
+      call phys_register()
+
       ! Initialize constituent information
       !    This will set the total number of constituents and the
       !    number of advected constituents.

src/data/generate_registry_data.py

@@ -1783,16 +1783,7 @@ def gen_registry(registry_file, dycore, outdir, indent,
                                     logger, schema_path=schema_dir,
                                     error_on_noxmllint=error_on_no_validate)
     except CCPPError as ccpperr:
-        cemsg = f"{ccpperr}".split('\n', maxsplit=1)[0]
-        if cemsg[0:12] == 'Execution of':
-            xstart = cemsg.find("'")
-            if xstart >= 0:
-                xend = cemsg[xstart + 1:].find("'") + xstart + 1
-                emsg += '\n' + cemsg[xstart + 1:xend]
-            # end if (else, just keep original message)
-        elif cemsg[0:18] == 'validate_xml_file:':
-            emsg += "\n" + cemsg
-        # end if
+        emsg += f"\n{ccpperr}"
         file_ok = False
     # end try
     if not file_ok:

src/data/registry.xml

@@ -154,7 +154,7 @@
               allocatable="pointer">
       <long_name>inverse exner function w.r.t. surface pressure (ps/p)^(R/cp)</long_name>
       <dimensions>horizontal_dimension vertical_layer_dimension</dimensions>
-      <ic_file_input_names>exner state_exner</ic_file_input_names>
+      <ic_file_input_names>exner state_exner inverse_exner_function_wrt_surface_pressure</ic_file_input_names>
     </variable>
     <variable local_name="zm"
               standard_name="geopotential_height_wrt_surface"

src/dynamics/mpas/dyn_comp.F90

@@ -812,7 +812,6 @@ subroutine mark_variable_as_initialized()
         call mark_as_initialized('eastward_wind')
         call mark_as_initialized('geopotential_height_wrt_surface')
         call mark_as_initialized('geopotential_height_wrt_surface_at_interface')
-        call mark_as_initialized('reciprocal_of_dimensionless_exner_function_wrt_surface_air_pressure')
         call mark_as_initialized('lagrangian_tendency_of_air_pressure')
         call mark_as_initialized('ln_air_pressure')
         call mark_as_initialized('ln_air_pressure_at_interface')
@@ -821,6 +820,7 @@ subroutine mark_variable_as_initialized()
         call mark_as_initialized('northward_wind')
         call mark_as_initialized('reciprocal_of_air_pressure_thickness')
         call mark_as_initialized('reciprocal_of_air_pressure_thickness_of_dry_air')
+        call mark_as_initialized('reciprocal_of_dimensionless_exner_function_wrt_surface_air_pressure')
         call mark_as_initialized('surface_air_pressure')
         call mark_as_initialized('surface_geopotential')
         call mark_as_initialized('surface_pressure_of_dry_air')

src/dynamics/mpas/dyn_grid.F90

@@ -5,13 +5,14 @@ module dyn_grid
     use cam_grid_support, only: cam_grid_register, cam_grid_attribute_register, &
                                 horiz_coord_t, horiz_coord_create, &
                                 max_hcoordname_len
+    use cam_history_support, only: add_vert_coord
     use cam_initfiles, only: initial_file_get_id
     use cam_map_utils, only: kind_imap => imap
     use dyn_comp, only: dyn_debug_print, mpas_dynamical_core, &
         ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels, &
         ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max, &
         sphere_radius
-    use dynconst, only: constant_pi => pi, rad_to_deg, dynconst_init
+    use dynconst, only: constant_p0 => pref, constant_pi => pi, rad_to_deg, dynconst_init
     use physics_column_type, only: kind_pcol, physics_column_t
     use physics_grid, only: phys_decomp, phys_grid_init
     use ref_pres, only: ref_pres_init
@@ -107,7 +108,7 @@ subroutine model_grid_init()
             call endrun('Numbers of vertical layers mismatch', subname, __LINE__)
         end if
 
-        ! Initialize reference pressure.
+        ! Initialize reference pressure for use by physics.
         call dyn_debug_print('Calling init_reference_pressure')
 
         call init_reference_pressure()
@@ -123,7 +124,7 @@ subroutine model_grid_init()
         call define_cam_grid()
     end subroutine model_grid_init
 
-    !> Initialize reference pressure by computing necessary variables and calling `ref_pres_init`.
+    !> Initialize reference pressure for use by physics.
     !> (KCW, 2024-03-25)
     subroutine init_reference_pressure()
         character(*), parameter :: subname = 'dyn_grid::init_reference_pressure'
@@ -139,10 +140,16 @@ subroutine init_reference_pressure()
         ! `zw` denotes zeta at w-wind levels (i.e., at layer interfaces).
         ! `dzw` denotes the delta/difference between `zw`.
         ! `rdzw` denotes the reciprocal of `dzw`.
-        real(kind_r8), allocatable :: zu(:), zw(:), dzw(:)
+        real(kind_r8), allocatable :: dzw(:)
         real(kind_r8), pointer :: rdzw(:)
+        real(kind_r8), pointer :: zu(:) ! CANNOT be safely deallocated because `add_vert_coord`
+                                        ! just uses pointers to point at it internally.
+        real(kind_r8), pointer :: zw(:) ! CANNOT be safely deallocated because `add_vert_coord`
+                                        ! just uses pointers to point at it internally.
 
         nullify(rdzw)
+        nullify(zu)
+        nullify(zw)
 
         ! Compute reference height.
         call mpas_dynamical_core % get_variable_pointer(rdzw, 'mesh', 'rdzw')
@@ -169,11 +176,19 @@ subroutine init_reference_pressure()
             zu(k) = 0.5_kind_r8 * (zw(k + 1) + zw(k))
         end do
 
+        deallocate(dzw)
+
+        ! Register zeta coordinates with history.
+        call add_vert_coord('ilev', pverp, 'Height (zeta) level at layer interfaces', 'm', zw, &
+            positive='up')
+        call add_vert_coord('lev', pver, 'Height (zeta) level at layer midpoints', 'm', zu, &
+            positive='up')
+
         ! Compute reference pressure from reference height.
         allocate(p_ref_int(pverp), stat=ierr)
         call check_allocate(ierr, subname, 'p_ref_int(pverp)', 'dyn_grid', __LINE__)
 
-        call std_atm_pres(zw, p_ref_int)
+        call std_atm_pres(zw, p_ref_int, user_specified_ps=constant_p0)
 
         allocate(p_ref_mid(pver), stat=ierr)
         call check_allocate(ierr, subname, 'p_ref_mid(pver)', 'dyn_grid', __LINE__)
@@ -201,6 +216,12 @@ subroutine init_reference_pressure()
             ' |  ' // stringify([p_ref_int(k) / 100.0_kind_r8]))
 
         call ref_pres_init(p_ref_int, p_ref_mid, num_pure_p_lev)
+
+        deallocate(p_ref_int)
+        deallocate(p_ref_mid)
+
+        nullify(zu)
+        nullify(zw)
     end subroutine init_reference_pressure
 
     !> Initialize physics grid in terms of dynamics decomposition.
@@ -209,7 +230,7 @@ end subroutine init_reference_pressure
     subroutine init_physics_grid()
         character(*), parameter :: subname = 'dyn_grid::init_physics_grid'
         character(max_hcoordname_len), allocatable :: dyn_attribute_name(:)
-        integer :: hdim1_d, hdim2_d
+        integer :: hdim1_d, hdim2_d ! First and second horizontal dimensions of physics grid.
         integer :: i
         integer :: ierr
         integer, pointer :: indextocellid(:)  ! Global indexes of cell centers.
@@ -274,6 +295,9 @@ subroutine init_physics_grid()
         call check_allocate(ierr, subname, 'dyn_attribute_name(0)', 'dyn_grid', __LINE__)
 
         call phys_grid_init(hdim1_d, hdim2_d, 'mpas', dyn_column, 'mpas_cell', dyn_attribute_name)
+
+        deallocate(dyn_column)
+        deallocate(dyn_attribute_name)
     end subroutine init_physics_grid
 
     !> This subroutine defines and registers four variants of dynamics grids in terms of dynamics decomposition.
@@ -298,25 +322,25 @@ subroutine define_cam_grid()
         real(kind_r8), pointer :: lonedge(:)   ! Edge node longitudes (radians).
         real(kind_r8), pointer :: lonvertex(:) ! Vertex node longitudes (radians).
 
-        ! Global grid indexes. CAN be safely deallocated because its values are copied into
-        ! `cam_grid_attribute_*_t` and `horiz_coord_t`.
+        ! Global grid indexes. CAN be safely deallocated because its values are copied internally by
+        ! `cam_grid_attribute_register` and `horiz_coord_create`.
         ! `kind_imap` is an integer kind of `PIO_OFFSET_KIND`.
         integer(kind_imap), pointer :: global_grid_index(:)
-        ! Global grid maps. CANNOT be safely deallocated because `cam_filemap_t`
-        ! just uses pointers to point at it.
+        ! Global grid maps. CANNOT be safely deallocated because `cam_grid_register`
+        ! just uses pointers to point at it internally.
         ! `kind_imap` is an integer kind of `PIO_OFFSET_KIND`.
         integer(kind_imap), pointer :: global_grid_map(:, :)
-        ! Cell areas (square meters). CANNOT be safely deallocated because `cam_grid_attribute_*_t`
-        ! just uses pointers to point at it.
+        ! Cell areas (square meters). CANNOT be safely deallocated because `cam_grid_attribute_register`
+        ! just uses pointers to point at it internally.
         real(kind_r8), pointer :: cell_area(:)
-        ! Cell weights normalized to unity. CANNOT be safely deallocated because `cam_grid_attribute_*_t`
-        ! just uses pointers to point at it.
+        ! Cell weights normalized to unity. CANNOT be safely deallocated because `cam_grid_attribute_register`
+        ! just uses pointers to point at it internally.
         real(kind_r8), pointer :: cell_weight(:)
-        ! Latitude coordinates. CANNOT be safely deallocated because `cam_grid_t`
-        ! just uses pointers to point at it.
+        ! Latitude coordinates. CANNOT be safely deallocated because `cam_grid_register`
+        ! just uses pointers to point at it internally.
         type(horiz_coord_t), pointer :: lat_coord
-        ! Longitude coordinates. CANNOT be safely deallocated because `cam_grid_t`
-        ! just uses pointers to point at it.
+        ! Longitude coordinates. CANNOT be safely deallocated because `cam_grid_register`
+        ! just uses pointers to point at it internally.
         type(horiz_coord_t), pointer :: lon_coord
 
         nullify(indextocellid, indextoedgeid, indextovertexid)

src/physics/utils/phys_comp.F90

@@ -10,6 +10,7 @@ module phys_comp
    private
 
    public :: phys_readnl
+   public :: phys_register
    public :: phys_init
    public :: phys_timestep_init
    public :: phys_run1
@@ -129,26 +130,14 @@ subroutine phys_readnl(nlfilename)
 
    end subroutine phys_readnl
 
-   subroutine phys_init()
-      use cam_abortutils,       only: endrun
-      use physics_grid,         only: columns_on_task
-      use vert_coord,           only: pver, pverp
-      use cam_thermo,           only: cam_thermo_init
-      use physics_types,        only: allocate_physics_types_fields
-      use cam_ccpp_cap,         only: cam_ccpp_physics_initialize
+   subroutine phys_register()
+      use cam_ccpp_cap,         only: cam_ccpp_physics_register
       use cam_ccpp_cap,         only: ccpp_physics_suite_part_list
+      use cam_abortutils,       only: endrun
 
       ! Local variables
       integer                    :: i_group
 
-      call cam_thermo_init(columns_on_task, pver, pverp)
-
-      call allocate_physics_types_fields(columns_on_task, pver, pverp,        &
-           set_init_val_in=.true., reallocate_in=.false.)
-      call cam_ccpp_physics_initialize(phys_suite_name)
-      if (errcode /= 0) then
-         call endrun('cam_ccpp_physics_initialize: '//trim(errmsg))
-      end if
       call ccpp_physics_suite_part_list(phys_suite_name, suite_parts,         &
            errmsg, errcode)
       if (errcode /= 0) then
@@ -158,11 +147,35 @@ subroutine phys_init()
       ! Confirm that the suite parts are as expected
       do i_group = 1, size(suite_parts)
          if (.not. any(suite_parts(i_group) == suite_parts_expect)) then
-            write(errmsg, *) 'phys_init: SDF suite groups MUST be ',             &
+            write(errmsg, *) 'phys_register: SDF suite groups MUST be ',      &
                 'ONLY `physics_before_coupler` and/or `physics_after_coupler`'
             call endrun(errmsg)
          end if
       end do
+      ! Call CCPP register phase
+      call cam_ccpp_physics_register(phys_suite_name)
+      if (errcode /= 0) then
+         call endrun('cam_ccpp_physics_register: '//trim(errmsg))
+      end if
+
+   end subroutine phys_register
+
+   subroutine phys_init()
+      use cam_abortutils,       only: endrun
+      use physics_grid,         only: columns_on_task
+      use vert_coord,           only: pver, pverp
+      use cam_thermo,           only: cam_thermo_init
+      use physics_types,        only: allocate_physics_types_fields
+      use cam_ccpp_cap,         only: cam_ccpp_physics_initialize
+
+      call cam_thermo_init(columns_on_task, pver, pverp)
+
+      call allocate_physics_types_fields(columns_on_task, pver, pverp,        &
+           set_init_val_in=.true., reallocate_in=.false.)
+      call cam_ccpp_physics_initialize(phys_suite_name)
+      if (errcode /= 0) then
+         call endrun('cam_ccpp_physics_initialize: '//trim(errmsg))
+      end if
 
    end subroutine phys_init
 

src/utils/std_atm_profile.F90

@@ -52,16 +52,27 @@ module std_atm_profile
 CONTAINS
 !=========================================================================================
 
-subroutine std_atm_pres(height, pstd)
+subroutine std_atm_pres(height, pstd, user_specified_ps)
 
    ! arguments
-   real(r8), intent(in)  :: height(:) ! height above sea level in meters
-   real(r8), intent(out) :: pstd(:)   ! std pressure in Pa
+   real(r8),           intent(in)  :: height(:) ! height above sea level in meters
+   real(r8),           intent(out) :: pstd(:)   ! std pressure in Pa
+   real(r8), optional, intent(in)  :: user_specified_ps
+
+   integer  :: i, ii, k, nlev
+   real(r8) :: pb_local(nreg)
 
-   integer :: i, ii, k, nlev
    character(len=*), parameter :: routine = 'std_atm_pres'
    !----------------------------------------------------------------------------
 
+   ! Initialize local standard pressure values array
+   pb_local = pb
+
+   ! Set new surface pressure value if provided by the caller
+   if (present(user_specified_ps)) then
+      pb_local(1) = user_specified_ps
+   end if
+
    nlev = size(height)
    do k = 1, nlev
       if (height(k) < 0.0_r8) then
@@ -78,13 +89,12 @@ subroutine std_atm_pres(height, pstd)
       end if
 
       if (lb(ii) /= 0._r8) then
-         pstd(k) = pb(ii) * ( tb(ii) / (tb(ii) + lb(ii)*(height(k) - hb(ii)) ) )**(c1/lb(ii))
+         pstd(k) = pb_local(ii) * ( tb(ii) / (tb(ii) + lb(ii)*(height(k) - hb(ii)) ) )**(c1/lb(ii))
       else
-         pstd(k) = pb(ii) * exp( -c1*(height(k) - hb(ii))/tb(ii) )
+         pstd(k) = pb_local(ii) * exp( -c1*(height(k) - hb(ii))/tb(ii) )
       end if
 
    end do
-
 end subroutine std_atm_pres
 
 !=========================================================================================