Initial commit from Polytechnique Montreal

1 files changed, 402 insertions, 0 deletions

diff --git a/Donjon/src/PCRDRV.f b/Donjon/src/PCRDRV.f
new file mode 100644
index 0000000..70468e2
--- /dev/null
+++ b/Donjon/src/PCRDRV.f
@@ -0,0 +1,402 @@
+*DECK PCRDRV
+      SUBROUTINE PCRDRV(LCUBIC,NMIX,IMPX,NCAL,ITER,MAXNIS,TERP,NISO,
+     1 HISO,CONC,LMIXC,XS_CALC)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute TERP factors for PMAXS file interpolation. Use user-defined
+* global and local parameters.
+*
+*Copyright:
+* Copyright (C) 2019 Ecole Polytechnique de Montreal
+*
+*Author(s): 
+* A. Hebert
+*
+*Parameters: input
+* LCUBIC  =.TRUE.: cubic Ceschino interpolation; =.FALSE: linear
+*         Lagrange interpolation.
+* NMIX    maximum number of material mixtures in the microlib.
+* IMPX    print parameter (equal to zero for no print).
+* NCAL    number of elementary calculations in the PMAXS file.
+*
+*Parameters: output
+* ITER    completion flag (=0: all over; =1: use another PMAXS file;
+*         =2 use another L_MAP + PMAXS file).
+* MAXNIS  maximum value of NISO(I) in user data.
+* TERP    interpolation factors.
+* NISO    number of user-selected isotopes.
+* HISO    name of the user-selected isotopes.
+* CONC    user-defined number density of the user-selected isotopes.
+* LMIXC   flag set to .true. for fuel-map mixtures to process.
+* XS_CALC pointers towards PMAXS elementary calculations.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+      USE PCRDATA
+      IMPLICIT NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER, PARAMETER::MAXISD=200
+      INTEGER NMIX,IMPX,NCAL,ITER,MAXNIS,NISO(NMIX),HISO(2,NMIX,MAXISD)
+      REAL TERP(NCAL,NMIX),CONC(NMIX,MAXISD)
+      LOGICAL LCUBIC,LMIXC(NMIX)
+      TYPE(XSBLOCK_ITEM) XS_CALC(NCAL)
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER, PARAMETER::IOUT=6
+      INTEGER, PARAMETER::MAXLIN=50
+      INTEGER, PARAMETER::MAXPAR=50
+      INTEGER, PARAMETER::MAXVAL=200
+      REAL, PARAMETER::REPS=1.0E-4
+      REAL FLOTT, SUM
+      INTEGER I0, IBM, ICAL, INDIC, IPAR, ITYPE, I, JBM, J, NCOMLI,
+     & NITMA, NPAR, IBRA, IBSET, II, IND, INDELT, NBURN, NNV
+      CHARACTER TEXT12*12,PARKEY(MAXPAR)*12,HSMG*131,COMMEN(MAXLIN)*80,
+     1 RECNAM*12,HCUBIC*12
+      INTEGER NVALUE(MAXPAR),MUPLET(MAXPAR),MUTYPE(MAXPAR)
+      DOUBLE PRECISION DFLOTT
+      REAL VALR(MAXPAR,2),VREAL(MAXVAL,MAXPAR)
+      LOGICAL LCUB2(MAXPAR)
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: MUBASE
+      LOGICAL, ALLOCATABLE, DIMENSION(:) :: LDELTA
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(LDELTA(NMIX))
+*----
+*  RECOVER TABLE-OF-CONTENT INFORMATION FOR THE PMAXS FILE. THE I-TH
+*  PMAXS FILE INFORMATION CORRESPONDS TO POINTERS bran_i and PMAX.
+*----
+      NPAR=bran_i%Nstat_var
+      NVALUE(:NPAR)=0
+      DO IPAR=1,bran_i%Nstat_var
+        PARKEY(IPAR)=bran_i%var_nam(IPAR)
+      ENDDO
+      IF(PMAX%NBset.GT.0) THEN
+        NPAR=NPAR+1
+        PARKEY(NPAR)='B'
+        NVALUE(NPAR)=PMAX%Bset(1)%NBURN
+        NNV=NVALUE(NPAR)
+        VREAL(:NNV,NPAR)=REAL(PMAX%Bset(1)%burns(:NNV))
+      ENDIF
+      IF(NPAR.GT.MAXPAR) CALL XABORT('PCRDRV: MAXPAR OVERFLOW.')
+      IF(NHST.NE.1) CALL XABORT('PCRDRV: MULTIPLE HISTORY CASE NOT IMP'
+     1 //'LEMENTED.')
+      NCOMLI=6
+      COMMEN(:6)=hcomment(:6)
+      DO IBRA=1,NBRA
+        DO IPAR=1,bran_i%Nstat_var
+          FLOTT=REAL(bran_i%state(IPAR,IBRA))
+          IF(NVALUE(IPAR).EQ.0) THEN
+            NVALUE(IPAR)=1
+            VREAL(1,IPAR)=FLOTT
+          ELSE
+            DO I=1,NVALUE(IPAR)
+              IF(FLOTT.EQ.VREAL(I,IPAR)) THEN
+                GO TO 10
+              ELSE IF(FLOTT.LT.VREAL(I,IPAR)) THEN
+                DO J=NVALUE(IPAR),I,-1
+                  VREAL(J+1,IPAR)=VREAL(J,IPAR)
+                ENDDO
+                VREAL(I,IPAR)=FLOTT
+                NVALUE(IPAR)=NVALUE(IPAR)+1
+                GO TO 10
+              ENDIF
+            ENDDO
+            IF(FLOTT.GT.VREAL(NVALUE(IPAR),IPAR)) THEN
+              NVALUE(IPAR)=NVALUE(IPAR)+1
+              VREAL(NVALUE(IPAR),IPAR)=FLOTT
+            ENDIF
+          ENDIF
+   10     CONTINUE
+        ENDDO
+      ENDDO
+      IF((IMPX.GT.0).AND.(bran_i%Nstat_var.GT.0))THEN
+        DO IPAR=1,NPAR
+          WRITE(RECNAM,'(''pval'',I8.8)') IPAR
+          WRITE(IOUT,'(13H PCRDRV: KEY=,A,18H TABULATED POINTS=,
+     1    1P,6E12.4/(43X,6E12.4))') PARKEY(IPAR),(VREAL(I,IPAR),I=1,
+     2    NVALUE(IPAR))
+        ENDDO
+      ENDIF
+*----
+*  PRINT PMAXS FILE AND FUELMAP STATISTICS
+*----
+      IF(IMPX.GT.0) THEN
+        WRITE(IOUT,'(43H PCRDRV: NUMBER OF CALCULATIONS IN PMAXS FI,
+     1  3HLE=,I6)') NCAL
+        WRITE(IOUT,'(43H PCRDRV: NUMBER OF MATERIAL MIXTURES IN FUE,
+     1  6HL MAP=,I6)') NMIX
+        WRITE(IOUT,'(43H PCRDRV: NUMBER OF LOCAL VARIABLES INCLUDIN,
+     1  9HG BURNUP=,I6)') NPAR
+        WRITE(IOUT,'(28H PCRDRV: PMAXS FILE COMMENTS,60(1H-))')
+        WRITE(IOUT,'(1X,A)') (COMMEN(I),I=1,NCOMLI)
+        WRITE(IOUT,'(9H PCRDRV: ,79(1H-))')
+      ENDIF
+*----
+*  SCAN THE PMAXS FILE INFORMATION TO RECOVER THE MUPLET DATABASE
+*----
+      IF(IMPX.GT.5) THEN
+        WRITE(IOUT,'(24H PCRDRV: MUPLET DATABASE/12H CALCULATION,4X,
+     1  6HMUPLET)')
+        WRITE(IOUT,'(16X,20A4)') PARKEY(:NPAR)
+      ENDIF
+      ALLOCATE(MUBASE(NPAR,NCAL))
+      ICAL=0
+      DO IBRA=1,NBRA
+        INDELT=0
+        DO IPAR=1,NPAR
+          IF(bran_i%state_nam(IBRA).EQ.PARKEY(IPAR)) THEN
+            INDELT=IPAR
+            CYCLE
+          ENDIF
+        ENDDO
+        IBSET=PMAX%BRANCH(IBRA,1)%IBSET
+        NBURN=PMAX%Bset(IBSET)%NBURN
+        DO IPAR=1,bran_i%Nstat_var
+          FLOTT=REAL(bran_i%state(IPAR,IBRA))
+          IND=0
+          DO I=1,NVALUE(IPAR)
+            IF(FLOTT.EQ.VREAL(I,IPAR)) THEN
+               IND=I
+               EXIT
+            ENDIF
+          ENDDO
+          IF(IND.EQ.0) THEN
+            CALL XABORT('PCRDRV: MUPLET ALGORITHM FAILURE.')
+          ELSE
+            MUPLET(IPAR)=IND
+          ENDIF
+        ENDDO
+        IF((NBURN.EQ.PMAX%Bset(1)%NBURN).OR.(NBURN.EQ.1)) THEN
+          DO I=1,NBURN
+            MUPLET(bran_i%Nstat_var+1)=I
+            II=ICAL+I
+            MUBASE(:bran_i%Nstat_var+1,II)=MUPLET(:bran_i%Nstat_var+1)
+            XS_CALC(ICAL+I)%IBURN=I
+            XS_CALC(ICAL+I)%XS=>PMAX%BRANCH(IBRA,1)%XS(I)
+            XS_CALC(ICAL+I)%TIV=>PMAX%TIVB(1)%TIV(I)
+            IF(INDELT.GT.0) THEN
+              XS_CALC(ICAL+I)%DELTA=bran_i%state(INDELT,IBRA)-
+     1        bran_i%state(INDELT,1)
+            ELSE
+              XS_CALC(ICAL+I)%DELTA=0.0
+            ENDIF
+          ENDDO
+        ELSE
+          CALL XABORT('PCRDRV: INVALID VALUE OF NBURN.')
+        ENDIF
+        IF(IMPX.GT.5) THEN
+          DO I=ICAL+1,ICAL+NBURN
+            WRITE(IOUT,'(I8,2X,A2,2X,20I4/(14X,20I4))') I,
+     1      bran_i%state_nam(IBRA),MUBASE(:NPAR,I)
+          ENDDO
+        ENDIF
+        ICAL=ICAL+NBURN
+      ENDDO !IBRA
+      IF(ICAL.NE.NCAL) CALL XABORT('PCRDRV: MUPLET ALGORITHM FAILURE.')
+*----
+*  READ (INTERP_DATA) AND SET VALR PARAMETERS CORRESPONDING TO THE
+*  INTERPOLATION POINT. FILL MUPLET FOR PARAMETERS SET WITHOUT
+*  INTERPOLATION.
+*----
+      NISO(:NMIX)=0
+      TERP(:NCAL,:NMIX)=0.0
+      LMIXC(:NMIX)=.FALSE.
+*----
+*  READ (INTERP_DATA) AND SET VALR PARAMETERS CORRESPONDING TO THE
+*  INTERPOLATION POINT. FILL MUPLET FOR PARAMETERS SET WITHOUT
+*  INTERPOLATION.
+*----
+      IBM=0
+      MAXNIS=0
+      NISO(:NMIX)=0
+      LDELTA(:NMIX)=.FALSE.
+   20 CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+      IF(INDIC.NE.3) CALL XABORT('PCRDRV: CHARACTER DATA EXPECTED.')
+   30 IF(TEXT12.EQ.'MIX') THEN
+         MUPLET(:NPAR)=0
+         MUTYPE(:NPAR)=0
+         VALR(:NPAR,1)=0.0
+         VALR(:NPAR,2)=0.0
+         LCUB2(:NPAR)=LCUBIC
+         CALL REDGET(INDIC,IBM,FLOTT,TEXT12,DFLOTT)
+         IF(INDIC.NE.1) CALL XABORT('PCRDRV: INTEGER DATA EXPECTED.')
+         IF(IBM.GT.NMIX) CALL XABORT('PCRDRV: NMIX OVERFLOW.')
+         LMIXC(IBM)=.TRUE.
+         CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+         IF(INDIC.NE.3) CALL XABORT('PCRDRV: CHARACTER DATA EXPECTED.')
+         GO TO 30
+      ELSE IF(TEXT12.EQ.'MICRO') THEN
+         IF(IBM.EQ.0) CALL XABORT('PCRDRV: MIX NOT SET (1).')
+         CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+         IF(INDIC.NE.3) CALL XABORT('PCRDRV: CHARACTER DATA EXPECTED.')
+   50    IF(TEXT12.EQ.'ENDMIX') THEN
+            GO TO 30
+         ELSE
+            NISO(IBM)=NISO(IBM)+1
+            IF(NISO(IBM).GT.MAXISD) CALL XABORT('PCRDRV: MAXISD OVERFL'
+     1      //'OW.')
+            MAXNIS=MAX(MAXNIS,NISO(IBM))
+            READ(TEXT12,'(2A4)') (HISO(I0,IBM,NISO(IBM)),I0=1,2)
+            CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+            IF(INDIC.EQ.2) THEN
+               CONC(IBM,NISO(IBM))=FLOTT
+            ELSE
+               CALL XABORT('PCRDRV: INVALID HISO DATA.')
+            ENDIF
+            CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+            IF(INDIC.NE.3) CALL XABORT('PCRDRV: CHARACTER DATA EXPECTE'
+     1      //'D.')
+            GO TO 50
+         ENDIF
+      ELSE IF((TEXT12.EQ.'SET').OR.(TEXT12.EQ.'DELTA')) THEN
+         IF(IBM.EQ.0) CALL XABORT('PCRDRV: MIX NOT SET (2).')
+         ITYPE=0
+         IF(TEXT12.EQ.'SET') THEN
+            ITYPE=1
+         ELSE IF(TEXT12.EQ.'DELTA') THEN
+            ITYPE=2
+            LDELTA(IBM)=.TRUE.
+         ENDIF
+         CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+         IF(INDIC.NE.3) CALL XABORT('PCRDRV: CHARACTER DATA EXPECTED.')
+         IF((TEXT12.EQ.'LINEAR').OR.(TEXT12.EQ.'CUBIC')) THEN
+            HCUBIC=TEXT12
+            CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+         ELSE
+            HCUBIC=' '
+         ENDIF
+         IF(INDIC.NE.3) CALL XABORT('PCRDRV: CHARACTER DATA EXPECTED.')
+         DO 60 I=1,NPAR
+         IF(TEXT12.EQ.PARKEY(I)) THEN
+            IPAR=I
+            GO TO 70
+         ENDIF
+   60    CONTINUE
+         CALL XABORT('PCRDRV: PARAMETER '//TEXT12//' NOT FOUND.')
+   70    IF(HCUBIC.EQ.'LINEAR') THEN
+            LCUB2(IPAR)=.FALSE.
+         ELSE IF(HCUBIC.EQ.'CUBIC') THEN
+            LCUB2(IPAR)=.TRUE.
+         ENDIF
+         CALL REDGET(INDIC,NITMA,VALR(IPAR,1),TEXT12,DFLOTT)
+         IF(INDIC.NE.2) CALL XABORT('PCRDRV: REAL DATA EXPECTED.')
+         VALR(IPAR,2)=VALR(IPAR,1)
+         CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+         IF(INDIC.EQ.2) THEN
+            VALR(IPAR,2)=FLOTT
+            CALL REDGET(INDIC,NITMA,FLOTT,TEXT12,DFLOTT)
+         ENDIF
+         IF(VALR(IPAR,1).EQ.VALR(IPAR,2)) THEN
+            DO 80 J=1,NVALUE(IPAR)
+            IF(ABS(VALR(IPAR,1)-VREAL(J,IPAR)).LE.REPS*
+     1      ABS(VREAL(J,IPAR)))THEN
+               MUPLET(IPAR)=J
+               IF(ITYPE.NE.1) MUPLET(IPAR)=-1
+               MUTYPE(IPAR)=ITYPE
+               GO TO 30
+            ENDIF
+   80       CONTINUE
+         ENDIF
+         IF(VALR(IPAR,1).LT.VREAL(1,IPAR)) THEN
+            WRITE(HSMG,'(23HPCRDRV: REAL PARAMETER ,A,10H WITH VALU,
+     1      1HE,1P,E12.4,25H IS OUTSIDE THE DOMAIN (<,E12.4,1H))')
+     2      PARKEY(IPAR),VALR(IPAR,1),VREAL(1,IPAR)
+            CALL XABORT(HSMG)
+         ELSE IF(VALR(IPAR,2).GT.VREAL(NVALUE(IPAR),IPAR)) THEN
+            WRITE(HSMG,'(23HPCRDRV: REAL PARAMETER ,A,10H WITH VALU,
+     1      1HE,1P,E12.4,25H IS OUTSIDE THE DOMAIN (>,E12.4,1H))')
+     2      PARKEY(IPAR),VALR(IPAR,1),VREAL(NVALUE(IPAR),IPAR)
+            CALL XABORT(HSMG)
+         ELSE IF(VALR(IPAR,1).GT.VALR(IPAR,2)) THEN
+            WRITE(HSMG,'(23HPCRDRV: REAL PARAMETER ,A,9H IS DEFIN,
+     1      7HED WITH,1P,E12.4,2H >,E12.4,1H.)') PARKEY(IPAR),
+     2      VALR(IPAR,1),VALR(IPAR,2)
+            CALL XABORT(HSMG)
+         ENDIF
+         MUPLET(IPAR)=-1
+         MUTYPE(IPAR)=ITYPE
+         GO TO 30
+      ELSE IF(TEXT12.EQ.'ENDMIX') THEN
+*----
+*  COMPUTE THE TERP FACTORS USING TABLE-OF-CONTENT INFORMATION.
+*----
+         IF(IMPX.GT.0) THEN
+           DO IPAR=1,NPAR
+             IF(LCUB2(IPAR)) THEN
+               WRITE(IOUT,'(26H PCRDRV: GLOBAL PARAMETER:,A12,5H ->CU,
+     1         18HBIC INTERPOLATION.)') PARKEY(IPAR)
+             ELSE
+               WRITE(IOUT,'(26H PCRDRV: GLOBAL PARAMETER:,A12,5H ->LI,
+     1         19HNEAR INTERPOLATION.)') PARKEY(IPAR)
+             ENDIF
+           ENDDO
+         ENDIF
+         IF(IBM.GT.NMIX) CALL XABORT('PCRDRV: MIX OVERFLOW (MICROLIB).')
+         IF(NCAL.EQ.1) THEN
+           TERP(1,IBM)=1.0
+         ELSE
+           CALL PCRTRP(LCUB2,IMPX,NPAR,NCAL,NVALUE,MUPLET,MUTYPE,VALR,
+     1     0.0,MUBASE,VREAL,TERP(1,IBM))
+         ENDIF
+         IBM=0
+      ELSE IF((TEXT12.EQ.'PMAXS').OR.(TEXT12.EQ.'TABLE').OR.
+     1   (TEXT12.EQ.';')) THEN
+*----
+*  CHECK TERP FACTORS AND RETURN
+*----
+         IF(TEXT12.EQ.';') ITER=0
+         IF(TEXT12.EQ.'PMAXS') ITER=1
+         IF(TEXT12.EQ.'TABLE') ITER=2
+         DO 150 IBM=1,NMIX
+         IF(.NOT.LMIXC(IBM)) GO TO 150
+         IF(NISO(IBM).GT.MAXNIS) CALL XABORT('PCRDRV: MAXNIS OVERFLOW.')
+         IF(LDELTA(IBM)) THEN
+            SUM=0.0
+         ELSE
+            SUM=1.0
+         ENDIF
+         DO 140 ICAL=1,NCAL
+         SUM=SUM-TERP(ICAL,IBM)
+  140    CONTINUE
+         IF(ABS(SUM).GT.1.0E-4) THEN
+            WRITE(HSMG,'(43HPCRDRV: INVALID INTERPOLATION FACTORS IN MI,
+     1      5HXTURE,I4,1H.)') IBM
+            CALL XABORT(HSMG)
+         ENDIF
+  150    CONTINUE
+         GO TO 160
+      ELSE
+         CALL XABORT('PCRDRV: '//TEXT12//' IS AN INVALID KEYWORD.')
+      ENDIF
+      GO TO 20
+*----
+*  PRINT INTERPOLATION (TERP) FACTORS
+*----
+  160 IF(IMPX.GT.2) THEN
+        WRITE(IOUT,'(/30H PCRDRV: INTERPOLATION FACTORS)')
+        DO ICAL=1,NCAL
+          DO IBM=1,NMIX
+            IF(TERP(ICAL,IBM).NE.0.0) THEN
+              WRITE(IOUT,170) ICAL,(TERP(ICAL,JBM),JBM=1,NMIX)
+              EXIT
+            ENDIF
+          ENDDO
+        ENDDO
+      ENDIF
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(MUBASE,LDELTA)
+      RETURN
+  170 FORMAT(6H CALC=,I8,6H TERP=,1P,8E13.5/(20X,8E13.5))
+      END