Initial commit from Polytechnique Montreal

1 files changed, 335 insertions, 0 deletions

diff --git a/Donjon/src/PCRMIC.f b/Donjon/src/PCRMIC.f
new file mode 100644
index 0000000..9219e44
--- /dev/null
+++ b/Donjon/src/PCRMIC.f
@@ -0,0 +1,335 @@
+*DECK PCRMIC
+      SUBROUTINE PCRMIC(MAXNIS,MAXISO,IPLIB,IACCS,NMIX,NGRP,IMPX,
+     1 NCAL,TERP,NISO,HISO,CONC,LMIXC,XS_CALC,B2)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Build the microlib by scanning the NCAL elementary calculations from
+* PMAXS file and weighting them with TERP factors.
+*
+*Copyright:
+* Copyright (C) 2019 Ecole Polytechnique de Montreal
+*
+*Author(s): 
+* A. Hebert
+*
+*Parameters: input
+* MAXNIS  maximum value of NISO(I) in user data.
+* MAXISO  maximum allocated space for output microlib TOC information.
+* IPLIB   address of the output microlib LCM object.
+* IACCS   =0 microlib is created; =1 ... is updated.
+* NMIX    maximum number of material mixtures in the microlib.
+* NGRP    number of energy groups.
+* IMPX    print parameter (equal to zero for no print).
+* NCAL    number of elementary calculations in the PMAXS file.
+* 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.
+* B2      buckling
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+      USE PCRDATA
+      IMPLICIT NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPLIB
+      INTEGER MAXNIS,MAXISO,IACCS,NMIX,NGRP,IMPX,NCAL,NISO(NMIX),
+     1 HISO(2,NMIX,MAXNIS)
+      REAL TERP(NCAL,NMIX),CONC(NMIX,MAXNIS),B2
+      LOGICAL LMIXC(NMIX)
+      TYPE(XSBLOCK_ITEM) XS_CALC(NCAL)
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER, PARAMETER::IOUT=6
+      INTEGER, PARAMETER::MAXED=50
+      INTEGER, PARAMETER::NSTATE=40
+      INTEGER I0, IBM, ICAL, IED1, IED2, IGR, ISO, ITRANC, KSO1, I,
+     & JSO, KSO, NBISO1, NBISO2, NCOMB2, NCOMB, NDEL, NDEPL, NED1,
+     & NED2, NL, ITSTMP, MAXMIX, NBISO
+      REAL WEIGHT,TMPDAY(3)
+      CHARACTER TEXT12*12,HNAME*12,HVECT1(MAXED)*8,HHISO*8,TEXT8*8,
+     & HVECT2(MAXED)*8
+      INTEGER ISTATE(NSTATE)
+      TYPE(C_PTR) IPTMP,JPTMP,KPTMP,JPLIB,KPLIB
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: IMIX2,ITOD2,ISTY2,MILVO,
+     & IMICR
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: HUSE1,HUSE2,HNAM2
+      REAL, ALLOCATABLE, DIMENSION(:) :: DENS2,ENER,GAR1,GAR2
+      LOGICAL, ALLOCATABLE, DIMENSION(:) :: MASKL
+      INTEGER, POINTER, DIMENSION(:) :: ISONA,ISOMI
+      REAL, POINTER, DIMENSION(:) :: DENIS
+      TYPE(C_PTR) ISONA_PTR,ISOMI_PTR,DENIS_PTR
+      TYPE(C_PTR), ALLOCATABLE, DIMENSION(:) :: IPLIST
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(HUSE1(3,MAXISO),IMIX2(MAXISO),ITOD2(MAXISO),
+     & ISTY2(MAXISO),HUSE2(3,MAXISO),HNAM2(3,MAXISO),MILVO(NMIX),
+     & IPLIST(MAXISO))
+      ALLOCATE(DENS2(MAXISO),ENER(NGRP+1))
+*----
+*  MICROLIB INITIALIZATION
+*----
+      ITRANC=0
+      DENS2(:MAXISO)=0.0
+      IMIX2(:MAXISO)=0
+      ITOD2(:MAXISO)=0
+      ISTY2(:MAXISO)=3
+      IPLIST(:MAXISO)=C_NULL_PTR
+      IF(IACCS.EQ.0) THEN
+         NBISO2=0
+         NCOMB2=0
+         NED2=0
+         TEXT12='L_LIBRARY'
+         CALL LCMPTC(IPLIB,'SIGNATURE',12,TEXT12)
+      ELSE
+         CALL LCMGET(IPLIB,'STATE-VECTOR',ISTATE)
+         IF(ISTATE(1).NE.NMIX) CALL XABORT('PCRMIC: INVALID NUMBER OF '
+     1   //'MATERIAL MIXTURES IN THE MICROLIB.')
+         IF(ISTATE(3).NE.NGRP) CALL XABORT('PCRMIC: INVALID NUMBER OF '
+     1   //'ENERGY GROUPS IN THE MICROLIB.')
+         NBISO2=ISTATE(2)
+         NCOMB2=ISTATE(12)
+         IF(NBISO2.GT.MAXISO) CALL XABORT('PCRMIC: MAXISO OVERFLOW(1).')
+         NED2=ISTATE(13)
+         IF(NED2.GT.MAXED) CALL XABORT('PCRMIC: MAXED OVERFLOW.')
+         CALL LCMGET(IPLIB,'ISOTOPESUSED',HUSE2)
+         CALL LCMGET(IPLIB,'ISOTOPERNAME',HNAM2)
+         CALL LCMGET(IPLIB,'ISOTOPESDENS',DENS2)
+         CALL LCMGET(IPLIB,'ISOTOPESMIX',IMIX2)
+         CALL LCMGET(IPLIB,'ISOTOPESTODO',ITOD2)
+         CALL LCMGET(IPLIB,'ISOTOPESTYPE',ISTY2)
+         IF(NED2.GT.0) CALL LCMGTC(IPLIB,'ADDXSNAME-P0',8,NED2,HVECT2)
+         CALL LCMGET(IPLIB,'ENERGY',ENER)
+      ENDIF
+*----
+*  LOOP OVER MICROLIB MIXTURES
+*----
+      MILVO(:NMIX)=0
+      NCOMB=0
+      DO 190 IBM=1,NMIX
+      IF(.NOT.LMIXC(IBM)) GO TO 190
+      IF(NISO(IBM).GT.MAXNIS) CALL XABORT('PCRMIC: MAXNIS OVERFLOW.')
+*----
+*  FIND THE VALUE OF NBISO1 AND HUSE1 IN MIXTURE IBM
+*----
+      NBISO1=1
+      TEXT12='*MAC*RES'
+      READ(TEXT12,'(3A4)') (HUSE1(I,1),I=1,3)
+      IF(NXST.GT.4) THEN
+        NBISO1=3
+        TEXT12='Xe135'
+        READ(TEXT12,'(3A4)') (HUSE1(I,2),I=1,3)
+        TEXT12='Sm149'
+        READ(TEXT12,'(3A4)') (HUSE1(I,3),I=1,3)
+      ENDIF
+*----
+*  LOOP OVER ELEMENTARY CALCULATIONS
+*----
+      CALL LCMOP(IPTMP,'*CALCULATIONS*',0,1,0)
+      JPTMP=LCMLID(IPTMP,'CALCULATIONS',NCAL)
+      DO 70 ICAL=1,NCAL
+      WEIGHT=TERP(ICAL,IBM)
+      IF(WEIGHT.EQ.0.0) GO TO 70
+      KPTMP=LCMDIL(JPTMP,ICAL)
+      CALL PCRONE(IMPX,ICAL,KPTMP,NCAL,NGRP,XS_CALC)
+      IF(IMPX.GT.0) THEN
+        WRITE(IOUT,'(33H PCRMIC: PMAXS ACCESS FOR MIXTURE,I8,6H AND C,
+     1  10HALCULATION,I8,9H. WEIGHT=,1P,E12.4)') IBM,ICAL,WEIGHT
+        IF(IMPX.GT.50) CALL LCMLIB(KPTMP)
+      ENDIF
+      CALL LCMGET(KPTMP,'STATE-VECTOR',ISTATE)
+      IF(ISTATE(1).NE.1) CALL XABORT('PCRMIC: INVALID NUMBER OF MATERI'
+     1 //'AL MIXTURES IN THE PMAXS FILE.')
+      IF(ISTATE(2).NE.NBISO1) CALL XABORT('PCRMIC: INVALID NBISO1.')
+      IF(ISTATE(3).NE.NGRP) CALL XABORT('PCRMIC: INVALID NUMBER OF ENE'
+     1 //'RGY GROUPS IN THE COMPO.')
+      NL=ISTATE(4)
+      ITRANC=ISTATE(5)
+      NDEPL=0
+      NED1=ISTATE(13)
+      NDEL=ISTATE(19)
+      IF(NED1.GT.MAXED) CALL XABORT('PCRMIC: MAXED OVERFLOW.')
+      IF(NED1.GT.0) CALL LCMGTC(KPTMP,'ADDXSNAME-P0',8,NED1,HVECT1)
+      CALL LCMGET(KPTMP,'ENERGY',ENER)
+      DO 30 IED1=1,NED1
+      DO 20 IED2=1,NED2
+      IF(HVECT1(IED1).EQ.HVECT2(IED2)) GO TO 30
+   20 CONTINUE
+      NED2=NED2+1
+      HVECT2(NED2)=HVECT1(IED1)
+   30 CONTINUE
+      CALL LCMGPD(KPTMP,'ISOTOPESUSED',ISONA_PTR)
+      CALL LCMGPD(KPTMP,'ISOTOPESDENS',DENIS_PTR)
+      CALL C_F_POINTER(ISONA_PTR,ISONA,(/ NBISO1 /))
+      CALL C_F_POINTER(DENIS_PTR,DENIS,(/ NBISO1 /))
+      DO 60 ISO=1,NBISO1
+      WRITE(TEXT8,'(2A4)') (ISONA(3*(ISO-1)+I0),I0=1,2)
+      IF(TEXT8.EQ.'*MAC*RES') THEN
+        DENIS(ISO)=1.0
+      ELSE
+        KSO1=0
+        DO 40 KSO=1,NISO(IBM) ! user-selected isotope
+        WRITE(HHISO,'(2A4)') (HISO(I0,IBM,KSO),I0=1,2)
+        IF(TEXT8.EQ.HHISO) THEN
+          KSO1=KSO
+          GO TO 50
+        ENDIF
+   40   CONTINUE
+   50   IF(KSO1.GT.0) DENIS(ISO)=CONC(IBM,KSO1)
+      ENDIF
+   60 CONTINUE
+      CALL LCMPPD(KPTMP,'ISOTOPESDENS',NBISO1,2,DENIS_PTR)
+   70 CONTINUE
+*----
+*  SELECT MICROLIB ISOTOPES CORRESPONDING TO PMAXS ISOTOPES
+*----
+      DO 90 ISO=1,NBISO1 ! PMAXS isotope
+      WRITE(TEXT12,'(2A4)') (HUSE1(I,ISO),I=1,2)
+      NBISO2=NBISO2+1
+      IF(NBISO2.GT.MAXISO) THEN
+        WRITE(IOUT,'(/16H PCRMIC: NBISO2=,I6,8H MAXISO=,I6)') NBISO2,
+     1  MAXISO
+        CALL XABORT('PCRMIC: MAXISO OVERFLOW(2).')
+      ENDIF
+      READ(TEXT12,'(3A4)') (HUSE2(I0,NBISO2),I0=1,3)
+      DO 80 I0=1,3
+      HNAM2(I0,NBISO2)=HUSE1(I0,ISO)
+   80 CONTINUE
+      IMIX2(NBISO2)=IBM
+      DENS2(NBISO2)=DENIS(ISO)
+   90 CONTINUE
+      ALLOCATE(IMICR(NBISO1))
+      IMICR(:NBISO1)=0
+      DO 130 ISO=1,NBISO2 ! microlib isotope
+      IF(IMIX2(ISO).NE.IBM) GO TO 130
+      DO 120 JSO=1,NBISO1 ! PMAXS isotope
+      IF((HUSE1(1,JSO).EQ.HUSE2(1,ISO)).AND.(HUSE1(2,JSO).EQ.
+     1 HUSE2(2,ISO))) THEN
+        IMICR(JSO)=ISO
+        GO TO 130
+      ENDIF
+  120 CONTINUE
+      WRITE(TEXT12,'(3A4)') (HUSE2(I0,ISO),I0=1,3)
+      CALL XABORT('PCRMIC: UNABLE TO FIND '//TEXT12//'.')
+  130 CONTINUE
+*----
+*  PROCESS ISOTOPE DIRECTORIES FOR MICROLIB MIXTURE IBM
+*----
+      DO 180 JSO=1,NBISO1 ! PMAXS isotope
+      ISO=IMICR(JSO) ! microlib isotope
+      IF(ISO.EQ.0) GO TO 180
+      WRITE(HNAME,'(3A4)') (HUSE1(I0,JSO),I0=1,3)
+      CALL LCMOP(KPLIB,'*ISOTOPE*',0,1,0)
+      IPLIST(ISO)=KPLIB ! set isot ISO
+      CALL PCRISO(KPLIB,JPTMP,HNAME,JSO,NCAL,NGRP,NL,NED2,HVECT2,NDEL,
+     1 IMPX,TERP(1,IBM))
+  180 CONTINUE
+      DEALLOCATE(IMICR)
+      CALL LCMCL(IPTMP,2)
+  190 CONTINUE
+*----
+*  END OF LOOP OVER MICROLIB MIXTURES
+*----
+*----
+*  CREATE ISOTOPE LIST DIRECTORY IN MICROLIB
+*----
+      JPLIB=LCMLID(IPLIB,'ISOTOPESLIST',NBISO2)
+      DO 195 ISO=1,NBISO2 ! microlib isotope
+      IF(C_ASSOCIATED(IPLIST(ISO))) THEN
+        KPLIB=LCMDIL(JPLIB,ISO) ! step up isot ISO
+        CALL LCMEQU(IPLIST(ISO),KPLIB)
+        CALL LCMCL(IPLIST(ISO),2)
+        IPLIST(ISO)=C_NULL_PTR
+      ENDIF
+  195 CONTINUE
+*----
+*  MICROLIB FINALIZATION
+*----
+      ISTATE(:NSTATE)=0
+      ISTATE(1)=NMIX
+      ISTATE(2)=NBISO2
+      ISTATE(3)=NGRP
+      ISTATE(4)=NL
+      ISTATE(5)=ITRANC
+      ISTATE(7)=1
+      ISTATE(11)=NDEPL
+      ISTATE(12)=NCOMB+NCOMB2
+      ISTATE(13)=NED2
+      ISTATE(14)=NMIX
+      ISTATE(18)=1
+      ISTATE(19)=NDEL
+      ISTATE(22)=MAXISO/NMIX
+      IF(NBISO2.EQ.0) CALL XABORT('PCRMIC: NBISO2=0.')
+      CALL LCMPUT(IPLIB,'STATE-VECTOR',NSTATE,1,ISTATE)
+      CALL LCMPUT(IPLIB,'ISOTOPESUSED',3*NBISO2,3,HUSE2)
+      CALL LCMPUT(IPLIB,'ISOTOPERNAME',3*NBISO2,3,HNAM2)
+      CALL LCMPUT(IPLIB,'ISOTOPESDENS',NBISO2,2,DENS2)
+      CALL LCMPUT(IPLIB,'ISOTOPESMIX',NBISO2,1,IMIX2)
+      IF(NED2.GT.0) CALL LCMPTC(IPLIB,'ADDXSNAME-P0',8,NED2,HVECT2)
+      CALL LCMPUT(IPLIB,'ISOTOPESTODO',NBISO2,1,ITOD2)
+      CALL LCMPUT(IPLIB,'ISOTOPESTYPE',NBISO2,1,ISTY2)
+      CALL LCMPUT(IPLIB,'ENERGY',NGRP+1,2,ENER)
+      IF(IMPX.GT.5) CALL LCMLIB(IPLIB)
+      IACCS=1
+*----
+*  BUILD EMBEDDED MACROLIB
+*----
+      CALL LCMGET(IPLIB,'STATE-VECTOR',ISTATE)
+      MAXMIX=ISTATE(1)
+      IF(MAXMIX.NE.NMIX) CALL XABORT('PCRMIC: INVALID NMIX.')
+      NBISO=ISTATE(2)
+      ALLOCATE(MASKL(NGRP))
+      CALL LCMGPD(IPLIB,'ISOTOPESUSED',ISONA_PTR)
+      CALL LCMGPD(IPLIB,'ISOTOPESMIX',ISOMI_PTR)
+      CALL LCMGPD(IPLIB,'ISOTOPESDENS',DENIS_PTR)
+      CALL C_F_POINTER(ISONA_PTR,ISONA,(/ NBISO /))
+      CALL C_F_POINTER(ISOMI_PTR,ISOMI,(/ NBISO /))
+      CALL C_F_POINTER(DENIS_PTR,DENIS,(/ NBISO /))
+      MASKL(:NGRP)=.TRUE.
+      ITSTMP=0
+      TMPDAY(1)=0.0
+      TMPDAY(2)=0.0
+      TMPDAY(3)=0.0
+      CALL LIBMIX(IPLIB,MAXMIX,NGRP,NBISO,ISONA,ISOMI,DENIS,LMIXC,MASKL,
+     1 ITSTMP,TMPDAY)
+      DEALLOCATE(MASKL)
+*----
+*  INCLUDE LEAKAGE IN THE MACROLIB (USED ONLY FOR NON-REGRESSION TESTS)
+*----
+      IF(B2.NE.0.0) THEN
+        IF(IMPX.GT.0) WRITE(6,'(/34H PCRMIC: INCLUDE LEAKAGE IN THE MA,
+     1  11HCROLIB (B2=,1P,E12.5,2H).)') B2
+        CALL LCMSIX(IPLIB,'MACROLIB',1)
+        JPLIB=LCMGID(IPLIB,'GROUP')
+        ALLOCATE(GAR1(NMIX),GAR2(NMIX))
+        DO 210 IGR=1,NGRP
+          KPLIB=LCMGIL(JPLIB,IGR)
+          CALL LCMGET(KPLIB,'NTOT0',GAR1)
+          CALL LCMGET(KPLIB,'DIFF',GAR2)
+          DO 200 IBM=1,NMIX
+            IF(LMIXC(IBM)) GAR1(IBM)=GAR1(IBM)+B2*GAR2(IBM)
+  200     CONTINUE
+          CALL LCMPUT(KPLIB,'NTOT0',NMIX,2,GAR1)
+  210   CONTINUE
+        DEALLOCATE(GAR2,GAR1)
+        CALL LCMSIX(IPLIB,' ',2)
+      ENDIF
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(ENER,DENS2)
+      DEALLOCATE(IPLIST,MILVO,HNAM2,HUSE2,ISTY2,ITOD2,IMIX2,HUSE1)
+      RETURN
+      END