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