diff --git a/components/eam/bld/build-namelist b/components/eam/bld/build-namelist
index 4660e6b8d704..95416cb5e4ac 100755
--- a/components/eam/bld/build-namelist
+++ b/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.');
 
 
diff --git a/components/eam/bld/namelist_files/namelist_definition.xml b/components/eam/bld/namelist_files/namelist_definition.xml
index 0106c418ded0..258c3ef00e27 100644
--- a/components/eam/bld/namelist_files/namelist_definition.xml
+++ b/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.
diff --git a/components/eam/src/dynamics/se/dp_coupling.F90 b/components/eam/src/dynamics/se/dp_coupling.F90
index 5655948e2f39..149fff8b1c4c 100644
--- a/components/eam/src/dynamics/se/dp_coupling.F90
+++ b/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
         !-----------------------------------------------------------------------
diff --git a/components/eam/src/dynamics/se/gravity_waves_sources.F90 b/components/eam/src/dynamics/se/gravity_waves_sources.F90
index 16985e8083e0..c6fcc76dc1fe 100644
--- a/components/eam/src/dynamics/se/gravity_waves_sources.F90
+++ b/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
@@ -97,32 +105,45 @@ subroutine compute_frontogenesis(frontgf,frontga,tl,elem,hybrid,nets,nete,nphys)
   !   gradth = grad(theta)
   !   C = ( gradth dot grad ) U
   ! 
+  ! For more details on the frontal GW scheme and fronotogenesis source calculation:
+  !   Charron, M., and E. Manzini, 2002: Gravity Waves from Fronts:
+  !   Parameterization and Middle Atmosphere Response in a General
+  !   Circulation Model. J. Atmos. Sci., 59, 923–941
+  !
   ! 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 and Wandi Yu, Oct 2025:
+  !   added pressure gradient correction options for consistency with original C&M paper
   ! 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    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
+#ifdef MODEL_THETA_L
+    use element_ops,        only: get_phi
+#else
+    use element_ops,        only: get_phi, get_phi_i
+#endif
     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)
   
@@ -130,32 +151,151 @@ subroutine compute_frontogenesis(frontgf,frontga,tl,elem,hybrid,nets,nete,nphys)
     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 geopotential
+    real(kind=real_kind) :: zmid(np,np,nlev)            ! mid-point geopotential
+    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
+#ifdef MODEL_THETA_L
+          call get_phi(elem(ie), zmid, zint, hvcoord, tl)
+#else
+          call get_phi(elem(ie), zmid, hvcoord, tl, -1)
+          call get_phi_i(elem(ie), zint, hvcoord, tl, -1)
+#endif
+          ! 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 +332,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
diff --git a/components/eam/src/physics/cam/phys_control.F90 b/components/eam/src/physics/cam/phys_control.F90
index 501daee3855d..31734d9f636c 100644
--- a/components/eam/src/physics/cam/phys_control.F90
+++ b/components/eam/src/physics/cam/phys_control.F90
@@ -174,6 +174,9 @@ module phys_control
 ! Convective
 logical, public, protected :: use_gw_convect = .false.
 
+! Frontogenesis gradient correction options
+logical, public, protected :: use_fgf_pgrad_correction = .false.
+logical, public, protected :: use_fgf_zgrad_correction = .false.
 !GW energy fix
 logical, public, protected :: use_gw_energy_fix = .false.
 
@@ -258,7 +261,7 @@ subroutine phys_ctl_readnl(nlfile)
       use_qqflx_fixer, & 
       print_fixer_message, & 
       use_hetfrz_classnuc, use_gw_oro, use_gw_front, use_gw_convect, &
-      use_gw_energy_fix, &
+      use_fgf_pgrad_correction,use_fgf_zgrad_correction, use_gw_energy_fix, &
       use_od_ls,use_od_bl,use_od_ss,use_od_fd,&
       cld_macmic_num_steps, micro_do_icesupersat, &
       fix_g1_err_ndrop, ssalt_tuning, resus_fix, convproc_do_aer, &
@@ -378,6 +381,8 @@ subroutine phys_ctl_readnl(nlfile)
    call mpibcast(use_hetfrz_classnuc,             1 , mpilog,  0, mpicom)
    call mpibcast(use_gw_oro,                      1 , mpilog,  0, mpicom)
    call mpibcast(use_gw_front,                    1 , mpilog,  0, mpicom)
+   call mpibcast(use_fgf_pgrad_correction,        1 , mpilog,  0, mpicom)
+   call mpibcast(use_fgf_zgrad_correction,        1 , mpilog,  0, mpicom)
    call mpibcast(use_gw_convect,                  1 , mpilog,  0, mpicom)
    call mpibcast(use_gw_energy_fix,               1 , mpilog,  0, mpicom)
    call mpibcast(use_od_ls,                       1 , mpilog,  0, mpicom)
@@ -429,6 +434,12 @@ subroutine phys_ctl_readnl(nlfile)
    ! Defaults for PBL and microphysics are set in build-namelist.  Check here that
    ! values have been set to guard against problems with hand edited namelists.
 
+   ! check frontal gwd gradient correction options
+   if (use_fgf_pgrad_correction .and. use_fgf_zgrad_correction) then
+      write(iulog,*)'phys_setopts: use_fgf_pgrad_correction and use_fgf_zgrad_correction cannot both be true'
+      call endrun('phys_setopts: use_fgf_pgrad_correction and use_fgf_zgrad_correction cannot both be true')
+   end if
+
    ! WACCM-X run option set in build-namelist. Check for valid values
    if (.not. (waccmx_opt == 'ionosphere' .or. waccmx_opt == 'neutral' .or. waccmx_opt == 'off')) then
       write(iulog,*)'waccm: illegal value of waccmx_opt:', waccmx_opt
diff --git a/components/homme/src/preqx/share/element_ops.F90 b/components/homme/src/preqx/share/element_ops.F90
index 01e31e5ea362..85ca6be5bb61 100644
--- a/components/homme/src/preqx/share/element_ops.F90
+++ b/components/homme/src/preqx/share/element_ops.F90
@@ -221,6 +221,24 @@ subroutine get_temperature(elem,temperature,hvcoord,nt)
   
   end subroutine get_temperature
 
+  !_____________________________________________________________________
+  subroutine get_hydro_pressure_i(p_i,dp,hvcoord)
+  !
+  implicit none
+
+  real (kind=real_kind), intent(out)  :: p_i(np,np,nlevp)
+  real (kind=real_kind), intent(in)   :: dp(np,np,nlev)
+  type (hvcoord_t),     intent(in)    :: hvcoord                      ! hybrid vertical coordinate struct
+
+  integer :: k
+
+  p_i(:,:,1)=hvcoord%hyai(1)*hvcoord%ps0
+  do k=1,nlev  ! SCAN
+     p_i(:,:,k+1)=p_i(:,:,k) + dp(:,:,k)
+  enddo
+
+  end subroutine get_hydro_pressure_i
+
   !_____________________________________________________________________
   subroutine copy_state(elem,nin,nout)
   implicit none
diff --git a/components/homme/src/theta-l/share/element_ops.F90 b/components/homme/src/theta-l/share/element_ops.F90
index 853d3be3085c..033c95609e97 100644
--- a/components/homme/src/theta-l/share/element_ops.F90
+++ b/components/homme/src/theta-l/share/element_ops.F90
@@ -43,6 +43,7 @@
 !  get_pottemp()
 !  get_dpnh_dp()
 !  get_hydro_pressure()
+!  get_hydro_pressure_i()
 !  get_nonhydro_pressure()
 !  get_phi()
 !  get_cp_star()
@@ -68,7 +69,7 @@ module element_ops
   type(elem_state_t), dimension(:), allocatable :: state0 ! storage for save_initial_state routine
 
   public get_field, get_field_i, get_state, get_pottemp
-  public get_temperature, get_phi, get_R_star, get_hydro_pressure
+  public get_temperature, get_phi, get_R_star, get_hydro_pressure, get_hydro_pressure_i
   public set_thermostate, set_state, set_state_i, set_elem_state
   public set_forcing_rayleigh_friction
   public initialize_reference_states, set_theta_ref
@@ -307,6 +308,24 @@ subroutine get_hydro_pressure(p,dp,hvcoord)
 
   end subroutine get_hydro_pressure
 
+  !_____________________________________________________________________
+  subroutine get_hydro_pressure_i(p_i,dp,hvcoord)
+  !
+  implicit none
+
+  real (kind=real_kind), intent(out)  :: p_i(np,np,nlevp)
+  real (kind=real_kind), intent(in)   :: dp(np,np,nlev)
+  type (hvcoord_t),     intent(in)    :: hvcoord                      ! hybrid vertical coordinate struct
+
+  integer :: k
+
+  p_i(:,:,1)=hvcoord%hyai(1)*hvcoord%ps0
+  do k=1,nlev  ! SCAN
+     p_i(:,:,k+1)=p_i(:,:,k) + dp(:,:,k)
+  enddo
+
+  end subroutine get_hydro_pressure_i
+
 
   !_____________________________________________________________________
   subroutine get_nonhydro_pressure(elem,pnh,exner,hvcoord,nt)