*******************
Primitive Equations
*******************

At resolutions larger than 1 km (more marginally above 100 m), The ocean is a fluid that can 
be described as a good approximation by the primitive equations (PE). The PE equations are 
simplifications from the Navier-Stokes equations made from scale considerations, along 
with a nonlinear equation of state, which couples the two active tracers (temperature and salinity):

* Hydrostatic hypothesis: the vertical momentum equation is reduced to a balance between 
  the vertical pressure gradient and the buoyancy force (non-hydrostatic processes such as 
  convection must be parametrized)

* Boussinesq hypothesis: density variations are neglected except in their contribution 
  to the buoyancy force

* Incompressibility hypothesis (stemming from the former): the three-dimensional divergence 
  of the velocity vector is assumed to be zero.

* Spherical earth approximation: the geopotential surfaces are assumed to be spheres so 
  that gravity (local vertical) is parallel to the earth's radius
 
* Thin-shell approximation: the ocean depth is neglected compared to the earth's radius

* Turbulent closure hypothesis: the turbulent fluxes (which represent the effect of 
  small-scale processes on the large scale) are expressed in terms of large-scale features

By default (#undef NBQ), CROCO solves the primitive equations. However, it also 
has the ability to relax the first 3 hypotheses (#define NBQ). When SOLVE3D is 
not activated, CROCO can be used as a classical shallow water model.
 


Equations in Cartesian coordinates
----------------------------------

* The momentum balance in zonal x and meridional y directions, written in terms of 
  grid-scale (resolved) and subgrid-scale velocity components: 

.. math::

  \frac{\partial u}{\partial t} 
              + \vec{\bf \nabla} . \left ( \vec{\textbf v} u \right ) - f v   = 
              - \frac{\partial \phi}{\partial x}  + \mathcal{F}_u +  \mathcal{D}_u

  \frac{\partial v}{\partial t} 
              + \vec{\bf \nabla} . \left ( \vec{\textbf v} v \right ) + f u = 
              - \frac{\partial \phi}{\partial y}  + \mathcal{F}_v +  \mathcal{D}_v

  Turbulent closure schemes are applied to parametrized subgrid-scale vertical fluxes.

* The time evolution of a scalar concentration field, :math:`C(x,y,z,t)` (e.g., salinity, 
  temperature, or nutrients), is governed by the advective-diffusive equation :

.. math::

  \frac{\partial C}{\partial t} 
               + \vec{\bf \nabla} . \left ( \vec{\textbf v} C \right ) = 
               \mathcal{F}_C +  \mathcal{D}_C


* The equation of state is given by :

.. math::

  \rho = \rho(T,S,P)


* In the Boussinesq approximation, density variations are neglected in the momentum 
  equations except in their contribution to the buoyancy force in the vertical momentum 
  equation. Under the hydrostatic approximation, it is further assumed that the vertical 
  pressure gradient balances the buoyancy force :

.. math::

  \frac{\partial \phi}{\partial z} = - \frac{\rho g}{\rho_0}

* The final equation expresses the continuity equation. For an incompressible 
  fluid (Boussinesq approximation):

.. math:: 

   \vec{\bf \nabla} . \vec{\textbf v} = 
            \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}  
          + \frac{\partial w}{\partial z} = 0


The variables used are :

:math:`\mathcal{D}_u, \mathcal{D}_v, \mathcal{D}_C` : diffusive terms

:math:`\mathcal{F}_u, \mathcal{F}_v, \mathcal{F}_C` : forcing terms

:math:`f(x,y)` : Traditional Coriolis parameter :math:`2 \Omega sin \phi`

:math:`g` : acceleration of gravity

:math:`\phi(x,y,z,t)` : dynamic pressure :math:`\phi=P/\rho_0`, with P the total pressure

:math:`\rho_0+\rho(x,y,z,t)` : total in situ density

:math:`u,v,w` : the (x,y,z) components of vector velocity :math:`\vec{\textbf v}`



Equations in terrain-following coordinates
------------------------------------------

We first introduce a generalized stretched vertical coordinate system (:math:`s`), 
which sets the variable bottom flat at :math:`z=-h(x,y)`. :math:`s` spans the range 
from -1 (bottom) to 0 (surface), and the transformation rules are: 

.. math::

  \begin{aligned}
  {\left(\frac{\partial}{\partial x}\right)_{z}=\left(\frac{\partial}{\partial x}\right)_{s}-\left(\frac{1}{H_{z}}\right)\left(\frac{\partial z}{\partial x}\right)_{s} \frac{\partial}{\partial s}} \\ 
  {\left(\frac{\partial}{\partial y}\right)_{z}=\left(\frac{\partial}{\partial y}\right)_{s}-\left(\frac{1}{H_{z}}\right)\left(\frac{\partial z}{\partial y}\right)_{s} \frac{\partial}{\partial s}} \\ 
  {\frac{\partial}{\partial z}=\left(\frac{\partial s}{\partial z}\right) \frac{\partial}{\partial s}=\frac{1}{H_{z}} \frac{\partial}{\partial s}} \\ 
 \text{where } {H_{z} \equiv \frac{\partial z}{\partial s}}
  \end{aligned}

The vertical velocity in :math:`s` coordinate is:

.. math::

 \begin{aligned} \Omega(x, y, s, t) &=\frac{1}{H_{z}}\left[w-(1+s) \frac{\partial \zeta}{\partial t}-u \frac{\partial z}{\partial x}-v \frac{\partial z}{\partial y}\right] \\ 
  w &=\frac{\partial z}{\partial t}+u \frac{\partial z}{\partial x}+v \frac{\partial z}{\partial y}+\Omega H_{z} \end{aligned}

:math:`\Omega=0` at both surface and bottom. 

Next, the requirement for a laterally variable grid resolution can also be met for 
suitably smooth domains by introducing an appropriate orthogonal coordinate 
transformation in the horizontal. Let the new coordinates be :math:`\xi(x,y)` 
and :math:`\eta(x,y)` where the relationship of horizontal arc length to the differential 
distance is given by:

.. math::

 \qquad(d s)_{\xi}=\left(\frac{1}{m}\right) d \xi

 \qquad(d s)_{\eta}=\left(\frac{1}{n}\right) d \eta


Here :math:`m(\xi,\eta)` and :math:`n(\xi,\eta)` are the scale factors which relate the 
differential distances  :math:`(\Delta \xi, \Delta \eta)` to the actual (physical) arc lengths.

.. math::

  \begin{aligned} \frac{\partial}{\partial t}\left(\frac{H_{z} u}{m n}\right)+\frac{\partial}{\partial \xi}\left(\frac{H_{z} u^{2}}{n}\right)+\frac{\partial}{\partial \eta}\left(\frac{H_{z} u v}{m}\right)+\frac{\partial}{\partial s}\left(\frac{H_{z} u \Omega}{m n}\right) \\-\left\{\left(\frac{f}{m n}\right)+v \frac{\partial}{\partial \xi}\left(\frac{1}{n}\right)-u \frac{\partial}{\partial \eta}\left(\frac{1}{m}\right)\right\} H_{z} v=\\-\left(\frac{H_{z}}{n}\right)\left(\frac{\partial \phi}{\partial \xi}+\frac{g \rho}{\rho_{o}} \frac{\partial z}{\partial \xi}+g \frac{\partial \zeta}{\partial \xi}\right)+\frac{H_{z}}{m n}\left(\mathcal{F}_{u}+\mathcal{D}_{u}\right) \\ \frac{\partial}{\partial t}\left(\frac{H_{z} v}{m n}\right)+\frac{\partial}{\partial \xi}\left(\frac{H_{z} u v}{n}\right)+\frac{\partial}{\partial \eta}\left(\frac{H_{z} v^{2}}{m}\right)+\frac{\partial}{\partial s}\left(\frac{H_{z} v \Omega}{m n}\right) \\+\left\{\left(\frac{f}{m n}\right)+v \frac{\partial}{\partial \xi}\left(\frac{1}{n}\right)-u \frac{\partial}{\partial \eta}\left(\frac{1}{m}\right)\right\} H_{z} u=\\-\left(\frac{H_{z}}{m}\right)\left(\frac{\partial \phi}{\partial \eta}+\frac{g \rho}{\rho_{o}} \frac{\partial z}{\partial \eta}+g \frac{\partial \zeta}{\partial \eta}\right)+\frac{H_{z}}{m n}\left(\mathcal{F}_{v}+\mathcal{D}_{v}\right) \\
 \frac{\partial}{\partial t}\left(\frac{H_{z} T}{m n}\right)+\frac{\partial}{\partial \xi}\left(\frac{H_{z} u T}{n}\right)+\frac{\partial}{\partial \eta}\left(\frac{H_{z} v T}{m}\right)+\frac{\partial}{\partial s}\left(\frac{H_{z} \Omega T}{m n}\right) &=\frac{H_{z}}{m n}\left(\mathcal{F}_{T}+\mathcal{D}_{T}\right) \\ \frac{\partial}{\partial t}\left(\frac{H_{z} S}{m n}\right)+\frac{\partial}{\partial \xi}\left(\frac{H_{z} u S}{n}\right)+\frac{\partial}{\partial \eta}\left(\frac{H_{z} v S}{m}\right)+\frac{\partial}{\partial s}\left(\frac{H_{z} \Omega S}{m n}\right) &=\frac{H_{z}}{m n}\left(\mathcal{F}_{S}+\mathcal{D}_{S}\right) \\ \rho &=\rho(T, S, P) \\ \frac{\partial \phi}{\partial s} &=-\left(\frac{g H_{z} \rho}{\rho_{o}}\right) \\ \frac{\partial}{\partial t}\left(\frac{H_{z}}{m n}\right)+\frac{\partial}{\partial \xi}\left(\frac{H_{z} u}{n}\right)+\frac{\partial}{\partial \eta}\left(\frac{H_{z} v}{m}\right)+\frac{\partial}{\partial s}\left(\frac{H_{z} \Omega}{m n}\right) &=0 \end{aligned}