darwin3 |
|
|
|
|||
Warning, /doc/algorithm/algorithm.rst is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit d4a066fa on 2025-09-10 18:05:35 UTC4f2617d475 Jeff*0001 .. _discret_algorithm: 0002 0003 Discretization and Algorithm 0004 **************************** 0005 0006 This chapter lays out the numerical schemes that are employed in the 0007 core MITgcm algorithm. Whenever possible links are made to actual 0008 program code in the MITgcm implementation. The chapter begins with a 0009 discussion of the temporal discretization used in MITgcm. This 0010 discussion is followed by sections that describe the spatial 0011 discretization. The schemes employed for momentum terms are described 0012 first, afterwards the schemes that apply to passive and dynamically 0013 active tracers are described. 0014 0015 Notation 0016 ======== 0017 0018 Because of the particularity of the vertical direction in stratified 0019 fluid context, in this chapter, the vector notations are mostly used for 0020 the horizontal component: the horizontal part of a vector is simply 0bad585a21 Navi*0021 written :math:`\vec{\bf v}` (instead of :math:`{{\bf v}_h}` or 4f2617d475 Jeff*0022 :math:`\vec{\mathbf{v}}_{h}` in chapter 1) and a 3D vector is simply 0023 written :math:`\vec{v}` (instead of :math:`\vec{\mathbf{v}}` in chapter 0024 1). 0025 0026 The notations we use to describe the discrete formulation of the model 0027 are summarized as follows. 0028 0029 General notation: 0030 0031 | :math:`\Delta x, \Delta y, \Delta r` grid spacing in X, Y, R directions 0032 | 0033 | :math:`A_c,A_w,A_s,A_{\zeta}` : horizontal area of a grid cell surrounding :math:`\theta,u,v,\zeta` point 0034 | 0035 | :math:`{\cal V}_u , {\cal V}_v , {\cal V}_w , {\cal V}_\theta` : Volume of the grid box surrounding :math:`u,v,w,\theta` point 0036 | 0037 | :math:`i,j,k` : current index relative to X, Y, R directions 0038 | 0039 0040 Basic operators: 0041 0042 | :math:`\delta_i` : 0043 :math:`\delta_i \Phi = \Phi_{i+1/2} - \Phi_{i-1/2}` 0044 | 0045 | :math:`~^{-i}` : 0046 :math:`\overline{\Phi}^i = ( \Phi_{i+1/2} + \Phi_{i-1/2} ) / 2` 0047 | 0048 | :math:`\delta_x` : 0049 :math:`\delta_x \Phi = \frac{1}{\Delta x} \delta_i \Phi` 0050 | 0bad585a21 Navi*0051 | :math:`\overline{ \nabla }` = horizontal gradient operator : 0052 :math:`\overline{ \nabla } \Phi = \{ \delta_x \Phi , \delta_y \Phi \}` 4f2617d475 Jeff*0053 | 0bad585a21 Navi*0054 | :math:`\overline{ \nabla } \cdot` = horizontal divergence operator : ab47de63dc Mart*0055 :math:`\overline{ \nabla }\cdot \vec{\mathrm{{\bf f}}} = 0056 \dfrac{1}{\cal A} \{ \delta_i \Delta y \, \mathrm{f}_x 4f2617d475 Jeff*0057 + \delta_j \Delta x \, \mathrm{f}_y \}` 0058 | 0059 | :math:`\overline{\nabla}^2` = horizontal Laplacian operator : ab47de63dc Mart*0060 :math:`\overline{\nabla}^2 \Phi = 0bad585a21 Navi*0061 \overline{ \nabla } \cdot \overline{ \nabla } \Phi` 4f2617d475 Jeff*0062 | 0063 68aadaa6bd Phob*0064 .. _time_stepping: 0065 4f2617d475 Jeff*0066 Time-stepping 0067 ============= 0068 0069 The equations of motion integrated by the model involve four prognostic 0070 equations for flow, :math:`u` and :math:`v`, temperature, 0071 :math:`\theta`, and salt/moisture, :math:`S`, and three diagnostic 0072 equations for vertical flow, :math:`w`, density/buoyancy, 0bad585a21 Navi*0073 :math:`\rho`/:math:`b`, and pressure/geo-potential, :math:`\phi_{\rm hyd}`. 4f2617d475 Jeff*0074 In addition, the surface pressure or height may by described by either a 0075 prognostic or diagnostic equation and if non-hydrostatics terms are 0076 included then a diagnostic equation for non-hydrostatic pressure is also 0077 solved. The combination of prognostic and diagnostic equations requires 0078 a model algorithm that can march forward prognostic variables while 0079 satisfying constraints imposed by diagnostic equations. 0080 0081 Since the model comes in several flavors and formulation, it would be 0082 confusing to present the model algorithm exactly as written into code 0083 along with all the switches and optional terms. Instead, we present the 0084 algorithm for each of the basic formulations which are: 0085 0086 #. the semi-implicit pressure method for hydrostatic equations with a 0087 rigid-lid, variables co-located in time and with Adams-Bashforth ab47de63dc Mart*0088 time-stepping; 4f2617d475 Jeff*0089 0090 #. as 1 but with an implicit linear free-surface; 0091 0092 #. as 1 or 2 but with variables staggered in time; 0093 0094 #. as 1 or 2 but with non-hydrostatic terms included; 0095 0096 #. as 2 or 3 but with non-linear free-surface. 0097 0098 In all the above configurations it is also possible to substitute the 0099 Adams-Bashforth with an alternative time-stepping scheme for terms 0100 evaluated explicitly in time. Since the over-arching algorithm is 0101 independent of the particular time-stepping scheme chosen we will 0102 describe first the over-arching algorithm, known as the pressure method, 0103 with a rigid-lid model in :numref:`press_meth_rigid`. This 0104 algorithm is essentially unchanged, apart for some coefficients, when 0105 the rigid lid assumption is replaced with a linearized implicit 0106 free-surface, described in :numref:`press_meth_linear`. These two flavors of the 0107 pressure-method encompass all formulations of the model as it exists 0108 today. The integration of explicit in time terms is out-lined in 0109 :numref:`adams-bashforth` and put into the context of the overall algorithm 0110 in :numref:`adams-bashforth-sync` and :numref:`adams-bashforth-staggered`. 0111 Inclusion of non-hydrostatic terms 0112 requires applying the pressure method in three dimensions instead of two 0113 and this algorithm modification is described in 0114 :numref:`non-hydrostatic`. Finally, the free-surface equation may be treated 0115 more exactly, including non-linear terms, and this is described in 0116 :numref:`nonlinear-freesurface`. 0117 0118 .. _press_meth_rigid: ab47de63dc Mart*0119 4f2617d475 Jeff*0120 Pressure method with rigid-lid 0121 ============================== 0122 0123 The horizontal momentum and continuity equations for the ocean 0124 (:eq:`eq-ocean-mom` and :eq:`eq-ocean-cont`), or for the atmosphere 0125 (:eq:`atmos-mom` and :eq:`atmos-cont`), can be summarized by: 0126 0127 .. math:: 0128 0129 \begin{aligned} 0bad585a21 Navi*0130 \partial_t u + g \partial_x \eta & = G_u \\ 0131 \partial_t v + g \partial_y \eta & = G_v \\ 0132 \partial_x u + \partial_y v + \partial_z w & = 0\end{aligned} 4f2617d475 Jeff*0133 0134 where we are adopting the oceanic notation for brevity. All terms in 0135 the momentum equations, except for surface pressure gradient, are 0136 encapsulated in the :math:`G` vector. The continuity equation, when 0137 integrated over the fluid depth, :math:`H`, and with the rigid-lid/no 0138 normal flow boundary conditions applied, becomes: 0139 0140 .. math:: 0141 \partial_x H \widehat{u} + \partial_y H \widehat{v} = 0 0142 :label: rigid-lid-continuity 0143 0144 Here, :math:`H\widehat{u} = \int_H u dz` is the depth integral of 0145 :math:`u`, similarly for :math:`H\widehat{v}`. The rigid-lid 0146 approximation sets :math:`w=0` at the lid so that it does not move but 0147 allows a pressure to be exerted on the fluid by the lid. The horizontal 0148 momentum equations and vertically integrated continuity equation are be 0149 discretized in time and space as follows: 0150 0151 .. math:: 0152 u^{n+1} + \Delta t g \partial_x \eta^{n+1} 0153 = u^{n} + \Delta t G_u^{(n+1/2)} 0154 :label: discrete-time-u 0155 0156 .. math:: 0157 v^{n+1} + \Delta t g \partial_y \eta^{n+1} 0158 = v^{n} + \Delta t G_v^{(n+1/2)} 0159 :label: discrete-time-v 0160 0161 .. math:: 0162 \partial_x H \widehat{u^{n+1}} 0163 + \partial_y H \widehat{v^{n+1}} = 0 0164 :label: discrete-time-cont-rigid-lid 0165 0166 As written here, terms on the LHS all involve time level :math:`n+1` 0167 and are referred to as implicit; the implicit backward time stepping 0168 scheme is being used. All other terms in the RHS are explicit in time. 0169 The thermodynamic quantities are integrated forward in time in parallel 0170 with the flow and will be discussed later. For the purposes of 0171 describing the pressure method it suffices to say that the hydrostatic 0172 pressure gradient is explicit and so can be included in the vector 0173 :math:`G`. 0174 0bad585a21 Navi*0175 Substituting the two momentum equations into the depth-integrated 4f2617d475 Jeff*0176 continuity equation eliminates :math:`u^{n+1}` and :math:`v^{n+1}` 0177 yielding an elliptic equation for :math:`\eta^{n+1}`. Equations 0178 :eq:`discrete-time-u`, :eq:`discrete-time-v` and 0179 :eq:`discrete-time-cont-rigid-lid` can then be re-arranged as follows: 0180 0181 .. math:: u^{*} = u^{n} + \Delta t G_u^{(n+1/2)} 0182 :label: ustar-rigid-lid 0183 0184 .. math:: v^{*} = v^{n} + \Delta t G_v^{(n+1/2)} 0185 :label: vstar-rigid-lid 0186 0187 .. math:: \partial_x \Delta t g H \partial_x \eta^{n+1} 0188 + \partial_y \Delta t g H \partial_y \eta^{n+1} 0189 = \partial_x H \widehat{u^{*}} ab47de63dc Mart*0190 + \partial_y H \widehat{v^{*}} 4f2617d475 Jeff*0191 :label: elliptic 0192 0193 .. math:: u^{n+1} = u^{*} - \Delta t g \partial_x \eta^{n+1} 0194 :label: un+1-rigid-lid 0195 0196 .. math:: v^{n+1} = v^{*} - \Delta t g \partial_y \eta^{n+1} 0197 :label: vn+1-rigid-lid 0198 0199 Equations :eq:`ustar-rigid-lid` to :eq:`vn+1-rigid-lid`, solved 0200 sequentially, represent the pressure method algorithm used in the model. 0201 The essence of the pressure method lies in the fact that any explicit 0202 prediction for the flow would lead to a divergence flow field so a 0203 pressure field must be found that keeps the flow non-divergent over each 0204 step of the integration. The particular location in time of the pressure 0205 field is somewhat ambiguous; in :numref:`pressure-method-rigid-lid` we 0206 depicted as co-located with the future flow field (time level 0207 :math:`n+1`) but it could equally have been drawn as staggered in time 0208 with the flow. 0209 0210 .. figure:: figs/pressure-method-rigid-lid.* 0211 :width: 70% 0212 :align: center 0213 :alt: pressure-method-rigid-lid 0214 :name: pressure-method-rigid-lid 0215 0216 A schematic of the evolution in time of the pressure method algorithm. A prediction for the flow variables at time level :math:`n+1` is made based only on the explicit terms, :math:`G^{(n+^1/_2)}`, and denoted :math:`u^*`, :math:`v^*`. Next, a pressure field is found such that :math:`u^{n+1}`, :math:`v^{n+1}` will be non-divergent. Conceptually, the :math:`*` quantities exist at time level :math:`n+1` but they are intermediate and only temporary. 0217 0218 0219 The correspondence to the code is as follows: 0220 0221 - the prognostic phase, equations :eq:`ustar-rigid-lid` and 0222 :eq:`vstar-rigid-lid`, stepping forward :math:`u^n` and :math:`v^n` to 0223 :math:`u^{*}` and :math:`v^{*}` is coded in :filelink:`timestep.F <model/src/timestep.F>` 0224 0225 - the vertical integration, :math:`H \widehat{u^*}` and :math:`H 0226 \widehat{v^*}`, divergence and inversion of the elliptic operator in 0227 equation :eq:`elliptic` is coded in :filelink:`solve_for_pressure.F <model/src/solve_for_pressure.F>` 0228 0229 - finally, the new flow field at time level :math:`n+1` given by 0230 equations :eq:`un+1-rigid-lid` and :eq:`vn+1-rigid-lid` is calculated 0231 in :filelink:`correction_step.F <model/src/correction_step.F>` 0232 0233 The calling tree for these routines is as follows: 0234 0235 .. _call-tree-press-meth: 0236 0237 .. admonition:: Pressure method calling tree 0238 :class: note 0239 0240 | :filelink:`FORWARD\_STEP <model/src/forward_step.F>` 0241 | :math:`\phantom{W}` :filelink:`DYNAMICS <model/src/dynamics.F>` 0242 | :math:`\phantom{WW}` :filelink:`TIMESTEP <model/src/timestep.F>` :math:`\phantom{xxxxxxxxxxxxxxxxxxxxxx}` :math:`u^*,v^*` :eq:`ustar-rigid-lid` , :eq:`vstar-rigid-lid` 0243 | :math:`\phantom{W}` :filelink:`SOLVE\_FOR\_PRESSURE <model/src/solve_for_pressure.F>` 0244 | :math:`\phantom{WW}` :filelink:`CALC\_DIV\_GHAT <model/src/calc_div_ghat.F>` :math:`\phantom{xxxxxxxxxxxxxxxx}` :math:`H\widehat{u^*},H\widehat{v^*}` :eq:`elliptic` 0245 | :math:`\phantom{WW}` :filelink:`CG2D <model/src/cg2d.F>` :math:`\phantom{xxxxxxxxxxxxxxxxxxxxxxxxxx}` :math:`\eta^{n+1}` :eq:`elliptic` 0246 | :math:`\phantom{W}` :filelink:`MOMENTUM\_CORRECTION\_STEP <model/src/momentum_correction_step.F>` 0247 | :math:`\phantom{WW}` :filelink:`CALC\_GRAD\_PHI\_SURF <model/src/calc_grad_phi_surf.F>` :math:`\phantom{xxxxxxxxxx}` :math:`\nabla \eta^{n+1}` 0248 | :math:`\phantom{WW}` :filelink:`CORRECTION\_STEP <model/src/correction_step.F>` :math:`\phantom{xxxxxxxxxxxxw}` :math:`u^{n+1},v^{n+1}` :eq:`un+1-rigid-lid` , :eq:`vn+1-rigid-lid` 0249 0250 In general, the horizontal momentum time-stepping can contain some terms 0251 that are treated implicitly in time, such as the vertical viscosity when 0bad585a21 Navi*0252 using the backward time-stepping scheme (:varlink:`implicitViscosity` ``=.TRUE.``). The method used to solve 4f2617d475 Jeff*0253 those implicit terms is provided in :numref:`implicit-backward-stepping`, and modifies equations 0254 :eq:`discrete-time-u` and :eq:`discrete-time-v` to give: 0255 0256 .. math:: 0257 0258 \begin{aligned} 0259 u^{n+1} - \Delta t \partial_z A_v \partial_z u^{n+1} 0bad585a21 Navi*0260 + \Delta t g \partial_x \eta^{n+1} & = u^{n} + \Delta t G_u^{(n+1/2)} 4f2617d475 Jeff*0261 \\ 0262 v^{n+1} - \Delta t \partial_z A_v \partial_z v^{n+1} 0bad585a21 Navi*0263 + \Delta t g \partial_y \eta^{n+1} & = v^{n} + \Delta t G_v^{(n+1/2)}\end{aligned} 4f2617d475 Jeff*0264 0265 0266 .. _press_meth_linear: 0267 0268 Pressure method with implicit linear free-surface 0269 ================================================= 0270 0271 The rigid-lid approximation filters out external gravity waves 0272 subsequently modifying the dispersion relation of barotropic Rossby 0273 waves. The discrete form of the elliptic equation has some zero 0274 eigenvalues which makes it a potentially tricky or inefficient problem 0275 to solve. 0276 0277 The rigid-lid approximation can be easily replaced by a linearization of 0278 the free-surface equation which can be written: 0279 0280 .. math:: 94151a9b18 Jeff*0281 \partial_t \eta + \partial_x H \widehat{u} + \partial_y H \widehat{v} = {\mathcal{P-E+R}} 4f2617d475 Jeff*0282 :label: linear-free-surface=P-E 0283 0bad585a21 Navi*0284 which differs from the depth-integrated continuity equation with 0285 rigid-lid :eq:`rigid-lid-continuity` by the time-dependent term and 4f2617d475 Jeff*0286 fresh-water source term. 0287 0288 Equation :eq:`discrete-time-cont-rigid-lid` in the rigid-lid pressure 0289 method is then replaced by the time discretization of 0290 :eq:`linear-free-surface=P-E` which is: 0291 0292 .. math:: 0293 \eta^{n+1} 0294 + \Delta t \partial_x H \widehat{u^{n+1}} 0295 + \Delta t \partial_y H \widehat{v^{n+1}} 94151a9b18 Jeff*0296 = \eta^{n} + \Delta t ( {\mathcal{P-E}}) 4f2617d475 Jeff*0297 :label: discrete-time-backward-free-surface 0298 0299 where the use of flow at time level :math:`n+1` makes the method 0300 implicit and backward in time. This is the preferred scheme since it 0301 still filters the fast unresolved wave motions by damping them. A d4a066fa68 Jean*0302 centered scheme, such as Crank-Nicolson (see 4f2617d475 Jeff*0303 :numref:`crank-nicolson_baro`), would alias the energy of the fast modes onto 0304 slower modes of motion. 0305 0306 As for the rigid-lid pressure method, equations :eq:`discrete-time-u`, 0307 :eq:`discrete-time-v` and :eq:`discrete-time-backward-free-surface` can be 0308 re-arranged as follows: 0309 0310 .. math:: 0311 u^{*} = u^{n} + \Delta t G_u^{(n+1/2)} 0312 :label: ustar-backward-free-surface 0313 0314 .. math:: 0315 v^{*} = v^{n} + \Delta t G_v^{(n+1/2)} 0316 :label: vstar-backward-free-surface 0317 0318 .. math:: 0bad585a21 Navi*0319 \eta^* = \epsilon_{\rm fs} ( \eta^{n} + \Delta t ({\mathcal{P-E}}) ) 4f2617d475 Jeff*0320 - \Delta t ( \partial_x H \widehat{u^{*}} 0321 + \partial_y H \widehat{v^{*}} ) 0322 :label: etastar-backward-free-surface 0323 0324 .. math:: 0325 \partial_x g H \partial_x \eta^{n+1} 0326 + \partial_y g H \partial_y \eta^{n+1} 0bad585a21 Navi*0327 - \frac{\epsilon_{\rm fs} \eta^{n+1}}{\Delta t^2} = 4f2617d475 Jeff*0328 - \frac{\eta^*}{\Delta t^2} 0329 :label: elliptic-backward-free-surface 0330 0331 .. math:: 0332 u^{n+1} = u^{*} - \Delta t g \partial_x \eta^{n+1} 0333 :label: un+1-backward-free-surface 0334 0335 .. math:: 0336 v^{n+1} = v^{*} - \Delta t g \partial_y \eta^{n+1} 0337 :label: vn+1-backward-free-surface 0338 0339 Equations :eq:`ustar-backward-free-surface` 0340 to :eq:`vn+1-backward-free-surface`, solved sequentially, represent the 0341 pressure method algorithm with a backward implicit, linearized free 0342 surface. The method is still formerly a pressure method because in the 0343 limit of large :math:`\Delta t` the rigid-lid method is recovered. 0344 However, the implicit treatment of the free-surface allows the flow to 0345 be divergent and for the surface pressure/elevation to respond on a 0346 finite time-scale (as opposed to instantly). To recover the rigid-lid dcaaa42497 Jeff*0347 formulation, we use a switch-like variable, 0bad585a21 Navi*0348 :math:`\epsilon_{\rm fs}` (:varlink:`freesurfFac`), which selects between the free-surface and 0349 rigid-lid; :math:`\epsilon_{\rm fs}=1` allows the free-surface to evolve; 0350 :math:`\epsilon_{\rm fs}=0` imposes the rigid-lid. The evolution in time and 4f2617d475 Jeff*0351 location of variables is exactly as it was for the rigid-lid model so 0352 that :numref:`pressure-method-rigid-lid` is still applicable. 0353 Similarly, the calling sequence, given :ref:`here <call-tree-press-meth>`, is as for the pressure-method. 0354 0355 .. _adams-bashforth: 0356 0357 Explicit time-stepping: Adams-Bashforth 0358 ======================================= 0359 0360 In describing the the pressure method above we deferred describing the 0361 time discretization of the explicit terms. We have historically used the 838416a165 Jeff*0362 quasi-second order Adams-Bashforth method (AB-II) for all explicit terms in both 4f2617d475 Jeff*0363 the momentum and tracer equations. This is still the default mode of 0364 operation but it is now possible to use alternate schemes for tracers ab47de63dc Mart*0365 (see :numref:`tracer_eqns`), or a 3rd order Adams-Bashforth method (AB-III). 838416a165 Jeff*0366 In the previous sections, we summarized an explicit scheme as: 4f2617d475 Jeff*0367 0368 .. math:: 0369 \tau^{*} = \tau^{n} + \Delta t G_\tau^{(n+1/2)} 0370 :label: taustar 0371 0372 where :math:`\tau` could be any prognostic variable (:math:`u`, 0373 :math:`v`, :math:`\theta` or :math:`S`) and :math:`\tau^*` is an 0374 explicit estimate of :math:`\tau^{n+1}` and would be exact if not for 0375 implicit-in-time terms. The parenthesis about :math:`n+1/2` indicates 838416a165 Jeff*0376 that the term is explicit and extrapolated forward in time. Below we describe 0377 in more detail the AB-II and AB-III schemes. 0378 0379 Adams-Bashforth II 0380 ------------------ 0381 0382 The quasi-second order Adams-Bashforth scheme is formulated as follows: 4f2617d475 Jeff*0383 0384 .. math:: 0bad585a21 Navi*0385 G_\tau^{(n+1/2)} = ( 3/2 + \epsilon_{\rm AB}) G_\tau^n 0386 - ( 1/2 + \epsilon_{\rm AB}) G_\tau^{n-1} 4f2617d475 Jeff*0387 :label: adams-bashforth2 0388 0389 This is a linear extrapolation, forward in time, to 0bad585a21 Navi*0390 :math:`t=(n+1/2+{\epsilon_{\rm AB}})\Delta t`. An extrapolation to the 4f2617d475 Jeff*0391 mid-point in time, :math:`t=(n+1/2)\Delta t`, corresponding to 0bad585a21 Navi*0392 :math:`\epsilon_{\rm AB}=0`, would be second order accurate but is weakly 4f2617d475 Jeff*0393 unstable for oscillatory terms. A small but finite value for 0bad585a21 Navi*0394 :math:`\epsilon_{\rm AB}` stabilizes the method. Strictly speaking, damping 4f2617d475 Jeff*0395 terms such as diffusion and dissipation, and fixed terms (forcing), do 0396 not need to be inside the Adams-Bashforth extrapolation. However, in the 0397 current code, it is simpler to include these terms and this can be 0398 justified if the flow and forcing evolves smoothly. Problems can, and 0399 do, arise when forcing or motions are high frequency and this 0400 corresponds to a reduced stability compared to a simple forward 0401 time-stepping of such terms. The model offers the possibility to leave 0402 terms outside the Adams-Bashforth extrapolation, by turning off the logical flag :varlink:`forcing_In_AB` 0bad585a21 Navi*0403 (parameter file ``data``, namelist ``PARM01``, default value = ``.TRUE.``) and then setting :varlink:`tracForcingOutAB` 4f2617d475 Jeff*0404 (default=0), :varlink:`momForcingOutAB` (default=0), and :varlink:`momDissip_In_AB` (parameter file ``data``, namelist ``PARM01``, 0405 default value = TRUE), respectively for the tracer terms, momentum forcing terms, and the dissipation terms. 0406 0407 A stability analysis for an oscillation equation should be given at this 0408 point. 0409 0410 A stability analysis for a relaxation equation should be given at this 0411 point. 0412 0413 .. figure:: figs/oscil+damp_AB2.* 0414 :width: 80% 0415 :align: center 0416 :alt: stability_analysis 0417 :name: oscil+damp_AB2 0418 0bad585a21 Navi*0419 Oscillatory and damping response of quasi-second order Adams-Bashforth scheme for different values of the :math:`\epsilon _{\rm AB}` 0420 parameter (0.0, 0.1, 0.25, from top to bottom) The analytical solution (in black), the physical mode (in blue) and the numerical 0421 mode (in red) are represented with a CFL step of 0.1. The left column represents the oscillatory response on the complex plane 0422 for CFL ranging from 0.1 up to 0.9. The right column represents the damping response amplitude (y-axis) function of the CFL (x-axis). 4f2617d475 Jeff*0423 838416a165 Jeff*0424 Adams-Bashforth III 0425 ------------------- 0426 0427 The 3rd order Adams-Bashforth time stepping (AB-III) provides several 0428 advantages (see, e.g., Durran 1991 :cite:`durran:91`) compared to the 0429 default quasi-second order Adams-Bashforth method: 0430 0431 - higher accuracy; 0432 0433 - stable with a longer time-step; 0434 0435 - no additional computation (just requires the storage of one 0436 additional time level). 0437 0438 The 3rd order Adams-Bashforth can be used to extrapolate 0439 forward in time the tendency (replacing :eq:`adams-bashforth2`) 0440 as: 0441 0442 .. math:: 0bad585a21 Navi*0443 G_\tau^{(n+1/2)} = ( 1 + \alpha_{\rm AB} + \beta_{\rm AB}) G_\tau^n 0444 - ( \alpha_{\rm AB} + 2 \beta_{\rm AB}) G_\tau^{n-1} 0445 + \beta_{\rm AB} G_\tau^{n-2} 838416a165 Jeff*0446 :label: adams-bashforth3 0447 0448 3rd order accuracy is obtained with 0bad585a21 Navi*0449 :math:`(\alpha_{\rm AB},\,\beta_{\rm AB}) = (1/2,\,5/12)`. Note that selecting 0450 :math:`(\alpha_{\rm AB},\,\beta_{\rm AB}) = (1/2+\epsilon_{AB},\,0)` one 838416a165 Jeff*0451 recovers AB-II. The AB-III time stepping improves the 0452 stability limit for an oscillatory problem like advection or Coriolis. 0453 As seen from :numref:`ab3_oscill_response`, it remains stable up to a 0454 CFL of 0.72, compared to only 0.50 with AB-II and 0bad585a21 Navi*0455 :math:`\epsilon_{\rm AB} = 0.1`. It is interesting to note that the 838416a165 Jeff*0456 stability limit can be further extended up to a CFL of 0.786 for an 0457 oscillatory problem (see :numref:`ab3_oscill_response`) using 0bad585a21 Navi*0458 :math:`(\alpha_{\rm AB},\,\beta_{\rm AB}) = (0.5,\,0.2811)` but then the scheme 838416a165 Jeff*0459 is only second order accurate. 0460 0461 However, the behavior of the AB-III for a damping problem (like diffusion) 0462 is less favorable, since the stability limit is reduced to 0.54 only 0bad585a21 Navi*0463 (and 0.64 with :math:`\beta_{\rm AB} = 0.2811`) compared to 1.0 (and 0.9 with 0464 :math:`\epsilon_{\rm AB} = 0.1`) with the AB-II (see 838416a165 Jeff*0465 :numref:`ab3_damp_response`). 0466 0467 A way to enable the use of a longer time step is to keep the dissipation 0468 terms outside the AB extrapolation (setting :varlink:`momDissip_In_AB` to ``.FALSE.`` 0469 in main parameter file ``data``, namelist ``PARM03``, thus returning to 0470 a simple forward time-stepping for dissipation, and to use AB-III only for 0471 advection and Coriolis terms. 0472 0473 The AB-III time stepping is activated by defining the option ``#define`` 0474 :varlink:`ALLOW_ADAMSBASHFORTH_3` in :filelink:`CPP_OPTIONS.h <model/inc/CPP_OPTIONS.h>`. The parameters 0bad585a21 Navi*0475 :math:`\alpha_{\rm AB},\beta_{\rm AB}` can be set from the main parameter file 838416a165 Jeff*0476 ``data`` (namelist ``PARM03``) and their default values correspond to 0477 the 3rd order Adams-Bashforth. A simple example is provided in 0478 :filelink:`verification/advect_xy/input.ab3_c4`. 0479 0480 AB-III is not yet available for the vertical momentum equation 0481 (non-hydrostatic) nor for passive tracers. 0482 0483 .. figure:: figs/stab_AB3_oscil.* 0484 :width: 80% 0485 :align: center 0486 :alt: ab3_stability_analysis 0487 :name: ab3_oscill_response 0488 0bad585a21 Navi*0489 Oscillatory response of third order Adams-Bashforth scheme for different values of the :math:`(\alpha_{\rm AB},\,\beta_{\rm AB})` parameters. 838416a165 Jeff*0490 The analytical solution (in black), the physical mode (in blue) and the numerical mode (in red) are represented with a CFL step of 0.1. 0491 0492 .. figure:: figs/stab_AB3_dampR.* 0493 :width: 80% 0494 :align: center 0495 :alt: ab3_damping_analysis 0496 :name: ab3_damp_response 0497 0bad585a21 Navi*0498 Damping response of third order Adams-Bashforth scheme for different values of the :math:`(\alpha_{\rm AB},\,\beta_{\rm AB})` parameters. 838416a165 Jeff*0499 The analytical solution (in black), the physical mode (in blue) and the numerical mode (in red) are represented with a CFL step of 0.1. 0500 0501 4f2617d475 Jeff*0502 .. _implicit-backward-stepping: 0503 0504 Implicit time-stepping: backward method 0505 ======================================= 0506 0507 Vertical diffusion and viscosity can be treated implicitly in time using 0508 the backward method which is an intrinsic scheme. Recently, the option 0509 to treat the vertical advection implicitly has been added, but not yet 0510 tested; therefore, the description hereafter is limited to diffusion and 0511 viscosity. For tracers, the time discretized equation is: 0512 0513 .. math:: 0514 \tau^{n+1} - \Delta t \partial_r \kappa_v \partial_r \tau^{n+1} = \tau^{n} + \Delta t G_\tau^{(n+1/2)} 0515 :label: implicit-diffusion 0516 0517 where :math:`G_\tau^{(n+1/2)}` is the remaining explicit terms 0518 extrapolated using the Adams-Bashforth method as described above. 0519 Equation :eq:`implicit-diffusion` can be split split into: 0520 0521 .. math:: 0522 \tau^* = \tau^{n} + \Delta t G_\tau^{(n+1/2)} 0523 :label: taustar-implicit 0524 0525 .. math:: 0526 \tau^{n+1} = {\cal L}_\tau^{-1} ( \tau^* ) 0527 :label: tau-n+1-implicit 0528 0529 where :math:`{\cal L}_\tau^{-1}` is the inverse of the operator 0530 0531 .. math:: {\cal L}_\tau = \left[ 1 + \Delta t \partial_r \kappa_v \partial_r \right] 0532 0533 Equation :eq:`taustar-implicit` looks exactly as :eq:`taustar` while 0534 :eq:`tau-n+1-implicit` involves an operator or matrix inversion. By 0535 re-arranging :eq:`implicit-diffusion` in this way we have cast the method 0536 as an explicit prediction step and an implicit step allowing the latter 0537 to be inserted into the over all algorithm with minimal interference. 0538 0539 The calling sequence for stepping forward a tracer variable such as temperature with implicit diffusion is 0540 as follows: 0541 0542 .. _adams-bash-calltree: 0543 0544 .. admonition:: Adams-Bashforth calling tree 0545 :class: note 0546 0547 | :filelink:`FORWARD\_STEP <model/src/forward_step.F>` 0548 | :math:`\phantom{W}` :filelink:`THERMODYNAMICS <model/src/thermodynamics.F>` 0549 | :math:`\phantom{WW}` :filelink:`TEMP\_INTEGRATE <model/src/temp_integrate.F>` 0550 | :math:`\phantom{WWW}` :filelink:`GAD\_CALC\_RHS <pkg/generic_advdiff/gad_calc_rhs.F>` :math:`\phantom{xxxxxxxxxw}` :math:`G_\theta^n = G_\theta( u, \theta^n)` 0551 | :math:`\phantom{WWW}` either 0552 | :math:`\phantom{WWWW}` :filelink:`EXTERNAL\_FORCING <model/src/external_forcing.F>` :math:`\phantom{xxxx}` :math:`G_\theta^n = G_\theta^n + {\cal Q}` 0553 | :math:`\phantom{WWWW}` :filelink:`ADAMS\_BASHFORTH2 <model/src/adams_bashforth2.F>` :math:`\phantom{xxi}` :math:`G_\theta^{(n+1/2)}` :eq:`adams-bashforth2` 0554 | :math:`\phantom{WWW}` or 0555 | :math:`\phantom{WWWW}` :filelink:`EXTERNAL\_FORCING <model/src/external_forcing.F>` :math:`\phantom{xxxx}` :math:`G_\theta^{(n+1/2)} = G_\theta^{(n+1/2)} + {\cal Q}` 0556 | :math:`\phantom{WW}` :filelink:`TIMESTEP\_TRACER <model/src/timestep_tracer.F>` :math:`\phantom{xxxxxxxxxx}` :math:`\tau^*` :eq:`taustar` 0557 | :math:`\phantom{WW}` :filelink:`IMPLDIFF <model/src/impldiff.F>` :math:`\phantom{xxxxxxxxxxxxxxxxxw}` :math:`\tau^{(n+1)}` :eq:`tau-n+1-implicit` 0558 0559 In order to fit within the pressure method, the implicit viscosity must 0560 not alter the barotropic flow. In other words, it can only redistribute 0561 momentum in the vertical. The upshot of this is that although vertical 0562 viscosity may be backward implicit and unconditionally stable, no-slip 0563 boundary conditions may not be made implicit and are thus cast as a an 0564 explicit drag term. 0565 0566 .. _adams-bashforth-sync: 0567 0568 Synchronous time-stepping: variables co-located in time 0569 ======================================================= 0570 0571 .. figure:: figs/adams-bashforth-sync.* 0572 :width: 70% 0573 :align: center 0574 :alt: adams-bash-sync 0575 :name: adams-bash-sync 0576 0577 A schematic of the explicit Adams-Bashforth and implicit time-stepping phases of the algorithm. All prognostic variables are co-located in time. Explicit tendencies are evaluated at time level :math:`n` as a function of the state at that time level (dotted arrow). The explicit tendency from the previous time level, :math:`n-1`, is used to extrapolate tendencies to :math:`n+1/2` (dashed arrow). This extrapolated tendency allows variables to be stably integrated forward-in-time to render an estimate (:math:`*` -variables) at the :math:`n+1` time level (solid arc-arrow). The operator :math:`{\cal L}` formed from implicit-in-time terms is solved to yield the state variables at time level :math:`n+1`. 0578 0579 0580 The Adams-Bashforth extrapolation of explicit tendencies fits neatly 0581 into the pressure method algorithm when all state variables are 0582 co-located in time. The algorithm can be represented by the sequential solution of the 0583 follow equations: 0584 ab47de63dc Mart*0585 .. math:: 4f2617d475 Jeff*0586 G_{\theta,S}^{n} = G_{\theta,S} ( u^n, \theta^n, S^n ) 0587 :label: Gt-n-sync 0588 0589 .. math:: 0590 G_{\theta,S}^{(n+1/2)} = (3/2+\epsilon_{AB}) G_{\theta,S}^{n}-(1/2+\epsilon_{AB}) G_{\theta,S}^{n-1} 0591 :label: Gt-n+5-sync 0592 0593 .. math:: 0594 (\theta^*,S^*) = (\theta^{n},S^{n}) + \Delta t G_{\theta,S}^{(n+1/2)} 0595 :label: tstar-sync 0596 0597 .. math:: 0598 (\theta^{n+1},S^{n+1}) = {\cal L}^{-1}_{\theta,S} (\theta^*,S^*) 0599 :label: t-n+1-sync 0600 0601 .. math:: 0bad585a21 Navi*0602 \phi^n_{\rm hyd} = \int b(\theta^n,S^n) dr 4f2617d475 Jeff*0603 :label: phi-hyd-sync 0604 0605 .. math:: 0bad585a21 Navi*0606 \vec{\bf G}_{\vec{\bf v}}^{n} = \vec{\bf G}_{\vec{\bf v}} ( \vec{\bf v}^n, \phi^n_{\rm hyd} ) 4f2617d475 Jeff*0607 :label: Gv-n-sync 0608 0609 .. math:: 0610 \vec{\bf G}_{\vec{\bf v}}^{(n+1/2)} = (3/2 + \epsilon_{AB} ) \vec{\bf G}_{\vec{\bf v}}^{n} - (1/2 + \epsilon_{AB} ) \vec{\bf G}_{\vec{\bf v}}^{n-1} 0611 :label: Gv-n+5-sync 0612 0613 .. math:: 0614 \vec{\bf v}^{*} = \vec{\bf v}^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}}^{(n+1/2)} 0615 :label: vstar-sync 0616 0617 .. math:: 0618 \vec{\bf v}^{**} = {\cal L}_{\vec{\bf v}}^{-1} ( \vec{\bf v}^* ) 0619 :label: vstarstar-sync 0620 0621 .. math:: 0bad585a21 Navi*0622 \eta^* = \epsilon_{\rm fs} \left( \eta^{n} + \Delta t ({\mathcal{P-E}}) \right)- \Delta t 0623 \nabla \cdot H \widehat{ \vec{\bf v}^{**} } 4f2617d475 Jeff*0624 :label: nstar-sync 0625 0626 .. math:: 0bad585a21 Navi*0627 \nabla \cdot g H \nabla \eta^{n+1} - \frac{\epsilon_{\rm fs} \eta^{n+1}}{\Delta t^2} ~ = ~ - \frac{\eta^*}{\Delta t^2} 4f2617d475 Jeff*0628 :label: elliptic-sync 0629 0630 .. math:: 0bad585a21 Navi*0631 \vec{\bf v}^{n+1} = \vec{\bf v}^{**} - \Delta t g \nabla \eta^{n+1} 4f2617d475 Jeff*0632 :label: v-n+1-sync 0633 0634 :numref:`adams-bash-sync` illustrates the location of variables 0635 in time and evolution of the algorithm with time. The Adams-Bashforth 0636 extrapolation of the tracer tendencies is illustrated by the dashed 0637 arrow, the prediction at :math:`n+1` is indicated by the solid arc. 0638 Inversion of the implicit terms, :math:`{\cal 0639 L}^{-1}_{\theta,S}`, then yields the new tracer fields at :math:`n+1`. 0640 All these operations are carried out in subroutine :filelink:`THERMODYNAMICS <model/src/thermodynamics.F>` and 0641 subsidiaries, which correspond to equations :eq:`Gt-n-sync` to 0642 :eq:`t-n+1-sync`. Similarly illustrated is the Adams-Bashforth 0643 extrapolation of accelerations, stepping forward and solving of implicit 0644 viscosity and surface pressure gradient terms, corresponding to 0645 equations :eq:`Gv-n-sync` to :eq:`v-n+1-sync`. These operations are 0646 carried out in subroutines :filelink:`DYNAMICS <model/src/dynamics.F>`, 0647 :filelink:`SOLVE\_FOR\_PRESSURE <model/src/solve_for_pressure.F>` and 0648 :filelink:`MOMENTUM\_CORRECTION\_STEP <model/src/momentum_correction_step.F>`. 0649 This, then, represents an entire algorithm 0650 for stepping forward the model one time-step. The corresponding calling 0651 tree for the overall synchronous algorithm using 0652 Adams-Bashforth time-stepping is given below. The place where the model geometry 0653 :varlink:`hFac` factors) is updated is added here but is only relevant 0654 for the non-linear free-surface algorithm. 0655 For completeness, the external forcing, 0656 ocean and atmospheric physics have been added, although they are mainly optional. 0657 0658 .. admonition:: Synchronous Adams-Bashforth calling tree 0659 :class: note 0660 0661 | :filelink:`FORWARD\_STEP <model/src/forward_step.F>` 0662 | :math:`\phantom{WWW}` :filelink:`EXTERNAL\_FIELDS\_LOAD <model/src/external_fields_load.F>` 0663 | :math:`\phantom{WWW}` :filelink:`DO\_ATMOSPHERIC\_PHYS <model/src/do_atmospheric_phys.F>` 0664 | :math:`\phantom{WWW}` :filelink:`DO\_OCEANIC\_PHYS <model/src/do_oceanic_phys.F>` 0665 | :math:`\phantom{WW}` :filelink:`THERMODYNAMICS <model/src/thermodynamics.F>` 0666 | :math:`\phantom{WWW}` :filelink:`CALC\_GT <model/src/calc_gt.F>` 0667 | :math:`\phantom{WWWW}` :filelink:`GAD\_CALC\_RHS <pkg/generic_advdiff/gad_calc_rhs.F>` :math:`\phantom{xxxxxxxxxxxxxlwww}` :math:`G_\theta^n = G_\theta( u, \theta^n )` :eq:`Gt-n-sync` 0668 | :math:`\phantom{WWWW}` :filelink:`EXTERNAL\_FORCING <model/src/external_forcing.F>` :math:`\phantom{xxxxxxxxxxlww}` :math:`G_\theta^n = G_\theta^n + {\cal Q}` 0669 | :math:`\phantom{WWWW}` :filelink:`ADAMS\_BASHFORTH2 <model/src/adams_bashforth2.F>` :math:`\phantom{xxxxxxxxxxxw}` :math:`G_\theta^{(n+1/2)}` :eq:`Gt-n+5-sync` 0670 | :math:`\phantom{WWW}` :filelink:`TIMESTEP\_TRACER <model/src/timestep_tracer.F>` :math:`\phantom{xxxxxxxxxxxxxxxww}` :math:`\theta^*` :eq:`tstar-sync` 0671 | :math:`\phantom{WWW}` :filelink:`IMPLDIFF <model/src/impldiff.F>` :math:`\phantom{xxxxxxxxxxxxxxxxxxxxxvwww}` :math:`\theta^{(n+1)}` :eq:`t-n+1-sync` 0672 | :math:`\phantom{WW}` :filelink:`DYNAMICS <model/src/dynamics.F>` 0bad585a21 Navi*0673 | :math:`\phantom{WWW}` :filelink:`CALC\_PHI\_HYD <model/src/calc_phi_hyd.F>` :math:`\phantom{xxxxxxxxxxxxxxxxxxxxi}` :math:`\phi_{\rm hyd}^n` :eq:`phi-hyd-sync` 4f2617d475 Jeff*0674 | :math:`\phantom{WWW}` :filelink:`MOM\_FLUXFORM <pkg/mom_fluxform/mom_fluxform.F>` or :filelink:`MOM\_VECINV <pkg/mom_vecinv/mom_vecinv.F>` :math:`\phantom{xxi}` :math:`G_{\vec{\bf v}}^n` :eq:`Gv-n-sync` 0675 | :math:`\phantom{WWW}` :filelink:`TIMESTEP <model/src/timestep.F>` :math:`\phantom{xxxxxxxxxxxxxxxxxxxxxxxxxx}` :math:`\vec{\bf v}^*` :eq:`Gv-n+5-sync`, :eq:`vstar-sync` 0676 | :math:`\phantom{WWW}` :filelink:`IMPLDIFF <model/src/impldiff.F>` :math:`\phantom{xxxxxxxxxxxxxxxxxxxxxxxxlw}` :math:`\vec{\bf v}^{**}` :eq:`vstarstar-sync` 0677 | :math:`\phantom{WW}` :filelink:`UPDATE\_R\_STAR <model/src/update_r_star.F>` or :filelink:`UPDATE\_SURF\_DR <model/src/update_surf_dr.F>` (NonLin-FS only) 0678 | :math:`\phantom{WW}` :filelink:`SOLVE\_FOR\_PRESSURE <model/src/solve_for_pressure.F>` 0679 | :math:`\phantom{WWW}` :filelink:`CALC\_DIV\_GHAT <model/src/calc_div_ghat.F>` :math:`\phantom{xxxxxxxxxxxxxxxxxxxx}` :math:`\eta^*` :eq:`nstar-sync` 0680 | :math:`\phantom{WWW}` :filelink:`CG2D <model/src/cg2d.F>` :math:`\phantom{xxxxxxxxxxxxxxxxxxxxxxxxxxxxxi}` :math:`\eta^{n+1}` :eq:`elliptic-sync` 0681 | :math:`\phantom{WW}` :filelink:`MOMENTUM\_CORRECTION\_STEP <model/src/momentum_correction_step.F>` 0682 | :math:`\phantom{WWW}` :filelink:`CALC\_GRAD\_PHI\_SURF <model/src/calc_grad_phi_surf.F>` :math:`\phantom{xxxxxxxxxxxxxx}` :math:`\nabla \eta^{n+1}` 0683 | :math:`\phantom{WWW}` :filelink:`CORRECTION\_STEP <model/src/correction_step.F>` :math:`\phantom{xxxxxxxxxxxxxxxxw}` :math:`u^{n+1},v^{n+1}` :eq:`v-n+1-sync` 0684 | :math:`\phantom{WW}` :filelink:`TRACERS\_CORRECTION\_STEP <model/src/tracers_correction_step.F>` 0685 | :math:`\phantom{WWW}` :filelink:`CYCLE\_TRACER <model/src/cycle_tracer.F>` :math:`\phantom{xxxxxxxxxxxxxxxxxxxxx}` :math:`\theta^{n+1}` 0686 | :math:`\phantom{WWW}` :filelink:`SHAP_FILT_APPLY_TS <pkg/shap_filt/shap_filt_apply_ts.F>` or :filelink:`ZONAL_FILT_APPLY_TS <pkg/zonal_filt/zonal_filt_apply_ts.F>` 0687 | :math:`\phantom{WWW}` :filelink:`CONVECTIVE_ADJUSTMENT <model/src/convective_adjustment.F>` 0688 0689 0690 .. _adams-bashforth-staggered: 0691 0692 Staggered baroclinic time-stepping 0693 ================================== 0694 0695 .. figure:: figs/adams-bashforth-staggered.* 0696 :width: 80% 0697 :align: center 0698 :alt: adams-bash-staggered 0699 :name: adams-bash-staggered 0700 0bad585a21 Navi*0701 A schematic of the explicit Adams-Bashforth and implicit time-stepping phases of the algorithm but with staggering in time of thermodynamic variables with the flow. Explicit momentum tendencies are evaluated at time level :math:`n-1/2` as a function of the flow field at that time level :math:`n-1/2`. The explicit tendency from the previous time level, :math:`n-3/2`, is used to extrapolate tendencies to :math:`n` (dashed arrow). The hydrostatic pressure/geo-potential :math:`\phi _{\rm hyd}` is evaluated directly at time level :math:`n` (vertical arrows) and used with the extrapolated tendencies to step forward the flow variables from :math:`n-1/2` to :math:`n+1/2` (solid arc-arrow). The implicit-in-time operator :math:`{\cal L}_{\bf u,v}` (vertical arrows) is then applied to the previous estimation of the the flow field (:math:`*` -variables) and yields to the two velocity components :math:`u,v` at time level :math:`n+1/2`. These are then used to calculate the advection term (dashed arc-arrow) of the thermo-dynamics tendencies at time step :math:`n`. The extrapolated thermodynamics tendency, from time level :math:`n-1` and :math:`n` to :math:`n+1/2`, allows thermodynamic variables to be stably integrated forward-in-time (solid arc-arrow) up to time level :math:`n+1`. 4f2617d475 Jeff*0702 0703 For well-stratified problems, internal gravity waves may be the limiting 0704 process for determining a stable time-step. In the circumstance, it is 0705 more efficient to stagger in time the thermodynamic variables with the 0706 flow variables. :numref:`adams-bash-staggered` illustrates the 0707 staggering and algorithm. The key difference between this and 0708 :numref:`adams-bash-sync` is that the thermodynamic variables are 0709 solved after the dynamics, using the recently updated flow field. This 0710 essentially allows the gravity wave terms to leap-frog in time giving 0711 second order accuracy and more stability. 0712 0713 The essential change in the staggered algorithm is that the 0714 thermodynamics solver is delayed from half a time step, allowing the use 0715 of the most recent velocities to compute the advection terms. Once the 0716 thermodynamics fields are updated, the hydrostatic pressure is computed 0717 to step forward the dynamics. Note that the pressure gradient must also 0718 be taken out of the Adams-Bashforth extrapolation. Also, retaining the 0719 integer time-levels, :math:`n` and :math:`n+1`, does not give a user the 0720 sense of where variables are located in time. Instead, we re-write the 0721 entire algorithm, :eq:`Gt-n-sync` to :eq:`v-n+1-sync`, annotating the 0722 position in time of variables appropriately: 0723 0724 .. math:: 0bad585a21 Navi*0725 \phi^{n}_{\rm hyd} = \int b(\theta^{n},S^{n}) dr 4f2617d475 Jeff*0726 :label: phi-hyd-staggered 0727 0728 .. math:: 0729 \vec{\bf G}_{\vec{\bf v}}^{n-1/2} = \vec{\bf G}_{\vec{\bf v}} ( \vec{\bf v}^{n-1/2} ) 0730 :label: Gv-n-staggered 0731 0732 .. math:: 0733 \vec{\bf G}_{\vec{\bf v}}^{(n)} = (3/2 + \epsilon_{AB} ) \vec{\bf G}_{\vec{\bf v}}^{n-1/2} - (1/2 + \epsilon_{AB} ) \vec{\bf G}_{\vec{\bf v}}^{n-3/2} 0734 :label: Gv-n+5-staggered 0735 0736 .. math:: 0bad585a21 Navi*0737 \vec{\bf v}^{*} = \vec{\bf v}^{n-1/2} + \Delta t \left( \vec{\bf G}_{\vec{\bf v}}^{(n)} - \nabla \phi_{\rm hyd}^{n} \right) 4f2617d475 Jeff*0738 :label: vstar-staggered 0739 0740 .. math:: 0741 \vec{\bf v}^{**} = {\cal L}_{\vec{\bf v}}^{-1} ( \vec{\bf v}^* ) 0742 :label: vstarstar-staggered 0743 0744 .. math:: 0bad585a21 Navi*0745 \eta^* = \epsilon_{\rm fs} \left( \eta^{n-1/2} + \Delta t ({\mathcal{P-E}})^n \right)- \Delta t 0746 \nabla \cdot H \widehat{ \vec{\bf v}^{**} } 4f2617d475 Jeff*0747 :label: nstar-staggered 0748 0749 .. math:: 0bad585a21 Navi*0750 \nabla \cdot g H \nabla \eta^{n+1/2} - \frac{\epsilon_{\rm fs} \eta^{n+1/2}}{\Delta t^2} 0751 = - \frac{\eta^*}{\Delta t^2} 4f2617d475 Jeff*0752 :label: elliptic-staggered 0753 0754 .. math:: 0bad585a21 Navi*0755 \vec{\bf v}^{n+1/2} = \vec{\bf v}^{**} - \Delta t g \nabla \eta^{n+1/2} 4f2617d475 Jeff*0756 :label: v-n+1-staggered 0757 0758 .. math:: 0759 G_{\theta,S}^{n} = G_{\theta,S} ( u^{n+1/2}, \theta^{n}, S^{n} ) 0760 :label: Gt-n-staggered 0761 0762 .. math:: 0763 G_{\theta,S}^{(n+1/2)} = (3/2+\epsilon_{AB}) G_{\theta,S}^{n}-(1/2+\epsilon_{AB}) G_{\theta,S}^{n-1} 0764 :label: Gt-n+5-staggered 0765 0766 .. math:: 0767 (\theta^*,S^*) = (\theta^{n},S^{n}) + \Delta t G_{\theta,S}^{(n+1/2)} 0768 :label: tstar-staggered 0769 0770 .. math:: 0771 (\theta^{n+1},S^{n+1}) = {\cal L}^{-1}_{\theta,S} (\theta^*,S^*) 0772 :label: t-n+1-staggered 0773 0774 The corresponding calling tree is given below. The staggered algorithm is 0bad585a21 Navi*0775 activated with the run-time flag :varlink:`staggerTimeStep` ``=.TRUE.`` in 4f2617d475 Jeff*0776 parameter file ``data``, namelist ``PARM01``. 0777 0778 .. admonition:: Staggered Adams-Bashforth calling tree 0779 :class: note 0780 0781 | :filelink:`FORWARD\_STEP <model/src/forward_step.F>` 0782 | :math:`\phantom{WWW}` :filelink:`EXTERNAL\_FIELDS\_LOAD <model/src/external_fields_load.F>` 0783 | :math:`\phantom{WWW}` :filelink:`DO\_ATMOSPHERIC\_PHYS <model/src/do_atmospheric_phys.F>` 0784 | :math:`\phantom{WWW}` :filelink:`DO\_OCEANIC\_PHYS <model/src/do_oceanic_phys.F>` 0785 | :math:`\phantom{WW}` :filelink:`DYNAMICS <model/src/dynamics.F>` 0bad585a21 Navi*0786 | :math:`\phantom{WWW}` :filelink:`CALC\_PHI\_HYD <model/src/calc_phi_hyd.F>` :math:`\phantom{xxxxxxxxxxxxxxxxxxxxi}` :math:`\phi_{\rm hyd}^n` :eq:`phi-hyd-staggered` 4f2617d475 Jeff*0787 | :math:`\phantom{WWW}` :filelink:`MOM\_FLUXFORM <pkg/mom_fluxform/mom_fluxform.F>` or :filelink:`MOM\_VECINV <pkg/mom_vecinv/mom_vecinv.F>` :math:`\phantom{xxi}` :math:`G_{\vec{\bf v}}^{n-1/2}` :eq:`Gv-n-staggered` 0788 | :math:`\phantom{WWW}` :filelink:`TIMESTEP <model/src/timestep.F>` :math:`\phantom{xxxxxxxxxxxxxxxxxxxxxxxxxx}` :math:`\vec{\bf v}^*` :eq:`Gv-n+5-staggered`, :eq:`vstar-staggered` 0789 | :math:`\phantom{WWW}` :filelink:`IMPLDIFF <model/src/impldiff.F>` :math:`\phantom{xxxxxxxxxxxxxxxxxxxxxxxxlw}` :math:`\vec{\bf v}^{**}` :eq:`vstarstar-staggered` 0790 | :math:`\phantom{WW}` :filelink:`UPDATE\_R\_STAR <model/src/update_r_star.F>` or :filelink:`UPDATE\_SURF\_DR <model/src/update_surf_dr.F>` (NonLin-FS only) 0791 | :math:`\phantom{WW}` :filelink:`SOLVE\_FOR\_PRESSURE <model/src/solve_for_pressure.F>` 0792 | :math:`\phantom{WWW}` :filelink:`CALC\_DIV\_GHAT <model/src/calc_div_ghat.F>` :math:`\phantom{xxxxxxxxxxxxxxxxxxxx}` :math:`\eta^*` :eq:`nstar-staggered` 0793 | :math:`\phantom{WWW}` :filelink:`CG2D <model/src/cg2d.F>` :math:`\phantom{xxxxxxxxxxxxxxxxxxxxxxxxxxxxxi}` :math:`\eta^{n+1/2}` :eq:`elliptic-staggered` 0794 | :math:`\phantom{WW}` :filelink:`MOMENTUM\_CORRECTION\_STEP <model/src/momentum_correction_step.F>` 0795 | :math:`\phantom{WWW}` :filelink:`CALC\_GRAD\_PHI\_SURF <model/src/calc_grad_phi_surf.F>` :math:`\phantom{xxxxxxxxxxxxxx}` :math:`\nabla \eta^{n+1/2}` 0796 | :math:`\phantom{WWW}` :filelink:`CORRECTION\_STEP <model/src/correction_step.F>` :math:`\phantom{xxxxxxxxxxxxxxxxw}` :math:`u^{n+1/2},v^{n+1/2}` :eq:`v-n+1-staggered` 0797 | :math:`\phantom{WW}` :filelink:`THERMODYNAMICS <model/src/thermodynamics.F>` 0798 | :math:`\phantom{WWW}` :filelink:`CALC\_GT <model/src/calc_gt.F>` 0799 | :math:`\phantom{WWWW}` :filelink:`GAD\_CALC\_RHS <pkg/generic_advdiff/gad_calc_rhs.F>` :math:`\phantom{xxxxxxxxxxxxxlwww}` :math:`G_\theta^n = G_\theta( u, \theta^n )` :eq:`Gt-n-staggered` 0800 | :math:`\phantom{WWWW}` :filelink:`EXTERNAL\_FORCING <model/src/external_forcing.F>` :math:`\phantom{xxxxxxxxxxlww}` :math:`G_\theta^n = G_\theta^n + {\cal Q}` 0801 | :math:`\phantom{WWWW}` :filelink:`ADAMS\_BASHFORTH2 <model/src/adams_bashforth2.F>` :math:`\phantom{xxxxxxxxxxxw}` :math:`G_\theta^{(n+1/2)}` :eq:`Gt-n+5-staggered` 0802 | :math:`\phantom{WWW}` :filelink:`TIMESTEP\_TRACER <model/src/timestep_tracer.F>` :math:`\phantom{xxxxxxxxxxxxxxxww}` :math:`\theta^*` :eq:`tstar-staggered` 0803 | :math:`\phantom{WWW}` :filelink:`IMPLDIFF <model/src/impldiff.F>` :math:`\phantom{xxxxxxxxxxxxxxxxxxxxxvwww}` :math:`\theta^{(n+1)}` :eq:`t-n+1-staggered` 0804 | :math:`\phantom{WW}` :filelink:`TRACERS\_CORRECTION\_STEP <model/src/tracers_correction_step.F>` 0805 | :math:`\phantom{WWW}` :filelink:`CYCLE\_TRACER <model/src/cycle_tracer.F>` :math:`\phantom{xxxxxxxxxxxxxxxxxxxxx}` :math:`\theta^{n+1}` 0806 | :math:`\phantom{WWW}` :filelink:`SHAP_FILT_APPLY_TS <pkg/shap_filt/shap_filt_apply_ts.F>` or :filelink:`ZONAL_FILT_APPLY_TS <pkg/zonal_filt/zonal_filt_apply_ts.F>` 0807 | :math:`\phantom{WWW}` :filelink:`CONVECTIVE_ADJUSTMENT <model/src/convective_adjustment.F>` 0808 0809 The only difficulty with this approach is apparent in equation 0810 :eq:`Gt-n-staggered` and illustrated by the dotted arrow connecting 0811 :math:`u,v^{n+1/2}` with :math:`G_\theta^{n}`. The flow used to advect 0812 tracers around is not naturally located in time. This could be avoided 0813 by applying the Adams-Bashforth extrapolation to the tracer field itself** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 815.
0814 and advecting that around but this approach is not yet available. We’re 0815 not aware of any detrimental effect of this feature. The difficulty lies 0816 mainly in interpretation of what time-level variables and terms 0817 correspond to. 0818 0819 .. _non-hydrostatic: 0820 0821 Non-hydrostatic formulation 0822 =========================== 0823 0824 0825 The non-hydrostatic formulation re-introduces the full vertical momentum 0826 equation and requires the solution of a 3-D elliptic equations for 0827 non-hydrostatic pressure perturbation. We still integrate vertically for 0828 the hydrostatic pressure and solve a 2-D elliptic equation for the 0829 surface pressure/elevation for this reduces the amount of work needed to 0830 solve for the non-hydrostatic pressure. 0831 0832 The momentum equations are discretized in time as follows: 0833 0834 .. math:: 0bad585a21 Navi*0835 \frac{1}{\Delta t} u^{n+1} + g \partial_x \eta^{n+1} + \partial_x \phi_{\rm nh}^{n+1} 4f2617d475 Jeff*0836 = \frac{1}{\Delta t} u^{n} + G_u^{(n+1/2)} 0837 :label: discrete-time-u-nh 0838 0839 .. math:: 0bad585a21 Navi*0840 \frac{1}{\Delta t} v^{n+1} + g \partial_y \eta^{n+1} + \partial_y \phi_{\rm nh}^{n+1} 4f2617d475 Jeff*0841 = \frac{1}{\Delta t} v^{n} + G_v^{(n+1/2)} 0842 :label: discrete-time-v-nh 0843 0844 .. math:: 0bad585a21 Navi*0845 \frac{1}{\Delta t} w^{n+1} + \partial_r \phi_{\rm nh}^{n+1} 4f2617d475 Jeff*0846 = \frac{1}{\Delta t} w^{n} + G_w^{(n+1/2)} 0847 :label: discrete-time-w-nh 0848 0849 which must satisfy the discrete-in-time depth integrated continuity, 0850 equation :eq:`discrete-time-backward-free-surface` and the local 0851 continuity equation 0852 0853 .. math:: 0854 \partial_x u^{n+1} + \partial_y v^{n+1} + \partial_r w^{n+1} = 0 0855 :label: non-divergence-nh 0856 0857 As before, the explicit predictions for momentum are consolidated as: 0858 0859 .. math:: 0860 0861 \begin{aligned} 0bad585a21 Navi*0862 u^* & = u^n + \Delta t G_u^{(n+1/2)} \\ 0863 v^* & = v^n + \Delta t G_v^{(n+1/2)} \\ 0864 w^* & = w^n + \Delta t G_w^{(n+1/2)}\end{aligned} 4f2617d475 Jeff*0865 0866 but this time we introduce an intermediate step by splitting the 0867 tendency of the flow as follows: 0868 0869 .. math:: 0870 0871 \begin{aligned} 0bad585a21 Navi*0872 u^{n+1} = u^{**} - \Delta t \partial_x \phi_{\rm nh}^{n+1} 4f2617d475 Jeff*0873 & & 0874 u^{**} = u^{*} - \Delta t g \partial_x \eta^{n+1} \\ 0bad585a21 Navi*0875 v^{n+1} = v^{**} - \Delta t \partial_y \phi_{\rm nh}^{n+1} 4f2617d475 Jeff*0876 & & 0877 v^{**} = v^{*} - \Delta t g \partial_y \eta^{n+1}\end{aligned} 0878 0879 Substituting into the depth integrated continuity 0bad585a21 Navi*0880 :eq:`discrete-time-backward-free-surface` gives 4f2617d475 Jeff*0881 0882 .. math:: 0bad585a21 Navi*0883 \partial_x H \partial_x \left( g \eta^{n+1} + \widehat{\phi}_{\rm nh}^{n+1} \right) 0884 + \partial_y H \partial_y \left( g \eta^{n+1} + \widehat{\phi}_{\rm nh}^{n+1} \right) 0885 - \frac{\epsilon_{\rm fs}\eta^{n+1}}{\Delta t^2} 4f2617d475 Jeff*0886 = - \frac{\eta^*}{\Delta t^2} 0887 :label: substituting-in-cont 0888 0889 which is approximated by equation :eq:`elliptic-backward-free-surface` 0bad585a21 Navi*0890 on the basis that i) :math:`\phi_{\rm nh}^{n+1}` is not yet known and ii) 0891 :math:`\nabla \widehat{\phi}_{\rm nh} \ll g \nabla \eta`. 0892 If :eq:`elliptic-backward-free-surface` is solved 0893 accurately then the implication is that :math:`\widehat{\phi}_{\rm nh} 4f2617d475 Jeff*0894 \approx 0` so that the non-hydrostatic pressure field does not drive 0895 barotropic motion. 0896 0897 The flow must satisfy non-divergence (equation :eq:`non-divergence-nh`) 0898 locally, as well as depth integrated, and this constraint is used to 0bad585a21 Navi*0899 form a 3-D elliptic equations for :math:`\phi_{\rm nh}^{n+1}`: 4f2617d475 Jeff*0900 0901 .. math:: 0bad585a21 Navi*0902 \partial_{xx} \phi_{\rm nh}^{n+1} + \partial_{yy} \phi_{\rm nh}^{n+1} + 0903 \partial_{rr} \phi_{\rm nh}^{n+1} = 4f2617d475 Jeff*0904 \partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*} 0905 :label: elliptic-pnh 0906 0907 The entire algorithm can be summarized as the sequential solution of the 0908 following equations: 0909 0910 .. math:: 0911 u^{*} = u^{n} + \Delta t G_u^{(n+1/2)} 0912 :label: ustar-nh 0913 ab47de63dc Mart*0914 .. math:: 4f2617d475 Jeff*0915 v^{*} = v^{n} + \Delta t G_v^{(n+1/2)} 0916 :label: vstar-nh ab47de63dc Mart*0917 4f2617d475 Jeff*0918 .. math:: 0919 w^{*} = w^{n} + \Delta t G_w^{(n+1/2)} 0920 :label: wstar-nh 0921 0922 .. math:: 0bad585a21 Navi*0923 \eta^* ~ = ~ \epsilon_{\rm fs} \left( \eta^{n} + \Delta t ({\mathcal{P-E}}) \right) 4f2617d475 Jeff*0924 - \Delta t \left( \partial_x H \widehat{u^{*}} 0925 + \partial_y H \widehat{v^{*}} \right) 0926 :label: etastar-nh 0927 0928 .. math:: 0929 \partial_x g H \partial_x \eta^{n+1} 0930 + \partial_y g H \partial_y \eta^{n+1} 0bad585a21 Navi*0931 - \frac{\epsilon_{\rm fs} \eta^{n+1}}{\Delta t^2} 4f2617d475 Jeff*0932 ~ = ~ - \frac{\eta^*}{\Delta t^2} 0933 :label: elliptic-nh 0934 0935 .. math:: 0936 u^{**} = u^{*} - \Delta t g \partial_x \eta^{n+1} 0937 :label: unx-nh 0938 0939 .. math:: 0940 v^{**} = v^{*} - \Delta t g \partial_y \eta^{n+1} 0941 :label: vnx-nh 0942 0943 .. math:: 0bad585a21 Navi*0944 \partial_{xx} \phi_{\rm nh}^{n+1} + \partial_{yy} \phi_{\rm nh}^{n+1} + 0945 \partial_{rr} \phi_{\rm nh}^{n+1} = 4f2617d475 Jeff*0946 \partial_x u^{**} + \partial_y v^{**} + \partial_r w^{*} 0947 :label: phi-nh 0948 0949 .. math:: 0bad585a21 Navi*0950 u^{n+1} = u^{**} - \Delta t \partial_x \phi_{\rm nh}^{n+1} 4f2617d475 Jeff*0951 :label: un+1-nh 0952 0953 .. math:: 0bad585a21 Navi*0954 v^{n+1} = v^{**} - \Delta t \partial_y \phi_{\rm nh}^{n+1} 4f2617d475 Jeff*0955 :label: vn+1-nh 0956 0957 .. math:: 0958 \partial_r w^{n+1} = - \partial_x u^{n+1} - \partial_y v^{n+1} 0959 :label: wn+1-nh 0960 0961 where the last equation is solved by vertically integrating for 0962 :math:`w^{n+1}`. 0963 0964 Variants on the Free Surface 0965 ============================ 0966 0967 We now describe the various formulations of the free-surface that d4a066fa68 Jean*0968 include non-linear forms, implicit in time using Crank-Nicolson, 4f2617d475 Jeff*** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 970.
0969 explicit and [one day] split-explicit. First, we’ll reiterate the 0970 underlying algorithm but this time using the notation consistent with 0971 the more general vertical coordinate :math:`r`. The elliptic equation 0972 for free-surface coordinate (units of :math:`r`), corresponding to 0973 :eq:`discrete-time-backward-free-surface`, and assuming no 0bad585a21 Navi*0974 non-hydrostatic effects (:math:`\epsilon_{\rm nh} = 0`) is: 4f2617d475 Jeff*0975 0976 .. math:: 0bad585a21 Navi*0977 \epsilon_{\rm fs} {\eta}^{n+1} - 0978 \nabla _h \cdot \Delta t^2 (R_o-R_{\rm fixed}) \nabla _h b_s 4f2617d475 Jeff*0979 {\eta}^{n+1} = {\eta}^* 0980 :label: eq-solve2D 0981 0982 where 0983 0984 .. math:: 0bad585a21 Navi*0985 {\eta}^* = \epsilon_{\rm fs} \: {\eta}^{n} - 0986 \Delta t \nabla _h \cdot \int_{R_{\rm fixed}}^{R_o} \vec{\bf v}^* dr 0987 \: + \: \epsilon_{\rm fw} \Delta t ({\mathcal{P-E}})^{n} 4f2617d475 Jeff*0988 :label: eq-solve2D_rhs 0989 0990 .. admonition:: S/R :filelink:`SOLVE_FOR_PRESSURE <model/src/solve_for_pressure.F>` 0991 :class: note 0992 0993 | :math:`u^*` : :varlink:`gU` ( :filelink:`DYNVARS.h <model/inc/DYNVARS.h>` ) 0994 | :math:`v^*` : :varlink:`gV` ( :filelink:`DYNVARS.h <model/inc/DYNVARS.h>` ) 0995 | :math:`{\eta}^*` : :varlink:`cg2d_b` ( :filelink:`SOLVE_FOR_PRESSURE.h <model/inc/SOLVE_FOR_PRESSURE.h>` ) 0996 | :math:`{\eta}^{n+1}` : :varlink:`etaN` ( :filelink:`DYNVARS.h <model/inc/DYNVARS.h>` ) 0997 0998 0999 Once :math:`{\eta}^{n+1}` has been found, substituting into 1000 :eq:`discrete-time-u`, :eq:`discrete-time-v` yields 1001 :math:`\vec{\bf v}^{n+1}` if the model is hydrostatic 0bad585a21 Navi*1002 (:math:`\epsilon_{\rm nh}=0`): 4f2617d475 Jeff*1003 1004 .. math:: 1005 \vec{\bf v}^{n+1} = \vec{\bf v}^{*} 0bad585a21 Navi*1006 - \Delta t \nabla _h b_s {\eta}^{n+1} 4f2617d475 Jeff*1007 1008 This is known as the correction step. However, when the model is 0bad585a21 Navi*1009 non-hydrostatic (:math:`\epsilon_{\rm nh}=1`) we need an additional step and 1010 an additional equation for :math:`\phi'_{\rm nh}`. This is obtained by 4f2617d475 Jeff*1011 substituting :eq:`discrete-time-u-nh`, :eq:`discrete-time-v-nh` and 1012 :eq:`discrete-time-w-nh` into continuity: 1013 1014 .. math:: 0bad585a21 Navi*1015 [ \nabla _h^2 + \partial_{rr} ] {\phi'_{\rm nh}}^{n+1} ab47de63dc Mart*1016 = \frac{1}{\Delta t} 0bad585a21 Navi*1017 \nabla _h \cdot \vec{\bf v}^{**} + \partial_r \dot{r}^* ab47de63dc Mart*1018 :label: sub-u-v-w-in-cont 4f2617d475 Jeff*1019 1020 where 1021 0bad585a21 Navi*1022 .. math:: \vec{\bf v}^{**} = \vec{\bf v}^* - \Delta t \nabla _h b_s {\eta}^{n+1} 4f2617d475 Jeff*1023 1024 Note that :math:`\eta^{n+1}` is also used to update the second RHS term 1025 :math:`\partial_r \dot{r}^*` since the vertical velocity at the surface 0bad585a21 Navi*1026 (:math:`\dot{r}_{\rm surf}`) is evaluated as 4f2617d475 Jeff*1027 :math:`(\eta^{n+1} - \eta^n) / \Delta t`. 1028 1029 Finally, the horizontal velocities at the new time level are found by: 1030 1031 .. math:: 1032 \vec{\bf v}^{n+1} = \vec{\bf v}^{**} 0bad585a21 Navi*1033 - \epsilon_{\rm nh} \Delta t \nabla _h {\phi'_{\rm nh}}^{n+1} ab47de63dc Mart*1034 :label: v-new-time-lev 4f2617d475 Jeff*1035 1036 and the vertical velocity is found by integrating the continuity 1037 equation vertically. Note that, for the convenience of the restart 1038 procedure, the vertical integration of the continuity equation has been 1039 moved to the beginning of the time step (instead of at the end), without 1040 any consequence on the solution. 1041 1042 .. _correction_step_sr_in-out: 1043 1044 .. admonition:: S/R :filelink:`CORRECTION_STEP <model/src/correction_step.F>` 1045 :class: note 1046 1047 | :math:`{\eta}^{n+1}` : :varlink:`etaN` ( :filelink:`DYNVARS.h <model/inc/DYNVARS.h>` ) 0bad585a21 Navi*1048 | :math:`{\phi}^{n+1}_{\rm nh}` : :varlink:`phi_nh` ( :filelink:`NH_VARS.h <model/inc/NH_VARS.h>` ) 4f2617d475 Jeff*1049 | :math:`u^*` : :varlink:`gU` ( :filelink:`DYNVARS.h <model/inc/DYNVARS.h>` ) 1050 | :math:`v^*` : :varlink:`gV` ( :filelink:`DYNVARS.h <model/inc/DYNVARS.h>` ) 1051 | :math:`u^{n+1}` : :varlink:`uVel` ( :filelink:`DYNVARS.h <model/inc/DYNVARS.h>` ) 1052 | :math:`v^{n+1}` : :varlink:`vVel` ( :filelink:`DYNVARS.h <model/inc/DYNVARS.h>` ) 1053 1054 1055 Regarding the implementation of the surface pressure solver, all 1056 computation are done within the routine 1057 :filelink:`SOLVE_FOR_PRESSURE <model/src/solve_for_pressure.F>` and its 1058 dependent calls. The standard method to solve the 2D elliptic problem 1059 :eq:`eq-solve2D` uses the conjugate gradient method 1060 (routine :filelink:`CG2D <model/src/cg2d.F>`); the 1061 solver matrix and conjugate gradient operator are only function of the 1062 discretized domain and are therefore evaluated separately, before the 1063 time iteration loop, within :filelink:`INI_CG2D <model/src/ini_cg2d.F>`. The computation of the RHS 1064 :math:`\eta^*` is partly done in :filelink:`CALC_DIV_GHAT <model/src/calc_div_ghat.F>` and in 1065 :filelink:`SOLVE\_FOR\_PRESSURE <model/src/solve_for_pressure.F>`. 1066 1067 The same method is applied for the non hydrostatic part, using a 1068 conjugate gradient 3D solver (:filelink:`CG3D <model/src/cg3d.F>`) that is initialized in 1069 :filelink:`INI_CG3D <model/src/ini_cg3d.F>`. The RHS terms of 2D and 3D problems are computed together 1070 at the same point in the code. 1071 1072 .. toctree:: 1073 crank-nicol.rst 1074 nonlinear-freesurf.rst 1075 1c87316fba Jeff*1076 .. _spatial_discret_dyn_eq: 4f2617d475 Jeff*1077 1078 Spatial discretization of the dynamical equations 1079 ================================================= 1080 1081 Spatial discretization is carried out using the finite volume method. 1082 This amounts to a grid-point method (namely second-order centered finite 1083 difference) in the fluid interior but allows boundaries to intersect a 1084 regular grid allowing a more accurate representation of the position of 1085 the boundary. We treat the horizontal and vertical directions as 1086 separable and differently. 1087 1088 .. toctree:: 1089 finitevol-meth.rst 1090 c-grid.rst 1091 horiz-grid.rst 1092 vert-grid.rst ab47de63dc Mart*1093 4f2617d475 Jeff*1094 Continuity and horizontal pressure gradient term 1095 ================================================= 1096 1097** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 1099.
1098 The core algorithm is based on the “C grid” discretization of the 1099 continuity equation which can be summarized as: 1100 1101 .. math:: 0bad585a21 Navi*1102 \partial_t u + \frac{1}{\Delta x_c} \delta_i \left. \frac{ \partial \Phi}{\partial r}\right|_{s} \eta + \frac{\epsilon_{\rm nh}}{\Delta x_c} \delta_i \Phi_{\rm nh}' = G_u - \frac{1}{\Delta x_c} \delta_i \Phi_h' 4f2617d475 Jeff*1103 :label: discrete-momu 1104 1105 .. math:: 0bad585a21 Navi*1106 \partial_t v + \frac{1}{\Delta y_c} \delta_j \left. \frac{ \partial \Phi}{\partial r}\right|_{s} \eta + \frac{\epsilon_{\rm nh}}{\Delta y_c} \delta_j \Phi_{\rm nh}' = G_v - \frac{1}{\Delta y_c} \delta_j \Phi_h' 4f2617d475 Jeff*1107 :label: discrete-momv 1108 1109 .. math:: 0bad585a21 Navi*1110 \epsilon_{\rm nh} \left( \partial_t w + \frac{1}{\Delta r_c} \delta_k \Phi_{\rm nh}' \right) = \epsilon_{\rm nh} G_w + \overline{b}^k - \frac{1}{\Delta r_c} \delta_k \Phi_{h}' 4f2617d475 Jeff*1111 :label: discrete-momw 1112 1113 .. math:: 1114 \delta_i \Delta y_g \Delta r_f h_w u + 1115 \delta_j \Delta x_g \Delta r_f h_s v + 94151a9b18 Jeff*1116 \delta_k {\cal A}_c w = {\cal A}_c \delta_k (\mathcal{P-E})_{r=0} 4f2617d475 Jeff*1117 :label: discrete-continuity 1118 1119 where the continuity equation has been most naturally discretized by 1120 staggering the three components of velocity as shown in 1121 :numref:`cgrid3d`. The grid lengths :math:`\Delta x_c` and 1122 :math:`\Delta y_c` are the lengths between tracer points (cell centers). 1123 The grid lengths :math:`\Delta x_g`, :math:`\Delta y_g` are the grid 1124 lengths between cell corners. :math:`\Delta r_f` and :math:`\Delta r_c` 1125 are the distance (in units of :math:`r`) between level interfaces 1126 (w-level) and level centers (tracer level). The surface area presented 1127 in the vertical is denoted :math:`{\cal 1128 A}_c`. The factors :math:`h_w` and :math:`h_s` are non-dimensional 1129 fractions (between 0 and 1) that represent the fraction cell depth that** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 1131.
1130 is “open” for fluid flow. 1131 1132 The last equation, the discrete continuity equation, can be summed in 1133 the vertical to yield the free-surface equation: 1134 1135 .. math:: 1136 {\cal A}_c \partial_t \eta + \delta_i \sum_k \Delta y_g \Delta r_f h_w 1137 u + \delta_j \sum_k \Delta x_g \Delta r_f h_s v = {\cal 94151a9b18 Jeff*1138 A}_c(\mathcal{P-E})_{r=0} 4f2617d475 Jeff*1139 :label: discrete-freesurface 1140 94151a9b18 Jeff*1141 The source term :math:`\mathcal{P-E}` on the rhs of continuity accounts for the 4f2617d475 Jeff*1142 local addition of volume due to excess precipitation and run-off over 1143 evaporation and only enters the top-level of the ocean model. 1144 1145 Hydrostatic balance 1146 =================== 1147 1148 The vertical momentum equation has the hydrostatic or quasi-hydrostatic 1149 balance on the right hand side. This discretization guarantees that the 1150 conversion of potential to kinetic energy as derived from the buoyancy 1151 equation exactly matches the form derived from the pressure gradient 1152 terms when forming the kinetic energy equation. 1153 1154 In the ocean, using z-coordinates, the hydrostatic balance terms are 1155 discretized: 1156 1157 .. math:: 0bad585a21 Navi*1158 \epsilon_{\rm nh} \partial_t w 4f2617d475 Jeff*1159 + g \overline{\rho'}^k + \frac{1}{\Delta z} \delta_k \Phi_h' = \ldots 1160 :label: discrete_hydro_ocean 1161 1162 In the atmosphere, using p-coordinates, hydrostatic balance is 1163 discretized: 1164 1165 .. math:: 1166 \overline{\theta'}^k + \frac{1}{\Delta \Pi} \delta_k \Phi_h' = 0 1167 :label: discrete_hydro_atmos 1168 1169 where :math:`\Delta \Pi` is the difference in Exner function between 1170 the pressure points. The non-hydrostatic equations are not available in 1171 the atmosphere. 1172 1173 The difference in approach between ocean and atmosphere occurs because 1174 of the direct use of the ideal gas equation in forming the potential 1175 energy conversion term :math:`\alpha \omega`. Because of the different 1176 representation of hydrostatic balance between 1177 ocean and atmosphere there is no elegant way to represent both systems 1178 using an arbitrary coordinate. 1179 1180 The integration for hydrostatic pressure is made in the positive 1181 :math:`r` direction (increasing k-index). For the ocean, this is from 1182 the free-surface down and for the atmosphere this is from the ground up. 1183 1184 The calculations are made in the subroutine :filelink:`CALC_PHI_HYD <model/src/calc_phi_hyd.F>`. Inside 1185 this routine, one of other of the atmospheric/oceanic form is selected 1186 based on the string variable :varlink:`buoyancyRelation`. 1187 68aadaa6bd Phob*1188 .. _flux-form_momentum_equations: 4f2617d475 Jeff*1189 1190 Flux-form momentum equations 1191 ============================ 1192 1193 1194 The original finite volume model was based on the Eulerian flux form 1195 momentum equations. This is the default though the vector invariant form 1196 is optionally available (and recommended in some cases). 1197** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 1199.
1198 The “G’s” (our colloquial name for all terms on rhs!) are broken into 1199 the various advective, Coriolis, horizontal dissipation, vertical 1200 dissipation and metric forces: 1201 1202 .. math:: 0bad585a21 Navi*1203 G_u = G_u^{\rm adv} + G_u^{\rm Cor} + G_u^{h- \rm diss} + G_u^{v- \rm diss} + 1204 G_u^{\rm metric} + G_u^{\rm nh-metric} 4f2617d475 Jeff*1205 :label: gsplit_momu 1206 1207 .. math:: 0bad585a21 Navi*1208 G_v = G_v^{\rm adv} + G_v^{\rm Cor} + G_v^{h- \rm diss} + G_v^{v- \rm diss} + 1209 G_v^{\rm metric} + G_v^{\rm nh-metric} 4f2617d475 Jeff*1210 :label: gsplit_momv 1211 1212 .. math:: 0bad585a21 Navi*1213 G_w = G_w^{\rm adv} + G_w^{\rm Cor} + G_w^{h- \rm diss} + G_w^{v- \rm diss} + 1214 G_w^{\rm metric} + G_w^{\rm nh-metric} 4f2617d475 Jeff*1215 :label: gsplit_momw 1216 0bad585a21 Navi*1217 In the hydrostatic limit, :math:`G_w=0` and :math:`\epsilon_{\rm nh}=0`, 4f2617d475 Jeff*1218 reducing the vertical momentum to hydrostatic balance. 1219 1220 These terms are calculated in routines called from subroutine 1221 :filelink:`MOM_FLUXFORM <pkg/mom_fluxform/mom_fluxform.F>` and collected into the global arrays :varlink:`gU`, :varlink:`gV`, and 1222 :varlink:`gW`. 1223 1224 .. admonition:: S/R :filelink:`MOM_FLUXFORM <pkg/mom_fluxform/mom_fluxform.F>` 1225 :class: note 1226 1227 | :math:`G_u` : :varlink:`gU` ( :filelink:`DYNVARS.h <model/inc/DYNVARS.h>` ) 1228 | :math:`G_v` : :varlink:`gV` ( :filelink:`DYNVARS.h <model/inc/DYNVARS.h>` ) 1229 | :math:`G_w` : :varlink:`gW` ( :filelink:`NH_VARS.h <model/inc/NH_VARS.h>` ) ab47de63dc Mart*1230 4f2617d475 Jeff*1231 1232 Advection of momentum 1233 --------------------- 1234 1235 The advective operator is second order accurate in space: 1236 1237 .. math:: 0bad585a21 Navi*1238 {\cal A}_w \Delta r_f h_w G_u^{\rm adv} = 4f2617d475 Jeff*1239 \delta_i \overline{ U }^i \overline{ u }^i 1240 + \delta_j \overline{ V }^i \overline{ u }^j 1241 + \delta_k \overline{ W }^i \overline{ u }^k 1242 :label: discrete-momadvu 1243 1244 .. math:: 0bad585a21 Navi*1245 {\cal A}_s \Delta r_f h_s G_v^{\rm adv} = 4f2617d475 Jeff*1246 \delta_i \overline{ U }^j \overline{ v }^i 1247 + \delta_j \overline{ V }^j \overline{ v }^j ab47de63dc Mart*1248 + \delta_k \overline{ W }^j \overline{ v }^k 4f2617d475 Jeff*1249 :label: discrete-momadvv 1250 1251 .. math:: 0bad585a21 Navi*1252 {\cal A}_c \Delta r_c G_w^{\rm adv} = 4f2617d475 Jeff*1253 \delta_i \overline{ U }^k \overline{ w }^i 1254 + \delta_j \overline{ V }^k \overline{ w }^j 1255 + \delta_k \overline{ W }^k \overline{ w }^k 1256 :label: discrete-momadvw 1257 1258 and because of the flux form does not contribute to the global budget 1259 of linear momentum. The quantities :math:`U`, :math:`V` and :math:`W` 1260 are volume fluxes defined: 1261 1262 .. math:: 1263 U = \Delta y_g \Delta r_f h_w u 1264 :label: utrans 1265 1266 .. math:: 1267 V = \Delta x_g \Delta r_f h_s v 1268 :label: vtrans 1269 1270 .. math:: 1271 W = {\cal A}_c w 1272 :label: rtrans 1273 1274 The advection of momentum takes the same form as the advection of 1275 tracers but by a translated advective flow. Consequently, the 1276 conservation of second moments, derived for tracers later, applies to 1277 :math:`u^2` and :math:`v^2` and :math:`w^2` so that advection of 1278 momentum correctly conserves kinetic energy. 1279 1280 .. admonition:: S/R :filelink:`MOM_U_ADV_UU <pkg/mom_fluxform/mom_u_adv_uu.F>`, :filelink:`MOM_U_ADV_VU <pkg/mom_fluxform/mom_u_adv_vu.F>`, :filelink:`MOM_U_ADV_WU <pkg/mom_fluxform/mom_u_adv_wu.F>` 1281 :class: note 1282 1283 | :math:`uu, vu, wu` : :varlink:`fZon`, :varlink:`fMer`, :varlink:`fVerUkp` ( local to :filelink:`MOM_FLUXFORM.F <pkg/mom_fluxform/mom_fluxform.F>` ) 1284 1285 .. admonition:: S/R :filelink:`MOM_V_ADV_UV <pkg/mom_fluxform/mom_v_adv_uv.F>`, :filelink:`MOM_V_ADV_VV <pkg/mom_fluxform/mom_v_adv_vv.F>`, :filelink:`MOM_V_ADV_WV <pkg/mom_fluxform/mom_v_adv_wv.F>` 1286 :class: note 1287 1288 | :math:`uv, vv, wv` : :varlink:`fZon`, :varlink:`fMer`, :varlink:`fVerVkp` ( local to :filelink:`MOM_FLUXFORM.F <pkg/mom_fluxform/mom_fluxform.F>` ) 1289 9c8516d9da Jeff*1290 .. _fluxform_cor_terms: 1291 4f2617d475 Jeff*1292 Coriolis terms 1293 -------------- 1294** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 1296.
1295 The “pure C grid” Coriolis terms (i.e. in absence of C-D scheme) are 1296 discretized: 1297 1298 .. math:: 0bad585a21 Navi*1299 {\cal A}_w \Delta r_f h_w G_u^{\rm Cor} = 4f2617d475 Jeff*1300 \overline{ f {\cal A}_c \Delta r_f h_c \overline{ v }^j }^i 0bad585a21 Navi*1301 - \epsilon_{\rm nh} \overline{ f' {\cal A}_c \Delta r_f h_c \overline{ w }^k }^i 4f2617d475 Jeff*1302 :label: cdscheme_gu 1303 1304 .. math:: 0bad585a21 Navi*1305 {\cal A}_s \Delta r_f h_s G_v^{\rm Cor} = 4f2617d475 Jeff*1306 - \overline{ f {\cal A}_c \Delta r_f h_c \overline{ u }^i }^j 1307 :label: cdscheme_gv 1308 1309 .. math:: 0bad585a21 Navi*1310 {\cal A}_c \Delta r_c G_w^{\rm Cor} = 1311 \epsilon_{\rm nh} \overline{ f' {\cal A}_c \Delta r_f h_c \overline{ u }^i }^k 4f2617d475 Jeff*1312 :label: cdscheme_gw 1313 1314 where the Coriolis parameters :math:`f` and :math:`f'` are defined: 1315 1316 .. math:: 1317 1318 \begin{aligned} 0bad585a21 Navi*1319 f & = 2 \Omega \sin{\varphi} \\ 1320 f' & = 2 \Omega \cos{\varphi}\end{aligned} 4f2617d475 Jeff*1321 1322 where :math:`\varphi` is geographic latitude when using spherical 1323 geometry, otherwise the :math:`\beta`-plane definition is used: 1324 1325 .. math:: 1326 1327 \begin{aligned} 0bad585a21 Navi*1328 f & = f_o + \beta y \\ 1329 f' & = 0\end{aligned} 4f2617d475 Jeff*1330 1331 This discretization globally conserves kinetic energy. It should be 1332 noted that despite the use of this discretization in former 1333 publications, all calculations to date have used the following different 1334 discretization: 1335 1336 .. math:: 0bad585a21 Navi*1337 G_u^{\rm Cor} = f_u \overline{ v }^{ji} 1338 - \epsilon_{\rm nh} f_u' \overline{ w }^{ik} 4f2617d475 Jeff*1339 :label: gu_cor 1340 1341 .. math:: 0bad585a21 Navi*1342 G_v^{\rm Cor} = - f_v \overline{ u }^{ij} 4f2617d475 Jeff*1343 :label: gv_cor 1344 1345 .. math:: 0bad585a21 Navi*1346 G_w^{\rm Cor} = \epsilon_{\rm nh} f_w' \overline{ u }^{ik} 4f2617d475 Jeff*1347 :label: gw_cor 1348 1349 where the subscripts on :math:`f` and :math:`f'` indicate evaluation of 1350 the Coriolis parameters at the appropriate points in space. The above 1351 discretization does *not* conserve anything, especially energy, but for 1352 historical reasons is the default for the code. A flag controls this 7843dde2de jm-c 1353 discretization: set run-time integer :varlink:`selectCoriScheme` to two (=2) 1354 (which otherwise defaults to zero) 1355 to select the energy-conserving conserving form :eq:`cdscheme_gu`, :eq:`cdscheme_gv`, and :eq:`cdscheme_gw` above. 1356 4f2617d475 Jeff*1357 1358 .. admonition:: S/R :filelink:`CD_CODE_SCHEME <pkg/cd_code/cd_code_scheme.F>`, :filelink:`MOM_U_CORIOLIS <pkg/mom_fluxform/mom_u_coriolis.F>`, :filelink:`MOM_V_CORIOLIS <pkg/mom_fluxform/mom_v_coriolis.F>` 1359 :class: note 1360 0bad585a21 Navi*1361 | :math:`G_u^{\rm Cor}, G_v^{\rm Cor}` : :varlink:`cF` ( local to :filelink:`MOM_FLUXFORM.F <pkg/mom_fluxform/mom_fluxform.F>` ) 4f2617d475 Jeff*1362 1363 Curvature metric terms 1364 ---------------------- 1365 1366 The most commonly used coordinate system on the sphere is the geographic 1367 system :math:`(\lambda,\varphi)`. The curvilinear nature of these** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 1369.
1368 coordinates on the sphere lead to some “metric” terms in the component 1369 momentum equations. Under the thin-atmosphere and hydrostatic 1370 approximations these terms are discretized: 1371 1372 .. math:: 0bad585a21 Navi*1373 {\cal A}_w \Delta r_f h_w G_u^{\rm metric} = 4f2617d475 Jeff*1374 \overline{ \frac{ \overline{u}^i }{a} \tan{\varphi} {\cal A}_c \Delta r_f h_c \overline{ v }^j }^i 1375 :label: gu_metric ab47de63dc Mart*1376 4f2617d475 Jeff*1377 .. math:: 0bad585a21 Navi*1378 {\cal A}_s \Delta r_f h_s G_v^{\rm metric} = 4f2617d475 Jeff*1379 - \overline{ \frac{ \overline{u}^i }{a} \tan{\varphi} {\cal A}_c \Delta r_f h_c \overline{ u }^i }^j \\ 1380 :label: gv_metric ab47de63dc Mart*1381 4f2617d475 Jeff*1382 .. math:: 0bad585a21 Navi*1383 G_w^{\rm metric} = 0 4f2617d475 Jeff*1384 :label: gw_metric ab47de63dc Mart*1385 4f2617d475 Jeff*1386 where :math:`a` is the radius of the planet (sphericity is assumed) or 1387 the radial distance of the particle (i.e. a function of height). It is 1388 easy to see that this discretization satisfies all the properties of the 1389 discrete Coriolis terms since the metric factor :math:`\frac{u}{a} 1390 \tan{\varphi}` can be viewed as a modification of the vertical Coriolis 1391 parameter: :math:`f \rightarrow f+\frac{u}{a} \tan{\varphi}`. 1392 1393 However, as for the Coriolis terms, a non-energy conserving form has 1394 exclusively been used to date: 1395 1396 .. math:: 1397 1398 \begin{aligned} 0bad585a21 Navi*1399 G_u^{\rm metric} & = \frac{u \overline{v}^{ij} }{a} \tan{\varphi} \\ 1400 G_v^{\rm metric} & = \frac{ \overline{u}^{ij} \overline{u}^{ij}}{a} \tan{\varphi}\end{aligned} 4f2617d475 Jeff*1401 1402 where :math:`\tan{\varphi}` is evaluated at the :math:`u` and :math:`v` 1403 points respectively. 1404 1405 .. admonition:: S/R :filelink:`MOM_U_METRIC_SPHERE <pkg/mom_fluxform/mom_u_metric_sphere.F>`, :filelink:`MOM_V_METRIC_SPHERE <pkg/mom_fluxform/mom_v_metric_sphere.F>` 1406 :class: note 1407 0bad585a21 Navi*1408 | :math:`G_u^{\rm metric}, G_v^{\rm metric}` : :varlink:`mT` ( local to :filelink:`MOM_FLUXFORM.F <pkg/mom_fluxform/mom_fluxform.F>` ) 4f2617d475 Jeff*1409 dcaaa42497 Jeff*1410 .. _non_hyd_metric_terms: 4f2617d475 Jeff*1411 1412 Non-hydrostatic metric terms 1413 ---------------------------- 1414 1415 For the non-hydrostatic equations, dropping the thin-atmosphere 1416 approximation re-introduces metric terms involving :math:`w` which are 1417 required to conserve angular momentum: 1418 1419 .. math:: 0bad585a21 Navi*1420 {\cal A}_w \Delta r_f h_w G_u^{\rm metric} = 4f2617d475 Jeff*1421 - \overline{ \frac{ \overline{u}^i \overline{w}^k }{a} {\cal A}_c \Delta r_f h_c }^i 1422 :label: gu_metricnh ab47de63dc Mart*1423 4f2617d475 Jeff*1424 .. math:: 0bad585a21 Navi*1425 {\cal A}_s \Delta r_f h_s G_v^{\rm metric} = 4f2617d475 Jeff*1426 - \overline{ \frac{ \overline{v}^j \overline{w}^k }{a} {\cal A}_c \Delta r_f h_c}^j 1427 :label: gv_metricnh ab47de63dc Mart*1428 4f2617d475 Jeff*1429 .. math:: 0bad585a21 Navi*1430 {\cal A}_c \Delta r_c G_w^{\rm metric} = 4f2617d475 Jeff*1431 \overline{ \frac{ {\overline{u}^i}^2 + {\overline{v}^j}^2}{a} {\cal A}_c \Delta r_f h_c }^k 1432 :label: wv_metricnh 1433 1434 Because we are always consistent, even if consistently wrong, we have, 1435 in the past, used a different discretization in the model which is: 1436 1437 .. math:: 1438 1439 \begin{aligned} ab47de63dc Mart*1440 G_u^{\rm metric} & = 4f2617d475 Jeff*1441 - \frac{u}{a} \overline{w}^{ik} \\ ab47de63dc Mart*1442 G_v^{\rm metric} & = 4f2617d475 Jeff*1443 - \frac{v}{a} \overline{w}^{jk} \\ ab47de63dc Mart*1444 G_w^{\rm metric} & = 4f2617d475 Jeff*1445 \frac{1}{a} ( {\overline{u}^{ik}}^2 + {\overline{v}^{jk}}^2 )\end{aligned} 1446 1447 .. admonition:: S/R :filelink:`MOM_U_METRIC_NH <pkg/mom_common/mom_u_metric_nh.F>`, :filelink:`MOM_V_METRIC_NH <pkg/mom_common/mom_v_metric_nh.F>` 1448 :class: note 1449 0bad585a21 Navi*1450 | :math:`G_u^{\rm metric}, G_v^{\rm metric}` : :varlink:`mT` ( local to :filelink:`MOM_FLUXFORM.F <pkg/mom_fluxform/mom_fluxform.F>` ) 4f2617d475 Jeff*1451 9c8516d9da Jeff*1452 .. _fluxform_lat_dissip: 4f2617d475 Jeff*1453 1454 Lateral dissipation 1455 ------------------- 1456 1457 Historically, we have represented the SGS Reynolds stresses as simply 1458 down gradient momentum fluxes, ignoring constraints on the stress tensor 1459 such as symmetry. 1460 1461 .. math:: 0bad585a21 Navi*1462 {\cal A}_w \Delta r_f h_w G_u^{h- \rm diss} = 4f2617d475 Jeff*1463 \delta_i \Delta y_f \Delta r_f h_c \tau_{11} 1464 + \delta_j \Delta x_v \Delta r_f h_\zeta \tau_{12} 1465 :label: gu_h-diss 1466 1467 .. math:: 0bad585a21 Navi*1468 {\cal A}_s \Delta r_f h_s G_v^{h- \rm diss} = 4f2617d475 Jeff*1469 \delta_i \Delta y_u \Delta r_f h_\zeta \tau_{21} 1470 + \delta_j \Delta x_f \Delta r_f h_c \tau_{22} 1471 :label: gv_h-diss 1472 1473 1474 The lateral viscous stresses are discretized: 1475 1476 .. math:: 1477 \tau_{11} = A_h c_{11\Delta}(\varphi) \frac{1}{\Delta x_f} \delta_i u 1478 -A_4 c_{11\Delta^2}(\varphi) \frac{1}{\Delta x_f} \delta_i \nabla^2 u 1479 :label: tau11 1480 1481 .. math:: 1482 \tau_{12} = A_h c_{12\Delta}(\varphi) \frac{1}{\Delta y_u} \delta_j u 1483 -A_4 c_{12\Delta^2}(\varphi)\frac{1}{\Delta y_u} \delta_j \nabla^2 u 1484 :label: tau12 1485 1486 .. math:: 1487 \tau_{21} = A_h c_{21\Delta}(\varphi) \frac{1}{\Delta x_v} \delta_i v 1488 -A_4 c_{21\Delta^2}(\varphi) \frac{1}{\Delta x_v} \delta_i \nabla^2 v 1489 :label: tau21 1490 1491 .. math:: 1492 \tau_{22} = A_h c_{22\Delta}(\varphi) \frac{1}{\Delta y_f} \delta_j v 1493 -A_4 c_{22\Delta^2}(\varphi) \frac{1}{\Delta y_f} \delta_j \nabla^2 v 1494 :label: tau22 1495 025afa4065 Jeff*1496 where the non-dimensional factors :math:`c_{lm\Delta^n}(\varphi), \{l,m,n\} \in \{1,2\}`** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 1498.
1497 define the “cosine” scaling with latitude which can be applied 4f2617d475 Jeff*1498 in various ad-hoc ways. For instance, :math:`c_{11\Delta} = 1499 c_{21\Delta} = (\cos{\varphi})^{3/2}`, 1500 :math:`c_{12\Delta}=c_{22\Delta}=1` would represent the anisotropic** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 1502.
1501 cosine scaling typically used on the “lat-lon” grid for Laplacian 1502 viscosity. 1503 1504 It should be noted that despite the ad-hoc nature of the scaling, some 1505 scaling must be done since on a lat-lon grid the converging meridians 1506 make it very unlikely that a stable viscosity parameter exists across 1507 the entire model domain. 1508 1509 The Laplacian viscosity coefficient, :math:`A_h` (:varlink:`viscAh`), has units 1510 of :math:`m^2 s^{-1}`. The bi-harmonic viscosity coefficient, 1511 :math:`A_4` (:varlink:`viscA4`), has units of :math:`m^4 s^{-1}`. 1512 1513 .. admonition:: S/R :filelink:`MOM_U_XVISCFLUX <pkg/mom_fluxform/mom_u_xviscflux.F>`, :filelink:`MOM_U_YVISCFLUX <pkg/mom_fluxform/mom_u_yviscflux.F>` 1514 :class: note 1515 1516 | :math:`\tau_{11}, \tau_{12}` : :varlink:`vF`, :varlink:`v4F` ( local to :filelink:`MOM_FLUXFORM.F <pkg/mom_fluxform/mom_fluxform.F>` ) 1517 1518 .. admonition:: S/R :filelink:`MOM_V_XVISCFLUX <pkg/mom_fluxform/mom_v_xviscflux.F>`, :filelink:`MOM_V_YVISCFLUX <pkg/mom_fluxform/mom_v_yviscflux.F>` 1519 :class: note 1520 1521 | :math:`\tau_{21}, \tau_{22}` : :varlink:`vF`, :varlink:`v4F` ( local to :filelink:`MOM_FLUXFORM.F <pkg/mom_fluxform/mom_fluxform.F>` ) 1522 1523 Two types of lateral boundary condition exist for the lateral viscous 1524 terms, no-slip and free-slip. 1525 1526 The free-slip condition is most convenient to code since it is 1527 equivalent to zero-stress on boundaries. Simple masking of the stress 1528 components sets them to zero. The fractional open stress is properly 1529 handled using the lopped cells. 1530 1531 The no-slip condition defines the normal gradient of a tangential flow 1532 such that the flow is zero on the boundary. Rather than modify the** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 1534.
1533 stresses by using complicated functions of the masks and “ghost” points 1534 (see Adcroft and Marshall (1998) :cite:`adcroft:98`) we add the boundary stresses as an 1535 additional source term in cells next to solid boundaries. This has the** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 1537.
1536 advantage of being able to cope with “thin walls” and also makes the 1537 interior stress calculation (code) independent of the boundary** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 1539.
1538 conditions. The “body” force takes the form: 1539 1540 .. math:: 0bad585a21 Navi*1541 G_u^{\rm side-drag} = 4f2617d475 Jeff*1542 \frac{4}{\Delta z_f} \overline{ (1-h_\zeta) \frac{\Delta x_v}{\Delta y_u} }^j 1543 \left( A_h c_{12\Delta}(\varphi) u - A_4 c_{12\Delta^2}(\varphi) \nabla^2 u \right) 1544 :label: gu_sidedrag 1545 1546 .. math:: 0bad585a21 Navi*1547 G_v^{\rm side-drag} = 4f2617d475 Jeff*1548 \frac{4}{\Delta z_f} \overline{ (1-h_\zeta) \frac{\Delta y_u}{\Delta x_v} }^i 1549 \left( A_h c_{21\Delta}(\varphi) v - A_4 c_{21\Delta^2}(\varphi) \nabla^2 v \right) 1550 :label: gv_sidedrag 1551 1552 In fact, the above discretization is not quite complete because it 1553 assumes that the bathymetry at velocity points is deeper than at 1554 neighboring vorticity points, e.g. :math:`1-h_w < 1-h_\zeta` 1555 1556 .. admonition:: S/R :filelink:`MOM_U_SIDEDRAG <pkg/mom_common/mom_u_sidedrag.F>`, :filelink:`MOM_V_SIDEDRAG <pkg/mom_common/mom_v_sidedrag.F>` 1557 :class: note 1558 0bad585a21 Navi*1559 | :math:`G_u^{\rm side-drag}, G_v^{\rm side-drag}` : :varlink:`vF` ( local to :filelink:`MOM_FLUXFORM.F <pkg/mom_fluxform/mom_fluxform.F>` ) 4f2617d475 Jeff*1560 1561 Vertical dissipation 1562 -------------------- 1563 1564 Vertical viscosity terms are discretized with only partial adherence to 1565 the variable grid lengths introduced by the finite volume formulation. 1566 This reduces the formal accuracy of these terms to just first order but 1567 only next to boundaries; exactly where other terms appear such as linear 1568 and quadratic bottom drag. 1569 1570 .. math:: 0bad585a21 Navi*1571 G_u^{v- \rm diss} = 4f2617d475 Jeff*1572 \frac{1}{\Delta r_f h_w} \delta_k \tau_{13} 1573 :label: gu_u-diss 1574 1575 .. math:: 0bad585a21 Navi*1576 G_v^{v- \rm diss} = 4f2617d475 Jeff*1577 \frac{1}{\Delta r_f h_s} \delta_k \tau_{23} 1578 :label: gv_v-diss 1579 1580 .. math:: 0bad585a21 Navi*1581 G_w^{v- \rm diss} = \epsilon_{\rm nh} 4f2617d475 Jeff*1582 \frac{1}{\Delta r_f h_d} \delta_k \tau_{33} 1583 :label: gw_w-diss 1584 1585 represents the general discrete form of the vertical dissipation terms. 1586 1587 In the interior the vertical stresses are discretized: 1588 1589 .. math:: 1590 1591 \begin{aligned} 0bad585a21 Navi*1592 \tau_{13} & = A_v \frac{1}{\Delta r_c} \delta_k u \\ 1593 \tau_{23} & = A_v \frac{1}{\Delta r_c} \delta_k v \\ 1594 \tau_{33} & = A_v \frac{1}{\Delta r_f} \delta_k w\end{aligned} 4f2617d475 Jeff*1595 1596 It should be noted that in the non-hydrostatic form, the stress tensor 1597 is even less consistent than for the hydrostatic (see Wajsowicz (1993) 1598 :cite:`wajsowicz:93`). It is well known how to do this 1599 properly (see Griffies and Hallberg (2000) :cite:`griffies:00`) and is on the list of** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 1601.
1600 to-do’s. 1601 1602 .. admonition:: S/R :filelink:`MOM_U_RVISCFLUX <pkg/mom_common/mom_u_rviscflux.F>`, :filelink:`MOM_V_RVISCFLUX <pkg/mom_common/mom_v_rviscflux.F>` 1603 :class: note 1604 1605 | :math:`\tau_{13}` : :varlink:`fVrUp`, :varlink:`fVrDw` ( local to :filelink:`MOM_FLUXFORM.F <pkg/mom_fluxform/mom_fluxform.F>` ) 1606 | :math:`\tau_{23}` : :varlink:`fVrUp`, :varlink:`fVrDw` ( local to :filelink:`MOM_FLUXFORM.F <pkg/mom_fluxform/mom_fluxform.F>` ) 1607 1608 As for the lateral viscous terms, the free-slip condition is equivalent 1609 to simply setting the stress to zero on boundaries. The no-slip ab47de63dc Mart*1610 condition is implemented as an additional term acting in conjunction with the 4f2617d475 Jeff*1611 interior and free-slip stresses. Bottom drag represents additional 1612 friction, in addition to that imposed by the no-slip condition at the 1613 bottom. The drag is cast as a stress expressed as a linear or quadratic 1614 function of the mean flow in the layer above the topography: 1615 1616 .. math:: 0bad585a21 Navi*1617 \tau_{13}^{\rm bottom-drag} = 4f2617d475 Jeff*1618 \left( 1619 2 A_v \frac{1}{\Delta r_c} 1620 + r_b 0bad585a21 Navi*1621 + C_d \sqrt{ \overline{2 \mathrm{KE}}^i } 4f2617d475 Jeff*1622 \right) u 1623 :label: tau13 1624 1625 .. math:: 0bad585a21 Navi*1626 \tau_{23}^{\rm bottom-drag} = 4f2617d475 Jeff*1627 \left( 1628 2 A_v \frac{1}{\Delta r_c} 1629 + r_b 0bad585a21 Navi*1630 + C_d \sqrt{ \overline{2 \mathrm{KE}}^j } 4f2617d475 Jeff*1631 \right) v 1632 :label: tau23 1633 1634 where these terms are only evaluated immediately above topography. ab47de63dc Mart*1635 :math:`r_b` (:varlink:`bottomDragLinear`) has units of m s\ :sup:`-1` and a 1636 typical value of the order 0.0002 m s\ :sup:`-1`. :math:`C_d` 4f2617d475 Jeff*1637 (:varlink:`bottomDragQuadratic`) is dimensionless with typical values in the** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 1639.
1638 range 0.001–0.003. 1639 ab47de63dc Mart*1640 After defining :varlink:`ALLOW_BOTTOMDRAG_ROUGHNESS` in 1641 :filelink:`MOM_COMMON_OPTIONS.h <pkg/mom_common/MOM_COMMON_OPTIONS.h>`, 1642 the quadratic drag coefficient can be 1643 computed by the logarithmic law of the wall: 1644 1645 .. math:: 1646 u(z) = \left(\frac{\tau}{\rho}\right)^{\frac{1}{2}}\frac{1}{0.4} 1647 \ln{\frac{z+z_r}{z_r}} 1648 :label: logLawWallu 1649 1650 where :math:`z_r` is the roughness length (runtime parameter 1651 :varlink:`zRoughBot`). Here, :math:`z` is the height from the seafloor and 1652 :math:`\tau` is the bottom stress (and stress in the log-layer). The velocity 1653 is computed at the center of the bottom cell :math:`z_b=\frac{1}{2}\Delta r_f 1654 h_w`, so stress on the bottom cell is :math:`\tau / \rho = C_d u_b^2`, where 1655 :math:`u_b = u(z_b)` is the bottom cell velocity and: 1656 1657 .. math:: 1658 C_d = \left(\frac{0.4}{ 1659 \ln{\frac{\frac{1}{2} \Delta r_f h_w + z_{r} }{z_{r}}}}\right)^2 1660 :label: logLawWallcd 1661 1662 This formulation assumes that the bottommost cell is sufficiently thin that it 1663 is in a constant-stress log layer. A value of :varlink:`zRoughBot` of 0.01 m 1664 yields a quadratic drag coefficient of 0.0022 for :math:`\Delta r_f =` 100 m, 1665 or a quadratic drag coefficient of 0.0041 for :math:`\Delta r_f =` 10 m. 1666 1667 For :math:`z_r = 0`, :math:`C_d` defaults to the value of 1668 :varlink:`bottomDragQuadratic`. 1669 1670 .. admonition:: S/R :filelink:`MOM_U_BOTTOMDRAG 1671 <pkg/mom_common/mom_u_bottomdrag.F>`, 1672 :filelink:`MOM_V_BOTTOMDRAG 1673 <pkg/mom_common/mom_v_bottomdrag.F>` 4f2617d475 Jeff*1674 :class: note 1675 0bad585a21 Navi*1676 | :math:`\tau_{13}^{\rm bottom-drag} / \Delta r_f , \tau_{23}^{\rm bottom-drag} / \Delta r_f` : :varlink:`vF` ( local to :filelink:`MOM_FLUXFORM.F <pkg/mom_fluxform/mom_fluxform.F>` ) 4f2617d475 Jeff*1677 1678 Derivation of discrete energy conservation 1679 ------------------------------------------ 1680 1681 These discrete equations conserve kinetic plus potential energy using 1682 the following definitions: 1683 1684 .. math:: 0bad585a21 Navi*1685 \mathrm{KE} = \frac{1}{2} \left( \overline{ u^2 }^i + \overline{ v^2 }^j + 1686 \epsilon_{\rm nh} \overline{ w^2 }^k \right) 4f2617d475 Jeff*1687 :label: KE_discrete 1688 9ce7d74115 Jeff*1689 .. _mom_diagnostics: 1690 4f2617d475 Jeff*1691 Mom Diagnostics 1692 --------------- 1693 1694 :: 1695 1696 1697 ------------------------------------------------------------------------ ab47de63dc Mart*1698 <-Name->|Levs|<-parsing code->|<-- Units -->|<- Tile (max=80c) 4f2617d475 Jeff*1699 ------------------------------------------------------------------------ 1700 VISCAHZ | 15 |SZ MR |m^2/s |Harmonic Visc Coefficient (m2/s) (Zeta Pt) 1701 VISCA4Z | 15 |SZ MR |m^4/s |Biharmonic Visc Coefficient (m4/s) (Zeta Pt) 1702 VISCAHD | 15 |SM MR |m^2/s |Harmonic Viscosity Coefficient (m2/s) (Div Pt) 1703 VISCA4D | 15 |SM MR |m^4/s |Biharmonic Viscosity Coefficient (m4/s) (Div Pt) 1704 VAHZMAX | 15 |SZ MR |m^2/s |CFL-MAX Harm Visc Coefficient (m2/s) (Zeta Pt) 1705 VA4ZMAX | 15 |SZ MR |m^4/s |CFL-MAX Biharm Visc Coefficient (m4/s) (Zeta Pt) 1706 VAHDMAX | 15 |SM MR |m^2/s |CFL-MAX Harm Visc Coefficient (m2/s) (Div Pt) 1707 VA4DMAX | 15 |SM MR |m^4/s |CFL-MAX Biharm Visc Coefficient (m4/s) (Div Pt) 1708 VAHZMIN | 15 |SZ MR |m^2/s |RE-MIN Harm Visc Coefficient (m2/s) (Zeta Pt) 1709 VA4ZMIN | 15 |SZ MR |m^4/s |RE-MIN Biharm Visc Coefficient (m4/s) (Zeta Pt) 1710 VAHDMIN | 15 |SM MR |m^2/s |RE-MIN Harm Visc Coefficient (m2/s) (Div Pt) 1711 VA4DMIN | 15 |SM MR |m^4/s |RE-MIN Biharm Visc Coefficient (m4/s) (Div Pt) 1712 VAHZLTH | 15 |SZ MR |m^2/s |Leith Harm Visc Coefficient (m2/s) (Zeta Pt) 1713 VA4ZLTH | 15 |SZ MR |m^4/s |Leith Biharm Visc Coefficient (m4/s) (Zeta Pt) 1714 VAHDLTH | 15 |SM MR |m^2/s |Leith Harm Visc Coefficient (m2/s) (Div Pt) 1715 VA4DLTH | 15 |SM MR |m^4/s |Leith Biharm Visc Coefficient (m4/s) (Div Pt) 1716 VAHZLTHD| 15 |SZ MR |m^2/s |LeithD Harm Visc Coefficient (m2/s) (Zeta Pt) 1717 VA4ZLTHD| 15 |SZ MR |m^4/s |LeithD Biharm Visc Coefficient (m4/s) (Zeta Pt) 1718 VAHDLTHD| 15 |SM MR |m^2/s |LeithD Harm Visc Coefficient (m2/s) (Div Pt) 1719 VA4DLTHD| 15 |SM MR |m^4/s |LeithD Biharm Visc Coefficient (m4/s) (Div Pt) 1720 VAHZSMAG| 15 |SZ MR |m^2/s |Smagorinsky Harm Visc Coefficient (m2/s) (Zeta Pt) 1721 VA4ZSMAG| 15 |SZ MR |m^4/s |Smagorinsky Biharm Visc Coeff. (m4/s) (Zeta Pt) 1722 VAHDSMAG| 15 |SM MR |m^2/s |Smagorinsky Harm Visc Coefficient (m2/s) (Div Pt) 1723 VA4DSMAG| 15 |SM MR |m^4/s |Smagorinsky Biharm Visc Coeff. (m4/s) (Div Pt) 1724 momKE | 15 |SM MR |m^2/s^2 |Kinetic Energy (in momentum Eq.) 1725 momHDiv | 15 |SM MR |s^-1 |Horizontal Divergence (in momentum Eq.) 1726 momVort3| 15 |SZ MR |s^-1 |3rd component (vertical) of Vorticity 1727 Strain | 15 |SZ MR |s^-1 |Horizontal Strain of Horizontal Velocities 1728 Tension | 15 |SM MR |s^-1 |Horizontal Tension of Horizontal Velocities 1729 UBotDrag| 15 |UU 129MR |m/s^2 |U momentum tendency from Bottom Drag 1730 VBotDrag| 15 |VV 128MR |m/s^2 |V momentum tendency from Bottom Drag 1731 USidDrag| 15 |UU 131MR |m/s^2 |U momentum tendency from Side Drag 1732 VSidDrag| 15 |VV 130MR |m/s^2 |V momentum tendency from Side Drag 1733 Um_Diss | 15 |UU 133MR |m/s^2 |U momentum tendency from Dissipation 1734 Vm_Diss | 15 |VV 132MR |m/s^2 |V momentum tendency from Dissipation 1735 Um_Advec| 15 |UU 135MR |m/s^2 |U momentum tendency from Advection terms 1736 Vm_Advec| 15 |VV 134MR |m/s^2 |V momentum tendency from Advection terms 1737 Um_Cori | 15 |UU 137MR |m/s^2 |U momentum tendency from Coriolis term 1738 Vm_Cori | 15 |VV 136MR |m/s^2 |V momentum tendency from Coriolis term 1739 Um_Ext | 15 |UU 137MR |m/s^2 |U momentum tendency from external forcing 1740 Vm_Ext | 15 |VV 138MR |m/s^2 |V momentum tendency from external forcing 1741 Um_AdvZ3| 15 |UU 141MR |m/s^2 |U momentum tendency from Vorticity Advection 1742 Vm_AdvZ3| 15 |VV 140MR |m/s^2 |V momentum tendency from Vorticity Advection 1743 Um_AdvRe| 15 |UU 143MR |m/s^2 |U momentum tendency from vertical Advection (Explicit part) 1744 Vm_AdvRe| 15 |VV 142MR |m/s^2 |V momentum tendency from vertical Advection (Explicit part) 1745 ADVx_Um | 15 |UM 145MR |m^4/s^2 |Zonal Advective Flux of U momentum 1746 ADVy_Um | 15 |VZ 144MR |m^4/s^2 |Meridional Advective Flux of U momentum 1747 ADVrE_Um| 15 |WU LR |m^4/s^2 |Vertical Advective Flux of U momentum (Explicit part) 1748 ADVx_Vm | 15 |UZ 148MR |m^4/s^2 |Zonal Advective Flux of V momentum 1749 ADVy_Vm | 15 |VM 147MR |m^4/s^2 |Meridional Advective Flux of V momentum 1750 ADVrE_Vm| 15 |WV LR |m^4/s^2 |Vertical Advective Flux of V momentum (Explicit part) 1751 VISCx_Um| 15 |UM 151MR |m^4/s^2 |Zonal Viscous Flux of U momentum 1752 VISCy_Um| 15 |VZ 150MR |m^4/s^2 |Meridional Viscous Flux of U momentum 1753 VISrE_Um| 15 |WU LR |m^4/s^2 |Vertical Viscous Flux of U momentum (Explicit part) 1754 VISrI_Um| 15 |WU LR |m^4/s^2 |Vertical Viscous Flux of U momentum (Implicit part) 1755 VISCx_Vm| 15 |UZ 155MR |m^4/s^2 |Zonal Viscous Flux of V momentum 1756 VISCy_Vm| 15 |VM 154MR |m^4/s^2 |Meridional Viscous Flux of V momentum 1757 VISrE_Vm| 15 |WV LR |m^4/s^2 |Vertical Viscous Flux of V momentum (Explicit part) 1758 VISrI_Vm| 15 |WV LR |m^4/s^2 |Vertical Viscous Flux of V momentum (Implicit part) 1759 d67096e55c Jeff*1760 .. _vec_invar_mom_eqs: 4f2617d475 Jeff*1761 1762 Vector invariant momentum equations 1763 =================================== 1764 1765 The finite volume method lends itself to describing the continuity and 1766 tracer equations in curvilinear coordinate systems. However, in 1767 curvilinear coordinates many new metric terms appear in the momentum 1768 equations (written in Lagrangian or flux-form) making generalization far 1769 from elegant. Fortunately, an alternative form of the equations, the 1770 vector invariant equations are exactly that; invariant under coordinate 1771 transformations so that they can be applied uniformly in any orthogonal 1772 curvilinear coordinate system such as spherical coordinates, boundary 1773 following or the conformal spherical cube system. 1774 1775 The non-hydrostatic vector invariant equations read: 1776 1777 .. math:: 0bad585a21 Navi*1778 \partial_t \vec{\bf v} + ( 2\vec{\boldsymbol{\Omega}} + \vec{\boldsymbol{\zeta}}) \times \vec{\bf v} 1779 - b \hat{\bf r} 1780 + \nabla B = \nabla \cdot \vec{\boldsymbol{\tau}} 4f2617d475 Jeff*1781 :label: vect_invar_nh 1782 1783 which describe motions in any orthogonal curvilinear coordinate system. ab47de63dc Mart*1784 Here, :math:`B` is the Bernoulli function and :math:`\vec{\boldsymbol{\zeta}}= \nabla 0bad585a21 Navi*1785 \times \vec{\bf v}` is the vorticity vector. We can take advantage of the 4f2617d475 Jeff*1786 elegance of these equations when discretizing them and use the discrete 1787 definitions of the grad, curl and divergence operators to satisfy 1788 constraints. We can also consider the analogy to forming derived 1789 equations, such as the vorticity equation, and examine how the 1790 discretization can be adjusted to give suitable vorticity advection 1791 among other things. 1792 1793 The underlying algorithm is the same as for the flux form equations. All** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 1795.
1794 that has changed is the contents of the “G’s”. For the time-being, only 1795 the hydrostatic terms have been coded but we will indicate the points 1796 where non-hydrostatic contributions will enter: 1797 1798 .. math:: 1799 G_u = G_u^{fv} + G_u^{\zeta_3 v} + G_u^{\zeta_2 w} + G_u^{\partial_x B} 0bad585a21 Navi*1800 + G_u^{\partial_z \tau^x} + G_u^{h- \rm diss} + G_u^{v- \rm diss} 4f2617d475 Jeff*1801 :label: gu_vecinv 1802 1803 .. math:: 1804 G_v = G_v^{fu} + G_v^{\zeta_3 u} + G_v^{\zeta_1 w} + G_v^{\partial_y B} 0bad585a21 Navi*1805 + G_v^{\partial_z \tau^y} + G_v^{h- \rm diss} + G_v^{v- \rm diss} 4f2617d475 Jeff*1806 :label: gv_vecinv 1807 1808 .. math:: 1809 G_w = G_w^{fu} + G_w^{\zeta_1 v} + G_w^{\zeta_2 u} + G_w^{\partial_z B} 0bad585a21 Navi*1810 + G_w^{h- \rm diss} + G_w^{v- \rm diss} 4f2617d475 Jeff*1811 :label: gw_vecinv 1812 1813 .. admonition:: S/R :filelink:`MOM_VECINV <pkg/mom_vecinv/mom_vecinv.F>` 1814 :class: note 1815 1816 | :math:`G_u` : :varlink:`gU` ( :filelink:`DYNVARS.h <model/inc/DYNVARS.h>` ) 1817 | :math:`G_v` : :varlink:`gV` ( :filelink:`DYNVARS.h <model/inc/DYNVARS.h>` ) 1818 | :math:`G_w` : :varlink:`gW` ( :filelink:`NH_VARS.h <model/inc/NH_VARS.h>` ) 1819 1820 Relative vorticity 1821 ------------------ 1822 1823 The vertical component of relative vorticity is explicitly calculated 1824 and use in the discretization. The particular form is crucial for 1825 numerical stability; alternative definitions break the conservation 1826 properties of the discrete equations. 1827 1828 Relative vorticity is defined: 1829 1830 .. math:: 1831 \zeta_3 = \frac{\Gamma}{A_\zeta} 1832 = \frac{1}{{\cal A}_\zeta} ( \delta_i \Delta y_c v - \delta_j \Delta x_c u ) 1833 :label: zeta3 1834 1835 where :math:`{\cal A}_\zeta` is the area of the vorticity cell 1836 presented in the vertical and :math:`\Gamma` is the circulation about 1837 that cell. 1838 1839 .. admonition:: S/R :filelink:`MOM_CALC_RELVORT3 <pkg/mom_common/mom_calc_relvort3.F>` 1840 :class: note 1841 1842 | :math:`\zeta_3` : :varlink:`vort3` ( local to :filelink:`MOM_VECINV.F <pkg/mom_vecinv/mom_vecinv.F>` ) 1843 1844 Kinetic energy 1845 -------------- 1846 0bad585a21 Navi*1847 The kinetic energy, denoted :math:`\mathrm{KE}`, is defined: 4f2617d475 Jeff*1848 1849 .. math:: ab47de63dc Mart*1850 \mathrm{KE} = \frac{1}{2} ( \overline{ u^2 }^i + \overline{ v^2 }^j 0bad585a21 Navi*1851 + \epsilon_{\rm nh} \overline{ w^2 }^k ) 4f2617d475 Jeff*1852 :label: KE_vecinv 1853 1854 .. admonition:: S/R :filelink:`MOM_CALC_KE <pkg/mom_common/mom_calc_KE.F>` 1855 :class: note 1856 0bad585a21 Navi*1857 | :math:`\mathrm{KE}` : :varlink:`KE` ( local to :filelink:`MOM_VECINV.F <pkg/mom_vecinv/mom_vecinv.F>` ) 4f2617d475 Jeff*1858 1859 Coriolis terms 1860 -------------- 1861 1862 The potential enstrophy conserving form of the linear Coriolis terms are 1863 written: 1864 1865 .. math:: 1866 G_u^{fv} = \frac{1}{\Delta x_c} 1867 \overline{ \frac{f}{h_\zeta} }^j \overline{ \overline{ \Delta x_g h_s v }^j }^i 1868 :label: gu_fv 1869 1870 .. math:: 1871 G_v^{fu} = - \frac{1}{\Delta y_c} 1872 \overline{ \frac{f}{h_\zeta} }^i \overline{ \overline{ \Delta y_g h_w u }^i }^j 1873 :label: gv_fv 1874 1875 Here, the Coriolis parameter :math:`f` is defined at vorticity (corner) 1876 points. 1877 1878 The potential enstrophy conserving form of the non-linear Coriolis terms 1879 are written: 1880 1881 .. math:: 1882 G_u^{\zeta_3 v} = \frac{1}{\Delta x_c} 1883 \overline{ \frac{\zeta_3}{h_\zeta} }^j \overline{ \overline{ \Delta x_g h_s v }^j }^i 1884 :label: gu_zeta3 1885 1886 .. math:: 1887 G_v^{\zeta_3 u} = - \frac{1}{\Delta y_c} 1888 \overline{ \frac{\zeta_3}{h_\zeta} }^i \overline{ \overline{ \Delta y_g h_w u }^i }^j 1889 :label: gv_zeta3 1890 1891 The Coriolis terms can also be evaluated together and expressed in terms 1892 of absolute vorticity :math:`f+\zeta_3`. The potential enstrophy 1893 conserving form using the absolute vorticity is written: 1894 1895 .. math:: 1896 G_u^{fv} + G_u^{\zeta_3 v} = \frac{1}{\Delta x_c} 1897 \overline{ \frac{f + \zeta_3}{h_\zeta} }^j \overline{ \overline{ \Delta x_g h_s v }^j }^i 1898 :label: gu_together 1899 1900 .. math:: 1901 G_v^{fu} + G_v^{\zeta_3 u} = - \frac{1}{\Delta y_c} 1902 \overline{ \frac{f + \zeta_3}{h_\zeta} }^i \overline{ \overline{ \Delta y_g h_w u }^i }^j 1903 :label: gv_together 1904 1905 1906 The distinction between using absolute vorticity or relative vorticity 1907 is useful when constructing higher order advection schemes; monotone 1908 advection of relative vorticity behaves differently to monotone 1909 advection of absolute vorticity. Currently the choice of 1910 relative/absolute vorticity, centered/upwind/high order advection is 1911 available only through commented subroutine calls. 1912 1913 .. admonition:: S/R :filelink:`MOM_VI_CORIOLIS <pkg/mom_vecinv/mom_vi_coriolis.F>`, :filelink:`MOM_VI_U_CORIOLIS <pkg/mom_vecinv/mom_vi_u_coriolis.F>`, :filelink:`MOM_VI_V_CORIOLIS <pkg/mom_vecinv/mom_vi_v_coriolis.F>` 1914 :class: note 1915 1916 | :math:`G_u^{fv} , G_u^{\zeta_3 v}` : :varlink:`uCf` ( local to :filelink:`MOM_VECINV.F <pkg/mom_vecinv/mom_vecinv.F>` ) 1917 | :math:`G_v^{fu} , G_v^{\zeta_3 u}` : :varlink:`vCf` ( local to :filelink:`MOM_VECINV.F <pkg/mom_vecinv/mom_vecinv.F>` ) 1918 1919 Shear terms 1920 ----------- 1921 1922 The shear terms (:math:`\zeta_2w` and :math:`\zeta_1w`) are are 1923 discretized to guarantee that no spurious generation of kinetic energy 1924 is possible; the horizontal gradient of Bernoulli function has to be 1925 consistent with the vertical advection of shear: 1926 1927 .. math:: 1928 G_u^{\zeta_2 w} = \frac{1}{ {\cal A}_w \Delta r_f h_w } \overline{ 0bad585a21 Navi*1929 \overline{ {\cal A}_c w }^i ( \delta_k u - \epsilon_{\rm nh} \delta_i w ) }^k 4f2617d475 Jeff*1930 :label: gu_zeta2w 1931 1932 .. math:: 1933 G_v^{\zeta_1 w} = \frac{1}{ {\cal A}_s \Delta r_f h_s } \overline{ 0bad585a21 Navi*1934 \overline{ {\cal A}_c w }^j ( \delta_k v - \epsilon_{\rm nh} \delta_j w ) }^k 4f2617d475 Jeff*1935 :label: gv_zeta1w 1936 1937 .. admonition:: S/R :filelink:`MOM_VI_U_VERTSHEAR <pkg/mom_vecinv/mom_vi_u_vertshear.F>`, :filelink:`MOM_VI_V_VERTSHEAR <pkg/mom_vecinv/mom_vi_v_vertshear.F>` 1938 :class: note 1939 1940 | :math:`G_u^{\zeta_2 w}` : :varlink:`uCf` ( local to :filelink:`MOM_VECINV.F <pkg/mom_vecinv/mom_vecinv.F>` ) 1941 | :math:`G_v^{\zeta_1 w}` : :varlink:`vCf` ( local to :filelink:`MOM_VECINV.F <pkg/mom_vecinv/mom_vecinv.F>` ) 1942 1943 1944 Gradient of Bernoulli function 1945 ------------------------------ 1946 1947 .. math:: 1948 G_u^{\partial_x B} = 0bad585a21 Navi*1949 \frac{1}{\Delta x_c} \delta_i ( \phi' + \mathrm{KE} ) 4f2617d475 Jeff*1950 :label: gu_dBdx 1951 1952 .. math:: 1953 G_v^{\partial_y B} = 0bad585a21 Navi*1954 \frac{1}{\Delta x_y} \delta_j ( \phi' + \mathrm{KE} ) 4f2617d475 Jeff*1955 :label: gv_dBdy 1956 1957 .. admonition:: S/R :filelink:`MOM_VI_U_GRAD_KE <pkg/mom_vecinv/mom_vi_u_grad_ke.F>`, :filelink:`MOM_VI_V_GRAD_KE <pkg/mom_vecinv/mom_vi_v_grad_ke.F>` 1958 :class: note 1959 0bad585a21 Navi*1960 | :math:`G_u^{\partial_x \mathrm{KE}}` : :varlink:`uCf` ( local to :filelink:`MOM_VECINV.F <pkg/mom_vecinv/mom_vecinv.F>` ) 1961 | :math:`G_v^{\partial_y \mathrm{KE}}` : :varlink:`vCf` ( local to :filelink:`MOM_VECINV.F <pkg/mom_vecinv/mom_vecinv.F>` ) 4f2617d475 Jeff*1962 1963 1964 Horizontal divergence 1965 --------------------- 1966 1967 The horizontal divergence, a complimentary quantity to relative 1968 vorticity, is used in parameterizing the Reynolds stresses and is 1969 discretized: 1970 1971 .. math:: 1972 D = \frac{1}{{\cal A}_c h_c} ( 1973 \delta_i \Delta y_g h_w u 1974 + \delta_j \Delta x_g h_s v ) 1975 :label: horiz_div)vecinv 1976 1977 .. admonition:: S/R :filelink:`MOM_CALC_KE <pkg/mom_common/mom_calc_ke.F>` 1978 :class: note 1979 1980 | :math:`D` : :varlink:`hDiv` ( local to :filelink:`MOM_VECINV.F <pkg/mom_vecinv/mom_vecinv.F>` ) 1981 1982 Horizontal dissipation 1983 ---------------------- 1984 1985 The following discretization of horizontal dissipation conserves 1986 potential vorticity (thickness weighted relative vorticity) and 1987 divergence and dissipates energy, enstrophy and divergence squared: 1988 1989 .. math:: 0bad585a21 Navi*1990 G_u^{h- \rm diss} = 4f2617d475 Jeff*1991 \frac{1}{\Delta x_c} \delta_i ( A_D D - A_{D4} D^*) 1992 - \frac{1}{\Delta y_u h_w} \delta_j h_\zeta ( A_\zeta \zeta - A_{\zeta4} \zeta^* ) 1993 :label: gu_h-dissip 1994 1995 .. math:: 0bad585a21 Navi*1996 G_v^{h- \rm diss} = 4f2617d475 Jeff*1997 \frac{1}{\Delta x_v h_s} \delta_i h_\zeta ( A_\zeta \zeta - A_\zeta \zeta^* ) 1998 + \frac{1}{\Delta y_c} \delta_j ( A_D D - A_{D4} D^* ) 1999 :label: gv_h-dissip 2000 2001 where 2002 2003 .. math:: 2004 2005 \begin{aligned} 0bad585a21 Navi*2006 D^* & = \frac{1}{{\cal A}_c h_c} ( 4f2617d475 Jeff*2007 \delta_i \Delta y_g h_w \nabla^2 u 2008 + \delta_j \Delta x_g h_s \nabla^2 v ) \\ 0bad585a21 Navi*2009 \zeta^* & = \frac{1}{{\cal A}_\zeta} ( 4f2617d475 Jeff*2010 \delta_i \Delta y_c \nabla^2 v 2011 - \delta_j \Delta x_c \nabla^2 u )\end{aligned} 2012 2013 .. admonition:: S/R :filelink:`MOM_VI_HDISSIP <pkg/mom_vecinv/mom_vi_hdissip.F>` 2014 :class: note 2015 2016 | :math:`G_u^{h-dissip}` : :varlink:`uDissip` ( local to :filelink:`MOM_VI_HDISSIP.F <pkg/mom_vecinv/mom_vi_hdissip.F>` ) 2017 | :math:`G_v^{h-dissip}` : :varlink:`vDissip` ( local to :filelink:`MOM_VI_HDISSIP.F <pkg/mom_vecinv/mom_vi_hdissip.F>` ) 2018 2019 2020 Vertical dissipation 2021 -------------------- 2022 2023 Currently, this is exactly the same code as the flux form equations. 2024 2025 .. math:: 0bad585a21 Navi*2026 G_u^{v- \rm diss} = 4f2617d475 Jeff*2027 \frac{1}{\Delta r_f h_w} \delta_k \tau_{13} 2028 :label: gu_v-dissip 2029 2030 .. math:: 0bad585a21 Navi*2031 G_v^{v- \rm diss} = 4f2617d475 Jeff*2032 \frac{1}{\Delta r_f h_s} \delta_k \tau_{23} 2033 :label: gv_v-dissip 2034 2035 represents the general discrete form of the vertical dissipation terms. 2036 2037 In the interior the vertical stresses are discretized: 2038 2039 .. math:: 2040 2041 \begin{aligned} 0bad585a21 Navi*2042 \tau_{13} & = A_v \frac{1}{\Delta r_c} \delta_k u \\ 2043 \tau_{23} & = A_v \frac{1}{\Delta r_c} \delta_k v\end{aligned} 4f2617d475 Jeff*2044 2045 .. admonition:: S/R :filelink:`MOM_U_RVISCFLUX <pkg/mom_common/mom_u_rviscflux.F>`, :filelink:`MOM_V_RVISCFLUX <pkg/mom_common/mom_u_rviscflux.F>` 2046 :class: note 2047 2048 | :math:`\tau_{13}, \tau_{23}` : :varlink:`vrf` ( local to :filelink:`MOM_VECINV.F <pkg/mom_vecinv/mom_vecinv.F>` ) 2049 2050 .. _tracer_eqns: 2051 2052 Tracer equations 2053 ================ 2054 2055 The basic discretization used for the tracer equations is the second 2056 order piece-wise constant finite volume form of the forced 2057 advection-diffusion equations. There are many alternatives to second 2058 order method for advection and alternative parameterizations for the 2059 sub-grid scale processes. The Gent-McWilliams eddy parameterization, KPP 2060 mixing scheme and PV flux parameterization are all dealt with in 2061 separate sections. The basic discretization of the advection-diffusion 2062 part of the tracer equations and the various advection schemes will be 2063 described here. 2064 68aadaa6bd Phob*2065 .. _sub_tracer_eqns_ab: 2066 4f2617d475 Jeff*2067 Time-stepping of tracers: ABII 2068 ------------------------------ 2069 2070 The default advection scheme is the centered second order method which 2071 requires a second order or quasi-second order time-stepping scheme to be 2072 stable. Historically this has been the quasi-second order 2073 Adams-Bashforth method (ABII) and applied to all terms. For an arbitrary 2074 tracer, :math:`\tau`, the forced advection-diffusion equation reads: 2075 0bad585a21 Navi*2076 .. math:: \partial_t \tau + G_{\rm adv}^\tau = G_{\rm diff}^\tau + G_{\rm forc}^\tau 4f2617d475 Jeff*2077 :label: trac_forced_advdiff 2078 0bad585a21 Navi*2079 where :math:`G_{\rm adv}^\tau`, :math:`G_{\rm diff}^\tau` and 2080 :math:`G_{\rm forc}^\tau` are the tendencies due to advection, diffusion and 4f2617d475 Jeff*2081 forcing, respectively, namely: 2082 2083 .. math:: 0bad585a21 Navi*2084 G_{\rm adv}^\tau = \partial_x (u \tau) + \partial_y (v \tau) + \partial_r (w \tau) 2085 - \tau \nabla \cdot {\bf v} 4f2617d475 Jeff*2086 :label: g_adv-tau 2087 2088 .. math:: 0bad585a21 Navi*2089 G_{\rm diff}^\tau = \nabla \cdot \left ( {\bf K} \nabla \tau \right ) 4f2617d475 Jeff*2090 :label: g_diff-tau 2091 2092 and the forcing can be some arbitrary function of state, time and 2093 space. 2094 0bad585a21 Navi*2095 The term, :math:`\tau \nabla \cdot {\bf v}`, is required to retain local 4f2617d475 Jeff*2096 conservation in conjunction with the linear implicit free-surface. It 2097 only affects the surface layer since the flow is non-divergent 2098 everywhere else. This term is therefore referred to as the surface 2099 correction term. Global conservation is not possible using the flux-form 2100 (as here) and a linearized free-surface 2101 (Griffies and Hallberg (2000) :cite:`griffies:00` , Campin et al. (2004) :cite:`cam:04`). 2102 2103 The continuity equation can be recovered by setting 0bad585a21 Navi*2104 :math:`G_{\rm diff}=G_{\rm forc}=0` and :math:`\tau=1`. 4f2617d475 Jeff*2105 2106 The driver routine that calls the routines to calculate tendencies are 2107 :filelink:`CALC_GT <model/src/calc_gt.F>` and :filelink:`CALC_GS <model/src/calc_gs.F>` for temperature and salt (moisture), 2108 respectively. These in turn call a generic advection diffusion routine 2109 :filelink:`GAD_CALC_RHS <pkg/generic_advdiff/gad_calc_rhs.F>` that is called with the flow field and relevant 2110 tracer as arguments and returns the collective tendency due to advection 2111 and diffusion. Forcing is add subsequently in :filelink:`CALC_GT <model/src/calc_gt.F>` 2112 or :filelink:`CALC_GS <model/src/calc_gs.F>` to the same tendency array. 2113 2114 .. admonition:: S/R :filelink:`GAD_CALC_RHS <pkg/generic_advdiff/gad_calc_rhs.F>` 2115 :class: note 2116 2117 | :math:`\tau` : :varlink:`tau` ( argument ) 2118 | :math:`G^{(n)}` : :varlink:`gTracer` ( argument ) 2119 | :math:`F_r` : :varlink:`fVerT` ( argument ) 2120 2121 The space and time discretization are treated separately (method of 2122 lines). Tendencies are calculated at time levels :math:`n` and 2123 :math:`n-1` and extrapolated to :math:`n+1/2` using the Adams-Bashforth 2124 method: 2125 2126 .. math:: ab47de63dc Mart*2127 G^{(n+1/2)} = 0bad585a21 Navi*2128 (\tfrac{3}{2} + \epsilon) G^{(n)} - (\tfrac{1}{2} + \epsilon) G^{(n-1)} 4f2617d475 Jeff*2129 :label: g_nponehalf 2130 0bad585a21 Navi*2131 where :math:`G^{(n)} = G_{\rm adv}^\tau + G_{\rm diff}^\tau + G_{\rm src}^\tau` at 4f2617d475 Jeff*2132 time step :math:`n`. The tendency at :math:`n-1` is not re-calculated 2133 but rather the tendency at :math:`n` is stored in a global array for 2134 later re-use. 2135 2136 .. admonition:: S/R :filelink:`ADAMS_BASHFORTH2 <model/src/gad_calc_rhs.F>` 2137 :class: note 2138 2139 | :math:`G^{(n+1/2)}` : :varlink:`gTracer` ( argument on exit ) 2140 | :math:`G^{(n)}` : :varlink:`gTracer` ( argument on entry ) 2141 | :math:`G^{(n-1)}` : :varlink:`gTrNm1` ( argument ) 2142 | :math:`\epsilon` : :varlink:`ABeps` ( :filelink:`PARAMS.h <model/inc/PARAMS.h>` ) 2143 2144 The tracers are stepped forward in time using the extrapolated tendency: 2145 2146 .. math:: \tau^{(n+1)} = \tau^{(n)} + \Delta t G^{(n+1/2)} 2147 :label: tau_np1 2148 2149 .. admonition:: S/R :filelink:`TIMESTEP_TRACER <model/src/timestep_tracer.F>` 2150 :class: note 2151 2152 | :math:`\tau^{(n+1)}` : :varlink:`gTracer` ( argument on exit ) 2153 | :math:`\tau^{(n)}` : :varlink:`tracer` ( argument on entry ) 2154 | :math:`G^{(n+1/2)}` : :varlink:`gTracer` ( argument ) 2155 | :math:`\Delta t` : :varlink:`deltaTtracer` ( :filelink:`PARAMS.h <model/inc/PARAMS.h>` ) 2156 2157 Strictly speaking the ABII scheme should be applied only to the 2158 advection terms. However, this scheme is only used in conjunction with 2159 the standard second, third and fourth order advection schemes. Selection 2160 of any other advection scheme disables Adams-Bashforth for tracers so 2161 that explicit diffusion and forcing use the forward method. 2162 68aadaa6bd Phob*2163 .. _advection_schemes: 4f2617d475 Jeff*2164 68aadaa6bd Phob*2165 Advection schemes 2166 ================= 4f2617d475 Jeff*2167 ab47de63dc Mart*2168 .. toctree:: 68aadaa6bd Phob*2169 adv-schemes.rst 4f2617d475 Jeff*2170 025afa4065 Jeff*2171 .. _shapiro_filter: 4f2617d475 Jeff*2172 2173 Shapiro Filter 2174 ============== 2175 2176 The Shapiro filter (Shapiro 1970) :cite:`shapiro:70` is a high order 2177 horizontal filter that efficiently remove small scale grid noise without 2178 affecting the physical structures of a field. It is applied at the end 2179 of the time step on both velocity and tracer fields. 2180 2181 Three different space operators are considered here (S1,S2 and S4). They 2182 differ essentially by the sequence of derivative in both X and Y 2183 directions. Consequently they show different damping response function 2184 specially in the diagonal directions X+Y and X-Y. 2185 2186 Space derivatives can be computed in the real space, taking into account 2187 the grid spacing. Alternatively, a pure computational filter can be 2188 defined, using pure numerical differences and ignoring grid spacing. 2189 This later form is stable whatever the grid is, and therefore specially 2190 useful for highly anisotropic grid such as spherical coordinate grid. A 2191 damping time-scale parameter :math:`\tau_{shap}` defines the strength of 2192 the filter damping. 2193 2194 The three computational filter operators are : 2195 2196 .. math:: 0bad585a21 Navi*2197 \begin{aligned} 2198 & \mathrm{S1c:}\hspace{2cm} 2199 [1 - 1/2 \frac{\Delta t}{\tau_{\rm Shap}} ab47de63dc Mart*2200 \{ (\frac{1}{4}\delta_{ii})^n 0bad585a21 Navi*2201 + (\frac{1}{4}\delta_{jj})^n \} ]\\ 2202 & \mathrm{S2c:}\hspace{2cm} ab47de63dc Mart*2203 [1 - \frac{\Delta t}{\tau_{\rm Shap}} 0bad585a21 Navi*2204 \{ \frac{1}{8} (\delta_{ii} + \delta_{jj}) \}^n]\\ 2205 & \mathrm{S4c:}\hspace{2cm} 2206 [1 - \frac{\Delta t}{\tau_{\rm Shap}} (\frac{1}{4}\delta_{ii})^n] ab47de63dc Mart*2207 [1 - \frac{\Delta t}{\tau_{\rm Shap}} (\frac{1}{4}\delta_{jj})^n]\end{aligned} 4f2617d475 Jeff*2208 2209 In addition, the S2 operator can easily be extended to a physical space 2210 filter: 2211 2212 .. math:: 2213 2214 \mathrm{S2g:}\hspace{2cm} 0bad585a21 Navi*2215 [1 - \frac{\Delta t}{\tau_{\rm Shap}} 2216 \{ \frac{L_{\rm Shap}^2}{8} \overline{\nabla}^2 \}^n] 4f2617d475 Jeff*2217 2218 with the Laplacian operator :math:`\overline{\nabla}^2` and a length 0bad585a21 Navi*2219 scale parameter :math:`L_{\rm Shap}`. The stability of this S2g filter 2220 requires :math:`L_{\rm Shap} < \mathrm{Min}^{(\rm Global)}(\Delta x,\Delta y)`. 4f2617d475 Jeff*2221 9ce7d74115 Jeff*2222 .. _shapiro_diagnostics: 2223 4f2617d475 Jeff*2224 SHAP Diagnostics 2225 ---------------- 2226 2227 :: 2228 2229 -------------------------------------------------------------- ab47de63dc Mart*2230 <-Name->|Levs|parsing code|<-Units->|<- Tile (max=80c) 4f2617d475 Jeff*2231 -------------------------------------------------------------- 2232 SHAP_dT | 5 |SM MR |K/s |Temperature Tendency due to Shapiro Filter 2233 SHAP_dS | 5 |SM MR |g/kg/s |Specific Humidity Tendency due to Shapiro Filter 2234 SHAP_dU | 5 |UU 148MR |m/s^2 |Zonal Wind Tendency due to Shapiro Filter 2235 SHAP_dV | 5 |VV 147MR |m/s^2 |Meridional Wind Tendency due to Shapiro Filter 2236 9c8516d9da Jeff*2237 .. _nonlinear_vis_schemes: 4f2617d475 Jeff*2238 2239 Nonlinear Viscosities for Large Eddy Simulation 2240 =============================================== 2241 2242 In Large Eddy Simulations (LES), a turbulent closure needs to be 2243 provided that accounts for the effects of subgridscale motions on the 2244 large scale. With sufficiently powerful computers, we could resolve the 2245 entire flow down to the molecular viscosity scales 2246 (:math:`L_{\nu}\approx 1 \rm cm`). Current computation allows perhaps 2247 four decades to be resolved, so the largest problem computationally 2248 feasible would be about 10m. Most oceanographic problems are much larger 2249 in scale, so some form of LES is required, where only the largest scales 2250 of motion are resolved, and the subgridscale effects on the 2251 large-scale are parameterized. 2252 2253 To formalize this process, we can introduce a filter over the 2254 subgridscale L: :math:`u_\alpha\rightarrow \overline{u_\alpha}` and L: 2255 :math:`b\rightarrow \overline{b}`. This filter has some intrinsic length and time 2256 scales, and we assume that the flow at that scale can be characterized 2257 with a single velocity scale (:math:`V`) and vertical buoyancy gradient 2258 (:math:`N^2`). The filtered equations of motion in a local Mercator 2259 projection about the gridpoint in question (see Appendix for notation 2260 and details of approximation) are: 2261 2262 .. math:: 2263 {\frac{{ \overline{D} {{\tilde {\overline{u}}}}}} {{\overline{Dt}}}} - \frac{{{\tilde {\overline{v}}}} 2264 \sin\theta}{{\rm Ro}\sin\theta_0} + \frac{{M_{Ro}}}{{\rm Ro}} \frac{\partial{\overline{\pi}}}{\partial{x}} 2265 = -\left({\overline{\frac{D{\tilde u}}{Dt} }} - {\frac{{\overline{D} {{\tilde {\overline{u}}}}}}{{\overline{Dt}}} }\right) 2266 +\frac{\nabla^2{{\tilde {\overline{u}}}}}{{\rm Re}} 2267 :label: mercat 2268 2269 .. math:: 2270 {\frac{{ \overline{D} {{\tilde {\overline{v}}}}}} {{\overline{Dt}}}} - \frac{{{\tilde {\overline{u}}}} 2271 \sin\theta}{{\rm Ro}\sin\theta_0} + \frac{{M_{Ro}}}{{\rm Ro}} \frac{\partial{\overline{\pi}}}{\partial{y}} 2272 = -\left({\overline{\frac{D{\tilde v}}{Dt} }} - {\frac{{\overline{D} {{\tilde {\overline{v}}}}}}{{\overline{Dt}}} }\right) 2273 +\frac{\nabla^2{{\tilde {\overline{v}}}}}{{\rm Re}} 2274 :label: mercat_v 2275 2276 .. math:: 2277 \frac{{\overline{D} \overline w}}{{\overline{Dt}}} + \frac{ \frac{\partial{\overline{\pi}}}{\partial{z}} - \overline b}{{\rm Fr}^2\lambda^2} 2278 = -\left(\overline{\frac{D{w}}{Dt}} - \frac{{\overline{D} \overline w}}{{\overline{Dt}}}\right) 0bad585a21 Navi*2279 +\frac{\nabla^2 \overline w}{{\rm Re}} 4f2617d475 Jeff*2280 :label: mercat_w 2281 2282 .. math:: 2283 \frac{{\overline{D} \bar b}}{{\overline{Dt}}} + \overline w = 2284 -\left(\overline{\frac{D{b}}{Dt}} - \frac{{\overline{D} \bar b}}{{\overline{Dt}}}\right) 0bad585a21 Navi*2285 +\frac{\nabla^2 \overline b}{\Pr{\rm Re}} 4f2617d475 Jeff*2286 :label: mercat_b 2287 2288 .. math:: 2289 \mu^2\left({\frac{\partial{\tilde {\overline{u}}}}{\partial{x}}} + {\frac{\partial{\tilde {\overline{v}}}}{\partial{y}}} \right) 2290 + {\frac{\partial{\overline w}}{\partial{z}}} = 0 2291 :label: cont_bfk 2292 2293 Tildes denote multiplication by :math:`\cos\theta/\cos\theta_0` to 2294 account for converging meridians. 2295 2296 The ocean is usually turbulent, and an operational definition of** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 2298.
2297 turbulence is that the terms in parentheses (the ’eddy’ terms) on the 2298 right of :eq:`mercat` - :eq:`mercat_b`) are of comparable magnitude to the terms on the 2299 left-hand side. The terms proportional to the inverse of , instead, are 2300 many orders of magnitude smaller than all of the other terms in 2301 virtually every oceanic application. 2302 2303 Eddy Viscosity 2304 -------------- 2305** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 2307.
2306 A turbulent closure provides an approximation to the ’eddy’ terms on the 2307 right of the preceding equations. The simplest form of LES is just to 2308 increase the viscosity and diffusivity until the viscous and diffusive 2309 scales are resolved. That is, we approximate :eq:`mercat` - :eq:`mercat_b`: 2310 2311 .. math:: 2312 \left({\overline{\frac{D{\tilde u}}{Dt} }} - {\frac{{\overline{D} {{\tilde {\overline{u}}}}}}{{\overline{Dt}}} }\right) 2313 \approx\frac{\nabla^2_h{{\tilde {\overline{u}}}}}{{\rm Re}_h} 2314 +\frac{{\frac{\partial^2{{\tilde {\overline{u}}}}}{{\partial{z}}^2}}}{{\rm Re}_v} 2315 :label: eddyvisc_u 2316 2317 .. math:: 2318 \left({\overline{\frac{D{\tilde v}}{Dt} }} - {\frac{{\overline{D} {{\tilde {\overline{v}}}}}}{{\overline{Dt}}} }\right) 2319 \approx\frac{\nabla^2_h{{\tilde {\overline{v}}}}}{{\rm Re}_h} 2320 +\frac{{\frac{\partial^2{{\tilde {\overline{v}}}}}{{\partial{z}}^2}}}{{\rm Re}_v} 2321 :label: eddyvisc_v 2322 2323 .. math:: 2324 \left(\overline{\frac{D{w}}{Dt}} - \frac{{\overline{D} \overline w}}{{\overline{Dt}}}\right) 2325 \approx\frac{\nabla^2_h \overline w}{{\rm Re}_h} 2326 +\frac{{\frac{\partial^2{\overline w}}{{\partial{z}}^2}}}{{\rm Re}_v} 2327 :label: eddyvisc_w 2328 2329 .. math:: 2330 \left(\overline{\frac{D{b}}{Dt}} - \frac{{\overline{D} \bar b}}{{\overline{Dt}}}\right) 2331 \approx\frac{\nabla^2_h \overline b}{\Pr{\rm Re}_h} 0bad585a21 Navi*2332 +\frac{{\frac{\partial^2{\overline b}}{{\partial{z}}^2}}}{\Pr{\rm Re}_v} 4f2617d475 Jeff*2333 :label: eddyvisc_b 2334 2335 Reynolds-Number Limited Eddy Viscosity 2336 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 2337 2338 One way of ensuring that the gridscale is sufficiently viscous (i.e., 2339 resolved) is to choose the eddy viscosity :math:`A_h` so that the 2340 gridscale horizontal Reynolds number based on this eddy viscosity, 2341 :math:`{\rm Re}_h`, is O(1). That is, if the gridscale is to be 2342 viscous, then the viscosity should be chosen to make the viscous terms 2343 as large as the advective ones. Bryan et al. (1975) 2344 :cite:`bryan:75` notes that a computational mode is 2345 squelched by using :math:`{\rm Re}_h<`\ 2. 2346 2347 MITgcm users can select horizontal eddy viscosities based on 2348 :math:`{\rm Re}_h` using two methods. 1) The user may estimate the 2349 velocity scale expected from the calculation and grid spacing and set 2350 :varlink:`viscAh` to satisfy :math:`{\rm Re}_h<2`. 2) The user may use 2351 :varlink:`viscAhReMax`, which ensures that the viscosity is always chosen so that 2352 :math:`{\rm Re}_h<` :varlink:`viscAhReMax`. This last option should be used with 2353 caution, however, since it effectively implies that viscous terms are 2354 fixed in magnitude relative to advective terms. While it may be a useful 2355 method for specifying a minimum viscosity with little effort, tests 2356 Bryan et al. (1975) :cite:`bryan:75` have shown that setting :varlink:`viscAhReMax` =2 2357 often tends to increase the viscosity substantially over other more** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 2359.
2358 ’physical’ parameterizations below, especially in regions where 2359 gradients of velocity are small (and thus turbulence may be weak), so 2360 perhaps a more liberal value should be used, e.g. :varlink:`viscAhReMax` =10. 2361 2362 While it is certainly necessary that viscosity be active at the 2363 gridscale, the wavelength where dissipation of energy or enstrophy 2364 occurs is not necessarily :math:`L=A_h/U`. In fact, it is by ensuring 2365 that either the dissipation of energy in a 3-d turbulent cascade 2366 (Smagorinsky) or dissipation of enstrophy in a 2-d turbulent cascade 2367 (Leith) is resolved that these parameterizations derive their physical 2368 meaning. 2369 2370 Vertical Eddy Viscosities 2371 ~~~~~~~~~~~~~~~~~~~~~~~~~ 2372 2373 Vertical eddy viscosities are often chosen in a more subjective way, as 2374 model stability is not usually as sensitive to vertical viscosity.** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 2376.
2375 Usually the ’observed’ value from finescale measurements is used 2376 (e.g. :varlink:`viscAr`\ :math:`\approx1\times10^{-4} m^2/s`). However, 2377 Smagorinsky (1993) :cite:`smag:93` notes that the Smagorinsky 2378 parameterization of isotropic turbulence implies a value of the vertical 2379 viscosity as well as the horizontal viscosity (see below). 2380 2381 Smagorinsky Viscosity 2382 ~~~~~~~~~~~~~~~~~~~~~ 2383 2384 Some suggest (see Smagorinsky 1963 :cite:`smag:63`; Smagorinsky 1993 :cite:`smag:93`) choosing a viscosity 2385 that *depends on the resolved motions*. Thus, the overall viscous 2386 operator has a nonlinear dependence on velocity. Smagorinsky chose his** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 2388.
2387 form of viscosity by considering Kolmogorov’s ideas about the energy 2388 spectrum of 3-d isotropic turbulence. 2389 2390 Kolmogorov supposed that energy is injected into the flow at** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 2392.
2391 large scales (small :math:`k`) and is ’cascaded’ or transferred 2392 conservatively by nonlinear processes to smaller and smaller scales 2393 until it is dissipated near the viscous scale. By setting the energy 2394 flux through a particular wavenumber :math:`k`, :math:`\epsilon`, to be 2395 a constant in :math:`k`, there is only one combination of viscosity and 2396 energy flux that has the units of length, the Kolmogorov wavelength. It 2397 is :math:`L_\epsilon(\nu)\propto\pi\epsilon^{-1/4}\nu^{3/4}` (the 2398 :math:`\pi` stems from conversion from wavenumber to wavelength). To 2399 ensure that this viscous scale is resolved in a numerical model, the 2400 gridscale should be decreased until :math:`L_\epsilon(\nu)>L` (so-called 2401 Direct Numerical Simulation, or DNS). Alternatively, an eddy viscosity 2402 can be used and the corresponding Kolmogorov length can be made larger 2403 than the gridscale, 2404 :math:`L_\epsilon(A_h)\propto\pi\epsilon^{-1/4}A_h^{3/4}` (for Large 2405 Eddy Simulation or LES). 2406 2407 There are two methods of ensuring that the Kolmogorov length is resolved 2408 in MITgcm. 1) The user can estimate the flux of energy through spectral 2409 space for a given simulation and adjust grid spacing or :varlink:`viscAh` to ensure 2410 that :math:`L_\epsilon(A_h)>L`; 2) The user may use the approach of 2411 Smagorinsky with :varlink:`viscC2Smag`, which estimates the energy flux at every 2412 grid point, and adjusts the viscosity accordingly. 2413 2414 Smagorinsky formed the energy equation from the momentum equations by 2415 dotting them with velocity. There are some complications when using the 2416 hydrostatic approximation as described by Smagorinsky (1993) 2417 :cite:`smag:93`. The positive definite energy 2418 dissipation by horizontal viscosity in a hydrostatic flow is 2419 :math:`\nu D^2`, where D is the deformation rate at the viscous scale.** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 2421.
2420 According to Kolmogorov’s theory, this should be a good approximation to 2421 the energy flux at any wavenumber :math:`\epsilon\approx\nu D^2`. 2422 Kolmogorov and Smagorinsky noted that using an eddy viscosity that 2423 exceeds the molecular value :math:`\nu` should ensure that the energy 2424 flux through viscous scale set by the eddy viscosity is the same as it 2425 would have been had we resolved all the way to the true viscous scale. 2426 That is, :math:`\epsilon\approx 0bad585a21 Navi*2427 A_{h \rm Smag} \overline D^2`. If we use this approximation to estimate the 4f2617d475 Jeff*2428 Kolmogorov viscous length, then 2429 2430 .. math:: 0bad585a21 Navi*2431 L_\epsilon(A_{h \rm Smag})\propto\pi\epsilon^{-1/4}A_{h \rm Smag}^{3/4}\approx\pi(A_{h \rm Smag} 2432 \overline D^2)^{-1/4}A_{h\rm Smag}^{3/4} = \pi A_{h\rm Smag}^{1/2}\overline D^{-1/2} 4f2617d475 Jeff*2433 :label: kolm_visc_len 2434 0bad585a21 Navi*2435 To make :math:`L_\epsilon(A_{h\rm Smag})` scale with the gridscale, then 4f2617d475 Jeff*2436 0bad585a21 Navi*2437 .. math:: A_{h\rm Smag} = \left(\frac{{\sf viscC2Smag}}{\pi}\right)^2L^2|\overline D| 4f2617d475 Jeff*2438 :label: AhSmag 2439 2440 Where the deformation rate appropriate for hydrostatic flows with 2441 shallow-water scaling is 2442 2443 .. math:: 2444 |\overline D|=\sqrt{\left({\frac{\partial{\overline {\tilde u}}}{\partial{x}}} - {\frac{\partial{\overline {\tilde v}}}{\partial{y}}}\right)^2 2445 + \left({\frac{\partial{\overline {\tilde u}}}{\partial{y}}} + {\frac{\partial{\overline {\tilde v}}}{\partial{x}}}\right)^2} 2446 :label: Dbar_mag 2447 2448 The coefficient :varlink:`viscC2Smag` is what an MITgcm user sets, and it replaces 2449 the proportionality in the Kolmogorov length with an equality. Others 2450 (Griffies and Hallberg, 2000 :cite:`griffies:00`) suggest values of :varlink:`viscC2Smag` from 2.2 to 2451 4 for oceanic problems. Smagorinsky (1993) :cite:`smag:93` 2452 shows that values from 0.2 to 0.9 have been used in atmospheric modeling. 2453 2454 Smagorinsky (1993) :cite:`smag:93` shows that a corresponding 2455 vertical viscosity should be used: 2456 2457 .. math:: 0bad585a21 Navi*2458 A_{v \rm Smag} = \left(\frac{{\sf viscC2Smag}}{\pi}\right)^2H^2 4f2617d475 Jeff*2459 \sqrt{\left({\frac{\partial{\overline {\tilde u}}}{\partial{z}}}\right)^2 2460 + \left({\frac{\partial{\overline {\tilde v}}}{\partial{z}}}\right)^2} 2461 :label: A_vsmag 2462 2463 This vertical viscosity is currently not implemented in MITgcm. 2464 9c8516d9da Jeff*2465 .. _leith_viscosity: 2466 4f2617d475 Jeff*2467 Leith Viscosity 2468 ~~~~~~~~~~~~~~~ 2469 9c8516d9da Jeff*2470 Leith (1968, 1996) :cite:`leith:68` :cite:`leith:96` notes that 2-D turbulence is 2471 quite different from 3-D. In 2-D turbulence, energy cascades 4f2617d475 Jeff*2472 to larger scales, so there is no concern about resolving the scales of 2473 energy dissipation. Instead, another quantity, enstrophy, (which is the 9c8516d9da Jeff*2474 vertical component of vorticity squared) is conserved in 2-D turbulence, 4f2617d475 Jeff*2475 and it cascades to smaller scales where it is dissipated. 2476 2477 Following a similar argument to that above about energy flux, the 2478 enstrophy flux is estimated to be equal to the positive-definite 0bad585a21 Navi*2479 gridscale dissipation rate of enstrophy :math:`\eta\approx A_{h \rm Leith} 4f2617d475 Jeff*2480 |\nabla\overline \omega_3|^2`. By dimensional analysis, the 0bad585a21 Navi*2481 enstrophy-dissipation scale is :math:`L_\eta(A_{h \rm Leith})\propto\pi 2482 A_{h \rm Leith}^{1/2}\eta^{-1/6}`. Thus, the Leith-estimated length scale of 4f2617d475 Jeff*2483 enstrophy-dissipation and the resulting eddy viscosity are 2484 2485 .. math:: 0bad585a21 Navi*2486 L_\eta(A_{h \rm Leith})\propto\pi A_{h \rm Leith}^{1/2}\eta^{-1/6} 2487 = \pi A_{h \rm Leith}^{1/3}| \nabla \overline \omega_3|^{-1/3} 4f2617d475 Jeff*2488 :label: L_eta 2489 2490 .. math:: 0bad585a21 Navi*2491 A_{h \rm Leith} = \left(\frac{{\sf viscC2Leith}}{\pi}\right)^3L^3| \nabla \overline\omega_3| 4f2617d475 Jeff*2492 :label: Ah_Leith 2493 2494 .. math:: 0bad585a21 Navi*2495 | \nabla \omega_3| \equiv \sqrt{\left[{\frac{\partial{\ }}{\partial{x}}} 4f2617d475 Jeff*2496 \left({\frac{\partial{\overline {\tilde v}}}{\partial{x}}} - {\frac{\partial{\overline {\tilde u}}}{\partial{y}}}\right)\right]^2 2497 + \left[{\frac{\partial{\ }}{\partial{y}}}\left({\frac{\partial{\overline {\tilde v}}}{\partial{x}}} 2498 - {\frac{\partial{\overline {\tilde u}}}{\partial{y}}}\right)\right]^2} 2499 :label: Leith3 2500 9c8516d9da Jeff*2501 The runtime flag :varlink:`useFullLeith` controls whether or not to calculate the full gradients for the Leith viscosity (.TRUE.) 2502 or to use an approximation (.FALSE.). The only reason to set :varlink:`useFullLeith` = .FALSE. is if your simulation fails when 2503 computing the gradients. This can occur when using the cubed sphere and other complex grids. f59d76b0dd Ed D*2504 4f2617d475 Jeff*2505 Modified Leith Viscosity 2506 ~~~~~~~~~~~~~~~~~~~~~~~~ 2507 2508 The argument above for the Leith viscosity parameterization uses 2509 concepts from purely 2-dimensional turbulence, where the horizontal flow 2510 field is assumed to be non-divergent. However, oceanic flows are only 2511 quasi-two dimensional. While the barotropic flow, or the flow within 2512 isopycnal layers may behave nearly as two-dimensional turbulence, there 2513 is a possibility that these flows will be divergent. In a 2514 high-resolution numerical model, these flows may be substantially 2515 divergent near the grid scale, and in fact, numerical instabilities 2516 exist which are only horizontally divergent and have little vertical 2517 vorticity. This causes a difficulty with the Leith viscosity, which can 2518 only respond to buildup of vorticity at the grid scale. 2519 2520 MITgcm offers two options for dealing with this problem. 1) The 2521 Smagorinsky viscosity can be used instead of Leith, or in conjunction 2522 with Leith -- a purely divergent flow does cause an increase in Smagorinsky 2523 viscosity; 2) The :varlink:`viscC2LeithD` parameter can be set. This is a damping 2524 specifically targeting purely divergent instabilities near the 2525 gridscale. The combined viscosity has the form: 2526 2527 .. math:: 0bad585a21 Navi*2528 A_{h \rm Leith} = L^3\sqrt{\left(\frac{{\sf viscC2Leith}}{\pi}\right)^6 2529 | \nabla \overline \omega_3|^2 + \left(\frac{{\sf viscC2LeithD}}{\pi}\right)^6 2530 | \nabla ( \nabla \cdot \overline {\tilde u}_h)|^2} 4f2617d475 Jeff*2531 :label: Ah_Leithcomb 2532 2533 .. math:: 0bad585a21 Navi*2534 | \nabla ( \nabla \cdot \overline {\tilde u}_h)| \equiv 4f2617d475 Jeff*2535 \sqrt{\left[{\frac{\partial{\ }}{\partial{x}}}\left({\frac{\partial{\overline {\tilde u}}}{\partial{x}}} 2536 + {\frac{\partial{\overline {\tilde v}}}{\partial{y}}}\right)\right]^2 2537 + \left[{\frac{\partial{\ }}{\partial{y}}}\left({\frac{\partial{\overline {\tilde u}}}{\partial{x}}} 2538 + {\frac{\partial{\overline {\tilde v}}}{\partial{y}}}\right)\right]^2} 2539 :label: Leith_mod 2540 2541 Whether there is any physical rationale for this correction is unclear, 2542 but the numerical consequences are good. The divergence 2543 in flows with the grid scale larger or comparable to the Rossby radius 2544 is typically much smaller than the vorticity, so this adjustment only 2545 rarely adjusts the viscosity if :varlink:`viscC2LeithD` = :varlink:`viscC2Leith`. 2546 However, the rare regions where this 2547 viscosity acts are often the locations for the largest vales of vertical 2548 velocity in the domain. Since the CFL condition on vertical velocity is 2549 often what sets the maximum timestep, this viscosity may substantially 2550 increase the allowable timestep without severely compromising the verity 2551 of the simulation. Tests have shown that in some calculations, a 2552 timestep three times larger was allowed when :varlink:`viscC2LeithD` = :varlink:`viscC2Leith`. 2553 f59d76b0dd Ed D*2554 2555 Quasi-Geostrophic Leith Viscosity 2556 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 2557 2558 A variant of Leith viscosity can be derived for quasi-geostrophic dynamics. 2559 This leads to a slightly different equation for the viscosity that includes 2560 a contribution from quasigeostrophic vortex stretching (Bachman et al. 2017 :cite:`bachman:17`). 2561 The viscosity is given by 2562 2563 .. math:: ab47de63dc Mart*2564 \nu_{*} = \left(\frac{\Lambda \Delta s}{\pi}\right)^{3} \left| \nabla _{h}(f\mathbf{\hat{z}}) + 0bad585a21 Navi*2565 \nabla _{h}( \nabla \times \mathbf{v}_{h*}) + \partial_{z}\frac{f}{N^{2}} \nabla _{h} b \right| f59d76b0dd Ed D*2566 :label: bachman2017_eq39 2567 2568 where :math:`\Lambda` is a tunable parameter of :math:`\mathcal{O}(1)`, 2569 :math:`\Delta s = \sqrt{\Delta x \Delta y}` is the grid scale, :math:`f\mathbf{\hat{z}}` 2570 is the vertical component of the Coriolis parameter, :math:`\mathbf{v}_{h*}` is the horizontal velocity, 2571 :math:`N^{2}` is the Brunt-Väisälä frequency, and :math:`b` is the buoyancy. 2572 2573 However, the viscosity given by :eq:`bachman2017_eq39` does not constrain purely 2574 divergent motions. As such, a small :math:`\mathcal{O}(\epsilon)` correction is added 2575 2576 .. math:: ab47de63dc Mart*2577 \nu_{*} = \left(\frac{\Lambda \Delta s}{\pi}\right)^{3} 2578 \sqrt{\left| \nabla _{h}(f\mathbf{\hat{z}}) + \nabla _{h}( \nabla \times \mathbf{v}_{h*}) + 0bad585a21 Navi*2579 \partial_{z} \frac{f}{N^{2}} \nabla _{h} b\right|^{2} + | \nabla ( \nabla \cdot \mathbf{v}_{h})|^{2}} f59d76b0dd Ed D*2580 :label: bachman2017_eq40 2581 2582 This form is, however, numerically awkward; as the Brunt-Väisälä Frequency becomes very small 2583 in regions of weak or vanishing stratification, the vortex stretching term becomes very large. 2584 The resulting large viscosities can lead to numerical instabilities. Bachman et al. (2017) :cite:`bachman:17` 2585 present two limiting forms for the viscosity based on flow parameters such as :math:`Fr_{*}`, 2586 the Froude number, and :math:`Ro_{*}`, the Rossby number. The second of which, 2587 2588 .. math:: 2589 \begin{aligned} 2590 \nu_{*} = & \left(\frac{\Lambda \Delta s}{\pi}\right)^{3} \\ 0bad585a21 Navi*2591 & \sqrt{\min \left( \left| \nabla _{h}q_{2*} + \partial_{z} \frac{f^{2}}{N^{2}} \nabla _{h} b \right|, 2592 \left( 1 + \frac{Fr_{*}^{2}}{Ro_{*}^{2}} + Fr_{*}^{4}\right) | \nabla _{h}q_{2*}|\right)^{2} + | \nabla ( \nabla \cdot \mathbf{v}_{h}) |^{2}} f59d76b0dd Ed D*2593 \end{aligned} 2594 :label: bachman2017_eq56 2595 2596 2597 has been implemented and is active when ``#define`` :varlink:`ALLOW_LEITH_QG` is included 2598 in a copy of :filelink:`MOM_COMMON_OPTIONS.h<pkg/mom_common/MOM_COMMON_OPTIONS.h>` in 2599 a code mods directory (specified through :ref:`-mods <mods_option>` command 2600 line option in :filelink:`genmake2 <tools/genmake2>`). 2601 2602 LeithQG viscosity is designed to work best in simulations that resolve some mesoscale features. 2603 In simulations that are too coarse to permit eddies or fine enough to resolve submesoscale features, 2604 it should fail gracefully. The non-dimensional parameter :varlink:`viscC2LeithQG` corresponds to 2605 :math:`\Lambda` in the above equations and scales the viscosity; the recommended value is 1. 2606 2607 There is no reason to use the quasi-geostrophic form of Leith at the same time as either 2608 standard Leith or modified Leith. Therefore, the model will not run if non-zero values have 2609 been set for these coefficients; the model will stop during the configuration check. 2610 LeithQG can be used regardless of the setting for :varlink:`useFullLeith`. Just as for the 2611 other forms of Leith viscosity, this flag determines whether or not the full gradients are used. 2612 The simplified gradients were originally intended for use on complex grids, but have been 2613 shown to produce better kinetic energy spectra even on very straightforward grids. 2614 2615 To add the LeithQG viscosity to the GMRedi coefficient, as was done in some of the simulations 2616 in Bachman et al. (2017) :cite:`bachman:17`, ``#define`` :varlink:`ALLOW_LEITH_QG` must be specified, 2617 as described above. In addition to this, the compile-time flag :varlink:`ALLOW_GM_LEITH_QG` 2618 must also be defined in a (``-mods``) copy of :filelink:`GMREDI_OPTIONS.h<pkg/gmredi/GMREDI_OPTIONS.h>` 2619 when the model is compiled, and the runtime parameter :varlink:`GM_useLeithQG` set to .TRUE. in ``data.gmredi``. 2620 This will use the value of :varlink:`viscC2LeithQG` specified in the ``data`` input file to compute the coefficient. 2621 9c8516d9da Jeff*2622 .. _CFL_constraint_visc: f59d76b0dd Ed D*2623 4f2617d475 Jeff*** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 2625.
2624 Courant–Freidrichs–Lewy Constraint on Viscosity 2625 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 2626 2627 Whatever viscosities are used in the model, the choice is constrained by** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 2629.
2628 gridscale and timestep by the Courant–Freidrichs–Lewy (CFL) constraint 2629 on stability: 2630 2631 .. math:: 2632 \begin{aligned} 2633 A_h & < \frac{L^2}{4\Delta t} \\ 2634 A_4 & \le \frac{L^4}{32\Delta t}\end{aligned} 2635 2636 The viscosities may be automatically limited to be no greater than 2637 these values in MITgcm by specifying :varlink:`viscAhGridMax` :math:`<1` and 2638 :varlink:`viscA4GridMax` :math:`<1`. Similarly-scaled minimum values of 2639 viscosities are provided by :varlink:`viscAhGridMin` and :varlink:`viscA4GridMin`, which if 2640 used, should be set to values :math:`\ll 1`. :math:`L` is roughly the 2641 gridscale (see below). 2642 2643 Following Griffies and Hallberg (2000) :cite:`griffies:00`, we note that there is a 2644 factor of :math:`\Delta x^2/8` difference between the harmonic and biharmonic viscosities. Thus, 2645 whenever a non-dimensional harmonic coefficient is used in the MITgcm 2646 (e.g. :varlink:`viscAhGridMax` :math:`<1`), the biharmonic equivalent is scaled 2647 so that the same non-dimensional value can be used (e.g. :varlink:`viscA4GridMax` :math:`<1`). 2648 2649 Biharmonic Viscosity 2650 ~~~~~~~~~~~~~~~~~~~~ 2651 2652 Holland (1978) :cite:`holland:78` suggested that eddy viscosities ought to be 2653 focused on the dynamics at the grid scale, as larger motions would be** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 2655.
2654 ’resolved’. To enhance the scale selectivity of the viscous operator, he 2655 suggested a biharmonic eddy viscosity instead of a harmonic (or Laplacian) viscosity: 2656 2657 .. math:: 2658 \left({\overline{\frac{D{\tilde u}}{Dt} }} - {\frac{{\overline{D} {{\tilde {\overline{u}}}}}}{{\overline{Dt}}} }\right) 2659 \approx \frac{-\nabla^4_h{{\tilde {\overline{u}}}}}{{\rm Re}_4} 2660 + \frac{{\frac{\partial^2{{\tilde {\overline{u}}}}}{{\partial{z}}^2}}}{{\rm Re}_v} 2661 :label: bieddyvisc_u 2662 2663 2664 .. math:: 2665 \left({\overline{\frac{D{\tilde v}}{Dt} }} - {\frac{{\overline{D} {{\tilde {\overline{v}}}}}}{{\overline{Dt}}} }\right) 2666 \approx \frac{-\nabla^4_h{{\tilde {\overline{v}}}}}{{\rm Re}_4} 0bad585a21 Navi*2667 + \frac{{\frac{\partial^2{{\tilde {\overline{v}}}}}{{\partial{z}}^2}}}{{\rm Re}_v} 4f2617d475 Jeff*2668 :label: bieddyvisc_v 2669 2670 2671 .. math:: 2672 \left(\overline{\frac{D{w}}{Dt}} - \frac{{\overline{D} \overline w}}{{\overline{Dt}}}\right) 0bad585a21 Navi*2673 \approx\frac{-\nabla^4_h\overline w}{{\rm Re}_4} + \frac{{\frac{\partial^2{\overline w}}{{\partial{z}}^2}}}{{\rm Re}_v} 4f2617d475 Jeff*2674 :label: bieddyvisc_w 2675 2676 .. math:: 2677 \left(\overline{\frac{D{b}}{Dt}} - \frac{{\overline{D} \bar b}}{{\overline{Dt}}}\right) 2678 \approx \frac{-\nabla^4_h \overline b}{\Pr{\rm Re}_4} 0bad585a21 Navi*2679 +\frac{{\frac{\partial^2{\overline b}}{{\partial{z}}^2}}}{\Pr{\rm Re}_v} 4f2617d475 Jeff*2680 :label: bieddyvisc_b 2681 2682 2683 Griffies and Hallberg (2000) :cite:`griffies:00` propose that if one scales the 2684 biharmonic viscosity by stability considerations, then the biharmonic 2685 viscous terms will be similarly active to harmonic viscous terms at the 2686 gridscale of the model, but much less active on larger scale motions. 2687 Similarly, a biharmonic diffusivity can be used for less diffusive 2688 flows. 2689 2690 In practice, biharmonic viscosity and diffusivity allow a less viscous, 2691 yet numerically stable, simulation than harmonic viscosity and 2692 diffusivity. However, there is no physical rationale for such operators 2693 being of leading order, and more boundary conditions must be specified 2694 than for the harmonic operators. If one considers the approximations of 2695 :eq:`eddyvisc_u` - :eq:`eddyvisc_b` and :eq:`bieddyvisc_u` - :eq:`bieddyvisc_b` 2696 to be terms in the Taylor series 2697 expansions of the eddy terms as functions of the large-scale gradient, 2698 then one can argue that both harmonic and biharmonic terms would occur 2699 in the series, and the only question is the choice of coefficients. 2700 Using biharmonic viscosity alone implies that one zeros the first 2701 non-vanishing term in the Taylor series, which is unsupported by any 2702 fluid theory or observation. 2703 2704 Nonetheless, MITgcm supports a plethora of biharmonic viscosities and 2705 diffusivities, which are controlled with parameters named similarly to 2706 the harmonic viscosities and diffusivities with the substitution 2707 h :math:`\rightarrow 4` in the MITgcm parameter name. MITgcm also supports biharmonic Leith and 2708 Smagorinsky viscosities: 2709 2710 .. math:: 0bad585a21 Navi*2711 A_{4 \rm Smag} = \left(\frac{{\sf viscC4Smag}}{\pi}\right)^2\frac{L^4}{8}|D| 4f2617d475 Jeff*2712 :label: A4_Smag 2713 2714 .. math:: 0bad585a21 Navi*2715 A_{4 \rm Leith} = \frac{L^5}{8}\sqrt{\left(\frac{{\sf viscC4Leith}}{\pi}\right)^6 2716 | \nabla \overline \omega_3|^2 + \left(\frac{{\sf viscC4LeithD}}{\pi}\right)^6 2717 | \nabla ( \nabla \cdot \overline {\bf {\tilde u}}_h)|^2} 4f2617d475 Jeff*2718 :label: A4_Leith 2719 2720 However, it should be noted that unlike the harmonic forms, the 2721 biharmonic scaling does not easily relate to whether energy-dissipation 2722 or enstrophy-dissipation scales are resolved. If similar arguments are 2723 used to estimate these scales and scale them to the gridscale, the 2724 resulting biharmonic viscosities should be: 2725 2726 .. math:: 0bad585a21 Navi*2727 A_{4 \rm Smag} = \left(\frac{{\sf viscC4Smag}}{\pi}\right)^5L^5 4f2617d475 Jeff*2728 |\nabla^2\overline {\bf {\tilde u}}_h| 2729 :label: A4_Smag_alt 2730 2731 .. math:: 0bad585a21 Navi*2732 A_{4 \rm Leith} = L^6\sqrt{\left(\frac{{\sf viscC4Leith}}{\pi}\right)^{12} 4f2617d475 Jeff*2733 |\nabla^2 \overline \omega_3|^2 + \left(\frac{{\sf viscC4LeithD}}{\pi}\right)^{12} 0bad585a21 Navi*2734 |\nabla^2 ( \nabla \cdot \overline {\bf {\tilde u}}_h)|^2} 4f2617d475 Jeff*2735 :label: A4_Leith_alt 2736 2737 Thus, the biharmonic scaling suggested by Griffies and Hallberg (2000) 2738 :cite:`griffies:00` implies: 2739 2740 .. math:: 2741 \begin{aligned} 2742 |D| & \propto L|\nabla^2\overline {\bf {\tilde u}}_h|\\ 0bad585a21 Navi*2743 | \nabla \overline \omega_3| & \propto L|\nabla^2 \overline \omega_3|\end{aligned} 4f2617d475 Jeff*2744 2745 It is not at all clear that these assumptions ought to hold. Only the 2746 Griffies and Hallberg (2000) :cite:`griffies:00` forms are currently implemented in 2747 MITgcm. 2748 2749 Selection of Length Scale 2750 ~~~~~~~~~~~~~~~~~~~~~~~~~ 2751 2752 Above, the length scale of the grid has been denoted :math:`L`. However, 2753 in strongly anisotropic grids, :math:`L_x` and :math:`L_y` will be quite 2754 different in some locations. In that case, the CFL condition suggests 2755 that the minimum of :math:`L_x` and :math:`L_y` be used. On the other 2756 hand, other viscosities which involve whether a particular wavelength is** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 2758.
2757 ’resolved’ might be better suited to use the maximum of :math:`L_x` and 2758 :math:`L_y`. Currently, MITgcm uses :varlink:`useAreaViscLength` to select between dcaaa42497 Jeff*2759 two options. If false, the square root of the `harmonic mean <https://en.wikipedia.org/wiki/Harmonic_mean>`_ 2760 of :math:`L^2_x` and 4f2617d475 Jeff*2761 :math:`L^2_y` is used for all viscosities, which is closer to the 2762 minimum and occurs naturally in the CFL constraint. If :varlink:`useAreaViscLength` 2763 is true, then the square root of the area of the grid cell is used. 2764 2765 Mercator, Nondimensional Equations 2766 ---------------------------------- 2767 2768 The rotating, incompressible, Boussinesq equations of motion 2769 (Gill, 1982) :cite:`gill:82` on a sphere can be written in Mercator 2770 projection about a latitude :math:`\theta_0` and geopotential height 2771 :math:`z=r-r_0`. The nondimensional form of these equations is: 2772 2773 .. math:: 2774 {\rm Ro} \frac{D{\tilde u}}{Dt} - \frac{{\tilde v} 2775 \sin\theta}{\sin\theta_0}+M_{Ro}{\frac{\partial{\pi}}{\partial{x}}} 2776 + \frac{\lambda{\rm Fr}^2 M_{Ro}\cos \theta}{\mu\sin\theta_0} w 2777 = -\frac{{\rm Fr}^2 M_{Ro} {\tilde u} w}{r/H} 2778 + \frac{{\rm Ro} {\bf \hat x}\cdot\nabla^2{\bf u}}{{\rm Re}} 2779 :label: gill_u 2780 2781 .. math:: 2782 {\rm Ro} \frac{D{\tilde v}}{Dt} + 2783 \frac{{\tilde u}\sin\theta}{\sin\theta_0} + M_{Ro}{\frac{\partial{\pi}}{\partial{y}}} ab47de63dc Mart*2784 = -\frac{\mu{\rm Ro} \tan\theta({\tilde u}^2 + {\tilde v}^2)}{r/L} 4f2617d475 Jeff*2785 - \frac{{\rm Fr}^2M_{Ro} {\tilde v} w}{r/H} 2786 + \frac{{\rm Ro} {\bf \hat y}\cdot\nabla^2{\bf u}}{{\rm Re}} 2787 :label: gill_v 2788 2789 .. math:: 2790 {\rm Fr}^2\lambda^2\frac{D{w}}{Dt} - b + {\frac{\partial{\pi}}{\partial{z}}} 2791 -\frac{\lambda\cot \theta_0 {\tilde u}}{M_{Ro}} 2792 = \frac{\lambda\mu^2({\tilde u}^2+{\tilde v}^2)}{M_{Ro}(r/L)} 2793 + \frac{{\rm Fr}^2\lambda^2{\bf \hat z}\cdot\nabla^2{\bf u}}{{\rm Re}} 2794 :label: gill_w 2795 2796 .. math:: 2797 \frac{D{b}}{Dt} + w = \frac{\nabla^2 b}{\Pr{\rm Re}}\nonumber 2798 :label: gill_b 2799 2800 .. math:: 2801 \mu^2\left({\frac{\partial{\tilde u}}{\partial{x}}} + {\frac{\partial{\tilde v}}{\partial{y}}} \right)+{\frac{\partial{w}}{\partial{z}}} = 0 2802 :label: gill_cont 2803 2804 Where 2805 2806 .. math:: 2807 \mu\equiv\frac{\cos\theta_0}{\cos\theta},\ \ \ 2808 {\tilde u}=\frac{u^*}{V\mu},\ \ \ {\tilde v}=\frac{v^*}{V\mu} 2809 2810 .. math:: ab47de63dc Mart*2811 f_0\equiv2\Omega\sin\theta_0,\ \ \ 2812 %,\ \ \ \BFKDt\ \equiv \mu^2\left({\tilde u}\BFKpd x\ 2813 %+{\tilde v} \BFKpd y\ \right)+\frac{\rm Fr^2M_{Ro}}{\rm Ro} w\BFKpd z\ 2814 \frac{D}{Dt} \equiv \mu^2\left({\tilde u}\frac{\partial}{\partial x} 4f2617d475 Jeff*2815 +{\tilde v} \frac{\partial}{\partial y} \right) 2816 +\frac{\rm Fr^2M_{Ro}}{\rm Ro} w\frac{\partial}{\partial z} 2817 2818 .. math:: ab47de63dc Mart*2819 x\equiv \frac{r}{L} \phi \cos \theta_0, \ \ \ 4f2617d475 Jeff*2820 y\equiv \frac{r}{L} \int_{\theta_0}^\theta ab47de63dc Mart*2821 \frac{\cos \theta_0 {\,\rm d\theta}'}{\cos\theta'}, \ \ \ 4f2617d475 Jeff*2822 z\equiv \lambda\frac{r-r_0}{L} 2823 2824 .. math:: t^*=t \frac{L}{V},\ \ \ b^*= b\frac{V f_0M_{Ro}}{\lambda} 2825 2826 .. math:: ab47de63dc Mart*2827 \pi^* = \pi V f_0 LM_{Ro},\ \ \ 4f2617d475 Jeff*2828 w^* = w V \frac{{\rm Fr}^2 \lambda M_{Ro}}{\rm Ro} 2829 2830 .. math:: {\rm Ro} \equiv \frac{V}{f_0 L},\ \ \ M_{Ro}\equiv \max[1,\rm Ro] 2831 2832 .. math:: ab47de63dc Mart*2833 {\rm Fr} \equiv \frac{V}{N \lambda L}, \ \ \ 2834 {\rm Re} \equiv \frac{VL}{\nu}, \ \ \ 4f2617d475 Jeff*2835 {\rm Pr} \equiv \frac{\nu}{\kappa} 2836 2837 Dimensional variables are denoted by an asterisk where necessary. If we 2838 filter over a grid scale typical for ocean models: 2839 2840 | 1m :math:`< L <` 100km 2841 | 0.0001 :math:`< \lambda <` 1 2842 | 0.001m/s :math:`< V <` 1 m/s 2843 | :math:`f_0 <` 0.0001 s :sup:`-1` 2844 | 0.01 s :sup:`-1` :math:`< N <` 0.0001 s :sup:`-1` ab47de63dc Mart*2845 | 4f2617d475 Jeff*2846 2847 these equations are very well approximated by 2848 2849 .. math:: 2850 {\rm Ro}\frac{D{\tilde u}}{Dt} - \frac{{\tilde v} 2851 \sin\theta}{\sin\theta_0}+M_{Ro}{\frac{\partial{\pi}}{\partial{x}}} 2852 = -\frac{\lambda{\rm Fr}^2M_{Ro}\cos \theta}{\mu\sin\theta_0} w 2853 + \frac{{\rm Ro}\nabla^2{{\tilde u}}}{{\rm Re}} \\ 2854 :label: gill_u_filt 2855 2856 .. math:: 2857 {\rm Ro}\frac{D{\tilde v}}{Dt} + 2858 \frac{{\tilde u}\sin\theta}{\sin\theta_0}+M_{Ro}{\frac{\partial{\pi}}{\partial{y}}} 2859 = \frac{{\rm Ro}\nabla^2{{\tilde v}}}{{\rm Re}} \\ 2860 :label: gill_v_filt 2861 2862 .. math:: 2863 {\rm Fr}^2\lambda^2\frac{D{w}}{Dt} - b + {\frac{\partial{\pi}}{\partial{z}}} 2864 = \frac{\lambda\cot \theta_0 {\tilde u}}{M_{Ro}} ab47de63dc Mart*2865 +\frac{{\rm Fr}^2\lambda^2\nabla^2w}{{\rm Re}} 4f2617d475 Jeff*2866 :label: gill_w_filt 2867 2868 .. math:: ab47de63dc Mart*2869 \frac{D{b}}{Dt} + w = \frac{\nabla^2 b}{\Pr{\rm Re}} 4f2617d475 Jeff*2870 :label: gill_b_filt 2871 2872 2873 .. math:: 2874 \mu^2\left({\frac{\partial{\tilde u}}{\partial{x}}} + {\frac{\partial{\tilde v}}{\partial{y}}} \right)+{\frac{\partial{w}}{\partial{z}}} = 0 2875 :label: gill_cont_filt 2876 2877 .. math:: 2878 \nabla^2 \approx \left(\frac{\partial^2}{\partial x^2} 2879 +\frac{\partial^2}{\partial y^2} 2880 +\frac{\partial^2}{\lambda^2\partial z^2}\right) 2881 2882 Neglecting the non-frictional terms on the right-hand side is usually** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 2884.
2883 called the ’traditional’ approximation. It is appropriate, with either 2884 large aspect ratio or far from the tropics. This approximation is used 2885 here, as it does not affect the form of the eddy stresses which is the 2886 main topic. The frictional terms are preserved in this approximate form 2887 for later comparison with eddy stresses.
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated from https://github.com/darwinproject/darwin3 by the 2.3.7-MITgcm-0.1 LXR engine. The LXR team |
|