cam6_4_042: Fix reference pressures for MPAS

doc/ChangeLog

@@ -1,5 +1,51 @@
 ===============================================================
 
+Tag name: cam6_4_xxx
+Originator(s): pel, nusbaume
+Date: Oct 9, 2024
+One-line Summary: Fix reference pressures in MPAS
+Github PR URL: https://github.com/ESCOMP/CAM/pull/1169
+
+Purpose of changes (include the issue number and title text for each relevant GitHub issue):
+- https://github.com/ESCOMP/CAM/issues/1168
+
+Describe any changes made to build system: N/A
+
+Describe any changes made to the namelist: N/A
+
+List any changes to the defaults for the boundary datasets: N/A
+
+Describe any substantial timing or memory changes: N/A
+
+Code reviewed by: nusbaume, cacraigucar
+
+List all files eliminated: N/A
+
+List all files added and what they do: N/A
+
+List all existing files that have been modified, and describe the changes:
+
+M       src/dynamics/mpas/dyn_grid.F90
+M       src/utils/std_atm_profile.F90
+ - Use CAM reference pressure instead of standard atmosphere surface pressure.
+
+If there were any failures reported from running test_driver.sh on any test
+platform, and checkin with these failures has been OK'd by the gatekeeper,
+then copy the lines from the td.*.status files for the failed tests to the
+appropriate machine below.  All failed tests must be justified.
+
+derecho/intel/aux_cam:
+
+derecho/nvhpc/aux_cam:
+
+izumi/nag/aux_cam:
+
+izumi/gnu/aux_cam:
+
+Summarize any changes to answers: all MPAS tests will likely fail against baselines
+
+===============================================================
+
 Tag name: cam6_4_041
 Originator(s): fvitt
 Date: 16 Oct 2024

src/dynamics/mpas/dyn_grid.F90

@@ -105,6 +105,7 @@ subroutine dyn_grid_init()
    ! coupling code requires constructing global fields for the cell center
    ! grid which is used by the physics parameterizations.
 
+   use hycoef,              only: ps0
    use ref_pres,            only: ref_pres_init
    use std_atm_profile,     only: std_atm_pres
    use time_manager,        only: get_step_size
@@ -141,7 +142,7 @@ subroutine dyn_grid_init()
    ierr = pio_get_att(fh_ini, pio_global, 'sphere_radius', sphere_radius)
 
    ! Compute reference pressures from reference heights.
-   call std_atm_pres(zw, pref_edge)
+   call std_atm_pres(zw, pref_edge, user_specified_ps=ps0)
    pref_mid = (pref_edge(1:plev) + pref_edge(2:plevp)) * 0.5_r8
 
    num_pr_lev = 0
@@ -173,7 +174,7 @@ subroutine dyn_grid_init()
    ! centered grid is used by the physics parameterizations.  The physics
    ! decomposition of the cell centered grid is defined in phys_grid_init.
    call define_cam_grids()
-   
+
 end subroutine dyn_grid_init
 
 !=========================================================================================
@@ -380,7 +381,7 @@ integer function get_dyn_grid_parm(name) result(ival)
       ival = nCells_g
    else if(name == 'plev') then
       ival = plev
-   else	
+   else
       ival = -1
    end if
 
@@ -428,7 +429,7 @@ subroutine dyn_grid_get_elem_coords(ie, rlon, rlat, cdex )
    ! for a field on a dynamics decomposition.  The code in hbuf_accum_addlcltime appears
    ! to also assume that the field is on the physics grid since there is no argument
    ! passed to specify which dynamics grid the coordinates are for.
-   
+
    call endrun(subname//': not implemented for the MPAS grids')
 
 end subroutine dyn_grid_get_elem_coords
@@ -545,7 +546,7 @@ subroutine define_cam_grids()
    use cam_grid_support, only: horiz_coord_t, horiz_coord_create, iMap
    use cam_grid_support, only: cam_grid_register, cam_grid_attribute_register
    use shr_const_mod,    only: PI => SHR_CONST_PI
- 
+
    ! Local variables
    integer :: i, j
 
@@ -600,7 +601,7 @@ subroutine define_cam_grids()
           'degrees_north', 1, nCellsSolve, latCell(1:nCellsSolve)*rad2deg, map=gidx)
    lon_coord => horiz_coord_create('lonCell', 'nCells', nCells_g, 'longitude',     &
           'degrees_east', 1, nCellsSolve, lonCell(1:nCellsSolve)*rad2deg, map=gidx)
- 
+
    allocate(areaWeight(nCellsSolve), stat=ierr)
    if( ierr /= 0 ) call endrun(subname//':failed to allocate area_weight :'//int2str(__LINE__))
    call get_dyn_grid_info(hdim1_d, hdim2_d, num_levels, index_model_top_layer, index_surface_layer, unstructured, dyn_cols)
@@ -667,7 +668,7 @@ subroutine define_cam_grids()
           'degrees_north', 1, nEdgesSolve, latEdge(1:nEdgesSolve)*rad2deg, map=gidx)
    lon_coord => horiz_coord_create('lonEdge', 'nEdges', nEdges_g, 'longitude',     &
           'degrees_east', 1, nEdgesSolve, lonEdge(1:nEdgesSolve)*rad2deg, map=gidx)
- 
+
    ! Map for edge node grid
    allocate(grid_map(3, nEdgesSolve), stat=ierr)
    if( ierr /= 0 ) call endrun(subname//':failed to allocate grid_map array at line:'//int2str(__LINE__))
@@ -704,7 +705,7 @@ subroutine define_cam_grids()
           'degrees_north', 1, nVerticesSolve, latVertex(1:nVerticesSolve)*rad2deg, map=gidx)
    lon_coord => horiz_coord_create('lonVertex', 'nVertices', nVertices_g, 'longitude',     &
           'degrees_east', 1, nVerticesSolve, lonVertex(1:nVerticesSolve)*rad2deg, map=gidx)
- 
+
    ! Map for vertex node grid
    allocate(grid_map(3, nVerticesSolve), stat=ierr)
    if( ierr /= 0 ) call endrun(subname//':failed to allocate grid_map array at line:'//int2str(__LINE__))
@@ -723,7 +724,7 @@ subroutine define_cam_grids()
    nullify(grid_map)
    nullify(lat_coord)
    nullify(lon_coord)
-   
+
 end subroutine define_cam_grids
 
 end module dyn_grid

src/utils/std_atm_profile.F90

@@ -47,21 +47,33 @@ module std_atm_profile
 real(r8), parameter :: g0 = 9.80665_r8   ! gravitational acceleration (m/s^2)
 real(r8), parameter :: mw = 0.0289644_r8 ! molar mass of dry air (kg/mol)
 real(r8), parameter :: c1 = g0*mw/rg
-  
+
 !=========================================================================================
 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
-
-   integer :: i, ii, k, nlev
+   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
+   integer  :: ierr
+   real(r8) :: pb_local(nreg)
+
    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
@@ -76,33 +88,32 @@ subroutine std_atm_pres(height, pstd)
             end if
          end do find_region
       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 do
 end subroutine std_atm_pres
 
 !=========================================================================================
 
 subroutine std_atm_height(pstd, height)
-    
+
    ! arguments
    real(r8), intent(in)   :: pstd(:)   ! std pressure in Pa
    real(r8), intent(out)  :: height(:) ! height above sea level in meters
-    
+
    integer :: i, ii, k, nlev
    logical :: found_region
    character(len=*), parameter :: routine = 'std_atm_height'
    !----------------------------------------------------------------------------
-    
+
    nlev = size(height)
    do k = 1, nlev
-      
+
       if (pstd(k) <= pb(nreg)) then
          ii = nreg
       else if (pstd(k) > pb(1)) then
@@ -129,16 +140,16 @@ end subroutine std_atm_height
 !=========================================================================================
 
 subroutine std_atm_temp(height, temp)
-    
+
    ! arguments
    real(r8), intent(in)   :: height(:) ! std pressure in Pa
    real(r8), intent(out)  :: temp(:)   ! temperature
-    
+
    ! local vars
    integer :: i, ii, k, nlev
    character(len=*), parameter :: routine = 'std_atm_temp'
    !----------------------------------------------------------------------------
-    
+
    nlev = size(height)
    do k = 1, nlev
       if (height(k) < 0.0_r8) then
@@ -158,7 +169,7 @@ subroutine std_atm_temp(height, temp)
       else
          temp(k) = tb(ii)
       end if
-      
+
    end do
 
 end subroutine std_atm_temp