| 1 |
\documentclass[12pt]{article} |
| 2 |
|
| 3 |
\pagestyle{myheadings} |
| 4 |
\markright{Adcroft \& Hill \hspace*{3mm} |
| 5 |
Redi/GM/Griffies scheme in the MIT OGCM} |
| 6 |
|
| 7 |
\usepackage{html} |
| 8 |
|
| 9 |
\begin{document} |
| 10 |
|
| 11 |
\bodytext{bgcolor="#FFFFFFFF"} |
| 12 |
|
| 13 |
\begin{center} |
| 14 |
{\Large \bf Implementation of the GM/Redi/Griffies schemes in the MIT |
| 15 |
OGCM} |
| 16 |
|
| 17 |
\vspace*{4mm} |
| 18 |
{\large A.Adcroft and C.Hill, MIT} |
| 19 |
|
| 20 |
\vspace*{3mm} |
| 21 |
{\large July 2000} |
| 22 |
\end{center} |
| 23 |
|
| 24 |
\section{Overview} |
| 25 |
|
| 26 |
There are two parts to the Redi/GM parameterization of geostrophic |
| 27 |
eddies. The first aims to mix tracer properties along isentropes |
| 28 |
(neutral surfaces) by means of a diffusion operator oriented along the |
| 29 |
local isentropic surface (Redi). The second part, adiabatically |
| 30 |
re-arranges tracers through an advective flux where the advecting flow |
| 31 |
is a function of slope of the isentropic surfaces (GM). |
| 32 |
|
| 33 |
The first GCM implementation of the Redi scheme was by Cox 1987 in the |
| 34 |
GFDL ocean circulation model. The original approach failed to |
| 35 |
distinguish between isopycnals and surfaces of locally referenced |
| 36 |
potential density (now called neutral surfaces) which are proper |
| 37 |
isentropes for the ocean. As will be discussed later, it also appears |
| 38 |
that the Cox implementation is susceptible to a computational mode. |
| 39 |
Due to this mode, the Cox scheme requires a background lateral |
| 40 |
diffusion to be present to conserve the integrity of the model fields. |
| 41 |
|
| 42 |
The GM parameterization was then added to the GFDL code in the form of |
| 43 |
a non-divergent bolus velocity. The method defines two |
| 44 |
stream-functions expressed in terms of the isoneutral slopes subject |
| 45 |
to the boundary condition of zero value on upper and lower |
| 46 |
boundaries. The horizontal bolus velocities are then the vertical |
| 47 |
derivative of these functions. Here in lies a problem highlighted by |
| 48 |
Griffies et al., 1997: the bolus velocities involve multiple |
| 49 |
derivatives on the potential density field, which can consequently |
| 50 |
give rise to noise. Griffies et al. point out that the GM bolus fluxes |
| 51 |
can be identically written as a skew flux which involves fewer |
| 52 |
differential operators. Further, combining the skew flux formulation |
| 53 |
and Redi scheme, substantial cancellations take place to the point |
| 54 |
that the horizontal fluxes are unmodified from the lateral diffusion |
| 55 |
parameterization. |
| 56 |
|
| 57 |
\section{Redi scheme: Isopycnal diffusion} |
| 58 |
|
| 59 |
The Redi scheme diffuses tracers along isopycnals and introduces a |
| 60 |
term in the tendency (rhs) of such a tracer (here $\tau$) of the form: |
| 61 |
\begin{equation} |
| 62 |
\bf{\nabla} \cdot \kappa_\rho \bf{K}_{Redi} \bf{\nabla} \tau |
| 63 |
\end{equation} |
| 64 |
where $\kappa_\rho$ is the along isopycnal diffusivity and |
| 65 |
$\bf{K}_{Redi}$ is a rank 2 tensor that projects the gradient of |
| 66 |
$\tau$ onto the isopycnal surface. The unapproximated projection tensor is: |
| 67 |
\begin{equation} |
| 68 |
\bf{K}_{Redi} = \left( |
| 69 |
\begin{array}{ccc} |
| 70 |
1 + S_x& S_x S_y & S_x \\ |
| 71 |
S_x S_y & 1 + S_y & S_y \\ |
| 72 |
S_x & S_y & |S|^2 \\ |
| 73 |
\end{array} |
| 74 |
\right) |
| 75 |
\end{equation} |
| 76 |
Here, $S_x = -\partial_x \sigma / \partial_z \sigma$ and $S_y = |
| 77 |
-\partial_y \sigma / \partial_z \sigma$ are the components of the |
| 78 |
isoneutral slope. |
| 79 |
|
| 80 |
The first point to note is that a typical slope in the ocean interior |
| 81 |
is small, say of the order $10^{-4}$. A maximum slope might be of |
| 82 |
order $10^{-2}$ and only exceeds such in unstratified regions where |
| 83 |
the slope is ill defined. It is therefore justifiable, and |
| 84 |
customary, to make the small slope approximation, $|S| << 1$. The Redi |
| 85 |
projection tensor then becomes: |
| 86 |
\begin{equation} |
| 87 |
\bf{K}_{Redi} = \left( |
| 88 |
\begin{array}{ccc} |
| 89 |
1 & 0 & S_x \\ |
| 90 |
0 & 1 & S_y \\ |
| 91 |
S_x & S_y & |S|^2 \\ |
| 92 |
\end{array} |
| 93 |
\right) |
| 94 |
\end{equation} |
| 95 |
|
| 96 |
|
| 97 |
\section{GM parameterization} |
| 98 |
|
| 99 |
The GM parameterization aims to parameterise the ``advective'' or |
| 100 |
``transport'' effect of geostrophic eddies by means of a ``bolus'' |
| 101 |
velocity, $\bf{u}^*$. The divergence of this advective flux is added |
| 102 |
to the tracer tendency equation (on the rhs): |
| 103 |
\begin{equation} |
| 104 |
- \bf{\nabla} \cdot \tau \bf{u}^* |
| 105 |
\end{equation} |
| 106 |
|
| 107 |
The bolus velocity is defined as: |
| 108 |
\begin{eqnarray} |
| 109 |
u^* & = & - \partial_z F_x \\ |
| 110 |
v^* & = & - \partial_z F_y \\ |
| 111 |
w^* & = & \partial_x F_x + \partial_y F_y |
| 112 |
\end{eqnarray} |
| 113 |
where $F_x$ and $F_y$ are stream-functions with boundary conditions |
| 114 |
$F_x=F_y=0$ on upper and lower boundaries. The virtue of casting the |
| 115 |
bolus velocity in terms of these stream-functions is that they are |
| 116 |
automatically non-divergent ($\partial_x u^* + \partial_y v^* = - |
| 117 |
\partial_{xz} F_x - \partial_{yz} F_y = - \partial_z w^*$). $F_x$ and |
| 118 |
$F_y$ are specified in terms of the isoneutral slopes $S_x$ and $S_y$: |
| 119 |
\begin{eqnarray} |
| 120 |
F_x & = & \kappa_{GM} S_x \\ |
| 121 |
F_y & = & \kappa_{GM} S_y |
| 122 |
\end{eqnarray} |
| 123 |
This is the form of the GM parameterization as applied by Donabasaglu, |
| 124 |
1997, in MOM versions 1 and 2. |
| 125 |
|
| 126 |
\subsection{Griffies Skew Flux} |
| 127 |
|
| 128 |
Griffies notes that the discretisation of bolus velocities involves |
| 129 |
multiple layers of differencing and interpolation that potentially |
| 130 |
lead to noisy fields and computational modes. He pointed out that the |
| 131 |
bolus flux can be re-written in terms of a non-divergent flux and a |
| 132 |
skew-flux: |
| 133 |
\begin{eqnarray*} |
| 134 |
\bf{u}^* \tau |
| 135 |
& = & |
| 136 |
\left( \begin{array}{c} |
| 137 |
- \partial_z ( \kappa_{GM} S_x ) \tau \\ |
| 138 |
- \partial_z ( \kappa_{GM} S_y ) \tau \\ |
| 139 |
(\partial_x \kappa_{GM} S_x + \partial_y \kappa_{GM} S_y)\tau |
| 140 |
\end{array} \right) |
| 141 |
\\ |
| 142 |
& = & |
| 143 |
\left( \begin{array}{c} |
| 144 |
- \partial_z ( \kappa_{GM} S_x \tau) \\ |
| 145 |
- \partial_z ( \kappa_{GM} S_y \tau) \\ |
| 146 |
\partial_x ( \kappa_{GM} S_x \tau) + \partial_y ( \kappa_{GM} S_y) \tau) |
| 147 |
\end{array} \right) |
| 148 |
+ \left( \begin{array}{c} |
| 149 |
\kappa_{GM} S_x \partial_z \tau \\ |
| 150 |
\kappa_{GM} S_y \partial_z \tau \\ |
| 151 |
- \kappa_{GM} S_x \partial_x \tau - \kappa_{GM} S_y) \partial_y \tau |
| 152 |
\end{array} \right) |
| 153 |
\end{eqnarray*} |
| 154 |
The first vector is non-divergent and thus has no effect on the tracer |
| 155 |
field and can be dropped. The remaining flux can be written: |
| 156 |
\begin{equation} |
| 157 |
\bf{u}^* \tau = - \kappa_{GM} \bf{K}_{GM} \bf{\nabla} \tau |
| 158 |
\end{equation} |
| 159 |
where |
| 160 |
\begin{equation} |
| 161 |
\bf{K}_{GM} = |
| 162 |
\left( |
| 163 |
\begin{array}{ccc} |
| 164 |
0 & 0 & -S_x \\ |
| 165 |
0 & 0 & -S_y \\ |
| 166 |
S_x & S_y & 0 |
| 167 |
\end{array} |
| 168 |
\right) |
| 169 |
\end{equation} |
| 170 |
is an anti-symmetric tensor. |
| 171 |
|
| 172 |
This formulation of the GM parameterization involves fewer derivatives |
| 173 |
than the original and also involves only terms that already appear in |
| 174 |
the Redi mixing scheme. Indeed, a somewhat fortunate cancellation |
| 175 |
becomes apparent when we use the GM parameterization in conjunction |
| 176 |
with the Redi isoneutral mixing scheme: |
| 177 |
\begin{equation} |
| 178 |
\kappa_\rho \bf{K}_{Redi} \bf{\nabla} \tau |
| 179 |
- u^* \tau = |
| 180 |
( \kappa_\rho \bf{K}_{Redi} + \kappa_{GM} \bf{K}_{GM} ) \bf{\nabla} \tau |
| 181 |
\end{equation} |
| 182 |
In the instance that $\kappa_{GM} = \kappa_{\rho}$ then |
| 183 |
\begin{equation} |
| 184 |
\kappa_\rho \bf{K}_{Redi} + \kappa_{GM} \bf{K}_{GM} = |
| 185 |
\kappa_\rho |
| 186 |
\left( \begin{array}{ccc} |
| 187 |
1 & 0 & 0 \\ |
| 188 |
0 & 1 & 0 \\ |
| 189 |
2 S_x & 2 S_y & |S|^2 |
| 190 |
\end{array} |
| 191 |
\right) |
| 192 |
\end{equation} |
| 193 |
which differs from the variable laplacian diffusion tensor by only |
| 194 |
two non-zero elements in the $z$-row. |
| 195 |
|
| 196 |
\section{Variable $\kappa_{GM}$} |
| 197 |
|
| 198 |
Visbeck et al., 1996, suggest making the eddy coefficient, |
| 199 |
$\kappa_{GM}$, a function of the Eady growth rate, |
| 200 |
$|f|/\sqrt{Ri}$. The formula involves a non-dimensional constant, |
| 201 |
$\alpha$, and a length-scale $L$: |
| 202 |
\begin{displaymath} |
| 203 |
\kappa_{GM} = \alpha L^2 \overline{ \frac{|f|}{\sqrt{Ri}} }^z |
| 204 |
\end{displaymath} |
| 205 |
where the Eady growth rate has been depth averaged (indicated by the |
| 206 |
over-line). A local Richardson number is defined $Ri = N^2 / (\partial |
| 207 |
u/\partial z)^2$ which, when combined with thermal wind gives: |
| 208 |
\begin{displaymath} |
| 209 |
\frac{1}{Ri} = \frac{(\frac{\partial u}{\partial z})^2}{N^2} = |
| 210 |
\frac{ ( \frac{g}{f \rho_o} | {\bf \nabla} \sigma | )^2 }{N^2} = |
| 211 |
\frac{ M^4 }{ |f|^2 N^2 } |
| 212 |
\end{displaymath} |
| 213 |
where $M^2$ is defined $M^2 = \frac{g}{\rho_o} |{\bf \nabla} \sigma|$. |
| 214 |
Substituting into the formula for $\kappa_{GM}$ gives: |
| 215 |
\begin{displaymath} |
| 216 |
\kappa_{GM} = \alpha L^2 \overline{ \frac{M^2}{N} }^z = |
| 217 |
\alpha L^2 \overline{ \frac{M^2}{N^2} N }^z = |
| 218 |
\alpha L^2 \overline{ |S| N }^z |
| 219 |
\end{displaymath} |
| 220 |
|
| 221 |
|
| 222 |
\section{Tapering and stability} |
| 223 |
|
| 224 |
Experience with the GFDL model showed that the GM scheme has to be |
| 225 |
matched to the convective parameterization. This was originally |
| 226 |
expressed in connection with the introduction of the KPP boundary |
| 227 |
layer scheme (Large et al., 97) but infact, as subsequent experience |
| 228 |
with the MIT model has found, is necessary for any convective |
| 229 |
parameterization. |
| 230 |
|
| 231 |
Deep convection sites and the mixed layer are indicated by |
| 232 |
homogenized, unstable or nearly unstable stratification. The slopes |
| 233 |
become in such regions are infinite, very large with a sign reversal |
| 234 |
or simply very large. From a numerical point of view, large slopes |
| 235 |
lead to large variations in the tensor elements (implying large bolus |
| 236 |
flow) and can be numerically unstable. This was first reognized by |
| 237 |
Cox, 1987, who implemented ``slopes clipping'' in the isopycnal mixing |
| 238 |
tensor. Here, the slope magnitude is simply restricted by an upper limit: |
| 239 |
\begin{eqnarray} |
| 240 |
|\nabla \sigma| & = & \sqrt{ \sigma_x^2 + \sigma_y^2 } \\ |
| 241 |
S_{lim} & = & - \frac{|\nabla \sigma|}{ S_{max} } |
| 242 |
\;\;\;\;\;\;\;\; \mbox{where $S_{max}$ is a parameter} \\ |
| 243 |
\sigma_z^\star & = & min( \sigma_z , S_{lim} ) \\ |
| 244 |
{[s_x,s_y]} & = & - \frac{ [\sigma_x,\sigma_y] }{\sigma_z^\star} |
| 245 |
\end{eqnarray} |
| 246 |
Notice that this algorithm assumes stable stratification through the |
| 247 |
``min'' function. In the case where the fluid is well stratified ($\sigma_z < S_{lim}$) then the slopes evaluate to: |
| 248 |
\begin{equation} |
| 249 |
{[s_x,s_y]} = - \frac{ [\sigma_x,\sigma_y] }{\sigma_z} |
| 250 |
\end{equation} |
| 251 |
while in the limited regions ($\sigma_z > S_{lim}$) the slopes become: |
| 252 |
\begin{equation} |
| 253 |
{[s_x,s_y]} = \frac{ [\sigma_x,\sigma_y] }{|\nabla \sigma|/S_{max}} |
| 254 |
\end{equation} |
| 255 |
so that the slope magnitude is limited $\sqrt{s_x^2 + s_y^2} = |
| 256 |
S_{max}$. |
| 257 |
|
| 258 |
Even using slope clipping, it is normally the case that the vertical |
| 259 |
diffusion term (with coefficient $\kappa_\rho{\bf K}_{33} = \kappa_\rho |
| 260 |
S_{max}^2$) is large and must be time-stepped using an implicit |
| 261 |
procedure (see section on discretisation and code later). |
| 262 |
Fig. \ref{fig-mixedlayer} shows the mixed layer depth resulting from |
| 263 |
a) using the GM scheme with clipping and b) no GM scheme (horizontal |
| 264 |
diffusion). The classic result of dramatically reducd mixed layers is |
| 265 |
evident. Indeed, the deep convection sites to just one or two points |
| 266 |
each and are much shallower than we might prefer. |
| 267 |
|
| 268 |
This it turns out is due to the over zealous restratification due to |
| 269 |
the bolus transport parameterization. |
| 270 |
|
| 271 |
Gerdes et al., 1991 .... |
| 272 |
|
| 273 |
Danabasoglu and McWilliams, 1995 .... |
| 274 |
|
| 275 |
Large et al., 1997 .... |
| 276 |
|
| 277 |
Issue: should GM and Redi be tapered? cf convective paper |
| 278 |
|
| 279 |
Issue: is adiabtic important in these regions? changes discretization |
| 280 |
|
| 281 |
|
| 282 |
|
| 283 |
\begin{figure} |
| 284 |
\include{mixedlayer-cox.eps} |
| 285 |
\include{mixedlayer-diff.eps} |
| 286 |
\caption{Mixed layer depth using GM parameterization with a) Cox slope |
| 287 |
clipping and for comparison b) using horizontal constant diffusion.} |
| 288 |
\ref{fig-mixedlayer} |
| 289 |
\end{figure} |
| 290 |
|
| 291 |
\begin{figure} |
| 292 |
\include{slopelimits.eps} |
| 293 |
\caption{Effective slope as a function of ``true'' slope using a) Cox |
| 294 |
slope clipping, b) GKW91 limiting, c) DM95 limiting and d) LDD97 |
| 295 |
limiting.} |
| 296 |
\end{figure} |
| 297 |
|
| 298 |
|
| 299 |
\begin{figure} |
| 300 |
%\includegraphics{coxslope.eps} |
| 301 |
%\includegraphics{gkw91slope.eps} |
| 302 |
%\includegraphics{dm95slope.eps} |
| 303 |
%\includegraphics{ldd97slope.eps} |
| 304 |
\caption{Effective slope magnitude at 100~m depth evaluated using a) |
| 305 |
Cox slope clipping, b) GKW91 limiting, c) DM95 limiting and d) LDD97 |
| 306 |
limiting.} |
| 307 |
\end{figure} |
| 308 |
|
| 309 |
\section{Discretisation and code} |
| 310 |
|
| 311 |
The Gent-McWilliams-Redi parameterization is implemented through the |
| 312 |
package ``gmredi''. There are two necessary calls to ``gmredi'' |
| 313 |
routines other than initialization; 1) to calculate the slope tensor |
| 314 |
as a function of the current model state ({\bf gmredi\_calc\_tensor}) |
| 315 |
and 2) evaluation of the lateral and vertical fluxes due to gradients |
| 316 |
along isopycnals or bolus transport ({\bf gmredi\_xtransport}, {\bf |
| 317 |
gmredi\_ytransport} and {\bf gmredi-rtransport}). |
| 318 |
|
| 319 |
Each element of the tensor is discretised to be adiabatic and so that |
| 320 |
there would be no flux if the gmredi operator is applied to buoyancy. |
| 321 |
To acheive this we have to consider both these constraints for each |
| 322 |
row of the tensor, each row corresponding to a 'u', 'v' or 'w' point |
| 323 |
on the model grid. |
| 324 |
|
| 325 |
|
| 326 |
|
| 327 |
|
| 328 |
\vspace*{2cm} |
| 329 |
|
| 330 |
This is the old documentation..... |
| 331 |
|
| 332 |
|
| 333 |
The code that implements the Redi/GM/Griffies schemes involves an |
| 334 |
original core routine {\bf inc\_tracer()} that is used to calculate |
| 335 |
the tendency in the tracers (namely, salt and potential temperature) |
| 336 |
and a new routine {\bf RediTensor()} that calculates the tensor |
| 337 |
components and $\kappa_{GM}$. |
| 338 |
|
| 339 |
\subsection{subroutine RediTensor()} |
| 340 |
|
| 341 |
{\small |
| 342 |
\begin{verbatim} |
| 343 |
subroutine RediTensor(Temp,Salt,Kredigm,K31,K32,K33, nIter,DumpFlag) |
| 344 |
|---in--| |-------out-------| |
| 345 |
! Input |
| 346 |
real Temp(Nx,Ny,Nz) ! Potential temperature |
| 347 |
real Salt(Nx,Ny,Nz) ! Salinity |
| 348 |
! Output |
| 349 |
real Kredigm(Nx,Ny,Nz) ! Redi/GM eddy coefficient |
| 350 |
real K31(Nx,Ny,Nz) ! Redi/GM (3,1) tensor component |
| 351 |
real K32(Nx,Ny,Nz) ! Redi/GM (3,2) tensor component |
| 352 |
real K33(Nx,Ny,Nz) ! Redi/GM (3,3) tensor component |
| 353 |
! Auxiliary input |
| 354 |
integer nIter ! interation/time-step number |
| 355 |
logical DumpFlag ! flag to indicate routine should ``dump'' |
| 356 |
\end{verbatim} |
| 357 |
} |
| 358 |
|
| 359 |
The subroutine {\bf RediTensor()} is called from {\bf model()} with |
| 360 |
input arguments $T$ and $S$. It returns the 3D-arrays {\tt Kredigm}, |
| 361 |
{\t K31}, {\tt K32} and {\tt K33} which represent $\kappa_{GM}$ (at |
| 362 |
$T/S$ points) and the three components of the bottom row in the |
| 363 |
Redi/GM tensor; $2 S_x$, $2 S_y$ and $|S|^2$ respectively, all at $W$ |
| 364 |
points. |
| 365 |
|
| 366 |
The discretisations and algorithm within {\bf RediTensor()} are as |
| 367 |
follows. The routine first calculates the locally reference potential |
| 368 |
density $\sigma_\theta$ from $T$ and $S$ and calculates the potential |
| 369 |
density gradients in subroutine {\bf gradSigma()}: |
| 370 |
|
| 371 |
\centerline{\begin{tabular}{ccl} |
| 372 |
& & \\ |
| 373 |
Array & Grid-point & Definition \\ |
| 374 |
{\tt SigX} & U & |
| 375 |
$\sigma_x = \frac{1}{\Delta x} \delta_x \sigma|_{z(k)}$ |
| 376 |
\\ |
| 377 |
{\tt SigY} & V & |
| 378 |
$\sigma_y = \frac{1}{\Delta y} \delta_y \sigma|_{z(k)}$ |
| 379 |
\\ |
| 380 |
{\tt SigZ} & W & |
| 381 |
$\sigma_z = \frac{1}{\Delta z} |
| 382 |
[ \sigma|_{z(k)}(k-1/2) - \sigma|_{z(k)}(k+1/2) ]$ |
| 383 |
\\ |
| 384 |
\\ |
| 385 |
\end{tabular}} |
| 386 |
|
| 387 |
Note that $\sigma_z$ is the static stability because the potential |
| 388 |
densities are referenced to the same reference level ($W$-level). |
| 389 |
|
| 390 |
The next step calculates the three tensor components {\tt K13}, {\tt |
| 391 |
K23} and {\tt K33} in subroutine {\bf KtensorWface()}. First, the |
| 392 |
lateral gradients $\sigma_x$ and $\sigma_y$ are interpolated to the |
| 393 |
$W$ points and stored in intermediate variables: |
| 394 |
\begin{eqnarray*} |
| 395 |
\mbox{\tt Sx} & = & \overline{ \overline{ \sigma_x }^x }^z \\ |
| 396 |
\mbox{\tt Sy} & = & \overline{ \overline{ \sigma_y }^y }^z |
| 397 |
\end{eqnarray*} |
| 398 |
Next, the magnitude of ${\bf \nabla}_z \sigma$ is stored in an intermediate |
| 399 |
variable: |
| 400 |
\begin{displaymath} |
| 401 |
\mbox{\tt Sxy2} = \sqrt{ {\tt Sx}^2 + {\tt Sy}^2 } |
| 402 |
\end{displaymath} |
| 403 |
The stratification ($\sigma_z$) is ``checked'' such that the slope |
| 404 |
vector has magnitude less than or equal to {\tt Smax} and stored in |
| 405 |
an intermediate variable: |
| 406 |
\begin{displaymath} |
| 407 |
\mbox{\tt Sz} = \max ( \sigma_z , - \mbox{\tt Sxy2/Smax} ) |
| 408 |
\end{displaymath} |
| 409 |
This guarantees stability and at the same time retains the lateral |
| 410 |
orientation of the slope vector. The tensor components are then calculated: |
| 411 |
\begin{eqnarray*} |
| 412 |
\mbox{\tt K13} & = & -2 {\tt Sx/Sz} \\ |
| 413 |
\mbox{\tt K23} & = & -2 {\tt Sx/Sz} \\ |
| 414 |
\mbox{\tt K33} & = & ({\tt Sx/Sz})^2 + ({\tt Sy/Sz})^2 |
| 415 |
\end{eqnarray*} |
| 416 |
|
| 417 |
Finally, {\tt Kredigm} ($\kappa_{GM}$) is calculated in subroutine |
| 418 |
{\bf GMRediCoefficient()}. First, all the gradients are interpolated |
| 419 |
to the $T/S$ points and stored in intermediate variables: |
| 420 |
\begin{eqnarray*} |
| 421 |
\mbox{\tt Sx} & = & \overline{ \sigma_x }^x \\ |
| 422 |
\mbox{\tt Sy} & = & \overline{ \sigma_y }^y \\ |
| 423 |
\mbox{\tt Sz} & = & \overline{ \sigma_z }^z |
| 424 |
\end{eqnarray*} |
| 425 |
Again, a nominal stratification is found by ``check'' the magnitude of |
| 426 |
the slope vector but here is converted to a Brunt-Vasala frequency: |
| 427 |
\begin{eqnarray*} |
| 428 |
{\tt M2} & = & \sqrt{ {\tt Sx}^2 + {\tt Sy}^2} \\ |
| 429 |
{\tt N2} & = & - \frac{g}{\rho_o} \max ( {\tt Sz} , -{\tt M2 / Smax} |
| 430 |
\end{eqnarray*} |
| 431 |
The magnitude of the slope is then $|S| = {\tt M2}/{\tt N2}$. The Eady |
| 432 |
growth rate is defined as $|f|/\sqrt(Ri) = |S| N$ and is calculated |
| 433 |
as: |
| 434 |
\begin{displaymath} |
| 435 |
{\tt FrRi} = \frac{\tt M2}{\tt N2} ( - \frac{g}{\rho} {\tt Sz} ) |
| 436 |
\end{displaymath} |
| 437 |
The Eady growth rate is then averaged over the upper layers (about |
| 438 |
1100m) and $\kappa_{GM}$ specified from this 2D-variable: |
| 439 |
\begin{displaymath} |
| 440 |
{\tt Kredigm} = 0.02 * (200d3 **2) * {\tt FrRi} |
| 441 |
\end{displaymath} |
| 442 |
|
| 443 |
\subsection{subroutine inc\_tracer()} |
| 444 |
|
| 445 |
{\bf inc\-tracer()} is called from {\bf model()} and has {\em four |
| 446 |
new} arguments: |
| 447 |
\begin{verbatim} |
| 448 |
subroutine inc_tracer( ...,Kredigm,K31,K32,K33, ... ) |
| 449 |
real Kredigm(Nx,Ny,Nz) ! Eddy coefficient |
| 450 |
real K31(Nx,Ny,Nz) ! (3,1) tensor coefficient |
| 451 |
real K32(Nx,Ny,Nz) ! (3,2) tensor coefficient |
| 452 |
real K33(Nx,Ny,Nz) ! (3,3) tensor coefficient |
| 453 |
\end{verbatim} |
| 454 |
|
| 455 |
Within the routine, the lateral fluxes, {\tt fluxWest} and {\tt |
| 456 |
fluxSouth}, in the Redi/GM/Griffies scheme are very similar to the |
| 457 |
conventional horizontal diffusion terms except that the diffusion |
| 458 |
coefficient is a function of space and must be interpolated from the |
| 459 |
$T/S$ points: |
| 460 |
\begin{eqnarray*} |
| 461 |
{\tt fluxWest}(\tau) & = & \ldots + |
| 462 |
\overline{\tt Kredigm}^x \partial_x \tau \\ |
| 463 |
{\tt fluxSouth}(\tau) & = & \ldots + |
| 464 |
\overline{\tt Kredigm}^y \partial_y \tau |
| 465 |
\end{eqnarray*} |
| 466 |
|
| 467 |
The Redi/GM/Griffies scheme adds three terms to the vertical flux |
| 468 |
({\tt fluxUpper}) in the tracer equation. It is discretise simply: |
| 469 |
\begin{displaymath} |
| 470 |
{\tt fluxUpper}(\tau) = \ldots + \overline{\tt Kredigm}^z |
| 471 |
\left( |
| 472 |
{\tt K13} \overline{\partial_x \tau}^{xz} + |
| 473 |
{\tt K23} \overline{\partial_y \tau}^{yz} + |
| 474 |
{\tt K33} \partial_z \tau |
| 475 |
\right) |
| 476 |
\end{displaymath} |
| 477 |
On boundaries, {\tt fluxUpper} is set to zero. |
| 478 |
|
| 479 |
\section{Outstanding issues?} |
| 480 |
|
| 481 |
A second innovation discussed by Griffies concerns the nonlinearity in |
| 482 |
the equation of state and its implication for the evaluation of the |
| 483 |
isoneutral slopes. He argues that ${\bf \nabla} \sigma$ should be |
| 484 |
written ${\bf \nabla} \sigma = \rho_o ( \beta {\bf \nabla} S - \alpha |
| 485 |
{\bf \nabla} \theta )$ so that, for example, the zonal slope becomes: |
| 486 |
\begin{displaymath} |
| 487 |
S_x = \frac{ \alpha \partial_x \theta - \beta \partial_x S }{ \beta |
| 488 |
\partial_z S - \alpha \partial_z \theta } |
| 489 |
\end{displaymath} |
| 490 |
where $\alpha=\frac{1}{\rho_o} \frac{\partial \rho}{\partial \theta}$ |
| 491 |
and $\beta=\frac{1}{\rho_o} \frac{\partial \rho}{\partial S}$ are |
| 492 |
locally evaluated. He also advocates a very complicated interpolation |
| 493 |
scheme which he calls the ``triads'' scheme. Both of these innovations |
| 494 |
are designed to overcome computational modes and at the same time make |
| 495 |
the scheme {\em more} adiabatic. |
| 496 |
|
| 497 |
Compared to the current implementation out-lined above, the Griffies |
| 498 |
innovations requires two non-linear equations to replace the equation |
| 499 |
of state, twice as many differences because of the explicit |
| 500 |
expressions with $\theta$ and $S$, and a much wider stencil (ie. much |
| 501 |
more computation) for the ``triads'' interpolation. |
| 502 |
|
| 503 |
We have not implemented these aspects of Griffies schemes on the basis |
| 504 |
that the additional computational costs and code complexity, so far, |
| 505 |
outweigh the immediate benefits. |
| 506 |
|
| 507 |
\end{document} |
| 508 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph