summaryrefslogtreecommitdiff
path: root/Dragon/src/SNFKC2.F
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/SNFKC2.F')
-rw-r--r--Dragon/src/SNFKC2.F562
1 files changed, 562 insertions, 0 deletions
diff --git a/Dragon/src/SNFKC2.F b/Dragon/src/SNFKC2.F
new file mode 100644
index 0000000..c95de6b
--- /dev/null
+++ b/Dragon/src/SNFKC2.F
@@ -0,0 +1,562 @@
+*DECK SNFKC2
+ SUBROUTINE SNFKC2(NKBA,NUN,NGEFF,IMPX,INCONV,NGIND,LX,LY,IELEM,
+ 1 NM,NMX,NMY,NMAT,NPQ,NSCT,MAT,VOL,TOTAL,NCODE,ZCODE,QEXT,LFIXUP,
+ 2 DU,DE,W,MRM,MRMY,DB,DA,MN,DN,WX,WY,CST,ISADPT,ISBS,NBS,ISBSM,BS,
+ 3 MAXL,FUNKNO)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Perform one inner iteration for solving SN equations in 2D Cartesian
+* geometry for the HODD method. KBA-like multithreading. Albedo
+* boundary conditions. Boltzmann (BTE) discretization.
+*
+*Copyright:
+* Copyright (C) 2025 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
+* NKBA number of macrocells along each axis
+* 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.
+* NM number of moments in space and energy for flux 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.
+* NCODE boundary condition indices.
+* ZCODE albedos.
+* 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.
+* CST constants for the polynomial approximations.
+* ISADPT flag to enable/disable adaptive flux calculations.
+*
+*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 NKBA,NUN,NGEFF,IMPX,NGIND(NGEFF),LX,LY,IELEM,NM,NMX,NMY,
+ 1 NMAT,NPQ,NSCT,MAT(LX,LY),NCODE(4),MRM(NPQ),MRMY(NPQ),ISBS,NBS,
+ 2 ISBSM(4*ISBS,NPQ*ISBS,NGEFF*ISBS),MAXL
+ LOGICAL INCONV(NGEFF)
+ REAL VOL(LX,LY),TOTAL(0:NMAT,NGEFF),ZCODE(4),QEXT(NUN,NGEFF),
+ 1 DU(NPQ),DE(NPQ),W(NPQ),DB(LX,NPQ),DA(LX,LY,NPQ),
+ 2 FUNKNO(NUN,NGEFF),BS(MAXL*ISBS,NBS*ISBS),WX(IELEM+1),
+ 3 WY(IELEM+1),CST(IELEM),MN(NPQ,NSCT),DN(NSCT,NPQ)
+ LOGICAL LFIXUP,ISADPT(2)
+*----
+* LOCAL VARIABLES
+*----
+ INTEGER NPQD(4),IIND(4),P
+ REAL WX0(IELEM+1),WY0(IELEM+1)
+ DOUBLE PRECISION Q(NM),Q2(NM,NM+1),V
+ PARAMETER(IUNOUT=6,RLOG=1.0E-8,PI=3.141592654)
+ LOGICAL ISFIX(2)
+*----
+* ALLOCATABLE ARRAYS
+*----
+ INTEGER, ALLOCATABLE, DIMENSION(:,:) :: INDANG
+ DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: FLUX
+ DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:,:,:) :: FLUX_G
+ DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:,:) :: XNI,XNJ
+*----
+* SCRATCH STORAGE ALLOCATION
+*----
+ ALLOCATE(INDANG(NPQ,4))
+ ALLOCATE(FLUX(NM,NSCT))
+ ALLOCATE(FLUX_G(NM,NSCT,LX,LY,NGEFF))
+ ALLOCATE(XNI(IELEM,LY,NPQ,NGEFF),XNJ(IELEM,LX,NPQ,NGEFF))
+*----
+* LENGTH OF FUNKNO COMPONENTS (IN ORDER)
+*----
+ LFLX=NM*LX*LY*NSCT
+ LXNI=NMX*LY*NPQ
+ LXNJ=NMY*LX*NPQ
+*----
+* NUMBER OF MACROCELLS (MACRO*)
+* NUMBER OF LZ LAYERS IN EACH MACROCELL (NCELL*)
+*----
+ MACROX=NKBA
+ MACROY=NKBA
+ NCELLX=1+(LX-1)/MACROX
+ NCELLY=1+(LY-1)/MACROY
+*----
+* SET OCTANT SWAPPING ORDER.
+*----
+ NPQD(:4)=0
+ INDANG(:NPQ,:4)=0
+ IIND(:)=0
+ DO M=1,NPQ
+ VU=DU(M)
+ VE=DE(M)
+ IF(W(M).EQ.0) CYCLE
+ 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
+ WX0=WX
+ WY0=WY
+
+ DO 240 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
+
+*----
+* KBA DIAGONAL LOOP
+*----
+ XNI(:IELEM,:LY,:NPQ,:NGEFF)=0.0D0
+ XNJ(:IELEM,:LX,:NPQ,:NGEFF)=0.0D0
+ DO 230 IDI=1,MACROX+MACROY-1
+
+*----
+* MAIN SWAPPING LOOPS FOR SN FLUX CALCULATION
+*----
+
+*$OMP PARALLEL DO
+*$OMP+ PRIVATE(ITID,FLUX,M,IG,Q,Q2,IOF,IER,II,JJ,IEL,I,J,P)
+*$OMP+ PRIVATE(IPQD,IMX,IMY,IBM,SIGMA,V,ISFIX,IX,JX,IY,JY,JCEL)
+*$OMP+ FIRSTPRIVATE(WX,WY,WX0,WY0) SHARED(FUNKNO,XNI,XNJ)
+*$OMP+ REDUCTION(+:FLUX_G) COLLAPSE(3)
+
+ ! LOOP FOR GROUPS TO EXECUTE IN PARALLEL
+ DO 220 IG=1,NGEFF
+
+ ! LOOP OVER ALL DIRECTIONS
+ DO 210 IPQD=1,NPQD(IND)
+
+ ! LOOP OVER MACROCELLS IN WAVEFRONT
+ DO 200 ICEL=MAX(1,IDI-MACROY+1),MIN(MACROX,IDI)
+
+ IF(.NOT.INCONV(IG)) GO TO 200
+ M=INDANG(IPQD,IND)
+ IF(W(M).EQ.0.0) GO TO 200
+
+ ! 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
+
+*----
+* LOOP OVER X- AND Y-DIRECTED AXES.
+*----
+ JCEL=IDI-ICEL+1
+
+ ! X-BOUNDARIES CONDITIONS
+ IF(ICEL.EQ.1) THEN
+ DO IMY=1,MIN(NCELLY,LY-(JCEL-1)*NCELLY)
+ J=(JCEL-1)*NCELLY+IMY
+ IF((IND.EQ.3).OR.(IND.EQ.4)) J=LY+1-J
+ 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,IPQD,IG)=FUNKNO(LFLX+IOF,IG)*ZCODE(1)
+ ELSE
+ XNI(IEL,J,IPQD,IG)=FUNKNO(LFLX+IOF,IG)*ZCODE(2)
+ ENDIF
+ ENDDO
+
+ ! X-BOUNDARIES FIXED SOURCES
+ IF(ISBS.EQ.1) THEN
+ IF((IND.EQ.2.OR.IND.EQ.3).AND.ISBSM(2,M,IG).NE.0) THEN
+ XNI(1,J,IPQD,IG)=XNI(1,J,IPQD,IG)+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,IPQD,IG)=XNI(1,J,IPQD,IG)+BS(J,ISBSM(1,M,IG))
+ ENDIF
+ ENDIF
+ ENDDO
+ ENDIF
+
+ ! Y-BOUNDARIES CONDITIONS
+ IF(JCEL.EQ.1) THEN
+ DO IMX=1,MIN(NCELLX,LX-(ICEL-1)*NCELLX)
+ I=(ICEL-1)*NCELLX+IMX
+ IF((IND.EQ.2).OR.(IND.EQ.3)) I=LX+1-I
+ DO IEL=1,NMY
+ IOF=(M-1)*NMY*LX+(I-1)*NMY+IEL
+ IF((IND.EQ.1).OR.(IND.EQ.2)) THEN
+ XNJ(IEL,I,IPQD,IG)=FUNKNO(LFLX+LXNI+IOF,IG)*ZCODE(3)
+ ELSE
+ XNJ(IEL,I,IPQD,IG)=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,I,IPQD,IG)=XNJ(1,I,IPQD,IG)+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,I,IPQD,IG)=XNJ(1,I,IPQD,IG)+BS(I,ISBSM(3,M,IG))
+ ENDIF
+ ENDIF
+ ENDDO
+ ENDIF
+
+ ! X-AXIS LOOP
+ DO 190 IMX=1,MIN(NCELLX,LX-(ICEL-1)*NCELLX)
+ I=(ICEL-1)*NCELLX+IMX
+ IF((IND.EQ.2).OR.(IND.EQ.3)) I=LX+1-I
+
+ ! Y-AXIS LOOP
+ DO 180 IMY=1,MIN(NCELLY,LY-(JCEL-1)*NCELLY)
+ J=(JCEL-1)*NCELLY+IMY
+ IF((IND.EQ.3).OR.(IND.EQ.4)) J=LY+1-J
+
+ ! INITIALIZE FLUXES
+ FLUX(:NM,:NSCT)=0.0D0
+
+ ! IF(ITID.EQ.1)THEN
+ ! WRITE(*,*) 'I,J', JCEL, NCELLY, IMY
+ ! WRITE(*,*) IND, IPQD, I, J
+ ! WRITE(*,*) XNI(:,J,IPQD,IG)
+ ! WRITE(*,*) XNJ(:,I,IPQD,IG)
+ ! ENDIF
+
+ ! DATA
+ IBM=MAT(I,J)
+ IF(IBM.EQ.0) GO TO 200
+ SIGMA=TOTAL(IBM,IG)
+ V=VOL(I,J)
+
+ ! 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
+
+ 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
+ II=IELEM*(IY-1)+IX
+ JJ=IELEM*(JY-1)+JX
+
+ ! DIAGONAL TERMS
+ IF(II.EQ.JJ) THEN
+ Q2(II,JJ)=SIGMA*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
+ ! X-SPACE TERMS
+ IF(IY.EQ.JY) 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
+ ! Y-SPACE TERMS
+ ELSEIF(IX.EQ.JX) 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
+ ! X-SPACE TERMS
+ IF(IY.EQ.JY) THEN
+ IF(MOD(IX+JX,2).EQ.1) THEN
+ Q2(II,JJ)=CST(IX)*CST(JX)*(WX(JX+1)-2)*DA(I,J,M)
+ ELSE
+ Q2(II,JJ)=CST(IX)*CST(JX)*WX(JX+1)*ABS(DA(I,J,M))
+ ENDIF
+ ! Y-SPACE TERMS
+ ELSEIF(IX.EQ.JX) THEN
+ IF(MOD(IY+JY,2).EQ.1) THEN
+ Q2(II,JJ)=CST(IY)*CST(JY)*(WY(JY+1)-2)*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
+
+ ! FLUX SOURCE VECTOR
+ DO IY=1,IELEM
+ DO IX=1,IELEM
+ II=IELEM*(IY-1)+IX
+ Q2(II,NM+1)=Q(II)*V
+ ! 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(IY,J,IPQD,IG)*ABS(DA(I,J,M))
+ ELSE
+ Q2(II,NM+1)=Q2(II,NM+1)-CST(IX)*(1+WX(1))
+ 1 *XNI(IY,J,IPQD,IG)*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(IX,I,IPQD,IG)*ABS(DB(I,M))
+ ELSE
+ Q2(II,NM+1)=Q2(II,NM+1)-CST(IY)*(1+WY(1))
+ 1 *XNJ(IX,I,IPQD,IG)*DB(I,M)
+ ENDIF
+ ENDDO
+ ENDDO
+
+ CALL ALSBD(NM,1,Q2,IER,NM)
+ IF(IER.NE.0) CALL XABORT('SNFKC2: SINGULAR MATRIX.')
+
+ ! IF(ITID.EQ.1)THEN
+ ! WRITE(*,*) 'I,J'
+ ! WRITE(*,*) Q2(:,NM+1)
+ ! WRITE(*,*) ' '
+ ! WRITE(*,*) ' '
+ ! WRITE(*,*) ' '
+ ! WRITE(*,*) ' '
+ ! ENDIF
+
+ ! ADAPTIVE CORRECTION OF WEIGHTING PARAMETERS
+ IF(ANY(ISADPT)) THEN
+ IF(ISADPT(1)) THEN
+ CALL SNADPT(IELEM,NM,IELEM,Q2(1:IELEM:1,NM+1),
+ 1 XNI(:NMX,J,IPQD,IG),1.0,WX,ISFIX(1))
+ ELSE
+ ISFIX(1)=.TRUE.
+ ENDIF
+ IF(ISADPT(2)) THEN
+ CALL SNADPT(IELEM,NM,IELEM,Q2(1:NM:IELEM,NM+1),
+ 1 XNJ(:NMY,I,IPQD,IG),1.0,WY,ISFIX(2))
+ ELSE
+ ISFIX(2)=.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,IPQD,IG)=WX(1)*XNI(:NMX,J,IPQD,IG)
+ XNJ(:NMY,I,IPQD,IG)=WY(1)*XNJ(:NMY,I,IPQD,IG)
+ DO IY=1,IELEM
+ DO IX=1,IELEM
+ II=IELEM*(IY-1)+IX
+ ! X-SPACE
+ IF(MOD(IX,2).EQ.1) THEN
+ XNI(IY,J,IPQD,IG)=XNI(IY,J,IPQD,IG)+CST(IX)*WX(IX+1)
+ 1 *Q2(II,NM+1)
+ ELSE
+ XNI(IY,J,IPQD,IG)=XNI(IY,J,IPQD,IG)+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(IX,I,IPQD,IG)=XNJ(IX,I,IPQD,IG)+CST(IY)*WY(IY+1)
+ 1 *Q2(II,NM+1)
+ ELSE
+ XNJ(IX,I,IPQD,IG)=XNJ(IX,I,IPQD,IG)+CST(IY)*WY(IY+1)
+ 1 *Q2(II,NM+1)*SIGN(1.0,DB(I,M))
+ ENDIF
+ ENDDO
+ ENDDO
+ IF(IELEM.EQ.1.AND.LFIXUP.AND.(XNI(1,J,IPQD,IG).LE.RLOG))
+ 1 XNI(1,J,IPQD,IG)=0.0
+ IF(IELEM.EQ.1.AND.LFIXUP.AND.(XNJ(1,I,IPQD,IG).LE.RLOG))
+ 1 XNJ(1,I,IPQD,IG)=0.0
+ WX=WX0
+ WY=WY0
+
+ ! SAVE LEGENDRE MOMENT OF THE FLUX
+ DO P=1,NSCT
+ DO IEL=1,NM
+ FLUX(IEL,P)=FLUX(IEL,P)+Q2(IEL,NM+1)*DN(P,M)
+ ENDDO
+ ENDDO
+
+ ! SAVE X-BOUNDARY CONDITIONS
+ IF((ICEL-1)*NCELLX+IMX.EQ.LX) THEN
+ DO IEL=1,NMX
+ IOF=(M-1)*NMX*LY+(J-1)*NMX+IEL
+ FUNKNO(LFLX+IOF,IG)=REAL(XNI(IEL,J,IPQD,IG))
+ ENDDO
+ ENDIF
+
+ ! SAVE Y-BOUNDARY CONDITIONS
+ IF((JCEL-1)*NCELLY+IMY.EQ.LY) THEN
+ DO IEL=1,NMY
+ IOF=(M-1)*NMY*LX+(I-1)*NMY+IEL
+ FUNKNO(LFLX+LXNI+IOF,IG)=REAL(XNJ(IEL,I,IPQD,IG))
+ ENDDO
+ ENDIF
+
+ ! SAVE FLUX INFORMATION
+ FLUX_G(:,:,I,J,IG)=FLUX_G(:,:,I,J,IG)+FLUX(:,:)
+
+ 180 CONTINUE ! END OF Y-LOOP
+ 190 CONTINUE ! END OF X-LOOP
+
+ 200 CONTINUE ! END OF MACROCELL LOOP
+ 210 CONTINUE ! END OF DIRECTION LOOP
+ 220 CONTINUE ! END OF ENERGY LOOP
+*$OMP END PARALLEL DO
+ 230 CONTINUE ! END OF WAVEFRONT LOOP
+ 240 CONTINUE ! END OF OCTANT LOOP
+
+ ! SAVE FLUX INFORMATION
+ DO 250 IG=1,NGEFF
+ IF(.NOT.INCONV(IG)) GO TO 250
+ FUNKNO(:LFLX,IG)=
+ 1 RESHAPE(REAL(FLUX_G(:NM,:NSCT,:LX,:LY,IG)),(/ LFLX /) )
+ 250 CONTINUE
+
+ ! CALL XABORT('SNFKC2: testing')
+
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ DEALLOCATE(XNI,XNJ,FLUX_G,FLUX,INDANG)
+ RETURN
+ 400 FORMAT(16H SNFKC2: thread=,I8,12H --->(group=,I4,7H angle=,I4,1H))
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 13:22:23 +0000