Calculation of the stratification


   subroutine stratification(nlev,buoy_method,dt,cnpar,nub,gamB)

This routine computes the mean potential density, $ {\langle \rho \rangle} $, the mean potential buoyancy, $ B$, defined in (33), and the mean buoyancy frequency,

$\displaystyle N^2 = - \dfrac{g}{\rho_0} \dfrac{\partial {\rho}}{\partial {z}} = \dfrac{\partial {B}}{\partial {z}} \; , \quad$ (40)

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, (244), according to

$\displaystyle {\langle \rho \rangle} = \hat{\rho} (\Theta,S,P_R) \; , \quad$ (41)

where $ \Theta$ denotes the mean potential temperature, $ S$ the mean salinity and $ P_R$ the mean reference pressure. The buoyancy frequency defined in (40) can be decomposed into contributions due to potential temperature and salinity stratification,

$\displaystyle N^2 = N_\Theta^2 + N_S^2 \; , \quad$ (42)

where we introduced the quantities

$\displaystyle N_\Theta^2 = - \dfrac{g}{\rho_0} \dfrac{\partial {\rho}}{\partial...
...{S} = g \alpha(\Theta,S,P_R) \dfrac{\partial {\Theta}}{\partial {z}} \; , \quad$ (43)

with the thermal expansion coefficient defined in (246), and

$\displaystyle N_S^2 = - \dfrac{g}{\rho_0} \dfrac{\partial {\rho}}{\partial {z}}...
...\Theta} = - g \beta(\Theta,S,P_R) \dfrac{\partial {S}}{\partial {z}} \; , \quad$ (44)

with the saline contraction coefficient defined in (247). 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 section 8.8, the model is truely based on potential temperature and density. Otherwise, the model is based on in-situ quantities.

Alternatively to the procedure outlined above, depending on the values of the parameter buoy_method, the buoyancy may be calculated by means of the transport equation (34). This equation then replaces the computation of $ \Theta$ and $ S$ and is only recommended for idealized studies.


   use meanflow,   only: h,S,T,buoy,rho
   use meanflow,   only: NN,NNT,NNS
   use meanflow,   only: gravity,rho_0
   use eqstate,    only: eqstate1
   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)
   Original author(s): Karsten Bolding & Hans Burchard

Karsten Bolding 2012-12-28