Fortran: Module Interface kpp: the KPP-turbulence model (Source File: kpp.F90)

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 (43). It assumes that the turbulent flux is the sum of a down-gradient flux and a non-local contribution,

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

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

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

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

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

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

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

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

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

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, (214) can be re-expressed as

$$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]$$ (215)

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 (212), if we require that

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

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 (212), and setting it equal to the $z$-derivative of $\nu_{ti}^\phi$. For the surface layer this results in

$$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$$ (217)

where we used the relation

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

if the motion of the free surface is ignored.

The derivative of $w_\phi$ appearing in (217) can be evaluted with the help of the formulae given in section 4.8.6. As discussed in section 4.8.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 (225) is given by (227), and it is easy to show that

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

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 (219), one can re-write (217) as

$$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$$ (220)

valid only for the surface boundary layer. For the bottom boundary layer, the minus sign in (218) disappears, with the consequence that the minus sign in (220) 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 (211) 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,

$$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$$ (221)

defined by Large et al. (1994), reaches the critical value $Ri_c$. The subscript "r" in (221) 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 (221) by linear interpolation.

To check the boundary layer limit according to the condition

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

two methods have been implemented in GOTM. The first method simply evaluates (222) with a linear interpolation scheme. The second method is activated if the pre-processor macro ${\tt KPP\_IP\_FC}$ is defined. Then, the condition (222) is reformulated as

$$F_c(z_{bl}) = Ri_ (z_{bl}) - Ri_c Ri_ (z_{bl}) = 0 \quad .$$ (223)

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
    $Log: kpp.F90,v $
    Revision 1.6  2011-04-05 14:53:31  jorn
    removed reference to eqstate from kpp
    Revision 1.5  2010-09-17 12:53:52  jorn
    extensive code clean-up to ensure proper initialization and clean-up of all variables
    Revision 1.4  2007-01-06 11:49:15  kbk
    namelist file extension changed .inp --> .nml
    Revision 1.3  2005/11/15 11:35:02  lars
    documentation finish for print
    Revision 1.2  2005/07/21 10:20:00  lars
    polished documentation
    Revision 1.1  2005/06/27 10:54:33  kbk
    new files needed

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