stratification.F90 6.8 KB
!$Id: stratification.F90,v 1.7 2005-07-18 08:54:33 lars Exp$
#include"cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: Calculation of the stratification\label{sec:stratification}
!
! !INTERFACE:
subroutine stratification(nlev,buoy_method,dt,cnpar,nub,gamB)
!
! !DESCRIPTION:
! This routine computes the mean potential density, $\mean{\rho}$, the mean
! potential buoyancy, $B$, defined in \eq{DefBuoyancy}, and the mean buoyancy
! frequency,
!  \begin{equation}
!    \label{DefBuoyancyFrequency}
!     N^2 = - \dfrac{g}{\rho_0} \partder{\rho}{z} = \partder{B}{z}
!     \comma
!  \end{equation}
! which is based on potential density or buoyancy such that for $N^2=0$, the entropy
! is constant in the whole water column and mixing does not work against buoyancy
! forces. If GOTM used as a turbulence library in your own three-dimensional model,
! you have to insure that the $N^2$ computed by you, and passed to the turbulence
! routines in GOTM, is consistent with the concept of potential density and your
! equation of state.
!
! The mean potential density is evaluated from the equation of state, \eq{DefEOS},
! according to
!  \begin{equation}
!    \label{DefPotentialDensity}
!     \mean{\rho} = \hat{\rho} (\Theta,S,P_R)
!     \comma
!  \end{equation}
!  where $\Theta$ denotes the mean potential temperature, $S$ the mean salinity
!  and $P_R$ the mean reference pressure. The buoyancy frequency defined in
! \eq{DefBuoyancyFrequency} can be decomposed into contributions due to
!  potential temperature and salinity stratification,
!  \begin{equation}
!    \label{NDecompostionA}
!     N^2 = N_\Theta^2 + N_S^2
!     \comma
!  \end{equation}
!  where we introduced the quantities
!  \begin{equation}
!    \label{NNT}
!     N_\Theta^2  = - \dfrac{g}{\rho_0} \partder{\rho}{z} \Big|_{S}
!                 = g \alpha(\Theta,S,P_R) \partder{\Theta}{z}
!     \comma
!  \end{equation}
!  with the thermal expansion coefficient defined in \eq{eosAlpha}, and
!  \begin{equation}
!    \label{NNS}
!     N_S^2  = - \dfrac{g}{\rho_0} \partder{\rho}{z} \Big|_{\Theta}
!                 = - g \beta(\Theta,S,P_R) \partder{S}{z}
!  \comma
!  \end{equation}
!  with the saline contraction coefficient defined in \eq{eosBeta}. It is important
!  to note that in the actual code the reference pressure, $P_R$, has been replaced by
!  the (approximate) hydrostatic pressure. Only if this dependence is replaced by
!  the constant reference pressure at the surface in the equation of state,
!  see \sect{sec:eqstate}, the model is truely based on potential temperature and density.
!  Otherwise,  the model is based on \emph{in-situ} quantities.
!
!  Alternatively to the procedure outlined above, depending on the values of the
!  parameter {\tt buoy\_method}, the buoyancy may be calculated by means of the
!  transport equation \eq{bEq}. This equation then replaces the computation of $\Theta$
!  and $S$ and is only recommended for idealized studies.
!
! !USES:
use meanflow,   only: h,S,T,buoy,rho
use meanflow,   only: NN,NNT,NNS
use meanflow,   only: gravity,rho_0
use eqstate,    only: eqstate1
IMPLICIT NONE
!
! !INPUT PARAMETERS:

!  number of vertical layers
integer,  intent(in)                :: nlev

!  method to compute buoyancy
integer,  intent(in)                :: buoy_method

!   time step (s)
REALTYPE, intent(in)                :: dt

!  numerical "implicitness" parameter
REALTYPE, intent(in)                :: cnpar

!  diffusivity of buoyancy (m^2/s)
REALTYPE, intent(in)                :: nub(0:nlev)

!  non-local buoyancy flux (m^2/s^3)
REALTYPE, intent(in)                :: gamb(0:nlev)

!
! !OUTPUT PARAMETERS:

!
! !REVISION HISTORY:
!  Original author(s): Karsten Bolding & Hans Burchard
!
!  $Log: stratification.F90,v$
!  Revision 1.7  2005-07-18 08:54:33  lars
!  changed docu for html compliance
!
!  Revision 1.6  2005/06/27 13:44:07  kbk
!  modified + removed traling blanks
!
!  Revision 1.5  2003/03/28 09:20:35  kbk
!
!  Revision 1.4  2003/03/28 08:56:56  kbk
!  removed tabs
!
!  Revision 1.3  2003/03/10 08:50:07  gotm
!  Improved documentation and cleaned up code
!
!  Revision 1.2  2001/11/18 11:50:37  gotm
!  Cleaned
!
!  Revision 1.1.1.1  2001/02/12 15:55:57  gotm
!  initial import into CVS
!
!EOP
!
! !LOCAL VARIABLES:
integer                   :: i
REALTYPE                  :: buoyp,buoym,Sface,Tface
REALTYPE                  :: zCenter,zFace,dz
REALTYPE                  :: pFace,pCenter
integer,  parameter       :: USEEQSTATE=1
REALTYPE, parameter       :: oneTenth=_ONE_/10.
!
!-----------------------------------------------------------------------
!BOC
if (buoy_method .EQ. USEEQSTATE) then

!     initialize parameters for uppermost grid box
zFace   = _ZERO_
zCenter = 0.5*h(nlev)
pCenter = oneTenth*zCenter

buoy(nlev) = eqstate1(S(nlev),T(nlev),pCenter,gravity,rho_0)
rho(nlev)  = rho_0 - rho_0/gravity*buoy(nlev)

do i=nlev-1,1,-1

! compute distances between centers
dz     = 0.5*(h(i)+h(i+1))

! compute local depths
zFace   = zFace   + h(i+1)
zCenter = zCenter + dz

! compute approximate local pressure
pFace   = oneTenth*zFace
pCenter = oneTenth*zCenter

! interpolate T and S from centers to faces
Sface  = ( S(i+1)*h(i) + S(i)*h(i+1) ) / ( h(i+1) + h(i) )
Tface  = ( T(i+1)*h(i) + T(i)*h(i+1) ) / ( h(i+1) + h(i) )

! T contribution to buoyancy frequency
buoyp  = eqstate1(Sface,T(i+1),pFace,gravity,rho_0)
buoym  = eqstate1(Sface,T(i  ),pFace,gravity,rho_0)
NNT(i) = (buoyp-buoym)/dz

! S contribution to buoyancy frequency
buoyp  = eqstate1(S(i+1),Tface,pFace,gravity,rho_0)
buoym  = eqstate1(S(i  ),Tface,pFace,gravity,rho_0)
NNS(i) = (buoyp-buoym)/dz

! total buoyancy frequency is the sum
NN(i) = NNT(i) + NNS(i)

! compute buoyancy and density
buoy(i) = eqstate1(S(i),T(i),pCenter,gravity,rho_0)
rho(i)  = rho_0 - rho_0/gravity*buoy(i)
end do

else

! for some idealized cases, compute buoyancy from
! prognostic equation
call buoyancy(nlev,dt,cnpar,nub,gamb)

! compute density and buoyancy frequency
rho(nlev)  = rho_0 - rho_0/gravity*buoy(nlev)

do i=nlev-1,1,-1
dz     = 0.5*(h(i)+h(i+1))
NN(i)  = (buoy(i+1)-buoy(i))/dz
rho(i) = rho_0 - rho_0/gravity*buoy(i)
end do
end if

! update boundary values
NN(0)    = _ZERO_
NN(nlev) = _ZERO_

return
end subroutine stratification
!EOC

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