COUPLE, November 20, 1997 Version November 20, 1997 Documentation by Richard B. Evans Woods Hole Oceanographic Institution 404 Bigelow Bldg., MS #9, Woods Hole, MA 02543, USA (508) 289-3513 and Science Applications International Corporation 21 Montauk Avenue, Suite 201, New London, CT 06320, USA (860) 442-0564 INTRODUCTION: COUPLE is a program that performs a stepwise coupled mode [1] calculation of the two dimensional underwater acoustic field due a vertical array of harmonic point or line sources in cylindrical or plane geometry respectively. The environment may be range dependent, consisting of a series of horizontally stratified fluid filled regions. Each of the regions consists of a water layer and a sediment or bottom layer that are assumed to have a piecewise linear depth variation in the index of refraction squared and a piecewise exponential depth variation in the density. The index of refracton squared and density are continuous in the water layer and in sediment layer, but they may be discontinuous at the water sediment interface. The index of refraction squared may be complex in both the water and bottom sediment layers to allow for attenuation. The modal calculation in each of the regions is accomplished with a search or Galerkin procedure. If the problem consist of a homogeneous water layer over a homogeneous bottom sediment layer and only the water depth varies with range, then the complex eigenvalues can be found using a search procedure based on Newton's method. In the more general case with inhomogeneous water and bottom sediment layers, range dependent material properties and range dependent bathymetry, the Galerkin matrix method is used. The Galerkin solution is constructed on a basis set from a simple background model consisting of a homogeneous water layer over a homogeneous bottom sediment layer, which retains the density discontinuity at the water sediment interface. The Galerkin procedure, with a finite basis set, yields a complex matrix eigenvalue problem. The complex eigenvalues are found using a Jacobi type algorithm designed for complex symmetric matrices. This routine does not order the eigenvalues by increaing real part. A description of the Galerkin method for the constant density water layer and homogeneous sediment layer case can be found in [2]. When density variations are present, the eigenvalue problem becomes a generalized eigenvalue problem and the Cholesky factorizaton is applied before the Jacobi type algorithm. The range dependence is described by a series of ranges and water depths along with the associated sound speed, density and attenuation profiles. A pressure release boundary is imposed at the surface of the water and at the greatest depth in the bottom sediment. Although the water depth varies with range, the total depth of the water and sediment is fixed. A decoupling algorithm [3] is used to solve the matrix two-point boundary value problem in range that occurs in the stepwise coupled mode formulation. This eliminates the numerical instability that occurs with other methods such as shooting. The purpose of this document is to describe the inputs and outputs of COUPLE. There is also a brief description of how to compile, link and run the FORTRAN source code that comprises the program. I. INPUTS: Block 0 (1 record): TITLE ------- Block 1 (1 record): HB,FREQ,ROHW,ROHB,DBPWL ------- HB = Total thickness of water and sediment layer (m). FREQ = Frequency (Hz). ROHW = Density of the water layer (g/cm**3). ROHB = Density of the sediment layer (g/cm**3). DBPWL = Attenuation in the sediment layer (dB/wavelength). Block 2 (1 record): CB,CW ------- CB = Sound speed in the sediment (m/s). CW = Sound speed in the water (m/s). Block 3 (1 record): M,NATEN,IEIG,NDPTH,IGEOM,IREFL ------- M = Number of contributing modes (Maximum=MX=190). NATEN = Number of incremental parts of DBPWL to use in the search for eigenvalues. IEIG = Flag used to indicate the source of the estimates of complex eigenvalues in region 2 through N in the search mode. IEIG=0 - The real eigenvalue problem is used and the attenuation is introduced, based on NATEN. IEIG=1 - The complex eigenvalues of the previous region are used to start. NDPTH = Number of incremental parts of the depth step to use with IEIG=1 in the search for complex eigenvalues. IGEOM = Flag used to indicate cylindrical or plane geometry. IGEOM=0 - Cylindrical geometry is used with point sources. IGEOM=1 - Plane geometry is used with line sources. IREFL = Flag to control the source condition at range zero. IREFL=0 - Perfectly cylindrically symmetric geometry is assumed. IREFL=1 - The source is backed by an absorber at negative ranges in cylindrical or plane geometry. Block 4 (1 record): NSDEP,ZSMIN,ZSINC ------- NSDEP = Number of source depths (Maximum=MD=400). ZSMIN = Minimum source depth (m). ZSINC = Increment in source depth (m) Block 5 (1 record): NDEP,ZMIN,ZINC ------- NDEP= Number of receiver depths in the pressure calculation (Maximum=MD=400). ZMIN = Minimum receiver depth (m). ZINC = Increment between the receiver depths (m). Block 6 (1 record): RMIN,RMAX,RINC ------- RMIN = Minimum range in the pressure calculation (km). RMAX = Maximum range in the pressure calculation (km). RINC = Increment between ranges in the pressure calculation (km). Block 7 (1 record): N,IPRT,IOUT,MAMP ------- N = Number of environmental regions, headed by range depth pairs, to follow (Maximum=NRX=2002). IPRT = Print flag for details of the calculation such as eigenvalues, coupling matrices and coefficients of outgoing and ingoing waves. IOUT = Flag that controls the matching condition used at the vertical interfaces in the range marching scheme. IOUT=0 - The program computes the fully elliptic two way solution that matches both the pressure and radial particle velocity. IOUT=1 - The one way pressure matching solution is found. IOUT=-1 - A single scatter approximation is used. MAMP = Flag to control generation of a file of mode amplitudes. ******** The following input block is repeated for J=1,N ******** Block 8(J) (1+NPW(J)+NPB(J) records): ---------- First record: RANG(J),IRLIN(J),DPTH(J),NPW(J),NPB(J),IREIG(J) RANG(J) = Range to the end of the J th region (km). IRLIN (J) = Number of additional regions generated in subroutine BATHY to represent a linearly sloping bottom between RANG(J) and RANG(J+1). DPTH (J) = Depth of the water in the J th region (m). NPW(J) = Number of points in the water sound speed, attenuations, density profile in the J th region. NPW(J)=0 - The water is assumed to be homogeneous with sound speed CW and the search procedure is used. NPW(J)>1 - The water column may be nonhomogeneous, a sound speed profile must be provided and the Galerkin procedure is used. NPB(J) = Number of points in the bottom sound speed, attenuation, density profile in the J th region. NPB(J) must follow the same two cases as NPW(J). NPB(J)=0 - The bottom is homogeneous with sound speed CB, density RHOB and attenuation DBPWL. NPB(J)>1 - The bottom may be nonhomogeneous and a bottom sound speed, attenuation, density profile must be provided. IREIG (J) = Sequence number for storage or retrieval of eigenvalues on a direct access file for use with perodic media only. IREIG(J)=0 - Omit storage on direct access file IREIG(J)>0 - Write eigenvalues on direct access file at record IREIG(J) IREIG(J)<0 - Read eigenvalues from direct access file at record ABS(IREIG(J)) Next NPW(J) records: DEPW(J,K),SVPW(J,K),DBPWLW(J,K),RHOW(J,K), K=1,NPW(J) DEPW(J,K) = Depth (m) of sound speed in the water (DEP(J,1)=0.0). SVPW(J,K) = Sound speed in the water (m/s). DBPWLW(J,K) = Attenuation in the water (dB/wavelength). RHOW(J,K) = Density in the water (g/cm**3). Last NPB(J) records: DEPB(J,L),SVPB(J,L),DBPWLB(J,L),RHOB(J,L), L=1,NPB(J) DEPB(J,L) = Depth (m) of sound speed in the bottom measured from the water surface (DEPB(J,1)=DEP(J,NPW(J))). SVPB(J,L) = Sound speed in the bottom (m/s). DBPWLB(J,L) = Attenuation in the bottom (dB/wavelength). RHOB(J,L) = Density in the bottom (g/cm**3). II. SUGGESTIONS AND COMMENTS: Block 1: ------- HB should be large enough to provide space for an absorber in the deepest part of the problem. The artificial absorber, where the attenuation increases significantly with depth, is often included at the bottom of the physical sediment. ROHW, ROHB and DBPWL are used in search mode and ignored in the Galerkin mode. Block 2: ------- CB is always nonzero. It is used to compute the index of refraction everywhere. CB is the bottom sound speed in the search mode and it is the bottom sound speed for the background two layer problem in every region in the Galerkin procedure. CW is used in the search mode and is ignored in the Galerkin mode. Block 3: ------- M should be chosen to be significantly larger than the number of propagating modes to represent high angles and to help in the expansion of the modes in the Galerkin procedure. Convergence should be checked by increasing M until the change in the calculated field is insignificant. NATEN, IEIG and NDPTH are only used in the search mode. IEIG=0 is always assumed in the first region. If IEIG=1, then NDPTH is used when redundant eigenvalues are found, which means the one was missed (See subroutine SEARCH). NATEN, IEIG and NDPTH can be set to zero, in the Galerkin mode, to avoid confusion. IREFL=0 yields the source condition described in [1 and 2] for perfectly cylindrically symmetric geometry. In this case, there exists standing waves, in range, in the region containing the source. This is a multiple scattering effect. In perfectly cylindrically symmetric geometry, the alternative IREFL=1 is nonphysical, but it is consistent with the single scatter approximation. IREFL=1 is always used in plane geometry. Block 4 ------- The sources may be in both the water layer and/or the bottom sediment layer. If NSDEP>1, then a vertical array of in-phase point or line sources projecting a horizontal beam is used. The field is found by a superpostion of the fields due to the NSDEP sources. The sum of sources is divided by NSDEP. See the main program COUPLE. Block 6 ------- RMIN can be zero, with IOUT=0, in the calculation of the scattered field that is regular at the source. The other fields are irregular at range zero. RMAX can be larger than the end of the last environmental region RNG(N). Block 7 ------- If IPRT>0, then details of the calculation are printed for the first, last and every IPRT region in between. If IOUT=-1, an approximation to the single scatter approximation due to Porter, Jensen and Ferla [4] is used. See subroutine PROPM for the implementation. If MAMP>0, then the mode amplitudes are written out in CPRES every MAMP range steps. This file can get very big so pick MAMP carefully. Block 8(J) : ---------- In long flat regions with IOUT=0, it may be necessary to include zero height bathymetry changes to activate the decoupling algorithm to normalize the solution and accurately compute ingoing and outgoing waves. RNG(J) is the range to the end of the J th region so RNG(1)>0. RNG(N) has no significance, since the last region extends indefinitely. If IRLIN(J)>0 then interpolation and extrapolation procedures are used. As the water depth increases (decreases) the sound speed, attenuation, density profile in the water will be extrapolated (interpolated). As the water depth increases (decreases),the bottom sound speed, attenuation, density profile is shifted down (up) and the interpolation (extrapolated) takes place in deepest part of sediment profile at HB. NPW(J) controls the procedure used to find eigenvalues. If NPW(J)=0, then the water is assumed to be homogeneous with sound speed CW. No depth sound speed pairs are input. The eigenvalues are found in subroutine SEARCH based on NATEN, NDPTH and IEIG. The analytic forms of the depth functions are used. If NPW(J)=2 or more, then the water column may be nonhomogeneous. A sound speed, attenuation, density profile must be input. The eigenvalues are found by the Galerkin method in subroutine GALRKN based on the Jacobi type matrix eigenvalue algorithm. The depth functions are constructed from the matrix eigenvectors. NPB(J) follows the same two cases as NPW(J). If NPW(J)=2 or more, then NPB(J)=2 or more and a bottom sound speed, attenuation, density profile must be provided. If NPW(J)=0 is then it is necesary that NPB(J)=0 and no bottom profile may be input. IREIG(J) is used with perodic media to avoid recomputation of identical sets of eigenvalues using subroutine SEARCH or GALRKN. It is not valid when IRLIN(J) is nonzero because it is too hard to anticipate the sequence numbers needed. It is not valid with a mixture of calls to subroutines SEARCH and GALRKN to make it easy to keep track of the records written on direct access file. RHOW(J,K) overrides the value of ROHW (Note the different spelling). DBPWLB(J,L) and RHOB(J,L) override the values of ROHB and DBPWL. III. EXAMPLES: The following four examples are chosen because the calculated transmission loss appears in the references [1-5] and they may be used to test the code. The Galerkin procedure is less prone to problems, but input for the search procedure is given, in some examples, for contrast. The results from both procedures should be essentially the same. Example 1. Kutschale Problem [1]. Kutschale Problem, Search procedure, 70 Modes 5000. 25.0 1.0 1.15 0.5 !HB,FREQ,ROHW,ROHB,DBPWL 1704.5 1500. !CB,CW 70 20 1 0 0 0 !M,NATEN,IEIG,NDPTH,IGEOM,IREFL 1 25.0 0.0 !NSDEP,ZSMIN,ZSINC 1 25.0 0.0 !NDEP,ZMIN,ZINC 0.05 5.1 .05 !RMIN,RMAX,RINC 2 1 0 0 !N,IPRT,IOUT,MAMP 2.50 0 32.5 0 0 0 !RANGE,IRLIN,DPTH,NPW,NPB,IREIG 5.10 0 32.5 0 0 0 Kutschale Problem, Galerkin procedure, 25 Modes 1000. 25.0 0.0 0.00 0.0 !HB,FREQ,ROHW,ROHB,DBPWL 1704.5 00.00 !CB,CW 25 0 0 0 0 0 !M,NATEN,IEIG,NDPTH,IGEOM,IREFL 1 25.0 0.0 !NSDEP,ZSMIN,ZSINC 1 25.0 0.0 !NDEP,ZMIN,ZINC 0.05 5.1 .05 !RMIN,RMAX,RINC 2 1 0 0 !N,IPRT,IOUT,MAMP 2.50 0 32.5 2 3 0 !RANGE,IRLIN,DPTH,NP,NPB,IREIG 0.0 1500.0 0.0 1.00 !DEPW,SVPW,DBPWLW,RHOW 32.5 1500.0 0.0 1.00 32.5 1704.5 0.5 1.15 !DEPB,SVPB,DBPWLB,RHOB 250.0 1704.5 0.5 1.15 1000.0 1704.5 2.5 1.15 5.10 0 32.5 2 3 0 0.0 1500.0 0.0 1.00 !DEPW,SVPW,DBPWLW,RHOW 32.5 1500.0 0.0 1.00 32.5 1704.5 0.5 1.15 !DEPB,SVPB,DBPWLB,RHOB 250.0 1704.5 0.5 1.15 1000.0 1704.5 2.5 1.15 Example 2. Square wave 5-10 km patch of roughness [2 and 5]. The published transmission loss plots are under sampled, at a range increment of .05 km, in both references. The .01 km range increment used here shows a good deal more structure. Patch of roughness, Square wave corrugation, Search procedure, 40 modes 1500. 25. 1.0 2.5 0.5 !HB,FREQ,ROHW,ROHB,DBPWL 1704.5 1500.0 !CB,CW 40 20 0 0 0 0 !M,NATEN,IEIG,NDPTH,IGEOM,IREFL 1 18.0 0.0 !NSDEP,ZSMIN,ZSINC 1 50.0 0.0 !NDEP,ZMIN,ZINC 0.00 18.1 .010 !RMIN,RMAX,RINC 101 10 0 0 !N,IPRT,IOUT,MAMP 5.00 0 100. 0 0 1 !RANGE,IRLIN,DPTH,NPW,NPB,IREIG 5.05 0 90. 0 0 2 5.10 0 100. 0 0 -1 5.15 0 90. 0 0 -2 (Repeat the last two records 47 times increasing both ranges by .1 km in each cycle. The -1,-2 sequence eliminates recomputation of eigevalues.) 9.90 0 100. 0 0 -1 9.95 0 90. 0 0 -2 20.00 0 100. 0 0 -1 Example 3. PE Workshop II, Test Case 3b [3]. PE Workshop II, Test Case 3b , Galerkin procedure, 28 modes 400. 250.0 0.0 0.0 0.0 !HB,FREQ,ROHW,ROHB,DBPWL 1590.0 00.00 !CB,CW 28 0 0 0 0 0 !M,NATEN,IEIG,NDPTH,IGEOM,IREFL 1 99.5 0.0 !NSDEP,ZSMIN,ZSINC 1 99.5 0.0 !NDEP,ZMIN,ZINC 4.9 10.1 .01 !RMIN,RMAX,RINC 2 1 0 0 !N,IPRT,IOUT,MAMP 5.00 0 100.0 2 2 0 !RANGE,IRLIN,DPTH,NPW,NPB,IREIG 0.0 1500.0 0.0 1.0 !DEPW,SVPW,DBPWLW,RHOW 100.0 1500.0 0.0 1.0 100.0 1590.0 0.5 1.2 !DEPB,SVPB,DBPWLB,RHOB 400.0 1590.0 0.5 1.2 5.10 0 100.0 2 2 0 0.0 1500.0 0.0 1.0 !DEPW,SVPW,DBPWLW,RHOW 100.0 1500.0 0.0 1.0 100.0 1590.0 0.5 1.2 !DEPB,SVPB,DBPWLB,RHOB 400.0 1590.0 0.5 1.2 Example 4. Acoustical Society of America (ASA) Benchmark Wedge [4]. ASA Benchmark Wedge, Galerkin procedure, 45 modes, 200 steps 1500. 25.0 1.0 1.50 0.0 !HB,FREQ,ROHW,ROHB,DBPWL 1700.0 00. !CB,CW 45 0 0 0 0 0 !M,NATEN,IEIG,NDPTH,IGEOM,IREFL 1 100.0 0.0 !NSDEP,ZSMIN,ZSINC 2 30.0 120.0 !NDEP,ZMIN,ZINC 0.01 4.01 .01 !RMIN,RMAX,RINC 2 50 0 0 !N,IPRT,IOUT,MAMP 0.01 200 200.0 2 3 0 !RANGE,IRLIN,DPTH,NPW,NPB,IREIG 0.0 1500.0 .0 1.00 !DEPW,SVPW,DBPWLW,RHOW 200.0 1500.0 .0 1.00 200.0 1700.0 .5 1.50 !DEPB,SVPB,DBPWLB,RHOB 1000.0 1700.0 .5 1.50 1500.0 1700.0 2.5 1.50 4.000 0 0.01 2 3 0 0.0 1500.0 .0 1.00 !DEPW,SVPW,DBPWLW,RHOW 0.01 1500.0 .0 1.00 0.01 1700.0 .5 1.50 !DEPB,SVPB,DBPWLB,RHOB 1000.0 1700.0 .5 1.50 1500.0 1700.0 2.5 1.50 ASA Benchmark Wedge, Search procedure, 90 modes, 200 steps 3000. 25.0 1.0 1.50 0.5 !HB,FREQ,ROHW,ROHB,DBPWL 1700.0 1500.0 !CB,CW 90 30 1 20 0 0 !M,NATEN,IEIG,NDPTH,IGEOM,IREFL 1 100.0 0.0 !NSDEP,ZSMIN,ZSINC 2 30.0 120.0 !NDEP,ZMIN,ZINC 0.01 4.01 .01 !RMIN,RMAX,RINC 2 50 0 0 !N,IPRT,IOUT,MAMP 0.01 200 200.0 0 0 0 !RANGE,IRLIN,DPTH,NPW,NPB,IREIG 4.000 0 0.01 0 0 0 IV. INPUT AND OUTPUT FILES: COUPLE.DAT = Free format input data as described above. COUPLE.LOG = Formatted log file (may appear on screen). COUPLE.PRT = Formatted print file. COUPLE.TL = Formatted transmission loss. COUPLEO.TL = Formatted outgoing transmission loss if IOUT=0 COUPLEI.TL = Formatted ingoing transmission loss if IOUT=0 COUPLES.TL = Formatted scattered transmission loss in region 1 if IOUT=0 COUPLE.CPR = Unformatted real and imaginary parts of the complex pressure. COUPLEO.CPR = Unformatted outgoing complex pressure if IOUT=0 COUPLEI.CPR = Unformatted ingoing complex pressure if IOUT=0 COUPLES.CPR = Unformatted scattered complex pressure in region 1 if IOUT=0 COUPLE.AMP = Unformatted mode amplitudes when MAMP>0 V. COMPILING, LINKING AND RUNNING: Windows NT or 95, Command window (DOS): There is a batch file called "FALLSUBS.BAT" that compiles all the subroutines (70) using the Lahey Fortran 90 compiler. The Fortran source code has names of the form *.FOR. The compiler option 'vax' is used along with targeting for your your processor. To compile the subroutines type "FALLSUBS.BAT" at the DOS prompt in command window. There is a batch file "FLCOUPLE.BAT" that compiles the main program COUPLE and links it with the subroutines that are listed in the file "COUPLE.LNK" read by the linker. To compile and link the main program COUPLE, type "FLCOUPLE.BAT" at the DOS prompt in the command window. There is a batch file called "COUPLE.BAT" that imports the name of the *.dat input file, runs the program and handles the output file naming. To run COUPLE create an input file called "name.dat" and type "COUPLE.BAT name" at the DOS prompt in the command widow. When the program finishes, some or all of the output files listed above are created with "name" in place of COUPLE and with the various extensions. The string "name" should be no more than 7 characters for backward compatibiitly with DOS, where the maximun is 8. The format of the output files is established in the subroutines HEADER and CPRES. Unix: There is a unix make file called "makefile" that lists all the subroutines. To compile and link type "make couple" at the unix prompt. There is a executable script file called "couple.bat" that runs the program and handles the input and output file naming. To run couple perform the following steps. First create an input file, as described above, with a unique name like "newrun.dat", Next use the vi editor to replace the "oldrun" from the last run with the "newrun" everywhere in couple.bat. This can be done with the vi command "%s/oldrun/newrun/g" at the ":" prompt in vi. Exit vi. Lastly run couple in the batch mode by typing the three lines "batch " and "couple.bat " and "control d " at the unix prompt. It is necessary to open a file in the main program, say COUPLE.LOG, to except the output written with WRITE(* that normally comes to the screen in a DOS window. It is also be necessary to change timming calls in the main program COUPLE and the date call in HEADER, which are system dependent. REFERENCES: 1. R. B. Evans, "A coupled mode solution for acoustic propagation in a waveguide with stepwise depth variations of a penetrable bottom," J. Acoust. Soc. Am. 74, 188-195 (1983). 2. R. B. Evans and K. E. Gilbert, "Acoustic propagation in an refracting ocean waveguide with an irregular interface," Comp. Maths. Appl. 11, 795-805 (1985). 3. R. B. Evans, "The decoupling of stepwise coupled modes," J. Acoust. Soc. Am. 80, 1414-1418 (1986). 4. M. B. Porter, F. B. Jensen and C. M. Ferla, "The problem of energy conservation on one way models," J. Acoust. Soc. Am. 89, 1058-1067 (1991). 5. R. B. Evans and K. E. Gilbert, "The perodic extension of stepwise coupled modes," J. Acoust. Soc. Am. 77, 983-988 (1985).