* *----------------------------------------------------------------------* * * * 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