Module kpp: the KPP-turbulence model


INTERFACE:

   module kpp
DESCRIPTION:

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,

$\displaystyle {\langle w'\phi' \rangle} = - \nu_t^\phi \dfrac{\partial { {\langle \phi \rangle} }}{\partial {z}} + \tilde{\Gamma}_\phi \; , \quad$ (212)

where the super- or subscript $ \phi$ is a placeholder for the symbols $ m$, $ h$, and $ s$, 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 KPP_SALINITY has been defined in cppdefs.h.

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

$\displaystyle \tilde{\Gamma}_\phi = \nu_t^\phi \gamma_\phi \; , \quad$ (213)

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

$\displaystyle \nu_t^\phi = h w_\phi(\sigma) G(\sigma) \; , \quad$ (214)

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

$\displaystyle \sigma = \dfrac{d}{h} \; , \quad$ (215)

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

$\displaystyle G(\sigma) = a_0 + a_1 \sigma + a_2 \sigma^2 + a_3 \sigma^3 \quad .$ (216)

Physical arguments discussed in Large et al. (1994) require $ a_0=0$, $ a_1=1$. The remaining two parameters $ a_2$ and $ a_3$ may be re-expressed in terms of the value of $ G$ and its derivative, $ G'$, at the edge of the boundary layer, $ \sigma=1$. Then, (216) can be re-expressed as

$\displaystyle G(\sigma) = \sigma \Big[ 1 + \sigma \Big( \left( \sigma - 2 \righ...
... \left( 3 - 2 \sigma \right) G(1) + \left( \sigma - 1 \right) G'(1) \Big) \Big]$ (217)

Apart from the boundary layer diffusivities, the KPP model also computes "interior" diffusivities, here denoted by $ \nu_{ti}^{\phi}$. The function $ G$ 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

$\displaystyle G(1) = \dfrac{1}{h w_\phi(1)} \nu_{ti}^\phi(z_{bl}) \; , \quad$ (218)

where $ z_{bl}$ denotes the vertical coordinate of the surface (or bottom) boundary layer.

A condition for the continuity of the derivatives of $ \nu_t^\phi$ and $ \nu_{ti}^\phi$ can be obtained by carrying out the derivative with respect to $ z$ of (214), and setting it equal to the $ z$-derivative of $ \nu_{ti}^\phi$. For the surface layer this results in

$\displaystyle G'(1) = - \dfrac{G(1)}{w(1)} \dfrac{\partial {w}}{\partial {\sigm...
...\dfrac{\partial {\nu_{ti}^\phi}}{\partial {z}} \Big\vert _{z=z_{bl}} \; , \quad$ (219)

where we used the relation

$\displaystyle \dfrac{\partial {}}{\partial {z}} = - \frac{1}{h} \dfrac{\partial {}}{\partial {\sigma}} \; , \quad$ (220)

if the motion of the free surface is ignored.

The derivative of $ w_\phi$ 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 $ \sigma=1$, the derivative of $ w_\phi$ is different from zero only for stably stratified flows. Then, the non-dimensional function $ \Phi_\phi$ appearing in (227) is given by (229), and it is easy to show that

$\displaystyle \dfrac{\partial {w}}{\partial {\sigma}} \Big\vert _{\sigma=1} = - 5 h w_\phi(1) \dfrac{B_f}{u_*^4} \; , \quad$ (221)

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 $ \kappa$ in this relation.

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

$\displaystyle G'(1) = \dfrac{B_f}{u_*^4} \nu_{ti}^\phi \Big\vert _{z=z_{bl}} - ...
...\dfrac{\partial {\nu_{ti}^\phi}}{\partial {z}} \Big\vert _{z=z_{bl}} \; , \quad$ (222)

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 KPP_CLIP_GS is defined in cppdef.h, the slope of $ G$ is set to zero for negative slopes. For stably stratified flows with a stabilizing buoyancy flux, this limiter breaks the continuity of the first derivatives.

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, $ z_{bl}$, corresponds to the position where the bulk Richardson number,

$\displaystyle Ri_b = \dfrac{(B_r - B(z_{bl})) d} {\vert {\boldsymbol U}_r - {\boldsymbol U}(z_{bl}) \vert^2 + V_t^2(z_{bl})} \; , \quad$ (223)

defined by Large et al. (1994), reaches the critical value $ Ri_c$. 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 $ {\tt KPP\_TWOPOINT\_REF}$ is defined) or by simply identifying the neighbouring center-value with the value at the face. The "turbulent velocity shear", $ V_t$, is computed as described by Large et al. (1994). The value of $ z_{bl}$ is then found from (223) by linear interpolation.

To check the boundary layer limit according to the condition

$\displaystyle Ri_b(z_{bl}) = \dfrac{Ri_\text{top}(z_{bl})}{Ri_\text{bot}(z_{bl})} = Ri_c \; , \quad$ (224)

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 $ {\tt KPP\_IP\_FC}$ is defined. Then, the condition (224) is reformulated as

$\displaystyle F_c(z_{bl}) = Ri_$top$\displaystyle (z_{bl}) - Ri_c Ri_$bot$\displaystyle (z_{bl}) = 0 \quad .$ (225)

The position where the function $ F_c$ 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 MEMBER FUNCTIONS:
  public init_kpp, do_kpp, clean_kpp
PUBLIC DATA MEMBERS:
  z-position of surface boundary layer depth
  REALTYPE, public                 :: zsbl
 
  z-position of bottom boundary layer depth
  REALTYPE, public                 :: zbbl
DEFINED PARAMETERS:
   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
REVISION HISTORY:
   Original author(s): Lars Umlauf
LOCAL VARIABLES:
   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



Subsections
Karsten Bolding 2012-12-28