Initial commit from Polytechnique Montreal

1 files changed, 313 insertions, 0 deletions

diff --git a/Dragon/src/LIBDEP.F b/Dragon/src/LIBDEP.F
new file mode 100644
index 0000000..e1b539a
--- /dev/null
+++ b/Dragon/src/LIBDEP.F
@@ -0,0 +1,313 @@
+*DECK LIBDEP
+      SUBROUTINE LIBDEP(IPLIB,IMPX,NDEPL)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Read the information related to the depletion calculation.
+*
+*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): A. Hebert and G. Marleau
+*
+*Parameters: input
+* IPLIB   pointer to the internal microscopic cross section library
+*         (L_LIBRARY signature).
+* IMPX    print flag.
+*
+*Parameters: output
+* NDEPL   number of depleting isotopes.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+#if defined(HDF5_LIB)
+      USE hdf5_wrap
+#endif /* defined(HDF5_LIB) */
+      IMPLICIT NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPLIB
+      INTEGER IMPX,NDEPL
+*----
+*  LOCAL PARAMETERS
+*----
+      TYPE(C_PTR)      IPDRL
+      INTEGER          IOUT,NSTATE,MAXR,INDIC,NEL,IEVOT,NITMA,NDFI,
+     >                 NDFP,NHEAVY,NLIGHT,NOTHER,NSTABL,NREAC,NPAR,
+     >                 ITEXT4,I,J,ISTA,ILONG,ITYLCM,NBESP
+      REAL             FLOTT
+      PARAMETER       (IOUT=6,NSTATE=40,MAXR=12)
+      DOUBLE PRECISION DBLINP
+      CHARACTER        NMDEPL(MAXR)*8,TEXT4*4,HSMG*131,CFILNA*64,
+     >                 HHLIB*8,TEXT12*12,NAMLCM*12,NAMMY*12
+      LOGICAL          EMPTY,LCM,LEXIST
+      INTEGER          ISTATE(NSTATE)
+#if defined(HDF5_LIB)
+      CHARACTER        CFILNA1*64,CFILNA2*64
+#endif /* defined(HDF5_LIB) */
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: INAM,IZAE,HREAC,IDR,KPAR,
+     > ITNAM,ITZEA,MATNO,KPAX
+      REAL, ALLOCATABLE, DIMENSION(:) :: RER,RRD,BPAR,YIELD,BPAX,ENER
+#if defined(HDF5_LIB)
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: DIMS_AP
+#endif /* defined(HDF5_LIB) */
+*----
+*  DATA STATEMENTS
+*----
+      SAVE      NMDEPL
+      DATA      NMDEPL/'DECAY   ','NFTOT   ','NG      ','N2N     ',
+     >                 'N3N     ','N4N     ','NA      ','NP      ',
+     >                 'N2A     ','NNP     ','ND      ','NT      '/
+*----
+*  READ INFORMATION AVAILABLE ON INPUT
+*----
+      CALL REDGET(INDIC,NEL,FLOTT,TEXT4,DBLINP)
+      IEVOT=-99
+      NBESP=1
+      IF(INDIC.EQ.1) THEN
+         IEVOT=0
+      ELSE IF((INDIC.EQ.3).AND.(TEXT4.EQ.'LIB:')) THEN
+        CALL REDGET(INDIC,NITMA,FLOTT,HHLIB,DBLINP)
+        IF(INDIC.NE.3) THEN
+          CALL XABORT('LIBDEP: CHARACTER LIBRARY NAME REQUIRED.')
+        ELSE IF((HHLIB.NE.'DRAGON  ') .AND. (HHLIB.NE.'WIMSAECL') .AND.
+     >          (HHLIB.NE.'WIMSD4  ') .AND. (HHLIB.NE.'WIMSE   ') .AND.
+     >          (HHLIB.NE.'APLIB2  ') .AND. (HHLIB.NE.'APLIB3  ') .AND.
+     >          (HHLIB.NE.'NDAS    ') .AND. (HHLIB.NE.'APXSM   ') ) THEN
+          WRITE(HSMG,'(30HLIBDEP: INVALID EVOL LIB TYPE ,A8)') HHLIB
+          CALL XABORT(HSMG)
+        ENDIF
+        CALL REDGET(INDIC,NITMA,FLOTT,TEXT4,DBLINP)
+        IF((INDIC.NE.3).OR.(TEXT4.NE.'FIL:'))
+     >    CALL XABORT('LIBDEP: FIL: EXPECTED.')
+        CFILNA=' '
+        CALL REDGET(INDIC,NITMA,FLOTT,CFILNA,DBLINP)
+        IF(INDIC.NE.3) CALL XABORT('LIBDEP: CHARACTER DATA EXPECTED.')
+        IF(HHLIB.EQ.'DRAGON') THEN
+          TEXT12=CFILNA(:12)
+          CALL LCMINF(IPLIB,NAMLCM,NAMMY,EMPTY,ILONG,LCM)
+          IF(TEXT12.EQ.NAMLCM) THEN
+            IPDRL=IPLIB
+          ELSE
+            INQUIRE(FILE=TRIM(TEXT12),EXIST=LEXIST)
+            IF(.NOT.LEXIST) THEN
+              WRITE(HSMG,'(17HLIBDEP: XSM FILE ,A,14H DOESNT EXIST.)')
+     >        TRIM(TEXT12)
+              CALL XABORT(HSMG)
+            ENDIF
+            CALL LCMOP(IPDRL,TEXT12,2,2,0)
+          ENDIF
+          CALL LCMLEN(IPDRL,'DEPL-CHAIN',ILONG,ITYLCM)
+          IF(ILONG.EQ.0) THEN
+             CALL XABORT('LIBDEP: NO BURNUP DATA ON DRAGLIB NAMED '//
+     >       TEXT12//'.')
+          ENDIF
+          CALL LCMSIX(IPDRL,'DEPL-CHAIN',1)
+          CALL LCMGET(IPDRL,'STATE-VECTOR',ISTATE)
+          NDEPL=ISTATE(1)
+          NDFI=ISTATE(2)
+          NDFP=ISTATE(3)
+          NHEAVY=ISTATE(4)
+          NLIGHT=ISTATE(5)
+          NOTHER=ISTATE(6)
+          NSTABL=ISTATE(7)
+          NREAC=ISTATE(8)
+          NPAR=ISTATE(9)
+          NBESP=MAX(1,ISTATE(10))
+          ALLOCATE(INAM(3*NDEPL),IZAE(NDEPL),HREAC(2*NREAC),
+     1    IDR(NREAC*NDEPL),RER(NREAC*NDEPL),RRD(NDEPL),KPAR(NPAR*NDEPL),
+     2    BPAR(NPAR*NDEPL))
+          IF(NDFP.GT.0) ALLOCATE(YIELD(NBESP*NDFI*NDFP))
+          CALL LCMGET(IPDRL,'STATE-VECTOR',ISTATE)
+          CALL LCMGET(IPDRL,'ISOTOPESDEPL',INAM)
+          CALL LCMGET(IPDRL,'CHARGEWEIGHT',IZAE)
+          CALL LCMGET(IPDRL,'DEPLETE-IDEN',HREAC)
+          CALL LCMGET(IPDRL,'DEPLETE-REAC',IDR)
+          CALL LCMGET(IPDRL,'DEPLETE-ENER',RER)
+          CALL LCMGET(IPDRL,'DEPLETE-DECA',RRD)
+          CALL LCMGET(IPDRL,'PRODUCE-REAC',KPAR)
+          CALL LCMGET(IPDRL,'PRODUCE-RATE',BPAR)
+          IF(NDFI*NDFP.GT.0) CALL LCMGET(IPDRL,'FISSIONYIELD',YIELD)
+          CALL LCMSIX(IPDRL,' ',2)
+          IF(TEXT12.NE.NAMLCM) CALL LCMCL(IPDRL,1)
+          GO TO 20
+        ELSE IF(HHLIB.EQ.'WIMSAECL') THEN
+          CALL LIBEWI(CFILNA,NEL)
+          IEVOT=2
+        ELSE IF(HHLIB.EQ.'WIMSD4') THEN
+          CALL LIBENI(CFILNA,NEL)
+          IEVOT=3
+        ELSE IF(HHLIB.EQ.'APLIB2') THEN
+          CALL LIBEAI(CFILNA,NEL)
+          IEVOT=4
+        ELSE IF(HHLIB.EQ.'NDAS') THEN
+          CALL LIBND5(CFILNA,NEL)
+          IEVOT=5
+        ELSE IF(HHLIB.EQ.'APXSM') THEN
+          CALL LIBXS1(CFILNA,NEL)
+          IEVOT=6
+        ELSE IF(HHLIB.EQ.'WIMSE') THEN
+          CALL LIBENI(CFILNA,NEL)
+          IEVOT=7
+        ELSE IF(HHLIB.EQ.'APLIB3') THEN
+#if defined(HDF5_LIB)
+          I = INDEX(CFILNA, ":")
+          IF(I.EQ.0) THEN
+            CFILNA1=CFILNA
+            CFILNA2=" "
+          ELSE
+            CFILNA1=CFILNA(:I-1)
+            CFILNA2=CFILNA(I+1:)
+          ENDIF
+          CALL hdf5_open_file(CFILNA1, IPDRL, .TRUE.)
+          CALL hdf5_read_data(IPDRL, "Head/nbIs", NEL)
+          CALL hdf5_close_file(IPDRL)
+          IF(CFILNA2.NE.' ') THEN
+            CALL hdf5_open_file(CFILNA2, IPDRL, .TRUE.)
+            CALL hdf5_get_shape(IPDRL,"/Yields/YieldEnMshInMeV",DIMS_AP)
+            CALL hdf5_close_file(IPDRL)
+            NBESP=DIMS_AP(1)-1
+            DEALLOCATE(DIMS_AP)
+          ENDIF
+          IEVOT=8
+#else
+          CALL XABORT('LIBDEP: THE HDF5 API IS NOT AVAILABLE(1).')
+#endif /* defined(HDF5_LIB) */
+        ENDIF
+      ELSE
+        CALL XABORT('LIBDEP: INVALID KEY WORD.')
+      ENDIF
+      IF(IEVOT.EQ.0.OR.IEVOT.GT.1) THEN
+*----
+*  ALLOCATE/INITIALIZE WORK VECTORS FOR WIMS-AECL, WIMSD4
+*  AND INPUT FILE
+*----
+        ALLOCATE(ENER(NBESP+1),ITNAM(3*NEL),ITZEA(NEL),MATNO(NEL),
+     1  KPAX((NEL+MAXR)*NEL),BPAX((NEL+MAXR)*NEL*NBESP))
+        TEXT4='    '
+        READ(TEXT4,'(A4)') ITEXT4
+        ITNAM(:3*NEL)=ITEXT4
+        ITZEA(:NEL)=0
+        MATNO(:NEL)=0
+        KPAX(:(NEL+MAXR)*NEL)=0
+        BPAX(:(NEL+MAXR)*NEL*NBESP)=0.0
+        IF(IEVOT.EQ.0) THEN
+          CALL LIBEIR(MAXR,NEL,NMDEPL,ITNAM,ITZEA,KPAX,BPAX)
+        ELSE IF(IEVOT.EQ.2) THEN
+          CALL LIBEWR(CFILNA,MAXR,NEL,ITNAM,KPAX,BPAX)
+        ELSE IF(IEVOT.EQ.3) THEN
+          CALL LIBENR(CFILNA,4,MAXR,NEL,ITNAM,KPAX,BPAX)
+        ELSE IF(IEVOT.EQ.4) THEN
+          CALL LIBEAR(CFILNA,MAXR,NEL,NMDEPL,ITNAM,ITZEA,KPAX,BPAX)
+        ELSE IF(IEVOT.EQ.5) THEN
+          CALL LIBND6(CFILNA,MAXR,NEL,ITNAM,KPAX,BPAX)
+        ELSE IF(IEVOT.EQ.6) THEN
+          CALL LIBXS2(CFILNA,MAXR,NEL,NMDEPL,ITNAM,ITZEA,KPAX,BPAX)
+        ELSE IF(IEVOT.EQ.7) THEN
+          CALL LIBENR(CFILNA,5,MAXR,NEL,ITNAM,KPAX,BPAX)
+        ELSE IF(IEVOT.EQ.8) THEN
+#if defined(HDF5_LIB)
+          CALL LIBE3R(CFILNA1,CFILNA2,MAXR,NEL,NBESP,IMPX,ITNAM,ITZEA,
+     1    KPAX,BPAX,ENER)
+#else
+          CALL XABORT('LIBDEP: THE HDF5 API IS NOT AVAILABLE(2).')
+#endif /* defined(HDF5_LIB) */
+        ENDIF
+        CALL LIBWET(MAXR,NEL,NBESP,NSTATE,NMDEPL,ITNAM,ISTATE,MATNO,
+     >  KPAX,BPAX)
+        NDEPL=ISTATE(1)
+        NDFI=ISTATE(2)
+        NDFP=ISTATE(3)
+        NHEAVY=ISTATE(4)
+        NLIGHT=ISTATE(5)
+        NOTHER=ISTATE(6)
+        NSTABL=ISTATE(7)
+        NREAC=ISTATE(8)
+        NPAR=ISTATE(9)
+        NBESP=ISTATE(10)
+      ENDIF
+*----
+*  ALLOCATE DECAY CHAIN
+*----
+      NDEPL=MAX(NDEPL,1)
+      NDFI=MAX(NDFI,1)
+      NDFP=MAX(NDFP,1)
+      ALLOCATE(INAM(3*NDEPL),IZAE(NDEPL),IDR(NREAC*NDEPL),
+     1 RER(NREAC*NDEPL),RRD(NDEPL),KPAR(NPAR*NDEPL),BPAR(NPAR*NDEPL),
+     2 YIELD(NDFI*NDFP*NBESP))
+*----
+*  SET DECAY CHAIN
+*----
+      CALL LIBWED(MAXR,NEL,NBESP,NDEPL,NDFI,NDFP,NHEAVY,NLIGHT,NOTHER,
+     >            NREAC,NPAR,ITNAM,ITZEA,MATNO,KPAX,BPAX,INAM,IZAE,IDR,
+     >            RER,RRD,KPAR,BPAR,YIELD)
+*----
+*  RELEASE WORK VECTORS FOR WIMS-AECL, WIMS-NEA, DRAGLIB
+*  AND INPUT FILE
+*----
+      DEALLOCATE(BPAX,KPAX,MATNO,ITZEA,ITNAM)
+*----
+*  SELECT USED DEPLETION REACTION NAMES
+*----
+      ALLOCATE(HREAC(2*NREAC))
+      DO 10 I=1,NREAC
+        READ(NMDEPL(I),'(2A4)') (HREAC(2*(I-1)+J),J=1,2)
+   10 CONTINUE
+*----
+*  PRINT DECAY CHAIN IF REQUIRED
+*----
+   20 CALL LIBEPR(IMPX,NBESP,NDEPL,NSTABL,NDFI,NDFP,NREAC,NPAR,INAM,
+     >            HREAC,IDR,RER,RRD,KPAR,BPAR,YIELD,IZAE)
+*----
+*  SAVE CHAIN
+*----
+      CALL LCMSIX(IPLIB,'DEPL-CHAIN',1)
+      CALL LCMPUT(IPLIB,'STATE-VECTOR',NSTATE,1,ISTATE)
+      CALL LCMPUT(IPLIB,'ISOTOPESDEPL',3*NDEPL,3,INAM)
+      CALL LCMPUT(IPLIB,'CHARGEWEIGHT',NDEPL,1,IZAE)
+      CALL LCMPUT(IPLIB,'DEPLETE-IDEN',2*NREAC,3,HREAC)
+      CALL LCMPUT(IPLIB,'DEPLETE-REAC',NREAC*NDEPL,1,IDR)
+      CALL LCMPUT(IPLIB,'DEPLETE-ENER',NREAC*NDEPL,2,RER)
+      CALL LCMPUT(IPLIB,'DEPLETE-DECA',NDEPL,2,RRD)
+      CALL LCMPUT(IPLIB,'PRODUCE-REAC',NPAR*NDEPL,1,KPAR)
+      CALL LCMPUT(IPLIB,'PRODUCE-RATE',NPAR*NDEPL,2,BPAR)
+      IF(NDFI*NDFP.GT.0) THEN
+        CALL LCMPUT(IPLIB,'FISSIONYIELD',NDFI*NDFP*NBESP,2,YIELD)
+        IF(NBESP.GT.1) CALL LCMPUT(IPLIB,'ENERGY-YIELD',NBESP+1,2,ENER)
+      ENDIF
+      CALL LCMSIX(IPLIB,' ',2)
+      IF(IMPX.GE.2) WRITE(IOUT,6000) (ISTATE(ISTA),ISTA=1,10)
+*----
+*  RELEASE DECAY CHAIN
+*----
+      DEALLOCATE(HREAC)
+      IF(NDFP.GT.0) DEALLOCATE(YIELD)
+      DEALLOCATE(BPAR,KPAR,RER,RRD,IDR,IZAE,INAM)
+      IF(IEVOT.GT.1) DEALLOCATE(ENER)
+      RETURN
+*----
+*  FORMAT
+*----
+ 6000 FORMAT(/' STATE-VECTOR FOR 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)'/
+     >  ' NBESP  ',I6,' (NUMBER OF ENERGY-DEPENDENT FISSION YIELD MAT',
+     >  'RICES)'/)
+      END