C ----------------------------------------------------------------------
C			Program EIGEN
C
C
C  2-4 INPUT TEST DATA FILES, SHOULD HAVE NAMES 'DINP',TAPE1' & 'RYSFIL'
C  5 OUTPUT OF EIGEN PROGRAM FOR ELNINO INPUT
C   ONE WAY TO RUN PROGRAM IS BUILD LIBRARY FROM 'CORE' ROUTINES,
C   COMBINE WITH THE GRAPHICS WRITE ROUTINES AND RUN WITH THE EIGEN
C   ROUTINES. ALL THREE DATA FILES MUST BE PRESENT.
C ----------------------------------------------------------------------
c.....Some History:
c     As with HARPO, PROFILE and CONPLT, EIGEN was written with pre-FORTRAN 77
c     FORTRAN for a Cyber NOS environment.  Many modifications have been made
c     in order to port this code to an IBM PC-compatible DOS machine.  The
c     original code has usually been commented out rather than deleted when
c     possible.  These comments can usually be identified by one of the
c     following:
c
c     1. "ccccc"
c     2. "c" followed by one or more "*"
c     3. "crb"
c
c     In addition, new comments have usually been written using lowercase.
c     New source code also has been usually been written in lowercase, although
c     many changes in the small routines which write to NDEVTXT e.g. MARKER,
c     are in uppercase.
c
c
      PROGRAM EIGEN
	
	implicit double precision (a-h,o-z)
 
      REAL*8 AZIMUTH,ELEVA,RANGE
      LOGICAL XLOR,WSIGAP,INGAP,HADEIS
C
C     COMMON DECK "HDR" INSERTED HERE
      CHARACTER*10 INITID*80,DAT,TOD
      COMMON/HDR/SEC
      COMMON/HDRC/INITID,DAT,TOD
C     COMMON DECK "WWE" INSERTED HERE
      PARAMETER (NWARSZ=1000)
      character*80 id
      common /ident/ id
      COMMON/WW/dum(10),MAXW,W(NWARSZ)
      REAL*8 NRSLMI
      character*132 rbtemp
      EQUIVALENCE (W(80),APOG) ,(W(401),PTYP) ,(W(403),ELMIN)
     1   ,(W(404),ELMAX) ,(W(405),ELVINC) ,(W(406),RGINC)
     2   ,(W(407),HOPLME) ,(W(408),NRSLME)
     3   ,(W(409),GAPLO), (W(410),GAPHI), (W(411),SYMRED)
     4   ,(W(428),GRPMIN) ,(W(429),GRPMAX) ,(W(430),GRPVINC)
     5   ,(W(432),HOPLMG) ,(W(433),NRSLMG)
     6   ,(W(434),RFSPEED) ,(W(435),TMFT) ,(W(436),RMAX)
     7   ,(W(453),RCRNG) ,(W(454),SLPCUT)
     8   ,(W(456),HOPLMI) ,(W(457),NRSLMI) ,(W(458),ESET)
      REAL*8 MAXSTP,MAXERR,INTYP,LLAT,LLON
      EQUIVALENCE (EARTHR,W(1)),(RAY,W(2)),(XMTRH,W(3)),(TLAT,W(4)),
     1 (TLON,W(5)),(OW,W(6)),(FBEG,W(7)),(FEND,W(8)),(FSTEP,W(9)),
     2 (AZ1,W(10)),(AZBEG,W(11)),(AZEND,W(12)),(AZSTEP,W(13)),
     3 (BETA,W(14)),(ELBEG,W(15)),(ELEND,W(16)),(ELSTEP,W(17)),
     8   (RUNSUP,W(18)),(RCVRH,W(20)),
     4 (ONLY,W(21)),(HOP,W(22)),(MAXSTP,W(23)),(PLAT,W(24)),(PLON,W(25))
     5,(HMAX,W(26)),(RAYFNC,W(29)),(EXTINC,W(33)),
     6 (HMIN,W(27)),(RGMAX,W(28)),
     8 (INTYP,W(41)),(MAXERR,W(42)),(ERATIO,W(43)),
     6 (STEP1,W(44)),(STPMAX,W(45)),(STPMIN,W(46)),(FACTR,W(47)),
     7 (SKIP,W(71)),(RAYSET,W(72)),(PRTSRP,W(74)),(HITLET,W(75))
     9 ,(BINRAY,W(76)),(PAGLN,W(77)),(PLT,W(81)),(PFACTR,W(82)),
     1 (LLAT,W(83)),(LLON,W(84)),(RLAT,W(85)),(RLON,W(86))
     2,(TIC,W(87)),(HB,W(88)),(HT,W(89)),(TICV,W(96))
      EQUIVALENCE (W(153),REFC)
C     COMMON DECK "CONST" INSERTED HERE
      COMMON/PCONST/CREF,RGAS,GAMMA
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10
C     COMMON DECK "RAYDEV" INSERTED HERE
C     DEVICE ASSIGNED TO RAYTRC INPUT FILE
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT
 
C     COMMON DECK "CRAYCOM" INSERTED HERE
      COMMON/RAYCOM/NRAYS,MXRYS,MXHOP,MXXHOP,NRSLIM,ICODE
      COMMON/RAYCOM/NHOPS(333),MNHOPS(333)
C     COMMON DECK "CRAYDAT" INSERTED HERE
      COMMON/RAYDAT/AZIMUTH(333),ELEVA(333),RANGE(10000)
      COMMON/RAYDAT/GRPTIM(10000)
      COMMON/RAYDAT/NELMTS
      COMMON/GRPC/GRPN,GRPX
C
      INTEGER NCODES(3),PLTCOD
      CHARACTER INITRY*80
      CHARACTER*3 TESTID,PRVID
      PARAMETER (NHSZ=154,NWSSZ=400)
      REAL*8 WH(NHSZ),WS(NWSSZ)
      REAL*8 EL(333),RG(333),TT(333),GTIM(333)
C     COMMON DECK "CODES" INSERTED HERE
      COMMON/CODES/ICODRE,ICODRT,ICODEI
C
      CHARACTER TDID*30
      COMMON/CDATIM/TDID
	character cdate*10,ctime*10

C.....Added the following for input file OPENing.
      LOGICAL found
      CHARACTER*80 fname

C
      EQUIVALENCE (NCODES,ICODRE)
C
C     COMMON DECK "CEFILES" INSERTED HERE
      COMMON/EFILES/INFILE,IOTFIL,LOGFIL,NCNTRL,NRYSFL
C
      DATA NPLTS,NPPLT,NEIGS/3*0/
      DATA NTRAYS,INPLT,UNDEF/0,1,1.E25/
      DATA PRVID,LMXRYS,NPR/' ',0,0/
C
C+
C+      GENERAL INITIALIZATION OF SUPPORT SYSTEM AND PARAMETERS
C+
      CALL SYSINI
      CALL CONBLK
C
      INFILE=1
      IOTFIL=3
      NDTMP=4
      NDEVTMP=4
ccccc      NDEVGRP=10
      LOGFIL=7
      NCNTRL=8
      NRYSFL=9
C*      NPUNAM=12
      ICODRE=16
      ICODRT=32
      ICODEI=64
      NELMTS=10000
C+
C+I     PRIMARY FILE INITIALIZATION
C+
ccccc
      INQUIRE (FILE = 'RYSFIL', EXIST = found)
      if (found) then
           OPEN(UNIT=NRYSFL,FILE='RYSFIL',STATUS='OLD',IOSTAT=IST)
           IF (IST.LT.0) THEN
	        STOP 'Unexpected end of input rayset file during OPEN.'
           ELSE IF (IST.GT.0) THEN
                PRINT 444,IST
444	        FORMAT ('ERROR OPENING INPUT FILE.  ERROR:  ',I3)
	        STOP 'Check error code in FORTRAN Reference Manual.'
           ENDIF
      else
           WRITE (*,*)' No RYSFIL file exists for EIGEN program.'
30         WRITE (*,222)
222        FORMAT ('  Enter filename on next line >>>>>>')
           READ (*,333) fname
333        FORMAT (A)
           INQUIRE (FILE = fname, EXIST = found)
           if (.not. found) then
                WRITE (*,*)' Could not find file: ',fname
                go to 30
           endif
           OPEN(UNIT=NRYSFL,FILE= fname,STATUS='OLD',IOSTAT=IST)
           IF (IST.LT.0) THEN
	        STOP 'Unexpected end of input rayset file during OPEN.'
           ELSE IF (IST.GT.0) THEN
                PRINT 444,IST
	        STOP 'Check error code in FORTRAN Reference Manual.'
           ENDIF
           WRITE (*,*)' Using file: ',fname(1:len(fname))
      endif

ccccc
ccccc      OPEN(NRYSFL,FILE='RYSFIL',STATUS='OLD',ERR=2500)
      REWIND NRYSFL
C+
C+I     READ ORIGINAL DATE-TIME INFO
C+
      READ(NRYSFL,320,END=10000) CDATE,CTIME
320   FORMAT(T80,A10,2x,a10)
10000 IF(CDATE.EQ.' ') then
ccccc           CALL DATE(CDATE)
           continue
      endif
      IF(CTIME.NE.' ') GO TO 310
C+
C+D1    NON-NULL INFO?
C+
C  IN DEFAULT OF NOT FINDING A TIME VALUE, USE THE FILE NAME
C  ASSUMED TO BE IN FILE 'PUNAM'.
C+
C+R11   (NO) ATTEMPT TO READ ALTERNATE FILE WITH NEW DATE-TIME INFO
C+
C Until I discover what file PUNAM shoud be, CALL TIME...
ccccc	call time(ctime)
 
C*      OPEN(NPUNAM,FILE='PUNAM',STATUS='OLD',ERR=2600)
C*      REWIND NPUNAM
C*      READ(NPUNAM,315,END=310  ) CTIME
315   FORMAT(BZ,2(A10,2x))
C
C+
C+R12   (YES)
C+
C+    IDENTIFY OURSELVES TO 'READW' FOR SPECIAL MODE
C+    FOR WHICH ONLY TABULAR DATA IS WRITTEN TO TEMP FILE
C+    AND TABULAR ARRAYS ARE NOT NEEDED
C+
310   NFRMAT=-2
C
C+
C+I     OPEN THE FOLLOWING FILES
C+              'TAPE5'         GRAPHICS
C+              'LOGOUT'        SUMMARY INTERPOLATION INFO
C+              'DINP'          CONTROL INPUT FILE
C+              'TAPE1'         ORIGINAL RAY TRACE INPUT FILE
C+              'DOUTP'         OUTPUT TO BE USED AS INPUT TO RAY TRACE
C+              'TAPE4'         TEMPORARY WORK FILE
C+
ccccc      CALL OPNURP(NDEVGRP,'TAPE5')
      CALL OPNREP(LOGFIL,'LOGOUT')
c..........The following OPEN is performed in CONBLK instead of here.
ccccc      OPEN(NCNTRL,FILE='DINP',STATUS='OLD',ERR=2700)
ccccc      NRYIND=NCNTRL
ccccc      REWIND NCNTRL
      OPEN(INFILE,FILE='TAPE1',STATUS='OLD',IOSTAT=ios)
      if (ios .ne. 0) then
           print*,'Error opening TAPE1.  Error = ',ios
           stop
      endif
c..........DOUTP is now attached to logical unit 43.  All references to
c          logical unit IOTFIL have been commented out.
ccccc      CALL OPNREP(IOTFIL,'DOUTP')
C
 998  IF(MAXW.GT.500) then
crb   IF(MAXW.GT.500) then
           print*,'ERROR:  MAXW > 500'
           stop
      endif !!!!!CALL STOPIT('MAXW.GT.500')
ccccc      READ(NRYIND,10,END=10001) INITID
 
10    FORMAT(A)
 
      call stdint
 
ccccc 10001 WRITE(IOTFIL,10) INITID
      DO 15 I=1,NWSSZ
15    WS(I)=0.0
C
C+
C+      INITIALIZE W-ARRAY AND SET DEFAULT RUN PARAMETER VALUES
C+
ccccc      CALL SETW
C
      HOPLMI=30
      HOPLME=30
      HOPLMG=30
      SLPCUT=.1
      RFSPEED=1.5
      RGINC=5.
      RMAX=30.
      RCRNG=26.5
      NRSLMI=100.
      NRSLME=100.
      NRSLMG=100.
      TMFT=1000.
      LRUN=0
C
C LOOP UNTIL ALL RUNSETS ARE READ
C  Begin Main DO loop --------------------------------------------------
 
      DO 650 NRUN=1,10000
C
C+
C+    SWITCH TO CONTROL FILE FOR W-ARRAY INPUT
C+
      NRYIND=NCNTRL
      IF(READW(0).NE.0.) GO TO 670
C
      WRITE(TESTID,'(A3)') ID(1:3)
      WRITE(TDID,330) ID(1:3),CDATE,CTIME
330   FORMAT(' ',A3,2(A10,1x))
C
C+    DETERMINE VALIDITY OF FUNCTION CODES SPECIFIED
      PRINT 10,TDID
Csp      ICODE=ND2B(INT(RAYFNC))
ccccc      ICODE=ND2B(jidint(RAYFNC))
      icode = int(rayfnc)
C
C
      PLTCOD=IOR(ICODRE,ICODRT)
      IALL=IOR(PLTCOD,ICODEI)
C
C     IF W(29) IS ZERO DO ALL FUNCTIONS
C*      IF(AND(ICODE,IALL).EQ.0) ICODE=IALL
      IF(IAND(ICODE,IALL).EQ.0) ICODE=IALL
C
      PRINT 111, nrun,rayfnc
111   format(1x,'PROCESSING RUNSET #',I2,' FOR RAY FUNCTION ',f14.6)
C
      N=0
      DO 200 I=1,NHSZ
           IF(I.EQ.100) N=401-I
200        WH(I)=W(I+N)
C+
C+    USE MAXIMUM HOP LIMIT OF 3 CODES TO DETERMINE RAY
C+    DIMENSION OF RAY-HOP ARRAY.
C+
csp      MXHOP=AMAX1(HOPLMI,HOPLME,HOPLMG)
      MXHOP=dMAX1(HOPLMI,HOPLME,HOPLMG)
      MXXHOP=MXHOP
      MXRYS=MIN0(333,NELMTS/MXHOP)
C
C+
C+      SWITCH TO ORIGINAL DINP FILE FOR W-ARRAY INPUT
C+
      NRYIND=INFILE
      REWIND NRYIND
      READ(NRYIND,10,END=10002) INITRY
10002 REWIND NRYSFL
C+
C+      READ TO CURRENT RUN-SET OR THAT SPECIFIED IN 'ESET' PARAM
C+
      NSET=NRUN
      IF(ESET.NE.0.0) NSET=ESET
C
      NRYSET=0
      DO 400 NS=1,NSET
           IF(READW(0).ne.0.) then
                PRINT 410,TESTID
410             FORMAT('CANNOT LOCATE RUN ID ',A3,' IN RUNSET DATA')
                STOP
           else
               IF(REFC.GT.0.) CREF=REFC
           endif
C
C+
C+         COMPUTE CUMULATIVE COUNT OF RAYSETS PRODUCED AND
C+         USE FOR CURRENT RAYSET NUMBER
C+
           IF(RAYSET.NE.0.0) NRYSET=NRYSET+1
400   CONTINUE
C
425   DO 20 I=1,NELMTS
20          RANGE(I)=UNDEF
C
C OVERRIDE WITH CONTROL FILE VALUES
      N=0
      DO 220 I=1,NHSZ
      IF(I.EQ.100) N=401-I
220   W(I+N)=WH(I)
      PRINT *,'PROCESSING RAY SET #',NRYSET
C+
C+    READ RAYSET VALUES INTO ARRAYS
C+
 
      IF(GETRAYS(NRYSFL,NRYSET).NE.0.0) GO TO 670
      PRVID=TESTID
      LMXRYS=MXRYS
C
230   PRINT 235, NRAYS,sngl(AZIMUTH(1))
C+
C+    REPORT NUMBER OF RAYS READ AND AZIMUTH
C+
      WRITE(LOGFIL,240) ID,NRAYS,sngl(AZIMUTH(1))
235   FORMAT(1x,I4,' RAYS PROCESSED, AZIMUTH=',F9.3)
240   FORMAT(' ',A80/I4,' RAYS PROCESSED, AZIMUTH=',F9.3)
      IF(NRAYS.LE.0) GO TO 650
C
C
C FORCE AN OVERLAY MODE PLOT
      PLT=-dABS(PLT)
      ELSTEP=0.
      RAYFNC=0.
C
C     SKIP THOSE WITH FEWER THAN 2 RAYS. CHECK HEIGHT
C     AND LIMIT TO VALUES GREATER THAN A MINIMUM.
C     SUBMIT VECTORS TO PLOTTING PROGRAM AND
C
      GRPN=1.E25
      GRPX=-GRPN
C
C     DEFINE MIDDLE AND WIDTH OF THE GAP IN WHICH TO DO POINT PLOTTING
C+
C+    SET UP POINT PLOT FEATURE IF GAP IS SPECIFIED.
C+
      ELMID=.5*DEGS*(GAPLO+GAPHI)
      ELWID=.5*DEGS*(GAPHI-GAPLO)
C   REQUIRE A NON-ZERO RANGE FOR ANY EFFECT OF THIS FEATURE
      IF(ELWID.LE.0.0) ELWID=-1.
C
C+
C+    LOOP THROUGH 3 POSSIBLE PROGRAM MODES, PERFORM THOSE
C+    REQUESTED BY RAYFNC W-ELEMENT
C+
 
C  Begin loop over 3 program modes -------------------------------------
 
      DO 640 MODE=1,3
      NCODE=IAND(ICODE,NCODES(MODE))
      IF(NCODE.EQ.0) GO TO 640
      IF(NCODE.EQ.ICODRE) THEN
         PRINT *,'******** RANGE-ELEVATION PLOTS ********'
         MXHOP=HOPLME
         NRSLIM=NRSLME
      ELSEIF(NCODE.EQ.ICODRT) THEN
         PRINT *,'******** RANGE-DELAY PLOTS ********'
         MXHOP=HOPLMG
         NRSLIM=NRSLMG
      ELSEIF(NCODE.EQ.ICODEI) THEN
         PRINT *,' **EIGEN-RAY ESTIMATION'
         MXHOP=HOPLMI
         NRSLIM=NRSLMI
      ENDIF
C
      MXXHOP=MXHOP
      MXRYS=MIN0(333,NELMTS/MXHOP)
ccccc      PRINT *,' MAXIMUM NUMBER OF RAYS/HOPS= ',MXRYS,'/',MXXHOP
      NPPLT=0
C
      NBRKS=0
C+
C+    LOOP THRU SET OF HOPS
C+
C Begin loop over maximum hops -----------------------------------------
 
 
      DO 625 NHOP=1,MXHOP
         NCRAY=0
         INGAP=.FALSE.
         WSIGAP=.FALSE.
C
430      NPTS=0
         IF(NCRAY.GT.0) NBRKS=NBRKS+1
450      IF(.NOT.XLOR(INGAP,WSIGAP)) NCRAY=NCRAY+1
 
            RR=0.0
            IF(NCRAY.GT.NRAYS) then
                 GO TO 550
            endif
            IF(NHOP.GT.NHOPS(NCRAY).AND.NPTS.LT.1) then
                 GO TO 450
            endif
            if (nhops(ncray) .eq.9999) then
                 goto 500
            endif
            LNHOP=(NCRAY-1)*MXXHOP+NHOP
            RR=RANGE(LNHOP)
 
            IF(RR.EQ.UNDEF) then
                 GO TO 500
            endif
C+
C+       DETERMINE IF THIS POINT IS IN THE GAP
C+
            ELV=ELEVA(NCRAY)
            WSIGAP=INGAP
 
            INGAP=dABS(ELV-ELMID).LT.ELWID
C+
C+       IF IT WASNT AND NOW IS OR VICA VERSA,
C+          WE HAVE A PLOT MODE CHANGE AND MUST EXIT HERE
C+
            IF(XLOR(INGAP,WSIGAP))then
               GO TO 500
            endif
C
            NPTS=NPTS+1
            EL(NPTS)=ELV
            RG(NPTS)=RR
            GTIM(NPTS)=GRPTIM(LNHOP)
            T=(GRPTIM(LNHOP)/CREF-RR/RFSPEED)*TMFT
            TT(NPTS)=T
csp            GRPN=AMIN1(T,GRPN)
csp            GRPX=AMAX1(T,GRPX)
            GRPN=dMIN1(T,GRPN)
            GRPX=dMAX1(T,GRPX)
            GO TO 450
C  MOVE PAST ANY GROUPS OF UNDEFINED RANGES
500      IF(NPTS.LT.1) then
            GO TO 430
         endif
550      IF(NPTS.LT.1) then
            GO TO 625
         endif
C
      IF(NCODE.NE.ICODEI) GO TO 480
C
25    FORMAT(I3,G14.7)
C
C+D2  IS A NEW RUN SET COMMENCING?
C
      IF(NRUN.GT.LRUN) THEN
C+R21
C+D3     WERE THERE EIGEN RAYS PRODUCED LAST RUNSET?
         IF(HADEIS) THEN
C+R31       IF SO TERMINATE THE PLOTS
ccccc            IF(NRUN.GT.1) WRITE(IOTFIL,25) 81,sngl(DABS(PLT)),0
         ELSE
C+R32       OTHERWISE CANCEL THE LAST RUN SET(FOR THE RAY TRACE PROGRAM
C+          BY SETTING W(29) TO NON-ZERO BUT NON-SPECIFIC CODE(.1)
C+          INSURE THIS IS ONLY DONE ONCE BY TESTING WS(29).
C
C NO write to DOUTP.DAT
C***            IF(WS(29).NE.0.1) WRITE(IOTFIL,25) 29,.1,0
C           NOTIFY THE W-ARRAY OF SAME SO NEXT RUN SET WILL
C           GET A RESET W(29).
            WS(29)=.1
         ENDIF
C+       RESET THE 'HAD EIGEN RAYS' FLAG
         HADEIS=.FALSE.
C+       WRITE ALL W-ARRAY DATA AND TABLULAR DATA TO DINP FILE
C+       ONLY DIFFERENCES FROM PREVIOUS RUN SET ARE WRITTEN OUT
C
ccccc         CALL WRITW(IOTFIL,NDTMP,ID,W,WS,NWSSZ)
ccccc         WRITE(IOTFIL,25) 81,sngl(PLT)
crb   This is where eigen creates an output file "DOUTP"
crb   for use in rerunning HARPO to produce desired eigenrays
crb
c.....NOTE:  ron brent refers to "runsets" as "cases". jah
c
crb   It is crucial that cases be deliminated by the line
crb     "      ***********END OF RUNSET 1*************  "
crb     "      ***********END OF RUNSET 2************* " etc.
crb   in order for the program to work.
crb
crb   First one rewinds TAPE1
      rewind infile
crb   Open the file DOUTP
      call opnrep(43,'doutp')
crb   RBCOUNT is a variable that keeps track of which case one has
crb   read.  W(458) is the case number that EIGEN processes.  Once
crb   one has hit that case the write to DOUTP ceases.
crb
      rbcount=0
      do 255 i=1,5000
           read(infile,999,err=256)rbtemp
           if( rbtemp(25:27) .eq. '***' ) rbcount=rbcount+1
crb........If we have just read the case EIGEN processes one exits the loop
           if ( rbcount .eq. w(458) ) goto 256
crb   .....If we have not reaches the case number of interest we have to
crb   .....eliminate flags that run previuos cases.  This is done by NOT
crb   .....writing the "0" flag as well as the followinf ID line
           if( rbtemp(25:27) .eq. '***' ) then
                read(infile,999,err=256)
                goto 255
           ENDIF
c..........Added this IF statement to ensure that ID is written to the 2nd
c          line of DOUTP.  ID is a string containing the line describing the
c          data set e.g. "EL1-1    EL NINO PERTURBATION MODEL".  It is taken
c          from the appropriate runset in TAPE1.  The RB changes did not
c          write the correct ID to DOUTP. jah.

           if (i.eq.2) then
                write(43,257) id
           else
                write(43,999)rbtemp
           endif
 255  continue

 999  format(a132)
crb   This is where various flags are written to DOUTP

 256     write(43,25) 81,sngl(plt)
         write(43,25) 17,0.
         write(43,25) 29,0.
         LRUN=NRUN
         NPR=0
      ENDIF
257   format (a80)

C+R22
      N=NPR
C+    DO INTERPOLATION OF EIGEN RAYS
ccccc 475   CALL VINTERP(IOTFIL,EL,RG,GTIM,NPTS,NHOP,NPR)

c.....NOTE:  VINTERP is now called with '43' instead of IOTFIL in order
c            to take advantage of Ron Brent's modifications. jah.

475   CALL VINTERP(43,EL,RG,GTIM,NPTS,NHOP,NPR)
crb   write(43,25) 0
      NEIGS=NEIGS+NPR-N
      IF(NPR.NE.0) HADEIS=.TRUE.
      GO TO 600
C
C  Plot rays
 
480   IF(INPLT.NE.0) CALL PLOTINI(.3,SYMRED)
      INPLT=0
C HOP NO 1 IS THE FIRST CURVE OF THE PLOT
      IF(NPPLT.EQ.0) CALL PLTTLT(NCODE)
      MARK=0
C+    SET POINT MODE IF IN GAP
      IF(WSIGAP.OR.NPTS.EQ.1) MARK=3
C+    DO R-EL PLOT
      IF(NCODE.EQ.ICODRE)
     1   CALL KURVE(EL,RG,NPTS,MARK,DEGS*ELMIN,DEGS*ELMAX)
C+    DO DELAY-RANGE PLOT
      IF(NCODE.EQ.ICODRT)
     1  CALL KURVE(TT,RG,NPTS,MARK,GRPMIN*TMFT,GRPMAX*TMFT)
C
600   NPPLT=NPPLT+NPTS
 
C.....Get next ray for the current hop
      IF(RR.EQ.UNDEF .OR. XLOR(INGAP,WSIGAP) .or.
     &     nhops(ncray).eq.9999) GO TO 430
 
625   CONTINUE
C ----------------------- End loop over maximum hops -------------------
 
C
      NTRAYS=NTRAYS+NRAYS
ccccc      IF(NCODE.EQ.ICODRT) PRINT *,' GRPN, GRPX=',GRPN,GRPX
ccccc      IF(NCODE.EQ.ICODRE)
ccccc     1     PRINT *,' ELMIN, ELMAX=',DEGS*ELMIN,DEGS*ELMAX
C
      IF(NCODE.EQ.ICODEI) GO TO 640
      CALL endpl(0)
 
      NPLTS=NPLTS+1
 
640   CONTINUE
C ------------------- End loop over 3 program modes --------------------
 
C
ccccc      IF(NBRKS.GT.0) PRINT *,' NUMBER OF BREAKS IN CURVES= ',NBRKS
C
650   CONTINUE
 
C ---------------- End of Main loop over runsets -----------------------
 
C
crb670     write(43,25) 0
670   IF(NPLTS.GT.0) CALL DONEPL
      IF(NCODE.EQ.ICODEI) THEN
         PRINT *,' NUMBER OF EIGENRAYS GENERATED=',NEIGS
         PRINT *,'   EIGENRAY RANGE=',RCRNG,' AT HEIGHT=',RCVRH
C
ccccc         IF(NEIGS.GT.0) THEN
ccccc            WRITE(IOTFIL,25) 81,sngl(dABS(PLT)),0
ccccc         ELSE
ccccc            WRITE(IOTFIL,25) 29,.1,0
ccccc         ENDIF
      ENDIF
      PRINT *,' NUMBER OF RUNSETS READ=',NRUN-1
C
      write(43,25) 0
      WRITE(LOGFIL,680) NPPLT
      PRINT 680,NPPLT
      IF(NPPLT.LE.0) then
          PRINT *,' BE SURE THAT W(29) IS 10000,100000 OR 1000000'
      else
          print*,' Normal completion.'
          close (ndevtmp,status='DELETE')
      endif

680   FORMAT(1x,I4,' RAY POINTS PROCESSED')

ccccc      STOP
C
ccccc2500  PRINT *,' RAY SET FILE(RYSFIL) NOT PRESENT'
ccccc      STOP
C
C*2600  PRINT *,' FILE NAME FILE(PUNAM) NOT PRESENT'
C*      STOP
C
ccccc2700  PRINT *,' INPUT FILE DINP NOT PRESENT'
ccccc      STOP
C
ccccc2800  PRINT *,' INPUT FILE TAPE1 NOT PRESENT'
ccccc      STOP
C
      END
 
C -------------- End of EIGEN ------------------------------------------
 
C ----------------------------------------------------------------------
C			Function XLOR
C
C ----------------------------------------------------------------------
 
      LOGICAL FUNCTION XLOR(XL1,XL2)
      LOGICAL XL1,XL2
C
      XLOR=(XL1.AND..NOT.XL2).OR.(.NOT.XL1.AND.XL2)
      END
 
C ----------------------------------------------------------------------
C			Subroutine WRITW
C
c.....NOTE:  This routine is no longer used. Its function has been replaced
c            by Ron Brent's modifications. jah.
C ----------------------------------------------------------------------
 
      SUBROUTINE WRITW(NDEV,NDEVTMP,ID,W,WS,MAXW)
C
      implicit double precision (a-h,o-z)
 
      character id*80,line*80
      REAL*8 W(MAXW),WS(MAXW)
C*      REAL LINE(8)
C
C
      REWIND NDEVTMP
      DO 200 I=1,1000
           READ(NDEVTMP,190,END=10000) LINE
C*190   FORMAT(BZ,8A10)
C     SKIP THE ID
           WRITE(NDEV,190) LINE
200   CONTINUE
190   FORMAT(BZ,A80)
C
10000 REWIND NDEVTMP
C
      DO 100 I=1,MAXW
           IF(W(I).NE.WS(I)) WRITE(NDEV,20) I,sngl(W(I))
100        WS(I)=W(I)
20    FORMAT(I3,G14.7)
      END
C********************************************************************
      SUBROUTINE MARKER(MK)
C********************************************************************
C     SPECIFIES MARKER SYMBOL TO BE USED BY CURVE ROUTINE(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C     COMMON DECK "RAYDEV" INSERTED HERE
C     DEVICE ASSIGNED TO RAYTRC INPUT FILE
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT
 
      WRITE(NDEVTXT,100) 28,MK,0
100   FORMAT (3I2)
      END

