Merge pull request #1467 from PeterHjortLauritzen/frontogenesis-fv-fix

cam6_4_150: fix FV frontogenesis bug and diagnostic PS/PSDRY bug in SE log file

Author: cacraigucar

doc/ChangeLog

@@ -1,3 +1,82 @@
+
+===============================================================
+
+Tag name: cam6_4_150
+Originator(s): Jack Chen, Peter Lauritzen, cacraigucar
+Date: Feb 9, 2026
+One-line Summary: fix frontogenesis function bug for FV
+Github PR URL: https://github.com/ESCOMP/CAM/pull/1467
+
+Purpose of changes (include the issue number and title text for each relevant GitHub issue):
+- Gradient computations for the fronotogenesis function were originally computed over terrain-following coordinate.
+  Correction terms are added so gradients are on isobaric coordinates.
+- diagnostic PS and PSDRY computation in CAM-SE dycore written to log file is missing mass over model top.
+
+
+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: 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/fv/gravity_waves_sources.F90
+   M       src/dynamics/se/dycore/prim_state_mod.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:
+
+  FAIL ERI_D_Ln18.ne16pg3_ne16pg3_mt232.FHIST_C4.derecho_intel.cam-outfrq3s_eri
+     - pre-existing failure
+     - ERI test discovered a bug in CICE:  See: https://github.com/ESCOMP/CESM_CICE/issues/34
+
+  ERP_Lh12.f19_f19_mg17.FW4madSD.derecho_intel.cam-outfrq3h (Overall: DIFF) details:
+  ERP_Ln9.f09_f09_mg17.FCSD_HCO.derecho_intel.cam-outfrq9s (Overall: DIFF) details:
+  ERP_Ln9.f19_f19_mg17.FWsc1850.derecho_intel.cam-outfrq9s (Overall: DIFF) details:
+  ERS_Ln9.f09_f09_mg17.FX2000.derecho_intel.cam-outfrq9s (Overall: DIFF) details:
+  ERS_Ln9.f19_f19_mg17.FXSD.derecho_intel.cam-outfrq9s (Overall: DIFF) details:
+  SMS_D_Ln9.f09_f09_mg17.FCts2nudged.derecho_intel.cam-outfrq9s_leapday (Overall: DIFF) details:
+  SMS_D_Ln9.f09_f09_mg17.FCvbsxHIST.derecho_intel.cam-outfrq9s (Overall: DIFF) details:
+  SMS_D_Ln9.f19_f19_mg17.FCARMA2000climo.derecho_intel.cam-outfrq9s (Overall: DIFF) details:
+  SMS_D_Ln9.f19_f19_mg17.FWma2000climo.derecho_intel.cam-outfrq9s (Overall: DIFF) details:
+  SMS_D_Ln9.f19_f19_mg17.FWma2000climo.derecho_intel.cam-outfrq9s_waccm_ma_mam4 (Overall: DIFF) details:
+  SMS_D_Ln9.f19_f19_mg17.QPC2000climo.derecho_intel.cam-outfrq3s_usecase (Overall: DIFF) details:
+  SMS_Ld1.f09_f09_mg17.FCHIST_GC.derecho_intel.cam-outfrq1d (Overall: DIFF) details:
+  SMS_Ld1.f09_f09_mg17.FW2000climo.derecho_intel.cam-outfrq1d (Overall: DIFF) details:
+  SMS_Lh12.f09_f09_mg17.FCSD_HCO.derecho_intel.cam-outfrq3h (Overall: DIFF) details:
+  SMS_Ln9.f09_f09_mg17.FW1850.derecho_intel.cam-reduced_hist3s (Overall: DIFF) details:
+      - Answer changes for FV expected
+
+derecho/nvhpc/aux_cam: all pass
+
+izumi/nag/aux_cam:
+  ERC_D_Ln9.f10_f10_mg37.QPWmaC6.izumi_nag.cam-outfrq3s (Overall: DIFF) details:
+     - Answer changes for FV expected
+
+izumi/gnu/aux_cam:
+  ERC_D_Ln9.f10_f10_mg37.QPC4.izumi_gnu.cam-outfrq3s_diags (Overall: DIFF) details:
+  SMS_D_Ln9.f10_f10_mg37.FWmaHIST.izumi_gnu.cam-outfrq9s_mee_fluxes (Overall: DIFF) details:
+  SMS_D_Ln9.f10_f10_mg37.QPWmaC4.izumi_gnu.cam-outfrq9s_apmee (Overall: DIFF) details:
+  SMS_D_Ln9.f19_f19_mt232.FW4madSD.izumi_gnu.cam-outfrq9s (Overall: DIFF) details:
+     - Answer changes for FV expected
+
+Summarize any changes to answers:
+  FV answers differences are climate changing - ACOM will need to retune at a future date
+
+===============================================================
 ===============================================================
 
 Tag name: cam6_4_149

src/dynamics/fv/gravity_waves_sources.F90

@@ -58,6 +58,11 @@ subroutine gws_src_fnct(grid, u3,v3,pt, q3, pe, frontgf, frontga)
       real(r8) ::  ug(grid%ifirstxy-1:grid%ilastxy+1,plev,grid%jfirstxy-1:grid%jlastxy+1)  ! U-wind ghosted
       real(r8) ::  vg(grid%ifirstxy-1:grid%ilastxy+1,plev,grid%jfirstxy-1:grid%jlastxy+1)  ! V-wind ghosted
 
+      real(r8) :: pmg(grid%ifirstxy-1:grid%ilastxy+1,plev,grid%jfirstxy-1:grid%jlastxy+1)  ! mid-point pressure ghosted
+      real(r8) :: dptdp(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy)
+      real(r8) :: dudp(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy)
+      real(r8) :: dvdp(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy)
+
       real(r8) :: tglat                               ! tangent-latitude
       integer  :: i,j,k
       integer  :: im, ip
@@ -100,6 +105,12 @@ subroutine gws_src_fnct(grid, u3,v3,pt, q3, pe, frontgf, frontga)
       call ghost_array(grid, u3, ug)
       call ghost_array(grid, v3, vg)
 
+      call ghost_array(grid, pm, pmg)
+
+      call compute_vertical_derivative(grid,pe,pm,ptemp,dptdp)
+      call compute_vertical_derivative(grid,pe,pm,u3,dudp)
+      call compute_vertical_derivative(grid,pe,pm,v3,dvdp)
+
       !$omp parallel do private (i,j,k)
       do k=1, plev
          do i=beglonxy, endlonxy
@@ -108,15 +119,32 @@ subroutine gws_src_fnct(grid, u3,v3,pt, q3, pe, frontgf, frontga)
                ! Pot. Temperature
                pty(i,k,j) = ( ptg(i,k,j+1) - ptg(i,k,j-1) ) / (2._r8 * grid%dp)
                pty(i,k,j) = pty(i,k,j) / aearth
+               !
+               ! Topography/isobaric correction term:
+               ! The horizontal temperature gradient for the frontogenesis function is initially
+               ! computed in terrain-following coordinates. Correction terms are then applied to
+               ! convert it to an isobaric (pressure) surface reference. The same correction is
+               ! performed for the wind component terms.
+               !
+               pty(i,k,j) = pty(i,k,j) - dptdp(i,k,j) * &
+                            ( pmg(i,k,j+1) - pmg(i,k,j-1) ) / (2._r8 * grid%dp) / aearth
 
                ! U-wind
                uy(i,k,j) = ( ug(i,k,j+1) - ug(i,k,j-1) ) / (2._r8 * grid%dp)
                uy(i,k,j) = uy(i,k,j) / aearth
 
+               ! Topography correction term
+               uy(i,k,j) = uy(i,k,j) - dudp(i,k,j) * &
+                           ( pmg(i,k,j+1) - pmg(i,k,j-1) ) / (2._r8 * grid%dp) / aearth
+
                ! V-wind
                vy(i,k,j) = ( vg(i,k,j+1) - vg(i,k,j-1) ) / (2._r8 * grid%dp)
                vy(i,k,j) = vy(i,k,j) / aearth
 
+               ! Topography correction term
+               vy(i,k,j) = vy(i,k,j) - dvdp(i,k,j) * &
+                           ( pmg(i,k,j+1) - pmg(i,k,j-1) ) / (2._r8 * grid%dp) / aearth
+
             end do
          end do
       end do
@@ -135,14 +163,26 @@ subroutine gws_src_fnct(grid, u3,v3,pt, q3, pe, frontgf, frontga)
                ptx(i,k,j) = ( ptg(ip,k,j) - ptg(im,k,j) ) / (2._r8 * grid%dl)
                ptx(i,k,j) = ptx(i,k,j) / (aearth * (grid%cosp(j)+1.e-3_r8))
 
+               ! Topography correction
+               ptx(i,k,j) = ptx(i,k,j) - dptdp(i,k,j) * &
+                            ( pmg(ip,k,j) - pmg(im,k,j) )/(2._r8 * grid%dl)/(aearth * (grid%cosp(j)+1.e-3_r8))
+
                ! U-wind
                ux(i,k,j) = ( ug(ip,k,j) - ug(im,k,j) ) / (2._r8 *grid%dl)
+               ! Topography correction
                ux(i,k,j) = ux(i,k,j) / (aearth * (grid%cosp(j)+1.e-3_r8))
 
+               ! Topography correction
+               ux(i,k,j) = ux(i,k,j) - dudp(i,k,j) * &
+                           ( pmg(ip,k,j) - pmg(im,k,j) )/(2._r8 * grid%dl)/(aearth * (grid%cosp(j)+1.e-3_r8))
+
                ! V-wind
                vx(i,k,j) = ( vg(ip,k,j) - vg(im,k,j) ) / (2._r8 *grid%dl)
                vx(i,k,j) = vx(i,k,j) / (aearth * (grid%cosp(j)+1.e-3_r8))
 
+               ! Topography correction
+               vx(i,k,j) = vx(i,k,j) - dvdp(i,k,j) * &
+                           ( pmg(ip,k,j) - pmg(im,k,j) )/(2._r8 * grid%dl)/(aearth * (grid%cosp(j)+1.e-3_r8))
             end do
          end do
       end do
@@ -180,6 +220,66 @@ subroutine gws_src_fnct(grid, u3,v3,pt, q3, pe, frontgf, frontga)
 
     end subroutine gws_src_fnct
 
+  subroutine compute_vertical_derivative(grid,pint,pmid,data,ddata_dp)
+    !---------------------------------------------------------------------------
+    use dynamics_vars, only: t_fvdycore_grid
+
+    type (t_fvdycore_grid), intent(in) :: grid    ! grid for XY decomp
+    real(r8),   intent(in ) :: pint(grid%ifirstxy:grid%ilastxy,plevp,grid%jfirstxy:grid%jlastxy)         ! interface pressure
+    real(r8),   intent(in ) :: pmid(grid%ifirstxy:grid%ilastxy   ,plev,grid%jfirstxy:grid%jlastxy)       ! mid-point pressure
+    real(r8),   intent(in ) :: data(grid%ifirstxy:grid%ilastxy   ,plev,grid%jfirstxy:grid%jlastxy)
+    real(r8),   intent(out) :: ddata_dp(grid%ifirstxy:grid%ilastxy   ,plev,grid%jfirstxy:grid%jlastxy)
+    !---------------------------------------------------------------------------
+    integer :: i,j,k
+    real(r8) :: pint_above ! pressure interpolated to interface above the current k mid-point
+    real(r8) :: pint_below ! pressure interpolated to interface below the current k mid-point
+    real(r8) :: dint_above ! data interpolated to interface above the current k mid-point
+    real(r8) :: dint_below ! data interpolated to interface below the current k mid-point
+    integer :: beglatxy, endlatxy, beglonxy, endlonxy
+    !---------------------------------------------------------------------------
+    beglatxy = grid%jfirstxy
+    endlatxy = grid%jlastxy
+    beglonxy = grid%ifirstxy
+    endlonxy = grid%ilastxy
+
+    !$omp parallel do private (i,j,k,pint_above,pint_below,dint_above,dint_below)
+    do j = beglatxy, endlatxy
+         do k = 2, plev-1
+            do i = beglonxy, endlonxy
+              pint_above = pint(i,k,j)
+              pint_below = pint(i,k+1,j)
+              dint_above = ( data(i,k-1,j) + data(i,k,j) ) / 2.0_r8
+              dint_below = ( data(i,k+1,j) + data(i,k,j) ) / 2.0_r8
+              ddata_dp(i,k,j) = ( dint_above - dint_below ) / ( pint_above - pint_below )
+            end do
+         end do
+    end do
+
+    k = 1
+    !$omp parallel do private (i,j,pint_above,pint_below,dint_above,dint_below)
+    do j = beglatxy, endlatxy
+         do i = beglonxy, endlonxy
+              pint_above = pmid(i,k,j)
+              pint_below = pint(i,k+1,j)
+              dint_above = data(i,k,j)
+              dint_below = ( data(i,k+1,j) + data(i,k,j) ) / 2.0_r8
+              ddata_dp(i,k,j) = ( dint_above - dint_below ) / ( pint_above - pint_below )
+         end do
+    end do
+
+    k = plev
+    !$omp parallel do private (i,j,pint_above,pint_below,dint_above,dint_below)
+    do j = beglatxy, endlatxy
+         do i = beglonxy, endlonxy
+                pint_above = pint(i,k,j)
+                pint_below = pmid(i,k,j)
+                dint_above = ( data(i,k-1,j) + data(i,k,j) ) / 2.0_r8
+                dint_below = data(i,k,j)
+                ddata_dp(i,k,j) = ( dint_above - dint_below ) / ( pint_above - pint_below )
+         end do
+    end do
+  end subroutine compute_vertical_derivative
+
     subroutine ghost_array(grid, x, xg)
 
       ! subroutine added by Francis Vitt -- 7 July 2008

src/dynamics/se/dycore/prim_state_mod.F90

@@ -27,6 +27,7 @@ subroutine prim_printstate(elem, tl,hybrid,nets,nete, fvm, omega_cn)
     use se_dyn_time_mod,        only: tstep
     use control_mod,            only: rsplit, qsplit
     use perf_mod,       only: t_startf, t_stopf
+    use hycoef,                 only: hyai, ps0
     type (element_t),             intent(inout) :: elem(:)
     type (TimeLevel_t), target,   intent(in)    :: tl
     type (hybrid_t),              intent(in)    :: hybrid
@@ -62,7 +63,7 @@ subroutine prim_printstate(elem, tl,hybrid,nets,nete, fvm, omega_cn)
     ! moist surface pressure
     if (use_cslam) then
       do ie=nets,nete
-        moist_ps_fvm(:,:,ie)=SUM(fvm(ie)%dp_fvm(1:nc,1:nc,:),DIM=3)
+        moist_ps_fvm(:,:,ie)=SUM(fvm(ie)%dp_fvm(1:nc,1:nc,:),DIM=3)+hyai(1)*ps0
         do q=dry_air_species_num+1,thermodynamic_active_species_num
           m_cnst = thermodynamic_active_species_idx(q)
           do k=1,nlev
@@ -134,8 +135,8 @@ subroutine prim_printstate(elem, tl,hybrid,nets,nete, fvm, omega_cn)
         max_local(ie,5)  = 0.0_r8
       end if
       if (use_cslam) then
-        min_local(ie,6) = MINVAL(SUM(fvm(ie)%dp_fvm(1:nc,1:nc,:),DIM=3))
-        max_local(ie,6) = MAXVAL(SUM(fvm(ie)%dp_fvm(1:nc,1:nc,:),DIM=3))
+        min_local(ie,6) = MINVAL(SUM(fvm(ie)%dp_fvm(1:nc,1:nc,:),DIM=3))+hyai(1)*ps0
+        max_local(ie,6) = MAXVAL(SUM(fvm(ie)%dp_fvm(1:nc,1:nc,:),DIM=3))+hyai(1)*ps0
         min_local(ie,7) = MINVAL(moist_ps_fvm(:,:,ie))
         max_local(ie,7) = MINVAL(moist_ps_fvm(:,:,ie))
         min_local(ie,8)  = MINVAL(elem(ie)%state%psdry(:,:))
@@ -207,7 +208,7 @@ subroutine prim_printstate(elem, tl,hybrid,nets,nete, fvm, omega_cn)
           tmp_fvm(:,:,q,ie) = SUM(fvm(ie)%c(1:nc,1:nc,:,q)*fvm(ie)%dp_fvm(1:nc,1:nc,:),DIM=3)
         end do
         q=statediag_numtrac+1
-        tmp_fvm(:,:,q,ie) = SUM(fvm(ie)%dp_fvm(1:nc,1:nc,:),DIM=3)
+        tmp_fvm(:,:,q,ie) = SUM(fvm(ie)%dp_fvm(1:nc,1:nc,:),DIM=3)+hyai(1)*ps0
         q=statediag_numtrac+2
         tmp_fvm(:,:,q,ie) = moist_ps_fvm(:,:,ie)
       end do
@@ -291,7 +292,9 @@ subroutine prim_printstate(elem, tl,hybrid,nets,nete, fvm, omega_cn)
 101 format (A12,A23,A23,A23,A23)
 
 #ifdef waccm_debug
-    call prim_printstate_cslam_gamma(elem, tl,hybrid,nets,nete, fvm)
+    if (use_cslam) then
+       call prim_printstate_cslam_gamma(elem, tl,hybrid,nets,nete, fvm)
+    end if
 #endif
     call prim_printstate_U(elem, tl,hybrid,nets,nete, fvm)
   end subroutine prim_printstate