Initial commit from Polytechnique Montreal

1 files changed, 185 insertions, 0 deletions

diff --git a/Dragon/src/LIBND3.f b/Dragon/src/LIBND3.f
new file mode 100644
index 0000000..c2d2b62
--- /dev/null
+++ b/Dragon/src/LIBND3.f
@@ -0,0 +1,185 @@
+*DECK LIBND3
+      SUBROUTINE LIBND3(NGF,NGFR,NGRO,NBDIL,SN,SB,DILUS,DELTA,NF,XA,XS,
+     > XF,XN,GAR1,SCAT,GAR2,WT0)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Interpolate dilution-tabulated information, perform 
+* Livolant-Jeanpierre normalization and compute self-shielded 
+* cross sections at specific dilution.
+*
+*Copyright:
+* Copyright (C) 2006 Ecole Polytechnique de Montreal
+*
+*Author(s): A. Hebert
+*
+*Parameters: input
+* NGF     number of fast groups without self-shielding.
+* NGFR    number of fast and resonance groups.
+* NGRO    number of energy groups.
+* NBDIL   number of dilutions.
+* SN      actual dilution of the nuclide.
+* SB      dilution of the nuclide, as used in Livolant-Jeanpierre
+*         normalization.
+* DILUS   tabulation points in dilution.
+* DELTA   lethargy widths.
+* NF      flag set to 3 if fission information is present.
+* XA      tabulated absorption effective reaction rates.
+* XS      tabulated scattering effective reaction rates.
+* XF      tabulated nu*fission effective reaction rates.
+* XN      tabulated NJOY flux.
+* GAR1    infinite dilution cross sections.
+* SCAT    infinite dilution P0 differential scattering cross sections.
+*
+*Parameters: output
+* GAR2    interpolated self-shielded cross sections.
+* WT0     NJOY flux.
+*
+*Reference:
+* Copyright (C) from NDAS Atomic Energy of Canada Limited utility (2006)
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT NONE
+*----
+*  Subroutine arguments
+*----
+      INTEGER NGF,NGFR,NGRO,NBDIL,NF
+      REAL SN(NGRO),SB(NGRO),DILUS(NBDIL),DELTA(NGRO),
+     1 XA(NGFR-NGF,NBDIL),XS(NGFR-NGF,NBDIL),XF(NGFR-NGF,NBDIL),
+     2 XN(NGFR-NGF,NBDIL),GAR1(NGRO,6),SCAT(NGRO,NGRO),GAR2(NGRO,6),
+     3 WT0(NGRO)
+*----
+*  Local variables
+*----
+      REAL WW,ZNGAR,SSFACT,AUX,XN3
+      INTEGER I,IG,IG2
+      CHARACTER HSMG*131
+      LOGICAL LCUBIC
+      PARAMETER(LCUBIC=.TRUE.)
+      REAL, ALLOCATABLE, DIMENSION(:) :: TERPD,DD,XA2,XS2,XF2,WK
+      LOGICAL, ALLOCATABLE, DIMENSION(:) :: LINF
+*----
+*  Scratch storage allocation
+*----
+      ALLOCATE(TERPD(NBDIL),LINF(NGRO))
+*----
+*  Dilution interpolation
+*----
+      ALLOCATE(DD(NBDIL))
+      DO I=1,NBDIL
+        DD(I)=LOG10(DILUS(I))
+      ENDDO
+      ALLOCATE(XA2(NGRO),XS2(NGRO),XF2(NGRO))
+      DO IG=1,NGF
+        XA2(IG)=GAR1(IG,2)
+        XS2(IG)=GAR1(IG,5)
+        XF2(IG)=GAR1(IG,4)
+        WT0(IG)=1.0
+      ENDDO
+      ALLOCATE(WK(3*NBDIL))
+      DO IG=NGF+1,NGFR
+        TERPD(:NBDIL)=0.0
+        LINF(IG)=SN(IG).EQ.DILUS(NBDIL)
+        DO I=1,NBDIL
+          IF(ABS(SN(IG)-DILUS(I)).LE.1.0E-5*ABS(SN(IG))) THEN
+            TERPD(I)=1.0
+            GO TO 10
+          ENDIF
+        ENDDO
+        IF((NBDIL.EQ.1).OR.(SN(IG).GE.DILUS(NBDIL))) THEN
+*         No interpolation above infinite dilution
+          TERPD(NBDIL)=1.0
+        ELSE IF((NBDIL.EQ.2).OR.(SN(IG).GE.DILUS(NBDIL-1))) THEN
+*         One over SN interpolation near infinite dilution
+          LINF(IG)=.TRUE.
+          TERPD(NBDIL-1)=DILUS(NBDIL-1)/SN(IG)
+          TERPD(NBDIL)=1.0-DILUS(NBDIL-1)/SN(IG)
+        ELSE
+*         Perform Ceschino cubic interpolation
+          CALL ALTERP(LCUBIC,NBDIL-1,DD,LOG10(SN(IG)),.FALSE.,
+     >    TERPD,WK)
+        ENDIF
+   10   XA2(IG)=0.0
+        XS2(IG)=0.0
+        XF2(IG)=0.0
+        WT0(IG)=0.0
+        DO I=1,NBDIL
+          WW=TERPD(I)
+          IF(ABS(WW).GT.1.0E-6) THEN
+            XA2(IG)=XA2(IG)+WW*XA(IG-NGF,I)
+            XS2(IG)=XS2(IG)+WW*XS(IG-NGF,I)
+            IF(NF.EQ.3) XF2(IG)=XF2(IG)+WW*XF(IG-NGF,I)
+            WT0(IG)=WT0(IG)+WW*XN(IG-NGF,I)
+          ENDIF
+        ENDDO
+      ENDDO
+      DEALLOCATE(WK,DD)
+*----
+*  Livolant-Jeanpierre normalization
+*----
+      DO IG=NGF+1,NGFR
+        IF(SB(IG).NE.SN(IG)) THEN
+          ZNGAR=-(XA2(IG)+XS2(IG))
+          DO IG2=1,IG
+            SSFACT=XS2(IG2)/GAR1(IG2,5)
+            ZNGAR=ZNGAR+SCAT(IG,IG2)*SSFACT*DELTA(IG2)/DELTA(IG)
+          ENDDO
+          XN3=(WT0(IG)-1.0)*SN(IG)
+          IF(LINF(IG)) THEN
+*           Use an interpolated value near infinite dilution
+            AUX=(DILUS(NBDIL-1)/SB(IG))**2
+            XN3=AUX*XN3+(1.0-AUX)*ZNGAR
+            XN3=1.0+XN3/SB(IG)
+          ELSE
+            XN3=1.0+XN3/SB(IG)
+          ENDIF
+          IF((XN3.LE.0.0).OR.(XN3.GT.2.0)) THEN
+            WRITE (HSMG,100) XN3,IG,SB(IG),SN(IG)
+            CALL XABORT(HSMG)
+          ELSE IF(XN3.GT.1.2) THEN
+            WRITE (HSMG,100) XN3,IG,SB(IG),SN(IG)
+            WRITE(6,'(1X,A)') HSMG
+          ENDIF
+          WT0(IG)=XN3
+        ENDIF
+      ENDDO
+*----
+*  Divide effective reaction rates by NJOY flux for obtaining
+*  self-shielded cross sections
+*----
+      DO I=1,6
+        DO IG=1,NGRO
+          GAR2(IG,I)=GAR1(IG,2)
+        ENDDO
+      ENDDO
+      DO IG=NGF+1,NGFR
+*       Absorption xs
+        GAR2(IG,2)=XA2(IG)/WT0(IG)
+*       P0 scattering xs
+        GAR2(IG,5)=XS2(IG)/WT0(IG)
+*       nu*fission xs
+        GAR2(IG,4)=XF2(IG)/WT0(IG)
+*       Transport-corrected total xs
+        GAR2(IG,1)=GAR1(IG,1)*(GAR2(IG,2)+GAR2(IG,5))/(GAR1(IG,2)+
+     >  GAR1(IG,5))
+*       Fission xs
+        IF((NF.EQ.3).AND.(XF2(IG).EQ.0.0)) THEN
+          GAR2(IG,3)=GAR1(IG,3)
+        ELSE IF(NF.EQ.3) THEN
+          GAR2(IG,3)=GAR1(IG,3)*GAR2(IG,4)/GAR1(IG,4)
+        ENDIF
+*       P1 scattering xs
+        GAR2(IG,6)=GAR1(IG,6)*GAR2(IG,5)/GAR1(IG,5)
+      ENDDO
+      DEALLOCATE(XF2,XS2,XA2)
+*----
+*  Scratch storage deallocation
+*----
+      DEALLOCATE(LINF,TERPD)
+      RETURN
+*
+  100 FORMAT(37H LIBND3: Invalid value of NJOY flux (,1P,E11.3,
+     1 10H) in group,I4,11H. Dilution=,E11.3,2H (,E11.3,2H).)
+      END