*
*----------------------------------------------------------------------*
*                                                                      *
*  Package of FLUKA user files for the explicite calculation of dose   *
*  rates from induced radioactivity.                                   *
*                                                                      *
*  This version dated 25.10.06 by Stefan Roesler.                      *
*                                                                      *
*----------------------------------------------------------------------*
*
*                             USRICALL
*                      (first and second step)
*
*
*  SDUM = TCOOL, TCOOLM, TCOOLH, TCOOLD, TCOOLY 
*  (first step)
*
*    WHAT(1-6)  Cooling time in seconds, minutes, hours, days, years
*               
*    Note: The card can be given repeatedly. 
*          Maximum number of cooling times = 8.
*    
*
*  SDUM = DUMPING
*  (first step)
*
*    WHAT(1)    = 0  gamma ray, positron, and electron emitter 
*               = 1  gamma ray emitter only
*               = 2  positron emitter only
*               = 3  electron emitter only
*               = 4  gamma ray and positron emitter only
*
*
*  SDUM = DUMPREG
*  (first step)
*
*    WHAT(1)    Lower bound of region index for isotope dumping
*
*    WHAT(2)    Upper bound of region index for isotope dumping
*
*    WHAT(3-4)  As WHAT(1-2), here further range of region indices
*
*    WHAT(5-6)  As WHAT(1-2), here further range of region indices
*
*    Note: The card can be given repeatedly. 
*
*
*  SDUM = SAMPREG
*  (second step)
*
*    WHAT(1)    Lower bound of region index for isotope sampling
*
*    WHAT(2)    Upper bound of region index for isotope sampling
*
*    WHAT(3-4)  As WHAT(1-2), here further range of region indices
*
*    WHAT(5-6)  As WHAT(1-2), here further range of region indices
*
*    Note: The card can be given repeatedly. 
*
*
*  SDUM = OUTPUT
*  (first and second step)
*
*    WHAT(1)    < 0: binary isotope dump file
*               Default: ASCII format
*
*  
*  SDUM = BIASING
*  (first step)
*
*    WHAT(1)  number of entries in dump file used for adjusting biasing
*             factors
*             Default: 200
*
*    WHAT(2)  maximum fraction of abundance of certain isotope within
*             WHAT(2) entries of the dump files above which biasing is
*             applied
*             Default: 0.2
*
*    WHAT(3)  Minimum weight for isotope dumping biasing
*             Default: 0.01
*    
*
*  SDUM = BIASREG
*  (first step)
*
*    WHAT(1)    Lower bound of region index for biasing in dumping
*
*    WHAT(2)    Upper bound of region index for biasing in dumping
*
*    WHAT(3-4)  As WHAT(1-2), here further range of region indices
*
*    WHAT(5-6)  As WHAT(1-2), here further range of region indices
*
*    Note: The card can be given repeatedly. 
*
*
*  SDUM = WCOOL
*  (second step)
*
*    WHAT(1-6)  Weight factors taken from output of first calculation
*               step corresponding to cooling times as defined in first
*               step
*
*----------------------------------------------------------------------*
*
*                                 SOURCE
*                              (second step)
*
*
*    WHAT(1)    Index of cooling time for which dose rates have to be
*               computed
*
*    WHAT(2)    = 0  gamma ray, positron, and electron emission
*               = 1  gamma ray emission only
*               = 2  positron emission only
*               = 3  electron emission only
*               = 4  gamma ray and positron emission only
*
*----------------------------------------------------------------------*
*$ CREATE USRINI.FOR
*COPY USRINI
*
*=== usrini ===========================================================*
*
      SUBROUTINE USRINI ( WHAT, SDUM )
*
*----------------------------------------------------------------------*
*                                                                      *
*  This routine is part of the package of user files for the           *
*  calculation of dose rates from induced activity.                    *
*                                                                      *
*  Initialization of histograms. Input of cooling times (in calculation*
*  of isotopes, 1st step) and weights (in calculation of dose rates,   *
*  2nd step) via USRICALL.                                             *
*                                                                      *
*  This version dated 01.08.06 by Stefan Roesler.                      *
*                                                                      *
*----------------------------------------------------------------------*
*
      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
      DIMENSION WHAT (6)
      CHARACTER SDUM*8

      LOGICAL LFIRST

      PARAMETER (LUNISO = 31)

      PARAMETER (MXCOOL = 8,
     &           IZMAX  = 100,
     &           IAMAX  = 257)
      COMMON /ISOSTA/ YLDISO(MXCOOL,2,IZMAX,IAMAX),
     &                WGTISO(MXCOOL,2,IZMAX,IAMAX),
     &                ICTISO(MXCOOL,2,IZMAX,IAMAX),
     &                IJISO(5,2,IZMAX,IAMAX),NTOISO(MXXRGN,MXCOOL)
      COMMON /COOLIN/ TCOOL(MXCOOL),WGTCOO(MXCOOL),KCOOL
      COMMON /SAMPLE/ VOLSAM,RHOSAM
      LOGICAL LBIAS,LBINRY,LREGDR,LNA24
      COMMON /DOSPAR/ FBIAS,WMINBI,FNA24,KBIAS,KDUMP,
     &                LREGDR(MXXRGN),LBIAS(MXXRGN),LNA24(MXXRGN),
     &                LBINRY

      DATA LFIRST /.TRUE./
*
*  Don't change the following line:
      LUSRIN = .TRUE.
*
* Initialization of histograms
*
      IF (LFIRST) THEN
         KCOOL = 0
         MCOOL = 0
         DO 1 ICOOL=1,MXCOOL
            TCOOL(ICOOL)  = ZERZER
            DO 2 IZ=1,IZMAX
               DO 3 IA=1,IAMAX
                  YLDISO(ICOOL,1,IZ,IA) = ZERZER
                  YLDISO(ICOOL,2,IZ,IA) = ZERZER
                  WGTISO(ICOOL,1,IZ,IA) = ZERZER
                  WGTISO(ICOOL,2,IZ,IA) = ZERZER
                  ICTISO(ICOOL,1,IZ,IA) = 0
                  ICTISO(ICOOL,2,IZ,IA) = 0
                  IF (ICOOL.LE.5) THEN
                     IJISO(ICOOL,1,IZ,IA)  = 0
                     IJISO(ICOOL,2,IZ,IA)  = 0
                  ENDIF
    3          CONTINUE
    2       CONTINUE
    1    CONTINUE
         DO 6 IREG=1,MXXRGN
            DO 13 ICOOL=1,MXCOOL
               NTOISO(IREG,ICOOL) = 0
   13       CONTINUE
            LREGDR(IREG) = .FALSE.
            LBIAS(IREG)  = .FALSE.
            LNA24(IREG)  = .FALSE.
    6    CONTINUE
         VOLSAM = ONEONE
         RHOSAM = ONEONE
         LBINRY = .FALSE.
         KDUMP  = 0
         FBIAS  = 0.2D0
         KBIAS  = 200
         WMINBI = 0.01D0
         FNA24  = 1.44927D0
         LFIRST = .FALSE.
      ENDIF
*
* Input of cooling times (1st step)
*
      IF (SDUM(1:5).EQ.'TCOOL') THEN
         FAC = ONEONE
         IF (SDUM(6:6).EQ.'M') THEN
            FAC = 60.0D0
         ELSEIF (SDUM(6:6).EQ.'H') THEN
            FAC = 3600.0D0
         ELSEIF (SDUM(6:6).EQ.'D') THEN
            FAC = 86400.0D0
         ELSEIF (SDUM(6:6).EQ.'Y') THEN
            FAC = 31536000.0D0
         ENDIF
         DO 4 I=1,6
            IF (WHAT(I).GT.ZERZER) THEN
               KCOOL = KCOOL+1
               IF (KCOOL.GT.MXCOOL) STOP ' KCOOL > MXCOOL !'
               TCOOL(KCOOL) = FAC*WHAT(I)
            ENDIF
    4    CONTINUE
      ENDIF
*
* Input of weights (2nd step)
*
      IF (SDUM(1:5).EQ.'WCOOL') THEN
         DO 5 I=1,6
            IF (WHAT(I).GT.ZERZER) THEN
               MCOOL = MCOOL+1
               IF (MCOOL.GT.MXCOOL) STOP ' MCOOL > MXCOOL !'
               WGTCOO(MCOOL) = WHAT(I)
            ENDIF
    5    CONTINUE
      ENDIF
*
* Type of emitter (1st step)
*
      IF (SDUM(1:7).EQ.'DUMPING') THEN
         KDUMP = INT(WHAT(1))
         IF ((KDUMP.LT.0).OR.(KDUMP.GT.4)) THEN
            WRITE(LUNOUT,*) ' Invalid parameter KDUMP ! ',KDUMP
            STOP
         ENDIF
      ENDIF
*
* Input of dumping (1st step) or sampling region numbers (2nd step)
*
      IF ((SDUM(1:7).EQ.'DUMPREG').OR.(SDUM(1:7).EQ.'SAMPREG')) THEN
         I0 = INT(WHAT(1))
         I1 = INT(WHAT(2))
         IF ((I0.GE.1).AND.(I1.GE.1).AND.(I1.GE.I0)) THEN
            DO 7 IREG=I0,I1
               LREGDR(IREG) = .TRUE.
    7       CONTINUE
         ENDIF
         I0 = INT(WHAT(3))
         I1 = INT(WHAT(4))
         IF ((I0.GE.1).AND.(I1.GE.1).AND.(I1.GE.I0)) THEN
            DO 8 IREG=I0,I1
               LREGDR(IREG) = .TRUE.
    8       CONTINUE
         ENDIF
         I0 = INT(WHAT(5))
         I1 = INT(WHAT(6))
         IF ((I0.GE.1).AND.(I1.GE.1).AND.(I1.GE.I0)) THEN
            DO 9 IREG=I0,I1
               LREGDR(IREG) = .TRUE.
    9       CONTINUE
         ENDIF
      ENDIF
*
* Material-specific information (volume, density)
*
      IF (SDUM(1:8).EQ.'MATERIAL') THEN
         VOLSAM = WHAT(1)
         RHOSAM = WHAT(2)
      ENDIF
*
* Output-specific parameters
*
      IF (SDUM(1:6).EQ.'OUTPUT') THEN
         IF (WHAT(1).LT.ZERZER) LBINRY = .TRUE.
      ENDIF
*
* Biasing in dumping isotopes (1st step)
*
*  1) Parameters
*
      IF (SDUM(1:7).EQ.'BIASING') THEN
         IF (WHAT(1).GT.ZERZER) THEN
            KBIAS = INT(WHAT(1))
         ENDIF
         IF (WHAT(2).GT.ZERZER) THEN
            FBIAS = WHAT(2)
         ENDIF
         IF (WHAT(3).GT.ZERZER) THEN
            WMINBI = WHAT(3)
         ENDIF
      ENDIF
*
*  2) Region numbers
*
      IF (SDUM(1:7).EQ.'BIASREG') THEN
         I0 = INT(WHAT(1))
         I1 = INT(WHAT(2))
         IF ((I0.GE.1).AND.(I1.GE.1).AND.(I1.GE.I0)) THEN
            DO 10 IREG=I0,I1
               LBIAS(IREG) = .TRUE.
   10       CONTINUE
         ENDIF
         I0 = INT(WHAT(3))
         I1 = INT(WHAT(4))
         IF ((I0.GE.1).AND.(I1.GE.1).AND.(I1.GE.I0)) THEN
            DO 11 IREG=I0,I1
               LBIAS(IREG) = .TRUE.
   11       CONTINUE
         ENDIF
         I0 = INT(WHAT(5))
         I1 = INT(WHAT(6))
         IF ((I0.GE.1).AND.(I1.GE.1).AND.(I1.GE.I0)) THEN
            DO 12 IREG=I0,I1
               LBIAS(IREG) = .TRUE.
   12       CONTINUE
         ENDIF
      ENDIF
*
* Additional factor to scale for underestimation in Al-->Na24 production
*
      IF (SDUM(1:7).EQ.'AL-NA24') THEN
         I0 = INT(WHAT(1))
         IF ((I0.GT.0).AND.(I0.LE.MXXRGN)) LNA24(I0) = .TRUE.
         I0 = INT(WHAT(2))
         IF ((I0.GT.0).AND.(I0.LE.MXXRGN)) LNA24(I0) = .TRUE.
         I0 = INT(WHAT(3))
         IF ((I0.GT.0).AND.(I0.LE.MXXRGN)) LNA24(I0) = .TRUE.
         I0 = INT(WHAT(4))
         IF ((I0.GT.0).AND.(I0.LE.MXXRGN)) LNA24(I0) = .TRUE.
         I0 = INT(WHAT(5))
         IF ((I0.GT.0).AND.(I0.LE.MXXRGN)) LNA24(I0) = .TRUE.
         I0 = INT(WHAT(6))
         IF ((I0.GT.0).AND.(I0.LE.MXXRGN)) LNA24(I0) = .TRUE.
      ENDIF

      RETURN
      END

*$ CREATE USROUT.FOR
*COPY USROUT
*
*=== usrout ===========================================================*
*
      SUBROUTINE USROUT ( WHAT, SDUM )
*
*----------------------------------------------------------------------*
*                                                                      *
*  This routine is part of the package of user files for the           *
*  calculation of dose rates from induced activity.                    *
*                                                                      *
*  Output of sampling statistics and histograms.                       *
*                                                                      *
*  This version dated 25.10.06 by Stefan Roesler.                      *
*                                                                      *
*----------------------------------------------------------------------*
*
      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
      INCLUDE '(CASLIM)'
*
      DIMENSION WHAT (6)
      CHARACTER SDUM*8

      PARAMETER (LUNISO = 31)

      PARAMETER (MXCOOL = 8,
     &           IZMAX  = 100,
     &           IAMAX  = 257)
      COMMON /ISOSTA/ YLDISO(MXCOOL,2,IZMAX,IAMAX),
     &                WGTISO(MXCOOL,2,IZMAX,IAMAX),
     &                ICTISO(MXCOOL,2,IZMAX,IAMAX),
     &                IJISO(5,2,IZMAX,IAMAX),NTOISO(MXXRGN,MXCOOL)
      COMMON /SAMPLE/ VOLSAM,RHOSAM
      COMMON /COOLIN/ TCOOL(MXCOOL),WGTCOO(MXCOOL),KCOOL

      DIMENSION DOSEPS(MXCOOL),NTOT(MXCOOL)
*
C     WRITE(LUNOUT,*)' Number of primary particles (requested): ',NCASES
C     WRITE(LUNOUT,*)'                             (sampled)  : ',NCASE
      WRITE(LUNOUT,*)' Total number of isotope entries in isodump.out: '
      DO 4 I=1,MXCOOL
         NTOT(I) = 0
         DO 7 J=1,MXXRGN
            IF (NTOISO(J,I).GT.0) THEN
               NTOT(I) = NTOT(I)+NTOISO(J,I)
               WRITE(LUNOUT,1000) I,J,NTOISO(J,I)
 1000          FORMAT(42X,I3,I5,I15)
            ENDIF
    7    CONTINUE
    4 CONTINUE

      WRITE(LUNOUT,*)' Parameters for USRICALL:'
      WRITE(LUNOUT,*)'    - all regions '
      NLINES = MXCOOL/6
      IF (NLINES.EQ.0) THEN
         WRITE(LUNOUT,1004) (DBLE(NTOT(I))/DBLE(NCASE),I=1,MXCOOL)
 1004    FORMAT(10X,1P,6E10.3)
      ELSE
         DO 6 ILINE=1,NLINES
            I0 = 1+(ILINE-1)*6
            I1 = I0+5
            WRITE(LUNOUT,1004) (DBLE(NTOT(I))/DBLE(NCASE),I=I0,I1)
    6    CONTINUE
         NLEFT = KCOOL-6*NLINES
         IF (NLEFT.GT.0) THEN
            I0 = I1+1
            I1 = I0+NLEFT-1
            WRITE(LUNOUT,1004) (DBLE(NTOT(I))/DBLE(NCASE),I=I0,I1)
         ENDIF
      ENDIF
      DO 8 J=1,MXXRGN
         IF (NTOISO(J,1).GT.0) THEN
            WRITE(LUNOUT,*)'    - regions no. ',J
            IF (NLINES.EQ.0) THEN
               WRITE(LUNOUT,1004) 
     &            (DBLE(NTOISO(J,I))/DBLE(NCASE),I=1,MXCOOL)
            ELSE
               DO 9 ILINE=1,NLINES
                  I0 = 1+(ILINE-1)*6
                  I1 = I0+5
                  WRITE(LUNOUT,1004) 
     &               (DBLE(NTOISO(J,I))/DBLE(NCASE),I=I0,I1)
    9          CONTINUE
               NLEFT = KCOOL-6*NLINES
               IF (NLEFT.GT.0) THEN
                  I0 = I1+1
                  I1 = I0+NLEFT-1
                  WRITE(LUNOUT,1004) 
     &               (DBLE(NTOISO(J,I))/DBLE(NCASE),I=I0,I1)
               ENDIF
            ENDIF
         ENDIF
    8 CONTINUE

      WRITE(LUNOUT,*) ' Material properties used for normalization ',
     &                '(volume,density):', VOLSAM,RHOSAM
      DO 1 ICOOL=1,MXCOOL
         DOSEPS(ICOOL) = 0.0D0
         WRITE(LUNOUT,1001) ICOOL
 1001    FORMAT(/,'Cooling time no.',I2)
         DO 2 IZ=1,IZMAX
            DO 3 IA=1,IAMAX
               DOSEPS(ICOOL) = DOSEPS(ICOOL)
     &            +WGTISO(ICOOL,1,IZ,IA)+WGTISO(ICOOL,2,IZ,IA)
               IF (YLDISO(ICOOL,1,IZ,IA).GT.ZERZER)
     &            WRITE(LUNOUT,1002) ' 1',IZ,IA,
     &               YLDISO(ICOOL,1,IZ,IA)/DBLE(NCASES)/VOLSAM/RHOSAM,
     &               WGTISO(ICOOL,1,IZ,IA)/
     &                  DBLE(MAX(ICTISO(ICOOL,1,IZ,IA),1)),
     &               ICTISO(ICOOL,1,IZ,IA),
     &               (IJISO(I,1,IZ,IA),I=1,5)
               IF (YLDISO(ICOOL,2,IZ,IA).GT.ZERZER)
     &            WRITE(LUNOUT,1002) ' 2',IZ,IA,
     &               YLDISO(ICOOL,2,IZ,IA)/DBLE(NCASES)/VOLSAM/RHOSAM,
     &               WGTISO(ICOOL,2,IZ,IA)/
     &                  DBLE(MAX(ICTISO(ICOOL,2,IZ,IA),1)),
     &               ICTISO(ICOOL,2,IZ,IA),
     &               (IJISO(I,2,IZ,IA),I=1,5)
 1002          FORMAT(A2,2I6,1P,2E15.4,I6,5I6)
    3       CONTINUE
    2    CONTINUE
         DOSEPS(ICOOL) = DOSEPS(ICOOL)/DBLE(NCASES)
    1 CONTINUE
      WRITE(LUNOUT,*) 'Dose estimate (point source, nSv/h)for r=10cm:'
      FAC = 1.D9
      DO 5 ICOOL=1,MXCOOL
         WRITE(LUNOUT,1003) ICOOL,FAC*1.D-8/7.0D0*DOSEPS(ICOOL)/100.0D0
 1003    FORMAT(I4,1P,E15.5)
    5 CONTINUE
*
      RETURN
      END

*$ CREATE USRRNC.FOR
*COPY USRRNC
*
*=== Usrrnc ===========================================================*
*
      SUBROUTINE USRRNC (ICRES,IBRES,ISRES,XA,YA,ZA,MREG,RULL,ICALL)

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
*                                                                      *
*  This routine is part of the package of user files for the           *
*  calculation of dose rates from induced activity.                    *
*                                                                      *
*  This routine dumps information on radioactive isotope production    *
*  into external file 'isodump.out'.                                   *
*                                                                      *
*  This version dated 19.01.03 by Stefan Roesler.                      *
*                                                                      *
*  Based on the default geoden routine:                                *
*                                                                      *
*----------------------------------------------------------------------*
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 2005-2005      by    Alfredo Ferrari & Paola Sala  *
*     All Rights Reserved.                                             *
*                                                                      *
*                                                                      *
*     USeR Residual NuClei:                                            *
*                                                                      *
*     Created on   06 april 2005   by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*     Last change on 06-apr-05     by    Alfredo Ferrari               *
*                                                                      *
*                                                                      *
*----------------------------------------------------------------------*
*
      LOGICAL LFIRST, LDUMP
 
      PARAMETER (LUNISO=31)

      PARAMETER (MXCOOL = 8,
     &           MAXISO = 20,
     &           MXBIAS = 20,
     &           IZMAX  = 100,
     &           IAMAX  = 257)
      COMMON /ISOSTA/ YLDISO(MXCOOL,2,IZMAX,IAMAX),
     &                WGTISO(MXCOOL,2,IZMAX,IAMAX),
     &                ICTISO(MXCOOL,2,IZMAX,IAMAX),
     &                IJISO(5,2,IZMAX,IAMAX),NTOISO(MXXRGN,MXCOOL)
      COMMON /COOLIN/ TCOOL(MXCOOL),WGTCOO(MXCOOL),KCOOL
      LOGICAL LBIAS,LBINRY,LREGDR,LNA24
      COMMON /DOSPAR/ FBIAS,WMINBI,FNA24,KBIAS,KDUMP,
     &                LREGDR(MXXRGN),LBIAS(MXXRGN),LNA24(MXXRGN),
     &                LBINRY
*
* flag to indicate if decay data have already been dumped
      DIMENSION ISODEC(2,IZMAX,IAMAX)
*
      DIMENSION IZTMP(IZMAX,IAMAX,MXCOOL,MAXISO),
     &          IATMP(IZMAX,IAMAX,MXCOOL,MAXISO),
     &          ISTMP(IZMAX,IAMAX,MXCOOL,MAXISO),
     &          YLDTMP(IZMAX,IAMAX,MXCOOL,MAXISO),
     &          NISTMP(IZMAX,IAMAX,MXCOOL),KTMP(IZMAX,IAMAX,MXCOOL),
     &          ICBIAS(2,IZMAX,IAMAX),WBIAS(MXBIAS)
*
* variable to temporarily store usrsuw output before dumping it
      DIMENSION IZDEC(MAXISO),IADEC(MAXISO),ISDEC(MAXISO),
     &          YLDDEC(MAXISO)
*
* interface variable for usrsuw
      DIMENSION JZDEC(MXCOOL,MAXISO),JADEC(MXCOOL,MAXISO),
     &          JSDEC(MXCOOL,MAXISO),ZLDDEC(MXCOOL,MAXISO),KISO(MXCOOL)
*
      DATA LFIRST /.TRUE./

      IF (LREGDR(MREG)) THEN
*
* open dump file
*
         IF (LFIRST) THEN
            IF (LBINRY) THEN
               OPEN (LUNISO,FILE='isodump.out',STATUS='UNKNOWN',
     &                                       FORM='UNFORMATTED')
            ELSE
               CALL OAUXFI('isodump.out',LUNISO,'UNKNOWN',IERR)
            ENDIF
            IF (IERR.GT.0) STOP 'Cannot open isotope dump file !'
            WRITE(LUNOUT,*) 
     &        ' GEORSN: Isotopes dumped for the following cooling times'
            WRITE(LUNOUT,*) (TCOOL(ICOOL),ICOOL=1,KCOOL)
            IF (LBINRY) THEN
               WRITE(LUNISO) KDUMP,KCOOL
            ELSE
               WRITE(LUNISO,'(2I3)') KDUMP,KCOOL
            ENDIF
            DO 51 ICOOL=1,KCOOL
               IF (LBINRY) THEN
                  WRITE(LUNISO) TCOOL(ICOOL)
               ELSE
                  WRITE(LUNISO,'(A3,I3,1P,E16.7)') 
     &               'tc ',ICOOL,TCOOL(ICOOL)
               ENDIF
   51       CONTINUE
            DO 16 I=1,IAMAX
               DO 17 J=1,IZMAX
                  ISODEC(1,J,I) = 0
                  ISODEC(2,J,I) = 0
                  ICBIAS(1,J,I) = 0
                  ICBIAS(2,J,I) = 0
                  DO 18 ICOOL=1,MXCOOL
                     NISTMP(J,I,ICOOL) = 0
                     KTMP(J,I,ICOOL)   = 0
                     DO 19 ISO=1,MAXISO
                        IZTMP(J,I,ICOOL,ISO)  = 0
                        IATMP(J,I,ICOOL,ISO)  = 0
                        ISTMP(J,I,ICOOL,ISO)  = 0
                        YLDTMP(J,I,ICOOL,ISO) = ZERZER
   19                CONTINUE
   18             CONTINUE
   17          CONTINUE
   16       CONTINUE
            DO 52 I=1,MXBIAS
               WBIAS(I) = ONEONE
   52       CONTINUE
            NBIAS  = 0
            NWBIAS = 0
            IF (LBIAS(MREG)) THEN
               WRITE(LUNOUT,*) ' GEORSN: biased dumping of isotopes ',
     &                         'requested with parameters '
               WRITE(LUNOUT,*) '         Fbias = ',FBIAS,
     &                         '   Kbias = ',KBIAS,'  Wminbi = ',WMINBI
               IF (KDUMP.EQ.ONEONE) THEN
                  WRITE(LUNOUT,*) ' GEORSN: gamma ray emitters only !'
               ELSEIF (KDUMP.EQ.TWOTWO) THEN
                  WRITE(LUNOUT,*) ' GEORSN: positron emitters only !'
               ELSEIF (KDUMP.EQ.THRTHR) THEN
                  WRITE(LUNOUT,*) ' GEORSN: electron emitters only !'
               ELSEIF (KDUMP.EQ.FOUFOU) THEN
                  WRITE(LUNOUT,*) 
     &               ' GEORSN: gamma ray and positron emitters only !'
               ENDIF
            ENDIF
            LFIRST = .FALSE.
         ENDIF
*
         INRES = IBRES-ICRES
         NMZ   = INRES-ICRES
         IF ((IBRES.GT.0).AND.(NMZ.GE.-4)) THEN
*
            IF (IJ.EQ.1) THEN
               JJ = 1
            ELSEIF (IJ.EQ.8) THEN
               JJ = 2
            ELSEIF (IJ.EQ.13) THEN
               JJ = 3
            ELSEIF (IJ.EQ.14) THEN
               JJ = 4
            ELSE
               JJ = 5
            ENDIF
*
* calculate intensities for the various cooling times
*
            JISO  = 0
            NCOOL = 0
            DO 10 ICOOL=1,KCOOL
               WEI  = ONEONE
               NISO = 0
               KISO(ICOOL) = 0
               IF (KTMP(ICRES,IBRES,ICOOL).EQ.0) THEN
                  CALL USRSUW_EVO(ICRES,IBRES,WEI,TCOOL(ICOOL),
     &                            NISO,IZDEC,IADEC,ISDEC,YLDDEC)
                  KTMP(ICRES,IBRES,ICOOL)   = 1
                  NISTMP(ICRES,IBRES,ICOOL) = NISO
                  IF (NISO.GT.0) THEN
                     DO 200 ISO=1,NISO
                        IZTMP(ICRES,IBRES,ICOOL,ISO)  = IZDEC(ISO)
                        IATMP(ICRES,IBRES,ICOOL,ISO)  = IADEC(ISO)
                        ISTMP(ICRES,IBRES,ICOOL,ISO)  = ISDEC(ISO)
                        YLDTMP(ICRES,IBRES,ICOOL,ISO) = YLDDEC(ISO)
  200                CONTINUE
                  ENDIF
               ELSE
                  NISO = NISTMP(ICRES,IBRES,ICOOL)
                  DO 201 ISO=1,MAXISO
                     IZDEC(ISO)  = IZTMP(ICRES,IBRES,ICOOL,ISO)
                     IADEC(ISO)  = IATMP(ICRES,IBRES,ICOOL,ISO)
                     ISDEC(ISO)  = ISTMP(ICRES,IBRES,ICOOL,ISO)
                     YLDDEC(ISO) = YLDTMP(ICRES,IBRES,ICOOL,ISO)
  201             CONTINUE
               ENDIF
               IF (NISO.GT.0) THEN
                  IF (NISO.GT.JISO) JISO = NISO
                  KISO(ICOOL) = NISO
                  NCOOL = NCOOL+1
                  DO 11 ISO=1,NISO
                     JZDEC(ICOOL,ISO)  = IZDEC(ISO)
                     JADEC(ICOOL,ISO)  = IADEC(ISO)
                     JSDEC(ICOOL,ISO)  = ISDEC(ISO)
                     ZLDDEC(ICOOL,ISO) = YLDDEC(ISO)
   11             CONTINUE
               ENDIF
   10       CONTINUE
*
* dump intensities
*
*   Niso, Kiso(Icool) number of daugher isotopes for this cooling time
*   Jiso              max. number of daugher isotopes for all cooling times
*   Ncool             number of cooling times with non-zero yields
*   Jzdec(Icool,Iso), Jadec(Icool,Iso)
*                     mass number and charge for this cooling time and 
*                     daughter isotope
*   Jsdec(Icool,Iso)  isomer flag for this cooling time and daughter isotope
*   Zlddec(Icool,Iso) yield for this cooling time and daughter isotope
*
*   Isodec(Is,Iz,Ia)  flag to indicate if isotope data have already been
*                     dumped
*   Mcool = Ncool     isotope decay data follow
*         = -Ncool    isotope decay data have been dumped previously
*
   12       CONTINUE
            IF (JISO.GT.0) THEN
*   only one decay product
               IF (JISO.EQ.1) THEN
*
* artificial biasing for Cu64
C                 IF ((JZDEC(1,1).EQ.29).AND.(JADEC(1,1).EQ.64)) THEN
C                    R64CU = 0.1D0
C                    RKILL = FLRNDM(R64CU)
C                    IF (RKILL.GT.R64CU) GOTO 20
C                    RULL = RULL/R64CU
C                 ENDIF
* biasing of number of entries in dump file
                  IF (LBIAS(MREG)) THEN
                     IF (NBIAS.EQ.KBIAS) THEN
                        write(lunout,*) ' Checking biasing factors '
                        DO 53 I=1,IZMAX
                           DO 54 J=I,IAMAX
                              FRAC = DBLE(ICBIAS(1,I,J))/DBLE(NBIAS)
                              IF (FRAC.GT.FBIAS) THEN
                                 IF (ICBIAS(2,I,J).EQ.0) THEN
                                    NWBIAS = NWBIAS+1
                                    ICBIAS(2,I,J) = NWBIAS
                                 ENDIF
                                 WBIAS(ICBIAS(2,I,J)) = 
     &                              MAX(WBIAS(ICBIAS(2,I,J))/TWOTWO,
     &                                                       WMINBI)
                                 write(lunout,*) I,J,ICBIAS(2,I,J),
     &                                    WBIAS(ICBIAS(2,I,J)),FRAC
                              ENDIF
                              ICBIAS(1,I,J) = 0
   54                      CONTINUE
   53                   CONTINUE
                        NBIAS = 0
                     ENDIF
                     IF (ICBIAS(2,JZDEC(1,1),JADEC(1,1)).GT.0) THEN
                        FRAC = FLRNDM(RULL)
                        IF (FRAC.GT.WBIAS(ICBIAS(2,JZDEC(1,1),
     &                                             JADEC(1,1)))) GOTO 20
                        RULL = RULL
     &                         /WBIAS(ICBIAS(2,JZDEC(1,1),JADEC(1,1)))
                     ENDIF
                     ICBIAS(1,JZDEC(1,1),JADEC(1,1)) = 
     &                  ICBIAS(1,JZDEC(1,1),JADEC(1,1))+1
                     NBIAS = NBIAS+1
                  ENDIF
* apply biasing factor for Na24 production from Al
C                 IF ((JZDEC(1,1).EQ.11).AND.(JADEC(1,1).EQ.24).AND.
C    &                                             LNA24(MREG)) THEN
C                    RULL = RULL*FNA24
C                 ENDIF
*
                  MCOOL = -NCOOL
                  IF (ISODEC(JSDEC(1,1),JZDEC(1,1),JADEC(1,1)).EQ.0)
     &               MCOOL = NCOOL
                  IF (LBINRY) THEN
                     IFLAG = 1
                     WRITE(LUNISO) IFLAG
                     WRITE(LUNISO) MCOOL,JZDEC(1,1),JADEC(1,1),
     &                             JSDEC(1,1),MREG,XA,YA,ZA,RULL
                  ELSE
                     WRITE(LUNISO,1001)'#1',MCOOL,JZDEC(1,1),JADEC(1,1),
     &                                 JSDEC(1,1),MREG,XA,YA,ZA,RULL
 1001                FORMAT(A2,I3,2I4,I2,I5,1P,4E13.5)
                  ENDIF
                  IF ((JZDEC(1,1).LT.0).OR.(JZDEC(1,1).GT.IZMAX).OR.
     &                (JADEC(1,1).LT.0).OR.(JADEC(1,1).GT.IAMAX)) THEN
                     WRITE(LUNOUT,*) ' GEORSN: Z or A out of range ! ',
     &                               JZDEC(1,1),JADEC(1,1),JSDEC(1,1)
                     GOTO 20
                  ENDIF
                  IJISO(JJ,JSDEC(1,1),JZDEC(1,1),JADEC(1,1)) = 
     &               IJISO(JJ,JSDEC(1,1),JZDEC(1,1),JADEC(1,1))+1
                  DO 5 ICOOL=1,NCOOL
                     NTOISO(MREG,ICOOL) = NTOISO(MREG,ICOOL)+1
                     YLDISO(ICOOL,JSDEC(1,1),JZDEC(1,1),JADEC(1,1)) = 
     &                  YLDISO(ICOOL,JSDEC(1,1),JZDEC(1,1),JADEC(1,1))+
     &                  RULL*ZLDDEC(ICOOL,1)
                     WGTISO(ICOOL,JSDEC(1,1),JZDEC(1,1),JADEC(1,1)) = 
     &                  WGTISO(ICOOL,JSDEC(1,1),JZDEC(1,1),JADEC(1,1))+
     &                  RULL
                     ICTISO(ICOOL,JSDEC(1,1),JZDEC(1,1),JADEC(1,1)) = 
     &                  ICTISO(ICOOL,JSDEC(1,1),JZDEC(1,1),JADEC(1,1))+1
    5             CONTINUE
                  IF (ISODEC(JSDEC(1,1),JZDEC(1,1),JADEC(1,1)).EQ.0)THEN
                     NLINES = NCOOL/6
                     IF (NLINES.EQ.0) THEN
                        IF (LBINRY) THEN
                           WRITE(LUNISO) 
     &                        (ZLDDEC(ICOOL,1),ICOOL=1,NCOOL)
                        ELSE
                           WRITE(LUNISO,1002) 
     &                        (ZLDDEC(ICOOL,1),ICOOL=1,NCOOL)
 1002                      FORMAT(1P,6E12.5)
                        ENDIF
                     ELSE
                        DO 13 ILINE=1,NLINES
                           I0 = 1+(ILINE-1)*6
                           I1 = I0+5
                           IF (LBINRY) THEN
                              WRITE(LUNISO) 
     &                           (ZLDDEC(ICOOL,1),ICOOL=I0,I1)
                           ELSE
                              WRITE(LUNISO,1002) 
     &                           (ZLDDEC(ICOOL,1),ICOOL=I0,I1)
                           ENDIF
   13                   CONTINUE
                        NLEFT = NCOOL-6*NLINES
                        IF (NLEFT.GT.0) THEN
                           I0 = I1+1
                           I1 = I0+NLEFT-1
                           IF (LBINRY) THEN
                              WRITE(LUNISO) 
     &                           (ZLDDEC(ICOOL,1),ICOOL=I0,I1)
                           ELSE
                              WRITE(LUNISO,1002) 
     &                           (ZLDDEC(ICOOL,1),ICOOL=I0,I1)
                           ENDIF
                        ENDIF
                     ENDIF
                     ISODEC(JSDEC(1,1),JZDEC(1,1),JADEC(1,1)) = 1
                  ENDIF
*   several decay products
               ELSE
                  MCOOL = -NCOOL
                  LDUMP = .FALSE.
                  DO 60 ICOOL=1,MXCOOL
                     DO 61 ISO=1,KISO(ICOOL)
                        IZ = JZDEC(ICOOL,ISO)
                        IA = JADEC(ICOOL,ISO)
                        IS = JSDEC(ICOOL,ISO)
                        IF ((IZ.NE.ICRES).OR.(IA.NE.IBRES)) THEN
C                          WRITE(LUNOUT,*) ' GEORSN: decay data ',
C    &                        'cannot be dumped !',ICRES,IBRES,IZ,IA
                           LDUMP = .TRUE.
                           MCOOL = NCOOL
                           GOTO 62
                        ENDIF
                        IF (ISODEC(IS,IZ,IA).EQ.0) THEN
                           LDUMP = .TRUE.
                           ISODEC(IS,IZ,IA) = 1
                           MCOOL = NCOOL
                        ENDIF
   61                CONTINUE
   60             CONTINUE
**for the moment dump anyway
C                 LDUMP = .TRUE.
C                 MCOOL = NCOOL
**
   62             CONTINUE
                  IF (LBINRY) THEN
                     IFLAG = 2
                     WRITE(LUNISO) IFLAG
                     WRITE(LUNISO) MCOOL,ICRES,IBRES,MREG,XA,YA,ZA,RULL
                  ELSE
                     WRITE(LUNISO,1004)
     &                  '#2',MCOOL,ICRES,IBRES,MREG,XA,YA,ZA,RULL
 1004                FORMAT(A2,I3,2I4,I5,1P,4E13.5)
                  ENDIF
                  DO 63 ICOOL=1,MXCOOL
                     IF (KISO(ICOOL).GT.0) THEN
                        NTOISO(MREG,ICOOL) = NTOISO(MREG,ICOOL)
     &                                       +KISO(ICOOL)
                        DO 6 ISO=1,KISO(ICOOL)
                           IZ = JZDEC(ICOOL,ISO)
                           IA = JADEC(ICOOL,ISO)
                           IS = JSDEC(ICOOL,ISO)
                           YLDISO(ICOOL,IS,IZ,IA) = 
     &                        YLDISO(ICOOL,IS,IZ,IA)+
     &                        RULL*ZLDDEC(ICOOL,ISO)
                           WGTISO(ICOOL,IS,IZ,IA) = 
     &                        WGTISO(ICOOL,IS,IZ,IA)+RULL
                           ICTISO(ICOOL,IS,IZ,IA) = 
     &                        ICTISO(ICOOL,IS,IZ,IA)+1
                           IF (ICOOL.EQ.1) THEN
                              IJISO(JJ,IS,IZ,IA) = 
     &                           IJISO(JJ,IS,IZ,IA)+1
                           ENDIF
    6                   CONTINUE
                     ENDIF
   63             CONTINUE
                  IF (LDUMP) THEN
                     IF (LBINRY) THEN
                        WRITE(LUNISO)(KISO(ICOOL),ICOOL=1,MXCOOL)
                     ELSE
                        WRITE(LUNISO,'(100I3)')
     &                     (KISO(ICOOL),ICOOL=1,MXCOOL)
                     ENDIF
                     DO 14 ICOOL=1,MXCOOL
                        IF (KISO(ICOOL).GT.0) THEN
                           NLINES = KISO(ICOOL)/3
                           IF (NLINES.EQ.0) THEN
                              IF (LBINRY) THEN
                                 WRITE(LUNISO) 
     &                              ICOOL,(JZDEC(ICOOL,ISO),
     &                              JADEC(ICOOL,ISO),JSDEC(ICOOL,ISO),
     &                              ZLDDEC(ICOOL,ISO),ISO=1,KISO(ICOOL))
                              ELSE
                                 WRITE(LUNISO,1003) 
     &                              ICOOL,(JZDEC(ICOOL,ISO),
     &                              JADEC(ICOOL,ISO),JSDEC(ICOOL,ISO),
     &                              ZLDDEC(ICOOL,ISO),ISO=1,KISO(ICOOL))
 1003                            FORMAT(I3,3(2I4,I2,E12.5))
                              ENDIF
                           ELSE
                              DO 15 ILINE=1,NLINES
                                 I0 = 1+(ILINE-1)*3
                                 I1 = I0+2
                                 IF (LBINRY) THEN
                                    WRITE(LUNISO)
     &                                ICOOL,(JZDEC(ICOOL,ISO),
     &                                JADEC(ICOOL,ISO),JSDEC(ICOOL,ISO),
     &                                ZLDDEC(ICOOL,ISO),ISO=I0,I1)
                                 ELSE
                                    WRITE(LUNISO,1003)
     &                                ICOOL,(JZDEC(ICOOL,ISO),
     &                                JADEC(ICOOL,ISO),JSDEC(ICOOL,ISO),
     &                                ZLDDEC(ICOOL,ISO),ISO=I0,I1)
                                 ENDIF
   15                         CONTINUE
                              NLEFT = KISO(ICOOL)-3*NLINES
                              IF (NLEFT.GT.0) THEN
                                 I0 = I1+1
                                 I1 = I0+NLEFT-1
                                 IF (LBINRY) THEN
                                    WRITE(LUNISO)
     &                                ICOOL,(JZDEC(ICOOL,ISO),
     &                                JADEC(ICOOL,ISO),JSDEC(ICOOL,ISO),
     &                                ZLDDEC(ICOOL,ISO),ISO=I0,I1)
                                 ELSE
                                    WRITE(LUNISO,1003)
     &                                ICOOL,(JZDEC(ICOOL,ISO),
     &                                JADEC(ICOOL,ISO),JSDEC(ICOOL,ISO),
     &                                ZLDDEC(ICOOL,ISO),ISO=I0,I1)
                                 ENDIF
                              ENDIF
                           ENDIF
                        ENDIF
   14                CONTINUE
c                 ELSE
c                    DO 8 ICOOL=1,MXCOOL
c                       IF (KISO(ICOOL).GT.0) THEN
c                          NTOISO(MREG,ICOOL) = NTOISO(MREG,ICOOL)
c    &                                          +KISO(ICOOL)
c                          WRITE(LUNISO,1005) ICOOL,KISO(ICOOL),
c    &                        (JZDEC(ICOOL,ISO),JADEC(ICOOL,ISO),
c    &                         JSDEC(ICOOL,ISO),ISO=1,KISO(ICOOL))
c1005                      FORMAT(2I3,100(2I4,I2))
c                       ENDIF
c   8                CONTINUE
                  ENDIF
               ENDIF
            ENDIF

         ELSE
            IF (IBRES.LE.0) THEN
               WRITE(LUNOUT,*) ' GEORSN: IBRES < 0 ! ',
     &                         ICRES,IBRES,IJ,ICALL,RULL
            ELSE
               WRITE(LUNOUT,*) ' GEORSN: INRES-ICRES < -5 ! ',
     &                         IBRES,ICRES,INRES
            ENDIF
         ENDIF
* d, 3-H, 3-He and 4-He (no gamma emitters)
C        IF (ICALL.LT.0) THEN
C           WRITE(LUNOUT,*) ' GEORSN: ICALL < 0 !'
C        ENDIF
      ENDIF
   20 CONTINUE

      RETURN
*=== End of subroutine Usrrnc =========================================*
      END

**sr modified from stand-alone postprocessor
C     PROGRAM USRSUWEV

*$ CREATE USRSUW_EVO.FOR
*COPY USRSUW_EVO
*
*=== usrsuw_evo =======================================================*
*
      SUBROUTINE USRSUW_EVO(IZEVO,IAEVO,WEI,TDCY,
     &                      NISO,IZDEC,IADEC,ISODEC,YLDDEC)

************************************************************************
* Yield of daughter isotopes which are gamma emitters for certain      *
* colling time of the mother isotope.                                  *
*                                                                      *
* Input:    IZEVO,IAEVO charge/mass number of mother isotope           *
*           WEI         weight of  mother isotop                       *
*           TDCY        cooling time in seconds                        *
*                                                                      *
* Output:   NISO        number of daughter isotopes                    *
*           IZDEC,IADEC charge/mass number of daughter isotopes        *
*                       (isomer if Iadec>260)                          *
*           ISODEC      isomer flag of daughter isotopes               *
*           YLDDEC      yield of daughter isotopes                     *
*                                                                      *
* This version dated 29.7.06   is written by S. Roesler.               *
************************************************************************
**

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1995-2006      by    Alfredo Ferrari & Paola Sala  *
*     All Rights Reserved.                                             *
*                                                                      *
*                                                                      *
*     USer Residual scoring SUm and Weighting and EVolution:           *
*                                                                      *
*     Created on    12 july 1995   by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*     Last change on 09-apr-06     by    Alfredo Ferrari               *
*                                                                      *
*                                                                      *
*----------------------------------------------------------------------*
*
      INCLUDE '(RDCYCM)'
      INCLUDE '(USRSNC)'

**sr modified from stand-alone postprocessor
      LOGICAL LFIRST
**
*
      CHARACTER FILE*80,FILE1*80,RUNTIT*80,RUNTIM*32,CHSTAT*10,CHERR*3,
     &          CHISMR*10,YES*1
*
**sr modified from stand-alone postprocessor
C     PARAMETER ( MXDUMM = 500000 )
      PARAMETER ( MXDUMM =2000000)
**
      PARAMETER ( NMXGRP = 1000 )
      PARAMETER ( NMXEBN = 1000 )
      PARAMETER ( MXITHR = 1200 )
      PARAMETER ( MXIDAU =  110 )
      PARAMETER ( MXACTV =  100 * 100 * MXIDAU )
*
      DIMENSION GMSTOR (MXDUMM), GDSTOR (MXDUMM), GBSTOR (MXDUMM),
     &          GESTOR (MXDUMM), GCSTOR (MXDUMM),
     &          TOTTOT (MXRSNC), TOTERR (MXRSNC),
     &          BIRHLP (MXIRRD), TIRHLP (MXIRRD),
     &          KAURSN (MXRSNC), JX (MXRSNC), ACCUMA (NMXEBN+NMXGRP),
     &          ACCUMZ (NMXEBN+NMXGRP), RCCUMA (NMXEBN+NMXGRP),
     &          RCCUMZ (NMXEBN+NMXGRP), TDECAY (MXRSNC),
     &          KZURSN (MXRSNC), IABEG (100), IAEND (100)
      REAL WEIPRI, BEAINT, TIMEIR, TDECAY, VURHLP, GMSTOR, TOTA,
     &     TOTB  , RCCUMA, RCCUMZ, BIRHLP, TIRHLP
      DIMENSION IAD   ( MXIDAU ), IZD ( MXIDAU ), ISD ( MXIDAU ),
     &          ACTIV ( MXACTV ), BRCTOT ( MXIDAU )
*
      DIMENSION IATHR  (MXITHR), IZTHR  (MXITHR), ACTHR  (MXITHR),
     &          ERTHR  (MXITHR), IAIZTH (MXITHR), ICORTH (MXITHR),
     &          IMTHR  (MXITHR)
      DIMENSION T12DDA (4), BRDDAU (4), KKA (4), KKZ (4)
**sr modified from stand-alone postprocessor
      DOUBLE PRECISION WEI,TDCY
**
*
      CHARACTER CHSYM*2
      LOGICAL LISMRS, LSTATI, LCHECK, LPATCH, LFIEND, LEVOLN, LOFFEV,
     &        LOLDIR
      LOGICAL LPRINT
**sr
**sr modified from stand-alone postprocessor
C     PARAMETER (MXIRRA =33500, MXCOOL = 1, MXFILE = 4000, MXDET = 50)
      PARAMETER (MXIRRA = 5000, MXCOOL = 1, MXFILE = 170, MXDET = 1)
**
      LOGICAL LCYCLE
      DOUBLE PRECISION PRTSEC,TIRRAD,TDECY,ADDNUC,ADDERR,TOTWGT,WEIGHT
      CHARACTER CINFLE*80,CRUNTT*80,CRUNTM*32
      COMMON /IRRPRP/ PRTSEC,TIRRAD,TDECY(MXCOOL,MXIRRA),TOTWGT,
     &                WEIGHT(MXIRRA),NIRRAD,NCOOL,LCYCLE,
     &                CINFLE(MXFILE),NFILE,CRUNTT,CRUNTM,NPRIM
      DIMENSION ADDNUC(MXDET,100,520),ADDERR(MXDET,100,520)
**sr modified from stand-alone postprocessor
      PARAMETER (MAXISO=20)
      DOUBLE PRECISION YLDDEC
      DIMENSION IZDEC(MAXISO),IADEC(MAXISO),ISODEC(MAXISO),
     &          YLDDEC(MAXISO)
**
**
*
      EXTERNAL TIM1O2, BDNOPT
      DIMENSION IASTLN (100), RESNUC (100,260), RESERR (100,260)
      DATA IASTLN /
     &     1,       4,       7,       9,      11,      12,
     &    14,      16,      19,      20,      23,      24,
     &    27,      28,      31,      32,      35,      40,
     &    39,      40,      45,      48,      51,      52,
     &    55,      56,      59,      59,      64,      65,
     &    70,      73,      75,      79,      80,      84,
     &    86,      88,      89,      91,      93,      96,
     &    97,     101,     103,     107,     108,     113,
     &   115,     118,     122,     128,     127,     131,
     &   133,     137,     139,     140,     141,     144,
     &   145,     150,     152,     157,     159,     163,
     &   165,     167,     169,     173,     175,     179,
     &   181,     184,     186,     190,     192,     195,
     &   197,     201,     204,     207,     209,     209,
     &   210,     222,     223,     226,     227,     232,
     &   231,     238,     237,     244,     243,     247,
     &   247,     251,     254,     257 /
*
**sr modified from stand-alone postprocessor
      DATA LFIRST /.TRUE./
**

      CALL CMSPPR
*
*     Note that all addresses are for zero index (KBURSN,KAURSN)
*
      DO 5555 I = 1, MXDUMM
         GMSTOR (I) = 0.E+00
         GBSTOR (I) = ZERZER
         GDSTOR (I) = ZERZER
         GCSTOR (I) = ZERZER
         GESTOR (I) = ZERZER
 5555 CONTINUE
      DO 4444 I = 1, MXRSNC
         TOTTOT (I) = ZERZER
         TOTERR (I) = ZERZER
         NISMRS (I) = 0
 4444 CONTINUE
      CHERR  = '+/-'
      NCTOT  = 0
      WCTOT  = ZERZER
      MBATCH = 0
**sr
      NFILE = 0
**
      LPRINT = .FALSE.
      LSTATI = .FALSE.
      LISMRS = .FALSE.
      LEVOLN = .FALSE.
      K0     = 0
      KK0    = 0
      KI0    = 0
      KI1    = 0
      KMAX   = 0
      NRNDC  = 0
      LQOLD  = 0
**sr modified from stand-alone postprocessor
C     WRITE (*,*)' Do you want to patch isomers (def=yes)?'
C     READ  (*,'(A)') YES
      YES = 'y'
**
      LPATCH = YES .NE. 'N' .AND. YES .NE. 'n'
**sr modified from stand-alone postprocessor
C     WRITE (*,*) ' Multiplicative factor (=part/s):'
C     READ  (*,*) FACMUL
      FACMUL = 1.0D0
C     CALL ISOINI
**
**sr modified from stand-alone postprocessor
      KFILE = 0
**
1111  CONTINUE
      LFIEND = .FALSE.
**sr modified from stand-alone postprocessor
      KFILE = KFILE+1
C     WRITE(*,'('' Type the input file: '',$)')
C     READ (*,'(A)')FILE
C     LQ     = LNNBLN (FILE)
C     IF (LQ .LE. 0) GO TO 2000
      IF (KFILE.GT.1) GO TO 2000
C     OPEN (UNIT=1,FILE=FILE,STATUS='OLD',FORM='UNFORMATTED',ERR=1111)
**
      IRECRD = 0
**sr modified from stand-alone postprocessor
C     READ (1) RUNTIT, RUNTIM, WEIPRI, NCASE
**
      IRECRD = IRECRD + 1
**sr modified from stand-alone postprocessor
C     WRITE (*,*) RUNTIT
C     WRITE (*,*) RUNTIM
C     WRITE (*,*) WEIPRI
C     WRITE (*,*) ABS (NCASE)
      WEIPRI = 1.0D0
      NCASE  = 1.0D0
**
**sr
      NFILE  = NFILE+1
      IF (NFILE.GT.MXFILE) STOP ' nfile > mxfile !'
      CINFLE(NFILE) = FILE
      IF (NFILE.EQ.1) THEN
         NPRIM  = NCASE
         TOTWGT = WEIPRI
         CRUNTT = RUNTIT
         CRUNTM = RUNTIM
      ENDIF
**
**sr modified from stand-alone postprocessor
C     LOFFEV = INDEX ( RUNTIM, 'Time Evolution' ) .GT. 0
C     LEVOLN = LOFFEV .OR. NCASE .LT. 0
**
      LOLDIR = NCASE  .GT. 0
      NCASE  = ABS (NCASE)
      IF ( LEVOLN ) THEN
         IF ( LOLDIR ) THEN
            READ (1) BEAINT, TIMEIR
            IRECRD = IRECRD + 1
            NIRRDT = 1
            BIRRDT (1) = BEAINT
            TIRRDT (1) = TIMEIR
         ELSE
            READ (1) NIRRDT, ( BIRHLP(IR), TIRHLP(IR), IR = 1, NIRRDT )
            IRECRD = IRECRD + 1
            BEAINT = ZERZER
            TPREV  = ZERZER
            DO IR = 1, NIRRDT
               BIRRDT (IR) = BIRHLP (IR)
               TIRRDT (IR) = TIRHLP (IR)
               BEAINT = BEAINT + BIRRDT (IR) * ( TIRRDT (IR) - TPREV )
               TPREV  = TIRRDT (IR)
            END DO
            TIMEIR = TIRRDT (NIRRDT)
            BEAINT = BEAINT / TIMEIR
         END IF
         WRITE (*,*) ' Intensity:', BEAINT, ' pr/s'
         WRITE (*,*) ' T_irrad. :', TIMEIR, ' s'
         WRITE (*,*)
     &   ' Are you sure you want to evolve an evolution file??'
      END IF
      NCTOT  = NCTOT  + NCASE
      WCTOT  = WCTOT  + DBLE(WEIPRI)
      MBATCH = MBATCH + 1
      KLAST  = 0
      KMAX   = 0
      NRNDC  = 0
      MRNDC  = 0
      NRNMN  = 0
*  +-------------------------------------------------------------------*
*  |                            Read from binary file the results!!!
      DO 6600 IX = 1,1000
         NRN = IX - NRNMN
**sr modified from stand-alone postprocessor
C        READ (1,ERR=6350,END=1000) CHSTAT,IJX
         CHSTAT = '          '
**
         IRECRD = IRECRD + 1
         IF ( CHSTAT .EQ. 'STATISTICS') THEN
            LSTATI = .TRUE.
            NRNDC  = MRNDC
            IF ( .NOT. LEVOLN ) GO TO 1000
            NRNMN = NRNMN + 1
            DO J = 1, NRNDC
               READ (1)
               READ (1)
               READ (1)
               READ (1)
               READ (1)
               READ (1)
               IRECRD = IRECRD + 6
               IF ( LISMRS ) THEN
                  READ (1)
                  IRECRD = IRECRD + 1
               END IF
            END DO
            MRNDC = 0
            GO TO 6600
         END IF
         IRECRD = IRECRD - 1
 6350    CONTINUE
*        BACKSPACE 1
**sr modified from stand-alone postprocessor
C        REWIND 1
C        DO IR = 1, IRECRD
C           READ (1)
C        END DO
C        READ (1,ERR=6400,END=1000)
C    &                MRN, TIURSN(NRN), ITURSN(NRN), NRURSN(NRN),
C    &                VURHLP, IMRHGH (NRN), IZRHGH (NRN), NMZMNM
         MRN = 1
         TIURSN(NRN) = '          '
         ITURSN(NRN) = 1
         NRURSN(NRN) = 1
         VURHLP      = 1.0D0
         NMZMNM      = -5
         IMRHGH(NRN) = IAEVO-2*IZEVO+6-NMZMNM
         IZRHGH(NRN) = IZEVO+4
         INEVO       = IAEVO-IZEVO
         IMEVO       = INEVO-IZEVO-NMZMNM
**
         IRECRD = IRECRD + 1
         VURSNC (NRN) = VURHLP
**sr modified from stand-alone postprocessor
C        IF ( LEVOLN ) READ (1) TDECAY (NRN)
C        IF ( LEVOLN ) IRECRD = IRECRD + 1
**
         GO TO 6500
 6400    CONTINUE
*        BACKSPACE 1
         REWIND 1
         DO IR = 1, IRECRD
            READ (1)
         END DO
         READ (1,ERR=6450) CHSTAT,IJX
         IRECRD = IRECRD + 1
         IF ( CHSTAT .EQ. 'STATISTICS' ) THEN
            LSTATI = .TRUE.
            NRNDC  = MRNDC
            IF ( .NOT. LEVOLN ) GO TO 1000
            NRNMN = NRNMN + 1
            DO J = 1, NRNDC
               READ (1)
               READ (1)
               READ (1)
               READ (1)
               READ (1)
               READ (1)
               IRECRD = IRECRD + 6
               IF ( LISMRS ) THEN
                  READ (1)
                  IRECRD = IRECRD + 1
               END IF
            END DO
            MRNDC = 0
            GO TO 6600
         END IF
 6450    CONTINUE
         STOP ' 2ND-RECORD'
 6500    CONTINUE
         JX (IX-NRNMN) = NRN
         NRNMX   = NRN
         MRNDC   = MRNDC + 1
***************************************************************
*********** Note KBURSN is used for 0 address *****************
***************************************************************
         KBURSN (NRN) = KLAST
         K0     = KLAST + 1
         K1     = IMRHGH (NRN) * IZRHGH (NRN) + K0 - 1
         KLAST  = K1
**sr modified from stand-alone postprocessor
C        READ (1) (GMSTOR(J), J = K0, K1)
         IDX = IZEVO+IZRHGH(NRN)*(IMEVO-1)
         IF (IDX.LE.0) THEN
            WRITE(11,*) ' IZEVO,IAEVO,IMEVO,IZRHGH(NRN) ',
     &                   IZEVO,IAEVO,IMEVO,IZRHGH(NRN)
            STOP 'USRSUW: IDX .le. 0 !'
         ENDIF
         GMSTOR(IDX) = WEI
**
         IRECRD = IRECRD + 1
         KMAX   = MAX ( KMAX, K1 )
**sr modified from stand-alone postprocessor
C        READ (1,ERR=6530,END=6525) CHISMR,NISMRS(NRN)
C        IRECRD = IRECRD + 1
C        IF ( CHISMR .EQ. 'ISOMERS:  ' ) THEN
C           KI0   = K1 + 1
C           KI1   = K1 + NISMRS (NRN)
C           KLAST = KI1
C           KMAX  = MAX ( KMAX, KI1 )
C           IF ( KI1 .GE. KI0 ) THEN
C              READ (1) (GMSTOR(J), J = KI0, KI1)
C           ELSE
C              READ (1)
C           END IF
C           IRECRD = IRECRD + 1
C           LISMRS = .TRUE.
C           LPATCH = .FALSE.
C           GO TO 6535
C        ELSE
*           WRITE (*,*) ' CHISMR  ', CHISMR
C           NISMRS(NRN) = 0
C           IRECRD = IRECRD - 1
C           GO TO 6530
C        END IF
**
 6525    CONTINUE
         LFIEND = .TRUE.
 6530    CONTINUE
*        BACKSPACE 1
**sr modified from stand-alone postprocessor
C        REWIND 1
C        DO IR = 1, IRECRD
C           READ (1)
C        END DO
**
         KI1    = K1
*+++++++++++++++ Patch to add the isomers:
*  |  +----------------------------------------------------------------*
*  |  |
         IF ( LPATCH ) THEN
            LISMRS = .TRUE.
            IAMXMM = 2 * IZRHGH (NRN) + IMRHGH (NRN) + NMZMNM
            NISMRS (NRN) = ISMHLP (IAMXMM)
            KI0    = K1 + 1
            KI1    = K1 + NISMRS (NRN)
            KLAST  = KI1
            KMAX   = MAX ( KMAX, KI1 )
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            DO IM = 1, IMRHGH (NRN)
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |
               DO IZ = 1, IZRHGH (NRN)
                  IN  = IM + IZ + NMZMNM
                  IA  = IN + IZ
                  IZN = IZ + IZRHGH (NRN) * ( IM - 1 ) + KBURSN (NRN)
*  |  |  |  |  +-------------------------------------------------------*
*  |  |  |  |  |  No unstable nucleus for A=<2
                  IF ( IA .GT. 2 .AND. IZ .LE. IA .AND. GMSTOR (IZN)
     &                 .GT. 0.E+00 ) THEN
                     ISOMER = -1
                     T12CUR = TIM1O2 ( IA, IZ, KKA, KKZ, T12DDA, BRDDAU,
     &                                 ISOMER )
*  |  |  |  |  |  avoid to "boost" 207m-Pb production, since it
*  |  |  |  |  |  is mostly produced by thermal capture on the ground
*  |  |  |  |  |  state (??? to be checked)
                     IF ( IA .EQ. 207 .AND. IZ .EQ. 82 ) ISOMER = 0
                     ISOMER = ABS (ISOMER)
                     DO IS = 1, ISOMER
                        KS = 0
                        CALL ISMRCH ( IA, IZ, IS, KS, ADUMMY, BDUMMY,
     &                                JDUMMY, IDUMMY )
                        IF ( KS .LE. 0 ) WRITE(*,*)' IA,IZ,IS,KS',
     &                                               IA,IZ,IS,KS
                        IF ( KS .LE. 0 ) STOP ' STOP:KS-ISOMER'
                        JZN = KI0 - 1 + KS
                        IF ( JZN .GT. KLAST ) STOP ' STOP:JZN-KLAST'
                        GMSTOR (JZN) = GMSTOR (IZN) / DBLE (ISOMER+1)
                     END DO
                     GMSTOR (IZN) = GMSTOR (IZN) / DBLE (ISOMER+1)
                  END IF
*  |  |  |  |  |
*  |  |  |  |  +-------------------------------------------------------*
               END DO
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
            END DO
*  |  |  |
*  |  |  +-------------------------------------------------------------*
**sr modified from stand-alone postprocessor
C           WRITE (*,*)' *** Isomer patch completed ***',NISMRS(NRN),NRN
**
*           WRITE (*,*)' *** K0,K1,KI0,KI1',K0,K1,KI0,KI1,' ***'
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |  With the new decay channels going on isomeric states we need
*  |  |  the isomer storage allocation anyway:
         ELSE IF ( .NOT. LISMRS ) THEN
            LISMRS = .TRUE.
            IAMXMM = 2 * IZRHGH (NRN) + IMRHGH (NRN) + NMZMNM
            NISMRS (NRN) = ISMHLP (IAMXMM)
            KI0    = K1 + 1
            KI1    = K1 + NISMRS (NRN)
            KLAST  = KI1
            KMAX   = MAX ( KMAX, KI1 )
            DO KS = KI0, KI1
               GMSTOR (KS) = 0.E+00
            END DO
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
 6535    CONTINUE
         DO J = K0, KI1
            GDSTOR (J) = GDSTOR (J) + DBLE (GMSTOR(J))    * DBLE (NCASE)
            GBSTOR (J) = GBSTOR (J) + DBLE (GMSTOR(J))**2 * DBLE (NCASE)
         END DO
         KMAX = MAX ( KMAX, KI1 )
         IF ( LFIEND ) GO TO 1000
 6600 CONTINUE
*  |
*  +-------------------------------------------------------------------*
 1000 CONTINUE
      IX = NRNMX
      IF ( NRNDC .EQ. 0 ) NRNDC = NRNMX
**sr modified from stand-alone postprocessor
C     WRITE (*,*) ' Nrnmx:', NRNMX
**
*  +-------------------------------------------------------------------*
*  |  First run: allocate the room for A-mass yields:
      IF ( MBATCH .EQ. 1 ) THEN
         DO IIX = 1, IX
            NRN = JX (IIX)
            KAURSN (NRN) = KLAST
            KLAST = KLAST + IMRHGH (NRN) + 2 * IZRHGH (NRN) + NMZMNM
            KZURSN (NRN) = KLAST
            KLAST = KLAST + IZRHGH (NRN)
            KMAX  = MAX ( KMAX, KLAST )
         END DO
         FILE1 = FILE
         LQOLD = LQ
      END IF
*  |
*  +-------------------------------------------------------------------*
      CLOSE (UNIT=1)
      KLAST = 0
*  +-------------------------------------------------------------------*
*  |  Loop over the scorings:
      DO 7000 IIX = 1, IX
         NRN    = JX (IIX)
         ACCUMP = ZERZER
*  |  +----------------------------------------------------------------*
*  |  |
         DO IG  = 1, IMRHGH (NRN) + 2 * IZRHGH (NRN) + NMZMNM
            ACCUMA (IG) = ZERZER
            ACCUMZ (IG) = ZERZER
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
         AFACT = DBLE (NCASE)
         K0    = KLAST + 1
         K1    = IMRHGH (NRN) * IZRHGH (NRN) + K0 - 1
         KLAST = K1
*  |  +----------------------------------------------------------------*
*  |  |
         DO 6540 IM = 1, IMRHGH (NRN)
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            DO 6520 IZ = 1, IZRHGH (NRN)
               IN     = IM + IZ + NMZMNM
               IA     = IN + IZ
               IZN    = IZ + IZRHGH (NRN) * ( IM - 1 ) + KBURSN (NRN)
               KK0    = IZN
               KMAX   = MAX ( KMAX, KK0 )
               DSCORE = DBLE ( GMSTOR (IZN) )
               IF ( IA .LT. 1 ) GO TO 6520
               ACCUMP = ACCUMP + DSCORE
               ACCUMA (IA) = ACCUMA (IA) + DSCORE
               ACCUMZ (IZ) = ACCUMZ (IZ) + DSCORE
 6520       CONTINUE
*  |  |  |
*  |  |  +-------------------------------------------------------------*
 6540    CONTINUE
*  |  |
*  |  +----------------------------------------------------------------*
         KI0   = K1 + 1
         K1I   = K1 + NISMRS (NRN)
         KLAST = KI1
*  |  +----------------------------------------------------------------*
*  |  |  Loop over possible isomers:
         DO 6550 KS = 1, NISMRS (NRN)
            IZN = KS - 1 + KI0
            IF ( GMSTOR (IZN) .GT. 0.E+00 ) THEN
               CALL ISMRCH ( IA, IZ, IS, KS, ADUMMY, BDUMMY, JDUMMY,
     &                       IDUMMY )
               DSCORE = DBLE ( GMSTOR (IZN) )
               IF ( IA .LT. 1 ) GO TO 6550
               ACCUMP = ACCUMP + DSCORE
               ACCUMA (IA) = ACCUMA (IA) + DSCORE
               ACCUMZ (IZ) = ACCUMZ (IZ) + DSCORE
            END IF
 6550    CONTINUE
*  |  |
*  |  +----------------------------------------------------------------*
         TOTTOT (NRN) = TOTTOT (NRN) + ACCUMP * AFACT
         TOTERR (NRN) = TOTERR (NRN) + ACCUMP * ACCUMP * AFACT
*  |  +----------------------------------------------------------------*
*  |  |
         DO IA = 1, IMRHGH (NRN) + 2 * IZRHGH (NRN) + NMZMNM
            KK0 = KAURSN (NRN) + IA
            GDSTOR (KK0) = GDSTOR (KK0) + ACCUMA (IA) * AFACT
            GBSTOR (KK0) = GBSTOR (KK0) + ACCUMA (IA) * ACCUMA (IA)
     &                   * AFACT
            KMAX = MAX ( KMAX, KK0 )
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
*  |  +----------------------------------------------------------------*
*  |  |
         DO IZ = 1, IZRHGH (NRN)
            KK0 = KZURSN (NRN) + IZ
            GDSTOR (KK0) = GDSTOR (KK0) + ACCUMZ (IZ) * AFACT
            GBSTOR (KK0) = GBSTOR (KK0) + ACCUMZ (IZ) * ACCUMZ (IZ)
     &                   * AFACT
            KMAX = MAX ( KMAX, KK0 )
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
         KMAX = MAX ( KMAX, KK0 )
 7000 CONTINUE
*  |
*  +-------------------------------------------------------------------*
**sr bug patch
      DO J = K0, KI1
         GMSTOR(J) = ZERZER
      END DO
**
      GO TO 1111
 2000 CONTINUE
*  +-------------------------------------------------------------------*
*  |
      DO J = 1, KMAX
         GDSTOR (J) = GDSTOR (J) / DBLE (NCTOT)
         GBSTOR (J) = GBSTOR (J) / DBLE (NCTOT)
      END DO
*  |
*  +-------------------------------------------------------------------*
*  +-------------------------------------------------------------------*
*  |
      DO IIX = 1, IX
         NRN = JX (IIX)
         TOTTOT (NRN) = TOTTOT (NRN) / DBLE (NCTOT)
         TOTERR (NRN) = TOTERR (NRN) / DBLE (NCTOT)
      END DO
*  |
*  +-------------------------------------------------------------------*
*  +-------------------------------------------------------------------*
*  |  Real sum of Mbatch runs:
      IF ( MBATCH - 1 .GT. 0 ) THEN
         DO J = 1, KMAX
            GBSTOR(J)=SQRT(GBSTOR(J))
            GBSTOR(J)=(GBSTOR(J)-GDSTOR(J))*(GBSTOR(J)+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
            NRN = JX (IIX)
            TOTERR(NRN)=SQRT(TOTERR(NRN))
            TOTERR(NRN)=(TOTERR(NRN)-TOTTOT(NRN))
     &                 *(TOTERR(NRN)+TOTTOT(NRN))/DBLE(MBATCH-1)
            IF (ABS(TOTERR(NRN)).LT.1.D-12*TOTTOT(NRN)**2)
     &              TOTERR(NRN)=ZERZER
            IF (TOTTOT(NRN) .GT. ZERZER) THEN
               TOTERR(NRN)=MIN(SQRT(TOTERR(NRN))/TOTTOT(NRN),0.99D+00)
            END IF
         END DO
         LSTATI = .FALSE.
*  |
*  +-------------------------------------------------------------------*
*  |  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)
         IF ( LEVOLN ) READ (11)
         KLAST = 0
         DO JJX = 1, IX / NRNDC
            DO JIX = 1, NRNDC
               READ (11)
               IF ( LEVOLN ) READ (11)
               READ (11)
               IF ( LISMRS .AND. .NOT. LPATCH ) THEN
                  READ (11)
                  READ (11)
               END IF
            END DO
            READ (11) CHSTAT,IJX
            IF ( CHSTAT .NE. 'STATISTICS' .OR. IJX .NE. NRNDC )
     &         STOP 'STATI'
            DO 7400 JIX = 1, NRNDC
               IIX = (JJX-1)*NRNDC + JIX
               NRN = JX (IIX)
               READ (11) TOTA, TOTB
               IF ( ABS (TOTA-TOTTOT(NRN)) .GT. 1.D-05*TOTA )
     &            WRITE (*,*)' *** Error Tota:',TOTA,TOTTOT(NRN)
               TOTERR (NRN) = TOTB
               KK0 = 1
               KK1 = IMRHGH (NRN) + 2 * IZRHGH (NRN) + NMZMNM
               READ (11) ( RCCUMA (J), J = KK0, KK1 )
               IF ( .NOT. LEVOLN ) THEN
                  DO 7100 J = KK0, KK1
                     IF ( ABS (RCCUMA(J)-GDSTOR(J+KAURSN(NRN))) .GT.
     &                    1.D-05*RCCUMA(J) )
     &               WRITE (*,*)' *** Error Gdstor:',RCCUMA(J),
     &                            GDSTOR(J+KAURSN(NRN))
 7100             CONTINUE
               END IF
               READ (11) ( RCCUMA (J), J = KK0, KK1 )
               DO 7150 J = KK0, KK1
                  GBSTOR (J+KAURSN(NRN)) = RCCUMA (J)
 7150          CONTINUE
               KK0 = 1
               KK1 = IZRHGH (NRN)
               READ (11) ( RCCUMZ (J), J = KK0, KK1 )
               IF ( .NOT. LEVOLN ) THEN
                  DO 7200 J = KK0, KK1
                     IF ( ABS (RCCUMZ(J)-GDSTOR(J+KZURSN(NRN))) .GT.
     &                    1.D-05*RCCUMZ(J) )
     &               WRITE (*,*)' *** Error Gdstor:',RCCUMZ(J),
     &                            GDSTOR(J+KZURSN(NRN))
 7200             CONTINUE
               END IF
               READ (11) ( RCCUMZ (J), J = KK0, KK1 )
               DO 7250 J = KK0, KK1
                  GBSTOR (J+KZURSN(NRN)) = RCCUMZ (J)
 7250          CONTINUE
               K0    = KLAST + 1
               K1    = IMRHGH (NRN) * IZRHGH (NRN) + K0 - 1
               KLAST = K1
               READ (11) ( GMSTOR(J), J = K0, K1 )
               DO 7300 J = K0, K1
                  GBSTOR (J) = DBLE (GMSTOR(J))
 7300          CONTINUE
               IF ( LISMRS .AND. .NOT. LPATCH ) THEN
                  KI0   = KLAST + 1
                  KI1   = KLAST + NISMRS (NRN)
                  KLAST = KI1
                  IF ( NISMRS (NRN) .GT. 0 ) THEN
                     READ (11) ( GMSTOR(J), J = KI0, KI1 )
                  ELSE
                     READ (11)
                  END IF
                  DO 7350 J = KI0, KI1
                     GBSTOR (J) = DBLE (GMSTOR(J))
 7350             CONTINUE
               ELSE IF ( LISMRS ) THEN
                  KI0   = KLAST + 1
                  KI1   = KLAST + NISMRS (NRN)
                  KLAST = KI1
                  DO KS = 1, NISMRS (NRN)
                     J = KI0 + KS - 1
                     CALL ISMRCH ( IA, IZ, IS, KS, ADUMMY, BDUMMY,
     &                             JDUMMY, IDUMMY )
                     IM = IA - 2 * IZ - NMZMNM
                     IF ( IA .GE. 1 ) THEN
                        IZN = K0 - 1 + ( IM - 1 ) * IZRHGH (NRN) + IZ
                        GBSTOR (J) = DBLE (GBSTOR(IZN))
                     END IF
                  END DO
               END IF
 7400       CONTINUE
         END DO
         KLAST = 0
         CLOSE ( UNIT=11 )
*  |
*  +-------------------------------------------------------------------*
*  |  Only one plain run, no possibility of making statistics:
      ELSE
         LSTATI = .FALSE.
         DO J =1, KMAX
            IF ( GDSTOR (J) .GT. ZERZER ) THEN
               GBSTOR (J) = 0.99D+00
            END IF
         END DO
         DO IIX = 1, IX
            NRN = JX (IIX)
            IF ( TOTTOT (NRN) .GT. ZERZER ) THEN
               TOTERR (NRN) = 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 ) THEN
         IF ( LEVOLN ) THEN
            RUNTIM='   ***** Time Evolution *****'
         ELSE
            RUNTIM='      ***** Sum file *****'
         END IF
**sr modified from stand-alone postprocessor
C        WRITE(*,*)' Type the output file name:'
C        READ (*,'(A)') FILE
C        LQ = LNNBLN (FILE)
**
* Start_VAX_seq
*        OPEN (UNIT=1,FILE=FILE,STATUS='NEW',FORM='UNFORMATTED')
* End_VAX_seq
**sr commented
* Start_UNIX_seq
c        OPEN (UNIT=1,FILE=FILE,STATUS='UNKNOWN',FORM='UNFORMATTED')
* End_UNIX_seq
c        IF ( LEVOLN ) THEN
c           WRITE (1) RUNTIT, RUNTIM, SNGL (WCTOT),-NCTOT
c           WRITE (1) NIRRDT, ( SNGL (BIRRDT(IR)), SNGL (TIRRDT(IR)),
c    &                          IR = 1, NIRRDT )
c        ELSE
c           WRITE (1) RUNTIT, RUNTIM, SNGL (WCTOT), NCTOT
c        END IF
**
*  |
*  +-------------------------------------------------------------------*
*  |
      ELSE
         FILE = FILE1
         LQ   = LQOLD
      END IF
*  |
*  +-------------------------------------------------------------------*
      KLAST  = 0
*  +-------------------------------------------------------------------*
*  |  Loop over scorings:
      DO 6800 JJ=1,IX
         NRN   = JX(JJ)
**sr commented
c        IF ( .NOT. LSTATI )
c    &   WRITE (1)    NRN, TIURSN (NRN), ITURSN (NRN), NRURSN (NRN),
c    &          SNGL (VURSNC (NRN)), IMRHGH(NRN), IZRHGH(NRN), NMZMNM
c        IF ( LEVOLN .AND. .NOT. LSTATI ) WRITE (1) TDECAY (NRN)
**
         K0    = KLAST  + 1
         K1    = IMRHGH (NRN) * IZRHGH (NRN) + K0 - 1
         KLAST = K1
**sr commented
c        IF ( .NOT. LSTATI )
c    &   WRITE (1) (SNGL(GDSTOR(J)), J = K0, K1)
**
         IF ( .NOT. LSTATI .AND. LISMRS ) THEN
**sr commented
c           WRITE (1) 'ISOMERS:  ', NISMRS (NRN)
**
            KI0   = KLAST + 1
            KI1   = KLAST + NISMRS (NRN)
            KLAST = KI1
**sr commented
c           IF ( KI1 .GE. KI0 ) THEN
c              WRITE (1) (SNGL(GDSTOR(J)), J = KI0, KI1)
c           ELSE
c              WRITE (1)
c           END IF
**
         END IF
 6800 CONTINUE
*  |
*  +-------------------------------------------------------------------*
*  Now write the statistical data, first a record to flag them
**sr commented
c     IF ( .NOT. LSTATI ) WRITE (1) 'STATISTICS',IX
**
*  Then the errors concerning the binning data
      KLAST  = 0
      INDX   = INDLST ( FILE, '.' )
      IF ( INDX .GT. 0 ) THEN
* Start_VAX_seq
*        OPEN (UNIT=3,FILE=FILE(1:INDX-1)//'_SUM.LIS',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
* End_VAX_seq
* Start_UNIX_seq
**sr commented
c        OPEN (UNIT=3,FILE=FILE(1:INDX-1)//'_sum.lis',STATUS='UNKNOWN')
**
* End_UNIX_seq
      ELSE
* Start_VAX_seq
*        OPEN (UNIT=3,FILE=FILE(1:LQ)//'_SUM.LIS',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
* End_VAX_seq
* Start_UNIX_seq
**sr commented
c        OPEN (UNIT=3,FILE=FILE(1:LQ)//'_sum.lis',STATUS='UNKNOWN')
**
* End_UNIX_seq
      END IF
**sr commented
c     WRITE (3,*)' **** ',FILE(1:LQ),' ****'
c     WRITE (3,*)
c     WRITE (3,*)
c     WRITE (3,*)'   Total primaries run:',NCTOT
c     WRITE (3,*)
c     WRITE (3,*)'   Total weight of the primaries run:',SNGL(WCTOT)
c     WRITE (3,*)
**
*  +-------------------------------------------------------------------*
*  |  Loop over scorings:
      DO 6850 JJ=1,IX
         NRN = JX(JJ)
**sr commented
c        WRITE (1) SNGL(TOTTOT(NRN)),SNGL(TOTERR(NRN))
c        IF ( LEVOLN ) THEN
c           WRITE (3,*)
c           WRITE (3,*) ' Beam intensity  :',BEAINT,' pr/s'
c           WRITE (3,*)
c           WRITE (3,*) ' Irradiation time:',TIMEIR,' s'
c           WRITE (3,*)
c           WRITE (3,*) ' Decay time      :',TDECAY(NRN),' s'
c        END IF
c        WRITE (3,*)
c        WRITE (3,*) ' Detector n: ', NRN,'  ',TIURSN (NRN)
c        WRITE (3,*) '    (Region',NRURSN(NRN),' Volume:  ',
c    &                     SNGL(VURSNC(NRN)),' cmc,'
c        WRITE (3,*) '     distr. type  :',ITURSN(NRN),'   ,'
c        WRITE (3,*) '     Z_max:',IZRHGH(NRN),', N-Z_max:',
c    &                     IMRHGH(NRN)+NMZMNM,', N-Z_min:',NMZMNM+1,
c    &               ')'
c        WRITE (3,*)
c        WRITE (3,*)'    Tot. response (n/cmc/pr)', SNGL(TOTTOT(NRN)),
c    &              '  +/-',SNGL(TOTERR(NRN))*100.0,' %'
c        WRITE (3,*)'    ( --> Nuclei/pr         ',
c    &                      SNGL(TOTTOT(NRN)*VURSNC(NRN)),
c    &              '  +/-',SNGL(TOTERR(NRN))*100.0,' % )'
c        WRITE (3,*)
c        WRITE (3,*)
c    &         ' **** Isotope Yield as a function of Mass Number ****',
c    &         ' ****            (nuclei / cmc / pr)             ****'
c        WRITE (3,*)
c        WRITE (3,*)
c        WRITE (3,*) ' A_min:',2,' -  A_max:',IMRHGH(NRN)+2*IZRHGH(NRN)
c    &               + NMZMNM
c        WRITE (3,*)
**
         KK0 = KAURSN (NRN) + 1
         KK1 = KAURSN (NRN) + IMRHGH (NRN) + 2 * IZRHGH (NRN) + NMZMNM
         DO J = KK1, KK0, -1
            IA = J - KK0 + 1
**sr commented
c           IF ( GDSTOR(J) .GT. AZRZRZ .AND. IA .GE. 1 )
c    &      WRITE (3,*) ' A: ', IA, SNGL (GDSTOR(J)),' +/- ',
c    &                  SNGL (GBSTOR(J)*100.D+00),' % '
**
         END DO
**sr commented
c        IF ( .NOT. LSTATI ) THEN
c           WRITE (1) ( SNGL (GDSTOR(J)), J = KK0, KK1 )
c           WRITE (1) ( SNGL (GBSTOR(J)), J = KK0, KK1 )
c        END IF
c        WRITE (3,*)
c        WRITE (3,*)
c    &       ' **** Isotope Yield as a function of Atomic Number ****',
c    &       ' ****             (nuclei / cmc / pr)              ****'
c        WRITE (3,*)
c        WRITE (3,*)
c        WRITE (3,*) ' Z_min:',1,' -  Z_max:',IZRHGH(NRN)
c        WRITE (3,*)
**
         KK0 = KZURSN (NRN) + 1
         KK1 = KZURSN (NRN) + IZRHGH (NRN)
         DO J = KK1, KK0, -1
            IZ = J - KK0 + 1
**sr commented
c           IF ( GDSTOR(J) .GT. AZRZRZ )
c    &      WRITE (3,*) ' Z: ', IZ, SNGL (GDSTOR(J)),' +/- ',
c    &                  SNGL (GBSTOR(J)*100.D+00),' % '
**
         END DO
         IF ( .NOT. LSTATI ) THEN
**sr commented
c           WRITE (1) ( SNGL (GDSTOR(J)), J = KK0, KK1 )
c           WRITE (1) ( SNGL (GBSTOR(J)), J = KK0, KK1 )
**
         END IF
         K0    = KLAST + 1
         K1    = IMRHGH (NRN) * IZRHGH (NRN) + K0 - 1
         KLAST = K1
**sr commented
c        IF ( .NOT. LSTATI )
c    &   WRITE (1) (SNGL(GBSTOR(J)), J = K0, K1)
c        WRITE (3,*)
c        WRITE (3,*)
c        WRITE (3,*)' Residual nuclei distribution:'
c        WRITE (3,*)
c    &         ' **** Residual nuclei distribution  ****',
c    &         ' ****    (nuclei / cmc / pr)        ****'
c        WRITE (3,*)
c        WRITE (3,*)
**
         DO IA = 1, 260
            DO IZ = 1, 100
               RESNUC (IZ,IA) = ZERZER
               RESERR (IZ,IA) = ZERZER
**sr
               DO IDET=1,MXDET
                  ADDNUC (IDET,IZ,IA) = 0.E+00
                  ADDERR (IDET,IZ,IA) = 0.E+00
                  ADDNUC (IDET,IZ,IA+260) = 0.E+00
                  ADDERR (IDET,IZ,IA+260) = 0.E+00
               ENDDO
**
               IABEG  (IZ) = 500
               IAEND  (IZ) = 0
            END DO
         END DO
         DO IM = 1, IMRHGH (NRN)
            DO IZ  = 1, IZRHGH (NRN)
               IN  = IM + IZ + NMZMNM
               IA  = IN + IZ
               IF ( IA .GE. 1 ) THEN
                  IZN = K0 - 1 + ( IM - 1 ) * IZRHGH (NRN) + IZ
                  RESNUC (IZ,IA) = GDSTOR (IZN)
                  IF ( ABS (GDSTOR (IZN)) .GT. 0.E+00 ) THEN
                     IABEG (IZ) = MIN ( IABEG (IZ), IA )
                     IAEND (IZ) = MAX ( IAEND (IZ), IA )
                     RESERR (IZ,IA) = GBSTOR (IZN) * 100.D+00
                  END IF
               END IF
            END DO
         END DO
         IASTRT = IMRHGH (NRN) + 2 * IZRHGH (NRN) + NMZMNM
         IZ0    = IZRHGH (NRN) + 1
 9220    CONTINUE
         IZ1 = IZ0 - 1
         IZ0 = MAX ( 1, IZ1 - 10 )
*        IAS1= IASTLN (MIN(MAX(IZ1,1),100))
*        IAS0= IASTLN (MIN(MAX(IZ0,1),100))
*        IA1 = MIN ( IASTRT, IAS1 + MAX ( MIN (IZ1/5,10),  5 ) )
*        IA0 = MAX ( IZ0   , IAS0 - MAX ( MIN (IZ1/3,20), 10 ) )
*        IA1 = MIN ( IA1, IMRHGH (NRN) + 2 * IZ1 + NMZMNM )
*        IA0 = MAX ( IA0, 1 + 2 * IZ0 + NMZMNM )
*        IA0 = MAX ( IA0, 1 )
         IA1 = 0
         IA0 = 500
         DO IZ = IZ0, IZ1
            IA1 = MAX ( IA1, IAEND (IZ) )
            IA0 = MIN ( IA0, IABEG (IZ) )
         END DO
**sr commented
c        WRITE (3,9240) ( IZ, IZ = IZ0, IZ1 )
**
         DO 9230 IA = IA1, IA0, - 1
            LCHECK = .FALSE.
            DO IZ = IZ0, IZ1
               LCHECK = LCHECK .OR. RESNUC (IZ,IA) .GE. AZRZRZ
            END DO
            IF ( LCHECK ) THEN
**sr commented
c           WRITE (3,9250) IA, ( RESNUC (IZ,IA), IZ = IZ0, IZ1 )
c           WRITE (3,9260)     ( CHERR, RESERR (IZ,IA), IZ = IZ0, IZ1 )
**
            END IF
 9230    CONTINUE
         IF ( IZ0 .GT. 1 ) GO TO 9220
**sr
C        CALL ANALYS(RESNUC,RESERR,VURSNC(NRN),NRN,0)
**
**sr commented
c        WRITE (3,*)
**
         IF ( LISMRS ) THEN
            KI0   = KLAST + 1
            KI1   = KLAST + NISMRS (NRN)
            KLAST = KI1
            IF ( .NOT. LSTATI ) THEN
               IF ( NISMRS (NRN) .GT. 0 ) THEN
**sr commented
c                 WRITE (1) (SNGL(GBSTOR(J)),J=KI0,KI1)
**
               ELSE
**sr commented
c                 WRITE (1)
**
               END IF
            END IF
            IF ( NISMRS (NRN) .GT. 0 ) THEN
**sr commented
c              WRITE (3,*) ' **** Isomers (nuclei / cmc / pr)     ****'
c              WRITE (3,*) '     A        Z       mth'
**
            END IF
            DO KS = 1, NISMRS (NRN)
               IZN = KI0 + KS - 1
               IF ( GDSTOR(IZN) .GT. AZRZRZ ) THEN
                  CALL ISMRCH ( IA, IZ, IS, KS, ADUMMY, BDUMMY, JDUMMY,
     &                          IDUMMY )
**sr commented
c                 WRITE (3,9270) IA, IZ, IS, SNGL(GDSTOR(IZN)),
c    &                               SNGL(GBSTOR(IZN)*100.D+0)
**
               END IF
            END DO
**sr commented
c           WRITE (3,*)
**
         END IF
 6850 CONTINUE
*  |
*  +-------------------------------------------------------------------*
      IF ( .NOT. LSTATI ) CLOSE ( UNIT = 1 )
**sr commented
c     CLOSE ( UNIT = 3 )
**
**sr
**sr modified from stand-alone postprocessor
C     WRITE (*,*)
C    &  ' Do you want to calculate an irradiation cycle (def=yes)?'
C     READ  (*,'(A)') YES
      YES = 'y'
**
      LCYCLE = YES .NE. 'N' .AND. YES .NE. 'n'
      ICOOL  = 0
      IF (LCYCLE) THEN
**sr modified from stand-alone postprocessor
C        OPEN(50,FILE='irrcyc.out',STATUS='UNKNOWN')
C        OPEN(51,FILE='irrcyc.lst',STATUS='UNKNOWN')
C        WRITE(50,3001) CRUNTT(1:59),CRUNTM(8:16),CRUNTM(24:32),
C    &                  NPRIM,TOTWGT
**
 3001    FORMAT(/,1X,'Title of the run : ',A,/,1X,'Date and Time    : ',
     &          2A,/,/,1X,'Number of primaries       : ',I7,/,
     &          1X,'Total weight of primaries : ',F10.2,/)
**sr modified from stand-alone postprocessor
C        WRITE(50,3002)
C3002    FORMAT(1X,'FLUKA output read from file(s) : ')
C        DO 211 II=1,NFILE
C           WRITE(50,'(1X,A)') CINFLE(II)
C 211    CONTINUE
C        CALL IRRCYC
         CALL IRRCYC(TDCY)
**
         TIMEIR = TIRRAD
         FACMUL = FACMUL*PRTSEC
**sr added. to be checked!!!
         TIMEDC = 0.D+00
**
      ELSE
         TIMEDC = 0.D+00
**sr modified from stand-alone postprocessor
C        WRITE (*,*)' Irradiation time (s):'
C        READ  (*,*)  TIMEIR
**
      ENDIF
  210 CONTINUE
      ICOOL = ICOOL+1
      IRRAD = 0
      DO IA = 1, 260
         DO IZ = 1, 100
            DO IDET=1,MXDET
               ADDNUC (IDET,IZ,IA) = 0.E+00
               ADDERR (IDET,IZ,IA) = 0.E+00
               ADDNUC (IDET,IZ,IA+260) = 0.E+00
               ADDERR (IDET,IZ,IA+260) = 0.E+00
            ENDDO
         END DO
      END DO
c     TIMEDC = 0.D+00
c     WRITE (*,*)' Irradiation time (s):'
c     READ  (*,*)  TIMEIR
**
***** here start irradiation build up
**sr commented
c     TIMEDC = ZERZER
c     WRITE (*,*)' Irradiation time (s):'
c     READ  (*,*)  TIMEIR
**
      NIRRDT = 1
      NTRDCY = 1
      TIRRDT   (0)    = ZERZER
      CIRRDT   (0)    = ZERZER
      TIRRDT (NIRRDT) = TIMEIR
      BIRRDT (NIRRDT) = FACMUL
      CIRRDT (NIRRDT) = BIRRDT (NIRRDT) * TIRRDT (NIRRDT)
      TIRDCY (NTRDCY) = TIMEDC
      LEVOLN = .TRUE.
      BEAINT = FACMUL
      RUNTIM ='   ***** Time Evolution *****'
**sr
      IF ((.NOT.LCYCLE).OR.(LCYCLE.AND.(ICOOL.EQ.1))) THEN
**
**sr modified from stand-alone postprocessor
C     WRITE(*,*)' Type the evolution output file name:'
C     READ (*,'(A)') FILE
**
      LQ     = LNNBLN (FILE)
* Start_VAX_seq
*     OPEN (UNIT=1,FILE=FILE,STATUS='NEW',FORM='UNFORMATTED')
* End_VAX_seq
* Start_UNIX_seq
**sr commented
c     OPEN (UNIT=1,FILE=FILE,STATUS='UNKNOWN',FORM='UNFORMATTED')
**
* End_UNIX_seq
**sr commented
c     WRITE (1) RUNTIT, RUNTIM, SNGL(WCTOT), -NCTOT
c     WRITE (1) NIRRDT, ( SNGL (BIRRDT(IR)), SNGL (TIRRDT(IR)),
c    &                    IR = 1,NIRRDT )
**
      INDX   = INDEX ( FILE,'.')
      IF ( INDX .GT. 0 ) THEN
* Start_VAX_seq
*        OPEN (UNIT=3,FILE=FILE(1:INDX-1)//'_SUM.LIS',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
*        OPEN (UNIT=66,FILE=FILE(1:INDX-1)//'_RES.LIS',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
*        OPEN (UNIT=67,FILE=FILE(1:INDX-1)//'_TOT.LIS',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
* End_VAX_seq
* Start_UNIX_seq
**sr commented
c        OPEN (UNIT=3,FILE=FILE(1:INDX-1)//'_sum.lis',STATUS='UNKNOWN')
c        OPEN (UNIT=66,FILE=FILE(1:INDX-1)//'_res.lis',STATUS='UNKNOWN')
c        OPEN (UNIT=67,FILE=FILE(1:INDX-1)//'_tot.lis',STATUS='UNKNOWN')
**
* End_UNIX_seq
      ELSE
* Start_VAX_seq
*        OPEN (UNIT=3,FILE=FILE(1:LQ)//'_SUM.LIS',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
*        OPEN (UNIT=66,FILE=FILE(1:LQ)//'_RES.LIS',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
*        OPEN (UNIT=67,FILE=FILE(1:LQ)//'_TOT.LIS',STATUS='NEW',
*    &         CARRIAGECONTROL='LIST')
* End_VAX_seq
* Start_UNIX_seq
**sr commented
c        OPEN (UNIT=3,FILE=FILE(1:LQ)//'_sum.lis',STATUS='UNKNOWN')
c        OPEN (UNIT=66,FILE=FILE(1:LQ)//'_res.lis',STATUS='UNKNOWN')
c        OPEN (UNIT=67,FILE=FILE(1:LQ)//'_tot.lis',STATUS='UNKNOWN')
**
* End_UNIX_seq
      END IF
**sr
      ENDIF
**
**sr commented
c     WRITE (3,*)' **** ',FILE(1:LQ),' ****'
c     WRITE (3,*)
c     WRITE (3,*)
c     WRITE (3,*)'   Total primaries run:',NCTOT
c     WRITE (3,*)
c     WRITE (3,*)'   Total weight of the primaries run:',SNGL(WCTOT)
c     WRITE (3,*)
**************
c     WRITE (66,*)' **** ',FILE(1:LQ),' ****'
c     WRITE (66,*)
c     WRITE (66,*)
c     WRITE (66,*)'   Total primaries run:',NCTOT
c     WRITE (66,*)
c     WRITE (66,*)'   Total weight of the primaries run:',SNGL(WCTOT)
c     WRITE (66,*)
**************
**************
c     WRITE (67,*)' **** ',FILE(1:LQ),' ****'
c     WRITE (67,*) NRNMX
**
**************
***** HERE ASK FOR DECAY TIMES
 9999 CONTINUE
**sr new version: load here the internal array of decay times
      IF (LCYCLE) THEN
         IRRAD = IRRAD+1
         IF (IRRAD.GT.NIRRAD) THEN
            TIMEDC = 0.0D0
         ELSE
            TIMEDC = TDECY(ICOOL,IRRAD)
         ENDIF
      ELSE
**sr modified from stand-alone postprocessor
C        WRITE (*,*)' Decay time (s):'
C        READ  (*,*) TIMEDC
**
      ENDIF
c     WRITE (*,*)' Decay time (s):'
c     READ  (*,*) TIMEDC
**
      IF ( TIMEDC .GT. ZERZER ) THEN 
         NTRDCY = NTRDCY + 1
         TIRDCY (NTRDCY) = TIMEDC
         GO TO 9999      
      END IF
      NXIDAU = MXACTV / MAX ( NIRRDT, 1 ) / MAX ( NTRDCY, 1 )
      NXIDAU = MIN ( NXIDAU, MXIDAU )
      MXCRIR = MAX ( NIRRDT, 1 )
      MXCRDT = MAX ( NTRDCY, 1 )
      MXCRDA = NXIDAU
      KDIM  = KMAX
      KKMAX = KMAX * NTRDCY
*  +-------------------------------------------------------------------*
*  |  Set up the starting condition at T_ir (only for radioactive
*  |  nuclei):
      if (kkmax.gt.MXDUMM) then
         write(lunout,*) 'KKMAX ',KKMAX
         write(lunout,*) 'IZEVO,IAEVO,WEI,TDCY ',IZEVO,IAEVO,WEI,TDCY
         stop
      endif
      DO J = 1, KKMAX
         GESTOR (J) = GDSTOR (J)
         GCSTOR (J) = GBSTOR (J)
         GDSTOR (J) = ZERZER
         GBSTOR (J) = ZERZER
      END DO
*  |
*  +-------------------------------------------------------------------*
      KLAST  = 0
*  +-------------------------------------------------------------------*
*  |  Loop on the residual nuclei sets:
      DO 8300 JJ = 1, IX
         NRN   = JX(JJ)
         K0    = KLAST  + 1
         K1    = IMRHGH (NRN) * IZRHGH (NRN) + K0 - 1
         KI0   = K1 + 1
         KLAST = K1
*  |  +----------------------------------------------------------------*
*  |  |  Loop on N-Z-NMZmin:
         DO 8888 IM = 1, IMRHGH (NRN)
*  |  |  +-------------------------------------------------------------*
*  |  |  |  Loop on Z:
            DO 8000 IZ  = 1, IZRHGH (NRN)
               IN  = IM + IZ + NMZMNM
               IA  = IN + IZ
*               LPRINT = IA .EQ. 141
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |  Possible isotope:
               IF ( IA .GE. 1 .AND. IZ .LE. IA ) THEN
                  IZN    = K0 - 1 + ( IM - 1 ) * IZRHGH (NRN) + IZ
                  IF ( GESTOR (IZN) .LE. AZRZRZ ) GO TO 8000
                  ISOMER = 0
                  CALL FLLACT ( IA, IZ, ISOMER, NDAUGH, IAD, IZD, 
     &                          ISD, ACTIV, BRCTOT, MXCRIR, MXCRDT, 
     &                          MXCRDA, MXIDAU )
                  IF ( MXCRIR .NE. NIRRDT .OR. MXCRDT .NE. NTRDCY .OR.
     &                 MXCRDA .NE. NXIDAU )
     &               CALL FLABRT ( 'USRSUWEV', 'MXCRIR,MXCRDT,MXCRDA' )
                  WEIISO = ONEONE
*  |  |  |  |  +-------------------------------------------------------*
*  |  |  |  |  |  Not stable:
                  IF ( NDAUGH .GT. 0 ) THEN
                     ORI  = GESTOR (IZN)
                     EORI = GCSTOR (IZN)
*  |  |  |  |  |  +----------------------------------------------------*
*  |  |  |  |  |  |  Loop on daughters 
                     DO 8050 IDAUGH = 1, NDAUGH
                        KZ  = IZD (IDAUGH)
                        KA  = IAD (IDAUGH)
                        ISM = ISD (IDAUGH)
*  |  |  |  |  |  |  Possible if of negligible weight:
                        IF ( KA .LE. 0 ) GO TO 8050
                        KN  = KA - KZ
*  |  |  |  |  |  |  +--------------------------------------------------*
*  |  |  |  |  |  |  |  Normal (ground state) isotope:
                        IF ( ISM .LE. 0 ) THEN
                           KM = KA - 2 * KZ - NMZMNM
                           IF ( KM .LE. 0 .OR. KM .GT. IMRHGH (NRN) )
     &                          GO TO 8050
                           IF ( KZ .LE. 0 .OR. KZ .GT. IZRHGH (NRN) )
     &                          GO TO 8050
                           KZN = K0 - 1 + ( KM - 1 ) * IZRHGH (NRN) + KZ
*  |  |  |  |  |  |  |
*  |  |  |  |  |  |  +--------------------------------------------------*
*  |  |  |  |  |  |  |  Isomer:
                        ELSE
                           KSM = 0
                           CALL ISMRCH ( KA, KZ, ISM, KSM, ADUMMY, 
     &                                   BDUMMY, JDUMMY, IDUMMY )
                           IF ( KSM .GT. NISMRS (NRN) ) GO TO 8050
                           KZN = KI0 - 1 + KSM
                        END IF
*  |  |  |  |  |  |  |
*  |  |  |  |  |  |  +--------------------------------------------------*
*  |  |  |  |  |  |  +--------------------------------------------------*
*  |  |  |  |  |  |  |
                        DO 8700 ITR = 1, NIRRDT
*                          ORIT = ORI / ( TIRRDT (ITR) - TIRRDT (ITR-1) )
                           ORIT = ORI * BIRRDT (ITR)
*  |  |  |  |  |  |  |  +-----------------------------------------------*
*  |  |  |  |  |  |  |  |
                           DO 8701 ITD = 1, NTRDCY
*  |  |  |  |  |  |  |  |  Activ = Activ(Nirrdt,Ntrdcy,Nxidau)
                              IRIDIA = ITR + NIRRDT * ( ITD - 1 )
     &                               + NIRRDT * NTRDCY * ( IDAUGH - 1 )
                              DCYFCT = ACTIV ( IRIDIA )
                              KIND   = KZN + (ITD-1) * KDIM
                              GDSTOR (KIND) = GDSTOR (KIND) + WEIISO
     &                                      * ORIT * DCYFCT
                              GBSTOR (KIND) = GBSTOR (KIND) + ( WEIISO
     &                                      * EORI * ORIT * DCYFCT )**2
 8701                      CONTINUE
*  |  |  |  |  |  |  |  |
*  |  |  |  |  |  |  |  +-----------------------------------------------*
 8700                   CONTINUE
*  |  |  |  |  |  |  |
*  |  |  |  |  |  |  +--------------------------------------------------*
 8050                CONTINUE
*  |  |  |  |  |  |
*  |  |  |  |  |  +----------------------------------------------------*
                  END IF
*  |  |  |  |  |  
*  |  |  |  |  +-------------------------------------------------------*
               END IF
*  |  |  |  | 
*  |  |  |  +----------------------------------------------------------*
 8000       CONTINUE 
*  |  |  |
*  |  |  +-------------------------------------------------------------*
 8888    CONTINUE
*  |  |
*  |  +----------------------------------------------------------------*
***********  The same for isomers
         KI0   = K1  + 1
         KI1   = KI0 + NISMRS (NRN) - 1
         KLAST = KI1
*  |  +----------------------------------------------------------------*
*  |  |  Isomers:
         IF ( LISMRS .AND. NISMRS (NRN) .GT. 0 ) THEN
*  |  |  +-------------------------------------------------------------*
*  |  |  |  Loop on ( possible) isomers of the current residual
*  |  |  |  nucleus:
            DO 8889 KS = 1, NISMRS (NRN)
               IDAUGH = -1000
               IZN    = KI0 + KS - 1
               IF ( GESTOR (IZN) .LE. AZRZRZ ) GO TO 8889
               WEIISO = ONEONE
               CALL ISMRCH ( IA, IZ, IS, KS, ADUMMY, BDUMMY, JDUMMY,
     &                       IDUMMY )
               JSOMER = IS
               CALL FLLACT ( IA, IZ, JSOMER, NDAUGH, IAD, IZD, ISD,
     &                       ACTIV, BRCTOT, MXCRIR, MXCRDT, MXCRDA, 
     &                       MXIDAU )
               IF ( MXCRIR .NE. NIRRDT .OR. MXCRDT .NE. NTRDCY .OR.
     &              MXCRDA .NE. NXIDAU )
     &            CALL FLABRT ( 'USRSUWEV', 'MXCRIR,MXCRDT,MXCRDA' )
*  |  |  |  +----------------------------------------------------------*
*  |  |  |  |  Not stable:
               IF ( NDAUGH .GT. 0 ) THEN
                  ORI  = GESTOR (IZN)
                  EORI = GCSTOR (IZN)
*  |  |  |  |  +-------------------------------------------------------*
*  |  |  |  |  |  Loop on daughters 
                  DO 8250 IDAUGH = 1, NDAUGH
                     KZ  = IZD (IDAUGH)
                     KA  = IAD (IDAUGH)
                     ISM = ISD (IDAUGH)
*  |  |  |  |  |  Possible if of negligible weight:
                     IF ( KA .LE. 0 ) GO TO 8250
                     KN  = KA - KZ
*  |  |  |  |  |  +----------------------------------------------------*
*  |  |  |  |  |  |  Normal (ground state) isotope:
                     IF ( ISM .LE. 0 ) THEN
                        KM  = KA - 2 * KZ - NMZMNM
                        IF ( KM .LE. 0 .OR. KM .GT. IMRHGH (NRN) )
     &                       GO TO 8250
                        IF ( KZ .LE. 0 .OR. KZ .GT. IZRHGH (NRN) )
     &                       GO TO 8250
                        KZN = K0 - 1 + ( KM - 1 ) * IZRHGH (NRN) + KZ
*  |  |  |  |  |  |
*  |  |  |  |  |  +----------------------------------------------------*
*  |  |  |  |  |  |  Isomer:
                     ELSE
                        KSM = 0
                        CALL ISMRCH ( KA, KZ, ISM, KSM, ADUMMY, 
     &                       BDUMMY, JDUMMY, IDUMMY )
                        IF ( KSM .GT. NISMRS (NRN) ) GO TO 8250
                        KZN = KI0 - 1 + KSM
                     END IF
*  |  |  |  |  |  |
*  |  |  |  |  |  +----------------------------------------------------*
*  |  |  |  |  |  +----------------------------------------------------*
*  |  |  |  |  |  |
                     DO 8702 ITR = 1, NIRRDT
*                       ORIT = ORI / ( TIRRDT (ITR) - TIRRDT (ITR-1) )
                        ORIT = ORI * BIRRDT (ITR)
*  |  |  |  |  |  |  +-------------------------------------------------*
*  |  |  |  |  |  |  |
                        DO 8703 ITD = 1, NTRDCY
*  |  |  |  |  |  |  |  Activ = Activ(Nirrdt,Ntrdcy,Nxidau)
                           IRIDIA = ITR + NIRRDT * ( ITD - 1 )
     &                            + NIRRDT * NTRDCY * ( IDAUGH - 1 )
                           DCYFCT = ACTIV ( IRIDIA )
                           KIND   = KZN + (ITD-1) * KDIM
                           GDSTOR (KIND) = GDSTOR (KIND) + WEIISO
     &                                   * ORIT * DCYFCT
                           GBSTOR (KIND) = GBSTOR (KIND) + ( WEIISO
     &                                   * EORI * ORIT * DCYFCT )**2
 8703                   CONTINUE
*  |  |  |  |  |  |  |
*  |  |  |  |  |  |  +-------------------------------------------------*
 8702                CONTINUE
*  |  |  |  |  |  |
*  |  |  |  |  |  +----------------------------------------------------*
*  |  |  |  |  |
*  |  |  |  |  +-------------------------------------------------------*
 8250             CONTINUE
               ELSE
                  WRITE (*,*)' ??? Isomer',IA,IZ,' with no daughters'
               END IF
*  |  |  |  |
*  |  |  |  +----------------------------------------------------------*
 8889       CONTINUE 
*  |  |  |
*  |  |  +-------------------------------------------------------------*
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
 8300 CONTINUE
*  |
*  +-------------------------------------------------------------------*
****  Now prepare for output 
*  +-------------------------------------------------------------------*
*  |  Loop on decay times
      DO 9000 ITD = 1, NTRDCY
         KOFFS = ( ITD - 1 ) * KDIM
         KLAST = KOFFS
         close (unit=99)
*  |  +----------------------------------------------------------------*
*  |  |  Loop on detectors for totals and errors
         DO 9001 JJ = 1, IX
            NRN   = JX (JJ)
            K0    = KLAST  + 1
            K1    = IMRHGH (NRN) * IZRHGH (NRN) + K0 - 1
            KLAST = K1
            IF ( LISMRS ) THEN
               KI0   = KLAST + 1
               KI1   = KLAST + NISMRS (NRN)
               KLAST = KI1
            ELSE
               KI1   = K1
            END IF
            TOTTOT (NRN) = ZERZER
            TOTERR (NRN) = ZERZER
*  |  |  +-------------------------------------------------------------*
*  |  |  |
            DO J = K0, KI1
               IF ( GDSTOR (J) .LT. AZRZRZ ) GDSTOR (J) = ZERZER
               TOTTOT (NRN) = TOTTOT (NRN) + GDSTOR (J)
               TOTERR (NRN) = TOTERR (NRN) + GBSTOR (J)
            END DO
*  |  |  |
*  |  |  +-------------------------------------------------------------*
            TOTERR (NRN) = SQRT (TOTERR (NRN))
            IF ( TOTTOT (NRN) .GT. ZERZER ) TOTERR (NRN) = TOTERR (NRN)
     &                                                   / TOTTOT (NRN)
 9001    CONTINUE
*  |  |
*  |  +----------------------------------------------------------------*
*  |  +----------------------------------------------------------------*
*  |  |  Compute errors:
         DO J = KOFFS + 1, KOFFS + KDIM
            GBSTOR (J) = SQRT (GBSTOR(J))
            IF ( GDSTOR (J) .GT. ZERZER )
     &         GBSTOR (J) = GBSTOR (J) / GDSTOR (J)
         END DO
*  |  |
*  |  +----------------------------------------------------------------*
         KLAST = KOFFS
*  |  +----------------------------------------------------------------*
*  |  |  Now write down:
         DO 7850 JJ=1,IX
            NRN   = JX(JJ)
            K0    = KLAST + 1
            K1    = IMRHGH (NRN) * IZRHGH (NRN) + K0 - 1
            KLAST = K1
**sr commented
c           WRITE (3,*)
c           WRITE (3,*) ' Beam intensity  :',BEAINT,' pr/s'
c           WRITE (3,*)
c           WRITE (3,*) ' Irradiation time:',SNGL(TIRRDT(NIRRDT)),' s'
c           WRITE (3,*)
c           WRITE (3,*) ' Decay time      :',SNGL(TIRDCY(ITD)),   ' s'
c           WRITE (3,*)
c           WRITE (3,*) ' Detector n: ', NRN,'  ',TIURSN (NRN)
c           WRITE (3,*) '    (Region',NRURSN(NRN),' Volume:  ',
c    &           SNGL(VURSNC(NRN)),' cmc,'
c           WRITE (3,*) '     distr. type  :',ITURSN(NRN),'   ,'
c           WRITE (3,*) '     Z_max:',IZRHGH(NRN),', N-Z_max:',
c    &           IMRHGH(NRN)+NMZMNM,', N-Z_min:',NMZMNM+1,
c    &           ')'
c           WRITE (3,*)
c           WRITE (3,*)'    Tot. response (Bq/cmc)',
c    &           SNGL(TOTTOT(NRN)),
c    &           '  +/-',SNGL(TOTERR(NRN))*100.0,' %'
c           WRITE (3,*)'    ( --> Bq              ',
c    &           SNGL(TOTTOT(NRN)*VURSNC(NRN)),
c    &           '  +/-',SNGL(TOTERR(NRN))*100.0,' % )'
c           WRITE (3,*)
c           WRITE (3,*)
c    &           ' **** Residual nuclei distribution  ****',
c    &           ' ****           (Bq/cmc)            ****'
c           WRITE (3,*)
c           WRITE (3,*)
**
            THRESH = TOTTOT(NRN) * 0.001
*************************
**sr commented
c           WRITE (66,*)
c           WRITE (66,*)
c           WRITE (66,*) ' Beam intensity  :',BEAINT,' pr/s'
c           WRITE (66,*)
c           WRITE (66,*) ' Irradiation time:',SNGL(TIRRDT(NIRRDT)),' s'
c           WRITE (66,*)
c           WRITE (66,*) ' Decay time      :',SNGL(TIRDCY(ITD)),   ' s'
c           WRITE (66,*)
c           WRITE (66,*) ' Detector n: ', NRN,'  ',TIURSN (NRN)
c           WRITE (66,*) '    (Region',NRURSN(NRN),' Volume:  ',
c    &           SNGL(VURSNC(NRN)),' cmc,'
c           WRITE (66,*) '     distr. type  :',ITURSN(NRN),'   ,'
c           WRITE (66,*) '     Z_max:',IZRHGH(NRN),', N-Z_max:',
c    &           IMRHGH(NRN)+NMZMNM,', N-Z_min:',NMZMNM+1,
c    &           ')'
c           WRITE (66,*)
c           WRITE (66,*)'    Tot. response (Bq/cmc)',
c    &           SNGL(TOTTOT(NRN)),
c    &           '  +/-',SNGL(TOTERR(NRN))*100.0,' %'
c           WRITE (66,*)'    ( --> Bq              ',
c    &           SNGL(TOTTOT(NRN)*VURSNC(NRN)),
c    &           '  +/-',SNGL(TOTERR(NRN))*100.0,' % )'
c           WRITE (66,*)
c           WRITE (66,*)
c    &           ' **** Residual nuclei distribution  ****',
c    &           ' ****           (Bq/cmc)            ****'
c           WRITE (66,*)
c           WRITE (66,9340)
c           WRITE (66,*)
**
*************
*************
**sr commented
c           WRITE (67,*) NRN, SNGL(TIRRDT(NIRRDT)), SNGL(TIRDCY(ITD))
c           WRITE (67,*) SNGL(TOTTOT(NRN)), SNGL(TOTERR(NRN))
**
*************
            DO IA = 1, 260
               DO IZ = 1, 100
                  RESNUC (IZ,IA) = ZERZER
                  RESERR (IZ,IA) = ZERZER
                  IABEG  (IZ) = 500
                  IAEND  (IZ) = 0
               END DO
            END DO
            ITHR  = 0
            DO IM = 1, IMRHGH (NRN)
               DO IZ  = 1, IZRHGH (NRN)
                  IN  = IM + IZ + NMZMNM
                  IA  = IN + IZ
                  IF ( IA .GE. 1 ) THEN
                     IZN = K0 - 1 + ( IM - 1 ) * IZRHGH (NRN) + IZ
                     RESNUC (IZ,IA) = GDSTOR (IZN)
                     IF ( ABS (GDSTOR (IZN)) .GT. ZERZER ) THEN
                        IABEG (IZ) = MIN ( IABEG (IZ), IA )
                        IAEND (IZ) = MAX ( IAEND (IZ), IA )
                        RESERR (IZ,IA) = GBSTOR (IZN) * 100.D+00
                        IF ( RESNUC (IZ,IA). GT. THRESH ) THEN
                           ITHR = ITHR + 1
                           IATHR (ITHR) = IA
                           IZTHR (ITHR) = IZ
                           IMTHR (ITHR) = 0
                           ACTHR (ITHR) = RESNUC (IZ,IA)
                           ERTHR (ITHR) = RESERR (IZ,IA)
                        END IF
                     END IF
                  END IF
               END DO
            END DO
*
            DO JTHR = 1, ITHR
               IAIZTH (JTHR) = IATHR (JTHR) * 1000 + IZTHR (JTHR) * 10
     &              + IMTHR (JTHR)
               ICORTH (JTHR) = JTHR
            END DO
            CALL IORDIN ( IAIZTH, ICORTH, ITHR )
            DO KTHR = 1, ITHR
               JTHR = ICORTH (KTHR)
               IF ( ICHSYZ ( IZTHR (JTHR), CHSYM ) .NE. IZTHR (JTHR) )
     &              STOP 'ICHSYZ'
               ISOMER = IMTHR (JTHR)
               T12CUR = TIM1O2 ( IATHR (JTHR), IZTHR (JTHR), KKA, KKZ,
     &              T12DDA, BRDDAU, ISOMER )
**sr commented
c              WRITE (66,9350) IATHR (JTHR), CHSYM, IZTHR (JTHR),
c    &              ACTHR (JTHR), ERTHR (JTHR), T12CUR
**
            END DO
            ITHR   = 0
*     
            IASTRT = IMRHGH (NRN) + 2 * IZRHGH (NRN) + NMZMNM
            IZ0    = IZRHGH (NRN) + 1
 8220       CONTINUE
            IZ1 = IZ0 - 1
            IZ0 = MAX ( 1, IZ1 - 10 )
*        IAS1= IASTLN (MIN(MAX(IZ1,1),100))
*        IAS0= IASTLN (MIN(MAX(IZ0,1),100))
*        IA1 = MIN ( IASTRT, IAS1 + MAX ( MIN (IZ1/5,10),  5 ) )
*        IA0 = MAX ( IZ0   , IAS0 - MAX ( MIN (IZ1/3,20), 10 ) )
*        IA1 = MIN ( IA1, IMRHGH (NRN) + 2 * IZ1 + NMZMNM )
*        IA0 = MAX ( IA0, 1 + 2 * IZ0 + NMZMNM )
*        IA0 = MAX ( IA0, 2 )
            IA1 = 0
            IA0 = 500
            DO IZ = IZ0, IZ1
               IA1 = MAX ( IA1, IAEND (IZ) )
               IA0 = MIN ( IA0, IABEG (IZ) )
            END DO
**sr commented
c           WRITE (3,9240) ( IZ, IZ = IZ0, IZ1 )
**
            DO 8230 IA = IA1, IA0, - 1
               LCHECK = .FALSE.
               DO IZ = IZ0, IZ1
                  LCHECK = LCHECK .OR. RESNUC (IZ,IA) .GE. AZRZRZ
               END DO
               IF ( LCHECK ) THEN
**sr commented
c                 WRITE (3,9250) IA, ( RESNUC (IZ,IA), IZ = IZ0, IZ1 )
c                 WRITE (3,9260) 
c    &                        ( CHERR, RESERR (IZ,IA), IZ = IZ0, IZ1 )
**
**sr
               IF (TIRDCY(ITD).GT.1.0E-38) THEN
                  if (irrad.le.0) stop '  IRRAD .LE. 0 !'
                  DO 200 IZ = IZ0, IZ1
                     ADDNUC(NRN,IZ,IA) = ADDNUC(NRN,IZ,IA)+
     &                                   WEIGHT(ITD-1)*RESNUC(IZ,IA)
                     ADDERR(NRN,IZ,IA) = ADDERR(NRN,IZ,IA)+
     &                RESERR(IZ,IA)/100.0E0*WEIGHT(ITD-1)*RESNUC(IZ,IA)
  200             CONTINUE
               ENDIF
**
               END IF
 8230       CONTINUE
            IF ( IZ0 .GT. 1 ) GO TO 8220
**sr
C        CALL ANALYS(RESNUC,RESERR,VURSNC(NRN),NRN,1)
**
**sr commented
c           WRITE (3,*)
***********
c           WRITE (66,*)
**
            ITHR = 0
***********
            IF ( LISMRS ) THEN
               KI0   = KLAST + 1
               KI1   = KLAST + NISMRS (NRN)
               KLAST = KI1
               IF ( NISMRS (NRN) .GT. 0 ) THEN
**sr commented
c                 WRITE (3,*) ' **** Isomers (Bq/cmc)     ****'
c                 WRITE (3,*) '     A        Z       mth'
***********
c                 WRITE(66,*) ' **** Isomers (Bq/cmc)     ****'
c                 WRITE(66,*)
c                 WRITE(66,9360)
c                 WRITE(66,*)
**
***********
               END IF
               DO KS = 1, NISMRS (NRN)
                  IZN = KI0 + KS - 1
                  IF ( GDSTOR(IZN) .GT. AZRZRZ ) THEN
                     CALL ISMRCH ( IA, IZ, IS, KS, ADUMMY, BDUMMY, 
     &                             JDUMMY, IDUMMY )
**sr commented
c                    WRITE (3,9270) IA, IZ, IS, SNGL(GDSTOR(IZN)),
c    &                              SNGL(GBSTOR(IZN)*100.D+00)
**
***********
**sr
                     IF (TIRDCY(ITD).GT.1.0E-38) THEN
                        ADDNUC(NRN,IZ,IA+260) = ADDNUC(NRN,IZ,IA+260)
     &                                      +WEIGHT(ITD-1)*GDSTOR(IZN)
                        ADDERR(NRN,IZ,IA+260) = ADDERR(NRN,IZ,IA+260)+
     &                     GBSTOR(IZN)*WEIGHT(ITD-1)*GDSTOR(IZN)
                     ENDIF
**
                     IF ( GDSTOR(IZN) .GT. THRESH ) THEN
                        ITHR = ITHR + 1
                        IZTHR (ITHR) = IZ
                        IATHR (ITHR) = IA
                        IMTHR (ITHR) = IS
                        ACTHR (ITHR) = GDSTOR (IZN)
                        ERTHR (ITHR) = GBSTOR (IZN) * 100.D+00
                     END IF
***********
                  END IF
               END DO
*
               DO JTHR = 1, ITHR
                  IAIZTH (JTHR) = IATHR (JTHR) * 10000
     &                 + IZTHR (JTHR) * 10 + IMTHR (JTHR)
                  ICORTH (JTHR) = JTHR
               END DO
               CALL IORDIN ( IAIZTH, ICORTH, ITHR )
               DO KTHR = 1, ITHR
                  JTHR = ICORTH (KTHR)
                  IF ( ICHSYZ ( IZTHR (JTHR), CHSYM ) 
     &                 .NE. IZTHR (JTHR) )
     &                 STOP 'ICHSYZ'
                  ISOMER = IMTHR (JTHR)
                  T12CUR = TIM1O2 ( IATHR (JTHR), IZTHR (JTHR), KKA,
     &                 KKZ, T12DDA, BRDDAU, ISOMER )
**sr commented
c                 WRITE (66,9370) IATHR (JTHR), CHSYM, IZTHR (JTHR),
c    &                            IMTHR (JTHR), ACTHR (JTHR),
c    &                            ERTHR (JTHR), T12CUR
**
               END DO
               ITHR = 0
*     
**sr commented
c              WRITE (3,*)
***********
c              WRITE (66,*)
**
***********
            END IF
*  |  Write on binary file the results:
**sr commented
c           WRITE (1)    NRN, TIURSN(NRN), ITURSN(NRN), NRURSN(NRN),
c    &              SNGL(VURSNC (NRN)), IMRHGH(NRN), IZRHGH(NRN), NMZMNM
c           IF ( LEVOLN ) WRITE (1) SNGL (TIRDCY(ITD))
c           WRITE (1) (SNGL(GDSTOR(J)), J = K0, K1)
c           IF ( LISMRS ) THEN
c              WRITE (1) 'ISOMERS:  ', NISMRS (NRN)
c              IF ( KI1 .GE. KI0 ) THEN
c                 WRITE (1) (SNGL(GDSTOR(J)), J = KI0, KI1)
c              ELSE
c                 WRITE (1)
c              END IF
c           END IF
**
 7850    CONTINUE
*  |  |
*  |  +----------------------------------------------------------------*
         KLAST = KOFFS
*  Now write the statistical data, first a record to flag them
**sr commented
c        WRITE (1) 'STATISTICS',IX
**
*  |  +----------------------------------------------------------------*
*  |  |
         DO 7860 JJ=1,IX
            NRN   = JX(JJ)
            K0    = KLAST + 1
            K1    = IMRHGH (NRN) * IZRHGH (NRN) + K0 - 1
            KLAST = K1
**sr commented
c           WRITE (1) SNGL (TOTTOT(NRN)), SNGL (TOTERR(NRN))
**
*     |  These data are zeroes for an evolution file:
            KK0 = KAURSN (NRN) + 1
            KK1 = KAURSN (NRN) + IMRHGH (NRN) + 2 * IZRHGH (NRN) +NMZMNM
**sr commented
c           WRITE (1) ( SNGL (GDSTOR(J)), J = KK0, KK1 )
c           WRITE (1) ( SNGL (GBSTOR(J)), J = KK0, KK1 )
**
*     |  These data are zeroes for an evolution file:
            KK0 = KZURSN (NRN) + 1
            KK1 = KZURSN (NRN) + IZRHGH (NRN)
**sr commented
c           WRITE (1) ( SNGL (GDSTOR(J)), J = KK0, KK1 )
c           WRITE (1) ( SNGL (GBSTOR(J)), J = KK0, KK1 )
**
*     |  Errors on the A/Z matrix:
**sr commented
c           WRITE (1) ( SNGL (GBSTOR(J)), J = K0 , K1  )
**
            IF ( LISMRS ) THEN
               KI0   = KLAST + 1
               KI1   = KLAST + NISMRS (NRN)
               KLAST = KI1
               IF ( KI1 .GE. KI0 ) THEN
**sr commented
c                 WRITE (1) ( SNGL (GBSTOR(J)), J = KI0, KI1 )
**
               ELSE
**sr commented
c                 WRITE (1)
**
               END IF
            END IF
 7860    CONTINUE
*  |  |
*  |  +----------------------------------------------------------------*
 9000 CONTINUE
*  |
*  +-------------------------------------------------------------------*
**sr
      IF (LCYCLE) THEN
         DO 204 JJ=1,IX
            NRN   = JX(JJ)
**sr modified from stand-alone postprocessor
C           WRITE(50,3000) NRN,TIURSN(NRN),NRURSN(NRN),VURSNC(NRN),
C    &                     ITURSN(NRN)
**
 3000       FORMAT(/,1X,'-------------------------------------------',
     &             '-----------------------------------',/,/,1X,
     &             'Isotope yield for detector no.',I3,' (',A,') :',/,
     &             1X,'(Region ',I3,'  Volume: ',E10.4,' cmc,',
     &             ' distribution type: ',I1,')',/)
            IASTRT = IMRHGH (NRN) + 2 * IZRHGH (NRN) + NMZMIN
            IZ0    = IZRHGH (NRN) + 1
  203       CONTINUE
            IZ1 = IZ0 - 1
            IZ0 = MAX ( 1, IZ1 - 10 )
            IA1 = 0
            IA0 = 500
**sr 3.2.00 bug-fix
C           DO IZ = IZ0, IZ1
C              IA1 = MAX ( IA1, IAEND (IZ) )
C              IA0 = MIN ( IA0, IABEG (IZ) )
C           END DO
            IA0 = 1
            IA1 = 260
**
C           WRITE (50,9240) ( IZ, IZ = IZ0, IZ1 )
**sr modified from stand-alone postprocessor
C           DO 201 IA = IA1, IA0, - 1
C              LCHECK = .FALSE.
C              DO IZ = IZ0, IZ1
C                 LCHECK = LCHECK .OR. ADDNUC (NRN,IZ,IA) .GE. 1.E-38
C              END DO
C              IF ( LCHECK ) THEN
C              DO 202 IZ = IZ0, IZ1
C                 IF (ADDNUC(NRN,IZ,IA).GE.1.E-38) THEN
C                    ADDERR(NRN,IZ,IA) = 100.0E0*ADDERR(NRN,IZ,IA)/
C    &                                   ADDNUC(NRN,IZ,IA)
C                 ELSE
C                    ADDERR(NRN,IZ,IA) = 0.0E0
C                 ENDIF
C 202          CONTINUE
C              END IF
C              LCHECK = .FALSE.
C              DO IZ = IZ0, IZ1
C                 LCHECK = LCHECK.OR.ADDNUC(NRN,IZ,IA+260).GE.1.E-38
C              END DO
C              IF ( LCHECK ) THEN
C              DO 205 IZ = IZ0, IZ1
C                 IF (ADDNUC(NRN,IZ,IA+260).GE.1.E-38) THEN
C                    ADDERR(NRN,IZ,IA+260) =
C    &               100.0E0*ADDERR(NRN,IZ,IA+260)/ADDNUC(NRN,IZ,IA+260)
C                 ELSE
C                    ADDERR(NRN,IZ,IA+260) = 0.0E0
C                 ENDIF
C 205          CONTINUE
CC             WRITE (50,9250) IA, ( ADDNUC (NRN,IZ,IA), IZ = IZ0, IZ1 )
CC             WRITE (50,9260) (CHERR,ADDERR(NRN,IZ,IA),IZ=IZ0,IZ1)
C              END IF
C 201       CONTINUE
**
            IF ( IZ0 .GT. 1 ) GO TO 203
            IF (ICOOL.LT.NCOOL) GOTO 210
**sr modified from stand-alone postprocessor
C           CALL PRTISO(ADDNUC,ADDERR,NRN,50)
            CALL PRTISO_NEW(IZEVO,IAEVO,ADDNUC,
     &                      NISO,IZDEC,IADEC,ISODEC,YLDDEC)
**
  204    CONTINUE
      ENDIF
**
      CLOSE (UNIT= 1)
      CLOSE (UNIT= 3)
      CLOSE (UNIT=66)
      CLOSE (UNIT=67)
*
 9240 FORMAT ( /, 1X,' A \\ Z',2X,I3,3X,10(5X,I3,3X) )
 9250 FORMAT ( 1X,I3,11(3X,1P,E8.2) )
 9260 FORMAT ( 1X,3X,11(2X,A3,F4.1,' %') )
 9270 FORMAT(1X,I8,I8,I8,4X,1P,E8.2,' +/- ',0P,F4.1,' %')
 9340 FORMAT(7X,'A ',2X,'Sym.',2X,6X,'Z ',4X,' Bq/cmc  ',5X,'Err.',2X,
     &       4X,'  T 1/2 (s)')
 9350 FORMAT(1X,I8,3X,A2,3X,I8,4X,1P,E9.3,' +/- ',0P,F4.1,' %',
     &       4X,1P,E9.3)
 9360 FORMAT(7X,'A ',2X,'Sym.',2X,6X,'Z ',' mth',4X,' Bq/cmc  ',
     &       5X,'Err.',2X,4X,'  T 1/2 (s)')
 9370 FORMAT(1X,I8,3X,A2,3X,I8,I4,4X,1P,E9.3,' +/- ',0P,F4.1,' %',
     &       4X,1P,E9.3)
**sr modified from stand-alone postprocessor
C     STOP
**
      END
*
*===irrcyc=============================================================*
*
      SUBROUTINE IRRCYC(TDCY)

************************************************************************
* Preprocessor for USRSUW to generate irradiation cycles.              *
* Last change 22.4.99 by S. Roesler                                    *
************************************************************************

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE

      PARAMETER (LUNINP = 10)
      PARAMETER (ZERO = 0.0D0)

      CHARACTER CLINE*78,CTIME*1,CWHAT*77,CTUNIT*1

      PARAMETER (MXIRRA = 5000, MXCOOL = 1, MXFILE = 170)
      LOGICAL LCYCLE
      CHARACTER CINFLE*80,CRUNTT*80,CRUNTM*32
      COMMON /IRRPRP/ PRTSEC,TIRRAD,TDECY(MXCOOL,MXIRRA),TOTWGT,
     &                WEIGHT(MXIRRA),NIRRAD,NCOOL,LCYCLE,
     &                CINFLE(MXFILE),NFILE,CRUNTT,CRUNTM,NPRIM

      DIMENSION CTUNIT(5),NSECON(5),TCOOL(MXCOOL)
      DATA NSECON /  1, 60,3600,86400,31536000/
      DATA CTUNIT /'s','m', 'h',  'd',     'y'/

      CALL OAUXFI('irrcyc.inp',LUNINP,'UNKNOWN',IERR)

      PRTSEC = ZERO
      TIRRAD = ZERO
      TGAP   = ZERO
      DO 1 I=1,MXCOOL
         TCOOL(I) = ZERO
    1 CONTINUE
      DO 12 I=1,MXIRRA
         WEIGHT(I) = 1.0D0
   12 CONTINUE
      NIRRAD = 0
      NCOOL  = 0

* number of beam particles per second
      READ(LUNINP,*) PRTSEC
      IF (PRTSEC.LE.ZERO) THEN
         WRITE(*,*) 
     &     ' Warning! Number of beam particles per second < 0: ',PRTSEC
         STOP
      ENDIF

* irradiation time
      READ(LUNINP,'(A78)') CLINE
      READ(CLINE,'(A1,A77)') CTIME,CWHAT
      READ(CWHAT,*) TIRRAD
      DO 2 ISECON=1,5
         IF (CTIME.EQ.CTUNIT(ISECON)) GOTO 3
    2 CONTINUE
      WRITE(*,*) ' Warning! Time unit not recognized: ',CTIME
      STOP
    3 CONTINUE
      TIRRAD = TIRRAD*NSECON(ISECON)
      IF (TIRRAD.LE.ZERO) THEN
         WRITE(*,*) 
     &     ' Warning! Inconsistent irradiation time: ',TIRRAD
         STOP
      ENDIF

* cooling time in between irradiations
      READ(LUNINP,'(A78)') CLINE
      READ(CLINE,'(A1,A77)') CTIME,CWHAT
      READ(CWHAT,*) TGAP
      DO 4 ISECON=1,5
         IF (CTIME.EQ.CTUNIT(ISECON)) GOTO 5
    4 CONTINUE
      WRITE(*,*) ' Warning! Time unit not recognized: ',CTIME
      STOP
    5 CONTINUE
      TGAP = TGAP*NSECON(ISECON)
      IF (TGAP.LT.ZERO) THEN
         WRITE(*,*) 
     &     ' Warning! Inconsistent gap-time: ',TGAP
         STOP
      ENDIF

* number of irradiations
      READ(LUNINP,*) NIRRAD
      IF (NIRRAD.EQ.0) THEN
         WRITE(*,*) ' Warning! Number of irradiations = 0: ',NIRRAD
         STOP
      ELSEIF (NIRRAD.LT.0) THEN
         NIRRAD = -NIRRAD
         DO 13 I=1,NIRRAD
            READ(LUNINP,*) WEIGHT(NIRRAD+1-I)
   13    CONTINUE
      ENDIF
      IF (NIRRAD.GT.MXIRRA) THEN
         WRITE(*,*) ' Warning! Too many irradiations: ',NIRRAD,MXIRRA
         STOP
      ENDIF

* cooling time after the end of the last irradiation
C   6 CONTINUE
C     READ(LUNINP,'(A78)',END=9000) CLINE
C     READ(CLINE,'(A1,A77)') CTIME,CWHAT
C     NCOOL = NCOOL+1
C     IF (NCOOL.GT.MXCOOL) THEN
C        WRITE(*,*) ' Warning! Too many cooling times: ',NCOOL,MXCOOL
C        STOP
C     ENDIF
C     READ(CWHAT,*) TCOOL(NCOOL)
C     DO 7 ISECON=1,5
C        IF (CTIME.EQ.CTUNIT(ISECON)) GOTO 8
C   7 CONTINUE
C     WRITE(*,*) ' Warning! Time unit not recognized: ',CTIME
C     STOP
C   8 CONTINUE
C     TCOOL(NCOOL) = TCOOL(NCOOL)*NSECON(ISECON)
C     IF (TCOOL(NCOOL).LE.ZERO) THEN
C        WRITE(*,*) 
C    &     ' Warning! Inconsistent cooling time: ',TCOOL(NCOOL)
C        STOP
C     ENDIF
C     GOTO 6
      NCOOL = NCOOL+1
      TCOOL(NCOOL) = TDCY
      
 9000 CONTINUE
      CLOSE(LUNINP)

C     WRITE(50,1005)
 1005 FORMAT(/,1X,'Irradiation properties :',/)
C     WRITE(50,1000) PRTSEC
C     WRITE(51,'(1X,E10.4)') PRTSEC
 1000 FORMAT(9X,'Number of beam particles per second:',2X,1P,E10.4)
C     WRITE(50,1001) TIRRAD
 1001 FORMAT(24X,'Irradiation time (s):',2X,1P,E10.4)
C     WRITE(50,1002) TGAP
 1002 FORMAT(32X,'Gap time (s):',2X,1P,E10.4)
C     WRITE(50,1003) NIRRAD
 1003 FORMAT(22X,'Number of irradiations:',2X,I10)
C     DO 9 I=1,NCOOL
C        WRITE(50,1004) TCOOL(I)
C        WRITE(51,'(1X,1P,E10.4)') TCOOL(I)
C   9 CONTINUE
 1004 FORMAT(28X,'Cooling time (s):',2X,1P,E10.4)

      DO 10 I=1,NCOOL
         DO 11 J=1,NIRRAD
            TDECY(I,J) = DBLE(J-1)*(TIRRAD+TGAP)+TCOOL(I)
   11    CONTINUE
   10 CONTINUE
      
      RETURN
      END
*
*===prtiso_new=========================================================*
*
      SUBROUTINE PRTISO_NEW(IZEVO,IAEVO,YIELD,
     &                      NISO,IZDEC,IADEC,ISODEC,YLDDEC)

************************************************************************
* Yield of daughter isotopes which are gamma emitters.                 *
*                                                                      *
* Input:    IZEVO,IAEVO charge/mass number of mother isotope           *
*           YIELD       array of yield of daughter isotopes            *
*                                                                      *
* Output:   NISO        number of daughter isotopes                    *
*           IZDEC,IADEC charge/mass number of daughter isotopes        *
*                       (isomer if Iadec>260)                          *
*           ISODEC      isomer flag of daughter isotopes               *
*           YLDDEC      yield of daughter isotopes                     *
*                                                                      *
* This version dated 10.1.03   is written by S. Roesler.               *
************************************************************************

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'

      PARAMETER (MXDET = 1)
      DIMENSION YIELD(MXDET,100,520)

      PARAMETER (MAXISO=20)
      DIMENSION IZDEC(MAXISO),IADEC(MAXISO),ISODEC(MAXISO),
     &          YLDDEC(MAXISO)

      LOGICAL LBIAS,LBINRY,LREGDR,LNA24
      COMMON /DOSPAR/ FBIAS,WMINBI,FNA24,KBIAS,KDUMP,
     &                LREGDR(MXXRGN),LBIAS(MXXRGN),LNA24(MXXRGN),
     &                LBINRY

      NISO = 0
      DO 1 IZ=1,IZEVO+1
         DO 2 IA=IZ,IAEVO

            IF (YIELD(1,IZ,IA).GE.1.0D-38) THEN
               FLAG   = -1.0D0
               WEIGT1 = 0.0D0
               WEIGT2 = 0.0D0
               WEIGT3 = 0.0D0
               ISO    = 1
               IF (IA.GT.260) ISO = 2
               IF ((KDUMP.EQ.0).OR.(KDUMP.EQ.1).OR.(KDUMP.EQ.4))
     &            CALL READGL(IZ,IA,ISO,FLAG,WEIGT1)
               IF ((KDUMP.EQ.0).OR.(KDUMP.EQ.2).OR.(KDUMP.EQ.4))
     &            CALL READPL(IZ,IA,ISO,FLAG,WEIGT2)
               IF ((KDUMP.EQ.0).OR.(KDUMP.EQ.3))
     &            CALL READEL(IZ,IA,ISO,FLAG,WEIGT3)
*
               IF ((WEIGT1.GT.0.0D0).OR.(WEIGT2.GT.0.0D0)
     &                              .OR.(WEIGT3.GT.0.0D0)) THEN
                  NISO = NISO+1
                  IF (NISO.GT.MAXISO) STOP ' PRTISO: NISO > MAXISO !'
                  IZDEC(NISO)  = IZ
                  IADEC(NISO)  = IA
                  ISODEC(NISO) = ISO
                  YLDDEC(NISO) = YIELD(1,IZ,IA)
C                 WRITE(50,*) 
C    &               ' isotope: ',IZ,IA,ISO,YIELD(1,IZ,IA)
               ENDIF
            ENDIF

            IF (YIELD(1,IZ,IA+260).GE.1.0D-38) THEN
* check for gamma ray emitter
               FLAG   = -1.0D0
               WEIGT1 = 0.0D0
               WEIGT2 = 0.0D0
               WEIGT3 = 0.0D0
               ISO    = 2
               IF ((KDUMP.EQ.0).OR.(KDUMP.EQ.1).OR.(KDUMP.EQ.4))
     &            CALL READGL(IZ,IA,ISO,FLAG,WEIGT1)
               IF ((KDUMP.EQ.0).OR.(KDUMP.EQ.2).OR.(KDUMP.EQ.4))
     &            CALL READPL(IZ,IA,ISO,FLAG,WEIGT2)
               IF ((KDUMP.EQ.0).OR.(KDUMP.EQ.3))
     &            CALL READEL(IZ,IA,ISO,FLAG,WEIGT3)
*
               IF ((WEIGT1.GT.0.0D0).OR.(WEIGT2.GT.0.0D0)
     &                              .OR.(WEIGT3.GT.0.0D0)) THEN
                  NISO = NISO+1
                  IF (NISO.GT.MAXISO) STOP ' PRTISO: NISO > MAXISO !'
                  IZDEC(NISO)  = IZ
                  IADEC(NISO)  = IA
                  ISODEC(NISO) = ISO
                  YLDDEC(NISO) = YIELD(1,IZ,IA+260)
C                 WRITE(50,*) 
C    &               ' isotope: ',IZ,IA,ISO,YIELD(1,IZ,IA+260)
               ENDIF
            ENDIF

    2    CONTINUE
    1 CONTINUE

      RETURN
      END
*
*=== bspect ===========================================================*
*
      SUBROUTINE BSPECT(KZ,KA,KISO,EP,WEIGHT)
*
*----------------------------------------------------------------------*
*                                                                      *
*  This routine is part of the package of user files for the           *
*  calculation of dose rates from induced activity.                    *
*                                                                      *
*  Sampling from energy spectra for beta radiation. The spectra are    *
*  read from ICRP38.BET of the Nucdecay-program.                       *
*                                                                      *
*  This version dated 23.01.03 by Stefan Roesler.                      *
*                                                                      *
*----------------------------------------------------------------------*
*
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE

      LOGICAL LFIRST
      CHARACTER*78 CLINE

      PARAMETER (LUNBET = 32)

      DIMENSION EBIN(150)
      DATA EBIN /
     & 0.0E+00,1.0E-04,1.1E-04,1.2E-04,1.3E-04,1.4E-04,1.5E-04,1.6E-04,
     & 1.8E-04,2.0E-04,2.2E-04,2.4E-04,2.6E-04,2.8E-04,3.0E-04,3.2E-04,
     & 3.6E-04,4.0E-04,4.5E-04,5.0E-04,5.5E-04,6.0E-04,6.5E-04,7.0E-04,
     & 7.5E-04,8.0E-04,8.5E-04,9.0E-04,1.0E-03,1.1E-03,1.2E-03,1.3E-03,
     & 1.4E-03,1.5E-03,1.6E-03,1.8E-03,2.0E-03,2.2E-03,2.4E-03,2.6E-03,
     & 2.8E-03,3.0E-03,3.2E-03,3.6E-03,4.0E-03,4.5E-03,5.0E-03,5.5E-03,
     & 6.0E-03,6.5E-03,7.0E-03,7.5E-03,8.0E-03,8.5E-03,9.0E-03,1.0E-02,
     & 1.1E-02,1.2E-02,1.3E-02,1.4E-02,1.5E-02,1.6E-02,1.8E-02,2.0E-02,
     & 2.2E-02,2.4E-02,2.6E-02,2.8E-02,3.0E-02,3.2E-02,3.6E-02,4.0E-02,
     & 4.5E-02,5.0E-02,5.5E-02,6.0E-02,6.5E-02,7.0E-02,7.5E-02,8.0E-02,
     & 8.5E-02,9.0E-02,1.0E-01,1.1E-01,1.2E-01,1.3E-01,1.4E-01,1.5E-01,
     & 1.6E-01,1.8E-01,2.0E-01,2.2E-01,2.4E-01,2.6E-01,2.8E-01,3.0E-01,
     & 3.2E-01,3.6E-01,4.0E-01,4.5E-01,5.0E-01,5.5E-01,6.0E-01,6.5E-01,
     & 7.0E-01,7.5E-01,8.0E-01,8.5E-01,9.0E-01,1.0E+00,1.1E+00,1.2E+00,
     & 1.3E+00,1.4E+00,1.5E+00,1.6E+00,1.8E+00,2.0E+00,2.2E+00,2.4E+00,
     & 2.6E+00,2.8E+00,3.0E+00,3.2E+00,3.6E+00,4.0E+00,4.5E+00,5.0E+00,
     & 5.5E+00,6.0E+00,6.5E+00,7.0E+00,7.5E+00,8.0E+00,8.5E+00,9.0E+00,
     & 1.0E+01,1.1E+01,1.2E+01,1.3E+01,1.4E+01,1.5E+01,1.6E+01,1.8E+01,
     & 2.0E+01,2.2E+01,2.4E+01,2.6E+01,2.8E+01,3.0E+01/
*
      CHARACTER CISO*2,CISOIN*2
      DIMENSION CISO(109)
      DATA CISO /
     &    'H ','He','Li','Be','B ','C ','N ','O ','F ','Ne','Na',
     &    'Mg','Al','Si','P ','S ','Cl','Ar','K ','Ca','Sc','Ti',
     &    'V ','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As',
     &    'Se','Br','Kr','Rb','Sr','Y ','Zr','Nb','Mo','Tc','Ru',
     &    'Rh','Pd','Ag','Cd','In','Sn','Sb','Te','I ','Xe','Cs',
     &    'Ba','La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy',
     &    'Ho','Er','Tm','Yb','Lu','Hf','Ta','W ','Re','Os','Ir',
     &    'Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra',
     &    'Ac','Th','Pa','U ','Np','Pu','Am','Cm','Bk','Cf','Es',
     &    'Fm','Md','No','Lr','Rf','Ha','Sg','Ns','Hs','Mt'/

      PARAMETER (IZMAX  = 100,
     &           IAMAX  = 257,
     &           NEBIN  = 100,
     &           IEMAX  = 135,
     &           ISOMAX = 546)
      DIMENSION IPOINT(2,IZMAX,IAMAX),
     &          EBET(ISOMAX,NEBIN),SPCBET(ISOMAX,NEBIN),
     &          SPECTR(IEMAX)

      DATA LFIRST /.TRUE./
*
*----------------------------------------------------------------------*
* Initializations, read beta spectra
*
      IF (LFIRST) THEN
         DO 1 IZ=1,IZMAX
            DO 2 IA=1,IAMAX
               IPOINT(1,IZ,IA) = 0
               IPOINT(2,IZ,IA) = 0
    2       CONTINUE
    1    CONTINUE
         DO 3 IS=1,ISOMAX
            DO 4 IE=1,NEBIN
               EBET(IS,IE)   = 0.0D0
               SPCBET(IS,IE) = 0.0D0
    4       CONTINUE
    3    CONTINUE
         DO 5 IE=1,IEMAX
            SPECTR(IE) = 0.0D0
    5    CONTINUE
         ISO   = 0
         CALL OAUXFI('ICRP38.BET',LUNBET,'OLD',IERR)
*
   10    CONTINUE
         READ(LUNBET,'(A78)',END=21) CLINE
*
*   determine charge and mass
*
         IF (CLINE(2:2).EQ.'-') THEN
            CISOIN(1:2) = CLINE(1:1)//' '
            IF ((CLINE(5:5).EQ.'m').OR.(CLINE(5:5).EQ.'a')
     &                             .OR.(CLINE(5:5).EQ.'b')) THEN
               IS = 1
               IF (CLINE(5:5).EQ.'m') IS = 2
               READ(CLINE,'(2X,I2)') IA
            ELSEIF ((CLINE(6:6).EQ.'m').OR.(CLINE(6:6).EQ.'a')
     &                                 .OR.(CLINE(6:6).EQ.'b')) THEN
               IS = 1
               IF (CLINE(6:6).EQ.'m') IS = 2
               READ(CLINE,'(2X,I3)') IA
            ELSE
               IS = 1
               READ(CLINE,'(2X,I3)') IA
            ENDIF
         ELSEIF (CLINE(3:3).EQ.'-') THEN
            CISOIN = CLINE(1:2)
            IF ((CLINE(6:6).EQ.'m').OR.(CLINE(6:6).EQ.'a')
     &                             .OR.(CLINE(6:6).EQ.'b')) THEN
               IS = 1
               IF (CLINE(6:6).EQ.'m') IS = 2
               READ(CLINE,'(3X,I2)') IA
            ELSEIF ((CLINE(7:7).EQ.'m').OR.(CLINE(7:7).EQ.'a')
     &                                 .OR.(CLINE(7:7).EQ.'b')) THEN
               IS = 1
               IF (CLINE(7:7).EQ.'m') IS = 2
               READ(CLINE,'(3X,I3)') IA
            ELSE
               IS = 1
               READ(CLINE,'(3X,I3)') IA
            ENDIF
         ELSE
            WRITE(6,*) ' Inconsistent record : '
            WRITE(6,*) CLINE
            STOP
         ENDIF
         DO 11 I=1,109
            IF (CISOIN.EQ.CISO(I)) THEN
               IZ = I
               GOTO 12
            ENDIF
   11    CONTINUE
   12    CONTINUE
*
*   read spectra and normalize to unity
*
         READ(LUNBET,*) EMAX
         ISO = ISO+1
         IF (ISO.GT.ISOMAX) STOP ' Iso > Isomax !'
         IPOINT(IS,IZ,IA) = ISO
         IE     = 0
         ISTART = 0
C        write(90,'(A2,3I4)') CISOIN,IZ,IA,IS
C        write(91,'(A2,3I4)') CISOIN,IZ,IA,IS
C        write(92,'(A2,3I4)') CISOIN,IZ,IA,IS
   14    CONTINUE
         IE = IE+1
         IF (IE.GT.IEMAX) STOP ' Ie > Iemax ! '
         READ(LUNBET,'(E9.3)') SPECTR(IE)
         IF (SPECTR(IE).NE.0.0D0) THEN
            ISTART = 1
C           write(90,'(2E15.5)') EBIN(IE),SPECTR(IE)
         ENDIF
         IF ((SPECTR(IE).EQ.0.0D0).AND.(ISTART.EQ.1)) THEN
            PEMI  = 0.0D0
            EAVG  = 0.0D0
            DO 15 I=2,IE
               IF (I.EQ.IE) THEN
                  IF (EMAX.LT.EBIN(I-1)) GOTO 15
                  E = (EBIN(I-1)+EMAX)/2.0D0
               ELSE
                  E = (EBIN(I-1)+EBIN(I))/2.0D0
               ENDIF
               SBIN = (SPECTR(I-1)+SPECTR(I))/2.0D0
               PEMI = PEMI+(EBIN(I)-EBIN(I-1))*SBIN
               EAVG = EAVG+E*(EBIN(I)-EBIN(I-1))*SBIN
   15       CONTINUE
            DO 16 I=1,IE
               SPECTR(I) = SPECTR(I)/PEMI
   16       CONTINUE
*
*   re-bin spectra
*
            DEBIN = EMAX/DBLE(NEBIN)
            PEMI1 = 0.0D0
            DO 17 I=1,NEBIN
               EBET(ISO,I) = DBLE(I)*DEBIN
               E = EBET(ISO,I)-DEBIN/2.0D0
               DO 18 J=1,150
                  IF (EBIN(J).GT.E) THEN
                     I0 = J-1
                     I1 = J
                     GOTO 19
                  ENDIF
   18          CONTINUE
   19          CONTINUE
               FRAC = (E-EBIN(I0))/(EBIN(I1)-EBIN(I0))
               SPCBET(ISO,I) = SPECTR(I0)+FRAC*(SPECTR(I1)-SPECTR(I0))
               PEMI1 = PEMI1+DEBIN*SPCBET(ISO,I)
C              write(91,'(2E15.5)') E,SPCBET(ISO,I)*PEMI
   17       CONTINUE
*
*   normalize to unity
*
            DO 20 I=1,NEBIN
               SPCBET(ISO,I) = SPCBET(ISO,I)/PEMI1
C              write(92,'(3E15.5)') 
C    &            (DBLE(I)-0.5D0)*DEBIN,SPCBET(ISO,I),DEBIN
   20       CONTINUE

            GOTO 10
         ELSE
            GOTO 14
         ENDIF
      
   21    CONTINUE

         LFIRST = .FALSE.
      ENDIF
*
*----------------------------------------------------------------------*
* Sample energy and calculate weight for given isotope
*
      WEIGHT = 0.0D0
      EP     = 0.0D0
      ISO    = IPOINT(KISO,KZ,KA)
      IF (ISO.EQ.0) RETURN

      EP     = EBET(ISO,NEBIN)*FLRNDM(EMAX)
      DEBIN  = EBET(ISO,NEBIN)/DBLE(NEBIN)
      IBIN   = INT(EP/DEBIN)+1
      IF (IBIN.GT.NEBIN) IBIN = NEBIN
      WEIGHT = DBLE(NEBIN)*DEBIN*SPCBET(ISO,IBIN)
      EP     = 1.D-3*EP

      END

*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
*
*===prtiso=============================================================*
*
      SUBROUTINE PRTISO(YIELD,ERROR,IDET,LUNOUT)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE

      INCLUDE '(DIMPAR)'
      INCLUDE '(ISOTOP)'

      PARAMETER (ZERO = 0.0D0, TWOTWO = 2.0D0, FOUFOU = 4.0D0,
     &           ANGLGB = 5.0D-16, AINFNT = 1.0D+30)

      PARAMETER (EFRAC=0.01D0)

      CHARACTER CCISO*1,CISO*2,CERR*3,CTUNIT*1,CDCY*2,CMETA*1

      PARAMETER (MXDET = 1)

      DIMENSION CISO(109),YIELD(MXDET,100,520),ERROR(MXDET,100,520),
     &          TOTYLD(2,MXDET),TOTERR(2,MXDET),IDXIS1(5000),
     &          IDXIS2(5000),IDXIS3(5000),IDXIS4(5000),
     &          T12(5000),IZISO(5000),IAISO(5000),YLD(5000),
     &          CDCY(5000),EDEPOS(5000),NENERG(100,520)
      DIMENSION KA(4),KZ(4),BRDAUG(4),T12DAU(4)
      DIMENSION CTUNIT(5),NSECON(5)

      CHARACTER*4 CDECAY,CTEMP
      DIMENSION GAMENE(5000),GAMPRO(5000),CDECAY(5000),GAMEFF(5000)

      DIMENSION EGAM(50000),ETOGAM(50000),AGAM(50000),
     &          IZGAM(50000),IAGAM(50000),IDXIS5(50000)

      DATA NSECON /  1, 60,3600,86400,31536000/
      DATA CTUNIT /'s','m', 'h',  'd',     'y'/

      DATA CISO /
     &    ' H','He','Li','Be',' B',' C',' N',' O',' F','Ne','Na',
     &    'Mg','Al','Si',' P',' S','Cl','Ar',' K','Ca','Sc','Ti',
     &    ' V','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As',
     &    'Se','Br','Kr','Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru',
     &    'Rh','Pd','Ag','Cd','In','Sn','Sb','Te',' I','Xe','Cs',
     &    'Ba','La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy',
     &    'Ho','Er','Tm','Yb','Lu','Hf','Ta',' W','Re','Os','Ir',
     &    'Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra',
     &    'Ac','Th','Pa',' U','Np','Pu','Am','Cm','Bk','Cf','Es',
     &    'Fm','Md','No','Lr','Rf','Ha','Sg','Ns','Hs','Mt'/
      DATA CERR /'+/-'/

      DO 100 I=1,5000
         IDXIS1(I) = 0
         IDXIS2(I) = 0
         IDXIS3(I) = 0
         IDXIS4(I) = 0
         T12(I) = ZERO
         YLD(I) = ZERO
  100 CONTINUE
      EDEPTO = ZERO
      NISO = 0
      NGAM = 0
      EGAMTO = ZERO

      WRITE(LUNOUT,1000) 
 1000 FORMAT(/,/,6X,'Z   A   Tot. Yield',20X,'T1/2',8X,'N_gamma',
     &       /,15X,'(Bq/cmc)',/)
      TOTYLD(1,IDET) = 0.0D0
      TOTYLD(2,IDET) = 0.0D0
      TOTERR(1,IDET) = 0.0D0
      TOTERR(2,IDET) = 0.0D0
      DO 1 IZ=1,100
         DO 2 IA=1,520
            IIA = IA
            IF (IA.GT.260) IIA = IA-260
            IF (YIELD(IDET,IZ,IA).GE.1.0D-38) THEN
               TOTYLD(1,IDET) = TOTYLD(1,IDET)+YIELD(IDET,IZ,IA)
               TOTERR(1,IDET) = TOTERR(1,IDET)+ERROR(IDET,IZ,IA)
     &                        *YIELD(IDET,IZ,IA)/100.0D0
               NISO = NISO+1
               IDXIS1(NISO) = NISO
               IDXIS2(NISO) = NISO
               IDXIS4(NISO) = NISO
               IZISO(NISO)  = IZ
               IAISO(NISO)  = IA
               YLD(NISO)    = YIELD(IDET,IZ,IA)
* check gamma-transition data
               write(*,*) ' checking g-transition table for Z,A = ',
     &                    IZ,IA
               CALL EGAMMA(IZ,IA,GAMENE,GAMPRO,THLIFE,CDECAY,
     &                                         NENERG(IZ,IA))
               IF (NENERG(IZ,IA).GT.0) THEN
                  DO 20 K=1,5
                     T = THLIFE/DBLE(NSECON(K))
                     KK = K
                     IF (T.LT.100.0D0) GOTO 21
   20             CONTINUE
   21             CONTINUE
                  THLIFE = T
                  ITUNIT = KK
                  IF ((ITUNIT.EQ.5).AND.(THLIFE.LT.1.0D0)) THEN
                     THLIFE = THLIFE*DBLE(NSECON(5))/DBLE(NSECON(4))
                     ITUNIT = 4
                  ENDIF
                  GAMTOT = 0.0D0
                  DO 23 N=1,NENERG(IZ,IA)
                     GAMEFF(N) = GAMENE(N)*ABS(GAMPRO(N))/100.0D0
                     GAMTOT = GAMTOT+GAMEFF(N)
                     IDXIS3(N) = N
   23             CONTINUE
                  CALL SORT(GAMEFF,IDXIS3,NENERG(IZ,IA),1,
     &                                    NENERG(IZ,IA),2)
                  TOTYLD(2,IDET) = TOTYLD(2,IDET)+YIELD(IDET,IZ,IA)
                  TOTERR(2,IDET) = TOTERR(2,IDET)+ERROR(IDET,IZ,IA)
     &                           *YIELD(IDET,IZ,IA)/100.0D0
               ENDIF
               write(*,*) ' ',NENERG(IZ,IA),' energies found'
* half-life
               JSOMER = 0
**sr
C              IF (IA.GT.260) JSOMER = -1
               IF (IA.GT.260) JSOMER = 1
**
               IF ( IZ .LT. IIA ) THEN
                  T12(NISO) = TIM1O2(IIA,IZ,KA,KZ,T12DAU,BRDAUG,JSOMER)
               ELSE
                  T12(NISO) = 0.0D0
               ENDIF
               DO 3 K=1,5
                  THALF = T12(NISO)/DBLE(NSECON(K))
                  KK = K
                  IF (THALF.LT.100.0D0) GOTO 4
    3          CONTINUE
    4          CONTINUE
               IF ((KK.EQ.5).AND.(THALF.LT.1.0D0)) THEN
                  THALF = THALF*DBLE(NSECON(5))/DBLE(NSECON(4))
                  KK = 4
               ENDIF
 
               IACURR = IIA
               IZCURR = IZ
               ACURR  = IIA
               ZCURR  = IZ
               IF ( IZ .LT. IIA ) THEN
                  AMORIG = ENKNOW ( ACURR, ZCURR, IZFLG )
               ELSE
                  write(*,*) ' iz < iia : ',IZ,IIA
                  IZFLG  = 0
                  AMORIG = AINFNT
               ENDIF
               CDCY(NISO) = '  '

* check for alpha emitter
               IZALPH = IZ - 2
               ZCURR  = IZALPH
               IAALPH = IIA - 4
               AALPHA = IAALPH
               IF ( IZALPH .GE. 0 .AND. IZALPH .LT. IAALPH ) THEN
                  AMMALP = ENKNOW ( FOUFOU, TWOTWO, IZFLG )
                  AMALRS = ENKNOW ( AALPHA,  ZCURR, IZFLG )
               ELSE
                  IZFLG  = 0
                  AMMALP = AINFNT
                  AMALRS = AINFNT
               END IF
               QALPHA = AMORIG - AMMALP - AMALRS
               IF ( QALPHA .GE. ANGLGB ) CDCY(NISO) = 'al'
 
* check for beta- emitter
               ACURR  = IIA
               ZCURR  = IZ
               IZBETM = IZ + 1
               ZCURR  = IZBETM
               IF ( IZBETM .GE. 0 .AND. IZBETM .LE. IIA ) THEN
                  AMBETM = ENKNOW ( ACURR, ZCURR, IZFLG )
               ELSE
                  IZFLG  = 0
                  AMBETM = AINFNT
               END IF
               QBETAM = AMORIG - AMBETM
               IF ( QBETAM .GE. ANGLGB ) CDCY(NISO) = 'b-'

* check for beta+ emitter
               ACURR  = IIA
               ZCURR  = IZ
               IZBPEC = IZ - 1
               ZCURR  = IZBPEC
               IF ( IZBPEC .GE. 0 .AND. IZBPEC .LE. IIA ) THEN
                  AMBPEC = ENKNOW ( ACURR, ZCURR, IZFLG )
               ELSE
                  IZFLG  = 0
                  AMBPEC = AINFNT
               END IF
               QBEPEC = AMORIG - AMBPEC
               IF ( QBEPEC .GE. ANGLGB ) CDCY(NISO) = 'b+'

               IF (NENERG(IZ,IA).GT.0) THEN
                  IF (IA.LE.260) THEN
                    IF (THALF.LT.1.E12) THEN
                       WRITE(LUNOUT,1004)
     &                          CISO(IZ),IZ,IIA,YIELD(IDET,IZ,IA),
     &                          CERR,ERROR(IDET,IZ,IA),THALF,CTUNIT(KK),
     &                          CDCY(NISO),NENERG(IZ,IA),THLIFE,
     &                          CTUNIT(ITUNIT)
                       WRITE(51,1204) 
     &                          CISO(IZ),IZ,IIA,YIELD(IDET,IZ,IA),
     &                          ERROR(IDET,IZ,IA),THALF,CTUNIT(KK),
     &                          THLIFE,CTUNIT(ITUNIT)
                    ELSE
                       WRITE(LUNOUT,1104)
     &                          CISO(IZ),IZ,IIA,YIELD(IDET,IZ,IA),
     &                          CERR,ERROR(IDET,IZ,IA),THALF,CTUNIT(KK),
     &                          CDCY(NISO),NENERG(IZ,IA),THLIFE,
     &                          CTUNIT(ITUNIT)
                       WRITE(51,1304) 
     &                          CISO(IZ),IZ,IIA,YIELD(IDET,IZ,IA),
     &                          ERROR(IDET,IZ,IA),THALF,CTUNIT(KK),
     &                          THLIFE,CTUNIT(ITUNIT)
                    ENDIF
                  ELSE
                    IF (THALF.LT.1.E12) THEN
                       WRITE(LUNOUT,1017)
     &                          'm'//CISO(IZ),IZ,IIA,YIELD(IDET,IZ,IA),
     &                          CERR,ERROR(IDET,IZ,IA),THALF,CTUNIT(KK),
     &                          CDCY(NISO),NENERG(IZ,IA),THLIFE,
     &                          CTUNIT(ITUNIT)
                       WRITE(51,1217) 
     &                          'm'//CISO(IZ),IZ,IIA,YIELD(IDET,IZ,IA),
     &                          ERROR(IDET,IZ,IA),THALF,CTUNIT(KK),
     &                          THLIFE,CTUNIT(ITUNIT)
                    ELSE
                       WRITE(LUNOUT,1117)
     &                          'm'//CISO(IZ),IZ,IIA,YIELD(IDET,IZ,IA),
     &                          CERR,ERROR(IDET,IZ,IA),THALF,CTUNIT(KK),
     &                          CDCY(NISO),NENERG(IZ,IA),THLIFE,
     &                          CTUNIT(ITUNIT)
                       WRITE(51,1317) 
     &                          'm'//CISO(IZ),IZ,IIA,YIELD(IDET,IZ,IA),
     &                          ERROR(IDET,IZ,IA),THALF,CTUNIT(KK),
     &                          THLIFE,CTUNIT(ITUNIT)
                    ENDIF
                  ENDIF
                  WRITE(LUNOUT,*)
                  DO 22 N=1,NENERG(IZ,IA)
                     JJ = IDXIS3(N)
                     CTEMP = CDECAY(JJ)
                     CMETA = ' '
                     IF (GAMPRO(JJ).LT.0.0D0) CMETA = 'm'
                     IF (CTEMP(1:3).EQ.'e+b') CTEMP(1:3) = 'e b'
                     WRITE(LUNOUT,1006)GAMENE(JJ),ABS(GAMPRO(JJ)),CTEMP,
     &                                 CMETA
 1006                FORMAT(36X,'E_g =',F9.3,' keV  I_g =',
     &                      F10.5,2X,A4,1X,A1)
                     IF (ABS(GAMPRO(JJ)).GT.ZERO) THEN
                        NGAM = NGAM+1
                        IF (NGAM.GT.50000) STOP 'NGAM > 50000 !!'
                        EGAM(NGAM)   = GAMENE(JJ)
                        AGAM(NGAM)   = YIELD(IDET,IZ,IA)*ABS(GAMPRO(JJ))
     &                                                          /100.0D0
                        ETOGAM(NGAM) = EGAM(NGAM)*AGAM(NGAM)
                        IZGAM(NGAM)  = IZ
                        IAGAM(NGAM)  = IA
                        IDXIS5(NGAM) = NGAM
                        EGAMTO = EGAMTO+ETOGAM(NGAM)
                     ENDIF
   22             CONTINUE
                  WRITE(LUNOUT,1007) GAMTOT
 1007             FORMAT(42X,'E_g        =',F10.3,' keV/decay')
                  WRITE(LUNOUT,1008) GAMTOT*YIELD(IDET,IZ,IA)
 1008             FORMAT(42X,'E_g(total) =',E12.4,' keV')
                  WRITE(LUNOUT,*)
                  EDEPOS(NISO) = GAMTOT*YIELD(IDET,IZ,IA)
                  EDEPTO = EDEPTO+EDEPOS(NISO)
               ELSE
                  IF (IA.LE.260) THEN
                    IF (THALF.LT.1.E12) THEN
                       WRITE(LUNOUT,1001)
     &                          CISO(IZ),IZ,IIA,YIELD(IDET,IZ,IA),
     &                          CERR,ERROR(IDET,IZ,IA),THALF,CTUNIT(KK),
     &                          CDCY(NISO)
                    ELSE
                       WRITE(LUNOUT,1101)
     &                          CISO(IZ),IZ,IIA,YIELD(IDET,IZ,IA),
     &                          CERR,ERROR(IDET,IZ,IA),THALF,CTUNIT(KK),
     &                          CDCY(NISO)
                    ENDIF
                  ELSE
                    IF (THALF.LT.1.E12) THEN
                       WRITE(LUNOUT,1016)
     &                          'm'//CISO(IZ),IZ,IIA,YIELD(IDET,IZ,IA),
     &                          CERR,ERROR(IDET,IZ,IA),THALF,CTUNIT(KK),
     &                          CDCY(NISO)
                    ELSE
                       WRITE(LUNOUT,1116)
     &                          'm'//CISO(IZ),IZ,IIA,YIELD(IDET,IZ,IA),
     &                          CERR,ERROR(IDET,IZ,IA),THALF,CTUNIT(KK),
     &                          CDCY(NISO)
                    ENDIF
 1016               FORMAT(A3,2I4,E12.4,2X,A3,F5.1,' %',F15.2,A1,2X,
     &                                                              A2)
 1116               FORMAT(A3,2I4,E12.4,2X,A3,F5.1,' %',E15.5,A1,
     &                                                           2X,A2)
                  ENDIF
                  EDEPOS(NISO) = ZERO
               ENDIF
 1004          FORMAT(1X,A2,2I4,E12.4,2X,A3,F5.1,' %',F15.2,A1,2X,A2,
     &                I5,F15.2,A1)
 1104          FORMAT(1X,A2,2I4,E12.4,2X,A3,F5.1,' %',E15.5,A1,2X,A2,
     &                I5,F15.2,A1)
 1204          FORMAT(1X,A2,2I4,E12.4,2X,F5.1,2X,F15.2,A1,2X,
     &                F15.2,A1)
 1304          FORMAT(1X,A2,2I4,E12.4,2X,F5.1,2X,E15.5,A1,2X,
     &                E15.5,A1)
 1017          FORMAT(A3,2I4,E12.4,2X,A3,F5.1,' %',F15.2,A1,2X,A2,
     &                I5,F15.2,A1)
 1117          FORMAT(A3,2I4,E12.4,2X,A3,F5.1,' %',E15.5,A1,2X,A2,
     &                I5,F15.2,A1)
 1217          FORMAT(A3,2I4,E12.4,2X,F5.1,2X,F15.2,A1,2X,
     &                F15.2,A1)
 1317          FORMAT(A3,2I4,E12.4,2X,F5.1,2X,E15.5,A1,2X,
     &                E15.5,A1)
 1001          FORMAT(1X,A2,2I4,E12.4,2X,A3,F5.1,' %',F15.2,A1,2X,A2)
 1101          FORMAT(1X,A2,2I4,E12.4,2X,A3,F5.1,' %',E15.5,A1,2X,A2)

            ENDIF
    2    CONTINUE
    1 CONTINUE

      IF (TOTYLD(1,IDET).GE.1.0D-38) THEN
         TOTERR(1,IDET) = TOTERR(1,IDET)/TOTYLD(1,IDET)*100.0D0
      ELSE
         TOTERR(1,IDET) = 0.0D0
      ENDIF
      IF (TOTYLD(2,IDET).GE.1.0D-38) THEN
         TOTERR(2,IDET) = TOTERR(2,IDET)/TOTYLD(2,IDET)*100.0D0
      ELSE
         TOTERR(2,IDET) = 0.0D0
      ENDIF
      WRITE(LUNOUT,1002) TOTYLD(1,IDET),CERR,TOTERR(1,IDET),EDEPTO,
     &                   TOTYLD(2,IDET),CERR,TOTERR(2,IDET)
 1002 FORMAT(/,7X,'sum:',E12.4,2X,A3,F5.1,' %',7X,'E_total = ',E12.4,
     &       ' keV',/,11X,E12.4,2X,A3,F5.1,' %')

      CALL SORT(EDEPOS,IDXIS4,NISO,1,NISO,2)
      WRITE(LUNOUT,1010) 
 1010 FORMAT(/,/,6X,'Z   A   Tot. Yield',20X,'T1/2',8X,'E_g(total)',
     &       /,15X,'(Bq/cmc)',/)
      DO 11 I=1,NISO
         IF (EDEPOS(I).GT.ZERO) THEN
            IX = IDXIS4(I)
            IZ = IZISO(IX)
            IA = IAISO(IX)
            IIA = IA
            IF (IA.GT.260) IIA = IA-260
            DO 12 K=1,5
               JSOMER = 0
               IF (IA.GT.260) JSOMER = -1
               IF ( IZ .LT. IIA ) THEN
                  THALF = TIM1O2(IIA,IZ,KA,KZ,T12DAU,BRDAUG,JSOMER)
               ELSE
                  THALF = 0.0D0
               ENDIF
               THALF = THALF/DBLE(NSECON(K))
               KK = K
               IF (THALF.LT.100.0D0) GOTO 13
   12       CONTINUE
   13       CONTINUE
            IF ((KK.EQ.5).AND.(THALF.LT.1.0D0)) THEN
               THALF = THALF*DBLE(NSECON(5))/DBLE(NSECON(4))
               KK = 4
            ENDIF
            FRAC = EDEPOS(I)/EDEPTO*100.0D0
            CCISO = ' '
            IF (IA.GT.260) CCISO = 'm'
            IF (THALF.LT.1.E12) THEN
               WRITE(LUNOUT,1009)
     &                     CCISO,CISO(IZ),IZ,IIA,YIELD(IDET,IZ,IA),
     &                     CERR,ERROR(IDET,IZ,IA),THALF,CTUNIT(KK),
     &                     EDEPOS(I),FRAC
            ELSE
               WRITE(LUNOUT,1109)
     &                     CCISO,CISO(IZ),IZ,IIA,YIELD(IDET,IZ,IA),
     &                     CERR,ERROR(IDET,IZ,IA),THALF,CTUNIT(KK),
     &                     EDEPOS(I),FRAC
            ENDIF
 1009       FORMAT(A1,A2,2I4,E12.4,2X,A3,F5.1,' %',F15.2,A1,3X,E12.4,
     &                2X,'(',F7.3,' %)')
 1109       FORMAT(A1,A2,2I4,E12.4,2X,A3,F5.1,' %',E15.5,A1,3X,E12.4,
     &                2X,'(',F7.3,' %)')
         ENDIF
   11 CONTINUE


      CALL SORT(T12,IDXIS1,NISO,1,NISO,1)
      WRITE(LUNOUT,1000) 
      DO 5 I=1,NISO
         IX = IDXIS1(I)
         IZ = IZISO(IX)
         IA = IAISO(IX)
         IIA = IA
         IF (IA.GT.260) IIA = IA-260
         DO 6 K=1,5
            THALF = T12(I)/DBLE(NSECON(K))
            KK = K
            IF (THALF.LT.100.0D0) GOTO 7
    6    CONTINUE
    7    CONTINUE
         IF ((KK.EQ.5).AND.(THALF.LT.1.0D0)) THEN
            THALF = THALF*DBLE(NSECON(5))/DBLE(NSECON(4))
            KK = 4
         ENDIF
         CCISO = ' '
         IF (IA.GT.260) CCISO = 'm'
         IF (THALF.LT.1.E12) THEN
            WRITE(LUNOUT,1005)
     &                     CCISO,CISO(IZ),IZ,IIA,YIELD(IDET,IZ,IA),
     &                     CERR,ERROR(IDET,IZ,IA),THALF,CTUNIT(KK),
     &                     CDCY(IX),NENERG(IZ,IA)
         ELSE
            WRITE(LUNOUT,1105)
     &                     CCISO,CISO(IZ),IZ,IIA,YIELD(IDET,IZ,IA),
     &                     CERR,ERROR(IDET,IZ,IA),THALF,CTUNIT(KK),
     &                     CDCY(IX),NENERG(IZ,IA)
         ENDIF
 1005    FORMAT(A1,A2,2I4,E12.4,2X,A3,F5.1,' %',F15.2,A1,2X,A2,I5)
 1105    FORMAT(A1,A2,2I4,E12.4,2X,A3,F5.1,' %',E15.5,A1,2X,A2,I5)
    5 CONTINUE

      CALL SORT(YLD,IDXIS2,NISO,1,NISO,2)
      WRITE(LUNOUT,1000) 
      DO 8 I=1,NISO
         IX = IDXIS2(I)
         IZ = IZISO(IX)
         IA = IAISO(IX)
         IIA = IA
         IF (IA.GT.260) IIA = IA-260
         DO 9 K=1,5
            JSOMER = 0
            IF (IA.GT.260) JSOMER = -1
            IF ( IZ .LT. IIA ) THEN
               THALF = TIM1O2(IIA,IZ,KA,KZ,T12DAU,BRDAUG,JSOMER)
            ELSE
               THALF = 0.0D0
            ENDIF
            THALF = THALF/DBLE(NSECON(K))
            KK = K
            IF (THALF.LT.100.0D0) GOTO 10
    9    CONTINUE
   10    CONTINUE
         IF ((KK.EQ.5).AND.(THALF.LT.1.0D0)) THEN
            THALF = THALF*DBLE(NSECON(5))/DBLE(NSECON(4))
            KK = 4
         ENDIF
         FRAC = YIELD(IDET,IZ,IA)/TOTYLD(1,IDET)*100.0D0
         CCISO = ' '
         IF (IA.GT.260) CCISO = 'm'
         IF (THALF.LT.1.E12) THEN
            WRITE(LUNOUT,1003)
     &                     CCISO,CISO(IZ),IZ,IIA,YIELD(IDET,IZ,IA),
     &                     CERR,ERROR(IDET,IZ,IA),THALF,CTUNIT(KK),
     &                     CDCY(IX),NENERG(IZ,IA),FRAC
         ELSE
            WRITE(LUNOUT,1103)
     &                     CCISO,CISO(IZ),IZ,IIA,YIELD(IDET,IZ,IA),
     &                     CERR,ERROR(IDET,IZ,IA),THALF,CTUNIT(KK),
     &                     CDCY(IX),NENERG(IZ,IA),FRAC
         ENDIF
 1003    FORMAT(A1,A2,2I4,E12.4,2X,A3,F5.1,' %',F15.2,A1,2X,A2,
     &          I5,3X,'(',F7.3,' %)')
 1103    FORMAT(A1,A2,2I4,E12.4,2X,A3,F5.1,' %',E15.5,A1,2X,A2,
     &          I5,3X,'(',F7.3,' %)')
    8 CONTINUE

      CALL SORT(ETOGAM,IDXIS5,NGAM,1,NGAM,2)
      AGAMTO = ZERO
      AGAMSE = ZERO
      DO 14 I=1,NGAM
         IX = IDXIS5(I)
         AGAMTO = AGAMTO+AGAM(IX)
         IF (ETOGAM(I)/EGAMTO.GT.EFRAC) THEN
            AGAMSE = AGAMSE+AGAM(IX)
         ENDIF
   14 CONTINUE
      WRITE(LUNOUT,1013) 
 1013 FORMAT(/,/,' FIASCO input values:',/)
      WRITE(LUNOUT,1014) AGAMTO,EFRAC*100.0D0,AGAMSE
 1014 FORMAT(1X,'Total activity (gamma decays):',E12.4,' Bq/cmc',
     &       /,/,1X,'Total activity A (gamma decays) contributing',
     &       ' more than ',F5.1,'% to the',/,2X,'total emitted gamma',
     &       ' energy E:',E12.4,' Bq/cmc')
      WRITE(LUNOUT,1011) 
 1011 FORMAT(/,/,'   E_gamma         Z   A   Tot. Yield   E_tot',
     &       '    Fract. in E    Fract. in A',/,
     &       '    (keV)                   (Bq/cmc)    (keV)',
     &       '        (%)            (%)',/)
      DO 15 I=1,NGAM
         IF (ETOGAM(I)/EGAMTO.GT.EFRAC) THEN
            IX = IDXIS5(I)
            IF (IAGAM(IX).GT.260) THEN
               IIA = IAGAM(IX)-260
               CCISO = 'm'
            ELSE
               IIA = IAGAM(IX)
               CCISO = ' '
            ENDIF
            WRITE(LUNOUT,1012) EGAM(IX),CCISO,CISO(IZGAM(IX)),
     &                      IZGAM(IX),IIA,AGAM(IX),ETOGAM(I),
     &                      ETOGAM(I)/EGAMTO*100.0D0,
     &                      AGAM(IX)/AGAMSE*100.0D0
 1012       FORMAT(1X,F9.3,4X,A1,A2,2I4,E12.4,3X,F7.1,4X,F7.3,
     &                                                   8X,F7.3)
         ENDIF
   15 CONTINUE

      RETURN
      END
*
*===egamma=============================================================*
*
      SUBROUTINE EGAMMA(IZ,IIA,GAMENE,GAMPRO,THLIFE,CDECAY,NENERG)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE

      CHARACTER CLINE*78,CENERGY*14,CPROB*14,CDMODE*14,CLIFET*19,
     &          CISO*2,CISOTO*12,CDIG*1

      CHARACTER*4 CDECAY
      DIMENSION GAMENE(5000),GAMPRO(5000),CDECAY(5000)

      DIMENSION NSECON(5),CDIG(10)
      DATA NSECON / 1,60,3600,86400,31536000/
      DIMENSION CISO(109)
      DATA CISO /
     &    'H ','He','Li','Be','B ','C ','N ','O ','F ','Ne','Na',
     &    'Mg','Al','Si','P ','S ','Cl','Ar','K ','Ca','Sc','Ti',
     &    'V ','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As',
     &    'Se','Br','Kr','Rb','Sr','Y ','Zr','Nb','Mo','Tc','Ru',
     &    'Rh','Pd','Ag','Cd','In','Sn','Sb','Te','I ','Xe','Cs',
     &    'Ba','La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy',
     &    'Ho','Er','Tm','Yb','Lu','Hf','Ta','W ','Re','Os','Ir',
     &    'Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra',
     &    'Ac','Th','Pa','U ','Np','Pu','Am','Cm','Bk','Cf','Es',
     &    'Fm','Md','No','Lr','Rf','Ha','Sg','Ns','Hs','Mt'/
      DATA CDIG /'1','2','3','4','5','6','7','8','9','0'/

      JA = IIA
      IF (IIA.GT.260) THEN
         IA = IIA-260
      ELSE
         IA = IIA
      ENDIF

    1 CONTINUE
      NENERG = 0
      CALL OAUXFI('gamma.dat',10,'OLD',IERR)

      NISO = 0
   10 CONTINUE

      READ(10,'(A78)',END=9999) CLINE
      NISO = NISO+1
      READ(CLINE,1000) CENERGY,CPROB,CDMODE,CLIFET,CISOTO
 1000 FORMAT(3A14,A19,A12)

* extract isotope charge and mass
      I0 = 0
      DO 20 K=1,12
         I0 = I0+1
         IF (CISOTO(I0:I0).NE.' ') GOTO 21
   20 CONTINUE
      STOP ' digit not found!'
   21 CONTINUE
      I1 = I0
      DO 22 K=1,12
         I1 = I1+1
         IFOUND = 0
         DO 23 I=1,10
            IF (CISOTO(I1:I1).EQ.CDIG(I)) IFOUND = 1
   23    CONTINUE
         IF (IFOUND.EQ.0) GOTO 24
   22 CONTINUE
   24 CONTINUE
      I1 = I1-1
      READ(CISOTO(I0:I1),*,ERR=9998) IMASS
      ICHAR = 0
      IMETA = 1
      IF (CISOTO(I1+1:I1+2).EQ.'m2') THEN
         I1 = I1+2
         IMETA = -1
      ELSEIF (CISOTO(I1+1:I1+1).EQ.'m') THEN
         I1 = I1+1
         IMETA = -1
      ENDIF
      DO 11 K=1,109
         IF (CISOTO(I1+1:I1+2).EQ.CISO(K)) THEN
            ICHAR = K
            GOTO 12
         ENDIF
   11 CONTINUE
   12 CONTINUE
      IF ((ICHAR.EQ.0).OR.(IMASS.LE.0).OR.(IMASS.GT.260)) THEN
         WRITE(*,*) ' Isotope not found: ',CISOTO,ICHAR,IMASS,NISO
         WRITE(*,*) '                    ',JA,IMETA
         IF ((JA.GT.260).AND.(IMETA.EQ.-1)) THEN
            JA = JA-260
            GOTO 1
         ENDIF
         STOP
      ENDIF

      IF ((IMASS-IA).GT.70) GOTO 9999

      IF (((JA.LE.260).AND.(IMETA.EQ.-1)).OR.
     &    ((JA.GT.260).AND.(IMETA.NE.-1))) GOTO 10

      IF ((ICHAR.EQ.IZ).AND.(IMASS.EQ.IA)) THEN

         NENERG = NENERG+1

* gamma energy and transition probability
         READ(CENERGY,*,ERR=9998) E
         IF (CPROB.NE.'              ') READ(CPROB,*,ERR=9998) P
         GAMENE(NENERG) = E
         GAMPRO(NENERG) = P*DBLE(IMETA)

* decay mode
         I0 = 1
         IF (CDMODE.NE.'              ') THEN
            DO 13 K=1,14
               IF (CDMODE(I0:I0).EQ.' ') THEN
                  I0 = I0+1
               ELSE
                  GOTO 14
               ENDIF
   13       CONTINUE
   14       CONTINUE
            CDECAY(NENERG) = CDMODE(I0:I0+3)
         ENDIF

* half life
         IUNIT = -1
         IF (CLIFET.NE.'                   ') THEN
            READ(CLIFET,*,ERR=9998) T,IUNIT
         ENDIF
         IF ((IUNIT.GE.0).AND.(IUNIT.LE.4)) THEN
            THLIFE = T*DBLE(NSECON(IUNIT+1))
         ELSE
            WRITE(*,*) ' life-time unit not found: ',ICHAR,IMASS
            STOP
         ENDIF

      ENDIF

      GOTO 10

 9998 CONTINUE
      WRITE(*,*) ' inconsistent format ',NISO

 9999 CONTINUE
      CLOSE(10)

      RETURN
      END
*
*===sort===============================================================*
*
      SUBROUTINE SORT(A,IDX,N,I0,I1,MODE)

************************************************************************
* This subroutine sorts entries in A in increasing/decreasing order    *
* of A(i).                                                             *
*              MODE  = 1     increasing in A(i=1..N)                   *
*                    = 2     decreasing in A(i=1..N)                   *
* This version dated 21.04.95 is revised by S. Roesler                 *
************************************************************************

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE

      DIMENSION A(N),IDX(N)

      M = I1
   10 CONTINUE
      M = I1-1
      IF (M.LE.0) RETURN
      L = 0
      DO 20 I=I0,M
         J = I+1
         IF (MODE.EQ.1) THEN
            IF (A(I).LE.A(J)) GOTO 20
         ELSE
            IF (A(I).GE.A(J)) GOTO 20
         ENDIF
         B    = A(I)
         A(I) = A(J)
         A(J) = B
         IX     = IDX(I)
         IDX(I) = IDX(J)
         IDX(J) = IX
         L = 1
   20 CONTINUE
      IF (L.EQ.1) GOTO 10

      RETURN
      END
*
*===analys=============================================================*
*
      SUBROUTINE ANALYS(RESNUC,RESERR,VOL,IDET,MODE)

************************************************************************
* Analysis routine for residual nuclei output.                         *
* This version dated 16.01.01 is written by S. Roesler                 *
************************************************************************

C     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE

      LOGICAL LFIRST

      DIMENSION RESNUC (100,260), RESERR (100,260)

      DATA LFIRST /.TRUE./

      IF (LFIRST) THEN
         OPEN(99,FILE='analysis.out',STATUS='UNKNOWN')
         LFIRST = .FALSE.
      ENDIF

      WRITE(99,1000) IDET,VOL,RESNUC(1,3),RESERR(1,3),MODE
 1000 FORMAT(I5,F12.0,1P,2(E12.2),I5)

      RETURN
      END
      SUBROUTINE TEST

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
      SAVE

      DIMENSION RBIN(100),EBIN(100)
      DATA RBIN,EBIN /100*0.0D0,100*0.0D0/

      DIMENSION JZ(10),JA(10),JISO(10)
      DATA   JZ /  1,  4,  6, 21, 21, 23, 23, 80, 83, 18/
      DATA   JA /  3,  7, 11, 44, 44, 49, 52,192,207, 39/
      DATA JISO /  1,  1,  1,  1,  2,  1,  1,  1,  1,  1/

      EG     = -ONEONE
      WEIGHT = ZERZER
      DO 10 I=1,10
C        CALL READGL(JZ(I),JA(I),JISO(I),EG,WEIGHT,IBIN)
         CALL READGL(JZ(I),JA(I),JISO(I),EG,WEIGHT)
         WRITE(*,*) JZ(I),JA(I),JISO(I),EG,WEIGHT
   10 CONTINUE
      STOP

      IZ  = 80
      IA  = 192
      ISO = 1
      EG  = ZERZER
      WEIGHT = ZERZER

      NEVENT = 1000000

      IBIN = 1
      DO 1 I=1,NEVENT
C        CALL READGL(IZ,IA,ISO,EG,WEIGHT,IBIN)
         CALL READGL(IZ,IA,ISO,EG,WEIGHT)
         RBIN(IBIN) = RBIN(IBIN)+1.0D0
         EBIN(IBIN) = EBIN(IBIN)+EG
    1 CONTINUE
      DO 2 I=1,100
         EBIN(I) = EBIN(I)/MAX(ONEONE,RBIN(I))
         RBIN(I) = RBIN(I)/DBLE(NEVENT)*WEIGHT*100.0D0
         WRITE(*,*) I,EBIN(I),RBIN(I)
    2 CONTINUE
      END
*
*===readgl=============================================================*
*
C     SUBROUTINE READGL(KZ,KA,KISO,EG,WEIGHT,IBIN)
      SUBROUTINE READGL(KZ,KA,KISO,EG,WEIGHT)

************************************************************************
* Sampling of gamma-ray energy for the radioactive decay of a given    *
* gamma-emitting isotope.                                              *
*                                                                      *
* Input:    KZ, KA     charge/mass number of isotope                   *
*           KISO       flag for isomeric state                         *
* Output:   EG         energy of sampled gamma ray (GeV)               *
*           WEIGHT     weigth of the gamma ray                         *
*                      =total emission probability ( sum[Ig] )         *
*                                                                      *
* If input EG < 0: just check if this is a gamma ray emitter without   *
*                  sampling. If yes: WEIGHT = 1, otherwise WEIGHT = 0  *
*                                                                      *
* This version dated 5.8.02    is written by S. Roesler.               *
************************************************************************

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
      SAVE

      LOGICAL LFIRST

      PARAMETER (LUNGAL = 30,
     &           LOUT   = 6)
C    &           LOUT   = LUNOUT)

      PARAMETER (IZMAX = 100,
     &           IAMAX = 257,
     &           NMZMIN = -5,
     &           NMZMAX = 157,
     &           NGLMAX = 68617)

      DIMENSION GAMENE(NGLMAX),GAMPRO(NGLMAX),
     &          IDXGL(2,IZMAX,NMZMIN:NMZMAX),NGL(2,IZMAX,NMZMIN:NMZMAX),
     &          GAMTOT(2,IZMAX,NMZMIN:NMZMAX)

      DATA LFIRST /.TRUE./
*
*----------------------------------------------------------------------*
* Initializations
*
      IF (LFIRST) THEN
*
*   zero arrays
         DO 1 I=1,NGLMAX
            GAMENE(I) = ZERZER
            GAMPRO(I) = ZERZER
    1    CONTINUE
         DO 2 IZ=1,IZMAX
            DO 3 NMZ=NMZMIN,NMZMAX
               IDXGL(1,IZ,NMZ)  = 0
               IDXGL(2,IZ,NMZ)  = 0
               NGL(1,IZ,NMZ)    = 0
               NGL(2,IZ,NMZ)    = 0
               GAMTOT(1,IZ,NMZ) = ZERZER
               GAMTOT(2,IZ,NMZ) = ZERZER
    3       CONTINUE
    2    CONTINUE
*
*   open data file
         CALL OAUXFI('glines.dat',LUNGAL,'OLD',IERR)
         IF (IERR.GT.0) STOP 'Cannot open gamma library !'
*
*   read data sets
*
*     Iz, Ia    mass number and charge of isotope
*     In, Nmz   number of neutrons, neutron excess
*     Iso       isomer flag (=1 ground state, =2 isomer)
*
*     Nl, Ngl(Iso,Iz,Nmz)       number of emitted gamma lines
*     Prob, Gamtot(Iso,Iz,Nmz)  total emission probability ( sum[Ig] )
*     Idxgl(Iso,Iz,Nmz)         pointer to the first gamma line in
*                               Gamene and Gampro
*     Gamene(i)                 gamma energy in keV
*     Gampro(i)                 cumulative emission probability
*     Etot                      energy emitted per decay ( sum[Eg*Ig] )
*
         KMAX = 0
    4    CONTINUE
         READ(LUNGAL,1000,END=10) IZ,IA,NL,ETOT,PROB
 1000    FORMAT(3I4,1P,2E15.5)
         K0  = KMAX+1
         ISO = 1
         IF (IA.GT.260) THEN
            ISO = 2
            IA  = IA-260
         ENDIF
         IN  = IA-IZ
         NMZ = IN-IZ
         IF ((NMZ.LT.NMZMIN).OR.(NMZ.GT.NMZMAX)) THEN
            WRITE(LOUT,*) ' N-Z out of range ! ',IZ,IA
            STOP
         ENDIF
         IDXGL(ISO,IZ,NMZ)  = K0
         NGL(ISO,IZ,NMZ)    = NL
         GAMTOT(ISO,IZ,NMZ) = PROB
         NLINES = NL/10
         IF (NLINES.EQ.0) THEN
            I0 = K0
            I1 = I0+NL-1
            READ(LUNGAL,1001) (GAMENE(I),I=I0,I1)
 1001       FORMAT(10F7.1)
         ELSE
            DO 5 ILINE=1,NLINES
               I0 = K0+(ILINE-1)*10
               I1 = I0+9
               READ(LUNGAL,1001) (GAMENE(I),I=I0,I1)
    5       CONTINUE
            NLEFT = NL-10*NLINES
            IF (NLEFT.GT.0) THEN
               I0 = I1+1
               I1 = I0+NLEFT-1
               READ(LUNGAL,1001) (GAMENE(I),I=I0,I1)
            ENDIF
         ENDIF
         IF (NLINES.EQ.0) THEN
            I0 = K0
            I1 = I0+NL-1
            READ(LUNGAL,1002) (GAMPRO(I),I=I0,I1)
 1002       FORMAT(1P,10E9.2)
         ELSE
            DO 6 ILINE=1,NLINES
               I0 = K0+(ILINE-1)*10
               I1 = I0+9
               READ(LUNGAL,1002) (GAMPRO(I),I=I0,I1)
    6       CONTINUE
            NLEFT = NL-10*NLINES
            IF (NLEFT.GT.0) THEN
               I0 = I1+1
               I1 = I0+NLEFT-1
               READ(LUNGAL,1002) (GAMPRO(I),I=I0,I1)
            ENDIF
         ENDIF
         KMAX = I1
         GOTO 4
   10    CONTINUE

         WRITE(LOUT,*) ' READGL: Gamma-line data base read.'
         CLOSE(LUNGAL)
         LFIRST = .FALSE.

***temporary
C        CALL FL48UT (ISRM48,ISEED1,ISEED2)
C        CALL FL48IN (54217137,ISEED1,ISEED2)
***

      ENDIF
*
*----------------------------------------------------------------------*
* Check if this is a gamma ray emitter (without sampling)
*
      WEIGHT = ZERZER
      IF (EG.LT.ZERZER) THEN
         KN  = KA-KZ
         NMZ = KN-KZ
         IF ((NMZ.LT.NMZMIN).OR.(NMZ.GT.NMZMAX)) THEN
            WRITE(LOUT,*) ' READGL: N-Z out of range ! ',NMZ
            RETURN
         ENDIF
         IF ((KISO.LT.1).OR.(KISO.GT.2)) THEN
            WRITE(LOUT,*) ' READGL: Inconsistent isomer flag ! ',KISO
            RETURN
         ENDIF
         IF (GAMTOT(KISO,KZ,NMZ).GT.ZERZER) WEIGHT = ONEONE
         RETURN
      ELSE
*
*----------------------------------------------------------------------*
* Sample gamma ray energy for given isotope
*
         EG  = ZERZER
         KN  = KA-KZ
         NMZ = KN-KZ
         IF ((NMZ.LT.NMZMIN).OR.(NMZ.GT.NMZMAX)) THEN
            WRITE(LOUT,*) ' READGL: N-Z out of range ! ',NMZ
            RETURN
         ENDIF
         IF ((KISO.LT.1).OR.(KISO.GT.2)) THEN
            WRITE(LOUT,*) ' READGL: Inconsistent isomer flag ! ',KISO
            RETURN
         ENDIF
         IF (GAMTOT(KISO,KZ,NMZ).LE.ZERZER) THEN
            WEIGHT = ZERZER
C           WRITE(LOUT,*) ' READGL: This is not a gamma ray emitter ! '
            RETURN
         ENDIF
         WEIGHT = GAMTOT(KISO,KZ,NMZ)
         IF (WEIGHT.GT.ZERZER) THEN
            K0 = IDXGL(KISO,KZ,NMZ)
            K1 = IDXGL(KISO,KZ,NMZ)+NGL(KISO,KZ,NMZ)-1
            FRAC = FLRNDM(WEIGHT)
            DO 20 I=K0,K1
               IF (FRAC.LT.(ONEONE-GAMPRO(I))) THEN
                  EG   = 1.0D-6*GAMENE(I)
C                 IBIN = I-IDXGL(KISO,KZ,NMZ)+1
                  RETURN
               ENDIF
   20       CONTINUE
            WRITE(LOUT,*) 
     &         ' READGL: Inconsistent emission probabilities !',
     &         KZ,KA,KISO,GAMTOT(KISO,KZ,NMZ),NGL(KISO,KZ,NMZ),
     &         GAMPRO(K1)
            STOP
         ENDIF
      ENDIF
*
*----------------------------------------------------------------------*
* Online consistency check
*
C 100 CONTINUE
C     WRITE(LOUT,*) ' IZ, IA = '
C     READ(*,*) IZ,IA
C     WRITE(LOUT,*) ' Isomer = '
C     READ(*,*) ISO
C     IN  = IA-IZ
C     NMZ = IN-IZ
C     IF ((NMZ.LT.NMZMIN).OR.(NMZ.GT.NMZMAX)) THEN
C        WRITE(LOUT,*) ' N-Z out of range: ',NMZ
C        GOTO 100
C     ENDIF
C     IF ((ISO.LT.1).OR.(ISO.GT.2)) GOTO 100

C     I0   = IDXGL(ISO,IZ,NMZ)
C     NL   = NGL(ISO,IZ,NMZ)
C     PROB = GAMTOT(ISO,IZ,NMZ)
C     WRITE(LOUT,*) ' Number of lines:            ',NL
C     WRITE(LOUT,*) ' Total emission probability: ',PROB
C     DO 21 I=I0,I0+NL-1
C        IF (I.EQ.I0) THEN
C           WRITE(LOUT,'(3X,F7.1,F10.5)') 
C    &         GAMENE(I),(ONEONE-GAMPRO(I))*PROB*100.0D0
C           WRITE(LOUT,'(3X,F7.1,1P,E15.5)') 
C    &         GAMENE(I),(ONEONE-GAMPRO(I))
C        ELSE
C           WRITE(LOUT,'(3X,F7.1,F10.5)') 
C    &         GAMENE(I),(GAMPRO(I-1)-GAMPRO(I))*PROB*100.0D0
C           WRITE(LOUT,'(3X,F7.1,1P,E15.5)')
C    &         GAMENE(I),(ONEONE-GAMPRO(I))
C        ENDIF
C  21 CONTINUE

C     GOTO 100

      RETURN
      END
*
*===readpl=============================================================*
*
      SUBROUTINE READPL(KZ,KA,KISO,EP,WEIGHT)

************************************************************************
* Sampling of positron energy for the radioactive decay of a given     *
* positron-emitting isotope.                                           *
*                                                                      *
* Input:    KZ, KA     charge/mass number of isotope                   *
*           KISO       flag for isomeric state                         *
* Output:   EP         energy of sampled positron (GeV)                *
*           WEIGHT     weigth of the positron                          *
*                      =total emission probability ( sum[Ip] )         *
*                                                                      *
* If input EP < 0: just check if this is a positron emitter without    *
*                  sampling. If yes: WEIGHT = 1, otherwise WEIGHT = 0  *
*                                                                      *
* This version dated 17.1.03   is written by S. Roesler.               *
************************************************************************

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
      SAVE

      LOGICAL LFIRST

      PARAMETER (LUNPOL = 30,
     &           LOUT   = 6)
C    &           LOUT   = LUNOUT)

      PARAMETER (IZMAX = 100,
     &           IAMAX = 257,
     &           NMZMIN = -5,
     &           NMZMAX = 157,
     &           NPLMAX = 68617)

      DIMENSION POSENE(NPLMAX),POSPRO(NPLMAX),
     &          IDXPL(2,IZMAX,NMZMIN:NMZMAX),NPL(2,IZMAX,NMZMIN:NMZMAX),
     &          POSTOT(2,IZMAX,NMZMIN:NMZMAX)

      DATA LFIRST /.TRUE./
*
*----------------------------------------------------------------------*
* Initializations
*
      IF (LFIRST) THEN
*
*   zero arrays
         DO 1 I=1,NPLMAX
            POSENE(I) = ZERZER
            POSPRO(I) = ZERZER
    1    CONTINUE
         DO 2 IZ=1,IZMAX
            DO 3 NMZ=NMZMIN,NMZMAX
               IDXPL(1,IZ,NMZ)  = 0
               IDXPL(2,IZ,NMZ)  = 0
               NPL(1,IZ,NMZ)    = 0
               NPL(2,IZ,NMZ)    = 0
               POSTOT(1,IZ,NMZ) = ZERZER
               POSTOT(2,IZ,NMZ) = ZERZER
    3       CONTINUE
    2    CONTINUE
*
*   open data file
         CALL OAUXFI('plines.dat',LUNPOL,'OLD',IERR)
         IF (IERR.GT.0) STOP 'Cannot open positron library !'
*
*   read data sets
*
*     Iz, Ia    mass number and charge of isotope
*     In, Nmz   number of neutrons, neutron excess
*     Iso       isomer flag (=1 ground state, =2 isomer)
*
*     Nl, Ngl(Iso,Iz,Nmz)       number of emitted positron lines
*     Prob, Postot(Iso,Iz,Nmz)  total emission probability ( sum[Ig] )
*     Idxgl(Iso,Iz,Nmz)         pointer to the first positron line in
*                               Posene and Pospro
*     Posene(i)                 positron energy in keV
*     Pospro(i)                 cumulative emission probability
*     Etot                      energy emitted per decay ( sum[Eg*Ig] )
*
         KMAX = 0
    4    CONTINUE
         READ(LUNPOL,1000,END=10) IZ,IA,NL,ETOT,PROB
 1000    FORMAT(3I4,1P,2E15.5)
         K0  = KMAX+1
         ISO = 1
         IF (IA.GT.260) THEN
            ISO = 2
            IA  = IA-260
         ENDIF
         IN  = IA-IZ
         NMZ = IN-IZ
         IF ((NMZ.LT.NMZMIN).OR.(NMZ.GT.NMZMAX)) THEN
            WRITE(LOUT,*) ' N-Z out of range ! ',IZ,IA
            STOP
         ENDIF
         IDXPL(ISO,IZ,NMZ)  = K0
         NPL(ISO,IZ,NMZ)    = NL
         POSTOT(ISO,IZ,NMZ) = PROB
         NLINES = NL/10
         IF (NLINES.EQ.0) THEN
            I0 = K0
            I1 = I0+NL-1
            READ(LUNPOL,1001) (POSENE(I),I=I0,I1)
 1001       FORMAT(10F7.1)
         ELSE
            DO 5 ILINE=1,NLINES
               I0 = K0+(ILINE-1)*10
               I1 = I0+9
               READ(LUNPOL,1001) (POSENE(I),I=I0,I1)
    5       CONTINUE
            NLEFT = NL-10*NLINES
            IF (NLEFT.GT.0) THEN
               I0 = I1+1
               I1 = I0+NLEFT-1
               READ(LUNPOL,1001) (POSENE(I),I=I0,I1)
            ENDIF
         ENDIF
         IF (NLINES.EQ.0) THEN
            I0 = K0
            I1 = I0+NL-1
            READ(LUNPOL,1002) (POSPRO(I),I=I0,I1)
 1002       FORMAT(1P,10E9.2)
         ELSE
            DO 6 ILINE=1,NLINES
               I0 = K0+(ILINE-1)*10
               I1 = I0+9
               READ(LUNPOL,1002) (POSPRO(I),I=I0,I1)
    6       CONTINUE
            NLEFT = NL-10*NLINES
            IF (NLEFT.GT.0) THEN
               I0 = I1+1
               I1 = I0+NLEFT-1
               READ(LUNPOL,1002) (POSPRO(I),I=I0,I1)
            ENDIF
         ENDIF
         KMAX = I1
         GOTO 4
   10    CONTINUE

         WRITE(LOUT,*) ' READPL: Positron-line data base read.'
         CLOSE(LUNPOL)
         LFIRST = .FALSE.

***temporary
C        CALL FL48UT (ISRM48,ISEED1,ISEED2)
C        CALL FL48IN (54217137,ISEED1,ISEED2)
***

      ENDIF
*
*----------------------------------------------------------------------*
* Check if this is a positron emitter (without sampling)
*
      WEIGHT = ZERZER
      IF (EP.LT.ZERZER) THEN
         KN  = KA-KZ
         NMZ = KN-KZ
         IF ((NMZ.LT.NMZMIN).OR.(NMZ.GT.NMZMAX)) THEN
            WRITE(LOUT,*) ' READPL: N-Z out of range ! ',NMZ
            RETURN
         ENDIF
         IF ((KISO.LT.1).OR.(KISO.GT.2)) THEN
            WRITE(LOUT,*) ' READPL: Inconsistent isomer flag ! ',KISO
            RETURN
         ENDIF
         IF (POSTOT(KISO,KZ,NMZ).GT.ZERZER) WEIGHT = ONEONE
         RETURN
      ELSE
*
*----------------------------------------------------------------------*
* Sample positron energy for given isotope
*
         EP  = ZERZER
         KN  = KA-KZ
         NMZ = KN-KZ
         IF ((NMZ.LT.NMZMIN).OR.(NMZ.GT.NMZMAX)) THEN
            WRITE(LOUT,*) ' READPL: N-Z out of range ! ',NMZ
            RETURN
         ENDIF
         IF ((KISO.LT.1).OR.(KISO.GT.2)) THEN
            WRITE(LOUT,*) ' READPL: Inconsistent isomer flag ! ',KISO
            RETURN
         ENDIF
         IF (POSTOT(KISO,KZ,NMZ).LE.ZERZER) THEN
            WEIGHT = ZERZER
C           WRITE(LOUT,*) ' READPL: This is not a positron emitter ! '
            RETURN
         ENDIF
         WEIGHT = POSTOT(KISO,KZ,NMZ)
         IF (WEIGHT.GT.ZERZER) THEN
            K0 = IDXPL(KISO,KZ,NMZ)
            K1 = IDXPL(KISO,KZ,NMZ)+NPL(KISO,KZ,NMZ)-1
            FRAC = FLRNDM(WEIGHT)
            DO 20 I=K0,K1
               IF (FRAC.LT.(ONEONE-POSPRO(I))) THEN
                  EP   = 1.0D-6*POSENE(I)
C                 IBIN = I-IDXPL(KISO,KZ,NMZ)+1
                  RETURN
               ENDIF
   20       CONTINUE
            WRITE(LOUT,*) 
     &         ' READPL: Inconsistent emission probabilities !',
     &         KZ,KA,KISO,POSTOT(KISO,KZ,NMZ),NPL(KISO,KZ,NMZ),
     &         POSPRO(K1)
            STOP
         ENDIF
      ENDIF
*
*----------------------------------------------------------------------*
* Online consistency check
*
C 100 CONTINUE
C     WRITE(LOUT,*) ' IZ, IA = '
C     READ(*,*) IZ,IA
C     WRITE(LOUT,*) ' Isomer = '
C     READ(*,*) ISO
C     IN  = IA-IZ
C     NMZ = IN-IZ
C     IF ((NMZ.LT.NMZMIN).OR.(NMZ.GT.NMZMAX)) THEN
C        WRITE(LOUT,*) ' N-Z out of range: ',NMZ
C        GOTO 100
C     ENDIF
C     IF ((ISO.LT.1).OR.(ISO.GT.2)) GOTO 100

C     I0   = IDXPL(ISO,IZ,NMZ)
C     NL   = NPL(ISO,IZ,NMZ)
C     PROB = POSTOT(ISO,IZ,NMZ)
C     WRITE(LOUT,*) ' Number of lines:            ',NL
C     WRITE(LOUT,*) ' Total emission probability: ',PROB
C     DO 21 I=I0,I0+NL-1
C        IF (I.EQ.I0) THEN
C           WRITE(LOUT,'(3X,F7.1,F10.5)') 
C    &         POSENE(I),(ONEONE-POSPRO(I))*PROB*100.0D0
C           WRITE(LOUT,'(3X,F7.1,1P,E15.5)') 
C    &         POSENE(I),(ONEONE-POSPRO(I))
C        ELSE
C           WRITE(LOUT,'(3X,F7.1,F10.5)') 
C    &         POSENE(I),(POSPRO(I-1)-POSPRO(I))*PROB*100.0D0
C           WRITE(LOUT,'(3X,F7.1,1P,E15.5)')
C    &         POSENE(I),(ONEONE-POSPRO(I))
C        ENDIF
C  21 CONTINUE

C     GOTO 100

      RETURN
      END
*
*===readel=============================================================*
*
      SUBROUTINE READEL(KZ,KA,KISO,EE,WEIGHT)

************************************************************************
* Sampling of electron energy for the radioactive decay of a given     *
* electron-emitting isotope.                                           *
*                                                                      *
* Input:    KZ, KA     charge/mass number of isotope                   *
*           KISO       flag for isomeric state                         *
* Output:   EE         energy of sampled electron (GeV)                *
*           WEIGHT     weigth of the electron                          *
*                      =total emission probability ( sum[Ip] )         *
*                                                                      *
* If input EE < 0: just check if this is a electron emitter without    *
*                  sampling. If yes: WEIGHT = 1, otherwise WEIGHT = 0  *
*                                                                      *
* This version dated 27.1.03   is written by S. Roesler.               *
************************************************************************

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
      SAVE

      LOGICAL LFIRST

      PARAMETER (LUNELL = 30,
     &           LOUT   = 6)
C    &           LOUT   = LUNOUT)

      PARAMETER (IZMAX = 100,
     &           IAMAX = 257,
     &           NMZMIN = -5,
     &           NMZMAX = 157,
     &           NELMAX = 68617)

      DIMENSION ELCENE(NELMAX),ELCPRO(NELMAX),
     &          IDXEL(2,IZMAX,NMZMIN:NMZMAX),NEL(2,IZMAX,NMZMIN:NMZMAX),
     &          ELCTOT(2,IZMAX,NMZMIN:NMZMAX)

      DATA LFIRST /.TRUE./
*
*----------------------------------------------------------------------*
* Initializations
*
      IF (LFIRST) THEN
*
*   zero arrays
         DO 1 I=1,NELMAX
            ELCENE(I) = ZERZER
            ELCPRO(I) = ZERZER
    1    CONTINUE
         DO 2 IZ=1,IZMAX
            DO 3 NMZ=NMZMIN,NMZMAX
               IDXEL(1,IZ,NMZ)  = 0
               IDXEL(2,IZ,NMZ)  = 0
               NEL(1,IZ,NMZ)    = 0
               NEL(2,IZ,NMZ)    = 0
               ELCTOT(1,IZ,NMZ) = ZERZER
               ELCTOT(2,IZ,NMZ) = ZERZER
    3       CONTINUE
    2    CONTINUE
*
*   open data file
         CALL OAUXFI('eiclines38.dat',LUNELL,'OLD',IERR)
         IF (IERR.GT.0) STOP 'Cannot open electron library !'
*
*   read data sets
*
*     Iz, Ia    mass number and charge of isotope
*     In, Nmz   number of neutrons, neutron excess
*     Iso       isomer flag (=1 ground state, =2 isomer)
*
*     Nl, Ngl(Iso,Iz,Nmz)       number of emitted electron lines
*     Prob, Postot(Iso,Iz,Nmz)  total emission probability ( sum[Ig] )
*     Idxgl(Iso,Iz,Nmz)         pointer to the first electron line in
*                               Posene and Pospro
*     Posene(i)                 electron energy in keV
*     Pospro(i)                 cumulative emission probability
*     Etot                      energy emitted per decay ( sum[Eg*Ig] )
*
         KMAX = 0
    4    CONTINUE
         READ(LUNELL,1000,END=10) IZ,IA,NL,ETOT,PROB
 1000    FORMAT(3I4,1P,2E15.5)
         K0  = KMAX+1
         ISO = 1
         IF (IA.GT.260) THEN
            ISO = 2
            IA  = IA-260
         ENDIF
         IN  = IA-IZ
         NMZ = IN-IZ
         IF ((NMZ.LT.NMZMIN).OR.(NMZ.GT.NMZMAX)) THEN
            WRITE(LOUT,*) ' N-Z out of range ! ',IZ,IA
            STOP
         ENDIF
         IDXEL(ISO,IZ,NMZ)  = K0
         NEL(ISO,IZ,NMZ)    = NL
         ELCTOT(ISO,IZ,NMZ) = PROB
         NLINES = NL/10
         IF (NLINES.EQ.0) THEN
            I0 = K0
            I1 = I0+NL-1
            READ(LUNELL,1001) (ELCENE(I),I=I0,I1)
 1001       FORMAT(10F7.1)
         ELSE
            DO 5 ILINE=1,NLINES
               I0 = K0+(ILINE-1)*10
               I1 = I0+9
               READ(LUNELL,1001) (ELCENE(I),I=I0,I1)
    5       CONTINUE
            NLEFT = NL-10*NLINES
            IF (NLEFT.GT.0) THEN
               I0 = I1+1
               I1 = I0+NLEFT-1
               READ(LUNELL,1001) (ELCENE(I),I=I0,I1)
            ENDIF
         ENDIF
         IF (NLINES.EQ.0) THEN
            I0 = K0
            I1 = I0+NL-1
            READ(LUNELL,1002) (ELCPRO(I),I=I0,I1)
 1002       FORMAT(1P,10E9.2)
         ELSE
            DO 6 ILINE=1,NLINES
               I0 = K0+(ILINE-1)*10
               I1 = I0+9
               READ(LUNELL,1002) (ELCPRO(I),I=I0,I1)
    6       CONTINUE
            NLEFT = NL-10*NLINES
            IF (NLEFT.GT.0) THEN
               I0 = I1+1
               I1 = I0+NLEFT-1
               READ(LUNELL,1002) (ELCPRO(I),I=I0,I1)
            ENDIF
         ENDIF
         KMAX = I1
         GOTO 4
   10    CONTINUE

         WRITE(LOUT,*) ' READEL: Electron-line data base read.'
         CLOSE(LUNELL)
         LFIRST = .FALSE.

***temporary
C        CALL FL48UT (ISRM48,ISEED1,ISEED2)
C        CALL FL48IN (54217137,ISEED1,ISEED2)
***

      ENDIF
*
*----------------------------------------------------------------------*
* Check if this is a electron emitter (without sampling)
*
      WEIGHT = ZERZER
      IF (EE.LT.ZERZER) THEN
         KN  = KA-KZ
         NMZ = KN-KZ
         IF ((NMZ.LT.NMZMIN).OR.(NMZ.GT.NMZMAX)) THEN
            WRITE(LOUT,*) ' READEL: N-Z out of range ! ',NMZ
            RETURN
         ENDIF
         IF ((KISO.LT.1).OR.(KISO.GT.2)) THEN
            WRITE(LOUT,*) ' READEL: Inconsistent isomer flag ! ',KISO
            RETURN
         ENDIF
         IF (ELCTOT(KISO,KZ,NMZ).GT.ZERZER) WEIGHT = ONEONE
         RETURN
      ELSE
*
*----------------------------------------------------------------------*
* Sample electron energy for given isotope
*
         EE  = ZERZER
         KN  = KA-KZ
         NMZ = KN-KZ
         IF ((NMZ.LT.NMZMIN).OR.(NMZ.GT.NMZMAX)) THEN
            WRITE(LOUT,*) ' READEL: N-Z out of range ! ',NMZ
            RETURN
         ENDIF
         IF ((KISO.LT.1).OR.(KISO.GT.2)) THEN
            WRITE(LOUT,*) ' READEL: Inconsistent isomer flag ! ',KISO
            RETURN
         ENDIF
         IF (ELCTOT(KISO,KZ,NMZ).LE.ZERZER) THEN
            WEIGHT = ZERZER
C           WRITE(LOUT,*) ' READEL: This is not a electron emitter ! '
            RETURN
         ENDIF
         WEIGHT = ELCTOT(KISO,KZ,NMZ)
         IF (WEIGHT.GT.ZERZER) THEN
            K0 = IDXEL(KISO,KZ,NMZ)
            K1 = IDXEL(KISO,KZ,NMZ)+NEL(KISO,KZ,NMZ)-1
            FRAC = FLRNDM(WEIGHT)
            DO 20 I=K0,K1
               IF (FRAC.LT.(ONEONE-ELCPRO(I))) THEN
                  EE   = 1.0D-6*ELCENE(I)
C                 IBIN = I-IDXEL(KISO,KZ,NMZ)+1
                  RETURN
               ENDIF
   20       CONTINUE
            WRITE(LOUT,*) 
     &         ' READEL: Inconsistent emission probabilities !',
     &         KZ,KA,KISO,ELCTOT(KISO,KZ,NMZ),NEL(KISO,KZ,NMZ),
     &         ELCPRO(K1)
            STOP
         ENDIF
      ENDIF
*
*----------------------------------------------------------------------*
* Online consistency check
*
C 100 CONTINUE
C     WRITE(LOUT,*) ' IZ, IA = '
C     READ(*,*) IZ,IA
C     WRITE(LOUT,*) ' Isomer = '
C     READ(*,*) ISO
C     IN  = IA-IZ
C     NMZ = IN-IZ
C     IF ((NMZ.LT.NMZMIN).OR.(NMZ.GT.NMZMAX)) THEN
C        WRITE(LOUT,*) ' N-Z out of range: ',NMZ
C        GOTO 100
C     ENDIF
C     IF ((ISO.LT.1).OR.(ISO.GT.2)) GOTO 100

C     I0   = IDXEL(ISO,IZ,NMZ)
C     NL   = NEL(ISO,IZ,NMZ)
C     PROB = ELCTOT(ISO,IZ,NMZ)
C     WRITE(LOUT,*) ' Number of lines:            ',NL
C     WRITE(LOUT,*) ' Total emission probability: ',PROB
C     DO 21 I=I0,I0+NL-1
C        IF (I.EQ.I0) THEN
C           WRITE(LOUT,'(3X,F7.1,F10.5)') 
C    &         ELCENE(I),(ONEONE-ELCPRO(I))*PROB*100.0D0
C           WRITE(LOUT,'(3X,F7.1,1P,E15.5)') 
C    &         ELCENE(I),(ONEONE-ELCPRO(I))
C        ELSE
C           WRITE(LOUT,'(3X,F7.1,F10.5)') 
C    &         ELCENE(I),(ELCPRO(I-1)-ELCPRO(I))*PROB*100.0D0
C           WRITE(LOUT,'(3X,F7.1,1P,E15.5)')
C    &         ELCENE(I),(ONEONE-ELCPRO(I))
C        ENDIF
C  21 CONTINUE

C     GOTO 100

      RETURN
      END
