13 |
part of the tracer equations and the various advection schemes will be |
part of the tracer equations and the various advection schemes will be |
14 |
described here. |
described here. |
15 |
|
|
16 |
|
\subsection{Time-stepping of tracers: ABII} |
17 |
|
|
18 |
|
The default advection scheme is the centered second order method which |
19 |
|
requires a second order or quasi-second order time-stepping scheme to |
20 |
|
be stable. Historically this has been the quasi-second order |
21 |
|
Adams-Bashforth method (ABII) and applied to all terms. For an |
22 |
|
arbitrary tracer, $\tau$, the forced advection-diffusion equation |
23 |
|
reads: |
24 |
|
\begin{equation} |
25 |
|
\partial_t \tau + G_{adv}^\tau = G_{diff}^\tau + G_{forc}^\tau |
26 |
|
\end{equation} |
27 |
|
where $G_{adv}^\tau$, $G_{diff}^\tau$ and $G_{forc}^\tau$ are the |
28 |
|
tendencies due to advection, diffusion and forcing, respectively, |
29 |
|
namely: |
30 |
|
\begin{eqnarray} |
31 |
|
G_{adv}^\tau & = & \partial_x u \tau + \partial_y v \tau + \partial_r w \tau |
32 |
|
- \tau \nabla \cdot {\bf v} \\ |
33 |
|
G_{diff}^\tau & = & \nabla \cdot {\bf K} \nabla \tau |
34 |
|
\end{eqnarray} |
35 |
|
and the forcing can be some arbitrary function of state, time and |
36 |
|
space. |
37 |
|
|
38 |
|
The term, $\tau \nabla \cdot {\bf v}$, is required to retain local |
39 |
|
conservation in conjunction with the linear implicit free-surface. It |
40 |
|
only affects the surface layer since the flow is non-divergent |
41 |
|
everywhere else. This term is therefore referred to as the surface |
42 |
|
correction term. Global conservation is not possible using the |
43 |
|
flux-form (as here) and a linearized free-surface |
44 |
|
(\cite{Griffies00,Campin02}). |
45 |
|
|
46 |
|
The continuity equation can be recovered by setting |
47 |
|
$G_{diff}=G_{forc}=0$ and $\tau=1$. |
48 |
|
|
49 |
|
The driver routine that calls the routines to calculate tendancies are |
50 |
|
{\em S/R CALC\_GT} and {\em S/R CALC\_GS} for temperature and salt |
51 |
|
(moisture), respectively. These in turn call a generic advection |
52 |
|
diffusion routine {\em S/R GAD\_CALC\_RHS} that is called with the |
53 |
|
flow field and relevent tracer as arguments and returns the collective |
54 |
|
tendancy due to advection and diffusion. Forcing is add subsequently |
55 |
|
in {\em S/R CALC\_GT} or {\em S/R CALC\_GS} to the same tendancy |
56 |
|
array. |
57 |
|
|
58 |
|
\fbox{ \begin{minipage}{4.75in} |
59 |
|
{\em S/R GAD\_CALC\_RHS} ({\em pkg/generic\_advdiff/gad\_calc\_rhs.F}) |
60 |
|
|
61 |
|
$\tau$: {\bf tracer} (argument) |
62 |
|
|
63 |
|
$G^{(n)}$: {\bf gTracer} (argument) |
64 |
|
|
65 |
|
$F_r$: {\bf fVerT} (argument) |
66 |
|
|
67 |
|
\end{minipage} } |
68 |
|
|
69 |
|
|
70 |
|
The space and time discretizations are treated seperately (method of |
71 |
|
lines). The Adams-Bashforth time discretization reads: |
72 |
|
\marginpar{$\epsilon$: {\bf AB\_eps}} |
73 |
|
\marginpar{$\Delta t$: {\bf deltaTtracer}} |
74 |
|
\begin{equation} |
75 |
|
\tau^{(n+1)} = \tau^{(n)} + \Delta t \left( |
76 |
|
(\frac{3}{2} + \epsilon) G^{(n)} - (\frac{1}{2} + \epsilon) G^{(n-1)} |
77 |
|
\right) |
78 |
|
\end{equation} |
79 |
|
where $G^{(n)} = G_{adv}^\tau + G_{diff}^\tau + G_{src}^\tau$ at time |
80 |
|
step $n$. |
81 |
|
|
82 |
|
Strictly speaking the ABII scheme should be applied only to the |
83 |
|
advection terms. However, this scheme is only used in conjuction with |
84 |
|
the standard second, third and fourth order advection |
85 |
|
schemes. Selection of any other advection scheme disables |
86 |
|
Adams-Bashforth for tracers so that explicit diffusion and forcing use |
87 |
|
the forward method. |
88 |
|
|
89 |
|
\fbox{ \begin{minipage}{4.75in} |
90 |
|
{\em S/R TIMESTEP\_TRACER} ({\em model/src/timestep\_tracer.F}) |
91 |
|
|
92 |
|
$\tau$: {\bf tracer} (argument) |
93 |
|
|
94 |
|
$G^{(n)}$: {\bf gTracer} (argument) |
95 |
|
|
96 |
|
$G^{(n-1)}$: {\bf gTrNm1} (argument) |
97 |
|
|
98 |
|
$\Delta t$: {\bf deltaTtracer} (PARAMS.h) |
99 |
|
|
100 |
|
\end{minipage} } |
101 |
|
|
102 |
|
\begin{figure} |
103 |
|
\resizebox{5.5in}{!}{\includegraphics{part2/advect-1d-lo.eps}} |
104 |
|
\caption{ |
105 |
|
Comparison of 1-D advection schemes. Courant number is 0.05 with 60 |
106 |
|
points and solutions are shown for T=1 (one complete period). |
107 |
|
a) Shows the upwind biased schemes; first order upwind, DST3, |
108 |
|
third order upwind and second order upwind. |
109 |
|
b) Shows the centered schemes; Lax-Wendroff, DST4, centered second order, |
110 |
|
centered fourth order and finite volume fourth order. |
111 |
|
c) Shows the second order flux limiters: minmod, Superbee, |
112 |
|
MC limiter and the van Leer limiter. |
113 |
|
d) Shows the DST3 method with flux limiters due to Sweby with |
114 |
|
$\mu=1$, $\mu=c/(1-c)$ and a fourth order DST method with Sweby limiter, |
115 |
|
$\mu=c/(1-c)$. |
116 |
|
\label{fig:advect-1d-lo} |
117 |
|
} |
118 |
|
\end{figure} |
119 |
|
|
120 |
|
\begin{figure} |
121 |
|
\resizebox{5.5in}{!}{\includegraphics{part2/advect-1d-hi.eps}} |
122 |
|
\caption{ |
123 |
|
Comparison of 1-D advection schemes. Courant number is 0.89 with 60 |
124 |
|
points and solutions are shown for T=1 (one complete period). |
125 |
|
a) Shows the upwind biased schemes; first order upwind and DST3. |
126 |
|
Third order upwind and second order upwind are unstable at this Courant number. |
127 |
|
b) Shows the centered schemes; Lax-Wendroff, DST4. Centered second order, |
128 |
|
centered fourth order and finite volume fourth order and unstable at this |
129 |
|
Courant number. |
130 |
|
c) Shows the second order flux limiters: minmod, Superbee, |
131 |
|
MC limiter and the van Leer limiter. |
132 |
|
d) Shows the DST3 method with flux limiters due to Sweby with |
133 |
|
$\mu=1$, $\mu=c/(1-c)$ and a fourth order DST method with Sweby limiter, |
134 |
|
$\mu=c/(1-c)$. |
135 |
|
\label{fig:advect-1d-hi} |
136 |
|
} |
137 |
|
\end{figure} |
138 |
|
|
139 |
|
\section{Linear advection schemes} |
140 |
|
|
141 |
|
The advection schemes known as centered second order, centered fourth |
142 |
|
order, first order upwind and upwind biased third order are known as |
143 |
|
linear advection schemes because the coefficient for interpolation of |
144 |
|
the advected tracer are linear and a function only of the flow, not |
145 |
|
the tracer field it self. We discuss these first since they are most |
146 |
|
commonly used in the field and most familiar. |
147 |
|
|
148 |
\subsection{Centered second order advection-diffusion} |
\subsection{Centered second order advection-diffusion} |
149 |
|
|
150 |
The basic discretization, centered second order, is the default. It is |
The basic discretization, centered second order, is the default. It is |
151 |
designed to be consistant with the continuity equation to facilitate |
designed to be consistant with the continuity equation to facilitate |
152 |
conservation properties analogous to the continuum: |
conservation properties analogous to the continuum. However, centered |
153 |
|
second order advection is notoriously noisey and must be used in |
154 |
|
conjuction with some finite amount of diffusion to produce a sensible |
155 |
|
solution. |
156 |
|
|
157 |
|
The advection operator is discretized: |
158 |
\begin{equation} |
\begin{equation} |
159 |
{\cal A}_c \Delta r_f h_c \partial_\theta |
{\cal A}_c \Delta r_f h_c G_{adv}^\tau = |
160 |
+ \delta_i F_x |
\delta_i F_x + \delta_j F_y + \delta_k F_r |
|
+ \delta_j F_y |
|
|
+ \delta_k F_r |
|
|
= {\cal A}_c \Delta r_f h_c {\cal S}_\theta |
|
|
+ \theta {\cal A}_c \delta_k (P-E)_{r=0} |
|
161 |
\end{equation} |
\end{equation} |
162 |
where the area integrated fluxes are given by: |
where the area integrated fluxes are given by: |
163 |
\begin{eqnarray} |
\begin{eqnarray} |
164 |
F_x & = & U \overline{ \theta }^i |
F_x & = & U \overline{ \tau }^i \\ |
165 |
- \kappa_h \frac{\Delta y_g \Delta r_f h_w}{\Delta x_c} \delta_i \theta \\ |
F_y & = & V \overline{ \tau }^j \\ |
166 |
F_y & = & V \overline{ \theta }^j |
F_r & = & W \overline{ \tau }^k |
|
- \kappa_h \frac{\Delta x_g \Delta r_f h_s}{\Delta y_c} \delta_j \theta \\ |
|
|
F_r & = & W \overline{ \theta }^k |
|
|
- \kappa_v \frac{\Delta x_g \Delta y_g}{\Delta r_c} \delta_k \theta |
|
167 |
\end{eqnarray} |
\end{eqnarray} |
168 |
The quantities $U$, $V$ and $W$ are volume fluxes defined: |
The quantities $U$, $V$ and $W$ are volume fluxes defined: |
169 |
\marginpar{$U$: {\bf uTrans} } |
\marginpar{$U$: {\bf uTrans} } |
174 |
V & = & \Delta x_g \Delta r_f h_s v \\ |
V & = & \Delta x_g \Delta r_f h_s v \\ |
175 |
W & = & {\cal A}_c w |
W & = & {\cal A}_c w |
176 |
\end{eqnarray} |
\end{eqnarray} |
|
${\cal S}$ represents the ``parameterized'' SGS processes and physics |
|
|
and sources associated with the tracer. For instance, potential |
|
|
temperature equation in the ocean has is forced by surface and |
|
|
partially penetrating heat fluxes: |
|
|
\begin{equation} |
|
|
{\cal A}_c \Delta r_f h_c {\cal S}_\theta = \frac{1}{c_p \rho_o} \delta_k {\cal A}_c {\cal Q} |
|
|
\end{equation} |
|
|
while the salt equation has no real sources, ${\cal S}=0$, which |
|
|
leaves just the $P-E$ term. |
|
|
|
|
|
The continuity equation can be recovered by setting ${\cal Q}=0$, $\kappa_h = \kappa_v = 0$ and |
|
|
$\theta=1$. The term $\theta (P-E)_{r=0}$ is required to retain local |
|
|
conservation of $\theta$. Global conservation is not possible using |
|
|
the flux-form (as here) and a linearized free-surface |
|
|
(\cite{Griffies00,Campin02}). |
|
177 |
|
|
178 |
|
For non-divergent flow, this discretization can be shown to conserve |
179 |
|
the tracer both locally and globally and to globally conserve tracer |
180 |
|
variance, $\tau^2$. The proof is given in \cite{Adcroft95,Adcroft97}. |
181 |
|
|
182 |
|
\fbox{ \begin{minipage}{4.75in} |
183 |
|
{\em S/R GAD\_C2\_ADV\_X} ({\em gad\_c2\_adv\_x.F}) |
184 |
|
|
185 |
|
$F_x$: {\bf uT} (argument) |
186 |
|
|
187 |
|
$U$: {\bf uTrans} (argument) |
188 |
|
|
189 |
|
$\tau$: {\bf tracer} (argument) |
190 |
|
|
191 |
|
{\em S/R GAD\_C2\_ADV\_Y} ({\em gad\_c2\_adv\_y.F}) |
192 |
|
|
193 |
|
$F_y$: {\bf vT} (argument) |
194 |
|
|
195 |
|
$V$: {\bf vTrans} (argument) |
196 |
|
|
197 |
|
$\tau$: {\bf tracer} (argument) |
198 |
|
|
199 |
|
{\em S/R GAD\_C2\_ADV\_R} ({\em gad\_c2\_adv\_r.F}) |
200 |
|
|
201 |
|
$F_r$: {\bf wT} (argument) |
202 |
|
|
203 |
|
$W$: {\bf rTrans} (argument) |
204 |
|
|
205 |
|
$\tau$: {\bf tracer} (argument) |
206 |
|
|
207 |
|
\end{minipage} } |
208 |
|
|
209 |
|
|
210 |
|
\subsection{Third order upwind bias advection} |
211 |
|
|
212 |
|
Upwind biased third order advection offers a relatively good |
213 |
|
compromise between accuracy and smoothness. It is not a ``positive'' |
214 |
|
scheme meaning false extrema are permitted but the amplitude of such |
215 |
|
are significantly reduced over the centered second order method. |
216 |
|
|
217 |
|
The third order upwind fluxes are discretized: |
218 |
|
\begin{eqnarray} |
219 |
|
F_x & = & U \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^i |
220 |
|
+ \frac{1}{2} |U| \delta_i \frac{1}{6} \delta_{ii} \tau \\ |
221 |
|
F_y & = & V \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^j |
222 |
|
+ \frac{1}{2} |V| \delta_j \frac{1}{6} \delta_{jj} \tau \\ |
223 |
|
F_r & = & W \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^k |
224 |
|
+ \frac{1}{2} |W| \delta_k \frac{1}{6} \delta_{kk} \tau |
225 |
|
\end{eqnarray} |
226 |
|
|
227 |
|
At boundaries, $\delta_{\hat{n}} \tau$ is set to zero allowing |
228 |
|
$\delta_{nn}$ to be evaluated. We are currently examing the accuracy |
229 |
|
of this boundary condition and the effect on the solution. |
230 |
|
|
231 |
|
|
232 |
|
\subsection{Centered fourth order advection} |
233 |
|
|
234 |
|
Centered fourth order advection is formally the most accurate scheme |
235 |
|
we have implemented and can be used to great effect in high resolution |
236 |
|
simultation where dynamical scales are well resolved. However, the |
237 |
|
scheme is noisey like the centered second order method and so must be |
238 |
|
used with some finite amount of diffusion. Bi-harmonic is recommended |
239 |
|
since it is more scale selective and less likely to diffuse away the |
240 |
|
well resolved gradient the fourth order scheme worked so hard to |
241 |
|
create. |
242 |
|
|
243 |
|
The centered fourth order fluxes are discretized: |
244 |
|
\begin{eqnarray} |
245 |
|
F_x & = & U \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^i \\ |
246 |
|
F_y & = & V \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^j \\ |
247 |
|
F_r & = & W \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^k |
248 |
|
\end{eqnarray} |
249 |
|
|
250 |
|
As for the third order scheme, the best discretization near boundaries |
251 |
|
is under investigation but currenlty $\delta_i \tau=0$ on a boundary. |
252 |
|
|
253 |
|
\subsection{First order upwind advection} |
254 |
|
|
255 |
|
Although the upwind scheme is the underlying scheme for the robust or |
256 |
|
non-linear methods given later, we haven't actually supplied this |
257 |
|
method for general use. It would be very diffusive and it is unlikely |
258 |
|
that it could ever produce more useful results than the positive |
259 |
|
higher order schemes. |
260 |
|
|
261 |
|
Upwind bias is introduced into many schemes using the {\em abs} |
262 |
|
function and is allows the first order upwind flux to be written: |
263 |
|
\begin{eqnarray} |
264 |
|
F_x & = & U \overline{ \tau }^i - \frac{1}{2} |U| \delta_i \tau \\ |
265 |
|
F_y & = & V \overline{ \tau }^j - \frac{1}{2} |V| \delta_j \tau \\ |
266 |
|
F_r & = & W \overline{ \tau }^k - \frac{1}{2} |W| \delta_k \tau |
267 |
|
\end{eqnarray} |
268 |
|
|
269 |
|
If for some reason, the above method is required, then the second |
270 |
|
order flux limiter scheme described later reduces to the above scheme |
271 |
|
if the limiter is set to zero. |
272 |
|
|
273 |
|
|
274 |
|
\section{Non-linear advection schemes} |
275 |
|
|
276 |
|
Non-linear advection schemes invoke non-linear interpolation and are |
277 |
|
widely used in computational fluid dynamics (non-linear does not refer |
278 |
|
to the non-linearity of the advection operator). The flux limited |
279 |
|
advection schemes belong to the class of finite volume methods which |
280 |
|
neatly ties into the spatial discretization of the model. |
281 |
|
|
282 |
|
When employing the flux limited schemes, first order upwind or |
283 |
|
direct-space-time method the time-stepping is switched to forward in |
284 |
|
time. |
285 |
|
|
286 |
|
\subsection{Second order flux limiters} |
287 |
|
|
288 |
|
The second order flux limiter method can be cast in several ways but |
289 |
|
is generally expressed in terms of other flux approximations. For |
290 |
|
example, in terms of a first order upwind flux and second order |
291 |
|
Lax-Wendroff flux, the limited flux is given as: |
292 |
|
\begin{equation} |
293 |
|
F = F_1 + \psi(r) F_{LW} |
294 |
|
\end{equation} |
295 |
|
where $\psi(r)$ is the limiter function, |
296 |
|
\begin{equation} |
297 |
|
F_1 = u \overline{\tau}^i - \frac{1}{2} |u| \delta_i \tau |
298 |
|
\end{equation} |
299 |
|
is the upwind flux, |
300 |
|
\begin{equation} |
301 |
|
F_{LW} = F_1 + \frac{|u|}{2} (1-c) \delta_i \tau |
302 |
|
\end{equation} |
303 |
|
is the Lax-Wendroff flux and $c = \frac{u \Delta t}{\Delta x}$ is the |
304 |
|
Courant (CFL) number. |
305 |
|
|
306 |
|
The limiter function, $\psi(r)$, takes the slope ratio |
307 |
|
\begin{eqnarray} |
308 |
|
r = \frac{ \tau_{i-1} - \tau_{i-2} }{ \tau_{i} - \tau_{i-1} } & \forall & u > 0 |
309 |
|
\\ |
310 |
|
r = \frac{ \tau_{i+1} - \tau_{i} }{ \tau_{i} - \tau_{i-1} } & \forall & u < 0 |
311 |
|
\end{eqnarray} |
312 |
|
as it's argument. There are many choices of limiter function but we |
313 |
|
only provide the Superbee limiter \cite{Roe85}: |
314 |
|
\begin{equation} |
315 |
|
\psi(r) = \max[0,\min[1,2r],\min[2,r]] |
316 |
|
\end{equation} |
317 |
|
|
318 |
|
|
319 |
|
\subsection{Third order direct space time} |
320 |
|
|
321 |
|
The direct-space-time method deals with space and time discretization |
322 |
|
together (other methods that treat space and time seperately are known |
323 |
|
collectively as the ``Method of Lines''). The Lax-Wendroff scheme |
324 |
|
falls into this category; it adds sufficient diffusion to a second |
325 |
|
order flux that the forward-in-time method is stable. The upwind |
326 |
|
biased third order DST scheme is: |
327 |
|
\begin{eqnarray} |
328 |
|
F = u \left( \tau_{i-1} |
329 |
|
+ d_0 (\tau_{i}-\tau_{i-1}) + d_1 (\tau_{i-1}-\tau_{i-2}) \right) |
330 |
|
& \forall & u > 0 \\ |
331 |
|
F = u \left( \tau_{i} |
332 |
|
- d_0 (\tau_{i}-\tau_{i-1}) - d_1 (\tau_{i+1}-\tau_{i}) \right) |
333 |
|
& \forall & u < 0 |
334 |
|
\end{eqnarray} |
335 |
|
where |
336 |
|
\begin{eqnarray} |
337 |
|
d_1 & = & \frac{1}{6} ( 2 - |c| ) ( 1 - |c| ) \\ |
338 |
|
d_2 & = & \frac{1}{6} ( 1 - |c| ) ( 1 + |c| ) |
339 |
|
\end{eqnarray} |
340 |
|
The coefficients $d_0$ and $d_1$ approach $1/3$ and $1/6$ repectively |
341 |
|
as the Courant number, $c$, vanishes. In this limit, the conventional |
342 |
|
third order upwind method is recovered. For finite Courant number, the |
343 |
|
deviations from the linear method are analogous to the diffusion added |
344 |
|
to centered second order advection in the Lax-Wendroff scheme. |
345 |
|
|
346 |
|
The DST3 method described above must be used in a forward-in-time |
347 |
|
manner and is stable for $0 \le |c| \le 1$. Although the scheme |
348 |
|
appears to be forward-in-time, it is in fact second order in time and |
349 |
|
the accuracy increases with the Courant number! For low Courant |
350 |
|
number, DST3 produces very similar results (indistinguishable in |
351 |
|
Fig.~\ref{fig:advect-1d-lo}) to the linear third order method but for |
352 |
|
large Courant number, where the linear upwind third order method is |
353 |
|
unstable, the scheme is extremely accurate |
354 |
|
(Fig.~\ref{fig:advect-1d-hi}) with only minor overshoots. |
355 |
|
|
356 |
|
\subsection{Third order direct space time with flux limiting} |
357 |
|
|
358 |
|
The overshoots in the DST3 method can be controlled with a flux limiter. |
359 |
|
The limited flux is written: |
360 |
|
\begin{equation} |
361 |
|
F = |
362 |
|
\frac{1}{2}(u+|u|)\left( \tau_{i-1} + \psi(r^+)(\tau_{i} - \tau_{i-1} )\right) |
363 |
|
+ |
364 |
|
\frac{1}{2}(u-|u|)\left( \tau_{i-1} + \psi(r^-)(\tau_{i} - \tau_{i-1} )\right) |
365 |
|
\end{equation} |
366 |
|
where |
367 |
|
\begin{eqnarray} |
368 |
|
r^+ & = & \frac{\tau_{i-1} - \tau_{i-2}}{\tau_{i} - \tau_{i-1}} \\ |
369 |
|
r^- & = & \frac{\tau_{i+1} - \tau_{i}}{\tau_{i} - \tau_{i-1}} |
370 |
|
\end{eqnarray} |
371 |
|
and the limiter is the Sweby limiter: |
372 |
|
\begin{equation} |
373 |
|
\psi(r) = \max[0, \min[\min(1,d_0+d_1r],\frac{1-c}{c}r ]] |
374 |
|
\end{equation} |
375 |
|
|
376 |
|
\subsection{Multi-dimensional advection} |
377 |
|
|
378 |
|
In many of the aforementioned advection schemes the behaviour in |
379 |
|
multiple dimensions is not necessarily as good as the one dimensional |
380 |
|
behaviour. For instance, a shape preserving monotonic scheme in one |
381 |
|
dimension can have severe shape distortion in two dimensions if the |
382 |
|
two components of horizontal fluxes are treated independently. There |
383 |
|
is a large body of literature on the subject dealing with this problem |
384 |
|
and among the fixes are operator and flux splitting methods, corner |
385 |
|
flux methods and more. We have adopted a variant on the standard |
386 |
|
splitting methods that allows the flux calculations to be implemented |
387 |
|
as if in one dimension: |
388 |
|
\begin{eqnarray} |
389 |
|
\tau^{n+1/3} & = & \tau^{n} |
390 |
|
- \Delta t \left( \frac{1}{\Delta x} \delta_i F^x(\tau^{n}) |
391 |
|
+ \tau^{n} \frac{1}{\Delta x} \delta_i u \right) \\ |
392 |
|
\tau^{n+2/3} & = & \tau^{n} |
393 |
|
- \Delta t \left( \frac{1}{\Delta y} \delta_j F^y(\tau^{n+1/3}) |
394 |
|
+ \tau^{n} \frac{1}{\Delta y} \delta_i v \right) \\ |
395 |
|
\tau^{n+3/3} & = & \tau^{n} |
396 |
|
- \Delta t \left( \frac{1}{\Delta r} \delta_k F^x(\tau^{n+2/3}) |
397 |
|
+ \tau^{n} \frac{1}{\Delta r} \delta_i w \right) |
398 |
|
\end{eqnarray} |
399 |
|
|
400 |
|
In order to incorporate this method into the general model algorithm, |
401 |
|
we compute the effective tendancy rather than update the tracer so |
402 |
|
that other terms such as diffusion are using the $n$ time-level and |
403 |
|
not the updated $n+3/3$ quantities: |
404 |
|
\begin{equation} |
405 |
|
G^{n+1/2}_{adv} = \frac{1}{\Delta t} ( \tau^{n+3/3} - \tau^{n} ) |
406 |
|
\end{equation} |
407 |
|
So that the over all time-stepping looks likes: |
408 |
|
\begin{equation} |
409 |
|
\tau^{n+1} = \tau^{n} + \Delta t \left( G^{n+1/2}_{adv} + G_{diff}(\tau^{n}) + G^{n}_{forcing} \right) |
410 |
|
\end{equation} |