diff options
Diffstat (limited to 'Dragon/src/EDILUM.f')
| -rw-r--r-- | Dragon/src/EDILUM.f | 433 |
1 files changed, 433 insertions, 0 deletions
diff --git a/Dragon/src/EDILUM.f b/Dragon/src/EDILUM.f new file mode 100644 index 0000000..bbf8c36 --- /dev/null +++ b/Dragon/src/EDILUM.f @@ -0,0 +1,433 @@ +*DECK EDILUM + SUBROUTINE EDILUM(IPRINT,IPEDIT,MAXFP,NBISO,NBFISS,NBDPF,NSUPS, + & NREAC,NFATH,NBCH,HICH,HISO,MYLIST,HREAC,IDREAC,DENER,DDECA, + & IPREAC,PRATE,YIELD,LISO,NBFISS2,NBFPCH2) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Complete and lump the burnup chain from NBISO to NBCH isotopes. +* Write the lumped chain on the LCM object. Based on subroutine +* dralum.f in dragr module. +* +*Copyright: +* Copyright (C) 2007 Ecole Polytechnique de Montreal +* This library is free software; you can redistribute it and/or +* modify it under the terms of the GNU Lesser General Public +* License as published by the Free Software Foundation; either +* version 2.1 of the License, or (at your option) any later version +* +*Author(s): A. Hebert +* +*Parameters: input +* IPRINT print parameter. +* IPEDIT pointer to the edition LCM object. +* MAXFP second dimension of array 'YIELD'. +* NBISO number of depleting nuclides before lumping. +* NBFISS number of fissile isotopes with fission yields. +* NBDPF number of fission products before lumping. +* NSUPS number of stable isotopes producing energy. +* NREAC number of depleting reactions including radioactive decay +* NFATH maximum number of parent isotopes leading to the production of +* an isotope in the depletion chain before lumping. +* NBCH number of depleting nuclides after lumping. +* HICH names of remaining isotopes after lumping. +* HISO 'ISOTOPESDEPL' record before lumping. +* MYLIST 'CHARGEWEIGHT' record before lumping. +* HREAC 'DEPLETE-IDEN' record before lumping. +* IDREAC 'DEPLETE-REAC' record before lumping. +* DENER 'DEPLETE-ENER' record before lumping. +* DDECA 'DEPLETE-DECA' record before lumping. +* IPREAC 'PRODUCE-REAC' record before lumping. +* PRATE 'PRODUCE-RATE' record before lumping. +* YIELD 'FISSIONYIELD' record before lumping. +* LISO =.true. when all isotopes are kept separately during +* the merging. +* NBFISS2 new number of fissile isotopes with fission yields. +* NBFPCH2 new maximum number of fission products after lumping. +* +*----------------------------------------------------------------------- +* + USE GANLIB +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPEDIT + INTEGER IPRINT,MAXFP,NBISO,NBFISS,NBDPF,NSUPS,NREAC,NFATH, + & MYLIST(NBISO),HICH(3,NBCH),HISO(3,NBISO),NBCH,HREAC(2,NREAC), + & IDREAC(NREAC,NBISO),IPREAC(NFATH,NBISO),NBFISS2,NBFPCH2 + REAL DENER(NREAC,NBISO),DDECA(NBISO),PRATE(NFATH,NBISO), + & YIELD(NBFISS,MAXFP) + LOGICAL LISO +*---- +* LOCAL VARIABLES +*---- + PARAMETER (NSTATE=40,MAXIT=20,NSYSO=6) + CHARACTER HNAME*8,TEXT4*4 + INTEGER ISTATE(NSTATE),MIX(NBCH),IMIX + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: JPREAC,JDREAC,IPOS,HHHH + REAL, ALLOCATABLE, DIMENSION(:) :: EYIEL,DDDD + REAL, ALLOCATABLE, DIMENSION(:,:) :: RRATE,EENER +*---- +* SCRATCH STORAGE ALLOCATION +*---- + PARAMETER (MAXFAT=25) + ALLOCATE(JPREAC(MAXFAT,NBCH),JDREAC(NREAC,NBCH),IPOS(NBCH,3), + & HHHH(3,NBCH)) + ALLOCATE(RRATE(MAXFAT,NBCH),EENER(NREAC,NBCH), + & EYIEL(NBFISS2*NBFPCH2),DDDD(NBCH)) +*---- +* FIND THE POSITION OF THE LUMPED ISOTOPES IN THE COMPLETE CHAIN +*---- + J0OLD=0 + J0=0 + DO ISO=1,NBCH + DO JSO=1,NBISO + J0=JSO + IF((HICH(1,ISO).EQ.HISO(1,JSO)).AND. + & (HICH(2,ISO).EQ.HISO(2,JSO))) GO TO 10 + ENDDO + WRITE(NSYSO,'(/35H EDILUM: LIST OF AVAILABLE ISOTOPES)') + DO JSO=1,NBISO + WRITE(NSYSO,'(2X,I5,5H --> ,2A4)') JSO,HISO(1,JSO),HISO(2,JSO) + ENDDO + WRITE(HNAME,'(2A4)') HICH(1,ISO),HICH(2,ISO) + CALL XABORT('EDILUM: UNABLE TO FIND '//HNAME) + 10 IF(.NOT.LISO) THEN + IF(J0.LE.J0OLD) THEN + WRITE(HNAME,'(2A4)') HISO(1,J0),HISO(2,J0) + CALL XABORT('EDILUM: ISOTOPIC DATA NOT SORTED:'//HNAME) + ENDIF + ENDIF + J0OLD=J0 + IPOS(ISO,1)=J0 + ENDDO + DO I=1,NBCH + IMIX=1 + IF(I.EQ.1) GOTO 15 + DO J=1,I-1 + IF(IPOS(I,1).EQ.IPOS(J,1)) IMIX=IMIX+1 + ENDDO + 15 MIX(I)=IMIX + ENDDO +*---- +* REMOVE THE LUMPED FATHERS THAT ARE NEUTRON INDUCED REACTIONS. +*---- + DO ISO=1,NBISO + DO IFATH=1,NFATH + IF(IPREAC(IFATH,ISO).EQ.0) GO TO 20 + IF(MOD(IPREAC(IFATH,ISO),100).EQ.1) GO TO 20 + JND=IPREAC(IFATH,ISO)/100 + DO J=1,NBCH + IF(IPOS(J,1).EQ.JND) GO TO 20 + ENDDO + IPREAC(IFATH,ISO)=0 + PRATE(IFATH,ISO)=0.0 + 20 CONTINUE + ENDDO + ENDDO +*---- +* LUMP IDREAC, DENER, IPREAC AND PRATE. +*---- + ITER=0 + 40 ITER=ITER+1 + IF(ITER.GT.MAXIT) CALL XABORT('EDILUM: TOO MANY ITERATIONS.') + NLUMP=0 + DO ISO=1,NBCH + IND=IPOS(ISO,1) + DO IFATH=1,NFATH + IF(IPREAC(IFATH,IND).EQ.0) GO TO 50 + IF(MOD(IPREAC(IFATH,IND),100).NE.1) GO TO 50 + JND=IPREAC(IFATH,IND)/100 + IF(MYLIST(JND).EQ.0) GO TO 50 + DO J=1,NBCH + IF(IPOS(J,1).EQ.JND) GO TO 50 + ENDDO + NLUMP=NLUMP+1 ! ISOTOPE JND IS LUMPED + DO IDA=1,NBISO + DO IFA=1,NFATH + IPGAR=IPREAC(IFA,IDA) + IF((IPGAR/100.EQ.JND).AND.(MOD(IPGAR,100).EQ.1)) THEN + IF(MYLIST(IDA).EQ.0) GO TO 50 + ENDIF + ENDDO + ENDDO + DO IDA=1,NBISO + DO IFA=1,NFATH + IPGAR=IPREAC(IFA,IDA) + IF((IPGAR/100.EQ.JND).AND.(MOD(IPGAR,100).EQ.1)) THEN + IF(IDA.EQ.JND) CALL XABORT('EDILUM: BUG.') + PRGAR=PRATE(IFA,IDA) + DO IM=IFA,NFATH-1 + IPREAC(IM,IDA)=IPREAC(IM+1,IDA) + PRATE(IM,IDA)=PRATE(IM+1,IDA) + ENDDO + IPREAC(NFATH,IDA)=0 + PRATE(NFATH,IDA)=0.0 + DO JFATH=1,NFATH + IF(IPREAC(JFATH,JND).EQ.0) GO TO 45 + IM=NFATH+1 + DO K=NFATH,1,-1 + IF(IPREAC(K,IDA).EQ.IPREAC(JFATH,JND)) THEN + PRATE(K,IDA)=PRATE(K,IDA)+PRGAR*PRATE(JFATH,JND) + GO TO 44 + ENDIF + IF(IPREAC(K,IDA).EQ.0) IM=K + ENDDO + IF(IM.GT.NFATH) CALL XABORT('EDILUM: NFATH OVERFLOW.') + IPREAC(IM,IDA)=IPREAC(JFATH,JND) + PRATE(IM,IDA)=PRGAR*PRATE(JFATH,JND) + 44 CONTINUE + ENDDO + 45 IF(MOD(IDREAC(2,JND),100).EQ.5) THEN + JFP=IDREAC(2,JND)/100 + IF(MOD(IDREAC(2,IDA),100).EQ.5) THEN + IFP=IDREAC(2,IDA)/100 + ELSE + NBDPF=NBDPF+1 + IF(NBDPF.GT.MAXFP) THEN + CALL XABORT('EDILUM: MAXFP OVERFLOW.') + ENDIF + IFP=NBDPF + YIELD(:NBFISS,IFP)=0.0 + ENDIF + DO IFI=1,NBFISS + YIELD(IFI,IFP)=YIELD(IFI,IFP)+YIELD(IFI,JFP)*PRGAR + ENDDO + IDREAC(2,IDA)=IFP*100+5 + ENDIF + ENDIF + ENDDO + ENDDO + DO JFATH=1,NFATH + IF(IPREAC(JFATH,JND).GT.0) THEN + KT=MOD(IPREAC(JFATH,JND),100) + KND=IPREAC(JFATH,JND)/100 + DENER(KT,KND)=DENER(KT,KND)+PRATE(JFATH,JND)*DENER(1,JND) + ENDIF + IPREAC(JFATH,JND)=0 + PRATE(JFATH,JND)=0.0 + ENDDO + YMAX=0.0 + IF(MOD(IDREAC(2,JND),100).EQ.5) THEN + JFP=IDREAC(2,JND)/100 + DO KSO=1,NBISO + IF(MOD(IDREAC(2,KSO),100).EQ.4) THEN + IFI=IDREAC(2,KSO)/100 + DENER(2,KSO)=DENER(2,KSO)+YIELD(IFI,JFP)*DENER(1,JND) + ENDIF + ENDDO + DO IFI=1,NBFISS + YMAX=MAX(YMAX,ABS(YIELD(IFI,JFP))) + YIELD(IFI,JFP)=0.0 + ENDDO + DENER(2,JND)=0.0 + IDREAC(2,JND)=0 + ENDIF + DENER(1,JND)=0.0 + IDREAC(1,JND)=0 + HALF=1.0E8*LOG(2.0)/DDECA(JND)/86400.0 + IF(DDECA(JND).EQ.0.0) THEN + WRITE(HNAME,'(2A4)') HISO(1,JND),HISO(2,JND) + IF(YMAX.GT.1.0E-2) THEN + WRITE(NSYSO,6020) TRIM(HNAME),HALF,YMAX*100.0 + CALL XABORT('EDILUM: ISOTOPE '//HNAME//' SHOULD NOT BE L' + & //'UMPED.(1)') + ENDIF + IF(IPRINT.GT.2) WRITE(NSYSO,6020) HNAME,HALF,YMAX*100.0 + ELSE IF((HALF.GT.30.0).AND.(HALF.LT.999999.99)) THEN + WRITE(HNAME,'(2A4)') HISO(1,JND),HISO(2,JND) + IF(YMAX.GT.1.0E-2) THEN + WRITE(NSYSO,6020) TRIM(HNAME),HALF,YMAX*100.0 + CALL XABORT('EDILUM: ISOTOPE '//HNAME//' SHOULD NOT BE L' + & //'UMPED.(2)') + ENDIF + IF(IPRINT.GT.2) WRITE(NSYSO,6020) HNAME,HALF,YMAX*100.0 + ELSE IF(HALF.GT.30.0) THEN + WRITE(HNAME,'(2A4)') HISO(1,JND),HISO(2,JND) + IF(YMAX.GT.1.0E-2) THEN + WRITE(NSYSO,6020) TRIM(HNAME),HALF,YMAX*100.0 + CALL XABORT('EDILUM: ISOTOPE '//HNAME//' SHOULD NOT BE L' + & //'UMPED.(3)') + ENDIF + IF(IPRINT.GT.2) WRITE(NSYSO,6020) HNAME,HALF,YMAX*100.0 + ENDIF + DDECA(JND)=0.0 + MYLIST(JND)=0 + 50 CONTINUE + ENDDO + ENDDO + IF(IPRINT.GT.2) WRITE(NSYSO,'('' ......... NLUMP='',I5)') NLUMP + IF(NLUMP.GT.0) GO TO 40 +*---- +* WRITE VECTORS 'PRODUCE-REAC' AND 'PRODUCE-RATE' TO THE LCM OBJECT +*---- + DO ISO=1,NBCH + DO IFATH=1,MAXFAT + JPREAC(IFATH,ISO)=0 + RRATE(IFATH,ISO)=0.0 + ENDDO + IND=IPOS(ISO,1) + NN=0 + DO IFATH=1,NFATH + IF(IPREAC(IFATH,IND).NE.0) THEN + DO J=1,IFATH-1 + IF(IPREAC(IFATH,IND).EQ.IPREAC(J,IND)) THEN + JND1=IPREAC(IFATH,IND)/100 + JND2=IPREAC(J,IND)/100 + WRITE(NSYSO,'(/27H EDILUM: DUPLICATE FATHERS:,2A4, + & 1X,2A4)') HISO(1,JND1),HISO(2,JND1),HISO(1,JND2), + & HISO(2,JND2) + WRITE(HNAME,'(2A4)') HISO(1,IND),HISO(2,IND) + CALL XABORT('EDILUM: DUPLICATE FATHERS FOR '//HNAME) + ENDIF + ENDDO + DO I=1,NBCH + JSO=I + IF((IPOS(I,1).EQ.IPREAC(IFATH,IND)/100).AND. + & (MIX(I).EQ.MIX(ISO)))GO TO 70 + ENDDO + JSO=-1 + 70 IF(JSO.EQ.-1) THEN + JND=IPREAC(IFATH,IND)/100 + IF(IPRINT.GT.2) THEN + WRITE(NSYSO,'(/24H EDILUM: UNKNOWN FATHER ,2A4,5H FOR , + & 2A4)') HISO(1,JND),HISO(2,JND),HISO(1,IND),HISO(2,IND) + ENDIF + ELSE + NN=NN+1 + IF(NN.GT.MAXFAT) THEN + WRITE(TEXT4,'(I4)') NN + CALL XABORT('EDILUM: MAXFAT OVERFLOW NN='//TEXT4) + ENDIF + JPREAC(NN,ISO)=100*JSO+MOD(IPREAC(IFATH,IND),100) + RRATE(NN,ISO)=PRATE(IFATH,IND) + ENDIF + ENDIF + ENDDO + ENDDO + CALL LCMPUT(IPEDIT,'PRODUCE-REAC',MAXFAT*NBCH,1,JPREAC) + CALL LCMPUT(IPEDIT,'PRODUCE-RATE',MAXFAT*NBCH,2,RRATE) +*---- +* WRITE THE ISOTOPE ASCII NAMES ON LCM OBJECT +*---- + CALL LCMPUT(IPEDIT,'DEPLETE-IDEN',2*NREAC,3,HREAC(1,1)) +*---- +* WRITE THE LUMPED FISSION YIELD MATRIX TO THE LCM OBJECT +*---- + IBFI=0 + IBFP=0 + DO ISO=1,NBCH + IND=IPOS(ISO,1) + IPOS(ISO,2)=0 + IPOS(ISO,3)=0 + IF(MOD(IDREAC(2,IND),100).EQ.4) THEN + IBFI=IBFI+1 + IF(IBFI.GT.NBFISS2) CALL XABORT('EDILUM: NBFISS OVERFLOW.') + IPOS(ISO,2)=IBFI + ELSE IF(MOD(IDREAC(2,IND),100).EQ.5) THEN + IBFP=IBFP+1 + IF(IBFP.GT.NBFPCH2) CALL XABORT('EDILUM: NBFPCH2 OVERFLOW.') + IPOS(ISO,3)=IBFP + ENDIF + ENDDO + DO ISO=1,NBCH + IFI=IPOS(ISO,2) + IF(IFI.GT.0) THEN + DO JSO=1,NBCH + IFP=IPOS(JSO,3) + IF(IFP.GT.0) THEN + IND=IPOS(ISO,1) + JND=IPOS(JSO,1) + IF(MIX(ISO).NE.MIX(JSO)) THEN + EYIEL((IFP-1)*IBFI+IFI)=0.0 + ELSE + EYIEL((IFP-1)*IBFI+IFI)=YIELD(IDREAC(2,IND)/100, + & IDREAC(2,JND)/100) + ENDIF + ENDIF + ENDDO + ENDIF + ENDDO + IF(IBFI*IBFP.GT.0) CALL LCMPUT(IPEDIT,'FISSIONYIELD',IBFI*IBFP,2, + & EYIEL) +*---- +* WRITE VECTORS 'DEPLETE-REAC' AND 'DEPLETE-ENER' TO THE LCM OBJECT +*---- + DO ISO=1,NBCH + IND=IPOS(ISO,1) + DO I=1,NREAC + IF(IDREAC(I,IND)/100.GT.0) THEN + KREAC=MOD(IDREAC(I,IND),100) + IF((KREAC.LE.0).OR.(KREAC.GT.5)) THEN + CALL XABORT('EDILUM: INVALID REACTION.') + ENDIF + ENDIF + IF((I.EQ.2).AND.(MOD(IDREAC(I,IND),100).EQ.4)) THEN + IF(IPOS(ISO,2).EQ.0) CALL XABORT('EDILUM: FAILURE 1.') + JDREAC(I,ISO)=IPOS(ISO,2)*100+4 + IF(IPRINT.GT.2) THEN + WRITE(NSYSO,6010) (HISO(I0,IPOS(ISO,1)),I0=1,3) + ENDIF + ELSE IF((I.EQ.2).AND.(MOD(IDREAC(I,IND),100).EQ.5)) THEN + IF(IPOS(ISO,3).EQ.0) CALL XABORT('EDILUM: FAILURE 2.') + JDREAC(I,ISO)=IPOS(ISO,3)*100+5 + ELSE + JDREAC(I,ISO)=IDREAC(I,IND) + ENDIF + EENER(I,ISO)=DENER(I,IND) + ENDDO + ENDDO + CALL LCMPUT(IPEDIT,'DEPLETE-REAC',NREAC*NBCH,1,JDREAC) + CALL LCMPUT(IPEDIT,'DEPLETE-ENER',NREAC*NBCH,2,EENER) +*---- +* WRITE VECTORS 'CHARGEWEIGHT', 'DEPLETE-DECA', 'ISOTOPESDEPL' AND +* 'STATE-VECTOR' TO THE LCM OBJECT +*---- + NBHEAV=0 + NSUPS2=0 + DO ISO=1,NBCH + IF(IPOS(ISO,1).GT.NBISO-NSUPS) NSUPS2=NSUPS2+1 + HHHH(1,ISO)=HISO(1,IPOS(ISO,1)) + HHHH(2,ISO)=HISO(2,IPOS(ISO,1)) + HHHH(3,ISO)=HISO(3,IPOS(ISO,1)) + DDDD(ISO)=DDECA(IPOS(ISO,1)) + IPOS(ISO,1)=MYLIST(IPOS(ISO,1)) + ENDDO + CALL LCMPUT(IPEDIT,'ISOTOPESDEPL',3*NBCH,3,HHHH) + CALL LCMPUT(IPEDIT,'CHARGEWEIGHT',NBCH,1,IPOS(1,1)) + CALL LCMPUT(IPEDIT,'DEPLETE-DECA',NBCH,2,DDDD) + ISTATE(:NSTATE)=0 + ISTATE(1)=NBCH + ISTATE(2)=IBFI + ISTATE(3)=IBFP + ISTATE(4)=NBHEAV + ISTATE(5)=NBCH-NBHEAV + ISTATE(7)=NSUPS2 + ISTATE(8)=NREAC + ISTATE(9)=MAXFAT + CALL LCMPUT(IPEDIT,'STATE-VECTOR',NSTATE,1,ISTATE) + IF(IPRINT.GT.2) WRITE(NSYSO,6000) (ISTATE(ISTA),ISTA=1,9) +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(DDDD,EYIEL,EENER,RRATE) + DEALLOCATE(HHHH,IPOS,JDREAC,JPREAC) + RETURN +*---- +* FORMAT +*---- + 6000 FORMAT(/' STATE-VECTOR FOR LUMPED DEPLETION CHAIN'/' -------'/ + > ' NDEPL ',I6,' (NUMBER OF DEPLETING ISOTOPES)'/ + > ' NDFI ',I6,' (NUMBER OF DIRECT FISSILE ISOTOPES)'/ + > ' NDFP ',I6,' (NUMBER OF DIRECT FISSION PRODUCT)'/ + > ' NHEAVY ',I6,' (NUMBER OF HEAVY ISOTOPES)'/ + > ' NLIGHT ',I6,' (NUMBER OF FISSION PRODUCTS)'/ + > ' NOTHER ',I6,' (NUMBER OF OTHER ISOTOPES)'/ + > ' NSTABL ',I6,' (NUMBER OF STABLE ISOTOPES PRODUCING ENERGY)'/ + > ' NREAC ',I6,' (MAXIMUM NUMBER OF DEPLETION REACTIONS)'/ + > ' NPAR ',I6,' (MAXIMUM NUMBER OF PARENT REACTIONS)'/) + 6010 FORMAT(45H EDILUM: FISSILE ISOTOPE WITH FISSION YIELD: ,3A4) + 6020 FORMAT(18H EDILUM: ISOTOPE ',A,30H' IS LUMPED AND HAS A HALF-LIF, + > 4HE OF,1P,E12.4,25H DAYS. MAX FISSION YIELD=,1P,E8.1,2H%.) + END |