1 | % ================================================================ |
---|
2 | % Chapter Ñ Vertical Ocean Physics (ZDF) |
---|
3 | % ================================================================ |
---|
4 | \chapter{Vertical Ocean Physics (ZDF)} |
---|
5 | \label{ZDF} |
---|
6 | \minitoc |
---|
7 | |
---|
8 | %gm% Add here a small introduction to ZDF and naming of the different physics (similar to what have been written for TRA and DYN. |
---|
9 | |
---|
10 | |
---|
11 | \newpage |
---|
12 | $\ $\newline % force a new ligne |
---|
13 | |
---|
14 | |
---|
15 | % ================================================================ |
---|
16 | % Vertical Mixing |
---|
17 | % ================================================================ |
---|
18 | \section{Vertical Mixing} |
---|
19 | \label{ZDF_zdf} |
---|
20 | |
---|
21 | The discrete form of the ocean subgrid scale physics has been presented in |
---|
22 | \S\ref{TRA_zdf} and \S\ref{DYN_zdf}. At the surface and bottom boundaries, |
---|
23 | the turbulent fluxes of momentum, heat and salt have to be defined. At the |
---|
24 | surface they are prescribed from the surface forcing (see Chap.~\ref{SBC}), |
---|
25 | while at the bottom they are set to zero for heat and salt, unless a geothermal |
---|
26 | flux forcing is prescribed as a bottom boundary condition ($i.e.$ \key{trabbl} |
---|
27 | defined, see \S\ref{TRA_bbc}), and specified through a bottom friction |
---|
28 | parameterisation for momentum (see \S\ref{ZDF_bfr}). |
---|
29 | |
---|
30 | In this section we briefly discuss the various choices offered to compute |
---|
31 | the vertical eddy viscosity and diffusivity coefficients, $A_u^{vm}$ , |
---|
32 | $A_v^{vm}$ and $A^{vT}$ ($A^{vS}$), defined at $uw$-, $vw$- and $w$- |
---|
33 | points, respectively (see \S\ref{TRA_zdf} and \S\ref{DYN_zdf}). These |
---|
34 | coefficients can be assumed to be either constant, or a function of the local |
---|
35 | Richardson number, or computed from a turbulent closure model (either |
---|
36 | TKE or KPP formulation). The computation of these coefficients is initialized |
---|
37 | in the \mdl{zdfini} module and performed in the \mdl{zdfric}, \mdl{zdftke} or |
---|
38 | \mdl{zdfkpp} modules. The trends due to the vertical momentum and tracer |
---|
39 | diffusion, including the surface forcing, are computed and added to the |
---|
40 | general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively. |
---|
41 | These trends can be computed using either a forward time stepping scheme |
---|
42 | (namelist parameter \np{np\_zdfexp}=true) or a backward time stepping |
---|
43 | scheme (\np{np\_zdfexp}=false) depending on the magnitude of the mixing |
---|
44 | coefficients, and thus of the formulation used (see \S\ref{DOM_nxt}). |
---|
45 | |
---|
46 | % ------------------------------------------------------------------------------------------------------------- |
---|
47 | % Constant |
---|
48 | % ------------------------------------------------------------------------------------------------------------- |
---|
49 | \subsection{Constant (\key{zdfcst})} |
---|
50 | \label{ZDF_cst} |
---|
51 | %--------------------------------------------namzdf--------------------------------------------------------- |
---|
52 | \namdisplay{namzdf} |
---|
53 | %-------------------------------------------------------------------------------------------------------------- |
---|
54 | |
---|
55 | When \key{zdfcst} is defined, the momentum and tracer vertical eddy coefficients |
---|
56 | are set to constant values over the whole ocean. This is the crudest way to define |
---|
57 | the vertical ocean physics. It is recommended that this option is only used in |
---|
58 | process studies, not in basin scale simulations. Typical values used in this case are: |
---|
59 | \begin{align*} |
---|
60 | A_u^{vm} = A_v^{vm} &= 1.2\ 10^{-4}~m^2.s^{-1} \\ |
---|
61 | \\ |
---|
62 | A^{vT} = A^{vS} &= 1.2\ 10^{-5}~m^2.s^{-1} |
---|
63 | \end{align*} |
---|
64 | |
---|
65 | These values are set through the \np{avm0} and \np{avt0} namelist parameters. |
---|
66 | In all cases, do not use values smaller that those associated with the molecular |
---|
67 | viscosity and diffusivity, that is $\sim10^{-6}~m^2.s^{-1}$ for momentum, |
---|
68 | $\sim10^{-7}~m^2.s^{-1}$ for temperature and $\sim10^{-9}~m^2.s^{-1}$ for salinity. |
---|
69 | |
---|
70 | |
---|
71 | % ------------------------------------------------------------------------------------------------------------- |
---|
72 | % Richardson Number Dependent |
---|
73 | % ------------------------------------------------------------------------------------------------------------- |
---|
74 | \subsection{Richardson Number Dependent (\key{zdfric})} |
---|
75 | \label{ZDF_ric} |
---|
76 | |
---|
77 | %--------------------------------------------namric--------------------------------------------------------- |
---|
78 | \namdisplay{namzdf_ric} |
---|
79 | %-------------------------------------------------------------------------------------------------------------- |
---|
80 | |
---|
81 | When \key{zdfric} is defined, a local Richardson number dependent formulation |
---|
82 | for the vertical momentum and tracer eddy coefficients is set. The vertical mixing |
---|
83 | coefficients are diagnosed from the large scale variables computed by the model. |
---|
84 | \textit{In situ} measurements have been used to link vertical turbulent activity to |
---|
85 | large scale ocean structures. The hypothesis of a mixing mainly maintained by the |
---|
86 | growth of Kelvin-Helmholtz like instabilities leads to a dependency between the |
---|
87 | vertical eddy coefficients and the local Richardson number ($i.e.$ the |
---|
88 | ratio of stratification to vertical shear). Following \citet{Pacanowski_Philander_JPO81}, the following |
---|
89 | formulation has been implemented: |
---|
90 | \begin{equation} \label{Eq_zdfric} |
---|
91 | \left\{ \begin{aligned} |
---|
92 | A^{vT} &= \frac {A_{ric}^{vT}}{\left( 1+a \; Ri \right)^n} + A_b^{vT} \\ |
---|
93 | \\ |
---|
94 | A^{vm} &= \frac{A^{vT} }{\left( 1+ a \;Ri \right) } + A_b^{vm} |
---|
95 | \end{aligned} \right. |
---|
96 | \end{equation} |
---|
97 | where $Ri = N^2 / \left(\partial_z \textbf{U}_h \right)^2$ is the local Richardson |
---|
98 | number, $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}), |
---|
99 | $A_b^{vT} $ and $A_b^{vm}$ are the constant background values set as in the |
---|
100 | constant case (see \S\ref{ZDF_cst}), and $A_{ric}^{vT} = 10^{-4}~m^2.s^{-1}$ |
---|
101 | is the maximum value that can be reached by the coefficient when $Ri\leq 0$, |
---|
102 | $a=5$ and $n=2$. The last three values can be modified by setting the |
---|
103 | \np{avmri}, \np{alp} and \np{nric} namelist parameters, respectively. |
---|
104 | |
---|
105 | % ------------------------------------------------------------------------------------------------------------- |
---|
106 | % TKE Turbulent Closure Scheme |
---|
107 | % ------------------------------------------------------------------------------------------------------------- |
---|
108 | \subsection{TKE Turbulent Closure Scheme (\key{zdftke})} |
---|
109 | \label{ZDF_tke} |
---|
110 | |
---|
111 | %--------------------------------------------namzdf_tke-------------------------------------------------- |
---|
112 | \namdisplay{namzdf_tke} |
---|
113 | %-------------------------------------------------------------------------------------------------------------- |
---|
114 | |
---|
115 | The vertical eddy viscosity and diffusivity coefficients are computed from a TKE |
---|
116 | turbulent closure model based on a prognostic equation for $\bar{e}$, the turbulent |
---|
117 | kinetic energy, and a closure assumption for the turbulent length scales. This |
---|
118 | turbulent closure model has been developed by \citet{Bougeault1989} in the |
---|
119 | atmospheric case, adapted by \citet{Gaspar1990} for the oceanic case, and |
---|
120 | embedded in OPA, the ancestor of NEMO, by \citet{Blanke1993} for equatorial Atlantic |
---|
121 | simulations. Since then, significant modifications have been introduced by |
---|
122 | \citet{Madec1998} in both the implementation and the formulation of the mixing |
---|
123 | length scale. The time evolution of $\bar{e}$ is the result of the production of |
---|
124 | $\bar{e}$ through vertical shear, its destruction through stratification, its vertical |
---|
125 | diffusion, and its dissipation of \citet{Kolmogorov1942} type: |
---|
126 | \begin{equation} \label{Eq_zdftke_e} |
---|
127 | \frac{\partial \bar{e}}{\partial t} = |
---|
128 | \frac{A^{vm}}{{e_3}^2 }\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2 |
---|
129 | +\left( {\frac{\partial v}{\partial k}} \right)^2} \right] |
---|
130 | -A^{vT}\,N^2 |
---|
131 | +\frac{1}{e_3} \;\frac{\partial }{\partial k}\left[ {\frac{A^{vm}}{e_3 } |
---|
132 | \;\frac{\partial \bar{e}}{\partial k}} \right] |
---|
133 | - c_\epsilon \;\frac{\bar {e}^{3/2}}{l_\epsilon } |
---|
134 | \end{equation} |
---|
135 | \begin{equation} \label{Eq_zdftke_kz} |
---|
136 | \begin{split} |
---|
137 | A^{vm} &= C_k\ l_k\ \sqrt {\bar{e}\; } \\ |
---|
138 | A^{vT} &= A^{vm} / P_{rt} |
---|
139 | \end{split} |
---|
140 | \end{equation} |
---|
141 | where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}), |
---|
142 | $l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales, |
---|
143 | $P_{rt}$ is the Prandtl number. The constants $C_k = 0.1$ and |
---|
144 | $C_\epsilon = \sqrt {2} /2$ $\approx 0.7$ are designed to deal with vertical mixing |
---|
145 | at any depth \citep{Gaspar1990}. They are set through namelist parameters |
---|
146 | \np{nn\_ediff} and \np{nn\_ediss}. $P_{rt}$ can be set to unity or, following |
---|
147 | \citet{Blanke1993}, be a function of the local Richardson number, $R_i$: |
---|
148 | \begin{align*} \label{Eq_prt} |
---|
149 | P_{rt} = \begin{cases} |
---|
150 | \ \ \ 1 & \text{if $\ R_i \leq 0.2$} \\ |
---|
151 | 5\,R_i & \text{if $\ 0.2 \leq R_i \leq 2$} \\ |
---|
152 | \ \ 10 & \text{if $\ 2 \leq R_i$} |
---|
153 | \end{cases} |
---|
154 | \end{align*} |
---|
155 | The choice of $P_{rt}$ is controlled by the \np{nn\_pdl} namelist parameter. |
---|
156 | |
---|
157 | For computational efficiency, the original formulation of the turbulent length |
---|
158 | scales proposed by \citet{Gaspar1990} has been simplified. Four formulations |
---|
159 | are proposed, the choice of which is controlled by the \np{nn\_mxl} namelist |
---|
160 | parameter. The first two are based on the following first order approximation |
---|
161 | \citep{Blanke1993}: |
---|
162 | \begin{equation} \label{Eq_tke_mxl0_1} |
---|
163 | l_k = l_\epsilon = \sqrt {2 \bar{e}\; } / N |
---|
164 | \end{equation} |
---|
165 | which is valid in a stable stratified region with constant values of the Brunt- |
---|
166 | Vais\"{a}l\"{a} frequency. The resulting length scale is bounded by the distance |
---|
167 | to the surface or to the bottom (\np{nn\_mxl}=0) or by the local vertical scale factor |
---|
168 | (\np{nn\_mxl}=1). \citet{Blanke1993} notice that this simplification has two major |
---|
169 | drawbacks: it makes no sense for locally unstable stratification and the |
---|
170 | computation no longer uses all the information contained in the vertical density |
---|
171 | profile. To overcome these drawbacks, \citet{Madec1998} introduces the |
---|
172 | \np{nn\_mxl}=2 or 3 cases, which add an extra assumption concerning the vertical |
---|
173 | gradient of the computed length scale. So, the length scales are first evaluated |
---|
174 | as in \eqref{Eq_tke_mxl0_1} and then bounded such that: |
---|
175 | \begin{equation} \label{Eq_tke_mxl_constraint} |
---|
176 | \frac{1}{e_3 }\left| {\frac{\partial l}{\partial k}} \right| \leq 1 |
---|
177 | \qquad \text{with }\ l = l_k = l_\epsilon |
---|
178 | \end{equation} |
---|
179 | \eqref{Eq_tke_mxl_constraint} means that the vertical variations of the length |
---|
180 | scale cannot be larger than the variations of depth. It provides a better |
---|
181 | approximation of the \citet{Gaspar1990} formulation while being much less |
---|
182 | time consuming. In particular, it allows the length scale to be limited not only |
---|
183 | by the distance to the surface or to the ocean bottom but also by the distance |
---|
184 | to a strongly stratified portion of the water column such as the thermocline |
---|
185 | (Fig.~\ref{Fig_mixing_length}). In order to impose the \eqref{Eq_tke_mxl_constraint} |
---|
186 | constraint, we introduce two additional length scales: $l_{up}$ and $l_{dwn}$, |
---|
187 | the upward and downward length scales, and evaluate the dissipation and |
---|
188 | mixing length scales as (and note that here we use numerical indexing): |
---|
189 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
190 | \begin{figure}[!t] \label{Fig_mixing_length} \begin{center} |
---|
191 | \includegraphics[width=1.00\textwidth]{./TexFiles/Figures/Fig_mixing_length.pdf} |
---|
192 | \caption {Illustration of the mixing length computation. } |
---|
193 | \end{center} |
---|
194 | \end{figure} |
---|
195 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
196 | \begin{equation} \label{Eq_tke_mxl2} |
---|
197 | \begin{aligned} |
---|
198 | l_{up\ \ }^{(k)} &= \min \left( l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3t}^{(k)}\ \ \ \; \right) |
---|
199 | \quad &\text{ from $k=1$ to $jpk$ }\ \\ |
---|
200 | l_{dwn}^{(k)} &= \min \left( l^{(k)} \ , \ l_{dwn}^{(k-1)} + e_{3t}^{(k-1)} \right) |
---|
201 | \quad &\text{ from $k=jpk$ to $1$ }\ \\ |
---|
202 | \end{aligned} |
---|
203 | \end{equation} |
---|
204 | where $l^{(k)}$ is computed using \eqref{Eq_tke_mxl0_1}, |
---|
205 | $i.e.$ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$. |
---|
206 | |
---|
207 | In the \np{nmxl}=2 case, the dissipation and mixing length scales take the same |
---|
208 | value: $ l_k= l_\epsilon = \min \left(\ l_{up} \;,\; l_{dwn}\ \right)$, while in the |
---|
209 | \np{nmxl}=2 case, the dissipation and mixing turbulent length scales are give |
---|
210 | as in \citet{Gaspar1990}: |
---|
211 | \begin{equation} \label{Eq_tke_mxl_gaspar} |
---|
212 | \begin{aligned} |
---|
213 | & l_k = \sqrt{\ l_{up} \ \ l_{dwn}\ } \\ |
---|
214 | & l_\epsilon = \min \left(\ l_{up} \;,\; l_{dwn}\ \right) |
---|
215 | \end{aligned} |
---|
216 | \end{equation} |
---|
217 | |
---|
218 | At the sea surface the value of $\bar{e}$ is prescribed from the wind |
---|
219 | stress field: $\bar{e}=rn\_ebb\;\left| \tau \right|$ (\np{rn\_ebb}=60 by default) |
---|
220 | with a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist |
---|
221 | parameters). Its value at the bottom of the ocean is assumed to be |
---|
222 | equal to the value of the level just above. The time integration of the |
---|
223 | $\bar{e}$ equation may formally lead to negative values because the |
---|
224 | numerical scheme does not ensure its positivity. To overcome this |
---|
225 | problem, a cut-off in the minimum value of $\bar{e}$ is used (\np{rn\_emin} |
---|
226 | namelist parameter). Following \citet{Gaspar1990}, the cut-off value is set |
---|
227 | to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$. This allows the subsequent formulations |
---|
228 | to match that of \citet{Gargett1984} for the diffusion in the thermocline and |
---|
229 | deep ocean : $(A^{vT} = 10^{-3} / N)$. |
---|
230 | In addition, a cut-off is applied on $A^{vm}$ and $A^{vT}$ to avoid numerical |
---|
231 | instabilities associated with too weak vertical diffusion. They must be |
---|
232 | specified at least larger than the molecular values, and are set through |
---|
233 | \np{avm0} and \np{avt0} (namelist parameters). |
---|
234 | |
---|
235 | % ------------------------------------------------------------------------------------------------------------- |
---|
236 | % TKE Turbulent Closure Scheme : new organization to energetic considerations |
---|
237 | % ------------------------------------------------------------------------------------------------------------- |
---|
238 | \subsection{TKE Turbulent Closure Scheme - time integration (\key{zdftke} and (\key{zdftke2})} |
---|
239 | \label{ZDF_tke2} |
---|
240 | |
---|
241 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
242 | \begin{figure}[!t] \label{Fig_TKE_time_scheme} \begin{center} |
---|
243 | \includegraphics[width=1.00\textwidth]{./TexFiles/Figures/Fig_ZDF_TKE_time_scheme.pdf} |
---|
244 | \caption {Illustration of the TKE time integration and its links to the momentum and tracer time integration. } |
---|
245 | \end{center} |
---|
246 | \end{figure} |
---|
247 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
248 | |
---|
249 | The production of turbulence by vertical shear (the first term of the right hand side |
---|
250 | of \eqref{Eq_zdftke_e}) should balance the loss of kinetic energy associated with |
---|
251 | the vertical momentum diffusion (first line in \eqref{Eq_PE_zdf}). |
---|
252 | The total loss of kinetic energy (in 1D for the demonstration) |
---|
253 | due to this term is obtained by multiply this quantity by $u^n$ and verticaly integrating: |
---|
254 | |
---|
255 | \begin{equation} \label{Eq_energ1} |
---|
256 | \int_{k_b}^{k_s} {u^t \frac{1}{e_3} |
---|
257 | \frac{\partial } |
---|
258 | {\partial k} \left( \frac{A^{vm}}{e_3} |
---|
259 | \frac{\partial{u^{t+1}}} |
---|
260 | {\partial k } \right) \; e_3 \; dk } |
---|
261 | = \left[ \frac{u^t}{e_3} A^{vm} \frac{\partial{u^{t+1}}}{\partial k} \right]_{k_b}^{k_s} |
---|
262 | - \int_{k_b}^{k_s}{\frac{A^{vm}}{{e_3}}\frac{\partial{u^t}}{\partial k}\frac{\partial{u^{t+1}}}{\partial k}} \ dk |
---|
263 | \end{equation} |
---|
264 | |
---|
265 | The first term of the right hand side of \eqref{Eq_energ1} represents the kinetic |
---|
266 | energy transfer at the surface (atmospheric forcing) and at the bottom (friction effect). |
---|
267 | The second term is always negative and have to balance the term of \eqref{Eq_zdftke_e} |
---|
268 | previously identified. |
---|
269 | |
---|
270 | The sink term (possibly a source term in statically unstable situations) of turbulence |
---|
271 | by buoyancy (second term of the right hand side of \eqref{Eq_zdftke_e}) must balance |
---|
272 | the source of potential energy associated with the vertical diffusion |
---|
273 | in the density equation (second line in \eqref{Eq_PE_zdf}). The source of potential |
---|
274 | energy (in 1D for the demonstration) due to this term is obtained by multiply this quantity |
---|
275 | by $gz{\rho_r}^{-1}$ and verticaly integrating: |
---|
276 | |
---|
277 | \begin{equation} \label{Eq_energ2} |
---|
278 | \begin{aligned} |
---|
279 | \int_{k_b}^{k_s}{\frac{g\;z}{e_3} \frac{\partial }{\partial k} |
---|
280 | \left( \frac{A^{vT}}{e_3} |
---|
281 | \frac{\partial{\rho^{t+1}}}{\partial k} \right)} \; e_3 \; dk |
---|
282 | =\left[ g\;z \frac{A^{vT}}{e_3} |
---|
283 | \frac{\partial{\rho^{t+1}}}{\partial k} \right]_{k_b}^{k_s} |
---|
284 | - \int_{k_b}^{k_s}{ \frac{A^{vT}}{e_3} g \frac{\partial{\rho^{t+1}}}{\partial k}} \; dk\\ |
---|
285 | \\ |
---|
286 | = - \left[ z\,A^{vT} {N_{t+1}}^2 \right]_{k_b}^{k_s} |
---|
287 | + \int_{k_b}^{k_s}{ \rho^{t+1} \; A^{vT}{N_{t+1}}^2 \; e_3 \; dk }\\ |
---|
288 | \end{aligned} |
---|
289 | \end{equation} |
---|
290 | where $N^2_{t+1}$ is the Brunt-Vaissala frequency at $t+1$ |
---|
291 | and noting that $\frac{\partial z}{\partial k} = e_3$. |
---|
292 | |
---|
293 | The first term is always zero because the Brunt Vaissala frequency is set to zero at the |
---|
294 | surface and the bottom. The second term is of opposite sign than the buoyancy term |
---|
295 | identified previously. |
---|
296 | |
---|
297 | Under these energetic considerations, \eqref{Eq_zdftke_e} have to be written like this |
---|
298 | to be consistant: |
---|
299 | |
---|
300 | \begin{equation} \label{Eq_zdftke_ene} |
---|
301 | \frac{\partial \bar{e}}{\partial t} = |
---|
302 | \frac{A^{vm}}{{e_3}^2 }\;\left[ {\left( {\frac{\partial u^a}{\partial k}} \right)\left( {\frac{\partial u^n}{\partial k}} \right) |
---|
303 | +\left( {\frac{\partial v^a}{\partial k}} \right)\left( {\frac{\partial v^n}{\partial k}} \right) |
---|
304 | } \right]-A^{vT}\,{N_n}^2+... |
---|
305 | \end{equation} |
---|
306 | |
---|
307 | Note that during a time step, the equation \eqref{Eq_zdftke_e} is resolved before those |
---|
308 | of momentum and density. So, the indice "a" (after) become "n" (now) and the indice "n" |
---|
309 | (now) become "b" (before). |
---|
310 | |
---|
311 | Moreover, the vertical shear have to be multiply by the appropriate viscosity for numerical |
---|
312 | stability. Thus, the vertical shear at U-point have to be multiply by the viscosity avmu and |
---|
313 | the vertical shear at V-point have to be multiply by the viscosity avmv. Next, these two |
---|
314 | quantities are averaged to obtain a production term by vertical shear at W-point : |
---|
315 | |
---|
316 | \begin{equation} \label{Eq_zdftke_ene2} |
---|
317 | \frac{\partial \bar{e}}{\partial t} = |
---|
318 | \frac{1}{{e_3}^2 }\;\left[ { |
---|
319 | A^{vmu}\left({\frac{\partial u^a}{\partial k}} \right)\left( {\frac{\partial u^n}{\partial k}} \right) |
---|
320 | +A^{vmv}\left({\frac{\partial v^a}{\partial k}} \right)\left( {\frac{\partial v^n}{\partial k}} \right) |
---|
321 | } \right]-A^{vT}\,{N_n}^2+... |
---|
322 | \end{equation} |
---|
323 | |
---|
324 | The TKE equation is resolved before the mixing length, the viscosity and diffusivity. Two tabs |
---|
325 | are then declared : dissl (dissipation length) (Remark : it's not only the dissipation lenght, |
---|
326 | it's the root of the TKE divided by the dissipation lenght) and avmt (viscosity at the points T) |
---|
327 | used for the vertical diffusion of the TKE. |
---|
328 | |
---|
329 | This new organization needs also a reorganization of the routine step.F90 (controled by |
---|
330 | the key \key{ene\_cons}). The bigger change is the estimation of the Brunt-Vaissala |
---|
331 | frequency at "n" instead of "b". Moreover for energetic considerations, the call of tranxt.F90 |
---|
332 | is done after the density update. On the contrary, the density is updated with scalars fields |
---|
333 | filtered by the Asselin filter. |
---|
334 | |
---|
335 | This new organisation of the routine zdftke force to save five three dimensionnal tabs in |
---|
336 | the restart : avmu, avmv, avt, avmt and dissl are needed to calculate $e_n$. At the end |
---|
337 | of the run (time step = nitend), the alternative is to save only $e_n$ estimated at the |
---|
338 | following time step (nitend+1). The next run using this restart file, the mixing length |
---|
339 | and turbulents coefficients are directly calculated using $e_n$. It is the same thing |
---|
340 | for the intermediate restart. |
---|
341 | |
---|
342 | |
---|
343 | %%GM for figure of the time scheme: |
---|
344 | \begin{equation} |
---|
345 | \rho^{t+1} \; A^{vT}{N_{t+1}}^2 \; dk \\ |
---|
346 | \end{equation} |
---|
347 | |
---|
348 | %%end GM |
---|
349 | |
---|
350 | |
---|
351 | % ------------------------------------------------------------------------------------------------------------- |
---|
352 | % K Profile Parametrisation (KPP) |
---|
353 | % ------------------------------------------------------------------------------------------------------------- |
---|
354 | \subsection{K Profile Parametrisation (KPP) (\key{zdfkpp}) } |
---|
355 | \label{ZDF_kpp} |
---|
356 | |
---|
357 | %--------------------------------------------namkpp-------------------------------------------------------- |
---|
358 | \namdisplay{namzdf_kpp} |
---|
359 | %-------------------------------------------------------------------------------------------------------------- |
---|
360 | |
---|
361 | The KKP scheme has been implemented by J. Chanut ... |
---|
362 | |
---|
363 | \colorbox{yellow}{Add a description of KPP here.} |
---|
364 | |
---|
365 | |
---|
366 | % ================================================================ |
---|
367 | % Convection |
---|
368 | % ================================================================ |
---|
369 | \section{Convection} |
---|
370 | \label{ZDF_conv} |
---|
371 | |
---|
372 | %--------------------------------------------namzdf-------------------------------------------------------- |
---|
373 | \namdisplay{namzdf} |
---|
374 | %-------------------------------------------------------------------------------------------------------------- |
---|
375 | |
---|
376 | Static instabilities (i.e. light potential densities under heavy ones) may |
---|
377 | occur at particular ocean grid points. In nature, convective processes |
---|
378 | quickly re-establish the static stability of the water column. These |
---|
379 | processes have been removed from the model via the hydrostatic |
---|
380 | assumption so they must be parameterized. Three parameterisations |
---|
381 | are available to deal with convective processes: a non-penetrative |
---|
382 | convective adjustment or an enhanced vertical diffusion, or/and the |
---|
383 | use of a turbulent closure scheme. |
---|
384 | |
---|
385 | % ------------------------------------------------------------------------------------------------------------- |
---|
386 | % Non-Penetrative Convective Adjustment |
---|
387 | % ------------------------------------------------------------------------------------------------------------- |
---|
388 | \subsection [Non-Penetrative Convective Adjustment (\np{ln\_tranpc}) ] |
---|
389 | {Non-Penetrative Convective Adjustment (\np{ln\_tranpc}=.true.) } |
---|
390 | \label{ZDF_npc} |
---|
391 | |
---|
392 | %--------------------------------------------namzdf-------------------------------------------------------- |
---|
393 | \namdisplay{namzdf} |
---|
394 | %-------------------------------------------------------------------------------------------------------------- |
---|
395 | |
---|
396 | |
---|
397 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
398 | \begin{figure}[!htb] \label{Fig_npc} \begin{center} |
---|
399 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_npc.pdf} |
---|
400 | \caption {Example of an unstable density profile treated by the non penetrative |
---|
401 | convective adjustment algorithm. $1^{st}$ step: the initial profile is checked from |
---|
402 | the surface to the bottom. It is found to be unstable between levels 3 and 4. |
---|
403 | They are mixed. The resulting $\rho$ is still larger than $\rho$(5): levels 3 to 5 |
---|
404 | are mixed. The resulting $\rho$ is still larger than $\rho$(6): levels 3 to 6 are |
---|
405 | mixed. The $1^{st}$ step ends since the density profile is then stable below |
---|
406 | the level 3. $2^{nd}$ step: the new $\rho$ profile is checked following the same |
---|
407 | procedure as in $1^{st}$ step: levels 2 to 5 are mixed. The new density profile |
---|
408 | is checked. It is found stable: end of algorithm.} |
---|
409 | \end{center} \end{figure} |
---|
410 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
411 | |
---|
412 | The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}=true. |
---|
413 | It is applied at each \np{nnpc1} time step and mixes downwards instantaneously |
---|
414 | the statically unstable portion of the water column, but only until the density |
---|
415 | structure becomes neutrally stable ($i.e.$ until the mixed portion of the water |
---|
416 | column has \textit{exactly} the density of the water just below) \citep{Madec_al_JPO91}. |
---|
417 | The associated algorithm is an iterative process used in the following way |
---|
418 | (Fig. \ref{Fig_npc}): starting from the top of the ocean, the first instability is |
---|
419 | found. Assume in the following that the instability is located between levels |
---|
420 | $k$ and $k+1$. The potential temperature and salinity in the two levels are |
---|
421 | vertically mixed, conserving the heat and salt contents of the water column. |
---|
422 | The new density is then computed by a linear approximation. If the new |
---|
423 | density profile is still unstable between levels $k+1$ and $k+2$, levels $k$, |
---|
424 | $k+1$ and $k+2$ are then mixed. This process is repeated until stability is |
---|
425 | established below the level $k$ (the mixing process can go down to the |
---|
426 | ocean bottom). The algorithm is repeated to check if the density profile |
---|
427 | between level $k-1$ and $k$ is unstable and/or if there is no deeper instability. |
---|
428 | |
---|
429 | This algorithm is significantly different from mixing statically unstable levels |
---|
430 | two by two. The latter procedure cannot converge with a finite number |
---|
431 | of iterations for some vertical profiles while the algorithm used in \NEMO |
---|
432 | converges for any profile in a number of iterations which is less than the |
---|
433 | number of vertical levels. This property is of paramount importance as |
---|
434 | pointed out by \citet{Killworth1989}: it avoids the existence of permanent |
---|
435 | and unrealistic static instabilities at the sea surface. This non-penetrative |
---|
436 | convective algorithm has been proved successful in studies of the deep |
---|
437 | water formation in the north-western Mediterranean Sea |
---|
438 | \citep{Madec_al_JPO91, Madec_al_DAO91, Madec_Crepon_Bk91}. |
---|
439 | |
---|
440 | Note that in the current implementation of this algorithm presents several |
---|
441 | limitations. First, potential density referenced to the sea surface is used to |
---|
442 | check whether the density profile is stable or not. This is a strong |
---|
443 | simplification which leads to large errors for realistic ocean simulations. |
---|
444 | Indeed, many water masses of the world ocean, especially Antarctic Bottom |
---|
445 | Water, are unstable when represented in surface-referenced potential density. |
---|
446 | The scheme will erroneously mix them up. Second, the mixing of potential |
---|
447 | density is assumed to be linear. This assures the convergence of the algorithm |
---|
448 | even when the equation of state is non-linear. Small static instabilities can thus |
---|
449 | persist due to cabbeling: they will be treated at the next time step. |
---|
450 | Third, temperature and salinity, and thus density, are mixed, but the |
---|
451 | corresponding velocity fields remain unchanged. When using a Richardson |
---|
452 | Number dependent eddy viscosity, the mixing of momentum is done through |
---|
453 | the vertical diffusion: after a static adjustment, the Richardson Number is zero |
---|
454 | and thus the eddy viscosity coefficient is at a maximum. When this convective |
---|
455 | adjustment algorithm is used with constant vertical eddy viscosity, spurious |
---|
456 | solutions can occur since the vertical momentum diffusion remains small even |
---|
457 | after a static adjustment. In that case, we recommend the addition of momentum |
---|
458 | mixing in a manner that mimics the mixing in temperature and salinity |
---|
459 | \citep{Speich_PhD92, Speich_al_JPO96}. |
---|
460 | |
---|
461 | % ------------------------------------------------------------------------------------------------------------- |
---|
462 | % Enhanced Vertical Diffusion |
---|
463 | % ------------------------------------------------------------------------------------------------------------- |
---|
464 | \subsection [Enhanced Vertical Diffusion (\np{ln\_zdfevd})] |
---|
465 | {Enhanced Vertical Diffusion (\np{ln\_zdfevd}=.true.)} |
---|
466 | \label{ZDF_evd} |
---|
467 | |
---|
468 | %--------------------------------------------namzdf-------------------------------------------------------- |
---|
469 | \namdisplay{namzdf} |
---|
470 | %-------------------------------------------------------------------------------------------------------------- |
---|
471 | |
---|
472 | The enhanced vertical diffusion parameterisation is used when \np{ln\_zdfevd}=true. |
---|
473 | In this case, the vertical eddy mixing coefficients are assigned very large values |
---|
474 | (a typical value is $10\;m^2s^{-1})$ in regions where the stratification is unstable |
---|
475 | ($i.e.$ when the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{Lazar_PhD97, Lazar_al_JPO99}. |
---|
476 | This is done either on tracers only (\np{n\_evdm}=0) or on both momentum and |
---|
477 | tracers (\np{n\_evdm}=1). |
---|
478 | |
---|
479 | In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and |
---|
480 | if \np{n\_evdm}=1, the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ |
---|
481 | values also, are set equal to the namelist parameter \np{avevd}. A typical value |
---|
482 | for $avevd$ is between 1 and $100~m^2.s^{-1}$. This parameterisation of |
---|
483 | convective processes is less time consuming than the convective adjustment |
---|
484 | algorithm presented above when mixing both tracers and momentum in the |
---|
485 | case of static instabilities. It requires the use of an implicit time stepping on |
---|
486 | vertical diffusion terms (i.e. \np{ln\_zdfexp}=false). |
---|
487 | |
---|
488 | % ------------------------------------------------------------------------------------------------------------- |
---|
489 | % Turbulent Closure Scheme |
---|
490 | % ------------------------------------------------------------------------------------------------------------- |
---|
491 | \subsection{Turbulent Closure Scheme (\key{zdftke})} |
---|
492 | \label{ZDF_tcs} |
---|
493 | |
---|
494 | The TKE turbulent closure scheme presented in \S\ref{ZDF_tke} and used |
---|
495 | when the \key{zdftke} is defined, in theory solves the problem of statically |
---|
496 | unstable density profiles. In such a case, the term corresponding to the |
---|
497 | destruction of turbulent kinetic energy through stratification in \eqref{Eq_zdftke_e} |
---|
498 | becomes a source term, since $N^2$ is negative. It results in large values of |
---|
499 | $A_T^{vT}$ and $A_T^{vT}$, and also the four neighbouring |
---|
500 | $A_u^{vm} {and}\;A_v^{vm}$ (up to $1\;m^2s^{-1})$. These large values |
---|
501 | restore the static stability of the water column in a way similar to that of the |
---|
502 | enhanced vertical diffusion parameterisation (\S\ref{ZDF_evd}). However, |
---|
503 | in the vicinity of the sea surface (first ocean layer), the eddy coefficients |
---|
504 | computed by the turbulent closure scheme do not usually exceed $10^{-2}m.s^{-1}$, |
---|
505 | because the mixing length scale is bounded by the distance to the sea surface |
---|
506 | (see \S\ref{ZDF_tke}). It can thus be useful to combine the enhanced vertical |
---|
507 | diffusion with the turbulent closure scheme, $i.e.$ setting the \np{ln\_zdfnpc} |
---|
508 | namelist parameter to true and defining the \key{zdftke} CPP key all together. |
---|
509 | |
---|
510 | The KPP turbulent closure scheme already includes enhanced vertical diffusion |
---|
511 | in the case of convection, as governed by the variables $bvsqcon$ and $difcon$ |
---|
512 | found in \mdl{zdfkpp}, therefore \np{np\_zdfevd} should not be used with the KPP |
---|
513 | scheme. %gm% + one word on non local flux with KPP scheme trakpp.F90 module... |
---|
514 | |
---|
515 | % ================================================================ |
---|
516 | % Double Diffusion Mixing |
---|
517 | % ================================================================ |
---|
518 | \section [Double Diffusion Mixing (\textit{zdfddm} - \key{zdfddm})] |
---|
519 | {Double Diffusion Mixing (\mdl{zdfddm} module - \key{zdfddm})} |
---|
520 | \label{ZDF_ddm} |
---|
521 | |
---|
522 | %-------------------------------------------namddm-------------------------------------------------------- |
---|
523 | \namdisplay{namzdf_ddm} |
---|
524 | %-------------------------------------------------------------------------------------------------------------- |
---|
525 | |
---|
526 | Double diffusion occurs when relatively warm, salty water overlies cooler, fresher |
---|
527 | water, or vice versa. The former condition leads to salt fingering and the latter |
---|
528 | to diffusive convection. Double-diffusive phenomena contribute to diapycnal |
---|
529 | mixing in extensive regions of the ocean. \citet{Merryfield1999} include a |
---|
530 | parameterisation of such phenomena in a global ocean model and show that |
---|
531 | it leads to relatively minor changes in circulation but exerts significant regional |
---|
532 | influences on temperature and salinity. |
---|
533 | |
---|
534 | Diapycnal mixing of S and T are described by diapycnal diffusion coefficients |
---|
535 | \begin{align*} % \label{Eq_zdfddm_Kz} |
---|
536 | &A^{vT} = A_o^{vT}+A_f^{vT}+A_d^{vT} \\ |
---|
537 | &A^{vS} = A_o^{vS}+A_f^{vS}+A_d^{vS} |
---|
538 | \end{align*} |
---|
539 | where subscript $f$ represents mixing by salt fingering, $d$ by diffusive convection, |
---|
540 | and $o$ by processes other than double diffusion. The rates of double-diffusive mixing depend on the buoyancy ratio $R_\rho = \alpha \partial_z T / \beta \partial_z S$, |
---|
541 | where $\alpha$ and $\beta$ are coefficients of thermal expansion and saline |
---|
542 | contraction (see \S\ref{TRA_eos}). To represent mixing of $S$ and $T$ by salt |
---|
543 | fingering, we adopt the diapycnal diffusivities suggested by Schmitt (1981): |
---|
544 | \begin{align} \label{Eq_zdfddm_f} |
---|
545 | A_f^{vS} &= \begin{cases} |
---|
546 | \frac{A^{\ast v}}{1+(R_\rho / R_c)^n } &\text{if $R_\rho > 1$ and $N^2>0$ } \\ |
---|
547 | 0 &\text{otherwise} |
---|
548 | \end{cases} |
---|
549 | \\ \label{Eq_zdfddm_f_T} |
---|
550 | A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho |
---|
551 | \end{align} |
---|
552 | |
---|
553 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
554 | \begin{figure}[!t] \label{Fig_zdfddm} \begin{center} |
---|
555 | \includegraphics[width=0.99\textwidth]{./TexFiles/Figures/Fig_zdfddm.pdf} |
---|
556 | \caption {From \citet{Merryfield1999} : (a) Diapycnal diffusivities $A_f^{vT}$ |
---|
557 | and $A_f^{vS}$ for temperature and salt in regions of salt fingering. Heavy |
---|
558 | curves denote $A^{\ast v} = 10^{-3}~m^2.s^{-1}$ and thin curves |
---|
559 | $A^{\ast v} = 10^{-4}~m^2.s^{-1}$ ; (b) diapycnal diffusivities $A_d^{vT}$ and |
---|
560 | $A_d^{vS}$ for temperature and salt in regions of diffusive convection. Heavy |
---|
561 | curves denote the Federov parameterisation and thin curves the Kelley |
---|
562 | parameterisation. The latter is not implemented in \NEMO. } |
---|
563 | \end{center} \end{figure} |
---|
564 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
565 | |
---|
566 | The factor 0.7 in \eqref{Eq_zdfddm_f_T} reflects the measured ratio |
---|
567 | $\alpha F_T /\beta F_S \approx 0.7$ of buoyancy flux of heat to buoyancy |
---|
568 | flux of salt ($e.g.$, \citet{McDougall_Taylor_JMR84}). Following \citet{Merryfield1999}, |
---|
569 | we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$. |
---|
570 | |
---|
571 | To represent mixing of S and T by diffusive layering, the diapycnal diffusivities suggested by Federov (1988) is used: |
---|
572 | \begin{align} \label{Eq_zdfddm_d} |
---|
573 | A_d^{vT} &= \begin{cases} |
---|
574 | 1.3635 \, \exp{\left( 4.6\, \exp{ \left[ -0.54\,( R_{\rho}^{-1} - 1 ) \right] } \right)} |
---|
575 | &\text{if $0<R_\rho < 1$ and $N^2>0$ } \\ |
---|
576 | 0 &\text{otherwise} |
---|
577 | \end{cases} |
---|
578 | \\ \label{Eq_zdfddm_d_S} |
---|
579 | A_d^{vS} &= \begin{cases} |
---|
580 | A_d^{vT}\ \left( 1.85\,R_{\rho} - 0.85 \right) |
---|
581 | &\text{if $0.5 \leq R_\rho<1$ and $N^2>0$ } \\ |
---|
582 | A_d^{vT} \ 0.15 \ R_\rho |
---|
583 | &\text{if $\ \ 0 < R_\rho<0.5$ and $N^2>0$ } \\ |
---|
584 | 0 &\text{otherwise} |
---|
585 | \end{cases} |
---|
586 | \end{align} |
---|
587 | |
---|
588 | The dependencies of \eqref{Eq_zdfddm_f} to \eqref{Eq_zdfddm_d_S} on $R_\rho$ |
---|
589 | are illustrated in Fig.~\ref{Fig_zdfddm}. Implementing this requires computing |
---|
590 | $R_\rho$ at each grid point on every time step. This is done in \mdl{eosbn2} at the |
---|
591 | same time as $N^2$ is computed. This avoids duplication in the computation of |
---|
592 | $\alpha$ and $\beta$ (which is usually quite expensive). |
---|
593 | |
---|
594 | % ================================================================ |
---|
595 | % Bottom Friction |
---|
596 | % ================================================================ |
---|
597 | \section [Bottom Friction (\textit{zdfbfr})] |
---|
598 | {Bottom Friction (\mdl{zdfbfr} module)} |
---|
599 | \label{ZDF_bfr} |
---|
600 | |
---|
601 | %--------------------------------------------nambfr-------------------------------------------------------- |
---|
602 | \namdisplay{namzdf_bfr} |
---|
603 | %-------------------------------------------------------------------------------------------------------------- |
---|
604 | |
---|
605 | Both the surface momentum flux (wind stress) and the bottom momentum |
---|
606 | flux (bottom friction) enter the equations as a condition on the vertical |
---|
607 | diffusive flux. For the bottom boundary layer, one has: |
---|
608 | \begin{equation} \label{Eq_zdfbfr_flux} |
---|
609 | A^{vm} \left( \partial \textbf{U}_h / \partial z \right) = \textbf{F}_h |
---|
610 | \end{equation} |
---|
611 | where $\textbf{F}_h$ is supposed to represent the horizontal momentum flux |
---|
612 | outside the logarithmic turbulent boundary layer (thickness of the order of |
---|
613 | 1~m in the ocean). How $\textbf{F}_h$ influences the interior depends on the |
---|
614 | vertical resolution of the model near the bottom relative to the Ekman layer |
---|
615 | depth. For example, in order to obtain an Ekman layer depth |
---|
616 | $d = \sqrt{2\;A^{vm}} / f = 50$~m, one needs a vertical diffusion coefficient |
---|
617 | $A^{vm} = 0.125$~m$^2$s$^{-1}$ (for a Coriolis frequency |
---|
618 | $f = 10^{-4}$~m$^2$s$^{-1}$). With a background diffusion coefficient |
---|
619 | $A^{vm} = 10^{-4}$~m$^2$s$^{-1}$, the Ekman layer depth is only 1.4~m. |
---|
620 | When the vertical mixing coefficient is this small, using a flux condition is |
---|
621 | equivalent to entering the viscous forces (either wind stress or bottom friction) |
---|
622 | as a body force over the depth of the top or bottom model layer. To illustrate |
---|
623 | this, consider the equation for $u$ at $k$, the last ocean level: |
---|
624 | \begin{equation} \label{Eq_zdfbfr_flux2} |
---|
625 | \frac{\partial u \; (k)}{\partial t} = \frac{1}{e_{3u}} \left[ A^{vm} \; (k) \frac{U \; (k-1) - U \; (k)}{e_{3uw} \; (k-1)} - F_u \right] \approx - \frac{F_u}{e_{3u}} |
---|
626 | \end{equation} |
---|
627 | For example, if the bottom layer thickness is 200~m, the Ekman transport will |
---|
628 | be distributed over that depth. On the other hand, if the vertical resolution |
---|
629 | is high (1~m or less) and a turbulent closure model is used, the turbulent |
---|
630 | Ekman layer will be represented explicitly by the model. However, the |
---|
631 | logarithmic layer is never represented in current primitive equation model |
---|
632 | applications: it is \emph{necessary} to parameterize the flux $\textbf{F}_h $. |
---|
633 | Two choices are available in \NEMO: a linear and a quadratic bottom friction. |
---|
634 | Note that in both cases, the rotation between the interior velocity and the |
---|
635 | bottom friction is neglected in the present release of \NEMO. |
---|
636 | |
---|
637 | % ------------------------------------------------------------------------------------------------------------- |
---|
638 | % Linear Bottom Friction |
---|
639 | % ------------------------------------------------------------------------------------------------------------- |
---|
640 | \subsection{Linear Bottom Friction (\np{nbotfr} = 1) } |
---|
641 | \label{ZDF_bfr_linear} |
---|
642 | |
---|
643 | The linear bottom friction parameterisation assumes that the bottom friction |
---|
644 | is proportional to the interior velocity (i.e. the velocity of the last model level): |
---|
645 | \begin{equation} \label{Eq_zdfbfr_linear} |
---|
646 | \textbf{F}_h = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b |
---|
647 | \end{equation} |
---|
648 | where $\textbf{U}_h^b$ is the horizontal velocity vector of the bottom ocean |
---|
649 | layer and $r$ is a friction coefficient expressed in m.s$^{-1}$. This coefficient |
---|
650 | is generally estimated by setting a typical decay time $\tau$ in the deep ocean, |
---|
651 | and setting $r = H / \tau$, where $H$ is the ocean depth. Commonly accepted |
---|
652 | values of $\tau$ are of the order of 100 to 200 days \citep{Weatherly_JMR84}. |
---|
653 | A value $\tau^{-1} = 10^{-7}$~s$^{-1}$ equivalent to 115 days, is usually used |
---|
654 | in quasi-geostrophic models. One may consider the linear friction as an |
---|
655 | approximation of quadratic friction, $r \approx 2\;C_D\;U_{av}$ (\citet{Gill1982}, |
---|
656 | Eq. 9.6.6). For example, with a drag coefficient $C_D = 0.002$, a typical speed |
---|
657 | of tidal currents of $U_{av} =0.1$~m.s$^{-1}$, and assuming an ocean depth |
---|
658 | $H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m.s$^{-1}$. |
---|
659 | This is the default value used in \NEMO. It corresponds to a decay time scale |
---|
660 | of 115~days. It can be changed by specifying \np{bfric1} (namelist parameter). |
---|
661 | |
---|
662 | In the code, the bottom friction is imposed by updating the value of the |
---|
663 | vertical eddy coefficient at the bottom level. Indeed, the discrete formulation |
---|
664 | of (\ref{Eq_zdfbfr_linear}) at the last ocean $T-$level, using the fact that |
---|
665 | $\textbf {U}_h =0$ below the ocean floor, leads to |
---|
666 | \begin{equation} \label{Eq_zdfbfr_linKz} |
---|
667 | \begin{split} |
---|
668 | A_u^{vm} &= r\;e_{3uw}\\ |
---|
669 | A_v^{vm} &= r\;e_{3vw}\\ |
---|
670 | \end{split} |
---|
671 | \end{equation} |
---|
672 | |
---|
673 | This update is done in \mdl{zdfbfr} when \np{nbotfr}=1. The value of $r$ |
---|
674 | used is \np{bfric1}. Setting \np{nbotfr}=3 is equivalent to setting $r=0$ and |
---|
675 | leads to a free-slip bottom boundary condition. Setting \np{nbotfr}=0 sets |
---|
676 | $r=2\;A_{vb}^{\rm {\bf U}} $, where $A_{vb}^{\rm {\bf U}} $ is the background |
---|
677 | vertical eddy coefficient, and a no-slip boundary condition is imposed. |
---|
678 | Note that this latter choice generally leads to an underestimation of the |
---|
679 | bottom friction: for example with a deepest level thickness of $200~m$ |
---|
680 | and $A_{vb}^{\rm {\bf U}} =10^{-4}$m$^2$.s$^{-1}$, the friction coefficient |
---|
681 | is only $r=10^{-6}$m.s$^{-1}$. |
---|
682 | |
---|
683 | % ------------------------------------------------------------------------------------------------------------- |
---|
684 | % Non-Linear Bottom Friction |
---|
685 | % ------------------------------------------------------------------------------------------------------------- |
---|
686 | \subsection{Non-Linear Bottom Friction (\np{nbotfr} = 2)} |
---|
687 | \label{ZDF_bfr_nonlinear} |
---|
688 | |
---|
689 | The non-linear bottom friction parameterisation assumes that the bottom |
---|
690 | friction is quadratic: |
---|
691 | \begin{equation} \label{Eq_zdfbfr_nonlinear} |
---|
692 | \textbf {F}_h = \frac{A^{vm}}{e_3 }\frac{\partial \textbf {U}_h |
---|
693 | }{\partial k}=C_D \;\sqrt {u_b ^2+v_b ^2+e_b } \;\; \textbf {U}_h^b |
---|
694 | \end{equation} |
---|
695 | |
---|
696 | with $\textbf{U}_h^b = (u_b\;,\;v_b)$ the horizontal interior velocity ($i.e.$ |
---|
697 | the horizontal velocity of the bottom ocean layer), $C_D$ a drag coefficient, |
---|
698 | and $e_b $ a bottom turbulent kinetic energy due to tides, internal waves |
---|
699 | breaking and other short time scale currents. A typical value of the drag |
---|
700 | coefficient is $C_D = 10^{-3} $. As an example, the CME experiment |
---|
701 | \citep{Treguier_JGR92} uses $C_D = 10^{-3}$ and $e_b = 2.5\;10^{-3}$m$^2$.s$^{-2}$, |
---|
702 | while the FRAM experiment \citep{Killworth1992} uses $e_b =0$ |
---|
703 | and $e_b =2.5\;\;10^{-3}$m$^2$.s$^{-2}$. The FRAM choices have been |
---|
704 | set as default values (\np{bfric2} and \np{bfeb2} namelist parameters). |
---|
705 | |
---|
706 | As for the linear case, the bottom friction is imposed in the code by |
---|
707 | updating the value of the vertical eddy coefficient at the bottom level: |
---|
708 | \begin{equation} \label{Eq_zdfbfr_nonlinKz} |
---|
709 | \begin{split} |
---|
710 | A_u^{vm} &=C_D\; e_{3uw} \left[ u^2 + \left(\bar{\bar{v}}^{i+1,j}\right)^2 + e_b \right]^ |
---|
711 | {1/2}\\ |
---|
712 | A_v^{vm} &=C_D\; e_{3uw} \left[ \left(\bar{\bar{u}}^{i,j+1}\right)^2 + v^2 + e_b \right]^{1/2}\\ |
---|
713 | \end{split} |
---|
714 | \end{equation} |
---|
715 | |
---|
716 | This update is done in \mdl{zdfbfr}. The coefficients that control the strength of the |
---|
717 | non-linear bottom friction are initialized as namelist parameters: $C_D$= \np{bfri2}, |
---|
718 | and $e_b$ =\np{bfeb2}. |
---|
719 | |
---|
720 | % ================================================================ |
---|