Replaced **0.5 to sqrt() in calc_sqg_struct

Author: Wendazhang33

src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90

@@ -561,8 +561,6 @@ subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE)
    real,                                      intent(in)    :: dt !< Time increment [T ~> s]
   type(VarMix_CS),                           intent(inout) :: CS !< Variable mixing control struct
   type(MEKE_type),                           intent(in)     :: MEKE !< MEKE struct
-  !  type(ocean_OBC_type),                      pointer       :: OBC !< Open
-!  boundaries control structure.
 
   ! Local variables
   real, dimension(SZI_(G), SZJ_(G),SZK_(GV)+1) :: &
@@ -576,7 +574,6 @@ subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE)
   real, dimension(SZI_(G), SZJ_(G)) :: f  ! Absolute value of the Coriolis parameter at h point [T-1 ~> s-1]
   real :: N2  ! Positive buoyancy frequency square or zero [L2 Z-2 T-2 ~> s-2]
   real :: dzc  ! Spacing between two adjacent layers in stretched vertical coordinate [m]
-!  real, dimension(SZK_(GV)) :: zs  ! Stretched vertical coordinate [m]
   real, dimension(SZI_(G), SZJ_(G)) :: Le  ! Eddy length scale [m]
   integer :: i, j, k, is, ie, js, je, nz
 
@@ -596,28 +593,23 @@ subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE)
   if (allocated(MEKE%Le)) then
     do j=js,je ; do i=is,ie
       Le(i,j) = MEKE%Le(i,j)
-      f(i,j) = max(0.25*abs(G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) + &
-                       G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1)), 1.0e-8/US%T_to_s)
+      f(i,j) = max(0.25 * abs((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + &
+                       (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))), 1.0e-8 * US%s_to_T)
     enddo ; enddo
   else
     do j=js,je ; do i=is,ie
       Le(i,j) = sqrt(G%areaT(i,j))
-      f(i,j) = max(0.25*abs(G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) + &
-                       G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1)), 1.0e-8*US%T_to_s)
+      f(i,j) = max(0.25 * abs((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + &
+                       (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))), 1.0e-8 * US%s_to_T)
     enddo ; enddo
   endif
   do k=2,nz ; do j=js,je ; do i=is,ie
-    N2 = max(0.25*(N2_u(I-1,j,k) + N2_u(I,j,k) + N2_v(i,J-1,k) + N2_v(i,J,k)),0.0)
-    dzc = 0.25*(dzu(I-1,j,k) + dzu(I,j,k) + dzv(i,J-1,k) + dzv(i,J,k)) * &
-            N2**0.5/f(i,j)
-!    dzs = -N2**0.5/f(i,j)*dzc
-    CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1)*exp(-CS%sqg_expo*(dzc/Le(i,j)))
+    N2 = max(0.25 * ((N2_u(I-1,j,k) + N2_u(I,j,k)) + (N2_v(i,J-1,k) + N2_v(i,J,k))), 0.0)
+    dzc = 0.25 * ((dzu(I-1,j,k) + dzu(I,j,k)) + (dzv(i,J-1,k) + dzv(i,J,k))) 
+    CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1) * &
+            exp(-CS%sqg_expo * (dzc * sqrt(N2)/(f(i,j) * Le(i,j))))
   enddo ; enddo ; enddo
 
-  if (CS%debug) then
-    call hchksum(CS%sqg_struct, 'sqg_struct', G%HI)
-  endif
-
 
   if (query_averaging_enabled(CS%diag)) then
     if (CS%id_sqg_struct > 0) call post_data(CS%id_sqg_struct, CS%sqg_struct, CS%diag)