| 1 |
% $Header$ |
% $Header$ |
| 2 |
% $Name$ |
% $Name$ |
| 3 |
|
|
| 4 |
%\section{Time Integration} |
|
| 5 |
|
|
| 6 |
\subsection{Non-linear free surface} |
\subsection{Non-linear free surface} |
| 7 |
|
|
| 8 |
Recently, 2 options have added to the model |
Recently, two options have been added to the model (and have not yet |
| 9 |
(and therefore, have not yet been extensively tested) |
been extensively tested) that concern the free surface formulation. |
| 10 |
that concern the free surface formulation. |
|
| 11 |
|
|
|
%------------------------------------------ |
|
| 12 |
\subsubsection{Non-uniform linear-relation for the surface potential} |
\subsubsection{Non-uniform linear-relation for the surface potential} |
| 13 |
|
|
| 14 |
The linear relation between |
The linear relation between surface pressure/geo-potential |
| 15 |
surface pressure / geo- potential ($\Phi_{surf}$) |
($\Phi_{surf}$) and surface displacement ($\eta$) could be considered |
| 16 |
and surface displacement ($\eta$) |
to be a constant ($b_s=$ constant) |
| 17 |
has been considered as uniform ($b_{surf} = Constant$) |
\marginpar{add a reference to part.1 here} |
| 18 |
but is in fact |
but is in fact dependent on the position ($x,y,r$) |
| 19 |
dependent on the position |
since we linearize: |
| 20 |
since we have |
$$\Phi_{surf}=\int_{R_o}^{R_o+\eta} b dr \simeq b_s \eta |
| 21 |
$\Phi_{surf}=\int_{R_o}^{R_o+\eta} b dr \simeq b_{surf} \eta$ |
~\mathrm{with}~ b_s = b(\theta,S,r)_{r=R_o} |
| 22 |
|
\simeq b_s(\theta_{ref}(R_o),S_{ref}(R_o),R_o)$$ |
| 23 |
For the ocean, $b_{surf} = g \rho_{surf} / \rho_o \simeq g$ |
Note that, for convenience, the effect on $b_s$ of the local surface |
| 24 |
is a fairly good approximation since the relative difference |
$\theta,S$ are not considered here, but are incorporated in to |
| 25 |
in surface density are usually small. |
$\Phi'_{hyd}$. |
| 26 |
For the atmosphere, because of topographic effects, |
|
| 27 |
the reference surface pressure $R_o$ has large spacial differences |
For the ocean, $b_s = g \rho_{surf} / \rho_o \simeq g$ is a very good |
| 28 |
that are responsible for significant $b_surf$ variations |
approximation since the relative difference in surface density are |
| 29 |
(from 0.8 to 1.2 $[]$). |
usually small and only due to local $\theta,S$ gradients (because the |
| 30 |
|
upper surface, $R_o = 0$, is essentially flat). Therefore, they can |
| 31 |
Practically, in an oceanic configuration or |
easily be incorporated in $\Phi'_{hyd}$. |
| 32 |
when the default value (TRUE) of the parameter {\bf uniformLin\_PhiSurf} |
|
| 33 |
is used ) |
For the atmosphere, however, because of topographic effects, the |
| 34 |
then $b_{surf}$ is simply set to $g$ for the ocean |
reference surface pressure $R_o$ has large spacial variations that |
| 35 |
and $1.$ for the atmosphere.\\ |
are responsible for significant $b_s$ variations (from 0.8 to 1.2 |
| 36 |
Turning {\bf uniformLin\_PhiSurf} to "FALSE", allows to |
$[m^3/kg]$). For this reason, we use a non-uniform linear coefficient |
| 37 |
evaluate $b_{surf}$ from the reference Theta,Q vertical profile |
$b_s$. |
| 38 |
({\it S/R INI\_LINEAR\_PHISURF}) |
|
| 39 |
according to the reference surface pressure $P_o$ ($=R_o$): |
In practice, in an oceanic configuration or when the default value |
| 40 |
$b_{surf} = c_p \kappa (P_o / Pc)^{(\kappa - 1)} \theta_{ref}(P_o)$ |
(TRUE) of the parameter {\bf uniformLin\_PhiSurf} is used, then $b_s$ |
| 41 |
|
is simply set to $g$ for the ocean and $1.$ for the atmosphere. |
| 42 |
|
Turning {\bf uniformLin\_PhiSurf} to "FALSE", tells the code to |
| 43 |
|
evaluate $b_s$ from the reference vertical profile $\theta_{ref}$ |
| 44 |
|
({\it S/R INI\_LINEAR\_PHISURF}) according to the reference surface |
| 45 |
|
pressure $P_o$ ($=R_o$): $b_s = c_p \kappa (P_o / Pc)^{(\kappa - 1)} |
| 46 |
|
\theta_{ref}(P_o)$ |
| 47 |
|
|
| 48 |
|
|
|
%------------------------------------------ |
|
| 49 |
\subsubsection{Free surface effect on column total thickness |
\subsubsection{Free surface effect on column total thickness |
| 50 |
(Non-linear free surface)} |
(Non-linear free surface)} |
| 51 |
|
|
| 52 |
The total thickness of the fluid column is |
The total thickness of the fluid column is $r_{surf} - R_{fixed} = |
| 53 |
$r_{surf} - R_{min} = \eta + R_o - R_{min}$ |
\eta + R_o - R_{fixed}$ In the linear free surface approximation |
| 54 |
In the linear free surface approximation |
(detailed before), only the fixed part of it ($R_o - R_{fixed})$ is |
| 55 |
(detailed before), only the fixed part of |
considered when we integrate the continuity equation or compute tracer |
| 56 |
it ($R_o - R_{min})$ is considered when we integrate the |
and momentum advection term. |
| 57 |
continuity equation or compute tracer and momentum advection term. |
|
| 58 |
|
This approximation is dropped when using the non-linear free surface |
| 59 |
This approximation is dropped when using |
formulation. Here we discuss sections the barotropic part. In |
| 60 |
the non linear free surface formulation. |
sections \ref{sect:freesurf-tracer-advection} and |
| 61 |
Details follow here after for the barotropic part |
\ref{sect:freesurf-momentum-advection} we consider the baroclinic |
| 62 |
and in the 2 next subsection for the baroclinic |
component. |
| 63 |
part. |
|
| 64 |
|
|
| 65 |
%------------------------------------------ |
The continuous form of the model equations remains unchanged, except |
| 66 |
% Non Linear Barotropic part |
for the 2D continuity equation (\ref{eq-tCsC-eta}) which is now |
| 67 |
|
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-cont-2D}) that is now integrated |
|
|
from $R_{min}(x,y)$ up to $r_{surf}=R_o+\eta$ : |
|
| 68 |
|
|
| 69 |
\begin{displaymath} |
\begin{displaymath} |
| 70 |
\epsilon_{fs} \partial_t \eta = |
\epsilon_{fs} \partial_t \eta = |
| 71 |
\left. \dot{r} \right|_{r=r_{surf}} + \epsilon_{fw} (P-E) = |
\left. \dot{r} \right|_{r=r_{surf}} + \epsilon_{fw} (P-E) = |
| 72 |
- {\bf \nabla}_r \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 |
| 73 |
+ \epsilon_{fw} (P-E) |
+ \epsilon_{fw} (P-E) |
| 74 |
\end{displaymath} |
\end{displaymath} |
| 75 |
|
|
| 76 |
Since $\eta$ has a direct effect on the horizontal |
Since $\eta$ has a direct effect on the horizontal velocity (through |
| 77 |
velocity (through $\nabla_r \Phi_{surf}$), this |
$\nabla_h \Phi_{surf}$), this adds a non-linear term to the free |
| 78 |
add a non linear term to the free surface equation. |
surface equation. Several options for the time discretization of this |
| 79 |
|
non-linear part have been tested. |
| 80 |
Regarding the time discretization of this non linear part, |
|
| 81 |
several option are now being tested: |
If the column thickness is evaluated at time step $n$, and with |
| 82 |
|
implicit treatment of the surface potential gradient, equations |
| 83 |
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{solve_2d} \& \ref{solve_2d_rhs}) |
|
|
become: |
|
| 84 |
\begin{eqnarray*} |
\begin{eqnarray*} |
| 85 |
\epsilon_{fs} {\eta}^{n+1} - |
\epsilon_{fs} {\eta}^{n+1} - |
| 86 |
{\bf \nabla}_r \cdot \Delta t^2 (\eta^{n}+R_o-R_{min}) |
{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n}+R_o-R_{fixed}) |
| 87 |
{\bf \nabla}_r B_o {\eta}^{n+1} |
{\bf \nabla}_h b_s {\eta}^{n+1} |
| 88 |
= {\eta}^* |
= {\eta}^* |
|
%\label{solve_2d} |
|
| 89 |
\end{eqnarray*} |
\end{eqnarray*} |
| 90 |
where |
where |
| 91 |
\begin{eqnarray*} |
\begin{eqnarray*} |
| 92 |
{\eta}^* = \epsilon_{fs} \: {\eta}^{n} - |
{\eta}^* = \epsilon_{fs} \: {\eta}^{n} - |
| 93 |
\Delta t {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o+\eta} \vec{\bf v}^* dr |
\Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o+\eta^n} \vec{\bf v}^* dr |
| 94 |
\: + \: \epsilon_{fw} \Delta_t (P-E)^{n} |
\: + \: \epsilon_{fw} \Delta_t (P-E)^{n} |
|
%\label{solve_2d_rhs} |
|
| 95 |
\end{eqnarray*} |
\end{eqnarray*} |
| 96 |
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. |
| 97 |
|
|
| 98 |
Alternatively, the non-linear contribution can be evaluated fully |
Alternatively, the non-linear contribution can be evaluated fully |
| 99 |
explicitly: |
explicitly: |
| 100 |
\begin{eqnarray*} |
\begin{eqnarray*} |
| 101 |
\epsilon_{fs} {\eta}^{n+1} - |
\epsilon_{fs} {\eta}^{n+1} - |
| 102 |
{\bf \nabla}_r \cdot \Delta t^2 (R_o-R_{min}) |
{\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{fixed}) |
| 103 |
{\bf \nabla}_r B_o {\eta}^{n+1} |
{\bf \nabla}_h b_s {\eta}^{n+1} |
| 104 |
= {\eta}^* |
= {\eta}^* |
| 105 |
+{\bf \nabla}_r \cdot \Delta t^2 (\eta^{n}) |
+{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n}) |
| 106 |
{\bf \nabla}_r B_o {\eta}^{n} |
{\bf \nabla}_h b_s {\eta}^{n} |
| 107 |
\end{eqnarray*} |
\end{eqnarray*} |
| 108 |
This formulation allows to keep the initial solver matrix |
This formulation allows one to keep the initial solver matrix |
| 109 |
since the non-linear free surface only affects the RHS. |
unchanged though throughout the integration, since the non-linear free |
| 110 |
|
surface only affects the RHS. |
| 111 |
|
|
| 112 |
|
Finally, another option is a "linearized" formulation where the total |
| 113 |
|
column thickness appears only in the integral term of the RHS |
| 114 |
|
(\ref{eq-solve2D_rhs}) but not directly in the equation |
| 115 |
|
(\ref{eq-solve2D}). |
| 116 |
|
|
|
An other option is a "linearized" formulation where the |
|
|
total column thickness appears only in the integral term of |
|
|
the RHS (\ref{solve_2d_rhs}) but not directly in |
|
|
the equation (\ref{solve_2d}). |
|
| 117 |
|
|
|
%------------------------------------------ |
|
| 118 |
\subsubsection{Free surface effect on the surface level thickness |
\subsubsection{Free surface effect on the surface level thickness |
| 119 |
(Non-linear free surface): Tracer advection} |
(Non-linear free surface): Tracer advection} |
| 120 |
|
\label{sect:freesurf-tracer-advection} |
| 121 |
|
|
| 122 |
To ensure a global tracer conservation (i.e., the total amount) |
To ensure global tracer conservation (i.e., the total amount) as well |
| 123 |
as well as the local one (see tracer section for more details), |
as local conservation, the change in the surface level thickness must |
| 124 |
the change in the surface level thickness must be consistent with |
be consistent with the way the continuity equation is integrated, both |
| 125 |
the way the continuity equation is integrated, both in |
in the barotropic part (to find $\eta$) and baroclinic part (to find |
| 126 |
in the barotropic part (to find $\eta$) and baroclinic part |
$w = \dot{r}$). |
|
(to find $\dot{r}$). |
|
| 127 |
|
|
| 128 |
To illustrate this, let's consider the shallow water model, |
To illustrate this, consider the shallow water model, with uniform |
| 129 |
with uniform cartesian horizontal grid: |
Cartesian horizontal grid: |
| 130 |
$$ |
$$ |
| 131 |
\partial_t h + \nabla \cdot h \vec{\bf v} = 0 |
\partial_t h + \nabla \cdot h \vec{\bf v} = 0 |
| 132 |
$$ |
$$ |
| 135 |
$$ |
$$ |
| 136 |
\partial_t (h \theta) + \nabla \cdot ( h \theta \vec{\bf v})= 0 |
\partial_t (h \theta) + \nabla \cdot ( h \theta \vec{\bf v})= 0 |
| 137 |
$$ |
$$ |
| 138 |
Using the implicit (non-linear) free surface described before, we have: |
Using the implicit (non-linear) free surface described above (section |
| 139 |
|
\ref{sect:??}, we have: |
| 140 |
\begin{eqnarray*} |
\begin{eqnarray*} |
| 141 |
h^{n+1} = h^{n} - \Delta_t \nabla (h^n \, \vec{\bf v}^{n+1} ) \\ |
h^{n+1} = h^{n} - \Delta_t \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) \\ |
| 142 |
\end{eqnarray*} |
\end{eqnarray*} |
| 143 |
The discretized form of the tracer equation must use the same |
The discretized form of the tracer equation must adopt the same |
| 144 |
"geometry" to compute the tracer fluxes, that is, the same value of |
``form'' in the computation of tracer fluxes, that is, the same value |
| 145 |
h, as the continuity equation did: |
of $h$, as used in the continuity equation: |
| 146 |
\begin{eqnarray*} |
\begin{eqnarray*} |
| 147 |
h^{n+1} \, \theta^{n+1} = h^n \, \theta^n |
h^{n+1} \, \theta^{n+1} = h^n \, \theta^n |
| 148 |
- \Delta_t \nabla (h^n \, \theta^n \, \vec{\bf v}^{n+1}) |
- \Delta_t \nabla \cdot (h^n \, \theta^n \, \vec{\bf v}^{n+1}) |
| 149 |
\end{eqnarray*} |
\end{eqnarray*} |
| 150 |
|
|
| 151 |
In order to deal with the Adams-Bashforth time stepping, |
For Adams-Bashforth time-stepping, we implement this scheme slightly |
| 152 |
we implement this scheme slightly differently: |
differently from the linear free-surface method, using two steps: the |
| 153 |
the variation of the water column thickness (from |
variation of the water column thickness (from $h^n$ to $h^{n+1}$) is |
| 154 |
$h^n$ to $h^{n+1}$) |
not incorporated directly into the tracer equation. Instead, the |
| 155 |
is not incorporated directly to the tracer equation. |
model uses the $G_\theta$ terms (first step) as in the linear free |
| 156 |
Instead, |
surface formulation (with the "{\it surface correction}" turned "on", |
| 157 |
the model continues to evaluate the $G_\theta$ term |
see tracer section): |
|
exactly as it did with the Linear free surface formulation, |
|
|
with the "{\it surface correction}" turned "on" (see tracer section): |
|
| 158 |
$$ |
$$ |
| 159 |
G_\theta = \left(- \nabla (h^n \, \theta^n \, \vec{\bf v}^{n+1}) |
G_\theta^n = \left(- \nabla \cdot (h^n \, \theta^n \, \vec{\bf v}^{n+1}) |
| 160 |
+ w_{surf}^{n+1} \theta^n \right) / h^n |
- \dot{r}_{surf}^{n+1} \theta^n \right) / h^n |
| 161 |
$$ |
$$ |
| 162 |
Then after, |
Then, in a second step, the thickness variation (expansion/reduction) |
| 163 |
thickness variation (expansion/reduction) is taken into account : |
is taken into account: |
| 164 |
$$ |
$$ |
| 165 |
\theta^{n+1} = (\theta^n + \Delta_t G_\theta^{(n+1)}) \frac{h^{n+1}}{h^n} |
\theta^{n+1} = \theta^n + \Delta_t \frac{h^n}{h^{n+1}} G_\theta^{(n+1/2)} |
| 166 |
$$ |
$$ |
| 167 |
Note that with a simple forward time step (no Adams-Bashforth), |
Note that with a simple forward time step (no Adams-Bashforth), |
| 168 |
since |
since |
| 169 |
$ |
$ |
| 170 |
w_{surf} = - \nabla (h^n \, \vec{\bf v}^{n+1} ) = (h^{n+1} - h^{n}) |
\dot{r}_{surf}^{n+1} |
| 171 |
|
= - \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) = (h^{n+1} - h^{n}) |
| 172 |
/ \Delta_t |
/ \Delta_t |
| 173 |
$ |
$ |
| 174 |
those 2 formulations are equivalent. |
these two formulations are equivalent. |
| 175 |
|
|
| 176 |
|
Implementation in the MITgcm is as follows. The model ``geometry'' |
| 177 |
|
(here the {\bf hFacC,W,S}) is updated just before entering {\it |
| 178 |
|
SOLVE\_FOR\_PRESSURE}, using the current $\eta$ field. Then, at the |
| 179 |
|
end of the time step, the variables are advanced in time, so that |
| 180 |
|
$\eta^n$ replaces $\eta^{n-1}$. At the next time step, the tracer |
| 181 |
|
tendency ($G_\theta$) is computed using the same geometry, which is |
| 182 |
|
now consistent with $\eta^{n-1}$. Finally, in S/R {\it |
| 183 |
|
TIMESTEP\_TRACER}, the expansion/reduction ratio is applied to the |
| 184 |
|
surface level to compute the new tracer field. |
| 185 |
|
|
|
The implementation in the MITgcm follows this scheme. |
|
|
The model "geometry" (here the {\bf hFacC,W,S}) is updated |
|
|
at the beginning of {\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. |
|
| 186 |
|
|
|
%------------------------------------------ |
|
| 187 |
\subsubsection{Free surface effect on the surface level thickness |
\subsubsection{Free surface effect on the surface level thickness |
| 188 |
(Non-linear free surface): Momentum advection} |
(Non-linear free surface): Momentum advection} |
| 189 |
|
\label{sect:freesurf-momentum-advection} |
| 190 |
|
|
| 191 |
Regarding momentum advection, |
Regarding momentum advection, |
| 192 |
the vector invariant formulation is similar to the |
the vector invariant formulation is similar to the |
| 197 |
at one given place (like describe before) is sufficient. |
at one given place (like describe before) is sufficient. |
| 198 |
|
|
| 199 |
With the flux form formulation, advection of momentum |
With the flux form formulation, advection of momentum |
| 200 |
can be treated exactly as the tracer advection is, |
can be treated exactly as the tracer advection is. |
| 201 |
Here the expansion/reduction factors ($hFacW^{n+1}/hFacW^n$ for $u$ |
Here the expansion/reduction factors ($hFacW^{n+1}/hFacW^n$ for $u$ |
| 202 |
and $hFacS^{n+1}/hFacS^n$ for $v$) are simply applied in the |
and $hFacS^{n+1}/hFacS^n$ for $v$) are simply applied in the |
| 203 |
subroutine {\it TIMESTEP}. |
subroutine {\it TIMESTEP}. |
Parent Directory
| 
Revision Log
| 
Revision Graph
| 
Patch