      PROGRAM COUPLE
C
C     THIS PROGRAM IS THE DRIVER FOR A COUPLED NORMAL MODE
C     CALCULATION BASED ON A MODEL CONSISTING OF A HOMOGENEOUS
C     WATER LAYER OVER A HOMOGENEOUS BOTTOM SEDIMENT LAYER. IT IS
C     EXTENDED, USING THE GALERKIN METHOD, TO ALLOW A PIECEWISE 
C     LINEAR VARIATIONS OF THE INDEX OF REFRACTION SQUARED IN THE 
C     WATER AND SEDIMENT LAYERS AS WELL AS EXPONENTIAL VARIATIONS
C     OF THE SEDIMENT DENSITY.
C
C     A PRESSURE RELEASE BOUNDARY IS IMPOSED AT THE SURFACE OF THE
C     WATER AND AT THE GREATEST DEPTH IN THE BOTTOM. THE WATER DEPTH
C     MAY VARY WITH RANGE BUT THE TOTAL DEPTH OF WATER AND SEDIMENT 
C     IS FIXED.
C
C     A DECOUPLING ALGORITHM IS USED TO SOLVE THE MATRIX TWO-POINT
C     BOUNDARY VALUE PROBLEM THAT OCCURS IN THE STEPWISE COUPLED
C     MODE FORMULATION. THIS ELIMINATES THE NUMERICAL PROBLEM
C     INHERENT IN THE SHOOTING METHOD THAT WAS USED IN AN EARLIER 
C     VERSION OF THE PROGRAM.
C
C     IMPLICIT REAL*8 (A-H,O-Z) ADDED TO MAIN AND ALL SUBROUTINES
C     TO ELIMINATE MIXED MODE PROBLEMS (R. B. EVANS, SAIC, NEW 
C     LONDON, CT, 10/23/92)
C
C     MODIFICATIONS TO INCORPORATE A COMPLEX PIECEWISE LINEAR 
C     INDEX OF REFRACTION SQUARED IN THE BOTTOM.
C     (R. B. EVANS, SAIC, NEW LONDON, CT, 3/29/94)
C
C     MODIFICATIONS TO INCORPORATE AN EXPONENTIAL DENSITY PROFILE
C     IN THE BOTTOM. (R. B. EVANS, SAIC, NEW LONDON, CT, 2/21/95)
C
C     MODIFICATIONS TO INCORPORATE MULTIPLE POINT SOURECES FOR 
C     DIRECTIONALITY (R. B. EVANS, SAIC, NEW LONDON, CT, 3/3/95)
C
C     THE MAIN VARIABLES ARE AS FOLLOWS.
C
C        INPUT-
C          HB=TOTAL DEPTH OF WATER AND BOTTOM LAYER (M)
C          FREQ=FREQUENCY (HZ)
C          ROHW=DENSITY IN THE WATER (G/CM**3)
C          ROHB=DENSITY IN THE BOTTOM IN THE SEARCH MODE (G/CM**3)
C          DBPWL=ATTENUATION IN THE BOTTOM (DB/WAVE LENGTH)
C                IN THE SEARCH MODE WHEN THE BOTTOM IS ISOSPEED.
C                IT IS IGNORED IN THE GALERKIN MODE.
C
C          CB=SOUND SPEED IN THE BOTTOM (M/SEC). CB IS USED
C             TO COMPUTE THE INDEX OF REFRACTION EVERY WHERE.
C             IT IS THE BOTTOM SOUND SPEED FOR THE BACKGROUND
C             TWO LAYER PROBLEM IN THE GALERKIN MODE.
C          CW=SOUND SPEED IN THE WATER (M/SEC) IN THE SEARCH
C             MODE WHEN THE WATER IS ISOSPEED. IT IS IGNORED
C             IN THE GALERKIN MODE. 
C
C          M=NUMBER OF CONTRIBUTING MODES (MAX.=MX=190)
C          NATEN=NUMBER OF PARTS OF ATEN TO USE IN INCREMENTAL (ON ATEN)
C                SEARCH FOR EIGEN VALUES (SEE SEARCH)
C          IEIG=FLAG TO INDICATE SOURCE OF ESTIMATES OF COMPLEX EIGEN
C               VALUES IN REGION 2 THROUGH N. IF IEIG=0 THE REAL PROB.
C               IS USED (SEE SEARCH)
C          NDPTH=NUMBER OF PARTS OF DEPTH STEP TO USE WITH IEIG=1 IN
C                INCREMENTAL (ON DEPTH) SEARCH FOR EIGEN VALUES. NDPTH
C                IS ONLY USED IF MULTIPLE EIGEN VALUES
C                ARE INCOUNTERED (SEE SEARCH)
C          IGEOM=FLAG WHICH WHEN SET TO ONE (1) CAUSES THE CODE TO DO 
C                THE CALCULATION IN PLANE GEOMETRY WITH A LINE SOURCE.
C          NEWSBC=FLAG WHICH WHEN SET TO ONE (1) INVOKES A RADIATION
C                 BOUNDARY CONDITON AT THE SOURCE THAT IS COMPATIBLE 
C                 WITH THE SINGLE SCATTER APPROXOMATION IN CYLINDRICAL
C                 GEOMETRY AND IS AWAYS USED WITH PLANE GEOMETRY.
C
C          NSDEP=NUMBER OF SOURCE DEPTHS (MAX.=MD=400)
C          ZSMIN=MIN. SOURCE DEPTH (M)
C          ZSINC=INC. IN SOURCE DEPTH (M)
C
C          NDEP=NUMBER OF RECEIVER DEPTHS FOR PRESS. CAL. 
C               (MAX.=MD=400)
C          ZMIN=MIN. RECEIVER DEPTH FOR PRESS. CAL.(M)
C          ZINC=INC. IN RECIEVER DEPTH FOP PRESS. CAL. (M)
C
C          RMIN=MIN. RANGE FOR PRESSURE CALCULATION (KM)
C          RMAX=MAX. RANGE FOR PRESSURE CALCULATION (KM)
C          RINC=INC. IN RANGE FOR PRESSURE CALCULATION (KM)
C
C          N=NUMBER OF INPUT RANGE DEPTH PAIRS (RANG,DPTH)
C          IPRT=PRINT FLAG FOR DETAILS OF CALCULATION SUCH AS EIGEN
C               VALUES AND COEFFICIENTS OF OUTGOING AND BACKSCATTERED
C               WAVES. IF IPRT .GE. 0 THE RESULTS ARE PRINTED OUT
C               FOR THE FIRST AND THE LAST REGION AND EVERY IPRT
C               REGION IN BETWEEN.
C          ISOUT=FLAG FLAG WHICH WHEN SET TO ZERO (0) COMPUTES THE
C                FULLY ELLIPTIC TWO WAY SOLUTION THAT MATCHES BOTH 
C                PRESSURE AND RADIAL PARTICAL VELOCITY. IF IOUT IS 
C                SET TO ONE (1), THE ONE-WAY PRESSURE MATCHING 
C                SOLUTION IS FOUND. IF IOUT IS SET TO MINUS ONE (-1),
C                AN APPROXIMATION TO THE SINGLE SCATTER APPROXIMATION
C                DUE TO PORTER AND JENSEN IS USED (SEE PROPM).
C          MAMP=FLAG TO GENERATE A FILE OF MODE AMPLITUDES ON UNIT11.
C               IF MAMP .GT. 0 THEN THE MODE AMPLITUDES ARE WRITTEN 
C               OUT IN CPRES EVERY MAMP RANGE STEPS. THIS FILE CAN GET 
C               VERY BIG SO PICK MAMP CAREFULLY. 
C
C               ** ALL THE REMAINING INPUTS ARE REPEATED FOR J=1,N **
C
C          RANG=RANGE TO END OF J TH REGION (KM)
C          DPTH=DEPTH OF WATER IN J TH REGION (M)
C          IRLIN=NUMBER OF ADDITIONAL REGIONS GENERATED TO REPRESENT
C                A LINEARLY SLOPING BOTTOM BETWEEN RANG(J) AND
C                RANG(J+1). SEE SUBROUTINE BATHY.
C          NP=NUMBER OF POINTS IN THE SOUND SPEED PROFILE IN THE J TH
C             REGION. IF NP .LE. 1 (IMAT=0) THEN THE WATER IS ASSUMED
C             TO BE HOMOGENEOUS (CW). NO DEPTH SOUND SPEED PAIRS ARE 
C             INPUT. THE EIGENVALUES ARE FOUND IN SUBROUTINE SEARCH
C             BASED ON NATEN, NDPTH AND IEIG. THE ANALYTIC FORM OF 
C             THE DEPTH FUNCTIONS IS USED. IF NP .GE. 2 (IMAT=1) THEN 
C             THE WATER COLUMN MAY BE NONHOMOGENEOUS. A SOUND SPEED
C             PROFILE MUST BE INPUT. THE EIGENVALUES ARE FOUND BY THE
C             GALERKIN MATRIX METHOD IN SUBROUTINE GALRKN (NATEN,NDPTH
C             AND IEIG ARE IGNORED). THE DEPTH FUNCTIONS ARE CONSTRUCTED
C             FROM THE MATRIX EIGENVECTORS.
C          NPB=NUMBER OF POINTS IN THE BOTTOM SOUND SPEED, ATTENUATION,
C             DENSITY PROFILE IN THE J TH REGION. IF NP .GE. 2 THEN NPB 
C             MUST BE .GE. 2 AND A BOTTOM SOUND SPEED, ATTENUATION, 
C             DENSITY PROFILE MUST BE PROVIDED. IF NP IS .LT. 2 THEN 
C             NO BOTTOM PROFILE MAY BE INPUT.
C          IREIG=SEQUENCE NUMBER FOR STORAGE OR RETRIEVAL OF
C                EIGENVALUES ON THE DIRECT ACCESS FILE UNIT=18.
C                IREIG(J) .EQ. 0 : OMIT STORAGE ON UNIT=18
C                IREIG(J) .GT. 0 : WRITE EIGENVALUES ON UNIT=18 AT
C                RECORD DETERMINED BY IREIG(J) 
C                IREIG(J) .LT. 0 : READ EIGENVALUES FROM UNIT=18 AT
C                RECORD DETERMINED BY ABS(IREIG(J))
C                THIS OPTION IS USED TO AVOID RECOMPUTATION OF
C                IDENTICAL EIGENVALUES USING SEARCH AND GALRKN (QR).
C                IT IS NOT VALID WHEN IRLIN .NE. ZERO BECAUSE IT IS
C                TOO HARD TO ANTICIPATE THE SEQUENCE NUMBERS NEEDED.
C                IT IS NOT VALID WITH A MIXTURE OF SEARCH AND GALRKN
C                TO MAKE IT EASY TO KEEP TRACK OF THE NUMBER OF 
C                RECORDS WRITEN ON UNIT=18.
C
C          DEP=DEPTH (M) OF SOUND SPEED IN THE WATER (DEP(1)=0.0)
C          SVP=SOUND SPEED IN THE WATER (M/SEC)
C
C          DEPB=DEPTH (M) OF SOUND SPEED IN THE BOTTOM (MEASURED FROM
C               THE WATER SURFACE AND DEPB(1)=DEP(NP))
C          SVPB=SOUND SPEED IN THE BOTTOM (M/SEC)
C          DBPWLB=ATTENUATION IN THE BOTTOM (DB/WAVE LENGTH). THIS
C                 INPUT OVERIDES THE DBPWL IN THE FIRST RECORD.
C          RHOB= DENSITY IN THE BOTTOM (G/CM**3)
C
C        INTERMEDIATE-
C
C          ATEN=IMAGINARY PART OF SQUARE OF INDEX OF REFRACTION IN
C                THE BOTTOM. SEE SUBROUTINE INPUT AND ESPI. USED
C                ONLY IN THE SEARCH MODE WHEN THE BOTOM IS ISOSPEED.
C          ZZ=REAL EIGENVALUES
C          CW=AVERAGE SOUND SPEED IN WATER
C          IMAT=FLAG WHICH CONTROLS THE METHOD FOR FINDING COMPLEX 
C               EIGENVALUES AND HOW THE DEPTH FUNCTIONS ARE CONSTRUCTED.
C               IMAT IS SET IN SUBROUTINE INPUT BASED ON NP.
C          BSOUT=INTEGER FLAG WHICH IS SET TO ABS(IOUT) IN INPUT. IF 
C                SET TO ONE (1) IT TURNS OFF CALCULATION OF BACKSCATTER
C                AND YIELDS AN OUTGOING SOLUTION WHICH IS EITHER THE
C                THE ONE-WAY OR SINGLE SCATTER APPROXIMATION DEPENDING
C                ON IOUT. IF BSOUT IS EQUAL TO ZERO (0), THE TWO-WAY 
C                CALCULATION IS DONE. BSOUT IS CALLED IOUT IN ALL 
C                SUBROUTINES EXCEPT INPUT AND PROPM.
C          LH=TNE WATER DEPTH HW IS BETWEEN DEP(LH) AND DEP(LH+1).
C              IF LH=NP THEN EXTRAPOLATION IS DONE.
C          LHB=THE BOTTOM DEPTH HB IS BETWEEN DEPB(LHB) AND DEPB(LHB+1).
C              IF LHB=NPB THEN EXTRAPOLATION IS DONE.
C          GRAD=GRADIENTS IN THE REPRESENTATION OF THE REAL INDEX OF 
C               REFRACTION SQUARED IN THE WATER. SEE SUBROUTINE RFACIN.
C          FINT=INTERCEPTS IN THE REPRESENTATION OF THE REAL INDEX OF 
C               REFRACTION SQUARED IN THE WATER. SEE SUBROUTINE RFACIN.
C          AVG=AVERAGE OF THE SQUARE OF THE INDEX OF REFRACTION
C              IN THE WATER IN THE CURRENT DEPTH REGION. SEE SUBROUTINE
C              AVGIN. 
C          AVGW=STORED VALUES OF AVG FOR USED IN FIELD CALUCULATION.
C          CGRADB=GRADIENTS IN THE REPRESENTATION OF THE COMPLEX
C                 INDEX OF REFRACTION SQUARED IN THE BOTTOM. 
C                 COMPUTED IN CFACIN AND STORED ON LUSVP.
C          CINTB=INTERCEPTS IN THE REPRESENTATION OF THE COMPLEX 
C                INDEX OF REFRACTION SQUARED IN THE BOTTOM.
C                COMPUTED IN CFACIN AND STORED ON LUSVP.
C          AVGB=THE REAL PART OF THE SQUARE OF THE BACKGROUND INDEX
C               OF REFRACTION IN THE BOTTOM. THIS IS CURRENTLY SET TO
C               ONE FOR ALL REGIONS IN BOTDEPS AND ADDED TO THE 
C               NEGATIVE THE TRUE INDEX OF REFRACTION IN AMAT.
C          GRHOB=LOGARITHMIC GRADIENTS OF THE DENSITY IN THE BOTTOM
C                (1/M). COMPUTED IN SUBROUTINE INPUT AND STORED ON
C                LUSVP.
C          ROHBL=DENSITY OF BOTTOM ON LEFT SIDE OF VERTICAL INTERFACE
C                IN CRSCUP THAT WAS USED TO FIND THE BASIC EIGENVALUES
C                AND EIGENFUNCTIONS (IMAT=1).
C          ROHBR=DENSITY OF BOTTOM ON RIGHT SIDE OF VERTICAL INTERFACE
C                IN CRSCUP THAT WAS USED TO FIND THE BASIC EIGENVALUES
C                AND EIGENFUNCTIONS (IMAT=1).
C          BEGVL=BASIC SET OF EIGENVALUES ON LEFT (IMAT=1)
C          BEGVR=BASIC SET OF EIGENVALUES ON RIGHT (IMAT=1)
C          BMATL=MATRIX OF NORMALIZED EIGENVECTORS CORRESPONDING
C                TO BEGVL (COEFFICIENTS FOR THE DEPTH FUNCTIONS)
C          BMATR=MATRIX OF NORMALIZED EIGENVECTORS CORRESPONDING
C                TO BEGVR (COEFFICIENTS FOR THE DEPTH FUNCTIONS)
C          EGVL=EIGEN VALUE ON THE LEFT
C          EGVR=EIGEN VALUE ON THE RIGHT
C          FNORM=NORMALIZATION FACTOR
C          CLR=CROSS COUPLING (LEFT TO RIGHT)
C          CRL=CROSS COUPLING (RIGHT TO LEFT)
C          R=PROPAGATOR MATRIX   
C          T=UNITARY TRANSFORMATION MATRIX, R(J)*T(J)=T(J+1)*U(J)
C          U=UPPER TRIANGULAR MATRIX
C          Y1=DECOUPLED OUTGOING SOLUTION
C          Y3=DECOUPLED INGOING SOLUTION
C          C1=COEFFICIENT OF OUTGOING WAVE IN REGION 1
C          A=COEFICIENT OF OUTGOING WAVE
C          B=COEFICIENT OF INGOING WAVE
C          ZSD=SOURCE DEPTHS IN SOURCE ARRAY
C          ZS=MIDPOINT OF SOURCE ARRAY WRITTEN IN HEADER OF
C             OUTPUT FILES
C          ZR=RECIEVER DEPTHS
C          CDF=DEPTH FUNCTIONS AT OUTPUT DEPTHS
C          RANGE,DEPTH=SEQUENCE GENERATED IN SUBROUTINE BATHY  
C          NREG=NUMBER OF REGIONS
C
C        COMPUTED(IN CPRES)-
C
C          P=COMPLEX PRESSURE
C          P0=SCATTERED COMPLEX PRESSURE IN REGION 1
C             IF BSOUT .NE. 1
C          P1=OUTGOING COMPLEX PRESSURE IF BSOUT .NE. 1
C          P2=INGOING COMPLEX PRESSURE IF BSOUT .NE. 1
C          TL=TRANSMISION LOSS
C          TL0=SCATTERED TRANSMISION LOSS IN REGION 1
C              IF BSOUT .NE. 1
C          TL1=OUTGOING TRANSMISION LOSS IF BSOUT .NE. 1
C          TL2=INGOING TRANSMISION LOSS IF BSOUT .NE. 1
C
C     THE FILES USED ARE AS FOLLOWS.
C
C        INPUT-
C
C          UNIT5=INPUT
C
C        INTERMEDIATE-
C
C          UNIT21=DECOUPLED OUTGOING SOLUTION Y1
C          UNIT2=TOP TWO BLOCK MATRICES (UA & UB) IN UPPER 
C                TRIANGULAR MATRIX U
C          UNIT3=DECOUPLED INGOING SOLUTION Y3
C          UNIT4=UNITARY TRANSFORMATION MATRIX
C          UNIT8=EIGENVALUES AND NORMALIZATION FACTORS
C                ALSO BASIC EIGENVALUES AND NORMALIZED
C                EIGENVECTORS WHEN IMAT=1       
C          UNIT18=UNIQUE SEQUENCE OF EIGENVALUES, NORMALIZATION 
C                 FACTORS, BASIC EIGENVALUES AND NORMALIZED 
C                 EIGENVECTORS SPECIFIED BY IREIG (DIRECT ACCESS).
C          LUSVP=30, TABLES INDEX OF REFRACTION SQUARED IN THE WATER
C                AND COMPLEX INDEX OF REFRACTION SQUARED, DENSITY
C                AND LOGARITHMIC DENSITY GRADIENTS IN THE BOTTOM
C                (WRITTEN IN SUBROUTINE INPUT AND READ IN THE MAIN
C                 PROGRAM COUPLE) 
C
C        OUTPUT-
C
C          UNIT6= FREE FORMAT LOGFILE
C          UNIT9= FORMATTED PRINT FILE
C          UNIT7= REAL AND IMAGINARY PARTS OF COMPLEX PRESSURE
C          UNIT10= TRANSMISION LOSS
C          UNIT11= MODE AMPLITUDES WHEN MAMP .GT. 0
C          UNIT12= OUTGOING TRANSMISION LOSS IF BSOUT .NE. 1
C          UNIT13= INGOING TRANSMISION LOSS IF BSOUT .NE. 1
C          UNIT14= OUTGOING COMPLEX PRESSUE IF BSOUT .NE. 1
C          UNIT15= INGOING COMPLEX PRESSURE IF BSOUT .NE. 1
C          UNIT16= SCATTERED COMPLEX PRESSURE IN REGION 1
C                  IF BSOUT .NE. 1
C          UNIT17= SCATTERED TRANSMISSION LOSS IN REGION 1
C                  BSOUT .NE. 1
C
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(MX=190,MXT2=2*MX,MXMX=MX*MX,MD=400,NRX=2002,
     & MXT2Q=MXT2*MXT2)
C
      DIMENSION IPS(MX),INTH(MX)
      DIMENSION RANG(NRX),DPTH(NRX),RANGE(NRX),DEPTH(NRX)
      DIMENSION IRLIN(NRX),IREIG(NRX)
      DIMENSION DEP(MX),SVP(MX),GRAD(MX),FINT(MX)
      DIMENSION DEPB(MX),SVPB(MX),DBPWLB(MX),RHOB(MX),GRHOB(MX)
      DIMENSION DEPBL(MX),RHOBL(MX),GRHOBL(MX)
      DIMENSION DEPBR(MX),RHOBR(MX),GRHOBR(MX)
      DIMENSION AVGW(NRX),RHOBTOP(NRX)
      DIMENSION ZSD(MD),ZR(MD),TL(MD),TL0(MD),TL1(MD),TL2(MD)
C
      REAL*8 ZZ(MX)
      REAL*8 RSPD(MX,MX),VMAT(MX,MX)
      COMPLEX*16 CINTB(NRX),CGRADB(NRX)
      COMPLEX*16 FNORM(MX)
      COMPLEX*16 EGVL(MX),EGVR(MX),BEGVL(MX),BEGVR(MX)
      COMPLEX*16 CDSRT,BDF,BDF1
      COMPLEX*16 CRL(MX,MX),CLR(MX,MX),BMATL(MX,MX),BMATR(MX,MX)
      COMPLEX*16 R(MXT2,MXT2),T(MXT2,MXT2),U(MXT2,MXT2)
      COMPLEX*16 WRK2(MXT2,MXT2),WRK1(MX,MX),BUFF2(MXT2Q)
      COMPLEX*16 Y1(MX,MX),Y3(MX,MX)
      COMPLEX*16 T3(MX,MX),T4(MX,MX)
      COMPLEX*16 UB(MX,MX),UC(MX,MX),UE(MX,MX)
      COMPLEX*16 C1(MX),A(MX),B(MX),WRK(MX)
      COMPLEX*16 CDF(MX,MD)
      COMPLEX*16 P(MD),P0(MD),P1(MD),P2(MD)
      COMPLEX*16 CI,ONE,FKH,FKVS,FKVW
      COMPLEX*16 CXP,ZERO
      COMPLEX*16 HK01,SUM
      COMPLEX*16 BUFF1(MXMX), ALPHA(MX), FAC(MX), RAT(MX),
     1           BETA(MX),ALPHAR(MX),BETAR(MX),FNORMR(MX),
     2           ALPHAL(MX),BETAL(MX),FNORML(MX),FKHA(MX)
C
      INTEGER BSOUT
      CHARACTER*3 NEWSVP(NRX)
      CHARACTER*10 TITLE(8)
C
      COMMON /ACCUR/ CONST1
      COMMON /BLKEVN/ HB,CW,CB,FKW,FKB,ROHW,ROHB,ATEN
      COMMON /BLKCON/ CI,ONE,ZERO
      COMMON /FLAGS/ IGEOM,NEWSBC
C
      EQUIVALENCE (WRK2(1,1),BUFF2(1))
C     EQUIVALENCE (IRLIN(1),IPS(1))
      EQUIVALENCE (CRL(1,1),UB(1,1)),(CLR(1,1),UC(1,1))
      EQUIVALENCE (CRL(1,1),T3(1,1)),(CLR(1,1),T4(1,1))
      EQUIVALENCE (BMATL(1,1),UE(1,1))
      EQUIVALENCE (R(1,1),U(1,1))
C     EQUIVALENCE (T(1,1),CDF(1,1))
      EQUIVALENCE (BUFF1(1),ALPHA(1)),       (BUFF1(MX+1),FAC(1)),
     1            (BUFF1(2*MX+1),RAT(1)),    (BUFF1(3*MX+1),BETA(1)),
     2            (BUFF1(1),ALPHAR(1)),      (BUFF1(MX+1),BETAR(1)),
     3            (BUFF1(2*MX+1),FNORMR(1)), (BUFF1(3*MX+1),ALPHAL(1)),
     4            (BUFF1(4*MX+1),BETAL(1)),  (BUFF1(5*MX+1),FNORML(1)),
     5            (BUFF1(6*MX+1),FKHA(1)) 
C
C     USE ON FPS 164 TO START TIMER
C
C     CALL SYS$CLTIME
C
C     USE WITH SVS TO START TIMER
C
C      CALL GETTIM(IH,IM,IS,IHU)
C      STARTS=(60.0*IH+IM)*60.0+IS+.01*IHU
C
C     USE WITH LAHEY TO START TIMER
C
C      CALL TIMER(IHU)
C      STARTS=.01*IHU
C
      ISTARTS=time()
C
      PI=2.*ACOS(0.0)
      CI=DCMPLX(0.0,1.0)
      ONE=DCMPLX(1.0,0.0)
      ZERO=DCMPLX(0.0,0.0)
C
C     OPEN INPUT AND OUTPUT FILES
C
      OPEN(5,FILE='COUPLE.DAT',FORM='FORMATTED') 
      OPEN(6,FILE='COUPLE.LOG',FORM='FORMATTED') 
      OPEN(9,FILE='COUPLE.PRT',FORM='FORMATTED')
      LUSVP=30
      OPEN(LUSVP,STATUS='SCRATCH',FILE='LUSVP.DAT',FORM='FORMATTED')
C
      WRITE(*,10)
      WRITE(9,10)
   10 FORMAT(1X,/,'  ********  ',/,
     & ' THIS PC VERSION OF COUPLE IS RANGE DEPENDENT ',/,
     & ' HAS PLANE OR CYLINDRICAL GEOMETRY OPTION ',/,
     & ' AND STANDARD OR NEW INITIAL SOURCE BOUNDARY CONDITIONS ',/,
     & '  ********  ',/)
C
C     GET INPUT DATA 
C
      CALL INPUT(TITLE,FREQ,M,N,IPRT,IOUT,BSOUT,NSDEP,NDEP,
     1           NATEN,NDPTH,IEIG,IMAT,RANGE,DEPTH,NREG,MAMP,
     2           ZS,ZSD,RMIN,RMAX,RINC,ZR,NDINC,MD,NRX,
     3           DEP,SVP,GRAD,FINT,RANG,DPTH,IRLIN,IREIG,MX,
     4           DEPB,SVPB,DBPWLB,CGRADB,CINTB,RHOB,GRHOB,
     5           NEWSVP,LUSVP)
C
C
      H1=DEPTH(1)
      RKM1=RANGE(1)
C
      IF(IMAT .EQ. 1)   THEN 
      REWIND LUSVP
      READ(LUSVP,*) NP,NPB
      DO 11 KK=1,NP
      READ(LUSVP,*) DEP(KK),FINT(KK),GRAD(KK)
   11 CONTINUE
      DO 12 KK=1,NPB
      READ(LUSVP,*) DEPB(KK),CINTB(KK),CGRADB(KK),RHOB(KK),GRHOB(KK)
   12 CONTINUE
      CALL AVGIN(DEP,GRAD,FINT,NP,MX,H1,AVG,LH)
      CALL BOTDEPS(H1,DEPB,CGRADB,CINTB,NPB,AVGB,LHB,MX)
      CW=CB/SQRT(AVG)
      FKW=2.0*PI*FREQ/CW
      AVGW(1)=AVG
      RHOBTOP(1)=RHOB(1)
      ROHB=RHOB(1)
C
C     CHECK FOR MODE CUTOFF, AVOID IF NECESSARY BY INCREASING
C     THE WATER DEPTH BY 1/1000 OF A WAVE LENGTH, AND REDEFINE
C     THE AVERAGE INDEX OF REFRACTION IN THE WATER AND THE DEPTHS
C     AND INTERCEPTS IN THE BOTTOM
C
        CUTOFF=RCF(H1,1.0D0)
        IF(DABS(CUTOFF) .LT. 1.0D-12) THEN
        DEPINC=2.0D0*PI/(1000.0D0*FKW)
        WRITE(*,*) ' *** WARNING *** '
        WRITE(*,*) 'MODE CUTOFF AT RANGE=',RANGE(1),' KM'
        WRITE(*,*) 'DEPTH INCREASED BY ',DEPINC,' M'
        WRITE(9,*) ' '
        WRITE(9,*) ' *** WARNING *** '
        WRITE(9,*) 'MODE CUTOFF AT RANGE=',RANGE(1),' KM'
        WRITE(9,*) 'DEPTH INCREASED BY ',DEPINC,' M'
        DEPTH(1)=DEPTH(1)+DEPINC
        H1=DEPTH(1)
        CALL AVGIN(DEP,GRAD,FINT,NP,MX,H1,AVG,LH)
        CALL BOTDEPS(H1,DEPB,CGRADB,CINTB,NPB,AVGB,LHB,MX)
        CW=CB/SQRT(AVG)
        FKW=2.0*PI*FREQ/CW
        AVGW(1)=AVG
        END IF
C
      NPBL=NPB
      LHBL=LHB
      DO 13 KK=1,NPBL
      DEPBL(KK)=DEPB(KK)
      RHOBL(KK)=RHOB(KK)
      GRHOBL(KK)=GRHOB(KK)
   13 CONTINUE
C
C      WRITE(9,*) 'NPBL,LHBL= ',NPBL,LHBL
C      WRITE(*,*) 'NPBL,LHBL= ',NPBL,LHBL
C
      END IF
CCC
      FKWL=FKW
      ROHBL=ROHB
CCC
C         WRITE(9,*) 'CW,CB= ',CW,CB
C         WRITE(9,*) 'NP,LH= ',NP,LH
C         DO 15 KK=1,NP
C         WRITE(9,14) DEP(KK),FINT(KK)+GRAD(KK)*DEP(KK),GRAD(KK)
C   14    FORMAT(1X,F10.4,2X,F10.4,2X,E12.4)
C   15    CONTINUE
C         WRITE(9,14) H1,FINT(LH)+GRAD(LH)*H1,GRAD(LH)
C         WRITE(9,*) 'NPB,LHB= ',NPB,LHB
C         DO 17 KK=1,NPB
C         WRITE(9,16) DEPB(KK),CINTB(KK)+CGRADB(KK)*DEPB(KK),CGRADB(KK)
C   16    FORMAT(1X,F10.4,2X,F10.4,2X,E12.4,2X,E12.4,2X,E12.4)
C   17    CONTINUE
C         WRITE(9,16) HB,CINTB(LHB)+CGRADB(LHB)*HB,CGRADB(LHB)
CC
C     OPEN SCRATCH AND OTHER FILES
C
      CALL FILE(MAMP,BSOUT,MX,MXT2)
C
C     WRITE HEADERS ON THE OUTPUT FILES
C
      CALL HEADER(TITLE,FREQ,ZS,NDEP,NDINC,ZR,NREG,RANGE,DEPTH,MD,
     1            NRX,MAMP,BSOUT)
C
C     FIND THE FIRST SET OF EIGENVALUES AND NORMALIZATION
C     FACTORS
C
      IF(IMAT .EQ. 0)   THEN
      CALL SEARCH(1,RKM1,H1,H1,M,ZZ,EGVL,FNORM,
     &            EGVL,NATEN,NDPTH,IEIG,IPRT,MX,WRK)
      ELSE
      CALL GALRKN(H1,M,ZZ,EGVL,FNORM,
     1            BEGVL,BMATL,VMAT,MX,WRK,WRK1,RSPD,CLR,CRL,
     2            IPS,INTH,DEP,GRAD,FINT,NP,AVG,LH,
     3            DEPB,RHOB,GRHOB,CGRADB,CINTB,NPB,AVGB,LHB)
      END IF
C
      IF(IPRT .LE. 0) GO TO 50
C
      WRITE(9,20)
   20 FORMAT(/)
      WRITE(*,25) 1,RANGE(1),DEPTH(1)
      WRITE(9,25) 1,RANGE(1),DEPTH(1)
   25 FORMAT(1X,26HTHE EIGENVALUES IN REGION ,I5,11H ENDING AT ,F10.3,
     1       16H  KM WITH DEPTH ,F10.2,18H  M ARE AS FOLLOWS//)
      WRITE(9,30)
   30 FORMAT(2X,3HNO.,6X,10HEIGENVALUE,9X,20HNORMALIZATION FACTOR,
     1       8X,15HHORIZ. WAVE NO.,8X,22HVERT. WAVE NO. IN BOT.,
     2        4X,21HVERT. WAVE NO. IN H2O/)
C
      DO 40 I=1,M
      FKH=FKB*CDSRT(ONE-EGVL(I))
      IF(REAL(FKH) .LT. 0.0)   THEN
      WRITE(9,*) ' ERROR IN FIRST REGION :'
      WRITE(9,*) ' NEGATIVE WAVE NUMBER FOR MODE NO. ',I
      WRITE(9,*) ' WAVE NUMBER : ',FKH
      END IF
      FKVW=FKB*CDSRT((FKW/FKB)**2+EGVL(I)-ONE)
      FKVS=FKB*CDSRT(EGVL(I)+CI*ATEN)
      WRITE(9,35) I,EGVL(I),FNORM(I),FKH,FKVS,FKVW
   35 FORMAT(1X,I3,1P,2(1X,2E11.3),3(1X,2E13.5))
   40 CONTINUE
C
C
   50 CONTINUE
C
      IF(IREIG(1) .GT. 0) THEN
C
C     STORE UNIQUE SEQUENCE OF EIGENVALUES ECT., FOR RETRIEVAL.
C
        IF(IMAT .GE. 1)  THEN
        ISTART=(IREIG(1)-1)*(M+3)
        WRITE(18,REC=ISTART+1) (BEGVL(I),I=1,M)
        DO 52 I=1,M
        WRITE(18,REC=ISTART+1+I) (BMATL(II,I),II=1,M)
   52   CONTINUE
        WRITE(18,REC=ISTART+M+2) (EGVL(I),I=1,M)
        WRITE(18,REC=ISTART+M+3) (FNORM(I),I=1,M)
        ELSE
        ISTART=(IREIG(1)-1)*2
        WRITE(18,REC=ISTART+1) (EGVL(I),I=1,M)
        WRITE(18,REC=ISTART+2) (FNORM(I),I=1,M)
        END IF
C
      END IF
C
      IF(BSOUT .EQ. 1) GO TO 55
      IF(IMAT .GE. 1) WRITE(8) (BEGVL(I),I=1,M)
      IF(IMAT .GE. 1) WRITE(8) ((BMATL(II,I),II=1,M),I=1,M)
      WRITE(8) (EGVL(I),I=1,M)
      WRITE(8) (FNORM(I),I=1,M)
   55 CONTINUE
C
      RNG1=1000.*RANGE(1)
      DO 60 I=1,M
      C1(I)=ZERO
   60 CONTINUE
C
C     STORE THE VECTOR d IN C1. IT IS ASSUMED THAT ALL THE
C     SOURCES ARE IN THE WATER LAYER WITH DENSITY ROHW
C
      IF(IMAT .LT. 1) THEN
C
C     USE ONLY BASIC EIGEFUNCTIONS
C
      DO 70 I=1,M
C
      FKH=FKB*CDSRT(ONE-EGVL(I))
      FKVW=FKB*CDSRT((FKW/FKB)**2+EGVL(I)-ONE)
      FKVS=FKB*CDSRT(EGVL(I)+CI*ATEN)
      IF(IGEOM .EQ. 1) THEN
      FAC(I)=CDEXP(CI*RNG1*FKH)*CI/(ROHW*2.0*FKH)
      ELSE
      FAC(I)=HK01(RNG1*FKH)*CI/(ROHW*4.0)
      END IF
C
C     SUM OVER SOURCE DEPTHS
C
      DO 65 J=1,NSDEP
      BDF1=BDF(ZSD(J),H1,FKVW,FKVS,FNORM(I))
      C1(I)=C1(I)+BDF1*FAC(I)
   65 CONTINUE
      C1(I)=C1(I)/NSDEP
C
   70 CONTINUE
C
      ELSE
C
C     USE DEPTH FUNCTIONS CONSTRUCTED FROM THE BASIC
C     EIGENFUNCTIONS AND THE MATRIX EIGENVECTORS IN
C     THE ARRAY BMATL 
C
      DO 80 I=1,M
      ALPHAL(I)=FKB*CDSRT((FKW/FKB)**2+BEGVL(I)-ONE)
      BETAL(I)=FKB*CDSRT(BEGVL(I))
      FKH=FKB*CDSRT(ONE-EGVL(I))
      IF(IGEOM .EQ. 1) THEN
      FAC(I)=CDEXP(CI*RNG1*FKH)*CI/(ROHW*2.0*FKH)
      ELSE
      FAC(I)=HK01(RNG1*FKH)*CI/(ROHW*4.0)
      END IF
   80 CONTINUE
C
      DO 100 I=1,M
C
C     SUM OVER SOURCE DEPTHS
C
      DO 90 J=1,NSDEP
C
C     SUM OVER BASIC EIGENFUNCTIONS
C
      SUM=ZERO
      DO 85 II=1,M
      SUM=SUM+BMATL(II,I)*BDF(ZSD(J),H1,ALPHAL(II),
     1                        BETAL(II),FNORM(II))
   85 CONTINUE
C
      C1(I)=C1(I)+SUM*FAC(I)
   90 CONTINUE
      C1(I)=C1(I)/NSDEP
C
  100 CONTINUE
C
      END IF
C
      IF(BSOUT .NE. 1) GO TO 200
C
C     FIND THE OUTGOING COEFFICIENTS IN THE FIRST REGION (BSOUT=1)
C
      WRITE(9,105)
  105 FORMAT(//)
      WRITE(*,110) BSOUT
      WRITE(9,110) BSOUT
  110 FORMAT(1X,'BACKSCATTER WILL BE NEGLECTED AND ONLY AN',
     1       ' OUTGOING SOLUTION IS FOUND,','BSOUT= ',I5,/)   
C    
      IF (IOUT .EQ. 1) THEN  
      WRITE(*,112) IOUT
      WRITE(9,112) IOUT
  112 FORMAT(1X,'ONE-WAY SOLUTION IS USED: IOUT= ',I5,/)   
      ELSE IF (IOUT .EQ. -1) THEN
      WRITE(*,113) IOUT
      WRITE(9,113) IOUT
  113 FORMAT(1X,'SINGLE SCATTER APPROXIMATION IS USED: IOUT= ',I5,///)
      END IF
C
      DO 120 I=1,M
      A(I)=C1(I)
      B(I)=ZERO
  120 CONTINUE
C
C      UNORMALIZED IN AND OUTGOING COEFICIENTS ARE PRINTED
C      BELOW AND IN ALL SUBSQUENT PRINTS.
C
C
      IF(IPRT .LE. 0) GO TO 150
C
      WRITE(*,130) 1,RANGE(1),DEPTH(1)
      WRITE(9,130) 1,RANGE(1),DEPTH(1)
  130 FORMAT(1X,27HTHE COEFFICIENTS IN REGION ,I5,11H ENDING AT ,F10.3,
     1       16H  KM WITH DEPTH ,F10.2,18H  M ARE AS FOLLOWS///
     2       10H  MODE NO.,11X,8HOUTGOING,12X,13HOUTGOING (DB)//)
C
      DO 140 I=1,M
C
      ATL=200.
      AABS=CDABS(A(I))
      IF(AABS .GT. 1.0D-10) ATL=-20.*DLOG10(AABS)
C
      WRITE(9,135) I,A(I),ATL
  135 FORMAT(1X,I5,5X,1X,2E12.4,5X,F8.2)
  140 CONTINUE
C
  150 CONTINUE
C
C     THIS INITIALIZES THE RANGE LOOP FOR THE PRESSURE CALCULATION
C     (BSOUT=1)
C
      RKM=RMIN
C
      IF(RANGE(1) .LT. RMIN) GO TO 200
C
      R2KM=RANGE(1)
      R1KM=RANGE(1)
C
      JBOT=0
C
C     COMPUTE THE DEPTH FUNCTIONS AND PRESSURE IN THE FIRST REGION
C     (BSOUT=1)
C
      CALL CDFUN(M,EGVL,FNORM,H1,ZR,NDEP,CDF,MX,MD,
     1           IMAT,BEGVL,BMATL,ALPHA,FAC,RAT,BETA)
C
C
      CALL CPRES(M,RINC,RMAX,R1KM,R2KM,RKM,EGVL,A,B,
     1            NDEP,NDINC,CDF,JBOT,MX,MD,P,P0,P1,P2,
     2            TL,TL0,TL1,TL2,BSOUT,MAMP)
C
C
  200 CONTINUE
C
C
C ******** START OUTGOING LOOP TO FIND EIGENVALUES,COMPUTE *********
C ******** PRESSURE (BSOUT=1) AND GENERATE R, T, U AND Y1   *********
C
      NREM1=NREG-1
      DO 1000 J=1,NREM1
      H1=DEPTH(J)
      H2=DEPTH(J+1)
      RKMJP1=RANGE(J+1)
C
C     FIND OR RETRIEVE EIGENVALUES AND NORMALIZATION FACTORS
C
      IF(IMAT .EQ. 0)   THEN
        IF(IREIG(J+1) .LT. 0) THEN
        ISTART=(-IREIG(J+1)-1)*2
        READ(18,REC=ISTART+1) (EGVR(I),I=1,M)
        READ(18,REC=ISTART+2) (FNORM(I),I=1,M)
        ELSE
        CALL SEARCH(J+1,RKMJP1,H1,H2,M,ZZ,EGVL,FNORM,
     &              EGVR,NATEN,NDPTH,IEIG,IPRT,MX,WRK)
        END IF
      ELSE
       IF(NEWSVP(J+1) .EQ. 'NEW')   THEN
         READ(LUSVP,*) NP,NPB
         DO 215 KK=1,NP
         READ(LUSVP,*) DEP(KK),FINT(KK),GRAD(KK)
  215    CONTINUE
         DO 216 KK=1,NPB
         READ(LUSVP,*) DEPB(KK),CINTB(KK),CGRADB(KK),RHOB(KK),GRHOB(KK)
  216    CONTINUE
       END IF
       CALL AVGIN(DEP,GRAD,FINT,NP,MX,H2,AVG,LH)
       CALL BOTDEPS(H2,DEPB,CGRADB,CINTB,NPB,AVGB,LHB,MX)
       CW=CB/SQRT(AVG)
       FKW=2.0*PI*FREQ/CW
       AVGW(J+1)=AVG
       RHOBTOP(J+1)=RHOB(1)
       ROHB=RHOB(1)
C
C     CHECK FOR MODE CUTOFF, AVOID IF NECESSARY BY INCREASING
C     THE WATER DEPTH BY 1/1000 OF A WAVE LENGTH, AND REDEFINE
C     THE AVERAGE INDEX OF REFRACTION IN THE WATER AND THE DEPTHS
C     AND INTERCEPTS IN THE BOTTOM
C
        CUTOFF=RCF(H2,1.0D0)
        IF(DABS(CUTOFF) .LT. 1.0D-12) THEN
        DEPINC=2.0D0*PI/(1000.0D0*FKW)
        WRITE(*,*) ' *** WARNING *** '
        WRITE(*,*) 'MODE CUTOFF AT RANGE=',RANGE(J+1),' KM'
        WRITE(*,*) 'DEPTH INCREASED BY ',DEPINC,' M'
        WRITE(9,*) ' '
        WRITE(9,*) ' *** WARNING *** '
        WRITE(9,*) 'MODE CUTOFF AT RANGE=',RANGE(J+1),' KM'
        WRITE(9,*) 'DEPTH INCREASED BY ',DEPINC,' M'
        DEPTH(J+1)=DEPTH(J+1)+DEPINC
        H2=DEPTH(J+1)
        CALL AVGIN(DEP,GRAD,FINT,NP,MX,H2,AVG,LH)
        CALL BOTDEPS(H2,DEPB,CGRADB,CINTB,NPB,AVGB,LHB,MX)
        CW=CB/SQRT(AVG)
        FKW=2.0*PI*FREQ/CW
        AVGW(J+1)=AVG
        END IF
CC       
C         WRITE(9,*) 'CW,CB= ',CW,CB
C         WRITE(9,*) 'NP,LH= ',NP,LH
C         DO 218 KK=1,NP
C         WRITE(9,217) DEP(KK),FINT(KK)+GRAD(KK)*DEP(KK),GRAD(KK)
C  217    FORMAT(1X,F10.4,2X,F10.4,2X,E12.4)
C  218    CONTINUE
C         WRITE(9,217) H2,FINT(LH)+GRAD(LH)*H2,GRAD(LH)
C         WRITE(9,*) 'NPB,LHB= ',NPB,LHB
C         DO 229 KK=1,NPB
C         WRITE(9,227) DEPB(KK),CINTB(KK)+CGRADB(KK)*DEPB(KK),CGRADB(KK)
C  227    FORMAT(1X,F10.4,2X,F10.4,2X,E12.4,2X,E12.4,2X,E12.4)
C  229    CONTINUE
C         WRITE(9,227) HB,CINTB(LHB)+CGRADB(LHB)*HB,CGRADB(LHB)
CC
        IF(IREIG(J+1) .LT. 0) THEN
        ISTART=(-IREIG(J+1)-1)*(M+3)
        READ(18,REC=ISTART+1) (BEGVR(I),I=1,M)
        DO 242 I=1,M
        READ(18,REC=ISTART+1+I) (BMATR(II,I),II=1,M)
  242   CONTINUE
        READ(18,REC=ISTART+M+2) (EGVR(I),I=1,M)
        READ(18,REC=ISTART+M+3) (FNORM(I),I=1,M)
        ELSE
        CALL GALRKN(H2,M,ZZ,EGVR,FNORM,
     1             BEGVR,BMATR,VMAT,MX,WRK,WRK1,RSPD,CLR,CRL,
     2             IPS,INTH,DEP,GRAD,FINT,NP,AVG,LH,
     3             DEPB,RHOB,GRHOB,CGRADB,CINTB,NPB,AVGB,LHB)
        END IF
C
      NPBR=NPB
      LHBR=LHB
      DO 253 KK=1,NPBR
      DEPBR(KK)=DEPB(KK)
      RHOBR(KK)=RHOB(KK)
      GRHOBR(KK)=GRHOB(KK)
  253 CONTINUE
C
      END IF
C
C      MOVED FOLLOWING LINE HERE ON 8/20/91. IT WAS ONLY DONE
C      IN THE GALERKIN METHOD BUT WAS ALSO NEEDED IN THE SEARCH 
C      METHOD. R. B. EVANS (SAIC, NEW LONDON, CT)
CCC
      FKWR=FKW
      ROHBR=ROHB
CCC
C
      DO 280 I=1,M
      FKH=FKB*CDSRT(ONE-EGVR(I))
      IF(REAL(FKH) .LT. 0.0)   THEN
      WRITE(9,*) ' ERROR IN REGION NO :',J+1
      WRITE(9,*) ' NEGATIVE WAVE NUMBER FOR MODE NO. ',I
      WRITE(9,*) ' WAVE NUMBER : ',FKH
      END IF
  280 CONTINUE
      JP1=J+1
      IF(IPRT .LE. 0) GO TO 330
      IF(JP1 .EQ. NREG) GO TO 290
      IF(MOD(JP1,IPRT) .NE. 0) GO TO 330
  290 CONTINUE
C
      WRITE(9,20)
      WRITE(*,25) JP1,RANGE(JP1),DEPTH(JP1)
      WRITE(9,25) JP1,RANGE(JP1),DEPTH(JP1)
      WRITE(9,30)
C
      DO 300 I=1,M
      FKH=FKB*CDSRT(ONE-EGVR(I))
      FKVS=FKB*CDSRT(EGVR(I)+CI*ATEN)
      FKVW=FKB*CDSRT((FKW/FKB)**2+EGVR(I)-ONE)
      WRITE(9,35) I,EGVR(I),FNORM(I),FKH,FKVS,FKVW
  300 CONTINUE
C
  330 CONTINUE
C
C     STORE UNIQUE SEQUENCE OF EIGENVALUES ECT., FOR RETRIEVAL
C
      IF(IREIG(J+1) .GT. 0) THEN
        IF(IMAT .GE. 1) THEN
        ISTART=(IREIG(J+1)-1)*(M+3)
        WRITE(18,REC=ISTART+1) (BEGVR(I),I=1,M)
        DO 342 I=1,M
        WRITE(18,REC=ISTART+1+I) (BMATR(II,I),II=1,M)
  342   CONTINUE
        WRITE(18,REC=ISTART+M+2) (EGVR(I),I=1,M)
        WRITE(18,REC=ISTART+M+3) (FNORM(I),I=1,M)
        ELSE
        ISTART=(IREIG(J+1)-1)*2
        WRITE(18,REC=ISTART+1) (EGVR(I),I=1,M)
        WRITE(18,REC=ISTART+2) (FNORM(I),I=1,M)
        END IF
      END IF
C
      IF(BSOUT .EQ. 1) GO TO 350
      IF(IMAT .GE. 1) WRITE(8) (BEGVR(I),I=1,M)
      IF(IMAT .GE. 1) WRITE(8) ((BMATR(II,I),II=1,M),I=1,M)
      WRITE(8) (EGVR(I),I=1,M)
      WRITE(8) (FNORM(I),I=1,M)
  350 CONTINUE
C
C     COMPUTE THE CROSS COUPLING MATRICES
C
      IF(IMAT .EQ. 0)   THEN
      CALL CRSCUP(H1,H2,M,EGVL,EGVR,CLR,CRL,MX,
     &            IMAT,BMATL,BMATR,WRK1,
     &            ALPHAR,BETAR,FNORMR,
     &            ALPHAL,BETAL,FNORML,FKWL,FKWR,ROHBL,ROHBR)
C      
C
      ELSE
      CALL CRSCUPD(H1,H2,M,BEGVL,BEGVR,CLR,CRL,MX,
     1            BMATL,BMATR,WRK1,ALPHAR,BETAR,FNORMR,
     2            ALPHAL,BETAL,FNORML,FKWL,FKWR,ROHBL,ROHBR,
     3            NPBL,DEPBL,RHOBL,GRHOBL,LHBL,
     4            NPBR,DEPBR,RHOBR,GRHOBR,LHBR)
C      
C
      END IF
C
C      WRITE(9,*) 'ROHBL,ROHBR= ',ROHBL,ROHBR
C
      IF(IPRT .EQ. 1) CALL PRTMAT(M,M,'CLR ',CLR,MX) 
      IF(IPRT .EQ. 1) CALL PRTMAT(M,M,'CRL ',CRL,MX) 
C
      RNG1=RANGE(J)*1000.
      IF(J .GT. 1)  RNG1=RANGE(J-1)*1000.
      RNG2=RANGE(J)*1000.
      MT2=M*2
C
C     COMPUTE THE PROPAGATOR MATRIX R
C
      CALL PROPM(M,RNG1,RNG2,EGVL,EGVR,CLR,CRL,R,MX,MXT2,IOUT)
C
C
      IF(BSOUT .NE. 1) GO TO 500
C
C     GENERATE THE OUTGOING COEFFICIENTS (BSOUT=1)
C
      DO 370 II=1,M
      SUM=ZERO
      DO 360 I=1,M
      SUM=SUM+R(II+M,I+M)*A(I)
  360 CONTINUE
      B(II)=SUM
  370 CONTINUE
C
      DO 380 I=1,M
      A(I)=B(I)
      B(I)=ZERO
  380 CONTINUE  
C
      IF(IPRT .LE. 0) GO TO 400
      IF(JP1 .EQ. NREG) GO TO 385
      IF(MOD(JP1,IPRT) .NE. 0) GO TO 400
  385 CONTINUE
C
      WRITE(9,105)
      WRITE(*,130) JP1,RANGE(JP1),DEPTH(JP1)
      WRITE(9,130) JP1,RANGE(JP1),DEPTH(JP1)
C
      DO 390 I=1,M
C
      ATL=200.
      AABS=CDABS(A(I))
      IF(AABS .GT. 1.0D-10) ATL=-20.*DLOG10(AABS)
      WRITE(9,135) I,A(I),ATL
  390 CONTINUE
C
  400 CONTINUE
C
      IF((RANGE(J+1) .LT. RMIN) .OR. (RANGE(J) .GE. RMAX)) GO TO 900
      IF(RKM .GT. RANGE(J+1)) GO TO 900
C
      R1KM=RANGE(J)
      R2KM=RANGE(J+1)
      JBOT=0
C
C     CALCULATE DEPTH FUNCTION AND PRESSURE (BSOUT=1)
C
C
      CALL CDFUN(M,EGVR,FNORM,H2,ZR,NDEP,CDF,MX,MD,
     1           IMAT,BEGVR,BMATR,ALPHA,FAC,RAT,BETA)
C
C
      CALL CPRES(M,RINC,RMAX,R1KM,R2KM,RKM,EGVR,A,B,
     1            NDEP,NDINC,CDF,JBOT,MX,MD,P,P0,P1,P2,
     2            TL,TL0,TL1,TL2,BSOUT,MAMP)
C
C
      GO TO 900
C 
  500 CONTINUE  
C
      IF(J .GT. 1) GO TO 580
C
      DO 550 I=1,MT2
      DO 540 II=1,MT2
      WRK2(II,I)=R(II,I)
  540 CONTINUE
  550 CONTINUE
C
      GO TO 600
C
  580 CONTINUE
C
C     STORE THE PRODUCT R*T IN WRK2 (WRK2=R*T)
C
      CALL MMULT(R,T,WRK2,MT2,MXT2) 
C
C
  600 CONTINUE
C
C     DO MODIFIED GRAM-SCHMIDT DECOMPOSITION OF (WRK2=T*U) 
C
C
      CALL MGS(WRK2,T,U,MT2,MXT2)
C
C     THE MATRICES T AND  U ARE PARTITIONED AS FOLLOWS:
C
C         T4 T3              UB  UC
C     T =         AND   U = 
C         T2 T1              0   UE
C
C     DEFINE UE
C
      DO 650 I=1,M
      DO 640 II=1,M
      UE(II,I)=U(M+II,M+I)
  640 CONTINUE
  650 CONTINUE
C
C     STORE UB, UC AND T
C
      WRITE(2) ((U(II,I),II=1,I),I=1,M)
      WRITE(2) ((U(II,M+I),II=1,M),I=1,M)
      WRITE(4) ((T(II,I),II=1,MT2),I=1,MT2)
C
      IF(J .GT. 1) GO TO 820
C
      DO 800 I=1,M
      DO 700 II=1,M
C
      Y1(II,I)=UE(II,I)
C
  700 CONTINUE
  800 CONTINUE
C
      GO TO 850
C
  820 CONTINUE
C
C     ADVANCE Y1 (Y1=UE*Y1)
C
      CALL MMULT(UE,Y1,WRK1,M,MX) 
C
C
      DO 840 I=1,M
      DO 830 II=1,I
      Y1(II,I)=WRK1(II,I)
  830 CONTINUE
  840 CONTINUE
C
  850 CONTINUE
C
      WRITE(21) ((Y1(II,I),II=1,I),I=1,M)
C
  900 CONTINUE
C
      DO 920 I=1,M
      BEGVL(I)=BEGVR(I)
      DO 910 II=1,M
      BMATL(II,I)=BMATR(II,I)
  910 CONTINUE
  920 CONTINUE
C
      DO 950 I=1,M
      EGVL(I)=EGVR(I)
  950 CONTINUE
C
      IF(IMAT .NE. 0) THEN
      NPBL=NPBR
      LHBL=LHBR
      DO 953 KK=1,NPBL
      DEPBL(KK)=DEPBR(KK)
      RHOBL(KK)=RHOBR(KK)
      GRHOBL(KK)=GRHOBR(KK)
  953 CONTINUE
      END IF
CCC
      FKWL=FKWR
      ROHBL=ROHBR
CCC
C 
 1000 CONTINUE
C
C ***** IF BSOUT=1,THE OUTGOING SOLUTION IS HAS BEEN OBTAINED *****
C *****               THE CALCULATION IN OVER                *****
C
      IF(BSOUT .EQ. 1) GO TO 8000
C
C     THE MATRIX Y1, IN THE LAST REGION, IS NOT READ FROM UNIT21
C     IN THE LOOP TO GENERATE Y3, BUT IS USED IMEDIATELY IN THE 
C     SOLVE FOR Y3
C
      ENDFILE 21
      BACKSPACE 21
      BACKSPACE 21
      ENDFILE 2
      BACKSPACE 2
C
      ENDFILE 4
      REWIND 4
C
      ENDFILE 8
      REWIND 8
C
C     SET UP MATRIX EQUATION FOR Y3 IN THE LAST REGION. T3 IS 
C     MULTIPLIED BY -1 FOR CONVENIENCE. 
C
      DO 1200 I=1,M
      DO 1100 II=1,M
      T4(II,I)=T(II,I)
      T3(II,I)=-T(II,M+I)
 1100 CONTINUE
 1200 CONTINUE
C
C************ SOLVE MATRIX EQUATION FOR Y3 ***************
C
      CALL MMULT(T3,Y1,Y3,M,MX)
      IF(IPRT .GE. 1) CALL PRTMAT(M,M,'T4  ',T4,MX) 
      CALL CLINEQS(M,MX,M,MX,T4,Y3,WRK1,WRK,EGVR,IPS,IER)
C
      IFLAG=3
C
C     ABORT NOT COMMENTED OUT
C
      IF(IER .NE. 0) CALL ABORT(IFLAG)
C
      WRITE(9,105)
      IF(IER .EQ. 1)   THEN
      WRITE(*,1261) 
      WRITE(9,1261) 
      ELSE
      WRITE(*,1260) IER
      WRITE(9,1260) IER
      END IF
 1260 FORMAT(1X,4HIER=,I5///)
 1261 FORMAT(1X,///,
     & ' ***************************************************** ',
     & //,' WARNING :  FAILURE TO IMPROVE BY IMPRUV IN CLINEQS ',/, 
     &    '            AS ITERATION DID NOT CONVERGE IN SOLVE ',//,
     & ' *********************************************** ',///)
C
      WRITE(3) ((Y3(II,I),II=1,M),I=1,M)  
C
C************* START IN GOING LOOP TO GENERATE Y3 **************
C
C
      DO 2000 JJ=1,NREM1
      J=NREM1-JJ+1
C
      BACKSPACE 2
      BACKSPACE 2
      READ(2) ((UB(II,I),II=1,I),I=1,M)
      READ(2) ((UC(II,I),II=1,M),I=1,M)
      BACKSPACE 2
      BACKSPACE 2
C
      IF(J .GT. 1) GO TO 1500
C
      DO 1400 I=1,M
      DO 1300 II=1,M
      WRK1(II,I)=Y3(II,I)-UC(II,I)
 1300 CONTINUE
 1400 CONTINUE
      GO TO 1900
C
 1500 CONTINUE
C
      BACKSPACE 21
      READ(21) ((Y1(II,I),II=1,I),I=1,M)
      BACKSPACE 21
C
C
      CALL MMULT(UC,Y1,WRK1,M,MX)
C
C
      DO 1800 I=1,M
      DO 1700 II=1,M
      WRK1(II,I)=Y3(II,I)-WRK1(II,I)
 1700 CONTINUE
 1800 CONTINUE
C
C
 1900 CONTINUE
C
C
      CALL BACKSUB(UB,Y3,WRK1,M,MX)
C
C
      IF(J .GT. 1) WRITE(3) ((Y3(II,I),II=1,M),I=1,M)
C
 2000 CONTINUE
C
      ENDFILE 3
      BACKSPACE 3
C 
      REWIND 21
C
C     THE LAST Y3 (J=1) COMPUTED IS THE REFLECTION MATRIX FOR
C     REGION 2 THROUGH NREG. THE MATRIX Y3 IS USED IN THE
C     SET UP OF THE MATRIX EQUATION FOR THE VECTOR C1 THAT 
C     SATISFIES THE SOURCE CONDITION IN THE FIRST REGION.
C     THE MATRIX Y3 IS ALSO USED IN TO GENERATE THE INCOMING
C     COEFFICIENTS IN THE FIRST REGION. 
C
      IF(IMAT .GE. 1) READ(8) (BEGVL(I),I=1,M)
      IF(IMAT .GE. 1) READ(8) ((BMATL(II,I),II=1,M),I=1,M)
      READ(8) (EGVL(I),I=1,M)
      READ(8) (FNORM(I),I=1,M)
C
      H1=DEPTH(1)
C
      IF(IMAT .NE. 0)   THEN
      AVG=AVGW(1)
      ROHB=RHOBTOP(1)
      CW=CB/SQRT(AVG)
      FKW=2.0*PI*FREQ/CW
      END IF
C
      RNG1=RANGE(1)*1000.
C
C     STORE THE MATRIX I-D*Y3 IN T4
C
      IF(NEWSBC .GT. 0)   THEN
      WRITE(*,2100)
      WRITE(9,2100)
 2100 FORMAT(1X,/,
     & ' *** NEW INITIAL SOURCE BOUNDARY CONDITIONS ARE COMPUTED ***',/)
      END IF
C
      DO 2500 II=1,M
C
      FKH=FKB*CDSRT(ONE-EGVL(II))
      IF(-DIMAG(2.*FKH*RNG1) .LT. -88.) GO TO 2150
      IF(IGEOM .EQ. 1) THEN
      CXP=CDEXP(2.*CI*FKH*RNG1)
      ELSE
      CXP=-CI*CDEXP(2.*CI*FKH*RNG1)
      END IF
      GO TO 2200
 2150 CONTINUE
      CXP=ZERO
 2200 CONTINUE
      IF(NEWSBC .EQ. 0)   THEN
      DO 2400 I=1,M
      T4(II,I)=-CXP*Y3(II,I)
      IF(I .EQ. II) T4(II,I)=ONE+T4(II,I)
 2400 CONTINUE
      ELSE
      DO 2450 I=1,M
      T4(II,I)=ZERO
      IF(I .EQ. II) T4(II,I)=ONE
 2450 CONTINUE
      END IF
C
 2500 CONTINUE
C
C 
      IF(IPRT .GE. 1) CALL PRTMAT(M,M,'REFL',Y3,MX) 
      IF(IPRT .GE. 1) CALL PRTMAT(M,M,'I-DY',T4,MX) 
C
C ********** SOLVE MATRIX EQUATION FOR C1 *****************
C
      CALL CLINEQS(M,MX,1,1,T4,C1,WRK1,WRK,EGVR,IPS,IER)
C     
C
      IFLAG=3
C
C     ABORT NOT COMMENTED OUT
C
      IF(IER .NE. 0) CALL ABORT(IFLAG)
C
      WRITE(9,105)
      IF(IER .EQ. 1)   THEN
      WRITE(*,1261) 
      WRITE(9,1261) 
      ELSE
      WRITE(*,1260) IER
      WRITE(9,1260) IER
      END IF
C
      WRITE(*,3690) 1,RANGE(1),DEPTH(1)
      WRITE(9,3690) 1,RANGE(1),DEPTH(1)
 3690 FORMAT(1X,27HTHE COEFFICIENTS IN REGION ,I5,11H ENDING AT ,F10.3,
     1       16H  KM WITH DEPTH ,F10.2,18H  M ARE AS FOLLOWS///
     2       10H  MODE NO.,11X,8HOUTGOING,14X,11HBACKSCATTER,
     3       12X,13HOUTGOING (DB),4X,16HBACKSCATTER (DB),3X,
     4       10HTOTAL (DB)//)
C
C     GENERATE THE COEFFICIENTS IN THE FIRST REGION
C
      DO 4000 II=1,M
      A(II)=C1(II)
      B(II)=C1(II)
 4000 CONTINUE
C
C     B=Y3*C1
C
      CALL CVMULT(M,Y3,B,MX,WRK) 
C
C
      DO 4100 I=1,M
C
      ATL=200.
      BTL=200.
      TTL=200.
      AABS=CDABS(A(I))
      BABS=CDABS(B(I))
      TABS=CDABS(A(I)+B(I))
      IF(AABS .GT. 1.0D-10) ATL=-20.*DLOG10(AABS)
      IF(BABS .GT. 1.0D-10) BTL=-20.*DLOG10(BABS)
      IF(TABS .GT. 1.0D-10) TTL=-20.*DLOG10(TABS)
C
      WRITE(9,4050) I,A(I),B(I),ATL,BTL,TTL
 4050 FORMAT(1X,I5,5X,2(1X,2E12.4),5X,F8.2,2(10X,F8.2))
 4100 CONTINUE
C
C     THIS INITIALIZES THE RANGE LOOP FOR THE PRESSURE CALCULATION
C     (BSOUT .NE. 1)
C
      RKM=RMIN
C
      IF(RANGE(1) .LT. RMIN) GO TO 4200
C
      R2KM=RANGE(1)
      R1KM=RANGE(1)
C
      JBOT=0
C
C     COMPUTE THE DEPTH FUNCTIONS AND PRESSURE IN THE FIRST REGION
C     (BSOUT .NE. 1)
C
      CALL CDFUN(M,EGVL,FNORM,H1,ZR,NDEP,CDF,MX,MD,
     1           IMAT,BEGVL,BMATL,ALPHA,FAC,RAT,BETA)
C
C
      CALL CPRES(M,RINC,RMAX,R1KM,R2KM,RKM,EGVL,A,B,
     1            NDEP,NDINC,CDF,JBOT,MX,MD,P,P0,P1,P2,
     2            TL,TL0,TL1,TL2,BSOUT,MAMP)
C
C
 4200 CONTINUE
C
C ************** LOOP FOR PRESSURE CALCULATION *********************
C **************         (BSOUT .NE. 1)         *********************
C
C
      DO 7000 J=1,NREM1
C
      IF(IMAT .GE. 1) READ(8) (BEGVR(I),I=1,M)
      IF(IMAT .GE. 1) READ(8) ((BMATR(II,I),II=1,M),I=1,M)
      READ(8) (EGVR(I),I=1,M)
      READ(8) (FNORM(I),I=1,M)
C
      H1=DEPTH(J)
      H2=DEPTH(J+1)
C
      IF(IMAT .EQ. 1) THEN
      AVG=AVGW(J+1)
      ROHB=RHOBTOP(J+1)
      CW=CB/SQRT(AVG)
      FKW=2.0*PI*FREQ/CW
      END IF
C
      READ(4) ((T(II,I),II=1,MT2),I=1,MT2)
C
 6500 CONTINUE
C 
      READ(21) ((Y1(II,I),II=1,I),I=1,M)
C
      BACKSPACE 3
      READ(3) ((Y3(II,I),II=1,M),I=1,M)
      BACKSPACE 3
C
C     GENERATE COEFFICIENTS A AND B
C
      CALL GENAB(T,Y1,Y3,C1,A,B,WRK1,CLR,M,MX,MXT2)
C
C
      JP1=J+1
      IF(IPRT .LE. 0) GO TO 6850
      IF(JP1 .EQ. NREG) GO TO 6780
      IF(MOD(JP1,IPRT) .NE. 0) GO TO 6850
 6780 CONTINUE
C
      WRITE(9,105)
      WRITE(*,3690) JP1,RANGE(JP1),DEPTH(JP1)
      WRITE(9,3690) JP1,RANGE(JP1),DEPTH(JP1)
C
      DO 6800 I=1,M
C
      ATL=200.
      BTL=200.
      TTL=200.
      AABS=CDABS(A(I))
      BABS=CDABS(B(I))
      TABS=CDABS(A(I)+B(I))
      IF(AABS .GT. 1.0D-10) ATL=-20.*DLOG10(AABS)
      IF(BABS .GT. 1.0D-10) BTL=-20.*DLOG10(BABS)
      IF(TABS .GT. 1.0D-10) TTL=-20.*DLOG10(TABS)
C
      WRITE(9,4050) I,A(I),B(I),ATL,BTL,TTL
 6800 CONTINUE
      WRITE(9,105)
C
 6850 CONTINUE
C
      IF(J .LT. NREM1) GO TO 6880
C
      BMAX=0.
      DO 6860 I=1,M
      BMAX=DMAX1(BMAX,CDABS(B(I)))
      B(I)=ZERO
 6860 CONTINUE
C
      WRITE (6,105)
      WRITE(*,6870) BMAX
      WRITE(9,6870) BMAX
 6870 FORMAT(1X,'THE MAXIMUM MAGNITUDE OF THE BACKSCATTER COEFFICIENTS',
     1 ' IN THE LAST REGION WAS',/,E12.4,25H BEFORE BEING SET TO ZERO/)
C
 6880 CONTINUE
C
      IF((RANGE(J+1) .LT. RMIN) .OR. (RANGE(J) .GE. RMAX)) GO TO 6900
      IF(RKM .GT. RANGE(J+1)) GO TO 6900
C
      R1KM=RANGE(J)
      R2KM=RANGE(J+1)
      JBOT=0
C
C     CALCULATE DEPTH FUNCTION AND PRESSURE (BSOUT .NE. 1)
C
      CALL CDFUN(M,EGVR,FNORM,H2,ZR,NDEP,CDF,MX,MD,
     1           IMAT,BEGVR,BMATR,ALPHA,FAC,RAT,BETA)
C
C
      CALL CPRES(M,RINC,RMAX,R1KM,R2KM,RKM,EGVR,A,B,
     1            NDEP,NDINC,CDF,JBOT,MX,MD,P,P0,P1,P2,
     2            TL,TL0,TL1,TL2,BSOUT,MAMP)
C
C
 6900 CONTINUE
C
      DO 6920 I=1,M
      BEGVL(I)=BEGVR(I)
      DO 6910 II=1,M
      BMATL(II,I)=BMATR(II,I)
 6910 CONTINUE
 6920 CONTINUE
C
      DO 6950 I=1,M
      EGVL(I)=EGVR(I)
 6950 CONTINUE
C
 7000 CONTINUE
 8000 CONTINUE
C
      ENDFILE 7
      ENDFILE 10
      IF(MAMP .GT. 0) ENDFILE 11
      IF(BSOUT .NE. 1) THEN
      ENDFILE 12
      ENDFILE 13
      ENDFILE 14
      ENDFILE 15
      ENDFILE 16
      ENDFILE 17
      END IF
C
      REWIND 7
      REWIND 10
      REWIND LUSVP
      IF(MAMP .GT. 0) REWIND 11
      IF(BSOUT .NE. 1) THEN
      REWIND 12
      REWIND 13
      REWIND 14
      REWIND 15
      REWIND 16
      REWIND 17
      END IF
C
C     PRINT TRANSMISSION LOSS
C
      WRITE(*,*) 'CALL PRTTL'
      WRITE(9,*) 'CALL PRTTL'
      CALL PRTTL(TITLE,FREQ,ZS,ZR,TL)
C
C     CLOSE THE SCRATCH FILES
C
      CLOSE(21)
      CLOSE(2)
      CLOSE(3)
      CLOSE(4)
      CLOSE(8)
      CLOSE(18)
      CLOSE(LUSVP)
C 
C     FIND ELAPSED TIME ON FPS 164
C
C     CALL SYS$RDTIME(T1)
C
C     FIND ELAPSED TIME IN SVS
C
C      CALL GETTIM(IH,IM,IS,IHU)
C      ENDS=(60.0*IH+IM)*60.0+IS+.01*IHU
C
C      CALL TIMER(IHU)
C      ENDS=.01*IHU
C      ELAPSE=ENDS-STARTS
C
      IENDS=time()
      ELAPSE=IENDS-ISTARTS
C
      WRITE(*,8100) ELAPSE
      WRITE(9,8100) ELAPSE
 8100 FORMAT(1X,'TOTAL ELAPSED TIME= ',F10.3,'  SECONDS')    
C
C
      STOP
      END
