Initial commit from Polytechnique Montreal

1 files changed, 312 insertions, 0 deletions

diff --git a/Dragon/src/LIBWET.f b/Dragon/src/LIBWET.f
new file mode 100644
index 0000000..f30cd99
--- /dev/null
+++ b/Dragon/src/LIBWET.f
@@ -0,0 +1,312 @@
+*DECK LIBWET
+      SUBROUTINE LIBWET(MAXR,NEL,NBESP,NSTATE,NMDEPL,ITNAM,ISTATE,
+     >                  MATNO,KPAX,BPAX)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Identify the depleting isotopes by type and reorder them (recompute
+* the KPAX and BPAX matrices).
+*
+*Copyright:
+* Copyright (C) 2002 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): G. Marleau
+*
+*Parameters: input
+* MAXR    number of reaction types.
+* NEL     number of isotopes on library.
+* NBESP   number of energy-dependent fission yield matrices.
+* NSTATE  number of parameters in the state vector.
+* NMDEPL  names of reactions:
+*           NMDEPL(1)='DECAY'; NMDEPL(2)='NFTOT';
+*           NMDEPL(3)='NG'   ; NMDEPL(4)='N2N';
+*           etc.
+* ITNAM   reactive isotope names in chain.
+*
+*Parameters: output
+* ISTATE  state vector containing the library parameters.
+* MATNO   reaction material indices.
+*
+*Parameters: input/output
+* KPAX    complete reaction type matrix.
+* BPAX    complete branching ratio matrix.
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER      MAXR,NEL,NBESP,NSTATE,ITNAM(3,NEL),ISTATE(NSTATE),
+     >             MATNO(NEL),KPAX(NEL+MAXR,NEL)
+      CHARACTER    NMDEPL(MAXR)*8
+      REAL         BPAX(NBESP,NEL+MAXR,NEL)
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER      NDFP,JSO,ISO,ITR,NUNTIL,NHEAVY,NDFI,IUNTIL,NEW,
+     >             NLIGHT,NOTHER,NIFP,KPAX0,KPAX1,KPAX2,IEL,JEL,
+     >             KREAC,II,NPAR,ICOUNT,NREAC,NSTABL
+      REAL         SUMFI
+      CHARACTER    TEXT12*12
+      LOGICAL      LOGPF
+*----
+*  INTERNAL PARAMETERS
+*----
+      INTEGER      KDECAY,KFISSP,KCAPTU,KFISSI,KHEAT
+      PARAMETER   (KDECAY=1,KFISSP=2,KCAPTU=3,KFISSI=6,KHEAT=9)
+*----
+*  COMPUTE NUMBER OF DIRECT FISSION PRODUCT
+*  NUMBER DIRECT FISSION PRODUCT
+*
+*  SET MATNO TO LIGHT (-KFISSP) FOR FISSION PRODUCT
+*  SET MATNO TO HEAVY (-KFISSI) FOR FISSILE ISOTOPES
+*  SET MATNO TO STABLE (-KHEAT) FOR STABLE ISOTOPES PRODUCING ENERGY
+*  SET MATNO TO 0 TO UNUSED ISOTOPES
+*----
+      NDFI=0
+      NDFP=0
+      DO 100 JSO=NEL,1,-1
+        MATNO(JSO)=0
+        KPAX0=KPAX(NEL+KDECAY,JSO)
+        KPAX1=KPAX(NEL+KFISSP,JSO)
+        KPAX2=KPAX(NEL+KCAPTU,JSO)
+        IF((KPAX0.EQ.-9999).OR.(KPAX1.EQ.-9999).OR.(KPAX2.EQ.-9999))
+     >  THEN
+*         JSO IS A STABLE ISOTOPE
+          MATNO(JSO)=-KHEAT
+          DO 222 ITR=1,MAXR
+            IF(BPAX(1,NEL+ITR,JSO).GT.0.0) THEN
+              KPAX(NEL+ITR,JSO)=2
+              KPAX(JSO,:NEL)=0
+            ENDIF
+ 222      CONTINUE
+          GO TO 100
+        ELSE IF(KPAX1.LT.0) THEN
+*         JSO IS A FISSION PRODUCT
+          MATNO(JSO)=-KFISSP
+          NIFP=0
+          DO 101 ISO=NEL,1,-1
+            IF((KPAX(JSO,ISO).EQ.KFISSP).AND.
+     >         (BPAX(1,JSO,ISO).GT.0.0)) THEN
+              NIFP=NIFP+1
+            ENDIF
+ 101      CONTINUE
+          IF(NIFP.GT.0) THEN
+            NDFP=NDFP+1
+            KPAX(NEL+KFISSP,JSO)=5+100*NDFP
+          ELSE
+            KPAX(NEL+KFISSP,JSO)=5
+          ENDIF
+        ELSE IF(KPAX1.GT.0) THEN
+*         JSO IS A FISSILE ISOTOPE
+          MATNO(JSO)=-KFISSI
+          SUMFI=0.0
+          DO 221 ISO=1,NEL
+            IF(KPAX(ISO,JSO).EQ.KFISSP) SUMFI=SUMFI+BPAX(1,ISO,JSO)
+ 221      CONTINUE
+          IF(SUMFI.GT.0.0) THEN
+            NDFI=NDFI+1
+            KPAX(NEL+KFISSP,JSO)=4+100*NDFI
+          ELSE
+            KPAX(NEL+KFISSP,JSO)=3
+          ENDIF
+        ENDIF
+ 100  CONTINUE
+*----
+*  CHECK IF THE DEPLETION CHAIN IS COHERENT
+*----
+      DO 200 IEL=1,NEL
+        DO 210 JEL=1,NEL
+          KREAC=KPAX(JEL,IEL)
+          IF(KREAC.EQ.KFISSP) THEN
+            IF(MOD(KPAX(NEL+KFISSP,JEL),100).NE.5) THEN
+              WRITE(TEXT12,'(3A4)') (ITNAM(II,JEL),II=1,3)
+              CALL XABORT('LIBWET: SON '//TEXT12//' IS NOT A FISSION '//
+     >                    'PRODUCT')
+            ENDIF
+            IF((KPAX(NEL+KFISSP,IEL).NE.3).AND.
+     >         (MOD(KPAX(NEL+KFISSP,IEL),100).NE.4)) THEN
+              WRITE(TEXT12,'(3A4)') (ITNAM(II,IEL),II=1,3)
+              CALL XABORT('LIBWET: PARENT '//TEXT12//' IS NOT A FISSI'//
+     >                    'LE ISOTOPE')
+            ENDIF
+          ELSE IF(KREAC.NE.0) THEN
+            IF(KPAX(NEL+KREAC,IEL).EQ.0) THEN
+              WRITE(TEXT12,'(3A4)') (ITNAM(II,IEL),II=1,3)
+              CALL XABORT('LIBWET: PARENT '//TEXT12//' DOES NOT DEPL'//
+     >                    'ETE VIA REACTION '//NMDEPL(KREAC))
+            ENDIF
+          ENDIF
+ 210    CONTINUE
+ 200  CONTINUE
+*----
+*  SET MATNO TO OTHER (-KDECAY) FOR ISOTOPES WITH DAUGHTER OR FATHER
+*  AND FOR ISOTOPES WITH DECAY
+*----
+      DO 112 ISO=1,NEL
+        IF(MATNO(ISO).EQ.0) THEN
+          IF(KPAX(NEL+KDECAY,ISO).GT.0 .OR.
+     >       KPAX(NEL+KCAPTU,ISO).GT.0) THEN
+            MATNO(ISO)=-KDECAY
+            GO TO 115
+          ENDIF
+          DO 113 JSO=1,NEL
+            IF(KPAX(JSO,ISO).GT.0) THEN
+              MATNO(ISO)=-KDECAY
+              GO TO 115
+            ELSE IF(KPAX(ISO,JSO).GT.0) THEN
+              MATNO(ISO)=-KDECAY
+              GO TO 115
+            ENDIF
+ 113      CONTINUE
+        ENDIF
+ 115    CONTINUE
+ 112  CONTINUE
+*----
+*  SET MATNO TO STABLE (-KHEAT) FOR OTHER ISOTOPES PRODUCING ENERGY
+*----
+      DO 212 ISO=1,NEL
+        IF(MATNO(ISO).EQ.0) THEN
+          DO 213 ITR=2,MAXR
+            IF(BPAX(1,NEL+ITR,ISO).NE.0.0) THEN
+              MATNO(ISO)=-KHEAT
+              GO TO 215
+            ENDIF
+ 213      CONTINUE
+        ENDIF
+ 215    CONTINUE
+ 212  CONTINUE
+*----
+*  COMPUTE THE MAXIMUM NUMBER OF PARENTS FOR AN ISOTOPE
+*----
+      NPAR=0
+      DO 116 ISO=1,NEL
+        ICOUNT=0
+        DO 114 JSO=1,NEL
+          LOGPF=((MATNO(ISO).EQ.-KFISSP).AND.(MATNO(JSO).EQ.-KFISSI))
+          IF((BPAX(1,ISO,JSO).NE.0.0).AND.(.NOT.LOGPF)) ICOUNT=ICOUNT+1
+ 114    CONTINUE
+        NPAR=MAX(NPAR,ICOUNT)
+ 116  CONTINUE
+*----
+*  COMPUTE THE MAXIMUM NUMBER OF DEPLETION REACTIONS
+*----
+      NREAC=4
+      DO 118 ISO=1,NEL
+        DO 117 ITR=1,MAXR
+          IF(KPAX(NEL+ITR,ISO).NE.0) NREAC=MAX(NREAC,ITR)
+ 117    CONTINUE
+ 118  CONTINUE
+*----
+*  SET MATNO TO HEAVY (-KFISSI) FOR ISOTOPES RESULTING FROM DECAY
+*  OR CAPTURE OF HEAVY ISOTOPE UNTIL ALL HEAVY ISOTOPES IDENTIFIED
+*----
+      NUNTIL=NEL
+      NHEAVY=0
+      DO 120 IUNTIL=1,NUNTIL
+        NEW=0
+        DO 121 ISO=NEL,1,-1
+          IF(MATNO(ISO).EQ.-KFISSI) THEN
+            NHEAVY=NHEAVY+1
+            DO 122 JSO=NEL,1,-1
+              IF(MATNO(JSO).NE.-KFISSI) THEN
+                IF((KPAX(JSO,ISO).NE.0).AND.(KPAX(JSO,ISO).NE.KFISSP))
+     >          THEN
+                  NEW=NEW+1
+                  MATNO(JSO)=-KFISSI
+                ENDIF
+              ENDIF
+ 122        CONTINUE
+          ENDIF
+ 121    CONTINUE
+        IF(NEW.EQ.0) GO TO 125
+        NHEAVY=0
+ 120  CONTINUE
+ 125  CONTINUE
+*----
+*  SET MATNO TO LIGHT (-KFISSP) FOR ISOTOPES RESULTING FROM DECAY
+*  OR CAPTURE OF LIGHT ISOTOPE UNTIL ALL LIGHT ISOTOPES IDENTIFIED
+*----
+      NUNTIL=NUNTIL-NHEAVY
+      NLIGHT=0
+      DO 130 IUNTIL=1,NUNTIL
+        NEW=0
+        DO 131 ISO=NEL,1,-1
+          IF(MATNO(ISO).EQ.-KFISSP) THEN
+            NLIGHT=NLIGHT+1
+            DO 132 JSO=NEL,1,-1
+              IF(MATNO(JSO).NE.-KFISSI.AND.
+     >           MATNO(JSO).NE.-KFISSP) THEN
+                IF(KPAX(JSO,ISO).NE.0) THEN
+                  NEW=NEW+1
+                  MATNO(JSO)=-KFISSP
+                ENDIF
+              ENDIF
+ 132        CONTINUE
+          ENDIF
+ 131    CONTINUE
+        IF(NEW.EQ.0) GO TO 135
+        NLIGHT=0
+ 130  CONTINUE
+ 135  CONTINUE
+*----
+*  SET MATNO TO OTHER (-KDECAY) FOR ISOTOPES RESULTING FROM DECAY
+*  OR CAPTURE OF OTHER ISOTOPE UNTIL ALL OTHER ISOTOPES IDENTIFIED
+*----
+      NUNTIL=NUNTIL-NLIGHT
+      NOTHER=0
+      DO 140 IUNTIL=1,NUNTIL
+        NEW=0
+        DO 141 ISO=NEL,1,-1
+          IF(MATNO(ISO).EQ.-KDECAY) THEN
+            NOTHER=NOTHER+1
+            DO 142 JSO=NEL,1,-1
+              IF(MATNO(JSO).NE.-KFISSI.AND.
+     >           MATNO(JSO).NE.-KFISSP.AND.
+     >           MATNO(JSO).NE.-KDECAY) THEN
+                IF(KPAX(JSO,ISO).NE.0) THEN
+                  NEW=NEW+1
+                  MATNO(JSO)=-KDECAY
+                ENDIF
+              ENDIF
+ 142        CONTINUE
+          ENDIF
+ 141    CONTINUE
+        IF(NEW.EQ.0) GO TO 145
+        NOTHER=0
+ 140  CONTINUE
+ 145  CONTINUE
+*----
+*  COUNT THE NUMBER OF STABLE ISOTOPES AND SET TO ZERO THEIR RADIOACTIVE
+*  DECAY CONSTANT.
+*----
+      NSTABL=0
+      DO 150 ISO=1,NEL
+        IF(MATNO(ISO).EQ.-KHEAT) THEN
+           NSTABL=NSTABL+1
+           BPAX(:,NEL+KDECAY,ISO)=0.0
+        ENDIF
+ 150  CONTINUE
+*
+      ISTATE(:NSTATE)=0
+      ISTATE(1)=NHEAVY+NLIGHT+NOTHER+NSTABL
+      ISTATE(2)=NDFI
+      ISTATE(3)=NDFP
+      ISTATE(4)=NHEAVY
+      ISTATE(5)=NLIGHT
+      ISTATE(6)=NOTHER
+      ISTATE(7)=NSTABL
+      ISTATE(8)=NREAC
+      ISTATE(9)=NPAR
+      ISTATE(10)=NBESP
+*----
+*  RETURN
+*----
+      RETURN
+      END