summaryrefslogtreecommitdiff
path: root/Dragon/src/ASMDRV.f
blob: 27e7296962d2c888a170ed8664aea52b02269dd1 (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
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
*DECK ASMDRV
      SUBROUTINE ASMDRV(IPSYS,IPTRK,IPMACR,IFTRAK,CDOOR,IPRNTP,NGROUP,
     >                  NBMIX,NREGIO,NANI,NANIST,NW,MATCOD,VOLUME,
     >                  LEAKSW,ITRANC,LDIFF,IBFP,TITRE,ITPIJ,LNORM,
     >                  IPHASE,ISTRM,KNORM,NALBP)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Dragon assembly and pij phases.
*
*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): G. Marleau
*
*Parameters: input
* IPSYS   pointer to the pij LCM object.
* IPTRK   pointer to the tracking LCM object.
* IPMACR  pointer to the macrolib LCM object.
* IFTRAK  file unit number for tracks.
* CDOOR   name of the pij calculation door.
* IPRNTP  print option for pij calculations.
* NGROUP  number of groups treated.
* NBMIX   number of mixtures considered.
* NREGIO  number of regions considered .
* NANI    number of Legendre orders for scattering cross sections.
* NANIST  number of Legendre orders for scattering cross sections
*         if streaming leakage is present.
* NW      type of weighting for P1 cross section info (=0: P0; =1: P1).
* MATCOD  mixture code in each region.
* VOLUME  volume of each region.
* LEAKSW  leakage switch.
* ITRANC  type of transport correction.
* LDIFF   diffusion coefficient switch.
* IBFP    Fokker-Planck solution (=0: off; =1/2: on).
* TITRE   execution title.
* ITPIJ   type of collision probability available:
*         =1 scatt mod pij (wij);
*         =2 stand. pij;
*         =3 scatt mod pij+pijk (wij,wijk);
*         =4 stand. pij+pijk.
* LNORM   switch for removing leakage from collision probabilities and
*         keeping the pis information.
* IPHASE  type of assembly (=1 for ass and 2 for pij).
* ISTRM   type of streaming effect:
*         =1 no streaming effect;
*         =2 isotropic streaming effect;
*         =3 anisotropic streaming effect.
* KNORM   type of pij normalization.
* NALBP   number of physical albedos.
*
*-----------------------------------------------------------------------
*
      USE GANLIB
*----
*  SUBROUTINE ARGUMENTS
*----
      CHARACTER   CDOOR*12,TITRE*72,TEXT12*12
      LOGICAL     LEAKSW,LDIFF,LNORM
      TYPE(C_PTR) IPSYS,IPTRK,IPMACR
      INTEGER     IFTRAK,IPRNTP,NGROUP,NBMIX,NREGIO,NANI,NANIST,NW,
     >            MATCOD(NREGIO),ITRANC,IBFP,ITPIJ,IPHASE,ISTRM,KNORM,
     >            NALBP
      REAL        VOLUME(NREGIO)
*----
*  LOCAL VARIABLES
*----
      PARAMETER  (IUNOUT=6,ILCMUP=1,ILCMDN=2)
      LOGICAL     LTRANC
      CHARACTER   HSMG*130,CM*2
      INTEGER, ALLOCATABLE, DIMENSION(:) :: NPSYS
      TYPE(C_PTR) JPSYS,KPSYS,JPMACR,KPMACR
      REAL, ALLOCATABLE, DIMENSION(:) :: TEMP,ENERGY
      REAL, ALLOCATABLE, DIMENSION(:,:) :: XSSCOR,XSDIFF,ALBP,EMOMTR
      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: XSSIGT,XSSIGW,ESTOPW
*----
*  SCRATCH STORAGE ALLOCATION
*----
      ALLOCATE(NPSYS(NGROUP))
      ALLOCATE(XSSIGT(0:NBMIX,NGROUP,NW+1),XSSCOR(NBMIX,NGROUP),
     >         XSSIGW(0:NBMIX,NANIST,NGROUP),XSDIFF(0:NBMIX,NGROUP),
     >         ESTOPW(0:NBMIX,2,NGROUP),ENERGY(NGROUP+1),
     >         EMOMTR(0:NBMIX,NGROUP))
      ALLOCATE(ALBP(NALBP,NGROUP))
*----
*  CHECK MIXTURE INDICES.
*----
      DO 10 I=1,NREGIO
      IF(MATCOD(I).GT.NBMIX) THEN
         WRITE (HSMG,5090) NBMIX
         CALL XABORT(HSMG)
      ENDIF
   10 CONTINUE
*----
*  RECOVER PHYSICAL ALBEDOS.
*----
      IF(NALBP.GT.0) THEN
        CALL LCMLEN(IPMACR,'ALBEDO',ILCMLN,ITYLCM)
        IF(ILCMLN.EQ.NALBP*NGROUP) THEN
          CALL LCMGET(IPMACR,'ALBEDO',ALBP)
        ELSE
          CALL LCMLIB(IPMACR)
          CALL XABORT('ASMDRV: READ ERROR ON LCM RECORD= ALBEDO')
        ENDIF
      ENDIF
*----
*  RECOVER ENERGY MESH VALUES.
*----
      IF(IBFP.GT.0) CALL LCMGET(IPMACR,'ENERGY',ENERGY)
*----
*  READ X-SECTIONS AND COMPUTE TRANSPORT CORRECTED X-SECTIONS.
*----
      IF(IPRNTP.GE.1) THEN
        IF(IPHASE.EQ.1) THEN
          WRITE(IUNOUT,6200) CDOOR
        ELSE
          WRITE(IUNOUT,6201) CDOOR
        ENDIF
        IF(ITRANC.NE.0) WRITE(IUNOUT,6101) ITRANC
      ENDIF
      IF(IPRNTP.GE.2) THEN
        WRITE(IUNOUT,6000)
        WRITE(IUNOUT,6001) (IREGIO,VOLUME(IREGIO),MATCOD(IREGIO),
     >                      IREGIO=1,NREGIO)
      ENDIF
      CALL LCMLEN(IPMACR,'GROUP',ILON,ITYLCM)
      IF(ILON.NE.NGROUP) CALL XABORT('ASMDRV: INVALID MACROLIB.')
      JPMACR=LCMGID(IPMACR,'GROUP')
      JPSYS=LCMLID(IPSYS,'GROUP',NGROUP)
      ITRAN2=0
      DO 60 IGR=1,NGROUP
        KPMACR=LCMGIL(JPMACR,IGR)
        DO 20 IW=1,MIN(NW+1,10)
        WRITE(TEXT12,'(4HNTOT,I1)') IW-1
        CALL LCMLEN(KPMACR,TEXT12,ILCMLN,ITYLCM)
        IF(ILCMLN.GT.0) THEN
          CALL LCMGET(KPMACR,TEXT12,XSSIGT(1,IGR,IW))
          XSSIGT(0,IGR,IW)=0.0
        ELSE IF(IW.EQ.1) THEN
          CALL LCMLIB(KPMACR)
          CALL XABORT('ASMDRV: READ ERROR ON LCM RECORD= TOTAL')
        ELSE
          CALL LCMGET(KPMACR,'NTOT0',XSSIGT(1,IGR,IW))
          XSSIGT(0,IGR,IW)=0.0
        ENDIF
   20   CONTINUE
        DO 30 IL=1,NANIST
        WRITE(CM,'(I2.2)') IL-1
        CALL LCMLEN(KPMACR,'SIGW'//CM,ILCMLN,ITYLCM)
        IF(ILCMLN.GT.0) THEN
          CALL LCMGET(KPMACR,'SIGW'//CM,XSSIGW(1,IL,IGR))
          XSSIGW(0,IL,IGR)=0.0
        ELSE IF(IL.EQ.1) THEN
          CALL LCMLIB(KPMACR)
          CALL XABORT('ASMDRV: READ ERROR ON LCM RECORD= SIGW'//CM)
        ELSE
          XSSIGW(0:NBMIX,IL,IGR)=0.0
        ENDIF
   30   CONTINUE
        CALL LCMLEN(KPMACR,'TRANC',ILCMLN,ITYLCM)
        LTRANC=ILCMLN.GT.0
        IF((ITRANC.NE.0).AND.LTRANC) THEN
*         TRANSPORT CORRECTION (INCLUDE THE LEAKAGE CORRECTION).
          ITRAN2=ITRANC
          CALL LCMGET(KPMACR,'TRANC',XSSCOR(1,IGR))
        ELSE
          ITRAN2=0
          XSSCOR(:NBMIX,IGR)=0.0
        ENDIF
        IF(ITRAN2.NE.0) THEN
*         INCLUDE TRANSPORT CORRECTION.
          DO 40 IMAT=1,NBMIX
            DELTA=XSSIGT(IMAT,IGR,1)-XSSIGW(IMAT,1,IGR)
            XSSIGT(IMAT,IGR,1)=XSSIGT(IMAT,IGR,1)-XSSCOR(IMAT,IGR)
            XSSIGW(IMAT,1,IGR)=XSSIGT(IMAT,IGR,1)-DELTA
*           Tibere is using transport-corrected XS for the second
*           equation. Scattering reduction must be performed with
*           transport-corrected SIGS1 values.
            IF(ISTRM.EQ.3) THEN
              XSSIGW(IMAT,2,IGR)=XSSIGW(IMAT,2,IGR)-XSSCOR(IMAT,IGR)
            ENDIF
   40     CONTINUE
        ENDIF
*
        IF(NW.GT.0) THEN
*         PERFORM A P0_TOTAL LEAKAGE CORRECTION.
          DO 55 IW=2,MIN(NANIST,NW+1)
          DO 50 IMAT=1,NBMIX
          DELTA=XSSIGT(IMAT,IGR,1)-XSSIGT(IMAT,IGR,IW)
          IF((ITRAN2.NE.0).AND.(DELTA.NE.0.0)) THEN
            CALL XABORT('ASMDRV: CANNOT PERFORM BOTH TRANSPORT AND LEA'
     >      //'KAGE CORRECTIONS.')
          ENDIF
          XSSIGT(IMAT,IGR,IW)=XSSIGT(IMAT,IGR,1)
          XSSIGW(IMAT,IW,IGR)=XSSIGW(IMAT,IW,IGR)+DELTA
   50     CONTINUE
   55     CONTINUE
        ENDIF
*
        IF(IPRNTP.GE.3) THEN
          WRITE(IUNOUT,6002)  IGR
          WRITE(IUNOUT,6003) (IMIX,XSSIGT(IMIX,IGR,1),XSSIGW(IMIX,1,IGR)
     >    ,IMIX=1,NBMIX)
        ENDIF
        IF(LDIFF) THEN
*         INCLUDE DIFFUSION COEFFICIENTS.
          CALL LCMGET(KPMACR,'DIFF',XSDIFF(1,IGR))
          XSDIFF(0,IGR)=1.0E10
        ENDIF
        IF(IBFP.GT.0) THEN
*         INCLUDE RESTRICTED STOPPING POWER.
          ALLOCATE(TEMP(2*NBMIX))
          CALL LCMGET(KPMACR,'ESTOPW',TEMP)
          ESTOPW(0,:2,IGR)=0.0
          ESTOPW(1:NBMIX,1,IGR)=TEMP(:NBMIX)
          ESTOPW(1:NBMIX,2,IGR)=TEMP(NBMIX+1:)
          DEALLOCATE(TEMP)
*         INCLUDE RESTRICTED MOMENTUM TRANSFER.
          EMOMTR(0,IGR)=0.0
          CALL LCMLEN(KPMACR,'EMOMTR',ILCMLN,ITYLCM)
          IF(ILCMLN.NE.0) THEN
            ALLOCATE(TEMP(NBMIX))
            CALL LCMGET(KPMACR,'EMOMTR',TEMP)
            EMOMTR(1:NBMIX,IGR)=TEMP(:NBMIX)
            DEALLOCATE(TEMP)
          ELSE
            EMOMTR(1:NBMIX,IGR)=0.0
          ENDIF
        ENDIF
   60 CONTINUE
*----
*  COMPUTE THE CP OR RESPONSE MATRIX INFORMATION FOR THE SOLUTION OF
*  THE BALANCE EQUATION.
*----
      IPIJK=1
      IF(ISTRM.EQ.3) IPIJK=4
      DO 70 IGR=1,NGROUP
      NPSYS(IGR)=IGR
      KPSYS=LCMDIL(JPSYS,IGR)
      CALL LCMPUT(KPSYS,'DRAGON-TXSC',NBMIX+1,2,XSSIGT(0,IGR,1))
      IF(NW.GT.0) THEN
         CALL LCMPUT(KPSYS,'DRAGON-T1XSC',NBMIX+1,2,XSSIGT(0,IGR,2))
      ENDIF
      CALL LCMPUT(KPSYS,'DRAGON-S0XSC',(NBMIX+1)*NANI,2,XSSIGW(0,1,IGR))
      IF(LDIFF) CALL LCMPUT(KPSYS,'DRAGON-DIFF',NBMIX+1,2,XSDIFF(0,IGR))
      IF(IBFP.GT.0) THEN
        CALL LCMPUT(KPSYS,'DRAGON-ESTOP',(NBMIX+1)*2,2,ESTOPW(0,1,IGR))
        CALL LCMPUT(KPSYS,'DRAGON-EMOMT',NBMIX+1,2,EMOMTR(0,IGR))
        DELTAE=(ENERGY(IGR)-ENERGY(IGR+1))/1.0E6
        CALL LCMPUT(KPSYS,'DRAGON-DELTE',1,2,DELTAE)
        IF(IGR.EQ.NGROUP) THEN
          CALL LCMPUT(KPSYS,'DRAGON-ISLG',1,1,1)
        ELSE
          CALL LCMPUT(KPSYS,'DRAGON-ISLG',1,1,0)
        ENDIF
      ENDIF
      IF(NALBP.GT.0) CALL LCMPUT(KPSYS,'ALBEDO',NALBP,2,ALBP(1,IGR))
   70 CONTINUE
      IF(IBFP.GT.0) THEN
        CALL LCMPUT(IPSYS,'ECUTOFF',1,2,ENERGY(NGROUP+1)/1.0E6)
      ENDIF
      IF(IPHASE.EQ.2) THEN
         CALL DOORPV(CDOOR,JPSYS,NPSYS,IPTRK,IFTRAK,IPRNTP,NGROUP,
     >   NREGIO,NBMIX,NANI,MATCOD,VOLUME,KNORM,IPIJK,LEAKSW,ITPIJ,
     >   LNORM,TITRE,NALBP)
      ELSE
         CALL DOORAV(CDOOR,JPSYS,NPSYS,IPTRK,IFTRAK,IPRNTP,NGROUP,
     >   NREGIO,NBMIX,NANI,NW,MATCOD,VOLUME,KNORM,LEAKSW,TITRE,NALBP,
     >   ISTRM)
      ENDIF
*----
*  COMPUTE THE P1 CP OR RESPONSE MATRIX INFORMATION FOR THE ECCO
*  ISOTROPIC STREAMING MODEL.
*----
      IF(ISTRM.EQ.2) THEN
         CALL LCMSIX(IPSYS,'STREAMING',ILCMUP)
         JPSYS=LCMLID(IPSYS,'GROUP',NGROUP)
         IF(ITRAN2.NE.0) THEN
*          REMOVE TRANSPORT CORRECTION.
           DO 85 IGR=1,NGROUP
           DO 80 IMAT=1,NBMIX
             XSSIGT(IMAT,IGR,1)=XSSIGT(IMAT,IGR,1)+XSSCOR(IMAT,IGR)
             XSSIGW(IMAT,1,IGR)=XSSIGW(IMAT,1,IGR)+XSSCOR(IMAT,IGR)
   80      CONTINUE
   85      CONTINUE
         ENDIF
         IF(NANIST.LE.1) CALL XABORT('ASMDRV: MISSING P1 XS INFO.')
         DO 90 IGR=1,NGROUP
         NPSYS(IGR)=IGR
         KPSYS=LCMDIL(JPSYS,IGR)
         CALL LCMPUT(KPSYS,'DRAGON-TXSC',NBMIX+1,2,XSSIGT(0,IGR,1))
         CALL LCMPUT(KPSYS,'DRAGON-S0XSC',(NBMIX+1)*(NANIST-1),2,
     >   XSSIGW(0,2,IGR))
         IF(LDIFF) CALL LCMPUT(KPSYS,'DRAGON-DIFF',NBMIX+1,2,
     >   XSDIFF(0,IGR))
   90    CONTINUE
         IPIJK=1
         IF(IPHASE.EQ.2) THEN
            CALL DOORPV(CDOOR,JPSYS,NPSYS,IPTRK,IFTRAK,IPRNTP,NGROUP,
     >      NREGIO,NBMIX,NANIST-1,MATCOD,VOLUME,KNORM,IPIJK,LEAKSW,
     >      ITPIJ,LNORM,TITRE,NALBP)
         ELSE
            CALL DOORAV(CDOOR,JPSYS,NPSYS,IPTRK,IFTRAK,IPRNTP,NGROUP,
     >      NREGIO,NBMIX,NANIST-1,NW,MATCOD,VOLUME,KNORM,LEAKSW,TITRE,
     >      NALBP,ISTRM)
         ENDIF
         CALL LCMSIX(IPSYS,' ',ILCMDN)
      ENDIF
*----
*  COMPUTE RESPONSE MATRIX INFORMATION FOR THE TIBERE
*  ANISOTROPIC STREAMING MODEL.
*----
      IF((ISTRM.EQ.3).AND.(IPHASE.EQ.1)) THEN
         CALL LCMSIX(IPSYS,'STREAMING',ILCMUP)
         JPSYS=LCMLID(IPSYS,'GROUP',NGROUP)
         IF(NANIST.LE.1) CALL XABORT('ASMDRV: MISSING P1 XS INFO.')
         DO 100 IGR=1,NGROUP
         NPSYS(IGR)=IGR
         KPSYS=LCMDIL(JPSYS,IGR)
         CALL LCMPUT(KPSYS,'DRAGON-TXSC',NBMIX+1,2,XSSIGT(0,IGR,1))
         CALL LCMPUT(KPSYS,'DRAGON-S0XSC',(NBMIX+1)*(NANIST-1),2,
     >   XSSIGW(0,2,IGR))
         IF(LDIFF) CALL LCMPUT(KPSYS,'DRAGON-DIFF',NBMIX+1,2,
     >   XSDIFF(0,IGR))
  100    CONTINUE
         CALL DOORAV(CDOOR,JPSYS,NPSYS,IPTRK,IFTRAK,IPRNTP,NGROUP,
     >   NREGIO,NBMIX,NANIST-1,NW,MATCOD,VOLUME,KNORM,LEAKSW,TITRE,
     >   NALBP,ISTRM)
         CALL LCMSIX(IPSYS,' ',ILCMDN)
      ENDIF
      IF(LNORM) LEAKSW=.FALSE.
*----
*  SCRATCH STORAGE DEALLOCATION
*----
      DEALLOCATE(ALBP)
      DEALLOCATE(ENERGY,ESTOPW,EMOMTR,XSDIFF,XSSCOR,XSSIGW,XSSIGT)
      DEALLOCATE(NPSYS)
      RETURN
*----
*  FORMATS
*----
 5090 FORMAT(32HASMDRV: INVALID VALUE OF NBMIX (,I5,2H).)
 6000 FORMAT(//30X,' EDITION REGION/VOLUME/MIXTURE '//
     >3(5X,'REGION',5X,'VOLUME  ',5X,'MIXTURE')/)
 6001 FORMAT(1P,3(5X,I4,4X,E12.5,4X,I4,4X))
 6002 FORMAT(//30X,' G R O U P : ',I5//31X,
     >'TOTAL AND WITHIN-GROUP MACROSCOPIC CROSS SECTIONS PER MIXTURE '/)
 6003 FORMAT(3(1X,'MIXTURE',4X,'NTOT0',11X,'SIGW',3X)/
     >1P,3(1X,I4,3X,E12.5,3X,E12.5))
 6101 FORMAT(//' USE TRANSPORT CORRECTED CROSS-SECTIONS (ITRANC=',I4,
     >' )')
 6200 FORMAT(//' COMPUTATION OF DRAGON RESPONSE MATRICES BY DOOR =',
     >3X,A12)
 6201 FORMAT(//' COMPUTATION OF DRAGON COMPLETE CP BY DOOR =',
     >3X,A12)
      END