summaryrefslogtreecommitdiff
path: root/Donjon/src/SCRISO.f
blob: dc381e85837b906f84ee9ce8309896569446edb8 (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
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
*DECK SCRISO
      SUBROUTINE SCRISO(IPLIB,NREA,NGRP,NL,NPRC,NOMREA,NWT0,XS,SIGS,
     > SS2D,TAUXFI,LXS,LAMB,CHIRS,BETAR,INVELS,INAME,LSTRD,LPURE,ILUPS,
     > ITRANC,IFISS)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Store an isotopic data recovered from a Saphyb into a Microlib.
*
*Copyright:
* Copyright (C) 2012 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
* IPLIB   address of the output microlib LCM object
* NREA    number of reactions in the Saphyb object
* NGRP    number of energy groups
* NL      maximum Legendre order (NL=1 is for isotropic scattering)
* NPRC    number of delayed neutron precursor groups
* NOMREA  names of reactions in the Saphyb
* NWT0    average flux
* XS      cross sections per reaction
* SIGS    scattering cross sections
* SS2D    complete scattering matrix
* TAUXFI  interpolated fission rate
* LXS     existence flag of each reaction
* LAMB    decay constants of the delayed neutron precursor groups
* CHIRS   delayed neutron emission spectrums
* BETAR   delayed neutron fractions
* INVELS  group-average of the inverse neutron velocity
* INAME   name of the isotope.
* LSTRD   flag set to .true. if B2=0.0.
* LPURE   =.true. if the interpolation is a pure linear interpolation 
*         with TERP factors.
* ILUPS   up-scattering removing flag (=1 to remove up-scattering from
*         output cross-sections).
*
*Parameters: output
* ITRANC  transport correction flag
* IFISS   fission flag
*
*-----------------------------------------------------------------------
*
      USE GANLIB
      IMPLICIT NONE
*----
*  SUBROUTINE ARGUMENTS
*----
      TYPE(C_PTR) IPLIB
      INTEGER NREA,NGRP,NL,NPRC,INAME(2),ILUPS,ITRANC,IFISS
      REAL NWT0(NGRP),XS(NGRP,NREA),SIGS(NGRP,NL),SS2D(NGRP,NGRP,NL),
     > TAUXFI,LAMB(NPRC),CHIRS(NGRP,NPRC),BETAR(NPRC),INVELS(NGRP)
      LOGICAL LXS(NREA),LSTRD,LPURE
      CHARACTER NOMREA(NREA)*12
*----
*  LOCAL VARIABLES
*----
      INTEGER I0, IGFROM, IGMAX, IGMIN, IGR, JGR, IGTO, ILEG, IPRC,
     & IREA, NXSCMP, IL, IRENT0,IRENT1
      REAL FF,CSCAT
      LOGICAL LDIFF,LZERO
      CHARACTER TEXT12*12
      CHARACTER HCM(0:10)*2,NAMLEG*2
      INTEGER, ALLOCATABLE, DIMENSION(:) :: ITYPRO,NJJ,IJJ
      REAL, ALLOCATABLE, DIMENSION(:) :: STRD,WRK,XSSCMP
      DATA HCM  /'00','01','02','03','04','05','06','07','08','09','10'/
*----
*  SCRATCH STORAGE ALLOCATION
*----
      ALLOCATE(STRD(NGRP))
*----
*  UP-SCATTERING CORRECTION
*----
      IF(ILUPS.EQ.1) THEN
        IRENT0=0
        IRENT1=0
        DO IREA=1,NREA
          IF(NOMREA(IREA).EQ.'TOTALE') IRENT0=IREA
          IF(NOMREA(IREA).EQ.'TOTALE P1') IRENT1=IREA
        ENDDO
        IF(IRENT0.EQ.0) CALL XABORT('SCRISO: MISSING NTOT0.')
        DO JGR=2,NGRP
          DO IGR=1,JGR-1 ! IGR < JGR
            CSCAT=SS2D(IGR,JGR,1)
            FF=NWT0(JGR)/NWT0(IGR)
            XS(IGR,IRENT0)=XS(IGR,IRENT0)-CSCAT*FF
            XS(JGR,IRENT0)=XS(JGR,IRENT0)-CSCAT
            IF((IRENT1.GT.0).AND.(NL.GT.1)) THEN
              CSCAT=SS2D(IGR,JGR,2)
              XS(IGR,IRENT1)=XS(IGR,IRENT1)-CSCAT*FF
              XS(JGR,IRENT1)=XS(JGR,IRENT1)-CSCAT
            ENDIF
            DO IL=1,NL
              CSCAT=SS2D(IGR,JGR,IL)
              SIGS(IGR,IL)=SIGS(IGR,IL)-CSCAT*FF
              SIGS(JGR,IL)=SIGS(JGR,IL)-CSCAT
              SS2D(JGR,IGR,IL)=SS2D(JGR,IGR,IL)-CSCAT*FF
              SS2D(IGR,JGR,IL)=0.0
            ENDDO
          ENDDO
        ENDDO
      ENDIF
*----
*  BUILD MICROLIB
*----
      WRITE(TEXT12,'(2A4)') (INAME(I0),I0=1,2)
      CALL LCMPTC(IPLIB,'ALIAS',12,TEXT12)
      CALL LCMPUT(IPLIB,'NWT0',NGRP,2,NWT0)
      IF(NPRC.GT.0) THEN
        CALL LCMPUT(IPLIB,'LAMBDA-D',NPRC,2,LAMB)
        CALL LCMPUT(IPLIB,'OVERV',NGRP,2,INVELS)
      ENDIF
      ITRANC=0
      IFISS=0
      LDIFF=.FALSE.
      STRD(:NGRP)=0.0
      DO IREA=1,NREA
        IF(.NOT.LXS(IREA)) CYCLE
        LZERO=.TRUE.
        DO IGR=1,NGRP
          LZERO=LZERO.AND.(XS(IGR,IREA).EQ.0.0)
        ENDDO
        IF(LZERO) CYCLE
        IF(NOMREA(IREA).EQ.'TOTALE') THEN
          IF(LSTRD) THEN
            DO IGR=1,NGRP
              STRD(IGR)=STRD(IGR)+XS(IGR,IREA)
            ENDDO
          ENDIF
          CALL LCMPUT(IPLIB,'NTOT0',NGRP,2,XS(1,IREA))
        ELSE IF(NOMREA(IREA).EQ.'TOTALE P1') THEN
          CALL LCMPUT(IPLIB,'NTOT1',NGRP,2,XS(1,IREA))
        ELSE IF(NOMREA(IREA).EQ.'EXCESS') THEN
*         correct scattering XS with excess XS
          DO IGR=1,NGRP
            SIGS(IGR,1)=SIGS(IGR,1)+XS(IGR,IREA)
          ENDDO
          CALL LCMPUT(IPLIB,'N2N',NGRP,2,XS(1,IREA))
        ELSE IF(NOMREA(IREA).EQ.'FISSION') THEN
          CALL LCMPUT(IPLIB,'NFTOT',NGRP,2,XS(1,IREA))
        ELSE IF(NOMREA(IREA).EQ.'ABSORPTION') THEN
          CALL LCMPUT(IPLIB,'NG',NGRP,2,XS(1,IREA))
        ELSE IF(NOMREA(IREA).EQ.'SPECTRE') THEN
          IF(.NOT.LPURE) THEN
            DO IGR=1,NGRP
              IF(XS(IGR,IREA).NE.0.0) THEN
                XS(IGR,IREA)=XS(IGR,IREA)/TAUXFI
              ENDIF
            ENDDO
          ENDIF
          CALL LCMPUT(IPLIB,'CHI',NGRP,2,XS(1,IREA))
          DO IPRC=1,NPRC
            WRITE(TEXT12,'(A3,I2.2)') 'CHI',IPRC
            CALL LCMPUT(IPLIB,TEXT12,NGRP,2,CHIRS(1,IPRC))
          ENDDO
        ELSE IF(NOMREA(IREA).EQ.'NU*FISSION') THEN
          IFISS=1
          CALL LCMPUT(IPLIB,'NUSIGF',NGRP,2,XS(1,IREA))
          IF(NPRC.GT.0) THEN
            ALLOCATE(WRK(NGRP))
            DO IPRC=1,NPRC
              DO IGR=1,NGRP
                WRK(IGR)=XS(IGR,IREA)*BETAR(IPRC)
              ENDDO
              WRITE(TEXT12,'(A6,I2.2)') 'NUSIGF',IPRC
              CALL LCMPUT(IPLIB,TEXT12,NGRP,2,WRK)
            ENDDO
            DEALLOCATE(WRK)
          ENDIF
        ELSE IF(NOMREA(IREA).EQ.'ENERGIE') THEN
          ALLOCATE(WRK(NGRP))
          DO IGR=1,NGRP
            WRK(IGR)=XS(IGR,IREA)*1.0E6 ! convert MeV to eV
          ENDDO
          CALL LCMPUT(IPLIB,'H-FACTOR',NGRP,2,WRK)
          DEALLOCATE(WRK)
        ELSE IF(NOMREA(IREA).EQ.'SELF') THEN
          CALL LCMPUT(IPLIB,'SIGW00',NGRP,2,XS(1,IREA))
        ELSE IF(NOMREA(IREA).EQ.'TRANSP-CORR') THEN
          ITRANC=2
          IF(LSTRD) THEN
            DO IGR=1,NGRP
              STRD(IGR)=STRD(IGR)-XS(IGR,IREA)
            ENDDO
          ENDIF
          CALL LCMPUT(IPLIB,'TRANC',NGRP,2,XS(1,IREA))
        ELSE IF(NOMREA(IREA).EQ.'FUITES') THEN
          LDIFF=LSTRD
          IF(.NOT.LSTRD) THEN
            DO IGR=1,NGRP
              LDIFF=LDIFF.OR.(XS(IGR,IREA).NE.0.0)
              STRD(IGR)=XS(IGR,IREA)
            ENDDO
          ENDIF
        ELSE IF(NOMREA(IREA).EQ.'DIFFUSION') THEN
          CYCLE
        ELSE IF(NOMREA(IREA).EQ.'TRANSFERT') THEN
          CYCLE
        ELSE
          CALL LCMPUT(IPLIB,NOMREA(IREA),NGRP,2,XS(1,IREA))
        ENDIF
      ENDDO
      IF(LSTRD) THEN
        IF((ITRANC.EQ.0).AND.(NL.GT.1)) THEN
*         Apollo-type transport correction
          DO IGR=1,NGRP
            STRD(IGR)=STRD(IGR)-SIGS(IGR,2)
          ENDDO
        ENDIF
      ELSE
        DO IGR=1,NGRP
          STRD(IGR)=1.0/(3.0*STRD(IGR))
        ENDDO
      ENDIF
      IF((ITRANC.EQ.0).AND.(NL.GT.1)) THEN
*       Apollo-type transport correction
        ITRANC=2
        CALL LCMPUT(IPLIB,'TRANC',NGRP,2,SIGS(1,2))
      ENDIF
      IF(LDIFF.OR.LSTRD) CALL LCMPUT(IPLIB,'STRD',NGRP,2,STRD)
*----
*  SAVE SCATTERING VECTORS AND MATRICES (DO NOT USE XDRLGS TO SAVE CPU
*  TIME)
*----
      ALLOCATE(NJJ(NGRP),IJJ(NGRP),XSSCMP(NGRP*NGRP),ITYPRO(NL))
      DO ILEG=1,NL
        IF(ILEG.LE.11) THEN
          NAMLEG=HCM(ILEG-1)
        ELSE
          WRITE(NAMLEG,'(I2.2)') ILEG-1
        ENDIF
        CALL LCMPUT(IPLIB,'SIGS'//NAMLEG,NGRP,2,SIGS(1,ILEG))
        NXSCMP=0
        DO IGTO=1,NGRP
          IGMIN=IGTO
          IGMAX=IGTO
          DO IGFROM=1,NGRP
            IF(SS2D(IGTO,IGFROM,ILEG).NE.0.0) THEN
              IGMIN=MIN(IGMIN,IGFROM)
              IGMAX=MAX(IGMAX,IGFROM)
            ENDIF
          ENDDO
          IJJ(IGTO)=IGMAX
          NJJ(IGTO)=IGMAX-IGMIN+1
          DO IGFROM=IGMAX,IGMIN,-1
            NXSCMP=NXSCMP+1
            XSSCMP(NXSCMP)=SS2D(IGTO,IGFROM,ILEG)
          ENDDO
        ENDDO
        CALL LCMPUT(IPLIB,'NJJS'//NAMLEG,NGRP,1,NJJ)
        CALL LCMPUT(IPLIB,'IJJS'//NAMLEG,NGRP,1,IJJ)
        CALL LCMPUT(IPLIB,'SCAT'//NAMLEG,NXSCMP,2,XSSCMP)
        ITYPRO(ILEG)=1
      ENDDO
      CALL LCMPUT(IPLIB,'SCAT-SAVED',NL,1,ITYPRO)
      DEALLOCATE(ITYPRO,XSSCMP,IJJ,NJJ)
*----
*  SCRATCH STORAGE DEALLOCATION
*----
      DEALLOCATE(STRD)
      RETURN
      END