C*********************************************************************
C***  HEADER  ********************************************************
C
C  GBOOK - A SMALL HISTOGRAM PACKAGE
C  WRITTEN DECEMBER 1979, LAST CHANGED MARCH 1983
C  AUTHOR: TORBJORN SJOSTRAND, DEPT. OF THEORETICAL PHYSICS,
C  UNIVERSITY OF LUND, SOLVEGATAN 14A, S-223 62 LUND, SWEDEN
C  PRESENT ADDRESS: FNAL, THEORY, TEL. 312-840-3753
C  PLEASE REPORT ANY ERRORS TO THE AUTHOR
C
C***  INTRODUCTION  **************************************************
C
C  GBOOK IS A SMALL SUBROUTINE PACKAGE FOR GETTING ONE- OR TWO-
C  DIMENSIONAL HISTOGRAMS ON AN ORDINARY LINE PRINTER OR TERMINAL.
C  IT CAN BE USED IN A MANNER VERY SIMILAR TO HBOOK, BUT HAS BEEN
C  WRITTEN COMPLETELY INDEPENDENTLY. THE CODE IS FULLY IN FORTRAN77,
C  WITH A MINIMUM OF MACHINE DEPENDENCE. THE USAGE CAN BE DIVIDED INTO
C  FOUR STEPS: BOOKING, FILLING, EDITING AND PRINTING. ALL SUBROUTINES
C  HAVE NAMES SIX CHARACTERS LONG AND BEGINNING WITH G.
C
C***  BOOKING  *******************************************************
C
C  THERE ARE NMAX HISTOGRAMS AT THE DISPOSAL OF THE USER, EACH GIVEN BY
C  A NUMBER BETWEEN 1 AND NMAX. BEFORE A HISTOGRAM CAN BE USED, SPACE
C  MUST BE RESERVED FOR IT. NMAX IS A COMPILE-TIME PARAMETER ADDED BY
C  MIKE SEYMOUR.
C
C  CALL GBOOK1(ID,TITLE,NX,XL,XU)
C  PURPOSE: BOOK A ONE-DIMENSIONAL HISTOGRAM.
C  ID: HISTOGRAM NUMBER, INTEGER BETWEEN 1 AND NMAX.
C  TITLE: HISTOGRAM TITLE, CAN BE GIVEN EITHER AS A CHARACTER STRING
C     OF AT MOST 60 CHARACTERS OR AS A CHARACTER*60 VARIABLE.
C  NX: NUMBER OF BINS (IN X DIRECTION) IN HISTOGRAM, INTEGER BETWEEN
C     1 AND 100.
C  XL, XU: LOWER AND UPPER BOUND, RESPECTIVELY, ON THE (X) RANGE
C     COVERED BY THE HISTOGRAM.
C
C  CALL GBOOK2(ID,TITLE,NX,XL,XU,NY,YL,YU)
C  PURPOSE: BOOK A TWO-DIMENSIONAL HISTOGRAM.
C  ID: HISTOGRAM NUMBER, INTEGER BETWEEN 1 AND NMAX.
C  TITLE: HISTOGRAM TITLE, SEE GBOOK1.
C  NX: NUMBER OF BINS IN X DIRECTION, INTEGER BETWEEN 1 AND 50.
C  XL, XU: LOWER AND UPPER BOUND ON X RANGE OF HISTOGRAM.
C  NY: NUMBER OF BINS IN Y DIRECTION, ARBITRARY POSITIVE INTEGER.
C  YL, YU: LOWER AND UPPER BOUND OF Y RANGE OF HISTOGRAM.
C
C***  FILLING  *******************************************************
C
C  FOR BOOKED HISTOGRAMS WEIGHTS CAN BE FILLED AT GIVEN COORDINATES.
C
C  CALL GFILL1(ID,X,W)
C  PURPOSE: FILL IN A ONE-DIMENSIONAL HISTOGRAM.
C  ID: HISTOGRAM NUMBER.
C  X: X COORDINATE OF POINT.
C  W: WEIGHT TO BE ADDED IN THIS POINT.
C
C  CALL GFILL2(ID,X,Y,W)
C  PURPOSE: FILL IN A TWO-DIMENSIONAL HISTOGRAM.
C  ID: HISTOGRAM NUMBER.
C  X: X COORDINATE OF POINT.
C  Y: Y COORDINATE OF POINT.
C  W: WEIGHT TO BE ADDED IN THIS POINT.
C
C***  EDITING  *******************************************************
C
C  FOR EDITING OF HISTOGRAMS BEFORE PRINTOUT TWO ROUTINES ARE AVAILABLE.
C
C  CALL GSCALE(ID,F)
C  PURPOSE: RESCALE THE CONTENTS OF A HISTOGRAM.
C  ID: HISTOGRAM NUMBER.
C  F: RESCALING FACTOR, I.E. FACTOR THAT ALL BIN CONTENTS (INCLUDING
C     OVERFLOW ETC.) ARE MULTIPLIED BY.
C  REMARK: A TYPICAL RESCALING FACTOR FOR A ONE-DIMENSIONAL HISTOGRAM
C     COULD BE F = 1/(BIN SIZE * NUMBER OF EVENTS) =
C     = NX/(XU-XL) * 1/(NUMBER OF EVENTS).
C
C  CALL GOPERA(ID1,OPER,ID2,ID3,F1,F2)
C  PURPOSE: THIS IS A GENERAL PURPOSE ROUTINE FOR EDITING ONE OR SEVERAL
C     HISTOGRAMS, WHICH ALL ARE ASSUMED TO HAVE THE SAME NUMBER OF
C     BINS. OPERATIONS ARE CARRIED OUT BIN BY BIN, INCLUDING OVERFLOW
C     BINS ETC.
C  OPER: GIVES THE TYPE OF OPERATION TO BE CARRIED OUT, A ONE-CHARACTER
C     STRING OR A CHARACTER*1 VARIABLE.
C  OPER= '+', '-', '*', '/': ADD, SUBTRACT, MULTIPLY OR DIVIDE THE
C     CONTENTS IN ID1 AND ID2 AND PUT THE RESULT IN ID3. F1 AND F2, IF
C     NOT 1., GIVE FACTORS BY WHICH THE ID1 AND ID2 BIN CONTENTS ARE
C     MULTIPLIED BEFORE THE INDICATED OPERATION. (DIVISION WITH A
C     VANISHING BIN CONTENT WILL GIVE 0.)
C  OPER= 'A', 'S', 'L': FOR 'S' THE SQUARE ROOT OF THE CONTENT IN ID1
C     IS TAKEN (RESULT 0 FOR NEGATIVE BIN CONTENTS) AND FOR 'L' THE
C     10-LOGARITHM IS TAKEN (A NONPOSITIVE BIN CONTENT IS BEFORE THAT
C     REPLACED BY 0.8 TIMES THE SMALLEST POSITIVE BIN CONTENT).
C     THEREAFTER, IN ALL THREE CASES, THE CONTENT IS MULTIPLIED BY F1
C     AND ADDED WITH F2, AND THE RESULT IS PLACED IN ID3. THUS ID2
C     IS DUMMY IN THESE CASES.
C  OPER= 'M': INTENDED FOR STATISTICAL ANALYSIS, BIN-BY-BIN MEAN AND
C     STANDARD DEVIATION OF A VARIABLE, ASSUMING THAT ID1 CONTAINS
C     ACCUMULATED WEIGHTS, ID2 ACCUMULATED WEIGHT*VARIABLE AND
C     ID3 ACCUMULATED WEIGHT*VARIABLE-SQUARED. AFTERWARDS ID2 WILL
C     CONTAIN THE MEAN VALUES (=ID2/ID1) AND ID3 THE STANDARD
C     DEVIATIONS (=SQRT(ID3/ID1-(ID2/ID1)**2)). IN THE END, F1
C     MULTIPLIES ID1 (FOR NORMALIZATION PURPOSES), WHILE F2 IS DUMMY.
C
C***  PRINTING  ******************************************************
C
C  AT PRINTING AXES ARE CHOSEN SUCH THAT HISTOGRAMS ARE SUPPOSED
C  TO FIT INTO ONE PAGE. FOR ONE-DIMENSIONAL HISTOGRAMS NUMBERS
C  SMALLER IN MAGNITUDE THAN 10**(-10) ARE CONSIDERED TO BE 0,
C  OTHERWISE SCALES ARE CHOSEN AUTOMATICALLY FOR MAXIMUM RESOLUTION.
C  TWO-DIMENSIONAL HISTOGRAMS ARE STRONGLY ORIENTED TOWARDS HAVING
C  WEIGHTS IN THE ORDER OF OR BIGGER THAN UNITY : SIGNS 1 - 9 ARE
C  CHOSEN IN A LINEAR SCALE UP TO WEIGHTS 9.5, WHEREAS A - Z WILL
C  EITHER BE CHOSEN IN A LINEAR SCALE WITH STEP 1 OR, IF MAXIMUM
C  WEIGHT IS LARGER THAN 36.5, IN A LOGARITHMICALLY EVEN SCALE.
C  NEGATIVE NUMBERS ARE ALLOWED BOTH IN ONE- AND TWO-DIMENSIONAL
C  HISTOGRAMS, AND ARE PRECEDED BY A - SIGN.
C
C  CALL GCLEAR
C  PURPOSE: PRINT OUT ALL HISTOGRAMS THAT HAVE BEEN FILLED, AND
C     RESET THEM THEREAFTER TO 0.
C
C  CALL GPRINT(ID)
C  PURPOSE: PRINT OUT A SINGLE HISTOGRAM.
C  ID: HISTOGRAM TO BE PRINTED.
C
C  CALL GRESET(ID)
C  PURPOSE: RESET ALL BIN CONTENTS, INCLUDING OVERFLOW ETC., TO 0.
C  ID: HISTOGRAM TO BE RESET.
C
C***  COMMON BLOCK AND SPACE REQUIREMENTS  ***************************
C
C  A COMMONBLOCK
C  PARAMETER (NSIZE=200000,NMAX=2000)
C  COMMON /GBOOK/ A(NSIZE)
C  IS USED TO STORE HISTOGRAM INFORMATION. THE HISTOGRAM INDEX TAKES
C  NMAX+2 POSITIONS. EACH BOOKED ONE- (TWO-) DIMENSIONAL HISTOGRAM TAKES
C  AN ADDITIONAL 38+NX (38+NX*NY)  POSITIONS. THE PROGRAM HAS TO BE
C  RECOMPILED WITH CHANGED COMMONBLOCK IF MORE IS REQUIRED.
C
C***  END  DESCRIPTION  **********************************************
      SUBROUTINE IDATE_MOD(iyr, imo, iday)
C     Drop-in replacement for old 'idate' subroutine
C     Returns year, month, day as integers

      INTEGER iyr, imo, iday
      INTEGER values(8)

      CALL DATE_AND_TIME(values=values)

      iyr = values(1)
      imo = values(2)
      iday = values(3)

      RETURN
      END
C*********************************************************************
      SUBROUTINE GBOOK1(ID,TITLE,NX,XL,XU)
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NSIZE=200000,NMAX=2000)
      COMMON /GBOOK/ A(NSIZE)
      CHARACTER TITLE*(*),TITFX*60
      EQUIVALENCE (REQ,IEQ)
      IF (ID.GT.NMAX .OR. A(1)+A(2)+38+NX.GT.NSIZE) THEN
         WRITE (6,200) ID,NMAX,INT(A(1)+A(2)+38+NX+0.5),NSIZE
         RETURN
      ENDIF
      A(ID+2)=A(1)+A(2)
      A(2)=A(2)+38+NX
      IS=A(ID+2)+0.5
      A(IS+1)=NX
      A(IS+2)=XL
      A(IS+3)=XU
      A(IS+4)=(XU-XL)/NX
      A(IS+5)=1
      CALL GRESET(ID)
      TITFX=TITLE//' '
      DO 100 IT=1,20
      IEQ=256**2*ICHAR(TITFX(3*IT-2:3*IT-2))+256*ICHAR(TITFX(3*IT-1:
     &3*IT-1))+ICHAR(TITFX(3*IT:3*IT))
  100 A(IS+18+NX+IT)=REQ
      RETURN
  200 FORMAT (' ERROR: Too much space requested in GBOOK1!'/
     &  ' Requested   ID=',I4,',    Maximum   ID=',I4/
     &  ' Requested space=',I6,', Maximum space=',I6/
     &  ' Recompile with larger NMAX and/or NSIZE')
      END
C*********************************************************************
      SUBROUTINE GBOOK2(ID,TITLE,NX,XL,XU,NY,YL,YU)
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NSIZE=200000,NMAX=2000)
      COMMON /GBOOK/ A(NSIZE)
      CHARACTER TITLE*(*),TITFX*60
      EQUIVALENCE (REQ,IEQ)
      IF (ID.GT.NMAX .OR. A(1)+A(2)+38+NX*NY.GT.NSIZE) THEN
         WRITE (6,200) ID,NMAX,INT(A(1)+A(2)+38+NX*NY+0.5),NSIZE
         RETURN
      ENDIF
      A(ID+2)=A(1)+A(2)
      A(2)=A(2)+38+NX*NY
      IS=A(ID+2)+0.5
      A(IS+1)=NX
      A(IS+2)=XL
      A(IS+3)=XU
      A(IS+4)=(XU-XL)/NX
      A(IS+5)=NY
      A(IS+6)=YL
      A(IS+7)=YU
      A(IS+8)=(YU-YL)/NY
      CALL GRESET(ID)
      TITFX=TITLE//' '
      DO 100 IT=1,20
      IEQ=256**2*ICHAR(TITFX(3*IT-2:3*IT-2))+256*ICHAR(TITFX(3*IT-1:
     &3*IT-1))+ICHAR(TITFX(3*IT:3*IT))
  100 A(IS+18+NX*NY+IT)=REQ
      RETURN
  200 FORMAT (' ERROR: Too much space requested in GBOOK2!'/
     &  ' Requested   ID=',I4,',    Maximum   ID=',I4/
     &  ' Requested space=',I6,', Maximum space=',I6/
     &  ' Recompile with larger NMAX and/or NSIZE')
      END
C*********************************************************************
      SUBROUTINE GFILL1(ID,X,W)
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NSIZE=200000,NMAX=2000)
      COMMON /GBOOK/ A(NSIZE)
      IF (ID.GT.NMAX) RETURN
      IS=A(ID+2)+0.5
      A(IS+9)=A(IS+9)+1.
      IOX=2
      IF(X.LT.A(IS+2)) IOX=1
      IF(X.GE.A(IS+3)) IOX=3
      A(IS+12+IOX)=A(IS+12+IOX)+W
      IF(IOX.NE.2) RETURN
      IX=(X-A(IS+2))/A(IS+4)
      A(IS+19+IX)=A(IS+19+IX)+W
      RETURN
      END
C*********************************************************************
      SUBROUTINE GFILL2(ID,X,Y,W)
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NSIZE=200000,NMAX=2000)
      COMMON /GBOOK/ A(NSIZE)
      IF (ID.GT.NMAX) RETURN
      IS=A(ID+2)+0.5
      A(IS+9)=A(IS+9)+1.
      IOX=2
      IF(X.LT.A(IS+2)) IOX=1
      IF(X.GE.A(IS+3)) IOX=3
      IOY=2
      IF(Y.LT.A(IS+6)) IOY=1
      IF(Y.GE.A(IS+7)) IOY=3
      A(IS+6+3*IOY+IOX)=A(IS+6+3*IOY+IOX)+W
      IF(IOX.NE.2.OR.IOY.NE.2) RETURN
      IX=(X-A(IS+2))/A(IS+4)
      IY=(Y-A(IS+6))/A(IS+8)
      IC=INT(A(IS+1)+0.5)*IY+IX
      A(IS+19+IC)=A(IS+19+IC)+W
      RETURN
      END
C*********************************************************************
      SUBROUTINE GSCALE(ID,F)
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NSIZE=200000,NMAX=2000)
      COMMON /GBOOK/ A(NSIZE)
      IF (ID.GT.NMAX) RETURN
      IS=A(ID+2)+0.5
      DO 100 IC=IS+10,IS+18+INT(A(IS+1)+0.5)*INT(A(IS+5)+0.5)
  100 A(IC)=F*A(IC)
      RETURN
      END
C*********************************************************************
      SUBROUTINE GOPERA(ID1,OPER,ID2,ID3,F1,F2)
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NSIZE=200000,NMAX=2000)
      COMMON /GBOOK/ A(NSIZE)
      CHARACTER OPER*1
      IF (ID1.GT.NMAX .OR. ID2.GT.NMAX .OR. ID3.GT.NMAX) RETURN
      IS1=A(ID1+2)+0.5
      IS2=A(ID2+2)+0.5
      IS3=A(ID3+2)+0.5
      NC=INT(A(IS3+1)+0.5)*INT(A(IS3+5)+0.5)
      IF(OPER.EQ.'+'.OR.OPER.EQ.'-'.OR.OPER.EQ.'*'.OR.OPER.EQ.'/')
     &A(IS3+9)=A(IS1+9)+A(IS2+9)
      IF(OPER.EQ.'A'.OR.OPER.EQ.'S'.OR.OPER.EQ.'L') A(IS3+9)=A(IS1+9)
      IF(OPER.EQ.'+') THEN
      DO 100 IC=10,18+NC
  100 A(IS3+IC)=F1*A(IS1+IC)+F2*A(IS2+IC)
      ELSEIF(OPER.EQ.'-') THEN
      DO 110 IC=10,18+NC
  110 A(IS3+IC)=F1*A(IS1+IC)-F2*A(IS2+IC)
      ELSEIF(OPER.EQ.'*') THEN
      DO 120 IC=10,18+NC
  120 A(IS3+IC)=F1*A(IS1+IC)*F2*A(IS2+IC)
      ELSEIF(OPER.EQ.'/') THEN
      DO 130 IC=10,18+NC
      FA2=F2*A(IS2+IC)
      IF(ABS(FA2).LE.1E-10) A(IS3+IC)=0.
  130 IF(ABS(FA2).GT.1E-10) A(IS3+IC)=F1*A(IS1+IC)/FA2
      ELSEIF(OPER.EQ.'A') THEN
      DO 140 IC=10,18+NC
  140 A(IS3+IC)=F1*A(IS1+IC)+F2
      ELSEIF(OPER.EQ.'S') THEN
      ZERO=0
      DO 150 IC=10,18+NC
  150 A(IS3+IC)=F1*SQRT(MAX(ZERO,A(IS1+IC)))+F2
      ELSEIF(OPER.EQ.'L') THEN
      ZMIN=1E30
      DO 160 IC=19,18+NC
  160 IF(A(IS1+IC).LT.ZMIN.AND.A(IS1+IC).GT.1E-20) ZMIN=0.8*A(IS1+IC)
      DO 170 IC=10,18+NC
  170 A(IS3+IC)=F1*LOG10(MAX(A(IS1+IC),ZMIN))+F2
      ELSEIF(OPER.EQ.'M') THEN
      DO 180 IC=10,18+NC
      IF(ABS(A(IS1+IC)).LE.1E-10) A(IS2+IC)=0.
      IF(ABS(A(IS1+IC)).GT.1E-10) A(IS2+IC)=A(IS2+IC)/A(IS1+IC)
      IF(ID3.NE.0.AND.ABS(A(IS1+IC)).LE.1E-10) A(IS3+IC)=0.
      ZERO=0
      IF(ID3.NE.0.AND.ABS(A(IS1+IC)).GT.1E-10) A(IS3+IC)=
     &SQRT(MAX(A(IS3+IC)/A(IS1+IC)-A(IS2+IC)**2,ZERO))
  180 A(IS1+IC)=F1*A(IS1+IC)
      ENDIF
      RETURN
      END
C*********************************************************************
      SUBROUTINE GCLEAR
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NSIZE=200000,NMAX=2000)
      COMMON /GBOOK/ A(NSIZE)
      IF (ID.GT.NMAX) RETURN
      DO 100 ID=1,INT(A(1)+0.5)
      IS=A(ID+2)+0.5
      IF(IS.EQ.0.OR.A(IS+9).LT.0.5) GOTO 100
      CALL GPRINT(ID)
      CALL GRESET(ID)
  100 CONTINUE
      RETURN
      END
C*********************************************************************
      SUBROUTINE GPRINT(ID)
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER GLAST
      PARAMETER (NSIZE=200000,NMAX=2000)
      COMMON /GBOOK/ A(NSIZE)
      CHARACTER TITLE*60, CTIME*8, OUT*108, CHA(40)*1
      DIMENSION IROW(100), IFRA(100), DYAC(10), EV(20)
      EQUIVALENCE (REQ,IEQ)
      DATA DYAC/.04,.05,.06,.08,.10,.12,.15,.20,.25,.30/, LIN/41/
      DATA CHA/' ','0','1','2','3','4','5','6','7','8','9','A','B',
     &'C','D','E','F','G','H','I','J','K','L','M','N','O','P','Q',
     &'R','S','T','U','V','W','X','Y','Z','*','-','!'/
      IF (ID.GT.NMAX) RETURN
      IS=A(ID+2)+0.5
      IF(A(IS+9).LT.0.5) WRITE(6,400) ID
      IF(A(IS+9).LT.0.5) RETURN
      NX=INT(A(IS+1)+0.5)
      NY=INT(A(IS+5)+0.5)
      DO 100 IT=1,20
      EV(IT)=0.
      REQ=A(IS+18+NX*NY+IT)
  100 TITLE(3*IT-2:3*IT)=CHAR(IEQ/256**2)//CHAR(MOD(IEQ,256**2)/256)
     &//CHAR(MOD(IEQ,256))
      CALL IDATE_MOD(IMON,IDAY,IYEAR)
      CALL DATE_AND_TIME(CTIME)
      WRITE(6,410) ID, TITLE, IYEAR, IMON, IDAY, CTIME(1:5)
      IF(NY.EQ.1) THEN
      YMIN=A(IS+19)
      YMAX=A(IS+19)
      DO 120 IX=IS+20,IS+18+NX
      IF(A(IX).LT.YMIN) YMIN=A(IX)
  120 IF(A(IX).GT.YMAX) YMAX=A(IX)
      IF(YMAX-YMIN.GT.LIN*DYAC(1)*1E-9) THEN
      IF(YMIN.GT.0..AND.YMIN.LT.YMAX/10.) YMIN=0.
      IF(YMAX.LT.0..AND.YMAX.GT.YMIN/10.) YMAX=0.
      IPOT=INT(LOG10(YMAX-YMIN)+10.)-10
      IF(YMAX-YMIN.LT.LIN*DYAC(1)*10.**IPOT) IPOT=IPOT-1
      IF(YMAX-YMIN.GT.LIN*DYAC(10)*10.**IPOT) IPOT=IPOT+1
      DELY=DYAC(1)
      DO 130 IDEL=1,9
  130 IF(YMAX-YMIN.GE.LIN*DYAC(IDEL)*10.**IPOT) DELY=DYAC(IDEL+1)
      DY=DELY*10.**IPOT
      DO 140 IX=1,NX
      CTA=ABS(A(IS+18+IX))/DY
      IROW(IX)=SIGN(CTA+0.95,A(IS+18+IX))
  140 IFRA(IX)=10.*(CTA+1.25-AINT(CTA+0.95))
      IRMI=SIGN(ABS(YMIN)/DY+0.95,YMIN)
      IRMA=SIGN(ABS(YMAX)/DY+0.95,YMAX)
      DO 160 IR=IRMA,IRMI,-1
      IF(IR.EQ.0) GOTO 160
      OUT=' '
      DO 150 IX=1,NX
      IF(IR.EQ.IROW(IX)) OUT(IX:IX)=CHA(IFRA(IX)+23*(IFRA(IX)/12))
  150 IF(IR*(IROW(IX)-IR).GT.0) OUT(IX:IX)=CHA(35)
      WRITE(6,420) IR*DELY, IPOT, OUT(1:GLAST(OUT))
  160 CONTINUE
      IPOT=INT(LOG10(MAX(YMAX,-YMIN))+10.0001)-10
      OUT=' '
      DO 170 IX=1,NX
      IF(A(IS+18+IX).LT.-10.**(IPOT-4)) OUT(IX:IX)=CHA(39)
  170 IROW(IX)=10.**(3-IPOT)*ABS(A(IS+18+IX))+0.5
      WRITE(6,430) OUT(1:GLAST(OUT))
      DO 190 IR=4,1,-1
      DO 180 IX=1,NX
  180 OUT(IX:IX)=CHA(2+MOD(IROW(IX),10**IR)/10**(IR-1))
  190 WRITE(6,440) IPOT+IR-4, OUT(1:GLAST(OUT))
      IPOT=INT(LOG10(MAX(-A(IS+2),A(IS+3)-A(IS+4)))+10.0001)-10
      OUT=' '
      DO 200 IX=1,NX
      IF(A(IS+2)+(IX-1)*A(IS+4).LT.-10**(IPOT-3)) OUT(IX:IX)=CHA(39)
  200 IROW(IX)=10.**(2-IPOT)*ABS(A(IS+2)+(IX-1)*A(IS+4))+0.5
      WRITE(6,450) OUT(1:GLAST(OUT))
      DO 220 IR=3,1,-1
      DO 210 IX=1,NX
  210 OUT(IX:IX)=CHA(2+MOD(IROW(IX),10**IR)/10**(IR-1))
  220 WRITE(6,440) IPOT+IR-3, OUT(1:GLAST(OUT))
      ENDIF
      DO 230 IX=1,NX
      CTA=ABS(A(IS+18+IX))
      X=A(IS+2)+(IX-0.5)*A(IS+4)
      EV(1)=EV(1)+CTA
      EV(2)=EV(2)+CTA*X
  230 EV(3)=EV(3)+CTA*X**2
      SMALL=1E-20
      XMEAN=EV(2)/MAX(EV(1),SMALL)
      ZERO=0
      XRMS=SQRT(MAX(ZERO,EV(3)/MAX(EV(1),SMALL)-XMEAN**2))
      WRITE(6,460) INT(A(IS+9)+0.5),XMEAN,A(IS+13),A(IS+2),A(IS+14),
     &XRMS,A(IS+15),A(IS+3)
      ELSE
      ZMAX=35.4
      DO 240 IC=IS+19,IS+NX*NY+18
  240 IF(ABS(A(IC)).GT.ZMAX) ZMAX=1.0001*ABS(A(IC))
      ZFAC=26./LOG(ZMAX/9.5)
      DO 290 IY=1,MAX(NY+4,38)
      OUT=' '
      IF(IY.EQ.1.OR.IY.EQ.NY+2) THEN
      OUT(4:4)=CHA(38)
      DO 250 IX=5,2*NX+5
  250 OUT(IX:IX)=CHA(39)
      OUT(2*NX+6:2*NX+6)=CHA(38)
      ELSEIF(IY.EQ.NY+3) THEN
      DO 260 IX=6,2*NX+4,2
  260 OUT(IX:IX)=CHA(2+(IX-3)/20)
      ELSEIF(IY.EQ.NY+4) THEN
      DO 270 IX=6,2*NX+4,2
  270 OUT(IX:IX)=CHA(IX/2-10*((IX-3)/20))
      ELSEIF(IY.LE.NY+1) THEN
      OUT(1:1)=CHA(2+(NY+2-IY)/10)
      OUT(2:2)=CHA(4+NY-IY-10*((NY+2-IY)/10))
      OUT(4:4)=CHA(40)
      DO 280 IX=6,2*NX+4,2
      CT=A(IS+16+NX*(NY+1-IY)+IX/2)
      IF(CT.LE.-0.05) OUT(IX-1:IX-1)=CHA(39)
      CTA=ABS(CT)
      IF(CTA.GE.0.1) OUT(IX:IX)=CHA(2)
      IF(CTA.GE.0.5.AND.CTA.LT.36.) OUT(IX:IX)=CHA(2+INT(CTA+0.5))
  280 IF(ZMAX.GT.35.5.AND.CTA.GE.9.5) OUT(IX:IX)=
     &CHA(12+INT(ZFAC*LOG(CTA/9.5)))
      OUT(2*NX+6:2*NX+6)=CHA(40)
      ENDIF
      IF(IY.GE.2.AND.IY.LE.38) THEN
      WMIN=0.1
      IF(IY.GE.3) WMIN=IY-2.5
      IF(ZMAX.GT.35.5.AND.IY.GE.12) WMIN=9.5*EXP((IY-12)/ZFAC)
      ENDIF
      IF(IY.EQ.1) WRITE(6,470) OUT(1:GLAST(OUT))
      IF(IY.GT.1.AND.IY.LE.38)WRITE(6,480)CHA(IY),WMIN,OUT(1:GLAST(OUT))
  290 IF(IY.GE.39) WRITE(6,490) OUT(1:GLAST(OUT))
      DO 300 IY=1,NY
      Y=A(IS+6)+(IY-0.5)*A(IS+8)
      DO 300 IX=1,NX
      X=A(IS+2)+(IX-0.5)*A(IS+4)
      CTA=ABS(A(IS+18+NX*(IY-1)+IX))
      EV(1)=EV(1)+CTA
      EV(2)=EV(2)+CTA*X
      EV(3)=EV(3)+CTA*X**2
      EV(4)=EV(4)+CTA*Y
      EV(5)=EV(5)+CTA*Y**2
  300 EV(6)=EV(6)+CTA*X*Y
      SMALL=1E-20
      XMEAN=EV(2)/MAX(EV(1),SMALL)
      ZERO=0
      XRMS=SQRT(MAX(ZERO,EV(3)/MAX(EV(1),SMALL)-XMEAN**2))
      YMEAN=EV(4)/MAX(EV(1),SMALL)
      ZERO=0
      YRMS=SQRT(MAX(ZERO,EV(5)/MAX(EV(1),SMALL)-YMEAN**2))
      XYCOR=(EV(6)/MAX(EV(1),SMALL)-XMEAN*YMEAN)/MAX(SMALL,XRMS*YRMS)
      WRITE(6,500) INT(A(IS+9)+0.5),(A(IS+J), J=16,18), XMEAN, A(IS+2),
     &XRMS, A(IS+3), (A(IS+J), J=13,15), YMEAN, A(IS+6), YRMS, A(IS+7),
     &(A(IS+J), J=10,12), XYCOR
      ENDIF
      RETURN
  400 FORMAT(/5X,'HISTOGRAM NO',I4,' : NO ENTRIES')
  410 FORMAT('1'/5X,'HISTOGRAM NO',I4,4X,A60,'19',I2,'-',I2,'-',I2,
     &1X,A5/)
  420 FORMAT(3X,F7.2,'*10**',I2,4X,A)
  430 FORMAT(/9X,'CONTENTS',4X,A)
  440 FORMAT(10X,'*10**',I2,4X,A)
  450 FORMAT(/9X,'LOW EDGE',4X,A)
  460 FORMAT(/5X,'ENTRIES  =',I10,1P,6X,'MEAN =',E10.3,6X,'UNDERFLOW ='
     &,E10.3,6X,'LOW EDGE  =',E10.3/5X,'ALL CHAN =',E10.3,6X,
     &'RMS  =',E10.3,6X,'OVERFLOW  =',E10.3,6X,'HIGH EDGE =',E10.3)
  470 FORMAT(8X,'SCALE',6X,A)
  480 FORMAT(6X,A1,1X,F7.1,'-',3X,A)
  490 FORMAT(19X,A)
  500 FORMAT(/5X,'ENTRIES =',I10,1P,17X,E10.3,' ? ',E10.3,' ? ',E10.3,
     &7X,'XMEAN =',E10.3/5X,'XMIN    =',E10.3,17X,36('-'),7X,'XRMS  =',
     &E10.3/5X,'XMAX    =',E10.3,7X,'OVERFLOW',2X,E10.3,' ? ',E10.3,
     &' ? ',E10.3,7X,'YMEAN =',E10.3/5X,'YMIN    =',E10.3,17X,36('-'),
     &7X,'YRMS  =',E10.3/5X,'YMAX    =',E10.3,17X,E10.3,' ? ',E10.3,
     &' ? ',E10.3,7X,'XYCOR =',E10.3)
      END
C*********************************************************************
      INTEGER FUNCTION GLAST(STRING)
C  LOOK FOR LAST NON-SPACE IN STRING
      CHARACTER STRING*108
      I=108
 10   IF (STRING(I:I).EQ.' ') THEN
        I=I-1
        IF (I.GT.1) GOTO 10
      ENDIF
      GLAST=I
      END
C*********************************************************************
      SUBROUTINE GRESET(ID)
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NSIZE=200000,NMAX=2000)
      COMMON /GBOOK/ A(NSIZE)
      IF (ID.GT.NMAX) RETURN
      IS=A(ID+2)+0.5
      DO 110 IC=IS+9,IS+18+INT(A(IS+1)+0.5)*INT(A(IS+5)+0.5)
  110 A(IC)=0.
      RETURN
      END
C*********************************************************************
      BLOCK DATA GDATA
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NSIZE=200000,NMAX=2000)
      COMMON /GBOOK/ A(NSIZE)
      DATA (A(J), J=1,NMAX+2)/NMAX,2,NMAX*0/
      END
C***  END GBOOK  *****************************************************
C-----------------------------------------------------------------------
C  SOME NEW ROUTINES WRITTEN BY MIKE SEYMOUR:
C  GSCAT2(ID,X,Y,N)
C     EXACTLY THE SAME AS GFILL2(ID,X,Y,1.0) EXCEPT THAT ONLY THE
C     FIRST N CALLS ARE USED. VERY USEFUL FOR MAKING SCATTER PLOTS
C     THAT ALL HAVE THE SAME NUMBER OF ENTRIES.
C
C  GAREA(ID,AREA)
C     RESCALE HISTOGRAM SO THAT AREA UNDER IT BECOMES AREA
C
C  GACCUM(ID)
C     REPLACE HISTOGRAM BY A CUMULATIVE SUM OF ITS ENTRIES
C-----------------------------------------------------------------------
      SUBROUTINE GSCAT2(ID,X,Y,N)
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NSIZE=200000,NMAX=2000)
      COMMON /GBOOK/ A(NSIZE)
      IF (ID.GT.NMAX) RETURN
      ONE=1
      IS=A(ID+2)+0.5
      IF (INT(A(IS+9)+0.5).LT.N) CALL GFILL2(ID,X,Y,ONE)
      END
C-----------------------------------------------------------------------
      SUBROUTINE GAREA(ID,AREA)
C
C  SCALES HISTOGRAM ID SO THAT THE AREA UNDER IT BECOMES AREA
C
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NSIZE=200000,NMAX=2000)
      COMMON /GBOOK/ A(NSIZE)
      IF (ID.GT.NMAX) RETURN
      IS=NINT(A(ID+2))
      BNWDTH=A(IS+4)
      AREAOL=A(IS+14)
      IF (AREAOL*BNWDTH.EQ.0.0) RETURN
      AREASC=AREA/(AREAOL*BNWDTH)
      CALL GSCALE (ID,AREASC)
      END
C-----------------------------------------------------------------------
      SUBROUTINE GACCUM(ID)
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NSIZE=200000,NMAX=400)
      COMMON /GBOOK/ A(NSIZE)
      IF (ID.GT.NMAX) RETURN
      TOTAL=0.0
      IS=NINT(A(ID+2))
      IF (IS.EQ.0) RETURN
      DO 50 IX=1,NINT(A(IS+1))
        TOTAL=TOTAL+A(IS+18+IX)
        A(IS+18+IX)=TOTAL
 50   CONTINUE
      END
C-----------------------------------------------------------------------
      SUBROUTINE GTOPDR(ID,NEW,HIST,LOG)
C
C  M SEYMOUR'S GBOOK->TOPDRAWER INTERFACE
C  WRITES HISTOGRAM WITH TOPDRAWER PACKAGE
C
C  NEW  IS .TRUE.  TO PUT HISTOGRAM ON A NEW PAGE
C          .FALSE. TO PUT IT ON TOP OF THE LAST ONE
C
C  HIST IS .TRUE.  FOR HISTOGRAM
C          .FALSE. FOR GRAPH
C
C  LOG  IS .TRUE.  FOR LOGARITHMIC Y-SCALE
C          .FALSE. FOR LINEAR Y-SCALE
C
C  .TRUE.  MAY BE ABBREVIATED TO 1 WHEN CALLING, AND
C  .FALSE. TO 0
C
C  IF LOG IS SELECTED, THEN ALL ZERO OR NEGATIVE BINS ARE SET TO
C  ONE THOUSANTH OF THE LOWEST OTHER BIN
C
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NSIZE=200000,NMAX=2000)
      COMMON /GBOOK/ A(NSIZE)
      LOGICAL NEW,HIST,LOG
      CALL GTOPER(ID,NEW,HIST,LOG,0)
      END
C-----------------------------------------------------------------------
      SUBROUTINE GTOPER(ID,NEW,HIST,LOG,IDERR)
C
C  IDENTICAL TO GTOPDR EXCEPT THAT IF IDERR IS NON-ZERO IT IS USED FOR
C  ERROR BARS FOR 1-DIM HISTOGRAMS
C
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NSIZE=200000,NMAX=2000)
      COMMON /GBOOK/ A(NSIZE)
      LOGICAL NEW,HIST,LOG
      IF (ID.GT.NMAX) RETURN
      IS=NINT(A(ID+2))
      IF (IS .EQ. 0) RETURN
      IF (A(IS+14) .EQ. 0.0) THEN
         WRITE (6,1000) ID
      ELSE
         WRITE (6,1010) ID
         IF (A(IS+5) .LT. 1.5) THEN
            CALL GTOP1(ID,NEW,HIST,LOG,IDERR)
         ELSE
            CALL GTOP2(ID,NEW)
         ENDIF
      ENDIF
 1000 FORMAT (/5X,'HISTOGRAM NO',I4,' : NO ENTRIES FOUND')
 1010 FORMAT (/5X,'HISTOGRAM NO',I4,' : OUTPUTTING TO TOPDRAWER FILE')
      END
C-----------------------------------------------------------------------
      SUBROUTINE GTOP1(ID,NEW,HIST,LOG,IDERR)
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NSIZE=200000,NMAX=2000)
      COMMON /GBOOK/ A(NSIZE)
      LOGICAL NEW,HIST,LOG, DIV
      CHARACTER*8 TYPE(4),PLOT
      COMMON /GTCOM/ TYPE,NTYPE
      PARAMETER (NXMAX=1000)
      DIMENSION X(NXMAX),Y(NXMAX),E(NXMAX)
      LOGICAL OPEN
      CHARACTER TITLE*60
      COMMON /GFILE/ OPEN, TITLE
      IF (ID.GT.NMAX) RETURN
      IF (NEW) THEN
            CALL GCLOSEF
            CALL GTITLE(ID,TITLE,LENTIT)
            CALL GOPENF
      ENDIF
      WEIGHT=0.0
      AVGE=0.0
      VARIAN=0.0
      IS=NINT(A(ID+2))
      NX=NINT(A(IS+1))
      IF (NX.GT.NXMAX) THEN
        WRITE (6,2000) NXMAX
        NX=NXMAX
      ENDIF
      ISERR=0
      IF (IDERR.NE.0) THEN
         ISERR=NINT(A(IDERR+2))
         IF (ISERR .NE. 0) THEN
            IF (A(ISERR+5) .GT. 1.5) THEN
               WRITE (6,1000) IDERR
               ISERR=0
            ELSEIF (NX.NE.NINT(A(ISERR+1))) THEN
               WRITE (6,1010) IDERR
               ISERR=0
            ELSE
               WRITE (6,1020) IDERR
            ENDIF
         ENDIF
      ENDIF
      YMIN=1D30
      DO 50  IX=1, NX
 50      IF (A(IS+18+IX).GT.0.0) YMIN=MIN(YMIN,A(IS+18+IX))
      DIV=.FALSE.
      DO 100 IX=1, NX
         X(IX)=(IX-0.5)*A(IS+4)+A(IS+2)
         Y(IX)=A(IS+18+IX)
         E(IX)=0
         IF (ISERR.NE.0) E(IX)=A(ISERR+18+IX)
         WEIGHT=WEIGHT+Y(IX)
         AVGE=AVGE+X(IX)*Y(IX)
         VARIAN=VARIAN+X(IX)**2*Y(IX)
         IF (LOG .AND. Y(IX).LE.0.0) THEN
            WRITE (6,2010) Y(IX),YMIN/1000
            Y(IX)=YMIN/1000
            DIV=.TRUE.
         ENDIF
 100  CONTINUE
      IF (WEIGHT.EQ.0) WEIGHT=1
      AVGE=AVGE/WEIGHT
      VARIAN=VARIAN/WEIGHT-AVGE**2
C  SET UP PAGE...
      IF (NEW) THEN
         WRITE (21,*) 'NEW FRAME'
         WRITE (21,*) 'SET WINDOW X 4.9 9.5 Y 2 9'
         WRITE (21,*) 'SET FONT DUPLEX'
         IF (LOG) WRITE (21,*) 'SET SCALE Y LOG'
         IF (DIV) WRITE (21,*) 'SET LIMIT YMIN',YMIN/10
         WRITE (21,*) 'TITLE TOP ''',TITLE(1:LENTIT),''''
         IF (NINT(A(IS+2)*10).EQ.0.AND.NINT(A(IS+3)*10).EQ.10) THEN
           WRITE (21,*) 'SET LIMITS X 0 1'
         ELSE
           WRITE (21,*) 'SET LIMITS X',SNGL(A(IS+2)),SNGL(A(IS+3))
         ENDIF
         WRITE (21,*) 'SET PATT 0.02 0.08'
         WRITE (21,*) 'SET ORDER X Y DY'
         NTYPE=0
      ENDIF
      NTYPE=MOD(NTYPE,4)+1
C  PLOT...
      DO 200 I=1, NX
         IF ((ABS(X(I)).GE.1E-3.AND.ABS(X(I)).LT.1E5.OR.X(I).EQ.0).AND.
     &       (ABS(Y(I)).GE.1E-5.AND.ABS(Y(I)).LT.1E3.OR.Y(I).EQ.0)) THEN
            WRITE (21,'(F14.5,2F14.7)') X(I),Y(I),E(I)
         ELSE
            WRITE (21,'(3(1PE14.4))') X(I),Y(I),E(I)
         ENDIF
 200  CONTINUE
      PLOT=' '
      IF (ISERR.NE.0) PLOT=' ; PLOT'
      IF (HIST) THEN
         WRITE (21,*) 'HIST ',TYPE(NTYPE),PLOT
      ELSE
         WRITE (21,*) 'JOIN ',TYPE(NTYPE),PLOT
      ENDIF
      WRITE (6,1100) A(IS+13)
      WRITE (6,1110) A(IS+14)
      WRITE (6,1120) A(IS+15)
      WRITE (6,1130) AVGE
      IF (VARIAN.GE.0.0) THEN
         WRITE (6,1140) SQRT(VARIAN)
      ELSE
         WRITE (6,1150) VARIAN
      ENDIF
 1000 FORMAT (5X,'HISTOGRAM NO',I4,' : SHOULD BE 1-DIM BUT ISN''T')
 1010 FORMAT (5X,'HISTOGRAM NO',I4,' : SHOULD BE SAME SIZE BUT ISN''T')
 1020 FORMAT (5X,'HISTOGRAM NO',I4,' : BEING USED FOR ERROR BARS')
 1100 FORMAT (5X,'Total underflow=',G10.3)
 1110 FORMAT (5X,'Total entry    =',G10.3)
 1120 FORMAT (5X,'Total overflow =',G10.3)
 1130 FORMAT (5X,'Mean value        =',G10.3)
 1140 FORMAT (5X,'Standard deviation=',G10.3)
 1150 FORMAT (5X,'Variance          =',G10.3)
 2000 FORMAT (5X,'WARNING: Histogram too big, using first',I4,' bins')
 2010 FORMAT (5X,'WARNING:',G10.3,' cannot be logged. Now using',G10.3)
      END
C-----------------------------------------------------------------------
      SUBROUTINE GTOP2(ID,NEW)
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NSIZE=200000,NMAX=2000)
      COMMON /GBOOK/ A(NSIZE)
      LOGICAL NEW
      PARAMETER (NZ=250)
      DIMENSION Z(NZ,NZ)
      LOGICAL OPEN
      CHARACTER TITLE*60
      COMMON /GFILE/ OPEN,TITLE
      REAL RAN
      DATA ISEED/12345/
      IF (ID.GT.NMAX) RETURN
      IF (NEW) THEN
            CALL GCLOSEF
            CALL GTITLE(ID,TITLE,LENTIT)
            CALL GOPENF
      ENDIF
      WEIGHT=0.0
      XAVGE=0.0
      XVARIN=0.0
      YAVGE=0.0
      YVARIN=0.0
      ZMAX=0
      IS=NINT(A(ID+2))
      NX=NINT(A(IS+1))
      NY=NINT(A(IS+5))
      IF (NX.GT.NZ) THEN
        WRITE (6,2000) NZ
        NX=NZ
      ENDIF
      IF (NY.GT.NZ) THEN
        WRITE (6,2000) NZ
        NY=NZ
      ENDIF
      DO 100 IX=1, NX
      DO 100 IY=1, NY
         X=(IX-1)*A(IS+4)+A(IS+2)
         Y=(IY-1)*A(IS+8)+A(IS+6)
         Z(IX,IY)=A(IS+18+NX*(IY-1)+IX)
         SMALL=1E-20
         IF (ABS(Z(IX,IY)).LT.SMALL) Z(IX,IY)=0.0
         IF (Z(IX,IY).GT.ZMAX) ZMAX=Z(IX,IY)
         WEIGHT=WEIGHT+Z(IX,IY)
         XAVGE=XAVGE+X*Z(IX,IY)
         XVARIN=XVARIN+X**2*Z(IX,IY)
         YAVGE=YAVGE+Y*Z(IX,IY)
         YVARIN=YVARIN+Y**2*Z(IX,IY)
 100  CONTINUE
      IF (ZMAX.GE.1 .AND. ZMAX.LE.20) THEN
         ZMAX=20
      ELSE
         WRITE (6,2010) ZMAX/20
      ENDIF
      XAVGE=XAVGE/WEIGHT
      XVARIN=XVARIN/WEIGHT-XAVGE**2
      YAVGE=YAVGE/WEIGHT
      YVARIN=YVARIN/WEIGHT-YAVGE**2
C  SET UP PAGE...
      IF (NEW) THEN
        WRITE (21,*) 'NEW FRAME'
        IF (A(IS+2).EQ.A(IS+6) .AND. A(IS+3).EQ.A(IS+7))
     &  WRITE (21,*) 'SET WINDOW X 3.7 10.7 Y 2 9'
        WRITE (21,*) 'SET FONT DUPLEX'
        WRITE (21,*) 'TITLE TOP ''',TITLE(1:LENTIT),''''
        WRITE (21,*) 'SET LIMITS X',SNGL(A(IS+2)),SNGL(A(IS+3)),
     &    ' Y',SNGL(A(IS+6)),SNGL(A(IS+7))
      ENDIF
C  PLOT...
      NN=0
      DO 210 IX=1, NX
      DO 210 IY=1, NY
         DO 200 IZ=1, NINT(20*Z(IX,IY)/ZMAX)
            X=(IX-RAN(ISEED))*A(IS+4)+A(IS+2)
            Y=(IY-RAN(ISEED))*A(IS+8)+A(IS+6)
            IF ((ABS(X).GE.1E-5.AND.ABS(X).LT.1E3.OR.X.EQ.0).AND.
     &          (ABS(Y).GE.1E-5.AND.ABS(Y).LT.1E3.OR.Y.EQ.0)) THEN
               WRITE (21,'(2F14.7)') X,Y
            ELSE
               WRITE (21,'(2(1PE14.4))') X,Y
            ENDIF
            NN=NN+1
 200     CONTINUE
         IF (NN.GT.500) THEN
            WRITE (21,*) 'PLOT'
            NN=0
          ENDIF
 210  CONTINUE
      IF (NN.GT.0) WRITE (21,*) 'PLOT'
      WRITE (6,1000)
      WRITE (6,1030) (A(IS+J), J=16,18)
      WRITE (6,1020) (A(IS+J), J=13,15)
      WRITE (6,1010) (A(IS+J), J=10,12)
      WRITE (6,1040) XAVGE
      IF (XVARIN.GT.0) WRITE (6,1050) SQRT(XVARIN)
      WRITE (6,1060) YAVGE
      IF (YVARIN.GT.0) WRITE (6,1070) SQRT(YVARIN)
 1000 FORMAT (20X,'   X under',5X,' X on grid',5X,'    X over')
 1010 FORMAT (5X,'   Y under',3(5X,G10.3))
 1020 FORMAT (5X,' Y on grid',3(5X,G10.3))
 1030 FORMAT (5X,'    Y over',3(5X,G10.3))
 1040 FORMAT (5X,'Mean x value        =',G10.3)
 1050 FORMAT (5X,'X standard deviation=',G10.3)
 1060 FORMAT (5X,'Mean y value        =',G10.3)
 1070 FORMAT (5X,'Y standard deviation=',G10.3)
 2000 FORMAT (5X,'WARNING: Histogram too big, using first',I4,' bins')
 2010 FORMAT (5X,'WARNING: Dividing by ',G10.3)
      END
C-----------------------------------------------------------------------
      SUBROUTINE GTOP
C---OUTPUT ALL HISTOGRAMS TO TOPDRAWER
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NSIZE=200000,NMAX=2000)
      COMMON /GBOOK/ A(NSIZE)
      DO 100 ID=1,INT(A(1)+0.5)
      IS=A(ID+2)+0.5
      IF(IS.EQ.0.OR.A(IS+9).LT.0.5) GOTO 100
      CALL GTOPDR(ID,.TRUE.,.TRUE.,.FALSE.)
  100 CONTINUE
      END
C-----------------------------------------------------------------------
      SUBROUTINE GTITLE(ID,TITLE,LENTIT)
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NSIZE=200000,NMAX=2000)
      COMMON /GBOOK/ A(NSIZE)
      CHARACTER TITLE*60
      EQUIVALENCE (REQ,IEQ)
      IF (ID.GT.NMAX) RETURN
      IS=NINT(A(ID+2))
      NX=NINT(A(IS+1))
      NY=NINT(A(IS+5))
      DO 100 IT=1, 20
         REQ=A(IS+18+NX*NY+IT)
         TITLE(3*IT-2:3*IT)=CHAR(IEQ/256**2)//CHAR(MOD(IEQ,256**2)
     &      /256)//CHAR(MOD(IEQ,256))
 100  CONTINUE
      LENTIT=60
 200  IF (LENTIT.GT.1 .AND. TITLE(LENTIT:LENTIT).EQ.' ') THEN
         LENTIT=LENTIT-1
         GOTO 200
      ENDIF
      END
C-----------------------------------------------------------------------
      SUBROUTINE GOPENF
C---OPEN TopDrawer FILE
      LOGICAL OPEN
      CHARACTER TITLE * 60
      CHARACTER FNAME * 74
      CHARACTER*40 OPATH
      COMMON /GFILE/ OPEN, TITLE
      COMMON /FPATH/ OPATH
      DATA OPEN/.FALSE./
      IF (OPEN) RETURN
      OPEN=.TRUE.
      N=0
      FNAME=TRIM(OPATH)//'gtopdraw_'//TRIM(TITLE)//'.top'
      OPEN (21,FILE=TRIM(FNAME),STATUS='NEW',ERR=10)
      WRITE (6,*) 'Using file "'//TRIM(FNAME)//'"'
      RETURN
 10   N=N+1
      IF (N.GE.100) STOP
      WRITE (6,*) 'Could not open file "',TRIM(FNAME),'"'
      WRITE (FNAME,'(A,I2.2,A4)') TRIM(OPATH)//'gtopdraw_'
     $      //TRIM(TITLE)//'_',N,'.top'
      OPEN (21,FILE=TRIM(FNAME),STATUS='NEW',ERR=10)
      WRITE (6,*) 'Using file "',TRIM(FNAME),'" instead'
      END
C-----------------------------------------------------------------------
      SUBROUTINE GCLOSEF
C---CLOSE TopDrawer FILE
      LOGICAL OPEN
      CHARACTER TITLE * 60
      COMMON /GFILE/ OPEN,TITLE
      IF (.NOT.OPEN) RETURN
      CLOSE(21)
      OPEN = .FALSE.
      END
C-----------------------------------------------------------------------
      BLOCK DATA GTDAT
      IMPLICIT INTEGER (I-N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER*8 TYPE(4)
      COMMON /GTCOM/ TYPE,NTYPE
      DATA TYPE,NTYPE/'        ','DASH    ','DOT-DASH','PATTERN ',0/
      END
C-----------------------------------------------------------------------
      double precision function ran(iseed)
      implicit none
      integer iseed
      double precision r(1)
      call rangen(1,r)
      iseed=iseed+1
      ran=r(1)
      end
C-----------------------------------------------------------------------