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