+Replaced ML_USE_OMEGA with ML_OMEGA_FRAC

  Replaced the logical test to toggle between using the local value of f or
omega in the rotational scaling for the decay of turbulence.  ML_USE_OMEGA=true
is replicated by ML_OMEGA_FRAC=1.0, and ML_USE_OMEGA=false (the default) is the
same as ML_OMEGA_FRAC=0.0 (the default), but now intermediate values can be
used as well.  Due to clever design, all test cases are bitwise identical, but
run-time parameter names and MOM_parameter_doc files have changed.

src/parameterizations/vertical/MOM_bulk_mixed_layer.F90

@@ -118,9 +118,9 @@ module MOM_bulk_mixed_layer
   integer :: ML_presort_nz_conv_adj ! If ML_resort is true, do convective
                              ! adjustment on this many layers (starting from the
                              ! top) before sorting the remaining layers.
-  logical :: use_omega       !   If true, use the absolute rotation rate instead
-                             ! of the vertical component of rotation when
-                             ! setting the decay scale for turbulence.
+  real    :: omega_frac      !   When setting the decay scale for turbulence, use
+                             ! this fraction of the absolute rotation rate blended
+                             ! with the local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).
   logical :: correct_absorption ! If true, the depth at which penetrating
                              ! shortwave radiation is absorbed is corrected by
                              ! moving some of the heating upward in the water
@@ -1368,7 +1368,7 @@ subroutine find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA,
   is = G%isc ; ie = G%iec ; nz = G%ke
   diag_wt = dt * Idt_diag
 
-  if (CS%use_omega) absf = 2.0*CS%omega
+  if (CS%omega_frac >= 1.0) absf = 2.0*CS%omega
   do i=is,ie
     U_Star = fluxes%ustar(i,j)
     if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then
@@ -1378,9 +1378,12 @@ subroutine find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA,
     endif
 
     if (U_Star < CS%ustar_min) U_Star = CS%ustar_min
-    if (.not.CS%use_omega) &
+    if (CS%omega_frac < 1.0) then
       absf = 0.25*((abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I-1,J-1))) + &
                    (abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I-1,J))))
+      if (CS%omega_frac > 0.0) &
+        absf = sqrt(CS%omega_frac*4.0*CS%omega**2 + (1.0-CS%omega_frac)*absf**2)
+    endif
     absf_Ustar = absf / U_Star
     Idecay_len_TKE(i) = (absf_Ustar * CS%TKE_decay) * GV%H_to_m
 
@@ -3355,8 +3358,9 @@ subroutine bulkmixedlayer_init(Time, G, GV, param_file, diag, CS)
 ! This include declares and sets the variable "version".
 #include "version_variable.h"
   character(len=40)  :: mod = "MOM_mixed_layer"  ! This module's name.
+  real :: omega_frac_dflt
   integer :: isd, ied, jsd, jed
-  logical :: use_temperature
+  logical :: use_temperature, use_omega
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
 
   if (associated(CS)) then
@@ -3446,10 +3450,20 @@ subroutine bulkmixedlayer_init(Time, G, GV, param_file, diag, CS)
   call get_param(param_file, mod, "OMEGA",CS%omega, &
                  "The rotation rate of the earth.", units="s-1", &
                  default=7.2921e-5)
-  call get_param(param_file, mod, "ML_USE_OMEGA", CS%use_omega, &
+  call get_param(param_file, mod, "ML_USE_OMEGA", use_omega, &
                  "If true, use the absolute rotation rate instead of the \n"//&
                  "vertical component of rotation when setting the decay \n"//&
-                 "scale for turbulence.", default=.false.)
+                 "scale for turbulence.", default=.false., do_not_log=.true.)
+  omega_frac_dflt = 0.0
+  if (use_omega) then
+    call MOM_error(WARNING, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
+    omega_frac_dflt = 1.0
+  endif
+  call get_param(param_file, mod, "ML_OMEGA_FRAC", CS%omega_frac, &
+                 "When setting the decay scale for turbulence, use this \n"//&
+                 "fraction of the absolute rotation rate blended with the \n"//&
+                 "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
+                 units="nondim", default=omega_frac_dflt)
   call get_param(param_file, mod, "ML_RESORT", CS%ML_resort, &
                  "If true, resort the topmost layers by potential density \n"//&
                  "before the mixed layer calculations.", default=.false.)

src/parameterizations/vertical/MOM_energetic_PBL.F90

@@ -99,9 +99,9 @@ module MOM_energetic_PBL
                              ! problems, in m s-1.  If the value is small enough,
                              ! this should not affect the solution.
   real    :: omega           !   The Earth's rotation rate, in s-1.
-  logical :: use_omega       !   If true, use the absolute rotation rate instead
-                             ! of the vertical component of rotation when
-                             ! setting the decay scale for turbulence.
+  real    :: omega_frac      !   When setting the decay scale for turbulence, use
+                             ! this fraction of the absolute rotation rate blended
+                             ! with the local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).
   real    :: wstar_ustar_coef ! A ratio relating the efficiency with which
                              ! convectively released energy is converted to a
                              ! turbulent velocity, relative to mechanically
@@ -443,10 +443,12 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, &
       endif
       if (U_Star < CS%ustar_min) U_Star = CS%ustar_min
 
-      if (CS%use_omega) then ; absf(i) = 2.0*CS%omega
+      if (CS%omega_frac >= 1.0) then ; absf(i) = 2.0*CS%omega
       else
         absf(i) = 0.25*((abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I-1,J-1))) + &
                      (abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I-1,J))))
+        if (CS%omega_frac > 0.0) &
+          absf = sqrt(CS%omega_frac*4.0*CS%omega**2 + (1.0-CS%omega_frac)*absf**2)
       endif
 
       mech_TKE(i) = (dt*CS%mstar*GV%Rho0)*((U_Star**3))
@@ -1117,8 +1119,9 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS)
 ! This include declares and sets the variable "version".
 #include "version_variable.h"
   character(len=40)  :: mod = "MOM_energetic_PBL"  ! This module's name.
+  real :: omega_frac_dflt
   integer :: isd, ied, jsd, jed
-  logical :: use_temperature
+  logical :: use_temperature, use_omega
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
 
   if (associated(CS)) then
@@ -1157,10 +1160,20 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS)
   call get_param(param_file, mod, "OMEGA",CS%omega, &
                  "The rotation rate of the earth.", units="s-1", &
                  default=7.2921e-5)
-  call get_param(param_file, mod, "ML_USE_OMEGA", CS%use_omega, &
+  call get_param(param_file, mod, "ML_USE_OMEGA", use_omega, &
                  "If true, use the absolute rotation rate instead of the \n"//&
                  "vertical component of rotation when setting the decay \n"//&
-                 "scale for turbulence.", default=.false.)
+                 "scale for turbulence.", default=.false., do_not_log=.true.)
+  omega_frac_dflt = 0.0
+  if (use_omega) then
+    call MOM_error(WARNING, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
+    omega_frac_dflt = 1.0
+  endif
+  call get_param(param_file, mod, "ML_OMEGA_FRAC", CS%omega_frac, &
+                 "When setting the decay scale for turbulence, use this \n"//&
+                 "fraction of the absolute rotation rate blended with the \n"//&
+                 "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
+                 units="nondim", default=omega_frac_dflt)
   call get_param(param_file, mod, "WSTAR_USTAR_COEF", CS%wstar_ustar_coef, &
                  "A ratio relating the efficiency with which convectively\n"//&
                  "released energy is converted to a turbulent velocity,\n"//&

src/parameterizations/vertical/MOM_set_diffusivity.F90

@@ -255,6 +255,9 @@ module MOM_set_diffusivity
   logical :: ML_use_omega     ! If true, use absolute rotation rate instead
                               ! of the vertical component of rotation when
                               ! setting the decay scale for mixed layer turbulence.
+  real    :: ML_omega_frac    !   When setting the decay scale for turbulence, use
+                              ! this fraction of the absolute rotation rate blended
+                              ! with the local value of f, as f^2 ~= (1-of)*f^2 + of*4*omega^2.
   logical :: user_change_diff ! If true, call user-defined code to change diffusivity.
   logical :: useKappaShear    ! If true, use the kappa_shear module to find the
                               ! shear-driven diapycnal diffusivity.
@@ -1791,11 +1794,13 @@ subroutine add_MLrad_diffusivity(h, fluxes, j, G, GV, CS, Kd, TKE_to_Kd, Kd_int)
   do k=1,kml ; do i=is,ie ; h_ml(i) = h_ml(i) + GV%H_to_m*h(i,j,k) ; enddo ; enddo
 
   do i=is,ie ; if (do_i(i)) then
-    if (CS%ML_use_omega) then
+    if (CS%ML_omega_frac >= 1.0) then
       f_sq = 4.0*Omega2
     else
       f_sq = 0.25*((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + &
                    (G%CoriolisBu(I,J-1)**2 + G%CoriolisBu(I-1,J)**2))
+      if (CS%ML_omega_frac > 0.0) &
+        f_sq = CS%ML_omega_frac*4.0*Omega2 + (1.0-CS%ML_omega_frac)*f_sq
     endif
 
     ustar_sq = max(fluxes%ustar(i,j), CS%ustar_min)**2
@@ -2509,13 +2514,14 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp
 
   real :: decay_length, utide, zbot, hamp
   type(vardesc) :: vd
-  logical :: read_tideamp
+  logical :: read_tideamp, ML_use_omega
 ! This include declares and sets the variable "version".
 #include "version_variable.h"
   character(len=40)  :: mod = "MOM_set_diffusivity"  ! This module's name.
   character(len=20)  :: tmpstr
   character(len=200) :: filename, tideamp_file, h2_file, Niku_TKE_input_file
   real :: Niku_scale ! local variable for scaling the Nikurashin TKE flux data
+  real :: omega_frac_dflt
   integer :: i, j, is, ie, js, je
 
   if (associated(CS)) then
@@ -2587,11 +2593,20 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp
     call get_param(param_file, mod, "TKE_DECAY", CS%TKE_decay, &
                  "The ratio of the natural Ekman depth to the TKE decay scale.", &
                  units="nondim", default=2.5)
-    call get_param(param_file, mod, "ML_USE_OMEGA", CS%ML_use_omega, &
+    call get_param(param_file, mod, "ML_USE_OMEGA", ML_use_omega, &
                  "If true, use the absolute rotation rate instead of the \n"//&
                  "vertical component of rotation when setting the decay \n"//&
-                 "scale for turbulence in the mixed layer. \n"//&
-                 "This is only used if ML_RADIATION is true.", default=.false.)
+                 "scale for turbulence.", default=.false., do_not_log=.true.)
+    omega_frac_dflt = 0.0
+    if (ML_use_omega) then
+      call MOM_error(WARNING, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
+      omega_frac_dflt = 1.0
+    endif
+    call get_param(param_file, mod, "ML_OMEGA_FRAC", CS%ML_omega_frac, &
+                   "When setting the decay scale for turbulence, use this \n"//&
+                   "fraction of the absolute rotation rate blended with the \n"//&
+                   "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
+                   units="nondim", default=omega_frac_dflt)
   endif
 
   call get_param(param_file, mod, "BOTTOMDRAGLAW", CS%bottomdraglaw, &

src/parameterizations/vertical/MOM_set_viscosity.F90

@@ -116,9 +116,9 @@ module MOM_set_visc
                             ! this should not affect the solution.
   real    :: TKE_decay      ! The ratio of the natural Ekman depth to the TKE
                             ! decay scale, nondimensional.
-  logical :: use_omega      !   If true, use the absolute rotation rate instead
-                            ! of the vertical component of rotation when
-                            ! setting the decay scale for turbulence.
+  real    :: omega_frac     !   When setting the decay scale for turbulence, use
+                            ! this fraction of the absolute rotation rate blended
+                            ! with the local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).
   logical :: debug          ! If true, write verbose checksums for debugging purposes.
   type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the
                             ! timing of diagnostic output.
@@ -1064,8 +1064,11 @@ subroutine set_viscous_ML(u, v, h, tv, fluxes, visc, dt, G, GV, CS)
           vhtot(I) = 0.25 * dt_Rho0 * ((fluxes%tauy(i,J) + fluxes%tauy(i+1,J-1)) + &
                                        (fluxes%tauy(i,J-1) + fluxes%tauy(i+1,J)))
 
-          if (CS%use_omega) then ; absf = 2.0*CS%omega
-          else ; absf = 0.5*(abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I,J-1))) ; endif
+          if (CS%omega_frac >= 1.0) then ; absf = 2.0*CS%omega ; else
+            absf = 0.5*(abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I,J-1)))
+            if (CS%omega_frac > 0.0) &
+              absf = sqrt(CS%omega_frac*4.0*CS%omega**2 + (1.0-CS%omega_frac)*absf**2)
+          endif
           U_Star = max(CS%ustar_min, 0.5 * (fluxes%ustar(i,j) + fluxes%ustar(i+1,j)))
           Idecay_len_TKE(I) = ((absf / U_Star) * CS%TKE_decay) * H_to_m
         endif
@@ -1304,8 +1307,12 @@ subroutine set_viscous_ML(u, v, h, tv, fluxes, visc, dt, G, GV, CS)
           uhtot(i) = 0.25 * dt_Rho0 * ((fluxes%taux(I,j) + fluxes%tauy(I-1,j+1)) + &
                                        (fluxes%taux(I-1,j) + fluxes%tauy(I,j+1)))
 
-         if (CS%use_omega) then ; absf = 2.0*CS%omega
-         else ; absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J))) ; endif
+         if (CS%omega_frac >= 1.0) then ; absf = 2.0*CS%omega ; else
+           absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J)))
+           if (CS%omega_frac > 0.0) &
+             absf = sqrt(CS%omega_frac*4.0*CS%omega**2 + (1.0-CS%omega_frac)*absf**2)        
+         endif
+
          U_Star = max(CS%ustar_min, 0.5 * (fluxes%ustar(i,j) + fluxes%ustar(i,j+1)))
          Idecay_len_TKE(i) = ((absf / U_Star) * CS%TKE_decay) * H_to_m
 
@@ -1617,8 +1624,9 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS)
 !                 for this module
   real    :: Csmag_chan_dflt, smag_const1, TKE_decay_dflt, bulk_Ri_ML_dflt
   real    :: Kv_background
+  real    :: omega_frac_dflt
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
-  logical :: use_kappa_shear, adiabatic, differential_diffusion
+  logical :: use_kappa_shear, adiabatic, differential_diffusion, use_omega
 ! This include declares and sets the variable "version".
 #include "version_variable.h"
   character(len=40)  :: mod = "MOM_set_visc"  ! This module's name.
@@ -1695,10 +1703,20 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS)
                  "mixed layer viscosity.  By default, \n"//&
                  "TKE_DECAY_VISC = TKE_DECAY or 0.", units="nondim", &
                  default=TKE_decay_dflt)
-    call get_param(param_file, mod, "ML_USE_OMEGA", CS%use_omega, &
+    call get_param(param_file, mod, "ML_USE_OMEGA", use_omega, &
                  "If true, use the absolute rotation rate instead of the \n"//&
                  "vertical component of rotation when setting the decay \n"//&
-                 "scale for turbulence.", default=.false.)
+                   "scale for turbulence.", default=.false., do_not_log=.true.)
+    omega_frac_dflt = 0.0
+    if (use_omega) then
+      call MOM_error(WARNING, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
+      omega_frac_dflt = 1.0
+    endif
+    call get_param(param_file, mod, "ML_OMEGA_FRAC", CS%omega_frac, &
+                   "When setting the decay scale for turbulence, use this \n"//&
+                   "fraction of the absolute rotation rate blended with the \n"//&
+                   "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
+                   units="nondim", default=omega_frac_dflt)
     call get_param(param_file, mod, "OMEGA", CS%omega, &
                  "The rotation rate of the earth.", units="s-1", &
                  default=7.2921e-5)