summaryrefslogtreecommitdiff
path: root/Dragon/src/LIBTE2.f
blob: 77a10b74ef5e83913aee4a9763f198125239b606 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
*DECK LIBTE2
      SUBROUTINE LIBTE2 (NGRO,NSUBM,TMIX,SMIX,TN,SN,TERP)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute the Lagrange interpolation factors (TERP) for temperature and
* dilution interpolation of cross sections. TRANSX 2.0 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)
      REAL, ALLOCATABLE, DIMENSION(:) :: WORK
*----
*  LOCAL VARIABLES
*----
      PARAMETER (IOUT=6)
      DOUBLE PRECISION CHECK
*
      IF(NSUBM.EQ.1) THEN
         DO 10 I=1,NGRO
         TERP(I,1)=1.0
   10    CONTINUE
         RETURN
      ENDIF
*----
*  SCRATCH STORAGE ALLOCATION
*----
      ALLOCATE(WORK(NSUBM))
*
      DO 115 ISUBM=1,NSUBM
      TMAT=TMIX(ISUBM)
      SMAT=SMIX(ISUBM)
*----
*  FIND TEMPERATURES VALUES
*----
      NTEMP=0
      DO 20 JSUBM=1,NSUBM
      IF(SMIX(JSUBM).LT.0.9E10) GO TO 20
      NTEMP=NTEMP+1
      WORK(NTEMP)=TMIX(JSUBM)
   20 CONTINUE
      TERPT=0.0
*----
*  LIMIT TEMPERATURE INTERPOLATION TO LINEAR IS TN A GRID TEMPERATURE?
*----
      DO 30 ITMP=1,NTEMP
      TT=WORK(ITMP)
      IF(ABS(TN-TT).LT.1.E-3*TN+1.E-3) THEN
         IF(ABS(TN-TMAT).LT.1.E-3*TMAT+1.E-3) TERPT=1.0
         GO TO 70
      ENDIF
   30 CONTINUE
*----
*  IF TEMP OUT OF RANGE USE ENDPOINTS
*----
      IF((TN.GT.WORK(NTEMP).AND.ABS(TMAT-WORK(NTEMP)).LT.1.E-3*TMAT
     > +1.E-3).OR.(TN.LT.WORK(1).AND.ABS(TMAT-WORK(1)).LT.
     > 1.E-3*TMAT+1.E-3)) THEN
         TERPT=1.0
         GO TO 70
      ENDIF
*----
*  FIND BRACKETING TEMPS
*----
      IF(NTEMP.EQ.1) THEN
         TERPT=1.0
         GO TO 70
      ENDIF
      DO 40 ITMP=1,NTEMP-1
      IF(WORK(ITMP).LT.TN.AND.WORK(ITMP+1).GT.TN) THEN
         ILOW=ITMP
         IHIGH=ITMP+1
         IF(ABS(TMAT-WORK(ITMP)).LT.1.E-3*TMAT+1.E-3.OR.ABS(TMAT-
     >    WORK(ITMP+1)).LT.1.E-3*TMAT+1.E-3) GO TO 50
      ENDIF
   40 CONTINUE
      GO TO 70
   50 TERPT=1.0
      DO 60 ITMP=ILOW,IHIGH
      TT=WORK(ITMP)
      IF(ABS(TMAT-TT).LT.1.E-3*TMAT+1.E-3) GO TO 60
      TERPT=TERPT*(TN-TT)/(TMAT-TT)
      IF(ABS(TERPT).LT.1.E-3) GO TO 70
   60 CONTINUE
*
*----
*  FIND SIGMA-ZERO VALUES
*----
   70 NTEMP=0
      NSIGZ=0
      DO 80 JSUBM=1,NSUBM
      IF(SMIX(JSUBM).GE.0.9E10) NTEMP=NTEMP+1
      IF(NTEMP.GT.1) GO TO 80
      NSIGZ=NSIGZ+1
      WORK(NSIGZ)=SMIX(JSUBM)
   80 CONTINUE
*----
*  FIND TERP FACTOR FOR SIGMA-ZERO
*----
      DO 110 JJ=1,NGRO
      TERPS=0.0
      IF((SN(JJ).GE.WORK(1)).OR.(NSIGZ.EQ.1)) THEN
         IF(SMAT.EQ.WORK(1)) TERPS=1.0
      ELSE IF(SN(JJ).LE.WORK(NSIGZ)) THEN
         IF(SMAT.EQ.WORK(NSIGZ)) TERPS=1.0
      ELSE IF((SN(JJ).GT.WORK(2)).OR.(NSIGZ.EQ.2)) THEN
         IF(SMAT.EQ.WORK(1)) TERPS=1.0-WORK(2)/SN(JJ)
         IF(SMAT.EQ.WORK(2)) TERPS=WORK(2)/SN(JJ)
      ELSE
         DO 90 I=2,NSIGZ-1
         IF(SN(JJ).LT.WORK(I+1)) GO TO 90
         IF(SMAT.EQ.WORK(I+1)) TERPS=(ALOG10(WORK(I))-ALOG10(SN(JJ)))
     X    /(ALOG10(WORK(I))-ALOG10(WORK(I+1)))
         IF(SMAT.EQ.WORK(I)) TERPS=(ALOG10(SN(JJ))-ALOG10(WORK(I+1)))
     X    /(ALOG10(WORK(I))-ALOG10(WORK(I+1)))
         GO TO 100
   90    CONTINUE
      ENDIF
  100 TERP(JJ,ISUBM)=TERPT*TERPS
  110 CONTINUE
  115 CONTINUE
*----
*  CHECK FOR CONSISTENCY OF THE TERP FACTORS.
*----
      DO 130 JJ=1,NGRO
      CHECK=0.0D0
      DO 120 ISUBM=1,NSUBM
      CHECK=CHECK+TERP(JJ,ISUBM)
  120 CONTINUE
      IF(ABS(CHECK-1.0D0).GT.5.0D-3) THEN
         WRITE (IOUT,200) JJ,CHECK,(TERP(JJ,ISUBM),ISUBM=1,NSUBM)
         CALL XABORT('LIBTE2: INTERPOLATION FAILURE.')
      ENDIF
  130 CONTINUE
*----
*  SCRATCH STORAGE DEALLOCATION
*----
      DEALLOCATE(WORK)
      RETURN
*
  200 FORMAT (/51H LIBTE2: 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