/[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.2 by jmc, Fri Aug 17 18:38:10 2001 UTC revision 1.8 by jmc, Wed Oct 13 18:56:52 2004 UTC
# Line 1  Line 1 
1  % $Header$  % $Header$
2  % $Name$  % $Name$
3    
4  %\section{Time Integration}  
5    
6  \subsection{Non-linear free surface}  \subsection{Non-linear free surface}
7    \label{sect:nonlinear-freesurface}
8    
9  Recently, 2 options have added to the model  Recently, new options have been added to the model
 (and therefore, have not yet been extensively tested)  
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{Non-uniform linear-relation for the surface potential}
14    
15  The linear relation between  The linear relation between surface pressure/geo-potential
16  surface pressure / geo- potential ($\Phi_{surf}$)  ($\Phi_{surf}$) and surface displacement ($\eta$) could be considered
17  and surface displacement ($\eta$)  to be a constant ($b_s=$ constant)
18  has been considered as uniform ($b_s =$ Constant)  \marginpar{add a reference to part.1 here}
19  but is in fact  but is in fact dependent on the position ($x,y,r$)
 dependent on the position ($x,y,r$)  
20  since we linearize:  since we linearize:
21  $$\Phi_{surf}=\int_{R_o}^{R_o+\eta} b dr \simeq b_s \eta  $$\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}  ~\mathrm{with}~ b_s = b(\theta,S,r)_{r=R_o}
23  \simeq b_s(\theta_{ref}(R_o),S_{ref}(R_o),R_o)$$  \simeq b_s(\theta_{ref}(R_o),S_{ref}(R_o),R_o)$$
24  Note that, for convinience, the effect of the local  Note that, for convenience, the effect on $b_s$ of the local surface
25  surface $\theta,S$ on $b_s$  $\theta,S$ are not considered here, but are incorporated in to
26  are not considered here, but incorporated in $\Phi'_{hyd}$.  $\Phi'_{hyd}$.
27    
28  For the ocean, $b_s = g \rho_{surf} / \rho_o \simeq g$  For the ocean, $b_s = g \rho_{surf} / \rho_o \simeq g$ is a very good
29  is a fairly good approximation since the relative difference  approximation since the relative difference in surface density are
30  in surface density are usually small and only due to  usually small and only due to local $\theta,S$ gradients (because the
31  local $\theta,S$ gradient (because $R_o = 0$);  upper surface, $R_o = 0$, is essentially flat). Therefore, they can
32  Therefore, they can easely be incorporated in $\Phi'_{hyd}$.  easily be incorporated in $\Phi'_{hyd}$.
33    
34  For the atmosphere, because of topographic effects,  For the atmosphere, however, because of topographic effects, the
35  the reference surface pressure $R_o$ has large spacial differences  reference surface pressure $R_o$ has large spatial variations that
36  that are responsible for significant $b_s$ variations  are responsible for significant $b_s$ variations (from 0.8 to 1.2
37  (from 0.8 to 1.2 $[m^3/kg]$). For this reason,  $[m^3/kg]$). For this reason, we use a non-uniform linear coefficient
38  we use a non-uniform linear coefficient $b_s$.  $b_s$.
39    
40  Practically, in an oceanic configuration or  In practice, in an oceanic configuration or when the default value
41  when the default value (TRUE) of the parameter  (TRUE) of the parameter {\bf uniformLin\_PhiSurf} is used, then $b_s$
42  {\bf uniformLin\_PhiSurf} is used  is simply set to $g$ for the ocean and $1.$ for the atmosphere.
43  then $b_s$ is simply set to $g$ for the ocean  Turning {\bf uniformLin\_PhiSurf} to "FALSE", tells the code to
44  and $1.$ for the atmosphere.\\  evaluate $b_s$ from the reference vertical profile $\theta_{ref}$
45  Turning {\bf uniformLin\_PhiSurf} to "FALSE", allows to  ({\it S/R INI\_LINEAR\_PHISURF}) according to the reference surface
46  evaluate $b_s$ from the reference vertical profile $\theta_{ref}$  pressure $P_o$ ($=R_o$): $b_s = c_p \kappa (P_o / Pc)^{(\kappa - 1)}
47  ({\it S/R INI\_LINEAR\_PHISURF})  \theta_{ref}(P_o)$
48  according to the reference surface pressure $P_o$ ($=R_o$):  
 $b_s = c_p \kappa (P_o / Pc)^{(\kappa - 1)} \theta_{ref}(P_o)$  
49    
 %------------------------------------------  
50  \subsubsection{Free surface effect on column total thickness  \subsubsection{Free surface effect on column total thickness
51  (Non-linear free surface)}  (Non-linear free surface)}
52    
53  The total thickness of the fluid column is  The total thickness of the fluid column is $r_{surf} - R_{fixed} =
54  $r_{surf} - R_{min} = \eta + R_o - R_{min}$  \eta + R_o - R_{fixed}$. In most applications, the free surface
55  In the linear free surface approximation  displacements are small compared to the total thickness
56  (detailed before), only the fixed part of  $\eta << H_o = R_o - R_{fixed}$.
57  it ($R_o - R_{min})$ is considered when we integrate the  In the previous sections and in older version of the model,
58  continuity equation or compute tracer and momentum advection term.  the linearized free-surface approximation was made, assuming
59    $r_{surf} - R_{fixed} \simeq H_o$ when the horizontal transport is
60  This approximation is dropped when using  computed, either in the continuity equation or in tracer and momentum
61  the non-linear free surface formulation.  advection terms.
62  Details follow here after for the barotropic part  This approximation is dropped when using the non-linear free surface
63  and in the 2 next subsections for the baroclinic  formulation and the total thickness, including the time varying part
64  part.  $\eta$, is consisdered when computing horizontal transport.
65    Implications for the barotropic part are presented hereafter.
66  %------------------------------------------  In sections \ref{sect:freesurf-tracer-advection} and
67  % Non-Linear Barotropic part  \ref{sect:freesurf-momentum-advection}, consequences for tracer
68    and momentum are brifly discussed. a more detailed description
69  The continuous form of the model equations remains  is available in \cite{campin:02}.
70  unchanged, except for the 2D continuity equation  
71  (\ref{eq-tCsC-eta}) that is now integrated  
72  from $R_{min}(x,y)$ up to $r_{surf}=R_o+\eta$ :  In the non-linear formulation, the continuous form of the model equations
73    remains unchanged, except for the 2D continuity equation
74    (\ref{eq:discrete-time-backward-free-surface}) which is now
75    integrated from $R_{fixed}(x,y)$ up to $r_{surf}=R_o+\eta$ :
76    
77  \begin{displaymath}  \begin{displaymath}
78  \epsilon_{fs} \partial_t \eta =  \epsilon_{fs} \partial_t \eta =
79  \left. \dot{r} \right|_{r=r_{surf}} + \epsilon_{fw} (P-E) =  \left. \dot{r} \right|_{r=r_{surf}} + \epsilon_{fw} (P-E) =
80  - {\bf \nabla}_h \cdot \int_{R_{min}}^{R_o+\eta} \vec{\bf v} dr  - {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o+\eta} \vec{\bf v} dr
81  + \epsilon_{fw} (P-E)  + \epsilon_{fw} (P-E)
82  \end{displaymath}  \end{displaymath}
83    
84  Since $\eta$ has a direct effect on the horizontal  Since $\eta$ has a direct effect on the horizontal velocity (through
85  velocity (through $\nabla_h \Phi_{surf}$), this  $\nabla_h \Phi_{surf}$), this adds a non-linear term to the free
86  adds a non-linear term to the free surface equation.  surface equation. Several options for the time discretization of this
87    non-linear part can be considered, as detailed below.
88  Regarding the time discretization of this non-linear part,  
89  several options are now being tested:  If the column thickness is evaluated at time step $n$, and with
90    implicit treatment of the surface potential gradient, equations
91  With the column thickness evaluated at time step $n$,  (\ref{eq-solve2D} and \ref{eq-solve2D_rhs}) becomes:
 and the surface potential gradient still implicit,  
 equation (\ref{eq-solve2D} \& \ref{eq-solve2D_rhs})  
 become:  
92  \begin{eqnarray*}  \begin{eqnarray*}
93  \epsilon_{fs} {\eta}^{n+1} -  \epsilon_{fs} {\eta}^{n+1} -
94  {\bf \nabla}_h \cdot \Delta t^2 (\eta^{n}+R_o-R_{min})  {\bf \nabla}_h \cdot \Delta t^2 (\eta^{n}+R_o-R_{fixed})
95  {\bf \nabla}_h b_s {\eta}^{n+1}  {\bf \nabla}_h b_s {\eta}^{n+1}
96  = {\eta}^*  = {\eta}^*
 %\label{solve_2d}  
97  \end{eqnarray*}  \end{eqnarray*}
98  where  where
99  \begin{eqnarray*}  \begin{eqnarray*}
100  {\eta}^* = \epsilon_{fs} \: {\eta}^{n} -  {\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
101  \Delta t {\bf \nabla}_h \cdot \int_{R_{min}}^{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
102  \: + \: \epsilon_{fw} \Delta_t (P-E)^{n}  \: + \: \epsilon_{fw} \Delta_t (P-E)^{n}
 %\label{solve_2d_rhs}  
103  \end{eqnarray*}  \end{eqnarray*}
104  This method requires to update the solver matrix at each time step.  This method requires us to update the solver matrix at each time step.
105    
106  Alternatively, the non-linear contribution can be evaluated fully  Alternatively, the non-linear contribution can be evaluated fully
107  explicitly:  explicitly:
108  \begin{eqnarray*}  \begin{eqnarray*}
109  \epsilon_{fs} {\eta}^{n+1} -  \epsilon_{fs} {\eta}^{n+1} -
110  {\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{min})  {\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{fixed})
111  {\bf \nabla}_h b_s {\eta}^{n+1}  {\bf \nabla}_h b_s {\eta}^{n+1}
112  = {\eta}^*  = {\eta}^*
113  +{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n})  +{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n})
114  {\bf \nabla}_h b_s {\eta}^{n}  {\bf \nabla}_h b_s {\eta}^{n}
115  \end{eqnarray*}  \end{eqnarray*}
116  This formulation allows to keep the initial solver matrix  This formulation allows one to keep the initial solver matrix
117  since the non-linear free surface only affects the RHS.  unchanged though throughout the integration, since the non-linear free
118    surface only affects the RHS.
119    
120    Finally, another option is a "linearized" formulation where the total
121    column thickness appears only in the integral term of the RHS
122    (\ref{eq-solve2D_rhs}) but not directly in the equation
123    (\ref{eq-solve2D}).
124    
125    Those different options (see tab.?? for the one still available)
126    have been tested and show litle differences. However, we recommand
127    the use of the most precise method (the 1rst one) since the
128    computation cost involved in the solver matrix update are negligeable.
129    
130    \begin{center}
131     \begin{tabular}[htb]{|l|c|l|}
132       \hline
133       parameter & value & description \\
134       \hline
135                       & -1 & linear free-surface, restart from a pickup file \\
136                       &    & produced with \#undef EXACT\_CONSERV code\\
137       \cline{2-3}
138                       & 0 & Linear free-Surface \\
139       \cline{2-3}
140        nonlinFreeSurf & 4 & Non-linear free-surface \\
141       \cline{2-3}
142                       & 3 & same as 4 but neglecting
143                               $\int_{R_o}^{R_o+\eta} b' dr $ in $\Phi'_{hyd}$ \\
144       \cline{2-3}
145                       & 2 & same as 3 but do not update cg2d solver matrix \\
146       \cline{2-3}
147                      & 1 & same as 2 but treat momentum as in Linear FS \\
148       \hline
149                      & 0 & do not use $r*$ vertical coordinate (= default)\\
150       \cline{2-3}
151        select\_rStar & 2 & use $r^*$ vertical coordinate \\
152       \cline{2-3}
153                      & 1 & same as 2 but without the contribution of the\\
154                      &   & slope of the coordinate in $\nabla \Phi$ \\
155       \hline
156      \end{tabular}
157    \end{center}
158    
 An other option is a "linearized" formulation where the  
 total column thickness appears only in the integral term of  
 the RHS (\ref{eq-solve2D_rhs}) but not directly in  
 the equation (\ref{eq-solve2D}).  
159    
 %------------------------------------------  
160  \subsubsection{Free surface effect on the surface level thickness  \subsubsection{Free surface effect on the surface level thickness
161  (Non-linear free surface): Tracer advection}  (Non-linear free surface): Tracer advection}
162    \label{sect:freesurf-tracer-advection}
163    
164  To ensure a global tracer conservation (i.e., the total amount)  To ensure global tracer conservation (i.e., the total amount) as well
165  as well as the local one (see tracer section for more details),  as local conservation, the change in the surface level thickness must
166  the change in the surface level thickness must be consistent with  be consistent with the way the continuity equation is integrated, both
167  the way the continuity equation is integrated, both in  in the barotropic part (to find $\eta$) and baroclinic part (to find
168  in the barotropic part (to find $\eta$) and baroclinic part  $w = \dot{r}$).
 (to find $w = \dot{r}$).  
169    
170  To illustrate this, let's consider the shallow water model,  To illustrate this, consider the shallow water model, with uniform
171  with uniform cartesian horizontal grid:  Cartesian horizontal grid:
172  $$  $$
173  \partial_t h + \nabla \cdot h \vec{\bf v} = 0  \partial_t h + \nabla \cdot h \vec{\bf v} = 0
174  $$  $$
# Line 147  To conserve the tracer $\theta$ we have Line 177  To conserve the tracer $\theta$ we have
177  $$  $$
178  \partial_t (h \theta) + \nabla \cdot ( h \theta \vec{\bf v})= 0  \partial_t (h \theta) + \nabla \cdot ( h \theta \vec{\bf v})= 0
179  $$  $$
180  Using the implicit (non-linear) free surface described before, we have:  Using the implicit (non-linear) free surface described above (section
181    \ref{sect:pressure-method-linear-backward}) we have:
182  \begin{eqnarray*}  \begin{eqnarray*}
183  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} ) \\
184  \end{eqnarray*}  \end{eqnarray*}
185  The discretized form of the tracer equation must use the same  The discretized form of the tracer equation must adopt the same
186  "geometry" to compute the tracer fluxes, that is, the same value of  ``form'' in the computation of tracer fluxes, that is, the same value
187  h, as the continuity equation did:  of $h$, as used in the continuity equation:
188  \begin{eqnarray*}  \begin{eqnarray*}
189  h^{n+1} \, \theta^{n+1} = h^n \, \theta^n  h^{n+1} \, \theta^{n+1} = h^n \, \theta^n
190          - \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})
191  \end{eqnarray*}  \end{eqnarray*}
192    
193  In order to deal with the Adams-Bashforth time stepping,  The use of a 3 time-levels timestepping scheme such as the Adams-Bashforth
194  we implement this scheme slightly differently, in two step:  make the conservation less straitforward.
195  the variation of the water column thickness (from  The current implementation with the Adams-Bashforth time-stepping
196  $h^n$ to $h^{n+1}$)  provides an exact local conservation and prevents any drift in
197  is not incorporated directly to the tracer equation.  the global tracer content (\cite{campin:02}).
198  Instead,  Compared to the linear free-surface method, an additional step is required:
199  the model continues to evaluate the $G_\theta$ term (first step)  the variation of the water column thickness (from $h^n$ to $h^{n+1}$) is
200  as it use to do with the Linear free surface formulation  not incorporated directly into the tracer equation.  Instead, the
201  (with the "{\it surface correction}" turned "on", see tracer section):  model uses the $G_\theta$ terms (first step) as in the linear free
202    surface formulation (with the "{\it surface correction}" turned "on",
203    see tracer section):
204  $$  $$
205  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})
206           - \dot{r}_{surf}^{n+1} \theta^n \right) / h^n           - \dot{r}_{surf}^{n+1} \theta^n \right) / h^n
207  $$  $$
208  Then in a second step,  Then, in a second step, the thickness variation (expansion/reduction)
209  thickness variation (expansion/reduction) is taken into account :  is taken into account:
210  $$  $$
211  \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}} G_\theta^{(n+1/2)}
212  $$  $$
# Line 184  $ Line 217  $
217  = - \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) = (h^{n+1} - h^{n})  = - \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) = (h^{n+1} - h^{n})
218  / \Delta_t  / \Delta_t
219  $  $
220  those 2 formulations are equivalent.  these two formulations are equivalent.
221    
222    Implementation in the MITgcm is as follows.  The model ``geometry''
223    (here the {\bf hFacC,W,S}) is updated just before entering {\it
224    SOLVE\_FOR\_PRESSURE}, using the current $\eta$ field.  Then, at the
225    end of the time step, the variables are advanced in time, so that
226    $\eta^n$ replaces $\eta^{n-1}$.  At the next time step, the tracer
227    tendency ($G_\theta$) is computed using the same geometry, which is
228    now consistent with $\eta^{n-1}$.  Finally, in S/R {\it
229    TIMESTEP\_TRACER}, the expansion/reduction ratio is applied to the
230    surface level to compute the new tracer field.
231    
 The implementation in the MITgcm follows this scheme.  
 The model "geometry" (here the {\bf hFacC,W,S}) is updated  
 just before entering {\it SOLVE\_FOR\_PRESSURE},  
 according to the current $\eta$ field.  
 Then, at the end of the time step, the variables are  
 advanced in time, so that $\eta^n$ becomes $\eta^{n-1}$.  
 At the next time step, the tracer tendency ($G_\theta$) is computed  
 using the same geometry, that is now consistent with  
 $\eta^{n-1}$.  
 Finally, in S/R {\it TIMESTEP\_TRACER}, the expansion/reduction  
 ratio is applied to the surface level to compute the new tracer field.  
232    
 %------------------------------------------  
233  \subsubsection{Free surface effect on the surface level thickness  \subsubsection{Free surface effect on the surface level thickness
234  (Non-linear free surface): Momentum advection}      (Non-linear free surface): Momentum advection}    
235    \label{sect:freesurf-momentum-advection}
236    
237    With the flux form formulation, advection of momentum
238    can be treated exactly as the tracer advection is.
239    Here the expansion/reduction factors ($hFacW^{n+1}/hFacW^n$ for $u$
240    and $hFacS^{n+1}/hFacS^n$ for $v$) are simply applied in the
241    subroutine {\it TIMESTEP}.
242    
243  Regarding momentum advection,  Regarding momentum advection,
244  the vector invariant formulation is similar to the  the vector invariant formulation is similar to the
# Line 210  free surface effect on the surface level Line 248  free surface effect on the surface level
248  Updating the {\bf hFacC,W,S} and the {\bf recip\_hFac}(s)  Updating the {\bf hFacC,W,S} and the {\bf recip\_hFac}(s)
249  at one given place (like describe before) is sufficient.  at one given place (like describe before) is sufficient.
250    
251  With the flux form formulation, advection of momentum  \subsubsection{Non-linear free surface and vertical resolution}
252  can be treated exactly as the tracer advection is.  \label{sect:nonlin-freesurf-dz_surf}
253  Here the expansion/reduction factors ($hFacW^{n+1}/hFacW^n$ for $u$  
254  and $hFacS^{n+1}/hFacS^n$ for $v$) are simply applied in the  When the amplitude of the free-surface variations becomes
255  subroutine {\it TIMESTEP}.  as large as the vertical resolution near the surface,
256    the surface layer thickness can decrease to nearly zero or
257    can even vanishe completly.
258    This later possibility has not been implemented, and a
259    minimum relative thickness is imposed ({\bf hFacInf},
260    parameter file {\em data}, namelist {\em PARM01}) to prevent
261    numerical instabilities caused by very thin surface level.
262    
263    A better atlternative to the vanishing level problem has been
264    found and implemented recently, and rely on a different
265    vertical coordinate $r^*$~:
266    The time variation ot the total column thickness becomes
267    part of the r* coordinate motion, as in a $\sigma_{z},\sigma_{p}$
268    model, but the fixed part related to topography is treated
269    as in a height or pressure coordinate model.
270    A complete description is given in \cite{adcroft:04}.
271    
272    The time-stepping implementation of the $r^*$ coordinate is
273    identical to the non-linear free-surface in $r$ coordinate,
274    and differences appear only in the spacial discretisation.
275    \marginpar{needs a subsection ref. here}
276    

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

  ViewVC Help
Powered by ViewVC 1.1.22