*
*=== Usxrea ===========================================================*
*
      PROGRAM USXREA
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 2003-2019:  CERN & INFN                            *
*     Copyright (C) 2020-    :  CERN                                   *
*     All Rights Reserved.                                             *
*                                                                      *
*     USer boundary X (crossing) scoring REAding:                      *
*                                                                      *
*     Created on 03 November 1989  by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*     Modified on  27 June 2021                                        *
*                                                                      *
*----------------------------------------------------------------------*
*
      CHARACTER FILE*232,RUNTIT*80,RUNTIM*32,CHSTAT*10,CHCTOT*15
*
*=== usrbdx ===========================================================*
*
*----------------------------------------------------------------------*
*     Module USRBDX:                                                   *
*     A. Ferrari : user defined boundary crossing scoring              *
*                                                                      *
*     Up to MXUSBX user defined bdrx are allowed                       *
*            lusbdx = logical flag, .true. if at least 1 user defined  *
*                     bdrx is used                                     *
*            nusrbx = number of user defined bdrx used                 *
*            itusbx = type of binning: 1 = linear energy, linear angle *
*                     2 = linear energy, logarithmic angle, -1 = loga- *
*                     rithmic energy, linear angle, -2 = logarithmic   *
*                     energy, logarithmic angle                        *
*            idusbx = distribution to be scored: the usual values are  *
*                     allowed                                          *
*            nr1usx = first region                                     *
*            nr2usx = second region                                    *
*            ausbdx = area (cm**2) of the detector                     *
*            lwusbx = one way if false, two ways if true               *
*            lfusbx = current if false, fluence if true                *
*            llnusx = no low energy neutron scoring if false, yes if   *
*                     true                                             *
*            titusx = bdrx name                                        *
*            ipusbx = logical unit to print the results on: formatted  *
*                     if > 0, unformatted if < 0                       *
*            igmusx = maximum low energy neutron group to be scored    *
*            kbusbx = initial location in blank common of the consi-   *
*                     dered bdrx (real*8 address)                      *
*                     I T   I S   F O R  A D D R E S S  OF  1 !!!!!    *
*            nebxbn = number of energy intervals                       *
*            nabxbn = number of angular intervals                      *
*     ebxlow/ebxhgh = minimum and maximum energies                     *
*     abxlow/abxhgh = minimum and maximum angle (steradian)            *
*            debxbn = energy bin width                                 *
*            dabxbn = angular (steradian) bin width                    *
*                                                                      *
*----------------------------------------------------------------------*
*
      PARAMETER ( MXUSBX = 1100 )
      LOGICAL LUSBDX, LFUSBX, LWUSBX, LLNUSX
      CHARACTER*10 TITUSX
      COMMON /UUSRBX/ EBXLOW(MXUSBX), EBXHGH(MXUSBX), ABXLOW(MXUSBX),
     &                ABXHGH(MXUSBX), DEBXBN(MXUSBX), DABXBN(MXUSBX),
     &                AUSBDX(MXUSBX),
     &                NEBXBN(MXUSBX), NABXBN(MXUSBX), NR1USX(MXUSBX),
     &                NR2USX(MXUSBX), ITUSBX(MXUSBX), IDUSBX(MXUSBX),
     &                KBUSBX(MXUSBX), IPUSBX(MXUSBX), IGMUSX(MXUSBX),
     &                LFUSBX(MXUSBX), LWUSBX(MXUSBX), LLNUSX(MXUSBX),
     &                NUSRBX, LUSBDX
      COMMON /UUSXCH/ TITUSX(MXUSBX)

      PARAMETER ( PIPIPI = 3.141592653589793D+00 )
      PARAMETER ( TWOPIP = 2.D+00 * PIPIPI )

      REAL,DIMENSION(:),ALLOCATABLE,TARGET :: GMSTOR
      REAL,DIMENSION(:),ALLOCATABLE :: GBSTOR, OMGMAX, ENGMAX, EPGMAX
      REAL,DIMENSION(:),POINTER :: GDSTOR

      DIMENSION TOTTOT (MXUSBX), TOTERR (MXUSBX), KFUSBX (MXUSBX),
     &          JX (MXUSBX)

      DOUBLE PRECISION DOMEGA

      EQUIVALENCE (ANNX,NNX)
*
*     Note that all addresses are for zero index (Kbusbx,Kfusbx)
*
      MAXSIZE = 0
      MNEBXBN = 0
      MIGMUSX = 0
      MNABXBN = 0

      WRITE(*,'('' Type the input file: '',$)')
      READ (*,'(A)')FILE
      LQ     = LNNBLN(FILE)
      IF (LQ .LE. 0) GO TO 2300
      OPEN (UNIT=1,FILE=FILE,STATUS='OLD',FORM='UNFORMATTED')
      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,ERR=120) RUNTIT,RUNTIM,WEIPRI,NCASE
      IRECRD = IRECRD + 1
      WRITE(*,*)RUNTIT
      WRITE(*,*)RUNTIM
      WRITE(*,*)WEIPRI
      WRITE(*,*)NCASE
      MCASE  = 0
      NBATCH = 1
      GO TO 200
 120  CONTINUE
      REWIND 1
      IRECRD = 0
      READ (1,ERR=150) RUNTIT,RUNTIM,NCASE
      IRECRD = IRECRD + 1
      WRITE(*,*)RUNTIT
      WRITE(*,*)RUNTIM
      WRITE(*,*)NCASE
      WEIPRI = REAL(NCASE)
      MCASE  = 0
      NBATCH = 1
      GO TO 200
 150  CONTINUE
      REWIND 1
      IRECRD = 0
      WRITE(*,*)' How many particle for this run?'
      READ(*,*)NCASE
      WEIPRI = REAL(NCASE)
      MCASE  = 0
      NBATCH = 1
 200  CONTINUE
      NCTOT  = NCASE
      MCTOT  = MCASE
      IF ( NCTOT .GE. 1000000000 ) THEN
         NCTOT = NCTOT - 1000000000
         MCTOT = MCTOT + 1
      END IF
      WCTOT  = DBLE (WEIPRI)
      MBATCH = NBATCH

*  +-------------------------------------------------------------------*
*  | First pass to determine array size needed
      IF (MAXSIZE .EQ. 0) THEN
         IRKEEP = IRECRD

         DO IX = 1, MXUSBX
            NX = IX
            READ (1,ERR=310,END=1000)
     &                   MX, TITUSX(NX), ITUSBX(NX), IDUSBX(NX),
     &                   NR1USX(NX), NR2USX(NX), AUSBDX (NX),
     &                   LWUSBX(NX), LFUSBX(NX), LLNUSX (NX),
     &                   EBXLOW(NX), EBXHGH(NX),
     &                   NEBXBN(NX), DEBXBN(NX),
     &                   ABXLOW(NX), ABXHGH(NX),
     &                   NABXBN(NX), DABXBN(NX)
            IRECRD = IRECRD + 1
            IF ( LLNUSX (NX) ) THEN
               READ (1) IGMUSX(NX)
               IRECRD = IRECRD + 1
            ELSE
               IGMUSX (NX) = 0
            END IF
            GO TO 350

 310        CONTINUE
*           BACKSPACE 1
            REWIND 1
            DO IR = 1, IRECRD
               READ (1)
            END DO
            READ (1,ERR=320,END=1000)
     &                   MX, TITUSX(NX), ITUSBX(NX), IDUSBX(NX),
     &                   NR1USX(NX), NR2USX(NX), AUSBDX (NX),
     &                   LWUSBX(NX), LFUSBX(NX),
     &                   EBXLOW(NX), EBXHGH(NX),
     &                   NEBXBN(NX), DEBXBN(NX),
     &                   ABXLOW(NX), ABXHGH(NX),
     &                   NABXBN(NX), DABXBN(NX)
            IRECRD = IRECRD + 1
            IGMUSX (NX) = 0
            GO TO 350

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

 330        CONTINUE
            STOP 'UNKNOWN RECORD TYPE'

 350        CONTINUE

            READ (1)
            IRECRD = IRECRD + 1

            MAXSIZE = MAXSIZE + ( 2 + NABXBN(NX) ) *
     &                   ( NEBXBN(NX) + IGMUSX(NX) )
            MNEBXBN = MAX ( MNEBXBN, NEBXBN(NX) )
            MIGMUSX = MAX ( MIGMUSX, IGMUSX(NX) )
            MNABXBN = MAX ( MNABXBN, NABXBN(NX) )

         END DO
 390     CONTINUE

         ALLOCATE( GMSTOR(MAXSIZE) )
         ALLOCATE( GBSTOR(MAXSIZE) )
         GDSTOR => GMSTOR

         ALLOCATE( ENGMAX(MIGMUSX + 1) )
         ALLOCATE( EPGMAX(MNEBXBN + MIGMUSX + 1) )
         ALLOCATE( OMGMAX(0:MNABXBN) )

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

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

      KLAST  = 0
      KMAX   = 0
*  +-------------------------------------------------------------------*
*  |                            Read from binary file the results!!!
      DO 6600 IX = 1, MXUSBX
         NX = IX
         READ (1,ERR=6400,END=1000)
     &                MX, TITUSX(NX), ITUSBX(NX), IDUSBX(NX),
     &                NR1USX(NX), NR2USX(NX), AUSBDX (NX),
     &                LWUSBX(NX), LFUSBX(NX), LLNUSX (NX),
     &                EBXLOW(NX), EBXHGH(NX),
     &                NEBXBN(NX), DEBXBN(NX),
     &                ABXLOW(NX), ABXHGH(NX),
     &                NABXBN(NX), DABXBN(NX)
         IRECRD = IRECRD + 1
         WRITE (*,*) NX,' ',TITUSX(NX)
         IF ( LLNUSX (NX) ) THEN
            READ (1) IGMUSX(NX), (ENGMAX (IG), IG=1,IGMUSX(NX)+1)
            IRECRD = IRECRD + 1
         ELSE
            IGMUSX (NX) = 0
         END IF
         GO TO 6500
 6400    CONTINUE
*        BACKSPACE 1
         REWIND 1
         DO IR = 1, IRECRD
            READ (1)
         END DO
         READ (1,ERR=6450,END=1000)
     &                MX, TITUSX(NX), ITUSBX(NX), IDUSBX(NX),
     &                NR1USX(NX), NR2USX(NX), AUSBDX (NX),
     &                LWUSBX(NX), LFUSBX(NX),
     &                EBXLOW(NX), EBXHGH(NX),
     &                NEBXBN(NX), DEBXBN(NX),
     &                ABXLOW(NX), ABXHGH(NX),
     &                NABXBN(NX), DABXBN(NX)
         IRECRD = IRECRD + 1
         LLNUSX (NX) = .FALSE.
         IGMUSX (NX) = 0
         GO TO 6500
 6450    CONTINUE
*        BACKSPACE 1
         REWIND 1
         DO IR = 1, IRECRD
            READ (1)
         END DO
         READ (1,ERR=6460,END=1000)CHSTAT,ANNX
         IRECRD = IRECRD + 1
         IF (CHSTAT .EQ. 'STATISTICS') GO TO 2000
 6460    CONTINUE
         STOP 'UNKNOWN RECORD TYPE'
 6500    CONTINUE
         JX (IX) = MX
*        JX (IX) = NX
         KBUSBX (NX) = KLAST
         K0 = KLAST + 1
         K1 = NEBXBN (NX) * NABXBN (NX) + K0 - 1
         K2 = K1 + IGMUSX (NX) * NABXBN (NX)
         KLAST = K2
         KMAX = MAX ( KMAX, K2 )
         READ (1) (GMSTOR(J), J = K0, K2)
         IRECRD = IRECRD + 1
 6600 CONTINUE
*  |
*  +-------------------------------------------------------------------*
 1000 CONTINUE
      WRITE (*,*) ' END OF FILE'
      IX = IX-1
      GO TO 2300
 2000 CONTINUE
      IX = IX-1
      INDX = INDLST ( FILE, '.' )
* Start_VAX_seq
*     IF ( INDX .GT. 0 ) THEN
*        OPEN (UNIT=3,FILE=FILE(1:INDX-1)//'_SUM.LIS',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
*     ELSE
*        OPEN (UNIT=3,FILE=FILE(1:LQ)//'_SUM.LIS',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
*     END IF
* End_VAX_seq
* Start_UNIX_seq
      IF ( INDX .GT. 0 ) THEN
         OPEN (UNIT=3,FILE=FILE(1:INDX-1)//'_sum.lis',STATUS='UNKNOWN')
      ELSE
         OPEN (UNIT=3,FILE=FILE(1:LQ)//'_sum.lis',STATUS='UNKNOWN')
      END IF
* End_UNIX_seq
      WRITE (3,*)' **** ',FILE(1:LQ),' ****'
      WRITE (3,*)
      WRITE (3,*)' ',RUNTIT
      WRITE (3,*)
      WRITE (CHCTOT,'(F15.0)') DBLE (NCTOT) + DBLE (MCTOT) * 1.D+09
      WRITE (3,'(A,2X,A14)')'  Total primaries run:',CHCTOT
      WRITE (3,*)'  Total weight of the primaries run:', WCTOT
      WRITE (3,*)
      DOMEGA = 0.D+00
      GFACT  = 0.E+00
*  +-------------------------------------------------------------------*
*  |
      DO 6850 JJ=1,IX
         NX=JJ
         MX=JX(JJ)
         KFUSBX (NX) = KLAST
         KLAST = KLAST + 2 * ( NEBXBN (NX) + IGMUSX (NX) )
         READ (1) TOTTOT(NX),TOTERR(NX)
         READ (1) NEBXBN(NX),IGMUSX(NX),EBXLOW(NX),(EPGMAX(IG),
     &            IG=1,NEBXBN(NX)+IGMUSX(NX)+1)
         KK0 = KFUSBX (NX) + 1
         KK2 = KFUSBX (NX) + NEBXBN (NX) + IGMUSX (NX)
         READ (1) ( GMSTOR(J), J = KK0, KK2 )
         READ (1) ( GBSTOR(J), J = KK0, KK2 )
         KK0 = KK2 + 1
         KK2 = KK2 + NEBXBN (NX) + IGMUSX (NX)
         READ (1) ( GMSTOR(J), J = KK0, KK2 )
         READ (1) ( GBSTOR(J), J = KK0, KK2 )
         K0 = KBUSBX (NX) + 1
         K1 = NEBXBN (NX) * NABXBN (NX) + K0 - 1
         K2 = K1 + IGMUSX (NX) * NABXBN (NX)
         READ (1) ( GBSTOR(J), J = K0, K2)
         KK0 = KFUSBX (NX) + 1
         KK1 = KFUSBX (NX) + NEBXBN (NX)
         KK2 = KFUSBX (NX) + NEBXBN (NX) + IGMUSX (NX)
         WRITE (3,*)
         WRITE (3,*) ' Detector n: ', NX,'(',MX,')  ',TITUSX (NX)
         WRITE (3,*) '    (Area:         ',AUSBDX(NX),' cmq,'
         WRITE (3,*) '     distr. scored:',IDUSBX(NX),'   ,'
         WRITE (3,*) '     from reg.',NR1USX(NX),' to ',NR2USX(NX),','
         IF ( .NOT. LWUSBX(NX) ) THEN
            WRITE (3,*) '     one way scoring,'
         ELSE
            WRITE (3,*) '     two ways scoring,'
         END IF
         IF ( .NOT. LFUSBX(NX) ) THEN
            WRITE (3,*) '     current scoring)'
         ELSE
            WRITE (3,*) '     fluence scoring scoring)'
         END IF
         WRITE (3,*)
         WRITE (3,10000)'    Tot. resp. (Part/cmq/pr)',
     &                   TOTTOT(NX),'  +/-',
     &                   TOTERR(NX)*100.0,' %'
         WRITE (3,10000)'    ( -->      (Part/pr)    ',
     &                   TOTTOT(NX)*AUSBDX(NX),'  +/-',
     &                   TOTERR(NX)*100.0,' % )'
         WRITE (3,*)
         WRITE (3,*)
     &         ' **** Different. Fluxes as a function of energy ****',
     &         ' ****      (integrated over solid angle)        ****'
         WRITE (3,*)
         WRITE (3,*)
         WRITE (3,*) ' Energy boundaries (GeV):'
         WRITE (3,*)
         WRITE (3,10002) ( EPGMAX (J), J = 1, NEBXBN (NX) )
         WRITE (3,10003) ' Lowest boundary   (GeV):',
     &                     EPGMAX (NEBXBN(NX)+1)
         WRITE (3,*)
         WRITE (3,*) ' Flux (Part/GeV/cmq/pr):'
         WRITE (3,*)
         WRITE (3,10001) ( GDSTOR(J),' +/- ',
     &                     GBSTOR(J)*100.0,' % ',
     &                     J = KK0, KK1, 1 )
         IF ( LLNUSX (NX) ) THEN
            WRITE (3,*)
            WRITE (3,*) ' Energy boundaries (GeV):'
            WRITE (3,*)
            WRITE (3,10002) ( EPGMAX (J), J = NEBXBN (NX) + 1,
     &                        NEBXBN (NX) + IGMUSX (NX) )
            WRITE (3,10003) ' Lowest boundary   (GeV):',
     &                        EPGMAX (NEBXBN(NX)+IGMUSX(NX)+1)
            WRITE (3,*)
            WRITE (3,*) ' Flux (Part/GeV/cmq/pr):'
            WRITE (3,*)
            WRITE (3,10001) ( GDSTOR(J),' +/- ',
     &                        GBSTOR(J)*100.0,' % ',
     &                        J = KK1+1, KK2 )
         END IF
         KK0 = KFUSBX (NX) + NEBXBN (NX) + IGMUSX (NX) + 1
         KK1 = KFUSBX (NX) + 2 * NEBXBN (NX) + IGMUSX (NX)
         KK2 = KFUSBX (NX) + 2 * ( NEBXBN (NX) + IGMUSX (NX) )
         WRITE (3,*)
         WRITE (3,*)
         WRITE (3,*)
     &         ' **** Cumulative Fluxes as a function of energy ****',
     &         ' ****      (integrated over solid angle)        ****'
         WRITE (3,*)
         WRITE (3,*)
         WRITE (3,*) ' Energy boundaries (GeV):'
         WRITE (3,*)
         WRITE (3,10002) ( EPGMAX (J), J = 1, NEBXBN (NX) )
         WRITE (3,10003) ' Lowest boundary   (GeV):',
     &                     EPGMAX (NEBXBN(NX)+1)
         WRITE (3,*)
         WRITE (3,*) ' Cumul. Flux (Part/cmq/pr):'
         WRITE (3,*)
         WRITE (3,10001) ( GDSTOR(J),' +/- ',
     &                     GBSTOR(J)*100.0,' % ',
     &                     J = KK0, KK1, 1 )
         IF ( LLNUSX (NX) ) THEN
            WRITE (3,*)
            WRITE (3,*) ' Energy boundaries (GeV):'
            WRITE (3,*)
            WRITE (3,10002) ( EPGMAX (J), J = NEBXBN (NX) + 1,
     &                        NEBXBN (NX) + IGMUSX (NX) )
            WRITE (3,10003) ' Lowest boundary   (GeV):',
     &                        EPGMAX (NEBXBN(NX)+IGMUSX(NX)+1)
            WRITE (3,*)
            WRITE (3,*) ' Cumul. Flux (Part/cmq/pr):'
            WRITE (3,*)
            WRITE (3,10001) ( GDSTOR(J),' +/- ',
     &                        GBSTOR(J)*100.0,' % ' ,
     &                        J = KK1+1, KK2 )
         END IF
         IF ( NABXBN (NX) .LE. 1 .AND. ABXLOW (NX) .LT. 1.E-05
     &        .AND. ABXHGH (NX) .GT. 6.28 ) GO TO 6850
*  |  +----------------------------------------------------------------*
*  |  |
         IF ( ABS ( ITUSBX (NX) ) .GE. 2 ) THEN
            OMEGA0 = ABXLOW (NX)
            GFACT  = DABXBN (NX)
            OMGMAX (0) = 0.0
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |
         ELSE
            OMEGA0 = ABXLOW (NX)
            DOMEGA = DBLE (DABXBN (NX))
            OMGMAX (0) = OMEGA0
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
*  |  +----------------------------------------------------------------*
*  |  |
         IF ( ABS (ITUSBX (NX)) .GE. 2 ) THEN
            OMEOLD = 0.0
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |
         ELSE
            OMEOLD = OMEGA0
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
*  |  +----------------------------------------------------------------*
*  |  |
         DO IA = 1, NABXBN (NX)
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            IF ( ABS (ITUSBX (NX)) .GE. 2 ) THEN
               OMEGAA = OMEGA0 * GFACT**(IA - 1)
               DOMEGA = DBLE (OMEGAA - OMEOLD)
*  |  |  |
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            ELSE
               OMEGAA = OMEGA0 + DOMEGA * FLOAT (IA)
            END IF
*  |  |  |
*  |  |  +-------------------------------------------------------------*
            OMEOLD = OMEGAA
            OMGMAX (IA) = OMEGAA
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
         WRITE (3,*)
         WRITE (3,*)
     &         ' **** Double diff. Fluxes as a function of energy ****'
         WRITE (3,*)
         WRITE (3,*)
         WRITE (3,10003) ' Solid angle minimum value (sr):',
     &                     OMGMAX (0)
         WRITE (3,*)
         WRITE (3,*) ' Solid angle upper boundaries (sr):'
         WRITE (3,*)
         WRITE (3,10002) ( OMGMAX (J), J = 1, NABXBN (NX) )
         WRITE (3,*)
         WRITE (3,10003) ' Angular minimum value (deg.):',
     &                     SNGL ( ACOS ( MAX ( 1.D+00 - OMGMAX (0)
     &                          / TWOPIP, - 1.D+00 ) )
     &                          * 180.D+00 / PIPIPI )
         WRITE (3,*)
         WRITE (3,*) ' Angular upper boundaries (deg.):'
         WRITE (3,*)
         WRITE (3,10002) ( SNGL ( ACOS ( MAX ( 1.D+00 - OMGMAX (J)
     &                          / TWOPIP, - 1.D+00 ) )
     &                          * 180.D+00 / PIPIPI ),
     &                          J = 1, NABXBN (NX) )
         WRITE (3,*)
         DO 6750 JE = 1, NEBXBN (NX)
            WRITE (3,*)
            WRITE (3,10004) ' Energy interval (GeV):', EPGMAX (JE),
     &                        EPGMAX (JE+1)
            IE = NEBXBN (NX) - JE + 1
            WRITE (3,*)
            WRITE (3,*) ' Flux (Part/sr/GeV/cmq/pr):'
            WRITE (3,*)
            WRITE (3,10001) ( GDSTOR(IE+NEBXBN(NX)*(IA-1)+K0-1),
     &                        ' +/- ',
     &                        GBSTOR(IE+NEBXBN(NX)*(IA-1)+K0-1)
     &                             *100.0,' % ',
     &                        IA = 1, NABXBN (NX) )
            WRITE (3,*)
            WRITE (3,*) ' Flux (Part/deg/GeV/cmq/pr):'
            WRITE (3,*)
            WRITE (3,10001) ( GDSTOR(IE+NEBXBN(NX)*(IA-1)+K0-1)
     &                         *(OMGMAX(IA)-OMGMAX(IA-1))
     &                    / ( ACOS ( MAX ( 1.D+00 - OMGMAX(IA) / TWOPIP,
     &                        - 1.D+00 ) ) * 180.D+00 / PIPIPI -
     &                      ACOS ( MAX ( 1.D+00 - OMGMAX(IA-1) / TWOPIP,
     &                        - 1.D+00 ) ) * 180.D+00 / PIPIPI ),
     &                        ' +/- ',
     &                        GBSTOR(IE+NEBXBN(NX)*(IA-1)+K0-1)
     &                             *100.0,' % ',
     &                        IA = 1, NABXBN (NX) )
 6750    CONTINUE
         IF ( LLNUSX (NX) ) THEN
            DO 6780 IE = NEBXBN(NX) + 1, NEBXBN (NX) + IGMUSX (NX)
               WRITE (3,*)
               WRITE (3,10004) ' Energy interval (GeV):', EPGMAX (IE),
     &                           EPGMAX (IE+1)
               IG  = IE - NEBXBN (NX)
               WRITE (3,*)
               WRITE (3,*) ' Flux (Part/sr/GeV/cmq/pr):'
               WRITE (3,*)
               WRITE (3,10001) ( GDSTOR(IG+IGMUSX(NX)*(IA-1)+K1),
     &                           ' +/- ',
     &                           GBSTOR(IG+IGMUSX(NX)*(IA-1)+K1)
     &                                 * 100.0,' % ',
     &                           IA = 1, NABXBN (NX) )
               WRITE (3,*)
               WRITE (3,*) ' Flux (Part/deg/GeV/cmq/pr):'
               WRITE (3,*)
               WRITE (3,10001) ( GDSTOR(IG+IGMUSX(NX)*(IA-1)+K1)
     &                         *(OMGMAX(IA)-OMGMAX(IA-1))
     &                    / ( ACOS ( MAX ( 1.D+00 - OMGMAX(IA) / TWOPIP,
     &                        - 1.D+00 ) ) * 180.D+00 / PIPIPI -
     &                      ACOS ( MAX ( 1.D+00 - OMGMAX(IA-1) / TWOPIP,
     &                        - 1.D+00 ) ) * 180.D+00 / PIPIPI ),
     &                        ' +/- ',
     &                        GBSTOR(IG+IGMUSX(NX)*(IA-1)+K1)
     &                              * 100.0,' % ',
     &                        IA = 1, NABXBN (NX) )
 6780       CONTINUE
         END IF
 6850 CONTINUE
*  |
*  +-------------------------------------------------------------------*

      DEALLOCATE(GMSTOR)
      DEALLOCATE(GBSTOR)
      GDSTOR => NULL()
      DEALLOCATE(ENGMAX)
      DEALLOCATE(EPGMAX)
      DEALLOCATE(OMGMAX)

 2300 CONTINUE
      CLOSE (UNIT=1)
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))
10004 FORMAT ((1X,A,1X,1P,G14.7,1X,G14.7))
      STOP
      END

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

*
      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
