Initial commit from Polytechnique Montreal

1 files changed, 142 insertions, 0 deletions

diff --git a/Dragon/src/LIBMOM.f b/Dragon/src/LIBMOM.f
new file mode 100644
index 0000000..94f472b
--- /dev/null
+++ b/Dragon/src/LIBMOM.f
@@ -0,0 +1,142 @@
+*DECK LIBMOM
+      SUBROUTINE LIBMOM(NFS,NDIL,NPAR,DELTA,SIGTF,SIGSF,SIGFF,NOR,
+     1 SIGERD,MOMT,MOMP,SEFFER)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute a set of total (SIGT**N) and partial (SIGA*(SIGT**N)) moments
+* and a set of reference self-shielded flux and cross sections.
+*
+*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
+*
+*Parameters: input
+* NFS     number of fine energy groups.
+* NDIL    number of dilutions.
+* NPAR    number of partial cross sections (NPAR=0, 1 or 2).
+* DELTA   lethargy widths of the fine groups.
+* SIGTF   microscopic principal x-sections in the fine groups.
+* SIGSF   microscopic partial x-sections 1 in the fine groups.
+* SIGFF   microscopic partial x-sections 2 in the fine groups.
+* NOR     related to the number of moments to preserve.
+*         For the total moments: -NOR+1 <= N <= NOR.
+*         For the partial moments: -NOR/2 <= N <= (NOR-1)/2.
+* SIGERD  dilutions used to compute SEFFER.
+*
+*Parameters: output
+* MOMT    total moments.
+* MOMP    partial moments in absorption.
+* SEFFER  Bondarenko self-shielded flux and cross sections at each
+*         dilution.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER NFS,NDIL,NPAR,NOR
+      REAL DELTA(NFS),SIGTF(NFS),SIGSF(NFS),SIGFF(NFS),SIGERD(NDIL),
+     1 SEFFER(NPAR+2,NDIL)
+      DOUBLE PRECISION MOMT(2*NOR),MOMP(NOR,NPAR)
+*----
+*  LOCAL VARIABLES
+*----
+      DOUBLE PRECISION T,DEL
+*
+      DEL=0.0D0
+      DO 10 I=1,2*NOR
+      MOMT(I)=0.0D0
+   10 CONTINUE
+      DO 21 I=1,NOR
+      DO 20 J=1,NPAR
+      MOMP(I,J)=0.0D0
+   20 CONTINUE
+   21 CONTINUE
+      DO 26 IDIL=1,NDIL
+      DO 25 IPAR=1,NPAR+2
+      SEFFER(IPAR,IDIL)=0.0
+   25 CONTINUE
+   26 CONTINUE
+*
+      DO 70 IGF=1,NFS
+        DELF=DELTA(IGF)
+        SIGT=MAX(0.001,SIGTF(IGF))
+        DEL=DEL+DELF
+        T=DELF
+        DO 30 INOR=NOR,2*NOR
+        MOMT(INOR)=MOMT(INOR)+T
+        T=T*SIGT
+   30   CONTINUE
+        T=DELF/SIGT
+        DO 40 INOR=NOR-1,1,-1
+        MOMT(INOR)=MOMT(INOR)+T
+        T=T/SIGT
+   40   CONTINUE
+        DO 45 IDIL=1,NDIL
+        T=SIGERD(IDIL)*DELF
+        SEFFER(1,IDIL)=SEFFER(1,IDIL)+REAL(T)/(SIGERD(IDIL)+SIGT)
+        SEFFER(2,IDIL)=SEFFER(2,IDIL)+REAL(T)*SIGT/(SIGERD(IDIL)+SIGT)
+   45   CONTINUE
+        IF(NPAR.GT.0) THEN
+          SIGS=MAX(1.E-9,SIGSF(IGF))
+          T=DELF*SIGS
+          DO 50 INOR=NOR/2+1,NOR
+          MOMP(INOR,1)=MOMP(INOR,1)+T
+          T=T*SIGT
+   50     CONTINUE
+          T=DELF*SIGS/SIGT
+          DO 60 INOR=NOR/2,1,-1
+          MOMP(INOR,1)=MOMP(INOR,1)+T
+          T=T/SIGT
+   60     CONTINUE
+          DO 65 IDIL=1,NDIL
+          T=SIGERD(IDIL)*DELF
+          SEFFER(3,IDIL)=SEFFER(3,IDIL)+REAL(T)*SIGS/(SIGERD(IDIL)+SIGT)
+   65     CONTINUE
+          IF(NPAR.EQ.2) THEN
+            SIGF=MAX(1.E-9,SIGFF(IGF))
+            T=DELF*SIGF
+            DO 500 INOR=NOR/2+1,NOR
+            MOMP(INOR,2)=MOMP(INOR,2)+T
+            T=T*SIGT
+  500       CONTINUE
+            T=DELF*SIGF/SIGT
+            DO 600 INOR=NOR/2,1,-1
+            MOMP(INOR,2)=MOMP(INOR,2)+T
+            T=T/SIGT
+  600       CONTINUE
+            DO 650 IDIL=1,NDIL
+            T=SIGERD(IDIL)*DELF
+            SEFFER(4,IDIL)=SEFFER(4,IDIL)+REAL(T)*SIGF/(SIGERD(IDIL)+
+     1      SIGT)
+  650       CONTINUE
+          ENDIF
+        ENDIF
+   70 CONTINUE
+*
+      IF(DEL.EQ.0.0) CALL XABORT('LIBMOM: ALGORITHM FAILURE.')
+      DO 80 INOR=1,2*NOR
+      MOMT(INOR)=MOMT(INOR)/DEL
+   80 CONTINUE
+      IF(NPAR.GT.0) THEN
+         DO 95 INOR=1,NOR
+         DO 90 IPAR=1,NPAR
+         MOMP(INOR,IPAR)=MOMP(INOR,IPAR)/REAL(DEL)
+   90    CONTINUE
+   95    CONTINUE
+      ENDIF
+      DO 110 IDIL=1,NDIL
+      DO 100 IPAR=2,NPAR+2
+      SEFFER(IPAR,IDIL)=SEFFER(IPAR,IDIL)/SEFFER(1,IDIL)
+  100 CONTINUE
+      SEFFER(1,IDIL)=SEFFER(1,IDIL)/REAL(DEL)
+  110 CONTINUE
+      RETURN
+      END