F. Parameterization__of__Mesoscale__Eddies__________________________
In addition to the diffusion of tracers along isopycnals (Redi, 1982; Cox, 1987), the
parameterization of Gent & McWilliams (1990) includes an additional tracer advection.
The new equation for potential temperature, , is
t + L+ () = R(AI ; ) + (Asio Vz )z ; (F1)
with
L+ (G) = ___1_____acosOE[(u + u* )G] + [(v + v* )GcosOE]OE + [(w + w* )G]z : (F2)
Here, u, v, and w are the longitudinal (), latitudinal (OE), and vertical (z) Eulerian-mean
velocity components, respectively, t is time, a is the mean radius of the earth. Also, R is
the diffusion operator along isopycnals (see below; see also Redi, 1982; Cox, 1987; Gent &
McWilliams, 1990), AI is the isopycnal tracer diffusion coefficient, and u* , v* , and w* are
the longitudinal, latitudinal, and vertical components of the non-divergent, eddy-induced
transport velocity, respectively. This transport velocity is defined by
u* = AIT_D____acosOE%_%; v* = AIT_D____%OE_ ;
z z a %z z (F3)
w* = - ____1____acosOE AIT_D____acosOE%_%+ AIT_D____%OEcosOE___ ;
z a %z OE
where % is local potential density, and AIT D is the isopycnal thickness diffusivity. In
(F1), Ao V is the diapycnal diffusion coefficient, and the superscript si indicates that the
small background values are greatly enhanced in regions of static instability. Finally, G
represents a generic variable in (F2).
Note that in (F2) advection is by the effective transport velocity which is the sum of
the Eulerian-mean and eddy-induced transport velocities. The same conservation equation,
(F1), is used also for S. The discretization of the isopycnal transport parameterization
formulae is now given.
Consistent with the GFDL ocean general circulation model (Bryan, 1969; Cox, 1984),
all of the continuous equations of the isopycnal transport parameterization are discretized
on the Arakawa B-grid. In addition, our implementation of the parameterization is a
variant of that of Cox (1987) and Pacanowski et al. (1991, 1993). In this section, we use
the generic shorthand notation for spatial differencing, i.e.,
ffi = __1___[Gi+1=2 - Gi-1=2 ]: (F4)
i
Here, the subscript i denotes the discrete spatial location at which the gradient is evaluated
and 1=2 refers to one-half grid point offset with respect to the point i. , OE, and z
F-1
are the grid spacings in the longitudinal, latitudinal, and vertical directions, respectively.
Similarly, we apply the following generic shorthand notation for spatial averaging,
___ 1
G = __2[Gi+1=2 + Gi-1=2 ]: (F5)
Because the leap-frog scheme used in the model is unconditionally unstable if the
diffusion terms are evaluated at the central time level, the isopycnal diffusion and eddy-
induced transport terms are evaluated at the previous time level. In the following discrete
relations, we omit the temporal index for clarity.
A. Isopycnal Diffusion
The diffusion operator R, given in (F1), is defined by
R(AI ; G) r3D . [AI K . r3D G]; (F6)
where r3D is the three-dimensional gradient operator in spherical coordinates and K is a
three-dimensional second-rank tensor,
T
K SI SS. S : (F7)
Here, I is the 2 x 2 identity matrix, T denotes transpose, and
S -r%=%z : (F8)
In the derivation of (F7), the small slope approximation is used, i. e., |S| <<1.
Local potential densities, %, are computed relative to some selected reference depths,
differing from Cox (1987). This approach expedites the parameterization implementation
in a general coordinate model. In the current 2 x 2 and 3 x 3 models, we use 6 and 8
reference depth (pressure) levels. For example, 150, 650, 1500, 2500, 3500, and 4500m
are used as reference levels in the 3 x 3 model. Thus, % is referenced to a depth at most
500m away. The test cases with more reference depths showed no significant changes in
the solutions.
Whenever a new set of local potential density gradients is computed, the slopes of
the isopycnal surfaces are checked to prevent numerical instability (Cox, 1987). For this
purpose, we consider two approaches. The first method, already implemented in the GFDL
model, is based on changing the value of %z if the slope is larger than a permissible
maximum slope, Smax , e. g., q ____________
%2 + %2OE
check = - _______________S; (F9)
max
if %z > check =) %z = check: (F10)
F-2
Because the slopes of the isopycnals are limited by Smax , slopes larger than this produce
spurious diapycnal diffusion. This can be avoided to some extent by increasing Smax at the
expense of smaller time steps. In our present solutions, we use Smax = 0:01, allowing the
inclusion of high slopes encountered in frontal regions. This first approach is the default
when the ifdef option isopycmix is defined.
As stated above, a major disadvantage of the maximum slope approach is that it will
produce diapycnal diffusion for slopes larger than Smax . To remedy this problem, the
second method, while still allowing mixing along the computed slopes, actually tapers AI
and AIT D to zero. To activate this method, the ifdef option isopycmixspatialvar must
be defined in addition to the isopycmix option. The tapering is accomplished by the
application of two functions,
AI =) f1 f2 AI ; AIT D =) f1 f2 AIT D : (F11)
Here, f1 is a hyperbolic tangent function,
f1 = 0:5 1 + tanh -|S|_+_Sc______S ; (F12)
d
which diminishes the diffusion coefficients as |S| gets large. In (F12), Sc is the slope at
which f1 = 0:5, and Sd is the half length of the interval in which f1 changes with a steep
slope. (F12) is applied in the interval 0 -|S| -Smax . Outside of this interval, f1 is set
to zero. The second function, f2 , is designed to turn off the isopycnal mixing scheme near
the ocean surface, thus allowing the specified vertical mixing processes to function fully.
It is based on the depth of the water parcel, the local Rossby deformation radius, and the
local potential density slope. The details of this function can be found in the isopyc.F
subroutine. We strongly recommend the use of the second approach especially for models
with high vertical resolution near the ocean surface. Note that in both (F9) and (F12), the
negative sign is used to denote that the isopycnal transport parameterization is applicable
only in stably stratified regions.
The (1-3) component of the tensor, K13 , is obtained at the center of the eastern face
of the tracer grid box. Consequently, the associated local potential density gradients are
obtained there using
T
% = m_____affi [%]; (F13)
%OE = 1_affiOE[__%;OE]; (F14)
%z = ffiz [__%;z ]: (F15)
In (F13), mT is secOE evaluated at the tracer grid points.
After checking the slopes of the isopycnal surfaces, K13 is computed from
K13 = - %___%: (F16)
z
F-3
Note that (F14) is exclusively used during the slope check procedure, and it is set to zero
if any one of the neighboring tracer grid points is land.
The (2-3) component of the tensor, K23 , is obtained at the center of the northern face
of the tracer grid box. Hence, the required potential density gradients are computed at
these locations. The discrete forms are
u __
% = m____affi [% ;OE ]; (F17)
%OE = 1_affiOE[%]; (F18)
%z = ffiz [__%OE;z]; (F19)
where mu is secOE computed at the velocity grid points. After slope check, K23 is computed
from %
K23 = - _OE_%: (F20)
z
(F17) is used only for slope check, and % = 0 if any one of the neighboring tracer grid
points is land.
K13 and K23 are subject to the zero flux condition at horizontal boundaries. In order
to compute the vertical gradients of % at the ocean surface and bottom, an extrapolative
estimator is used to obtain %.
K31 , K32 , and K33 are computed at the center of the vertical face of the tracer grid
box. The discrete forms for the local potential density gradients are
T ______________u;ffiz [%]
% = m_____a________________T; (F21)
______________uOE;z
%OE = 1_aOE__ffiOE[%]_____T; (F22)
OE
%z = ffiz [%]: (F23)
Here, the superscripts u and T refer to the velocity and tracer grid points, respectively.
Note that in case of equal mesh spacing in the horizontal directions, the use of u ffi or
OEu ffiOE is equivalent to computing averages of the gradients. After slope check, we use
K31 = - %___%; (F24)
z
K32 = - %OE_%; (F25)
z
%2 + %2OE
K33 = ____________%2: (F26)
z
For K31 , K32 , and K33 components of the tensor, in addition to the horizontal boundaries,
we impose the zero flux condition at the top and bottom boundaries.
F-4
Finally, (F6) is discretized as follows,
T mT ___;z
R(AI ; G) = m_____affi AI ______affi [G] + AI K13 ffiz [G ]
T 1 1 ___OE;z
+ m_____affiOE AI _______mufafiOE[G] + _____muAI K23 ffiz [G ] (F27)
T _______________u;z _______________uOE;z
+ ffiz m_____aAI K31 __ffi__[G]_________T+ 1_ AI K32 OE__ffiOE[G]_______T+ AI K33 ffiz [G] :
a OE
Equation (F27) is subject to the no tracer flux condition at all boundaries. The necessary
vertical gradients at the first level are computed assuming that the surface tracer values
are the same as those at the first level. For the bottom level, an extrapolative estimator
is used to compute the tracer values at the ocean floor.
B. Eddy-Induced Transport Velocity and Advection
The eddy-induced transport velocity components, u* , v* , and w* are computed at the
centers of the eastern, northern, and vertical faces of the tracer grid box, respectively. The
discrete forms for u* and v* are
______________z
u* = -ffiz [AIT D K13 ]; (F28)
______________z
v* = -ffiz [AIT D K23 ]: (F29)
The vertical component of the eddy-induced transport velocity is directly obtained by
integration of the continuity equation for this velocity. Note that the application of (F28)
and (F29) at the top and bottom levels requires the values of the tensor elements at the
ocean top and bottom. The requirement that the vertical eddy-induced velocity is zero
at the ocean top and bottom shows that AIT D K13 and AIT D K23 must be zero at these
surfaces.
Finally, the eddy-induced advection terms in (F2) are discretized as
T ___ v* ___GOE ___z
L(G) = m_____affi (u* G ) + ffiOE ________mu + ffiz (w* G ): (F 30)
REMARKS:
1. In order to eliminate any horizontal diffusion of tracers, the horizontal diffusivity
coefficient, ah, must be set to zero.
2. AI and AIT D have the same spatial dependency.
3. Because the current implementation uses an extrapolative formula near the ocean
surface and bottom, the values in the bottom topography array, kmt, must be 2 in
the active ocean.
F-5
4. If the ifdef option isopycmix is specified, the ifdef option skipland must not be used.
5. Although it has the same option name as the Redi (1982) and Cox (1987) isopycnal
mixing parameterization in MOM 1.0 and MOM 1.1, it is NOT the Redi-Cox version.
6. AI and AIT D may have vertical variation (fzisop in the model). The default is no
depth variation. To utilize this capability, the user must insert a function into the 2586
do-loop in ocean.F. The vertical variation is ineffective when the isopycmixspatialvar
option is defined.
Bryan, K., 1969: A numerical method for the study of the circulation of the world ocean.
J. Computational Phys., 4, 347-376.
Cox, M. D., 1984: A primitive equation, 3-dimensional model of the ocean. GFDL Ocean
Group Tech. Rep. No. 1, GFDL/Princeton University, Princeton.
Cox, M. D., 1987: Isopycnal diffusion in a z-coordinate ocean model. Ocean Modelling,
74, 1-5.
Gent, P. R. & J. C. McWilliams, 1990: Isopycnal mixing in ocean circulation models. J.
Phys. Oceanogr., 20, 150-155.
Pacanowski, R. C., K. Dixon, & A. Rosati, 1991, 1993: The GFDL modular ocean model
users guide. GFDL Ocean Group Tech. Rep. No. 2, GFDL, Princeton.
Redi, M. H., 1982: Oceanic isopycnal mixing by coordinate rotation. J. Phys. Oceanogr.,
12, 1154-1158.
F-6