/[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.9 by jmc, Sun Oct 17 04:14:21 2004 UTC revision 1.15 by jmc, Tue Feb 26 21:24:12 2013 UTC
# Line 4  Line 4 
4    
5    
6  \subsection{Non-linear free-surface}  \subsection{Non-linear free-surface}
7  \label{sect:nonlinear-freesurface}  \label{sec: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{pressure/geo-potential and free surface}  \subsubsection{pressure/geo-potential and free surface}
14    \label{sec:phi-freesurface}
15    
16  For the atmosphere, since $\phi = \phi_{topo} - \int^p_{p_s} \alpha dp$,  For the atmosphere, since $\phi = \phi_{topo} - \int^p_{p_s} \alpha dp$,
17  subtracting the reference state defined in section  subtracting the reference state defined in section
18  \ref{sec:hpe-p-geo-potential-split}~:\\  \ref{sec:hpe-p-geo-potential-split}~:\\
19  $$  $$
20  \phi_o = \phi_{topo} - \int^p_{p_o} \alpha_o dp    \phi_o = \phi_{topo} - \int^p_{p_o} \alpha_o dp
21  \hspace{5mm}\mathrm{with}\hspace{3mm} \phi_o(p_o)=\phi_{topo}  \hspace{5mm}\mathrm{with}\hspace{3mm} \phi_o(p_o)=\phi_{topo}
22  $$  $$
23  we get:  we get:
# Line 27  For the ocean, the reference state is si Line 28  For the ocean, the reference state is si
28  on $z$ ($b_o=g$) and the surface reference position is uniformly $z=0$ ($R_o=0$),  on $z$ ($b_o=g$) and the surface reference position is uniformly $z=0$ ($R_o=0$),
29  and the same subtraction leads to a similar relation.  and the same subtraction leads to a similar relation.
30  For both fluid, using the isomorphic notations, we can write:  For both fluid, using the isomorphic notations, we can write:
31  $$  $$
32  \phi' = \int^{r_{surf}}_r b~ dr - \int^{R_o}_r b_o dr  \phi' = \int^{r_{surf}}_r b~ dr - \int^{R_o}_r b_o dr
33  $$  $$
34  \begin{eqnarray}  and re-write as:
35  \mathrm{and~re~write:}\hspace{10mm}  \begin{equation}
36  \phi' = \int^{r_{surf}}_{R_o} b~ dr & + & \int^{R_o}_r (b - b_o) dr  \phi' = \int^{r_{surf}}_{R_o} b~ dr + \int^{R_o}_r (b - b_o) dr
37  \label{eq:split-phi-Ro} \\  \label{eq:split-phi-Ro}
38  \mathrm{or:}\hspace{10mm}  \end{equation}
39  \phi' = \int^{r_{surf}}_{R_o} b_o dr & + & \int^{r_{surf}}_r (b - b_o) dr  or:
40    \begin{equation}
41    \phi' = \int^{r_{surf}}_{R_o} b_o dr + \int^{r_{surf}}_r (b - b_o) dr
42  \label{eq:split-phi-bo}  \label{eq:split-phi-bo}
43  \end{eqnarray}  \end{equation}
44    
45  In section \ref{sec:finding_the_pressure_field}, following eq.\ref{eq:split-phi-Ro},  In section \ref{sec:finding_the_pressure_field}, following eq.\ref{eq:split-phi-Ro},
46  the pressure/geo-potential $\phi'$ has been separated into surface ($\phi_s$),  the pressure/geo-potential $\phi'$ has been separated into surface ($\phi_s$),
47  and hydrostatic anomaly ($\phi'_{hyd}$).  and hydrostatic anomaly ($\phi'_{hyd}$).
48  In this section, the split between $\phi_s$ and $\phi'_{hyd}$ is  In this section, the split between $\phi_s$ and $\phi'_{hyd}$ is
49  made according to equation \ref{eq:split-phi-bo}. This slightly  made according to equation \ref{eq:split-phi-bo}. This slightly
50  different definition reflects the actual implementation in the code  different definition reflects the actual implementation in the code
51  and is valid for both linear and non-linear  and is valid for both linear and non-linear
52  free-surface formulation, in both r-coordinate and r*-coordinate.  free-surface formulation, in both r-coordinate and r*-coordinate.
53    
54  Because the linear free-surface approximation ignore the tracer content  Because the linear free-surface approximation ignore the tracer content
55  of the fluid parcel between $R_o$ and $r_{surf}=R_o+\eta$,  of the fluid parcel between $R_o$ and $r_{surf}=R_o+\eta$,
56  for consistency reasons, this part is also neglected in $\phi'_{hyd}$~:  for consistency reasons, this part is also neglected in $\phi'_{hyd}$~:
57  $$  $$
58  \phi'_{hyd} = \int^{r_{surf}}_r (b - b_o) dr \simeq \int^{R_o}_r (b - b_o) dr  \phi'_{hyd} = \int^{r_{surf}}_r (b - b_o) dr \simeq \int^{R_o}_r (b - b_o) dr
59  $$  $$
60  Note that in this case, the two definitions of $\phi_s$ and $\phi'_{hyd}$  Note that in this case, the two definitions of $\phi_s$ and $\phi'_{hyd}$
61  from equation \ref{eq:split-phi-Ro} and \ref{eq:split-phi-bo} converge toward  from equation \ref{eq:split-phi-Ro} and \ref{eq:split-phi-bo} converge toward
62  the same (approximated) expressions: $\phi_s = \int^{r_{surf}}_{R_o} b_o dr$  the same (approximated) expressions: $\phi_s = \int^{r_{surf}}_{R_o} b_o dr$
63  and $\phi'_{hyd}=\int^{R_o}_r b' dr$.\\  and $\phi'_{hyd}=\int^{R_o}_r b' dr$.\\
64  On the contrary, the unapproximated formulation ("non-linear free-surface",  On the contrary, the unapproximated formulation ("non-linear free-surface",
65  see the next section) retains the full expression:  see the next section) retains the full expression:
# Line 64  $\phi'_{hyd} = \int^{r_{surf}}_r (b - b_ Line 67  $\phi'_{hyd} = \int^{r_{surf}}_r (b - b_
67  This is obtained by selecting {\bf nonlinFreeSurf}=4 in parameter  This is obtained by selecting {\bf nonlinFreeSurf}=4 in parameter
68  file {\em data}.\\  file {\em data}.\\
69    
70  Regarding the surface potential:  Regarding the surface potential:
71  $$\phi_s = \int_{R_o}^{R_o+\eta} b_o dr = b_s \eta  $$\phi_s = \int_{R_o}^{R_o+\eta} b_o dr = b_s \eta
72  \hspace{5mm}\mathrm{with}\hspace{5mm}  \hspace{5mm}\mathrm{with}\hspace{5mm}
73  b_s = \frac{1}{\eta} \int_{R_o}^{R_o+\eta} b_o dr $$  b_s = \frac{1}{\eta} \int_{R_o}^{R_o+\eta} b_o dr $$
74  $b_s \simeq b_o(R_o)$ is an excellent approximation (better than  $b_s \simeq b_o(R_o)$ is an excellent approximation (better than
75  the usual numerical truncation, since generally $|\eta|$ is smaller  the usual numerical truncation, since generally $|\eta|$ is smaller
76  than the vertical grid increment).  than the vertical grid increment).
77    
78  For the ocean, $\phi_s = g \eta$ and $b_s = g$ is uniform.  For the ocean, $\phi_s = g \eta$ and $b_s = g$ is uniform.
79  For the atmosphere, however, because of topographic effects, the  For the atmosphere, however, because of topographic effects, the
80  reference surface pressure $R_o=p_o$ has large spatial variations that  reference surface pressure $R_o=p_o$ has large spatial variations that
81  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
82  $[m^3/kg]$). For this reason, when {\bf uniformLin\_PhiSurf} {\em=.FALSE.}  $[m^3/kg]$). For this reason, when {\bf uniformLin\_PhiSurf} {\em=.FALSE.}
83  (parameter file {\em data}, namelist {\em PARAM01})  (parameter file {\em data}, namelist {\em PARAM01})
84  a non-uniform linear coefficient $b_s$ is used and computed  a non-uniform linear coefficient $b_s$ is used and computed
85  ({\it S/R INI\_LINEAR\_PHISURF}) according to the reference surface  ({\it S/R INI\_LINEAR\_PHISURF}) according to the reference surface
86  pressure $p_o$:  pressure $p_o$:
87  $b_s = b_o(R_o) = c_p \kappa (p_o / P^o_{SL})^{(\kappa - 1)} \theta_{ref}(p_o)$.  $b_s = b_o(R_o) = c_p \kappa (p_o / P^o_{SL})^{(\kappa - 1)} \theta_{ref}(p_o)$.
88  with $P^o_{SL}$ the mean sea-level pressure.  with $P^o_{SL}$ the mean sea-level pressure.
# Line 89  with $P^o_{SL}$ the mean sea-level press Line 92  with $P^o_{SL}$ the mean sea-level press
92  (Non-linear free-surface)}  (Non-linear free-surface)}
93    
94  The total thickness of the fluid column is $r_{surf} - R_{fixed} =  The total thickness of the fluid column is $r_{surf} - R_{fixed} =
95  \eta + R_o - R_{fixed}$. In most applications, the free surface  \eta + R_o - R_{fixed}$. In most applications, the free surface
96  displacements are small compared to the total thickness  displacements are small compared to the total thickness
97  $\eta \ll H_o = R_o - R_{fixed}$.  $\eta \ll H_o = R_o - R_{fixed}$.
98  In the previous sections and in older version of the model,  In the previous sections and in older version of the model,
99  the linearized free-surface approximation was made, assuming  the linearized free-surface approximation was made, assuming
100  $r_{surf} - R_{fixed} \simeq H_o$ when computing horizontal transports,  $r_{surf} - R_{fixed} \simeq H_o$ when computing horizontal transports,
101  either in the continuity equation or in tracer and momentum  either in the continuity equation or in tracer and momentum
102  advection terms.  advection terms.
103  This approximation is dropped when using the non-linear free-surface  This approximation is dropped when using the non-linear free-surface
104  formulation and the total thickness, including the time varying part  formulation and the total thickness, including the time varying part
105  $\eta$, is considered when computing horizontal transports.  $\eta$, is considered when computing horizontal transports.
106  Implications for the barotropic part are presented hereafter.  Implications for the barotropic part are presented hereafter.
107  In section \ref{sect:freesurf-tracer-advection} consequences for  In section \ref{sec:freesurf-tracer-advection} consequences for
108  tracer conservation is briefly discussed (more details can be  tracer conservation is briefly discussed (more details can be
109  found in \cite{campin:02})~; the general time-stepping is presented  found in \cite{campin:02})~; the general time-stepping is presented
110  in section \ref{sect:nonlin-freesurf-timestepping} with some  in section \ref{sec:nonlin-freesurf-timestepping} with some
111  limitations regarding the vertical resolution in section  limitations regarding the vertical resolution in section
112  \ref{sect:nonlin-freesurf-dz_surf}.  \ref{sec:nonlin-freesurf-dz_surf}.
113    
114  In the non-linear formulation, the continuous form of the model  In the non-linear formulation, the continuous form of the model
115  equations remains unchanged, except for the 2D continuity equation  equations remains unchanged, except for the 2D continuity equation
116  (\ref{eq:discrete-time-backward-free-surface}) which is now  (\ref{eq:discrete-time-backward-free-surface}) which is now
117  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$ :
118    
# Line 139  where Line 142  where
142  {\eta}^* = \epsilon_{fs} \: {\eta}^{n} -  {\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
143  \Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o+\eta^n} \vec{\bf v}^* dr  \Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o+\eta^n} \vec{\bf v}^* dr
144  \: + \: \epsilon_{fw} \Delta_t (P-E)^{n}  \: + \: \epsilon_{fw} \Delta_t (P-E)^{n}
145  \end{eqnarray*}  \end{eqnarray*}
146  This method requires us to update the solver matrix at each time step.  This method requires us to update the solver matrix at each time step.
147    
148  Alternatively, the non-linear contribution can be evaluated fully  Alternatively, the non-linear contribution can be evaluated fully
# Line 151  explicitly: Line 154  explicitly:
154  = {\eta}^*  = {\eta}^*
155  +{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n})  +{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n})
156  {\bf \nabla}_h b_s {\eta}^{n}  {\bf \nabla}_h b_s {\eta}^{n}
157  \end{eqnarray*}  \end{eqnarray*}
158  This formulation allows one to keep the initial solver matrix  This formulation allows one to keep the initial solver matrix
159  unchanged though throughout the integration, since the non-linear free  unchanged though throughout the integration, since the non-linear free
160  surface only affects the RHS.  surface only affects the RHS.
# Line 163  column thickness appears only in the int Line 166  column thickness appears only in the int
166    
167  Those different options (see Table \ref{tab:nonLinFreeSurf_flags})  Those different options (see Table \ref{tab:nonLinFreeSurf_flags})
168  have been tested and show little differences. However, we recommend  have been tested and show little differences. However, we recommend
169  the use of the most precise method (the 1rst one) since the  the use of the most precise method (the 1rst one) since the
170  computation cost involved in the solver matrix update is negligible.  computation cost involved in the solver matrix update is negligible.
171    
172  \begin{table}[htb]  \begin{table}[htb]
# Line 171  computation cost involved in the solver Line 174  computation cost involved in the solver
174  \centering  \centering
175   \begin{tabular}[htb]{|l|c|l|}   \begin{tabular}[htb]{|l|c|l|}
176     \hline     \hline
177     parameter & value & description \\     parameter & value & description \\
178     \hline     \hline
179                     & -1 & linear free-surface, restart from a pickup file \\                     & -1 & linear free-surface, restart from a pickup file \\
180                     &    & produced with \#undef EXACT\_CONSERV code\\                     &    & produced with \#undef EXACT\_CONSERV code\\
181     \cline{2-3}     \cline{2-3}
182                     & 0 & Linear free-surface \\                     & 0 & Linear free-surface \\
183     \cline{2-3}     \cline{2-3}
184      nonlinFreeSurf & 4 & Non-linear free-surface \\      nonlinFreeSurf & 4 & Non-linear free-surface \\
185     \cline{2-3}     \cline{2-3}
# Line 202  computation cost involved in the solver Line 205  computation cost involved in the solver
205    
206    
207  \subsubsection{Tracer conservation with non-linear free-surface}  \subsubsection{Tracer conservation with non-linear free-surface}
208  \label{sect:freesurf-tracer-advection}  \label{sec:freesurf-tracer-advection}
209    
210  To ensure global tracer conservation (i.e., the total amount) as well  To ensure global tracer conservation (i.e., the total amount) as well
211  as local conservation, the change in the surface level thickness must  as local conservation, the change in the surface level thickness must
# Line 222  $$ Line 225  $$
225    = P \theta_{\mathrm{rain}}    = P \theta_{\mathrm{rain}}
226  $$  $$
227  Using the implicit (non-linear) free surface described above (section  Using the implicit (non-linear) free surface described above (section
228  \ref{sect:pressure-method-linear-backward}) we have:  \ref{sec:pressure-method-linear-backward}) we have:
229  \begin{eqnarray*}  \begin{eqnarray*}
230  h^{n+1} = h^{n} - \Delta t \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) + \Delta t P \\  h^{n+1} = h^{n} - \Delta t \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) + \Delta t P \\
231  \end{eqnarray*}  \end{eqnarray*}
# Line 230  The discretized form of the tracer equat Line 233  The discretized form of the tracer equat
233  ``form'' in the computation of tracer fluxes, that is, the same value  ``form'' in the computation of tracer fluxes, that is, the same value
234  of $h$, as used in the continuity equation:  of $h$, as used in the continuity equation:
235  \begin{eqnarray*}  \begin{eqnarray*}
236  h^{n+1} \, \theta^{n+1} = h^n \, \theta^n  h^{n+1} \, \theta^{n+1} = h^n \, \theta^n
237          - \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})
238          + \Delta t P \theta_{rain}          + \Delta t P \theta_{rain}
239  \end{eqnarray*}  \end{eqnarray*}
# Line 247  model uses the $G_\theta$ terms (first s Line 250  model uses the $G_\theta$ terms (first s
250  surface formulation (with the "{\it surface correction}" turned "on",  surface formulation (with the "{\it surface correction}" turned "on",
251  see tracer section):  see tracer section):
252  $$  $$
253  G_\theta^n = \left(- \nabla \cdot (h^n \, \theta^n \, \vec{\bf v}^{n+1})  G_\theta^n = \left(- \nabla \cdot (h^n \, \theta^n \, \vec{\bf v}^{n+1})
254           - \dot{r}_{surf}^{n+1} \theta^n \right) / h^n           - \dot{r}_{surf}^{n+1} \theta^n \right) / h^n
255  $$  $$
256  Then, in a second step, the thickness variation (expansion/reduction)  Then, in a second step, the thickness variation (expansion/reduction)
257  is taken into account:  is taken into account:
258  $$  $$
259  \theta^{n+1} = \theta^n + \Delta t \frac{h^n}{h^{n+1}}  \theta^{n+1} = \theta^n + \Delta t \frac{h^n}{h^{n+1}}
260     \left( G_\theta^{(n+1/2)} + P (\theta_{\mathrm{rain}} - \theta^n )/h^n \right)     \left( G_\theta^{(n+1/2)} + P (\theta_{\mathrm{rain}} - \theta^n )/h^n \right)
261  %\theta^{n+1} = \theta^n + \frac{\Delta t}{h^{n+1}}  %\theta^{n+1} = \theta^n + \frac{\Delta t}{h^{n+1}}
262  %   \left( h^n G_\theta^{(n+1/2)} + P (\theta_{\mathrm{rain}} - \theta^n ) \right)  %   \left( h^n G_\theta^{(n+1/2)} + P (\theta_{\mathrm{rain}} - \theta^n ) \right)
263  $$  $$
264  Note that with a simple forward time step (no Adams-Bashforth),  Note that with a simple forward time step (no Adams-Bashforth),
265  these two formulations are equivalent,    these two formulations are equivalent,
266  since  since
267  $  $
268  (h^{n+1} - h^{n})/ \Delta t =  (h^{n+1} - h^{n})/ \Delta t =
269  P - \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) = P + \dot{r}_{surf}^{n+1}  P - \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) = P + \dot{r}_{surf}^{n+1}
270  $  $
271    
272  \subsubsection{Time stepping implementation of the  \subsubsection{Time stepping implementation of the
273  non-linear free-surface}      non-linear free-surface}
274  \label{sect:nonlin-freesurf-timestepping}  \label{sec:nonlin-freesurf-timestepping}
275    
276  The grid cell thickness was hold constant with the linear  The grid cell thickness was hold constant with the linear
277  free-surface~; with the non-linear free-surface, it is now varying  free-surface~; with the non-linear free-surface, it is now varying
278  in time, at least at the surface level.  in time, at least at the surface level.
279  This implies some modifications of the general algorithm described  This implies some modifications of the general algorithm described
280  earlier in sections \ref{sect:adams-bashforth-sync} and  earlier in sections \ref{sec:adams-bashforth-sync} and
281  \ref{sect:adams-bashforth-staggered}.  \ref{sec:adams-bashforth-staggered}.
282    
283  A simplified version of the staggered in time, non-linear  A simplified version of the staggered in time, non-linear
284  free-surface algorithm is detailed hereafter, and can be compared  free-surface algorithm is detailed hereafter, and can be compared
285  to the equivalent linear free-surface case (eq. \ref{eq:Gv-n-staggered}  to the equivalent linear free-surface case (eq. \ref{eq:Gv-n-staggered}
286  to \ref{eq:t-n+1-staggered}) and can also be easily transposed  to \ref{eq:t-n+1-staggered}) and can also be easily transposed
287  to the synchronous time-stepping case.  to the synchronous time-stepping case.
# Line 291  $h^n$ and $dh^n$ are the column and grid Line 294  $h^n$ and $dh^n$ are the column and grid
294  \begin{eqnarray}  \begin{eqnarray}
295  \phi^{n}_{hyd} & = & \int b(\theta^{n},S^{n},r) dr  \phi^{n}_{hyd} & = & \int b(\theta^{n},S^{n},r) dr
296  \label{eq:phi-hyd-nlfs} \\  \label{eq:phi-hyd-nlfs} \\
297  \vec{\bf G}_{\vec{\bf v}}^{n-1/2}\hspace{-2mm} & = &  \vec{\bf G}_{\vec{\bf v}}^{n-1/2}\hspace{-2mm} & = &
298  \vec{\bf G}_{\vec{\bf v}} (dh^{n-1},\vec{\bf v}^{n-1/2})  \vec{\bf G}_{\vec{\bf v}} (dh^{n-1},\vec{\bf v}^{n-1/2})
299  \hspace{+2mm};\hspace{+2mm}  \hspace{+2mm};\hspace{+2mm}
300  \vec{\bf G}_{\vec{\bf v}}^{(n)} =    \vec{\bf G}_{\vec{\bf v}}^{(n)} =
301     \frac{3}{2} \vec{\bf G}_{\vec{\bf v}}^{n-1/2}     \frac{3}{2} \vec{\bf G}_{\vec{\bf v}}^{n-1/2}
302  -  \frac{1}{2} \vec{\bf G}_{\vec{\bf v}}^{n-3/2}  -  \frac{1}{2} \vec{\bf G}_{\vec{\bf v}}^{n-3/2}
303  \label{eq:Gv-n-nlfs} \\  \label{eq:Gv-n-nlfs} \\
304  %\vec{\bf G}_{\vec{\bf v}}^{(n)} & = &  %\vec{\bf G}_{\vec{\bf v}}^{(n)} & = &
305  %   \frac{3}{2} \vec{\bf G}_{\vec{\bf v}}^{n-1/2}  %   \frac{3}{2} \vec{\bf G}_{\vec{\bf v}}^{n-1/2}
306  %-  \frac{1}{2} \vec{\bf G}_{\vec{\bf v}}^{n-3/2}  %-  \frac{1}{2} \vec{\bf G}_{\vec{\bf v}}^{n-3/2}
307  %\label{eq:Gv-n+5-nlfs} \\  %\label{eq:Gv-n+5-nlfs} \\
308  %\vec{\bf v}^{*} & = & \vec{\bf v}^{n-1/2} + \frac{\Delta t}{dh^{n}} \left(  %\vec{\bf v}^{*} & = & \vec{\bf v}^{n-1/2} + \frac{\Delta t}{dh^{n}} \left(
# Line 307  $h^n$ and $dh^n$ are the column and grid Line 310  $h^n$ and $dh^n$ are the column and grid
310  \vec{\bf v}^{*} & = & \vec{\bf v}^{n-1/2} + \Delta t \frac{dh^{n-1}}{dh^{n}} \left(  \vec{\bf v}^{*} & = & \vec{\bf v}^{n-1/2} + \Delta t \frac{dh^{n-1}}{dh^{n}} \left(
311  \vec{\bf G}_{\vec{\bf v}}^{(n)} + F_{\vec{\bf v}}^{n}/dh^{n-1} \right)  \vec{\bf G}_{\vec{\bf v}}^{(n)} + F_{\vec{\bf v}}^{n}/dh^{n-1} \right)
312  - \Delta t \nabla \phi_{hyd}^{n}  - \Delta t \nabla \phi_{hyd}^{n}
313  \label{eq:vstar-nlfs} \\  \label{eq:vstar-nlfs}
314    \mathrm{update}\hspace{-4mm} & & \hspace{-4mm}\mathrm{  \end{eqnarray}
315    model~geometry~:~{\bf hFac}}(dh^n)\nonumber \\  \hspace{3cm}$\longrightarrow$~~{\it update model~geometry~:~}${\bf hFac}(dh^n)$\\
316  \eta^{n+1/2} \hspace{-2mm} & = &  \begin{eqnarray}
317    \eta^{n+1/2} \hspace{-2mm} & = &
318  \eta^{n-1/2} + \Delta t P^{n+1/2} - \Delta t  \eta^{n-1/2} + \Delta t P^{n+1/2} - \Delta t
319    \nabla \cdot \int \vec{\bf v}^{n+1/2} dh^{n} \nonumber \\    \nabla \cdot \int \vec{\bf v}^{n+1/2} dh^{n} \nonumber \\
320               & = & \eta^{n-1/2} + \Delta t P^{n+1/2} - \Delta t               & = & \eta^{n-1/2} + \Delta t P^{n+1/2} - \Delta t
321    \nabla \cdot \int \!\!\! \left( \vec{\bf v}^* - g \Delta t \nabla \eta^{n+1/2} \right) dh^{n}    \nabla \cdot \int \!\!\! \left( \vec{\bf v}^* - g \Delta t \nabla \eta^{n+1/2} \right) dh^{n}
322  \label{eq:nstar-nlfs} \\  \label{eq:nstar-nlfs} \\
323  \vec{\bf v}^{n+1/2}\hspace{-2mm} & = &  \vec{\bf v}^{n+1/2}\hspace{-2mm} & = &
324  \vec{\bf v}^{*} - g \Delta t \nabla \eta^{n+1/2}  \vec{\bf v}^{*} - g \Delta t \nabla \eta^{n+1/2}
325  \label{eq:v-n+1-nlfs} \\  \label{eq:v-n+1-nlfs} \\
326  h^{n+1} & = & h^{n} + \Delta t P^{n+1/2} - \Delta t  h^{n+1} & = & h^{n} + \Delta t P^{n+1/2} - \Delta t
# Line 326  G_{\theta}^{n} & = & G_{\theta} ( dh^{n} Line 330  G_{\theta}^{n} & = & G_{\theta} ( dh^{n}
330  \hspace{+2mm};\hspace{+2mm}  \hspace{+2mm};\hspace{+2mm}
331  G_{\theta}^{(n+1/2)} = \frac{3}{2} G_{\theta}^{n} - \frac{1}{2} G_{\theta}^{n-1}  G_{\theta}^{(n+1/2)} = \frac{3}{2} G_{\theta}^{n} - \frac{1}{2} G_{\theta}^{n-1}
332  \label{eq:Gt-n-nlfs} \\  \label{eq:Gt-n-nlfs} \\
333  %\theta^{n+1} & = &\theta^{n} + \frac{\Delta t}{dh^{n+1}} \left( dh^n  %\theta^{n+1} & = &\theta^{n} + \frac{\Delta t}{dh^{n+1}} \left( dh^n
334  %G_{\theta}^{(n+1/2)} + Q^{n+1/2} + P^{n+1/2} (\theta_{\mathrm{rain}}-\theta^n) \right)  %G_{\theta}^{(n+1/2)} + Q^{n+1/2} + P^{n+1/2} (\theta_{\mathrm{rain}}-\theta^n) \right)
335  \theta^{n+1} & = &\theta^{n} + \Delta t \frac{dh^n}{dh^{n+1}} \left(  \theta^{n+1} & = &\theta^{n} + \Delta t \frac{dh^n}{dh^{n+1}} \left(
336  G_{\theta}^{(n+1/2)}  G_{\theta}^{(n+1/2)}
337  +( P^{n+1/2} (\theta_{\mathrm{rain}}-\theta^n) + Q^{n+1/2})/dh^n \right)  +( P^{n+1/2} (\theta_{\mathrm{rain}}-\theta^n) + Q^{n+1/2})/dh^n \right)
338  \nonumber \\  \nonumber \\
339  & & \label{eq:t-n+1-nlfs}  & & \label{eq:t-n+1-nlfs}
340  \end{eqnarray}  \end{eqnarray}
341  %-------------------------------------------------------------  %-------------------------------------------------------------
342  Two steps have been added to linear free-surface algorithm  Two steps have been added to linear free-surface algorithm
343  (eq. \ref{eq:Gv-n-staggered} to \ref{eq:t-n+1-staggered}):  (eq. \ref{eq:Gv-n-staggered} to \ref{eq:t-n+1-staggered}):
344  Firstly, the model ``geometry''  Firstly, the model ``geometry''
345  (here the {\bf hFacC,W,S}) is updated just before entering {\it  (here the {\bf hFacC,W,S}) is updated just before entering {\it
346  SOLVE\_FOR\_PRESSURE}, using the current $dh^{n}$ field.  SOLVE\_FOR\_PRESSURE}, using the current $dh^{n}$ field.
347  Secondly, the vertically integrated continuity equation  Secondly, the vertically integrated continuity equation
348  (eq.\ref{eq:h-n+1-nlfs}) has been added ({\bf exactConserv}{\em =TRUE},  (eq.\ref{eq:h-n+1-nlfs}) has been added ({\bf exactConserv}{\em =TRUE},
349  in parameter file {\em data}, namelist {\em PARM01})  in parameter file {\em data}, namelist {\em PARM01})
350  just before computing the vertical velocity, in subroutine  just before computing the vertical velocity, in subroutine
351  {\em INTEGR\_CONTINUITY}. This ensures that tracer and continuity equation  {\em INTEGR\_CONTINUITY}.
352  discretization  a Although this equation might appear  %This ensures that tracer and continuity equation discretization a
353  redundant with eq.\ref{eq:nstar-nlfs}, the integrated column  Although this equation might appear redundant with eq.\ref{eq:nstar-nlfs},
354  thickness $h^{n+1}$ can be different from $\eta^{n+1/2} + H$~:  the integrated column thickness $h^{n+1}$ will be different from
355    $\eta^{n+1/2} + H$~ in the following cases:
356  \begin{itemize}  \begin{itemize}
357  \item when Crank-Nickelson time-stepping is used (see section  \item when Crank-Nickelson time-stepping is used (see section
358  \ref{sect:freesurf-CrankNick}).  \ref{sec:freesurf-CrankNick}).
359  \item when filters are applied to the flow field, after  \item when filters are applied to the flow field, after
360  (\ref{eq:v-n+1-nlfs}) and alter the divergence of the flow.  (\ref{eq:v-n+1-nlfs}) and alter the divergence of the flow.
361  \item when the solver does not iterate until convergence~;  \item when the solver does not iterate until convergence~;
# Line 362  In this staggered time-stepping algorith Line 367  In this staggered time-stepping algorith
367  are computed using $dh^{n-1}$ geometry factors.  are computed using $dh^{n-1}$ geometry factors.
368  (eq.\ref{eq:Gv-n-nlfs}) and then rescaled in subroutine {\it TIMESTEP},  (eq.\ref{eq:Gv-n-nlfs}) and then rescaled in subroutine {\it TIMESTEP},
369  (eq.\ref{eq:vstar-nlfs}), similarly to tracer tendencies (see section  (eq.\ref{eq:vstar-nlfs}), similarly to tracer tendencies (see section
370  \ref{sect:freesurf-tracer-advection}).  \ref{sec:freesurf-tracer-advection}).
371  The tracers are stepped forward later, using the recently updated  The tracers are stepped forward later, using the recently updated
372  flow field ${\bf v}^{n+1/2}$ and the corresponding model geometry  flow field ${\bf v}^{n+1/2}$ and the corresponding model geometry
373  $dh^{n}$ to compute the tendencies (eq.\ref{eq:Gt-n-nlfs});  $dh^{n}$ to compute the tendencies (eq.\ref{eq:Gt-n-nlfs});
# Line 370  Then the tendencies are rescaled by $dh^ Line 375  Then the tendencies are rescaled by $dh^
375  the new tracers values $(\theta,S)^{n+1}$ (eq.\ref{eq:t-n+1-nlfs},  the new tracers values $(\theta,S)^{n+1}$ (eq.\ref{eq:t-n+1-nlfs},
376  in subroutine {\em CALC\_GT, CALC\_GS}).  in subroutine {\em CALC\_GT, CALC\_GS}).
377    
378  Note that the fresh-water input is added in a consistent way in the  Note that the fresh-water input is added in a consistent way in the
379  continuity equation and in the tracer equation, taking into account  continuity equation and in the tracer equation, taking into account
380  the fresh-water temperature $\theta_{\mathrm{rain}}$.  the fresh-water temperature $\theta_{\mathrm{rain}}$.
381    
382  Regarding the restart procedure,  Regarding the restart procedure,
383  two 2.D fields $h^{n-1}$ and $(h^n-h^{n-1})/\Delta t$  two 2.D fields $h^{n-1}$ and $(h^n-h^{n-1})/\Delta t$
384  in addition to the standard  in addition to the standard
385  state variables and tendencies ($\eta^{n-1/2}$, ${\bf v}^{n-1/2}$,  state variables and tendencies ($\eta^{n-1/2}$, ${\bf v}^{n-1/2}$,
386  $\theta^n$, $S^n$, ${\bf G}_{\bf v}^{n-3/2}$, $G_{\theta,S}^{n-1}$)  $\theta^n$, $S^n$, ${\bf G}_{\bf v}^{n-3/2}$, $G_{\theta,S}^{n-1}$)
# Line 384  The model restarts reading this "{\em pi Line 389  The model restarts reading this "{\em pi
389  then update the model geometry according to $h^{n-1}$,  then update the model geometry according to $h^{n-1}$,
390  and compute $h^n$ and the vertical velocity  and compute $h^n$ and the vertical velocity
391  %$h^n=h^{n-1} + \Delta t [(h^n-h^{n-1})/\Delta t]$,  %$h^n=h^{n-1} + \Delta t [(h^n-h^{n-1})/\Delta t]$,
392  before starting the main calling sequence (eq.\ref{eq:phi-hyd-nlfs}  before starting the main calling sequence (eq.\ref{eq:phi-hyd-nlfs}
393  to \ref{eq:t-n+1-nlfs}, {\em S/R FORWARD\_STEP}).  to \ref{eq:t-n+1-nlfs}, {\em S/R FORWARD\_STEP}).
394  \\  \\
395    
# Line 395  $h^{n+1} -H_o$: {\bf etaH} ({\em DYNVARS Line 400  $h^{n+1} -H_o$: {\bf etaH} ({\em DYNVARS
400    
401  $h^{n} -H_o$: {\bf etaHnm1} ({\em SURFACE.h})  $h^{n} -H_o$: {\bf etaHnm1} ({\em SURFACE.h})
402    
403  $h^{n+1}-h^{n}/\Delta t$: {\bf dEtaHdt} ({\em SURFACE.h})  $(h^{n+1}-h^{n})/\Delta t$: {\bf dEtaHdt} ({\em SURFACE.h})
404    
405  \end{minipage} }  \end{minipage} }
406    
407  \subsubsection{Non-linear free-surface and vertical resolution}  \subsubsection{Non-linear free-surface and vertical resolution}
408  \label{sect:nonlin-freesurf-dz_surf}  \label{sec:nonlin-freesurf-dz_surf}
409    
410  When the amplitude of the free-surface variations becomes  When the amplitude of the free-surface variations becomes
411  as large as the vertical resolution near the surface,  as large as the vertical resolution near the surface,
412  the surface layer thickness can decrease to nearly zero or  the surface layer thickness can decrease to nearly zero or
413  can even vanish completely.  can even vanish completely.
414  This later possibility has not been implemented, and a  This later possibility has not been implemented, and a
415  minimum relative thickness is imposed ({\bf hFacInf},  minimum relative thickness is imposed ({\bf hFacInf},
416  parameter file {\em data}, namelist {\em PARM01}) to prevent  parameter file {\em data}, namelist {\em PARM01}) to prevent
417  numerical instabilities caused by very thin surface level.  numerical instabilities caused by very thin surface level.
418    
419  A better alternative to the vanishing level problem has been  A better alternative to the vanishing level problem has been
420  found and implemented recently, and rely on a different  found and implemented recently, and rely on a different
421  vertical coordinate $r^*$~:  vertical coordinate $r^*$~:
422  The time variation ot the total column thickness becomes  The time variation ot the total column thickness becomes
423  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}$
424  model, but the fixed part related to topography is treated  model, but the fixed part related to topography is treated
425  as in a height or pressure coordinate model.  as in a height or pressure coordinate model.
426  A complete description is given in \cite{adcroft:04}.  A complete description is given in \cite{adcroft:04a}.
427    
428  The time-stepping implementation of the $r^*$ coordinate is  The time-stepping implementation of the $r^*$ coordinate is
429  identical to the non-linear free-surface in $r$ coordinate,  identical to the non-linear free-surface in $r$ coordinate,

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

  ViewVC Help
Powered by ViewVC 1.1.22