*DECK SNFLUX SUBROUTINE SNFLUX(KPSYS,INCONV,NGIND,IPTRK,IMPX,NGEFF,NREG, 1 NBMIX,NUN,MAT,VOL,KEYFLX,FUNKNO,SUNKNO,ITER,NBS,KPSOU1,KPSOU2, 2 FLUXC) * *----------------------------------------------------------------------- * *Purpose: * Solve a single non-accelerated scattering iteration of the N-group * transport equation for fluxes using the discrete ordinates (SN) * method. * *Copyright: * Copyright (C) 2007 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): A. Hebert * *Parameters: input * KPSYS pointer to the assembly matrices. KPSYS is an array of * directories. * INCONV energy group convergence flag (set to .FALSE. if converged). * NGIND energy group indices assign to the NGEFF set. * IPTRK pointer to the tracking (L_TRACK signature). * IMPX print flag (equal to zero for no print). * NGEFF number of energy groups processed in parallel. * NREG total number of regions for which specific values of the * neutron flux and reactions rates are required. * NBMIX number of mixtures. * NUN total number of unknowns in vectors SUNKNO and FUNKNO. * MAT index-number of the mixture type assigned to each volume. * VOL volumes. * KEYFLX position of averaged flux elements in FUNKNO vector. * SUNKNO input source vector. * ITER number of previous calls to SNFLUX. * NBS * KPSOU1 * KPSOU2 * *Parameters: input/output * FUNKNO unknown vector. * FLUXC flux at the cutoff energy. * *----------------------------------------------------------------------- * USE GANLIB *---- * SUBROUTINE ARGUMENTS *---- TYPE(C_PTR) KPSYS(NGEFF),IPTRK,KPSOU1(NGEFF),KPSOU2(NGEFF) INTEGER NGEFF,NGIND(NGEFF),IMPX,NREG,NBMIX,NUN,MAT(NREG), 1 KEYFLX(NREG),ITER,NBS(NGEFF) LOGICAL INCONV(NGEFF) REAL VOL(NREG),FUNKNO(NUN,NGEFF),SUNKNO(NUN,NGEFF) REAL,OPTIONAL :: FLUXC(NREG) *---- * LOCAL VARIABLES *---- PARAMETER (IUNOUT=6,NSTATE=40) INTEGER IPAR(NSTATE),NCODE(6),BSINFO(2),EELEM,ESCHM,P,Q, > LOZSWP(3,6) REAL ZCODE(6),SIDE,LAMBDA0 LOGICAL LFIXUP,LDSA,LSHOOT,LADPT(4) *---- * ALLOCATABLE ARRAYS *--- INTEGER, ALLOCATABLE, DIMENSION(:) :: KEYSPN INTEGER, ALLOCATABLE, DIMENSION(:,:) :: COORDMAP INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: ISBSM REAL, ALLOCATABLE, DIMENSION(:,:) :: QEXT,OLD,SGAR,BS,EMOMTR, 1 SIGMATR,MTR REAL, ALLOCATABLE, DIMENSION(:,:,:) :: SGAS,ESTOPW * TYPE(C_PTR) U_PTR,W_PTR,PL_PTR,JOP_PTR,UPQ_PTR,WPQ_PTR,ALPHA_PTR, 1 PLZ_PTR,SURF_PTR,XXX_PTR,DU_PTR,DE_PTR,MRM_PTR,MRMX_PTR,MRMY_PTR, 2 MRMZ_PTR,DB_PTR,DA_PTR,DAL_PTR,DZ_PTR,DC_PTR,WX_PTR,WE_PTR, 3 CST_PTR,MN_PTR,DN_PTR,IL_PTR INTEGER, POINTER, DIMENSION(:) :: JOP,MRM,MRMX,MRMY,MRMZ,ISLG,IL REAL, POINTER, DIMENSION(:) :: U,W,PL,UPQ,WPQ,ALPHA,PLZ,SURF,DU, 1 DE,XXX,DB,DA,DAL,DZ,DC,DELTAE,WX,WE,CST,MN,DN *---- * SCRATCH STORAGE ALLOCATION *---- ALLOCATE(QEXT(NUN,NGEFF),KEYSPN(NREG),DELTAE(NGEFF),ISLG(NGEFF)) *---- * RECOVER SN SPECIFIC PARAMETERS. *---- CALL LCMGET(IPTRK,'STATE-VECTOR',IPAR) ITYPE=IPAR(6) NSCT=IPAR(7) IELEM=IPAR(8) ISCHM=IPAR(10) L4=IPAR(11) LX=IPAR(12) LY=IPAR(13) LZ=IPAR(14) NLF=IPAR(15) ISCAT=IPAR(16) LFIXUP=IPAR(18).EQ.1 LDSA=IPAR(19).EQ.1 NSDSA=IPAR(21) ISPLH=IPAR(26) NKBA=IPAR(28) IGAV=IPAR(29) LSHOOT=.TRUE. IF(IPAR(30).EQ.0) LSHOOT=.FALSE. IBFP=IPAR(31) EELEM=IPAR(35) ESCHM=IPAR(36) LADPT(1:4)=.FALSE. IF(ESCHM.EQ.3) LADPT(1)=.TRUE. IF(ISCHM.EQ.3) LADPT(2:4)=.TRUE. *---- * TEST FOR DSA and SAVE OLD FLUX. *---- LDSA=.FALSE. IF(MOD(ITER,(NSDSA+1))==0) LDSA=.TRUE. IF(LDSA)THEN ALLOCATE(OLD(NUN,NGEFF)) OLD(:NUN,:NGEFF)=FUNKNO(:NUN,:NGEFF) ENDIF *---- * RECOVER TOTAL AND WITHIN-GROUP SCATTERING MULTIGROUP CROSS SECTIONS. *---- ALLOCATE(SGAR(0:NBMIX,NGEFF),SGAS(0:NBMIX,ISCAT,NGEFF)) ALLOCATE(ESTOPW(0:NBMIX,2,NGEFF),EMOMTR(0:NBMIX,NGEFF)) NANI=1 DO 10 II=1,NGEFF IF(.NOT.INCONV(II)) GO TO 10 CALL LCMLEN(KPSYS(II),'DRAGON-TXSC',ILONG,ITYLCM) IF(ILONG.NE.NBMIX+1) CALL XABORT('SNFLUX: INVALID TXSC LENGTH.') CALL LCMLEN(KPSYS(II),'DRAGON-S0XSC',ILONG,ITYLCM) NANI=MAX(NANI,ILONG/(NBMIX+1)) IF(NANI.GT.ISCAT) CALL XABORT('SNFLUX: INVALID S0XSC LENGTH.') CALL LCMGET(KPSYS(II),'DRAGON-TXSC',SGAR(0,II)) CALL LCMGET(KPSYS(II),'DRAGON-S0XSC',SGAS(0,1,II)) *---- * TEST FOR FOKKER-PLANCK TREATMENT. *---- IF(IBFP.GT.0) THEN CALL LCMLEN(KPSYS(II),'DRAGON-ESTOP',IFP,ITYLCM) IF(IFP.NE.(NBMIX+1)*2) CALL XABORT('SNFLUX: INVALID ESTOPW LEN' 1 //'GTH.') CALL LCMGET(KPSYS(II),'DRAGON-ESTOP',ESTOPW(0,1,II)) CALL LCMGET(KPSYS(II),'DRAGON-DELTE',DELTAE(II)) CALL LCMGET(KPSYS(II),'DRAGON-ISLG',ISLG(II)) IF(ISLG(II).EQ.1) FLUXC(:)=0.0 CALL LCMLEN(KPSYS(II),'DRAGON-EMOMT',IFP,ITYLCM) IF(IFP.NE.NBMIX+1) CALL XABORT('SNFLUX: INVALID EMOMTR LENGTH.') CALL LCMGET(KPSYS(II),'DRAGON-EMOMT',EMOMTR(0,II)) ENDIF *---- * PRINT ZEROTH MOMENT OF SOURCES. *---- IF(IMPX.GT.3) THEN WRITE(IUNOUT,500) NGIND(II) WRITE(IUNOUT,'(1P,6(5X,E15.7))') (SUNKNO(KEYFLX(I),II),I=1,NREG) ENDIF 10 CONTINUE *---- * RECOVER INFORMATION ABOUT BOUNDARY SOURCES *---- ISBS=0 MAXL=0 IF(SUM(NBS).NE.0) THEN ISBS=1 IF(ITYPE.EQ.2) THEN MAXL=1 ALLOCATE(ISBSM(2,NLF,NGEFF)) ELSE IF(ITYPE.EQ.5) THEN MAXL=MAX(LX,LY) CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) ALLOCATE(ISBSM(4,NPQ,NGEFF)) ELSE IF(ITYPE.EQ.7) THEN MAXL=MAX(LX*LY,LX*LZ,LY*LZ) CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) ALLOCATE(ISBSM(6,NPQ,NGEFF)) ELSE CALL XABORT('SNFLUX: BOUNDARY SOURCE NOT IMPLEMENTED FOR THAT' 1 //' GEOMETRY') ALLOCATE(ISBSM(0,0,0),BS(0,0)) ENDIF ALLOCATE(BS(MAXL,SUM(NBS))) BS=0.0 ISBSM=0 JJ=1 DO 11 II=1,NGEFF IF(NBS(II).NE.0) THEN DO N=1,NBS(II) CALL LCMGDL(KPSOU1(II),N,BS(1,JJ)) CALL LCMGDL(KPSOU2(II),N,BSINFO) ISBSM(BSINFO(1),BSINFO(2),II)=JJ JJ=JJ+1 ENDDO ENDIF 11 CONTINUE ELSE ALLOCATE(ISBSM(0,0,0),BS(0,0)) ENDIF *---- * COMPUTE THE FLUX. *---- IF((ITYPE.EQ.2).AND.(IBFP.EQ.0)) THEN *------------ * 1D SLAB *------------ ! EXTRACTING PARAMETERS IF(IELEM*NSCT*NREG.GT.NUN) CALL XABORT('SNFLUX: QEXT OVERFLOW' 1 //'(1)') CALL LCMGPD(IPTRK,'U',U_PTR) CALL LCMGPD(IPTRK,'WX',WX_PTR) CALL LCMGPD(IPTRK,'CST',CST_PTR) CALL LCMGPD(IPTRK,'MN',MN_PTR) CALL LCMGPD(IPTRK,'DN',DN_PTR) CALL LCMGPD(IPTRK,'IL',IL_PTR) CALL C_F_POINTER(U_PTR,U,(/ NLF /)) CALL C_F_POINTER(WX_PTR,WX,(/ IELEM+1 /)) CALL C_F_POINTER(CST_PTR,CST,(/ IELEM /)) CALL C_F_POINTER(MN_PTR,MN,(/ NLF*NSCT /)) CALL C_F_POINTER(DN_PTR,DN,(/ NLF*NSCT /)) CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) CALL LCMGET(IPTRK,'NCODE',NCODE) CALL LCMGET(IPTRK,'ZCODE',ZCODE) ! LOOP FOR GROUPS DO 40 II=1,NGEFF IF(.NOT.INCONV(II)) GO TO 40 IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN/1D-slab' ! FIXED VOLUMIC SOURCES QEXT(:NUN,II)=SUNKNO(:NUN,II) ! IN-GROUP SCATTERING DO 30 IR=1,NREG IBM=MAT(IR) IF(IBM.EQ.0) GO TO 30 DO 20 P=1,NSCT IF(IL(P).GT.NANI-1) GO TO 20 DO 15 IEL=1,IELEM IND=(IR-1)*NSCT*IELEM+IELEM*(P-1)+IEL QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,IL(P)+1,II)*FUNKNO(IND,II) 15 CONTINUE 20 CONTINUE 30 CONTINUE ! ONE-SPEED FLUX CALCULATION CALL SNFBC1(NREG,NBMIX,IELEM,NLF,NSCT,U,MAT,VOL,SGAR(0,II), 1 NCODE,ZCODE,QEXT(1,II),LFIXUP,LSHOOT,ISBS,SUM(NBS), 2 ISBSM(1,1,II),BS,WX,CST,LADPT(2),NUN,FUNKNO(1,II),MN,DN) 40 CONTINUE ! END OF ENERGY LOOP ELSE IF(ITYPE.EQ.2) THEN *------------ * 1D SLAB BOLTZMANN-FOKKER-PLANCK *------------ ! EXTRACTING PARAMETERS NM=IELEM*EELEM IF((IELEM*NLF+NM*NSCT)*NREG.GT.NUN) THEN CALL XABORT('SNFLUX: QEXT OVERFLOW(1a)') ENDIF IOF=NM*NSCT*NREG+1 CALL LCMGPD(IPTRK,'U',U_PTR) CALL LCMGPD(IPTRK,'WX',WX_PTR) CALL LCMGPD(IPTRK,'WE',WE_PTR) CALL LCMGPD(IPTRK,'CST',CST_PTR) CALL LCMGPD(IPTRK,'MN',MN_PTR) CALL LCMGPD(IPTRK,'DN',DN_PTR) CALL LCMGPD(IPTRK,'IL',IL_PTR) CALL C_F_POINTER(U_PTR,U,(/ NLF /)) CALL C_F_POINTER(WX_PTR,WX,(/ IELEM+1 /)) CALL C_F_POINTER(WE_PTR,WE,(/ EELEM+1 /)) CALL C_F_POINTER(CST_PTR,CST,(/ MAX(IELEM,EELEM) /)) CALL C_F_POINTER(MN_PTR,MN,(/ NLF*NSCT /)) CALL C_F_POINTER(DN_PTR,DN,(/ NLF*NSCT /)) CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) CALL LCMGET(IPTRK,'NCODE',NCODE) CALL LCMGET(IPTRK,'ZCODE',ZCODE) ! ANGULAR FOKKER-PLANCK OPERATOR ALLOCATE(SIGMATR(NSCT,NSCT),MTR(NLF,NLF)) SIGMATR=0.0 DO P=1,NSCT SIGMATR(P,P)=-IL(P)*(IL(P)+1) ENDDO MTR=MATMUL(MATMUL(RESHAPE(MN,[NLF,NSCT]),SIGMATR), 1 RESHAPE(DN,[NSCT,NLF])) LAMBDA0=0.0 DO M=1,NLF IF(MTR(M,M).LT.-LAMBDA0) LAMBDA0=-MTR(M,M) ENDDO DO M=1,NLF MTR(M,M)=MTR(M,M)+LAMBDA0 ENDDO SIGMATR=MATMUL(MATMUL(RESHAPE(DN,[NSCT,NLF]),MTR), 1 RESHAPE(MN,[NLF,NSCT])) DEALLOCATE(MTR) DO II=1,NGEFF DO IBM=1,NBMIX SGAR(IBM,II)=SGAR(IBM,II)+LAMBDA0*EMOMTR(IBM,II) ENDDO ENDDO ! LOOP FOR GROUPS DO 80 II=1,NGEFF IF(.NOT.INCONV(II)) GO TO 80 IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN-BFP/1D-slab' ! FIXED VOLUMIC SOURCES QEXT(:NUN,II)=SUNKNO(:NUN,II) ! IN-GROUP SCATTERING DO 70 IR=1,NREG IBM=MAT(IR) IF(IBM.EQ.0) GO TO 70 DO 60 P=1,NSCT IF(IL(P).GT.NANI-1) GO TO 60 DO 50 IEL=1,NM IND=(IR-1)*NSCT*NM+NM*(P-1)+IEL QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,IL(P)+1,II)*FUNKNO(IND,II) DO Q=1,NSCT JND=(IR-1)*NSCT*NM+NM*(Q-1)+IEL QEXT(IND,II)=QEXT(IND,II)+EMOMTR(IBM,II)*SIGMATR(P,Q) 1 *FUNKNO(JND,II) ENDDO 50 CONTINUE 60 CONTINUE 70 CONTINUE ! ONE-SPEED FLUX CALCULATION CALL SNFE1D(NREG,NBMIX,IELEM,EELEM,NM,NLF,NSCT,U,MAT,VOL, 1 SGAR(0,II),ESTOPW(0,1,II),NCODE,ZCODE,DELTAE(II),QEXT(1,II), 2 LFIXUP,LSHOOT,FUNKNO(1,II),ISBS,SUM(NBS),ISBSM(1,1,II),BS,WX,WE, 3 CST,LADPT(1:2),IBFP,NUN,MN,DN) ! SOLUTION FLUX AT THE CUTOFF ENERGY IF(ISLG(II).EQ.1) THEN LXNI=EELEM*NLF IF(LSHOOT) LXNI=0 DO 71 IR=1,NREG IBM=MAT(IR) FLUXC(IR)=0.0 DO 72 IP=1,NLF IND=NM*NSCT*NREG+LXNI+(IR-1)*NLF*IELEM+(IP-1)*IELEM+1 JND=(IP-1)*NSCT+1 FLUXC(IR)=FLUXC(IR)+FUNKNO(IND,II)*DN(JND) 72 CONTINUE 71 CONTINUE ENDIF 80 CONTINUE ! END OF ENERGY LOOP DEALLOCATE(SIGMATR) ELSE IF(ITYPE.EQ.3) THEN *------------ * 1D TUBE/CYLINDRICAL *------------ ! EXTRACTING PARAMETERS IF(IELEM.NE.1) CALL XABORT('SNFLUX: DIAM 0 EXPECTED(1).') IF(NSCT*NREG.GT.NUN) CALL XABORT('SNFLUX: QEXT OVERFLOW(2)') M2=NLF/2 CALL LCMLEN(IPTRK,'UPQ',NPQ,ITYLCM) CALL LCMGPD(IPTRK,'JOP',JOP_PTR) CALL LCMGPD(IPTRK,'U',U_PTR) CALL LCMGPD(IPTRK,'UPQ',UPQ_PTR) CALL LCMGPD(IPTRK,'WPQ',WPQ_PTR) CALL LCMGPD(IPTRK,'ALPHA',ALPHA_PTR) CALL LCMGPD(IPTRK,'PLZ',PLZ_PTR) CALL LCMGPD(IPTRK,'PL',PL_PTR) CALL LCMGPD(IPTRK,'SURF',SURF_PTR) CALL LCMGET(IPTRK,'ZCODE',ZCODE) CALL C_F_POINTER(JOP_PTR,JOP,(/ M2 /)) CALL C_F_POINTER(U_PTR,U,(/ NPQ /)) CALL C_F_POINTER(UPQ_PTR,UPQ,(/ NPQ /)) CALL C_F_POINTER(WPQ_PTR,WPQ,(/ NPQ /)) CALL C_F_POINTER(ALPHA_PTR,ALPHA,(/ NPQ /)) CALL C_F_POINTER(PLZ_PTR,PLZ,(/ NSCT*M2 /)) CALL C_F_POINTER(PL_PTR,PL,(/ NSCT*NPQ /)) CALL C_F_POINTER(SURF_PTR,SURF,(/ NREG+1 /)) ! LOOP FOR GROUPS DO 120 II=1,NGEFF IF(.NOT.INCONV(II)) GO TO 120 IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN/1D-cyl' ! FIXED VOLUMIC SOURCES QEXT(:NUN,II)=SUNKNO(:NUN,II) ! IN-GROUP SCATTERING DO 110 IR=1,NREG IBM=MAT(IR) IF(IBM.EQ.0) GO TO 110 IOF=0 DO 100 P=0,NANI-1 DO 90 IM=0,P IF(MOD(P+IM,2).EQ.1) GO TO 90 IOF=IOF+1 IND=(IR-1)*NSCT+IOF QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,P+1,II)*FUNKNO(IND,II) 90 CONTINUE 100 CONTINUE 110 CONTINUE ! ONE-SPEED FLUX CALCULATION CURR=0.0 CALL SNFT1C(NREG,NBMIX,M2,NPQ,ISCAT,NSCT,JOP,U,UPQ,WPQ,ALPHA, 1 PLZ,PL,MAT,VOL,SURF,SGAR(0,II),IGAV,QEXT(1,II),LFIXUP,CURR, 2 FUNKNO(1,II)) IF(ZCODE(2).NE.0.0) THEN CA=CURR CB=1.0 CALL SNFT1C(NREG,NBMIX,M2,NPQ,ISCAT,NSCT,JOP,U,UPQ,WPQ,ALPHA, 1 PLZ,PL,MAT,VOL,SURF,SGAR(0,II),IGAV,QEXT(1,II),LFIXUP,CB, 2 FUNKNO(1,II)) CURR=ZCODE(2)*CA/(1.0+ZCODE(2)*(CA-CB)) CALL SNFT1C(NREG,NBMIX,M2,NPQ,ISCAT,NSCT,JOP,U,UPQ,WPQ,ALPHA, 1 PLZ,PL,MAT,VOL,SURF,SGAR(0,II),IGAV,QEXT(1,II),LFIXUP,CURR, 2 FUNKNO(1,II)) ENDIF 120 CONTINUE ! END OF ENERGY LOOP ELSE IF(ITYPE.EQ.4) THEN *------------ * 1D SPHERE *------------ ! EXTRACTING PARAMETERS IF(IELEM.NE.1) CALL XABORT('SNFLUX: DIAM 0 EXPECTED(2).') NSCT=ISCAT IF(NSCT*NREG.GT.NUN) CALL XABORT('SNFLUX: QEXT OVERFLOW(3)') CALL LCMGPD(IPTRK,'U',U_PTR) CALL LCMGPD(IPTRK,'W',W_PTR) CALL LCMGPD(IPTRK,'ALPHA',ALPHA_PTR) CALL LCMGPD(IPTRK,'PLZ',PLZ_PTR) CALL LCMGPD(IPTRK,'PL',PL_PTR) CALL LCMGPD(IPTRK,'SURF',SURF_PTR) CALL LCMGPD(IPTRK,'XXX',XXX_PTR) CALL LCMGET(IPTRK,'ZCODE',ZCODE) CALL C_F_POINTER(U_PTR,U,(/ NLF /)) CALL C_F_POINTER(W_PTR,W,(/ NLF /)) CALL C_F_POINTER(ALPHA_PTR,ALPHA,(/ NLF /)) CALL C_F_POINTER(PLZ_PTR,PLZ,(/ NSCT /)) CALL C_F_POINTER(PL_PTR,PL,(/ NSCT*NLF /)) CALL C_F_POINTER(SURF_PTR,SURF,(/ NREG+1 /)) CALL C_F_POINTER(XXX_PTR,XXX,(/ NREG+1 /)) ! LOOP FOR GROUPS DO 150 II=1,NGEFF IF(.NOT.INCONV(II)) GO TO 150 IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN/1D-sph' ! FIXED VOLUMIC SOURCES QEXT(:NUN,II)=SUNKNO(:NUN,II) ! IN-GROUP SCATTERING DO 140 IR=1,NREG IBM=MAT(IR) IF(IBM.EQ.0) GO TO 140 DO 130 P=0,NANI-1 IND=(IR-1)*NSCT+P+1 QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,P+1,II)*FUNKNO(IND,II) 130 CONTINUE 140 CONTINUE ! ONE-SPEED FLUX CALCULATION CURR=0.0 CALL SNFT1S(NREG,NBMIX,NLF,NSCT,U,W,ALPHA,PLZ,PL,MAT,VOL,SURF, 1 XXX,SGAR(0,II),IGAV,QEXT(1,II),LFIXUP,CURR,FUNKNO(1,II)) IF(ZCODE(2).NE.0.0) THEN CA=CURR CB=1.0 CALL SNFT1S(NREG,NBMIX,NLF,NSCT,U,W,ALPHA,PLZ,PL,MAT,VOL, 1 SURF,XXX,SGAR(0,II),IGAV,QEXT(1,II),LFIXUP,CB,FUNKNO(1,II)) CURR=ZCODE(2)*CA/(1.0+ZCODE(2)*(CA-CB)) CALL SNFT1S(NREG,NBMIX,NLF,NSCT,U,W,ALPHA,PLZ,PL,MAT,VOL, 1 SURF,XXX,SGAR(0,II),IGAV,QEXT(1,II),LFIXUP,CURR,FUNKNO(1,II)) ENDIF 150 CONTINUE ! END OF ENERGY LOOP ELSE IF((ITYPE.EQ.5).AND.(IBFP.EQ.0)) THEN *------------ * 2D CARTESIAN *------------ ! EXTRACTING PARAMETERS NM=IELEM**2 NMX=IELEM IF(NM*NSCT*NREG.GT.NUN) CALL XABORT('SNFLUX: QEXT OVE' 1 //'RFLOW(4a)') CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) CALL LCMGPD(IPTRK,'DU',DU_PTR) CALL LCMGPD(IPTRK,'DE',DE_PTR) CALL LCMGPD(IPTRK,'W',W_PTR) CALL LCMGPD(IPTRK,'MRM',MRM_PTR) CALL LCMGPD(IPTRK,'MRMY',MRMY_PTR) CALL LCMGPD(IPTRK,'DB',DB_PTR) CALL LCMGPD(IPTRK,'DA',DA_PTR) CALL LCMGPD(IPTRK,'WX',WX_PTR) CALL LCMGPD(IPTRK,'CST',CST_PTR) CALL LCMGPD(IPTRK,'IL',IL_PTR) CALL LCMGPD(IPTRK,'MN',MN_PTR) CALL LCMGPD(IPTRK,'DN',DN_PTR) CALL C_F_POINTER(DU_PTR,DU,(/ NPQ /)) CALL C_F_POINTER(DE_PTR,DE,(/ NPQ /)) CALL C_F_POINTER(W_PTR,W,(/ NPQ /)) CALL C_F_POINTER(MRM_PTR,MRM,(/ NPQ /)) CALL C_F_POINTER(MRMY_PTR,MRMY,(/ NPQ /)) CALL C_F_POINTER(DB_PTR,DB,(/ LX*NPQ /)) CALL C_F_POINTER(DA_PTR,DA,(/ LX*LY*NPQ /)) CALL C_F_POINTER(WX_PTR,WX,(/ IELEM+1 /)) CALL C_F_POINTER(CST_PTR,CST,(/ IELEM /)) CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) CALL C_F_POINTER(MN_PTR,MN,(/ NPQ*NSCT /)) CALL C_F_POINTER(DN_PTR,DN,(/ NPQ*NSCT /)) CALL LCMGET(IPTRK,'NCODE',NCODE) CALL LCMGET(IPTRK,'ZCODE',ZCODE) ! LOOP FOR GROUPS DO 200 II=1,NGEFF IF(.NOT.INCONV(II)) GO TO 200 IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN/2D-car' ! FIXED VOLUMIC SOURCES QEXT(:NUN,II)=SUNKNO(:NUN,II) ! IN-GROUP SCATTERING DO 190 IR=1,NREG IBM=MAT(IR) IF(IBM.EQ.0) GO TO 190 IOF=0 DO 180 P=1,NSCT IF(IL(P).GT.NANI-1) GO TO 180 DO 160 IEL=1,NM IND=(IR-1)*NSCT*NM+(P-1)*NM+IEL QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,IL(P)+1,II)*FUNKNO(IND,II) 160 CONTINUE 180 CONTINUE 190 CONTINUE 200 CONTINUE ! END OF ENERGY LOOP IF(NKBA.EQ.0)THEN CALL SNFBC2(NUN,NGEFF,IMPX,INCONV,NGIND,LX,LY,IELEM,NM,NMX, 1 NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR,NCODE,ZCODE,QEXT,LFIXUP, 2 DU,DE,W,MRM,MRMY,DB,DA,MN,DN,WX,WX,CST,LADPT(2:3), 3 ISBS,SUM(NBS),ISBSM,BS,MAXL,FUNKNO) ELSE IF(NKBA.GT.0)THEN IF(ISBS.EQ.1) CALL XABORT('SNFLUX: BOUNDARY SOURCES NOT YET ' 1 //'VERIFIED WITH KBA ALGORITHM. YET TO BE VALIDATED. COMMENT ' 2 //'CALL TO XABORT, RECOMPILE AND PROCEED WITH CAUTION.') IF(ISCHM.EQ.3) CALL XABORT('SNFLUX: AWD METHOD NOT YET ' 1 //'VERIFIED WITH KBA ALGORITHM. YET TO BE VALIDATED. COMMENT ' 2 //'CALL TO XABORT, RECOMPILE AND PROCEED WITH CAUTION.') CALL SNFKC2(NKBA,NUN,NGEFF,IMPX,INCONV,NGIND,LX,LY,IELEM,NM, 1 NMX,NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR,NCODE,ZCODE,QEXT,LFIXUP, 2 DU,DE,W,MRM,MRMY,DB,DA,MN,DN,WX,WX,CST,LADPT(2:3), 3 ISBS,SUM(NBS),ISBSM,BS,MAXL,FUNKNO) ENDIF ELSE IF(ITYPE.EQ.5) THEN *------------ * 2D CARTESIAN BOLTZMANN-FOKKER-PLANCK *------------ ! EXTRACTING PARAMETERS NM=IELEM**2*EELEM NME=IELEM**2 NMX=IELEM*EELEM CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) IF((NSCT*NM+NPQ*NME)*NREG.GT.NUN) THEN CALL XABORT('SNFLUX: QEXT OVERFLOW(4b)') ENDIF CALL LCMGPD(IPTRK,'DU',DU_PTR) CALL LCMGPD(IPTRK,'DE',DE_PTR) CALL LCMGPD(IPTRK,'W',W_PTR) CALL LCMGPD(IPTRK,'MRM',MRM_PTR) CALL LCMGPD(IPTRK,'MRMY',MRMY_PTR) CALL LCMGPD(IPTRK,'DB',DB_PTR) CALL LCMGPD(IPTRK,'DA',DA_PTR) CALL LCMGPD(IPTRK,'WX',WX_PTR) CALL LCMGPD(IPTRK,'WE',WE_PTR) CALL LCMGPD(IPTRK,'CST',CST_PTR) CALL LCMGPD(IPTRK,'IL',IL_PTR) CALL LCMGPD(IPTRK,'MN',MN_PTR) CALL LCMGPD(IPTRK,'DN',DN_PTR) CALL C_F_POINTER(DU_PTR,DU,(/ NPQ /)) CALL C_F_POINTER(DE_PTR,DE,(/ NPQ /)) CALL C_F_POINTER(W_PTR,W,(/ NPQ /)) CALL C_F_POINTER(MRM_PTR,MRM,(/ NPQ /)) CALL C_F_POINTER(MRMY_PTR,MRMY,(/ NPQ /)) CALL C_F_POINTER(DB_PTR,DB,(/ LX*NPQ /)) CALL C_F_POINTER(DA_PTR,DA,(/ LX*LY*NPQ /)) CALL C_F_POINTER(WX_PTR,WX,(/ IELEM+1 /)) CALL C_F_POINTER(WE_PTR,WE,(/ EELEM+1 /)) CALL C_F_POINTER(CST_PTR,CST,(/ MAX(IELEM,EELEM) /)) CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) CALL C_F_POINTER(MN_PTR,MN,(/ NPQ*NSCT /)) CALL C_F_POINTER(DN_PTR,DN,(/ NPQ*NSCT /)) CALL LCMGET(IPTRK,'NCODE',NCODE) CALL LCMGET(IPTRK,'ZCODE',ZCODE) ! ANGULAR FOKKER-PLANCK OPERATOR ALLOCATE(SIGMATR(NSCT,NSCT),MTR(NPQ,NPQ)) SIGMATR=0.0 DO P=1,NSCT SIGMATR(P,P)=-IL(P)*(IL(P)+1) ENDDO MTR=MATMUL(MATMUL(RESHAPE(MN,[NPQ,NSCT]),SIGMATR), 1 RESHAPE(DN,[NSCT,NPQ])) LAMBDA0=0.0 DO M=1,NPQ IF(MTR(M,M).LT.-LAMBDA0) LAMBDA0=-MTR(M,M) ENDDO DO M=1,NPQ MTR(M,M)=MTR(M,M)+LAMBDA0 ENDDO SIGMATR=MATMUL(MATMUL(RESHAPE(DN,[NSCT, NPQ]),MTR), 1 RESHAPE(MN,[NPQ,NSCT])) DEALLOCATE(MTR) DO II=1,NGEFF DO IBM=1,NBMIX SGAR(IBM,II)=SGAR(IBM,II)+LAMBDA0*EMOMTR(IBM,II) ENDDO ENDDO ! LOOP FOR GROUPS DO 250 II=1,NGEFF IF(.NOT.INCONV(II)) GO TO 250 IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN-BFP/2D-car' ! FIXED VOLUMIC SOURCES QEXT(:NUN,II)=SUNKNO(:NUN,II) ! IN-GROUP SCATTERING DO 240 IR=1,NREG IBM=MAT(IR) IF(IBM.EQ.0) GO TO 240 IOF=0 DO 230 P=1,NSCT IF(IL(P).GT.NANI-1) GO TO 230 DO 210 IEL=1,NM IND=(IR-1)*NSCT*NM+(P-1)*NM+IEL QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,IL(P)+1,II)*FUNKNO(IND,II) DO Q=1,NSCT JND=(IR-1)*NSCT*NM+(Q-1)*NM+IEL QEXT(IND,II)=QEXT(IND,II)+EMOMTR(IBM,II)*SIGMATR(P,Q) 1 *FUNKNO(JND,II) ENDDO 210 CONTINUE 230 CONTINUE 240 CONTINUE 250 CONTINUE ! END OF ENERGY LOOP DEALLOCATE(SIGMATR) ! FLUX CALCULATION CALL SNFE2D(NUN,NGEFF,IMPX,INCONV,NGIND,LX,LY,IELEM, 1 EELEM,NM,NME,NMX,NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR,ESTOPW, 2 NCODE,ZCODE,DELTAE,QEXT,LFIXUP,DU,DE,W,MRM,MRMY,DB,DA,FUNKNO, 3 ISLG,FLUXC,ISBS,SUM(NBS),ISBSM,BS,MAXL,WX,WX,WE,CST,LADPT(1:3), 4 IBFP,MN,DN) ELSE IF(ITYPE.EQ.6) THEN *------------ * TUBE 2D (R-Z) *------------ ! EXTRACTING PARAMETERS IF(IELEM.NE.1) CALL XABORT('SNFLUX: DIAM 0 EXPECTED(3).') IF(NSCT*NREG.GT.NUN) CALL XABORT('SNFLUX: QEXT OVERFLOW(5)') CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) CALL LCMGPD(IPTRK,'DU',DU_PTR) CALL LCMGPD(IPTRK,'DE',DE_PTR) CALL LCMGPD(IPTRK,'W',W_PTR) CALL LCMGPD(IPTRK,'MRM',MRM_PTR) CALL LCMGPD(IPTRK,'MRMY',MRMY_PTR) CALL LCMGPD(IPTRK,'DB',DB_PTR) CALL LCMGPD(IPTRK,'DA',DA_PTR) CALL LCMGPD(IPTRK,'DAL',DAL_PTR) CALL LCMGPD(IPTRK,'IL',IL_PTR) CALL LCMGPD(IPTRK,'MN',MN_PTR) CALL LCMGPD(IPTRK,'DN',DN_PTR) CALL C_F_POINTER(DU_PTR,DU,(/ NPQ /)) CALL C_F_POINTER(DE_PTR,DE,(/ NPQ /)) CALL C_F_POINTER(W_PTR,W,(/ NPQ /)) CALL C_F_POINTER(MRM_PTR,MRM,(/ NPQ /)) CALL C_F_POINTER(MRMY_PTR,MRMY,(/ NPQ /)) CALL C_F_POINTER(DB_PTR,DB,(/ LX*NPQ /)) CALL C_F_POINTER(DA_PTR,DA,(/ LX*LY*NPQ /)) CALL C_F_POINTER(DAL_PTR,DAL,(/ LX*LY*NPQ /)) CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) CALL C_F_POINTER(MN_PTR,MN,(/ NPQ*NSCT /)) CALL C_F_POINTER(DN_PTR,DN,(/ NPQ*NSCT /)) CALL LCMGET(IPTRK,'NCODE',NCODE) CALL LCMGET(IPTRK,'ZCODE',ZCODE) ! LOOP FOR GROUPS DO 290 II=1,NGEFF IF(.NOT.INCONV(II)) GO TO 290 IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN/2D-rz' ! FIXED VOLUMIC SOURCES QEXT(:NUN,II)=SUNKNO(:NUN,II) ! IN-GROUP SCATTERING DO 280 IR=1,NREG IBM=MAT(IR) IF(IBM.EQ.0) GO TO 280 IOF=0 DO 270 P=1,NSCT IF(IL(P).GT.MIN(ISCAT,NANI)-1) GO TO 270 IND=(IR-1)*NSCT+(P-1)*IELEM*IELEM+1 QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,IL(P)+1,II)*FUNKNO(IND,II) 270 CONTINUE 280 CONTINUE ! ONE-SPEED FLUX CALCULATION CALL SNFC12(LX,LY,NBMIX,NPQ,NSCT,MAT,VOL,SGAR(0,II),NCODE,ZCODE, 1 QEXT(1,II),LFIXUP,DU,DE,W,MRM,MRMY,DB,DA,DAL,FUNKNO(1,II), 2 FUNKNO(L4+1,II),FUNKNO(L4+IELEM*LY*NPQ+1,II),MN,DN) 290 CONTINUE ! END OF ENERGY GROUP ELSE IF((ITYPE.EQ.7).AND.(IBFP.EQ.0)) THEN *---- * 3D CARTESIAN CASE *---- ! EXTRACTING PARAMETERS NM=IELEM**3 NMX=IELEM**2 IF(NM*NSCT*NREG.GT.NUN) CALL XABORT('SNFLUX: QE' 1 //'XT OVERFLOW(6a)') CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) CALL LCMGPD(IPTRK,'DU',DU_PTR) CALL LCMGPD(IPTRK,'DE',DE_PTR) CALL LCMGPD(IPTRK,'DZ',DZ_PTR) CALL LCMGPD(IPTRK,'W',W_PTR) CALL LCMGPD(IPTRK,'MRMX',MRMX_PTR) CALL LCMGPD(IPTRK,'MRMY',MRMY_PTR) CALL LCMGPD(IPTRK,'MRMZ',MRMZ_PTR) CALL LCMGPD(IPTRK,'DC',DC_PTR) CALL LCMGPD(IPTRK,'DB',DB_PTR) CALL LCMGPD(IPTRK,'DA',DA_PTR) CALL LCMGPD(IPTRK,'WX',WX_PTR) CALL LCMGPD(IPTRK,'CST',CST_PTR) CALL LCMGPD(IPTRK,'IL',IL_PTR) CALL LCMGPD(IPTRK,'MN',MN_PTR) CALL LCMGPD(IPTRK,'DN',DN_PTR) CALL C_F_POINTER(DU_PTR,DU,(/ NPQ /)) CALL C_F_POINTER(DE_PTR,DE,(/ NPQ /)) CALL C_F_POINTER(DZ_PTR,DZ,(/ NPQ /)) CALL C_F_POINTER(W_PTR,W,(/ NPQ /)) CALL C_F_POINTER(MRMX_PTR,MRMX,(/ NPQ /)) CALL C_F_POINTER(MRMY_PTR,MRMY,(/ NPQ /)) CALL C_F_POINTER(MRMZ_PTR,MRMZ,(/ NPQ /)) CALL C_F_POINTER(DC_PTR,DC,(/ LX*LY*NPQ /)) CALL C_F_POINTER(DB_PTR,DB,(/ LX*LY*NPQ /)) CALL C_F_POINTER(DA_PTR,DA,(/ LX*LY*NPQ /)) CALL C_F_POINTER(WX_PTR,WX,(/ IELEM+1 /)) CALL C_F_POINTER(CST_PTR,CST,(/ IELEM /)) CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) CALL C_F_POINTER(MN_PTR,MN,(/ NPQ*NSCT /)) CALL C_F_POINTER(DN_PTR,DN,(/ NPQ*NSCT /)) CALL LCMGET(IPTRK,'NCODE',NCODE) CALL LCMGET(IPTRK,'ZCODE',ZCODE) ! LOOP FOR GROUPS DO 340 II=1,NGEFF IF(.NOT.INCONV(II)) GO TO 340 IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN/3D-car' ! FIXED VOLUMIC SOURCES QEXT(:NUN,II)=SUNKNO(:NUN,II) ! IN-GROUP SCATTERING DO 330 IR=1,NREG IBM=MAT(IR) IF(IBM.EQ.0) GO TO 330 IOF=0 DO 320 P=1,NSCT IF(IL(P).GT.NANI-1) GO TO 320 DO 300 IEL=1,NM IND=(IR-1)*NSCT*NM+(P-1)*NM+IEL QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,IL(P)+1,II)*FUNKNO(IND,II) 300 CONTINUE 320 CONTINUE 330 CONTINUE 340 CONTINUE ! END OF ENERGY LOOP ! FLUX CALCULATION IF(NKBA.EQ.0)THEN CALL SNFBC3(NUN,NGEFF,IMPX,INCONV,NGIND,LX,LY,LZ,IELEM,NM,NMX, 1 NMX,NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR,NCODE,ZCODE,QEXT,LFIXUP, 2 DU,DE,DZ,W,MRMX,MRMY,MRMZ,DC,DB,DA,MN,DN,WX,WX,WX,CST, 3 LADPT(2:4),ISBS,SUM(NBS),ISBSM,BS,MAXL,FUNKNO) ELSE IF(NKBA.GT.0)THEN IF(ISBS.EQ.1) CALL XABORT('SNFLUX: BOUNDARY SOURCES NOT YET ' 1 //'VERIFIED WITH KBA ALGORITHM. YET TO BE VALIDATED. COMMENT ' 2 //'CALL TO XABORT, RECOMPILE AND PROCEED WITH CAUTION.') IF(ISCHM.EQ.3) CALL XABORT('SNFLUX: AWD METHOD NOT YET ' 1 //'VERIFIED WITH KBA ALGORITHM. YET TO BE VALIDATED. COMMENT ' 2 //'CALL TO XABORT, RECOMPILE AND PROCEED WITH CAUTION.') CALL SNFKC3(NKBA,NUN,NGEFF,IMPX,INCONV,NGIND,LX,LY,LZ,IELEM, 1 NM,NMX,NMX,NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR,NCODE,ZCODE,QEXT, 2 LFIXUP,DU,DE,DZ,W,MRMX,MRMY,MRMZ,DC,DB,DA,MN,DN,WX,WX,WX,CST, 3 LADPT(2:4),ISBS,SUM(NBS),ISBSM,BS,MAXL,FUNKNO) ENDIF ELSE IF(ITYPE.EQ.7) THEN *------------ * 3D CARTESIAN BOLTZMANN-FOKKER-PLANCK *------------ ! EXTRACTING PARAMETERS NM=IELEM**3*EELEM NME=IELEM**3 NMX=IELEM**2*EELEM CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) IF((NSCT*NM+NPQ*NME)*NREG.GT.NUN) THEN CALL XABORT('SNFLUX: QEXT OVERFLOW(6b)') ENDIF CALL LCMGPD(IPTRK,'DU',DU_PTR) CALL LCMGPD(IPTRK,'DE',DE_PTR) CALL LCMGPD(IPTRK,'DZ',DZ_PTR) CALL LCMGPD(IPTRK,'W',W_PTR) CALL LCMGPD(IPTRK,'MRMX',MRMX_PTR) CALL LCMGPD(IPTRK,'MRMY',MRMY_PTR) CALL LCMGPD(IPTRK,'MRMZ',MRMZ_PTR) CALL LCMGPD(IPTRK,'DC',DC_PTR) CALL LCMGPD(IPTRK,'DB',DB_PTR) CALL LCMGPD(IPTRK,'DA',DA_PTR) CALL LCMGPD(IPTRK,'WX',WX_PTR) CALL LCMGPD(IPTRK,'WE',WE_PTR) CALL LCMGPD(IPTRK,'CST',CST_PTR) CALL LCMGPD(IPTRK,'IL',IL_PTR) CALL LCMGPD(IPTRK,'MN',MN_PTR) CALL LCMGPD(IPTRK,'DN',DN_PTR) CALL C_F_POINTER(DU_PTR,DU,(/ NPQ /)) CALL C_F_POINTER(DE_PTR,DE,(/ NPQ /)) CALL C_F_POINTER(DZ_PTR,DZ,(/ NPQ /)) CALL C_F_POINTER(W_PTR,W,(/ NPQ /)) CALL C_F_POINTER(MRMX_PTR,MRMX,(/ NPQ /)) CALL C_F_POINTER(MRMY_PTR,MRMY,(/ NPQ /)) CALL C_F_POINTER(MRMZ_PTR,MRMZ,(/ NPQ /)) CALL C_F_POINTER(DC_PTR,DC,(/ LX*LY*NPQ /)) CALL C_F_POINTER(DB_PTR,DB,(/ LX*LY*NPQ /)) CALL C_F_POINTER(DA_PTR,DA,(/ LX*LY*NPQ /)) CALL C_F_POINTER(WX_PTR,WX,(/ IELEM+1 /)) CALL C_F_POINTER(WE_PTR,WE,(/ EELEM+1 /)) CALL C_F_POINTER(CST_PTR,CST,(/ MAX(IELEM,EELEM) /)) CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) CALL C_F_POINTER(MN_PTR,MN,(/ NPQ*NSCT /)) CALL C_F_POINTER(DN_PTR,DN,(/ NPQ*NSCT /)) CALL LCMGET(IPTRK,'NCODE',NCODE) CALL LCMGET(IPTRK,'ZCODE',ZCODE) ! ANGULAR FOKKER-PLANCK OPERATOR ALLOCATE(SIGMATR(NSCT,NSCT),MTR(NPQ,NPQ)) SIGMATR=0.0 DO P=1,NSCT SIGMATR(P,P)=-IL(P)*(IL(P)+1) ENDDO MTR=MATMUL(MATMUL(RESHAPE(MN,[NPQ, NSCT]),SIGMATR), 1 RESHAPE(DN,[NSCT,NPQ])) LAMBDA0=0.0 DO M=1,NPQ IF(MTR(M,M).LT.-LAMBDA0) LAMBDA0=-MTR(M,M) ENDDO DO M=1,NPQ MTR(M,M)=MTR(M,M)+LAMBDA0 ENDDO SIGMATR=MATMUL(MATMUL(RESHAPE(DN,[NSCT,NPQ]),MTR), 1 RESHAPE(MN,[NPQ,NSCT])) DEALLOCATE(MTR) DO II=1,NGEFF DO IBM=1,NBMIX SGAR(IBM,II)=SGAR(IBM,II)+LAMBDA0*EMOMTR(IBM,II) ENDDO ENDDO ! LOOP FOR GROUPS DO 390 II=1,NGEFF IF(.NOT.INCONV(II)) GO TO 390 IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN-BFP/3D-car' ! FIXED VOLUMIC SOURCES QEXT(:NUN,II)=SUNKNO(:NUN,II) ! IN-GROUP SCATTERING DO 380 IR=1,NREG IBM=MAT(IR) IF(IBM.EQ.0) GO TO 380 IOF=0 DO 370 P=1,NSCT IF(IL(P).GT.NANI-1) GO TO 370 DO 350 IEL=1,NM IND=(IR-1)*NSCT*NM+(P-1)*NM+IEL QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,IL(P)+1,II)*FUNKNO(IND,II) DO Q=1,NSCT JND=(IR-1)*NSCT*NM+(Q-1)*NM+IEL QEXT(IND,II)=QEXT(IND,II)+EMOMTR(IBM,II)*SIGMATR(P,Q) 1 *FUNKNO(JND,II) ENDDO 350 CONTINUE 370 CONTINUE 380 CONTINUE 390 CONTINUE ! END OF ENERGY LOOP DEALLOCATE(SIGMATR) ! FLUX CALCULATION CALL SNFE3D(NUN,NGEFF,IMPX,INCONV,NGIND,LX,LY,LZ, 1 IELEM,EELEM,NM,NME,NMX,NMX,NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR, 2 ESTOPW,NCODE,ZCODE,DELTAE,QEXT,LFIXUP,DU,DE,DZ,W,MRMX,MRMY,MRMZ, 3 DC,DB,DA,FUNKNO,ISLG,FLUXC,ISBS,SUM(NBS),ISBSM,BS,MAXL, 4 WX,WX,WX,WE,CST,LADPT,IBFP,MN,DN) ELSE IF(ITYPE.EQ.8) THEN *------------ * 2D HEXAGONAL *------------ ! EXTRACTING PARAMETERS NM=IELEM**2 NMX=IELEM IF(NM*NSCT*NREG.GT.NUN) CALL XABORT('SNFLUX: QEXT OVE' 1 //'RFLOW(7)') CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) CALL LCMGPD(IPTRK,'DU',DU_PTR) CALL LCMGPD(IPTRK,'DE',DE_PTR) CALL LCMGPD(IPTRK,'W',W_PTR) CALL LCMGPD(IPTRK,'DB',DB_PTR) CALL LCMGPD(IPTRK,'DA',DA_PTR) CALL LCMGPD(IPTRK,'WX',WX_PTR) CALL LCMGPD(IPTRK,'CST',CST_PTR) CALL LCMGPD(IPTRK,'IL',IL_PTR) CALL LCMGPD(IPTRK,'MN',MN_PTR) CALL LCMGPD(IPTRK,'DN',DN_PTR) CALL C_F_POINTER(DU_PTR,DU,(/ NPQ /)) CALL C_F_POINTER(DE_PTR,DE,(/ NPQ /)) CALL C_F_POINTER(W_PTR,W,(/ NPQ /)) CALL C_F_POINTER(DB_PTR,DB,(/ LX*NPQ /)) CALL C_F_POINTER(DA_PTR,DA,(/ LX*LY*NPQ /)) CALL C_F_POINTER(WX_PTR,WX,(/ IELEM+1 /)) CALL C_F_POINTER(CST_PTR,CST,(/ IELEM /)) CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) CALL C_F_POINTER(MN_PTR,MN,(/ NPQ*NSCT /)) CALL C_F_POINTER(DN_PTR,DN,(/ NPQ*NSCT /)) CALL LCMGET(IPTRK,'SIDE',SIDE) CALL LCMGET(IPTRK,'LOZSWP',LOZSWP) NHEX=LX/(3*ISPLH**2) ALLOCATE(COORDMAP(3,NHEX)) CALL LCMGET(IPTRK,'COORDMAP',COORDMAP) ! LOOP FOR GROUPS DO 440 II=1,NGEFF IF(.NOT.INCONV(II)) GO TO 440 IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN/2D-hex' ! FIXED VOLUMIC SOURCES QEXT(:NUN,II)=SUNKNO(:NUN,II) ! IN-GROUP SCATTERING DO 430 IR=1,NREG IBM=MAT(IR) IF(IBM.EQ.0) GO TO 430 IOF=0 DO 420 P=1,NSCT IF(IL(P).GT.NANI-1) GO TO 420 DO 400 IEL=1,NM IND=(((IR-1)*NSCT+(P-1))*NM)+IEL QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,IL(P)+1,II)*FUNKNO(IND,II) 400 CONTINUE 420 CONTINUE 430 CONTINUE 440 CONTINUE ! END OF ENERGY LOOP ! FLUX CALCULATION IF(NKBA.EQ.0)THEN CALL SNFBH2(NUN,NGEFF,IMPX,INCONV,NGIND,NHEX,ISPLH,SIDE, 1 IELEM,NM,NMX,NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR,QEXT,LFIXUP,DU, 2 DE,W,DB,DA,MN,DN,WX,WX,CST,LADPT(2:3),LOZSWP,COORDMAP, 3 FUNKNO) ELSE IF(NKBA.GE.1)THEN IF(ISCHM.EQ.3) CALL XABORT('SNFLUX: AWD METHOD NOT YET ' 1 //'VERIFIED WITH KBA ALGORITHM. YET TO BE VALIDATED. COMMENT ' 2 //'CALL TO XABORT, RECOMPILE AND PROCEED WITH CAUTION.') CALL SNFKH2(NUN,NGEFF,IMPX,INCONV,NGIND,NHEX,ISPLH,SIDE, 1 IELEM,NM,NMX,NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR,QEXT,LFIXUP,DU, 2 DE,W,DB,DA,MN,DN,WX,WX,CST,LADPT(2:3),LOZSWP,COORDMAP, 3 FUNKNO) ENDIF DEALLOCATE(COORDMAP) ELSE IF(ITYPE.EQ.9) THEN *------------ * 3D HEXAGONAL *------------ ! EXTRACTING PARAMETERS NM=IELEM**3 NMX=IELEM**2 IF(IELEM*IELEM*IELEM*NSCT*NREG.GT.NUN) CALL XABORT('SNFLUX: QE' 1 //'XT OVERFLOW(8)') CALL LCMLEN(IPTRK,'DU',NPQ,ITYLCM) CALL LCMGPD(IPTRK,'DU',DU_PTR) CALL LCMGPD(IPTRK,'DE',DE_PTR) CALL LCMGPD(IPTRK,'DZ',DZ_PTR) CALL LCMGPD(IPTRK,'W',W_PTR) CALL LCMGPD(IPTRK,'MRMZ',MRMZ_PTR) CALL LCMGPD(IPTRK,'DC',DC_PTR) CALL LCMGPD(IPTRK,'DB',DB_PTR) CALL LCMGPD(IPTRK,'DA',DA_PTR) CALL LCMGPD(IPTRK,'WX',WX_PTR) CALL LCMGPD(IPTRK,'CST',CST_PTR) CALL LCMGPD(IPTRK,'IL',IL_PTR) CALL LCMGPD(IPTRK,'MN',MN_PTR) CALL LCMGPD(IPTRK,'DN',DN_PTR) CALL C_F_POINTER(DU_PTR,DU,(/ NPQ /)) CALL C_F_POINTER(DE_PTR,DE,(/ NPQ /)) CALL C_F_POINTER(DZ_PTR,DZ,(/ NPQ /)) CALL C_F_POINTER(W_PTR,W,(/ NPQ /)) CALL C_F_POINTER(MRMZ_PTR,MRMZ,(/ NPQ /)) CALL C_F_POINTER(DC_PTR,DC,(/ LX*LY*NPQ /)) CALL C_F_POINTER(DB_PTR,DB,(/ LX*LY*NPQ /)) CALL C_F_POINTER(DA_PTR,DA,(/ LX*LY*NPQ /)) CALL C_F_POINTER(IL_PTR,IL,(/ NSCT /)) CALL C_F_POINTER(WX_PTR,WX,(/ IELEM+1 /)) CALL C_F_POINTER(CST_PTR,CST,(/ IELEM /)) CALL C_F_POINTER(MN_PTR,MN,(/ NPQ*NSCT /)) CALL C_F_POINTER(DN_PTR,DN,(/ NPQ*NSCT /)) CALL LCMGET(IPTRK,'NCODE',NCODE) CALL LCMGET(IPTRK,'ZCODE',ZCODE) CALL LCMGET(IPTRK,'SIDE',SIDE) CALL LCMGET(IPTRK,'LOZSWP',LOZSWP) ! Number of hexagons in one plane only NHEX=LX/(3*ISPLH**2) ALLOCATE(COORDMAP(3,NHEX)) CALL LCMGET(IPTRK,'COORDMAP',COORDMAP) ! LOOP FOR GROUPS DO 490 II=1,NGEFF IF(.NOT.INCONV(II)) GO TO 490 IF(IMPX.GT.2) WRITE(IUNOUT,510) NGIND(II),'SN/3D-hex' ! FIXED VOLUMIC SOURCES QEXT(:NUN,II)=SUNKNO(:NUN,II) ! IN-GROUP SCATTERING DO 480 IR=1,NREG IBM=MAT(IR) IF(IBM.EQ.0) GO TO 480 IOF=0 DO 470 P=1,NSCT IF(IL(P).GT.NANI-1) GO TO 470 DO 450 IEL=1,IELEM*IELEM*IELEM IND=(((IR-1)*NSCT+(P-1))*IELEM**3)+IEL QEXT(IND,II)=QEXT(IND,II)+SGAS(IBM,IL(P)+1,II)*FUNKNO(IND,II) 450 CONTINUE 470 CONTINUE 480 CONTINUE 490 CONTINUE ! END OF ENERGY LOOP ! FLUX CALCULATION IF(NKBA.EQ.0)THEN CALL SNFBH3(NUN,NGEFF,IMPX,INCONV,NGIND,NHEX,LZ,ISPLH,SIDE, 1 IELEM,NM,NMX,NMX,NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR,NCODE,ZCODE, 2 QEXT,LFIXUP,DU,DE,DZ,W,MRMZ,DC,DB,DA,MN,DN,WX,WX,WX,CST, 3 LADPT(2:4),LOZSWP,COORDMAP,FUNKNO) ELSE IF(NKBA.GE.1)THEN IF(ISCHM.EQ.3) CALL XABORT('SNFLUX: AWD METHOD NOT YET ' 1 //'VERIFIED WITH KBA ALGORITHM. YET TO BE VALIDATED. COMMENT ' 2 //'CALL TO XABORT, RECOMPILE AND PROCEED WITH CAUTION.') CALL SNFKH3(NKBA,NUN,NGEFF,IMPX,INCONV,NGIND,NHEX,LZ,ISPLH, 1 SIDE,IELEM,NM,NMX,NMX,NMX,NBMIX,NPQ,NSCT,MAT,VOL,SGAR,NCODE, 2 ZCODE,QEXT,LFIXUP,DU,DE,DZ,W,MRMZ,DC,DB,DA,MN,DN,WX,WX,WX, 3 CST,LADPT(2:4),LOZSWP,COORDMAP,FUNKNO) ENDIF DEALLOCATE(COORDMAP) ELSE CALL XABORT('SNFLUX: TYPE OF DISCRETIZATION NOT IMPLEMENTED.') ENDIF DEALLOCATE(ESTOPW,SGAS,SGAR) *---- * PRINT COMPLETE UNKNOWN VECTOR. *---- DO 495 II=1,NGEFF IF(.NOT.INCONV(II)) GO TO 495 IF(IMPX.GT.5) THEN WRITE(IUNOUT,520) NGIND(II) WRITE(IUNOUT,'(1P,4(5X,E15.7))') (FUNKNO(:,II)) ENDIF 495 CONTINUE *---- * DIFFUSION SYNTHETIC ACCELERATION. *---- IF(LDSA) THEN ISOLVSA=IPAR(33) CALL LCMSIX(IPTRK,'DSA',1) CALL LCMGET(IPTRK,'KEYFLX',KEYSPN) CALL LCMGET(IPTRK,'STATE-VECTOR',IPAR) IF(NREG.NE.IPAR(1)) CALL XABORT('SNFLUX: INVALID NREG (DSA).') NUNSA=IPAR(2) ITYPE=IPAR(6) IF((ITYPE.EQ.7).OR.(ITYPE.EQ.9)) THEN IELEMSA=IPAR(9) ELSE IF(ISOLVSA.EQ.1)THEN IELEMSA=IPAR(8) ELSEIF(ISOLVSA.EQ.2)THEN IELEMSA=IPAR(9) ENDIF ENDIF IMPY=MAX(0,IMPX-1) CALL LCMSIX(IPTRK,' ',2) CALL SNDSA(KPSYS,INCONV,NGIND,IPTRK,IMPY,NGEFF,NREG,NBMIX, 1 NUN,MAT,VOL,KEYFLX,KEYSPN,NUNSA,IELEMSA,ZCODE,OLD,FUNKNO, 2 NHEX) DEALLOCATE(OLD) ENDIF *---- * SCRATCH STORAGE DEALLOCATION *---- DEALLOCATE(DELTAE,KEYSPN,QEXT,BS,ISBSM) RETURN * 500 FORMAT(//41H SNFLUX: N E U T R O N S O U R C E S (,I5,3H ):) 510 FORMAT(/25H SNFLUX: PROCESSING GROUP,I5,6H WITH ,A,1H.) 520 FORMAT(//41H SNFLUX: A F T E R T R A N S P O R T (,I5,3H ):) END