ATMOCEAN.460777 1998/01/01 A GLOBAL ATMOSPHERE-OCEAN MODEL FOR CLIMATE STUDIES: MODEL FORMULATION Gary L. Russell NASA / Goddard Space Flight Center Institute for Space Studies 2880 Broadway, New York, NY 10025 Abstract A new coupled atmosphere-ocean model has been developed for climate predictions at decade to century time scales. Atmospheric source terms are similar to those of Hansen et al. [1983]. Atmospheric dynamic equations for mass and momentum are solved using Arakawa and Lamb's [1977] C grid scheme. Advection of potential enthalpy and water vapor uses the linear upstream scheme [Russell and Lerner, 1981]. The new global ocean model conserves mass as opposed to volume and allows for divergent flow and a free surface. Ocean mass and momentum equations are solved using a new C grid scheme. The linear upstream scheme is used for advection of potential enthalpy and salt. A version of the k-profile parameterization (KPP) mixing sceme [Large et al., 1994] is used for ocean vertical mixing. The atmosphere-ocean model uses realistic land distribution and realistic solid ground topography. Insolation has proper seasonal and diurnal cycles. Atmospheric and ocean surface fluxes are of opposite sign and are applied synchronously. Flux corrections are not used. Except for partial strength alternating binomial filters which are applied to the momentum components in the atmosphere and oceans, there is no explicit horizontal diffusion. Climate simulations of the coupled model at 4 x 5 horizontal resolution with 9 vertical layers for the atmosphere and 13 layers for the ocean is discussed by Russell et al. [2000]. Table of Contents Abstract Table of Contents 1. Introduction 2. The grid 3. Prognostic variables 4. The time scheme 5. Atmospheric mass and momentum equations 6. Atmospheric polar filter 7. Ocean mass and momentum equations excluding the pressure gradient force 8. Ocean pressure gradient force 9. Ocean polar filter 10. Alternating binomial filter on the momentum components 11. Advection of tracers 12. Moist convection and large scale condensation 13. Radiation 14. Surface fluxes 15. River flow 16. Ocean vertical mixing 17. Ocean bottom friction 18. Ocean straits 19. Sea ice formation, melting and thermal diffusion 20. Sea ice advection 21. Problems and future improvements 22. Conclusion 23. Acknowledgements Appendix: Ocean functions and the equation of state References 1. Introduction Global atmospheric general circulation models (AGCMs) have been used extensively to simulate the Earth's present climate with reasonable success. AGCMs also have been used to estimate climate changes that might occur by adding greenhouse gases to the atmosphere. Because the poleward heat transport by the ocean circulation is comparable to that by the atmospheric circulation, climate change simulations with AGCMs need some technique to incorporate ocean heat transport. Two such techniques and their combination are discussed. In technique 1, Russell et al. [1985] ran the Goddard Institute for Space Studies AGCM with climatologically specified ocean surface temperatures and sea ice distributions, and for each grid box derived the daily converged ocean heat transport that is necessary to maintain the current climate. Hansen et al. [1984, 1988] incorporated these converged ocean heat transports with the same AGCM to estimate the climate response to doubling of atmospheric carbon dioxide and the climate response to three greenhouse gas scenarios in transient simulations. They were able to estimate three principal feedback factors that drive the global temperature change: the water vapor feedback factor, the snow-ice albedo feedback factor, and the cloud feedback factor. The use of technique 1 prevented them from estimating the ocean circulation feedback factor, which may be of the same order. In technique 2, Washington and Meehl [1989] combined a version of the NCAR atmospheric Community Climate Model with the ocean model developed by Bryan [1969] and modified by Semtner [1974, 1987]. The ocean circulation of this coupled model responds to climate change. The present model, a synchronously coupled atmosphere-ocean model but with radically different dynamics and physics, also incorporates the ocean circulation feedback factor. The transport of heat by the oceans in the present model is accomplished by dynamical advection and not by diffusion. Both of these atmosphere-ocean models develop a climate drift in the deep ocean and significant differences from climatological ocean surface temperatures at a few grid boxes. Manabe et al. [1991] and Sausen et al. [1988] essentially combined techniques 1 and 2. They separately built coupled atmosphere-ocean models as in technique 2 but added surface flux correction terms. Technique 1 essentially uses flux correction terms with an ocean model with no circulation. Manabe et al. and Sausen et al. show annual flux correction terms of comparable magnitude to the total flux used by Hansen et al. [1984]. In theory, technique 2 is superior to technique 1 because it allows the ocean circulation to respond to climate changes during a simulation. In practice, the existence of climate drift or the necessity of using flux correction terms indicate model deficiencies. A principal task for current model development is to improve on these atmosphere-ocean models. The coupled model described in this paper does not solve these problems, but it does introduce a few improvements and new insights on such models. 2. The grid The global atmosphere-ocean model described in this paper is a grid point model that simulates the planet Earth. Grid boxes are defined by a regular latitude-longitude grid in the horizontal and by layers in the vertical. The vertical coordinate, m (kg/m), increases downward and at any level is the mass per unit area above the level; this includes the masses of air, ice and water. Although the local pressure P (Pa) equals mg, m is used as the vertical coordinate instead of P because mass is a conserved quantity and g (m/s), the Earth's effective gravitational acceleration, has variations in both space and time although not so in the model at present. (Do not confuse the m coordinate with the abbreviation for meter which is often used in the parenthetical units that follow a quantity.) In the atmosphere, each column uses the same number of vertical layers. There is a constant level, MT (kg/m), above which atmospheric dynamics do not occur. The three dimensional prognostic variable MA (kg/m) is the mass thickness per unit area of a layer. The vertical sum of MA is the atmospheric column mass from MT to the surface. The model can accept layers of constant MA's in the upper part of the atmosphere and a sigma coordinate system in the lower part. The standard resolution of 4 in latitude and 5 in longitude uses only a 9 layer sigma coordinate system whose constant mass fractions counted from the surface upward are: .051335, .082135, .137577, .174538, .164271, .138604, .107803, .082135 and .061602. In the ocean, the number of grid layers depends on horizontal location but is constant in time. At the end of the dynamical subroutines, the ratio of the sea water masses between two layers in a column depends only on the two layer numbers and not on horizontal location or time. Each ocean column for the standard resolution uses 13 or fewer layers, counted downward, and the ratio of an ocean grid box mass divided by the grid box mass below it is 2/3 at the end of the dynamics. Atmospheric and ocean topography is realistic and is measured from mean sea level. The ocean surface height varies in space and time. The realistic bottom topography under oceans is adjusted for the standard resolution to be -24(1.5-1) meters where Z is the number of layers in the ocean column. Because of the vertical column mass ratios and adjusted bottom topography, the ocean layer thicknesses for the standard resolution are approximately 12, 18, 27, 40, 61, 91, 137, 205, 308, 461, 692, 1038 and 1557 meters. Each horizontal grid square is either all land or all ocean. Each land square has fractions that are ground, glacial ice, open lake, or lake ice. Each ocean square may contain sea ice, whose horizontal fraction and single thickness varies with location and time according to conservation of mass and energy. Figure 1 shows the number of grid boxes for each ocean column for the standard resolution. The Central American isthmus and the Malaysian peninsula and the islands south of it have been filled in and a northwest passage between North America and Greenland has been opened. Large islands have been retained. 3. Prognostic variables The atmospheric prognostic variables are three dimensional arrays and are listed in Table 1. Evaporation or precipitation of water does not affect the atmosphere's mass, but it does affect the ground reservoirs or the ocean's mass. The only effect of water vapor mass on the atmospheric model is to measure latent energy. This inconsistency is retained in order to treat the atmosphere as a single substance ideal gas which allows potential enthalpy to be advected as a conserved tracer in the atmosphere. The definition of atmospheric potential specific enthalpy, G (J/kg), is G = CT(P/P) (3.1) where C (J/kgK) is the specific heat capacity of dry air at constant pressure, T (K) is temperature, P is the constant mean sea level pressure, P is pressure, x = R/C is the exponent of the exner function, and R (J/kgK) is the gas constant for dry air. The ocean prognostic variables are three dimensional arrays and are listed in Table 2. Ocean entropy, which is conserved under adiabatic advection, would be the ideal prognostic variable for heat were advection properties the only consideration. But the ocean equation of state, which is used to calculate the specific volume from the prognostic variables, is not known accurately for entropy. For that reason potential enthalpy, for which the equation of state can be calculated accurately, is chosen as the prognostic variable for heat. The determination of the equation of state from the prognostic variables is shown in the Appendix. Potential specific enthalpy, G, of a parcel at depth in the ocean is defined to be the specific enthalpy of the parcel at standard atmospheric pressure, P, if the parcel were raised adiabatically to P. Advecting G as a conserved tracer is not strictly correct; for example two parcels at the same depth, when mixed, should conserve their enthalpy. In the ocean model, those parcels are raised adiabatically to P, are mixed conserving enthalpy at P, and then lowered adiabatically to the original depth. In spite of this inaccuracy, G is a better prognostic variable than would be potential temperature because it properly accounts for variations in salinity and specific heat capacity. Mass, potential enthalpy, water vapor, and salt are defined as integrals over the three dimensional grid boxes defined in section 2. For the latter three quantities, the gradients with respect to mass are defined in each direction. The horizontal velocity components are defined on the C grid as shown in Figure 2. Thus UA and UO are defined at the edge between two adjacent west-east mass boxes and VA and VO are defined between two adjacent south-north mass boxes. Ocean velocity components between a land and an ocean box or between two land boxes are set to zero. The prognostic variables for the land reservoirs and for sea ice and lake ice reservoirs are shown in Table 3. For the standard resolution, these reservoirs contain two vertical layers for mass and heat. Mass of snow contains its own prognostic variable, but the heat content of snow is lumped with the heat of the upper layer. Lake temperature and lake ice mass are not prognostic variables, but are specified from climatology. Lake and river mass above the sill depth is prognostic. 4. The time scheme The prognostic variables are given at the start of a simulation, and are updated iteratively by the following processes in order: Moist convection and large scale condensation of water vapor and calculation of clouds Update of ground and ice reservoirs due to precipitation, and calculation of runoff Radiation which uses the calculated clouds Calculation of momentum stress, evaporation, and sensible and latent heat fluxes, and their application to atmospheric variables Update of ground and ice reservoirs due to evaporation and surface flux heating, and melting of sea ice River flow Update of surface ocean variables due to precipitation, radiation, momentum stress, evaporation, and surface fluxes, and freezing of ocean water Ocean convection and vertical diffusion Ocean bottom friction Acceleration of strait flow by pressure gradient force Convection and vertical diffusion in straits Bottom friction in straits Advection of mass, potential enthalpy, and salt through straits Atmospheric dynamics of mass and momentum, and advection of potential enthalpy and water vapor Alternating binomial filter on atmospheric momentum components Ocean dynamics of mass and momentum Ocean advection of potential enthalpy and salt Alternating binomial filter on ocean momentum components The atmospheric dynamics of mass, momentum and tracers (potential enthalpy and water vapor) contains its own sub-loop with a leap frog time scheme. After the prognostic variables have been updated by the source terms, forward and backward steps are performed on the dynamic terms to advance the odd solution one dynamic time step which is half of the leap frog time step. The leap frog scheme is applied to the dynamic terms of the solutions, ending with a final even solution. The tracers are advected only for the even solutions. The calculation of the pressure gradient force for the even leap frog steps of momentum uses a time average of potential enthalpy at the beginning and end of the even step. For the ocean model, potential specific enthalpy and salinity are kept constant while a leap frog time scheme is applied to the dynamic terms of ocean mass and momentum. After mass and momentum have been integrated, the summed mass fluxes from the even steps are used to advect ocean potential enthalpy and salt. For the 4 by 5 standard resolution, the iterative loop at the beginning of this section is executed once for each simulated hour. The dynamic time step for both the atmosphere and ocean is 7.5 minutes. The atmosphere-ocean model integrates 4 simulated years each real day on an SGI Origin 3000 single processor. The computer time percentages for the individual components are listed in Table 4. 5. Atmospheric mass and momentum equations After discretizing the vertical coordinate, the form of the primitive equations on the spherical grid used in the atmospheric model are: M 1 Mu Mvcos + + + Wê - Wë = 0, (5.1) t acos Mu 1 Muu Mvucos + + + (Wu)ê - (Wu)ë t acos tan M (+CT) M(P/P) G - f + uMv + - = 0, (5.2) a acos acos Mv 1 Muv Mvvcos + + + (Wv)ê - (Wv)ë t acos tan M (+CT) M(P/P) G + f + uMu + - = 0, (5.3) a a a where M is mass per unit area, u is zonal velocity, v is meridional velocity, t is time, is longitude, is latitude, a is the radius of the Earth, f is the Coriolis parameter 2׎sin, is geopotential, and G is defined by (3.1) below which C, T, P, P and R are defined. W (kg/ms) is the upward vertical mass flux which is located at discrete layer edges such as Z+. W depends on the horizontal mass convergence and is defined so that M's in different layers maintain their specified mass in the upper atmosphere or their specified sigma fractions in the lower atmosphere. Arakawa and Lamb's [1977] second order C grid scheme is used to integrate the discretized mass equation and the nonlinear terms of the the discretized momentum equation. In the finite difference equations, M becomes MA, u becomes UA, and v becomes VA. Arakawa and Lamb's [1977] choice for the Coriolis and metric terms is not used because it does not conserve zonal angular momentum even if the corner fluxes in the nonlinear terms are eliminated. The new form of the Coriolis and metric terms in (5.2) would exactly conserve zonal angular momentum were the corner fluxes in the nonlinear terms eliminated. The Coriolis and metric terms for the meridional momentum equation are designed to conserve kinetic energy in concert with the new zonal momentum equation. This improvement was tested in shallow water equations models on a spherical grid. For the Rossby-Haurwitz wave 6 initial conditions, the model with the old form of the Coriolis and metric terms produced significant eddy kinetic energy near the poles which was not present in the new model nor in higher resolution simulations of either model. It is possible to design forms for the Coriolis and metric terms that would conserve zonal angular momentum even with the corner fluxes included in the nonlinear terms, but it is inadvisable because the meridional momentum equation would have a term with v in the denominator in order to conserve kinetic energy. Because the atmosphere is considered to be dry and R and C are constants, many different forms for the pressure gradient force can be used. For example: dP + d = d(+CT) + (P/P)͎dG (5.4) where is the specific volume of dry air. The form for the pressure gradient force chosen here has two terms: a spatial derivative of geopotential plus specific enthalpy and the exner function multiplied by a spatial derivative of potential specific enthalpy. The first term represents the gradient of specific energy and the second term represents the deviation of the pressure gradient force due to the vertical coordinate. The second term would vanish if potential specific enthalpy were the vertical coordinate and each layer contained the same value of potential specific enthalpy for all horizontal columns. These choices for the vertical coordinate and pressure gradient force are ideal for Arakawa's potential enstrophy conserving schemes [Arakawa and Lamb, 1981; Takano and Wurtele, 1982] because they would again conserve potential enstrophy as they do in the shallow water equations. 6. Atmospheric polar filter The atmospheric polar filter used in this model is similar to that described by Arakawa and Lamb [1977] and is identical to that described by Hansen et al. [1983]. The purpose of the filter is to allow a longer dynamical time step (7.5 minutes for 4 x 5 horizontal resolution) than would be required to satisfy the Courant-Friedrichs- Lewy condition. Both the zonal mass flux and the zonal pressure gradient force are spectrally analyzed in the east-west direction. Spectral coefficients of each field are multiplied by a smoothing factor S which decreases with wave number and with proximity to the poles according to the formula: x S = (6.1) y sin(n/2) where x is the east-west grid distance for the particular latitude, y is the constant north-south grid distance, n is the particular wave number, and is 2 divided by the number of east-west grid boxes. If S is greater than 1, it is reset to 1 before multiplying the spectral coefficients. The modified spectral coefficients are then transformed back into grid point values. 7. Ocean mass and momentum equation excluding the pressure gradient force After discretizing the vertical coordinate, the form of the primitive equations on the spherical grid used in the ocean model are: M 1 Mu Mvcos + + + Wê - Wë = 0, (7.1) t acos Mu 1 Muu Mvucos + + + (Wu)ê - (Wu)ë t acos tan M P M - f + uMv + + = 0, (7.2) a acos acos Mv 1 Muv Mvvcos + + + (Wv)ê - (Wv)ë t acos tan M P M + f + uMu + + = 0. (7.3) a a a Most variables have the same meanings as those used in section 5 except (m/kg) is the specific volume of sea water. For the ocean, W (kg/ms) is the downward vertical mass flux which is located at discrete layer edges such as Z+. W depends on the horizontal mass convergence and is defined so that M's in different layers tend to their specified ratios, namely that M of a layer is 1.5 times that of the layer above for the standard vertical resolution. In the finite difference equations, M becomes MO, u becomes UO, and v becomes VO. In an earlier version of this ocean model, the nonlinear advection terms for the momentum equation were taken from Arakawa and Lamb [1977]. In that version, the mass associated with a velocity point came from six mass boxes, some of which may have been land boxes having no ocean mass. This seemed to cause serious problems only in the Arctic Ocean. We note that Arakawa and Lamb's scheme conserved enstrophy under nondivergent flow and energy away from the coasts, but it did not conserve zonal angular momentum by advection. In the current dynamical scheme, only the two mass boxes on either side of a velocity point are used to calculate the momentum. The north-south flux of zonal momentum flows into the corner boxes as well as to the boxes directly north or south of a momentum box in a way similar to that in Arakawa and Lamb's [1977] scheme. However, the east-west flux of zonal momentum flows only into the boxes directly east or west of a momentum box. For meridional momentum, all orientations are switched. This new scheme conserves energy away from the coasts but not enstrophy. The new Coriolis and metric terms for zonal momentum are designed to conserve zonal angular momentum exactly by advection away from the coasts. The Coriolis and metric terms for meridional momentum are designed not to generate kinetic energy in concert with the zonal momentum. Conservation of zonal angular momentum near the poles improves the flow there. 8. Ocean pressure gradient force The pressure gradient force accelerates the velocity between two mass grid boxes at the same layer. The force in a layer has two terms: the gradient of the mean pressure multiplied by the average mean volume of the two boxes and the gradient of the mean geopotential multiplied by the average mass of the two boxes. The horizontal gradients of potential enthalpy and salt are not used in calculating the above mean quantities, but the vertical gradients are used. The mean vertical coordinate, m, of a grid box is the average of the top edge m and the bottom edge m. m is the difference of m at the two edges. The mean pressure of a grid box is P(m) = mg, where g (m/s) is the uniform effective gravity. To calculate the mean specific volume of a grid box we assume that (G(m),S(m),P(m)) = (m) fits a quadratic polynomial in m. The mean value of the quadratic (m) from m-m to m+m is equal to [(m-x) + (m+x)] where x = m/12 (see Figure 3). The mean volume of a grid box is equal to the mass multiplied by the mean specific volume which is calculated using quadratic precision using two evaluations of (G,S,P) per grid box at the vertical coordinates m-x and m+x. This shows how the first term of the pressure gradient term is modeled. The second term is modeled as follows. We again assume that (m) fits a quadratic in each grid box. The height, h (meters), at the column bottom is specified by the bottom topography. Integrating upwards, the change in h from the bottom layer edge to any level in a layer is equal to the integral of (m)dm. The change in h to the top layer edge uses the two evaluations of mentioned in the previous paragraph and is calculated with quadratic precision. The top edge value of h of a layer becomes the bottom edge value of h for the layer above. The mean geopotential, (m/s), of the layer is the mass weighted throughout the layer which is a double integral of (m) and is equal to _ 1 m+m = (m) dm m m-m 1 m+m m+m = (m+m) + g (n) dn dm m m-m m = (m+m) + [(m-x)(3-3) + (m+x)(3+3)] mg/12 (8.1) where is equal to gh. The mean is calculated with quadratic precision using the same two evaluations of (G,S,P). Given an arbitrary quadratic polynomial in an interval [-m,m], the abcissas of the intersection points between the polynomial and the least square fit line that fits the polynomial in the interval are the points -x and x mentioned above. These points do not depend on the coefficients of the polynomial. Conversely, there are no two other points such that the mean value of an arbitrary quadratic polynomial can be derived from evaluating the polynomial at the two points. This mathematical result is even more elegant because the double integral in (8.1) can be calculated from an evaluation of the polynomial at the points -x and x. 9. Ocean polar filter For 4 by 5 resolution, the 7 minute time step is chosen as the largest convenient time step such that gravity waves of wave length 2y, or 8 at the equator, do not cause the numerical solution of the momentum equation to diverge. The short gravity waves that are resolved at middle and high latitudes in the zonal direction presumably could cause the numerical solution to diverge. To prevent this, the zonal pressure gradient force and the zonal velocity used in the zonal mass flux are spectrally analyzed, and coefficients of waves shorter than 2y are multiplied by the wave length divided by 2y. This technique is used in the model of Arakawa and Lamb [1977] and other grid point atmospheric models. Ocean models are different because separate ocean basins should be filtered and are filtered separately. For an individual ocean with N grid boxes from coast to coast at a given latitude, the two quantities to be filtered are defined on the N-1 edges. The quantities are spectrally analyzed on 2N periodic values which include zeroes on each coast and N-1 reflected values of opposite sign. Again, the spectral coefficients of waves shorter than 2y are multiplied by the wave length divided by 2y. Since each filtered value is a linear combination of the N-1 unfiltered values, the coefficients of the linear combinations are calculated once offline. These coefficients depend only on N and latitude. 10. Alternating binomial filter on the momentum components C grid schemes have a strong tendency to develop horizontal grid alternating patterns in their velocity fields in regions of mass variations, i.e. near coast lines and mountains. Alternating patterns are even more noticeable in the vertical mass flux field which is derived from the horizontal mass convergence. The alternating pattern is an artifact of the numerical scheme and not of the physical flow. It is efficiently suppressed by a high order alternating binomial filter which selectively damps the alternating pattern of a field. The filter also damps other short waves but not as strongly. In the atmospheric model, a one dimensional eighth order alternating binomial filter is applied with full strength to the zonal and meridional momentum componets every hour in the east-west direction. For the zonal component: __ __ 1 8 Y+8 __ MAUA(I) MAUA(I) - (-1) MAUA(I+Y) (10.1) 4 -8 16 where MA is the average of the mass boxes associated with the velocity point, UA is the prognostic velocity, I is an arbitrary longitude index, Y is the summation index, and (Y+8 16) is the binomial coefficient "Y+8 choose 16". In the ocean model, a one dimensional eighth order alternating binomial filter is modified to properly handle coast lines and the poles. The filter is applied to the UO and VO fields once every four hours at one quarter strength in the east-west direction and once every sixth hour at one sixth strength in the north-south direction. (A quarter strength filter means that the final changes in the fields are multiplied by one quarter.) The reduction of ocean kinetic energy by the alternating binomial filter is approximately one fifth of the production by the atmospheric stress. 11. Advection of tracers A significant difference between this atmosphere-ocean model and other climate models is the use of the linear upstream scheme [Russell and Lerner, 1981] for the advection of potential enthalpy, water vapor, and salt. For each of these tracers, each grid box predicts the mean tracer and the tracer gradients in each of the three directions with respect to mass. The tracer distribution defined by the mean and the gradients have discontinuities at the grid box edges. Advection is performed one direction at a time. Figure 4 shows how the linear upstream scheme performs advection in one direction. If the gradients were set to zero at the beginning of each step, the linear upstream scheme would reduce to the standard upstream scheme. Russell and Lerner [1981] showed that the linear upstream scheme is more accurate than a fourth order tracer advection scheme when the mass of the fluid is the same in all grid boxes. When there are mass variations, the linear upstream scheme is decidedly superior to the second order or fourth order schemes. Mass variations occur due to topography variations in sigma coordinate general circulation models, due to the convergence of meridians toward the pole in latitude- longitude grids, and due to nonuniform vertical layering as in this atmosphere-ocean model. Tracer concentration is defined as the tracer mass divided by the fluid mass. The natural coordinates for a tracer are the mass coordinates of the fluid with respect to which the tracer is defined. In a three dimensional model, the mass coordinates are functions of the three spatial coordinates and time whereas a tracer is a function only of the three mass coordinates. When mass variations exist, the tracer is defined on an irregular grid of the mass coordinates; and the second order, fourth order, and spectral tracer advection schemes become inaccurate because they expect the tracer to be defined on a regular grid. This is demonstrated in Russell and Lerner [1981]. Mass variations are not a problem for upstream schemes because they make no inherent assumptions that a tracer is defined on a regular grid. The use of center difference or spectral schemes for tracer advection causes a slowly increasing systematic error that may be unnoticed in a tropospheric model where convection overwhelms and corrects the temperature and humidity fields, but the unstableness is deadly to a stratospheric or ocean model. The normal remedy is the application of strong diffusion. The linear upstream scheme on the other hand is stable, produces smooth fields, and is only weakly diffusive. 12. Moist convection and large scale condensation Using the mean and gradients, potential enthalpy and water vapor are calculated for each horizontal quarter cell for moist convection and for each eighth cell for large scale condensation. Moist convection is performed each hour, but large scale condensation is applied only once every fifth hour, just prior to calling the radiation subroutines. The Model does not have the ability to store condensate in the atmosphere; condensate is reevaporated or precipitated each time step. Cloud optical depths, used by radiation, are determined from the amount of condensate leaving each layer and the local temperature. Only one of the four quarter columns is used to determine cloud optical depths that radiation applies to whole grid cells. For moist convection, a plume is formed with one fourth of the mass of a quarter cell of layer 1, the lowest atmospheric layer. The plume temporarily rises to the layer above it where water vapor in excess of saturation is condensed. If there is no condensation, nor if the plume, with the weight of the condensate, is less buoyant than that of the environment air, then the plume is abandoned. If the plume is more buoyant, then the condensate is removed, and the plume temporarily rises to the next higher layer, continuing the process. Subsidence of the environment occurs using the vertical gradients of heat and water vapor. Condensate may reevaporate as it falls. After this process is completed, another plume is formed with one fourth of the mass of a quarter cell of layer 2, and the process is repeated. Initial plume layers continue up to the second highest layer in the troposphere. The whole process is repeated for all horizontal quarter columns and for all grid columns. Note that plumes do not entrain adjacent air from their environment as they are rising, and that plumes are slightly under saturated after condensation. For large scale condensation and starting from the top, if water vapor of an eighth of a cell exceeds saturation, then vapor is condensed until the concentration is slightly below saturation. The condensate falls into the eighth cell below it, where it may be fully or partially reevaporated. Super saturation conditions are checked for each eighth cell. 13. Radiation Atmospheric and surface radiation is similar to that described by Hansen et al. [1983], except for the calculation of optical depths due to clouds. Each wave length for insolation is attributed to one of 6 discontinuous spectral domains dependent on the absorption line's strength. Thermal radiation uses 25 spectral domains which are again apportioned due to absorption line strength. Radiation calculations are performed for each spectral domain after which the vertical radiative fluxes are summed. All significant gases, cloud particles and background aerosols are used. Solar albedoes and thermal atmospheric heating rates are calculated on whole grid cells and are maintained for five hours. Incident sun light at the top of the atmosphere, which is calculated half hourly, is multiplied by the albedoes to determine solar radiation heating rates. The upward thermal radiation from the surface is recalculated each time the ground temperature changes. Cloud optical depths depend on temperature between 0 and -30C, on the square root of the condensate amount leaving a grid cell, on the square root of the air mass in the cell, and on whether the clouds are due to moist convection or large scale condensation. The formula for the cloud optical depth for each layer is: / [0.075 MAMCC + 0.30 MALSC] when T > 0C OD = < (13.1) \ [0.025 MAMCC + 0.01 MALSC] when T < -30C where MA is the air mass defined in section 3, MCC (kg/m) is the moist convective condensate leaving the layer, LSC (kg/m) is the large scale condensate leaving the layer, and T is the layer temperature. When the temperature is between 0 and -30C, the quantities inside the square roots are linearized with respect to temperature. Because radiative fluxes are not proportional to optical depths, the optical depths from each horizontal quarter column are not averaged but are cyclically used every fourth call to the radiation routines. The Model's cloud diagnostics do not affect radiation but are determined from the cloud optical depths. If the vertically integrated cloud optical depth is less than 1 at a given time, then the cloud diagnostics accumulate nothing. If the depth is greater than 1, then the cloud diagnostics are incremented in the level where optical depth above it is 1. Low clouds are diagnosed if that level is in the first three layers (on average below 720 mb); high clouds are diagnosed for layers 6 through 9 (on average above 390 mb). Although the Model calculates cloud optical depths simultaneously at all levels, Model diagnostics report only the highest clouds at each time step in order to compare with satellite observations. Absorbed solar radiation on the open ocean surface penetrates through the first three layers of the ocean according to the formula of Paulson and Simpson [1977]. The absorbed insolation affects both the mean and the vertical gradients of potential enthalpy in the ocean. 14. Surface fluxes The atmosphere-ocean model calculates the following fluxes at the surface: Fresh water flux from evaporation and precipitation Vertical heat flux including solar radiation, thermal radiation, latent heating, sensible heating, and precipitation heat flux Momentum surface stress As mentioned in section 2, each land grid square may contain horizontal fractions of ground, glacial ice, open lake, and lake ice whereas each ocean grid square may contain horizontal fractions of open ocean and sea ice. The surface fluxes are calculated separately for each surface fraction and are applied to the atmosphere and the appropriate subsurface reservoir. The ocean receives these surface fluxes without corrections every hour for the standard resolution. Horizontal gradients of potential enthalpy and water vapor are ignored during the calculations. Prior to applying the surface fluxes, the surface air speed, WS, is assumed to be 3/4 of the first layer air speed, and the surface air temperature and humidity are 1/3 of the way down the vertical gradient from the first layer center to the ground. The Richardson number, Ri, between the ground and surface values is calculated from which the drag coefficient, CD, and cross isobar angle, , are determined: CD = CDN / (0.5 + 2Ri), (14.1) 0.09375 - 0.03125 / (1 + 4Ri) = 2 , (14.2) 1 + 0.3WS where CDN is the constant neutral drag coefficient that depends on topography and vegetation. If Ri is negative, it is reset to 0 before CD and are evaluated. Classical drag laws are used to calculate momentum, heat and moisture transfers using the same CD. In Hansen et al. [1983], the initial surface fluxes and their derivatives with respect to ground temperature were evaluated after which an implicit time scheme was used to calculate the final surface fluxes. In the present model, derivatives of the fluxes with respect to surface air temperature and humidity are also evaluated, and a multi-variable implicit time scheme is used to determine the final surface fluxes. This simplified surface parameterization with a better implicit time scheme is a significant improvement over the complex coding in Hansen et al. The ocean-sea ice heat flux is identical to that in Hansen et al. [1983]. The first ocean layer receives the atmospheric surface fluxes on the open ocean and the ocean-sea ice fluxes under sea ice. Because there is no ocean mixed layer model, the vertical gradients of potential enthalpy and salt are set to zero for the ocean's first layer. 15. River flow The atmosphere-ocean model uses river flow to conserve water and to return continental runoff to the oceans. The river direction file and an offline version of the river routing model are described by Miller et al. [1994]. For each continental grid box, the river direction file indicates which of the eight adjacent or diagonal boxes is downstream from the given box, except that some continental boxes have no outlet. The river mass flux leaving a grid box to its downstream neighbor is F (kg/s) = MRu/d (15.1) where MR (kg) is the river and lake mass above the sill depth in the box and d (m) is the mean distance between a grid box and its downstream neighbor. u (m/s) is an effective flow speed that depends on the topography gradient, i: u = 0.35 i/i (15.2) where i = 0.00005 is the reference topography gradient. u is limited to vary between 0.15 and 5 m/s. If MR is negative, then F is set to zero. Source runoff or unabsorbed precipitation and snow melt in a grid box immediately increases the river and lake mass of the box. Precipitation and evaporation on the lake fraction of a continental grid box also increases or decreases the river and lake mass of the box. Substantial evaporation on the lake fraction can cause MR to become negative. River and lake mass above the sill depth is passed downstream to lower rivers and lakes and eventually into the ocean. Water mass and energy are conserved globally by river transport, although lake energy by itself is not conserved. The coefficients of u were optimized by Miller et al. [1994] so that the fresh water flux and energy flux from rivers are fed into the proper coastal ocean boxes with the proper seasonal timing. 16. Ocean vertical mixing Ocean vertical mixing is similar to that described by Large et al. [1994]. 17. Ocean bottom friction The ocean velocity of the lowest layer of each ocean column is reduced in order to simulate bottom friction. The bottom drag vector is (kg/ms) = CD (u,v) (u,v) (17.1) where CD is a dimensionless drag coefficient, (kg/m) is density, and (u,v) is velocity. In the ocean model, CD is set to 1 kg/m. The destruction of kinetic energy by the bottom drag is two orders of magnitude smaller than the production by the atmospheric stress. 18. Ocean straits Prognostic variables for straits, listed in Table 5, are arrays with one vertical dimension. Mass of water in a strait is nondivergent and is not prognostic. The mass flux is dynamically accelerated only by the pressure gradient force between the whole ocean grid boxes at each end. The calculation is similar to that between normal grid boxes except that the distance in the denominator is equal to the length of the strait plus the distances from the straits' ends to the respective grid box centers. Advection of potential enthalpy and salt through a strait follows the basic principles of the linear upstream scheme but is more complicated. The exchange of water between the end of a strait and its whole ocean grid box does not occur at a single edge of the grid box, but occurs in the interior of the grid box. Convection and vertical diffusion are applied to the two halves of each strait in the same way they are applied to the quarter boxes of whole ocean grid boxes. Bottom friction with the same coefficient is applied in the lowest layer in each strait. Sources and sinks are not allowed to occur in a strait. Parameters for the 12 straits for the standard resolution are listed in Table 6. 19. Sea ice formation, melting, and thermal diffusion The atmosphere-ocean model predicts sea ice by conserving mass and energy and by calculating their fluxes from above and below. When surface fluxes cause the first layer of the ocean to cool below the freezing point, the ocean stays at the freezing point and the required amount of sea ice is formed with a 0.5 m thickness. If sea ice draws sufficient heat from the ocean below to cool it to below freezing, then the ocean stays at the freezing point and the sea ice thickens. Sea ice contains no salt, so when sea ice forms or thickens, the salt concentration of the first ocean layer increases. When surface fluxes, principally insolation, cause the sea ice temperature to rise above 0C, the temperature stays at 0C and snow or sea ice is melted which runs into the ocean below. If sea ice becomes thinner than 0.5 m, then it is contracted horizontally so that it remains at 0.5 m thick. If the temperature of the first ocean layer rises above 0C, then sea ice is melted vertically and horizontally drawing the necessary energy from the ocean which cools back to 0C. 20. Advection of sea ice 21. Problems and future improvements This section describes a list of needed improvements to the atmosphere-ocean model which are currently being developed or which are considered important. The atmospheric model with climatological ocean surface temperatures and sea ice distributions is not in radiation balance. For the standard resolution, the atmospheric model adds annually and globally 7.4 W/m into the oceans. This may be the most significant cause of the global climate drift experienced by the coupled model. The principal reason for this imbalance may be that the model's ocean surface albedo is too low. Over the tropical ocean, the International Satellite Cloud Climatology Project data [Rossow and Zhang, 1994; Zhang et al., 1994] shows a surface albedo of about 8% whereas the model's albedo is less than 5%. As mentioned in section 11, center difference or spectral schemes for advection are inaccurate because grid boxes are unequally distributed in the mass coordinate system. This statement also applies to momentum and its related functions such as angular momentum or absolute vorticity times mass. Replacing the center difference advection scheme for momentum by a linear or quadratic upstream type scheme should significantly improve both the atmospheric and ocean models. This is not a straightforward task because of the pressure gradient force, the nonlinearity of momentum, and other factors, but progress is being made. Implementation of this change should allow reduction or elimination of the alternating binomial filters that are applied to the momentum components. The 0.8% of the Earth's surface covered by continental lakes use a climatologically specified lake surface temperature and lake ice fraction. Lakes currently do not respond to climate change. For completeness, lake parameters need to be predicted in the future. As mentioned in section 3, water vapor is not included in the atmospheric mass and does not affect the specific heat capacity of air. Implementation of this effect in an atmospheric model reduced the difference between the surface air temperature and the ocean surface temperature, which was then in much better agreement with observations. This coding will be incorporated into the atmosphere-ocean model. The model subroutines that calculate condensation and reevaporation of atmospheric water vapor are currently under scrutiny. Del Genio and Yao [1993] have improved the convective mass flux calculation and have added the capability of allowing convective downdrafts in some atmospheric model simulations. DelGenio et al. [1993] has also added prognostic variables for cloud condensate which may be retained in the atmosphere, reevaporated later, or precipitated to the ground. These improvements will be added to the atmosphere-ocean model. The Earth's gravitational acceleration is a global vertical constant in the present version of the atmosphere-ocean model, although the model was coded to allow it to vary. Tides will be added to the coupled model to test their effect and to compare with observations. Incorporation of tides may increase the ocean vertical mixing near coast lines and in the straits. At present, the radiation subroutines cannot accept vertical gradients of temperature, humidity or clouds that are calculated by other parts of the model, and the subroutines cannot return the vertical gradient of heating rates within each layer. Changing the radiation routines will not be done soon. As mentioned in the Appendix, (G,S,P) is interpolated from table values. The range of S in the table is from 0 to .04 (kg salt/kg sea water). This is because the range of S in the input functions, e.g. (T,S,P), is from 0 to .04. However, salinity often exceeds .04 in the Persian Gulf grid cell for the current climate, and may exceed that limit in other areas if the climate changes. The range of S in the table for (G,S,P) should be increased. The real ocean has numerous island chains and underwater mountain ridges which effectively prohibit the mass transport at some depth from one body of water to another, e.g., the Caribbean archipelago which prevents the excange of deep water between the Caribbean Sea and the Atlantic Ocean, and the Indonesian archipelago extending from Malaysia which restricts the flow of water between the Pacific and Indian Oceans. Such ridges may be narrow even though they may be very long. A finite resolution ocean model, which bases its mean bottom topography for each grid box on the area weighted topography within the box, may not detect the effect of such ridges. One possible solution is to raise the bottom topography at such locations which was done for the ocean model's standard resolution. Another solution (which was not implemented) is to reduce the mass flux between adjacent ocean boxes if there is a ridge between them. An analogous problem also occurs in the atmospheric model. 22. Conclusion This paper has described the mathematical formulations of a new atmosphere-ocean model. Many aspects of this model are unique or unusual including the linear upstream scheme for the advection of potential enthalpy, water vapor, and salt, a new C grid dynamical scheme for ocean momentum, continental river flow that feeds the ocean, the calculation of strait flow in the ocean, conservation of ocean mass, and prediction of the ocean surface height. The use of the linear upstream scheme for advection of the tracer quantities is the single most significant advantage of this coupled model over other climate models. The linear upstream scheme is stable and accurate and does not require external diffusion. Where other schemes would use a single prognostic variable for a quantity, the linear upsrtream scheme uses 4, the mean quantity and its three directional gradients. This implies that the effective resolution for tracer quantites is finer than the advertised resolution. The 4 x 5 standard resolution of this atmosphere-ocean model has been used for two simulations for the Intergovernmental Panel on Climate Change's 1995 Scientific Assessment of Climate Change [1996]. The first was a 97 year integration with fixed atmospheric composition, and the second was a 74 year integration with a compounded 1% per annum increase of CO. Diagnostics from the last ten years of each simulation have been delivered to the German Climate Computer Center in Hamburg which maintains a database of similar atmosphere- ocean model simulations from other groups. The difference in surface air temperature for the last ten years between the 1% CO experiment and the control was 1.42C. The seasonal simulation of the ocean for years 56 to 63 of the control simulation mentioned above have been analyzed by Miller and Russell [1994]. The key results of their investigation were: The global root mean square error between the coupled model's and climatological sea surface temperature is 2.4C in February and 2.6C in August. The model simulates the major features of the present climate system including the major ocean gyres. The global meridional ocean heat transport is 1.4 PW at 16N which is consistent with observations and other models. The Antarctic Circumpolar Current through the Drake Passage is 125 Sv, and the Gulf Stream transport at 32N is 44 Sv. The meridional mass stream function in the North Atlantic has a maximum of 11.7 Sv at 28N and 900 m. The development of a complex atmosphere-ocean model is a continuing process. The model described by this paper may not be in the middle of this process, but it is certainly not near either end. 23. Acknowledgements A complex atmosphere-ocean model cannot be built by one person. It requires contributions, large and small, from many individuals. James R. Miller of Rutgers University must be cited as a most significant contributor to both the model's development and the writing of this mauscript. Numerous other people, who, wittingly or not, contributed to the design of this atmosphere-ocean model, are listed in alphabetical order: Frank Abramopoulos, James Bishop, Kirk Bryan, Guilherme Caliri, Mark Cane, Mark Chandler, Anthony Del Genio, Keith Dixon, Inez Fung, Dale Haidvogel, James Hansen, Xingjian Jiang, Andrew Lacis, Jean Lerner, Jochem Marotzke, Douglas Martinson, Ernst Maier-Reimer, Paulo Malanotte-Rizzoli, David Rind, Anthony Rosati, Cynthia Rosenzweig, William Rossow, Reto Ruedy, Gavin Schmidt, Suki Sheth, Peter Stone, Robert Suozzo, and Hanspeter Zinn. Model development support was provided by the NASA Climate Program Office, by the Earth Observing System interdisciplinary funding, and by the United States Environmental Protection Agency office of Policy, Planning and Evaluation, Global Climate Change Division. Appendix: Ocean Functions and the Equation of State The following variables are defined for the ocean model: C (J/kgC) = specific heat capacity at constant pressure; G (J/kg) = potential specific enthalpy; g (m/s) = Earth's vertical gravitational acceleration; H (J/kg) = specific enthalpy; m (kg/m) = previously defined vertical coordinate; P (Pa) = pressure; P (Pa) = global mean atmospheric pressure at sea level; S (1) = salinity; T (C) = temperature; (m/kg) = specific volume; (m/kg) = potential specific volume; (C/Pa) = adiabatic lapse rate; (C) = potential temperature. G and S can be calculated at any location in the ocean from m and the prognostic variables. The purpose of this section is to derive P, H, T, and from G, S and m or G, S and P. Remember that P = mg. (T,S,P), C(T,S,P) and (T,S,P) are provided by Fofonoff and Millard [1983], but also see Fofonoff [1985]. H(T,S,P) is derived as follows: H/P = (T,S,P), (A.1) T/P = (T,S,P) (A.2) are integrated simultaneously from P to P at constant entropy and salinity to obtain H(T,S,P) - H(,S,P). P is the reference pressure for potential temperature so that = T at P. Next evaluate H(,S,P) - H(,0,P) at constant entropy and pressure which is provided by Millero and Leung [1976]. Finally, the differential equation H/T = C(T,0,P) (A.3) is integrated from to 0C at constant salinity and pressure to obtain H(,0,P) - H(0,0,P). The reference value for specific enthalpy is set by H(0,0,P) = 0. Thus H(T,S,P) is determined by the above formulas. (T,S,P) and T(,S,P) are derived by integrating (A.2) from P to P. The G of a parcel is defined to be the H of the parcel at P were it raised adiabatically. Thus, G(T,S,P) = H((T,S,P),S,P) = G(,S). (G,S) is derived by inverting the first argument of G(,S) which is possible because G and are in a one-to-one relationship at constant S and P. T, H, and can now be derived as functions of G, S and P. T(G,S,P) = T((G,S),S,P) = T(,S,P); H(G,S,P) = H(T((G,S),S,P),S,P); and (G,S,P) = (T((G,S),S,P),S,P). (G,S) is defined to be (G,S,P). The evaluation of (G,S,P), called the equation of state, is performed by linear interpolation within a look-up table in the ocean model's computer program. The increments in the table are 4000 (J/kg) for G, .001 for S, and 210 (Pa) for P. References Abramopoulos F, Rosenzweig C, and Choudhury B. 1988. Improved ground hydrology calculations for global climate models (GCMs): Soil water movement and evapotranspiration. J Clim 1: 921-941. Arakawa A and Lamb VR. 1977. Computational design of the basic dynamical processes of the UCLA general circulation model. Methods in Computational Physics, Vol. 17, Academic Press, 337 pp. Arakawa A and Lamb VR. 1981. A potential enstrophy and energy conserving scheme for the shallow water equations. Mon Wea Rev 109: 18-36. Bryan K. 1969. A numerical method for the study of the circulation of the world ocean. J Comput Phys 3: 347-376. DelGenio AD and Yao M-S. 1993. Efficient cumulus parameterization for long-term climate studies: the GISS scheme. In "The Representation of Cumulus Convection in numerical models", Emanuel KA and Raymond DA, eds. Amer Meteor Soc, Meteorological Monographs 24: 181-184. DelGenio AD, Yao M-S, and Wendell CE. 1993. GCM feedback sensitivity to interactive cloud water budget parameterization. In Preprints, Fourth Symposium on Global Change Studies, Anaheim, Amer Meteor Soc: 176-181. Fofonoff NP. 1985. Physical properties of seawater: A new salinity scale and the equation of state for seawater. J Geophys Res 90: 3332-3342. Fofonoff NP and Millard RC. 1983. Computation of fundamental properties of seawater. UNESCO, Technical Papers in Marine Science 44. Hansen J, Fung I, Lacis A, Lebedeff S, Rind D, Ruedy R, and Russell G. 1988. Global climate changes as forecast by the GISS 3-D model. J Geophys Res 92: 14739-14760. Hansen J, Lacis A, Rind D, Russell G, Stone P, Fung I, Ruedy R, and Lerner J. 1984. Climate sensitivity: Analysis of feedback mechanisms. In "Climate Processes and Climate Sensitivity", Hansen JE and Takahashi T, eds. Geophysical Monograph 29, AGU, Washington DC: 130-163. Hansen J, Russell G, Rind D, Stone P, Lacis A, Lebedeff S, Ruedy R, and Travis L. 1983. Efficient three-dimensional global models for climate studies: Models I and II. Mon Wea Rev 111: 609-662. Intergovernmental Panel on Climate Change. 1996. "Climate Change 1995", Houghton JT, Meira Filho LG, Chandler BA, Harris N, Kattenberg A, and Maskell K, eds. Cambridge University Press, 572 pp. Large WG, McWilliams JC, and Doney S.C. 1994. Oceaninc vertical mixing: Review and a model with non-local boundary parameterization. Rev Geophys 32: 363-403. Manabe S, Stouffer RJ, Spelman MJ, and Bryan K. 1991. Transient responses of a coupled ocean-atmosphere model to gradual changes of atmospheric CO; Part I: Annual mean response. J Clim 4: 785-818. Miller JR, Russell GR, and Caliri G. 1994. Continental scale river flow in climate models. J Clim 7: 914-926. Millero FJ and Leung WH. 1976. The thermodynamics of seawater at one atmosphere. Amer J Sci 276. Paulson CA and Simpson JJ. 1977. Irradiance measurements in the upper ocean. J Appl Oceanog 7: 952-956. Rossow WB and Zhang Y-C. 1994. Calculation of surface and top of atmosphere radiative fluxes from physical quantities derived from ISCCP datasets, 2. validation and first results. J Geophys Res 100: 1167-1197. Russell GL and Lerner JA. 1981. A new finite-differencing scheme for the tracer transport equation. J Appl Meteorol 20: 1483-1498. Russell GL, Miller JR, and Rind D. 1995. A coupled atmosphere-ocean model for transient climate change studies. Atmos-Ocean 33: 683-730. Russell GL, Miller JR, Rind D, Ruedy RA, Schmidt GA, and Sheth S. 2000. Comparison of model and observed regional temperature changes during the past 40 years. J Geophys Res 105(D11): 14891-14898. Russell GL, Miller JR, and Tsang L-C. 1985. Seasonal oceanic heat transports computed from an atmospheric model. Dyn Atmos Oceans 9: 253-271. Sausen R, Barthel K, and Hasselmann K. 1988. Coupled ocean-atmosphere models with flux correction. Clim Dyn 2: 145-164. Semtner AJ. 1974. An oceanic general circulation model with bottom topography. Tech Rep No 9, Dept Meteor, University of California, Los Angeles, 99 pp. Semtner AJ. 1987. A numerical study of sea ice and ocean circulation in the Arctic. J Phys Oceanogr 17: 1077-1099. Takano K and Wurtele MG. 1982. A fourth order energy and potential enstrophy conserving difference scheme. AFGL-TR-82-0205 (NTIS AD-A126626) Air Force Geophysics Laboratory, Hanscom AFB, MA 01731, 85 pp. Washington WM and Meehl GA. 1989. Climate sensitivity due to increased CO: Experiments with a coupled atmosphere and ocean general circulation model. Clim Dyn 4: 1-38. Zhang Y-C, Rossow WB, and Lacis AA. 1994. Calculation of surface and top of atmosphere radiative fluxes from physical quantities based on ISCCP datasets, 1. method and sensitivity to input uncertainties. J Geophys Res 100: 1149-1165. Table 1. Atmospheric prognostic variables. MA (kg/m) = mass per unit area in each atmospheric box UA (m/s) = eastward velocity component VA (m/s) = northward velocity component H0M (J) = potential enthalpy in each atmospheric box HXM (J) = west to east gradient of potential enthalpy HYM (J) = south to north gradient of potential enthalpy HZM (J) = upward vertical gradient of potential enthalpy Q0M (kg) = water vapor in each atmospheric box QXM (kg) = west to east gradient of water vapor QYM (kg) = south to north gradient of water vapor QZM (kg) = upward vertical gradient of water vapor Table 2. Ocean prognostic variables. MO (kg/m) = mass per unit area in each ocean box UO (m/s) = eastward velocity component VO (m/s) = northward velocity component G0M (J) = potential enthalpy in each ocean box GXM (J) = west to east gradient of potential enthalpy GYM (J) = south to north gradient of potential enthalpy GZM (J) = downward vertical gradient of potential enthalpy S0M (kg) = salt in each ocean box SXM (kg) = west to east gradient of salt SYM (kg) = south to north gradient of salt SZM (kg) = downward vertical gradient of salt Table 3. Prognostic variables for land reservoirs, for sea ice or lake ice reservoirs, and for lake and river mass MGL (kg/m) = Mass of ground liquid water MGS (kg/m) = Mass of ground solid ice and snow HGR (C) = Temperature of ground MGI (kg/m) = Snow mass above glacial ice HGI (J/m) = Heat, including latent energy, of glacial ice RSI = Horizontal ratio of sea ice to ocean MSI (kg/m) = Sea ice mass and snow mass above sea or lake ice HSI (J/m) = Heat, including latent energy, of sea or lake ice MR (kg) = Lake and river mass above the sill depth Table 4. Computer time percentages for different processes. Atmospheric dynamics 33.7 Ocean dynamics 26.2 Moist convection and condensation 10.6 Radiation 12.1 Surface fluxes 2.7 Ocean source terms 13.3 Diagnostic subroutines 1.3 Writing information to disk .1 Table 5. Prognostic variables for straits. MUST (kg/s) = mass flux in and out of the strait G0MST (J) = potential enthalpy in the strait GXMST (J) = horizontal gradient of potential enthalpy GZMST (J) = downward vertical gradient of potential enthalpy S0MST (kg) = salt in the strait SXMST (kg) = horizontal gradient of salt SZMST (kg) = downward vertical gradient of salt Table 6. Name, location and parameters in meters for the 12 straits for the 4 by 5 standard resolution. DISTPG is equal to Length plus the distances from the strait ends to the grid box centers. Strait From To Depth Width Length DISTPG Fury & Hecla 86W,73N 96W,67N 30 20000 489935 926089 Nares 71W,79N 64W,81N 158 50000 146155 510600 Gibraltar 6W,35N 1E,37N 158 25000 497073 1066702 English 1W,51N 1E,53N 30 35000 56153 561589 Kattegat 9E,59N 16W,57N 30 35000 56153 561589 Bosporous 29E,39N 31E,41N 30 35000 56153 561589 Red Sea 39E,21N 41E,19N 249 250000 343282 686499 Bab-al-Mandab 44E,17N 46E,15N 249 25000 69570 695673 Hormuz 54E,25N 56E,19N 30 50000 543944 1088514 Malacca 99E, 5N 106E, 3N 57 50000 1041765 1610458 Korea 129E,35N 131E,37N 98 170000 253184 632921 Soya 139E,43N 141E,45N 30 40000 119698 598492 Figure 1. Number of ocean vertical layers for each horizontal grid column for the standard 4 x 5 resolution. An A corresponds to 10 ocean layers, B = 11, C = 12, D = 13. Figure 2. Location of the u and v velocity components (arrows) with respect to the mass field (dots) for the C grid near the North Pole. Figure 3. A quadratic polynomial (m) (solid curve) and its least square fit line (dashed line) between m-m and m+m. The mean value of (m) in the interval is equal to [(m-x) + (m+x)] where x = m/12. x does not depend on the coefficients of . Figure 4. Diagram of the linear upstream scheme in one dimension (from Russell and Lerner [1981]). A is the mass of air entering grid box i from the left and A' is the mass of air leaving grid box i to the right. The heavy solid lines show the piecewise linear distribution of tracer concentration at the beginning of a time step. The dashed line segment is the least square fit line of tracer concentration in the interval (-M/2-A, M/2-A'). It determines the linear distribution in grid box i at the end of the time step.