Attempts to include Rossby radius

NoraLoose · Apr 4, 2024 · 36e57c5 · 36e57c5
1 parent 4470c6d
commit 36e57c5
Showing 1 changed file with 50 additions and 12 deletions.
diff --git a/src/parameterizations/lateral/MOM_Zanna_Bolton.F90 b/src/parameterizations/lateral/MOM_Zanna_Bolton.F90
@@ -58,13 +58,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]
@@ -251,7 +257,13 @@ subroutine ZB_2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020)
                  "\t 1 - memory for eddy forcing\n" //&
                  "\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)
@@ -305,6 +317,12 @@ subroutine ZB_2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020)
   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
 
@@ -704,7 +722,7 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, eddy)
   type(verticalGrid_type),  intent(in)  :: GV            !< The ocean's vertical grid structure.
   real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
                          intent(inout)    :: u              !< The zonal velocity [L T-1 ~> m s-1].
-  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
+  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), 
                          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].
@@ -714,9 +732,11 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, eddy)
   ! 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]
+  real :: Rd                                          !< first Rossby radius of deformation in baroclinic DG [L ~> m]
   integer :: is, ie, js, je, nz, Isq, Ieq, Jsq, Jeq
   integer :: i, j, k
   character(len=160) :: mesg  ! The text of an error message
@@ -761,17 +781,35 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, eddy)
   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)) &
+
+      ! 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)) )
-      D = sqrt(CS%sh_xx(i,j,k) * CS%sh_xx(i,j,k) + sh_xy_h * sh_xy_h)
+          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
+            subroundoff_Cor = 1e-30 * US%T_to_s
+            ICoriolis = 1. / (0.25 * ((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) &
+                          + (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))) + subroundoff_Cor)
+            Rd = sqrt(GV%g_prime(2) * h(i,j,1) * h(i,j,2) / (h(i,j,1)+h(i,j,2))) * ICoriolis
+            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 = Rd / sqrt(G%areaT(i,j)) / (D + Rd / sqrt(G%areaT(i,j)) /CS%ZB_max_memory)
+          endif
+
+      endif
 
       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 +826,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