1 | \documentclass[../main/NEMO_manual]{subfiles} |
---|
2 | |
---|
3 | \begin{document} |
---|
4 | \chapter{Ocean Dynamics (DYN)} |
---|
5 | \label{chap:DYN} |
---|
6 | |
---|
7 | \chaptertoc |
---|
8 | |
---|
9 | Using the representation described in \autoref{chap:DOM}, |
---|
10 | several semi-discrete space forms of the dynamical equations are available depending on |
---|
11 | the vertical coordinate used and on the conservation properties of the vorticity term. |
---|
12 | In all the equations presented here, the masking has been omitted for simplicity. |
---|
13 | One must be aware that all the quantities are masked fields and |
---|
14 | that each time an average or difference operator is used, the resulting field is multiplied by a mask. |
---|
15 | |
---|
16 | The prognostic ocean dynamics equation can be summarized as follows: |
---|
17 | \[ |
---|
18 | \text{NXT} = \dbinom {\text{VOR} + \text{KEG} + \text {ZAD} } |
---|
19 | {\text{COR} + \text{ADV} } |
---|
20 | + \text{HPG} + \text{SPG} + \text{LDF} + \text{ZDF} |
---|
21 | \] |
---|
22 | NXT stands for next, referring to the time-stepping. |
---|
23 | The first group of terms on the rhs of this equation corresponds to the Coriolis and advection terms that |
---|
24 | are decomposed into either a vorticity part (VOR), a kinetic energy part (KEG) and |
---|
25 | a vertical advection part (ZAD) in the vector invariant formulation, |
---|
26 | or a Coriolis and advection part (COR+ADV) in the flux formulation. |
---|
27 | The terms following these are the pressure gradient contributions |
---|
28 | (HPG, Hydrostatic Pressure Gradient, and SPG, Surface Pressure Gradient); |
---|
29 | and contributions from lateral diffusion (LDF) and vertical diffusion (ZDF), |
---|
30 | which are added to the rhs in the \mdl{dynldf} and \mdl{dynzdf} modules. |
---|
31 | The vertical diffusion term includes the surface and bottom stresses. |
---|
32 | The external forcings and parameterisations require complex inputs |
---|
33 | (surface wind stress calculation using bulk formulae, estimation of mixing coefficients) |
---|
34 | that are carried out in modules SBC, LDF and ZDF and are described in |
---|
35 | \autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively. |
---|
36 | |
---|
37 | In the present chapter we also describe the diagnostic equations used to compute the horizontal divergence, |
---|
38 | curl of the velocities (\emph{divcur} module) and the vertical velocity (\emph{wzvmod} module). |
---|
39 | |
---|
40 | The different options available to the user are managed by namelist variables. |
---|
41 | For term \textit{ttt} in the momentum equations, the logical namelist variables are \textit{ln\_dynttt\_xxx}, |
---|
42 | where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme. |
---|
43 | %If a CPP key is used for this term its name is \key{ttt}. |
---|
44 | The corresponding code can be found in the \textit{dynttt\_xxx} module in the DYN directory, |
---|
45 | and it is usually computed in the \textit{dyn\_ttt\_xxx} subroutine. |
---|
46 | |
---|
47 | The user has the option of extracting and outputting each tendency term from the 3D momentum equations |
---|
48 | (\texttt{trddyn?} defined), as described in \autoref{chap:MISC}. |
---|
49 | Furthermore, the tendency terms associated with the 2D barotropic vorticity balance (when \texttt{trdvor?} is defined) |
---|
50 | can be derived from the 3D terms. |
---|
51 | %%% |
---|
52 | \gmcomment{STEVEN: not quite sure I've got the sense of the last sentence. does |
---|
53 | MISC correspond to "extracting tendency terms" or "vorticity balance"?} |
---|
54 | |
---|
55 | \section{Sea surface height and diagnostic variables ($\eta$, $\zeta$, $\chi$, $w$)} |
---|
56 | \label{sec:DYN_divcur_wzv} |
---|
57 | |
---|
58 | %-------------------------------------------------------------------------------------------------------------- |
---|
59 | % Horizontal divergence and relative vorticity |
---|
60 | %-------------------------------------------------------------------------------------------------------------- |
---|
61 | \subsection[Horizontal divergence and relative vorticity (\textit{divcur.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{divcur})} |
---|
62 | \label{subsec:DYN_divcur} |
---|
63 | |
---|
64 | The vorticity is defined at an $f$-point (\ie\ corner point) as follows: |
---|
65 | \begin{equation} |
---|
66 | \label{eq:DYN_divcur_cur} |
---|
67 | \zeta =\frac{1}{e_{1f}\,e_{2f} }\left( {\;\delta_{i+1/2} \left[ {e_{2v}\;v} \right] |
---|
68 | -\delta_{j+1/2} \left[ {e_{1u}\;u} \right]\;} \right) |
---|
69 | \end{equation} |
---|
70 | |
---|
71 | The horizontal divergence is defined at a $T$-point. |
---|
72 | It is given by: |
---|
73 | \[ |
---|
74 | % \label{eq:DYN_divcur_div} |
---|
75 | \chi =\frac{1}{e_{1t}\,e_{2t}\,e_{3t} } |
---|
76 | \left( {\delta_i \left[ {e_{2u}\,e_{3u}\,u} \right] |
---|
77 | +\delta_j \left[ {e_{1v}\,e_{3v}\,v} \right]} \right) |
---|
78 | \] |
---|
79 | |
---|
80 | Note that although the vorticity has the same discrete expression in $z$- and $s$-coordinates, |
---|
81 | its physical meaning is not identical. |
---|
82 | $\zeta$ is a pseudo vorticity along $s$-surfaces |
---|
83 | (only pseudo because $(u,v)$ are still defined along geopotential surfaces, |
---|
84 | but are not necessarily defined at the same depth). |
---|
85 | |
---|
86 | The vorticity and divergence at the \textit{before} step are used in the computation of |
---|
87 | the horizontal diffusion of momentum. |
---|
88 | Note that because they have been calculated prior to the Asselin filtering of the \textit{before} velocities, |
---|
89 | the \textit{before} vorticity and divergence arrays must be included in the restart file to |
---|
90 | ensure perfect restartability. |
---|
91 | The vorticity and divergence at the \textit{now} time step are used for the computation of |
---|
92 | the nonlinear advection and of the vertical velocity respectively. |
---|
93 | |
---|
94 | %-------------------------------------------------------------------------------------------------------------- |
---|
95 | % Sea Surface Height evolution |
---|
96 | %-------------------------------------------------------------------------------------------------------------- |
---|
97 | \subsection[Horizontal divergence and relative vorticity (\textit{sshwzv.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})} |
---|
98 | \label{subsec:DYN_sshwzv} |
---|
99 | |
---|
100 | The sea surface height is given by: |
---|
101 | \begin{equation} |
---|
102 | \label{eq:DYN_spg_ssh} |
---|
103 | \begin{aligned} |
---|
104 | \frac{\partial \eta }{\partial t} |
---|
105 | &\equiv \frac{1}{e_{1t} e_{2t} }\sum\limits_k { \left\{ \delta_i \left[ {e_{2u}\,e_{3u}\;u} \right] |
---|
106 | +\delta_j \left[ {e_{1v}\,e_{3v}\;v} \right] \right\} } |
---|
107 | - \frac{\textit{emp}}{\rho_w } \\ |
---|
108 | &\equiv \sum\limits_k {\chi \ e_{3t}} - \frac{\textit{emp}}{\rho_w } |
---|
109 | \end{aligned} |
---|
110 | \end{equation} |
---|
111 | where \textit{emp} is the surface freshwater budget (evaporation minus precipitation), |
---|
112 | expressed in Kg/m$^2$/s (which is equal to mm/s), |
---|
113 | and $\rho_w$=1,035~Kg/m$^3$ is the reference density of sea water (Boussinesq approximation). |
---|
114 | If river runoff is expressed as a surface freshwater flux (see \autoref{chap:SBC}) then |
---|
115 | \textit{emp} can be written as the evaporation minus precipitation, minus the river runoff. |
---|
116 | The sea-surface height is evaluated using exactly the same time stepping scheme as |
---|
117 | the tracer equation \autoref{eq:TRA_nxt}: |
---|
118 | a leapfrog scheme in combination with an Asselin time filter, |
---|
119 | \ie\ the velocity appearing in \autoref{eq:DYN_spg_ssh} is centred in time (\textit{now} velocity). |
---|
120 | This is of paramount importance. |
---|
121 | Replacing $T$ by the number $1$ in the tracer equation and summing over the water column must lead to |
---|
122 | the sea surface height equation otherwise tracer content will not be conserved |
---|
123 | \citep{griffies.pacanowski.ea_MWR01, leclair.madec_OM09}. |
---|
124 | |
---|
125 | The vertical velocity is computed by an upward integration of the horizontal divergence starting at the bottom, |
---|
126 | taking into account the change of the thickness of the levels: |
---|
127 | \begin{equation} |
---|
128 | \label{eq:DYN_wzv} |
---|
129 | \left\{ |
---|
130 | \begin{aligned} |
---|
131 | &\left. w \right|_{k_b-1/2} \quad= 0 \qquad \text{where } k_b \text{ is the level just above the sea floor } \\ |
---|
132 | &\left. w \right|_{k+1/2} = \left. w \right|_{k-1/2} + \left. e_{3t} \right|_{k}\; \left. \chi \right|_k |
---|
133 | - \frac{1} {2 \rdt} \left( \left. e_{3t}^{t+1}\right|_{k} - \left. e_{3t}^{t-1}\right|_{k}\right) |
---|
134 | \end{aligned} |
---|
135 | \right. |
---|
136 | \end{equation} |
---|
137 | |
---|
138 | In the case of a non-linear free surface (\texttt{vvl?}), the top vertical velocity is $-\textit{emp}/\rho_w$, |
---|
139 | as changes in the divergence of the barotropic transport are absorbed into the change of the level thicknesses, |
---|
140 | re-orientated downward. |
---|
141 | \gmcomment{not sure of this... to be modified with the change in emp setting} |
---|
142 | In the case of a linear free surface, the time derivative in \autoref{eq:DYN_wzv} disappears. |
---|
143 | The upper boundary condition applies at a fixed level $z=0$. |
---|
144 | The top vertical velocity is thus equal to the divergence of the barotropic transport |
---|
145 | (\ie\ the first term in the right-hand-side of \autoref{eq:DYN_spg_ssh}). |
---|
146 | |
---|
147 | Note also that whereas the vertical velocity has the same discrete expression in $z$- and $s$-coordinates, |
---|
148 | its physical meaning is not the same: |
---|
149 | in the second case, $w$ is the velocity normal to the $s$-surfaces. |
---|
150 | Note also that the $k$-axis is re-orientated downwards in the \fortran\ code compared to |
---|
151 | the indexing used in the semi-discrete equations such as \autoref{eq:DYN_wzv} |
---|
152 | (see \autoref{subsec:DOM_Num_Index_vertical}). |
---|
153 | |
---|
154 | \section{Coriolis and advection: vector invariant form} |
---|
155 | \label{sec:DYN_adv_cor_vect} |
---|
156 | %-----------------------------------------nam_dynadv---------------------------------------------------- |
---|
157 | |
---|
158 | \begin{listing} |
---|
159 | \nlst{namdyn_adv} |
---|
160 | \caption{\forcode{&namdyn_adv}} |
---|
161 | \label{lst:namdyn_adv} |
---|
162 | \end{listing} |
---|
163 | %------------------------------------------------------------------------------------------------------------- |
---|
164 | |
---|
165 | The vector invariant form of the momentum equations is the one most often used in |
---|
166 | applications of the \NEMO\ ocean model. |
---|
167 | The flux form option (see next section) has been present since version $2$. |
---|
168 | Options are defined through the \nam{dyn_adv}{dyn\_adv} namelist variables Coriolis and |
---|
169 | momentum advection terms are evaluated using a leapfrog scheme, |
---|
170 | \ie\ the velocity appearing in these expressions is centred in time (\textit{now} velocity). |
---|
171 | At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied following |
---|
172 | \autoref{chap:LBC}. |
---|
173 | |
---|
174 | \subsection[Vorticity term (\textit{dynvor.F90})]{Vorticity term (\protect\mdl{dynvor})} |
---|
175 | \label{subsec:DYN_vor} |
---|
176 | %------------------------------------------nam_dynvor---------------------------------------------------- |
---|
177 | |
---|
178 | \begin{listing} |
---|
179 | \nlst{namdyn_vor} |
---|
180 | \caption{\forcode{&namdyn_vor}} |
---|
181 | \label{lst:namdyn_vor} |
---|
182 | \end{listing} |
---|
183 | %------------------------------------------------------------------------------------------------------------- |
---|
184 | |
---|
185 | Options are defined through the \nam{dyn_vor}{dyn\_vor} namelist variables. |
---|
186 | Four discretisations of the vorticity term (\texttt{ln\_dynvor\_xxx}\forcode{=.true.}) are available: |
---|
187 | conserving potential enstrophy of horizontally non-divergent flow (ENS scheme); |
---|
188 | conserving horizontal kinetic energy (ENE scheme); |
---|
189 | conserving potential enstrophy for the relative vorticity term and |
---|
190 | horizontal kinetic energy for the planetary vorticity term (MIX scheme); |
---|
191 | or conserving both the potential enstrophy of horizontally non-divergent flow and horizontal kinetic energy |
---|
192 | (EEN scheme) (see \autoref{subsec:INVARIANTS_vorEEN}). |
---|
193 | In the case of ENS, ENE or MIX schemes the land sea mask may be slightly modified to ensure the consistency of |
---|
194 | vorticity term with analytical equations (\np[=.true.]{ln_dynvor_con}{ln\_dynvor\_con}). |
---|
195 | The vorticity terms are all computed in dedicated routines that can be found in the \mdl{dynvor} module. |
---|
196 | |
---|
197 | %------------------------------------------------------------- |
---|
198 | % enstrophy conserving scheme |
---|
199 | %------------------------------------------------------------- |
---|
200 | \subsubsection[Enstrophy conserving scheme (\forcode{ln_dynvor_ens})]{Enstrophy conserving scheme (\protect\np{ln_dynvor_ens}{ln\_dynvor\_ens})} |
---|
201 | \label{subsec:DYN_vor_ens} |
---|
202 | |
---|
203 | In the enstrophy conserving case (ENS scheme), |
---|
204 | the discrete formulation of the vorticity term provides a global conservation of the enstrophy |
---|
205 | ($ [ (\zeta +f ) / e_{3f} ]^2 $ in $s$-coordinates) for a horizontally non-divergent flow (\ie\ $\chi$=$0$), |
---|
206 | but does not conserve the total kinetic energy. |
---|
207 | It is given by: |
---|
208 | \begin{equation} |
---|
209 | \label{eq:DYN_vor_ens} |
---|
210 | \left\{ |
---|
211 | \begin{aligned} |
---|
212 | {+\frac{1}{e_{1u} } } & {\overline {\left( { \frac{\zeta +f}{e_{3f} }} \right)} }^{\,i} |
---|
213 | & {\overline{\overline {\left( {e_{1v}\,e_{3v}\;v} \right)}} }^{\,i, j+1/2} \\ |
---|
214 | {- \frac{1}{e_{2v} } } & {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)} }^{\,j} |
---|
215 | & {\overline{\overline {\left( {e_{2u}\,e_{3u}\;u} \right)}} }^{\,i+1/2, j} |
---|
216 | \end{aligned} |
---|
217 | \right. |
---|
218 | \end{equation} |
---|
219 | |
---|
220 | %------------------------------------------------------------- |
---|
221 | % energy conserving scheme |
---|
222 | %------------------------------------------------------------- |
---|
223 | \subsubsection[Energy conserving scheme (\forcode{ln_dynvor_ene})]{Energy conserving scheme (\protect\np{ln_dynvor_ene}{ln\_dynvor\_ene})} |
---|
224 | \label{subsec:DYN_vor_ene} |
---|
225 | |
---|
226 | The kinetic energy conserving scheme (ENE scheme) conserves the global kinetic energy but not the global enstrophy. |
---|
227 | It is given by: |
---|
228 | \begin{equation} |
---|
229 | \label{eq:DYN_vor_ene} |
---|
230 | \left\{ |
---|
231 | \begin{aligned} |
---|
232 | {+\frac{1}{e_{1u}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right) |
---|
233 | \; \overline {\left( {e_{1v}\,e_{3v}\;v} \right)} ^{\,i+1/2}} }^{\,j} } \\ |
---|
234 | {- \frac{1}{e_{2v}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right) |
---|
235 | \; \overline {\left( {e_{2u}\,e_{3u}\;u} \right)} ^{\,j+1/2}} }^{\,i} } |
---|
236 | \end{aligned} |
---|
237 | \right. |
---|
238 | \end{equation} |
---|
239 | |
---|
240 | %------------------------------------------------------------- |
---|
241 | % mix energy/enstrophy conserving scheme |
---|
242 | %------------------------------------------------------------- |
---|
243 | \subsubsection[Mixed energy/enstrophy conserving scheme (\forcode{ln_dynvor_mix})]{Mixed energy/enstrophy conserving scheme (\protect\np{ln_dynvor_mix}{ln\_dynvor\_mix})} |
---|
244 | \label{subsec:DYN_vor_mix} |
---|
245 | |
---|
246 | For the mixed energy/enstrophy conserving scheme (MIX scheme), a mixture of the two previous schemes is used. |
---|
247 | It consists of the ENS scheme (\autoref{eq:DYN_vor_ens}) for the relative vorticity term, |
---|
248 | and of the ENE scheme (\autoref{eq:DYN_vor_ene}) applied to the planetary vorticity term. |
---|
249 | \[ |
---|
250 | % \label{eq:DYN_vor_mix} |
---|
251 | \left\{ { |
---|
252 | \begin{aligned} |
---|
253 | {+\frac{1}{e_{1u} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^{\,i} |
---|
254 | \; {\overline{\overline {\left( {e_{1v}\,e_{3v}\;v} \right)}} }^{\,i,j+1/2} -\frac{1}{e_{1u} } |
---|
255 | \; {\overline {\left( {\frac{f}{e_{3f} }} \right) |
---|
256 | \;\overline {\left( {e_{1v}\,e_{3v}\;v} \right)} ^{\,i+1/2}} }^{\,j} } \\ |
---|
257 | {-\frac{1}{e_{2v} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^j |
---|
258 | \; {\overline{\overline {\left( {e_{2u}\,e_{3u}\;u} \right)}} }^{\,i+1/2,j} +\frac{1}{e_{2v} } |
---|
259 | \; {\overline {\left( {\frac{f}{e_{3f} }} \right) |
---|
260 | \;\overline {\left( {e_{2u}\,e_{3u}\;u} \right)} ^{\,j+1/2}} }^{\,i} } \hfill |
---|
261 | \end{aligned} |
---|
262 | } \right. |
---|
263 | \] |
---|
264 | |
---|
265 | %------------------------------------------------------------- |
---|
266 | % energy and enstrophy conserving scheme |
---|
267 | %------------------------------------------------------------- |
---|
268 | \subsubsection[Energy and enstrophy conserving scheme (\forcode{ln_dynvor_een})]{Energy and enstrophy conserving scheme (\protect\np{ln_dynvor_een}{ln\_dynvor\_een})} |
---|
269 | \label{subsec:DYN_vor_een} |
---|
270 | |
---|
271 | In both the ENS and ENE schemes, |
---|
272 | it is apparent that the combination of $i$ and $j$ averages of the velocity allows for |
---|
273 | the presence of grid point oscillation structures that will be invisible to the operator. |
---|
274 | These structures are \textit{computational modes} that will be at least partly damped by |
---|
275 | the momentum diffusion operator (\ie\ the subgrid-scale advection), but not by the resolved advection term. |
---|
276 | The ENS and ENE schemes therefore do not contribute to dump any grid point noise in the horizontal velocity field. |
---|
277 | Such noise would result in more noise in the vertical velocity field, an undesirable feature. |
---|
278 | This is a well-known characteristic of $C$-grid discretization where |
---|
279 | $u$ and $v$ are located at different grid points, |
---|
280 | a price worth paying to avoid a double averaging in the pressure gradient term as in the $B$-grid. |
---|
281 | \gmcomment{ To circumvent this, Adcroft (ADD REF HERE) |
---|
282 | Nevertheless, this technique strongly distort the phase and group velocity of Rossby waves....} |
---|
283 | |
---|
284 | A very nice solution to the problem of double averaging was proposed by \citet{arakawa.hsu_MWR90}. |
---|
285 | The idea is to get rid of the double averaging by considering triad combinations of vorticity. |
---|
286 | It is noteworthy that this solution is conceptually quite similar to the one proposed by |
---|
287 | \citep{griffies.gnanadesikan.ea_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:INVARIANTS}). |
---|
288 | |
---|
289 | The \citet{arakawa.hsu_MWR90} vorticity advection scheme for a single layer is modified |
---|
290 | for spherical coordinates as described by \citet{arakawa.lamb_MWR81} to obtain the EEN scheme. |
---|
291 | First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point: |
---|
292 | \[ |
---|
293 | % \label{eq:DYN_pot_vor} |
---|
294 | q = \frac{\zeta +f} {e_{3f} } |
---|
295 | \] |
---|
296 | where the relative vorticity is defined by (\autoref{eq:DYN_divcur_cur}), |
---|
297 | the Coriolis parameter is given by $f=2 \,\Omega \;\sin \varphi _f $ and the layer thickness at $f$-points is: |
---|
298 | \begin{equation} |
---|
299 | \label{eq:DYN_een_e3f} |
---|
300 | e_{3f} = \overline{\overline {e_{3t} }} ^{\,i+1/2,j+1/2} |
---|
301 | \end{equation} |
---|
302 | |
---|
303 | \begin{figure}[!ht] |
---|
304 | \centering |
---|
305 | \includegraphics[width=0.66\textwidth]{Fig_DYN_een_triad} |
---|
306 | \caption[Triads used in the energy and enstrophy conserving scheme (EEN)]{ |
---|
307 | Triads used in the energy and enstrophy conserving scheme (EEN) for |
---|
308 | $u$-component (upper panel) and $v$-component (lower panel).} |
---|
309 | \label{fig:DYN_een_triad} |
---|
310 | \end{figure} |
---|
311 | % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
312 | |
---|
313 | A key point in \autoref{eq:DYN_een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made. |
---|
314 | It uses the sum of masked t-point vertical scale factor divided either by the sum of the four t-point masks |
---|
315 | (\np[=1]{nn_een_e3f}{nn\_een\_e3f}), or just by $4$ (\np[=.true.]{nn_een_e3f}{nn\_een\_e3f}). |
---|
316 | The latter case preserves the continuity of $e_{3f}$ when one or more of the neighbouring $e_{3t}$ tends to zero and |
---|
317 | extends by continuity the value of $e_{3f}$ into the land areas. |
---|
318 | This case introduces a sub-grid-scale topography at f-points |
---|
319 | (with a systematic reduction of $e_{3f}$ when a model level intercept the bathymetry) |
---|
320 | that tends to reinforce the topostrophy of the flow |
---|
321 | (\ie\ the tendency of the flow to follow the isobaths) \citep{penduff.le-sommer.ea_OS07}. |
---|
322 | |
---|
323 | Next, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as |
---|
324 | the following triad combinations of the neighbouring potential vorticities defined at f-points |
---|
325 | (\autoref{fig:DYN_een_triad}): |
---|
326 | \begin{equation} |
---|
327 | \label{eq:DYN_Q_triads} |
---|
328 | _i^j \mathbb{Q}^{i_p}_{j_p} |
---|
329 | = \frac{1}{12} \ \left( q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p} \right) |
---|
330 | \end{equation} |
---|
331 | where the indices $i_p$ and $k_p$ take the values: $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$. |
---|
332 | |
---|
333 | Finally, the vorticity terms are represented as: |
---|
334 | \begin{equation} |
---|
335 | \label{eq:DYN_vor_een} |
---|
336 | \left\{ { |
---|
337 | \begin{aligned} |
---|
338 | +q\,e_3 \, v &\equiv +\frac{1}{e_{1u} } \sum_{\substack{i_p,\,k_p}} |
---|
339 | {^{i+1/2-i_p}_j} \mathbb{Q}^{i_p}_{j_p} \left( e_{1v}\,e_{3v} \;v \right)^{i+1/2-i_p}_{j+j_p} \\ |
---|
340 | - q\,e_3 \, u &\equiv -\frac{1}{e_{2v} } \sum_{\substack{i_p,\,k_p}} |
---|
341 | {^i_{j+1/2-j_p}} \mathbb{Q}^{i_p}_{j_p} \left( e_{2u}\,e_{3u} \;u \right)^{i+i_p}_{j+1/2-j_p} \\ |
---|
342 | \end{aligned} |
---|
343 | } \right. |
---|
344 | \end{equation} |
---|
345 | |
---|
346 | This EEN scheme in fact combines the conservation properties of the ENS and ENE schemes. |
---|
347 | It conserves both total energy and potential enstrophy in the limit of horizontally nondivergent flow |
---|
348 | (\ie\ $\chi$=$0$) (see \autoref{subsec:INVARIANTS_vorEEN}). |
---|
349 | Applied to a realistic ocean configuration, it has been shown that it leads to a significant reduction of |
---|
350 | the noise in the vertical velocity field \citep{le-sommer.penduff.ea_OM09}. |
---|
351 | Furthermore, used in combination with a partial steps representation of bottom topography, |
---|
352 | it improves the interaction between current and topography, |
---|
353 | leading to a larger topostrophy of the flow \citep{barnier.madec.ea_OD06, penduff.le-sommer.ea_OS07}. |
---|
354 | |
---|
355 | %-------------------------------------------------------------------------------------------------------------- |
---|
356 | % Kinetic Energy Gradient term |
---|
357 | %-------------------------------------------------------------------------------------------------------------- |
---|
358 | \subsection[Kinetic energy gradient term (\textit{dynkeg.F90})]{Kinetic energy gradient term (\protect\mdl{dynkeg})} |
---|
359 | \label{subsec:DYN_keg} |
---|
360 | |
---|
361 | As demonstrated in \autoref{apdx:INVARIANTS}, |
---|
362 | there is a single discrete formulation of the kinetic energy gradient term that, |
---|
363 | together with the formulation chosen for the vertical advection (see below), |
---|
364 | conserves the total kinetic energy: |
---|
365 | \[ |
---|
366 | % \label{eq:DYN_keg} |
---|
367 | \left\{ |
---|
368 | \begin{aligned} |
---|
369 | -\frac{1}{2 \; e_{1u} } & \ \delta_{i+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right] \\ |
---|
370 | -\frac{1}{2 \; e_{2v} } & \ \delta_{j+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right] |
---|
371 | \end{aligned} |
---|
372 | \right. |
---|
373 | \] |
---|
374 | |
---|
375 | %-------------------------------------------------------------------------------------------------------------- |
---|
376 | % Vertical advection term |
---|
377 | %-------------------------------------------------------------------------------------------------------------- |
---|
378 | \subsection[Vertical advection term (\textit{dynzad.F90})]{Vertical advection term (\protect\mdl{dynzad})} |
---|
379 | \label{subsec:DYN_zad} |
---|
380 | |
---|
381 | The discrete formulation of the vertical advection, t |
---|
382 | ogether with the formulation chosen for the gradient of kinetic energy (KE) term, |
---|
383 | conserves the total kinetic energy. |
---|
384 | Indeed, the change of KE due to the vertical advection is exactly balanced by |
---|
385 | the change of KE due to the gradient of KE (see \autoref{apdx:INVARIANTS}). |
---|
386 | \[ |
---|
387 | % \label{eq:DYN_zad} |
---|
388 | \left\{ |
---|
389 | \begin{aligned} |
---|
390 | -\frac{1} {e_{1u}\,e_{2u}\,e_{3u}} &\ \overline{\ \overline{ e_{1t}\,e_{2t}\;w } ^{\,i+1/2} \;\delta_{k+1/2} \left[ u \right]\ }^{\,k} \\ |
---|
391 | -\frac{1} {e_{1v}\,e_{2v}\,e_{3v}} &\ \overline{\ \overline{ e_{1t}\,e_{2t}\;w } ^{\,j+1/2} \;\delta_{k+1/2} \left[ u \right]\ }^{\,k} |
---|
392 | \end{aligned} |
---|
393 | \right. |
---|
394 | \] |
---|
395 | When \np[=.true.]{ln_dynzad_zts}{ln\_dynzad\_zts}, |
---|
396 | a split-explicit time stepping with 5 sub-timesteps is used on the vertical advection term. |
---|
397 | This option can be useful when the value of the timestep is limited by vertical advection \citep{lemarie.debreu.ea_OM15}. |
---|
398 | Note that in this case, |
---|
399 | a similar split-explicit time stepping should be used on vertical advection of tracer to ensure a better stability, |
---|
400 | an option which is only available with a TVD scheme (see \np{ln_traadv_tvd_zts}{ln\_traadv\_tvd\_zts} in \autoref{subsec:TRA_adv_tvd}). |
---|
401 | |
---|
402 | \section{Coriolis and advection: flux form} |
---|
403 | \label{sec:DYN_adv_cor_flux} |
---|
404 | %------------------------------------------nam_dynadv---------------------------------------------------- |
---|
405 | |
---|
406 | %------------------------------------------------------------------------------------------------------------- |
---|
407 | |
---|
408 | Options are defined through the \nam{dyn_adv}{dyn\_adv} namelist variables. |
---|
409 | In the flux form (as in the vector invariant form), |
---|
410 | the Coriolis and momentum advection terms are evaluated using a leapfrog scheme, |
---|
411 | \ie\ the velocity appearing in their expressions is centred in time (\textit{now} velocity). |
---|
412 | At the lateral boundaries either free slip, |
---|
413 | no slip or partial slip boundary conditions are applied following \autoref{chap:LBC}. |
---|
414 | |
---|
415 | %-------------------------------------------------------------------------------------------------------------- |
---|
416 | % Coriolis plus curvature metric terms |
---|
417 | %-------------------------------------------------------------------------------------------------------------- |
---|
418 | \subsection[Coriolis plus curvature metric terms (\textit{dynvor.F90})]{Coriolis plus curvature metric terms (\protect\mdl{dynvor})} |
---|
419 | \label{subsec:DYN_cor_flux} |
---|
420 | |
---|
421 | In flux form, the vorticity term reduces to a Coriolis term in which the Coriolis parameter has been modified to account for the "metric" term. |
---|
422 | This altered Coriolis parameter is thus discretised at $f$-points. |
---|
423 | It is given by: |
---|
424 | \begin{multline*} |
---|
425 | % \label{eq:DYN_cor_metric} |
---|
426 | f+\frac{1}{e_1 e_2 }\left( {v\frac{\partial e_2 }{\partial i} - u\frac{\partial e_1 }{\partial j}} \right) \\ |
---|
427 | \equiv f + \frac{1}{e_{1f} e_{2f} } \left( { \ \overline v ^{i+1/2}\delta_{i+1/2} \left[ {e_{2u} } \right] |
---|
428 | - \overline u ^{j+1/2}\delta_{j+1/2} \left[ {e_{1u} } \right] } \ \right) |
---|
429 | \end{multline*} |
---|
430 | |
---|
431 | Any of the (\autoref{eq:DYN_vor_ens}), (\autoref{eq:DYN_vor_ene}) and (\autoref{eq:DYN_vor_een}) schemes can be used to |
---|
432 | compute the product of the Coriolis parameter and the vorticity. |
---|
433 | However, the energy-conserving scheme (\autoref{eq:DYN_vor_een}) has exclusively been used to date. |
---|
434 | This term is evaluated using a leapfrog scheme, \ie\ the velocity is centred in time (\textit{now} velocity). |
---|
435 | |
---|
436 | %-------------------------------------------------------------------------------------------------------------- |
---|
437 | % Flux form Advection term |
---|
438 | %-------------------------------------------------------------------------------------------------------------- |
---|
439 | \subsection[Flux form advection term (\textit{dynadv.F90})]{Flux form advection term (\protect\mdl{dynadv})} |
---|
440 | \label{subsec:DYN_adv_flux} |
---|
441 | |
---|
442 | The discrete expression of the advection term is given by: |
---|
443 | \[ |
---|
444 | % \label{eq:DYN_adv} |
---|
445 | \left\{ |
---|
446 | \begin{aligned} |
---|
447 | \frac{1}{e_{1u}\,e_{2u}\,e_{3u}} |
---|
448 | \left( \delta_{i+1/2} \left[ \overline{e_{2u}\,e_{3u}\;u }^{i } \ u_t \right] |
---|
449 | + \delta_{j } \left[ \overline{e_{1u}\,e_{3u}\;v }^{i+1/2} \ u_f \right] \right. \ \; \\ |
---|
450 | \left. + \delta_{k } \left[ \overline{e_{1w}\,e_{2w}\;w}^{i+1/2} \ u_{uw} \right] \right) \\ |
---|
451 | \\ |
---|
452 | \frac{1}{e_{1v}\,e_{2v}\,e_{3v}} |
---|
453 | \left( \delta_{i } \left[ \overline{e_{2u}\,e_{3u }\;u }^{j+1/2} \ v_f \right] |
---|
454 | + \delta_{j+1/2} \left[ \overline{e_{1u}\,e_{3u }\;v }^{i } \ v_t \right] \right. \ \, \, \\ |
---|
455 | \left. + \delta_{k } \left[ \overline{e_{1w}\,e_{2w}\;w}^{j+1/2} \ v_{vw} \right] \right) \\ |
---|
456 | \end{aligned} |
---|
457 | \right. |
---|
458 | \] |
---|
459 | |
---|
460 | Two advection schemes are available: |
---|
461 | a $2^{nd}$ order centered finite difference scheme, CEN2, |
---|
462 | or a $3^{rd}$ order upstream biased scheme, UBS. |
---|
463 | The latter is described in \citet{shchepetkin.mcwilliams_OM05}. |
---|
464 | The schemes are selected using the namelist logicals \np{ln_dynadv_cen2}{ln\_dynadv\_cen2} and \np{ln_dynadv_ubs}{ln\_dynadv\_ubs}. |
---|
465 | In flux form, the schemes differ by the choice of a space and time interpolation to define the value of |
---|
466 | $u$ and $v$ at the centre of each face of $u$- and $v$-cells, \ie\ at the $T$-, $f$-, |
---|
467 | and $uw$-points for $u$ and at the $f$-, $T$- and $vw$-points for $v$. |
---|
468 | |
---|
469 | %------------------------------------------------------------- |
---|
470 | % 2nd order centred scheme |
---|
471 | %------------------------------------------------------------- |
---|
472 | \subsubsection[CEN2: $2^{nd}$ order centred scheme (\forcode{ln_dynadv_cen2})]{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln_dynadv_cen2}{ln\_dynadv\_cen2})} |
---|
473 | \label{subsec:DYN_adv_cen2} |
---|
474 | |
---|
475 | In the centered $2^{nd}$ order formulation, the velocity is evaluated as the mean of the two neighbouring points: |
---|
476 | \begin{equation} |
---|
477 | \label{eq:DYN_adv_cen2} |
---|
478 | \left\{ |
---|
479 | \begin{aligned} |
---|
480 | u_T^{cen2} &=\overline u^{i } \quad & u_F^{cen2} &=\overline u^{j+1/2} \quad & u_{uw}^{cen2} &=\overline u^{k+1/2} \\ |
---|
481 | v_F^{cen2} &=\overline v ^{i+1/2} \quad & v_F^{cen2} &=\overline v^j \quad & v_{vw}^{cen2} &=\overline v ^{k+1/2} \\ |
---|
482 | \end{aligned} |
---|
483 | \right. |
---|
484 | \end{equation} |
---|
485 | |
---|
486 | The scheme is non diffusive (\ie\ conserves the kinetic energy) but dispersive (\ie\ it may create false extrema). |
---|
487 | It is therefore notoriously noisy and must be used in conjunction with an explicit diffusion operator to |
---|
488 | produce a sensible solution. |
---|
489 | The associated time-stepping is performed using a leapfrog scheme in conjunction with an Asselin time-filter, |
---|
490 | so $u$ and $v$ are the \emph{now} velocities. |
---|
491 | |
---|
492 | %------------------------------------------------------------- |
---|
493 | % UBS scheme |
---|
494 | %------------------------------------------------------------- |
---|
495 | \subsubsection[UBS: Upstream Biased Scheme (\forcode{ln_dynadv_ubs})]{UBS: Upstream Biased Scheme (\protect\np{ln_dynadv_ubs}{ln\_dynadv\_ubs})} |
---|
496 | \label{subsec:DYN_adv_ubs} |
---|
497 | |
---|
498 | The UBS advection scheme is an upstream biased third order scheme based on |
---|
499 | an upstream-biased parabolic interpolation. |
---|
500 | For example, the evaluation of $u_T^{ubs} $ is done as follows: |
---|
501 | \begin{equation} |
---|
502 | \label{eq:DYN_adv_ubs} |
---|
503 | u_T^{ubs} =\overline u ^i-\;\frac{1}{6} |
---|
504 | \begin{cases} |
---|
505 | u"_{i-1/2}& \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i \geqslant 0$ } \\ |
---|
506 | u"_{i+1/2}& \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i < 0$ } |
---|
507 | \end{cases} |
---|
508 | \end{equation} |
---|
509 | where $u"_{i+1/2} =\delta_{i+1/2} \left[ {\delta_i \left[ u \right]} \right]$. |
---|
510 | This results in a dissipatively dominant (\ie\ hyper-diffusive) truncation error |
---|
511 | \citep{shchepetkin.mcwilliams_OM05}. |
---|
512 | The overall performance of the advection scheme is similar to that reported in \citet{farrow.stevens_JPO95}. |
---|
513 | It is a relatively good compromise between accuracy and smoothness. |
---|
514 | It is not a \emph{positive} scheme, meaning that false extrema are permitted. |
---|
515 | But the amplitudes of the false extrema are significantly reduced over those in the centred second order method. |
---|
516 | As the scheme already includes a diffusion component, it can be used without explicit lateral diffusion on momentum |
---|
517 | (\ie\ \np[=]{ln_dynldf_lap}{ln\_dynldf\_lap}\np[=.false.]{ln_dynldf_bilap}{ln\_dynldf\_bilap}), |
---|
518 | and it is recommended to do so. |
---|
519 | |
---|
520 | The UBS scheme is not used in all directions. |
---|
521 | In the vertical, the centred $2^{nd}$ order evaluation of the advection is preferred, \ie\ $u_{uw}^{ubs}$ and |
---|
522 | $u_{vw}^{ubs}$ in \autoref{eq:DYN_adv_cen2} are used. |
---|
523 | UBS is diffusive and is associated with vertical mixing of momentum. \gmcomment{ gm pursue the |
---|
524 | sentence:Since vertical mixing of momentum is a source term of the TKE equation... } |
---|
525 | |
---|
526 | For stability reasons, the first term in (\autoref{eq:DYN_adv_ubs}), |
---|
527 | which corresponds to a second order centred scheme, is evaluated using the \textit{now} velocity (centred in time), |
---|
528 | while the second term, which is the diffusion part of the scheme, |
---|
529 | is evaluated using the \textit{before} velocity (forward in time). |
---|
530 | This is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the Quick advection scheme. |
---|
531 | |
---|
532 | Note that the UBS and QUICK (Quadratic Upstream Interpolation for Convective Kinematics) schemes only differ by |
---|
533 | one coefficient. |
---|
534 | Replacing $1/6$ by $1/8$ in (\autoref{eq:DYN_adv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}. |
---|
535 | This option is not available through a namelist parameter, since the $1/6$ coefficient is hard coded. |
---|
536 | Nevertheless it is quite easy to make the substitution in the \mdl{dynadv\_ubs} module and obtain a QUICK scheme. |
---|
537 | |
---|
538 | Note also that in the current version of \mdl{dynadv\_ubs}, |
---|
539 | there is also the possibility of using a $4^{th}$ order evaluation of the advective velocity as in ROMS. |
---|
540 | This is an error and should be suppressed soon. |
---|
541 | %%% |
---|
542 | \gmcomment{action : this have to be done} |
---|
543 | %%% |
---|
544 | |
---|
545 | \section[Hydrostatic pressure gradient (\textit{dynhpg.F90})]{Hydrostatic pressure gradient (\protect\mdl{dynhpg})} |
---|
546 | \label{sec:DYN_hpg} |
---|
547 | %------------------------------------------nam_dynhpg--------------------------------------------------- |
---|
548 | |
---|
549 | \begin{listing} |
---|
550 | \nlst{namdyn_hpg} |
---|
551 | \caption{\forcode{&namdyn_hpg}} |
---|
552 | \label{lst:namdyn_hpg} |
---|
553 | \end{listing} |
---|
554 | %------------------------------------------------------------------------------------------------------------- |
---|
555 | |
---|
556 | Options are defined through the \nam{dyn_hpg}{dyn\_hpg} namelist variables. |
---|
557 | The key distinction between the different algorithms used for |
---|
558 | the hydrostatic pressure gradient is the vertical coordinate used, |
---|
559 | since HPG is a \emph{horizontal} pressure gradient, \ie\ computed along geopotential surfaces. |
---|
560 | As a result, any tilt of the surface of the computational levels will require a specific treatment to |
---|
561 | compute the hydrostatic pressure gradient. |
---|
562 | |
---|
563 | The hydrostatic pressure gradient term is evaluated either using a leapfrog scheme, |
---|
564 | \ie\ the density appearing in its expression is centred in time (\emph{now} $\rho$), |
---|
565 | or a semi-implcit scheme. |
---|
566 | At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied. |
---|
567 | |
---|
568 | %-------------------------------------------------------------------------------------------------------------- |
---|
569 | % z-coordinate with full step |
---|
570 | %-------------------------------------------------------------------------------------------------------------- |
---|
571 | \subsection[Full step $Z$-coordinate (\forcode{ln_dynhpg_zco})]{Full step $Z$-coordinate (\protect\np{ln_dynhpg_zco}{ln\_dynhpg\_zco})} |
---|
572 | \label{subsec:DYN_hpg_zco} |
---|
573 | |
---|
574 | The hydrostatic pressure can be obtained by integrating the hydrostatic equation vertically from the surface. |
---|
575 | However, the pressure is large at great depth while its horizontal gradient is several orders of magnitude smaller. |
---|
576 | This may lead to large truncation errors in the pressure gradient terms. |
---|
577 | Thus, the two horizontal components of the hydrostatic pressure gradient are computed directly as follows: |
---|
578 | |
---|
579 | for $k=km$ (surface layer, $jk=1$ in the code) |
---|
580 | \begin{equation} |
---|
581 | \label{eq:DYN_hpg_zco_surf} |
---|
582 | \left\{ |
---|
583 | \begin{aligned} |
---|
584 | \left. \delta_{i+1/2} \left[ p^h \right] \right|_{k=km} |
---|
585 | &= \frac{1}{2} g \ \left. \delta_{i+1/2} \left[ e_{3w} \ \rho \right] \right|_{k=km} \\ |
---|
586 | \left. \delta_{j+1/2} \left[ p^h \right] \right|_{k=km} |
---|
587 | &= \frac{1}{2} g \ \left. \delta_{j+1/2} \left[ e_{3w} \ \rho \right] \right|_{k=km} \\ |
---|
588 | \end{aligned} |
---|
589 | \right. |
---|
590 | \end{equation} |
---|
591 | |
---|
592 | for $1<k<km$ (interior layer) |
---|
593 | \begin{equation} |
---|
594 | \label{eq:DYN_hpg_zco} |
---|
595 | \left\{ |
---|
596 | \begin{aligned} |
---|
597 | \left. \delta_{i+1/2} \left[ p^h \right] \right|_{k} |
---|
598 | &= \left. \delta_{i+1/2} \left[ p^h \right] \right|_{k-1} |
---|
599 | + \frac{1}{2}\;g\; \left. \delta_{i+1/2} \left[ e_{3w} \ \overline {\rho}^{k+1/2} \right] \right|_{k} \\ |
---|
600 | \left. \delta_{j+1/2} \left[ p^h \right] \right|_{k} |
---|
601 | &= \left. \delta_{j+1/2} \left[ p^h \right] \right|_{k-1} |
---|
602 | + \frac{1}{2}\;g\; \left. \delta_{j+1/2} \left[ e_{3w} \ \overline {\rho}^{k+1/2} \right] \right|_{k} \\ |
---|
603 | \end{aligned} |
---|
604 | \right. |
---|
605 | \end{equation} |
---|
606 | |
---|
607 | Note that the $1/2$ factor in (\autoref{eq:DYN_hpg_zco_surf}) is adequate because of the definition of $e_{3w}$ as |
---|
608 | the vertical derivative of the scale factor at the surface level ($z=0$). |
---|
609 | Note also that in case of variable volume level (\texttt{vvl?} defined), |
---|
610 | the surface pressure gradient is included in \autoref{eq:DYN_hpg_zco_surf} and |
---|
611 | \autoref{eq:DYN_hpg_zco} through the space and time variations of the vertical scale factor $e_{3w}$. |
---|
612 | |
---|
613 | %-------------------------------------------------------------------------------------------------------------- |
---|
614 | % z-coordinate with partial step |
---|
615 | %-------------------------------------------------------------------------------------------------------------- |
---|
616 | \subsection[Partial step $Z$-coordinate (\forcode{ln_dynhpg_zps})]{Partial step $Z$-coordinate (\protect\np{ln_dynhpg_zps}{ln\_dynhpg\_zps})} |
---|
617 | \label{subsec:DYN_hpg_zps} |
---|
618 | |
---|
619 | With partial bottom cells, tracers in horizontally adjacent cells generally live at different depths. |
---|
620 | Before taking horizontal gradients between these tracer points, |
---|
621 | a linear interpolation is used to approximate the deeper tracer as if |
---|
622 | it actually lived at the depth of the shallower tracer point. |
---|
623 | |
---|
624 | Apart from this modification, |
---|
625 | the horizontal hydrostatic pressure gradient evaluated in the $z$-coordinate with partial step is exactly as in |
---|
626 | the pure $z$-coordinate case. |
---|
627 | As explained in detail in section \autoref{sec:TRA_zpshde}, |
---|
628 | the nonlinearity of pressure effects in the equation of state is such that |
---|
629 | it is better to interpolate temperature and salinity vertically before computing the density. |
---|
630 | Horizontal gradients of temperature and salinity are needed for the TRA modules, |
---|
631 | which is the reason why the horizontal gradients of density at the deepest model level are computed in |
---|
632 | module \mdl{zpsdhe} located in the TRA directory and described in \autoref{sec:TRA_zpshde}. |
---|
633 | |
---|
634 | %-------------------------------------------------------------------------------------------------------------- |
---|
635 | % s- and s-z-coordinates |
---|
636 | %-------------------------------------------------------------------------------------------------------------- |
---|
637 | \subsection{$S$- and $Z$-$S$-coordinates} |
---|
638 | \label{subsec:DYN_hpg_sco} |
---|
639 | |
---|
640 | Pressure gradient formulations in an $s$-coordinate have been the subject of a vast number of papers |
---|
641 | (\eg, \citet{song_MWR98, shchepetkin.mcwilliams_OM05}). |
---|
642 | A number of different pressure gradient options are coded but the ROMS-like, |
---|
643 | density Jacobian with cubic polynomial method is currently disabled whilst known bugs are under investigation. |
---|
644 | |
---|
645 | $\bullet$ Traditional coding (see for example \citet{madec.delecluse.ea_JPO96}: (\np[=.true.]{ln_dynhpg_sco}{ln\_dynhpg\_sco}) |
---|
646 | \begin{equation} |
---|
647 | \label{eq:DYN_hpg_sco} |
---|
648 | \left\{ |
---|
649 | \begin{aligned} |
---|
650 | - \frac{1} {\rho_o \, e_{1u}} \; \delta_{i+1/2} \left[ p^h \right] |
---|
651 | + \frac{g\; \overline {\rho}^{i+1/2}} {\rho_o \, e_{1u}} \; \delta_{i+1/2} \left[ z_t \right] \\ |
---|
652 | - \frac{1} {\rho_o \, e_{2v}} \; \delta_{j+1/2} \left[ p^h \right] |
---|
653 | + \frac{g\; \overline {\rho}^{j+1/2}} {\rho_o \, e_{2v}} \; \delta_{j+1/2} \left[ z_t \right] \\ |
---|
654 | \end{aligned} |
---|
655 | \right. |
---|
656 | \end{equation} |
---|
657 | |
---|
658 | Where the first term is the pressure gradient along coordinates, |
---|
659 | computed as in \autoref{eq:DYN_hpg_zco_surf} - \autoref{eq:DYN_hpg_zco}, |
---|
660 | and $z_T$ is the depth of the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point |
---|
661 | ($e_{3w}$). |
---|
662 | |
---|
663 | $\bullet$ Traditional coding with adaptation for ice shelf cavities (\np[=.true.]{ln_dynhpg_isf}{ln\_dynhpg\_isf}). |
---|
664 | This scheme need the activation of ice shelf cavities (\np[=.true.]{ln_isfcav}{ln\_isfcav}). |
---|
665 | |
---|
666 | $\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np[=.true.]{ln_dynhpg_prj}{ln\_dynhpg\_prj}) |
---|
667 | |
---|
668 | $\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{shchepetkin.mcwilliams_OM05} |
---|
669 | (\np[=.true.]{ln_dynhpg_djc}{ln\_dynhpg\_djc}) (currently disabled; under development) |
---|
670 | |
---|
671 | Note that expression \autoref{eq:DYN_hpg_sco} is commonly used when the variable volume formulation is activated |
---|
672 | (\texttt{vvl?}) because in that case, even with a flat bottom, |
---|
673 | the coordinate surfaces are not horizontal but follow the free surface \citep{levier.treguier.ea_rpt07}. |
---|
674 | The pressure jacobian scheme (\np[=.true.]{ln_dynhpg_prj}{ln\_dynhpg\_prj}) is available as |
---|
675 | an improved option to \np[=.true.]{ln_dynhpg_sco}{ln\_dynhpg\_sco} when \texttt{vvl?} is active. |
---|
676 | The pressure Jacobian scheme uses a constrained cubic spline to |
---|
677 | reconstruct the density profile across the water column. |
---|
678 | This method maintains the monotonicity between the density nodes. |
---|
679 | The pressure can be calculated by analytical integration of the density profile and |
---|
680 | a pressure Jacobian method is used to solve the horizontal pressure gradient. |
---|
681 | This method can provide a more accurate calculation of the horizontal pressure gradient than the standard scheme. |
---|
682 | |
---|
683 | \subsection{Ice shelf cavity} |
---|
684 | \label{subsec:DYN_hpg_isf} |
---|
685 | |
---|
686 | Beneath an ice shelf, the total pressure gradient is the sum of the pressure gradient due to the ice shelf load and |
---|
687 | the pressure gradient due to the ocean load (\np[=.true.]{ln_dynhpg_isf}{ln\_dynhpg\_isf}).\\ |
---|
688 | |
---|
689 | The main hypothesis to compute the ice shelf load is that the ice shelf is in an isostatic equilibrium. |
---|
690 | The top pressure is computed integrating from surface to the base of the ice shelf a reference density profile |
---|
691 | (prescribed as density of a water at 34.4 PSU and -1.9\deg{C}) and |
---|
692 | corresponds to the water replaced by the ice shelf. |
---|
693 | This top pressure is constant over time. |
---|
694 | A detailed description of this method is described in \citet{losch_JGR08}.\\ |
---|
695 | |
---|
696 | The pressure gradient due to ocean load is computed using the expression \autoref{eq:DYN_hpg_sco} described in |
---|
697 | \autoref{subsec:DYN_hpg_sco}. |
---|
698 | |
---|
699 | %-------------------------------------------------------------------------------------------------------------- |
---|
700 | % Time-scheme |
---|
701 | %-------------------------------------------------------------------------------------------------------------- |
---|
702 | \subsection[Time-scheme (\forcode{ln_dynhpg_imp})]{Time-scheme (\protect\np{ln_dynhpg_imp}{ln\_dynhpg\_imp})} |
---|
703 | \label{subsec:DYN_hpg_imp} |
---|
704 | |
---|
705 | The default time differencing scheme used for the horizontal pressure gradient is a leapfrog scheme and |
---|
706 | therefore the density used in all discrete expressions given above is the \textit{now} density, |
---|
707 | computed from the \textit{now} temperature and salinity. |
---|
708 | In some specific cases |
---|
709 | (usually high resolution simulations over an ocean domain which includes weakly stratified regions) |
---|
710 | the physical phenomenon that controls the time-step is internal gravity waves (IGWs). |
---|
711 | A semi-implicit scheme for doubling the stability limit associated with IGWs can be used |
---|
712 | \citep{brown.campana_MWR78, maltrud.smith.ea_JGR98}. |
---|
713 | It involves the evaluation of the hydrostatic pressure gradient as |
---|
714 | an average over the three time levels $t-\rdt$, $t$, and $t+\rdt$ |
---|
715 | (\ie\ \textit{before}, \textit{now} and \textit{after} time-steps), |
---|
716 | rather than at the central time level $t$ only, as in the standard leapfrog scheme. |
---|
717 | |
---|
718 | $\bullet$ leapfrog scheme (\np[=.true.]{ln_dynhpg_imp}{ln\_dynhpg\_imp}): |
---|
719 | |
---|
720 | \begin{equation} |
---|
721 | \label{eq:DYN_hpg_lf} |
---|
722 | \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \; |
---|
723 | -\frac{1}{\rho_o \,e_{1u} }\delta_{i+1/2} \left[ {p_h^t } \right] |
---|
724 | \end{equation} |
---|
725 | |
---|
726 | $\bullet$ semi-implicit scheme (\np[=.true.]{ln_dynhpg_imp}{ln\_dynhpg\_imp}): |
---|
727 | \begin{equation} |
---|
728 | \label{eq:DYN_hpg_imp} |
---|
729 | \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \; |
---|
730 | -\frac{1}{4\,\rho_o \,e_{1u} } \delta_{i+1/2} \left[ p_h^{t+\rdt} +2\,p_h^t +p_h^{t-\rdt} \right] |
---|
731 | \end{equation} |
---|
732 | |
---|
733 | The semi-implicit time scheme \autoref{eq:DYN_hpg_imp} is made possible without |
---|
734 | significant additional computation since the density can be updated to time level $t+\rdt$ before |
---|
735 | computing the horizontal hydrostatic pressure gradient. |
---|
736 | It can be easily shown that the stability limit associated with the hydrostatic pressure gradient doubles using |
---|
737 | \autoref{eq:DYN_hpg_imp} compared to that using the standard leapfrog scheme \autoref{eq:DYN_hpg_lf}. |
---|
738 | Note that \autoref{eq:DYN_hpg_imp} is equivalent to applying a time filter to the pressure gradient to |
---|
739 | eliminate high frequency IGWs. |
---|
740 | Obviously, when using \autoref{eq:DYN_hpg_imp}, |
---|
741 | the doubling of the time-step is achievable only if no other factors control the time-step, |
---|
742 | such as the stability limits associated with advection or diffusion. |
---|
743 | |
---|
744 | In practice, the semi-implicit scheme is used when \np[=.true.]{ln_dynhpg_imp}{ln\_dynhpg\_imp}. |
---|
745 | In this case, we choose to apply the time filter to temperature and salinity used in the equation of state, |
---|
746 | instead of applying it to the hydrostatic pressure or to the density, |
---|
747 | so that no additional storage array has to be defined. |
---|
748 | The density used to compute the hydrostatic pressure gradient (whatever the formulation) is evaluated as follows: |
---|
749 | \[ |
---|
750 | % \label{eq:DYN_rho_flt} |
---|
751 | \rho^t = \rho( \widetilde{T},\widetilde {S},z_t) |
---|
752 | \quad \text{with} \quad |
---|
753 | \widetilde{X} = 1 / 4 \left( X^{t+\rdt} +2 \,X^t + X^{t-\rdt} \right) |
---|
754 | \] |
---|
755 | |
---|
756 | Note that in the semi-implicit case, it is necessary to save the filtered density, |
---|
757 | an extra three-dimensional field, in the restart file to restart the model with exact reproducibility. |
---|
758 | This option is controlled by \np{nn_dynhpg_rst}{nn\_dynhpg\_rst}, a namelist parameter. |
---|
759 | |
---|
760 | \section[Surface pressure gradient (\textit{dynspg.F90})]{Surface pressure gradient (\protect\mdl{dynspg})} |
---|
761 | \label{sec:DYN_spg} |
---|
762 | %-----------------------------------------nam_dynspg---------------------------------------------------- |
---|
763 | |
---|
764 | \begin{listing} |
---|
765 | \nlst{namdyn_spg} |
---|
766 | \caption{\forcode{&namdyn_spg}} |
---|
767 | \label{lst:namdyn_spg} |
---|
768 | \end{listing} |
---|
769 | %------------------------------------------------------------------------------------------------------------ |
---|
770 | |
---|
771 | Options are defined through the \nam{dyn_spg}{dyn\_spg} namelist variables. |
---|
772 | The surface pressure gradient term is related to the representation of the free surface (\autoref{sec:MB_hor_pg}). |
---|
773 | The main distinction is between the fixed volume case (linear free surface) and |
---|
774 | the variable volume case (nonlinear free surface, \texttt{vvl?} is defined). |
---|
775 | In the linear free surface case (\autoref{subsec:MB_free_surface}) |
---|
776 | the vertical scale factors $e_{3}$ are fixed in time, |
---|
777 | while they are time-dependent in the nonlinear case (\autoref{subsec:MB_free_surface}). |
---|
778 | With both linear and nonlinear free surface, external gravity waves are allowed in the equations, |
---|
779 | which imposes a very small time step when an explicit time stepping is used. |
---|
780 | Two methods are proposed to allow a longer time step for the three-dimensional equations: |
---|
781 | the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:MB_flt?}), |
---|
782 | and the split-explicit free surface described below. |
---|
783 | The extra term introduced in the filtered method is calculated implicitly, |
---|
784 | so that the update of the next velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}. |
---|
785 | |
---|
786 | The form of the surface pressure gradient term depends on how the user wants to |
---|
787 | handle the fast external gravity waves that are a solution of the analytical equation (\autoref{sec:MB_hor_pg}). |
---|
788 | Three formulations are available, all controlled by a CPP key (ln\_dynspg\_xxx): |
---|
789 | an explicit formulation which requires a small time step; |
---|
790 | a filtered free surface formulation which allows a larger time step by |
---|
791 | adding a filtering term into the momentum equation; |
---|
792 | and a split-explicit free surface formulation, described below, which also allows a larger time step. |
---|
793 | |
---|
794 | The extra term introduced in the filtered method is calculated implicitly, so that a solver is used to compute it. |
---|
795 | As a consequence the update of the $next$ velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}. |
---|
796 | |
---|
797 | %-------------------------------------------------------------------------------------------------------------- |
---|
798 | % Explicit free surface formulation |
---|
799 | %-------------------------------------------------------------------------------------------------------------- |
---|
800 | \subsection[Explicit free surface (\forcode{ln_dynspg_exp})]{Explicit free surface (\protect\np{ln_dynspg_exp}{ln\_dynspg\_exp})} |
---|
801 | \label{subsec:DYN_spg_exp} |
---|
802 | |
---|
803 | In the explicit free surface formulation (\np{ln_dynspg_exp}{ln\_dynspg\_exp} set to true), |
---|
804 | the model time step is chosen to be small enough to resolve the external gravity waves |
---|
805 | (typically a few tens of seconds). |
---|
806 | The surface pressure gradient, evaluated using a leap-frog scheme (\ie\ centered in time), |
---|
807 | is thus simply given by : |
---|
808 | \begin{equation} |
---|
809 | \label{eq:DYN_spg_exp} |
---|
810 | \left\{ |
---|
811 | \begin{aligned} |
---|
812 | - \frac{1}{e_{1u}\,\rho_o} \; \delta_{i+1/2} \left[ \,\rho \,\eta\, \right] \\ |
---|
813 | - \frac{1}{e_{2v}\,\rho_o} \; \delta_{j+1/2} \left[ \,\rho \,\eta\, \right] |
---|
814 | \end{aligned} |
---|
815 | \right. |
---|
816 | \end{equation} |
---|
817 | |
---|
818 | Note that in the non-linear free surface case (\ie\ \texttt{vvl?} defined), |
---|
819 | the surface pressure gradient is already included in the momentum tendency through |
---|
820 | the level thickness variation allowed in the computation of the hydrostatic pressure gradient. |
---|
821 | Thus, nothing is done in the \mdl{dynspg\_exp} module. |
---|
822 | |
---|
823 | %-------------------------------------------------------------------------------------------------------------- |
---|
824 | % Split-explict free surface formulation |
---|
825 | %-------------------------------------------------------------------------------------------------------------- |
---|
826 | \subsection[Split-explicit free surface (\forcode{ln_dynspg_ts})]{Split-explicit free surface (\protect\np{ln_dynspg_ts}{ln\_dynspg\_ts})} |
---|
827 | \label{subsec:DYN_spg_ts} |
---|
828 | %------------------------------------------namsplit----------------------------------------------------------- |
---|
829 | % |
---|
830 | %\nlst{namsplit} |
---|
831 | %------------------------------------------------------------------------------------------------------------- |
---|
832 | |
---|
833 | The split-explicit free surface formulation used in \NEMO\ (\np{ln_dynspg_ts}{ln\_dynspg\_ts} set to true), |
---|
834 | also called the time-splitting formulation, follows the one proposed by \citet{shchepetkin.mcwilliams_OM05}. |
---|
835 | The general idea is to solve the free surface equation and the associated barotropic velocity equations with |
---|
836 | a smaller time step than $\rdt$, the time step used for the three dimensional prognostic variables |
---|
837 | (\autoref{fig:DYN_spg_ts}). |
---|
838 | The size of the small time step, $\rdt_e$ (the external mode or barotropic time step) is provided through |
---|
839 | the \np{nn_baro}{nn\_baro} namelist parameter as: $\rdt_e = \rdt / nn\_baro$. |
---|
840 | This parameter can be optionally defined automatically (\np[=.true.]{ln_bt_nn_auto}{ln\_bt\_nn\_auto}) considering that |
---|
841 | the stability of the barotropic system is essentially controled by external waves propagation. |
---|
842 | Maximum Courant number is in that case time independent, and easily computed online from the input bathymetry. |
---|
843 | Therefore, $\rdt_e$ is adjusted so that the Maximum allowed Courant number is smaller than \np{rn_bt_cmax}{rn\_bt\_cmax}. |
---|
844 | |
---|
845 | %%% |
---|
846 | The barotropic mode solves the following equations: |
---|
847 | % \begin{subequations} |
---|
848 | % \label{eq:DYN_BT} |
---|
849 | \begin{equation} |
---|
850 | \label{eq:DYN_BT_dyn} |
---|
851 | \frac{\partial {\mathrm \overline{{\mathbf U}}_h} }{\partial t}= |
---|
852 | -f\;{\mathrm {\mathbf k}}\times {\mathrm \overline{{\mathbf U}}_h} |
---|
853 | -g\nabla _h \eta -\frac{c_b^{\textbf U}}{H+\eta} \mathrm {\overline{{\mathbf U}}_h} + \mathrm {\overline{\mathbf G}} |
---|
854 | \end{equation} |
---|
855 | \[ |
---|
856 | % \label{eq:DYN_BT_ssh} |
---|
857 | \frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\mathrm{\mathbf \overline{U}}}_h \,} \right]+P-E |
---|
858 | \] |
---|
859 | % \end{subequations} |
---|
860 | where $\mathrm {\overline{\mathbf G}}$ is a forcing term held constant, containing coupling term between modes, |
---|
861 | surface atmospheric forcing as well as slowly varying barotropic terms not explicitly computed to gain efficiency. |
---|
862 | The third term on the right hand side of \autoref{eq:DYN_BT_dyn} represents the bottom stress |
---|
863 | (see section \autoref{sec:ZDF_drg}), explicitly accounted for at each barotropic iteration. |
---|
864 | Temporal discretization of the system above follows a three-time step Generalized Forward Backward algorithm |
---|
865 | detailed in \citet{shchepetkin.mcwilliams_OM05}. |
---|
866 | AB3-AM4 coefficients used in \NEMO\ follow the second-order accurate, |
---|
867 | "multi-purpose" stability compromise as defined in \citet{shchepetkin.mcwilliams_ibk09} |
---|
868 | (see their figure 12, lower left). |
---|
869 | |
---|
870 | %> > > > > > > > > > > > > > > > > > > > > > > > > > > > |
---|
871 | \begin{figure}[!t] |
---|
872 | \centering |
---|
873 | \includegraphics[width=0.66\textwidth]{Fig_DYN_dynspg_ts} |
---|
874 | \caption[Split-explicit time stepping scheme for the external and internal modes]{ |
---|
875 | Schematic of the split-explicit time stepping scheme for the external and internal modes. |
---|
876 | Time increases to the right. |
---|
877 | In this particular exemple, |
---|
878 | a boxcar averaging window over \np{nn_baro}{nn\_baro} barotropic time steps is used |
---|
879 | (\np[=1]{nn_bt_flt}{nn\_bt\_flt}) and \np[=5]{nn_baro}{nn\_baro}. |
---|
880 | Internal mode time steps (which are also the model time steps) are denoted by |
---|
881 | $t-\rdt$, $t$ and $t+\rdt$. |
---|
882 | Variables with $k$ superscript refer to instantaneous barotropic variables, |
---|
883 | $< >$ and $<< >>$ operator refer to time filtered variables using respectively primary |
---|
884 | (red vertical bars) and secondary weights (blue vertical bars). |
---|
885 | The former are used to obtain time filtered quantities at $t+\rdt$ while |
---|
886 | the latter are used to obtain time averaged transports to advect tracers. |
---|
887 | a) Forward time integration: |
---|
888 | \protect\np[=.true.]{ln_bt_fw}{ln\_bt\_fw}, \protect\np[=.true.]{ln_bt_av}{ln\_bt\_av}. |
---|
889 | b) Centred time integration: |
---|
890 | \protect\np[=.false.]{ln_bt_fw}{ln\_bt\_fw}, \protect\np[=.true.]{ln_bt_av}{ln\_bt\_av}. |
---|
891 | c) Forward time integration with no time filtering (POM-like scheme): |
---|
892 | \protect\np[=.true.]{ln_bt_fw}{ln\_bt\_fw}, \protect\np[=.false.]{ln_bt_av}{ln\_bt\_av}.} |
---|
893 | \label{fig:DYN_spg_ts} |
---|
894 | \end{figure} |
---|
895 | %> > > > > > > > > > > > > > > > > > > > > > > > > > > > |
---|
896 | |
---|
897 | In the default case (\np[=.true.]{ln_bt_fw}{ln\_bt\_fw}), |
---|
898 | the external mode is integrated between \textit{now} and \textit{after} baroclinic time-steps |
---|
899 | (\autoref{fig:DYN_spg_ts}a). |
---|
900 | To avoid aliasing of fast barotropic motions into three dimensional equations, |
---|
901 | time filtering is eventually applied on barotropic quantities (\np[=.true.]{ln_bt_av}{ln\_bt\_av}). |
---|
902 | In that case, the integration is extended slightly beyond \textit{after} time step to |
---|
903 | provide time filtered quantities. |
---|
904 | These are used for the subsequent initialization of the barotropic mode in the following baroclinic step. |
---|
905 | Since external mode equations written at baroclinic time steps finally follow a forward time stepping scheme, |
---|
906 | asselin filtering is not applied to barotropic quantities.\\ |
---|
907 | Alternatively, one can choose to integrate barotropic equations starting from \textit{before} time step |
---|
908 | (\np[=.false.]{ln_bt_fw}{ln\_bt\_fw}). |
---|
909 | Although more computationaly expensive ( \np{nn_baro}{nn\_baro} additional iterations are indeed necessary), |
---|
910 | the baroclinic to barotropic forcing term given at \textit{now} time step become centred in |
---|
911 | the middle of the integration window. |
---|
912 | It can easily be shown that this property removes part of splitting errors between modes, |
---|
913 | which increases the overall numerical robustness. |
---|
914 | %references to Patrick Marsaleix' work here. Also work done by SHOM group. |
---|
915 | |
---|
916 | %%% |
---|
917 | |
---|
918 | As far as tracer conservation is concerned, |
---|
919 | barotropic velocities used to advect tracers must also be updated at \textit{now} time step. |
---|
920 | This implies to change the traditional order of computations in \NEMO: |
---|
921 | most of momentum trends (including the barotropic mode calculation) updated first, tracers' after. |
---|
922 | This \textit{de facto} makes semi-implicit hydrostatic pressure gradient |
---|
923 | (see section \autoref{subsec:DYN_hpg_imp}) |
---|
924 | and time splitting not compatible. |
---|
925 | Advective barotropic velocities are obtained by using a secondary set of filtering weights, |
---|
926 | uniquely defined from the filter coefficients used for the time averaging (\citet{shchepetkin.mcwilliams_OM05}). |
---|
927 | Consistency between the time averaged continuity equation and the time stepping of tracers is here the key to |
---|
928 | obtain exact conservation. |
---|
929 | |
---|
930 | %%% |
---|
931 | |
---|
932 | One can eventually choose to feedback instantaneous values by not using any time filter |
---|
933 | (\np[=.false.]{ln_bt_av}{ln\_bt\_av}). |
---|
934 | In that case, external mode equations are continuous in time, |
---|
935 | \ie\ they are not re-initialized when starting a new sub-stepping sequence. |
---|
936 | This is the method used so far in the POM model, the stability being maintained by |
---|
937 | refreshing at (almost) each barotropic time step advection and horizontal diffusion terms. |
---|
938 | Since the latter terms have not been added in \NEMO\ for computational efficiency, |
---|
939 | removing time filtering is not recommended except for debugging purposes. |
---|
940 | This may be used for instance to appreciate the damping effect of the standard formulation on |
---|
941 | external gravity waves in idealized or weakly non-linear cases. |
---|
942 | Although the damping is lower than for the filtered free surface, |
---|
943 | it is still significant as shown by \citet{levier.treguier.ea_rpt07} in the case of an analytical barotropic Kelvin wave. |
---|
944 | |
---|
945 | %>>>>>=============== |
---|
946 | \gmcomment{ %%% copy from griffies Book |
---|
947 | |
---|
948 | \textbf{title: Time stepping the barotropic system } |
---|
949 | |
---|
950 | Assume knowledge of the full velocity and tracer fields at baroclinic time $\tau$. |
---|
951 | Hence, we can update the surface height and vertically integrated velocity with a leap-frog scheme using |
---|
952 | the small barotropic time step $\rdt$. |
---|
953 | We have |
---|
954 | |
---|
955 | \[ |
---|
956 | % \label{eq:DYN_spg_ts_eta} |
---|
957 | \eta^{(b)}(\tau,t_{n+1}) - \eta^{(b)}(\tau,t_{n+1}) (\tau,t_{n-1}) |
---|
958 | = 2 \rdt \left[-\nabla \cdot \textbf{U}^{(b)}(\tau,t_n) + \text{EMP}_w(\tau) \right] |
---|
959 | \] |
---|
960 | \begin{multline*} |
---|
961 | % \label{eq:DYN_spg_ts_u} |
---|
962 | \textbf{U}^{(b)}(\tau,t_{n+1}) - \textbf{U}^{(b)}(\tau,t_{n-1}) \\ |
---|
963 | = 2\rdt \left[ - f \textbf{k} \times \textbf{U}^{(b)}(\tau,t_{n}) |
---|
964 | - H(\tau) \nabla p_s^{(b)}(\tau,t_{n}) +\textbf{M}(\tau) \right] |
---|
965 | \end{multline*} |
---|
966 | \ |
---|
967 | |
---|
968 | In these equations, araised (b) denotes values of surface height and vertically integrated velocity updated with |
---|
969 | the barotropic time steps. |
---|
970 | The $\tau$ time label on $\eta^{(b)}$ and $U^{(b)}$ denotes the baroclinic time at which |
---|
971 | the vertically integrated forcing $\textbf{M}(\tau)$ |
---|
972 | (note that this forcing includes the surface freshwater forcing), |
---|
973 | the tracer fields, the freshwater flux $\text{EMP}_w(\tau)$, |
---|
974 | and total depth of the ocean $H(\tau)$ are held for the duration of the barotropic time stepping over |
---|
975 | a single cycle. |
---|
976 | This is also the time that sets the barotropic time steps via |
---|
977 | \[ |
---|
978 | % \label{eq:DYN_spg_ts_t} |
---|
979 | t_n=\tau+n\rdt |
---|
980 | \] |
---|
981 | with $n$ an integer. |
---|
982 | The density scaled surface pressure is evaluated via |
---|
983 | \[ |
---|
984 | % \label{eq:DYN_spg_ts_ps} |
---|
985 | p_s^{(b)}(\tau,t_{n}) = |
---|
986 | \begin{cases} |
---|
987 | g \;\eta_s^{(b)}(\tau,t_{n}) \;\rho(\tau)_{k=1}) / \rho_o & \text{non-linear case} \\ |
---|
988 | g \;\eta_s^{(b)}(\tau,t_{n}) & \text{linear case} |
---|
989 | \end{cases} |
---|
990 | \] |
---|
991 | To get started, we assume the following initial conditions |
---|
992 | \[ |
---|
993 | % \label{eq:DYN_spg_ts_eta} |
---|
994 | \begin{split} |
---|
995 | \eta^{(b)}(\tau,t_{n=0}) &= \overline{\eta^{(b)}(\tau)} \\ |
---|
996 | \eta^{(b)}(\tau,t_{n=1}) &= \eta^{(b)}(\tau,t_{n=0}) + \rdt \ \text{RHS}_{n=0} |
---|
997 | \end{split} |
---|
998 | \] |
---|
999 | with |
---|
1000 | \[ |
---|
1001 | % \label{eq:DYN_spg_ts_etaF} |
---|
1002 | \overline{\eta^{(b)}(\tau)} = \frac{1}{N+1} \sum\limits_{n=0}^N \eta^{(b)}(\tau-\rdt,t_{n}) |
---|
1003 | \] |
---|
1004 | the time averaged surface height taken from the previous barotropic cycle. |
---|
1005 | Likewise, |
---|
1006 | \[ |
---|
1007 | % \label{eq:DYN_spg_ts_u} |
---|
1008 | \textbf{U}^{(b)}(\tau,t_{n=0}) = \overline{\textbf{U}^{(b)}(\tau)} \\ \\ |
---|
1009 | \textbf{U}(\tau,t_{n=1}) = \textbf{U}^{(b)}(\tau,t_{n=0}) + \rdt \ \text{RHS}_{n=0} |
---|
1010 | \] |
---|
1011 | with |
---|
1012 | \[ |
---|
1013 | % \label{eq:DYN_spg_ts_u} |
---|
1014 | \overline{\textbf{U}^{(b)}(\tau)} = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau-\rdt,t_{n}) |
---|
1015 | \] |
---|
1016 | the time averaged vertically integrated transport. |
---|
1017 | Notably, there is no Robert-Asselin time filter used in the barotropic portion of the integration. |
---|
1018 | |
---|
1019 | Upon reaching $t_{n=N} = \tau + 2\rdt \tau$ , |
---|
1020 | the vertically integrated velocity is time averaged to produce the updated vertically integrated velocity at |
---|
1021 | baroclinic time $\tau + \rdt \tau$ |
---|
1022 | \[ |
---|
1023 | % \label{eq:DYN_spg_ts_u} |
---|
1024 | \textbf{U}(\tau+\rdt) = \overline{\textbf{U}^{(b)}(\tau+\rdt)} = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau,t_{n}) |
---|
1025 | \] |
---|
1026 | The surface height on the new baroclinic time step is then determined via a baroclinic leap-frog using |
---|
1027 | the following form |
---|
1028 | |
---|
1029 | \begin{equation} |
---|
1030 | \label{eq:DYN_spg_ts_ssh} |
---|
1031 | \eta(\tau+\Delta) - \eta^{F}(\tau-\Delta) = 2\rdt \ \left[ - \nabla \cdot \textbf{U}(\tau) + \text{EMP}_w \right] |
---|
1032 | \end{equation} |
---|
1033 | |
---|
1034 | The use of this "big-leap-frog" scheme for the surface height ensures compatibility between |
---|
1035 | the mass/volume budgets and the tracer budgets. |
---|
1036 | More discussion of this point is provided in Chapter 10 (see in particular Section 10.2). |
---|
1037 | |
---|
1038 | In general, some form of time filter is needed to maintain integrity of the surface height field due to |
---|
1039 | the leap-frog splitting mode in equation \autoref{eq:DYN_spg_ts_ssh}. |
---|
1040 | We have tried various forms of such filtering, |
---|
1041 | with the following method discussed in \cite{griffies.pacanowski.ea_MWR01} chosen due to |
---|
1042 | its stability and reasonably good maintenance of tracer conservation properties (see ??). |
---|
1043 | |
---|
1044 | \begin{equation} |
---|
1045 | \label{eq:DYN_spg_ts_sshf} |
---|
1046 | \eta^{F}(\tau-\Delta) = \overline{\eta^{(b)}(\tau)} |
---|
1047 | \end{equation} |
---|
1048 | Another approach tried was |
---|
1049 | |
---|
1050 | \[ |
---|
1051 | % \label{eq:DYN_spg_ts_sshf2} |
---|
1052 | \eta^{F}(\tau-\Delta) = \eta(\tau) |
---|
1053 | + (\alpha/2) \left[\overline{\eta^{(b)}}(\tau+\rdt) |
---|
1054 | + \overline{\eta^{(b)}}(\tau-\rdt) -2 \;\eta(\tau) \right] |
---|
1055 | \] |
---|
1056 | |
---|
1057 | which is useful since it isolates all the time filtering aspects into the term multiplied by $\alpha$. |
---|
1058 | This isolation allows for an easy check that tracer conservation is exact when |
---|
1059 | eliminating tracer and surface height time filtering (see ?? for more complete discussion). |
---|
1060 | However, in the general case with a non-zero $\alpha$, |
---|
1061 | the filter \autoref{eq:DYN_spg_ts_sshf} was found to be more conservative, and so is recommended. |
---|
1062 | |
---|
1063 | } %%end gm comment (copy of griffies book) |
---|
1064 | |
---|
1065 | %>>>>>=============== |
---|
1066 | |
---|
1067 | %-------------------------------------------------------------------------------------------------------------- |
---|
1068 | % Filtered free surface formulation |
---|
1069 | %-------------------------------------------------------------------------------------------------------------- |
---|
1070 | \subsection{Filtered free surface (\forcode{dynspg_flt?})} |
---|
1071 | \label{subsec:DYN_spg_fltp} |
---|
1072 | |
---|
1073 | The filtered formulation follows the \citet{roullet.madec_JGR00} implementation. |
---|
1074 | The extra term introduced in the equations (see \autoref{subsec:MB_free_surface}) is solved implicitly. |
---|
1075 | The elliptic solvers available in the code are documented in \autoref{chap:MISC}. |
---|
1076 | |
---|
1077 | %% gm %%======>>>> given here the discrete eqs provided to the solver |
---|
1078 | \gmcomment{ %%% copy from chap-model basics |
---|
1079 | \[ |
---|
1080 | % \label{eq:DYN_spg_flt} |
---|
1081 | \frac{\partial {\mathrm {\mathbf U}}_h }{\partial t}= {\mathrm {\mathbf M}} |
---|
1082 | - g \nabla \left( \tilde{\rho} \ \eta \right) |
---|
1083 | - g \ T_c \nabla \left( \widetilde{\rho} \ \partial_t \eta \right) |
---|
1084 | \] |
---|
1085 | where $T_c$, is a parameter with dimensions of time which characterizes the force, |
---|
1086 | $\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, |
---|
1087 | and $\mathrm {\mathbf M}$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient, |
---|
1088 | non-linear and viscous terms in \autoref{eq:MB_dyn}. |
---|
1089 | } %end gmcomment |
---|
1090 | |
---|
1091 | Note that in the linear free surface formulation (\texttt{vvl?} not defined), |
---|
1092 | the ocean depth is time-independent and so is the matrix to be inverted. |
---|
1093 | It is computed once and for all and applies to all ocean time steps. |
---|
1094 | |
---|
1095 | \section[Lateral diffusion term and operators (\textit{dynldf.F90})]{Lateral diffusion term and operators (\protect\mdl{dynldf})} |
---|
1096 | \label{sec:DYN_ldf} |
---|
1097 | %------------------------------------------nam_dynldf---------------------------------------------------- |
---|
1098 | |
---|
1099 | \begin{listing} |
---|
1100 | \nlst{namdyn_ldf} |
---|
1101 | \caption{\forcode{&namdyn_ldf}} |
---|
1102 | \label{lst:namdyn_ldf} |
---|
1103 | \end{listing} |
---|
1104 | %------------------------------------------------------------------------------------------------------------- |
---|
1105 | |
---|
1106 | Options are defined through the \nam{dyn_ldf}{dyn\_ldf} namelist variables. |
---|
1107 | The options available for lateral diffusion are to use either laplacian (rotated or not) or biharmonic operators. |
---|
1108 | The coefficients may be constant or spatially variable; |
---|
1109 | the description of the coefficients is found in the chapter on lateral physics (\autoref{chap:LDF}). |
---|
1110 | The lateral diffusion of momentum is evaluated using a forward scheme, |
---|
1111 | \ie\ the velocity appearing in its expression is the \textit{before} velocity in time, |
---|
1112 | except for the pure vertical component that appears when a tensor of rotation is used. |
---|
1113 | This latter term is solved implicitly together with the vertical diffusion term (see \autoref{chap:TD}). |
---|
1114 | |
---|
1115 | At the lateral boundaries either free slip, |
---|
1116 | no slip or partial slip boundary conditions are applied according to the user's choice (see \autoref{chap:LBC}). |
---|
1117 | |
---|
1118 | \gmcomment{ |
---|
1119 | Hyperviscous operators are frequently used in the simulation of turbulent flows to |
---|
1120 | control the dissipation of unresolved small scale features. |
---|
1121 | Their primary role is to provide strong dissipation at the smallest scale supported by |
---|
1122 | the grid while minimizing the impact on the larger scale features. |
---|
1123 | Hyperviscous operators are thus designed to be more scale selective than the traditional, |
---|
1124 | physically motivated Laplace operator. |
---|
1125 | In finite difference methods, |
---|
1126 | the biharmonic operator is frequently the method of choice to achieve this scale selective dissipation since |
---|
1127 | its damping time (\ie\ its spin down time) scale like $\lambda^{-4}$ for disturbances of wavelength $\lambda$ |
---|
1128 | (so that short waves damped more rapidelly than long ones), |
---|
1129 | whereas the Laplace operator damping time scales only like $\lambda^{-2}$. |
---|
1130 | } |
---|
1131 | |
---|
1132 | % Vertical diffusion term |
---|
1133 | % External Forcing |
---|
1134 | % Wetting and drying |
---|
1135 | % Time evolution term |
---|