Changeset 11422
- Timestamp:
- 2019-08-08T15:40:47+02:00 (3 years ago)
- Location:
- NEMO/branches/2019/fix_vvl_ticket1791
- Files:
-
- 33 deleted
- 139 edited
- 32 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/fix_vvl_ticket1791/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg
r10075 r11422 22 22 cn_exp = "Pacific" ! experience name 23 23 nn_it000 = 1 ! first time step 24 nn_itend = 75! last time step (std 5475)24 nn_itend = 80 ! last time step (std 5475) 25 25 nn_istate = 0 ! output the initial state (1) or not (0) 26 26 / … … 30 30 ln_linssh = .false. ! =T linear free surface ==>> model level are fixed in time 31 31 ! 32 rn_rdt = 5 760. ! time step for the dynamics and tracer32 rn_rdt = 5400. ! time step for the dynamics and tracer 33 33 / 34 34 !----------------------------------------------------------------------- -
NEMO/branches/2019/fix_vvl_ticket1791/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg
r10072 r11422 22 22 cn_exp = "Nordic" ! experience name 23 23 nn_it000 = 1 ! first time step 24 nn_itend = 3 00 ! last time step24 nn_itend = 320 ! last time step 25 25 nn_istate = 0 ! output the initial state (1) or not (0) 26 26 ln_clobber = .true. ! clobber (overwrite) an existing file … … 31 31 ln_linssh = .false. ! =T linear free surface ==>> model level are fixed in time 32 32 ! 33 rn_rdt = 1 440. ! time step for the dynamics (and tracer if nn_acc=0)33 rn_rdt = 1350. ! time step for the dynamics (and tracer if nn_acc=0) 34 34 / 35 35 !----------------------------------------------------------------------- -
NEMO/branches/2019/fix_vvl_ticket1791/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg
r10072 r11422 22 22 cn_exp = "Nordic" ! experience name 23 23 nn_it000 = 1 ! first time step 24 nn_itend = 9 00 ! last time step24 nn_itend = 960 ! last time step 25 25 nn_istate = 0 ! output the initial state (1) or not (0) 26 26 ln_clobber = .true. ! clobber (overwrite) an existing file … … 31 31 ln_linssh = .false. ! =T linear free surface ==>> model level are fixed in time 32 32 ! 33 rn_rdt = 4 80. ! time step for the dynamics (and tracer if nn_acc=0)33 rn_rdt = 450. ! time step for the dynamics (and tracer if nn_acc=0) 34 34 / 35 35 !----------------------------------------------------------------------- -
NEMO/branches/2019/fix_vvl_ticket1791/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg
r10075 r11422 22 22 cn_exp = "ORCA2" ! experience name 23 23 nn_it000 = 1 ! first time step 24 nn_itend = 75! last time step (std 5475)24 nn_itend = 80 ! last time step (std 5475) 25 25 nn_istate = 0 ! output the initial state (1) or not (0) 26 26 / … … 30 30 ln_linssh = .false. ! =T linear free surface ==>> model level are fixed in time 31 31 ! 32 rn_rdt = 5 760. ! time step for the dynamics and tracer32 rn_rdt = 5400. ! time step for the dynamics and tracer 33 33 / 34 34 !----------------------------------------------------------------------- -
NEMO/branches/2019/fix_vvl_ticket1791/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
r10190 r11422 22 22 cn_exp = "ORCA2" ! experience name 23 23 nn_it000 = 1 ! first time step 24 nn_itend = 5 475! last time step (std 5475)24 nn_itend = 5840 ! last time step (std 5475) 25 25 nn_istate = 0 ! output the initial state (1) or not (0) 26 26 / … … 28 28 &namdom ! time and space domain 29 29 !----------------------------------------------------------------------- 30 rn_rdt = 5 760. ! time step for the dynamics and tracer30 rn_rdt = 5400. ! time step for the dynamics and tracer 31 31 / 32 32 !----------------------------------------------------------------------- … … 75 75 &namsbc ! Surface Boundary Condition manager (default: NO selection) 76 76 !----------------------------------------------------------------------- 77 nn_fsbc = 3! frequency of SBC module call77 nn_fsbc = 4 ! frequency of SBC module call 78 78 ! (also = the frequency of sea-ice & iceberg model call) 79 79 ! Type of air-sea fluxes -
NEMO/branches/2019/fix_vvl_ticket1791/cfgs/SPITZ12/EXPREF/namelist_ice_cfg
r10911 r11422 102 102 &namdia ! Diagnostics 103 103 !------------------------------------------------------------------------------ 104 ln_icediachk = .true. ! check online the heat, mass & salt budgets (T) or not (F) 104 105 / -
NEMO/branches/2019/fix_vvl_ticket1791/doc
- Property svn:ignore deleted
-
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex
- Property svn:ignore
-
old new 1 *.aux 2 *.bbl 3 *.blg 4 *.dvi 5 *.fdb* 6 *.fls 7 *.idx 8 *.ilg 9 *.ind 10 *.log 11 *.maf 12 *.mtc* 13 *.out 14 *.pdf 15 *.toc 16 _minted-* 1 figures
-
- Property svn:ignore
-
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/HTML_latex2html.sh
r10146 r11422 1 1 #!/bin/bash 2 2 3 ./inc/clean.sh4 ./inc/build.sh3 #./inc/clean.sh 4 #./inc/build.sh 5 5 6 cd tex_main 7 sed -i -e 's#\\documentclass#%\\documentclass#' -e '/{document}/ s/^/%/' \ 8 ../tex_sub/*.tex 9 sed -i 's#\\subfile{#\\include{#g' \ 10 NEMO_manual.tex 11 latex2html -local_icons -no_footnode -split 4 -link 2 -mkdir -dir ../html_LaTeX2HTML \ 12 $* \ 13 NEMO_manual 14 sed -i -e 's#%\\documentclass#\\documentclass#' -e '/{document}/ s/^%//' \ 15 ../tex_sub/*.tex 16 sed -i 's#\\include{#\\subfile{#g' \ 17 NEMO_manual.tex 6 sed -i -e 's#utf8#latin1#' \ 7 -e 's#\[outputdir=../build\]{minted}#\[\]{minted}#' \ 8 -e '/graphicspath/ s#{../#{../../#g' \ 9 global/packages.tex 10 11 cd ./NEMO/main 12 sed -i -e 's#\\documentclass#%\\documentclass#' -e '/{document}/ s#^#%#' ../subfiles/*.tex 13 sed -i 's#\\subfile{#\\input{#' chapters.tex appendices.tex 14 15 #latex2html -noimages -local_icons -no_footnode -split 4 -link 2 -dir ../html_LaTeX2HTML $* NEMO_manual 16 latex2html -debug -noreuse -init_file ../../l2hconf.pm -local_icons -dir ../build/html NEMO_manual 17 18 sed -i -e 's#%\\documentclass#\\documentclass#' -e '/{document}/ s#^%##' ../subfiles/*.tex 19 sed -i 's#\\input{#\\subfile{#' chapters.tex appendices.tex 18 20 cd - 19 21 22 sed -i -e 's#latin1#utf8#' \ 23 -e 's#\[\]{minted}#\[outputdir=../build\]{minted}#' \ 24 -e '/graphicspath/ s#{../../#{../#g' \ 25 global/packages.tex 26 20 27 exit 0 -
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO
- Property svn:ignore deleted
-
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/main
- Property svn:ignore
-
old new 1 *.aux2 *.bbl3 *.blg4 *.dvi5 *.fdb*6 *.fls7 *.idx8 1 *.ilg 9 2 *.ind 10 *.log11 *.maf12 *.mtc*13 *.out14 *.pdf15 *.toc16 _minted-*
-
-
Property
svn:externals
set to
^/utils/dev/bibtool.rsc bibtool.rsc
- Property svn:ignore
-
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/annex_A.tex
r10442 r11422 10 10 11 11 \minitoc 12 13 \vfill 14 \begin{figure}[b] 15 \subsubsection*{Changes record} 16 \begin{tabular}{l||l|m{0.65\linewidth}} 17 Release & Author & Modifications \\ 18 {\em 4.0} & {\em Mike Bell} & {\em review} \\ 19 {\em 3.x} & {\em Gurvan Madec} & {\em original} \\ 20 \end{tabular} 21 \end{figure} 22 12 23 13 24 \newpage … … 29 40 \begin{equation} 30 41 \label{apdx:A_s_slope} 31 \sigma_1 =\frac{1}{e_1 } \;\left. {\frac{\partial z}{\partial i}} \right|_s42 \sigma_1 =\frac{1}{e_1 } \; \left. {\frac{\partial z}{\partial i}} \right|_s 32 43 \quad \text{and} \quad 33 \sigma_2 =\frac{1}{e_2 }\;\left. {\frac{\partial z}{\partial j}} \right|_s 34 \end{equation} 35 36 The chain rule to establish the model equations in the curvilinear $s-$coordinate system is: 44 \sigma_2 =\frac{1}{e_2 } \; \left. {\frac{\partial z}{\partial j}} \right|_s . 45 \end{equation} 46 47 The model fields (e.g. pressure $p$) can be viewed as functions of $(i,j,z,t)$ (e.g. $p(i,j,z,t)$) or as 48 functions of $(i,j,s,t)$ (e.g. $p(i,j,s,t)$). The symbol $\bullet$ will be used to represent any one of 49 these fields. Any ``infinitesimal'' change in $\bullet$ can be written in two forms: 50 \begin{equation} 51 \label{apdx:A_s_infin_changes} 52 \begin{aligned} 53 & \delta \bullet = \delta i \left. \frac{ \partial \bullet }{\partial i} \right|_{j,s,t} 54 + \delta j \left. \frac{ \partial \bullet }{\partial i} \right|_{i,s,t} 55 + \delta s \left. \frac{ \partial \bullet }{\partial s} \right|_{i,j,t} 56 + \delta t \left. \frac{ \partial \bullet }{\partial t} \right|_{i,j,s} , \\ 57 & \delta \bullet = \delta i \left. \frac{ \partial \bullet }{\partial i} \right|_{j,z,t} 58 + \delta j \left. \frac{ \partial \bullet }{\partial i} \right|_{i,z,t} 59 + \delta z \left. \frac{ \partial \bullet }{\partial z} \right|_{i,j,t} 60 + \delta t \left. \frac{ \partial \bullet }{\partial t} \right|_{i,j,z} . 61 \end{aligned} 62 \end{equation} 63 Using the first form and considering a change $\delta i$ with $j, z$ and $t$ held constant, shows that 64 \begin{equation} 65 \label{apdx:A_s_chain_rule} 66 \left. {\frac{\partial \bullet }{\partial i}} \right|_{j,z,t} = 67 \left. {\frac{\partial \bullet }{\partial i}} \right|_{j,s,t} 68 + \left. {\frac{\partial s }{\partial i}} \right|_{j,z,t} \; 69 \left. {\frac{\partial \bullet }{\partial s}} \right|_{i,j,t} . 70 \end{equation} 71 The term $\left. \partial s / \partial i \right|_{j,z,t}$ can be related to the slope of constant $s$ surfaces, 72 (\autoref{apdx:A_s_slope}), by applying the second of (\autoref{apdx:A_s_infin_changes}) with $\bullet$ set to 73 $s$ and $j, t$ held constant 74 \begin{equation} 75 \label{apdx:a_delta_s} 76 \delta s|_{j,t} = 77 \delta i \left. \frac{ \partial s }{\partial i} \right|_{j,z,t} 78 + \delta z \left. \frac{ \partial s }{\partial z} \right|_{i,j,t} . 79 \end{equation} 80 Choosing to look at a direction in the $(i,z)$ plane in which $\delta s = 0$ and using 81 (\autoref{apdx:A_s_slope}) we obtain 82 \begin{equation} 83 \left. \frac{ \partial s }{\partial i} \right|_{j,z,t} = 84 - \left. \frac{ \partial z }{\partial i} \right|_{j,s,t} \; 85 \left. \frac{ \partial s }{\partial z} \right|_{i,j,t} 86 = - \frac{e_1 }{e_3 }\sigma_1 . 87 \label{apdx:a_ds_di_z} 88 \end{equation} 89 Another identity, similar in form to (\autoref{apdx:a_ds_di_z}), can be derived 90 by choosing $\bullet$ to be $s$ and using the second form of (\autoref{apdx:A_s_infin_changes}) to consider 91 changes in which $i , j$ and $s$ are constant. This shows that 92 \begin{equation} 93 \label{apdx:A_w_in_s} 94 w_s = \left. \frac{ \partial z }{\partial t} \right|_{i,j,s} = 95 - \left. \frac{ \partial z }{\partial s} \right|_{i,j,t} 96 \left. \frac{ \partial s }{\partial t} \right|_{i,j,z} 97 = - e_3 \left. \frac{ \partial s }{\partial t} \right|_{i,j,z} . 98 \end{equation} 99 100 In what follows, for brevity, indication of the constancy of the $i, j$ and $t$ indices is 101 usually omitted. Using the arguments outlined above one can show that the chain rules needed to establish 102 the model equations in the curvilinear $s-$coordinate system are: 37 103 \begin{equation} 38 104 \label{apdx:A_s_chain_rule} 39 105 \begin{aligned} 40 106 &\left. {\frac{\partial \bullet }{\partial t}} \right|_z = 41 \left. {\frac{\partial \bullet }{\partial t}} \right|_s 42 -\frac{\partial \bullet }{\partial s}\;\frac{\partial s}{\partial t}\\107 \left. {\frac{\partial \bullet }{\partial t}} \right|_s 108 + \frac{\partial \bullet }{\partial s}\; \frac{\partial s}{\partial t} , \\ 43 109 &\left. {\frac{\partial \bullet }{\partial i}} \right|_z = 44 110 \left. {\frac{\partial \bullet }{\partial i}} \right|_s 45 -\frac{\partial \bullet }{\partial s}\;\frac{\partial s}{\partial i}=46 \left. {\frac{\partial \bullet }{\partial i}} \right|_s 47 -\frac{e_1 }{e_3 }\sigma_1 \frac{\partial \bullet }{\partial s} \\111 +\frac{\partial \bullet }{\partial s}\; \frac{\partial s}{\partial i}= 112 \left. {\frac{\partial \bullet }{\partial i}} \right|_s 113 -\frac{e_1 }{e_3 }\sigma_1 \frac{\partial \bullet }{\partial s} , \\ 48 114 &\left. {\frac{\partial \bullet }{\partial j}} \right|_z = 49 \left. {\frac{\partial \bullet }{\partial j}} \right|_s 50 -\frac{\partial \bullet }{\partial s}\;\frac{\partial s}{\partial j}=51 \left. {\frac{\partial \bullet }{\partial j}} \right|_s 52 - \frac{e_2 }{e_3 }\sigma_2 \frac{\partial \bullet }{\partial s} \\53 &\;\frac{\partial \bullet }{\partial z} \;\; = \frac{1}{e_3 }\frac{\partial \bullet }{\partial s} 115 \left. {\frac{\partial \bullet }{\partial j}} \right|_s 116 + \frac{\partial \bullet }{\partial s}\;\frac{\partial s}{\partial j}= 117 \left. {\frac{\partial \bullet }{\partial j}} \right|_s 118 - \frac{e_2 }{e_3 }\sigma_2 \frac{\partial \bullet }{\partial s} , \\ 119 &\;\frac{\partial \bullet }{\partial z} \;\; = \frac{1}{e_3 }\frac{\partial \bullet }{\partial s} . 54 120 \end{aligned} 55 121 \end{equation} 56 122 57 In particular applying the time derivative chain rule to $z$ provides the expression for $w_s$,58 the vertical velocity of the $s-$surfaces referenced to a fix z-coordinate:59 \begin{equation}60 \label{apdx:A_w_in_s}61 w_s = \left. \frac{\partial z }{\partial t} \right|_s62 = \frac{\partial z}{\partial s} \; \frac{\partial s}{\partial t}63 = e_3 \, \frac{\partial s}{\partial t}64 \end{equation}65 123 66 124 % ================================================================ … … 79 137 { 80 138 \begin{array}{*{20}l} 81 \nabla \cdot {\ rm {\bf U}}139 \nabla \cdot {\mathrm {\mathbf U}} 82 140 &= \frac{1}{e_1 \,e_2 } \left[ \left. {\frac{\partial (e_2 \,u)}{\partial i}} \right|_z 83 141 +\left. {\frac{\partial(e_1 \,v)}{\partial j}} \right|_z \right] … … 115 173 $, it becomes:} 116 174 % 117 \nabla \cdot {\ rm {\bf U}}175 \nabla \cdot {\mathrm {\mathbf U}} 118 176 & = \frac{1}{e_1 \,e_2 \,e_3 } \left[ 119 177 \left. \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s … … 131 189 \end{subequations} 132 190 133 Here, $w$ is the vertical velocity relative to the $z-$coordinate system. 134 Introducing the dia-surface velocity component, 135 $\omega $, defined as the volume flux across the moving $s$-surfaces per unit horizontal area: 191 Here, $w$ is the vertical velocity relative to the $z-$coordinate system. 192 Using the first form of (\autoref{apdx:A_s_infin_changes}) 193 and the definitions (\autoref{apdx:A_s_slope}) and (\autoref{apdx:A_w_in_s}) for $\sigma_1$, $\sigma_2$ and $w_s$, 194 one can show that the vertical velocity, $w_p$ of a point 195 moving with the horizontal velocity of the fluid along an $s$ surface is given by 196 \begin{equation} 197 \label{apdx:A_w_p} 198 \begin{split} 199 w_p = & \left. \frac{ \partial z }{\partial t} \right|_s 200 + \frac{u}{e_1} \left. \frac{ \partial z }{\partial i} \right|_s 201 + \frac{v}{e_2} \left. \frac{ \partial z }{\partial j} \right|_s \\ 202 = & w_s + u \sigma_1 + v \sigma_2 . 203 \end{split} 204 \end{equation} 205 The vertical velocity across this surface is denoted by 136 206 \begin{equation} 137 207 \label{apdx:A_w_s} 138 \omega = w - w_s - \sigma_1 \,u - \sigma_2 \,v \\ 139 \end{equation} 140 with $w_s$ given by \autoref{apdx:A_w_in_s}, 141 we obtain the expression for the divergence of the velocity in the curvilinear $s-$coordinate system: 142 \begin{subequations} 143 \begin{align*} 144 { 145 \begin{array}{*{20}l} 146 \nabla \cdot {\rm {\bf U}} 147 &= \frac{1}{e_1 \,e_2 \,e_3 } \left[ 208 \omega = w - w_p = w - ( w_s + \sigma_1 \,u + \sigma_2 \,v ) . 209 \end{equation} 210 Hence 211 \begin{equation} 212 \frac{1}{e_3 } \frac{\partial}{\partial s} \left[ w - u\;\sigma_1 - v\;\sigma_2 \right] = 213 \frac{1}{e_3 } \frac{\partial}{\partial s} \left[ \omega + w_s \right] = 214 \frac{1}{e_3 } \left[ \frac{\partial \omega}{\partial s} 215 + \left. \frac{ \partial }{\partial t} \right|_s \frac{\partial z}{\partial s} \right] = 216 \frac{1}{e_3 } \frac{\partial \omega}{\partial s} + \frac{1}{e_3 } \left. \frac{ \partial e_3}{\partial t} . \right|_s 217 \end{equation} 218 219 Using (\autoref{apdx:A_w_s}) in our expression for $\nabla \cdot {\mathrm {\mathbf U}}$ we obtain 220 our final expression for the divergence of the velocity in the curvilinear $s-$coordinate system: 221 \begin{equation} 222 \nabla \cdot {\mathrm {\mathbf U}} = 223 \frac{1}{e_1 \,e_2 \,e_3 } \left[ 148 224 \left. \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s 149 225 +\left. \frac{\partial (e_1 \,e_3 \,v)}{\partial j} \right|_s \right] 150 226 + \frac{1}{e_3 } \frac{\partial \omega }{\partial s} 151 + \frac{1}{e_3 } \frac{\partial w_s }{\partial s} \\ \\ 152 &= \frac{1}{e_1 \,e_2 \,e_3 } \left[ 153 \left. \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s 154 +\left. \frac{\partial (e_1 \,e_3 \,v)}{\partial j} \right|_s \right] 155 + \frac{1}{e_3 } \frac{\partial \omega }{\partial s} 156 + \frac{1}{e_3 } \frac{\partial}{\partial s} \left( e_3 \; \frac{\partial s}{\partial t} \right) \\ \\ 157 &= \frac{1}{e_1 \,e_2 \,e_3 } \left[ 158 \left. \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s 159 +\left. \frac{\partial (e_1 \,e_3 \,v)}{\partial j} \right|_s \right] 160 + \frac{1}{e_3 } \frac{\partial \omega }{\partial s} 161 + \frac{\partial}{\partial s} \frac{\partial s}{\partial t} 162 + \frac{1}{e_3 } \frac{\partial s}{\partial t} \frac{\partial e_3}{\partial s} \\ \\ 163 &= \frac{1}{e_1 \,e_2 \,e_3 } \left[ 164 \left. \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s 165 +\left. \frac{\partial (e_1 \,e_3 \,v)}{\partial j} \right|_s \right] 166 + \frac{1}{e_3 } \frac{\partial \omega }{\partial s} 167 + \frac{1}{e_3 } \frac{\partial e_3}{\partial t} 168 \end{array} 169 } 170 \end{align*} 171 \end{subequations} 227 + \frac{1}{e_3 } \left. \frac{\partial e_3}{\partial t} \right|_s . 228 \end{equation} 172 229 173 230 As a result, the continuity equation \autoref{eq:PE_continuity} in the $s-$coordinates is: … … 178 235 {\left. {\frac{\partial (e_2 \,e_3 \,u)}{\partial i}} \right|_s 179 236 + \left. {\frac{\partial (e_1 \,e_3 \,v)}{\partial j}} \right|_s } \right] 180 +\frac{1}{e_3 }\frac{\partial \omega }{\partial s} = 0 181 \end{equation} 182 A additional term has appeared that takeinto account237 +\frac{1}{e_3 }\frac{\partial \omega }{\partial s} = 0 . 238 \end{equation} 239 An additional term has appeared that takes into account 183 240 the contribution of the time variation of the vertical coordinate to the volume budget. 184 241 … … 210 267 + w \;\frac{\partial u}{\partial z} \\ \\ 211 268 &= \left. {\frac{\partial u }{\partial t}} \right|_z 212 - \left. \zeta \right|_z v 213 + \frac{1}{e_1 \,e_2 }\left[ { \left.{ \frac{\partial (e_2 \,v)}{\partial i} }\right|_z 269 - \frac{1}{e_1 \,e_2 }\left[ { \left.{ \frac{\partial (e_2 \,v)}{\partial i} }\right|_z 214 270 -\left.{ \frac{\partial (e_1 \,u)}{\partial j} }\right|_z } \right] \; v 215 271 + \frac{1}{2e_1} \left.{ \frac{\partial (u^2+v^2)}{\partial i} } \right|_z … … 230 286 } \\ \\ 231 287 &= \left. {\frac{\partial u }{\partial t}} \right|_z 232 +\left. \zeta \right|_s \;v288 - \left. \zeta \right|_s \;v 233 289 + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s \\ 234 290 &\qquad \qquad \qquad \quad 235 291 + \frac{w}{e_3 } \;\frac{\partial u}{\partial s} 236 -\left[ {\frac{\sigma_1 }{e_3 }\frac{\partial v}{\partial s}292 + \left[ {\frac{\sigma_1 }{e_3 }\frac{\partial v}{\partial s} 237 293 - \frac{\sigma_2 }{e_3 }\frac{\partial u}{\partial s}} \right]\;v 238 294 - \frac{\sigma_1 }{2e_3 }\frac{\partial (u^2+v^2)}{\partial s} \\ \\ 239 295 &= \left. {\frac{\partial u }{\partial t}} \right|_z 240 +\left. \zeta \right|_s \;v296 - \left. \zeta \right|_s \;v 241 297 + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s \\ 242 298 &\qquad \qquad \qquad \quad … … 245 301 - \sigma_1 u\frac{\partial u}{\partial s} - \sigma_1 v\frac{\partial v}{\partial s}} \right] \\ \\ 246 302 &= \left. {\frac{\partial u }{\partial t}} \right|_z 247 +\left. \zeta \right|_s \;v303 - \left. \zeta \right|_s \;v 248 304 + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s 249 305 + \frac{1}{e_3} \left[ w - \sigma_2 v - \sigma_1 u \right] 250 \; \frac{\partial u}{\partial s} \\306 \; \frac{\partial u}{\partial s} . \\ 251 307 % 252 \intertext{Introducing $\omega$, the dia- a-surface velocity given by (\autoref{apdx:A_w_s}) }308 \intertext{Introducing $\omega$, the dia-s-surface velocity given by (\autoref{apdx:A_w_s}) } 253 309 % 254 310 &= \left. {\frac{\partial u }{\partial t}} \right|_z 255 +\left. \zeta \right|_s \;v311 - \left. \zeta \right|_s \;v 256 312 + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s 257 + \frac{1}{e_3 } \left( \omega -w_s \right) \frac{\partial u}{\partial s} \\313 + \frac{1}{e_3 } \left( \omega + w_s \right) \frac{\partial u}{\partial s} \\ 258 314 \end{array} 259 315 } … … 266 322 { 267 323 \begin{array}{*{20}l} 268 w_s\;\frac{\partial u}{\partial s}269 = \frac{\partial s}{\partial t}\; \frac{\partial u }{\partial s}270 = \left. {\frac{\partial u }{\partial t}} \right|_s - \left. {\frac{\partial u }{\partial t}} \right|_z \ quad ,324 \frac{w_s}{e_3} \;\frac{\partial u}{\partial s} 325 = - \left. \frac{\partial s}{\partial t} \right|_z \; \frac{\partial u }{\partial s} 326 = \left. {\frac{\partial u }{\partial t}} \right|_s - \left. {\frac{\partial u }{\partial t}} \right|_z \ . 271 327 \end{array} 272 328 } 273 329 \] 274 leads to the $s-$coordinate formulation of the total $z-$coordinate time derivative,330 This leads to the $s-$coordinate formulation of the total $z-$coordinate time derivative, 275 331 \ie the total $s-$coordinate time derivative : 276 332 \begin{align} … … 278 334 \left. \frac{D u}{D t} \right|_s 279 335 = \left. {\frac{\partial u }{\partial t}} \right|_s 280 +\left. \zeta \right|_s \;v336 - \left. \zeta \right|_s \;v 281 337 + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s 282 + \frac{1}{e_3 } \omega \;\frac{\partial u}{\partial s} 338 + \frac{1}{e_3 } \omega \;\frac{\partial u}{\partial s} . 283 339 \end{align} 284 340 Therefore, the vector invariant form of the total time derivative has exactly the same mathematical form in … … 306 362 + \frac{1}{e_3} \frac{\partial \omega}{\partial s} \right] \\ \\ 307 363 &&- \frac{v}{e_1 e_2 }\left( v \;\frac{\partial e_2 }{\partial i} 308 -u \;\frac{\partial e_1 }{\partial j} \right) \\364 -u \;\frac{\partial e_1 }{\partial j} \right) . \\ 309 365 \end{array} 310 366 } … … 328 384 - e_2 u \;\frac{\partial e_3 }{\partial i} 329 385 - e_1 v \;\frac{\partial e_3 }{\partial j} \right) 330 -\frac{1}{e_3} \frac{\partial \omega}{\partial s} \right] \\ \\386 + \frac{1}{e_3} \frac{\partial \omega}{\partial s} \right] \\ \\ 331 387 && - \frac{v}{e_1 e_2 }\left( v \;\frac{\partial e_2 }{\partial i} 332 388 -u \;\frac{\partial e_1 }{\partial j} \right) \\ \\ … … 337 393 && - \,u \left[ \frac{1}{e_1 e_2 e_3} \left( \frac{\partial(e_2 e_3 \, u)}{\partial i} 338 394 + \frac{\partial(e_1 e_3 \, v)}{\partial j} \right) 339 -\frac{1}{e_3} \frac{\partial \omega}{\partial s} \right]395 + \frac{1}{e_3} \frac{\partial \omega}{\partial s} \right] 340 396 - \frac{v}{e_1 e_2 }\left( v \;\frac{\partial e_2 }{\partial i} 341 -u \;\frac{\partial e_1 }{\partial j} \right) \\397 -u \;\frac{\partial e_1 }{\partial j} \right) . \\ 342 398 % 343 399 \intertext {Introducing a more compact form for the divergence of the momentum fluxes, … … 346 402 % 347 403 &= \left. {\frac{\partial u }{\partial t}} \right|_s 348 &+ \left. \nabla \cdot \left( {{\ rm {\bf U}}\,u} \right) \right|_s404 &+ \left. \nabla \cdot \left( {{\mathrm {\mathbf U}}\,u} \right) \right|_s 349 405 + \,u \frac{1}{e_3 } \frac{\partial e_3}{\partial t} 350 406 - \frac{v}{e_1 e_2 }\left( v \;\frac{\partial e_2 }{\partial i} … … 359 415 \label{apdx:A_sco_Dt_flux} 360 416 \left. \frac{D u}{D t} \right|_s = \frac{1}{e_3} \left. \frac{\partial ( e_3\,u)}{\partial t} \right|_s 361 + \left. \nabla \cdot \left( {{\ rm {\bf U}}\,u} \right) \right|_s417 + \left. \nabla \cdot \left( {{\mathrm {\mathbf U}}\,u} \right) \right|_s 362 418 - \frac{v}{e_1 e_2 }\left( v \;\frac{\partial e_2 }{\partial i} 363 -u \;\frac{\partial e_1 }{\partial j} \right) 419 -u \;\frac{\partial e_1 }{\partial j} \right). 364 420 \end{flalign} 365 421 which is the total time derivative expressed in the curvilinear $s-$coordinate system. … … 377 433 & =-\frac{1}{\rho_o e_1 }\left[ {\left. {\frac{\partial p}{\partial i}} \right|_s -\frac{e_1 }{e_3 }\sigma_1 \frac{\partial p}{\partial s}} \right] \\ 378 434 & =-\frac{1}{\rho_o \,e_1 }\left. {\frac{\partial p}{\partial i}} \right|_s +\frac{\sigma_1 }{\rho_o \,e_3 }\left( {-g\;\rho \;e_3 } \right) \\ 379 &=-\frac{1}{\rho_o \,e_1 }\left. {\frac{\partial p}{\partial i}} \right|_s -\frac{g\;\rho }{\rho_o }\sigma_1 435 &=-\frac{1}{\rho_o \,e_1 }\left. {\frac{\partial p}{\partial i}} \right|_s -\frac{g\;\rho }{\rho_o }\sigma_1 . 380 436 \end{split} 381 437 \] 382 438 Applying similar manipulation to the second component and 383 replacing $\sigma_1$ and $\sigma_2$ by their expression \autoref{apdx:A_s_slope}, it comes:439 replacing $\sigma_1$ and $\sigma_2$ by their expression \autoref{apdx:A_s_slope}, it becomes: 384 440 \begin{equation} 385 441 \label{apdx:A_grad_p_1} … … 391 447 -\frac{1}{\rho_o \, e_2 }\left. {\frac{\partial p}{\partial j}} \right|_z 392 448 &=-\frac{1}{\rho_o \,e_2 } \left( \left. {\frac{\partial p}{\partial j}} \right|_s 393 + g\;\rho \;\left. {\frac{\partial z}{\partial j}} \right|_s \right) \\449 + g\;\rho \;\left. {\frac{\partial z}{\partial j}} \right|_s \right) . \\ 394 450 \end{split} 395 451 \end{equation} … … 399 455 400 456 As in $z$-coordinate, 401 the horizontal pressure gradient can be split in two parts following \citet{ Marsaleix_al_OM08}.457 the horizontal pressure gradient can be split in two parts following \citet{marsaleix.auclair.ea_OM08}. 402 458 Let defined a density anomaly, $d$, by $d=(\rho - \rho_o)/ \rho_o$, 403 459 and a hydrostatic pressure anomaly, $p_h'$, by $p_h'= g \; \int_z^\eta d \; e_3 \; dk$. … … 405 461 \[ 406 462 \begin{split} 407 p &= g\; \int_z^\eta \rho \; e_3 \; dk = g\; \int_z^\eta \ left( \rho_o \,d + 1 \right) \; e_3 \; dk \\408 &= g \, \rho_o \; \int_z^\eta d \; e_3 \; dk + g \, \int_z^\eta e_3 \; dk463 p &= g\; \int_z^\eta \rho \; e_3 \; dk = g\; \int_z^\eta \rho_o \left( d + 1 \right) \; e_3 \; dk \\ 464 &= g \, \rho_o \; \int_z^\eta d \; e_3 \; dk + \rho_o g \, \int_z^\eta e_3 \; dk . 409 465 \end{split} 410 466 \] … … 412 468 \begin{equation} 413 469 \label{apdx:A_pressure} 414 p = \rho_o \; p_h' + g \, ( z + \eta)470 p = \rho_o \; p_h' + \rho_o \, g \, ( \eta - z ) 415 471 \end{equation} 416 472 and the hydrostatic pressure balance expressed in terms of $p_h'$ and $d$ is: 417 473 \[ 418 \frac{\partial p_h'}{\partial k} = - d \, g \, e_3 474 \frac{\partial p_h'}{\partial k} = - d \, g \, e_3 . 419 475 \] 420 476 421 477 Substituing \autoref{apdx:A_pressure} in \autoref{apdx:A_grad_p_1} and 422 using the definition of the density anomaly it comes theexpression in two parts:478 using the definition of the density anomaly it becomes an expression in two parts: 423 479 \begin{equation} 424 480 \label{apdx:A_grad_p_2} … … 426 482 -\frac{1}{\rho_o \, e_1 } \left. {\frac{\partial p}{\partial i}} \right|_z 427 483 &=-\frac{1}{e_1 } \left( \left. {\frac{\partial p_h'}{\partial i}} \right|_s 428 + g\; d \;\left. {\frac{\partial z}{\partial i}} \right|_s \right) - \frac{g}{e_1 } \frac{\partial \eta}{\partial i} \\484 + g\; d \;\left. {\frac{\partial z}{\partial i}} \right|_s \right) - \frac{g}{e_1 } \frac{\partial \eta}{\partial i} , \\ 429 485 % 430 486 -\frac{1}{\rho_o \, e_2 }\left. {\frac{\partial p}{\partial j}} \right|_z 431 487 &=-\frac{1}{e_2 } \left( \left. {\frac{\partial p_h'}{\partial j}} \right|_s 432 + g\; d \;\left. {\frac{\partial z}{\partial j}} \right|_s \right) - \frac{g}{e_2 } \frac{\partial \eta}{\partial j} \\488 + g\; d \;\left. {\frac{\partial z}{\partial j}} \right|_s \right) - \frac{g}{e_2 } \frac{\partial \eta}{\partial j} . \\ 433 489 \end{split} 434 490 \end{equation} … … 463 519 - \frac{1}{e_1 } \left( \frac{\partial p_h'}{\partial i} + g\; d \; \frac{\partial z}{\partial i} \right) 464 520 - \frac{g}{e_1 } \frac{\partial \eta}{\partial i} 465 + D_u^{\vect{U}} + F_u^{\vect{U}} 521 + D_u^{\vect{U}} + F_u^{\vect{U}} , 466 522 \end{multline} 467 523 \begin{multline} … … 473 529 - \frac{1}{e_2 } \left( \frac{\partial p_h'}{\partial j} + g\; d \; \frac{\partial z}{\partial j} \right) 474 530 - \frac{g}{e_2 } \frac{\partial \eta}{\partial j} 475 + D_v^{\vect{U}} + F_v^{\vect{U}} 531 + D_v^{\vect{U}} + F_v^{\vect{U}} . 476 532 \end{multline} 477 533 \end{subequations} … … 483 539 \label{apdx:A_PE_dyn_flux_u} 484 540 \frac{1}{e_3} \frac{\partial \left( e_3\,u \right) }{\partial t} = 485 \nabla \cdot \left( {{\rm {\bf U}}\,u} \right)541 - \nabla \cdot \left( {{\mathrm {\mathbf U}}\,u} \right) 486 542 + \left\{ {f + \frac{1}{e_1 e_2 }\left( v \;\frac{\partial e_2 }{\partial i} 487 543 -u \;\frac{\partial e_1 }{\partial j} \right)} \right\} \,v \\ 488 544 - \frac{1}{e_1 } \left( \frac{\partial p_h'}{\partial i} + g\; d \; \frac{\partial z}{\partial i} \right) 489 545 - \frac{g}{e_1 } \frac{\partial \eta}{\partial i} 490 + D_u^{\vect{U}} + F_u^{\vect{U}} 546 + D_u^{\vect{U}} + F_u^{\vect{U}} , 491 547 \end{multline} 492 548 \begin{multline} 493 549 \label{apdx:A_dyn_flux_v} 494 550 \frac{1}{e_3}\frac{\partial \left( e_3\,v \right) }{\partial t}= 495 - \nabla \cdot \left( {{\ rm {\bf U}}\,v} \right)496 +\left\{ {f + \frac{1}{e_1 e_2 }\left( v \;\frac{\partial e_2 }{\partial i}551 - \nabla \cdot \left( {{\mathrm {\mathbf U}}\,v} \right) 552 - \left\{ {f + \frac{1}{e_1 e_2 }\left( v \;\frac{\partial e_2 }{\partial i} 497 553 -u \;\frac{\partial e_1 }{\partial j} \right)} \right\} \,u \\ 498 554 - \frac{1}{e_2 } \left( \frac{\partial p_h'}{\partial j} + g\; d \; \frac{\partial z}{\partial j} \right) 499 555 - \frac{g}{e_2 } \frac{\partial \eta}{\partial j} 500 + D_v^{\vect{U}} + F_v^{\vect{U}} 556 + D_v^{\vect{U}} + F_v^{\vect{U}} . 501 557 \end{multline} 502 558 \end{subequations} … … 505 561 \begin{equation} 506 562 \label{apdx:A_dyn_zph} 507 \frac{\partial p_h'}{\partial k} = - d \, g \, e_3 563 \frac{\partial p_h'}{\partial k} = - d \, g \, e_3 . 508 564 \end{equation} 509 565 … … 531 587 \left[ \frac{\partial }{\partial i} \left( {e_2 \,e_3 \;Tu} \right) 532 588 + \frac{\partial }{\partial j} \left( {e_1 \,e_3 \;Tv} \right) \right] \\ 533 +\frac{1}{e_3} \frac{\partial }{\partial k} \left( Tw \right)589 - \frac{1}{e_3} \frac{\partial }{\partial k} \left( Tw \right) 534 590 + D^{T} +F^{T} 535 591 \end{multline} 536 592 537 The expression for the advection term is a straight consequence of ( A.4),593 The expression for the advection term is a straight consequence of (\autoref{apdx:A_sco_Continuity}), 538 594 the expression of the 3D divergence in the $s-$coordinates established above. 539 595 -
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/annex_B.tex
r10442 r11422 53 53 { 54 54 \begin{array}{*{20}l} 55 D^T=& \frac{1}{e_1\,e_2\,e_3 }\;\left[ {\ \ \ \ e_2\,e_3\,A^{lT} \;\left. 56 {\frac{\partial }{\partial i}\left( {\frac{1}{e_1}\;\left. {\frac{\partial T}{\partial i}} \right|_s -\frac{\sigma_1 }{e_3 }\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\ 57 &\qquad \quad \ \ \ +e_1\,e_3\,A^{lT} \;\left. {\frac{\partial }{\partial j}\left( {\frac{1}{e_2 }\;\left. {\frac{\partial T}{\partial j}} \right|_s -\frac{\sigma_2 }{e_3 }\;\frac{\partial T}{\partial s}} \right)} \right|_s \\ 58 &\qquad \quad \ \ \ + e_1\,e_2\,A^{lT} \;\frac{\partial }{\partial s}\left( {-\frac{\sigma_1 }{e_1 }\;\left. {\frac{\partial T}{\partial i}} \right|_s -\frac{\sigma_2 }{e_2 }\;\left. {\frac{\partial T}{\partial j}} \right|_s } \right. 59 \left. {\left. {+\left( {\varepsilon +\sigma_1^2+\sigma_2 ^2} \right)\;\frac{1}{e_3 }\;\frac{\partial T}{\partial s}} \right)\;\;} \right] 55 D^T= \frac{1}{e_1\,e_2\,e_3 } & \left\{ \quad \quad \frac{\partial }{\partial i} \left. \left[ e_2\,e_3 \, A^{lT} 56 \left( \ \frac{1}{e_1}\; \left. \frac{\partial T}{\partial i} \right|_s 57 -\frac{\sigma_1 }{e_3 } \; \frac{\partial T}{\partial s} \right) \right] \right|_s \right. \\ 58 & \quad \ + \ \left. \frac{\partial }{\partial j} \left. \left[ e_1\,e_3 \, A^{lT} 59 \left( \ \frac{1}{e_2 }\; \left. \frac{\partial T}{\partial j} \right|_s 60 -\frac{\sigma_2 }{e_3 } \; \frac{\partial T}{\partial s} \right) \right] \right|_s \right. \\ 61 & \quad \ + \ \left. e_1\,e_2\, \frac{\partial }{\partial s} \left[ A^{lT} \; \left( 62 -\frac{\sigma_1 }{e_1 } \; \left. \frac{\partial T}{\partial i} \right|_s 63 -\frac{\sigma_2 }{e_2 } \; \left. \frac{\partial T}{\partial j} \right|_s 64 +\left( \varepsilon +\sigma_1^2+\sigma_2 ^2 \right) \; \frac{1}{e_3 } \; \frac{\partial T}{\partial s} \right) \; \right] \; \right\} . 60 65 \end{array} 61 66 } 62 67 \end{align*} 63 68 64 Equation\autoref{apdx:B2} is obtained from \autoref{apdx:B1} without any additional assumption.69 \autoref{apdx:B2} is obtained from \autoref{apdx:B1} without any additional assumption. 65 70 Indeed, for the special case $k=z$ and thus $e_3 =1$, 66 71 we introduce an arbitrary vertical coordinate $s = s (i,j,z)$ as in \autoref{apdx:A} and … … 90 95 { 91 96 \begin{array}{*{20}l} 92 \intertext{Noting that $\frac{1}{e_1} \left. \frac{\partial e_3 }{\partial i} \right|_s = \frac{\partial \sigma_1 }{\partial s}$, itbecomes:}97 \intertext{Noting that $\frac{1}{e_1} \left. \frac{\partial e_3 }{\partial i} \right|_s = \frac{\partial \sigma_1 }{\partial s}$, this becomes:} 93 98 % 94 & =\frac{1}{e_1\,e_2\,e_3 }\left[ {\left. {\;\;\;\frac{\partial }{\partial i}\left( {\frac{e_2\,e_3 }{e_1}\,A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)} \right|_s \left. -\, {e_3 \frac{\partial }{\partial i}\left( {\frac{e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\99 D^T & =\frac{1}{e_1\,e_2\,e_3 }\left[ {\left. {\;\;\;\frac{\partial }{\partial i}\left( {\frac{e_2\,e_3 }{e_1}\,A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)} \right|_s \left. -\, {e_3 \frac{\partial }{\partial i}\left( {\frac{e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\ 95 100 & \qquad \qquad \quad -e_2 A^{lT}\;\frac{\partial \sigma_1 }{\partial s}\left. {\frac{\partial T}{\partial i}} \right|_s -e_1 \,\sigma_1 \frac{\partial }{\partial s}\left( {\frac{e_2 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right) \\ 96 & \qquad \qquad \quad\shoveright{ \left. { +e_1 \,\sigma_1 \frac{\partial }{\partial s}\left( {\frac{e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right)+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 }{e_3 }A^{vT}\;\frac{\partial T}{\partial z}} \right)\;\;\;} \right] }\\101 & \qquad \qquad \quad\shoveright{ \left. { +e_1 \,\sigma_1 \frac{\partial }{\partial s}\left( {\frac{e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right)+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 }{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right] }\\ 97 102 \\ 98 103 &=\frac{1}{e_1 \,e_2 \,e_3 } \left[ {\left. {\;\;\;\frac{\partial }{\partial i} \left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)} \right|_s \left. {-\frac{\partial }{\partial i}\left( {e_2 \,\sigma_1 A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\ 99 104 & \qquad \qquad \quad \left. {+\frac{e_2 \,\sigma_1 }{e_3}A^{lT}\;\frac{\partial T}{\partial s} \;\frac{\partial e_3 }{\partial i}} \right|_s -e_2 A^{lT}\;\frac{\partial \sigma_1 }{\partial s}\left. {\frac{\partial T}{\partial i}} \right|_s \\ 100 105 & \qquad \qquad \quad-e_2 \,\sigma_1 \frac{\partial}{\partial s}\left( {A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 \,\sigma_1 ^2}{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right) \\ 101 & \qquad \qquad \quad\shoveright{ \left. {-\frac{\partial \left( {e_1 \,e_2 \,\sigma_1 } \right)}{\partial s} \left( {\frac{\sigma_1 }{e_3}A^{lT}\;\frac{\partial T}{\partial s}} \right) + \frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 }{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right]} 106 & \qquad \qquad \quad\shoveright{ \left. {-\frac{\partial \left( {e_1 \,e_2 \,\sigma_1 } \right)}{\partial s} \left( {\frac{\sigma_1 }{e_3}A^{lT}\;\frac{\partial T}{\partial s}} \right) + \frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 }{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right]} . 102 107 \end{array} 103 108 } \\ … … 105 110 \begin{array}{*{20}l} 106 111 % 107 \intertext{ using the same remark as just above, itbecomes:}112 \intertext{Using the same remark as just above, $D^T$ becomes:} 108 113 % 109 &= \frac{1}{e_1 \,e_2 \,e_3 } \left[ {\left. {\;\;\;\frac{\partial }{\partial i} \left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s -e_2 \,\sigma_1 A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right.\;\;\; \\114 D^T &= \frac{1}{e_1 \,e_2 \,e_3 } \left[ {\left. {\;\;\;\frac{\partial }{\partial i} \left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s -e_2 \,\sigma_1 A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right.\;\;\; \\ 110 115 & \qquad \qquad \quad+\frac{e_1 \,e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}\;\frac{\partial \sigma_1 }{\partial s} - \frac {\sigma_1 }{e_3} A^{lT} \;\frac{\partial \left( {e_1 \,e_2 \,\sigma_1 } \right)}{\partial s}\;\frac{\partial T}{\partial s} \\ 111 116 & \qquad \qquad \quad-e_2 \left( {A^{lT}\;\frac{\partial \sigma_1 }{\partial s}\left. {\frac{\partial T}{\partial i}} \right|_s +\frac{\partial }{\partial s}\left( {\sigma_1 A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)-\frac{\partial \sigma_1 }{\partial s}\;A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right) \\ 112 & \qquad \qquad \quad\shoveright{\left. {+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 \,\sigma_1 ^2}{e_3 }A^{lT}\;\frac{\partial T}{\partial s}+\frac{e_1 \,e_2}{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right] }117 & \qquad \qquad \quad\shoveright{\left. {+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 \,\sigma_1 ^2}{e_3 }A^{lT}\;\frac{\partial T}{\partial s}+\frac{e_1 \,e_2}{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right] . } 113 118 \end{array} 114 119 } \\ … … 117 122 % 118 123 \intertext{Since the horizontal scale factors do not depend on the vertical coordinate, 119 the last term of the first line and the first term of the lastline cancel, while120 the second line reduces to a single vertical derivative, so it becomes:}124 the two terms on the second line cancel, while 125 the third line reduces to a single vertical derivative, so it becomes:} 121 126 % 122 & =\frac{1}{e_1 \,e_2 \,e_3 }\left[ {\left. {\;\;\;\frac{\partial }{\partial i}\left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s -e_2 \,\sigma_1 \,A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\127 D^T & =\frac{1}{e_1 \,e_2 \,e_3 }\left[ {\left. {\;\;\;\frac{\partial }{\partial i}\left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s -e_2 \,\sigma_1 \,A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\ 123 128 & \qquad \qquad \quad \shoveright{ \left. {+\frac{\partial }{\partial s}\left( {-e_2 \,\sigma_1 \,A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s +A^{lT}\frac{e_1 \,e_2 }{e_3 }\;\left( {\varepsilon +\sigma_1 ^2} \right)\frac{\partial T}{\partial s}} \right)\;\;\;} \right]} \\ 124 129 % 125 \intertext{ in other words, the horizontal/vertical Laplacian operator in the ($i$,$s$) plane takes the following form:}130 \intertext{In other words, the horizontal/vertical Laplacian operator in the ($i$,$s$) plane takes the following form:} 126 131 \end{array} 127 132 } \\ … … 162 167 the ($i$,$j$,$k$) curvilinear coordinate system in which 163 168 the equations of the ocean circulation model are formulated, 164 takes the following form \citep{ Redi_JPO82}:169 takes the following form \citep{redi_JPO82}: 165 170 166 171 \begin{equation} … … 169 174 \left[ {{ 170 175 \begin{array}{*{20}c} 171 {1+a_ 1 ^2} \hfill & {-a_1 a_2 } \hfill & {-a_1} \hfill \\172 {-a_1 a_2 } \hfill & {1+a_2 ^2} \hfill & {-a_2} \hfill \\173 {-a_1 } \hfill & {-a_2} \hfill & {\varepsilon +a_1 ^2+a_2 ^2} \hfill \\176 {1+a_2 ^2 +\varepsilon a_1 ^2} \hfill & {-a_1 a_2 (1-\varepsilon)} \hfill & {-a_1 (1-\varepsilon) } \hfill \\ 177 {-a_1 a_2 (1-\varepsilon) } \hfill & {1+a_1 ^2 +\varepsilon a_2 ^2} \hfill & {-a_2 (1-\varepsilon)} \hfill \\ 178 {-a_1 (1-\varepsilon)} \hfill & {-a_2 (1-\varepsilon)} \hfill & {\varepsilon +a_1 ^2+a_2 ^2} \hfill \\ 174 179 \end{array} 175 180 }} \right] 176 181 \end{equation} 177 where ($a_1$, $a_2$) are the isopycnal slopes in ($\textbf{i}$, $\textbf{j}$) directions, relative to geopotentials: 182 where ($a_1$, $a_2$) are $(-1) \times$ the isopycnal slopes in 183 ($\textbf{i}$, $\textbf{j}$) directions, relative to geopotentials (or 184 equivalently the slopes of the geopotential surfaces in the isopycnal 185 coordinate framework): 178 186 \[ 179 187 a_1 =\frac{e_3 }{e_1 }\left( {\frac{\partial \rho }{\partial i}} \right)\left( {\frac{\partial \rho }{\partial k}} \right)^{-1} … … 182 190 \right)\left( {\frac{\partial \rho }{\partial k}} \right)^{-1} 183 191 \] 184 185 In practice, isopycnal slopes are generally less than $10^{-2}$ in the ocean, 186 so $\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{Cox1987}: 192 and, as before, $\epsilon = A^{vT} / A^{lT}$. 193 194 In practice, $\epsilon$ is small and isopycnal slopes are generally less than $10^{-2}$ in the ocean, 195 so $\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{cox_OM87}. Keeping leading order terms\footnote{Apart from the (1,0) 196 and (0,1) elements which are set to zero. See \citet{griffies_bk04}, section 14.1.4.1 for a discussion of this point.}: 187 197 \begin{subequations} 188 198 \label{apdx:B4} … … 226 236 The isopycnal diffusion operator \autoref{apdx:B4}, 227 237 \autoref{apdx:B_ldfiso} conserves tracer quantity and dissipates its square. 228 The demonstration of the first property is trivial as \autoref{apdx:B4} is the divergence of fluxes. 229 Let us demonstrate the second one:238 As \autoref{apdx:B4} is the divergence of a flux, the demonstration of the first property is trivial, providing that the flux normal to the boundary is zero 239 (as it is when $A_h$ is zero at the boundary). Let us demonstrate the second one: 230 240 \[ 231 241 \iiint\limits_D T\;\nabla .\left( {\textbf{A}}_{\textbf{I}} \nabla T \right)\,dv … … 236 246 { 237 247 \begin{array}{*{20}l} 238 \nabla T\;.\left( {{\ rm {\bf A}}_{\rm {\bf I}} \nabla T}248 \nabla T\;.\left( {{\mathrm {\mathbf A}}_{\mathrm {\mathbf I}} \nabla T} 239 249 \right)&=A^{lT}\left[ {\left( {\frac{\partial T}{\partial i}} \right)^2-2a_1 240 250 \frac{\partial T}{\partial i}\frac{\partial T}{\partial k}+\left( … … 246 256 j}-a_2 \frac{\partial T}{\partial k}} \right)^2} 247 257 +\varepsilon \left(\frac{\partial T}{\partial k}\right) ^2\right] \\ 248 & \geq 0 258 & \geq 0 . 249 259 \end{array} 250 260 } … … 275 285 \end{equation} 276 286 277 where ($r_1$, $r_2$) are the isopycnal slopes in ($\textbf{i}$, $\textbf{j}$) directions, 278 relative to $s$-coordinate surfaces: 287 where ($r_1$, $r_2$) are $(-1)\times$ the isopycnal slopes in ($\textbf{i}$, $\textbf{j}$) directions, 288 relative to $s$-coordinate surfaces (or equivalently the slopes of the 289 $s$-coordinate surfaces in the isopycnal coordinate framework): 279 290 \[ 280 291 r_1 =\frac{e_3 }{e_1 }\left( {\frac{\partial \rho }{\partial i}} \right)\left( {\frac{\partial \rho }{\partial s}} \right)^{-1} … … 284 295 \] 285 296 286 To prove \autoref{apdx:B 5} by direct re-expression of \autoref{apdx:B_ldfiso} is straightforward, but laborious.297 To prove \autoref{apdx:B_ldfiso_s} by direct re-expression of \autoref{apdx:B_ldfiso} is straightforward, but laborious. 287 298 An easier way is first to note (by reversing the derivation of \autoref{sec:B_2} from \autoref{sec:B_1} ) that 288 299 the weak-slope operator may be \emph{exactly} reexpressed in non-orthogonal $i,j,\rho$-coordinates as … … 306 317 the (rotated,orthogonal) isoneutral axes to the non-orthogonal $i,j,\rho$-coordinates. 307 318 The further transformation into $i,j,s$-coordinates is exact, whatever the steepness of the $s$-surfaces, 308 in the same way as the transformation of horizontal/vertical Laplacian diffusion in $z$-coordinates ,319 in the same way as the transformation of horizontal/vertical Laplacian diffusion in $z$-coordinates in 309 320 \autoref{sec:B_1} onto $s$-coordinates is exact, however steep the $s$-surfaces. 310 321 … … 316 327 \label{sec:B_3} 317 328 318 The second order momentum diffusion operator (Laplacian) in the $z$-coordinateis found by329 The second order momentum diffusion operator (Laplacian) in $z$-coordinates is found by 319 330 applying \autoref{eq:PE_lap_vector}, the expression for the Laplacian of a vector, 320 331 to the horizontal velocity vector: … … 361 372 \end{align*} 362 373 Using \autoref{eq:PE_div}, the definition of the horizontal divergence, 363 the third compon ant of the second vector is obviously zero and thus :374 the third component of the second vector is obviously zero and thus : 364 375 \[ 365 \Delta {\textbf{U}}_h = \nabla _h \left( \chi \right) - \nabla _h \times \left( \zeta \ right) + \frac {1}{e_3 } \frac {\partial }{\partial k} \left( {\frac {1}{e_3 } \frac{\partial {\textbf{ U}}_h }{\partial k}} \right)376 \Delta {\textbf{U}}_h = \nabla _h \left( \chi \right) - \nabla _h \times \left( \zeta \textbf{k} \right) + \frac {1}{e_3 } \frac {\partial }{\partial k} \left( {\frac {1}{e_3 } \frac{\partial {\textbf{ U}}_h }{\partial k}} \right) . 366 377 \] 367 378 … … 379 390 - \nabla _h \times \left( {A^{lm}\;\zeta \;{\textbf{k}}} \right) 380 391 + \frac{1}{e_3 }\frac{\partial }{\partial k}\left( {\frac{A^{vm}\;}{e_3 } 381 \frac{\partial {\ rm {\bf U}}_h }{\partial k}} \right)\\392 \frac{\partial {\mathrm {\mathbf U}}_h }{\partial k}} \right) , \\ 382 393 \end{equation} 383 394 that is, in expanded form: … … 386 397 & = \frac{1}{e_1} \frac{\partial \left( {A^{lm}\chi } \right)}{\partial i} 387 398 -\frac{1}{e_2} \frac{\partial \left( {A^{lm}\zeta } \right)}{\partial j} 388 +\frac{1}{e_3} \frac{\partial u}{\partial k}\\399 +\frac{1}{e_3} \frac{\partial }{\partial k} \left( \frac{A^{vm}}{e_3} \frac{\partial u}{\partial k} \right) , \\ 389 400 D^{\textbf{U}}_v 390 401 & = \frac{1}{e_2 }\frac{\partial \left( {A^{lm}\chi } \right)}{\partial j} 391 402 +\frac{1}{e_1 }\frac{\partial \left( {A^{lm}\zeta } \right)}{\partial i} 392 +\frac{1}{e_3} \frac{\partial v}{\partial k}403 +\frac{1}{e_3} \frac{\partial }{\partial k} \left( \frac{A^{vm}}{e_3} \frac{\partial v}{\partial k} \right) . 393 404 \end{align*} 394 405 -
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/annex_E.tex
r10442 r11422 49 49 50 50 This results in a dissipatively dominant (\ie hyper-diffusive) truncation error 51 \citep{ Shchepetkin_McWilliams_OM05}.52 The overall performance of the advection scheme is similar to that reported in \cite{ Farrow1995}.51 \citep{shchepetkin.mcwilliams_OM05}. 52 The overall performance of the advection scheme is similar to that reported in \cite{farrow.stevens_JPO95}. 53 53 It is a relatively good compromise between accuracy and smoothness. 54 54 It is not a \emph{positive} scheme meaning false extrema are permitted but … … 65 65 the second term which is the diffusive part of the scheme, is evaluated using the \textit{before} velocity 66 66 (forward in time). 67 This is discussed by \citet{ Webb_al_JAOT98} in the context of the Quick advection scheme.67 This is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the Quick advection scheme. 68 68 UBS and QUICK schemes only differ by one coefficient. 69 Substituting 1/6 with 1/8 in (\autoref{eq:tra_adv_ubs}) leads to the QUICK advection scheme \citep{ Webb_al_JAOT98}.69 Substituting 1/6 with 1/8 in (\autoref{eq:tra_adv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}. 70 70 This option is not available through a namelist parameter, since the 1/6 coefficient is hard coded. 71 71 Nevertheless it is quite easy to make the substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme. … … 80 80 $\tau_w^{ubs}$ will be evaluated using either \textit{(a)} a centered $2^{nd}$ order scheme, 81 81 or \textit{(b)} a TVD scheme, or \textit{(c)} an interpolation based on conservative parabolic splines following 82 \citet{ Shchepetkin_McWilliams_OM05} implementation of UBS in ROMS, or \textit{(d)} an UBS.82 \citet{shchepetkin.mcwilliams_OM05} implementation of UBS in ROMS, or \textit{(d)} an UBS. 83 83 The $3^{rd}$ case has dispersion properties similar to an eight-order accurate conventional scheme. 84 84 … … 255 255 \subsection{Griffies iso-neutral diffusion operator} 256 256 257 Let try to define a scheme that get its inspiration from the \citet{ Griffies_al_JPO98} scheme,257 Let try to define a scheme that get its inspiration from the \citet{griffies.gnanadesikan.ea_JPO98} scheme, 258 258 but is formulated within the \NEMO framework 259 259 (\ie using scale factors rather than grid-size and having a position of $T$-points that … … 272 272 Nevertheless, this technique works fine for $T$ and $S$ as they are active tracers 273 273 (\ie they enter the computation of density), but it does not work for a passive tracer. 274 \citep{ Griffies_al_JPO98} introduce a different way to discretise the off-diagonal terms that274 \citep{griffies.gnanadesikan.ea_JPO98} introduce a different way to discretise the off-diagonal terms that 275 275 nicely solve the problem. 276 276 The idea is to get rid of combinations of an averaged in one direction combined with … … 308 308 \begin{figure}[!ht] 309 309 \begin{center} 310 \includegraphics[width= 0.70\textwidth]{Fig_ISO_triad}310 \includegraphics[width=\textwidth]{Fig_ISO_triad} 311 311 \caption{ 312 312 \protect\label{fig:ISO_triad} … … 508 508 \] 509 509 510 \citep{ Griffies_JPO98} introduces another way to implement the eddy induced advection, the so-called skew form.510 \citep{griffies_JPO98} introduces another way to implement the eddy induced advection, the so-called skew form. 511 511 It is based on a transformation of the advective fluxes using the non-divergent nature of the eddy induced velocity. 512 512 For example in the (\textbf{i},\textbf{k}) plane, the tracer advective fluxes can be transformed as follows: … … 574 574 The horizontal component reduces to the one use for an horizontal laplacian operator and 575 575 the vertical one keeps the same complexity, but not more. 576 This property has been used to reduce the computational time \citep{ Griffies_JPO98},576 This property has been used to reduce the computational time \citep{griffies_JPO98}, 577 577 but it is not of practical use as usually $A \neq A_e$. 578 578 Nevertheless this property can be used to choose a discret form of \autoref{eq:eiv_skew_continuous} which -
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/annex_iso.tex
r10442 r11422 4 4 \newcommand{\rML}[1][i]{\ensuremath{_{\mathrm{ML}\,#1}}} 5 5 \newcommand{\rMLt}[1][i]{\tilde{r}_{\mathrm{ML}\,#1}} 6 \newcommand{\triad}[6][]{\ensuremath{{}_{#2}^{#3}{\mathbb{#4}_{#1}}_{#5}^{\,#6}}} 6 %% Move to ../../global/new_cmds.tex to avoid error with \listoffigures 7 %\newcommand{\triad}[6][]{\ensuremath{{}_{#2}^{#3}{\mathbb{#4}_{#1}}_{#5}^{\,#6}} 7 8 \newcommand{\triadd}[5]{\ensuremath{{}_{#1}^{#2}{\mathbb{#3}}_{#4}^{\,#5}}} 8 9 \newcommand{\triadt}[5]{\ensuremath{{}_{#1}^{#2}{\tilde{\mathbb{#3}}}_{#4}^{\,#5}}} … … 52 53 the vertical skew flux is further reduced to ensure no vertical buoyancy flux, 53 54 giving an almost pure horizontal diffusive tracer flux within the mixed layer. 54 This is similar to the tapering suggested by \citet{ Gerdes1991}. See \autoref{subsec:Gerdes-taper}55 This is similar to the tapering suggested by \citet{gerdes.koberle.ea_CD91}. See \autoref{subsec:Gerdes-taper} 55 56 \item[\np{ln\_botmix\_triad}] 56 57 See \autoref{sec:iso_bdry}. … … 71 72 \label{sec:iso} 72 73 73 We have implemented into \NEMO a scheme inspired by \citet{ Griffies_al_JPO98},74 We have implemented into \NEMO a scheme inspired by \citet{griffies.gnanadesikan.ea_JPO98}, 74 75 but formulated within the \NEMO framework, using scale factors rather than grid-sizes. 75 76 … … 194 195 \subsection{Expression of the skew-flux in terms of triad slopes} 195 196 196 \citep{ Griffies_al_JPO98} introduce a different discretization of the off-diagonal terms that197 \citep{griffies.gnanadesikan.ea_JPO98} introduce a different discretization of the off-diagonal terms that 197 198 nicely solves the problem. 198 199 % Instead of multiplying the mean slope calculated at the $u$-point by … … 201 202 \begin{figure}[tb] 202 203 \begin{center} 203 \includegraphics[width= 1.05\textwidth]{Fig_GRIFF_triad_fluxes}204 \includegraphics[width=\textwidth]{Fig_GRIFF_triad_fluxes} 204 205 \caption{ 205 206 \protect\label{fig:ISO_triad} … … 265 266 \begin{figure}[tb] 266 267 \begin{center} 267 \includegraphics[width= 0.80\textwidth]{Fig_GRIFF_qcells}268 \includegraphics[width=\textwidth]{Fig_GRIFF_qcells} 268 269 \caption{ 269 270 \protect\label{fig:qcells} … … 473 474 474 475 To complete the discretization we now need only specify the triad volumes $_i^k\mathbb{V}_{i_p}^{k_p}$. 475 \citet{ Griffies_al_JPO98} identifies these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter cells,476 \citet{griffies.gnanadesikan.ea_JPO98} identifies these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter cells, 476 477 defined in terms of the distances between $T$, $u$,$f$ and $w$-points. 477 478 This is the natural discretization of \autoref{eq:cts-var}. … … 658 659 \begin{figure}[h] 659 660 \begin{center} 660 \includegraphics[width= 0.60\textwidth]{Fig_GRIFF_bdry_triads}661 \includegraphics[width=\textwidth]{Fig_GRIFF_bdry_triads} 661 662 \caption{ 662 663 \protect\label{fig:bdry_triads} … … 685 686 As discussed in \autoref{subsec:LDF_slp_iso}, 686 687 iso-neutral slopes relative to geopotentials must be bounded everywhere, 687 both for consistency with the small-slope approximation and for numerical stability \citep{ Cox1987, Griffies_Bk04}.688 both for consistency with the small-slope approximation and for numerical stability \citep{cox_OM87, griffies_bk04}. 688 689 The bound chosen in \NEMO is applied to each component of the slope separately and 689 690 has a value of $1/100$ in the ocean interior. … … 732 733 \[ 733 734 % \label{eq:iso_tensor_ML} 734 D^{lT}=\nabla {\ rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad735 D^{lT}=\nabla {\mathrm {\mathbf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 735 736 \mbox{with}\quad \;\;\Re =\left( {{ 736 737 \begin{array}{*{20}c} … … 829 830 (\eg the green triad $i_p=1/2,k_p=-1/2$) are tapered to the appropriate basal triad.} 830 831 % } 831 \includegraphics[width= 0.60\textwidth]{Fig_GRIFF_MLB_triads}832 \includegraphics[width=\textwidth]{Fig_GRIFF_MLB_triads} 832 833 \end{figure} 833 834 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 847 848 \[ 848 849 % \label{eq:iso_tensor_ML2} 849 D^{lT}=\nabla {\ rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad850 D^{lT}=\nabla {\mathrm {\mathbf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 850 851 \mbox{with}\quad \;\;\Re =\left( {{ 851 852 \begin{array}{*{20}c} … … 859 860 \footnote{ 860 861 To ensure good behaviour where horizontal density gradients are weak, 861 we in fact follow \citet{ Gerdes1991} and862 we in fact follow \citet{gerdes.koberle.ea_CD91} and 862 863 set $\rML^*=\mathrm{sgn}(\tilde{r}_i)\min(|\rMLt^2/\tilde{r}_i|,|\tilde{r}_i|)-\sigma_i$. 863 864 } … … 865 866 This approach is similar to multiplying the iso-neutral diffusion coefficient by 866 867 $\tilde{r}_{\mathrm{max}}^{-2}\tilde{r}_i^{-2}$ for steep slopes, 867 as suggested by \citet{ Gerdes1991} (see also \citet{Griffies_Bk04}).868 as suggested by \citet{gerdes.koberle.ea_CD91} (see also \citet{griffies_bk04}). 868 869 Again it is applied separately to each triad $_i^k\mathbb{R}_{i_p}^{k_p}$ 869 870 … … 925 926 926 927 However, when \np{ln\_traldf\_triad} is set true, 927 \NEMO instead implements eddy induced advection according to the so-called skew form \citep{ Griffies_JPO98}.928 \NEMO instead implements eddy induced advection according to the so-called skew form \citep{griffies_JPO98}. 928 929 It is based on a transformation of the advective fluxes using the non-divergent nature of the eddy induced velocity. 929 930 For example in the (\textbf{i},\textbf{k}) plane, … … 1139 1140 it is equivalent to a horizontal eiv (eddy-induced velocity) that is uniform within the mixed layer 1140 1141 \autoref{eq:eiv_v}. 1141 This ensures that the eiv velocities do not restratify the mixed layer \citep{ Treguier1997,Danabasoglu_al_2008}.1142 This ensures that the eiv velocities do not restratify the mixed layer \citep{treguier.held.ea_JPO97,danabasoglu.ferrari.ea_JC08}. 1142 1143 Equivantly, in terms of the skew-flux formulation we use here, 1143 1144 the linear slope tapering within the mixed-layer gives a linearly varying vertical flux, … … 1153 1154 $uw$ (integer +1/2 $i$, integer $j$, integer +1/2 $k$) and $vw$ (integer $i$, integer +1/2 $j$, integer +1/2 $k$) 1154 1155 points (see Table \autoref{tab:cell}) respectively. 1155 We follow \citep{ Griffies_Bk04} and calculate the streamfunction at a given $uw$-point from1156 We follow \citep{griffies_bk04} and calculate the streamfunction at a given $uw$-point from 1156 1157 the surrounding four triads according to: 1157 1158 \[ -
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_ASM.tex
r10442 r11422 8 8 \label{chap:ASM} 9 9 10 Authors: D. Lea, M. Martin, K. Mogensen, A. Weaver, ... % do we keep 10 \minitoc 11 11 12 \minitoc 12 \vfill 13 \begin{figure}[b] 14 \subsubsection*{Changes record} 15 \begin{tabular}{l||l|m{0.65\linewidth}} 16 Release & Author & Modifications \\ 17 {\em 4.0} & {\em D. J. Lea} & {\em \NEMO 4.0 updates} \\ 18 {\em 3.4} & {\em D. J. Lea, M. Martin, K. Mogensen, A. Weaver} & {\em Initial version} \\ 19 \end{tabular} 20 \end{figure} 13 21 14 22 \newpage 15 23 16 24 The ASM code adds the functionality to apply increments to the model variables: temperature, salinity, 17 sea surface height, velocity and sea ice concentration. 25 sea surface height, velocity and sea ice concentration. 18 26 These are read into the model from a NetCDF file which may be produced by separate data assimilation code. 19 27 The code can also output model background fields which are used as an input to data assimilation code. … … 37 45 it may be preferable to introduce the increment gradually into the ocean model in order to 38 46 minimize spurious adjustment processes. 39 This technique is referred to as Incremental Analysis Updates (IAU) \citep{ Bloom_al_MWR96}.47 This technique is referred to as Incremental Analysis Updates (IAU) \citep{bloom.takacs.ea_MWR96}. 40 48 IAU is a common technique used with 3D assimilation methods such as 3D-Var or OI. 41 49 IAU is used when \np{ln\_asmiau} is set to true. 42 50 43 With IAU, the model state trajectory ${\ bf x}$ in the assimilation window ($t_{0} \leq t_{i} \leq t_{N}$)51 With IAU, the model state trajectory ${\mathbf x}$ in the assimilation window ($t_{0} \leq t_{i} \leq t_{N}$) 44 52 is corrected by adding the analysis increments for temperature, salinity, horizontal velocity and SSH as 45 53 additional tendency terms to the prognostic equations: 46 54 \begin{align*} 47 55 % \label{eq:wa_traj_iau} 48 {\ bf x}^{a}(t_{i}) = M(t_{i}, t_{0})[{\bf x}^{b}(t_{0})] \; + \; F_{i} \delta \tilde{\bf x}^{a}56 {\mathbf x}^{a}(t_{i}) = M(t_{i}, t_{0})[{\mathbf x}^{b}(t_{0})] \; + \; F_{i} \delta \tilde{\mathbf x}^{a} 49 57 \end{align*} 50 where $F_{i}$ is a weighting function for applying the increments $\delta\tilde{\ bf x}^{a}$ defined such that58 where $F_{i}$ is a weighting function for applying the increments $\delta\tilde{\mathbf x}^{a}$ defined such that 51 59 $\sum_{i=1}^{N} F_{i}=1$. 52 ${\ bf x}^b$ denotes the model initial state and ${\bf x}^a$ is the model state after the increments are applied.60 ${\mathbf x}^b$ denotes the model initial state and ${\mathbf x}^a$ is the model state after the increments are applied. 53 61 To control the adjustment time of the model to the increment, 54 62 the increment can be applied over an arbitrary sub-window, $t_{m} \leq t_{i} \leq t_{n}$, … … 56 64 Typically the increments are spread evenly over the full window. 57 65 In addition, two different weighting functions have been implemented. 58 The first function employs constant weights,66 The first function (namelist option \np{niaufn} = 0) employs constant weights, 59 67 \begin{align} 60 68 \label{eq:F1_i} … … 62 70 =\left\{ 63 71 \begin{array}{ll} 64 0 & {\ rm if} \; \; \; t_{i} < t_{m} \\65 1/M & {\ rm if} \; \; \; t_{m} < t_{i} \leq t_{n} \\66 0 & {\ rm if} \; \; \; t_{i} > t_{n}72 0 & {\mathrm if} \; \; \; t_{i} < t_{m} \\ 73 1/M & {\mathrm if} \; \; \; t_{m} < t_{i} \leq t_{n} \\ 74 0 & {\mathrm if} \; \; \; t_{i} > t_{n} 67 75 \end{array} 68 \right. 76 \right. 69 77 \end{align} 70 78 where $M = m-n$. 71 The second function employs peaked hat-like weights in order to give maximum weight in the centre of the sub-window,79 The second function (namelist option \np{niaufn} = 1) employs peaked hat-like weights in order to give maximum weight in the centre of the sub-window, 72 80 with the weighting reduced linearly to a small value at the window end-points: 73 81 \begin{align} … … 76 84 =\left\{ 77 85 \begin{array}{ll} 78 0 & {\ rm if} \; \; \; t_{i} < t_{m} \\79 \alpha \, i & {\ rm if} \; \; \; t_{m} \leq t_{i} \leq t_{M/2} \\80 \alpha \, (M - i +1) & {\ rm if} \; \; \; t_{M/2} < t_{i} \leq t_{n} \\81 0 & {\ rm if} \; \; \; t_{i} > t_{n}86 0 & {\mathrm if} \; \; \; t_{i} < t_{m} \\ 87 \alpha \, i & {\mathrm if} \; \; \; t_{m} \leq t_{i} \leq t_{M/2} \\ 88 \alpha \, (M - i +1) & {\mathrm if} \; \; \; t_{M/2} < t_{i} \leq t_{n} \\ 89 0 & {\mathrm if} \; \; \; t_{i} > t_{n} 82 90 \end{array} 83 91 \right. 84 92 \end{align} 85 where $\alpha^{-1} = \sum_{i=1}^{M/2} 2i$ and $M$ is assumed to be even. 93 where $\alpha^{-1} = \sum_{i=1}^{M/2} 2i$ and $M$ is assumed to be even. 86 94 The weights described by \autoref{eq:F2_i} provide a smoother transition of the analysis trajectory from 87 95 one assimilation cycle to the next than that described by \autoref{eq:F1_i}. … … 92 100 \label{sec:ASM_div_dmp} 93 101 94 The velocity increments may be initialized by the iterative application of a divergence damping operator. 95 In iteration step $n$ new estimates of velocity increments $u^{n}_I$ and $v^{n}_I$ are updated by: 102 It is quite challenging for data assimilation systems to provide non-divergent velocity increments. 103 Applying divergent velocity increments will likely cause spurious vertical velocities in the model. This section describes a method to take velocity increments provided to \NEMO ($u^0_I$ and $v^0_I$) and adjust them by the iterative application of a divergence damping operator. The method is also described in \citet{dobricic.pinardi.ea_OS07}. 104 105 In iteration step $n$ (starting at $n=1$) new estimates of velocity increments $u^{n}_I$ and $v^{n}_I$ are updated by: 106 96 107 \begin{equation} 97 108 \label{eq:asm_dmp} … … 105 116 \right., 106 117 \end{equation} 107 where 118 119 where the divergence is defined as 120 108 121 \[ 109 122 % \label{eq:asm_div} … … 112 125 +\delta_j \left[ {e_{1v}\,e_{3v}\,v^{n-1}_I} \right]} \right). 113 126 \] 114 By the application of \autoref{eq:asm_dmp} and \autoref{eq:asm_dmp} the divergence is filtered in each iteration, 127 128 By the application of \autoref{eq:asm_dmp} the divergence is filtered in each iteration, 115 129 and the vorticity is left unchanged. 116 130 In the presence of coastal boundaries with zero velocity increments perpendicular to the coast … … 118 132 This type of the initialisation reduces the vertical velocity magnitude and 119 133 alleviates the problem of the excessive unphysical vertical mixing in the first steps of the model integration 120 \citep{ Talagrand_JAS72, Dobricic_al_OS07}.134 \citep{talagrand_JAS72, dobricic.pinardi.ea_OS07}. 121 135 Diffusion coefficients are defined as $A_D = \alpha e_{1t} e_{2t}$, where $\alpha = 0.2$. 122 136 The divergence damping is activated by assigning to \np{nn\_divdmp} in the \textit{nam\_asminc} namelist 123 137 a value greater than zero. 124 By choosing this value to be of the order of 100 the increments in 125 the vertical velocity will be significantly reduced. 138 This specifies the number of iterations of the divergence damping. Setting a value of the order of 100 will result in a significant reduction in the vertical velocity induced by the increments. 126 139 127 140 … … 131 144 \label{sec:ASM_details} 132 145 133 Here we show an example \ngn{nam asm} namelist and the header of an example assimilation increments file on146 Here we show an example \ngn{nam\_asminc} namelist and the header of an example assimilation increments file on 134 147 the ORCA2 grid. 135 148 136 %------------------------------------------nam asm-----------------------------------------------------149 %------------------------------------------nam_asminc----------------------------------------------------- 137 150 % 138 151 \nlst{nam_asminc} -
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_CONFIG.tex
r10442 r11422 18 18 \label{sec:CFG_intro} 19 19 20 The purpose of this part of the manual is to introduce the \NEMO reference configurations.20 The purpose of this part of the manual is to introduce the NEMO reference configurations. 21 21 These configurations are offered as means to explore various numerical and physical options, 22 22 thus allowing the user to verify that the code is performing in a manner consistent with that we are running. … … 24 24 The reference configurations also provide a sense for some of the options available in the code, 25 25 though by no means are all options exercised in the reference configurations. 26 Configuration is defined manually through the \textit{namcfg} namelist variables. 26 27 27 28 %------------------------------------------namcfg---------------------------------------------------- … … 33 34 % 1D model configuration 34 35 % ================================================================ 35 \section{C1D: 1D Water column model (\protect\key{c1d}) } 36 \section[C1D: 1D Water column model (\texttt{\textbf{key\_c1d}})] 37 {C1D: 1D Water column model (\protect\key{c1d})} 36 38 \label{sec:CFG_c1d} 37 39 38 BE careful: to be re-written according to suppression of jpizoom and jpjzoom !!!! 39 40 The 1D model option simulates a stand alone water column within the 3D \NEMO system. 40 The 1D model option simulates a stand alone water column within the 3D NEMO system. 41 41 It can be applied to the ocean alone or to the ocean-ice system and can include passive tracers or a biogeochemical model. 42 42 It is set up by defining the position of the 1D water column in the grid 43 (see \textit{ CONFIG/SHARED/namelist\_ref}).43 (see \textit{cfgs/SHARED/namelist\_ref}). 44 44 The 1D model is a very useful tool 45 45 \textit{(a)} to learn about the physics and numerical treatment of vertical mixing processes; … … 50 50 \textit{(d)} to produce extra diagnostics, without the large memory requirement of the full 3D model. 51 51 52 The methodology is based on the use of the zoom functionality over the smallest possible domain: 53 a 3x3 domain centered on the grid point of interest, with some extra routines. 54 There is no need to define a new mesh, bathymetry, initial state or forcing, 55 since the 1D model will use those of the configuration it is a zoom of. 56 The chosen grid point is set in \textit{\ngn{namcfg}} namelist by 57 setting the \np{jpizoom} and \np{jpjzoom} parameters to the indices of the location of the chosen grid point. 52 The methodology is based on the configuration of the smallest possible domain: 53 a 3x3 domain with 75 vertical levels. 58 54 59 55 The 1D model has some specifies. First, all the horizontal derivatives are assumed to be zero, 60 56 and second, the two components of the velocity are moved on a $T$-point. 61 Therefore, defining \key{c1d} changes five mainthings in the code behaviour:57 Therefore, defining \key{c1d} changes some things in the code behaviour: 62 58 \begin{description} 63 59 \item[(1)] 64 the lateral boundary condition routine (\rou{lbc\_lnk}) set the value of the central column of65 the 3x3 domain is imposed over the whole domain;66 \item[(3)]67 a call to \rou{lbc\_lnk} is systematically done when reading input data (\ie in \mdl{iom});68 \item[(3)]69 60 a simplified \rou{stp} routine is used (\rou{stp\_c1d}, see \mdl{step\_c1d} module) in which 70 61 both lateral tendancy terms and lateral physics are not called; 71 \item[( 4)]62 \item[(2)] 72 63 the vertical velocity is zero 73 64 (so far, no attempt at introducing a Ekman pumping velocity has been made); 74 \item[( 5)]65 \item[(3)] 75 66 a simplified treatment of the Coriolis term is performed as $U$- and $V$-points are the same 76 67 (see \mdl{dyncor\_c1d}). 77 68 \end{description} 78 All the relevant \textit{\_c1d} modules can be found in the NEMOGCM/NEMO/OPA\_SRC/C1D directory of79 the \NEMO distribution.69 All the relevant \textit{\_c1d} modules can be found in the src/OCE/C1D directory of 70 the NEMO distribution. 80 71 81 72 % to be added: a test case on the yearlong Ocean Weather Station (OWS) Papa dataset of Martin (1985) … … 88 79 89 80 The ORCA family is a series of global ocean configurations that are run together with 90 the LIM sea-ice model (ORCA-LIM) and possibly with PISCES biogeochemical model (ORCA-LIM-PISCES), 91 using various resolutions. 92 An appropriate namelist is available in \path{CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg} for ORCA2. 81 the SI3 model (ORCA-ICE) and possibly with PISCES biogeochemical model (ORCA-ICE-PISCES). 82 An appropriate namelist is available in \path{cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg} for ORCA2. 93 83 The domain of ORCA2 configuration is defined in \ifile{ORCA\_R2\_zps\_domcfg} file, 94 this file is available in tar file in the wiki of NEMO: \\ 95 https://forge.ipsl.jussieu.fr/nemo/wiki/Users/ReferenceConfigurations/ORCA2\_LIM3\_PISCES \\ 84 this file is available in tar file on the NEMO community zenodo platform: \\ 85 https://doi.org/10.5281/zenodo.2640723 86 96 87 In this namelist\_cfg the name of domain input file is set in \ngn{namcfg} block of namelist. 97 88 … … 99 90 \begin{figure}[!t] 100 91 \begin{center} 101 \includegraphics[width= 0.98\textwidth]{Fig_ORCA_NH_mesh}92 \includegraphics[width=\textwidth]{Fig_ORCA_NH_mesh} 102 93 \caption{ 103 94 \protect\label{fig:MISC_ORCA_msh} … … 106 97 The two "north pole" are the foci of a series of embedded ellipses (blue curves) which 107 98 are determined analytically and form the i-lines of the ORCA mesh (pseudo latitudes). 108 Then, following \citet{ Madec_Imbard_CD96}, the normal to the series of ellipses (red curves) is computed which99 Then, following \citet{madec.imbard_CD96}, the normal to the series of ellipses (red curves) is computed which 109 100 provides the j-lines of the mesh (pseudo longitudes). 110 101 } … … 119 110 \label{subsec:CFG_orca_grid} 120 111 121 The ORCA grid is a tripolar is based on the semi-analytical method of \citet{Madec_Imbard_CD96}.112 The ORCA grid is a tripolar grid based on the semi-analytical method of \citet{madec.imbard_CD96}. 122 113 It allows to construct a global orthogonal curvilinear ocean mesh which has no singularity point inside 123 114 the computational domain since two north mesh poles are introduced and placed on lands. … … 131 122 \begin{figure}[!tbp] 132 123 \begin{center} 133 \includegraphics[width= 1.0\textwidth]{Fig_ORCA_NH_msh05_e1_e2}134 \includegraphics[width= 0.80\textwidth]{Fig_ORCA_aniso}124 \includegraphics[width=\textwidth]{Fig_ORCA_NH_msh05_e1_e2} 125 \includegraphics[width=\textwidth]{Fig_ORCA_aniso} 135 126 \caption { 136 127 \protect\label{fig:MISC_ORCA_e1e2} … … 158 149 159 150 % ------------------------------------------------------------------------------------------------------------- 160 % ORCA- LIM(-PISCES) configurations151 % ORCA-ICE(-PISCES) configurations 161 152 % ------------------------------------------------------------------------------------------------------------- 162 153 \subsection{ORCA pre-defined resolution} … … 199 190 The ORCA\_R2 configuration has the following specificity: starting from a 2\deg~ORCA mesh, 200 191 local mesh refinements were applied to the Mediterranean, Red, Black and Caspian Seas, 201 so that the resolution is 1\deg \time 1\degthere.192 so that the resolution is 1\deg~ there. 202 193 A local transformation were also applied with in the Tropics in order to refine the meridional resolution up to 203 0.5\deg at the Equator.194 0.5\deg~ at the Equator. 204 195 205 196 The ORCA\_R1 configuration has only a local tropical transformation to refine the meridional resolution up to … … 211 202 For ORCA\_R1 and R025, setting the configuration key to 75 allows to use 75 vertical levels, otherwise 46 are used. 212 203 In the other ORCA configurations, 31 levels are used 213 (see \autoref{tab:orca_zgr} %\sfcomment{HERE I need to put new table for ORCA2 values} and \autoref{fig:zgr}).214 215 Only the ORCA\_R2 is provided with all its input files in the \NEMO distribution.216 It is very similar to that used as part of the climate model developed at IPSL for the 4th IPCC assessment of217 climate change (Marti et al., 2009).218 It is also the basis for the \NEMO contribution to the Coordinate Ocean-ice Reference Experiments (COREs)219 documented in \citet{Griffies_al_OM09}.204 (see \autoref{tab:orca_zgr}). %\sfcomment{HERE I need to put new table for ORCA2 values} and \autoref{fig:zgr}). 205 206 Only the ORCA\_R2 is provided with all its input files in the NEMO distribution. 207 %It is very similar to that used as part of the climate model developed at IPSL for the 4th IPCC assessment of 208 %climate change (Marti et al., 2009). 209 %It is also the basis for the \NEMO contribution to the Coordinate Ocean-ice Reference Experiments (COREs) 210 %documented in \citet{griffies.biastoch.ea_OM09}. 220 211 221 212 This version of ORCA\_R2 has 31 levels in the vertical, with the highest resolution (10m) in the upper 150m 222 213 (see \autoref{tab:orca_zgr} and \autoref{fig:zgr}). 223 214 The bottom topography and the coastlines are derived from the global atlas of Smith and Sandwell (1997). 224 The default forcing uses the boundary forcing from \citet{ Large_Yeager_Rep04} (see \autoref{subsec:SBC_blk_core}),215 The default forcing uses the boundary forcing from \citet{large.yeager_rpt04} (see \autoref{subsec:SBC_blk_core}), 225 216 which was developed for the purpose of running global coupled ocean-ice simulations without 226 217 an interactive atmosphere. 227 This \citet{ Large_Yeager_Rep04} dataset is available through218 This \citet{large.yeager_rpt04} dataset is available through 228 219 the \href{http://nomads.gfdl.noaa.gov/nomads/forms/mom4/CORE.html}{GFDL web site}. 229 The "normal year" of \citet{ Large_Yeager_Rep04} has been chosen of the \NEMO distribution since release v3.3.230 231 ORCA\_R2 pre-defined configuration can also be run with an AGRIF zoom over the Agulhas current area232 (\key{agrif} defined) and, by setting the appropriate variables, see \path{CONFIG/SHARED/namelist_ref}. 220 The "normal year" of \citet{large.yeager_rpt04} has been chosen of the NEMO distribution since release v3.3. 221 222 ORCA\_R2 pre-defined configuration can also be run with multiply online nested zooms (\ie with AGRIF, \key{agrif} defined). This is available as the AGRIF\_DEMO configuration that can be found in the \path{cfgs/AGRIF_DEMO/} directory. 223 233 224 A regional Arctic or peri-Antarctic configuration is extracted from an ORCA\_R2 or R05 configurations using 234 225 sponge layers at open boundaries. … … 237 228 % GYRE family: double gyre basin 238 229 % ------------------------------------------------------------------------------------------------------------- 239 \section{GYRE family: double gyre basin 230 \section{GYRE family: double gyre basin} 240 231 \label{sec:CFG_gyre} 241 232 242 The GYRE configuration \citep{ Levy_al_OM10} has been built to233 The GYRE configuration \citep{levy.klein.ea_OM10} has been built to 243 234 simulate the seasonal cycle of a double-gyre box model. 244 It consists in an idealized domain similar to that used in the studies of \citet{ Drijfhout_JPO94} and245 \citet{ Hazeleger_Drijfhout_JPO98, Hazeleger_Drijfhout_JPO99, Hazeleger_Drijfhout_JGR00, Hazeleger_Drijfhout_JPO00},235 It consists in an idealized domain similar to that used in the studies of \citet{drijfhout_JPO94} and 236 \citet{hazeleger.drijfhout_JPO98, hazeleger.drijfhout_JPO99, hazeleger.drijfhout_JGR00, hazeleger.drijfhout_JPO00}, 246 237 over which an analytical seasonal forcing is applied. 247 238 This allows to investigate the spontaneous generation of a large number of interacting, transient mesoscale eddies 248 239 and their contribution to the large scale circulation. 249 240 241 The GYRE configuration run together with the PISCES biogeochemical model (GYRE-PISCES). 250 242 The domain geometry is a closed rectangular basin on the $\beta$-plane centred at $\sim$ 30\deg{N} and 251 243 rotated by 45\deg, 3180~km long, 2120~km wide and 4~km deep (\autoref{fig:MISC_strait_hand}). … … 253 245 The configuration is meant to represent an idealized North Atlantic or North Pacific basin. 254 246 The circulation is forced by analytical profiles of wind and buoyancy fluxes. 255 The applied forcings vary seasonally in a sinusoidal manner between winter and summer extrema \citep{ Levy_al_OM10}.247 The applied forcings vary seasonally in a sinusoidal manner between winter and summer extrema \citep{levy.klein.ea_OM10}. 256 248 The wind stress is zonal and its curl changes sign at 22\deg{N} and 36\deg{N}. 257 249 It forces a subpolar gyre in the north, a subtropical gyre in the wider part of the domain and … … 266 258 The GYRE configuration is set like an analytical configuration. 267 259 Through \np{ln\_read\_cfg}\forcode{ = .false.} in \textit{namcfg} namelist defined in 268 the reference configuration \path{ CONFIG/GYRE/EXP00/namelist_cfg}260 the reference configuration \path{cfgs/GYRE_PISCES/EXPREF/namelist_cfg} 269 261 analytical definition of grid in GYRE is done in usrdef\_hrg, usrdef\_zgr routines. 270 262 Its horizontal resolution (and thus the size of the domain) is determined by 271 263 setting \np{nn\_GYRE} in \ngn{namusr\_def}: \\ 264 272 265 \np{jpiglo} $= 30 \times$ \np{nn\_GYRE} + 2 \\ 266 273 267 \np{jpjglo} $= 20 \times$ \np{nn\_GYRE} + 2 \\ 268 274 269 Obviously, the namelist parameters have to be adjusted to the chosen resolution, 275 see the Configurations pages on the NEMO web site ( Using NEMO\/Configurations).270 see the Configurations pages on the NEMO web site (NEMO Configurations). 276 271 In the vertical, GYRE uses the default 30 ocean levels (\jp{jpk}\forcode{ = 31}) (\autoref{fig:zgr}). 277 272 … … 281 276 even though the physical integrity of the solution can be compromised. 282 277 Benchmark is activate via \np{ln\_bench}\forcode{ = .true.} in \ngn{namusr\_def} in 283 namelist \path{ CONFIG/GYRE/EXP00/namelist_cfg}.278 namelist \path{cfgs/GYRE_PISCES/EXPREF/namelist_cfg}. 284 279 285 280 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 286 281 \begin{figure}[!t] 287 282 \begin{center} 288 \includegraphics[width= 1.0\textwidth]{Fig_GYRE}283 \includegraphics[width=\textwidth]{Fig_GYRE} 289 284 \caption{ 290 285 \protect\label{fig:GYRE} 291 286 Snapshot of relative vorticity at the surface of the model domain in GYRE R9, R27 and R54. 292 From \citet{ Levy_al_OM10}.287 From \citet{levy.klein.ea_OM10}. 293 288 } 294 289 \end{center} … … 304 299 The AMM, Atlantic Margins Model, is a regional model covering the Northwest European Shelf domain on 305 300 a regular lat-lon grid at approximately 12km horizontal resolution. 306 The appropriate \textit{\&namcfg} namelist is available in \textit{ CONFIG/AMM12/EXP00/namelist\_cfg}.301 The appropriate \textit{\&namcfg} namelist is available in \textit{cfgs/AMM12/EXPREF/namelist\_cfg}. 307 302 It is used to build the correct dimensions of the AMM domain. 308 303 309 304 This configuration tests several features of NEMO functionality specific to the shelf seas. 310 In particular, the AMM uses $S$-coordinates in the vertical rather than $z$-coordinates and 311 is forced with tidal lateral boundary conditions using a flather boundary condition from the BDY module. 312 The AMM configuration uses the GLS (\key{zdfgls}) turbulence scheme, 313 the VVL non-linear free surface(\key{vvl}) and time-splitting (\key{dynspg\_ts}). 305 In particular, the AMM uses $s$-coordinates in the vertical rather than $z$-coordinates and 306 is forced with tidal lateral boundary conditions using a Flather boundary condition from the BDY module. 307 Also specific to the AMM configuration is the use of the GLS turbulence scheme (\np{ln\_zdfgls} \forcode{= .true.}). 314 308 315 309 In addition to the tidal boundary condition the model may also take open boundary conditions from -
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_DIA.tex
r10509 r11422 10 10 \minitoc 11 11 12 \vfill 13 \begin{figure}[b] 14 \subsubsection*{Changes record} 15 \begin{tabular}{l||l|m{0.65\linewidth}} 16 Release & Author & Modifications \\ 17 {\em 4.0} & {\em Mirek Andrejczuk, Massimiliano Drudi} & {\em } \\ 18 {\em } & {\em Dorotea Iovino, Nicolas Martin} & {\em } \\ 19 {\em 3.6} & {\em Gurvan Madec, Sebastien Masson } & {\em } \\ 20 {\em 3.4} & {\em Gurvan Madec, Rachid Benshila, Andrew Coward } & {\em } \\ 21 {\em } & {\em Christian Ethe, Sebastien Masson } & {\em } \\ 22 \end{tabular} 23 \end{figure} 24 12 25 \newpage 13 26 … … 15 28 % Old Model Output 16 29 % ================================================================ 17 \section{ Old model output (default)}30 \section{Model output} 18 31 \label{sec:DIA_io_old} 19 32 … … 25 38 the same run performed in one step. 26 39 It should be noted that this requires that the restart file contains two consecutive time steps for 27 all the prognostic variables, and that it is saved in the same binary format as the one used by the computer that 28 is to read it (in particular, 32 bits binary IEEE format must not be used for this file). 40 all the prognostic variables. 29 41 30 42 The output listing and file(s) are predefined but should be checked and eventually adapted to the user's needs. … … 38 50 the writing work is distributed over the processors in massively parallel computing. 39 51 A complete description of the use of this I/O server is presented in the next section. 40 41 By default, \key{iomput} is not defined,42 NEMO produces NetCDF with the old IOIPSL library which has been kept for compatibility and its easy installation.43 However, the IOIPSL library is quite inefficient on parallel machines and, since version 3.2,44 many diagnostic options have been added presuming the use of \key{iomput}.45 The usefulness of the default IOIPSL-based option is expected to reduce with each new release.46 If \key{iomput} is not defined, output files and content are defined in the \mdl{diawri} module and47 contain mean (or instantaneous if \key{diainstant} is defined) values over a regular period of48 nn\_write time-steps (namelist parameter).49 52 50 53 %\gmcomment{ % start of gmcomment … … 91 94 in a very easy way. 92 95 All details of iomput functionalities are listed in the following subsections. 93 Examples of the XML files that control the outputs can be found in: \path{NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef.xml}, 94 \path{NEMOGCM/CONFIG/SHARED/field_def.xml} and \path{NEMOGCM/CONFIG/SHARED/domain_def.xml}. \\ 96 Examples of the XML files that control the outputs can be found in: 97 \path{cfgs/ORCA2_ICE_PISCES/EXPREF/iodef.xml}, 98 \path{cfgs/SHARED/field_def_nemo-oce.xml}, 99 \path{cfgs/SHARED/field_def_nemo-pisces.xml}, 100 \path{cfgs/SHARED/field_def_nemo-ice.xml} and \path{cfgs/SHARED/domain_def_nemo.xml}. \\ 95 101 96 102 The second functionality targets output performance when running in parallel (\key{mpp\_mpi}). … … 101 107 102 108 In version 3.6, the iom\_put interface depends on 103 an external code called \href{https://forge.ipsl.jussieu.fr/ioserver/browser/XIOS/branchs/xios- 1.0}{XIOS-1.0}109 an external code called \href{https://forge.ipsl.jussieu.fr/ioserver/browser/XIOS/branchs/xios-2.5}{XIOS-2.5} 104 110 (use of revision 618 or higher is required). 105 111 This new IO server can take advantage of the parallel I/O functionality of NetCDF4 to … … 168 174 \xmlline|<variable id="using_server" type="bool"></variable>| 169 175 170 The {\tt using\_server} setting determines whether or not the server will be used in \textit{attached mode}171 (as a library) [{\tt > false <}] or in \textit{detached mode}172 (as an external executable on N additional, dedicated cpus) [{\tt > true <}].176 The {\ttfamily using\_server} setting determines whether or not the server will be used in \textit{attached mode} 177 (as a library) [{\ttfamily> false <}] or in \textit{detached mode} 178 (as an external executable on N additional, dedicated cpus) [{\ttfamily > true <}]. 173 179 The \textit{attached mode} is simpler to use but much less efficient for massively parallel applications. 174 180 The type of each file can be either ''multiple\_file'' or ''one\_file''. … … 207 213 \subsubsection{Control of XIOS: the context in iodef.xml} 208 214 209 As well as the {\tt using\_server} flag, other controls on the use of XIOS are set in the XIOS context in iodef.xml.215 As well as the {\ttfamily using\_server} flag, other controls on the use of XIOS are set in the XIOS context in iodef.xml. 210 216 See the XML basics section below for more details on XML syntax and rules. 211 217 … … 257 263 See the installation guide on the \href{http://forge.ipsl.jussieu.fr/ioserver/wiki}{XIOS} wiki for help and guidance. 258 264 NEMO will need to link to the compiled XIOS library. 259 The \href{https://forge.ipsl.jussieu.fr/nemo/ wiki/Users/ModelInterfacing/InputsOutputs#Inputs-OutputsusingXIOS}260 { XIOS with NEMO} guide provides an example illustration of how this can be achieved.265 The \href{https://forge.ipsl.jussieu.fr/nemo/chrome/site/doc/NEMO/guide/html/install.html#extract-and-install-xios} 266 {Extract and install XIOS} guide provides an example illustration of how this can be achieved. 261 267 262 268 \subsubsection{Add your own outputs} … … 269 275 \begin{enumerate} 270 276 \item[1.] 271 in NEMO code, add a \forcode{CALL iom \_put( 'identifier', array )} where you want to output a 2D or 3D array.277 in NEMO code, add a \forcode{CALL iom_put( 'identifier', array )} where you want to output a 2D or 3D array. 272 278 \item[2.] 273 279 If necessary, add \forcode{USE iom ! I/O manager library} to the list of used modules in … … 442 448 \xmlline|<context src="./nemo_def.xml" />| 443 449 444 \noindent In NEMO, by default, the field and domain definition is done in 2 separate files: 445 \path{NEMOGCM/CONFIG/SHARED/field_def.xml} and \path{NEMOGCM/CONFIG/SHARED/domain_def.xml} that 450 \noindent In NEMO, by default, the field definition is done in 3 separate files ( 451 \path{cfgs/SHARED/field_def_nemo-oce.xml}, 452 \path{cfgs/SHARED/field_def_nemo-pisces.xml} and 453 \path{cfgs/SHARED/field_def_nemo-ice.xml} ) and the domain definition is done in another file ( \path{cfgs/SHARED/domain_def_nemo.xml} ) 454 that 446 455 are included in the main iodef.xml file through the following commands: 447 456 \begin{xmllines} 448 <field_definition src="./field_def.xml" /> 449 <domain_definition src="./domain_def.xml" /> 457 <context id="nemo" src="./context_nemo.xml"/> 450 458 \end{xmllines} 451 459 … … 508 516 509 517 Secondly, the group can be used to replace a list of elements. 510 Several examples of groups of fields are proposed at the end of the file \path{CONFIG/SHARED/field_def.xml}. 518 Several examples of groups of fields are proposed at the end of the XML field files ( 519 \path{cfgs/SHARED/field_def_nemo-oce.xml}, 520 \path{cfgs/SHARED/field_def_nemo-pisces.xml} and 521 \path{cfgs/SHARED/field_def_nemo-ice.xml} ) . 511 522 For example, a short list of the usual variables related to the U grid: 512 523 … … 514 525 <field_group id="groupU" > 515 526 <field field_ref="uoce" /> 516 <field field_ref="s uoce" />527 <field field_ref="ssu" /> 517 528 <field field_ref="utau" /> 518 529 </field_group> … … 538 549 the tag family domain. 539 550 It must therefore be done in the domain part of the XML file. 540 For example, in \path{ CONFIG/SHARED/domain_def.xml}, we provide the following example of a definition of551 For example, in \path{cfgs/SHARED/domain_def.xml}, we provide the following example of a definition of 541 552 a 5 by 5 box with the bottom left corner at point (10,10). 542 553 543 554 \begin{xmllines} 544 <domain _group id="grid_T">545 < domain id="myzoom" zoom_ibegin="10" zoom_jbegin="10" zoom_ni="5" zoom_nj="5" />555 <domain id="myzoomT" domain_ref="grid_T"> 556 <zoom_domain ibegin="10" jbegin="10" ni="5" nj="5" /> 546 557 \end{xmllines} 547 558 … … 551 562 \begin{xmllines} 552 563 <file id="myfile_vzoom" output_freq="1d" > 553 <field field_ref="toce" domain_ref="myzoom "/>564 <field field_ref="toce" domain_ref="myzoomT"/> 554 565 </file> 555 566 \end{xmllines} … … 576 587 \subsubsection{Define vertical zooms} 577 588 578 Vertical zooms are defined through the attributs zoom\_begin and zoom\_ endof the tag family axis.589 Vertical zooms are defined through the attributs zoom\_begin and zoom\_n of the tag family axis. 579 590 It must therefore be done in the axis part of the XML file. 580 For example, in \path{NEMOGCM/CONFIG/ORCA2_LIM/iodef_demo.xml}, we provide the following example: 581 582 \begin{xmllines} 583 <axis_group id="deptht" long_name="Vertical T levels" unit="m" positive="down" > 584 <axis id="deptht" /> 585 <axis id="deptht_myzoom" zoom_begin="1" zoom_end="10" /> 591 For example, in \path{cfgs/ORCA2_ICE_PISCES/EXPREF/iodef_demo.xml}, we provide the following example: 592 593 \begin{xmllines} 594 <axis_definition> 595 <axis id="deptht" long_name="Vertical T levels" unit="m" positive="down" /> 596 <axis id="deptht_zoom" azix_ref="deptht" > 597 <zoom_axis zoom_begin="1" zoom_n="10" /> 598 </axis> 586 599 \end{xmllines} 587 600 … … 765 778 \end{xmllines} 766 779 767 Note that, then the code is crashing, writting real4 variables forces a numerical conve ction from780 Note that, then the code is crashing, writting real4 variables forces a numerical conversion from 768 781 real8 to real4 which will create an internal error in NetCDF and will avoid the creation of the output files. 769 782 Forcing double precision outputs with prec="8" (for example in the field\_definition) will avoid this problem. … … 938 951 \hline 939 952 \end{tabularx} 940 \caption{Field tags ("\tt {field\_*}")}953 \caption{Field tags ("\ttfamily{field\_*}")} 941 954 \end{table} 942 955 … … 974 987 \hline 975 988 \end{tabularx} 976 \caption{File tags ("\tt {file\_*}")}989 \caption{File tags ("\ttfamily{file\_*}")} 977 990 \end{table} 978 991 … … 1007 1020 \hline 1008 1021 \end{tabularx} 1009 \caption{Axis tags ("\tt {axis\_*}")}1022 \caption{Axis tags ("\ttfamily{axis\_*}")} 1010 1023 \end{table} 1011 1024 … … 1040 1053 \hline 1041 1054 \end{tabularx} 1042 \caption{Domain tags ("\tt {domain\_*)}"}1055 \caption{Domain tags ("\ttfamily{domain\_*)}"} 1043 1056 \end{table} 1044 1057 … … 1073 1086 \hline 1074 1087 \end{tabularx} 1075 \caption{Grid tags ("\tt {grid\_*}")}1088 \caption{Grid tags ("\ttfamily{grid\_*}")} 1076 1089 \end{table} 1077 1090 … … 1114 1127 \hline 1115 1128 \end{tabularx} 1116 \caption{Reference attributes ("\tt {*\_ref}")}1129 \caption{Reference attributes ("\ttfamily{*\_ref}")} 1117 1130 \end{table} 1118 1131 … … 1150 1163 \hline 1151 1164 \end{tabularx} 1152 \caption{Domain attributes ("\tt {zoom\_*}")}1165 \caption{Domain attributes ("\ttfamily{zoom\_*}")} 1153 1166 \end{table} 1154 1167 … … 1318 1331 \subsection{CF metadata standard compliance} 1319 1332 1320 Output from the XIOS -1.0IO server is compliant with1333 Output from the XIOS IO server is compliant with 1321 1334 \href{http://cfconventions.org/Data/cf-conventions/cf-conventions-1.5/build/cf-conventions.html}{version 1.5} of 1322 1335 the CF metadata standard. … … 1332 1345 % NetCDF4 support 1333 1346 % ================================================================ 1334 \section{NetCDF4 support (\protect\key{netcdf4})} 1347 \section[NetCDF4 support (\texttt{\textbf{key\_netcdf4}})] 1348 {NetCDF4 support (\protect\key{netcdf4})} 1335 1349 \label{sec:DIA_nc4} 1336 1350 … … 1340 1354 Chunking and compression can lead to significant reductions in file sizes for a small runtime overhead. 1341 1355 For a fuller discussion on chunking and other performance issues the reader is referred to 1342 the NetCDF4 documentation found \href{http ://www.unidata.ucar.edu/software/netcdf/docs/netcdf.html#Chunking}{here}.1356 the NetCDF4 documentation found \href{https://www.unidata.ucar.edu/software/netcdf/docs/netcdf_perf_chunking.html}{here}. 1343 1357 1344 1358 The new features are only available when the code has been linked with a NetCDF4 library … … 1389 1403 \end{forlines} 1390 1404 1391 \noindent for a standard ORCA2\_LIM configuration gives chunksizes of {\small\tt 46x38x1} respectively in1392 the mono-processor case (\ie global domain of {\small\tt 182x149x31}).1405 \noindent for a standard ORCA2\_LIM configuration gives chunksizes of {\small\ttfamily 46x38x1} respectively in 1406 the mono-processor case (\ie global domain of {\small\ttfamily 182x149x31}). 1393 1407 An illustration of the potential space savings that NetCDF4 chunking and compression provides is given in 1394 1408 table \autoref{tab:NC4} which compares the results of two short runs of the ORCA2\_LIM reference configuration with … … 1450 1464 % Tracer/Dynamics Trends 1451 1465 % ------------------------------------------------------------------------------------------------------------- 1452 \section{Tracer/Dynamics trends (\protect\ngn{namtrd})} 1466 \section[Tracer/Dynamics trends (\texttt{namtrd})] 1467 {Tracer/Dynamics trends (\protect\ngn{namtrd})} 1453 1468 \label{sec:DIA_trd} 1454 1469 … … 1462 1477 (\ie at the end of each $dyn\cdots.F90$ and/or $tra\cdots.F90$ routines). 1463 1478 This capability is controlled by options offered in \ngn{namtrd} namelist. 1464 Note that the output are done with xIOS, and therefore the \key{IOM} is required.1479 Note that the output are done with XIOS, and therefore the \key{iomput} is required. 1465 1480 1466 1481 What is done depends on the \ngn{namtrd} logical set to \forcode{.true.}: … … 1488 1503 1489 1504 Note that the mixed layer tendency diagnostic can also be used on biogeochemical models via 1490 the \key{trdtrc} and \key{trdm ld\_trc} CPP keys.1505 the \key{trdtrc} and \key{trdmxl\_trc} CPP keys. 1491 1506 1492 1507 \textbf{Note that} in the current version (v3.6), many changes has been introduced but not fully tested. 1493 1508 In particular, options associated with \np{ln\_dyn\_mxl}, \np{ln\_vor\_trd}, and \np{ln\_tra\_mxl} are not working, 1494 and none of the options have been tested with variable volume (\ie \ key{vvl} defined).1509 and none of the options have been tested with variable volume (\ie \np{ln\_linssh}\forcode{ = .true.}). 1495 1510 1496 1511 % ------------------------------------------------------------------------------------------------------------- 1497 1512 % On-line Floats trajectories 1498 1513 % ------------------------------------------------------------------------------------------------------------- 1499 \section{FLO: On-Line Floats trajectories (\protect\key{floats})} 1514 \section[FLO: On-Line Floats trajectories (\texttt{\textbf{key\_floats}})] 1515 {FLO: On-Line Floats trajectories (\protect\key{floats})} 1500 1516 \label{sec:FLO} 1501 1517 %--------------------------------------------namflo------------------------------------------------------- … … 1506 1522 The on-line computation of floats advected either by the three dimensional velocity field or constraint to 1507 1523 remain at a given depth ($w = 0$ in the computation) have been introduced in the system during the CLIPPER project. 1508 Options are defined by \ngn{namflo} namelis variables.1509 The algorithm used is based either on the work of \cite{ Blanke_Raynaud_JPO97} (default option),1524 Options are defined by \ngn{namflo} namelist variables. 1525 The algorithm used is based either on the work of \cite{blanke.raynaud_JPO97} (default option), 1510 1526 or on a $4^th$ Runge-Hutta algorithm (\np{ln\_flork4}\forcode{ = .true.}). 1511 Note that the \cite{ Blanke_Raynaud_JPO97} algorithm have the advantage of providing trajectories which1527 Note that the \cite{blanke.raynaud_JPO97} algorithm have the advantage of providing trajectories which 1512 1528 are consistent with the numeric of the code, so that the trajectories never intercept the bathymetry. 1513 1529 … … 1519 1535 In case of Ariane convention, input filename is \np{init\_float\_ariane}. 1520 1536 Its format is: \\ 1521 {\scriptsize \texttt{I J K nisobfl itrash itrash}}1537 {\scriptsize \texttt{I J K nisobfl itrash}} 1522 1538 1523 1539 \noindent with: … … 1577 1593 In that case, output filename is trajec\_float. 1578 1594 1579 Another possiblity of writing format is Netcdf (\np{ln\_flo\_ascii}\forcode{ = .false.}). 1580 There are 2 possibilities: 1581 1582 - if (\key{iomput}) is used, outputs are selected in iodef.xml. 1595 Another possiblity of writing format is Netcdf (\np{ln\_flo\_ascii}\forcode{ = .false.}) with 1596 \key{iomput} and outputs selected in iodef.xml. 1583 1597 Here it is an example of specification to put in files description section: 1584 1598 … … 1597 1611 \end{xmllines} 1598 1612 1599 - if (\key{iomput}) is not used, a file called \ifile{trajec\_float} will be created by IOIPSL library.1600 1601 See also \href{http://stockage.univ-brest.fr/~grima/Ariane/}{here} the web site describing the off-line use of1602 this marvellous diagnostic tool.1603 1613 1604 1614 % ------------------------------------------------------------------------------------------------------------- 1605 1615 % Harmonic analysis of tidal constituents 1606 1616 % ------------------------------------------------------------------------------------------------------------- 1607 \section{Harmonic analysis of tidal constituents (\protect\key{diaharm}) } 1617 \section[Harmonic analysis of tidal constituents (\texttt{\textbf{key\_diaharm}})] 1618 {Harmonic analysis of tidal constituents (\protect\key{diaharm})} 1608 1619 \label{sec:DIA_diag_harm} 1609 1620 1610 %------------------------------------------nam dia_harm----------------------------------------------------1621 %------------------------------------------nam_diaharm---------------------------------------------------- 1611 1622 % 1612 1623 \nlst{nam_diaharm} … … 1616 1627 This on-line Harmonic analysis is actived with \key{diaharm}. 1617 1628 1618 Some parameters are available in namelist \ngn{nam dia\_harm}:1629 Some parameters are available in namelist \ngn{nam\_diaharm}: 1619 1630 1620 1631 - \np{nit000\_han} is the first time step used for harmonic analysis … … 1652 1663 % Sections transports 1653 1664 % ------------------------------------------------------------------------------------------------------------- 1654 \section{Transports across sections (\protect\key{diadct}) } 1665 \section[Transports across sections (\texttt{\textbf{key\_diadct}})] 1666 {Transports across sections (\protect\key{diadct})} 1655 1667 \label{sec:DIA_diag_dct} 1656 1668 … … 1664 1676 1665 1677 Each section is defined by the coordinates of its 2 extremities. 1666 The pathways between them are contructed using tools which can be found in \texttt{ NEMOGCM/TOOLS/SECTIONS\_DIADCT}1667 and are written in a binary file \texttt{section\_ijglobal.diadct \_ORCA2\_LIM} which is later read in by1678 The pathways between them are contructed using tools which can be found in \texttt{tools/SECTIONS\_DIADCT} 1679 and are written in a binary file \texttt{section\_ijglobal.diadct} which is later read in by 1668 1680 NEMO to compute on-line transports. 1669 1681 … … 1684 1696 \subsubsection{Creating a binary file containing the pathway of each section} 1685 1697 1686 In \texttt{ NEMOGCM/TOOLS/SECTIONS\_DIADCT/run},1698 In \texttt{tools/SECTIONS\_DIADCT/run}, 1687 1699 the file \textit{ {list\_sections.ascii\_global}} contains a list of all the sections that are to be computed 1688 1700 (this list of sections is based on MERSEA project metrics). … … 1733 1745 1734 1746 The script \texttt{job.ksh} computes the pathway for each section and creates a binary file 1735 \texttt{section\_ijglobal.diadct \_ORCA2\_LIM} which is read by NEMO. \\1747 \texttt{section\_ijglobal.diadct} which is read by NEMO. \\ 1736 1748 1737 1749 It is possible to use this tools for new configuations: \texttt{job.ksh} has to be updated with … … 1809 1821 The steric effect is therefore not explicitely represented. 1810 1822 This approximation does not represent a serious error with respect to the flow field calculated by the model 1811 \citep{ Greatbatch_JGR94}, but extra attention is required when investigating sea level,1823 \citep{greatbatch_JGR94}, but extra attention is required when investigating sea level, 1812 1824 as steric changes are an important contribution to local changes in sea level on seasonal and climatic time scales. 1813 1825 This is especially true for investigation into sea level rise due to global warming. 1814 1826 1815 1827 Fortunately, the steric contribution to the sea level consists of a spatially uniform component that 1816 can be diagnosed by considering the mass budget of the world ocean \citep{ Greatbatch_JGR94}.1828 can be diagnosed by considering the mass budget of the world ocean \citep{greatbatch_JGR94}. 1817 1829 In order to better understand how global mean sea level evolves and thus how the steric sea level can be diagnosed, 1818 1830 we compare, in the following, the non-Boussinesq and Boussinesq cases. … … 1888 1900 the ocean surface, not by changes in mean mass of the ocean: the steric effect is missing in a Boussinesq fluid. 1889 1901 1890 Nevertheless, following \citep{ Greatbatch_JGR94}, the steric effect on the volume can be diagnosed by1902 Nevertheless, following \citep{greatbatch_JGR94}, the steric effect on the volume can be diagnosed by 1891 1903 considering the mass budget of the ocean. 1892 1904 The apparent changes in $\mathcal{M}$, mass of the ocean, which are not induced by surface mass flux 1893 1905 must be compensated by a spatially uniform change in the mean sea level due to expansion/contraction of the ocean 1894 \citep{ Greatbatch_JGR94}.1906 \citep{greatbatch_JGR94}. 1895 1907 In others words, the Boussinesq mass, $\mathcal{M}_o$, can be related to $\mathcal{M}$, 1896 1908 the total mass of the ocean seen by the Boussinesq model, via the steric contribution to the sea level, … … 1924 1936 This value is a sensible choice for the reference density used in a Boussinesq ocean climate model since, 1925 1937 with the exception of only a small percentage of the ocean, density in the World Ocean varies by no more than 1926 2$\%$ from this value (\cite{ Gill1982}, page 47).1938 2$\%$ from this value (\cite{gill_bk82}, page 47). 1927 1939 1928 1940 Second, we have assumed here that the total ocean surface, $\mathcal{A}$, … … 1931 1943 1932 1944 Third, the discretisation of \autoref{eq:steric_Bq} depends on the type of free surface which is considered. 1933 In the non linear free surface case, \ie \ key{vvl} defined, it is given by1945 In the non linear free surface case, \ie \np{ln\_linssh}\forcode{ = .true.}, it is given by 1934 1946 1935 1947 \[ … … 1954 1966 so that there are no associated ocean currents. 1955 1967 Hence, the dynamically relevant sea level is the effective sea level, 1956 \ie the sea level as if sea ice (and snow) were converted to liquid seawater \citep{ Campin_al_OM08}.1968 \ie the sea level as if sea ice (and snow) were converted to liquid seawater \citep{campin.marshall.ea_OM08}. 1957 1969 However, in the current version of \NEMO the sea-ice is levitating above the ocean without mass exchanges between 1958 1970 ice and ocean. … … 1970 1982 where $S_o$ and $p_o$ are the initial salinity and pressure, respectively. 1971 1983 1972 Both steric and thermosteric sea level are computed in \mdl{diaar5} which needs the \key{diaar5} defined to 1973 be called. 1984 Both steric and thermosteric sea level are computed in \mdl{diaar5}. 1974 1985 1975 1986 % ------------------------------------------------------------------------------------------------------------- 1976 1987 % Other Diagnostics 1977 1988 % ------------------------------------------------------------------------------------------------------------- 1978 \section{Other diagnostics (\protect\key{diahth}, \protect\key{diaar5})} 1989 \section[Other diagnostics] 1990 {Other diagnostics} 1979 1991 \label{sec:DIA_diag_others} 1980 1992 … … 1982 1994 The available ready-to-add diagnostics modules can be found in directory DIA. 1983 1995 1984 \subsection{Depth of various quantities (\protect\mdl{diahth})} 1996 \subsection[Depth of various quantities (\textit{diahth.F90})] 1997 {Depth of various quantities (\protect\mdl{diahth})} 1985 1998 1986 1999 Among the available diagnostics the following ones are obtained when defining the \key{diahth} CPP key: 1987 2000 1988 - the mixed layer depth (based on a density criterion \citep{de _Boyer_Montegut_al_JGR04}) (\mdl{diahth})2001 - the mixed layer depth (based on a density criterion \citep{de-boyer-montegut.madec.ea_JGR04}) (\mdl{diahth}) 1989 2002 1990 2003 - the turbocline depth (based on a turbulent mixing coefficient criterion) (\mdl{diahth}) … … 1994 2007 - the depth of the thermocline (maximum of the vertical temperature gradient) (\mdl{diahth}) 1995 2008 1996 % -----------------------------------------------------------1997 % Poleward heat and salt transports1998 % -----------------------------------------------------------1999 2000 \subsection{Poleward heat and salt transports (\protect\mdl{diaptr})}2001 2002 %------------------------------------------namptr-----------------------------------------2003 2004 \nlst{namptr}2005 %-----------------------------------------------------------------------------------------2006 2007 The poleward heat and salt transports, their advective and diffusive component,2008 and the meriodional stream function can be computed on-line in \mdl{diaptr} \np{ln\_diaptr} to true2009 (see the \textit{\ngn{namptr} } namelist below).2010 When \np{ln\_subbas}\forcode{ = .true.}, transports and stream function are computed for the Atlantic, Indian,2011 Pacific and Indo-Pacific Oceans (defined north of 30\deg{S}) as well as for the World Ocean.2012 The sub-basin decomposition requires an input file (\ifile{subbasins}) which contains three 2D mask arrays,2013 the Indo-Pacific mask been deduced from the sum of the Indian and Pacific mask (\autoref{fig:mask_subasins}).2014 2009 2015 2010 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 2016 2011 \begin{figure}[!t] 2017 2012 \begin{center} 2018 \includegraphics[width= 1.0\textwidth]{Fig_mask_subasins}2013 \includegraphics[width=\textwidth]{Fig_mask_subasins} 2019 2014 \caption{ 2020 2015 \protect\label{fig:mask_subasins} … … 2032 2027 % CMIP specific diagnostics 2033 2028 % ----------------------------------------------------------- 2034 \subsection{CMIP specific diagnostics (\protect\mdl{diaar5})} 2035 2036 A series of diagnostics has been added in the \mdl{diaar5}. 2037 They corresponds to outputs that are required for AR5 simulations (CMIP5) 2029 \subsection[CMIP specific diagnostics (\textit{diaar5.F90}, \textit{diaptr.F90})] 2030 {CMIP specific diagnostics (\protect\mdl{diaar5})} 2031 2032 A series of diagnostics has been added in the \mdl{diaar5} and \mdl{diaptr}. 2033 In \mdl{diaar5} they correspond to outputs that are required for AR5 simulations (CMIP5) 2038 2034 (see also \autoref{sec:DIA_steric} for one of them). 2039 Activating those outputs requires to define the \key{diaar5} CPP key. 2035 The module \mdl{diaar5} is active when one of the following outputs is required : global total volume (voltot), global mean ssh (sshtot), global total mass (masstot), global mean temperature (temptot), global mean ssh steric (sshsteric), global mean ssh thermosteric (sshthster), global mean salinity (saltot), sea water pressure at sea floor (botpres), dynamic sea surface height (sshdyn). 2036 2037 In \mdl{diaptr} when \np{ln\_diaptr}\forcode{ = .true.} 2038 (see the \textit{\ngn{namptr} } namelist below) can be computed on-line the poleward heat and salt transports, their advective and diffusive component, and the meriodional stream function . 2039 When \np{ln\_subbas}\forcode{ = .true.}, transports and stream function are computed for the Atlantic, Indian, 2040 Pacific and Indo-Pacific Oceans (defined north of 30\deg{S}) as well as for the World Ocean. 2041 The sub-basin decomposition requires an input file (\ifile{subbasins}) which contains three 2D mask arrays, 2042 the Indo-Pacific mask been deduced from the sum of the Indian and Pacific mask (\autoref{fig:mask_subasins}). 2043 2044 %------------------------------------------namptr----------------------------------------- 2045 2046 \nlst{namptr} 2047 %----------------------------------------------------------------------------------------- 2040 2048 2041 2049 % ----------------------------------------------------------- -
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_DIU.tex
r10442 r11422 33 33 (\ie from the temperature of the top few model levels) or from some other source. 34 34 It must be noted that both the cool skin and warm layer models produce estimates of the change in temperature 35 ($\Delta T_{\ rm{cs}}$ and $\Delta T_{\rm{wl}}$) and35 ($\Delta T_{\mathrm{cs}}$ and $\Delta T_{\mathrm{wl}}$) and 36 36 both must be added to a foundation SST to obtain the true skin temperature. 37 37 … … 60 60 %=============================================================== 61 61 62 The warm layer is calculated using the model of \citet{ Takaya_al_JGR10} (TAKAYA10 model hereafter).62 The warm layer is calculated using the model of \citet{takaya.bidlot.ea_JGR10} (TAKAYA10 model hereafter). 63 63 This is a simple flux based model that is defined by the equations 64 64 \begin{align} 65 \frac{\partial{\Delta T_{\ rm{wl}}}}{\partial{t}}&=&\frac{Q(\nu+1)}{D_T\rho_w c_p65 \frac{\partial{\Delta T_{\mathrm{wl}}}}{\partial{t}}&=&\frac{Q(\nu+1)}{D_T\rho_w c_p 66 66 \nu}-\frac{(\nu+1)ku^*_{w}f(L_a)\Delta T}{D_T\Phi\!\left(\frac{D_T}{L}\right)} \mbox{,} 67 67 \label{eq:ecmwf1} \\ 68 68 L&=&\frac{\rho_w c_p u^{*^3}_{w}}{\kappa g \alpha_w Q }\mbox{,}\label{eq:ecmwf2} 69 69 \end{align} 70 where $\Delta T_{\ rm{wl}}$ is the temperature difference between the top of the warm layer and the depth $D_T=3$\,m at which there is assumed to be no diurnal signal.70 where $\Delta T_{\mathrm{wl}}$ is the temperature difference between the top of the warm layer and the depth $D_T=3$\,m at which there is assumed to be no diurnal signal. 71 71 In equation (\autoref{eq:ecmwf1}) $\alpha_w=2\times10^{-4}$ is the thermal expansion coefficient of water, 72 72 $\kappa=0.4$ is von K\'{a}rm\'{a}n's constant, $c_p$ is the heat capacity at constant pressure of sea water, 73 73 $\rho_w$ is the water density, and $L$ is the Monin-Obukhov length. 74 74 The tunable variable $\nu$ is a shape parameter that defines the expected subskin temperature profile via 75 $T(z) = T(0) - \left( \frac{z}{D_T} \right)^\nu \Delta T_{\ rm{wl}}$,75 $T(z) = T(0) - \left( \frac{z}{D_T} \right)^\nu \Delta T_{\mathrm{wl}}$, 76 76 where $T$ is the absolute temperature and $z\le D_T$ is the depth below the top of the warm layer. 77 77 The influence of wind on TAKAYA10 comes through the magnitude of the friction velocity of the water $u^*_{w}$, … … 82 82 the diurnal layer, \ie 83 83 \[ 84 Q = Q_{\ rm{sol}} + Q_{\rm{lw}} + Q_{\rm{h}}\mbox{,}84 Q = Q_{\mathrm{sol}} + Q_{\mathrm{lw}} + Q_{\mathrm{h}}\mbox{,} 85 85 % \label{eq:e_flux_eqn} 86 86 \] 87 where $Q_{\ rm{h}}$ is the sensible and latent heat flux, $Q_{\rm{lw}}$ is the long wave flux,88 and $Q_{\ rm{sol}}$ is the solar flux absorbed within the diurnal warm layer.89 For $Q_{\ rm{sol}}$ the 9 term representation of \citet{Gentemann_al_JGR09} is used.87 where $Q_{\mathrm{h}}$ is the sensible and latent heat flux, $Q_{\mathrm{lw}}$ is the long wave flux, 88 and $Q_{\mathrm{sol}}$ is the solar flux absorbed within the diurnal warm layer. 89 For $Q_{\mathrm{sol}}$ the 9 term representation of \citet{gentemann.minnett.ea_JGR09} is used. 90 90 In equation \autoref{eq:ecmwf1} the function $f(L_a)=\max(1,L_a^{\frac{2}{3}})$, 91 91 where $L_a=0.3$\footnote{ … … 118 118 %=============================================================== 119 119 120 The cool skin is modelled using the framework of \citet{ Saunders_JAS82} who used a formulation of the near surface temperature difference based upon the heat flux and the friction velocity $u^*_{w}$.121 As the cool skin is so thin (~1\,mm) we ignore the solar flux component to the heat flux and the Saunders equation for the cool skin temperature difference $\Delta T_{\ rm{cs}}$ becomes120 The cool skin is modelled using the framework of \citet{saunders_JAS67} who used a formulation of the near surface temperature difference based upon the heat flux and the friction velocity $u^*_{w}$. 121 As the cool skin is so thin (~1\,mm) we ignore the solar flux component to the heat flux and the Saunders equation for the cool skin temperature difference $\Delta T_{\mathrm{cs}}$ becomes 122 122 \[ 123 123 % \label{eq:sunders_eqn} 124 \Delta T_{\ rm{cs}}=\frac{Q_{\rm{ns}}\delta}{k_t} \mbox{,}124 \Delta T_{\mathrm{cs}}=\frac{Q_{\mathrm{ns}}\delta}{k_t} \mbox{,} 125 125 \] 126 where $Q_{\ rm{ns}}$ is the, usually negative, non-solar heat flux into the ocean and126 where $Q_{\mathrm{ns}}$ is the, usually negative, non-solar heat flux into the ocean and 127 127 $k_t$ is the thermal conductivity of sea water. 128 128 $\delta$ is the thickness of the skin layer and is given by … … 132 132 \end{equation} 133 133 where $\mu$ is the kinematic viscosity of sea water and $\lambda$ is a constant of proportionality which 134 \citet{ Saunders_JAS82} suggested varied between 5 and 10.134 \citet{saunders_JAS67} suggested varied between 5 and 10. 135 135 136 The value of $\lambda$ used in equation (\autoref{eq:sunders_thick_eqn}) is that of \citet{ Artale_al_JGR02},137 which is shown in \citet{ Tu_Tsuang_GRL05} to outperform a number of other parametrisations at136 The value of $\lambda$ used in equation (\autoref{eq:sunders_thick_eqn}) is that of \citet{artale.iudicone.ea_JGR02}, 137 which is shown in \citet{tu.tsuang_GRL05} to outperform a number of other parametrisations at 138 138 both low and high wind speeds. 139 139 Specifically, -
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_DOM.tex
r10502 r11422 18 18 % - domclo: closed sea and lakes.... management of closea sea area : specific to global configuration, both forced and coupled 19 19 20 \vfill 21 \begin{figure}[b] 22 \subsubsection*{Changes record} 23 \begin{tabular}{m{0.08\linewidth}||m{0.32\linewidth}|m{0.6\linewidth}} 24 Release & Author(s) & Modifications \\ 25 \hline 26 {\em 4.0} & {\em Simon M{\"u}ller \& Andrew Coward} & {\em Compatibility changes for v4.0. Major simplication has moved many of the options to external domain configuration tools. For now this information has been retained in \autoref{apdx:DOMAINcfg} } \\ 27 {\em 3.x} & {\em Sebastien Masson, Gurvan Madec \& Rashid Benshila } & {\em } \\ 28 \end{tabular} 29 \end{figure} 30 20 31 \newpage 21 32 22 33 Having defined the continuous equations in \autoref{chap:PE} and chosen a time discretization \autoref{chap:STP}, 23 we need to choose a discretization on a grid, and numerical algorithms.34 we need to choose a grid for spatial discretization and related numerical algorithms. 24 35 In the present chapter, we provide a general description of the staggered grid used in \NEMO, 25 and other information relevant to the main directory routines as well as the DOM (DOMain) directory.36 and other relevant information about the DOM (DOMain) source-code modules . 26 37 27 38 % ================================================================ … … 40 51 \begin{figure}[!tb] 41 52 \begin{center} 42 \includegraphics[ ]{Fig_cell}53 \includegraphics[width=\textwidth]{Fig_cell} 43 54 \caption{ 44 55 \protect\label{fig:cell} … … 55 66 The numerical techniques used to solve the Primitive Equations in this model are based on the traditional, 56 67 centred second-order finite difference approximation. 57 Special attention has been given to the homogeneity of the solution in the three spa cedirections.68 Special attention has been given to the homogeneity of the solution in the three spatial directions. 58 69 The arrangement of variables is the same in all directions. 59 70 It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector points $(u, v, w)$ defined in 60 71 the centre of each face of the cells (\autoref{fig:cell}). 61 72 This is the generalisation to three dimensions of the well-known ``C'' grid in Arakawa's classification 62 \citep{ Mesinger_Arakawa_Bk76}.73 \citep{mesinger.arakawa_bk76}. 63 74 The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge and 64 75 the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points. … … 71 82 Each scale factor is defined as the local analytical value provided by \autoref{eq:scale_factors}. 72 83 As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and 73 $\pd[]{z}$ are evaluated i na uniform mesh with a grid size of unity.84 $\pd[]{z}$ are evaluated is a uniform mesh with a grid size of unity. 74 85 Discrete partial derivatives are formulated by the traditional, centred second order finite difference approximation 75 86 while the scale factors are chosen equal to their local analytical value. … … 79 90 the continuous properties (see \autoref{apdx:C}). 80 91 A similar, related remark can be made about the domain size: 81 when needed, an area, volume, or the total ocean depth must be evaluated as the sum of the relevant scale factors92 when needed, an area, volume, or the total ocean depth must be evaluated as the product or sum of the relevant scale factors 82 93 (see \autoref{eq:DOM_bar} in the next section). 83 94 … … 87 98 \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|} 88 99 \hline 89 T& $i $ & $j $ & $k $ \\100 t & $i $ & $j $ & $k $ \\ 90 101 \hline 91 102 u & $i + 1/2$ & $j $ & $k $ \\ … … 107 118 \protect\label{tab:cell} 108 119 Location of grid-points as a function of integer or integer and a half value of the column, line or level. 109 This indexing is only used for the writing of the semi -discrete equation .110 In the code, the indexing uses integer values only and has a reverse direction in the vertical120 This indexing is only used for the writing of the semi -discrete equations. 121 In the code, the indexing uses integer values only and is positive downwards in the vertical with $k=1$ at the surface. 111 122 (see \autoref{subsec:DOM_Num_Index}) 112 123 } 113 124 \end{center} 114 125 \end{table} 126 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 127 128 Note that the definition of the scale factors 129 (\ie as the analytical first derivative of the transformation that 130 results in $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) 131 is specific to the \NEMO model \citep{marti.madec.ea_JGR92}. 132 As an example, a scale factor in the $i$ direction is defined locally at a $t$-point, 133 whereas many other models on a C grid choose to define such a scale factor as 134 the distance between the $u$-points on each side of the $t$-point. 135 Relying on an analytical transformation has two advantages: 136 firstly, there is no ambiguity in the scale factors appearing in the discrete equations, 137 since they are first introduced in the continuous equations; 138 secondly, analytical transformations encourage good practice by the definition of smoothly varying grids 139 (rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{treguier.dukowicz.ea_JGR96}. 140 An example of the effect of such a choice is shown in \autoref{fig:zgr_e3}. 141 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 142 \begin{figure}[!t] 143 \begin{center} 144 \includegraphics[width=\textwidth]{Fig_zgr_e3} 145 \caption{ 146 \protect\label{fig:zgr_e3} 147 Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, 148 and (b) analytically derived grid-point position and scale factors. 149 For both grids here, the same $w$-point depth has been chosen but 150 in (a) the $t$-points are set half way between $w$-points while 151 in (b) they are defined from an analytical function: 152 $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$. 153 Note the resulting difference between the value of the grid-size $\Delta_k$ and 154 those of the scale factor $e_k$. 155 } 156 \end{center} 157 \end{figure} 115 158 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 116 159 … … 132 175 Following \autoref{eq:PE_grad} and \autoref{eq:PE_lap}, the gradient of a variable $q$ defined at 133 176 a $t$-point has its three components defined at $u$-, $v$- and $w$-points while 134 its Laplacian is defined at $t$-point.177 its Laplacian is defined at the $t$-point. 135 178 These operators have the following discrete forms in the curvilinear $s$-coordinates system: 136 179 \[ … … 171 214 \end{equation} 172 215 173 The vertical average over the whole water column denoted by an overbar becomes for a quantity $q$ which174 is a masked field (i.e. equal to zero inside solid area):216 The vertical average over the whole water column is denoted by an overbar and is for 217 a masked field $q$ (\ie a quantity that is equal to zero inside solid areas): 175 218 \begin{equation} 176 219 \label{eq:DOM_bar} … … 178 221 \end{equation} 179 222 where $H_q$ is the ocean depth, which is the masked sum of the vertical scale factors at $q$ points, 180 $k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $ k^o$ refers to a summation over223 $k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $\sum \limits_k$ refers to a summation over 181 224 all grid points of the same type in the direction indicated by the subscript (here $k$). 182 225 … … 193 236 vector points $(u,v,w)$. 194 237 195 Let $a$ and $b$ be two fields defined on the mesh, with value zero inside continental area.196 Using integration by parts it can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$)238 Let $a$ and $b$ be two fields defined on the mesh, with a value of zero inside continental areas. 239 It can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) 197 240 are skew-symmetric linear operators, and further that the averaging operators $\overline{\cdots}^{\, i}$, 198 241 $\overline{\cdots}^{\, j}$ and $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie … … 218 261 \begin{figure}[!tb] 219 262 \begin{center} 220 \includegraphics[ ]{Fig_index_hor}263 \includegraphics[width=\textwidth]{Fig_index_hor} 221 264 \caption{ 222 265 \protect\label{fig:index_hor} … … 228 271 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 229 272 230 The array representation used in the \fortran code requires an integer indexing while231 the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with the use of232 integer values for $t$-points and both integer and integer and a half values for all the other points.233 Therefore a specific integer indexing must bedefined for points other than $t$-points273 The array representation used in the \fortran code requires an integer indexing. 274 However, the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with the use of 275 integer values for $t$-points only while all the other points involve integer and a half values. 276 Therefore, a specific integer indexing has been defined for points other than $t$-points 234 277 (\ie velocity and vorticity grid-points). 235 Furthermore, the direction of the vertical indexing has been changed so that the surface level isat $k = 1$.278 Furthermore, the direction of the vertical indexing has been reversed and the surface level set at $k = 1$. 236 279 237 280 % ----------------------------------- … … 253 296 \label{subsec:DOM_Num_Index_vertical} 254 297 255 In the vertical, the chosen indexing requires special attention since the $k$-axis is re-orientated downwardin256 the \fortran code compared to the indexingused in the semi -discrete equations and298 In the vertical, the chosen indexing requires special attention since the direction of the $k$-axis in 299 the \fortran code is the reverse of that used in the semi -discrete equations and 257 300 given in \autoref{subsec:DOM_cell}. 258 The sea surface corresponds to the $w$-level $k = 1$ which is the same index as$t$-level just below301 The sea surface corresponds to the $w$-level $k = 1$, which is the same index as the $t$-level just below 259 302 (\autoref{fig:index_vert}). 260 The last $w$-level ($k = jpk$) either corresponds to the ocean floor or is inside the bathymetry while 261 the last $t$-level is always inside the bathymetry (\autoref{fig:index_vert}). 262 Note that for an increasing $k$ index, a $w$-point and the $t$-point just below have the same $k$ index, 263 in opposition to what is done in the horizontal plane where 264 it is the $t$-point and the nearest velocity points in the direction of the horizontal axis that 265 have the same $i$ or $j$ index 303 The last $w$-level ($k = jpk$) either corresponds to or is below the ocean floor while 304 the last $t$-level is always outside the ocean domain (\autoref{fig:index_vert}). 305 Note that a $w$-point and the directly underlaying $t$-point have a common $k$ index (\ie $t$-points and their 306 nearest $w$-point neighbour in negative index direction), in contrast to the indexing on the horizontal plane where 307 the $t$-point has the same index as the nearest velocity points in the positive direction of the respective horizontal axis index 266 308 (compare the dashed area in \autoref{fig:index_hor} and \autoref{fig:index_vert}). 267 309 Since the scale factors are chosen to be strictly positive, 268 a \textit{minus sign} appears in the \fortran code \textit{beforeall the vertical derivatives} of269 the discrete equations given in this documentation.310 a \textit{minus sign} is included in the \fortran implementations of \textit{all the vertical derivatives} of 311 the discrete equations given in this manual in order to accommodate the opposing vertical index directions in implementation and documentation. 270 312 271 313 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 272 314 \begin{figure}[!pt] 273 315 \begin{center} 274 \includegraphics[ ]{Fig_index_vert}316 \includegraphics[width=\textwidth]{Fig_index_vert} 275 317 \caption{ 276 318 \protect\label{fig:index_vert} 277 319 Vertical integer indexing used in the \fortran code. 278 Note that the $k$-axis is orient ated downward.279 The dashed area indicates the cell in which variables contained in arrays have the same$k$-index.320 Note that the $k$-axis is oriented downward. 321 The dashed area indicates the cell in which variables contained in arrays have a common $k$-index. 280 322 } 281 323 \end{center} … … 283 325 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 284 326 327 % ------------------------------------------------------------------------------------------------------------- 328 % Domain configuration 329 % ------------------------------------------------------------------------------------------------------------- 330 \section{Spatial domain configuration} 331 \label{subsec:DOM_config} 332 333 \nlst{namcfg} 334 335 Two typical methods are available to specify the spatial domain 336 configuration; they can be selected using parameter \np{ln\_read\_cfg} 337 parameter in namelist \ngn{namcfg}. 338 339 If \np{ln\_read\_cfg} is set to \forcode{.true.}, the domain-specific parameters 340 and fields are read from a netCDF input file, whose name (without its .nc 341 suffix) can be specified as the value of the \np{cn\_domcfg} parameter in 342 namelist \ngn{namcfg}. 343 344 If \np{ln\_read\_cfg} is set to \forcode{.false.}, the domain-specific 345 parameters and fields can be provided (\eg analytically computed) by subroutines 346 \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}. These subroutines can be supplied in 347 the \path{MY_SRC} directory of the configuration, and default versions that 348 configure the spatial domain for the GYRE reference configuration are present in 349 the \path{src/OCE/USR} directory. 350 351 In version 4.0 there are no longer any options for reading complex bathmetries and 352 performing a vertical discretization at run-time. Whilst it is occasionally convenient 353 to have a common bathymetry file and, for example, to run similar models with and 354 without partial bottom boxes and/or sigma-coordinates, supporting such choices leads to 355 overly complex code. Worse still is the difficulty of ensuring the model configurations 356 intended to be identical are indeed so when the model domain itself can be altered by runtime 357 selections. The code previously used to perform vertical discretization has be incorporated 358 into an external tool (\path{tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMAINcfg}. 359 360 The next subsections summarise the parameter and fields related to the 361 configuration of the whole model domain. These represent the minimum information 362 that must be provided either via the \np{cn\_domcfg} file or set by code 363 inserted into user-supplied versions of the \mdl{usrdef\_*} subroutines. The 364 requirements are presented in three sections: the domain size 365 (\autoref{subsec:DOM_size}), the horizontal mesh 366 (\autoref{subsec:DOM_hgr}), and the vertical grid 367 (\autoref{subsec:DOM_zgr}). 368 285 369 % ----------------------------------- 286 370 % Domain Size 287 371 % ----------------------------------- 288 \subs ubsection{Domain size}372 \subsection{Domain size} 289 373 \label{subsec:DOM_size} 290 374 291 The total size of the computational domain is set by the parameters \np{jpiglo}, 292 \np{jpjglo} and \np{jpkglo} in the $i$, $j$ and $k$ directions respectively. 293 Parameters $jpi$ and $jpj$ refer to the size of each processor subdomain when 294 the code is run in parallel using domain decomposition (\key{mpp\_mpi} defined, 295 see \autoref{sec:LBC_mpp}). 296 297 % ================================================================ 298 % Domain: List of fields needed 299 % ================================================================ 300 \section{Needed fields} 301 \label{sec:DOM_fields} 302 The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that 303 gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 304 The grid-points are located at integer or integer and a half values of as indicated in \autoref{tab:cell}. 305 The associated scale factors are defined using the analytical first derivative of the transformation 306 \autoref{eq:scale_factors}. 307 Necessary fields for configuration definition are: 308 309 \begin{itemize} 310 \item 311 Geographic position: 312 longitude with \texttt{glamt}, \texttt{glamu}, \texttt{glamv}, \texttt{glamf} and 313 latitude with \texttt{gphit}, \texttt{gphiu}, \texttt{gphiv}, \texttt{gphif} 314 (all respectively at T, U, V and F point) 315 \item 316 Coriolis parameter (if domain not on the sphere): \texttt{ff\_f} and \texttt{ff\_t} 317 (at T and F point) 318 \item 319 Scale factors: 320 \texttt{e1t}, \texttt{e1u}, \texttt{e1v} and \texttt{e1f} (on i direction), 321 \texttt{e2t}, \texttt{e2u}, \texttt{e2v} and \texttt{e2f} (on j direction) and 322 \texttt{ie1e2u\_v}, \texttt{e1e2u}, \texttt{e1e2v}. \\ 323 \texttt{e1e2u}, \texttt{e1e2v} are u and v surfaces (if gridsize reduction in some straits), 324 \texttt{ie1e2u\_v} is to flag set u and v surfaces are neither read nor computed. 325 \end{itemize} 326 327 These fields can be read in an domain input file which name is setted in \np{cn\_domcfg} parameter specified in 328 \ngn{namcfg}. 329 330 \nlst{namcfg} 331 332 Or they can be defined in an analytical way in \path{MY_SRC} directory of the configuration. 333 For Reference Configurations of NEMO input domain files are supplied by NEMO System Team. 334 For analytical definition of input fields two routines are supplied: \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}. 335 They are an example of GYRE configuration parameters, and they are available in \path{src/OCE/USR} directory, 336 they provide the horizontal and vertical mesh. 337 % ------------------------------------------------------------------------------------------------------------- 338 % Needed fields 339 % ------------------------------------------------------------------------------------------------------------- 340 %\subsection{List of needed fields to build DOMAIN} 341 %\label{subsec:DOM_fields_list} 342 375 The total size of the computational domain is set by the parameters 376 \np{jpiglo}, \np{jpjglo} and \np{jpkglo} for the $i$, $j$ and $k$ 377 directions, respectively. Note, that the variables \forcode{jpi} and \forcode{jpj} 378 refer to the size of each processor subdomain when the code is run in 379 parallel using domain decomposition (\key{mpp\_mpi} defined, see 380 \autoref{sec:LBC_mpp}). 381 382 The name of the configuration is set through parameter \np{cn\_cfg}, 383 and the nominal resolution through parameter \np{nn\_cfg} (unless in 384 the input file both of variables \forcode{ORCA} and \forcode{ORCA_index} 385 are present, in which case \np{cn\_cfg} and \np{nn\_cfg} are set from these 386 values accordingly). 387 388 The global lateral boundary condition type is selected from 8 options 389 using parameter \np{jperio}. See \autoref{sec:LBC_jperio} for 390 details on the available options and the corresponding values for 391 \np{jperio}. 343 392 344 393 % ================================================================ 345 394 % Domain: Horizontal Grid (mesh) 346 395 % ================================================================ 347 \section{Horizontal grid mesh (\protect\mdl{domhgr})} 348 \label{sec:DOM_hgr} 349 350 % ------------------------------------------------------------------------------------------------------------- 351 % Coordinates and scale factors 352 % ------------------------------------------------------------------------------------------------------------- 353 \subsection{Coordinates and scale factors} 354 \label{subsec:DOM_hgr_coord_e} 355 356 The ocean mesh (\ie the position of all the scalar and vector points) is defined by 357 the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 358 The grid-points are located at integer or integer and a half values of as indicated in \autoref{tab:cell}. 359 The associated scale factors are defined using the analytical first derivative of the transformation 360 \autoref{eq:scale_factors}. 361 These definitions are done in two modules, \mdl{domhgr} and \mdl{domzgr}, 362 which provide the horizontal and vertical meshes, respectively. 363 This section deals with the horizontal mesh parameters. 364 365 In a horizontal plane, the location of all the model grid points is defined from 366 the analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a function of $(i,j)$. 367 The horizontal scale factors are calculated using \autoref{eq:scale_factors}. 368 For example, when the longitude and latitude are function of a single value 369 ($i$ and $j$, respectively) (geographical configuration of the mesh), 370 the horizontal mesh definition reduces to define the wanted $\lambda(i)$, $\varphi(j)$, 371 and their derivatives $\lambda'(i) \ \varphi'(j)$ in the \mdl{domhgr} module. 372 The model computes the grid-point positions and scale factors in the horizontal plane as follows: 373 \begin{align*} 374 \lambda_t &\equiv \text{glamt} = \lambda (i ) 375 &\varphi_t &\equiv \text{gphit} = \varphi (j ) \\ 376 \lambda_u &\equiv \text{glamu} = \lambda (i + 1/2) 377 &\varphi_u &\equiv \text{gphiu} = \varphi (j ) \\ 378 \lambda_v &\equiv \text{glamv} = \lambda (i ) 379 &\varphi_v &\equiv \text{gphiv} = \varphi (j + 1/2) \\ 380 \lambda_f &\equiv \text{glamf} = \lambda (i + 1/2) 381 &\varphi_f &\equiv \text{gphif} = \varphi (j + 1/2) \\ 382 e_{1t} &\equiv \text{e1t} = r_a |\lambda'(i ) \; \cos\varphi(j ) | 383 &e_{2t} &\equiv \text{e2t} = r_a |\varphi'(j ) | \\ 384 e_{1u} &\equiv \text{e1t} = r_a |\lambda'(i + 1/2) \; \cos\varphi(j ) | 385 &e_{2u} &\equiv \text{e2t} = r_a |\varphi'(j ) | \\ 386 e_{1v} &\equiv \text{e1t} = r_a |\lambda'(i ) \; \cos\varphi(j + 1/2) | 387 &e_{2v} &\equiv \text{e2t} = r_a |\varphi'(j + 1/2) | \\ 388 e_{1f} &\equiv \text{e1t} = r_a |\lambda'(i + 1/2) \; \cos\varphi(j + 1/2) | 389 &e_{2f} &\equiv \text{e2t} = r_a |\varphi'(j + 1/2) | 390 \end{align*} 391 where the last letter of each computational name indicates the grid point considered and 392 $r_a$ is the earth radius (defined in \mdl{phycst} along with all universal constants). 393 Note that the horizontal position of and scale factors at $w$-points are exactly equal to those of $t$-points, 394 thus no specific arrays are defined at $w$-points. 395 396 Note that the definition of the scale factors 397 (\ie as the analytical first derivative of the transformation that 398 gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) 399 is specific to the \NEMO model \citep{Marti_al_JGR92}. 400 As an example, $e_{1t}$ is defined locally at a $t$-point, 401 whereas many other models on a C grid choose to define such a scale factor as 402 the distance between the $U$-points on each side of the $t$-point. 403 Relying on an analytical transformation has two advantages: 404 firstly, there is no ambiguity in the scale factors appearing in the discrete equations, 405 since they are first introduced in the continuous equations; 406 secondly, analytical transformations encourage good practice by the definition of smoothly varying grids 407 (rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{Treguier1996}. 408 An example of the effect of such a choice is shown in \autoref{fig:zgr_e3}. 409 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 410 \begin{figure}[!t] 411 \begin{center} 412 \includegraphics[]{Fig_zgr_e3} 413 \caption{ 414 \protect\label{fig:zgr_e3} 415 Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, 416 and (b) analytically derived grid-point position and scale factors. 417 For both grids here, the same $w$-point depth has been chosen but 418 in (a) the $t$-points are set half way between $w$-points while 419 in (b) they are defined from an analytical function: 420 $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$. 421 Note the resulting difference between the value of the grid-size $\Delta_k$ and 422 those of the scale factor $e_k$. 423 } 424 \end{center} 425 \end{figure} 426 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 427 428 % ------------------------------------------------------------------------------------------------------------- 429 % Choice of horizontal grid 430 % ------------------------------------------------------------------------------------------------------------- 431 \subsection{Choice of horizontal grid} 432 \label{subsec:DOM_hgr_msh_choice} 433 434 % ------------------------------------------------------------------------------------------------------------- 435 % Grid files 436 % ------------------------------------------------------------------------------------------------------------- 437 \subsection{Output grid files} 438 \label{subsec:DOM_hgr_files} 439 440 All the arrays relating to a particular ocean model configuration (grid-point position, scale factors, masks) 441 can be saved in files if \np{nn\_msh} $\not = 0$ (namelist variable in \ngn{namdom}). 442 This can be particularly useful for plots and off-line diagnostics. 443 In some cases, the user may choose to make a local modification of a scale factor in the code. 444 This is the case in global configurations when restricting the width of a specific strait 445 (usually a one-grid-point strait that happens to be too wide due to insufficient model resolution). 446 An example is Gibraltar Strait in the ORCA2 configuration. 447 When such modifications are done, 448 the output grid written when \np{nn\_msh} $\not = 0$ is no more equal to the input grid. 396 \subsection{Horizontal grid mesh (\protect\mdl{domhgr})} 397 \label{subsec:DOM_hgr} 398 399 % ================================================================ 400 % Domain: List of hgr-related fields needed 401 % ================================================================ 402 \subsubsection{Required fields} 403 \label{sec:DOM_hgr_fields} 404 The explicit specification of a range of mesh-related fields are required for the definition of a configuration. These include: 405 406 \begin{Verbatim}[fontsize=\tiny] 407 int jpiglo, jpjglo, jpkglo /* global domain sizes */ 408 int jperio /* lateral global domain b.c. */ 409 double glamt, glamu, glamv, glamf /* geographic longitude (t,u,v and f points respectively) */ 410 double gphit, gphiu, gphiv, gphif /* geographic latitude */ 411 double e1t, e1u, e1v, e1f /* horizontal scale factors */ 412 double e2t, e2u, e2v, e2f /* horizontal scale factors */ 413 \end{Verbatim} 414 415 The values of the geographic longitude and latitude arrays at indices $i,j$ correspond to the analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a function of $(i,j)$, evaluated at the values as specified in Table \autoref{tab:cell} for the respective grid-point position. The calculation of the values of the horizontal scale factor arrays in general additionally involves partial derivatives of $\lambda$ and $\varphi$ with respect to $i$ and $j$, evaluated for the same arguments as $\lambda$ and $\varphi$. 416 417 \subsubsection{Optional fields} 418 \begin{Verbatim}[fontsize=\tiny] 419 /* Optional: */ 420 int ORCA, ORCA_index /* configuration name, configuration resolution */ 421 double e1e2u, e1e2v /* U and V surfaces (if grid size reduction in some straits) */ 422 double ff_f, ff_t /* Coriolis parameter (if not on the sphere) */ 423 \end{Verbatim} 424 425 NEMO can support the local reduction of key strait widths by altering individual values of 426 e2u or e1v at the appropriate locations. This is particularly useful for locations such as 427 Gibraltar or Indonesian Throughflow pinch-points (see \autoref{sec:MISC_strait} for 428 illustrated examples). The key is to reduce the faces of $T$-cell (\ie change the value of 429 the horizontal scale factors at $u$- or $v$-point) but not the volume of the cells. Doing 430 otherwise can lead to numerical instability issues. In normal operation the surface areas 431 are computed from $\texttt{e1u} * \texttt{e2u}$ and $\texttt{e1v} * \texttt{e2v}$ but in 432 cases where a gridsize reduction is required, the unaltered surface areas at $u$ and $v$ 433 grid points (\texttt{e1e2u} and \texttt{e1e2v}, respectively) must be read or pre-computed 434 in \mdl{usrdef\_hgr}. If these arrays are present in the \np{cn\_domcfg} file they are 435 read and the internal computation is suppressed. Versions of \mdl{usrdef\_hgr} which set 436 their own values of \texttt{e1e2u} and \texttt{e1e2v} should set the surface-area 437 computation flag: \texttt{ie1e2u\_v} to a non-zero value to suppress their re-computation. 438 439 \smallskip 440 Similar logic applies to the other optional fields: \texttt{ff\_f} and \texttt{ff\_t} 441 which can be used to provide the Coriolis parameter at F- and T-points respectively if the 442 mesh is not on a sphere. If present these fields will be read and used and the normal 443 calculation ($2*\Omega*\sin(\varphi)$) suppressed. Versions of \mdl{usrdef\_hgr} which set 444 their own values of \texttt{ff\_f} and \texttt{ff\_t} should set the Coriolis computation 445 flag: \texttt{iff} to a non-zero value to suppress their re-computation. 446 447 Note that longitudes, latitudes, and scale factors at $w$ points are exactly 448 equal to those of $t$ points, thus no specific arrays are defined at $w$ points. 449 449 450 450 451 % ================================================================ 451 452 % Domain: Vertical Grid (domzgr) 452 453 % ================================================================ 453 \section{Vertical grid (\protect\mdl{domzgr})} 454 \label{sec:DOM_zgr} 455 %-----------------------------------------nam_zgr & namdom------------------------------------------- 456 % 457 %\nlst{namzgr} 458 459 \nlst{namdom} 454 \subsection[Vertical grid (\textit{domzgr.F90})] 455 {Vertical grid (\protect\mdl{domzgr})} 456 \label{subsec:DOM_zgr} 457 %-----------------------------------------namdom------------------------------------------- 458 \nlst{namdom} 460 459 %------------------------------------------------------------------------------------------------------------- 461 460 462 Variables are defined through the \ngn{namzgr} and \ngn{namdom} namelists.463 461 In the vertical, the model mesh is determined by four things: 464 (1) the bathymetry given in meters; 465 (2) the number of levels of the model (\jp{jpk}); 466 (3) the analytical transformation $z(i,j,k)$ and the vertical scale factors (derivatives of the transformation); and 467 (4) the masking system, \ie the number of wet model levels at each 468 $(i,j)$ column of points. 462 \begin{enumerate} 463 \item the bathymetry given in meters; 464 \item the number of levels of the model (\jp{jpk}); 465 \item the analytical transformation $z(i,j,k)$ and the vertical scale factors (derivatives of the transformation); and 466 \item the masking system, \ie the number of wet model levels at each 467 $(i,j)$ location of the horizontal grid. 468 \end{enumerate} 469 469 470 470 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 471 471 \begin{figure}[!tb] 472 472 \begin{center} 473 \includegraphics[ ]{Fig_z_zps_s_sps}473 \includegraphics[width=\textwidth]{Fig_z_zps_s_sps} 474 474 \caption{ 475 475 \protect\label{fig:z_zps_s_sps} … … 480 480 (d) hybrid $s-z$ coordinate, 481 481 (e) hybrid $s-z$ coordinate with partial step, and 482 (f) same as (e) but in the non-linear free surface (\protect\np{ln\_linssh} ~\forcode{= .false.}).482 (f) same as (e) but in the non-linear free surface (\protect\np{ln\_linssh}\forcode{ = .false.}). 483 483 Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e). 484 484 } … … 487 487 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 488 488 489 The choice of a vertical coordinate, even if it is made through \ngn{namzgr} namelist parameters, 490 must be done once of all at the beginning of an experiment. 491 It is not intended as an option which can be enabled or disabled in the middle of an experiment. 492 Three main choices are offered (\autoref{fig:z_zps_s_sps}): 493 $z$-coordinate with full step bathymetry (\np{ln\_zco}~\forcode{= .true.}), 494 $z$-coordinate with partial step bathymetry (\np{ln\_zps}~\forcode{= .true.}), 495 or generalized, $s$-coordinate (\np{ln\_sco}~\forcode{= .true.}). 496 Hybridation of the three main coordinates are available: 497 $s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps} and \autoref{fig:z_zps_s_sps}). 498 By default a non-linear free surface is used: the coordinate follow the time-variation of the free surface so that 499 the transformation is time dependent: $z(i,j,k,t)$ (\autoref{fig:z_zps_s_sps}). 500 When a linear free surface is assumed (\np{ln\_linssh}~\forcode{= .true.}), 501 the vertical coordinate are fixed in time, but the seawater can move up and down across the $z_0$ surface 502 (in other words, the top of the ocean in not a rigid-lid). 503 The last choice in terms of vertical coordinate concerns the presence (or not) in 504 the model domain of ocean cavities beneath ice shelves. 505 Setting \np{ln\_isfcav} to true allows to manage ocean cavities, otherwise they are filled in. 506 This option is currently only available in $z$- or $zps$-coordinate, 507 and partial step are also applied at the ocean/ice shelf interface. 508 509 Contrary to the horizontal grid, the vertical grid is computed in the code and no provision is made for 510 reading it from a file. 511 The only input file is the bathymetry (in meters) (\ifile{bathy\_meter}) 512 \footnote{ 513 N.B. in full step $z$-coordinate, a \ifile{bathy\_level} file can replace the \ifile{bathy\_meter} file, 514 so that the computation of the number of wet ocean point in each water column is by-passed}. 515 If \np{ln\_isfcav}~\forcode{= .true.}, an extra file input file (\ifile{isf\_draft\_meter}) describing 516 the ice shelf draft (in meters) is needed. 517 518 After reading the bathymetry, the algorithm for vertical grid definition differs between the different options: 519 \begin{description} 520 \item[\textit{zco}] 521 set a reference coordinate transformation $z_0(k)$, and set $z(i,j,k,t) = z_0(k)$. 522 \item[\textit{zps}] 523 set a reference coordinate transformation $z_0(k)$, and calculate the thickness of the deepest level at 524 each $(i,j)$ point using the bathymetry, to obtain the final three-dimensional depth and scale factor arrays. 525 \item[\textit{sco}] 526 smooth the bathymetry to fulfill the hydrostatic consistency criteria and 527 set the three-dimensional transformation. 528 \item[\textit{s-z} and \textit{s-zps}] 529 smooth the bathymetry to fulfill the hydrostatic consistency criteria and 530 set the three-dimensional transformation $z(i,j,k)$, 531 and possibly introduce masking of extra land points to better fit the original bathymetry file. 532 \end{description} 533 %%% 534 \gmcomment{ add the description of the smoothing: envelop topography...} 535 %%% 536 537 Unless a linear free surface is used (\np{ln\_linssh}~\forcode{= .false.}), 538 the arrays describing the grid point depths and vertical scale factors are three set of 539 three dimensional arrays $(i,j,k)$ defined at \textit{before}, \textit{now} and \textit{after} time step. 540 The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively. 541 They are updated at each model time step using a fixed reference coordinate system which 542 computer names have a $\_0$ suffix. 543 When the linear free surface option is used (\np{ln\_linssh}~\forcode{= .true.}), \textit{before}, 544 \textit{now} and \textit{after} arrays are simply set one for all to their reference counterpart. 545 546 % ------------------------------------------------------------------------------------------------------------- 547 % Meter Bathymetry 548 % ------------------------------------------------------------------------------------------------------------- 549 \subsection{Meter bathymetry} 550 \label{subsec:DOM_bathy} 551 552 Three options are possible for defining the bathymetry, according to the namelist variable \np{nn\_bathy} 553 (found in \ngn{namdom} namelist): 554 \begin{description} 555 \item[\np{nn\_bathy}~\forcode{= 0}]: 556 a flat-bottom domain is defined. 557 The total depth $z_w (jpk)$ is given by the coordinate transformation. 558 The domain can either be a closed basin or a periodic channel depending on the parameter \np{jperio}. 559 \item[\np{nn\_bathy}~\forcode{= -1}]: 560 a domain with a bump of topography one third of the domain width at the central latitude. 561 This is meant for the "EEL-R5" configuration, a periodic or open boundary channel with a seamount. 562 \item[\np{nn\_bathy}~\forcode{= 1}]: 563 read a bathymetry and ice shelf draft (if needed). 564 The \ifile{bathy\_meter} file (Netcdf format) provides the ocean depth (positive, in meters) at 565 each grid point of the model grid. 566 The bathymetry is usually built by interpolating a standard bathymetry product (\eg ETOPO2) onto 567 the horizontal ocean mesh. 568 Defining the bathymetry also defines the coastline: where the bathymetry is zero, 569 no model levels are defined (all levels are masked). 570 571 The \ifile{isfdraft\_meter} file (Netcdf format) provides the ice shelf draft (positive, in meters) at 572 each grid point of the model grid. 573 This file is only needed if \np{ln\_isfcav}~\forcode{= .true.}. 574 Defining the ice shelf draft will also define the ice shelf edge and the grounding line position. 575 \end{description} 576 577 When a global ocean is coupled to an atmospheric model it is better to represent all large water bodies 578 (\eg great lakes, Caspian sea...) even if the model resolution does not allow their communication with 579 the rest of the ocean. 580 This is unnecessary when the ocean is forced by fixed atmospheric conditions, 581 so these seas can be removed from the ocean domain. 582 The user has the option to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}), 583 but the code has to be adapted to the user's configuration. 584 585 % ------------------------------------------------------------------------------------------------------------- 586 % z-coordinate and reference coordinate transformation 587 % ------------------------------------------------------------------------------------------------------------- 588 \subsection[$Z$-coordinate (\protect\np{ln\_zco}~\forcode{= .true.}) and ref. coordinate] 589 {$Z$-coordinate (\protect\np{ln\_zco}~\forcode{= .true.}) and reference coordinate} 590 \label{subsec:DOM_zco} 591 592 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 593 \begin{figure}[!tb] 594 \begin{center} 595 \includegraphics[]{Fig_zgr} 596 \caption{ 597 \protect\label{fig:zgr} 598 Default vertical mesh for ORCA2: 30 ocean levels (L30). 599 Vertical level functions for (a) T-point depth and (b) the associated scale factor as computed from 600 \autoref{eq:DOM_zgr_ana_1} using \autoref{eq:DOM_zgr_coef} in $z$-coordinate. 601 } 602 \end{center} 603 \end{figure} 604 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 605 606 The reference coordinate transformation $z_0(k)$ defines the arrays $gdept_0$ and $gdepw_0$ for $t$- and $w$-points, 607 respectively. 608 As indicated on \autoref{fig:index_vert} \jp{jpk} is the number of $w$-levels. 609 $gdepw_0(1)$ is the ocean surface. 610 There are at most \jp{jpk}-1 $t$-points inside the ocean, 611 the additional $t$-point at $jk = jpk$ is below the sea floor and is not used. 612 The vertical location of $w$- and $t$-levels is defined from the analytic expression of the depth $z_0(k)$ whose 613 analytical derivative with respect to $k$ provides the vertical scale factors. 614 The user must provide the analytical expression of both $z_0$ and its first derivative with respect to $k$. 615 This is done in routine \mdl{domzgr} through statement functions, 616 using parameters provided in the \ngn{namcfg} namelist. 617 618 It is possible to define a simple regular vertical grid by giving zero stretching (\np{ppacr}~\forcode{= 0}). 619 In that case, the parameters \jp{jpk} (number of $w$-levels) and 620 \np{pphmax} (total ocean depth in meters) fully define the grid. 621 622 For climate-related studies it is often desirable to concentrate the vertical resolution near the ocean surface. 623 The following function is proposed as a standard for a $z$-coordinate (with either full or partial steps): 624 \begin{gather} 625 \label{eq:DOM_zgr_ana_1} 626 z_0 (k) = h_{sur} - h_0 \; k - \; h_1 \; \log \big[ \cosh ((k - h_{th}) / h_{cr}) \big] \\ 627 e_3^0(k) = \lt| - h_0 - h_1 \; \tanh \big[ (k - h_{th}) / h_{cr} \big] \rt| 628 \end{gather} 629 where $k = 1$ to \jp{jpk} for $w$-levels and $k = 1$ to $k = 1$ for $T-$levels. 630 Such an expression allows us to define a nearly uniform vertical location of levels at the ocean top and bottom with 631 a smooth hyperbolic tangent transition in between (\autoref{fig:zgr}). 632 633 If the ice shelf cavities are opened (\np{ln\_isfcav}~\forcode{= .true.}), the definition of $z_0$ is the same. 634 However, definition of $e_3^0$ at $t$- and $w$-points is respectively changed to: 635 \begin{equation} 636 \label{eq:DOM_zgr_ana_2} 637 \begin{split} 638 e_3^T(k) &= z_W (k + 1) - z_W (k ) \\ 639 e_3^W(k) &= z_T (k ) - z_T (k - 1) 640 \end{split} 641 \end{equation} 642 This formulation decrease the self-generated circulation into the ice shelf cavity 643 (which can, in extreme case, leads to blow up).\\ 644 645 The most used vertical grid for ORCA2 has $10~m$ ($500~m$) resolution in the surface (bottom) layers and 646 a depth which varies from 0 at the sea surface to a minimum of $-5000~m$. 647 This leads to the following conditions: 648 \begin{equation} 649 \label{eq:DOM_zgr_coef} 650 \begin{array}{ll} 651 e_3 (1 + 1/2) = 10. & z(1 ) = 0. \\ 652 e_3 (jpk - 1/2) = 500. & z(jpk) = -5000. 653 \end{array} 654 \end{equation} 655 656 With the choice of the stretching $h_{cr} = 3$ and the number of levels \jp{jpk}~$= 31$, 657 the four coefficients $h_{sur}$, $h_0$, $h_1$, and $h_{th}$ in 658 \autoref{eq:DOM_zgr_ana_2} have been determined such that 659 \autoref{eq:DOM_zgr_coef} is satisfied, through an optimisation procedure using a bisection method. 660 For the first standard ORCA2 vertical grid this led to the following values: 661 $h_{sur} = 4762.96$, $h_0 = 255.58, h_1 = 245.5813$, and $h_{th} = 21.43336$. 662 The resulting depths and scale factors as a function of the model levels are shown in 663 \autoref{fig:zgr} and given in \autoref{tab:orca_zgr}. 664 Those values correspond to the parameters \np{ppsur}, \np{ppa0}, \np{ppa1}, \np{ppkth} in \ngn{namcfg} namelist. 665 666 Rather than entering parameters $h_{sur}$, $h_0$, and $h_1$ directly, it is possible to recalculate them. 667 In that case the user sets \np{ppsur}~$=$~\np{ppa0}~$=$~\np{ppa1}~$= 999999$., 668 in \ngn{namcfg} namelist, and specifies instead the four following parameters: 489 The choice of a vertical coordinate is made when setting up the configuration; 490 it is not intended to be an option which can be changed in the middle of an 491 experiment. The one exception to this statement being the choice of linear or 492 non-linear free surface. In v4.0 the linear free surface option is implemented 493 as a special case of the non-linear free surface. This is computationally 494 wasteful since it uses the structures for time-varying 3D metrics for fields 495 that (in the linear free surface case) are fixed. However, the linear 496 free-surface is rarely used and implementing it this way means a single configuration 497 file can support both options. 498 499 By default a non-linear free surface is used (\np{ln\_linssh} set to \forcode{ = 500 .false.} in \ngn{namdom}): the coordinate follow the time-variation of the free 501 surface so that the transformation is time dependent: $z(i,j,k,t)$ 502 (\eg \autoref{fig:z_zps_s_sps}f). When a linear free surface is assumed 503 (\np{ln\_linssh} set to \forcode{ = .true.} in \ngn{namdom}), the vertical 504 coordinates are fixed in time, but the seawater can move up and down across the 505 $z_0$ surface (in other words, the top of the ocean in not a rigid lid). 506 507 Note that settings: \np{ln\_zco}, \np{ln\_zps}, \np{ln\_sco} and \np{ln\_isfcav} mentioned 508 in the following sections appear to be namelist options but they are no longer truly 509 namelist options for NEMO. Their value is written to and read from the domain configuration file 510 and they should be treated as fixed parameters for a particular configuration. They are 511 namelist options for the \forcode{DOMAINcfg} tool that can be used to build the 512 configuration file and serve both to provide a record of the choices made whilst building the 513 configuration and to trigger appropriate code blocks within NEMO. 514 These values should not be altered in the \np{cn\_domcfg} file. 515 516 \medskip 517 The decision on these choices must be made when the \np{cn\_domcfg} file is constructed. 518 Three main choices are offered (\autoref{fig:z_zps_s_sps}a-c): 519 669 520 \begin{itemize} 670 \item 671 \np{ppacr}~$= h_{cr}$: stretching factor (nondimensional). 672 The larger \np{ppacr}, the smaller the stretching. 673 Values from $3$ to $10$ are usual. 674 \item 675 \np{ppkth}~$= h_{th}$: is approximately the model level at which maximum stretching occurs 676 (nondimensional, usually of order 1/2 or 2/3 of \jp{jpk}) 677 \item 678 \np{ppdzmin}: minimum thickness for the top layer (in meters). 679 \item 680 \np{pphmax}: total depth of the ocean (meters). 521 \item $z$-coordinate with full step bathymetry (\np{ln\_zco}\forcode{ = .true.}), 522 \item $z$-coordinate with partial step ($zps$) bathymetry (\np{ln\_zps}\forcode{ = .true.}), 523 \item Generalized, $s$-coordinate (\np{ln\_sco}\forcode{ = .true.}). 681 524 \end{itemize} 682 As an example, for the $45$ layers used in the DRAKKAR configuration those parameters are: 683 \jp{jpk}~$= 46$, \np{ppacr}~$= 9$, \np{ppkth}~$= 23.563$, \np{ppdzmin}~$= 6~m$, \np{pphmax}~$= 5750~m$. 684 685 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 686 \begin{table} 687 \begin{center} 688 \begin{tabular}{c||r|r|r|r} 689 \hline 690 \textbf{LEVEL} & \textbf{gdept\_1d} & \textbf{gdepw\_1d} & \textbf{e3t\_1d } & \textbf{e3w\_1d} \\ 691 \hline 692 1 & \textbf{ 5.00} & 0.00 & \textbf{ 10.00} & 10.00 \\ 693 \hline 694 2 & \textbf{ 15.00} & 10.00 & \textbf{ 10.00} & 10.00 \\ 695 \hline 696 3 & \textbf{ 25.00} & 20.00 & \textbf{ 10.00} & 10.00 \\ 697 \hline 698 4 & \textbf{ 35.01} & 30.00 & \textbf{ 10.01} & 10.00 \\ 699 \hline 700 5 & \textbf{ 45.01} & 40.01 & \textbf{ 10.01} & 10.01 \\ 701 \hline 702 6 & \textbf{ 55.03} & 50.02 & \textbf{ 10.02} & 10.02 \\ 703 \hline 704 7 & \textbf{ 65.06} & 60.04 & \textbf{ 10.04} & 10.03 \\ 705 \hline 706 8 & \textbf{ 75.13} & 70.09 & \textbf{ 10.09} & 10.06 \\ 707 \hline 708 9 & \textbf{ 85.25} & 80.18 & \textbf{ 10.17} & 10.12 \\ 709 \hline 710 10 & \textbf{ 95.49} & 90.35 & \textbf{ 10.33} & 10.24 \\ 711 \hline 712 11 & \textbf{ 105.97} & 100.69 & \textbf{ 10.65} & 10.47 \\ 713 \hline 714 12 & \textbf{ 116.90} & 111.36 & \textbf{ 11.27} & 10.91 \\ 715 \hline 716 13 & \textbf{ 128.70} & 122.65 & \textbf{ 12.47} & 11.77 \\ 717 \hline 718 14 & \textbf{ 142.20} & 135.16 & \textbf{ 14.78} & 13.43 \\ 719 \hline 720 15 & \textbf{ 158.96} & 150.03 & \textbf{ 19.23} & 16.65 \\ 721 \hline 722 16 & \textbf{ 181.96} & 169.42 & \textbf{ 27.66} & 22.78 \\ 723 \hline 724 17 & \textbf{ 216.65} & 197.37 & \textbf{ 43.26} & 34.30 \\ 725 \hline 726 18 & \textbf{ 272.48} & 241.13 & \textbf{ 70.88} & 55.21 \\ 727 \hline 728 19 & \textbf{ 364.30} & 312.74 & \textbf{ 116.11} & 90.99 \\ 729 \hline 730 20 & \textbf{ 511.53} & 429.72 & \textbf{ 181.55} & 146.43 \\ 731 \hline 732 21 & \textbf{ 732.20} & 611.89 & \textbf{ 261.03} & 220.35 \\ 733 \hline 734 22 & \textbf{ 1033.22} & 872.87 & \textbf{ 339.39} & 301.42 \\ 735 \hline 736 23 & \textbf{ 1405.70} & 1211.59 & \textbf{ 402.26} & 373.31 \\ 737 \hline 738 24 & \textbf{ 1830.89} & 1612.98 & \textbf{ 444.87} & 426.00 \\ 739 \hline 740 25 & \textbf{ 2289.77} & 2057.13 & \textbf{ 470.55} & 459.47 \\ 741 \hline 742 26 & \textbf{ 2768.24} & 2527.22 & \textbf{ 484.95} & 478.83 \\ 743 \hline 744 27 & \textbf{ 3257.48} & 3011.90 & \textbf{ 492.70} & 489.44 \\ 745 \hline 746 28 & \textbf{ 3752.44} & 3504.46 & \textbf{ 496.78} & 495.07 \\ 747 \hline 748 29 & \textbf{ 4250.40} & 4001.16 & \textbf{ 498.90} & 498.02 \\ 749 \hline 750 30 & \textbf{ 4749.91} & 4500.02 & \textbf{ 500.00} & 499.54 \\ 751 \hline 752 31 & \textbf{ 5250.23} & 5000.00 & \textbf{ 500.56} & 500.33 \\ 753 \hline 754 \end{tabular} 755 \end{center} 756 \caption{ 757 \protect\label{tab:orca_zgr} 758 Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration as computed from 759 \autoref{eq:DOM_zgr_ana_2} using the coefficients given in \autoref{eq:DOM_zgr_coef} 760 } 761 \end{table} 762 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 763 764 % ------------------------------------------------------------------------------------------------------------- 765 % z-coordinate with partial step 766 % ------------------------------------------------------------------------------------------------------------- 767 \subsection{$Z$-coordinate with partial step (\protect\np{ln\_zps}~\forcode{= .true.})} 768 \label{subsec:DOM_zps} 769 %--------------------------------------------namdom------------------------------------------------------- 770 771 \nlst{namdom} 772 %-------------------------------------------------------------------------------------------------------------- 773 774 In $z$-coordinate partial step, 775 the depths of the model levels are defined by the reference analytical function $z_0(k)$ as described in 776 the previous section, \textit{except} in the bottom layer. 777 The thickness of the bottom layer is allowed to vary as a function of geographical location $(\lambda,\varphi)$ to 778 allow a better representation of the bathymetry, especially in the case of small slopes 779 (where the bathymetry varies by less than one level thickness from one grid point to the next). 780 The reference layer thicknesses $e_{3t}^0$ have been defined in the absence of bathymetry. 781 With partial steps, layers from 1 to \jp{jpk}-2 can have a thickness smaller than $e_{3t}(jk)$. 782 The model deepest layer (\jp{jpk}-1) is allowed to have either a smaller or larger thickness than $e_{3t}(jpk)$: 783 the maximum thickness allowed is $2*e_{3t}(jpk - 1)$. 784 This has to be kept in mind when specifying values in \ngn{namdom} namelist, 785 as the maximum depth \np{pphmax} in partial steps: 786 for example, with \np{pphmax}~$= 5750~m$ for the DRAKKAR 45 layer grid, 787 the maximum ocean depth allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk - 1)$ being $250~m$). 788 Two variables in the namdom namelist are used to define the partial step vertical grid. 789 The mimimum water thickness (in meters) allowed for a cell partially filled with bathymetry at level jk is 790 the minimum of \np{rn\_e3zps\_min} (thickness in meters, usually $20~m$) or $e_{3t}(jk)*$\np{rn\_e3zps\_rat} 791 (a fraction, usually 10\%, of the default thickness $e_{3t}(jk)$). 792 793 \gmcomment{ \colorbox{yellow}{Add a figure here of pstep especially at last ocean level } } 794 795 % ------------------------------------------------------------------------------------------------------------- 796 % s-coordinate 797 % ------------------------------------------------------------------------------------------------------------- 798 \subsection{$S$-coordinate (\protect\np{ln\_sco}~\forcode{= .true.})} 799 \label{subsec:DOM_sco} 800 %------------------------------------------nam_zgr_sco--------------------------------------------------- 801 % 802 %\nlst{namzgr_sco} 803 %-------------------------------------------------------------------------------------------------------------- 804 Options are defined in \ngn{namzgr\_sco}. 805 In $s$-coordinate (\np{ln\_sco}~\forcode{= .true.}), the depth and thickness of the model levels are defined from 806 the product of a depth field and either a stretching function or its derivative, respectively: 807 808 \begin{align*} 809 % \label{eq:DOM_sco_ana} 810 z(k) &= h(i,j) \; z_0 (k) \\ 811 e_3(k) &= h(i,j) \; z_0'(k) 812 \end{align*} 813 814 where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $t$-point location in the horizontal and 815 $z_0(k)$ is a function which varies from $0$ at the sea surface to $1$ at the ocean bottom. 816 The depth field $h$ is not necessary the ocean depth, 817 since a mixed step-like and bottom-following representation of the topography can be used 818 (\autoref{fig:z_zps_s_sps}) or an envelop bathymetry can be defined (\autoref{fig:z_zps_s_sps}). 819 The namelist parameter \np{rn\_rmax} determines the slope at which 820 the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate. 821 The coordinate can also be hybridised by specifying \np{rn\_sbot\_min} and \np{rn\_sbot\_max} as 822 the minimum and maximum depths at which the terrain-following vertical coordinate is calculated. 823 824 Options for stretching the coordinate are provided as examples, 825 but care must be taken to ensure that the vertical stretch used is appropriate for the application. 826 827 The original default NEMO s-coordinate stretching is available if neither of the other options are specified as true 828 (\np{ln\_s\_SH94}~\forcode{= .false.} and \np{ln\_s\_SF12}~\forcode{= .false.}). 829 This uses a depth independent $\tanh$ function for the stretching \citep{Madec_al_JPO96}: 830 831 \[ 832 z = s_{min} + C (s) (H - s_{min}) 833 % \label{eq:SH94_1} 834 \] 835 836 where $s_{min}$ is the depth at which the $s$-coordinate stretching starts and 837 allows a $z$-coordinate to placed on top of the stretched coordinate, 838 and $z$ is the depth (negative down from the asea surface). 839 \begin{gather*} 840 s = - \frac{k}{n - 1} \quad \text{and} \quad 0 \leq k \leq n - 1 841 % \label{eq:DOM_s} 842 \\ 843 % \label{eq:DOM_sco_function} 844 C(s) = \frac{[\tanh(\theta \, (s + b)) - \tanh(\theta \, b)]}{2 \; \sinh(\theta)} 845 \end{gather*} 846 847 A stretching function, 848 modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_s\_SH94}~\forcode{= .true.}), 849 is also available and is more commonly used for shelf seas modelling: 850 851 \[ 852 C(s) = (1 - b) \frac{\sinh(\theta s)}{\sinh(\theta)} 853 + b \frac{\tanh \lt[ \theta \lt(s + \frac{1}{2} \rt) \rt] - \tanh \lt( \frac{\theta}{2} \rt)} 854 { 2 \tanh \lt( \frac{\theta}{2} \rt)} 855 % \label{eq:SH94_2} 856 \] 857 858 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 859 \begin{figure}[!ht] 860 \begin{center} 861 \includegraphics[]{Fig_sco_function} 862 \caption{ 863 \protect\label{fig:sco_function} 864 Examples of the stretching function applied to a seamount; 865 from left to right: surface, surface and bottom, and bottom intensified resolutions 866 } 867 \end{center} 868 \end{figure} 869 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 870 871 where $H_c$ is the critical depth (\np{rn\_hc}) at which the coordinate transitions from pure $\sigma$ to 872 the stretched coordinate, and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) are the surface and 873 bottom control parameters such that $0 \leqslant \theta \leqslant 20$, and $0 \leqslant b \leqslant 1$. 874 $b$ has been designed to allow surface and/or bottom increase of the vertical resolution 875 (\autoref{fig:sco_function}). 876 877 Another example has been provided at version 3.5 (\np{ln\_s\_SF12}) that allows a fixed surface resolution in 878 an analytical terrain-following stretching \citet{Siddorn_Furner_OM12}. 879 In this case the a stretching function $\gamma$ is defined such that: 880 881 \begin{equation} 882 z = - \gamma h \quad \text{with} \quad 0 \leq \gamma \leq 1 883 % \label{eq:z} 884 \end{equation} 885 886 The function is defined with respect to $\sigma$, the unstretched terrain-following coordinate: 887 888 \begin{gather*} 889 % \label{eq:DOM_gamma_deriv} 890 \gamma = A \lt( \sigma - \frac{1}{2} (\sigma^2 + f (\sigma)) \rt) 891 + B \lt( \sigma^3 - f (\sigma) \rt) + f (\sigma) \\ 892 \intertext{Where:} 893 % \label{eq:DOM_gamma} 894 f(\sigma) = (\alpha + 2) \sigma^{\alpha + 1} - (\alpha + 1) \sigma^{\alpha + 2} 895 \quad \text{and} \quad \sigma = \frac{k}{n - 1} 896 \end{gather*} 897 898 This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of 899 the user prescribed stretching parameter $\alpha$ (\np{rn\_alpha}) that stretches towards 900 the surface ($\alpha > 1.0$) or the bottom ($\alpha < 1.0$) and 901 user prescribed surface (\np{rn\_zs}) and bottom depths. 902 The bottom cell depth in this example is given as a function of water depth: 903 904 \[ 905 % \label{eq:DOM_zb} 906 Z_b = h a + b 907 \] 908 909 where the namelist parameters \np{rn\_zb\_a} and \np{rn\_zb\_b} are $a$ and $b$ respectively. 910 911 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 912 \begin{figure}[!ht] 913 \includegraphics[]{Fig_DOM_compare_coordinates_surface} 914 \caption{ 915 A comparison of the \citet{Song_Haidvogel_JCP94} $S$-coordinate (solid lines), 916 a 50 level $Z$-coordinate (contoured surfaces) and 917 the \citet{Siddorn_Furner_OM12} $S$-coordinate (dashed lines) in the surface $100~m$ for 918 a idealised bathymetry that goes from $50~m$ to $5500~m$ depth. 919 For clarity every third coordinate surface is shown. 920 } 921 \label{fig:fig_compare_coordinates_surface} 922 \end{figure} 923 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 924 925 This gives a smooth analytical stretching in computational space that is constrained to 926 given specified surface and bottom grid cell thicknesses in real space. 927 This is not to be confused with the hybrid schemes that 928 superimpose geopotential coordinates on terrain following coordinates thus 929 creating a non-analytical vertical coordinate that 930 therefore may suffer from large gradients in the vertical resolutions. 931 This stretching is less straightforward to implement than the \citet{Song_Haidvogel_JCP94} stretching, 932 but has the advantage of resolving diurnal processes in deep water and has generally flatter slopes. 933 934 As with the \citet{Song_Haidvogel_JCP94} stretching the stretch is only applied at depths greater than 935 the critical depth $h_c$. 936 In this example two options are available in depths shallower than $h_c$, 937 with pure sigma being applied if the \np{ln\_sigcrit} is true and pure z-coordinates if it is false 938 (the z-coordinate being equal to the depths of the stretched coordinate at $h_c$). 939 940 Minimising the horizontal slope of the vertical coordinate is important in terrain-following systems as 941 large slopes lead to hydrostatic consistency. 942 A hydrostatic consistency parameter diagnostic following \citet{Haney1991} has been implemented, 943 and is output as part of the model mesh file at the start of the run. 944 945 % ------------------------------------------------------------------------------------------------------------- 946 % z*- or s*-coordinate 947 % ------------------------------------------------------------------------------------------------------------- 948 \subsection{\zstar- or \sstar-coordinate (\protect\np{ln\_linssh}~\forcode{= .false.})} 949 \label{subsec:DOM_zgr_star} 950 951 This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO web site. 952 953 %gm% key advantage: minimise the diffusion/dispertion associated with advection in response to high frequency surface disturbances 525 526 Additionally, hybrid combinations of the three main coordinates are available: 527 $s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps}d and \autoref{fig:z_zps_s_sps}e). 528 529 A further choice related to vertical coordinate concerns the presence (or not) of ocean 530 cavities beneath ice shelves within the model domain. A setting of \np{ln\_isfcav} as 531 \forcode{.true.} indicates that the domain contains ocean cavities, otherwise the top, 532 wet layer of the ocean will always be at the ocean surface. This option is currently only 533 available for $z$- or $zps$-coordinates. In the latter case, partial steps are also applied 534 at the ocean/ice shelf interface. 535 536 Within the model, the arrays describing the grid point depths and vertical scale factors 537 are three set of three dimensional arrays $(i,j,k)$ defined at \textit{before}, 538 \textit{now} and \textit{after} time step. The time at which they are defined is 539 indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively. They are updated at each 540 model time step. The initial fixed reference coordinate system is held in variable names 541 with a $\_0$ suffix. When the linear free surface option is used 542 (\np{ln\_linssh}\forcode{ = .true.}), \textit{before}, \textit{now} and \textit{after} 543 arrays are initially set to their reference counterpart and remain fixed. 544 545 \subsubsection{Required fields} 546 \label{sec:DOM_zgr_fields} 547 The explicit specification of a range of fields related to the vertical grid are required for the definition of a configuration. These include: 548 549 \begin{Verbatim}[fontsize=\tiny] 550 int ln_zco, ln_zps, ln_sco /* flags for z-coord, z-coord with partial steps and s-coord */ 551 int ln_isfcav /* flag for ice shelf cavities */ 552 double e3t_1d, e3w_1d /* reference vertical scale factors at T and W points */ 553 double e3t_0, e3u_0, e3v_0, e3f_0, e3w_0 /* vertical scale factors 3D coordinate at T,U,V,F and W points */ 554 double e3uw_0, e3vw_0 /* vertical scale factors 3D coordinate at UW and VW points */ 555 int bottom_level, top_level /* last wet T-points, 1st wet T-points (for ice shelf cavities) */ 556 /* For reference: */ 557 float bathy_metry /* bathymetry used in setting top and bottom levels */ 558 \end{Verbatim} 559 560 This set of vertical metrics is sufficient to describe the initial depth and thickness of 561 every gridcell in the model regardless of the choice of vertical coordinate. With constant 562 z-levels, e3 metrics will be uniform across each horizontal level. In the partial step 563 case each e3 at the \np{bottom\_level} (and, possibly, \np{top\_level} if ice cavities are 564 present) may vary from its horizontal neighbours. And, in s-coordinates, variations can 565 occur throughout the water column. With the non-linear free-surface, all the coordinates 566 behave more like the s-coordinate in that variations occurr throughout the water column 567 with displacements related to the sea surface height. These variations are typically much 568 smaller than those arising from bottom fitted coordinates. The values for vertical metrics 569 supplied in the domain configuration file can be considered as those arising from a flat 570 sea surface with zero elevation. 571 572 The \np{bottom\_level} and \np{top\_level} 2D arrays define the \np{bottom\_level} and top 573 wet levels in each grid column. Without ice cavities, \np{top\_level} is essentially a land 574 mask (0 on land; 1 everywhere else). With ice cavities, \np{top\_level} determines the 575 first wet point below the overlying ice shelf. 576 577 954 578 955 579 % ------------------------------------------------------------------------------------------------------------- 956 580 % level bathymetry and mask 957 581 % ------------------------------------------------------------------------------------------------------------- 958 \subs ection{Level bathymetry and mask}582 \subsubsection{Level bathymetry and mask} 959 583 \label{subsec:DOM_msk} 960 584 961 Whatever the vertical coordinate used, the model offers the possibility of representing the bottom topography with 962 steps that follow the face of the model cells (step like topography) \citep{Madec_al_JPO96}. 963 The distribution of the steps in the horizontal is defined in a 2D integer array, mbathy, which 964 gives the number of ocean levels (\ie those that are not masked) at each $t$-point. 965 mbathy is computed from the meter bathymetry using the definiton of gdept as the number of $t$-points which 966 gdept $\leq$ bathy. 967 968 Modifications of the model bathymetry are performed in the \textit{bat\_ctl} routine (see \mdl{domzgr} module) after 969 mbathy is computed. 970 Isolated grid points that do not communicate with another ocean point at the same level are eliminated. 971 972 As for the representation of bathymetry, a 2D integer array, misfdep, is created. 973 misfdep defines the level of the first wet $t$-point. 974 All the cells between $k = 1$ and $misfdep(i,j) - 1$ are masked. 975 By default, $misfdep(:,:) = 1$ and no cells are masked. 976 977 In case of ice shelf cavities, modifications of the model bathymetry and ice shelf draft into 978 the cavities are performed in the \textit{zgr\_isf} routine. 979 The compatibility between ice shelf draft and bathymetry is checked. 980 All the locations where the isf cavity is thinnest than \np{rn\_isfhmin} meters are grounded (\ie masked). 981 If only one cell on the water column is opened at $t$-, $u$- or $v$-points, 982 the bathymetry or the ice shelf draft is dug to fit this constrain. 983 If the incompatibility is too strong (need to dig more than 1 cell), the cell is masked. 984 985 From the \textit{mbathy} and \textit{misfdep} array, the mask fields are defined as follows: 585 586 From \np{top\_level} and \np{bottom\_level} fields, the mask fields are defined as follows: 986 587 \begin{alignat*}{2} 987 588 tmask(i,j,k) &= & & 988 589 \begin{cases} 989 0 &\text{if $ k < misfdep(i,j)$} \\990 1 &\text{if $ misfdep(i,j) \leq k \leq mbathy(i,j)$} \\991 0 &\text{if $ k > mbathy(i,j)$}590 0 &\text{if $ k < top\_level(i,j)$} \\ 591 1 &\text{if $bottom\_level(i,j) \leq k \leq top\_level(i,j)$} \\ 592 0 &\text{if $ k > bottom\_level(i,j)$} 992 593 \end{cases} 993 594 \\ … … 1005 606 exactly in the same way as for the bottom boundary. 1006 607 1007 The specification of closed lateral boundaries requires that at least 1008 the first and last rows and columns of the \textit{mbathy} array are set to zero. 1009 In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to 1010 the second one and its first column equal to the last but one (and so too the mask arrays) 1011 (see \autoref{fig:LBC_jperio}). 608 %% The specification of closed lateral boundaries requires that at least 609 %% the first and last rows and columns of the \textit{mbathy} array are set to zero. 610 %% In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to 611 %% the second one and its first column equal to the last but one (and so too the mask arrays) 612 %% (see \autoref{fig:LBC_jperio}). 613 614 615 %------------------------------------------------------------------------------------------------- 616 % Closed seas 617 %------------------------------------------------------------------------------------------------- 618 \subsection{Closed seas} \label{subsec:DOM_closea} 619 620 When a global ocean is coupled to an atmospheric model it is better to represent all large 621 water bodies (\eg great lakes, Caspian sea...) even if the model resolution does not allow 622 their communication with the rest of the ocean. This is unnecessary when the ocean is 623 forced by fixed atmospheric conditions, so these seas can be removed from the ocean 624 domain. The user has the option to set the bathymetry in closed seas to zero (see 625 \autoref{sec:MISC_closea}) and to optionally decide on the fate of any freshwater 626 imbalance over the area. The options are explained in \autoref{sec:MISC_closea} but it 627 should be noted here that a successful use of these options requires appropriate mask 628 fields to be present in the domain configuration file. Among the possibilities are: 629 630 \begin{Verbatim}[fontsize=\tiny] 631 int closea_mask /* non-zero values in closed sea areas for optional masking */ 632 int closea_mask_rnf /* non-zero values in closed sea areas with runoff locations (precip only) */ 633 int closea_mask_emp /* non-zero values in closed sea areas with runoff locations (total emp) */ 634 \end{Verbatim} 635 636 % ------------------------------------------------------------------------------------------------------------- 637 % Grid files 638 % ------------------------------------------------------------------------------------------------------------- 639 \subsection{Output grid files} 640 \label{subsec:DOM_meshmask} 641 642 \nlst{namcfg} 643 644 Most of the arrays relating to a particular ocean model configuration dicussed in this 645 chapter (grid-point position, scale factors) can be saved in a file if namelist parameter 646 \np{ln\_write\_cfg} (namelist \ngn{namcfg}) is set to \forcode{.true.}; the output 647 filename is set thorugh parameter \np{cn\_domcfg\_out}. This is only really useful 648 if the fields are computed in subroutines \mdl{usrdef\_hgr} or \mdl{usrdef\_zgr} and 649 checking or confirmation is required. 650 651 \nlst{namdom} 652 653 Alternatively, all the arrays relating to a particular ocean model configuration 654 (grid-point position, scale factors, depths and masks) can be saved in a file called 655 \texttt{mesh\_mask} if namelist parameter \np{ln\_meshmask} (namelist \ngn{namdom}) is set 656 to \forcode{.true.}. This file contains additional fields that can be useful for 657 post-processing applications 1012 658 1013 659 % ================================================================ 1014 660 % Domain: Initial State (dtatsd & istate) 1015 661 % ================================================================ 1016 \section{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})} 662 \section[Initial state (\textit{istate.F90} and \textit{dtatsd.F90})] 663 {Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})} 1017 664 \label{sec:DTA_tsd} 1018 665 %-----------------------------------------namtsd------------------------------------------- 1019 1020 666 \nlst{namtsd} 1021 667 %------------------------------------------------------------------------------------------ 1022 668 1023 Options are defined in \ngn{namtsd}. 1024 By default, the ocean start from rest (the velocity field is set to zero) and the initialization of temperature and 1025 salinity fields is controlled through the \np{ln\_tsd\_ini} namelist parameter. 669 Basic initial state options are defined in \ngn{namtsd}. By default, the ocean starts 670 from rest (the velocity field is set to zero) and the initialization of temperature and 671 salinity fields is controlled through the \np{ln\_tsd\_init} namelist parameter. 672 1026 673 \begin{description} 1027 \item[\np{ln\_tsd\_init} ~\forcode{= .true.}]1028 use a T and S input files that can be given on the model grid itself or on their native input data grid.1029 In the latter case,1030 the data will be interpolated on-the-fly both in the horizontal and the vertical to the model grid1031 (see \autoref{subsec:SBC_iof}).1032 The information relative to the input files are given in the \np{sn\_tem} and \np{sn\_sal} structures.1033 The computation is done in the \mdl{dtatsd} module. 1034 \item[\np{ln\_tsd\_init}~\forcode{= .false.}] 1035 use constant salinity value of $35.5~psu$ and an analytical profile of temperature1036 (typical of the tropical ocean), see \rou{istate\_t\_s} subroutine called from \mdl{istate} module.674 \item[\np{ln\_tsd\_init}\forcode{= .true.}] 675 Use T and S input files that can be given on the model grid itself or on their native 676 input data grids. In the latter case, the data will be interpolated on-the-fly both in 677 the horizontal and the vertical to the model grid (see \autoref{subsec:SBC_iof}). The 678 information relating to the input files are specified in the \np{sn\_tem} and 679 \np{sn\_sal} structures. The computation is done in the \mdl{dtatsd} module. 680 \item[\np{ln\_tsd\_init}\forcode{= .false.}] 681 Initial values for T and S are set via a user supplied \rou{usr\_def\_istate} routine 682 contained in \mdl{userdef\_istate}. The default version sets horizontally uniform T and 683 profiles as used in the GYRE configuration (see \autoref{sec:CFG_gyre}). 1037 684 \end{description} 1038 685 -
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_DYN.tex
r10499 r11422 65 65 % Horizontal divergence and relative vorticity 66 66 %-------------------------------------------------------------------------------------------------------------- 67 \subsection{Horizontal divergence and relative vorticity (\protect\mdl{divcur})} 67 \subsection[Horizontal divergence and relative vorticity (\textit{divcur.F90})] 68 {Horizontal divergence and relative vorticity (\protect\mdl{divcur})} 68 69 \label{subsec:DYN_divcur} 69 70 … … 101 102 % Sea Surface Height evolution 102 103 %-------------------------------------------------------------------------------------------------------------- 103 \subsection{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})} 104 \subsection[Horizontal divergence and relative vorticity (\textit{sshwzv.F90})] 105 {Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})} 104 106 \label{subsec:DYN_sshwzv} 105 107 … … 127 129 Replacing $T$ by the number $1$ in the tracer equation and summing over the water column must lead to 128 130 the sea surface height equation otherwise tracer content will not be conserved 129 \citep{ Griffies_al_MWR01, Leclair_Madec_OM09}.131 \citep{griffies.pacanowski.ea_MWR01, leclair.madec_OM09}. 130 132 131 133 The vertical velocity is computed by an upward integration of the horizontal divergence starting at the bottom, … … 181 183 % Vorticity term 182 184 % ------------------------------------------------------------------------------------------------------------- 183 \subsection{Vorticity term (\protect\mdl{dynvor})} 185 \subsection[Vorticity term (\textit{dynvor.F90})] 186 {Vorticity term (\protect\mdl{dynvor})} 184 187 \label{subsec:DYN_vor} 185 188 %------------------------------------------nam_dynvor---------------------------------------------------- … … 203 206 % enstrophy conserving scheme 204 207 %------------------------------------------------------------- 205 \subsubsection{Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})} 208 \subsubsection[Enstrophy conserving scheme (\forcode{ln_dynvor_ens = .true.})] 209 {Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})} 206 210 \label{subsec:DYN_vor_ens} 207 211 … … 226 230 % energy conserving scheme 227 231 %------------------------------------------------------------- 228 \subsubsection{Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})} 232 \subsubsection[Energy conserving scheme (\forcode{ln_dynvor_ene = .true.})] 233 {Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})} 229 234 \label{subsec:DYN_vor_ene} 230 235 … … 246 251 % mix energy/enstrophy conserving scheme 247 252 %------------------------------------------------------------- 248 \subsubsection{Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{ = .true.}) } 253 \subsubsection[Mixed energy/enstrophy conserving scheme (\forcode{ln_dynvor_mix = .true.})] 254 {Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{ = .true.})} 249 255 \label{subsec:DYN_vor_mix} 250 256 … … 271 277 % energy and enstrophy conserving scheme 272 278 %------------------------------------------------------------- 273 \subsubsection{Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.}) } 279 \subsubsection[Energy and enstrophy conserving scheme (\forcode{ln_dynvor_een = .true.})] 280 {Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.})} 274 281 \label{subsec:DYN_vor_een} 275 282 … … 287 294 Nevertheless, this technique strongly distort the phase and group velocity of Rossby waves....} 288 295 289 A very nice solution to the problem of double averaging was proposed by \citet{ Arakawa_Hsu_MWR90}.296 A very nice solution to the problem of double averaging was proposed by \citet{arakawa.hsu_MWR90}. 290 297 The idea is to get rid of the double averaging by considering triad combinations of vorticity. 291 298 It is noteworthy that this solution is conceptually quite similar to the one proposed by 292 \citep{ Griffies_al_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:C}).293 294 The \citet{ Arakawa_Hsu_MWR90} vorticity advection scheme for a single layer is modified295 for spherical coordinates as described by \citet{ Arakawa_Lamb_MWR81} to obtain the EEN scheme.299 \citep{griffies.gnanadesikan.ea_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:C}). 300 301 The \citet{arakawa.hsu_MWR90} vorticity advection scheme for a single layer is modified 302 for spherical coordinates as described by \citet{arakawa.lamb_MWR81} to obtain the EEN scheme. 296 303 First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point: 297 304 \[ … … 309 316 \begin{figure}[!ht] 310 317 \begin{center} 311 \includegraphics[width= 0.70\textwidth]{Fig_DYN_een_triad}318 \includegraphics[width=\textwidth]{Fig_DYN_een_triad} 312 319 \caption{ 313 320 \protect\label{fig:DYN_een_triad} … … 327 334 (with a systematic reduction of $e_{3f}$ when a model level intercept the bathymetry) 328 335 that tends to reinforce the topostrophy of the flow 329 (\ie the tendency of the flow to follow the isobaths) \citep{ Penduff_al_OS07}.336 (\ie the tendency of the flow to follow the isobaths) \citep{penduff.le-sommer.ea_OS07}. 330 337 331 338 Next, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as … … 356 363 (\ie $\chi$=$0$) (see \autoref{subsec:C_vorEEN}). 357 364 Applied to a realistic ocean configuration, it has been shown that it leads to a significant reduction of 358 the noise in the vertical velocity field \citep{ Le_Sommer_al_OM09}.365 the noise in the vertical velocity field \citep{le-sommer.penduff.ea_OM09}. 359 366 Furthermore, used in combination with a partial steps representation of bottom topography, 360 367 it improves the interaction between current and topography, 361 leading to a larger topostrophy of the flow \citep{ Barnier_al_OD06, Penduff_al_OS07}.368 leading to a larger topostrophy of the flow \citep{barnier.madec.ea_OD06, penduff.le-sommer.ea_OS07}. 362 369 363 370 %-------------------------------------------------------------------------------------------------------------- 364 371 % Kinetic Energy Gradient term 365 372 %-------------------------------------------------------------------------------------------------------------- 366 \subsection{Kinetic energy gradient term (\protect\mdl{dynkeg})} 373 \subsection[Kinetic energy gradient term (\textit{dynkeg.F90})] 374 {Kinetic energy gradient term (\protect\mdl{dynkeg})} 367 375 \label{subsec:DYN_keg} 368 376 … … 384 392 % Vertical advection term 385 393 %-------------------------------------------------------------------------------------------------------------- 386 \subsection{Vertical advection term (\protect\mdl{dynzad}) } 394 \subsection[Vertical advection term (\textit{dynzad.F90})] 395 {Vertical advection term (\protect\mdl{dynzad})} 387 396 \label{subsec:DYN_zad} 388 397 … … 403 412 When \np{ln\_dynzad\_zts}\forcode{ = .true.}, 404 413 a split-explicit time stepping with 5 sub-timesteps is used on the vertical advection term. 405 This option can be useful when the value of the timestep is limited by vertical advection \citep{ Lemarie_OM2015}.414 This option can be useful when the value of the timestep is limited by vertical advection \citep{lemarie.debreu.ea_OM15}. 406 415 Note that in this case, 407 416 a similar split-explicit time stepping should be used on vertical advection of tracer to ensure a better stability, … … 430 439 % Coriolis plus curvature metric terms 431 440 %-------------------------------------------------------------------------------------------------------------- 432 \subsection{Coriolis plus curvature metric terms (\protect\mdl{dynvor}) } 441 \subsection[Coriolis plus curvature metric terms (\textit{dynvor.F90})] 442 {Coriolis plus curvature metric terms (\protect\mdl{dynvor})} 433 443 \label{subsec:DYN_cor_flux} 434 444 … … 451 461 % Flux form Advection term 452 462 %-------------------------------------------------------------------------------------------------------------- 453 \subsection{Flux form advection term (\protect\mdl{dynadv}) } 463 \subsection[Flux form advection term (\textit{dynadv.F90})] 464 {Flux form advection term (\protect\mdl{dynadv})} 454 465 \label{subsec:DYN_adv_flux} 455 466 … … 475 486 a $2^{nd}$ order centered finite difference scheme, CEN2, 476 487 or a $3^{rd}$ order upstream biased scheme, UBS. 477 The latter is described in \citet{ Shchepetkin_McWilliams_OM05}.488 The latter is described in \citet{shchepetkin.mcwilliams_OM05}. 478 489 The schemes are selected using the namelist logicals \np{ln\_dynadv\_cen2} and \np{ln\_dynadv\_ubs}. 479 490 In flux form, the schemes differ by the choice of a space and time interpolation to define the value of … … 484 495 % 2nd order centred scheme 485 496 %------------------------------------------------------------- 486 \subsubsection{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{ = .true.})} 497 \subsubsection[CEN2: $2^{nd}$ order centred scheme (\forcode{ln_dynadv_cen2 = .true.})] 498 {CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{ = .true.})} 487 499 \label{subsec:DYN_adv_cen2} 488 500 … … 507 519 % UBS scheme 508 520 %------------------------------------------------------------- 509 \subsubsection{UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{ = .true.})} 521 \subsubsection[UBS: Upstream Biased Scheme (\forcode{ln_dynadv_ubs = .true.})] 522 {UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{ = .true.})} 510 523 \label{subsec:DYN_adv_ubs} 511 524 … … 523 536 where $u"_{i+1/2} =\delta_{i+1/2} \left[ {\delta_i \left[ u \right]} \right]$. 524 537 This results in a dissipatively dominant (\ie hyper-diffusive) truncation error 525 \citep{ Shchepetkin_McWilliams_OM05}.526 The overall performance of the advection scheme is similar to that reported in \citet{ Farrow1995}.538 \citep{shchepetkin.mcwilliams_OM05}. 539 The overall performance of the advection scheme is similar to that reported in \citet{farrow.stevens_JPO95}. 527 540 It is a relatively good compromise between accuracy and smoothness. 528 541 It is not a \emph{positive} scheme, meaning that false extrema are permitted. … … 542 555 while the second term, which is the diffusion part of the scheme, 543 556 is evaluated using the \textit{before} velocity (forward in time). 544 This is discussed by \citet{ Webb_al_JAOT98} in the context of the Quick advection scheme.557 This is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the Quick advection scheme. 545 558 546 559 Note that the UBS and QUICK (Quadratic Upstream Interpolation for Convective Kinematics) schemes only differ by 547 560 one coefficient. 548 Replacing $1/6$ by $1/8$ in (\autoref{eq:dynadv_ubs}) leads to the QUICK advection scheme \citep{ Webb_al_JAOT98}.561 Replacing $1/6$ by $1/8$ in (\autoref{eq:dynadv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}. 549 562 This option is not available through a namelist parameter, since the $1/6$ coefficient is hard coded. 550 563 Nevertheless it is quite easy to make the substitution in the \mdl{dynadv\_ubs} module and obtain a QUICK scheme. … … 560 573 % Hydrostatic pressure gradient term 561 574 % ================================================================ 562 \section{Hydrostatic pressure gradient (\protect\mdl{dynhpg})} 575 \section[Hydrostatic pressure gradient (\textit{dynhpg.F90})] 576 {Hydrostatic pressure gradient (\protect\mdl{dynhpg})} 563 577 \label{sec:DYN_hpg} 564 578 %------------------------------------------nam_dynhpg--------------------------------------------------- … … 582 596 % z-coordinate with full step 583 597 %-------------------------------------------------------------------------------------------------------------- 584 \subsection{Full step $Z$-coordinate (\protect\np{ln\_dynhpg\_zco}\forcode{ = .true.})} 598 \subsection[Full step $Z$-coordinate (\forcode{ln_dynhpg_zco = .true.})] 599 {Full step $Z$-coordinate (\protect\np{ln\_dynhpg\_zco}\forcode{ = .true.})} 585 600 \label{subsec:DYN_hpg_zco} 586 601 … … 627 642 % z-coordinate with partial step 628 643 %-------------------------------------------------------------------------------------------------------------- 629 \subsection{Partial step $Z$-coordinate (\protect\np{ln\_dynhpg\_zps}\forcode{ = .true.})} 644 \subsection[Partial step $Z$-coordinate (\forcode{ln_dynhpg_zps = .true.})] 645 {Partial step $Z$-coordinate (\protect\np{ln\_dynhpg\_zps}\forcode{ = .true.})} 630 646 \label{subsec:DYN_hpg_zps} 631 647 … … 652 668 653 669 Pressure gradient formulations in an $s$-coordinate have been the subject of a vast number of papers 654 (\eg, \citet{ Song1998, Shchepetkin_McWilliams_OM05}).670 (\eg, \citet{song_MWR98, shchepetkin.mcwilliams_OM05}). 655 671 A number of different pressure gradient options are coded but the ROMS-like, 656 672 density Jacobian with cubic polynomial method is currently disabled whilst known bugs are under investigation. 657 673 658 $\bullet$ Traditional coding (see for example \citet{ Madec_al_JPO96}: (\np{ln\_dynhpg\_sco}\forcode{ = .true.})674 $\bullet$ Traditional coding (see for example \citet{madec.delecluse.ea_JPO96}: (\np{ln\_dynhpg\_sco}\forcode{ = .true.}) 659 675 \begin{equation} 660 676 \label{eq:dynhpg_sco} … … 679 695 $\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np{ln\_dynhpg\_prj}\forcode{ = .true.}) 680 696 681 $\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{ Shchepetkin_McWilliams_OM05}697 $\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{shchepetkin.mcwilliams_OM05} 682 698 (\np{ln\_dynhpg\_djc}\forcode{ = .true.}) (currently disabled; under development) 683 699 684 700 Note that expression \autoref{eq:dynhpg_sco} is commonly used when the variable volume formulation is activated 685 701 (\key{vvl}) because in that case, even with a flat bottom, 686 the coordinate surfaces are not horizontal but follow the free surface \citep{ Levier2007}.702 the coordinate surfaces are not horizontal but follow the free surface \citep{levier.treguier.ea_rpt07}. 687 703 The pressure jacobian scheme (\np{ln\_dynhpg\_prj}\forcode{ = .true.}) is available as 688 704 an improved option to \np{ln\_dynhpg\_sco}\forcode{ = .true.} when \key{vvl} is active. … … 704 720 corresponds to the water replaced by the ice shelf. 705 721 This top pressure is constant over time. 706 A detailed description of this method is described in \citet{ Losch2008}.\\722 A detailed description of this method is described in \citet{losch_JGR08}.\\ 707 723 708 724 The pressure gradient due to ocean load is computed using the expression \autoref{eq:dynhpg_sco} described in … … 712 728 % Time-scheme 713 729 %-------------------------------------------------------------------------------------------------------------- 714 \subsection{Time-scheme (\protect\np{ln\_dynhpg\_imp}\forcode{ = .true./.false.})} 730 \subsection[Time-scheme (\forcode{ln_dynhpg_imp = .{true,false}.})] 731 {Time-scheme (\protect\np{ln\_dynhpg\_imp}\forcode{ = .\{true,false\}}.)} 715 732 \label{subsec:DYN_hpg_imp} 716 733 … … 722 739 the physical phenomenon that controls the time-step is internal gravity waves (IGWs). 723 740 A semi-implicit scheme for doubling the stability limit associated with IGWs can be used 724 \citep{ Brown_Campana_MWR78, Maltrud1998}.741 \citep{brown.campana_MWR78, maltrud.smith.ea_JGR98}. 725 742 It involves the evaluation of the hydrostatic pressure gradient as 726 743 an average over the three time levels $t-\rdt$, $t$, and $t+\rdt$ … … 773 790 % Surface Pressure Gradient 774 791 % ================================================================ 775 \section{Surface pressure gradient (\protect\mdl{dynspg})} 792 \section[Surface pressure gradient (\textit{dynspg.F90})] 793 {Surface pressure gradient (\protect\mdl{dynspg})} 776 794 \label{sec:DYN_spg} 777 795 %-----------------------------------------nam_dynspg---------------------------------------------------- … … 790 808 which imposes a very small time step when an explicit time stepping is used. 791 809 Two methods are proposed to allow a longer time step for the three-dimensional equations: 792 the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:PE_flt }),810 the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:PE_flt?}), 793 811 and the split-explicit free surface described below. 794 812 The extra term introduced in the filtered method is calculated implicitly, … … 811 829 % Explicit free surface formulation 812 830 %-------------------------------------------------------------------------------------------------------------- 813 \subsection{Explicit free surface (\protect\key{dynspg\_exp})} 831 \subsection[Explicit free surface (\texttt{\textbf{key\_dynspg\_exp}})] 832 {Explicit free surface (\protect\key{dynspg\_exp})} 814 833 \label{subsec:DYN_spg_exp} 815 834 … … 837 856 % Split-explict free surface formulation 838 857 %-------------------------------------------------------------------------------------------------------------- 839 \subsection{Split-explicit free surface (\protect\key{dynspg\_ts})} 858 \subsection[Split-explicit free surface (\texttt{\textbf{key\_dynspg\_ts}})] 859 {Split-explicit free surface (\protect\key{dynspg\_ts})} 840 860 \label{subsec:DYN_spg_ts} 841 861 %------------------------------------------namsplit----------------------------------------------------------- … … 845 865 846 866 The split-explicit free surface formulation used in \NEMO (\key{dynspg\_ts} defined), 847 also called the time-splitting formulation, follows the one proposed by \citet{ Shchepetkin_McWilliams_OM05}.867 also called the time-splitting formulation, follows the one proposed by \citet{shchepetkin.mcwilliams_OM05}. 848 868 The general idea is to solve the free surface equation and the associated barotropic velocity equations with 849 869 a smaller time step than $\rdt$, the time step used for the three dimensional prognostic variables … … 862 882 \begin{equation} 863 883 \label{eq:BT_dyn} 864 \frac{\partial {\ rm \overline{{\bf U}}_h} }{\partial t}=865 -f\;{\ rm {\bf k}}\times {\rm \overline{{\bf U}}_h}866 -g\nabla _h \eta -\frac{c_b^{\textbf U}}{H+\eta} \ rm {\overline{{\bf U}}_h} + \rm {\overline{\bf G}}884 \frac{\partial {\mathrm \overline{{\mathbf U}}_h} }{\partial t}= 885 -f\;{\mathrm {\mathbf k}}\times {\mathrm \overline{{\mathbf U}}_h} 886 -g\nabla _h \eta -\frac{c_b^{\textbf U}}{H+\eta} \mathrm {\overline{{\mathbf U}}_h} + \mathrm {\overline{\mathbf G}} 867 887 \end{equation} 868 888 \[ 869 889 % \label{eq:BT_ssh} 870 \frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\ rm{\bf \overline{U}}}_h \,} \right]+P-E890 \frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\mathrm{\mathbf \overline{U}}}_h \,} \right]+P-E 871 891 \] 872 892 % \end{subequations} 873 where $\ rm {\overline{\bf G}}$ is a forcing term held constant, containing coupling term between modes,893 where $\mathrm {\overline{\mathbf G}}$ is a forcing term held constant, containing coupling term between modes, 874 894 surface atmospheric forcing as well as slowly varying barotropic terms not explicitly computed to gain efficiency. 875 895 The third term on the right hand side of \autoref{eq:BT_dyn} represents the bottom stress 876 896 (see section \autoref{sec:ZDF_bfr}), explicitly accounted for at each barotropic iteration. 877 897 Temporal discretization of the system above follows a three-time step Generalized Forward Backward algorithm 878 detailed in \citet{ Shchepetkin_McWilliams_OM05}.898 detailed in \citet{shchepetkin.mcwilliams_OM05}. 879 899 AB3-AM4 coefficients used in \NEMO follow the second-order accurate, 880 "multi-purpose" stability compromise as defined in \citet{ Shchepetkin_McWilliams_Bk08}900 "multi-purpose" stability compromise as defined in \citet{shchepetkin.mcwilliams_ibk09} 881 901 (see their figure 12, lower left). 882 902 … … 884 904 \begin{figure}[!t] 885 905 \begin{center} 886 \includegraphics[width= 0.7\textwidth]{Fig_DYN_dynspg_ts}906 \includegraphics[width=\textwidth]{Fig_DYN_dynspg_ts} 887 907 \caption{ 888 908 \protect\label{fig:DYN_dynspg_ts} … … 936 956 and time splitting not compatible. 937 957 Advective barotropic velocities are obtained by using a secondary set of filtering weights, 938 uniquely defined from the filter coefficients used for the time averaging (\citet{ Shchepetkin_McWilliams_OM05}).958 uniquely defined from the filter coefficients used for the time averaging (\citet{shchepetkin.mcwilliams_OM05}). 939 959 Consistency between the time averaged continuity equation and the time stepping of tracers is here the key to 940 960 obtain exact conservation. … … 953 973 external gravity waves in idealized or weakly non-linear cases. 954 974 Although the damping is lower than for the filtered free surface, 955 it is still significant as shown by \citet{ Levier2007} in the case of an analytical barotropic Kelvin wave.975 it is still significant as shown by \citet{levier.treguier.ea_rpt07} in the case of an analytical barotropic Kelvin wave. 956 976 957 977 %>>>>>=============== … … 1051 1071 the leap-frog splitting mode in equation \autoref{eq:DYN_spg_ts_ssh}. 1052 1072 We have tried various forms of such filtering, 1053 with the following method discussed in \cite{ Griffies_al_MWR01} chosen due to1073 with the following method discussed in \cite{griffies.pacanowski.ea_MWR01} chosen due to 1054 1074 its stability and reasonably good maintenance of tracer conservation properties (see ??). 1055 1075 … … 1081 1101 % Filtered free surface formulation 1082 1102 %-------------------------------------------------------------------------------------------------------------- 1083 \subsection{Filtered free surface (\protect\key{dynspg\_flt})} 1103 \subsection[Filtered free surface (\texttt{\textbf{key\_dynspg\_flt}})] 1104 {Filtered free surface (\protect\key{dynspg\_flt})} 1084 1105 \label{subsec:DYN_spg_fltp} 1085 1106 1086 The filtered formulation follows the \citet{ Roullet_Madec_JGR00} implementation.1107 The filtered formulation follows the \citet{roullet.madec_JGR00} implementation. 1087 1108 The extra term introduced in the equations (see \autoref{subsec:PE_free_surface}) is solved implicitly. 1088 1109 The elliptic solvers available in the code are documented in \autoref{chap:MISC}. … … 1092 1113 \[ 1093 1114 % \label{eq:spg_flt} 1094 \frac{\partial {\ rm {\bf U}}_h }{\partial t}= {\rm {\bf M}}1115 \frac{\partial {\mathrm {\mathbf U}}_h }{\partial t}= {\mathrm {\mathbf M}} 1095 1116 - g \nabla \left( \tilde{\rho} \ \eta \right) 1096 1117 - g \ T_c \nabla \left( \widetilde{\rho} \ \partial_t \eta \right) … … 1098 1119 where $T_c$, is a parameter with dimensions of time which characterizes the force, 1099 1120 $\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, 1100 and $\ rm {\bf M}$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient,1121 and $\mathrm {\mathbf M}$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient, 1101 1122 non-linear and viscous terms in \autoref{eq:PE_dyn}. 1102 1123 } %end gmcomment … … 1109 1130 % Lateral diffusion term 1110 1131 % ================================================================ 1111 \section{Lateral diffusion term and operators (\protect\mdl{dynldf})} 1132 \section[Lateral diffusion term and operators (\textit{dynldf.F90})] 1133 {Lateral diffusion term and operators (\protect\mdl{dynldf})} 1112 1134 \label{sec:DYN_ldf} 1113 1135 %------------------------------------------nam_dynldf---------------------------------------------------- … … 1143 1165 1144 1166 % ================================================================ 1145 \subsection[Iso-level laplacian (\ protect\np{ln\_dynldf\_lap}\forcode{= .true.})]1146 1167 \subsection[Iso-level laplacian (\forcode{ln_dynldf_lap = .true.})] 1168 {Iso-level laplacian operator (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})} 1147 1169 \label{subsec:DYN_ldf_lap} 1148 1170 … … 1152 1174 \left\{ 1153 1175 \begin{aligned} 1154 D_u^{l{\ rm {\bf U}}} =\frac{1}{e_{1u} }\delta_{i+1/2} \left[ {A_T^{lm}1176 D_u^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{1u} }\delta_{i+1/2} \left[ {A_T^{lm} 1155 1177 \;\chi } \right]-\frac{1}{e_{2u} {\kern 1pt}e_{3u} }\delta_j \left[ 1156 1178 {A_f^{lm} \;e_{3f} \zeta } \right] \\ \\ 1157 D_v^{l{\ rm {\bf U}}} =\frac{1}{e_{2v} }\delta_{j+1/2} \left[ {A_T^{lm}1179 D_v^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{2v} }\delta_{j+1/2} \left[ {A_T^{lm} 1158 1180 \;\chi } \right]+\frac{1}{e_{1v} {\kern 1pt}e_{3v} }\delta_i \left[ 1159 1181 {A_f^{lm} \;e_{3f} \zeta } \right] … … 1169 1191 % Rotated laplacian operator 1170 1192 %-------------------------------------------------------------------------------------------------------------- 1171 \subsection[Rotated laplacian (\ protect\np{ln\_dynldf\_iso}\forcode{= .true.})]1172 1193 \subsection[Rotated laplacian (\forcode{ln_dynldf_iso = .true.})] 1194 {Rotated laplacian operator (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})} 1173 1195 \label{subsec:DYN_ldf_iso} 1174 1196 … … 1228 1250 % Iso-level bilaplacian operator 1229 1251 %-------------------------------------------------------------------------------------------------------------- 1230 \subsection[Iso-level bilaplacian (\ protect\np{ln\_dynldf\_bilap}\forcode{= .true.})]1231 1252 \subsection[Iso-level bilaplacian (\forcode{ln_dynldf_bilap = .true.})] 1253 {Iso-level bilaplacian operator (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})} 1232 1254 \label{subsec:DYN_ldf_bilap} 1233 1255 … … 1243 1265 % Vertical diffusion term 1244 1266 % ================================================================ 1245 \section{Vertical diffusion term (\protect\mdl{dynzdf})} 1267 \section[Vertical diffusion term (\textit{dynzdf.F90})] 1268 {Vertical diffusion term (\protect\mdl{dynzdf})} 1246 1269 \label{sec:DYN_zdf} 1247 1270 %----------------------------------------------namzdf------------------------------------------------------ … … 1326 1349 There are two main options for wetting and drying code (wd): 1327 1350 (a) an iterative limiter (il) and (b) a directional limiter (dl). 1328 The directional limiter is based on the scheme developed by \cite{ WarnerEtal13} for RO1351 The directional limiter is based on the scheme developed by \cite{warner.defne.ea_CG13} for RO 1329 1352 MS 1330 which was in turn based on ideas developed for POM by \cite{ Oey06}. The iterative1353 which was in turn based on ideas developed for POM by \cite{oey_OM06}. The iterative 1331 1354 limiter is a new scheme. The iterative limiter is activated by setting $\mathrm{ln\_wd\_il} = \mathrm{.true.}$ 1332 1355 and $\mathrm{ln\_wd\_dl} = \mathrm{.false.}$. The directional limiter is activated … … 1372 1395 % Iterative limiters 1373 1396 %----------------------------------------------------------------------------------------- 1374 \subsection [Directional limiter (\textit{wet\_dry})]1375 1397 \subsection[Directional limiter (\textit{wet\_dry.F90})] 1398 {Directional limiter (\mdl{wet\_dry})} 1376 1399 \label{subsec:DYN_wd_directional_limiter} 1377 1400 The principal idea of the directional limiter is that … … 1400 1423 1401 1424 1402 \cite{ WarnerEtal13} state that in their scheme the velocity masks at the cell faces for the baroclinic1425 \cite{warner.defne.ea_CG13} state that in their scheme the velocity masks at the cell faces for the baroclinic 1403 1426 timesteps are set to 0 or 1 depending on whether the average of the masks over the barotropic sub-steps is respectively less than 1404 1427 or greater than 0.5. That scheme does not conserve tracers in integrations started from constant tracer … … 1412 1435 %----------------------------------------------------------------------------------------- 1413 1436 1414 \subsection [Iterative limiter (\textit{wet\_dry})]1415 1437 \subsection[Iterative limiter (\textit{wet\_dry.F90})] 1438 {Iterative limiter (\mdl{wet\_dry})} 1416 1439 \label{subsec:DYN_wd_iterative_limiter} 1417 1440 1418 \subsubsection [Iterative flux limiter (\textit{wet\_dry})]1419 1441 \subsubsection[Iterative flux limiter (\textit{wet\_dry.F90})] 1442 {Iterative flux limiter (\mdl{wet\_dry})} 1420 1443 \label{subsubsec:DYN_wd_il_spg_limiter} 1421 1444 … … 1494 1517 \end{equation} 1495 1518 1496 Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\it [Q: Why is1519 Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\itshape [Q: Why is 1497 1520 this necessary/desirable?]}. Substituting from (\ref{dyn_wd_continuity_coef}) gives an 1498 1521 expression for the coefficient needed to multiply the outward flux at this cell in order … … 1522 1545 % Surface pressure gradients 1523 1546 %---------------------------------------------------------------------------------------- 1524 \subsubsection [Modification of surface pressure gradients (\textit{dynhpg})]1525 1547 \subsubsection[Modification of surface pressure gradients (\textit{dynhpg.F90})] 1548 {Modification of surface pressure gradients (\mdl{dynhpg})} 1526 1549 \label{subsubsec:DYN_wd_il_spg} 1527 1550 … … 1541 1564 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1542 1565 \begin{figure}[!ht] \begin{center} 1543 \includegraphics[width= 0.8\textwidth]{Fig_WAD_dynhpg}1566 \includegraphics[width=\textwidth]{Fig_WAD_dynhpg} 1544 1567 \caption{ \label{Fig_WAD_dynhpg} 1545 1568 Illustrations of the three possible combinations of the logical variables controlling the … … 1588 1611 conditions. 1589 1612 1590 \subsubsection [Additional considerations (\textit{usrdef\_zgr})]1591 1613 \subsubsection[Additional considerations (\textit{usrdef\_zgr.F90})] 1614 {Additional considerations (\mdl{usrdef\_zgr})} 1592 1615 \label{subsubsec:WAD_additional} 1593 1616 … … 1603 1626 % The WAD test cases 1604 1627 %---------------------------------------------------------------------------------------- 1605 \subsection [The WAD test cases (\textit{usrdef\_zgr})]1606 1628 \subsection[The WAD test cases (\textit{usrdef\_zgr.F90})] 1629 {The WAD test cases (\mdl{usrdef\_zgr})} 1607 1630 \label{WAD_test_cases} 1608 1631 … … 1614 1637 % Time evolution term 1615 1638 % ================================================================ 1616 \section{Time evolution term (\protect\mdl{dynnxt})} 1639 \section[Time evolution term (\textit{dynnxt.F90})] 1640 {Time evolution term (\protect\mdl{dynnxt})} 1617 1641 \label{sec:DYN_nxt} 1618 1642 -
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_LBC.tex
r10614 r11422 17 17 % Boundary Condition at the Coast 18 18 % ================================================================ 19 \section{Boundary condition at the coast (\protect\np{rn\_shlat})} 19 \section[Boundary condition at the coast (\texttt{rn\_shlat})] 20 {Boundary condition at the coast (\protect\np{rn\_shlat})} 20 21 \label{sec:LBC_coast} 21 22 %--------------------------------------------nam_lbc------------------------------------------------------- … … 56 57 \begin{figure}[!t] 57 58 \begin{center} 58 \includegraphics[width= 0.90\textwidth]{Fig_LBC_uv}59 \includegraphics[width=\textwidth]{Fig_LBC_uv} 59 60 \caption{ 60 61 \protect\label{fig:LBC_uv} … … 85 86 \begin{figure}[!p] 86 87 \begin{center} 87 \includegraphics[width= 0.90\textwidth]{Fig_LBC_shlat}88 \includegraphics[width=\textwidth]{Fig_LBC_shlat} 88 89 \caption{ 89 90 \protect\label{fig:LBC_shlat} … … 147 148 % Boundary Condition around the Model Domain 148 149 % ================================================================ 149 \section{Model domain boundary condition (\protect\np{jperio})} 150 \section[Model domain boundary condition (\texttt{jperio})] 151 {Model domain boundary condition (\protect\np{jperio})} 150 152 \label{sec:LBC_jperio} 151 153 … … 158 160 % Closed, cyclic (\np{jperio}\forcode{ = 0..2}) 159 161 % ------------------------------------------------------------------------------------------------------------- 160 \subsection{Closed, cyclic (\protect\np{jperio}\forcode{= [0127]})} 162 \subsection[Closed, cyclic (\forcode{jperio = [0127]})] 163 {Closed, cyclic (\protect\np{jperio}\forcode{ = [0127]})} 161 164 \label{subsec:LBC_jperio012} 162 165 … … 194 197 \begin{figure}[!t] 195 198 \begin{center} 196 \includegraphics[width= 1.0\textwidth]{Fig_LBC_jperio}199 \includegraphics[width=\textwidth]{Fig_LBC_jperio} 197 200 \caption{ 198 201 \protect\label{fig:LBC_jperio} … … 206 209 % North fold (\textit{jperio = 3 }to $6)$ 207 210 % ------------------------------------------------------------------------------------------------------------- 208 \subsection{North-fold (\protect\np{jperio}\forcode{ = 3..6})} 211 \subsection[North-fold (\forcode{jperio = [3-6]})] 212 {North-fold (\protect\np{jperio}\forcode{ = [3-6]})} 209 213 \label{subsec:LBC_north_fold} 210 214 … … 218 222 \begin{figure}[!t] 219 223 \begin{center} 220 \includegraphics[width= 0.90\textwidth]{Fig_North_Fold_T}224 \includegraphics[width=\textwidth]{Fig_North_Fold_T} 221 225 \caption{ 222 226 \protect\label{fig:North_Fold_T} … … 232 236 % Exchange with neighbouring processors 233 237 % ==================================================================== 234 \section{Exchange with neighbouring processors (\protect\mdl{lbclnk}, \protect\mdl{lib\_mpp})} 238 \section[Exchange with neighbouring processors (\textit{lbclnk.F90}, \textit{lib\_mpp.F90})] 239 {Exchange with neighbouring processors (\protect\mdl{lbclnk}, \protect\mdl{lib\_mpp})} 235 240 \label{sec:LBC_mpp} 236 241 … … 280 285 \begin{figure}[!t] 281 286 \begin{center} 282 \includegraphics[width= 0.90\textwidth]{Fig_mpp}287 \includegraphics[width=\textwidth]{Fig_mpp} 283 288 \caption{ 284 289 \protect\label{fig:mpp} … … 360 365 \begin{figure}[!ht] 361 366 \begin{center} 362 \includegraphics[width= 0.90\textwidth]{Fig_mppini2}367 \includegraphics[width=\textwidth]{Fig_mppini2} 363 368 \caption { 364 369 \protect\label{fig:mppini2} … … 395 400 396 401 The BDY module was modelled on the OBC module (see NEMO 3.4) and shares many features and 397 a similar coding structure \citep{ Chanut2005}.402 a similar coding structure \citep{chanut_rpt05}. 398 403 The specification of the location of the open boundary is completely flexible and 399 404 allows for example the open boundary to follow an isobath or other irregular contour. … … 475 480 \label{subsec:BDY_FRS_scheme} 476 481 477 The Flow Relaxation Scheme (FRS) \citep{ Davies_QJRMS76,Engerdahl_Tel95},482 The Flow Relaxation Scheme (FRS) \citep{davies_QJRMS76,engedahl_T95}, 478 483 applies a simple relaxation of the model fields to externally-specified values over 479 484 a zone next to the edge of the model domain. … … 514 519 \label{subsec:BDY_flather_scheme} 515 520 516 The \citet{ Flather_JPO94} scheme is a radiation condition on the normal,521 The \citet{flather_JPO94} scheme is a radiation condition on the normal, 517 522 depth-mean transport across the open boundary. 518 523 It takes the form … … 535 540 \label{subsec:BDY_orlanski_scheme} 536 541 537 The Orlanski scheme is based on the algorithm described by \citep{ Marchesiello2001}, hereafter MMS.542 The Orlanski scheme is based on the algorithm described by \citep{marchesiello.mcwilliams.ea_OM01}, hereafter MMS. 538 543 539 544 The adaptive Orlanski condition solves a wave plus relaxation equation at the boundary: … … 636 641 \begin{figure}[!t] 637 642 \begin{center} 638 \includegraphics[width= 1.0\textwidth]{Fig_LBC_bdy_geom}643 \includegraphics[width=\textwidth]{Fig_LBC_bdy_geom} 639 644 \caption { 640 645 \protect\label{fig:LBC_bdy_geom} … … 670 675 These restrictions mean that data files used with versions of the 671 676 model prior to Version 3.4 may not work with Version 3.4 onwards. 672 A \fortran utility {\it bdy\_reorder} exists in the TOOLS directory which677 A \fortran utility {\itshape bdy\_reorder} exists in the TOOLS directory which 673 678 will re-order the data in old BDY data files. 674 679 … … 676 681 \begin{figure}[!t] 677 682 \begin{center} 678 \includegraphics[width= 1.0\textwidth]{Fig_LBC_nc_header}683 \includegraphics[width=\textwidth]{Fig_LBC_nc_header} 679 684 \caption { 680 685 \protect\label{fig:LBC_nc_header} -
NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_LDF.tex
r10442 r11422 23 23 These three aspects of the lateral diffusion are set through namelist parameters 24 24 (see the \ngn{nam\_traldf} and \ngn{nam\_dynldf} below). 25 Note that this chapter describes the standard implementation of iso-neutral tracer mixing ,26 and Griffies's implementation, which is used if \np{traldf\_grif}\forcode{ = .true.},27 is described in Appdx\autoref{apdx:triad}28 29 %-----------------------------------nam _traldf - nam_dynldf--------------------------------------------25 Note that this chapter describes the standard implementation of iso-neutral tracer mixing. 26 Griffies's implementation, which is used if \np{ln\_traldf\_triad}\forcode{ = .true.}, 27 is described in \autoref{apdx:triad} 28 29 %-----------------------------------namtra_ldf - namdyn_ldf-------------------------------------------- 30 30 31 31 \nlst{namtra_ldf} … … 34 34 %-------------------------------------------------------------------------------------------------------------- 35 35 36 % ================================================================ 37 % Lateral Mixing Operator 38 % ================================================================ 39 \section[Lateral mixing operators] 40 {Lateral mixing operators} 41 \label{sec:LDF_op} 42 We remind here the different lateral mixing operators that can be used. Further details can be found in \autoref{subsec:TRA_ldf_op} and \autoref{sec:DYN_ldf}. 43 44 \subsection[No lateral mixing (\forcode{ln_traldf_OFF}, \forcode{ln_dynldf_OFF})] 45 {No lateral mixing (\protect\np{ln\_traldf\_OFF}, \protect\np{ln\_dynldf\_OFF})} 46 47 It is possible to run without explicit lateral diffusion on tracers (\protect\np{ln\_traldf\_OFF}\forcode{ = .true.}) and/or 48 momentum (\protect\np{ln\_dynldf\_OFF}\forcode{ = .true.}). The latter option is even recommended if using the 49 UBS advection scheme on momentum (\np{ln\_dynadv\_ubs}\forcode{ = .true.}, 50 see \autoref{subsec:DYN_adv_ubs}) and can be useful for testing purposes. 51 52 \subsection[Laplacian mixing (\forcode{ln_traldf_lap}, \forcode{ln_dynldf_lap})] 53 {Laplacian mixing (\protect\np{ln\_traldf\_lap}, \protect\np{ln\_dynldf\_lap})} 54 Setting \protect\np{ln\_traldf\_lap}\forcode{ = .true.} and/or \protect\np{ln\_dynldf\_lap}\forcode{ = .true.} enables 55 a second order diffusion on tracers and momentum respectively. Note that in \NEMO 4, one can not combine 56 Laplacian and Bilaplacian operators for the same variable. 57 58 \subsection[Bilaplacian mixing (\forcode{ln_traldf_blp}, \forcode{ln_dynldf_blp})] 59 {Bilaplacian mixing (\protect\np{ln\_traldf\_blp}, \protect\np{ln\_dynldf\_blp})} 60 Setting \protect\np{ln\_traldf\_blp}\forcode{ = .true.} and/or \protect\np{ln\_dynldf\_blp}\forcode{ = .true.} enables 61 a fourth order diffusion on tracers and momentum respectively. It is implemented by calling the above Laplacian operator twice. 62 We stress again that from \NEMO 4, the simultaneous use Laplacian and Bilaplacian operators is not allowed. 36 63 37 64 % ================================================================ 38 65 % Direction of lateral Mixing 39 66 % ================================================================ 40 \section{Direction of lateral mixing (\protect\mdl{ldfslp})} 67 \section[Direction of lateral mixing (\textit{ldfslp.F90})] 68 {Direction of lateral mixing (\protect\mdl{ldfslp})} 41 69 \label{sec:LDF_slp} 42 70 … … 44 72 \gmcomment{ 45 73 we should emphasize here that the implementation is a rather old one. 46 Better work can be achieved by using \citet{ Griffies_al_JPO98, Griffies_Bk04} iso-neutral scheme.74 Better work can be achieved by using \citet{griffies.gnanadesikan.ea_JPO98, griffies_bk04} iso-neutral scheme. 47 75 } 48 76 … … 87 115 %gm% caution I'm not sure the simplification was a good idea! 88 116 89 These slopes are computed once in \rou{ldf slp\_init} when \np{ln\_sco}\forcode{ = .true.}rue,117 These slopes are computed once in \rou{ldf\_slp\_init} when \np{ln\_sco}\forcode{ = .true.}, 90 118 and either \np{ln\_traldf\_hor}\forcode{ = .true.} or \np{ln\_dynldf\_hor}\forcode{ = .true.}. 91 119 … … 119 147 %In practice, \autoref{eq:ldfslp_iso} is of little help in evaluating the neutral surface slopes. Indeed, for an unsimplified equation of state, the density has a strong dependancy on pressure (here approximated as the depth), therefore applying \autoref{eq:ldfslp_iso} using the $in situ$ density, $\rho$, computed at T-points leads to a flattening of slopes as the depth increases. This is due to the strong increase of the $in situ$ density with depth. 120 148 121 %By definition, neutral surfaces are tangent to the local $in situ$ density \citep{ McDougall1987}, therefore in \autoref{eq:ldfslp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters).149 %By definition, neutral surfaces are tangent to the local $in situ$ density \citep{mcdougall_JPO87}, therefore in \autoref{eq:ldfslp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters). 122 150 123 151 %In the $z$-coordinate, the derivative of the \autoref{eq:ldfslp_iso} numerator is evaluated at the same depth \nocite{as what?} ($T$-level, which is the same as the $u$- and $v$-levels), so the $in situ$ density can be used for its evaluation. … … 135 163 thus the $in situ$ density can be used. 136 164 This is not the case for the vertical derivatives: $\delta_{k+1/2}[\rho]$ is replaced by $-\rho N^2/g$, 137 where $N^2$ is the local Brunt-Vais\"{a}l\"{a} frequency evaluated following \citet{ McDougall1987}165 where $N^2$ is the local Brunt-Vais\"{a}l\"{a} frequency evaluated following \citet{mcdougall_JPO87} 138 166 (see \autoref{subsec:TRA_bn2}). 139 167 … … 144 172 \item[$s$- or hybrid $s$-$z$- coordinate: ] 145 173 in the current release of \NEMO, iso-neutral mixing is only employed for $s$-coordinates if 146 the Griffies scheme is used (\np{ traldf\_grif}\forcode{ = .true.};147 see Appdx\autoref{apdx:triad}).174 the Griffies scheme is used (\np{ln\_traldf\_triad}\forcode{ = .true.}; 175 see \autoref{apdx:triad}). 148 176 In other words, iso-neutral mixing will only be accurately represented with a linear equation of state 149 (\np{ nn\_eos}\forcode{ = 1..2}).177 (\np{ln\_seos}\forcode{ = .true.}). 150 178 In the case of a "true" equation of state, the evaluation of $i$ and $j$ derivatives in \autoref{eq:ldfslp_iso} 151 179 will include a pressure dependent part, leading to the wrong evaluation of the neutral slopes. … … 154 182 Note: The solution for $s$-coordinate passes trough the use of different (and better) expression for 155 183 the constraint on iso-neutral fluxes. 156 Following \citet{ Griffies_Bk04}, instead of specifying directly that there is a zero neutral diffusive flux of184 Following \citet{griffies_bk04}, instead of specifying directly that there is a zero neutral diffusive flux of 157 185 locally referenced potential density, we stay in the $T$-$S$ plane and consider the balance between 158 186 the neutral direction diffusive fluxes of potential temperature and salinity: … … 197 225 198 226 This implementation is a rather old one. 199 It is similar to the one proposed by Cox [1987], except for the background horizontal diffusion.200 Indeed, the Coximplementation of isopycnal diffusion in GFDL-type models requires227 It is similar to the one proposed by \citet{cox_OM87}, except for the background horizontal diffusion. 228 Indeed, the \citet{cox_OM87} implementation of isopycnal diffusion in GFDL-type models requires 201 229 a minimum background horizontal diffusion for numerical stability reasons. 202 230 To overcome this problem, several techniques have been proposed in which the numerical schemes of 203 the ocean model are modified \citep{ Weaver_Eby_JPO97, Griffies_al_JPO98}.204 Griffies's scheme is now available in \NEMO if \np{ traldf\_grif\_iso} is set true; see Appdx\autoref{apdx:triad}.205 Here, another strategy is presented \citep{ Lazar_PhD97}:231 the ocean model are modified \citep{weaver.eby_JPO97, griffies.gnanadesikan.ea_JPO98}. 232 Griffies's scheme is now available in \NEMO if \np{ln\_traldf\_triad}=\forcode{= .true.}; see \autoref{apdx:triad}. 233 Here, another strategy is presented \citep{lazar_phd97}: 206 234 a local filtering of the iso-neutral slopes (made on 9 grid-points) prevents the development of 207 235 grid point noise generated by the iso-neutral diffusion operator (\autoref{fig:LDF_ZDF1}). … … 212 240 213 241 Nevertheless, this iso-neutral operator does not ensure that variance cannot increase, 214 contrary to the \citet{ Griffies_al_JPO98} operator which has that property.242 contrary to the \citet{griffies.gnanadesikan.ea_JPO98} operator which has that property. 215 243 216 244 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 217 245 \begin{figure}[!ht] 218 246 \begin{center} 219 \includegraphics[width= 0.70\textwidth]{Fig_LDF_ZDF1}247 \includegraphics[width=\textwidth]{Fig_LDF_ZDF1} 220 248 \caption { 221 249 \protect\label{fig:LDF_ZDF1} … … 235 263 236 264 237 % In addition and also for numerical stability reasons \citep{ Cox1987, Griffies_Bk04},265 % In addition and also for numerical stability reasons \citep{cox_OM87, griffies_bk04}, 238 266 % the slopes are bounded by $1/100$ everywhere. This limit is decreasing linearly 239 267 % to zero fom $70$ meters depth and the surface (the fact that the eddies "feel" the 240 268 % surface motivates this flattening of isopycnals near the surface). 241 269 242 For numerical stability reasons \citep{ Cox1987, Griffies_Bk04}, the slopes must also be bounded by243 $1/100$everywhere.270 For numerical stability reasons \citep{cox_OM87, griffies_bk04}, the slopes must also be bounded by 271 the namelist scalar \np{rn\_slpmax} (usually $1/100$) everywhere. 244 272 This constraint is applied in a piecewise linear fashion, increasing from zero at the surface to 245 273 $1/100$ at $70$ metres and thereafter decreasing to zero at the bottom of the ocean 246 274 (the fact that the eddies "feel" the surface motivates this flattening of isopycnals near the surface). 275 \colorbox{yellow}{The way slopes are tapered has be checked. Not sure that this is still what is actually done.} 247 276 248 277 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 249 278 \begin{figure}[!ht] 250 279 \begin{center} 251 \includegraphics[width= 0.70\textwidth]{Fig_eiv_slp}280 \includegraphics[width=\textwidth]{Fig_eiv_slp} 252 281 \caption{ 253 282 \protect\label{fig:eiv_slp} … … 297 326 (see \autoref{sec:LBC_coast}). 298 327 299 300 % ================================================================301 % Lateral Mixing Operator302 % ================================================================303 \section{Lateral mixing operators (\protect\mdl{traldf}, \protect\mdl{traldf}) }304 \label{sec:LDF_op}305 306 307 328 308 329 % ================================================================ 309 330 % Lateral Mixing Coefficients 310 331 % ================================================================ 311 \section{Lateral mixing coefficient (\protect\mdl{ldftra}, \protect\mdl{ldfdyn}) } 332 \section[Lateral mixing coefficient (\forcode{nn_aht_ijk_t}, \forcode{nn_ahm_ijk_t})] 333 {Lateral mixing coefficient (\protect\np{nn\_aht\_ijk\_t}, \protect\np{nn\_ahm\_ijk\_t})} 312 334 \label{sec:LDF_coef} 313 335 314 Introducing a space variation in the lateral eddy mixing coefficients changes the model core memory requirement, 315 adding up to four extra three-dimensional arrays for the geopotential or isopycnal second order operator applied to 316 momentum. 317 Six CPP keys control the space variation of eddy coefficients: three for momentum and three for tracer. 318 The three choices allow: 319 a space variation in the three space directions (\key{traldf\_c3d}, \key{dynldf\_c3d}), 320 in the horizontal plane (\key{traldf\_c2d}, \key{dynldf\_c2d}), 321 or in the vertical only (\key{traldf\_c1d}, \key{dynldf\_c1d}). 322 The default option is a constant value over the whole ocean on both momentum and tracers. 323 324 The number of additional arrays that have to be defined and the gridpoint position at which 325 they are defined depend on both the space variation chosen and the type of operator used. 326 The resulting eddy viscosity and diffusivity coefficients can be a function of more than one variable. 327 Changes in the computer code when switching from one option to another have been minimized by 328 introducing the eddy coefficients as statement functions 329 (include file \textit{ldftra\_substitute.h90} and \textit{ldfdyn\_substitute.h90}). 330 The functions are replaced by their actual meaning during the preprocessing step (CPP). 331 The specification of the space variation of the coefficient is made in \mdl{ldftra} and \mdl{ldfdyn}, 332 or more precisely in include files \textit{traldf\_cNd.h90} and \textit{dynldf\_cNd.h90}, with N=1, 2 or 3. 333 The user can modify these include files as he/she wishes. 334 The way the mixing coefficient are set in the reference version can be briefly described as follows: 335 336 \subsubsection{Constant mixing coefficients (default option)} 337 When none of the \key{dynldf\_...} and \key{traldf\_...} keys are defined, 338 a constant value is used over the whole ocean for momentum and tracers, 339 which is specified through the \np{rn\_ahm0} and \np{rn\_aht0} namelist parameters. 340 341 \subsubsection{Vertically varying mixing coefficients (\protect\key{traldf\_c1d} and \key{dynldf\_c1d})} 342 The 1D option is only available when using the $z$-coordinate with full step. 343 Indeed in all the other types of vertical coordinate, 344 the depth is a 3D function of (\textbf{i},\textbf{j},\textbf{k}) and therefore, 345 introducing depth-dependent mixing coefficients will require 3D arrays. 346 In the 1D option, a hyperbolic variation of the lateral mixing coefficient is introduced in which 347 the surface value is \np{rn\_aht0} (\np{rn\_ahm0}), the bottom value is 1/4 of the surface value, 348 and the transition takes place around z=300~m with a width of 300~m 349 (\ie both the depth and the width of the inflection point are set to 300~m). 350 This profile is hard coded in file \textit{traldf\_c1d.h90}, but can be easily modified by users. 351 352 \subsubsection{Horizontally varying mixing coefficients (\protect\key{traldf\_c2d} and \protect\key{dynldf\_c2d})} 353 By default the horizontal variation of the eddy coefficient depends on the local mesh size and 336 The specification of the space variation of the coefficient is made in modules \mdl{ldftra} and \mdl{ldfdyn}. 337 The way the mixing coefficients are set in the reference version can be described as follows: 338 339 \subsection[Mixing coefficients read from file (\forcode{nn_aht_ijk_t = -20, -30}, \forcode{nn_ahm_ijk_t = -20,-30})] 340 { Mixing coefficients read from file (\protect\np{nn\_aht\_ijk\_t}\forcode{ = -20, -30}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = -20, -30})} 341 342 Mixing coefficients can be read from file if a particular geographical variation is needed. For example, in the ORCA2 global ocean model, 343 the laplacian viscosity operator uses $A^l$~= 4.10$^4$ m$^2$/s poleward of 20$^{\circ}$ north and south and 344 decreases linearly to $A^l$~= 2.10$^3$ m$^2$/s at the equator \citep{madec.delecluse.ea_JPO96, delecluse.madec_icol99}. 345 Similar modified horizontal variations can be found with the Antarctic or Arctic sub-domain options of ORCA2 and ORCA05. 346 The provided fields can either be 2d (\np{nn\_aht\_ijk\_t}\forcode{ = -20}, \np{nn\_ahm\_ijk\_t}\forcode{ = -20}) or 3d (\np{nn\_aht\_ijk\_t}\forcode{ = -30}, \np{nn\_ahm\_ijk\_t}\forcode{ = -30}). They must be given at U, V points for tracers and T, F points for momentum (see \autoref{tab:LDF_files}). 347 348 %-------------------------------------------------TABLE--------------------------------------------------- 349 \begin{table}[htb] 350 \begin{center} 351 \begin{tabular}{|l|l|l|l|} 352 \hline 353 Namelist parameter & Input filename & dimensions & variable names \\ \hline 354 \np{nn\_ahm\_ijk\_t}\forcode{ = -20} & \forcode{eddy_viscosity_2D.nc } & $(i,j)$ & \forcode{ahmt_2d, ahmf_2d} \\ \hline 355 \np{nn\_aht\_ijk\_t}\forcode{ = -20} & \forcode{eddy_diffusivity_2D.nc } & $(i,j)$ & \forcode{ahtu_2d, ahtv_2d} \\ \hline 356 \np{nn\_ahm\_ijk\_t}\forcode{ = -30} & \forcode{eddy_viscosity_3D.nc } & $(i,j,k)$ & \forcode{ahmt_3d, ahmf_3d} \\ \hline 357 \np{nn\_aht\_ijk\_t}\forcode{ = -30} & \forcode{eddy_diffusivity_3D.nc } & $(i,j,k)$ & \forcode{ahtu_3d, ahtv_3d} \\ \hline 358 \end{tabular} 359 \caption{ 360 \protect\label{tab:LDF_files} 361 Description of expected input files if mixing coefficients are read from NetCDF files. 362 } 363 \end{center} 364 \end{table} 365 %-------------------------------------------------------------------------------------------------------------- 366 367 \subsection[Constant mixing coefficients (\forcode{nn_aht_ijk_t = 0}, \forcode{nn_ahm_ijk_t = 0})] 368 { Constant mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 0}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 0})} 369 370 If constant, mixing coefficients are set thanks to a velocity and a length scales ($U_{scl}$, $L_{scl}$) such that: 371 372 \begin{equation} 373 \label{eq:constantah} 374 A_o^l = \left\{ 375 \begin{aligned} 376 & \frac{1}{2} U_{scl} L_{scl} & \text{for laplacian operator } \\ 377 & \frac{1}{12} U_{scl} L_{scl}^3 & \text{for bilaplacian operator } 378 \end{aligned} 379 \right. 380 \end{equation} 381 382 $U_{scl}$ and $L_{scl}$ are given by the namelist parameters \np{rn\_Ud}, \np{rn\_Uv}, \np{rn\_Ld} and \np{rn\_Lv}. 383 384 \subsection[Vertically varying mixing coefficients (\forcode{nn_aht_ijk_t = 10}, \forcode{nn_ahm_ijk_t = 10})] 385 {Vertically varying mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 10}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 10})} 386 387 In the vertically varying case, a hyperbolic variation of the lateral mixing coefficient is introduced in which 388 the surface value is given by \autoref{eq:constantah}, the bottom value is 1/4 of the surface value, 389 and the transition takes place around z=500~m with a width of 200~m. 390 This profile is hard coded in module \mdl{ldfc1d\_c2d}, but can be easily modified by users. 391 392 \subsection[Mesh size dependent mixing coefficients (\forcode{nn_aht_ijk_t = 20}, \forcode{nn_ahm_ijk_t = 20})] 393 {Mesh size dependent mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 20}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 20})} 394 395 In that case, the horizontal variation of the eddy coefficient depends on the local mesh size and 354 396 the type of operator used: 355 397 \begin{equation} … … 357 399 A_l = \left\{ 358 400 \begin{aligned} 359 & \frac{ \max(e_1,e_2)}{e_{max}} A_o^l& \text{for laplacian operator } \\360 & \frac{ \max(e_1,e_2)^{3}}{e_{max}^{3}} A_o^l& \text{for bilaplacian operator }401 & \frac{1}{2} U_{scl} \max(e_1,e_2) & \text{for laplacian operator } \\ 402 & \frac{1}{12} U_{scl} \max(e_1,e_2)^{3} & \text{for bilaplacian operator } 361 403 \end{aligned} 362 404 \right. 363 405 \end{equation} 364 where $e_{max}$ is the maximum of $e_1$ and $e_2$ taken over the whole masked ocean domain, 365 and $A_o^l$ is the \np{rn\_ahm0} (momentum) or \np{rn\_aht0} (tracer) namelist parameter. 406 where $U_{scl}$ is a user defined velocity scale (\np{rn\_Ud}, \np{rn\_Uv}). 366 407 This variation is intended to reflect the lesser need for subgrid scale eddy mixing where 367 408 the grid size is smaller in the domain. 368 It was introduced in the context of the DYNAMO modelling project \citep{ Willebrand_al_PO01}.369 Note that such a grid scale dependance of mixing coefficients significantly increase the range of stability of370 model configurations presenting large changes in grid pacing such as global ocean models.409 It was introduced in the context of the DYNAMO modelling project \citep{willebrand.barnier.ea_PO01}. 410 Note that such a grid scale dependance of mixing coefficients significantly increases the range of stability of 411 model configurations presenting large changes in grid spacing such as global ocean models. 371 412 Indeed, in such a case, a constant mixing coefficient can lead to a blow up of the model due to 372 413 large coefficient compare to the smallest grid size (see \autoref{sec:STP_forward_imp}), 373 414 especially when using a bilaplacian operator. 374 415 375 Other formulations can be introduced by the user for a given configuration. 376 For example, in the ORCA2 global ocean model (see Configurations), 377 the laplacian viscosity operator uses \np{rn\_ahm0}~= 4.10$^4$ m$^2$/s poleward of 20$^{\circ}$ north and south and 378 decreases linearly to \np{rn\_aht0}~= 2.10$^3$ m$^2$/s at the equator \citep{Madec_al_JPO96, Delecluse_Madec_Bk00}. 379 This modification can be found in routine \rou{ldf\_dyn\_c2d\_orca} defined in \mdl{ldfdyn\_c2d}. 380 Similar modified horizontal variations can be found with the Antarctic or Arctic sub-domain options of 381 ORCA2 and ORCA05 (see \&namcfg namelist). 382 383 \subsubsection{Space varying mixing coefficients (\protect\key{traldf\_c3d} and \key{dynldf\_c3d})} 384 385 The 3D space variation of the mixing coefficient is simply the combination of the 1D and 2D cases, 416 \colorbox{yellow}{CASE \np{nn\_aht\_ijk\_t} = 21 to be added} 417 418 \subsection[Mesh size and depth dependent mixing coefficients (\forcode{nn_aht_ijk_t = 30}, \forcode{nn_ahm_ijk_t = 30})] 419 {Mesh size and depth dependent mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 30}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 30})} 420 421 The 3D space variation of the mixing coefficient is simply the combination of the 1D and 2D cases above, 386 422 \ie a hyperbolic tangent variation with depth associated with a grid size dependence of 387 423 the magnitude of the coefficient. 388 424 389 \subsubsection{Space and time varying mixing coefficients} 390 391 There is no default specification of space and time varying mixing coefficient. 392 The only case available is specific to the ORCA2 and ORCA05 global ocean configurations. 393 It provides only a tracer mixing coefficient for eddy induced velocity (ORCA2) or both iso-neutral and 394 eddy induced velocity (ORCA05) that depends on the local growth rate of baroclinic instability. 395 This specification is actually used when an ORCA key and both \key{traldf\_eiv} and \key{traldf\_c2d} are defined. 425 \subsection[Velocity dependent mixing coefficients (\forcode{nn_aht_ijk_t = 31}, \forcode{nn_ahm_ijk_t = 31})] 426 {Flow dependent mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 31}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 31})} 427 In that case, the eddy coefficient is proportional to the local velocity magnitude so that the Reynolds number $Re = \lvert U \rvert e / A_l$ is constant (and here hardcoded to $12$): 428 \colorbox{yellow}{JC comment: The Reynolds is effectively set to 12 in the code for both operators but shouldn't it be 2 for Laplacian ?} 429 430 \begin{equation} 431 \label{eq:flowah} 432 A_l = \left\{ 433 \begin{aligned} 434 & \frac{1}{12} \lvert U \rvert e & \text{for laplacian operator } \\ 435 & \frac{1}{12} \lvert U \rvert e^3 & \text{for bilaplacian operator } 436 \end{aligned} 437 \right. 438 \end{equation} 439 440 \subsection[Deformation rate dependent viscosities (\forcode{nn_ahm_ijk_t = 32})] 441 {Deformation rate dependent viscosities (\protect\np{nn\_ahm\_ijk\_t}\forcode{ = 32})} 442 443 This option refers to the \citep{smagorinsky_MW63} scheme which is here implemented for momentum only. Smagorinsky chose as a 444 characteristic time scale $T_{smag}$ the deformation rate and for the lengthscale $L_{smag}$ the maximum wavenumber possible on the horizontal grid, e.g.: 445 446 \begin{equation} 447 \label{eq:smag1} 448 \begin{split} 449 T_{smag}^{-1} & = \sqrt{\left( \partial_x u - \partial_y v\right)^2 + \left( \partial_y u + \partial_x v\right)^2 } \\ 450 L_{smag} & = \frac{1}{\pi}\frac{e_1 e_2}{e_1 + e_2} 451 \end{split} 452 \end{equation} 453 454 Introducing a user defined constant $C$ (given in the namelist as \np{rn\_csmc}), one can deduce the mixing coefficients as follows: 455 456 \begin{equation} 457 \label{eq:smag2} 458 A_{smag} = \left\{ 459 \begin{aligned} 460 & C^2 T_{smag}^{-1} L_{smag}^2 & \text{for laplacian operator } \\ 461 & \frac{C^2}{8} T_{smag}^{-1} L_{smag}^4 & \text{for bilaplacian operator } 462 \end{aligned} 463 \right. 464 \end{equation} 465 466 For stability reasons, upper and lower limits are applied on the resulting coefficient (see \autoref{sec:STP_forward_imp}) so that: 467 \begin{equation} 468 \label{eq:smag3} 469 \begin{aligned} 470 & C_{min} \frac{1}{2} \lvert U \rvert e < A_{smag} < C_{max} \frac{e^2}{ 8\rdt} & \text{for laplacian operator } \\ 471 & C_{min} \frac{1}{12} \lvert U \rvert e^3 < A_{smag} < C_{max} \frac{e^4}{64 \rdt} & \text{for bilaplacian operator } 472 \end{aligned} 473 \end{equation} 474 475 where $C_{min}$ and $C_{max}$ are adimensional namelist parameters given by \np{rn\_minfac} and \np{rn\_maxfac} respectively. 476 477 \subsection{About space and time varying mixing coefficients} 396 478 397 479 The following points are relevant when the eddy coefficient varies spatially: … … 406 488 (\autoref{sec:dynldf_properties}). 407 489 408 (3) for isopycnal diffusion on momentum or tracers, an additional purely horizontal background diffusion with409 uniform coefficient can be added by setting a non zero value of \np{rn\_ahmb0} or \np{rn\_ahtb0},410 a background horizontal eddy viscosity or diffusivity coefficient411 (namelist parameters whose default values are $0$).412 However, the technique used to compute the isopycnal slopes is intended to get rid of such a background diffusion,413 since it introduces spurious diapycnal diffusion (see \autoref{sec:LDF_slp}).414 415 (4) when an eddy induced advection term is used (\key{traldf\_eiv}),416 $A^{eiv}$, the eddy induced coefficient has to be defined.417 Its space variations are controlled by the same CPP variable as for the eddy diffusivity coefficient418 (\ie \key{traldf\_cNd}).419 420 (5) the eddy coefficient associated with a biharmonic operator must be set to a \emph{negative} value.421 422 (6) it is possible to use both the laplacian and biharmonic operators concurrently.423 424 (7) it is possible to run without explicit lateral diffusion on momentum425 (\np{ln\_dynldf\_lap}\forcode{ = .?.}\np{ln\_dynldf\_bilap}\forcode{ = .false.}).426 This is recommended when using the UBS advection scheme on momentum (\np{ln\_dynadv\_ubs}\forcode{ = .true.},427 see \autoref{subsec:DYN_adv_ubs}) and can be useful for testing purposes.428 429 490 % ================================================================ 430 491 % Eddy Induced Mixing 431 492 % ================================================================ 432 \section{Eddy induced velocity (\protect\mdl{traadv\_eiv}, \protect\mdl{ldfeiv})} 493 \section[Eddy induced velocity (\forcode{ln_ldfeiv = .true.})] 494 {Eddy induced velocity (\protect\np{ln\_ldfeiv}\forcode{ = .true.})} 495 433 496 \label{sec:LDF_eiv} 497 498 %--------------------------------------------namtra_eiv--------------------------------------------------- 499 500 \nlst{namtra_eiv} 501 502 %-------------------------------------------------------------------------------------------------------------- 503 434 504 435 505 %%gm from Triad appendix&