/[MITgcm]/manual/s_phys_pkgs/text/gmredi.tex
ViewVC logotype

Contents of /manual/s_phys_pkgs/text/gmredi.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.16 - (show annotations) (download) (as text)
Mon Jun 27 02:08:35 2011 UTC (14 years, 1 month ago) by dfer
Branch: MAIN
Changes since 1.15: +82 -45 lines
File MIME type: application/x-tex
Fix a bit the GM/Redi section.

1 \subsection{GMREDI: Gent-McWilliams/Redi SGS Eddy Parameterization}
2 \label{sec:pkg:gmredi}
3 \begin{rawhtml}
4 <!-- CMIREDIR:gmredi: -->
5 \end{rawhtml}
6
7 There are two parts to the Redi/GM parameterization of geostrophic
8 eddies. The first, the Redi scheme \citep{re82}, aims to mix tracer properties along isentropes
9 (neutral surfaces) by means of a diffusion operator oriented along the
10 local isentropic surface. The second part, GM \citep{gen-mcw:90,gen-eta:95}, adiabatically
11 re-arranges tracers through an advective flux where the advecting flow
12 is a function of slope of the isentropic surfaces.
13
14 The first GCM implementation of the Redi scheme was by \cite{Cox87} in the
15 GFDL ocean circulation model. The original approach failed to
16 distinguish between isopycnals and surfaces of locally referenced
17 potential density (now called neutral surfaces) which are proper
18 isentropes for the ocean. As will be discussed later, it also appears
19 that the Cox implementation is susceptible to a computational mode.
20 Due to this mode, the Cox scheme requires a background lateral
21 diffusion to be present to conserve the integrity of the model fields.
22
23 The GM parameterization was then added to the GFDL code in the form of
24 a non-divergent bolus velocity. The method defines two
25 stream-functions expressed in terms of the isoneutral slopes subject
26 to the boundary condition of zero value on upper and lower
27 boundaries. The horizontal bolus velocities are then the vertical
28 derivative of these functions. Here in lies a problem highlighted by
29 \cite{gretal:98}: the bolus velocities involve multiple
30 derivatives on the potential density field, which can consequently
31 give rise to noise. Griffies et al. point out that the GM bolus fluxes
32 can be identically written as a skew flux which involves fewer
33 differential operators. Further, combining the skew flux formulation
34 and Redi scheme, substantial cancellations take place to the point
35 that the horizontal fluxes are unmodified from the lateral diffusion
36 parameterization.
37
38 \subsubsection{Redi scheme: Isopycnal diffusion}
39
40 The Redi scheme diffuses tracers along isopycnals and introduces a
41 term in the tendency (rhs) of such a tracer (here $\tau$) of the form:
42 \begin{equation}
43 \bf{\nabla} \cdot \kappa_\rho \bf{K}_{Redi} \bf{\nabla} \tau
44 \end{equation}
45 where $\kappa_\rho$ is the along isopycnal diffusivity and
46 $\bf{K}_{Redi}$ is a rank 2 tensor that projects the gradient of
47 $\tau$ onto the isopycnal surface. The unapproximated projection tensor is:
48 \begin{equation}
49 \bf{K}_{Redi} = \left(
50 \begin{array}{ccc}
51 1 + S_x& S_x S_y & S_x \\
52 S_x S_y & 1 + S_y & S_y \\
53 S_x & S_y & |S|^2 \\
54 \end{array}
55 \right)
56 \end{equation}
57 Here, $S_x = -\partial_x \sigma / \partial_z \sigma$ and $S_y =
58 -\partial_y \sigma / \partial_z \sigma$ are the components of the
59 isoneutral slope.
60
61 The first point to note is that a typical slope in the ocean interior
62 is small, say of the order $10^{-4}$. A maximum slope might be of
63 order $10^{-2}$ and only exceeds such in unstratified regions where
64 the slope is ill defined. It is therefore justifiable, and
65 customary, to make the small slope approximation, $|S| << 1$. The Redi
66 projection tensor then becomes:
67 \begin{equation}
68 \bf{K}_{Redi} = \left(
69 \begin{array}{ccc}
70 1 & 0 & S_x \\
71 0 & 1 & S_y \\
72 S_x & S_y & |S|^2 \\
73 \end{array}
74 \right)
75 \end{equation}
76
77
78 \subsubsection{GM parameterization}
79
80 The GM parameterization aims to represent the ``advective'' or
81 ``transport'' effect of geostrophic eddies by means of a ``bolus''
82 velocity, $\bf{u}^\star$. The divergence of this advective flux is added
83 to the tracer tendency equation (on the rhs):
84 \begin{equation}
85 - \bf{\nabla} \cdot \tau \bf{u}^\star
86 \end{equation}
87
88 The bolus velocity $\bf{u}^\star$ is defined as the rotational of a
89 streamfunction $\bf{F}^\star$=$(F_x^\star,F_y^\star,0)$:
90 \begin{equation}
91 \bf{u}^\star = \nabla \times \bf{F}^\star =
92 \left( \begin{array}{c}
93 - \partial_z F_y^\star \\
94 + \partial_z F_x^\star \\
95 \partial_x F_y^\star - \partial_y F_x^\star
96 \end{array} \right),
97 \end{equation}
98 and thus is automatically non-divergent. In the GM parameterization, the streamfunction is
99 specified in terms of the isoneutral slopes $S_x$ and $S_y$:
100 \begin{eqnarray}
101 F_x^\star & = & -\kappa_{GM} S_y \\
102 F_y^\star & = & \kappa_{GM} S_x
103 \end{eqnarray}
104 with boundary conditions $F_x^\star=F_y^\star=0$ on upper and lower boundaries.
105 In the end, the bolus transport in the GM parameterization is given by:
106 \begin{equation}
107 \bf{u}^\star = \left(
108 \begin{array}{c}
109 u^\star \\
110 v^\star \\
111 w^\star
112 \end{array}
113 \right) = \left(
114 \begin{array}{c}
115 - \partial_z (\kappa_{GM} S_x) \\
116 - \partial_z (\kappa_{GM} S_y) \\
117 \partial_x (\kappa_{GM} S_x) + \partial_y (\kappa_{GM} S_y)
118 \end{array}
119 \right)
120 \end{equation}
121
122 This is the form of the GM parameterization as applied by Donabasaglu,
123 1997, in MOM versions 1 and 2.
124
125 Note that in the MITgcm, the variables containing the GM bolus streamfunction are:
126 \begin{equation}
127 \left(
128 \begin{array}{c}
129 GM\_PsiX \\
130 GM\_PsiY
131 \end{array}
132 \right) = \left(
133 \begin{array}{c}
134 \kappa_{GM} S_x \\
135 \kappa_{GM} S_y
136 \end{array}
137 \right)= \left(
138 \begin{array}{c}
139 F_y^\star \\
140 -F_x^\star
141 \end{array}
142 \right).
143 \end{equation}
144
145 \subsubsection{Griffies Skew Flux}
146
147 \cite{gr:98} notes that the discretisation of bolus velocities involves
148 multiple layers of differencing and interpolation that potentially
149 lead to noisy fields and computational modes. He pointed out that the
150 bolus flux can be re-written in terms of a non-divergent flux and a
151 skew-flux:
152 \begin{eqnarray*}
153 \bf{u}^\star \tau
154 & = &
155 \left( \begin{array}{c}
156 - \partial_z ( \kappa_{GM} S_x ) \tau \\
157 - \partial_z ( \kappa_{GM} S_y ) \tau \\
158 (\partial_x \kappa_{GM} S_x + \partial_y \kappa_{GM} S_y)\tau
159 \end{array} \right)
160 \\
161 & = &
162 \left( \begin{array}{c}
163 - \partial_z ( \kappa_{GM} S_x \tau) \\
164 - \partial_z ( \kappa_{GM} S_y \tau) \\
165 \partial_x ( \kappa_{GM} S_x \tau) + \partial_y ( \kappa_{GM} S_y \tau)
166 \end{array} \right)
167 + \left( \begin{array}{c}
168 \kappa_{GM} S_x \partial_z \tau \\
169 \kappa_{GM} S_y \partial_z \tau \\
170 - \kappa_{GM} S_x \partial_x \tau - \kappa_{GM} S_y \partial_y \tau
171 \end{array} \right)
172 \end{eqnarray*}
173 The first vector is non-divergent and thus has no effect on the tracer
174 field and can be dropped. The remaining flux can be written:
175 \begin{equation}
176 \bf{u}^\star \tau = - \kappa_{GM} \bf{K}_{GM} \bf{\nabla} \tau
177 \end{equation}
178 where
179 \begin{equation}
180 \bf{K}_{GM} =
181 \left(
182 \begin{array}{ccc}
183 0 & 0 & -S_x \\
184 0 & 0 & -S_y \\
185 S_x & S_y & 0
186 \end{array}
187 \right)
188 \end{equation}
189 is an anti-symmetric tensor.
190
191 This formulation of the GM parameterization involves fewer derivatives
192 than the original and also involves only terms that already appear in
193 the Redi mixing scheme. Indeed, a somewhat fortunate cancellation
194 becomes apparent when we use the GM parameterization in conjunction
195 with the Redi isoneutral mixing scheme:
196 \begin{equation}
197 \kappa_\rho \bf{K}_{Redi} \bf{\nabla} \tau
198 - u^\star \tau =
199 ( \kappa_\rho \bf{K}_{Redi} + \kappa_{GM} \bf{K}_{GM} ) \bf{\nabla} \tau
200 \end{equation}
201 In the instance that $\kappa_{GM} = \kappa_{\rho}$ then
202 \begin{equation}
203 \kappa_\rho \bf{K}_{Redi} + \kappa_{GM} \bf{K}_{GM} =
204 \kappa_\rho
205 \left( \begin{array}{ccc}
206 1 & 0 & 0 \\
207 0 & 1 & 0 \\
208 2 S_x & 2 S_y & |S|^2
209 \end{array}
210 \right)
211 \end{equation}
212 which differs from the variable Laplacian diffusion tensor by only
213 two non-zero elements in the $z$-row.
214
215 \fbox{ \begin{minipage}{4.75in}
216 {\em S/R GMREDI\_CALC\_TENSOR} ({\em pkg/gmredi/gmredi\_calc\_tensor.F})
217
218 $\sigma_x$: {\bf SlopeX} (argument on entry)
219
220 $\sigma_y$: {\bf SlopeY} (argument on entry)
221
222 $\sigma_z$: {\bf SlopeY} (argument)
223
224 $S_x$: {\bf SlopeX} (argument on exit)
225
226 $S_y$: {\bf SlopeY} (argument on exit)
227
228 \end{minipage} }
229
230
231
232 \subsubsection{Variable $\kappa_{GM}$}
233
234 \cite{visbeck:97} suggest making the eddy coefficient,
235 $\kappa_{GM}$, a function of the Eady growth rate,
236 $|f|/\sqrt{Ri}$. The formula involves a non-dimensional constant,
237 $\alpha$, and a length-scale $L$:
238 \begin{displaymath}
239 \kappa_{GM} = \alpha L^2 \overline{ \frac{|f|}{\sqrt{Ri}} }^z
240 \end{displaymath}
241 where the Eady growth rate has been depth averaged (indicated by the
242 over-line). A local Richardson number is defined $Ri = N^2 / (\partial
243 u/\partial z)^2$ which, when combined with thermal wind gives:
244 \begin{displaymath}
245 \frac{1}{Ri} = \frac{(\frac{\partial u}{\partial z})^2}{N^2} =
246 \frac{ ( \frac{g}{f \rho_o} | {\bf \nabla} \sigma | )^2 }{N^2} =
247 \frac{ M^4 }{ |f|^2 N^2 }
248 \end{displaymath}
249 where $M^2$ is defined $M^2 = \frac{g}{\rho_o} |{\bf \nabla} \sigma|$.
250 Substituting into the formula for $\kappa_{GM}$ gives:
251 \begin{displaymath}
252 \kappa_{GM} = \alpha L^2 \overline{ \frac{M^2}{N} }^z =
253 \alpha L^2 \overline{ \frac{M^2}{N^2} N }^z =
254 \alpha L^2 \overline{ |S| N }^z
255 \end{displaymath}
256
257
258 \subsubsection{Tapering and stability}
259
260 Experience with the GFDL model showed that the GM scheme has to be
261 matched to the convective parameterization. This was originally
262 expressed in connection with the introduction of the KPP boundary
263 layer scheme \citep{lar-eta:94} but in fact, as subsequent experience
264 with the MIT model has found, is necessary for any convective
265 parameterization.
266
267 \fbox{ \begin{minipage}{4.75in}
268 {\em S/R GMREDI\_SLOPE\_LIMIT} ({\em
269 pkg/gmredi/gmredi\_slope\_limit.F})
270
271 $\sigma_x, s_x$: {\bf SlopeX} (argument)
272
273 $\sigma_y, s_y$: {\bf SlopeY} (argument)
274
275 $\sigma_z$: {\bf dSigmadRReal} (argument)
276
277 $z_\sigma^{*}$: {\bf dRdSigmaLtd} (argument)
278
279 \end{minipage} }
280
281 \begin{figure}
282 \begin{center}
283 \resizebox{5.0in}{3.0in}{\includegraphics{s_phys_pkgs/figs/tapers.eps}}
284 \end{center}
285 \caption{Taper functions used in GKW91 and DM95.}
286 \label{fig:tapers}
287 \end{figure}
288
289 \begin{figure}
290 \begin{center}
291 \resizebox{5.0in}{3.0in}{\includegraphics{s_phys_pkgs/figs/effective_slopes.eps}}
292 \end{center}
293 \caption{Effective slope as a function of ``true'' slope using Cox
294 slope clipping, GKW91 limiting and DM95 limiting.}
295 \label{fig:effective_slopes}
296 \end{figure}
297
298
299 \subsubsection*{Slope clipping}
300
301 Deep convection sites and the mixed layer are indicated by
302 homogenized, unstable or nearly unstable stratification. The slopes in
303 such regions can be either infinite, very large with a sign reversal
304 or simply very large. From a numerical point of view, large slopes
305 lead to large variations in the tensor elements (implying large bolus
306 flow) and can be numerically unstable. This was first recognized by
307 \cite{Cox87} who implemented ``slope clipping'' in the isopycnal mixing
308 tensor. Here, the slope magnitude is simply restricted by an upper
309 limit:
310 \begin{eqnarray}
311 |\nabla \sigma| & = & \sqrt{ \sigma_x^2 + \sigma_y^2 } \\
312 S_{lim} & = & - \frac{|\nabla \sigma|}{ S_{max} }
313 \;\;\;\;\;\;\;\; \mbox{where $S_{max}$ is a parameter} \\
314 \sigma_z^\star & = & \min( \sigma_z , S_{lim} ) \\
315 {[s_x,s_y]} & = & - \frac{ [\sigma_x,\sigma_y] }{\sigma_z^\star}
316 \end{eqnarray}
317 Notice that this algorithm assumes stable stratification through the
318 ``min'' function. In the case where the fluid is well stratified ($\sigma_z < S_{lim}$) then the slopes evaluate to:
319 \begin{equation}
320 {[s_x,s_y]} = - \frac{ [\sigma_x,\sigma_y] }{\sigma_z}
321 \end{equation}
322 while in the limited regions ($\sigma_z > S_{lim}$) the slopes become:
323 \begin{equation}
324 {[s_x,s_y]} = \frac{ [\sigma_x,\sigma_y] }{|\nabla \sigma|/S_{max}}
325 \end{equation}
326 so that the slope magnitude is limited $\sqrt{s_x^2 + s_y^2} =
327 S_{max}$.
328
329 The slope clipping scheme is activated in the model by setting {\bf
330 GM\_tap\-er\_scheme = 'clipping'} in {\em data.gmredi}.
331
332 Even using slope clipping, it is normally the case that the vertical
333 diffusion term (with coefficient $\kappa_\rho{\bf K}_{33} =
334 \kappa_\rho S_{max}^2$) is large and must be time-stepped using an
335 implicit procedure (see section on discretisation and code later).
336 Fig. \ref{fig-mixedlayer} shows the mixed layer depth resulting from
337 a) using the GM scheme with clipping and b) no GM scheme (horizontal
338 diffusion). The classic result of dramatically reduced mixed layers is
339 evident. Indeed, the deep convection sites to just one or two points
340 each and are much shallower than we might prefer. This, it turns out,
341 is due to the over zealous re-stratification due to the bolus transport
342 parameterization. Limiting the slopes also breaks the adiabatic nature
343 of the GM/Redi parameterization, re-introducing diabatic fluxes in
344 regions where the limiting is in effect.
345
346 \subsubsection*{Tapering: Gerdes, Koberle and Willebrand, Clim. Dyn. 1991}
347
348 The tapering scheme used in \cite{gkw:91}
349 addressed two issues with the clipping method: the introduction of
350 large vertical fluxes in addition to convective adjustment fluxes is
351 avoided by tapering the GM/Redi slopes back to zero in
352 low-stratification regions; the adjustment of slopes is replaced by a
353 tapering of the entire GM/Redi tensor. This means the direction of
354 fluxes is unaffected as the amplitude is scaled.
355
356 The scheme inserts a tapering function, $f_1(S)$, in front of the
357 GM/Redi tensor:
358 \begin{equation}
359 f_1(S) = \min \left[ 1, \left( \frac{S_{max}}{|S|}\right)^2 \right]
360 \end{equation}
361 where $S_{max}$ is the maximum slope you want allowed. Where the
362 slopes, $|S|<S_{max}$ then $f_1(S) = 1$ and the tensor is un-tapered
363 but where $|S| \ge S_{max}$ then $f_1(S)$ scales down the tensor so
364 that the effective vertical diffusivity term $\kappa f_1(S) |S|^2 =
365 \kappa S_{max}^2$.
366
367 The GKW91 tapering scheme is activated in the model by setting {\bf
368 GM\_tap\-er\_scheme = 'gkw91'} in {\em data.gmredi}.
369
370 \subsubsection*{Tapering: Danabasoglu and McWilliams, J. Clim. 1995}
371
372 The tapering scheme used by \cite{dm:95} followed a similar procedure but used a different
373 tapering function, $f_1(S)$:
374 \begin{equation}
375 f_1(S) = \frac{1}{2} \left( 1+\tanh \left[ \frac{S_c - |S|}{S_d} \right] \right)
376 \end{equation}
377 where $S_c = 0.004$ is a cut-off slope and $S_d=0.001$ is a scale over
378 which the slopes are smoothly tapered. Functionally, the operates in
379 the same way as the GKW91 scheme but has a substantially lower
380 cut-off, turning off the GM/Redi SGS parameterization for weaker
381 slopes.
382
383 The DM95 tapering scheme is activated in the model by setting {\bf
384 GM\_tap\-er\_scheme = 'dm95'} in {\em data.gmredi}.
385
386 \subsubsection*{Tapering: Large, Danabasoglu and Doney, JPO 1997}
387
388 The tapering used in \cite{ldd:97} is based on the
389 DM95 tapering scheme, but also tapers the scheme with an additional
390 function of height, $f_2(z)$, so that the GM/Redi SGS fluxes are
391 reduced near the surface:
392 \begin{equation}
393 f_2(z) = \frac{1}{2} \left( 1 + \sin(\pi \frac{z}{D} - \frac{\pi}{2})\right)
394 \end{equation}
395 where $D = L_\rho |S|$ is a depth-scale and $L_\rho=c/f$ with
396 $c=2$~m~s$^{-1}$. This tapering with height was introduced to fix
397 some spurious interaction with the mixed-layer KPP parameterization.
398
399 The LDD97 tapering scheme is activated in the model by setting {\bf
400 GM\_tap\-er\_scheme = 'ldd97'} in {\em data.gmredi}.
401
402
403
404
405 \begin{figure}
406 \begin{center}
407 %\includegraphics{mixedlayer-cox.eps}
408 %\includegraphics{mixedlayer-diff.eps}
409 Figure missing.
410 \end{center}
411 \caption{Mixed layer depth using GM parameterization with a) Cox slope
412 clipping and for comparison b) using horizontal constant diffusion.}
413 \label{fig-mixedlayer}
414 \end{figure}
415
416 \subsubsection{Package Reference}
417 \label{sec:pkg:gmredi:diagnostics}
418
419 {\footnotesize
420 \begin{verbatim}
421 ------------------------------------------------------------------------
422 <-Name->|Levs|<-parsing code->|<-- Units -->|<- Tile (max=80c)
423 ------------------------------------------------------------------------
424 GM_VisbK| 1 |SM P M1 |m^2/s |Mixing coefficient from Visbeck etal parameterization
425 GM_Kux | 15 |UU P 177MR |m^2/s |K_11 element (U.point, X.dir) of GM-Redi tensor
426 GM_Kvy | 15 |VV P 176MR |m^2/s |K_22 element (V.point, Y.dir) of GM-Redi tensor
427 GM_Kuz | 15 |UU 179MR |m^2/s |K_13 element (U.point, Z.dir) of GM-Redi tensor
428 GM_Kvz | 15 |VV 178MR |m^2/s |K_23 element (V.point, Z.dir) of GM-Redi tensor
429 GM_Kwx | 15 |UM 181LR |m^2/s |K_31 element (W.point, X.dir) of GM-Redi tensor
430 GM_Kwy | 15 |VM 180LR |m^2/s |K_32 element (W.point, Y.dir) of GM-Redi tensor
431 GM_Kwz | 15 |WM P LR |m^2/s |K_33 element (W.point, Z.dir) of GM-Redi tensor
432 GM_PsiX | 15 |UU 184LR |m^2/s |GM Bolus transport stream-function : X component
433 GM_PsiY | 15 |VV 183LR |m^2/s |GM Bolus transport stream-function : Y component
434 GM_KuzTz| 15 |UU 186MR |degC.m^3/s |Redi Off-diagonal Tempetature flux: X component
435 GM_KvzTz| 15 |VV 185MR |degC.m^3/s |Redi Off-diagonal Tempetature flux: Y component
436 \end{verbatim}
437 }
438
439 \subsubsection{Experiments and tutorials that use gmredi}
440 \label{sec:pkg:gmredi:experiments}
441
442 \begin{itemize}
443 \item{Global Ocean tutorial, in tutorial\_global\_oce\_latlon verification directory,
444 described in section \ref{sec:eg-global} }
445 \item{ Front Relax experiment, in front\_relax verification directory.}
446 \item{ Ideal 2D Ocean experiment, in ideal\_2D\_oce verification directory.}
447 \end{itemize}
448
449 % DO NOT EDIT HERE
450
451
452

  ViewVC Help
Powered by ViewVC 1.1.22