Merge pull request #1415 from PeterHjortLauritzen/frontogenesis-bug

cam6_4_126: Fix frontogenesis bug (SE dycore only)

Author: cacraigucar

.gitmodules

@@ -107,13 +107,6 @@
         fxtag = clubb_4ncar_20240605_73d60f6_gpufixes_posinf
         fxDONOTUSEurl = https://github.com/larson-group/clubb_release
 
-[submodule "cism"]
-path = components/cism
-url = https://github.com/ESCOMP/CISM-wrapper
-fxtag = cismwrap_2_2_010
-fxrequired = ToplevelRequired
-fxDONOTUSEurl = https://github.com/ESCOMP/CISM-wrapper
-
 [submodule "rtm"]
 path = components/rtm
 url = https://github.com/ESCOMP/RTM

cime_config/config_compsets.xml

@@ -423,7 +423,7 @@
   </compset>
     <compset>
     <alias>FCHIST_GC</alias>
-    <lname>HIST_CAM60%GEOSCHEM%HEMCO_CLM50%SP_CICE%PRES_DOCN%DOM_MOSART_CISM2%NOEVOLVE_SWAV</lname>
+    <lname>HIST_CAM60%GEOSCHEM%HEMCO_CLM50%SP_CICE%PRES_DOCN%DOM_MOSART_SGLC_SWAV</lname>
   </compset>
 
   <compset>

components/cism

@@ -1,5 +1,97 @@
 ===============================================================
 
+Tag name: cam6_4_126
+Originator(s):  Walter Hannah, Jack Chen, Peter Lauritzen, Cheryl Craig
+Date: Oct 16, 2025
+One-line Summary: Fix frontogenesis function bug and Change CISM2 to SGLC in FCHIST_GC compset
+Github PR URL: https://github.com/ESCOMP/pull/1415 and 1417
+
+Purpose of changes (include the issue number and title text for each relevant GitHub issue):
+  - Fix bug in frontogenesis function computation: https://github.com/ESCOMP/CAM/issues/1414
+  - Change CISM2 to SGLC in FCHIST_GC compset: https://github.com/ESCOMP/CAM/issues/1416 and 1412
+
+Describe any changes made to build system: Using SGLC for FCHIST_GC compset
+
+Describe any changes made to the namelist: N/A
+
+List any changes to the defaults for the bound N/Aary datasets:
+
+Describe any substantial timing or memory changes: N/A
+
+Code reviewed by: Peter Lauritzen, Jack Chen, cacraigucar
+
+List all files eliminated:
+D       components/cism
+           - CISM external removed from CAM checkouts
+
+List all files added and what they do: N/A
+
+List all existing files that have been modified, and describe the changes:
+M       .gitmodules
+           - CISM external removed from CAM checkouts
+
+M       cime_config/config_compsets.xml
+           - FCHIST_GC compset now uses data ice
+
+M       src/dynamics/se/gravity_waves_sources.F90
+           - bug fixes for frontogenesis
+
+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:
+  ERC_D_Ln9.ne30pg2_ne30pg2_mt232.QPC7.derecho_intel.cam-outfrq9s (Overall: DIFF) details:
+  ERC_D_Ln9.ne30pg3_ne30pg3_mt232.F1850C_LTso.derecho_intel.cam-outfrq9s (Overall: DIFF) details:
+  ERI_D_Ln18.ne30pg3_ne30pg3_mt232.FHISTC_LTso.derecho_intel.cam-outfrq3s_eri (Overall: DIFF) details:
+  ERP_D_Ln9.ne30pg3_ne30pg3_mt232.F1850C_MTso.derecho_intel.cam-outfrq9s (Overall: DIFF) details:
+  ERP_Ld3.ne16pg3_ne16pg3_mg17.FHISTC_WAt1ma.derecho_intel.cam-reduced_hist1d (Overall: DIFF) details:
+  ERP_Ld3.ne30pg3_ne30pg3_mt232.FHISTC_MTt4s.derecho_intel.cam-outfrq1d_aoa (Overall: DIFF) details:
+  ERP_Ln9.ne30pg3_ne30pg3_mg17.FCnudged.derecho_intel.cam-outfrq9s (Overall: DIFF) details  :
+  ERP_Ln9.ne30pg3_ne30pg3_mg17.FHISTC_WAma.derecho_intel.cam-outfrq9s (Overall: DIFF) details:
+  ERS_Ln9.ne30pg3_ne30pg3_mg17.FHISTC_WXma.derecho_intel.cam-outfrq9s_ctem (Overall: DIFF) details:
+  SMS_C2_D_Ln9.ne16pg3_ne16pg3_mg17.FHISTC_WXma.derecho_intel.cam-outfrq9s (Overall: DIFF) details:
+  SMS_D_Ln9.ne16_ne16_mg17.QPX2000.derecho_intel.cam-outfrq9s_amie (Overall: DIFF) details:
+  SMS_D_Ln9.ne30pg3_ne30pg3_mt232.FHISTC_MTso.derecho_intel.cam-outfrq9s (Overall: DIFF) details:
+  SMS_D_Ln9_P1280x1.ne0CONUSne30x8_ne0CONUSne30x8_mt12.FCHIST.derecho_intel.cam-outfrq9s (Overall: DIFF) details:
+  SMS_D_Ln9_P1280x1.ne30pg3_ne30pg3_mt232.FHISTC_MTt1s.derecho_intel.cam-outfrq9s_Leung_dust (Overall: DIFF) details:
+  SMS_Ld1.ne30pg3_ne30pg3_mg17.FC2010climo.derecho_intel.cam-outfrq1d (Overall: DIFF) details:
+  SMS_Ln9.ne30pg3_ne30pg3_mg17.FW2000climo.derecho_intel.cam-outfrq9s_rrtmgp (Overall: DIFF) details:
+    - All CAM7 and WACCM tests with SE dycore should have answer changes.
+
+  SMS_Ld1.f09_f09_mg17.FCHIST_GC.derecho_intel.cam-outfrq1d (Overall: DIFF) details:
+    - Answer and namelist changes due to changing CISM to be data ice
+
+derecho/nvhpc/aux_cam:
+  ERS_Ln9.ne30pg3_ne30pg3_mt232.FHISTC_LTso.derecho_nvhpc.cam-outfrq9s_gpu_default (Overall: DIFF) details:
+    - All CAM7 and WACCM tests with SE dycore should have answer changes.
+
+izumi/nag/aux_cam:
+  SMS_D_Ln3.ne5pg3_ne5pg3_mg37.QPX2000.izumi_nag.cam-outfrq3s (Overall: DIFF) details:
+  SMS_D_Ln6.ne5_ne5_mg37.QPWmaC4.izumi_nag.cam-outfrq3s_physgrid_tem (Overall: DIFF) details:
+    - All CAM7 and WACCM tests with SE dycore should have answer changes.
+
+
+izumi/gnu/aux_cam: all BFB
+
+Summarize any changes to answers, i.e.,
+- what code configurations:
+    - All CAM7 and WACCM tests with SE dycore should have answer changes due to using the frontal scheme
+
+- what platforms/compilers:
+    - all
+
+- nature of change (roundoff; larger than roundoff but same climate; new
+  climate): climate changing
+
+If this tag changes climate describe the run(s) done to evaluate the new
+climate in enough detail that it(they) could be reproduced, i.e.,
+    - Answers were scientifically verified by Peter and Julio
+
+===============================================================
+===============================================================
+
 Tag name: cam6_4_125
 Originator(s): eaton
 Date: 9 October 2025
@@ -92,7 +184,7 @@ src/physics/cam/macrop_driver.F90
 src/physics/cam/micro_pumas_cam.F90
 src/physics/cam7/micro_pumas_cam.F90
 . remove pbuf fields am_evp_st, evprain_st, evpsnow_st.  these are only used by
-  unicon. 
+  unicon.
 
 src/physics/cam/physpkg.F90
 . remove call to unicon_cam_org_diags

doc/ChangeLog

@@ -247,6 +247,8 @@ subroutine compute_frontogenesis(frontgf,frontga,tl,tlq,elem,ederiv,hybrid,nets,
   ! Integrated into gravity_waves_sources module, several arguments made global
   !  to prevent repeated allocation/initialization
   !
+  ! Frontogenesis function correction by Walter Hannah, Mark Taylor, and Jack Chen. October 2025
+  !
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     use physconst,      only: cappa
     use air_composition,only: dry_air_species_num, thermodynamic_active_species_num
@@ -273,14 +275,46 @@ subroutine compute_frontogenesis(frontgf,frontga,tl,tlq,elem,ederiv,hybrid,nets,
     real(r8) :: frontga_gll(np,np,nlev,nets:nete)
     integer  :: k,kptr,i,j,ie,component,h,nq,m_cnst
     real(r8) :: gradth(np,np,2,nlev,nets:nete) ! grad(theta)
-    real(r8) :: p(np,np)                       ! pressure at mid points
-    real(r8) :: pint(np,np)                    ! pressure at interface points
-    real(r8) :: theta(np,np)                   ! potential temperature at mid points
+
+    real(r8) :: p(np,np,nlev)               ! pressure at mid points
+    real(r8) :: pint(np,np,nlev+1)          ! pressure at interface points
+    real(r8) :: gradp(np,np,2)              ! grad(pressure)
+    real(r8) :: theta(np,np,nlev)           ! potential temperature at mid points
+    real(r8) :: dtheta_dp(np,np,nlev)       ! d(theta)/dp    for eta to pressure surface correction
+    real(r8) :: dum_grad(np,np,2)           ! ?
+    real(r8) :: dum_cart(np,np,3,nlev)      ! d/dp of ?
+    real(r8) :: ddp_dum_cart(np,np,3,nlev)  ! ?
     real(r8) :: C(np,np,2), sum_water(np,np)
 
+    !  By Mark Taylor
+    !  For a vector velocity "v", a tensor "grad(v)", and a vector "grad(theta)",
+    !  this loop computes the vector "grad(theta)*grad(v)"
+    !
+    !  Representing the tensor "grad(v)" in spherical coordinates is difficult.  This routine
+    !  avoids this by computing a mathematically equivalent form using a mixture of
+    !  Cartesian and spherical coordinates
+    !
+    !  This routine is a modified version of derivative_mod.F90:ugradv_sphere() in that the
+    !  grad(v) term is modified to compute grad_p(v) - the gradient on p-surfaces expressed
+    !  in terms of the gradient on model surfaces and a vertical pressure gradient.
+    !
+    !  First, v is represented in cartesian coordinates  v(c) for c=1,2,3
+    !  For each v(c), we compute its gradient on p-surfaces via:
+    !     grad(v(c)) - d(v(c))/dz grad(p)
+    !  Each of these gradients is represented in *spherical* coordinates (i=1,2)
+    !
+    !  We then dot each of these vectors with grad(theta).  This dot product is computed
+    !  in spherical coordinates.  The end result is dum_cart(c), for c=1,2,3
+    !  These three scalars are the three Cartesian coefficients of
+    !  the vector "grad(theta)*grad(v)"
+    !
+    !  This Cartesian vector is then transformed back to spherical coordinates
+    !
+
     do ie=nets,nete
       ! pressure at model top
-      pint(:,:) = hvcoord%hyai(1)*hvcoord%ps0
+      pint(:,:,1) = hvcoord%hyai(1)*hvcoord%ps0
+
       do k=1,nlev
         ! moist pressure at mid points
         sum_water(:,:) = 1.0_r8
@@ -291,15 +325,49 @@ subroutine compute_frontogenesis(frontgf,frontga,tl,tlq,elem,ederiv,hybrid,nets,
           !
           sum_water(:,:) = sum_water(:,:) + elem(ie)%state%Qdp(:,:,k,m_cnst,tlq)/elem(ie)%state%dp3d(:,:,k,tl)
         end do
-        p(:,:) = pint(:,:) + 0.5_r8*sum_water(:,:)*elem(ie)%state%dp3d(:,:,k,tl)
+        p(:,:,k) = pint(:,:,k) + 0.5_r8*sum_water(:,:)*elem(ie)%state%dp3d(:,:,k,tl)
         ! moist pressure at interface for next iteration
-        pint(:,:) = pint(:,:)+elem(ie)%state%dp3d(:,:,k,tl)
+        pint(:,:,k+1) = pint(:,:,k)+elem(ie)%state%dp3d(:,:,k,tl)
         !
-        theta(:,:) = elem(ie)%state%T(:,:,k,tl)*(psurf_ref / p(:,:))**cappa
-        ! gradth(:,:,:,k,ie) = gradient_sphere(theta,ederiv,elem(ie)%Dinv)
-        call gradient_sphere(theta,ederiv,elem(ie)%Dinv,gradth(:,:,:,k,ie))
-        ! compute C = (grad(theta) dot grad ) u
-        C(:,:,:) = ugradv_sphere(gradth(:,:,:,k,ie), elem(ie)%state%v(:,:,:,k,tl),ederiv,elem(ie))
+        theta(:,:,k) = elem(ie)%state%T(:,:,k,tl)*(psurf_ref / p(:,:,k))**cappa
+      end do
+
+      call compute_vertical_derivative(pint,p,theta,dtheta_dp)
+
+      do k=1,nlev
+        call gradient_sphere(theta(:,:,k),ederiv,elem(ie)%Dinv,gradth(:,:,:,k,ie))
+
+        call gradient_sphere(p(:,:,k),ederiv,elem(ie)%Dinv,gradp)
+
+        do component=1,2
+          gradth(:,:,component,k,ie) = gradth(:,:,component,k,ie) - dtheta_dp(:,:,k) * gradp(:,:,component)
+        end do
+      end do
+
+      do k=1,nlev
+        do component=1,3
+          dum_cart(:,:,component,k) = sum( elem(ie)%vec_sphere2cart(:,:,component,:) * elem(ie)%state%v(:,:,:,k,tl),3 )
+        end do
+      end do
+
+      do component=1,3
+        call compute_vertical_derivative(pint,p,dum_cart(:,:,component,:),ddp_dum_cart(:,:,component,:))
+      end do
+      do k=1,nlev
+        call gradient_sphere(p(:,:,k),ederiv,elem(ie)%Dinv,gradp)
+
+        do component=1,3
+          call gradient_sphere(dum_cart(:,:,component,k),ederiv,elem(ie)%Dinv,dum_grad)
+          do i=1,2
+            dum_grad(:,:,i) = dum_grad(:,:,i) - ddp_dum_cart(:,:,component,k) * gradp(:,:,i)
+          end do
+          dum_cart(:,:,component,k) = sum( gradth(:,:,:,k,ie) * dum_grad , 3 )
+        end do
+
+        do i=1,2
+          C(:,:,i) = sum(dum_cart(:,:,:,k)*elem(ie)%vec_sphere2cart(:,:,:,i), 3)
+        end do
+
         ! gradth dot C
         frontgf_gll(:,:,k,ie) = -( C(:,:,1)*gradth(:,:,1,k,ie) +  C(:,:,2)*gradth(:,:,2,k,ie)  )
         ! apply mass matrix
@@ -351,5 +419,38 @@ subroutine compute_frontogenesis(frontgf,frontga,tl,tlq,elem,ederiv,hybrid,nets,
     enddo
   end subroutine compute_frontogenesis
 
+  subroutine compute_vertical_derivative(pint,pmid,data,ddata_dp)
+    !---------------------------------------------------------------------------
+    real(r8),   intent(in ) :: pint(np,np,nlev+1)
+    real(r8),   intent(in ) :: pmid(np,np,nlev)
+    real(r8),   intent(in ) :: data(np,np,nlev)
+    real(r8),   intent(out) :: ddata_dp(np,np,nlev)
+    !---------------------------------------------------------------------------
+    integer :: k
+    real(r8) :: pint_above(np,np) ! pressure interpolated to interface above the current k mid-point
+    real(r8) :: pint_below(np,np) ! pressure interpolated to interface below the current k mid-point
+    real(r8) :: dint_above(np,np) ! data interpolated to interface above the current k mid-point
+    real(r8) :: dint_below(np,np) ! data interpolated to interface below the current k mid-point
+    !---------------------------------------------------------------------------
+    do k = 1,nlev
+      if (k==1) then
+        pint_above = pmid(:,:,k)
+        pint_below = pint(:,:,k+1)
+        dint_above = data(:,:,k)
+        dint_below = ( data(:,:,k+1) + data(:,:,k) ) / 2.0_r8
+      elseif (k==nlev) then
+        pint_above = pint(:,:,k)
+        pint_below = pmid(:,:,k)
+        dint_above = ( data(:,:,k-1) + data(:,:,k) ) / 2.0_r8
+        dint_below = data(:,:,k)
+      else
+        pint_above = pint(:,:,k)
+        pint_below = pint(:,:,k+1)
+        dint_above = ( data(:,:,k-1) + data(:,:,k) ) / 2.0_r8
+        dint_below = ( data(:,:,k+1) + data(:,:,k) ) / 2.0_r8
+      end if
+      ddata_dp(:,:,k) = ( dint_above - dint_below ) / ( pint_above - pint_below )
+    end do
+  end subroutine compute_vertical_derivative
 
 end module gravity_waves_sources