*DECK SNFE2D SUBROUTINE SNFE2D(NUN,NGEFF,IMPX,INCONV,NGIND,LX,LY,IELEM, 1 EELEM,NM,NME,NMX,NMY,NMAT,NPQ,NSCT,MAT,VOL,TOTAL,ESTOPW, 2 NCODE,ZCODE,DELTAE,QEXT,LFIXUP,DU,DE,W,MRM,MRMY,DB,DA,FUNKNO, 3 ISLG,FLUXC,ISBS,NBS,ISBSM,BS,MAXL,WX,WY,WE,CST,ISADPT,IBFP,MN,DN) * *----------------------------------------------------------------------- * *Purpose: * Perform one inner iteration for solving SN equations in 2D Cartesian * geometry for the HODD method. Energy-angle multithreading. Albedo * boundary conditions. Boltzmann-Fokker-Planck (BFP) discretization. * *Copyright: * Copyright (C) 2021 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, A. A. Calloo and C. Bienvenue * *Parameters: input * NUN total number of unknowns in vector FUNKNO. * NGEFF number of energy groups processed in parallel. * IMPX print flag (equal to zero for no print). * INCONV energy group convergence flag (set to .FALSE. if converged). * NGIND energy group indices assign to the NGEFF set. * LX number of meshes along X axis. * LY number of meshes along Y axis. * IELEM measure of order of the spatial approximation polynomial: * =1 constant - default for HODD; * =2 linear - default for DG; * >3 higher orders. * EELEM measure of order of the energy approximation polynomial: * =1 constant - default for HODD; * =2 linear - default for DG; * >3 higher orders. * NM number of moments in space and energy for flux components * NME number of moments for energy boundaries components * NMX number of moments for X axis boundaries components * NMY number of moments for Y axis boundaries components * NMAT number of material mixtures. * NPQ number of SN directions in four octants (including zero-weight * directions). * NSCT maximum number of spherical harmonics moments of the flux. * MAT material mixture index in each region. * VOL volumes of each region. * TOTAL macroscopic total cross sections. * ESTOPW stopping power. * NCODE boundary condition indices. * ZCODE albedos. * DELTAE energy group width in MeV. * QEXT Legendre components of the fixed source. * LFIXUP flag to enable negative flux fixup. * DU first direction cosines ($\\mu$). * DE second direction cosines ($\\eta$). * W weights. * MRM quadrature index. * MRMY quadrature index. * DB diamond-scheme parameter. * DA diamond-scheme parameter. * MN moment-to-discrete matrix. * DN discrete-to-moment matrix. * ISBS flag to indicate the presence or not of boundary fixed * sources. * NBS number of boundary fixed sources. * ISBSM flag array to indicate the presence or not of boundary fixed * source in each unit surface. * BS boundary source array with their intensities. * MAXL maximum size of boundary source array. * WX spatial X axis closure relation weighting factors. * WY spatial Y axis closure relation weighting factors. * WE energy closure relation weighting factors. * CST constants for the polynomial approximations. * ISADPTX flag to enable/disable adaptive X axis flux calculations. * ISADPTY flag to enable/disable adaptive Y axis flux calculations. * ISADPTE flag to enable/disable adaptive energy flux calculations. * IBFP type of energy proparation relation: * =1 Galerkin type; * =2 heuristic Przybylski and Ligou type. * *Parameters: input/output * FUNKNO Legendre components of the flux and boundary fluxes. * FLUXC flux at the cutoff energy. * *----------------------------------------------------------------------- #if defined(_OPENMP) USE omp_lib #endif * *---- * SUBROUTINE ARGUMENTS *---- INTEGER NUN,NGEFF,IMPX,NGIND(NGEFF),LX,LY,IELEM,EELEM, 1 NM,NME,NMX,NMY,NMAT,NPQ, 2 NSCT,MAT(LX,LY),NCODE(4),MRM(NPQ),MRMY(NPQ),ISLG(NGEFF),ISBS, 3 NBS,ISBSM(4*ISBS,NPQ*ISBS,NGEFF*ISBS),MAXL LOGICAL INCONV(NGEFF) REAL VOL(LX,LY),TOTAL(0:NMAT,NGEFF),ESTOPW(0:NMAT,2,NGEFF), 1 ZCODE(4),DELTAE(NGEFF),QEXT(NUN,NGEFF),DU(NPQ),DE(NPQ),W(NPQ), 2 DB(LX,NPQ),DA(LX,LY,NPQ),FUNKNO(NUN,NGEFF),FLUXC(LX,LY), 3 BS(MAXL*ISBS,NBS*ISBS),WX(IELEM+1),WY(IELEM+1),WE(EELEM+1), 4 CST(MAX(IELEM,EELEM)),MN(NPQ,NSCT),DN(NSCT,NPQ) LOGICAL LFIXUP,ISADPT(3) *---- * LOCAL VARIABLES *---- INTEGER NPQD(4),IIND(4),P REAL BM,BP,TB,WX0(IELEM+1),WY0(IELEM+1),WE0(EELEM+1) DOUBLE PRECISION Q(NM),Q2(NM,NM+1),FEP(NME), 1 XNJ(NMY),V PARAMETER(IUNOUT=6,RLOG=1.0E-8,PI=3.141592654) LOGICAL ISFIX(3) *---- * ALLOCATABLE ARRAYS *---- INTEGER, ALLOCATABLE, DIMENSION(:,:) :: INDANG DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:,:) :: FLUX,FLUX0 DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:,:,:) :: FLUX_G, 1 FLUX0_G DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: XNI *---- * SCRATCH STORAGE ALLOCATION *---- ALLOCATE(INDANG(NPQ,4)) ALLOCATE(XNI(NMX,LY),FLUX(NM,NSCT,LX,LY), 1 FLUX0(NME,NPQ,LX,LY)) ALLOCATE(FLUX_G(NM,NSCT,LX,LY,NGEFF), 1 FLUX0_G(NME,NPQ,LX,LY,NGEFF)) *---- * LENGTH OF FUNKNO COMPONENTS (IN ORDER) *---- LFLX=NM*LX*LY*NSCT LXNI=NMX*LY*NPQ LXNJ=NMY*LX*NPQ LFEP=NME*LX*LY*NPQ *---- * SET OCTANT SWAPPING ORDER. *---- NPQD(:4)=0 INDANG(:NPQ,:4)=0 DO M=1,NPQ VU=DU(M) VE=DE(M) IF((VU.GE.0.0).AND.(VE.GE.0.0)) THEN IND=1 JND=4 ELSE IF((VU.LE.0.0).AND.(VE.GE.0.0)) THEN IND=2 JND=3 ELSE IF((VU.LE.0.0).AND.(VE.LE.0.0)) THEN IND=3 JND=1 ELSE IND=4 JND=2 ENDIF IIND(JND)=IND NPQD(IND)=NPQD(IND)+1 INDANG(NPQD(IND),IND)=M ENDDO *---- * MAIN LOOP OVER OCTANTS. *---- FLUX_G(:NM,:NSCT,:LX,:LY,:NGEFF)=0.0D0 FLUX0_G(:NME,:NPQ,:LX,:LY,:NGEFF)=0.0D0 WE0=WE WX0=WX WY0=WY DO 190 JND=1,4 IND=IIND(JND) *---- * PRELIMINARY LOOPS FOR SETTING BOUNDARY CONDITIONS. *---- *$OMP PARALLEL DO *$OMP+ PRIVATE(M,IG,VU,VE,M1,IOF,JOF,IEL,I,J,IPQD) *$OMP+ SHARED(FUNKNO) COLLAPSE(2) DO 70 IG=1,NGEFF DO 60 IPQD=1,NPQD(IND) IF(.NOT.INCONV(IG)) GO TO 60 M=INDANG(IPQD,IND) VU=DU(M) VE=DE(M) ! X-BOUNDARY IF(VU.GT.0.0)THEN M1=MRM(M) IF((NCODE(1).NE.4))THEN DO IEL=1,NMX DO J=1,LY IOF=((M-1)*LY+(J-1))*NMX+IEL JOF=((M1-1)*LY+(J-1))*NMX+IEL FUNKNO(LFLX+IOF,IG)=FUNKNO(LFLX+JOF,IG) ENDDO ENDDO ENDIF ELSEIF(VU.LT.0.0)THEN M1=MRM(M) IF((NCODE(2).NE.4))THEN DO IEL=1,NMX DO J=1,LY IOF=((M-1)*LY+(J-1))*NMX+IEL JOF=((M1-1)*LY+(J-1))*NMX+IEL FUNKNO(LFLX+IOF,IG)=FUNKNO(LFLX+JOF,IG) ENDDO ENDDO ENDIF ENDIF ! Y-BOUNDARY IF(VE.GT.0.0)THEN M1=MRMY(M) IF((NCODE(3).NE.4))THEN DO IEL=1,NMY DO I=1,LX IOF=((M-1)*LX+(I-1))*NMY+IEL JOF=((M1-1)*LX+(I-1))*NMY+IEL FUNKNO(LFLX+LXNI+IOF,IG)= > FUNKNO(LFLX+LXNI+JOF,IG) ENDDO ENDDO ENDIF ELSEIF(VE.LT.0.0)THEN M1=MRMY(M) IF((NCODE(4).NE.4))THEN DO IEL=1,NMY DO I=1,LX IOF=((M-1)*LX+(I-1))*NMY+IEL JOF=((M1-1)*LX+(I-1))*NMY+IEL FUNKNO(LFLX+LXNI+IOF,IG)= > FUNKNO(LFLX+LXNI+JOF,IG) ENDDO ENDDO ENDIF ENDIF 60 CONTINUE 70 CONTINUE *$OMP END PARALLEL DO *---- * MAIN SWAPPING LOOPS FOR SN FLUX CALCULATION *---- *$OMP PARALLEL DO *$OMP+ PRIVATE(ITID,FLUX,M,IG,XNI,XNJ,Q,Q2,IOF,IER,II,JJ,IEL,I,J,L) *$OMP+ PRIVATE(FEP,FLUX0,BM,BP,IIE,IIX,IIY,IE,IX,IY) *$OMP+ PRIVATE(ISFIX,JX,JY,JE,TB,V,SIGMA,IBM,J0,I0,IPQD) *$OMP+ FIRSTPRIVATE(WE,WX,WY,WE0,WX0,WY0) SHARED(FUNKNO) *$OMP+ REDUCTION(+:FLUX_G,FLUX0_G,FLUXC) COLLAPSE(2) ! LOOP FOR GROUPS TO EXECUTE IN PARALLEL DO 180 IG=1,NGEFF ! LOOP OVER ALL DIRECTIONS DO 170 IPQD=1,NPQD(IND) IF(.NOT.INCONV(IG)) GO TO 170 M=INDANG(IPQD,IND) IF(W(M).EQ.0.0) GO TO 170 ! GET AND PRINT THREAD NUMBER #if defined(_OPENMP) ITID=omp_get_thread_num() #else ITID=0 #endif IF(IMPX.GT.5) WRITE(IUNOUT,400) ITID,NGIND(IG),IPQD ! INITIALIZE FLUXES FLUX(:NM,:NSCT,:LX,:LY)=0.0D0 FLUX0(:NME,:NPQ,:LX,:LY)=0.0D0 *---- * LOOP OVER X- AND Y-DIRECTED AXES. *---- ! X-AXIS LOOP DO 155 I0=1,LX I=I0 IF((IND.EQ.2).OR.(IND.EQ.3)) I=LX+1-I ! Y-BOUNDARIES CONDITIONS XNJ=0.0 DO IEL=1,NMY IOF=(M-1)*NMY*LX+(I-1)*NMY+IEL IF((IND.EQ.1).OR.(IND.EQ.2)) THEN XNJ(IEL)=FUNKNO(LFLX+LXNI+IOF,IG)*ZCODE(3) ELSE XNJ(IEL)=FUNKNO(LFLX+LXNI+IOF,IG)*ZCODE(4) ENDIF ENDDO ! Y-BOUNDARIES FIXED SOURCES IF(ISBS.EQ.1) THEN IF((IND.EQ.3.OR.IND.EQ.4).AND.ISBSM(4,M,IG).NE.0) THEN XNJ(1)=XNJ(1)+BS(I,ISBSM(4,M,IG)) ELSE IF((IND.EQ.1.OR.IND.EQ.2).AND.ISBSM(3,M,IG).NE.0) THEN XNJ(1)=XNJ(1)+BS(I,ISBSM(3,M,IG)) ENDIF ENDIF ! Y-AXIS LOOP DO 140 J0=1,LY J=J0 IF((IND.EQ.3).OR.(IND.EQ.4)) J=LY+1-J ! X-BOUNDARIES CONDITIONS IF(I0.EQ.1) THEN XNI(:NMX,J)=0.0 DO IEL=1,NMX IOF=(M-1)*NMX*LY+(J-1)*NMX+IEL IF((IND.EQ.1).OR.(IND.EQ.4)) THEN XNI(IEL,J)=FUNKNO(LFLX+IOF,IG)*ZCODE(1) ELSE XNI(IEL,J)=FUNKNO(LFLX+IOF,IG)*ZCODE(2) ENDIF ENDDO ENDIF ! X-BOUNDARIES FIXED SOURCES IF(ISBS.EQ.1.AND.I0.EQ.1) THEN IF((IND.EQ.2.OR.IND.EQ.3).AND.ISBSM(2,M,IG).NE.0) THEN XNI(1,J)=XNI(1,J)+BS(J,ISBSM(2,M,IG)) ELSE IF((IND.EQ.1.OR.IND.EQ.4).AND.ISBSM(1,M,IG).NE.0) THEN XNI(1,J)=XNI(1,J)+BS(J,ISBSM(1,M,IG)) ENDIF ENDIF ! DATA IBM=MAT(I,J) IF(IBM.EQ.0) GO TO 140 SIGMA=TOTAL(IBM,IG) BM=ESTOPW(IBM,1,IG)/DELTAE(IG) BP=ESTOPW(IBM,2,IG)/DELTAE(IG) V=VOL(I,J) ! TYPE OF ENERGY PROPAGATION FACTOR IF(IBFP.EQ.1) THEN ! GALERKIN TYPE TB=BM/BP WE(1)=WE(1)*TB WE(2:EELEM+1)=(WE(2:EELEM+1)-1)*TB+1 ELSE ! PRZYBYLSKI AND LIGOU TYPE TB=1.0 ENDIF ! SOURCE DENSITY TERM DO IEL=1,NM Q(IEL)=0.0D0 DO P=1,NSCT IOF=((J-1)*LX*NSCT+(I-1)*NSCT+(P-1))*NM+IEL Q(IEL)=Q(IEL)+QEXT(IOF,IG)*MN(M,P) ENDDO ENDDO ! ENERGY GROUP UPPER BOUNDARY INCIDENT FLUX DO IEL=1,NME IOF=((J-1)*LX*NPQ+(I-1)*NPQ+(M-1))*NME+IEL FEP(IEL)=QEXT(LFLX+LXNI+LXNJ+IOF,IG) ENDDO ISFIX=.FALSE. DO WHILE (.NOT.ALL(ISFIX)) ! LOOP FOR ADAPTIVE CALCULATION ! FLUX MOMENT COEFFICIENTS MATRIX Q2(:NM,:NM+1)=0.0D0 DO IY=1,IELEM DO JY=1,IELEM DO IX=1,IELEM DO JX=1,IELEM DO IE=1,EELEM DO JE=1,EELEM II=IELEM*EELEM*(IY-1)+EELEM*(IX-1)+IE JJ=IELEM*EELEM*(JY-1)+EELEM*(JX-1)+JE ! DIAGONAL TERMS IF(II.EQ.JJ) THEN Q2(II,JJ)=(SIGMA+CST(IE)**2*WE(JE+1)*BP+(IE-1)*(BM-BP))*V 1 +CST(IX)**2*WX(JX+1)*ABS(DA(I,J,M)) 2 +CST(IY)**2*WY(JY+1)*ABS(DB(I,M)) ! UPPER DIAGONAL TERMS ELSEIF(II.LT.JJ) THEN IF(IY.EQ.JY) THEN ! ENERGY TERMS IF(IX.EQ.JX) THEN IF(MOD(IE+JE,2).EQ.1) THEN Q2(II,JJ)=-CST(IE)*CST(JE)*WE(JE+1)*BP*V ELSE Q2(II,JJ)=CST(IE)*CST(JE)*WE(JE+1)*BP*V ENDIF ! X-SPACE TERMS ELSEIF(IE.EQ.JE) THEN IF(MOD(IX+JX,2).EQ.1) THEN Q2(II,JJ)=CST(IX)*CST(JX)*WX(JX+1)*DA(I,J,M) ELSE Q2(II,JJ)=CST(IX)*CST(JX)*WX(JX+1)*ABS(DA(I,J,M)) ENDIF ENDIF ! Y-SPACE TERMS ELSEIF(IX.EQ.JX.AND.IE.EQ.JE) THEN IF(MOD(IY+JY,2).EQ.1) THEN Q2(II,JJ)=CST(IY)*CST(JY)*WY(JY+1)*DB(I,M) ELSE Q2(II,JJ)=CST(IY)*CST(JY)*WY(JY+1)*ABS(DB(I,M)) ENDIF ENDIF ! UNDER DIAGONAL TERMS ELSE IF(IY.EQ.JY) THEN ! ENERGY TERMS IF(IX.EQ.JX) THEN IF(MOD(IE+JE,2).EQ.1) THEN Q2(II,JJ)=-CST(IE)*CST(JE)*(WE(JE+1)*BP-BM-BP)*V ELSE Q2(II,JJ)=CST(IE)*CST(JE)*(WE(JE+1)*BP+BM-BP)*V ENDIF ! X-SPACE TERMS ELSEIF(IE.EQ.JE) THEN IF(MOD(IX+JX,2).EQ.1) THEN Q2(II,JJ)=CST(IX)*CST(JX)*(WX(JX+1)-2.0D0)*DA(I,J,M) ELSE Q2(II,JJ)=CST(IX)*CST(JX)*WX(JX+1)*ABS(DA(I,J,M)) ENDIF ENDIF ! Y-SPACE TERMS ELSEIF(IX.EQ.JX.AND.IE.EQ.JE) THEN IF(MOD(IY+JY,2).EQ.1) THEN Q2(II,JJ)=CST(IY)*CST(JY)*(WY(JY+1)-2.0D0)*DB(I,M) ELSE Q2(II,JJ)=CST(IY)*CST(JY)*WY(JY+1)*ABS(DB(I,M)) ENDIF ENDIF ENDIF ENDDO ENDDO ENDDO ENDDO ENDDO ENDDO ! FLUX SOURCE VECTOR DO IY=1,IELEM DO IX=1,IELEM DO IE=1,EELEM II=IELEM*EELEM*(IY-1)+EELEM*(IX-1)+IE IIE=IELEM*(IY-1)+IX IIX=EELEM*(IY-1)+IE IIY=EELEM*(IX-1)+IE Q2(II,NM+1)=Q(II)*V ! ENERGY TERMS IF(MOD(IE,2).EQ.1) THEN Q2(II,NM+1)=Q2(II,NM+1)+CST(IE)*(BM-WE(1)*BP)*FEP(IIE)*V ELSE Q2(II,NM+1)=Q2(II,NM+1)+CST(IE)*(BM+WE(1)*BP)*FEP(IIE)*V ENDIF ! X-SPACE TERMS IF(MOD(IX,2).EQ.1) THEN Q2(II,NM+1)=Q2(II,NM+1)+CST(IX)*(1-WX(1)) 1 *XNI(IIX,J)*ABS(DA(I,J,M)) ELSE Q2(II,NM+1)=Q2(II,NM+1)-CST(IX)*(1+WX(1)) 1 *XNI(IIX,J)*DA(I,J,M) ENDIF ! Y-SPACE TERMS IF(MOD(IY,2).EQ.1) THEN Q2(II,NM+1)=Q2(II,NM+1)+CST(IY)*(1-WY(1)) 1 *XNJ(IIY)*ABS(DB(I,M)) ELSE Q2(II,NM+1)=Q2(II,NM+1)-CST(IY)*(1+WY(1)) 1 *XNJ(IIY)*DB(I,M) ENDIF ENDDO ENDDO ENDDO CALL ALSBD(NM,1,Q2,IER,NM) IF(IER.NE.0) CALL XABORT('SNFE2D: SINGULAR MATRIX.') ! ADAPTIVE CORRECTION OF WEIGHTING PARAMETERS IF(ANY(ISADPT)) THEN IF(ISADPT(1)) THEN CALL SNADPT(EELEM,NM,IELEM**2,Q2(1:EELEM:1,NM+1), 1 FEP,TB,WE,ISFIX(1)) ELSE ISFIX(1)=.TRUE. ENDIF IF(ISADPT(2)) THEN CALL SNADPT(IELEM,NM,EELEM*IELEM,Q2(1:IELEM*EELEM:IELEM,NM+1), 1 XNI(:NMX,J),1.0,WX,ISFIX(2)) ELSE ISFIX(2)=.TRUE. ENDIF IF(ISADPT(3)) THEN CALL SNADPT(IELEM,NM,EELEM*IELEM,Q2(1:NM:IELEM*EELEM,NM+1), 1 XNJ,1.0,WY,ISFIX(3)) ELSE ISFIX(3)=.TRUE. ENDIF ELSE ISFIX=.TRUE. ENDIF END DO ! END OF ADAPTIVE LOOP ! CLOSURE RELATIONS IF(IELEM.EQ.1.AND.LFIXUP.AND.(Q2(1,2).LE.RLOG)) Q2(1,2)=0.0 XNI(:NMX,J)=WX(1)*XNI(:NMX,J) XNJ(:NMY)=WY(1)*XNJ(:NMY) FEP(:NME)=WE(1)*FEP(:NME) DO IY=1,IELEM DO IX=1,IELEM DO IE=1,EELEM II=IELEM*EELEM*(IY-1)+EELEM*(IX-1)+IE IIE=IELEM*(IY-1)+IX IIX=EELEM*(IY-1)+IE IIY=EELEM*(IX-1)+IE ! ENERGY IF(MOD(IE,2).EQ.1) THEN FEP(IIE)=FEP(IIE)+CST(IE)*WE(IE+1)*Q2(II,NM+1) ELSE FEP(IIE)=FEP(IIE)-CST(IE)*WE(IE+1)*Q2(II,NM+1) ENDIF ! X-SPACE IF(MOD(IX,2).EQ.1) THEN XNI(IIX,J)=XNI(IIX,J)+CST(IX)*WX(IX+1) 1 *Q2(II,NM+1) ELSE XNI(IIX,J)=XNI(IIX,J)+CST(IX)*WX(IX+1) 1 *Q2(II,NM+1)*SIGN(1.0,DA(I,J,M)) ENDIF ! Y-SPACE IF(MOD(IY,2).EQ.1) THEN XNJ(IIY)=XNJ(IIY)+CST(IY)*WY(IY+1) 1 *Q2(II,NM+1) ELSE XNJ(IIY)=XNJ(IIY)+CST(IY)*WY(IY+1) 1 *Q2(II,NM+1)*SIGN(1.0,DB(I,M)) ENDIF ENDDO ENDDO ENDDO IF(IELEM.EQ.1.AND.LFIXUP.AND.(XNI(1,J).LE.RLOG)) XNI(1,J)=0.0 IF(IELEM.EQ.1.AND.LFIXUP.AND.(XNJ(1).LE.RLOG)) XNJ(1)=0.0 WE=WE0 WX=WX0 WY=WY0 ! SAVE ENERGY GROUP LOWER BOUNDARY OUTGOING FLUX FLUX0(:NME,M,I,J)=REAL(FEP(:NME))/DELTAE(IG) ! SAVE LAST GROUP LOWER BOUNDARY FLUX IF(ISLG(IG).EQ.1) THEN FLUXC(I,J)=FLUXC(I,J)+REAL(FLUX0(1,M,I,J))*DN(1,M) ENDIF ! SAVE LEGENDRE MOMENT OF THE FLUX DO P=1,NSCT DO IEL=1,NM FLUX(IEL,P,I,J)=FLUX(IEL,P,I,J)+Q2(IEL,NM+1)*DN(P,M) ENDDO ENDDO 140 CONTINUE ! END OF Y-LOOP ! SAVE Y-BOUNDARY CONDITIONS DO IEL=1,NMY IOF=(M-1)*NMY*LX+(I-1)*NMY+IEL FUNKNO(LFLX+LXNI+IOF,IG)=REAL(XNJ(IEL)) ENDDO 155 CONTINUE ! END OF X-LOOP ! SAVE BOUNDARY CONDITIONS DO J=1,LY DO IEL=1,NMX IOF=(M-1)*NMX*LY+(J-1)*NMX+IEL FUNKNO(LFLX+IOF,IG)=REAL(XNI(IEL,J)) ENDDO ENDDO ! SAVE FLUX INFORMATION FLUX_G(:,:,:,:,IG)=FLUX_G(:,:,:,:,IG)+FLUX(:,:,:,:) FLUX0_G(:,:,:,:,IG)=FLUX0_G(:,:,:,:,IG)+FLUX0(:,:,:,:) 170 CONTINUE ! END OF DIRECTION LOOP 180 CONTINUE ! END OF ENERGY LOOP *$OMP END PARALLEL DO 190 CONTINUE ! END OF OCTANT LOOP ! SAVE FLUX INFORMATION DO 200 IG=1,NGEFF IF(.NOT.INCONV(IG)) GO TO 200 FUNKNO(:LFLX,IG)= 1 RESHAPE(REAL(FLUX_G(:NM,:NSCT,:LX,:LY,IG)), 2 (/ LFLX /) ) FUNKNO(LFLX+LXNI+LXNJ+1:LFLX+LXNI+LXNJ+LFEP,IG)= 1 RESHAPE(REAL(FLUX0_G(:NME,:NPQ,:LX,:LY,IG)), (/ LFEP /) ) 200 CONTINUE *---- * SCRATCH STORAGE DEALLOCATION *---- DEALLOCATE(XNI,FLUX0_G,FLUX_G,FLUX0,FLUX,INDANG) RETURN 400 FORMAT(16H SNFP12: thread=,I8,12H --->(group=,I4,7H angle=,I4,1H)) END