Merge pull request #1304 from PeterHjortLauritzen/MPAS-pmid

cam6_4_092: MPAS dp_coupling bug fixes

Author: cacraigucar

doc/ChangeLog

@@ -1,3 +1,57 @@
+Tag name: cam6_4_092
+Originator(s): pel, cacraig
+Date: May 14, 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: nusbaume, kuanchihwang
+
+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/dp_coupling.F90
+          - Changes as described above
+
+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:
+  SMS_D_Ln9_P1280x1.ne0CONUSne30x8_ne0CONUSne30x8_mt12.FCHIST.derecho_intel.cam-outfrq9s (Overall: FAIL) details:
+    FAIL SMS_D_Ln9_P1280x1.ne0CONUSne30x8_ne0CONUSne30x8_mt12.FCHIST.derecho_intel.cam-outfrq9s SETUP
+      - Pre-existing failures due to build-namelist error requiring CLM/CTSM external update
+
+  ERS_D_Ln9.mpasa480_mpasa480.F2000climo.derecho_intel.cam-outfrq9s_mpasa480_pertlim (Overall: DIFF) details:
+  ERS_D_Ln9.mpasa480_mpasa480.FHISTC_LTso.derecho_intel.cam-outfrq9s_mpasa480 (Overall: DIFF) details:
+  ERS_Ln9_P288x1.mpasa120_mpasa120.F2000climo.derecho_intel.cam-outfrq9s_mpasa120 (Overall: DIFF) details:
+     - Answer changes due to MPAS changes
+
+derecho/nvhpc/aux_cam: all BFB
+
+izumi/nag/aux_cam:  all BFB
+
+izumi/gnu/aux_cam:
+  ERS_D_Ln9.mpasa480_mpasa480.FHISTC_LTso.izumi_gnu.cam-outfrq9s_mpasa480 (Overall: DIFF) details:
+  ERS_Ln9_P24x1.mpasa480_mpasa480.F2000climo.izumi_gnu.cam-outfrq9s_mpasa480 (Overall: DIFF) details:
+     - Answer changes due to MPAS changes
+
 ===============================================================
 
 Tag name: cam6_4_091
@@ -346,7 +400,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
@@ -458,7 +458,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
@@ -475,8 +475,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'
@@ -545,13 +543,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
@@ -568,6 +559,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,pverp) / 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)
@@ -877,6 +872,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, dpdry_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
@@ -925,8 +922,16 @@ 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
+        dpdry_epsilon = dpdry(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) - dpdry_epsilon), pintdry(k + 1, iCell) + dpdry_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)