/[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.7 by adcroft, Tue Nov 13 20:27:54 2001 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, two options have been added to the model (and have not yet
10    been extensively tested) that concern the free surface formulation.
11    
 Recently, 2 options have added to the model  
 (and therefore, have not yet been extensively tested)  
 that concern the free surface formulation.  
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 the linear free surface approximation
55  In the linear free surface approximation  (detailed before), only the fixed part of it ($R_o - R_{fixed})$ is
56  (detailed before), only the fixed part of  considered when we integrate the continuity equation or compute tracer
57  it ($R_o - R_{min})$ is considered when we integrate the  and momentum advection term.
58  continuity equation or compute tracer and momentum advection term.  
59    This approximation is dropped when using the non-linear free surface
60  This approximation is dropped when using  formulation.  Here we discuss sections the barotropic part. In
61  the non-linear free surface formulation.  sections \ref{sect:freesurf-tracer-advection} and
62  Details follow here after for the barotropic part  \ref{sect:freesurf-momentum-advection} we consider the baroclinic
63  and in the 2 next subsections for the baroclinic  component.
64  part.  
65    
66  %------------------------------------------  The continuous form of the model equations remains unchanged, except
67  % Non-Linear Barotropic part  for the 2D continuity equation (\ref{eq:discrete-time-backward-free-surface}) which is now
68    integrated from $R_{fixed}(x,y)$ up to $r_{surf}=R_o+\eta$ :
 The continuous form of the model equations remains  
 unchanged, except for the 2D continuity equation  
 (\ref{eq-tCsC-eta}) that is now integrated  
 from $R_{min}(x,y)$ up to $r_{surf}=R_o+\eta$ :  
69    
70  \begin{displaymath}  \begin{displaymath}
71  \epsilon_{fs} \partial_t \eta =  \epsilon_{fs} \partial_t \eta =
72  \left. \dot{r} \right|_{r=r_{surf}} + \epsilon_{fw} (P-E) =  \left. \dot{r} \right|_{r=r_{surf}} + \epsilon_{fw} (P-E) =
73  - {\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
74  + \epsilon_{fw} (P-E)  + \epsilon_{fw} (P-E)
75  \end{displaymath}  \end{displaymath}
76    
77  Since $\eta$ has a direct effect on the horizontal  Since $\eta$ has a direct effect on the horizontal velocity (through
78  velocity (through $\nabla_h \Phi_{surf}$), this  $\nabla_h \Phi_{surf}$), this adds a non-linear term to the free
79  adds a non-linear term to the free surface equation.  surface equation. Several options for the time discretization of this
80    non-linear part have been tested.
81  Regarding the time discretization of this non-linear part,  
82  several options are now being tested:  If the column thickness is evaluated at time step $n$, and with
83    implicit treatment of the surface potential gradient, equations
84  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:  
85  \begin{eqnarray*}  \begin{eqnarray*}
86  \epsilon_{fs} {\eta}^{n+1} -  \epsilon_{fs} {\eta}^{n+1} -
87  {\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})
88  {\bf \nabla}_h b_s {\eta}^{n+1}  {\bf \nabla}_h b_s {\eta}^{n+1}
89  = {\eta}^*  = {\eta}^*
 %\label{solve_2d}  
90  \end{eqnarray*}  \end{eqnarray*}
91  where  where
92  \begin{eqnarray*}  \begin{eqnarray*}
93  {\eta}^* = \epsilon_{fs} \: {\eta}^{n} -  {\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
94  \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
95  \: + \: \epsilon_{fw} \Delta_t (P-E)^{n}  \: + \: \epsilon_{fw} \Delta_t (P-E)^{n}
 %\label{solve_2d_rhs}  
96  \end{eqnarray*}  \end{eqnarray*}
97  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.
98    
99  Alternatively, the non-linear contribution can be evaluated fully  Alternatively, the non-linear contribution can be evaluated fully
100  explicitly:  explicitly:
101  \begin{eqnarray*}  \begin{eqnarray*}
102  \epsilon_{fs} {\eta}^{n+1} -  \epsilon_{fs} {\eta}^{n+1} -
103  {\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{min})  {\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{fixed})
104  {\bf \nabla}_h b_s {\eta}^{n+1}  {\bf \nabla}_h b_s {\eta}^{n+1}
105  = {\eta}^*  = {\eta}^*
106  +{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n})  +{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n})
107  {\bf \nabla}_h b_s {\eta}^{n}  {\bf \nabla}_h b_s {\eta}^{n}
108  \end{eqnarray*}  \end{eqnarray*}
109  This formulation allows to keep the initial solver matrix  This formulation allows one to keep the initial solver matrix
110  since the non-linear free surface only affects the RHS.  unchanged though throughout the integration, since the non-linear free
111    surface only affects the RHS.
112    
113    Finally, another option is a "linearized" formulation where the total
114    column thickness appears only in the integral term of the RHS
115    (\ref{eq-solve2D_rhs}) but not directly in the equation
116    (\ref{eq-solve2D}).
117    
 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}).  
118    
 %------------------------------------------  
119  \subsubsection{Free surface effect on the surface level thickness  \subsubsection{Free surface effect on the surface level thickness
120  (Non-linear free surface): Tracer advection}  (Non-linear free surface): Tracer advection}
121    \label{sect:freesurf-tracer-advection}
122    
123  To ensure a global tracer conservation (i.e., the total amount)  To ensure global tracer conservation (i.e., the total amount) as well
124  as well as the local one (see tracer section for more details),  as local conservation, the change in the surface level thickness must
125  the change in the surface level thickness must be consistent with  be consistent with the way the continuity equation is integrated, both
126  the way the continuity equation is integrated, both in  in the barotropic part (to find $\eta$) and baroclinic part (to find
127  in the barotropic part (to find $\eta$) and baroclinic part  $w = \dot{r}$).
 (to find $w = \dot{r}$).  
128    
129  To illustrate this, let's consider the shallow water model,  To illustrate this, consider the shallow water model, with uniform
130  with uniform cartesian horizontal grid:  Cartesian horizontal grid:
131  $$  $$
132  \partial_t h + \nabla \cdot h \vec{\bf v} = 0  \partial_t h + \nabla \cdot h \vec{\bf v} = 0
133  $$  $$
# Line 147  To conserve the tracer $\theta$ we have Line 136  To conserve the tracer $\theta$ we have
136  $$  $$
137  \partial_t (h \theta) + \nabla \cdot ( h \theta \vec{\bf v})= 0  \partial_t (h \theta) + \nabla \cdot ( h \theta \vec{\bf v})= 0
138  $$  $$
139  Using the implicit (non-linear) free surface described before, we have:  Using the implicit (non-linear) free surface described above (section
140    \ref{sect:pressure-method-linear-backward}) we have:
141  \begin{eqnarray*}  \begin{eqnarray*}
142  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} ) \\
143  \end{eqnarray*}  \end{eqnarray*}
144  The discretized form of the tracer equation must use the same  The discretized form of the tracer equation must adopt the same
145  "geometry" to compute the tracer fluxes, that is, the same value of  ``form'' in the computation of tracer fluxes, that is, the same value
146  h, as the continuity equation did:  of $h$, as used in the continuity equation:
147  \begin{eqnarray*}  \begin{eqnarray*}
148  h^{n+1} \, \theta^{n+1} = h^n \, \theta^n  h^{n+1} \, \theta^{n+1} = h^n \, \theta^n
149          - \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})
150  \end{eqnarray*}  \end{eqnarray*}
151    
152  In order to deal with the Adams-Bashforth time stepping,  For Adams-Bashforth time-stepping, we implement this scheme slightly
153  we implement this scheme slightly differently, in two step:  differently from the linear free-surface method, using two steps: the
154  the variation of the water column thickness (from  variation of the water column thickness (from $h^n$ to $h^{n+1}$) is
155  $h^n$ to $h^{n+1}$)  not incorporated directly into the tracer equation.  Instead, the
156  is not incorporated directly to the tracer equation.  model uses the $G_\theta$ terms (first step) as in the linear free
157  Instead,  surface formulation (with the "{\it surface correction}" turned "on",
158  the model continues to evaluate the $G_\theta$ term (first step)  see tracer section):
 as it use to do with the Linear free surface formulation  
 (with the "{\it surface correction}" turned "on", see tracer section):  
159  $$  $$
160  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})
161           - \dot{r}_{surf}^{n+1} \theta^n \right) / h^n           - \dot{r}_{surf}^{n+1} \theta^n \right) / h^n
162  $$  $$
163  Then in a second step,  Then, in a second step, the thickness variation (expansion/reduction)
164  thickness variation (expansion/reduction) is taken into account :  is taken into account:
165  $$  $$
166  \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)}
167  $$  $$
# Line 184  $ Line 172  $
172  = - \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})
173  / \Delta_t  / \Delta_t
174  $  $
175  those 2 formulations are equivalent.  these two formulations are equivalent.
176    
177    Implementation in the MITgcm is as follows.  The model ``geometry''
178    (here the {\bf hFacC,W,S}) is updated just before entering {\it
179    SOLVE\_FOR\_PRESSURE}, using the current $\eta$ field.  Then, at the
180    end of the time step, the variables are advanced in time, so that
181    $\eta^n$ replaces $\eta^{n-1}$.  At the next time step, the tracer
182    tendency ($G_\theta$) is computed using the same geometry, which is
183    now consistent with $\eta^{n-1}$.  Finally, in S/R {\it
184    TIMESTEP\_TRACER}, the expansion/reduction ratio is applied to the
185    surface level to compute the new tracer field.
186    
 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.  
187    
 %------------------------------------------  
188  \subsubsection{Free surface effect on the surface level thickness  \subsubsection{Free surface effect on the surface level thickness
189  (Non-linear free surface): Momentum advection}      (Non-linear free surface): Momentum advection}    
190    \label{sect:freesurf-momentum-advection}
191    
192  Regarding momentum advection,  Regarding momentum advection,
193  the vector invariant formulation is similar to the  the vector invariant formulation is similar to the
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22