OS2IFD: A microcomputer implementation of the parabolic equation for predicting underwater sound propagation

Computers & GeosciencesVol. 17, No. 6, pp. 731-757, 1991 Printed in Great Britain

0098-3004/91 $3.00+ 0.00 PergamonPress pie

OS2IFD: A MICROCOMPUTER IMPLEMENTATION OF THE PARABOLIC EQUATION FOR PREDICTING U N D E R W A T E R S O U N D PROPAGATION J. S. ROBERTSON,1W. L. SIEGMANN,2 and M. J. JACOBSON2 Department of Mathematical Sciences,United States Military Academy, West Point, NY 10996-1786and 2Department of Mathematical Sciences, Rensselaer Polytechnic Institute, Troy, NY 12180-3590, U.S.A.

(Received 6 October 1989; accepted 17 April 1990; received for publication 14 March 1991) Abstract--Underwater acoustic signals are powerful scientifictools for probing large regions of the ocean which would otherwise be inaccessible. Oceanographers and other scientists use acoustic energy as a mechanism to examine the structure of ocean regions for a variety of purposes. Predicting the behavior of sound in different types of ocean environments is an extraordinarily difficult problem which has been studied intensely for many decades. The parabolic approximation, introduced to the oceanographic community more than a decade ago has proven to be a powerful and effectiveocean acoustics propagation model. Whereas these models were once exclusively run on mainframe computers, the advent of fast microcomputer chips, together with operating systems that can exploit the powerful features of these chips, now makes personal computers an attractive tool for performing many propagation prediction computations. This paper describes a full-featured version of one such widely used underwater acoustic propagation model which runs on PCs under the OS/2 operating system.

Key Words: Ocean acoustics, Sound propagation, Numerical modeling, Parabolic approximation, Finite differences.

INTRODUCTION

THE PARABOLIC APPROXIMATION

Underwater acoustic signals are powerful scientific tools for probing large regions of the ocean which would otherwise be inaccessible. Oceanographers and other scientists use acoustic energy as a mechanism to examine the structure of ocean regions for a variety of purposes. Predicting the behavior of sound in different types of ocean environments is an extraordinarily difficult problem which has been studied intensely for many decades. Of the many modeling developments that have emerged from this effort, the parabolic approximation, introduced to the oceanographic community more than a decade ago (Tappert, 1977) has proven to be a powerful and effective ocean acoustics propagation model. One numerical implementation of the parabolic approximation which has seen broad use throughout the ocean science community is IFD, an implicit finitedifference code which solves the parabolic approximation (Lee and McDaniel, 1987). It has been run on a variety of large computing machines, and numerous research papers using IFD to study the environmental effects of ocean variations on underwater acoustic transmissions have been published. The advent of fast microcomputer chips, together with operating systems that can exploit the powerful features of these chips, now makes personal computers (PCs) an attractive tool for performing many propagation prediction computations. The purpose of this paper is to describe a full-featured version of IFD which runs on PCs under the OS/2 operating system.

Parabolic approximations (or parabolic equations) are obtained typically from the Helmholtz equation

cAoeo ,7/~--A

731

V2p(r,O,z) + k~n2(r,O,z)p(r,O,z) = 0,

(!)

which governs the acoustic pressure field p in a steady, quiescent medium resulting from an harmonic point source of frequency f. In Equation (1), ko = 2nf/co is the wave number, n = Co/c(r,O,z) is the index of refraction, c(r,O,z) is the sound speed, and Co is a reference sound speed. In the far-field p is assumed to have the form given by

p = u(r,O,z)H~ol)(kor),

(2)

where H~)(ko r) is the Hankel function of the first type of order zero. Equation (2), assumed to hold for sufficiently large values of r, supposes that there are only outward traveling waves, that is no backscattering. The quantity u(r,O,z) is a slowly varying function of position which modulates the Hankel function. If Equation (2) is substituted into Equation (1), and appropriate conditions are satisfied, the following "standard PE" can be shown to hold: 2/k0u~+ u= + k2o(n2 - 1)u = 0.

(3)

Because no 0-derivatives are in Equation (3), we suppress the 0-dependence of u and write u(r,z). Details on the derivation of Equation (3) are given in Tappert (1977), Robertson, Siegmann, and Jacobson (1985), and elsewhere.

732

J.S. ROBERTSON,W. L. SIEGMANN,and M. J. JACOBSON

There are several numerical algorithms available for solving Equation (3), together with an appropriate set of boundary and initial conditions. In particular, one using implicit finite-differencing of the partial derivatives has been distributed widely in the under water acoustics community. It utilizes a Crank-Nicolson scheme, together with a fast tridiagonal solver (Carnahan, Luther, and Wilkes, 1969), to march the solution in range with a formal order of accuracy of O[(Ar)2 + (Az)4], where Ar and Az are the range and depth step sizes, respectively. This method is popular for parabolic systems because of its accuracy and absolute stability. One advantage to this approach is its particularly convenient ability to model horizontal as well as irregular interfaces between layers with different acoustic properties. Additional details of this algorithm are given in Lee and McDaniel (1987). IMPLEMENTATION OVERVIEW As noted, the underlying mathematical and numerical theory of the parabolic approximation and the implicit finite-difference scheme used to solve it have been discussed thoroughly elsewhere (Tappert, 1977; Robertson, Siegmann, and Jacobson, 1985; Lee and McDaniel, 1987, and references therein). The mainframe implementation, known as IFD, contains a number of powerful and flexible features which permit its application to a broad variety of ocean acoustic environments. F o r instance, an input stream provides a mechanism for the user to prescribe source data, specify selected numerical parameters such as mesh size, describe many types of complex propagation environments, and select a variety of output options. I F D implements pressurerelease and rigid boundary conditions on flat or sloping boundaries, and matches the wave field from the water column into any number of subbottom layers with different densities, sound speeds, and volume attenuations. It also contains a feature to apply an artificial absorbing layer (beneath subbottom layers) which is used to enforce a pressurerelease bottom boundary in many applications. We have preserved the flexibility and power of the mainframe version; at the same time, we have implemented a number of improvements and enhancements intended to improve the portability and ease of use of this model. Here we will discuss changes made to the mainframe version of I F D and outline its application to underwater acoustic problems on a personal computer. The complete source code is contained in three files, each of which is listed in the Appendix. The main program is contained in OS2IFD.FOR. The remaining pair are "include" files: C O M M O N . I F D contains the common declarations for the main program, and D O C . I F D serves as a useful documentation reference to the variables, flow-of-control, and important input-output considerations.

Note first that double precision real (REAL*8) and complex (COMPLEX*I6) variables and functions are used throughout the program. We have determined that accuracy attained in this way is worth the additional computational expense entailed by the extended precision. On the other hand, we decided to use 16-bit integers (INTEGER*2) because their use can yield noticeable improvements in running time and because longer integers are never necessary. Ample use also is made of P A R A M E T E R statements which make subsequent modifications easier to implement and provide the compiler with optimizing information. Most of the code has been rearranged to exploit the I F - T H E N - E L S E control structures available in F O R T R A N 77. Common blocks have been rcconfigured so that large arrays are separated, allowing the compiler to generate more efficient code. A number of nonportable F O R T R A N language constructs, present in the mainframe version of IFD, have been replaced with code conforming to accepted FORTRAN-77 standards. Array sizes can be controlled through the include file C O M M O N . I F D . The size for the main arrays, 214 elements, should be adequate for most problems. We have tested problems that used larger arrays (2 ~6elements); nevertheless, we recommend retaining the sizes listed in C O M M O N . I F D , because it is not efficient to require the code to manage unnecessarily large data structures. In order to facilitate the use of OS2IFD with batch files, an important file naming feature has been added. At the beginning of program execution, the code reads four file names from a file termed N A M E S . I F D . These file names, consisting of up to 32 characters, can include path names and drive names as well. The first line of N A M E S . I F D lists the name of the input file, which contains all the execution information required by OS2IFD. The second file name is used to hold a summary of the input information read by the code, and serves to record a number of optional outputs, such as the starting field, range information, etc., that may be requested by the user via the input file. The third file contains the complex-valued solution to the PE, whereas the fourth file contains the propagation loss in decibels. The complex solution output file should only be generated if it is needed specifically, for example for checking the conformance of an accuracy benchmark (Robertson, 1989). This file will remain empty and be deleted at program termination by setting either the variables W D R or WDZ, contained in the input file, to a value of zero. Additional details are contained in the Appendix. OS2IFD contains no facility for the graphical visualization of output. Several modestly priced commercial plotting packages provide ample ability to examine propagation predictions on the screen, as well as preparing high-quality hardcopy output on plotters, printers, and other peripheral devices.

Predicting underwater sound propagation Table 2. Data summary file for Example 1

Table 1. Input runstream for Example I 500.0, 50.0, 0.0, 0, 0.0, 250.0, 500, 0, 3, 0 25000.0, 5.0, 0.0, 0.0, 50.0, 50.0, 0, 0, 0 0.0, 100.0 25000.0, 100.0

OS/2-IFDPE Solution--Double Precision Version Gaussian Starting Field Frq = 500.0Hz ZS = 50.0 m CO = 1500.0m/sec RA = 0.0m ZA = 250.0m N = 500 DR = 5.0 m WDR -0.0 m PDR = 50.0 m DZ = 0.5000m WDZ 0.5000m PDZ = 50.0000m ITYPEB = 3 ITYPES = 0 RMAX = 25000.0 m

- 1.0, - 1.0

0.0 0 2 1.0 0.0, 1500.0 100.0, 1500.0 200.0, 1.2, 1.0 100.0, 1550.0 200.0, 1550.0 I00.0, 1.0, -

=

SYSTEM CONSIDERATIONS The implementation described here is targeted for microcomputers based on the Intel 80286/386 ® architecture and which run the OS/2t operating system. OS/2 is a descendant of MS-DOS, and is rapidly seeing a growth of support throughout the software industry. It is a particularly attractive operating system for scientific applications because it preserves the appearance and style of MS-DOS, supports many MS-DOS applications through a feature known as the DOS compatibility box, and, perhaps most importantly, permits data structures and code segments which exceed 640 kilobytes (kbyte) in size. This last feature, together with chips that run at clock speeds in excess of 25 MHz, makes feasible the porting to these smaller machines of many codes which heretofore ran only on minicomputers or bigger machines. The examples which follow were executed on an IBM PS/2 System 70 Model A21 with an 80387 math coprocessor running at a clock speed of 25 MHz. The operating system version was OS/2 Extended Edition 1.1. The machine was equipped with 16 megabytes (Mbyte) of random-access memory (RAM). The performance numbers to be reported are sensitive to the hardware configuration. For example, a configuration with only 4 Mbyte of RAM, or one with an 80286/287 combination, would execute more slowly. The program performance also is sensitive to presence of other concurrent processes. That is, if the multitasking ability of OS/2 is being used to accomplish other tasks while OS2IFD is running, there will be a predictable performance degradation. The compiler used was Microsoft FORTRAN 4.1, a full FORTRAN-77 implementation. The only system-dependent calls utilized were system clock reads which were used to estimate the number of CPU seconds used per run. This call is noted in the listing each time it appears. Once these are deleted, the source code should compile on any FORTRAN-77 compiler. Finally, we note that OS2IFD can be compiled to run under MS-DOS, although allowable array sizes are typically too small to be of much practical use owing to the memory constraints of that

733

Layer Max Depth (m) Density (g/cm**3) Att (dB/wl) 1

100.000000

! .000000

- 1.000000

2 200.000000 ! .200000 1.000000 3 250.000000 1.200000 1.000000 End of run. OS2IFD terminated after 708.0CPU seconds. operating system. Nevertheless, users may care to experiment, and it may be possible to perform useful but limited propagation predictions under that operating system. EXAMPLES The following two example problems serve as speed benchmarks (Robertson, 1989). Both problems are drawn from a recent monograph (Lee and McDaniel, 1987). The input files are listed in tables, and a sample illustration of the solution generated for each example is discussed briefly. The first example models an isospeed, shallowwater sound channel. The input data file is given line-by-line in Table 1. In this example a continuous wave (cw) signal at 500 Hz at a depth of 50 m is propagated in an underwater sound channel to a range of 25 km. The water column is 100 m deep and overlies a fluid bottom. Both the density and sound speed are discontinuous across the interface: the density increases from 1.0 to 1.2 g cm -3 whereas the sound speed jumps from 1500 to 1550ms -1. The bottom layer also has a mild attenuation factor of

-40 -5O~ v"o

~O D S21F' BE~NCHMAE R,X KAMPL1E

--60 -70

_z

-80 -90

-10C

5

10

15

20

25

R~OE (KM)

?Trademark of IBM Corp. :~Tradcmark of Microsoft Corp.

Figure 1. R e l a t i v e intensity v s r a n g e for E x a m p l e R u n n i n g time = 708 sec.

I.

J.S. ROBERTSON, W. L. SIEGMANN, and M. J. JACOBsoN

734

Table 3. Input runstream for Example 2

-60

100.0, 30.0, 0.0, 0, 0.0, 1200.0, 600, 0, 3, 0 20000.0, 2.0, 0.0, 0.0, 50.0, 90.0, 0, 0, 0 0.0, 240.0 20000.0, 240.0

-65

-

1.0,

0.0 0 2 240.0, 0.0, 120.0, 240.0, 512.0, 240.0, 512.0,

- 1.0

O$21FD BENCHMARKEXAMPLE2

-70

~

75

~j -80

1.0, 0.0 1500.0 1498,0 1500.0 2.1, 0.0 1505.0 1505.0

-85 -90

i

0

,

i

10

15

20

RANGE (KM)

Figure 2. Relative intensity vs range for Example 2. Running time = 2208 sec.

1.0 dB/wavelength. Table 2 contains a listing of the output summary file.This file is useful for verifying the correctness of information supplied in the input file,and, when certain options are selected, contains detailed information on the sound speed profiles and interface-bottom topographies. Figure I illustrates one way to visualize the output from O S 2 I F D for this model ocean channel. Shown in this figure is relative intensity (the negative of propagation loss) plotted as a function of range for a point receiver that is 50 m deep. The running time for this example on the system described earlier was 708 sec, comparable in our experience to running times encountered on small minicomputer systems (e.g. a Prime 850). Note that the figure shows the intensity at only one receiver depth. In fact, the output file from this one run contains a substantial amount of additional data at other receiver depths. The second example models a strong horizontal interface in a moderately deep channel with an inhomogeneous sound speed. Table 3 contains the input data filefor this example, whereas Table 4 contains the data summary. The source frequency is I00 H z and the

source is located at depth 30 m. The range o f the sound channel is 20 km. The bottom layer contains a strong density discontinuity which increases from 1.0 to 2.1 g cm -3 across the interface, but no attenuation. This particular problem uses a smaller range step size than Example 1 and the channel is substantially deeper. As a result, the running time is naturally longer. Figure 2 depicts relative intensity vs range when the receiver is 90 m deep. The running time here was 2208 sec, about three time.s longer than Example 1. We emphasize that each run produces a substantial amount of information, so that the ocean scientist obtains much useful modeling information from each run.

Table 4. Data summary file for Example 2

Acknowledgments--This work was supported by Code

OS/2-IFDPE Solution--Double Precision Version Gaussian Starting Field Frq ZS CO RA ZA N DR WDR PDR DZ WDZ PDZ ITYPEB ITYPES RMAX

= 100.0Hz = 30.0 m = 1499.0 m/see = 0.0m = 1200.0 m = 600 = 2.0 m = 0.0 m = 50.0 m = 2.0000 m = 2.0000 m = 90.0000 m = 3 -0 = 20000.0 m

Layer Max Depth (m) Density (g/cm**3) Att (dB/wl) 1 240.000000 ! .000000 0.000000 2 512.000000 2.100000 0.000000 3 1200.000000 2.100000 0.000000 End of run. OS2IFD terminated after 2208.0 CPU seconds.

SUMMARY

An implementation of a widely used implicit finitedifference model for predicting ocean acoustic propagation presented here has been tested for, and is intended to run under, the 0S/2 operating system. The program, termed OS2IFD, enables oceanographers and other ocean scientists to perform a broad variety of ocean acoustic propagation predictions on microcomputers.

1125OA, Office of Naval Research, and the Signatures, Sensors, and Signal Processing Technology Office, U.S. Army Laboratory Command. REFERENCES

Carnahan, B., Luther, H. A., and Wilkes, J. O., 1969, Applied numerical methods: John Wiley & Sons, New York, 604 p. Lee, D., and McDaniel, S. T., 1987, Ocean acoustic propagation by finitedifference methods: Comp. & Maths. Appls., v. 14, no. 5, p. 305--423. Robertson, J. S., 1989, A classificationscheme for computational ocean acoustic benchmarks: Appl. Acoust., v. 27, no. I, p. 65-68. Robertson, J. S., Siegmann, W. L., and Jacobson, M. J. 1985, Current and current shear effectsin the parabolic approximation for underwater sound channels: Jour. Acoust. Soc. America, v. 77, no. 5, p. 1768-1780. Tappert, F. D., 1977, The parabolic approximation method, /n Keller, J. B., and Papadakis, J. S. ~Is., Wave propagation and underwater acoustics: Springer, Berlin, p. 224--287.

Predicting

underwater

sound

propagation

735

APPENDIX

Source Code for DOC.IFD ******************************************************************

* Implicit Finite-Difference Method for Solving the * Parabolic Equation : U =A*U+B*U ; where *

ZZ

R

*

* *

A=I*.5*K0*(N*N-I)

*

;

AND

B=I*.S/K0

*

****************************************************************** * Based on the IFD model by * D. Lee and G. Botseas, Code 3342 * Naval Underwater Systems Center * New London, Connecticut 06320

* * * *

******************************************************************

* * * * * * *

Double Precision Version developed for Microsoft Fortran 4.1 running under OS/2 v. I.I J.S. Robertson Department of Mathematics U.S. Military Academy by West Point, New York 10996-1786 W.L. Siegmann and M.J. Jacobson * Department of Mathematical Sciences * Rensselaer Polytechnic Institute * Troy, New York 12180-3590

by

* * * * * * * *

* *

******************************************************************

* This source code is available on diskette. Contact J.S. Robertson at the above address. ******************************************************************

*

*

*

* NOTE: This code uses one non-standard system call, GE~TIM. * This subroutine is called to estimate the number of CPU * seconds used during progrum execution. This call must be * deleted or changed in order to compile under a different * compiler.

* * * * *

******************************************************************

---Variables **

ACOFX

** **

ACOFY ALPHA

**

ATT

,

** **

BCOF BETA

** **

BOTX BDTY

**

BTA

**

CO

**

CSVP

** **

DEG

** **

DR1

**

DR

DZ DZZ

*

** EPS ** FRQ

and Constants---

- Coefficient 'A' at bottom - at advanced range - Coefficient 'A' at bottom - at present range - Volume attenuation - dB/meter Attenuation coefficient for artificial absorbing layer at bottom or surface - Coefficient 'B' - Range independent - Array - attenuation in layers - db/wavelength - Complex pressure at bottom at advanced range ra+dr Complex pressure at bottom at present range ra Array - partial solution of system of equations Reference speed of sound - meters/sec - Array - sound velocity - meters/sec - Conversion factor - degrees/radiun (180.DO / PI) - Range step - meters - Last dr used in routine diag - Depth increment of solution - meters - Depth increment for adjusting layer depths in Sloping bottom - Parameter - I.ODO-IO - Frequency - hz

(1)

(k r) 0 0 ** HNKL - External function - computes hankel function ** IBOT Bottom depth print flag = 0 - do not print bottom depths = I - print bottom depths * * IDIAG ~ Diagonal update flag ** IHNK - Hankel function flag = 0 - hankel function not used. lO*log(r) added to solution. = I - starting field divided by hankel function. solution multiplied by hankel function before computing propagation loss. * * ILYR - Index for arrays beta, zlyr, and rho ** IPZ - Every ipz'th value in depth is printed ** ISF Starting field flag * = 0 - program generates gaussian starting field * at range = 0.0. see subroutine sfield. ** ,

HNK

-

Hankel function H

* * • * *

• • * * *

• • * •

*

* * *

, • • * * * *

* • * *

* , *

* • • * *

* *

736

J . S . ROBERTSON, W . L. S1EOMANN, a n d M . J. JACOB,TON

C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C

user supplies starting field, see subroutine ufield. ISFLD - Starting field print flag = 0 - do not print starting field = I - print starting field ISVP - Svp print flag = 0 - do not print sound velocity profile = I - print sound velocity profile ITRK - Index for array track ITYPEB- Type of bottom 0 - Rigid bottom - program supplies bottom condition - user supplies bottom condition I - Not rigid - see subroutine bcon - absorbing layer introduced 2 - Not rigid - follows contour of bottom 3 - Not rigid - absorbing layer introduced - but - bottom of absorbing layer kept flat ITYPES- Type of surface = 0 - pressure release, scon sets sury and surx = 0.0 = .he. 0 - user inserts code in scon to compute sur =

C

**

** * *

** ** * *

*

$

i

-

And surx

- Index increment of receiver solutions to be written On disk ** I X S V P Array of pointers which point to entries in csvp an

**

INZ

Zsvp

**

KSVP

-

* * * *

* *

* * * *

* * * * *

* * * * * * *

* *

* *

Ixsvp(1) points to bottom depth and speed in layer i * Ixsvp(2) points to bottom depth and speed in layer 2 * Etc. -* Svp profile flag * Profile is read from input file - see subroutine SVP * User supplies profile 1 - see subroutine usvp * User supplies profile 2 - see subroutine usvp *

* *

0 I

-

*

2

-

1~

-

*

-

* * * *

*

**

** , **

** * **

** **

User supplies profile n - see subroutine usvp Integer variable that is local to a do loop, Added for execution speed. The underscore is is replaced by a digit. MLYR Temporary - number of layers involved in specific Calculation MM Index - m m ÷ l points to first value of artificial Absorbing Layer in array u M X L Y R - P a r a m e t e r - maximum number of layers - max dimension of arrays MXN Parameter - maximum dimension of c, x, y, r12 a n d u Arrays MXSVP - Parameter - maximum dimension of arrays csvp and zsvp MXTRK - Parameter - maximum dimension of array track LOOP_

N

N1 NA NDR NIU NLYR

NOLD N0U NPU NqU NSVP NWSVP NZ DLDR PDR PDZ PI PL R1 R2 RI2 RA

Number of equi-spaced points in u Includes bottom point - does not include surface Point Diagonal elements nl thru n will be computed Number of points in absorbing layer - Parameter - Range step file number - Parameter - input unit number - Number of layers - Number of receiver depths on entry to routine crnk - Parameter - output unit number - Parameter - printer unit number - Parameter - propagation loss file number - Number of points in csvp and zsvp - Number of points in layer I svp N u m b e r of solution depths to be written on disk - Range increment at start of problem - R a n g e increment at which solution is printed - meters - Depth increment at which solution is printed - meters Parameter - The value of pi Propagation loss - dB range at which bottom depth is available - meter Next range at which bottom depth is available - meter - Array of density ratios Horizontal range of starting field from source -

* * *

* * * * *

* *

* * • * * * * • * * • , • * * • *

* * * * *

* •

Predicting underwater sound propagation C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C

* , *

RA+DR ** KHO ** RMAX -

**

**

RSVP

-

**

SURX

-

**

SURY

-

TF_.~P * * THETA **

* * TRACK

-

** U ,

** *

W D R

-

**

NDZ

-

* * *

**

X

**

XK0

-

**

XPR

-

**

XWR

-

**

Y

**

Z1

-

**

Z2 ZA

-

**

-

*

**

ZERO

-

**

ZLYR

-

** ZI ** ZS **

-

Z S V P

-

Meters Ra is set to 0.0 if starting field is gaussian, ra is Incremented as solution is marched out in range. Range to which solution is to be advanced - meters Array - density in layer Maximum range of solution - meters Next range at which svp is available - meters Complex pressure at surface at advanced range ra+dr Complex pressure at surface at present range ra Temporary variable - complex Slope of bottom - radians .eq.O -- flat bottom .gt.O -- shallow to deep .it.O -- deep to shallow 2 dim. array - range and depth of water - meters Array - complex acoustic pressure field at mth range step R a n g e step at which solution is written on disk Meters Depth increment at which solution is written on disk Disk in meters Wdz should be selected so that plot program does not Interpolate between widely spaced receivers. Array - main diagonal of matrix at advanced range Reference wave number Range at which solution is printed - meters Range at which solution is written on disk - meters Array - main diagonal of matrix at present range Depth of water at range rl - meters Depth of water at range r2 - meters Depth of field at range ra - meters Initial depth of starting field at range ra is As follows: If itypeb = 0 or I, za is maximum depth oT Bottom-most sediment layer at initial range of Starting field, if itypeb = 2 or 3, za is maximum Depth of artificial absorbing layer at initial Range of starting field, program inserts layer. Rho and beta are obtained from layer above. Speed is bottom-most speed from layer above. As solution progresses, Za is updated if ocean bottom not flat. if itypeb = S Bottom of absorbing layer remains flat. Parameter - the value 0.0 Array - depth to bottom of layer - meters Depth of receiver 'i' - meters Source depth - meters Array - depth of sound velocity - meters

737

, * *

* • *

* * * * * * * *

* * * * *

* , * • * , * *

* • • , * ,

* * * * *

* *

, * * , , , , ,

****************************************************************** * Input ******************************************************************

,

* * * * *

* * • * ,

In order to generate maximum flexibility, such as running a chain of batch jobs, input and output file names are contained in a file called NAMES.IFD. Under 0S/2 V 1.1, this file must be in the directory from which the program is being executed.

*

Line i:

Data filename

attached to NIU

* * * *

Line 2:

Output filename I attached to NPU, which will contain general program information, time stamp, a summary of input data, program options, and.error messages.

* • * ,

* *

Line 3:

Output filename 2 attached to NOU, will contain the solution vector u(n).

* ,

* *

Line 4:

Output filename the propagation

* ,

S attached to NQU, will contain loss.

* This file is read from Unit 8. If this unit is not available * on your system, then you must change the code in * subroutine READER.

,

* , ,

738

J.S. P,OBERTSON, W. L. SIEGMANN, and M. J. JACOBSON

* *

I 2 3 4

unit : : : :

* *

Line N

:

Line

N+I:

-1,-1

* *

Line Line Line Line Line Line

N+2: N+3: N+4: N+5: N+6: N+7:

RSVP KSVP NLYR ZLYB(I) ,B/i0(I),BETA(I) ZSVP(1) ,CSVP(1) ZSVP(2) ,CSVP(2)

* *

* *

* *

* *

Input Line Line Line Line

number ffi NIU FRQ,ZS,C0,ISF,BA,ZA,N,IHNK,ITYPEB,ITYPES RMAX,DR,WDR,WDZ,PDR,PDZ,ISFLD,ISVP,IBOT RI,ZI ** R2,Z2 * Bottom profile ** Range, Water Depth (meters)

################# ## ## ## # ## Repeat # # for # each #

*

##

*

__

N+M:

* * * *

*

Line

* *

**

* * *

*

* *

ZSVP(J)

*

* * * *

Repeat * for * each * pro* file *

#

*

#

*

#

#

*

##

#

*

_

,CSVP(J)

# # # # ## ## ## ##

layer

#################

*

****************************************************************** *** Quick reference and notes for Line input *** Units: Meters and meters/sec except as noted *****~*****~**~*****************~***************~**************** * F R Q = Frequency (HZ) * ZS = Source depth * CO = Reference sound speed. If C0 = 0.0, CO is set to * average speed in first layer. * ISF ffi starting field flag. 0 ffi gaussian, i ffi user field. * if isf = 0, ra is set to zero. * RA = horizontal range from source to starting field. * ra is set to 0.0 if isf = 0. * ZA = depth of starting field at range ra. if za = 0.0, za * set * to max depth of bottom layer in first profile, if * ITYPEB = * *

2

* * * *

* * *

* *

* * * *

* * * *

* *

N

or

3

and

za

= 0.0,

za

is

set

to

(4/3)*max

depth

of

bottom layer, if itypeb ffi 2 or 3 and za not zero, the artificial bottom layer is extended to za meters provided that za is greater than or equal to max depth of bottom layer in first profile. ffi number of equispaced receivers in starting field, if =

O,

n is set so that the receiver depth increment is less than or equal to I/4 wavelength, if n is greater than mT~, n is set to mxn. IHNK = hankel function flag. ihnk ffi 0, don't use h a z e l function.

* * * *

* *

* * * * * ,

* , * *

* * * *

• * *

* ,

* . * ,

ihnk = 1, divide starting field by hankel function, the* multiply the solution field by hankel function before * computing propagation loss. if starting field is * gaussian, , ihnk should be set to 0. if starting field is elliptic,* ihnk should be set to i. * ITYPEB = type of bottom processing. = 0 - rigid bottom, program supplies bottom condition. = I - user supplies bottom condition, user ,rites subroutine bcon = 2 - artificial absorbing bottom introduced, follows contour of water/bottom interface. ffi 3 - artificial absorbing bottom introduced, bottom of layer kept flat. ITYPES = type of surface = 0 pressure release, scon sets sury and surx = 0.0 not 0 - user inserts code in scon to compute sury and - surx

Predicting underwater

* *

C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C

* *

* * * * *

* * * * * *

* * * *

* * * * * *

RMAX = maximum range of solution DR = range step. if dr = O, dr is set to I/2 wavelength. if bottom of problem is not flat, dr is recomputed so that max depth is either incremented or decremented by dz. solution is computed every dr meters. WDR = range step at which solution is written on disk. if w not O, a n output disk file is assigned. rounded to nearest dr. WDZ = depth increment at which solution is written on disk. rounded to nearest dz. PDR = range step at which solution is printed. rounded to nearest dr. PDZ = depth increment at which solution is printed. rounded to nearest dz. ISFLD = 0 - don't print starting field. = I - print starting field. ISVP = 0 - don't print sound velocity profile. = 1 - print sound velocity profile. IBOT = 0 - don't print bottom depths. = I - print bottom depths. R1,Z1 R2,Z2 RSVP KSVP

= = = =

*

=

* *

=

* * NLYB = * * ZLYR(I) * RHfl(I) * BETA(I) *

* * * * * *

sound propagation

ZSVP CSVP ZSVP(1) CSVP(1) ZSVP(J) CSVP(J)

* * * *

bottom profile, first range and depth of water. etc. (-I,-i) marks the end of the bottom profile. range o f svp svp flag. 0 - svp in input runstream not zero - profile (Lines n+4 thru n+m) is supplied by user. user writes subroutine usvp. ksvp may be used in computed goto statement t o transfer control in usvp. number of layers, if itypeb = 2 or 3, program inserts an artificial layer and increments nlyr by I. = max depth of layer i in profile density in layer i ( g / c m * * 3 ) = attenuation in layer i (db/wavelenEth) if beta(i) is negative, attenuation is computed. = depth array for sound speed = speed array for sound speed = depth of water to top of layer i = speed of sound at top of layer i = depth of water t o bottom of layer i = speed of sound at bottom of layer i if only one svp inputted, it is used thru entire problem.

if more than one svp inputted, last svp is used thru remainder of problem.

********%******************%********$*****************************

Source Code for COMMON.IFD PARAMETER(PI = 3. 141592653589793D0) PARAMETER (MXN = 4096, MXLYR = 128, MXTRK = MXLYR, MXSVP = 128) COMPLEX*16 ACOFX, ACDFY, BCOF, BDTX, BOTY, BTA, HNK, HNKL, SURX, 1 SURY, TEMP, U, X, Y

DIMENSION BTA(MXN), U(MXN), X(MXN), Y(MXN), R12(MXN), 1

BETA(MXLYR), CSVP(MXLYR),

2

TRACK(MXLYR,2),

IXSVP(MXLYR), KHO(MXLYR),

ZLYR(MXLYR),

ZSVP(MXLYR)

COMMON /IFDCM0/ 1 ACOFX, ACOFY, BCOF, BOTX, BOTY, SUP/, SUBY, 2 ALPHA, CO, DR, DR1, DZ, FRq, BA, RSVP, THETA, XKO, ZA, ZP, 3

ZS, IHNK, ISF, ITYPEB, ITYPES, KSVP, N, NI, NLYR, NSVP,

4

NWSVP, DMI, DM2 /IFDCMI/ BTA /IFDCM2/ U /IFDCM3/ X /IFDCM4/ Y /IFDCMS/ R12 /IFDCM6/ BETA, CSVP, IXSVP, KHO, TRACK, ZLYR, ZSVP

COMMON COMMON COMMON COMMON COMMON COMMON

739

* * * *

* * * * *

* * * * * *

* * * *

* * * * * * * * *

* * * * * * *

* * * * * * * • * *

740

J.S. ROBERTSON, W. L. SIEGMANN, and M. J. JACOBSON

Source Code for OS2IFD.FOR PROGRAM OS2IFD IMPLICIT REAL*8 (A - H, O - Z), INTEGER*2 (I - M) $INCLUDE : 'DOC. IFD ' C PARAMETER (NIU = 7, NQU = 8, NOU = 9, NPU = i0) $ I N C L U D E : 'COMMON. IFD ' CHA/RAC~*32 DATA, SUM~Y, VECTOR, PROPLS CHARACTER*7 STWORD INTEGER ITIME(4) COMMON /FILES/ DATA, SUMMRY, VECTOR, PROPLS PARAMETER (DEG - 180.D0 / PI, EPS = 1.0D-10, ZERO = O.ODO) C C *** READ INPUT PARAMETERS CALL GETTIM(ITIME(1), ITIME(2), ITIME(S), ITIME(4)) CPUTIM = 3600.0D0 * ITIME(1) + 60.ODO * ITIME(2) * ITIME(3) CALL READER OPEN (NIU, FILE-DATA, STAI~JS='OLD') OPEN (I~PU, FILE'SUM~Y, STATUS-'UNKNOWN') OPEN (NOU, FILE-VECTOR, STATUS='UNKNOWN') OPEN (NQU, FILE=PROPLS, STATUS" UNKNOWN' ) C REWIND (NIU)

REWIND(NPU) REWIND (NOU)

REWIND(NQU) C READ(NIU, *)FRO, ZS, CO, ISF, RA, ZA, I N, IHNK, ITYPEB, ITYPES BEAD(NIU, *)BMAX, DR, WDR, WDZ, PDR, PDZ, ISFLD, 1 ISVP, IBOT C C C C

C C

C

C C

C

C

C C

*** IF GAUSSIAN STARTING FIELD, BA MUST BE O. IF (ISF .Eq. O) IRA = 0.0

*** READ BOTTOM PROFILE - RANGE,DEPTH DO 70 I = I, MXTRK RFAD(NIU, *) TRACK( I , 1) ,TRACK( I , 2) ITRK= I *** END OF PROFILE? IF (TRACK(I,I) .LT. 0.0) GO TO 80 *** NO 70 CONTINUE *** ERI~3R? 80 IF (ITRK .EO. 1 .OR. ITRK .GT. MXI~RK) WRITE (NPU,120) MXTRK 120 FORMAT(IX, 'Error. Bottom profile missing or too many points. ', 1 ' Max i s ', I5) *** GO TO EXIT POINT GO TO 9 9 9 9 ENDIF *** NO. EXTEND LAST DEPTH BEYOND MAX RANGE. TRACK(ITRK,1) - 1.0D38 TRACK(ITRK,2) - ~RACK(ITRK - 1 , 2 ) ITRK - 1 R2 = TRACK(ITRK,1) Z2 = TRACK(ITRK,2) *** FIND BOTTOM SEGMENT WHICH CONTAINS STARTING RANGE 9 0 R 1 = R2 Zl = Z2 ITRK = ITRK + 1 IF (ITRK .GT. MXII~K) THEN WRITE (NPU,120) MA"FRK *** GO TO EXIT POINT GO TO 9 9 9 9 ENDIF R2 - TRACK(ITHK,I) Z2 = TRACK(ITRK,2) *** ADVANCE TRACK IF NECESSARY SO THAT RI.LE.RA.LT.R2 IF (RA .GE. R2) GO TO 90 *** R1 MUST BE LE RA IF (R1 .GT. RA)

WRITE (NPU, 100)

Predicting underwater sound propagation I00

FORMAT (IX, 'Depth of bottom at initial range missing ~) *** GO TO EXIT POINT GO TO 9999 ENDIF *** COMPUTE SLOPE 0F BOTTOM THETA = DATAN2(Z2 - Z 1 , R 2 - RI)

C C 140

148 150 160

* * * READ RANGE OF SVP READ(NIU, *)RSVP *** SVP BEYOND START RANGE? IF (RSVP .GT. RA) THEN NSVP = 0 GO TO 148 ENDIF *** NO. DETERMINE IF SVP IN RUNSTREAM OR SUPPLIED BY SUBROUTINE USVP. READ(NIU, *)KSVP IF (KSVP .EQ. O) THEN CALL SVP(NLYR, ZLYR, RHO, BETA, IXSVP, NSVP, ZSVP, CSVP) ELSE CALL USVP ENDIF * * * ERROR IN SVP? IF (NSVP .EQ. O) THEN *** YES. WRITE (NPU,160) RA FORMAT (1X,~SVP missing or input error. Range = ',F8.1,' m.') *** GO TO EXIT POINT GO TO 9 9 9 9 ENDIF *** SAVE POINTER TO LAST SVP VALUE IN FIRST LAYER. NWSVP = IXSVP(1)

*** IF CO NOT SPECIFIED, SET CO TO AVERAGE SPEED OF FIRST PROFILE *** IN FIRST LAYER AT INITIAL RANGE IF (DABS(CO) .LT. EPS) THEN DO 180 I = 2, NWSVP C0 = CO + (ZSVP(I) - ZSVP(I - 1)) * (CSVP(I - Z) I + 0.SDO * (CSVP(I) - CSVP(I - I))) 180 CONTINUE CO = C0 / zsvp(m4svp) ENDIF *** COMPUTE REFERENCE WAVE NUMBER XKO = 2.0D0 * PI * FRQ / C0 *** COMPUTE ATTENUATION - SACLANT MEMO SM-121 (JENSEN + FERLA) *** MODIFIED AS FOLLOWS: *** IF INPUTrED BETA IS LT 0.0, ALPHA IS COMPUTED IN DB/METER AND USED FOR BETA ALPHA = FI~Q * FRQ * 1 (7.0D-3 + (0.155D0 * 1.7D0) / 2 (1.7D0 * 1.7DO + FRQ * FRQ * I.OD-6)) * 1.0D-09 * * * ADJUST LAYER DEPTHS IN CASE BOTTOM SLOPES AND RA.NE.RSVP *** ASSUMES LAYERS ARE PARALLEL AND FOLLOW BOTTOM CONTOUR. * * * LAYER DEPTHS ARE ENTERED WITH SVP. DO 200 ILYR = 1, NLYR ZLYR(ILYR) = ZLYR(ILYR) + (RA - RSVP) * DTAN(THETA) 2 0 0 CONTINUE *** GET RANGE OF NEXT SVP READ(NIU, *, END = 210)RSVP GO TO 220 *** ONLY ONE SVP - SET RSVP LARGE SO SAME SVP USED FOR WHOLE PROBLEM 210 RSVP - 1 . 0 D 3 8 220 CONTINUE *** IF STARTING FIELD IS BEYOND NEXT SVP, GO BACK AND GET NEXT SVP IF (RSVP .LE. RA) GO TO 140 * * * I F ITYPEB = 2 0R 3 , AND ZA • O, EXTEND BOTTOM BY SETTING ZA TO * * * 4 / 3 MAX DEPTH

741

J.S. ROBERTSON,W. L. SIEGMANN, and M. J. JACOBSON

742

*** SEE NORDA TECH NOTE 12, JAN 78, H. K. BROCK *** IF ITYPEB - 2 OR 3 AND ZA NOT O, EXTEND BOTTOM TO ZA METERS IF (ITYPEB .EQ. 2 .OR. ITYPEB .EQ. 3) THEN C *** EXTEND BOTTOM TO ZA METERS IF ZA NOT ZERO C ~ * EXTEND BOTTOM 4/3 MAX DEPTH IF ZA = ZERO IF (DABS(ZA) .LT. EPS) ZA = 4.0DO * ZLYR(NLYR) / 3.0DO IF (ZA .LT. ZLYR(NLYR)) THEN C ~e~ USER ATTEMPTED TO EXTEND DEPTH IN NEGATIVE DIRECTION WRITE (NPU,230) 230 FORMAT(IX, 'Error. ZA reset to max depth of bottom layer. ') ZA - ZLYR(NLYR) C * * * INSERT PARAMETERS FOR EXTENDED BOTTOM IN APPROPRIATE C ARRAYS ENDIF NLYR = NLYR + 1 C *** STORE DEPTH OF ARTIFICIAL ABSORBING LAYER ZLYR(NLYR) = ZA C *** USE SAME DENSITY AND ATTENUATION AS LAYER ABOVE REOCELYe) - SeO(m..YS - 1) BETA(NLYR) - BETA(ELYB, - 1) C *** USE BOTTOM SP.EED OF LAYER ABOVE NSVP - NSVP + 1 IXSVP(ELYR) - NSVP ZSVPCNSVP) n ZLYR(NLY~)

C C

CSVP(NSVP)

C C C ~-

-

CSVP(NSVP

-

1)

ELSE *e,w I F BOTTOM NOT EXTENDED AND ZA=O, SET ZA TO MAX DEPTH INPUTTED ~*~ WITH FIRST SVP IF (DABS(ZA) .LT. EPS) ZA = ZLYR(NLYR) ENDIF

*** IF N NOT SPECIFIED, SET N SO THAT DZ .LE. 1/4 WAVELENGTH IF (N .EQ. 0) N = INT((4.0D0 * ZA t FRq / C0)) + 1 IF (N .GT. MXN) THEN WRITE (NPU,270) MXN FORMAT (' Error. N too large. N reset to ', I5, '.') 270 N=MXN ENDIF *** COMPUTE RECEIVER DEPTH INCP/24ENT - DZ MAY BE SUCH THAT RECEIVERS DO ~** NOT LIE EXACTLY ON LAYER INTERFACES DZ = ZA / DBLE(FLOAT(N))

C C C C

*** IF BOTTOM IS FLAT, RANGE STEP MAY BE SPECIFIED *** IF RANGE STEP NOT SPECIFIED, SET IT EQUAL TO i/2 WAVELENGTH IF (DABS(DR) .LT. EPS) THEN DR = O.BDO * CO / FRQ ENDIF *** IF BOTTOM IS NOT FLAT, COMPUTE RANGE STEP * * * IT MAY BE NECESSARY TO REDUCE DZ SUCH THAT DR IS SMALL ENOUGH *** TO PRODUCE ACCURATE RESULTS IF (ITYPEB.NE.3.AND.DABS(THETA).GT.EPS) DR = ABS(DZ/TAN(THETA)) ~e~ COMPUTE DEPTH INCREMENT FOR ADJUSTING LAYERS DZZ = DR * DTAN(THETA)

C C

e** GET STARTING FIELD I F ( I S F .Eq. 0) THEN CALL SFIELD ELSE CALL UFIELD ERDIF

IF (IHNK .NE. O) THEN *** DIVIDE STARTING FIELD BY HANKEL FUNCTION IF REQUESTED BY USER HNK = HNKL(XKO * BA) DO 290 I - I, N g(I) - U(I) / lINK CONTINUE 290 ENDIF

***

COMPUTE DEPTH WRITE INCREMENT TO NEAREST DZ

Predicting underwater sound propagation

743

IF (WDZ .LT. DZ) WDZ = DZ IWZ = WDZ / DZ + 0.SDO WDZ = IWZ * DZ *** COMPUTE DEPTH PRINT INCREMENT TO NEAREST DZ IPZ = PDZ / DZ + 0.5DO IF (PDZ .GT. 0.0 .AND. IPZ .EQ. O) IPZ = 1 PDZ = IPZ * DZ *** PRINT PROBLEM PARAMETERS WRITE (NPU,310) DATE, TIME 310 FORMAT (//IX, 'OS/2-1FDPE Solution - Double Precision Version ', 1

A16,

lX,

A16)

IF (ISF .EQ. O) THEN WRITE (NPU,320) FORMAT(1X, 'Gaussian Starting Field', /) ELSE WRITE (NPU,330) FORMAT (IX, 'User Starting Field', /) ENDIF IF (IHNK .NE. O) WRITE (NPU,340) FORMAT (IX, 'Starting filed divided by Hankel Function') IF (ITYPES .NE. O) WRITE (NPU,350) FORMAT (1X, 'User Surface Condition') IF (ITYPEB .EO. I) WRITE (NPU,360) FORMAT (IX, 'User Bottom Condition') WRITE (NPU,370) FP~, ZS, CO, RA, ZA, N FORMAT (IX, 'Frq = ', F7.1, ' Hz', /, 1X, 'ZS = ', FT.I,

320 330 340 350 360 370

1

' m',

2 3

'BA

/,

1X,

'CO

= ',

1;7.1,

' m/sec',

/,

1X,

= ', F7.1, ' m', /, lX, 'ZA = ', FT.I, ' m', /, 1X, 'N = ', I5) WRITE (NPU,380) DR, WDR, PDR, DZ, WDZ, PDZ, ITYPEB, ITYPES, RMAX 380 FORMAT(IX, 'DR = ', F7 1, ' m', /, lX, 'WDR = ', F7.1, 1 ' m', /, 1X, 'PDR = ', F7.1, ' m', /, lX, 'DZ - ', 2 F9.4, ' m', /, lX, 'WDZ = ', F9.4, ' m', /, IX, 3 'PDZ = ', F9.4, ' m', /, IX, 'ITYPEB ', II, /, IX, 4 'ITYPES = ', If, /, IX, 'RMAX -', F8.1, ' m', /, /o 5 2X, 'Layer Max Depth(m) Density(g/cm,,3) Art(dR/el)' 6

,

/)

DO 400 ILYR = 1, NLYR WRITE (NPU,390) ILYR, ZLYR(ILYR), RHO(ILYR), BETA(ILYR)

390 FORMAT(2X, I3, 3X, F13.6, 5X, F13.6, 5X, F13.6) 400 CONTINUE *** PRINT BOTTOM DEPTHS IF REQUESTED IF (IBOT .NE. O) THEN TH = THETA * DEG WRITE (NPU,610) RI, Zl, R2, Z2, TH ENDIF IF (ISFLD .NE. O) THEN

420

* * * PRINT STARTING FIELD WRITE (NPU,420) FORMAT (/, IX, 'Starting Field') DO 440 I = 1, N Zl = I * DZ

WRITE (NPU,430i I , ZI, U(I) 430 440

FORMAT (IX, I4, 3X, FIO.2, 3X, ' ( ' , CONTINUE ENDIF

E12,S, 2X, E12.5,

C C

IF (ISVP .NE. O) THEN *** PRINT SVP WRITE (NPU,460) BA 460 FORMAT (//IX, 'Sound Velocity Profile at Range ', F9.2, I ' meters',/) DD 480 I ffi I, NSVP WRITE (NPU,470) I, ZSVP(I), CSVP(I) 470 FORMAT (1X, I4, 3X, F8.1, 3X, F8.1) 480 CONTINUE ENDIF **~ COMPUTE 'B' Coefficient

BCOF = DCMPLX(0.0D0, 0.5D0 / XK0) C

' )')

744

J.S. ROSER'rsoN, W. L. SmOMA~, and M. J. JACOBSON *** ASSIGN OUTPUT FILE IF OUTPUT REQUESTED IF (WDR .GT. ZERO) THEN WRITE (NOU,510) (U(1), I= IWZ, N, IWZ) FORMAT (' ~, F12.10, IX, F12.10) 510 ENDIF *** INITIALIZE RANGE VARIABLE AT WHICH SOLUTION IS TO BE PRINTED XPR = RA + PDR *** INITIALIZE RANGE VARIABLE AT WHICH SOLUTION IS TO BE RECORDED XWR = RA ÷ ',/DR IF (DABS(WDR) .LT. EPS) XWR = RA + RMAX ÷ 1.0

*** SAVE RANGE STEP OLDR = DR *** INITIALIZE PARAMS FOR DIAG THEN COMPUTE MAIN DIAGONALS X AND Y N1 = 1 DR1 = DR *** COMPUTE X DIAGONAL *** COMPUTE Y DIAGONAL CALL DIAG CALL DIAG MAIN LOOP STARTS HERE * * * < .......... SOLUTION WILL BE ADVANCED FROM RANGE RA TO RANGE RA+DR *** IMPLEMENTED AS A WHILE LOOP 530 IF (RA .LT. RMAX) THEN IDIAG = 0 IF (ITYPEB .F~. 3) GO TO 540 IF (DABS(THETA) .GT. EPS) GO TO 550 540 IF (DABS(DR - OLDR) .GT. EPS) IDIAG = I DR = OLDR 550 CONTINUE *** ***

560

*** DOES SVP CHANGE BEFORE BOTTOM PROFILE? IF (RSVP .GE. R2) GO TO 560 *** YES. DOES SVP CHANGE BEFORE NEXT SOLUTION RANGE? IF (RA + DR .LE. RSVP) GO TO 650 *** YES. ADJUST DR SO THAT SOLUTION ADVANCES TO RSVP DR = R S V P - RA GO TO 620 *** BOTTOM PROFILE CHANGES BEFORE OR AT SAME RANGE AS SVP CONTINUE *** DETERMINE IF BOTTOM DEPTHS ARE TO BE UPDATED *** FIRST PASS - RA WILL BE .LT. R2 *** UPDATE BOTTOM DEPTHS? IF (RA .GE. R2) GO TO 570 *** NO. *** DETERMINE IF RANGE STEP TOO LARGE IF (RA + DR .LE. R2) GO TO 650 *** RANGE STEP IS TOO LARGE - RESET DR - ADV SOLUTION TO R2 DR = R2 - RA GO TO 620

570

*** UPDATE BOTTOM DEPTHS CONTINUE RI = R2 Z1 = Z2 ITRK = ITRK + 1 R2 = TRACK(ITRK,1) Z2 = TRACK(ITRK,2) *** TWO DEPTHS AT SAME RANGE INDICATE VERTICAL DISCONTINUITY. *** ADVANCE TRACK FORWARD. IF (DABS(R1 - R2) .LT. EPS) GO TO 570 *** RESTORE DR DR = OLDR *** COMPUTE SLOPE OF BOTTOM THETA = DATAN2(Z2 - ZI, R2 - RI)

C C C

*** IF BOTTOM OF _ARTIFICIAL LAYER IS FLAT, DO NOT RECOMPUTE DR.

Predicting underwater sound propagation *** IF BOTTOM IS NOT FLAT, COMPUTE NEW RANGE STEP IF (DABS(THETA) .GT. EPS .AND. ITYPEB .WE. 3) THEN DR = DABS(DZ / DTAN(THETA)) ENDIF *** IF RANGE STEP TOO LARGE, PRINT WARNING MESSAGE. RECOMPUTE DR. IF (RA + DR .GT. R2) THEN WRITE (NPU,S90) FORMAT (1X, ' P ~ e step too large for bottom irregularities') DR = R2 - RA ENDIF IF (RA + DR .GT. RSVP) THEN DR = RSVP - RA ENDIF

590

610 1 2 3 620

630

*** PRINT BOTTOM DEPTHS IF (IBOT .WE. O) THEN TH m THETA * DEG WRITE (MPU,BI0) RI, Zl, R2, Z2, TH FORMAT (//1X, 'Bottom Depths ', //, IX, 'RI = ', F8.1, ' m ' , /,iX, 'Z1 ,, ', F8.1, ' M', /, IX, 'R2 = ', F8.1, ' m',/, IX, 'Z2 = ', F8.1, ' m ' , / , 1X, 'Theta = ' , F 8 . 1 , ' D e g ' ) ENDIF CONTINUE DZZ = DR * DTAN(THETA) WRITE (NPU,630) DR, RA FORMAT (iX, 'Range step = ', F7.2, ' m at range ', F9.2, ' m') *** ADVANCE RANGE ONE STEP RA = RA + DR IDIAG = 1 GO TO 660 CONTINUE RA = RA + DR CONTINUE *** READ NEW SVP PROFILE FLAG? IF(R~.GE.RSVP) BY~AD(NIU, *)KSVP *** I F KSVP NOT O, SUBROUTINE USVP IS CALLED FOR SVP. USER IS *** RESPONSIBLE FOR ADJUSTING DEPTHS OF LAYERS WHEN BOTTOM *** IS NOT FLAT. IF (KSVP .NE. o) THEN CALL USVP IDIAG -- 1 ELSE

650 660

*** NEW SVP AT ADVANCED RANGE? IF (RA .GE. RSVP) GO TO 720 *** NO. IS BOTTOM FLAT?. IF (DABS(THETA) .LT. EPS) GO TO 790 MLYR = NLYR

IF (ITYPEB .EQ. 3) MLYR = N L Y R - I *** UPDATE DEPTHS OF LAYERS * * * ASSUMES LAYERS ARE PARALLEL AND FOLLOW BOTTOM CONTOUR *** IF ITYPEB = 3 , BOTTOM OF ARTIFICIAL LAYER REMAINS FLAT. DO 690 ILYR = 1, MLYR ZLYR(ILYR) = ZLYR(ILYR) + DZZ CONTINUE N1 = ZLYR(1) / DZ NI=I ZA = ZLYR(NLYR) I F (ITYPEB .EO. S .AND. ZLYR(NLYR) .LT. ZLYR(MLYR)) WRITE

690 C

I 700 I C C C

(NPU,700)

FORMAT (1X, 'Err. Depth of artificial layer less than', ' layer above') *** ADJUST D E I ~ S OF PROFILES IN SLOPING LAYERS *** ASSUMES SAME SVP IN LAYERS NNSVP = IXSVP(1) NNWSVP = NWSVP + 1 DO 710 I " NNNSVP, NSVP ZSVP(I) = ZSVP(I) + DZZ

745

746

J.S. ROBERTSON,W. L. SIEGMANN,and M. J. JACOBSON 710

CONTINUE GO TO 800

C C

*** GET NEXT SVP CONTINUE CALL SVP(NLYR, ZLYR, RHO, BETA, IXSVP, NSVP, ZSVP, CSVP) IDIAG = I ENDIF * * * ERROR DETECTED? IF (NSVP .EQ. O) THEN WRITE (NPU, 160)RA GO TO 9999 ENDIF IF (ITYPEB .NE. S) ZA = ZA + DZZ *** NO. READ RANGE OF NEXT PROFILE? IF(RA.GE.RSVP) THEN READ(NIU, *, END = 7S4)RSVP GO TO 736 RSVP = 1.0D38 CONTINUE ENDIF *** ARTIFICIAL ABSORBING LAYER? IF (ITYPEB .EQ. 2 .OR. ITYPEB .EQ. S) THEN * * * YES. UPDATE DENSITY, ATTEN AND SPEED. IF (ZA .GE. ZLYR(NLYR)) THEN NLYR = NLYR + 1 ZLYR(NLYR) = ZA RHO(NLYR) = RHO(NLYR - 1) BETA(NLYR) -- BETA(NLYR - 1) NSVP = NSVP + I IXSVP(NLYR) = NSVP ZSVP(NSVP) = ZLYR(NLYR) CSVP(NSVP) = C S V P ( N S V P - i) ENDIF ENDIF

720

C

C

734 736 C C

C C

*** PRINT SVP IF (ISVP .NE. 01 THEN WRITE (NPU,460) PA DO 770 I = I, NSVP WRITE (NPU,470) I, ZSVP(I), CSVP(I) CONTINUE ENDIF

770 C

*** IF NEW DR AND/OR SVP - U P D A T E IF (IDIAG .EQ. O) GO TO 810

C

790

X AND

Y

DIAGONALS

NI=I

800

CALL DIAG DR1 = DR CONTINUE

810 C C

*** ADVANCE SOLUTION ONE STEP FORWARD NOLD = N CALL CP,NK

C C

*** APPLY ABSORPTION IF ITYPEB = 2 OR 3 IF (ITYPEB .EQ. 2 .OR. ITYPEB .EQ. 3) THEN MM = ZLYR(NLYR - 1) / DZ NA=N-MM IF (NA .LT. O) THEN WRITE (NPU,820) PA FORMAT (1X, 'Err in thickness of absorbing layer at ', 820 F8.1, ' m') 1 GO TO EXIT POINT C GO TO 9 9 9 9 ENDIF

C

1

840

IF (NA .GT. 01 THEN *** SEE AESD PE MODEL BY BROCK - NORDA TECH NOTE 12 - JAN 78 DO 840 I = 1 , NA ATT = DEXP(-.OIDO * DR * DEXP(-DBLE(FLOAT(I - NA) / (FLOAT(NA) / 3 . 0 D 0 ) ) * * 2 ) ) LOOP1 = MM + I U(LOOP1) -- U(LOOP1) * A F T CONTINUE ENDIF

Predicting underwater sound propagation

747

ENDIF

C C

***

860

C C 880

C C

TERMINATE BUN IF N HAS BEEN DECREMENTED TO 5 - ARBITRARY IF (N .LE. 5) THEN WRITE (NPU,860) FORMAT (IX, 'Program terminated - only five points remain') GO TD 980 ENDIF *** TERMINATE RUN IF N HAS, BEEN INCREMENTED BEYOND MAXIMUM IF (N .GT. MXN) THEN WRITE (NPU,880) BA FORMAT (IX, 'Max N exceeded a t range ', FI0.2, ' m.') GO TO 980 ENDIF *** DETERMINE IF NEW POINT ADDED IF (N .GT. NOLD) THEN N1 = NOLO CALL DIAG DR1 = DR ENDIF

C C C

*** IF SOLUTION IS TO BE WRITTEN ON DISK, *** WRITE SELECTED PRESSURE FIELD AT RANGE KA IF (XWR .LE. BA .AND. WDR .GT. ZERO) THEN WHITE (NOU, 5101 (U(I), I= IWZ, N, IWZ)

C C C

*** DETEBMINE NEXT RANGE AT WHICH TO WHITE SOLUTION ON DISK XWH = X W R + WDR ENDIF

C C

*** DErERMINE IF SOLUTION IS TO BE PRINTED IF (XPR .LE. RA .AND. IPZ .NE. 0 .AND. DABS(PDR)

C C

.GT. EPS) THEN

*** COMPUTE HANKEL FUNCTION IF (IHNK .NE. 01 HNK = HNKL(XKO * KA)

C C

*** COMPUTE AND PRINT PROPAGATION LOSS AT EACH IPZ'TH DEPTH DO 960 I = IPZ, N, IPZ ZI = I * DZ IF (IHNK .EQ. 11 THEN PL = CDABS(U(I) *. HNK) ELSE PL = CDABS(U(I)) ENDIF IF (PL .LE. O.ODO) THEN PL = 999.9D0 ELSE PL = -20.ODO * DLOGIO(PL) IF (IHNK .EQ. 01 PL = PL + IO.ODO * DLOGIO(RA) ENDIF

C WRITE(NQU, 9591 -PL FORMAT(F7.21 CONTINUE

959 960 C C

C C

*** DETERMINE NEXT RANGE AT WHICH TO PRINT SOLUTION XPR = XPR + PDR ENDIF *** THIS IS THE BOTTOM OF THE MAIN LOOP (WHILE) GO TO 5 3 0 ENDIF

C

C

< ..........

***

TERMINATE RUN CALL GETTIM(ITIME(1), ITIME(2), ITIME(31, ITIME(4)) CPUTIM - 3600.0D0 * ITIME(1) + 60.DO * ITIME(2) + I ITIME(3) - CPUTIM 980 WHITE (NPU,990) CPUTIM 990 FORMAT (/, IX, 'End of run. OS2IFD terminated after ', 1 F12.1, ' CPU s e c o n d s . ' ) C THIS IS THE EXIT POINT 9999 CONTINUE

CAOEO 1716-,-B

748

J.S. ROSERTSON, W. L. SmGMANN, and M. J. JACOSSON CLOSE (NIU) ENDFILE(NPU) CLOSE (NPU) IF((WDR .GE. 0.OD0 .AND. WDR .LE. 0.0D0) .OR. I (WDZ .GE. 0.OD0 .AND. WDZ .LE. 0.0D0)) %'HEN STWORD = 'DELETE' ELSE STWORD = 'KEEP ' ENDIF ENDFILE(NOU) CLOSE (NOU, STATUS.--STWORD) ENDFILE(NQU) CLOSE (NQU) STOP ' OS2IFD Normal termination.'

END SUBROUTINE CRNK IMPLICIT REAL*8 (A - H, 0 - Z), INTEGER*2 (I - M) C C

******$**************************************$**********$******$**

****************************************************************** * THIS ROUTINE USES THE CRANK-NICOLSON METHOD FOR SOLVING THE *

C C

* SYSTEM OF EQUATIONS. THE SOLUTION IS ADVANCED FROM RANGE RA-DR *

C

* TO RANGE RA WHERE RA IS THE ADVANCED RANGE.

C C C

*** SYSTEM OF EQUATIONS IS AS FOLLOWS:

C

* X(1)

C

*

C

*

C

***

*

ALL VALUES ARE AT ADVANCED RANGE -R12(1)

-1 0

X(2) -i

0 -R12(2) X(N)

*

* U(1)

*

* *

* U(2) * U(N)

* *

***

***

* SURX * -

***

* 0 * BOTX

* *

***

***

***

***

=

C ALL VALUES ARE AT PRESENT RANGE

C C

***

C

* Y(1) * ÷1 * 0

C C C C

***

+R12(1)

Y(2)

0

*

+R12(2) *

+1

Y(N)

*

***

***

* U(1) * * U(2) * * U(N) *

* SURY * +

*

0

*

* BOTY *

* $$$$$$*$*$*$$$*$$$$***$$****$$$$$$$$$$$$$*$$**$$$*$$*$$*$*$*$***$$

C

C

D

- ARRAY - RIGHT HAND SIDE

C

RA

- RANGE TO WHICH SOLUTION IS TO BE ADVANCED - M ~ I ~ S

C C

TA

C C

C C C C C C C

C C

TB

- RA IS INCREMEWrED BY DR PRIOR TO ENTERING THIS ROUTINE - TEMPORARY VARIABLE USED IN MANIPULATION OF BOTTOM POINT. - TEMPORARY VARIABLE USED IN MANIPULATION OF BOTTOM POINT.

*** SEE MAIN PROGRAM FOR OTHER DEFINITIONS *** UNDEFINED VARIABLES ARE TEMPORARY VARIABLES *** SUBROUTINE CRNK RETURNS: * N - NUMBER OF EQUI-SPACED POINTS IN U AT RANGE RA * U - ARRAY - COMPLEX ACOUSTIC PRESSURE FIELD AT RANGE RA - INCLUDES BOTTOM POINT - DOES NOT INCLUDE SURFACE POINT *** SUBROUTINE TRID IS CALLED TO SOLVE THE TRIDIAGONAL MATRIX

C

$INCLUDE:'COMMON.IFD' DIMENSION D(MXN) COMMON ICRNKER/ D COMPLEX*16 D, SUM, A, EX, P, Q, XX, SI, S2, CA, CB, CC, CD, TA, TB PARAMETER (EPS = I.OD-IO) C C *** CALL SCON FOR SURFACE CONDITION CALL SCON C C *** COMPUTE RIGHT HAND SIDE D(1) THROUGH D(N-1) C *** D(N) IS COMPUTED LATER D(1)

= Y(1)

* U(1)

+ R12(1)

* U(2)

+ SURY + SURX

DO I0 1 = 2, N - 1 10 D(I) = U(I - I) ÷ Y(I) * U(I) + RI2(I) * U(I + 1)

Predicting underwater sound propagation

* * * BOTTOM TYPE IB = ITYPEB + 1 GO TO ( 2 0 , 5 0 , 9 0 ,

I00),

749

IB

*** BOTTOM IS RIGID - ITYPEB = 0 20 CONTINUE IF (THETA .GT. 0.0) GO TO 40 IF (DABS(THETA) .LT. EPS) THEN *** RIGID BOTTOM IS FLAT BOTY = U(N) D(N) = U(N - I ) + Y(N) * U(N) + BOTY TA = (O.ODO, O.ODO) *** BOTX=U(N) AT ADVANCED RANGE TH = (I.ODO, O.ODO) CALL TRID(N, N, TA, TB) RETURN

*** RIGID BOTTOM SLOPES UPWARD - DELETE A POINT - USE 2ND ORDER ODE. ENDIF BOTY = U(N) N=N-1 D(N) = U(N - I) + Y(N) * U(N) + BOTY COT = DCOS(THETA) / DSIN(THETA) P = -COT / BCOF q = (ACOFX + DCMPLX(O.ODO, XKO)) / BCOF XX = C D S O R T ( P * P - 4.0DO * {~) $1 = (-P + XX) 1 2.0DO $2 = (-P - XX) / 2.0DO CA = (I.0DO - $I * DZ) / (-XX * DZ) CB = 1 . 0 D O / (XX * DZ) CC = 1.0DO - CA CD = CB TA = -CD * CDEXP(SI * DZ) + CB * CDEXP(S2 * DZ) TB = CC * CDEXP(S1 * DZ) + CA * CDEXP(S2 * DZ) CALL TRID(N, N, TA, TH) RETURN

40 CONTINUE *** RIGID BOTTOM SLOPES DOWNWARD - USE 2ND ORDER O.D.E. COT = DCOS(THETA) / DSIN(THETA) P = -COT / BCOF = (ACOFY + DCMPLX(O.ODO,XKO)) / BCOF XX = CDSQRT(P * P - 4 . 0 D O * Q) $1 = (-P + XX) / 2.0DO $2 = (-P - XX) / 2.0DO CA = (I.0DO - $1 * DZ) I (-XX * DZ) CB = 1 . O D O / (XX * DZ) CC = 1.0DO - CA CD=CB TA = -CD * CDEXP(Sl * DZ) + CB * CDEXP(S2 * DZ) TB = C( * CDEXP(S1 * DZ) + CA * CDEXP(S2 * DZ) C C C

MODIFY LOWER AND Y DIAGONALS IN ROW N BY TA AND TH COMPUTE D(N) D(N) -- (1.0DO + TA) * U(N - 1) ÷ (Y(N) ÷ TB) * U(N) *** ***

*** COMPUTE TA AND TB FOR LOWER AND X DIAGONALS IN R0N N O = (ACOFX + DCMPLX(0.0DO,XK0)) / BCOF XX = CDS@RT(P * P - 4.0DO * Q) $I = (-P + XX) / 2.0DO $2 = ( - P - XX) / 2.0DO CA = (I.0DO - $1 * DZ) / ('XX * DZ) CB = 1 . 0 D O / ( X X * DZ) CC = I . O D O - CA CD=CB TA = -CD * CDEXP(SI * DZ) + CB * CDEXP(S2 * DZ) TB = CC * CDEXP(SI * DZ) ÷ CA * CDEXP(S2 * DZ) CALL TRID(N, N, TA, TB)

ADD A POINT CB = ( ( U ( N ) - UCN - 1 ) ) / D Z CA = U ( N ) - CB ***

- $1 * UCN)) /

(-XX)

- THEN

750

J.S. ROBERTSON, W. L. SIEOMANN, and M. J. JACOBSON N=N+ I U(N) = CA * CDEXP(SI*DZ) RETURN

+ CB * CDEXP(S2

* DZ)

50 CONTINUE *** BOTTOM CONDITION SUPPLIED BY USER - ITYPEB = I IF (DABS(THETA) .GT. EPS) GO TO 80 IF (DABS(THETA) .GT. EPS) GO TO 80 *** FLAT BOTTOM *** USER SUPPLIES BOTY ON BOTTOM INTERFACE AT PRESENT RANGE *** USER SUPPLIES BOTX ON BOTTOM INTERFACE AT ADVANCED RANGE 60 CALL BCON TA = 0.ODO TB= O.ODO N=N-1 D(N) = U(N - I) + Y(N) * U(N) + BOTY + BOTX CALL TRID(N, N, TA, TB) N=N+I U(N) = BOTX RETURN 70 CONTINUE *** LAYER SLOPES DOWN *** USER SUPPLIES ; BOTY - ONE DZ IN BOTTOM AT PRESENT RANGE BOTX - ON BOTTOM INTERFACE AT ADVANCED RANGE CALL BCON TA = O.ODO T B = O.ODO D(N) = U(N - 1) + Y(N) * U(N) + BOTY + BOTX CALL TRID(N, N, TA, TB) N - - N + I

U(N)

=

BOTX

RETURN 80 CONTINUE *** LAYER SLOPES UP *** USER. SUPPLIES : BOTX - ONE DZ BELOW BOTTOM INTERFACE AT ADVANCED RANGE CALL BCON N=N-1 D(N) = U(N - 1) + Y(N) * U(N) + BOTY + BOTX TA = O.ODO TB = O.ODO CALL TRID(N, N, TA, TB) RETURN C C C

*** ITYPEB = 2 -- ARTIFICIAL ABSORBING 90 CONTINUE IF (DABS(THETA) .GT. EPS) GO TO 110

LAYER INTRODUCED

*** BOTTOM OF ABSORBING LAYER IS FLAT *** ITYPEB = 2 OR 3 -- ARTIFICIAL ABSORBING i00 U(N) = O.ODO TA = O.ODO TB = O.ODO CALL TRID(N - 1, N, TA, TB) RETURN

110

*** ITYPEB = 2 -- ARTIFICIAL ABSORBING IF (THETA .GT. 0.0D0) GO TO 120

LAYER NOT FLAT

*** BOTTOM OF ABSORBING LAYER SLOPES UP TA = 0.0D0 TB = 0.OD0 U(N - 1 ) = 0.OD0 U(N) = 0.0D0 CALL T R I D ( N - 2, N - I, TA, TB) *** DELETE A POINT N=N-I RETURN

LAYER

Predicting underwater sound propagation

751

120 CONTINUE *** BOTTOM OF ABSORBING LAYER SLOPES DOWN D(N) = U(N - 1) TA = 0.0D0 TB = O.ODO CALL TRIO(N, N, TA, TB) *** ADD A POINT N=N+I U(N)

=

O.ODO

RETURN END SUBROUTINE TRID(L, M, TA, TB) IMPLICIT REAL*8 (A - H, 0 - Z), INTEGER*2

(I - M)

*** SPECIALIZED VERSION OF METHOD PRESENTED ON PG 442 OF CARNAHAN *** ET AL FOR SOLVING A TRIDIAGONAL MATRIX. *** TRIO RETURNS THE SOLUTION FIELD IN U *** LOWER DIAGONAL = -I *** UPPER DIAGONAL = -RI2 *** BTA IS PARTIAL SOLUTION COMPUTED IN ROUTINE DIAG *** D CONTAINS R.H.S. *** GAMMA - TEMPORARY STORAGE *** TA - TEMPORARY VARIABLE USED IN MANIPULATION OF BOTTOM *** POINT *** TB - TEMPORARY VARIABLE USED IN MANIPULATION OF BOTTOM *** POINT ****************************************************************** $INCLUDE : 'COMMON. IFD ' COMMON /TRIDI/ GAMMA COMMON /CRNKER/ D DIMENSION D(MXN), GAMMA(MXN) COMPLEX*16 GAMMA, D, TA, TB C *** X AND LOWER DIAGONAL IN ROW L MUST BE MODIFIED BY TB AND TA BTA(L) = (X(L) - TB) - ((I.0D0 + TA)*RI2(L - i)) / BTA(L - I) GAMMA(1) = D(1) / BTA(1) DO 10 I = 2, L GAMMA(I) = (D(I) + GAMMA(I - I)) / BTA(I) 10 CONTINUE GAMMA(L) = G A ~ ( L ) + (TA * GAMMA(L - I)) / BTA(L) C *** IF L IS LESS THAN M, U(M) IS KNOWN. IF (L .EQ. M) U(M) = GAMMA(M) ML00P = M - 1 DO 20 I = I, MI..OOP MM = MLOOP I + 1 U(MM) = GAMMA(MM) ÷ RI2(MM) * U(MM + I) / BTA(MM) 20 CONTINUE RETURN END SUBROUTINE DIAG IMPLICIT REAL*8 ( A - H, 0 - Z), INTEGER*2 (I - M) C ***********~***~************~*************~*********************** C ****************************************************************** C * COMPUTES RANGE AND DEPTH DEPENDENT MAIN DIAGONALS OF MATRICES C * DIAG IS CALLED WHENEVER BOTTOM DEPTH OR SVP CHANGES C ****************************************************************** C ****************************************************************** C C *** SEE MAIN PROGRAM FOR DEFINITIONS C *** SUBROUTINE DIAG RETURNS: C ACOFX - Coefficient 'A' AT BOTTOM - AT RANGE RA C ACOFY - Coefficient 'A ~ AT BOTTOM - AT RANGE RA-DR C BTA - ARRAY - PARTIAL SOLUTION OF SYSTEM OF EQUATIONS C X - ARRAY - MAIN DIAGONAL OF MATRIX AT RANGE RA C Y - ARRAY - MAIN DIAGONAL OF MATRIX AT RANGE RA-DR C ****************************************************************** C $INCLUDE : ' COMMON. IFD ' COMPLEX*16 B, XN1, XN2, EYE PARAMETER (EYE = (O.ODO, 1.0DO)) PARAMETER (EPS = I.OD-IO) INTEGER IF/ C

• *

752

J.S. ROmmTSON, W. L. S1EGMANN, and M. J. JACOBSON *** COMPUTE NEW X AND Y DIAGONALS

C

C C

*** GET INDICES, DENSITY, AND ATTENUATION FOR FIRST LAYER IEX - 0 ILYR - 1 L=I M = IXSVP(1) m = ~o(1) BETA1 = BETA(1) DO I0 I = NI, N *** TRANSFORM X INTO Y Y(I) = -X(I) - EYE * 2.0DO * DZ * DZ * XKO * 1 (1.0DO + R12(1)) * (I.ODO / DRI + 1.0DO / DR) 10 CONTINUE DO 120 I = NI, N *** ZI IS RECEIVER DEPTH ZI = I * DZ *** IS RECEIVER IN THIS LAYER? 20 IF (ZI .LE. ZLYR(ILYR)) GO TO 30

C

* * * NO. ALL LAYERS CHECKED?

C

***

IF (ILYR .EQ. NLYR) GO TO 50 NO. SET INDICES, DENSITY, AND ATTENUATION FOR NEXT LAYER. ILYR = ILYR ÷ 1 L=M+I M = IXSVP(ILTR) R1 ,, RH0(ILYR) BETA1 = B E T A ( I L Y R )

C C

C

C

GO TO 20 *** DEPTH ZI IS IN THIS LAYER. 30 CONTINUE *** DETERMINE SOUND SPEED CI AT DEPTH ZI DO 40 J = L, M IF (ZI .GT. ZSVP(J)) GO TO 40 INTERPOLATE CI = CSVP(J - 1) + ( C S V P ( J ) 1 ) / ( Z S V P ( J ) - ZSVP(J - 1 ) ) GO TO 90 40 CONTINUE L=J * ** EXTRAPOLATE 50 CONTINUE ***

IF (DABS(ZSVP(N) - ZSVP(N-1)) CI = CSVP(M) GO TO 70 CI = CSVP(M -

60 i

70 80 90

- CSVP(J -

1))

* (ZI -

ZSVP(J - 1)

.GT. EPS) GO TO 60

1) + (CSVP(M) - CSVP(M -

1))

* ( Z I - ZSVP(M -

/ (ZSVP(H) - ZSVP(M - I)) IF (IEX .EQ. 1) GO TO 90 WRITE (NPU,80) FORMAT (lX, 'Warning. Extrapolation of SVP performed. ') IEX=I L=J

*** SAVE SPEED IN MEDIUM 1 Cl = CI C *** IS RECEIVER ON OR WITHIN 1 DZ OF INTERFACE? IF (ZLYR(ILYR) - ZI .GE. DZ) GO TO 100 C *** YES. IS THIS BOTTOM LAYER? IF (ILYR .EQ. NLYR) GO TO 100 C *** NO. GET PARAMETERS FOR MEDIUM 2 R2 = RHO(ILYR ÷ 1) BETA2 = BETA(ILYR + 1) C2 - CSVP(IXSVP(ILYR) ÷ I) GO TO 110 100 CONTINUE C * * * NOT AN INTERFACE. USE MEDIUM 1 PARAMETERS. C2 - Cl R2 = R1 C

BETA2 ,, BETA1 * * * COMPUTE DENSITY RATIO 110 R 1 2 ( I ) = RI l R2 C * * * THE NEXT FEW LINES COMPUTE THE X DIAGONAL

C

XN = CO / C1

IF (BETAI .LT. O.ODO) BETA1 = ALPHA * C1 / FRO XN = XN * XN XNI = XN * BETAI / 27.28752TD0

XN1

=

DCMPLX(XN,XNI)

1))

Predicting underwater sound propagation

753

XN = CO / C2 IF (BETA2 .LT. O.ODO) BETA2 ffi ALPHA * C2 / FRO XN ffi X N * XN XNI ," XN * BETA2 / 27.287527 XN2 = DCMPLX(XN,XNI) B = (-XK0 * 0.5D0 * ((XNI - 1.0D0) + Rl2(I) * (XN2 - 1.0D0))) XN = 2.0D0 * (I.0D0 + RI2(1)) / DR X(I) = DZ * DZ * XK0 * (B - DCMPLX(0.0D0,XN)) + 1.0D0 + RI2(I) 120 CONTINUE *** COMPUTE BTA(NI) THROUGH BTA(N) * * * ARRAY BTA CONTAINS PARTIAL SOLUTION OF SYSTEM OF EOUATIONS IF (NI .Eq. 1) G0 TO 130 M = NI GO TO 140 130 BTA(1) = X(1) M

=

2

140 DO 150 1 = M, N LOOP1 = I - I B T A ( I ) ffi X ( I ) - R12(LOOP1) / BTA(LOOP1) 150 CONTINUE *** COMPUTE 'A' Coefficients AT BOTTOM ACOFY ffi ACOFX ACOFX = DCMPLX(O.ODO, O.SDO) * XKO * (XN2 - 1.0DO) RETURN END SUBROUTINE SCON IMPLICIT REAL*8 (A - H, D - Z), INTEGER*2 (I - M) ***************************************************************** *** SURFACE CONDITION SUBROUTINE *** I F ITYPES ffi O, SCDN SETS SURY AND SURX = 0 . 0 . *** I F ITYPES NOT 0 , THE USER MUST SUPPLY SURY AND SURX. *** SEE MAIN PROGRAM FOR DEFINITIONS

$INCLUDE:'COMMON.IFD' IF (ITYPES .Zq. O) THEN C C * * * PRESSURE RELEASE SURFACE SURY = DCMPLX(O.ODO, O.ODO) SURX ffi DCMPLX(O.ODO, O.ODO) * * * USER SURFACE CONDITION ELSE *** USER ADDS SURFACE CONDITION HERE ENDIF RETURN END COMPLEX*16 FUNCTION HNKL(X) IMPLICIT REAL*8 (A - H, 0 - Z), INTEGER*2 * *

HANKEL FUNCTION HO(1) - POLYNOMIAL APPROXIMATION HANDBOOK OF MATH FUNCTIONS - N.B.S. - NOV 1967

REAL*8 JO PARAMETER(PI f f i 3 . 1 4 1 5 9 2 6 5 3 5 8 9 7 9 3 D O ) PARAMETER ( * A1 " -2.2499997D0, * A2 = 1.2656208D0, * A 3 = -0.3163866D0, * A4 = 0 . 0 ! ~ 7 9 D 0 , * A5 = - 0 . 0 0 3 9 4 4 4 D 0 , * A 6 = 0.0002100D0

*

(I - M)

) PARAMETER ( * B1 ffi 0.36746691D0, * B 2 = 0.60559366D0, * B3 = -0.7435038D0, * B4 = 0.25300117D0, * B5 = -0.04261214D0, * B6 = 0.00427916D0, * B 7 ffi -0.00024846D0 * )

* *

754

J . S . ROBERTSON, W. L. SIEGMANN, a n d M . J. JACOBSON

PARAMETER * * * * * * *

Cl C2 C3 C4 C5 C6 C7

= = = = = = =

PARAMETER D1 = D2 = D3 = D4 = D5 = D6 = D7 =

( 0.79788456D0, -0.00000077D0, -0.00552740D0, -0.00009512D0, 0.00137237D0, -0.00072805D0, 0.00014476D0 ( 0.78539816D0, -0.04166397D0, -0.00003954D0,

0.00262573D0, -0.00054125D0, -0.00029333D0, 0.00013558D0

IF (X .LE. 3 . 0 D O )

THEN

Y = X * X / 9.0DO JO = I . O D O + Y * (AI + Y * (A2 + Y * (A3 + Y * (A4 + Y * (A5 + Y * ( A S ) ) ) ) ) ) YO = 2.0DO * DLOG(O.SDO * X) * JO / PI + B1 + Y * (B2 1 2

Y

*

(B3 + Y *

Y * (BT)))))) HNKL = DCMPLX(JO, ELSE

(B4 + Y *

(B5 + Y *

(B6 +

YO)

Y = 3.0DO / X FO = Cl + Y * (C2 + Y * (C3 + Y * (C4 + I Y * (C5 + Y * (C6 + Y * (C7)))))) TO -- X - D1 + Y * (D2 ÷ Y * (D3 + Y * (D4 + 1 Y * (D5 + Y * (D6 + Y * (D7)))))) HNKL -- FO * CDEXP(DCMPLX(O.ODO, TO)) / DSQRT(X) ENDIF RETURN END SUBROUTINE SVP(NLYR, ZLYR, RHO, BETA, IXSVP, NSVP, ZSVP, CSVP) IMPLICIT REAL*8 (A - H, 0 - Z), INTEGER*2 (I - M) ****************************************************************** • ** SOUND VELOCITY PROFILE SUBROUTINE • ** CALLING P R O G R A M SUPPLIES: NOTHING • ** SVP SUBROUTINE RETURNS: • N L Y R - NUMBER OF LAYERS. LAYER 1 IS WATER. OTHERS ARE SEDIMENT • ZLYR - ARRAY - DEPTH OF EACH LAYER. FIRST IS DEPTH OF WATER. • RHO - ARRAY - DENSITY OF EACH LAYER. GRAMS/CUBIC CM • BETA - ARRAY - ATTENUATION IN EACH LAYER. DE/WAVELENGTH • IXSVP- ARRAY - CONTAINS POINTERS. POINTS TO LAST VALUE OF SVP • IN CORRESPONDING LAYER. SVP IS STORED IN A R ~ Y S ZSVP • AND CSVP. IXSVP(1) POINTS TO LAST SVP P O I N T IN WAT?_A. • IXSVP(NLYR) POINTS TO LAST SVP POINT IN BOTTOM-MOST • LAYER. • N S V P - NUMBER OF POINTS IN ZSVP AND CSVP. ZSVP AND CSVP • CONTAIN THE PROFILES FOR ALL LAYERS. • ZSVP - ARRAY - SVP DEPTHS - METERS • CSVP - ARRAY - SOUND SPEED - METERS/SEC

PARAMETER PARANETER PARAMETER DIMENSION 1

( N I U = 7 , NOU = 8 , NOU = 9 , NPU = 1 0 ) (NXN = 1024, MXLYR = 64, MXTRK = MXLYR, MXSVP = 64) (MXDM -- MXLYR) ZLYR(bIXDM), RHO(MXDM), BETA(MXDM), IXSVP(MXDM), ZSVP(MXDM), CSVP(MXDM)

NSVP = 0 *** READ N U M B E R OF LAYERS READ(NIU, *, END = 9999)NLYR *** F I R S T LAYER IS WATER. OTHERS ARE SEDIMENT. DO 20 I = I, N L Y R *** READ DEPTH OF LAYER, DENSITY AND ATTF.A'UATION. READ(NIU,*) ZLYR(I) ,PJ{O(I) ,BETA(I)

Predicting underwater sound propagation *** READ PROFILE. I0 READ(NIU, *)ZV,CV NSVP = NSVP + 1 ZSVP(NSVP) = ZV CSVP(NSVP) = CV IF (ZV .LT. ZLYR(I)) IF (ZV .GT. ZLYR(I)) IXSVP(I) = NSVP 20 CONTINUE RETURN C C

GO T0 I0 G0 T0 30

*** ERROR EXIT 30 NSVP = 0 RETURN 9999 STOP ' Input end-of-file on Unit 7 in Subroutine SVP.' END SUBROUTINE SFIELD IMPLICIT REAL*8 (A - H, 0 - Z), INTEGER*2 (I - M) *************************************************************** *** fiAUSSIAN STARTING FIELD - SEE NORDA TECH NOTE 12 BY H.K.BROCK *** SFIELD IS CALLED IF INPUT PARAMETER ISF = 0 *************************************************************** *** CALLING ROUTINE SUPPLIES: * FRq - FREQUENCY IN HZ * CO - REFEI~F~CE SOUND SPEED - METERS/SEC * ZS - DEPTH OF SOURCE IN METERS. * N - NUMBER OF POINTS IN ARRAY U * DZ - DEPTH INCREMENT - METERS *** SFIELD SUBROUTINE SUPPLIES: * U - COMPLEX STARTING FIELD

$INCLUDE : 'COMMON. IFD ' C C THE FIELD IS DEFINED AS A GAUSSIAN BEAM AT RANGE = 0. C LOCAL VARIABLES - GA GAUSSIAN AMPLITUDE GW = 2.0D0 / XK0 GA = DSQRT(GW) / GW C WRITE(*, *)'JSR FIRST BIG LOOP' , GA, GW D0 10 I = I, N C WRITE(*, *) 'LOOP COUNT = ', I ZM -- I * DZ PR = GAUSS(GA,ZM,ZS,GW) - GAUSS(GA,-ZM,ZS,GW) U(I) = DCMPLX(PR, 0.0D0) I0 CONTINUE RETURN END DOUBLE PRECISION FUNCTION GAUSS(GA, Z, GD, GW) IMPLICIT REAL*8 (A - H, 0 - Z), INTEGER*2 (I - M) C INPUT - GA GAUSSIAN AMPLITUDE C OUTPUT - GAUSS = GA * DEXP(-((Z - GD) / GW)**2) C TEMPORARY VARIABLE - V V = (Z - GD) / GW v -- - ( v , v ) IF(V .LT. -200.ODO) THEN GAUSS = O.ODO ELSE GAUSS = GA * DEXP(V) ENDIF RETURN END SUBROUTINE USVP IMPLICIT REAL*8 C C C C C C C C C C

(A - H, 0 - Z), INTEGER*2

(I - M)

*** U S E R SOUND VELOCITY PROFILE SUBROUTINE * SUBROUTINE USVP IS CALLED EACH DR IN RANGE AS LONG AS * RSVP IS N0T ZERO. RSVP MAY BE USED BY USER T0 TRANSFER CONTROL * IN THIS SUBROUTINE. USER INSERTS LOGIC TO CLEAR KSVP * WHEN USVP IS NO LONGER NEEDED. IF RSVP NOT CLEARED BY USER, * USVP IS CALLED EACH STEP IN RANGE UNTIL BA = NEXT RSVP. ****************************************************************** *** USVP SUBROUTINE RETURNS: " N L Y R - N U M B E R OF LAYERS. L A Y E R I IS WATER. OTHERS ARE SEDIMENT

755

756

J.S. ROBERTSON,W. L. SIEGMANN, and M. J. JACOBSON C C C C C C C

* *

ZLYR - ARRAY - DEPTH OF EACH LAYER. FIRST IS DEPTH OF WATER. RHO - ARRAY - DENSITY OF EACH LAYER. GRAMS/CUBIC CM * B E T A - ARRAY - ATTENUATION IN EACH LAYER. DR/WAVELENGTH * IXSVP - ARRAY - CONTAINS POINTERS. POINTS TO LAST VALUE OF SVP * IN CORRESPONDING LAYER. SVP IS STORED IN ARRAYS ZSVP * AND CSVP. IXSVP(1) POINTS TO LAST SVP POINT IN WATER. * NSVP - NUMBER OF POINTS IN ZSVP AND CSVP. ZSVP AND CSVP * CONTAIN %'HE PROFILES FOR ALL LAYERS. * ZSVP - A R ~ Y - SVP DEPTHS - METERS * CSVP - ARRAY - SOUND SPEED - METERS/SEC * KSVP - AS DESCRIBED ABOVE. ******************************************************************

$INCLUDE : 'COMMON. IFD ' C IF (RSVP .EQ. I) THEN C C *** IF KSVP=I, CONTROL IS TRANSFERRED HERE. USER LOADS C NLYR,ZLYR(I),RBO(I) ,BETA(I), AND IXSVP(I) WHERE I=I,NLYR. C USER ALSO LOADS NSVP,ZSVP(I), AND CSVP(I) WHERE I--1,NSVP. C KSVP MAY BE ALTERED DEPENDING ON USER LOGIC. C C *** USER SUPPLIES SVP C C *** USER INSERTS CODE HERE IF DESIRED ELSEIF (RSVP .EQ. 2) THEN C *** USER INSERTS CODE HERE IF DESIRED C ELSEIF (ESVP .EQ. S) THEN C *** USER INSERTS CODE HERE IF DESIRED C ELSEIF (ESVP .EQ. 4) THEN C *** USER INSERTS CODE HERE IF DESIRED ELSE

C C C C C

C C

C

NSVP -- 0 ENDIF RETURN END SUBROUTINE BCON IMPLICIT REAL*8 (A - H, 0 - Z), INTEGER*2 (I - M) ***************************************************************** *** USER PREPARED BOTTOM CONDITION SUBROUTINE * B C O N IS C A L L E D IF I N P U T P A R A M E T E R I T Y P E B = 1 * SEE MAIN PROGRAM FOR DEFINITIONS *************************************************%%***************

*** SUBROUTINE RETURNS: * BOTY,BOTX *****************************************************************

C $INCLUDE : 'COMMON. IFD ' C

IF (THETA) 10, 2 0 , 30 C C

C

C C

C

*** THETA LESS THAN 0 . 0 .

BOTTOM SLOPES UP.

I0 CONTINUE BOTY = U(N) B O T X ............... RETURN *** THETA = 0.0. BOTTOM IS FLAT. 20 CONTINUE BOTY = U(N) BOTX................ RETURN

C C

C C

*** THETA GREATER THAN 0 . 0 , 30 CONTINUE BOTY= .............. B O T X = . ............. RETURN END

B O T T O M SLOPES DOWN.

SUBROUTINE UFIELD IMPLICIT REAL*8 (A - H, 0 - Z), INTEGER*2

(I - M)

Predicting underwater sound propagation

757

***************************************************************

*** USER STARTING FIELD *** USER WRITES THIS SUBROUTINE IF GAUSSIAN FIELD NOT DESIRED *** UFIELD IS CALLED IF INPUT PARAMETER ISF IS NOT ZERO ***************************************************************

*** UFIELD SUBROUTINE SUPPLIES: U - COMPLEX STARTING FIELD ****.22.**********************************.2.*.2.**.2222.*.2.**

$INCLUDE:'COMMON.IFD' C C *** STARTING FIELD GENERATED BY USER DO I0 I = i, N C U(I) . . . . . 10 CONTINUE RETURN END SUBROUTINE READER IMPLICIT REAL*8 (A - H, 0 - Z), INTEGEN*2 (I - M) *$*2*****************************************2****2*$****2***** C C * This subroutine sets up the input a n d output file names. * C * User must set up a data file named NAMES.IFD with the * C * followining format: * C 2 L I N E I: filenamel - the data file * C * LINE 2:filename2 - output (formatted) containing a * C * of model data. , C * LINE 3:filename3 - output (formatted) containing the * C * the solution vector. , C * LINE 4 : f i l e n a m e 4 - output (unformatted) containing the * C * propngaton loss. , * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C CHABACTER*32 DATA, SUMMRY, VECTOR, PROPLS COMMON /FILES/ DATA, SUMMRY, VECTOR, PROPLS OPEN (8, FILE = 'NAMES.IFD', I STATUS = 'OLD', IOSTAT = I F ~ ) REWIND(8) IF (IERR .NE. O) THEN WRITE(*, *)' IEBR = ', I E ~ CLOSE ( 8 ) WRITE(*, 10)

FORMAT(' The file containing filenames cannot b e opened. ,/' ProEram execution aborted') STOP ' 2.. Program aborted in subroutine Reader. ***' ENDIF READ (8, 20) DATA, SUMMRY, VECTOR, PROPLS 20 FOB24AT(A32, /, AS2, /, A32, /, A32) CLOSE (8) RETURN END I0

1

'