Add SQG vertical structure to eddy diffusivities (#738)

Added SQG vertical structure in Varmix to provide vertical profile for diffusivities - added function calc_sqg_struct in MOM_lateral_mixing_coeffs to compute sqg_struct - added sqg_expo to set the exponent of sqg_struct - to use sqg_struct for the backscatter, set BS_use_sqg=true, sqg_expo>0., and BS_EBT_power=0. - if SQG_USE_MEKE=True, use the eddy length scale from MEKE to compute sqg_struct - added eddy length scale Le in MEKE if SQG_USE_MEKE=True - added MEKE%Le into restart file if SQG_USE_MEKE=True - added MEKE in Varmix - registered N2_u and N2_v diagnostics when SQG_EXPO>0 (cherry picked from commit 6d3df0541c33d6f6d1f9fcb695f1a1eb961ec1b3) * Compute vertical structures for khth, khtr, backscatter, and kdgl90 all in VarMix - Vertical structures including khth_struct, khtr_struct, BS_struct, and kdgl90_struct are now computed in VarMix - Each diffusivity/viscosity have two vertical structure options, equivalent barotropic (EBT) and surface quasigeostrophic (SQG) mode structures - KHTH_USE_EBT_STRUCT, KHTR_USE_EBT_STRUCT, KDGL90_USE_EBT_STRUCT and BS_EBT_POWER parameters, which already existed, still control whether to use the EBT structure for khth, khtr, kdgl90, and backscatter, respectively - Added KHTH_USE_SQG_STRUCT, KHTR_USE_SQG_STRUCT, KDGL90_USE_SQG_STRUCT and BS_USE_SQG parameters to control whether to use the SQG structure for khth, khtr, kdgl90, and backscatter, respectively - If neither EBT nor SQG is called, no vertical structure will be used for that diffusivity/viscosity - An error will be called if both EBT and SQG structures are called for the same diffusivity/viscosity - Added dt as an input of calc_resoln_function. dt is needed for calc_sqg_struct called in calc_resoln_function Bug fixes - Changed ()**0.5 to sqrt in sqg_struct, which solves the dimension error - Added if statement for using BS_struct in MOM_hor_visc - Added if statement for CS%sqg_struct<=0. - Changed the subroundoff of Coriolis parameter to 1e-40 in calc_sqg_struct - Changed parameter BS_USE_SQG to BS_USE_SQG_STRUCT --------- Co-authored-by: Alistair Adcroft <[email protected]>
NOAA-GFDL · Nov 13, 2024 · 7a9adbc · 7a9adbc
1 parent 31a4d8b
commit 7a9adbc
Show file tree

Hide file tree

Showing 9 changed files with 379 additions and 78 deletions.
diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
@@ -751,7 +751,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
 
     if (CS%VarMix%use_variable_mixing) then
       call enable_averages(cycle_time, Time_start + real_to_time(US%T_to_s*cycle_time), CS%diag)
-      call calc_resoln_function(h, CS%tv, G, GV, US, CS%VarMix)
+      call calc_resoln_function(h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, dt)
       call calc_depth_function(G, CS%VarMix)
       call disable_averaging(CS%diag)
     endif
@@ -1899,7 +1899,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS
         if (.not. skip_diffusion) then
           if (CS%VarMix%use_variable_mixing) then
             call pass_var(CS%h, G%Domain)
-            call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix)
+            call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, dt_offline)
             call calc_depth_function(G, CS%VarMix)
             call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC)
           endif
@@ -1926,7 +1926,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS
         if (.not. skip_diffusion) then
           if (CS%VarMix%use_variable_mixing) then
             call pass_var(CS%h, G%Domain)
-            call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix)
+            call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, dt_offline)
             call calc_depth_function(G, CS%VarMix)
             call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC)
           endif

diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90
@@ -124,6 +124,7 @@ module MOM_MEKE
   logical :: debug      !< If true, write out checksums of data for debugging
   integer :: eke_src !< Enum specifying whether EKE is stepped forward prognostically (default),
                      !! read in from a file, or inferred via a neural network
+  logical :: sqg_use_MEKE !< If True, use MEKE%Le for the SQG vertical structure.
   type(diag_ctrl), pointer :: diag => NULL() !< A type that regulates diagnostics output
   !>@{ Diagnostic handles
   integer :: id_MEKE = -1, id_Ue = -1, id_Kh = -1, id_src = -1
@@ -400,6 +401,13 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
       call hchksum(LmixScale, 'MEKE LmixScale', G%HI, unscale=US%L_to_m)
     endif
 
+    if (allocated(MEKE%Le)) then
+      !$OMP parallel do default(shared)
+      do j=js,je ; do i=is,ie
+        MEKE%Le(i,j) = LmixScale(i,j)
+      enddo ; enddo
+    endif
+
     ! Aggregate sources of MEKE (background, frictional and GM)
     !$OMP parallel do default(shared)
     do j=js,je ; do i=is,ie
@@ -757,7 +765,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
     enddo ; enddo
   endif
 
-  if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au)) then
+  if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au) &
+      .or. allocated(MEKE%Le)) then
     call cpu_clock_begin(CS%id_clock_pass)
     call do_group_pass(CS%pass_Kh, G%Domain)
     call cpu_clock_end(CS%id_clock_pass)
@@ -1425,6 +1434,9 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME
                  "computing beta in the expression of Rhines scale. Use 1 if full "//&
                  "topographic beta effect is considered; use 0 if it's completely ignored.", &
                  units="nondim", default=0.0)
+  call get_param(param_file, mdl, "SQG_USE_MEKE", CS%sqg_use_MEKE, &
+                 "If true, the eddy scale of MEKE is used for the SQG vertical structure ",&
+                 default=.false.)
 
   ! Nonlocal module parameters
   call get_param(param_file, mdl, "CDRAG", cdrag, &
@@ -1531,13 +1543,20 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME
 
   CS%id_clock_pass = cpu_clock_id('(Ocean continuity halo updates)', grain=CLOCK_ROUTINE)
 
+
   ! Detect whether this instance of MEKE_init() is at the beginning of a run
   ! or after a restart. If at the beginning, we will initialize MEKE to a local
   ! equilibrium.
   CS%initialize = .not.query_initialized(MEKE%MEKE, "MEKE", restart_CS)
   if (coldStart) CS%initialize = .false.
   if (CS%initialize) call MOM_error(WARNING, &
                        "MEKE_init: Initializing MEKE with a local equilibrium balance.")
+  if (.not.query_initialized(MEKE%Le, "MEKE_Le", restart_CS) .and. allocated(MEKE%Le)) then
+    !$OMP parallel do default(shared)
+    do j=js,je ; do i=is,ie
+      MEKE%Le(i,j) = sqrt(G%areaT(i,j))
+    enddo ; enddo
+  endif
 
   ! Set up group passes.  In the case of a restart, these fields need a halo update now.
   if (allocated(MEKE%MEKE)) then
@@ -1548,8 +1567,10 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME
   if (allocated(MEKE%Kh)) call create_group_pass(CS%pass_Kh, MEKE%Kh, G%Domain)
   if (allocated(MEKE%Ku)) call create_group_pass(CS%pass_Kh, MEKE%Ku, G%Domain)
   if (allocated(MEKE%Au)) call create_group_pass(CS%pass_Kh, MEKE%Au, G%Domain)
+  if (allocated(MEKE%Le)) call create_group_pass(CS%pass_Kh, MEKE%Le, G%Domain)
 
-  if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au)) &
+  if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au) &
+      .or. allocated(MEKE%Le)) &
     call do_group_pass(CS%pass_Kh, G%Domain)
 
 end function MEKE_init
@@ -1839,6 +1860,7 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS)
   real :: MEKE_KHCoeff, MEKE_viscCoeff_Ku, MEKE_viscCoeff_Au  ! Coefficients for various terms [nondim]
   logical :: Use_KH_in_MEKE
   logical :: useMEKE
+  logical :: sqg_use_MEKE
   integer :: isd, ied, jsd, jed
 
 ! Determine whether this module will be used
@@ -1853,6 +1875,7 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS)
   MEKE_viscCoeff_Ku = 0. ; call read_param(param_file,"MEKE_VISCOSITY_COEFF_KU",MEKE_viscCoeff_Ku)
   MEKE_viscCoeff_Au = 0. ; call read_param(param_file,"MEKE_VISCOSITY_COEFF_AU",MEKE_viscCoeff_Au)
   Use_KH_in_MEKE = .false. ; call read_param(param_file,"USE_KH_IN_MEKE", Use_KH_in_MEKE)
+  sqg_use_MEKE = .false. ; call read_param(param_file,"SQG_USE_MEKE", sqg_use_MEKE)
 
   if (.not. useMEKE) return
 
@@ -1884,6 +1907,12 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS)
              longname="Lateral viscosity from Mesoscale Eddy Kinetic Energy", &
              units="m2 s-1", conversion=US%L_to_m**2*US%s_to_T)
   endif
+  if (sqg_use_MEKE) then
+    allocate(MEKE%Le(isd:ied,jsd:jed), source=0.0)
+    call register_restart_field(MEKE%Le, "MEKE_Le", .false., restart_CS, &
+             longname="Eddy length scale from Mesoscale Eddy Kinetic Energy", &
+             units="m", conversion=US%L_to_m)
+  endif
   if (Use_Kh_in_MEKE) then
     allocate(MEKE%Kh_diff(isd:ied,jsd:jed), source=0.0)
     call register_restart_field(MEKE%Kh_diff, "MEKE_Kh_diff", .false., restart_CS, &
@@ -1918,6 +1947,7 @@ subroutine MEKE_end(MEKE)
   if (allocated(MEKE%mom_src_bh)) deallocate(MEKE%mom_src_bh)
   if (allocated(MEKE%GM_src)) deallocate(MEKE%GM_src)
   if (allocated(MEKE%MEKE)) deallocate(MEKE%MEKE)
+  if (allocated(MEKE%Le)) deallocate(MEKE%Le)
 end subroutine MEKE_end
 
 !> \namespace mom_meke

diff --git a/src/parameterizations/lateral/MOM_MEKE_types.F90 b/src/parameterizations/lateral/MOM_MEKE_types.F90
@@ -24,6 +24,7 @@ module MOM_MEKE_types
                                     !! backscatter from unresolved eddies (see Jansen and Held, 2014).
   real, allocatable :: Au(:,:)      !< The MEKE-derived lateral biharmonic viscosity
                                     !! coefficient [L4 T-1 ~> m4 s-1].
+  real, allocatable :: Le(:,:)      !< Eddy length scale [L m]
 
   ! Parameters
   real :: KhTh_fac = 1.0 !< Multiplier to map Kh(MEKE) to KhTh [nondim]

diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90
@@ -446,6 +446,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
   logical :: use_MEKE_Ku
   logical :: use_MEKE_Au
   logical :: use_cont_huv
+  logical :: use_kh_struct
   integer :: is_vort, ie_vort, js_vort, je_vort  ! Loop ranges for vorticity terms
   integer :: is_Kh, ie_Kh, js_Kh, je_Kh  ! Loop ranges for thickness point viscosities
   integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
@@ -502,6 +503,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
   if (CS%id_FrictWorkIntz > 0) find_FrictWork = .true.
 
   if (allocated(MEKE%mom_src)) find_FrictWork = .true.
+  use_kh_struct = allocated(VarMix%BS_struct)
   backscat_subround = 0.0
   if (find_FrictWork .and. allocated(MEKE%mom_src) .and. (MEKE%backscatter_Ro_c > 0.0) .and. &
       (MEKE%backscatter_Ro_Pow /= 0.0)) &
@@ -663,7 +665,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
   !$OMP   is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, &
   !$OMP   is_vort, ie_vort, js_vort, je_vort, &
   !$OMP   is_Kh, ie_Kh, js_Kh, je_Kh, &
-  !$OMP   apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, &
+  !$OMP   apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, use_kh_struct, &
   !$OMP   use_MEKE_Ku, use_MEKE_Au, u_smooth, v_smooth, use_cont_huv, slope_x, slope_y, dz, &
   !$OMP   backscat_subround, GME_effic_h, GME_effic_q, &
   !$OMP   h_neglect, h_neglect3, inv_PI3, inv_PI6, &
@@ -1181,14 +1183,26 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
 
       if (use_MEKE_Ku .and. .not. CS%EY24_EBT_BS) then
         ! *Add* the MEKE contribution (which might be negative)
-        if (CS%res_scale_MEKE) then
-          do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
-            Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%Res_fn_h(i,j) * VarMix%BS_struct(i,j,k)
-          enddo ; enddo
+        if (use_kh_struct) then
+          if (CS%res_scale_MEKE) then
+            do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
+              Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%Res_fn_h(i,j) * VarMix%BS_struct(i,j,k)
+            enddo ; enddo
+          else
+            do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
+              Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%BS_struct(i,j,k)
+            enddo ; enddo
+          endif
         else
-          do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
-            Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%BS_struct(i,j,k)
-          enddo ; enddo
+          if (CS%res_scale_MEKE) then
+            do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
+              Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%Res_fn_h(i,j)
+            enddo ; enddo
+          else
+            do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
+              Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j)
+            enddo ; enddo
+          endif
         endif
       endif
 
@@ -1443,7 +1457,11 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
         if (visc_limit_h_flag(i,j,k) > 0) then
           Kh_BS(i,j) = 0.
         else
-          Kh_BS(i,j) = MEKE%Ku(i,j) * VarMix%BS_struct(i,j,k)
+          if (use_kh_struct) then
+            Kh_BS(i,j) = MEKE%Ku(i,j) * VarMix%BS_struct(i,j,k)
+          else
+            Kh_BS(i,j) = MEKE%Ku(i,j)
+          endif
         endif
       enddo ; enddo
 
@@ -1618,10 +1636,17 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
 
         if (use_MEKE_Ku .and. .not. CS%EY24_EBT_BS) then
           ! *Add* the MEKE contribution (might be negative)
-          Kh(I,J) = Kh(I,J) + 0.25*( ((MEKE%Ku(i,j)*VarMix%BS_struct(i,j,k)) + &
-                                     (MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + &
-                                     ((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + &
-                                     (MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) ) * meke_res_fn
+          if (use_kh_struct) then
+            Kh(I,J) = Kh(I,J) + 0.25*( ((MEKE%Ku(i,j)*VarMix%BS_struct(i,j,k)) + &
+                                       (MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + &
+                                       ((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + &
+                                       (MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) ) * meke_res_fn
+          else
+            Kh(I,J) = Kh(I,J) + 0.25*( (MEKE%Ku(i,j) + &
+                                       MEKE%Ku(i+1,j+1)) + &
+                                       (MEKE%Ku(i+1,j) + &
+                                        MEKE%Ku(i,j+1)) ) * meke_res_fn
+          endif
         endif
 
         if (CS%anisotropic) &
@@ -1789,10 +1814,17 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
         if (visc_limit_q_flag(I,J,k) > 0) then
           Kh_BS(I,J) = 0.
         else
-          Kh_BS(I,J) = 0.25*( ((MEKE%Ku(i,j)*VarMix%BS_struct(i,j,k)) + &
-                       (MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + &
-                       ((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + &
-                       (MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) )
+          if (use_kh_struct) then
+            Kh_BS(I,J) = 0.25*( ((MEKE%Ku(i,j)*VarMix%BS_struct(i,j,k)) + &
+                         (MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + &
+                         ((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + &
+                         (MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) )
+          else
+            Kh_BS(I,J) = 0.25*( (MEKE%Ku(i,j) + &
+                         MEKE%Ku(i+1,j+1)) + &
+                         (MEKE%Ku(i+1,j) + &
+                         MEKE%Ku(i,j+1)) )
+          endif
         endif
       enddo ; enddo