| 1 |
% $Header$ |
% $Header$ |
| 2 |
% $Name$ |
% $Name$ |
| 3 |
|
|
|
%\section{Time Integration} |
|
| 4 |
|
|
|
\subsection{Non-linear free surface} |
|
| 5 |
|
|
| 6 |
Recently, 2 options have added to the model |
\subsection{Non-linear free-surface} |
| 7 |
(and therefore, have not yet been extensively tested) |
\label{sec:nonlinear-freesurface} |
| 8 |
|
|
| 9 |
|
Recently, new options have been added to the model |
| 10 |
that concern the free surface formulation. |
that concern the free surface formulation. |
| 11 |
|
|
|
%------------------------------------------ |
|
|
\subsubsection{Non-uniform linear-relation for the surface potential} |
|
| 12 |
|
|
| 13 |
The linear relation between |
\subsubsection{pressure/geo-potential and free surface} |
| 14 |
surface pressure / geo- potential ($\Phi_{surf}$) |
\label{sec:phi-freesurface} |
| 15 |
and surface displacement ($\eta$) |
|
| 16 |
has been considered as uniform ($b_s =$ Constant) |
For the atmosphere, since $\phi = \phi_{topo} - \int^p_{p_s} \alpha dp$, |
| 17 |
but is in fact |
subtracting the reference state defined in section |
| 18 |
dependent on the position ($x,y,r$) |
\ref{sec:hpe-p-geo-potential-split}~:\\ |
| 19 |
since we linearize: |
$$ |
| 20 |
$$\Phi_{surf}=\int_{R_o}^{R_o+\eta} b dr \simeq b_s \eta |
\phi_o = \phi_{topo} - \int^p_{p_o} \alpha_o dp |
| 21 |
~\mathrm{with}~ b_s = b(\theta,S,r)_{r=R_o} |
\hspace{5mm}\mathrm{with}\hspace{3mm} \phi_o(p_o)=\phi_{topo} |
| 22 |
\simeq b_s(\theta_{ref}(R_o),S_{ref}(R_o),R_o)$$ |
$$ |
| 23 |
Note that, for convinience, the effect of the local |
we get: |
| 24 |
surface $\theta,S$ on $b_s$ |
$$ |
| 25 |
are not considered here, but incorporated in $\Phi'_{hyd}$. |
\phi' = \phi - \phi_o = \int^{p_s}_p \alpha dp - \int^{p_o}_p \alpha_o dp |
| 26 |
|
$$ |
| 27 |
For the ocean, $b_s = g \rho_{surf} / \rho_o \simeq g$ |
For the ocean, the reference state is simpler since $\rho_c$ does not dependent |
| 28 |
is a fairly good approximation since the relative difference |
on $z$ ($b_o=g$) and the surface reference position is uniformly $z=0$ ($R_o=0$), |
| 29 |
in surface density are usually small and only due to |
and the same subtraction leads to a similar relation. |
| 30 |
local $\theta,S$ gradient (because $R_o = 0$); |
For both fluid, using the isomorphic notations, we can write: |
| 31 |
Therefore, they can easely be incorporated in $\Phi'_{hyd}$. |
$$ |
| 32 |
|
\phi' = \int^{r_{surf}}_r b~ dr - \int^{R_o}_r b_o dr |
| 33 |
For the atmosphere, because of topographic effects, |
$$ |
| 34 |
the reference surface pressure $R_o$ has large spacial differences |
and re-write as: |
| 35 |
that are responsible for significant $b_s$ variations |
\begin{equation} |
| 36 |
(from 0.8 to 1.2 $[m^3/kg]$). For this reason, |
\phi' = \int^{r_{surf}}_{R_o} b~ dr + \int^{R_o}_r (b - b_o) dr |
| 37 |
we use a non-uniform linear coefficient $b_s$. |
\label{eq:split-phi-Ro} |
| 38 |
|
\end{equation} |
| 39 |
Practically, in an oceanic configuration or |
or: |
| 40 |
when the default value (TRUE) of the parameter |
\begin{equation} |
| 41 |
{\bf uniformLin\_PhiSurf} is used |
\phi' = \int^{r_{surf}}_{R_o} b_o dr + \int^{r_{surf}}_r (b - b_o) dr |
| 42 |
then $b_s$ is simply set to $g$ for the ocean |
\label{eq:split-phi-bo} |
| 43 |
and $1.$ for the atmosphere.\\ |
\end{equation} |
| 44 |
Turning {\bf uniformLin\_PhiSurf} to "FALSE", allows to |
|
| 45 |
evaluate $b_s$ from the reference vertical profile $\theta_{ref}$ |
In section \ref{sec:finding_the_pressure_field}, following eq.\ref{eq:split-phi-Ro}, |
| 46 |
({\it S/R INI\_LINEAR\_PHISURF}) |
the pressure/geo-potential $\phi'$ has been separated into surface ($\phi_s$), |
| 47 |
according to the reference surface pressure $P_o$ ($=R_o$): |
and hydrostatic anomaly ($\phi'_{hyd}$). |
| 48 |
$b_s = c_p \kappa (P_o / Pc)^{(\kappa - 1)} \theta_{ref}(P_o)$ |
In this section, the split between $\phi_s$ and $\phi'_{hyd}$ is |
| 49 |
|
made according to equation \ref{eq:split-phi-bo}. This slightly |
| 50 |
|
different definition reflects the actual implementation in the code |
| 51 |
|
and is valid for both linear and non-linear |
| 52 |
|
free-surface formulation, in both r-coordinate and r*-coordinate. |
| 53 |
|
|
| 54 |
|
Because the linear free-surface approximation ignore the tracer content |
| 55 |
|
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}$~: |
| 57 |
|
$$ |
| 58 |
|
\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}$ |
| 61 |
|
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$ |
| 63 |
|
and $\phi'_{hyd}=\int^{R_o}_r b' dr$.\\ |
| 64 |
|
On the contrary, the unapproximated formulation ("non-linear free-surface", |
| 65 |
|
see the next section) retains the full expression: |
| 66 |
|
$\phi'_{hyd} = \int^{r_{surf}}_r (b - b_o) dr $~. |
| 67 |
|
This is obtained by selecting {\bf nonlinFreeSurf}=4 in parameter |
| 68 |
|
file {\em data}.\\ |
| 69 |
|
|
| 70 |
|
Regarding the surface potential: |
| 71 |
|
$$\phi_s = \int_{R_o}^{R_o+\eta} b_o dr = b_s \eta |
| 72 |
|
\hspace{5mm}\mathrm{with}\hspace{5mm} |
| 73 |
|
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 |
| 75 |
|
the usual numerical truncation, since generally $|\eta|$ is smaller |
| 76 |
|
than the vertical grid increment). |
| 77 |
|
|
| 78 |
|
For the ocean, $\phi_s = g \eta$ and $b_s = g$ is uniform. |
| 79 |
|
For the atmosphere, however, because of topographic effects, the |
| 80 |
|
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 |
| 82 |
|
$[m^3/kg]$). For this reason, when {\bf uniformLin\_PhiSurf} {\em=.FALSE.} |
| 83 |
|
(parameter file {\em data}, namelist {\em PARAM01}) |
| 84 |
|
a non-uniform linear coefficient $b_s$ is used and computed |
| 85 |
|
({\it S/R INI\_LINEAR\_PHISURF}) according to the reference surface |
| 86 |
|
pressure $p_o$: |
| 87 |
|
$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. |
| 89 |
|
|
| 90 |
|
|
|
%------------------------------------------ |
|
| 91 |
\subsubsection{Free surface effect on column total thickness |
\subsubsection{Free surface effect on column total thickness |
| 92 |
(Non-linear free surface)} |
(Non-linear free-surface)} |
| 93 |
|
|
| 94 |
The total thickness of the fluid column is |
The total thickness of the fluid column is $r_{surf} - R_{fixed} = |
| 95 |
$r_{surf} - R_{min} = \eta + R_o - R_{min}$ |
\eta + R_o - R_{fixed}$. In most applications, the free surface |
| 96 |
In the linear free surface approximation |
displacements are small compared to the total thickness |
| 97 |
(detailed before), only the fixed part of |
$\eta \ll H_o = R_o - R_{fixed}$. |
| 98 |
it ($R_o - R_{min})$ is considered when we integrate the |
In the previous sections and in older version of the model, |
| 99 |
continuity equation or compute tracer and momentum advection term. |
the linearized free-surface approximation was made, assuming |
| 100 |
|
$r_{surf} - R_{fixed} \simeq H_o$ when computing horizontal transports, |
| 101 |
This approximation is dropped when using |
either in the continuity equation or in tracer and momentum |
| 102 |
the non-linear free surface formulation. |
advection terms. |
| 103 |
Details follow here after for the barotropic part |
This approximation is dropped when using the non-linear free-surface |
| 104 |
and in the 2 next subsections for the baroclinic |
formulation and the total thickness, including the time varying part |
| 105 |
part. |
$\eta$, is considered when computing horizontal transports. |
| 106 |
|
Implications for the barotropic part are presented hereafter. |
| 107 |
%------------------------------------------ |
In section \ref{sec:freesurf-tracer-advection} consequences for |
| 108 |
% Non-Linear Barotropic part |
tracer conservation is briefly discussed (more details can be |
| 109 |
|
found in \cite{campin:02})~; the general time-stepping is presented |
| 110 |
The continuous form of the model equations remains |
in section \ref{sec:nonlin-freesurf-timestepping} with some |
| 111 |
unchanged, except for the 2D continuity equation |
limitations regarding the vertical resolution in section |
| 112 |
(\ref{eq-tCsC-eta}) that is now integrated |
\ref{sec:nonlin-freesurf-dz_surf}. |
| 113 |
from $R_{min}(x,y)$ up to $r_{surf}=R_o+\eta$ : |
|
| 114 |
|
In the non-linear formulation, the continuous form of the model |
| 115 |
|
equations remains unchanged, except for the 2D continuity equation |
| 116 |
|
(\ref{eq:discrete-time-backward-free-surface}) which is now |
| 117 |
|
integrated from $R_{fixed}(x,y)$ up to $r_{surf}=R_o+\eta$ : |
| 118 |
|
|
| 119 |
\begin{displaymath} |
\begin{displaymath} |
| 120 |
\epsilon_{fs} \partial_t \eta = |
\epsilon_{fs} \partial_t \eta = |
| 121 |
\left. \dot{r} \right|_{r=r_{surf}} + \epsilon_{fw} (P-E) = |
\left. \dot{r} \right|_{r=r_{surf}} + \epsilon_{fw} (P-E) = |
| 122 |
- {\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 |
| 123 |
+ \epsilon_{fw} (P-E) |
+ \epsilon_{fw} (P-E) |
| 124 |
\end{displaymath} |
\end{displaymath} |
| 125 |
|
|
| 126 |
Since $\eta$ has a direct effect on the horizontal |
Since $\eta$ has a direct effect on the horizontal velocity (through |
| 127 |
velocity (through $\nabla_h \Phi_{surf}$), this |
$\nabla_h \Phi_{surf}$), this adds a non-linear term to the free |
| 128 |
adds a non-linear term to the free surface equation. |
surface equation. Several options for the time discretization of this |
| 129 |
|
non-linear part can be considered, as detailed below. |
| 130 |
Regarding the time discretization of this non-linear part, |
|
| 131 |
several options are now being tested: |
If the column thickness is evaluated at time step $n$, and with |
| 132 |
|
implicit treatment of the surface potential gradient, equations |
| 133 |
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: |
|
| 134 |
\begin{eqnarray*} |
\begin{eqnarray*} |
| 135 |
\epsilon_{fs} {\eta}^{n+1} - |
\epsilon_{fs} {\eta}^{n+1} - |
| 136 |
{\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}) |
| 137 |
{\bf \nabla}_h b_s {\eta}^{n+1} |
{\bf \nabla}_h b_s {\eta}^{n+1} |
| 138 |
= {\eta}^* |
= {\eta}^* |
|
%\label{solve_2d} |
|
| 139 |
\end{eqnarray*} |
\end{eqnarray*} |
| 140 |
where |
where |
| 141 |
\begin{eqnarray*} |
\begin{eqnarray*} |
| 142 |
{\eta}^* = \epsilon_{fs} \: {\eta}^{n} - |
{\eta}^* = \epsilon_{fs} \: {\eta}^{n} - |
| 143 |
\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 |
| 144 |
\: + \: \epsilon_{fw} \Delta_t (P-E)^{n} |
\: + \: \epsilon_{fw} \Delta_t (P-E)^{n} |
| 145 |
%\label{solve_2d_rhs} |
\end{eqnarray*} |
| 146 |
\end{eqnarray*} |
This method requires us to update the solver matrix at each time step. |
|
This method requires 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 |
| 149 |
explicitly: |
explicitly: |
| 150 |
\begin{eqnarray*} |
\begin{eqnarray*} |
| 151 |
\epsilon_{fs} {\eta}^{n+1} - |
\epsilon_{fs} {\eta}^{n+1} - |
| 152 |
{\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{min}) |
{\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{fixed}) |
| 153 |
{\bf \nabla}_h b_s {\eta}^{n+1} |
{\bf \nabla}_h b_s {\eta}^{n+1} |
| 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 to keep the initial solver matrix |
This formulation allows one to keep the initial solver matrix |
| 159 |
since the non-linear free surface only affects the RHS. |
unchanged though throughout the integration, since the non-linear free |
| 160 |
|
surface only affects the RHS. |
| 161 |
An other option is a "linearized" formulation where the |
|
| 162 |
total column thickness appears only in the integral term of |
Finally, another option is a "linearized" formulation where the total |
| 163 |
the RHS (\ref{eq-solve2D_rhs}) but not directly in |
column thickness appears only in the integral term of the RHS |
| 164 |
the equation (\ref{eq-solve2D}). |
(\ref{eq-solve2D_rhs}) but not directly in the equation |
| 165 |
|
(\ref{eq-solve2D}). |
| 166 |
%------------------------------------------ |
|
| 167 |
\subsubsection{Free surface effect on the surface level thickness |
Those different options (see Table \ref{tab:nonLinFreeSurf_flags}) |
| 168 |
(Non-linear free surface): Tracer advection} |
have been tested and show little differences. However, we recommend |
| 169 |
|
the use of the most precise method (the 1rst one) since the |
| 170 |
To ensure a global tracer conservation (i.e., the total amount) |
computation cost involved in the solver matrix update is negligible. |
| 171 |
as well as the local one (see tracer section for more details), |
|
| 172 |
the change in the surface level thickness must be consistent with |
\begin{table}[htb] |
| 173 |
the way the continuity equation is integrated, both in |
%\begin{center} |
| 174 |
in the barotropic part (to find $\eta$) and baroclinic part |
\centering |
| 175 |
(to find $w = \dot{r}$). |
\begin{tabular}[htb]{|l|c|l|} |
| 176 |
|
\hline |
| 177 |
|
parameter & value & description \\ |
| 178 |
|
\hline |
| 179 |
|
& -1 & linear free-surface, restart from a pickup file \\ |
| 180 |
|
& & produced with \#undef EXACT\_CONSERV code\\ |
| 181 |
|
\cline{2-3} |
| 182 |
|
& 0 & Linear free-surface \\ |
| 183 |
|
\cline{2-3} |
| 184 |
|
nonlinFreeSurf & 4 & Non-linear free-surface \\ |
| 185 |
|
\cline{2-3} |
| 186 |
|
& 3 & same as 4 but neglecting |
| 187 |
|
$\int_{R_o}^{R_o+\eta} b' dr $ in $\Phi'_{hyd}$ \\ |
| 188 |
|
\cline{2-3} |
| 189 |
|
& 2 & same as 3 but do not update cg2d solver matrix \\ |
| 190 |
|
\cline{2-3} |
| 191 |
|
& 1 & same as 2 but treat momentum as in Linear FS \\ |
| 192 |
|
\hline |
| 193 |
|
& 0 & do not use $r*$ vertical coordinate (= default)\\ |
| 194 |
|
\cline{2-3} |
| 195 |
|
select\_rStar & 2 & use $r^*$ vertical coordinate \\ |
| 196 |
|
\cline{2-3} |
| 197 |
|
& 1 & same as 2 but without the contribution of the\\ |
| 198 |
|
& & slope of the coordinate in $\nabla \Phi$ \\ |
| 199 |
|
\hline |
| 200 |
|
\end{tabular} |
| 201 |
|
\caption{Non-linear free-surface flags} |
| 202 |
|
\label{tab:nonLinFreeSurf_flags} |
| 203 |
|
%\end{center} |
| 204 |
|
\end{table} |
| 205 |
|
|
| 206 |
|
|
| 207 |
|
\subsubsection{Tracer conservation with non-linear free-surface} |
| 208 |
|
\label{sec:freesurf-tracer-advection} |
| 209 |
|
|
| 210 |
|
To ensure global tracer conservation (i.e., the total amount) as well |
| 211 |
|
as local conservation, the change in the surface level thickness must |
| 212 |
|
be consistent with the way the continuity equation is integrated, both |
| 213 |
|
in the barotropic part (to find $\eta$) and baroclinic part (to find |
| 214 |
|
$w = \dot{r}$). |
| 215 |
|
|
| 216 |
To illustrate this, let's consider the shallow water model, |
To illustrate this, consider the shallow water model, with a source |
| 217 |
with uniform cartesian horizontal grid: |
of fresh water (P): |
| 218 |
$$ |
$$ |
| 219 |
\partial_t h + \nabla \cdot h \vec{\bf v} = 0 |
\partial_t h + \nabla \cdot h \vec{\bf v} = P |
| 220 |
$$ |
$$ |
| 221 |
where $h$ is the total thickness of the water column. |
where $h$ is the total thickness of the water column. |
| 222 |
To conserve the tracer $\theta$ we have to discretize: |
To conserve the tracer $\theta$ we have to discretize: |
| 223 |
$$ |
$$ |
| 224 |
\partial_t (h \theta) + \nabla \cdot ( h \theta \vec{\bf v})= 0 |
\partial_t (h \theta) + \nabla \cdot ( h \theta \vec{\bf v}) |
| 225 |
|
= P \theta_{\mathrm{rain}} |
| 226 |
$$ |
$$ |
| 227 |
Using the implicit (non-linear) free surface described before, we have: |
Using the implicit (non-linear) free surface described above (section |
| 228 |
|
\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} ) \\ |
h^{n+1} = h^{n} - \Delta t \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) + \Delta t P \\ |
| 231 |
\end{eqnarray*} |
\end{eqnarray*} |
| 232 |
The discretized form of the tracer equation must use the same |
The discretized form of the tracer equation must adopt the same |
| 233 |
"geometry" to compute the tracer fluxes, that is, the same value of |
``form'' in the computation of tracer fluxes, that is, the same value |
| 234 |
h, as the continuity equation did: |
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} |
| 239 |
\end{eqnarray*} |
\end{eqnarray*} |
| 240 |
|
|
| 241 |
In order to deal with the Adams-Bashforth time stepping, |
The use of a 3 time-levels time-stepping scheme such as the Adams-Bashforth |
| 242 |
we implement this scheme slightly differently, in two step: |
make the conservation sightly tricky. |
| 243 |
the variation of the water column thickness (from |
The current implementation with the Adams-Bashforth time-stepping |
| 244 |
$h^n$ to $h^{n+1}$) |
provides an exact local conservation and prevents any drift in |
| 245 |
is not incorporated directly to the tracer equation. |
the global tracer content (\cite{campin:02}). |
| 246 |
Instead, |
Compared to the linear free-surface method, an additional step is required: |
| 247 |
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 |
| 248 |
as it use to do with the Linear free surface formulation |
not incorporated directly into the tracer equation. Instead, the |
| 249 |
(with the "{\it surface correction}" turned "on", see tracer section): |
model uses the $G_\theta$ terms (first step) as in the linear free |
| 250 |
|
surface formulation (with the "{\it surface correction}" turned "on", |
| 251 |
|
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, |
Then, in a second step, the thickness variation (expansion/reduction) |
| 257 |
thickness variation (expansion/reduction) is taken into account : |
is taken into account: |
| 258 |
$$ |
$$ |
| 259 |
\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}} |
| 260 |
|
\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}} |
| 262 |
|
% \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, |
| 266 |
since |
since |
| 267 |
$ |
$ |
| 268 |
\dot{r}_{surf}^{n+1} |
(h^{n+1} - h^{n})/ \Delta t = |
| 269 |
= - \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 |
|
| 270 |
$ |
$ |
|
those 2 formulations are equivalent. |
|
| 271 |
|
|
| 272 |
The implementation in the MITgcm follows this scheme. |
\subsubsection{Time stepping implementation of the |
| 273 |
The model "geometry" (here the {\bf hFacC,W,S}) is updated |
non-linear free-surface} |
| 274 |
just before entering {\it SOLVE\_FOR\_PRESSURE}, |
\label{sec:nonlin-freesurf-timestepping} |
| 275 |
according to the current $\eta$ field. |
|
| 276 |
Then, at the end of the time step, the variables are |
The grid cell thickness was hold constant with the linear |
| 277 |
advanced in time, so that $\eta^n$ becomes $\eta^{n-1}$. |
free-surface~; with the non-linear free-surface, it is now varying |
| 278 |
At the next time step, the tracer tendency ($G_\theta$) is computed |
in time, at least at the surface level. |
| 279 |
using the same geometry, that is now consistent with |
This implies some modifications of the general algorithm described |
| 280 |
$\eta^{n-1}$. |
earlier in sections \ref{sec:adams-bashforth-sync} and |
| 281 |
Finally, in S/R {\it TIMESTEP\_TRACER}, the expansion/reduction |
\ref{sec:adams-bashforth-staggered}. |
| 282 |
ratio is applied to the surface level to compute the new tracer field. |
|
| 283 |
|
A simplified version of the staggered in time, non-linear |
| 284 |
%------------------------------------------ |
free-surface algorithm is detailed hereafter, and can be compared |
| 285 |
\subsubsection{Free surface effect on the surface level thickness |
to the equivalent linear free-surface case (eq. \ref{eq:Gv-n-staggered} |
| 286 |
(Non-linear free surface): Momentum advection} |
to \ref{eq:t-n+1-staggered}) and can also be easily transposed |
| 287 |
|
to the synchronous time-stepping case. |
| 288 |
Regarding momentum advection, |
Among the simplifications, salinity equation, implicit operator |
| 289 |
the vector invariant formulation is similar to the |
and detailed elliptic equation are omitted. Surface forcing is |
| 290 |
advective form (as opposed to the flux form) and therefore |
explicitly written as fluxes of temperature, fresh water and |
| 291 |
does not need specific modification to include the |
momentum, $Q^{n+1/2}, P^{n+1/2}, F_{\bf v}^n$ respectively. |
| 292 |
free surface effect on the surface level thickness. |
$h^n$ and $dh^n$ are the column and grid box thickness in r-coordinate. |
| 293 |
Updating the {\bf hFacC,W,S} and the {\bf recip\_hFac}(s) |
%------------------------------------------------------------- |
| 294 |
at one given place (like describe before) is sufficient. |
\begin{eqnarray} |
| 295 |
|
\phi^{n}_{hyd} & = & \int b(\theta^{n},S^{n},r) dr |
| 296 |
With the flux form formulation, advection of momentum |
\label{eq:phi-hyd-nlfs} \\ |
| 297 |
can be treated exactly as the tracer advection is. |
\vec{\bf G}_{\vec{\bf v}}^{n-1/2}\hspace{-2mm} & = & |
| 298 |
Here the expansion/reduction factors ($hFacW^{n+1}/hFacW^n$ for $u$ |
\vec{\bf G}_{\vec{\bf v}} (dh^{n-1},\vec{\bf v}^{n-1/2}) |
| 299 |
and $hFacS^{n+1}/hFacS^n$ for $v$) are simply applied in the |
\hspace{+2mm};\hspace{+2mm} |
| 300 |
subroutine {\it TIMESTEP}. |
\vec{\bf G}_{\vec{\bf v}}^{(n)} = |
| 301 |
|
\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} |
| 303 |
|
\label{eq:Gv-n-nlfs} \\ |
| 304 |
|
%\vec{\bf G}_{\vec{\bf v}}^{(n)} & = & |
| 305 |
|
% \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} |
| 307 |
|
%\label{eq:Gv-n+5-nlfs} \\ |
| 308 |
|
%\vec{\bf v}^{*} & = & \vec{\bf v}^{n-1/2} + \frac{\Delta t}{dh^{n}} \left( |
| 309 |
|
%dh^{n-1}\vec{\bf G}_{\vec{\bf v}}^{(n)} + F_{\vec{\bf v}}^{n} \right) |
| 310 |
|
\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) |
| 312 |
|
- \Delta t \nabla \phi_{hyd}^{n} |
| 313 |
|
\label{eq:vstar-nlfs} |
| 314 |
|
\end{eqnarray} |
| 315 |
|
\hspace{3cm}$\longrightarrow$~~{\it update model~geometry~:~}${\bf hFac}(dh^n)$\\ |
| 316 |
|
\begin{eqnarray} |
| 317 |
|
\eta^{n+1/2} \hspace{-2mm} & = & |
| 318 |
|
\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 \\ |
| 320 |
|
& = & \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} |
| 322 |
|
\label{eq:nstar-nlfs} \\ |
| 323 |
|
\vec{\bf v}^{n+1/2}\hspace{-2mm} & = & |
| 324 |
|
\vec{\bf v}^{*} - g \Delta t \nabla \eta^{n+1/2} |
| 325 |
|
\label{eq:v-n+1-nlfs} \\ |
| 326 |
|
h^{n+1} & = & h^{n} + \Delta t P^{n+1/2} - \Delta t |
| 327 |
|
\nabla \cdot \int \vec{\bf v}^{n+1/2} dh^{n} |
| 328 |
|
\label{eq:h-n+1-nlfs} \\ |
| 329 |
|
G_{\theta}^{n} & = & G_{\theta} ( dh^{n}, u^{n+1/2}, \theta^{n} ) |
| 330 |
|
\hspace{+2mm};\hspace{+2mm} |
| 331 |
|
G_{\theta}^{(n+1/2)} = \frac{3}{2} G_{\theta}^{n} - \frac{1}{2} G_{\theta}^{n-1} |
| 332 |
|
\label{eq:Gt-n-nlfs} \\ |
| 333 |
|
%\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) |
| 335 |
|
\theta^{n+1} & = &\theta^{n} + \Delta t \frac{dh^n}{dh^{n+1}} \left( |
| 336 |
|
G_{\theta}^{(n+1/2)} |
| 337 |
|
+( P^{n+1/2} (\theta_{\mathrm{rain}}-\theta^n) + Q^{n+1/2})/dh^n \right) |
| 338 |
|
\nonumber \\ |
| 339 |
|
& & \label{eq:t-n+1-nlfs} |
| 340 |
|
\end{eqnarray} |
| 341 |
|
%------------------------------------------------------------- |
| 342 |
|
Two steps have been added to linear free-surface algorithm |
| 343 |
|
(eq. \ref{eq:Gv-n-staggered} to \ref{eq:t-n+1-staggered}): |
| 344 |
|
Firstly, the model ``geometry'' |
| 345 |
|
(here the {\bf hFacC,W,S}) is updated just before entering {\it |
| 346 |
|
SOLVE\_FOR\_PRESSURE}, using the current $dh^{n}$ field. |
| 347 |
|
Secondly, the vertically integrated continuity equation |
| 348 |
|
(eq.\ref{eq:h-n+1-nlfs}) has been added ({\bf exactConserv}{\em =TRUE}, |
| 349 |
|
in parameter file {\em data}, namelist {\em PARM01}) |
| 350 |
|
just before computing the vertical velocity, in subroutine |
| 351 |
|
{\em INTEGR\_CONTINUITY}. |
| 352 |
|
%This ensures that tracer and continuity equation discretization a |
| 353 |
|
Although this equation might appear redundant with eq.\ref{eq:nstar-nlfs}, |
| 354 |
|
the integrated column thickness $h^{n+1}$ will be different from |
| 355 |
|
$\eta^{n+1/2} + H$~ in the following cases: |
| 356 |
|
\begin{itemize} |
| 357 |
|
\item when Crank-Nickelson time-stepping is used (see section |
| 358 |
|
\ref{sec:freesurf-CrankNick}). |
| 359 |
|
\item when filters are applied to the flow field, after |
| 360 |
|
(\ref{eq:v-n+1-nlfs}) and alter the divergence of the flow. |
| 361 |
|
\item when the solver does not iterate until convergence~; |
| 362 |
|
for example, because a too large residual target was set |
| 363 |
|
({\bf cg2dTargetResidual}, parameter file {\em data}, namelist |
| 364 |
|
{\em PARM02}). |
| 365 |
|
\end{itemize}\noindent |
| 366 |
|
In this staggered time-stepping algorithm, the momentum tendencies |
| 367 |
|
are computed using $dh^{n-1}$ geometry factors. |
| 368 |
|
(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 |
| 370 |
|
\ref{sec:freesurf-tracer-advection}). |
| 371 |
|
The tracers are stepped forward later, using the recently updated |
| 372 |
|
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}); |
| 374 |
|
Then the tendencies are rescaled by $dh^n/dh^{n+1}$ to derive |
| 375 |
|
the new tracers values $(\theta,S)^{n+1}$ (eq.\ref{eq:t-n+1-nlfs}, |
| 376 |
|
in subroutine {\em CALC\_GT, CALC\_GS}). |
| 377 |
|
|
| 378 |
|
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 |
| 380 |
|
the fresh-water temperature $\theta_{\mathrm{rain}}$. |
| 381 |
|
|
| 382 |
|
Regarding the restart procedure, |
| 383 |
|
two 2.D fields $h^{n-1}$ and $(h^n-h^{n-1})/\Delta t$ |
| 384 |
|
in addition to the standard |
| 385 |
|
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}$) |
| 387 |
|
are stored in a "{\em pickup}" file. |
| 388 |
|
The model restarts reading this "{\em pickup}" file, |
| 389 |
|
then update the model geometry according to $h^{n-1}$, |
| 390 |
|
and compute $h^n$ and the vertical velocity |
| 391 |
|
%$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} |
| 393 |
|
to \ref{eq:t-n+1-nlfs}, {\em S/R FORWARD\_STEP}). |
| 394 |
|
\\ |
| 395 |
|
|
| 396 |
|
\fbox{ \begin{minipage}{4.75in} |
| 397 |
|
{\em INTEGR\_CONTINUITY} ({\em integr\_continuity.F}) |
| 398 |
|
|
| 399 |
|
$h^{n+1} -H_o$: {\bf etaH} ({\em DYNVARS.h}) |
| 400 |
|
|
| 401 |
|
$h^{n} -H_o$: {\bf etaHnm1} ({\em SURFACE.h}) |
| 402 |
|
|
| 403 |
|
$(h^{n+1}-h^{n})/\Delta t$: {\bf dEtaHdt} ({\em SURFACE.h}) |
| 404 |
|
|
| 405 |
|
\end{minipage} } |
| 406 |
|
|
| 407 |
|
\subsubsection{Non-linear free-surface and vertical resolution} |
| 408 |
|
\label{sec:nonlin-freesurf-dz_surf} |
| 409 |
|
|
| 410 |
|
When the amplitude of the free-surface variations becomes |
| 411 |
|
as large as the vertical resolution near the surface, |
| 412 |
|
the surface layer thickness can decrease to nearly zero or |
| 413 |
|
can even vanish completely. |
| 414 |
|
This later possibility has not been implemented, and a |
| 415 |
|
minimum relative thickness is imposed ({\bf hFacInf}, |
| 416 |
|
parameter file {\em data}, namelist {\em PARM01}) to prevent |
| 417 |
|
numerical instabilities caused by very thin surface level. |
| 418 |
|
|
| 419 |
|
A better alternative to the vanishing level problem has been |
| 420 |
|
found and implemented recently, and rely on a different |
| 421 |
|
vertical coordinate $r^*$~: |
| 422 |
|
The time variation ot the total column thickness becomes |
| 423 |
|
part of the r* coordinate motion, as in a $\sigma_{z},\sigma_{p}$ |
| 424 |
|
model, but the fixed part related to topography is treated |
| 425 |
|
as in a height or pressure coordinate model. |
| 426 |
|
A complete description is given in \cite{adcroft:04a}. |
| 427 |
|
|
| 428 |
|
The time-stepping implementation of the $r^*$ coordinate is |
| 429 |
|
identical to the non-linear free-surface in $r$ coordinate, |
| 430 |
|
and differences appear only in the spacial discretization. |
| 431 |
|
\marginpar{needs a subsection ref. here} |
| 432 |
|
|
Parent Directory
| 
Revision Log
| 
Revision Graph
| 
Patch