C file EIGEN.FOR
C ----------------------------------------------------------------------!-------
Contains EIGEN, RETRANSMIT
C ----------------------------------------------------------------------!-------
C 29 Oct 2009 - Added more documentation to doutp3.
C 14 Aug 2009 - Added W(459) as a phase to subtract from interpolated phases
C 14 Aug 2009 - Changed W(453) from rcvr range to rcvr central earth angle
C 13 Aug 2009 - Added interpolation of phase, Absorption, and Doppler shift
C 13 Aug 2009 - Restored output of time, phase, absorp, Doppler to doutp2
C  4 May 2009 - replaced REAL*8 with DOUBLE PRECISION
C 19 Mar 2008 - Added equivalence plt with w(81) to subroutine retransmit
C 13 Mar 2008 - Corrected length in common block /RAYDAT2/
C 13 Mar 2008 - Corrected length in common block /RAYCOM/
C 13 Mar 2008 - Corrected length in common block /HDR/
C 13 Mar 2008 - Corrected length in common block /PCONST/
C 29 Nov 2007 - Fixed to not print spreading for end ray
C 23 Nov 2007 - Added output to doutp4 and doutp5.
C 23 Nov 2007 - Added a space in last column of doutp3
C 20 Sep 2007 - Ensured that the argument of log10 is not zero
C 22 May 2007 - Format change in printing
C 18 May 2007 - Normalized spreading calculation by elevation angle step
C 16 May 2007 - Output ID to doutp3
C  9 May 2007 - Improved the format for particle speed.
C  8 May 2007 - Corrected multiple frequencies in runsets
C  7 May 2007 - Corrected spreading loss formula for vertical incidence
C  4 May 2007 - Allow multiple frequencies in runsets
C  3 May 2007 - Using runsets other than the first now works correctly
C  3 May 2007 - Added printing of frequency, source amplitude
C  3 May 2007 - Changed output format
C  2 May 2007 - Added calculation of particle speed
C  2 May 2007 - Changed output format
C  1 May 2007 - Added spreading loss.  This calculation is valid only when 
C  1 May 2007 -  there is azimuthal symmetry of the medium (the atmosphere) 
C  1 May 2007 -  about the transmitter.  For practical purposes, that usually 
C  1 May 2007 -  means that the atmospheric parameters depend only on height.  
C  1 May 2007 - Although the calculation of particle speed includes the factor 
C  1 May 2007 -  of the square root of density, it neglects the effect of the 
C  1 May 2007 -  polarization relations (that is, the differences in 
C  1 May 2007 -  perturbations in pressure, temperature, density, or wind 
C  1 May 2007 -  other than the square root of the density).
C 30 Apr 2007 - Format change
C 27 Apr 2007 - Changed what is printed
C 26 Apr 2007 - Corrected reading of date and time
C 24 Apr 2007 - Enlarge common block /RAYDAT/
C 24 Mar 2006 - Moved XLOR to S7
C 24 Mar 2006 - Moved WRITW, MARKER to S8
C 24 Mar 2006 - Changed COMMON block /RAYDEV/
C 13 September 2005 - Receiver lat. and long. calcuation enabled
C 25 August 2005 - Added lower-case variables to RAYDAT common block
C 25 August 2005 - Added Subroutine RETRANSMIT
C 23 August 2005 - Shortened DOUTP output file to 80 characters
C ----------------------------------------------------------------------!-------
C			Program EIGEN
C ----------------------------------------------------------------------!-------
C PROGRAM TO PLOT RANGE VERSUS ELEVATION ANGLE OF TRANSMISSION AND 
C RANGE VERSUS TRAVEL TIME AND TO CALCULATE EIGENRAYS
C ----------------------------------------------------------------------!-------
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
C ----------------------------------------------------------------------!-------
      PROGRAM EIGEN
C ----------------------------------------------------------------------!-------
	
	implicit double precision (a-h,o-z)
 
      DOUBLE PRECISION AZIMUTH,ELEVA,RANGE
      LOGICAL XLOR,WSIGAP,INGAP,HADEIS
C
C     COMMON DECK "HDR" INSERTED HERE
      CHARACTER*10 INITID*80,DAT,TOD
      COMMON/HDR/SEC,stsec                            !Rev 7/13/90 ! 13 Mar 2008
      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)
      DOUBLE PRECISION NRSLMI
C     character*132 rbtemp
      character*80 rbtemp                             ! 15 September 2005
      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),rcCEA) ,(W(454),SLPCUT)
     8   ,(W(456),HOPLMI) ,(W(457),NRSLMI) ,(W(458),ESET)
     9   ,(W(459),PHASE0)                                         ! 14 Aug 2009
      DOUBLE PRECISION 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
     1 ,MO2,MN2,MO,MN,MH2O,MO3                        ! 9 Mar 2007 ! 13 Mar 2008
      DOUBLE PRECISION MO2,MN2,MO,MN,MH2O,MO3                      ! 13 Mar 2008
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10
C     COMMON DECK "RAYDEV" INSERTED HERE
C     DEVICE ASSIGNED TO RAYTRC INPUT FILE
C     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT         ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
 
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
      integer hopmax                                               ! 24 Apr 2007
      DOUBLE PRECISION logdensT, logdensR                          ! 24 Apr 2007
      character*1 Tevent, event                                    ! 24 Apr 2007
      COMMON/RAYDAT/ nwndabs(333), xmtrht(333),                    ! 24 Apr 2007
     1 xmtrlat(333), xmtrlong(333), rcvrht(333), omega(333),       !  3 May 2007
     2 AZIMUTH(333), ELEVA(333), polarTR(333), polarTI(333),       ! 24 Apr 2007
     3 hopmax(333), logdensT(333),                                 ! 24 Apr 2007
     4 recvht(10000), RANGE(10000), dazimuth(10000),               ! 24 Apr 2007
     5 rcvrAZ(10000), rcvrEL(10000),                               ! 24 Apr 2007
     6 GRPTIM(10000), phase(10000), absorb(10000),                 ! 24 Apr 2007
     7 Doppler(10000), polarRR(10000), polarRI(10000),             ! 24 Apr 2007
     8 nhopin(10000), logdensR(10000), NELMTS                      ! 24 Apr 2007
      
      CHARACTER FID*3                                              ! 13 Mar 2008
      COMMON/RAYDAT2/ FID, event(10000)                            ! 24 Apr 2007

c      COMMON/RAYDAT/xmtrlat(333),xmtrlong(333),AZIMUTH(333),ELEVA(333)
c      COMMON/RAYDAT/recvht(10000),dazimuth(10000),RANGE(10000)
c      COMMON/RAYDAT/GRPTIM(10000),NELMTS,nhopin(10000)
C     !25 August 2005 - added lower-case variables to RAYDAT common block

      COMMON/GRPC/GRPN,GRPX
C
      INTEGER NCODES(3),PLTCOD
      CHARACTER INITRY*80
      CHARACTER*3 TESTID,PRVID
      PARAMETER (NHSZ=154,NWSSZ=400)
      DOUBLE PRECISION WH(NHSZ),WS(NWSSZ)
      DOUBLE PRECISION EL(333),RG(333),TT(333),GTIM(333)
     1 ,ptim(333),abs(333),Dop(333)                               ! 13 Aug 2009
C     COMMON DECK "CODES" INSERTED HERE
      COMMON/CODES/ICODRE,ICODRT,ICODEI
C
      CHARACTER TDID*30
      COMMON/CDATIM/TDID
        character cdate*11,ctime*10                              ! 30 Mar 2006

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
c320   FORMAT(T80,A11,2x,a10)                                   ! 30 Mar 2006
320   FORMAT(T81,A10,2x,a10)                                       ! 26 Apr 2007
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
C     IF(READW(0).NE.0.) GO TO 670                      ! 23 August 2005
      IF(READW(xxx).NE.0.d0) GO TO 670                  ! 23 August 2005

C
      WRITE(TESTID,'(A3)') ID(1:3)
      WRITE(TDID,330) ID(1:3),CDATE,CTIME
330   FORMAT(' ',A3,A11,1x,A10)                      ! 30 Mar 2006

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+
      ICODE29=ND2B(INT(RAYFNC))                                     ! 3 May 2007
c           IF(RAYSET.NE.0.0) NRYSET=NRYSET+1
         IF(RAYSET.NE.0.0.and.iand(icode29,4).ne.0) NRYSET=NRYSET+1 ! 3 May 2007
c      print * ,'rayfnc = ', rayfnc, ', nryset = ', nryset, ', 
c     1 icode29 = ',icode29
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)
c            print*,'eleva = ',eleva
c            print*,'elv = ',elv
            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
c            print*,'el =',el
            RG(NPTS)=RR
            GTIM(NPTS)=GRPTIM(LNHOP)
            PTIM(NPTS)=PHASE(LNHOP)                               ! 13 Aug 2009
            abs(NPTS)=absorb(LNHOP)                               ! 13 Aug 2009
            Dop(NPTS)=Doppler(LNHOP)                              ! 13 Aug 2009
            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 NOW 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 runset 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
c      print *, 'rbcount = ', rbcount
c      print *, rbtemp
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 following 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

C999  format(a132)
 999  format(a80)                                ! 23 August 2005
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,ptim,abs,dop,NPTS,NHOP,NPR,omega)!13-14 Aug2009
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
          call retransmit
          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
C ----------------------------------------------------------------------!-------
      END PROGRAM EIGEN
C ----------------------------------------------------------------------!-------
 
C -------------- End of EIGEN ------------------------------------------
C 13 September 2005 - Receiver lat. and long. calcuation enabled
C 25 August 2005 - Added lower-case variables to RAYDAT common block
C 25 August 2005 - Added subroutine RETRANSMIT
C ----------------------------------------------------------------------!-------
C                       Subroutine RETRANSMIT
C
C ----------------------------------------------------------------------!-------
      SUBROUTINE RETRANSMIT
C ----------------------------------------------------------------------!-------
C Allows the eigenray program use the ends of rays
C as virtual transmitters to calculate Fresnel zones. 
C The eigenray program produces a file named doutp2, that can be
C used as input to the ray tracing program.  
C ----------------------------------------------------------------------!-------

      implicit double precision (a-h,o-z)
 
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10
C     COMMON DECK "CRAYDAT" INSERTED HERE
      integer hopmax                                               ! 24 Apr 2007
      DOUBLE PRECISION logdensT, logdensR                                    ! 24 Apr 2007
      character*1 Tevent, event                                    ! 24 Apr 2007
      COMMON/RAYDAT/ nwndabs(333), xmtrht(333),                    ! 24 Apr 2007
     1 xmtrlat(333), xmtrlong(333), rcvrht(333), omega(333),       !  3 May 2007
     2 AZIMUTH(333), ELEVA(333), polarTR(333), polarTI(333),       ! 24 Apr 2007
     3 hopmax(333), logdensT(333),                                 ! 24 Apr 2007
     4 recvht(10000), RANGE(10000), dazimuth(10000),               ! 24 Apr 2007
     5 rcvrAZ(10000), rcvrEL(10000),                               ! 24 Apr 2007
     6 GRPTIM(10000), phase(10000), absorb(10000),                 ! 24 Apr 2007
     7 Doppler(10000), polarRR(10000), polarRI(10000),             ! 24 Apr 2007
     8 nhopin(10000), logdensR(10000), NELMTS                      ! 24 Apr 2007
      
      CHARACTER FID*3                                              ! 13 Mar 2008
      COMMON/RAYDAT2/ FID, event(10000)                            ! 24 Apr 2007

c      COMMON/RAYDAT/xmtrlat(333),xmtrlong(333),AZIMUTH(333),ELEVA(333)
c      COMMON/RAYDAT/recvht(10000),dazimuth(10000),RANGE(10000)
c      COMMON/RAYDAT/GRPTIM(10000),NELMTS,nhopin(10000)

      COMMON/RAYCOM/NRAYS,MXRYS,MXHOP,MXXHOP,NRSLIM,ICODE
      COMMON/RAYCOM/NHOPS(333),MNHOPS(333)                         ! 13 Mar 2008
 
      PARAMETER (NWARSZ=1000)
      CHARACTER*80 ID
      COMMON/IDENT/ID
      COMMON/WW/dum(10),MAXW,W(NWARSZ)
      EQUIVALENCE (EARTHR,W(1)),(HOP,W(22))
     1 ,(PLT,W(81))                                               ! 19 Mar 2008
     2 ,(ampT,w(478)),(rT,w(479))                                 !  1 May 2007
      DIMENSION RECLAT(10000),RECLONG(10000)
     1 ,spreadxmtr(333),spreadrcvr(10000),partspeed(10000)        !  2 May 2007
      DIMENSION THT(5000),TLAT(5000),TLONG(5000)
     1 , Tfreq(5000)                                              !  8 May 2007
     2 , GTIM(5000),PTIM(5000),abs(5000),Dop(5000)                ! 13 Aug 2009
      CHARACTER*80 rbtemp
 
      DATA MAXHOP/0/
 
c      PRINT*,'SUBROUTINE RETRANSMIT REACHED'

      INFILE=1
      REWIND INFILE
      CALL OPNREP(73,'doutp2')
      CALL OPNREP(74,'doutp3')                                      ! 3 May 2007
      CALL OPNREP(75,'doutp4')                                     ! 23 Nov 2007
      CALL OPNREP(76,'doutp5')                                     ! 23 Nov 2007

      write(unit=74,fmt=819) ampT, rT                               ! 3 May 2007
819   FORMAT(
     1 '   The first table in this file gives ',                    ! 1 May 2007
     2 'information about the ray at the transmitter.  ',/,'The ',  ! 1 May 2007
     3 'second table in this file gives information about the ',    ! 1 May 2007
     4 'ray at the receiver.  The ',/,'ray numbers in each ',       ! 1 May 2007
     5 'table correspond, to give information about the same ray.', ! 1 May 2007
     6 /,'"Xmtr spreading" is proportional to the cross ',          ! 1 May 2007
     7 'sectional area of a ray tube near the ',/,                  ! 1 May 2007
     8 'transmitter, converted to dB.  "Rcvr spreading" is ',       ! 1 May 2007
     9 'the same for the receiver end of ',/,'the ray.',/,          ! 1 May 2007
     1 '   This calculation of spreading loss is valid only when ', ! 1 May 2007
     2 'there is azimuthal symmetry ',/,'of the medium ',           ! 1 May 2007
     3 '(the atmosphere) about the transmitter.  For practical ',   ! 1 May 2007
     4 'purposes, that ',/,'usually means that the atmospheric '    ! 1 May 2007
     5 'parameters depend only on height.  Although the ',/,        ! 1 May 2007
     6 'calculation of particle speed includes the factor of ',     ! 1 May 2007
     7 'the square root of density, it '/,'neglects the effect of ',! 1 May 2007
     8 'the polarization relations (that is, the differences in ',  ! 1 May 2007
     9 /,'perturbations in pressure, temperature, density, or ',    ! 1 May 2007
     1 'wind other than the square root ',/,'of the density).',     ! 1 May 2007
     2 /,'   The first table is also output to a file doutp4, so ',! 29 Oct 2009
     3 'it can be read into a '/, 'spreadsheet program.  ',        ! 29 Oct 2009
     4 'The second table is also output to a file doutp5, so ',    ! 29 Oct 2009
     5 'it can be '/, 'read into a spreadsheet program.  ',        ! 29 Oct 2009
     6 //,'   Source amplitude = ' ,f9.3, ' m/sec at ',             ! 3 May 2007
     7 f9.3,' km from source',/)                                    ! 3 May 2007

      WRITE(unit=74,fmt=257) id                                    ! 16 May 2007

      WRITE(UNIT=74,FMT=820)                                        ! 3 May 2007
      WRITE(UNIT=75,FMT=820)                                       ! 23 Nov 2007
820   FORMAT( 13x,'Xmtr', 8x,'Xmtr', 4x,'frequency' 2x,'Azimuth'
     1, 4x,'Elevation', 7x, 'log', 5x, '  Xmtr   '                 !  1 May 2007
     2 ,/, 2x,'Ray', 9x,'Lat.', 7x,'Long.', 7x,'Hz', 7x,'deg'      !  3 May 2007
     3 , 9x,'deg', 8x, 'density', 3x, 'spreading',/,82x,'dB')      !  3 May 2007
      DO 710 i=1,MXRYS
         sinbetaT=dsin(rad*eleva(i))                               !  1 May 2007
c           WRITE(74,'(A,I3)') 'i = ',i
c           WRITE(74,'(A,I3)') 'hopmax(i+1) = ',hopmax(i+1)
         IF(i.ne.1.and.i.ne.MXRYS.and.sinbetaT.ne.0.d0             !  1 May 2007
     1      .and.hopmax(i+1).ne.0                                  ! 29 Nov 2007
     2       ) then                                                !  1 May 2007
            cosbeta1=dcos(rad*eleva(i-1))                          !  1 May 2007
            cosbeta2=dcos(rad*eleva(i+1))                          !  1 May 2007
            IF(rT.ne.0.d0                                          ! 20 Sep 2007
     1      .and.(dabs(cosbeta1)+dabs(cosbeta2)).ne.0.d0           ! 20 Sep 2007
     2      .and.(cosbeta1-cosbeta2).ne.0.d0                       ! 20 Sep 2007
     3      .and.(eleva(i-1)-eleva(i+1)).ne.0.d0)                  ! 20 Sep 2007
     4      spreadxmtr(i)=10.d0*DLOG10(Dabs(rT**2                  !  1 May 2007
     1      *(dabs(cosbeta1)+dabs(cosbeta2))                       !  7 May 2007
     2      *(cosbeta1-cosbeta2)/sinbetaT                          !  1 May 2007
     3      /(eleva(i-1)-eleva(i+1))                               ! 18 May 2007
     4      /(eleva(i-1)-eleva(i+1))                               ! 18 May 2007
     5                                   ))                        !  1 May 2007
         endif                                                     !  1 May 2007
         MAXSTEP=0
         DO 700 m=1,MXXHOP
            MAXSTEP=MAX(MAXSTEP,NHOPIN((i-1)*MXXHOP+m))
c           WRITE(74,'(A,I3)') 'MAXSTEP = ',MAXSTEP
700   CONTINUE
c           WRITE(74,'(A,I3)') 'maxstep = ',maxstep
         IF (MAXSTEP.EQ.0) GO TO 710                               !  1 May 2007
      freq=omega(i)/PIT2                                           !  3 May 2007
c      write(74,"(1x, 'mxrys = ',i5)") mxrys  
         IF(i.ne.1.and.i.ne.mxrys.and.sinbetaT.ne.0.d0             !  1 May 2007
     1      .and.hopmax(i+1).ne.0                                  ! 29 Nov 2007
     2       ) then                                                !  1 May 2007
c           WRITE(74,'(A,I3)') 'i = ',i
      WRITE(74,800) i,xmtrlat(i),xmtrlong(i),
     1 freq,                                                       !  4 May 2007
     2 azimuth(i),eleva(i)
     3 , logdensT(i), spreadxmtr(i)                                !  1 May 2007
      WRITE(75,800) i,xmtrlat(i),xmtrlong(i),                      ! 23 Nov 2007
     1 freq,                                                       
     2 azimuth(i),eleva(i)
     3 , logdensT(i), spreadxmtr(i)                                
       else                                                        !  1 May 2007
c           WRITE(74,'(A,I3)') 'i = ',i
      WRITE(74,800) i,xmtrlat(i),xmtrlong(i),
     1 freq,                                                       ! 4 May 2007
     2 azimuth(i),eleva(i)
     3 , logdensT(i)                                               ! 27 Apr 2007
      WRITE(75,800) i,xmtrlat(i),xmtrlong(i),                      ! 23 Nov 2007
     1 freq,                                                       
     2 azimuth(i),eleva(i)
     3 , logdensT(i)                                               
      endif                                                        !  1 May 2007
c           WRITE(74,'(A,I3)') 'i = ',i
710   CONTINUE
800   FORMAT(I5,7F12.3)                                            !  1 May 2007
      WRITE(74,'(1X)')
 
      WRITE(UNIT=74,FMT=830)                                        ! 3 May 2007
      WRITE(UNIT=76,FMT=830)                                       ! 23 Nov 2007
c830   FORMAT('Ray',1x,'Hop',2x,'Ht.',5x,'DAz.',6x,'Range',9x'Time',
c     1 10x,'RLat.',9x,'RLong.')

830   FORMAT( 8x,'Height', 3x,'DAz.', 6x,'Range', 2x'Time',
     1 5x,'elev.', 2x,'absorp.', 2x,'log', 5x,'Rcvr'                ! 1 May 2007
     2 , 3x, 'particle'                                             ! 2 May 2007
     3 ,/, 'Ray', 1x,'Hop', 2x,'km', 6x,'deg', 7x,'km', 5x'Sec',    ! 3 May 2007
     4 6x,'deg', 6x,'dB', 3x,'density', 1x,'spreading'              ! 3 May 2007
     5 , 1x, 'speed',/,67x,'dB', 5x, 'm/sec')                       ! 3 May 2007
      MAXI=0
      n=0

c           WRITE(74,'(A,I3)') 'MXRYS = ',MXRYS
      DO 720 k=1,MXRYS

         MAXHOP=0
c           WRITE(74,'(A,I3)') 'MXXHOP = ',MXXHOP
         DO 730 j=1,MXXHOP
            MAXHOP=MAX(MAXHOP,NHOPIN((k-1)*MXXHOP+j))
c           if(maxhop.ne.0) WRITE(74,'(A,I3)') 'MAXHOP = ',MAXHOP
c           if (NHOPIN((k-1)*MXXHOP+j).ne.0) 
c     1     WRITE(74,'(A,I3)') 'NHOPIN = ',NHOPIN((k-1)*MXXHOP+j)
c           WRITE(74,'(A,I3)') 'NHOPIN = ',NHOPIN((k-1)*MXXHOP+j)
c           WRITE(74,'(A,I3)') 'i = ',i
730   CONTINUE
         DO 720 NN=1,MAXHOP
         i=(k-1)*MXXHOP+NN
         i1=(k-2)*MXXHOP+NN                                        !  1 May 2007
         i2=(k)*MXXHOP+NN                                          !  1 May 2007
         THETAT=PID2-RAD*xmtrlat(k)
c      print *, 'range (',i,') = ', range(i)
c      print *, 'earthr = ', earthr
         D=RANGE(i)/EARTHR
c      print *, 'D = ', D
         SINT=DSIN(THETAT)
         SIND=DSIN(D)
         COST=DCOS(THETAT)
         COSD=DCOS(D)
         CALL BETTERCEA(THETAT,D,THETAR,COST,COSD,COSR,SINT,SIND,SINR,
     1                 (AZIMUTH(k)-DAZIMUTH(i))*RAD)
         SINA=DSIN(RAD*(AZIMUTH(k)-DAZIMUTH(i)))
         RECLAT(i)=DEGS*(PID2-THETAR)
         RECLONG(i)=XMTRLONG(k)+DEGS*(SIN(SIND*SINA/SINR))
         N=N+1
         THT(n)=RECVHT(i)
         TLAT(n)=RECLAT(i)
         TLONG(n)=RECLONG(i)
         GTIM(n)=GRPTIM(i)                                        ! 13 Aug 2009
         PTIM(n)=PHASE(i)                                         ! 13 Aug 2009
         abs(n)=absorb(i)                                         ! 13 Aug 2009
         Dop(n)=Doppler(i)                                        ! 13 Aug 2009
         Tfreq(n)=omega(k)/PIT2                                    ! 8 May 2007
         MAXI=n
c      write(74,"(1x, 'mxrys = ',i5)") mxrys
c           WRITE(74,'(A,I3)') 'i1 = ',i1
c           WRITE(74,'(A,I3)') 'i2 = ',i2
c           WRITE(74,'(A,G)') 'range(i2) = ',range(i2)
      if(k.ne.1.and.k.ne.MXRYS
     1   .and.dsin(rad*rcvrEL(i)).ne.0.d0                          ! 20 Sep 2007
     2   .and.(range(i1)+range(i2)).ne.0.d0                        ! 20 Sep 2007
     3   .and.(range(i1)*dcos(rad*dazimuth(i1))                    ! 20 Sep 2007
     4   -range(i2)*dcos(rad*dazimuth(i2)) )  .ne.0.d0             ! 20 Sep 2007
     5   .and.(eleva(k-1)-eleva(k+1)).ne.0.d0                      ! 20 Sep 2007
     6   .and.nhopin(i2).ne.0                                      ! 29 Nov 2007
     7   ) then                                                    !  1 May 2007
         spreadrcvr(i)=10.d0*DLOG10(dabs(dsin(rad*rcvrEL(i))       !  1 May 2007
     1   *(range(i1)+range(i2))*                                   !  1 May 2007
     2   (range(i1)*dcos(rad*dazimuth(i1))                         !  7 May 2007
     3   -range(i2)*dcos(rad*dazimuth(i2)) )                       !  7 May 2007
     4      /(eleva(k-1)-eleva(k+1))                               ! 18 May 2007
     5      /(eleva(k-1)-eleva(k+1))                               ! 18 May 2007
     6                                       ))                    !  7 May 2007
c           WRITE(74,'(A,I3)') 'then k = ',k
         partspeed(i)=ampT*10.d0**(-(absorb(i)+spreadrcvr(i)
     1   -spreadxmtr(k)+10.d0*logdensR(i)-10.d0*logdensT(k))/20.d0) ! 2 May 2007
      WRITE(74,810) k,nhopin(i),recvht(i),dazimuth(i),range(i),
     1  grptim(i), rcvrEL(i), absorb(i), logdensR(i), spreadrcvr(i) ! 1 May 2007
     2  ,partspeed(i)                                               ! 2 May 2007
      WRITE(76,810) k,nhopin(i),recvht(i),dazimuth(i),range(i),    ! 23 Nov 2007
     1  grptim(i), rcvrEL(i), absorb(i), logdensR(i), spreadrcvr(i) 
     2  ,partspeed(i)                                               
      else                                                         !  1 May 2007
c           WRITE(74,'(A,I3)') 'else k = ',k
      WRITE(74,810) k,nhopin(i),recvht(i),dazimuth(i),range(i),
     1  grptim(i), rcvrEL(i), absorb(i), logdensR(i)               ! 27 Apr 2007
c     1              grptim(i),reclat(i),reclong(i)
      WRITE(76,810) k,nhopin(i),recvht(i),dazimuth(i),range(i),    ! 23 Nov 2007
     1  grptim(i), rcvrEL(i), absorb(i), logdensR(i)               
c     1              grptim(i),reclat(i),reclong(i)
      endif                                                        !  1 May 2007
c           WRITE(74,'(A,I3)') 'k = ',k
c      write(74,"(1x, 'mxrys = ',i5)") mxrys
720      CONTINUE
c810   FORMAT(2I3,2G9.2,6G11.4)                                    ! 30 Apr 2007
c810   FORMAT(2I3,2F8.2,F9.2,2F8.2,F9.2,F7.3, F8.2,G10.3)          !  2 May 2007
c810   FORMAT(2I3,2F8.2,F10.3,2F8.2,F9.2,F7.3, F8.2,G10.3)         ! 22 May 2007
810   FORMAT(2I3,2F8.2,F10.3,2F8.2,F9.2,F7.3, F8.2,1x,G10.3)       ! 23 Nov 2007

Crb   RBCOUNT is a variable that keeps track of which runset 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.
c      PRINT*,'MAXI = ',MAXI

      rbcount=0
      do 255 i=1,5000
           read(infile,999,err=256) rbtemp
           if (rbtemp(25:27).eq.'***') rbcount=rbcount+1
c      print *, 'rbcount = ', rbcount
c      print *, rbtemp
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 following 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(73,257) id
           ELSE
                WRITE(73,999) rbtemp
           ENDIF
255   CONTINUE
 
999   FORMAT(A80)
Crb   This is where various flags are written to DOUTP2

256   WRITE(73,25) 81,sngl(plt)
      WRITE(73,25) 29,100.
25    FORMAT(I3,G14.7)
26    FORMAT(I3,T25,'***')                                  ! 13 Aug 2009
257   FORMAT (A80)
258   FORMAT (A5)

      WRITE(73,89)                                                  ! 4 May 2007
      DO 900 i=1,MAXI-1
      WRITE(73,85) THT(i)
      WRITE(73,86) TLAT(i)
      WRITE(73,87) TLONG(i)
      WRITE(73,88) Tfreq(i)                                        ! 8 May 2007
      WRITE(73,90) GTIM(i)                                        ! 13 Aug 2009
      WRITE(73,91) PTIM(i)                                        ! 13 Aug 2009
      WRITE(73,92) abs(i)                                         ! 13 Aug 2009
      WRITE(73,93) Dop(i)                                         ! 13 Aug 2009
      WRITE(73,26) 0                                              ! 13 Aug 2009
C      WRITE(73,258) ID(1:5)
      WRITE(73,257) ID                                            ! 13 Aug 2009
900   CONTINUE
85    FORMAT('  3',T4,G14.7,T25,'TRANSMITTER HEIGHT ABOUT MSL')
86    FORMAT('  4',T4,G14.7,T18,'AN DG',T25,'N. TRANSMITTER LATITUDE')
87    FORMAT('  5',T4,G14.7,T18,'AN DG',T25,'E. TRANSMITTER LONGITUDE')
88    FORMAT('  7',T4,G14.7,T18,'FQ HZ',T25,'INITIAL FREQUENCY, HZ')! 4 May 2007
89    FORMAT('  9',T9,'0.0',T18,'FQ HZ',T25,'DO NOT STEP FREQUENCY')! 4 May 2007
90    FORMAT(' 56',T4,G14.7,T18,'     ',T25,'GROUP TIME, SEC')     ! 13 Aug 2009
91    FORMAT(' 57',T4,G14.7,T18,'     ',T25,'PHASE TIME, SEC')     ! 13 Aug 2009
92    FORMAT(' 58',T4,G14.7,T18,'     ',T25,'ABSORPTION, dB')     ! 13 Aug 2009
93    FORMAT(' 59',T4,G14.7,T18,'     ',T25,'Doppler shift, Hz')  ! 13 Aug 2009
      WRITE(73,85) THT(MAXI)
      WRITE(73,86) TLAT(MAXI)
      WRITE(73,87) TLONG(MAXI)
      WRITE(73,90) GTIM(MAXi)                                     ! 13 Aug 2009
      WRITE(73,91) PTIM(MAXi)                                     ! 13 Aug 2009
      WRITE(73,92) abs(MAXi)                                      ! 13 Aug 2009
      WRITE(73,93) Dop(MAXi)                                      ! 13 Aug 2009
      WRITE(73,26) 0                                              ! 13 Aug 2009

      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE RETRANSMIT
C ----------------------------------------------------------------------!-------

