*
*=== Gplevbin =========================================================*
*
      MODULE GPLMOD

         REAL,DIMENSION(:),ALLOCATABLE :: FLUMAT, FLEMAT
         DOUBLE PRECISION,DIMENSION(:,:),ALLOCATABLE :: GMRZTA, VLRZTA,
     &                                                  GERZTA

      END MODULE GPLMOD

      PROGRAM GPLEVBIN
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 2003-2019:  CERN & INFN                            *
*                                                                      *
*                   2020-    :  CERN                                   *
*                                                                      *
*     All Rights Reserved.                                             *
*                                                                      *
*----------------------------------------------------------------------*
*
      USE GPLMOD

      PARAMETER ( NM     = 50000 )
*     PARAMETER ( NMMIN  = 5000  )
*     PARAMETER ( EPS    = 0.25  )
*     PARAMETER ( EPS2   = 0.0625)
      CHARACTER YES*1,FILE1*250,FILE3*250,FILE4*250,
     &          CARD*80,TITPLO*80,FILCMP*80
*     CHARACTER CHELP*20
      DIMENSION XX(NM),YY(NM),XXX(NM),YYY(NM)
      DOUBLE PRECISION XPOINT (NM), YPOINT (NM), ZPOINT (NM), WHAT (6)
      DOUBLE PRECISION WX0, WY0, WZ0, WX1, WY1, WZ1, WXX, WXY, WXZ,
     &                 WYX, WYY, WYZ, WXAXLN, WYAXLN
*     COMMON /LGCURL/NUM,IVET(MXXXXX)
      LOGICAL LCMPFL, LCMPOU, LINSID, LINPRE, LSMPOI, LGEOPL
      LOGICAL LRHOFL, LROTA , LAXSWP, LGAXSW
      LOGICAL LRPHI

      DOUBLE PRECISION X, Y, XM, DDX, DXBIN, DYBIN

      CHARACTER SDUM*8
*
      CALL CMSPPR
      CALL RTMTXI
*     NUM   = MXXXXX
      FACTX = 1.E+00
      FACTY = 1.E+00
      GGXM  = 2047.0
      GGYM  = 2047.0
      WX0   = 0.D+00
      WY0   = 0.D+00
      WZ0   = 0.D+00
      WXX   = 0.D+00
      WXY   = 0.D+00
      WXZ   = 0.D+00
      WYX   = 0.D+00
      WYY   = 0.D+00
      WYZ   = 0.D+00
      LSMPOI=.FALSE.
      LGEOPL=.FALSE.
      LCMPFL=.FALSE.
      LCMPOU=.FALSE.
      LROTA =.FALSE.
      WRITE(*,*)' PLOTGEOM file (def. PLOTGEOM.STORE):'
      READ(*,'(A50)')FILE1
      IF ( FILE1 .EQ. '      ' ) FILE1 = 'PLOTGEOM.STORE'
      WRITE (*,'(A)') ' Swap plotgeom/bin axis (def=n)?'
      READ  (*,'(A)') YES
      LAXSWP = YES .EQ. 'Y' .OR. YES .EQ. 'y'
      LGAXSW = LAXSWP
      IF ( LAXSWP ) THEN
         WRITE (*,'(A)')
     &' Swap: both geo and bin axis (def=0), geo only (1) bin only (2)?'
         READ (*,'(A)') YES
         IFSWAP = 0
         IF ( YES .NE. ' ' ) READ (YES,*) IFSWAP
         IF ( IFSWAP .LT. 0 .OR. IFSWAP .GT. 2 ) IFSWAP = 0
         LAXSWP = MOD (IFSWAP,2) .EQ. 0
         LGAXSW = IFSWAP .LE. 1
      END IF
*
      OPEN ( UNIT=99, FILE=FILE1,FORM='UNFORMATTED',
     &       STATUS='OLD', ERR=133 )
*
      LGEOPL=.TRUE.
*  New for (geometry) rotations:
      WRITE (*,*)' Rotate geometry ? (y,n), def n'
      READ (*,'(A)')YES
      IF ( YES .EQ. 'y' .OR. YES .EQ. 'Y' ) THEN
         LROTA =.TRUE.
         WRITE(*,'(A)')' Enter the 6 rot-defi what''s'
         READ (*,*) (WHAT(I),I=1,6)
         SDUM=' '
         CALL SETROT (WHAT,SDUM)
         CALL PRNROT
         KROTAT = 1
      ELSE
         KROTAT = 0
         LROTA=.FALSE.
      END IF
*
*     WRITE(*,*)' Type the "X" axis expansion factor (def=1):'
*     READ (*,'(A)') CHELP
*     IF ( CHELP .NE. '  ' ) READ (CHELP,*) FACTX
*     WRITE(*,*) FACTX
*     WRITE(*,*)' Type the "Y" axis expansion factor (def=1):'
*     READ (*,'(A)') CHELP
*     IF ( CHELP .NE. '  ' ) READ (CHELP,*) FACTY
*     WRITE(*,*) FACTY
      OPEN (UNIT=99,FILE=FILE1,FORM='UNFORMATTED',
     1      STATUS='OLD')
      READ(99)TITPLO
      IF ( TITPLO .EQ. ' ***COMPRESSED***COMPRESSED*** ' ) THEN
         WRITE (*,*) ' !!! Compressed input file !!! '
         LCMPFL = .TRUE.
         LCMPOU = .FALSE.
         READ (99) TITPLO
      ELSE
         LCMPFL = .FALSE.
         WRITE (*,*)' (Possible) compressed output file (def='' ''=no)'
         READ (*,'(A)')FILCMP
         IF ( FILCMP .NE. '   ' ) THEN
            LCMPOU = .TRUE.
* Start_VAX_seq
*           OPEN ( UNIT = 98, FILE=FILCMP, STATUS='NEW',
*    &             FORM ='UNFORMATTED' )
* End_VAX_seq
* Start_UNIX_seq
            OPEN ( UNIT = 98, FILE=FILCMP, STATUS='UNKNOWN',
     &             FORM ='UNFORMATTED' )
* End_UNIX_seq
            WRITE(98)' ***COMPRESSED***COMPRESSED*** '
            WRITE(98) TITPLO
         ELSE
            LCMPOU=.FALSE.
         END IF
      END IF
      WRITE(*,'(A80)')TITPLO
      READ(99)X0,Y0,Z0,X1,Y1,Z1,XU,YU,ZU,XV,YV,ZV,XAXLEN,YAXLEN
      IF (LCMPOU)
     &WRITE(98)X0,Y0,Z0,X1,Y1,Z1,XU,YU,ZU,XV,YV,ZV,XAXLEN,YAXLEN
* X0,Y0,Z0 is the left-bottom corner, X1,Y1,Z1 is the right-top one
* Xu,Yu,Zu is the bottom-top axis versor (Y axis just to make
* confusion...), the length is Yaxlen
* Xv,Yv,Zv is the left-right axis versor (X axis), the length is Xaxlen
* Origins of the R (X) and T (Y) axis:
      T00=X0*XU+Y0*YU+Z0*ZU
      R00=X0*XV+Y0*YV+Z0*ZV
      EPSY  = 0.99999E+00 * YAXLEN / 2000.E+00
      EPSX  = 0.99999E+00 * XAXLEN / 2000.E+00
      DXEPS = 0.13E+00 * EPSX
      DYEPS = 0.13E+00 * EPSY
*+++++ Debug:
      WRITE(*,*)' R00,T00',R00,T00
*-----
      WX0=X0
      WY0=Y0
      WZ0=Z0
      WX1=X1
      WY1=Y1
      WZ1=Z1
      WYX=XU
      WYY=YU
      WYZ=ZU
      WXX=XV
      WXY=YV
      WXZ=ZV
      WXAXLN =XAXLEN
      WYAXLN =YAXLEN
* X0,Y0,Z0 is now the left-top corner:
      X0 =WX0+WYAXLN*WYX
      Y0 =WY0+WYAXLN*WYY
      Z0 =WZ0+WYAXLN*WYZ
* X1,Y1,Z1 is now the right-bottom corner:
      X1 =WX0+WXAXLN*WXX
      Y1 =WY0+WXAXLN*WXY
      Z1 =WZ0+WXAXLN*WXZ
* Xu,Yu,Zu is now the top-bottom axis versor:
      XU =-XU
      YU =-YU
      ZU =-ZU
      IF ( LROTA ) THEN
         XPOINT (1) = X0
         YPOINT (1) = Y0
         ZPOINT (1) = Z0
         XPOINT (2) = X1
         YPOINT (2) = Y1
         ZPOINT (2) = Z1
         NPOINT = 2
         CALL DOTRSF ( NPOINT, XPOINT, YPOINT, ZPOINT, KROTAT )
         X0 = XPOINT (1)
         Y0 = YPOINT (1)
         Z0 = ZPOINT (1)
         X1 = XPOINT (2)
         Y1 = YPOINT (2)
         Z1 = ZPOINT (2)
         XPOINT (1) = XU
         YPOINT (1) = YU
         ZPOINT (1) = ZU
         XPOINT (2) = XV
         YPOINT (2) = YV
         ZPOINT (2) = ZV
         NPOINT = 2
         CALL DORTNO ( NPOINT, XPOINT, YPOINT, ZPOINT, KROTAT )
         XU = XPOINT (1)
         YU = YPOINT (1)
         ZU = ZPOINT (1)
         XV = XPOINT (2)
         YV = YPOINT (2)
         ZV = ZPOINT (2)
      END IF
*+++++ Debug:
      WRITE(*,*)X0,Y0,Z0
      WRITE(*,*)X1,Y1,Z1
      WRITE(*,*)XU,YU,ZU
      WRITE(*,*)XV,YV,ZV
*-----
* End reshaffling of input positions/directions:
      T0=-(X0*XU+Y0*YU+Z0*ZU)
      T1=-(X1*XU+Y1*YU+Z1*ZU)
      R0=X0*XV+Y0*YV+Z0*ZV
      R1=X1*XV+Y1*YV+Z1*ZV
      XUMIN = MIN (R0,R1)
      XUMAX = MAX (R0,R1)
      YUMIN = MIN (T0,T1)
      YUMAX = MAX (T0,T1)
      IF ( LGAXSW ) THEN
         WRITE (*,*) YUMIN, YUMAX
         WRITE (*,*) ' X_min - X_max (def. provided)'
         READ  (*,'(A)') CARD
         IF ( CARD .NE. ' ' ) READ (CARD,*) YUMIN,YUMAX
         WRITE (*,*) XUMIN, XUMAX
         WRITE (*,*) ' Y_min - Y_max (def. provided)'
         READ  (*,'(A)') CARD
         IF ( CARD .NE. ' ' ) READ (CARD,*) XUMIN,XUMAX
      ELSE
         WRITE (*,*) XUMIN, XUMAX
         WRITE (*,*) ' X_min - X_max (def. provided)'
         READ  (*,'(A)') CARD
         IF ( CARD .NE. ' ' ) READ (CARD,*) XUMIN,XUMAX
         WRITE (*,*) YUMIN, YUMAX
         WRITE (*,*) ' Y_min - Y_max (def. provided)'
         READ  (*,'(A)') CARD
         IF ( CARD .NE. ' ' ) READ (CARD,*) YUMIN,YUMAX
      END IF
      DX = XUMAX-XUMIN
      DY = YUMAX-YUMIN
*++++++++++
      XWMIN=XUMIN-0.05E+00*DX
      YWMIN=YUMIN-0.05E+00*DY
      DX = XUMAX - XWMIN
      DY = YUMAX - YWMIN
*----------
      G0XM=200.
      G0YM=200.
      GMM=1947.
      IF (FACTX*DX .GE. FACTY*DY) THEN
         GXM=GMM
         GYM=DY*FACTY/DX/FACTX*(GXM-G0XM)+G0XM
      ELSE
         GYM=GMM
         GXM=DX*FACTX/DY/FACTY*(GYM-G0YM)+G0YM
      END IF
      G1XM = G0XM + (XUMIN-XWMIN) * (GXM-G0XM) / (XUMAX-XWMIN)
      G1YM = G0YM + (YUMIN-YWMIN) * (GYM-G0YM) / (YUMAX-YWMIN)
      IF ( LAXSWP .NEQV. LGAXSW ) THEN
         ZTALOW = YUMIN
         ZTAHGH = YUMAX
         RADLOW = XUMIN
         RADHGH = XUMAX
      ELSE
         ZTALOW = XUMIN
         ZTAHGH = XUMAX
         RADLOW = YUMIN
         RADHGH = YUMAX
      END IF
      ZMIN = 0.E+00
      ZMAX = 0.E+00
 133  CONTINUE
      DO 50 ID=1,100
      IPEN=1
      WRITE(*,*)' Bin file:'
      READ(*,'(A)')FILE3
      IF (FILE3 .EQ. '  ' .AND. ID .GT. 1 ) GO TO 60
      IF (FILE3 .EQ. '  ') GO TO 55
      WRITE(*,*)
     & ' Density file (name) or density value (pos.) or norm. (neg.):'
      READ(*,'(A)')FILE4
      OPEN (UNIT=11,FILE=FILE3,STATUS='OLD',FORM='UNFORMATTED'
* Start_VAX_seq
*    &     ,READONLY
* End_VAX_seq
     &     )
      IF (FILE4 .EQ. '  ') THEN
         LRHOFL = .FALSE.
         RHOFIX = 1.E+00
      ELSE
         RHOFIX =-1.E+00
         LRHOFL = .TRUE.
         OPEN (UNIT=12,FILE=FILE4,STATUS='OLD',FORM='UNFORMATTED',
* Start_VAX_seq
*    &         READONLY,
* End_VAX_seq
     &         ERR=53)
         GO TO 54
   53    CONTINUE
         READ(FILE4,*)RHOFIX
         IF ( RHOFIX .LT. 0.E+00 ) THEN
            LRHOFL = .FALSE.
            RHOFIX =-RHOFIX
         END IF
   54    CONTINUE
      END IF
      RHOTH=1.E-30
      WRITE(*,*)' Threshold density (def. provided):'
      READ (*,'(A)') CARD
      IF ( CARD .NE. ' ' ) READ (CARD,*) RHOTH
      IF (RHOTH .LE. 0.E+00) RHOTH=5.E-3
      WRITE (*,*) ' Rho_threshold: ', RHOTH
      WRITE (*,*) ' Which binning?'
      READ  (*,*) NBG
*     NBG=1
      CALL UBDOS2(NZTABN,NRADBN,ZTALOW,ZTAHGH,LGEOPL,
     &            RADLOW,RADHGH,RHOTH,NBG,LRHOFL,RHOFIX,LRPHI)
      CLOSE (UNIT=11)
      IF (LRHOFL) CLOSE (UNIT=12)
* Start_VAX_seq
*     OPEN ( UNIT=93, FILE='gplevh.dat', STATUS ='NEW' )
* End_VAX_seq
* Start_UNIX_seq
      OPEN ( UNIT=93, FILE='gplevh.dat', STATUS ='UNKNOWN' )
* End_UNIX_seq
*
      DXBIN = (DBLE(ZTAHGH)-DBLE(ZTALOW))/DBLE(NZTABN)
      ZMIN  = +1.E+30
      ZMAX  = -1.E+30
*  +-------------------------------------------------------------------*
*  | 2D histogram
      IF (NRADBN .GT. 1) THEN

         DYBIN = (DBLE(RADHGH)-DBLE(RADLOW))/DBLE(NRADBN)

*        If LRPHI find multiples of lines to dump for simulating a R-Phi plot
         IF ( LRPHI .AND. NZTABN.LT.32) THEN
             MULT = INT(32.D+00 / DBLE(NZTABN))   ! +/- mod/2 from the 32
             DDX  = DXBIN / DBLE(MULT)
         ELSE
             MULT = 1
             DDX  = 0.D+00
         ENDIF

         DO IE=1,NZTABN
            X=DBLE(ZTALOW)+(DBLE(IE)-1.0D+00) * DXBIN

            DO IM=1,MULT
               IF ( IM .EQ. 1 ) THEN
                   XM = X
               ELSE
                   XM = X + DBLE(IM-1)*DDX
               END IF

               DO IA=1,NRADBN
                 Y    = DBLE(RADLOW)+(DBLE(IA)-1.0D+00)*DYBIN
                 IDOS = (IA-1)*NZTABN + IE
                 IF ( LAXSWP ) THEN
                    WRITE (93,'(1P,2E18.10,2E13.5)')
     &                     Y, XM, FLUMAT(IDOS), FLEMAT(IDOS)*100.0
                 ELSE
                    WRITE (93,'(1P,2E18.10,2E13.5)')
     &                     XM, Y, FLUMAT(IDOS), FLEMAT(IDOS)*100.0
                 END IF

                 IF(ABS(FLUMAT(IDOS)).LT.ZMIN.AND.ABS(FLUMAT(IDOS))
     &              .GT.2.E-38) ZMIN=ABS(FLUMAT(IDOS))
                 IF(ABS(FLUMAT(IDOS)).GT.ZMAX) ZMAX=ABS(FLUMAT(IDOS))
               END DO

               IDOS = (NRADBN-1)*NZTABN + IE
               IF ( LAXSWP ) THEN
                  WRITE (93,'(1P,2E18.10,2E13.5)')
     &                   Y+DYBIN, XM, FLUMAT(IDOS), FLEMAT(IDOS)*100.0
               ELSE
                  WRITE (93,'(1P,2E18.10,2E13.5)')
     &                   XM, Y+DYBIN, FLUMAT(IDOS), FLEMAT(IDOS)*100.0
               END IF
               WRITE (93,*)
            END DO
         END DO

         X = X + DXBIN
         DO IA=1,NRADBN
            Y    = DBLE(RADLOW)+(DBLE(IA)-1.0D+00)*DYBIN
            IDOS = (IA-1)*NZTABN + NZTABN
            IF ( LAXSWP ) THEN
               WRITE (93,'(1P,2E18.10,2E13.5)')
     &                Y, X, FLUMAT(IDOS), FLEMAT(IDOS)*100.0
            ELSE
               WRITE (93,'(1P,2E18.10,2E13.5)')
     &                X, Y, FLUMAT(IDOS), FLEMAT(IDOS)*100.0
            END IF
         END DO

         IDOS = (NRADBN-1)*NZTABN + NZTABN
         IF ( LAXSWP ) THEN
            WRITE (93,'(1P,2E18.10,2E13.5)')
     &             Y+DYBIN, X, FLUMAT(IDOS), FLEMAT(IDOS)*100.0
         ELSE
            WRITE (93,'(1P,2E18.10,2E13.5)')
     &             X, Y+DYBIN, FLUMAT(IDOS), FLEMAT(IDOS)*100.0
         END IF
*  |
*  +-------------------------------------------------------------------*
*  | 1D histogram
      ELSE
         WRITE (93,'(A)') '# Detector n: 1D Projection'
         DO IE=1,NZTABN
              X = DBLE(ZTALOW)+(DBLE(IE)-1.0D+00)*DXBIN
              WRITE (93,'(1P,2E18.10,2E13.5)')
     &               X, X+DXBIN, FLUMAT(IE), FLEMAT(IE)*100.0
              IF(ABS(FLUMAT(IE)).LT.ZMIN.AND.ABS(FLUMAT(IE))
     &           .GT.2.E-38) ZMIN=ABS(FLUMAT(IE))
              IF(ABS(FLUMAT(IE)).GT.ZMAX) ZMAX=ABS(FLUMAT(IE))
         END DO
         WRITE (93,'(1P,2E18.10,2E13.5)')
     &          X+DXBIN, X+DXBIN, FLUMAT(NZTABN), FLEMAT(NZTABN)*100.0
      END IF
*  |
*  +-------------------------------------------------------------------*
      WRITE (*,*) ZMIN, ZMAX
   55 CONTINUE
      CLOSE (UNIT = 93)
      IF (ID .EQ. 1) THEN
* Start_VAX_seq
*        OPEN(UNIT=89,FILE='gplevh.npo',STATUS='NEW')
*        OPEN(UNIT=90,FILE='gplevh.poi',STATUS='NEW')
*        OPEN(UNIT=91,FILE='gplevh.lim',STATUS='NEW')
* End_VAX_seq
* Start_UNIX_seq
         OPEN(UNIT=89,FILE='gplevh.npo',STATUS='UNKNOWN')
         OPEN(UNIT=90,FILE='gplevh.poi',STATUS='UNKNOWN')
         OPEN(UNIT=91,FILE='gplevh.lim',STATUS='UNKNOWN')
* End_UNIX_seq
         IF ( .NOT. LGEOPL ) THEN
            IF ( LAXSWP ) THEN
               YUMIN = ZTALOW
               XUMIN = RADLOW
               YUMAX = ZTAHGH
               XUMAX = RADHGH
            ELSE
               XUMIN = ZTALOW
               YUMIN = RADLOW
               XUMAX = ZTAHGH
               YUMAX = RADHGH
            END IF
         END IF
         IF ( LGAXSW ) THEN
            WRITE(91,*)YUMIN,XUMIN
            WRITE(91,*)YUMAX,XUMAX
         ELSE
            WRITE(91,*)XUMIN,YUMIN
            WRITE(91,*)XUMAX,YUMAX
         END IF
         WRITE(91,*)ZTALOW,RADLOW
         WRITE(91,*)ZTAHGH,RADHGH
         WRITE(91,*)ZMIN,ZMAX
         CLOSE(91)
         LSMPOI=.FALSE.
*  Plot now the geometry:
         IF ( LGEOPL ) THEN
         DO 3000 II=1,100
*  Read the number of worms on tape:
            READ (99,ERR=3600,END=3600) NWORMS, NDUMMY
            IF ( NWORMS.EQ.-1 .AND. NDUMMY.EQ.-1 ) GO TO 3500
            IF ( LCMPOU ) WRITE(98) NWORMS, NDUMMY
            DO 2000 IW = 1, NWORMS
               READ (99) IWORM,IDUMMY,NPOINT
               IF ( IWORM .NE. IW ) STOP 'IWORM-IW'
               READ (99) (XX(J),YY(J),J=1,NPOINT)
               IF (.NOT.LCMPFL)  CALL SHRINK (XX,YY,NPOINT,DXEPS,DYEPS)
               IF ( LCMPOU ) THEN
                  WRITE(98)IWORM,IDUMMY,NPOINT
                  WRITE(98)(XX(J),YY(J),J=1,NPOINT)
               END IF
               DO 1500 IP=1,NPOINT
                  XX (IP) = XX (IP) + R00
                  YY (IP) = YY (IP) + T00
 1500          CONTINUE
               IF ( LROTA ) THEN
                  DO J=1,NPOINT
                     XPOINT (J)= WX0 + ( XX (J) - R00 ) * WXX
     &                               + ( YY (J) - T00 ) * WYX
                     YPOINT (J)= WY0 + ( XX (J) - R00 ) * WXY
     &                               + ( YY (J) - T00 ) * WYY
                     ZPOINT (J)= WZ0 + ( XX (J) - R00 ) * WXZ
     &                               + ( YY (J) - T00 ) * WYZ
                  END DO
                  CALL DOTRSF ( NPOINT, XPOINT, YPOINT, ZPOINT, KROTAT )
                  DO J=1,NPOINT
                     XX (J) = XPOINT (J) * WXX + YPOINT (J) * WXY
     &                      + ZPOINT (J) * WXZ
                     YY (J) = XPOINT (J) * WYX + YPOINT (J) * WYY
     &                      + ZPOINT (J) * WYZ
                  END DO
               END IF
               IF ( XUMIN .LT. 0.E+00 ) THEN
                  XZMIN = 1.000002E+00 * XUMIN
               ELSE
                  XZMIN = 0.999998E+00 * XUMIN
               END IF
               IF ( YUMIN .LT. 0.E+00 ) THEN
                  YZMIN = 1.000002E+00 * YUMIN
               ELSE
                  YZMIN = 0.999998E+00 * YUMIN
               END IF
               IF ( XUMAX .GT. 0.E+00 ) THEN
                  XZMAX = 1.000002E+00 * XUMAX
               ELSE
                  XZMAX = 0.999998E+00 * XUMAX
               END IF
               IF ( YUMAX .GT. 0.E+00 ) THEN
                  YZMAX = 1.000002E+00 * YUMAX
               ELSE
                  YZMAX = 0.999998E+00 * YUMAX
               END IF
               NPPPP = 0
               XXP   = XX (1)
               YYP   = YY (1)
               LINSID=XXP.GE.XZMIN.AND.XXP.LE.XZMAX.AND.YYP.GE.YZMIN
     &                .AND.YYP.LE.YZMAX
               IF ( LINSID ) THEN
                  NPPPP = NPPPP + 1
                  XXX (NPPPP) = XXP
                  YYY (NPPPP) = YYP
               END IF
               DO 1600 IP=2,NPOINT
                  XXO   = XXP
                  YYO   = YYP
                  LINPRE=LINSID
                  XXP   = XX(IP)
                  YYP   = YY(IP)
                  LINSID=XXP.GE.XZMIN.AND.XXP.LE.XZMAX.AND.YYP.GE.YZMIN
     &                   .AND.YYP.LE.YZMAX
                  IF ( LINPRE .AND. LINSID ) THEN
                     NPPPP = NPPPP + 1
                     XXX (NPPPP) = XXP
                     YYY (NPPPP) = YYP
                  ELSE IF ( .NOT. LINPRE .AND. LINSID ) THEN
                     DDXX = XXP - XXO
                     DDX0 = 0.E+00
                     IF ( XXO .LT. XZMIN ) THEN
                        DDX0 = XZMIN - XXO
                     ELSE IF ( XXO .GT. XZMAX ) THEN
                        DDX0 = XZMAX - XXO
                     END IF
                     DDYY = YYP - YYO
                     DDY0 = 0.E+00
                     IF ( YYO .LT. YZMIN ) THEN
                        DDY0 = YZMIN - YYO
                     ELSE IF ( YYO .GT. YZMAX ) THEN
                        DDY0 = YZMAX - YYO
                     END IF
                     IF ( ABS (DDXX) .GT. 1.E-10 .AND. ABS (DDYY)
     &                    .GT. 1.E-10 ) THEN
                        DDYX = DDYY / DDXX * DDX0
                     ELSE IF ( ABS (DDXX) .GT. 1.E-10 ) THEN
                        DDYX = 1.E+10
                     ELSE
                        DDYX = 0.E+00
                     END IF
                     IF ( ABS (DDYX) .GT. ABS (DDY0) ) THEN
                        XXA = XXO + DDX0
                        YYA = YYO + DDYY / DDXX * DDX0
                     ELSE
                        YYA = YYO + DDY0
                        XXA = XXO + DDXX / DDYY * DDY0
                     END IF
                     NPPPP = NPPPP + 1
                     XXX (NPPPP) = XXA
                     YYY (NPPPP) = YYA
                     NPPPP = NPPPP + 1
                     XXX (NPPPP) = XXP
                     YYY (NPPPP) = YYP
                  ELSE IF ( LINPRE .AND. .NOT. LINSID ) THEN
                     DDXX = XXO - XXP
                     DDX1 = 0.E+00
                     IF ( XXP .LT. XZMIN ) THEN
                        DDX1 = XZMIN - XXP
                     ELSE IF ( XXP .GT. XZMAX ) THEN
                        DDX1 = XZMAX - XXP
                     END IF
                     DDYY = YYO - YYP
                     DDY1 = 0.E+00
                     IF ( YYP .LT. YZMIN ) THEN
                        DDY1 = YZMIN - YYP
                     ELSE IF ( YYP .GT. YZMAX ) THEN
                        DDY1 = YZMAX - YYP
                     END IF
                     IF ( ABS (DDXX) .GT. 1.E-10 .AND. ABS (DDYY)
     &                    .GT. 1.E-10 ) THEN
                        DDYX = DDYY / DDXX * DDX1
                     ELSE IF ( ABS (DDXX) .GT. 1.E-10 ) THEN
                        DDYX = 1.E+10
                     ELSE
                        DDYX = 0.E+00
                     END IF
                     IF ( ABS (DDYX) .GT. ABS (DDY1) ) THEN
                        XXA = XXP + DDX1
                        YYA = YYP + DDYY / DDXX * DDX1
                     ELSE
                        YYA = YYP + DDY1
                        XXA = XXP + DDXX / DDYY * DDY1
                     END IF
                     NPPPP = NPPPP + 1
                     XXX (NPPPP) = XXA
                     YYY (NPPPP) = YYA
                  ELSE IF ( .NOT. LINPRE .AND. .NOT. LINSID ) THEN
                     DDXX = XXO - XXP
                     DDX1 = 0.E+00
                     IF ( XXP .LT. XZMIN .AND. XXO .LT. XZMIN ) THEN
                     ELSE IF ( XXP .GT. XZMAX .AND. XXO .GT. XZMAX )THEN
                     ELSE IF ( YYP .LT. YZMIN .AND. YYO .LT. YZMIN )THEN
                     ELSE IF ( YYP .GT. YZMAX .AND. YYO .GT. YZMAX )THEN
                     ELSE
                        DDYY = YYO - YYP
                        IF ( ABS (DDXX) .LT. 1.E-07 * ABS (XZMAX-XZMIN)
     &                     ) THEN
                           XXA  = XXO
                           YYA  = MAX ( YZMIN, MIN ( YYP, YYO ) )
                           NPPPP = NPPPP + 1
                           XXX (NPPPP) = XXA
                           YYY (NPPPP) = YYA
                           XXA  = XXP
                           YYA  = MIN ( YZMAX, MAX ( YYP, YYO ) )
                           NPPPP = NPPPP + 1
                           XXX (NPPPP) = XXA
                           YYY (NPPPP) = YYA
                        ELSE IF ( ABS (DDYY) .LT. 1.E-07
     &                     * ABS (YZMAX-YZMIN) ) THEN
                           YYA  = YYO
                           XXA  = MAX ( XZMIN, MIN ( XXP, XXO ) )
                           NPPPP = NPPPP + 1
                           XXX (NPPPP) = XXA
                           YYY (NPPPP) = YYA
                           YYA  = YYP
                           XXA  = MIN ( XZMAX, MAX ( XXP, XXO ) )
                           NPPPP = NPPPP + 1
                           XXX (NPPPP) = XXA
                           YYY (NPPPP) = YYA
                        ELSE
                           IPPP  = 0
                           EMME  = DDYY / DDXX
                           XXA   = MAX ( XZMIN, MIN ( XXP, XXO ) )
                           YYA   = EMME * ( XXA - XXP ) + YYP
                           IF ( YYA .GE. YZMIN .AND. YYA .LE. YZMAX )
     &                        THEN
                              IPPP  = IPPP  + 1
                              NPPPP = NPPPP + 1
                              XXX (NPPPP) = XXA
                              YYY (NPPPP) = YYA
                           END IF
                           XXA   = MIN ( XZMAX, MAX ( XXP, XXO ) )
                           YYA   = EMME * ( XXA - XXP ) + YYP
                           IF ( YYA .GE. YZMIN .AND. YYA .LE. YZMAX )
     &                        THEN
                              IPPP  = IPPP  + 1
                              NPPPP = NPPPP + 1
                              XXX (NPPPP) = XXA
                              YYY (NPPPP) = YYA
                           END IF
                           EMME  = DDXX / DDYY
                           YYA   = MIN ( YZMAX, MAX ( YYP, YYO ) )
                           XXA   = EMME * ( YYA - YYP ) + XXP
                           IF ( XXA .GE. XZMIN .AND. XXA .LE. XZMAX
     &                         .AND. IPPP .LT. 2 ) THEN
                              IPPP  = IPPP  + 1
                              NPPPP = NPPPP + 1
                              XXX (NPPPP) = XXA
                              YYY (NPPPP) = YYA
                           END IF
                           YYA   = MAX ( YZMIN, MIN ( YYP, YYO ) )
                           XXA   = EMME * ( YYA - YYP ) + XXP
                           IF ( XXA .GE. XZMIN .AND. XXA .LE. XZMAX
     &                         .AND. IPPP .LT. 2 ) THEN
                              IPPP  = IPPP  + 1
                              NPPPP = NPPPP + 1
                              XXX (NPPPP) = XXA
                              YYY (NPPPP) = YYA
                           END IF
                        END IF
                     END IF
                  END IF
 1600          CONTINUE
               IF ( NPPPP .GT. 0 ) WRITE(89,*)NPPPP
               LSMPOI = LSMPOI .OR. NPPPP .GT. 0
               IF ( LGAXSW ) THEN
                  DO J=1,NPPPP
                     WRITE(90,*)YYY(J),XXX(J)
                  END DO
               ELSE
                  DO J=1,NPPPP
                     WRITE(90,*)XXX(J),YYY(J)
                  END DO
               END IF
               WRITE(90,*)
 2000       CONTINUE
 3000    CONTINUE
 3500 CONTINUE
      WRITE(*,*) 'Materials and magnetic field information found.'
* Start_UNIX_seq
      OPEN  ( UNIT = 92, FILE = 'plotgeom.his', STATUS='OLD', ERR=3550 )
      CLOSE ( UNIT = 92, STATUS = 'DELETE' )
 3550 CONTINUE
* End_UNIX_seq
* Start_VAX_seq
*     OPEN ( UNIT=93, FILE='plotgeom.dat', STATUS = 'NEW' )
* End_VAX_seq
* Start_UNIX_seq
      OPEN ( UNIT=93, FILE='plotgeom.dat', STATUS ='UNKNOWN' )
* End_UNIX_seq
      READ(99,ERR=3600,END=3600) IMXPFX, IMXPFY
      IF ( LGAXSW ) THEN
         XMIN=YUMIN
         XMAX=YUMAX
         YMIN=XUMIN
         YMAX=XUMAX
      ELSE
         XMIN=XUMIN
         XMAX=XUMAX
         YMIN=YUMIN
         YMAX=YUMAX
      END IF
      DXSTEP = (XMAX-XMIN)/FLOAT(IMXPFX)
      DYSTEP = (YMAX-YMIN)/FLOAT(IMXPFY)
      DO IY=1,IMXPFY
         Y=YMIN+(FLOAT(IY)-0.5)*DYSTEP
         DO IX=1,IMXPFX
            X=XMIN+(FLOAT(IX)-0.5)*DXSTEP
            READ(99) NREG,MAT,LAT,BX,BY,B
            IF ( LGAXSW ) THEN
               XN = Y
               YN = X
            ELSE
               XN = X
               YN = Y
            END IF
            WRITE (93,*) XN,YN,NREG,MAT,LAT,BX,BY,B
         END DO
         WRITE (93,*)
      END DO
      IC = 0
 3600 CONTINUE
      END IF
      IF ( .NOT. LSMPOI ) THEN
         WRITE (89,*) 0
         WRITE (90,*) 0.E+00, 0.E+00
      END IF
      CLOSE (UNIT=99)
      IF ( LCMPOU ) CLOSE (UNIT=98)
      END IF
      IF (FILE3 .EQ. '  ') GO TO 60
*  End plot of the geometry
   50 CONTINUE
   60 CONTINUE
*9999 CONTINUE
      CLOSE(89)
      CLOSE(90)

      DEALLOCATE(FLUMAT)
      DEALLOCATE(FLEMAT)

      END

      SUBROUTINE SHRINK (XX,YY,NPOINT,DXEPS,DYEPS)

      DIMENSION XX(*), YY(*)
      IF ( NPOINT .LE. 3 ) RETURN
      XSCALE = DXEPS
      YSCALE = DYEPS
      DUEPS  = DXEPS / XSCALE
      DVEPS  = DYEPS / YSCALE
      DDEPS  = SQRT ( DUEPS**2 + DYEPS**2 )
      IPOINT = 1
      IPP    = 1
      ULAST  = XX (1) / XSCALE
      VLAST  = YY (1) / YSCALE
      DO 1000 IP = 2, NPOINT - 1
         UU = XX (IP) / XSCALE
         VV = YY (IP) / YSCALE
         DU = ABS ( UU - ULAST )
         DV = ABS ( VV - VLAST )
         IF ( ( DU .GE. DUEPS .OR. DV .GE. DVEPS ) .AND. IPOINT .EQ. 1 )
     &     THEN
            IPOINT = IPOINT + 1
            XX (IPOINT) = XX (IP)
            YY (IPOINT) = YY (IP)
            ULSM1 = ULAST
            VLSM1 = VLAST
            ULAST = UU
            VLAST = VV
            IPP   = IP
            IF ( ABS (ULAST-ULSM1) .LE. 1.0E-02*DUEPS ) THEN
               EMME  = 1.E+31
            ELSE IF ( ABS (VLAST-VLSM1) .LE. 1.0E-02*DVEPS ) THEN
               EMME  = 0.E+00
            ELSE
               EMME  = ( VLAST - VLSM1 ) / ( ULAST - ULSM1 )
               QU    = VLSM1 - EMME * ULSM1
            END IF
* The line is : y = m x + q
         ELSE IF ( DU .GE. DUEPS .OR. DV .GE. DVEPS ) THEN
* x=const line:
            IF ( ABS (EMME) .GT. 1.E+30 ) THEN
               DIST = ABS (UU-ULAST)
            ELSE IF ( ABS (EMME) .LT. 1.E-30 ) THEN
               DIST = ABS (VV-VLAST)
            ELSE
               DIST = ABS ( QU + EMME * UU - VV )
     &              / SQRT ( 1.E+00 + EMME**2 )
            END IF
            IF ( DIST .GE. DDEPS ) THEN
               IF ( IPP .LT. IP - 1 ) THEN
                  IPOINT = IPOINT + 1
                  XX (IPOINT) = XX (IP-1)
                  YY (IPOINT) = YY (IP-1)
                  ULSM1 = ULAST
                  VLSM1 = VLAST
                  ULAST = XX (IPOINT) / XSCALE
                  VLAST = YY (IPOINT) / YSCALE
                  IPP   = IP - 1
               END IF
               IPOINT = IPOINT + 1
               XX (IPOINT) = XX (IP)
               YY (IPOINT) = YY (IP)
               ULSM1 = ULAST
               VLSM1 = VLAST
               ULAST = UU
               VLAST = VV
               IPP   = IP
               IF ( ABS (ULAST-ULSM1) .LE. 1.0E-02*DUEPS ) THEN
                  EMME  = 1.E+31
               ELSE IF ( ABS (VLAST-VLSM1) .LE. 1.0E-02*DVEPS ) THEN
                  EMME  = 0.E+00
               ELSE
                  EMME  = ( VLAST - VLSM1 ) / ( ULAST - ULSM1 )
                  QU    = VLSM1 - EMME * ULSM1
               END IF
            ELSE
               IPP   = IP
               ULAST = UU
               VLAST = VV
               XX (IPOINT) = XX (IP)
               YY (IPOINT) = YY (IP)
* Do not recompute m,q:
*              IF ( ABS (ULAST-ULSM1) .LE. 1.0E-02*DUEPS ) THEN
*                 EMME  = 1.E+31
*              ELSE IF ( ABS (VLAST-VLSM1) .LE. 1.0E-02*DVEPS ) THEN
*                 EMME  = 0.E+00
*              ELSE
*                 EMME  = ( VLAST - VLSM1 ) / ( ULAST - ULSM1 )
*                 QU    = VLSM1 - EMME * ULSM1
*              END IF
            END IF
         END IF
 1000 CONTINUE
      IPOINT = IPOINT + 1
      XX (IPOINT) = XX (NPOINT)
      YY (IPOINT) = YY (NPOINT)
      NPOINT = IPOINT
      RETURN
      END

      SUBROUTINE UBDOS2 (NZTABN,NRADBN,SZTLOW,SZTHGH,LGEOPL,SRDLOW,
     &                   SRDHGH,RHOTH,NBG,LRHOFL,RHOFIX,LRPHI)

      USE GPLMOD

      CHARACTER RUNTIT*80,RUNTIM*32,CHSTAT*10
*
*=== usrbin ==========================================================*
*
*---------------------------------------------------------------------*
*     Module USRBIN:                                                  *
*     A. Ferrari & A. Fasso': user defined binnings                   *
*                                                                     *
*     Up to MXUSBN user defined binnings are allowed                  *
*            lusbin = logical flag, .true. if at least 1 user defined *
*                     binning is used                                 *
*            lusevt = logical flag, .true. if at least 1 user defined *
*                     binning event by event is used                  *
*            lustkb = logical flag, .true. if at least 1 user defined *
*                     track-length binning is used. Track length bin- *
*                     ning are recognized as binnings where accurate  *
*                     deposition along the track is requested but     *
*                     for generalized particles other than 208/211    *
*            nusrbn = number of user defined binnings used            *
*            itusbn = type of binning: 0 = cartesian, .ne. 0 = RZ     *
*            idusbn = distribution to be scored: the usual values are *
*                     allowed.                                        *
*            titusb = binning name                                    *
*            ipusbn = logical unit to print the results on: formatted *
*                     if > 0, unformatted if < 0                      *
*            kbusbn = initial location in blank common of the consi-  *
*                     dered binning (real*8 address)                  *
*            nxbin  = number of x (r for RZ) intervals                *
*            nybin  = number of y (1 for RZ) intervals                *
*            nzbin  = number of z intervals                           *
*         xlow/high = minimum and maximum x (r for RZ)                *
*         ylow/high = minimum and maximum y (meaningless for RZ)      *
*         zlow/high = minimum and maximum z                           *
*            dxusbn = x (r) bin width                                 *
*            dyusbn = y bin width                                     *
*            dzusbn = z bin width                                     *
*            tcusbn = time cut-off (seconds) for this binning         *
*            bkusbn = 1st Birk's law parameter for this binning       *
*                     (meaningful only for energy scoring)            *
*            b2usbn = 2nd Birk's law parameter for this binning       *
*                     (meaningful only for energy scoring)            *
*            levtbn = logical flag for binning to be printed at the   *
*                     end of each event only                          *
*            lntzer = logical flag for printing only non zero cells   *
*            ltrkbn = logical flag for flagging track-length binnings *
*                                                                     *
*---------------------------------------------------------------------*
*
      PARAMETER ( MXUSBN = 400 )
      LOGICAL LUSBIN, LEVTBN, LNTZER, LUSEVT, LUSTKB, LTRKBN
      CHARACTER*10 TITUSB
      COMMON /UUSRBN/ XLOW  (MXUSBN), XHIGH (MXUSBN), YLOW  (MXUSBN),
     &                YHIGH (MXUSBN), ZLOW  (MXUSBN), ZHIGH (MXUSBN),
     &                DXUSBN(MXUSBN), DYUSBN(MXUSBN), DZUSBN(MXUSBN),
     &                TCUSBN(MXUSBN), BKUSBN(MXUSBN), B2USBN(MXUSBN),
     &                NXBIN (MXUSBN), NYBIN (MXUSBN), NZBIN (MXUSBN),
     &                ITUSBN(MXUSBN), IDUSBN(MXUSBN), KBUSBN(MXUSBN),
     &                IPUSBN(MXUSBN), LEVTBN(MXUSBN), LNTZER(MXUSBN),
     &                LTRKBN(MXUSBN), NUSRBN, LUSBIN, LUSEVT, LUSTKB
      COMMON /USRCH/  TITUSB(MXUSBN)
      PARAMETER ( PIPIPI = 3.141592653589793D+00 )
*     PARAMETER ( MEMORY = 40000000 )
*     DIMENSION GMSTOR (MEMORY), GESTOR (MEMORY), RHO (MEMORY),
      DIMENSION JB (MXUSBN), ACCDIS (MXUSBN)
      REAL,DIMENSION(:),ALLOCATABLE :: GMSTOR, GESTOR, RHO
      DOUBLE PRECISION  ACCDIS, ACCUDS, DOSMAX, DOSMIN, DOSINT,
     &  DOSECO, VOLUCO, ZTAZTA, ZTALOW, ZTAHGH,
     &  DZTABN, RADLOW, RADHGH, DRADBN, XXXMIN, XXXMAX,
     &  YYYMIN, YYYMAX, ZZZMIN, ZZZMAX, XCOLOW, XCOHGH, DXCOBN,
     &  YCOLOW, YCOHGH, DYCOBN, ZCOLOW, ZCOHGH, DZCOBN,
     &  WEIGXT, WEIGX2, WEIGYT, WEIGY2, WEIGZT, WEIGZ2, XCSXCS, YUIYUI,
     &  PARVOL, VOLCHK, VOLPAR, VOLPA0, PARVHL, DOSERR
      LOGICAL LCYL, LFIRST, LRHOFL, LGEOPL, LRADBN, LSTAT, LRPHI
      CHARACTER CARD*80
      PARAMETER ( GEVGRA = 1.60219E-07 )
* For the TP: 1E34 adopted as peak luminosity
      PARAMETER ( PLUMIN = 1.0    E+00 )
*     PARAMETER ( PLUMIN = 1.0    E+34 )
**    PARAMETER ( PLUMIN = 1.7    E+34 )
      PARAMETER ( ALUMIN = 1.0    E+00 )
*     PARAMETER ( ALUMIN = 1.0    E+34 )
      PARAMETER ( SIGIN  = 1.0    E+00 )
*     PARAMETER ( SIGIN  = 8.0    E-26 )
      PARAMETER ( TIMEYR = 1.0    E+00 )
*     PARAMETER ( TIMEYR = 1.0    E+07 )
      PARAMETER ( TIMKHZ = 1.0    E+00 )
*     PARAMETER ( TIMKHZ = 1.0    E-03 )
*
*     DO I = 1, MEMORY
*        GMSTOR (I) = 0.D+00
*        GESTOR (I) = 0.D+00
*        RHO    (I) = 0.E+00
*     END DO
*
      IF ( .NOT. LRHOFL ) THEN
         IOPT = 1
         ULUMIN = PLUMIN
         TIMUSD = TIMKHZ
*        ULUMIN = ALUMIN
*        TIMUSD = TIMEYR
      ELSE
         ULUMIN = ALUMIN
         TIMUSD = TIMEYR
      END IF
      SITILU = SIGIN*ULUMIN*TIMUSD
      LFIRST = .TRUE.
      IB0    = 0
      LSTAT  = .FALSE.
      IRECRD = 0
      READ (11,ERR=120) RUNTIT,RUNTIM,WEIPRI,NCASE,MCASE,NBATCH
      IRECRD = IRECRD + 1
      WRITE(*,*)RUNTIT
      WRITE(*,*)RUNTIM
      WRITE(*,*)WEIPRI
      WRITE(*,*)MCASE, NCASE
      WRITE(*,*)NBATCH
      GO TO 150

100   CONTINUE
      IRECRD = 0
      REWIND 11
      READ (11,ERR=100) RUNTIT,RUNTIM,WEIPRI,NCASE,NBATCH
      IRECRD = IRECRD + 1
      MCASE  = 0
      WRITE(*,*)RUNTIT
      WRITE(*,*)RUNTIM
      WRITE(*,*)WEIPRI
      WRITE(*,*)NCASE
      WRITE(*,*)NBATCH
      GO TO 150

120   CONTINUE
      IRECRD = 0
      REWIND 11
      READ (11,ERR=100) RUNTIT,RUNTIM,WEIPRI,NCASE
      IRECRD = IRECRD + 1
      NBATCH = 1
      MCASE  = 0
      WRITE(*,*)RUNTIT
      WRITE(*,*)RUNTIM
      WRITE(*,*)WEIPRI
      WRITE(*,*)NCASE

150   CONTINUE
      ACCUTT = 0.D0
      ACCUDS = 0.D0
*     KLAST  = 0
      LCYL   = .FALSE.
      NBLAST = 0
*  +-------------------------------------------------------------------*
*  | Read Binnings
      DO 1000 IB=IB0+1,1000
         NB=IB
         READ (11,ERR=400,END=2000) MB, TITUSB(NB), ITUSBN(NB),
     &                   IDUSBN(NB),
     &                   XLOW(NB),XHIGH(NB),NXBIN(NB),
     &                   DXUSBN(NB),YLOW(NB),
     &                   YHIGH(NB),NYBIN(NB),DYUSBN(NB),
     &                   ZLOW(NB),ZHIGH(NB),NZBIN(NB),
     &                   DZUSBN(NB),LNTZER(NB),BKUSBN(NB),
     &                   B2USBN(NB),TCUSBN(NB)
         IRECRD = IRECRD + 1
         GO TO 500
 400     CONTINUE
*        BACKSPACE 11
         REWIND 11
         DO IR = 1, IRECRD
            READ (11)
         END DO
         READ (11,ERR=450,END=2000) MB, TITUSB(NB), ITUSBN(NB),
     &                     IDUSBN(NB),
     &                     XLOW(NB),XHIGH(NB),NXBIN(NB),
     &                     DXUSBN(NB),YLOW(NB),
     &                     YHIGH(NB),NYBIN(NB),DYUSBN(NB),
     &                     ZLOW(NB),ZHIGH(NB),NZBIN(NB),
     &                     DZUSBN(NB)
         IRECRD = IRECRD + 1
         LNTZER(NB)=.FALSE.
         TCUSBN(NB)=1.E+38
         BKUSBN(NB)=0.E+00
         B2USBN(NB)=0.E+00
         GO TO 500
 450     CONTINUE
*        BACKSPACE 11
         REWIND 11
         DO IR = 1, IRECRD
            READ (11)
         END DO
         READ (11,ERR=460) CHSTAT,IJX
         IRECRD = IRECRD + 1
         IF ( CHSTAT .EQ. 'STATISTICS' ) THEN
            LSTAT = .TRUE.
            GO TO 2000
         END IF
 460     CONTINUE
         STOP 'UNKNOWN RECORD TYPE'
 500     CONTINUE
         NBLAST = NB
         WRITE(*,*)' NB',NB,' ',TITUSB(NB)
*  |  +----------------------------------------------------------------*
*  |  | Requested bin ask user for limits
         IF ( NB .EQ. NBG ) THEN
            IF ( MOD (ITUSBN(NB),10) .EQ. 1 .OR.
     &           MOD (ITUSBN(NB),10) .EQ. 7 ) THEN
               LRADBN = .TRUE.
            ELSE
               LRADBN = .FALSE.
            ENDIF
            IRAXIS = 0
            IZAXIS = 0
            WRITE (*,*) XLOW (NB), XHIGH(NB), NXBIN (NB)
            IF ( LRADBN ) THEN
               WRITE (*,*)' Rmin, Rmax (def. provided)'
            ELSE
               WRITE (*,*)' Xmin, Xmax (def. provided)'
            ENDIF
            XCOLOW = XLOW (NB)
            XCOHGH = XHIGH(NB)
            READ (*,'(A)') CARD
            IF ( CARD .NE. ' ' ) READ (CARD,*) XCOLOW, XCOHGH
            IF ( LRADBN ) THEN
               WRITE (*,*)' NR (def. N_R_max)'
            ELSE
               WRITE (*,*)' Nx (def. N_x_max)'
            END IF
            NXCOBN = NXBIN (NB)
            READ (*,'(A)') CARD
            IF ( CARD .NE. ' ' ) READ (CARD,*) NXCOBN
            DXCOBN = ( XCOHGH - XCOLOW ) / DBLE (NXCOBN)
            IF ( NXCOBN .GT. 1 ) THEN
               IRAXIS = 1
               RADLOW = XCOLOW
               RADHGH = XCOHGH
               NRADBN = NXCOBN
               DRADBN = DXCOBN
            END IF
            WRITE (*,*) YLOW (NB), YHIGH(NB), NYBIN (NB)
            IF ( LRADBN ) THEN
               WRITE (*,*)' PHImin, PHImax (def. provided)'
            ELSE
               WRITE (*,*)' Ymin, Ymax (def. provided)'
            ENDIF
            YCOLOW = YLOW (NB)
            YCOHGH = YHIGH(NB)
            READ (*,'(A)') CARD
            IF ( CARD .NE. ' ' ) READ (CARD,*) YCOLOW, YCOHGH
            IF ( LRADBN ) THEN
               WRITE (*,*)' NPHI (def. N_PHI_max)'
            ELSE
               WRITE (*,*)' Ny (def. N_y_max)'
            ENDIF
            NYCOBN = NYBIN (NB)
            READ (*,'(A)') CARD
            IF ( CARD .NE. ' ' ) READ (CARD,*) NYCOBN
            DYCOBN = ( YCOHGH - YCOLOW ) / DBLE (NYCOBN)
            IF ( NYCOBN .GT. 1 .AND. IRAXIS .LE. 0 ) THEN
               IRAXIS = 2
               RADLOW = YCOLOW
               RADHGH = YCOHGH
               NRADBN = NYCOBN
               DRADBN = DYCOBN
            ELSE IF ( NYCOBN .GT. 1 .AND. IRAXIS .GT. 0 ) THEN
               IZAXIS = 2
               ZTALOW = YCOLOW
               ZTAHGH = YCOHGH
               NZTABN = NYCOBN
               DZTABN = DYCOBN
            END IF
            WRITE (*,*) ZLOW (NB), ZHIGH(NB), NZBIN (NB)
            WRITE (*,*)' Zmin, Zmax (def. provided)'
            ZCOLOW = ZLOW (NB)
            ZCOHGH = ZHIGH(NB)
            READ (*,'(A)') CARD
            IF ( CARD .NE. ' ' ) READ (CARD,*) ZCOLOW, ZCOHGH
            WRITE (*,*)' Nz (def. N_z_max)'
            NZCOBN = NZBIN (NB)
            READ (*,'(A)') CARD
            IF ( CARD .NE. ' ' ) READ (CARD,*) NZCOBN
            DZCOBN = ( ZCOHGH - ZCOLOW ) / DBLE (NZCOBN)
            LRPHI = LRADBN .AND. NZCOBN .EQ. 1
            IF ( NZCOBN .GT. 1 .AND. IZAXIS .LE. 0 ) THEN
               IZAXIS = 3
               ZTALOW = ZCOLOW
               ZTAHGH = ZCOHGH
               NZTABN = NZCOBN
               DZTABN = DZCOBN
            ELSE IF ( NZCOBN .LE. 1 .AND. IZAXIS .GT. 0 ) THEN
            END IF
* One dimensional scans
*  Z axis:
            IF (NXCOBN .EQ. 1 .AND. NYCOBN .EQ. 1) THEN
               IRAXIS = 0
               IZAXIS = 3
               ZTALOW = ZCOLOW
               ZTAHGH = ZCOHGH
               NZTABN = NZCOBN
               DZTABN = DZCOBN
               NRADBN = 1
*  Y axis:
            ELSE IF (NXCOBN .EQ. 1 .AND. NZCOBN .EQ. 1) THEN
               IRAXIS = 0
               IZAXIS = 2
               ZTALOW = YCOLOW
               ZTAHGH = YCOHGH
               NZTABN = NYCOBN
               DZTABN = DYCOBN
               NRADBN = 1
*  X axis:
            ELSE IF (NYCOBN .EQ. 1 .AND. NZCOBN .EQ. 1) THEN
               IRAXIS = 0
               IZAXIS = 1
               ZTALOW = XCOLOW
               ZTAHGH = XCOHGH
               NZTABN = NXCOBN
               DZTABN = DXCOBN
               NRADBN = 1
            ELSE IF (NXCOBN .EQ. 1 .AND. NYCOBN .EQ. 1.AND.
     &               NZCOBN .EQ. 1) THEN
               STOP ' ABSURD BINNING'
            END IF
            WRITE (*,*)' IRAXIS,IZAXIS',IRAXIS,IZAXIS
            WRITE (*,*)' ZTALOW,ZTAHGH,DZTABN',ZTALOW,ZTAHGH,DZTABN
            WRITE (*,*)' RADLOW,RADHGH,DRADBN',RADLOW,RADHGH,DRADBN

            ALLOCATE(VLRZTA(NRADBN,NZTABN))
            ALLOCATE(GMRZTA(NRADBN,NZTABN))
            ALLOCATE(GERZTA(NRADBN,NZTABN))

            VLRZTA = 0.D+00
            GMRZTA = 0.D+00
            GERZTA = 0.D+00

         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         JB(IB)     = NB
         ACCDIS(NB) = 0.D+00
*         K0         = KLAST + 1
         IRECRD = IRECRD + 1
         IF ( NB .NE. NBG ) THEN
             KBUSBN(NB) = 0
             READ (11)        ! skip data
         ELSE
             K0 = 1
             K1 = NXBIN(NB) * NYBIN(NB) * NZBIN(NB) + K0 -1
             ALLOCATE(GMSTOR(K1))
             ALLOCATE(GESTOR(K1))
*             KLAST      = K1
             KBUSBN(NB) = K0
*  |  Read data
             READ (11) (GMSTOR(J), J = K0, K1)
         END IF
*  |  +----------------------------------------------------------------*
*  |  |
         IF ( NB .EQ. NBG .AND. LFIRST ) THEN
            LFIRST =.FALSE.
            KK0 = 1
            KK1 = K1 - K0 + 1
            ALLOCATE(RHO(KK1))
            IF ( LRHOFL .AND. RHOFIX .LT. 0.E+00 ) THEN
               READ(12)
               READ(12)
               READ(12) (RHO(J), J = KK0, KK1)
               CLOSE (UNIT=12)
            ELSE
               DO J = KK0,KK1
                  RHO (J) = RHOFIX
               END DO
            END IF
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
 1000 CONTINUE
*  |
*  +-------------------------------------------------------------------*
 2000 CONTINUE
*  +-------------------------------------------------------------------*
*  | Read statistics
      IF ( LSTAT ) THEN
*         KLAST  = 0
*  |  +----------------------------------------------------------------*
*  |  | Read errors
         DO IB = 1, NBLAST
            NB = IB
            IRECRD = IRECRD + 1
            IF ( NB .NE. NBG ) THEN
                READ (11)
            ELSE
                K0 = 1
                K1 = NXBIN(NB) * NYBIN(NB) * NZBIN(NB) + K0 -1
*               KLAST      = K1
                READ (11) (GESTOR(J), J = K0, K1)
            END IF
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
      END IF
*  |
*  +-------------------------------------------------------------------*
      NB = NBG
      IF ( LRADBN ) THEN
         PARVO0 = 0.5D+00 * DYUSBN(NB) * DZUSBN(NB) * DXUSBN(NB)**2
         PARVHL = 2.D+00 * XLOW(NB) / DXUSBN (NB)
      ELSE
         PARVOL = DXUSBN(NB) * DYUSBN(NB) * DZUSBN(NB)
      END IF
*  +-------------------------------------------------------------------*
*  |
      DO 7540 IX = 1,NXBIN(NB)
*  |  +----------------------------------------------------------------*
*  |  |
         IF ( LRADBN ) THEN
            PARVOL = PARVO0 * ( DBLE ( 2*IX - 1 ) + PARVHL )
            IF ( PARVOL .LE. 0.0001 ) THEN
               WRITE(*,*)'PARVOL : PARVO0,PARVHL ,ix '
               WRITE(*,*)PARVOL,PARVO0,PARVHL ,ix
            END IF
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
*  |  +----------------------------------------------------------------*
*  |  |
         DO 7530 IY = 0,NYBIN(NB)-1
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            DO 7520 IZ = 0,NZBIN(NB)-1
               IXYZ0 = ( IZ * NYBIN (NB) + IY ) * NXBIN (NB)
     &                  + IX
               IXYZ  = IXYZ0 + KBUSBN (NB) - 1
               XXX   = XLOW(NB)+(FLOAT(IX)-0.5E+00)*DXUSBN(NB)
               YYY   = YLOW(NB)+(FLOAT(IY)+0.5E+00)*DYUSBN(NB)
               ZZZ   = ZLOW(NB)+(FLOAT(IZ)+0.5E+00)*DZUSBN(NB)
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |
               IF ( ABS (GMSTOR (IXYZ)) .GT. 0.E+00 ) THEN
                  IF ( RHO (IXYZ0) .LE. 1.E-05 .AND. LRHOFL ) THEN
                     WRITE (77,*)' RHO,IXYZ0,X,Y,Z,GMSTOR,NB',
     &                  RHO(IXYZ0),IXYZ0,XXX,YYY,ZZZ,GMSTOR(IXYZ),NB
                     ACCDIS(NB) = ACCDIS(NB) + GMSTOR (IXYZ) * PARVOL
                     GMSTOR (IXYZ) = 0.E+00
                     GESTOR (IXYZ) = 0.E+00
                     DOSECO = 0.D+00
                     DOSERR = 0.D+00
                     VOLUCO = 0.D+00
                  ELSE IF ( RHO (IXYZ0) .LE. RHOTH ) THEN
                     ACCDIS(NB) = ACCDIS(NB) + GMSTOR (IXYZ) * PARVOL
                     GMSTOR (IXYZ) = 0.E+00
                     DOSECO = 0.D+00
                     DOSERR = 0.D+00
                     VOLUCO = 0.D+00
                  ELSE
                     GMSTOR (IXYZ) = GMSTOR (IXYZ) / RHO (IXYZ0)
     &                             * SITILU
                     IF ( LRHOFL ) THEN
                        GMSTOR (IXYZ) = GMSTOR (IXYZ) * GEVGRA
                     END IF
                     DOSECO = DBLE(GMSTOR(IXYZ)) * PARVOL
     &                         * RHO (IXYZ0)
                     DOSERR =  (DBLE(GESTOR(IXYZ)) * DOSECO)**2
                     IF ( ABS (DOSECO) .LT. 1.D-70 ) THEN
     &                     DOSECO = SIGN ( 1.D-70, DOSECO )
                     VOLUCO = DBLE (PARVOL*RHO(IXYZ0))
                  END IF
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |
               ELSE
                  IF ( RHO (IXYZ0) .GT. 1.E-05 ) THEN
                     GMSTOR (IXYZ) = 0.E+00
                     GESTOR (IXYZ) = 0.E+00
                     DOSECO = 1.D-70
                     DOSERR = 0.D+00
                     VOLUCO = DBLE (PARVOL*RHO(IXYZ0))
                  ELSE
                     GMSTOR (IXYZ) = 0.E+00
                     GESTOR (IXYZ) = 0.E+00
                     DOSECO = 0.D+00
                     DOSERR = 0.D+00
                     VOLUCO = 0.D+00
                  END IF
               END IF
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |
               IF ( ABS (DOSECO) .GT. 0.D+00 ) THEN
                  IF ( ABS (DOSECO) .LT. 1.D-60 ) DOSECO = 0.D+00
                  RHOVOL = 0.D+00
                  VOLCHK = 0.D+00
*                 LPOINT=1
                  LPOINT=3
*                 LPOINT=10
*                 MPOINT=1
                  MPOINT=3
*                 MPOINT=10
*                 NPOINT=1
                  NPOINT=3
*                 NPOINT=10
                  WEIGXT=1.D+00/DBLE(LPOINT)
                  WEIGX2=0.5D+00*WEIGXT
                  WEIGYT=1.D+00/DBLE(MPOINT)
                  WEIGY2=0.5D+00*WEIGYT
                  WEIGZT=1.D+00/DBLE(NPOINT)
                  WEIGZ2=0.5D+00*WEIGZT
                  IF ( LRADBN ) THEN
                     VOLPA0 = WEIGZT * DZUSBN(NB) * WEIGY2
     &                      * DYUSBN (NB)
                  ELSE
                     VOLPA0 = WEIGZT * WEIGXT * WEIGYT * PARVOL
                  END IF
*  |  |  |  |  +-------------------------------------------------------*
*  |  |  |  |  |
                  DO 7510 JX=1,LPOINT
                     XXXMIN = XLOW(NB)+(DBLE(IX-1)+WEIGXT*DBLE(JX-1))
     &                      * DXUSBN(NB)
                     XXXMAX = XLOW(NB)+(DBLE(IX-1)+WEIGXT*DBLE(JX))
     &                      * DXUSBN(NB)
                     IF ( LRADBN ) THEN
                        VOLPAR = VOLPA0 * ( XXXMAX - XXXMIN )
     &                         * ( XXXMAX + XXXMIN )
                     ELSE
                        VOLPAR = VOLPA0
                     END IF

*  |  |  |  |  |  +----------------------------------------------------*
*  |  |  |  |  |  |
                     DO 7508 JY=1,MPOINT
                        YYYMIN = YLOW(NB)+(DBLE(IY)
     &                         + WEIGYT*DBLE(JY-1)) * DYUSBN(NB)
                        YYYMAX = YLOW(NB)+(DBLE(IY)
     &                         + WEIGYT*DBLE(JY)) * DYUSBN(NB)
*  |  |  |  |  |  |  +-------------------------------------------------*
*  |  |  |  |  |  |  |
                        DO 7505 JZ=1,NPOINT
                           ZZZMIN = ZLOW(NB)+(DBLE(IZ)
     &                            + WEIGZT*DBLE(JZ-1)) * DZUSBN(NB)
                           ZZZMAX = ZLOW(NB)+(DBLE(IZ)
     &                            + WEIGZT*DBLE(JZ)) * DZUSBN(NB)
                           XCSXCS = 0.5D+00 * ( XXXMIN + XXXMAX )
                           YUIYUI = 0.5D+00 * ( YYYMIN + YYYMAX )
                           ZTAZTA = 0.5D+00 * ( ZZZMIN + ZZZMAX )
                           VOLCHK = VOLCHK + VOLPAR
                           IF ( XCSXCS .LT. XCOLOW .OR. XCSXCS .GE.
     &                          XCOHGH .OR. YUIYUI .LT. YCOLOW .OR.
     &                          YUIYUI .GE. YCOHGH .OR. ZTAZTA .LT.
     &                          ZCOLOW .OR. ZTAZTA .GE. ZCOHGH )
     &                          GO TO 7505
                           IF ( IRAXIS .EQ. 0 ) THEN
                              IRAD = 1
                           ELSE IF ( IRAXIS .EQ. 1 ) THEN
                              IRAD = INT ( ( XCSXCS - RADLOW )
     &                               / DRADBN ) + 1
                           ELSE IF ( IRAXIS .EQ. 2 ) THEN
                              IRAD = INT ( ( YUIYUI - RADLOW )
     &                               / DRADBN ) + 1
                           ELSE
                              STOP ' IRAXIS??'
                           END IF
                           IF ( IZAXIS .EQ. 1 ) THEN
                              IZTA = INT ( ( XCSXCS - ZTALOW )
     &                             / DZTABN ) + 1
                           ELSE IF ( IZAXIS .EQ. 2 ) THEN
                              IZTA = INT ( ( YUIYUI - ZTALOW )
     &                             / DZTABN ) + 1
                           ELSE IF ( IZAXIS .EQ. 3 ) THEN
                              IZTA = INT ( ( ZTAZTA - ZTALOW )
     &                             / DZTABN ) + 1
                           ELSE
                              STOP ' IZAXIS??'
                           END IF
                           IF ( IRAD .LE. 0 ) GO TO 7505
                           IF ( IZTA .LE. 0 ) GO TO 7505
                           IRAD = MIN ( IRAD, NRADBN )
                           IZTA = MIN ( IZTA, NZTABN )
                           GMRZTA (IRAD,IZTA) = GMRZTA (IRAD,IZTA)
     &                            + DOSECO * VOLPAR / PARVOL
                           GERZTA (IRAD,IZTA) = GERZTA (IRAD,IZTA)
     &                            + DOSERR * VOLPAR / PARVOL
                           VLRZTA (IRAD,IZTA) = VLRZTA (IRAD,IZTA)
     &                            + VOLUCO * VOLPAR / PARVOL
7505                    CONTINUE
*  |  |  |  |  |  |  |
*  |  |  |  |  |  |  +-------------------------------------------------*
7508                 CONTINUE
*  |  |  |  |  |  |
*  |  |  |  |  |  +----------------------------------------------------*
7510              CONTINUE
*  |  |  |  |  |
*  |  |  |  |  +-------------------------------------------------------*
                  IF (ABS(PARVOL-VOLCHK) .GT. 1.D-05*PARVOL)
     &               WRITE(*,*)'  IX,IY,IZ,NB,PARVOL,VOLCHK',
     &                            IX,IY,IZ,NB,PARVOL,VOLCHK
               END IF
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
7520        CONTINUE
*  |  |  |
*  |  |  +-------------------------------------------------------------*
7530     CONTINUE
*  |  |
*  |  +----------------------------------------------------------------*
7540  CONTINUE
*  |
*  +-------------------------------------------------------------------*
      DOSMAX = -1.D+70
      DOSMIN = +1.D+70
      DOSINT =  0.D+00

      ALLOCATE(FLUMAT(NRADBN*NZTABN))
      ALLOCATE(FLEMAT(NRADBN*NZTABN))

      DO IE=1,NZTABN
         DO IA=1,NRADBN
            DOSINT = DOSINT + GMRZTA (IA,IE)
            IF ( ABS (GMRZTA (IA,IE)) .GT. 0.D+00 ) THEN
               GERZTA (IA,IE) = SQRT (GERZTA(IA,IE)) / GMRZTA (IA,IE)
               GMRZTA (IA,IE) = GMRZTA (IA,IE) / VLRZTA (IA,IE)
            END IF
            IDOS   = (IA-1)*NZTABN + IE
*  Introduce here the ABS!!
            FLUMAT (IDOS) = GMRZTA (IA,IE)
            FLEMAT (IDOS) = GERZTA (IA,IE)
            IF ( DBLE (ABS(FLUMAT (IDOS))) .LT. 0.999999D+38 )
     &      DOSMAX = MAX ( DOSMAX, DBLE (ABS(FLUMAT (IDOS))) )
            DOSMIN = MIN ( DOSMIN, DBLE (ABS(FLUMAT (IDOS))) )
         END DO
      END DO
      WRITE (*,*) ' Dosmax, Dosmin: ', DOSMAX, DOSMIN
      WRITE (*,*) ' Binning integral: ', DOSINT
      DO IE=1,NZTABN
         DO IA=1,NRADBN
            IDOS = (IA-1)*NZTABN + IE
            IF ( DBLE (ABS(FLUMAT(IDOS))) .GT. 0.D+00 .AND.
     &              DBLE (ABS(FLUMAT(IDOS))) .LT. 1.D-25*DOSMAX ) THEN
               WRITE (*,*)' *** Zeroed!',FLUMAT(IDOS),SNGL(DOSMAX),
     &                                      IE,IA
               FLUMAT (IDOS) = 0.E+00
            END IF
         END DO
      END DO
* Start_VAX_seq
*     OPEN(UNIT=33,FILE='DOSLEV.DAT',STATUS='NEW')
* End_VAX_seq
* Start_UNIX_seq
      OPEN(UNIT=33,FILE='doslev.dat',STATUS='UNKNOWN')
* End_UNIX_seq
      WRITE(33,8700)((FLUMAT((IA-1)*NZTABN+IE),IA=1,NRADBN),IE=1,NZTABN)
8700  FORMAT (1(5X,1P,5(1X,E11.4)))
      CLOSE(UNIT=33)
      IB0=IB-1
      SRDLOW=RADLOW
      SRDHGH=RADHGH
      SZTLOW=ZTALOW
      SZTHGH=ZTAHGH

      DEALLOCATE(GMSTOR)
      DEALLOCATE(GESTOR)
      DEALLOCATE(RHO)

      DEALLOCATE(VLRZTA)
      DEALLOCATE(GMRZTA)
      DEALLOCATE(GERZTA)

      RETURN
      END
