## 5.6 Discretization     continued...

Scalar conservation equations (such as for turbulence K.E.) can also be reduced to finite difference form in a similar manner:

Similar equations can be written down for , S and . The momentum equation (70) can be written as:

Therefore the finite difference form of the momentum equation is

Notice that a pseudo-Coriolis parameter can be defined in an orthogonal curvilinear coordinate system:

accounts for the acceleration terms arising out of the changes in grid size from grid point to grid point. The finite difference form of the momentum equation (71) is

Note that for convenience and algebraic clarity we have avoided writing down the terms on the R.H.S. of equations (103) and (105).

The external mode (barotropic) equations (90) to (92) are also cast in similar finite difference forms, which, for brevity, we are not going to write down. The resulting external and internal mode finite difference equations are solved in a split-mode method, as discussed earlier, to obtain solutions for the prognostic quantities as a function of time. The external mode is solved using explicit differencing techniques, while the internal mode equations are solved implicitly in the vertical direction, but explicitly in the horizontal. The finite difference equations for the internal mode are cast in a tri-diagonal matrix form, which are then solved using well-known techniques such as the Thomas tri-diagonal algorithm (see chapter 6.2.2 volume 1 of [14]). Close synchronization between the two modes, external and internal, is assured by passing onto the external mode, from the internal mode, the bottom friction and the component of the pressure gradient that is due to the density field. In turn, the external mode supplies the internal mode with the sea surface elevation and vertically averaged currents that it needs. This interaction takes place at each internal time step. The procedure then is to start with one external mode (barotropic) calculation and then pass on the needed information to the internal (baroclinic) mode and take one time step with it. The internal mode then passes on the needed information to the external mode and the whole process is repeated until the end of the calculation is reached.

As mentioned in earlier sections, the data to initialize these models come from historical archives, such as the National Oceanographic Data Center (NODC), and forcing data for the model can be from historical archives at National Oceanic and Atmospheric Administration (NOAA) or Numerical Weather Prediction (NWP) Centers around the world, such as the National Meteorological Center (NMC) in the U.S. or European Center for Medium range Weather Forecasting (ECMWF) in Europe. Excellent bathymetric data for the world, including the oceans, at about 9 km resolution, called ETOP05, is available from the National Center for Atmospheric Research (NCAR). From these archives, it is possible to derive the bottom topography, as well as the initial 3-D temperature and salinity fields at all model grid points. This involves considerable processing of the data bases to remove data errors and other uncertainties and obtain meaningful fields to initialize the model. The same applies to deriving the momentum and buoyancy flux fields at the ocean surface from NWP or NOAA products to drive the model. For the sake of brevity and clarity, we will not detail these procedures. Suffice it to say that data preparation (and analyses of the model output) is a very significant fraction of any ocean modeling effort.

It is important to realize that we have described only the dynamical oceanographic component of the model. In many practical applications such as investigating the fate of man-made carbon dioxide, a large fraction of which is absorbed and sequestered in the global oceans, we need to consider and solve equations for the chemical and biological variables also. Once again, this involves writing down the governing partial differential equations (PDEs) very much similar to the governing equations for and S (equations (38) and (39)) for the quantities of interest, casting them in finite difference form as described above, and solving them numerically. For example, for investigating the inorganic carbon cycle in the ocean, where the carbon dioxide is absorbed at the sub-polar/polar latitudes (and expelled in warmer tropics) and converted to the carbonate, bicarbonate and dissolved gas forms, it is necessary to solve for the total carbon and the pH of the water parcel. The governing equations for these are similar to the one for salinity, except that the right hand side now involves source and sink terms. Also an equation involving chemical reaction rates that results in an appropriate balance between the three forms of dissolved carbon needs to be solved. Nevertheless, the same methodology can be used to explore how the oceans absorb the excess carbon mankind is introducing into the atmosphere by fossil-fuel burning. This problem is of obvious importance to the future climate our children will live in, as a result of our massive uncontrolled tinkering with our global environment by profligate use of fossil carbon that took nature millions of years to create.

The carbon cycle in the oceans involves another component, the organic cycle. Phytoplankton (microscopic plants) in the upper layers of the ocean take up the dissolved carbon dioxide within the photic zone (the upper layers where there is sufficient light for photosynthesis), and are in turn consumed by zooplanktonic organisms, whose carbonate shells rain down onto the ocean floor upon their death and remain undissolved and hence sequestered from the atmosphere for thousands of years. Thus the amount of carbon sequestered by the ocean and thus not contributing to an immediate global warming, is a strong function of the biological productivity of the upper layers of the ocean. However, planktonic organisms need for growth nutrients such as nitrates and phosphates, which have to come from the deeper ocean. Thus both light and nutrients are needed for biological productivity, and it is only a fraction of the oceanic upper layers (a few per cent) around the globe that are responsible for more than 90% of the biological productivity and the food chain so crucial to the growing population on the globe. These are the so-called upwelling regions, mostly around coasts where the winds blow such that the ocean waters are expelled from the coast and deeper, nutrient-rich waters come up to replace them. If conditions are right, as for example the spring--summer time, biological productivity literally explodes.

To account for the organic cycle, it is necessary to solve for the chemical and biological variables involved, such as the nutrient concentrations, and the phytoplankton and zooplankton densities in the water column, in addition to solving for the circulation. Once again it is a matter of solving the appropriate PDEs cast in finite difference form as above. The only difference, of course is that these equations involve complex source/sink terms on the right hand side that need to be appropriately parameterized, in view of our knowledge of how the planktonic organisms behave. Thus while in theory at least, it is rather straightforward to solve for the physical, chemical and biological state of our oceans, practically speaking, the solution becomes extremely complex, time-consuming and computationally very expensive. It is not unusual for a relatively high resolution (50 km grid size) physical global ocean model itself to take several hundred hours on a modern computer such as the CRAY YMP, just to do a year's worth of model simulation. It is easy to imagine therefore a full physical--chemical--biological model of the global oceans requiring the next generation of computers, the massively parallel ones of the CM-5 type. In fact, understanding and predicting the state of the global oceans is one of the ``grand challenge'' problems that require the power of the Teraflop massively parallel computers that are expected to become available early in the next century.

(See exercise 7.)