*DECK NAPPPR SUBROUTINE NAPPPR(IPMAP,IPTRK,IPFLU,IPMTX,IPMAC,NSTATE) * *----------------------------------------------------------------------- * *Purpose: * Perform the Pin Power Reconstruction for core with * heterogeneous mixture * *Copyright: * Copyright (C) 2014 Ecole Polytechnique de Montreal. * *Author(s): * R. Chambon (EPM) and R. Nguyen Van Ho (URANUS) * *Parameters: input/output * IPMAP LCM object address of Map. * IPTRK LCM object address of Tracking. * IPFLU LCM object address of Flux. * IPMTX LCM object address of Matex. * IPMAC LCM object address of Macrolib of the fuel. * NSTATE length of the state vector * *----------------------------------------------------------------------- * USE GANLIB IMPLICIT NONE *---- * SUBROUTINE ARGUMENTS *---- INTEGER NSTATE TYPE(C_PTR) IPMAP,IPTRK,IPFLU,IPMTX,IPMAC *---- * LOCAL VARIABLES *---- INTEGER IOUT,NGPT REAL REPS PARAMETER (REPS=1.0E-4,IOUT=6,NGPT=2) TYPE(C_PTR) JPFLU,JPMAP,KPMAP INTEGER INDIC,NITMA,LENGTH,NBPIN CHARACTER TEXT*12 REAL FLOT DOUBLE PRECISION DFLOT INTEGER ISTATE(NSTATE),IMPX,IMETH INTEGER NXP,NYP,NXD,NYD,NZD,NAX,NAY, 1 NASS,NCOMB,NG,NASS2,NREG,NXM,NYM,NZM, 2 NXT,NYT,NZT,NXDA,NYDA,NZDA,NCH,NZASS,NPIN,IFX, 3 NUN,IEL,NMIX,NAMIX,NGFF CHARACTER LABEL*8 CHARACTER TFDINF*12 INTEGER I,J,K,IP,JP,I1,I2,J1,J2,K1,K2,IASS,ICH,IG,IM,JM,ID,JD,KM, 1 IAX,JAX,IGP,JGP,KGP,IMIX,IPIN,ICHX,IDIM,LC,L4,MAXKN,MKN,ITYLCM, 2 ITYPE REAL POW,FACT,POWTOT,POWASS,DX,DY,DZ,FPD,FQ,PMAX, 1 HOTPINPOW,PINPOW,FXY,VTOT REAL ZGKSIX(NGPT),ZGKSIY(NGPT),ZGKSIZ(NGPT),WGKSIX(NGPT), 1 WGKSIY(NGPT),WGKSIZ(NGPT),X(NGPT),Y(NGPT),Z(NGPT), 2 FLUGP(NGPT,NGPT,NGPT) REAL E(25) LOGICAL LSPX,LSPY,LSPZ,LCH,LPOW,LNOINT,LDEBUG *---- * ALLOCATABLE ARRAYS *---- CHARACTER*4, ALLOCATABLE, DIMENSION(:) :: NAMX,NAMY INTEGER, ALLOCATABLE, DIMENSION(:) :: NBAX,IBAX,BMIXP,AZONE, 1 ACOMB,KN INTEGER, ALLOCATABLE, DIMENSION(:,:) :: CODEA INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: KEYFLX,BMIX,MAT REAL, ALLOCATABLE, DIMENSION(:) :: MXP,MYP,MXD,MYD,MZD,MXM, 1 MYM,MZM,MXDA,MYDA,MZDA,FLXD,VOLM,FXYZ,PLINMAXZ,FXYASS, 2 FACTASS,PWASS,PWASS2 REAL, ALLOCATABLE, DIMENSION(:,:) :: FXTD,FYTD,BUNDPW REAL, ALLOCATABLE, DIMENSION(:,:,:) :: HFA,HFM,FDINFM, 1 FTINFM,AXPOW,VPIN REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: FLXDA,VOL REAL, ALLOCATABLE, DIMENSION(:,:,:,:,:) :: FLXP,HF,FDINF,FTINF * IMPX=0 FACT=1.0 LSPX=.FALSE. LSPY=.FALSE. LSPZ=.FALSE. LPOW=.FALSE. LNOINT=.FALSE. NZASS=0 NPIN=0 NBPIN=0 IFX=0 POW=1.0 FQ=0.0 FXY=0.0 HOTPINPOW=0.0 PINPOW=0.0 PMAX=0.0 VTOT=0.0 LDEBUG=.false. * Read mandatory keywords if(LDEBUG)write(6,*) 'NAPPPR begin debug' * [EDIT] PPR CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT) IF(INDIC.NE.3) CALL XABORT('NAPPPR: character data expected.') IF(TEXT.EQ.'EDIT') THEN CALL REDGET(INDIC,IMPX,FLOT,TEXT,DFLOT) IF(INDIC.NE.1) CALL XABORT('NAPPPR: inteGEr data expected.') CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT) IF(INDIC.NE.3) CALL XABORT('NAPPPR: character data expected.') ENDIF IF(TEXT.NE.'PPR') CALL XABORT('NAPPPR: ''PPR'' keyword '// 1 'expected.') !* NPIN ! CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT) ! IF(INDIC.NE.3) CALL XABORT('NAPPPR: character data expected.') ! IF(TEXT.NE.'NPIN') CALL XABORT('NAPPPR: ''NPIN'' keyword '// ! 1 'expected.') ! CALL REDGET(INDIC,NPIN,FLOT,TEXT,DFLOT) ! IF(INDIC.NE.1) CALL XABORT('NAPPPR: NPIN inteGEr expected.') ! NXP=NPIN ! NYP=NPIN * NZASS CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT) IF(INDIC.NE.3) CALL XABORT('NAPPPR: character data expected.') IF(TEXT.NE.'NZASS') CALL XABORT('NAPPPR: ''NZASS'' keyword '// 1 'expected.') CALL REDGET(INDIC,NZASS,FLOT,TEXT,DFLOT) IF(INDIC.NE.1) CALL XABORT('NAPPPR: NZASS inteGEr expected.') C* SPIN + SGAP C CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT) C IF(INDIC.NE.3) CALL XABORT('NAPPPR: character data expected.') C IF(TEXT.NE.'DIM') CALL XABORT('NAPPPR: ''NZASS'' keyword '// C 1 'expected.') C CALL REDGET(INDIC,NITMA,SPIN,TEXT,DFLOT) C IF(INDIC.NE.2) CALL XABORT('NAPPPR: SPIN real expected.') C CALL REDGET(INDIC,NITMA,SGAP,TEXT,DFLOT) C IF(INDIC.NE.2) CALL XABORT('NAPPPR: SGAP real expected.') C * GEt core GEometry description in matex IF(IMPX.GE.100)WRITE(6,*) 'debug:GEt matex info' CALL LCMGET(IPMTX,'STATE-VECTOR',ISTATE) NG=ISTATE(1) NXD=ISTATE(8) NYD=ISTATE(9) NZD=ISTATE(10) ALLOCATE(MXD(NXD+1),MYD(NYD+1),MZD(NZD+1)) CALL LCMGET(IPMTX,'MESHX',MXD) CALL LCMGET(IPMTX,'MESHY',MYD) CALL LCMGET(IPMTX,'MESHZ',MZD) * GEt KEYFLX IF(IMPX.GE.100)WRITE(6,*) 'debug:GEt track info' CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE) NREG=ISTATE(1) NUN=ISTATE(2) ITYPE=ISTATE(6) IEL=ISTATE(9) L4=ISTATE(11) ICHX=ISTATE(12) NXT=ISTATE(14) NYT=ISTATE(15) NZT=ISTATE(16) IDIM=1 IF((ITYPE.EQ.5).OR.(ITYPE.EQ.6).OR.(ITYPE.EQ.8)) IDIM=2 IF((ITYPE.EQ.7).OR.(ITYPE.EQ.9)) IDIM=3 IF((NXD.NE.NXT).OR.(NYD.NE.NYT).OR.(NZD.NE.NZT)) CALL XABORT 1 ('NAPPPR: dimension do not match between MATEX and TRACKING') ALLOCATE(KEYFLX(NXT,NYT,NZT),MAT(NXT,NYT,NZT)) CALL LCMGET(IPTRK,'KEYFLX',KEYFLX) CALL LCMGET(IPTRK,'MATCOD',MAT) ALLOCATE(FLXD(NUN)) * GEt assembly GEometry in map IF(IMPX.GE.100)WRITE(6,*) 'debug:GEt map info' CALL LCMGET(IPMAP,'STATE-VECTOR',ISTATE) NCH=ISTATE(2) NASS=ISTATE(14) NAX=ISTATE(15) NAY=ISTATE(16) ALLOCATE(AZONE(NCH)) ALLOCATE(NAMX(NAX),NAMY(NAY)) CALL LCMGET(IPMAP,'A-ZONE',AZONE) CALL LCMGTC(IPMAP,'AXNAME',4,NAX,NAMX) CALL LCMGTC(IPMAP,'AYNAME',4,NAY,NAMY) CALL LCMSIX(IPMAP,'GEOMAP',1) CALL LCMGET(IPMAP,'STATE-VECTOR',ISTATE) NXM=ISTATE(3) NYM=ISTATE(4) NZM=ISTATE(5) ALLOCATE(MXM(NXM+1),MYM(NYM+1),MZM(NZM+1)) ALLOCATE(NBAX(NAY),IBAX(NAY)) ALLOCATE(BMIX(NXM,NYM,NZM)) CALL LCMGET(IPMAP,'MESHX',MXM) CALL LCMGET(IPMAP,'MESHY',MYM) CALL LCMGET(IPMAP,'MESHZ',MZM) CALL LCMLEN(IPMAP,'A-NX',LENGTH,INDIC) IF(LENGTH.NE.NAY) CALL XABORT('NAPPPR: Number of assembly along' 1 //'Y direction do not match between MAP and embedded GEometry') CALL LCMGET(IPMAP,'A-NX',NBAX) CALL LCMGET(IPMAP,'A-IBX',IBAX) CALL LCMSIX(IPMAP,'GEOMAP',2) CALL LCMGET(IPMAP,'BMIX',BMIX) C* GEt data in pin by pin assembly GEometry C IF(IMPX.GE.100)WRITE(6,*) 'debug:GEt map pinBypin info' C CALL LCMGET(IPMPP,'STATE-VECTOR',ISTATE) C NCHP=ISTATE(2) C NASSP=ISTATE(14) C* total number of fuel bundles = tot. nb. of .XS C NXS=ISTATE(9) C IF(NASS.NE.NASSP)CALL XABORT('NAPPPR: number of assembly do not ' C 1 //'match between unfolded GEometries') C ALLOCATE(AZONEP(NCHP)) C CALL LCMGET(IPMPP,'A-ZONE',AZONEP) C CALL LCMSIX(IPMPP,'GEOMAP',1) C CALL LCMGET(IPMPP,'STATE-VECTOR',ISTATE) C NXMP=ISTATE(3) C NYMP=ISTATE(4) C NZMP=ISTATE(5) C CALL LCMSIX(IPMPP,'GEOMAP',2) C ALLOCATE(BMIXP(NXMP,NYMP,NZMP)) C CALL LCMGET(IPMPP,'BMIX',BMIXP) IF(IMPX.GE.5) THEN WRITE(6,*) 'MATEX dimension (het):',NXD,NYD,NZD WRITE(6,*) 'TRACKING dimension(het):',NXT,NYT,NZT WRITE(6,*) 'MAP dimension (het):',NXM,NYM,NZM ENDIF * Read remaining input file NCOMB=0 ALLOCATE(ACOMB(NASS)) IF(IMPX.GE.100)WRITE(6,*) 'debug: beg read input' IMETH=0 5 CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT) IF(TEXT.EQ.'METH') THEN CALL REDGET(INDIC,NITMA,FLOT,TEXT,DFLOT) IF(INDIC.NE.3) CALL XABORT('NAPPPR: character data expected.') IF(TEXT.EQ.'GPPR') THEN IMETH=1 CALL REDGET(INDIC,IFX,FLOT,TEXT,DFLOT) IF(INDIC.NE.1) CALL XABORT('NAPPPR: inteGEr data expected.') WRITE(TFDINF,500) IFX ELSE CALL XABORT('NAPPPR: '//TEXT//' is a wrong method keyword.') ENDIF GOTO 5 ELSEIF(TEXT.EQ.'POWER') THEN LPOW=.TRUE. CALL REDGET(INDIC,NITMA,POW,TEXT,DFLOT) IF(INDIC.NE.2) CALL XABORT('NAPPPR: POWER real expected.') GOTO 5 ELSEIF(TEXT.EQ.';') THEN GOTO 50 ELSE CALL XABORT('NAPPPR: '//TEXT//' is a wrong keyword.') ENDIF *----------------------------- 50 CONTINUE IF(IMPX.GE.100)WRITE(6,*) 'debug: computation begin' * Compute mesh X and Y for a pin-by-pin assembly CALL LCMGET(IPMAC,'STATE-VECTOR',ISTATE) NMIX=ISTATE(2) NAMIX=NMIX/NASS/NZASS IF(IMPX.GE.1) WRITE(6,*) 'Number of Mix per assembly per plane'// 1 ' NAMIX = ',NAMIX NGFF=ISTATE(16) IF(NGFF.EQ.0) CALL XABORT('NAPPPR: NGFF.NE.0 expected.') CALL LCMSIX(IPMAC,'GFF',1) CALL LCMSIX(IPMAC,'GFF-GEOM',1) CALL LCMGET(IPMAC,'STATE-VECTOR',ISTATE) NXP=ISTATE(3) NYP=ISTATE(4) NPIN=NXP ALLOCATE(MXP(NXP+1),MYP(NYP+1)) CALL LCMGET(IPMAC,'MESHX',MXP) CALL LCMGET(IPMAC,'MESHY',MYP) DO I=2,NXP+1 MXP(I)=MXP(I)-MXP(1) ENDDO MXP(1)=0.0 DO I=2,NYP+1 MYP(I)=MYP(I)-MYP(1) ENDDO MYP(1)=0.0 CALL LCMSIX(IPMAC,'GFF-GEOM',2) CALL LCMSIX(IPMAC,'GFF',2) * Compute IX-,IX+,IY-,IY+,IZ-,IZ+ for each assembly in core GEometry IF(IMPX.GE.100) WRITE(6,*) 'debug PPR:IX-,IX+,IY-,IY+,IZ-,IZ+' ALLOCATE(CODEA(NASS,6)) CODEA(:NASS,:6)=0 ICH=0 I1=0 I2=0 NASS2=0 DO IASS=1,NASS CODEA(IASS,1)=NXD+1 CODEA(IASS,2)=0 CODEA(IASS,3)=NYD+1 CODEA(IASS,4)=0 CODEA(IASS,5)=NZD+1 CODEA(IASS,6)=0 ENDDO DO 150 JM=1,NYM DO 130 IM=1,NXM LCH=.TRUE. IASS=0 DO 100 KM=1,NZM IF(BMIX(IM,JM,KM).NE.0) THEN IF(LCH) THEN ICH=ICH+1 LCH=.FALSE. IASS=AZONE(ICH) NASS2=MAX(NASS2,IASS) DO I=1,NXD+1 IF(MXD(I).EQ.MXM(IM)) I1=I IF(MXD(I).EQ.MXM(IM+1)) I2=I ENDDO CODEA(IASS,1)=MIN(I1,CODEA(IASS,1)) CODEA(IASS,2)=MAX(I2,CODEA(IASS,2)) DO I=1,NYD+1 IF(MYD(I).EQ.MYM(JM)) I1=I IF(MYD(I).EQ.MYM(JM+1)) I2=I ENDDO CODEA(IASS,3)=MIN(I1,CODEA(IASS,3)) CODEA(IASS,4)=MAX(I2,CODEA(IASS,4)) DO I=1,NZD+1 IF(MZD(I).EQ.MZM(KM)) I1=I IF(MZD(I).EQ.MZM(KM+1)) I2=I ENDDO CODEA(IASS,5)=MIN(I1,CODEA(IASS,5)) CODEA(IASS,6)=MAX(I2,CODEA(IASS,6)) ELSE DO I=2,NZD+1 IF(MZD(I).EQ.MZM(KM+1)) I2=I ENDDO CODEA(IASS,6)=MAX(I2,CODEA(IASS,6)) ENDIF ENDIF 100 CONTINUE 130 CONTINUE 150 CONTINUE IF(IMPX.GE.10) THEN WRITE(6,*) 'Position of the assemblies in the core' WRITE(6,*) 'IX-,IX+,IY-,IY+,IZ-,IZ+' do iass=1,nass WRITE(6,*) 'Assembly #',iass,':',(CODEA(iass,i),i=1,6) ENDDO ENDIF IF(NASS2.NE.NASS)CALL XABORT('NAPPPR: number of assembly do not' 1 //' match: NASS2.NE.NASS') * For all assembly perform PPR ALLOCATE(FLXP(NXP,NYP,NZASS,NG,NASS)) ALLOCATE(AXPOW(NXP,NYP,NASS)) ALLOCATE(VPIN(NXP,NYP,NASS)) ALLOCATE(FXYASS(NASS)) ALLOCATE(FACTASS(NASS)) ALLOCATE(PWASS(NASS),PWASS2(NASS)) ALLOCATE(BUNDPW(NASS,NZASS)) IF(.NOT.LPOW) THEN CALL LCMGET(IPMAP,'BUND-PW',BUNDPW) ENDIF DO IASS=1,NASS PWASS(IASS)=0.0 DO IP=1,NXP DO JP=1,NYP AXPOW(IP,JP,IASS)=0.0 VPIN(IP,JP,IASS)=0.0 ENDDO ENDDO FXYASS(IASS)=0.0 IF(.NOT.LPOW) THEN DO K=1,NZASS PWASS(IASS)=PWASS(IASS)+BUNDPW(IASS,K) ENDDO ENDIF ENDDO DO IASS=1,NASS * GEt flux at core GEometry level for assembly only I1=CODEA(IASS,1) I2=CODEA(IASS,2) J1=CODEA(IASS,3) J2=CODEA(IASS,4) K1=CODEA(IASS,5) K2=CODEA(IASS,6) NXDA=I2-I1 NYDA=J2-J1 NZDA=K2-K1 ALLOCATE(FLXDA(NXDA,NYDA,NZDA,NG)) FLXDA(:NXDA,:NYDA,:NZDA,:NG)=0.0 IF(NZDA.NE.NZASS) CALL XABORT('NAPPPR: incoherent number of mesh' 1 //' in Z direction for an assembly: NZDA.NE.NZASS') JPFLU=LCMGID(IPFLU,'FLUX') IF((LNOINT).OR.(IMPX.GE.0)) THEN DO IG=1,NG CALL LCMGDL(JPFLU,IG,FLXD) DO I=I1,I2-1 DO J=J1,J2-1 DO K=K1,K2-1 FLXDA(I-I1+1,J-J1+1,K-K1+1,IG)=FLXD(KEYFLX(I,J,K)) ENDDO C end K ENDDO C end J ENDDO C end I ENDDO C end IG ENDIF ALLOCATE(MXDA(NXDA+1),MYDA(NYDA+1),MZDA(NZDA+1)) DO I=I1,I2 MXDA(I-I1+1)=MXD(I)-MXD(I1)+MXP(1) ENDDO IF(ABS(MXDA(NXDA+1)-MXDA(1)-MXP(NXP+1)+MXP(1)).GT.0.0001) THEN WRITE(6,*) 'Assembly Transport and Core meshX do not match:'// 1 'Transport=',MXP(NXP+1)-MXP(1),'Core=',MXDA(NXDA+1)-MXDA(1) CALL XABORT('Sizes do not match') ENDIF DO J=J1,J2 MYDA(J-J1+1)=MYD(J)-MYD(J1)+MYP(1) ENDDO IF(ABS(MYDA(NYDA+1)-MYDA(1)-MYP(NYP+1)+MYP(1)).GT.0.0001) THEN WRITE(6,*) 'Assembly Transport and Core meshY do not match:'// 1 'Transport=',MYP(NYP+1)-MYP(1),'Core=',MYDA(NYDA+1)-MYDA(1) CALL XABORT('Sizes do not match') ENDIF DO K=K1,K2 MZDA(K-K1+1)=MZD(K) ENDDO IF(IMPX.GE.10) THEN WRITE(6,*) 'Coarse Flux and mesh at assembly level' WRITE(6,*) 'Mesh X:',(MXDA(I),I=1,NXDA+1) WRITE(6,*) 'Mesh Y:',(MYDA(I),I=1,NYDA+1) WRITE(6,*) 'Mesh Z:',(MZDA(I),I=1,NZDA+1) WRITE(6,*) 'Flux:' DO IG=1,NG WRITE(6,*) 'Group #',IG DO K=1,NZDA WRITE(6,*) 'Plan #',K DO J=1,NYDA WRITE(6,*) (FLXDA(I,J,K,IG),I=1,NXDA) ENDDO ENDDO ENDDO ENDIF * project flux at assembly level ALLOCATE(FXTD(NXP,NXDA),FYTD(NYP,NYDA)) FXTD(:NXP,:NXDA)=0.0 FYTD(:NYP,:NYDA)=0.0 * compute fraction of the transport volumes occupied by diffusion volumes CALL NAPFTD(NXP,MXP,NXDA,MXDA,FXTD) CALL NAPFTD(NYP,MYP,NYDA,MYDA,FYTD) ! DO IP=1,NXP ! DXP=MXP(IP+1)-MXP(IP) ! DO ID=1,NXDA ! IF((MXDA(ID).LE.MXP(IP)).AND.(MXDA(ID+1).GE.MXP(IP+1))) THEN ! FXTD(IP,ID)=1.0 ! ELSEIF ((MXDA(ID).LE.MXP(IP)).AND.(MXDA(ID+1).GT.MXP(IP))) THEN ! FXTD(IP,ID)=(MXDA(ID+1)-MXP(IP))/DXP ! ELSEIF ((MXDA(ID).GE.MXP(IP)).AND. ! 1 (MXDA(ID+1).LE.MXP(IP+1))) THEN ! FXTD(IP,ID)=(MXDA(ID+1)-MXDA(ID))/DXP ! ELSEIF ((MXDA(ID).LT.MXP(IP+1)).AND. ! 1 (MXDA(ID+1).GE.MXP(IP+1))) THEN ! FXTD(IP,ID)=(MXP(IP+1)-MXDA(ID))/DXP ! ENDIF ! ENDDO ! ENDDO * ! DO IP=1,NYP ! DYP=MYP(IP+1)-MYP(IP) ! DO ID=1,NYDA ! IF((MYDA(ID).LE.MYP(IP)).AND.(MYDA(ID+1).GE.MYP(IP+1))) THEN ! FYTD(IP,ID)=1.0 ! ELSEIF ((MYDA(ID).LE.MYP(IP)).AND.(MYDA(ID+1).GT.MYP(IP))) THEN ! FYTD(IP,ID)=(MYDA(ID+1)-MYP(IP))/DYP ! ELSEIF ((MYDA(ID).GE.MYP(IP)).AND. ! 1 (MYDA(ID+1).LE.MYP(IP+1))) THEN ! FYTD(IP,ID)=(MYDA(ID+1)-MYDA(ID))/DYP ! ELSEIF ((MYDA(ID).LT.MYP(IP+1)).AND. ! 1 (MYDA(ID+1).GE.MYP(IP+1))) THEN ! FYTD(IP,ID)=(MYP(IP+1)-MYDA(ID))/DYP ! ENDIF ! ENDDO ! ENDDO ! adds up all fluxes if(LDEBUG)write(6,*)'NXP,NYP',NXP,NYP DO IG=1,NG IF(.NOT.LNOINT) CALL LCMGDL(JPFLU,IG,FLXD) DO K=1,NZASS DO IP=1,NXP DO JP=1,NYP FLXP(IP,JP,K,IG,IASS)=0.0 DO ID=1,NXDA DO JD=1,NYDA IF(LNOINT) THEN * No interpolation: use averaGE flux FLXP(IP,JP,K,IG,IASS)=FLXP(IP,JP,K,IG,IASS) 1 +FLXDA(ID,JD,K,IG)*FXTD(IP,ID)*FYTD(JP,JD) * Interpolate flux with polynomial representation * (only if pin and macro region have a non-nul intersection) ELSEIF(FXTD(IP,ID)*FYTD(JP,JD).NE.0.0) THEN * indent removed * compute gauss points and weights CALL ALGPT(NGPT,MAX(MXP(IP),MXDA(ID)),MIN(MXP(IP+1),MXDA(ID+1)), 1 ZGKSIX,WGKSIX) DX=MIN(MXP(IP+1),MXDA(ID+1))-MAX(MXP(IP),MXDA(ID)) CALL ALGPT(NGPT,MAX(MYP(JP),MYDA(JD)),MIN(MYP(JP+1),MYDA(JD+1)), 1 ZGKSIY,WGKSIY) DY=MIN(MYP(JP+1),MYDA(JD+1))-MAX(MYP(JP),MYDA(JD)) CALL ALGPT(NGPT,MZDA(K),MZDA(K+1),ZGKSIZ,WGKSIZ) DZ=MZDA(K+1)-MZDA(K) IF(IMPX.GE.10) 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 WRITE(6,*) 'Gauss point ZGWG:',(ZGKSIZ(I),I=1,NGPT), 1 (WGKSIZ(I),I=1,NGPT),'DZ',DZ ENDIF * interpolate flux FPD=0.0 DO IGP=1,NGPT X(IGP)=MXD(I1)-MXP(1)+ZGKSIX(IGP) ENDDO DO JGP=1,NGPT Y(JGP)=MYD(J1)-MYP(1)+ZGKSIY(JGP) ENDDO DO KGP=1,NGPT Z(KGP)=ZGKSIZ(KGP) ENDDO IF(IMPX.GE.10) then WRITE(6,*) 'Gauss point X:',(X(I),I=1,NGPT) WRITE(6,*) 'Gauss point Y:',(Y(I),I=1,NGPT) WRITE(6,*) 'Gauss point Z:',(Z(I),I=1,NGPT) ENDIF IF(ICHX.EQ.1) THEN * Variational collocation method CALL LCMLEN(IPTRK,'KN',MAXKN,ITYLCM) MKN=MAXKN/(NXD*NYD*NZD) 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 VALUE2(LC,MKN,NXD,NYD,NZD,L4,X,Y,Z,MXD,MYD,MZD, 1 FLXD,MAT,KN,NGPT,NGPT,NGPT,E,FLUGP) DEALLOCATE(KN) ELSE IF(ICHX.EQ.2) THEN * Raviart-Thomas finite element method CALL VALUE4(IEL,NUN,NXD,NYD,NZD,X,Y,Z,MXD,MYD,MZD, 1 FLXD,MAT,KEYFLX,NGPT,NGPT,NGPT,FLUGP) ELSE IF(ICHX.EQ.3) THEN * Nodal collocation method (MCFD) CALL VALUE1(IDIM,NXD,NYD,NZD,L4,X,Y,Z,MXD,MYD,MZD, 1 FLXD,MAT,IEL,NGPT,NGPT,NGPT,FLUGP) ELSE CALL XABORT('NAPPPR: INTERPOLATION NOT IMPLEMENTED.') ENDIF IF(IMPX.GE.10) then WRITE(6,*) 'Gauss flux values:' DO KGP=1,NGPT WRITE(6,*) 'KGP=:',KGP DO JGP=1,NGPT WRITE(6,*) (FLUGP(IGP,JGP,KGP),IGP=1,NGPT) ENDDO ENDDO ENDIF * integrate flux (gauss method) DO IGP=1,NGPT DO JGP=1,NGPT DO KGP=1,NGPT FPD=FPD+FLUGP(IGP,JGP,KGP)*WGKSIX(IGP)*WGKSIY(JGP)*WGKSIZ(KGP) ENDDO ENDDO ENDDO * GEt averaGE flux FPD=FPD/DX/DY/DZ if(LDEBUG)write(6,*)'FLXP,FPD,FXTD,FYTD',FLXP(IP,JP,K,IG,IASS), 1 FPD,FXTD(IP,ID),FYTD(JP,JD) FLXP(IP,JP,K,IG,IASS)=FLXP(IP,JP,K,IG,IASS) 1 +FPD*FXTD(IP,ID)*FYTD(JP,JD) if(LDEBUG)write(6,*)'FLXP after',FLXP(IP,JP,K,IG,IASS) * indent back ENDIF ENDDO ENDDO ENDDO ENDDO ENDDO ENDDO * DEALLOCATE(FXTD,FYTD) DEALLOCATE(MXDA,MYDA,MZDA) DEALLOCATE(FLXDA) IF(IMPX.GE.100)WRITE(6,*) 'debug PPR:projection flux for one ' 1 //'assem end' ! end of DO IASS=1,NASS ENDDO IF(IMPX.GE.100)WRITE(6,*) 'debug PPR:projection flux for all ' 1 //'assem end' * GPPR IF(IMETH.EQ.1) THEN IF(IMPX.GE.100)WRITE(6,*) 'debug PPR:',TFDINF * GEt Volume, phi^t,inf_p and phi^d,inf_m,p from macrolib of fuel * Note: if homoGEneous (normal PPR), m=1 ALLOCATE(VOLM(NGFF),HFM(NMIX,NGFF,NG), 1 FTINFM(NMIX,NGFF,NG),FDINFM(NMIX,NGFF,NG)) ALLOCATE(VOL(NPIN,NPIN,NZASS,NASS)) ALLOCATE(HF(NPIN,NPIN,NZASS,NG,NASS)) ALLOCATE(FTINF(NPIN,NPIN,NZASS,NG,NASS)) ALLOCATE(FDINF(NPIN,NPIN,NZASS,NG,NASS)) ALLOCATE(BMIXP(NPIN*NPIN)) VOL(:NPIN,:NPIN,:NZASS,:NASS)=0.0 HF(:NPIN,:NPIN,:NZASS,:NG,:NASS)=0.0 FTINF(:NPIN,:NPIN,:NZASS,:NG,:NASS)=0.0 FDINF(:NPIN,:NPIN,:NZASS,:NG,:NASS)=0.0 if(LDEBUG) call LCMLIB(IPMAC) CALL LCMSIX(IPMAC,'GFF',1) if(LDEBUG) call LCMLIB(IPMAC) CALL LCMGET(IPMAC,'VOLUME',VOLM) CALL LCMGET(IPMAC,'H-FACTOR',HFM) CALL LCMGET(IPMAC,'NWT0',FTINFM) CALL LCMGET(IPMAC,TFDINF,FDINFM) CALL LCMSIX(IPMAC,'GFF-GEOM',1) CALL LCMGET(IPMAC,'MIX',BMIXP) CALL LCMSIX(IPMAC,'GFF-GEOM',2) CALL LCMSIX(IPMAC,'GFF',2) DO IG=1,NG DO IASS=1,NASS K1=CODEA(IASS,5) DO K=1,NZASS ! NAMIX = 1 for homogeneous assembly ! > 1 for heterogeneous assembly ! Note that all values of HFM are identical ! for all the mix in a specific assembly IMIX=(IASS-1+(K-1)*NASS)*NAMIX+1 DO J=1,NPIN DO I=1,NPIN IPIN=I+(J-1)*NPIN HF(I,J,K,IG,IASS)=HFM(IMIX,BMIXP(IPIN),IG) FTINF(I,J,K,IG,IASS)=FTINFM(IMIX,BMIXP(IPIN),IG) FDINF(I,J,K,IG,IASS)=FDINFM(IMIX,BMIXP(IPIN),IG) VOL(I,J,K,IASS)=(MXP(I+1)-MXP(I))*(MYP(J+1)-MYP(J)) 3 *(MZM(K1+K)-MZM(K1+K-1)) ENDDO ENDDO ENDDO ! end of DO IASS=1,NASS ENDDO ! end of DO IG=1,NG ENDDO IF(IMPX.GE.6) then DO iass=1,nass WRITE(6,*) 'XS for assembly #',IASS DO k=1,nzass WRITE(6,*) 'Plane #',K DO ig=1,ng WRITE(6,*) 'group #',ig WRITE(6,*) 'HF #' DO j=1,npin WRITE(6,*) (HF(I,J,K,ig,iass),I=1,NPIN) ENDDO WRITE(6,*) 'FTINF #' DO j=1,npin WRITE(6,*) (FTINF(I,J,K,ig,iass),I=1,NPIN) ENDDO WRITE(6,*) 'FLXP #' DO j=1,npin WRITE(6,*) (FLXP(I,J,K,ig,iass),I=1,NPIN) ENDDO WRITE(6,*) 'FDINF #' DO j=1,npin WRITE(6,*) (FDINF(I,J,K,ig,iass),I=1,NPIN) ENDDO ! end of do ig=1,ng ENDDO WRITE(6,*) 'VOL #' DO j=1,npin WRITE(6,*) (VOL(I,J,K,iass),I=1,NPIN) ENDDO ENDDO ENDDO ENDIF * Print and save reaction rates IF(IMPX.GE.100)WRITE(6,*) 'debug: Print and save reaction rates' POWTOT=0.0 DO IASS=1,NASS PWASS2(IASS)=0.0 K1=CODEA(IASS,5) DO K=1,NZASS DO J=1,NPIN DO I=1,NPIN VTOT=VTOT+VOL(I,J,K,IASS) DO IG=1,NG POWTOT=POWTOT+HF(I,J,K,IG,IASS)*FTINF(I,J,K,IG,IASS) 1 *FLXP(I,J,K,IG,IASS)/FDINF(I,J,K,IG,IASS) 2 *VOL(I,J,K,IASS) PWASS2(IASS)=PWASS2(IASS) 1 +HF(I,J,K,IG,IASS)*FTINF(I,J,K,IG,IASS) 1 *FLXP(I,J,K,IG,IASS)/FDINF(I,J,K,IG,IASS) 2 *VOL(I,J,K,IASS) ENDDO ENDDO ENDDO ENDDO ENDDO IF(IMPX.GE.2) WRITE(6,*) 'POWTOT:',POWTOT IF(LPOW) THEN DO IASS=1,NASS FACTASS(IASS)=POW/POWTOT ENDDO ELSE DO IASS=1,NASS FACTASS(IASS)=PWASS(IASS)/PWASS2(IASS) ENDDO ENDIF IF(IMPX.GE.2) WRITE(6,*) 'FACTASS:',(FACTASS(I),I=1,NASS) ALLOCATE(HFA(NPIN,NPIN,NZASS)) ALLOCATE(FXYZ(NZASS)) ALLOCATE(PLINMAXZ(NZASS)) DO K=1,NZASS FXYZ(K)=0.0 PLINMAXZ(K)=0.0 ENDDO JPMAP=LCMLID(IPMAP,'ASSEMBLY',NASS) IAX=0 JAX=1 DO IASS=1,NASS K1=CODEA(IASS,5) IAX=IAX+1 IF(IAX.GT.NBAX(JAX)) THEN IAX=1 JAX=JAX+1 ENDIF WRITE(LABEL,'(A4,A4)') NAMY(JAX),NAMX(IBAX(JAX)+IAX-1) IF(IMPX.GE.5) THEN WRITE(6,*) 'Reaction rates for assembly #',IASS,' Label:', 1 LABEL ENDIF DO K=1,NZASS IF(IMPX.GE.5) WRITE(6,*) 'Plane #',K DO J=1,NPIN DO I=1,NPIN HFA(I,J,K)=0.0 DO IG=1,NG HFA(I,J,K)=HFA(I,J,K)+HF(I,J,K,IG,IASS)*FTINF(I,J,K,IG,IASS) 1 *FLXP(I,J,K,IG,IASS)/FDINF(I,J,K,IG,IASS) 2 *FACTASS(IASS) 2 *VOL(I,J,K,IASS) IF((PLINMAXZ(K)*(MZM(K1+K)-MZM(K1+K-1))).LT.HFA(I,J,K)) THEN PLINMAXZ(K)=HFA(I,J,K)/(MZM(K1+K)-MZM(K1+K-1)) ENDIF ! end of DO IG=1,NG ENDDO ! end of I=1,NPIN ENDDO IF(IMPX.GE.5) WRITE(6,*) (HFA(I,J,K),I=1,NPIN) ! end of J=1,NPIN ENDDO ! end of DO K=1,NZASS ENDDO * KPMAP=LCMDIL(JPMAP,IASS) CALL LCMPTC(KPMAP,'LABEL',8,LABEL) CALL LCMPUT(KPMAP,'PIN-POWER',NPIN*NPIN*NZASS,2,HFA) CALL LCMPUT(KPMAP,'FLUX',NPIN*NPIN*NZASS*NG,2, 1 FLXP(1,1,1,1,IASS)) POWASS=0.0 DO I=1,NPIN DO J=1,NPIN DO K=1,NZASS POWASS=POWASS+HFA(I,J,K)!power of the assembly iass VPIN(I,J,IASS)=VPIN(I,J,IASS)+VOL(I,J,K,IASS) !pin volume ENDDO ENDDO ENDDO DO I=1,NPIN DO J=1,NPIN DO K=1,NZASS AXPOW(I,J,IASS)=HFA(I,J,K) 2 +AXPOW(I,J,IASS) !AXPOW:axially integrated pin power per pin !normalized to the pin mean power IF(PMAX.LT.HFA(I,J,K)) THEN PMAX=HFA(I,J,K)!maximal 3D local power ENDIF ENDDO AXPOW(I,J,IASS)=AXPOW(I,J,IASS)/(POWASS/NPIN/NPIN) ENDDO ENDDO * IF(IMPX.GE.2) WRITE(6,*) 'Power of assembly #',IASS,":",POWASS DO I=1,NPIN DO J=1,NPIN IF(IMPX.GE.2) THEN WRITE(6,*) 'AXPOW for assembly #',IASS NBPIN=NBPIN+1 WRITE(6,*) 'ASS:',IASS,'PIN #',NBPIN,":",AXPOW(I,J,IASS) ENDIF PINPOW=AXPOW(I,J,IASS)*VPIN(I,J,IASS) IF(HOTPINPOW.LT.PINPOW) THEN HOTPINPOW=PINPOW !power of the hot pin normalized to the pin mean power ENDIF IF(FXYASS(IASS).LT.AXPOW(I,J,IASS)) THEN FXYASS(IASS)=AXPOW(I,J,IASS) ENDIF ENDDO ENDDO NBPIN=0 * IF(IMPX.GE.1) THEN WRITE(6,*) 'Fxy for assembly #',IASS,":",FXYASS(IASS) ENDIF CALL LCMPUT(KPMAP,'ASS-POWER',1,2,POWASS) ! end of DO IASS=1,NASS ENDDO ! end of IF(IMETH.EQ.1) THEN ENDIF * FQ=PMAX FXY=HOTPINPOW * IF(IMPX.GE.0) THEN WRITE(6,*) 'FQ=',FQ WRITE(6,*) 'FXY=',FXY DO K=1,NZASS FXYZ(K)=PLINMAXZ(K) IF(IMPX.GE.0) WRITE(6,*) 'Plane #',K,'---> FXYZ(Z)=',FXYZ(K) ENDDO ENDIF CALL LCMPUT(IPMAP,'FQ',1,2,FQ) CALL LCMPUT(IPMAP,'FXY',1,2,FXY) CALL LCMPUT(IPMAP,'FXYZ',NZASS,2,FXYZ) CALL LCMPUT(IPMAP,'FXYASS',IASS,2,FXYASS) CALL LCMGET(IPMAP,'STATE-VECTOR',ISTATE) ISTATE(17)=NZASS CALL LCMPUT(IPMAP,'STATE-VECTOR',NSTATE,1,ISTATE) ! IF(IMPX.GE.100)WRITE(6,*) 'debug: beging deallacate' DEALLOCATE(FLXP,FLXD) DEALLOCATE(MXP,MYP) DEALLOCATE(CODEA) IF(IMETH.EQ.1) THEN DEALLOCATE(VOLM,HFM,FTINFM,FDINFM) DEALLOCATE(VOL,HF,FTINF,FDINF,AXPOW,FXYASS) DEALLOCATE(HFA,FXYZ,PLINMAXZ,VPIN) DEALLOCATE(FACTASS,PWASS,PWASS2) DEALLOCATE(BUNDPW) ENDIF DEALLOCATE(ACOMB) DEALLOCATE(AZONE) DEALLOCATE(NAMX,NAMY) DEALLOCATE(BMIX,BMIXP) DEALLOCATE(MXM,MYM,MZM) DEALLOCATE(NBAX,IBAX) DEALLOCATE(KEYFLX,MAT) DEALLOCATE(MXD,MYD,MZD) RETURN 500 FORMAT(5HFINF_,I3.3,4H ) END