Merge branch 'whannah/eam/test-frontogenesis-correction' into next (PR #7797)

Frontogenesis function pressure gradient correction

This implements a missing correction for the frontogenesis function (FGF) used for the frontal gravity wave drag scheme. The FGF gradient calculations are done on model levels for convenience, but the original papers indicate a simple correction should be used to effectively put these gradients on height or pressure levels. However, going back to the original implementation in CESM with FV dynamics this correction was never implemented, and this mistake was carried through to E3SM with the SE dynamics.

The correction addresses problematic biases of excessive gravity wave drag over areas of steep topography. It is controlled by a new namelist variable "use_fgf_pgrad_correction" that is off by default to preserve the behavior of official v3 release version.

[BFB]

* whannah/eam/test-frontogenesis-correction:
  comment fix
  fix namelist variable description
  misc fixes
  misc fixes
  fix white space
  update compute_frontogenesis for zgrad correction
  add get_hydro_pressure_interface subroutine
  add use_fgf_zgrad_correction namelist variable
  remove pmid_hydro
  update pressure calculation
  update comments
  switch convention for clarity
  add namelist control for frontogenesis pgrad correction
  fix typos
  add comments describing tensor vector contraction code
  bug fix
  add experimental pressure gradient correction for frontogensis function calculation

Author: singhbalwinder

components/eam/bld/build-namelist

@@ -4137,6 +4137,9 @@ if ($nl->get_value('use_gw_convect') =~ /$TRUE/io) {
     add_default($nl, 'gw_drag_file');
 }
 
+# this correction is disabled in the released version E3SMv3 and prior versions
+add_default($nl, 'use_fgf_pgrad_correction', 'val'=>'.false.');
+add_default($nl, 'use_fgf_zgrad_correction', 'val'=>'.false.');
 add_default($nl, 'use_gw_energy_fix'    , 'val'=>'.true.');
 
 

components/eam/bld/namelist_files/namelist_definition.xml

@@ -1072,6 +1072,18 @@ Whether or not to enable gravity waves produced by convection.
 Default: set by build-namelist.
 </entry>
 
+<entry id="use_fgf_pgrad_correction" type="logical" category="gw_drag"
+       group="phys_ctl_nl" valid_values="" >
+Flag to enable frontogenesis pressure gradient correction for frontal GWD
+Default: set by build-namelist.
+</entry>
+
+<entry id="use_fgf_zgrad_correction" type="logical" category="gw_drag"
+       group="phys_ctl_nl" valid_values="" >
+Flag to enable frontogenesis geopotential gradient correction for frontal GWD
+Default: set by build-namelist.
+</entry>
+
 <entry id="use_gw_energy_fix" type="logical" category="gw_drag"
        group="phys_ctl_nl" valid_values="" >
 Whether or not to enable GWD brute-force energy fix.

components/eam/src/dynamics/se/dp_coupling.F90

@@ -38,7 +38,7 @@ subroutine d_p_coupling(phys_state, phys_tend,  pbuf2d, dyn_out)
     use cam_abortutils,          only: endrun
     use gravity_waves_sources,   only: gws_src_fnct
     use dyn_comp,                only: frontgf_idx, frontga_idx, hvcoord
-    use phys_control,            only: use_gw_front
+    use phys_control,            only: use_gw_front, use_fgf_pgrad_correction, use_fgf_zgrad_correction
     use dyn_comp,                only: dom_mt
     use gllfvremap_mod,          only: gfr_dyn_to_fv_phys
 
@@ -107,7 +107,10 @@ subroutine d_p_coupling(phys_state, phys_tend,  pbuf2d, dyn_out)
       elem => dyn_out%elem
       tl_f = TimeLevel%n0  ! time split physics (with forward-in-time RK)
 
-      if (use_gw_front) call gws_src_fnct(elem, tl_f, nphys, frontgf, frontga)
+      if (use_gw_front) call gws_src_fnct(elem, tl_f, nphys, &
+                                          use_fgf_pgrad_correction, &
+                                          use_fgf_zgrad_correction, &
+                                          frontgf, frontga)
 
       if (fv_nphys > 0) then
         !-----------------------------------------------------------------------

components/eam/src/dynamics/se/gravity_waves_sources.F90

@@ -1,6 +1,6 @@
 module gravity_waves_sources
   use derivative_mod, only : derivative_t
-  use dimensions_mod, only : np,nlev
+  use dimensions_mod, only : np,nlev,nlevp
   use edgetype_mod, only       : EdgeBuffer_t
   use element_mod, only    : element_t
   use hybrid_mod, only     : hybrid_t
@@ -38,7 +38,7 @@ subroutine gws_init(elem)
 
   end subroutine gws_init
   !-------------------------------------------------------------------------------------------------
-  subroutine gws_src_fnct(elem, tl, nphys, frontgf, frontga)
+  subroutine gws_src_fnct(elem, tl, nphys, use_fgf_pgrad_correction, use_fgf_zgrad_correction, frontgf, frontga)
     use dimensions_mod, only  : npsq, nelemd
     use dof_mod, only         : UniquePoints
     use dyn_comp, only        : dom_mt
@@ -51,6 +51,8 @@ subroutine gws_src_fnct(elem, tl, nphys, frontgf, frontga)
     type (element_t),      intent(inout), dimension(:) :: elem
     integer,               intent(in   ) :: tl
     integer,               intent(in   ) :: nphys
+    logical,               intent(in   ) :: use_fgf_pgrad_correction
+    logical,               intent(in   ) :: use_fgf_zgrad_correction
     real (kind=real_kind), intent(out  ) :: frontgf(nphys*nphys,pver,nelemd)
     real (kind=real_kind), intent(out  ) :: frontga(nphys*nphys,pver,nelemd)
 
@@ -67,7 +69,10 @@ subroutine gws_src_fnct(elem, tl, nphys, frontgf, frontga)
     hybrid = hybrid_create(par,ithr,hthreads)
     allocate(frontgf_thr(nphys,nphys,nlev,nets:nete))
     allocate(frontga_thr(nphys,nphys,nlev,nets:nete))
-    call compute_frontogenesis(frontgf_thr,frontga_thr,tl,elem,hybrid,nets,nete,nphys)
+    call compute_frontogenesis( frontgf_thr, frontga_thr, tl, &
+                                use_fgf_pgrad_correction, &
+                                use_fgf_zgrad_correction, &
+                                elem, hybrid, nets, nete, nphys )
     if (fv_nphys>0) then
       do ie = nets,nete
         do k = 1,nlev
@@ -88,7 +93,10 @@ subroutine gws_src_fnct(elem, tl, nphys, frontgf, frontga)
 
   end subroutine gws_src_fnct
   !-------------------------------------------------------------------------------------------------
-  subroutine compute_frontogenesis(frontgf,frontga,tl,elem,hybrid,nets,nete,nphys)
+  subroutine compute_frontogenesis( frontgf, frontga, tl, &
+                                    use_fgf_pgrad_correction, &
+                                    use_fgf_zgrad_correction, &
+                                    elem, hybrid, nets, nete, nphys )
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! compute frontogenesis function F
   !   F =  -gradth dot C
@@ -99,63 +107,180 @@ subroutine compute_frontogenesis(frontgf,frontga,tl,elem,hybrid,nets,nete,nphys)
   ! 
   ! Original by Mark Taylor, July 2011
   ! Change by Santos, 10 Aug 2011:
-  ! Integrated into gravity_waves_sources module, several arguments made global
-  !  to prevent repeated allocation/initialization
+  !   Integrated into gravity_waves_sources module, several arguments made global
+  !   to prevent repeated allocation/initialization
   ! Change by Aaron Donahue, April 2017:
   !   Fixed bug where boundary information was called for processors not associated
   !   with dynamics when dyn_npes<npes
+  ! Change by Walter Hannah, Oct 2025:
+  !   added pressure gradient correction
   ! 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    use physical_constants, only: kappa
+    use physical_constants, only: kappa, g
     use derivative_mod,     only: gradient_sphere, ugradv_sphere
-    use edge_mod,           only : edge_g, edgevpack_nlyr, edgevunpack_nlyr
+    use edge_mod,           only: edge_g, edgevpack_nlyr, edgevunpack_nlyr
     use bndry_mod,          only: bndry_exchangev
     use dyn_comp,           only: hvcoord
     use spmd_utils,         only: iam 
     use parallel_mod,       only: par 
-    use element_ops,        only : get_temperature
+    use element_ops,        only: get_temperature, get_hydro_pressure_i, get_phi
     use dyn_grid,           only: fv_nphys
-    use prim_driver_mod, only     : deriv1
-    use element_ops,        only : get_temperature
-    use gllfvremap_mod, only  : gfr_g2f_scalar, gfr_g2f_vector
+    use prim_driver_mod,    only: deriv1
+    use gllfvremap_mod,     only: gfr_g2f_scalar, gfr_g2f_vector
     implicit none
     type(hybrid_t),        intent(in   ) :: hybrid
     type(element_t),target,intent(inout) :: elem(:)
     integer,               intent(in   ) :: nets,nete,nphys
     integer,               intent(in   ) :: tl ! timelevel to use
+    logical,               intent(in   ) :: use_fgf_pgrad_correction
+    logical,               intent(in   ) :: use_fgf_zgrad_correction
     real(kind=real_kind),  intent(out  ) :: frontgf(nphys,nphys,nlev,nets:nete)
     real(kind=real_kind),  intent(out  ) :: frontga(nphys,nphys,nlev,nets:nete)
 
     ! local
     integer :: k,kptr,i,j,ie,component
     real(kind=real_kind) :: frontgf_gll(np,np,nlev,nets:nete)
     real(kind=real_kind) :: gradth_gll(np,np,2,nlev,nets:nete)  ! grad(theta)
-    real(kind=real_kind) :: p(np,np)        ! pressure at mid points
-    real(kind=real_kind) :: theta(np,np)    ! potential temperature at mid points
-    real(kind=real_kind) :: temperature(np,np,nlev)  ! Temperature
+    real(kind=real_kind) :: zint(np,np,nlevp)           ! interface altitude
+    real(kind=real_kind) :: zmid(np,np,nlev)            ! mid-point altitude
+    real(kind=real_kind) :: pint(np,np,nlevp)           ! interface hydrostatic pressure
+    real(kind=real_kind) :: pmid(np,np,nlev)            ! mid-point hydrostatic pressure
+    real(kind=real_kind) :: temperature(np,np,nlev)     ! Temperature
     real(kind=real_kind) :: C(np,np,2), wf1(nphys*nphys,nlev), wf2(nphys*nphys,nlev)
+
+    real(kind=real_kind) :: theta(np,np,nlev)           ! potential temperature at mid points
+    real(kind=real_kind) :: dum_grad(np,np,2)           ! temporary variable for spherical gradient calcualtions
+    real(kind=real_kind) :: dum_cart(np,np,3,nlev)      ! temporary variable for cartesian gradient calcualtions
+
+    ! variables needed for eta to pressure surface correction
+    real(kind=real_kind) :: grad_vert_gll(np,np,2,nlev) ! grad(vertical coordinate) - either pressure or geopotential
+    real(kind=real_kind) :: theta_dvert(np,np,nlev)     ! d(theta)/dp for eta to pressure surface correction
+    real(kind=real_kind) :: dum_cart_dvert(np,np,3,nlev)! vertical pressure derivative of dum_cart
+
+    !---------------------------------------------------------------------------
+    !
+    ! 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_z(v) (or grad_p(v))
+    ! - the gradient on z-surfaces expressed in terms of the gradient on model
+    ! surfaces and a vertical geopotential gradient.
+    !
+    ! The old version only computed gradients on model surfaces, which creates
+    ! issues around topograpy. This is address with use_fgf_pgrad_correction=.true.
+    !
+    ! First, v is represented in cartesian coordinates  v(c) for c=1,2,3
+    ! For each v(c), we compute its gradient on z (or p) surfaces via:
+    !    grad(v(c)) - d(v(c))/dz grad(z)
+    ! 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
-      do k = 1,nlev
-        ! pressure at mid points
-        p(:,:) = hvcoord%hyam(k)*hvcoord%ps0 + hvcoord%hybm(k)*elem(ie)%state%ps_v(:,:,tl)
-        ! potential temperature: theta = T (p/p0)^kappa
+
+      if (use_fgf_pgrad_correction .or. use_fgf_zgrad_correction) then
+
+        ! compute pressure at interfaces and mid-points
+        call get_hydro_pressure_i(pint,elem(ie)%state%dp3d(:,:,:,tl),hvcoord)
+
         call get_temperature(elem(ie),temperature,hvcoord,tl)
-        theta(:,:) = temperature(:,:,k)*(psurf_ref / p(:,:))**kappa
-        gradth_gll(:,:,:,k,ie) = gradient_sphere(theta,deriv1,elem(ie)%Dinv)
-        ! compute C = (grad(theta) dot grad ) u
-        C(:,:,:) = ugradv_sphere(gradth_gll(:,:,:,k,ie), elem(ie)%state%v(:,:,:,k,tl),deriv1,elem(ie))
-        ! gradth_gll dot C
-        frontgf_gll(:,:,k,ie) = -( C(:,:,1)*gradth_gll(:,:,1,k,ie) + C(:,:,2)*gradth_gll(:,:,2,k,ie)  )
-        ! apply mass matrix
-        gradth_gll(:,:,1,k,ie) = gradth_gll(:,:,1,k,ie)*elem(ie)%spheremp(:,:)
-        gradth_gll(:,:,2,k,ie) = gradth_gll(:,:,2,k,ie)*elem(ie)%spheremp(:,:)
-        frontgf_gll(:,:,k,ie)  = frontgf_gll(:,:,k,ie)*elem(ie)%spheremp(:,:)
-      end do ! k
+
+        ! calculate pmid and potential temperature: theta = T (p/p0)^kappa
+        do k = 1,nlev
+          pmid(:,:,k) = pint(:,:,k) + elem(ie)%state%dp3d(:,:,k,tl)/2
+          theta(:,:,k) = temperature(:,:,k)*(psurf_ref / pmid(:,:,k))**kappa
+        end do
+
+        if (use_fgf_pgrad_correction) then
+          ! compute d(theta)/dp
+          call compute_vertical_derivative(pint,theta,theta_dvert)
+        end if
+
+        if (use_fgf_zgrad_correction) then
+          ! compute geopotential
+          call get_phi(elem(ie), zmid, zint, hvcoord, tl)
+          ! compute d(theta)/dz
+          call compute_vertical_derivative(zint,theta,theta_dvert)
+        end if
+
+        do k = 1,nlev
+          gradth_gll(:,:,:,k,ie) = gradient_sphere(theta(:,:,k),deriv1,elem(ie)%Dinv)
+          if (use_fgf_pgrad_correction) grad_vert_gll(:,:,:,k) = gradient_sphere(pmid(:,:,k),deriv1,elem(ie)%Dinv)
+          if (use_fgf_zgrad_correction) grad_vert_gll(:,:,:,k) = gradient_sphere(zmid(:,:,k),deriv1,elem(ie)%Dinv)
+          do component=1,2
+            gradth_gll(:,:,component,k,ie) = gradth_gll(:,:,component,k,ie) - theta_dvert(:,:,k) * grad_vert_gll(:,:,component,k)
+          end do
+        end do
+
+        do k = 1,nlev
+          ! latlon -> cartesian - Summing along the third dimension is a sum over components for each point
+          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
+          if (use_fgf_pgrad_correction) call compute_vertical_derivative(pint,dum_cart(:,:,component,:),dum_cart_dvert(:,:,component,:))
+          if (use_fgf_zgrad_correction) call compute_vertical_derivative(zint,dum_cart(:,:,component,:),dum_cart_dvert(:,:,component,:))
+        end do
+
+        do k = 1,nlev
+          ! Do ugradv on the cartesian components - Dot u with the gradient of each component
+          do component=1,3
+            dum_grad(:,:,:) = gradient_sphere(dum_cart(:,:,component,k),deriv1,elem(ie)%Dinv)
+            do i=1,2
+              dum_grad(:,:,i) = dum_grad(:,:,i) - dum_cart_dvert(:,:,component,k) * grad_vert_gll(:,:,i,k)
+            end do
+            dum_cart(:,:,component,k) = sum( gradth_gll(:,:,:,k,ie) * dum_grad ,3)
+          enddo
+          ! cartesian -> latlon - vec_sphere2cart is its own pseudoinverse.
+          do i=1,2
+            C(:,:,i) = sum(dum_cart(:,:,:,k)*elem(ie)%vec_sphere2cart(:,:,:,i), 3)
+          end do
+          ! gradth_gll dot C
+          frontgf_gll(:,:,k,ie) = -( C(:,:,1)*gradth_gll(:,:,1,k,ie) + C(:,:,2)*gradth_gll(:,:,2,k,ie)  )
+          ! apply mass matrix
+          gradth_gll(:,:,1,k,ie) = gradth_gll(:,:,1,k,ie)*elem(ie)%spheremp(:,:)
+          gradth_gll(:,:,2,k,ie) = gradth_gll(:,:,2,k,ie)*elem(ie)%spheremp(:,:)
+          frontgf_gll(:,:,k,ie)  = frontgf_gll(:,:,k,ie)*elem(ie)%spheremp(:,:)
+        end do ! k
+
+      else ! .not. use_fgf_pgrad_correction
+
+        do k = 1,nlev
+          ! pressure at mid points - this pressure preserves the old behavior of E3SMv3 and prior
+          pmid(:,:,k) = hvcoord%hyam(k)*hvcoord%ps0 + hvcoord%hybm(k)*elem(ie)%state%ps_v(:,:,tl)
+          ! potential temperature: theta = T (p/p0)^kappa
+          call get_temperature(elem(ie),temperature,hvcoord,tl)
+          theta(:,:,k) = temperature(:,:,k)*(psurf_ref / pmid(:,:,k))**kappa
+          gradth_gll(:,:,:,k,ie) = gradient_sphere(theta(:,:,k),deriv1,elem(ie)%Dinv)
+          ! compute C = (grad(theta) dot grad ) u
+          C(:,:,:) = ugradv_sphere(gradth_gll(:,:,:,k,ie), elem(ie)%state%v(:,:,:,k,tl),deriv1,elem(ie))
+          ! gradth_gll dot C
+          frontgf_gll(:,:,k,ie) = -( C(:,:,1)*gradth_gll(:,:,1,k,ie) + C(:,:,2)*gradth_gll(:,:,2,k,ie)  )
+          ! apply mass matrix
+          gradth_gll(:,:,1,k,ie) = gradth_gll(:,:,1,k,ie)*elem(ie)%spheremp(:,:)
+          gradth_gll(:,:,2,k,ie) = gradth_gll(:,:,2,k,ie)*elem(ie)%spheremp(:,:)
+          frontgf_gll(:,:,k,ie)  = frontgf_gll(:,:,k,ie)*elem(ie)%spheremp(:,:)
+        end do ! k
+
+      end if ! use_fgf_pgrad_correction
+
       ! pack
-       call edgeVpack_nlyr(edge_g, elem(ie)%desc, frontgf_gll(:,:,:,ie),nlev,0,3*nlev)
-       call edgeVpack_nlyr(edge_g, elem(ie)%desc, gradth_gll(:,:,:,:,ie),2*nlev,nlev,3*nlev)
+      call edgeVpack_nlyr(edge_g, elem(ie)%desc, frontgf_gll(:,:,:,ie),nlev,0,3*nlev)
+      call edgeVpack_nlyr(edge_g, elem(ie)%desc, gradth_gll(:,:,:,:,ie),2*nlev,nlev,3*nlev)
 
     end do ! ie
 
@@ -192,4 +317,37 @@ subroutine compute_frontogenesis(frontgf,frontga,tl,elem,hybrid,nets,nete,nphys)
 
   end subroutine compute_frontogenesis
   !-------------------------------------------------------------------------------------------------
+  subroutine compute_vertical_derivative(vert_int, data_mid, ddata_dvert)
+    !---------------------------------------------------------------------------
+    real(kind=real_kind),   intent(in ) :: vert_int(np,np,nlev) ! vertical coord on interfaces (i.e. pint or zint)
+    real(kind=real_kind),   intent(in ) :: data_mid(np,np,nlev) ! input data in mid-points
+    real(kind=real_kind),   intent(out) :: ddata_dvert(np,np,nlev) ! vertical derivative of data
+    !---------------------------------------------------------------------------
+    integer :: k
+    real(kind=real_kind) :: vert_above(np,np) ! pressure interpolated to interface above the current k mid-point
+    real(kind=real_kind) :: vert_below(np,np) ! pressure interpolated to interface below the current k mid-point
+    real(kind=real_kind) :: data_above(np,np) ! data interpolated to interface above the current k mid-point
+    real(kind=real_kind) :: data_below(np,np) ! data interpolated to interface below the current k mid-point
+    !---------------------------------------------------------------------------
+    do k = 1,nlev
+      if (k==1) then
+        data_above = data_mid(:,:,k)
+        data_below = ( data_mid(:,:,k+1) + data_mid(:,:,k) ) / 2.0 ! interpolate to interface k+1
+        vert_above = ( vert_int(:,:,k+1) + vert_int(:,:,k) ) / 2.0 ! interpolate to mid-point k
+        vert_below = vert_int(:,:,k+1)
+      elseif (k==nlev) then
+        data_above = ( data_mid(:,:,k-1) + data_mid(:,:,k) ) / 2.0 ! interpolate to interface
+        data_below = data_mid(:,:,k)
+        vert_above = vert_int(:,:,k)
+        vert_below = ( vert_int(:,:,k+1) + vert_int(:,:,k) ) / 2.0 ! interpolate to mid-point k
+      else
+        data_above = ( data_mid(:,:,k-1) + data_mid(:,:,k) ) / 2.0 ! interpolate to interface k
+        data_below = ( data_mid(:,:,k+1) + data_mid(:,:,k) ) / 2.0 ! interpolate to interface k+1
+        vert_above = vert_int(:,:,k)
+        vert_below = vert_int(:,:,k+1)
+      end if
+      ddata_dvert(:,:,k) = ( data_above - data_below ) / ( vert_above - vert_below )
+    end do
+  end subroutine compute_vertical_derivative
+  !-------------------------------------------------------------------------------------------------
 end module gravity_waves_sources