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

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

  ViewVC Help
Powered by ViewVC 1.1.22