- Timestamp:
- 2010-04-12T16:59:59+02:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_ZDF.tex
r1225 r1831 7 7 8 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 \gmcomment{Steven remark (not taken into account : problem here with turbulent vs turbulence. I've changed "turbulent closure" to "turbulence closure", "turbulent mixing length" to "turbulence mixing length", but I've left "turbulent kinetic energy" alone - though I think it is an historical abberation! 10 Gurvan : I kept "turbulent closure etc "...} 9 10 11 \newpage 12 $\ $\newline % force a new ligne 11 13 12 14 … … 74 76 75 77 %--------------------------------------------namric--------------------------------------------------------- 76 \namdisplay{nam ric}78 \namdisplay{namzdf_ric} 77 79 %-------------------------------------------------------------------------------------------------------------- 78 80 … … 84 86 growth of Kelvin-Helmholtz like instabilities leads to a dependency between the 85 87 vertical eddy coefficients and the local Richardson number ($i.e.$ the 86 ratio of stratification to vertical shear). Following \citet{Pac Phil1981}, the following88 ratio of stratification to vertical shear). Following \citet{Pacanowski_Philander_JPO81}, the following 87 89 formulation has been implemented: 88 90 \begin{equation} \label{Eq_zdfric} … … 107 109 \label{ZDF_tke} 108 110 109 %--------------------------------------------nam tke---------------------------------------------------------110 \namdisplay{nam tke}111 %--------------------------------------------namzdf_tke-------------------------------------------------- 112 \namdisplay{namzdf_tke} 111 113 %-------------------------------------------------------------------------------------------------------------- 112 114 113 115 The vertical eddy viscosity and diffusivity coefficients are computed from a TKE 114 turbulent closure model based on a prognostic equation for $\bar 116 turbulent closure model based on a prognostic equation for $\bar{e}$, the turbulent 115 117 kinetic energy, and a closure assumption for the turbulent length scales. This 116 118 turbulent closure model has been developed by \citet{Bougeault1989} in the 117 119 atmospheric case, adapted by \citet{Gaspar1990} for the oceanic case, and 118 embedded in OPA by \citet{Blanke1993} for equatorial Atlantic simulations. Since119 then, significant modifications have been introduced by \citet{Madec1998} in both120 the implementation and the formulation of the mixing length scale. The time121 evolution of $\bar{e}$ is the result of the production of $\bar{e}$ through vertical122 shear, its destruction through stratification, its vertical diffusion, and its dissipation123 of \citet{Kolmogorov1942} type: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: 124 126 \begin{equation} \label{Eq_zdftke_e} 125 127 \frac{\partial \bar{e}}{\partial t} = 126 \frac{A^{vm}}{ e_3}\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2128 \frac{A^{vm}}{{e_3}^2 }\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2 127 129 +\left( {\frac{\partial v}{\partial k}} \right)^2} \right] 128 130 -A^{vT}\,N^2 … … 133 135 \begin{equation} \label{Eq_zdftke_kz} 134 136 \begin{split} 135 A^{vm} &= C_k\ l_k\ \sqrt {\bar{e} } \\137 A^{vm} &= C_k\ l_k\ \sqrt {\bar{e}\; } \\ 136 138 A^{vT} &= A^{vm} / P_{rt} 137 139 \end{split} … … 139 141 where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}), 140 142 $l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales, 141 $P_{rt} $ is the Prandtl number. The constants $C_k = \sqrt {2} /2$ and142 $C_\epsilon = 0.1$ are designed to deal with vertical mixing at any depth143 \citep{Gaspar1990}. They are set through namelist parameters \np{ediff}144 and \np{ediss}. $P_{rt} $ can be set to unity or, following \citet{Blanke1993},145 be a function of the local Richardson number, $R_i$: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$: 146 148 \begin{align*} \label{Eq_prt} 147 149 P_{rt} = \begin{cases} … … 151 153 \end{cases} 152 154 \end{align*} 153 Note that a horizontal Shapiro filter can optionally be applied to $R_i$. 154 However it is an obsolescent option that is not recommended. 155 The choice of $P_{rt} $ is controlled by the \np{npdl} namelist parameter. 155 The choice of $P_{rt}$ is controlled by the \np{nn\_pdl} namelist parameter. 156 156 157 157 For computational efficiency, the original formulation of the turbulent length 158 158 scales proposed by \citet{Gaspar1990} has been simplified. Four formulations 159 are proposed, the choice of which is controlled by the \np{n mxl} namelist159 are proposed, the choice of which is controlled by the \np{nn\_mxl} namelist 160 160 parameter. The first two are based on the following first order approximation 161 161 \citep{Blanke1993}: 162 162 \begin{equation} \label{Eq_tke_mxl0_1} 163 l_k = l_\epsilon = \sqrt {2 \bar e} / N164 \end{equation} 165 which is valid in a stable stratified region with constant values of the brunt-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 166 Vais\"{a}l\"{a} frequency. The resulting length scale is bounded by the distance 167 to the surface or to the bottom (\np{nmxl}=0) or by the local vertical scale factor (\np{nmxl}=1). \citet{Blanke1993} notice that this simplification has two major 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 168 169 drawbacks: it makes no sense for locally unstable stratification and the 169 170 computation no longer uses all the information contained in the vertical density 170 171 profile. To overcome these drawbacks, \citet{Madec1998} introduces the 171 \np{n mxl}=2 or 3 cases, which add an extra assumption concerning the vertical172 \np{nn\_mxl}=2 or 3 cases, which add an extra assumption concerning the vertical 172 173 gradient of the computed length scale. So, the length scales are first evaluated 173 174 as in \eqref{Eq_tke_mxl0_1} and then bounded such that: … … 185 186 constraint, we introduce two additional length scales: $l_{up}$ and $l_{dwn}$, 186 187 the upward and downward length scales, and evaluate the dissipation and 187 mixing turbulent length scales as (and note that here we use numerical 188 indexing): 188 mixing length scales as (and note that here we use numerical indexing): 189 189 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 190 190 \begin{figure}[!t] \label{Fig_mixing_length} \begin{center} … … 196 196 \begin{equation} \label{Eq_tke_mxl2} 197 197 \begin{aligned} 198 l_{up\ \ }^{(k)} &= \min \left( l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3 T}^{(k)}\ \ \ \; \right)198 l_{up\ \ }^{(k)} &= \min \left( l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3t}^{(k)}\ \ \ \; \right) 199 199 \quad &\text{ from $k=1$ to $jpk$ }\ \\ 200 l_{dwn}^{(k)} &= \min \left( l^{(k)} \ , \ l_{dwn}^{(k-1)} + e_{3 T}^{(k-1)} \right)200 l_{dwn}^{(k)} &= \min \left( l^{(k)} \ , \ l_{dwn}^{(k-1)} + e_{3t}^{(k-1)} \right) 201 201 \quad &\text{ from $k=jpk$ to $1$ }\ \\ 202 202 \end{aligned} 203 203 \end{equation} 204 204 where $l^{(k)}$ is computed using \eqref{Eq_tke_mxl0_1}, 205 $i.e.$ $l^{(k)} = \sqrt {2 \bar e^{(k)} / N^{(k)} }$.205 $i.e.$ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$. 206 206 207 207 In the \np{nmxl}=2 case, the dissipation and mixing length scales take the same 208 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 length scales are give209 \np{nmxl}=2 case, the dissipation and mixing turbulent length scales are give 210 210 as in \citet{Gaspar1990}: 211 211 \begin{equation} \label{Eq_tke_mxl_gaspar} … … 217 217 218 218 At the sea surface the value of $\bar{e}$ is prescribed from the wind 219 stress field: $\bar{e}= ebb\;\left| \tau \right|$ ($ebb=60$by default)220 with a minimal threshold of $emin0=10^{-4}~m^2.s^{-2}$ (namelist219 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 221 parameters). Its value at the bottom of the ocean is assumed to be 222 222 equal to the value of the level just above. The time integration of the 223 223 $\bar{e}$ equation may formally lead to negative values because the 224 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. Following 226 \citet{Gaspar1990}, the cut-off value is set to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$. 227 This allows the subsequent formulations to match that of\citet{Gargett1984} 228 for the diffusion in the thermocline and deep ocean : $(A^{vT} = 10^{-3} / N)$. 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)$. 229 230 In addition, a cut-off is applied on $A^{vm}$ and $A^{vT}$ to avoid numerical 230 231 instabilities associated with too weak vertical diffusion. They must be 231 232 specified at least larger than the molecular values, and are set through 232 \textit{avm0} and \textit{avt0} (\textbf{namelist} parameters). 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 233 350 234 351 % ------------------------------------------------------------------------------------------------------------- … … 239 356 240 357 %--------------------------------------------namkpp-------------------------------------------------------- 241 \namdisplay{namkpp} 242 %-------------------------------------------------------------------------------------------------------------- 243 244 The K-Profile Parametrization (KKP) developed by \cite{Large_al_RG94} has been 245 implemented in \NEMO by J. Chanut (PhD reference to be added here!). 358 \namdisplay{namzdf_kpp} 359 %-------------------------------------------------------------------------------------------------------------- 360 361 The KKP scheme has been implemented by J. Chanut ... 246 362 247 363 \colorbox{yellow}{Add a description of KPP here.} … … 274 390 \label{ZDF_npc} 275 391 276 %--------------------------------------------nam npc--------------------------------------------------------277 \namdisplay{nam npc}392 %--------------------------------------------namzdf-------------------------------------------------------- 393 \namdisplay{namzdf} 278 394 %-------------------------------------------------------------------------------------------------------------- 279 395 … … 298 414 the statically unstable portion of the water column, but only until the density 299 415 structure becomes neutrally stable ($i.e.$ until the mixed portion of the water 300 column has \textit{exactly} the density of the water just below) \citep{Madec 1991a}.416 column has \textit{exactly} the density of the water just below) \citep{Madec_al_JPO91}. 301 417 The associated algorithm is an iterative process used in the following way 302 418 (Fig. \ref{Fig_npc}): starting from the top of the ocean, the first instability is … … 320 436 convective algorithm has been proved successful in studies of the deep 321 437 water formation in the north-western Mediterranean Sea 322 \citep{Madec 1991a, Madec1991b, Madec1991c}.438 \citep{Madec_al_JPO91, Madec_al_DAO91, Madec_Crepon_Bk91}. 323 439 324 440 Note that in the current implementation of this algorithm presents several … … 341 457 after a static adjustment. In that case, we recommend the addition of momentum 342 458 mixing in a manner that mimics the mixing in temperature and salinity 343 \citep{Speich 1992, Speich1996}.459 \citep{Speich_PhD92, Speich_al_JPO96}. 344 460 345 461 % ------------------------------------------------------------------------------------------------------------- … … 357 473 In this case, the vertical eddy mixing coefficients are assigned very large values 358 474 (a typical value is $10\;m^2s^{-1})$ in regions where the stratification is unstable 359 ($i.e.$ when the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{Lazar 1997, Lazar1999}.475 ($i.e.$ when the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{Lazar_PhD97, Lazar_al_JPO99}. 360 476 This is done either on tracers only (\np{n\_evdm}=0) or on both momentum and 361 477 tracers (\np{n\_evdm}=1). … … 405 521 406 522 %-------------------------------------------namddm-------------------------------------------------------- 407 \namdisplay{nam ddm}523 \namdisplay{namzdf_ddm} 408 524 %-------------------------------------------------------------------------------------------------------------- 409 525 … … 422 538 \end{align*} 423 539 where subscript $f$ represents mixing by salt fingering, $d$ by diffusive convection, 424 and $o$ by processes other than double diffusion. The rates of double-diffusive mixing 425 depend on the buoyancy ratio $R_\rho = \alpha \partial_z T / \beta \partial_z S$, 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$, 426 541 where $\alpha$ and $\beta$ are coefficients of thermal expansion and saline 427 542 contraction (see \S\ref{TRA_eos}). To represent mixing of $S$ and $T$ by salt … … 454 569 we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$. 455 570 456 To represent mixing of S and T by diffusive layering, the diapycnal diffusivities suggested 457 by Federov (1988) is used: 571 To represent mixing of S and T by diffusive layering, the diapycnal diffusivities suggested by Federov (1988) is used: 458 572 \begin{align} \label{Eq_zdfddm_d} 459 573 A_d^{vT} &= \begin{cases} … … 486 600 487 601 %--------------------------------------------nambfr-------------------------------------------------------- 488 \namdisplay{nam bfr}602 \namdisplay{namzdf_bfr} 489 603 %-------------------------------------------------------------------------------------------------------------- 490 604 … … 536 650 is generally estimated by setting a typical decay time $\tau$ in the deep ocean, 537 651 and setting $r = H / \tau$, where $H$ is the ocean depth. Commonly accepted 538 values of $\tau$ are of the order of 100 to 200 days \citep{Weatherly 1984}.652 values of $\tau$ are of the order of 100 to 200 days \citep{Weatherly_JMR84}. 539 653 A value $\tau^{-1} = 10^{-7}$~s$^{-1}$ equivalent to 115 days, is usually used 540 654 in quasi-geostrophic models. One may consider the linear friction as an … … 585 699 breaking and other short time scale currents. A typical value of the drag 586 700 coefficient is $C_D = 10^{-3} $. As an example, the CME experiment 587 \citep{Treguier 1992} uses $C_D = 10^{-3}$ and $e_b = 2.5\;10^{-3}$m$^2$.s$^{-2}$,701 \citep{Treguier_JGR92} uses $C_D = 10^{-3}$ and $e_b = 2.5\;10^{-3}$m$^2$.s$^{-2}$, 588 702 while the FRAM experiment \citep{Killworth1992} uses $e_b =0$ 589 703 and $e_b =2.5\;\;10^{-3}$m$^2$.s$^{-2}$. The FRAM choices have been
Note: See TracChangeset
for help on using the changeset viewer.
