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/NAPPPR.f | |
Initial commit from Polytechnique Montreal
Diffstat (limited to 'Donjon/src/NAPPPR.f')
| -rw-r--r-- | Donjon/src/NAPPPR.f | 866 |
1 files changed, 866 insertions, 0 deletions
diff --git a/Donjon/src/NAPPPR.f b/Donjon/src/NAPPPR.f new file mode 100644 index 0000000..d1f891c --- /dev/null +++ b/Donjon/src/NAPPPR.f @@ -0,0 +1,866 @@ +*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
|