*$ CREATE USXSUW.FOR
*COPY USXSUW
*
*=== Usxsuw ===========================================================*
*
      PROGRAM USXSUW
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1990-2008      by    Alfredo Ferrari & Paola Sala  *
*     All Rights Reserved.                                             *
*                                                                      *
*----------------------------------------------------------------------*
*
      CHARACTER FILE*80,RUNTIT*80,RUNTIM*32,CHSTAT*10,CHCTOT*15
*
*=== usrbdx ==========================================================*
*
*---------------------------------------------------------------------*
*     Module USRBDX:                                                  *
*     A. Ferrari : user defined boundary crossing scoring             *
*          Last change A. Ferrari 16-may-1990                         *
*                                                                     *
*                                                                     *
*     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 = 700 )
      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 ( MXDUMM = 300000 )
      PARAMETER ( NMXGRP = 2000 )
      PARAMETER ( NMXEBN = 2200 )
      PARAMETER ( PIPIPI = 3.141592653589793D+00 )
      PARAMETER ( TWOPIP = 2.D+00 * PIPIPI )
      DIMENSION GMSTOR (MXDUMM), GDSTOR (MXDUMM), GBSTOR (MXDUMM),
     &          TOTTOT (MXUSBX), TOTERR (MXUSBX),
     &          KFUSBX (MXUSBX), JX (MXUSBX), ENGMAX (NMXGRP+1),
     &          EPGMAX (NMXEBN), ACCUMG (NMXEBN+NMXGRP),
     &          ACUMCU (NMXEBN+NMXGRP), OMGMAX (0:NMXEBN)
      DOUBLE PRECISION GDSTOR, ACCUMP, GBSTOR, ACCUMG, COSTHE, AFACT,
     &                 TOTTOT, TOTERR, DSCORE, DSCOR2, DELTAE, DOMEGA,
     &                 ACUMCU, ACCUMT, ACCUME, WCTOT , GBHELP, TOTHLP
*
*     Note that all addresses are for zero index (Kbusbx,Kfusbx)
*
      DO 5555 I = 1, MXDUMM
         GBSTOR (I) = 0.D+00
         GDSTOR (I) = 0.D+00
 5555 CONTINUE
      DO 4444 I = 1, MXUSBX
         TOTTOT (I) = 0.D+00
         TOTERR (I) = 0.D+00
 4444 CONTINUE
      NCTOT  = 0
      MCTOT  = 0
      WCTOT  = 0.D+00
      MBATCH = 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)
      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  = 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
      KLAST  = 0
      KMAX   = 0
*  +-------------------------------------------------------------------*
*  |                            Read from binary file the results!!!
      DO 6600 IX = 1,1000
         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
         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) CHSTAT,IJX
         IRECRD = IRECRD + 1
         IF ( CHSTAT .EQ. 'STATISTICS') GO TO 1000
 6460    CONTINUE
         STOP 'UNKNOWN RECORD TYPE'
 6500    CONTINUE
         JX (IX) = NX
*        JX (IX) = MX
***************************************************************
*********** Note Kbusbx is used for 0 address *****************
***************************************************************
         KBUSBX (NX) = KLAST
         K0 = KLAST + 1
         K1 = NEBXBN (NX) * NABXBN (NX) + K0 - 1
         K2 = K1 + IGMUSX (NX) * NABXBN (NX)
         KLAST = K2
         READ (1) (GMSTOR(J), J = K0, K2)
         IRECRD = IRECRD +1
         DO J = K0, K2
            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 )
 6600 CONTINUE
*  |
*  +-------------------------------------------------------------------*
 1000 CONTINUE
      IX = IX-1
      CLOSE (UNIT=1)
      IF ( MBATCH .EQ. 1 ) THEN
         DO IIX = 1, IX
            NX = IIX
            MX = JX (IIX)
            KFUSBX (NX) = KLAST
            KLAST = KLAST + 2 * ( NEBXBN (NX) + IGMUSX (NX) )
         END DO
      END IF
      KLAST = 0
*  +-------------------------------------------------------------------*
*  |
      DO 7000 IIX = 1, IX
         NX = IIX
         MX = JX (IIX)
         ACCUMP = 0.D+00
         DO IG = 1, NEBXBN (NX) + IGMUSX (NX)
            ACCUMG (IG) = 0.D+00
            ACUMCU (IG) = 0.D+00
         END DO
         AFACT  = DBLE (NCASE) + DBLE (MCASE) * 1.D+09
         K0     = KLAST + 1
         K1     = NEBXBN (NX) * NABXBN (NX) + K0 - 1
         K2     = K1 + IGMUSX (NX) * NABXBN (NX)
         KLAST  = K2
         ENERG0 = EBXLOW (NX)
*  |  +----------------------------------------------------------------*
*  |  |
         IF ( ITUSBX (NX) .LT. 0 ) THEN
            EFACT  = DEBXBN (NX)
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |
         ELSE
            DELTAE = DBLE (DEBXBN (NX))
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         ENEOLD = ENERG0
*  |  +----------------------------------------------------------------*
*  |  |
         IF ( ABS ( ITUSBX (NX) ) .GE. 2 ) THEN
            OMEGA0 = ABXLOW (NX)
            GFACT  = DABXBN (NX)
            XGMIN  = 0.3 * LOG10 (OMEGA0)
            XGMAX  = LOG10 ( ABXHGH (NX) )
            XGMIN  = ACOS(-10.0**XGMIN/TWOPIP+1.0)*180.0/PIPIPI
            XGMAX  = ACOS(MAX(-10.0**XGMAX/TWOPIP+1.0,-1.0))
     &             * 180.0/PIPIPI
            XGMIN  = 0.0
            OMGMAX (0) = 0.0
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |
         ELSE
            OMEGA0 = ABXLOW (NX)
            DOMEGA = DBLE (DABXBN (NX))
            XGMIN = OMEGA0
            XGMAX = ABXHGH (NX)
            XGMIN  = ACOS(-XGMIN/TWOPIP+1.0)*180.0/PIPIPI
            XGMAX  = ACOS(-XGMAX/TWOPIP+1.0)*180.0/PIPIPI
            OMGMAX (0) = OMEGA0
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         YGMAX = -1.E+30
         YGMIN = 1.E+30
*  |  +----------------------------------------------------------------*
*  |  |
         DO 6540 IE = 1, NEBXBN (NX) + IGMUSX (NX)
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            IF ( IE .LE. NEBXBN (NX) ) THEN
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |
               IF ( ITUSBX (NX) .LT. 0 ) THEN
                  ENERGE = ENERG0 * EFACT**IE
                  DELTAE = DBLE (ENERGE - ENEOLD)
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |
               ELSE
                  ENERGE = EBXLOW (NX) + DEBXBN (NX) * FLOAT (IE)
               END IF
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
*  |  |  |
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            ELSE
               IF ( IE-NEBXBN (NX) .EQ. 1) ENEOLD = ENGMAX (1)
               ENERGE = ENGMAX (IE-NEBXBN(NX)+1)
               DELTAE = DBLE(ENEOLD - ENERGE)
            END IF
*  |  |  |
*  |  |  +-------------------------------------------------------------*
            ENEOLD = ENERGE

            IF (DELTAE.LE. 0.D+00)WRITE(*,*)' **** DELTAE',IE,NX

*  |  |  +-------------------------------------------------------------*
*  |  |  |
            IF ( ABS (ITUSBX (NX)) .GE. 2 ) THEN
               OMEOLD = 0.0
*  |  |  |
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            ELSE
               OMEOLD = OMEGA0
            END IF
*  |  |  |
*  |  |  +-------------------------------------------------------------*
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            DO 6520 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

               IF (DOMEGA.LE. 0.D+00)WRITE(*,*)' **** DOMEGA',IE,NX

*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |
               IF ( IE .LE. NEBXBN (NX) ) THEN
                  IEA = IE + NEBXBN (NX) * ( IA - 1 ) + K0 - 1
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |
               ELSE
                  IG  = IE - NEBXBN (NX)
                  IEA = IG + IGMUSX (NX) * ( IA - 1 ) + K1
               END IF
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
               DSCORE = DBLE ( GMSTOR (IEA) ) * ( DELTAE * DOMEGA )
               DSCOR2 = DBLE ( GMSTOR (IEA) ) * DOMEGA
               ACCUMP = ACCUMP + DSCORE
               ACUMCU (IE) = ACUMCU (IE) + DSCORE
               ACCUMG (IE) = ACCUMG (IE) + DSCOR2
               YGMAX = MAX ( YGMAX, GMSTOR (IEA) )
               YGMIN = MIN ( YGMIN, GMSTOR (IEA) )
6520        CONTINUE
*  |  |  |
*  |  |  +-------------------------------------------------------------*
6540     CONTINUE
*  |  |
*  |  +----------------------------------------------------------------*
         TOTTOT (NX) = TOTTOT (NX) + ACCUMP * AFACT
         TOTERR (NX) = TOTERR (NX) + ACCUMP * ACCUMP * AFACT

         ACCUMT=0.D+00
         DO IE = 1,NEBXBN(NX)+IGMUSX(NX)
            ACCUMT = ACCUMT + ACUMCU (IE)
         END DO
         IF ( ABS ( ACCUMT - ACCUMP ) .GT.
     &        1.D-13 * ACCUMP ) WRITE (*,*)' *** ERROR 0', NX,
     &        ACCUMT, ACCUMP
         ACCUME=0.D+00
         DO IE = 1,NEBXBN(NX)
            ACCUME = ACCUME + ACUMCU (IE)
         END DO

*  |  +----------------------------------------------------------------*
*  |  |
         DO IE = NEBXBN (NX) - 1, 1, -1
            ACUMCU (IE) = ACUMCU (IE) + ACUMCU (IE+1)
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
         ACUMCU (NEBXBN(NX)+1) = ACUMCU (NEBXBN(NX)+1) + ACUMCU (1)
*  |  +----------------------------------------------------------------*
*  |  |
         DO IE = NEBXBN (NX) + 2, NEBXBN (NX) + IGMUSX (NX)
            ACUMCU (IE) = ACUMCU (IE) + ACUMCU (IE-1)
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
         IF ( LLNUSX (NX) .AND.
     &      ABS ( ACUMCU ( NEBXBN (NX) + IGMUSX (NX) ) - ACCUMP ) .GT.
     &      1.D-13 * ACCUMP ) WRITE (*,*)' *** ERROR', NX,
     &      ACUMCU ( NEBXBN (NX) + IGMUSX (NX) ), ACCUMP
*  |  +----------------------------------------------------------------*
*  |  |
         DO IE = 1, NEBXBN (NX) + IGMUSX (NX)
            KK0 = KFUSBX (NX) + IE
            GDSTOR (KK0) = GDSTOR (KK0) + ACCUMG (IE) * AFACT
            GBSTOR (KK0) = GBSTOR (KK0) + ACCUMG (IE) * ACCUMG (IE) *
     &                     AFACT
            KK0 = KK0 + NEBXBN (NX) + IGMUSX (NX)
            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
         NX = IIX
         MX = JX (IIX)
         TOTTOT (NX) = TOTTOT (NX)
     &               / ( DBLE (NCTOT) + DBLE (MCTOT) * 1.D+09 )
         TOTERR (NX) = TOTERR (NX)
     &               / ( DBLE (NCTOT) + DBLE (MCTOT) * 1.D+09 )
      END DO
      IF ( MBATCH - 1 .GT. 0 ) THEN
         DO J = 1, KMAX
            GBHELP     = SQRT ( GBSTOR (J) )
            GBSTOR (J) = ( ( GBHELP - GDSTOR (J) )
     &                   * ( GBHELP + GDSTOR (J) ) )
     &                 / DBLE (MBATCH-1)
            IF ( GBSTOR (J) .LT. 1.D-12 * GDSTOR (J)**2 )
     &           GBSTOR (J) = 0.D+00
            IF ( GDSTOR (J) .GT. 0.D+00 )
     &          GBSTOR (J) = MIN ( SQRT ( GBSTOR (J) ) / GDSTOR (J),
     &                             0.99D+00 )
         END DO
         DO IIX = 1, IX
            NX = IIX
            MX = JX (IIX)
            TOTHLP = SQRT ( TOTERR (NX) )
            TOTERR (NX) = ( ( TOTHLP - TOTTOT (NX) )
     &                    * ( TOTHLP + TOTTOT (NX) ) )
     &                  / DBLE (MBATCH-1)
            IF ( TOTERR (NX) .LT. 1.D-12 * TOTTOT (NX)**2 )
     &           TOTERR (NX) = 0.D+00
            IF ( TOTTOT (NX) .GT. 0.D+00 ) THEN
               TOTERR (NX) = MIN ( SQRT ( TOTERR (NX) ) / TOTTOT (NX),
     &                             0.99D+00 )
            END IF
         END DO
      ELSE
         DO J =1, KMAX
            IF ( GDSTOR(J) .GT. 0.D+00 ) THEN
               GBSTOR (J) = 0.99D+00
            END IF
         END DO
         DO IIX = 1, IX
            NX = IIX
            MX = JX (IIX)
            IF ( TOTTOT (NX) .GT. 0.D+00 ) THEN
               TOTERR (NX) = 0.99D+00
            END IF
         END DO
      END IF
      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
      OPEN (UNIT=1,FILE=FILE,STATUS='UNKNOWN',FORM='UNFORMATTED')
* End_UNIX_seq
      WRITE (1) RUNTIT,RUNTIM,SNGL(WCTOT),NCTOT,MCTOT,MBATCH
      KLAST = 0
*  +-------------------------------------------------------------------*
*  |
      DO 6800 JJ=1,IX
         NX=JJ
         MX=JX(JJ)
         WRITE (1)    NX, 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)
         IF ( LLNUSX (NX) ) THEN
            WRITE (1) IGMUSX(NX), (ENGMAX (IG), IG=1,IGMUSX(NX)+1)
         END IF
         K0 = KLAST + 1
         K1 = NEBXBN (NX) * NABXBN (NX) + K0 - 1
         K2 = K1 + IGMUSX (NX) * NABXBN (NX)
         KLAST = K2
         WRITE (1) (SNGL(GDSTOR(J)), J = K0, K2)
6800  CONTINUE
*  |
*  +-------------------------------------------------------------------*
*  Now write the statistical data, first a record to flag them
      WRITE (1) 'STATISTICS',IX
*  Then the errors concerning the binning data
      KLAST = 0
      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
* Start_VAX_seq
*     IF ( INDX .GT. 0 ) THEN
*        OPEN (UNIT=4,FILE=FILE(1:INDX-1)//'_TAB.LIS',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
*     ELSE
*        OPEN (UNIT=4,FILE=FILE(1:LQ)//'_TAB.LIS',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
*     END IF
* End_VAX_seq
* Start_UNIX_seq
      IF ( INDX .GT. 0 ) THEN
         OPEN (UNIT=4,FILE=FILE(1:INDX-1)//'_tab.lis',STATUS='UNKNOWN')
      ELSE
         OPEN (UNIT=4,FILE=FILE(1:LQ)//'_tab.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:',SNGL(WCTOT)
      WRITE (3,*)
*  +-------------------------------------------------------------------*
*  |
      DO 6850 JJ = 1, IX
         NX     = JJ
         MX     = JX (JJ)
*        WRITE (2,*)'  N: ',MX,'  ',TITUSX(NX)
         WRITE (1) SNGL(TOTTOT(NX)),SNGL(TOTERR(NX))
*        WRITE (2,*)'    Tot. response', SNGL(TOTTOT(NX)),
*    &              '  +/-',SNGL(TOTERR(NX))*100.0,' %'
         ENERG0 = EBXLOW (NX)
*  |  +----------------------------------------------------------------*
*  |  |
         IF ( ITUSBX (NX) .LT. 0 ) THEN
            EFACT  = DEBXBN (NX)
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |
         ELSE
            DELTAE = DBLE (DEBXBN (NX))
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         ENEOLD = ENERG0
*  |  +----------------------------------------------------------------*
*  |  |
         DO IE = 1, NEBXBN (NX)
            IG = NEBXBN (NX) - IE + 1
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            IF ( ITUSBX (NX) .LT. 0 ) THEN
               ENERGE = ENERG0 * EFACT**IE
               DELTAE = DBLE (ENERGE - ENEOLD)
*  |  |  |
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            ELSE
               ENERGE = EBXLOW (NX) + DEBXBN (NX) * FLOAT (IE)
            END IF
*  |  |  |
*  |  |  +-------------------------------------------------------------*
            EPGMAX (IG) = ENERGE
            ENEOLD = ENERGE
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
*  |  +----------------------------------------------------------------*
*  |  |
         DO IG = NEBXBN (NX) + 1, NEBXBN (NX) + IGMUSX (NX)
            EPGMAX (IG) = ENGMAX (IG-NEBXBN(NX))
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
*  |  +----------------------------------------------------------------*
*  |  |
         IF (LLNUSX (NX)) THEN
            EPGMAX (NEBXBN(NX)+IGMUSX(NX)+1) = ENGMAX (IGMUSX(NX)+1)
            IF (IDUSBX (NX) .NE. 8) EPGMAX (NEBXBN (NX)+1) = EBXLOW (NX)
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |
         ELSE
            EPGMAX (NEBXBN(NX)+1) = EBXLOW (NX)
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
         WRITE (1) NEBXBN(NX),IGMUSX(NX),EBXLOW(NX),(EPGMAX(IG),
     &             IG=1,NEBXBN(NX)+IGMUSX(NX)+1)
         KK0 = KFUSBX (NX) + 1
         KK1 = KFUSBX (NX) + NEBXBN (NX)
         KK2 = KFUSBX (NX) + NEBXBN (NX) + IGMUSX (NX)

*        ACCUMP=0.D+00
*        DO IE=1,NEBXBN(NX)
*           IG = NEBXBN (NX) - IE + 1
*           J=KFUSBX(NX)+IE
*           ACCUMP=ACCUMP+GDSTOR(J)*(EPGMAX(IG)-EPGMAX(IG+1))
*        END DO
*        WRITE(2,*) ACCUMP
*        DO IE=NEBXBN(NX)+1,NEBXBN(NX)+IGMUSX(NX)
*           IG = IE
*           J=KFUSBX(NX)+IE
*           ACCUMP=ACCUMP+GDSTOR(J)*(EPGMAX(IG)-EPGMAX(IG+1))
*        END DO
*        WRITE(2,*) ACCUMP

         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)',
     &                   SNGL(TOTTOT(NX)),'  +/-',
     &                   SNGL(TOTERR(NX))*100.0,' %'
         WRITE (3,10000)'    ( -->      (Part/pr)    ',
     &                   SNGL(TOTTOT(NX)*AUSBDX(NX)),'  +/-',
     &                   SNGL(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) ( SNGL (GDSTOR(J)),' +/- ',
     &                     SNGL (GBSTOR(J)*100.D+00),' % ',
     &                     J = KK1, KK0,-1 )
         IF ( LLNUSX (NX) ) THEN
            IF (IDUSBX (NX) .NE. 8)
     &         EPGMAX (NEBXBN (NX)+1) = ENGMAX (1)
            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) ( SNGL (GDSTOR(J)),' +/- ',
     &                        SNGL (GBSTOR(J)*100.D+00),' % ',
     &                        J = KK1+1, KK2 )
            IF (IDUSBX (NX) .NE. 8)
     &         EPGMAX (NEBXBN (NX)+1) = EBXLOW (NX)
         END IF
         WRITE (1) ( SNGL (GDSTOR(J)), J = KK1, KK0,-1 ),
     &             ( SNGL (GDSTOR(J)), J = KK1+1, KK2 )
         WRITE (1) ( SNGL (GBSTOR(J)), J = KK1, KK0,-1 ),
     &             ( SNGL (GBSTOR(J)), J = KK1+1, KK2 )
*****   here to write on "table" file
         IF (NX .GT. 1) THEN
*****       write 2 empty lines so gnuplot can use it with 'index #'
            WRITE (4,*)
            WRITE (4,*)
         END IF
         WRITE (4,*) '# Detector n: ', NX,'  ',TITUSX (NX),
     &                                '(integrated over solid angle)'
         WRITE (4,*) '# N. of energy intervals' ,NEBXBN(NX)+IGMUSX(NX)
         IF ( LLNUSX (NX) ) THEN
            IF (IDUSBX (NX) .NE. 8)
     &         EPGMAX (NEBXBN (NX)+1) = ENGMAX (1)
            KK = KK2 + 1
            DO 6001 J=NEBXBN(NX)+IGMUSX(NX), NEBXBN(NX)+1,-1
               KK = KK - 1
               WRITE (4,9501) EPGMAX (J+1), EPGMAX (J),
     &                   SNGL (GDSTOR(KK)), SNGL (GBSTOR(KK)*100.D+00)
 6001       CONTINUE
            IF (IDUSBX (NX) .NE. 8)
     &         EPGMAX (NEBXBN (NX)+1) = EBXLOW (NX)
         END IF
         KK=KK0-1
         DO 6002 J=NEBXBN(NX),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 = 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) ( SNGL (GDSTOR(J)),' +/- ',
     &                     SNGL (GBSTOR(J)*100.D+00),' % ',
     &                     J = KK1, KK0,-1 )
         IF ( LLNUSX (NX) ) THEN
            IF (IDUSBX (NX) .NE. 8)
     &         EPGMAX (NEBXBN (NX)+1) = ENGMAX (1)
            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) ( SNGL (GDSTOR(J)),' +/- ',
     &                        SNGL (GBSTOR(J)*100.D+00),' % ' ,
     &                        J = KK1+1, KK2 )
            IF (IDUSBX (NX) .NE. 8)
     &         EPGMAX (NEBXBN (NX)+1) = EBXLOW (NX)
         END IF
         WRITE (1) ( SNGL (GDSTOR(J)), J = KK1, KK0,-1 ),
     &             ( SNGL (GDSTOR(J)), J = KK1+1, KK2 )
         WRITE (1) ( SNGL (GBSTOR(J)), J = KK1, KK0,-1 ),
     &             ( SNGL (GBSTOR(J)), J = KK1+1, KK2 )
         K0 = KLAST + 1
         K1 = NEBXBN (NX) * NABXBN (NX) + K0 - 1
         K2 = K1 + IGMUSX (NX) * NABXBN (NX)
         KLAST = K2
         WRITE (1) (SNGL(GBSTOR(J)), J = K0, K2)
*  |
*  |     Commented out 06 Apr 2016
*  |     P. SC. : Consistancy, we always want the double differential
*  |              even in the simple 1-bin case
*  | 
c         IF ( NABXBN (NX) .LE. 1 .AND. ABXLOW (NX) .LT. 1.E-05
c     &        .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) ( SNGL (GDSTOR(IE+NEBXBN(NX)*(IA-1)+K0-1)),
     &                        ' +/- ',
     &                        SNGL (GBSTOR(IE+NEBXBN(NX)*(IA-1)+K0-1)
     &                             *100.D+00),' % ',
     &                        IA = 1, NABXBN (NX) )
            WRITE (3,*)
            WRITE (3,*) ' Flux (Part/deg/GeV/cmq/pr):'
            WRITE (3,*)
            WRITE (3,10001) ( SNGL (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 ) ),
     &                        ' +/- ',
     &                        SNGL (GBSTOR(IE+NEBXBN(NX)*(IA-1)+K0-1)
     &                             *100.D+00),' % ',
     &                        IA = 1, NABXBN (NX) )
 6750    CONTINUE
         IF ( LLNUSX (NX) ) THEN
            IF (IDUSBX (NX) .NE. 8)
     &         EPGMAX (NEBXBN (NX)+1) = ENGMAX (1)
            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) ( SNGL (GDSTOR(IG+IGMUSX(NX)*(IA-1)+K1)),
     &                           ' +/- ',
     &                           SNGL (GBSTOR(IG+IGMUSX(NX)*(IA-1)+K1)
     &                                *100.D+00),' % ',
     &                           IA = 1, NABXBN (NX) )
               WRITE (3,*)
               WRITE (3,*) ' Flux (Part/deg/GeV/cmq/pr):'
               WRITE (3,*)
               WRITE (3,10001) ( SNGL (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 ) ),
     &                        ' +/- ',
     &                        SNGL (GBSTOR(IG+IGMUSX(NX)*(IA-1)+K1)
     &                             *100.D+00),' % ',
     &                        IA = 1, NABXBN (NX) )
 6780       CONTINUE
            IF (IDUSBX (NX) .NE. 8)
     &         EPGMAX (NEBXBN (NX)+1) = EBXLOW (NX)
         END IF
*****   here to write on "table" file
         WRITE (4,*)
         WRITE (4,*)
         WRITE (4,*) '# double differential distributions'
         WRITE (4,*) '# number of solid angle intervals', NABXBN(NX)
*        WRITE (4,9502) ( OMGMAX (J-1), OMGMAX (J), J = 1, NABXBN (NX) )
         DO 7003 IA = 1, NABXBN (NX), 1
            WRITE (4,*) '# Block n: ',IA,' ',
     &                  OMGMAX (IA-1),' :',OMGMAX(IA)
            IF ( LLNUSX (NX) ) THEN
               IF (IDUSBX (NX) .NE. 8)
     &            EPGMAX (NEBXBN (NX)+1) = ENGMAX (1)
               DO 7001 JE = NEBXBN (NX) + IGMUSX (NX), NEBXBN(NX)+1,-1
                  IG  = JE - NEBXBN (NX)
                  EMAX = EPGMAX(JE)
                  WRITE (4,9504) EPGMAX (JE+1), EPGMAX(JE), !OMGMAX (IA-1),
     &                SNGL (GDSTOR(IG+IGMUSX(NX)*(IA-1)+K1)),
     &                SNGL (GBSTOR(IG+IGMUSX(NX)*(IA-1)+K1)*100.D+00)
 7001          CONTINUE
               IF (IDUSBX (NX) .NE. 8)
     &            EPGMAX (NEBXBN (NX)+1) = EBXLOW (NX)
           END IF
           DO 7002 JE = NEBXBN(NX), 1, -1
              IE = NEBXBN (NX) - JE + 1
              EMAX = EPGMAX(JE)
              WRITE (4,9504) EPGMAX (JE+1), EPGMAX(JE), !OMGMAX (IA-1),
     &            SNGL (GDSTOR(IE+NEBXBN(NX)*(IA-1)+K0-1)),
     &            SNGL (GBSTOR(IE+NEBXBN(NX)*(IA-1)+K0-1)*100.D+00)
 7002      CONTINUE
           WRITE (4,9504) EMAX, EMAX, 0.0, 0.0
*           WRITE (4,9504) EMAX, OMGMAX(IA-1), 0.0, 0.0
           WRITE (4,*)
 7003    CONTINUE

         IF ( LLNUSX (NX) ) THEN
            IF (IDUSBX (NX) .NE. 8)
     &         EPGMAX (NEBXBN (NX)+1) = ENGMAX (1)
            DO 7011 JE = NEBXBN (NX) + IGMUSX (NX), NEBXBN(NX)+1,-1
               EMAX = EPGMAX(JE)
               WRITE (4,9504) EPGMAX (JE+1), OMGMAX (NABXBN (NX)),
     &             0.0, 0.0
 7011       CONTINUE
            IF (IDUSBX (NX) .NE. 8)
     &         EPGMAX (NEBXBN (NX)+1) = EBXLOW (NX)
         END IF
         DO 7012 JE = NEBXBN(NX), 1, -1
            EMAX = EPGMAX(JE)
            WRITE (4,9504) EPGMAX (JE+1), OMGMAX (NABXBN (NX)),
     &             0.0, 0.0
 7012    CONTINUE
         WRITE (4,9504) EMAX, OMGMAX(NABXBN(NX)), 0.0, 0.0
 6850 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))
 9501 FORMAT (4(1P,1X,E10.3))
*9502 FORMAT (('# ',20(1X,1P,E10.3)))
 9502 FORMAT ((20(1X,1P,E10.3)))
*9503 FORMAT ((20(1X,1P,E10.3)))
 9503 FORMAT (2(1X,1P,E10.3),/,(20(1X,1P,E10.3)))
 9504 FORMAT (4(1X,1P,E10.3))
      STOP
      END

*$ CREATE LNNBLN.FOR
*COPY LNNBLN
*
*=== 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

