Module kpp: the KPP-turbulence model

*INTERFACE:*

module kpp

This implentation of the KPP turbulence parameterisation is based on the publications of Large et al. (1994) and Durksi et al. (2004). The general expression for the turbulent fluxes used in the KPP model is identical to that suggested in (45). It assumes that the turbulent flux is the sum of a down-gradient flux and a non-local contribution,

where the super- or subscript is a placeholder for the symbols , , and , indicating whether a quantity relates to momentum, heat, or salinity (or any other tracer), respectively. Note that turbulence parameters due to salinity stratification are updated only if the pre-processor macro

In the notation of the KPP model, the non-local flux is expressed as

where independent models are valid for and . The KPP model assumes that the turbulent diffusivity, , inside the surface or bottom boundary layer is determined by a relation of the form

where denotes the thickness of the boundary layer, computed according to the algorithm discussed below. The non-dimensional boundary layer coordinate is defined according to

where is the distance from the free surface (or the bottom boundary). The velocity scale, , in (214) is computed as described in section 4.35.6. The dimensionless shape function is a cubic polynomial,

Physical arguments discussed in Large et al. (1994) require , . The remaining two parameters and may be re-expressed in terms of the value of and its derivative, , at the edge of the boundary layer, . Then, (216) can be re-expressed as

Apart from the boundary layer diffusivities, the KPP model also computes "interior" diffusivities, here denoted by . The function and its derivative can be evaluted from the requirement that, at the edge of the boundary layer, the boundary layer diffusivity and its derivative correspond exactly to the interior diffusivity and its derivative, respectively.

Continuity of the boundary and interior diffusivites is obviously insured, see (214), if we require that

where denotes the vertical coordinate of the surface (or bottom) boundary layer.

A condition for the continuity of the derivatives of and can be obtained by carrying out the derivative with respect to of (214), and setting it equal to the -derivative of . For the surface layer this results in

where we used the relation

if the motion of the free surface is ignored.

The derivative of appearing in (219) can be evaluted with the help of the formulae given in section 4.35.6. As discussed in section 4.35.6, at , the derivative of is different from zero only for stably stratified flows. Then, the non-dimensional function appearing in (227) is given by (229), and it is easy to show that

valid for both, bottom and surface boundary layers. Note that in the original publication of Large et al. (1994), erroneously, there appears an additional factor in this relation.

With the help of (221), one can re-write (219) as

valid only for the surface boundary layer. For the bottom boundary layer, the minus sign in (220) disappears, with the consequence that the minus sign in (222) has to be replaced by a plus. Note that if the pre-processor macro

The non-local transport term defined in (213) is computed
as described in Large et al. (1994), if the pre-processor macro
`NONLOCAL` is defined. Otherwise, non-local transport is ignored.

The position of the surface boundary layer depth, , corresponds to the position where the bulk Richardson number,

defined by Large et al. (1994), reaches the critical value . The subscript "r" in (223) denotes a certain reference value of the buoyancy and velocity close to the surface. The choice of this reference value is not unique, and several possibilities have been implemented in numerical models. Presently, GOTM uses the uppermost grid point as the reference value. The bulk Richardson-number is then computed at the grid faces by linear interpolation of quantities defined at the centers (if is defined) or by simply identifying the neighbouring center-value with the value at the face. The "turbulent velocity shear", , is computed as described by Large et al. (1994). The value of is then found from (223) by linear interpolation.

To check the boundary layer limit according to the condition

two methods have been implemented in GOTM. The first method simply evaluates (224) with a linear interpolation scheme. The second method is activated if the pre-processor macro is defined. Then, the condition (224) is reformulated as

The position where the function changes sign is computed from linear interpolation. This method has been suggested in the ROMS code as the numerically more stable alternative. Clearly, all approaches are grid-depending, a difficulty that cannot be overcome with the KPP model.

Finally, provided `clip_mld=.true.` in `kpp.nml`, the boundary layer is cut
if it exceeds the Ekman or the Monin-Obukhov length scale, see Large et al. (1994).

*USES:*

use turbulence, only: num,nuh,nus use turbulence, only: gamu,gamv,gamh,gams use turbulence, only: Rig use turbulence, only: kappa ifdef EXTRA_OUTPUT use turbulence, only: turb1,turb2,turb3,turb4,turb5 endif use eqstate IMPLICIT NONE private

public init_kpp, do_kpp, clean_kpp

z-position of surface boundary layer depth REALTYPE, public :: zsbl z-position of bottom boundary layer depth REALTYPE, public :: zbbl

non-dimensional extent of the surface layer (epsilon=0.1) REALTYPE, parameter :: epsilon = 0.1 critical gradient Richardson number below which turbulent mixing occurs (Ri0=0.7) REALTYPE, parameter :: Ri0 = 0.7 value of double-diffusive density ratio where mixing goes to zero in salt fingering (Rrho0=1.9) REALTYPE, parameter :: Rrho0 = 1.9 buoancy frequency (1/s2) limit for convection (bvfcon=-2.0E-5) REALTYPE, parameter :: bvfcon = -2.0E-5 scaling factor for double diffusion of temperature in salt fingering case (fdd=0.7) REALTYPE, parameter :: fdd = 0.7 maximum interior convective viscosity and diffusivity due to shear instability (nu0c=0.01) REALTYPE, parameter :: nu0c = 0.01 maximum interior viscosity (m2/s) due to shear instability (nu0m=10.0E-4) REALTYPE, parameter :: nu0m = 10.0E-4 maximum interior diffusivity (m2/s) due to shear instability (nu0s=10.0E-4) REALTYPE, parameter :: nu0s = 10.0E-4 scaling factor for double diffusion in salt fingering (nu=1.5E-6) REALTYPE, parameter :: nu = 1.5E-6 scaling factor for double diffusion in salt fingering (nuf=10.0E-4) REALTYPE, parameter :: nuf = 10.0E-4 interior viscosity (m2/s) due to wave breaking (nuwm=1.0E-5) REALTYPE, parameter :: nuwm = 1.0E-5 interior diffusivity (m2/s) due to wave breaking (nuwm=1.0E-6) REALTYPE, parameter :: nuws = 1.0E-6 double diffusion constant for salinity in diffusive convection case (sdd1=0.15) REALTYPE, parameter :: sdd1 = 0.15 double diffusion constant for salinity in diffusive convection (sdd2=1.85) REALTYPE, parameter :: sdd2 = 1.85 double diffusion constant for salinity in diffusive convection (sdd3=0.85) REALTYPE, parameter :: sdd3 = 0.85 double diffusion constant for temperature in diffusive convection (tdd1=0.909) REALTYPE, parameter :: tdd1 = 0.909 double diffusion constant for temperature in diffusive convection (tdd2=4.6) REALTYPE, parameter :: tdd2 = 4.6 double diffusion constant for temperature in diffusive convection case (tdd3=0.54). REALTYPE, parameter :: tdd3 = 0.54 proportionality coefficient parameterizing nonlocal transport (Cstar=10.0) REALTYPE, parameter :: Cstar = 10.0 ratio of interior Brunt-Vaisala frequency to that at entrainment depth (Cv=1.5-1.6) REALTYPE, parameter :: Cv = 1.6 ratio of entrainment flux to surface buoyancy flux (betaT=-0.2) REALTYPE, parameter :: betaT = -0.2 constant for computation of Ekman scale (cekman=0.7) REALTYPE, parameter :: cekman = 0.7 constant for computation of Monin-Obukhov scale (cmonob = 1.0) REALTYPE, parameter :: cmonob = 1.0 coefficient of flux profile for momentum in their 1/3 power law regimes (am=1.26) REALTYPE, parameter :: am = 1.257 coefficient of flux profile for momentum in their 1/3 power law regimes (as=-28.86) REALTYPE, parameter :: as = -28.86 coefficient of flux profile for momentum in their 1/3 power law regimes (cm=8.38) REALTYPE, parameter :: cm = 8.38 coefficient of flux profile for momentum in their 1/3 power law regimes (cs=98.96) REALTYPE, parameter :: cs = 98.96 maximum stability parameter "zeta" value of the 1/3 power law regime of flux profile for momentum (zetam=-0.2) REALTYPE, parameter :: zetam = -0.2 maximum stability parameter "zeta" value of the 1/3 power law regime of flux profile for tracers (zetas=-1.0) REALTYPE, parameter :: zetas = -1.0

Original author(s): Lars Umlauf

proportionality coefficient for parameterizing non-local transport REALTYPE :: Cg coefficient from computation of turbulent shear velocity REALTYPE :: Vtc acceleration of gravity REALTYPE :: g reference density REALTYPE :: rho_0 g/rho_0 REALTYPE :: gorho0 critical bulk Richardson number REALTYPE :: Ric compute surface and bottom BBL logical :: kpp_sbl,kpp_bbl compute internal mixing logical :: kpp_interior use clipping of MLD at Ekman and Monin-Oboukhov scale logical :: clip_mld positions of grid faces and centers REALTYPE, dimension(:), allocatable :: z_w,z_r distance between centers REALTYPE, dimension(:), allocatable :: h_r integer :: ksblOld REALTYPE :: zsblOld