The dynamical wave forcing and associated diffusion in SOCRATES are estimated as a function of the zonal wind and derived according to linear wave theories. Although it is still necessary to parameterize wave processes, the wave components are able to respond to and interact with the mean field. Three types of waves are included in the model: planetary, gravity and tidal waves. Transience, dissipation, and breaking of these waves produces convergence of the wave fluxes and deposit momentum on the zonal wind field, which is reflected in the governing equations (Eq. 2). Wave diffusive mixing is also generated in the process and principally affects the temperature distribution (Eq. 1) and chemical transport (Eq. 79). Details of how the wave forcing and diffusion are derived is described next.
2.2.1 Planetary waves
Following Smith and Avery (1987), the vorticity and thermodynamic equations for the planetary wave take the forms:
The potential vorticity q is:
and the wave component of the meridional wind is:
where 
' is the geopotential disturbance, with prime referring to wave quantities. To represent wave dissipation from Newtonian cooling, Rayleigh friction and wave breaking, a linear damping term is added to the right hand side of equations (39) and (40), where 
is the Newtonian cooling rate, 
is the Rayleigh dissipation rate, and 
is the dissipation rate from planetary wave breaking (method to derive d will be described later). To derive the potential vorticity equation, w' is eliminated between (39) and (40), and substituted with (41) and (42) to form a single equation of the planetary wave geopotential '. For a more concise formulation of the potential vorticity equation, the following transformed coordinates are defined:
Note that the log pressure z here is a dimensionless altitude unit used only within the wave model which is different from the pressure altitude coordinate used throughout SOCRATES.
Furthermore, a new variable 
is introduced to define the geopotential ' at zonal wave number m, such that
Thus, the equation can be written as:
In the current version of the model, we are using s = 0.02, independent of altitude, and only planetary wave one is considered, with m = 1. Then applying implicit method for time integration, Equation (43) can be rewritten as:
Using a space-centered differencing scheme, Equation (44) is solved using a numerical technique developed by Lindzen and Kuo (1969). Forcing at the lower boundary (16 km) is specified according to climatological values of wave one geopotential height at 100 mb (Randel, 1987). From the wave stream function, the wave zonal wind u' and meridional wind v' can be derived. The planetary wave momentum flux divergence, or the Eliassen-Palm (EP) flux divergence 
F , is then calculated from
The planetary wave diffusion coefficient Kyy and the wave breaking dissipation rate d are derived according to methods developed by Garcia (1991). Garcia adopted barotropic instability as the criterion for planetary wave breaking, and assumed that there is dissipation when |q'y| 
qy, generally when the zonal mean potential vorticity meridional gradient becomes negative. In this case, the dissipation rate is calculated from the WKB approximation [Eq. 17 of Garcia (1991)]:
where cgy and cgz are the horizontal and vertical group velocities given by
The planetary wave diffusion coefficient Kyy coresponding to this dissipation rate is derived from
Evaluation of the dissipation rate requires knowledge of the phase speed (c), meridional wave number (l), and vertical wave number (m) for the group velocities. These quantities are likewise derived from Garcia (1991) [Eqs. (35) and (36)]:
where Im denotes the imaginary component.
2.2.2 Gravity waves
In the new version of the model, two options are available for the parameterization of the gravity wave breaking process. The default standard run uses the Lindzen-Holton gravity wave breaking scheme (Lindzen 1981, Holton 1982) for estimating the gravity wave forcing and vertical diffusion coefficient. In this case, the gravity wave momentum flux divergence is expressed as:
and the vertical eddy diffusion resulting from wave breaking is given by
where Pr is the Prandtl number and Ag is the amplitude coefficient. These parameters are unchanged from the previous versions of SOCRATES (Brasseur and Hitchman 1987, Brasseur et al. 1990).
A second option is to use the parameterization scheme developed primarily by Fritts and Lu (1993) (hereafter FL93), which will be described in more detail here.
The Fritts-Lu gravity wave forcing parameterization scheme uses a general gravity wave spectral formulation (Fritts and VanZandt 1993) provided by the observed gravity wave spectrum to deduce how the gravity wave energy flux responds to variation in the background environment. The variation of total gravity wave energy E0 with height is given by E0 
dz/HE, where the scale height for energy growth is specified as
with z in kilometers. Equation (50) is slightly different from Equation (28) of Fritts and VanZandt (1993), based on more recent observations of the total energy variations (Wilson et al. 1991a,b), and produces a more realistic estimate of the gravity wave forcing profile. The lower boundary condition for gravity wave energy forcing is currently set at 30 km, and the total energy at 30 km is specified as a function of latitude and season according to Hirota (1984). With the assumption that the total energy is comprised of two component spectra propagating east and west, the energy at the lower boundary is equally split into the east and west components, while their intrinsic phase speeds are calculated from the relation [Eq. (22) of Fritts and VanZandt (1993), E0/2]
where j denotes the east or west direction, and c*j is the intrinsic phase speed.
To determine how the total energy varies with height, with the atmospheric stratification and the local mean wind, the relationship expressed in Eq. (31) of FL93 ( E0 N1/2ez/HE ) is used to obtain
where the superscript z in this case denotes the altitude grid. Meanwhile, the energy spectra for each directional component ( component energy ) are assumed to vary with changes in the zonal wind in the vertical direction. When the vertical gradient of zonal wind is positive, the westerly phase speed c*w increases with height according to
with the corresponding easterly component energy Ee calculated from Equation (51). With Ee known, Ew (westerly component energy) can be calculated from E0 - Ee. Similar conditions are applied when the vertical gradient of the zonal wind is negative, in which case c*e increases with altitude.
Aside from the constraint of energy (or wave saturation) at large vertical wave number that arise from variation of wave energy with height and static stability as mentioned above, another saturation criteria set by FL93 is provided by the maximum anisotropy () allowed from wave propagation. Fritts and Lu assumed that no less than a fraction of the energy may be associated with any one component (Ej/E0 < 1 - ), which is motivated by the insensitivity of motions of small vertical wave number to mean wind changes and the tendency of instability to occur under large wind shear conditions. In the model, we have set a=0.25. When conditions are such that Ej/E0 < 1 - , than Ej is adjusted back to Ej = E0(1 - ).
With the variation with altitude of the two energy components calculated, the momentum fluxes and its vertical divergence can be estimated. The relationships between energy (Ej), the momentum flux (FP), and energy flux (FE) are expressed according to FL93
The gravity wave momentum flux divergence FG and the energy dissipation rate are calculated from
With the energy dissipation rate (
) known, the gravity wave diffusion coefficient (Kzz) is calculated according to the relationship
As discussed in Weinstock (1984), the value of the coefficient A (equivalent to the inverse of Prandtl no.) could be as small as 1/3 or even less. In the model, we adopt the value of A = 0.2.
2.2.3. Tidal waves.
The parameterization of tidal wave forcing and diffusion in the 2-D model is based on the work of Lindzen (1981). Lindzen showed that under the WKBJ approximation of the linear wave equation, the momentum flux divergence from tidal gravity wave can be expressed as
and the diffusion coefficient can be written as
From theoretical calculations, it can be shown that only the first diurnal mode of the tidal waves within 30o latitude is able to propagate up into the mesosphere, and break to generate forcing and turbulence. At 85 km altitude where the waves break, with the tidal wave phase speed of 465 m/s and k = l/6400 for the main diurnal propagating mode, FT = -16.3 m/s/day, and Kzz = 183 m2/s, according to Eqs. (59) and (60). Above 85 km, we let FT and Kzz decrease with altitude to take into account the rapid increase of the contribution from molecular viscosity and conductivity. At level z (km), above the breaking level, the momentum flux divergence FT and diffusion coefficient Kzz calculated at 85 km are multiplied by exp[ (85-z)/5]. Since waves are still expected to generate turbulence below the breaking level, a factor exp[(z-85)/5] is applied to FT and Kzz below the breaking level (85 km).
Because the tidal wave mode that propagates up to the mesosphere is restricted to latitudes -30o < lat < 30o, a factor of exp[-(lat/20o)2] is applied to FT and Kzz beyond the latitudes 30o N and 30o S. Within the latitude of 30o, FT and Kzz are assumed to be constant with latitude.
2.2.4. Tropospheric wave forcing and latent heating specification.
As mentioned in the Section 1, the lower boundary condition for the circulation is applied at 2 km. In order to solve the stream function equation (7), the tropospheric wave momentum flux divergence and the latent heat release rate that is particularly important in the troposphere, needs to be accounted for in the forcing term CF. In SOCRATES, we chose to specify the EP flux divergence and latent heating rate in the troposphere according to a seasonal monthly-mean climatology.
Based on the EP flux divergence climatology established in Randel (1992), we formulated the EP flux divergence as a function of latitude, height and season as following
If -
/2 < 
< -/6 or /6 < < /2 (mid to high latitudes):
If -/6 < < /6 (tropics):
EP flux divergence = 0.
The amplitude of Atrop varies with season and hemisphere:
/2 < <-/6/6 < < /2 /2 < <-/6/6 < < /2 /2 < <-/6/6 < < /2 The tropospheric latent heating rate due to condensation during convection events is specified to vary with latitude, height and season, based on Boville (1985) and Peixoto and Oort (1991). We use the following formulation for condensation heating rate:
(z-3)/12) exp(-2/100) (K/day) (62a)
(z+10)/20) cos(9/4*) (K/day) (62b)
(z+10)/20) cos(3/2*) (K/day) (62b)
The latent heating rate QLH calculated from Eq. (62) is finally processed with the Gaussian filter smoothing routine, with weights of 0.25, 0.5, and 0.25, in both y and z direction for five times.
2.2.5 The circulation global mass balance correction.
The meridional circulation estimated from the stream function equation (7) is not necessarily in perfect global mass balance, either from inaccuracies in the radiative or dynamical forcing scheme, or from omitted physical processes such as small-scale turbulent processes. Thus, corrections have to be made to ensure that as the model integration proceeds, the integrated zonally-averaged mass flux between poles along a pressure surface remains equal zero. The global mass balance requirement can be written as
which is the latitudinal integration of area-weighted vertical residual wind along the log pressure height surface. Since w* does not necessarily meet the condition as shown in Eq.(63), a correction is applied to w* such that
where is a function of height only. Thus, the correction factor can be obtained from
2.2.6 Molecular diffusivity parameterization.
The effect of molecular diffusion on the chemical species and thermal field becomes important above the mesopause. Since the model domain of SOCRATES extends up to 120 km, the effect of molecular diffusion effect is taken into account.
Following Banks and Kockarts (1973), the vertical flux of a minor constituent (i) relative to air from molecular diffusion can be written as:
with the molecular diffusion coefficient Km for a particular gas i expressed as:
where [n] is the number density of the minor consituent, H(i) = RT/(mig) is the scale height for trace species (i), T the thermal diffusion factor (-0.38 for H), mair is the molecular weight of air (g/mole) which is specified as a function of altitude, mi is the molecular weight of gas i (g/mole), and [M] is the number density of air (cm-3). To derive the diffusivity flux formulation of Eq. (66) consistent with the chemical transport equation that will be described in section 2.3 [Equation (79)], [n] in Eq. (66) is substituted by the mixing ratio Xi.:
where
Thus, if we define the vertical diffusive velocity by
the diffusivity flux can be expressed as
Since the time rate of change of mixing ratio expressed in flux form is
the effect of molecular diffusion on the chemical transport is given by
where Xi represents the zonally averaged mixing ratio of gas i.
The thermal conductivity coefficient KT associated with the molecular diffusion of heat in Equation (6) is adapted from the NCAR Thermospheric General Circulation Model (Alan Burns, private communication), in the form of
where X represents the mixing ratio. To convert KT to units of m2/s for the model, Equation (72) is multiplied by a factor of (mairg2H2)/(P0cpRTe-z/H).
2.2.7. Quasibiennial oscillation simulation
The quasibiennial oscillation (QBO) of the temperature, zonal wind, and meridional circulation in the equatorial stratosphere is simulated in the 2-D model according to Politowicz and Hitchman (1997). An analytical expression for the QBO radiative forcing is derived from a composite of observed zonal wind in the tropics, which is used to obtain the QBO forced circulation and the temperature oscillation in the stratosphere. Since the 2-D model does not explicitly deal with the zonal component of the momentum equation, using a QBO type radiative forcing in the stream function equation (Eq. 7) is a convenient way to simulate QBO processes for temperature and circulation without changes in the logistics of the model physics. The QBO forcing is specified and therefore does not respond to changes in zonal mean conditions. However, with this parameterization, the model can be used to examine various QBO radiative, transport, and chemical feedback effects on stratospheric chemistry.
From radiosonde time series at 15 tropical stations, Politowitz and Hitchman (1997) created a composite QBO for zonal wind at tropical latitudes between 16 to 40 km altitude with the formulation:
where hqbo is the QBO vertical extent of 24 km and t is the model time (in units of day with tqbo in units of days). The maximum zonal wind QBO perturbation (umax ) is 25 m/s and is centered over the equator. The perturbation decays exponentially in latitude with a Gaussian latitude width ( ) of 15o, with negative tails beyond 15o on either side of the equator with the factor . The perturbation varies sinusoidally in altitude, between 16 km and 40 km, with maximum vertical amplitude at zmax, (28 km). A node between easterlies and westerlies descends with time, with a descent period of tqbo=27 months (or 821 days) and a half width (zhfwd) of 27 km.
A analytical expression for the time rate of change of temperature is derived based on the equatorial beta plane thermal wind relationship:
Assuming T = T0exp(-2/w2) , where T0 is the equatorial temperature with w in radial units, the thermal wind relation can be written as
Substituting the Eq. 73 into Equation (76), the time rate of change of temperature as a function of time, latitude, and altitude can be written as
Although in principle, temperature QBO can be derived from Equation (75) knowing the zonal wind profile, the temperature is actually derived from dynamical and thermal equations through the QBO forced meridional circulation. One option is to include the time rate of temperature change from Equation (76) as a diabatic heat forcing in the stream function equation (7) to obtain a meridional circulation forced by the QBO temperature damping rate. Another option is to include in the stream function equation a direct QBO momentum forcing with the formulation:
where Aqbo is an adjustable parameter to optimize the results. The QBO forced circulation is then used to integrate in time the thermodynamic equation (1), which provides the temperature. Through this procedure, both the circulation and temperature can be affected by the QBO.