diff options
| author | stainer_t <thomas.stainer@oecd-nea.org> | 2025-09-08 13:48:49 +0200 |
|---|---|---|
| committer | stainer_t <thomas.stainer@oecd-nea.org> | 2025-09-08 13:48:49 +0200 |
| commit | 7dfcc480ba1e19bd3232349fc733caef94034292 (patch) | |
| tree | 03ee104eb8846d5cc1a981d267687a729185d3f3 /Donjon/src/NAPCPO.f | |
Initial commit from Polytechnique Montreal
Diffstat (limited to 'Donjon/src/NAPCPO.f')
| -rw-r--r-- | Donjon/src/NAPCPO.f | 602 |
1 files changed, 602 insertions, 0 deletions
diff --git a/Donjon/src/NAPCPO.f b/Donjon/src/NAPCPO.f new file mode 100644 index 0000000..74efe96 --- /dev/null +++ b/Donjon/src/NAPCPO.f @@ -0,0 +1,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
|