      PROGRAM TRIM82
      INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL
      PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7,
     &   DSTFIL = 37, IGNFIL = 41 )
      DIMENSION Z(2,50),V(2,50)
      COMMON /MIRROR/ FSRL, BRLT(200), BPST(200)
	DATA BRLT/200*1./, BPST/200*0./
      COMMON /INFO/ R0, R, OMEGA, ATT,ICOH,IRAN,IT2,IT3,IGNRAY,IFLMR
      COMMON /PROF/ RONE, N1, Z1(50), V1(50),
     1                RTWO, N2, Z2(50), V2(50)
      COMMON /LOUDNS/ NRCD, RCD(100), XINT(400)
	COMMON /TLPLOT/ RKM(1000), XTL(1000,8), ITYP(1000), IRCD(1000)
	DATA ITYP/1000*-1./
      COMMON /PRO/ IPROP, ICUR, REARTH, R1KM, DRKM, R2KM
      COMMON /PATTRN/ SD, ITBP, DATD, SORLEV
      COMMON /PIDEF/ PI, DEGREE, TWO PI
      COMMON /TLE/ ITITLE (18)
      COMMON/ABC/ NBRS, NSRS, ALIM
      COMMON /RANST/ RSKM, IREC, IFSK, SDS
	DATA XLG/1000./, YHGT/6./, XL/1./, YL/6.4/, MAXPLOT/1000/
	DATA IPROP/1/, ICUR/1/, IFLMR/1/, REARTH/6371221.3/
C*****   VAX CONVERSION OF TRIMAIN--CF. N R L REPORT # 7827, 12/16/74.
C   E B WRIGHT  CODE 5160  NAVAL RESEARCH LAB, WASHINGTON D C  20375.
C*****   TRY TO OBSERVE VARIABLE-NAME CONVENTION: LENGTHS WITHOUT 'KM'
C	OR SOME OTHER UNIT QUALIFIER ARE ALWAYS IN METERS.
      PI = 3.141592654
      TWO PI = PI + PI
      DEGREE = 180./PI
      READ( INPFIL, 900 )  ITITLE
900   FORMAT(20A4)
901	FORMAT(16I5)
902	FORMAT(10F8.0)
903	FORMAT(16F5.0)
907	FORMAT(' RECEIVER RANGES R1KM, DRKM, R2KM = ', 3F10.3)
908	FORMAT(' RECEIVER DEPTHS(M) = ', 10F10.3)
      WRITE( PRTFIL, 900) ITITLE
      READ( INPFIL, 902 ) SD, FHZ, SORLEV, DATD, SLDB, XATT, XBOT
	IFLMR = 1
	IF(SD.GT.200. .OR. SD*FHZ.GT.15000.) IFLMR = 0
	ITBP = 0
	IF(DATD.NE.0.) ITBP = 1
	IF(DATD.GT.90.) DATD = 0.
	READ( INPFIL, 901 ) ICOH, IRAN, IT2, IT3, IGNRAY, IRD, IRF, 
	1                            IPL, NSRS, NBRS, LIMDB
	INT = IRAN
	IF(ICOH.LE.0 .AND. IGNRAY.LE.0) GO TO 15
	IF(ICOH.GT.0 .AND. IRAN.LE.0) IRAN = ICOH
	IF(IGNRAY.GT.0 .AND. IRAN.LE.0) IRAN = IGNRAY
	IF(ICOH.GT.0) INT = MIN0(ICOH,IRAN)
	IF(IGNRAY.GT.0) INT = MIN0(IGNRAY,INT)
	IF(IGNRAY.GT.0 .AND. INT.EQ.0) INT = IGNRAY
15	IF(IT2.GT.0 .AND. INT.LE.0) INT = IT2
	IF(IT3.GT.0 .AND. INT.LE.0) INT = IT3
      FSRL = 10.**(-.1*SLDB)
      OMEGA = TWO PI*FHZ
      FSQ = FHZ*FHZ
	FKHZSQ = FSQ*1.E-6
      WRITE( PRTFIL, 902) SD, SORLEV, FHZ, DATD
	SDS = SD
      IF (ICUR.NE.0) SD = SD*(1.+.5*SD/REARTH)
      ATT=0.
C*****  RELAXATION FREQUENCY AT 4.4 C IS ABOUT 70 KHZ.  ATT IN DB/KM.
C	THE FIRST FORMULA IS FOR THORP ATTENUATION.
C	THE SECOND IS CYBULSKI/GOODMAN (?).
      IF(XATT.NE.0.) ATT = 0.10936*FKHZSQ/(1.+FKHZSQ) + 
     1                         	        50.*FKHZSQ/(5000. + FKHZSQ)
C      IF(XATT.EQ.2.) ATT = 0.00044*FKHZSQ + 50.*FKHZSQ/(5000. + FKHZSQ)
	WRITE( PRTFIL, 915) ATT
915	FORMAT(' VOLUME ATTENUATION IN DB/KM = ', E10.3)
      IF(NBRS.LE.0) NBRS = 9999
      IF(NSRS.LE.0) NSRS = 9999
	IF(LIMDB.LE.0) LIMDB=200
	WRITE( PRTFIL, 901) ICOH, IRAN, IT2, IT3, NSRS, NBRS, LIMDB
	ALIM = 10.**(-.1*LIMDB)
	NRC = 7
	NRCD = 0
	READ( INPFIL, 902 ) R1KM, DRKM, R2KM, (RCD(I),I=1,7)
	IF(R1KM.EQ.DRKM) R1KM = 0.
	DELR12KM = R2KM - R1KM
	R1 = 1000.*R1KM
	DR = 1000.*DRKM
	R2 = 1000.*R2KM
	IF(RCD(7).LE.0.) GO TO 20
	NRC = 17
	READ( INPFIL, 902 ) (RCD(I),I=8,17)
	IF(RCD(17).LE.0.) GO TO 20
	NRC = 27
	READ( INPFIL, 902 ) (RCD(I),I=18,27)
20	DO 30 I=1,NRC
	IF(RCD(I).LE.0.) GO TO 30
	NRCD = I
	IF(ICUR.NE.0) RCD(I) = RCD(I)*(1.+.5*RCD(I)/REARTH)
30	CONTINUE
31      WRITE( PRTFIL, 907) R1KM, DRKM, R2KM
	IF(NRCD.GT.0 .OR. INT.EQ.0) GO TO 32
	NRCD = 1
	RCD(1) = 1.
32	WRITE( PRTFIL, 908) (RCD(I),I=1,NRCD)
	CALL INITRY
	RBKM = 50000.
	IF(XBOT.GE.0.) GO TO 36
	WRITE( PRTFIL, *) 'BOTTOM IS PERFECTLY ABSORBING.'
	DO 35 I=1,200
	BRLT(I) = 0.
35	CONTINUE
	RBLKM = 50000.
	GO TO 40
36	IF(XBOT.GT.0.) GO TO 38	
	WRITE( PRTFIL, *) 'BOTTOM IS PERFECTLY REFLECTING.'
	RBLKM = 50000.
	GO TO 40
38	WRITE( PRTFIL, *) 'BOTTOM LOSS CLASS TO BE INPUT.'
	XBOT = ICOH
	CALL BOTLOSS(XBOT)
	RBLKM = XBOT
40	RBL = 1000.*RBLKM
      CALL NEWPRO(ZM)
      CALL NEWPRO(ZM)
      CALL CONNCT
      IF(ZMAX.EQ.0.) ZMAX=ZM
	R0 = 0.
	R = 0.
	IF(R1.NE.0.) GO TO 50
	IF(IRF.GT.0) CALL RAYOUT
	GO TO 95
50	IF(RTWO.GE.R1 .AND. RBL.GE.R1) GO TO 75
	R = AMIN1(RBL,RTWO)
	CALL ADVNCE
	R0 = R
	IF(R.LT.RBL) GO TO 55
	CALL NEWBLT(RBLKM)
	RBL = 1000.*RBLKM
55	IF(R.LT.RTWO) GO TO 50
	CALL NEWPRO(ZM)
	CALL CONNCT
	GO TO 50
75	R = R1
	CALL ADVNCE
	R0 = R
	IF(INT.GT.0) CALL INTSTY(0)
	IF(IRD.GT.0) CALL ZDIST
	IF(IRF.GT.0) CALL RAYOUT
	IF(R.NE.RBL) GO TO 90
	CALL NEWBLT(RBLKM)
90	IF(R.NE.RTWO) GO TO 95
	CALL NEWPRO(ZM)
	CALL CONNCT
95	NRS = (R2 - R1)/DR + .5
	NPLOTS = 0
	IF(IPL.EQ.0 .OR. NRS.LE.MAXPLOT) GO TO 105
	WRITE( PRTFIL, 986) NPLOTS
986	FORMAT(' TOO MANY RANGE STEPS FOR NPLOTS=', I6)
	STOP
105	NEXTI = INT
	NEXTRD = IRD
	NEXTRF = IRF
	NEXTBL = (RBLKM-R1KM)/DRKM + .5
	NEXTP = (RTWO-R1)/DR + .5
	DO 1000 NR=1,NRS
	R = R + DR
	CALL ADVNCE
	R0 = R
	IF(NR.NE.NEXTI) GO TO 200
	CALL INTSTY(NEXTI)
	NEXTI = NEXTI + INT
200	IF(NR.NE.NEXTRD) GO TO 300
	CALL ZDIST
	NEXTRD = NEXTRD + IRD	
300	IF(NR.NE.NEXTRF) GO TO 400
	CALL RAYOUT
	NEXTRF = NEXTRF + IRF
400	IF(NR.NE.NEXTP) GO TO 500
	CALL NEWPRO(ZM)
	CALL CONNCT
	NEXTP = (RTWO-R1)/DR + .5
500	IF(NR.NE.NEXTBL) GO TO 800
	CALL NEWBLT(RBLKM)
	NEXTBL = (RBLKM-R1KM)/DRKM + .5
800	IF(INT.GT.0) CALL INTPRT
1000	CONTINUE	
	IF(IPL.EQ.0) GO TO 73
	NPLOTS = (ICOH+IRAN+IT2+IT3)*NRCD
	XLG = 8.
	IF(NRS.GT.500) XLG = 12.
	YHGT = 6.
	DELRKM = 100.
	IF(DELR12KM.LT.675.) DELRKM = 50.
	IF(DELR12KM.LT.375.) DELRKM = 20.
	IF(DELR12KM.LT.175.) DELRKM = 10.
	IF(DELR12KM.LT.65.)  DELRKM = 5.
	IF(DELR12KM.LT.33.)  DELRKM = 2.
	IF(DELR12KM.LT.15.)  DELRKM = 1.
73	STOP
      END

      SUBROUTINE BOTLOSS(RBKM)
      INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL
      PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7,
     &   DSTFIL = 37, IGNFIL = 41 )
	DATA IENT/0/
      COMMON /PIDEF/ PI, DEGREE, TWO PI
      COMMON /INFO/ ROLD, RANGE, OMEGA, IDUMM(7)
      COMMON /MIRROR/ FSRL,BRLT(200),BPST(200)
      DIMENSION DB(50),DG(50),DTSTR(600),BR(91,4),BP(91,4),ISW(4),
     *   RBLKM(50),ICL(50), DTH(600)
      COMMON/ABC/NBRS,NSRS,ALIM
      EQUIVALENCE (DB,DTSTR(551)),(DG,DTSTR(501)),(ICL,DTSTR(451)),
     *(RBLKM,DTSTR(401)),(ISW,DTSTR(397)),(BR,DTSTR),(ICOH,DTSTR(396)), 
     *   (IE,DTSTR(395)),(IB,DTSTR(394)),(K,DTSTR(393)),
     *   (IL,DTSTR(392)),(IH,DTSTR(391)),(IC,DTSTR(390))
      EQUIVALENCE (N,DTSTR(389)),(T,DTSTR(388)),(F,DTSTR(387)),
     *   (PHI,DTSTR(386)),(DCB,DTSTR(385)),(RKM,DTSTR(384)),
     *   (RNKM,DTSTR(383))
      EQUIVALENCE (BP,DTSTR(183)), (DTSTR,DTH1)
      COMMON/BLDATA/DTH1(100),DTH2(100),DTH3(100),DTH4(100),DTH5(100),
     1   DTH6(100)
      DATA DTH1 /7.5, 11., 13.5, 15., 17.5, 20., 25., 27.5, 2*90.,
     2           7.5, 10., 11.5, 14., 15., 17.5, 20., 25., 29.5, 90.,
     3           7.5, 12., 15., 18., 20., 22.5, 25., 27., 2*90.,
     4           7.5, 10., 14.5, 17.5, 20., 22.5, 25., 28., 2*90.,
     5           2.5, 5.5, 8.0, 10., 15., 17.5, 20., 23., 2*90.,
     1           6.5, 8.6, 9.8, 10.8, 11., 11.7, 12.7, 3*13.,
     2           7.0, 8.8, 10., 10.6, 11.1, 11.9, 12.5, 13.4, 2*14.,
     3           8.0, 11., 12.8, 14., 14.8, 15.4, 15.8, 3*16.,
     4           9.0, 11., 13.5, 14.8, 15.6, 16.2, 16.6, 3*17.,
     5           10., 13., 15., 16., 17.8, 18.4, 18.7, 3*19./
      DATA DTH2 /2.5, 10., 12.5, 15., 20., 25., 30., 35., 40., 90.,
     2           2.5, 7.5, 12.5, 20., 25., 32.5, 37.5, 42.5, 2*90.,
     3           2.5, 5.0, 7.5, 12.5, 17.5, 22.5, 27.5, 32.5, 2*90.,
     4           2.5, 7.5, 10., 11.5, 13., 15., 20., 24., 2*90.,
     5           2.5, 5.0, 8.5, 10., 12.5, 15., 20., 22.5, 27.5, 90.,
     1           3.5, 5.4, 6.0, 6.6, 7.7, 8.6, 9.4, 9.8, 2*10.,
     2           6.5, 8.0, 9.5, 11., 11.9, 13., 13.4, 3*13.5,
     3           7.0, 9.4, 11.5, 14., 14.7, 15.1, 15.4, 3*15.5,
     4           9.0, 13.5, 15.2, 16., 16.6, 17., 17.4, 3*17.5,
     5           11., 13.6, 16., 17., 18.1, 18.8, 19.6, 19.8, 2*20./
      DATA DTH3 /18., 23., 30., 40., 50., 55., 4*90.,
     2           17., 20., 25., 30., 35., 45., 55., 3*90.,
     3           13., 25., 30., 35., 42.5, 45., 50., 3*90.,
     4           9.0, 12.5, 20., 22.5, 27., 5*90.,
     5           2.5, 4.5, 7.5, 10., 14.5, 17.5, 20., 23., 2*90.,
     1           0., 2.0, 4.0, 6.0, 7.4, 5*8.,
     2           0., 1.7, 3.9, 5.6, 6.8, 8.8, 4*10.,
     3           3.0, 6.6, 8.0, 9.1, 10.4, 10.7, 4*11.,
     4           7.0, 8.0, 12., 13.1, 6*14.,
     5           9.0, 10., 12., 13.9, 16., 17.1, 17.7, 3*18./
      DATA DTH4 /19., 25., 35., 45., 51., 5*90.,
     2           18., 25., 35., 45., 55., 5*90.,
     3           14., 20., 30., 40., 50., 53., 4*90.,
     4           8.0, 11., 20., 22.5, 26., 5*90.,
     5           2.5, 4.5, 7.5, 10., 14.5, 17.5, 20., 23., 2*90.,
     1           0., 1.5, 3.5, 5.1, 6*6.,
     2           0., 2.8, 5.6, 7.6, 6*9.,
     3           3.0, 5.0, 7.5, 9.4, 10.8, 5*11.,
     4           5.0, 6.1, 12.2, 13.2, 6*14.,
     5           9.0, 10., 12., 13.9, 16., 17.1, 17.7, 3*18./
      DATA DTH5 /18., 30., 40., 50., 55., 5*90.,
     2           17., 20., 30., 40., 50., 55., 4*90.,
     3           13., 20., 35., 45., 52.5, 5*90.,
     4           7.5, 11., 20., 25., 30., 34., 4*90.,
     5           2.5, 5.0, 7.5, 15., 17.5, 20., 22., 3*90.,
     1           0., 2.6, 4.4, 5.6, 6*6.,
     2           0., 1.3, 4.6, 7.0, 8.5, 5*9.,
     3           3.0, 5.3, 8.7, 10.3, 6*11.,
     4           3.0, 4.0, 10.4, 12.7, 13.7, 5*14.,
     5           8.0, 8.7, 10., 14.1, 15., 15.8, 4*16./
      DATA DTH6 /19., 35., 40., 52., 6*90.,
     2           19., 25., 35., 45., 55., 5*90.,
     3           11., 20., 25., 35., 45., 56., 4*90.,
     4           11., 20., 30., 40., 55., 5*90.,
     5            5., 7.5, 16.5, 17.5, 19., 20., 21., 22.5, 2*90.,
     1            0., 3.7, 4.5, 7*6.,
     2            0., 2.3, 4.6, 6.5, 6*8.,
     3            0., 3.0, 4.4, 6.7, 8.5, 5*10.,
     4            2., 5.2, 7.7, 9.8, 6*12.,
     5            4., 5.2, 12., 12.8, 13.4, 13.7, 13.9, 3*14./
      IF(IENT.EQ.1) GO TO 2
      DO 1 J=1,600
    1 DTH(J) = DTSTR(J)
      IENT = 1
      GO TO 3
    2 DO 4 J=1,600
    4 DTSTR(J) = DTH(J)
      IB=0
    3 FR=OMEGA/TWOPI
      IF(FR.LT.300.) GO TO 25
      IF(FR.GT.1500.) GO TO 10
      IF(FR.LT.750.) GO TO 5
      IQ=301
      GO TO 20
5     IQ=401
      GO TO 20
10    IF(FR.LT.2750.) GO TO 15
      IQ=101
      IF(FR.GT.6000.) IQ=1
      GO TO 20
15    IQ=201
20    DO 24 I=501,600
      DTSTR(I)=DTH1(IQ)
24    IQ=IQ+1
25    DO 30 I=1,500
30    DTSTR(I)=0.
      ICOH=RBKM
40    IB=IB+1
      READ( INPFIL, 900 ) RBLKM(IB), ICL(IB)
      IF(RBLKM(IB)) 41,41,42
41    RBLKM(IB)= 50000.
42    IF(ICL(IB)) 43,43,44
43    ICL(IB)=0
      GO TO 60
44    IF(ICL(IB)-9) 46,46,45
45    ICL(IB)=ICL(IB)/10
46    IF(ICL(IB)-5) 60,60,47
47    K=ICL(IB)-5
      IF(ISW(K)) 48,48,60
48    ISW(K)=1
      IL=1
49    IH=MIN0(91,IL+15)
 
      READ( INPFIL, 901 ) (BR(I,K), I=IL,IH)
      IL=IH + 1
      IF(BR(IH,K)) 51,51,49
50    IF(BR(IH,K)) 51,51,53
51    IH=IH-1
      IF(IH) 52,52,50
52    IH=1
53    DO 54 I=IH,90
54    BR(I+1,K)=BR(IH,K)
      IL=1
      IF(ICOH) 60,60,55
55    IH=MIN0(91,IL+15)
      READ( INPFIL, 901 ) (BP(I,K),I=IL,IH)
      IL=IH + 1
      IF(BP(IH,K)) 55,60,55
60    IF(RBLKM(IB).LT.40000.) GO TO 40
      DO 65 K=6,9
      IF(ISW(K-5)) 65,65,61
61    WRITE( PRTFIL, 902) K
      WRITE( PRTFIL, 903) (I,BR(I+1,K-5),BP(I+1,K-5),I=IE,90)
65    CONTINUE
      WRITE( PRTFIL, 900)
      IF(IB.LT.2) GO TO 71
      DO 70 I=2,IB
      RNKM =RBLKM(I-1)
      WRITE( PRTFIL, 904) RKM,RNKM,ICL(I-1)
70    RKM =RNKM
71    WRITE( PRTFIL, 905) RKM,ICL(IB)
      IE=1
75    IC=ICL(IE)
      IF(IC) 76,76,80
76    DO 77 I=1,200
      BRLT(I)=1.
77    BPST(I)=0.
      GO TO 100
80    DO 99 I=1,200
      T=DEGREE*ATAN(.01*(I-1))
      K=IC-5
      IF(K) 81,81,95
81    IL=10*IC
      IF(DG(IL-9).GE.T) GO TO 93
82    IF(DG(IL-8).GE.T) GO TO 92
      IL=IL+1
      GO TO 82
92    DCB=DB(IL-9)+(T-DG(IL-9))*(DB(IL-8)-DB(IL-9))/(DG(IL-8)-DG(IL-9))
      GO TO 94
93    DCB=DB(IL-9)
94    PHI=0.
      GO TO 98
95    N=T
      F=T-N
      DCB=BR(N+1,K)+F*(BR(N+2,K)-BR(N+1,K))
      PHI=0.
      IF(ICOH.NE.0) PHI=BP(N+1,K)+F*(BP(N+2,K)-BP(N+1,K))
98    BPST(I)=PHI
99    BRLT(I) = 10.**(-.1*DCB)
100   RBKM=RBLKM(IE)
      RETURN
      ENTRY NEWBLT(RBKM)
      IE=IE+1
      IF(IC-ICL(IE)) 75,100,75
900   FORMAT(F8.4,I2)
901   FORMAT(16F5.2)
902   FORMAT(38H0USER SUPPLIED BOTTOM LOSS TABLE CLASS  ,I5/,
     1   5(7X,'GA     DB      PHI '))
  903 FORMAT(5(I9,F9.3,F8.3))
904   FORMAT(5H FROM,F10.2,3H TO,F10.2,16H KM BOTTOM CLASS,I5)
905   FORMAT(5H FROM,F10.2,3H TO,3X,23HEND OF RUN BOTTOM CLASS,I5)
      END
 
      SUBROUTINE NEWPRO(ZMAX)
      INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL
      PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7,
     &   DSTFIL = 37, IGNFIL = 41 )
      COMMON /INFO/ ROLD, RANGE, OMEGA,ATT,ICOH,IRAN,IT2,IT3,IGNRAY,IFLMR
      COMMON /PROF/ R1,N1,Z1(50),V1(50),R2,N2,Z2(50),V2(50)
      COMMON/INPUT/ RBOTKM(250),ZBOT(250),ZIL(2,50),VIL(2,50),
     * NN(2), ICON(2,50)
      COMMON /TLE/ ITITLE (18)
      COMMON /PRO/ IPROP, ICUR, REARTH, R1KM, DRKM, R2KM
      DIMENSION TITLE(10),PM(3),LIN (101)
      COMMON/ABC/NBRS,NSRS,ALIM
      COMMON /RANST/ RSKM,IDUM,IDUMM,SDS
      DATA NBPT/0/
      DATA IFEOF /0/
      DATA LIN/101*1H  /
      DATA IENT/0/
      DATA IBLANK /'    '/ , IVEE/'V   '/ , ISTAR/'*   '/
      DATA RAKM/0./, RBKM/0./
  952 IF (IENT.EQ.1) GO TO 950
      NBPT=0
      IFEOF=0
      IENT=1
      N2=0
      R2 = 0.
      N1=0
      R1 = 0.
      DO 951 J=1,101
  951 LIN(J)=IBLANK
  950 IF(NBPT.GT.0) GO TO 20
C        FIRST ENTRY.  READ IN BOTTOM TRACK
      IL=1
2     IH=IL+4
      READ( INPFIL, 900 ) (RBOTKM(I), ZBOT(I), I=IL,IH)
      IF(RBOTKM(IH).LE.0.) GO TO 3
      IL=IL+5
      IF(IH.LT.250) GO TO 2
      IH = 251
3     IH=IH-1
      IF(RBOTKM(IH).LE.0.) GO TO 3
      NBPT=IH
      IF (ICUR.NE.0) GO TO 1
      DO 5 I=1,NBPT
    5 ZBOT(I)=ZBOT(I)*(1.+.5*ZBOT(I)/REARTH)
    1 CONTINUE
  611 FORMAT (20A4,/,(10F8.3))
  911 FORMAT (10F8.3)
      ZMAX=0.
      DO 4 I=1,NBPT
	XI = (RBOTKM(I)-R1KM)/DRKM
        II = XI + .5
	IF(XI.GE.0.) RBOTKM(I) = R1KM + II*DRKM
      IF(ZBOT(I).GT.ZMAX) ZMAX=ZBOT(I)
4     CONTINUE
      WRITE( PRTFIL, 901)(I,RBOTKM(I),ZBOT(I),I=1,NBPT)
      IBPT=0
      DO 55 I=1,2
      DO 55 J=1,50
      ZIL(I,J) = 0.
      VIL(I,J) = 0.
   55 ICON(I,J)=0
      NN(1)=0
      NIB=0
      DO 56 J=1,50
      Z1(J) = 0.
      V1(J) = 0.
      Z2(J) = 0.
   56 V2(J) = 0.
      ASSIGN 6 TO IRET
      GO TO 70
C        AT 70 IS THE READIN AND CONNECT ROUTINE.  RETURN IS TO IRET
6     ASSIGN 20 TO IRET
      GO TO 70
20    IBPT=IBPT+1
      IF(IBPT.LE.NBPT) GO TO 21
      RBTKM=.001*R2+100.
      GO TO 22
21    RBTKM=RBOTKM(IBPT)
      ZBT=ZBOT(IBPT)
22    IF(RBTKM.LE.RBKM) GO TO 25
C        GET NEW INPUT PROFILE
      ASSIGN 23 TO IRET
      GO TO 70
   23 IF (IBPT.GT.NBPT .OR. RBOTKM(IBPT-1).EQ.RAKM) GO TO 25
      IBPT=IBPT-1
      ZBT=ZBT+(RAKM-RBTKM)*(ZBT-ZBOT(IBPT))/(RBTKM-RBOTKM(IBPT))
C        THIS IS THE BOTTOM DEPTH AT PROFILE RANGE
      RBTKM=RAKM
25    R1=R2
      N1=N2
C        MOVE PROFILE 2 TO PROFILE 1
      IF(N2.LT.0 .OR. N2.GT.50) GO TO 27
      DO 26 I=1,N2
      Z1(I)=Z2(I)
 26   V1(I)=V2(I)
   27 R2=1000.*RBTKM
      FA=(RBKM-RBTKM)/(RBKM-RAKM)
      FB=1.-FA
      N2=1
   29 IF (ICON(1,N2).GT.NN(1)) ICON (1,N2)=NN(1)
      IF (ICON(2,N2).GT.NN(2)) ICON (2,N2)=NN(2)
      Z=FA*ZIL(1,ICON(1,N2))+FB*ZIL(2,ICON(2,N2))
      V=FA*VIL(1,ICON(1,N2))+FB*VIL(2,ICON(2,N2))
      IF (N2.LE.1 .OR. Z.GE.Z2(N2-1)) GO TO 126
  125 Z=Z2(N2-1)
      V=V2(N2-1)
  126 IF(Z.GT.ZBT-1.) GO TO 35
      Z2(N2)=Z
      V2(N2)=V
      N2=N2+1
      IF(N2.LE.NIB) GO TO 29
C        EXTRAPOLATION TO BOTTOM
      Z2(N2)=ZBT
      M=N2-1
31    M=M-1
      Z=Z2(M)
      IF(ABS(Z-Z2(N2-1)).LT.Z*.00001) GO TO 31
      V=V2(M)
      V2(N2)=V2(M)  +(ZBT-Z)*(V2(N2-1)-V)/(Z2(N2-1)-Z)
      GO TO 40
35    Z2(N2)=ZBT
      V2(N2)=V2(N2-1)+(ZBT-Z2(N2-1))*(V-V2(N2-1))/(Z-Z2(N2-1))
C        REMOVE DUPLICATE POINTS
40    M=2
      IDC=0
   42 IF (Z2(M).GT.Z2(M-1)+1.) GO TO 45
      IDC=1
      N2=N2-1
      DO 44 I=M,N2
      Z2(I)=Z2(I+1)
44    V2(I)=V2(I+1)
45    M=M+1
      IF(M.LE.N2) GO TO 42
      IF (IDC.GT.0) GO TO 40
      IF (Z2(N2+1).EQ.0.) GO TO 46
      N5=N2+1
      DO 47 I=N5,50
      Z2(I) = 0.
   47 V2(I) = 0.
   46 	IF (IPROP.EQ.0.OR.IPROP.EQ.1) GO TO 54
	WRITE( PRTFIL, 902) R2
            WRITE( PRTFIL, 914)
      DO 7 J=1,N2
      NV=V2(J)-1455.5
      IF (NV.LT.4.OR.NV.GT.101) GO TO 48
      LIN(NV)= ISTAR
      WRITE( PRTFIL, 915) Z2(J),V2(J),LIN
      LIN(NV)= IBLANK
      GO TO 7
   48 WRITE( PRTFIL, 915) Z2(J),V2(J)
    7 CONTINUE
C54      IF (IPFL.EQ.2) CALL PRPLIT (RBKM,N3,INPF)
54      IF (RSKM.GE.RBOTKM(2).OR.IFEOF.EQ.1) RETURN
      WRITE( PRTFIL, 275)
      STOP
  275 FORMAT (5X,'PROGRAM ABORTED, TWO SOUND SPEED PROFILES INPUT BEFORE
     * SECOND BOTTOM POINT')
70    IF(NN(2).EQ.0) GO TO 75
      NB=NN(2)
      DO 71 I=1,NB
      ZIL(1,I)=ZIL(2,I)
71    VIL(1,I)=VIL(2,I)
      RAKM=RBKM
75    NN(1)=NN(2)
      IF(IFEOF.EQ.0) GO TO 76
      RBKM=RBKM+1.E6
      GO TO 82
   76 READ( INPFIL, 903 ) NCUR, RBKM, TITLE
	IF(NCUR.GE.0) GO TO 78
77    IFEOF=1
      RBKM=RAKM+1.E6
      GO TO 82
78    IL=1
      RSKM=RBKM
79    IH=IL+4
      READ( INPFIL, 900 ) (ZIL(2,I), VIL(2,I), I=IL,IH)
      IF(VIL(2,IH).LE.0.) GO TO 81
      IL=IL+5
      IF(IH.LT.50) GO TO 79
      IH = 51
81    IH=IH-1
      IF(VIL(2,IH).LE.0.) GO TO 81
      NN(2)=IH
      N3=IH
      INPF=0
      IF (IPROP.EQ.0) GO TO 74
      WRITE( PRTFIL, 904) RBKM,TITLE
      WRITE( PRTFIL, 914)
   74 DO 28 I=1,IH
      IF (VIL(2,I).GE.1400. .AND. VIL(2,I).LE.10000.)  GO TO 57
   58 WRITE( PRTFIL, 920) RBKM,ZIL(2,I),VIL(2,I)
  920 FORMAT(1H1,5X,'INPUT ERROR IN PROFILE AT RANGE',F10.4,5X,'VALUE AT
     * DEPTH',F10.3,'=',F10.3,' M/SEC')
      STOP
   57 IF(ZIL(2,I).NE.SDS) GO TO 257
      ZIL(2,I)=ZIL(2,I)+0.1
      WRITE( PRTFIL, 258) SDS,RBKM, ZIL(2,I)
  258 FORMAT (5X,'PROGRAM WARNING, INPUT PROFILE DEPTH',F10.3,1X,'=SOURC
     *E DEPTH, AT RANGE',F10.3,' KM',5X,'NEW PROFILE DEPTH=',F10.3,/)
  257 IF (I.EQ.1) GO TO 157
C*****   ISOVELOCITY LAYERS ELIMINATED 9/80 BECAUSE THEY CAN LEAD TO
C   ROTATION ANGLE OF 90 DEGREES FOR A TRIANGLE AND DIVISION BY THE
C   COSINE OF ROTATION ANGLE=CT IS USED IN RAY INTENSITY IN ADVNCE.
C*****   THE ASC PRECISION FOR REAL*4 NUMBERS BETWEEN 1024. AND 2048. IS
C   2048/16**6 = .000122.  USE .0004 AS A DISTINGUISHABLE DIFFERENCE.
      IF(VIL(2,I).EQ.VIL(2,I-1)) VIL(2,I) = VIL(2,I) + .0004
      IF (ZIL(2,I).LT.ZIL (2,I-1)) GO TO 58
  157 IF (IPROP.EQ.0) GO TO 28
      NV=VIL(2,I)-1455.5
      IF (NV.LT.4.OR.NV.GT.101) GO TO 38
      LIN(NV)= IVEE
      WRITE( PRTFIL, 915) ZIL(2,I),VIL(2,I),LIN
      LIN (NV)= IBLANK
      GO TO 28
   38 WRITE( PRTFIL, 915) ZIL(2,I),VIL(2,I)
   28 CONTINUE
      IF (NCUR.NE.0) GO TO 89
	WRITE( PRTFIL, *) 'CURVED-EARTH CORRECTION APPLIED TO PROFILE.'
      DO 80 I=1,IH
      VIL(2,I)=VIL(2,I)*(1.+ZIL(2,I)/REARTH)
   80 ZIL(2,I)=ZIL(2,I)*(1.+.5*ZIL(2,I)/REARTH)
C89      IF (IPFL.EQ.2) CALL PRPLIT (RBKM,N3,INPF)
89      IF(NN(1).EQ.0)              GO TO IRET,(6,20,23)
82    ICON(1,1)=1
      ICON(2,1)=1
      NA=2
      NB=2
      NIB=2
   83 IF (ZIL(2,NIB).EQ.ZIL(1,NIB)) GO TO 92
      IF(NB.LE.NN(2)) GO TO 84
      IP=2
      GO TO 88
84    IF(NA.LE.NN(1)) GO TO 85
      IP=1
      GO TO 88
85    DVMA=VIL(1,NA)-VIL(1,NA-1)
      DVLA=DVMA
      DVHA=DVMA
      IF(NA.GT.2) DVLA=VIL(1,NA-1)-VIL(1,NA-2)
      IF(NA.LT.NN(1)) DVHA=VIL(1,NA+1)-VIL(1,NA)
      DVMB=VIL(2,NB)-VIL(2,NB-1)
      DVLB=DVMB
      DVHB=DVMB
      IF(NB.GT.2) DVLB=VIL(2,NB-1)-VIL(2,NB-2)
      IF(NB.LT.NN(2)) DVHB=VIL(2,NB+1)-VIL(2,NB)
      PM(1)=ABS(ZIL(1,NA-1)-ZIL(2,NB))
      IF(DVLA*DVMB.GE.0..AND.DVMA*DVHB.GE.0.) PM(1)=PM(1)-250.
      PM(2)=ABS(ZIL(1,NA)-ZIL(2,NB-1))
      IF(DVMA*DVLB.GE.0..AND.DVHA*DVMB.GE.0.) PM(2)=PM(2)-250.
      PM(3)=ABS(ZIL(1,NA)-ZIL(2,NB))
      IF(DVMA*DVMB.GE.0..AND.DVHA*DVHB.GE.0.) PM(3)=PM(3)-250.
      BP=PM(1)
      IP=1
      DO 87 I=2,3
      IF(PM(I).GE.BP) GO TO 87
      IP=I
      BP=PM(I)
87    CONTINUE
88    GO TO (90,91,92),IP
90    ICON(1,NIB)=NA-1
      ICON(2,NIB)=NB
      NB=NB+1
      GO TO 93
91    ICON(1,NIB)=NA
      ICON(2,NIB)=NB-1
      NA=NA+1
      GO TO 93
92    ICON(1,NIB)=NA
      ICON(2,NIB)=NB
      NA=NA+1
      NB=NB+1
93    IF(NA.GT.NN(1).AND.NB.GT.NN(2)) GO TO IRET,(6,20,23)
      NIB=NIB+1
      IF (NIB.LE.50) GO TO 83
      WRITE( PRTFIL, 905)
      STOP
900   FORMAT(10F8.4)
  901 FORMAT('0LISTING OF BOTTOM TRACK'//5(9X, 'R(KM)  Z(M)')
     1                                  //(5(I6, F8.2, F6.0)))
902   FORMAT(30H0INTERPOLATED PROFILE AT RANGE,F10.0,2H M//)
903   FORMAT(I8,F8.0,10A4)
904   FORMAT('1INPUT PROFILE AT RANGE' ,F11.3,3H KM,5X,10A4//)
  914 FORMAT (3X,'DEPTH (M)',2X,'C(M/SEC)','1460',6X,'1470',6X,'1480',
     *6X,'1490',6X,'1500',6X,'1510',6X,'1520',6X,'1530',6X,'1540',6X,
     *'1550',6X,'1560',//)
  915 FORMAT (1X,F11.4,F10.4,101A1)
  905 FORMAT ( '0TOO MANY POINTS INTERPOLATED' )
      END
 
      SUBROUTINE CHANNL(ZR,T,ZTO,ZTU)
      INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL
      PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7,
     &   DSTFIL = 37, IGNFIL = 41 )
      COMMON /VELPRF/ N,Z(100),C(100)
      REAL*8 ARG,V,DR,CM,TG1,TG2,F,RCYCLE
      CALL VELCAL
20    DO 25 I=2,N
      IF(ZR.GT.Z(I)) GO TO 25
      IZ=I
      F=(ZR-Z(I-1))/(Z(I)-Z(I-1))
      CM=DSQRT(1.D0+T*T)*(F*C(I)+(1.D0-F)*C(I-1))
      GO TO 30
25    CONTINUE
      ZTU=-1.
      ZTO=-1.
      RETURN
30    IN=IZ-1
31    IF(C(IN).GT.CM) GO TO 35
      IN=IN-1
      IF(IN.NE.0) GO TO 31
      ZTO=0.
      GO TO 40
35    F=(CM-C(IN+1))/(C(IN)-C(IN+1))
      ZTO=F*Z(IN)+(1.-F)*Z(IN+1)
40    IN=IZ
41    IF(C(IN).GT.CM) GO TO 45
      IN=IN+1
      IF(IN.LE.N) GO TO 41
      ZTU=Z(N)
      RETURN
45    F=(CM-C(IN-1))/(C(IN)-C(IN-1))
      ZTU=F*Z(IN)+(1.-F)*Z(IN-1)
      RETURN
      ENTRY RCALC(ZR,T,ZTO,ZTU)
      IF(ZTO.GT.-.5) GO TO 50
      RCYCLE=1. D30
      RETURN
50    IN=1
      RCYCLE=0.
      TG2=0.D0
      DO 60 I=1,N
      IF(Z(I).GT.ZTO.AND.Z(I).LT.ZTU) GO TO 55
      IF(IN.EQ.1) GO TO 60
      RCYCLE = RCYCLE + 2.D0*(ZTU-Z(I-1))/TG1
      RETURN
   55 ARG = CM/C(I)
      ARG = ARG*ARG-1.D0
      TG1 = 0.D0
      IF(ARG.GT.0.D0) TG1 = DSQRT(ARG)
      IF(I.EQ.1) GO TO 59
      IF(IN.EQ.1) GO TO 56
      DR = 2.D0*(Z(I)-Z(I-1))/(TG1+TG2)
      GO TO 57
   56 DR = 1. D30
      IF(TG1.NE.0.D0) DR = 2.D0*(Z(I)-ZTO)/TG1
57    RCYCLE=RCYCLE+DR
59    TG2=TG1
      IN=2
60    CONTINUE
      RETURN
      ENTRY WDENS(ZR,T)
      IF(RCYCLE.GT.1.D10) GO TO 71
      IN=IZ
65    IF(Z(IN).GT.ZR) GO TO 70
      IF(Z(IN+1).GT.ZR) GO TO 75
      IN=IN+1
      IF(IN.GE.N) GO TO 71
      GO TO 65
70    IF(Z(IN-1).LE.ZR) GO TO 74
      IN=IN-1
      IF(IN.GE.1) GO TO 65
71    T=0.
      RETURN
74    IN=IN-1
75    V=C(IN)+(C(IN+1)-C(IN))*(ZR-Z(IN))/(Z(IN+1)-Z(IN))
      ARG = CM/V
      ARG = ARG*ARG - 1.D0
      T = 1. E+10
      IF(ARG.GT.0.D0) T = 1.D0/(RCYCLE*DSQRT(ARG))
      RETURN
      END
 
      SUBROUTINE ZDIST 
      INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL
      PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7,
     &   DSTFIL = 37, IGNFIL = 41 )
      COMMON/RAYS/NRAY,TH0(1000),TGAM(1000),ZZ(1000),SS(1000),DSUM(1000)
     1   ,TIME(1000),PHASE(1000),NBR(1000),NTU(1000),NSR(1000),NTO(1000)
      COMMON /START/ CT0(1000)
      COMMON /INFO/ ROLD, RANGE, OMEGA,ATT,ICOH,IRAN,IT2,IT3,IGNRAY,IFLMR
      COMMON /SCRACH/ DDD(182)
      COMMON /TLE/ ITITLE (18)
      DIMENSION LINE(51)
      COMMON /PIDEF/ PI, DEGREE, TWO PI
      EQUIVALENCE (LINE,DDD)
      COMMON /ABC/  NBRS, NSRS, ALIM
      DIMENSION IDEPTH(6)
      DATA IBLANK/' '/, IPLUS/'+'/, IMINUS/'-'/, ISTAR/'*'/
      DATA IBEE/'B'/, IESS/'S'/, IDOT/'.'/
    3 WRITE( DSTFIL, 2) ITITLE
    1 CALL CHANNL(1.,1000.,ZS,ZB)
      DO 5 I=1,51
    5 LINE(I) = IBLANK
      DO 51 I=1,51,10
      LINE(I) = IDOT
   51 CONTINUE
      RKM = .001*RANGE
      SCALE = .2*ZB
      DO 75 I=1,6
      IDEPTH(I) = (I-1)*SCALE + .5
   75 CONTINUE
      WRITE( DSTFIL, 900) RKM, IDEPTH, LINE
      DO 8 I=1,51
    8 LINE(I) = IBLANK
      DO 40 I=1,NRAY
   35 IF(SS(I).GT.ALIM .AND.NBR(I).LE.NBRS .AND.NSR(I).LE.NSRS) GO TO 12
   15 WRITE( DSTFIL, 901) I, TH0(I)
900   FORMAT('0RAY DEPTH DISTRIBUTION AT', F9.2,' KM.', 51X, 'DEPTH (M)'
     1   , /, 60X, 5I10, I9, /,
     2   '  NRAY THETA0 NBR NTU NSR NTO DEPTH  THETA  LEVEL    TIME',
     3   '   CMEAN   ', 51A1)
  901 FORMAT(1X, I5, F7.3, 10X, 'TERMINATED.')
      GO TO 40
   12 T=TIME(I)
      CMEAN = RANGE/T
      Z=ZZ(I)
      LINE(1)= IESS
      LINE(51)= IBEE
      CALL CHANNL(Z,TGAM(I),ZTO,ZTU)
      IZP = 50.*ZTU/ZB + 1.5
      IF(IZP.LT.1 .OR. IZP.GT.51) GO TO 45
      LINE(IZP)= IPLUS
   45 IZM = 50.*ZTO/ZB + 1.5
      IF(IZM.LT.1 .OR. IZM.GT.51) GO TO 55
      LINE(IZM)= IMINUS
   55 IZR = 50.*Z/ZB + 1.5
      IF(IZR.LT.1 .OR. IZR.GT.51) GO TO 60
      LINE(IZR)= ISTAR
   60 TGR = TGAM(I)
      STR = TGR/SQRT(1.+TGR*TGR)
      TH = DEGREE*ATAN(TGR)
      RTL = SS(I)
      IF(DSUM(I).NE.0) RTL = RTL*CT0(I)/(RANGE*STR*DSUM(I))
      RTL = 10.*ALOG10(ABS(RTL))
   25 WRITE( DSTFIL, 902) I, TH0(I), NBR(I), NTU(I), NSR(I), NTO(I),
     1                         Z, TH, RTL, T, CMEAN, LINE
  902 FORMAT(1X, I5, F7.3, 4I4, F6.0, 2F7.1, F9.3, F7.1, 3X, 51A1)
   30 IF (IZP.LE.51.AND.IZP.GE.1) LINE(IZP)= IBLANK
      IF (IZM.LE.51.AND.IZM.GE.1) LINE(IZM)= IBLANK
      IF (IZR.LE.51.AND.IZR.GE.1) LINE(IZR)= IBLANK
40    CONTINUE
    2 FORMAT ('1',20A4,/)
      RETURN
      END
 
      SUBROUTINE INITRY
      INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL
      PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7,
     &   DSTFIL = 37, IGNFIL = 41 )
      DIMENSION VBP(91)
      COMMON /PATTRN/ SD,ITBP,DATD,SORLEV
      COMMON /PIDEF/ PI, DEGREE, TWO PI
      COMMON/RAYS/NRAY,TH0(1000),TGAM(1000),ZZ(1000),SS(1000),DSUM(1000)
     1   ,TIME(1000),PHASE(1000),NBR(1000),NTU(1000),NSR(1000),NTO(1000)
      COMMON/START/CT0(1000)
      NRAY=0
      OG=1000.
      NFAN = 0
      WRITE( PRTFIL, 905)
  905 FORMAT('0 THE SPECIFIED ANGULAR APERTURES FOLLOW. WITHIN EACH FAN'
     1   ,' THE RAYS TRACED ARE CENTERED IN THE ANGULAR INCREMENTS.',/,
     2 '  FAN   DOWNLIMIT  UPLIMIT INCREMENT NRAYS    LEVEL     PHASE')
   10 NFAN = NFAN + 1
      READ( INPFIL, 900 ) GAMDD, DGAMD, GAMUD, SL, PH, XCONT
      NRAYS = (GAMUD - GAMDD)/DGAMD + .1
      WRITE( PRTFIL, 906) NFAN, GAMDD, GAMUD, DGAMD, NRAYS, SL, PH
  906 FORMAT(1H0, I4, 3F10.4, I5, 2F10.2)
      SL = 10.**(.1*(SL+SORLEV))
C     IFF=1
      IF(ABS(OG-GAMDD).LT..001) GO TO 20
      IF(NRAY.EQ.0) GO TO 15
C*****   THE RAY WITH ZERO INTENSITY GOES BETWEEN NON-CONTIGUOUS FANS.
      NRAY=NRAY+1
      ZZ(NRAY) = SD
      TIME(NRAY) = 0.
      NBR(NRAY) = 0
      NTU(NRAY) = 0
      NSR(NRAY) = 0
      NTO(NRAY) = 0
      PHASE(NRAY) = PH/DEGREE
      TH0(NRAY) = OG
      S1 = SIN(OG/DEGREE)
      TGAM(NRAY) = S1/SQRT(1. - S1*S1)
      SS(NRAY)=0.
      CT0(NRAY) = 0.
15    OG=GAMDD
      SOG=SIN(OG/DEGREE)
C     IFF=0
20    G=OG+DGAMD
      S=SIN(G/DEGREE)
      NRAY=NRAY+1
      IF(NRAY.GT.1000) GO TO 101
      TH0(NRAY) = G - .5*DGAMD
      S1 = SIN(TH0(NRAY)/DEGREE)
      CT0(NRAY) = SQRT(1. - S1*S1)
      TGAM(NRAY) = S1/CT0(NRAY)
      ZZ(NRAY)=SD
      TIME(NRAY)= 0.
      DSUM(NRAY) = 0.
      NBR(NRAY) = 0
      NTU(NRAY) = 0
      NSR(NRAY) = 0
      NTO(NRAY) = 0
      PHASE(NRAY)=PH/DEGREE
      SS(NRAY) = SL*ABS(S-SOG)
      CT0(NRAY) = CT0(NRAY)/SS(NRAY)
C
C     IFF=1
      OG=G
      SOG=S
      IF ( (G + .5*DGAMD).LT. GAMUD ) GO TO 20
      IF(XCONT.NE.0.) GO TO 10
23    IF(ITBP.EQ.0) GO TO 35
      IL=1
25    IH=IL+15
      IF(IH.GT.91) IH = 91
      READ( INPFIL, 901 ) (VBP(I),I=IL,IH)
      WRITE( PRTFIL, 907) (VBP(I),I=IL,IH)
907   FORMAT('0VBP ',16F5.1)
      IF(VBP(IH).EQ.0.) GO TO 27
      IL=IL+16
      GO TO 25
27    IH=IH-1
      IF(IH.EQ.1) GO TO 28
      IF(VBP(IH).EQ.0.) GO TO 27
28    DO 29 I=IH,90
29    VBP(I+1)=VBP(IH)
      DO 31 I=1,NRAY
      D=ABS(DEGREE*ATAN(TGAM(I))+DATD)
      N=D
      D=D-N
      D=(1.-D)*VBP(N+1)+D*VBP(N+2)
31    SS(I) = SS(I)*10.**(-.1*D)
      WRITE( PRTFIL, 902)
      WRITE( PRTFIL, 903) (I, TH0(I), SS(I), PHASE(I), I=1,NRAY)
   35 RETURN
  101 NRAY=1000
      WRITE( PRTFIL, 908)
100   IF(XCONT.EQ.0.) GO TO 23
      READ( INPFIL, 900 ) XCONT
      GO TO 100
900   FORMAT (10F8.0)
901   FORMAT(16F5.1)
902   FORMAT('0  INITIAL RAY ANGLE, SIGNAL LEVEL, AND PHASE'//)
  903 FORMAT(5(I6, F8.3, F7.4, F5.2))
  908 FORMAT ('0TOO MANY RAYS INPUT, FIRST 1000 RAYS WILL BE USED')
      END
 
      SUBROUTINE INTSTY(INT)
      INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL
      PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7,
     &   DSTFIL = 37, IGNFIL = 41 )
      COMMON /LOUDNS/ NRCD,RCD(100),TYPE1C(100),TYPE1R(100),
     *   TYPE2(100),TYPE3(100)
      COMMON /TLE/ ITITLE (18)
      COMMON/RAYS/NRAY,TH0(1000),TGAM(1000),ZZ(1000),SS(1000),DSUM(1000)
     1   ,TIME(1000),PHASE(1000),NBR(1000),NTU(1000),NSR(1000),NTO(1000)
      COMMON /INFO/ROLD,RANGE,OMEGA,ATT,ICOH,IRAN,IT2,IT3,IGNRAY,IFLMR  
      COMMON /PIDEF/ PI, DEGREE, TWO PI
      COMMON /VELPRF/NN,Z(100),C(100)
      COMMON /SCRACH/ DDD(182)
      DIMENSION ID(1000), IHP(100), SINT(100,4), SSL(100), QUAD(100),
     1                QUS(100), TYPS(100)
      EQUIVALENCE (QUAD,DDD)
      COMMON/ABC/NBRS,NSRS,ALIM
      DATA IBLANK/'    '/, ISTAR/'*   '/, ISTARS/'**  '/
      DATA IHP /100*0/, SINT /400*0./, ID /1000*0/, IENT/0/
      SEC(T)=FVA*SQRT(1.+T*T)
    3 FVA=1.
      IF(ATT.NE.0.) FVA = 10.**(-.0001*ATT*RANGE)
      CALL VELCAL
      WAVEL=TWOPI*C(1)/OMEGA
      AK=TWOPI/WAVEL
      FLM=1.
      DO 5 I=1,NRCD
      TYPE1C(I) = 0.
      TYPE1R(I) = 0.
      TYPE2(I) = 0.
      TYPE3(I) = 0.
    5 QUAD(I) = 0.
      IF(IRAN.EQ.0.AND.IT2.EQ.0) GO TO 60
      CALL CHANNL(1.,1000.,ZS,ZB)
      IF(IRAN.EQ.0) GO TO 35
      DO 20 I=2,NRAY
      NDTUO = IABS(NTO(I)-NTO(I-1)) + IABS(NTU(I)-NTU(I-1))
      NDSBR = IABS(NSR(I)-NSR(I-1)) + IABS(NBR(I)-NBR(I-1))
      ID(I) = NDSBR + NDTUO
      IF(ID(I).GT.3) ID(I)=3
      W=AMIN1(SS(I),SS(I-1))
      IF(W.LE.ALIM .OR. NBR(I).GT.NBRS .OR. NSR(I).GT.NSRS) ID(I)=3
20    CONTINUE
      DO 21 I=1,NRCD
      SSL(I) = 0.
      QUS(I) = 0.
      TYPS(I) = 0.
      IHP(I) = 0.
   21 CONTINUE
      IPR=0
      INW = ID(2)
      DO 30 I=2,NRAY
      IF(I+1.LE.NRAY) INXT = ID(I+1)
      IF(I+2.LE.NRAY) INXTT = ID(I+2)
      IF(INW.GE.2) GO TO 29
      IQ1=0
      IQ2=0
      IQ3=0
      Z1=ZZ(I-1)
      S1=SEC(TGAM(I-1))*SS(I-1)
      DS=SEC(TGAM(I))*SS(I)-S1
      T1=TIME(I-1)
      DT=TIME(I)-T1
      TG1=TGAM(I-1)
      DTG=TGAM(I)-TG1
      DZ=ZZ(I)-Z1
C**********   TRY THIS FOR THAT 1.-IN-1.E6 CHANCE ***********
	IF(DZ.EQ.0.) GO TO 29
C************************************************************
      DO 25 J=1,NRCD
      ZR=RCD(J)
      IF(ZR.GT.ZB) GO TO 25
      F=(ZR-Z1)/DZ
C     DZ1=ZZ(I+1)-ZZ(I)
C     F1=(ZR-ZZ(I))/DZ1
      IF(F.LT.-.5 .OR. F.GT.1.5) GO TO 25
      IF(INW.NE.0 .OR. F.GE.1. .OR. F.LE.0.) GO TO 24
C     IF (INW.EQ.0 .AND. F.LE.1. .AND. F.GE.0. .AND. IHP(J).EQ.0) 23,24
C*****   AT ST. 23, RAYS I & I-1 BRACKET THE RECEIVER AND HAVE SAME
C   TURNAROUND HISTORY, SO WE HAVE A QUALITY 1 EIGENRAY.
   23 IQ1=1
      GO TO 26
   24 IF(F.LE.1. .OR. INW.NE.0 .OR. INXT.EQ.INW .OR. IHP(J).NE.0)
     1   GO TO 32
C*****   AT ST. 31 RAYS I & I-1 HAVE SAME TURNAROUND HISTORY BUT DO NOT
C   BRACKET THE RECEIVER; BUT THEY COME CLOSE AND NO QUALITY 1 EIGENRAY
C   HAS BEEN FOUND BETWEEN THEM OR WILL BE FOUND BETWEEN THE NEXT PAIR.
C   WE ACCEPT THIS AS A QUALITY 2 EIGENRAY.  NOTE THAT FOR QUALITY 2 WE
C   HAVE 1. < (ZR - Z(I-1))/(Z(I)-Z(I-1)) < 1.5, SO THAT WE ARE LOOKING
C   FOR A RECEIVER LOCATION ALMOST REACHED BY THE "ADVANCING" RAY FAN.
   31 IQ2=1
      GO TO 26
   32 IF(F.GE.0. .OR. IPR.EQ.INW .OR. INW.NE.0) GO TO 34
C*****   AT ST. 33 RAYS I & I-1 HAVE SAME TURNAROUND HISTORY AND COME
C   CLOSE TO BRACKETING RECEIVER; NO QUALITY 1 OR 2 EIGENRAY HAS BEEN
C   FOUND BETWEEN THEM AND NO QUALITY 1 EIGENRAY WAS FOUND BETWEEN THE
C   PREVIOUS PAIR.  WE ACCEPT THIS AS A QUALITY 3 EIGENRAY.  NOTE THAT
C   THIS RECEIVER LOCATION HAD   -.5 < (ZR-Z(I-1))/(Z(I)-Z(I-1)) < 0. ,
C   SO THAT IT WAS MISSED IN THE WAKE OF THE RAY FAN.
   33 IQ3=1
      GO TO 26
C*****   AT ST. 34 NO EIGENRAY HAS BEEN FOUND BECAUSE RAYS I & I-1 HAVE
C   DIFFERENT TURNAROUND HISTORIES.
   34 IQ1=0
      IQ2=0
      IQ3=0
      GO TO 25
C*****  AT ST. 26 AN EIGENRAY HAS BEEN FOUND; INTERPOLATE OR EXTRAPOLATE
C   IN DEPTH ITS TRAJECTORY PARAMETERS FROM THOSE OF RAYS I & I-1.
   26 SL=(S1+F*DS)/(RANGE*ABS(DZ))
      IF (SL.LE.0.) GO TO 25
C*****   THE NEXT STATEMENT IS TO IGNORE THIS EIGENRAY IF THE RAYS ARE
C   VERY CLOSE TOGETHER NEAR THIS RECEIVER; I.E., IF AT A CAUSTIC.
      IF (ABS(DZ).LT..001) GO TO 25
      TG=TG1+F*DTG
      IF (IQ1.EQ.1) IHP(J)=1
      TH=DEGREE*ATAN(TG)
      IF(INW.GE.1.AND.INXT.GE.1) GO TO 25
      IF(IFLMR.NE.0) FLM=2.*SIN(ZR*AK*TG/SQRT(1.+TG*TG))**2
      T=T1+F*DT
      IF (IQ1.EQ.1 .AND. SINT(J,1).EQ.0.) SINT(J,4)=FLM
      IF (IQ1.EQ.1 .AND. SINT(J,1).EQ.0.) SINT(J,1)=SL
      IF (IQ2.EQ.1 .AND. SINT(J,1).EQ.0. .AND.
     1             SINT(J,2).EQ.0. .AND. SINT(J,3).EQ.0.) SINT(J,4)=FLM
      IF (IQ2.EQ.1 .AND. SINT(J,1).EQ.0. .AND.
     1              SINT(J,2).EQ.0. .AND. SINT(J,3).EQ.0.) SINT(J,2)=SL
      IF (IQ3.EQ.1 .AND. SINT(J,1).EQ.0. .AND.
     1          SINT(J,2).EQ.0. .AND. SINT(J,3).EQ.0.) SINT(J,4)=FLM
      IF (IQ3.EQ.1 .AND. SINT(J,1).EQ.0. .AND.
     1           SINT(J,2).EQ.0. .AND. SINT(J,3).EQ.0.) SINT(J,3)=SL
      IF(IGNRAY.EQ.0) GO TO 22
      IF (IENT.NE.0) GO TO 41
	WRITE( IGNFIL, 1899) ITITLE
 1899 FORMAT ('1 EIGENRAY SET     ',20A4)
      WRITE( IGNFIL, 899)
      IENT = 1
41      IF(SL.LT.ALIM) GO TO 22
      PR=10.*ALOG10(SL)
      IF (IQ3.EQ.1) IQUAL=3
      IF (IQ2.EQ.1) IQUAL=2
      IF (IQ1.EQ.1) IQUAL=1
      IREM=IBLANK
      IF (IQ3.EQ.1 .AND. SINT(J,3).NE.0.) IREM=ISTAR
      IF (SINT(J,3).NE.0. .AND. IQ1.EQ.1) IREM=ISTARS
      RMKM = .001*RANGE
      WRITE( IGNFIL, 900) I, TH0(I), NBR(I), NTU(I), NSR(I), NTO(I),
     1                RMKM, ZR, TH, T, PR, IQUAL, IREM
899   FORMAT('0NRAY THETA0  NBR  NTU  NSR  NTO      RANGE(km)  DEPTH(m)'
     1   ,'  THETA(deg)  TIME(sec)', 4X,'SL(dB)',7X,'IQUAL',4X,'REM')
900   FORMAT(I5, F7.3, 4I5, F15.3, 2F10.2, F13.5, F10.1, 6X, I5, 5X, A2)
 1878 FORMAT(I5, F9.3, 4I5, 3F8.2, F12.5, F10.1)
   22 FLM=SINT(J,4)
      IF (SINT(J,1).NE.0.) SL=SINT(J,1)
      IF (SINT(J,2).NE.0. .AND. SINT(J,1).EQ.0.) SL=SINT(J,2)
      IF (SINT(J,3).NE.0. .AND. SINT(J,1).EQ.0.) SL=SINT(J,3)
      IF (SINT(J,3).EQ.0.) GO TO 27
      IF (IQ1.NE.1 .OR. SINT(J,1).EQ.0.) GO TO 72
   71 TYPE1R(J)=SSL(J)
      IF (ICOH.EQ.0) GO TO 27
      TYPE1C(J)=TYPS(J)
      QUAD(J)=QUS(J)
      GO TO 27
   72 SSL(J)=TYPE1R(J)
   27 TYPE1R(J)=TYPE1R(J)+SL*FLM
      IF(ICOH.EQ.0) GO TO 25
      PR=SQRT(SL)
      P=OMEGA*T+PHASE(I)
      IF (IQ3.EQ.1 .AND. IQ1.EQ.0) QUS(J)=QUAD(J)
      QUAD(J)=QUAD(J)+PR*SIN(P)
      IF (IQ3.EQ.1 .AND. IQ1.EQ.0) TYPS(J)=TYPE1C(J)
      TYPE1C(J)=TYPE1C(J)+PR*COS(P)
25    IF (INXT.GE.1 .OR. IPR.NE.INW .AND. INW.NE.INXT)IHP(J)=0
29    IPR=INW
      INW=INXT
C     IF (INW.EQ.0 .OR. INW.LT.1 .AND. INXTT.LT.1) GO TO 30
      IF(INW.EQ.0) GO TO 30
      DO 28 LL=1,NRCD
      SINT(LL,1)=0.
      SINT(LL,2)=0.
      SINT(LL,3)=0.
      SINT(LL,4)=0.
   28 CONTINUE
30    CONTINUE
35    DO 37 I=1,NRCD
      TYPE1C(I)=TYPE1C(I)*TYPE1C(I) + QUAD(I)*QUAD(I)
   37 CONTINUE
      IF(IT2.EQ.0) GO TO 60
      SW = 0.
      SDZ = 0.
      DO 38 I=2,NRAY
      W = AMIN1(SS(I),SS(I-1))
      SDZ = SDZ + W*ABS(ZZ(I) - ZZ(I-1))
      SW = SW + W
   38 CONTINUE
      DZBAR = WAVEL
      IF(SW.GT.0.) DZBAR = AMAX1(DZBAR,SDZ/SW)
      DZM = ZB/SQRT(FLOAT(NRAY))
      IF(DZBAR.GT.DZM) DZBAR = DZM
C*****   ASC WORD FORMAT LIMITS MAGNITUDE OF EXPONENTIAL ARGUMENT.
      DZBAR = AMAX1(DZBAR,.016*ZB)
      DO 385 J=1,NRCD
      IF(RCD(J).LE.ZB) QUAD(J) = EXP(-RCD(J)/DZBAR)
  385 CONTINUE
      F=2.*RANGE*DZBAR
      EB=EXP(-ZB/DZBAR)
      DO 50 I=1,NRAY
      IF(ABS(SS(I)).LT.ALIM)   GO TO 50
      EZ=EXP(-ZZ(I)/DZBAR)
      SL=SEC(TGAM(I))*SS(I)/F
      SGK=TGAM(I)*AK/SQRT(1.+TGAM(I)*TGAM(I))
      SL=SL/(1.-.5*(EB/EZ+EZ))
      DO 40 J=1,NRCD
      IF(RCD(J).GT.ZB) GO TO 50
      ER=QUAD(J)
      IF(IFLMR.NE.0) FLM=2.*SIN(RCD(J)*SGK)**2
      IF(ER.GT.EZ) GO TO 39
      S=SL*ER/EZ
      GO TO 40
39    S=SL*EZ/ER
40    TYPE2(J)=TYPE2(J)+S*FLM
50    CONTINUE
60    IF(IT3.EQ.0) RETURN
      DO 70 I=1,NRAY
      CALL CHANNL(ZZ(I),TGAM(I),ZTO,ZTU)
      CALL RCALC(ZZ(I),TGAM(I),ZTO,ZTU)
      IF(IFLMR.NE.0) CALL WDENS(ZZ(I),CTRR)
      SL=FVA*SS(I)/RANGE
      DO 65 J=1,NRCD
      IF(RCD(J).GT.ZTU) GO TO 70
      IF(RCD(J).LT.ZTO) GO TO 65
      CALL WDENS(RCD(J),S)
      IF(IFLMR.EQ.0) GO TO 64
      TG=TGAM(I)*CTRR/S
      FLM=2.*SIN(AK*RCD(J)*TG/SQRT(1.+TG*TG))**2
C**************************  VAX FLP OVERFLOW W/TG=4.5E19.
      S=S*FLM
64    TYPE3(J)=TYPE3(J)+S*SL
65    CONTINUE
70    CONTINUE
      RETURN
      END
 
      DOUBLE PRECISION FUNCTION DELB(DR, TGRP, TGRPN)
      REAL*8 DR, T0, T1, TGRP, TGRPN, T0T1
      T0 = TGRP
      T1 = TGRPN
      T0T1 = T0*T1
      IF(T0T1.NE.0.D0) DELB = -(1. + 1./T0T1)
      IF(T0.EQ.0.D0) DELB = 1./(T1*T1) - 1.
      IF(T1.EQ.0.D0) DELB = 1./(T0*T0) - 1.
      DELB = DR*DELB
      RETURN
      END

      SUBROUTINE CONNCT
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /PROF/ R1,N1,Z1(50),V1(50),R2,N2,Z2(50),V2(50)
      COMMON /TRIANG/ AP(100,2),BP(100,2),AL(100),BL(100),ZZERO(100),
     *   RZERO(100),AA(100),BB(100),SST(100),CCT(100),NTRI
      REAL*4 R1, Z1, V1, R2, Z2, V2
      NA=1
      NB=1
      NTRI=0
      IADD=0
10    NTRI=NTRI+1
      IF(NB.EQ.N2) GO TO 12
      IF (IADD.EQ.1.OR.NA.EQ.N1) GO TO 30
12    RZERO(NTRI)=R2
      ZZERO(NTRI)=Z2(NB)
      AL(NTRI)=Z1(NA+1)
      BL(NTRI)=(Z2(NB)-Z1(NA+1))/(R2-R1)
      AA(NTRI)=1.D0/V2(NB)**2
      BZ=(1.D0/V1(NA)**2-1.D0/V1(NA+1)**2)/(Z1(NA)-Z1(NA+1))
      BR=(1.D0/V2(NB)**2-1.D0/V1(NA)**2-BZ*(Z2(NB)-Z1(NA)))/(R2-R1)
      NA=NA+1
      IADD=1
      GO TO 40
30    RZERO(NTRI)=R1
      ZZERO(NTRI)=Z1(NA)
      AL(NTRI)=Z1(NA)
      BL(NTRI)=(Z1(NA)-Z2(NB+1))/(R1-R2)
      AA(NTRI)=1.D0/V1(NA)**2
      BZ=(1.D0/V2(NB)**2-1.D0/V2(NB+1)**2)/(Z2(NB)-Z2(NB+1))
      BR=(1.D0/V1(NA)**2-1.D0/V2(NB)**2-BZ*(Z1(NA)-Z2(NB)))/(R1-R2)
      NB=NB+1
      IADD=0
   40 BB(NTRI) = DSIGN(DSQRT(BZ*BZ + BR*BR), BZ)
C     IF(BB(NTRI).EQ.0.) GO TO 41
      IF(DABS(BB(NTRI)).LT.1.D-20) GO TO 41
      SST(NTRI)=BR/BB(NTRI)
      CCT(NTRI)=BZ/BB(NTRI)
      GO TO 42
41    CCT(NTRI)=1.D0
      SST(NTRI)=0.
42    AL(NTRI)=AL(NTRI)-BL(NTRI)*R1
      TANT=SST(NTRI)/DMAX1(CCT(NTRI),1.D-30)
      BP(NTRI,2)=(BL(NTRI)+TANT)/(1.-BL(NTRI)*TANT)
      B=BL(NTRI-1)
      IF(NTRI.EQ.1) B=0
      BP(NTRI,1)=(B+TANT)/(1.D0-B*TANT)
      AP(NTRI,1) = 0.
      AP(NTRI,2) = 0.
      IF(NA.LT.N1.OR.NB.LT.N2) GO TO 10
      RETURN
      END

      SUBROUTINE RAYOUT
      INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL
      PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7,
     &   DSTFIL = 37, IGNFIL = 41 )
      COMMON/RAYS/NRAY,TH0(1000),TGAM(1000),ZZ(1000),SS(1000),DSUM(1000)
     1   ,TIME(1000),PHASE(1000),NBR(1000),NTU(1000),NSR(1000),NTO(1000)
      COMMON /INFO/ROLD,RANGE,OMEGA,ATT,ICOH,IRAN,IT2,IT3,IGNRAY,IFLMR
      COMMON /TRIANG/ AP(100,2),BP(100,2),AL(100),BL(100),ZZERO(100),
     *   RZERO(100),AA(100),BB(100),SST(100),CCT(100),NTRI
      REAL*8 AP,BP,AL,BL,ZZERO,RZERO,AA,BB,SST,CCT
      COMMON /TLE/ ITITLE (18)
      COMMON/ABC/NBRS, NSRS, ALIM
      COMMON /RANST/ RDUM,IREC,IFSK,SDUM
      DATA IFSK/-1/
      DATA IENT/0/
      ZB = AL(NTRI) + BL(NTRI)*RANGE
C      IF (IENT.EQ.1) GO TO 1
C      WRITE( RAYFIL )ITITLE, ALIM, NBRS, NSRS
C      IENT=1
C      IREC=1
C      IFSK=IFSK+1
    1 WRITE( RAYFIL ) NRAY, RANGE, ZB,
     &   (TGAM(I),I=1,NRAY), (ZZ(I),I=1,NRAY),
     1   (SS(I),I=1,NRAY), (TIME(I),I=1,NRAY), (PHASE(I),I=1,NRAY),
     2   (NBR(I),I=1,NRAY), (NTU(I),I=1,NRAY), (NSR(I),I=1,NRAY),
     3   (NTO(I),I=1,NRAY)
      IREC=IREC+1
      RETURN
      END
 
      SUBROUTINE VELCAL
      COMMON /INFO/ RDUM, RANGE, OMEGA, IDUMM(7)
      COMMON /VELPRF/ N,Z(100),C(100)
      COMMON /TRIANG/ AP(100,2),BP(100,2),AL(100),BL(100),ZZERO(100),
     *   RZERO(100),AA(100),BB(100),SST(100),CCT(100),NTRI
      REAL*8 AP,BP,AL,BL,ZZERO,RZERO,AA,BB,SST,CCT
      COMMON/ABC/NBRS,NSRS,ALIM
      DATA ROLD/-1. E+20/
    3 IF(ABS(RANGE-ROLD).LT.1.) GO TO 20
      ROLD=RANGE
      N=NTRI+1
      DO 10 I=1,N
      IF(I.EQ.1) GO TO 1
      J=I-1
      ZZ=AL(J)+BL(J)*RANGE
      GO TO 2
1      J=1
      ZZ=0.
2     ZP=CCT(J)*(ZZ-ZZERO(J))+SST(J)*(RANGE-RZERO(J))
      C(I)=1.D0/DSQRT(AA(J)+BB(J)*ZP)
10    Z(I)=ZZ
20    RETURN
      END
 
      SUBROUTINE ADVNCE
      IMPLICIT REAL*8(A-H,O-Z)
      INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL
      PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7,
     &   DSTFIL = 37, IGNFIL = 41 )
      COMMON /TRIANG/ AP(100,2),BP(100,2),AL(100),BL(100),ZZERO(100),
     *   RZERO(100), AA(100), BB(100), SST(100), CCT(100), NTRI
      COMMON/RAYS/NRAY,TH0(1000),TGAM(1000),ZZ(1000),SS(1000),DSUM(1000)
     1   ,TIME(1000),PHASE(1000),NBR(1000),NTU(1000),NSR(1000),NTO(1000)
      COMMON /START/ CT0(1000)
      COMMON /INFO/RSTART,RMAX,OMEGA,ATT,ICOH,IRAN,IT2,IT3,IGNRAY,IFLMR
      COMMON /PIDEF/ PI, DEGREE, TWO PI
      COMMON /MIRROR/ SL, BRLT(200), BPST(200)
      COMMON/ABC/ NBRS, NSRS, ALIM
      DIMENSION RPNEW(4), RNEW(4), ZNEW(4), ZPNEW(4), D(7), RTNEXT(1000)
      REAL*4 TH0, TGAM, ZZ, SS, TIME, PHASE, RSTART, RMAX, OMEGA, ATT
      REAL*4 DEGREE, PI, TWO PI, BRLT, BPST, ALIM, DSUM, CT0, SL
      DATA ONE MM /.001D0/, FIVE MM /.005D0/, RTNEXT /1000*1.D0/
      DATA DLIMIT /1.D-30/
      LOGICAL ONUP,ONLW
      DELT(DR,T,U) = DABS(DR)*DSQRT(CMIS)*(1.D0+(T*(T+U)+U*U)/3.D0)
C        LOOP OVER ALL THE RAYS
      DO 100 IRAY=1,NRAY
      IHIST = 0
      IF(NBRS.EQ.100000 .AND. NSRS.EQ.100000) GO TO 10
      IF(NBR(IRAY).GT.NBRS .OR. NSR(IRAY).GT.NSRS) SS(IRAY) = .1*ALIM
   10 IF(ABS(SS(IRAY)).LT.ALIM) GO TO 100
      ZR = ZZ(IRAY)
      RR = RSTART
      TGR = TGAM(IRAY)
      TRAY = TIME(IRAY)
      DRDTH0 = DSUM(IRAY)
      DO 20 I=1,NTRI
C        FIND CORRECT LAYER
      P = ZR-AL(I)-BL(I)*RR
      IF(P.GT.ONE MM) GO TO 20
      NTRR = I
      IF(DABS(P).GT.ONE MM .OR. TGR.GT.-BL(I)) GO TO 30
20    CONTINUE
      GO TO 86
30    CT = CCT(NTRR)
C      IF(ITLC.NE.9) GO TO 301
C      STR = TGR/DSQRT(1.+TGR*TGR)
C      RTL = 0.
C      IF(RR.GT.0. .AND. DRDTH0.NE.0.)
C     1   RTL = 10.*DLOG10(DABS(CT0(IRAY)*SS(IRAY)/(RR*DRDTH0*STR)))
C      WRITE( PRTFIL, 901) IRAY, NTRR, IHIST, RR, ZR, TGR, TRAY,
C     1   RTL, PHASE(IRAY), NBR(IRAY), NTU(IRAY), NSR(IRAY), NTO(IRAY)
  301 ST = SST(NTRR)
      Z0 = ZZERO(NTRR)
C        TRANSFORM TO PRIMED COORDINATES
      R0 = RZERO(NTRR)
      ZRP = CT*(ZR-Z0)+ST*(RR-R0)
      RRP = CT*(RR-R0)-ST*(ZR-Z0)
      CIS = AA(NTRR)+BB(NTRR)*ZRP
      TANT = ST/DMAX1(CT,DLIMIT)
      TGRP = (TANT - TGR)/(1. + TANT*TGR)
      CMIS = CIS/(1.+TGRP*TGRP)
C        CALCULATE PARABOLIC RAY PATH.  ALPHA IS THE CURVATURE
      ALPHA = .25*BB(NTRR)/CMIS
      TWO ALP = 2.*ALPHA
C        FIND INTERSECTIONS WITH UPPER AND LOWER LAYER BOUNDARIES
      DO 40 I=1,2
      C = AP(NTRR,I) + BP(NTRR,I)*RRP - ZRP
      P = BP(NTRR,I) - TGRP
      IF(DABS(ALPHA).GT.1.D-25) GO TO 303
C        RAY PATH IS LINEAR FOR THESE STATEMENTS
      DRP = -C/P
      K = 2*I-1
      RPNEW(K) = RRP+DRP
      ZPNEW(K) = ZRP + TGRP*DRP
      ZNEW(K) = CT*ZPNEW(K) - ST*RPNEW(K) + Z0
      RNEW(K) = CT*RPNEW(K) + ST*ZPNEW(K) + R0
      ZNEW(K+1) = 1.D+30
      ZPNEW(K+1) = 1.D+30
      RNEW(K+1) = 1.D+30
      RPNEW(K+1) = 1.D+30
      GO TO 40
C        SOLVE QUADRATIC EQUATION FOR PARABOLIC RAY PATH
303   DISC = P*P+4.*ALPHA*C
      IF(DISC.LT.0.) GO TO 32
      DISC = DSQRT(DISC)
      DO 31 J=1,2
      K = 2*I+J-2
      S = 2*J-3
      F = P+S*DISC
      DRP = F/TWO ALP
C        ITERATION TO IMPROVE ACCURACY FOR SMALL CURVATURE RAYS
      IF(DABS(F).LT..1*DABS(P)) DRP=(ALPHA*((ALPHA*DRP*DRP-C)/P)**2-C)/P
      RPNEW(K) = RRP+DRP
      ZPNEW(K) = ZRP+DRP*(TGRP+ALPHA*DRP)
      ZNEW(K) = CT*ZPNEW(K) - ST*RPNEW(K) + Z0
31    RNEW(K) = CT*RPNEW(K) + ST*ZPNEW(K) + R0
      GO TO 40
32    DO 33 J=1,2
      K = 2*I + J-2
33    RNEW(K) = -1.D+30
40    CONTINUE
      ONLW= DABS(ZR-AL(NTRR)-BL(NTRR)*RR).LT.FIVEMM.AND.TGR.GT.-BL(NTRR)
C        SELECT CORRECT INTERSECTION AS NEXT POSITION
      IF(NTRR.EQ.1) GO TO 405
      ONUP = DABS(AL(NTRR-1)+BL(NTRR-1)*RR-ZR).LT.FIVE MM
     *   .AND. TGR.LT.-BL(NTRR-1)
      GO TO 41
405   ONUP = DABS(ZR).LT.FIVE MM .AND. TGR.LT.0.
41    IUP = -1
      IDN = -1
      IF(.NOT.ONUP) GO TO 411
      IUP = 1
      IF(DABS(RR-RNEW(1)).LT.DABS(RR-RNEW(2))) IUP = 2
      IF(RNEW(IUP).LT.RR .OR. RNEW(IUP).GT.RMAX) IUP = 0
411   IF(.NOT.ONLW) GO TO 412
      IDN = 3
      IF(DABS(RR-RNEW(3)).LT.DABS(RR-RNEW(4))) IDN = 4
      IF(RNEW(IDN).LT.RR .OR. RNEW(IDN).GT.RMAX) IDN = 0
412   IF(IUP.GE.0) GO TO 414
      IUP = 0
      RTRY = RMAX
      DO 413 I=1,2
      IF(RNEW(I).LT.RR .OR. RNEW(I).GT.RTRY) GO TO 413
      RTRY = RNEW(I)
      IUP = I
413   CONTINUE
414   IF(IDN.GE.0) GO TO 416
      IDN = 0
      RTRY = RMAX
      DO 415 I=3,4
      IF(RNEW(I).LT.RR .OR. RNEW(I).GT.RTRY) GO TO 415
      IDN = I
      RTRY = RNEW(I)
415   CONTINUE
416   IF(IDN.EQ.0) GO TO 420
      IF(IUP.NE.0) GO TO 418
417   ITRY = IDN
      RTRY = RNEW(IDN)
      GO TO 429
418   IF(RNEW(IUP).GE.RNEW(IDN)) GO TO 417
419   ITRY = IUP
      RTRY = RNEW(IUP)
      GO TO 429
420   IF(IUP.EQ.0) GO TO 50
      GO TO 419
  429 DELR = RPNEW(ITRY) - RRP
      TGRPN = TGRP + TWO ALP*DELR
      TRAY = TRAY + DELT(DELR, TGRP, TGRPN)
      DRDTH0 = DRDTH0 + DELB(DELR, TGRP, TGRPN)*TGRP*RTNEXT(IRAY)/CT
      TRP = TGRP
      IF(DABS(TGRP).LT.DLIMIT) TRP = DSIGN(DLIMIT,TGRP)
      TRPN = TGRPN
      IF(DABS(TGRPN).LT.DLIMIT) TRPN = DSIGN(DLIMIT,TGRPN)
      RTNEXT(IRAY) = RTNEXT(IRAY)*TRP/TRPN
      TGRN = (TANT - TGRPN)/(1. + TGRPN*TANT)
      RR = RTRY
      ZR = ZNEW(ITRY)
C        CHECK FOR SURFACE AND BOTTOM HITS, TURNOVERS AND TURNUNDERS
      IF(TGR*TGRN.GE.0.) GO TO 44
      IF(TGR.GT.0.) GO TO 43
      NTU(IRAY) = NTU(IRAY) + 1
      GO TO 44
   43 NTO(IRAY) = NTO(IRAY) + 1
44    IF(ITRY.GT.2) GO TO 45
      NTRR = NTRR-1
      IF(NTRR.GT.0) GO TO 46
      NTRR = 1
      SS(IRAY) = SS(IRAY)*SL
      PHASE(IRAY) = PHASE(IRAY) + PI
      NSR(IRAY) = NSR(IRAY) + 1
      TGR = -TGRN
      GO TO 30
45    NTRR = NTRR + 1
      IF(NTRR.LE.NTRI) GO TO 46
      NTRR = NTRI
      NBR(IRAY) = NBR(IRAY) + 1
      TGRAZE = (TGRN + BL(NTRI))/(1.D0 - TGRN*BL(NTRI))
      S = 100.D0*DABS(TGRAZE) + 1.D0
      N = S
      S = S-N
      IF(N.GT.199) N = 199
      SS(IRAY) = SS(IRAY)*((1.-S)*BRLT(N) + S*BRLT(N + 1))
      PHASE(IRAY) = PHASE(IRAY) + ((1.-S)*BPST(N) + S*BPST(N + 1))
      TGR = -(TGRAZE + BL(NTRI))/(1. - TGRAZE*BL(NTRI))
      IF(TGR.LT.-BL(NTRI)) GO TO 86
      GO TO 30
46    TGR = TGRN
      GO TO 30
C        CALCULATE RAY INTERSECTION WITH VERTICAL BOUNDARY, R = RMAX
   50 IF(DABS(ST).LT.2.D-4) GO TO 60
      BV = -CT/ST
      AV = (RMAX-R0)/ST
      C = AV + RRP*BV - ZRP
      P = BV-TGRP
      IF(DABS(ALPHA).GT.1.D-25) GO TO 501
C        LINEAR RAYS
      RPNEW(1) = RRP-C/P
      RPNEW(2) = 1.D+30
      GO TO 505
C        PARABOLIC RAYS
501   DISC = P*P+4.D0*ALPHA*C
      IF(DISC.LT.0.) GO TO 86
      DISC = DSQRT(DISC)
      S = 1.
      DO 503 I=1,2
      F = P + S*DISC
      DRP = F/TWO ALP
      IF(DABS(F).LT..1*DABS(P)) DRP=(ALPHA*((ALPHA*DRP*DRP-C)/P)**2-C)/P
      RPNEW(I) = RRP+DRP
503   S = -S
  505 DO 51 I=1,2
      ZPNEW(I) = ZRP + (RPNEW(I)-RRP)*(TGRP + ALPHA*(RPNEW(I)-RRP))
51    ZNEW(I) = CT*ZPNEW(I) - ST*RPNEW(I) + Z0
      I = 2
C        SELECT CORRECT INTERSECTION
      IF(DABS(ZNEW(1)-ZR).LT.DABS(ZNEW(2)-ZR)) I = 1
      DELR = RPNEW(I) - RRP
      TGRPN = TGRP + TWO ALP*DELR
      TRAY = TRAY + DELT(DELR, TGRP, TGRPN)
      DRDTH0 = DRDTH0 + DELB(DELR, TGRP, TGRPN)*TGRP*RTNEXT(IRAY)/CT
      TRP = TGRP
      IF(DABS(TGRP).LT.DLIMIT) TRP = DSIGN(DLIMIT,TGRP)
      TRPN = TGRPN
      IF(DABS(TGRPN).LT.DLIMIT) TRPN = DSIGN(DLIMIT,TGRPN)
      RTNEXT(IRAY) = RTNEXT(IRAY)*TRP/TRPN
      TGRN = (TANT - TGRPN)/(1. + TGRPN*TANT)
C        CHECK FOR TURNOVERS AND TURNUNDERS
      IF(TGR*TGRN.GE.0.) GO TO 54
      IF(TGR.LT.0.) GO TO 53
      NTO(IRAY) = NTO(IRAY) + 1
      GO TO 54
   53 NTU(IRAY) = NTU(IRAY) + 1
54    TGAM(IRAY) = TGRN
      IF(NTRR.LE.1) GO TO 545
C        CHECK THAT RAY IS WITHIN THE PROPER LIMITS
      IF(ZNEW(I).LT.AL(NTRR-1) + BL(NTRR-1)*RMAX-.1D0) GO TO 86
545   IF(ZNEW(I).GT.AL(NTRR) + BL(NTRR)*RMAX + .1D0) GO TO 86
C        THROW IN VOLUME ATTENUATION
      IF(ATT.NE.0)
     1  SS(IRAY)=SS(IRAY)*10.**(-.0001*ATT*(TRAY-TIME(IRAY))/DSQRT(CIS))
   55 ZZ(IRAY) = ZNEW(I)
      TIME(IRAY) = TRAY
      DSUM(IRAY) = DRDTH0
      GO TO 100
60    IHIST = 60
C        WHEN SMALL OR ZERO HORIZONTAL GRADIENTS, EQUATION FOR SECTION
C        50 IS SINGULAR.  THIS IS AN ITERATIVE SOLUTION FOR SUCH A CASE
C        GUESS VALUE FOR ITERATION
      ZNW = ZR-TGR*(RMAX-RR)
   66 DO 70 I=1,7
      RPNW = CT*(RMAX-R0)-ST*(ZNW-Z0)
      DRP = RPNW-RRP
      ZPNW = ZRP + DRP*(TGRP + ALPHA*DRP)
      ZNW2 = CT*ZPNW-ST*RPNW+Z0
      D(I) = DABS(ZNW2-ZNW)
      ZNW = ZNW2
   70 CONTINUE
C        CHECK FOR CONVERGENCE
      IF(D(7).LT..01D0) GO TO 68
      IF(D(7).LT..5*D(3)) GO TO 66
      IHIST = 66
      GO TO 86
68    RPNEW(1) = RPNW
      RPNEW(2) = RPNW
      GO TO 505
C        PRINT ERROR MESSAGE
86    SS(IRAY) = 0.
      WRITE( PRTFIL, 986) IRAY, IHIST, RR, ZR, TGR, RMAX, ZNEW(I), TGRN
      GO TO 55
  100 CONTINUE
      RETURN
  901 FORMAT(' RAY', I4, ' IN LAYER', I4, ' IHIST=', I3,
     1   '  R=', F10.2, '  Z=', F8.2, '  TG=', F10.5,
     2   '  T=', F10.4, '  L=', F6.1, '  P=', F8.3, '  NCTR=', 4I4)
  986 FORMAT(' RAY', I5, ' TERMINATED; IHIST=', I2, '   RR, ZR, TGR =',
     1   3D11.3, '   NEW R, Z, TG =', 3D11.3)
      END

      SUBROUTINE INTPRT 
      INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL
      PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7,
     &   DSTFIL = 37, IGNFIL = 41 )
      COMMON /INFO/ RSTART,RMAX,OMEGA,ATT,ICOH,IRAN,IT2,IT3,IGNRAY,IFLMR
      COMMON/LIMITS/ R1 , DR , R2 
      COMMON/ABC/NBRS,NSRS,ALIM
      COMMON /TLE/ ITITLE (18)
      DATA IENT/0/
      DIMENSION IS(4), ITT(4)
      EQUIVALENCE (ICOH,ITT)
      COMMON /LOUDNS/ NRCD,RCD(100),XINT(100,4)
	COMMON /TLPLOT/ RKM(1000), XTL(1000,8), ITYP(1000), IRCD(1000)
      DATA IS/1HC, 1HR, 1H2, 1H3/, NR/0/
      DATA IBLANK/'    '/
40    IF(IENT.EQ.1) GO TO 41
      WRITE( PRTFIL, 900) ITITLE
900    FORMAT ('1',20A4)
      IENT=1
      WRITE( PRTFIL, 905) (RCD(I),I=1,NRCD)
      WRITE( PRTFIL, 908)
905   FORMAT('0RECEIVED LEVEL VS RANGE',19X,'AT DEPTHS'/
     *   '0   R(KM)  TYPE',13F9.3,/,15X,13F9.3)
908   FORMAT(1X)
41	NR = NR + 1
	RKM(NR) = .001*RMAX
      NPLOT=0
	NPR = 0
      DO 50 I=1,4
      IF(ITT(I).EQ.0) GO TO 50
      DO 45 J=1,NRCD
      XINT(J,I) = 10.*ALOG10(AMAX1(ALIM  ,XINT(J,I)))
	NPLOT = NPLOT + 1
	XTL(NR,NPLOT) = ABS(XINT(J,I))
	ITYP(NPLOT) = I-1
	IRCD(NPLOT) = RCD(J) + .5
   45 CONTINUE
      IF(NPR.EQ.0) GO TO 46
C	WRITE(INTFIL, 909) RKM(NR), (XINT(J,I),J=1,NRCD)
      WRITE( PRTFIL, 906) IS(I),(XINT(J,I),J=1,NRCD)
906   FORMAT(13X,A1,1X,13F9.1,/,15X,13F9.1)
      GO TO 50
46	NPR = 1
C	WRITE(INTFIL, 909) RKM(NR), (XINT(J,I),J=1,NRCD)
      WRITE( PRTFIL, 907) RKM(NR),IS(I),(XINT(J,I),J=1,NRCD)
907   FORMAT(F9.3,4X,A1,1X,13F9.1,/,15X,13F9.1)
909	FORMAT(10F8.2)
   50 CONTINUE
      RETURN
      END
