*DECK FLU2DR SUBROUTINE FLU2DR(IPRT,IPMACR,IPFLUX,IPSYS,IPTRK,IPFLUP,IPSOU, 1 IGPT,IFTRAK,CXDOOR,TITLE,NUNKNO,NREG,NSOUT,NANIS,NLF,NLIN,NFUNL, 2 NGRP,NMAT,NIFIS,LFORW,LEAKSW,MAXINR,EPSINR,MAXOUT,EPSUNK,EPSOUT, 3 NCPTL,NCPTA,ITYPEC,IPHASE,ITPIJ,ILEAK,OPTION,REFKEF,MATCOD, 4 KEYFLX,VOL,XSTRC,XSDIA,XSNUF,XSCHI,LREBAL,INITFL,NMERG,IMERG) * *----------------------------------------------------------------------- * *Purpose: * Fixed source problem or inverse power method for K-effective or * buckling iteration. Perform thermal iterations. * *Copyright: * Copyright (C) 2002 Ecole Polytechnique de Montreal * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public * License as published by the Free Software Foundation; either * version 2.1 of the License, or (at your option) any later version * *Author(s): R. Roy * *Parameters: input * IPRT print flag. * IPMACR pointer to the macrolib LCM object. * IPFLUX pointer to the flux LCM object. * IPSYS pointer to the system LCM object. * IPTRK pointer to the tracking LCM object. * IPFLUP pointer to the unperturbed flux LCM object (if ITYPEC=1). * IPSOU pointer to the fixed source LCM object (if ITYPEC=0 or 1). * IGPT index of the fixed source eigenvalue problem to solve. * IFTRAK tracking file unit number. * CXDOOR character name of the flux solution door. * TITLE title. * NUNKNO number of unknowns per energy group including spherical * harmonic terms, interface currents and fundamental * currents. * NREG number of regions. * NSOUT number of outer surfaces. * NANIS maximum cross section Legendre order in object IPMACR. * NLF number of Legendre orders for the flux. * NLIN number of polynomial components in flux spatial expansion. * NFUNL number of spherical harmonics components. * NGRP number of energy groups. * NMAT number of mixtures in the macrolib. * NIFIS number of fissile isotopes. * LFORW flag set to .false. to solve an adjoint problem. * LEAKSW leakage flag (=.true. if leakage is present on the outer * surface). * MAXINR maximum number of thermal iterations. * EPSINR thermal iterations epsilon. * MAXOUT maximum number of outer iterations. * EPSUNK outer iterations eigenvector epsilon. * EPSOUT outer iterations eigenvalue epsilon. * NCPTL number of free iterations in an acceleration cycle. * NCPTA number of accelerated iterations in an acceleration cycle. * ITYPEC type of flux evaluation: * =-2 Fourier analysis; * =-1 skip the flux calculation; * =0 fixed sources; * =1 fixed source eigenvalue problem (GPT type); * =2 fission sources/k effective convergence; * =3 fission sources/k effective convergence/ * db2 buckling evaluation; * =4 fission sources/db2 buckling convergence; * =5 b2 sources/db2 buckling convergence. * IPHASE type of flux solution door (1 for asm 2 for pij). * ITPIJ type of cp available: * =1 scatt mod pij (wij); * =2 stand. pij; * =3 scatt mod pij+pijk (wij,wijk); * =4 standard pij+pijk. * ILEAK method used to include DB2 effect: * =1 the scattering modified cp matrix is multiplied by PNLR; * =2 the reduced cp matrix is multiplied by PNL; * =3 sigs0-db2 approximation; * =4 albedo approximation; * =5 Todorova-type isotropic streaming model; * =6 Ecco-type isotropic streaming model; * >6 Tibere type anisotropic streaming model. * OPTION type of leakage coefficients: * 'LKRD' (recover leakage coefficients in Macrolib); * 'RHS' (recover leakage coefficients in RHS flux object); * 'B0' (B-0), 'P0' (P-0), 'B1' (B-1), * 'P1' (P-1), 'B0TR' (B-0 with transport correction) or 'P0TR' * (P-0 with transport correction). * REFKEF target effective multiplication factor. * MATCOD mixture indices. * KEYFLX index of L-th order flux components in unknown vector. * VOL volumes. * XSTRC transport-corrected macroscopic total cross sections. * XSDIA transport-corrected macroscopic within-group scattering cross * sections. * XSNUF nu*macroscopic fission cross sections. * XSCHI fission spectrum. * LREBAL thermal iteration rebalancing flag (=.true. if thermal * rebalancing required). * INITFL flux initialization flag (=0/1/2: uniform flux/LCM/DSA). * NMERG number of leakage zones. * IMERG leakage zone index in each material mixture zone. * *----------------------------------------------------------------------- * *---- * INTERNAL PARAMETERS: * SYBILF : SYBIL FLUX SOLUTION DOOR EXT ROUTINE * TRFICF : DEFAULT CP FLUX SOLUTION DOOR EXT ROUTINE * BIVAF : DEFAULT 2D DIFFUSION FLUX SOLUTION DOOR EXT ROUTINE * TRIVAF : DEFAULT 3D DIFFUSION FLUX SOLUTION DOOR EXT ROUTINE * PNF : DEFAULT PN/SPN FLUX SOLUTION DOOR EXT ROUTINE * SNF : DEFAULT SN FLUX SOLUTION DOOR EXT ROUTINE *---- * USE GANLIB USE DOORS_MOD *---- * SUBROUTINE ARGUMENTS *---- TYPE(C_PTR) IPMACR,IPFLUX,IPSYS,IPTRK,IPFLUP,IPSOU INTEGER IPRT,IGPT,IFTRAK,NUNKNO,NREG,NSOUT,NANIS,NLF,NLIN,NFUNL, 1 NGRP,NMAT,NIFIS,MAXINR,MAXOUT,NCPTL,NCPTA,ITYPEC,IPHASE,ITPIJ, 2 ILEAK,MATCOD(NREG),KEYFLX(NREG,NLIN,NFUNL),INITFL,NMERG, 3 IMERG(NMAT) REAL EPSINR,EPSUNK,EPSOUT,VOL(NREG),XSTRC(0:NMAT,NGRP), 1 XSDIA(0:NMAT,0:NANIS,NGRP),XSNUF(0:NMAT,NIFIS,NGRP), 2 XSCHI(0:NMAT,NIFIS,NGRP) CHARACTER CXDOOR*12,TITLE*72,OPTION*4,HLEAK*6 LOGICAL LFORW,LEAKSW,LREBAL,CFLI,CEXE DOUBLE PRECISION REFKEF *---- * LOCAL VARIABLES *---- PARAMETER(NSTATE=40,PI=3.141592654) TYPE(C_PTR) IPREB,J1,JPSOU,JPFLUX,JPMACR,KPMACR,JPSYS,KPSYS,IPSTR, 1 JPSTR,KPSTR,JPFLUP1,JPFLUP2,JPSOUR INTEGER JPAR(NSTATE),KEYSPN(NREG) CHARACTER CAN(0:19)*2,MESSIN*8,MESSOU*5,HTYPE(0:5)*4 INTEGER INDD(3) DOUBLE PRECISION AKEEP(8),FISOUR,OLDBIL,AKEFF,AKEFFO,AFLNOR, 1 BFLNOR,DDELN1,DDELD1,PROD,FLXIN LOGICAL LSCAL,LEXAC,REBFLG REAL ALBEDO(6),FLUXC(NREG),B2(4) * ************************************************************************ * * * ICHAR : COUNTER FOR NUM. OF OUTER ITERATIONS * * ICTOT : TOTAL NUMBER OF FLUX CALCULATIONS * * * ************************************************************************ * *---- * ALLOCATABLE ARRAYS *---- INTEGER, ALLOCATABLE, DIMENSION(:) :: IJJ,NJJ,IPOS,NPSYS,KEYCUR, 1 MATALB REAL, ALLOCATABLE, DIMENSION(:) :: DHOM,FXSOR,XSCAT,GAMMA,V,FL,DFL REAL, ALLOCATABLE, DIMENSION(:,:) :: DIFHET,SFNU REAL, ALLOCATABLE, DIMENSION(:,:,:) :: FLUX DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: XCSOU *---- * FOR NUMERICAL FOURIER ANALYSIS *---- REAL OMEGA,XLEN,EKKK,EVALRHO,SPECR REAL, ALLOCATABLE, DIMENSION(:) :: ARRAYRHO1 REAL, ALLOCATABLE, DIMENSION(:) :: XXX *---- * DATA STATEMENTS *---- SAVE CAN,HTYPE DATA CAN /'00','01','02','03','04','05','06','07','08','09', > '10','11','12','13','14','15','16','17','18','19'/ DATA (HTYPE(JJ),JJ=0,5)/'S ','P ',2*'K ','B ','L '/ *---- * SCRATCH STORAGE ALLOCATION * DHOM homogeneous leakage coefficients. * DIFHET heterogeneous leakage coefficients. * FLUX iteration flux: * FLUX(:,:,1) <=old outer; * FLUX(:,:,2) <=present outer; * FLUX(:,:,3) <=new outer; * FLUX(:,:,4) <=source outer; * FLUX(:,:,5) <=old inner; * FLUX(:,:,6) <=present inner; * FLUX(:,:,7) <=new inner; * FLUX(:,:,8) <=source inner. *---- ALLOCATE(IJJ(0:NMAT),NJJ(0:NMAT),IPOS(0:NMAT),NPSYS(NGRP)) ALLOCATE(FLUX(NUNKNO,NGRP,8),XSCAT(0:NMAT*NGRP),GAMMA(NGRP), 1 DHOM(NGRP),DIFHET(NMERG,NGRP),XCSOU(NGRP)) * REBFLG=.TRUE. IPREB=IPMACR * AKEEP(:8)=0.0D0 ICHAR=0 ICTOT=0 *---- * RECOVER INDEX FOR THE CURRENTS IN FLUX, NUMERICAL SURFACES, * ALBEDO IF NEEDED BY THE REBALANCING. *---- ICREB=0 NNN=0 INSB=0 IBFP=0 NDIM=0 NFOU=0 LX=0 ITYPE=0 IF(CXDOOR.EQ.'MCCG') THEN CALL LCMGET(IPTRK,'STATE-VECTOR',JPAR) NANIS_TRK=JPAR(6) NDIM=JPAR(16) INSB=JPAR(22) CALL LCMLEN(IPTRK,'KEYCUR$MCCG',ICREB,ITYLCM) IF(ICREB.GT.0) THEN CALL LCMLEN(IPTRK,'NZON$MCCG',ILONG,ITYLCM) NNN=ILONG-ICREB ALLOCATE(KEYCUR(NSOUT),V(ILONG),MATALB(ILONG)) CALL LCMGET(IPTRK,'KEYCUR$MCCG',KEYCUR) CALL LCMGET(IPTRK,'V$MCCG',V) CALL LCMGET(IPTRK,'ALBEDO',ALBEDO) CALL LCMGET(IPTRK,'NZON$MCCG',MATALB) ENDIF ELSE IF(CXDOOR.EQ.'SN') THEN CALL LCMGET(IPTRK,'STATE-VECTOR',JPAR) ITYPE=JPAR(6) NDIM=JPAR(9) LX=JPAR(12) LY=JPAR(13) INSB=JPAR(27) IBFP=JPAR(31) NFOU=JPAR(34) ENDIF *---- * SELECT THE CALCULATION DOORS FOR WHICH A GROUP-BY-GROUP SCALAR * PROCEDURE WILL BE USED. A VECTORIAL APPROACH WILL BE USED WITH * OTHER DOORS. LSCAL is true for CXDOOR = 'TRAFIC'. *---- LSCAL=(INSB.EQ.0) * CALL KDRCPU(CPU0) IF(ILEAK.LT.6) THEN INORM=1 ELSE IF(ILEAK.EQ.6) THEN INORM=2 ! Ecco ELSE IF(ILEAK.GE.7) THEN INORM=3 ! Tibere ENDIF LEXAC=.FALSE. AKEEP(5)=1.0D0 AKEEP(6)=1.0D0 AKEEP(7)=1.0D0 DIFHET(:NMERG,:NGRP)=0.0 GAMMA(:NGRP)=1.0 *---- * EXTERNAL FLUX(:,:,2) INITIALISATION AND FIXED-EXTERNAL SOURCE IN * FLUX(:,:,4) *---- IF(ITYPEC.GE.3) THEN CALL LCMGET(IPFLUX,'B2 B1HOM',B2(4)) IF(ILEAK.GE.7) CALL LCMGET(IPFLUX,'B2 HETE',B2) ELSE B2(:4)=0.0 ENDIF AKEFFO=0.0D0 IF(LFORW) THEN IF(ITYPEC.EQ.1) THEN J1=LCMGID(IPSOU,'DSOUR') JPSOU=LCMGIL(J1,IGPT) J1=LCMGID(IPFLUX,'DFLUX') JPFLUX=LCMGIL(J1,IGPT) ELSE IF(C_ASSOCIATED(IPSOU)) THEN J1=LCMGID(IPSOU,'DSOUR') JPSOU=LCMGIL(J1,1) ENDIF JPFLUX=LCMGID(IPFLUX,'FLUX') JPSOUR=LCMLID(IPFLUX,'SOUR',NGRP) ENDIF ELSE IF(ITYPEC.EQ.1) THEN J1=LCMGID(IPSOU,'ASOUR') JPSOU=LCMGIL(J1,IGPT) J1=LCMGID(IPFLUX,'ADFLUX') JPFLUX=LCMGIL(J1,IGPT) ELSE IF(C_ASSOCIATED(IPSOU)) THEN J1=LCMGID(IPSOU,'ASOUR') JPSOU=LCMGIL(J1,1) ENDIF JPFLUX=LCMGID(IPFLUX,'AFLUX') JPSOUR=LCMLID(IPFLUX,'SOUR',NGRP) ENDIF ENDIF ALLOCATE(FXSOR(0:NMAT)) JPMACR=LCMGID(IPMACR,'GROUP') DO 20 IG=1,NGRP FLUX(:NUNKNO,IG,2)=0.0 FLUX(:NUNKNO,IG,4)=0.0 CALL LCMLEL(JPFLUX,1,ILINIT,ITYLCM) IF(LFORW) THEN CALL LCMGDL(JPFLUX,IG,FLUX(1,IG,2)) ELSE CALL LCMGDL(JPFLUX,NGRP-IG+1,FLUX(1,IG,2)) ENDIF * IF(((ITYPEC.EQ.0).OR.(ITYPEC.EQ.-2)).AND. 1 (.NOT.C_ASSOCIATED(IPSOU))) THEN KPMACR=LCMGIL(JPMACR,IG) FXSOR(0)=0.0 CALL LCMGET(KPMACR,'FIXE',FXSOR(1)) CALL DOORS(CXDOOR,IPTRK,NMAT,0,NUNKNO,FXSOR, > FLUX(1,IG,4)) ELSE IF(((ITYPEC.EQ.0).OR.(ITYPEC.EQ.1).OR.(ITYPEC.EQ.-2)) 1 .AND.C_ASSOCIATED(IPSOU))THEN IF(LFORW) THEN CALL LCMGDL(JPSOU,IG,FLUX(1,IG,4)) ELSE CALL LCMGDL(JPSOU,NGRP-IG+1,FLUX(1,IG,4)) ENDIF ENDIF 20 CONTINUE DEALLOCATE(FXSOR) *------- * IF IMPORTED FLUX PRESENT FOR SN, REORDER FLUX. *------- IF((CXDOOR.EQ.'SN').AND.(INITFL.EQ.2)) THEN CALL LCMLEN(IPFLUX,'KEYFLX',ILINIT,ITYLCM) IF(ILINIT.NE.NREG) THEN WRITE(*,*) NREG, ILINIT CALL XABORT('FLU2DR: NUMBER OF REGIONS FROM SPN CALCULATION' 1 //' (OBTAINED FROM LENGTH OF KEYFLX) DOES NOT MATCH NUMBER ' 2 //'OF REGIONS IN SN CALCULATION. CHECK INPUT FILE FOR ' 3 //'POTENTIAL ERRORS.') ENDIF KEYSPN(:) = 0 CALL LCMGET(IPFLUX,'KEYFLX',KEYSPN) DO 25 IG=1,NGRP CALL SNEST(IPTRK,IPRT,NREG,NUNKNO,MATCOD,IG,KEYFLX,KEYSPN, 1 FLUX(:,IG,2)) 25 CONTINUE ENDIF *---- * FOURIER ANALYSIS FLUX INITIALISATION *---- ALLOCATE(ARRAYRHO1(NFOU**NDIM)) ARRAYRHO1(:)=0.0 IFACOUNT=-1 26 IFACOUNT=IFACOUNT+1 * IF(ITYPEC.EQ.-2) THEN IF(NFOU.EQ.0) > CALL XABORT('FLU2DR: NEED TO SPECIFY FOURIER ANALYSIS ' > //'KEYWORD NFOU IN TRACKING, AS WELL AS NUMBER OF ' > //'FREQUENCIES TO INVESTIGATE.') IF(CXDOOR.NE.'SN') > CALL XABORT('FLU2DR: FOURIER ANALYSIS ONLY MEANT FOR SN') IF(NGRP.NE.1) > CALL XABORT('FLU2DR: FOURIER ANALYSIS NOT MEANT FOR MULTI-' > //'GROUP PROBLEMS. CONSIDER ADDING THAT FUNCTIONALITY. ') IGR=1 FLUX(:,:,:) = 0.0 SUMXSTRC=0.0 DO IR=1,NREG SUMXSTRC = SUMXSTRC + XSTRC(MATCOD(IR),1) ENDDO AVXSTRC = (SUMXSTRC/NREG) ALLOCATE(XXX(LX+1)) CALL LCMGET(IPTRK,'XXX',XXX) CALL LCMGET(IPTRK,'XLEN',XLEN) OMEGA = (2*PI)/(XLEN*AVXSTRC) EKKK =(REAL(IFACOUNT)/NFOU) IF(ITYPE.EQ.2)THEN PARTX=0.0 DO IX=1,LX PARTX = XXX(IX) + ((XXX(IX+1)-XXX(IX))/2) IND=KEYFLX(IX,1,1) IF(IND.GT.0) > FLUX(IND,IGR,2) = COS(EKKK*OMEGA*AVXSTRC*PARTX) ENDDO ELSE CALL XABORT('FLU2DR: FOURIER ANALYSIS FOR GEOMETRIES OTHER ' > //'THAN CARTESIAN 1D NOT AVAILABLE.') ENDIF DEALLOCATE(XXX) ENDIF *---- * COMPUTE FIRST K-EFFECTIVE *---- IF((ITYPEC.EQ.0).OR.(ITYPEC.EQ.5).OR.(ITYPEC.EQ.-2)) THEN AKEFFO=1.0D0 AKEFF=1.0D0 AFLNOR=1.0D0 ELSE IF(ITYPEC.EQ.1) THEN CALL LCMGET(IPFLUP,'STATE-VECTOR',JPAR) IF(JPAR(6).GE.2) THEN CALL LCMGET(IPFLUP,'K-EFFECTIVE',RKEFF) CALL LCMGET(IPFLUP,'K-INFINITY',CUREIN) ENDIF AKEFF=RKEFF IF(JPAR(6).GE.3) THEN B2(:4)=0.0 CALL LCMGET(IPFLUP,'B2 B1HOM',B2(4)) ENDIF IF((JPAR(6).GT.2).AND.(JPAR(7).GE.6)) THEN CALL LCMGET(IPFLUP,'B2 HETE',B2) ENDIF IF((JPAR(6).GT.2).AND.(JPAR(7).GE.5)) THEN CALL LCMGET(IPFLUP,'GAMMA',GAMMA) ENDIF AKEFFO=AKEFF AKEEP(2)=AKEFF AFLNOR=1.0D0/RKEFF ELSE OLDBIL=0.0D0 CALL FLUKEF(IPRT,IPMACR,NGRP,NREG,NUNKNO,NMAT,NIFIS,NANIS, 1 MATCOD(1),VOL,KEYFLX(1,1,1),XSTRC,XSDIA,XSNUF,XSCHI,NMERG, 2 IMERG,DIFHET,FLUX(1,1,2),B2,ILEAK,LEAKSW,OLDBIL,AKEFF,AFLNOR) AKEFFO=AKEFF AKEEP(2)=AKEFF ENDIF B2VALO=B2(4) * NCTOT=NCPTA+NCPTL IF(NCPTA.EQ.0) THEN NCPTM=NCTOT+1 ELSE NCPTM=NCPTL ENDIF MESSOU=' ' IF(IPRT.GT.0) WRITE(6,1090) 0,1.0,EPSOUT,AKEFFO,B2(4) *---- * CALCULATION OF THE INITIAL LEAKAGE COEFFICIENTS *---- IF(ITYPEC.GT.2) THEN DIFHET(:NMERG,:NGRP)=0.0 IF(OPTION.EQ.'LKRD') THEN CALL LCMGET(IPMACR,'STATE-VECTOR',JPAR) IF(JPAR(2).NE.NMAT) THEN CALL XABORT('FLU2DR: INVALID NMAT IN THE MACROLIB.') ELSE IF(JPAR(9).NE.1) THEN CALL XABORT('FLU2DR: INVALID LEAKAGE IN THE MACROLIB.') ENDIF ALLOCATE(FL(NMAT),DFL(NMAT)) JPMACR=LCMGID(IPMACR,'GROUP') DO IG=1,NGRP KPMACR=LCMGIL(JPMACR,IG) CALL LCMLEN(KPMACR,'DIFF',ILONG,ITYLCM) IF(ILONG.EQ.0) CALL XABORT('FLU2DR: UNABLE TO RECOVER T' > //'HE DIFF RECORD IN THE MACROLIB OBJECT.') IF(NMERG.EQ.NMAT) THEN CALL LCMGET(KPMACR,'DIFF',DIFHET(1,IG)) ELSE CALL LCMGET(KPMACR,'DIFF',DFL) CALL LCMGET(KPMACR,'FLUX-INTG',FL) DO INM=1,NMERG DDELN1=0.D0 DDELD1=0.D0 DO IBM=1,NMAT IF(IMERG(IBM).EQ.INM) THEN DDELN1=DDELN1+FL(IBM)/DFL(IBM) DDELD1=DDELD1+FL(IBM) ENDIF ENDDO IF(DDELN1.EQ.0.D0) CALL XABORT('FLU2DR: DDELN1=0.') DIFHET(INM,IG)=REAL(DDELD1/DDELN1) ENDDO ENDIF ENDDO DEALLOCATE(DFL,FL) GAMMA(:NGRP)=1.0 ELSE IF(OPTION.EQ.'RHS') THEN CALL LCMLEN(IPFLUX,'DIFFHET',ILONG,ITYLCM) IF(ILONG.EQ.NMERG*NGRP) THEN IF(LFORW) THEN CALL LCMGET(IPFLUX,'DIFFHET',DIFHET) ELSE * Permute the diffusion coefficients if the LKRD option * is set for an adjoint calculation CALL LCMGET(IPFLUX,'DIFFHET',DIFHET) DO INM=1,NMERG GAMMA(:NGRP)=DIFHET(INM,:NGRP) DO IG=1,NGRP DIFHET(INM,IG)=GAMMA(NGRP-IG+1) ENDDO ENDDO ENDIF ELSE CALL XABORT('FLU2DR: UNABLE TO RECOVER THE DIFFHET RECO' > //'RD IN THE FLUX OBJECT.(1)') ENDIF GAMMA(:NGRP)=1.0 ELSE IF((LEAKSW).OR.(ILEAK.EQ.5)) THEN * Todorova heterogeneous leakage model. CALL FLULPN(IPMACR,NUNKNO,OPTION,'DIFF',NGRP,NREG,NMAT, 1 VOL,MATCOD,NMERG,IMERG,KEYFLX(1,1,1),FLUX(1,1,2),B2(4), 2 IPRT,DIFHET,DHOM) GAMMA(:NGRP)=1.0 ELSE * FUNDAMENTAL MODE CONDITION. IF(NMERG.NE.1) CALL XABORT('FLU2DR: ONE LEAKAGE ZONE EXPEC' 1 //'TED.(1)') CALL B1HOM(IPMACR,LEAKSW,NUNKNO,OPTION,'DIFF',NGRP,NREG, 1 NMAT,NIFIS,VOL,MATCOD,KEYFLX(1,1,1),FLUX(1,1,2),REFKEF, 2 IPRT,DIFHET(1,1),GAMMA,AKEFF,INORM,B2) ENDIF CALL LCMPUT(IPFLUX,'B2 B1HOM',1,2,B2(4)) CALL LCMPUT(IPFLUX,'DIFFHET',NMERG*NGRP,2,DIFHET) ENDIF * **** OUTER LOOP ****************************************** IGDEB=1 CFLI=.FALSE. CEXE=.FALSE. DO 400 IT=1,MAXOUT CALL KDRCPU(CPU1) MESSIN=' ' *---- * FISSION SOURCE CALCULATION IN FLUX(:,:,4) *---- IF(((ITYPEC.EQ.0).OR.(ITYPEC.EQ.-2)) 1 .AND.(.NOT.C_ASSOCIATED(IPSOU))) THEN ALLOCATE(FXSOR(0:NMAT)) JPMACR=LCMGID(IPMACR,'GROUP') FLUX(:NUNKNO,:NGRP,4)=0.0 DO IG=1,NGRP KPMACR=LCMGIL(JPMACR,IG) FXSOR(0)=0.0 CALL LCMGET(KPMACR,'FIXE',FXSOR(1)) CALL DOORS(CXDOOR,IPTRK,NMAT,0,NUNKNO,FXSOR, > FLUX(1,IG,4)) ENDDO DEALLOCATE(FXSOR) ELSE IF((ITYPEC.EQ.0).OR.(ITYPEC.EQ.-2)) THEN DO IG=1,NGRP IF(LFORW) THEN CALL LCMGDL(JPSOU,IG,FLUX(1,IG,4)) ELSE CALL LCMGDL(JPSOU,NGRP-IG+1,FLUX(1,IG,4)) ENDIF ENDDO ELSE IF(ITYPEC.EQ.1) THEN DO IG=1,NGRP IF(LFORW) THEN CALL LCMGDL(JPSOU,IG,FLUX(1,IG,4)) ELSE CALL LCMGDL(JPSOU,NGRP-IG+1,FLUX(1,IG,4)) ENDIF DO IUN=1,NUNKNO FLUX(IUN,IG,4)=-FLUX(IUN,IG,4) ENDDO ENDDO ELSE FLUX(:NUNKNO,:NGRP,4)=0.0 ENDIF IF(NIFIS.GT.0) THEN IF(CXDOOR.EQ.'BIVAC') THEN CALL BIVFIS(IPTRK,NREG,NMAT,NIFIS,NUNKNO,NGRP,MATCOD,VOL, > XSCHI,XSNUF,FLUX(1,1,2),FLUX(1,1,4)) ELSE IF(CXDOOR.EQ.'TRIVAC') THEN CALL TRIFIS(IPTRK,NREG,NMAT,NIFIS,NUNKNO,NGRP,MATCOD,VOL, > XSCHI,XSNUF,FLUX(1,1,2),FLUX(1,1,4)) ELSE ALLOCATE(FXSOR(NUNKNO)) DO IS=1,NIFIS FXSOR(:NUNKNO)=0.0 DO IG=1,NGRP CALL DOORS(CXDOOR,IPTRK,NMAT,0,NUNKNO,XSNUF(0,IS,IG), > FXSOR,FLUX(1,IG,2)) ENDDO DO IR=1,NREG IBM=MATCOD(IR) IF(IBM.EQ.0) CYCLE DO IE=1,NLIN IND=KEYFLX(IR,IE,1) IF(IND.EQ.0) CYCLE DO IG=1,NGRP FLUX(IND,IG,4)=FLUX(IND,IG,4)+XSCHI(IBM,IS,IG)* > FXSOR(IND) ENDDO ENDDO ENDDO ENDDO ! IS DEALLOCATE(FXSOR) ENDIF FLUX(:NUNKNO,:NGRP,4)=FLUX(:NUNKNO,:NGRP,4)*REAL(AFLNOR) ENDIF *---- * VOLUME-INTEGRATED SOURCE CALCULATION FOR USE IN NEUTRON BALANCE *---- DO IG=1,NGRP XCSOU(IG)=0.0D0 DO IR=1,NREG IND=KEYFLX(IR,1,1) IF((CXDOOR.EQ.'BIVAC').OR.(CXDOOR.EQ.'TRIVAC')) THEN ! volumes are already included in the sources IF(IND.GT.0) XCSOU(IG)=XCSOU(IG)+FLUX(IND,IG,4) ELSE IF(IND.GT.0) XCSOU(IG)=XCSOU(IG)+FLUX(IND,IG,4)*VOL(IR) ENDIF ENDDO ENDDO ISBS=-1 IF(C_ASSOCIATED(IPSOU)) CALL LCMLEN(IPSOU,'NBS',ISBS,ITYLCM) *---- * SET THE STARTING ENERGY GROUP *---- DO 40 IG=1,NGRP IGDEB=IG IF(XCSOU(IG).NE.0.0.OR.ISBS.NE.0) GO TO 50 40 CONTINUE *---- * DOWNLOAD FROM EXTERNAL FLUX(:,:,2) TO PRESENT INTERNAL FLUX(:,:,6) *---- 50 FLUX(:NUNKNO,:NGRP,6)=FLUX(:NUNKNO,:NGRP,2) * **** INNER LOOP ****************************************** DO 270 JT=1,MAXINR FLUX(:NUNKNO,IGDEB:NGRP,7)=FLUX(:NUNKNO,IGDEB:NGRP,6) FLUX(:NUNKNO,IGDEB:NGRP,8)=FLUX(:NUNKNO,IGDEB:NGRP,4) JPMACR=LCMGID(IPMACR,'GROUP') DO 140 IG=IGDEB,NGRP *---- * PROCESS SELF-SCATTERING REDUCTION IN INNER SOURCES. *---- IF((ITPIJ.EQ.2).OR.(ITPIJ.EQ.4)) THEN IF((CXDOOR.EQ.'BIVAC').OR.(CXDOOR.EQ.'TRIVAC')) THEN CALL XABORT('FLU2DR: SCATTERING REDUCTION IS MANDATORY.') ENDIF CALL DOORS(CXDOOR,IPTRK,NMAT,NANIS,NUNKNO,XSDIA(0,0,IG), > FLUX(1,IG,8),FLUX(1,IG,7)) ENDIF IF(ILEAK.EQ.6) THEN * ECCO ISOTROPIC STREAMING MODEL. CCLBD=0.0 IF((ITPIJ.EQ.1).OR.(ITPIJ.EQ.3).AND.(OPTION.EQ.'B1')) > CCLBD=1.0-GAMMA(IG) DO 75 IE=1,NLIN DO 70 IR=1,NREG IBM=MATCOD(IR) IND=NUNKNO/2+KEYFLX(IR,IE,1) IF(IND.EQ.NUNKNO/2) GO TO 70 IF(OPTION(2:2).EQ.'1') THEN * B1 OR P1 CASE. IF(ITPIJ.EQ.2) THEN FLUX(IND,IG,8)=FLUX(IND,IG,8)+XSDIA(IBM,1,IG)* > FLUX(IND,IG,7) ENDIF ELSE IF(ITPIJ.EQ.1) THEN * B0, P0, B0TR OR P0TR CASE. FLUX(IND,IG,8)=FLUX(IND,IG,8)-XSDIA(IBM,1,IG)* > FLUX(IND,IG,7)*GAMMA(IG) ENDIF FLUX(IND,IG,8)=FLUX(IND,IG,8)+CCLBD*XSDIA(IBM,1,IG)* > FLUX(IND,IG,7) 70 CONTINUE 75 CONTINUE ELSE IF(ILEAK.GE.7) THEN * TIBERE ANISOTROPIC STREAMING MODEL. CCLBD=0.0 IF((ITPIJ.EQ.3).AND.(OPTION.EQ.'B1')) CCLBD=1.0-GAMMA(IG) DO 86 IE=1,NLIN DO 85 IR=1,NREG IND0=KEYFLX(IR,IE,1) IF(IND0.EQ.0) GO TO 85 IBM=MATCOD(IR) INDD(1)=NUNKNO/4+IND0 INDD(2)=NUNKNO/2+IND0 INDD(3)=3*NUNKNO/4+IND0 DO 80 IDIR=1,3 IND=INDD(IDIR) IF(OPTION(2:2).EQ.'1') THEN * B1 OR P1 CASE. IF(ITPIJ.EQ.4) THEN FLUX(IND,IG,8)=FLUX(IND,IG,8)+XSDIA(IBM,1,IG)* > FLUX(IND,IG,7) ENDIF ELSE IF(ITPIJ.EQ.3) THEN * B0, P0, B0TR OR P0TR CASE. FLUX(IND,IG,8)=FLUX(IND,IG,8)-XSDIA(IBM,1,IG)* > FLUX(IND,IG,7)*GAMMA(IG) ENDIF FLUX(IND,IG,8)=FLUX(IND,IG,8)+CCLBD*XSDIA(IBM,1,IG)* > FLUX(IND,IG,7) 80 CONTINUE 85 CONTINUE 86 CONTINUE ENDIF *---- * COMPUTE INNER SOURCES ASSUMING SELF-SCATTERING REDUCTION. *---- IF(.NOT.LSCAL) THEN KPMACR=LCMGIL(JPMACR,IG) NUNK2=NUNKNO IF(ILEAK.EQ.6) NUNK2=NUNKNO/2 IF(ILEAK.GE.7) NUNK2=NUNKNO/4 IF((CXDOOR.EQ.'SN').AND.(IBFP.EQ.0)) THEN NUNK2=NUNKNO CALL SNSOUR(NUNKNO,IG,IPTRK,KPMACR,NANIS,NREG,NMAT,NUNK2, 1 NGRP,MATCOD,FLUX(1,1,7),FLUX(1,1,8)) ELSE IF(CXDOOR.EQ.'SN') THEN NUNK2=NUNKNO IPSTR=LCMGID(IPSYS,'STREAMING') JPSTR=LCMGID(IPSTR,'GROUP') KPSYS=LCMGIL(JPSTR,IG) CALL SNSBFP(IG,IPTRK,KPMACR,KPSYS,NANIS,NLF,NREG,NMAT, 1 NUNK2,NGRP,MATCOD,FLUX(1,1,7),FLUX(1,1,8)) ELSE HLEAK=' ' CALL FLUSOU(CXDOOR,HLEAK,NUNKNO,IG,IPTRK,KPMACR,NMAT,NANIS, 1 NUNK2,NGRP,FLUX(1,1,7),FLUX(1,1,8)) ENDIF IF((ILEAK.EQ.6).AND.(OPTION(2:2).EQ.'1')) THEN * ECCO ISOTROPIC STREAMING MODEL. HLEAK='ECCO ' CALL FLUSOU(CXDOOR,HLEAK,NUNKNO,IG,IPTRK,KPMACR,NMAT,NANIS, 1 NUNK2,NGRP,FLUX(1,1,7),FLUX(1,1,8)) ELSE IF(ILEAK.GE.7) THEN * TIBERE ANISOTROPIC STREAMING MODEL. HLEAK='TIBERE' CALL FLUSOU(CXDOOR,HLEAK,NUNKNO,IG,IPTRK,KPMACR,NMAT,NANIS, 1 NUNK2,NGRP,FLUX(1,1,7),FLUX(1,1,8)) ENDIF ENDIF 140 CONTINUE *---- * FLUX COMPUTATION *---- NPSYS(:NGRP)=0 DO 150 IG=IGDEB,NGRP NPSYS(IG)=IG 150 CONTINUE JPSTR=C_NULL_PTR IF(C_ASSOCIATED(IPSYS)) THEN JPSYS=LCMGID(IPSYS,'GROUP') IF(ILEAK.EQ.6.OR.((MOD(ILEAK,10).EQ.7).AND.(IPHASE.EQ.1))) THEN IPSTR=LCMGID(IPSYS,'STREAMING') JPSTR=LCMGID(IPSTR,'GROUP') ENDIF ENDIF IF((.NOT.LSCAL).AND.(ILEAK.EQ.0)) THEN IDIR=0 CALL DOORFV(CXDOOR,JPSYS,NPSYS,IPTRK,IFTRAK,IPRT,NGRP, 1 NMAT,IDIR,NREG,NUNKNO,IPHASE,LEXAC,MATCOD,VOL,KEYFLX,TITLE, 2 FLUX(1,1,8),FLUX(1,1,7),IPREB,IPSOU,REBFLG,FLUXC,EVALRHO) ELSE IF(.NOT.LSCAL) THEN CALL FLUDBV(CXDOOR,IPHASE,JPSYS,JPSTR,NPSYS,IPTRK,IFTRAK, 1 IPRT,NREG,NUNKNO,NFUNL,NGRP,NMAT,NANIS,LEXAC,MATCOD,VOL,KEYFLX, 2 TITLE,ILEAK,LEAKSW,XSTRC,XSDIA,B2,NMERG,IMERG,DIFHET,GAMMA, 3 FLUX(1,1,2),FLUX(1,1,8),FLUX(1,1,7),IPREB,IPSOU,REBFLG,FLUXC) ELSE * A GROUP-BY-GROUP SCALAR PROCEDURE IS BEEN USED. IF(.NOT.C_ASSOCIATED(IPSYS)) THEN CALL XABORT('FLU2DR: MISSING L_PIJ OBJECT.') ENDIF KPSTR=C_NULL_PTR DO 230 IG=IGDEB,NGRP IF(IPRT.GT.10) WRITE(6,'(/25H FLU2DR: PROCESSING GROUP,I5, > 6H WITH ,A,1H.)') IG,CXDOOR KPMACR=LCMGIL(JPMACR,IG) NUNK2=NUNKNO IF(ILEAK.EQ.6) NUNK2=NUNKNO/2 IF(ILEAK.GE.7) NUNK2=NUNKNO/4 IF(CXDOOR.EQ.'BIVAC') THEN CALL BIVSOU(NUNKNO,IG,IPTRK,KPMACR,NANIS,NREG,NMAT,NUNK2, 1 NGRP,MATCOD,VOL,FLUX(1,1,7),FLUX(1,1,8)) ELSE IF(CXDOOR.EQ.'TRIVAC') THEN CALL TRIVSO(NUNKNO,IG,IPTRK,KPMACR,NANIS,NREG,NMAT,NUNK2, 1 NGRP,MATCOD,VOL,FLUX(1,1,7),FLUX(1,1,8)) ELSE IF((CXDOOR.EQ.'SN').AND.(IBFP.EQ.0)) THEN CALL SNSOUR(NUNKNO,IG,IPTRK,KPMACR,NANIS,NREG,NMAT,NUNK2, 1 NGRP,MATCOD,FLUX(1,1,7),FLUX(1,1,8)) ELSE IF(CXDOOR.EQ.'SN') THEN JPSYS=LCMGID(IPSYS,'GROUP') KPSYS=LCMGIL(JPSYS,IG) CALL SNSBFP(IG,IPTRK,KPMACR,KPSYS,NANIS,NLF,NREG,NMAT, 1 NUNK2,NGRP,MATCOD,FLUX(1,1,7),FLUX(1,1,8)) ELSE HLEAK=' ' CALL FLUSOU(CXDOOR,HLEAK,NUNKNO,IG,IPTRK,KPMACR,NMAT,NANIS, 1 NUNK2,NGRP,FLUX(1,1,7),FLUX(1,1,8)) ENDIF IF((ILEAK.EQ.6).AND.(OPTION(2:2).EQ.'1')) THEN * ECCO ISOTROPIC STREAMING MODEL. HLEAK='ECCO ' CALL FLUSOU(CXDOOR,HLEAK,NUNKNO,IG,IPTRK,KPMACR,NMAT,NANIS, 1 NUNK2,NGRP,FLUX(1,1,7),FLUX(1,1,8)) ELSE IF(ILEAK.GE.7) THEN * TIBERE ANISOTROPIC STREAMING MODEL. HLEAK='TIBERE' CALL FLUSOU(CXDOOR,HLEAK,NUNKNO,IG,IPTRK,KPMACR,NMAT,NANIS, 1 NUNK2,NGRP,FLUX(1,1,7),FLUX(1,1,8)) ENDIF * NPSYS(:NGRP)=0 NPSYS(IG)=IG IF(ILEAK.EQ.0) THEN IDIR=0 CALL DOORFV(CXDOOR,JPSYS,NPSYS,IPTRK,IFTRAK,IPRT,NGRP, 1 NMAT,IDIR,NREG,NUNKNO,IPHASE,LEXAC,MATCOD,VOL,KEYFLX,TITLE, 2 FLUX(1,1,8),FLUX(1,1,7),IPREB,IPSOU,REBFLG,FLUXC,EVALRHO) ELSE CALL FLUDBV(CXDOOR,IPHASE,JPSYS,JPSTR,NPSYS,IPTRK,IFTRAK, 1 IPRT,NREG,NUNKNO,NFUNL,NGRP,NMAT,NANIS,LEXAC,MATCOD,VOL, 2 KEYFLX,TITLE,ILEAK,LEAKSW,XSTRC,XSDIA,B2,NMERG,IMERG,DIFHET, 3 GAMMA,FLUX(1,1,2),FLUX(1,1,8),FLUX(1,1,7),IPREB,IPSOU,REBFLG, 4 FLUXC) ENDIF 230 CONTINUE ENDIF IF(LREBAL.AND.(ITYPEC.NE.5)) THEN CALL FLUBAL(IPMACR,NGRP,ILEAK,NMAT,NREG,ICREB,NUNKNO,NANIS, 1 MATCOD,VOL,KEYFLX(1,1,1),XSTRC,XSDIA,XCSOU,IGDEB,B2,NMERG, 2 IMERG,DIFHET,KEYCUR,MATALB(NNN+1),ALBEDO,V(NNN+1),FLUX(1,1,7)) ENDIF *---- * ACCELERATING INNER ITERATIONS CYCLICALLY DEPENDING ON PARAM. *---- IF(MOD(JT-1,NCTOT).GE.NCPTM) THEN CALL FLU2AC(NGRP,NUNKNO,IGDEB,FLUX(1,1,5),AKEEP(5),ZMU) ELSE ZMU=1.0 ENDIF *---- * CALCULATING ERROR AND PREC BETWEEN PRESENT AND NEW FLUX FOR * EACH GROUP. RETAIN LARGEST ERROR BETWEEN ANY GROUP. *---- EINN=0.0 ICHAR=ICHAR+1 ICTOT=ICTOT+NGRP-IGDEB+1 IGDEBO=IGDEB DO 260 IG=IGDEBO,NGRP GINN=0.0 FINN=0.0 DO 240 IR=1,NREG IND=KEYFLX(IR,1,1) IF(IND.EQ.0) GO TO 240 GINN=MAX(GINN,ABS(FLUX(IND,IG,6)-FLUX(IND,IG,7))) FINN=MAX(FINN,ABS(FLUX(IND,IG,7))) 240 CONTINUE FLUX(:NUNKNO,IG,5)=FLUX(:NUNKNO,IG,6) FLUX(:NUNKNO,IG,6)=FLUX(:NUNKNO,IG,7) GINN=GINN/FINN IF((GINN.LT.EPSINR).AND.(IGDEB.EQ.IG)) THEN IGDEB=IGDEB+1 ELSEIF((IGDEB.EQ.IG).AND.(IG.LT.NGRP)) THEN ERRDEB1=GINN ENDIF EINN=MAX(EINN,GINN) 260 CONTINUE * ITERF=JT IF(IPRT.GT.0) WRITE(6,1080) JT,EINN,EPSINR,IGDEB,ZMU IF((IPRT.GT.0).AND.(IGDEB.GT.1).AND.(IGDEB.LE.NGRP)) THEN WRITE(6,1082) ERRDEB1 ENDIF IF(EINN.LT.EPSINR) THEN * thermal convergence is reached CFLI=CEXE GOTO 280 ENDIF * near convergence (eps < 10.0 criterion) a new outer iteration * is started IF((IGDEB.GT.1).AND.(EINN.LT.10.*EPSINR)) GOTO 281 270 CONTINUE MESSIN='*NOT*' **** END OF INNER LOOP ****************************************** * 281 MESSIN='*NEARLY*' 280 IF(LREBAL) THEN IF(LEAKSW) THEN IF(ICREB.EQ.0) THEN WRITE(6,*) ' *** INCOMPATIBILITY ON LEAKAGE SWITCH ***' CALL XABORT('FLU2DR: ERROR ON LEAKAGE SWITCH') ELSE IF(IPRT.GT.0) & WRITE(6,*) 'FLU2DR: LEAKAGE & ICREB -> REBALANCING ON' ENDIF ELSE IF(IPRT.GT.0) & WRITE(6,*) 'FLU2DR: NO LEAKAGE-> REBALANCING ON' ENDIF ELSE IF(IPRT.GT.0) WRITE(6,*) 'FLU2DR: LEAKAGE-> REBALANCING OFF' ENDIF CALL KDRCPU(CPU2) IF(IPRT.GT.0) WRITE(6,1040) CPU2-CPU1,'INTERNAL',MESSIN,ITERF *---- * PROMOTE FROM NEW INTERNAL FLUX(,,,7) TO NEW EXTERNAL FLUX(,,,3) *---- FLUX(:NUNKNO,:NGRP,3)=FLUX(:NUNKNO,:NGRP,7) FLUX(:NUNKNO,:NGRP,4)=FLUX(:NUNKNO,:NGRP,8) *---- * HOTELLING DEFLATION IN GPT CASES. *---- IF(ITYPEC.EQ.1) THEN JPFLUP1=LCMGID(IPFLUP,'FLUX') JPFLUP2=LCMGID(IPFLUP,'AFLUX') DDELN1=0.0D0 DDELD1=0.0D0 DO 300 IG=1,NGRP IF(LFORW) THEN CALL LCMGDL(JPFLUP1,IG,FLUX(1,IG,5)) ! EVECT CALL LCMGDL(JPFLUP2,IG,FLUX(1,IG,6)) ! ADECT ELSE CALL LCMGDL(JPFLUP2,NGRP-IG+1,FLUX(1,IG,5)) ! ADECT CALL LCMGDL(JPFLUP1,NGRP-IG+1,FLUX(1,IG,6)) ! EVECT ENDIF 300 CONTINUE FLUX(:NUNKNO,:NGRP,7)=0.0 IF(CXDOOR.EQ.'BIVAC') THEN CALL BIVFIS(IPTRK,NREG,NMAT,NIFIS,NUNKNO,NGRP,MATCOD,VOL, > XSCHI,XSNUF,FLUX(1,1,6),FLUX(1,1,7)) ELSE IF(CXDOOR.EQ.'TRIVAC') THEN CALL TRIFIS(IPTRK,NREG,NMAT,NIFIS,NUNKNO,NGRP,MATCOD,VOL, > XSCHI,XSNUF,FLUX(1,1,6),FLUX(1,1,7)) ELSE ALLOCATE(FXSOR(NUNKNO)) DO IS=1,NIFIS FXSOR(:NUNKNO)=0.0 DO IG=1,NGRP CALL DOORS(CXDOOR,IPTRK,NMAT,0,NUNKNO,XSNUF(0,IS,IG), > FXSOR,FLUX(1,IG,6)) ENDDO DO IR=1,NREG IBM=MATCOD(IR) IF(IBM.EQ.0) CYCLE DO IE=1,NLIN IND=KEYFLX(IR,IE,1) IF(IND.EQ.0) CYCLE DO IG=1,NGRP FLUX(IND,IG,7)=FLUX(IND,IG,7)+XSCHI(IBM,IS,IG)* > FXSOR(IND) ENDDO ENDDO ENDDO ENDDO ! IS DEALLOCATE(FXSOR) ENDIF * DO 335 IG=1,NGRP DO 330 IND=1,NUNKNO DDELN1=DDELN1+FLUX(IND,IG,7)*FLUX(IND,IG,3) DDELD1=DDELD1+FLUX(IND,IG,7)*FLUX(IND,IG,5) 330 CONTINUE 335 CONTINUE DO 345 IG=1,NGRP DO 340 IND=1,NUNKNO FLUX(IND,IG,3)=FLUX(IND,IG,3)-REAL(DDELN1/DDELD1)*FLUX(IND,IG,5) 340 CONTINUE 345 CONTINUE ENDIF * IF(ITYPEC.EQ.2) THEN * NO B-N LEAKAGE CALCULATION REQUIRED IF(B2(4).NE.0.0) CALL XABORT('FLU2DR: NON ZERO BUCKLING.') CALL FLUKEF(IPRT,IPMACR,NGRP,NREG,NUNKNO,NMAT,NIFIS,NANIS, 1 MATCOD(1),VOL,KEYFLX(1,1,1),XSTRC,XSDIA,XSNUF,XSCHI,NMERG, 2 IMERG,DIFHET,FLUX(1,1,3),B2,ILEAK,LEAKSW,OLDBIL,AKEFF,AFLNOR) ELSE IF(ITYPEC.GT.2) THEN * PERFORM LEAKAGE CALCULATION. CALL LCMLEN(IPFLUX,'DIFFHET',ILONG,ITYLCM) IF(ILONG.EQ.0) THEN CALL XABORT('FLU2DR: UNABLE TO RECOVER THE DIFFHET RECORD ' 1 //'IN THE FLUX OBJECT.(2)') ENDIF CALL LCMGET(IPFLUX,'DIFFHET',DIFHET) GAMMA(:NGRP)=1.0 IF(ILEAK.EQ.5) THEN * Todorova heterogeneous leakage model. CALL FLULPN(IPMACR,NUNKNO,OPTION,HTYPE(ITYPEC),NGRP,NREG, 1 NMAT,VOL,MATCOD,NMERG,IMERG,KEYFLX(1,1,1),FLUX(1,1,3),B2(4), 2 IPRT,DIFHET,DHOM) IF(.NOT.LEAKSW) THEN CALL B1HOM(IPMACR,LEAKSW,NUNKNO,'LKRD',HTYPE(ITYPEC),NGRP, 1 NREG,NMAT,NIFIS,VOL,MATCOD,KEYFLX(1,1,1),FLUX(1,1,3), 2 REFKEF,IPRT,DHOM(1),GAMMA,AKEFF,INORM,B2) GO TO 350 ENDIF ENDIF IF(LEAKSW) THEN IF(HTYPE(ITYPEC).NE.'K') THEN CALL XABORT('FLU2DR: TYPE K EXPECTED WITH OPEN GEOMETRY.') ENDIF JPMACR=LCMGID(IPMACR,'GROUP') ALLOCATE(SFNU(NMAT,NIFIS)) PROD=0.0D0 DO IGR=1,NGRP KPMACR=LCMGIL(JPMACR,IGR) SFNU(:NMAT,:NIFIS)=0.0 IF(NIFIS.GT.0) CALL LCMGET(KPMACR,'NUSIGF',SFNU) DO IBM=1,NMAT FLXIN=0.0D0 DO I=1,NREG IND=KEYFLX(I,1,1) IF((MATCOD(I).EQ.IBM).AND.(IND.GT.0)) THEN FLXIN=FLXIN+FLUX(IND,IGR,3)*VOL(I) ENDIF ENDDO DO NF=1,NIFIS PROD=PROD+SFNU(IBM,NF)*FLXIN ENDDO ENDDO ENDDO DEALLOCATE(SFNU) AKEFF=AKEFF*PROD/OLDBIL OLDBIL=PROD IF(IPRT.GT.0) WRITE (6,1150) B2(4),AKEFF ELSE * FUNDAMENTAL MODE CONDITION. IF(NMERG.NE.1) CALL XABORT('FLU2DR: ONE LEAKAGE ZONE EXPEC' 1 //'TED.(2)') CALL B1HOM(IPMACR,LEAKSW,NUNKNO,OPTION,HTYPE(ITYPEC),NGRP, 1 NREG,NMAT,NIFIS,VOL,MATCOD,KEYFLX(1,1,1),FLUX(1,1,3), 2 REFKEF,IPRT,DIFHET(1,1),GAMMA,AKEFF,INORM,B2) IF(ILEAK.GE.7) THEN * COMPUTE THE DIRECTIONNAL BUCKLING COMPONENTS FOR TIBERE. IHETL=ILEAK/10-1 IF(IHETL.GT.0) THEN CALL FLUBLN(IPMACR,IPRT,NGRP,NMAT,NREG,NUNKNO,NIFIS, 1 MATCOD,VOL,KEYFLX(1,1,1),FLUX(1,1,3),IHETL,REFKEF,B2) ENDIF ENDIF ENDIF 350 CALL LCMPUT(IPFLUX,'B2 B1HOM',1,2,B2(4)) CALL LCMPUT(IPFLUX,'DIFFHET',NMERG*NGRP,2,DIFHET) ENDIF IF(ITYPEC.GE.3) THEN IF(B2(4).EQ.0.0) THEN BFLNOR=1.0D0 ELSE BFLNOR=1.0D0/ABS(B2(4)) ENDIF EEXT=REAL(ABS(B2(4)-B2VALO)*BFLNOR) B2VALO=B2(4) ENDIF IEXTF=IT IF((ITYPEC.GT.1).AND.(ITYPEC.LT.5)) THEN IF(AKEFF.NE.0.0) AFLNOR=1.0D0/AKEFF EEXT=REAL(ABS(AKEFF-AKEFFO)/AKEFF) ELSE EEXT=0.0 ENDIF AKEEP(3)=AKEFF *---- * ACCELERATING INNER ITERATIONS CYCLICALLY DEPENDING ON PARAM. *---- IF(MOD(IT-1,NCTOT).GE.NCPTM) THEN CALL FLU2AC(NGRP,NUNKNO,1,FLUX(1,1,1),AKEEP(1),ZMU) ELSE ZMU=1.0 ENDIF * EINN=0.0 IF(IPRT.GT.0) WRITE(6,1090) IT,EEXT,EPSOUT,AKEFF,B2(4) IF(EEXT.LT.EPSOUT) THEN * COMPARE FLUX FOR OUTER ITERATIONS DO 370 IG=1,NGRP GINN=0.0 FINN=0.0 DO 360 IR=1,NREG IND=KEYFLX(IR,1,1) IF(IND.EQ.0) GO TO 360 GINN=MAX(GINN,ABS(FLUX(IND,IG,2)-FLUX(IND,IG,3))) FINN=MAX(FINN,ABS(FLUX(IND,IG,3))) 360 CONTINUE FLUX(:NUNKNO,IG,1)=FLUX(:NUNKNO,IG,2) FLUX(:NUNKNO,IG,2)=FLUX(:NUNKNO,IG,3) GINN=GINN/FINN EINN=MAX(EINN,GINN) 370 CONTINUE IF(IPRT.GT.0) WRITE(6,1100) IT,EINN,EPSUNK,AFLNOR,ZMU CEXE=.TRUE. ELSE FLUX(:NUNKNO,:NGRP,1)=FLUX(:NUNKNO,:NGRP,2) FLUX(:NUNKNO,:NGRP,2)=FLUX(:NUNKNO,:NGRP,3) IF(IPRT.GT.0) WRITE(6,1110) IT,AFLNOR,ZMU ENDIF IF((ITYPEC.GE.2).AND.(AKEFF.NE.0.0)) THEN AFLNOR=(AKEFF/AKEEP(3))*AFLNOR ENDIF AKEEP(1)=AKEEP(2) AKEEP(2)=AKEEP(3) *---- * UPDATE KEFF *---- AKEFFO=AKEFF IF((EEXT.LT.EPSOUT).AND.(EINN.LT.EPSUNK).AND.(IT.GE.2)) GO TO 410 400 CONTINUE WRITE(6,*) '*** FLU2DR: CONVERGENCE NOT REACHED ***' WRITE(6,*) '*** FLU2DR: CONVERGENCE NOT REACHED ***' WRITE(6,*) '*** FLU2DR: CONVERGENCE NOT REACHED ***' MESSOU='*NOT*' * **** CONVERGENCE REACHED ****************************************** 410 RKEFF=REAL(AKEFF) IF(IPRT.GE.3) THEN WRITE(6,1010) (IR,IR=1,NREG) ALLOCATE(FL(NREG)) DO 425 IG=1,NGRP WRITE(6,1070) IG FL(:NREG)=0.0 DO 420 IR=1,NREG IND=KEYFLX(IR,1,1) IF(IND.GT.0) FL(IR)=FLUX(IND,IG,3) 420 CONTINUE WRITE(6,1020) (FL(IR),IR=1,NREG) 425 CONTINUE DEALLOCATE(FL) ENDIF IF(IPRT.GE.4) THEN ALLOCATE(FL(NREG)) DO 445 IG=1,NGRP WRITE(6,1070) IG DO 440 IA=2,NFUNL FL(:NREG)=0.0 DO 430 IR=1,NREG IND=KEYFLX(IR,1,IA) IF(IND.GT.0) FL(IR)=FLUX(IND,IG,3) 430 CONTINUE WRITE(6,1030) IA,(FL(IR),IR=1,NREG) 440 CONTINUE 445 CONTINUE DEALLOCATE(FL) ENDIF *---- * COMPUTE K-INF *---- IF(ITYPEC.GE.2) THEN FISOUR=0.0D0 OLDBIL=0.0D0 DO 490 IG=1,NGRP DO 460 IR=1,NREG IND=KEYFLX(IR,1,1) IF(IND.EQ.0) GO TO 460 DO 450 IS=1,NIFIS FISOUR=FISOUR+XSNUF(MATCOD(IR),IS,IG)*FLUX(IND,IG,3)*VOL(IR) 450 CONTINUE OLDBIL=OLDBIL+XSTRC(MATCOD(IR),IG)*FLUX(IND,IG,3)*VOL(IR) 460 CONTINUE KPMACR=LCMGIL(JPMACR,IG) CALL LCMGET(KPMACR,'NJJS00',NJJ(1)) CALL LCMGET(KPMACR,'IJJS00',IJJ(1)) CALL LCMGET(KPMACR,'IPOS00',IPOS(1)) CALL LCMGET(KPMACR,'SCAT00',XSCAT(1)) DO 480 IR=1,NREG IBM=MATCOD(IR) IF(IBM.GT.0) THEN IND=KEYFLX(IR,1,1) JG=IJJ(IBM) DO 470 JND=1,NJJ(IBM) IF(JG.EQ.IG) THEN OLDBIL=OLDBIL-XSDIA(IBM,0,IG)*FLUX(IND,JG,3)*VOL(IR) ELSE OLDBIL=OLDBIL-XSCAT(IPOS(IBM)+JND-1)*FLUX(IND,JG,3)*VOL(IR) ENDIF JG=JG-1 470 CONTINUE ENDIF 480 CONTINUE 490 CONTINUE CUREIN=0.0 IF(FISOUR.NE.0.0) CUREIN=REAL(FISOUR/OLDBIL) * * FLUX NORMALIZATION TO KEFF. IF(ITYPEC.LT.5) THEN FLUX(:NUNKNO,:NGRP,3)=FLUX(:NUNKNO,:NGRP,3)*REAL(AKEFF/FISOUR) FLUX(:NUNKNO,:NGRP,4)=FLUX(:NUNKNO,:NGRP,4)*REAL(AKEFF/FISOUR) ENDIF ENDIF *---- * PRINT TIME TAKEN *---- CALL KDRCPU(CPU1) IF(IPRT.GE.1) WRITE(6,1040) CPU1-CPU0,'EXTERNAL',MESSOU,IEXTF *---- * FOURIER ANALYSIS: STORE E-VALUE AND FIND SPECTRAL RADIUS *---- IF(ITYPEC.EQ.-2) THEN WRITE(6,1130) ARRAYRHO1(IFACOUNT+1)=EVALRHO IF(IFACOUNT.LT.(NFOU-1)) GO TO 26 SPECR=MAXVAL(ARRAYRHO1) WRITE(6,1140) SPECR CALL LCMPUT(IPFLUX,'SPEC-RADIUS',1,2,SPECR) ENDIF DEALLOCATE(ARRAYRHO1) *---- * PRINT TRACKING INFORMATION *---- IF(IPRT.GE.1) THEN IF((ITYPEC.EQ.0).OR.(ITYPEC.EQ.-2)) THEN WRITE(6,1050) ICHAR,EEXT ELSE WRITE(6,1060) ICHAR,CUREIN,AKEFF,B2(4),EEXT ENDIF WRITE(6,1120) ICTOT ENDIF *---- * RELEASE ARRAYS *---- IF(CXDOOR.EQ.'MCCG') THEN IF(ICREB.GT.0) DEALLOCATE(MATALB,V,KEYCUR) ENDIF *---- * SAVE THE SOLUTION *---- DO 510 IG=1,NGRP IF(LFORW) THEN CALL LCMPDL(JPFLUX,IG,NUNKNO,2,FLUX(1,IG,3)) CALL LCMPDL(JPSOUR,IG,NUNKNO,2,FLUX(1,IG,4)) ELSE CALL LCMPDL(JPFLUX,NGRP-IG+1,NUNKNO,2,FLUX(1,IG,3)) CALL LCMPDL(JPSOUR,NGRP-IG+1,NUNKNO,2,FLUX(1,IG,4)) ENDIF 510 CONTINUE IF(C_ASSOCIATED(IPSOU)) THEN CALL LCMLEN(IPSOU,'NORM-FS',ILEN,ITYLCM) IF(ILEN.GT.0) THEN CALL LCMGET(IPSOU,'NORM-FS',ZNORM) CALL LCMPUT(IPFLUX,'NORM-FS',1,2,ZNORM) CALL LCMPUT(IPFLUX,'MATCOD',NREG,1,MATCOD) ENDIF ENDIF IF(IBFP.NE.0) THEN CALL LCMGET(IPSYS,'ECUTOFF',ECUTOFF) CALL LCMPUT(IPFLUX,'ECUTOFF',1,2,ECUTOFF) CALL LCMPUT(IPFLUX,'FLUXC',NREG,2,FLUXC) ENDIF IF(ITYPEC.GE.2) THEN CALL LCMPUT(IPFLUX,'K-EFFECTIVE',1,2,RKEFF) CALL LCMPUT(IPFLUX,'K-INFINITY',1,2,CUREIN) ENDIF IF(ITYPEC.GE.3) THEN CALL LCMPUT(IPFLUX,'B2 B1HOM',1,2,B2(4)) ENDIF IF((ITYPEC.GT.2).AND.(ILEAK.GE.7)) THEN CALL LCMPUT(IPFLUX,'B2 HETE',3,2,B2) ENDIF IF((ITYPEC.GT.2).AND.(ILEAK.GE.6)) THEN CALL LCMPUT(IPFLUX,'GAMMA',NGRP,2,GAMMA) ENDIF *---- * SCRATCH STORAGE DEALLOCATION *---- DEALLOCATE(XCSOU,DIFHET,DHOM,GAMMA,XSCAT,FLUX) DEALLOCATE(NPSYS,IPOS,NJJ,IJJ) RETURN * 1010 FORMAT (/28H FLUXES AVERAGED :/ 1 (9X,2H/=,I4,:,6X,2H/=,I4,:,6X,2H/=,I4,:,6X,2H/=,I4,:,6X,2H/=, 2 I4,:,6X,2H/=,I4,:,6X,2H/=,I4,:,6X,2H/=,I4,:,6X,2H/=,I4,:,6X, 3 2H/=,I4)) 1020 FORMAT (7H FLUX ,2H: ,1P,10E12.5/(9X,10E12.5)) 1030 FORMAT (5H CUR ,I2,2H: ,1P,10E12.5/(9X,10E12.5)) 1040 FORMAT (18H FLU2DR: CPU TIME=, F10.0,1X,A8,13H CONVERGENCE , 1 A8,14H REACHED AFTER ,I6,12H ITERATIONS. ) 1050 FORMAT (/20H ++ TRACKING CALLED=,I4,6H TIMES , 1 11H PRECISION=,E9.2) 1060 FORMAT (/20H ++ TRACKING CALLED=,I4,6H TIMES , 1 12H FINAL KINF=,1P,E13.6, 2 12H FINAL KEFF=,E13.6,4H B2=,E12.5, 3 11H PRECISION=,E9.2) 1070 FORMAT (/14H ENERGY GROUP ,I6) 1080 FORMAT (10X,3HIN(,I3,6H) FLX:,5H PRC=,1P,E9.2,5H TAR=,E9.2, 1 7H IGDEB=, I13,6H ACCE=,0P,F12.5) 1082 FORMAT (18X,28HFIRST UNCONVERGED GROUP PRC=,1P,E9.2) 1090 FORMAT (5H OUT(,I3,6H) EIG:,5H PRC=,1P,E9.2,5H TAR=,E9.2, 1 6H KEFF=,E13.6,6H BUCK=,E12.5) 1100 FORMAT (5H OUT(,I3,6H) FLX:,5H PRC=,1P,E9.2,5H TAR=,E9.2, 1 6H FNOR=,E13.6,6H ACCE=,0P,F12.5) 1110 FORMAT (5H OUT(,I3,6H) FLX:,28X,6H FNOR=,1P,E13.6,6H ACCE=, 1 0P,F12.5) 1120 FORMAT (38H ++ TOTAL NUMBER OF FLUX CALCULATIONS=,I10) 1130 FORMAT (24H CONVERGENCE NOT SOUGHT.) 1140 FORMAT (49H FLU2DR: SPECTRAL RADIUS FOR FOURIER ANALYSIS IS , 1 E13.6) 1150 FORMAT(/18H FLU2DR: BUCKLING=,1P,E13.5,15H K-EFFECTIVE =,E13.5) END