summaryrefslogtreecommitdiff
path: root/Donjon/src/FLPDRV.f
blob: c57c28be85b757c6170c800ba43915ec885c963c (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
*DECK FLPDRV
      SUBROUTINE FLPDRV(IPPOW,IPNFX,IPFLX,IPKIN,IPTRK,IPMTX,IPMAP,
     1 IPMAC,PTOT,LNEW,LMAP,JMOD,LFLX,LPOW,LRAT,IMPX,FSTH,LFSTH,LFLU,
     2 LBUN,LNRM)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Driver for the powers and fluxes computations.
*
*Copyright:
* Copyright (C) 2007 Ecole Polytechnique de Montreal.
*
*Author(s): 
* D. Sekki, M. Guyot
*
*Parameters: input/output
* IPPOW  pointer to power information.
* IPNFX  pointer to normalized flux information.
* IPFLX  pointer to flux information.
* IPKIN  pointer to kinetics information.
* IPTRK  pointer to tracking information.
* IPMTX  pointer to matex information.
* IPMAP  pointer to fuel-map information.
* IPMAC  pointer to macrolib information.
* PTOT   given total reactor power in mega-watts.
* LNEW   new total power flag (=.true. for new computation).
* LMAP   fuel-map printing on file flag (=.true. for print).
* JMOD   modification index for L_MAP object.
* LFLX   flux printing on file flag (=.true. for print).
* LPOW   power printing on file flag (=.true. for print).
* LRAT   flux-ratio printing on file flag (=.true. for print).
* IMPX   printing on screen index (=0 for no print).
* FSTH   thermal to fission power ratio
* LFSTH  =.true if the thermal fission ratio is specified
* LFLU   =.true. if an output flux is to be created
* LBUN   =.true. if the output flux is a flux per bundle
* LNRM   =.true. if the output flux is a flux per mesh-splitted element
*
*-----------------------------------------------------------------------
*
      USE GANLIB
*----
*  SUBROUTINE ARGUMENTS
*----
      TYPE(C_PTR) IPPOW,IPNFX,IPFLX,IPKIN,IPTRK,IPMTX,IPMAP,IPMAC
      INTEGER JMOD,IMPX
      LOGICAL LNEW,LMAP,LFLX,LPOW,LRAT,LFSTH,LFLU,LBUN,LNRM
      REAL PTOT,FSTH
*----
*  LOCAL VARIABLES
*----
      PARAMETER(NSTATE=40)
      INTEGER ISTATE(NSTATE),IGR,NUN,NGRP,IGEO
      DOUBLE PRECISION ZNRM,VTOT,POWR
      CHARACTER HSIGN*12
      TYPE(C_PTR) JPMAC,KPMAC,MPFLUX
      INTEGER, ALLOCATABLE, DIMENSION(:) :: MAT,IDL,FMIX
      REAL, ALLOCATABLE, DIMENSION(:) :: HFC,VECT,FLUX,VOL,FXYZ,RAT,
     1 PXYZ,VMAP,FMAP,POWC,POWB,VECNR
*----
*  CHECK THE TYPE OF FLUX SELECTED FOR OUTPUT
*----
      IF(LFLU)THEN
        IF(LNRM.AND.LBUN) CALL XABORT('@FLPDRV: KEYWORD NORM AND BUND '
     1    //'BOTH SELECTED.')
        IF((.NOT.LNRM).AND.(.NOT.LBUN)) THEN
          LNRM=.TRUE.
          WRITE(6,*) 'FLPDRV: default option for L_FLUX object is NORM.'
        ENDIF
      ENDIF
*----
*  RECOVER INFORMATION
*----
      ISTATE(:NSTATE)=0
      IF(C_ASSOCIATED(IPFLX)) THEN
*       L_FLUX object
        CALL LCMGET(IPFLX,'STATE-VECTOR',ISTATE)
        NGRP=ISTATE(1)
        NUN=ISTATE(2)
        CALL LCMGET(IPFLX,'K-EFFECTIVE',FKEFF)
      ELSE IF(C_ASSOCIATED(IPKIN)) THEN
*       L_KINET object
        CALL LCMGET(IPKIN,'STATE-VECTOR',ISTATE)
        NGRP=ISTATE(3)
        NUN=ISTATE(6)
        CALL LCMGET(IPKIN,'E-KEFF',FKEFF)
        CALL LCMGET(IPKIN,'E-POW',PTOT)
      ENDIF
      ISTATE(:NSTATE)=0
      CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE)
      NEL=ISTATE(1)
      IF(ISTATE(2).NE.NUN) CALL XABORT('@FLPDRV: INCOMPATIBLE L_TRACK '
     1 //'AND L_FLUX/L_KINET OBJECTS')
      LX=ISTATE(14)
      LY=ISTATE(15)
      LZ=ISTATE(16)
*----
*  RECOVER H-FACTOR
*----
      ISTATE(:NSTATE)=0
      IF(C_ASSOCIATED(IPMAC))THEN
        CALL LCMGET(IPMAC,'STATE-VECTOR',ISTATE)
        IF(ISTATE(1).NE.NGRP)CALL XABORT('@FLPDRV: INVALID NUM'
     1   //'BER OF ENERGY GROUPS IN MACROLIB.')
        NMIX=ISTATE(2)
        ALLOCATE(HFC(NMIX*NGRP))
        JPMAC=LCMGID(IPMAC,'GROUP')
        DO JGR=1,NGRP
          KPMAC=LCMGIL(JPMAC,JGR)
          CALL LCMLEN(KPMAC,'H-FACTOR',LENGT,ITYP)
          IF(LENGT.NE.NMIX)CALL XABORT('@FLPDRV: UNABLE TO FIND'
     1     //' H-FACTOR BLOCK DATA IN THE MACROLIB.')
          CALL LCMGET(KPMAC,'H-FACTOR',HFC((JGR-1)*NMIX+1))
        ENDDO
      ELSEIF(C_ASSOCIATED(IPMTX))THEN
        CALL LCMGET(IPMTX,'STATE-VECTOR',ISTATE)
        IF(ISTATE(1).NE.NGRP)CALL XABORT('@FLPDRV: INVALID NUM'
     1   //'BER OF ENERGY GROUPS IN MATEX.')
        NMIX=ISTATE(2)
        ALLOCATE(HFC(NMIX*NGRP))
        HFC(:NMIX*NGRP)=0.0
        CALL LCMGET(IPMTX,'H-FACTOR',HFC)
      ENDIF
*----
*  RECOVER FUELMAP AND MATEX INFORMATION
*----
      IF(C_ASSOCIATED(IPMAP).AND.C_ASSOCIATED(IPMTX))THEN
        CALL LCMGET(IPMAP,'STATE-VECTOR',ISTATE)
        NB=ISTATE(1)
        NCH=ISTATE(2)
        IGEO=ISTATE(12)
        IF((IGEO.NE.7).AND.(IGEO.NE.9))CALL XABORT('@FLPDRV: INVALID'
     1   //' GEOMETRY IN FUEL MAP : ONLY 3-D CARTESIAN OR 3-D HEXAGO'
     2   //'NAL GEOMETRIES AVAILABLE') 
        IF(ISTATE(4).NE.NGRP)CALL XABORT('@FLPDRV: INVALID NUM'
     1   //'BER OF ENERGY GROUPS IN FUEL MAP OR FLUX.')
        ISTATE(:NSTATE)=0
        CALL LCMGET(IPMTX,'STATE-VECTOR',ISTATE)
        NMIX=ISTATE(2)
        NTOT=ISTATE(5)
        IF(ISTATE(6).NE.IGEO)CALL XABORT('@FLPDRV: GEOMETRIES IN'
     1   //' MATEX AND FUEL MAP ARE DIFFERENT')
        IF(ISTATE(7).NE.NEL)CALL XABORT('@FLPDRV: INVALID TOTAL'
     1   //' NUMBER OF REGIONS IN FUEL MAP OR TRACK.')
      ENDIF
*----
*  FLUX NORMALIZATION
*----
      ALLOCATE(VECT(NUN*NGRP),FLUX(NEL*NGRP),MAT(NEL),VOL(NEL),IDL(NEL))
      IF(LNEW)THEN
*       NEW TOTAL REACTOR POWER
        CALL LCMLEN(IPPOW,'NORM',LENGT,ITYP)
        IF(LENGT.EQ.0)CALL XABORT('@FLPDRV: UNABLE TO F'
     1   //'IND FLUX NORMALIZATION FACTOR IN L_POWER.')
        CALL LCMGET(IPPOW,'NORM',ZNRM)
        CALL FLPTOT(IPFLX,IPKIN,IPTRK,NMIX,NGRP,NEL,NUN,VECT,FLUX,MAT,
     1  VOL,IDL,HFC,POWR,ZNRM,IMPX)
      ELSE
*       GIVEN TOTAL REACTOR POWER
        POWR=DBLE(PTOT*10**6)
        IF(PTOT.EQ.0.0)CALL XABORT('@FLPDRV: PTOT IS NOT DEFINED')
        CALL FLPNRM(IPFLX,IPKIN,IPTRK,NMIX,NGRP,NEL,NUN,VECT,FLUX,MAT,
     1  VOL,IDL,HFC,POWR,ZNRM,IMPX)
        CALL LCMPUT(IPPOW,'NORM',1,4,ZNRM)
      ENDIF
      DEALLOCATE(IDL)
      CALL LCMPUT(IPPOW,'FLUX',NEL*NGRP,2,FLUX)
      POWR=POWR/(10**6)
      CALL LCMPUT(IPPOW,'PTOT',1,4,POWR)
*----
*  WHOLE REACTOR
*----
      ALLOCATE(FXYZ(NEL*NGRP),RAT(NEL*(NGRP-1)))
*     FLUX DISTRIBUTION
      IF(IGEO.EQ.7) THEN
        CALL FLPFLX(NGRP,NEL,LX,LY,LZ,MAT,VOL,FLUX,FXYZ,RAT,VTOT,IMPX,
     1  LFLX,LRAT)
      ELSEIF(IGEO.EQ.9) THEN
        CALL FLPHFX(NGRP,NEL,LX,LZ,MAT,VOL,FLUX,FXYZ,RAT,VTOT,IMPX,
     1  LFLX,LRAT)
      ENDIF
      CALL LCMPUT(IPPOW,'VTOT',1,4,VTOT)
      CALL LCMPUT(IPPOW,'FLUX-DISTR',NEL*NGRP,2,FXYZ)
      IF(NGRP.GT.1) CALL LCMPUT(IPPOW,'FLUX-RATIO',NEL*(NGRP-1),2,RAT)
      DEALLOCATE(RAT,FXYZ)
*     POWER DISTRIBUTION
      ALLOCATE(PXYZ(NEL))
      IF(IGEO.EQ.7) THEN
        CALL FLPOWR(NMIX,NGRP,NEL,LX,LY,LZ,MAT,VOL,FLUX,HFC,PXYZ,VTOT,
     1  IMPX,LPOW)
      ELSEIF(IGEO.EQ.9) THEN
        CALL FLPHPW(NMIX,NGRP,NEL,LX,LZ,MAT,VOL,FLUX,HFC,PXYZ,VTOT,
     1  IMPX,LPOW)
      ENDIF
      CALL LCMPUT(IPPOW,'POWER-DISTR',NEL,2,PXYZ)
      DEALLOCATE(PXYZ)
      CALL LCMPUT(IPPOW,'K-EFFECTIVE',1,2,FKEFF)
*----
*  FUEL-MAP
*----
      IF((C_ASSOCIATED(IPMAP)).AND.(C_ASSOCIATED(IPMTX))) THEN
        ALLOCATE(FMIX(NCH*NB),VMAP(NCH*NB),FMAP(NCH*NB*NGRP))
        CALL LCMGET(IPMAP,'FLMIX',FMIX)
*       COMPUTE FLUXES
        CALL FLPFLB(IPMTX,NTOT,NGRP,NEL,NCH,NB,FLUX,VOL,FMIX,VMAP,FMAP,
     1  IMPX,LMAP)
        ALLOCATE(POWC(NCH),POWB(NCH*NB))
*       COMPUTE POWERS
        CALL FLPOWB(IPPOW,IPMAP,IPMTX,NMIX,NTOT,NGRP,NCH,NB,NEL,MAT,
     1  VOL,HFC,FLUX,POWB,POWC,IMPX,POWR,FSTH,LFSTH,FMIX,FMAP,IGEO)
        CALL LCMPUT(IPPOW,'POWER-CHAN',NCH,2,POWC)
        CALL LCMPUT(IPPOW,'POWER-BUND',NCH*NB,2,POWB)
        CALL LCMPUT(IPPOW,'FLUX',NEL*NGRP,2,FLUX)
        CALL LCMPUT(IPPOW,'VOLU-BUND',NCH*NB,2,VMAP)
        CALL LCMPUT(IPPOW,'FLUX-BUND',NCH*NB*NGRP,2,FMAP)
        CALL LCMPUT(IPPOW,'FLMIX',NCH*NB,1,FMIX)
        IF(JMOD.GE.1) THEN
          CALL LCMPUT(IPMAP,'TOT-PW',1,2,PTOT)
          CALL LCMPUT(IPMAP,'BUND-PW',NCH*NB,2,POWB)
          CALL LCMPUT(IPMAP,'FLUX-AV',NCH*NB*NGRP,2,FMAP)
          IF(JMOD.EQ.2) THEN
            CALL LCMPUT(IPMAP,'BUND-PW-INI',NCH*NB,2,POWB)
          ENDIF
          PTOT=0.0
          DO I=1,NCH*NB
            PTOT=PTOT+POWB(I)
          ENDDO
          PTOT=PTOT/1.0E3
          CALL LCMPUT(IPMAP,'REACTOR-PW',1,2,PTOT)
        ENDIF
        DEALLOCATE(POWB,POWC)
      ENDIF
      DEALLOCATE(VOL,MAT,FLUX,HFC)

      IF(LFLU) THEN
*----
*  STATE-VECTOR FOR L_FLUX OBJECT
*----
        IF(LNRM) THEN
          ISTATE(:NSTATE)=0
          IF(C_ASSOCIATED(IPFLX)) THEN
            CALL LCMGET(IPFLX,'STATE-VECTOR',ISTATE)
          ELSE IF(C_ASSOCIATED(IPKIN)) THEN
            ISTATE(1)=NGRP
            ISTATE(2)=NUN
          ENDIF
          CALL LCMPUT(IPNFX,'STATE-VECTOR',NSTATE,1,ISTATE)
          HSIGN='L_FLUX'
          CALL LCMPTC(IPNFX,'SIGNATURE',12,HSIGN)
          ALLOCATE(VECNR(NUN*NGRP))
          DO 10 I=1,NUN*NGRP
            VECNR(I)=VECT(I)*REAL(ZNRM)
   10     CONTINUE
          MPFLUX=LCMLID(IPNFX,'FLUX',NGRP)
          DO 20 IGR=1,NGRP
            IOFSET=(IGR-1)*NUN
            CALL LCMPDL(MPFLUX,IGR,NUN,2,VECNR(IOFSET+1))
   20     CONTINUE        
          DEALLOCATE(VECNR)
        ELSE
          MPFLUX=LCMLID(IPNFX,'FLUX',NGRP)
          DO 30 IGR=1,NGRP
              IOFSET=(IGR-1)*NB*NCH
              CALL LCMPDL(MPFLUX,IGR,NB*NCH,2,FMAP(IOFSET+1))
   30     CONTINUE 
        ENDIF
      ENDIF
      IF(C_ASSOCIATED(IPMAP)) DEALLOCATE(VMAP,FMAP,FMIX)
      DEALLOCATE(VECT)
*----
*  STATE-VECTOR FOR L_FLUX OBJECT
*----
      IF(LFLU) THEN
        ISTATE(:NSTATE)=0
        IF(C_ASSOCIATED(IPFLX)) THEN
          CALL LCMGET(IPFLX,'STATE-VECTOR',ISTATE)
        ELSE IF(C_ASSOCIATED(IPKIN)) THEN
          ISTATE(1)=NGRP
          ISTATE(2)=NUN
        ENDIF
        IF(LBUN) ISTATE(2)=NB*NCH
        CALL LCMPUT(IPNFX,'STATE-VECTOR',NSTATE,1,ISTATE)
        HSIGN='L_FLUX'
        CALL LCMPTC(IPNFX,'SIGNATURE',12,HSIGN)
      ENDIF
*----
*  STATE-VECTOR FOR L_POWER OBJECT
*----
      ISTATE(:NSTATE)=0
      ISTATE(1)=NGRP
      ISTATE(2)=NEL
      ISTATE(3)=LX
      ISTATE(4)=LY
      ISTATE(5)=LZ
      ISTATE(6)=NCH
      ISTATE(7)=NB
      ISTATE(8)=IGEO
      CALL LCMPUT(IPPOW,'STATE-VECTOR',NSTATE,1,ISTATE)
      HSIGN='L_POWER'
      CALL LCMPTC(IPPOW,'SIGNATURE',12,HSIGN)
      IF(IMPX.GT.1)CALL LCMLIB(IPPOW)
      RETURN
      END