Frontogenesis function pressure gradient correction

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/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

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

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)