s_algorithm/text/nonlin_frsurf.tex

% $Header: /u/gcmpack/manual/part2/nonlin_frsurf.tex,v 1.7 2001/11/13 20:27:54 adcroft 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{Non-uniform linear-relation for the surface potential}

The linear relation between surface pressure/geo-potential
($\Phi_{surf}$) and surface displacement ($\eta$) could be considered
to be a constant ($b_s=$ constant)
\marginpar{add a reference to part.1 here}
but is in fact dependent on the position ($x,y,r$)
since we linearize:
$$\Phi_{surf}=\int_{R_o}^{R_o+\eta} b dr \simeq b_s \eta
~\mathrm{with}~ b_s = b(\theta,S,r)_{r=R_o} 
\simeq b_s(\theta_{ref}(R_o),S_{ref}(R_o),R_o)$$
Note that, for convenience, the effect on $b_s$ of the local surface
$\theta,S$ are not considered here, but are incorporated in to
$\Phi'_{hyd}$.

For the ocean, $b_s = g \rho_{surf} / \rho_o \simeq g$ is a very good
approximation since the relative difference in surface density are
usually small and only due to local $\theta,S$ gradients (because the
upper surface, $R_o = 0$, is essentially flat). Therefore, they can
easily be incorporated in $\Phi'_{hyd}$.

For the atmosphere, however, because of topographic effects, the
reference surface pressure $R_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, we use a non-uniform linear coefficient
$b_s$.

In practice, in an oceanic configuration or when the default value
(TRUE) of the parameter {\bf uniformLin\_PhiSurf} is used, then $b_s$
is simply set to $g$ for the ocean and $1.$ for the atmosphere.
Turning {\bf uniformLin\_PhiSurf} to "FALSE", tells the code to
evaluate $b_s$ from the reference vertical profile $\theta_{ref}$
({\it S/R INI\_LINEAR\_PHISURF}) according to the reference surface
pressure $P_o$ ($=R_o$): $b_s = c_p \kappa (P_o / Pc)^{(\kappa - 1)}
\theta_{ref}(P_o)$


\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 << 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 the horizontal transport is 
computed, 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 consisdered when computing horizontal transport.
Implications for the barotropic part are presented hereafter.
In sections \ref{sect:freesurf-tracer-advection} and
\ref{sect:freesurf-momentum-advection}, consequences for tracer 
and momentum are brifly discussed. a more detailed description
is available in \cite{campin:02}.


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 tab.?? for the one still available)
have been tested and show litle differences. However, we recommand
the use of the most precise method (the 1rst one) since the 
computation cost involved in the solver matrix update are negligeable.

\begin{center}
 \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}
\end{center}


\subsubsection{Free surface effect on the surface level thickness
(Non-linear free surface): Tracer advection}
\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 uniform
Cartesian horizontal grid:
$$
\partial_t h + \nabla \cdot h \vec{\bf v} = 0
$$
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})= 0
$$
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} ) \\
\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})
\end{eqnarray*}

The use of a 3 time-levels timestepping scheme such as the Adams-Bashforth
make the conservation less straitforward.
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}} G_\theta^{(n+1/2)} 
$$
Note that with a simple forward time step (no Adams-Bashforth), 
since
$
\dot{r}_{surf}^{n+1} 
= - \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) = (h^{n+1} - h^{n})
/ \Delta_t
$
these two formulations are equivalent. 

Implementation in the MITgcm is as follows.  The model ``geometry''
(here the {\bf hFacC,W,S}) is updated just before entering {\it
SOLVE\_FOR\_PRESSURE}, using the current $\eta$ field.  Then, at the
end of the time step, the variables are advanced in time, so that
$\eta^n$ replaces $\eta^{n-1}$.  At the next time step, the tracer
tendency ($G_\theta$) is computed using the same geometry, which is
now consistent with $\eta^{n-1}$.  Finally, in S/R {\it
TIMESTEP\_TRACER}, the expansion/reduction ratio is applied to the
surface level to compute the new tracer field.


\subsubsection{Free surface effect on the surface level thickness
(Non-linear free surface): Momentum advection}     
\label{sect:freesurf-momentum-advection}

With the flux form formulation, advection of momentum
can be treated exactly as the tracer advection is.
Here the expansion/reduction factors ($hFacW^{n+1}/hFacW^n$ for $u$
and $hFacS^{n+1}/hFacS^n$ for $v$) are simply applied in the
subroutine {\it TIMESTEP}.

Regarding momentum advection,
the vector invariant formulation is similar to the
advective form (as opposed to the flux form) and therefore
does not need specific modification to include the 
free surface effect on the surface level thickness.
Updating the {\bf hFacC,W,S} and the {\bf recip\_hFac}(s) 
at one given place (like describe before) is sufficient.

\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 vanishe completly. 
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 atlternative 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:04}. 

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 discretisation.
\marginpar{needs a subsection ref. here}

1	jmc	1.8	% $Header: /u/gcmpack/manual/part2/nonlin_frsurf.tex,v 1.7 2001/11/13 20:27:54 adcroft Exp $
2	jmc	1.1	% $Name: $
3
4	adcroft	1.4
5	jmc	1.1
6			\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			\subsubsection{Non-uniform linear-relation for the surface potential}
14
15	adcroft	1.4	The linear relation between surface pressure/geo-potential
16			($\Phi_{surf}$) and surface displacement ($\eta$) could be considered
17			to be a constant ($b_s=$ constant)
18	jmc	1.3	\marginpar{add a reference to part.1 here}
19			but is in fact dependent on the position ($x,y,r$)
20	jmc	1.2	since we linearize:
21			$$\Phi_{surf}=\int_{R_o}^{R_o+\eta} b dr \simeq b_s \eta
22			~\mathrm{with}~ b_s = b(\theta,S,r)_{r=R_o}
23			\simeq b_s(\theta_{ref}(R_o),S_{ref}(R_o),R_o)$$
24	adcroft	1.4	Note that, for convenience, the effect on $b_s$ of the local surface
25			$\theta,S$ are not considered here, but are incorporated in to
26			$\Phi'_{hyd}$.
27
28			For the ocean, $b_s = g \rho_{surf} / \rho_o \simeq g$ is a very good
29			approximation since the relative difference in surface density are
30			usually small and only due to local $\theta,S$ gradients (because the
31			upper surface, $R_o = 0$, is essentially flat). Therefore, they can
32			easily be incorporated in $\Phi'_{hyd}$.
33
34			For the atmosphere, however, because of topographic effects, the
35	cnh	1.5	reference surface pressure $R_o$ has large spatial variations that
36	adcroft	1.4	are responsible for significant $b_s$ variations (from 0.8 to 1.2
37			$[m^3/kg]$). For this reason, we use a non-uniform linear coefficient
38			$b_s$.
39
40			In practice, in an oceanic configuration or when the default value
41			(TRUE) of the parameter {\bf uniformLin\_PhiSurf} is used, then $b_s$
42			is simply set to $g$ for the ocean and $1.$ for the atmosphere.
43			Turning {\bf uniformLin\_PhiSurf} to "FALSE", tells the code to
44			evaluate $b_s$ from the reference vertical profile $\theta_{ref}$
45			({\it S/R INI\_LINEAR\_PHISURF}) according to the reference surface
46			pressure $P_o$ ($=R_o$): $b_s = c_p \kappa (P_o / Pc)^{(\kappa - 1)}
47			\theta_{ref}(P_o)$
48
49	jmc	1.1
50			\subsubsection{Free surface effect on column total thickness
51			(Non-linear free surface)}
52
53	adcroft	1.4	The total thickness of the fluid column is $r_{surf} - R_{fixed} =
54	jmc	1.8	\eta + R_o - R_{fixed}$. In most applications, the free surface
55			displacements are small compared to the total thickness
56			$\eta << H_o = R_o - R_{fixed}$.
57			In the previous sections and in older version of the model,
58			the linearized free-surface approximation was made, assuming
59			$r_{surf} - R_{fixed} \simeq H_o$ when the horizontal transport is
60			computed, either in the continuity equation or in tracer and momentum
61			advection terms.
62	adcroft	1.4	This approximation is dropped when using the non-linear free surface
63	jmc	1.8	formulation and the total thickness, including the time varying part
64			$\eta$, is consisdered when computing horizontal transport.
65			Implications for the barotropic part are presented hereafter.
66			In sections \ref{sect:freesurf-tracer-advection} and
67			\ref{sect:freesurf-momentum-advection}, consequences for tracer
68			and momentum are brifly discussed. a more detailed description
69			is available in \cite{campin:02}.
70	adcroft	1.4
71
72	jmc	1.8	In the non-linear formulation, the continuous form of the model equations
73			remains unchanged, except for the 2D continuity equation
74			(\ref{eq:discrete-time-backward-free-surface}) which is now
75	adcroft	1.4	integrated from $R_{fixed}(x,y)$ up to $r_{surf}=R_o+\eta$ :
76	jmc	1.1
77			\begin{displaymath}
78			\epsilon_{fs} \partial_t \eta =
79			\left. \dot{r} \right\|_{r=r_{surf}} + \epsilon_{fw} (P-E) =
80	jmc	1.3	- {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o+\eta} \vec{\bf v} dr
81	jmc	1.1	+ \epsilon_{fw} (P-E)
82			\end{displaymath}
83
84	adcroft	1.4	Since $\eta$ has a direct effect on the horizontal velocity (through
85			$\nabla_h \Phi_{surf}$), this adds a non-linear term to the free
86			surface equation. Several options for the time discretization of this
87	jmc	1.8	non-linear part can be considered, as detailed below.
88	adcroft	1.4
89			If the column thickness is evaluated at time step $n$, and with
90			implicit treatment of the surface potential gradient, equations
91			(\ref{eq-solve2D} and \ref{eq-solve2D_rhs}) becomes:
92	jmc	1.1	\begin{eqnarray*}
93			\epsilon_{fs} {\eta}^{n+1} -
94	jmc	1.3	{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n}+R_o-R_{fixed})
95	jmc	1.2	{\bf \nabla}_h b_s {\eta}^{n+1}
96	jmc	1.1	= {\eta}^*
97			\end{eqnarray*}
98			where
99			\begin{eqnarray*}
100			{\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
101	jmc	1.3	\Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o+\eta^n} \vec{\bf v}^* dr
102	jmc	1.1	\: + \: \epsilon_{fw} \Delta_t (P-E)^{n}
103			\end{eqnarray*}
104	adcroft	1.4	This method requires us to update the solver matrix at each time step.
105	jmc	1.1
106			Alternatively, the non-linear contribution can be evaluated fully
107			explicitly:
108			\begin{eqnarray*}
109			\epsilon_{fs} {\eta}^{n+1} -
110	jmc	1.3	{\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{fixed})
111	jmc	1.2	{\bf \nabla}_h b_s {\eta}^{n+1}
112	jmc	1.1	= {\eta}^*
113	jmc	1.2	+{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n})
114			{\bf \nabla}_h b_s {\eta}^{n}
115	jmc	1.1	\end{eqnarray*}
116	adcroft	1.4	This formulation allows one to keep the initial solver matrix
117			unchanged though throughout the integration, since the non-linear free
118			surface only affects the RHS.
119
120			Finally, another option is a "linearized" formulation where the total
121			column thickness appears only in the integral term of the RHS
122			(\ref{eq-solve2D_rhs}) but not directly in the equation
123			(\ref{eq-solve2D}).
124	jmc	1.1
125	jmc	1.8	Those different options (see tab.?? for the one still available)
126			have been tested and show litle differences. However, we recommand
127			the use of the most precise method (the 1rst one) since the
128			computation cost involved in the solver matrix update are negligeable.
129
130			\begin{center}
131			\begin{tabular}[htb]{\|l\|c\|l\|}
132			\hline
133			parameter & value & description \\
134			\hline
135			& -1 & linear free-surface, restart from a pickup file \\
136			& & produced with \#undef EXACT\_CONSERV code\\
137			\cline{2-3}
138			& 0 & Linear free-Surface \\
139			\cline{2-3}
140			nonlinFreeSurf & 4 & Non-linear free-surface \\
141			\cline{2-3}
142			& 3 & same as 4 but neglecting
143			$\int_{R_o}^{R_o+\eta} b' dr $ in $\Phi'_{hyd}$ \\
144			\cline{2-3}
145			& 2 & same as 3 but do not update cg2d solver matrix \\
146			\cline{2-3}
147			& 1 & same as 2 but treat momentum as in Linear FS \\
148			\hline
149			& 0 & do not use $r*$ vertical coordinate (= default)\\
150			\cline{2-3}
151			select\_rStar & 2 & use $r^*$ vertical coordinate \\
152			\cline{2-3}
153			& 1 & same as 2 but without the contribution of the\\
154			& & slope of the coordinate in $\nabla \Phi$ \\
155			\hline
156			\end{tabular}
157			\end{center}
158
159	jmc	1.1
160			\subsubsection{Free surface effect on the surface level thickness
161			(Non-linear free surface): Tracer advection}
162	adcroft	1.4	\label{sect:freesurf-tracer-advection}
163	jmc	1.1
164	adcroft	1.4	To ensure global tracer conservation (i.e., the total amount) as well
165			as local conservation, the change in the surface level thickness must
166			be consistent with the way the continuity equation is integrated, both
167			in the barotropic part (to find $\eta$) and baroclinic part (to find
168			$w = \dot{r}$).
169	jmc	1.1
170	adcroft	1.4	To illustrate this, consider the shallow water model, with uniform
171			Cartesian horizontal grid:
172	jmc	1.1	$$
173			\partial_t h + \nabla \cdot h \vec{\bf v} = 0
174			$$
175			where $h$ is the total thickness of the water column.
176			To conserve the tracer $\theta$ we have to discretize:
177			$$
178			\partial_t (h \theta) + \nabla \cdot ( h \theta \vec{\bf v})= 0
179			$$
180	adcroft	1.4	Using the implicit (non-linear) free surface described above (section
181	adcroft	1.6	\ref{sect:pressure-method-linear-backward}) we have:
182	jmc	1.1	\begin{eqnarray*}
183	jmc	1.2	h^{n+1} = h^{n} - \Delta_t \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) \\
184	jmc	1.1	\end{eqnarray*}
185	adcroft	1.4	The discretized form of the tracer equation must adopt the same
186			``form'' in the computation of tracer fluxes, that is, the same value
187			of $h$, as used in the continuity equation:
188	jmc	1.1	\begin{eqnarray*}
189			h^{n+1} \, \theta^{n+1} = h^n \, \theta^n
190	jmc	1.2	- \Delta_t \nabla \cdot (h^n \, \theta^n \, \vec{\bf v}^{n+1})
191	jmc	1.1	\end{eqnarray*}
192
193	jmc	1.8	The use of a 3 time-levels timestepping scheme such as the Adams-Bashforth
194			make the conservation less straitforward.
195			The current implementation with the Adams-Bashforth time-stepping
196			provides an exact local conservation and prevents any drift in
197			the global tracer content (\cite{campin:02}).
198			Compared to the linear free-surface method, an additional step is required:
199			the variation of the water column thickness (from $h^n$ to $h^{n+1}$) is
200	adcroft	1.4	not incorporated directly into the tracer equation. Instead, the
201			model uses the $G_\theta$ terms (first step) as in the linear free
202			surface formulation (with the "{\it surface correction}" turned "on",
203			see tracer section):
204	jmc	1.1	$$
205	jmc	1.2	G_\theta^n = \left(- \nabla \cdot (h^n \, \theta^n \, \vec{\bf v}^{n+1})
206			- \dot{r}_{surf}^{n+1} \theta^n \right) / h^n
207	jmc	1.1	$$
208	adcroft	1.4	Then, in a second step, the thickness variation (expansion/reduction)
209			is taken into account:
210	jmc	1.1	$$
211	jmc	1.2	\theta^{n+1} = \theta^n + \Delta_t \frac{h^n}{h^{n+1}} G_\theta^{(n+1/2)}
212	jmc	1.1	$$
213			Note that with a simple forward time step (no Adams-Bashforth),
214			since
215			$
216	jmc	1.2	\dot{r}_{surf}^{n+1}
217			= - \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) = (h^{n+1} - h^{n})
218	jmc	1.1	/ \Delta_t
219			$
220	adcroft	1.4	these two formulations are equivalent.
221
222			Implementation in the MITgcm is as follows. The model ``geometry''
223			(here the {\bf hFacC,W,S}) is updated just before entering {\it
224			SOLVE\_FOR\_PRESSURE}, using the current $\eta$ field. Then, at the
225			end of the time step, the variables are advanced in time, so that
226			$\eta^n$ replaces $\eta^{n-1}$. At the next time step, the tracer
227			tendency ($G_\theta$) is computed using the same geometry, which is
228			now consistent with $\eta^{n-1}$. Finally, in S/R {\it
229			TIMESTEP\_TRACER}, the expansion/reduction ratio is applied to the
230			surface level to compute the new tracer field.
231	jmc	1.1
232
233			\subsubsection{Free surface effect on the surface level thickness
234			(Non-linear free surface): Momentum advection}
235	adcroft	1.4	\label{sect:freesurf-momentum-advection}
236	jmc	1.1
237	jmc	1.8	With the flux form formulation, advection of momentum
238			can be treated exactly as the tracer advection is.
239			Here the expansion/reduction factors ($hFacW^{n+1}/hFacW^n$ for $u$
240			and $hFacS^{n+1}/hFacS^n$ for $v$) are simply applied in the
241			subroutine {\it TIMESTEP}.
242
243	jmc	1.1	Regarding momentum advection,
244			the vector invariant formulation is similar to the
245			advective form (as opposed to the flux form) and therefore
246			does not need specific modification to include the
247			free surface effect on the surface level thickness.
248			Updating the {\bf hFacC,W,S} and the {\bf recip\_hFac}(s)
249			at one given place (like describe before) is sufficient.
250
251	jmc	1.8	\subsubsection{Non-linear free surface and vertical resolution}
252			\label{sect:nonlin-freesurf-dz_surf}
253
254			When the amplitude of the free-surface variations becomes
255			as large as the vertical resolution near the surface,
256			the surface layer thickness can decrease to nearly zero or
257			can even vanishe completly.
258			This later possibility has not been implemented, and a
259			minimum relative thickness is imposed ({\bf hFacInf},
260			parameter file {\em data}, namelist {\em PARM01}) to prevent
261			numerical instabilities caused by very thin surface level.
262
263			A better atlternative to the vanishing level problem has been
264			found and implemented recently, and rely on a different
265			vertical coordinate $r^*$~:
266			The time variation ot the total column thickness becomes
267			part of the r* coordinate motion, as in a $\sigma_{z},\sigma_{p}$
268			model, but the fixed part related to topography is treated
269			as in a height or pressure coordinate model.
270			A complete description is given in \cite{adcroft:04}.
271
272			The time-stepping implementation of the $r^*$ coordinate is
273			identical to the non-linear free-surface in $r$ coordinate,
274			and differences appear only in the spacial discretisation.
275			\marginpar{needs a subsection ref. here}
276	jmc	1.1