      PROGRAM TRIPLT
      INTEGER PLTFIL, RAYFIL
      PARAMETER ( PLTFIL = 1, RAYFIL = 21 )
      COMMON /BOTTOM/  RB(8000), ZB(8000)
C
      DIMENSION PLTITL(20), CODES(5000), NRPLOT(65), ZR(1000), TR(1000),
     1   SS(1000), TIME(1000), PHASE(1000), NBR(1000), NTU(1000),
     2   NSR(1000), NTO(1000), R(8000)
      COMMON /TRAJCT/ ZSTOR(36000), TSTOR(36000)
      DATA ZSTOR/36000*-1./
      COMMON /HITS/ NBHITS, NSHITS
      INTEGER*2 NBHIT(36000), NSHIT(36000)
      COMMON /GRAPH/ R1KM, DRKM, R2KM, Z1, DZ, Z2, XLG, YHGT
      DATA YHGT /6./, TENSE/.1/
      REAL RBKM(8000)
C
900   FORMAT(20A4)
901   FORMAT(16I5)
  902 FORMAT(10F8.0)
  903 FORMAT(' PLOT AXIS RANGES(KM), DEPTH(M), LGTH/HGT(IN)=', 5F8.1)
  904 FORMAT(' DB, NSR, NBR LIMITS, TITLE FLAG=', F8.1, 8I5)
905   FORMAT(5H0NRAY,20I5)
  906 FORMAT(1X, I6, ' RAYS;', I10, ' BLOCKS OF LENGTH=', I6)
C
C*****   THIS PROGRAM READS THE RAY-HISTORY FILE CREATED BY TRIMAIN,
C   FOR001.DAT, AND PLOTS SELECTED RAY PATHS IN THE RANGE-DEPTH PLANE.
C*****   PLOT IS FROM RANGES R1KM TO R2KM AND DEPTHS 0 TO ZMAX.
C   HEIGHT IS 8 INCHES, LENGTH IN INCHES=XLG.
C   ALIM=TRANSMISSION LOSS AT WHICH A RAY IS DROPPED FROM PLOT.
C
      READ( PLTFIL, * ) R1KM, DRKM, R2KM, ZMAX, XLG, ALIM
C
C*****   NSRS AND NBRS ARE SURFACE AND BOTTOM-HIT LIMITS AT WHICH
C   TO DROP A RAY FROM THE PLOT.  NRAYS=NUMBER OF RAYS TO PLOT.
C
      READ( PLTFIL, * ) NSRS, NBRS, NRAYS
C
      R1 = 1000.*R1KM
      DR = 1000.*DRKM
      R2 = 1000.*R2KM
      IF(ALIM.EQ.0.) ALIM = 200.
      WRITE(6,903) R1KM, R2KM, ZMAX, XLG, YHGT
      WRITE(6,904) ALIM, NSRS, NBRS
      SSMIN = 10.**(-.1*ALIM)
      IF(XLG.LE.0..OR.XLG.GT.120.) STOP
      IF(ZMAX.LE.0.) STOP
      Z1 = 0.
      DZ = 1000.
      IF(ZMAX.LT.4000.) DZ = 500.
      IF(ZMAX.LT.2000.) DZ = 250.
      NDZ = (ZMAX + DZ - 5.)/DZ
      Z2 = NDZ*DZ
C      IF(NRAYS.EQ.0) GO TO 5
      READ( PLTFIL, * ) (NRPLOT(I),I=1,NRAYS)
C      GO TO 26
C5      IL = 1
C6      IH = IL + 15      
C      READ( PLTFIL, 901, END=10 ) (NRPLOT(I),I=IL,IH)
C      GO TO 15
C10    IH=IL-1
C      GO TO 25
C   15 WRITE(6,905) (NRPLOT(I),I=IL,IH)
C      IF(NRPLOT(IH).EQ.0) GO TO 20
C      IL=IL+15
C      GO TO 6
C20    IH=IH-1
C      IF(NRPLOT(IH).NE.0) GO TO 25
C      IF(IH.GT.0) GO TO 20
C25    NRAYS=IH
26      IF(NRAYS.LE.0 .OR. NRAYS.GT.64) STOP
      LBLK = 36000/NRAYS
      NRS = LBLK
      MBLK=(NRS-1)/LBLK+1
      WRITE(6,906) NRAYS, MBLK, LBLK
      IF(MBLK.GT.1) STOP
      NR=NRS
      M=1
      IREC=0
      ZBMIN = 10000.
      DO 500 IR=1,LBLK
         IREC=IREC+1
         IFEF=0
   33    READ( RAYFIL, END=36) NRAY,R(IREC), ZB(IREC), (TR(I),I=1,NRAY),
     1     (ZR(I), I=1, NRAY), (SS(I), I=1, NRAY), (TIME(I), I=1, NRAY),
     2      (PHASE(I),I=1,NRAY), (NBR(I),I=1,NRAY), (NTU(I),I=1,NRAY),
     3      (NSR(I),I=1,NRAY), (NTO(I),I=1,NRAY)
         RB(IREC) = R(IREC)
         IF(IREC.LT.5) TYPE *, R(IREC), ZB(IREC), ZR(3)
         GO TO 40
36       IFEF=1
         IF(NRS.LT.IREC-1) GO TO 501
   37    NRS = IREC-1
         NR=NRS
         GO TO 501
   40    IF(R(IREC).LT.R1) GO TO 33
         IF(IREC.EQ.2 .AND. R(1).EQ.R(2)) GO TO 33
         IF(R(IREC).GT.R2) GO TO 37
         ZBMIN = AMIN1(ZBMIN,ZB(IREC))
         DO 200 I=1,NRAYS
            K=NRPLOT(I)
C           ------ CHECK THAT THE RAY ACTUALLY EXISTS
            IF ( K .GT. NRAY ) THEN
               TYPE *, '*** WARNING: RAY ',K,' DOES NOT EXIST ***'
               NRAYS = I-1
               GOTO 500
            ENDIF
            J=LBLK*(I-1)+IREC
            ZSTOR(J) = -1.
            IF(SS(K).LE.SSMIN .OR. NSR(K).GT.NSRS .OR.
     &         NBR(K).GT.NBRS) GO TO 200
            ZSTOR(J) = ZR(K)
            NBHIT(J) = NBR(K)
            NSHIT(J) = NSR(K)
            TSTOR(J) = -TR(K)
200      CONTINUE
500   CONTINUE
  501 RMAX=R(NR)
C
C*****   THESE DISSPLA CALLS DRAW THE AXES AND THE BOTTOM PROFILE.
C
      CALL COMPRS
      CALL COMPLX
C      CALL SHDCHR(0.,1,.005,1)
C      CALL BASALF('STANDARD')
      CALL NOBRDR
      CALL PAGE(XLG+4., YHGT+3.)
      CALL INTAXS
      CALL YAXANG(0.)
      CALL XTICKS(2)
      CALL YTICKS(2)
      CALL XNONUM
      CALL HEIGHT(.25)
      CALL TITLE( 'TRIMAIN RAY PATHS$', 100, 'Range (km)$', 100,
     &   'Depth (m)$', 100, XLG, YHGT ) 
      CALL FRAME
      CALL GRAF(R1KM, DRKM, R2KM, Z2, -DZ, Z1)
      CALL GRACE(0.)
      CALL NOCHEK
      IF(ZBMIN.GT.ZMAX) GO TO 550
      CALL THKCRV(.02)
C     ------ CONVERT TO KM
      DO 507 IR = 1, NR
         RBKM(IR) = RB(IR)/1000.0
  507 CONTINUE
      CALL CURVE(RB, ZB, NR, 0)
      CALL RESET('THKCRV')
C
C*****   PLOT ONE RAY PATH EACH PASS THROUGH LOOP. 
C
550   CONTINUE
      DO 1000 I=1,NRAYS
         L=(I-1)*LBLK+1
         L1 = L + NRS 
         DO 800 J=1,NRS
            JJ = L1 - J
            IF(ZSTOR(JJ).LT.0) GO TO 800
            NPTS = JJ
            GO TO 801
800      CONTINUE
801      NBHITS = NBHIT(NPTS)
         NSHITS = NSHIT(NPTS)
         NPTS = NPTS - L + 1
         CALL RAYPLT(R, NPTS, ZSTOR(L), TSTOR(L), NBHIT(L), NSHIT(L))
 1000 CONTINUE
      CALL RESET('XNONUM')
      CALL ENDPL(0)
      CALL DONEPL
C
      STOP
      END

      SUBROUTINE RAYPLT(R, N, Z, T, NBHIT, NSHIT)
      COMMON /BOTTOM/  RB(8000), ZB(8000)
      COMMON /HITS/ NBHITS, NSHITS
      DIMENSION R(8000), Z(2000), T(2000)
      INTEGER *2 NBHIT(2000), NSHIT(2000)
      REAL RKM(8000)
C
C*****   SUBROUTINE TO PLOT THE TRAJECTORY OF A SINGLE RAY.
C   THE  N  POINTS ON THE PATH ARE STORED IN ARRAYS R=RANGE(KM),
C   Z=DEPTH(M), T=TANGENT W/HORIZ(+DOWN). RANGES EQUALLY SPACED.
C   NBHIT,NSHIT=COUNTS OF BOTTOM AND SURFACE HITS AT EACH RANGE.  
C   THE BOTTOM CONTOUR IS IN RB(KM) AND ZB(M); R(I)=RB(I) AT ENTRY.
C   THIS ROUTINE COMPUTES APPROXIMATE LOCATION OF EACH SURFACE OR
C   BOTTOM HIT IN ORDER TO PLOT A SHARP REFLECTION, WITH SMOOTH
C   SEGMENT BETWEEN REFLECTIONS.
C
      N1 = 1
C
C*****   IF THE RAY NEVER HITS BOUNDARY, INTERPOLATE IT ALL IN ONE PIECE.
C
      IF(NBHITS.GT.0 .OR. NSHITS.GT.0) GO TO 50
      NPT = N-1
      GO TO 1001
C
C*****   OTHERWISE GO THROUGH THE LOOP TO FIND THE END OF THE SEGMENT.
C   THE DIRECTION OF THE RAY CHANGES CONTINUOUSLY WITHIN A SEGMENT.
C
50    N11 = N1 + 1
      NPT = 0
      DO 1000 I=N11,N
         NPT = NPT + 1
         DR = R(I) - R(I-1)
         IF(NSHIT(I).EQ.NSHIT(I-1)) GO TO 200
C
C*****   SURFACE HIT
C
         Z0 = (DR*T(I)*T(I-1) - Z(I)*T(I-1) + Z(I-1)*T(I))/
     1                            (T(I) - T(I-1))
         RHIT = R(I) + (Z0-Z(I))/T(I)
         ZHIT = 0.
         ZSAVE = Z(I)
         TSAVE = T(I)
         R(I) = RHIT
         Z(I) = ZHIT
         T(I) = 0.
         NSHIT(I-1) = NSHIT(I)
         GO TO 1001
200      IF(NBHIT(I).EQ.NBHIT(I-1)) GO TO 1000
C
C*****   BOTTOM HIT
C
         TB = (ZB(I)-ZB(I-1))/DR
         IF(TB.NE.0.) GO TO 210
C
C*****   PROJECT INCOMING & OUTGOING RAYS TO THEIR INTERSECTION AND
C      USE THE INTERSECTION RANGE AS THE RANGE OF THE BOTTOM HIT.
C
         Z0 = (DR*T(I)*T(I-1) - Z(I)*T(I-1) + Z(I-1)*T(I))/
     1                            (T(I) - T(I-1))
         DRHIT = (Z(I)-Z0)/T(I)
         GO TO 220
C
C*****   PROJECT INCOMING & OUTGOING RAYS ONTO THE BOTTOM AND USE THE
C      MEAN OF THE INTERSECTION RANGES AS THE RANGE OF THE BOTTOM HIT.
C
210      RH1 = (ZB(I-1) - Z(I-1))/(T(I-1) - TB)
         RH2 = (ZB(I) - Z(I))/(TB - T(I))
         DRHIT = .5*(DR + RH2 - RH1)
220      ZHIT = ZB(I) - TB*DRHIT
         RHIT = R(I) - DRHIT
         ZSAVE = Z(I)
         TSAVE = T(I)
         R(I) = RHIT
         Z(I) = ZHIT
         T(I) = 0.
         NBHIT(I-1) = NBHIT(I)
         GO TO 1001
1000  CONTINUE
C
C*****   NPTS IS NUMBER OF POINTS ON THE SEGMENT TO BE PLOTTED.
C   NOTE THAT AN EXTRA POINT IS INSERTED AT EACH BOUNDARY HIT, AND
C   THE ORIGINAL COORDINATES ARE RESTORED TO R AND Z FOR NEXT SEGMENT.
C
1001  NPTS = NPT + 1
C     ------ CONVERT TO KM
      DO 901 IR = N1, N1+NPTS-1
         RKM(IR) = R(IR)/1000.0
  901 CONTINUE
      CALL CURVE( RKM(N1), Z(N1), NPTS, 0 )
      R(N1) = RB(N1)
      N1 = N1 + NPT
      IF(N1.GE.N .AND. T(N1).NE.0.) GO TO 1010
C
C*****   IF LAST POINT OF RAY PATH HAS NOT BEEN PLOTTED, THERE MUST 
C   HAVE BEEN A SURFACE OR BOTTOM HIT CAUSING EXIT FROM LOOP.
C
      N2 = N1
      N1 = N2 - 1
      R(N1) = RHIT
      Z(N1) = ZHIT
      T(N1) = 0.
      R(N2) = RB(N2)
      Z(N2) = ZSAVE
      T(N2) = TSAVE
      GO TO 50
1010  N1 = MIN0(N1,N)
      R(N1) = RB(N1)
C
      RETURN
      END
