
        COUPLE, March 24,1995 Version, September 13, 1996 Documentation
     Richard B. Evans, 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 is assumed to have a
piecewise linear depth variation in the index of refraction squared, in the
water and sediment layers, as well as an exponential depth variation in the
sediment density. The index of refraction squared in the sediment may be
complex 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 complex eigenvalues are found using the QR matrix algorithm. A
description of the Galerkin method for the homogeneous sediment layer
case can be found in [2].

   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,NEWSBC
-------
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.
NEWSBC = Flag to control the source condition at range zero.
           NEWSBC=0 - Perfectly cylindrically symmetric geometry is assumed.
           NEWSBC=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+NP(J)+NPB(J) records):
----------

First record: RANG(J),IRLIN(J),DPTH(J),NP(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).
NP(J) = Number of points in the sound speed profile in the J th region.
        NP(J)=0 - The water is assumed to be homogeneous with sound speed CW
                  and the search procedure is used.
        NP(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 NP(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 NP(J) records: DEP(J,K),SVP(J,K), K=1,NP(J) 

DEP(J,K) = Depth (m) of sound speed in the water (DEP(J,1)=0.0).
SVP(J,K) = Sound speed in the water (m/s).

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,NP(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.
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 mode.
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. NDPTH is used if 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.

NEWSBC=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 
NEWSBC=1 is nonphysical, but it is  consistent with the single scatter
approximation. NEWSBC=1 is always used in plane geometry.

Block 4
-------
The sources are assumed to be in the water layer where the density is ROHW.
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 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.

NP(J) controls the procedure used to find eigenvalues.

If NP(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 NP(J)=2 or more, then the water column may be nonhomogeneous. A sound 
speed profile must be input. The eigenvalues are found by the Galerkin 
method in subroutine GALRKN based on the QR matrix eigenvalue algorithm.
The depth functions are constructed from the matrix eigenvectors.

NPB(J) follows the same two cases as NP(J). If NP(J)=2 or more, then 
NPB(J)=2 or more and a bottom sound speed, attenuation, density  profile
must be provided. If NP(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.

DBPWLB(J,L) and RHOB(J,L) override the values of ROHB and DBPWL.

III. EXAMPLE:

SUBBOTTOM DENSITY GRADIENT STEP, 2W, HB=400 M, 50 MODES
400.  100.0   1.0   0.00   0.00              ! HB,FREQ,ROHW,ROHB,DBPWL
1847.   0.00                                 ! CB,CW
50  0  0  0  0  0                            ! M,NATEN,IEIG,NDPTH,GEOM,IREFL
1   50.0    0.0                              ! NSDEP,ZSMIN,ZSINC
2   45.0   80.0                              ! NDEP,ZMIN,ZINC
0.00   4.1   .005                            ! RMIN,RMAX,RINC
2   1   0   0                                ! N,IPRT,IOUT,MAMP
2.00     0    73.   8  4   0                 ! RANGE,IRLIN,DPTH,NP,NPB,IREIG
  0.0    1520.0                              ! DEP,SVP
 15.0    1520.0
 20.0    1505.0
 30.0    1492.0
 40.0    1484.0
 50.0    1485.0
 60.0    1487.0
 73.0    1487.0
 73.0    1560.0      0.0818    1.86           ! DEPB,SVPB,DBPWLB,RHOB
 83.0    1610.0      0.1337    1.96
 93.0    1847.0      0.1826    2.13
400.0    1847.0      5.0000    2.13 
4.50     0    73.   8  3   0                  ! RANGE,IRLIN,DPTH,NP,NPB,IREIG
  0.0    1520.0                               ! DEP,SVP
 15.0    1520.0
 20.0    1505.0
 30.0    1492.0
 40.0    1484.0
 50.0    1485.0
 60.0    1487.0
 73.0    1487.0
 73.0    1847.0      0.1826    2.13           ! DEPB,SVPB,DBPWLB,RHOB
 93.0    1847.0      0.1826    2.13
400.0    1847.0      5.0000    2.13

IV. INPUT AND OUTPUT FILES:

COUPLE.DAT = Input data as described above.

COUPLE.LOG = Free format log file.
COUPLE.PRT = Formatted print file.
COUPLE.CPR = Real and imaginary parts of the complex pressure.
COUPLE.TL = Transmission loss.
COUPLE.AMP = Mode amplitudes when MAMP>0
COUPLEO.TL = Outgoing transmission loss if IOUT=0
COUPLEI.TL = Ingoing transmission loss if IOUT=0
COUPLEO.CPR = Outgoing complex pressure if IOUT=0
COUPLEI.CPR = Ingoing complex pressure if IOUT=0
COUPLES.CPR = Scattered complex pressure in region 1 if IOUT=0
COUPLES.TL = Scattered transmission loss in region 1 if IOUT=0

V. COMPILING, LINKING AND RUNNING:

  There is a unix make file called "makefile" that lists all the subrout-
ines. 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 <cr>" and  "couple.bat <cr>" and "control d <cr>" at the unix prompt.

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).
