/[MITgcm]/manual/s_algorithm/text/time_stepping.tex
ViewVC logotype

Diff of /manual/s_algorithm/text/time_stepping.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by adcroft, Wed Aug 8 16:15:16 2001 UTC revision 1.5 by jmc, Mon Sep 24 19:30:40 2001 UTC
# Line 8  at time $t = n \Delta t$ . We used the n Line 8  at time $t = n \Delta t$ . We used the n
8  when time interpolation is required to estimate the value of $\phi$  when time interpolation is required to estimate the value of $\phi$
9  at the time $n \Delta t$.  at the time $n \Delta t$.
10    
11  \section{Time Integration}  \section{Time integration}
12    
13  The discretization in time of the model equations (cf section I )  The discretization in time of the model equations (cf section I )
14  does not depend of the discretization in space of each  does not depend of the discretization in space of each
# Line 21  The continuous form of the model equatio Line 21  The continuous form of the model equatio
21    
22  \begin{eqnarray}  \begin{eqnarray}
23  \partial_t \theta & = & G_\theta  \partial_t \theta & = & G_\theta
24    \label{eq-tCsC-theta}
25  \\  \\
26  \partial_t S & = & G_s  \partial_t S & = & G_s
27    \label{eq-tCsC-salt}
28  \\  \\
29  b' & = & b'(\theta,S,r)  b' & = & b'(\theta,S,r)
30  \\  \\
31  \partial_r \phi'_{hyd} & = & -b'  \partial_r \phi'_{hyd} & = & -b'
32  \label{eq-r-split-hyd}  \label{eq-tCsC-hyd}
33  \\  \\
34  \partial_t \vec{\bf v}  \partial_t \vec{\bf v}
35  + {\bf \nabla}_r B_o \eta  + {\bf \nabla}_h b_s \eta
36  + \epsilon_{nh} {\bf \nabla}_r \phi'_{nh}  + \epsilon_{nh} {\bf \nabla}_h \phi'_{nh}
37  & = & \vec{\bf G}_{\vec{\bf v}}  & = & \vec{\bf G}_{\vec{\bf v}}
38  - {\bf \nabla}_r \phi'_{hyd}  - {\bf \nabla}_h \phi'_{hyd}
39  \label{eq-r-split-hmom}  \label{eq-tCsC-Hmom}
40  \\  \\
41  \epsilon_{nh} \frac {\partial{\dot{r}}}{\partial{t}}  \epsilon_{nh} \frac {\partial{\dot{r}}}{\partial{t}}
42  + \epsilon_{nh} \partial_r \phi'_{nh}  + \epsilon_{nh} \partial_r \phi'_{nh}
43  & = & \epsilon_{nh} G_{\dot{r}}  & = & \epsilon_{nh} G_{\dot{r}}
44  \label{eq-r-split-rdot}  \label{eq-tCsC-Vmom}
45  \\  \\
46  {\bf \nabla}_r \cdot \vec{\bf v} + \partial_r \dot{r}  {\bf \nabla}_h \cdot \vec{\bf v} + \partial_r \dot{r}
47  & = & 0  & = & 0
48  \label{eq-r-cont}  \label{eq-tCsC-cont}
49  \end{eqnarray}  \end{eqnarray}
50  where  where
51  \begin{eqnarray*}  \begin{eqnarray*}
52  G_\theta & = &  G_\theta & = &
53  - \vec{\bf v} \cdot {\bf \nabla}_r \theta + {\cal Q}_\theta  - \vec{\bf v} \cdot {\bf \nabla} \theta + {\cal Q}_\theta
54  \\  \\
55  G_S & = &  G_S & = &
56  - \vec{\bf v} \cdot {\bf \nabla}_r S + {\cal Q}_S  - \vec{\bf v} \cdot {\bf \nabla} S + {\cal Q}_S
57  \\  \\
58  \vec{\bf G}_{\vec{\bf v}}  \vec{\bf G}_{\vec{\bf v}}
59  & = &  & = &
60  - \vec{\bf v} \cdot {\bf \nabla}_r \vec{\bf v}  - \vec{\bf v} \cdot {\bf \nabla} \vec{\bf v}
61  - f \hat{\bf k} \wedge \vec{\bf v}  - f \hat{\bf k} \wedge \vec{\bf v}
62  + \vec{\cal F}_{\vec{\bf v}}  + \vec{\cal F}_{\vec{\bf v}}
63  \\  \\
64  G_{\dot{r}}  G_{\dot{r}}
65  & = &  & = &
66  - \vec{\bf v} \cdot {\bf \nabla}_r \dot{r}  - \vec{\bf v} \cdot {\bf \nabla} \dot{r}
67  + {\cal F}_{\dot{r}}  + {\cal F}_{\dot{r}}
68  \end{eqnarray*}  \end{eqnarray*}
69  The exact form of all the "{\it G}"s terms is described in the next  The exact form of all the "{\it G}"s terms is described in the next
# Line 70  all the RHS terms except the pressure / Line 72  all the RHS terms except the pressure /
72    
73  The switch $\epsilon_{nh}$ allows to activate the non hydrostatic  The switch $\epsilon_{nh}$ allows to activate the non hydrostatic
74  mode ($\epsilon_{nh}=1$) for the ocean model. Otherwise,  mode ($\epsilon_{nh}=1$) for the ocean model. Otherwise,
75  in the hydrostatic limit $\epsilon_{nh} = 0$ and the 3rd equation vanishes.  in the hydrostatic limit $\epsilon_{nh} = 0$
76    and equation \ref{eq-tCsC-Vmom} vanishes.
77    
78  The equation for $\eta$ is obtained by integrating the  The equation for $\eta$ is obtained by integrating the
79  continuity equation over the entire depth of the fluid,  continuity equation over the entire depth of the fluid,
80  from $R_{min}(x,y)$ up to $R_o(x,y)$  from $R_{fixed}(x,y)$ up to $R_o(x,y)$
81  (Linear free surface):  (Linear free surface):
82  \begin{displaymath}  \begin{eqnarray}
83  \epsilon_{fs} \partial_t \eta =  \epsilon_{fs} \partial_t \eta =
84  \left. \dot{r} \right|_{r=r_{surf}} + \epsilon_{fw} (P-E) =  \left. \dot{r} \right|_{r=r_{surf}} + \epsilon_{fw} (P-E) =
85  - {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o} \vec{\bf v} dr  - {\bf \nabla} \cdot \int_{R_{fixed}}^{R_o} \vec{\bf v} dr
86  + \epsilon_{fw} (P-E)  + \epsilon_{fw} (P-E)
87  \end{displaymath}  \label{eq-tCsC-eta}
88    \end{eqnarray}
89    
90  Where $\epsilon_{fs}$,$\epsilon_{fw}$ are two flags to  Where $\epsilon_{fs}$,$\epsilon_{fw}$ are two flags to
91  distinguish between a free-surface equation ($\epsilon_{fs}=1$)  distinguish between a free-surface equation ($\epsilon_{fs}=1$)
# Line 91  at the ocean surface (natural BC) ($\eps Line 95  at the ocean surface (natural BC) ($\eps
95  or not ($\epsilon_{fw} = 0$).  or not ($\epsilon_{fw} = 0$).
96    
97  The hydrostatic potential is found by  The hydrostatic potential is found by
98  integrating \ref{eq-r-split-hyd} with the boundary condition that  integrating \ref{eq-tCsC-hyd} with the boundary condition that
99  $\phi'_{hyd}(r=R_o) = 0$:  $\phi'_{hyd}(r=R_o) = 0$:
100  \begin{eqnarray*}  \begin{eqnarray*}
101  & &  & &
# Line 103  $\phi'_{hyd}(r=R_o) = 0$: Line 107  $\phi'_{hyd}(r=R_o) = 0$:
107  \phi'_{hyd}(x,y,r') = \int_{r'}^{R_o} b' dr  \phi'_{hyd}(x,y,r') = \int_{r'}^{R_o} b' dr
108  \end{eqnarray*}  \end{eqnarray*}
109    
110  \subsection{standard synchronous time stepping}  \subsection{General method}
111    
112    An overview of the general method is presented hereafter,
113    with explicit references to the Fortran code. This part
114    often refers to the discretized equations of the model
115    that are detailed in the following sections.
116    
117    The general algorithm consist in  a "predictor step" that computes
118    the forward tendencies ("G" terms") and all
119    the "first guess" values (star notation):
120    $\theta^*, S^*, \vec{\bf v}^*$ (and $\dot{r}^*$
121    in non-hydrostatic mode). This is done in the two routines
122    {\it THERMODYNAMICS} and {\it DYNAMICS}.
123    
124    Then the implicit terms that appear on the left hand side (LHS)
125    of equations \ref{eq-tDsC-theta} - \ref{eq-tDsC-cont},
126    are solved as follows:
127    Since implicit vertical diffusion and viscosity terms
128    are independent from the barotropic flow adjustment,
129    they are computed first, solving a 3 diagonal Nr x Nr linear system,
130    and incorporated at the end of the {\it THERMODYNAMICS} and
131    {\it DYNAMICS} routines.
132    Then the surface pressure and non hydrostatic pressure
133    are evaluated ({\it SOLVE\_FOR\_PRESSURE});
134    
135    Finally, within a "corrector step',
136    (routine {\it THE\_CORRECTION\_STEP})
137    the new values of $u,v,w,\theta,S$
138    are derived according to the above equations
139    (see details in II.1.3).
140    
141    At this point, the regular time step is over, but  
142    the correction step contains also other optional
143    adjustments such as convective adjustment algorithm, or filters
144    (zonal FFT filter, shapiro filter)
145    that applied on both momentum and tracer fields.
146    just before setting the $n+1$ new time step value.
147    
148    Since the pressure solver precision is of the order of
149    the "target residual" that could be lower than the
150    the computer truncation error, and also because some filters
151    might alter the divergence part of the flow field,
152    a final evaluation of the surface r anomaly $\eta^{n+1}$
153    is performed, according to \ref{eq-tDsC-eta} ({\it CALC\_EXACT\_ETA}).
154    This ensures a perfect volume conservation.
155    Note that there is no need for an equivalent Non-hydrostatic
156    "exact conservation" step, since W is already computed after
157    the filters applied.
158    
159    Regarding optional forcing terms (usually part of a "package"),
160    that account for a specific source or sink term (e.g.: condensation
161    as a sink of water vapor Q), they are generally incorporated
162    in the main algorithm as follows;
163    At the the beginning of the time step,
164    the additional tendencies are computed
165    as function of the present state (time step $n$) and external forcing ;
166    Then within the main part of model,
167    only those new tendencies are added to the model variables.
168    
169    [more details needed]\\
170    
171    The atmospheric physics follows this general scheme.
172    
173    [more about C\_grid, A\_grid conversion \& drag term]\\
174    
175    \subsection{Standard synchronous time stepping}
176    
177  In the standard formulation, the surface pressure is  In the standard formulation, the surface pressure is
178  evaluated at time step n+1 (implicit method).  evaluated at time step n+1 (implicit method).
# Line 112  as follows: Line 181  as follows:
181  \begin{eqnarray}  \begin{eqnarray}
182  \left[ 1 - \partial_r \kappa_v^\theta \partial_r \right]  \left[ 1 - \partial_r \kappa_v^\theta \partial_r \right]
183  \theta^{n+1} & = & \theta^*  \theta^{n+1} & = & \theta^*
184    \label{eq-tDsC-theta}
185  \\  \\
186  \left[ 1 - \partial_r \kappa_v^S \partial_r \right]  \left[ 1 - \partial_r \kappa_v^S \partial_r \right]
187  S^{n+1} & = & S^*  S^{n+1} & = & S^*
188    \label{eq-tDsC-salt}
189  \\  \\
190  %{b'}^{n} & = & b'(\theta^{n},S^{n},r)  %{b'}^{n} & = & b'(\theta^{n},S^{n},r)
191  %\partial_r {\phi'_{hyd}}^{n} & = & {-b'}^{n}  %\partial_r {\phi'_{hyd}}^{n} & = & {-b'}^{n}
192  %\\  %\\
193  {\phi'_{hyd}}^{n} & = & \int_{r'}^{R_o} b'(\theta^{n},S^{n},r) dr  {\phi'_{hyd}}^{n} & = & \int_{r'}^{R_o} b'(\theta^{n},S^{n},r) dr
194  \label{eq-rtd-hyd}  \label{eq-tDsC-hyd}
195  \\  \\
196  \vec{\bf v} ^{n+1}  \vec{\bf v} ^{n+1}
197  + \Delta t {\bf \nabla}_r B_o {\eta}^{n+1}  + \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
198  + \epsilon_{nh} \Delta t {\bf \nabla}_r {\phi'_{nh}}^{n+1}  + \epsilon_{nh} \Delta t {\bf \nabla} {\phi'_{nh}}^{n+1}
199  - \partial_r A_v \partial_r \vec{\bf v}^{n+1}  - \partial_r A_v \partial_r \vec{\bf v}^{n+1}
200  & = &  & = &
201  \vec{\bf v}^*  \vec{\bf v}^*
202  \label{eq-rtd-hmom}  \label{eq-tDsC-Hmom}
203  \\  \\
204  \epsilon_{fs} {\eta}^{n+1} + \Delta t  \epsilon_{fs} {\eta}^{n+1} + \Delta t
205  {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o} \vec{\bf v}^{n+1} dr  {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o} \vec{\bf v}^{n+1} dr
206  & = &  & = &
207      \epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta_t (P-E)^{n}      \epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta_t (P-E)^{n}
208  \nonumber  \nonumber
209  \\  \\
210  % = \epsilon_{fs} {\eta}^{n} & + & \epsilon_{fw} \Delta_t (P-E)^{n}  % = \epsilon_{fs} {\eta}^{n} & + & \epsilon_{fw} \Delta_t (P-E)^{n}
211  \label{eq-rtd-eta}  \label{eq-tDsC-eta}
212  \\  \\
213  \epsilon_{nh} \left( \dot{r} ^{n+1}  \epsilon_{nh} \left( \dot{r} ^{n+1}
214  + \Delta t \partial_r {\phi'_{nh}} ^{n+1}  + \Delta t \partial_r {\phi'_{nh}} ^{n+1}
215  \right)  \right)
216  & = & \epsilon_{nh} \dot{r}^*  & = & \epsilon_{nh} \dot{r}^*
217  \label{eq-rtd-rdot}  \label{eq-tDsC-Vmom}
218  \\  \\
219  {\bf \nabla}_r \cdot \vec{\bf v}^{n+1} + \partial_r \dot{r}^{n+1}  {\bf \nabla}_h \cdot \vec{\bf v}^{n+1} + \partial_r \dot{r}^{n+1}
220  & = & 0  & = & 0
221  \label{eq-rtd-cont}  \label{eq-tDsC-cont}
222  \end{eqnarray}  \end{eqnarray}
223  where  where
224  \begin{eqnarray}  \begin{eqnarray}
# Line 159  S ^{n} + \Delta t G_{S} ^{(n+1/2)} Line 230  S ^{n} + \Delta t G_{S} ^{(n+1/2)}
230  \\  \\
231  \vec{\bf v}^* & = &  \vec{\bf v}^* & = &
232  \vec{\bf v}^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}  \vec{\bf v}^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}
233  + \Delta t  {\bf \nabla}_r {\phi'_{hyd}}^{(n+1/2)}  + \Delta t  {\bf \nabla}_h {\phi'_{hyd}}^{(n+1/2)}
234  \\  \\
235  \dot{r}^* & = &  \dot{r}^* & = &
236  \dot{r} ^{n} + \Delta t G_{\dot{r}} ^{(n+1/2)}  \dot{r} ^{n} + \Delta t G_{\dot{r}} ^{(n+1/2)}
# Line 171  written separately here. Line 242  written separately here.
242    
243  To ensure a second order time discretization for both  To ensure a second order time discretization for both
244  momentum and tracer,  momentum and tracer,
245  The "G" terms are "extrapolated" forward in time  The "{\it G}" terms are "extrapolated" forward in time
246  (Adams Bashforth time stepping)  (Adams Bashforth time stepping)
247  from the values computed at time step $n$ and $n-1$  from the values computed at time step $n$ and $n-1$
248  to the time $n+1/2$:  to the time $n+1/2$:
# Line 183  In the standard non-stagger formulation, Line 254  In the standard non-stagger formulation,
254  the Adams-Bashforth time stepping is also applied to the  the Adams-Bashforth time stepping is also applied to the
255  hydrostatic (pressure / geo-) potential term $\nabla_h \Phi'_{hyd}$.  hydrostatic (pressure / geo-) potential term $\nabla_h \Phi'_{hyd}$.
256  Note that presently, this term is in fact incorporated to the  Note that presently, this term is in fact incorporated to the
257  $\vec{\bf G}_{\vec{\bf v}}$ arrays ({\it gU,gV}).  $\vec{\bf G}_{\vec{\bf v}}$ arrays ({\bf gU,gV}).
258    
259  \subsection{general method}  \subsection{Stagger baroclinic time stepping}
   
 The general algorithm consist in  a "predictor step" that computes  
 the forward tendencies ("G" terms") and all  
 the "first guess" values star notation):  
 $\vec{\bf v}^*, \theta^*, S^*$ (and $\dot{r}^*$  
 in non-hydrostatic mode). This is done in the routine {\it DYNAMICS}.  
   
 Then the implicit terms that appear here on the left hand side (LHS),  
 are solved as follows:  
 Since implicit vertical diffusion and viscosity terms  
 are independent from the barotropic flow adjustment,  
 they are computed first, solving a 3 diagonal Nr x Nr linear system,  
 and incorporated at the end of the {\it DYNAMICS} routines.  
 Then the surface pressure and non hydrostatic pressure  
 are evaluated ({\it SOLVE\_FOR\_PRESSURE});  
   
 Finally, within a "corrector step',  
 (routine {\it THE\_CORRECTION\_STEP})  
 the new values of $u,v,w,\theta,S$  
 are derived according to the above equations  
 (see details in II.1.3).  
   
 At this point, the regular time step is over, but    
 the correction step contains also other optional  
 adjustments such as convective adjustment algorithm, or filters  
 (zonal FFT filter, shapiro filter)  
 that applied on both momentum and tracer fields.  
 just before setting the $n+1$ new time step value.  
   
 Since the pressure solver precision is of the order of  
 the "target residual" that could be lower than the  
 the computer truncation error, and also because some filters  
 might alter the divergence part of the flow field,  
 a final evaluation of the surface r anomaly $\eta^{n+1}$  
 is performed, according to \ref{eq-rtd-eta} ({\it CALC\_EXACT\_ETA}).  
 This ensures a perfect volume conservation.  
 Note that there is no need for an equivalent Non-hydrostatic  
 "exact conservation" step, since W is already computed after  
 the filters applied.  
   
 optionnal forcing terms (package):\\  
 Regarding optional forcing terms (usually part of a "package"),  
 that a account for a specific source or sink term (e.g.: condensation  
 as a sink of water vapor Q), they are generally incorporated  
 in the main algorithm as follows;  
 At the the beginning of the time step,  
 the additionnal tendencies are computed  
 as function of the present state (time step $n$) and external forcing ;  
 Then within the main part of model,  
 only those new tendencies are added to the model variables.  
   
 [mode details needed]  
 The atmospheric physics follows this general scheme.  
   
 \subsection{stagger baroclinic time stepping}  
260    
261  An alternative is to evaluate $\phi'_{hyd}$ with the  An alternative is to evaluate $\phi'_{hyd}$ with the
262  new tracer fields, and step forward the momentum after.  new tracer fields, and step forward the momentum after.
263  This option is known as stagger baroclinic time stepping,  This option is known as stagger baroclinic time stepping,
264  since tracer and momentum are step forward in time one after the other.  since tracer and momentum are step forward in time one after the other.
265  It can be activated turning on a running flag parameter  It can be activated turning on a running flag parameter
266  {\it staggerTimeStep} in file "{\it data}").  {\bf staggerTimeStep} in file "{\it data}").
267    
268  The main advantage of this time stepping compared to a synchronous one,  The main advantage of this time stepping compared to a synchronous one,
269  is a better stability, specially regarding internal gravity waves,  is a better stability, specially regarding internal gravity waves,
# Line 284  And Line 300  And
300  %\\  %\\
301  %\partial_r {\phi'_{hyd}}^{n+1/2} & = & {-b'}^{n+1/2}  %\partial_r {\phi'_{hyd}}^{n+1/2} & = & {-b'}^{n+1/2}
302  {\phi'_{hyd}}^{n+1/2} & = & \int_{r'}^{R_o} b'(\theta^{n+1/2},S^{n+1/2},r) dr  {\phi'_{hyd}}^{n+1/2} & = & \int_{r'}^{R_o} b'(\theta^{n+1/2},S^{n+1/2},r) dr
303  %\label{eq-rtd-hyd}  %\label{eq-tDsC-hyd}
304  \\  \\
305  \vec{\bf v} ^{n+1}  \vec{\bf v} ^{n+1}
306  + \Delta t {\bf \nabla}_r B_o {\eta}^{n+1}  + \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
307  + \epsilon_{nh} \Delta t {\bf \nabla}_r {\phi'_{nh}}^{n+1}  + \epsilon_{nh} \Delta t {\bf \nabla}_h {\phi'_{nh}}^{n+1}
308  - \partial_r A_v \partial_r \vec{\bf v}^{n+1}  - \partial_r A_v \partial_r \vec{\bf v}^{n+1}
309  & = &  & = &
310  \vec{\bf v}^*  \vec{\bf v}^*
311  %\label{eq-rtd-hmom}  %\label{eq-tDsC-Hmom}
312  \\  \\
313  \epsilon_{fs} {\eta}^{n+1} + \Delta t  \epsilon_{fs} {\eta}^{n+1} + \Delta t
314  {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o} \vec{\bf v}^{n+1} dr  {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o} \vec{\bf v}^{n+1} dr
315  & = &  & = &
316  \epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta_t (P-E)^{n}  \epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta_t (P-E)^{n}
317  \\  \\
# Line 303  And Line 319  And
319  + \Delta t \partial_r {\phi'_{nh}} ^{n+1}  + \Delta t \partial_r {\phi'_{nh}} ^{n+1}
320  \right)  \right)
321  & = & \epsilon_{nh} \dot{r}^*  & = & \epsilon_{nh} \dot{r}^*
322  %\label{eq-rtd-rdot}  %\label{eq-tDsC-Vmom}
323  \\  \\
324  {\bf \nabla}_r \cdot \vec{\bf v}^{n+1} + \partial_r \dot{r}^{n+1}  {\bf \nabla}_h \cdot \vec{\bf v}^{n+1} + \partial_r \dot{r}^{n+1}
325  & = & 0  & = & 0
326  %\label{eq-rtd-cont}  %\label{eq-tDsC-cont}
327  \end{eqnarray*}  \end{eqnarray*}
328  with  with
329  \begin{eqnarray*}  \begin{eqnarray*}
330  \vec{\bf v}^* & = &  \vec{\bf v}^* & = &
331  \vec{\bf v} ^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}  \vec{\bf v} ^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}
332  + \Delta t  {\bf \nabla}_r {\phi'_{hyd}}^{n+1/2}  + \Delta t  {\bf \nabla}_h {\phi'_{hyd}}^{n+1/2}
333  \\  \\
334  \dot{r}^* & = &  \dot{r}^* & = &
335  \dot{r} ^{n} + \Delta t G_{\dot{r}} ^{(n+1/2)}  \dot{r} ^{n} + \Delta t G_{\dot{r}} ^{(n+1/2)}
# Line 321  with Line 337  with
337    
338  %---------------------------------------------------------------------  %---------------------------------------------------------------------
339    
340  \subsection{surface pressure}  \subsection{Surface pressure}
341    
342  Substituting \ref{eq-rtd-hmom} into \ref{eq-rtd-cont}, assuming  Substituting \ref{eq-tDsC-Hmom} into \ref{eq-tDsC-cont}, assuming
343  $\epsilon_{nh} = 0$ yields a Helmholtz equation for ${\eta}^{n+1}$:  $\epsilon_{nh} = 0$ yields a Helmholtz equation for ${\eta}^{n+1}$:
344  $$  \begin{eqnarray}
345  \epsilon_{fs} {\eta}^{n+1} -  \epsilon_{fs} {\eta}^{n+1} -
346  {\bf \nabla}_r \cdot \Delta t^2 (R_o-R_{min})  {\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{fixed})
347  {\bf \nabla}_r B_o {\eta}^{n+1}  {\bf \nabla}_h b_s {\eta}^{n+1}
348  = {\eta}^*  = {\eta}^*
349  \label{solve_2d}  \label{eq-solve2D}
350  $$  \end{eqnarray}
351  where  where
352  $$  \begin{eqnarray}
353  {\eta}^* = \epsilon_{fs} \: {\eta}^{n} -  {\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
354  \Delta t {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o} \vec{\bf v}^* dr  \Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o} \vec{\bf v}^* dr
355  \: + \: \epsilon_{fw} \Delta_t (P-E)^{n}  \: + \: \epsilon_{fw} \Delta_t (P-E)^{n}
356  $$  \label{eq-solve2D_rhs}
357    \end{eqnarray}
358    
359  Once ${\eta}^{n+1}$ has been found substituting into \ref{eq-rtd-hmom}  Once ${\eta}^{n+1}$ has been found substituting into \ref{eq-tDsC-Hmom}
360  would yield $\vec{\bf v}^{n+1}$ if the model is hydrostatic  would yield $\vec{\bf v}^{n+1}$ if the model is hydrostatic
361  ($\epsilon_{nh}=0$):  ($\epsilon_{nh}=0$):
362  $$  $$
363  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}
364  - \Delta t {\bf \nabla}_r B_o {\eta}^{n+1}  - \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
365  $$  $$
366    
367  This is known as the correction step. However, when the model is  This is known as the correction step. However, when the model is
368  non-hydrostatic ($\epsilon_{nh}=1$) we need an additional step and an  non-hydrostatic ($\epsilon_{nh}=1$) we need an additional step and an
369  additional equation for $\phi'_{nh}$. This is obtained by  additional equation for $\phi'_{nh}$. This is obtained by
370  substituting \ref{eq-rtd-hmom} and \ref{eq-rtd-rdot} into  substituting \ref{eq-tDsC-Hmom} and \ref{eq-tDsC-Vmom} into
371  \ref{eq-rtd-cont}:  \ref{eq-tDsC-cont}:
372  \begin{equation}  \begin{equation}
373  \left[ {\bf \nabla}_r^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1}  \left[ {\bf \nabla}_h^2 + \partial_{rr} \right] {\phi'_{nh}}^{n+1}
374  = \frac{1}{\Delta t} \left(  = \frac{1}{\Delta t} \left(
375  {\bf \nabla}_r \cdot \vec{\bf v}^{**} + \partial_r \dot{r}^* \right)  {\bf \nabla}_h \cdot \vec{\bf v}^{**} + \partial_r \dot{r}^* \right)
376  \end{equation}  \end{equation}
377  where  where
378  \begin{displaymath}  \begin{displaymath}
379  \vec{\bf v}^{**} = \vec{\bf v}^* - \Delta t {\bf \nabla}_r B_o {\eta}^{n+1}  \vec{\bf v}^{**} = \vec{\bf v}^* - \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
380  \end{displaymath}  \end{displaymath}
381  Note that $\eta^{n+1}$ is also used to update the second RHS term  Note that $\eta^{n+1}$ is also used to update the second RHS term
382  $\partial_r \dot{r}^* $ since  $\partial_r \dot{r}^* $ since
# Line 369  is evaluated as $(\eta^{n+1} - \eta^n) / Line 386  is evaluated as $(\eta^{n+1} - \eta^n) /
386  Finally, the horizontal velocities at the new time level are found by:  Finally, the horizontal velocities at the new time level are found by:
387  \begin{equation}  \begin{equation}
388  \vec{\bf v}^{n+1} = \vec{\bf v}^{**}  \vec{\bf v}^{n+1} = \vec{\bf v}^{**}
389  - \epsilon_{nh} \Delta t {\bf \nabla}_r {\phi'_{nh}}^{n+1}  - \epsilon_{nh} \Delta t {\bf \nabla}_h {\phi'_{nh}}^{n+1}
390  \end{equation}  \end{equation}
391  and the vertical velocity is found by integrating the continuity  and the vertical velocity is found by integrating the continuity
392  equation vertically.  equation vertically.
# Line 381  without any consequence on the solution. Line 398  without any consequence on the solution.
398  Regarding the implementation, all those computation are done  Regarding the implementation, all those computation are done
399  within the routine {\it SOLVE\_FOR\_PRESSURE} and its dependent  within the routine {\it SOLVE\_FOR\_PRESSURE} and its dependent
400  {\it CALL}s.  {\it CALL}s.
401  The standard method to solve the 2D elliptic problem (\ref{solve_2d})  The standard method to solve the 2D elliptic problem (\ref{eq-solve2D})
402  uses the conjugate gradient method (routine {\it CG2D}); The  uses the conjugate gradient method (routine {\it CG2D}); The
403  the solver matrix and conjugate gradient operator are only function  the solver matrix and conjugate gradient operator are only function
404  of the discretized domain and are therefore evaluated separately,  of the discretized domain and are therefore evaluated separately,
# Line 414  In the code, $\beta,\gamma$ are defined Line 431  In the code, $\beta,\gamma$ are defined
431  {\it implicSurfPress}, {\it implicDiv2DFlow}. They are read from  {\it implicSurfPress}, {\it implicDiv2DFlow}. They are read from
432  the main data file "{\it data}" and are set by default to 1,1.  the main data file "{\it data}" and are set by default to 1,1.
433    
434  Equations \ref{eq-rtd-hmom} and \ref{eq-rtd-eta} are modified as follows:  Equations \ref{eq-tDsC-Hmom} and \ref{eq-tDsC-eta} are modified as follows:
435  $$  $$
436  \frac{ \vec{\bf v}^{n+1} }{ \Delta t }  \frac{ \vec{\bf v}^{n+1} }{ \Delta t }
437  + {\bf \nabla}_r B_o [ \beta {\eta}^{n+1} + (1-\beta) {\eta}^{n} ]  + {\bf \nabla}_h b_s [ \beta {\eta}^{n+1} + (1-\beta) {\eta}^{n} ]
438  + \epsilon_{nh} {\bf \nabla}_r {\phi'_{nh}}^{n+1}  + \epsilon_{nh} {\bf \nabla}_h {\phi'_{nh}}^{n+1}
439   = \frac{ \vec{\bf v}^* }{ \Delta t }   = \frac{ \vec{\bf v}^* }{ \Delta t }
440  $$  $$
441  $$  $$
442  \epsilon_{fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t}  \epsilon_{fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t}
443  + {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o}  + {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o}
444  [ \gamma \vec{\bf v}^{n+1} + (1-\gamma) \vec{\bf v}^{n}] dr  [ \gamma \vec{\bf v}^{n+1} + (1-\gamma) \vec{\bf v}^{n}] dr
445  = \epsilon_{fw} (P-E)  = \epsilon_{fw} (P-E)
446  $$  $$
# Line 431  where: Line 448  where:
448  \begin{eqnarray*}  \begin{eqnarray*}
449  \vec{\bf v}^* & = &  \vec{\bf v}^* & = &
450  \vec{\bf v} ^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}  \vec{\bf v} ^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}
451  + (\beta-1) \Delta t {\bf \nabla}_r B_o {\eta}^{n}  + (\beta-1) \Delta t {\bf \nabla}_h b_s {\eta}^{n}
452  + \Delta t {\bf \nabla}_r {\phi'_{hyd}}^{(n+1/2)}  + \Delta t {\bf \nabla}_h {\phi'_{hyd}}^{(n+1/2)}
453  \\  \\
454  {\eta}^* & = &  {\eta}^* & = &
455  \epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta t (P-E)  \epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta t (P-E)
456  - \Delta t {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o}  - \Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o}
457  [ \gamma \vec{\bf v}^* + (1-\gamma) \vec{\bf v}^{n}] dr  [ \gamma \vec{\bf v}^* + (1-\gamma) \vec{\bf v}^{n}] dr
458  \end{eqnarray*}  \end{eqnarray*}
459  \\  \\
# Line 444  In the hydrostatic case ($\epsilon_{nh}= Line 461  In the hydrostatic case ($\epsilon_{nh}=
461  this allow to find ${\eta}^{n+1}$, according to:  this allow to find ${\eta}^{n+1}$, according to:
462  $$  $$
463  \epsilon_{fs} {\eta}^{n+1} -  \epsilon_{fs} {\eta}^{n+1} -
464  {\bf \nabla}_r \cdot \beta\gamma \Delta t^2 B_o (R_o - R_{min})  {\bf \nabla}_h \cdot \beta\gamma \Delta t^2 b_s (R_o - R_{fixed})
465  {\bf \nabla}_r {\eta}^{n+1}  {\bf \nabla}_h {\eta}^{n+1}
466  = {\eta}^*  = {\eta}^*
467  $$  $$
468  and then to compute (correction step):  and then to compute (correction step):
469  $$  $$
470  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}
471  - \beta \Delta t {\bf \nabla}_r B_o {\eta}^{n+1}  - \beta \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
472  $$  $$
473    
474  The non-hydrostatic part is solved as described previously.  The non-hydrostatic part is solved as described previously.
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22