#include"cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !MODULE: kpp: the KPP-turbulence model \label{sec:kpp}
!
! !INTERFACE:
   module kpp
!
! !DESCRIPTION:
! This implentation of the KPP turbulence parameterisation is based on the
! publications of \cite{Largeetal94} and \cite{Durskietal2004}.
! The general expression for the turbulent fluxes used in the KPP model is identical to
! that suggested in \eq{fluxes}. It assumes that the turbulent flux is the sum of a
! down-gradient flux and a non-local contribution,
! \begin{equation}
!   \label{kppFluxes}
!   \mean{w'\phi'} = - \nu_t^\phi  \partder{\mean{\phi}}{z}  + \tilde{\Gamma}_\phi
!  \comma
! \end{equation}
! 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 {\tt KPP\_SALINITY} has been defined in {\tt cppdefs.h}.
!
! In the notation of the KPP model, the non-local flux is expressed as
! \begin{equation}
!   \label{kppNonlocal}
!   \tilde{\Gamma}_\phi = \nu_t^\phi \gamma_\phi
!  \comma
! \end{equation}
! 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
! \begin{equation}
!   \label{Kx}
!   \nu_t^\phi = h w_\phi(\sigma) G(\sigma)
!  \comma
! \end{equation}
! 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
! \begin{equation}
!   \label{defSigma}
!   \sigma = \dfrac{d}{h}
!  \comma
! \end{equation}
! where $d$ is the distance from the free surface (or the bottom boundary).
! The velocity scale, $w_\phi$, in \eq{Kx} is computed as described in \sect{sec:wscale}. The
! dimensionless shape function $G$ is a cubic polynomial,
! \begin{equation}
!   \label{GsigmaA}
!   G(\sigma) = a_0 + a_1 \sigma + a_2 \sigma^2 + a_3 \sigma^3
! \point
! \end{equation}
! Physical arguments discussed in \cite{Largeetal94} 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, \eq{GsigmaA} can be re-expressed as
! \begin{equation}
!   \label{GsigmaB}
!   G(\sigma) = \sigma \Big[
!   1 + \sigma \Big(
!     \left( \sigma - 2 \right)
!   + \left( 3 - 2 \sigma  \right) G(1)
!   + \left( \sigma - 1 \right) G'(1)
!   \Big) \Big]
! \end{equation}
! 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
! \eq{Kx}, if we require that
! \begin{equation}
!   \label{GOfOne}
!   G(1) = \dfrac{1}{h w_\phi(1)} \nu_{ti}^\phi(z_{bl})
!  \comma
! \end{equation}
! 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 \eq{Kx}, and
! setting it equal to the $z$-derivative of $\nu_{ti}^\phi$. For the surface layer
! this results in
! \begin{equation}
!   \label{GPrimeOfOne}
!   G'(1) = - \dfrac{G(1)}{w(1)} \partder{w}{\sigma} \Big|_{\sigma=1}
!           - \dfrac{1}{w(1)} \partder{\nu_{ti}^\phi}{z} \Big|_{z=z_{bl}}
! \comma
! \end{equation}
! where we used the relation
! \begin{equation}
!   \label{Dsigma}
!   \partder{}{z} = - \frac{1}{h} \partder{}{\sigma}
! \comma
! \end{equation}
! if the motion of the free surface is ignored.
!
! The derivative of $w_\phi$ appearing in \eq{GPrimeOfOne} can be evaluted with
! the help of the formulae given in \sect{sec:wscale}. As discussed in \sect{sec:wscale},
! 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 \eq{wscale} is given by \eq{PhiStable}, and it is easy to show that
! \begin{equation}
!   \label{dWdS}
!     \partder{w}{\sigma} \Big|_{\sigma=1} =
!    - 5 h w_\phi(1) \dfrac{B_f}{u_*^4}
! \comma
! \end{equation}
! valid for both, bottom and surface boundary layers.
! Note that in the original publication of \cite{Largeetal94}, erroneously,
! there appears an additional factor $\kappa$ in this relation.
!
! With the help of \eq{dWdS}, one can re-write \eq{GPrimeOfOne} as
! \begin{equation}
!   \label{GPrimeOfOneB}
!   G'(1) = \dfrac{B_f}{u_*^4} \nu_{ti}^\phi \Big|_{z=z_{bl}}
!           - \dfrac{1}{w(1)} \partder{\nu_{ti}^\phi}{z} \Big|_{z=z_{bl}}
! \comma
! \end{equation}
! valid only for the surface boundary layer. For the bottom boundary layer,
! the minus sign in \eq{Dsigma} disappears, with the consequence that
! the minus sign in \eq{GPrimeOfOneB} has to be replaced by a plus. Note that
! if the pre-processor macro {\tt KPP\_CLIP\_GS} is defined in {\tt 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 \eq{kppNonlocal} is computed
! as described in \cite{Largeetal94}, if the pre-processor macro
! {\tt 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,
! \begin{equation}
!   \label{Rib}
!   Ri_b = \dfrac{(B_r - B(z_{bl})) d}
!                {\magn{\V U_r - \V U(z_{bl})}^2 + V_t^2(z_{bl})}
! \comma
! \end{equation}
! defined by \cite{Largeetal94}, reaches the critical value $Ri_c$.
! The subscript "r" in \eq{Rib} 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 \cite{Largeetal94}. The value of $z_{bl}$ is then found from
! \eq{Rib} by linear interpolation.
!
! To check the boundary layer limit according to the condition
! \begin{equation}
!   \label{RiConditionA}
!   Ri_b(z_{bl})
!   = \dfrac{Ri_\text{top}(z_{bl})}{Ri_\text{bot}(z_{bl})} = Ri_c
! \comma
! \end{equation}
! two methods have been implemented in GOTM. The first method simply evaluates
! \eq{RiConditionA} with a linear interpolation scheme.
! The second method is activated if the pre-processor
! macro ${\tt KPP\_IP\_FC}$ is defined. Then, the condition \eq{RiConditionA}
! is reformulated as
! \begin{equation}
!   \label{RiConditionB}
!   F_c(z_{bl}) = Ri_\text{top}(z_{bl})  - Ri_c Ri_\text{bot}(z_{bl}) = 0
! \point
! \end{equation}
! 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 {\tt clip\_mld=.true.} in {\tt kpp.nml}, the boundary layer is cut
! if it exceeds the Ekman or the Monin-Obukhov length scale, see \cite{Largeetal94}.
!
! !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 20.0  2013/12/14 00:13:46  fms
!  Merged revision 1.1.2.1 onto trunk
!
!  Revision 1.1.2.1  2012/05/15 16:00:53  smg
!  initial cvs ci for these modules to mom5.
!  AUTHOR:Griffies
!  REVIEWERS:
!  TEST STATUS:
!  CHANGES PUBLIC INTERFACES?
!  CHANGES ANSWERS?
!
!  Revision 1.1.2.2.346.1  2012/04/23 20:30:29  smg
!  updated to the gotm-2012.03.09 CVS tag.
!  AUTHOR:Martin Schmidt
!  REVIEWERS:Griffies
!  TEST STATUS:
!  CHANGES PUBLIC INTERFACES?
!  CHANGES ANSWERS?
!
!  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

!
!
!EOP
!-----------------------------------------------------------------------

   contains

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Initialise the KPP module
!
! !INTERFACE:
   subroutine init_kpp(namlst,fn,nlev,h0,h,kpp_g,kpp_rho_0)
!
! !DESCRIPTION:
! This routine first reads the namelist {\tt kpp}, which has to be contained
! in a file with filename specified by the string {\tt fn} (typically called
! {\tt kpp.nml}). Since the {\tt kpp} module uses fields defined in the
! {\tt turbulence} module, it has to allocate dynamic memory for them.
! Apart from this, this routine reports the model settings and initialises a
! number of parameters needed later in the time loop.
!
! If you call the GOTM KPP routines from a three-dimensional model, make sure
! that this function is called \emph{after} the call to {\tt init\_turbulence()}.
! Also make sure that you pass the correct arguments.
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:

!  namelist reference
   integer,          intent(in)        :: namlst

!  filename containing namelist
   character(len=*), intent(in)        :: fn

!  number of grid cells
   integer,          intent(in)        :: nlev

!  bathymetry (m)
   REALTYPE,         intent(in)        :: h0

!  size of grid cells (m)
   REALTYPE,         intent(in)        :: h(0:nlev)

!  acceleration of gravity (m/s^2)
   REALTYPE,         intent(in)        :: kpp_g

!  reference density (kg/m^3)
   REALTYPE,         intent(in)        :: kpp_rho_0

! !REVISION HISTORY:
!  Original author(s): Lars Umlauf
!
!EOP
!
! !LOCAL VARIABLES:
   integer                             :: k
   integer                             :: rc

   namelist /kpp/                      kpp_sbl,kpp_bbl,kpp_interior,    &
                                       clip_mld,Ric
!
!-----------------------------------------------------------------------
!BOC

   LEVEL1 'init_kpp...'

   ! read the variables from the namelist file
   open(namlst,file=fn,status='old',action='read',err=80)

   LEVEL2 'reading kpp namelist...'

   read(namlst,nml=kpp,err=81)
   close (namlst)

   LEVEL2 'done.'

!  allocate memory for variables defined in other modules
!
   allocate(num(0:nlev),stat=rc)
   if (rc /= 0) stop 'init_turbulence: Error allocating (num)'
   num = _ZERO_

   allocate(nuh(0:nlev),stat=rc)
   if (rc /= 0) stop 'init_turbulence: Error allocating (nuh)'
   nuh = _ZERO_

   allocate(nus(0:nlev),stat=rc)
   if (rc /= 0) stop 'init_turbulence: Error allocating (nus)'
   nus = _ZERO_

   allocate(gamu(0:nlev),stat=rc)
   if (rc /= 0) stop 'init_turbulence: Error allocating (gamu)'
   gamu = _ZERO_

   allocate(gamv(0:nlev),stat=rc)
   if (rc /= 0) stop 'init_turbulence: Error allocating (gamv)'
   gamv = _ZERO_

   allocate(gamh(0:nlev),stat=rc)
   if (rc /= 0) stop 'init_turbulence: Error allocating (gamh)'
   gamh = _ZERO_

   allocate(gams(0:nlev),stat=rc)
   if (rc /= 0) stop 'init_turbulence: Error allocating (gams)'
   gams = _ZERO_

   allocate(Rig(0:nlev),stat=rc)
   if (rc /= 0) stop 'init_turbulence: Error allocating (Rig)'
   Rig = _ZERO_

   allocate(z_w(0:nlev),stat=rc)
   if (rc /= 0) stop 'init_turbulence: Error allocating (z_w)'
   z_w = _ZERO_

   allocate(z_r(0:nlev),stat=rc)
   if (rc /= 0) stop 'init_turbulence: Error allocating (z_r)'
   z_r = _ZERO_

   allocate(h_r(0:nlev),stat=rc)
   if (rc /= 0) stop 'init_turbulence: Error allocating (h_r)'
   h_r = _ZERO_

# ifdef EXTRA_OUTPUT

   allocate(turb1(0:nlev),stat=rc)
   if (rc /= 0) stop 'init_turbulence: Error allocating (turb1)'
   turb1 = _ZERO_

   allocate(turb2(0:nlev),stat=rc)
   if (rc /= 0) stop 'init_turbulence: Error allocating (turb2)'
   turb2 = _ZERO_

   allocate(turb3(0:nlev),stat=rc)
   if (rc /= 0) stop 'init_turbulence: Error allocating (turb3)'
   turb3 = _ZERO_

   allocate(turb4(0:nlev),stat=rc)
   if (rc /= 0) stop 'init_turbulence: Error allocating (turb4)'
   turb4 = _ZERO_

   allocate(turb5(0:nlev),stat=rc)
   if (rc /= 0) stop 'init_turbulence: Error allocating (turb5)'
   turb5 = _ZERO_

# endif


!  report model parameters

   LEVEL2 '--------------------------------------------------------'
   LEVEL3 'You are using the KPP turbulence model          '
   LEVEL3 'with the following specifications:                      '
   LEVEL3 '                                                        '
   if (kpp_interior) then
      LEVEL3 'Interior mixing algorithm                  - active -   '
# ifdef KPP_SHEAR
      LEVEL4 'KPP shear instability mixing           - active -   '
# else
      LEVEL4 'KPP shear instability mixing       - not active -   '
# endif
# ifdef KPP_INTERNAL_WAVE
      LEVEL4 'KPP internal wave mixing               - active -   '
# else
      LEVEL4 'KPP internal wave mixing           - not active -   '
# endif
# ifdef KPP_CONVEC
      LEVEL4 'KPP convective instability mixing      - active -   '
# else
      LEVEL4 'KPP convective instability mixing  - not active -   '
# endif
# ifdef KPP_DDMIX
      LEVEL4 'KPP double diffusive mixing            - active -   '
# else
      LEVEL4 'KPP double diffusive mixing        - not active -   '
# endif
   else
      LEVEL3 'Interior mixing algorithm              - not active -   '
   endif

   if (kpp_sbl) then
      LEVEL3 'Surface layer mixing algorithm             - active -   '
      if (clip_mld) then
         LEVEL4 'Clipping at Ekman/Oboukhov scale       - active -   '
      else
         LEVEL4 'Clipping at Ekman/Oboukhov scale   - not active -   '
      endif
# ifdef KPP_SALINITY
      LEVEL4 'Compute salinity fluxes                - active -   '
# else
      LEVEL4 'Compute salinity fluxes            - not active -   '
# endif
# ifdef NONLOCAL
      LEVEL4 'Nonlocal fluxes                        - active -   '
# else
      LEVEL4 'Nonlocal fluxes                    - not active -   '
# endif
# ifdef KPP_TWOPOINT_REF
      LEVEL4 'Ri_b from 2-point interpolation        - active -   '
# else
      LEVEL4 'Ri_b from 2-point interpolation    - not active -   '
# endif
# ifdef KPP_IP_FC
      LEVEL4 'F_c =0 criterion for SL-depth          - active -   '
# else
      LEVEL4 'Ri_b - Ri_c =0 criterion for SL-depth  - active -   '
# endif
# ifdef KPP_CLIP_GS
      LEVEL4 'Clipping G''(sigma) for matching        - active -   '
# else
      LEVEL4 'Clipping G''(sigma) for matching    - not active -   '
# endif

      LEVEL4 'Ri_c  = ', Ric

   else
      LEVEL3 'Surface layer mixing algorithm         - not active -   '
   endif

   if (kpp_bbl) then
      LEVEL3 'Bottom layer mixing algorithm              - active -   '
      LEVEL4 '(Same parameters as surface layer mixing)'

   else
      LEVEL3 'Bottom layer mixing algorithm          - not active -   '
   endif

   LEVEL2 '--------------------------------------------------------'



!  pre-compute coefficient for turbulent shear velocity
   Vtc=Cv*sqrt(-betaT)/(sqrt(cs*epsilon)*Ric*kappa*kappa)


!  pre-compute proportionality coefficient for
!  boundary layer non-local transport
   Cg=Cstar*kappa*(cs*kappa*epsilon)**(1.0/3.0)


!  set acceleration of gravity and reference density
   g        = kpp_g
   rho_0    = kpp_rho_0
   gorho0   = g/rho_0

   LEVEL1 'done.'

   return

80 FATAL 'I could not open "kpp.nml"'
   stop 'init_kpp'
81 FATAL 'I could not read "kpp" namelist'
   stop 'init_kpp'

 end subroutine init_kpp
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Loop over the KPP-algorithm
!
! !INTERFACE:
   subroutine do_kpp(nlev,h0,h,rho,u,v,NN,NNT,NNS,SS,u_taus,u_taub,  &
                     tFlux,btFlux,sFlux,bsFlux,tRad,bRad,f)
!
! !DESCRIPTION:
! Here, the time step for the KPP model is managed. If {\tt kpp\_interior=.true.}
! in {\tt kpp.nml}, the mixing algorithm for the computation of the interior
! diffusivities is called first. This algorithm is described in \sect{sec:kppInterior}.
! Then, if {\tt kpp\_sbl=.true.} and {\tt kpp\_bbl=.true.}, the algorithms
! for the surface and bottom boundary layer are called. They are described in
! \sect{sec:kppSurface} and \sect{sec:kppBottom}, respectively.
!
! If this routine is called from a three-dimensional code, it is essential to
! pass the correct arguments. The first 3 parameters relate to the numerical
! grid, discussed in \sect{SectionNumericsMean}. Note that {\tt h0} denotes
! the local bathymetry, i.e.\ the positive distance between the reference level
! $z=0$ and the bottom.

! The next three parameters denote the potential density, $\rho$,
! and the two mean velocity components, $U$ and $V$. The buoyancy frequency, $N^2$,
! and the different contributions to it, $N_\Theta^2$ and $N_S^2$, have to be computed
! from the potential density as discussed in \sect{sec:stratification}. The shear frequency,
! $M^2$, is defined in \eq{MSquared}. The vertical discretisation does not necessarly
! have to follow \eq{shearsquared}, since in the KPP model no TKE equation is solved
! and thus energy conservation is not an issue. All three-dimensional fields have to
! be interpolated "in a smart way" to the water column defined by GOTM. The corresponding
! interpolation schemes may be quite different for the different staggered grids,
! finite volume, and finite element approaches used in the horizontal. Therefore,
! we cannot offer a general recipe here.

! The bottom friction velocity is computed as described in \sect{sec:friction}. If this
! parameter is passed from a three-dimensional code, it has to be insured that the parameter
! $r$ in \eq{uStar} is computed consistently, see \eq{rParam}.
!
! All fluxes without exception are counted positive, if they enter the water body.
! Note that for consistency, the equations of state in GOTM cannot
! be used if the KPP routines are called from a 3-D model. Therefore,
! it is necessary to pass the temperature and salinity fluxes, as well as the
! corresponding buoyancy fluxes. The same applies
! to the radiative fluxes. The user is responsible for
! performing the flux conversions in the correct way. To get an idea
! have a look at \sect{sec:convertFluxes}.
!
! The last argument is the Coriolis parameter, $f$. It is only used for clippling the mixing
! depth at the Ekman depth.
!
!
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
!  number of grid cells
   integer                                       :: nlev

!  bathymetry (m)
   REALTYPE                                      :: h0

!  thickness of grid cells (m)
   REALTYPE                                      :: h(0:nlev)

!  potential density at grid centers (kg/m^3)
   REALTYPE                                      :: rho(0:nlev)

!  velocity components at grid centers (m/s)
   REALTYPE                                      :: u(0:nlev),v(0:nlev)

!  square of buoyancy frequency (1/s^2)
   REALTYPE                                      :: NN(0:nlev)

!  square of buoyancy frequency caused by
!  temperature and salinity stratification
   REALTYPE                                      :: NNT(0:nlev),NNS(0:nlev)

!  square of shear frequency (1/s^2)
   REALTYPE                                      :: SS(0:nlev)

!  surface and bottom friction velocities (m/s)
   REALTYPE                                      :: u_taus,u_taub

!  surface temperature flux (K m/s) and
!  salinity flux (psu m/s) (negative for loss)
   REALTYPE                                      :: tFlux,sFlux

!  surface buoyancy fluxes (m^2/s^3) due to
!  heat and salinity fluxes
   REALTYPE                                      :: btFlux,bsFlux

!  radiative flux [ I(z)/(rho Cp) ] (K m/s)
!  and associated buoyancy flux (m^2/s^3)
   REALTYPE                                      :: tRad(0:nlev),bRad(0:nlev)

!  Coriolis parameter (rad/s)
   REALTYPE                                      :: f


! !REVISION HISTORY:
!  Original author(s): Lars Umlauf
!
!EOP
!
! !LOCAL VARIABLES:
   REALTYPE                            :: cff
   integer                             :: k

!
!-----------------------------------------------------------------------
!BOC
!
!-----------------------------------------------------------------------
! Update model grid
!-----------------------------------------------------------------------

!  Compute distance between centers (between rho-points)
!  Note that h is the distance between faces (between w-points)
   do k=1,nlev-1
      h_r(k) = 0.5*(h(k)+ h(k+1))
   enddo

!  Compute position of interfaces (w-points)
   z_w(0) = - h0
   do k=1,nlev
      z_w(k) = z_w(k-1) + h(k)
   enddo

!  Compute position of centers (rho-points)
   z_r(1) = - h0 + 0.5*h(1)
   do k=2,nlev
      z_r(k) = z_r(k-1) + h_r(k-1)
   enddo



!-----------------------------------------------------------------------
! compute interior mixing
!-----------------------------------------------------------------------

   if (kpp_interior) then
      call interior(nlev,NN,NNT,NNS,SS)
   else
      num = _ZERO_
      nuh = _ZERO_
      nus = _ZERO_
   endif


!-----------------------------------------------------------------------
! compute surface boundary layer mixing
!-----------------------------------------------------------------------

   if (kpp_sbl) then
      call surface_layer(nlev,h0,h,rho,u,v,NN,u_taus,u_taub,            &
                         tFlux,btFlux,sFlux,bsFlux,tRad,bRad,f)
   endif

!-----------------------------------------------------------------------
! compute bottom boundary layer mixing
!-----------------------------------------------------------------------

   if (kpp_bbl) then
      call bottom_layer(nlev,h0,h,rho,u,v,NN,u_taus,u_taub,             &
                        _ZERO_,_ZERO_,_ZERO_,_ZERO_,tRad,bRad,f)
   endif



 end subroutine do_kpp
!EOC



!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Compute interior fluxes \label{sec:kppInterior}
!
! !INTERFACE:
   subroutine interior(nlev,NN,NNT,NNS,SS)
!
! !DESCRIPTION:
! Here, the interior diffusivities (defined as the diffusivities outside the
! surface and bottom boundary layers) are computed. The algorithms are identical
! to those suggested by \cite{Largeetal94}. For numerical efficiency, the
! algorithms for different physical processes are active only if certain
! pre-processor macros are defined in {\tt cppdefs.h}.
! \begin{itemize}
!  \item The shear instability algorithm is active if the macro
!        {\tt KPP\_SHEAR} is defined.
!  \item The internal wave algorithm is active if the macro
!        {\tt KPP\_INTERNAL\_WAVE} is defined.
!  \item The convective instability algorithm is active if the macro
!        {\tt KPP\_CONVEC} is defined.
!  \item The double-diffusion algorithm is active if the macro
!        {\tt KPP\_DDMIX} is defined. Note that in this case, the
!        macro {\tt SALINITY} has to be defined as well.
! \end{itemize}

!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:

!  number of grid cells
   integer                                       :: nlev

!  square of buoyancy frequency (1/s^2)
   REALTYPE                                      :: NN(0:nlev)

!  square of buoyancy frequencies caused by
!  temperature and salinity stratification
   REALTYPE                                      :: NNT(0:nlev),NNS(0:nlev)

!  square of shear frequency (1/s^2)
   REALTYPE                                      :: SS(0:nlev)


! !REVISION HISTORY:
!  Original author(s): Lars Umlauf
!
!EOP
!
! !LOCAL VARIABLES:
   REALTYPE , parameter       :: eps=1.0E-14

   integer                    :: i
   REALTYPE                   :: cff,shear2
   REALTYPE                   :: nu_sx,nu_sxc
   REALTYPE                   :: iwm,iws
   REALTYPE                   :: drhoT,drhoS,Rrho,nu_dds,nu_ddt

!
!-----------------------------------------------------------------------
!BOC
!
!-----------------------------------------------------------------------
! Compute gradient Richardson number
!-----------------------------------------------------------------------
!
   do i=1,nlev-1
      Rig(i) = NN(i)/(SS(i) + eps)
   enddo

   Rig(0)    = _ZERO_
   Rig(nlev) = _ZERO_
!
!-----------------------------------------------------------------------
!  Compute "interior" viscosities and diffusivities everywhere as
!  the superposition of three processes: local Richardson number
!  instability due to resolved vertical shear, internal wave
!  breaking, and double diffusion.
!-----------------------------------------------------------------------
!
   do i=1,nlev-1
!
!     Smooth gradient Richardson number
      Rig(i)=0.25*Rig(i-1) + 0.50*Rig(i) + 0.25*Rig(i+1)
!
!     Compute interior diffusivity due to shear instability mixing.
# ifdef KPP_SHEAR
      cff=min(_ONE_,max(_ZERO_,Rig(i))/Ri0)
      nu_sx  = _ONE_-cff*cff
      nu_sx  = nu_sx*nu_sx*nu_sx
!
!     The shear mixing should be also a function of the actual magnitude
!     of the shear, see Polzin (1996, JPO, 1409-1425).
      shear2 = SS(i)
      cff    = shear2*shear2/(shear2*shear2+16.0E-10)
      nu_sx  = cff*nu_sx
# else KPP_SHEAR
      nu_sx=_ZERO_
# endif KPP_SHEAR

#ifdef KPP_INTERNAL_WAVE
!
!      Compute interior diffusivity due to wave breaking
!
!      Version A, see Gargett and Holloway (1984)
!      cff  =  _ONE_/sqrt(max(NN(i),1.0d-7))
!      iwm  =  1.0E-6*cff
!      iws  =  1.0E-7*cff

!     Version B, see Large et al. (1994)
      iwm  =  nuwm
      iws  =  nuws
#else
      iwm  =  _ZERO_
      iws  =  _ZERO_
#endif KPP_INTERNAL_WAVE


# ifdef KPP_CONVEC
!     Compute interior convective diffusivity due to static instability
!     mixing
      cff    =  max(NN(i),bvfcon)
      cff    =  min(_ONE_,(bvfcon-cff)/bvfcon)
      nu_sxc =  _ONE_-cff*cff
      nu_sxc =  nu_sxc*nu_sxc*nu_sxc
# else KPP_CONVEC
      nu_sxc =  _ZERO_
# endif KPP_CONVEC
!
!     Sum contributions due to internal wave breaking, shear instability
!     and convective diffusivity due to shear instability.
      num(i)=iwm+nu0m*nu_sx+nu0c*nu_sxc
      nuh(i)=iws+nu0s*nu_sx+nu0c*nu_sxc
      nus(i)=nuh(i)
!
# ifdef KPP_DDMIX
!
!-----------------------------------------------------------------------
!  Compute double-diffusive mixing.  It can occur when vertical
!  gradient of density is stable but the vertical gradient of
!  salinity (salt figering) or temperature (diffusive convection)
!  is unstable.
!-----------------------------------------------------------------------
!
!     Compute double-diffusive density ratio, Rrho.
      drhoT = NNT(i)
      drhoS = NNS(i)
      Rrho= - drhoT/drhoS
!
!
!     Salt fingering case.
      if ((Rrho.gt._ONE_).and.(drhoS.gt._ZERO_)) then
!
!        Compute interior diffusivity for double diffusive mixing of
!        salinity.  Upper bound "Rrho" by "Rrho0"; (Rrho0=1.9, nuf=0.001).
         Rrho=min(Rrho,Rrho0)
         nu_dds=_ONE_-((Rrho-_ONE_)/(Rrho0-_ONE_))**2.0
         nu_dds=nuf*nu_dds*nu_dds*nu_dds
!
!        Compute interior diffusivity for double diffusive mixing
!        of temperature (fdd=0.7).
         nu_ddt=fdd*nu_dds
!
!
!     Diffusive convection case.
      elseif ((_ZERO_.lt.Rrho).and.(Rrho.lt._ONE_).and.(drhoS.lt._ZERO_)) then
!
!        Compute interior diffusivity for double diffusive mixing of
!        temperature (Marmorino and Caldwell, 1976); (nu=1.5e-6,
!        tdd1=0.909, tdd2=4.6, tdd3=0.54).
         nu_ddt=nu*tdd1*exp(tdd2*exp(-tdd3*((_ONE_/Rrho)-_ONE_)))

!        Compute interior diffusivity for double diffusive mixing
!        of salinity (sdd1=0.15, sdd2=1.85, sdd3=0.85).
!
         if (Rrho.lt.0.5) then
            nu_dds=nu_ddt*sdd1*Rrho
         else
            nu_dds=nu_ddt*(sdd2*Rrho-sdd3)
         endif
      else
         nu_ddt=_ZERO_
         nu_dds=_ZERO_
      endif
!
!     Add double diffusion contribution to temperature and salinity
!     mixing coefficients.
      nuh(i)=nuh(i)  + nu_ddt
      nus(i)=nuh(i)  + nu_dds

# endif KPP_DDMIX

   enddo ! loop over interior points

 end subroutine interior
!EOC


!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Compute turbulence in the surface layer \label{sec:kppSurface}
!
! !INTERFACE:
   subroutine surface_layer(nlev,h0,h,rho,u,v,NN,u_taus,u_taub,       &
                            tFlux,btFlux,sFlux,bsFlux,tRad,bRad,f)
!
! !DESCRIPTION:
! In this routine all computations related to turbulence in the surface layer
! are performed. The algorithms are described in \sect{sec:kpp}. Note that these
! algorithms are affected by some pre-processor macros defined in {\tt cppdefs.inp},
! and by the parameters set in {\tt kpp.nml}, see \sect{sec:kpp}.
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
!  number of grid cells
   integer                                       :: nlev

!  bathymetry (m)
   REALTYPE                                      :: h0

!  thickness of grid cells (m)
   REALTYPE                                      :: h(0:nlev)

!  potential density at grid centers (kg/m^3)
   REALTYPE                                      :: rho(0:nlev)

!  velocity components at grid centers (m/s)
   REALTYPE                                      :: u(0:nlev),v(0:nlev)

!  square of buoyancy frequency (1/s^2)
   REALTYPE                                      :: NN(0:nlev)

!  surface and bottom friction velocities (m/s)
   REALTYPE                                      :: u_taus,u_taub

!  surface temperature flux (K m/s) and
!  salinity flux (sal m/s) (negative for loss)
   REALTYPE                                      :: tFlux,sFlux

!  surface buoyancy fluxes (m^2/s^3) due to
!  heat and salinity fluxes
   REALTYPE                                      :: btFlux,bsFlux

!  radiative flux [ I(z)/(rho Cp) ] (K m/s)
!  and associated buoyancy flux (m^2/s^3)
   REALTYPE                                      :: tRad(0:nlev),bRad(0:nlev)

!  Coriolis parameter (rad/s)
   REALTYPE                                      :: f


!
! !REVISION HISTORY:
!  Original author(s): Lars Umlauf
!
!EOP
!
! !LOCAL VARIABLES:
   REALTYPE, parameter          :: eps      = 1.0E-10

   integer                      :: k,ksbl

   REALTYPE                     :: Bo
   REALTYPE                     :: Bfsfc
   REALTYPE                     :: tRadSrf,tRadSbl
   REALTYPE                     :: bRadSrf,bRadSbl
   REALTYPE                     :: Gm1
   REALTYPE                     :: Gt1
   REALTYPE                     :: Gs1
   REALTYPE                     :: dGm1dS
   REALTYPE                     :: dGt1dS
   REALTYPE                     :: dGs1dS
   REALTYPE                     :: f1
   REALTYPE                     :: sl_dpth,sl_z
   REALTYPE                     :: swdk
   REALTYPE                     :: wm
   REALTYPE                     :: ws
   REALTYPE                     :: zgrid,depth

   REALTYPE                     :: Gm, Gt, Gs, K_bl, Ribot, Ritop, Rk, Rref
   REALTYPE                     :: Uk, Uref, Ustarb, Vk, Vref
   REALTYPE                     :: dR,dU,dV
   REALTYPE                     :: a1, a2, a3
   REALTYPE                     :: cff,cff_up,cff_dn
   REALTYPE                     :: c1,c2,c3
   REALTYPE                     :: dK_bl, hekman, hmonob, sigma, zbl

   REALTYPE, dimension (0:nlev) :: Bflux
   REALTYPE, dimension (0:nlev) :: FC


!-----------------------------------------------------------------------
!BOC
!
!
!-----------------------------------------------------------------------
!  Get approximation of surface layer depth using "epsilon" and
!  boundary layer depth from previous time step.
!-----------------------------------------------------------------------
!
   sl_dpth = epsilon*(z_w(nlev)-zsbl)
   sl_z    = epsilon*zsbl
!
!-----------------------------------------------------------------------
!  Compute total buoyancy flux at W-points.
!  Bo is negative, if heat is lost or salinity is gained (convection).
!  It does not include the short wave radiative flux at the surface.
!-----------------------------------------------------------------------
!
    tRadSrf   =   tRad(nlev)
    bRadSrf   =   bRad(nlev)

!  surface buoyancy flux (negative for buoyancy loss)
   Bo         = btFlux + bsFlux

!  include effect of short wave radiation
!  prepare non-local fluxes
   do k = 0,nlev
      Bflux(k)  = Bo  + ( bRadSrf - bRad(k) )
# ifdef NONLOCAL
      cff       = _ONE_-(0.5+sign(0.5d0,Bflux(k)))
      gamh(k)   =  -cff*( tFlux + tRadSrf - tRad(k) )
#  ifdef KPP_SALINITY
      gams(k)   =  -cff*sFlux
#  endif
# endif

   enddo

!-----------------------------------------------------------------------
!  Compute potential density and velocity components surface reference
!  values.
!-----------------------------------------------------------------------

!  simply use uppermost value
   Rref = rho(nlev)
   Uref =   u(nlev)
   Vref =   v(nlev)


!-----------------------------------------------------------------------
!  Compute critical function, FC, for bulk Richardson number
!-----------------------------------------------------------------------

   FC(0   ) = _ZERO_
   FC(nlev) = _ZERO_

   do k=nlev-1,2,-1

      depth=z_w(nlev)-z_w(k)

      if (Bflux(k).lt._ZERO_) then
         sigma=min(sl_dpth,depth)
      else
         sigma=depth
      endif

      call wscale (Bflux(k),u_taus,sigma,wm,ws)


#ifdef KPP_TWOPOINT_REF

!     interpolate reference value at grid face "k"
!     from values at grid centers

      cff = _ONE_/h_r(k)

      dR  = cff*( rho(k+1) - rho(k) )
      dU  = cff*( u  (k+1) - u  (k) )
      dV  = cff*( v  (k+1) - v  (k) )

      cff = _ONE_/2.0

      Rk  = rho(k) + h(k)*cff*dR
      Uk  =   u(k) + h(k)*cff*dU
      Vk  =   v(k) + h(k)*cff*dV

#else
!     identify reference value at grid face "k"
!     with value at center above

      Rk = rho(k+1)
      Uk =   u(k+1)
      Vk =   v(k+1)

!!$      c1 =  0.6
!!$      c2 =  0.2
!!$      c3 =  0.2
!!$
!!$      Rk = c1*rho(k+1) + c2*rho(k) + c3*rho(k-1)
!!$      Uk = c1*  u(k+1) + c2*  u(k) + c3*  u(k-1)
!!$      Vk = c1*  v(k+1) + c2*  v(k) + c3*  v(k-1)
#endif

!     compute numerator and denominator of Ri_b
      Ritop = -gorho0*(Rref-Rk)*depth

      Ribot = (Uref-Uk)**2+(Vref-Vk)**2 +                               &
              Vtc*depth*ws*sqrt(abs(NN(k)))

# ifdef KPP_IP_FC
      FC(k) = Ritop-Ric*Ribot
# else
      FC(k) = Ritop/(Ribot+eps)
# endif

   enddo  ! inner grid faces



!-----------------------------------------------------------------------
! Linearly interpolate to find "zsbl" where Rib/Ric=1.
!-----------------------------------------------------------------------

   ksbl = 1         ! ksbl is index of cell containing zsbl
   zsbl = z_w(0)

# ifdef KPP_IP_FC
!  look for position of vanishing F_crit
   do k=nlev,2,-1
      if ((ksbl.eq.1).and.(FC(k-1).gt._ZERO_)) then
         zsbl = (z_w(k)*FC(k-1)-z_w(k-1)*FC(k))/(FC(k-1)-FC(k))
         ksbl = k
      endif
   enddo
# else
!  look for position of vanishing Ri_b
   do k=nlev,2,-1
      if ((ksbl.eq.1).and.((FC(k).lt.Ric).and.(FC(k-1).ge.Ric))) then
         zsbl = ((FC(k-1)-Ric)*z_w(k) -                                 &
                         (FC(k)-Ric)*z_w(k-1))/(FC(k-1)-FC(k))
         ksbl = k
      endif
   enddo
# endif



!-----------------------------------------------------------------------
!  Compute total buoyancy flux at surface boundary layer depth
!-----------------------------------------------------------------------

!  interpolate from interface values to zsbl
   bRadSbl = ( bRad(ksbl-1)*(z_w(ksbl)-zsbl) +                                &
               bRad(ksbl  )*(zsbl-z_w(ksbl-1) ) )/ h(ksbl)

   Bfsfc   = Bo + (bRadSrf - bRadsbl)


!-----------------------------------------------------------------------
!  Limit boundary layer depth by Ekman and Monin-Obukhov depths
!  (under neutral and stable conditions)
!-----------------------------------------------------------------------

   if (clip_mld) then
      if ((u_taus.gt._ZERO_).and.(Bfsfc.gt._ZERO_)) then
         hekman = cekman*u_taus/max(abs(f),eps)
         hmonob = cmonob*u_taus*u_taus*u_taus/max(kappa*Bfsfc,eps)
         zsbl   = (z_w(nlev)-min(hekman,hmonob,z_w(nlev)-zsbl))
      endif
   endif


   zsbl = min(zsbl,z_w(nlev))
   zsbl = max(zsbl,z_w(0   ))

!  find new boundary layer index "ksbl".
   ksbl=1
   do k=nlev,2,-1
      if ((ksbl.eq.1).and.(z_w(k-1).lt.zsbl)) then
         ksbl = k
      endif
   enddo


!-----------------------------------------------------------------------
!  Compute total buoyancy flux at surface boundary layer depth
!-----------------------------------------------------------------------


   bRadSbl = ( bRad(ksbl-1)*(z_w(ksbl)-zsbl) +                          &
               bRad(ksbl  )*(zsbl-z_w(ksbl-1) ) )/ h(ksbl)


   Bfsfc   = Bo + (bRadSrf - bRadSbl)

!-----------------------------------------------------------------------
!  Compute tubulent velocity scales (wm,ws) at "zsbl".
!-----------------------------------------------------------------------

   zbl     = z_w(nlev)-zsbl
   sl_dpth = epsilon*zbl

   if (Bfsfc.gt._ZERO_) then
      cff  = _ONE_
   else
      cff  = epsilon
   endif

   sigma=cff*(z_w(nlev)-zsbl)

   call wscale (Bfsfc,u_taus,sigma,wm,ws)


!-----------------------------------------------------------------------
!  Compute nondimensional shape function Gx(sigma) in terms of the
!  interior diffusivities at sigma=1 (Gm1, Gs1, Gt1) and its vertical
!  derivative evaluated "zsbl" via interpolation.
!-----------------------------------------------------------------------

!   original code with kappa-bug
!   f1 = 5.0*max(_ZERO_,Bfsfc)*kappa/(u_taus*u_taus*u_taus*u_taus+eps)

!  new code code without kappa-bug
   f1 = 5.0*max(_ZERO_,Bfsfc)/(u_taus*u_taus*u_taus*u_taus+eps)

   if (zsbl.gt.z_w(1)) then
      ! if boundary layer does not touch lowest grid box

!     abbreviation for "ksbl"
      k      = ksbl
      cff    = _ONE_/h(k)
      cff_dn = cff*(zsbl-z_w(k-1))
      cff_up = cff*(z_w(k)-zsbl  )

!  Compute nondimensional shape function for viscosity "Gm1" and its
!  vertical derivative "dGm1dS" evaluated at "zsbl".

      K_bl   = cff_dn*num(k)+cff_up*num(k-1)
      dK_bl  = cff*(num(k)-num(k-1))
      Gm1    = K_bl/(zbl*wm+eps)

# ifdef KPP_CLIP_GS
      dGm1dS = min(_ZERO_,-dK_bl/(wm+eps)-K_bl*f1)
# else
      dGm1dS = -dK_bl/(wm+eps) + K_bl*f1
# endif

!  Compute nondimensional shape function for diffusion of temperature
!  "Gt1" and its vertical derivative "dGt1dS" evaluated at "zsbl".
!
      K_bl   = cff_dn*nuh(k)+cff_up*nuh(k-1)
      dK_bl  = cff*(nuh(k)-nuh(k-1))
      Gt1    = K_bl/(zbl*ws+eps)

# ifdef KPP_CLIP_GS
      dGt1dS = min(_ZERO_,-dK_bl/(ws+eps)-K_bl*f1)
# else
      dGt1dS = -dK_bl/(ws+eps) + K_bl*f1
# endif

# ifdef KPP_SALINITY
!
!  Compute nondimensional shape function for diffusion of salinity
!  "Gs1" and its vertical derivative "dGs1dS" evaluated at "zsbl".
!
      K_bl   = cff_dn*nus(k)+cff_up*nus(k-1)
      dK_bl  = cff*(nus(k)-nus(k-1))
      Gs1    = K_bl/(zbl*ws+eps)

# ifdef KPP_CLIP_GS
      dGs1dS = min(_ZERO_,-dK_bl/(ws+eps)-K_bl*f1)
# else
      dGs1dS = -dK_bl/(ws+eps) + K_bl*f1
# endif

# endif


   else
!     If the surface boundary layer extends to the bottom, assume that
!     the neutral boundary layer similarity theory holds at the bottom.
!
      dK_bl  = kappa*u_taub
      K_bl   = dK_bl*(zsbl-z_w(0))

!     Compute nondimensional bottom shape function for diffusion of
!     momentum

      Gm1    = K_bl/(zbl*wm+eps)

# ifdef KPP_CLIP_GS
      dGm1dS = min(_ZERO_,-dK_bl/(wm+eps)-K_bl*f1)
# else
      dGm1dS = -dK_bl/(wm+eps) + K_bl*f1
# endif


!     Compute nondimensional bottom shape function for diffusion of
!     temperature.

      Gt1    = K_bl/(zbl*ws+eps)

# ifdef KPP_CLIP_GS
      dGs1dS = min(_ZERO_,-dK_bl/(ws+eps)-K_bl*f1)
# else
      dGs1dS = -dK_bl/(ws+eps) + K_bl*f1
# endif

!     Compute nondimensional bottom shape function for diffusion of
!     salinity.
# ifdef KPP_SALINITY
      Gs1    = Gt1
      dGs1dS = dGt1dS
# endif

   endif
!
!-----------------------------------------------------------------------
!  Compute surface boundary layer mixing coefficients
!-----------------------------------------------------------------------
!
!  loop over the inner interfaces
!  (the outermost are not needed)
   do k=1,nlev-1

      if (k.ge.ksbl) then    ! interface above cell containing zsbl

!        Compute turbulent velocity scales at vertical W-points.
         depth=z_w(nlev)-z_w(k)
         if (Bflux(k).lt._ZERO_) then
            sigma=min(sl_dpth,depth)
         else
            sigma=depth
         endif

         call wscale(Bflux(k),u_taus,sigma,wm, ws)
!
!        Set polynomial coefficients for shape function.
         sigma = depth/(zbl+eps)

         a1    = sigma-2.0
         a2    = 3.0-2.0*sigma
         a3    = sigma-1.0
!
!        Compute nondimesional shape functions.
         Gm = a1+a2*Gm1+a3*dGm1dS
         Gt = a1+a2*Gt1+a3*dGt1dS
# ifdef KPP_SALINITY
         Gs = a1+a2*Gs1+a3*dGs1dS
# endif

!        Compute boundary layer mixing coefficients, combine them
!        with interior mixing coefficients.
         num(k) = depth*wm*(_ONE_+sigma*Gm)
         nuh(k) = depth*ws*(_ONE_+sigma*Gt)
# ifdef KPP_SALINITY
         nus(k) = depth*ws*(_ONE_+sigma*Gs)
# endif
# ifdef NONLOCAL
!        Compute boundary layer nonlocal transport (m/s2).
         cff      = Cg*(_ONE_-(0.5+sign(0.5d0,Bflux(k))))/(zbl*ws+eps)
         gamh(k) = cff*nuh(k)*gamh(k)
#  ifdef KPP_SALINITY
         gams(k) = cff*nus(k)*gams(k)
#  endif
# endif


!     Set non-local transport terms to zero outside boundary  layer.
      else
# ifdef NONLOCAL
         gamh(k) = _ZERO_
#  ifdef KPP_SALINITY
         gams(k) = _ZERO_
#  endif
# endif
      endif


   enddo

!     not non-local fluxes at top and bottom
      gamh(0   ) = _ZERO_
      gams(0   ) = _ZERO_
      gamh(nlev) = _ZERO_
      gams(nlev) = _ZERO_


 end subroutine surface_layer
!EOC




!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Compute turbulence in the bottom layer \label{sec:kppBottom}
!
! !INTERFACE:
   subroutine bottom_layer(nlev,h0,h,rho,u,v,NN,u_taus,u_taub, &
                            tFlux,btFlux,sFlux,bsFlux,tRad,bRad,f)
!
! !DESCRIPTION:
! In this routine all computations related to turbulence in the bottom layer
! are performed. The algorithms are described in \sect{sec:kpp}. Note that these
! algorithms are affected by some pre-processor macros defined in {\tt cppdefs.inp},
! and by the parameters set in {\tt kpp.nml}, see \sect{sec:kpp}.

! The computation of the bulk Richardson number is slightly different from the
! surface boundary layer, since for the bottom boundary layer this quantity
! is defined as,
! \begin{equation}
!   \label{RibBottom}
!   Ri_b = \dfrac{(B(z_{bl})-B_r) d}
!                {\magn{\V U(z_{bl})-\V U_r}^2 + V_t^2(z_{bl})}
! \comma
! \end{equation}
! where $z_{bl}$ denotes the position of the edge of the bottom boundary layer.
!
! Also different from the surface layer computations is the absence of non-local
! fluxes.

! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:

!  number of grid cells
   integer                                       :: nlev

!  bathymetry (m)
   REALTYPE                                      :: h0

!  thickness of grid cells (m)
   REALTYPE                                      :: h(0:nlev)

!  potential density at grid centers (kg/m^3)
   REALTYPE                                      :: rho(0:nlev)

!  velocity components at grid centers (m/s)
   REALTYPE                                      :: u(0:nlev),v(0:nlev)

!  square of buoyancy frequency (1/s^2)
   REALTYPE                                      :: NN(0:nlev)

!  surface and bottom friction velocities (m/s)
   REALTYPE                                      :: u_taus,u_taub

!  bottom temperature flux (K m/s) and
!  salinity flux (sal m/s) (negative for loss)
   REALTYPE                                      :: tFlux,sFlux

!  bottom buoyancy fluxes (m^2/s^3) due to
!  heat and salinity fluxes
   REALTYPE                                      :: btFlux,bsFlux

!  radiative flux [ I(z)/(rho Cp) ] (K m/s)
!  and associated buoyancy flux (m^2/s^3)
   REALTYPE                                      :: tRad(0:nlev),bRad(0:nlev)

!  Coriolis parameter (rad/s)
   REALTYPE                                      :: f

!
! !REVISION HISTORY:
!  Original author(s): Lars Umlauf
!
!EOP
!
! !LOCAL VARIABLES:
   REALTYPE, parameter          :: eps      = 1.0E-10

   integer                      :: k,kbbl

   REALTYPE                     :: Bo
   REALTYPE                     :: Bfbot
   REALTYPE                     :: tRadBot,tRadBbl
   REALTYPE                     :: bRadBot,bRadBbl
   REALTYPE                     :: Gm1
   REALTYPE                     :: Gt1
   REALTYPE                     :: Gs1
   REALTYPE                     :: dGm1dS
   REALTYPE                     :: dGt1dS
   REALTYPE                     :: dGs1dS
   REALTYPE                     :: f1
   REALTYPE                     :: bl_dpth,bl_z
   REALTYPE                     :: swdk
   REALTYPE                     :: wm
   REALTYPE                     :: ws
   REALTYPE                     :: zgrid,depth

   REALTYPE                     :: Gm, Gt, Gs, K_bl, Ribot, Ritop, Rk, Rref
   REALTYPE                     :: Uk, Uref, Ustarb, Vk, Vref
   REALTYPE                     :: dR,dU,dV
   REALTYPE                     :: a1, a2, a3
   REALTYPE                     :: cff,cff_up, cff_dn
   REALTYPE                     :: dK_bl, hekman, hmonob, sigma, zbl

   REALTYPE, dimension (0:nlev) :: Bflux
   REALTYPE, dimension (0:nlev) :: FC


!-----------------------------------------------------------------------
!BOC
!
!
!-----------------------------------------------------------------------
!  Get approximation of bottom layer depth using "epsilon" and
!  boundary layer depth from previous time step.
!-----------------------------------------------------------------------
!
   bl_dpth = epsilon*(zbbl-z_w(0))
   bl_z    = epsilon*zbbl

!-----------------------------------------------------------------------
!  Compute total buoyancy flux (m2/s3) at W-points.
!-----------------------------------------------------------------------
!
    tRadBot   =   tRad(0)
    bRadBot   =   bRad(0)

!  bottom buoyancy flux
!  (negative for buoyancy gain)
   Bo         = - (btFlux + bsFlux)

!  include effect of short-wave radiation
   do k = 0,nlev
      Bflux(k)  = Bo + ( bRad(k) - bRadBot )
   enddo


!-----------------------------------------------------------------------
!  Compute potential density and velocity components bottom reference
!  values.
!-----------------------------------------------------------------------


!  simply use lowest value
   Rref = rho(1)
   Uref =   u(1)
   Vref =   v(1)


!-----------------------------------------------------------------------
!  Compute critical function, FC, for bulk Richardson number
!-----------------------------------------------------------------------

   FC(0   ) = _ZERO_
   FC(nlev) = _ZERO_

   do k=1,nlev-1

      depth=z_w(k)-z_w(0)

      if (Bflux(k).lt._ZERO_) then
         sigma=min(bl_dpth,depth)
      else
         sigma=depth
      endif

      call wscale (Bflux(k),u_taub,sigma,wm,ws)



#ifdef KPP_TWOPOINT_REF

!     interpolate reference value at grid face "k"
!     from values at grid centers

      cff = _ONE_/h_r(k)

      dR  = cff*( rho(k+1) - rho(k) )
      dU  = cff*( u  (k+1) - u  (k) )
      dV  = cff*( v  (k+1) - v  (k) )

      cff = _ONE_/2.0

      Rk  = rho(k) + h(k)*cff*dR
      Uk  =   u(k) + h(k)*cff*dU
      Vk  =   v(k) + h(k)*cff*dV

#else
!     identify reference value at grid face "k"
!     with value at center below
      Rk  = rho(k)
      Uk  =   u(k)
      Vk  =   v(k)
#endif

!     compute numerator and denominator of Ri_b
      Ritop = -gorho0*(Rk-Rref)*depth
      Ribot = (Uk-Uref)**2+(Vk-Vref)**2 +                              &
              Vtc*depth*ws*sqrt(abs(NN(k)) )

# ifdef KPP_IP_FC
      FC(k)=Ritop-Ric*Ribot
# else
      FC(k)=Ritop/(Ribot+eps)
# endif

   enddo ! inner grid faces



!-----------------------------------------------------------------------
! Linearly interpolate to find "zbbl" where Rib/Ric=1.
!-----------------------------------------------------------------------

   kbbl = nlev           ! kbbl is  index of cell containing zbbl
   zbbl = z_w(nlev)

# ifdef KPP_IP_FC
!  look for position of vanishing F_crit
   do k=1,nlev-1
      if ((kbbl.eq.nlev).and.(FC(k).gt._ZERO_)) then
         zbbl = (z_w(k)*FC(k-1)-z_w(k-1)*FC(k))/(FC(k-1)-FC(k))
         kbbl = k
      endif
   enddo
# else
!  look for position of vanishing Ri_b
   do k=1,nlev-1
      if ((kbbl.eq.nlev).and.((FC(k-1).lt.Ric).and.(FC(k).ge.Ric))) then
         zbbl = ( (Ric  -  FC(k-1) )*z_w(k  ) +                         &
                  (FC(k) - Ric     )*z_w(k-1) )/(FC(k)-FC(k-1))
         kbbl = k
      endif
   enddo
# endif


!-----------------------------------------------------------------------
!  Limit boundary layer thickness by Ekman scale
!-----------------------------------------------------------------------

   if (clip_mld) then
      if (u_taub.gt._ZERO_) then
         hekman = cekman*u_taub/max(abs(f),eps)
         zbbl   = min(z_w(0)+hekman,zbbl)
      endif
   endif

   zbbl = min(zbbl,z_w(nlev))
   zbbl = max(zbbl,z_w(0   ))

!  find new boundary layer index "kbbl".
   kbbl=nlev
   do k=1,nlev-1
      if ((kbbl.eq.nlev).and.(z_w(k).gt.zbbl)) then
         kbbl = k
      endif
   enddo

!-----------------------------------------------------------------------
!  Compute total buoyancy flux at bottom boundary layer depth
!-----------------------------------------------------------------------

   bRadBbl = ( bRad(kbbl-1)*(z_w(kbbl)-zbbl)  +                         &
               bRad(kbbl  )*(zbbl-z_w(kbbl-1) ) )/ h(kbbl)

   Bfbot   = Bo + (bRadBbl-bRadBot)


!-----------------------------------------------------------------------
!  Compute tubulent velocity scales (wm,ws) at "zbbl".
!-----------------------------------------------------------------------

   zbl     = zbbl-z_w(0)
   bl_dpth = epsilon*zbl

   if (Bfbot.gt._ZERO_) then
      cff  = _ONE_
   else
      cff  = epsilon
   endif

   sigma=cff*zbl

   call wscale (Bfbot,u_taub,sigma,wm,ws)


!-----------------------------------------------------------------------
!  Compute nondimensional shape function Gx(sigma) in terms of the
!  interior diffusivities at sigma=1 (Gm1, Gs1, Gt1) and its vertical
!  derivative evaluated "zbbl" via interpolation.
!-----------------------------------------------------------------------

!  original code with kappa-bug
!   f1 = 5.0*max(_ZERO_,Bfbot)*kappa/(u_taub*u_taub*u_taub*u_taub+eps)

!   new code without kappa-bug
   f1 = 5.0*max(_ZERO_,Bfbot)/(u_taub*u_taub*u_taub*u_taub+eps)

   k      = kbbl
   cff    = _ONE_/h(k)
   cff_dn = cff*(zbbl-z_w(k-1))
   cff_up = cff*(z_w(k)-zbbl  )

!  Compute nondimensional shape function for viscosity "Gm1" and its
!  vertical derivative "dGm1dS" evaluated at "zbbl".

   K_bl   =  cff_dn*num(k)+cff_up*num(k-1)
   dK_bl  =  cff*(num(k)-num(k-1))
   Gm1    =  K_bl/(zbl*wm+eps)

# ifdef KPP_CLIP_GS
   dGm1dS = min(_ZERO_,dK_bl/(wm+eps)+K_bl*f1)
# else
   dGm1dS = dK_bl/(wm+eps) + K_bl*f1
# endif

!
!  Compute nondimensional shape function for diffusion of temperature
!  "Gt1" and its vertical derivative "dGt1dS" evaluated at "zbbl".
!
   K_bl   =  cff_dn*nuh(k)+cff_up*nuh(k-1)
   dK_bl  =  cff*(nuh(k)-nuh(k-1))
   Gt1    =   K_bl/(zbl*ws+eps)


# ifdef KPP_CLIP_GS
   dGt1dS = min(_ZERO_,dK_bl/(ws+eps)+K_bl*f1)
# else
   dGt1dS = dK_bl/(ws+eps) + K_bl*f1
# endif

# ifdef KPP_SALINITY
!
!  Compute nondimensional shape function for diffusion of salinity
!  "Gs1" and its vertical derivative "dGs1dS" evaluated at "zbbl".
!
   K_bl   =  cff_dn*nus(k)+cff_up*nus(k-1)
   dK_bl  =  cff*(nus(k)-nus(k-1))
   Gs1    =   K_bl/(zbl*ws+eps)

# ifdef KPP_CLIP_GS
   dGs1dS = min(_ZERO_,dK_bl/(ws+eps)+K_bl*f1)
# else
   dGs1dS = dK_bl/(ws+eps) + K_bl*f1
# endif

# endif

!
!-----------------------------------------------------------------------
!  Compute bottom boundary layer mixing coefficients
!-----------------------------------------------------------------------
!
!  loop over the inner interfaces
!  (the outermost are not needed)
   do k=1,nlev-1

      if (k.lt.kbbl) then

!        Compute turbulent velocity scales at vertical W-points.
         depth=z_w(k)-z_w(0)
         if (Bflux(k).lt._ZERO_) then
            sigma=min(bl_dpth,depth)
         else
            sigma=depth
         endif

         call wscale(Bflux(k),u_taub,sigma,wm, ws)
!
!        Set polynomial coefficients for shape function.
         sigma = depth/(zbl+eps)

         a1    = sigma-2.0
         a2    = 3.0-2.0*sigma
         a3    = sigma-1.0
!
!        Compute nondimesional shape functions.

         Gm = a1+a2*Gm1+a3*dGm1dS
         Gt = a1+a2*Gt1+a3*dGt1dS
# ifdef KPP_SALINITY
         Gs = a1+a2*Gs1+a3*dGs1dS
# endif
!
!        Compute boundary layer mixing coefficients, combine them
!        with interior mixing coefficients.
         num(k) = depth*wm*(_ONE_+sigma*Gm)
         nuh(k) = depth*ws*(_ONE_+sigma*Gt)
# ifdef KPP_SALINITY
         nus(k) = depth*ws*(_ONE_+sigma*Gs)
# endif

      endif
   enddo



 end subroutine bottom_layer
!EOC


!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Compute the velocity scale\label{sec:wscale}
!
! !INTERFACE:
   subroutine wscale(Bfsfc,u_taus,d,wm,ws)
!
! !DESCRIPTION:
!  This routine computes the turbulent velocity scale for momentum
!  and tracer as a function of the turbulent friction velocity, $u_*$,
!  the "limited" distance, $d_\text{lim}$, and the total buoyancy
!  flux, $B_f$, according to
!  \begin{equation}
!    \label{wscale}
!     w_\phi = \dfrac{\kappa u_*}{\Phi_\phi (\zeta)}
!     \point
!  \end{equation}
! In this equation, $\Phi_\phi$ is a non-dimensional function of the
! stability parameter $\zeta=d_\text{lim}/L$, using the Monin-Obukhov
! length,
!  \begin{equation}
!    \label{MOLength}
!     L = \dfrac{u_*^3}{\kappa B_f}
!     \point
!  \end{equation}
! In stable situations, $B_f \ge 0$, the length scale $d_\text{lim}$ is just the distance
! from the surface or bottom, $d$. Then, the non-dimensional function is of the form
! \begin{equation}
!   \label{PhiStable}
!   \Phi_\phi = 1 + \zeta
! \comma
! \end{equation}
! and identical for momentum, heat, and tracers.
!
! In unstable situations, $B_f < 0$, the scale $d_\text{lim}$ corresponds
! to the distance from surface or bottom only until it reaches the end of
! the surface (or bottom) layer at $d=\epsilon h$. Then it stays constant
! at this maximum value.
!
! The different functional forms of $\Phi_\phi(\zeta)$ for unstable flows
! are discussed in \cite{Largeetal94}.


! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:

!  buoyancy flux (m^2/s^3)
   REALTYPE, intent(in)                :: Bfsfc

!  friction velocity (m/s)
   REALTYPE, intent(in)                :: u_taus

!  (limited) distance (m)
   REALTYPE, intent(in)                :: d
!
! !OUTPUT PARAMETERS:

!  velocity scale (m/s)
!  for momentum and tracer
   REALTYPE, intent(out)               :: wm, ws
!
! !REVISION HISTORY:
!  Original author(s): Lars Umlauf
!
!EOP
!
! !LOCAL VARIABLES:
   REALTYPE, parameter                 :: eps = 1.0E-20
   REALTYPE, parameter                 :: r3  = 1.0/3.0

   REALTYPE                            :: u_taus3, zetahat, zetapar
!
!-----------------------------------------------------------------------
!BOC
!
!  pre-compute some quantities for faster execution
   u_taus3=u_taus*u_taus*u_taus
   zetahat=kappa*d*Bfsfc
   zetapar=zetahat/(u_taus3+eps)

   if (zetahat.ge._ZERO_) then
!     stable or neutral water column
      wm=kappa*u_taus/(_ONE_+5.0*zetapar)
      ws=wm
   else
!     instable water column
      if (zetapar.gt.zetam) then
         wm=kappa*u_taus*(_ONE_-16.0*zetapar)**0.25
      else
         wm=kappa*(am*u_taus3-cm*zetahat)**r3
      endif
      if (zetapar.gt.zetas) then
         ws=kappa*u_taus*(_ONE_-16.0*zetapar)**0.5
      else
         ws=kappa*(as*u_taus3-cs*zetahat)**r3
      endif
   endif

 end subroutine wscale
!EOC

! !IROUTINE: Clean up the kpp module
!
! !INTERFACE:
   subroutine clean_kpp()
!
! !DESCRIPTION:
!  De-allocate all memory allocated in init\_kpp().
!
! !USES:
   IMPLICIT NONE
!
! !REVISION HISTORY:
!  Original author(s): Jorn Bruggeman
!
!EOP
!-----------------------------------------------------------------------
!BOC
   LEVEL1 'clean_kpp'

   LEVEL2 'de-allocating KPP memory ...'
   if (allocated(num)) deallocate(num)
   if (allocated(nuh)) deallocate(nuh)
   if (allocated(nus)) deallocate(nus)
   if (allocated(gamu)) deallocate(gamu)
   if (allocated(gamv)) deallocate(gamv)
   if (allocated(gamh)) deallocate(gamh)
   if (allocated(gams)) deallocate(gams)
   if (allocated(Rig)) deallocate(Rig)
   if (allocated(z_w)) deallocate(z_w)
   if (allocated(z_r)) deallocate(z_r)
   if (allocated(h_r)) deallocate(h_r)

# ifdef EXTRA_OUTPUT

   if (allocated(turb1)) deallocate(turb1)
   if (allocated(turb2)) deallocate(turb2)
   if (allocated(turb3)) deallocate(turb3)
   if (allocated(turb4)) deallocate(turb4)
   if (allocated(turb5)) deallocate(turb5)

# endif

   return
   end subroutine clean_kpp

!-----------------------------------------------------------------------


!-----------------------------------------------------------------------

 end module kpp

!-----------------------------------------------------------------------
! Copyright by the GOTM-team under the GNU Public License - www.gnu.org
!-----------------------------------------------------------------------

