Initial commit from Polytechnique Montreal

1 files changed, 124 insertions, 0 deletions

diff --git a/Dragon/src/LIBTER.f b/Dragon/src/LIBTER.f
new file mode 100644
index 0000000..8fe3641
--- /dev/null
+++ b/Dragon/src/LIBTER.f
@@ -0,0 +1,124 @@
+*DECK LIBTER
+      SUBROUTINE LIBTER (NGRO,NSUBM,TMIX,SMIX,TN,SN,TERP)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the Lagrange interpolation factors (TERP) for temperature and
+* dilution interpolation of cross sections. TRANSX CTR algorithm.
+*
+*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
+* NGRO    number of energy groups.
+* NSUBM   number of submaterials (number of temperature/dilution
+*         collocation points).
+* TMIX    temperature of each submaterial in the library.
+* SMIX    dilution of each submaterial in the library.
+*         The submaterials are ordered by decreasing dilution and
+*         then by increasing temperature.
+* TN      temperature of the isotope.
+* SN      dilution cross section in each energy group of the isotope.
+*         A value of 1.0E10 is used for infinite dilution.
+*
+*Parameters: output
+* TERP    Lagrange interpolation factors.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER NGRO,NSUBM
+      REAL TMIX(NSUBM),SMIX(NSUBM),TN,SN(NGRO),TERP(NGRO,NSUBM)
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER (IOUT=6)
+      DOUBLE PRECISION CHECK
+*----
+*  TEMPERATURE AND BACKGROUND LOOP
+*----
+      IF (NSUBM.EQ.1) THEN
+         DO 10 I=1,NGRO
+         TERP(I,1)=1.0
+   10    CONTINUE
+         RETURN
+      ENDIF
+      BREAK=0.0
+      SMIN=1.0E10
+      DO 20 ISUBM=1,NSUBM
+      ST=SMIX(ISUBM)
+      IF ((ST.LT.1.0E8).AND.(ST.GT.BREAK)) BREAK=ST
+      IF (ST.LT.SMIN) SMIN=ST
+   20 CONTINUE
+      DO 70 ISUBM=1,NSUBM
+      DO 30 JJ=1,NGRO
+      TERP(JJ,ISUBM)=0.0
+   30 CONTINUE
+      TMAT=TMIX(ISUBM)
+      SMAT=SMIX(ISUBM)
+*----
+*  COMPUTE TERP FACTORS
+*----
+      TERPT=1.0
+      DO 40 ISM=1,NSUBM
+      TT=TMIX(ISM)
+      ST=SMIX(ISM)
+      IF (ST.LE.0.99E10) GO TO 40
+      IF (ABS(TMAT-TT).LT.1.0E-5*TMAT+1.0E-5) GO TO 40
+      TERPT=TERPT*(TN-TT)/(TMAT-TT)
+      IF (ABS(TERPT).LT.1.0E-5) GO TO 70
+   40 CONTINUE
+      DO 60 JJ=1,NGRO
+      TERPS=TERPT
+      IF ((SN(JJ).LT.SMIN).AND.(SMAT.GT.1.01*SMIN)) GO TO 60
+      IF ((SN(JJ).LT.BREAK).AND.(SMAT.GT.1.01*BREAK)) GO TO 60
+      IF ((SN(JJ).GE.BREAK).AND.(SMAT.LT.0.99*BREAK)) GO TO 60
+      TLAST=-1.0
+      DO 50 ISM=1,NSUBM
+      TT=TMIX(ISM)
+      ST=SMIX(ISM)
+      IF (TLAST.LT.0.) TLAST=TT
+      IF (TT.NE.TLAST) GO TO 50
+      IF (ABS(SMAT-ST).LT.1.0E-5*SMAT) GO TO 50
+      IF ((SN(JJ).GE.SMIN).AND.(SN(JJ).LT.BREAK)) THEN
+         IF (ST.GT.1.01*BREAK) GO TO 50
+         TERPS=TERPS*LOG(SN(JJ)/ST)/LOG(SMAT/ST)
+      ELSE IF ((SN(JJ).GE.SMIN).AND.(SN(JJ).GE.BREAK)) THEN
+         IF (ST.LT.0.99*BREAK) GO TO 50
+         TERPS=TERPS*((ST/SN(JJ))-1.)/((ST/SMAT)-1.)
+      ELSE
+         IF (ST.GT.1.01*SMIN) GO TO 50
+         TERPS=TERPS*(SN(JJ)**2-ST**2)/(SMAT**2-ST**2)
+      ENDIF
+      IF (ABS(TERPS).LE.1.0E-5) GO TO 60
+   50 CONTINUE
+      TERP(JJ,ISUBM)=TERPS
+   60 CONTINUE
+   70 CONTINUE
+*----
+*  CHECK FOR CONSISTENCY OF THE TERP FACTORS.
+*----
+      DO 90 JJ=1,NGRO
+      CHECK=0.0D0
+      DO 80 ISUBM=1,NSUBM
+      CHECK=CHECK+TERP(JJ,ISUBM)
+   80 CONTINUE
+      IF (ABS(CHECK-1.0D0).GT.5.0D-3) THEN
+         WRITE (IOUT,100) JJ,CHECK,(TERP(JJ,ISUBM),ISUBM=1,NSUBM)
+         CALL XABORT('LIBTER: INTERPOLATION FAILURE.')
+      ENDIF
+   90 CONTINUE
+      RETURN
+  100 FORMAT (/51H LIBTER: INCONSISTENT LAGRANGE INTERPOLATION FACTOR,
+     1 9H IN GROUP,I4,8H. CHECK=,1P,E13.3/9H FACTORS=,1P,9E13.3/
+     2 (9X,1P,9E13.3))
+      END