Option to test Danabasoglu et al 2006 limitation on zeta in the stable buoyancy case with windstorms

[vanroekel]

This branch adds an option to revert a part of the computation of the 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 (see for example cvmix_kpp_compute_turbulent_scales_1d). The default
for this flag is false, so doing nothing will result in the current
CVMIX implementation.

[mnlevy1981]

Hi Luke,

Sorry to take so long to get back to you about this (I think I also owe you an email about testing in POP). I like what you have done, but would make one suggestion -- at line 1858, where you have

      if(CVmix_kpp_params_in%l_LMD_ws) then
        do kw=1,n_sigma
        ! compute scales at sigma if sigma < surf_layer_ext, otherwise compute
        ! at surf_layer_ext
                if(surf_buoy_force .ge. cvmix_zero) then
                        zeta(kw) = sigma_coord(kw) * OBL_depth *                      &
                                surf_buoy_force*vonkar/(surf_fric_vel**3)
                else
                        zeta(kw) = min(surf_layer_ext, sigma_coord(kw)) * OBL_depth *         &
                                surf_buoy_force*vonkar/(surf_fric_vel**3)
                endif
        end do
      else
        do kw=1,n_sigma
                zeta(kw) = min(surf_layer_ext, sigma_coord(kw)) * OBL_depth *         &
                                surf_buoy_force*vonkar/(surf_fric_vel**3)
        end do
      endif

I would prefer to see something like

zeta(kw) = sigma_loc(kw) * OBL_depth *         &
                                surf_buoy_force*vonkar/(surf_fric_vel**3)

where sigma_loc is determined based on logicals... something like

if ((surf_buoy_force .ge. cvmix_zero).and.CVmix_kpp_params_in%l_LMD_ws) then
  sigma_loc(:) = sigma_coord(:)
else
  sigma_loc(:) = min(surf_layer_ext, sigma_coord(:))
end if
zeta(:) = sigma_loc(:) * OBL_depth * surf_buoy_force*vonkar/(surf_fric_vel**3)

That's completely untested code, so I may have some dimension issues, but the main idea is that I don't like seeing zeta set by identical formulas in two different places (lines 1866/67 and lines 1872/73); introducing a local variable for the proper value of sigma to use (either sigma_coord or surface_layer_ext, depending on the logicals) and computing zeta in just a single place seems less error prone if something changes in the future.

A similar comment holds for the block of code at line 2000.

[vanroekel]

Hello Michael,
I think I've addressed what you are requesting. If not, let me know. Just for my own understanding, may I ask why you prefer to compute this sigma_coord_loc twice instead of zeta? Is it because the zeta is used many places in the code?

[mnlevy1981]

Two things:

1. Something weird seems to be going on with this pull request, merging as is looks like it'll wipe the Langmuir mixing options we added to CVMix?
2. You're getting closer to what I was asking about, but I'm trying to avoid duplicate lines of code while 2012 & 2019 still look very similar.

I'll make a new branch this weekend with your changes & my code-style recommendations (without removing existing functionality) and point you to it to make sure it looks okay. At this point I think we want to try to start a clean pull request (I'll leave this one open until the new one is available).

[mnlevy1981]

Closing this in favor of #62 -- @vanroekel can you please try this version of the code and make sure it behaves the way you expect? It has passed my brief tests, so I believe the default behavior has not changed.