s_algorithm/text/nonlin_frsurf.tex

% $Header: /u/gcmpack/mitgcmdoc/part2/nonlin_frsurf.tex,v 1.2 2001/08/17 18:38:10 jmc Exp $
% $Name:  $

%\section{Time Integration}

\subsection{Non-linear free surface}

Recently, 2 options have added to the model
(and therefore, have not yet been extensively tested)
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$)
has been considered as uniform ($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 convinience, the effect of the local 
surface $\theta,S$ on $b_s$
are not considered here, but incorporated in $\Phi'_{hyd}$.

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

For the atmosphere, because of topographic effects,
the reference surface pressure $R_o$ has large spacial differences
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$.

Practically, 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", allows 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 the linear free surface approximation
(detailed before), only the fixed part of
it ($R_o - R_{fixed})$ is considered when we integrate the 
continuity equation or compute tracer and momentum advection term.

This approximation is dropped when using 
the non-linear free surface formulation. 
Details follow here after for the barotropic part
and in the 2 next subsections for the baroclinic
part.

%------------------------------------------
% Non-Linear Barotropic part

The continuous form of the model equations remains 
unchanged, except for the 2D continuity equation
(\ref{eq-tCsC-eta}) that 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.

Regarding the time discretization of this non-linear part,
several options are now being tested:

With the column thickness evaluated at time step $n$,
and the surface potential gradient still implicit,
equation (\ref{eq-solve2D} \& \ref{eq-solve2D_rhs})
become:
\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}^*
%\label{solve_2d}
\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}
%\label{solve_2d_rhs}
\end{eqnarray*} 
This method requires 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 to keep the initial solver matrix
since the non-linear free surface only affects the RHS.

An other 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}).

%------------------------------------------
\subsubsection{Free surface effect on the surface level thickness
(Non-linear free surface): Tracer advection}

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

To illustrate this, let's 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 before, 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 use the same
"geometry" to compute the tracer fluxes, that is, the same value of
h, as the continuity equation did:
\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*}

In order to deal with the Adams-Bashforth time stepping,
we implement this scheme slightly differently, in two step:
the variation of the water column thickness (from
$h^n$ to $h^{n+1}$)
is not incorporated directly to the tracer equation.
Instead,
the model continues to evaluate the $G_\theta$ term (first step)
as it use to do with 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,
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
$
those 2 formulations are equivalent. 

The implementation in the MITgcm follows this scheme.
The model "geometry" (here the {\bf hFacC,W,S}) is updated
just before entering {\it SOLVE\_FOR\_PRESSURE},
according to the current $\eta$ field.
Then, at the end of the time step, the variables are
advanced in time, so that $\eta^n$ becomes $\eta^{n-1}$.
At the next time step, the tracer tendency ($G_\theta$) is computed 
using the same geometry, that 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}     

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.

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}.

1	jmc	1.3	% $Header: /u/gcmpack/mitgcmdoc/part2/nonlin_frsurf.tex,v 1.2 2001/08/17 18:38:10 jmc Exp $
2	jmc	1.1	% $Name: $
3
4			%\section{Time Integration}
5
6			\subsection{Non-linear free surface}
7
8			Recently, 2 options have added to the model
9			(and therefore, have not yet been extensively tested)
10			that concern the free surface formulation.
11
12			%------------------------------------------
13			\subsubsection{Non-uniform linear-relation for the surface potential}
14
15			The linear relation between
16			surface pressure / geo- potential ($\Phi_{surf}$)
17			and surface displacement ($\eta$)
18	jmc	1.2	has been considered as uniform ($b_s =$ Constant)
19	jmc	1.3	\marginpar{add a reference to part.1 here}
20			but is in fact dependent on the position ($x,y,r$)
21	jmc	1.2	since we linearize:
22			$$\Phi_{surf}=\int_{R_o}^{R_o+\eta} b dr \simeq b_s \eta
23			~\mathrm{with}~ b_s = b(\theta,S,r)_{r=R_o}
24			\simeq b_s(\theta_{ref}(R_o),S_{ref}(R_o),R_o)$$
25			Note that, for convinience, the effect of the local
26			surface $\theta,S$ on $b_s$
27			are not considered here, but incorporated in $\Phi'_{hyd}$.
28	jmc	1.1
29	jmc	1.2	For the ocean, $b_s = g \rho_{surf} / \rho_o \simeq g$
30	jmc	1.1	is a fairly good approximation since the relative difference
31	jmc	1.2	in surface density are usually small and only due to
32			local $\theta,S$ gradient (because $R_o = 0$);
33			Therefore, they can easely be incorporated in $\Phi'_{hyd}$.
34
35	jmc	1.1	For the atmosphere, because of topographic effects,
36			the reference surface pressure $R_o$ has large spacial differences
37	jmc	1.2	that are responsible for significant $b_s$ variations
38			(from 0.8 to 1.2 $[m^3/kg]$). For this reason,
39			we use a non-uniform linear coefficient $b_s$.
40	jmc	1.1
41			Practically, in an oceanic configuration or
42	jmc	1.2	when the default value (TRUE) of the parameter
43			{\bf uniformLin\_PhiSurf} is used
44			then $b_s$ is simply set to $g$ for the ocean
45	jmc	1.1	and $1.$ for the atmosphere.\\
46			Turning {\bf uniformLin\_PhiSurf} to "FALSE", allows to
47	jmc	1.2	evaluate $b_s$ from the reference vertical profile $\theta_{ref}$
48	jmc	1.1	({\it S/R INI\_LINEAR\_PHISURF})
49			according to the reference surface pressure $P_o$ ($=R_o$):
50	jmc	1.2	$b_s = c_p \kappa (P_o / Pc)^{(\kappa - 1)} \theta_{ref}(P_o)$
51	jmc	1.1
52			%------------------------------------------
53			\subsubsection{Free surface effect on column total thickness
54			(Non-linear free surface)}
55
56			The total thickness of the fluid column is
57	jmc	1.3	$r_{surf} - R_{fixed} = \eta + R_o - R_{fixed}$
58	jmc	1.1	In the linear free surface approximation
59			(detailed before), only the fixed part of
60	jmc	1.3	it ($R_o - R_{fixed})$ is considered when we integrate the
61	jmc	1.1	continuity equation or compute tracer and momentum advection term.
62
63			This approximation is dropped when using
64	jmc	1.2	the non-linear free surface formulation.
65	jmc	1.1	Details follow here after for the barotropic part
66	jmc	1.2	and in the 2 next subsections for the baroclinic
67	jmc	1.1	part.
68
69			%------------------------------------------
70	jmc	1.2	% Non-Linear Barotropic part
71	jmc	1.1
72			The continuous form of the model equations remains
73			unchanged, except for the 2D continuity equation
74	jmc	1.2	(\ref{eq-tCsC-eta}) that is now integrated
75	jmc	1.3	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			Since $\eta$ has a direct effect on the horizontal
85	jmc	1.2	velocity (through $\nabla_h \Phi_{surf}$), this
86			adds a non-linear term to the free surface equation.
87	jmc	1.1
88	jmc	1.2	Regarding the time discretization of this non-linear part,
89			several options are now being tested:
90	jmc	1.1
91			With the column thickness evaluated at time step $n$,
92			and the surface potential gradient still implicit,
93	jmc	1.2	equation (\ref{eq-solve2D} \& \ref{eq-solve2D_rhs})
94	jmc	1.1	become:
95			\begin{eqnarray*}
96			\epsilon_{fs} {\eta}^{n+1} -
97	jmc	1.3	{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n}+R_o-R_{fixed})
98	jmc	1.2	{\bf \nabla}_h b_s {\eta}^{n+1}
99	jmc	1.1	= {\eta}^*
100			%\label{solve_2d}
101			\end{eqnarray*}
102			where
103			\begin{eqnarray*}
104			{\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
105	jmc	1.3	\Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o+\eta^n} \vec{\bf v}^* dr
106	jmc	1.1	\: + \: \epsilon_{fw} \Delta_t (P-E)^{n}
107			%\label{solve_2d_rhs}
108			\end{eqnarray*}
109			This method requires to update the solver matrix at each time step.
110
111			Alternatively, the non-linear contribution can be evaluated fully
112			explicitly:
113			\begin{eqnarray*}
114			\epsilon_{fs} {\eta}^{n+1} -
115	jmc	1.3	{\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{fixed})
116	jmc	1.2	{\bf \nabla}_h b_s {\eta}^{n+1}
117	jmc	1.1	= {\eta}^*
118	jmc	1.2	+{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n})
119			{\bf \nabla}_h b_s {\eta}^{n}
120	jmc	1.1	\end{eqnarray*}
121			This formulation allows to keep the initial solver matrix
122			since the non-linear free surface only affects the RHS.
123
124			An other option is a "linearized" formulation where the
125			total column thickness appears only in the integral term of
126	jmc	1.2	the RHS (\ref{eq-solve2D_rhs}) but not directly in
127			the equation (\ref{eq-solve2D}).
128	jmc	1.1
129			%------------------------------------------
130			\subsubsection{Free surface effect on the surface level thickness
131			(Non-linear free surface): Tracer advection}
132
133			To ensure a global tracer conservation (i.e., the total amount)
134			as well as the local one (see tracer section for more details),
135			the change in the surface level thickness must be consistent with
136			the way the continuity equation is integrated, both in
137			in the barotropic part (to find $\eta$) and baroclinic part
138	jmc	1.2	(to find $w = \dot{r}$).
139	jmc	1.1
140			To illustrate this, let's consider the shallow water model,
141			with uniform cartesian horizontal grid:
142			$$
143			\partial_t h + \nabla \cdot h \vec{\bf v} = 0
144			$$
145			where $h$ is the total thickness of the water column.
146			To conserve the tracer $\theta$ we have to discretize:
147			$$
148			\partial_t (h \theta) + \nabla \cdot ( h \theta \vec{\bf v})= 0
149			$$
150			Using the implicit (non-linear) free surface described before, we have:
151			\begin{eqnarray*}
152	jmc	1.2	h^{n+1} = h^{n} - \Delta_t \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) \\
153	jmc	1.1	\end{eqnarray*}
154			The discretized form of the tracer equation must use the same
155			"geometry" to compute the tracer fluxes, that is, the same value of
156			h, as the continuity equation did:
157			\begin{eqnarray*}
158			h^{n+1} \, \theta^{n+1} = h^n \, \theta^n
159	jmc	1.2	- \Delta_t \nabla \cdot (h^n \, \theta^n \, \vec{\bf v}^{n+1})
160	jmc	1.1	\end{eqnarray*}
161
162			In order to deal with the Adams-Bashforth time stepping,
163	jmc	1.2	we implement this scheme slightly differently, in two step:
164	jmc	1.1	the variation of the water column thickness (from
165			$h^n$ to $h^{n+1}$)
166			is not incorporated directly to the tracer equation.
167			Instead,
168	jmc	1.2	the model continues to evaluate the $G_\theta$ term (first step)
169			as it use to do with the Linear free surface formulation
170			(with the "{\it surface correction}" turned "on", see tracer section):
171	jmc	1.1	$$
172	jmc	1.2	G_\theta^n = \left(- \nabla \cdot (h^n \, \theta^n \, \vec{\bf v}^{n+1})
173			- \dot{r}_{surf}^{n+1} \theta^n \right) / h^n
174	jmc	1.1	$$
175	jmc	1.2	Then in a second step,
176	jmc	1.1	thickness variation (expansion/reduction) is taken into account :
177			$$
178	jmc	1.2	\theta^{n+1} = \theta^n + \Delta_t \frac{h^n}{h^{n+1}} G_\theta^{(n+1/2)}
179	jmc	1.1	$$
180			Note that with a simple forward time step (no Adams-Bashforth),
181			since
182			$
183	jmc	1.2	\dot{r}_{surf}^{n+1}
184			= - \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) = (h^{n+1} - h^{n})
185	jmc	1.1	/ \Delta_t
186			$
187			those 2 formulations are equivalent.
188
189			The implementation in the MITgcm follows this scheme.
190			The model "geometry" (here the {\bf hFacC,W,S}) is updated
191	jmc	1.2	just before entering {\it SOLVE\_FOR\_PRESSURE},
192	jmc	1.1	according to the current $\eta$ field.
193			Then, at the end of the time step, the variables are
194			advanced in time, so that $\eta^n$ becomes $\eta^{n-1}$.
195			At the next time step, the tracer tendency ($G_\theta$) is computed
196			using the same geometry, that is now consistent with
197			$\eta^{n-1}$.
198			Finally, in S/R {\it TIMESTEP\_TRACER}, the expansion/reduction
199			ratio is applied to the surface level to compute the new tracer field.
200
201			%------------------------------------------
202			\subsubsection{Free surface effect on the surface level thickness
203			(Non-linear free surface): Momentum advection}
204
205			Regarding momentum advection,
206			the vector invariant formulation is similar to the
207			advective form (as opposed to the flux form) and therefore
208			does not need specific modification to include the
209			free surface effect on the surface level thickness.
210			Updating the {\bf hFacC,W,S} and the {\bf recip\_hFac}(s)
211			at one given place (like describe before) is sufficient.
212
213			With the flux form formulation, advection of momentum
214	jmc	1.2	can be treated exactly as the tracer advection is.
215	jmc	1.1	Here the expansion/reduction factors ($hFacW^{n+1}/hFacW^n$ for $u$
216			and $hFacS^{n+1}/hFacS^n$ for $v$) are simply applied in the
217			subroutine {\it TIMESTEP}.
218