*
*=== Usysuw ===========================================================*
*
      PROGRAM USYSUW
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 2003-2019:  CERN & INFN                            *
*     Copyright (C) 2020-    :  CERN                                   *
*     All Rights Reserved.                                             *
*                                                                      *
*     USer Yiled scoring SUm and Weighting:                            *
*                                                                      *
*     Created on  14 October 1994  by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*     Modified on  5 July 2021                                         *
*                                                                      *
*----------------------------------------------------------------------*
*
      CHARACTER FILE*232,FILE1*232,RUNTIT*80,RUNTIM*32,CHSTAT*10,
     &          CHCTOT*15
*
*=== usryld ===========================================================*
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 2003-2019:  CERN & INFN                            *
*     All Rights Reserved.                                             *
*                                                                      *
*     Module usryld:                                                   *
*                                                                      *
*     A. Ferrari : user defined particle yield scoring                 *
*                                                                      *
*     Up to MXUSYL user defined yield estimators are allowed           *
*            lusryl = logical flag, .true. if at least 1 user defined  *
*                     yield estimator is used                          *
*            lscuyl = logical flag, .true. if at least 1 user defined  *
*                     yield estimator for secondaries coming out from  *
*                     hadronic inelastic interactions (nr1=-1, nr2=-2) *
*                     is used                                          *
*            lusydf = logical flag, .true. if the momentum etc quanti- *
*                     ties required for cms scoring have been explici- *
*                     tly defined by the user                          *
*            nusryl = number of user defined yeld estimators           *
*            itusyl = type of binning: itusyl = ie + ia*100            *
*                     |ie| = 1 : kinetic  energy binning               *
*                     |ie| = 2 : total momentum  binning               *
*                     |ie| = 3 : lab.   rapidity binning               *
*                     |ie| = 4 : cms    rapidity binning               *
*                     |ie| = 5 : lab. pseudorap. binning               *
*                     |ie| = 6 : cms  pseudorap. binning               *
*                     |ie| = 7 : lab. x          binning               *
*                     |ie| = 8 : cms  Feynmann x binning               *
*                     |ie| = 9 : transverse mom. binning               *
*                     |ie| =10 : transverse mass binning               *
*                     |ie| =11 : lab. long. mom. binning               *
*                     |ie| =12 : cms  long. mom. binning               *
*                     |ie| =13 : total energy    binning               *
*                     |ie| =14 : lab. angle      binning               *
*                     |ie| =15 : cms  angle      binning               *
*                     |ie| =16 : p_t squared     binning               *
*                     |ie| =17 : lab. angle      binning               *
*                                with 1/[2pi sin(theta)] weight        *
*                     |ie| =18 : p_t             binning               *
*                                with 1/(2pi p_t) weight               *
*                     |ie| =19 : frac. lab mom.  binning               *
*                     |ie| =20 : trans. kin. en. binning               *
*                     |ie| =21 : excitation en.  binning               *
*                     |ie| =22 : particle charge binning               *
*                     |ie| =23 : particle LET    binning               *
*                     |ie| =24 : lab. angle(deg) binning               *
*                     |ie| =25 : cms  angle(deg) binning               *
*                     |ie| =26 : lab. kin.en./n  binning               *
*                     |ie| =27 : lab. mom/n      binning               *
*                     |ie| =28 : part.bar.charge binning               *
*                     |ie| =29 : 4 mom.tran. -t  binning               *
*                     |ie| =30 : cms X_F_long    binning               *
*                     |ie| =31 : M_x squared     binning               *
*                     |ie| =32 : M_x squared / s binning               *
*                     |ie| =33 : time (s)        binning               *
*                     |ie| =34 : lab. angle      binning               *
*                                with sin(theta) weight                *
*                     |ie| =35 : cms momentum    binning               *
*                     |ie| =36 : cms tot. energy binning               *
*                     |ie| =37 : lab. angle (deg)binning               *
*                                with 1/[2pi sin(theta)] weight        *
*                     |ie| =38 : lab. angle (deg)binning               *
*                                with sin(theta) weight                *
*                     ie > 0 --> linear, ie < 0 --> logarithmic        *
*                     ia has the same meaning but for the 2nd variable *
*            ixusyl = cross section kind, ixa + ixm * 100              *
*                     ixa = 1 : plain d2 sigma / d x1 d x2 where x1,x2 *
*                               are the first and second variables     *
*                     ixa = 2 : invariant cross section E d3 sigma/dp3 *
*                     ixa = 3 : plain d2 N / d x1 d x2 where x1,x2 are *
*                               the first and second variables         *
*                     ixa = 4 : d2 (x2 N) / d x1 d x2  where x1,x2 are *
*                               the first and second variables         *
*                     ixa = 5 : d2 (x1 N) / d x1 d x2  where x1,x2 are *
*                               the first and second variables         *
*                     ixa = 6 : d2 N / d x1 d x2 cos(theta) where x1,x2*
*                               are the first and second variables,    *
*                               and theta is the angle with the normal *
*                               to the crossed surface                 *
*                     ixa = 7 : d2 (x2 x2 N) / d x1 d x2  where x1,x2  *
*                               are the first and second variables     *
*                     ixa = 8 : d2 (x1 x1 N) / d x1 d x2  where x1,x2  *
*                               are the first and second variables     *
*                     ixa = 9 : d2 N / (x2 d x1 d x2) where x1, x2     *
*                               are the first and second variables     *
*                     ixa =10 : d2 N / (x1 d x1 d x2) where x1, x2     *
*                               are the first and second variables     *
*                     ixa =11 : d2 N / d (x1/x2) d x2) where x1, x2    *
*                               are the first and second variables     *
*                     ixa =16 : d2 (x2 N) / d x1 d x2 cos(theta) where *
*                               x1,x2 are the first and second varia-  *
*                               bles, and theta is the angle with the  *
*                               normal to the crossed surface          *
*                     ixa =26 : d2 (x1 N) / d x1 d x2 cos(theta) where *
*                               x1,x2 are the first and second varia-  *
*                               bles, and theta is the angle with the  *
*                               normal to the crossed surface          *
*                     ixa =51,52: forbidden                            *
*                     ixa >52 : as per ixa = ixa - 50, but d2 sigma ...*
*                               instead of d2 N..., eg ixa=53 is the   *
*                               same as ixa=1, ixa 59 is
*                               d2 sigma / (x2 d x1 d x2) etc          *
*                     ixm = material number for cross section or LET   *
*                           calculations                               *
*            idusyl = distribution to be scored: the usual values are  *
*                     allowed                                          *
*            jdusyl = secondary distribution to be scored: the usual   *
*                     values are allowed                               *
*            nr1uyl = 1st region (from...)                             *
*            nr2uyl = 2nd region (to  ...)                             *
*            usnryl = normalization factor (itusyl dependent)          *
*            sgusyl = normalization cross section (itusyl dependent)   *
*            llnuyl = no low energy neutron scoring if false, yes if   *
*                     true                                             *
*            tituyl = yield estimator name                             *
*            dequyl = yield name for dose equivalent conversions       *
*            ipusyl = logical unit to print the results on: formatted  *
*                     if > 0, unformatted if < 0                       *
*            igmuyl = maximum low energy neutron group to be scored    *
*            jmtuyl = material index for cross section or LET calcula- *
*                     tions                                            *
*            kbusyl = initial location in blank common of the consi-   *
*                     dered yield estimator (real*8 address)           *
*                     I T   I S   F O R  A D D R E S S  OF  1 !!!!!    *
*            neylbn = number of energy or other quantity intervals     *
*     eyllow/eylhgh = minimum and maximum energies (or other)          *
*     ayllow/aylhgh = minimum and maximum angles (or other)            *
*            deylbn = energy (or other) bin width                      *
*            pusryl = lab. momentum of projectile to be used to define *
*                    (possible) Lorentz transformations, Feynmann X    *
*                     etc.                                             *
*            pusryt = total lab. momentum of projectile (different     *
*                     from pusryl for ion projectiles only)            *
*            pcmuyl = CMS momentum of projectile to be used to define  *
*                    (possible) Lorentz transformations, Feynmann X    *
*                     etc.                                             *
*            eusryl = lab. total energy of projectile to be used to    *
*                     define (possible) Lorentz transformations,       *
*                     Feynmann X etc.                                  *
*            eusryt = lab. total energy of projectile (different       *
*                     from eusryl for ion projectiles only)            *
*            ijusyl = projectile identity                              *
*            jtusyl = target     identity                              *
*        u,v,wusryl = laboratory projectile direction                  *
*        u,v,wcmuyl = cms        projectile direction                  *
*    etx,ety,etzuyl = eta for Lorentz transformation                   *
*            gamuyl = gamma for Lorentz transformation                 *
*            sqsuyl = cms energy for Lorentz transformation            *
*     iausyl,ieusyl = ausiliary arrays where itusyl is decoded         *
*            kdcyuy = flag for radioactive decay scoring (-1=prompt +  *
*                     all radioactive decays, 0=prompt only, i=i_th    *
*                     decay time only)                                 *
*                                                                      *
*----------------------------------------------------------------------*
*
      PARAMETER ( MKUSYL =   38 )
      PARAMETER ( MXUSYL = 1000 )
      LOGICAL LUSRYL, LLNUYL, LSCUYL, LUSYDF
      CHARACTER*10 TITUYL, DEQUYL
      COMMON /UUSRYL/ EYLLOW(MXUSYL), EYLHGH(MXUSYL), DEYLBN(MXUSYL),
     &                USNRYL(MXUSYL), SGUSYL(MXUSYL), AYLLOW(MXUSYL),
     &                AYLHGH(MXUSYL), PUSRYL, PUSRYT, PCMUYL, EUSRYL,
     &                EUSRYT, UUSRYL, VUSRYL, WUSRYL, ETXUYL, ETYUYL,
     &                ETZUYL, GAMUYL, SQSUYL, UCMUYL, VCMUYL, WCMUYL,
     &                IJUSYL, JTUSYL,
     &                NEYLBN(MXUSYL), NR1UYL(MXUSYL), NR2UYL(MXUSYL),
     &                IXUSYL(MXUSYL), ITUSYL(MXUSYL), IDUSYL(MXUSYL),
     &                JDUSYL(MXUSYL), KBUSYL(MXUSYL), IPUSYL(MXUSYL),
     &                IGMUYL(MXUSYL), IEUSYL(MXUSYL), IAUSYL(MXUSYL),
     &                JMTUYL(MXUSYL), KDCYUY(MXUSYL), LLNUYL(MXUSYL),
     &                NUSRYL, LUSRYL, LSCUYL, LUSYDF
      COMMON /UUSYCH/ TITUYL(MXUSYL), DEQUYL(MXUSYL)

      COMMON /TCHELP/ DAYLBN(MXUSYL), NAYLBN(MXUSYL)

      PARAMETER ( PIPIPI = 3.141592653589793D+00 )
      PARAMETER ( HLFHLF = 0.5D+00 )
      PARAMETER ( ZERZER = 0.D+00  )
      PARAMETER ( ONEONE = 1.D+00  )
      PARAMETER ( TWOTWO = 2.D+00  )
      PARAMETER ( TWOPIP = 2.D+00 * PIPIPI )

      REAL,DIMENSION(:),ALLOCATABLE :: GMSTOR, GESTOR, ENGMAX, EPGMAX,
     &                                 EPGMCK, ACCUMH, ACCUEH
      DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE :: GDSTOR, GBSTOR,
     &                                             ACCUMG, ACUMCU
      DOUBLE PRECISION TOTTOT, TOTERR, ACCUMP, AFACT , DSCORE, DSCOR2,
     &                 DELTAE, EFACT , WCTOT , GBHELP, TOTHLP
      DIMENSION TOTTOT (MXUSYL), TOTERR (MXUSYL), KFUSYL (MXUSYL),
     &          JX (MXUSYL)    , TOTTHL (MXUSYL), TOTEHL (MXUSYL)
      LOGICAL LSTATI
      CHARACTER CHVARI(MKUSYL)*30, CHYILD(3,76)*49

      DATA CHVARI / ' Laboratory Kinetic Energy    ',
     &              ' Laboratory Total Momentum    ',
     &              ' Laboratory Rapidity          ',
     &              ' CMS rapidity                 ',
     &              ' Laboratory Pseudorapidity    ',
     &              ' CMS Pseudorapidity           ',
     &              ' Laboratory Feynmann X        ',
     &              ' CMS Feynmann X               ',
     &              ' Transverse Momentum          ',
     &              ' Transverse Mass              ',
     &              ' Lab. Longitudinal Momentum   ',
     &              ' CMS Longitudinal Momentum    ',
     &              ' Laboratory Total Energy      ',
     &              ' Laboratory Angle (radians)   ',
     &              ' CMS Angle (radians)          ',
     &              ' Squared transverse momentum  ',
     &              ' Weighted Lab. Angle (radians)',
     &              ' Weighted transverse momentum ',
     &              ' Fractional momentum          ',
     &              ' Transverse kinetic energy    ',
     &              ' Excitation energy            ',
     &              ' Particle charge              ',
     &              ' Particle LET                 ',
     &              ' Laboratory Angle (degrees)   ',
     &              ' CMS Angle (degrees)          ',
     &              ' Lab. Kinetic Energy / nucleon',
     &              ' Laboratory Momentum / nucleon',
     &              ' Particle baryon number       ',
     &              ' four momentum transfer -t    ',
     &              ' CMS Long. Feynmann X         ',
     &              ' Excited mass squared         ',
     &              ' Excited mass squared / s     ',
     &              ' Time                         ',
     &              ' Sin(th) Wei. Lab. Angle (rad)',
     &              ' CMS Total Momentum           ',
     &              ' CMS Total Energy             ',
     &              ' Weighted Lab. Angle (degrees)',
     &              ' Sin(th) Wei. Lab. Angle (deg)'/
      DATA ( ( CHYILD (J,I), J = 1, 3 ), I = 1, 5 )
     &            / ' plain d2 sigma / d x1 d x2 where x1, x2 are     ',
     &              ' the first and second variables                  ',
     &              '                                                 ',
     &              ' invariant cross section E d3 sigma / dp3, diff. ',
     &              ' in x1, x2, the first and second variables)      ',
     &              '                                                 ',
     &              ' plain d2 N / d x1 d x2 where x1, x2 are         ',
     &              ' the first and second variables                  ',
     &              '                                                 ',
     &              ' d2 (x2 N) / d x1 d x2 where x1, x2 are          ',
     &              ' the first and second variables                  ',
     &              '                                                 ',
     &              ' d2 (x1 N) / d x1 d x2 where x1, x2 are          ',
     &              ' the first and second variables                  ',
     &              '                                                 '/
      DATA ( ( CHYILD (J,I), J = 1, 3 ), I = 6, 10 )
     &            / ' d2 N / d x1 d x2 cos(theta) where x1,x2 are     ',
     &              ' the first and second variables, and theta is the',
     &              ' the angle with the normal to the crossed surface',
     &              ' d2 (x2 x2 N) / d x1 d x2 where x1, x2 are       ',
     &              ' the first and second variables                  ',
     &              '                                                 ',
     &              ' d2 (x1 x1 N) / d x1 d x2 where x1, x2 are       ',
     &              ' the first and second variables                  ',
     &              '                                                 ',
     &              ' d2 N / (x2 d x1 d x2) where x1, x2 are          ',
     &              ' the first and second variables                  ',
     &              '                                                 ',
     &              ' d2 N / (x1 d x1 d x2) where x1, x2 are          ',
     &              ' the first and second variables                  ',
     &              '                                                 '/
      DATA ( ( CHYILD (J,I), J = 1, 3 ), I = 11, 15 )
     &            / ' d2 N / d (x1/x2) d x2) where x1, x2 are         ',
     &              ' the first and second variables                  ',
     &              '                                                 ',
     &              ' d2 (x1 Sqrt[x1^2-x2^2] N) / d x1 d x2 where     ',
     &              ' x1, x2 are the first and second variables       ',
     &              '                                                 ',
     &              ' d2 ([x1^2+x2^2] N) / d x1 d x2) where x1,       ',
     &              ' x2 are the first and second variables           ',
     &              '                                                 ',
     &              ' d2 ([x1+x2^2] N) / (pi d x1 d x2) where x1,     ',
     &              ' x2 are the first and second variables           ',
     &              '                                                 ',
     &              ' d2 ([x1^2+x2] N) / (pi d x1 d x2) where x1,     ',
     &              ' x2 are the first and second variables           ',
     &              '                                                 '/
      DATA ( ( CHYILD (J,I), J = 1, 3 ), I = 16, 16 )
     &            / ' d2 (x2 N) / d x1 d x2 cos(theta) where x1, x2 ar',
     &              'e the first and second variables, and theta is th',
     &              'e angle with the normal to the crossed surface   '/
      DATA ( ( CHYILD (J,I), J = 1, 3 ), I = 26, 26 )
     &            / ' d2 (x1 N) / d x1 d x2 cos(theta) where x1, x2 ar',
     &              'e the first and second variables, and theta is th',
     &              'e angle with the normal to the crossed surface   '/
      DATA ( ( CHYILD (J,I), J = 1, 3 ), I = 51, 55 )
     &            / ' plain d2 sigma / d x1 d x2 where x1, x2 are     ',
     &              ' the first and second variables                  ',
     &              '                                                 ',
     &              ' invariant cross section E d3 sigma / dp3, diff. ',
     &              ' in x1, x2, the first and second variables)      ',
     &              '                                                 ',
     &              ' plain d2 sigma / d x1 d x2 where x1, x2 are     ',
     &              ' the first and second variables                  ',
     &              '                                                 ',
     &              ' d2 (x2 sigma) / d x1 d x2 where x1, x2 are      ',
     &              ' the first and second variables                  ',
     &              '                                                 ',
     &              ' d2 (x1 sigma) / d x1 d x2 where x1, x2 are      ',
     &              ' the first and second variables                  ',
     &              '                                                 '/
      DATA ( ( CHYILD (J,I), J = 1, 3 ), I = 56, 60 )
     &            / ' d2 sigma / d x1 d x2 cos(theta) where x1,x2 are ',
     &              ' the first and second variables, and theta is the',
     &              ' the angle with the normal to the crossed surface',
     &              ' d2 (x2 x2 N) / d x1 d x2 where x1, x2 are       ',
     &              ' the first and second variables                  ',
     &              '                                                 ',
     &              ' d2 (x1 x1 sigma) / d x1 d x2 where x1, x2 are   ',
     &              ' the first and second variables                  ',
     &              '                                                 ',
     &              ' d2 sigma / (x2 d x1 d x2) where x1, x2 are      ',
     &              ' the first and second variables                  ',
     &              '                                                 ',
     &              ' d2 sigma / (x1 d x1 d x2) where x1, x2 are      ',
     &              ' the first and second variables                  ',
     &              '                                                 '/
      DATA ( ( CHYILD (J,I), J = 1, 3 ), I = 61, 65 )
     &            / ' d2 sigma / d (x1/x2) d x2) where x1, x2 are     ',
     &              ' the first and second variables                  ',
     &              '                                                 ',
     &              ' d2 (x1 Sqrt[x1^2-x2^2] sigma) / d x1 d x2 where ',
     &              ' x1, x2 are the first and second variables       ',
     &              '                                                 ',
     &              ' d2 ([x1^2+x2^2] sigma) / d x1 d x2) where x1,   ',
     &              ' x2 are the first and second variables           ',
     &              '                                                 ',
     &              ' d2 ([x1+x2^2] sigma) / (pi d x1 d x2) where x1, ',
     &              ' x2 are the first and second variables           ',
     &              '                                                 ',
     &              ' d2 ([x1^2+x2] sigma) / (pi d x1 d x2) where x1, ',
     &              ' x2 are the first and second variables           ',
     &              '                                                 '/
      DATA ( ( CHYILD (J,I), J = 1, 3 ), I = 66, 66 )
     &            / ' d2 (x2 sig) / d x1 d x2 cos(theta) where x1, x2 ',
     &              'are the first and second variables, and theta is ',
     &              'the angle with the normal to the crossed surface '/
      DATA ( ( CHYILD (J,I), J = 1, 3 ), I = 76, 76 )
     &            / ' d2 (x1 sig) / d x1 d x2 cos(theta) where x1, x2 ',
     &              'are the first and second variables, and theta is ',
     &              'the angle with the normal to the crossed surface '/
*
*     In order to avoid computer nuisances
*
      KK0    =-1000000000
      KMAX   =-1000000000
      ENERGE =-1.E+30
      DELTAE =-1.D+30
      EFACT  =-1.D+30
*
*     Note that all addresses are for zero index (KBUSYL,KFUSYL)
*
      DO 4444 I = 1, MXUSYL
         TOTTOT (I) = ZERZER
         TOTERR (I) = ZERZER
 4444 CONTINUE
      NCTOT   = 0
      MCTOT   = 0
      WCTOT   = ZERZER
      MBATCH  = 0
      MAXSIZE = 0
      MNEYLBN = 0
      MIGMUYL = 0
      IFILE   = 0
      LSTATI  = .FALSE.
      LQOLD   = 0
 1111 CONTINUE
      WRITE(*,'('' Type the input file: '',$)')
      READ (*,'(A)')FILE
      LQ     = LNNBLN(FILE)
      IF (LQ .LE. 0) GO TO 2000
      OPEN (UNIT=1,FILE=FILE,STATUS='OLD',FORM='UNFORMATTED',ERR=1111)
      IFILE  = IFILE + 1
      IRECRD = 0
      READ (1,ERR=100) RUNTIT,RUNTIM,WEIPRI,NCASE,MCASE,NBATCH
      IRECRD = IRECRD + 1
      WRITE(*,*)RUNTIT
      WRITE(*,*)RUNTIM
      WRITE(*,*)WEIPRI
      WRITE(*,*)NCASE
      WRITE(*,*)MCASE
      WRITE(*,*)NBATCH
      GO TO 200
 100  CONTINUE
      REWIND 1
      IRECRD = 0
      READ (1,ERR=110) RUNTIT,RUNTIM,WEIPRI,NCASE,NBATCH
      IRECRD = IRECRD + 1
      WRITE(*,*)RUNTIT
      WRITE(*,*)RUNTIM
      WRITE(*,*)WEIPRI
      WRITE(*,*)NCASE
      MCASE  = 0
      GO TO 200
 110  CONTINUE
      REWIND 1
      IRECRD = 0
      READ (1) RUNTIT,RUNTIM,WEIPRI,NCASE
      IRECRD = IRECRD + 1
      WRITE(*,*)RUNTIT
      WRITE(*,*)RUNTIM
      WRITE(*,*)WEIPRI
      WRITE(*,*)NCASE
      MCASE  = 0
      NBATCH = 1
 200  CONTINUE
      READ (1) IJUSYL, JTUSYL, PUSRYL, SQSUYL, UUSRYL, VUSRYL, WUSRYL
      IRECRD = IRECRD + 1
      NCTOT  = NCTOT + NCASE
      MCTOT  = MCTOT + MCASE
      IF ( NCTOT .GE. 1000000000 ) THEN
         NCTOT = NCTOT - 1000000000
         MCTOT = MCTOT + 1
      END IF
      WCTOT  = WCTOT  + DBLE (WEIPRI)
      MBATCH = MBATCH + NBATCH

*  +-------------------------------------------------------------------*
*  | First pass to determine array size needed
      IF (MAXSIZE .EQ. 0) THEN
         IRKEEP = IRECRD
         DO IX = 1, MXUSYL
            NYL = IX
            READ (1,ERR=320,END=390)
     &                MYL, TITUYL(NYL), ITUSYL(NYL),
     &                IXUSYL(NYL), IDUSYL(NYL),
     &                NR1UYL(NYL), NR2UYL(NYL), USNRYL (NYL),
     &                SGUSYL(NYL), LLNUYL(NYL),
     &                EYLLOW(NYL), EYLHGH(NYL),
     &                NEYLBN(NYL), DEYLBN(NYL),
     &                AYLLOW(NYL), AYLHGH(NYL)
            IRECRD = IRECRD + 1
            IF ( LLNUYL (NYL) ) THEN
               READ (1) IGMUYL(NYL)
               IRECRD = IRECRD + 1
            ELSE
               IGMUYL (NYL) = 0
            END IF
            GO TO 350

 320        CONTINUE
*           BACKSPACE 1
            REWIND 1
            DO IR = 1, IRECRD
               READ (1)
            END DO
            READ (1,ERR=330) CHSTAT,IJX
            IRECRD = IRECRD + 1
            IF ( CHSTAT .EQ. 'STATISTICS') GO TO 390

 330        CONTINUE
            STOP 'UNKNOWN RECORD TYPE'

 350        CONTINUE

            NAYLBN(NYL) = 1

            READ (1)
            IRECRD = IRECRD + 1

            MAXSIZE = MAXSIZE + ( 2 + NAYLBN(NYL) ) *
     &                   ( NEYLBN(NYL) + IGMUYL(NYL) )
            MNEYLBN = MAX ( MNEYLBN, NEYLBN(NYL) )
            MIGMUYL = MAX ( MIGMUYL, IGMUYL(NYL) )

         END DO
 390     CONTINUE

         ALLOCATE( GMSTOR(MAXSIZE) )
         ALLOCATE( GDSTOR(MAXSIZE) )
         ALLOCATE( GBSTOR(MAXSIZE) )
         ALLOCATE( GESTOR(MAXSIZE) )

         GBSTOR = ZERZER
         GDSTOR = ZERZER
         GMSTOR = 0.E+00
         GESTOR = 0.E+00

         ALLOCATE( ENGMAX(MIGMUYL + 1) )
         ALLOCATE( EPGMAX(MNEYLBN + MIGMUYL + 1) )
         ALLOCATE( EPGMCK(MNEYLBN + MIGMUYL + 1) )
         ALLOCATE( ACCUMG(MNEYLBN + MIGMUYL + 1) )
         ALLOCATE( ACUMCU(MNEYLBN + MIGMUYL + 1) )
         ALLOCATE( ACCUMH(MNEYLBN + MIGMUYL + 1) )
         ALLOCATE( ACCUEH(MNEYLBN + MIGMUYL + 1) )

         REWIND 1
         DO IR = 1, IRKEEP
             READ (1)
         END DO
         IRECRD = IRKEEP

      END IF
*  |
*  +-------------------------------------------------------------------*

      KLAST  = 0
      KMAX   = 0
*  +-------------------------------------------------------------------*
*  |  No statistics inside (or old fashioned file):
      IF ( NBATCH .LE. 1 ) THEN
*  |  +----------------------------------------------------------------*
*  |  |                       Read from binary file the results!!!
         DO 3600 IX = 1, MXUSYL
            NYL = IX
            READ (1,ERR=3400,END=1000)
     &                MYL, TITUYL(NYL), ITUSYL(NYL),
     &                IXUSYL(NYL), IDUSYL(NYL),
     &                NR1UYL(NYL), NR2UYL(NYL), USNRYL (NYL),
     &                SGUSYL(NYL), LLNUYL(NYL),
     &                EYLLOW(NYL), EYLHGH(NYL),
     &                NEYLBN(NYL), DEYLBN(NYL),
     &                AYLLOW(NYL), AYLHGH(NYL)
            AAUSYL = 0.01D+00 * DBLE (ITUSYL(NYL))
            IAUSYL (NYL) = NINT (AAUSYL)
            IEUSYL (NYL) = ITUSYL (NYL) - IAUSYL (NYL) * 100
            IRECRD = IRECRD + 1
            IF ( LLNUYL (NYL) ) THEN
               READ (1) IGMUYL(NYL), (ENGMAX (IG), IG=1,IGMUYL(NYL)+1)
               IRECRD = IRECRD + 1
            ELSE
               IGMUYL (NYL) = 0
            END IF
            GO TO 3500
 3400       CONTINUE
*           BACKSPACE 1
            REWIND 1
            DO IR = 1, IRECRD
               READ (1)
            END DO
            READ (1,ERR=3460,END=1000) CHSTAT,IJX
            IRECRD = IRECRD + 1
            IF ( CHSTAT .EQ. 'STATISTICS') THEN
               LSTATI = .TRUE.
               GO TO 1000
            END IF
 3460       CONTINUE
            STOP 'UNKNOWN RECORD TYPE'
 3500       CONTINUE
******
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            IF ( IAUSYL (NYL) .EQ. 14 .OR. IAUSYL (NYL) .EQ. 15 .OR.
     &           IAUSYL (NYL) .EQ. 17 .OR. IAUSYL (NYL) .EQ. 34 ) THEN
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |  Small angle (exact!):
               IF ( AYLHGH (NYL) .LT. 0.3D+00 * PIPIPI ) THEN
                  TTLOW  = TAN (HLFHLF*AYLLOW(NYL))
                  TTHGH  = TAN (HLFHLF*AYLHGH(NYL))
                  DAYLBN (NYL) = TWOPIP * TWOTWO * ( TTHGH + TTLOW )
     &                      * ( TTHGH - TTLOW ) / ( ONEONE + TTLOW**2 )
     &                      / ( ONEONE + TTHGH**2 )
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |  Large angle:
               ELSE
                  DAYLBN (NYL) = TWOPIP * ( COS (DBLE(AYLLOW(NYL)))
     &                         - COS (DBLE(AYLHGH(NYL))) )
               END IF
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
*  |  |  |
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            ELSE
               DAYLBN (NYL) = AYLHGH (NYL) - AYLLOW (NYL)
            END IF
*  |  |  |
*  |  |  +-------------------------------------------------------------*
            NAYLBN(NYL) = 1
******
            JX (IX) = NYL
***************************************************************
*********** Note KBUSYL is used for 0 address *****************
***************************************************************
            KBUSYL (NYL) = KLAST
            K0 = KLAST + 1
            K1 = NEYLBN (NYL) + K0 - 1
            K2 = K1 + IGMUYL (NYL)
            KLAST = K2
            READ (1) (GMSTOR(J), J = K0, K2)
            IRECRD = IRECRD + 1
            DO J = K0, K2
*              if (gmstor(j).lt.0.e+00) write(*,*)' *** Nyl,k,gmstor',
*     &           nyl,j-k0+1,gmstor(j)
               GDSTOR (J) = GDSTOR (J) + DBLE (GMSTOR(J))
     &                    * ( DBLE (NCASE) + DBLE (MCASE) * 1.D+09 )
               GBSTOR (J) = GBSTOR (J) + DBLE (GMSTOR(J))**2
     &                    * ( DBLE (NCASE) + DBLE (MCASE) * 1.D+09 )
            END DO
            KMAX = MAX ( KMAX, K2 )
 3600    CONTINUE
 1000    CONTINUE
         IX = IX - 1
*  |  +----------------------------------------------------------------*
*  |  |
         IF ( MBATCH .EQ. 1 .OR. IFILE .EQ. 1 ) THEN
            DO IIX = 1, IX
               NYL = JX (IIX)
               KFUSYL (NYL) = KLAST
               KLAST = KLAST + 2 * ( NEYLBN (NYL) + IGMUYL (NYL) )
               KMAX  = MAX ( KMAX, KLAST )
            END DO
            FILE1 = FILE
            LQOLD = LQ
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         CLOSE (UNIT=1)
*  |
*  +-------------------------------------------------------------------*
*  |  Treat the cases with statistics inside (when NBATCH.GT.1)
      ELSE
*  |  +----------------------------------------------------------------*
*  |  |                       Read from binary file the results!!!
         DO 4400 IX = 1, MXUSYL
            NYL = IX
            READ (1,ERR=4600)
     &                MYL, TITUYL(NYL), ITUSYL(NYL),
     &                IXUSYL(NYL), IDUSYL(NYL),
     &                NR1UYL(NYL), NR2UYL(NYL), USNRYL (NYL),
     &                SGUSYL(NYL), LLNUYL(NYL),
     &                EYLLOW(NYL), EYLHGH(NYL),
     &                NEYLBN(NYL), DEYLBN(NYL),
     &                AYLLOW(NYL), AYLHGH(NYL)
            AAUSYL = 0.01D+00 * DBLE (ITUSYL(NYL))
            IAUSYL (NYL) = NINT (AAUSYL)
            IEUSYL (NYL) = ITUSYL (NYL) - IAUSYL (NYL) * 100
            IRECRD = IRECRD + 1
            IF ( LLNUYL (NYL) ) THEN
               READ (1) IGMUYL(NYL), (ENGMAX (IG), IG=1,IGMUYL(NYL)+1)
               IRECRD = IRECRD + 1
            ELSE
               IGMUYL (NYL) = 0
            END IF
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            IF ( IAUSYL (NYL) .EQ. 14 .OR. IAUSYL (NYL) .EQ. 15 .OR.
     &           IAUSYL (NYL) .EQ. 17 .OR. IAUSYL (NYL) .EQ. 34 ) THEN
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |  Small angle (exact!):
               IF ( AYLHGH (NYL) .LT. 0.3D+00 * PIPIPI ) THEN
                  TTLOW  = TAN (HLFHLF*AYLLOW(NYL))
                  TTHGH  = TAN (HLFHLF*AYLHGH(NYL))
                  DAYLBN (NYL) = TWOPIP * TWOTWO * ( TTHGH + TTLOW )
     &                      * ( TTHGH - TTLOW ) / ( ONEONE + TTLOW**2 )
     &                      / ( ONEONE + TTHGH**2 )
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |  Large angle:
               ELSE
                  DAYLBN (NYL) = TWOPIP * ( COS (DBLE(AYLLOW(NYL)))
     &                         - COS (DBLE(AYLHGH(NYL))) )
               END IF
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
*  |  |  |
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            ELSE
               DAYLBN (NYL) = AYLHGH (NYL) - AYLLOW (NYL)
            END IF
*  |  |  |
*  |  |  +-------------------------------------------------------------*
            NAYLBN(NYL)=1
******
            JX (IX) = NYL
***************************************************************
*********** Note KBUSYL is used for 0 address *****************
***************************************************************
            KBUSYL (NYL) = KLAST
            K0 = KLAST + 1
            K1 = NEYLBN (NYL) + K0 - 1
            K2 = K1 + IGMUYL (NYL)
            KLAST = K2
            READ (1) (GMSTOR(J), J = K0, K2)
            IRECRD = IRECRD + 1
            KMAX = MAX ( KMAX, K2 )
 4400    CONTINUE
 4600    CONTINUE
*        BACKSPACE 1
         REWIND 1
         DO IR = 1, IRECRD
            READ (1)
         END DO
         READ (1,ERR=9999) CHSTAT,IJX
         IRECRD = IRECRD + 1
         IF ( CHSTAT .NE. 'STATISTICS') GO TO 9999
         LSTATI = .TRUE.
         IX     = IX - 1
         IF ( IJX .NE. IX ) STOP 'IJX!=IX'
*  |  +----------------------------------------------------------------*
*  |  |
         IF ( IFILE .EQ. 1 ) THEN
            DO IIX = 1, IX
               NYL = JX (IIX)
               KFUSYL (NYL) = KLAST
               KLAST = KLAST + 2 * ( NEYLBN (NYL) + IGMUYL (NYL) )
               KMAX  = MAX ( KMAX, KLAST )
            END DO
            FILE1 = FILE
            LQOLD = LQ
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |
         ELSE
            NYL   = JX (IX)
            KLAST = KFUSYL (NYL) + 2 * ( NEYLBN (NYL) + IGMUYL (NYL) )
            KMAX  = MAX ( KMAX, KLAST )
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         KLAST = 0
*  |  +----------------------------------------------------------------*
*  |  |
         DO IIX = 1, IX
            NYL = JX (IIX)
            READ (1) TOTTHL (NYL), TOTEHL (NYL)
            IRECRD = IRECRD + 1
            TOTEHL (NYL) = ( TOTEHL (NYL) * TOTTHL (NYL) )**2
            TOTEHL (NYL) = TOTEHL (NYL) * DBLE (NBATCH-1)
     &                   + TOTTHL (NYL)**2
            TOTTOT (NYL) = TOTTOT (NYL) + DBLE (TOTTHL (NYL))
     &                   * ( DBLE (NCASE) + DBLE (MCASE) * 1.D+09 )
            TOTERR (NYL) = TOTERR (NYL) + DBLE (TOTEHL (NYL))
     &                   * ( DBLE (NCASE) + DBLE (MCASE) * 1.D+09 )
            READ (1) NETCCK,IGMUCK,ETCLCK,(EPGMCK(IG),
     &               IG=1,NETCCK+IGMUCK+1)
            IRECRD = IRECRD + 1
            IF ( NETCCK .NE. NEYLBN (NYL) )
     &         WRITE (*,*)' *** Error Netcck:',NETCCK,NEYLBN(NYL)
            IF ( IGMUCK .NE. IGMUYL (NYL) )
     &         WRITE (*,*)' *** Error Igmuck:',IGMUCK,IGMUYL(NYL)
            IF ( ABS (ETCLCK-EYLLOW(NYL)) .GT. 1.E-05*ABS(ETCLCK) )
     &         WRITE (*,*)' *** Error Etclck:',ETCLCK,EYLLOW(NYL)
            KK0 = 1
            KK1 = NEYLBN (NYL)
            KK2 = NEYLBN (NYL) + IGMUYL (NYL)
            READ (1) ( ACCUMH (J), J = KK1, KK0,-1 ),
     &               ( ACCUMH (J), J = KK1+1, KK2 )
            IRECRD = IRECRD + 1
            READ (1) ( ACCUEH (J), J = KK1, KK0,-1 ),
     &               ( ACCUEH (J), J = KK1+1, KK2 )
            IRECRD = IRECRD + 1
            DO J = KK0, KK2
               ACCUEH (J) = ( ACCUEH (J) * ACCUMH (J) )**2
               ACCUEH (J) = ACCUEH (J) * DBLE (NBATCH-1)
     &                    + ACCUMH (J)**2
               GDSTOR (KFUSYL(NYL)+J) = GDSTOR (KFUSYL(NYL)+J)
     &                                + DBLE (ACCUMH (J))
     &               * ( DBLE (NCASE) + DBLE (MCASE) * 1.D+09 )
               GBSTOR (KFUSYL(NYL)+J) = GBSTOR (KFUSYL(NYL)+J)
     &                                + DBLE (ACCUEH (J))
     &               * ( DBLE (NCASE) + DBLE (MCASE) * 1.D+09 )
            END DO
            KFUHLP = KFUSYL (NYL) + NEYLBN (NYL) + IGMUYL (NYL)
            READ (1) ( ACCUMH (J), J = KK1, KK0,-1 ),
     &               ( ACCUMH (J), J = KK1+1, KK2 )
            IRECRD = IRECRD + 1
            READ (1) ( ACCUEH (J), J = KK1, KK0,-1 ),
     &               ( ACCUEH (J), J = KK1+1, KK2 )
            IRECRD = IRECRD + 1
            DO J = KK0, KK2
               ACCUEH (J) = ( ACCUEH (J) * ACCUMH (J) )**2
               ACCUEH (J) = ACCUEH (J) * DBLE (NBATCH-1)
     &                    + ACCUMH (J)**2
               GDSTOR (KFUHLP+J) = GDSTOR (KFUHLP+J) + DBLE (ACCUMH (J))
     &               * ( DBLE (NCASE) + DBLE (MCASE) * 1.D+09 )
               GBSTOR (KFUHLP+J) = GBSTOR (KFUHLP+J) + DBLE (ACCUEH (J))
     &               * ( DBLE (NCASE) + DBLE (MCASE) * 1.D+09 )
            END DO
            K0 = KLAST + 1
            K1 = NEYLBN (NYL) + K0 - 1
            K2 = K1 + IGMUYL (NYL)
            KLAST = K2
            READ (1) ( GESTOR (J), J = K0, K2 )
            IRECRD = IRECRD + 1
            DO J = K0, K2
               GESTOR (J) = (GESTOR(J)*GMSTOR(J))**2
               GESTOR (J) = GESTOR(J) * DBLE(NBATCH-1) + GMSTOR(J)**2
               GDSTOR (J) = GDSTOR (J) + DBLE (GMSTOR(J))
     &                    * ( DBLE (NCASE) + DBLE (MCASE) * 1.D+09 )
               GBSTOR (J) = GBSTOR (J) + DBLE (GESTOR(J))
     &                    * ( DBLE (NCASE) + DBLE (MCASE) * 1.D+09 )
            END DO
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
         CLOSE (UNIT=1)
         GO TO 1111
      END IF
*  |
*  +-------------------------------------------------------------------*
      KLAST = 0
*  +-------------------------------------------------------------------*
*  |
      DO 7000 IIX = 1, IX
         NYL = JX (IIX)
         ACCUMP = ZERZER
         DO IG = 1, NEYLBN (NYL) + IGMUYL (NYL)
            ACCUMG (IG) = ZERZER
            ACUMCU (IG) = ZERZER
         END DO
         AFACT  = DBLE (NCASE) + DBLE (MCASE) * 1.D+09
         K0     = KLAST + 1
         K1     = NEYLBN (NYL) + K0 - 1
         K2     = K1 + IGMUYL (NYL)
         KLAST  = K2
         ENERG0 = EYLLOW (NYL)
*  |  +----------------------------------------------------------------*
*  |  |  Linear 1st variable binning
         IF ( IEUSYL (NYL) .GT. 0 ) THEN
            IF ( IEUSYL (NYL) .NE. 14 .AND. IEUSYL (NYL) .NE. 15 )
     &         DELTAE = DBLE (DEYLBN (NYL))
            EFACT  = DBLE (DEYLBN (NYL))
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |  Log 1st variable binning
         ELSE
            EFACT  = DBLE (DEYLBN (NYL))
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         ENEOLD = ENERG0
*  |  +----------------------------------------------------------------*
*  |  |
         DO 6540 IE = 1, NEYLBN (NYL) + IGMUYL (NYL)
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            IF ( IE .LE. NEYLBN (NYL) ) THEN
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |  Log 1st variable binning
               IF ( IEUSYL (NYL) .LT. 0 ) THEN
                  ENERGE = DBLE (ENERG0) * EFACT**IE
                  IF ( IEUSYL (NYL) .NE. 14 .AND.
     &                 IEUSYL (NYL) .NE. 15 .AND.
     &                 IEUSYL (NYL) .NE. 17 .AND.
     &                 IEUSYL (NYL) .NE. 34 ) THEN
                     DELTAE = DBLE (ENERGE - ENEOLD)
                  ELSE
                     DELTAE = TWOPIP * ( COS (DBLE(ENEOLD))
     &                      - COS (DBLE(ENERGE)) )
                  END IF
                  ENEOLD = ENERGE
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |  Linear angular 1st variable binning:
               ELSE IF ( IEUSYL (NYL) .EQ. 14 .OR.
     &                   IEUSYL (NYL) .EQ. 15 .OR.
     &                   IEUSYL (NYL) .EQ. 17 .OR.
     &                   IEUSYL (NYL) .EQ. 34 ) THEN
                  ENERGE = MIN ( DBLE (PIPIPI), DBLE (ENEOLD + EFACT) )
                  DELTAE = TWOPIP * ( COS (DBLE(ENEOLD))
     &                   - COS (DBLE(ENERGE)) )
                  ENEOLD = ENERGE
               END IF
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
*  |  |  |
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            ELSE
               IF ( IE-NEYLBN (NYL) .EQ. 1) ENEOLD = ENGMAX (1)
               ENERGE = ENGMAX (IE-NEYLBN(NYL)+1)
               DELTAE = DBLE(ENEOLD - ENERGE)
            END IF
*  |  |  |
*  |  |  +-------------------------------------------------------------*
            ENEOLD = ENERGE

            IF (DELTAE.LE. ZERZER) WRITE(*,*)' **** DELTAE',IE,NYL

*  |  |  +----------------------------------------------------------*
*  |  |  |
            IF ( IE .LE. NEYLBN (NYL) ) THEN
               IEA = IE + K0 - 1
*  |  |  |
*  |  |  +----------------------------------------------------------*
*  |  |  |
            ELSE
               IG  = IE - NEYLBN (NYL)
               IEA = IG + K1
            END IF
*  |  |  |
*  |  |  +----------------------------------------------------------*
            DSCORE = DBLE ( GMSTOR (IEA) ) * DELTAE
            DSCOR2 = DBLE ( GMSTOR (IEA) )
            ACCUMP = ACCUMP + DSCORE
            ACUMCU (IE) = ACUMCU (IE) + DSCORE
            ACCUMG (IE) = ACCUMG (IE) + DSCOR2
6540     CONTINUE
*  |  |
*  |  +----------------------------------------------------------------*
         TOTTOT (NYL) = TOTTOT (NYL) + ACCUMP * AFACT
         TOTERR (NYL) = TOTERR (NYL) + ACCUMP * ACCUMP * AFACT
*  |  +----------------------------------------------------------------*
*  |  |
         DO IE = NEYLBN (NYL) - 1, 1, -1
            ACUMCU (IE) = ACUMCU (IE) + ACUMCU (IE+1)
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
         IF ( NEYLBN(NYL) .GT. 0 )
     &   ACUMCU (NEYLBN(NYL)+1) = ACUMCU (NEYLBN(NYL)+1) + ACUMCU (1)
*  |  +----------------------------------------------------------------*
*  |  |
         DO IE = NEYLBN (NYL) + 2, NEYLBN (NYL) + IGMUYL (NYL)
            ACUMCU (IE) = ACUMCU (IE) + ACUMCU (IE-1)
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
         IF ( LLNUYL (NYL) .AND.
     &      ABS ( ACUMCU ( NEYLBN (NYL) + IGMUYL (NYL) ) - ACCUMP ) .GT.
     &      1.D-13 * ACCUMP ) WRITE (*,*)' *** ERROR', NYL,
     &      ACUMCU ( NEYLBN (NYL) + IGMUYL (NYL) ), ACCUMP
*  |  +----------------------------------------------------------------*
*  |  |
         DO IE = 1, NEYLBN (NYL) + IGMUYL (NYL)
            KK0 = KFUSYL (NYL) + IE
            GDSTOR (KK0) = GDSTOR (KK0) + ACCUMG (IE) * AFACT
            GBSTOR (KK0) = GBSTOR (KK0) + ACCUMG (IE) * ACCUMG (IE) *
     &                     AFACT
            KK0 = KK0 + NEYLBN (NYL) + IGMUYL (NYL)
            GDSTOR (KK0) = GDSTOR (KK0) + ACUMCU (IE) * AFACT
            GBSTOR (KK0) = GBSTOR (KK0) + ACUMCU (IE) * ACUMCU (IE) *
     &                     AFACT
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
         KMAX = MAX ( KMAX, KK0 )
7000  CONTINUE
*  |
*  +-------------------------------------------------------------------*
      GO TO 1111
2000  CONTINUE
*  +-------------------------------------------------------------------*
*  |
      DO J = 1, KMAX
         GDSTOR (J) = GDSTOR (J)
     &              / ( DBLE (NCTOT) + DBLE (MCTOT) * 1.D+09 )
         GBSTOR (J) = GBSTOR (J)
     &              / ( DBLE (NCTOT) + DBLE (MCTOT) * 1.D+09 )
      END DO
*  |
*  +-------------------------------------------------------------------*
*  +-------------------------------------------------------------------*
*  |
      DO IIX = 1, IX
         NYL = JX (IIX)
         TOTTOT (NYL) = TOTTOT (NYL)
     &                / ( DBLE (NCTOT) + DBLE (MCTOT) * 1.D+09 )
         TOTERR (NYL) = TOTERR (NYL)
     &                / ( DBLE (NCTOT) + DBLE (MCTOT) * 1.D+09 )
      END DO
*  |
*  +-------------------------------------------------------------------*
*  +-------------------------------------------------------------------*
*  |  Real sum of Mbatch runs:
      IF ( MBATCH - 1 .GT. 0 .AND. ( .NOT. LSTATI .OR. IFILE .GT. 1 )
     &    ) THEN
         DO J = 1, KMAX
            GBHELP     = SQRT (GBSTOR(J))
            GBSTOR (J) = ( ( GBHELP - GDSTOR (J) )
     &                   * ( GBHELP + GDSTOR (J) ) )
     &                 / DBLE (MBATCH-1)
            IF ( ABS (GBSTOR(J)) .LT. 1.D-12 * GDSTOR(J)**2 )
     &         GBSTOR (J) = ZERZER
            IF ( GDSTOR(J) .GT. ZERZER ) THEN
               GBSTOR (J) = MIN ( SQRT (GBSTOR(J)) / GDSTOR(J),
     &                            0.99D+00 )
            END IF
         END DO
         DO IIX = 1, IX
            NYL = JX (IIX)
            TOTHLP       = SQRT (TOTERR(NYL))
            TOTERR (NYL) = ( ( TOTHLP - TOTTOT (NYL) )
     &                     * ( TOTHLP + TOTTOT (NYL) ) )
     &                   / DBLE (MBATCH-1)
            IF ( ABS (TOTERR(NYL)) .LT. 1.D-12 * TOTTOT (NYL)**2 )
     &              TOTERR (NYL) = ZERZER
            IF ( TOTTOT (NYL) .GT. ZERZER ) THEN
               TOTERR (NYL) = MIN ( SQRT ( TOTERR(NYL) ) / TOTTOT(NYL),
     &                              0.99D+00 )
            END IF
         END DO
*  |
*  +-------------------------------------------------------------------*
*  |  Only one run, check if it is already a sum with statistics:
      ELSE IF ( LSTATI ) THEN
         OPEN ( UNIT=11, FILE=FILE1, STATUS='OLD', FORM='UNFORMATTED' )
         READ (11)
         READ (11)
         DO 7050 IIX = 1, IX
            NYL = JX (IIX)
            READ (11)
            IF ( LLNUYL (NYL) ) READ (11)
            READ (11)
 7050    CONTINUE
         READ (11) CHSTAT,IJX
         IF ( CHSTAT .NE. 'STATISTICS' .OR. IJX .NE. IX ) THEN
            WRITE (*,*) IJX, IX
            STOP 'STATI'
         END IF
         KLAST = 0
         DO 7400 IIX = 1, IX
            NYL = JX (IIX)
            READ (11) TOTA, TOTB
            IF ( ABS (TOTA-TOTTOT(NYL)) .GT. 1.D-05*TOTA )
     &         WRITE (*,*)' *** Error Tota:',TOTA,TOTTOT(NYL), NYL
            TOTERR (NYL) = TOTB
            READ (11) NETCCK,IGMUCK,ETCLCK,(EPGMCK(IG),
     &                IG=1,NETCCK+IGMUCK+1)
            IF ( NETCCK .NE. NEYLBN (NYL) )
     &         WRITE (*,*)' *** Error Netcck:',NETCCK,NEYLBN(NYL)
            IF ( IGMUCK .NE. IGMUYL (NYL) )
     &         WRITE (*,*)' *** Error Igmuck:',IGMUCK,IGMUYL(NYL)
            IF ( ABS (ETCLCK-EYLLOW(NYL)) .GT. 1.E-05*ABS(ETCLCK) )
     &         WRITE (*,*)' *** Error Etclck:',ETCLCK,EYLLOW(NYL)
            KK0 = 1
            KK1 = NEYLBN (NYL)
            KK2 = NEYLBN (NYL) + IGMUYL (NYL)
            READ (11) ( ACCUMH (J), J = KK1, KK0,-1 ),
     &                ( ACCUMH (J), J = KK1+1, KK2 )
            DO 7100 J = KK0, KK1
               IF ( ABS (ACCUMH(J)-GDSTOR(J+KFUSYL(NYL))) .GT.
     &              1.D-05*ACCUMH(J) )
     &         WRITE (*,*)' *** Error Gdstor:',ACCUMH(J),
     &                      GDSTOR(J+KFUSYL(NYL))
 7100       CONTINUE
            READ (11) ( ACCUMH (J), J = KK1, KK0,-1 ),
     &                ( ACCUMH (J), J = KK1+1, KK2 )
            DO 7150 J = KK0, KK2
               GBSTOR (J+KFUSYL(NYL)) = ACCUMH (J)
 7150       CONTINUE
            KFUHLP = KFUSYL (NYL) + NEYLBN (NYL) + IGMUYL (NYL)
            READ (11) ( ACCUMH (J), J = KK1, KK0,-1 ),
     &                ( ACCUMH (J), J = KK1+1, KK2 )
            DO 7200 J = KK0, KK2
               IF ( ABS (ACCUMH(J)-GDSTOR(J+KFUHLP)) .GT.
     &              1.D-05*ACCUMH(J) )
     &         WRITE (*,*)' *** Error Gdstor_2:',ACCUMH(J),
     &                      GDSTOR(J+KFUHLP)
 7200       CONTINUE
            READ (11) ( ACCUMH (J), J = KK1, KK0,-1 ),
     &                ( ACCUMH (J), J = KK1+1, KK2 )
            DO 7250 J = KK0, KK2
               GBSTOR (J+KFUHLP) = ACCUMH (J)
 7250       CONTINUE
            K0 = KLAST + 1
            K1 = NEYLBN (NYL) + K0 - 1
            K2 = K1 + IGMUYL (NYL)
            KLAST = K2
            READ (11) ( GMSTOR(J), J = K0, K2 )
            DO 7300 J = K0, K2
               GBSTOR (J) = DBLE (GMSTOR(J))
 7300       CONTINUE
 7400    CONTINUE
         KLAST = 0
         CLOSE ( UNIT=11 )
*  |
*  +-------------------------------------------------------------------*
*  |  Only one plain run, no possibility of making statistics:
      ELSE
         DO J = 1, KMAX
            IF (GDSTOR (J) .GT. ZERZER) THEN
               GBSTOR (J) = 0.99D+00
            END IF
         END DO
         DO IIX = 1, IX
            NYL = JX (IIX)
            IF ( TOTTOT (NYL) .GT. ZERZER ) THEN
               TOTERR (NYL) = 0.99D+00
            END IF
         END DO
      END IF
*  |
*  +-------------------------------------------------------------------*
*  +-------------------------------------------------------------------*
*  |  If 1 run already with statistics, no need to create a new output
      IF ( .NOT. LSTATI .OR. IFILE .GT. 1 ) THEN
         RUNTIM = '      ***** Sum file *****'
         WRITE(*,*)' Type the output file name:'
         READ (*,'(A)')FILE
         LQ     = LNNBLN(FILE)
* Start_VAX_seq
*        OPEN (UNIT=1,FILE=FILE,STATUS='NEW',FORM='UNFORMATTED')
* End_VAX_seq
* Start_UNIX_seq
         LQS = LLOCCR(FILE,'/') + 1
*  |  +----------------------------------------------------------------*
*  |  |  Be sure to add the default suffix if no one is specified:
         IF ( LQ .GE. LQS ) THEN
            LQD = INDEX (FILE(LQS:LQ),'.')
            IF ( LQD .EQ. LQ-LQS+1 ) THEN
               LQD = 0
               LQ  = LQ - 1
            END IF
            IF ( LQD .LE. 0 ) THEN
               FILE = FILE(1:LQ)//'.yie'
               LQ   = LQ+4
            END IF
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         OPEN (UNIT=1,FILE=FILE,STATUS='UNKNOWN',FORM='UNFORMATTED')
* End_UNIX_seq
         WRITE (1) RUNTIT,RUNTIM,SNGL(WCTOT),NCTOT,MCTOT,MBATCH
         WRITE (1) IJUSYL, JTUSYL, PUSRYL, SQSUYL, UUSRYL, VUSRYL,
     &             WUSRYL
*  |
*  +-------------------------------------------------------------------*
*  |
      ELSE
         FILE = FILE1
         LQ   = LQOLD
      END IF
*  |
*  +-------------------------------------------------------------------*
      KLAST = 0
*  +-------------------------------------------------------------------*
*  |
      DO 6800 JJ=1,IX
         NYL=JX(JJ)
         IF ( .NOT. LSTATI .OR. IFILE .GT. 1 )
     &   WRITE (1)    NYL, TITUYL(NYL), ITUSYL(NYL),
     &                IXUSYL(NYL), IDUSYL(NYL),
     &                NR1UYL(NYL), NR2UYL(NYL), USNRYL (NYL),
     &                SGUSYL(NYL), LLNUYL(NYL),
     &                EYLLOW(NYL), EYLHGH(NYL),
     &                NEYLBN(NYL), DEYLBN(NYL),
     &                AYLLOW(NYL), AYLHGH(NYL)
         IF ( LLNUYL (NYL) .AND. ( .NOT. LSTATI .OR. IFILE .GT. 1 ))THEN
            WRITE (1) IGMUYL(NYL), (ENGMAX (IG), IG=1,IGMUYL(NYL)+1)
         END IF
         K0 = KLAST + 1
         K1 = NEYLBN (NYL) + K0 - 1
         K2 = K1 + IGMUYL (NYL)
         KLAST = K2
         IF ( .NOT. LSTATI .OR. IFILE .GT. 1 )
     &      WRITE (1) (SNGL(GDSTOR(J)), J = K0, K2)
6800  CONTINUE
*  |
*  +-------------------------------------------------------------------*
*  Now write the statistical data, first a record to flag them
      IF ( .NOT. LSTATI .OR. IFILE .GT. 1 )
     &   WRITE (1) 'STATISTICS',IX
*  Then the errors concerning the binning data
      KLAST = 0
      INDX  = INDLST ( FILE, '.' )
      IF ( INDX .GT. 0 ) THEN
* Start_UNIX_seq
         OPEN (UNIT=3,FILE=FILE(1:INDX-1)//'_sum.lis',STATUS='UNKNOWN')
* End_UNIX_seq
* Start_VAX_seq
*        OPEN (UNIT=3,FILE=FILE(1:INDX-1)//'_sum.lis',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
* End_VAX_seq
      ELSE
* Start_UNIX_seq
         OPEN (UNIT=3,FILE=FILE(1:LQ)//'_sum.lis',STATUS='UNKNOWN')
* End_UNIX_seq
* Start_VAX_seq
*        OPEN (UNIT=3,FILE=FILE(1:LQ)//'_sum.lis',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
* End_VAX_seq
      END IF
      IF ( INDX .GT. 0 ) THEN
* Start_UNIX_seq
         OPEN (UNIT=4,FILE=FILE(1:INDX-1)//'_tab.lis',STATUS='UNKNOWN')
* End_UNIX_seq
* Start_VAX_seq
*        OPEN (UNIT=4,FILE=FILE(1:INDX-1)//'_tab.lis',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
* End_VAX_seq
      ELSE
* Start_UNIX_seq
         OPEN (UNIT=4,FILE=FILE(1:LQ)//'_tab.lis',STATUS='UNKNOWN')
* End_UNIX_seq
* Start_VAX_seq
*        OPEN (UNIT=4,FILE=FILE(1:LQ)//'_tab.lis',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
* End_VAX_seq
      END IF
      WRITE (3,*)' **** ',FILE(1:LQ),' ****'
      WRITE (3,*)
      WRITE (3,*)
      WRITE (CHCTOT,'(F15.0)') DBLE (NCTOT) + DBLE (MCTOT) * 1.D+09
      WRITE (3,'(A,2X,A14)')'  Total primaries run:',CHCTOT
      WRITE (3,*)
      WRITE (3,*)'   Total weight of the primaries run:',SNGL(WCTOT)
      WRITE (3,*)
*  +-------------------------------------------------------------------*
*  |
      DO 6850 JJ = 1, IX
         NYL    = JX (JJ)
         WRITE (1) SNGL(TOTTOT(NYL)),SNGL(TOTERR(NYL))
         ENERG0 = EYLLOW (NYL)
*  |  +----------------------------------------------------------------*
*  |  |
         IF ( IEUSYL (NYL) .LT. 0 ) THEN
            EFACT  = DBLE (DEYLBN (NYL))
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |
         ELSE
            DELTAE = DBLE (DEYLBN (NYL))
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         ENEOLD = ENERG0
*  |  +----------------------------------------------------------------*
*  |  |
         DO IE = 1, NEYLBN (NYL)
            IG = NEYLBN (NYL) - IE + 1
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            IF ( IEUSYL (NYL) .LT. 0 ) THEN
               ENERGE = ENERG0 * EFACT**IE
               DELTAE = DBLE (ENERGE - ENEOLD)
*  |  |  |
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            ELSE
               ENERGE = EYLLOW (NYL) + DEYLBN (NYL) * FLOAT (IE)
            END IF
*  |  |  |
*  |  |  +-------------------------------------------------------------*
            EPGMAX (IG) = ENERGE
            ENEOLD = ENERGE
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
*  |  +----------------------------------------------------------------*
*  |  |
         DO IG = NEYLBN (NYL) + 1, NEYLBN (NYL) + IGMUYL (NYL)
            EPGMAX (IG) = ENGMAX (IG-NEYLBN(NYL))
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
*  |  +----------------------------------------------------------------*
*  |  |
         IF (LLNUYL (NYL)) THEN
            EPGMAX (NEYLBN(NYL)+IGMUYL(NYL)+1) = ENGMAX (IGMUYL(NYL)+1)
            IF (IDUSYL(NYL) .NE. 8)EPGMAX (NEYLBN(NYL)+1) = EYLLOW (NYL)
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |
         ELSE
            EPGMAX (NEYLBN(NYL)+1) = EYLLOW (NYL)
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         IF ( .NOT. LSTATI .OR. IFILE .GT. 1 )
     &   WRITE (1) NEYLBN(NYL),IGMUYL(NYL),EYLLOW(NYL),(EPGMAX(IG),
     &             IG=1,NEYLBN(NYL)+IGMUYL(NYL)+1)
         KK0 = KFUSYL (NYL) + 1
         KK1 = KFUSYL (NYL) + NEYLBN (NYL)
         KK2 = KFUSYL (NYL) + NEYLBN (NYL) + IGMUYL (NYL)

         WRITE (3,*)
         WRITE (3,*) ' Detector n: ', NYL,'  ',TITUYL (NYL)
         WRITE (3,*) '    (User norm:       ',USNRYL(NYL),' '
         WRITE (3,*) '     sigma:           ',SGUSYL(NYL),' mb'
         WRITE (3,*) '     distr. scored:',IDUSYL(NYL),'   ,'
         WRITE (3,*) '     from reg.',NR1UYL(NYL),' to reg. ',
     &                     NR2UYL(NYL),' )'
         IF ( LLNUYL (NYL) ) WRITE (3,9400) IGMUYL(NYL)
*  |  +----------------------------------------------------------------*
*  |  |  Linear 1st variable binning
         IF ( IEUSYL (NYL) .GT. 0 ) THEN
            WRITE (3,8320) EYLLOW(NYL), EYLHGH(NYL),
     &                     NEYLBN(NYL), DEYLBN(NYL)
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |  Log 1st variable binning
         ELSE
            WRITE (3,8340) EYLLOW(NYL), EYLHGH(NYL),
     &                     NEYLBN(NYL), DEYLBN(NYL)
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         ENEOLD = ENERG0
         WRITE (3,8360) AYLLOW (NYL), AYLHGH (NYL)
         WRITE (3,8370) CHVARI (ABS(IEUSYL(NYL)))
         WRITE (3,8380) CHVARI (ABS(IAUSYL(NYL)))
         IXM = MOD ( IXUSYL(NYL), 100 )
*        JXM = IXM / 10
*        IXM = MOD ( IXM, 10 ) + 3 * MIN ( JXM, 1 ) + MAX ( JXM - 1, 0 )
         WRITE (3,8400) ( CHYILD(I,IXM), I = 1, 3 )
         IF ( IEUSYL (NYL) .EQ. 14 .OR. IEUSYL (NYL) .EQ. 15 .OR.
     &        IEUSYL (NYL) .EQ. 17 .OR. IEUSYL (NYL) .EQ. 34 )
     &      WRITE (3,8470)
         IF ( IAUSYL (NYL) .EQ. 14 .OR. IAUSYL (NYL) .EQ. 15 .OR.
     &        IAUSYL (NYL) .EQ. 17 .OR. IAUSYL (NYL) .EQ. 34 )
     &      WRITE (3,8480)
         WRITE (3,*)
         WRITE (3,10000)'    Tot. response (integrated over x1)',
     &                   SNGL(TOTTOT(NYL)),'  +/-',
     &                   SNGL(TOTERR(NYL))*100.0,' %'
         WRITE (3,10000)
     &            '    ( --> Tot. response without normaliz. ',
     &                   SNGL(TOTTOT(NYL)*USNRYL(NYL)),'  +/-',
     &                   SNGL(TOTERR(NYL))*100.0,' % )'
         WRITE (3,*)
         WRITE (3,*)
     &         ' **** Double differ. yield as a function of x1 ****',
     &         ' ****          (not integrated)                ****'
         WRITE (3,*)
         WRITE (3,*)
         WRITE (3,*) ' x1 boundaries:'
         WRITE (3,*)
         WRITE (3,10002) ( EPGMAX (J), J=1,NEYLBN(NYL) )
         WRITE (3,10003) ' Lowest x1 boundary      :',
     &                     EPGMAX (NEYLBN(NYL)+1)
         WRITE (3,*)
         WRITE (3,*) ' Yield (double differ.):'
         WRITE (3,*)
         WRITE (3,10001) ( SNGL (GDSTOR(J)),' +/- ',
     &                     SNGL (GBSTOR(J)*100.D+00),' % ',
     &                     J = KK1, KK0,-1 )
         IF ( LLNUYL (NYL) ) THEN
            IF (IDUSYL(NYL) .NE. 8)
     &         EPGMAX (NEYLBN(NYL)+1) = ENGMAX (1)
            WRITE (3,*)
            WRITE (3,*) ' x1 boundaries:'
            WRITE (3,*)
            WRITE (3,10002) ( EPGMAX (J), J=NEYLBN(NYL)+1,
     &                        NEYLBN(NYL)+IGMUYL(NYL) )
            WRITE (3,10003) ' Lowest x1 boundary      :',
     &                        EPGMAX (NEYLBN(NYL)+IGMUYL(NYL)+1)
            WRITE (3,*)
            WRITE (3,*) ' Yield (double differ.):'
            WRITE (3,*)
            WRITE (3,10001) ( SNGL (GDSTOR(J)),' +/- ',
     &                        SNGL (GBSTOR(J)*100.D+00),' % ',
     &                        J = KK1+1, KK2 )
            IF (IDUSYL(NYL) .NE. 8)
     &         EPGMAX (NEYLBN(NYL)+1) = EYLLOW (NYL)
         END IF
         IF ( .NOT. LSTATI .OR. IFILE .GT. 1 )
     &   WRITE (1) ( SNGL (GDSTOR(J)), J = KK1, KK0,-1 ),
     &             ( SNGL (GDSTOR(J)), J = KK1+1, KK2 )
         IF ( .NOT. LSTATI .OR. IFILE .GT. 1 )
     &   WRITE (1) ( SNGL (GBSTOR(J)), J = KK1, KK0,-1 ),
     &             ( SNGL (GBSTOR(J)), J = KK1+1, KK2 )
*****   here to write on "table" file
         IF (NYL .GT. 1) THEN
*****       write 2 empty lines so gnuplot can use it with 'index #'
            WRITE (4,*)
            WRITE (4,*)
         END IF
         WRITE (4,*) '# Detector n: ', NYL,'  ',TITUYL (NYL)
         WRITE (4,*) '# N. of x1 intervals' ,NEYLBN(NYL)+IGMUYL(NYL)
         IF ( LLNUYL (NYL) ) THEN
            IF (IDUSYL(NYL) .NE. 8)
     &         EPGMAX (NEYLBN(NYL)+1) = ENGMAX (1)
            KK = KK2 + 1
            DO 6001 J=NEYLBN(NYL)+IGMUYL(NYL), NEYLBN(NYL)+1,-1
               KK = KK - 1
               WRITE(4,9501) EPGMAX (J+1), EPGMAX (J),
     &              SNGL (GDSTOR(KK)), SNGL (GBSTOR(KK)*100.D+00)
            IF (IDUSYL(NYL) .NE. 8)
     &         EPGMAX (NEYLBN(NYL)+1) = EYLLOW (NYL)
 6001       CONTINUE
         END IF
         KK=KK0-1
         DO 6002 J=NEYLBN(NYL),1,-1
            KK = KK + 1
            WRITE(4,9501) EPGMAX (J+1), EPGMAX (J),
     &                    SNGL (GDSTOR(KK)), SNGL (GBSTOR(KK)*100.D+00)
 6002    CONTINUE
******  end write on tab file
         KK0 = KFUSYL (NYL) + NEYLBN (NYL) + IGMUYL (NYL) + 1
         KK1 = KFUSYL (NYL) + 2 * NEYLBN (NYL) + IGMUYL (NYL)
         KK2 = KFUSYL (NYL) + 2 * ( NEYLBN (NYL) + IGMUYL (NYL) )
         WRITE (3,*)
         WRITE (3,*)
         WRITE (3,*)
     &         ' **** Cumulative Yields as a function of x1 ****',
     &         ' ****   (integrated over x1, differ. in x2) ****'
         WRITE (3,*)
         WRITE (3,*)
         WRITE (3,*) ' x1 boundaries:'
         WRITE (3,*)
         WRITE (3,*)
         WRITE (3,10002) ( EPGMAX (J), J=1,NEYLBN(NYL) )
         WRITE (3,10003) ' Lowest x1 boundary      :',
     &                     EPGMAX (NEYLBN(NYL)+1)
         WRITE (3,*)
         WRITE (3,*) ' Cumul. Yield:'
         WRITE (3,*)
         WRITE (3,10001) ( SNGL (GDSTOR(J)),' +/- ',
     &                     SNGL (GBSTOR(J)*100.D+00),' % ',
     &                     J = KK1, KK0,-1 )
         IF ( LLNUYL (NYL) ) THEN
            IF (IDUSYL(NYL) .NE. 8)
     &         EPGMAX (NEYLBN(NYL)+1) = ENGMAX (1)
            WRITE (3,*)
            WRITE (3,*) ' x1 boundaries:'
            WRITE (3,*)
            WRITE (3,10002) ( EPGMAX (J), J=NEYLBN(NYL)+1,
     &                        NEYLBN(NYL)+IGMUYL(NYL) )
            WRITE (3,10003) ' Lowest x1 boundary      :',
     &                        EPGMAX (NEYLBN(NYL)+IGMUYL(NYL)+1)
            WRITE (3,*)
            WRITE (3,*)
            WRITE (3,*) ' Cumul. Yield:'
            WRITE (3,*)
            WRITE (3,10001) ( SNGL (GDSTOR(J)),' +/- ',
     &                        SNGL (GBSTOR(J)*100.D+00),' % ',
     &                        J = KK1+1, KK2 )
            IF (IDUSYL(NYL) .NE. 8)
     &         EPGMAX (NEYLBN(NYL)+1) = EYLLOW (NYL)
         END IF
         IF ( .NOT. LSTATI .OR. IFILE .GT. 1 )
     &   WRITE (1) ( SNGL (GDSTOR(J)), J = KK1, KK0,-1 ),
     &             ( SNGL (GDSTOR(J)), J = KK1+1, KK2 )
         IF ( .NOT. LSTATI .OR. IFILE .GT. 1 )
     &   WRITE (1) ( SNGL (GBSTOR(J)), J = KK1, KK0,-1 ),
     &             ( SNGL (GBSTOR(J)), J = KK1+1, KK2 )
         K0 = KLAST + 1
         K1 = NEYLBN (NYL) + K0 - 1
         K2 = K1 + IGMUYL (NYL)
         KLAST = K2
         IF ( .NOT. LSTATI .OR. IFILE .GT. 1 )
     &   WRITE (1) (SNGL(GBSTOR(J)), J = K0, K2)
6850  CONTINUE
*  |
*  +-------------------------------------------------------------------*
      IF ( .NOT. LSTATI .OR. IFILE .GT. 1 ) CLOSE (UNIT=1)
 8320 FORMAT (6X,'linear 1st variable (x1) binning from  ',1P,E11.4,
     &        ' to ',E11.4,' ',0P,I5,
     &        ' bins (',1P,E11.4,' wide)')
 8340 FORMAT (6X,'logar. 1st variable (x1) binning from ',1P,E11.4,
     &        ' to ',E11.4,' ',0P,I5,
     &        ' bins (ratio :',1P,E11.4,')')
 8360 FORMAT (6X,'2nd variable (x2) ranges from  ',1P,E11.4,
     &        ' to ',E11.4)
 8370 FORMAT (6X,'1st variable (x1) is:',A30)
 8380 FORMAT (6X,'2nd variable (x2) is:',A30)
 8400 FORMAT (6X,'The scored double differential yield is ',
     &        '(normalized per primary particle):',/,
     &        3(7X,A49,/))
 8470 FORMAT (6X,'1st variable (x1) differential is a solid angle one')
 8480 FORMAT (6X,'2nd variable (x2) differential is a solid angle one')
 9400 FORMAT (6X,'low energy neutrons scored from group 1 to group ',
     &        I5)
 9501 FORMAT (4(1P,1X,G14.7))
10000 FORMAT ((1X,A,1X,1P,G14.7,A,1X,G14.7,A))
10001 FORMAT ((2(1X,1P,G14.7,A,1X,G14.7,A)))
10002 FORMAT ((5(1X,1P,G14.7)))
10003 FORMAT ((1X,A,1X,1P,G14.7))

      DEALLOCATE(GMSTOR)
      DEALLOCATE(GDSTOR)
      DEALLOCATE(GBSTOR)
      DEALLOCATE(GESTOR)
      DEALLOCATE(ENGMAX)
      DEALLOCATE(EPGMAX)
      DEALLOCATE(EPGMCK)
      DEALLOCATE(ACCUMG)
      DEALLOCATE(ACUMCU)
      DEALLOCATE(ACCUMH)
      DEALLOCATE(ACCUEH)

      STOP
 9999 CONTINUE
      STOP 'UNKNOWN RECORD TYPE'
      END

*
*=== Lnnbln ===========================================================*
*
      INTEGER FUNCTION LNNBLN (CARD)

*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 2003-2019:  CERN & INFN                            *
*     All Rights Reserved.                                             *
*                                                                      *
*     Return the length of character variable CARD up to the last non  *
*     blank character                                                  *
*                                                                      *
*     Created on 21 September 1992 by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*----------------------------------------------------------------------*
*
      CHARACTER CARD*(*)
*
      LENGTH = LEN (CARD)
      DO 100 LQ = LENGTH, 1, -1
         IF ( CARD (LQ:LQ).NE.' ') GO TO 200
  100 CONTINUE
      LQ     = 0
  200 CONTINUE
      LNNBLN = LQ
      RETURN
*=== End of function Lnnbln ===========================================*
      END

*
*=== Lloccr ===========================================================*
*
      INTEGER FUNCTION LLOCCR ( CARD, CWSH )

*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 2003-2019:  CERN & INFN                            *
*     All Rights Reserved.                                             *
*                                                                      *
*     Return the length of character variable CARD up to the last      *
*     (included) occurrence of character CWSH (0 if not found)         *
*                                                                      *
*     Created on 21 September 1992 by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*----------------------------------------------------------------------*
*
      CHARACTER CARD*(*), CWSH*1
*
      LENGTH = LEN (CARD)
      DO 100 LQ = LENGTH, 1, -1
         IF ( CARD (LQ:LQ) .EQ. CWSH ) GO TO 200
  100 CONTINUE
      LQ     = 0
  200 CONTINUE
      LLOCCR = LQ
      RETURN
*=== End of function Lloccr ===========================================*
      END
