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