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

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

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

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

Legend:
Removed from v.1.8  
changed lines
  Added in v.1.12

  ViewVC Help
Powered by ViewVC 1.1.22