4. Numerical techniques



4.1 Thermodynamic and stream function equations (temp.f, strcal.f).

The potential temperature thermodynamic equation is solved by first applying standard 2nd order finite difference approximations to the spatial derivative.The solution and the resultant system of ODE's is approximated with a spatially split, fourth order, multi-step, backward difference scheme from Van Der Houwen (1980). In this implementation the predictor value is taken as the last known solution. The split, discrete ODE system is iterated exactly four times. When starting a calculation without using prior results the first three solutions are obtained using a simple backward Euler temporal treatment rather then the full fourth order backward formula. The resultant linear system is solved directly by a specialized banded matrix solver 'band_slv.f'.

The stream function equation, which is a 2nd-order partial differential diagnostic equation is discretized by a center finite differencing scheme. The Jacobian matrix system of equation is solved using the band sysem solver in the file 'band_slv.f'.

4.2 Chemical transport equation

The chemical transport equation (Eq. 79) can be written in a numerical operator form each representing contribution from advection, diffusion, and chemistry:

eq108

where Ladv is the advection operator, Ldiff is the diffusion operator, and Lchem is the chemical operator. We solve for the discrete analogue to Eq. 108 by:

eq109

where Phi1 = PhiSLT is the discrete analogue to Ladv, Phi2 = Phichem is the discrete analogue to Lchem , and Phi3 = Phidiffis the discrete analogue to Ldiff . The following sections detail the numerical schemes used by each discrete operator.

4.2.1 Advection: Semi-Langrangian transport (slt.f).

The solution for the transport of chemicals by advection is obtained with the semi-Langrangian technique (Robert, 1981). This method requires a backward tracing of air parcels utilizing the residual wind (v*, w*) information, with the assumption that the tracer mixing ratio value is conserved along the trajectory. A 1 day timestep is the default for the back tracing calculation in the standard model run. Since the point of departure of a tracer parcel generally does not coincide with the model grid point, interpolations of the residual wind and the tracer mixing ratio value are necessary. The interpolation is done in subroutine LE3. To initialize for the diurnal chemistry, only the long-lived species at noon time are advected.

4.2.2 Chemical loss and production (chem_resol.f)

The changes in the concentration of long-lived and intermediate species by photochemical sources are calculated by the Euler-backward iterative method:

eq110

where delta.uct is the time step. The method is iterated (for a maximum of 10 times) until it is converged within 0.1%.

A different numerical method is used for selected species ( O(3P), O3, OH, HO2, NO, NO2, NO3, H) during night time (night_sub.f). They belong in a special category due to large diurnal variations in their chemical lifetimes. These species which are in photochemical equilbrium in the presence of sunlight must be explicitly time integrated during the night. The system of non-linear chemical equations are transformed to a system of non-linear algebraic equations through an implicit Euler time scheme with Newton Raphson iteration:

eq111-133

where F'(X) is the Jacobian Matrix (dF/dX ) and it is the iteration number. The timesteps for the chemical equations varies according to the length of day and night, with 8 timesteps every 24 hours. For the timestep during daylight hours, the length of daylight is divided by 4; likewise, night timestep is the duration of the night divided by 4.

4.2.3 Diffusion (diffus.f)

The molecular and eddy diffusion component of transport is solved like the potential temperature by first applying standard finite difference analogues to the spatial derivatives. The resultant linear ODE system is solved through a straight forward application of backward Euler with modifications to incorporate splitting along the spatial coordinates (y and z). The combined splitting, backward Euler algorithm is iterated three times.