H. Numerical_and_Coupling_Improvements______ a. Third-Order Upwind Algorithm Third-order upwind schemes and variants of this scheme are discussed fully * *in Leonard (1979). A somewhat more complicated scheme has been put into the MOM code by the FRAM group, see Farrow and Stevens (1995). The third-order upwind scheme is swi* *tched on by setting the ifdef option upwind3. The finite-difference expression for the advection of tracers is given by ADV (i; j; k)-=___1____dxt(__uETE*- __uWTW*) - _1__(__vNTN*- __vSTS*) i cstj dytj - _1__dzt(wk-1 TT*- wk TB*) k where __uW; __uE; __vS; __vNare the T-grid cell face centered velocity componen* *ts. __u ui;jdyuj_+_ui;j-1dyuj-1_ __ __ E(i)= 2dytj uW (i) = uE(i - 1) __v (vi;jdxui+_vi-1;jdxui-1)csuj_ N(j)= 2dxticstj __v (vi;j-1dxui+_vi-1;j-1dxui-1)csuj-1_ S(j)= 2dxticstj and the T *quantities are the cell face centered tracer concentration at the ea* *st, west, north, south, top, and bottom faces. In the standard second order scheme of MOM* *, these are taken as the arithmetic average of the tracer concentrations in adjacent ce* *lls, e.g., TE*= 1_2(Ti+1;j+ Ti;j) In the third-order upwind formulation, the T *s are replaced by an upwind d* *irection biased three-point interpolation to the interface. In the zonal direction they * *are: ( __ flTi-1 + fi+ Ti+ ff+ Ti+1 uE > 0 TE*(i) = __ fi- Ti+ ff- Ti+1+ ffiTi+2 u E < 0 TW*(i) = TE*(i - 1) H-1 where the coefficients ff; fi; fl; ffi are given by ff+= _________dxti_(2dxti+_dxti-1)________(dxt i+ dxti+1) (dxt* *i+1+ dxti-1+ 2dxti) ff-= _____dxti(2dxti+1+_dxti+2)____(dxt i+ dxti+1) (dxti+1+ dx* *ti+2) fi+= ___dxti+1_(2dxti+_dxti-1)___(dxt i+ dxti+1) (dxti+ dxti-1) fi-= _______dxti+1(2dxti+1+_dxti+2)_______(dxt i+ dxti+1) (dxt* *i+ dxti+2+ 2dxti+1) fl= ____________-dxti*_dxti+1____________(dxt i+ dxti-1) (dxt* *i+1+ dxti-1+ 2dxti) ffi= _____________-dxti*_dxti+1_____________(dxt i+1+ dxti+2* *) (dxti+ dxti+2+ 2dxti+1) The interface tracer values TN*; TS*; TT*; and TB*are obtained in an analogous * *manner. Farrow, D.E., and D.P. Stevens, 1995: A new tracer advection scheme for Bryan a* *nd Cox type ocean general circulation models. J. Phys. Oceanogr., 25, 1731-1741. Leonard, B.P., 1979: A stable and accurate convective modeling procedure based* * on quadratic upstream interpolation. Comp. Meth. in Appl. Eng., 19, 59-98. b. Single Vertical Velocity The second improvement concerns the vertical velocity computation. In the o* *riginal scheme, two independent vertical velocity computations, one in the momentum and* * one in the tracer equations, are performed, resulting in two different vertical vel* *ocity fields, with the former typically much noisier than the latter. In the new method, the* * vertical velocity is computed only for the tracer equations. Consequently, this velocit* *y is then area-averaged to obtain the vertical velocity for the momentum equations, elimi* *nating the second computation. The related advection operator in the momentum equations is* * also modified accordingly to maintain the flux and kinetic energy conservation prope* *rties. This new treatment of the vertical velocity produces a unique and smooth vertical ve* *locity field with rather modest effects on global distributions. Further details can be foun* *d in Webb (1995), where the new scheme is described in Section 7. Webb, D.J., 1995: The vertical advection of momentum in Bryan-Cox-Semtner ocean general circulation models. J. Phys. Oceanogr., 25, 3186-3195. H-2 c. North Pole Implementation This section describes a consistent way to treat the north pole in B-grid m* *odels. The traditional finite difference approach breaks down because the zonal grid spaci* *ng goes to zero and there is no set of values to the north of the north pole. The approach* * described here treats the north pole as a circular grid point of radius dytjmt which carr* *ies values for all tracers and the transport streamfunction at its center. All grid points to * *the south of the pole are arranged in a standard B-grid, with internal velocity points offse* *t one half grid spacing to the north and east from tracer and transport streamfunction points. * *There is no velocity point which corresponds to the tracer point at the pole. This option i* *s activated by specifying the ifdef option northpole. The momentum equations may be solved at the row just to the south of the no* *rth pole (jmt-1) using the standard solution procedure outlined in Section C. It is assu* *med here that a Laplacian frictional parameterization is used, some considerations may be nec* *essary in order not to overflow arrays when using biharmonic or higher frictional paramet* *erizations. The value of the internal mode velocities at j=jmt is not defined, but this pos* *es no problem because the area of the interface between the jmt-1 and jmt velocity grid boxes* * is zero. All terms in the finite difference formulation which multiply ui;jmtor vi;jmtare al* *so multiplied by cstjmt (the cosine at the latitude of the the j=jmt T-point) which, for the * *north pole, is zero. The solution of the tracer equation at the north pole requires some special* * treatment. The finite difference form of the tracer evolution equation at the circular gri* *d may be written (neglecting friction) as himt-1X T n+1= T n-1+ 2t (vi;jmt-1dxui+ vi-1;jmt-1dxui-1) i=2 i (Ti;jmt+ Ti;jmt-1)csujmt-1=(dytjmtdxti) : This is simply the sum of all meridional fluxes into the circular grid cell fro* *m the regular grid to the south. All terms on the right hand side are known. Once solved, t* *he single new tracer value at the north pole can be filled into the j=jmt row of the trac* *er array so that the meridional fluxes through the interfaces with the interior grid can be* * calculated during the next time step. The transport streamfunction is located at the same points as the tracers, * *so this will also require some special treatment at the pole. We take a similar approach by * *vertically averaging the internal mode zonal momentum equation and computing its line inte* *gral along the velocity points at row j=jmt-1. The contribution due to the surface * *pressure gradient identically vanishes, leaving Z Z __u __0 td = u td: H-3 __ where u0tis the vertical average of the internal mode tendency and is longitud* *e. Using the definition of the transport streamfunction, this relation may be rewritten * *as Z Z __ - _1_aH_d = u0td; where a is the mean radius of the earth and H is the model depth at tracer poin* *ts. Hence, the barotropic circulation tendency around the pole determines the v* *alue of the streamfunction tendency ( _) at the pole, relative to the average interi* *or value. This equation may then be used to solve for the time rate of change of the tran* *sport streamfunction at the north pole (p_ole) as imt-1X 2dxt ____________i_________pole= i=2 Hi;jmt-1+ Hi-1;jmt-1 imt-1X 2dxt imt-1X ____________i_________i;jmt-1- dxtidyujmt-1__u0 :t i=2 Hi;jmt-1+ Hi-1;jmt-1 i=2 The forcing term on the right hand side is known and an estimate of the str* *eamfunction tendency at j=jmt-1 is obtained during the standard iterative solution for the * *streamfunction in the interior. Once the estimate for the tendency at the pole is found, the * *usual iterative solution procedure is resumed. As for the tracer arrays, the single * *values for the streamfunction tendency and the streamfunction, once solved, must be stored* * in the appropriate array so the interior solution can proceed. REMARK: The present implementation of the northpole ifdef option assumes equal * *grid intervals in the zonal direction when computing fluxes into the north pole grid* *box. d. Tapering of Horizontal and Isopycnal Mixing Coefficients Our previous experience with global models showed that the diffusive comput* *ational stability at high latitudes, rather than advective stability, is usually the li* *miting factor for the model time step. Consequently, in the present model, we allow the hori* *zontal eddy viscosity coefficients AMH , AI, and AITD to vary in the meridional dire* *ction. In particular, they are reduced at high latitudes. This treatment necessitates no * *changes in the scalar tracer equations. In contrast, the vector momentum equation requires* * additional horizontal friction terms so that pure solid body rotation does not produce any* * stress on the fluid. Consequently, we modified these viscous terms to follow Wajsowicz (1993)* *. Because the model uses the leap-frog time integration scheme, the viscous/diffusive com* *putational stability limits in the zonal direction in which the grid spacing becomes small* *er at high latitudes are given by AMH ___t_U_____2 1_; AI____t_T____2 1_: (x cosOE) 4 (x cosOE) 4 H-4 In the present discussion, we consider only this direction for simplicity, beca* *use it is the most restrictive. In the above equations, t U and t T are the maximum momentum * *and tracer time steps, respectively, x is the longitudinal grid spacing, and OE is * *latitude. The model determines the allowable AMH , AI, and AITD from these stability equatio* *ns, so that all these coefficients are kept at their full values as poleward as possible. I* *n addition, to prevent excessively small AMH , AI, and AITD which can produce instabilities d* *ue to lack of enough diffusion, we restrict the minimum values to be 20% and 6% of their s* *pecified values in the 3x3 and 2x2 models, respectively. In the 3x3 model, AMH tapering* * starts at 81.9ON, and it is kept constant at the 20% level starting at 87.3ON. AI and * *AITD satisfy the above constraint at all latitudes. In the 2x2 model, AMH , and AI and AITD* * taperings start at 84.6ON and 85.2ON, respectively. They are kept constant at the 6% leve* *l north of 88.8ON. No tapering is applied in the Southern Hemisphere in both model configu* *rations. This option is activated by setting the ifdef option taperhrzvd. Wajsowicz, R. C., 1993: A consistent formulation of the anisotropic stress tens* *or for use in models of the large scale ocean circulation. J. Comput. Phys., 105, 333-* *338. e. K. Bryan Acceleration Technique The model time integration is performed in two stages. First, the accelerat* *ion technique of Bryan (1984) is applied in which the momentum and tracer equations have vari* *able time steps, i. e., n+1 n-1 ff ^u___-_^u___2t= Fn-1;n; n+1 n-1 fl(z)[;_S]____-_[;_S]__2t= Gn-1;n: The temporal discretization is done using the leap-frog scheme, and n - 1, n, a* *nd n + 1 represent the previous, current, and next time levels, respectively. Also, ^uis* * the horizontal baroclinic velocity vector, t is the fundamental (surface tracer) time step, i* *s potential temperature, and S is salinity. F and G represent the rest of the terms in thes* *e equations, involving quantities from both the n - 1 and n time levels. ff and fl, which ca* *n vary in the vertical (z), are the time scale, or acceleration, factors. In the case of* * ff > fl, the momentum time step is smaller than the tracer time step by a factor of fl=ff. The model calender and all the forcing fields are based on the surface trac* *er time step (fl-1 t) by setting fl = 1 above 1058 m. The depth and timesteps given he* *re and below apply to the 3x3 model version. With this choice, the annual cycles of he* *at and salt forcing are consistently applied. Here, t = 440 min was chosen, well within the* * advective computational stability limit. The smaller current speeds at depth are taken ad* *vantage of by setting fl = 0:2 between 1058 and 2502 m and fl = 0:1 between 2502 and 5000 * *m depth. Consequently, the deep ocean tracer time step is about 3 days; ten times larger* * than at the surface. To illustrate a disadvantage of this depth acceleration, consider a v* *ertical heat flux, Q, between two adjacent model layers, k and k + 1. Integration of the equ* *ation above H-5 leads to a change in the total heat content of the two layers by an amount H = 2t Q (fl-1k - fl-1k+1) : Conservation of heat requires H = 0, which is violated whenever flk 6= flk+1. T* *herefore, vertical discontinuities in fl are restricted to depths where the fluxes are th* *ought to be small enough that tracers are well conserved. In the momentum equations, horizontal diffusion determines computational st* *ability. At high latitudes this factor limits the baroclinic velocity time step to about* * 44 min, ff = 10, at all depths. To be consistent a time step of ff-1t is also used in t* *he barotropic streamfunction equation. With this time step, the surface wind forcing is not c* *onsistently applied, because the annual cycle harmonic is applied 10 times per momentum yea* *r. As a result, some distortion of the seasonal cycle of ocean currents is expected. * * Analogous problems are not serious in the tracer fields, because their depth acceleration* * (fl 6= 1) is restricted to depths where the influence of the seasonal cycle on tracers is ne* *gligible. The acceleration technique is implemented by setting the values of dtts, dtuv, dtsf* *, and dtxcel in the model input data. In the second stage of the time integration, the accelerated results were u* *sed as initial conditions for a synchronous time integration in which equal time steps and no * *depth acceleration were applied, i. e., t = 44 min and ff = fl = 1. This integration* * needs to be long enough to remedy any erroneous behavior spawned by violating conservati* *on laws and by forcing mismatches in the acceleration stage, see Danabasoglu et al. (19* *96). Bryan, K., 1984: Accelerating the convergence to equilibrium of ocean-climate * *models. J. Phys. Oceanogr., 14, 666-673. Danabasoglu, G., J.C. McWilliams, and W.G. Large, 1996: Approach to equilibrium* * in accelerated global oceanic models. J. Climate, in press. f. Ocean Surface Slope The components of the ocean surface slope are computed exactly the same way* * as the diagnostic surface pressure gradient in the original MOM code. These slope comp* *onents are accumulated at each timestep and only time-averaged fields over the nysnc s* *teps in a coupling interval are passed to the flux coupler. The surface slope is nondimen* *sionalized from the surface pressure gradient by dividing by the gravitational acceleratio* *n. This calculation is activated by setting the ifdef option sflxpvmicetilt. g. Sea-Ice Formation and Melting This improvement is needed when the ocean model is coupled to a sea-ice mod* *el, and is activated by the ifdef option sflxiceflx. At the end of each timestep, n, a* *fter surface fluxes and model transports have been applied, the upper kmxice model layers, w* *here kmxice is a model input parameter, are scanned for temperatures below the freez* *ing point, H-6 Tf. At present Tf is set to the constant value of -1:8OC. However, the freezing* * point is a function of salinity and a more accurate equation for Tf, that will be implemen* *ted in the near future, is Tf = - :0575 SK + 0:001710523 S3=2K - 0:0002154996 S2K ; at the local salinity, SK . At each layer the potential mass per unit area of i* *ce formation (potice > 0) or ice melt (potice < 0) is computed as ! potice (k) = MAX (%Cp)o_Hk (Tf - Tk) ; qice (k - 1) ; f where (%Cp)o = 4.1 x 106 j=m3=K, is the product of ocean density and heat capac* *ity, and f = 334000 j=Kg is the latent heat of fusion. Any ice that forms at depth is as* *sumed to float towards the surface and this ice flux, qice(k) (defined positive downw* *ards), so qice(k) 0 is accumulated bottom to top as Xk qice (k) = -potice (k) ; kmxice where qice(kmxice-1) = 0, assumes no ice formation below the kmxice layer. As i* *ce floats to the surface, it can either partially or completely melt in upper layers whos* *e temperatures are above freezing. At each layer the temperature and salinity are both adjusted in accord with* * the ice formed or melted in the layer: Tk = Tk + ___f_____H potice(k) k(%Cp)o Sk = Sk + (So_-_Si)_%potice(k) : o Hk The ice flux is accumulated at each location over the number of timesteps i* *n a coupling interval, nsync, as nsyncX aqice(i; j) = qice(k = 1)n: n=1 If, at the end of a vertical scan, the upper layer temperature T1 remains g* *reater than freezing, (qice(k = 1) = 0), then the excess heat melts previously formed * *ice and aqice(i; j), T1 and S1 are adjusted accordingly. Thus over the nsync timesteps* *, ice is assumed to remain where it was formed. After nsync timesteps, the accumulated i* *ce flux aqice(i; j) can be passed through the Flux Coupler to the sea-ice model as an e* *quivalent downward heat flux (W=m2). qflux(i; j) = __-f_____nsyncatqice(i; j) where aqice includes any adjustment due to melting previously formed ice in the* * last timestep, and t is the model timestep. H-7 Thus, ocean temperature and salinity adjustments corresponding to the total* * ice formed (qflux > 0) are made prior to the ice being passed to the sea-ice model.* * However, warm surface model temperatures result in (qflux < 0), which is a potential to * *melt ice in the sea-ice model. Since the ocean does not know if sufficient ice is presen* *t at a given location, temperature and salinity adjustments are delayed until the appropriat* *e heat and freshwater fluxes are received back from the sea-ice model through the Flux Cou* *pler. Note, if sufficient ice is present in the absence of all other fluxes, the heat flux * *will be equal to the cooling needed to make T1 = Tf, when applied over the next nsync timesteps. H-8