s_algorithm/text/nonlin_frsurf.tex

% $Header: /u/gcmpack/manual/part2/nonlin_frsurf.tex,v 1.12 2006/04/04 20:19:00 jmc Exp $
% $Name:  $


\subsection{Non-linear free-surface}
\label{sect:nonlinear-freesurface}

Recently, new options have been added to the model
that concern the free surface formulation.


\subsubsection{pressure/geo-potential and free surface}
\label{sect:phi-freesurface}

For the atmosphere, since $\phi = \phi_{topo} - \int^p_{p_s} \alpha dp$,
subtracting the reference state defined in section
\ref{sec:hpe-p-geo-potential-split}~:\\
$$
\phi_o = \phi_{topo} - \int^p_{p_o} \alpha_o dp  
\hspace{5mm}\mathrm{with}\hspace{3mm} \phi_o(p_o)=\phi_{topo}
$$
we get:
$$
\phi' = \phi - \phi_o = \int^{p_s}_p \alpha dp - \int^{p_o}_p \alpha_o dp
$$
For the ocean, the reference state is simpler since $\rho_c$ does not dependent
on $z$ ($b_o=g$) and the surface reference position is uniformly $z=0$ ($R_o=0$),
and the same subtraction leads to a similar relation.
For both fluid, using the isomorphic notations, we can write:
$$ 
\phi' = \int^{r_{surf}}_r b~ dr - \int^{R_o}_r b_o dr
$$
and re-write as:
\begin{equation}
\phi' = \int^{r_{surf}}_{R_o} b~ dr + \int^{R_o}_r (b - b_o) dr
\label{eq:split-phi-Ro}
\end{equation}
or:
\begin{equation}
\phi' = \int^{r_{surf}}_{R_o} b_o dr + \int^{r_{surf}}_r (b - b_o) dr
\label{eq:split-phi-bo}
\end{equation}

In section \ref{sec:finding_the_pressure_field}, following eq.\ref{eq:split-phi-Ro},
the pressure/geo-potential $\phi'$ has been separated into surface ($\phi_s$),
and hydrostatic anomaly ($\phi'_{hyd}$).
In this section, the split between $\phi_s$ and $\phi'_{hyd}$ is 
made according to equation \ref{eq:split-phi-bo}. This slightly
different definition reflects the actual implementation in the code 
and is valid for both linear and non-linear
free-surface formulation, in both r-coordinate and r*-coordinate.

Because the linear free-surface approximation ignore the tracer content
of the fluid parcel between $R_o$ and $r_{surf}=R_o+\eta$, 
for consistency reasons, this part is also neglected in $\phi'_{hyd}$~:
$$
\phi'_{hyd} = \int^{r_{surf}}_r (b - b_o) dr \simeq \int^{R_o}_r (b - b_o) dr
$$
Note that in this case, the two definitions of $\phi_s$ and $\phi'_{hyd}$ 
from equation \ref{eq:split-phi-Ro} and \ref{eq:split-phi-bo} converge toward 
the same (approximated) expressions: $\phi_s = \int^{r_{surf}}_{R_o} b_o dr$ 
and $\phi'_{hyd}=\int^{R_o}_r b' dr$.\\
On the contrary, the unapproximated formulation ("non-linear free-surface",
see the next section) retains the full expression:
$\phi'_{hyd} = \int^{r_{surf}}_r (b - b_o) dr $~.
This is obtained by selecting {\bf nonlinFreeSurf}=4 in parameter
file {\em data}.\\

Regarding the surface potential: 
$$\phi_s = \int_{R_o}^{R_o+\eta} b_o dr = b_s \eta
\hspace{5mm}\mathrm{with}\hspace{5mm} 
b_s = \frac{1}{\eta} \int_{R_o}^{R_o+\eta} b_o dr $$
$b_s \simeq b_o(R_o)$ is an excellent approximation (better than 
the usual numerical truncation, since generally $|\eta|$ is smaller 
than the vertical grid increment).

For the ocean, $\phi_s = g \eta$ and $b_s = g$ is uniform.
For the atmosphere, however, because of topographic effects, the
reference surface pressure $R_o=p_o$ has large spatial variations that
are responsible for significant $b_s$ variations (from 0.8 to 1.2
$[m^3/kg]$). For this reason, when {\bf uniformLin\_PhiSurf} {\em=.FALSE.} 
(parameter file {\em data}, namelist {\em PARAM01})
a non-uniform linear coefficient $b_s$ is used and computed
({\it S/R INI\_LINEAR\_PHISURF}) according to the reference surface 
pressure $p_o$:
$b_s = b_o(R_o) = c_p \kappa (p_o / P^o_{SL})^{(\kappa - 1)} \theta_{ref}(p_o)$.
with $P^o_{SL}$ the mean sea-level pressure.


\subsubsection{Free surface effect on column total thickness
(Non-linear free-surface)}

The total thickness of the fluid column is $r_{surf} - R_{fixed} =
\eta + R_o - R_{fixed}$. In most applications, the free surface 
displacements are small compared to the total thickness
$\eta \ll H_o = R_o - R_{fixed}$. 
In the previous sections and in older version of the model,
the linearized free-surface approximation was made, assuming
$r_{surf} - R_{fixed} \simeq H_o$ when computing horizontal transports,
either in the continuity equation or in tracer and momentum 
advection terms.
This approximation is dropped when using the non-linear free-surface
formulation and the total thickness, including the time varying part
$\eta$, is considered when computing horizontal transports.
Implications for the barotropic part are presented hereafter.
In section \ref{sect:freesurf-tracer-advection} consequences for 
tracer conservation is briefly discussed (more details can be 
found in \cite{campin:02})~; the general time-stepping is presented 
in section \ref{sect:nonlin-freesurf-timestepping} with some 
limitations regarding the vertical resolution in section 
\ref{sect:nonlin-freesurf-dz_surf}.

In the non-linear formulation, the continuous form of the model 
equations remains unchanged, except for the 2D continuity equation 
(\ref{eq:discrete-time-backward-free-surface}) which is now
integrated from $R_{fixed}(x,y)$ up to $r_{surf}=R_o+\eta$ :

\begin{displaymath}
\epsilon_{fs} \partial_t \eta =
\left. \dot{r} \right|_{r=r_{surf}} + \epsilon_{fw} (P-E) =
- {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o+\eta} \vec{\bf v} dr
+ \epsilon_{fw} (P-E)
\end{displaymath}

Since $\eta$ has a direct effect on the horizontal velocity (through
$\nabla_h \Phi_{surf}$), this adds a non-linear term to the free
surface equation. Several options for the time discretization of this
non-linear part can be considered, as detailed below.

If the column thickness is evaluated at time step $n$, and with
implicit treatment of the surface potential gradient, equations
(\ref{eq-solve2D} and \ref{eq-solve2D_rhs}) becomes:
\begin{eqnarray*}
\epsilon_{fs} {\eta}^{n+1} -
{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n}+R_o-R_{fixed})
{\bf \nabla}_h b_s {\eta}^{n+1}
= {\eta}^*
\end{eqnarray*}
where
\begin{eqnarray*}
{\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
\Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o+\eta^n} \vec{\bf v}^* dr
\: + \: \epsilon_{fw} \Delta_t (P-E)^{n}
\end{eqnarray*} 
This method requires us to update the solver matrix at each time step.

Alternatively, the non-linear contribution can be evaluated fully
explicitly:
\begin{eqnarray*}
\epsilon_{fs} {\eta}^{n+1} -
{\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{fixed})
{\bf \nabla}_h b_s {\eta}^{n+1}
= {\eta}^*
+{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n})
{\bf \nabla}_h b_s {\eta}^{n}
\end{eqnarray*} 
This formulation allows one to keep the initial solver matrix
unchanged though throughout the integration, since the non-linear free
surface only affects the RHS.

Finally, another option is a "linearized" formulation where the total
column thickness appears only in the integral term of the RHS
(\ref{eq-solve2D_rhs}) but not directly in the equation
(\ref{eq-solve2D}).

Those different options (see Table \ref{tab:nonLinFreeSurf_flags})
have been tested and show little differences. However, we recommend
the use of the most precise method (the 1rst one) since the 
computation cost involved in the solver matrix update is negligible.

\begin{table}[htb]
%\begin{center}
\centering
 \begin{tabular}[htb]{|l|c|l|}
   \hline
   parameter & value & description \\ 
   \hline
                   & -1 & linear free-surface, restart from a pickup file \\
                   &    & produced with \#undef EXACT\_CONSERV code\\
   \cline{2-3}
                   & 0 & Linear free-surface \\ 
   \cline{2-3}
    nonlinFreeSurf & 4 & Non-linear free-surface \\
   \cline{2-3}
                   & 3 & same as 4 but neglecting
                           $\int_{R_o}^{R_o+\eta} b' dr $ in $\Phi'_{hyd}$ \\
   \cline{2-3}
                   & 2 & same as 3 but do not update cg2d solver matrix \\
   \cline{2-3}
                  & 1 & same as 2 but treat momentum as in Linear FS \\
   \hline
                  & 0 & do not use $r*$ vertical coordinate (= default)\\
   \cline{2-3}
    select\_rStar & 2 & use $r^*$ vertical coordinate \\
   \cline{2-3}
                  & 1 & same as 2 but without the contribution of the\\
                  &   & slope of the coordinate in $\nabla \Phi$ \\
   \hline
  \end{tabular}
  \caption{Non-linear free-surface flags}
  \label{tab:nonLinFreeSurf_flags}
%\end{center}
\end{table}


\subsubsection{Tracer conservation with non-linear free-surface}
\label{sect:freesurf-tracer-advection}

To ensure global tracer conservation (i.e., the total amount) as well
as local conservation, the change in the surface level thickness must
be consistent with the way the continuity equation is integrated, both
in the barotropic part (to find $\eta$) and baroclinic part (to find
$w = \dot{r}$).

To illustrate this, consider the shallow water model, with a source
of fresh water (P):
$$
\partial_t h + \nabla \cdot h \vec{\bf v} = P
$$
where $h$ is the total thickness of the water column.
To conserve the tracer $\theta$ we have to discretize:
$$
\partial_t (h \theta) + \nabla \cdot ( h \theta \vec{\bf v})
  = P \theta_{\mathrm{rain}}
$$
Using the implicit (non-linear) free surface described above (section
\ref{sect:pressure-method-linear-backward}) we have:
\begin{eqnarray*}
h^{n+1} = h^{n} - \Delta t \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) + \Delta t P \\
\end{eqnarray*}
The discretized form of the tracer equation must adopt the same
``form'' in the computation of tracer fluxes, that is, the same value
of $h$, as used in the continuity equation:
\begin{eqnarray*}
h^{n+1} \, \theta^{n+1} = h^n \, \theta^n 
        - \Delta t \nabla \cdot (h^n \, \theta^n \, \vec{\bf v}^{n+1})
        + \Delta t P \theta_{rain}
\end{eqnarray*}

The use of a 3 time-levels time-stepping scheme such as the Adams-Bashforth
make the conservation sightly tricky.
The current implementation with the Adams-Bashforth time-stepping
provides an exact local conservation and prevents any drift in
the global tracer content (\cite{campin:02}).
Compared to the linear free-surface method, an additional step is required:
the variation of the water column thickness (from $h^n$ to $h^{n+1}$) is
not incorporated directly into the tracer equation.  Instead, the
model uses the $G_\theta$ terms (first step) as in the linear free
surface formulation (with the "{\it surface correction}" turned "on",
see tracer section):
$$
G_\theta^n = \left(- \nabla \cdot (h^n \, \theta^n \, \vec{\bf v}^{n+1}) 
         - \dot{r}_{surf}^{n+1} \theta^n \right) / h^n
$$
Then, in a second step, the thickness variation (expansion/reduction)
is taken into account:
$$
\theta^{n+1} = \theta^n + \Delta t \frac{h^n}{h^{n+1}} 
   \left( G_\theta^{(n+1/2)} + P (\theta_{\mathrm{rain}} - \theta^n )/h^n \right)
%\theta^{n+1} = \theta^n + \frac{\Delta t}{h^{n+1}} 
%   \left( h^n G_\theta^{(n+1/2)} + P (\theta_{\mathrm{rain}} - \theta^n ) \right)
$$
Note that with a simple forward time step (no Adams-Bashforth), 
these two formulations are equivalent,  
since
$
(h^{n+1} - h^{n})/ \Delta t = 
P - \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) = P + \dot{r}_{surf}^{n+1}
$

\subsubsection{Time stepping implementation of the
non-linear free-surface}     
\label{sect:nonlin-freesurf-timestepping}

The grid cell thickness was hold constant with the linear
free-surface~; with the non-linear free-surface, it is now varying 
in time, at least at the surface level.
This implies some modifications of the general algorithm described 
earlier in sections \ref{sect:adams-bashforth-sync} and 
\ref{sect:adams-bashforth-staggered}.

A simplified version of the staggered in time, non-linear 
free-surface algorithm is detailed hereafter, and can be compared 
to the equivalent linear free-surface case (eq. \ref{eq:Gv-n-staggered}
to \ref{eq:t-n+1-staggered}) and can also be easily transposed
to the synchronous time-stepping case.
Among the simplifications, salinity equation, implicit operator
and detailed elliptic equation are omitted. Surface forcing is
explicitly written as fluxes of temperature, fresh water and
momentum, $Q^{n+1/2}, P^{n+1/2}, F_{\bf v}^n$ respectively.
$h^n$ and $dh^n$ are the column and grid box thickness in r-coordinate.
%-------------------------------------------------------------
\begin{eqnarray}
\phi^{n}_{hyd} & = & \int b(\theta^{n},S^{n},r) dr
\label{eq:phi-hyd-nlfs} \\
\vec{\bf G}_{\vec{\bf v}}^{n-1/2}\hspace{-2mm} & = & 
\vec{\bf G}_{\vec{\bf v}} (dh^{n-1},\vec{\bf v}^{n-1/2})
\hspace{+2mm};\hspace{+2mm}
\vec{\bf G}_{\vec{\bf v}}^{(n)} =  
   \frac{3}{2} \vec{\bf G}_{\vec{\bf v}}^{n-1/2} 
-  \frac{1}{2} \vec{\bf G}_{\vec{\bf v}}^{n-3/2}
\label{eq:Gv-n-nlfs} \\
%\vec{\bf G}_{\vec{\bf v}}^{(n)} & = & 
%   \frac{3}{2} \vec{\bf G}_{\vec{\bf v}}^{n-1/2} 
%-  \frac{1}{2} \vec{\bf G}_{\vec{\bf v}}^{n-3/2}
%\label{eq:Gv-n+5-nlfs} \\
%\vec{\bf v}^{*} & = & \vec{\bf v}^{n-1/2} + \frac{\Delta t}{dh^{n}} \left(
%dh^{n-1}\vec{\bf G}_{\vec{\bf v}}^{(n)} + F_{\vec{\bf v}}^{n} \right)
\vec{\bf v}^{*} & = & \vec{\bf v}^{n-1/2} + \Delta t \frac{dh^{n-1}}{dh^{n}} \left(
\vec{\bf G}_{\vec{\bf v}}^{(n)} + F_{\vec{\bf v}}^{n}/dh^{n-1} \right)
- \Delta t \nabla \phi_{hyd}^{n}
\label{eq:vstar-nlfs}
\end{eqnarray}
\hspace{3cm}$\longrightarrow$~~{\it update model~geometry~:~}${\bf hFac}(dh^n)$\\
\begin{eqnarray}
\eta^{n+1/2} \hspace{-2mm} & = & 
\eta^{n-1/2} + \Delta t P^{n+1/2} - \Delta t
  \nabla \cdot \int \vec{\bf v}^{n+1/2} dh^{n} \nonumber \\
             & = & \eta^{n-1/2} + \Delta t P^{n+1/2} - \Delta t
  \nabla \cdot \int \!\!\! \left( \vec{\bf v}^* - g \Delta t \nabla \eta^{n+1/2} \right) dh^{n}
\label{eq:nstar-nlfs} \\
\vec{\bf v}^{n+1/2}\hspace{-2mm} & = & 
\vec{\bf v}^{*} - g \Delta t \nabla \eta^{n+1/2}
\label{eq:v-n+1-nlfs} \\
h^{n+1} & = & h^{n} + \Delta t P^{n+1/2} - \Delta t
  \nabla \cdot \int \vec{\bf v}^{n+1/2} dh^{n}
\label{eq:h-n+1-nlfs} \\
G_{\theta}^{n} & = & G_{\theta} ( dh^{n}, u^{n+1/2}, \theta^{n} )
\hspace{+2mm};\hspace{+2mm}
G_{\theta}^{(n+1/2)} = \frac{3}{2} G_{\theta}^{n} - \frac{1}{2} G_{\theta}^{n-1}
\label{eq:Gt-n-nlfs} \\
%\theta^{n+1} & = &\theta^{n} + \frac{\Delta t}{dh^{n+1}} \left( dh^n 
%G_{\theta}^{(n+1/2)} + Q^{n+1/2} + P^{n+1/2} (\theta_{\mathrm{rain}}-\theta^n) \right)
\theta^{n+1} & = &\theta^{n} + \Delta t \frac{dh^n}{dh^{n+1}} \left( 
G_{\theta}^{(n+1/2)}
+( P^{n+1/2} (\theta_{\mathrm{rain}}-\theta^n) + Q^{n+1/2})/dh^n \right)
\nonumber \\
& & \label{eq:t-n+1-nlfs}
\end{eqnarray}
%-------------------------------------------------------------
Two steps have been added to linear free-surface algorithm 
(eq. \ref{eq:Gv-n-staggered} to \ref{eq:t-n+1-staggered}):
Firstly, the model ``geometry''
(here the {\bf hFacC,W,S}) is updated just before entering {\it
SOLVE\_FOR\_PRESSURE}, using the current $dh^{n}$ field.
Secondly, the vertically integrated continuity equation 
(eq.\ref{eq:h-n+1-nlfs}) has been added ({\bf exactConserv}{\em =TRUE},
in parameter file {\em data}, namelist {\em PARM01})
just before computing the vertical velocity, in subroutine
{\em INTEGR\_CONTINUITY}. This ensures that tracer and continuity equation
discretization  a Although this equation might appear 
redundant with eq.\ref{eq:nstar-nlfs}, the integrated column
thickness $h^{n+1}$ can be different from $\eta^{n+1/2} + H$~:
\begin{itemize}
\item when Crank-Nickelson time-stepping is used (see section
\ref{sect:freesurf-CrankNick}).
\item when filters are applied to the flow field, after
(\ref{eq:v-n+1-nlfs}) and alter the divergence of the flow.
\item when the solver does not iterate until convergence~;
 for example, because a too large residual target was set
 ({\bf cg2dTargetResidual}, parameter file {\em data}, namelist
 {\em PARM02}).
\end{itemize}\noindent
In this staggered time-stepping algorithm, the momentum tendencies
are computed using $dh^{n-1}$ geometry factors.
(eq.\ref{eq:Gv-n-nlfs}) and then rescaled in subroutine {\it TIMESTEP},
(eq.\ref{eq:vstar-nlfs}), similarly to tracer tendencies (see section
\ref{sect:freesurf-tracer-advection}).
The tracers are stepped forward later, using the recently updated
flow field ${\bf v}^{n+1/2}$ and the corresponding model geometry
$dh^{n}$ to compute the tendencies (eq.\ref{eq:Gt-n-nlfs});
Then the tendencies are rescaled by $dh^n/dh^{n+1}$ to derive
the new tracers values $(\theta,S)^{n+1}$ (eq.\ref{eq:t-n+1-nlfs},
in subroutine {\em CALC\_GT, CALC\_GS}).

Note that the fresh-water input is added in a consistent way in the 
continuity equation and in the tracer equation, taking into account
the fresh-water temperature $\theta_{\mathrm{rain}}$.

Regarding the restart procedure,
two 2.D fields $h^{n-1}$ and $(h^n-h^{n-1})/\Delta t$ 
in addition to the standard
state variables and tendencies ($\eta^{n-1/2}$, ${\bf v}^{n-1/2}$,
$\theta^n$, $S^n$, ${\bf G}_{\bf v}^{n-3/2}$, $G_{\theta,S}^{n-1}$)
are stored in a "{\em pickup}" file.
The model restarts reading this "{\em pickup}" file,
then update the model geometry according to $h^{n-1}$,
and compute $h^n$ and the vertical velocity
%$h^n=h^{n-1} + \Delta t [(h^n-h^{n-1})/\Delta t]$,
before starting the main calling sequence (eq.\ref{eq:phi-hyd-nlfs} 
to \ref{eq:t-n+1-nlfs}, {\em S/R FORWARD\_STEP}).
\\

\fbox{ \begin{minipage}{4.75in}
{\em INTEGR\_CONTINUITY} ({\em integr\_continuity.F})

$h^{n+1} -H_o$: {\bf etaH} ({\em DYNVARS.h})

$h^{n} -H_o$: {\bf etaHnm1} ({\em SURFACE.h})

$h^{n+1}-h^{n}/\Delta t$: {\bf dEtaHdt} ({\em SURFACE.h})

\end{minipage} }

\subsubsection{Non-linear free-surface and vertical resolution}
\label{sect:nonlin-freesurf-dz_surf}

When the amplitude of the free-surface variations becomes
as large as the vertical resolution near the surface,
the surface layer thickness can decrease to nearly zero or
can even vanish completely. 
This later possibility has not been implemented, and a 
minimum relative thickness is imposed ({\bf hFacInf}, 
parameter file {\em data}, namelist {\em PARM01}) to prevent 
numerical instabilities caused by very thin surface level.

A better alternative to the vanishing level problem has been 
found and implemented recently, and rely on a different 
vertical coordinate $r^*$~:
The time variation ot the total column thickness becomes
part of the r* coordinate motion, as in a $\sigma_{z},\sigma_{p}$
model, but the fixed part related to topography is treated
as in a height or pressure coordinate model.
A complete description is given in \cite{adcroft:04a}. 

The time-stepping implementation of the $r^*$ coordinate is
identical to the non-linear free-surface in $r$ coordinate,
and differences appear only in the spacial discretization.
\marginpar{needs a subsection ref. here}

1	edhill	1.13	% $Header: /u/gcmpack/manual/part2/nonlin_frsurf.tex,v 1.12 2006/04/04 20:19:00 jmc Exp $
2	jmc	1.1	% $Name: $
3
4	adcroft	1.4
5	jmc	1.1
6	jmc	1.9	\subsection{Non-linear free-surface}
7	adcroft	1.6	\label{sect:nonlinear-freesurface}
8	jmc	1.1
9	jmc	1.8	Recently, new options have been added to the model
10			that concern the free surface formulation.
11	adcroft	1.4
12	jmc	1.1
13	jmc	1.9	\subsubsection{pressure/geo-potential and free surface}
14	jmc	1.11	\label{sect:phi-freesurface}
15	jmc	1.1
16	jmc	1.9	For the atmosphere, since $\phi = \phi_{topo} - \int^p_{p_s} \alpha dp$,
17			subtracting the reference state defined in section
18			\ref{sec:hpe-p-geo-potential-split}~:\\
19			$$
20			\phi_o = \phi_{topo} - \int^p_{p_o} \alpha_o dp
21			\hspace{5mm}\mathrm{with}\hspace{3mm} \phi_o(p_o)=\phi_{topo}
22			$$
23			we get:
24			$$
25			\phi' = \phi - \phi_o = \int^{p_s}_p \alpha dp - \int^{p_o}_p \alpha_o dp
26			$$
27			For the ocean, the reference state is simpler since $\rho_c$ does not dependent
28			on $z$ ($b_o=g$) and the surface reference position is uniformly $z=0$ ($R_o=0$),
29			and the same subtraction leads to a similar relation.
30			For both fluid, using the isomorphic notations, we can write:
31			$$
32			\phi' = \int^{r_{surf}}_r b~ dr - \int^{R_o}_r b_o dr
33			$$
34	edhill	1.13	and re-write as:
35			\begin{equation}
36			\phi' = \int^{r_{surf}}_{R_o} b~ dr + \int^{R_o}_r (b - b_o) dr
37			\label{eq:split-phi-Ro}
38			\end{equation}
39			or:
40			\begin{equation}
41			\phi' = \int^{r_{surf}}_{R_o} b_o dr + \int^{r_{surf}}_r (b - b_o) dr
42	jmc	1.9	\label{eq:split-phi-bo}
43	edhill	1.13	\end{equation}
44	jmc	1.9
45			In section \ref{sec:finding_the_pressure_field}, following eq.\ref{eq:split-phi-Ro},
46			the pressure/geo-potential $\phi'$ has been separated into surface ($\phi_s$),
47			and hydrostatic anomaly ($\phi'_{hyd}$).
48			In this section, the split between $\phi_s$ and $\phi'_{hyd}$ is
49			made according to equation \ref{eq:split-phi-bo}. This slightly
50			different definition reflects the actual implementation in the code
51			and is valid for both linear and non-linear
52			free-surface formulation, in both r-coordinate and r*-coordinate.
53
54			Because the linear free-surface approximation ignore the tracer content
55			of the fluid parcel between $R_o$ and $r_{surf}=R_o+\eta$,
56			for consistency reasons, this part is also neglected in $\phi'_{hyd}$~:
57			$$
58			\phi'_{hyd} = \int^{r_{surf}}_r (b - b_o) dr \simeq \int^{R_o}_r (b - b_o) dr
59			$$
60			Note that in this case, the two definitions of $\phi_s$ and $\phi'_{hyd}$
61			from equation \ref{eq:split-phi-Ro} and \ref{eq:split-phi-bo} converge toward
62			the same (approximated) expressions: $\phi_s = \int^{r_{surf}}_{R_o} b_o dr$
63			and $\phi'_{hyd}=\int^{R_o}_r b' dr$.\\
64			On the contrary, the unapproximated formulation ("non-linear free-surface",
65			see the next section) retains the full expression:
66			$\phi'_{hyd} = \int^{r_{surf}}_r (b - b_o) dr $~.
67			This is obtained by selecting {\bf nonlinFreeSurf}=4 in parameter
68			file {\em data}.\\
69
70			Regarding the surface potential:
71			$$\phi_s = \int_{R_o}^{R_o+\eta} b_o dr = b_s \eta
72			\hspace{5mm}\mathrm{with}\hspace{5mm}
73			b_s = \frac{1}{\eta} \int_{R_o}^{R_o+\eta} b_o dr $$
74			$b_s \simeq b_o(R_o)$ is an excellent approximation (better than
75			the usual numerical truncation, since generally $\|\eta\|$ is smaller
76			than the vertical grid increment).
77	adcroft	1.4
78	jmc	1.9	For the ocean, $\phi_s = g \eta$ and $b_s = g$ is uniform.
79	adcroft	1.4	For the atmosphere, however, because of topographic effects, the
80	jmc	1.9	reference surface pressure $R_o=p_o$ has large spatial variations that
81	adcroft	1.4	are responsible for significant $b_s$ variations (from 0.8 to 1.2
82	jmc	1.9	$[m^3/kg]$). For this reason, when {\bf uniformLin\_PhiSurf} {\em=.FALSE.}
83			(parameter file {\em data}, namelist {\em PARAM01})
84			a non-uniform linear coefficient $b_s$ is used and computed
85			({\it S/R INI\_LINEAR\_PHISURF}) according to the reference surface
86			pressure $p_o$:
87			$b_s = b_o(R_o) = c_p \kappa (p_o / P^o_{SL})^{(\kappa - 1)} \theta_{ref}(p_o)$.
88			with $P^o_{SL}$ the mean sea-level pressure.
89	adcroft	1.4
90	jmc	1.1
91			\subsubsection{Free surface effect on column total thickness
92	jmc	1.9	(Non-linear free-surface)}
93	jmc	1.1
94	adcroft	1.4	The total thickness of the fluid column is $r_{surf} - R_{fixed} =
95	jmc	1.8	\eta + R_o - R_{fixed}$. In most applications, the free surface
96			displacements are small compared to the total thickness
97	jmc	1.9	$\eta \ll H_o = R_o - R_{fixed}$.
98	jmc	1.8	In the previous sections and in older version of the model,
99			the linearized free-surface approximation was made, assuming
100	jmc	1.9	$r_{surf} - R_{fixed} \simeq H_o$ when computing horizontal transports,
101			either in the continuity equation or in tracer and momentum
102	jmc	1.8	advection terms.
103	jmc	1.9	This approximation is dropped when using the non-linear free-surface
104	jmc	1.8	formulation and the total thickness, including the time varying part
105	jmc	1.9	$\eta$, is considered when computing horizontal transports.
106	jmc	1.8	Implications for the barotropic part are presented hereafter.
107	jmc	1.9	In section \ref{sect:freesurf-tracer-advection} consequences for
108			tracer conservation is briefly discussed (more details can be
109			found in \cite{campin:02})~; the general time-stepping is presented
110			in section \ref{sect:nonlin-freesurf-timestepping} with some
111			limitations regarding the vertical resolution in section
112			\ref{sect:nonlin-freesurf-dz_surf}.
113	adcroft	1.4
114	jmc	1.9	In the non-linear formulation, the continuous form of the model
115			equations remains unchanged, except for the 2D continuity equation
116	jmc	1.8	(\ref{eq:discrete-time-backward-free-surface}) which is now
117	adcroft	1.4	integrated from $R_{fixed}(x,y)$ up to $r_{surf}=R_o+\eta$ :
118	jmc	1.1
119			\begin{displaymath}
120			\epsilon_{fs} \partial_t \eta =
121			\left. \dot{r} \right\|_{r=r_{surf}} + \epsilon_{fw} (P-E) =
122	jmc	1.3	- {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o+\eta} \vec{\bf v} dr
123	jmc	1.1	+ \epsilon_{fw} (P-E)
124			\end{displaymath}
125
126	adcroft	1.4	Since $\eta$ has a direct effect on the horizontal velocity (through
127			$\nabla_h \Phi_{surf}$), this adds a non-linear term to the free
128			surface equation. Several options for the time discretization of this
129	jmc	1.8	non-linear part can be considered, as detailed below.
130	adcroft	1.4
131			If the column thickness is evaluated at time step $n$, and with
132			implicit treatment of the surface potential gradient, equations
133			(\ref{eq-solve2D} and \ref{eq-solve2D_rhs}) becomes:
134	jmc	1.1	\begin{eqnarray*}
135			\epsilon_{fs} {\eta}^{n+1} -
136	jmc	1.3	{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n}+R_o-R_{fixed})
137	jmc	1.2	{\bf \nabla}_h b_s {\eta}^{n+1}
138	jmc	1.1	= {\eta}^*
139			\end{eqnarray*}
140			where
141			\begin{eqnarray*}
142			{\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
143	jmc	1.3	\Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o+\eta^n} \vec{\bf v}^* dr
144	jmc	1.1	\: + \: \epsilon_{fw} \Delta_t (P-E)^{n}
145			\end{eqnarray*}
146	adcroft	1.4	This method requires us to update the solver matrix at each time step.
147	jmc	1.1
148			Alternatively, the non-linear contribution can be evaluated fully
149			explicitly:
150			\begin{eqnarray*}
151			\epsilon_{fs} {\eta}^{n+1} -
152	jmc	1.3	{\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{fixed})
153	jmc	1.2	{\bf \nabla}_h b_s {\eta}^{n+1}
154	jmc	1.1	= {\eta}^*
155	jmc	1.2	+{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n})
156			{\bf \nabla}_h b_s {\eta}^{n}
157	jmc	1.1	\end{eqnarray*}
158	adcroft	1.4	This formulation allows one to keep the initial solver matrix
159			unchanged though throughout the integration, since the non-linear free
160			surface only affects the RHS.
161
162			Finally, another option is a "linearized" formulation where the total
163			column thickness appears only in the integral term of the RHS
164			(\ref{eq-solve2D_rhs}) but not directly in the equation
165			(\ref{eq-solve2D}).
166	jmc	1.1
167	jmc	1.9	Those different options (see Table \ref{tab:nonLinFreeSurf_flags})
168			have been tested and show little differences. However, we recommend
169	jmc	1.8	the use of the most precise method (the 1rst one) since the
170	jmc	1.9	computation cost involved in the solver matrix update is negligible.
171	jmc	1.8
172	jmc	1.9	\begin{table}[htb]
173			%\begin{center}
174			\centering
175	jmc	1.8	\begin{tabular}[htb]{\|l\|c\|l\|}
176			\hline
177			parameter & value & description \\
178			\hline
179			& -1 & linear free-surface, restart from a pickup file \\
180			& & produced with \#undef EXACT\_CONSERV code\\
181			\cline{2-3}
182	jmc	1.9	& 0 & Linear free-surface \\
183	jmc	1.8	\cline{2-3}
184			nonlinFreeSurf & 4 & Non-linear free-surface \\
185			\cline{2-3}
186			& 3 & same as 4 but neglecting
187			$\int_{R_o}^{R_o+\eta} b' dr $ in $\Phi'_{hyd}$ \\
188			\cline{2-3}
189			& 2 & same as 3 but do not update cg2d solver matrix \\
190			\cline{2-3}
191			& 1 & same as 2 but treat momentum as in Linear FS \\
192			\hline
193			& 0 & do not use $r*$ vertical coordinate (= default)\\
194			\cline{2-3}
195			select\_rStar & 2 & use $r^*$ vertical coordinate \\
196			\cline{2-3}
197			& 1 & same as 2 but without the contribution of the\\
198			& & slope of the coordinate in $\nabla \Phi$ \\
199			\hline
200			\end{tabular}
201	jmc	1.9	\caption{Non-linear free-surface flags}
202			\label{tab:nonLinFreeSurf_flags}
203			%\end{center}
204			\end{table}
205	jmc	1.8
206	jmc	1.1
207	jmc	1.9	\subsubsection{Tracer conservation with non-linear free-surface}
208	adcroft	1.4	\label{sect:freesurf-tracer-advection}
209	jmc	1.1
210	adcroft	1.4	To ensure global tracer conservation (i.e., the total amount) as well
211			as local conservation, the change in the surface level thickness must
212			be consistent with the way the continuity equation is integrated, both
213			in the barotropic part (to find $\eta$) and baroclinic part (to find
214			$w = \dot{r}$).
215	jmc	1.1
216	jmc	1.9	To illustrate this, consider the shallow water model, with a source
217			of fresh water (P):
218	jmc	1.1	$$
219	jmc	1.9	\partial_t h + \nabla \cdot h \vec{\bf v} = P
220	jmc	1.1	$$
221			where $h$ is the total thickness of the water column.
222			To conserve the tracer $\theta$ we have to discretize:
223			$$
224	jmc	1.9	\partial_t (h \theta) + \nabla \cdot ( h \theta \vec{\bf v})
225			= P \theta_{\mathrm{rain}}
226	jmc	1.1	$$
227	adcroft	1.4	Using the implicit (non-linear) free surface described above (section
228	adcroft	1.6	\ref{sect:pressure-method-linear-backward}) we have:
229	jmc	1.1	\begin{eqnarray*}
230	jmc	1.9	h^{n+1} = h^{n} - \Delta t \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) + \Delta t P \\
231	jmc	1.1	\end{eqnarray*}
232	adcroft	1.4	The discretized form of the tracer equation must adopt the same
233			``form'' in the computation of tracer fluxes, that is, the same value
234			of $h$, as used in the continuity equation:
235	jmc	1.1	\begin{eqnarray*}
236			h^{n+1} \, \theta^{n+1} = h^n \, \theta^n
237	jmc	1.9	- \Delta t \nabla \cdot (h^n \, \theta^n \, \vec{\bf v}^{n+1})
238			+ \Delta t P \theta_{rain}
239	jmc	1.1	\end{eqnarray*}
240
241	jmc	1.9	The use of a 3 time-levels time-stepping scheme such as the Adams-Bashforth
242			make the conservation sightly tricky.
243	jmc	1.8	The current implementation with the Adams-Bashforth time-stepping
244			provides an exact local conservation and prevents any drift in
245			the global tracer content (\cite{campin:02}).
246			Compared to the linear free-surface method, an additional step is required:
247			the variation of the water column thickness (from $h^n$ to $h^{n+1}$) is
248	adcroft	1.4	not incorporated directly into the tracer equation. Instead, the
249			model uses the $G_\theta$ terms (first step) as in the linear free
250			surface formulation (with the "{\it surface correction}" turned "on",
251			see tracer section):
252	jmc	1.1	$$
253	jmc	1.2	G_\theta^n = \left(- \nabla \cdot (h^n \, \theta^n \, \vec{\bf v}^{n+1})
254			- \dot{r}_{surf}^{n+1} \theta^n \right) / h^n
255	jmc	1.1	$$
256	adcroft	1.4	Then, in a second step, the thickness variation (expansion/reduction)
257			is taken into account:
258	jmc	1.1	$$
259	jmc	1.9	\theta^{n+1} = \theta^n + \Delta t \frac{h^n}{h^{n+1}}
260			\left( G_\theta^{(n+1/2)} + P (\theta_{\mathrm{rain}} - \theta^n )/h^n \right)
261			%\theta^{n+1} = \theta^n + \frac{\Delta t}{h^{n+1}}
262			% \left( h^n G_\theta^{(n+1/2)} + P (\theta_{\mathrm{rain}} - \theta^n ) \right)
263	jmc	1.1	$$
264			Note that with a simple forward time step (no Adams-Bashforth),
265	jmc	1.9	these two formulations are equivalent,
266	jmc	1.1	since
267			$
268	jmc	1.9	(h^{n+1} - h^{n})/ \Delta t =
269			P - \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) = P + \dot{r}_{surf}^{n+1}
270	jmc	1.1	$
271	adcroft	1.4
272	jmc	1.9	\subsubsection{Time stepping implementation of the
273			non-linear free-surface}
274			\label{sect:nonlin-freesurf-timestepping}
275
276			The grid cell thickness was hold constant with the linear
277			free-surface~; with the non-linear free-surface, it is now varying
278			in time, at least at the surface level.
279			This implies some modifications of the general algorithm described
280			earlier in sections \ref{sect:adams-bashforth-sync} and
281			\ref{sect:adams-bashforth-staggered}.
282
283			A simplified version of the staggered in time, non-linear
284			free-surface algorithm is detailed hereafter, and can be compared
285			to the equivalent linear free-surface case (eq. \ref{eq:Gv-n-staggered}
286			to \ref{eq:t-n+1-staggered}) and can also be easily transposed
287			to the synchronous time-stepping case.
288			Among the simplifications, salinity equation, implicit operator
289			and detailed elliptic equation are omitted. Surface forcing is
290			explicitly written as fluxes of temperature, fresh water and
291			momentum, $Q^{n+1/2}, P^{n+1/2}, F_{\bf v}^n$ respectively.
292			$h^n$ and $dh^n$ are the column and grid box thickness in r-coordinate.
293			%-------------------------------------------------------------
294			\begin{eqnarray}
295			\phi^{n}_{hyd} & = & \int b(\theta^{n},S^{n},r) dr
296			\label{eq:phi-hyd-nlfs} \\
297			\vec{\bf G}_{\vec{\bf v}}^{n-1/2}\hspace{-2mm} & = &
298			\vec{\bf G}_{\vec{\bf v}} (dh^{n-1},\vec{\bf v}^{n-1/2})
299			\hspace{+2mm};\hspace{+2mm}
300			\vec{\bf G}_{\vec{\bf v}}^{(n)} =
301			\frac{3}{2} \vec{\bf G}_{\vec{\bf v}}^{n-1/2}
302			- \frac{1}{2} \vec{\bf G}_{\vec{\bf v}}^{n-3/2}
303			\label{eq:Gv-n-nlfs} \\
304			%\vec{\bf G}_{\vec{\bf v}}^{(n)} & = &
305			% \frac{3}{2} \vec{\bf G}_{\vec{\bf v}}^{n-1/2}
306			%- \frac{1}{2} \vec{\bf G}_{\vec{\bf v}}^{n-3/2}
307			%\label{eq:Gv-n+5-nlfs} \\
308			%\vec{\bf v}^{*} & = & \vec{\bf v}^{n-1/2} + \frac{\Delta t}{dh^{n}} \left(
309			%dh^{n-1}\vec{\bf G}_{\vec{\bf v}}^{(n)} + F_{\vec{\bf v}}^{n} \right)
310			\vec{\bf v}^{*} & = & \vec{\bf v}^{n-1/2} + \Delta t \frac{dh^{n-1}}{dh^{n}} \left(
311			\vec{\bf G}_{\vec{\bf v}}^{(n)} + F_{\vec{\bf v}}^{n}/dh^{n-1} \right)
312			- \Delta t \nabla \phi_{hyd}^{n}
313	jmc	1.12	\label{eq:vstar-nlfs}
314			\end{eqnarray}
315			\hspace{3cm}$\longrightarrow$~~{\it update model~geometry~:~}${\bf hFac}(dh^n)$\\
316			\begin{eqnarray}
317	jmc	1.9	\eta^{n+1/2} \hspace{-2mm} & = &
318			\eta^{n-1/2} + \Delta t P^{n+1/2} - \Delta t
319			\nabla \cdot \int \vec{\bf v}^{n+1/2} dh^{n} \nonumber \\
320			& = & \eta^{n-1/2} + \Delta t P^{n+1/2} - \Delta t
321			\nabla \cdot \int \!\!\! \left( \vec{\bf v}^* - g \Delta t \nabla \eta^{n+1/2} \right) dh^{n}
322			\label{eq:nstar-nlfs} \\
323			\vec{\bf v}^{n+1/2}\hspace{-2mm} & = &
324			\vec{\bf v}^{*} - g \Delta t \nabla \eta^{n+1/2}
325			\label{eq:v-n+1-nlfs} \\
326			h^{n+1} & = & h^{n} + \Delta t P^{n+1/2} - \Delta t
327			\nabla \cdot \int \vec{\bf v}^{n+1/2} dh^{n}
328			\label{eq:h-n+1-nlfs} \\
329			G_{\theta}^{n} & = & G_{\theta} ( dh^{n}, u^{n+1/2}, \theta^{n} )
330			\hspace{+2mm};\hspace{+2mm}
331			G_{\theta}^{(n+1/2)} = \frac{3}{2} G_{\theta}^{n} - \frac{1}{2} G_{\theta}^{n-1}
332			\label{eq:Gt-n-nlfs} \\
333			%\theta^{n+1} & = &\theta^{n} + \frac{\Delta t}{dh^{n+1}} \left( dh^n
334			%G_{\theta}^{(n+1/2)} + Q^{n+1/2} + P^{n+1/2} (\theta_{\mathrm{rain}}-\theta^n) \right)
335			\theta^{n+1} & = &\theta^{n} + \Delta t \frac{dh^n}{dh^{n+1}} \left(
336			G_{\theta}^{(n+1/2)}
337			+( P^{n+1/2} (\theta_{\mathrm{rain}}-\theta^n) + Q^{n+1/2})/dh^n \right)
338			\nonumber \\
339			& & \label{eq:t-n+1-nlfs}
340			\end{eqnarray}
341			%-------------------------------------------------------------
342			Two steps have been added to linear free-surface algorithm
343			(eq. \ref{eq:Gv-n-staggered} to \ref{eq:t-n+1-staggered}):
344			Firstly, the model ``geometry''
345	adcroft	1.4	(here the {\bf hFacC,W,S}) is updated just before entering {\it
346	jmc	1.9	SOLVE\_FOR\_PRESSURE}, using the current $dh^{n}$ field.
347			Secondly, the vertically integrated continuity equation
348			(eq.\ref{eq:h-n+1-nlfs}) has been added ({\bf exactConserv}{\em =TRUE},
349			in parameter file {\em data}, namelist {\em PARM01})
350			just before computing the vertical velocity, in subroutine
351			{\em INTEGR\_CONTINUITY}. This ensures that tracer and continuity equation
352			discretization a Although this equation might appear
353			redundant with eq.\ref{eq:nstar-nlfs}, the integrated column
354			thickness $h^{n+1}$ can be different from $\eta^{n+1/2} + H$~:
355			\begin{itemize}
356			\item when Crank-Nickelson time-stepping is used (see section
357			\ref{sect:freesurf-CrankNick}).
358			\item when filters are applied to the flow field, after
359			(\ref{eq:v-n+1-nlfs}) and alter the divergence of the flow.
360			\item when the solver does not iterate until convergence~;
361			for example, because a too large residual target was set
362			({\bf cg2dTargetResidual}, parameter file {\em data}, namelist
363			{\em PARM02}).
364			\end{itemize}\noindent
365			In this staggered time-stepping algorithm, the momentum tendencies
366			are computed using $dh^{n-1}$ geometry factors.
367			(eq.\ref{eq:Gv-n-nlfs}) and then rescaled in subroutine {\it TIMESTEP},
368			(eq.\ref{eq:vstar-nlfs}), similarly to tracer tendencies (see section
369			\ref{sect:freesurf-tracer-advection}).
370			The tracers are stepped forward later, using the recently updated
371			flow field ${\bf v}^{n+1/2}$ and the corresponding model geometry
372			$dh^{n}$ to compute the tendencies (eq.\ref{eq:Gt-n-nlfs});
373			Then the tendencies are rescaled by $dh^n/dh^{n+1}$ to derive
374			the new tracers values $(\theta,S)^{n+1}$ (eq.\ref{eq:t-n+1-nlfs},
375			in subroutine {\em CALC\_GT, CALC\_GS}).
376
377			Note that the fresh-water input is added in a consistent way in the
378			continuity equation and in the tracer equation, taking into account
379			the fresh-water temperature $\theta_{\mathrm{rain}}$.
380
381			Regarding the restart procedure,
382			two 2.D fields $h^{n-1}$ and $(h^n-h^{n-1})/\Delta t$
383			in addition to the standard
384			state variables and tendencies ($\eta^{n-1/2}$, ${\bf v}^{n-1/2}$,
385			$\theta^n$, $S^n$, ${\bf G}_{\bf v}^{n-3/2}$, $G_{\theta,S}^{n-1}$)
386			are stored in a "{\em pickup}" file.
387			The model restarts reading this "{\em pickup}" file,
388			then update the model geometry according to $h^{n-1}$,
389			and compute $h^n$ and the vertical velocity
390			%$h^n=h^{n-1} + \Delta t [(h^n-h^{n-1})/\Delta t]$,
391			before starting the main calling sequence (eq.\ref{eq:phi-hyd-nlfs}
392			to \ref{eq:t-n+1-nlfs}, {\em S/R FORWARD\_STEP}).
393			\\
394
395			\fbox{ \begin{minipage}{4.75in}
396			{\em INTEGR\_CONTINUITY} ({\em integr\_continuity.F})
397
398			$h^{n+1} -H_o$: {\bf etaH} ({\em DYNVARS.h})
399
400			$h^{n} -H_o$: {\bf etaHnm1} ({\em SURFACE.h})
401
402			$h^{n+1}-h^{n}/\Delta t$: {\bf dEtaHdt} ({\em SURFACE.h})
403
404			\end{minipage} }
405	jmc	1.1
406	jmc	1.9	\subsubsection{Non-linear free-surface and vertical resolution}
407	jmc	1.8	\label{sect:nonlin-freesurf-dz_surf}
408
409			When the amplitude of the free-surface variations becomes
410			as large as the vertical resolution near the surface,
411			the surface layer thickness can decrease to nearly zero or
412	jmc	1.9	can even vanish completely.
413	jmc	1.8	This later possibility has not been implemented, and a
414			minimum relative thickness is imposed ({\bf hFacInf},
415			parameter file {\em data}, namelist {\em PARM01}) to prevent
416			numerical instabilities caused by very thin surface level.
417
418	jmc	1.9	A better alternative to the vanishing level problem has been
419	jmc	1.8	found and implemented recently, and rely on a different
420			vertical coordinate $r^*$~:
421			The time variation ot the total column thickness becomes
422			part of the r* coordinate motion, as in a $\sigma_{z},\sigma_{p}$
423			model, but the fixed part related to topography is treated
424			as in a height or pressure coordinate model.
425	jmc	1.10	A complete description is given in \cite{adcroft:04a}.
426	jmc	1.8
427			The time-stepping implementation of the $r^*$ coordinate is
428			identical to the non-linear free-surface in $r$ coordinate,
429	jmc	1.9	and differences appear only in the spacial discretization.
430	jmc	1.8	\marginpar{needs a subsection ref. here}
431	jmc	1.1