cam6_4_092: MPAS dp_coupling bug fixes

doc/ChangeLog

@@ -1,3 +1,36 @@
+Tag name: cam6_4_???
+Originator(s): pel
+Date: ??, 2025
+One-line Summary: bugfixes in dp_coupling for MPAS
+Github PR URL: https://github.com/ESCOMP/CAM/pull/
+
+This tag updates the following externals:
+
+Purpose of changes (include the issue number and title text for each relevant GitHub issue):
+   fix issue with mid-level pressure not being bounded by interface pressures in MPAS physics-dynamics coupling: https://github.com/ESCOMP/CAM/issues/1268
+   fix Exner bug in MPAS physics-dynamics coupling: https://github.com/ESCOMP/CAM/issues/753
+
+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: ??
+
+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       cime_config/config_archive.xml
+          - Changes for rpointer files
+
+Could result in answer changes for all MPAS regression tests with full physics and
+possibly with simpler physics, however, this seems to only be a problem at high resolution
+
 ===============================================================
 
 Tag name: cam6_4_090
@@ -172,7 +205,7 @@ derecho/intel/aux_cam:
   ERP_Ln9.f09_f09_mg17.FCSD_HCO.derecho_intel.cam-outfrq9s (Overall: FAIL)
   SMS_Ld1.f09_f09_mg17.FCHIST_GC.derecho_intel.cam-outfrq1d (Overall: DIFF)
   SMS_Lh12.f09_f09_mg17.FCSD_HCO.derecho_intel.cam-outfrq3h (Overall: DIFF)
-      - pre-existing failures due to HEMCO not having reproducible results (issues #1018 and #856)
+	- pre-existing failures due to HEMCO not having reproducible results (issues #1018 and #856)
 
   SMS_D_Ln9_P1280x1.ne0CONUSne30x8_ne0CONUSne30x8_mt12.FCHIST.derecho_intel.cam-outfrq9s (Overall: FAIL)
       - pre-existing failure due to build-namelist error requiring CLM/CTSM external update

src/dynamics/mpas/dp_coupling.F90

@@ -8,7 +8,7 @@ module dp_coupling
 use pmgrid,         only: plev
 use ppgrid,         only: begchunk, endchunk, pcols, pver, pverp
 use constituents,   only: pcnst, cnst_type
-use physconst,      only: gravit, cappa, zvir
+use physconst,      only: gravit, zvir
 use air_composition,only: cpairv
 use air_composition,only: dry_air_species_num
 use dyn_comp,       only: dyn_export_t, dyn_import_t
@@ -417,7 +417,7 @@ subroutine derived_phys(phys_state, phys_tend, pbuf2d)
    use shr_vmath_mod,   only: shr_vmath_log
    use phys_control,    only: waccmx_is
    use cam_thermo,      only: cam_thermo_dry_air_update, cam_thermo_water_update
-   use air_composition, only: rairv, dry_air_species_num
+   use air_composition, only: rairv, dry_air_species_num, cappav
    use qneg_module,     only: qneg3
    use shr_const_mod,   only: shr_const_rwv
    use constituents,    only: qmin
@@ -434,8 +434,6 @@ subroutine derived_phys(phys_state, phys_tend, pbuf2d)
    real(r8) :: factor(pcols,pver)
    real(r8) :: zvirv(pcols,pver)
 
-   real(r8), parameter :: pref = 1.e5_r8 ! reference pressure (Pa)
-
    type(physics_buffer_desc), pointer :: pbuf_chnk(:)
 
    character(len=*), parameter :: subname = 'dp_coupling::derived_phys'
@@ -504,13 +502,6 @@ subroutine derived_phys(phys_state, phys_tend, pbuf2d)
                             phys_state(lchnk)%lnpmid(:ncol,k), ncol)
       end do
 
-      do k = 1, pver
-         phys_state(lchnk)%exner(:ncol,k) = (pref / phys_state(lchnk)%pmid(:ncol,k))**cappa
-      end do
-
-
-
-
       if (dry_air_species_num>0) then
         !------------------------------------------------------------
         ! Apply limiters to mixing ratios of major species
@@ -527,6 +518,10 @@ subroutine derived_phys(phys_state, phys_tend, pbuf2d)
       else
         zvirv(:,:) = zvir
       endif
+      do k = 1, pver
+         phys_state(lchnk)%exner(:ncol,k) = (phys_state(lchnk)%pint(:ncol,pver+1) / phys_state(lchnk)%pmid(:ncol,k))**cappav(:ncol,k,lchnk)
+      end do
+
       !
       ! update cp_dycore in module air_composition.
       ! (note: at this point q is dry)
@@ -836,6 +831,8 @@ subroutine hydrostatic_pressure(nCells, nVertLevels, qsize, index_qv, zz, zgrid,
    real(r8), dimension(nVertLevels+1,nCells) :: pint  ! hydrostatic pressure at interface
    real(r8) :: sum_water
    real(r8) :: pk,rhok,rhodryk,thetavk,kap1,kap2,tvk,tk
+   real(r8), parameter :: epsilon = 0.05_r8
+   real(r8) :: dp_epsilon
    !
    ! For each column, integrate downward from model top to compute dry hydrostatic pressure at layer
    ! midpoints and interfaces. The pressure averaged to layer midpoints should be consistent with
@@ -884,8 +881,15 @@ subroutine hydrostatic_pressure(nCells, nVertLevels, qsize, index_qv, zz, zgrid,
         pintdry(k,iCell) = pintdry(k+1,iCell)+dpdry(k)
         pmid(k,iCell)    = dp(k)   *rgas*tvk/(gravit*dz(k))
         pmiddry(k,iCell) = dpdry(k)*rgas*tk /(gravit*dz(k))
-      end do
-    end do
+        !
+        ! PMID is not necessarily bounded by the hydrostatic interface pressure.
+        ! (has been found to be an issue at ~3.75km resolution in surface layer)
+        !
+        dp_epsilon = dp(k) * epsilon
+        pmid   (k, iCell) = max(min(pmid   (k, iCell), pint   (k, iCell) - dp_epsilon), pint   (k + 1, iCell) + dp_epsilon)
+        pmiddry(k, iCell) = max(min(pmiddry(k, iCell), pintdry(k, iCell) - dp_epsilon), pintdry(k + 1, iCell) + dp_epsilon)
+     end do
+   end do
 end subroutine hydrostatic_pressure
 
 subroutine tot_energy_dyn(nCells, nVertLevels, qsize, index_qv, zz, zgrid, rho_zz, theta_m, q, ux,uy,outfld_name_suffix)