1 |
adcroft |
1.1 |
% $Header: $ |
2 |
|
|
% $Name: $ |
3 |
|
|
|
4 |
|
|
The convention used in this section is as follows: |
5 |
|
|
Time is "discretize" using a time step $\Delta t$ |
6 |
|
|
and $\Phi^n$ refers to the variable $\Phi$ |
7 |
|
|
at time $t = n \Delta t$ . We used the notation $\Phi^{(n)}$ |
8 |
|
|
when time interpolation is required to estimate the value of $\phi$ |
9 |
|
|
at the time $n \Delta t$. |
10 |
|
|
|
11 |
|
|
\section{Time Integration} |
12 |
|
|
|
13 |
|
|
The discretization in time of the model equations (cf section I ) |
14 |
|
|
does not depend of the discretization in space of each |
15 |
|
|
term, so that this section can be read independently. |
16 |
|
|
Also for this purpose, we will refers to the continuous |
17 |
|
|
space-derivative form of model equations, and focus on |
18 |
|
|
the time discretization. |
19 |
|
|
|
20 |
|
|
The continuous form of the model equations is: |
21 |
|
|
|
22 |
|
|
\begin{eqnarray} |
23 |
|
|
\partial_t \theta & = & G_\theta |
24 |
|
|
\\ |
25 |
|
|
\partial_t S & = & G_s |
26 |
|
|
\\ |
27 |
|
|
b' & = & b'(\theta,S,r) |
28 |
|
|
\\ |
29 |
|
|
\partial_r \phi'_{hyd} & = & -b' |
30 |
|
|
\label{eq-r-split-hyd} |
31 |
|
|
\\ |
32 |
|
|
\partial_t \vec{\bf v} |
33 |
|
|
+ {\bf \nabla}_r B_o \eta |
34 |
|
|
+ \epsilon_{nh} {\bf \nabla}_r \phi'_{nh} |
35 |
|
|
& = & \vec{\bf G}_{\vec{\bf v}} |
36 |
|
|
- {\bf \nabla}_r \phi'_{hyd} |
37 |
|
|
\label{eq-r-split-hmom} |
38 |
|
|
\\ |
39 |
|
|
\epsilon_{nh} \frac {\partial{\dot{r}}}{\partial{t}} |
40 |
|
|
+ \epsilon_{nh} \partial_r \phi'_{nh} |
41 |
|
|
& = & \epsilon_{nh} G_{\dot{r}} |
42 |
|
|
\label{eq-r-split-rdot} |
43 |
|
|
\\ |
44 |
|
|
{\bf \nabla}_r \cdot \vec{\bf v} + \partial_r \dot{r} |
45 |
|
|
& = & 0 |
46 |
|
|
\label{eq-r-cont} |
47 |
|
|
\end{eqnarray} |
48 |
|
|
where |
49 |
|
|
\begin{eqnarray*} |
50 |
|
|
G_\theta & = & |
51 |
|
|
- \vec{\bf v} \cdot {\bf \nabla}_r \theta + {\cal Q}_\theta |
52 |
|
|
\\ |
53 |
|
|
G_S & = & |
54 |
|
|
- \vec{\bf v} \cdot {\bf \nabla}_r S + {\cal Q}_S |
55 |
|
|
\\ |
56 |
|
|
\vec{\bf G}_{\vec{\bf v}} |
57 |
|
|
& = & |
58 |
|
|
- \vec{\bf v} \cdot {\bf \nabla}_r \vec{\bf v} |
59 |
|
|
- f \hat{\bf k} \wedge \vec{\bf v} |
60 |
|
|
+ \vec{\cal F}_{\vec{\bf v}} |
61 |
|
|
\\ |
62 |
|
|
G_{\dot{r}} |
63 |
|
|
& = & |
64 |
|
|
- \vec{\bf v} \cdot {\bf \nabla}_r \dot{r} |
65 |
|
|
+ {\cal F}_{\dot{r}} |
66 |
|
|
\end{eqnarray*} |
67 |
|
|
The exact form of all the "{\it G}"s terms is described in the next |
68 |
|
|
section (). Here its sufficient to mention that they contains |
69 |
|
|
all the RHS terms except the pressure / geo- potential terms. |
70 |
|
|
|
71 |
|
|
The switch $\epsilon_{nh}$ allows to activate the non hydrostatic |
72 |
|
|
mode ($\epsilon_{nh}=1$) for the ocean model. Otherwise, |
73 |
|
|
in the hydrostatic limit $\epsilon_{nh} = 0$ and the 3rd equation vanishes. |
74 |
|
|
|
75 |
|
|
The equation for $\eta$ is obtained by integrating the |
76 |
|
|
continuity equation over the entire depth of the fluid, |
77 |
|
|
from $R_{min}(x,y)$ up to $R_o(x,y)$ |
78 |
|
|
(Linear free surface): |
79 |
|
|
\begin{displaymath} |
80 |
|
|
\epsilon_{fs} \partial_t \eta = |
81 |
|
|
\left. \dot{r} \right|_{r=r_{surf}} + \epsilon_{fw} (P-E) = |
82 |
|
|
- {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o} \vec{\bf v} dr |
83 |
|
|
+ \epsilon_{fw} (P-E) |
84 |
|
|
\end{displaymath} |
85 |
|
|
|
86 |
|
|
Where $\epsilon_{fs}$,$\epsilon_{fw}$ are two flags to |
87 |
|
|
distinguish between a free-surface equation ($\epsilon_{fs}=1$) |
88 |
|
|
or the rigid-lid approximation ($\epsilon_{fs}=0$); |
89 |
|
|
and to distinguish when exchange of Fresh-Water is included |
90 |
|
|
at the ocean surface (natural BC) ($\epsilon_{fw} = 1$) |
91 |
|
|
or not ($\epsilon_{fw} = 0$). |
92 |
|
|
|
93 |
|
|
The hydrostatic potential is found by |
94 |
|
|
integrating \ref{eq-r-split-hyd} with the boundary condition that |
95 |
|
|
$\phi'_{hyd}(r=R_o) = 0$: |
96 |
|
|
\begin{eqnarray*} |
97 |
|
|
& & |
98 |
|
|
\int_{r'}^{R_o} \partial_r \phi'_{hyd} dr = |
99 |
|
|
\left[ \phi'_{hyd} \right]_{r'}^{R_o} = |
100 |
|
|
\int_{r'}^{R_o} - b' dr |
101 |
|
|
\\ |
102 |
|
|
\Rightarrow & & |
103 |
|
|
\phi'_{hyd}(x,y,r') = \int_{r'}^{R_o} b' dr |
104 |
|
|
\end{eqnarray*} |
105 |
|
|
|
106 |
|
|
\subsection{standard synchronous time stepping} |
107 |
|
|
|
108 |
|
|
In the standard formulation, the surface pressure is |
109 |
|
|
evaluated at time step n+1 (implicit method). |
110 |
|
|
The above set of equations is then discretized in time |
111 |
|
|
as follows: |
112 |
|
|
\begin{eqnarray} |
113 |
|
|
\left[ 1 - \partial_r \kappa_v^\theta \partial_r \right] |
114 |
|
|
\theta^{n+1} & = & \theta^* |
115 |
|
|
\\ |
116 |
|
|
\left[ 1 - \partial_r \kappa_v^S \partial_r \right] |
117 |
|
|
S^{n+1} & = & S^* |
118 |
|
|
\\ |
119 |
|
|
%{b'}^{n} & = & b'(\theta^{n},S^{n},r) |
120 |
|
|
%\partial_r {\phi'_{hyd}}^{n} & = & {-b'}^{n} |
121 |
|
|
%\\ |
122 |
|
|
{\phi'_{hyd}}^{n} & = & \int_{r'}^{R_o} b'(\theta^{n},S^{n},r) dr |
123 |
|
|
\label{eq-rtd-hyd} |
124 |
|
|
\\ |
125 |
|
|
\vec{\bf v} ^{n+1} |
126 |
|
|
+ \Delta t {\bf \nabla}_r B_o {\eta}^{n+1} |
127 |
|
|
+ \epsilon_{nh} \Delta t {\bf \nabla}_r {\phi'_{nh}}^{n+1} |
128 |
|
|
- \partial_r A_v \partial_r \vec{\bf v}^{n+1} |
129 |
|
|
& = & |
130 |
|
|
\vec{\bf v}^* |
131 |
|
|
\label{eq-rtd-hmom} |
132 |
|
|
\\ |
133 |
|
|
\epsilon_{fs} {\eta}^{n+1} + \Delta t |
134 |
|
|
{\bf \nabla}_r \cdot \int_{R_{min}}^{R_o} \vec{\bf v}^{n+1} dr |
135 |
|
|
& = & |
136 |
|
|
\epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta_t (P-E)^{n} |
137 |
|
|
\nonumber |
138 |
|
|
\\ |
139 |
|
|
% = \epsilon_{fs} {\eta}^{n} & + & \epsilon_{fw} \Delta_t (P-E)^{n} |
140 |
|
|
\label{eq-rtd-eta} |
141 |
|
|
\\ |
142 |
|
|
\epsilon_{nh} \left( \dot{r} ^{n+1} |
143 |
|
|
+ \Delta t \partial_r {\phi'_{nh}} ^{n+1} |
144 |
|
|
\right) |
145 |
|
|
& = & \epsilon_{nh} \dot{r}^* |
146 |
|
|
\label{eq-rtd-rdot} |
147 |
|
|
\\ |
148 |
|
|
{\bf \nabla}_r \cdot \vec{\bf v}^{n+1} + \partial_r \dot{r}^{n+1} |
149 |
|
|
& = & 0 |
150 |
|
|
\label{eq-rtd-cont} |
151 |
|
|
\end{eqnarray} |
152 |
|
|
where |
153 |
|
|
\begin{eqnarray} |
154 |
|
|
\theta^* & = & |
155 |
|
|
\theta ^{n} + \Delta t G_{\theta} ^{(n+1/2)} |
156 |
|
|
\\ |
157 |
|
|
S^* & = & |
158 |
|
|
S ^{n} + \Delta t G_{S} ^{(n+1/2)} |
159 |
|
|
\\ |
160 |
|
|
\vec{\bf v}^* & = & |
161 |
|
|
\vec{\bf v}^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)} |
162 |
|
|
+ \Delta t {\bf \nabla}_r {\phi'_{hyd}}^{(n+1/2)} |
163 |
|
|
\\ |
164 |
|
|
\dot{r}^* & = & |
165 |
|
|
\dot{r} ^{n} + \Delta t G_{\dot{r}} ^{(n+1/2)} |
166 |
|
|
\end{eqnarray} |
167 |
|
|
|
168 |
|
|
Note that implicit vertical terms (viscosity and diffusivity) are |
169 |
|
|
not considered as part of the "{\it G}" terms, but are |
170 |
|
|
written separately here. |
171 |
|
|
|
172 |
|
|
To ensure a second order time discretization for both |
173 |
|
|
momentum and tracer, |
174 |
|
|
The "G" terms are "extrapolated" forward in time |
175 |
|
|
(Adams Bashforth time stepping) |
176 |
|
|
from the values computed at time step $n$ and $n-1$ |
177 |
|
|
to the time $n+1/2$: |
178 |
|
|
$$G^{(n+1/2)} = G^n + (1/2+\epsilon_{AB}) (G^n - G^{n-1})$$ |
179 |
|
|
A small number for the parameter $\epsilon_{AB}$ is generally used |
180 |
|
|
to stabilize this time stepping. |
181 |
|
|
|
182 |
|
|
In the standard non-stagger formulation, |
183 |
|
|
the Adams-Bashforth time stepping is also applied to the |
184 |
|
|
hydrostatic (pressure / geo-) potential term $\nabla_h \Phi'_{hyd}$. |
185 |
|
|
Note that presently, this term is in fact incorporated to the |
186 |
|
|
$\vec{\bf G}_{\vec{\bf v}}$ arrays ({\it gU,gV}). |
187 |
|
|
|
188 |
|
|
\subsection{general method} |
189 |
|
|
|
190 |
|
|
The general algorithm consist in a "predictor step" that computes |
191 |
|
|
the forward tendencies ("G" terms") and all |
192 |
|
|
the "first guess" values star notation): |
193 |
|
|
$\vec{\bf v}^*, \theta^*, S^*$ (and $\dot{r}^*$ |
194 |
|
|
in non-hydrostatic mode). This is done in the routine {\it DYNAMICS}. |
195 |
|
|
|
196 |
|
|
Then the implicit terms that appear here on the left hand side (LHS), |
197 |
|
|
are solved as follows: |
198 |
|
|
Since implicit vertical diffusion and viscosity terms |
199 |
|
|
are independent from the barotropic flow adjustment, |
200 |
|
|
they are computed first, solving a 3 diagonal Nr x Nr linear system, |
201 |
|
|
and incorporated at the end of the {\it DYNAMICS} routines. |
202 |
|
|
Then the surface pressure and non hydrostatic pressure |
203 |
|
|
are evaluated ({\it SOLVE\_FOR\_PRESSURE}); |
204 |
|
|
|
205 |
|
|
Finally, within a "corrector step', |
206 |
|
|
(routine {\it THE\_CORRECTION\_STEP}) |
207 |
|
|
the new values of $u,v,w,\theta,S$ |
208 |
|
|
are derived according to the above equations |
209 |
|
|
(see details in II.1.3). |
210 |
|
|
|
211 |
|
|
At this point, the regular time step is over, but |
212 |
|
|
the correction step contains also other optional |
213 |
|
|
adjustments such as convective adjustment algorithm, or filters |
214 |
|
|
(zonal FFT filter, shapiro filter) |
215 |
|
|
that applied on both momentum and tracer fields. |
216 |
|
|
just before setting the $n+1$ new time step value. |
217 |
|
|
|
218 |
|
|
Since the pressure solver precision is of the order of |
219 |
|
|
the "target residual" that could be lower than the |
220 |
|
|
the computer truncation error, and also because some filters |
221 |
|
|
might alter the divergence part of the flow field, |
222 |
|
|
a final evaluation of the surface r anomaly $\eta^{n+1}$ |
223 |
|
|
is performed, according to \ref{eq-rtd-eta} ({\it CALC\_EXACT\_ETA}). |
224 |
|
|
This ensures a perfect volume conservation. |
225 |
|
|
Note that there is no need for an equivalent Non-hydrostatic |
226 |
|
|
"exact conservation" step, since W is already computed after |
227 |
|
|
the filters applied. |
228 |
|
|
|
229 |
|
|
optionnal forcing terms (package):\\ |
230 |
|
|
Regarding optional forcing terms (usually part of a "package"), |
231 |
|
|
that a account for a specific source or sink term (e.g.: condensation |
232 |
|
|
as a sink of water vapor Q), they are generally incorporated |
233 |
|
|
in the main algorithm as follows; |
234 |
|
|
At the the beginning of the time step, |
235 |
|
|
the additionnal tendencies are computed |
236 |
|
|
as function of the present state (time step $n$) and external forcing ; |
237 |
|
|
Then within the main part of model, |
238 |
|
|
only those new tendencies are added to the model variables. |
239 |
|
|
|
240 |
|
|
[mode details needed] |
241 |
|
|
The atmospheric physics follows this general scheme. |
242 |
|
|
|
243 |
|
|
\subsection{stagger baroclinic time stepping} |
244 |
|
|
|
245 |
|
|
An alternative is to evaluate $\phi'_{hyd}$ with the |
246 |
|
|
new tracer fields, and step forward the momentum after. |
247 |
|
|
This option is known as stagger baroclinic time stepping, |
248 |
|
|
since tracer and momentum are step forward in time one after the other. |
249 |
|
|
It can be activated turning on a running flag parameter |
250 |
|
|
{\it staggerTimeStep} in file "{\it data}"). |
251 |
|
|
|
252 |
|
|
The main advantage of this time stepping compared to a synchronous one, |
253 |
|
|
is a better stability, specially regarding internal gravity waves, |
254 |
|
|
and a very natural implementation of a 2nd order in time |
255 |
|
|
hydrostatic pressure / geo- potential term. |
256 |
|
|
In the other hand, a synchronous time step might be better |
257 |
|
|
for convection problems; Its also make simpler time dependent forcing |
258 |
|
|
and diagnostic implementation ; and allows a more efficient threading. |
259 |
|
|
|
260 |
|
|
Although the stagger time step does not affect deeply the |
261 |
|
|
structure of the code --- a switch allows to evaluate the |
262 |
|
|
hydrostatic pressure / geo- potential from new $\theta,S$ |
263 |
|
|
instead of the Adams-Bashforth estimation --- |
264 |
|
|
this affect the way the time discretization is presented : |
265 |
|
|
|
266 |
|
|
\begin{eqnarray*} |
267 |
|
|
\left[ 1 - \partial_r \kappa_v^\theta \partial_r \right] |
268 |
|
|
\theta^{n+1/2} & = & \theta^* |
269 |
|
|
\\ |
270 |
|
|
\left[ 1 - \partial_r \kappa_v^S \partial_r \right] |
271 |
|
|
S^{n+1/2} & = & S^* |
272 |
|
|
\end{eqnarray*} |
273 |
|
|
with |
274 |
|
|
\begin{eqnarray*} |
275 |
|
|
\theta^* & = & |
276 |
|
|
\theta ^{(n-1/2)} + \Delta t G_{\theta} ^{(n)} |
277 |
|
|
\\ |
278 |
|
|
S^* & = & |
279 |
|
|
S ^{(n-1/2)} + \Delta t G_{S} ^{(n)} |
280 |
|
|
\end{eqnarray*} |
281 |
|
|
And |
282 |
|
|
\begin{eqnarray*} |
283 |
|
|
%{b'}^{n+1/2} & = & b'(\theta^{n+1/2},S^{n+1/2},r) |
284 |
|
|
%\\ |
285 |
|
|
%\partial_r {\phi'_{hyd}}^{n+1/2} & = & {-b'}^{n+1/2} |
286 |
|
|
{\phi'_{hyd}}^{n+1/2} & = & \int_{r'}^{R_o} b'(\theta^{n+1/2},S^{n+1/2},r) dr |
287 |
|
|
%\label{eq-rtd-hyd} |
288 |
|
|
\\ |
289 |
|
|
\vec{\bf v} ^{n+1} |
290 |
|
|
+ \Delta t {\bf \nabla}_r B_o {\eta}^{n+1} |
291 |
|
|
+ \epsilon_{nh} \Delta t {\bf \nabla}_r {\phi'_{nh}}^{n+1} |
292 |
|
|
- \partial_r A_v \partial_r \vec{\bf v}^{n+1} |
293 |
|
|
& = & |
294 |
|
|
\vec{\bf v}^* |
295 |
|
|
%\label{eq-rtd-hmom} |
296 |
|
|
\\ |
297 |
|
|
\epsilon_{fs} {\eta}^{n+1} + \Delta t |
298 |
|
|
{\bf \nabla}_r \cdot \int_{R_{min}}^{R_o} \vec{\bf v}^{n+1} dr |
299 |
|
|
& = & |
300 |
|
|
\epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta_t (P-E)^{n} |
301 |
|
|
\\ |
302 |
|
|
\epsilon_{nh} \left( \dot{r} ^{n+1} |
303 |
|
|
+ \Delta t \partial_r {\phi'_{nh}} ^{n+1} |
304 |
|
|
\right) |
305 |
|
|
& = & \epsilon_{nh} \dot{r}^* |
306 |
|
|
%\label{eq-rtd-rdot} |
307 |
|
|
\\ |
308 |
|
|
{\bf \nabla}_r \cdot \vec{\bf v}^{n+1} + \partial_r \dot{r}^{n+1} |
309 |
|
|
& = & 0 |
310 |
|
|
%\label{eq-rtd-cont} |
311 |
|
|
\end{eqnarray*} |
312 |
|
|
with |
313 |
|
|
\begin{eqnarray*} |
314 |
|
|
\vec{\bf v}^* & = & |
315 |
|
|
\vec{\bf v} ^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)} |
316 |
|
|
+ \Delta t {\bf \nabla}_r {\phi'_{hyd}}^{n+1/2} |
317 |
|
|
\\ |
318 |
|
|
\dot{r}^* & = & |
319 |
|
|
\dot{r} ^{n} + \Delta t G_{\dot{r}} ^{(n+1/2)} |
320 |
|
|
\end{eqnarray*} |
321 |
|
|
|
322 |
|
|
%--------------------------------------------------------------------- |
323 |
|
|
|
324 |
|
|
\subsection{surface pressure} |
325 |
|
|
|
326 |
|
|
Substituting \ref{eq-rtd-hmom} into \ref{eq-rtd-cont}, assuming |
327 |
|
|
$\epsilon_{nh} = 0$ yields a Helmholtz equation for ${\eta}^{n+1}$: |
328 |
|
|
$$ |
329 |
|
|
\epsilon_{fs} {\eta}^{n+1} - |
330 |
|
|
{\bf \nabla}_r \cdot \Delta t^2 (R_o-R_{min}) |
331 |
|
|
{\bf \nabla}_r B_o {\eta}^{n+1} |
332 |
|
|
= {\eta}^* |
333 |
|
|
\label{solve_2d} |
334 |
|
|
$$ |
335 |
|
|
where |
336 |
|
|
$$ |
337 |
|
|
{\eta}^* = \epsilon_{fs} \: {\eta}^{n} - |
338 |
|
|
\Delta t {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o} \vec{\bf v}^* dr |
339 |
|
|
\: + \: \epsilon_{fw} \Delta_t (P-E)^{n} |
340 |
|
|
$$ |
341 |
|
|
|
342 |
|
|
Once ${\eta}^{n+1}$ has been found substituting into \ref{eq-rtd-hmom} |
343 |
|
|
would yield $\vec{\bf v}^{n+1}$ if the model is hydrostatic |
344 |
|
|
($\epsilon_{nh}=0$): |
345 |
|
|
$$ |
346 |
|
|
\vec{\bf v}^{n+1} = \vec{\bf v}^{*} |
347 |
|
|
- \Delta t {\bf \nabla}_r B_o {\eta}^{n+1} |
348 |
|
|
$$ |
349 |
|
|
|
350 |
|
|
This is known as the correction step. However, when the model is |
351 |
|
|
non-hydrostatic ($\epsilon_{nh}=1$) we need an additional step and an |
352 |
|
|
additional equation for $\phi'_{nh}$. This is obtained by |
353 |
|
|
substituting \ref{eq-rtd-hmom} and \ref{eq-rtd-rdot} into |
354 |
|
|
\ref{eq-rtd-cont}: |
355 |
|
|
\begin{equation} |
356 |
|
|
\left[ {\bf \nabla}_r^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1} |
357 |
|
|
= \frac{1}{\Delta t} \left( |
358 |
|
|
{\bf \nabla}_r \cdot \vec{\bf v}^{**} + \partial_r \dot{r}^* \right) |
359 |
|
|
\end{equation} |
360 |
|
|
where |
361 |
|
|
\begin{displaymath} |
362 |
|
|
\vec{\bf v}^{**} = \vec{\bf v}^* - \Delta t {\bf \nabla}_r B_o {\eta}^{n+1} |
363 |
|
|
\end{displaymath} |
364 |
|
|
Note that $\eta^{n+1}$ is also used to update the second RHS term |
365 |
|
|
$\partial_r \dot{r}^* $ since |
366 |
|
|
the vertical velocity at the surface ($\dot{r}_{surf}$) |
367 |
|
|
is evaluated as $(\eta^{n+1} - \eta^n) / \Delta t$. |
368 |
|
|
|
369 |
|
|
Finally, the horizontal velocities at the new time level are found by: |
370 |
|
|
\begin{equation} |
371 |
|
|
\vec{\bf v}^{n+1} = \vec{\bf v}^{**} |
372 |
|
|
- \epsilon_{nh} \Delta t {\bf \nabla}_r {\phi'_{nh}}^{n+1} |
373 |
|
|
\end{equation} |
374 |
|
|
and the vertical velocity is found by integrating the continuity |
375 |
|
|
equation vertically. |
376 |
|
|
Note that, for convenience regarding the restart procedure, |
377 |
|
|
the integration of the continuity equation has been |
378 |
|
|
moved at the beginning of the time step (instead of at the end), |
379 |
|
|
without any consequence on the solution. |
380 |
|
|
|
381 |
|
|
Regarding the implementation, all those computation are done |
382 |
|
|
within the routine {\it SOLVE\_FOR\_PRESSURE} and its dependent |
383 |
|
|
{\it CALL}s. |
384 |
|
|
The standard method to solve the 2D elliptic problem (\ref{solve_2d}) |
385 |
|
|
uses the conjugate gradient method (routine {\it CG2D}); The |
386 |
|
|
the solver matrix and conjugate gradient operator are only function |
387 |
|
|
of the discretized domain and are therefore evaluated separately, |
388 |
|
|
before the time iteration loop, within {\it INI\_CG2D}. |
389 |
|
|
The computation of the RHS $\eta^*$ is partly |
390 |
|
|
done in {\it CALC\_DIV\_GHAT} and in {\it SOLVE\_FOR\_PRESSURE}. |
391 |
|
|
|
392 |
|
|
The same method is applied for the non hydrostatic part, using |
393 |
|
|
a conjugate gradient 3D solver ({\it CG3D}) that is initialized |
394 |
|
|
in {\it INI\_CG3D}. The RHS terms of 2D and 3D problems |
395 |
|
|
are computed together, within the same part of the code. |
396 |
|
|
|
397 |
|
|
\newpage |
398 |
|
|
%----------------------------------------------------------------------------------- |
399 |
|
|
\subsection{Crank-Nickelson barotropic time stepping} |
400 |
|
|
|
401 |
|
|
The full implicit time stepping described previously is unconditionally stable |
402 |
|
|
but damps the fast gravity waves, resulting in a loss of |
403 |
|
|
gravity potential energy. |
404 |
|
|
The modification presented hereafter allows to combine an implicit part |
405 |
|
|
($\beta,\gamma$) and an explicit part ($1-\beta,1-\gamma$) for the surface |
406 |
|
|
pressure gradient ($\beta$) and for the barotropic flow divergence ($\gamma$). |
407 |
|
|
\\ |
408 |
|
|
For instance, $\beta=\gamma=1$ is the previous fully implicit scheme; |
409 |
|
|
$\beta=\gamma=1/2$ is the non damping (energy conserving), unconditionally |
410 |
|
|
stable, Crank-Nickelson scheme; $(\beta,\gamma)=(1,0)$ or $=(0,1)$ |
411 |
|
|
corresponds to the forward - backward scheme that conserves energy but is |
412 |
|
|
only stable for small time steps.\\ |
413 |
|
|
In the code, $\beta,\gamma$ are defined as parameters, respectively |
414 |
|
|
{\it implicSurfPress}, {\it implicDiv2DFlow}. They are read from |
415 |
|
|
the main data file "{\it data}" and are set by default to 1,1. |
416 |
|
|
|
417 |
|
|
Equations \ref{eq-rtd-hmom} and \ref{eq-rtd-eta} are modified as follows: |
418 |
|
|
$$ |
419 |
|
|
\frac{ \vec{\bf v}^{n+1} }{ \Delta t } |
420 |
|
|
+ {\bf \nabla}_r B_o [ \beta {\eta}^{n+1} + (1-\beta) {\eta}^{n} ] |
421 |
|
|
+ \epsilon_{nh} {\bf \nabla}_r {\phi'_{nh}}^{n+1} |
422 |
|
|
= \frac{ \vec{\bf v}^* }{ \Delta t } |
423 |
|
|
$$ |
424 |
|
|
$$ |
425 |
|
|
\epsilon_{fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t} |
426 |
|
|
+ {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o} |
427 |
|
|
[ \gamma \vec{\bf v}^{n+1} + (1-\gamma) \vec{\bf v}^{n}] dr |
428 |
|
|
= \epsilon_{fw} (P-E) |
429 |
|
|
$$ |
430 |
|
|
where: |
431 |
|
|
\begin{eqnarray*} |
432 |
|
|
\vec{\bf v}^* & = & |
433 |
|
|
\vec{\bf v} ^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)} |
434 |
|
|
+ (\beta-1) \Delta t {\bf \nabla}_r B_o {\eta}^{n} |
435 |
|
|
+ \Delta t {\bf \nabla}_r {\phi'_{hyd}}^{(n+1/2)} |
436 |
|
|
\\ |
437 |
|
|
{\eta}^* & = & |
438 |
|
|
\epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta t (P-E) |
439 |
|
|
- \Delta t {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o} |
440 |
|
|
[ \gamma \vec{\bf v}^* + (1-\gamma) \vec{\bf v}^{n}] dr |
441 |
|
|
\end{eqnarray*} |
442 |
|
|
\\ |
443 |
|
|
In the hydrostatic case ($\epsilon_{nh}=0$), |
444 |
|
|
this allow to find ${\eta}^{n+1}$, according to: |
445 |
|
|
$$ |
446 |
|
|
\epsilon_{fs} {\eta}^{n+1} - |
447 |
|
|
{\bf \nabla}_r \cdot \beta\gamma \Delta t^2 B_o (R_o - R_{min}) |
448 |
|
|
{\bf \nabla}_r {\eta}^{n+1} |
449 |
|
|
= {\eta}^* |
450 |
|
|
$$ |
451 |
|
|
and then to compute (correction step): |
452 |
|
|
$$ |
453 |
|
|
\vec{\bf v}^{n+1} = \vec{\bf v}^{*} |
454 |
|
|
- \beta \Delta t {\bf \nabla}_r B_o {\eta}^{n+1} |
455 |
|
|
$$ |
456 |
|
|
|
457 |
|
|
The non-hydrostatic part is solved as described previously. |
458 |
|
|
\\ \\ |
459 |
|
|
N.B: |
460 |
|
|
\\ |
461 |
|
|
a) The non-hydrostatic part of the code has not yet been |
462 |
|
|
updated, %since it falls out of the purpose of this test, |
463 |
|
|
so that this option cannot be used with $(\beta,\gamma) \neq (1,1)$. |
464 |
|
|
\\ |
465 |
|
|
b) To remind, the stability criteria with the Crank-Nickelson time stepping |
466 |
|
|
for the pure linear gravity wave problem in cartesian coordinate is: |
467 |
|
|
\\ |
468 |
|
|
$\star$~ $\beta + \gamma < 1$ : unstable |
469 |
|
|
\\ |
470 |
|
|
$\star$~ $\beta \geq 1/2$ and $ \gamma \geq 1/2$ : stable |
471 |
|
|
\\ |
472 |
|
|
$\star$~ $\beta + \gamma \geq 1$ : stable if |
473 |
|
|
%, for all wave length $(k\Delta x,l\Delta y)$ |
474 |
|
|
$$ |
475 |
|
|
c_{max}^2 (\beta - 1/2)(\gamma - 1/2) + 1 \geq 0 |
476 |
|
|
$$ |
477 |
|
|
$$ |
478 |
|
|
\mbox{with }~ |
479 |
|
|
%c^2 = 2 g H {\Delta t}^2 |
480 |
|
|
%(\frac{1-cos 2 \pi / k}{\Delta x^2} |
481 |
|
|
%+\frac{1-cos 2 \pi / l}{\Delta y^2}) |
482 |
|
|
%$$ |
483 |
|
|
%Practically, the most stringent condition is obtained with $k=l=2$ : |
484 |
|
|
%$$ |
485 |
|
|
c_{max} = 2 \Delta t \: \sqrt{g H} \: |
486 |
|
|
\sqrt{ \frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} } |
487 |
|
|
$$ |