/[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.28 by jmc, Mon Aug 30 23:09:19 2010 UTC revision 1.34 by jmc, Mon Sep 25 18:11:06 2017 UTC
# Line 90  temporary.} Line 90  temporary.}
90    
91  \begin{figure}  \begin{figure}
92  \begin{center} \fbox{ \begin{minipage}{4.5in} \begin{tabbing}  \begin{center} \fbox{ \begin{minipage}{4.5in} \begin{tabbing}
93  aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill  aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill
94  \filelink{FORWARD\_STEP}{model-src-forward_step.F} \\  \filelink{FORWARD\_STEP}{model-src-forward_step.F} \\
95  \> DYNAMICS \\  \> DYNAMICS \\
96  \>\> TIMESTEP \` $u^*$,$v^*$ (\ref{eq:ustar-rigid-lid},\ref{eq:vstar-rigid-lid}) \\  \>\> TIMESTEP \` $u^*$,$v^*$ (\ref{eq:ustar-rigid-lid},\ref{eq:vstar-rigid-lid}) \\
# Line 101  aaa \= aaa \= aaa \= aaa \= aaa \= aaa \ Line 101  aaa \= aaa \= aaa \= aaa \= aaa \= aaa \
101  \>\> CALC\_GRAD\_PHI\_SURF \` $\nabla \eta^{n+1}$ \\  \>\> CALC\_GRAD\_PHI\_SURF \` $\nabla \eta^{n+1}$ \\
102  \>\> CORRECTION\_STEP \` $u^{n+1}$,$v^{n+1}$ (\ref{eq:un+1-rigid-lid},\ref{eq:vn+1-rigid-lid})  \>\> CORRECTION\_STEP \` $u^{n+1}$,$v^{n+1}$ (\ref{eq:un+1-rigid-lid},\ref{eq:vn+1-rigid-lid})
103  \end{tabbing} \end{minipage} } \end{center}  \end{tabbing} \end{minipage} } \end{center}
104  \caption{Calling tree for the pressure method algorithm  \caption{Calling tree for the pressure method algorithm
105    (\filelink{FORWARD\_STEP}{model-src-forward_step.F})}    (\filelink{FORWARD\_STEP}{model-src-forward_step.F})}
106  \label{fig:call-tree-pressure-method}  \label{fig:call-tree-pressure-method}
107  \end{figure}  \end{figure}
# Line 188  stepping forward $u^n$ and $v^n$ to $u^{ Line 188  stepping forward $u^n$ and $v^n$ to $u^{
188  \item  \item
189  the vertical integration, $H \widehat{u^*}$ and $H  the vertical integration, $H \widehat{u^*}$ and $H
190  \widehat{v^*}$, divergence and inversion of the elliptic operator in  \widehat{v^*}$, divergence and inversion of the elliptic operator in
191  equation \ref{eq:elliptic} is coded in  equation \ref{eq:elliptic} is coded in
192  \filelink{SOLVE\_FOR\_PRESSURE()}{model-src-solve_for_pressure.F}  \filelink{SOLVE\_FOR\_PRESSURE()}{model-src-solve_for_pressure.F}
193  \item  \item
194  finally, the new flow field at time level $n+1$ given by equations  finally, the new flow field at time level $n+1$ given by equations
195  \ref{eq:un+1-rigid-lid} and \ref{eq:vn+1-rigid-lid} is calculated in  \ref{eq:un+1-rigid-lid} and \ref{eq:vn+1-rigid-lid} is calculated in
196  \filelink{CORRECTION\_STEP()}{model-src-correction_step.F}.  \filelink{CORRECTION\_STEP()}{model-src-correction_step.F}.
197  \end{itemize}  \end{itemize}
198  The calling tree for these routines is given in  The calling tree for these routines is given in
# Line 200  Fig.~\ref{fig:call-tree-pressure-method} Line 200  Fig.~\ref{fig:call-tree-pressure-method}
200    
201    
202  %\paragraph{Need to discuss implicit viscosity somewhere:}  %\paragraph{Need to discuss implicit viscosity somewhere:}
203  In general, the horizontal momentum time-stepping can contain some terms  In general, the horizontal momentum time-stepping can contain some terms
204  that are treated implicitly in time,  that are treated implicitly in time,
205  such as the vertical viscosity when using the backward time-stepping scheme  such as the vertical viscosity when using the backward time-stepping scheme
206  (\varlink{implicitViscosity}{implicitViscosity} {\it =.TRUE.}).  (\varlink{implicitViscosity}{implicitViscosity} {\it =.TRUE.}).
207  The method used to solve those implicit terms is provided in  The method used to solve those implicit terms is provided in
208  section \ref{sec:implicit-backward-stepping}, and modifies  section \ref{sec:implicit-backward-stepping}, and modifies
209  equations \ref{eq:discrete-time-u} and \ref{eq:discrete-time-v} to  equations \ref{eq:discrete-time-u} and \ref{eq:discrete-time-v} to
210  give:  give:
# Line 263  As for the rigid-lid pressure method, eq Line 263  As for the rigid-lid pressure method, eq
263  \begin{eqnarray}  \begin{eqnarray}
264  u^{*} & = & u^{n} + \Delta t G_u^{(n+1/2)} \label{eq:ustar-backward-free-surface} \\  u^{*} & = & u^{n} + \Delta t G_u^{(n+1/2)} \label{eq:ustar-backward-free-surface} \\
265  v^{*} & = & v^{n} + \Delta t G_v^{(n+1/2)} \label{eq:vstar-backward-free-surface} \\  v^{*} & = & v^{n} + \Delta t G_v^{(n+1/2)} \label{eq:vstar-backward-free-surface} \\
266  \eta^* & = & \epsilon_{fs} \left( \eta^{n} + \Delta t (P-E) \right)- \Delta t  \eta^* & = & \epsilon_{fs} \left( \eta^{n} + \Delta t (P-E) \right)
267    \partial_x H \widehat{u^{*}}           - \Delta t \left( \partial_x H \widehat{u^{*}}
268  + \partial_y H \widehat{v^{*}}                           + \partial_y H \widehat{v^{*}} \right)
269  \\  \\
270    \partial_x g H \partial_x \eta^{n+1}    \partial_x g H \partial_x \eta^{n+1}
271  & + & \partial_y g H \partial_y \eta^{n+1}  & + & \partial_y g H \partial_y \eta^{n+1}
272   - \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2}   - \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2}
273   =   =
274  - \frac{\eta^*}{\Delta t^2}  - \frac{\eta^*}{\Delta t^2}
275  \label{eq:elliptic-backward-free-surface}  \label{eq:elliptic-backward-free-surface}
276  \\  \\
# Line 286  recovered. However, the implicit treatme Line 286  recovered. However, the implicit treatme
286  the flow to be divergent and for the surface pressure/elevation to  the flow to be divergent and for the surface pressure/elevation to
287  respond on a finite time-scale (as opposed to instantly). To recover  respond on a finite time-scale (as opposed to instantly). To recover
288  the rigid-lid formulation, we introduced a switch-like parameter,  the rigid-lid formulation, we introduced a switch-like parameter,
289  $\epsilon_{fs}$ (\varlink{freesurfFac}{freesurfFac}),  $\epsilon_{fs}$ (\varlink{freesurfFac}{freesurfFac}),
290  which selects between the free-surface and rigid-lid;  which selects between the free-surface and rigid-lid;
291  $\epsilon_{fs}=1$ allows the free-surface to evolve; $\epsilon_{fs}=0$  $\epsilon_{fs}=1$ allows the free-surface to evolve; $\epsilon_{fs}=0$
292  imposes the rigid-lid. The evolution in time and location of variables  imposes the rigid-lid. The evolution in time and location of variables
# Line 312  tracers (see section \ref{sec:tracer-adv Line 312  tracers (see section \ref{sec:tracer-adv
312    
313  \begin{figure}  \begin{figure}
314  \begin{center} \fbox{ \begin{minipage}{4.5in} \begin{tabbing}  \begin{center} \fbox{ \begin{minipage}{4.5in} \begin{tabbing}
315  aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill  aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill
316  FORWARD\_STEP \\  FORWARD\_STEP \\
317  \> THERMODYNAMICS \\  \> THERMODYNAMICS \\
318  \>\> CALC\_GT \\  \>\> CALC\_GT \\
# Line 358  simpler to include these terms and this Line 358  simpler to include these terms and this
358  and forcing evolves smoothly. Problems can, and do, arise when forcing  and forcing evolves smoothly. Problems can, and do, arise when forcing
359  or motions are high frequency and this corresponds to a reduced  or motions are high frequency and this corresponds to a reduced
360  stability compared to a simple forward time-stepping of such terms.  stability compared to a simple forward time-stepping of such terms.
361  The model offers the possibility to leave the tracer and momentum  The model offers the possibility to leave the tracer and momentum
362  forcing terms and the dissipation terms outside the  forcing terms and the dissipation terms outside the
363  Adams-Bashforth extrapolation, by turning off the logical flags  Adams-Bashforth extrapolation, by turning off the logical flags
364  \varlink{forcing\_In\_AB}{forcing_In_AB}  \varlink{forcing\_In\_AB}{forcing_In_AB}
365  (parameter file {\em data}, namelist {\em PARM01}, default value = True).  (parameter file {\em data}, namelist {\em PARM01}, default value = True).
366  (\varlink{tracForcingOutAB}{tracForcingOutAB}, default=0,  (\varlink{tracForcingOutAB}{tracForcingOutAB}, default=0,
# Line 403  The right column represents the damping Line 403  The right column represents the damping
403  \end{rawhtml}  \end{rawhtml}
404    
405  Vertical diffusion and viscosity can be treated implicitly in time  Vertical diffusion and viscosity can be treated implicitly in time
406  using the backward method which is an intrinsic scheme.  using the backward method which is an intrinsic scheme.
407  Recently, the option to treat the vertical advection  Recently, the option to treat the vertical advection
408  implicitly has been added, but not yet tested; therefore,  implicitly has been added, but not yet tested; therefore,
409  the description hereafter is limited to diffusion and viscosity.  the description hereafter is limited to diffusion and viscosity.
410  For tracers,  For tracers,
411  the time discretized equation is:  the time discretized equation is:
# Line 425  using the Adams-Bashforth method as desc Line 425  using the Adams-Bashforth method as desc
425  \end{eqnarray}  \end{eqnarray}
426  where ${\cal L}_\tau^{-1}$ is the inverse of the operator  where ${\cal L}_\tau^{-1}$ is the inverse of the operator
427  \begin{equation}  \begin{equation}
428  {\cal L} = \left[ 1 + \Delta t \partial_r \kappa_v \partial_r \right]  {\cal L}_\tau = \left[ 1 + \Delta t \partial_r \kappa_v \partial_r \right]
429  \end{equation}  \end{equation}
430  Equation \ref{eq:taustar-implicit} looks exactly as \ref{eq:taustar}  Equation \ref{eq:taustar-implicit} looks exactly as \ref{eq:taustar}
431  while \ref{eq:tau-n+1-implicit} involves an operator or matrix  while \ref{eq:tau-n+1-implicit} involves an operator or matrix
# Line 470  is solved to yield the state variables a Line 470  is solved to yield the state variables a
470    
471  \begin{figure}  \begin{figure}
472  \begin{center} \fbox{ \begin{minipage}{4.7in} \begin{tabbing}  \begin{center} \fbox{ \begin{minipage}{4.7in} \begin{tabbing}
473  aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill  aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill
474  FORWARD\_STEP \\  FORWARD\_STEP \\
475  \>\> EXTERNAL\_FIELDS\_LOAD\\  \>\> EXTERNAL\_FIELDS\_LOAD\\
476  \>\> DO\_ATMOSPHERIC\_PHYS \\  \>\> DO\_ATMOSPHERIC\_PHYS \\
# Line 501  FORWARD\_STEP \\ Line 501  FORWARD\_STEP \\
501  \end{tabbing} \end{minipage} } \end{center}  \end{tabbing} \end{minipage} } \end{center}
502  \caption{  \caption{
503  Calling tree for the overall synchronous algorithm using  Calling tree for the overall synchronous algorithm using
504  Adams-Bashforth time-stepping.  Adams-Bashforth time-stepping.
505  The place where the model geometry  The place where the model geometry
506  ({\bf hFac} factors) is updated is added here but is only relevant  ({\bf hFac} factors) is updated is added here but is only relevant
507  for the non-linear free-surface algorithm.  for the non-linear free-surface algorithm.
508  For completeness, the external forcing,  For completeness, the external forcing,
509  ocean and atmospheric physics have been added, although they are mainly  ocean and atmospheric physics have been added, although they are mainly
510  optional}  optional}
511  \label{fig:call-tree-adams-bashforth-sync}  \label{fig:call-tree-adams-bashforth-sync}
512  \end{figure}  \end{figure}
# Line 536  G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil Line 536  G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil
536  \label{eq:vstar-sync} \\  \label{eq:vstar-sync} \\
537  \vec{\bf v}^{**} & = & {\cal L}_{\vec{\bf v}}^{-1} ( \vec{\bf v}^* )  \vec{\bf v}^{**} & = & {\cal L}_{\vec{\bf v}}^{-1} ( \vec{\bf v}^* )
538  \label{eq:vstarstar-sync} \\  \label{eq:vstarstar-sync} \\
539  \eta^* & = & \epsilon_{fs} \left( \eta^{n} + \Delta t (P-E) \right)- \Delta t  \eta^* & = & \epsilon_{fs} \left( \eta^{n} + \Delta t (P-E) \right)- \Delta t
540    \nabla \cdot H \widehat{ \vec{\bf v}^{**} }    \nabla \cdot H \widehat{ \vec{\bf v}^{**} }
541  \label{eq:nstar-sync} \\  \label{eq:nstar-sync} \\
542  \nabla \cdot g H \nabla \eta^{n+1} & - & \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2}  \nabla \cdot g H \nabla \eta^{n+1} & - & \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2}
543  ~ = ~ - \frac{\eta^*}{\Delta t^2}  ~ = ~ - \frac{\eta^*}{\Delta t^2}
544  \label{eq:elliptic-sync} \\  \label{eq:elliptic-sync} \\
545  \vec{\bf v}^{n+1} & = & \vec{\bf v}^{*} - \Delta t g \nabla \eta^{n+1}  \vec{\bf v}^{n+1} & = & \vec{\bf v}^{**} - \Delta t g \nabla \eta^{n+1}
546  \label{eq:v-n+1-sync}  \label{eq:v-n+1-sync}
547  \end{eqnarray}  \end{eqnarray}
548  Fig.~\ref{fig:adams-bashforth-sync} illustrates the location of  Fig.~\ref{fig:adams-bashforth-sync} illustrates the location of
# Line 577  time-step. The corresponding calling tre Line 577  time-step. The corresponding calling tre
577  \caption{  \caption{
578  A schematic of the explicit Adams-Bashforth and implicit time-stepping  A schematic of the explicit Adams-Bashforth and implicit time-stepping
579  phases of the algorithm but with staggering in time of thermodynamic  phases of the algorithm but with staggering in time of thermodynamic
580  variables with the flow.  variables with the flow.
581  Explicit momentum tendencies are evaluated at time level $n-1/2$ as a  Explicit momentum tendencies are evaluated at time level $n-1/2$ as a
582  function of the flow field at that time level $n-1/2$.  function of the flow field at that time level $n-1/2$.
583  The explicit tendency from the previous time level, $n-3/2$, is used to  The explicit tendency from the previous time level, $n-3/2$, is used to
584  extrapolate tendencies to $n$ (dashed arrow).  extrapolate tendencies to $n$ (dashed arrow).
585  The hydrostatic pressure/geo-potential $\phi_{hyd}$ is evaluated directly  The hydrostatic pressure/geo-potential $\phi_{hyd}$ is evaluated directly
586  at time level $n$ (vertical arrows) and used with the extrapolated tendencies  at time level $n$ (vertical arrows) and used with the extrapolated tendencies
587  to step forward the flow variables from $n-1/2$ to $n+1/2$ (solid arc-arrow).  to step forward the flow variables from $n-1/2$ to $n+1/2$ (solid arc-arrow).
588  The implicit-in-time operator ${\cal L}_{\bf u,v}$ (vertical arrows) is  The implicit-in-time operator ${\cal L}_{\bf u,v}$ (vertical arrows) is
589  then applied to the previous estimation of the the flow field ($*$-variables)  then applied to the previous estimation of the the flow field ($*$-variables)
590  and yields to the two velocity components $u,v$ at time level $n+1/2$.  and yields to the two velocity components $u,v$ at time level $n+1/2$.
591  These are then used to calculate the advection term (dashed arc-arrow)  These are then used to calculate the advection term (dashed arc-arrow)
592  of the thermo-dynamics tendencies at time step $n$.  of the thermo-dynamics tendencies at time step $n$.
593  The extrapolated thermodynamics tendency, from time level $n-1$ and $n$  The extrapolated thermodynamics tendency, from time level $n-1$ and $n$
594  to $n+1/2$, allows thermodynamic variables to be stably integrated  to $n+1/2$, allows thermodynamic variables to be stably integrated
595  forward-in-time (solid arc-arrow) up to time level $n+1$.  forward-in-time (solid arc-arrow) up to time level $n+1$.
596  }  }
597  \label{fig:adams-bashforth-staggered}  \label{fig:adams-bashforth-staggered}
# Line 603  circumstance, it is more efficient to st Line 603  circumstance, it is more efficient to st
603  thermodynamic variables with the flow  thermodynamic variables with the flow
604  variables. Fig.~\ref{fig:adams-bashforth-staggered} illustrates the  variables. Fig.~\ref{fig:adams-bashforth-staggered} illustrates the
605  staggering and algorithm. The key difference between this and  staggering and algorithm. The key difference between this and
606  Fig.~\ref{fig:adams-bashforth-sync} is that the thermodynamic variables  Fig.~\ref{fig:adams-bashforth-sync} is that the thermodynamic variables
607  are solved after the dynamics, using the recently updated flow field.  are solved after the dynamics, using the recently updated flow field.
608  This essentially allows the gravity wave terms to leap-frog in  This essentially allows the gravity wave terms to leap-frog in
609  time giving second order accuracy and more stability.  time giving second order accuracy and more stability.
610    
611  The essential change in the staggered algorithm is that the  The essential change in the staggered algorithm is that the
612  thermodynamics solver is delayed from half a time step,  thermodynamics solver is delayed from half a time step,
613  allowing the use of the most recent velocities to compute  allowing the use of the most recent velocities to compute
614  the advection terms. Once the thermodynamics fields are  the advection terms. Once the thermodynamics fields are
615  updated, the hydrostatic pressure is computed  updated, the hydrostatic pressure is computed
616  to step forwrad the dynamics.  to step forward the dynamics.
617  Note that the pressure gradient must also be taken out of the  Note that the pressure gradient must also be taken out of the
618  Adams-Bashforth extrapolation. Also, retaining the integer time-levels,  Adams-Bashforth extrapolation. Also, retaining the integer time-levels,
619  $n$ and $n+1$, does not give a user the sense of where variables are  $n$ and $n+1$, does not give a user the sense of where variables are
# Line 631  position in time of variables appropriat Line 631  position in time of variables appropriat
631  \label{eq:vstar-staggered} \\  \label{eq:vstar-staggered} \\
632  \vec{\bf v}^{**} & = & {\cal L}_{\vec{\bf v}}^{-1} ( \vec{\bf v}^* )  \vec{\bf v}^{**} & = & {\cal L}_{\vec{\bf v}}^{-1} ( \vec{\bf v}^* )
633  \label{eq:vstarstar-staggered} \\  \label{eq:vstarstar-staggered} \\
634  \eta^* & = & \epsilon_{fs} \left( \eta^{n-1/2} + \Delta t (P-E)^n \right)- \Delta t  \eta^* & = & \epsilon_{fs} \left( \eta^{n-1/2} + \Delta t (P-E)^n \right)- \Delta t
635    \nabla \cdot H \widehat{ \vec{\bf v}^{**} }    \nabla \cdot H \widehat{ \vec{\bf v}^{**} }
636  \label{eq:nstar-staggered} \\  \label{eq:nstar-staggered} \\
637  \nabla \cdot g H \nabla \eta^{n+1/2} & - & \frac{\epsilon_{fs} \eta^{n+1/2}}{\Delta t^2}  \nabla \cdot g H \nabla \eta^{n+1/2} & - & \frac{\epsilon_{fs} \eta^{n+1/2}}{\Delta t^2}
638  ~ = ~ - \frac{\eta^*}{\Delta t^2}  ~ = ~ - \frac{\eta^*}{\Delta t^2}
639  \label{eq:elliptic-staggered} \\  \label{eq:elliptic-staggered} \\
640  \vec{\bf v}^{n+1/2} & = & \vec{\bf v}^{*} - \Delta t g \nabla \eta^{n+1/2}  \vec{\bf v}^{n+1/2} & = & \vec{\bf v}^{**} - \Delta t g \nabla \eta^{n+1/2}
641  \label{eq:v-n+1-staggered} \\  \label{eq:v-n+1-staggered} \\
642  G_{\theta,S}^{n} & = & G_{\theta,S} ( u^{n+1/2}, \theta^{n}, S^{n} )  G_{\theta,S}^{n} & = & G_{\theta,S} ( u^{n+1/2}, \theta^{n}, S^{n} )
643  \label{eq:Gt-n-staggered} \\  \label{eq:Gt-n-staggered} \\
# Line 650  G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil Line 650  G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil
650  \end{eqnarray}  \end{eqnarray}
651  The corresponding calling tree is given in  The corresponding calling tree is given in
652  \ref{fig:call-tree-adams-bashforth-staggered}.  \ref{fig:call-tree-adams-bashforth-staggered}.
653  The staggered algorithm is activated with the run-time flag  The staggered algorithm is activated with the run-time flag
654  {\bf staggerTimeStep}{\em=.TRUE.} in parameter file {\em data},  {\bf staggerTimeStep}{\em=.TRUE.} in parameter file {\em data},
655  namelist {\em PARM01}.  namelist {\em PARM01}.
656    
657  \begin{figure}  \begin{figure}
658  \begin{center} \fbox{ \begin{minipage}{4.7in} \begin{tabbing}  \begin{center} \fbox{ \begin{minipage}{4.7in} \begin{tabbing}
659  aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill  aaa \= aaa \= aaa \= aaa \= aaa \= aaa \kill
660  FORWARD\_STEP \\  FORWARD\_STEP \\
661  \>\> EXTERNAL\_FIELDS\_LOAD\\  \>\> EXTERNAL\_FIELDS\_LOAD\\
662  \>\> DO\_ATMOSPHERIC\_PHYS \\  \>\> DO\_ATMOSPHERIC\_PHYS \\
663  \>\> DO\_OCEANIC\_PHYS \\  \>\> DO\_OCEANIC\_PHYS \\
664  \> DYNAMICS \\  \> DYNAMICS \\
665  \>\> CALC\_PHI\_HYD \` $\phi_{hyd}^n$ (\ref{eq:phi-hyd-staggered}) \\  \>\> CALC\_PHI\_HYD \` $\phi_{hyd}^n$ (\ref{eq:phi-hyd-staggered}) \\
666  \>\> MOM\_FLUXFORM or MOM\_VECINV \` $G_{\vec{\bf v}}^{n-1/2}$  \>\> MOM\_FLUXFORM or MOM\_VECINV \` $G_{\vec{\bf v}}^{n-1/2}$
667      (\ref{eq:Gv-n-staggered})\\      (\ref{eq:Gv-n-staggered})\\
668  \>\> TIMESTEP \` $\vec{\bf v}^*$ (\ref{eq:Gv-n+5-staggered},  \>\> TIMESTEP \` $\vec{\bf v}^*$ (\ref{eq:Gv-n+5-staggered},
669                                    \ref{eq:vstar-staggered}) \\                                    \ref{eq:vstar-staggered}) \\
670  \>\> IMPLDIFF \` $\vec{\bf v}^{**}$ (\ref{eq:vstarstar-staggered}) \\  \>\> IMPLDIFF \` $\vec{\bf v}^{**}$ (\ref{eq:vstarstar-staggered}) \\
671  \> UPDATE\_R\_STAR or UPDATE\_SURF\_DR \` (NonLin-FS only)\\  \> UPDATE\_R\_STAR or UPDATE\_SURF\_DR \` (NonLin-FS only)\\
# Line 677  FORWARD\_STEP \\ Line 677  FORWARD\_STEP \\
677  \>\> CORRECTION\_STEP \` $u^{n+1/2}$,$v^{n+1/2}$ (\ref{eq:v-n+1-staggered})\\  \>\> CORRECTION\_STEP \` $u^{n+1/2}$,$v^{n+1/2}$ (\ref{eq:v-n+1-staggered})\\
678  \> THERMODYNAMICS \\  \> THERMODYNAMICS \\
679  \>\> CALC\_GT \\  \>\> CALC\_GT \\
680  \>\>\> GAD\_CALC\_RHS \` $G_\theta^n = G_\theta( u, \theta^n )$  \>\>\> GAD\_CALC\_RHS \` $G_\theta^n = G_\theta( u, \theta^n )$
681       (\ref{eq:Gt-n-staggered})\\       (\ref{eq:Gt-n-staggered})\\
682  \>\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\  \>\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\
683  \>\>\> ADAMS\_BASHFORTH2 \` $G_\theta^{(n+1/2)}$ (\ref{eq:Gt-n+5-staggered}) \\  \>\>\> ADAMS\_BASHFORTH2 \` $G_\theta^{(n+1/2)}$ (\ref{eq:Gt-n+5-staggered}) \\
# Line 690  FORWARD\_STEP \\ Line 690  FORWARD\_STEP \\
690  \end{tabbing} \end{minipage} } \end{center}  \end{tabbing} \end{minipage} } \end{center}
691  \caption{  \caption{
692  Calling tree for the overall staggered algorithm using  Calling tree for the overall staggered algorithm using
693  Adams-Bashforth time-stepping.  Adams-Bashforth time-stepping.
694  The place where the model geometry  The place where the model geometry
695  ({\bf hFac} factors) is updated is added here but is only relevant  ({\bf hFac} factors) is updated is added here but is only relevant
696  for the non-linear free-surface algorithm.  for the non-linear free-surface algorithm.
697  }  }
698  \label{fig:call-tree-adams-bashforth-staggered}  \label{fig:call-tree-adams-bashforth-staggered}
# Line 718  time-level variables and terms correspon Line 718  time-level variables and terms correspon
718  The non-hydrostatic formulation re-introduces the full vertical  The non-hydrostatic formulation re-introduces the full vertical
719  momentum equation and requires the solution of a 3-D elliptic  momentum equation and requires the solution of a 3-D elliptic
720  equations for non-hydrostatic pressure perturbation. We still  equations for non-hydrostatic pressure perturbation. We still
721  intergrate vertically for the hydrostatic pressure and solve a 2-D  integrate vertically for the hydrostatic pressure and solve a 2-D
722  elliptic equation for the surface pressure/elevation for this reduces  elliptic equation for the surface pressure/elevation for this reduces
723  the amount of work needed to solve for the non-hydrostatic pressure.  the amount of work needed to solve for the non-hydrostatic pressure.
724    
# Line 759  Substituting into the depth integrated c Line 759  Substituting into the depth integrated c
759  \partial_x H \partial_x \left( g \eta^{n+1} + \widehat{\phi}_{nh}^{n+1} \right)  \partial_x H \partial_x \left( g \eta^{n+1} + \widehat{\phi}_{nh}^{n+1} \right)
760  +  +
761  \partial_y H \partial_y \left( g \eta^{n+1} + \widehat{\phi}_{nh}^{n+1} \right)  \partial_y H \partial_y \left( g \eta^{n+1} + \widehat{\phi}_{nh}^{n+1} \right)
762   - \frac{\epsilon_{fs}\eta^*}{\Delta t^2}   - \frac{\epsilon_{fs}\eta^{n+1}}{\Delta t^2}
763  = - \frac{\eta^*}{\Delta t^2}  = - \frac{\eta^*}{\Delta t^2}
764  \end{equation}  \end{equation}
765  which is approximated by equation  which is approximated by equation
# Line 767  which is approximated by equation Line 767  which is approximated by equation
767  $\phi_{nh}^{n+1}$ is not yet known and ii) $\nabla \widehat{\phi}_{nh}  $\phi_{nh}^{n+1}$ is not yet known and ii) $\nabla \widehat{\phi}_{nh}
768  << g \nabla \eta$. If \ref{eq:elliptic-backward-free-surface} is  << g \nabla \eta$. If \ref{eq:elliptic-backward-free-surface} is
769  solved accurately then the implication is that $\widehat{\phi}_{nh}  solved accurately then the implication is that $\widehat{\phi}_{nh}
770  \approx 0$ so that thet non-hydrostatic pressure field does not drive  \approx 0$ so that the non-hydrostatic pressure field does not drive
771  barotropic motion.  barotropic motion.
772    
773  The flow must satisfy non-divergence  The flow must satisfy non-divergence
# Line 775  The flow must satisfy non-divergence Line 775  The flow must satisfy non-divergence
775  integrated, and this constraint is used to form a 3-D elliptic  integrated, and this constraint is used to form a 3-D elliptic
776  equations for $\phi_{nh}^{n+1}$:  equations for $\phi_{nh}^{n+1}$:
777  \begin{equation}  \begin{equation}
778  \partial_{xx} \phi_{nh}^{n+1} + \partial_{yy} \phi_{nh}^{n+1} +  \partial_{xx} \phi_{nh}^{n+1} + \partial_{yy} \phi_{nh}^{n+1} +
779  \partial_{rr} \phi_{nh}^{n+1} =  \partial_{rr} \phi_{nh}^{n+1} = \left(
780  \partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*}  \partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*}
781    \right) / \Delta t
782  \end{equation}  \end{equation}
783    
784  The entire algorithm can be summarized as the sequential solution of  The entire algorithm can be summarized as the sequential solution of
# Line 787  u^{*} & = & u^{n} + \Delta t G_u^{(n+1/2 Line 788  u^{*} & = & u^{n} + \Delta t G_u^{(n+1/2
788  v^{*} & = & v^{n} + \Delta t G_v^{(n+1/2)} \label{eq:vstar-nh} \\  v^{*} & = & v^{n} + \Delta t G_v^{(n+1/2)} \label{eq:vstar-nh} \\
789  w^{*} & = & w^{n} + \Delta t G_w^{(n+1/2)} \label{eq:wstar-nh} \\  w^{*} & = & w^{n} + \Delta t G_w^{(n+1/2)} \label{eq:wstar-nh} \\
790  \eta^* ~ = ~ \epsilon_{fs} \left( \eta^{n} + \Delta t (P-E) \right)  \eta^* ~ = ~ \epsilon_{fs} \left( \eta^{n} + \Delta t (P-E) \right)
791  & - & \Delta t  & - & \Delta t \left( \partial_x H \widehat{u^{*}}
792    \partial_x H \widehat{u^{*}}                      + \partial_y H \widehat{v^{*}} \right)
 + \partial_y H \widehat{v^{*}}  
793  \\  \\
794    \partial_x g H \partial_x \eta^{n+1}    \partial_x g H \partial_x \eta^{n+1}
795  + \partial_y g H \partial_y \eta^{n+1}  + \partial_y g H \partial_y \eta^{n+1}
# Line 800  w^{*} & = & w^{n} + \Delta t G_w^{(n+1/2 Line 800  w^{*} & = & w^{n} + \Delta t G_w^{(n+1/2
800  \\  \\
801  u^{**} & = & u^{*} - \Delta t g \partial_x \eta^{n+1} \label{eq:unx-nh}\\  u^{**} & = & u^{*} - \Delta t g \partial_x \eta^{n+1} \label{eq:unx-nh}\\
802  v^{**} & = & v^{*} - \Delta t g \partial_y \eta^{n+1} \label{eq:vnx-nh}\\  v^{**} & = & v^{*} - \Delta t g \partial_y \eta^{n+1} \label{eq:vnx-nh}\\
803  \partial_{xx} \phi_{nh}^{n+1} + \partial_{yy} \phi_{nh}^{n+1} +  \partial_{xx} \phi_{nh}^{n+1} + \partial_{yy} \phi_{nh}^{n+1} +
804  \partial_{rr} \phi_{nh}^{n+1} & = &  \partial_{rr} \phi_{nh}^{n+1} & = & \left(
805  \partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*}  \label{eq:phi-nh}\\  \partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*}
806    \right) / \Delta t  \label{eq:phi-nh}\\
807  u^{n+1} & = & u^{**} - \Delta t \partial_x \phi_{nh}^{n+1} \label{eq:un+1-nh}\\  u^{n+1} & = & u^{**} - \Delta t \partial_x \phi_{nh}^{n+1} \label{eq:un+1-nh}\\
808  v^{n+1} & = & v^{**} - \Delta t \partial_y \phi_{nh}^{n+1} \label{eq:vn+1-nh}\\  v^{n+1} & = & v^{**} - \Delta t \partial_y \phi_{nh}^{n+1} \label{eq:vn+1-nh}\\
809  \partial_r w^{n+1} & = & - \partial_x u^{n+1} - \partial_y v^{n+1}  \partial_r w^{n+1} & = & - \partial_x u^{n+1} - \partial_y v^{n+1}
810  \end{eqnarray}  \end{eqnarray}
811  where the last equation is solved by vertically integrating for  where the last equation is solved by vertically integrating for
812  $w^{n+1}$.  $w^{n+1}$.
   
   
813    
814  \section{Variants on the Free Surface}  \section{Variants on the Free Surface}
815  \label{sec:free-surface}  \label{sec:free-surface}
# Line 833  where Line 832  where
832  \begin{eqnarray}  \begin{eqnarray}
833  {\eta}^* = \epsilon_{fs} \: {\eta}^{n} -  {\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
834  \Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o} \vec{\bf v}^* dr  \Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o} \vec{\bf v}^* dr
835  \: + \: \epsilon_{fw} \Delta t (P-E)^{n}  \: + \: \epsilon_{fw} \Delta t (P-E)^{n}
836  \label{eq-solve2D_rhs}  \label{eq-solve2D_rhs}
837  \end{eqnarray}  \end{eqnarray}
838    
# Line 852  $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h) Line 851  $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h)
851    
852    
853  Once ${\eta}^{n+1}$ has been found, substituting into  Once ${\eta}^{n+1}$ has been found, substituting into
854  \ref{eq:discrete-time-u}, \ref{eq:discrete-time-v} yields $\vec{\bf v}^{n+1}$  \ref{eq:discrete-time-u}, \ref{eq:discrete-time-v} yields $\vec{\bf v}^{n+1}$
855  if the model is hydrostatic ($\epsilon_{nh}=0$):  if the model is hydrostatic ($\epsilon_{nh}=0$):
856  $$  $$
857  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}
# Line 875  where Line 874  where
874  \end{displaymath}  \end{displaymath}
875  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
876  $\partial_r \dot{r}^* $ since  $\partial_r \dot{r}^* $ since
877  the vertical velocity at the surface ($\dot{r}_{surf}$)  the vertical velocity at the surface ($\dot{r}_{surf}$)
878  is evaluated as $(\eta^{n+1} - \eta^n) / \Delta t$.  is evaluated as $(\eta^{n+1} - \eta^n) / \Delta t$.
879    
880  Finally, the horizontal velocities at the new time level are found by:  Finally, the horizontal velocities at the new time level are found by:
# Line 925  at the same point in the code. Line 924  at the same point in the code.
924    
925    
926    
927  \subsection{Crank-Nickelson barotropic time stepping}  \subsection{Crank-Nicolson barotropic time stepping}
928  \label{sec:freesurf-CrankNick}  \label{sec:freesurf-CrankNick}
929    
930  The full implicit time stepping described previously is  The full implicit time stepping described previously is
# Line 937  for the barotropic flow divergence ($\ga Line 936  for the barotropic flow divergence ($\ga
936  \\  \\
937  For instance, $\beta=\gamma=1$ is the previous fully implicit scheme;  For instance, $\beta=\gamma=1$ is the previous fully implicit scheme;
938  $\beta=\gamma=1/2$ is the non damping (energy conserving), unconditionally  $\beta=\gamma=1/2$ is the non damping (energy conserving), unconditionally
939  stable, Crank-Nickelson scheme; $(\beta,\gamma)=(1,0)$ or $=(0,1)$  stable, Crank-Nicolson scheme; $(\beta,\gamma)=(1,0)$ or $=(0,1)$
940  corresponds to the forward - backward scheme that conserves energy but is  corresponds to the forward - backward scheme that conserves energy but is
941  only stable for small time steps.\\  only stable for small time steps.\\
942  In the code, $\beta,\gamma$ are defined as parameters, respectively  In the code, $\beta,\gamma$ are defined as parameters, respectively
943  {\bf implicSurfPress}, {\bf implicDiv2DFlow}. They are read from  {\bf implicSurfPress}, {\bf implicDiv2DFlow}. They are read from
944  the main parameter file "{\em data}" and are set by default to 1,1.  the main parameter file "{\em data}" (namelist {\em PARM01})
945    and are set by default to 1,1.
946    
947  Equations \ref{eq:ustar-backward-free-surface} --  Equations \ref{eq:ustar-backward-free-surface} --
948  \ref{eq:vn+1-backward-free-surface} are modified as follows:  \ref{eq:vn+1-backward-free-surface} are modified as follows:
949  \begin{eqnarray*}  \begin{eqnarray*}
950  \frac{ \vec{\bf v}^{n+1} }{ \Delta t }  \frac{ \vec{\bf v}^{n+1} }{ \Delta t }
951  + {\bf \nabla}_h b_s [ \beta {\eta}^{n+1} + (1-\beta) {\eta}^{n} ]  + {\bf \nabla}_h b_s [ \beta {\eta}^{n+1} + (1-\beta) {\eta}^{n} ]
952  + \epsilon_{nh} {\bf \nabla}_h {\phi'_{nh}}^{n+1}  + \epsilon_{nh} {\bf \nabla}_h {\phi'_{nh}}^{n+1}
953   = \frac{ \vec{\bf v}^* }{ \Delta t }   = \frac{ \vec{\bf v}^{n} }{ \Delta t }
954     + \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}
955     + {\bf \nabla}_h {\phi'_{hyd}}^{(n+1/2)}
956  \end{eqnarray*}  \end{eqnarray*}
957  \begin{eqnarray}  \begin{eqnarray}
958  \epsilon_{fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t}  \epsilon_{fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t}
959  + {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o}  + {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o}
960  [ \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
961  = \epsilon_{fw} (P-E)  = \epsilon_{fw} (P-E)
962  \label{eq:eta-n+1-CrankNick}  \label{eq:eta-n+1-CrankNick}
963  \end{eqnarray}  \end{eqnarray}
964  where:  We set
965  \begin{eqnarray*}  \begin{eqnarray*}
966  \vec{\bf v}^* & = &  \vec{\bf v}^* & = &
967  \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)}
# Line 967  where: Line 969  where:
969  + \Delta t {\bf \nabla}_h {\phi'_{hyd}}^{(n+1/2)}  + \Delta t {\bf \nabla}_h {\phi'_{hyd}}^{(n+1/2)}
970  \\  \\
971  {\eta}^* & = &  {\eta}^* & = &
972  \epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta t (P-E)  \epsilon_{fs} {\eta}^{n} + \epsilon_{fw} \Delta t (P-E)
973  - \Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o}  - \Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o}
974  [ \gamma \vec{\bf v}^* + (1-\gamma) \vec{\bf v}^{n}] dr  [ \gamma \vec{\bf v}^* + (1-\gamma) \vec{\bf v}^{n}] dr
975  \end{eqnarray*}  \end{eqnarray*}
976  \\  \\
# Line 979  $$ Line 981  $$
981  {\bf \nabla}_h \cdot \beta\gamma \Delta t^2 b_s (R_o - R_{fixed})  {\bf \nabla}_h \cdot \beta\gamma \Delta t^2 b_s (R_o - R_{fixed})
982  {\bf \nabla}_h {\eta}^{n+1}  {\bf \nabla}_h {\eta}^{n+1}
983  = {\eta}^*  = {\eta}^*
984  $$  $$
985  and then to compute ({\em CORRECTION\_STEP}):  and then to compute ({\em CORRECTION\_STEP}):
986  $$  $$
987  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}
988  - \beta \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}  - \beta \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
989  $$  $$
990    
991  %The non-hydrostatic part is solved as described previously.  %The non-hydrostatic part is solved as described previously.
992    
993  \noindent  \noindent
994  Notes:  Notes:
995  \begin{enumerate}  \begin{enumerate}
996  \item The RHS term of equation \ref{eq:eta-n+1-CrankNick}  \item The RHS term of equation \ref{eq:eta-n+1-CrankNick}
997  corresponds the contribution of fresh water flux (P-E)  corresponds the contribution of fresh water flux (P-E)
998  to the free-surface variations ($\epsilon_{fw}=1$,  to the free-surface variations ($\epsilon_{fw}=1$,
999  {\bf useRealFreshWater}{\em=TRUE} in parameter file {\em data}).  {\bf useRealFreshWater}{\em=TRUE} in parameter file {\em data}).
1000  In order to remain consistent with the tracer equation, specially in  In order to remain consistent with the tracer equation, specially in
1001  the non-linear free-surface formulation, this term is also  the non-linear free-surface formulation, this term is also
1002  affected by the Crank-Nickelson time stepping. The RHS reads:  affected by the Crank-Nicolson time stepping. The RHS reads:
1003  $\epsilon_{fw} ( \gamma (P-E)^{n+1/2} + (1-\gamma) (P-E)^{n-1/2} )$  $\epsilon_{fw} ( \gamma (P-E)^{n+1/2} + (1-\gamma) (P-E)^{n-1/2} )$
1004  \item The non-hydrostatic part of the code has not yet been  %\item The non-hydrostatic part of the code has not yet been
1005  updated, and therefore cannot be used with $(\beta,\gamma) \neq (1,1)$.  %updated, and therefore cannot be used with $(\beta,\gamma) \neq (1,1)$.
1006  \item The stability criteria with Crank-Nickelson time stepping  \item The stability criteria with Crank-Nicolson time stepping
1007  for the pure linear gravity wave problem in cartesian coordinates is:  for the pure linear gravity wave problem in cartesian coordinates is:
1008  \begin{itemize}  \begin{itemize}
1009  \item $\beta + \gamma < 1$ : unstable  \item $\beta + \gamma < 1$ : unstable
1010  \item $\beta \geq 1/2$ and $ \gamma \geq 1/2$ : stable  \item $\beta \geq 1/2$ and $ \gamma \geq 1/2$ : stable
1011  \item $\beta + \gamma \geq 1$ : stable if  \item $\beta + \gamma \geq 1$ : stable if
1012  $$  $$
1013  c_{max}^2 (\beta - 1/2)(\gamma - 1/2) + 1 \geq 0  c_{max}^2 (\beta - 1/2)(\gamma - 1/2) + 1 \geq 0
1014  $$  $$
1015  $$  $$
1016  \mbox{with }~  \mbox{with }~
1017  %c^2 = 2 g H {\Delta t}^2  %c^2 = 2 g H {\Delta t}^2
1018  %(\frac{1-cos 2 \pi / k}{\Delta x^2}  %(\frac{1-cos 2 \pi / k}{\Delta x^2}
1019  %+\frac{1-cos 2 \pi / l}{\Delta y^2})  %+\frac{1-cos 2 \pi / l}{\Delta y^2})
1020  %$$  %$$
1021  %Practically, the most stringent condition is obtained with $k=l=2$ :  %Practically, the most stringent condition is obtained with $k=l=2$ :
1022  %$$  %$$
1023  c_{max} =  2 \Delta t \: \sqrt{g H} \:  c_{max} =  2 \Delta t \: \sqrt{g H} \:
1024  \sqrt{ \frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} }  \sqrt{ \frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} }
1025  $$  $$
1026  \end{itemize}  \end{itemize}
1027    \item A similar mixed forward/backward time-stepping is also available for
1028    the non-hydrostatic algorithm,
1029    with a fraction $\beta_{nh}$ ($ 0 < \beta_{nh} \leq 1$)
1030    of the non-hydrostatic pressure gradient being evaluated at time step $n+1$
1031    (backward in time) and the remaining part ($1 - \beta_{nh}$) being evaluated
1032    at time step $n$ (forward in time).
1033    The run-time parameter {\bf implicitNHPress} corresponding to the implicit
1034    fraction $\beta_{nh}$ of the non-hydrostatic pressure is set by default to
1035    the implicit fraction $\beta$ of surface pressure ({\bf implicSurfPress}),
1036    but can also be specified independently (in main parameter file {\em data},
1037    namelist {\em PARM01}).
1038  \end{enumerate}  \end{enumerate}
   

Legend:
Removed from v.1.28  
changed lines
  Added in v.1.34

  ViewVC Help
Powered by ViewVC 1.1.22