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
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
|
*DECK NAPCPO
SUBROUTINE NAPCPO(IPCPO,IPTRK,IPFLU,NSTATE)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Construct an 'enriched' multicompo with additional information
* needed by Pin Power Reconstruction.
*
*Copyright:
* Copyright (C) 2014 Ecole Polytechnique de Montreal.
*
*Author(s):
* R. Chambon
*
*Parameters: input/output
* IPCPO LCM object address of Multicompo.
* IPTRK LCM object address of Tracking.
* IPFLU LCM object address of Flux.
* NSTATE length of the state vector
*
*-----------------------------------------------------------------------
*
USE GANLIB
IMPLICIT NONE
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER NSTATE
TYPE(C_PTR) IPCPO,IPTRK,IPFLU
*----
* LOCAL VARIABLES
*----
INTEGER NCRCAL,NGPT
INTEGER IOUT,MAXPAR,MAXLIN,MAXVAL,MAXADD,MAXIFX
REAL REPS
PARAMETER (REPS=1.0E-4,IOUT=6,MAXPAR=50,MAXLIN=50,MAXVAL=200,
1 MAXADD=10,MAXIFX=5,NGPT=2)
CHARACTER PARKEY(MAXPAR)*12,PARFMT(MAXPAR)*8,RECNAM*12,
1 COMMEN(MAXLIN)*80,PARKEL(MAXPAR)*12,VALH(MAXPAR)*12
TYPE(C_PTR) JPCPO,KPCPO,LPCPO,JPMIC,JPFLU
CHARACTER TEXT*12,HSMG*131,DIRHOM*12,VCHAR(MAXVAL)*12,HVECT*8
INTEGER ISTATE(NSTATE),NPAR,NLOC,IMPX,IEL,NFDI,FINF(MAXIFX)
INTEGER IPAR,IBMOLD,IFX
REAL FLOT
DOUBLE PRECISION DFLOT
INTEGER VALI(MAXPAR),NVALUE(MAXPAR),VINTE(MAXVAL),
1 MUPLET(2*MAXPAR)
REAL VALR(2*MAXPAR,2),VREAL(MAXVAL),NVPO(2),PTR,PDF,PDF2
REAL ZGKSIX(NGPT),ZGKSIY(NGPT),WGKSIX(NGPT), WGKSIY(NGPT),
1 FLUGP(NGPT,NGPT),FPD,DX,DY
INTEGER NGFF,NXP,NYP,ITYPGP,NMIXP,NMIL,NG,NCOMLI,MAXNVP,STYPP,
1 NMCAL
INTEGER I,J,ICAL,INDIC,ITYLCM,IX,IY,LENGTH,NITMA,IREG,IREGP,IG,
1 IMIXP,ID,JD,IGP,JGP,IP,JP,J1,ICHX,IDIM,LC,L4,MAXKN,MKN
INTEGER NREG,NUN,NXD,NYD,ITYPGD,NREGP
REAL E(25)
LOGICAL LNOINT,FLAG
*----
* ALLOCATABLE ARRAYS
*----
INTEGER, ALLOCATABLE, DIMENSION(:) :: MIXP
REAL, ALLOCATABLE, DIMENSION(:) :: MXP,MYP,KN
REAL, ALLOCATABLE, DIMENSION(:) :: MXD,MYD,XX,YY
INTEGER, ALLOCATABLE, DIMENSION(:) :: KEYFLX,MATCOD,IXPD,
1 IYPD,JDEBAR,JARBVA
REAL, ALLOCATABLE, DIMENSION(:) :: FLXD
REAL, ALLOCATABLE, DIMENSION(:,:) :: FLXP,FT,FXTD,FYTD
LOGICAL, ALLOCATABLE, DIMENSION(:) :: FLAGMX
*
IMPX=0
CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT)
IF(INDIC.NE.3) CALL XABORT('NAPCPO: character data expected.')
IF(TEXT.EQ.'EDIT') THEN
CALL REDGET(INDIC,IMPX,FLOT,TEXT,DFLOT)
IF(INDIC.NE.1) CALL XABORT('NAPCPO: integer data expected.')
CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT)
IF(INDIC.NE.3) CALL XABORT('NAPCPO: character data expected.')
ENDIF
IF(TEXT.NE.'PROJECTION') CALL XABORT('NAPCPO: ''PROJECTION'' '//
1 'EXPECTED.')
CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT)
IF(TEXT.NE.'STEP') CALL XABORT('NAPCPO: ''STEP'' '//
1 'EXPECTED.')
CALL REDGET(INDIC,NITMA,FLOT,DIRHOM,DFLOT)
IF(INDIC.NE.3) CALL XABORT('NAPCPO: character data expected.')
CALL LCMSIX(IPCPO,DIRHOM,1)
*
IFX=1
LNOINT=.FALSE.
FINF(:MAXIFX)=-1
*----
* RECOVER TABLE-OF-CONTENT INFORMATION FOR THE COMPO.
*----
CALL LCMGET(IPCPO,'STATE-VECTOR',ISTATE)
NMIL=ISTATE(1)
NG=ISTATE(2)
NMCAL=ISTATE(4)
NPAR=ISTATE(5)
NLOC=ISTATE(6)
NCOMLI=ISTATE(10)
NGFF=ISTATE(14)
IF(NGFF.EQ.0) CALL XABORT('NAPCPO: NO GFF INFO IN MULTICOMPO.')
CALL LCMGTC(IPCPO,'COMMENT',80,NCOMLI,COMMEN)
IF(NPAR.GT.0)THEN
CALL LCMSIX(IPCPO,'GLOBAL',1)
CALL LCMGTC(IPCPO,'PARKEY',12,NPAR,PARKEY)
CALL LCMGTC(IPCPO,'PARFMT',8,NPAR,PARFMT)
CALL LCMGET(IPCPO,'NVALUE',NVALUE)
IF(IMPX.GT.10)THEN
DO IPAR=1,NPAR
WRITE(RECNAM,'(''pval'',I8.8)') IPAR
IF(PARFMT(IPAR).EQ.'INTEGER') THEN
CALL LCMGET(IPCPO,RECNAM,VINTE)
WRITE(IOUT,'(13H NAPCPO: KEY=,A,18H TABULATED POINTS=,
1 1P,6I12/(43X,6I12))') PARKEY(IPAR),(VINTE(I),I=1,
2 NVALUE(IPAR))
ELSE IF(PARFMT(IPAR).EQ.'REAL') THEN
CALL LCMGET(IPCPO,RECNAM,VREAL)
WRITE(IOUT,'(13H NAPCPO: KEY=,A,18H TABULATED POINTS=,
1 1P,6E12.4/(43X,6E12.4))') PARKEY(IPAR),(VREAL(I),I=1,
2 NVALUE(IPAR))
ELSE IF(PARFMT(IPAR).EQ.'STRING') THEN
CALL LCMGTC(IPCPO,RECNAM,12,NVALUE(IPAR),VCHAR)
WRITE(IOUT,'(13H NAPCPO: KEY=,A,18H TABULATED POINTS=,
1 1P,6A12/(43X,6A12))') PARKEY(IPAR),(VCHAR(I),I=1,
2 NVALUE(IPAR))
ENDIF
ENDDO
ENDIF
CALL LCMSIX(IPCPO,' ',2)
ENDIF
IF(NLOC.GT.0)THEN
CALL LCMSIX(IPCPO,'LOCAL',1)
CALL LCMGTC(IPCPO,'PARKEY',12,NLOC,PARKEL)
CALL LCMSIX(IPCPO,' ',2)
JPCPO=LCMGID(IPCPO,'MIXTURES')
DO IBMOLD=1,NMIL
KPCPO=LCMGIL(JPCPO,IBMOLD)
LPCPO=LCMGID(KPCPO,'TREE')
CALL LCMGET(LPCPO,'NVALUE',NVALUE)
IF(IMPX.GT.10)THEN
WRITE(IOUT,'(17H NAPCPO: MIXTURE=,I6)') IBMOLD
DO IPAR=1,NLOC
WRITE(RECNAM,'(''pval'',I8.8)') IPAR
CALL LCMGET(LPCPO,RECNAM,VREAL)
WRITE(IOUT,'(13H NAPCPO: KEY=,A,18H TABULATED POINTS=,
1 1P,6E12.4/(43X,6E12.4))') PARKEL(IPAR),(VREAL(I),I=1,
2 NVALUE(IPAR))
ENDDO
ENDIF
ENDDO
ENDIF
IF(IMPX.GT.10)WRITE(IOUT,'(1X,A)') (COMMEN(I),I=1,NCOMLI)
*----
* READ (INTERP_DATA) AND SET VALI, VALR AND VALH PARAMETERS
* CORRESPONDING TO THE INTERPOLATION POINT. FILL MUPLET FOR
* PARAMETERS.
*----
MUPLET(:NPAR+NLOC)=0
1020 CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT)
IF(INDIC.NE.3) CALL XABORT('NAPCPO: character data expected.')
IF(TEXT.EQ.'SET') THEN
CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT)
IF(INDIC.NE.3) CALL XABORT('NAPCPO: character data expected.')
DO 50 I=1,NPAR
IF(TEXT.EQ.PARKEY(I)) THEN
IPAR=I
GO TO 60
ENDIF
50 CONTINUE
GO TO 100
60 LPCPO=LCMGID(IPCPO,'GLOBAL')
CALL LCMGET(LPCPO,'NVALUE',NVALUE)
IF(NVALUE(IPAR).GT.MAXVAL) CALL XABORT('NAPCPO: MAXVAL OVERFL'
1 //'OW.')
WRITE(RECNAM,'(''pval'',I8.8)') IPAR
CALL LCMLEN(LPCPO,RECNAM,LENGTH,ITYLCM)
IF(LENGTH.EQ.0) THEN
WRITE(HSMG,'(25HNAPCPO: GLOBAL PARAMETER ,A,9H NOT SET.)')
1 PARKEY(IPAR)
CALL XABORT(HSMG)
ENDIF
IF(PARFMT(IPAR).EQ.'INTEGER') THEN
CALL REDGET(INDIC,VALI(IPAR),FLOT,TEXT,DFLOT)
IF(INDIC.NE.1) CALL XABORT('NAPCPO: integer data expected.')
CALL LCMGET(LPCPO,RECNAM,VINTE)
DO J=1,NVALUE(IPAR)
IF(VALI(IPAR).EQ.VINTE(J)) THEN
MUPLET(IPAR)=J
* MUTYPE(IPAR)=ITYPGD
GO TO 1020
ENDIF
ENDDO
WRITE(HSMG,'(26HNAPCPO: INTEGER PARAMETER ,A,9H WITH VAL,
1 2HUE,I5,29H NOT FOUND IN COMPO DATABASE.)') PARKEY(IPAR),
2 VALI(IPAR)
CALL XABORT(HSMG)
ELSE IF(PARFMT(IPAR).EQ.'REAL') THEN
CALL REDGET(INDIC,NITMA,VALR(IPAR,1),TEXT,DFLOT)
IF(INDIC.NE.2) CALL XABORT('NAPCPO: real data expected.')
VALR(IPAR,2)=VALR(IPAR,1)
! CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT)
CALL LCMGET(LPCPO,RECNAM,VREAL)
DO J=1,NVALUE(IPAR)
IF(ABS(VALR(IPAR,1)-VREAL(J)).LE.REPS*ABS(VREAL(J))) THEN
MUPLET(IPAR)=J
* MUTYPE(IPAR)=ITYPGD
GO TO 1020
ENDIF
ENDDO
WRITE(HSMG,'(23HNAPCPO: REAL PARAMETER ,A,9H WITH VAL,
1 2HUE,I5,29H NOT FOUND IN COMPO DATABASE.)') PARKEY(IPAR),
2 VALR(IPAR,1)
CALL XABORT(HSMG)
ELSE IF(PARFMT(IPAR).EQ.'STRING') THEN
CALL REDGET(INDIC,NITMA,FLOT,VALH(IPAR),DFLOT)
IF(INDIC.NE.3) CALL XABORT('NAPCPO: STRING DATA EXPECTED.')
CALL LCMGTC(LPCPO,RECNAM,12,NVALUE(IPAR),VCHAR)
DO J=1,NVALUE(IPAR)
IF(VALH(IPAR).EQ.VCHAR(J)) THEN
MUPLET(IPAR)=J
* MUTYPE(IPAR)=ITYPGD
GO TO 1020
ENDIF
ENDDO
WRITE(HSMG,'(25HNAPCPO: STRING PARAMETER ,A,10H WITH VALU,
1 2HE ,A12,29H NOT FOUND IN COMPO DATABASE.)') PARKEY(IPAR),
2 VALH(IPAR)
CALL XABORT(HSMG)
ENDIF
100 DO 110 I=1,NLOC
IF(TEXT.EQ.PARKEL(I)) THEN
IPAR=NPAR+I
GO TO 120
ENDIF
110 CONTINUE
CALL XABORT('NAPCPO: PARAMETER '//TEXT//' NOT FOUND.')
120 JPCPO=LCMGID(IPCPO,'MIXTURES')
IBMOLD=1
KPCPO=LCMGIL(JPCPO,IBMOLD)
LPCPO=LCMGID(KPCPO,'TREE')
CALL LCMGET(LPCPO,'NVALUE',NVALUE)
CALL REDGET(INDIC,NITMA,VALR(IPAR,1),TEXT,DFLOT)
IF(INDIC.NE.2) CALL XABORT('NAPCPO: real data expected.')
VALR(IPAR,2)=VALR(IPAR,1)
WRITE(RECNAM,'(''pval'',I8.8)') IPAR-NPAR
CALL LCMLEN(LPCPO,RECNAM,LENGTH,ITYLCM)
IF(LENGTH.EQ.0) THEN
WRITE(HSMG,'(24HNAPCPO: LOCAL PARAMETER ,A,9H NOT SET.)')
1 PARKEL(IPAR-NPAR)
CALL XABORT(HSMG)
ELSE IF(LENGTH.GT.MAXVAL) THEN
CALL XABORT('NAPCPO: MAXVAL OVERFLOW.')
ENDIF
CALL LCMGET(LPCPO,RECNAM,VREAL)
DO J=1,NVALUE(IPAR-NPAR)
IF(ABS(VALR(IPAR,1)-VREAL(J)).LE.REPS*ABS(VREAL(J))) THEN
MUPLET(IPAR)=J
* MUTYPE(IPAR)=ITYPGD
GO TO 1020
ENDIF
ENDDO
WRITE(HSMG,'(26HNAPCPO: INTEGER PARAMETER ,A,9H WITH VAL,
1 2HUE,I5,29H NOT FOUND IN COMPO DATABASE.)') PARKEY(IPAR),
2 VALI(IPAR)
CALL XABORT(HSMG)
ELSEIF(TEXT.EQ.'IFX') THEN
CALL REDGET(INDIC,IFX,FLOT,TEXT,DFLOT)
IF(INDIC.NE.1) CALL XABORT('NAPCPO: integer data expected.')
GO TO 1020
ELSEIF(TEXT.EQ.'NOINTP') THEN
LNOINT=.TRUE.
GO TO 1020
ELSEIF(TEXT.EQ.'INTERP') THEN
LNOINT=.FALSE.
GO TO 1020
ELSEIF(TEXT.EQ.';') THEN
GOTO 200
ENDIF
CALL XABORT('NAPCPO: '//TEXT//' is a wrong keyword')
*
200 CONTINUE
JPCPO=LCMGID(IPCPO,'MIXTURES')
IBMOLD=1
KPCPO=LCMGIL(JPCPO,IBMOLD)
LPCPO=LCMGID(KPCPO,'TREE')
CALL LCMGET(LPCPO,'NVP',NVPO)
CALL LCMLEN(LPCPO,'ARBVAL',MAXNVP,ITYLCM)
IF(NVPO(1).GT.MAXNVP) CALL XABORT('NAPCPO: NVP OVERFLOW.')
ALLOCATE(JDEBAR(MAXNVP+1),JARBVA(MAXNVP))
CALL LCMGET(LPCPO,'DEBARB',JDEBAR)
CALL LCMGET(LPCPO,'ARBVAL',JARBVA)
IF(IMPX.GE.20) THEN
WRITE(6,*) 'MUPLET: ',(MUPLET(I),I=1,NPAR+NLOC)
ENDIF
ICAL=NCRCAL(1,NVPO(1),NPAR+NLOC,JDEBAR,JARBVA,MUPLET)
IF(IMPX.GE.2) THEN
WRITE(6,*) 'Performing projection for calculation: ',ICAL
ENDIF
*
LPCPO=LCMGID(KPCPO,'CALCULATIONS')
JPMIC=LCMGIL(LPCPO,ICAL)
CALL LCMGET(JPMIC,'STATE-VECTOR',ISTATE)
CALL LCMSIX(JPMIC,'MACROLIB ',1)
CALL LCMSIX(JPMIC,'GFF ',1)
CALL LCMSIX(JPMIC,'GFF-GEOM ',1)
C get dimension in geometry from L_MULTICOMPO
CALL LCMGET(JPMIC,'STATE-VECTOR',ISTATE)
ITYPGP=ISTATE(1)
STYPP=ISTATE(11)
IF(ITYPGP.NE.5) CALL XABORT('NAPCPO: CAR2D geometry type '
1 //'expected in L_MULTICOMPO.')
IF(STYPP.NE.0) CALL XABORT('NAPCPO: No split in geometry expected'
1 //' in L_MULTICOMPO.')
NXP=ISTATE(3)
NYP=ISTATE(4)
NREGP=ISTATE(6)
NMIXP=ISTATE(7)
IF(NMIXP.NE.NGFF) CALL XABORT('NAPCPO: INVALID GFF-GEOM.')
ALLOCATE(MXP(NXP+1),MYP(NYP+1))
ALLOCATE(IXPD(NXP+1),IYPD(NYP+1))
ALLOCATE(MIXP(NREGP))
CALL LCMGET(JPMIC,'MESHX',MXP)
CALL LCMGET(JPMIC,'MESHY',MYP)
CALL LCMGET(JPMIC,'MIX',MIXP)
CALL LCMSIX(JPMIC,'GFF-GEOM ',2)
IXPD(:NXP+1)=0
IYPD(:NYP+1)=0
C get dimension in geometry from L_TRACK
CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE)
NREG=ISTATE(1)
NUN=ISTATE(2)
ITYPGD=ISTATE(6)
IEL=ISTATE(9)
IF(ITYPGD.NE.5) CALL XABORT('NAPCPO: CAR2D geometry type expected'
1 //' in L_TRACKING.')
IEL=ISTATE(9)
L4=ISTATE(11)
ICHX=ISTATE(12)
NXD=ISTATE(14)
NYD=ISTATE(15)
IDIM=2
IF(NREG.NE.NXD*NYD) CALL XABORT('NAPCPO: No Splitting allowed in '
1 //'CAR2D geometry type from L_TRACK.')
C compute X and Y mesh from L_TRACK
ALLOCATE(MXD(NXD+1),MYD(NYD+1))
ALLOCATE(XX(NREG),YY(NREG))
CALL LCMGET(IPTRK,'XX',XX)
CALL LCMGET(IPTRK,'YY',YY)
MXD(1)=MXP(1)
DO I=1,NXD
MXD(I+1)=MXD(I)+XX(I)
ENDDO
MYD(1)=MYP(1)
DO I=1,NYD
MYD(I+1)=MYD(I)+YY((I-1)*NXD+1)
ENDDO
if(IMPX.ge.10) then
WRITE(6,*) 'Respective mesh (Diffusion vs. Transport):'
WRITE(6,*) ' X direction :'
WRITE(6,*) 'MXD:',(MXD(I),I=1,NXD+1)
WRITE(6,*) 'MXP:',(MXP(I),I=1,NXP+1)
WRITE(6,*) ' Y direction :'
WRITE(6,*) 'MYD:',(MYD(I),I=1,NYD+1)
WRITE(6,*) 'MYP:',(MYP(I),I=1,NYP+1)
endif
IF((ABS(MXD(NXD+1)-MXP(NXP+1)).GE.1E-3).OR.
1 (ABS(MXD(NXD+1)-MXP(NXP+1)).GE.1E-3)) CALL XABORT('NAPCPO: '
2 //'Diffusion and transport geometries total size mismach')
ALLOCATE(FXTD(NXP,NXD),FYTD(NYP,NYD))
FXTD(:NXP,:NXD)=0.0
FYTD(:NYP,:NYD)=0.0
CALL NAPFTD(NXP,MXP,NXD,MXD,FXTD)
CALL NAPFTD(NYP,MYP,NYD,MYD,FYTD)
IF(LNOINT) THEN
C verify that both meshes match
J1=1
DO I=2,NXD+1
FLAG=.TRUE.
DO J=J1,NXP+1
IF(MXP(J).LT.MXD(I)) THEN
IXPD(J)=I-1
ENDIF
IF(ABS(MXD(I)-MXP(J)).LE.ABS(1E-5*MXP(J))) THEN
FLAG=.FALSE.
IXPD(J)=I
J1=J+1
ENDIF
ENDDO
IF(FLAG) CALL XABORT('NAPCPO: a X mesh in L_TRACK does not '
1 //'match the CAR2D geometry imbedded in L_MULTICOMPO.')
ENDDO
J1=1
DO I=2,NYD+1
FLAG=.TRUE.
DO J=J1,NYP+1
IF(MYP(J).LT.MYD(I)) THEN
IYPD(J)=I-1
ENDIF
IF(ABS(MYD(I)-MYP(J)).LE.ABS(1E-5*MYP(J))) THEN
FLAG=.FALSE.
IYPD(J)=I
J1=J+1
ENDIF
ENDDO
IF(FLAG) CALL XABORT('NAPCPO: a Y mesh in L_TRACK does not '
1 //'match the CAR2D geometry imbedded in L_MULTICOMPO.')
ENDDO
ENDIF
C project flux
ALLOCATE(KEYFLX(NREG),MATCOD(NREG))
CALL LCMGET(IPTRK,'KEYFLX',KEYFLX)
CALL LCMGET(IPTRK,'MATCOD',MATCOD)
ALLOCATE(FLXD(NUN),FLXP(NMIXP,NG))
ALLOCATE(FLAGMX(NMIXP))
JPFLU=LCMGID(IPFLU,'FLUX')
DO IG=1,NG
CALL LCMGDL(JPFLU,IG,FLXD)
DO IP=1,NXP
DO JP=1,NYP
IREGP=IP+(JP-1)*NXP
IF(IREGP.GT.NREGP) CALL XABORT('NAPCPO: NREGP OVERFLOW(1).')
IMIXP=MIXP(IREGP)
FLXP(IMIXP,IG)=0.0
IF(LNOINT) THEN
* integrated projected flux FLXP
IREG=IXPD(IP)+(IYPD(JP)-1)*NXD
FLXP(IMIXP,IG)=FLXD(KEYFLX(IREG))
ELSE
* interpolated projected flux FLXP
DO ID=1,NXD
DO JD=1,NYD
IF(FXTD(IP,ID)*FYTD(JP,JD).NE.0.0) THEN
* -----
CALL ALGPT(NGPT,MAX(MXP(IP),MXD(ID)),MIN(MXP(IP+1),MXD(ID+1)),
1 ZGKSIX,WGKSIX)
DX=MIN(MXP(IP+1),MXD(ID+1))-MAX(MXP(IP),MXD(ID))
CALL ALGPT(NGPT,MAX(MYP(JP),MYD(JD)),MIN(MYP(JP+1),MYD(JD+1)),
1 ZGKSIY,WGKSIY)
DY=MIN(MYP(JP+1),MYD(JD+1))-MAX(MYP(JP),MYD(JD))
IF(IMPX.GE.5) THEN
WRITE(6,*) 'IP,JP:',IP,JP,FXTD(IP,ID),'ID,JD:',ID,JD,FYTD(JP,JD)
WRITE(6,*) 'Gauss point ZGWG:',(ZGKSIX(I),I=1,NGPT),
1 (WGKSIX(I),I=1,NGPT),'DX',DX
WRITE(6,*) 'Gauss point ZGWG:',(ZGKSIY(I),I=1,NGPT),
1 (WGKSIY(I),I=1,NGPT),'DY',DY
ENDIF
FPD=0.0
* interpolate flux
IF(ICHX.EQ.1) THEN
* Variational collocation method
CALL LCMLEN(IPTRK,'KN',MAXKN,ITYLCM)
MKN=MAXKN/(NXD*NYD)
ALLOCATE(KN(MAXKN))
CALL LCMGET(IPTRK,'KN',KN)
CALL LCMSIX(IPTRK,'BIVCOL',1)
CALL LCMLEN(IPTRK,'T',LC,ITYLCM)
CALL LCMGET(IPTRK,'E',E)
CALL LCMSIX(IPTRK,' ',2)
CALL VALU2B(LC,MKN,NXD,NYD,L4,ZGKSIX,ZGKSIY,MXD,MYD,FLXD,MATCOD,
1 KN,NGPT,NGPT,E,FLUGP)
ELSE IF(ICHX.EQ.2) THEN
* Raviart-Thomas finite element method
CALL VALU4B(IEL,NUN,NXD,NYD,ZGKSIX,ZGKSIY,MXD,MYD,FLXD,MATCOD,
1 KEYFLX,NGPT,NGPT,FLUGP)
ELSE IF(ICHX.EQ.3) THEN
* Nodal collocation method (MCFD)
CALL VALU1B(IDIM,NXD,NYD,L4,ZGKSIX,ZGKSIY,MXD,MYD,FLXD,MATCOD,
1 IEL,NGPT,NGPT,FLUGP)
ELSE
CALL XABORT('NAPCPO: INTERPOLATION NOT IMPLEMENTED.')
ENDIF
IF(IMPX.GE.5) THEN
WRITE(6,*) 'Gauss flux values:'
do JGP=1,NGPT
WRITE(6,*) (FLUGP(IGP,JGP),IGP=1,NGPT)
ENDDO
ENDIF
* integrate flux (gauss method)
DO IGP=1,NGPT
DO JGP=1,NGPT
FPD=FPD+FLUGP(IGP,JGP)*WGKSIX(IGP)*WGKSIY(JGP)
ENDDO
ENDDO
* get average flux
FPD=FPD/DX/DY
FLXP(IMIXP,IG)=FLXP(IMIXP,IG)+FPD*FXTD(IP,ID)*FYTD(JP,JD)
* -----
ENDIF
ENDDO
ENDDO
ENDIF
ENDDO
ENDDO
ENDDO
C flux normalization
C get data from transport calculations
ALLOCATE(FT(NMIXP,NG))
CALL LCMGET(JPMIC,'NWT0',FT)
C group by group
DO IG=1,NG
C compute average flux from transport calculations
PTR=0.0
IREGP=0
DO IY=1,NYP
DO IX=1,NXP
IREGP=IREGP+1
IF(IREGP.GT.NREGP) CALL XABORT('NAPCPO: NREGP OVERFLOW(2).')
IMIXP=MIXP(IREGP)
PTR=PTR+FT(IMIXP,IG)*(MXP(IX+1)-MXP(IX))*(MYP(IY+1)-MYP(IY))
ENDDO
ENDDO
C compute average flux with projected diffusion flux
PDF=0.0
IREGP=0
DO IY=1,NYP
DO IX=1,NXP
IREGP=IREGP+1
IF(IREGP.GT.NREGP) CALL XABORT('NAPCPO: NREGP OVERFLOW(3).')
IMIXP=MIXP(IREGP)
PDF=PDF+FLXP(IMIXP,IG)*(MXP(IX+1)-MXP(IX))*(MYP(IY+1)-MYP(IY))
ENDDO
ENDDO
C renormalize flux
DO IMIXP=1,NMIXP
FLXP(IMIXP,IG)=FLXP(IMIXP,IG)/PDF*PTR
ENDDO
C
IF(IMPX.GT.5) THEN
PDF2=0.0
IREGP=0
DO IY=1,NYP
DO IX=1,NXP
IREGP=IREGP+1
IF(IREGP.GT.NREGP) CALL XABORT('NAPCPO: NREGP OVERFLOW(4).')
IMIXP=MIXP(IREGP)
PDF2=PDF2+FLXP(IMIXP,IG)
1 *(MXP(IX+1)-MXP(IX))*(MYP(IY+1)-MYP(IY))
ENDDO
ENDDO
WRITE(6,*)'NAPCPO: transport power:',PTR
WRITE(6,*)'NAPCPO: diffusion power (before normalization):',PDF
WRITE(6,*)'NAPCPO: diffusion power (after normalization):',PDF2
IREGP=0
WRITE(6,*) 'NAPCPO: FLXP/FT: group #',IG
DO IY=1,NYP
WRITE(6,*) (FLXP(MIXP(IREGP+I),IG)
1 /FT(MIXP(IREGP+I),IG),I=1,NXP)
IREGP=IREGP+NYP
ENDDO
ENDIF
C verify that all mixtures have a projected flux
DO IMIXP=1,NMIXP
IF(FLXP(IMIXP,IG).EQ.0.0) THEN
WRITE(HSMG,'(42HNAPCPO: no projected flux for mixture and ,
1 6Hgroup=,2I6,1H.)') IMIXP,IG
CALL XABORT(HSMG)
ENDIF
ENDDO
C end DO IG=1,NG
ENDDO
C save projected flux in L_MULTICOMPO for each original mixture
DO I=1,NMIL
JPCPO=LCMGID(IPCPO,'MIXTURES')
KPCPO=LCMGIL(JPCPO,I)
LPCPO=LCMGID(KPCPO,'CALCULATIONS')
JPMIC=LCMGIL(LPCPO,ICAL)
CALL LCMSIX(JPMIC,'MACROLIB ',1)
CALL LCMSIX(JPMIC,'GFF ',1)
CALL LCMLEN(JPMIC,'FINF_NUMBER ',NFDI,ITYLCM)
IF(NFDI+1.GT.MAXIFX) CALL XABORT('NAPCPO: MAXIFX OVERFLOW.')
IF(NFDI.GT.0) CALL LCMGET(JPMIC,'FINF_NUMBER ',FINF)
FINF(NFDI+1)=IFX
WRITE(HVECT,500) IFX
CALL LCMPUT(JPMIC,'FINF_NUMBER ',NFDI+1,1,FINF)
IF(IMPX.GE.10) THEN
WRITE(6,'(17H NAPCPO: MIXTURE=,I5,8H RECORD ,A8,1H=)') I,HVECT
DO IG=1,NG
WRITE(6,'(7H GROUP=,I5/(1X,1P,12E13.4))') IG,
1 (FLXP(IMIXP,IG),IMIXP=1,NMIXP)
ENDDO
ENDIF
CALL LCMPUT(JPMIC,HVECT,NMIXP*NG,2,FLXP)
CALL LCMSIX(JPMIC,'GFF ',2)
CALL LCMSIX(JPMIC,'*MAC*RES ',2)
ENDDO
DEALLOCATE(FT)
DEALLOCATE(FLAGMX)
DEALLOCATE(FLXD,FLXP)
DEALLOCATE(FXTD,FYTD)
DEALLOCATE(KEYFLX,MATCOD)
DEALLOCATE(MXD,MYD)
DEALLOCATE(XX,YY)
DEALLOCATE(MIXP)
DEALLOCATE(MXP,MYP)
DEALLOCATE(IXPD,IYPD)
DEALLOCATE(JDEBAR,JARBVA)
RETURN
*
500 FORMAT(5HFINF_,I3.3)
END
|
