A NUMERICAL SOLUTION TO THE RADIAL EQUATION OF THE TIDAL WAVE PROPAGATION

The tidal vave function y(x) is a solution to an inhoraogeneous, linear, second-order differential equation with variable coefficient. Numerical values for the height-dependence terms, In the observed tides, have been utilized in finding y(x) as a solution to an initial-value problem. Complex Fast Fourier Transform technique is also used to obtain the solution in a complex form. Based on a realistic temperature structure, the atmosphere below 110 km has been divided into layers with distinct characteristics, and thus the technique of propagation in stratified media has been applied. The reduced homogeneous equation assumes the form of Helmholtz equation and with initial conditions the general solution is obtained, HIRAHARE TRIESTE August 1981 * To be submitted for publication. + To be presented at XV Biennial Fluid Dynamics Sympos-iura, to Be held in Poland 6-12 September 1981. ** Permanent address: Applied Math. Dept., Faculty of Science, Ain Shams Univ. Cairo, Egypt. In the development of the tidal theory(Siebert, 1961 and Chapman and Lindzen, 1970), the assumption is made such that the dependence of the tides on height, Co-latitude and longitude rf> are separable. The latitude-dependences of the tidal fields are expressed in terms of the Hough functions (H ) ( 6 ) , vhlcfi. are solutions of Laplace's tidal equation. It is an eigenvalue-eigenfunction problem and has been the subject of extensive study during the Ia3t decade CKato, 19-66; Lindzen, 1966; Longuett-Higgins, 1967 and Makarious, 1968). The vertical structure of the tidal wave propagation, on the other band, is expressed in terms of the wave function y(x), which is the solution of the Radial Wave equation. For a given mode of oscillation, the variable coefficient, in that equation, depends solely on the vertical temperature distribution. Previous investigators tPekerls, 1937; Jacchia and Kopal, 1952; Butler and Small, 1963 and Lindzen; 1967) adopted,however, simplified vertical temperature structure in order to render the mathematical treatment more tractable. Consequently classical tidal theory has limitations, arising out of the assumptions on which it is based, which need to be recognized when making close comparison with the observational data (Groves, 1977). With a realistic temperature structure, however, no closed form solution to the vertical wave equation exists and the solution has to be approached numerically. This is the goal of the present study. Analytical expressions for the general solution, in conformity with the observed tidal oscillations In the atmosphere, have been obtained. The numerical solutions for the oscillatory and trapped modes of tidal wave propagation have also been presented and discussed.

Numerical values for the height-dependence terms, In the observed tides, have been utilized in finding y(x) as a solution to an initial-value problem.
Complex Fast Fourier Transform technique is also used to obtain the solution in a complex form.
Based on a realistic temperature structure, the atmosphere -below 110 km -has been divided into layers with distinct characteristics, and thus the technique of propagation in stratified media has been applied.The Cairo, Egypt.
In the development of the tidal theory- (Siebert, 1961 andChapman andLindzen, 1970), the assumption is made such that the dependence of the tides on height, Co-latitude and longitude rf> are separable.
The latitude-dependences of the tidal fields are expressed in terms of the Hough functions (H)(6), vhlcfi.are solutions of Laplace's tidal equation.
The vertical structure of the tidal wave propagation, on the other band, is expressed in terms of the wave function y(x), which is the solution of the Radial Wave equation.For a given mode of oscillation, the variable coefficient, in that equation, depends solely on the vertical temperature distribution.Previous investigators tPekerls, 1937; Jacchia and Kopal, 1952;Butler and Small, 1963 and Lindzen;1967) adopted,however, simplified vertical temperature structure in order to render the mathematical treatment more tractable.
Consequently classical tidal theory has limitations, arising out of the assumptions on which it is based, which need to be recognized when making close comparison with the observational data (Groves, 1977).With a realistic temperature structure, however, no closed form solution to the vertical wave equation exists and the solution has to be approached numerically.This is the goal of the present study.
Analytical expressions for the general solution, in conformity with the observed tidal oscillations In the atmosphere, have been obtained.The numerical solutions for the oscillatory and trapped modes of tidal wave propagation have also been presented and discussed. -2-2.

THE RADIAL EQUATION
The Baste equations of the theory of atmospheric tides are given in Siebert Cl9^l) and the derivations of the bas-ic equation* need not, therefore he repeated.The wave function y tx) is expressed In terms of the velocity divergence x (x) and the rate of heating J Cx) per unit mass as where x, the reduced height, is defined by f&7 (2) hn is the equivalent depth of the mode, it is the eigenvalue obtained through the solution of Laplace's tidal equation, wfiicB.depends only on tHe period of the mode.
In order to evaluate the reduced height tEq.2) realistic temperature data are utilized.These data are presented "by Cole et al (1965) for 0 < Z < 30 Fa and t>y Groves (l97l) for 30 < Z < 120 Km -at 5 Ki n intervals.
The scale height H(x) has been obtained, by using quadratic interpolating formula, at x = 0[0.l)l7.By least-squares approximation H{x) has been found to he represented fairly accurately by a fourth degree polynomial in the region 0 < x < 15.6.At x > 15.6, a constant temperature gradient is assumed, i.e.  la.
As a second-order differential equation, two boundary conditions for y*(x) must be satisfied.At high levels we shall require, first, that y be bounded as x •+ » and secondly, that the flow of energy be in the upwara direction.Following Wilkes [19^9) and according to Jacchia and Kopal (1952), this requirement should be inforced not at infinity, but at a finite altitude Z*, beyond vhich the highly rarefied atmosphere can he expected to float more or less passively on the less dense strata.Z* is taken, in the present study, * The subscript n will be dropped thereafter for simplicity In writing and we deal with each mode independently.
-h-to correspond to 110.Kf c and to conform \itn.the toundness of y-It la taken As a lover boundary condition, ±t is generally assumed that at the earth's surface, the vertical omponent of the velocity must vanish (Siebert,196l); this implies that: The conventional tidal fields may all 'be conviently expressed in terms 'of the y's and (HYS:

VERTICAL DEFEHDENCEE OF THE TIDAL FIELDS
The beat tidal observations are those for surface pressure, since they are obtained for a long period of time.The amplitudes and phases of the solar diurnal and semidiurnal oscillations have been calculated by  and Haurwitz (1956) respectively and are presented in terms of the associated Legendre functions as: -5- using the lover boundary condition (Eq.9), Tne coefficients y^' n have been previously obtained through the solution of Laplace's tidal equation.(Makarious, 1968).Therefore the tvo initial values y(o) and y'(0) have been evaluated at x ™ 0.
The vertical dependence terms in the horizontal vind fields (Eqa.10, 11) may be expressed as: ( 37 --r Cl6) This term had been'previously evaluated (Makarious, 1979) in analysing 6 years of horizontal vlnd data and the results were obtained at Z =• 3.5 (2.5) 60 Kin.
By Lagrangian quadratic interpolation, and using the lower boundary condition (q.), values of z(x) at x = 0 (O.l) 8.1t are evaluated.
The solution y(x) is expressed as a complex variable function, whose real and Imaginary parts are evaluated Independently.
The initial^yalue pro'Slem for the inhoroogeneoua.first order D.E is to determine ytcl that sfitiafi.es-Cl6"}K4tn.t&e initial value ylfl) as given in tl5).Multiplying Both, sides By the integrating factor exp(-x) and integrate "Eiotb aides of the resulting equation we get: ClT) The integral on the right-hand side has been evaluated using Simpson's rule of integration with the spacings Jx = 0.1.
By inspection, it was found that z(x) can "be approximated fairly accurately by Complex Fourier Series expansion in the form: (18) The coefficients C are ottained for the two regions 0 < x < h.Q and h.B < x < 9.6 using Complex Fast Fourier Transform Technique.The two regions are chosen this way as they are characterized by the different temperature behaviours and consequently different wavelengths for each mode.
The periodic values of z{x), expressed by Eq.Cl8), are therefore included in the integral of ClT) and solution y(x) is also obtained.The two obtained solutions are found to be comparable,within differences ±0-5ô n the average, reflecting the validity of the approximation (l8) and the numerical stability of the solution.The initial-value problem is considered to be well-set in the domain 0< x < 9.6 and thus the solution y(x) obtained in (IT) is the only solution for each given y(0) at x = 0.
The inhmogeneous term ztx) (Eq.16) is the forcing function and corresponds to the applied tidal force.Differentiating both sides of (l6) ve get: -T-Cl9) In comparing U?) vith. the inacmogeneous D.Eq.C3l, the source function QU) is expressed analytically, in terms of tne solution y(x), as: h.

THE REDUCED HOMOGENEOUS EQUATION
If QCx) = 0, the inhomogeneous second-order differential equation ( 3) reduces to the homogeneous form: An examination of this reduced equation ( 21) shows that there is a simple analogy with the propagation of plane waves in a medium of varying refractive index Cwilfces, 19 1 »9).This state of affairs is described by the Helmholtz equation, which, in general, is a typical one for oscillating systemB.
In order to obtain exact solutions to the Helmholtz equation, the variable coefficient u (x) can be replaced by a collection of continuous functions which look like acceptable approximation to u U).This type of functions will allow us to investigate effects of the continuous stratification by means of simple solutions.Insight into the properties of such a media may be obtained by the use of a linear law for y 2 (x) in the Helmholtz equation. -8-

M w
In a preyiou?investigation*, the region neloy JC = 15.6 has teen unequally divided into six auBregions, In each, of vftich..

C22)
The present study extends-the previous one 'By dividing tne region x < 15-6 into equally spaced subregions-.The results for P are given in Table It for the different modes of diurnal and semidiurnal oscillations.It can be inferred that the atmosphere is divided into k distinct regions with regard to the characteristics of its refr.ictiieindex, namely (i) 0 < x < It,8, Cii) U.8 < x < 9.6, (ill) 9.6 < x < 12.0 ana (iv) 12.0 < x < 16.8.For the region 15.6 < x < 16.8, the reduced homogeneous equation assumes the Bessel form equation, and for other layers it assumes Stoke's form.
The homogeneous solution has then been obtained in terms of Bessel functions in the form: C23) The initial-value problem for the second-order homogeneous linear differential equation ( 21) is to determine the general solution y(x) given by ( 23) that satisfies the initial conditions yc(x ) and y '(x ).Since the two solutions y and We*,,) = are linearly Independent, then we get the Wronskian: In order to get y tx ) and y ' (x ) ire call for the condition at CO CO XQ = 9.6; J is assumes to vanish at that freight, and thereafter, based on observations (Butler andSmall, 1963 andLlndzen, 1967).Therefore at x » 9.6, the general solution to the reduced equation ( 21) is equivalent to that obtained by ( 20) in setting Q = 0. Thus we get: Substituting in Eq. ( 25), ve can easily obtain the constants C^ and C2 and hence the homogeneous solution can then be evaluated numerically.This completes the method of numerical solution to the differential Equation (Eq.3).
Therefore we get for the two constants expressions: and (Eq.23) the following To be published. -9-5.

RESULTS AM) DISCUSSION
In this paper we presented a numerical method of solving linear inhomogeneous second-order differential equation vlth a variable coefficient, -10-aa applied to the vertical tidal*-w,ave propagation.Ttro solutions, are obtained, a particular solution that ts-derrved from the oBaerved tidal motion and a homogeneous solution using the characteristic* of the atmospheric structure.
As an illustration to the vertical Behaviour of the solutions, the results for two modes of the diurnal oscillation tl, -l) and tl, l) and for the main semidiurnal mode (2, 2) are presented xn Table 2.
It can be inferred from Table 2, that the particular solution grows exponentially vtth x and has greater values for the trapped mode (l, -l) than.the oscillating one, reflecting the effeetivennss of the former mode to the excitation.Although the (l,l) iiode is a propagating one, the constancy in phase is due to its large wavelength which is about U3 Kin in the troposphere and 25 Km in the stratasphere.For the main semidiurnal mode (2, 2), v (x) n ia almost zero through most of the stratasphere; i.e. extremely long vertical wavelength ( V150 Km).Thus not only does this mode receive the bulk of semidiurnal oscillation, but it must also respond to the excitation with particular efficiency.
With respect to the homogeneous solution, the function y(x) has been found to be non-oscillatory and the amplitude is exponentially increasing for 0 < x < k,B and decreasing above; in the latter region the energy is trapped; a feature associated with the negative mode.For the oscillating mode the function y(x) is in general oscillating and thus it is very sensitive to the variation in temperature.
To conclude, with such a method of solution, a general solution to the radial equation can be expressed analytically.It will be interesting to show hov this solution can be used to evaluate the thermal drive for the atmospheric tides.This will be the subject of the forthcoming work..73 .
Values obtained from the integration of z(x) as obtained from the observatIons.•* The subscript (s, n) stands for the type of oscillation (s » SL for migrating solar wave) and n is the mode number, (negative values stand for modes with negative equivalent depths). -ll*- Theoretical Physics, Trieste, Italy ABSTRACT The tidal vave function y(x) is a solution to an inhoraogeneous, linear, second-order differential equation with variable coefficient.

+
reduced homogeneous equation assumes the form of Helmholtz equation and with initial conditions the general solution is obtained, To be presented at XV Biennial Fluid Dynamics Sympos-iura, to Be held in Poland 6-12 September 1981.** Permanent address: Applied Math.Dept., Faculty of Science, Ain Shams Univ.
= RT(z)/g and k = 2/7 = {y -l)/y-The wave function yQCx) is a solution to the Radial Wave Equation, k + 0.25)/exp(3.9).Values of the approximated scale height §{x) as obtained from the assumed temperature model are given in Table

ACKHQNLEDGBtENT
The author would liRe to thank Professor AMua Salara, the International Atomic Energy Agency and UNESCO for hospitality at the International Centre for Theoretical Physics, Trieste.Thanks are due to the Centre's members of staff far their cooperation and for use of the computer facilities.He would also like to thank Frof. A. H. Cook of Cavendish Laboratory and Prof. D. Tritton of Newcastle University for fruitful discussions concerning the problem.-11--12-

TABLE 1 .
THE MODEL AIMOgPHSHE

TABLE 2 .
TEE MJtERICAL SOLUTION TO THE