Improvements from Luke Van Roekel

This branch adds an option to revert a part of the computation of the similarity
similarity functions in the turbulent scales routine to the original Large et al.
(1994) formulation.  In particular this applies to the stable buoyancy
forcing with wind stress case.  In KPP, the similarity function in this
regime is given as

    1 + 5*zeta

where zeta == sigma * OBL_depth / Monin_obukhov_length and sigma ==
depth/OBL_depth

In Large et al. (1994), zeta is allowed to vary from 0 to 1.  This is
mostly an assumption that scalings in the surface layer continue through
the boundary layer.  However, Appendix B of Large et al. (1994)
reference some observational support for this.

In Danabasoglu et al. (2006), Appendix A, zeta is confined to run
between zero and epsilon, where epsilon == surface layer extent /
OBL_depth (usually taken as 0.1).  This was done to increase the
unresolved velocity shear in the bulk richardson number calculation (see
Equations A1 and A2 of Danabasoglu et al (2006)).

In tests conducted against LES (I am using the NCAR LES) forced by a
constant wind stress and positive buoyancy forcing, the corresponding
SCM result without the limitation on zeta is closer to LES.

In this branch you can set l_LMD_ws in the cvmix_kpp_init routine. If
this is set to true, the limitaton on zeta is removed in stable buoyancy
forcing conditions, assuming the windstress is not zero.  The default
for this flag is false, so doing nothing will result in the current
CVMIX implementation.

src/shared/cvmix_kpp.F90

@@ -174,6 +174,14 @@ module cvmix_kpp
       logical        :: lenhanced_entr ! True => enhance entrainment by adding
                                        !  Stokes shear to the unresolved 
                                        !  vertial shear
+      logical        :: l_LMD_ws       ! flag to use original Large et al. (1994)
+                                       ! equations for computing turbulent scales
+                                       ! rather than the updated methodology in
+                                       ! Danabasoglu et al. (2006). The latter
+                                       ! limits sigma to be < surf_layer_extent
+                                       ! when computing turbulent scales while
+                                       ! the former only imposes this restriction
+                                       ! in unstable regimes.
   end type cvmix_kpp_params_type
 
 !EOP
@@ -192,7 +200,7 @@ subroutine cvmix_init_kpp(ri_crit, minOBLdepth, maxOBLdepth, minVtsqr,      &
                             Cv, interp_type, interp_type2, MatchTechnique,    &
                             old_vals, lEkman, lMonOb, lnoDGat1,               &
                             lenhanced_diff, lnonzero_surf_nonlocal,           &
-                            llangmuirEF, lenhanced_entr,                      &
+                            llangmuirEF, lenhanced_entr, l_LMD_ws,            &
                             CVmix_kpp_params_user)
 
 ! !DESCRIPTION:
@@ -225,7 +233,8 @@ subroutine cvmix_init_kpp(ri_crit, minOBLdepth, maxOBLdepth, minVtsqr,      &
                                               lenhanced_diff,                 &
                                               llangmuirEF,                    &
                                               lenhanced_entr,                 &
-                                              lnonzero_surf_nonlocal
+                                              lnonzero_surf_nonlocal,         &
+                                              l_LMD_ws
 
 ! !OUTPUT PARAMETERS:
     type(cvmix_kpp_params_type), intent(inout), target, optional ::           &
@@ -501,6 +510,15 @@ subroutine cvmix_init_kpp(ri_crit, minOBLdepth, maxOBLdepth, minVtsqr,      &
     end if
     call cvmix_put_kpp('nonlocal_coeff',nonlocal_coeff,CVmix_kpp_params_user)
 
+    ! By default, use sigma construction from Danabasoglu et al. when computing
+    ! turbulent scales. Set l_LMD_ws = .true. to use Large et al. construction.
+    if (present(l_LMD_ws)) then
+      call cvmix_put_kpp('l_LMD_ws', l_LMD_ws,                    &
+                         CVmix_kpp_params_user)
+    else
+      call cvmix_put_kpp('l_LMD_ws', .false., CVmix_kpp_params_user)
+    end if
+
 !EOC
 
   end subroutine cvmix_init_kpp
@@ -1168,6 +1186,8 @@ subroutine cvmix_put_kpp_logical(varname, val, CVmix_kpp_params_user)
         CVmix_kpp_params_out%llangmuirEF = val
       case ('lenhanced_entr')
         CVmix_kpp_params_out%lenhanced_entr = val
+      case ('l_LMD_ws')
+        CVmix_kpp_params_out%l_LMD_ws = val
       case DEFAULT
         print*, "ERROR: ", trim(varname), " is not a boolean variable!"
         stop 1
@@ -1864,8 +1884,8 @@ subroutine cvmix_kpp_compute_turbulent_scales_1d_sigma(sigma_coord,         &
 
     ! Local variables
     integer :: n_sigma, kw
-    logical :: compute_wm, compute_ws
-    real(cvmix_r8), dimension(size(sigma_coord)) :: zeta
+    logical :: compute_wm, compute_ws, l_LMD_ws
+    real(cvmix_r8), dimension(size(sigma_coord)) :: zeta, sigma_loc
     real(cvmix_r8) :: vonkar, surf_layer_ext
     type(cvmix_kpp_params_type), pointer :: CVmix_kpp_params_in
 
@@ -1878,16 +1898,21 @@ subroutine cvmix_kpp_compute_turbulent_scales_1d_sigma(sigma_coord,         &
 
     compute_wm = present(w_m)
     compute_ws = present(w_s)
-    vonkar = cvmix_get_kpp_real('vonkarman', CVmix_kpp_params_in)
-    surf_layer_ext = cvmix_get_kpp_real('surf_layer_ext', CVmix_kpp_params_in)
+
+    l_LMD_ws       = CVmix_kpp_params_in%l_LMD_ws
+    vonkar         = CVmix_kpp_params_in%vonkarman
+    surf_layer_ext = CVmix_kpp_params_in%surf_layer_ext
 
     if (surf_fric_vel.ne.cvmix_zero) then
-      do kw=1,n_sigma
-        ! compute scales at sigma if sigma < surf_layer_ext, otherwise compute
-        ! at surf_layer_ext
-        zeta(kw) = min(surf_layer_ext, sigma_coord(kw)) * OBL_depth *         &
-                   surf_buoy_force*vonkar/(surf_fric_vel**3)
-      end do
+      if ((surf_buoy_force.ge.cvmix_zero) .and. l_LMD_ws) then
+        sigma_loc(:) = sigma_coord(:)
+      else
+        sigma_loc(:) = min(surf_layer_ext, sigma_coord(:))
+      end if
+      ! compute scales at sigma if sigma < surf_layer_ext, otherwise compute
+      ! at surf_layer_ext
+      zeta(:) = sigma_loc(:) * OBL_depth * surf_buoy_force * vonkar /         &
+                (surf_fric_vel**3)
 
       if (compute_wm) then
         w_m(1) = compute_phi_inv(zeta(1), CVmix_kpp_params_in, lphi_m=.true.)*&
@@ -2029,8 +2054,8 @@ subroutine cvmix_kpp_compute_turbulent_scales_1d_OBL(sigma_coord,           &
 
     ! Local variables
     integer :: n_sigma, kw
-    logical :: compute_wm, compute_ws
-    real(cvmix_r8), dimension(size(surf_buoy_force)) :: zeta
+    logical :: compute_wm, compute_ws, l_LMD_ws
+    real(cvmix_r8), dimension(size(surf_buoy_force)) :: zeta, sigma_loc
     real(cvmix_r8) :: vonkar, surf_layer_ext
     type(cvmix_kpp_params_type), pointer :: CVmix_kpp_params_in
 
@@ -2043,16 +2068,20 @@ subroutine cvmix_kpp_compute_turbulent_scales_1d_OBL(sigma_coord,           &
 
     compute_wm = present(w_m)
     compute_ws = present(w_s)
-    vonkar = cvmix_get_kpp_real('vonkarman', CVmix_kpp_params_in)
-    surf_layer_ext = cvmix_get_kpp_real('surf_layer_ext', CVmix_kpp_params_in)
+
+    l_LMD_ws       = CVmix_kpp_params_in%l_LMD_ws
+    vonkar         = CVmix_kpp_params_in%vonkarman
+    surf_layer_ext = CVmix_kpp_params_in%surf_layer_ext
 
     if (surf_fric_vel.ne.cvmix_zero) then
-      do kw=1,n_sigma
-        ! compute scales at sigma if sigma < surf_layer_ext, otherwise compute
-        ! at surf_layer_ext
-        zeta(kw) = min(surf_layer_ext, sigma_coord) * OBL_depth(kw) *         &
-                   surf_buoy_force(kw)*vonkar/(surf_fric_vel**3)
-      end do
+      sigma_loc = min(surf_layer_ext, sigma_coord)
+      if (l_LMD_ws) then
+        where (surf_buoy_force.ge.cvmix_zero)
+          sigma_loc = sigma_coord
+        end where
+      end if
+      zeta(:) = sigma_loc(:) * OBL_depth(:) * surf_buoy_force(:) * vonkar /   &
+                (surf_fric_vel**3)
 
       if (compute_wm) then
         w_m(1) = compute_phi_inv(zeta(1), CVmix_kpp_params_in, lphi_m=.true.)*&