Merge branch 'whannah/eam/test-frontogenesis-correction'(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:
  rename altitude => geopotential for clarity
  add ifdef MODEL_THETA_L exceptions for phi calculations in  frontogenesis function
  dd preqx version of get_hydro_pressure_i
  fix namelist variable description
  update compute_frontogenesis for zgrad correction
  add get_hydro_pressure_interface subroutine
  add use_fgf_zgrad_correction namelist variable
  update pressure calculation
  update comments
  switch convention for clarity
  add namelist control for frontogenesis pgrad correction
  add comments describing tensor vector contraction code
  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

@@ -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/eam/src/physics/cam/phys_control.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/preqx/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)