TWO METHODS FOR THE NUMERICAL CALCUUTION OF ACOUSTIC NORMAL MODES IN THE OCEAN Ki rk Eden Evans Library Naval Postgraduate Scfxjof Monterey, California 93940 NAVAL POSTGRAOOATE SCHOOL Monterey, Galiiornia T TWO METHODS FOR THE NUMERICAL CALCULATION OF ACOUSTIC NORMAL MODES IN THE OCEAN j by Kirk Eden Evans Th esis Advisors G.H. Jung & A. B. Coppens September 1973 Tl 57081 AfpAxtvzd ^OA. pabtic -teXertie; dlst/uhatcon imtaruXtd. Two Methods for the Numerical Calculation of Acoustic Normal Modes in the Ocean by Kirk Eden Evans Lieutenant, United States Navy B.A., Miami University, 1966 Submitted in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE IN OCEANOGRAPHY from the NAVAL POSTGRADUATE SCHOOL September 1973 Library Naval Post.-^raduate Scnool Monterey. California ^3940 ABSTRACT Three computer programs were written to find the eigenvalues and eigenfunctions of acoustic normal modes in the ocean. The programs used two different methods: an iterative finite difference scheme, and a method based upon the WKB approximation of quantum mechanics. The raethods assume a flat fluid bottom and are designed for any arbitrary sound speed profile. While the results of both the finite difference and the WKB methods agreed, the WKB method proved faster. \ TABLE OF CONTENTS I. INTRODUCTION 11 II. THEORY 15 A. WAVE THEORY 15 1. The Wave Equation 15 2. Boundary Conditions 20 3. The Nature of the Solution 24 4. Mode Parameters 30 a. Phase and Group Speed 30 b. Bottom Decay Coefficient 30 B. THE WENTZEL-KRAMERS-BRILLOUIN (WKB) APPROXIMATION 31 1, The Characteristic Equation 31 2. Turning Point Connection 37 a. Asymptotic Method 38 b. The Airy Phase Method 40 III. THE PROGRAM METHODS 42 A. NORMOl - A FINITE DIFFERENCE TECHNIQUE "42 B. N0RM04 AND N0RM05 - THE WKB APPROACH - 42 IV. THE BASIC ALGORITHMS 50 A. NORMOl 50 1, General Description •-- 50 2. The Subroutines 52 B. N0RM04 AND N0RM05 54 1. General Description 54 2. SEARCH Subroutine -- 55 3. WKB Subroutine 57 V. RESULTS 58 A. TEST CASE ONE - NEGATIVE GRADIENT 58 B. TEST CASE TWO - POSITIVE GRADIENT 62 C. TEST CASE THREE - SYMMETRIC WAVE DUCT - 66 D. TEST CASE FOUR - DEEP SOUND CHANNEL 72 E. TEST CASE FIVE - DOUBLE CHANNEL 77 VI. CONCLUSIONS 79 A. GENERAL 79 B. FUTURE PROGRAM 82 APPENDIX A NORMOl FINITE DIFFERENCE SCHEME 85 APPENDIX B N0RM04 CONNECTION FORMULAE 87 APPENDIX C N0RM05 CONNECTION FORMULAE 90 APPENDIX D INPUT DATA 93 FLOWCHARTS 95 COMPUTER PROGRAMS 108 BIB LIOGRAPHY 156 INITIAL DISTRIBUTION LIST 158 FORM DD 1473 160 LIST OF TABLES I. Results of Test Case One 61 II. Results of Test Case One - NORMO Results 63 III. Results of Test Case Two 65 IV. Results of Test Case Three 70 V. Results of Test Case Four 74 VI. Results of Test Case Five 77 VII. CPU Times for Various Runs 80 LIST OF FIGURES 1. Typical Deep Ocean Sound Profile 21 2. Deep Ocean Velocity Function Profile 22 3. Typical Mode Solutions 26 4. Successive Modes 29 5. A WKB Profile 36 6. The Turning Point Connection 39 7. Surface Mode Function Vs. Horizontal Wave Number 45 8. Multiple Sound Channels 47 9. Test Case One Sound Profile 59 10. Test Case One Modes 60 11. Test Case Two Sound Profile 64 12. Test Case Three Sound Profile 67 13. Test Case Three Eigenvalue Errors 68 14. Test Case Three Modes 69 15. Test Case Four Sound Profile 73 16. Test Case Five Profile 76 17. Double Channel Mode Comparison 78 18. Sound Speed Profile With "Levels" 83 19. Proposed Sequence of Mode Calculations 84 TABLE OF SYMBOLS The following symbols are used within the text of this thesis, but do not necessarily apply to the computer programs. The figures in parentheses refer to the equation which defines the symbol, or in which the syrabol is first mentioned. a - coefficient for particular solution of Z, (57). Ai - Airy phase solution. b - coefficient for particular solution of Zj^ (57). Bi - Airy phase solution. c - sound velocity. c-u - bottom sound velocity. ^min ~ i^inimum sound velocity in the sound velocity profile. C - coefficient for particular solution of Z^ (24). C - group speed (40) . C - phase speed (39) . D - coefficient for particular solution of Z-^ (24). e - 2.718281828459045. E - energy level, or horizontal wave function (17). F(z) - square of the vertical wavenumber (92). h - depth grid spacing. H - depth from upper (lower) turning point to surface (bottom) (102). 1 r (1' 2) J-|.o Hankel functions of the first and second kind. k - wave number (14). ^max " naaximum possible real wave number (16). kj. - horizontal wave number (14). K-r - eigenvalue (horizontal wave number), ^m ° k„ - vertical wave number (14). Kg - bottom imaginary vertical wave nuinber (24). K - imaginary vertical wave number at water side of bottom (107). K_ - imaginary vertical wave number (49). Li - integral of imaginary wave number across a barrier (86). m - mode number. N - maximum number of discrete modes available (3 3). p - acoustic pressure. r . - horizontal range. R - horizontal displacement potential (10). Ri - reflection coefficient for a barrier (85). S(z) - an integral function of depth (34). Si(z) - vertical v/ave number integral for Zi (65). 80(2) - vertical wave number integral for Z^ (64). t - time . V(z) - velocity function (16). V^ - bottom value of the velocity function. ' '..>... .-t^ 8 z - depth. Zi^ - bottom depth. ZrpT - lower turning point depth (64). ZrpjT - upper turning point depth (66). zj - top depth of a barrier (86). Z2 - bottom depth of a barrier (86). Z - depth dependent vertical displacement potential (13). Z-i - bottom vertical displacement potential (22). Zq - vertical displacement potential at water side of bottora boundary. Z__ - eigenfunction for a particular mode, m. Zj - WKB solution in the "real" region (57). Z2 - WKB solution in the "imaginary" region (59). 0^ - mode attenuation coefficient (45). Q - bottom volume attenuation coefficient (45). AS - difference function (80). _5rt) - an undefined function of depth (34). 01. - phase at lower turning point (63). 0u - phase at upper turning point (67). 01 - phase at top of barrier (115). 6i - phase at bottom of barrier (115). ^^ - complex wave number in the bottom (43). £^T ~ complex horizontal wave number (44). V - mode normalization integral (41). f - an undefined function of depth (35). 'rf - 3. 141592653589793 p - density. p, - bottom density. p^ - ocean water density. Q. - velocity weighted mode normalization integral (42). ^ - angle made by ray with the horizontal (grazing). $ - displacement potential (2). uX - angular frequency. 10 I. INTRODUCTION Since the Second World War considerable interest has developed, in both naval and scientific communities, in the prediction of sound propagation within the ocean. For the scientist, propagation informa- tion can be a vital part in investigating the acoustic, physical and sometinaes chemical properties of the ocean. For the naval commu- nity, propagation is a vital element in the prediction of acoustic sen- sor performance. The prediction of such performance is not limited to the design and development of sensor systems, but is increasingly important in the operational employment of such sensors and selec- tion of the most fruitful tactics. Accordingly, there has been considerable development within the naval community of numerical acoustic propagation models. These models have, for the most part, been based upon the theory of ray acoustics. The techniques employed have developed considerable sophistication in order to deal with ray limitations with caustics, frequency effects, and interference. However, for prediction of prop- agation over great ranges and at low frequencies, ray techniques have two serious limitations: (1) The long ranges involved require long computation times. Each ray must be "traced" over many successive short time steps to simulate travel of the wavefront. Usually more than one hundred such rays are traced. Thus, as ranges increase, the computation time 11 increases proportionately. Such long coraputation times make the routine operational use of ray trace programs covering paths of thousands of kilometers too expensive to be practical. (2) Ray trace programs largely ignore such frequency dependent effects as diffraction and interference. The assumptions of ray a- coustics require that the wavelength of the acoustic signal be short in comparison to the depth and velocity gradient scale. As lower frequency signals are considered, the wavelength becomes an appre- ciable fraction of the depth scale, and so violates a basic validity assumption for ray acoustics. Because of these limitations there has been increased interest in normal mode theory as a basis for long range propagation models. The theory of normal mode propagation is not new and has had a num- ber of shallow water applications [Officer 1958, Bucker and Morris 1965]. A. O. Williams (1970) has outlined a method whereby the normal mode technique could be applied to deep ocean acoustics. As a part of any such technique, both the characteristic horizontal wave- nTimber (eigenvalue) and the mode shape (eigenfunction) must be solved for each mode. To do this, three basic methods have been used: (1) The wave equation has been solved in closed form [Tolstoy and Clay 1966, Williams and Home 1967]. In order to do this, the sound velocity profile must be fitted to an assumed analytic function. Such functions have taken the form of Epstein profiles [Bucker and 12 Morris 1967], profiles with a constant gradient of the reciprocal of the sound velocity or sound velocity squared [Williams and Home 1967, Tolstoy and Clay 1966]. The requirement that the sound veloc- ity profile be fitted to an arbitrary function places severe limitations on the flexibility of any model. (2) The wave equation has been vertically integrated while meet- ing appropriate boundary conditions. This method has the advantage of being adaptable to any arbitrary sound velocity profile. Kanabis (1972) and Newman and Ingenito (1972) have developed two such shal- low w^ater normal mode models. The first program presented with this thesis, NORMOl, is based upon these two shallow water models. (3) The Wentzel-Kramers -Brillouin (WKB) approximation of quantum mechanics provides an analytic solution to the wave equation. This solution is singular at depths which correspond to ray vertex depths. However, the WKB approxinaation may lend itself to solu- tion of the eigenvalue, after which the eigenfunction may be solved by an integration similar to the previous naethod. The last two pro- grams presented with this thesis, N0RM04 and N0RM05, use this method. The long term goal of this research is to develop a deep ocean normal mode acoustic propagation model of sufficient computational speed to be considered for operational use. In order to achieve this goal, a rapid method of solving for the eigenfunction must first be developed. Based upon this need the immediate objectives of this thesis are: 13 (1) To compare the results of the programs employing methods (2) and (3) above with analytic solutions, the published results of other programs, and each other. (2) Determine whether the WKB method offers promise of iin- provement in terms of time versus accuracy. (3) To outline future improvements to either method based upon the preliminary results. 14 II. THEORY A. WAVE THEORY In this section we apply a mathennatical separation of variables to the wave equation; then we impose boundary conditions on the re- sulting separated vertical equation. The results of this raathematical procedure offer us a differential equation and boundary conditions capable of numerical solution. We will close this section with a dis- cussion of some general properties of such a solution. The following discussion is based upon the treatments given by A. O. Williams (1970), and I. Tolstoy and C. S. Clay (1966). For a more complete descriptive treatment the reader is referred to the discussion given by Williams. 1. The Wave Equation The simple scalar wave equation for small amplitude waves IS Y72X _ 1 'S'^ (1) where <^ is the displacement potential. Here ^, the displacement of a fluid particle from its rest position, is represented by ^- VO, (2) 15 and the acoustic pressure, p, is given by Thus, if we can express ^ as a function of space, we can then find the resulting acoustic pressure and sound propagation loss due to spreading and refraction. If we use a cylindrical coordinate system and assume ^ is a function of range, depth, and time, we can rewrite (1) as Let us further assume that $ is separable in terms of r, z, and t, such that ^(r,z,t) -. R(T-)Z(H)e'"^ (5) where <*> is the source angular frequency, Z(z) is the vertical dis placement potential function, and R(r) is the radial displacement potential function. We immediately notice that b'§ _ _ W^ $ . (6) Using equations (5) and (6) we rewrite equation (4) as (7) 16 This expression has a range dependent term, a depth depen- dent term and a term with c, the sound velocity. If we could assume that c is independent of either depth or range, the expression is sep- arable into depth and range differential equations. The obvious choice is to assume that c varies only with depth. This assumption, termed horizontal stratification, although not true over great distance, is commonly made in oceanography and ocean acoustics. Actually, for our purposes it will be sufficient if c is horizontally stratified in the vicinity of our solution, and further if any horizontal variation of c is small with respect to the vertical variation. Making such an assump- tion, we continue by separating (7) into range and depth differential equations 2 ^^ F« " ^^ • (9) Solutions to equation (8) include the Hankel functions of first and second kind representing cylindrical waves, For distances such that kj.r is much greater than one, this function can be approximated by its asymptotic form --^i^^'^'''--''- 17 We now have an assumed time dependence and a solution to the range dependent function of SB . Thus, at long range from the source we have for $ J^'kr'^ ' (12) 2 Here kj. , an eigenvalue which appeared in the separation process as a mathematical constant, will be seen to be a horizontal wave number, We w^ill consider only positive values of k , -which represent outgoing wavefronts. It remains for us to find a method of solving for Z(z) in 2 terms of the separation parameter, k . We can rewrite (9) as f5-ff-A?)H>0. (13) This equation is the separated, space form of the wave equation, or Helmholtz equation, for the vertical wave component. The expression in parentheses represents a vertical wave number, k , and is related to the (true) wave number, k, and horizontal wave number, k , by r ■ Jk^=^^.Jkl^M^ . (14) Equation (13) with a pair of associated homogeneous boundary condi- tions forms a Sturm-Liouville problem, which is solved in terms of a set of eigenvalues and corresponding eigenfunctions. It now remains for us to define and employ boundary conditions associated with the ocean vertical sound profile. Before so doing, however, we will di- gress a moment to introduce a notational convenience. 18 Since the value of c does not vary greatly (perhaps five per- cent) throughout the depth profile, the square of the vertical wave number represents a small difference between two nearly equal num- bers Jb^^c^^Jk^, (15) In order to conserve accuracy and facilitate computational speed, a notational convention is borrowed from quantum mechanics. The 2 square of the (true) wave number, k , is subtracted from an arbitrary constant of the same order of magnitude (here the maximum possible wave number) to form a "potential function. " YCH)=i#' - 4'=J.^ - 4' . (16) The square of the horizontal wave number (our separation constant) is also subtracted from the maximum wave number to form a quantity corresponding to the "energy level" of quantum mechanics, and our new eigenvalue E-AL.--tr. (17) With this notation, the Helmholtz equation (13) becomes II -[e-Y(.)]Z=0 , (18) a form similar to Schroedinger 's equation. By this notation we hope to facilitate the adaptation of the results of quantum mechanics to the problem at hand. 19 2. Boundary Conditions We continue by developing appropriate vertical boundary conditions. Figure 1 represents a typical deep ocean sound profile. Figure 2 is the same profile represented in terms of V(z). Note that at some intermediate depth the value of sound velocity reaches a min- imum at the sound channel axis. Note also that the function V(z) forms a "potential well" with depth, representing the deep sound channel. Our domain of interest is bounded by the sea surface and a flat bottom. The air -water boundary at the surface approximates a free surface at which the stresses vanish. In a fluid this boundary corre- sponds to Z(o) = 0. (19) If the limit of the second derivative of the vertical displacement func- tion is finite at the surface, the stipulation that V(z) dis continuously approaches infinity at the surface has the effect of making Z{z) vanish at the surface. Thus, in terms of our quantum mechanics analogy, the free surface boundary condition corresponds to a perfectly rigid wall of the potential well. The ocean bottom provides a boundary more difficult to specify. The ocean bottom has a complex, often layered, horizontal- ly varying structure. In addition, little is known of the specific 20 TYPICAL DEEP OCEAN SOUND PROFILE I 0 0 ill 0 0 0 i • > I y SOUND SPEED mefers /sec I006 . j M>00- [ 3000- \ hOgO . jrooo- \ ■fcOTTOM M Figure 1 21 TYPICAL DEEP OCEAN VELOCITY FUNCTION lOOO ZOOO 3000 4000 5OO0 — 4— 0.06 __i o.oa 1— ROTTOM velocity function -2 (meter ) V(z) Figure 2 22 acoustic properties of much of the bottom required tomoclel fully the bottom boundary conditions. In the face of such complexities we will approximate the bottom with a relatively simple model. We will assume that the bottom is represented by a horao- geneous fluid of infinite depth, and constant velocity. Across a fluid- fluid interface we require that both pressure and vertical particle motion be continuous. From (3) and (6), continuity of pressure re- quires e.Zo = ^bZb (20) at the bottom. From (2) we find that continuity of vertical particle motion or displacement at the bottom requires We combine equations (20) and(21) to form a single bottom boundary condition foZo ^-t ebZb a* * (22) Since we have specified the bottom to be of constant velocity and density, a known solution to equation (18) for the bottom region is Zb- c.e^«^ ^ c.e-^«^ , (23) where Ko, a real constant, is the imaginary bottom wave number, 23 Kb =i Vb - E ■ . (24) If the amount of energy represented by the bottom solution is to be finite, the limit of the integral of Z, from the bottom boundary to infinity must in turn be finite. This condition requires that the expo nential growth term of Zi be suppressed. Thus, (23) becomes b = Cze , (25) and (22) becomes ZoTI-'W^B" ^ ♦ (26) We have now specified two boundary conditions which will enable us to find particular solutions to the vertical differential wave equation (18). Before proceeding with the technique of finding such solutions, let us discuss some of their properties. 3. The Nature of the Solution Before we find the solutions, Z(z), to the differential equa- tion (18) and boundary conditions (equations 19 and 26), it will be beneficial to consider the form and properties of the solutions w^e expect. a. In the previous section we saw^ that the solution in the bottom region is of the form of a decaying exponential (equation 25), This solution requires that the quantity Kg in equation 2 5 must be real. This requirement in turn implies that 24 E0, (29) or kj.<^ . (30) min Thus, we have an upper and lov/er bound for E and k , V^>E>0 (31) and ^ > K>^ . c . c (32) min b b. Within these limits the boundary value problem is so con- strained that it can only be solved in terms of a finite number of eigenvalues, E , and corresponding functions Z (z). This eigen- function, properly normalized, forms what is termed a normal mode. At sufficient range a vertical sound pressure profile can be approxi- mated by a linear combination of such eigenfvinctions, tn'i 25 TYPICAL MODE SOLUTION eigenlunciion t> imaginary region real region imaginary region V(z) Figure 3 26 c. At depths such that E^^ V(z) the eigenfunction Z^^(z) is oscillatory in character. We can represent such a function by where \{z) and S (z) are undefined functions of depth. We shall refer to such depths as the oscillatory or "real" region (a reference to the fact that the vertical wave number, k , is a real number). At depths such that E < V(z) the eigenfunction Z (z) becomes quasi- exponential in character, with a form we can represent by where, again, ^ (z) and S(z) are undefined functions. We shall refer to such depths as the exponential or "imaginary" region. The depth at which the sign of [e -V(z) changes is termed a turning point. At the turning point the character of the eigen- function changes from oscillatory to quasi -exponential (Fig. 3). d. Each mode, corresponds to one or more rays which traces paths w^ithin the oscillatory region between turning points, surface, or bottom boundaries. The angle the ray makes with the horizontal is defined by Cos (Z^Cz) = CCz) 4^ . (36) / CO From this we see that the mode turning point corresponds to the ray vertex depth (that depth at which a ray reaches its nnaximum depth excursion). When the oscillatory region of a mode is bounded by the 27 free surface boundary or the bottom, the corresponding ray is sur- face or bottoin reflected, respectively. e. The lower boundary placed on kj. represented by equation (28) has a more farailiar and perhaps more satisfying physical meaning. Substituting equation (36) into (28) we have Cos d > C(£) , (37) At the bottom this is the expression for the critical angle. Thus, we are limiting outselves to consider only those modes whose corres- ponding ray strikes the bottom at a grazing angle less than the criti- cal angle. Such modes are widely called "unattenuated modes." Also, there exists an infinite set of modes such that "^ T^ (38) However, because of bottom reflection loss, those modes tend to attenuate rapidly with range. For propagation problems at a con- siderable range, the "attenuated modes" contribute an insignificant amount to the acoustic pressure field, and are ignored. f. Given the oscillatoiry nature of the eigenfunction Z , we see that the eigenfunctions for each successive mode must change sign one additional time (Figure 4). Each sign change will be referred to as a mode crossing. Thus, corresponding to each eigenvalue E , there corresponds an eigenfunction Z with m-1 mode crossings. 28 SUCCESSIVE MODES MODE 1 MODE 2 MODE 3 0 crossings 2 crossings 29 4. Mode Parameters We should now introduce two sets of parameters, derived elsewhere, which are used in norraal mode calculations. a. Phase and Group Speed Phase speed describes the horizontal speed of the wave front represented by a mode. As the speed of advance of a "wave- front of constant phase, the phase speed is given by Cp = ^ • (39) Group speed is the rate of energy transport in the hori- zontal direction. Tolstoy and Clay (1966) give an expression for group speed based upon a theorem by Biot (1957). This method cal- culates group speed by Cg=i.^ (40) where V is the normalization integral. (41) and -d ccif (42) b. Bottom Decay Coefficient Kornhauser and Raney (1955) give an expression which calculates the effect of a bottom absorption on the mode amplitude. Assume the wave number in the bottom is complex and given by 30 K-% *^e (43) where /9 is a bottom attenuation coefficient. Then the horizontal wave number must also be complex and given by Ar = hv + i. c^ . (44) If & is small with respect to — , then the ratio of o(. to (^ for a given mode is approximated by % f-^|^/z'<.>d.. (45) i.- boLbom B. THE WENTZEL-KRAMERS-BRILLOUIN (WKB) APPROXIMATION The WKB approximation is a method borrowed from quantum mechanics which will enable us to approximate the eigenvalues, E , of equation (18). In the WKB approximation we assume a form for the mode solution which is valid at depths removed from the turning points. We then develop a characteristic equation for the eigenvalue, based upon the assumed solution. This section follows in many re- spects the descriptions given by Tolstoy and Clay (1966) and Lauvstad (1971). 1. The Characteristic Equation Based upon the sign of [E - V(z)] the vertical wave equation (18) applies to two domains, which we have termed the "imaginary" region and "real" region. Let us rephrase equation (18) into two separate equations for the two domains. So doing, we have 31 z; -^ )^i z. = 0 and for E > V(E), (46) Z^ -K/ 2, =0 for E: < VcE). (47) Here we have defined K- [e - V(.>] Vz and when E > V(z);, (48) Kz - [Yu) - E]' when E < \Jm. (49) Now let us assume a form for the solution to equation (46), similar to that which we have previously postulated in equation (34). The assumed solution is Zu) = 5cH) e .5(2) (50) By substituting equation (50) into (46) we have as a real part S'%) -5ci.)C5'fiE)^-J,/)=0, (51) and an imaginary part 3cz)S'^Ce) + 25'z>S'2) = 0 . (52) This equation immediately yields an expression for ^ (z), 5(H) - A/SV/ . (53) 32 Now let us simplify equation (51) with an assumption. Assume that the first term of equation (51) is negligible with respect to the others, or K 5" fz) 6 w) «1 (54) With this assumption we obtain as an expression for S(z) S'cH)" = Jk^ , or upon integrating il(E) = tj z We now have as the assumed solution Zi (z) = X^^ [a- Cos Si (H) + b Sin Si{e)J, (55) (56) (57) or, in a more convenient form, HifH) = i,-'''^[aSLn(5.(H) + Q,)] (58) Similarly we obtain as a corresponding solution for the "imaginary" region, Z^u,-E;''''[C^^'^'^De-''^'^]^ (59) where, Sp(z) is defined by S2 CH) =y Kzda: . (60) 33 Having ■written the two solutions, let us take a closer look at the assumption implied by equation (54). With equations (53) and (55) we can rewrite (54) for the "real" region as 1_ «^ » (61) and for the "imaginary" region, _1_ ir''''''\«'' (6z) This assumption requires that the order of magnitude of the vertical wave number, and thus the sound speed, vary slowly with respect to depth. This is generally the case in the ocean, where the total sound speed variation with depth is on the order of five percent. We also see that our solutions lose any validity as z approaches a turning point. Recalling the analogy of the turning point to the ray vertex depth, we see that as a wave approaches the turning point its orienta- tion becomes horizontal. Near a turning point the vertical wavenum- bers, k^ and K^, approach zero and the conditions of equations (61) and (62) are no longer satisfied. Thus the solutions represented by equations (58) and (59) become singular at the turning point location. The WKB method does not, at this point, appear satisfactory for the computation of the mode profile or eigenfvmction Z . How- ever, if the WKB approximation can be used to obtain an accurate estimate of the eigenvalues, E , a finite difference scheme can then be used to calculate the eigenfunction, Z^^. 34 Consider a typical V(z) profile, with upper and lower turning points, as well as bottom and pressure release boundaries (Figure 5). Remember that we have specified a derivative boundary condition at the bottom. This boundary condition determines the eigenfunction (except for a constant of multiplication) in the "imaginary" region ^TL'^'^b* "^^ ^^® lower turning point, ztXj, the "imaginary" region solution, Z2(z), corresponds to a particular "real" region solution, Z, (z). We can designate such a solution by the use of an initial phase angle, 9 l, such that (63) where S^u^-JK^dz , (64) BOTTOM and Sjczj = ) .K<^'^- . (65) z 2rL In a similar manner we can stipulate that in order to satisfy the free surface boundary condition, Z(0) = 0, the "real" region solution, Zj, must correspond to a particular "imaginary" region solution, Z2. Again we can represent this as a phase angle, 9g, in the "real" region solution, Zj. Thus we have required the argument of Sin (Si (z) +9j_^) at the upper turning point (or surface) to be our specified value 35 A WKB PROFILE surface boundary condition velocity function £ obottom boundary condition Figure 5 36 SiCHru) *Ql =■ J M^c^z - Qu = Qs . (66) Ztl If we let Q-(x -"^ " ®s ^® ^^^ rewrite this equation as /J^jcIe -G,-0u= \[l-\/u)]%*Q^*Qu. =ir. (67) Since this equation describes the arguments for which the sine of the integrated vertical wave number is zero (a Dirichlet boundary condi- tion), the surface boundary condition will also be satisfied if r[E-V(H)]''^dE + Gl-^Gu =itiTr W-1,2,...,N. (68) Ztl Values of E satisfying this characteristic equation will be the eigen- values, E , for which we are searching. Note that this equation is also that given by Tolstoy and Clay (1966) as the characteristic equa- tion for stratified acoustic waveguides. Now that we have a characteristic equation, some method of evaluating Q^ and 0 is required. The WKB approximation provides a method for finding the eigenvalues, E , once 0 and 0, are evalu- o o m u L ated. 2. Turning Point Connection The most vexing problem in using the WKB approximation is connecting the two solutions (58) and (59) across the turning points. A number of techniques have been developed, and here we will consider two such techniques which we will term the asymptotic and Airy phase methods . 37 a. Asymptotic Method We have two corresponding solutions represented by equations (58) and (59) which correspond to the same function on either side of a turning point. As the turning point is approached let us re- quire that the limit of the first derivative divided by the mode value be continuous. This is in fact the same as the fluid-fluid boundary condition of equation (22). For the purpose of these turning point discussions we will designate the turning point depth, z— , equal to zero and corresponding to z„ of equations (56) and (60). Further, z will be positive towards the "real" region and negative towards the "imaginary" region (Figure 6). Thus, from equation (22) we have JLim Zt = ilm Zg . (69) By taking advantage of the approximation involved in equation (61), we get by substitution for the above E-Ov ^ 2-0- Ce*^''^-De-^^'*^ (70) Taking the limit we have b.= CohCGo) = D-C. (71) It should be noted that this approach is simpler than the Airy phase technique (described next) and may not be as accurate. It is included, 38 THE TURNING POINT CONNECTION A Figure 6 39 however, since it provides an interesting comparison with the Airy phase method in the computer programs. b. The Airy Phase Method What follows is an outline of a technique developed by V. R. Lauvstad (1971) of the SACLANT ASW Research Centre. Lauvstad assumes asymptotic solutions of the same form as equations (57) and (59) in regions remote from the turning points. The gradient of the square of the vertical wave number a- cross the turning point is assumed to be linear, such that [J^J'= C.2-.... (72) This reduces the vertical Helmholtz equation to One solution to this equation is the Airy phase function H = AAc(-S(z)) - BB.(-5cE)), ("^^^ where S(z), a depth integration function, is positive in the "real" region, and negative in the "imaginary" region. The asymptotic forms of the Airy phase are then compared with our corresponding assumed solutions. We can represent the asymptotic form of the Airy pl\ase for negative real argument (the "real" region) as Zi = ± [(A-B)SinSu) ^ (A + B)Cc.sS(£)] , (75) 40 and for positive real argument as By matching the above coefficients with those in our solutions we have a set of coefficient matching equations representing the required turn- ing point connection. C=^(oL-b) D=^(a.b). (78) 41 III. THE PROGRAM METHODS Here we will discuss the basic methods used by the three pro- grams (NORMOl, N0RM04, and N0RM05) presented as a part of this thesis. This section is intended as only a brief introduction to the programs and their characteristics. A following section will give an outline of the program algorithms. A. NORMOi - A FINITE DIFFERENCE TECHNIQUE NORMOl is essentially an adaptation of two shallow water pro- grams developed by Kanabis (1972) and Newman and Ingenito (1972). For a given value of the horizontal wave number, k^^or E, the pro- gram begins with the bottom boundary condition (equation 22), and steps through a finite difference scheme to the surface. The finite difference scheme, which is derived in appendix A, is from Kanabis (1972). The finite difference scheme is a third-order Newton forward' difference scheme which iterates both the eigenfunction and its deriv- ative. Using this difference scherae, a search is made for those values of k , or E, which most closely satisfy the surface boundary condition (equation 19). Such values are the sought eigenvalues. A shortcoming of this technique is that under certain circum- stances the eigenfunction at the surface, Z (0), tends to grow expo- nentially rather than to decay towards zero. Some mention of this phenomenon is made by both Kanabis (1972) and Newman and Igenito 42 (1972). Kanabis recommends decreasing the incremental search limits for the eigenvalue (which imposes a consequential increase in computation time) to overcome the phenomenon. Newman and Igenito set the mode value equal to zero at depths above the last mode cross- ing or minimum. It happens that the circumstances for this expo- nential growth, or degenerate solution, occur frequently in deep ocean acoustic profiles. Thus, a more detailed treatment of the phenomenon is required. To understand this degenerate solution let us consider a typical deep ocean sound profile and our general forms for the eigenfunction (equations 34 and 35). Consider a mode whose eigenvalue, E , inter- sects the velocity function, V(z), at some depth, Zj'tt, below the sur- face. This is a common (in fact, the usual) occurrence for the lower modes in the deep ocean, where a deep sound channel and thermocline usually exist. As noted by Schiff (1955) and Lauvstad (1971), if the value of S(z) in equation (34) is not exactly -ff^ at the turning point, the 4 coefficients C and D of equation (35) are both non-zero. Thus, the exponential growth term has some finite value. Given both sufficient depth between the turning point and the surface, plus sufficiently large values of [V(z)-E], the growth exponential will eventually be- come dominant. Thus, the solution tends to grow exponentially as the surface is approached. This degenerate solution can complicate the search for eigen- values if it is based solely upon the value of Z(0). With the 43 degenerate solution, the value of Z(0) will oscillate between extreme- ly large positive and negative numbers, even with only small changes in the horizontal wave number or E. For the lowest modes of a deep sound channel profile the values of Z(0) can so oscillate with incre- -14 mental changes in E of only one part in 10 . Thus, a method such as regula falsi (described later) can become inappropriate since a map of Z(0) as a function of E can appear discontinuous (Figure 7). In order to circumvent this difficulty, NORMOl uses a halving procedure based upon the number of mode crossings. When the sur- face boundary condition is met, the eigenvalue E is the largest 7 ' to j-^ to possible value of E such that the eigenfunction has m-1 mode cross- ings. Based on this property, the number of mode crossings is used to converge upon the mode eigenvalue. The method of successive halving, although not sophisticated, has the versatility to deal suc- cesfully with the degenerate solution and any arbitrary profile. Once an eigenvalue is found, the function Z is then checked for a degen- erate solution near the surface. If it is present a new profile is iterated in the upper "imaginary" region and matched to the lower profile. This procedure may cause a discontinuity in the derivative of Zjj^ at the turning point; however, its ill effects would usually be less serious than those of the degenerate solution. B. N0RM04 AND N0RM05 - THE WKB APPROACH In both N0RM04 and N0RM05 the vertical wave number is inte- grated between the upper turning point (or surface) and lower turning 44 SURFACE MODE FUNCTION Vs. HORIZONTAL WAVE NUMBER i-4ax«o'« Zfe) fMOOS. VALue W SO<^FACe U+io* ^ — + IO* T r /■ A^ /f I* » n T — I 1 — I j 1 1 1 r — -to*' -I 1— I r- kr, HORimONTAU WAVE Nt»M3fiA J 1 1 T^ — I 1 r— J^ ^ »--io'^ .* X X X y » ^ Figure 7 45 point (or bottom). The programs differ only in the turning point connection formulae. Both programs integrate the vertical wave number using the trapezoidal rule over the depth grid points. The eigenvalue E is converged upon using the regula falsi method. This method makes successive estimates of the zero value of a func- tion by the linear interpolation of a negative and positive value. The method is represented by El*i = E- + (E^-E-)(-AS-) , (79) (AS^ - AS-j where the subscripts + and - correspond to the values associated with the last positive and negative values of A S; and AS, the differ- ence function, is defined by /Ztu Jk^6z ^Gu-e.-'m'^- (80) Htv If the mode is not found after a set number of iterations, the regula falsi method is dropped for a halving process. The phase effect of the upper and lower turning points, 9^ and 0, , are evaluated using the asymptotic method for N0RM04 and the Airy phase method for N0RM05. The vertical wave number is integ- rated from the uppermost turning point (or surface, if surface re- flected) to the lowest turning point (or bottom, if bottom reflected). Both programs must then provide, in addition to the values of 0^ and Qj, the effect of any intermediate "inaaginary" region or barrier. This occurs when multiple sound channels are encountered (Figure 8) 46 MULTIPLE SOLND CHANNELS 1000- 2000- 3000- Ul « o o m T o CM O y SOUND SPEED mef er$ / see Figure 8, 47 Details of the derivation of the specific formulae for 9 , 9t , and the phase effect of an intermediate barrier are given in appendices B and C for N0RM04 and N0RM05 respectively. If the mode is surface reflected or bottom reflected, both pro- grams use the same formulae for the phase effect of the reflection, 0^ and 0T . Since the free surface requires Z(0) = 0, for the surface reflected wave, we have from equations (63) and (66) L «<,dz +6^ = 9^= tnTf (81) and thus 9.-0. (82) For the bottom reflected mode, we require the conditions of a fluid fluid boundary to be met (equation 22). Substituting equations (63) and (25) into (22) we have J,,CotC8.V--^KB , (83) or where k^ is evaluated at the water side of the bottom boundary. The occurrence of multiple sound channels or ducts makes the WKB method more complicated. If the intermediate "imaginary" region or barrier is of sufficient extent and magnitude, sound 48 propagation in each channel becomes independent of the other. When this is the case, the waves in one channel undergo an effectively total reflection by the barrier. Thus, complete sound ducting occurs entirely within that channel. Lauvstad (1971) gives the ratio of the amplitudes of the reflected and incident waves in a duct as Ri= e*- - e"*- (85^ e*- +■ e"^ where L is the integral of the imaginary vertical wave number across the barrier. (86) This result can be derived using either the asymptotic or Airy phase connection formulae. As L increases, Ri rapidly approaches unity, and near total reflection is achieved by the barrier. The WKB pro- grams N0RM04 and N0RM05 evaluate L,, and if it is of sufficient magnitude, the programs treat the mode propagation in each channel independently. The programs accommodate a maximum of one main sound channel (which contains the velocity minimum), and one upper and one low^er sound channel. 49 IV. THE BASIC ALGORITHMS A. NORMOl 1 . General Description NORMOl iterates the wave equation vertically from t?ie bottom to the surface using a finite difference scheme over a depth grid of up to 501 points. The program is written in IBM FORTRAN IV with subroutines used as major functional blocks. This modular programming is designed for both ease of modification and debugging. A description of the sequence of events follows. NORMOl is basically composed of three nested loops: a frequency loop, a mode loop, and an iteration loop. The input data are read using an input subroutine. The required input data are listed in appendix D. Once the data are read and appropriate con- stants computed, the sound velocity profile is linearly interpolated for all grid points. Here the frequency loop begins, the appropriate frequency is selected, and the velocity function V(z) is computed for the grid points. Next the maximum and minimum horizontal wave numbers and maximum number of modes are computed. If the fre- quency is below cut-off (the maximum number of modes is zero), a new frequency is selected. At this point the mode loop begins, and the following steps are repeated for the first to the highest requested or possible mode. 50 First, a method of successive halving is used to find two values of the horizontal wave number k which have m-1 and m mode crossings. The desired eigenvalue must lie betv/een these two values, Then, again using a halving process, the program converges upon the eigenvalue k . The mode is considered found when one of two con- ditions (regulated by an input parameter, lEX) is satisfied; IZcoil i |Z1„„.10"'", (87) or the horizontal wave number changes between iterations by an amount such that ^. < 10 -^^^^ . (88) Once the mode is found, the eigenfunction is checked for a degenerate solution (exponential growth near the surface). If a de- generate solution is present, a new function is iterated from the sur- face to the upper turning point, then joined w^ith the lower function. The mode is then normalized to its absolute maximum value, and the profile within the homogenous bottom is computed. The computed mode is plotted using a printer plot routine. A subroutine then com- putes the normalization integral represented by equation (41), and the group speed, phase speed and mode attenuation coefficient. At this point the mode and frequency loops end. Included in NORMOl is an output subroutine which provides a punched card format listing of the eigenvalues, mode parameters. 51 sound velocity profile and naode profile. The output subroutine was written to provide local input to CALCOMP plots and propagation loss prqgrams. 2. The Subroutines The following is a brief description of the major subroutines which comprise the NORMOl program. INPUT 1 provides the input data for the program. A user supplied page heading is read on the first card, followed by the run parameters, list of frequencies desired, and a maximum of 50 depth and sound velocity pairs. If the parameter representing the number of depth and velocity pairs, NUMV, is negative, all length measure- ments read are converted from feet to meters. This subroutine also prints the input data for reference. INTRPL linearly interpolates the input depth and sound velocity pairs to compute the grid point values of sound velocity. If required, this subroutine also extrapolates the last sound velocity to the bottom assuming an isothermal, isohaline gradient (0.017sec~ ). MAXMIN computes and checks the raaximum and minimum values of the horizontal wave number, and finds the number of modes represented by the minimum horizontal wave number (the maximum number of modes). If the maximum number of modes is zero, the frequency is below cut-off and the program selects a new frequency or terminates the calculations. If the maximum possible number of modes is below that requested, this parameter is changed and noted. 52 HALF searches for two values of the horizontal wave number with m (the mode being searched) and m-1 mode crossings. These values represent lower and upper bounds on k which are then suc- cessively narrowed. At the sarae time this subroutine "maps" the values of k versus mode number in order to facilitate the search for later modes. ITERAT, the heart of the program, iterates the eigenfunc- tion from the bottom to the surface using equations (97) and (101). In addition, this subroutine counts the number of mode crossings and finds the maximum absolute value of the eigenfunction. DNORMl normalizes the eigenfunction with respect to its maximum absolute value. BOTTOM creates an exponentially decaying bottom mode profile from the bottom to an input depth, PLTMAX. MODPLT and CZPLOT are printer plot subroutines which plot the mode shape and sound velocity profile. Both these subrou- tines use a Naval Postgraduate School utility printer plot subroutine, UTPLOT. FIXUP tests for a degenerate solution near the surface. If the degenerate solution exists, this subroutine applies the finite dif- ference iteration scheme starting from the surface downward to the upper turning point. This upper raode shape is then scaled to join the original lower function at the turning point. If required, a new value for the absolute maximum of the eigenfunction is then found. 53 INTEGR computes the mode normalization integral (equa- tion 41) and the integral represented by (equation 42). The phase velocity is computed; then the two integrals are used to compute the group velocity. Finally, the mode attenuation coefficient (equation 45) is calculated. B. N0RM04 AND N0RM05 1. General Description N0RM04 and N0RM05 are identical with the exception of the turning point connection formulae (appendices B and C explain these differences). Similar to NORMOl, they consist of three nested loops; a frequency loop, a mode loop, and an iterative loop represented by the subroutine SEARCH. In addition, the input, interpolation, plot- ting, maximum /minimum and mode integration functions are per- formed by subroutines similar to those used in NORMOl. The sequence of events in the WKB method programs is similar to NORMOl. N0RM04 and N0RM05 read the input data and then interpolate the sound velocity profile over a maximum of 1001 depth grid points. If requested, a sound velocity printer plot or punched card format sound profile is provided. The frequency loop then computes the velocity function V(z) and checks for maximum horizontal wave number, minimum horizontal wave number, fre- quency cut-off, and the maximum number of modes available. At this point the mode loop begins and the search for the eigenvalue is controlled by the SEARCH subroutine. 54 After the mode is found, the mode shape is computed using a finite difference scheme. The finite difference scheme is fitted be- tween the upper and lower turning points, then a solution which decays away from the turning points is calculated. After the mode shape has been calculated and normalized, the mode is integrated to provide the normalization integral, phase and group speed, and the mode attenu- ation coefficient. Finally, if requested, a printer plot of the mode shape is provided. At this point the mode and frequency loops end. Because of their importance to these programs, the SEARCH and WKB subroutines deserve more extensive description. 2. SEARCH Subroutine The SEARCH subroutine controls the iterative search for the eigenvalues k . The subroutine first makes an initial guess of an upper and lower bound for the eigenvalue. The WKB subroutine is then called to compute the difference function given by equation (80). If the difference function for the lower bound is not positive, the bound is successively relaxed until a positive value is found. The next guess for the eigenvalue is calculated using a regula falsi zero -finding technique (given by equation 79). The re- sulting incremental change in eigenvalue and the absolute value of the difference function is checked to test whether the desired accuracy has been achieved. The WKB subroutine is then called to find the difference function for the new trial eigenvalue. If the resulting dif- ference function is negative the upper bound is replaced; if it is 55 positive the lower bound is replaced. A new guess for k is then estimated, and this iterative process continues until one of two con- ditions is met. The mode is considered found when the difference function falls below an absolute limit, |AS| < 10 "^^^ ' (89) or the value of the horizontal wave number approaches machine accuracy, < 10 -14 (90) If after a set number of iterations, the regula falsi method has failed to converge upon a solution, it is replaced by a cruder and slower (but more versatile) halving process. If more than one sound channel exists, the program must determine if the two channels are independent of one another. If the L of equation (86) is of sufficient magnitude, the subroutine SEARCH computes the eigenvalues for each channel separately. Thus, after an eigenvalue has been found, two tests are made. First, it is deter- mined -whether an upper or lower (separated) sound channel exists. Then, if an upper or lower channel exists, it is tested to determine whether it will propagate an additional mode. If such a mode can be propagated, the SEARCH subroutine computes its next mode within that channel. 56 3. WKB Subroutine The purpose of the WKB subroutine is to integrate the verti- cal wave number between the upper and lower turning points. The subroutine consists of a loop which steps through the depth grid from the surface to the bottom. In any "imaginary" region the imaginary vertical wave number K is integrated for use with the turning point connection formulae. Each time the depth loop exits an "imaginary" region for a "real" region, either the upper turning point phase 0 or the phase effect of a barrier must be calculated. If the barrier is considered to be separated (if L in equation 86 is larger than a set value) the integration is stopped with the upper channel. When the depth loop reaches the bottom, the bottom phase, Q-r , is computed for either a bottom reflected or a decaying solution as appropriate. Finally, the difference function is calculated and program control returns to the calling program. 57 V. RESULTS All three programs were run with a variety of sound velocity profiles and the calculated eigenvalues examined. The profiles select- ed include the published results of other models and profiles with known analytic solutions. In all the test runs the input variable lEX was set at 10. Thus, the iterations were halted when the absolute relative value of the eigenfunction at the surface, Z(o) , or the dif- ference function, A S , (equation 80) achieved a value less than 10 Additionally, computation was halted if the eigenvalue was incremented less than 10 during a halving process (in other words, when the solution approached machine accuracy). A. TEST CASE ONE - NEGATIVE GRADIENT This test case is one given by Kanabis (1972). The profile is characterized by a simple negative velocity gradient in a shallow water channel (Figure 9). The profile has a fast fluid bottom of the same density as the water. The test results are computed for a fre- quency of 1000 Hz. Table 1 lists the eigenvalues given by Kanabis and computed by NORMOl, N0RM04 and N0RM05 in yd . The list includes a set computed by the Kanabis program and a set computed by an unpublished program by C. L. Bartberger and L. E. Ackler. Except for the higher nnodes of N0RM05, the results from all three 5 programs agree to the fourth decimal place or one part in 10 . 58 TEST CASE ONE SOUND PROFILE '3in I . ^-: i-S3 _J BOTTOM Figure 9 59 MODE 1 TEST CASE ONE MODES MODE 2 MODE 3 1 ■■■ T r 1 50 - - - 100 - - - 150 - \ 200 r \ h» \ 250 I ■\ - ) « 300 m •o X -\ ■y C^ - 5' J j^ •" y ( \ i 350 - y^^ x^ - y' ) bollom ( \ K 400 ' ' Figure 10 60 in 0 0 h- 00 r~ vD 00 o 00 in m r-H 00 o o^ (M r^ en .— ( 00 vO -^ o OO r 1 CO r^ in ^ CO r— 1 o o 00 vO 00 r^ r- r^ r- r^ r-- r^ vO sO ^x> r-- r- r-- r- r^ r- r^ r^ r-- r-- r- CO OO CO CO CO '^ o r- 00 r-- v£) 00 o r- in o in 1 — 1 s o o^ (M r- CO 1— 1 00 vO CO 00 (M »-H 00 r^ in ^ ro 1— 1 o o^ r^ xO Ph' 00 r-- r^ r-- r-- r- t^ r^ sO vO vO O r- r^ r^ r-- r^ r-- c^ r^ t^ r-- r-- CO CO CO CO CO w O in sO sD in r- o^ r^ in CO ^ 00 z ^ o o ro r^ CO o 00 vO CO 00 r-H i-H 00 r- in ^ CO 1 — 1 o O^ r^ vO 0 P^ 00 r^ r- r^ r- t^ r- r^ vO o vO w CO 0 CO co CO CO CO CO CO CO CO CO CO < u H W W H p^ h w o Ort in vO in in r^ o^ t^ in CO '^f 00 « w o o (M r- CO o 00 vO CO 00 f—i w I— H 00 r^ in ^ CO 1—1 o o r-- ■^ H 00 r- r^ r- t^ r-- r^ r^ vO vD vO P r^ t^ r^ r~ r- t^ r- r- t^ r- r^ HU • . . . . . . . . . . «< CO CO CO CO CO CO CO CO CO CO CO ^^ w < 73 U ^ o in ^ vO o^ vO CO ro ro r^ o o^ OJ r^ CO o 00 o CO 00 I-H I-H 00 r- in '^ CO t— 1 o o r^ vO 00 r^ r- r- r- r^ r-- r-- o vO O t^ r^ r-- • r- r^ r^ r^ r- r^ t^ CO CO CO CO CO CO CO CO u I g > W Q O ^^rocor;finsor*ooo^o^^ £3 O ••-I O X 61 Table 2 gives the NORMOl, N0RM04 and N0RM05 results to the accuracy computed and relative differences between NORMOl and the other two prograras. NORMOl required 56.23 sec of central processing unit (CPU) time and 412 iterations to arrive at the solu- tions. N0RM04 and N0RM05 required 23. 50 seconds with 83 itera- tions and 22. 74 seconds with 79 iterations respectively. The CPU times cited include the time required for all operations including in- put, sound speed profile manipulation and graph printing. B. TEST CASE TWO - POSITIVE GRADIENT The second test case is a shallow water example cited by Newman and Ingenito (1972). The profile is characterized by a slow, variable positive sound speed gradient over a fast, dense bottom (Figure 11). The computation frequency was 700 Hz. It should be noted that the accuracy parameter used by Newman and Ingenito was an absolute value of less than 0.001 for the surface eigenfunction Z(0). This roughly corresponds to a value of lEX between 3 and 4 in the NORMO programs. With the exception of N0RM05, all the programs agreed to at least the third decimal place (Table 3). Considering the accu- racy limit set for the Newman and Ingenito run, this agreement is about as good as may be expected. NORMOl required 27.43 seconds and 196 iterations to complete computations, while N0RM04 and N0RM05 required 14. 73 seconds for 77 iterations and 14. 92 seconds for 79 iterations, respectively. 62 o I— I Q IT) in CO CO ^4^ (M ^H 'ri* ,-^ O vO o — ' ^ •— I 1-4 X H fc < t^ w o H O O « o o o ro r^ 00 o CO 1— ( (M LD ro iD o r- r^ vO nO CO o 00 LD ro O 00 o (T- (M r^ m ^M 00 xO ^ o ro i—H 00 r- ID -^ ro 1 — 1 o O 00 vO 00 r- vD r- r-- t^ ro ID rj ro ro CO in o CO ro ro o r- ro r- 00 o ro O . — 1 00 00 sO o r^ r^ o vO 00 o r~ ID O^ T}^ o o o (VJ r^ ro r — 1 00 ID CM 00 ra nH 00 r^ ID '^ ro •— t O o^ r- sO 00 t^ vO r- r- t^ ro o en CO w o ^ m m o < ID H ro (M >— ( 00 (M ID r^ o r- O ro 1 — 1 vD o o ^ vO ID ^ r^ 00 •^ ID ro ro r- S o o ro r- ro o 00 vO ro 00 1— 1 •—1 00 i^ in -^ ro f-H O CT^ r^ o D^ 00 r~- xO 0 r^ r-- r^ z ro ro CO u m ro ro sD in W W Q O ^ ro ^ in vD 00 CT^ o ^ P U 63 TEST CASE TWO PROFILE i.-ia i.Gl \.^z ;3 2. '51 FXiO I BjTTjM Figure 11 64 o O m (Nj un rsj (M o^ ^ LTI -^t* 00 00 ro ro CO in o ■ — 1 00 '^ un 00 00 ^ sO r— I o o 00 r^ LT) O cr- 00 00 00 00 (M (VJ (NJ ro fvj ro o 0) CO (M o H W < U H W H Ph O W ^ iTi 00 ^ •vO LO sO o ■— 1 ^ CO I— 1 xO o (NJ 00 O 00 ir> 00 ro '^ r- vO (M ^ o 1 — ) o o^ 00 r- un o o 00 00 00 00 z (N3 ro ro ro ro ro o o ID 1— f 1 — 1 in '^ r- ro a^ sD r^ ^ ro ro CO 00 in ^ •O CO rt^ r- o ro -^ 1 — 1 o o^ 00 r-- in cr- o CO 00 00 00 ro ro ro ro ro ro o 0) CO CO u CO CO ro PQ < CO 0) 6 Z W g Z ^ w p o ^ ro ro xO CO ro ■— 1 o^ o in r^ ^ ro CO 00 sO in 00 CO ^ r^ ^^ ro '^ 1— t o o^ 00 r^ in o o 00 • 00 00 00 ro ro ro ro (M '^ ro ro o o (U en CD ^1 (d ^ 0) <— 1 a •1-4 cj H o N P •i-i PU 0 u K 65 C. TEST CASE THREE - SYMMETRIC WAVE DUCT A profile was constructed corresponding to a symmetric wave duct where the sound speed is given by = JL _ a^w^fzo-z] Z /-i c c, o (115) between two depths. The solution to this profile is given by Tolstoy and Clay (1966) as "" L Co*- J The actual profile dimensions, shown in figure 12, are taken from Tolstoy and Clay. The profile was tested at 30.0 and 60.0 hertz. The computed eigenvalues agree with the anlytic solutions to one part in 10 , with the exception of the higher modes of NORMOl. Table 4 shows comparisons betw^een the analytic solution and NORMOl, N0RM04 or N0RM05 computed values. Of significant interest is the fact that the errors appeared rela- tively stable for N0RM04 and N0RM05. The errors for N0RM04 are consistently 0. 45x10 to 0.69x10" meters ~ low for 30 Hz, and 0. lOxlO"^ to 0. 12x10"^ meters " low for 60 Hz. In contrast, the errors for NORMOl increased from 0. 39x10" meters " low to 12. 5x10"^ meters "^ high for 30 Hz, and from 0. SOxlo" to 2. 8x10'^ meters ~ '■ low for 60 Hz. Figure 14 displays a graph of the program errors. Figure 13 shows the first three modes. 66 TEST CASE THREE PROFILE f>liG — I Mcj ZT 1000. 0 m 1539 m/sec m/sec Z2= 2828.8m Cy- 1539 m/sec BOTTOM Figure 12 67 TEST CASE THREE MODES hOCE i hOC[- z Figure 13 68 20t. '10L SXhr o.ao TEST CASE THREE ERRORS o o N0RM01 * ^ N0RM04 and N0RM05 60 Hz 12.51 K 10^ 10 MODE -t^O " ■»<---- •~-*<' X- ' ---'<- x-----^ .-i2.0 ^^i-o Figure 14 69 I— < n 00 '* ^ o CO o o 00 o o o 1 o o 1 o 1 o 1 o 1 o 1 o 1 — 1— I 1 o iri 00 00 00 in 1 — 1 m in 00 in ro o o oo o vD o o LO o '^ o ^ o ro ^ 't CT) 1 — 1 00 vO in (^ fO CM 0 1— i i-H z o • o h O o 00 XTi IT) o f— 1 vO o Tt< ^t< ^ tn sD sO un sO l-H ^ Q X o o O 1 o O o 1 o o 1 ^ o^ CO fvj o^ r-- ON ro 00 CVJ -H o^ r^ vD -* Tt< ro N o o I o o 00 o w w X H W in < U H w W H O H P w W o in 00 00 00 in -H m ^ o CO ^J< 00 in (M O^ v£) CO o %D ro o^ c^ f) o in O ^ o ^ o ro o CO 00 ^ vO 't CO •—1 c^ 00 o in fsj <— t o CO ro t^ nO tf (\J CO Tj* rf o ■ — 1 i—t ro (VJ ;z; o o o O u, h sO o t-~ o •^ o (VJ r-- vO l-H 1 CO CO ^ in IT) (M o^ UO o . . , , • • >— 1 o o o o o o I-H v^ 00 CO r- in 00 vO 00 o ON in ;-l 00 [N- 00 00 CO 00 tN- 00 o CO in 0) r^ 00 r- 00 in CO r— 1 o 00 ^ 00 rC o l-H h- s.n vO [N- ON 00 in o in CO 00 i-H in o CO ^ 'Ct^ ■ — 1 [N- o o CO ON ^ 00 xO CO o r~ ^ o i^ -o CO o 00 o in o in ON Tt< o^ CO ON ^ 00 vO Tf^ CO >— 1 ON 00 v£> in 00 l-H ON CO CO CO CO 00 00 00 00 >^ r^ r- vO 00 r— 1 00 l-H 00 <— 1 00 f— ( 00 1 — 1 00 f-H 00 l-H 00 l-H U •-7 00 00 00 00 in vO 00 W P a w -H CO CO 70 b 3 o in ^ 00 00 O LO o r- o r-H ro ^ CO f— I in ^ 1-H o xD ro o r^ O ro 00 — ( 00 00 in o^ fx^ 1 o •—1 ro I— I o . — 1 o u. o s —1 f-H »— 1 1— 1 1 — 1 1 — r • 1— 1 r—{ Q X 1 1 1 1 1 1 1 00 o 71 D. TEST CASE FOUR - DEEP SOUND CHANNEL Test case four is a typical deep ocean sound channel. The pro- file is characterized by a deep and sharp sound channel and by a bottom of the saine density as the water (Figure 15), With the excep 4 tion of the first raode, the eigenvalues agree within one part in 10 (Table 5). Note that for the higher frequency, 60 Hz, the agreement is somewhat improved. It is interesting that the first mode should show a larger disagreement among the programs. This effect may be a result of the combination of the sharp sound channel axis and a grid spacing of 8 meters. 72 TEST CASE FOUR PROFILE c. 4 BOTTOM Figure 15 73 RESULTS OF TEST CASE FOUR Frequency: 30 hertz NORMOl N0RM04 xlO^ MODE 4 5 6 7 8 9 10 0. 12614128 0. 12610844 567736 537387 514437 492671 472443 453745 436162 419601 404106 Frequency: 60 hertz 1 0.25255354 0.25250819 4 5 6 7 8 9 10 196659 154658 118755 089177 063654 039561 016594 196827 154783 118986 089755 039343 016618 0.24994761 0.24994890 N0RM05 DIFF xlO^ 567933 - 1. 97 538187 8.00 513923 - 5. 14 492224 - 4.47 472294 - 1.49 453659 - 0.86 436012 - 1.50 419117 - 4.84 400638 -34. 68 32.84 0.12610844 -32.84 567933 - 1.79 538187 8.00 513923 - 5. 14 492224 - 4.47 472294 - 1.49 453660 - 0.85 436018 1.44 419155 - 4.46 402169 -19.37 ■45.35 0.25250819 -45.35 1. 68 1.25 2. 31 5.78 196827 154783 118986 1.68 1.25 2.31 063571 - 0.83 2. 18 0.24 089755 5.78 063571 - 0.83 039343 - 2. 11 016618 0.24 1.29 0.24994890 1.29 973743 974255 5. 12 974255 5. 12 TABLE V 74 E. TEST CASE FIVE - DOUBLE CHANNEL Test case five is a double sound channel of nearly symnaetric dimensions (Figure 16). This profile was used to test the multi-duct properties of the N0RM04 and N0RM05 programs. Note that for 30 Hz the differences between NORMOl and N0RM04 (or N0RM05) increase by almost an order of magnitude for modes 5 and 7 (Table 6) For these modes the barrier connection formula (equation 114) is employed and the ducts are considered "connected." The mode pro- files (Figure 17) for these modes also do not completely correspond. Modes 5 and 7 appear to be "upside-down" in N0RM04, with the dominant amplitude waves in the lower channel. This is due to the fact that N0RM04 and N0RM05 iterate the vertical wave number from the surface downward. Because of this, the phase at the upper turning point of the lower channel represents the phase angle for the third and fourth mode in the lower channel. Equation (114) remains correct, although here it has been incorrectly applied. 75 TEST CASE FIVE PROFILE scijNri (.Ei-aniT"- fXiQ i.SQ BOTTOM Figure 16 76 RESULTS OF TEST CASE FIVE MODE NORMOl N0RM04 DIFF xlO^ N0RM05 Frequency: 30 hertz 1 0. 12608464 0. 12604450 2 3 4 5 6 08403 04450 0.12553486 0.12554660 52366 16376 11170 53937 36343 09716 0.12485241 0.12490245 10 75394 59913 42282 74793 59806 42565 Frequency: 60 hertz 1 0.25248206 0.25242998 48204 42998 0.25179326 0.25180252 79260 31541 30922 80252 31134 31100 7 8 9 10 0.25087788 0.25088532 84379 50385 41349 84790 49918 40770 -39.53 04450 11. 76 0. 12554660 15. 71 199. 67 -14. 54 53937 50. 04 0. 12490245 - 6. 01 - 1.07 2.83 74793 59806 42568 •52.08 0.25242998 -52.06 42998 9.26 0.25180252 9.92 - 4.07 1.78 80252 31134 31100 7.44 0.25088532 4. 11 - 4. 67 5. 79 84790 49918 40770 TABLE VI 77 DIFF xlO^ ■40. 14 0. 12604450 -40. 14 ■39.53 11.76 15. 71 36343 199.67 09716 -14.54 50.04 • 6.01 ■ 1.07 2.86 ■52/08 ■52.06 9.26 9.92 • 4.07 1.78 7.44 4. 11 ■ 4. 67 • 5. 79 O w I— 1 < Pu :^ o u w Q o :^ w <; u w CQ O O Q 'igure 17 78 VI. CONCLUSIONS A. GENERAL As seen by table 7, the WKB method has performed consistently- faster than the iterative technique. The WKB programs (N0RM04 and N0RM05) found solutions with about one -half the Central Pro- cessing Unit (CPU) time required by the finite difference program (NORMOl). In addition, one should consider the method used for in- tegration of the vertical wave number in the WKB programs. The trapezoid rule was used which required the evaluation of a square root at each grid point. Although simple, this is a time consuming task, while it is not necessarily very accurate. If a faster integra- ting procedure w^ere substituted for the trapezoid rule, the time re- quired for each iteration could be reduced considerably. Thus, in addition to showing an improvement in computation time, the WKB method also offers the potential of greater computational speed through the use of sophisticated integrators. The WKB and finite difference methods appear to have similar accuracy in comparison with other programs. The accuracy attain- able is adequate to calculate useful modes for a propagation loss pro- gram. However, the errors involved in the differences between successive eigenvalues will not permit accurate coherent phasing of the modes beyond the first or second convergence zone . Coherent phasing at extremely long ranges is, in fact, attainable only with an 79 CPU TIMES FOR VARIOUS RUNS Test Case One Test Case Two Test Case Four Test Case Five NORMOl 56. 23 sec 27.43 Test Case Three 106. 35 105. 67 114.07 N0RM04 23.50 sec 14.73 38.42 40.23 44.53 N0RM05 22. 74 sec 14. 92 39.47 40.43 47.09 TOTAL 409. 75 sec 161.41 sec 164. 65 sec TABLE VII 80 analytic profile and solution — an ideal mathematical model. In fact, when one considers the accuracy of the available sound velocity data, it is difficult to imagine any model attempting to yield accurate phasing effects at any large distance. Since the position of a target is unknov/n, in practice it is more important that the amplitude and nature of phasing effects be know, than for the sharp interference peaks to be accurately predicted. Thus, accuracies of one part in 5 10 should enable a model both to take account of the phasing effects in the first convergence zone, and to give a general description of the effects of phase interference at longer ranges. Although the WKB programs do not correctly interpret the barrier effects, the problem appears to lie in the approach taken in the pro- gramming, rather than in the mathematical and physical relationships. Within any duct, the effects of any upper and lower ducts must be treated as a turning point phase effect similar to the top and bottom boundary conditions. The modes 5 and 7 of the N0RM04 30 Hz double channel case are unwitting examples of this. In N0RM04 and N0RM05 the "connected" ducts are considered as a single wave sys- tem. In fact, with any sized barrier, no matter how small, one must consider the wave system within each duct separately. In this case the turning point phase angle between the duct and barrier is analo- gous to the impedence due to the other duct felt by the wave system in its duct. 81 B. FUTURE PROGRAM It appears that a faster WKB method normal mode program can be written which will properly interpret a barrier. In order to pro- perly consider the multiple channel cases, a mapping procedure should be used before actual mode searches commence. The sound profile would be divided into "levels" (Figure 18), found by selecting all the velocity raaxima and minima. Each "level" represents a range of the possible eigenvalues, E . Within each "level" there is a single combination of ducts, barriers and boundaries. At the start of computations the program would determine the limits of each level. Then, at the beginning of the frequency loop, the program would de- termine the maximum and minimum number of modes for each duct within each level. Thus, equipped with a map of the profile charac- teristics, the program would compute eigenvalues for each duct and level in sequence (Figure 19). After all the eigenvalues are computed and sorted (some may be out of numerical order), the eigenfunctions can be computed quickly using the finite difference procedure of NORMOl. Combined with a faster integrating method, such a pro- gram should yield fast, accurate normal modes for an arbitrary sound speed profile. 82 SOUND VELOCITY PROFILE WITH LEVELS velocity funcfion Figure 18 83 SEQUENCE OF MODE CAECUEATIONS VELOCITY FUNCTION a Maximum Modes Region 1 b Maximum Modes Region 2 e Maximum Modes Region 3 d Maximuni Modes Region 4 Figure 19 84 APPENDIX A NQRMOl FINITE DIFFERENCE SCHEME In order to iterate equation (18) through the depth grid, a finite difference scheme is used which computes values of the eigenfunction, Z(z) and its derivative, Z'(z) at each grid point. This scheme was originally developed by the Arthur D. Little Company [Report No. C- 70673, 31 January 1969]. Z(z) is expanded into a Taylor series, 2! 3? where h is the depth grid spacing. From equation (18) we have . Z"rz) = - [E-\/ui Z(z) = - Pu)H(K) , (92) where F(z) represents the square of the vertical v/ave number. Using a forward difference for Z"'(z), 2'cH) = _Z^^(z,U) - Z^fz^ (93) h Substituting equation (92), we have (94) Combining equations (91), (92) and (94), z: 3?L h J * ^ ' h 85 we have on clearing, 2te.h)ri*l^F(2.h)l --Zrz)ri-h'F(2)"] ^ hZ^. (96) Placing the above result in a vector subscript notation we have the formula used to compute the succeeding eigenfunction, We must now compute the value of the derivative of the eigenfunction Z'(z) at the next grid point. In order to do this let us express the eigenfunction at our original grid point as a backwards Taylor series of the new grid point, z+h, V 3? (Vo) Substituting equations (92) and (93), we have 2. e Rearranging, we have, ZT.*v.) - l[z«-h)(l-f Rz.h)) - Z(E)(1 * I Fez.))] ^ (100) or in vector notation, 86 APPENDIX B NQRM04 CONNECTION FORMULAE N0RM04 requires the evaluation of the phase at the upper and lower turning points in addition to the phase effect of a barrier. For the purposes of this discussion we will align the origin of the depth axis with the turning point, designating the positive direction towards the "real" region. Consider the upper turning point case (Figure 5). The free sur- face boundary condition requires that the "imaginary" region solution disappear at the surface. From equation (59) this requires that Z(-H)=K^(-H)[ce^^^-"^ - De-^^^-^'^] =0. (102) Thus, we can define C and D (except for a constant of multiplication); (103) Applying equation (71), we have f°i-v-^> ::s..H> • (104) C.--e-^''-"' D-e^^'""' or 9a = Arc ban [Txt. h (JK^ d^)] (105) 87 Now, consider the bottom boundary condition and lower turning point (Figure 5). The bottom boundary condition (equation 26) re - quires ^Kb = Ks Ce!i^"^- De"^^^-"' fb Q^^zC-^) ^ l)^-Sz.i-H) ' (106) where K is defined as the imaginary vertical wave number evaluated at the water side of the bottom interface, Ks^=[V(-H)-e]=i,;-_c^^ . CC-H) :io7) Thus, the imaginary solution can be written by defining C and D such that (108) Now, by applying the turning point connection formulae (equation 71) we have (109) or 0L = Arc tan D^_ll D^ -1 (110) 88 where D^ is defined as (111) Now, consider the situation with an upper and lower channel separated by a barrier (Figure 8). Let the z origin be at the lower connection point, z , w^ith the positive direction downward. From equations (59), (69) and (71) we have for a connection at the upper point, Zi, OLi De-^ Ce^ bi De-"- - Ce^ (112) where L is as defined in equation (86). Letting A-^ equal the numer- ator and Bi the denominator, we have, except for a constant, C=^(ai-k,)e-" (113; Applying these expressions to equation (71), we have o-i (ai + bi)e^ -^ Cai- bt^e"*- or b^ (cLt ^-bi)e'- - (ai-bOe""- (114) [ Q^ = Arcba.1 r(ta.,ei-l)e-.(tan9.-l)e-^" (tcLnQi ♦De'- - (toLnQi-l)e ^ (115) 89 APPENDIX C N0RM05 CONNECTION FORMULAE For N0RM05 the development of the upper and lower "imaginary' region particular solutions by the application of the free surface and bottom boundary conditions is identical to that for N0RM04. The differences in the two programs lie in the application of the Airy phase connection formulae vice the asymptotic connection fornraulae. To obtain the upper turning point phase, we apply equation (103) to equation (77) "^-JT^^^ ^ (116) ^ (117) or 0u - Arc tan where D, is defined as (118) Dx= Ee^^^^"''^ . (119) Note that as H approaches zero (the turning point approaches the sur- face) the value of 0 in equation (119) approaches arctan (1/3) vice the phase we required for the free surface boundary condition, (0 = 0 ). 90 To derive the lower turning point phase, we apply the particular lower "imaginary" region solution (equation 108) to the Airy phase connection formulae (equation 77), t - J, [£ (eb Ks - Co Kb) e^-^-'^Lc?bKs* e^K^) e"^^'"^] ^ ( 1 2 i ) Dividing, we have ^ED, ■^ 1 0L = Arctafl 2Da- 1 (122) where D^ is defined as in equation (111). Note that this equation also does not approach the bottom reflected boundary condition (equation 84) as H approaches zero (the lower turning point approaches the bottom). Lauvstad (1971) derived an expression for the phase coefficients of a wave entering a barrier region in terms of the phase of the wave leaving the opposite side of the barrier: (123) 75" L. _J • By rearranging, we obtain 91 which agrees with equation (114), the N0RM04 result, (124) 0L "" Arc tan (ai -^^1)6^ ^ (ai - hi) e"*" (114) 92 APPENDIX D INPUT DATA The following format is used for the input data to all three programs, NORMOl, N0RM04 and N0RM05. VARIABLE CARD ONE HED (20) CARD TWO NUMV FORMAT MEANING RANGE CARD THREE 20A4 110 NMOD no lEX no N no lOUT no IDEV* 12 H F10.3 CB FIO. 3 DRO F10.3 DRB FIO. 3 HP LOT FIO. 3 HEADING. TITLE CARD Number of Depth/Sound Speed 2 Pairs. Negative if conver- sion from feet to meters de- ■ sired. 50 Number of Modes Desired Iteration limit. 1 - 100 1 - 14 Number of grid points desired 100 -500 100-1000= Type of output desired Output device for cardfile see below Bottom Depth Bottom Sound Speed Water density- Bottom density Maximum depth to be plotted HPLOT H -N0RM04 and N0RM05 only. 93 CARD FOUR NFREQ 110 Number of frequencies desired CARD FIVE (SIX) FREQ (NFREQ) 8F10.3 Frequencies desired CARDS SIX TO 5+NUMV (SEVEN to 6+NUMV) DEP(I) CC (I) FIO. 3 FIO. 3 Depth Sound Speed 1 - 20 DEP(I)> DEP(I-i; For NORMOl, lOUT = l gives a punched deck output of the sound profile and mode parameters and profile. For N0RM04 and N0RM05, the lOUT parameter gives the following types of output; lOUT 0 1 2 3 4 5 PRINTER GRAPH PUNCHED CARDS DISK OUTPUT YES YES YES no no no no YES no no YES no no no YES YES no no 94 NORMOl FLOWCHART READ AND CHECK INPUTl DATA / <7 LINEARLY INTERPOLATE SOUND PROFILE SOUND VELOC ITY GRAPH PUNCH SOUND VELOC4 ITY PROFILE FREQUENCY LOOP BEGINS IFREQ^l, NFREQ COMPUTE VELOCITY FUNCTION below cut-off freq £ CHECK FOR MAXIMUM ANDI MINIMUM NUM BER OF MODES 95 A o cy o MODE LOOP BEGINS M = 1, NMOD > 3ZL USE METHOD OF HALVES TO GET UPPER AND LOW ER BOUND OF EIGENVALUE WITH M-1 AND M MODE CROSSINGS C JSZ ITERATION LOOP BEGINS Qim ) sz COMPUTE NEXT GUESS BY HALVING JS INCREMENTAL CHANGE < DESIRED ACCURACY? no yes v- CALL ITERAT I ^ ITERATIO^ COUNTER ^ 96 m SURFACE MODE VALUE < ACCURACY DESIRED? no ves e^SZ. MDDE yes^_,,*^ '•^^^ \^ no CROSSINGS \7 P4 O s 1^ PM O o 1-^ 2 REPLACE UPPER BOUNDS 8 XZ REPLACE LOWER BOUNDS <: T J END OF ITERATION LOOP ) 10 CHECK FOR DEGENERATE SOLUTION AND FIXUP IF REQUIRED 11 4^ ( >DDE HAS BEEN FOUNDJ- I NORMALIZE MODE TO ONE ^ m 97 fi^ (14 8 p^ •-J 0 0 g hJ 0 g P o* ^1 pj Fm :^ COMPUTE A BOTTOM SOLUTION KT GRAPH THE MODE 5Z. COMPUTE NORMALIZATION INTEGRAL, GROUP AND PHASE SPEEDS 12 C ^ es ^> PUNCH MODE OUTPUT DECK END OF MODE LOOP 15 } { END OF FREQUENCY LOOP 3 XT' ( STOP J 98 N0RM04 AND N0RM05 FLOWCHART Zi READ AND CHECK INPUT DATA 71 V 1 LINEARLY INTERPOLATE THE SOUND VELOCITY 99 (h o o l-J >< u n a- i COMPUTE VELOCITY FUNCTION V(Z) below cut-off frequency P4 O s i CHECK FOR MAXIMUM AND MINIMUM MODES MODE LOOP i MODE = 1, NMODF I CALL SEARCH (SEARCH FOR EIGENVALUE) i CALL SHAPE4 C COMPUTE THE MODE SHAPE) > (h(h INTEGRATE MODE, COMPUTE PHASE AND GROUP SPEED, AND DECAY COEF. 100 m o o o >* o s iii PRINT MODE PARAMETERS C C OUTPUT? ^>^H> PUNCH MODE SHAPE DECK 1 >DDE LOOP ENDS i FREQUENCY LOOP ENDS 717 ( STOP ) } ) 101 SEARCH SUBROUTINE FLOWCHART SET FOR SECOND DUCT no (T^ SET SOME constants' MKE FIRST GUESS no ^ rO SET UPPER LIMIT 3 14 0 0 102 32 SET LIMITS FOR UPPER WELL RELAX LOWER LINTT RE-SET LOWER LIMIT UPPER DIFFERENCE ^ " 16 COMPUTE NEXT GUESS BY REGITA FALSI 17 CALL WKB I ITERATION COUNTER ^^ d yes GO TO TOP DUCT ^ 103 16 17 24 19 "A DIRECT HIT" REPLACE UPPER LIMITS REPLACE LOWER LIMITS FIND NEXT GUESS BY HALVING I 25 1_1 THE EIGEWALUE HAS BEEN FOUND no^X^ ITERiMIONS < ITMAX? ves ^„ 26 1 27 NEXT DUCT WILL BE CENTER RJlCORD LIMITS C I i NEXT DUCT WILL BE TOP RECORD LIMITS RETURN -) C I >>> C CC^MON /DBLL/ DVZ ( 501 ) , CH, DH2, DCbSG ,DRORB , CKAP, DUMAX CC^'MQN /DdLE2/ DU2 , DCB , DC:-1 1 N , Dw , CC i , DCL S t COMMON /D6LE3/ DUZ ( 50 I ) ,OZZ ( bO i ) CCNMON /DtiLE4/ DS , CS5CC , DPV EL , DC- V EL COMMON /PARAM/ N,NP1tIEX COMMON /SOUND/ Z(5(.n,C(50J CCI^MON /^lEkO/ H,ED(20) CCMMON /SINbi/ ZM,C8,CMIN,CMAX,PLT^AX,UM,FG CCMMGN /INPUT/ NUMV,NMOD CCMMON /DMAP/ OKI(LOO) CCMMON /DKMAP/ DKM IN ,DKMAX , DKL , OKU CC^f^ON /DSQUND/ DCC(50i) CCMMON /RHC/ RO,Re COMMON /BGTTO/ DEC AY, ZB ( 20 ) ,UB ( 20) CCNMON /F^L^/ FQV(20) COMMON /ACC/ CEPS,DIFFI C C«<<<<<<<«<<««<<«<<<«<<<«<«<<<<<< DATA >>>>>»>>»» C D/^TA DP2/6.282185306D0/ DATA NPRINT/6/ C C«<«<««<<««««<<<««<««<<«<<<< PROGRAM BEGINS >>> C C C WRITE (NPRINT,16) DC 1 1=1,100 1 DKU I) =-l.CDO CALL INPUTl (IOUT,NFREC} CALL INTRPL DCMIN2=DCMIN=^DCMIN CALL CZPLQT IF (lOUT.NE.l) GO TO 2 CALL P-^GFIL (CZZjDCCNPl) 2 CONTINUE NMODP=NMGD C C . FRECUENCY LOOP EEGINS C c CO 15 IF=1,NFREQ FC=FQV( IFJ DW = 0P2^=FQ D^2=DW^DW 108 c c- c c c c c c c c c c c c DKAP=DW/DCNIN c c c c c c c c c c c COMPUTE THE V(Z) FUNCTION DC 3 CCZ=DCC(J) CCZ2=CCZ*DCZ 2 DVZ(J)=DW2=^( (CCZ-DCMIN)*(DCZ + DCMIN))/(DCZ2=*CCMIN2) \tiRllE (NPRII^T,18) NNCO=NKCDP C/iLL MAXMIN (£13) FQ MCDE LOOP BEGINS DC 12 M=1,NM0D IT=1 CALL HALF (M,£;9) ITERATICN LOOP EEGINS DKN=(DKL+DKU}*0.5D0 DINC=DABS(CKU-DKL)*C.5DC/DKN IF (DINC-OIFFI) 10,5,5 5 CALL ITERAT CDKM,MSJ IT=IT+I DEPSIL = DUMAX^'^OEPS IF (DA5S(CUZ(NPU )-DEPSIL) 11,6,6 IF (MS-M) 7,8,8 DKU=DKI^ GC TO 4 CKL=DKN GC TO 4 MODE HAS BEEN FCLND S DKN=(DKL+DKU)*0.5DC IC CALL FIXUP (DKN) 11 DKNAX=OKN DK1(H)=DKN CALL DNCRMl ( DUZ , OUMAX ,NP1 J CALL BOTTOM (DKN,M) CALL MODPLT (M, 1) WRITE (NPRINT,17i IT CALL INTEGR {DKN,M) IF (lOUT.NE.i) GO TO 12 CALL OUTPUT ( DUZ , DKK ,M ,NP1 ) 12 CONTINUE 13 DC 14 11=1,100 14 DKK II )=-l.CDO 15 CC^TINUE — NCDE LOOP ENDS FREGUENCY LCCP ENDS 109 c c« c c c STOP <<<<<<<<<<<<<<««<<<<<<<<«<<<<«<< FORMAT ST/iTEMENTS >> it FORMAT (•! K.E. EVANS NORMAL MCDE E IGE NVA LLE • , -'SCLVERt VERSION 1, REVISION 6, CATED 13 ALGUST 73') 17 FCRMAT CO NUMBER OF ITERATIONS =', 15) 18 FCRMAT CO FREQUENCY : •, F6.1, • HZ . • ) END 110 StBROUTINE INPUTl ( IOUT-,NFREQ ) IMPLICIT REAL''^8(D) c * C * THIS SUBROUTINE INPUTS THE SCUNC AND PLN DATA, C * DISPLAYS IT, THEN SETS SOME CONSTANTS FCR C * LATER USE. C * c c C<<<<<<<<«<<<<<<<<<<<<<<<<<<«<<< CONNCN BLOCKS »>>»>>>» C COMMON /OBLL/ OVZ ( 501 ) , DH ,DH2 ,OCeSC ,DRGR 6 , CKAP, DU^AX COMMON /DBLE2/ D W2 , DOB , DCMI N , DW , DC 1 ,DC i SG CCNMON /HEAD/ HEC(20J CCHMON /SINGI/ ZM , CB ,CKI N ,C^'AX , PLTNAX,UM , FC CCNKON /SOUND/ Z(50),C{50) CCNMON /INPUT/ NUNV.NMOD COMMON /PARAM/ N,NPI,IEX CCNNON /DHHH/ DH2S IX , DH2D3 , OHHL F CCMMON /RHC/ R0,k6 CCMMON /FREW/ FgV(20) CC^^ON /ACC/ CEPS.DIFFI C C«<<<<<<<<<<<<<<<<«<<<<<<<«<<<<«<<< DATA >>>>>>>>»>>>» C CATA NREAD,NPRINT/5,6/ C C<<<<<<<<<«<<<«««<<<«<««<<«<<<< PROGRAM EEGINS >>>» C C CCNVRT=l.OEO C CI ■ READ IN TFE RUN PARAMETERS C RE/SD {NREAC,4) HED WRITE (NPRINT,5) HED R£^D (NREAD,6) NUMV,NJMOD, I EX, N, I CUT N=IABS(N) IF (N.GT.500) N=50C IF (N.LT.5G) N=LOO NNCD=IAES(NMOD} IF (NMGD.GT.IOO) NM0D=10G IF (NMQC.EQ.OI NM0C=10 IF (N.EC.OJ N=500 IEX=IABS( lEX) IF (IEX.GT.I4) IEX=14 IF ( lEX.EQ.O) IEX = 7 NFl=N+l RE/^D (NREAD,7) ZM , CB, RO ,RB, PLTMAX 2M=ABS(ZM) PLTMAX=AdS(PLTMAX) RE/iD (NREAD,6) NFREC NFFEQ=IABS(NFREQ) IF (NFREC.GT.20} NFPEQ=20 RE/iD (NREA0,7) { FG V ( I ) , I = i , NFRE C ) F£=F&V( U C C— ' r ^-CONVERT FPCM FEET TC METERS C IF REQUIRED C IF (NUMV) 1,2,2 eCNVPT=0.3048 NUMV=-NUMV ZM=ZM*COMVRT Ce=CB-CCNVRT PLTMAX=PLTMAX*CONVRT y^RITE (NPRINT,8) V^RITE {NPkI(MT,9) (FGV( I) ,I = 1,NFR-EC) 111 VnRITe {^PRI^iT,lO) nnociex ksPITE (NPRINT.Il) ZN'.PLTMAX WPITE (NPRINT,L2) POtRB WRITE (NPRINT,L3i CB,N C c C READ IN /»KC DISPLAY SCUND DATA WRITE (NPRINT,14) HED DC 3 I=l,NUMV READ (NRcAC,7) Z( I ),CI CI=CI^CCNVRT Z( I)=Z( IJ^CGNVRT C(I)=CI 3 WRITE (NPRINT,15) Z(I) ,CI C c Q COiMPUTE SCNE CONSTANTS C Dh = CBLE(ZMi/DFLOAT(i\) D2^'=D6LE(■ZM) CRCRb=DeLE(RO/RB) Ch2 = 0ri-'CH CCe=DBLE(Ce) CCBSQ=CCe^CCB DI-2SIX=DH2^0.1666 66b66 6666 667D0 D(-2D3=Dh 2*0.3 3333 33 33 3333 33300 DhhLF = Dh-«0.5D0 DEPS=l.COr^^M-IEX) CIFFI = DEPS-'iE FINE GRID CONFUTATION ^** C DCyiN=DBLE(CMIN) CCMN2=LCf^.IN^-DCMIT^ DC1=0BLE(C{1) ) CCiSQ=CCl^-CCl C2^.=0BLE(ZM) 1 = 1 CCI=DC1 113 DCIP1=CELE(C(2) ) C2I=DBLE(Z(U J C2IPl=CeLE(Z(2) ) D2DIFF=DZIF1-DZI CFLN=DFLGAT(N ) CZF=C.OCO C c C«<««<«««<<«««<<«<< FINE GRID LOOP BEGINS C CC 8 J=l,NPi C22( JJ=CZP IF (OZP-DZIPL) 7,7,6 C C *** RESET THE I PARAMETERS *'** C 6 1=1+1 D2I=DZIFI CCI=DCIPI CCIPI = CELE(C( I+U ) CZIPI=DELE(Z(H-l) ) C2CiFF=DZIPl-DZI C C «** COMPUTE THE VALUES FOR C(Z) *** C 7 DC2=(DCI^{DZIP1-DZP)+DCIP1*(CZP-CZI ))/DZDIFF CCC(J)=DCZ c C *** COMPUTE THE NEXT DEPTH *** C C2P = DZM='^DFLnAT(J iVDFLN 8 CONTINUE C C<<<<<<««<<<<<<««<<<<<<<< FINE GRID LOOP ENDS C C RETURN ENC 114 SLBHOUTINE MAXMIN (*) IMPLICIT REAL*8(D) C c * C * THIS SUBROUTINE CHECKS THE NAXIMUM ANC MINIMUM C * VALUES OF THE HCRIZCiMTAL WAVE NUMBER; ThEN C * DETERMINES THE NUMBER OF MODES REPRESENTED C * BY EACH. IF NC MODES ARE PRESENT FOR THE C * MINIMUM WAVE NUMBER THE FRECLENCY IS EELCW C * CUT-CFF. IF THE NLMdER OF MGCES REQUESTED IS C * ABOVE THAT REPRESENTED BY TFE MINIMUM r^AVE C * NUMBER, THE REQUEST IS MODIFIED ACCORDINGLY. C * c C C<<<<<<<«<<<<<<<<<«<«<«<<<<«<«<<< COMMON BLCCKS >>>»> C C CCMMON /DBLl/ DVZ ( 501 ) ,DH ,0H2 , DCBSC ,DRORB » CKAP , DUMAX CCMMON /DBLE2/ DW2 t CCB , CCM I N t OW , CC 1,DC LS«Q COMMON /DBLE3/ DU Z (501 ), DZZ ( 50 1 ) CCMMON /PARAM/ N»NP1,IEX CCMMON /HEAD/ HED(20) CCMMON /SINGl/ ZM , C3 ,CM IN ,CMAX » PLTNAX ,UM , FG CCMMON /CMAP/ DKL(IOO) CCMMON /DKMAP/ DKM IN ,OKMAX » CKL , CKU CCMMON /INPUT/ NUMV,NMGD CCMMON /SOUND/ Z(50),C{50) C C<<<<<<<<<<<<<<<<<<«<<<<<<<<«<<<«<<< DATA >>>>>>>>>>>>>» C DATA NREAD,NPRINT/5,6/ DATA MAXI,MINI/«MAXI« , 'MINI'/ C C<<<<<<<<<<<<<<«<<«<<<<<«<«««<<<< FROGRAN BEGINS »>» C C C **« COMPUTE THE MAXIMUM AND MINIMUM DK'S *** C CKMN = OW/DCB DKMAX = 0^^/DCMIN C C *** CHECK OUT THE OKMIN FOR MAX NO. OF MCCES *** C CALL ITERAT (DKMIN,MrS) IF (MS) 1,4,1 1 MCCMAX=MS DK1(MS)*DKMIN C MCCMIN=0 MS = 0 C C *** PRINT THE RESULTS *** C WRITE (NPRINT,5) HED WRITE (NPRINT,6) WRITE (NPRINT,7) M AXI , CKM AX , MOD MN WRITE (NPRINT,7) f^ INI , OKMI N , MOD MAX C C *** CHECK THE NUMBER CF MCDES PECUESTEC *'^* C IF (MCCNAX-NMOD) 2,3,3 2 NNCD=MOOMAX WRITE (NPRINT,3) WRITE (N PR I NT, 9) WRITE (NPRINT,10) WRITE (NPRINT,9) WRITE (NPRINT,8) WRITE (NPRINT,il) NMOD 115 2 RETURN *=*»? 6ELGW CUT-OFF FREQUENCY **'5= UPITE WRITE WRITE WRITE WRITE WRITE RETUkN (NPRINT,8) (NPRINT, 9) (NPRINT, LO) (NPRINT,9) (NPRINT,8) {NFRINT,12) 1 C c C<<<<<<<<<<«<<<<««<<<<<<<<<«<«<<<< FORMAT ST/iTEMENTS » C c c c 8 c 10 11 12 5 FCR 6 FCR 1 7 FCR 1 FCR FCR FCR FCR FOP NAT MAT MAT (•IS CN HORIZONTAL* , F12.8, 3X, 20A4i (•0', ///, • WAVE CO', T22t •MODES :' (T33,13( '* (T33, 1H-,11X,1H*) (T33t'=?= CAUTION *') (•O't T26, 'MODES RECUESTEC RESET TO : CO', /, T25, 'FREQUENCY RECUESTEC IS CUT-OFF', //, T20, 'JOB TERMINATED T25, 'BOUNDS NUMBER' ) A4, 'NUM K : }) ENC NAT MAT MAT MAT MAT • BELOW •NEW FREQUENCY SELECTED.') 14) OR 116 SLEROUTINE HALF (M,*) I^FLICIT REAL-8(C) C Q *:»:k3|t^:>;j?3h*«j;i5k**3{c**«3|. ****♦*******:«!************** ^: ******* C ^ C * THIS SUBRCUTINE SEARCHES FOR AN UPPER VALUE C * OF CK WITH M-l KODE CROSSINGS A,\D A LG^nER C * VALUE CF DK WITH M f^GUE CROSSINGS. AT THE C * SAME TIME IT FILLS IN A MAP CF DK VS. C * MODE CROSSINGS FOR LATER USE. C * Q 5{=***********'':**^*** *>'**** A **********■*-***** ************ C C«<<<<<<<<<<<«<<<<<<<<<«<<<<<>>>>>>>>>>>> CCNMCN /DBLE2/ DW2 » CCB t COM IN,DW ,CC I ,DC LSC CCNFON /PARAM/ N,NPI,IEX CCMNON /DMAP/ DKI(ICO) CCMMGN /DKMAP/ DKM IN, DKMAX , CKL , DKU CCM^ON /ACC/ DEPStCIFFI C C<<<<<<<<<<<<<<<<««<<<<<<<<<«<< PROGRAM BEGINS >»»»>>> C LCW=-l LUf=-l C C *'?'* CHECK THE MAP FOR ALREADY COMPUTED GLESES *** C c c c c c c DC 2 I=Mt 100 IF (DKl(I) ) 1,1,3 1 CONTINUE 2 CONTINUE DKL=DKMIN GC TO 4 3 DKL=DK1( 1) 4 CKU=DKMAX C C C *** FIND £ CHECK THE NEW MID DK VALUE *** C 5 Ci>>>» C C CCNMCN /DELI/ DVZ ( 501 ) t CH, DH2 , DCBSG ,DRORB , DKAP , DUMAX CCKMON /DBLE2/ 0W2 , CC B , COM 1 1\ » DW ,CC L »DC1S G CCMMON /DBLE3/ DUZ ( 50 I ) ,OZZ ( 50i ) CCNFCN /PARAM/ N,NF1,IEX CCF.MON /SINGl/ ZM,CB,C>aN,Cf^AX,PLTNAX,UM,FG CCMf/ON /CKf^AP/ DKMIi\|,CKMAXt CKL,CKL CC/^MON /OHHH/ DH2S IX t DHzD3 » CHHL F C c C«<<<<<<<<«<<<<<<<<<<<<«<<<<<<««<< DATA >>>>>>>>>>>>>>> C DATA CUPPER/C .1D6C/ DATA NPPIiNT/6/ C C<<<<<<<<«<<<<<«<<<<<<<<<<<<«<«<<<< PROGRAM EEGINS >>>» C c C . DEFINE SCNE CONSTANTS C DSTART=i.CC-10 CK2=DK-CK CE=(DKAP-OK)*{DKAP+CK) C C- r — ' SET UP RUN PARAMETERS C 1 CLVAX=OSTART CN>lRK=i.OD0 NS = 0 C C ' INITIAL VALUES C CUJ=DSTART CLFJ=DkCR6«DSCRT ( ( CK-CKMIN )* ( DK■^DK^'IN) )*DSTART CFJ=D£-CVZ(1) CLZ(1)=CUJ c C — ' ^ -— ^ — DEPTH LOOP BEGINS C c DC 8 J=2,NPI C C ITERATE THE FUNCTION C CFJPL=LE-DVZ{J) CENCM=l.0DC+DH2SIX*CFJPI DLJP1 = ( {I.0D0-DHzD3^CFj J>^:=CUJ + Dh=?CUPJ)/DENCN D^^( l.GGC-DH2C3^DFJPl )^:=CUPJ CE=(i:FJ-fDFJPl-DH2S IX^ C F J->>>» C CC^^'CN /DBLl/ CVZ(501 » ,CH,CI-2TDCeSCiDR0RBtCKAP,DUMAX CCMMON /DBLE2/ 0'/^2 , DCB . DCM I N , D W tCCi ,0C1 SG CCf^MON /DELE3/ DUZ ( 501 ) , DZZ ( 501 ) CC^.XON /SINGl/ Z.M,CB,CMIN,CNAX,PLTNAX,UM,FC CCKMON /BOTTO/ DECAY, ZB ( 20 J ,UB ( 2C ) C C«<<<<<<<<«<<<<««<<<<<<<<<<<<<«<<< FRQGRA^' EEGINS >>>>> C Q COMPUTE CEPTh VALUES C C HE=(ZM-PLTyAX J*0.5E-i DhB=DBL£(F6) C C J COMPUTE SOME CONSTANTS C CGAM=DSQRT( ( DK-DW/DCE ) '!' ( DK + DW/DCB ) ) DARG = OGAM--CHB IF (DtXFU+OARG) 3,3,1 1 CEX=CEXP(DARGi CV = DUZ( 1)-CKCR6 C C : COMPUTE THE PROFILE C c CC 2 1=1,20 IF (DV.LT.l. 000-30) GC TO 4 DV-DV'^DEX 2 UE(I)=SNGL{DV ) C RETURN C C- . . I^UL L PR 0 F I L E C 3 1=1 C 4 CC 5 J=I,2C 5 UE(J)=0.0 C RETURN ENC 124 SLERCUTINE INTEGR (CK,NCDE) IMPLICIT REAL^^8(D) c * C * THIS SUBROUTINE COf^PLTES THE NC RI^AL I Z/^T I CN C * INTEGRAL, AND THE SCLND VELOCITY NORMAL- C * IZATICN INTEGRAL. THE.M THE GRCUP AND PHASE C * VELOCITIES ARE CCMPUTED. C * C C C«<<<<<<<<<<<<<<<<«<<<<<<<<<<<<<<<<<< COMMON BLOCKS >>>>» c c CCN^.ON /PARAM/ N,NPI,IEX CCNMCN /SINGI/ ZM , CB ,CM I N , CNAX , PLTNAX ,UM , FC CCMMON /DEL!/ OV^Z ( 50 I ), LH , DH2 , DCdSG tORORE , CKAP , DUMAX CCNKON /DBLEZ/ DW2 , CCB , LCMl N , DW , CC 1, DCl S<* CCNMON /D6LE3/ DUZ ( 5C1 ) , DZ Z (50 1 ) CCK>!ON /DELE4/ DS t CSSCC , DP V EL , DGVEL CCf^CN /PHC/ RO,RB CC^'^'ON /OSCUND/ DCC(501) CCNNCN /BCTTO/ DEC AY, Zb ( 20 ) t U8 ( 2C) cc^^'ON /cknap/ dk^^in,ckmax,ckl,cku ECLIVALENCE (DUZ 1 , CLZ ( L ) ) c C<<<<<<<«<<<«<<<<<«<<<<<<<<<<<<<<<<<< DATA >>>>>>>>>>>>>> C CnA NREAC,NPRINT/5,6/ C C<<<<<<<<<<«<<<«<«<<<<<<<<<<<<«<<<< PRCGRA^< BEGINS >>>» C c C *** SET UP THE END POINTS FOR THE TRAPAZCIDAL RULE C DLZlSC^CUZd ) IF (DABS(DLZLSQ) .LT. l.CC-30) DUZlSC=O.ODG CLZlSQ=CUZiS(.«DUZlSC CS = C.5DC-(DUZ1S>^) DCISQ=CCC( 1)*DCC( I ) CLZ1SG=CLZ1SQ/(CC1SC) CSSCC=G.5D0«(DUZ1SG) C Q INTEGRATING LOOP BEGINS C C CC 1 1=2, N CLZI2=CLZ( I) IF (DAES(DUZI2) .LT.i.CD-30) GO TC 1 -. CLZI2 = CUZI2--i'CUZI2 CS=DS+0LZI2 CCCI=DCC(I) ClZI2=CyZI2/(CCCI*CCCII DSSCC=DSSCC+DUZI2 I CONTINUE C c C INTEGRATING LOOP ENDS C c I C *** ADD BOTTOM EXPONENTIAL DECAY *** C CA2 = (DK-CKMN)*(CK + CKMN) CE = CbLE(R6)=^DUZl-C.50C/DSQRT(DA2) CS=CELE(RCi*OH*DS+Ce CSSCC = DSSCC=^DBL£( RC}*DH + DB/CCBS C CECAY = CkN*CB/(CCB»=DK*DSJ C c C **♦ CONPUTE THE PHASE AND GROUP VELOCITIES *** C 125 c c CFVEL=Cl\/DK OGVEL=CS/(CFVEL^DSSCC) k^PITE (NPRINT,2) ^'GCE,DPVEL WRITE (NPRINT,3) NCDE,DGVEL RETURN C c C<<<<<<<<<<«<<«<<<<<<<<<<<<<<<<<«<<< FORMAT ST/iTEMENTS » C C 2 FCP^^AT COPHASE VELCCITY FOR MOCES 14, • IS •, i G15.7t • f'/SEC) 2 FCRr'iAT COGRGUP VELCCITY FOR MOCES 14, • IS •, 1 G15.7, • N/SEC) ENC 126 c c c c c c c c c c c c c c c« c SLEROUTINE MODPLT (M,JUMP) INFLiClT RbAL^8(D) * THIS SlibRQUTIKE PLCTS THE MCCE SHAPE CN THE PRINTER USING THE N.P.S. ROUTINE UTPLOT . IF JUMP=1, ALL POINTS IN THE MODE FUNCTION WILL EE PLCTTEu. IF JUNF IS NOT EQUAL TO I, EVERY NPl/IOO TH POINT WILL BE PLOTTED. <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<«<<< CCMi^iON BLCCKS >>>>>» CCNMCN /DBLl/ DVZCSOU ,CH,DH2,DCESC,CR0RB ,CKAP2,DUNAX CCNMGN /DbLE2/ DV^2 , CCB, CCM IN, DVv , CCl tOC iSw CCr^r'CN /D3LE3/ DUZ ( 5Ci ) i DZZ (bO i ) CCKNON /PARAM/ N,NF1,I£X CCNNON /HEAD/ HEC(20) CCNMON /SINGl/ ZM,Ce,CMIN,Cf^AX, PLTNAX,UM,FG CCN'MON /DMAP/ DKl(lOO) CCf^NON /FLCT6/ RANGEl^) CCNMON /INPUT/ NUNV,NMGD CCN.NON /BGTTO/ DEC AY, 26 ( 20 ) ,U& { 2C ) C €<<<<<<<<«<<<<<<<<«<<<<<<<<«<<<<«< PROGRAM BEGINS >>»>> C C C C c c c c c c c c c c c c c c c c CK=DK1(N) *** CHECK IF FULL PLOT OF ALL FCINTS DESIRED =*** IF (JUMP-1) 1,2,1 1 JL^P=NP1/100 2 JLf'P2 = JU^P*2 NLNB=NP1/JUMP *^* PRINT THE HEALINGS ** l^PITE (6,0) HED,Fe WRITE (6,7) N,D(<1(^) WRITE (6,8) *^* CONFUTE THE RANGE FACTORS FCP PLOT SIZE ^^'^ PANGEin = l.l RANGE(2)=-1.1 R/!NGE(3)=ZN **^ WILL BOTTOM PROFILE BE INCLICEC? IF (ZM-PLTMAX) 4,3,3 2 NCCCUR=0 RANGE(4)=G.G GC TO 5 **=? DRAW BOTTOM FBCFILE *** 4 VCCCUR=1 RANG£(4)=ZN-PLTMAX CALL UTFLCT (UB , ZE , 20 .RANGE , 1 , MCCCUR ) MCDCUR=3 **♦ DRAW OCEAN U(Z) PROFILE ^^^ ** 127 c c 5 CALL UTPLGT ( CUZ , CZZ , NUMB , RANGE , JUMP2,M0CCLR) FETURN FCPMAT FCR^'AT FCRMAT I • ENC (• I', CO', ( 'O' , VS DEPTH') 20A4, ' FREQUENCY : T25, '^^GDE', 14, ' //, T23, 'NORMALIZEC , F6. I, • HZ' } K =', F16.12) EIGENFUNCTICN ' 128 c c SLBROUTINE CUTPUT (CUZ , CK, MODE, NFTS ) Ih'PLICIT REAL-*8{D) CINENSICN CZZ(NPTS), DUZCNPTS), CCCINPTS) CCNf^ON /PHC/ RO,Ra CC/^MON /SINGI/ ZM,C8,CMIN,CPAXtPLTNAX,UM.»FC CCFN.GN /BGTTG/ DEC AY, ZB ( 20 ) ,UB ( 20 ) CCNMON /DBLE4/ D S ,D SSCC ,OPVEL , DGVEL CCNMON /HEAD/ HED(20) cc^^'0^ /wcrk/ extra(20) DATA NPRINT,NPUNCH/6,7/ 1 yvRITE (NPUNCH.IO) MODE ,FQ ,DK,DPVEL , DGVEL , CS , DECAY DC 2 1=1,20 2 EXTRA( I )=UB(21-n c c c c c c c WRITE (NPUNCH,L1} EXTRA IvRITE (KPUNCHtLl) CUZ WRITE (NPRINT,13) MGDE RETURN ENTRY PROFIL(CZZ,CCC,NPTS I WRITE (NP0NCH,7) HED NPTS20=NPTS+20 WRITE (NPUNCH,8) NFTS20, RO, BE, CM IN t CB, ZM , PLTNAX DC 3 1=1, 20 EXTRACI )=ZN-Ze(21-I) WRITE (NPUNCH,12) EXTRA C2M=DBL£{Zf^} DC 4 I=l,IMPTS CZZC I)=D2K'-CZZ{I ) WRITE (NPUNCH,12) DZZ DC 5 1=1, 2G EXTRACI )=Cfi WRITE (NPUNCH,12) EXTRA WRITE (NPUNCH,12) CCC WRITE (NPRINT,SJ C c DC 6 I=1,NFTS 6 DZZ( I}=DZM-DZZ( I) RETURN 7 FCRMAT (20A4) e FCRNAT (I10,2F10.5,4F10.2) 9 FCRMAT COSCUNO ViLCCITY DETAILEC PROFILE PUNCHED*) END 129 SLtRCUTINE CZFLQT IMPLICIT REAL=*=8(D) C C * THIS "SUBRCUTINE PLCTS THE SCLND VELCCITV C * PROFILE CN THE PRINTER PLOT USING C * The N.P.S. ROUTINE UTPLCT. C * C C<<<<<<<<<<«<<<<<<<<<<<<<<<<<<<<<<<<<<< COMMON 6LCCKS >>>» C CCyNiO.M /DBLI/ DVZ( 501) ,DH,DH2»DCBSC ,DR0RB,CK/iF2,DUNAX CCNMON /DeLE2/ D W2 , CCB , DCM IN , DW , DC 1 ,DC LS G CC^MOl\ /DBLE3/ DUZ ( 5C1 ) ,DZZ (50 I ) CCNMGN /PiiRAM/ N,NP1,IEX cc^^;oN /fead/ hec(20) CCf^MON /UCRK/ WRKL(20) CCNMON /SINGl/ ZM , CB , CM IN , CMAX, PLTMAX,UM , FG CCNf^CN /PLCT6/ ftANGE(4) CCMI^GN /DSCUND/ DCC(501) CCFNQN /6GTT0/ DEC AY, ZB ( 20 ) ,UB ( 2C ) C C«<<<<<<<<<<<<<<<«<<<<<<<<<<<<<«<<< PROGRAM EEGINS >>>>» C C *** V^RITE THE HEADINGS *** C C c c t^RITE (e,5) nED V^RITE (^»6) WFITE (6,7) IF (PLTMAX.LT.ZM) PLTMAX=ZM Fe=(ZM-FLTMAXi"0.5C-l C C -*- DEPTH INCREMENT LCCP BEGINS C DC 1 1 = 1, 2C WFKKI )=CB ZE(n=HE^-^fLOAT(I-l) 1 CONTINUE C Q ^ ,^ DEPTH LCCF ENDS C C *** COMPUTE THE RANGE FACTORS *=** C RANGE { 1 ) = AMAXUCB,CMAX} R^NG£(2 )=CMIN R/SNGE(3)=2M- -C C *** WILL BOTTOM PROFILE BE DRAWN? -** IF (ZM-PLTMAX) 3,2,2 NCDCUR=0 RANGE(4)=0.0 GC TO 4 3 MCCCUR=1 RANGE(4 J=ZN-PLTMAX C C =*** DRAW THE BOTTOM C(Z) PROFILE *** C CALL UTPLOT ( WRKl , ZE ,20 ,RANGE , 1 , MCDCUR I MCCCUR=3 C C *** DRAW THE OCEAN CCZ) PROFILE *** C C ^ CALL UTPLCT ( CCC , CZZ ,NP 1 , RANGE , 2 ,MCCCUR ) C RETURN 130 5 FCPNAT (MS 20A4 ) 6 FCKMAT ( 'O* , // , T28, 7 FCRMAT CO', //, Til, 1 •NETERS AECVE BOTTCN ENC 'FLOT CF C(Z) VS CEFTHM •BOTTOM AT ZERO. DEPTH IN • C(Z) IN KETERS/SEC ) I^PLICIT REAL*8 (A-H, V-Z ) , L0GICAL*4 ($) C c c c C«<<<<<<<<««<<<<«<<<<<<««>>>>>>>>>>>>» C CCNNCN CCNMON CCN^'lON CC^'MON CCNMON cc^MON / CI / / 02 / / D3 / / 04 / / C5 / / 06 / CCN^ON / 07 / CCMMGN CCf^MOiN CCM^ON CCi^tVON cc^^'QN/ CCNI^ON CCNiVON / OS / / 010/ / OK / /AWKB/ CKi-'AP/ /HEAD/ / 11 / CCMiMON / -12 / CCr'MON CCMMQN CCMMGN CC^iMO.M 13 10 LI PI VELG ZZZ( CC(5 ORn» CM AX H, CH DH2, OK,C DIES ZiMAX woce F(,V( OKI { H£D( NUMV NP2i LOMO JHI , JUPP NREA $GPA DPI , (102 1G21 Oi,D ORB, ,CMI ,HPL H80T Kf^IN T,DS ,css 2C) , Cl!\T 100) 20) tNMG ,NPT DE, I JLCW ER, J D,NP FH,i DPIC 1), DEPTH (1021) ), DVZ(lOCl) EP(50) DRORB N,CB,CBSQ GT,0hhLF,DH2SIX,DH2C3, T,DHB ,DKMAX,DKOLC,CKUf DKL, MAX tDE tCECAY,DPVEL,DGVEL,G/iy NP1,WQCMIN,CVZB FCtW,^\2 ,GIFF,OS,DF 0,NFREC,N^^CCF,N♦NP1» S,NCDE,IT ,NF2 HNGDE ,Kl^ELL,IWELL,NWELL. BCTT,ISHAF,JTPL, JTPU RINT,NPUNCh eCTPR, JCAPCS,$CEL 2,DPIC4,DPiRT2,D2PI C C ' C«<<<<<<<<<<<<<«<«<<<<<<<< PROGRAM EEGINS >>>>>>>»>>>>» C C WRITE (KPRINT,6) CALL INFUT2 CALL LINEAR IF ($GRAPH) CALL CZPLCT IF ($CEL.OR.$CARDSi CALL PRCFIL ( C EPTH , VELC , NPTS ) C c- C DC 5 IFREG=ltNFREw FREQUENCY LCCP BEGINS C c- c FC = FQV( IFREQ) W = FC>"-D2FI wcce=w/CB ^ccMIN = ^^/c^> fc FCR^AT CI NORMQA :_WKB_METHCD NOR^'AL MODE' 133 ELCCK C INPLICI CCMVUN CCNMQN CC^MON CCMMOiM I CCNMON CCMMO.N D/TA CATa LATA CATA D^iTA C^TA C4TA C^TA ENC 134 SLEROUTINE MAXMI3 (*) INFLICIT Rf:/iL*8 (A-H, V-Z ) , LGGICAL*4 ($) LAST CHANGED 3 AUG 73 C C c c c €<<<<<<<<<<<<<<<<<«<<<<<<<<<<<<<<<<<<< COMMON BLCCKS >>>>>> C CCNMON / C5 / CMAX,CHlN,C6,C.iSQ CCFy.ON / D? / DK,CK(^1K,0KMAX fDKCLCCKU, CKLt 2 DTEST,DSKAXiDE CCNNON / DK / F4.V{20) ,FC,U,tv2 CCf-'MON /AWKB/ DKN , D I NT ,C I FF , OS , CN CCf^MCN /HEAD/ HtC(20) CCMMON / II / NUMVTNMOC,NFkEC,NMCCF,N,NPl, 3 NP21,r.PTS,MGDE,IT,NF2 CCf^^ON / 10 / NREAC,i\PRL\T»NPUNCF CCN'MON / LL / $GRAPr,$BCTPR,iCAHCS,$CEL CCMiMON i ?l ^ DPI,DPI02,DPIC4,DPIf>>>> C c^=o.oco DKMN=^/CB CKNAX = W/CiXIN CKN = OK-MIN CALL WKfc CSMAX=DS/DPI NCENAX= ICIiNT(DSMAX} CKCLD=OKMAX WRITE (NPftINT,3) HED,FQ CK^ = OK.^:AX NCDMIN=C ksFITE (NPRira,4) CATA MAXI/'MAXI' C C- c c c c- /, MIM/'MIMV MAXI tOKMAXf f'GC^'l^ WRITE (NPi^INT,5) M IN I , DKM IN , MODM AX WRITE (NPRINT,?) •CHECK NUMBER MODES PEGU5STED NNCCF = i\NCC iF (MODNAX-KMQD) 1»2,2 1 CCNTINUE ■IF WE GET FERE WE MLST RESET THE NUMBER CF MODES REQUESTED KNCOF=MCCMAX WRITE WRITE WRITE WRITE WRITE WRITE (NPRINT,6) (NFRINT,7) (KPRINTtS) (NPRINT,?) ( NPRINTt6) (NPRINT, 9) NMCDF 2 IF (MOCNAX.NE.O) RETURN •IF WE GET CUT - OFF FERE WE ARE BELOW FREQUENCY WRITE (NPRINT, 6) WRITE (NPRINT,?) WRITE (NPRINT, 8) ^RITE (NPRINT, 7) WRITE (NPRINT, 6) WRITE (NPRINT, 10) RETURN 1 135 c 6 7 £ 9 IC FCRMAT FCFMAK FCRi^AT FCFMAT FCh^'AT FCRMAT FCRMAT FCBMAT 1 • 2 'IjENCV END •C (T33, (T33, {T33, CO', I HCR12CNTAL ^.AVE NUMBER*) ,14 F12.8,3> ,'NCDES CIS lOAe, 4X, •FRECUENCY : ', Fd.l) ,// ,T25, 'BCUNDS ON •t T22, A4, »MUM K i3( •*' ) ) ih-i, llXtlH*) «-^ CALTIGN •■*•.«) T26, 'MJUES REQUESTEC RESET (•O* ,/ ,T25, CUT-OFF' t//TT20i SELECTED' ) 'FRECUENCY RECUESTEO IS JOB TERMINATED OR NEW TC : ', I EELQW FREQ' , 4) 136 SLERCUTINE IMEGR INPLICIT REALMS (A-H, V-Z), L0GIC/^L-«*4 ($) C C LAST CHANGED 9 AUG 73 C --^^ CCNyJN / Dl / VELD(1021), DEPTH (1021) CCNMQN / C2 / ZZZ(1C21), DVZ(lUOl) CCPHON / D4 / DROtCRbjDkORB CC^MCN / C5 / CMAXtCiy IN,CE3,CeS& CCf^MOK / D6 / H,CH,hFLCT,DHhLF,Dh2SIX,DH2C3, 1 DHPHPnTTOHP. CCN^ON / C7 / DK,CKMINtCKiXAX,CKGLC,CKU,DKL, 2 OTEST,DSMAX,DE CCMMQN / DS / ZMAX,CSS,CeCAY»DPVEL,DGVEL,GA^ CCN/^ON / DIO/ k^OCE,l>LCNFl tWCCr'IN^CVZb CC^^'GN / Gift / FCV (20) ,FG,W»l<2 CCMMON / II / NUiXV,NMaDTNFR£C,Nf^CCF,N,.MPl , 3 NP21, NPTS, MODE, IT ,NP2 C C c c c c DSS=C.OCO C.SSCC = O.0D0 c c C <<<<<«<<<<<<<<<<<<<<<<<<<<<< INTEGRATING LOOP BEGINS C C CC 1 1=2, N C ZZZI2=ZZZ( I) IF (DAES(ZZZI2) .LT.i.OD-30) GO TO 1 ZZZI2=ZZZ12^ZZZI2 CSS=DSS+ZZZI2 VELCI = VELC(I } 2ZZI2=ZZZI2/(V£Lai*VEL0I) CSSCC=DSSCC+ZZZi2 1 CCNTINUE C C <<<««<<««<«<<<«<<<<«<«< INTEGRATING LCCP ENDS C C IF (DABS(ZZZ(NP1 ) ) .LT.l .OC-30) GC TO 3 A2=DSQRT(CVZ3-DE) e = ZZZ(NPl)vZZZ(NPiJ*DRaR6--^DRQ/( 2.0DG=^A2) 2 CSS = CRO'=Oh -CSS+B CSSCC=DSSCC^^DRO-OH + B/CBSQ • CECAY = W*=b/(CB*OK^DSS} C CPVEL=W/OK CGVEL=OSS/(DPVEL^CSSCC) RETORN 3 e=C.ODO GC TO 2 ENC 137 c c SLBROUTINE WK3 IMPLICIT RHAL-8 ( A-H , V-2 ) , LCCi C ^ L'-i-H {$) LAST CMANG50 : 7 SEPT 73 ASYMP7CTIC B . C C c C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<«<< CCMMON ELCCKS >>»>> C CCKMQN CC/^VUN / D2 / / D4 / / D6 / CCMMON / D7 / CCMMON CCM.^ ON CC^MON / DIG/ /AWKB/ / II / CC^MGN / 12 / CCMMON / 13 / CCMMON / LL / CCMMON / Pi / CCMMON / DTH/ CATA DEXPU /I ZZZ ORG rii C UH2 DK, DTE woe OKN NUM NP2 LOM JHi JUP $GR DPI DTH .702 2(1021), DVZdOOi) H,H ,He DKM ST, ,DI V,f\ liN ODE tJL pr-R APh ,0F ETL PLC OTT IN, OSM CCN :JT, MOD PTS , IH ow , ja , $3 IC2 kORB T,OHhLF,DH2SIX,DH2C3, HR ,DH3 OK.MAX,OKQLC ,DKU,OKL AX,DE ,FI,'aOCMIN,CVZB DIFF ,C5,DM i,NFf<£C;,NMCCF,N,-MPl, ., MODE, IT ,NF2 IMODEjKwELL tlVvELL ,NWELL, CTT,ISHAP,JT?L, JTFL OTPK,$CARCS,$CEL ,DPIC4,DFIFT2,D2PI DWELL /6.0/ C C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<«<<<<<< PROGRAM EEGINS >»» C CINT=O.CDO JTFL=I CE=(DKMAX-OKN)*(DKMAX+CKN) CKB=DSORT{ DVZ&-Ot) $TCP=. FALSE. IVvELL=I ISIG=-i CAPG=O.CCO CARGP=C.CDC C c c C c C C C c c c c C c c c •DEPTH LGCP BEGINS DC 21 J=1,NP1 CKZ2=DE-DVZ( J) •CHECK FCR CONDITIONS AND TYPE CF SCLUTIGN If (DKZ2) 1,2,3 1 IF (ISIG) 4,5,6 2 IF (ISIG) 7,21,9 3 IF (ISIG) 10,19,20 IMAGINARY REGION CAPGP=DSQRT(-DKZ2) C4PG=CARG+CARGP GC TO 21 ISIG=-1 GC TO 4 C/!RGP = DSCRT(-DKZ2) DX=CKZ/ (CKZ+DARGP) CINT=DINT.+ CKZ='(DX-i.OD0)-0.5DO^CF ISIG=-1 JTFL=J-1 CARG=( 2.0DC-DX)'^DARGP*C.5D0 GC TO 21 ZERO VERTICAL WAVE NUMBER 138 c c c c c c c c c c 10 u 12 12 1^ it 1? 18 19 2C /Lc 23 2A ie/5CK=-l ISIG=0 GC TO ^1 IE/^CK=I GC TO 8 •REAL REGION IS IF JT CK CX IF CT CI ST GC CT IF IF IV. JT CT CI GC e = CE Cc f^C EC DI IF Ci GC CI CI C/i GC IG=l ($TC PC = J Z = DSU = CARG KG=DA ($Tu HETU = NT = DT CP=.T TO 1 1 = DM0 (DAR (IWE PG=C. ELL=I Pb=J hETU = NT = DT TO L C3I.\( ccas( L = OEX fa = DE = CEL^'= = DEL- NTP=D (DIN NTP = D TO L NT = DI KT = OI RG=0. TO 2 F) GC TO 11 RT( F/( kG-» P) CAT FET PUE e C(D G=^D LL. ceo UEL DKZ2i LKZ+OARGF) DARGP^^(DX-I.CC0)*C.5DC GO TO 13 AN{DTANF(CARG*CH) ) U INT,D2PI ) h-CWELL ) EG.KWELL) GO L + 1 TO ZQ CP FE £ DT CT P( XP (A (A in TP IN fc ^T NT CC 1 104 TU 11 1 ) DA (- + e + 6 T- .G TP RG^DH) DARG --Dh) ) + Ccf^L+(A-B) DT1+DATAN2( AC,BO) E.DINT-CTtST) GO TO +DPI02 17 + {2.CO0-CX)-DKZ'^C.5D0*=DH 0 CKZ=DSQRT(CKZ2) ISIG=+i CINT=DINT+DKZ^DH IF ($TCF) GC TO 13 JTPC=J GC TO 12 CKZ=DSCRT(DKZ2) CINT=DINT+OKZ'»DH 21 CONTINUE IF (DKZ2) 22,27,25 CEPTH LCCF ENCS •DECAYING SCLLTICN CKS=DARGP eARG=(DA^G-C.5D0'i«CARGP )=feDH*2.0D0 IF (CAKC-OEXPUi 23,24,24 C^=(CKS-i')R0R6"-0KB)*CEXP(-DARG}/ ( CKS + DRORE*C KB ) CTFETL=CATAN((1.0CC+D2)/tl.GD0-D2)) GC TO 27 CThETL=CPI04 GC TO 27 139 c c JTFL=NP1 IF (DKB.tC.C.ODO) GO TO 26 CThETL = CATAN(DKZ/(CKB'f=CRGRB)) GC TO 21 BOTTOM REFLECTED 2e CTI-ETL = CFIC2 ■COMPUTE CIFFERENCE 21 C£ = CINT+QTI-ETL CINT=DINT-DTHETU CIFF = DS-DM-*DPI N^;ELL = IWELL FETURN 6 ItNELL=IV^ELL+l •GC TO 24 ENC 140 SLEROUTINE SEARCH I^FLICIT i^E^L ^8 (A-H, V-Z), L0GICAL*4 1$) LAST CHANGE : 15 AUGUST 1973 C c c c C<<<<<<<<«<<<<<<<<«<<<<<<<<<<<<<<<<<< COMMON BLCCKS >>»» C CCMMON / C7 / CCNi'^ON / CC^NGN/CK CC^MCN /A CCMMON / D6 / iVAP/ WKB/ I 1 / CCNI/QN / 12 / CCNMQN CC^MON CC/^HON cc^^>>» C C C C C c c c c FIX=1.400 IF (imol;e.gt.u go TG 1 KWELL=1 $LCh=. FALSE. $hl=. FALSE. CKCLD=DKMAX JLPPER= 1 Jhl = l JLCW=NP1 IFMCDE=C LCNCD£=C CKCEN=CKGLD CKC£N2=CKCEN jeCTT=NPl INITIAL SETTINGS ■SET UP FOR DM=DFLGAT(MODE-LaMODE-IHMGDE) IT = 0 IF (SLOVn} GC TO 11 CKGUES=(CKCLD-OK.M IN) /(DSMAX-DM+ 1.000) CKL=CKCLC-CKGUES*FIX IF (OKL.LT.DKMIN) CKL=CKMIN CJ GO TO 2 7 CKCEN2=CKCEN CKCEN=CKK JFJ=JTFL JLCW=JTPL IF (JTPU.LT.JUPPER) IF (KWELL.LT.NWELL) KtsELL=I IF (JTPL.GT.JEOTT) LONCCE=0 ISFAP=I RETURN •CENTER WELL IHf^CCE = 0 GC TO 28 •LPPER WELL CI'CLC = CKCEN JLPPER=JHI $hl=. FALSE. K^^ELL = 2 ISFAP=2 RETURN •LOWER WELL CI>>>» C CCf^i^GN CCKMCN CCNNON CCNKON / D2 / / C4 / / C5 / / ce / CCi^MGN / D7 / CC^.MON CCNMQN CCf^MON CCNMGN CCNMQN / C9 / / 010/ / DW / /AkvKB/ / 11 / CCNNCN / 12 / CCNFGN CCNMON CCNKON / 13 / / IG / / LI / / OTH/ ZZZ ORG CMA H, C DH2 OK, DTE ZMA UGC FQV DKN NUf NP2 LGM JHI JUP NRt $GR OTH ( 1C2U tCRE.D XtCf^lN H,hPLC ,HBGTT CKMN , ST,DSM X ,CSS, B ,hCCi\ (20, F ,CINT, V,N^GD 1,NPTS CCE, Ih , JLQW PER, JB AD,NPR APh, $6 ETU , DVZ(LOCl) RORb ,CB,CESQ T,DHHLF,DH2SIX,0H2Ca , ,0H3 DKNAX,DKCLi:,CKU,CKL, AX,DE C£CAY,CPVEL,DGVEL,G/iN . PI ,WQCN'JK,CVZB Q , W , W 2 CIFF,CS, Ct* ,NFR£C,NNCCF,N,NP1, ,MUCE, IT ,KP2 NODE, KW ELL, I WELL, NVv ELL, GTT, ISHAP,JTPL, JTPL IiNT,i\FUNCF GTPR,$CARDS,$CEL C C«<<<<<<<<<<<<<<<<<<<<<<<<<<«<<<<<<<< DATA >>>>>»>>>>»» C C;iTA CUPPER /9.99C 60/ <<<<<<<<<<<<«««<<<<<<<<<<«<<<<<< PROGRAN BEGINS >>>>> C c<< c CSTAR CE = (D IF (I IF (I $GC = . ISTAR ISTCP GC TO ITCP= GC TO ITCP= GC TO ITCF = GC TO ISTCP GC TO ISTCP CCNTI AVKZ = ZJST = F/iCTG ZPJST ZJ = ZJ 2FJ = Z T = K,v SF SF FA T = = N ( 1 4 JL 4 JU ( = J 7 = J NU DI DS R = = A ST PJ l.CCO /iX-CK) -(DKN-^X + CK) AP..EQ.3.A^C. IhiVCCe.GE.l ) ISFAP = 6 /^P.. EQ.l.ANC.IFMGDE.GE.l) ISFAP=5 LSE. JTPL PI 1,1,1,2,3,3], ISFAP Cl^ + 1 FFER 7,5 ,6,7,7,6) , IShAP FI-1 bCTT-1 E NT/(CFLQAT{ JTPL-JTPU )*Ch) Ii\(DThETU)^'^CSTART ZJSJ^^iUll ISTART)-DVZ( 1ST ART- 1 ) ) / ( A VKZ*4. DO=»'DH) VKZ*DCCS (DTh ETU )-^OST ART -FACTOR ST C c c c FJ=DE-DVZ( ISTARTJ JSTART=ISTART JSTCF=ISTCP CEK0M=i.0DC-DH2SIX*{DE-DVZ( JSTART-l) ) ZJ;a=( ( 1.0D0-DH2D3^FJ)>^ZJ-DH^-ZFJ)/CEN0M ZNAX=G.CD.O CC 10 J=JSTART, JSTCP ZZZ( J)=ZJ FJ=DE-DVZ{ J) ZJF1=(2.0D0-DF2*FJ)«ZJ-ZJMI ZjM = ZJ •DEPTH LOOP BEGINS 145 c c c c c c c c c c c /S2J = CABS(ZJ) IF (J-JTPL) 9,9,8 e IF (AZJ.LT.ZFiAX*CTEST) GO TC 11 IF (DAbS(ZJPl) .GT.AZ-J) GQ TC 11 S IF (AZJ.GT.DUPPER) Gu TC 29 IF (AZJ .GT.ZMAX) Z^AX^AZJ IC ZJ=ZJP1 IF (ISTCP.EC.NPl) GO TO 16 J5TART=JSTCP+1 GC TO U 11 JSTART=(J+jTPL)/2 12 ZJ=ZZZ( JSTART) CC 15 J=JSTART, JSTCP 2ZZ( v)=ZJ CKZ = DSWf>>>>> CCN/^GN / 04 / DRQ,CRB,DR0R3 CCN^ON / C5 / CMAX,CMI,\i,CBtCBSG CC^NGN / 06 / H,CH,HPLCT ,DHHLF,DH2SIXiDH2C3, 1 DH2,HP0TT,DHB cc^^'.c^l / o? / CK,CKMK,L;K;^Ax,OKCLC,CKU,DKLf 2 DTEST,OS/-*AX,0E CCKMCN / 09 / ZMAX,CSS,LtCAY,DPVEL,DGVEL,G/i^ CCNf^ON / DU / FQVdC) ,FC,WtW2 CCNMGN /AWK6/ DKN , 0 INT , C IFF ,DS , DN cc^^»» C C /^CDE OUTPLT SEGTICN C 1 WRITE (^PU^CH,6) NGCE,FG,CKN,CPV£L,DGVEL, CSS, DECAY WRITE (NPUNCH,7) DLZ C IF ( .NGT.SGRAPH) GC TO 2 WRITE (NPRlNTjS) f'CCE RETURN C 2 WRITE (NPRINTilOJ RETURN C Q SCUND PRCFILE OUTPUT C E^TPY PFCFIL (DZZ, CCC, NPTS) WRITE (NPUNGH,3) HEC WRITE (NFUNCh,^) NPTS , ORG , CRB, GN I^ , CB, H, HPLCT WRITE (NFLNCH,6> CZZ WRITE (NFUNCH.B) CCC WRITE (^FRI^T,5) RETURN C c C«<<<<<<««<<<<«<«<<<<<<<<<<<<«<<< FORMAT STATEMENTS » C 3 FCR^AT (10A8) ^ FCRMAT niCt2Fl0.5,4FL0.2) 5 FCPI^AT (•CSCUND VELOCITY DETAILEt PROFILE PINCHED') 6 FCRMaT (I5,F5.0,G2G.12,2F10.2,2G15.8) 7 FCRNAT (1CF8.5) 8 FCFNAT (10F8.1) S FCRMAT CO GLTPUT FOR f'ODE*, 14, • HAS BEEN PUNCHEC) IC FCRMAT (•+•, T118, 'OLTPLT PUNCHEC*) ENC 143 SLERGUTINE INPUT2 1^FLICIT REALMS (^-H,V-Z}, LOGICAL-4 ($) C C :- C L/5ST Ch/iNGElD 3 AUGbST 1973 Q c C<<<<<<<<<<<<<<<<<<«<<<<<<<<<<<<<<<<<< CCMMON BLCCKS >>»>> C CCNFGN / C3 / CC(5C) tCEP(50} CC^^'C^ / D4 / DRC tCSBfChOKb CCNMON' / C5 / Ci^'.A>(,C^IN,C6,GBSQ CCKNON / Dfc / H,Ch,HFLCT,DhhLF,Ch2SlXTDh2C3, 1 DH2,hEC7TtOHB CC^MG'N / D7 / DK,CKyiN,CKMAX,OKCLC,CKUfOKLt 2 DTtST,OSMAX,De CC^^f^.ON / DU / FQV(2C),FC,l-<,^2 CCM^CN /HEAD/ HEC(20) CCMMQN / IL / NUMV,NM0D,NFREG,NMCDF,N,NP1, 3 NP21,NPTStMGDE, iT,l\F2 CCNMCN / IC / NRtAC,NPRINT,NPUMCF CCNMCN / LI / $GRAPH,$BGTPk»3.GARCS,$CEL C C<<<<<<<<«<<<<<«<«<<<<<<<<««<«<<< CATA >>>>>>>>>>>>>» C CATA CC^VRT / 1 . OOOCCOOCCOCOOGOO CO/ c C«<<<<<<««<<<<<<«<<<<<<<<<<<<<«<<< PROGRAM BEGINS >>>» C C C -— ' ^ READ IN THE HEADING C AND RUN FARAMcTERS C RE^D (NREAC,9) HED ksRITE (NPRINT,13J FED READ (NREAD,10} NUNV, NiVCD , I EX ,N , ICLT , IDEV C C . . '^ TEST THE RUN PARAMETERS C IF {N.GT.IOOG) N=1G00 IF (NUMV.GT.5C) ^L^'\/ = 50 IF (NMOC.GT.LOO} NNCD= 100 IEX=IABS(IEX) IF (IEX.GT.14) IEX=14 DTEST=LC.ODO-^v(-lEX} IF (IQGT.LE.2) SGPA FH= .TRUE . IF ({ IOLT.Ew.2) .OR. (ICUT.EQ.3)) $CEL=.TRLE. IF (( lOUT.EQ. i) .OR. ( IGGT.EQ.4) ) $CAkDS= . TR L E . IF (lOEV.NE.O) NPLNGh=iCEV NF1=N+1 NF21=N+21 BEAD (NREADtll) H , CB , DPC ,DRB ,HPLCT C C ■■ READ IN THE FREQUENCIES C RE/5C (NREADtLO) NFPEw IF (NFREG .GT.20) NFREG = 20 READ (NRcAOtLl) ( Fu V ( I ) ,1 =1 ,NFR EG) FG = FQV( 1) C C . CCNVERT AND DISPLAY C THE RUN FARAyETERS C IF (NUMV) 1,2,2 CCNVRT=C.3C48 NLNV=-i\LNV h = F-'CCNVRT CE=Cd^CCnVST FPLCT = hPLCT='CONVRT WRITE (NFRINT,14) l^RITE (KPRINT,15) ( F3V ( I) , I = 1 , NFR EG) l^FITE {NPRINT,16) NNGD,IEX 149 c c c c c URITE (NPRINT,17) WFITE (^PHI^T,la} ^^RITL- (^PRIl\TtL9) IF (IGRAPh) WRITE IF (SCARCS) ^;PITE IF (iCEL) WRITE (NPRINT,20) IDEV WRITE (NPRINT,23) HED H,HPLCT ce ,N (NPRINT,22) (i\FPlNTt21J •READ IN THE DEPTHS AhD SOUND VELOCITIES C c c c CC 2 RE/C (NRtAClU CEFI,CCI D£PI = DEFr-CCNVRT CC I = CCI^CONVRT CC (I)=CCI CEP{ n^DEPI WRITE {NPKINT,12) DEPI,CCI CCNTINUE CF = h/DFLCAT(f>l ) DPCR6=LhC/CRE CF2 = DH-'CH CeSC = CB''CB Dh2S IX = Ch 2 -^0.1666 0666 6 6£66 7 CF2D2=DH2SIX^2.CDC CFFLF = t(-^^ 0.500 IF (DEPI-H) 4,6»5 NLNV=f\U^V+l CC(NUM\/)-{h-0EPI ) ^'0.016 500 + 001 l;ep(nomv)=h COMPUTE SCNE CONSTANTS 6 NFTS=NP2L IF IHPLCT-h) 7,24,8 7 hPLOT=H 24 5eCTPR= .FALSE. NFTS=NP1 RETURN 8 HECTT=hPLGT-H CFe=HBGTT^'0.5D-l RETURN C C C«<<<<<<<<<<<<<«<<<<<<<<<<<<<<<<«<<< FORMAT STATEMENTS >> C 9 IC LI 12 13 14 15 16 17 F I 18 F 1 19 F 20 21 22 23 1 CPNAT CRMAT CRMAT CR^AT CRMAT CRMAT CRMAT CRMAT T24 CRMAT T22 CRMAT • METE CRMAT 'NUiMBE CRMAT CRMAT CRMAT CRMAT •DE NC (10A3 (5II0 (dFlO 0-, 1', 0«, 0', u', ITE C , .«MA 0«- RS/SEC ( R ( ( ( ( PTH 0', ,// OS 0', ) ,12) .3) T25 lOA T34 T2 3 T24 RATI T27 X DE T18 • ,/ T22 ,T42 T42 T42 T42 LOa METE , Fc.l, T57, F6.1J 8) , 'RU , 'FR , "MO ON LI , 'dC PTh p •BC ,Ti7 •CU •GRC •GU ' FU •SO // FS* N PAR ECU EN DES R MIT : TTCM LGTT5 TTCM , • ^ u M TPLTS UP AN TPUT NCFEC Lf^D V , T33 T44, AMETE CIES EgUES OEPTH D :• , SOUND BER C ftcOU D PHA ON FI CARC ELGCI , 'Vc •SCUN 1, , //, MM • VETE : • .(T42,F6.1, • HZ* ,/) ) TEC :• , 14 ( -M 12, :•, F3. 1, f8.l, • NETEFS* ) VELGCITY :«,F8 F GRID POINTS ;' ESTED : HCRIZ WAVE SE VELCCIT lES' ,// J LE FT' ,12, 'FCCl') CUTPUT' ) TV AND MCCE GRAPHS ' ) LCCITY PkCFILE' ,// ,T22, C VELOCITY, K/SEC, /) RS',//, 1, ,15) 150 SLBROUTIKE MODPLT INPLICIT REAL-^a (A-H, V-Z ) , LOGICAL-^4 ($) CINENSICN RANGE (4) C c => C * THIS SUSROUTIKE PLCTS THE MCCE SHAPE CN THE PRINTE C * UTILIZING THE STANCARD (AT NPS) ROUTIiNE UTPLCT C '^ c * c c c CCNNON / Ci / VeLC(1021J, DEPTH (1021) CCf^MON / D2 / ZZZ(iC21), DVZ(L0CIJ CCM^uN / C5 / CMAX,CMlN,CB,CbSg CC^'^'0^ / D6 / H,Ch ,HFLCT ,DHHLF,CH2SIX»DH2C3, 1 0H2?HB0TT,DHB CCNMON / D7 / DK,i:Ky,IN,CKi^AX,DKCLC,CKUf DKL, 2 DTEST,DS^^AX,DE cc^MQN / ca / diffi CC^"^'ON / D9 / ZMAX, CSS, DECAY, DPVELjDGVEL, CAN CCI^NON/DKMAP/ DKUICO) CCNMGiN / OW / FQV(20) ,FQ,W,W2 CC/^MON /AWKB/ DKN' , D INT , CI FF , CS , CN CCFNON /HEAD/ HE0(2C) CCNMQN / II / NUMV,NNCC,NFREC;,NNCCF,N,NP1, 3 NP21,NPTS jMGDE,IT.NF2 CCMMON / 10 / NREAD,NPRfNT,NPUNCH WRITE (NPRINT,!) HEC,FC WRITE (NPRINT,2) f^.GCE,DK WRITE (NFRINT,3) CALL UTPLGT ( ZZZ , DEPTH , NPTS , RANGE ,2 , 0 , PBCTT ) WRITE (i\PRINT,4) CPVEL , CSS , DEC AY , CGVEL , I T , CIFF I RETURN ENTRY C2FLCT FECTT=SNGL(H) WRITE (KPRINT,5} HEC WRITE (NPRINT,6) WRITE (NPRINT,7) RANGE( 1 J^SKGL(Co) RAKGE(2) = SNGL(CiMIN) RANGE(4 ) = 0 .0 RANGE(3)=SKGL(HPLCT) CALL UTPLCT ( VE LC , CEPTH ,NPT S , RANGE ,2 ,0 , PE CTT } RANGE( 1 )=i.O RANGE(2)=-I.O RETURN 1 FCPMAT (M', 10A8, • FREQUENCY : «, F6 . I , • HZ') 2 FCRMAT CO', T25, V^QDE', 14, • : K =S Fid. 12) 3 FCPNAT CO', // , T2j, 'NORMALIZED t IGENFUNCT ICN ' , I • VS CEPTH «) 4 FCRMAT COPHASE VELGCITY :', F7.1, ' M/SEC. DSS ^ I l'GI57 1 4X, ' BOTTOM LEAKAGE CQEF :', G15.7, //, $ ' GROUP VELOCITY :', 2 F7.1, ' ^'/SEC. ITERATIONS :', 15, ICX, 3 'DIFFERENCE FUNCTION :', GU.E ) 5 FCPFAT CI', 10A8) e FORMAT CC', //, T2^, 'GRAPH OF SCLND VELOCITY', I ' VS DEPTH' ) 7 FCFNAT CO', TI9, 'CEPTH IN f^ETERS.', 1 ' SOUND VELOCITY IN .METERS/SEC) ENC 151 SLBROUTIKE LINEAR IMPLICIT REAL*8 (A-h, V-Z ) , L0G1CAL''=4 ($) C C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<«<<< CCMMON eiCCKS >>>>» c cc^^c^ CCMFON CCf^NON ccpyoN CCVMON CI D2 D3 C5 D6 CW II CCNMON / LL / VELCdOZl), CEPTH (1021) ZZZ( 1C2I) , OVZCIOCI) CC(5C),D£P(50 ) CMA)(,Cf'IN,CB,CeS(. H,Dh,hPLCT,UHHLF,0H2SIX,DH2C2, Dh2,hLGTT,CHB ZMAX,CSS,CECAY,CPVEL,CGVEL,GA^ FQV( 20) ,FQ, W, V^2 !\UMV,NMCD,NFRE(.»NNCCF,N,NP1, l\P21 ,KPTS,MCDE,ITTfNF2 $GRAPF,$bGTPk,$CARLS,iCEL C C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<«<<< FRQGRAI- BEGINS >»» C CNAX=C.CCQ CMN = GB CL=C.OCO CL=DEP{2) VL = CG( I ) VL=CG(2) CL^CU=CL-DU VLNVU=VL-VU 1=2 C c c c c c c c c c c c c CC ^ J = I , i^J P L CEFJ=Dh=^-CFLCAT(J-l ) DEPTH( J)=DEPJ IF (DEPJ-DL) 2,2, 1 1 l=I+i CL = DL VL = VL CL = CEP( I ) VL = GG( I ) CLNDU = L>L-DU VLNVU=VL-VU 2 VELOJ=(DEPJ-DU)*VL^'VL/DLMDU+VU IF (VELCJ.GT.GNAX) CMAX=V£LCJ IF (VELCJ-CMIN) 3,4,4 ■CEPTH LCCP BEGINS 3 CMN = VELCJ A VELG(J)=VELOJ IF ( .NCT.SBCTFR) GC TO 6 NF2=NP1+1 •CEPTH LCCF ENDS CC 5 J=l,2C VELO{NPI+J)=GB CEPTH(NF1.+ J) = F+DHB*CFLGAT( J-1) 5 CCNTINUE ASSIGN BCTTCM VELCCITY VALUES RETGRN ENC 152 SLfcRCUTINE WKB IMPLICIT P.EAL*8(A-h,V-Z),LQGICAL^4(S) C C : C L^ST CH/INGED 2 7 SEPT 73 AIRY PHASE B.C. C _ ^ C C<<<<<<<<<<<<<<<<<<«<<<<<<<<<<<<<«<<< COFMQN BLCCKS >>>>>> C CCr-FCN / C2 / ZZ2(iC2U, CVZ(iOOi) CCf^MCN / D4 / DROtCRB.DRORB CCNNGfN / L6 / H, CH ,HPLCT » DHHLF, Dl-2SiX, DH2C 2 , 1 DH2 ,hGGTT ,DHB CCMMON / D7 / DK,DKNIN,DK.MAX,OKCLC,DKU,DKL, 2 DTEST,CSMAX,CE CCMMON / 010/ waCB,y^CCNPl,WCCf>'IN,CVZB CCMiMQN /AWKB/ DK N , D Ii\T , 0 I FF , CS , CN CCf'f^ON / U / NUFV,N^GC,NFREQ,Nf/CCF,N,NPl, 3 NP21,KPTS,MGDE,n »KP2 CCMMON / 12 / LGKCCE, lHNOOE,KWELL,Ii/vtLL ,N^ELL, ^ JHltJLGW CGKMCN / 13 / JUPPER, JBGTTt liHAP tJTPLtJTPL CCf^NON / Li / SGRAPF, $6GTPR,fCAHCS,SX£L CCI^NON / PI / DPI,CFIC2,DPIC4,DFIRT2,D2PI CCMMON / DTH/ DTHETL CATA DEX^U/l.7L2/,CWELL/6.0/ C C«<<<<<<«<<<<<<<<<<<<<<<<<<<<<<<«<<< PROGRAM BEGINS »>» C CIM=O.CDC $TCP=. FALSE. JTFL=1 CE={DKMAX-DKN)-(DKNAX+DKN) CKe=DS(.HT{CVZB-CE) IWbLL=l ISIG=-l CARG=C.CCO CAPGP=C.COC C C r CEPTH LCCP BEGINS C c CC 21 J=1,NP1 CKZ2=DE-DVZ{ J) C C . . CFECK FCR CONDITIONS C AND TYPE CF SOLUTION C IF (DKZ2) 1,2,3 1 IF (ISIG) 4,5,6 2 IF (ISIG) 7,21,9 3 IF (ISIG) 10,19,20 C C . IMAGINARY REGION C 4 CARGP=DSCRT(-DKZ2) CAPG=CA8G.+ CARGP GC TO 21 5 ISIG=-1 153 GC TO ^ C c c c c c IC 11 13 1^ 15 16 17 IS 19 2C C C C C CyiRGP = CSGRT(-CK22J CX=DKZ/ (CKZ+DAPGP) C1M=LIM + DKZ*(DX-1 .OCO J'i^O . 5D0*Ch ISIG=-L JTFL=J-1 C/RG=(2.0CO-CX)^D/iRGP*C.5DO GC TO 21 ■ZERO VERTICAL WAVE NUMBER 7 I£ACK=-1 6 ISIG=G GC TO 2 1 S IEACK=1 GC TO 8 ■REAL REGICN IS IF JT CK CX C/ IF C2 CT CI $T GC CT IF IF CA IVs JT CT CI GC A = E = CE CE AC EC CI IF CI GC CI CI CA GC IG = 1 {5T0P) GO TO 11 FU=J Z = DSQ = CARG RG=CA ($TQ = 2.CC hETU = NT=DT CP=.T TO 1 1 = DMC (DAR (I WE PG=C. ELL=I Flj=J FETU= NT = DT TO 1 CSIN{ CCCS( L = CEX fL=CE = DEL = = DEL-^ NTP=D (DIN NTP=D TO 1 NT = DI NT=DI RG=0. TO 2 ST(CKZ2) P/(CKZ+DARGP) F.G^CAFGP'i^(DX-1.0C0}*0.5DC Pi GO TG 13 C-D£XP(2.GDC-CARG^DH ) CAT AN ( (D2-1 .CD0)/(D2+1.0CC)) F-tTC RUE. C(DINT,D2P1 ) G-'Ch-CWELL ) LL.EG.KWELL) ccc WELL+1 15, 1^, lA GC TO 28 CPIO- FETU e CTl) CTl) F(DAI XP(-I (A+E (A + e INT-DTl + DAT/iN2(AC,B0) TF.GE.CINT-CTEST) GC TO INTP- 6 NTP NT+12.0D0-CX)*DKZ*C.5D0*DH ceo 1 iRG=^DH) ■DARG-^Ch ) ) + CEiML"(A-B) )-DEML-(A-B) •DTl + DAT/iN2(ACt E.CINT-CTEST) '+0PI02 17 CKZ = DS(.RT(CKZ2) ISIG=+1 CINT=DINT.+ DKZ=*OH IF ($TCF) GO TO 13 JTPU=J GC TO 12 CKZ=CSQRT(DKZ2) CINT=DINT+CKZ*DH 21 CCNTINUE •DEPTH LCCF ENDS 154 c c c c c IF (DKZ2) 22,27,25 ■DECAYING SCLtTICN 12 CKS=DARGP C/iFG=(C/5p,G-C.500*DARGP )^DH«2.0DC IF (CARG-DEXPU) 2b,24,2A 2 3 C5=(DKS-DRCRb*^DKBl*0EXF(-DARG)-«2.0CO/(OKS+CPCPE=PDKe) DThETL=CATAN(( l.0D0+D3)/{ 1.0D0-D2)) GC TO 27 2^ CTFETL=DPIQ4 GC TO 27 JTFL=NFI IF (DKB.FC.O.ODG) GC TO 26 CTFETL=Ca*TAN(DKZ/ (CK8*CRQRB) ) GC TO 27 ECTTOM REFLECTED It CTf-ETL=CFIC2 COMPUTE DIFFERENCE 26 CS=CINT+DTFETL CINT=DINT-DThETU DIFF = DS-CF'>DPI MnELL = I i\ELL RETURN Il>ELL = IyNElL+l GC TO 2^ END 155 BIBLIOGRAPHY Arthur D. Little, Inc., Report No. C-70673, Design of an Array to Excite Individual Norinal Modes in the BIFI Range, 31 January 1969. Biot, M. A. , "General Theorems on Equivalence of Group Velocity and Energy Transport," The Physical Review, vol. 105, 1129, 1957. Bucker, H. P. and Morris, H. E. , "Normal Mode Intensity Calcula - tions for a Constant -Depth Shallow -Water Channel, " J. Acoustic Soc. Amer. , vol. 38, p. 1010, 1963. Bucker, H. P. and Morris, H. E. , "Epstein Normal-Mode Model of a Surface Duct, " J. Acoustic Soc. Amer. , vol. 41, p. 1475, 1967. Churchill, R. V. , Fourier Series and Boundary Value Problems, McGraw-Hill, 1941. Naval Underwater Systems Center Technical Report 4319, Computer Programs to Calculate Normal Mode Propagation and Applica- tions to the Analysis of Explosive Sound Data in the BIFI Range, by W. G. Kanabis, 6 November 197Z. Kornhauser, E. T. , and Raney, W. P. , "Attenuation in Shallow "Water Propagation Due to an Absorbing Bottom, " J. Acoustic Soc. Amer., vol. 27, p. 689, 1955. North Atlantic Treaty Organization, SACLANT ASW Research Centre Technical Report No. 203, An Expansion Technique for Render- ing the WKBJ Method Uniformly Valid, by V. R. Lauvstad, 15 December 1971. Naval Research Laboratory Memorandum Report 2381, A Normal Mode Computer Program for Calculating Sound Propagation in Shallow Water with an Arbitrary Velocity Profile, by A. V. Newman and F. Ingenito, January 1972. Officer, C. B. , Introduction to the Theory of Sound Transmission, McGraw-Hill, 1958. Schiff, L. I., Quantum Mechanics, McGraw-Hill, 1955. Tolstoy, I. and Clay, C. S,, Ocean Acoustics. McGraw-Hill, 1966. 156 Williams, A. O. , Jr., and Horne, W. , "Axial Focusing of Sound in the Sofar Channel," J. Acoustic Soc. Amor., vol. 41, p. 189, 1967, "Williams, A. O. , "Normal Mode Methods in the Propagation of Underwater wSound, " Underwater Acoustics, ed. by R. W. B. Stephens, Wiley-Interscience, 1970. 157 INITIAL DISTRIBUTION LIST No. Copies 1. Defense Documentation Center 2 Cameron Station Alexandria, Virginia 22314 2. Library, Code 0212 2 Naval Postgraduate School Monterey, California 93940 3. Professor Glenn H. Jung, Code 58Jg 3 Department of Oceanography Naval Postgraduate School Monterey, California 93940 4. Associate Professor Alan B. Coppens, Code blCz 3 Department of Physics and Chemistry Naval Postgraduate School Monterey, California 93940 5. Department of Oceanography 3 Naval Postgraduate School Monterey, California 93940 6. The Oceanographer of the Navy 1 The Madison Building 732 North Washington Street Alexandria, Virginia 22 314 7. LT. Kirk E. Evans, USN 3 SMC 2886 Naval Postgraduate School Monterey, California 93940 8. Dr. Robert E. Stevenson 1 Scientific Liaison Office Scripps Institute of Oceanography La Jolla, California 92037 9. Office of Naval Research, Code 480 1 Arlington, Virginia 10. NAVOCEANO Library 1 U.S. Naval Oceanographic Office Washington, D. C. 20390 158 No. Copies 11. LCDR Ernest Young 2 Office of Naval Research Liaison Officer Room 411, Herrmann Hall East Wing Naval Postgraduate School Monterey, California 93940 12. Associate Professor Warren W. Denner, Code 58Dw 1 Department of Oceanography- Naval Postgraduate School Monterey, California 93940 13. Assistant Professor Robert H. Bourke, Code 58Bf 1 Department of Oceanography Naval Postgraduate School Monterey, California 93940 159 UNCLASSIFIED SECURITY CUASSinCATlON OF THIS PAGE (Wh»n Dbih Enterad) REPORT DGCUMENTATIOW PAGE READ INSTRUCTIONS DErG"$E COMJ'LKTING FORM 1. REPORT NUMBER 2. GOVT ACCESSION NO. 3. RECIPIENT'S CAT ALOG NUMBER 4. TITLE (and Subtitle) Two Methods for the Numerical Calculation of Acoustic Normal Modes in the Ocean 5. TYPE OF REPORT & PERIOD COVERED Master's Thesis; September 1973 6. PERFORMING ORG. REPORT NUMBER 7. AUTHOR^*; Kirk Eden Evans a. CONTRACT OR GRANT NUMBERf*; 9. PERFORMING ORGANIZATION NAME AND ADDRESS Naval Postgraduate School Monterey, California 93940 10. PROGRAM ELEMENT, PROJECT, TASK AREA * WORK UNIT NUMBERS 11. CONTROLLING OFFICE NAME AND ADDRESS Naval Postgraduate School Monterey, California 93940 12. REPORT DATE September 1973 13. NUMBER OF PAGES 161 t4. MONITORING AGENCY NAME A ADDRESSC// (