stratification.F90 6.8 KB
Newer Older
dumoda01's avatar
dumoda01 committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208
!$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
!  added new copyright to files
!
!  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
!-----------------------------------------------------------------------