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	% $Header: /u/gcmpack/mitgcmdoc/part2/nonlin_frsurf.tex,v 1.2 2001/08/17 18:38:10 jmc Exp $
2	% $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	has been considered as uniform ($b_s =$ Constant)
19	\marginpar{add a reference to part.1 here}
20	but is in fact dependent on the position ($x,y,r$)
21	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
29	For the ocean, $b_s = g \rho_{surf} / \rho_o \simeq g$
30	is a fairly good approximation since the relative difference
31	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	For the atmosphere, because of topographic effects,
36	the reference surface pressure $R_o$ has large spacial differences
37	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
41	Practically, in an oceanic configuration or
42	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	and $1.$ for the atmosphere.\\
46	Turning {\bf uniformLin\_PhiSurf} to "FALSE", allows to
47	evaluate $b_s$ from the reference vertical profile $\theta_{ref}$
48	({\it S/R INI\_LINEAR\_PHISURF})
49	according to the reference surface pressure $P_o$ ($=R_o$):
50	$b_s = c_p \kappa (P_o / Pc)^{(\kappa - 1)} \theta_{ref}(P_o)$
51
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	$r_{surf} - R_{fixed} = \eta + R_o - R_{fixed}$
58	In the linear free surface approximation
59	(detailed before), only the fixed part of
60	it ($R_o - R_{fixed})$ is considered when we integrate the
61	continuity equation or compute tracer and momentum advection term.
62
63	This approximation is dropped when using
64	the non-linear free surface formulation.
65	Details follow here after for the barotropic part
66	and in the 2 next subsections for the baroclinic
67	part.
68
69	%------------------------------------------
70	% Non-Linear Barotropic part
71
72	The continuous form of the model equations remains
73	unchanged, except for the 2D continuity equation
74	(\ref{eq-tCsC-eta}) that is now integrated
75	from $R_{fixed}(x,y)$ up to $r_{surf}=R_o+\eta$ :
76
77	\begin{displaymath}
78	\epsilon_{fs} \partial_t \eta =
79	\left. \dot{r} \right\|_{r=r_{surf}} + \epsilon_{fw} (P-E) =
80	- {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o+\eta} \vec{\bf v} dr
81	+ \epsilon_{fw} (P-E)
82	\end{displaymath}
83
84	Since $\eta$ has a direct effect on the horizontal
85	velocity (through $\nabla_h \Phi_{surf}$), this
86	adds a non-linear term to the free surface equation.
87
88	Regarding the time discretization of this non-linear part,
89	several options are now being tested:
90
91	With the column thickness evaluated at time step $n$,
92	and the surface potential gradient still implicit,
93	equation (\ref{eq-solve2D} \& \ref{eq-solve2D_rhs})
94	become:
95	\begin{eqnarray*}
96	\epsilon_{fs} {\eta}^{n+1} -
97	{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n}+R_o-R_{fixed})
98	{\bf \nabla}_h b_s {\eta}^{n+1}
99	= {\eta}^*
100	%\label{solve_2d}
101	\end{eqnarray*}
102	where
103	\begin{eqnarray*}
104	{\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
105	\Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o+\eta^n} \vec{\bf v}^* dr
106	\: + \: \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	{\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{fixed})
116	{\bf \nabla}_h b_s {\eta}^{n+1}
117	= {\eta}^*
118	+{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n})
119	{\bf \nabla}_h b_s {\eta}^{n}
120	\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	the RHS (\ref{eq-solve2D_rhs}) but not directly in
127	the equation (\ref{eq-solve2D}).
128
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	(to find $w = \dot{r}$).
139
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	h^{n+1} = h^{n} - \Delta_t \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) \\
153	\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	- \Delta_t \nabla \cdot (h^n \, \theta^n \, \vec{\bf v}^{n+1})
160	\end{eqnarray*}
161
162	In order to deal with the Adams-Bashforth time stepping,
163	we implement this scheme slightly differently, in two step:
164	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	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	$$
172	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	$$
175	Then in a second step,
176	thickness variation (expansion/reduction) is taken into account :
177	$$
178	\theta^{n+1} = \theta^n + \Delta_t \frac{h^n}{h^{n+1}} G_\theta^{(n+1/2)}
179	$$
180	Note that with a simple forward time step (no Adams-Bashforth),
181	since
182	$
183	\dot{r}_{surf}^{n+1}
184	= - \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) = (h^{n+1} - h^{n})
185	/ \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	just before entering {\it SOLVE\_FOR\_PRESSURE},
192	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	can be treated exactly as the tracer advection is.
215	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