C file S8.FOR
C ----------------------------------------------------------------------!-------
C 23 Jun 2009 - Changed the default character size for DDTEXT calls to 1.
C 30 Apr 2009 - replaced REAL*8 with DOUBLE PRECISION
C 10 Sep 2008 - Made data types consistent.
C 13 Mar 2008 - Corrected length in common block /PCONST/
C 13 Mar 2008 - Brought common block /HDR/ up to date
C  7 Mar 2008 - changed ABS to CABS, SQRT to DSQRT
C  7 Mar 2008 - changed AMAX1 to DMAX1, AMOD to DMOD, SIN to DSIN, COS to DCOS
C  2 Nov 2006 - moved decimal plot labels to the left in PLOTLB
C  2 Nov 2006 - lengthened format for decimal plot labels in PLOTLB
C 30 Mar 2006 - Moved text location for plot header label
C 27 Mar 2006 - Fixed type for common block ddscale
C 24 Mar 2006 - Moved DDCON, PLOTXY, MX1ALF, MX2ALF, DDSC, DDST, FULANN, PLOTHB,
C               PLOTLB from contour1
C 24 Mar 2006 - Moved KURVE, PLOTINI, PLTTLT,CPLOT,YTRIM,F,CLIP,KODE from eigen2
C 24 Mar 2006 - Moved WRITW, MARKER from eigen
C 24 Mar 2006 - Moved GRACE, GRAFB, FULANN from eigen4
C 24 Mar 2006 - Moved COMPRS, CURVE, DONEPL, ENDPL, MX1ALF, MX2ALF, NCLOSE, 
C               NOBRDR, PAGE, PHYSOR from eigen3 and pcprof1
C ----------------------------------------------------------------------!-------
CONTAINS SUBROUTINE COMPRS, SUBROUTINE CURVE, SUBROUTINE DONEPL, 
CONTAINS SUBROUTINE ENDPL, SUBROUTINE MX1ALF, SUBROUTINE MX2ALF, 
CONTAINS FUNCTION NCLOSE, SUBROUTINE NOBRDR, SUBROUTINE PAGE, SUBROUTINE PHYSOR, 
CONTAINS SUBROUTINE GRACE, SUBROUTINE GRAFB, SUBROUTINE FULANN, 
C SUBROUTINE SCLPIC, SUBROUTINE SUPNEG, SUBROUTINE TITLE, SUBROUTINE TITLXY, 
C SUBROUTINE VMNMX, SUBROUTINE XTICKS, SUBROUTINE YTICKS, SUBROUTINE TITLTX, 
C SUBROUTINE WRITW, SUBROUTINE MARKER, SUBROUTINE KURVE, SUBROUTINE PLOTINI, 
C SUBROUTINE PLTTLT, FUNCTION F, SUBROUTINE CPLOT, SUBROUTINE YTRIM, 
C FUNCTION CLIP, FUNCTION KODE, SUBROUTINE DDCON, SUBROUTINE PLOTXY, 
C SUBROUTINE DDSC, SUBROUTINE DDST, SUBROUTINE PLOTHB, SUBROUTINE PLOTLB, 
C Subroutine Flush 
C ----------------------------------------------------------------------!-------
      SUBROUTINE COMPRS
C ----------------------------------------------------------------------!-------
C     INITIALIZE THE METACODE FILE(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
C     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT         ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
C     NDEVTXT is rewound here because this is the first graphics code
C     routine.
C
c      print *, ' nryind  = ', nryind
c      print *, ' NDEVtmp = ', NDEVtmp
c      print *, ' Nfrmat  = ', Nfrmat
c      print *, ' NDEVtxt = ', NDEVtxt
c      print *, ' NDEVbin = ', NDEVbin
      REWIND NDEVTXT
      WRITE(NDEVTXT,100) 20
100   FORMAT(I2)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE COMPRS
C----------------------------------------------------------------------------
C ----------------------------------------------------------------------!-------
      SUBROUTINE CURVE(X,Y,N,IC)
C ----------------------------------------------------------------------!-------
C     DRAWS A CURVE THROUGH A SET OF DATA POINTS WITH OR WITHOUT
C     MARKER SYMBOLS(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
C     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT         ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
      DOUBLE PRECISION X(N),Y(N)

      if (n.eq.0) return

      WRITE(NDEVTXT,100) 29,N,IC
      if (n.ge.2) then
           WRITE(NDEVTXT,200) (X(I),Y(I),X(I+1),Y(I+1),I=1,N-1,2)
C..........Handle the case of an odd number of data points.
           if (MOD(n,2).ne.0) WRITE(NDEVTXT,300) X(N),Y(N)
      else if (n.eq.1) then
           WRITE(NDEVTXT,300) X(1),Y(1)
      endif
100   FORMAT(I2,2I5)
200   FORMAT(4E14.6)
300   FORMAT(2E14.6)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE CURVE
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE DONEPL
C ----------------------------------------------------------------------!-------
C     TERMINATES PLOTTING SESSION(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
C     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT         ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
C
      WRITE(NDEVTXT,100) 31
100   FORMAT(I2)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE DONEPL
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE ENDPL(N)
C ----------------------------------------------------------------------!-------
C     TERMINATES PLOTTING ON CURRENT PAGE(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
C     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT         ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
C
      WRITE(NDEVTXT,100) 30,N
100   FORMAT(I2,I3)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE ENDPL
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE MX1ALF(T1,T2)
C ----------------------------------------------------------------------!-------
C     SPECIFY USE OF ALTERNATE CHARACTER SET NUMBER 1(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
C     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT         ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
      CHARACTER*10 T1,T2

      WRITE(NDEVTXT,100) 11, T1, T2
100   FORMAT(I2,TR1,A,1X,A1)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE MX1ALF
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE MX2ALF(T1,T2)
C ----------------------------------------------------------------------!-------
C     SPECIFY USE OF ALTERNATE CHARACTER SET NUMBER 2(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
C     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT         ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
      CHARACTER*10 T1,T2

      WRITE(NDEVTXT,200) 12,T1,T2
200   FORMAT(I2,TR1,A,1X,A1)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE MX2ALF
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      FUNCTION NCLOSE(AM,AX,AI,NTICKS,AXIS)
C ----------------------------------------------------------------------!-------
C
C     ADDITIONAL LOGIC TO INSURE THAT DIGITS ARE NOT PACKED
C     TOO CLOSE TOGETHER. RULE OF THUMB: KEEP NUMBER OF DIGITS PER
C     INCH LESS THAN 4.
C
      implicit double precision (a-h,o-z)

      NCLOSE=0
      IF(AI.LE.0.0) RETURN
C
      IF(AI .GT. 0.1*(AX-AM)) THEN
         NCLOSE=1
         RETURN
      ENDIF
C
      NDS=dlog10(dmax1(dabs(AM),dabs(AX)))
C
      NDS=MOD(NDS,5)+3
      IF(AM.LT.0.0) NDS=NDS+1
      NCLOSE=(NDS*(AX-AM)/AI/AXIS-1.0)/4.0+1
      IF(NCLOSE.LE.1) RETURN
      AI=NCLOSE*AI
      NTICKS=4
      RETURN
C ----------------------------------------------------------------------!-------
      END FUNCTION NCLOSE
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE NOBRDR
C ----------------------------------------------------------------------!-------
C     SUPPRESSES DRAWING THE PAGE BORDER(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
C     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT         ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006

      WRITE(NDEVTXT,100) 36
100   FORMAT(I2)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE NOBRDR
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE PAGE(X,Y)
C ----------------------------------------------------------------------!-------
C     DEFINE LIMITS OF PHYSICAL PAGE(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
C     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT         ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
      WRITE(NDEVTXT,100) 23,X,Y
100   FORMAT(I2,2F6.2)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE PAGE
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE PHYSOR(X,Y)
C     DEFINE LOCATION OF PHYSICAL ORIGIN OF PLOTTING AREA(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
C     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT         ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
      WRITE(NDEVTXT,100) 22,X,Y
100   FORMAT(I2,2F5.2)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE PHYSOR
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE GRACE(X)
C ----------------------------------------------------------------------!-------
      DOUBLE PRECISION x                                   ! 23 August 2005
C     DEFINE MARGIN AROUND PLOTTING AREA(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
C     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT         ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
C
      WRITE(NDEVTXT,100) 21,X
100   FORMAT(I2,F4.2)
      RETURN                                     ! 23 August 2005
C ----------------------------------------------------------------------!-------
      END SUBROUTINE GRACE
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE GRAFB(A,B,C,D,E,F)
C ----------------------------------------------------------------------!-------
C     ROUTINE WHICH ESTABLISHES BOXED X AND Y AXES(DISSPLA)
C     THIS IS A PSEUDO CALL WHICH USES SEVERAL OTHER DISSPLA CALLS
C     IT IS ASSUMED THAT A PREVIOUS CALL TO TITLE HAS BEEN MADE
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
      DOUBLE PRECISION a,b,c,d,e,f
C     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT         ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
C
      WRITE(NDEVTXT,100) 39,A,B,C,D,E,F
100   FORMAT(I2,6E13.7)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE GRAFB
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE FULANN
C ----------------------------------------------------------------------!-------
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
C ----------------------------------------------------------------------!-------
      ENTRY SETANN
C ----------------------------------------------------------------------!-------
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

C ----------------------------------------------------------------------!-------
      ENTRY ANNFIL(S,C)
C ----------------------------------------------------------------------!-------
      C=S
C ----------------------------------------------------------------------!-------
      END SUBROUTINE FULANN
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE SCLPIC(S)
C ----------------------------------------------------------------------!-------
C     SCALES CURVE MARKER SYMBOLS(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
C     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT         ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
C
      WRITE(NDEVTXT,100) 24,S
100   FORMAT(I2,F5.2)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE SCLPIC
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE SUPNEG(IX,IY)
C ----------------------------------------------------------------------!-------
C     SET COUNTERS FOR NUMBER OF NEGATIVE AXIS VALUES TO MAKE POSITIVE
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE

c     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVGRP,NDEVBIN,NDEVTXT ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
      WRITE(NDEVTXT,100) 35,IX,IY
100   FORMAT(I2,2I3)
C ----------------------------------------------------------------------!-------
      END  SUBROUTINE SUPNEG
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE TITLE(A,B,C,D,E,F,G,H) 
C ----------------------------------------------------------------------!-------
C     DEFINES TITLE, AXES LABELS AND LENGTHS OF PLOTS(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C

C     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT         ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
C*****THESE 3 ARRAYS NEED TO BE REAL.  IF DECLARED CHAR*40, THE DIRECT
C*****ACCESS READ IN TAPRD WILL CAUSE EXTRA CHARACTERS TO BE PRESENT?!

ccccc      REAL A(4),C(4),E(4)
      character*(*) a,c,e
      integer b,d,f
      DOUBLE PRECISION G,H
C
C*****CHANGED ENTRY TO INCLUDE THE PARAMETER LIST.  MICROSOFT WOULD NOT
C*****COMPILE WITHOUT IT.  JAH.
C
C ----------------------------------------------------------------------!-------
      ENTRY TITLEW (A,B,C,D,E,F,G,H)
C ----------------------------------------------------------------------!-------
      WRITE(NDEVTXT,100) 25,G,H
100   FORMAT(I2,2F6.2)
      CALL TITLTX (A,B)
      CALL TITLTX (C,D)
      CALL TITLTX (E,F)
C
C ----------------------------------------------------------------------!-------
      END SUBROUTINE TITLE
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE TITLXY(I,J,K,L,M,N,XAXIS,YAXIS)
C ----------------------------------------------------------------------!-------
C TITLES FOR X AND Y AXES 
C ----------------------------------------------------------------------!-------
      implicit double precision (a-h,o-z)

      CHARACTER I*(*),K*(*),M*(*),CD*80 
      PARAMETER (NWARSZ=1000) 
      COMMON/WW/ dum(10), MAXW,W(NWARSZ)
      EQUIVALENCE (HITLET,W(75))
C
      F=1.0
C.....Changed from 6.65 inches. Too big for 8.5x11
      XAXIS=5.0

C.....Changed from 6.8.
      YAXIS=5.11

      IF(J.GT.0) GO TO 50
      F=.8
      XAXIS=9.0
      YAXIS=6.0
C
C.....Sets origin of plot to 2.3", 1.66" from LLeft
50    CALL PHYSOR(2.3,1.66)
      CALL NOBRDR
      IF(HITLET.LE.0.0) CALL HEIGHT(F*0.4)
      CALL SCLPIC(1.4)
      CALL TITLEW(CD(I),J,CD(K),L,CD(M),N,XAXIS,YAXIS)
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE TITLXY
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE VMNMX(V,N,  VMN,VMX)
C ----------------------------------------------------------------------!-------
C FINDS MINIMUM AND MAXIMUM OF AN ARRAY 
C ----------------------------------------------------------------------!-------
      implicit double precision (a-h,o-z)

      DOUBLE PRECISION V(N)
      IF(N.LT.1) RETURN
      DO 100 I=1,N
      VMN=dMIN1(V(I),VMN)
100   VMX=dMAX1(V(I),VMX)
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE VMNMX
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE XTICKS(IX)
C ----------------------------------------------------------------------!-------
C     SPECIFIES NUMBER OF SUBTICKS PER MAJOR TICK MARK(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
c     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVGRP,NDEVBIN,NDEVTXT ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
      WRITE(NDEVTXT,100) 32,IX
100   FORMAT(I2,I3) 
C ----------------------------------------------------------------------!-------
      END SUBROUTINE XTICKS
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE YTICKS(IX)
C ----------------------------------------------------------------------!-------
C     SPECIFIES NUMBER OF SUBTICKS PER MAJOR TICK MARK(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
c     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVGRP,NDEVBIN,NDEVTXT ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
      WRITE(NDEVTXT,100) 33,IX
100   FORMAT(I2,I3)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE YTICKS
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE TITLTX (ACHAR,B)
C ----------------------------------------------------------------------!-------
C     THIS IS CALLED BY TITLE IN ORDER TO WRITE TO NDEVTXT.
C
C     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT         ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
      CHARACTER*(*) ACHAR
c      LOGICAL FOUND
      INTEGER B
C
c      NMBLK=0
c      FOUND = .FALSE.
      lena = len(achar)
c      DO 10 I=lena,1,-1
c         IF (.NOT.FOUND) THEN
c            IF(ACHAR(I:I).EQ.' ') NMBLK=NMBLK+1
c            IF(ACHAR(I:I).EQ.'$') FOUND = .TRUE.
c         ENDIF
c10    CONTINUE
c      LENA = 40-NMBLK
C
      WRITE(NDEVTXT,200) B,ACHAR(1:lena) !!!!!,I=1,LENA)
200   FORMAT(I5,1X,A)
C
C ----------------------------------------------------------------------!-------
      END SUBROUTINE TITLTX
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
C			Subroutine WRITW
C ----------------------------------------------------------------------!-------
C
c.....NOTE:  This routine is no longer used. Its function has been replaced
c            by Ron Brent's modifications. jah.
C ----------------------------------------------------------------------
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE WRITW(NDEV,NDEVTMP,ID,W,WS,MAXW)
C ----------------------------------------------------------------------!-------
C
      implicit double precision (a-h,o-z)
 
      character id*80,line*80
      DOUBLE PRECISION W(MAXW),WS(MAXW)
C*      REAL LINE(8)
C
C
c      print *, ' NDEVtmp = ', NDEVtmp
      REWIND NDEVTMP
      DO 200 I=1,1000
           READ(NDEVTMP,190,END=10000) LINE
C*190   FORMAT(BZ,8A10)
C     SKIP THE ID
           WRITE(NDEV,190) LINE
200   CONTINUE
190   FORMAT(BZ,A80)
C
c      print *, ' NDEVtmp = ', NDEVtmp
10000 REWIND NDEVTMP
c      print *, ' NDEVtmp = ', NDEVtmp
C
      DO 100 I=1,MAXW
           IF(W(I).NE.WS(I)) WRITE(NDEV,20) I,sngl(W(I))
100        WS(I)=W(I)
20    FORMAT(I3,G14.7)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE WRITW
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE MARKER(MK)
C ----------------------------------------------------------------------!-------
C     SPECIFIES MARKER SYMBOL TO BE USED BY CURVE ROUTINE(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C     COMMON DECK "RAYDEV" INSERTED HERE
C     DEVICE ASSIGNED TO RAYTRC INPUT FILE
C     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT         ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
 
      WRITE(NDEVTXT,100) 28,MK,0
100   FORMAT (3I2)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE MARKER
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
C			Subroutine KURVE
C ----------------------------------------------------------------------!-------
C
C ----------------------------------------------------------------------
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE KURVE(X,Y,NPTS,MARK,XMIN,XMAX)
C ----------------------------------------------------------------------!-------
 
	implicit double precision (a-h,o-z)
 
C+P   (KURVE) MODIFIED VRS OF DISSPLA 'CURVE'. DOES CLIPPING AND
C+    MARKER SELECTION.
      DOUBLE PRECISION XW(333),YW(333),X(1),Y(1)
C     COMMON DECK "WWE" INSERTED HERE
      PARAMETER (NWARSZ=1000)
      COMMON/WW/dum(10),MAXW,W(NWARSZ)
      equivalence (w(436),rmax)
c
      COMMON/CALTC/INI,FX,FY,XL,XR,YB,YT
C
      XL=XMIN
      XR=XMAX
      YB=0.0
      YT=RMAX
      NW=333
      CALL CPLOT(X,Y,NPTS,XW,YW,NW)
      IF(NW.LE.0) RETURN
C
      IF(MARK.NE.0) CALL MARKER(MARK)
      N=0
      IF(MARK.NE.0) N=-1
      CALL CURVE (XW,YW,NW,N)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE KURVE
C ----------------------------------------------------------------------!-------
 
C ----------------------------------------------------------------------!-------
C			Subroutine PLOTINI
C ----------------------------------------------------------------------!-------
C
C ----------------------------------------------------------------------!-------
       SUBROUTINE PLOTINI(DHITE,SYMRED)
C ----------------------------------------------------------------------!-------
	
	implicit double precision (a-h,o-z)
 
C+P   (PLOTTING) INITIALIZE PLOTTING FOR EIGEN PROGRAM
C
C     COMMON DECK "HDR" INSERTED HERE
      CHARACTER*10 INITID*80,DAT,TOD
      COMMON/HDR/SEC,stsec                            !Rev 7/13/90 ! 13 Mar 2008
      COMMON/HDRC/INITID,DAT,TOD
C     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT         ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
C
C.....Write a code to TXTOUT indicating that eigen plots are being made.
C
      write (ndevtxt,111) 40
111   format (i2)
 
      CALL COMPRS
C     CALL GRACE(0.,0.)                               ! 23 August 2005
      CALL GRACE(0.D0)                                ! 23 August 2005
      CALL PHYSOR(2.,1.3)
      CALL PAGE(12.,10.)
      CALL HEIGHT(DHITE)
      CALL SETANN
      S=SYMRED
C  DEFAULT TO 15% REDUCTION ON POINT SIZE OVER DISSPLAY'S DEFAULT
      IF(S.EQ.0.0) S=.15
      CALL SCLPIC(S)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE PLOTINI
C ----------------------------------------------------------------------!-------
 
C ----------------------------------------------------------------------!-------
C			Subroutine PLTTLT
C
C ----------------------------------------------------------------------!-------
 
      SUBROUTINE PLTTLT(NCODE)
C ----------------------------------------------------------------------!-------
	
	implicit double precision (a-h,o-z)
 
C+P   (PLTTLT) PUT TITLES ON PLOTS FOR EIGEN PROGRAM
C
C     COMMON DECK "WWE" INSERTED HERE
      PARAMETER (NWARSZ=1000)
	character id*80
      common /ident/id
      COMMON/WW/dum(10),MAXW,W(NWARSZ)
      DOUBLE PRECISION NRSLMI
      EQUIVALENCE (W(80),APOG) ,(W(401),PTYP) ,(W(403),ELMIN)
     1   ,(W(404),ELMAX) ,(W(405),ELVINC) ,(W(406),RGINC)
     2   ,(W(407),HOPLME) ,(W(408),NRSLME)
     3   ,(W(409),GAPLO), (W(410),GAPHI), (W(411),SYMRED)
     4   ,(W(428),GRPMIN) ,(W(429),GRPMAX) ,(W(430),GRPVINC)
     5   ,(W(432),HOPLMG) ,(W(433),NRSLMG)
     6   ,(W(434),RFSPEED) ,(W(435),TMFT) ,(W(436),RMAX)
     &   ,(W(437),rmin)
     7   ,(W(453),RCRNG) ,(W(454),SLPCUT)
     8   ,(W(456),HOPLMI) ,(W(457),NRSLMI) ,(W(458),ESET)
      COMMON/GRPC/GRPN,GRPX
C
C     COMMON DECK "CONST" INSERTED HERE
      COMMON/PCONST/CREF,RGAS,GAMMA
     1 ,MO2,MN2,MO,MN,MH2O,MO3                        ! 9 Mar 2007 ! 13 Mar 2008
      DOUBLE PRECISION MO2,MN2,MO,MN,MH2O,MO3                      ! 13 Mar 2008
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10
C     COMMON DECK "CODES" INSERTED HERE
      COMMON/CODES/ICODRE,ICODRT,ICODEI
C
C     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT         ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
 
      CHARACTER TDID*30
      COMMON/CDATIM/TDID
	character label*80	
	equivalence (label,ilabel)
C
      DATA NTCHRS/21/
C
C.....Changed 7.8 by 6.9 to 5.65 by 5.0 for plot dimensions.
 
C.....Write a code to TXTOUT indicating that eigen plots are being made.
C
      write (ndevtxt,111) 40
111   format (i2)
 
      IF(NCODE.EQ.ICODRE) CALL TITLE (TDID,NTCHRS
     1     ,'ELEV. ANGLE OF TRANSMISSION (#DEG!)$',100
     2     ,'RANGE (#KM!)$',100
     3     ,5.65d0,5.0d0)
 
      IF(NCODE.EQ.ICODRT.AND.TMFT.EQ.1.) CALL TITLE (TDID,NTCHRS
     1     ,'RELATIVE TRAVEL TIME (#SEC!)$',100
     2     ,'RANGE (#KM!)$',100
     3     ,5.65d0,5.0d0)
 
      IF(NCODE.EQ.ICODRT.AND.TMFT.NE.1.) CALL TITLE (TDID,NTCHRS
     1     ,'RELATIVE TRAVEL TIME (#MS!)$',100
     2     ,'RANGE (#KM!)$',100
     3     ,5.65d0,5.0d0)
 
 
CCCCC      CALL FRAME  !!!!!
      ELINC=DEGS*ELVINC
      IF(ELVINC.LE.0.0) ELINC=DEGS*(ELMAX-ELMIN)/10.0
      GRPINC=GRPVINC
      IF(GRPVINC.LE.0.0) GRPINC=(GRPMAX-GRPMIN)/10.0
 
C.....Add the following text to TXTOUT.  This is a subheading for the
C     plot.  Eigen plots are essentially identical to profile plots,
C     therefore there must be a subheading for PSGRAPH.
 
c      write (ndevtxt,100) 3, 0, 0
      write (ndevtxt,100) 3, 335, 1150                     ! 30 Mar 2006
100   format (i2,2i4)
      write (ndevtxt,200) 0,80
200   format (2i3)
      write (ndevtxt,300) id
300   format (a)
 
      IF(NCODE.EQ.ICODRE)
     1   CALL GRAFB(DEGS*ELMIN,ELINC,DEGS*ELMAX,rmin,RGINC,RMAX)
      IF(NCODE.EQ.ICODRT)
     1  CALL GRAFB(GRPMIN*TMFT,GRPINC*TMFT,GRPMAX*TMFT,rmin,RGINC, RMAX)
C
 
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE PLTTLT
C ----------------------------------------------------------------------!-------
 
C ----------------------------------------------------------------------!-------
C			Function F
C
C ----------------------------------------------------------------------!-------
 
 
      character FUNCTION F(C)
C ----------------------------------------------------------------------!-------
 
	implicit double precision (a-h,o-z)
 
C+P   MODIFY ANNOTATIONS FOR EIGEN PROGRAM DEPENDING ON ANNOTATION
C+    MODEL.
      CHARACTER C*(*)
C
      CALL ANNFIL(C,F)
      RETURN
C ----------------------------------------------------------------------!-------
      END FUNCTION F
C ----------------------------------------------------------------------!-------
 
C ----------------------------------------------------------------------!-------
C			Subroutine CPLOT
C
C ----------------------------------------------------------------------!-------
 
 
      SUBROUTINE CPLOT(X,Y,NP,XW,YW,NW)
C ----------------------------------------------------------------------!-------
 
	implicit double precision (a-h,o-z)
 
C+P   (CPLOT) CLIPPED PLOT LINE PROGRAM. USES STANDARD CLIPPING
C+    ALGORITHM(CLIP) TO DETERMINE BREAKPOINTS AND INSERTS ANY
C+    LINE SEGMENTS TO PRODUCE CONNECTED LINE.
C
      INTEGER CLIP
      LOGICAL OUTSIDE
      DOUBLE PRECISION X(*),Y(*),XW(*),YW(*)
C
      COMMON/CALTC/INI,FX,FY,XL,XR,YB,YT
      DATA XS,YS/2*0.0/
C
C     HANDLE SINGLES SPECIALLY
      IF(NP.EQ.1) GO TO 180
C
C     FIRST TRIM THE POINTS BASED ON Y VALUES
      NPP=NP+1
      NM=NP
      CALL YTRIM(Y,YT,NM)
      NM=NM-1
      N=1
      IF(NM.LT.1) GO TO 140
C
      DO 100 I=1,NM
      XI=X(I)
      YI=Y(I)
      XX=X(I+1)
      YY=Y(I+1)
      OUTSIDE=CLIP(XI,YI,XX,YY).NE.0
      IF(OUTSIDE) GO TO 100
      IF((XS.EQ.XI.AND.YS.EQ.YI).OR.N.EQ.1) GO TO 50
      IF(N.EQ.NW) then
ccccc      CALL STOPIT('LOGIC ER 1')
           stop 'n = nw' !!??????
       endif
C      WRITE(6,111) I,sngl(XI),sngl(YI),sngl(XS),sngl(YS),
C     1             sngl(XX),sngl(YY)
111   FORMAT(I5,10F8.2)
      XW(N)=XS
      YW(N)=YS
      N=N+1
50    XS=XX
      YS=YY
      IF(N.EQ.NW) then
ccccc      CALL STOPIT('LOGIC ER 1')
           stop 'n = nw' !!??????
      endif      !!!CALL STOPIT('LOGIC ER 2')
      XW(N)=XI
      YW(N)=YI
      N=N+1
100   CONTINUE
C
      IF(.NOT.OUTSIDE) GO TO 150
140   NW=N-1
      CALL YTRIM(YW,YT,NW)
      RETURN
150   XW(N)=XX
      YW(N)=YY
      NW=N
      CALL YTRIM(YW,YT,NW)
C
      RETURN
C
180   NW=0
      IF(X(1).GT.XR) RETURN
      IF(X(1).LT.XL) RETURN
      IF(Y(1).GT.YT) RETURN
      IF(Y(1).LT.YB) RETURN
      XW(1)=X(1)
      YW(1)=Y(1)
      NW=1
C ----------------------------------------------------------------------!-------
      END SUBROUTINE CPLOT
C ----------------------------------------------------------------------!-------
 
C ----------------------------------------------------------------------!-------
C			Subroutine YTRIM
C
C ----------------------------------------------------------------------!-------
 
      SUBROUTINE YTRIM(Y,YT,NP)
C ----------------------------------------------------------------------!-------
 
	implicit double precision (a-h,o-z)
 
C+P   (YTRIM) ELIMINATES SEGMENTS OF LINE FROM HIGH END WHOSE
C+    Y VALUES EXCEED GIVEN VALUE.
      DOUBLE PRECISION Y(*)
      IF(NP.LT.2) RETURN
      NPP=NP+1
      K=NPP
      DO 10 I=1,NP
         IF(Y(NPP-I).LT.YT) GO TO 15
10    K=K-1
15    NP=MIN0(K,NP)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE YTRIM
C ----------------------------------------------------------------------!-------
 
C ----------------------------------------------------------------------!-------
C			Function CLIP
C
C ----------------------------------------------------------------------!-------
 
      INTEGER FUNCTION CLIP(X1,Y1,X2,Y2)
C ----------------------------------------------------------------------!-------
 
	implicit double precision (a-h,o-z)
 
C+P   (CLIP) PERFORM STANDARD CLIPPING ON SINGLE LINE SEGMENTS
C+    (X1,Y1;X2,Y2) AGAINST WINDOW (XL,YB;XR,YT). 'CLIP' TAKES
C+    ON NONZERO VALUE IF LINE SEGMENT IS COMPLETELY OUTSIDE WINDOW.
C
      LOGICAL SWAPPED
      INTEGER C1,C2
      COMMON/CALTC/INI,FX,FY,XL,XR,YB,YT
C
      SWAPPED=.FALSE.
      C2=KODE(X2,Y2)
C
C     BEGIN "COMPUTE" LOOP
C
C     WE TREAT PT ONE AS SPECIAL
10    C1=KODE(X1,Y1)
      CLIP=IAND(C1,C2)
ccccc111   FORMAT(4F8.2,3(1X,O4),L5)  !!!!!not used!?
C     IF BOTH END PTS ARE OFF-SCREEN, SO IS LINE
      IF(CLIP.NE.0) GO TO 100
C     IF BOTH END PTS ARE ON-SCREEN SO IS THE LINE
      IF(C1.EQ.0 .AND. C2.EQ.0) GO TO 100
C
      IF(C1.NE.0) GO TO 15
C     AT LEAST ONE END POINT IS OFF-SCREEN, MAKE SURE PT ONE IS ONE
      SWAPPED=.NOT.SWAPPED
      I=C1
      C1=C2
      C2=I
      CALL SWAP2(X1,X2,Y1,Y2)
C
C     SUCCESSIVELY REPLACE LINE BY ITS PROJECTIONS ALONG PRIMARY DIVIDER
15    IF(IAND(C1,1).EQ.0) GO TO 20
        Y1=Y1+(Y2-Y1)*(XL-X1)/(X2-X1)
        X1=XL
        GO TO 10
20    IF(IAND(C1,8).EQ.0) GO TO 30
        Y1=Y1+(Y2-Y1)*(XR-X1)/(X2-X1)
        X1=XR
        GO TO 10
30    IF(IAND(C1,64).EQ.0) GO TO 40
        X1=X1+(X2-X1)*(YB-Y1)/(Y2-Y1)
        Y1=YB
        GO TO 10
40    IF(IAND(C1,512).EQ.0) GO TO 10
        X1=X1+(X2-X1)*(YT-Y1)/(Y2-Y1)
        Y1=YT
C
      GO TO 10
100   IF(SWAPPED) CALL SWAP2(X1,X2,Y1,Y2)
C
C ----------------------------------------------------------------------!-------
      END FUNCTION CLIP
C ----------------------------------------------------------------------!-------
 
C ----------------------------------------------------------------------!-------
C			Function KODE
C
C ----------------------------------------------------------------------!-------
 
      FUNCTION KODE(X,Y)
C ----------------------------------------------------------------------!-------
C VALUE IS ZERO IF (X,Y) WITHIN PLOTTING AREA
C IF OUTSIDE OF PLOTTING AREA, HAS OTHER VALUES THAT DEPEND 
C ON WHICH OCTANT OUTSIDE OF PLOTTING AREA
C THE VALUES ARE AS FOLLOWS:
C    513 | 512 | 520
C    ---------------
C     1  |  0  |  8 
C    ---------------
C     63 |  64 |  72
C ----------------------------------------------------------------------!-------
 
	implicit double precision (a-h,o-z)
 
C
      COMMON/CALTC/INI,FX,FY,XL,XR,YB,YT
C
      IF(X.GE.XL) GO TO 10
        K=1
        GO TO 100
10    IF(X.LE.XR) GO TO 50
        K=8
        GO TO 100
C
50    K=0
100   IF(Y.GE.YB) GO TO 110
        K=IOR(64,K)
      GO TO 200
110   IF(Y.LE.YT) GO TO 200
        K=IOR(512,K)
C
200   KODE=K
      RETURN
C ----------------------------------------------------------------------!-------
      END FUNCTION KODE
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE DDCON (X, Y, Z, CC, MX, MY, MZ, MC)
C ----------------------------------------------------------------------!-------

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

C      real qqq,qq,scx,scy
      DOUBLE PRECISION qqq,scx,scy                                 ! 27 Mar 2006

C      COMMON /DDSCALE /QQQ (8), SCX, SCY, QQ (4), ISCX, ISCY
      COMMON /DDSCALE /QQQ(4), MQQQ(4), SCX, SCY, MQQ(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      PRINT*,'ddcon 1: ix= ',ix,'iy = ',iy

      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
c      PRINT*,'ddcon 2: ix= ',ix,'iy = ',iy

C      CALL DDVC_ACUM (numpts,ixtmp,iytmp, outflag)
      CALL DDVC_ACUM (numpts,ixtmp,iytmp)                         ! 3 Mar 2006

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
C ----------------------------------------------------------------------!-------
      END SUBROUTINE DDCON 
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE PLOTXY(X,Y,NC)
C ----------------------------------------------------------------------!-------
C PLOT ONE POINT 
C ----------------------------------------------------------------------!-------
      implicit double precision (a-h,o-z)

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

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

C      real xmin,xmax,ymin,ymax,scx,scy
      REAL(KIND=8) :: xmin,xmax,ymin,ymax,scx,scy          ! 27 Mar 2006

      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

c      PRINT*,'plotxy: x = ',x,'y = ',y,'ddhix= ',ddhix,'ddhiy = ',ddhiy
      SCX=X
      SCY=Y
      CALL DDSC
      IX=ISCX
      IY=ISCY
c      PRINT*,'plotxy: ix= ',ix,'iy = ',iy

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)                  ! 3 Mar 2006
C           CALL DDVC_ACUM (numpts, ixtmp, iytmp, outflag)
      endif

      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE PLOTXY
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE PLTANH_test(X,Y,NC,PLOT)
C ----------------------------------------------------------------------!-------
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
     1 ,MO2,MN2,MO,MN,MH2O,MO3                        ! 9 Mar 2007 ! 13 Mar 2008
      DOUBLE PRECISION MO2,MN2,MO,MN,MH2O,MO3                      ! 13 Mar 2008
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10
      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.0d0)                                  ! 10 Sep 2008
      rkmtmp = MOD(rkmtick,1.0d0)                                  ! 10 Sep 2008

      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(DMOD(TIC*FT+.0001,1.d0)).GT..001) GO TO 60 ! 7 Mar 2008 10 Sep 2008
C
C.....Shift label to right if it has only 1 or 2 places.

      if (INT(v*f+SIGN(.5d0,v)) .ge. 0 .and. 
     &    INT(v*f+SIGN(.5d0,v)) .lt. 100)
     &     ix = ix + 5                                             ! 10 Sep 2008
      WRITE(ANNOT,50) INT(V*F+SIGN(.5d0,V))                        ! 10 Sep 2008
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.d0,H)                     ! 10 Sep 2008
      IOR=0
      IS=1                                                        ! 23 Jun 2009
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))
C ----------------------------------------------------------------------!-------
      END SUBROUTINE PLTANH_test
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE PLTHLB_test(X,Y,NC)
C ----------------------------------------------------------------------!-------
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)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE PLTHLB_test
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE PLTLB_test(X,PY,NCP)
C ----------------------------------------------------------------------!-------
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_test(X,PY,NCP)
      NC=0
C
10    Y=PY
C
      IF(NCP.LT.0) GO TO 200
      F=1.
      IF(DMAX1(DABS(RMAX-EARTHR),DABS(RMIN-EARTHR)).LT.1.) F=1000. ! 7 Mar 2008
      Y=INT((Y-EARTHR)*F)/F+EARTHR
      IF(NC.EQ.0) GO TO 100
C
      V=Y
C      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(I4)  ! 23 jun 05
      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
      IS=1                                                        ! 23 Jun 2009
      CALL DDTEXT(1,ANNOT)
C
100   CALL PLOT(X,Y,NC)
      RETURN
C
200   CALL PLOT(X,Y,1)
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE PLTLB_test
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------
      SUBROUTINE TIKLINE_test(XL1,YB,XL,YT1,TICV,TIKSZ,PLOT)
C ----------------------------------------------------------------------!-------

      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)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE TIKLINE_test
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------
       SUBROUTINE DDSC
C ----------------------------------------------------------------------!-------
C***E     DDSC
C***.
C PLOT ONE VECTOR FOR CONTOUR PLOTTING 
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)
C      real xm,scx, xmin,xmax
      REAL(KIND=8) :: xm,scx, xmin,xmax                    ! 27 Mar 2006
      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/

c      PRINT*,'ddsc enter: iscx = ',iscx(1),'iscy = ',iscx(2)
c      PRINT*,'ddsc enter: mscx = ',mscx(1),'mscy = ',mscx(2)
      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 (.5d0, Q)                                     ! 10 Sep 2008
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
c      PRINT*,'ddsc exit: iscx = ',iscx(1),'iscy = ',iscx(2)
c      PRINT*,'ddsc exit: mscx = ',mscx(1),'mscy = ',mscx(2)
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE DDSC
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE DDST
C ----------------------------------------------------------------------!-------
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
     1 ,MO2,MN2,MO,MN,MH2O,MO3                        ! 9 Mar 2007 ! 13 Mar 2008
      DOUBLE PRECISION MO2,MN2,MO,MN,MH2O,MO3                      ! 13 Mar 2008
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10
C     COMMON DECK "WWR" INSERTED HERE
      PARAMETER (NWARSZ=1000)
      COMMON/WW/dum(10),MAXW,W(NWARSZ)
      EQUIVALENCE (EARTHR,W(1)),(PFACTR,W(82)),(PLT,W(81))

C      real xmin,xmax,ymin,ymax,scx,scy
      REAL(KIND=8) :: xmin,xmax,ymin,ymax,scx,scy          ! 27 Mar 2006

      COMMON /DDSCALE/ XMIN,XMAX,YMIN,YMAX,MINX,MAXX,MINY,MAXY,SCX,SCY,
     1 NSCX,NXCY,MSCX,MSCY,ISCX,ISCY
C
C
      COSC(X,Y)=DSIN(Y)*STH0*DCOS(X-PH0)+DCOS(Y)*CTH0 ! 7 Mar 2008
      SCFLIP(X)=DSQRT(DMAX1(0.d0,1.-X*X)) ! 7 Mar 2008             ! 10 Sep 2008
      SINB(X,Y)=DSIN(X-PH0)*DSIN(Y)/SCFLIP(COSCC) ! 7 Mar 2008
      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*DCOS(ALPHA) ! 7 Mar 2008
      XMAX=DMAX1(YMAX*DSIN(ALPHA),(YMAX-YMIN)/2.) ! 7 Mar 2008
      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)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE DDST
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE PLOTHB(X,Y,NC)
C ----------------------------------------------------------------------!-------
C TICK LABELING 
C ----------------------------------------------------------------------!-------
      implicit double precision (a-h,o-z)
      EXTERNAL PLOTXY
C
c      PRINT*,'plothb: x = ',x,'y = ',y,'ddhix= ',ddhix,'ddhiy = ',ddhiy

      CALL PLTANH(X,Y,NC,PLOTXY)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE PLOTHB
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE PLOTLB(X,Y,NC)
C ----------------------------------------------------------------------!-------
C TICK LABELING 
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
c      PRINT*,'plotlb: x = ',x,'y = ',y,'ddhix= ',ddhix,'ddhiy = ',ddhiy

      iroom = 0            ! usually do not need more room ! 2 Nov 2006

      IF(NC.LE.0) GO TO 100
      F=1.
      IF(DMAX1(DABS(RMAX-EARTHR),DABS(RMIN-EARTHR)).LT.1.) F=1000. ! 7 Mar 2008
      IF(ABS(AINT(TICV*F)-TICV*F).GT..001) GO TO 60
C
      WRITE(ANNOT,50) INT(ABS(Y-EARTHR)*F)
50    FORMAT(I4)                     ! 23 June 05
      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

      iroom = 15   ! give more room for labels with fractions - 2 Nov 2006

C80    FORMAT(F6.2)
80    FORMAT(F7.2)                                            ! 2 Nov 2006
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.
C      ix = ix-6
      ix = ix-6 - iroom ! more room for labels with fractions ! 2 Nov 2006
c      PRINT*,'plotlb: ix= ',ix,'iy = ',iy

      IS=1                                                        ! 23 Jun 2009
      CALL DDTEXT(1,ANNOT)
C
100   CALL PLOTXY(X,Y,IABS(NC))
C ----------------------------------------------------------------------!-------
      END SUBROUTINE PLOTLB
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE PLTANOT_test(ID,FREQ,XMF,YMF,XL,YB,XR,YT
     1    ,DEGS,LLAT,LLON,RLAT,RLON,ALTLOW,ALTHI,PLT,tic,earthr,PLOT)
C ----------------------------------------------------------------------!-------
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)
      IS=1                                                        ! 23 Jun 2009
      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  !!!!!
C     iy = 80    !!!!! for VT display, converts to last row - 2
C     Above line commented June 23 2005
      IY=40
      WRITE(LABEL,1840) DEGS*LLON
      CALL DDTEXT(2,LABEL)
      WRITE(LABEL,1850) DEGS*LLAT
ccccc      IY = iy - 32
C     iy = 45    !!!!! for VT display, converts to last row - 1
C     Above line commented June 23 2005
      IY=5
      CALL DDTEXT(2,LABEL)
C
C NEXT THE RIGHT LATITUDE/Longitude ANNOTATION
      CALL PLOT(XR,YB,1)
      IX=IX-130
C     iy = 80    !!!!! for VT display, converts to last row - 2
C     Above line commented June 23 2005
      IY=40
ccccc      IY=IY-NSPY
      WRITE(LABEL,1840) DEGS*RLON
      CALL DDTEXT(2,LABEL)
      WRITE(LABEL,1850) DEGS*RLAT
C     iy = 45    !!!!! for VT display, converts to last row - 1
C     Above line commented June 23 2005
      IY=5
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.0d0)                                  ! 10 Sep 2008
      rkmtmp = MOD(rkmtick,1.0d0)                                  ! 10 Sep 2008

      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 labeling.
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
C ----------------------------------------------------------------------!-------
      END SUBROUTINE PLTANOT_test
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      Subroutine Flush (numpts,ixtmp,iytmp,bpxtmp,bpytmp)
C ----------------------------------------------------------------------!-------

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
C ----------------------------------------------------------------------!-------
      end Subroutine Flush
C ----------------------------------------------------------------------!-------
