/[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.11 - (show annotations) (download) (as text)
Wed Apr 5 03:35:15 2006 UTC (19 years, 3 months ago) by molod
Branch: MAIN
Changes since 1.10: +21 -0 lines
File MIME type: application/x-tex
Add list of package diagnostics as a subsubsection (working on other packages too)

1 \subsection{GMREDI: Gent/McWiliams/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 aims to mix tracer properties along isentropes
9 (neutral surfaces) by means of a diffusion operator oriented along the
10 local isentropic surface (Redi). The second part, adiabatically
11 re-arranges tracers through an advective flux where the advecting flow
12 is a function of slope of the isentropic surfaces (GM).
13
14 The first GCM implementation of the Redi scheme was by Cox 1987 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 Griffies et al., 1997: 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 parameterise the ``advective'' or
81 ``transport'' effect of geostrophic eddies by means of a ``bolus''
82 velocity, $\bf{u}^*$. 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}^*
86 \end{equation}
87
88 The bolus velocity is defined as:
89 \begin{eqnarray}
90 u^* & = & - \partial_z F_x \\
91 v^* & = & - \partial_z F_y \\
92 w^* & = & \partial_x F_x + \partial_y F_y
93 \end{eqnarray}
94 where $F_x$ and $F_y$ are stream-functions with boundary conditions
95 $F_x=F_y=0$ on upper and lower boundaries. The virtue of casting the
96 bolus velocity in terms of these stream-functions is that they are
97 automatically non-divergent ($\partial_x u^* + \partial_y v^* = -
98 \partial_{xz} F_x - \partial_{yz} F_y = - \partial_z w^*$). $F_x$ and
99 $F_y$ are specified in terms of the isoneutral slopes $S_x$ and $S_y$:
100 \begin{eqnarray}
101 F_x & = & \kappa_{GM} S_x \\
102 F_y & = & \kappa_{GM} S_y
103 \end{eqnarray}
104 This is the form of the GM parameterization as applied by Donabasaglu,
105 1997, in MOM versions 1 and 2.
106
107 \subsubsection{Griffies Skew Flux}
108
109 Griffies notes that the discretisation of bolus velocities involves
110 multiple layers of differencing and interpolation that potentially
111 lead to noisy fields and computational modes. He pointed out that the
112 bolus flux can be re-written in terms of a non-divergent flux and a
113 skew-flux:
114 \begin{eqnarray*}
115 \bf{u}^* \tau
116 & = &
117 \left( \begin{array}{c}
118 - \partial_z ( \kappa_{GM} S_x ) \tau \\
119 - \partial_z ( \kappa_{GM} S_y ) \tau \\
120 (\partial_x \kappa_{GM} S_x + \partial_y \kappa_{GM} S_y)\tau
121 \end{array} \right)
122 \\
123 & = &
124 \left( \begin{array}{c}
125 - \partial_z ( \kappa_{GM} S_x \tau) \\
126 - \partial_z ( \kappa_{GM} S_y \tau) \\
127 \partial_x ( \kappa_{GM} S_x \tau) + \partial_y ( \kappa_{GM} S_y) \tau)
128 \end{array} \right)
129 + \left( \begin{array}{c}
130 \kappa_{GM} S_x \partial_z \tau \\
131 \kappa_{GM} S_y \partial_z \tau \\
132 - \kappa_{GM} S_x \partial_x \tau - \kappa_{GM} S_y) \partial_y \tau
133 \end{array} \right)
134 \end{eqnarray*}
135 The first vector is non-divergent and thus has no effect on the tracer
136 field and can be dropped. The remaining flux can be written:
137 \begin{equation}
138 \bf{u}^* \tau = - \kappa_{GM} \bf{K}_{GM} \bf{\nabla} \tau
139 \end{equation}
140 where
141 \begin{equation}
142 \bf{K}_{GM} =
143 \left(
144 \begin{array}{ccc}
145 0 & 0 & -S_x \\
146 0 & 0 & -S_y \\
147 S_x & S_y & 0
148 \end{array}
149 \right)
150 \end{equation}
151 is an anti-symmetric tensor.
152
153 This formulation of the GM parameterization involves fewer derivatives
154 than the original and also involves only terms that already appear in
155 the Redi mixing scheme. Indeed, a somewhat fortunate cancellation
156 becomes apparent when we use the GM parameterization in conjunction
157 with the Redi isoneutral mixing scheme:
158 \begin{equation}
159 \kappa_\rho \bf{K}_{Redi} \bf{\nabla} \tau
160 - u^* \tau =
161 ( \kappa_\rho \bf{K}_{Redi} + \kappa_{GM} \bf{K}_{GM} ) \bf{\nabla} \tau
162 \end{equation}
163 In the instance that $\kappa_{GM} = \kappa_{\rho}$ then
164 \begin{equation}
165 \kappa_\rho \bf{K}_{Redi} + \kappa_{GM} \bf{K}_{GM} =
166 \kappa_\rho
167 \left( \begin{array}{ccc}
168 1 & 0 & 0 \\
169 0 & 1 & 0 \\
170 2 S_x & 2 S_y & |S|^2
171 \end{array}
172 \right)
173 \end{equation}
174 which differs from the variable Laplacian diffusion tensor by only
175 two non-zero elements in the $z$-row.
176
177 \fbox{ \begin{minipage}{4.75in}
178 {\em S/R GMREDI\_CALC\_TENSOR} ({\em pkg/gmredi/gmredi\_calc\_tensor.F})
179
180 $\sigma_x$: {\bf SlopeX} (argument on entry)
181
182 $\sigma_y$: {\bf SlopeY} (argument on entry)
183
184 $\sigma_z$: {\bf SlopeY} (argument)
185
186 $S_x$: {\bf SlopeX} (argument on exit)
187
188 $S_y$: {\bf SlopeY} (argument on exit)
189
190 \end{minipage} }
191
192
193
194 \subsubsection{Variable $\kappa_{GM}$}
195
196 Visbeck et al., 1996, suggest making the eddy coefficient,
197 $\kappa_{GM}$, a function of the Eady growth rate,
198 $|f|/\sqrt{Ri}$. The formula involves a non-dimensional constant,
199 $\alpha$, and a length-scale $L$:
200 \begin{displaymath}
201 \kappa_{GM} = \alpha L^2 \overline{ \frac{|f|}{\sqrt{Ri}} }^z
202 \end{displaymath}
203 where the Eady growth rate has been depth averaged (indicated by the
204 over-line). A local Richardson number is defined $Ri = N^2 / (\partial
205 u/\partial z)^2$ which, when combined with thermal wind gives:
206 \begin{displaymath}
207 \frac{1}{Ri} = \frac{(\frac{\partial u}{\partial z})^2}{N^2} =
208 \frac{ ( \frac{g}{f \rho_o} | {\bf \nabla} \sigma | )^2 }{N^2} =
209 \frac{ M^4 }{ |f|^2 N^2 }
210 \end{displaymath}
211 where $M^2$ is defined $M^2 = \frac{g}{\rho_o} |{\bf \nabla} \sigma|$.
212 Substituting into the formula for $\kappa_{GM}$ gives:
213 \begin{displaymath}
214 \kappa_{GM} = \alpha L^2 \overline{ \frac{M^2}{N} }^z =
215 \alpha L^2 \overline{ \frac{M^2}{N^2} N }^z =
216 \alpha L^2 \overline{ |S| N }^z
217 \end{displaymath}
218
219
220 \subsubsection{Tapering and stability}
221
222 Experience with the GFDL model showed that the GM scheme has to be
223 matched to the convective parameterization. This was originally
224 expressed in connection with the introduction of the KPP boundary
225 layer scheme (Large et al., 97) but in fact, as subsequent experience
226 with the MIT model has found, is necessary for any convective
227 parameterization.
228
229 \fbox{ \begin{minipage}{4.75in}
230 {\em S/R GMREDI\_SLOPE\_LIMIT} ({\em
231 pkg/gmredi/gmredi\_slope\_limit.F})
232
233 $\sigma_x, s_x$: {\bf SlopeX} (argument)
234
235 $\sigma_y, s_y$: {\bf SlopeY} (argument)
236
237 $\sigma_z$: {\bf dSigmadRReal} (argument)
238
239 $z_\sigma^{*}$: {\bf dRdSigmaLtd} (argument)
240
241 \end{minipage} }
242
243 \begin{figure}
244 \begin{center}
245 \resizebox{5.0in}{3.0in}{\includegraphics{part6/tapers.eps}}
246 \end{center}
247 \caption{Taper functions used in GKW99 and DM95.}
248 \label{fig:tapers}
249 \end{figure}
250
251 \begin{figure}
252 \begin{center}
253 \resizebox{5.0in}{3.0in}{\includegraphics{part6/effective_slopes.eps}}
254 \end{center}
255 \caption{Effective slope as a function of ``true'' slope using Cox
256 slope clipping, GKW91 limiting and DM95 limiting.}
257 \label{fig:effective_slopes}
258 \end{figure}
259
260
261 Slope clipping:
262
263 Deep convection sites and the mixed layer are indicated by
264 homogenized, unstable or nearly unstable stratification. The slopes in
265 such regions can be either infinite, very large with a sign reversal
266 or simply very large. From a numerical point of view, large slopes
267 lead to large variations in the tensor elements (implying large bolus
268 flow) and can be numerically unstable. This was first recognized by
269 Cox, 1987, who implemented ``slope clipping'' in the isopycnal mixing
270 tensor. Here, the slope magnitude is simply restricted by an upper
271 limit:
272 \begin{eqnarray}
273 |\nabla \sigma| & = & \sqrt{ \sigma_x^2 + \sigma_y^2 } \\
274 S_{lim} & = & - \frac{|\nabla \sigma|}{ S_{max} }
275 \;\;\;\;\;\;\;\; \mbox{where $S_{max}$ is a parameter} \\
276 \sigma_z^\star & = & \min( \sigma_z , S_{lim} ) \\
277 {[s_x,s_y]} & = & - \frac{ [\sigma_x,\sigma_y] }{\sigma_z^\star}
278 \end{eqnarray}
279 Notice that this algorithm assumes stable stratification through the
280 ``min'' function. In the case where the fluid is well stratified ($\sigma_z < S_{lim}$) then the slopes evaluate to:
281 \begin{equation}
282 {[s_x,s_y]} = - \frac{ [\sigma_x,\sigma_y] }{\sigma_z}
283 \end{equation}
284 while in the limited regions ($\sigma_z > S_{lim}$) the slopes become:
285 \begin{equation}
286 {[s_x,s_y]} = \frac{ [\sigma_x,\sigma_y] }{|\nabla \sigma|/S_{max}}
287 \end{equation}
288 so that the slope magnitude is limited $\sqrt{s_x^2 + s_y^2} =
289 S_{max}$.
290
291 The slope clipping scheme is activated in the model by setting {\bf
292 GM\_tap\-er\_scheme = 'clipping'} in {\em data.gmredi}.
293
294 Even using slope clipping, it is normally the case that the vertical
295 diffusion term (with coefficient $\kappa_\rho{\bf K}_{33} =
296 \kappa_\rho S_{max}^2$) is large and must be time-stepped using an
297 implicit procedure (see section on discretisation and code later).
298 Fig. \ref{fig-mixedlayer} shows the mixed layer depth resulting from
299 a) using the GM scheme with clipping and b) no GM scheme (horizontal
300 diffusion). The classic result of dramatically reduced mixed layers is
301 evident. Indeed, the deep convection sites to just one or two points
302 each and are much shallower than we might prefer. This, it turns out,
303 is due to the over zealous re-stratification due to the bolus transport
304 parameterization. Limiting the slopes also breaks the adiabatic nature
305 of the GM/Redi parameterization, re-introducing diabatic fluxes in
306 regions where the limiting is in effect.
307
308 Tapering: Gerdes, Koberle and Willebrand, Clim. Dyn. 1991:
309
310 The tapering scheme used in Gerdes et al., 1999, (\cite{gkw:99})
311 addressed two issues with the clipping method: the introduction of
312 large vertical fluxes in addition to convective adjustment fluxes is
313 avoided by tapering the GM/Redi slopes back to zero in
314 low-stratification regions; the adjustment of slopes is replaced by a
315 tapering of the entire GM/Redi tensor. This means the direction of
316 fluxes is unaffected as the amplitude is scaled.
317
318 The scheme inserts a tapering function, $f_1(S)$, in front of the
319 GM/Redi tensor:
320 \begin{equation}
321 f_1(S) = \min \left[ 1, \left( \frac{S_{max}}{|S|}\right)^2 \right]
322 \end{equation}
323 where $S_{max}$ is the maximum slope you want allowed. Where the
324 slopes, $|S|<S_{max}$ then $f_1(S) = 1$ and the tensor is un-tapered
325 but where $|S| \ge S_{max}$ then $f_1(S)$ scales down the tensor so
326 that the effective vertical diffusivity term $\kappa f_1(S) |S|^2 =
327 \kappa S_{max}^2$.
328
329 The GKW tapering scheme is activated in the model by setting {\bf
330 GM\_tap\-er\_scheme = 'gkw91'} in {\em data.gmredi}.
331
332 \subsubsection{Tapering: Danabasoglu and McWilliams, J. Clim. 1995}
333
334 The tapering scheme used by Danabasoglu and McWilliams, 1995,
335 \cite{dm:95}, followed a similar procedure but used a different
336 tapering function, $f_1(S)$:
337 \begin{equation}
338 f_1(S) = \frac{1}{2} \left( 1+\tanh \left[ \frac{S_c - |S|}{S_d} \right] \right)
339 \end{equation}
340 where $S_c = 0.004$ is a cut-off slope and $S_d=0.001$ is a scale over
341 which the slopes are smoothly tapered. Functionally, the operates in
342 the same way as the GKW91 scheme but has a substantially lower
343 cut-off, turning off the GM/Redi SGS parameterization for weaker
344 slopes.
345
346 The DM tapering scheme is activated in the model by setting {\bf
347 GM\_tap\-er\_scheme = 'dm95'} in {\em data.gmredi}.
348
349 \subsubsection{Tapering: Large, Danabasoglu and Doney, JPO 1997}
350
351 The tapering used in Large et al., 1997, \cite{ldd:97}, is based on the
352 DM95 tapering scheme, but also tapers the scheme with an additional
353 function of height, $f_2(z)$, so that the GM/Redi SGS fluxes are
354 reduced near the surface:
355 \begin{equation}
356 f_2(S) = \frac{1}{2} \left( 1 + \sin(\pi \frac{z}{D} - \pi/2)\right)
357 \end{equation}
358 where $D = L_\rho |S|$ is a depth-scale and $L_\rho=c/f$ with
359 $c=2$~m~s$^{-1}$. This tapering with height was introduced to fix
360 some spurious interaction with the mixed-layer KPP parameterization.
361
362 The LDD tapering scheme is activated in the model by setting {\bf
363 GM\_tap\-er\_scheme = 'ldd97'} in {\em data.gmredi}.
364
365
366
367
368 \begin{figure}
369 \begin{center}
370 %\includegraphics{mixedlayer-cox.eps}
371 %\includegraphics{mixedlayer-diff.eps}
372 Figure missing.
373 \end{center}
374 \caption{Mixed layer depth using GM parameterization with a) Cox slope
375 clipping and for comparison b) using horizontal constant diffusion.}
376 \label{fig-mixedlayer}
377 \end{figure}
378
379 \subsubsection{Package Reference}
380 \label{sec:pkg:gmredi:diagnostics}
381
382 \begin{verbatim}
383 ------------------------------------------------------------------------
384 <-Name->|Levs|<-parsing code->|<-- Units -->|<- Tile (max=80c)
385 ------------------------------------------------------------------------
386 GM_VisbK| 1 |SM P M1 |m^2/s |Mixing coefficient from Visbeck etal parameterization
387 GM_Kux | 15 |UU P 177MR |m^2/s |K_11 element (U.point, X.dir) of GM-Redi tensor
388 GM_Kvy | 15 |VV P 176MR |m^2/s |K_22 element (V.point, Y.dir) of GM-Redi tensor
389 GM_Kuz | 15 |UU 179MR |m^2/s |K_13 element (U.point, Z.dir) of GM-Redi tensor
390 GM_Kvz | 15 |VV 178MR |m^2/s |K_23 element (V.point, Z.dir) of GM-Redi tensor
391 GM_Kwx | 15 |UM 181LR |m^2/s |K_31 element (W.point, X.dir) of GM-Redi tensor
392 GM_Kwy | 15 |VM 180LR |m^2/s |K_32 element (W.point, Y.dir) of GM-Redi tensor
393 GM_Kwz | 15 |WM P LR |m^2/s |K_33 element (W.point, Z.dir) of GM-Redi tensor
394 GM_PsiX | 15 |UU 184LR |m^2/s |GM Bolus transport stream-function : X component
395 GM_PsiY | 15 |VV 183LR |m^2/s |GM Bolus transport stream-function : Y component
396 GM_KuzTz| 15 |UU 186MR |degC.m^3/s |Redi Off-diagonal Tempetature flux: X component
397 GM_KvzTz| 15 |VV 185MR |degC.m^3/s |Redi Off-diagonal Tempetature flux: Y component
398 \end{verbatim}
399
400 \subsubsection{Package Reference}
401 % DO NOT EDIT HERE
402
403
404

  ViewVC Help
Powered by ViewVC 1.1.22