From 9cc084acbafb5a5b7a773b0c49a0f3fdafcdea99 Mon Sep 17 00:00:00 2001
From: Kuan-Chih Wang <kuanchihw@ucar.edu>
Date: Sat, 9 Nov 2024 04:33:42 +0800
Subject: [PATCH] Fix initialization of reference pressure for physics when
 using MPAS dynamical core (#317)

### Tag name (required for release branches):

None

### Originator(s):

PeterHjortLauritzen, nusbaume, kuanchihwang

### Description (include the issue title, and the keyword ['closes',
'fixes', 'resolves'] followed by the issue number):

Presently, MPAS dynamical core is not initializing reference pressure
for physics in the same way as others. This PR backports the fix from
ESCOMP/CAM#1169.

To confirm the fix, observe log entries similar to the following in
`atm.log.<job-id>.<date>-<time>`. The reference pressure at surface now
starts from 1000 hPa. Previously, it starts from 1013.25 hPa.

```
dyn_debug_print (0): Reference layer information:
dyn_debug_print (0): ----- | -------------- | --------------
dyn_debug_print (0): Index |     Height (m) | Pressure (hPa)
dyn_debug_print (0): ----- |   44999.999819 |       1.431394
dyn_debug_print (0):     1 |   43953.595504 |       1.655836
dyn_debug_print (0): ----- |   42907.191189 |       1.880278
dyn_debug_print (0):     2 |   41877.531046 |       2.177235
dyn_debug_print (0): ----- |   40847.870904 |       2.474192
dyn_debug_print (0):     3 |   39835.232205 |       2.867591
dyn_debug_print (0): ----- |   38822.593507 |       3.260989

... (SNIPPED) ...

dyn_debug_print (0): ----- |    1292.109873 |     856.014954
dyn_debug_print (0):    30 |     997.820649 |     887.686997
dyn_debug_print (0): ----- |     703.531425 |     919.359041
dyn_debug_print (0):    31 |     476.271636 |     945.093654
dyn_debug_print (0): ----- |     249.011847 |     970.828267
dyn_debug_print (0):    32 |     124.505924 |     985.414133
dyn_debug_print (0): ----- |       0.000000 |    1000.000000
```

Closes #315.

### Describe any changes made to build system:

None

### Describe any changes made to the namelist:

None

### List any changes to the defaults for the input datasets (e.g.
boundary datasets):

None

### List all files eliminated and why:

None

### List all files added and what they do:

None

### List all existing files that have been modified, and describe the
changes:

* `M       src/dynamics/mpas/dyn_comp.F90`
  * Sort statements
* `M       src/dynamics/mpas/dyn_grid.F90`
* Properly call `std_atm_pres` to initialize reference pressure for
physics
  * Wire up history support
  * More explicit memory management
  * Update code comments
* `M       src/utils/std_atm_profile.F90`
  * Allow custom surface pressure in `std_atm_pres`
  * Remove unused variable in `std_atm_pres`

---------

Co-authored-by: Peter Hjort Lauritzen <pel@ucar.edu>
Co-authored-by: Jesse Nusbaumer <nusbaume@ucar.edu>
---
 src/dynamics/mpas/dyn_comp.F90 |  2 +-
 src/dynamics/mpas/dyn_grid.F90 | 60 ++++++++++++++++++++++++----------
 src/utils/std_atm_profile.F90  | 24 ++++++++++----
 3 files changed, 60 insertions(+), 26 deletions(-)

diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90
index 30e4b2a5..b1892f47 100644
--- a/src/dynamics/mpas/dyn_comp.F90
+++ b/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')
diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90
index 26b0e345..07277b1c 100644
--- a/src/dynamics/mpas/dyn_grid.F90
+++ b/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)
diff --git a/src/utils/std_atm_profile.F90 b/src/utils/std_atm_profile.F90
index f0821d1f..aff962e4 100644
--- a/src/utils/std_atm_profile.F90
+++ b/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
 
 !=========================================================================================