Attempt to use deformation radius for memory time scale

src/parameterizations/lateral/MOM_Zanna_Bolton.F90

@@ -16,6 +16,7 @@ module MOM_Zanna_Bolton
 use MOM_cpu_clock,     only : CLOCK_MODULE, CLOCK_ROUTINE
 use MOM_ANN,           only : ANN_init, ANN_apply, ANN_end, ANN_CS
 use MOM_error_handler, only : MOM_error, FATAL, WARNING
+use MOM_lateral_mixing_coeffs, only : VarMix_CS
 
 implicit none ; private
 
@@ -58,13 +59,19 @@ module MOM_Zanna_Bolton
                               !! 1 - memory for eddy forcing
                               !! 2 - memory for large scale flow
                               !! 3 - memory for both eddy forcing and large scale flow
+  integer   :: ZB_memory_dynamical_scale !< Select which dynamical length scale to infer memory time
+                              !! scale from:
+                              !! 0 - no memory
+                              !! 1 - deformation scale (old), which probably differs from option 1 by factor of sqrt(2)
+                              !! 2 - grid scale
+                              !! 3 - deformation radius 
   real      :: ZB_eddy_memory !< The eddy memory time scale [T  ~> s] 
   real      :: ZB_eddy_memory_const !< The non-dimensional factor for the eddy memory time scale
-                               !! that scales the deformation scale [nondim]
+                               !! that scales the dynamical time scale [nondim]
                                !! Typical range: 1-10
   real      :: ZB_large_scale_memory !< The memory time scale for the velocity gradients [T  ~> s] 
   real      :: ZB_large_scale_memory_const !< The non-dimensional factor for the memory time scale
-                               !! for the velocity gradients that scales the deformation scale [nondim]
+                               !! for the velocity gradients that scales the dynamical time scale [nondim]
                                !! Typical range: 1-10
   real      :: ZB_max_memory  !< The maximum eddy memory time scale; default: 1 year [T  ~> s]
   real      :: DT             !< The baroclinic model time step [T  ~> s]
@@ -76,8 +83,7 @@ module MOM_Zanna_Bolton
                             !! points including metric terms [T-1 ~> s-1]
       vort_xy,            & !< Vertical vorticity (dv/dx - du/dy) in q (CORNER)
                             !! points including metric terms [T-1 ~> s-1]
-      hq                   !< Thickness in CORNER points [H ~> m or kg m-2]
-
+      hq                    !< Thickness in CORNER points [H ~> m or kg m-2]
   real, dimension(:,:,:), allocatable :: &
           Txx,     & !< Subgrid stress xx component in h [L2 T-2 ~> m2 s-2]
           Tyy,     & !< Subgrid stress yy component in h [L2 T-2 ~> m2 s-2]
@@ -95,7 +101,8 @@ module MOM_Zanna_Bolton
 
   real, dimension(:,:), allocatable :: &
           kappa_h, & !< Scaling coefficient in h points [L2 ~> m2]
-          kappa_q    !< Scaling coefficient in q points [L2 ~> m2]
+          kappa_q, & !< Scaling coefficient in q points [L2 ~> m2]
+          Rd_dx_h    !< The deformation radius compared with the grid spacing [nondim].
 
   real, allocatable ::    &
         ICoriolis_h(:,:), &  !< Inverse Coriolis parameter at h points [T ~> s]
@@ -166,6 +173,7 @@ subroutine ZB_2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020)
   type(ZB2020_CS),         intent(inout) :: CS         !< ZB2020 control structure.
   logical,                 intent(out)   :: use_ZB2020 !< If true, turns on ZB scheme.
 
+
   real :: subroundoff_Cor     !> A negligible parameter which avoids division by zero
                               !! but small compared to Coriolis parameter [T-1 ~> s-1]
 
@@ -252,6 +260,13 @@ subroutine ZB_2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020)
                  "\t 2 - memory for large scale flow\n" //&
                  "\t 3 - memory for both eddy forcing and large scale flow", default=0)
 
+  call get_param(param_file, mdl, "ZB_MEMORY_DYNAMICAL_SCALE", CS%ZB_memory_dynamical_scale, &
+                 "Select which dynamical length scale to infer memory time scale from:\n" //&
+                 "\t 0 - no memory\n" //&
+                 "\t 1 - deformation scale (old), which probably differs from option 1 by factor of sqrt(2)\n" //&
+                 "\t 2 - grid scale\n" //&
+                 "\t 3 - deformation radius", default=0)
+
   call get_param(param_file, mdl, "ZB_EDDY_MEMORY", CS%ZB_eddy_memory, &
                  "The eddy memory time scale. "//&
                  "A value of zero means no memory.", units="s", default=0.0)
@@ -304,6 +319,13 @@ subroutine ZB_2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020)
               "MOM_Zanna_Bolton: only one of ZB_EDDY_MEMORY_CONST or ZB_EDDY_MEMORY can be larger than 0")
   if (CS%ZB_large_scale_memory_const > 0 .AND. CS%ZB_large_scale_memory > 0) call MOM_error(FATAL, &
               "MOM_Zanna_Bolton: only one of ZB_LARGE_SCALE_MEMORY_CONST or ZB_LARGE_SCALE_MEMORY can be larger than 0")
+
+  if (CS%ZB_eddy_memory_const > 0 .OR. CS%ZB_large_scale_memory_const > 0) then
+    if (CS%ZB_memory_dynamical_scale == 0) then
+       call MOM_error(FATAL, &
+            "MOM_Zanna_Bolton: ZB_MEMORY_DYNAMICAL_SCALE must be greater than 0 if ZB_EDDY_MEMORY_CONST > 0 or ZB_LARGE_SCALE_MEMORY_CONST > 0")
+    endif
+  endif 
 
   ! Register fields for output from this module.
   CS%diag => diag
@@ -408,6 +430,9 @@ subroutine ZB_2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020)
     endif
   endif
 
+  if (CS%ZB_memory_dynamical_scale == 3) then
+    allocate(CS%Rd_dx_h(SZI_(G),SZJ_(G))); CS%Rd_dx_h(:,:) = 0.
+  endif
 
 
   ! Precomputing the scaling coefficient
@@ -579,7 +604,7 @@ end subroutine ZB_copy_gradient_and_thickness
 !> Baroclinic Zanna-Bolton-2020 parameterization, see
 !! eq. 6 in https://laurezanna.github.io/files/Zanna-Bolton-2020.pdf
 subroutine Zanna_Bolton_2020(u, v, h, diffu, diffv, G, GV, CS, &
-                             dx2h, dy2h, dx2q, dy2q)
+                             dx2h, dy2h, dx2q, dy2q, VarMix)
   type(ocean_grid_type),         intent(inout)    :: G      !< The ocean's grid structure.
   type(verticalGrid_type),       intent(in)    :: GV     !< The ocean's vertical grid structure.
   type(ZB2020_CS),               intent(inout) :: CS     !< ZB2020 control structure.
@@ -605,6 +630,7 @@ subroutine Zanna_Bolton_2020(u, v, h, diffu, diffv, G, GV, CS, &
   real, dimension(SZIB_(G),SZJB_(G)),           &
                                  intent(in) :: dx2q, & !< dx^2 at q points [L2 ~> m2]
                                                dy2q    !< dy^2 at q points [L2 ~> m2]
+  type(VarMix_CS), target,                    intent(in)    :: VarMix !< Variable mixing coefficients
 
   real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &
     ZB2020u           !< Zonal acceleration due to convergence of
@@ -631,7 +657,7 @@ subroutine Zanna_Bolton_2020(u, v, h, diffu, diffv, G, GV, CS, &
 
   ! Introduce memory into velocity gradients if specified
   if (CS%ZB_memory_location == 2 .OR. CS%ZB_memory_location == 3) then
-    if (CS%ZB_memory_type > 0) call step_forward_ZB_memory(u, v, h, G, GV, CS, eddy=.false.)
+    if (CS%ZB_memory_type > 0) call step_forward_ZB_memory(u, v, h, G, GV, CS, VarMix, eddy=.false.)
   endif
 
   ! Compute the stress tensor given the
@@ -646,7 +672,7 @@ subroutine Zanna_Bolton_2020(u, v, h, diffu, diffv, G, GV, CS, &
   call filter_stress(G, GV, CS)
 
   if (CS%ZB_memory_location == 1 .OR. CS%ZB_memory_location == 3) then
-    if (CS%ZB_memory_type > 0) call step_forward_ZB_memory(u, v, h, G, GV, CS, eddy=.true.)
+    if (CS%ZB_memory_type > 0) call step_forward_ZB_memory(u, v, h, G, GV, CS, VarMix, eddy=.true.)
   endif
 
   ! Compute acceleration from a given stress tensor
@@ -698,7 +724,7 @@ end subroutine Zanna_Bolton_2020
 
 
 
-subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, eddy)
+subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, VarMix, eddy)
   type(ZB2020_CS),       intent(inout)  :: CS            !< ZB2020 control structure.
   type(ocean_grid_type), intent(inout)     :: G             !< The ocean's grid structure.
   type(verticalGrid_type),  intent(in)  :: GV            !< The ocean's vertical grid structure.
@@ -708,15 +734,17 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, eddy)
                          intent(inout)    :: v              !< The meridional velocity [L T-1 ~> m s-1].
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),  &
                                  intent(in) :: h         !< Layer thicknesses [H ~> m or kg m-2].
+  type(VarMix_CS), target,                    intent(in)    :: VarMix !< Variable mixing coefficients
   logical, intent(in) :: eddy                           !< If true, step forward memory on eddy forcing; 
                                                         ! if false, step forward memory on velocity gradients
 
   ! Local variables
 
   real :: sh_xy_h                                     !< Shearing strain interpolated to h point [T-1 ~> s-1]
+  real :: vort_xy_h                                   !< Relative vorticity interpolated to h point [T-1 ~> s-1]
   real :: sh_xx_q                                     !< Horizontal tension interpolated to q point [T-1 ~> s-1]
   real :: D                                           !< deformation rate [T-1 ~> s-1]
-  real :: tau                                         !< eddy memory time scale tau [T ~> s]
+  real :: tau, tau0                                   !< eddy memory time scale tau [T ~> s]
   integer :: is, ie, js, je, nz, Isq, Ieq, Jsq, Jeq
   integer :: i, j, k
   character(len=160) :: mesg  ! The text of an error message
@@ -759,19 +787,44 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, eddy)
         call advection(u, v, G, GV, CS, q=CS%vort_xy_with_memory)
       endif
   endif
+
+  if (CS%ZB_memory_dynamical_scale == 3) then
+      if (allocated(VarMix%Rd_dx_h) .AND. allocated(CS%Rd_dx_h)) then
+        do j=js-1,je+1 ; do i=is-1,ie+1
+          CS%Rd_dx_h(i,j) = VarMix%Rd_dx_h(i,j)
+        enddo ; enddo
+      endif
+  endif
 
   do k=1,nz
 
-    ! relaxation: this is the Eulerian part of the averaging
     do j=js-1,je+1 ; do i=is-1,ie+1
-      ! compute eddy memory time scale at cell center
-      sh_xy_h = 0.25 * ( (CS%sh_xy(I-1,J-1,k) + CS%sh_xy(I,J,k)) &
-                       + (CS%sh_xy(I-1,J,k) + CS%sh_xy(I,J-1,k)) )
-      D = sqrt(CS%sh_xx(i,j,k) * CS%sh_xx(i,j,k) + sh_xy_h * sh_xy_h)
 
+      ! compute dynamical time scale tau0
+      if (CS%ZB_memory_dynamical_scale > 0) then
+          ! compute eddy memory time scale at cell center
+          sh_xy_h = 0.25 * ( (CS%sh_xy(I-1,J-1,k) + CS%sh_xy(I,J,k)) &
+                       + (CS%sh_xy(I-1,J,k) + CS%sh_xy(I,J-1,k)) )
+          vort_xy_h = 0.25 * ( (CS%vort_xy(I-1,J-1,k) + CS%vort_xy(I,J,k)) &
+                       + (CS%vort_xy(I-1,J,k) + CS%vort_xy(I,J-1,k)) )
+          if (CS%ZB_memory_dynamical_scale == 1) then
+            D = sqrt(CS%sh_xx(i,j,k) * CS%sh_xx(i,j,k) + sh_xy_h * sh_xy_h)
+            tau0 = 1 / (D + 1.0/CS%ZB_max_memory)
+          elseif (CS%ZB_memory_dynamical_scale == 2) then
+            D = sqrt(CS%sh_xx(i,j,k) * CS%sh_xx(i,j,k) + sh_xy_h * sh_xy_h + vort_xy_h * vort_xy_h)
+            tau0 = 1 / (D + 1.0/CS%ZB_max_memory)
+          elseif (CS%ZB_memory_dynamical_scale == 3) then
+
+            D = sqrt(CS%sh_xx(i,j,k) * CS%sh_xx(i,j,k) + sh_xy_h * sh_xy_h + vort_xy_h * vort_xy_h)
+            tau0 = CS%Rd_dx_h(i,j) / (D + CS%Rd_dx_h(i,j)/CS%ZB_max_memory)
+          endif
+
+      endif
+
+      ! relaxation: this is the Eulerian part of the averaging
       if (eddy) then
         if (CS%ZB_eddy_memory_const > 0.0) then
-          tau = CS%ZB_eddy_memory_const / (D + 1.0/CS%ZB_max_memory)
+          tau = CS%ZB_eddy_memory_const * tau0
         else
           tau = CS%ZB_eddy_memory
         endif
@@ -788,7 +841,7 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, eddy)
 
       else
         if (CS%ZB_large_scale_memory_const > 0.0) then
-          tau = CS%ZB_large_scale_memory_const / (D + 1.0/CS%ZB_max_memory)
+          tau = CS%ZB_large_scale_memory_const * tau0
         else
           tau = CS%ZB_large_scale_memory
         endif

src/parameterizations/lateral/MOM_hor_visc.F90

@@ -1692,7 +1692,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
 
   if (CS%use_ZB2020) then
     call Zanna_Bolton_2020(u, v, h, diffu, diffv, G, GV, CS%ZB2020, &
-                           CS%dx2h, CS%dy2h, CS%dx2q, CS%dy2q)
+                           CS%dx2h, CS%dy2h, CS%dx2q, CS%dy2q, VarMix)
   endif
 
 end subroutine horizontal_viscosity

src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90

@@ -1109,6 +1109,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
   real :: absurdly_small_freq  ! A miniscule frequency that is used to avoid division by 0 [T-1 ~> s-1].  The
              ! default value is roughly (pi / (the age of the universe)).
   logical :: Gill_equatorial_Ld, use_FGNV_streamfn, use_MEKE, in_use
+  integer :: ZB_memory_dynamical_scale
   integer :: default_answer_date  ! The default setting for the various ANSWER_DATE flags.
   integer :: remap_answer_date    ! The vintage of the order of arithmetic and expressions to use
                                   ! for remapping.  Values below 20190101 recover the remapping
@@ -1180,6 +1181,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
                  "If true, the viscosity contribution from MEKE is scaled by "//&
                  "the resolution function.", default=.false., do_not_log=.true.) ! Logged elsewhere.
   if (.not.use_MEKE) Resoln_scaled_MEKE_visc = .false.
+
   call get_param(param_file, mdl, "RESOLN_USE_EBT", CS%Resoln_use_ebt, &
                  "If true, uses the equivalent barotropic wave speed instead "//&
                  "of first baroclinic wave for calculating the resolution fn.",&
@@ -1211,6 +1213,13 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
                  default=.false., do_not_log=.true.)
   CS%calculate_cg1 = CS%calculate_cg1 .or. use_FGNV_streamfn .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct
   CS%calculate_Rd_dx = CS%calculate_Rd_dx .or. use_MEKE
+  call get_param(param_file, mdl, "ZB_MEMORY_DYNAMICAL_SCALE", ZB_memory_dynamical_scale, &
+                 "Select which dynamical length scale to infer memory time scale from:\n" //&
+                 "\t 0 - no memory\n" //&
+                 "\t 1 - deformation scale (old), which probably differs from option 1 by factor of sqrt(2)\n" //&
+                 "\t 2 - grid scale\n" //&
+                 "\t 3 - deformation radius", default=0)
+  CS%calculate_Rd_dx = CS%calculate_Rd_dx .or. (ZB_memory_dynamical_scale == 3)
   ! Indicate whether to calculate the Eady growth rate
   CS%calculate_Eady_growth_rate = use_MEKE .or. (KhTr_Slope_Cff>0.) .or. (KhTh_Slope_Cff>0.)
   call get_param(param_file, mdl, "KHTR_PASSIVITY_COEFF", KhTr_passivity_coeff, &