Initial commit from Polytechnique Montreal

1 files changed, 369 insertions, 0 deletions

diff --git a/Dragon/src/BREMAC.f b/Dragon/src/BREMAC.f
new file mode 100644
index 0000000..eaba3e9
--- /dev/null
+++ b/Dragon/src/BREMAC.f
@@ -0,0 +1,369 @@
+*DECK BREMAC
+      SUBROUTINE BREMAC(NC,IPMAC2,NG,NL,LX1,NMIX1,NMIX2,IMIX,IMIX1,
+     1 IGAP,ILEAKS,IDF,IPRINT,ZKEFF,B2,VOL1,FLX1,DC1,TOT1,CHI1,SIGF1,
+     2 SCAT1,JXM,JXP,FHETXM,FHETXP,ADF1)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Recover macroscopic cross sections.
+*
+*Copyright:
+* Copyright (C) 2021 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
+* NC      number of sn macrolibs.
+* IPMAC2  pointer to the sn macrolib.
+* NG      number of energy groups.
+* NL      Legendre order of TOT1 and SCAT1 arrays (=1 for isotropic
+*         scattering in LAB).
+* LX1     number of nodes in the reflector model.
+* NMIX1   number of mixtures in the nodal calculation.
+* NMIX2   number of mixtures in the sn calculation after edition.
+* IMIX    mix index of each node in output data.
+* IMIX1   mix index of each node in sn editions.
+* IGAP    mix index of the right gap where the surface flux is
+*         recovered.
+* ILEAKS  type of leakage calculation (=0: no; =1: isotropic;
+*         =2: anisotropic).
+* IDF     discontinuity factor flag (=0: not used; =3: recovered).
+* IPRINT  print parameter
+*
+*Parameters: output
+* ZKEFF   effective multiplication factor.
+* B2      buckling.
+* VOL1    volumes.
+* FLX1    averaged fluxes
+* DC1     diffusion coefficients.
+* TOT1    total cross sections.
+* CHI1    fission spectra.
+* SIGF1   nu*fission cross sections.
+* SCAT1   scattering cross sections.
+* JXM     left boundary currents.
+* JXP     right boundary currents.
+* FHETXM  left boundary fluxes.
+* FHETXP  right boundary fluxes.
+* ADF1    assembly discontinuity factors.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER NC
+      TYPE(C_PTR) IPMAC2(NC)
+      INTEGER NG,NL,LX1,NMIX1,NMIX2,IMIX(LX1),IMIX1(LX1),IGAP(LX1),
+     1 ILEAKS,IDF,IPRINT
+      REAL ZKEFF(NC),B2(NC),VOL1(NMIX1,NC),FLX1(NMIX1,NG,NC),
+     1 DC1(NMIX1,NG,NC),TOT1(NMIX1,NG,NL,NC),CHI1(NMIX1,NG,NC),
+     2 SIGF1(NMIX1,NG,NC),SCAT1(NMIX1,NG,NG,NL,NC),JXM(NMIX1,NG,NC),
+     3 JXP(NMIX1,NG,NC),FHETXM(NMIX1,NG,NL,NC),FHETXP(NMIX1,NG,NL,NC),
+     4 ADF1(NMIX1,NG,NC)
+*----
+*  LOCAL VARIABLES
+*----
+      TYPE(C_PTR) JPMAC2,KPMAC2
+      DOUBLE PRECISION DSFIS
+      CHARACTER CM*2,HADF*8,TEXT12*12
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: IJJ,NJJ,IPOS,IMIX2
+      REAL, ALLOCATABLE, DIMENSION(:) :: VOL,WORK,SFIS,SFIS1
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: DC,CHI,SIGF
+      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: FLX,TOT
+      REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: SCAT
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: DCOU
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(VOL(NMIX2),FLX(NMIX2,NG,NL),TOT(NMIX2,NG,NL),
+     1 DC(NMIX2,NG),CHI(NMIX2,NG),SIGF(NMIX2,NG),SCAT(NMIX2,NG,NG,NL),
+     2 IMIX2(NMIX2),SFIS(NMIX2),SFIS1(NMIX1),DCOU(NMIX2+1,NG))
+      ALLOCATE(IJJ(NMIX2),NJJ(NMIX2),IPOS(NMIX2),WORK(NG*NMIX2))
+*----
+*  SET IMIX2
+*----
+      IMIX2(:)=0
+      DO J=1,LX1
+        IBM=IMIX1(J)
+        IBG=IGAP(J)
+        IF(IBM.GT.NMIX2) CALL XABORT('BREMAC: NMIX2 OVERFLOW(1).')
+        IF(IBG.GT.NMIX2) CALL XABORT('BREMAC: NMIX2 OVERFLOW(2).')
+        IF(IBM.GT.0) IMIX2(IBM)=IMIX(J)
+        IF((IBG.GT.0).AND.(J.LT.LX1)) IMIX2(IBG)=IMIX(J+1)
+      ENDDO
+*----
+*  LOOP OVER SN MACROLIBS
+*----
+      DO IC=1,NC
+*----
+*  RECOVER FLUX, MACROSCOPIC CROSS SECTIONS AND DIFFUSION COEFFICIENTS
+*----
+        CALL LCMGET(IPMAC2(IC),'VOLUME',VOL)
+        CALL LCMGET(IPMAC2(IC),'K-EFFECTIVE',ZKEFF(IC))
+        B2(IC)=0.0
+        IF(ILEAKS.GT.0) THEN
+          CALL LCMLEN(IPMAC2(IC),'B2  B1HOM',ILONG,ITYLCM)
+          IF(ILONG.EQ.1) CALL LCMGET(IPMAC2(IC),'B2  B1HOM',B2(IC))
+        ENDIF
+        JPMAC2=LCMGID(IPMAC2(IC),'GROUP')
+        SCAT(:,:,:,:)=0.0
+        DO IGR=1,NG
+          KPMAC2=LCMGIL(JPMAC2,IGR)
+          CALL LCMGET(KPMAC2,'FLUX-INTG',FLX(1,IGR,1))
+          CALL LCMGET(KPMAC2,'NTOT0',TOT(1,IGR,1))
+          IF(NL.GE.2) THEN
+            CALL LCMLEN(KPMAC2,'FLUX-INTG-P1',ILCMLN,ITYLCM)
+            IF(ILCMLN.GT.0) THEN
+              CALL LCMGET(KPMAC2,'FLUX-INTG-P1',FLX(1,IGR,2))
+            ELSE
+              FLX(:NMIX2,IGR,2)=FLX(:NMIX2,IGR,1)
+            ENDIF
+            CALL LCMLEN(KPMAC2,'NTOT1',ILCMLN,ITYLCM)
+            IF(ILCMLN.GT.0) THEN
+              CALL LCMGET(KPMAC2,'NTOT1',TOT(1,IGR,2))
+            ELSE
+              TOT(:NMIX2,IGR,2)=TOT(:NMIX2,IGR,1)
+            ENDIF
+            DO IL=3,NL
+              WRITE(TEXT12,'(11HFLUX-INTG-P,I1)') IL-1
+              CALL LCMLEN(KPMAC2,TEXT12,ILCMLN,ITYLCM)
+              IF(ILCMLN.GT.0) THEN
+                CALL LCMGET(KPMAC2,TEXT12,FLX(1,IGR,IL))
+              ELSE
+                FLX(:NMIX2,IGR,IL)=FLX(:NMIX2,IGR,IL-2)
+              ENDIF
+              WRITE(TEXT12,'(4HNTOT,I1)') IL-1
+              CALL LCMLEN(KPMAC2,TEXT12,ILCMLN,ITYLCM)
+              IF(ILCMLN.GT.0) THEN
+                CALL LCMGET(KPMAC2,TEXT12,TOT(1,IGR,IL))
+              ELSE
+                TOT(:NMIX2,IGR,IL)=TOT(:NMIX2,IGR,IL-2)
+              ENDIF
+            ENDDO
+          ENDIF
+          CALL LCMGET(KPMAC2,'DIFF',DC(1,IGR))
+          CALL LCMGET(KPMAC2,'CHI',CHI(1,IGR))
+          CALL LCMGET(KPMAC2,'NUSIGF',SIGF(1,IGR))
+          DO IL=1,NL
+            WRITE(CM,'(I2.2)') IL-1
+            CALL LCMLEN(KPMAC2,'IJJS'//CM,ILONG,ITYLCM)
+            IF(ILONG.EQ.0) CYCLE
+            CALL LCMGET(KPMAC2,'IJJS'//CM,IJJ)
+            CALL LCMGET(KPMAC2,'NJJS'//CM,NJJ)
+            CALL LCMGET(KPMAC2,'IPOS'//CM,IPOS)
+            CALL LCMGET(KPMAC2,'SCAT'//CM,WORK)
+            DO IBM=1,NMIX2
+              IPOSDE=IPOS(IBM)-1
+              DO JGR=IJJ(IBM),IJJ(IBM)-NJJ(IBM)+1,-1
+                IPOSDE=IPOSDE+1
+                SCAT(IBM,IGR,JGR,IL)=WORK(IPOSDE) ! IGR <-- JGR
+              ENDDO
+            ENDDO
+          ENDDO
+          DO IBM=1,NMIX2
+            FLX(IBM,IGR,:NL)=FLX(IBM,IGR,:NL)/VOL(IBM)
+          ENDDO
+        ENDDO
+*----
+*  COMPUTE NET CURRENTS BETWEEN NODES USING A BALANCE RELATION
+*----
+        DCOU(:NMIX2+1,:NG)=0.0D0
+        DO IBM=1,NMIX2
+          DSFIS=0.0D0
+          DO IGR=1,NG
+            DSFIS=DSFIS+SIGF(IBM,IGR)*FLX(IBM,IGR,1)
+          ENDDO
+          DSFIS=DSFIS/ZKEFF(IC)
+          DO IGR=1,NG
+            DCOU(IBM+1,IGR)=DCOU(IBM,IGR)+VOL(IBM)*(CHI(IBM,IGR)*DSFIS-
+     1      (TOT(IBM,IGR,1)+B2(IC)*DC(IBM,IGR))*FLX(IBM,IGR,1))
+            DO JGR=1,NG
+              DCOU(IBM+1,IGR)=DCOU(IBM+1,IGR)+VOL(IBM)*
+     1        SCAT(IBM,IGR,JGR,1)*FLX(IBM,JGR,1)
+            ENDDO
+          ENDDO
+        ENDDO
+*----
+*  NORMALIZE THE ODD FLUX UNKNOWNS BETWEEN NODES
+*----
+        IF(NL.GT.1) THEN
+          DO J=1,LX1
+            IBM=IMIX1(J)
+            IBG=IGAP(J)
+            IF((IBM.GT.0).AND.(IBG.GT.0)) THEN
+              DO IGR=1,NG
+                FACT=REAL(DCOU(IBM+1,IGR))/FLX(IBG,IGR,2)
+                DO IL=2,NL,2
+                  FLX(IBG,IGR,IL)=FLX(IBG,IGR,IL)*FACT
+                ENDDO
+              ENDDO
+            ENDIF
+          ENDDO
+        ENDIF
+*----
+*  HOMOGENIZATION OVER THE GAPS AND NODES
+*----
+        VOL1(:,IC)=0.0
+        SFIS1(:)=0.0
+        FLX1(:,:,IC)=0.0
+        TOT1(:,:,:,IC)=0.0
+        DC1(:,:,IC)=0.0
+        SIGF1(:,:,IC)=0.0
+        CHI1(:,:,IC)=0.0
+        SCAT1(:,:,:,:,IC)=0.0
+        DO IL=1,NL,2
+          FHETXM(:,:,IL,IC)=1.0
+          FHETXP(:,:,IL,IC)=1.0
+        ENDDO
+        DO IL=2,NL,2
+          FHETXM(:,:,IL,IC)=0.0
+          FHETXP(:,:,IL,IC)=0.0
+        ENDDO
+        JXM(:,:,IC)=0.0
+        JXP(:,:,IC)=0.0
+        ADF1(:,:,IC)=0.0
+        DO J=1,LX1
+          IBM=IMIX1(J)
+          IBG=IGAP(J)
+          IF(IBG.GT.0) THEN
+            DO IL=1,NL
+              IF(IMIX(J).GT.0) FHETXP(IMIX(J),:,IL,IC)=FLX(IBG,:,IL)
+              IF(J.LT.LX1) THEN
+                IF(IMIX(J+1).GT.0) THEN
+                  FHETXM(IMIX(J+1),:,IL,IC)=FLX(IBG,:,IL)
+                ENDIF
+              ENDIF
+            ENDDO
+          ENDIF
+          IF(IBM.GT.0) THEN
+            IF(IMIX(J).GT.0) THEN
+              JXM(IMIX(J),:NG,IC)=REAL(DCOU(MAX(1,IBM-1),:NG))
+              JXP(IMIX(J),:NG,IC)=REAL(DCOU(IBM+1,:NG))
+            ENDIF
+          ENDIF
+        ENDDO
+        DO IBM2=1,NMIX2
+          IBM=IMIX2(IBM2)
+          IF(IBM.EQ.0) CYCLE
+          VOL1(IBM,IC)=VOL1(IBM,IC)+VOL(IBM2)
+          SFIS(IBM2)=0.0
+          DO IGR=1,NG
+            SFIS(IBM2)=SFIS(IBM2)+VOL(IBM2)*FLX(IBM2,IGR,1)*
+     1      SIGF(IBM2,IGR)
+          ENDDO
+          SFIS1(IBM)=SFIS1(IBM)+SFIS(IBM2)
+          DO IGR=1,NG
+            FLX1(IBM,IGR,IC)=FLX1(IBM,IGR,IC)+VOL(IBM2)*FLX(IBM2,IGR,1)
+            DO IL=1,NL
+              TOT1(IBM,IGR,IL,IC)=TOT1(IBM,IGR,IL,IC)+VOL(IBM2)*
+     1        FLX(IBM2,IGR,1)*TOT(IBM2,IGR,IL)
+            ENDDO
+            DC1(IBM,IGR,IC)=DC1(IBM,IGR,IC)+VOL(IBM2)*FLX(IBM2,IGR,1)*
+     1      DC(IBM2,IGR)
+            SIGF1(IBM,IGR,IC)=SIGF1(IBM,IGR,IC)+VOL(IBM2)*
+     1      FLX(IBM2,IGR,1)*SIGF(IBM2,IGR)
+            CHI1(IBM,IGR,IC)=CHI1(IBM,IGR,IC)+SFIS(IBM2)*CHI(IBM2,IGR)
+            DO IL=1,NL
+              DO JGR=1,NG
+                SCAT1(IBM,IGR,JGR,IL,IC)=SCAT1(IBM,IGR,JGR,IL,IC)+
+     1          VOL(IBM2)*FLX(IBM2,JGR,1)*SCAT(IBM2,IGR,JGR,IL)
+              ENDDO
+            ENDDO
+          ENDDO
+        ENDDO
+        DO IBM=1,NMIX1
+          DO IGR=1,NG
+            DO IL=1,NL
+              TOT1(IBM,IGR,IL,IC)=TOT1(IBM,IGR,IL,IC)/FLX1(IBM,IGR,IC)
+            ENDDO
+            DC1(IBM,IGR,IC)=DC1(IBM,IGR,IC)/FLX1(IBM,IGR,IC)
+            IF(SFIS1(IBM).NE.0.0) CHI1(IBM,IGR,IC)=CHI1(IBM,IGR,IC)/
+     1      SFIS1(IBM)
+            SIGF1(IBM,IGR,IC)=SIGF1(IBM,IGR,IC)/FLX1(IBM,IGR,IC)
+            DO JGR=1,NG
+              DO IL=1,NL
+                SCAT1(IBM,IGR,JGR,IL,IC)=SCAT1(IBM,IGR,JGR,IL,IC)/
+     1          FLX1(IBM,JGR,IC)
+              ENDDO
+            ENDDO
+          ENDDO
+          DO IGR=1,NG
+            FLX1(IBM,IGR,IC)=FLX1(IBM,IGR,IC)/VOL1(IBM,IC)
+          ENDDO
+        ENDDO
+*----
+*  RECOVER ADF
+*----
+        IF(IDF.EQ.3) THEN
+          CALL LCMSIX(IPMAC2(IC),'ADF',1)
+          IF(IPRINT.GT.5) CALL LCMLIB(IPMAC2(IC))
+          CALL LCMLEN(IPMAC2(IC),'HADF',NTYPE,ITYLCM)
+          IF(NTYPE/2.NE.1) CALL XABORT('BREMAC: NTYPE=1 EXPECTED.')
+          CALL LCMGTC(IPMAC2(IC),'HADF',8,HADF)
+          CALL LCMLEN(IPMAC2(IC),HADF,ILONG,ITYLCM)
+          IF(ILONG.NE.NMIX1*NG) CALL XABORT('BREMAC: ADF OVERFLOW.')
+          CALL LCMGET(IPMAC2(IC),HADF,ADF1(1,1,IC))
+        ENDIF
+*----
+*  PRINT CROSS SECTIONS
+*----
+        IF(IPRINT.GT.1) THEN
+          WRITE(6,'(/36H BREMAC: CROSS SECTION FOR MACROLIB=,I5)') IC
+          WRITE(6,'(31H BREMAC: SCATTERING ANISOTROPY=,I5)') NL-1
+          WRITE(6,'(/6H KEFF=,1P E13.5,4H B2=,E13.5/)') ZKEFF(IC),B2(IC)
+          WRITE(6,10) 'IMIX',IMIX(:)
+          WRITE(6,20) 'VOL1',VOL1(:,IC)
+          WRITE(6,20) 'SFIS1',SFIS1(:)
+          DO IGR=1,NG
+            WRITE(6,'(/29H BREMAC: PROCESS ENERGY GROUP,I5)') IGR
+            WRITE(6,20) 'FLX1',FLX1(:,IGR,IC)
+            WRITE(6,20) 'TOT1-P0',TOT1(:,IGR,1,IC)
+            IF(NL.GE.2) WRITE(6,20) 'TOT1-P1',TOT1(:,IGR,2,IC)
+            WRITE(6,20) 'SIGR',TOT1(:,IGR,1,IC)-SCAT1(:,IGR,IGR,1,IC)
+            WRITE(6,20) 'DC1',DC1(:,IGR,IC)
+            WRITE(6,20) 'CHI1',CHI1(:,IGR,IC)
+            WRITE(6,20) 'SIGF1',SIGF1(:,IGR,IC)
+            DO JGR=1,NG
+              IF(IGR.EQ.JGR) THEN
+                WRITE(6,20) 'INSCAT1-P0',SCAT1(:,IGR,IGR,1,IC)
+                IF(NL.EQ.2) THEN
+                  WRITE(6,20) 'INSCAT1-P1',SCAT1(:,IGR,IGR,2,IC)
+                ENDIF
+              ELSE
+                WRITE(6,20) 'OUTSCAT1-P0',SCAT1(:,JGR,IGR,1,IC)
+                IF(NL.GE.2) THEN
+                  WRITE(6,20) 'OUTSCAT1-P1',SCAT1(:,JGR,IGR,2,IC)
+                ENDIF
+              ENDIF
+            ENDDO
+            WRITE(6,20) 'JXM',JXM(:,IGR,IC)
+            WRITE(6,20) 'JXP',JXP(:,IGR,IC)
+            DO IL=1,NL
+              WRITE(TEXT12,'(8HFHETXM-P,I1)') IL-1
+              WRITE(6,20) TEXT12(:9),FHETXM(:,IGR,IL,IC)
+              WRITE(TEXT12,'(8HFHETXP-P,I1)') IL-1
+              WRITE(6,20) TEXT12(:9),FHETXP(:,IGR,IL,IC)
+            ENDDO
+            IF(IDF.EQ.3) WRITE(6,20) 'ADF1',ADF1(:,IGR,IC)
+          ENDDO
+        ENDIF
+      ENDDO
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(WORK,IPOS,NJJ,IJJ,DCOU,SFIS1,SFIS,IMIX2,SCAT,SIGF,CHI,
+     1 DC,TOT,FLX,VOL)
+      RETURN
+*
+   10 FORMAT(1X,A12,10I13/(12X,10I13))
+   20 FORMAT(1X,A12,1P,10E13.5/(12X,10E13.5))
+      END