C**********************************************************************

      PROGRAM RAYTRC

C**********************************************************************

C     MAIN PROGRAM FOR THE RAY TRACING PACKAGE.

C          SETS THE INITIAL CONDITIONS FOR EACH RAY AND CALLS TRACE

C

C     COMMON DECK "RAYDEV" INSERTED HERE

C     DEVICE ASSIGNED TO RAYTRC INPUT FILE

      implicit real*8 (a-h,o-z)

      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVGRP,NDEVBIN

      LOGICAL WCHANGE,FIRST

      REAL*8 WS(400),kay2,kay2i,kr,kth,kph

     1,maxstp,maxerr,intyp,llat,llon

      real*8 modg

C     COMMON DECK "FILEC" INSERTED HERE

      COMMON /FILEC/NPLTDP

C     COMMON DECK "CERR" INSERTED HERE

      COMMON/ERR/NERG,NERR,NERT,NERP

C     COMMON DECK "GG" INSERTED HERE

      COMMON/GG/MODG(4)

      COMMON/GG/G,PGR,PGRR,PGRTH,PGRPH

      COMMON/GG/PGTH,PGPH,PGTHTH,PGPHPH,PGTHPH,GSELECT,GTIME

C     COMMON DECK "HDR" INSERTED HERE

      CHARACTER*10 INITID*80,DAT,TOD

      COMMON/HDR/SEC,stsec                             !Rev 7/13/90

      COMMON/HDRC/INITID,DAT,TOD

C     COMMON DECK "CONST" INSERTED HERE

      COMMON/PCONST/CREF,RGAS,GAMMA

      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10

C     COMMON DECK "FLAG" INSERTED HERE

      LOGICAL NEWWR,NEWWP,NEWTRC,PENET

      COMMON /FLG/ NTYP,NEWWR,NEWWP,NEWTRC,PENET,LINES,IHOP,HPUNCH

      COMMON/FLGP/NSET

C     COMMON DECK "RINPLEX" INSERTED HERE

      COMPLEX*16 PNP,POLAR,LPOLAR

      LOGICAL SPACE

      CHARACTER DISPM*6,modrin*64

      COMMON/RINPL/DISPM

      COMMON /RIN/ MODRIN,RAYNAME(2,3),TYPE(3),SPACE

      COMMON/RIN/OMEGMIN,OMEGMAX,KAY2,KAY2I

      COMMON/RIN/PNP(10),POLAR,LPOLAR,SGN

C     COMMON DECK "RK" INSERTED HERE

C     DEFINE SIZE REQUIRED FOR RAY STATE SAVE ARRAY

      PARAMETER (LRKAMS=87+2*100,NXRKMS=12+LRKAMS,MXEQPT=21)

      PARAMETER (NRKSAV=NXRKMS+MXEQPT-1)

      COMMON /RK/ NEQS,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART

C     COMMON DECK "RKAM" INSERTED HERE

      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)

C     COMMON DECK "WWR" INSERTED HERE

      COMMON/WW/DUM(10),MAXW,W(1000)

      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))

      EQUIVALENCE (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))

C

      CHARACTER TMP80*80

C

      DATA FIRST/.TRUE./

C

C   SET CONSTANTS , INITITALIZE FILES  AND MAKE SYSTEM CALLS

c%      CALL SYSINI

      CALL CONBLK

      CALL SETDSP

      CALL STDINI

      call syssec(stsec)               !Rev 7/13/90  get start time

C     INITIALIZE LINE COUNTS FOR CURRENT AND LAST HEADER

      LINES=0

      LHDRPG=0

C

      REWIND NRYIND

C  INITIALIZE RAYSET FILE

5     READ(NRYIND,'(A)',END=8) TMP80

      WRITE(9,'(A)') TMP80

      GO TO 5

C

8     WRITE(9,'(A)') '****SOFT EOF***+++++'

      REWIND NRYIND

C     READ IN USERS NAME AND TELEPHONE EXTENSION FOR IDENTIFYING

C     MICROFILM PLOTS.

C

C     REPOSITION DINP AFTER IDENTIFIER

      READ(NRYIND,'(A)',END=6000) TMP80

C

C********* READ W ARRAY AND PRINT NON-ZERO VALUES

      IRUN=0

      NSET=0

C

10    xyz= READW(xx)

      NSET=NSET+1

      ICODE=ND2B(IDINT(RAYFNC))

C PROCESS RAYPATH CALCULATIONS ONLY IF W(29) IS BEING USED

c*      IF((RAYFNC.NE.0.).AND.(AND(ICODE,4).EQ.0)) GO TO 10

C*  For VAX A Newhall  (AND to IAND)

      IF((RAYFNC.NE.0.d0).AND.(IAND(ICODE,4).EQ.0)) GO TO 10

C

C********* LET ROUTINES PRINTR AND RAYPLT KNOW THERE IS A NEW W ARRAY

      NEWWP=.TRUE.

      NEWWR=.TRUE.

      IRUN=IRUN+1

C

C     SET BINARY RAY FILE UNIT TO ZERO IF NO OUTPUT WANTED

      NPLTDP=0

      IF(BINRAY.NE.0.d0) NPLTDP=NDEVBIN

C

C******** INITIALIZE THE MODELS VIA 'DISPER'

      CALL IDISPER

C

      IF(FIRST.AND.PLT.NE.0.d0) CALL DDINIT(8,INITID)

      IF(PLT.NE.0.d0) FIRST=.FALSE.

C

      OW=0.d0

      BETA=0.d0

      AZ1=0.d0

C

      CALL HEADER1

C

C********* PRINT OUT THE CONTENTS OF THE 'W' ARRAY

C********* DETERMINE TRANSMITTER LOCATION IN COMPUTATIONAL COORDINATE

C********* SYSTEM (GEOMAGNETIC COORDINATES IF DIPOLE FIELD IS USED)

      SP=DSIN (PLAT)

      CP=DSIN (PID2-PLAT)

      SDPH=DSIN (TLON-PLON)

      CDPH=DSIN (PID2-(TLON-PLON))

      SL=DSIN (TLAT)

      CL=DSIN (PID2-TLAT)

      ALPHA=DATAN2(-SDPH*CP,-CDPH*CP*SL+SP*CL)

      TH0=DACOS (CDPH*CP*CL+SP*SL)

      PH0=DATAN2(SDPH*CL,CDPH*SP*CL-CP*SL)

C

      R=EARTHR

      TH=TH0

      PH=PH0

      CALL TOPOG

C

C     OBTAIN ABSOLUTE HEIGHT OF THE TRANSMITTER.

C     IF IT WAS SPECIFIED AS RELATIVE TO THE TERRAIN, THEN REMOVE

C**   THE FLAG VALUE 10**40 WHICH WAS ADDED AT INPUT.

C     THE FLAG VALUE 10**10 WHICH WAS ADDED AT INPUT.

C

      TMP=XMTRH

C**   IF(XMTRH .EQ. 1.E-40) XMTRH=0.0

      IF(XMTRH .EQ. 1.d-10) XMTRH=0.0

C**   IF(DABS(XMTRH) .GE. 1.E20 ) XMTRH=XMTRH*1.E-40

      IF(DABS(XMTRH) .GE. 1.d10 ) XMTRH=XMTRH*1.d-10

      IF(TMP.NE.XMTRH) XMTRH=XMTRH-G/PGR

C

C     CHECK THAT TRANSMITTER IS ABOVE TERRAIN.

C

      IF(-G/PGR .LE. XMTRH) GO TO 655

C

      PRINT *,'transmitter below terrain; run #',irun

      WRITE(3,640) IRUN

640   FORMAT('0****** TRANSMITTER BELOW TERRAIN. RUN ',I3

     1 ,' TERMINATED.'/'SEE W-ARRAY PRINTOUT. INPUT CONTINUES.'//)

      xyy= setovr(xx)

      zyx= printw(xx)

      GO TO 10

C

C     PREVENT RERUNNING SAME MODEL CASES

655   IF(IRUN.EQ.1) GO TO 12

      IF(RUNSUP.NE.0.d0 .AND. .NOT.WCHANGE(WS(1),W(1))) THEN

         PRINT *,'SUPPRESSING REDUNDANT RUNSET #',NSET

         GO TO 10

      ENDIF

C

12    PRINT *,'PROCESSING FOR RUN SET NUMBER',NSET

C

      CALL RMOVE(WS,W,400)

C

      xyz=setovr (xx)

      zyx= printw(xx)

      CALL IPRINTR

C******** INITIALIZE PRINT CONTROL PARAMTERS

      LINSPP=PAGLN

      LNPHPG=LINSPP/2

      IF(LNPHPG.LT.40) LNPHPG=LINSPP

C

C********* LOOP ON FREQUENCY, AZIMUTH ANGLE, AND ELEVATION ANGLE

      NFREQ=1

      IF (FSTEP.NE.0.d0) NFREQ=(FEND-FBEG)/FSTEP+1.5d0

      NAZ=1

      IF (AZSTEP.NE.0.d0) NAZ=(AZEND-AZBEG)/AZSTEP+1.5d0

      NBETA=1

      IF (ELSTEP.NE.0.d0) NBETA=(ELEND-ELBEG)/ELSTEP+1.5d0

      DO 50 NF=1,NFREQ

      OW=FBEG+(NF-1)*FSTEP

      DO 45 J=1,NAZ

      AZ1=AZBEG+(J-1)*AZSTEP

      AZA=AZ1*DEGS

      GAMMA1=PI-AZ1+ALPHA

      SGAMMA=DSIN (GAMMA1)

      CGAMMA=DSIN (PID2-GAMMA1)

      DO 40 I=1,NBETA

      BETA=ELBEG+(I-1)*ELSTEP

      EL=BETA*DEGS

      print *,'ELEVATION ANGLE=',sngl(EL),'  DEG'  !Rev 7/16/90

      CBETA=DSIN (PID2-BETA)

      R=EARTHR+XMTRH

      TH=TH0

      PH=PH0

      KR=DSIN (BETA)

      KTH=CBETA*CGAMMA

      KPH=CBETA*SGAMMA

      TPULSE=0.d0

      RSTART=1.d0

C********* THE FOLLOWING LINE NEEDED FOR RAY TRACING IN COMPLEX SPACE

      SGN=1.d0

C********* CALL MODELS

      CALL DISPER

C...............LINPG= # of lines printed so far on this page

C...............LHDRPG= # of pages times lines per page

C...............LINES= total # of lines printed, assuming full pages

C...................(except for this page)

C...............LINSPP= # of lines per page

C...............LNPHPG= # of lines per half page

      lhdrpg=lines/linspp*linspp ! added to get correct line count

      LINPG=LINES-LHDRPG

      IF(I.EQ.1  .OR.  LINPG.GE.LINSPP-20) THEN

C

C     PUT OUT SUBHEADERS FROM MEDIA AND PRINTR ROUTINES

         CALL HEADER2

         CALL PRNHD1(' ')

C

C     COMPUTE LINE COUNT OF THIS HEADER

         LHDRPG=LINES/LINSPP*LINSPP

      ELSEIF(LINPG.GE.LNPHPG-10 .AND. LINPG.LE.LNPHPG) THEN

C     PUT OUT PAGE FEED WITH SUBHEADER IF AT HALF PAGE

         CALL PRNHD2('1')

      ELSE

C     PUT OUT SUBHEADER

         CALL PRNHD2(' ')

      ENDIF

C

      IF (KAY2.GT.0.d0) GO TO 30

      WRITE(3,2900) OMEGMIN,OMEGMAX

 2900 FORMAT (58H0TRANSMITTER IN EVANESCENT REGION, TRANSMISSION IMPOSSI

     1BLE/20H0MINIMUM FREQUENCY =,E17.10,20H MAXIMUM FREQUENCY =,E17.10)

      GO TO 44

30    r999=renorm(kr,kay2,3)

c30    CALL RENORM(KR,KAY2,3)

      if(neqs-6 .gt. 0) then

        CALL CLEAR(RKVARS,NEQS-6)

      endif

      CALL TOPOG

      IF(NERG.GT.0) RKVARS(NERG)=G

      IF(NERR.GT.0) RKVARS(NERR)=PGR

      IF(NERT.GT.0) RKVARS(NERT)=PGTH

      IF(NERP.GT.0) RKVARS(NERP)=PGPH

C

C     CALCULATE ONE RAY PATH

      CALL TRACE

      OSEC=SEC

      CALL SYSSEC(SEC)

      DIFF=SEC-OSEC

C

C     ADD TO LINES COUNT FOR ELAPSED TIME REPORT

      LINES=LINES+2

      WRITE(3,3500)SNGL(DIFF)

3500  FORMAT (/T93,'THIS RAY CALCULATION TOOK ',F8.3,' SEC')

C

      IF (PENET.AND.ONLY.NE.0.d0.AND.IHOP.EQ.1) GO TO 44

   40 CONTINUE

44    IF(PLT.GT.0.d0.AND.(NAZ.LE.1.OR.NBETA.GT.1)) CALL ENDPLT

   45 CONTINUE

      IF(PLT.GT.0.d0.AND.NAZ.GT.1.AND.NBETA.LE.1)CALL ENDPLT

      IF(PENET.AND.ONLY.NE.0.d0.AND.IHOP.EQ.1.AND.NAZ.EQ.1.AND.

     1 NBETA.EQ.1) GO TO 55

   50 CONTINUE

  55  IF (RAYSET.NE.0.d0) WRITE(9,5000)

 5000 FORMAT (78X,1H-)

      GO TO 10

6000  PRINT *,'DINP EMPTY OR NOT FOUND'

      STOP

      END

 

 

 

 

 

C**********************************************************************

      SUBROUTINE DFSYS

C**********************************************************************

 

C     ALL ROUTINES WHICH RELY ON SYSTEM BASED FUNCTIONS, SUCH AS TIME

C     OF DAY OR MAY BENEFIT FROM SYSTEM AVAILABLE ROUTINES FOR BETTER

C     PERFORMANCE, SUCH AS 'RMOVE' HAVE BEEN COLLECTED IN THIS ROUTINE.

 

      implicit double precision (A-H,O-Z)

      CHARACTER DAT*(*),DATX*10

      CHARACTER TIM*(*),TIMX*10

      character pc*1,typ*20

      common/crmach/rmach(9)

      integer ix(*),iy(*)

      real*8 x(*),y(*)



C**********************************************************************

      ENTRY RMOVE(X,Y,N)

C**********************************************************************

C     MOVE 'N' COMPONENTS OF ARRAY 'Y' TO ARRAY 'X'.

C     THIS LOGIC WILL HANDLE OVERLAP CASES ONLY FOR X(M)=X(M+N1)

C     WITH N1>0, M=1,N.

c     replaces CYBER MOVLEV routine

      DO 10 I=1,N

10    X(I)=Y(I)

      RETURN

C**********************************************************************

      ENTRY IMOVE(IX,IY,N)

C**********************************************************************

C     SPECIAL ENTRY FOR INTEGER ARRAYS FOR MACHINES HAVING DIFFERENT

C     WORD SIZES FOR INTEGERS THAN REALS.

C     replaces CYBER BUILT-IN ROUTINE:

C     CALL MOVLEV(Y,X,N)

       DO 20 I=1,N

 20    IX(I)=IY(I)

      RETURN

 

 

C**********************************************************************

      ENTRY SYSSEC(SECS)

C**********************************************************************

c      ms fortran call to get system time for run times

      call gettim(ihr,imin,isec,i100th)

      xhr=float(ihr)

      xmin=float(imin)

      secs=float(isec)

      x100th=float(i100th)/100.d0

      secs=secs+x100th

      secs=secs+60.d0*xmin

      secs=secs+3600.d0*xhr

      RETURN

 

C**********************************************************************

      ENTRY SYSTIM(TIM)

C**********************************************************************

c     ms fortran system time call for time labels

      call gettim(ihr,imin,iday,i100th)

      write (tim,'(2x,i2,1hh,i2,1hm,2x)') ihr,imin

      RETURN

 

 

C**********************************************************************

      ENTRY SYSDAT(DAT)

C**********************************************************************

c     ms fortran call to get system date for labels

      call getdat(iyr,imon,iday)

      write (dat,'(i2,1h/,i2,1h/,i4)') imon,iday,iyr

      RETURN



C**********************************************************************

      ENTRY MORTEM

C**********************************************************************

c     can add system post-mortem dump here

      STOP 'POST MORTEM'

c

c

C**********************************************************************

      entry pcodes (pc,typ)

C**********************************************************************

c     insert printer codes for form feed (pc), typography, etc.

c     the ones defined here are for the HP LaserJet 2

c     using line printer font with 8 lpi. 

c     LaserJet 1 uses A or B cartridge and prints in landscape mode.

c     pc=^L will work for most printers

      pc=16#0C

      write (typ,'(a1,a4,a1,a8)')16#1B,'&l8D',16#1B,'(s16.66H'

      return

c

c

C**********************************************************************

      entry mach

C**********************************************************************

c     replaces the DFCNST subroutine,

c     since these are the only machine constants used by HARPO

      rmach(3)=2.d0**(-52)  ! smallest relative spacing

      rmach(6)=2147483647.d0  ! largest 4-byte integer

      return

c

      END

 

 

 

 

C**********************************************************************

      SUBROUTINE READW1(AB,NWOK,MSET,MX,NTBL,ITBL,FRMTBL,GP)

C**********************************************************************

C     HANDLES INPUT OF TABULAR DATA REQUIRED BY SOME MODELS.

C     COMMON DECK "RAYDEV" INSERTED HERE

C     DEVICE ASSIGNED TO RAYTRC INPUT FILE

 

      implicit double precision (a-h,o-z)

 

      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVGRP,NDEVBIN

C     COMMON DECK "CUCON" INSERTED HERE

      COMMON/UCONV/CNVV(5,4)

      CHARACTER PCV*3,CNVC*2

      COMMON/UCONC/PCV(5),CNVC(5,4)

C     COMMON DECK "CONST" INSERTED HERE

      COMMON/PCONST/CREF,RGAS,GAMMA

      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10

C     COMMON DECK "WWR" INSERTED HERE

      COMMON/WW/DUM(10),MAXW,W(1000)

      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))

      LOGICAL NWOK,AB,ucon,uucon

C

      CHARACTER LINE*80,ALPHA*3,FBUF*3

      PARAMETER (NVAR=6,ALPHA='  A',NACDE=-1,NFIELD=13)

      INTEGER IPV(NVAR),NTBL(10),ITBL(10),FRMTBL(10)

      CHARACTER*8 PU(NVAR)

      REAL*8 GP(*),V(NVAR),CONV(NVAR)

C

C     THE DATA STRUCTURE EXPECTED BY THIS PROGRAM FOR A DATA BLOCK IS

C     SPECIFIED FOR A GIVEN MODEL BY THE THREE ARRAYS NTBL,ITBL,FRMTBL.

C     THE TYPE OF DATA EXPECTED FOR A GIVEN FORMAT GROUP IS TAKEN FROM

C     ARRAY FRMTBL. IT MUST HAVE AN ENTRY FOR A GIVEN TYPE OF FORMAT

C     WHICH IS EITHER ZERO, IF THAT FORMAT IS NOT ALLOWED OR AN ENTRY

C     WHOSE VALUE IS EQUAL TO THAT FORMAT TYPE. CURRENTLY ONLY TWO

C     FORMAT TYPES ARE ALLOWED, 1 FOR ALPHA AND 2 FOR NUMERIC.

C     TO ALLOW MORE FLEXIBLITY IN THE INPUT FILE THE NUMERIC FORMAT

C     TYPE IN TURN ALLOWS FOR 3 DIFFERENT DATA FORMATS, SEE BELOW.

C

C     OFFSETS INTO THE GENERAL ARRAY 'GP' FOR EACH FORMAT ARE GIVEN

C     IN THE ARRAY 'NTBL'. THE SPACING BETWEEN DATA VALUES IN EACH LINE

C     READ IS GIVEN IN ARRAY 'ITBL'. THIS SCHEME DOES NOT ALLOW FOR X,Y

C     ,Z CYCLES BUT RATHER SEPARATE ARRAYS OF X VALUES, THEN Y VALUES,

C     AND THEN Z VALUES(ETC. UP TO 'NVAR' VARIABLES).

C

      NWOK=.TRUE.

C

C     INITIALIZE ANY ALPHANUMERIC ARRAYS TO BLANKS

      dummy=1h 

      IF(FRMTBL(1).EQ.1) CALL SET2(GP(NTBL(1)),dummy,NTBL(2)-NTBL(1))

C

C     BEGIN MULTI-FORMAT LOOP

5     READ(NRYIND,'(A,A)',END=200) FBUF,LINE

      IF(NDEVTMP.GT.0) WRITE(NDEVTMP,'(A,A)') FBUF,LINE

      IF(FBUF.EQ.ALPHA) THEN

         NFRM=NACDE

      ELSE

         READ(FBUF,10) NFRM

      ENDIF

10    FORMAT(BZ,I3,A)

      IF(NFRM .EQ. 0) RETURN

C

C     FORMAT #1 IS ASSIGNED TO ALPHANUMERIC INPUT AND LARGER NUMBERS

C     FOR NUMERIC INPUT.

C

C     CHECK NOW FOR NUMERICAL FORMATS.

C     THESE ARE PROVIDED FOR EASE OF INPUT OF TABULATED DATA

C     AND ARE PARTLY 'TRANSPARENT' TO THE MODEL INVOLVED. I.E.

C     THE USER HAS THE OPTION OF SPECIFING THE NUMBER OF DATA COLUMNS.

C     INPUT FORMAT NUMBER SPECIFIES THE NUMBER OF COLUMNS OF INPUT DATA

C     I.E. FOR VALUE 1 A SINGLE DATA COLUMN IS EXPECTED, FOR A VALUE

C     2 TWO DATA COLUMNS ARE EXPECTED, ETC.

C

C     DETERMINE FORMAT NUMBER

      IF(NFRM.EQ.NACDE) THEN

         IT=1

C     ALPHANUMERIC FORMATS MUST BE FIRST IN LIST

C     SINCE NO OFFSET SPECIFIER IS ALLOWED IN A FULL ALPHA LINE

         IF(FRMTBL(1).EQ.1) GO TO 100

      ELSE

         READ(LINE,'(I2)') IT

         IT=MAX0(IT,2)

         IF(IT.LE.MX) GO TO 100

      ENDIF

C

15    AB=.TRUE.

      PRINT 20, MSET,NFRM,MX

20    FORMAT(' FOR SET',I5,' FORMAT NUMBER',I5,' EXCEEDS LIMIT OF',I5)

C      ERROR IN INPUT DATA , STOP HERE

      STOP

C

C     THE 'NTBL' TABLE PROVIDES A LIST OF OFFSETS TO VARIABLE GROUPS IN

C     MODEL ARRAYS, 'ITBL' GIVES ELEMENT SEPARATION TO ALLOW FOR 2-

C     DIMENSIONAL ARRAYS AND 'FRMTBL' SPECIFIES THE INPUT FORMAT.

C

100   N1=NTBL(IT)

      N3=ITBL(IT)

      N2=NTBL(IT+1)-N3

C

      IF(NFRM.NE.NACDE) GO TO 110

C

      READ(LINE,1010)((GP(I+J-1),I=N1,N2,N3),J=1,N3)

1010  FORMAT(10A8)

      GO TO 5

C

110   IF(NFRM.LT.1 .OR. NFRM.GT.NVAR) GO TO 200

C     HANDLE NUMERIC FORMATS

C

C     N1 IS STARTING ELEMENT WHICH WILL BE AFTER THE DATA COUNT

C     VALUE(NEED NELS*N3+1 TOTAL ELEMENTS)

      N1=N1+1

      N=0

C

C     # ELEMENTS IS EQUAL TO FORMAT NUMBER

      NELS=NFRM

C     READ UNITS CONVERSION SPECIFICATION LINE

      READ(LINE,1018) (PU(I),I=1,3)

      IF(PU(2).NE.' ') READ(LINE,2024) DELIM

      IF(PU(2).EQ.' ') READ(LINE,1020) DELIM,(PU(I),I=1,NELS)

C

C  ALLOW FOR BLANK CONVERSION LINE PRODUCED BY EIGENRAY ROUTINE

      IF(DELIM.EQ.0) READ(LINE,2024) DELIM

C

C     READ NEXT LINE AND EXTRACT UNIT CONVERSION SPECS IF PRESENT

C     IN FLOATING FORMAT

      READ(NRYIND,'(A)') LINE

      CALL GTUNIT(LINE,PU,NFIELD,NELS, NCARY)

C     CARRY COUNT IS ZERO IF LINE CONTAINED UNITS SPECIFIERS

      IF(NCARY.EQ.0.AND.NDEVTMP.GT.0) WRITE(NDEVTMP,'(A)') LINE

C

1018  FORMAT(10(A,5X))

1020  FORMAT(BZ,3X,G10.3,8X,10A)

2024  FORMAT(T43,G10.3)

1025  FORMAT(BZ,10G13.6)

C

C     LOOK UP CONVERSION CONSTANTS

      DO 2030 I=1,NELS

2030  uucon= UCON(IPV(I),PU(I),CONV(I))

C

C     DELIMITOR TESTING REQUIRES THE INNER LOOP BE EXECUTED

C     AT LEAST ONCE

      NX=MAX0(1,N3)

      N2=MAX0(N2,N1)

C     BEGIN SINGLE DATA FORMAT LOOP

      DO 150 J=0,10000

         IF(NELS.GT.1) THEN

            IF(NCARY.EQ.0) READ(NRYIND,'(A)') LINE

            NCARY=0

            IF(NDEVTMP.GT.0) WRITE(NDEVTMP,'(A)') LINE

            READ(LINE,1025) (V(I),I=1,NELS)

         ENDIF

         K=0

         DO 145 I=N1+J,N2+J,NX

               K=K+1

               IF(NELS.EQ.1) THEN

                  IF(NCARY.EQ.0) READ(NRYIND,'(A)',END=10000) LINE

                  NCARY=0

                  IF(NDEVTMP.GT.0) WRITE(NDEVTMP,'(A)') LINE

C     FOR SINGLE COLUMN CASE READ NOW

                  READ(LINE,1020) V(K)

               ENDIF

               IF(V(K) .EQ. DELIM) GO TO 160

C

C           APPLY ANY SCALING CONVERSIONS

               IF(CONV(K).GE.0.d0) THEN

                  V(K)=V(K)*CONV(K)

               ELSEIF(V(K).NE.0.d0) THEN

                  V(K)=-CONV(K)/V(K)

      ENDIF

               IF(J.LT.N3) GP(I)=V(K)

145            N=N+1

150   CONTINUE

C

10000 IF(V(1).NE.DELIM) CALL RERROR('READW1','NO DELIM',FLOAT(MSET))

C

C     SAVE NUMBER OF VALUES AND ELEMENTS(AS FRACTIONAL PART)

160   IF(N1.GT.1) GP(N1-1)=N+NELS/10.

C

C     CONTINUE MULTI-FORMAT LOOP

      GO TO 5

C

200   CALL RERROR('READW1','BAD FORMAT',FLOAT(MSET))

C

      END

 

 

 

 

 

C**********************************************************************

      SUBROUTINE GTUNIT(LINE,PU,NFIELD,NVAR, LN)

C**********************************************************************

C     GTUNIT IS USED TO EXTRACT UNIT CONVERSION SPECIFIERS FROM TABULAR

C     DATA INPUT.

C     CONVERSION SPECIFICATIONS ARE TAKEN FROM 'LINE' INPUT STRING.

C     A FLOATING FORMAT FALLING WITHIN SUCCESSIVE FIELDS OF

C     LENGTH 'NFIELD' IS PROVIDED FOR.

C     OUTPUT IS TO UNIT SPECIFIER ARRAY 'PU' FOR A MAXIMUM OF 'NVAR'

C     FIELDS. FOR A SUCCESS/FAIL RETURN VALUE, 'LN' GIVES THE NUMBER

C     OF NON-BLANK CHARACTERS USEFUL FOR UNIT DECODING, ZERO IMPLIES

C     NO USEFUL UNIT INFORMATION WAS FOUND.

C

 

      implicit double precision (a-h,o-z)

 

      PARAMETER (MAXLN=130)

      CHARACTER*(*) LINE,PU(NVAR),TMP*(MAXLN),TMP1*(MAXLN)

      INTEGER STRIM

C

      LN=STRIM(LINE)

      IF(LN.LT.1) THEN

         DO 5 I=1,NVAR

5           PU(I)=' '

         RETURN

      ENDIF

C

      TMP=' '

      CALL SFILTR(LINE,TMP,' 0123456789+-.E')

      IF(STRIM(TMP).LT.1) RETURN

C     SIGNIFY WE CONSUMED THIS LINE

      LN=0

C

      TMP1=LINE

C     PLACE INTER-SPECIFIER DELIMITORS

      DO 10 I=NFIELD,NFIELD*NVAR,NFIELD

         IF(TMP1(I:I).NE.' ') THEN

            IF(TMP1(I+1:I+1).NE.' ') THEN

C     IF NO ROOM FOR A BLANK, ITS AN ERROR

               CALL RERROR('GTUNIT','FORMAT ER1',FLOAT(I))

            ELSE

               TMP1(I+1:I+1)='#'

            ENDIF

         ELSE

            TMP1(I:I)='#'

         ENDIF

10    CONTINUE

C

C     NOW SQUEEZE OUT THE BLANKS

      CALL SFILTR(TMP1,TMP,' ')

C     HOPEFULLY SOMETHING LEFT

C

      NL=STRIM(TMP)

      ILIN=1

      DO 20 I=1,NVAR

C     CHECK FOR EMPTY CONVERSIONS

         IF(TMP(ILIN:ILIN).EQ.'#') THEN

            ILIN=ILIN+1

            PU(I)=' '

            GO TO 20

         ENDIF

C     GET NEXT UNIT TYPE

         PU(I)(1:2)=TMP(ILIN:ILIN+1)

C     GET NEXT PHYSICAL UNIT(BLANKS HAVE BEEN REMOVED)

C     SHORT PHRASE IF A BLANK WAS REMOVED

         IF(TMP(ILIN+3:ILIN+3).NE.'#') THEN

            IF(TMP(ILIN+4:ILIN+4).NE.'#') THEN

C     IF INCORRECT GROUPING OF LETTERS, ITS AN ERROR

               CALL RERROR('GTUNIT','FORMAT ER2',FLOAT(ILIN))

            ELSE

C     ADD TWO LETTER UNIT

               PU(I)(3:5)=' '//TMP(ILIN+2:ILIN+3)

               ILIN=ILIN+5

            ENDIF

         ELSE

C     ADD ONE LETTER UNIT

            PU(I)(3:5)='  '//TMP(ILIN+2:ILIN+2)

            ILIN=ILIN+4

         ENDIF

20    CONTINUE

C

C     SHOULD HAVE USED UP ALL THE TEXT BY NOW

      IF(ILIN.LE.NL) CALL RERROR('GTUNIT','FORMAT ER3',FLOAT(ILIN))

      END

 

 

 

 

 

C**********************************************************************

      SUBROUTINE SREAD1(AB,NWOK,NW)

C**********************************************************************

 

      implicit double precision (a-h,o-z)

      logical ab,nwok

C     THIS IS THE FALLTHROUGH FEATURE FOR READW WHEN ENCOUNTERING

C     A NON STANDARD LABELED COMMON BLOCK. SREAD1 WILL ADD TABULAR

C     DATA TO THE DEFAULT GP ARRAY PGP AND INCLUDE INDEX OFFSET

C     VALUES TO A LOOKUP ARRAY 'LIST'.

C     COMMON DECK "CPROCFL" INSERTED HERE

      INTEGER PMX,PNTBL,PITBL,PFRMTBL,IDSP(10)

C     PARAMETER DECK "PGROUPS"

      PARAMETER (NCHPG1=11,NWPV=250,NSPGP=NCHPG1+2*NWPV+1)

      PARAMETER (MNGRP=9,MXGRP=69,MXLIST=MXGRP-MNGRP+2)

      COMMON/PROCFL/LIST(MXLIST)

      COMMON/PROCFL/PMX,PNTBL(10),PITBL(10),PFRMTBL(10),PGP(NSPGP)

      EQUIVALENCE (PGP,IDSP)

C

      CHARACTER ITOC*7

C

      INTEGER PPX    ,PQTBL(10),PLTBL(10),PIRMTBL(10)

C

      DATA PPX/2/

      DATA PQTBL/1,NCHPG1,NSPGP,7*0/

      DATA PLTBL/1,NWPV,8*0/

      DATA PIRMTBL/1,2,8*0/

C

      NWA=IABS(NW)-(MNGRP-2)

      IF(NWA.LT.2 .OR. NWA.GT.MXLIST) CALL STOPIT

     1  ( 'SREAD1 NW='//ITOC(NW)//',LIST(1)='//ITOC(LIST(1)) )

C

      CALL READW1(AB,NWOK,NW,PMX,PNTBL,PITBL,PFRMTBL,PGP)

C

C     UPDATE THE LIST, LIST(1)=MAXIMUM USED ELEMENT

      LIST(1)=MAX0(LIST(1),NWA)

      LIST(NWA)=PNTBL(1)

C

C     GET NEXT NUMERICAL INDEX FROM THE POINTER ARRAY

      NLN=PNTBL(2)

C

C     NUMBER OF VALUES IS INTEGRAL PART OF FIRST ELEMENT

      NPTS=PGP(NLN)

C     NUMBER OF ROWS IS KEPT AS FRACTIONAL PART

      NELS=(PGP(NLN)-NPTS)*10.+.5

      N2=NPTS/NELS

      NMX=PITBL(2)

C     MAKE ADJUSTMENTS TO LENGTHS AND OFFSETS

      PITBL(2)=NMX-N2-(NCHPG1-1)

      NUP=NPTS+NCHPG1

      PNTBL(1)=PNTBL(1)+NUP

      PNTBL(2)=PNTBL(2)+NUP

      NLN=NLN+1

C

C     LOOP TO MOVE NELS-1 ROWS OF PGP(NMX,NELS) TO PGP(N2,NELS)

C     THUS MAKING PGP A CONTIGUOUS ARRAY

      NTO=NLN+N2

      NFRM=NLN+NMX

      DO 100 I=2,NELS

      IF(NTO.GT.0 .AND. NFRM+N2-1.LE.NSPGP)

     1        CALL RMOVE(PGP(NTO),PGP(NFRM),N2)

      NTO=NTO+N2

100   NFRM=NFRM+NMX

      RETURN

C

C**********************************************************************

      ENTRY SETSRD

C**********************************************************************

c%      CALL ISET2(LIST,0,MXLIST)

      do 66 i=1,mxlist

66    list(i)=0

      PMX=PPX

      CALL IMOVE(PNTBL,PQTBL,10)

      CALL IMOVE(PITBL,PLTBL,10)

      CALL IMOVE(PFRMTBL,PIRMTBL,10)

      END

 

 

 

 

 

C**********************************************************************

      FUNCTION   READW (XX)

C**********************************************************************

C     READS W ARRAY MAKING ANY NECESSARY CONVERSION OF UNITS.

C     VARIABLE LENGTH 'TABULAR' DATA IS READ BY SUBROUTINE 'READW1'

C     WHICH IS CALLED FOR NEGATIVE W-ARRAY INDICES.

C

C     COMMON DECK "CONST" INSERTED HERE

 

      implicit double precision (a-h,o-z)

 

      COMMON/PCONST/CREF,RGAS,GAMMA

      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10

C     COMMON DECK "WWR" INSERTED HERE

      COMMON/WW/DUM(10),MAXW,W(1000)

      COMMON/IDENT/ID

      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))

C     COMMON DECK "B1" INSERTED HERE

      INTEGER UMX,UNTBL,UITBL,UFRMTBL,IDSU(20)

      COMMON/B1/UMX,UNTBL(10),UITBL(10),UFRMTBL(10),UGP(10)

      EQUIVALENCE (UGP,IDSU)

C     COMMON DECK "B2" INSERTED HERE

      INTEGER DUMX,DUNTBL,DUITBL,DUFRMTB,IDSDU(20)

      COMMON/B2/DUMX,DUNTBL(10),DUITBL(10),DUFRMTB(10),DUGP(10)

      EQUIVALENCE (DUGP,IDSDU)

C     COMMON DECK "B3" INSERTED HERE

      INTEGER CMX,CNTBL,CITBL,CFRMTBL,IDSC(20)

      COMMON/B3/CMX,CNTBL(10),CITBL(10),CFRMTBL(10),CGP(512)

      EQUIVALENCE (CGP,IDSC),(ANC,CGP(11))

C     COMMON DECK "B4" INSERTED HERE

      INTEGER DCMX,DCNTBL,DCITBL,DCFRMTB,IDSDC(20)

      COMMON/B4/DCMX,DCNTBL(10),DCITBL(10),DCFRMTB(10),DCGP(10)

      EQUIVALENCE (DCGP,IDSDC)

C     COMMON DECK "B5" INSERTED HERE

      INTEGER TMX,TNTBL,TITBL,TFRMTBL,IDST(20)

      COMMON/B5/TMX,TNTBL(10),TITBL(10),TFRMTBL(10),TGP(262)

      EQUIVALENCE (TGP,IDST),(ANT,TGP(11))

C     COMMON DECK "B6" INSERTED HERE

      INTEGER DTMX,DTNTBL,DTITBL,DTFRMTB,IDSDT(20)

      COMMON/B6/DTMX,DTNTBL(10),DTITBL(10),DTFRMTB(10),DTGP(10)

      EQUIVALENCE (DTGP,IDSDT)

C     COMMON DECK "B7" INSERTED HERE

      INTEGER MMX,MNTBL,MITBL,MFRMTBL,IDSM(20)

      REAL*8 MGP

      COMMON/B7/MMX,MNTBL(10),MITBL(10),MFRMTBL(10),MGP(10)

      EQUIVALENCE (MGP,IDSM)

C     COMMON DECK "B9" INSERTED HERE

      INTEGER GMX,GNTBL,GITBL,GFRMTBL,IDSG(20)

      COMMON/B9/GMX,GNTBL(10),GITBL(10),GFRMTBL(10),GGP(113)

      EQUIVALENCE (GGP,IDSG),(ANG,GGP(11))

C

C     COMMON DECK "B10" INSERTED HERE

      INTEGER DGMX,DGNTBL,DGITBL,DGFRMTB,IDSDG(20)

      COMMON/B10/DGMX,DGNTBL(10),DGITBL(10),DGFRMTB(10),DGGP(10)

      EQUIVALENCE (DGGP,IDSDG)

C

C     COMMON DECK "B8" INSERTED HERE

      INTEGER RMX,RNTBL,RITBL,RFRMTBL,IDSR(20)

      COMMON/B8/RMX,RNTBL(10),RITBL(10),RFRMTBL(10),RGP(10)

      EQUIVALENCE (RGP,IDSR)

C

C     COMMON DECK "CB11" INSERTED HERE

      INTEGER SMX,SNTBL,SITBL,SFRMTBL,IDSS(20)

      COMMON/B11/SMX,SNTBL(10),SITBL(10),SFRMTBL(10),SGP(11)

      EQUIVALENCE (SGP,IDSS),(ANS,SGP(11))

C     COMMON DECK "CB12" INSERTED HERE

      INTEGER DSMX,DSNTBL,DSITBL,DSFRMTB,IDSDS(20)

      COMMON/B12/DSMX,DSNTBL(10),DSITBL(10),DSFRMTB(10),DSGP(11)

      EQUIVALENCE (DSGP,IDSDS),(ANDS,DSGP(11))

C     COMMON DECK "CB17" INSERTED HERE

      INTEGER VMX,VNTBL,VITBL,VFRMTBL,IDSV(20)

      COMMON/B17/VMX,VNTBL(10),VITBL(10),VFRMTBL(10),VGP(53)

      EQUIVALENCE (VGP,IDSV),(ANV,VGP(11))

C     COMMON DECK "CB18" INSERTED HERE

      INTEGER DVMX,DVNTBL,DVITBL,DVFRMTB,IDSDV(20)

      COMMON/B18/DVMX,DVNTBL(10),DVITBL(10),DVFRMTB(10),DVGP(11)

      EQUIVALENCE (DVGP,IDSDV),(ANDV,DVGP(11))

C     COMMON DECK "CB19" INSERTED HERE

      INTEGER PRMX,PRNTBL,PRITBL,PRFRMTB,IDSPR(20)

      COMMON/B19/PRMX,PRNTBL(10),PRITBL(10),PRFRMTB(10),PRGP(11)

      EQUIVALENCE (PRGP,IDSPR),(ANP,PRGP(11))

C     COMMON DECK "CB20" INSERTED HERE

      INTEGER DPMX,DPNTBL,DPITBL,DPFRMTB,IDSDP(20)

      COMMON/B20/DPMX,DPNTBL(10),DPITBL(10),DPFRMTB(10),DPGP(11)

      common/crmach/rmach(9)

      COMMON/HDRC/INITID,DAT,TOD    

      common/hdr/sec,stsec                                 !Rev 7/13/90

      EQUIVALENCE (DPGP,IDSDP),(ANDP,DPGP(11)),(bigval,rmach(6))

      PARAMETER (MXCMTS=83)

      CHARACTER*38 WFRMT,WNOTES(2),bell*1,tod*10,face*1    !Rev 7/13/90

      character*2 numstg

      CHARACTER*80 LINEX

      CHARACTER*100 STMP1,STMP2

      INTEGER STRIM,WCOMTS(MXCMTS)

      CHARACTER*1 DEG,KM,NM,FEET,CYCLE,PER,MSKMS

C     COMMON DECK "RAYDEV" INSERTED HERE

C     DEVICE ASSIGNED TO RAYTRC INPUT FILE

      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVGRP,NDEVBIN

      LOGICAL NWOK,AB,uucon,ucon,setucon

C

C     COMMON DECK "FLAG" INSERTED HERE

      LOGICAL NEWWR,NEWWP,NEWTRC,PENET

      COMMON /FLG/ NTYP,NEWWR,NEWWP,NEWTRC,PENET,LINES,IHOP,HPUNCH

      COMMON/FLGP/NSET

c	Common deck ' CGRID2 ' inserted here		A. Newhall

	parameter (mx=7,my=7,mz=7)

	common/cgrid2/lxx7,nyy7,mzz7,xc7(7),yc7(7),zc7(7),cavg7(7)

	common/knots/zknot(11),yknot(11),xknot(11)

	common/coefs/cco(7),spleof(7,7),cbispl(7,7*7)

C

      CHARACTER PC*1,ID*80,typ*20

C

      INTEGER XMX,XNTBL(10),XITBL(10),XFRMTBL(10)

      REAL*8 XGP(11)

C

      DATA XMX/2/

      DATA XNTBL/1,11,12,7*0/

      DATA XITBL/1,9*0/

      DATA XFRMTBL/1,2,8*0/

C

      DATA WCOMTS,WNOTES/MXCMTS*1, ' ',' INPUT OVERRIDDEN'/

*C

      AB=.FALSE.

      READW    =0.d0

C

      IF(NDEVTMP.GT.0) REWIND NDEVTMP

C

      READ(NRYIND,'(a)',END=3) ID

      IF(NDEVTMP.GT.0) WRITE(NDEVTMP,'(a)') ID

      print *,id

C*      READ(NRYIND,1000,END=3) ID

C* 1000 FORMAT (BZ,10A8)

C*      IF(NDEVTMP.GT.0) WRITE(NDEVTMP,1000) ID

      GO TO 4

C

3     READW    =1.d0

      IF(NFRMAT.LT.0) RETURN

C

C%      PRINT *,'  END OF INPUT DATA '     !Rev 7/13/90

      CALL ENDPLT

33    IF(PLT.NE.0.d0) CALL DDEND

      bell=16#07                 !Rev 7/13/90 here thru STOP:

      face=16#02                 ! to display elapsed time and 

      call systim(tod)           ! time completed.

      call syssec(sec) 

      elsec=(sec-stsec)

      if(sec.lt.stsec)elsec=(sec-stsec+86400.)  !add 24 hr

      elmin=elsec/60

      elhrs=elmin/60

      write (*,886)int(elhrs),int(mod(elmin,60)),mod(elsec,60)

886   format(/'  END OF INPUT DATA...',22x,'ELAPSED TIME = ',

     1  i2,' HR ',i2,' MIN ',f4.1,' SEC')

      write (*,888)bell,face,face,face,tod,face,face,face  

888   format (//12x,a1,a1,1x,a1,1x,a1,

     1   '   HARPO RUN COMPLETED NORMALLY AT ',A8,2x,a1,1x,a1,1x,a1)

      STOP '  '

C

4     READ(NRYIND,'(A)',END=10001) STMP1

C

      READ(STMP1,1100) NW,WWW,LINEX

1100  FORMAT (BZ,I3,E14.7,A)

      IF(NDEVTMP.GT.0.AND.(NFRMAT.NE.-2.OR.NW.LT.0))

     1      WRITE(NDEVTMP,'(A)') STMP1

10001 IF (NW.EQ.0) GO TO 10

      IF(NW.GT.MAXW) GO TO 3400

C

      IF (NW.LE.MAXW .AND. NW.GT.0) GO TO 5

C     TABULAR INPUT DATA

      NWOK=.FALSE.

C 'OPEN' THE TEMP FILE

1150  FORMAT(I3,T18,A)

      NWP=-NW

      IF(NFRMAT.EQ.-2) THEN

C     USE DUMMY ARRAYS TO ABSORB TABULAR DATA

         CALL READW1(AB,NWOK,NW,XMX,XNTBL,XITBL,XFRMTBL,XGP)

      ELSEIF(NWP.EQ.1) THEN

         CALL READW1(AB,NWOK,NW,UMX,UNTBL,UITBL,UFRMTBL,UGP)

      ELSEIF(NWP.EQ.2) THEN

         CALL READW1(AB,NWOK,NW,DUMX,DUNTBL,DUITBL,DUFRMTB,DUGP)

      ELSEIF(NWP.EQ.3) THEN

         CALL READW1(AB,NWOK,NW,CMX,CNTBL,CITBL,CFRMTBL,CGP)

      ELSEIF(NWP.EQ.4) THEN

         CALL READW1(AB,NWOK,NW,DCMX,DCNTBL,DCITBL,DCFRMTB,DCGP)

      ELSEIF(NWP.EQ.5) THEN

            CALL READW1(AB,NWOK,NW,TMX,TNTBL,TITBL,TFRMTBL,TGP)

      ELSEIF(NWP.EQ.6) THEN

         CALL READW1(AB,NWOK,NW,DTMX,DTNTBL,DTITBL,DTFRMTB,DTGP)

      ELSEIF(NWP.EQ.7) THEN

         CALL READW1(AB,NWOK,NW,MMX,MNTBL,MITBL,MFRMTBL,MGP)

      ELSEIF(NWP.EQ.8) THEN

         CALL READW1(AB,NWOK,NW,RMX,RNTBL,RITBL,RFRMTBL,RGP)

      ELSEIF(NWP.EQ.9) THEN

         CALL READW1(AB,NWOK,NW,GMX,GNTBL,GITBL,GFRMTBL,GGP)

      ELSEIF(NWP.EQ.10) THEN

         CALL READW1(AB,NWOK,NW,DGMX,DGNTBL,DGITBL,DGFRMTB,DGGP)

      ELSEIF(NWP.EQ.11) THEN

         CALL READW1(AB,NWOK,NW,SMX,SNTBL,SITBL,SFRMTBL,SGP)

      ELSEIF(NWP.EQ.12) THEN

         CALL READW1(AB,NWOK,NW,DSMX,DSNTBL,DSITBL,DSFRMTB,DSGP)

      ELSEIF(NWP.EQ.17) THEN

         CALL READW1(AB,NWOK,NW,VMX,VNTBL,VITBL,VFRMTBL,VGP)

      ELSEIF(NWP.EQ.18) THEN

         CALL READW1(AB,NWOK,NW,DVMX,DVNTBL,DVITBL,DVFRMTB,DVGP)

      ELSEIF(NWP.EQ.19) THEN

         CALL READW1(AB,NWOK,NW,PRMX,PRNTBL,PRITBL,PRFRMTB,PRGP)

      ELSEIF(NWP.EQ.20) THEN

         CALL READW1(AB,NWOK,NW,DPMX,DPNTBL,DPITBL,DPFRMTB,DPGP)

      ELSE

         CALL SREAD1(AB,NWOK,NW)

      ENDIF

      IF(NWOK) GO TO 4

C

3400  WRITE(3,4000) NW,MAXW

4000  FORMAT(15H1THE SUBSCRIPT ,I3,  ' ON THE W-ARRAY INPUT IS OUT OF BO

     1UNDS.  ALLOWABLE VALUES ARE -8 THROUGH,',I4,'.' )

      STOP 4000

C

5     READ(LINEX,70) DEG,KM,NM,FEET,CYCLE,PER,MSKMS

70    FORMAT(7A)

      W(NW)=WWW

C

C     CHECK FOR A 'TERRAIN RELATIVE' HEIGHT SPECIFICATION.

C**   IF SO ADD FLAG VALUE 10**40 TO BE TESTED FOR LATER IN 'RAYTRC'

C     IF SO ADD FLAG VALUE 10**25 TO BE TESTED FOR LATER IN 'RAYTRC'

C

C**   IF(MSKMS.EQ.'T' .AND. WWW.EQ.0.0) W(NW)=1.E-40

      IF(MSKMS.EQ.'T' .AND. WWW.EQ.0.d0) W(NW)=1.d-10

      IF(WWW.EQ.0.d0) GO TO 4

C**   IF(MSKMS.EQ.'T') WWW=WWW*1.E40

      IF(MSKMS.EQ.'T') WWW=WWW*1.d10

C

C     CHECK FOR KEYWORD UNITS SPECIFICATION, IF SO PERFORM CONVERSION

C

      uucon=  UCON(KU,LINEX(:10),CONV)

      if (.not. uucon) go to 60

75    IF(CONV.GE.0.d0) THEN

         W(NW)=WWW*SNGL(CONV)

      ELSEIF(W(NW).NE.0.d0) THEN

         W(NW)=-CONV/W(NW)

      ENDIF

      GO TO 4

C

60    IF (DEG.EQ.'1') WWW=WWW*RAD

      IF (KM.EQ.'1') WWW=WWW/EARTHR

      IF (NM.EQ.'1') WWW=WWW*1.852d0

      IF (FEET.EQ.'1') WWW=WWW*3.048006096d-4

      IF(CYCLE.EQ.'1') WWW=WWW*PIT2

      IF(PER.EQ.'1') WWW=PIT2/WWW

      IF(MSKMS.EQ.'1') WWW=WWW*1.d-3

      W(NW)=WWW

      GO TO 4

C

10    IF(.NOT.AB) RETURN

      PRINT 1200

1200  FORMAT(' A DATA FORMATTING ERROR PREVENTS CONTINUED EXECUTION')

      STOP 4500

C

C**********************************************************************

      ENTRY SETW(XX)

C**********************************************************************

C     THIS ENTRY IS CALLED ONCE BEFORE THE FIRST RUN SET IS READ

C

C INITIALIZE SOME MATHEMATICAL CONSTANTS

      PI=4.d0*dATAN(1.d0)

      PIT2=2.d0*PI

      PID2=PI/2.d0

      DEGS=180.d0/PI

      RAD=PI/180.d0

      ALN10=dLOG(10.d0)

CC********* INITIALIZE SOME VARIABLES IN THE W ARRAY

      MAXW=1000

      CALL CLEAR(W,MAXW)

      PLON=0.d0

      PLAT=PID2

      EARTHR=6370.d0

      INTYP=3.d0

      MAXERR=1.d-4

      ERATIO=50.d0

      STEP1=1.d0

      STPMAX=100.d0

      STPMIN=1.d-8

      FACTR=5.d-1

      HITLET=.15d0

      MAXSTP=1000.d0

      HOP=1.d0

      HMAX=500.d0

      EXTINC=999.999d0

      PAGLN=66.d0

C

C%      CALL SETUCON

      uucon=setucon(x,y,z)

      RETURN

C

C**********************************************************************

      ENTRY SETOVR (XX)

C**********************************************************************

C     THIS ENTRY IS CALLED AFTER EACH RUN SET IS READ

C     PERFORM OVERRIDE ASSIGNMENTS TO W-ELEMENTS

      CALL OVERRD(MAXSTP,0.d0,1.d3,WCOMTS(23),2,1)

      CALL OVERRD(ERATIO,0.d0,5.d1,WCOMTS(43),2,1)

      CALL OVERRD(FACTR,0.d0,5.d-1,WCOMTS(47),2,1)

      CALL OVERRD(SKIP,0.d0,bigval,WCOMTS(71),2,1)

      CALL OVERRD(HITLET,0.d0,1.5d-1,WCOMTS(75),2,1)

      IF(PAGLN.LT.3.d1) PAGLN=0.d0

      CALL OVERRD(PAGLN,0.d0,6.6d1,WCOMTS(77),2,1)

      CALL OVERRD(PFACTR,0.d0,1.d0,WCOMTS(82),2,1)

      RETURN

C

C**********************************************************************

      ENTRY PRINTW (XX)

C**********************************************************************

C     THIS ENTRY IS CALLED AFTER EACH RUN SET IS READ

C     TO PRINT VALUES OF THE W-ARRAY IN A FORMAT WHICH SHOWS

C     FULL ACCURACY AND THE EFFECTS OF ANY CONVERSIONS OR OVERRIDES.

C

C     GO TO NEXT PAGE, PUT OUT COPY OF INPUT DATA FOR THIS RUN SET

C     ASSUMING WE HAD A TEMPORARY FILE TO USE

C

      IF(NDEVTMP.LE.0) GO TO 1065

      CALL NEWPAG(NPAG,idint(pagln),PC)

      LNSXPG=LINES+idint(pagln)

      CALL PUTHDR(3,PC,NPAG)

      CALL PUTDVR(3)

      CALL PUTKBK(3,1)

C%      CALL PUTKCT(3,'INPUT DATA FILE FOR RUN SET NUMBER'

C%     1   //NUMSTG(NSET,1,'(I5)'))

      WRITE (3,1055)NSET

1055  FORMAT (47X,'INPUT DATA FILE FOR RUNSET NUMBER ',I3)

      CALL PUTKBK(3,1)

      CALL PUTDVR(3)

      CALL PUTKBK(3,1)

      REWIND NDEVTMP

C

1060  READ(NDEVTMP,'(A)',END=1065) LINEX

      IF(LINES.GE.LNSXPG) THEN

         CALL NEWPAG(NPAG,idint(pagln),PC)

         LNSXPG=LINES+idint(pagln)

         CALL PUTHDR(3,PC,NPAG)

         CALL PUTKBK(3,1)

      ENDIF

      CALL PUTKST(3,' '//LINEX)

      GO TO 1060

1065  CONTINUE

C

C     GO TO NEXT PAGE, PUT OUT W-ELEMENTS SHOWING FULL PRECISION

C     USE A TWO COLUMN FORMAT

C

      CALL NEWPAG(NPAG,idint(pagln),PC)

      CALL PUTHDR(3,PC,NPAG)

      CALL PUTDVR(3)

      CALL PUTKBK(3,1)

      CALL PUTKCT(3,'INITIAL VALUES FOR THE W ARRAY')

      CALL PUTKCT(3,'ONLY NONZERO VALUES PRINTED'

     1  //' -- ALL ANGLES IN RADIANS')

      CALL PUTKBK(3,1)

      CALL PUTDVR(3)

      LINSPP=idint(pagln)-7

C

      WFRMT='(BZ,I3,E14.7,A)'

      IF(NFRMAT.EQ.0) WFRMT='(I4,E24.15,A)'

C

      NWE=0

C     ALLOW MAXIMUM OF 10 PAGES OF OUTPUT(1 IS ENOUGH)

      DO 18 NPGS=1,10

C     ALLOW 3 LINES FOR INSIDE HEADER(SEE BELOW)

         LINSPP=LINSPP-3

C     NEXT ELEMENT TO SCAN IS ONE MORE THAN LAST ONE

         LWE=NWE+1

C

C     INITIALIZE SECOND COLUMN START INDEX

         NXCL=0

C     INITIALIZE TOTAL ELEMENT COUNTER

         NELS=0

C     LOOP TO FIND BREAK POINTS FOR FIRST/SECOND COLUMNS

         DO 14 I=LWE,MAXW

            IF(W(I).EQ.0.d0) GO TO 14

            NWE=I

            NELS=NELS+1

            IF(NELS.EQ.LINSPP) NXCL=I

            IF(NELS.EQ.2*LINSPP) GO TO 15

14       CONTINUE

15    IF(NELS.EQ.0) GO TO 22

C

         IF(NPGS.GT.1) THEN

            CALL NEWPAG(NPAG,idint(pagln),PC)

            CALL PUTHDR(3,PC,NPAG)

            LINSPP=idint(pagln)-1

         ENDIF

C

C     INSERT 'INNER' HEADER

         CALL PUTKBK(3,2)

         STMP1=' N              W(N)'

         CALL PUTKST(3,STMP1(:65)//' '//STMP1)

C

       NX=NXCL-1

       IF(NX.LT.LWE) NX=NWE

       DO 17 NW=LWE,NX

       NCOMT=MIN0(NW,MXCMTS)

       IF(W(NW).NE.0.d0) THEN

         WRITE(STMP1,WFRMT) NW,W(NW),WNOTES(WCOMTS(NCOMT))

         STMP2=' '

         IF(NXCL.EQ.0) GO TO 23

         DO 19 I=NXCL,MAXW

         IF(W(I).NE.0.d0) THEN

           NCOMT=MIN0(I,MXCMTS)

           WRITE(STMP2,WFRMT) I,W(I),WNOTES(WCOMTS(NCOMT))

           NXCL=I+1

           GO TO 23

         ENDIF

19       CONTINUE

23       CALL PUTKST(3,STMP1(:65)//' '//STMP2)

       ENDIF

17     CONTINUE

      IF(NELS.EQ.0) GO TO 22

18    CONTINUE

C

22    CONTINUE

      END

 

 







