      FUNCTION ALCOSH(X)

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



C     COMPUTE LOG(COSH(X)) AND USE LARGE ARGUMENT APPROXIMATION

C     WHEN POSSIBLE.

      DATA ALOG2/.6931471806/

C

      IF(ABS(X).GT.50.0) GO TO 10

      EX=EXP(X)

      ALCOSH=ALOG((EX+1.0/EX)*.5)

      RETURN

10    ALCOSH=ABS(X)-ALOG2

      END

C-----------------------------------------------------------------------------

      SUBROUTINE ARCTIC(THMIN,THMAX,HEIGHT,TICY,PLOT)

C     DRAW RANGE AXIS IN RAY TRACE PLOT. INCLUDES ANY CURVILINEAR

C     PROJECTIONS PROVIDED IN DDGRAPH.



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



C     COMMON DECK "WWR" INSERTED HERE

      PARAMETER (NWARSZ=1000)

      COMMON/WW/dum(10),MAXW,W(NWARSZ)

      EQUIVALENCE (TIC,W(87)),(HB,W(88)),(HT,W(89)),(TICV,W(96))



C.....The following statement makes the DDSCALE common size compatible

C     with later routines which have QQQ(8) as the first 8 elements.



      real xmin,xmax,ymin,ymax,scx,scy

      COMMON /DDSCALE/ XMIN,XMAX,YMIN,YMAX,MINX,MAXX,MINY,MAXY,SCX,SCY,

     1                 NSCX,NXCY,MSCX,MSCY,ISCX,ISCY

C

      NTIC=2

      IF(TIC.NE.0.) NTIC = 1+(THMAX-THMIN)/TIC

      NLINE=MAX0(2,100/NTIC)

C

      TICN=TIC/(NLINE-2)

      DO 10 I=1,NTIC

           X=THMIN+(I-1)*TIC

           CALL PLOT(X,HEIGHT+TICY,1)

           DO 10 J=2,NLINE

                XJ=X+(J-2)*TICN

                IF(XJ.GT.THMAX) GO TO 15

10              CALL PLOT(XJ,HEIGHT,0)

C



15    CALL PLOT(THMAX,HEIGHT,0)

      CALL PLOT(THMAX,HEIGHT+TICY,1)

      RETURN

      END

C-----------------------------------------------------------------------

      SUBROUTINE SETTRC

C     INITIALIZE /TRAC/ LABELED COMMON



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



C     COMMON DECK "TRAC" INSERTED HERE

      LOGICAL GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY

      COMMON /TRAC/ GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY,SMT,OSMT

      COMMON/TRAC/ROLD(20),DROLD(20),TOLD,ZDOT,D2Z,RAD,RAD1

      GROUND=.FALSE.

      SURF=.FALSE.

      END

C-------------------------------------------------------------------------

      SUBROUTINE CONHEAD(CL1,CL2,IOFF,NMD,MODEL1,MODEL2,MODEL3

     1  ,MODEL4,labelcon,lencon, ocean)



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



C

C.....Added LABELCON which contains the contour max,min and interval



      CHARACTER*(*) CL1,CL2

      CHARACTER*80 labelcon

      CHARACTER*50 CT1,CT2

      INTEGER CTRIM

ccccc      INTEGER MODEL1(4),MODEL2(4),MODEL3(4),MODEL4(4)

      double precision MODEL1(4),MODEL2(4),MODEL3(4),MODEL4(4)

C

C     PUT HEADERS ON CONTOUR PLOTS.

C

C     COMMON DECK "WWR" INSERTED HERE

      character*80 id

      common /ident/ id



C.....Set true when at or below sea level.  Prevents use of MODEL3 which

C     is the molecular weight model which only applies to atmosphere contours.



      logical ocean



      PARAMETER (NWARSZ=1000)

      COMMON/WW/dum(10),MAXW,W(NWARSZ)

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

C     COMMON DECK "HDR" INSERTED HERE

      CHARACTER*10 INITID*80,DAT,TOD

      COMMON/HDR/SEC

      COMMON/HDRC/INITID,DAT,TOD

      COMMON/DD/IN,IOR,IT,IS,IC,ICC,IX,IY

      CHARACTER*80 LABEL

C

      IX=0

      IOR=0

      IT=0

      IS=2

      IY=IOFF

      IF(IOFF.NE.0) GO TO 10

      IY=1023

      CT1=CL1

      CT2=CL2

      L1=CTRIM(CT1)

      L2=CTRIM(CT2)

C.....Set length to 5, then DDTEXT will always have enough characters

C     5 --> 50 in DDTEXT.  Used to be (L1+L2)/2 but this leads to

C     truncation!!



      CALL DDTEXT(5 , CT1(:L1)//CT2)

CCCCC      IX=1090

CCCCC      CALL DDTEXT (2,DAT)  ! DAT is the DATE from SYSDAT

      IX=0

      DTIC=TIC*EARTHR

C

10    KLM=3

      IF(NMD.GE.3) KLM=4



C.....51 will cause the ID to be on line 2 of the screen.



      iy = iy-51

      WRITE(LABEL,1000) ID

      CALL DDTEXT (8,LABEL)



C.....Put LABELCON (contour max,min,interval) on next line.

      iy = iy-34



      call ddtext (lencon/10 +1,labelcon)



C.....Put "MODELS" on the same line near the end of the line.

      ix=1050

C

C     NMD is number of models = 4 for sound speed contours;

C                             = 2 for wind/current (and components);

C                             = 1 for terrain contours.

C     There are probably more compact ways of doing the following, but

C     this way is easy to read and understand.



C     34 puts each model on a different line. 51 skips a line.



      write (label,1000) 'MODELS USED:'

      call ddtext (8,LABEL)

      iy=iy-51

      write (label,1700) model1(1),model1(2)

      CALL DDTEXT (8,LABEL)

      write (label,1700) model1(3),model1(4)

      iy = iy-34

      CALL DDTEXT (8,LABEL)

C

C.....RETURN if only 1 model i.e. for terrain contours.



      if (nmd .eq. 1) RETURN



      write (label,1700) model2(1),model2(2)

      iy = iy-34

      CALL DDTEXT (8,LABEL)

      write (label,1700) model2(3),model2(4)

      iy = iy-34

      CALL DDTEXT (8,LABEL)



C.....RETURN if only doing 2 models i.e. wind/current contours.



      if (nmd .eq. 2) RETURN



C.....Only write the molec.wt. model if doing atmos. contours.

      if (.not. ocean) then

           write (label,1700) model3(1),model3(2)

           iy = iy-34

           CALL DDTEXT (8,LABEL)

           write (label,1700) model3(3),model3(4)

           iy = iy-34

           CALL DDTEXT (8,LABEL)

      endif

      write (label,1700) model4(1),model4(2)

      iy = iy-34

      CALL DDTEXT (8,LABEL)

      write (label,1700) model4(3),model4(4)

      iy = iy-34

      CALL DDTEXT (8,LABEL)



      RETURN



1000  FORMAT (A)

1550  FORMAT('HEIGHT=',F8.2,' KM.')

1600  FORMAT(4(2(A8,F5.1,1X),2X))  !changed 10 to 8

1700  FORMAT((A8,F5.1,1X),2X)  !testing

      END

C-----------------------------------------------------------------------------

      SUBROUTINE COPYVL(V,W)

C

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



      DOUBLE PRECISION V(3),W(3,3)

C

      DO 10 I=1,3

10    W(1,I)=V(I)

      END

      SUBROUTINE CROSSP(V,W,R)

C

C     COMPUTE CROSS PRODUCT OF VECTORS V AND W , RESULT IN R

C

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



      DOUBLE PRECISION V(3),W(3),R(3)

C

      R(1)=V(2)*W(3)-V(3)*W(2)

      R(2)=-(V(1)*W(3)-V(3)*W(1))

      R(3)=V(1)*W(2)-V(2)*W(1)

      END

      INTEGER FUNCTION CTRIM(C)

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



      CHARACTER*(*) C

C

      L=LEN(C)

      DO 100 I=L,1,-1

100   IF(C(I:I) .NE. ' ') GO TO 200

      CTRIM=1

      RETURN

200   CTRIM=I

      END

C------------------------------------------------------------------------------

      SUBROUTINE DDCON (X, Y, Z, CC, MX, MY, MZ, MC)



C        VERSION 3.6

C

C***E     DDCON - CONTOUR PLOTTING

C        L. DAVID LEWIS     (NOAA-SEL)

C        THIS SUBROUTINE IS AN EXTENSION OF SUBROUTINE QUICK

C        BY DAVE ROBERTSON (NCAR).

C***1     CALL DDCON(X,Y,Z,C,MX,MY,MZ,MC)

C

C        X IS ARRAY OF MX X-VALUES.

C        Y IS ARRAY OF MY Y-VALUES.

C        Z IS ARRAY OF MX*MY FUNCTION VALUES.

C        MZ IS THE DEFINING ROW DIMENSION OF Z.

C        (( Z(I+MZ*(J-1)) = FUNCTION(X(I),Y(J)), I=1,MX), J=1,MY)

C        Cc IS ARRAY OF MC CONTOUR VALUES.  DDCONTUR OPERATES MORE

C          EFFICIENTLY IF C IS IN ASCENDING ORDER.. C(K).LE.C(K+1)

C        PLOT SCALING VALUES ARE SPECIFIED IN COMMON /DDSCALE/

C        INTENSITY IS SPECIFIED IN COMMON /DD/

C        USER MUST PROVIDE FRAME ADVANCE.

C***.

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



      DIMENSION X (MX), Y (MY), Z (MZ, MY), CC (MC)



C         COMMON DECK DD - DDPLOT PARAMETERS

      COMMON /DD/ IN,IOR,IT,IS,IC,ICC,IX,IY

C         IN - INTENSITY (0,1) FOR VECTORS, POINTS, SYMBOLS.

C         IOR - SYMBOL ORIENTATION (0,1) (UPRIGHT, ROTATED 90 DEG CCW).

C         IT - (0,1) (ROMAN, ITALIC) SYMBOL STYLE.

C         IS - (0,1,2,3) SYMBOL SIZE.

C         IC - (0,1) SYMBOL CASE (UPPER, LOWER).

C         ICC - (0-63) CHARACTER CODE.

C         IX - X-COORDINATE (0-1023).

C         IY - Y-COORDINATE (0-1023).

C         END COMMON DECK DD.



C..... added the following common

      INTEGER REZFLG

      PARAMETER (REZFLG=2001)

      COMMON/DDREZ/DDHIX,DDHIY



      integer*2 numpts, ixtmp(2000), iytmp(2000), bpxtmp(2000),

     &          bpytmp(2000)



C.....OUTFLAG is set by DDVC_ACUM if a plot point coordinate is out of range.



      logical outflag



      real qqq,qq,scx,scy



      COMMON /DDSCALE /QQQ (8), SCX, SCY, QQ (4), ISCX, ISCY



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

C

C        INTERPOLATION FUNCTIONS.

      VX (CA, CB) = X1 + DX * CA / (CA - CB)

      VY (CA, CB) = Y1 + DY * CA / (CA - CB)

C

      IF (MX .LT. 2 .OR. MY .LT. 2 .OR. MC .LT. 1) RETURN

ccccc      IF(MX.GT.MZ) CALL SYSTEM(52,' MX.GT.MZ')



      outflag = .false.



      ASSIGN 3 TO JV

      ASSIGN 2 TO JASC

      IF(MC.EQ.1) GO TO 11

      DO 101 K = 2, MC

      IF(CC(K).LT.CC(K-1)) GO TO 11

101   CONTINUE

      ASSIGN 3 TO JASC

11      DO 2 K = 1, MC

      C = CC (K)

      DO 3 I=2,MX

      X1 = X (I - 1)

      X2 = X (I)

      DX = X2 - X1

      DO 3 J = 2, MY

      Y1 = Y (J - 1)

      Y2 = Y (J)

      DY = Y2 - Y1

      Z1 = Z (I - 1, J - 1)

      Z2 = Z (I - 1, J)

      Z3 = Z (I, J)

      Z4 = Z (I, J - 1)

      ASSIGN 149 TO JD

      C1 = C - Z1

      C2 = C - Z2

      C3 = C - Z3

      C4 = C - Z4

C

C        BINARY TREE SEARCH FOR CONTOUR-EDGE INTERSECTIONS.

      IF (C1 .LT. 0.) GO TO 121

      IF (C2 .LT. 0.) GO TO 131

      IF (C3 .LT. 0.) GO TO 141

      IF (C4 .LT. 0.) GO TO 146

C        THIS CONTOUR LEVEL IS ABOVE ALL FOUR Z-VALUES, SO WE CAN

C        SKIP THE REST OF THE C-ARRAY IF IT IS IN ASCENDING ORDER.

      GO TO JASC, (2, 3)

 121  IF (C2 .LT. 0.) GO TO 133

      IF (C3 .LT. 0.) GO TO 147

      IF (C4 .LT. 0.) GO TO 144

 122  VX1 = X1

      VY1 = VY (C1, C2)

      VX2 = VX (C1, C4)

      VY2 = Y1

      GO TO 19

 131  IF (C3 .LT. 0.) GO TO 143

      IF (C4 .LT. 0.) GO TO 148

 132  VX1 = X1

      VY1 = VY (C1, C2)

      VX2 = VX (C2, C3)

      VY2 = Y2

      GO TO 19

 133  IF (C3 .LT. 0.) GO TO 145

      IF (C4 .LT. 0.) GO TO 142

 134  VX1 = VX (C2, C3)

      VY1 = Y2

      VX2 = VX (C1, C4)

      VY2 = Y1

      GO TO 19

 141  IF (C4 .LT. 0.) GO TO 134

 142  VX1 = VX (C2, C3)

      VY1 = Y2

      VX2 = X2

      VY2 = VY (C4, C3)

      GO TO 19

 143  IF (C4 .LT. 0.) GO TO 122

 144  VX1 = X1

      VY1 = VY (C1, C2)

      VX2 = X2

      VY2 = VY (C4, C3)

      GO TO 19

 145  IF (C4 .LT. 0.) GO TO 3  !!!!was 2

 146  VX1 = X2

      VY1 = VY (C4, C3)

      VX2 = VX (C1, C4)

      VY2 = Y1

      GO TO 19

 147  IF (C4 .LT. 0.) GO TO 132

C

C        FOUR INTERSECTIONS.  EVALUATE HYPERBOLIC DISCRIMINANT.

 148  GO TO JD, (149, 150)

 149  D = Z1 - Z2 + Z3 - Z4

      E = (Z2 - Z1) * (Z4 - Z1)

      ASSIGN 150 TO JD

 150  Q = C1 * D + E

      IF (ABS (Q) .LE. .01 * ABS (D)) GO TO 152

      IF (Q .GT. 0.) GO TO 153

 151  ASSIGN 154 TO JV

      GO TO 132

 152  ASSIGN 155 TO JV

      GO TO 134

 153  ASSIGN 156 TO JV

      GO TO 122

 154  ASSIGN 3 TO JV  !!!! was 2

      GO TO 146

 155  ASSIGN 3 TO JV    !!!

      GO TO 144

 156  ASSIGN 3 TO JV !!!

      GO TO 142

C

C         THE NEXT 12 STATEMENTS PLOT A VECTOR FROM VX1,VY1 TO VX2,VY2.

C

19    SCX = VX1

      SCY = VY1

      CALL DDSC

      IX = ISCX

      IY = ISCY



c*********************************



      IF(IX.EQ.REZFLG) THEN

         IX=DDHIX

         IY=DDHIY

      ELSE

         DDHIX=IX

         DDHIY=IY

      ENDIF

C

      MNIX=MIN0(MNIX,IX)

      MXIX=MAX0(MXIX,IX)

      MNIY=MIN0(MNIY,IY)

      MXIY=MAX0(MXIY,IY)



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

C DDVC_ACUM increments NUMPTS.  But it does so, before putting values into

C IXTMP,IYTMP.  Therefore to have only one place incrementing NUMPTS, we

C have to handle the fact that BPXTMP is filled before NUMPTS has been

C incremented, by always adding 1 to the current value of NUMPTS.



           bpxtmp(numpts+1) = int(1.d1*ddhix)

           bpytmp(numpts+1) = int(1.d1*ddhiY)



ccccc           CALL DDBP



      SCX = VX2

      SCY = VY2

      CALL DDSC

      IX = ISCX

      IY = ISCY



      CALL DDVC_ACUM (numpts,ixtmp,iytmp, outflag)



C.....When 2000 points have been accumulated, write the moveto,linetos

C     to TXTOUT and reset NUMPTS to zero.



      if (numpts.ge.2000) then

           call flush (numpts,ixtmp,iytmp,bpxtmp,bpytmp)

      endif



      GO TO JV, (3, 154, 155, 156)

   3  CONTINUE



      if (numpts .gt. 0) call flush (numpts,ixtmp,iytmp,bpxtmp,bpytmp)

   2  CONTINUE





      if (outflag) then

           print*,'WARNING:  One or more points were out of bounds.'

           print*,'          Check input data file e.g.,'

           print*,'          expansion factor, geographical limits.'

      endif

      RETURN

      END

C---------------------------------------------------------------------------

      SUBROUTINE PLOTXY(X,Y,NC)



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



      integer*2 numpts, ixtmp(2000), iytmp(2000)



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



      real xmin,xmax,ymin,ymax,scx,scy



      COMMON /DDSCALE/ XMIN,XMAX,YMIN,YMAX,MINX,MAXX,MINY,MAXY,SCX,SCY,

     1 NSCX,NXCY,MSCX,MSCY,ISCX,ISCY

      COMMON/DD/IN,IOR,IT,IS,IC,ICC,IX,IY

      INTEGER REZFLG

      PARAMETER (REZFLG=2001)

      COMMON/DDREZ/DDHIX,DDHIY  !!!!!



      character*12 lastcmd   !!!!!

      logical outflag



      SCX=X

      SCY=Y

      CALL DDSC

      IX=ISCX

      IY=ISCY



C.....DDBP has been eliminated due to the following logic.



      IF(NC.NE.0) then



           IF(IX.EQ.REZFLG) THEN

              IX=DDHIX

              IY=DDHIY

           ELSE

              DDHIX=IX

              DDHIY=IY

           ENDIF



           if (numpts .gt. 0) then

                WRITE (ndevtxt,444) lastcmd



                WRITE (ndevtxt,99) numpts

                nremain = MOD(numpts,8)

	        if (numpts .ge. 8) then

	             do ii = 1, numpts-7, 8

		          WRITE (ndevtxt,111)

     &			    (ixtmp(j),iytmp(j),j=ii,ii+7)

	             enddo

	        endif



C.....	    if NUMPTS was not evenly divisible by 8, then the

C.....      remaining points are written to the file one per line.



	        if (nremain.gt.0) then

	            WRITE (ndevtxt,222)

     &  	      (ixtmp(k),iytmp(k),k=numpts-nremain+1,numpts)

	        endif

           endif

           WRITE (lastcmd,333) 1,int(1.d1*ddhix),int(1.d1*ddhiy)

333        format (i2,2i5)

           numpts = 0

      endif



99    FORMAT (I5)

111   FORMAT (8(2I5))

222   FORMAT (2I5)

444   FORMAT (A)

      IF(NC.EQ.0) then

           CALL DDVC_ACUM (numpts, ixtmp, iytmp, outflag)

      endif



      RETURN

      END

C---------------------------------------------------------------------------

      SUBROUTINE PLTANH(X,Y,NC,PLOT)

C

C     TICK LABELING ROUTINE FOR HORIZONTAL and vertical PLOT PROJECTIONS

C

C.....Vertical axes only for contour plots.



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



      external plot

      CHARACTER ANNOT*10

      COMMON/DD/IN,IOR,IT,IS,IC,ICC,IX,IY

      character*80 id

      common /ident/ id

C     COMMON DECK "WWR" INSERTED HERE

      PARAMETER (NWARSZ=1000)

      COMMON/WW/dum(10),MAXW,W(NWARSZ)

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

C     COMMON DECK "CONST" INSERTED HERE

      COMMON/PCONST/CREF,RGAS,GAMMA

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

      COMMON/LABCLT/PROJCT,THMIN,THMAX,RMIN,RMAX

C

      COMMON/RAYCON/MCONP

C

      DATA LNC/-100/

C

      IF(NC.LE.0 .OR. NC.GT.2) GO TO 100

C

C     NORMALIZE LETTER SIZE FACTOR TO .15 INCHES

      HLETF=HITLET/.15

C

      IF(LNC.NE.NC) LICM=-100

      LNC=NC

C

      CALL PLOT(X,Y,1)

      IX=IX-55*HLETF      !!!!! changed 80 to 55

      ICM=IX

      IF(NC.GT.1) THEN

         ICM=IY

         IX=IX-40

      ENDIF

C

C     INSURE THAT OVERLAPS OF ANNOTATIONS DO NOT OCCUR

      IF(IABS(ICM-LICM).LT.INT(80*HLETF)) GO TO 100

C

      IF(PROJCT.EQ.2.0) GO TO 25

      IF(HB.GE.0.0.AND.Y.GT.(RMIN+RMAX)/2) GO TO 100

      IF(HB.LT.0.0.AND.Y.LT.(RMIN+RMAX)/2) GO TO 100

C

25    continue

      degtick = tic*degs*1.e+5

      rkmtick = tic*earthr*1000

      degtmp = MOD(degtick,1)

      rkmtmp = MOD(rkmtick,1)



      F=DEGS

      TMP=(THMAX-THMIN)*DEGS

      IF(TMP.LT.10.) F = EARTHR



      if (ABS(degtmp-NINT(degtmp)) .lt. ABS(rkmtmp-NINT(rkmtmp))) then

           f = degs

      else

           f = earthr

      endif

      FT=F



      IF(PROJCT.EQ.2.0 .AND. MCONP.EQ.0) F=F/EARTHR

      V=X-THMIN

      IF(NC.GT.1) V=Y-(RMIN+RMAX)/2.

      IF(ABS(AMOD(TIC*FT+.0001,1.)).GT..001) GO TO 60

C

C.....Shift label to right if it has only 1 or 2 places.



      if (INT(v*f+SIGN(.5,v)) .ge. 0 .and. INT(v*f+SIGN(.5,v)) .lt. 100)

     &     ix = ix + 5

      WRITE(ANNOT,50) INT(V*F+SIGN(.5,V))

50    FORMAT(I5)

      GO TO 90

C

60    WRITE(ANNOT,80) V*F

80    FORMAT(F8.2)

C

90    LICM=ICM

C     ALLOW 8 RASTERS FOR THE LETTER HEIGHT

      H=HB

      IF(PROJCT.EQ.2.0) H=0.

      IF(NC.EQ.1) IY=IY-8*HLETF-dSIGN(52.,H)

      IOR=0

95    FORMAT(3G13.6,2I5,1X,A10)

99    FORMAT(2G13.6,5I5,1X,A10)

      CALL DDTEXT(1,ANNOT)

C

100   CALL PLOT(X,Y,MIN0(1,NC))

      END

C-----------------------------------------------------------------------

      SUBROUTINE PLTHLB(X,Y,NC)

C     HORIZONTAL TICK ANNOTATION ROUTINE FOR RAYPLOT.

C

C

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



      EXTERNAL PLOT

C

      CALL PLTANH(X,Y,NC,PLOT)

      END

C---------------------------------------------------------------------------

      SUBROUTINE PLTLB(X,PY,NCP)

C     PUT VERTICAL TICK ANNOTATIONS ON RAY PLOT

C     THIS IS A SPECIAL PURPOSE SUBSTITUTE FOR THE 'PLOT' ROUTINE

C     IT COMPUTES THE NEAREST ROUNDED TIC POSITIONS FOR VERTICAL LINES.

C     WHEN THE NC PARAMETER IS > 0 AN ANNOTATION IS GENERATED

C     AT THE TIC POSITION. THE ADDITIONAL NC VALUE < 0 ALLOWS FOR

C     SIMPLE PEN UP MOTION TO THE X,Y POSITION.

C



C*****When this routine calls PLOT, it is calling subroutine "PLOT".

C     Unlike many other routines which have PLOT passed in as an argument

C     and can be one of several other routines.



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



      CHARACTER ANNOT*10

C

      character*80 id

      common /ident/ id

C     COMMON DECK "WWR" INSERTED HERE

      PARAMETER (NWARSZ=1000)

      COMMON/WW/dum(10),MAXW,W(NWARSZ)

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

      COMMON /DD/ IINT,IOR,IT,IS,IC,ICC,IX,IY

      COMMON/LABCLT/PROJCT,THMIN,THMAX,RMIN,RMAX

C

      DATA LIY/-100/

C

      NC=NCP

      GO TO 10

C

      ENTRY PLTLBN(X,PY,NCP)

      NC=0

C

10    Y=PY

C

      IF(NCP.LT.0) GO TO 200

      F=1.

      IF(AMAX1(ABS(RMAX-EARTHR),ABS(RMIN-EARTHR)).LT.1.) F=1000.

      Y=INT((Y-EARTHR)*F)/F+EARTHR

      IF(NC.EQ.0) GO TO 100

C

      V=Y

ccccc      IF(PROJCT.EQ.3.) V=RMIN+(Y-RMIN)/PFACTR  ! original



c.....M. Jones correction, May 1991



      if (projct.eq.3.) v=rmin+(y-rmin)/(pfactr*

     &                 (rmin/(earthr+(rmin-earthr)*pfactr)))



      IF(ABS(AINT(TICV*F)-TICV*F).GT..001) GO TO 60

C

      WRITE(ANNOT,50) INT(ABS(V-EARTHR)*F)

50    FORMAT(I3)

      GO TO 90

C

60    WRITE(ANNOT,80) ABS(V-EARTHR)*F

80    FORMAT(F6.2)

C

90    CALL PLOT(X,PY,1)

      CALL PLOT(X,Y,1)

C     INSURE THAT OVERLAPS OF ANNOTATIONS DO NOT OCCUR

      IF(IABS(IY-LIY).LT.80) GO TO 100

      LIY=IY

      IX=IX-100

      IOR=0

      CALL DDTEXT(1,ANNOT)

C

100   CALL PLOT(X,Y,NC)

      RETURN

C

200   CALL PLOT(X,Y,1)

      RETURN

      END



C------------------------------------------------------------------------------

      SUBROUTINE SETUP

C     OBTAIN DATE-TIME, READ USER IDENTIFICATION AND COMPUTE CONSTANTS.

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



C     COMMON DECK "RAYDEV" INSERTED HERE

C     DEVICE ASSIGNED TO RAYTRC INPUT FILE

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

C

C     COMMON DECK "CONST" INSERTED HERE

      COMMON/PCONST/CREF,RGAS,GAMMA

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

C     COMMON DECK "HDR" INSERTED HERE

      CHARACTER*10 INITID*80,DAT,TOD

      COMMON/HDR/SEC

      COMMON/HDRC/INITID,DAT,TOD

      PARAMETER (NWARSZ=1000)

      COMMON/WW/dum(10),MAXW,W(NWARSZ)

C

ccccc      CALL SYSDAT(DAT)

ccccc      CALL SYSTIM(TOD)

ccccc      CALL SYSSEC(SEC)

C     READ IN USERS NAME AND TELEPHONE EXTENSION FOR IDENTIFYING MICROF

C     PLOTS.

ccccc      READ(NRYIND,1100,END=10000) INITID

ccccc1100  FORMAT(A)

ccccc10000 CALL DDINIT(8,INITID)



CONPLT.21

C********* INITIALIZE MATH CONSTANTS AND THE W ARRAY

      CALL SETW

      END

C------------------------------------------------------------------------------

      SUBROUTINE TIKLINE(XL1,YB,XL,YT1,TICV,TIKSZ,PLOT)



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



C     DRAWS STRAIGHT LINE WITH TICKS AT INTERVALS. Usually the vertical axis.

C     Called for projections 1,3,4.

C

      XDF=XL-XL1

      YDF=YT1-YB

      DST=SQRT(XDF*XDF+YDF*YDF)

      TICE=DST

      T=0.0

      IF(TICV.EQ.0.0) GO TO 50

      TICE=TICV

      T=TIKSZ/TICV

C

50    TICVX=TICE*XDF/DST

      TICVY=TICE*YDF/DST

      TX=TICVY*T

      TY=-TICVX*T

      NTIC=1+DST/TICE

      CALL PLOT(XL1,YB,-1)

      DO 100 I=0,NTIC-1

           X=XL1+I*TICVX

           Y=YB+I*TICVY

           CALL PLOT(X,Y,0)

           CALL PLOT(X+TX,Y+TY,0)

100        CALL PLOT(X,Y,1)

      CALL PLOT(XL,YT1,0)

      CALL PLOT(XL,YT1,1)

      END



      LOGICAL FUNCTION UCON(JC1,U,CONV)



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



      LOGICAL SETUCON

C     PROVIDES KEYWORD UNITS CONVERSION FOR W-ARRAY INPUT.

C     UNITS SPECIFICATION MUST BE IN THE FORM 'DV UN' WHERE DV IS THE

C     TYPE OF VARIABLE(SUCH AS AN FOR ANGLE) AND UN IS THE UNITS CHOICE

C     SUCH AS RD FOR RADIANS). RETURN VALUE IS THE CONVERSION FACTOR.

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

      PARAMETER (NWARSZ=1000)

      COMMON/WW/dum(10),MAXW,W(NWARSZ)

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

      CHARACTER U*(*),PU*3,PV*2

C

      DATA NPV,NPU/5,4/

      UCON=.FALSE.

      CONV=1.0

C

      READ(U,100) PU,PV

100   FORMAT(A,A)

C

      DO 2010 J=1,NPV

2010     IF(PU.EQ.PCV(J)) GO TO 2015

      JC1=4

      RETURN

2015  JC1=J

      DO 2020 K=1,NPU

2020     IF(PV.EQ.CNVC(J,K)) GO TO 2025

      JC1=4

      RETURN

2025  UCON=.TRUE.

      CONV=CNVV(JC1,K)

      IF(CONV.EQ.-1.0) CONV=1.0/EARTHR

      RETURN

C

C     INITIAL CONVERSION CONSTANTS

C-----------------------------------------------------------------------------

      ENTRY SETUCON !!!!!(JC1,U,CONV)

      SETUCON=.TRUE.

C

      CNVV(1,2) = RAD

C     USE TAG -1 FOR DYNAMIC USE OF EARTHR

      CNVV(1,3) = -1.0

      CNVV(2,2) = 1.0E-3

      CNVV(2,3) = 3.048006096E-4

      CNVV(2,4) = 1.852

      CNVV(3,2) = PIT2

      CNVV(3,3) = PIT2*1.0E3

      CNVV(3,4) = -PIT2

      CNVV(4,1) = 1.0/4.3429448

      END

c-------------------------------------------------------------------------

      SUBROUTINE IZPZ(THETL,PHIL,THETR,PHIR)

C

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



      DOUBLE PRECISION KPX,KPY,KPZ,IPX,IPY,IPZ

      DOUBLE PRECISION LX,LY,LZ

      DOUBLE PRECISION L(3),R(3),KP(3),IP(3),V(3),RS(3)

C

      COMMON/ZPZC/LX,LY,LZ,RX,RY,RZ,KPX,KPY,KPZ,IPX,IPY,IPZ

      COMMON/ZPZC/VX,VY,VZ,RSX,RSY,RSZ

      COMMON/ZPZC/BETA,A(3,3)

C

      EQUIVALENCE (L,LX),(R,RX),(KP,KPX),(IP,IPX),(V,VX)

      EQUIVALENCE (RS,RSX)

C

C

      S=SIN(THETL)

      LX=S*COS(PHIL)

      LY=S*SIN(PHIL)

      LZ=COS(THETL)

C

      S=SIN(THETR)

      RX=S*COS(PHIR)

      RY=S*SIN(PHIR)

      RZ=COS(THETR)

C

C     DETERMINE PRIMED Z AXIS BY CROSS PRODUCT OF L AND R

      CALL CROSSP(L,R,KP)

      CALL NORMV(KP)

      IF(ABS(ABS(KPZ)-1.0).GT.1.E-10) GO TO 10

C

C     ONLY A SIMPLE Z AXIS ROTATION IS REQUIRED, SKIP MATRIX COMPUTATION

C     MAKE IT A CLEAN 90 DEGREES, MINUS ANGLES FOR SOUTH POLE CASE

C

      KPZ=SIGN(1.0,KPZ)

      BETA=PHIL*KPZ

      RETURN

C

10    CALL CROSSP(KP,L,IP)

      CALL COPYVL(L,A)

      CALL COPYVL(IP,A(2,1))

      CALL COPYVL(KP,A(3,1))

C

      END

      SUBROUTINE ZPZ(THETAP,PHIP,THETA,PHI)

C

C     'ZPZ' CONVERTS THE SPHERICAL COORDINATES THETAP,PHIP

C     OF A POINT IN PRIMED SPACE TO COORDINATES THETA,PHI IN

C     UNPRIMED SPACE WHERE PRIMED SPACE IS LOCATED RELATIVE TO

C     UNPRIMED SPACE SUCH THAT POINTS THETL,PHIL AND THETR,PHIR

C     LIE ON THE EQUATOR. SEE IZPZ FOR INITIALIZATION. METHOD

C     USED IS TO CONVERT TO CARTESIAN COORDINATES, ROTATE TO UNPRIMED

C     SYSTEM AND CONVERT BACK TO SPHERICAL COORDINATES.

C

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

C

C     COMMON DECK "CONST" INSERTED HERE

      COMMON/PCONST/CREF,RGAS,GAMMA

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

C

      DOUBLE PRECISION KPX,KPY,KPZ,IPX,IPY,IPZ

      DOUBLE PRECISION LX,LY,LZ

C

      COMMON/ZPZC/LX,LY,LZ,RX,RY,RZ,KPX,KPY,KPZ,IPX,IPY,IPZ

      COMMON/ZPZC/VX,VY,VZ,RSX,RSY,RSZ

      COMMON/ZPZC/BETA,A(3,3)

C

      IF(ABS(KPZ).NE.1.0) GO TO 100

      THETA=THETAP

      IF(KPZ.LT.0.0) THETA=PI-THETAP

      PHI=PHIP+BETA

      RETURN

C

100   S=1.0

      IF(THETAP.NE.PID2) S=SIN(THETAP)

      VX=S*COS(PHIP)

      VY=S*SIN(PHIP)

      VZ=0.0

      IF(THETAP.NE.PID2) VZ=COS(THETAP)

C

      CALL MATRIXM(A,VX,RSX)

      THETA=ACOS(RSZ)

      PHI=ATAN2(RSY,RSX)

      IF(ABS(RSZ).EQ.1.0) PHI=0.0

      END

C------------------------------------------------------------------------------

      SUBROUTINE DDINIT(N,TEXT)

C----------------------------------------------------------------------

C     INITIALIZES PLOTTING PROCESS(DDPLOT)

C     WRITE PARAMETERS TO GRAPHICS CALLS FILE

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,ndevtxt,NDEVBIN

      CHARACTER TEXT*(*)

C

      DATA IDD/0/

10    FORMAT(5A10)

      IF(IDD.EQ.0) REWIND ndevtxt

      IDD=1

c%      WRITE(NDEVGRP) 0,0,0

c%      WRITE(NDEVGRP) SNGL(M),LEN(TEXT),TEXT

      write (ndevtxt,33) 0,n,len(text)

      write (ndevtxt,34) text

33    format (i2,2i3)

34    format (a80)

      END

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

      SUBROUTINE DDBP

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

C     SETS A VECTOR ORIGIN(DDPLOT)

C     WRITE PARAMETERS TO GRAPHICS CALLS FILE

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,ndevtxt,NDEVBIN

      COMMON/DD/IN,IOR,IT,IS,IC,ICC,IX,IY

      INTEGER REZFLG

      PARAMETER (REZFLG=2001)

      COMMON/DDREZ/DDHIX,DDHIY

      COMMON/KNKN/KNBP,KNVC,KNDT

      COMMON/DDLIM/MXIX,MXIY,MNIX,MNIY

C

C     USE PRE-ASSIGNED FLAG VALUE FOR IX TO DETERMINE USE OF

C     FLOATING POINT RASTER VALUES OF X AND Y USED IN DDBP AND DDVC

C

      IF(IX.EQ.REZFLG) THEN

         IX=DDHIX

         IY=DDHIY

      ELSE

         DDHIX=IX

         DDHIY=IY

      ENDIF

C

      MNIX=MIN0(MNIX,IX)

      MXIX=MAX0(MXIX,IX)

      MNIY=MIN0(MNIY,IY)

      MXIY=MAX0(MXIY,IY)

      IF(IX.LT.0.OR.IX.GT.1023.OR.IY.LT.0.OR.IY.GT.1023

     1  ) PRINT *,'DDBP knbp,ix,iy,ddhix,ddhiy',KNBP,IX,IY,ddhix,ddhiy

ccccc10    FORMAT(A10,3I5)

      KNBP=KNBP+1

      write (ndevtxt,33)1,int(1.d1*ddhix),int(1.d1*ddhiy)

33    format (i2,2i5)

      END

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

      SUBROUTINE DDVC

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

C     PLOTS A VECTOR(DDPLOT)

C     WRITE PARAMETERS TO GRAPHICS CALLS FILE

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,ndevtxt,NDEVBIN

      COMMON/DD/IN,IOR,IT,IS,IC,ICC,IX,IY

      INTEGER REZFLG

      PARAMETER (REZFLG=2001)

      COMMON/DDREZ/DDHIX,DDHIY

      COMMON/KNKN/KNBP,KNVC,KNDT

      COMMON/DDLIM/MXIX,MXIY,MNIX,MNIY



      integer*2 numpts, ixtmp(2000), iytmp(2000)



C.....Set OUTFLAG to true when a plotting coordinate is out of bounds

C     This can happen if the DINP file contains a large expansion factor

C     or a large geographical area.



      logical outflag

C

C     USE PRE-ASSIGNED FLAG VALUE FOR IX TO DETERMINE USE OF

C     FLOATING POINT RASTER VALUES OF X AND Y USED IN DDBP AND DDVC

C

      entry DDVC_ACUM (numpts, ixtmp, iytmp, outflag)

      IF(IX.EQ.REZFLG) THEN

         IX=DDHIX

         IY=DDHIY

      ELSE

         DDHIX=IX

         DDHIY=IY

      ENDIF

C

      MNIX=MIN0(MNIX,IX)

      MXIX=MAX0(MXIX,IX)

      MNIY=MIN0(MNIY,IY)

      MXIY=MAX0(MXIY,IY)

      IF(IX.LT.0.OR.IX.GT.1023.OR.IY.LT.0.OR.IY.GT.1023

     1  ) outflag = .true.  !!!!! PRINT *,'DDVC ',KNVC,IX,IY

ccccc10    FORMAT(A10,3I5)

      KNVC=KNVC+1

c%      WRITE(NDEVGRP) 2,SNGL(DDHIX),SNGL(DDHIY)

c      write(ndevtxt,33)2,int(1.d1*ddhix),int(1.d1*ddhiy)

c33    format (i2/2i5)

      numpts = numpts + 1



      ixtmp(numpts) = int(1.d1*ddhix)

      iytmp(numpts) = int(1.d1*ddhiy)



      END

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

      SUBROUTINE DDTEXT(N,TEXT)

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

C     WRITES A CHARACTER ARRAY PACKED A10(DDPLOT)

C     WRITE PARAMETERS TO GRAPHICS CALLS FILE

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,ndevtxt,NDEVBIN

      CHARACTER TEXT*(*)

      COMMON/DD/IN,IOR,IT,IS,IC,ICC,IX,IY

      COMMON/KNKN/KNBP,KNVC,KNDT

100   FORMAT(A10,2I5,I2,1X,8A10)

      KNDT=KNDT+1

c%      WRITE(NDEVGRP) 3,IX,IY

      write (ndevtxt,33)3,ix,iy

33    format (i2,2i4)

      LN=MIN0(N*10,LEN(TEXT))

      IF(LN.EQ.N*10) THEN

           write (ndevtxt,34)ior,n*10,(text(i:i),i=1,len(text))

34         format (2i3/80a1)

      ELSE

           write (ndevtxt,35)ior,n*10,(text(i:i),i=1,len(text))

     1     ,(' ',i=ln+1,n*10)

35         format (2i3/80a1/80a1)

      ENDIF

      END

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

      SUBROUTINE DDFR

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

C     ADVANCE ONE PLOTTING FRAME(DDPLOT)

C     WRITE PARAMETERS TO GRAPHICS CALLS FILE

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,ndevtxt,NDEVBIN

10    FORMAT(A10)

c%      WRITE(NDEVGRP) -2,0,0

      write (ndevtxt,33)-2

33    format (i2)

      END

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

      SUBROUTINE DDEND

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

C     EMPTIES PLOT BUFFER AND RELEASES PLOTTING COMMAND FILE(DDPLOT)

C     WRITE PARAMETERS TO GRAPHICS CALLS FILE

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,ndevtxt,NDEVBIN

      COMMON/KNKN/KNBP,KNVC,KNDT

      COMMON/DDLIM/MXIX,MXIY,MNIX,MNIY

10    FORMAT(3A10,5I5)

c%      WRITE(NDEVGRP) -1,0,0

      write (ndevtxt,33)-1

33    format (i2)

      END

C----------------------------------------------------------------------------

      SUBROUTINE CLEAR(A,N)

C     SET N ELEMENTS OF ARRAY A TO ZERO

      DOUBLE PRECISION A(N)

      IF(N.LE.0) RETURN

      DO    10 I=1,N

10    A(I)=0.0

      RETURN

      END

      SUBROUTINE MX1ALF(T1,T2)

C     SPECIFY USE OF ALTERNATE CHARACTER SET NUMBER 1(DISSPLA)

C     WRITE PARAMETERS TO GRAPHICS CALLS FILE

C

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

      CHARACTER*10 T1,T2



      WRITE(NDEVTXT,100) 11, T1, T2

100   FORMAT(I2,TR1,A,1X,A1)

      END

C----------------------------------------------------------------------------

      SUBROUTINE MX2ALF(T1,T2)

C     SPECIFY USE OF ALTERNATE CHARACTER SET NUMBER 2(DISSPLA)

C     WRITE PARAMETERS TO GRAPHICS CALLS FILE

C

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

      CHARACTER*10 T1,T2



      WRITE(NDEVTXT,200) 12,T1,T2

200   FORMAT(I2,TR1,A,1X,A1)

      END



      SUBROUTINE OPNREP(IUN,FNAME)

C     THIS OPEN OPERATION ALLOWS FOR AN EXISTING VERSION OF

C     A FILE. IF THE FILE EXISTS IT IS OVERWRITTEN.

C     THE LOGIC USED IS HOPEFULLY MACHINE INDEPENDENT.

      CHARACTER*(*) FNAME

      CALL ZAPFIL(IUN,FNAME)

      OPEN(IUN,FILE=FNAME,STATUS='NEW')

      REWIND IUN

      RETURN

      ENTRY OPNURP(IUN,FNAME)

      CALL ZAPFIL(IUN,FNAME)

      OPEN(IUN,FILE=FNAME,STATUS='NEW',FORM='FORMATTED')

      REWIND IUN

      END

C----------------------------------------------------------------------------

      SUBROUTINE ZAPFIL(IUN,FNAME)

C     TEST FOR EXISTENCE OF FILE AND DELETE

      CHARACTER FNAME*(*)

      LOGICAL EX

      INQUIRE(FILE=FNAME,EXIST=EX)

      IF(EX) THEN

         OPEN(UNIT=IUN,FILE=FNAME,STATUS='OLD')

         CLOSE(UNIT=IUN,STATUS='DELETE')

      ENDIF

      END

C----------------------------------------------------------------------------

      SUBROUTINE HEIGHT(H)

C     SETS REFERENCE CHARACTER HEIGHT IN INCHES(DISSPLA)

C     WRITE PARAMETERS TO GRAPHICS CALLS FILE

C

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



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

C

      WRITE(NDEVTXT,100) 13,H

100   FORMAT(I2,F4.2)

      END

C----------------------------------------------------------------------------

      FUNCTION ITEST(I)

C     USED TO PASS INTEGER VALUES THROUGH FOR VARIABLES TYPED REAL

C     (AS IN MIXED REAL/INTEGER ARRAYS)

      ITEST=I

      END

C------------------------------------------------------------------------

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

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

      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

               PRINT*,'Format error in units specifier of tabular data.'

               STOP

            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

               PRINT*,'Format error in units specifier of tabular data.'

               STOP

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

           PRINT*,'Format error in units specifier of tabular data.'

           STOP

      endif

      END

       SUBROUTINE DDSC

C***E     DDSC

C***.

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



      INTEGER REZFLG

      PARAMETER (REZFLG=2001)

      COMMON/DDREZ/DDHIX,DDHIY

      DOUBLE PRECISION DDHI(2),ISC

      EQUIVALENCE (DDHIX,DDHI)

      real xm,scx, xmin,xmax

      COMMON /DDSCALE /XM(4), MX(4), SCX(2), NSCX(2), MSCX(2), ISCX(2)

      DIMENSION XMIN (3), XMAX (3), MINX (3), MAXX (3)

      EQUIVALENCE (XMIN, XM(1)), (XMAX, XM (2)), (MINX, MX(1)),

     &            (MAXX, MX (2))

      DATA XM,MX/4*0.,4*0/

      CALL DDST

      L2 = 1

      DO 3 L1 = 1, 2

      MD = MAXX (L2) - MINX (L2)

      Q = MD * (SCX (L1) - XMIN (L2)) / (XMAX (L2) - XMIN (L2))

      ISC = Q + SIGN (.5, Q)

C      ISCX (L1) = ISC + MINX (L2)

      DDHI(L1)=ISC+MINX(L2)

      IF (ISC .LT. 0) GO TO 1

      IF (ISC .GT. MD) GO TO 1

      NSCX (L1) = 0

      MSCX (L1) = ISCX (L1)

      GO TO 3

   1  NSCX (L1) = ISC / MD

      MSCX (L1) = ISC - MD * NSCX (L1)

      IF (MSCX (L1) .GE. 0) GO TO 2

      NSCX (L1) = NSCX (L1) - 1

      MSCX (L1) = MSCX (L1) + MD

   2  MSCX (L1) = MSCX (L1) + MINX (L2)

   3  L2 = 3

      ISCX(1)=REZFLG

      RETURN

      END

C---------------------------------------------------------------------------

      SUBROUTINE DDST

C     CALLED BY DDGRAPH TO COMPUTE CURVILINEAR PROJECTIONS AT SCALED

C     POINT SCX,SCY. ENTRY SETDDST COMPUTES SCALING CONSTANTS AND BOUND

C     OF PLOT REGION.

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



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

      PARAMETER (NWARSZ=1000)

      COMMON/WW/dum(10),MAXW,W(NWARSZ)

      EQUIVALENCE (EARTHR,W(1)),(PFACTR,W(82)),(PLT,W(81))



      real xmin,xmax,ymin,ymax,scx,scy



      COMMON /DDSCALE/ XMIN,XMAX,YMIN,YMAX,MINX,MAXX,MINY,MAXY,SCX,SCY,

     1 NSCX,NXCY,MSCX,MSCY,ISCX,ISCY

C

C

      COSC(X,Y)=SIN(Y)*STH0*COS(X-PH0)+COS(Y)*CTH0

      SCFLIP(X)=SQRT(AMAX1(0.,1.-X*X))

      SINB(X,Y)=SIN(X-PH0)*SIN(Y)/SCFLIP(COSCC)

      SINBC(X)=SIGN(SCFLIP(COSBC),THETA0-X)

C

      IF(PLT.EQ.4.) RETURN



      IF(PLT.EQ.2.) GO TO 10

      THETA=SCX

      R=SCY

C.....Correction M. Jones, March,1991



ccccc      IF(PLT.EQ.3.) R=(SCY-earthr)*PFACTR+earthr



ccccc      IF(PLT.EQ.3.) R=(SCY-RMIN)*PFACTR + RMIN  ! original



c.....newest correction, May 1991



      if (plt.eq.3.) then

           r = rmin +

     &         (scy-rmin)*pfactr*(rmin/(earthr+(rmin-earthr)*pfactr))

      endif



       SCX=R*SIN(THETA-THETA0)

       SCY=R*COS(THETA-THETA0)

ccccc      if (plt .eq. 3.) scy = scy + (rmin-earthr)*(1.-pfactr)  march



      RETURN

C

10    SCY=PID2-SCY

      COSCC=COSC(SCX,SCY)

      COSBC=1.0

C INSURE AGAINST DIVIDE BY ZERO

      IF(ABS(COSCC-1.0).GT.1.E-14) COSBC=SINB(SCX,SCY)

      C=ACOS(COSCC)

      SCX=C*COSBC

      SCY=C*SINBC(SCY)

      RETURN

C

      ENTRY SETDDST

C

      IF(PLT.LT.4.) GO TO 15

C     PLT = 4

      YMAX=YMIN+(XMAX-XMIN)*EARTHR/PFACTR

      RETURN

C

15    IF(PLT.EQ.2.) GO TO 20

C     PLT = 1 OR 3

      RMIN=YMIN

      THETA0=.5*(XMIN+XMAX)

      ALPHA=.5*(XMAX-XMIN)

      YMIN=YMIN*COS(ALPHA)

      XMAX=AMAX1(YMAX*SIN(ALPHA),(YMAX-YMIN)/2.)

      XMIN=-XMAX

      YMAX=2.*XMAX

      IF(PLT.EQ.1.) YMAX=YMAX/PFACTR

      YMAX=YMIN+YMAX



      RETURN

C

20    YMIN=PID2-YMIN

      YMAX=PID2-YMAX

      THETA0=.5*(YMIN+YMAX)

      PH0=.5*(XMIN+XMAX)

      STH0=SIN(THETA0)

      CTH0=COS(THETA0)

      COSCC=COSC(XMIN,YMIN)

      COSBC=SINB(XMIN,YMIN)

      C=ACOS(COSCC)

      XMIN=C*COSBC

      YMIN=C*SINBC(YMIN)

      COSCC=COSC(XMAX,YMAX)

      COSBC=SINB(XMAX,YMAX)

      C=ACOS(COSCC)

      XMAX=C*COSBC

      YMAX=C*SINBC(YMAX)

      END

C----------------------------------------------------------------------------

      SUBROUTINE DFCNST(DRMACH,MAX)

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



      DOUBLE PRECISION DRMACH(MAX)

C     RETURN MAX MACHINE CONSTANTS(OR THE NUMBER AVAILABLE)

C     IN ARRAY DRMACH.



C     THIS MODULE PROVIDES CONSTANTS WHICH ARE EXPLICITLY DEPENDENT

C     ON MACHINE ARCHITECTURE AND SHOULD BE EXAMINED WHEN MAKING ANY

C     MIGRATION OF THE RAY TRACING ROUTINES.

C

C

      DOUBLE PRECISION RMACH(5)

C

C     SINGLE PRECISION MACHINE CONSTANTS

C      (SEE DIGITAL SIGNAL PROCESSING, IEEE PRESS 1979 P S-7)

C

C     THE CONSTANTS ARE DERIVED FROM THE FOLLOWING

C     FORM FOR FLOATING POINT NUMBERS:

C

C        SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )

C     FOR MOST MACHINES THE BASE B=2.

C     FOR THE CYBER 60-BIT MACHINES

C        T=48, EMIN=-974 AND EMAX=1070.

C     FOR THE PDP-11 AND VAX-11 MACHINES

C        T=24, EMIN=-127 AND EMAX=127

C

C     RMACH(1)  SMALLEST POSITIVE MAGNITUDE = B**(EMIN-1)

C     RMACH(2)  LARGEST MAGNITUDE           = B**EMAX*(1-B**(-T))

C     RMACH(3)  SMALLEST RELATIVE SPACING   = B**(-T)

C     RMACH(4)  LARGEST RELATIVE SPACING    = B**(1-T)

C     RMACH(5)                              = LOG10(B=2)

C

C     CONSTANTS HERE ARE FOR CDC 6000/7000/CYBER SERIES COMPUTERS

C     PLEASE SUBSTITUTE FOR OTHER MACHINES

C

C

C      INTEGER SMALL(2)

C      INTEGER LARGE(2)

C      INTEGER RIGHT(2)

C      INTEGER DIVER(2)

C      INTEGER LOG10(2)

CC

C      EQUIVALENCE (RMACH(1),SMALL(1))

C      EQUIVALENCE (RMACH(2),LARGE(1))

C      EQUIVALENCE (RMACH(3),RIGHT(1))

C      EQUIVALENCE (RMACH(4),DIVER(1))

C      EQUIVALENCE (RMACH(5),LOG10(1))

C

C     MACHINE CONSTANTS FOR CYBER 170/180 COMPUTERS

C

ccccc      DATA RMACH(1) / O"00014000000000000000" /

ccccc      DATA RMACH(2) / O"37767777777777777777" /

ccccc      DATA RMACH(3) / O"16404000000000000000" /

ccccc      DATA RMACH(4) / O"16414000000000000000" /

ccccc      DATA RMACH(5) / O"17164642023241175720" /

C

C***REFERENCES  FOX, P.A., HALL, A.D., SCHRYER, N.L, *FRAMEWORK FOR

C                 A PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHE-

C                 MATICAL SOFTWARE, VOL. 4, NO. 2, JUNE 1978,

C                 PP. 177-188.

C***ROUTINES CALLED  XERROR

C***END PROLOGUE  R1MACH

C     MACHINE CONSTANTS FOR VAX 11/780

C     (EXPRESSED IN INTEGER AND HEXADECIMAL)

C    ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS***

C    *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***

C

C     DATA SMALL(1), SMALL(2) /        128,           0 /

C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /

C     DATA RIGHT(1), RIGHT(2) /       9344,           0 /

C     DATA DIVER(1), DIVER(2) /       9472,           0 /

C     DATA LOG10(1), LOG10(2) /  546979738,  -805796613 /

C

C     DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 /

C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /

C     DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 /

C     DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 /

C     DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFF884FB /

C

C   MACHINE CONSTANTS FOR VAX 11/780 (G-FLOATING)

C     (EXPRESSED IN INTEGER AND HEXADECIMAL)

C    ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS***

C    *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***

C

C     DATA SMALL(1), SMALL(2) /         16,           0 /

C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /

C     DATA RIGHT(1), RIGHT(2) /      15552,           0 /

C     DATA DIVER(1), DIVER(2) /      15568,           0 /

C     DATA LOG10(1), LOG10(2) /  1142112243, 2046775455 /

C

C     DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 /

C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /

C     DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 /

C     DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 /

C     DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F /

C

C     MACHINE CONSTANTS FOR THE CRAY 1

C

C     DATA SMALL(1) / 201354000000000000000B /

C     DATA SMALL(2) / 000000000000000000000B /

C

C     DATA LARGE(1) / 577767777777777777777B /

C     DATA LARGE(2) / 000007777777777777774B /

C

C     DATA RIGHT(1) / 376434000000000000000B /

C     DATA RIGHT(2) / 000000000000000000000B /

C

C     DATA DIVER(1) / 376444000000000000000B /

C     DATA DIVER(2) / 000000000000000000000B /

C

C     DATA LOG10(1) / 377774642023241175717B /

C     DATA LOG10(2) / 000007571421742254654B /

C

ccccc      CALL RMOVE(DRMACH,RMACH,MIN0(MAX,5))

C

      END

C----------------------------------------------------------------------

      SUBROUTINE DRAWTKS(DTICH,DTICV,XLPR,XR,YBPR,YT,PLOT)

C

C     HORIZONTAL PLOT PROJECTION (W(81) = 2) SUPPORT ROUTINE.

C     DRAW BOUNDARY of PLOT AREA, ticks AND tick LABELS.

C

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



      EXTERNAL PLOT

C

      XL=XLPR

      YB=YBPR

      YMID=.5*(YB+YT)

      TICX=.01*(YT-YB)

      TICY=.01*(XR-XL)

      IF(DTICV.GT.0.) THEN

         YBP=YMID-dINT((YMID-YB)/DTICV)*DTICV

         NTICX=(XR-XL)/DTICH

      ELSE

         YBP=YB

         NTICX=1

      ENDIF

      NTICY=1

      IF(DTICV.GT.0.) then

           NTICY = (YT-YBP+DTICV)/DTICV

      endif

C&     PRINT *,"DRAWTKS ",XLPR,XR,YBPR,YBP,YT,DTICH,DTICV,NTICX,NTICY

      IOF=0

C

C     TWO PASSES: BOTTOM AND LEFT axes THEN TOP AND RIGHT axes

C

      DO 100 J=1,2



C..........Moveto the lower left corner of boundary.



           CALL PLOT(XL,YBPR,1)



C..........Make the linetos to form the boundary and lineto to form the

C          tick mark.  Then label the tick mark via PLTANH.



           DO 30 I=1,NTICY

                Y=YBP+(I-1)*DTICV

                CALL PLOT(XL,Y,0)

                CALL PLOT(XL+TICY,Y,0)

30              CALL PLTANH(XL,Y,IOF+2,PLOT)



C..........Make the vertical line from the LL corner to UL corner.

           CALL PLOT(XL,yT,0)  !!!!! was YT

C

C..........Moveto the LL (or UR when J=2) corner again.

           CALL PLOT(XLPR,YB,1)

C

           DO 40 I=1,NTICX

                X=XLPR+I*DTICH

                CALL PLOT(X,YB,0)

                CALL PLOT(X,YB+TICX,0)

                CALL PLTANH(X,YB,IOF+1,PLOT)  !!!!! was YB

40              CALL PLOT(X,YB,1)



C..........Make the horizontal line to close the distance from the last

C          tick mark to the right boundary.



                CALL PLOT(XR,YB,0)



C..........Now change things to do the right Y axis and the top X axis.



           YB=YT

           XL=XR

           TICX=-TICX

           TICY=-TICY

100        IOF=10

C

C.....Make the horizontal line across the middle of the plot.



      CALL PLOT(XLPR,YMID,1)

      CALL PLOT(XR,YMID,0)

      RETURN

      END

C--------------------------------------------------------------------------

      SUBROUTINE FULANN

C     ANNOTATION MODEL SUITED FOR PUBLICATION QUALITY LETTERING

C

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



      CHARACTER*(*) S,C

C

C     INITIALIZES PLOT IN PUBLICATION-QUALITY MODE(REQUIRES DISSPLA)

C

C     COMMON DECK "WWR" INSERTED HERE

      PARAMETER (NWARSZ=1000) 

      COMMON/WW/ dum(10), MAXW,W(NWARSZ)

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

C     COMMON DECK "ANNOT" INSERTED HERE 

ccccc      CHARACTER*10 ANOTES,HNOTES

      character*13 anotes

      character*32 hnotes

      COMMON/ANNCTL/LENA(4),LENHA(4)

      COMMON/ANNCTC/ANOTES(4),HNOTES(4)

C

      INTEGER LHNA(4) ,LHNHA(4)

ccccc      CHARACTER*10 AQOTES(2,4),HQOTES(4,3)

      character*13 aqotes(4)

      character*32 hqotes(4)

C

ccccc      DATA LHNA,AQOTES/4*2,'DEPTH (#M!',')','DEPTH (#KM','!)',

ccccc     2'HEIGHT (#M','!)','HEIGHT (#K','M!)'/

      DATA LHNA,AQOTES/4*2,

     &                  'DEPTH (#M!)  ',

     &                  'DEPTH (#KM!) ',

     &                  'HEIGHT (#M!) ',

     &                  'HEIGHT (#KM!)'/



ccccc      DATA LHNHA,HQOTES/3,2,4, 'RANGE AT SEA LEVEL (#KM!)',' '

ccccc     2       ,'RANGE (#DE','G!)',2*' '

ccccc     3       ,'CROSS RANG','E AT SEA L','EVEL (#KM!',')'/

      DATA LHNHA,HQOTES/3,2,4,4,

     &                  'RANGE AT SEA LEVEL (#KM!)      ',

     &                  'RANGE (#DEG!)                  ',

     &                  'CROSS RANGE AT SEA LEVEL (#KM!)',

     &                  'CROSS RANGE AT SEA LEVEL (#DEG!)'/



C

      ENTRY SETANN

C

      CALL SCMPLX

C

C*****  PUT TRAILING BLANKS ON THE ARGUMENTS SO THAT NO EXTRANEOUS

C       CHARACTERS MAY BE INTRODUCED.  JAH.

C       THE ARGUMENTS ARE CHARACTER*10 IN MX1ALF,MX2ALF.  ESSENTIAL

C       FOR THE CYBER OPERATION.

C

      CALL MX1ALF('STAND     ','!')

      CALL MX2ALF('L/CSTD    ','#')

      CALL HEIGHT(HITLET)

C

      CALL IMOVE(LENA,LHNA,4) 

      CALL IMOVE(LENHA,LHNHA,4)  !!!! added a 4th element

ccccc DO 10 I=1,2

         DO 10 J=1,4

10       ANOTES(J)=AQOTES(J)

      DO 20 I=1,4  !!!!! added a 4th element

ccccc    DO 20 J=1,3

20       HNOTES(i)=HQOTES(i)

C

      RETURN

C



      ENTRY ANNFIL(S,C)

      C=S

      END

C-------------------------------------------------------------------------------

      SUBROUTINE SCMPLX

C     SPECIFY USE OF SIMPLEX CHARACTER SET(DISSPLA)

C     WRITE PARAMETERS TO GRAPHICS CALLS FILE

C

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

C

      WRITE(NDEVTXT,100) 10

100   FORMAT(I2)

      END

C----------------------------------------------------------------------------

      SUBROUTINE SET2(A,V,N)

C     SETS N COMPONENTS OF VECTOR 'A' TO SINGLE VALUE V

C



      DOUBLE PRECISION a(n),v

      integer ia(n)



C

      DO 100 I=1,N

100	  A(I)=V

      RETURN



C----------------------------------------------------------------------------

      ENTRY ISET2 (IA,IV,N)

C

C     SPECIAL ENTRY FOR INTEGER ARRAYS FOR MACHINES HAVING DIFFERENT

C     WORD SIZES FOR INTEGERS THAN DOUBLE PRECISION.

C



      DO 200 I=1,N

200	 IA(I) = IV

      RETURN

      END

C----------------------------------------------------------------------------

      SUBROUTINE SETXY(PROJ,XMIN,YMIN,XMAX,YMAX)

C     PLOT INITIALIZATION; SETS PROJECTION PARAMETERS

C

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



      COMMON/LABCLT/PROJCT,THMIN,THMAX,RMIN,RMAX

C

C

C   INITIAL ANNOTATION MODEL

      CALL SETANN

C

      PROJCT=PROJ

      THMIN=XMIN

      RMIN=YMIN

      THMAX=XMAX

      RMAX=YMAX

      END

C----------------------------------------------------------------------------

      SUBROUTINE SFILL(STG,LN,C)

C     FILLS A STRING WITH N SPECIFIED CHARACTERS

      CHARACTER STG*(*),C*1

C

      DO 10 I=1,LN

10    STG(I:I)=C

      RETURN

      END

C-----------------------------------------------------------------------------

      SUBROUTINE SFILTR(C,S,KSET)

C     FILTERS EXTRANEOUS CHARACTERS FROM PLOT LABELS

      CHARACTER*(*) C,S,KSET

C

      LN=LEN(C)

      J=0 

      DO 10 I=1,LN

        IF(INDEX(KSET,C(I:I)).EQ.0) THEN

           J=J+1

           S(J:J)=C(I:I)

        ENDIF

10    CONTINUE

      RETURN

      END

C-----------------------------------------------------------------------

      SUBROUTINE SREAD1(AB,NWOK,NW)

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'.



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



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

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

           PRINT*,'ERROR: Invalid W-ARRAY Value.',NW

           PRINT*,' Value must be >= 2 and <= 62.'

           STOP

      endif

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

      ENTRY SETSRD

C***** CHANGED "ISET2" TO "SET2".  SEE SET2 SUBROUTINE.

C*****	CALL SET2(LIST,0,MXLIST)

C***** changed to a loop rather than writing an integer version of SET2

C***** SET2 is for DOUBLE PRECISION*8.



      do 10 ii = 1,mxlist

10	   list(ii) = 0



      PMX=PPX

      CALL IMOVE(PNTBL,PQTBL,10)

      CALL IMOVE(PITBL,PLTBL,10)

      CALL IMOVE(PFRMTBL,PIRMTBL,10)

      END

C---------------------------------------------------------------------------

      INTEGER FUNCTION STRIM(C)

C     DETERMINES POSITION OF LAST NONBLANK CHARACTER OF A STRING

      CHARACTER*(*) C

C

      L=LEN(C)

      DO 100 I=L,1,-1

100   IF(C(I:I) .NE. ' ') GO TO 200

      I=0 

200   STRIM=I

      END

C---------------------------------------------------------------------------

      SUBROUTINE IMOVE(JOUT,INPI,N)

C     THIS ROUTINE REPLACES THE ENTRY POINT THAT WAS PART OF DFSYS. JAH.

C     THIS IS MACHINE INDEPENDENT AND THUS DOES NOT NEED TO BE PART OF

C     OF DFSYS.



      INTEGER INPI(*), JOUT(*),N



      DO 10 I=1,N

10    JOUT(I) = INPI(I)

      END

C-----------------------------------------------------------------------

      SUBROUTINE RMOVE(ROUT,RINP,N)

C     THIS ROUTINE REPLACES THE ENTRY POINT THAT WAS PART OF DFSYS. JAH.

  

      DOUBLE PRECISION RINP(*), ROUT(*)

      INTEGER N

  

      DO 10 I=1,N

10    ROUT(I) = RINP(I)

      END 

C------------------------------------------------------------------------

      FUNCTION GCIRCLE(LLAT,LLON,RLAT,RLON)



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



      DOUBLE PRECISION LLAT,LLON

C

      CTHR=SIN(RLAT)

      STHR=COS(RLAT)

      CTHL=SIN(LLAT)

      STHL=COS(LLAT)

      CLR=CTHL*CTHR+STHL*STHR*COS(RLON-LLON)

      GCIRCLE=ACOS(CLR)

      END

C-----------------------------------------------------------------------------

      FUNCTION GET1(Z)

C     FUNCTIONALLY THE SAME AS 'GET' PROGRAM, SEE DOCUMENTATION THERE

C     NEEDED BECAUSE RECEIVER MODELS WILL CALL GET TO OBTAIN TERRAIN VA

C     BUT THEY THEMSELVES ARE CALLED VIA GET. HENCE HAVE RE-ENTRANCE

C     PROBLEM.

C

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



      DOUBLE PRECISION Z(*)

C     COMMON DECK "SS" INSERTED HERE

      DOUBLE PRECISION MODSURF

      COMMON/SS/ MODSURF(4)

      COMMON/SS/U,PUR,PURR,PURTH,PURPH

      COMMON/SS/PUTH,PUPH,PUTHTH,PUPHPH,PUTHPH,USELECT,UTIME

C     COMMON DECK "RKAM" INSERTED HERE

      DOUBLE PRECISION KR,KTH,KPH

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

C     COMMON DECK "FNDER" INSERTED HERE

      COMMON/FNDER/GY(20),NZ,NPZR,NPZRR,NPZRTH,NPZRPH,NPZTH

      COMMON/FNDER/NPZPH,NPZTHTH,NPZPHPH,NPZTHPH,NSELECT,NTIME

C     COMMON DECK "RR" INSERTED HERE

      DOUBLE PRECISION MODREC

      COMMON/RR/ MODREC(4)

      COMMON/RR/F,PFR,PFRR,PFRTH,PFRPH

      COMMON/RR/PFTH,PFPH,PFTHTH,PFPHPH,PFTHPH,FSELECT,FTIME

C     COMMON DECK "GG" INSERTED HERE

      DOUBLE PRECISION MODG

      COMMON/GG/MODG(4)

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

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

C     COMMON DECK "RKTIME" INSERTED HERE

      COMMON/CRKTIME/RKTIME

C

C     COMMON DECK "RMACH" INSERTED HERE

      COMMON/CRMACH/RMACH(5)

      EQUIVALENCE(RKTIME,IRKTIME)

C     COMMON DECK "WWR" INSERTED HERE

      PARAMETER (NWARSZ=1000)

      COMMON/WW/dum(10),MAXW,W(NWARSZ)

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

C

      COMMON/CGET/ZERO

C

      IENTRY=1

      GO TO 5

C

      ENTRY DOT1(Z)

      IENTRY=2

5     CONTINUE

      IF(ZERO.EQ.0.0) ZERO=EARTHR*RMACH(3)*2.0

C

      IF(ITEST(Z(NTIME)).EQ.IRKTIME) GO TO 10

C

      IF(Z(NSELECT).EQ.USELECT) CALL SURFACE

      IF(Z(NSELECT).EQ.FSELECT) CALL RECEVER

      IF(Z(NSELECT).EQ.GSELECT) CALL TOPOG

      Z(NTIME)=RKTIME

C     REMOVE MACHINE ROUND OFF NOISE FROM EXACT RECEIVER LOCATIONS

      IF(ABS(Z(NZ)).LE.ZERO) Z(NZ)=0.0

C

10    IF(IENTRY.NE.1) GO TO 20

      GET1=Z(NZ)

      RETURN

C

20    GET1=Z(NPZR)*DRDT(1)+Z(NPZTH)*DRDT(2)+Z(NPZPH)*DRDT(3)

      RETURN

      END

      FUNCTION GET(Z)

C     'GET' AND ENTRY 'DOT' PROVIDE A CONTROL METHOD FOR AVOIDING

C     REDUNDANT CALLS TO THE TERRAIN AND RECEIVER MODELS. THE VALUES

C     RETURNED ARE THE FUNCTION VALUES FOR 'F' OR 'G' OR THEIR DERIVATI

C     (VIA 'DOT' ENTRY). UNNECESSARY CALLS ARE ELIMINATED THROUGH USE

C     OF 'TIME OF CALL' VARIABLES WHICH ARE COMPARED WITH THE CURRENT

C     LAST CALL TIME MAINTAINED BY THE 'RKAM' PROGRAM. WHEN VALUES ARE

C     NOT CURRENT THEY ARE UPDATED BY CALLS TO THE APPROPRIATE ROUTINES

C

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



      DOUBLE PRECISION Z(*),PF(10),PG(10)

C     COMMON DECK "SS" INSERTED HERE

      DOUBLE PRECISION MODSURF

      COMMON/SS/ MODSURF(4)

      COMMON/SS/U,PUR,PURR,PURTH,PURPH

      COMMON/SS/PUTH,PUPH,PUTHTH,PUPHPH,PUTHPH,USELECT,UTIME

C     COMMON DECK "RKAM" INSERTED HERE

      DOUBLE PRECISION KR,KTH,KPH

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

C     COMMON DECK "FNDER" INSERTED HERE

      COMMON/FNDER/GY(20),NZ,NPZR,NPZRR,NPZRTH,NPZRPH,NPZTH

      COMMON/FNDER/NPZPH,NPZTHTH,NPZPHPH,NPZTHPH,NSELECT,NTIME

C     COMMON DECK "RR" INSERTED HERE

      DOUBLE PRECISION MODREC

      COMMON/RR/ MODREC(4)

      COMMON/RR/F,PFR,PFRR,PFRTH,PFRPH

      COMMON/RR/PFTH,PFPH,PFTHTH,PFPHPH,PFTHPH,FSELECT,FTIME

C     COMMON DECK "GG" INSERTED HERE

      DOUBLE PRECISION MODG

      COMMON/GG/MODG(4)

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

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

C     COMMON DECK "RKTIME" INSERTED HERE

      COMMON/CRKTIME/RKTIME

C

C     COMMON DECK "RMACH" INSERTED HERE

      COMMON/CRMACH/RMACH(5)

      EQUIVALENCE(RKTIME,IRKTIME)

C

C     COMMON DECK "WWR" INSERTED HERE

      PARAMETER (NWARSZ=1000)

      COMMON/WW/dum(10),MAXW,W(NWARSZ)

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

C

      COMMON/CGET/ZERO

C

      DOUBLE PRECISION PI(10),PJ(10)

      COMMON/TMP1/NC,NSZR    ,NSZRR    ,NSZRTH

     1 ,NSZRPH    ,NSZTH    ,NSZPH    ,NSZTHTH    ,NSZPHPH    ,NSZTHPH

     2 ,NVELECT    ,NWIME

C

      DATA PI,PJ/1.0,9*0.0,1.0,9*0.0/

      DATA NC,NSZR,NSZRR,NSZRTH,NSZRPH,NSZTH,NSZPH,NSZTHTH,NSZPHPH

     1   ,NSZTHPH,NVELECT,NWIME/1,2,3,4,5,6,7,8,9,10,11,12/

      DATA FWIME,GWIME,FVELECT,GVELECT/2*-1.0,8HRECEIVER ,7HTERRAIN /

      DATA UWIME,UVELECT/-1.0,7HSURFACE /

      DATA ZHRO/0.0/

C

      IENTRY=1

      GO TO 5

C

      ENTRY SETGET

C

      CALL RMOVE(PF,PI,10)

      CALL RMOVE(PG,PJ,10)

      CALL IMOVE(NZ,NC,12)

      FTIME=FWIME

      GTIME=GWIME

      FSELECT=FVELECT

      GSELECT=GVELECT

      UTIME=UWIME

      USELECT=UVELECT

      ZERO=ZHRO

C

      RETURN

C

      ENTRY DOT(Z)

      IENTRY=2

5     CONTINUE

C

C     THE FOLLOWING THREE ASSIGNMENT STATEMENTS ARE NEEDED TO REPAIR

C     ARRAY ELEMENTS OVERWRITTEN ON THE CRAY VERSION OF THIS PROGRAM.

C     THE ADDITIONAL ARRAY 'GY' HAS ALSO BEEN ADDED TO THE /FNDER/

C     LABELED COMMON TO ACT AS A BUFFER AREA.

C     NO CAUSE FOR THIS PROBLEM HAS YET BEEN FOUND. THESE STATEMENTS

C     ARE NOT NEEDED ON OTHER MACHINES AS FAR AS IS KNOWN.

      NZ=1

      NPZR=2

      NPZRR=3

      IF(ZERO.EQ.0.0) ZERO=EARTHR*RMACH(3)*2.0

C

      IF(ITEST(Z(NTIME)).EQ.IRKTIME) GO TO 10

C

      IF(Z(NSELECT).EQ.USELECT) CALL SURFACE

      IF(Z(NSELECT).EQ.FSELECT) CALL RECEVER

      IF(Z(NSELECT).EQ.GSELECT) CALL TOPOG

      Z(NTIME)=RKTIME

C     REMOVE MACHINE ROUND OFF NOISE FROM EXACT RECEIVER LOCATIONS

      IF(ABS(Z(NZ)).LE.ZERO) Z(NZ)=0.0

C

10    IF(IENTRY.NE.1) GO TO 20

      GET=Z(NZ)

      RETURN

C

20    GET=Z(NPZR)*DRDT(1)+Z(NPZTH)*DRDT(2)+Z(NPZPH)*DRDT(3)

      RETURN

      END

C----------------------------------------------------------------------------

      FUNCTION ITOC(N)

C     RETURN 7 CHARACTER STRING REPRESENTATION OF INTEGER N

C     IF NUMBER IS TOO LARGE OR SMALL USE FLOATING POINT FORMAT

      CHARACTER ITOC*7

      IF(N.LT.-9999.OR.N.GT.99999) GO TO 100

      ITOC=' '

      WRITE(ITOC,'(I7)',ERR=100) N

      RETURN

100   WRITE(ITOC,'(2PG7.0)') FLOAT(N)

      END



      SUBROUTINE LABPLT

C          LABEL THE CURRENT PLOT

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



      CHARACTER*80 LABEL,CHID

      CHARACTER*4 ANGRANG,ANOTE

C     COMMON DECK "LL" INSERTED HERE

      DOUBLE PRECISION MODL

      COMMON/LL/MODL(4),APH,APHPT,APHPR,APHPTH,APHPPH

C     COMMON DECK "SS" INSERTED HERE

      DOUBLE PRECISION MODSURF

      COMMON/SS/ MODSURF(4)

      COMMON/SS/U,PUR,PURR,PURTH,PURPH

      COMMON/SS/PUTH,PUPH,PUTHTH,PUPHPH,PUTHPH,USELECT,UTIME

C     COMMON DECK "GG" INSERTED HERE

      DOUBLE PRECISION MODG

      COMMON/GG/MODG(4)

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

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

C     COMMON DECK "RR" INSERTED HERE

      DOUBLE PRECISION MODREC

      COMMON/RR/ MODREC(4)

      COMMON/RR/F,PFR,PFRR,PFRTH,PFRPH

      COMMON/RR/PFTH,PFPH,PFTHTH,PFPHPH,PFTHPH,FSELECT,FTIME

C     COMMON DECK "UU" INSERTED HERE

      DOUBLE PRECISION MODU

      COMMON/UU/MODU(4)

     1  ,V  ,PVT  ,PVR  ,PVTH  ,PVPH

     2  ,VR ,PVRT ,PVRR ,PVRTH ,PVRPH

     3  ,VTH,PVTHT,PVTHR,PVTHTH,PVTHPH

     4  ,VPH,PVPHT,PVPHR,PVPHTH,PVPHPH

C     COMMON DECK "CC" INSERTED HERE

      DOUBLE PRECISION MODC

      COMMON/CC/MODC(4),CS,PCST,PCSR,PCSTH,PCSPH

C     COMMON DECK "TT" INSERTED HERE

      DOUBLE PRECISION MODT

      COMMON/TT/MODT(4), T,PTT,PTR,PTTH,PTPH

C     COMMON DECK "MM" INSERTED HERE

      DOUBLE PRECISION M,MODM

      COMMON/MM/MODM(4),M,PMT,PMR,PMTH,PMPH

C     COMMON DECK "HDR" INSERTED HERE

      CHARACTER*10 INITID*80,DAT,TOD

      COMMON/HDR/SEC

      COMMON/HDRC/INITID,DAT,TOD

C     COMMON DECK "RINPLEX" INSERTED HERE

      DOUBLE PRECISION KAY2,KAY2I

      COMPLEX PNP,POLAR,LPOLAR

      LOGICAL SPACE

      CHARACTER DISPM*6

      COMMON/RINPL/DISPM

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

      COMMON/RIN/OMEGMIN,OMEGMAX,KAY2,KAY2I

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

      COMMON /DD/ INT,IOR,IT,IS,IC,ICC,IX,IY

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

      character*80 id

      common /ident/ id

C     COMMON DECK "WW" INSERTED HERE

      PARAMETER (NWARSZ=1000)

      COMMON/WW/dum(10),MAXW,W(NWARSZ)

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

      DOUBLE PRECISION MMODEL,MFORM,MID

C

C     WIND        100-124

      EQUIVALENCE (W(100),UMODEL),(W(101),UFORM),(W(102),UID)

C

C     DELTA WIND        125-149

      EQUIVALENCE (W(125),DUMODEL),(W(126),DUFORM),(W(127),DUID)

C

C     SOUND SPEED 150-174

      EQUIVALENCE (W(150),CMODEL),(W(151),CFORM),(W(152),CID)

      EQUIVALENCE (W(153),REFC)

C

C     DELTA SOUND SPEED 175-199

      EQUIVALENCE (W(175),DCMODEL),(W(176),DCFORM),(W(177),DCID)

C

C     TEMPERATURE    200-224

      EQUIVALENCE (W(200),TMODEL),(W(201),TFORM),(W(202),TID)

C

C     DELTA TEMPERATURE    225-249

      EQUIVALENCE (W(225),DTMODEL),(W(226),DTFORM),(W(227),DTID)

C

C     MOLECULAR   250-274

      EQUIVALENCE (W(250),MMODEL),(W(251),MFORM),(W(252),MID)

C

C     RECEIVER HEIGHT   275-299

      EQUIVALENCE (W(275),RMODEL),(W(276),RFORM),(W(277),RID)

C

C     TOPOGRAPHY   300-324

      EQUIVALENCE (W(300),GMODEL),(W(301),GFORM),(W(302),GID)

C

C     DELTA TOPOGRAPHY   325-349

      EQUIVALENCE (W(325),GUMODEL),(W(326),GUFORM),(W(327),GUID)

C

C     ATMOSPHERIC SURFACE TOPOGRAPHY   350-374

      EQUIVALENCE (W(350),SMODEL),(W(351),SFORM),(W(352),SID)

C

C     ATMOSPHERIC SURFACE TOPOGRAPHY   375-399

      EQUIVALENCE (W(375),SUMODEL),(W(376),SUFORM),(W(377),SUID)

C     PLOT ENHANCEMENTS CONTROL PARAMETERS

C

      EQUIVALENCE (W(490),XFQMDL),(W(491),YFQMDL)

C     ABSORPTION        500-524

      EQUIVALENCE (W(500),AMODEL),(W(501),AFORM),(W(502),AID)

C

C     DELTA ABSORPTION        525-549

      EQUIVALENCE (W(525),DAMODEL),(W(526),DAFORM),(W(527),DAID)

C

C     PRESSURE        550-574

      EQUIVALENCE (W(550),PMODEL),(W(551),PFORM),(W(552),PID)

C

C     DELTA PRESSURE        575-599

      EQUIVALENCE (W(575),DPMODEL),(W(576),DPFORM),(W(577),DPID)

C

C

C     COMMON DECK "AA" INSERTED HERE

      DOUBLE PRECISION MODA

      DOUBLE PRECISION MU,MUPT,MUPR,MUPTH,MUPPH

      DOUBLE PRECISION KAP,KAPPT,KAPPR,KAPPTH,KAPPPH

      COMMON/AA/MODA(4),MU,MUPT,MUPR,MUPTH,MUPPH

      COMMON/AA/KAP,KAPPT,KAPPR,KAPPTH,KAPPPH

C     COMMON DECK "PP" INSERTED HERE

      DOUBLE PRECISION MODP

      COMMON/PP/MODP(4),P,PPT,PPR,PPTH,PPPH

C

      WRITE(LABEL,900) ID

900   FORMAT(10A8)

      CHID=LABEL

      LABEL=CHID(5:)

      IOR=0

      IT=0

      IS=2

      IX=0

      IY=1023

      CALL DDTEXT (8,LABEL)

      IX=1090

      IY=0

      CALL DDTEXT (2,DAT)

      IY=1023

      F=OW/PIT2

      NANGLE=10

      ANGRANG='EL ='

      ANOTE='AZ ='

      ANEL=0.

      IF(ELSTEP.NE.0.0) ANEL=(ELEND-ELBEG)/ELSTEP+1.5

      IF(ANEL.GT.1.0.OR.PLT.LT.0.0) GO TO 100

         ANGRANG='AZ ='

         ANOTE='EL ='

         NANGLE=14

100   WRITE(LABEL,1300)  CHID(1:3),F,ANOTE,W(NANGLE)*DEGS

1300  FORMAT ('MODEL = ',A,'  ,FREQ =',F9.3,' HZ, ',A,F7.3,' DEG')

      IX=0

      IY=IY-32

      CALL DDTEXT (7,LABEL)

C

C INDEX OF OPPOSITE ANGLE

      NANGLE=(10+14-NANGLE)+1

      NANGLE2=NANGLE+2

      WRITE(LABEL,1400) ANGRANG,(W(I)*DEGS,I=NANGLE,NANGLE2)

1400  FORMAT(A,F7.2,' DEG TO',F7.2,' DEG, STEP =',F7.2,' DEG')

      IY=IY-32

      CALL DDTEXT (7,LABEL)

C

      WRITE(LABEL,1500) XMTRH,TLAT*DEGS,TLON*DEGS

1500  FORMAT('XMTR HT =',F6.2,' KM ,LAT =',F6.2,' DEG, LONG ='

     1   ,F6.2,' DEG')

      IY=IY-32

      CALL DDTEXT (7,LABEL)

C

      IY=IY-32

      WRITE(LABEL,'(10A8)') MODRIN

      CALL DDTEXT(8,LABEL)

C

      IX=1050

      CALL DDTEXT(1,'MODELS')

      IY=IY-15

C     LOOP FOR 8 MODELS AND PERTURBATIONS

      DO 1700 K=1,12

      I=(K-1)/2+1

C     GENERATE ALTERNATING 1,2;3,4 SERIES FOR MODEL AND PERTURBATION

      J1=2*(K+1-I*2)+1

      J2=J1+1

      IF(I.EQ.1) WRITE(LABEL,1600) (MODU(J),J=J1,J2)

      IF(I.EQ.2) WRITE(LABEL,1600) (MODC(J),J=J1,J2)

C

      IF(I.EQ.3) WRITE(LABEL,1600) (MODG(J),J=J1,J2)

      IF(I.EQ.4) WRITE(LABEL,1600) (MODL(J),J=J1,J2)

      IF(I.EQ.5 .AND. J1.EQ.1)

     1     WRITE(LABEL,1600) (MODREC(J),J=J1,J2)

      IF(I.EQ.6) WRITE(LABEL,1600) (MODSURF(J),J=J1,J2)

1600  FORMAT(2(2(A8,2X,F5.1,1X),2X))

C

1610  IF(LABEL(1:1).EQ.' ') GO TO 1700

      IY=IY-32

      CALL DDTEXT(2,LABEL)

1700  LABEL(1:1)=' '

C

         END

      SUBROUTINE MATRIXM(A,V,R)

C

      DOUBLE PRECISION A(3,3),V(3),R(3)

C

      DO 10 I=1,3

      R(I)=0.0

      DO 10 J=1,3

10    R(I)=R(I)+A(J,I)*V(J)

      END

C-----------------------------------------------------------------------------

      FUNCTION ND2B(INDEC)

C

C     CONVERT A NUMERIC DECIMAL DIGIT STRING TO A BIT STRING

C     WITH EACH BIT SET BY A CORRESPONDING DECIMAL DIGIT.

C

      ND2B=0

      IF(INDEC.LE.0) RETURN

      M=INDEC

      MB=1

C

10    IF(MOD(M,10).NE.0) ND2B=ND2B+MB

      MB=MB*2

      M=M/10

      IF(M.GT.0) GO TO 10

      END

      SUBROUTINE NORMV(V)

C

C     COMPUTE NORMLIZED VECTOR COMPONENTS FOR VECTOR V

C

      DOUBLE PRECISION V(3)

C

      R=V(1)*V(1)+V(2)*V(2)+V(3)*V(3)

      IF(R.LT.1.E-10) RETURN

      R=1.0/SQRT(R)

      DO 10 I=1,3

10    V(I)=R*V(I)

C

      END

C-------------------------------------------------------------------------

      FUNCTION NUMSTG(V,N,FRM)

C     CONVERTS A NUMERIC VALUE TO A STRING

C     USING FORMAT PROVIDED IN CALL

      CHARACTER NUMSTG*80,FRM*(*)

      INTEGER V(N)

C

      NUMSTG=' '

      WRITE(NUMSTG,FRM) V

      RETURN

      END

C----------------------------------------------------------------------------

      SUBROUTINE OVERRD(VAR,TST,DEFT,NFLG,NFVLEQ,NFVLNE)

C

C

C     OVERRIDE SUPPORT ROUTINE

C     TEST 'VAR' AGAINST 'TST', IF EQUAL SET TO DEFAULT 'DEFT'

C     AND ALSO SET INTEGER FLAG 'NFLG' TO VALUE 'NFVLEQ' ELSE 'NFVLNE'



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

C

      IF(VAR.EQ.TST) THEN

         VAR=DEFT

         NFLG=NFVLEQ

      ELSE

         NFLG=NFVLNE

      ENDIF

      RETURN

      END

C-----------------------------------------------------------------------------

      SUBROUTINE PLOTHB(X,Y,NC)

C

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

      EXTERNAL PLOTXY

C

      CALL PLTANH(X,Y,NC,PLOTXY)

      END

C-----------------------------------------------------------------------------

      SUBROUTINE PLOTLB(X,Y,NC)

C

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



      CHARACTER ANNOT*10

      character*80 id

      common /ident/ id

C     COMMON DECK "WWR" INSERTED HERE

      PARAMETER (NWARSZ=1000)

      COMMON/WW/dum(10),MAXW,W(NWARSZ)



      EQUIVALENCE (EARTHR,W(1)),(TICV,W(96))



      COMMON/LABCLT/PROJCT,THMIN,THMAX,RMIN,RMAX

C

      COMMON/DD/IN,IOR,IT,IS,IC,ICC,IX,IY

C

      DATA LIY/-100/

C

      IF(NC.LE.0) GO TO 100

      F=1.

      IF(AMAX1(ABS(RMAX-EARTHR),ABS(RMIN-EARTHR)).LT.1.) F=1000.

      IF(ABS(AINT(TICV*F)-TICV*F).GT..001) GO TO 60

C

      WRITE(ANNOT,50) INT(ABS(Y-EARTHR)*F)

50    FORMAT(I3)

      GO TO 90

C

C.....Added handling of 0.0 value.  Don't want to output "0.0", instead "0"

C     Started at 5 positions farther to the left.



60    if (abs(y-earthr)*f .ne. 0.0) then

           ix = ix - 5

           WRITE(ANNOT,80) ABS(Y-EARTHR)*F

      else

           write(annot,50) 0

      endif



80    FORMAT(F6.2)

C

90    CALL PLOTXY(X-.1*(THMAX-THMIN),Y,1)

C     INSURE THAT OVERLAPS OF ANNOTATIONS DO NOT OCCUR

      IF(IABS(IY-LIY).LT.50) GO TO 100

      LIY=IY



C.....Move the labels to the left.  Aligns better with proj. 3 axis.

      ix = ix-6



      CALL DDTEXT(1,ANNOT)

C

100   CALL PLOTXY(X,Y,IABS(NC))

      END

C----------------------------------------------------------------------------

      SUBROUTINE PLTANOT(ID,FREQ,XMF,YMF,XL,YB,XR,YT

     1    ,DEGS,LLAT,LLON,RLAT,RLON,ALTLOW,ALTHI,PLT,tic,earthr,PLOT)

C     PUTS STANDARD ANNOTATIONS ON PLOTS

C

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



      character*80 id



      external plot



      COMMON/LABCLT/PROJCT,THMIN,THMAX,RMIN,RMAX

      DOUBLE PRECISION LLAT,LLON

      CHARACTER LABEL*80

      COMMON /DD/ INT,IOR,IT,IS,IC,ICC,IX,IY

C     COMMON DECK "ANNOT" INSERTED HERE

      character*13 anotes

      character*32 hnotes

      COMMON/ANNCTL/LENA(4),LENHA(4)

      COMMON/ANNCTC/ANOTES(4),HNOTES(4)

ccccc CHARACTER*10 ANOTES,HNOTES

ccccc COMMON/ANNCTL/LENA(4),LENHA(3)

ccccc COMMON/ANNCTC/ANOTES(2,4),HNOTES(4,3)

C

      IF(PLT.EQ.2. .OR. XMF.LE.0.0 .OR. YMF.LE.0.0) GO TO 45

      CALL PLOT(XL+XMF*(XR-XL),YB+YMF*(YT-YB),1)

      IF(FREQ.LE.0.0) GO TO 30

C

      WRITE(LABEL,25) FREQ

25    FORMAT('FREQ =',F9.3)

      CALL DDTEXT(2,LABEL)

      IY=IY+40

C

30    WRITE(LABEL,35) ID

35    FORMAT('MODEL = ',A3)

      CALL DDTEXT(2,LABEL)

C

45    NSPY=97

      IOR=0

C

C FIRST THE left lat/lon ANNOTATION

50    CALL PLOT(XL,YB,1)

      IX=IX-155

      IYsave=IY-NSPY  !!!!!

      iy = 80    !!!!! for VT display, converts to last row - 2

      WRITE(LABEL,1840) DEGS*LLON

      CALL DDTEXT(2,LABEL)

      WRITE(LABEL,1850) DEGS*LLAT

ccccc      IY = iy - 32

      iy = 45    !!!!! for VT display, converts to last row - 1

      CALL DDTEXT(2,LABEL)

C

C NEXT THE RIGHT LATITUDE/Longitude ANNOTATION

      CALL PLOT(XR,YB,1)

      IX=IX-130

      iy = 80    !!!!! for VT display, converts to last row - 2

ccccc      IY=IY-NSPY

      WRITE(LABEL,1840) DEGS*RLON

      CALL DDTEXT(2,LABEL)

      WRITE(LABEL,1850) DEGS*RLAT

      iy = 45    !!!!! for VT display, converts to last row - 1

ccccc      IY = iy - 32

      CALL DDTEXT(2,LABEL)

C

      iy = iysave

      XMID=(XL+XR)/2.0

      IF(PLT.EQ.2.0) GO TO 55

C

C PUT THE HORIZONTAL TICK LABEL

      Y=YT

      IF(ALTLOW.GE.0.0) Y=YB

      CALL PLOT(XMID,Y,1)

C

      IF(ALTLOW.LT.0.0) IY=IY+80

      IF(ALTLOW.GE.0.0) IY=IY-95

      GO TO 60

C

55    CALL PLOT(XMID,YB,1)

      if (plt .eq. 2.) iy = 80

ccccc      IY=IY-95

C

60    IX=IX-180   !!!!! changed 235 to 180

ccccc      iy = iy-7    !!!!! put the e.g. RANGE AT SEA LEVEL lower

      NOTEA=1

ccccc      TMP=(THMAX-THMIN)*DEGS

ccccc      IF(TMP.GT.10.) NOTEA=2

C

      degtick = tic*degs*1.e+5

      rkmtick = tic*earthr*1000

      degtmp = MOD(degtick,1)

      rkmtmp = MOD(rkmtick,1)



      if (ABS(degtmp-NINT(degtmp)) .lt. ABS(rkmtmp-NINT(rkmtmp))) then

           notea = 2

      endif



      CALL DDTEXT (LENHA(NOTEA),HNOTES(NOTEA))

C

C PUT THE VERTICAL TICK LABEL

      CALL PLOT(XL,(YB+YT)/2.0,1)

      IOR=1

C

C   HORIZONTAL PLOT:  PUT Y-AXIS ANNOTATION

C

      IF(PLT.eq.2.) then

           IX=IX-125

           IY=200

C..........Use the 4th element of HNOTES if we used the 2nd element for

C          the horizontal axis i.e. using DEGrees instead of KM for both

C          labels.

C

           if (notea .eq. 2) then

                notea = 4

           else

                notea = 3

           endif

           CALL DDTEXT (LENHA(notea),HNOTES(notea))

      else

C

C..........Vertical projection Y-axis labelling.

C

           IX=IX-125

           IY=IY-75

           NOTEA=1

           IF(ALTLOW.GE.0.0) NOTEA=3

           IF(ABS(ALTHI-ALTLOW).GE.1.) NOTEA=NOTEA+1

           CALL DDTEXT (LENA(NOTEA),ANOTES(NOTEA))

      endif

C

      IOR=0

1840  FORMAT(F7.2,' Deg E.')

1850  FORMAT(F7.2,' Deg N.')

C

      END

C------------------------------------------------------------------------------

      Subroutine Flush (numpts,ixtmp,iytmp,bpxtmp,bpytmp)



C.....This routine is called by DDCON after every contour level is processed

C     or

C     in the event that NUMPTS reaches 2000.



C     Flush the buffers and reset NUMPTS to zero.



C.....NOTE:  Multiply NUMPTS by 2 before writing to TXTOUT.  This

C            makes PSGRAPH more straightforward.  In the CONTOUR

C            program, number of points indicates a "moveto,lineto"

C            pair.  In PSGRAPH, the each moveto is a point and each

C            lineto is a point.



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



      integer*2 numpts, ixtmp(2000), iytmp(2000), bpxtmp(2000),

     &          bpytmp(2000)





      if (numpts.eq. 0) RETURN



      write(ndevtxt,33) 4

      write(ndevtxt,44) numpts*2

33    format(I2)

44    format(I4)

      do ii = 1,numpts-3,4

           write(ndevtxt,55) bpxtmp(ii),bpytmp(ii),ixtmp(ii),iytmp(ii),

     &             bpxtmp(ii+1),bpytmp(ii+1),ixtmp(ii+1),iytmp(ii+1),

     &             bpxtmp(ii+2),bpytmp(ii+2),ixtmp(ii+2),iytmp(ii+2),

     &             bpxtmp(ii+3),bpytmp(ii+3),ixtmp(ii+3),iytmp(ii+3)

55         format(16(I5))

      enddo



      if (MOD(numpts,4) .ne. 0) then

           numremain = MOD(numpts,4)

           do jj = numpts-numremain+1, numpts

                write(ndevtxt,66) bpxtmp(jj),bpytmp(jj),

     &                            ixtmp(jj),iytmp(jj)

           enddo

      endif

66    format (4(I5))



      numpts = 0



      return

      end

