1 |
adcroft |
1.7 |
% $Header: /u/gcmpack/mitgcmdoc/part2/nonlin_frsurf.tex,v 1.6 2001/11/13 20:13:54 adcroft Exp $ |
2 |
jmc |
1.1 |
% $Name: $ |
3 |
|
|
|
4 |
adcroft |
1.4 |
|
5 |
jmc |
1.1 |
|
6 |
|
|
\subsection{Non-linear free surface} |
7 |
adcroft |
1.6 |
\label{sect:nonlinear-freesurface} |
8 |
jmc |
1.1 |
|
9 |
adcroft |
1.4 |
Recently, two options have been added to the model (and have not yet |
10 |
|
|
been extensively tested) that concern the free surface formulation. |
11 |
|
|
|
12 |
jmc |
1.1 |
|
13 |
|
|
\subsubsection{Non-uniform linear-relation for the surface potential} |
14 |
|
|
|
15 |
adcroft |
1.4 |
The linear relation between surface pressure/geo-potential |
16 |
|
|
($\Phi_{surf}$) and surface displacement ($\eta$) could be considered |
17 |
|
|
to be a constant ($b_s=$ constant) |
18 |
jmc |
1.3 |
\marginpar{add a reference to part.1 here} |
19 |
|
|
but is in fact dependent on the position ($x,y,r$) |
20 |
jmc |
1.2 |
since we linearize: |
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} |
23 |
|
|
\simeq b_s(\theta_{ref}(R_o),S_{ref}(R_o),R_o)$$ |
24 |
adcroft |
1.4 |
Note that, for convenience, the effect on $b_s$ of the local surface |
25 |
|
|
$\theta,S$ are not considered here, but are incorporated in to |
26 |
|
|
$\Phi'_{hyd}$. |
27 |
|
|
|
28 |
|
|
For the ocean, $b_s = g \rho_{surf} / \rho_o \simeq g$ is a very good |
29 |
|
|
approximation since the relative difference in surface density are |
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 |
32 |
|
|
easily be incorporated in $\Phi'_{hyd}$. |
33 |
|
|
|
34 |
|
|
For the atmosphere, however, because of topographic effects, the |
35 |
cnh |
1.5 |
reference surface pressure $R_o$ has large spatial variations that |
36 |
adcroft |
1.4 |
are responsible for significant $b_s$ variations (from 0.8 to 1.2 |
37 |
|
|
$[m^3/kg]$). For this reason, we use a non-uniform linear coefficient |
38 |
|
|
$b_s$. |
39 |
|
|
|
40 |
|
|
In practice, in an oceanic configuration or when the default value |
41 |
|
|
(TRUE) of the parameter {\bf uniformLin\_PhiSurf} is used, then $b_s$ |
42 |
|
|
is simply set to $g$ for the ocean and $1.$ for the atmosphere. |
43 |
|
|
Turning {\bf uniformLin\_PhiSurf} to "FALSE", tells the code to |
44 |
|
|
evaluate $b_s$ from the reference vertical profile $\theta_{ref}$ |
45 |
|
|
({\it S/R INI\_LINEAR\_PHISURF}) according to the reference surface |
46 |
|
|
pressure $P_o$ ($=R_o$): $b_s = c_p \kappa (P_o / Pc)^{(\kappa - 1)} |
47 |
|
|
\theta_{ref}(P_o)$ |
48 |
|
|
|
49 |
jmc |
1.1 |
|
50 |
|
|
\subsubsection{Free surface effect on column total thickness |
51 |
|
|
(Non-linear free surface)} |
52 |
|
|
|
53 |
adcroft |
1.4 |
The total thickness of the fluid column is $r_{surf} - R_{fixed} = |
54 |
|
|
\eta + R_o - R_{fixed}$ In the linear free surface approximation |
55 |
|
|
(detailed before), only the fixed part of it ($R_o - R_{fixed})$ is |
56 |
|
|
considered when we integrate the continuity equation or compute tracer |
57 |
|
|
and momentum advection term. |
58 |
|
|
|
59 |
|
|
This approximation is dropped when using the non-linear free surface |
60 |
|
|
formulation. Here we discuss sections the barotropic part. In |
61 |
|
|
sections \ref{sect:freesurf-tracer-advection} and |
62 |
|
|
\ref{sect:freesurf-momentum-advection} we consider the baroclinic |
63 |
|
|
component. |
64 |
|
|
|
65 |
|
|
|
66 |
|
|
The continuous form of the model equations remains unchanged, except |
67 |
adcroft |
1.7 |
for the 2D continuity equation (\ref{eq:discrete-time-backward-free-surface}) which is now |
68 |
adcroft |
1.4 |
integrated from $R_{fixed}(x,y)$ up to $r_{surf}=R_o+\eta$ : |
69 |
jmc |
1.1 |
|
70 |
|
|
\begin{displaymath} |
71 |
|
|
\epsilon_{fs} \partial_t \eta = |
72 |
|
|
\left. \dot{r} \right|_{r=r_{surf}} + \epsilon_{fw} (P-E) = |
73 |
jmc |
1.3 |
- {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o+\eta} \vec{\bf v} dr |
74 |
jmc |
1.1 |
+ \epsilon_{fw} (P-E) |
75 |
|
|
\end{displaymath} |
76 |
|
|
|
77 |
adcroft |
1.4 |
Since $\eta$ has a direct effect on the horizontal velocity (through |
78 |
|
|
$\nabla_h \Phi_{surf}$), this adds a non-linear term to the free |
79 |
|
|
surface equation. Several options for the time discretization of this |
80 |
|
|
non-linear part have been tested. |
81 |
|
|
|
82 |
|
|
If the column thickness is evaluated at time step $n$, and with |
83 |
|
|
implicit treatment of the surface potential gradient, equations |
84 |
|
|
(\ref{eq-solve2D} and \ref{eq-solve2D_rhs}) becomes: |
85 |
jmc |
1.1 |
\begin{eqnarray*} |
86 |
|
|
\epsilon_{fs} {\eta}^{n+1} - |
87 |
jmc |
1.3 |
{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n}+R_o-R_{fixed}) |
88 |
jmc |
1.2 |
{\bf \nabla}_h b_s {\eta}^{n+1} |
89 |
jmc |
1.1 |
= {\eta}^* |
90 |
|
|
\end{eqnarray*} |
91 |
|
|
where |
92 |
|
|
\begin{eqnarray*} |
93 |
|
|
{\eta}^* = \epsilon_{fs} \: {\eta}^{n} - |
94 |
jmc |
1.3 |
\Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o+\eta^n} \vec{\bf v}^* dr |
95 |
jmc |
1.1 |
\: + \: \epsilon_{fw} \Delta_t (P-E)^{n} |
96 |
|
|
\end{eqnarray*} |
97 |
adcroft |
1.4 |
This method requires us to update the solver matrix at each time step. |
98 |
jmc |
1.1 |
|
99 |
|
|
Alternatively, the non-linear contribution can be evaluated fully |
100 |
|
|
explicitly: |
101 |
|
|
\begin{eqnarray*} |
102 |
|
|
\epsilon_{fs} {\eta}^{n+1} - |
103 |
jmc |
1.3 |
{\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{fixed}) |
104 |
jmc |
1.2 |
{\bf \nabla}_h b_s {\eta}^{n+1} |
105 |
jmc |
1.1 |
= {\eta}^* |
106 |
jmc |
1.2 |
+{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n}) |
107 |
|
|
{\bf \nabla}_h b_s {\eta}^{n} |
108 |
jmc |
1.1 |
\end{eqnarray*} |
109 |
adcroft |
1.4 |
This formulation allows one to keep the initial solver matrix |
110 |
|
|
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 |
jmc |
1.1 |
|
118 |
|
|
|
119 |
|
|
\subsubsection{Free surface effect on the surface level thickness |
120 |
|
|
(Non-linear free surface): Tracer advection} |
121 |
adcroft |
1.4 |
\label{sect:freesurf-tracer-advection} |
122 |
jmc |
1.1 |
|
123 |
adcroft |
1.4 |
To ensure global tracer conservation (i.e., the total amount) as well |
124 |
|
|
as local conservation, the change in the surface level thickness must |
125 |
|
|
be consistent with the way the continuity equation is integrated, both |
126 |
|
|
in the barotropic part (to find $\eta$) and baroclinic part (to find |
127 |
|
|
$w = \dot{r}$). |
128 |
jmc |
1.1 |
|
129 |
adcroft |
1.4 |
To illustrate this, consider the shallow water model, with uniform |
130 |
|
|
Cartesian horizontal grid: |
131 |
jmc |
1.1 |
$$ |
132 |
|
|
\partial_t h + \nabla \cdot h \vec{\bf v} = 0 |
133 |
|
|
$$ |
134 |
|
|
where $h$ is the total thickness of the water column. |
135 |
|
|
To conserve the tracer $\theta$ we have to discretize: |
136 |
|
|
$$ |
137 |
|
|
\partial_t (h \theta) + \nabla \cdot ( h \theta \vec{\bf v})= 0 |
138 |
|
|
$$ |
139 |
adcroft |
1.4 |
Using the implicit (non-linear) free surface described above (section |
140 |
adcroft |
1.6 |
\ref{sect:pressure-method-linear-backward}) we have: |
141 |
jmc |
1.1 |
\begin{eqnarray*} |
142 |
jmc |
1.2 |
h^{n+1} = h^{n} - \Delta_t \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) \\ |
143 |
jmc |
1.1 |
\end{eqnarray*} |
144 |
adcroft |
1.4 |
The discretized form of the tracer equation must adopt the same |
145 |
|
|
``form'' in the computation of tracer fluxes, that is, the same value |
146 |
|
|
of $h$, as used in the continuity equation: |
147 |
jmc |
1.1 |
\begin{eqnarray*} |
148 |
|
|
h^{n+1} \, \theta^{n+1} = h^n \, \theta^n |
149 |
jmc |
1.2 |
- \Delta_t \nabla \cdot (h^n \, \theta^n \, \vec{\bf v}^{n+1}) |
150 |
jmc |
1.1 |
\end{eqnarray*} |
151 |
|
|
|
152 |
adcroft |
1.4 |
For Adams-Bashforth time-stepping, we implement this scheme slightly |
153 |
|
|
differently from the linear free-surface method, using two steps: the |
154 |
|
|
variation of the water column thickness (from $h^n$ to $h^{n+1}$) is |
155 |
|
|
not incorporated directly into the tracer equation. Instead, the |
156 |
|
|
model uses the $G_\theta$ terms (first step) as in the linear free |
157 |
|
|
surface formulation (with the "{\it surface correction}" turned "on", |
158 |
|
|
see tracer section): |
159 |
jmc |
1.1 |
$$ |
160 |
jmc |
1.2 |
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 |
162 |
jmc |
1.1 |
$$ |
163 |
adcroft |
1.4 |
Then, in a second step, the thickness variation (expansion/reduction) |
164 |
|
|
is taken into account: |
165 |
jmc |
1.1 |
$$ |
166 |
jmc |
1.2 |
\theta^{n+1} = \theta^n + \Delta_t \frac{h^n}{h^{n+1}} G_\theta^{(n+1/2)} |
167 |
jmc |
1.1 |
$$ |
168 |
|
|
Note that with a simple forward time step (no Adams-Bashforth), |
169 |
|
|
since |
170 |
|
|
$ |
171 |
jmc |
1.2 |
\dot{r}_{surf}^{n+1} |
172 |
|
|
= - \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) = (h^{n+1} - h^{n}) |
173 |
jmc |
1.1 |
/ \Delta_t |
174 |
|
|
$ |
175 |
adcroft |
1.4 |
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 |
jmc |
1.1 |
|
187 |
|
|
|
188 |
|
|
\subsubsection{Free surface effect on the surface level thickness |
189 |
|
|
(Non-linear free surface): Momentum advection} |
190 |
adcroft |
1.4 |
\label{sect:freesurf-momentum-advection} |
191 |
jmc |
1.1 |
|
192 |
|
|
Regarding momentum advection, |
193 |
|
|
the vector invariant formulation is similar to the |
194 |
|
|
advective form (as opposed to the flux form) and therefore |
195 |
|
|
does not need specific modification to include the |
196 |
|
|
free surface effect on the surface level thickness. |
197 |
|
|
Updating the {\bf hFacC,W,S} and the {\bf recip\_hFac}(s) |
198 |
|
|
at one given place (like describe before) is sufficient. |
199 |
|
|
|
200 |
|
|
With the flux form formulation, advection of momentum |
201 |
jmc |
1.2 |
can be treated exactly as the tracer advection is. |
202 |
jmc |
1.1 |
Here the expansion/reduction factors ($hFacW^{n+1}/hFacW^n$ for $u$ |
203 |
|
|
and $hFacS^{n+1}/hFacS^n$ for $v$) are simply applied in the |
204 |
|
|
subroutine {\it TIMESTEP}. |
205 |
|
|
|