+Made the von Karman constant a runtime parameter

  Made the von Karman constant into a runtime parameter, specified with
VON_KARMAN_CONST or VON_KARMAN_BBL, replacing a hard-coded value in multiple
places.  By default, all answers are bitwise identical, but there are one or
two new entries in the MOM_parameter_doc files.

src/parameterizations/lateral/MOM_mixed_layer_restrat.F90

@@ -48,6 +48,7 @@ module MOM_mixed_layer_restrat
   logical :: MLE_use_PBL_MLD       !< If true, use the MLD provided by the PBL parameterization.
                                    !! if false, MLE will calculate a MLD based on a density difference
                                    !! based on the parameter MLE_DENSITY_DIFF.
+  real    :: vonKar                !< The von Karman constant as used for mixed layer viscosity [nomdim]
   real    :: MLE_MLD_decay_time    !< Time-scale to use in a running-mean when MLD is retreating [T ~> s].
   real    :: MLE_MLD_decay_time2   !< Time-scale to use in a running-mean when filtered MLD is retreating [T ~> s].
   real    :: MLE_density_diff      !< Density difference used in detecting mixed-layer depth [R ~> kg m-3].
@@ -189,6 +190,8 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
   real :: dh        ! Portion of the layer thickness that is in the mixed layer [H ~> m or kg m-2]
   real :: res_scaling_fac ! The resolution-dependent scaling factor [nondim]
   real :: I_LFront ! The inverse of the frontal length scale [L-1 ~> m-1]
+  real :: vonKar_x_pi2    ! A scaling constant that is approximately the von Karman constant times
+                          ! pi squared [nondim]
   logical :: line_is_empty, keep_going, res_upscale
   integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
   integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
@@ -200,6 +203,8 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
   covTS(:)=0.0 !!Functionality not implemented yet; in future, should be passed in tv
   varS(:)=0.0
 
+  vonKar_x_pi2 = CS%vonKar * 9.8696
+
   if (.not.associated(tv%eqn_of_state)) call MOM_error(FATAL, "MOM_mixedlayer_restrat: "// &
          "An equation of state must be used with this module.")
   if (.not. allocated(VarMix%Rd_dx_h) .and. CS%front_length > 0.) &
@@ -380,11 +385,10 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
           ( sqrt( 0.5 * ( G%dxCu(I,j)**2 + G%dyCu(I,j)**2 ) ) * I_LFront ) &
           * min( 1., 0.5*( VarMix%Rd_dx_h(i,j) + VarMix%Rd_dx_h(i+1,j) ) )
 
-    ! peak ML visc: u_star * 0.41 * (h_ml*u_star)/(absf*h_ml + 4.0*u_star)
+    ! peak ML visc: u_star * von_Karman * (h_ml*u_star)/(absf*h_ml + 4.0*u_star)
     ! momentum mixing rate: pi^2*visc/h_ml^2
-    ! 0.41 is the von Karmen constant, 9.8696 = pi^2.
     h_vel = 0.5*((htot_fast(i,j) + htot_fast(i+1,j)) + h_neglect) * GV%H_to_Z
-    mom_mixrate = (0.41*9.8696)*u_star**2 / &
+    mom_mixrate = vonKar_x_pi2*u_star**2 / &
                   (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)
     timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
     timescale = timescale * CS%ml_restrat_coef
@@ -393,7 +397,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
         (Rml_av_fast(i+1,j)-Rml_av_fast(i,j)) * (h_vel**2 * GV%Z_to_H)
     ! As above but using the slow filtered MLD
     h_vel = 0.5*((htot_slow(i,j) + htot_slow(i+1,j)) + h_neglect) * GV%H_to_Z
-    mom_mixrate = (0.41*9.8696)*u_star**2 / &
+    mom_mixrate = vonKar_x_pi2*u_star**2 / &
                   (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)
     timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
     timescale = timescale * CS%ml_restrat_coef2
@@ -456,11 +460,10 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
           ( sqrt( 0.5 * ( (G%dxCv(i,J))**2 + (G%dyCv(i,J))**2 ) ) * I_LFront ) &
           * min( 1., 0.5*( VarMix%Rd_dx_h(i,j) + VarMix%Rd_dx_h(i,j+1) ) )
 
-    ! peak ML visc: u_star * 0.41 * (h_ml*u_star)/(absf*h_ml + 4.0*u_star)
+    ! peak ML visc: u_star * von_Karman * (h_ml*u_star)/(absf*h_ml + 4.0*u_star)
     ! momentum mixing rate: pi^2*visc/h_ml^2
-    ! 0.41 is the von Karmen constant, 9.8696 = pi^2.
     h_vel = 0.5*((htot_fast(i,j) + htot_fast(i,j+1)) + h_neglect) * GV%H_to_Z
-    mom_mixrate = (0.41*9.8696)*u_star**2 / &
+    mom_mixrate = vonKar_x_pi2*u_star**2 / &
                   (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)
     timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
     timescale = timescale * CS%ml_restrat_coef
@@ -469,7 +472,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
         (Rml_av_fast(i,j+1)-Rml_av_fast(i,j)) * (h_vel**2 * GV%Z_to_H)
     ! As above but using the slow filtered MLD
     h_vel = 0.5*((htot_slow(i,j) + htot_slow(i,j+1)) + h_neglect) * GV%H_to_Z
-    mom_mixrate = (0.41*9.8696)*u_star**2 / &
+    mom_mixrate = vonKar_x_pi2*u_star**2 / &
                   (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)
     timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
     timescale = timescale * CS%ml_restrat_coef2
@@ -617,6 +620,8 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
   real :: h_vel           ! htot interpolated onto velocity points [Z ~> m]. (The units are not H.)
   real :: absf            ! absolute value of f, interpolated to velocity points [T-1 ~> s-1]
   real :: u_star          ! surface friction velocity, interpolated to velocity points [Z T-1 ~> m s-1].
+  real :: vonKar_x_pi2    ! A scaling constant that is approximately the von Karman constant times
+                          ! pi squared [nondim]
   real :: mom_mixrate     ! rate at which momentum is homogenized within mixed layer [T-1 ~> s-1]
   real :: timescale       ! mixing growth timescale [T ~> s]
   real :: h_min           ! The minimum layer thickness [H ~> m or kg m-2].  h_min could be 0.
@@ -653,6 +658,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
   uDml(:)    = 0.0 ; vDml(:) = 0.0
   I4dt       = 0.25 / dt
   g_Rho0     = GV%g_Earth / GV%Rho0
+  vonKar_x_pi2 = CS%vonKar * 9.8696
   use_EOS    = associated(tv%eqn_of_state)
   h_neglect  = GV%H_subroundoff
   dz_neglect = GV%H_subroundoff*GV%H_to_Z
@@ -701,10 +707,9 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
 
     u_star = 0.5*(forces%ustar(i,j) + forces%ustar(i+1,j))
     absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J)))
-    ! peak ML visc: u_star * 0.41 * (h_ml*u_star)/(absf*h_ml + 4.0*u_star)
+    ! peak ML visc: u_star * von_Karman * (h_ml*u_star)/(absf*h_ml + 4.0*u_star)
     ! momentum mixing rate: pi^2*visc/h_ml^2
-    ! 0.41 is the von Karmen constant, 9.8696 = pi^2.
-    mom_mixrate = (0.41*9.8696)*u_star**2 / &
+    mom_mixrate = vonKar_x_pi2*u_star**2 / &
                   (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)
     timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
 
@@ -748,10 +753,9 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
 
     u_star = 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1))
     absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J)))
-    ! peak ML visc: u_star * 0.41 * (h_ml*u_star)/(absf*h_ml + 4.0*u_star)
+    ! peak ML visc: u_star * von_Karman * (h_ml*u_star)/(absf*h_ml + 4.0*u_star)
     ! momentum mixing rate: pi^2*visc/h_ml^2
-    ! 0.41 is the von Karmen constant, 9.8696 = pi^2.
-    mom_mixrate = (0.41*9.8696)*u_star**2 / &
+    mom_mixrate = vonKar_x_pi2*u_star**2 / &
                   (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)
     timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
 
@@ -877,6 +881,9 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS,
   call get_param(param_file, mdl, "USE_STANLEY_ML", CS%use_stanley_ml, &
                  "If true, turn on Stanley SGS T variance parameterization "// &
                  "in ML restrat code.", default=.false.)
+  call get_param(param_file, mdl, 'VON_KARMAN_CONST', CS%vonKar, &
+                 'The value the von Karman constant as used for mixed layer viscosity.', &
+                 units='nondim', default=0.41)
   ! We use GV%nkml to distinguish between the old and new implementation of MLE.
   ! The old implementation only works for the layer model with nkml>0.
   if (GV%nkml==0) then

src/parameterizations/vertical/MOM_bulk_mixed_layer.F90

@@ -49,6 +49,7 @@ module MOM_bulk_mixed_layer
                              !! the mixed layer is converted to TKE [nondim].
   real    :: bulk_Ri_convective !< The efficiency with which convectively
                              !! released mean kinetic energy becomes TKE [nondim].
+  real    :: vonKar          !< The von Karman constant as used for mixed layer viscosity [nomdim]
   real    :: Hmix_min        !< The minimum mixed layer thickness [H ~> m or kg m-2].
   real    :: H_limit_fluxes  !< When the total ocean depth is less than this
                              !! value [H ~> m or kg m-2], scale away all surface forcing to
@@ -316,7 +317,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C
   real :: H_nbr ! A minimum thickness based on neighboring thicknesses [H ~> m or kg m-2].
 
   real :: absf_x_H  ! The absolute value of f times the mixed layer thickness [Z T-1 ~> m s-1].
-  real :: kU_star   ! Ustar times the Von Karmen constant [Z T-1 ~> m s-1].
+  real :: kU_star   ! Ustar times the Von Karman constant [Z T-1 ~> m s-1].
   real :: dt__diag  ! A recaled copy of dt_diag (if present) or dt [T ~> s].
   logical :: write_diags  ! If true, write out diagnostics with this step.
   logical :: reset_diags  ! If true, zero out the accumulated diagnostics.
@@ -618,12 +619,12 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C
       ! as the third piece will then optimally describe mixed layer
       ! restratification.  For nkml>=4 the whole strategy should be revisited.
       do i=is,ie
-        kU_star = 0.41*fluxes%ustar(i,j) ! Maybe could be replaced with u*+w*?
+        kU_star = CS%vonKar*fluxes%ustar(i,j) ! Maybe could be replaced with u*+w*?
         if (associated(fluxes%ustar_shelf) .and. &
             associated(fluxes%frac_shelf_h)) then
           if (fluxes%frac_shelf_h(i,j) > 0.0) &
             kU_star = (1.0 - fluxes%frac_shelf_h(i,j)) * kU_star + &
-                      fluxes%frac_shelf_h(i,j) * (0.41*fluxes%ustar_shelf(i,j))
+                      fluxes%frac_shelf_h(i,j) * (CS%vonKar*fluxes%ustar_shelf(i,j))
         endif
         absf_x_H = 0.25 * GV%H_to_Z * h(i,0) * &
             ((abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I-1,J-1))) + &
@@ -1344,11 +1345,11 @@ subroutine find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA,
 !  the equatorial areas.  Although it is not cast as a parameter, it should
 !  be considered an empirical parameter, and it might depend strongly on the
 !  number of sublayers in the mixed layer and their locations.
-!  The 0.41 is VonKarman's constant.  This equation assumes that small & large
-!  scales contribute to mixed layer deepening at similar rates, even though
-!  small scales are dissipated more rapidly (implying they are less efficient).
-!     Ih = 1.0/(16.0*0.41*U_star*dt)
-    Ih = GV%H_to_Z/(3.0*0.41*U_star*dt)
+!    This equation assumes that small & large scales contribute to mixed layer
+!  deepening at similar rates, even though small scales are dissipated more
+!  rapidly (implying they are less efficient).
+!     Ih = 1.0/(16.0*CS%vonKar*U_star*dt)
+    Ih = GV%H_to_Z/(3.0*CS%vonKar*U_star*dt)
     cMKE(1,i) = 4.0 * Ih ; cMKE(2,i) = (absf_Ustar*GV%H_to_Z) * Ih
 
     if (Idecay_len_TKE(i) > 0.0) then
@@ -3387,6 +3388,9 @@ subroutine bulkmixedlayer_init(Time, G, GV, US, param_file, diag, CS)
                  "kinetic energy is converted to turbulent kinetic "//&
                  "energy.  By default BULK_RI_CONVECTIVE=BULK_RI_ML.", &
                  units="nondim", default=CS%bulk_Ri_ML)
+  call get_param(param_file, mdl, 'VON_KARMAN_CONST', CS%vonKar, &
+                 'The value the von Karman constant as used for mixed layer viscosity.', &
+                 units='nondim', default=0.41)
   call get_param(param_file, mdl, "HMIX_MIN", CS%Hmix_min, &
                  "The minimum mixed layer depth if the mixed layer depth "//&
                  "is determined dynamically.", units="m", default=0.0, scale=GV%m_to_H, &

src/parameterizations/vertical/MOM_diabatic_aux.F90

@@ -153,8 +153,7 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo)
         pressure(i,1) = ps(i) + (0.5*H_to_RL2_T2)*h(i,j,1)
       enddo
       do k=2,nz ; do i=is,ie
-        pressure(i,k) = pressure(i,k-1) + &
-          (0.5*H_to_RL2_T2) * (h(i,j,k) + h(i,j,k-1))
+        pressure(i,k) = pressure(i,k-1) + (0.5*H_to_RL2_T2) * (h(i,j,k) + h(i,j,k-1))
       enddo ; enddo
     endif
 

src/parameterizations/vertical/MOM_diapyc_energy_req.F90

@@ -27,8 +27,9 @@ module MOM_diapyc_energy_req
 type, public :: diapyc_energy_req_CS ; private
   logical :: initialized = .false. !< A variable that is here because empty
                                !! structures are not permitted by some compilers.
-  real :: test_Kh_scaling      !< A scaling factor for the diapycnal diffusivity.
-  real :: ColHt_scaling        !< A scaling factor for the column height change correction term.
+  real :: test_Kh_scaling      !< A scaling factor for the diapycnal diffusivity [nondim]
+  real :: ColHt_scaling        !< A scaling factor for the column height change correction term [nondim]
+  real :: VonKar               !< The von Karman coefficient as used in this module [nondim]
   logical :: use_test_Kh_profile !< If true, use the internal test diffusivity profile in place of
                                !! any that might be passed in as an argument.
   type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
@@ -104,7 +105,7 @@ subroutine diapyc_energy_req_test(h_3d, dt, tv, G, GV, US, CS, Kd_int)
       do K=2,nz
         tmp1 = h_top(K) * h_bot(K) * GV%H_to_Z
         Kd(K) = CS%test_Kh_scaling *  &
-                ustar * 0.41 * (tmp1*ustar) / (absf*tmp1 + htot*ustar)
+                ustar * CS%VonKar * (tmp1*ustar) / (absf*tmp1 + htot*ustar)
       enddo
     endif
     may_print = is_root_PE() .and. (i==ie) .and. (j==je)
@@ -1292,10 +1293,12 @@ subroutine diapyc_energy_req_init(Time, G, GV, US, param_file, diag, CS)
   call get_param(param_file, mdl, "ENERGY_REQ_COL_HT_SCALING", CS%ColHt_scaling, &
                  "A scaling factor for the column height change correction "//&
                  "used in testing the energy requirements.", default=1.0, units="nondim")
-  call get_param(param_file, mdl, "ENERGY_REQ_USE_TEST_PROFILE", &
-                 CS%use_test_Kh_profile, &
+  call get_param(param_file, mdl, "ENERGY_REQ_USE_TEST_PROFILE", CS%use_test_Kh_profile, &
                  "If true, use the internal test diffusivity profile in "//&
                  "place of any that might be passed in as an argument.", default=.false.)
+  call get_param(param_file, mdl, 'VON_KARMAN_CONST', CS%vonKar, &
+                 'The value the von Karman constant as used for mixed layer viscosity.', &
+                 units='nondim', default=0.41)
 
   CS%id_ERt = register_diag_field('ocean_model', 'EnReqTest_ERt', diag%axesZi, Time, &
                  "Diffusivity Energy Requirements, top-down", &

src/parameterizations/vertical/MOM_energetic_PBL.F90

@@ -36,8 +36,7 @@ module MOM_energetic_PBL
   logical :: initialized = .false. !< True if this control structure has been initialized.
 
   !/ Constants
-  real    :: VonKar = 0.41   !< The von Karman coefficient.  This should be a runtime parameter,
-                             !! but because it is set to 0.4 at runtime in KPP it might change answers.
+  real    :: VonKar          !< The von Karman coefficient as used in the ePBL module [nondim]
   real    :: omega           !< The Earth's rotation rate [T-1 ~> s-1].
   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
@@ -1982,6 +1981,9 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS)
                  "A nondimensional scaling factor controlling the inhibition "//&
                  "of the diffusive length scale by rotation. Making this larger "//&
                  "decreases the PBL diffusivity.", units="nondim", default=1.0)
+  call get_param(param_file, mdl, 'VON_KARMAN_CONST', CS%vonKar, &
+                 'The value the von Karman constant as used for mixed layer viscosity.', &
+                 units='nondim', default=0.41)
   call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
                  "This sets the default value for the various _ANSWER_DATE parameters.", &
                  default=99991231)

src/parameterizations/vertical/MOM_set_diffusivity.F90

@@ -73,6 +73,8 @@ module MOM_set_diffusivity
                              !! added.
   logical :: use_LOTW_BBL_diffusivity !< If true, use simpler/less precise, BBL diffusivity.
   logical :: LOTW_BBL_use_omega !< If true, use simpler/less precise, BBL diffusivity.
+  real    :: Von_Karm        !< The von Karman constant as used in the BBL diffusivity calculation
+                             !! [nondim].  See (http://en.wikipedia.org/wiki/Von_Karman_constant)
   real    :: BBL_effic       !< efficiency with which the energy extracted
                              !! by bottom drag drives BBL diffusion [nondim]
   real    :: cdrag           !< quadratic drag coefficient [nondim]
@@ -1406,7 +1408,6 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Kd_int
   logical :: Rayleigh_drag ! Set to true if there are Rayleigh drag velocities defined in visc, on
                            ! the assumption that this extracted energy also drives diapycnal mixing.
   integer :: i, k, km1
-  real, parameter :: von_karm = 0.41 ! Von Karman constant (http://en.wikipedia.org/wiki/Von_Karman_constant)
   logical :: do_diag_Kd_BBL
 
   if (.not.(CS%bottomdraglaw .and. (CS%BBL_effic > 0.0))) return
@@ -1483,7 +1484,7 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Kd_int
       if ( ustar_D + absf * ( z_bot * D_minus_z ) == 0.) then
         Kd_wall = 0.
       else
-        Kd_wall = ((von_karm * ustar2) * (z_bot * D_minus_z)) &
+        Kd_wall = ((CS%von_karm * ustar2) * (z_bot * D_minus_z)) &
                   / (ustar_D + absf * (z_bot * D_minus_z))
       endif
 
@@ -1995,6 +1996,7 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_
   ! This include declares and sets the variable "version".
 # include "version_variable.h"
   character(len=40)  :: mdl = "MOM_set_diffusivity"  ! This module's name.
+  real    :: vonKar          ! The von Karman constant as used for mixed layer viscosity [nomdim]
   real    :: omega_frac_dflt ! The default value for the fraction of the absolute rotation rate
                              ! that is used in place of the absolute value of the local Coriolis
                              ! parameter in the denominator of some expressions [nondim]
@@ -2152,6 +2154,12 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_
       call get_param(param_file, mdl, "LOTW_BBL_USE_OMEGA", CS%LOTW_BBL_use_omega, &
                  "If true, use the maximum of Omega and N for the TKE to diffusion "//&
                  "calculation. Otherwise, N is N.", default=.true.)
+      call get_param(param_file, mdl, 'VON_KARMAN_CONST', vonKar, &
+                 'The value the von Karman constant as used for mixed layer viscosity.', &
+                 units='nondim', default=0.41)
+      call get_param(param_file, mdl, 'VON_KARMAN_BBL', CS%von_Karm, &
+                 'The value the von Karman constant as used in calculating the BBL diffusivity.', &
+                 units='nondim', default=vonKar)
     endif
   else
     CS%use_LOTW_BBL_diffusivity = .false. ! This parameterization depends on a u* from viscous BBL