Initial commit from Polytechnique Montreal

1 files changed, 793 insertions, 0 deletions

diff --git a/Dragon/src/SNFBH3.F b/Dragon/src/SNFBH3.F
new file mode 100644
index 0000000..cc7b8d7
--- /dev/null
+++ b/Dragon/src/SNFBH3.F
@@ -0,0 +1,793 @@
+*DECK SNFBH3
+      SUBROUTINE SNFBH3(NUN,NGEFF,IMPX,INCONV,NGIND,NHEX,LZ,ISPLH,SIDE,
+     1 IELEM,NM,NMX,NMY,NMZ,NMAT,NPQ,NSCT,MAT,VOL,TOTAL,NCODE,ZCODE,
+     2 QEXT,LFIXUP,DU,DE,DZ,W,MRMZ,DC,DB,DA,MN,DN,WX,WY,WZ,CST,ISADPT,
+     3 LOZSWP,COORDMAP,FUNKNO)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Perform one inner iteration for solving SN equations in 3D Cartesian
+* geometry for the HODD method. Energy-angle multithreading. Albedo
+* boundary conditions on top/bottom, Void on sides. 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
+* 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.
+* NHEX     number of hexagons in X-Y plane.
+* ISPLH    splitting option for hexagons.
+* SIDE     side of an hexagon.
+* LZ       number of meshes along Z 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 for flux components
+* NMX      number of moments for X axis boundaries components
+* NMY      number of moments for Y axis boundaries components
+* NMZ      number of moments for Z axis boundaries components
+* NMAT     number of material mixtures.
+* NPQ      number of SN directions in height octants.
+* 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$).
+* DZ       third direction cosines ($\\xi$).
+* W        weights.
+* MRMZ     quadrature index.
+* DC       diamond-scheme parameter.
+* DB       diamond-scheme parameter.
+* DA       diamond-scheme parameter.
+* MN       moment-to-discrete matrix.
+* DN       discrete-to-moment matrix.
+* WX       spatial X axis closure relation weighting factors.
+* WY       spatial Y axis closure relation weighting factors.
+* WZ       spatial Z axis closure relation weighting factors.
+* CST      constants for the polynomial approximations.
+* ISADPT   flag to enable/disable adaptive flux calculations.
+* LOZSWP   lozenge sweep order depending on direction.
+* COORDMAP coordinate map - mapping the hexagons from the indices 
+*          within the DRAGON geometry to a Cartesian axial coordinate
+*          array (see redblobgames.com website).
+*
+*Parameters: input/output
+* FUNKNO  Legendre components of the flux and boundary fluxes.
+*
+*Comments:
+* 1. The direction of the axes I, J and D for the surface boundary 
+*    fluxes are shown in the diagram below. This means that 
+*    i) lozenge A has I- and D-boundaries (instead of I and J)
+*    i) lozenge B has I- and J-boundaries
+*    i) lozenge C has D- and J-boundaries (instead of I and J)
+*
+*                                  ^
+*                         j-axis   |
+*                                  |          ^
+*                              _________     /    d-axis
+*                             /       / \   /
+*                            /   B   /   \
+*                           /       /     \
+*                          (-------(   A   )
+*                           \       \     /
+*                            \  C    \   / 
+*                             \_______\_/   \
+*                                            \   i-axis
+*                                             ^
+*
+*-----------------------------------------------------------------------
+*
+#if defined(_OPENMP)
+      USE omp_lib
+#endif
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER NUN,NGEFF,IMPX,NGIND(NGEFF),NHEX,LZ,ISPLH,IELEM,NM,NMX,
+     1 NMY,NMZ,NMAT,NPQ,NSCT,MAT(ISPLH,ISPLH,3,NHEX,LZ),NCODE(6),
+     2 MRMZ(NPQ),LOZSWP(3,6),COORDMAP(3,NHEX)
+      LOGICAL INCONV(NGEFF)
+      REAL SIDE,VOL(ISPLH,ISPLH,3,NHEX,LZ),TOTAL(0:NMAT,NGEFF),ZCODE(6),
+     1 QEXT(NUN,NGEFF),DU(NPQ),DE(NPQ),DZ(NPQ),W(NPQ),
+     2 DC(ISPLH*ISPLH*3*NHEX,1,NPQ),DB(ISPLH*ISPLH*3*NHEX,LZ,NPQ),
+     3 DA(1,LZ,NPQ),FUNKNO(NUN,NGEFF),WX(IELEM+1),WY(IELEM+1),
+     4 WZ(IELEM+1),CST(IELEM),MN(NPQ,NSCT),DN(NSCT,NPQ)
+      LOGICAL LFIXUP,ISADPT(3)
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER :: NPQD(12),IIND(12),P,DCOORD
+      REAL :: JAC(2,2,3),MUH,ETAH,XIH,AAA,BBB,CCC,DDD,MUHTEMP,ETAHTEMP,
+     1 WX0(IELEM+1),WY0(IELEM+1),WZ0(IELEM+1)
+      DOUBLE PRECISION :: V,Q(NM),Q2(NM,NM+1),THETA,XNI(NMX),
+     > XNJ(NMY),XNK(NMZ)
+      PARAMETER(IUNOUT=6,RLOG=1.0E-8,PI=3.141592654)
+      LOGICAL ISFIX(3),LHEX(NHEX)
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: INDANG
+      INTEGER, ALLOCATABLE, DIMENSION(:,:,:,:,:,:) :: TMPMAT
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:,:,:) :: FLUX
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:,:,:,:) :: FLUX_G
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: TMPXNI,
+     > TMPXNJ, TMPXND
+      DOUBLE PRECISION,ALLOCATABLE, DIMENSION(:,:,:,:,:) :: TMPXNK
+*----
+*  MAP MATERIAL VALUES TO CARTESIAN AXIAL COORDINATE MAP
+*----
+      NRINGS=INT((SQRT(  REAL((4*NHEX-1)/3)  )+1.)/2.)
+      NCOL=2*NRINGS -1
+      ALLOCATE(TMPMAT(ISPLH,ISPLH,3,NCOL,NCOL,LZ))
+      TMPMAT(:,:,:,:,:,:) = -1
+      DO K=1,LZ
+        DO IHEX_XY=1,NHEX
+          TMPMAT(:,:,:,COORDMAP(1,IHEX_XY),COORDMAP(2,IHEX_XY),K) = 
+     >         MAT(:,:,:,IHEX_XY,K)
+        ENDDO
+      ENDDO
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(INDANG(NPQ,12))
+      ALLOCATE(FLUX(NM,NSCT,3*ISPLH**2,NHEX,LZ))
+      ALLOCATE(FLUX_G(NM,NSCT,3*ISPLH**2,NHEX,LZ,NGEFF))
+      ALLOCATE(TMPXNI(NMX,ISPLH,NCOL))
+      ALLOCATE(TMPXNJ(NMY,ISPLH,NCOL))
+      ALLOCATE(TMPXND(NMX,ISPLH,NCOL))
+      ALLOCATE(TMPXNK(NMZ,ISPLH,ISPLH,3,NHEX))
+*----
+*  CONSTRUCT JACOBIAN MATRIX FOR EACH LOZENGE
+*----
+      JAC = RESHAPE((/ 1., -SQRT(3.), 1., SQRT(3.), 2., 0., 1.,
+     >    SQRT(3.), 2., 0., -1., SQRT(3.) /), SHAPE(JAC))
+      JAC = (SIDE/2.)*JAC
+*----
+*  LENGTH OF FUNKNO COMPONENTS (IN ORDER)
+*----
+      LFLX=3*NM*(ISPLH**2)*NHEX*LZ*NSCT
+      L5=3*NMZ*(ISPLH**2)*NHEX
+*----
+*  SET DODECANT SWAPPING ORDER.
+*----
+      NPQD(:12)=0
+      INDANG(:NPQ,:12)=0
+      IIND(:12)=0
+      DO M=1,NPQ
+        VU=DU(M)
+        VE=DE(M)
+        VZ=DZ(M)
+        IF(W(M).EQ.0) CYCLE
+        THETA=0.0D0
+        IF(VE.GT.0.0)THEN
+          IF(VU.EQ.0.0)THEN
+            THETA = PI/2
+          ELSEIF(VU.GT.0.0)THEN
+            THETA = ATAN(ABS(VE/VU))
+          ELSEIF(VU.LT.0.0)THEN
+            THETA = PI - ATAN(ABS(VE/VU))
+          ENDIF
+        ELSEIF(VE.LT.0.0)THEN
+          IF(VU.EQ.0.0)THEN
+            THETA = 3*PI/2
+          ELSEIF(VU.LT.0.0)THEN
+            THETA = PI + ATAN(ABS(VE/VU))
+          ELSEIF(VU.GT.0.0)THEN
+            THETA = 2.*PI - ATAN(ABS(VE/VU))
+          ENDIF
+        ENDIF
+        ! UNFOLD DODECANTS
+        IND=0
+        IF(VZ.GE.0.0)THEN
+          IF((THETA.GT.0.0).AND.(THETA.LT.(PI/3.)))THEN
+            IND=1
+          ELSEIF((THETA.GT.(PI/3.)).AND.(THETA.LT.(2.*PI/3.)))THEN
+            IND=2
+          ELSEIF((THETA.GT.(2.*PI/3.)).AND.(THETA.LT.(PI)))THEN
+            IND=3
+          ELSEIF((THETA.GT.(PI)).AND.(THETA.LT.(4.*PI/3.)))THEN
+            IND=4
+          ELSEIF((THETA.GT.(4.*PI/3.)).AND.(THETA.LT.(5.*PI/3.)))THEN
+            IND=5
+          ELSEIF((THETA.GT.(5.*PI/3.)).AND.(THETA.LT.(2.*PI)))THEN
+            IND=6
+          ENDIF
+        ELSEIF(VZ.LT.0.0)THEN
+          IF((THETA.GT.0.0).AND.(THETA.LT.(PI/3.)))THEN
+            IND=7
+          ELSEIF((THETA.GT.(PI/3.)).AND.(THETA.LT.(2.*PI/3.)))THEN
+            IND=8
+          ELSEIF((THETA.GT.(2.*PI/3.)).AND.(THETA.LT.(PI)))THEN
+            IND=9
+          ELSEIF((THETA.GT.(PI)).AND.(THETA.LT.(4.*PI/3.)))THEN
+            IND=10
+          ELSEIF((THETA.GT.(4.*PI/3.)).AND.(THETA.LT.(5.*PI/3.)))THEN
+            IND=11
+          ELSEIF((THETA.GT.(5.*PI/3.)).AND.(THETA.LT.(2.*PI)))THEN
+            IND=12
+          ENDIF
+        ENDIF
+        ! Assume IIND(I)=I in hexagonal geometry
+        IIND(IND)=IND
+        NPQD(IND)=NPQD(IND)+1
+        INDANG(NPQD(IND),IND)=M
+      ENDDO
+*----
+*  MAIN LOOP OVER DODECANTS.
+*----
+
+      FLUX_G(:NM,:NSCT,:3*ISPLH**2,:NHEX,:LZ,:NGEFF)=0.0D0
+      WX0=WX
+      WY0=WY
+      WZ0=WZ
+
+      DO JND=1,12
+      IND=IIND(JND)
+      IND_XY=MOD(IND-1,6)+1
+      ! Needed because of S2 LS (8 dir. for 12 dodecants)
+      IF(IND.EQ.0) CYCLE
+*----
+*  PRELIMINARY LOOPS FOR SETTING BOUNDARY CONDITIONS.
+*----
+
+      IF((NCODE(5).NE.1).or.(NCODE(6).NE.1))THEN
+*$OMP  PARALLEL DO
+*$OMP+ PRIVATE(M,IG,VZ,M1,IOF,JOF,IPQD)
+*$OMP+ SHARED(FUNKNO) COLLAPSE(2)
+
+      DO IG=1,NGEFF
+      DO IPQD=1,NPQD(IND)
+      IF(.NOT.INCONV(IG)) CYCLE
+      M=INDANG(IPQD,IND)
+
+      VZ=DZ(M)
+      ! Z-BOUNDARY
+      IF(VZ.GT.0.0)THEN
+        M1=MRMZ(M)
+        IF(NCODE(5).NE.4)THEN
+          IOF=(M-1)*(L5)
+          JOF=(M1-1)*(L5)
+          FUNKNO(LFLX+IOF+1:LFLX+IOF+L5,IG)=
+     >      FUNKNO(LFLX+JOF+1:LFLX+JOF+L5,IG)
+        ENDIF
+      ELSEIF(VZ.LT.0.0)THEN
+        M1=MRMZ(M)
+        IF(NCODE(6).NE.4)THEN
+          IOF=(M-1)*(L5)
+          JOF=(M1-1)*(L5)
+          FUNKNO(LFLX+IOF+1:LFLX+IOF+L5,IG)=
+     >      FUNKNO(LFLX+JOF+1:LFLX+JOF+L5,IG)
+        ENDIF
+      ENDIF
+*
+      ENDDO
+      ENDDO
+
+*$OMP END PARALLEL DO
+      ENDIF
+
+      ! CALL XABORT('SNFBH3: testing 1 ')
+*----
+*  MAIN SWAPPING LOOPS FOR SN FLUX CALCULATION
+*  LOOP OVER ENERGY AND ANGLES
+*----
+
+*$OMP  PARALLEL DO 
+*$OMP+ PRIVATE(ITID,FLUX,M,IG,XNI,XNJ,XNK,Q,Q2,IOF,IER,II,JJ,I,J,K,K0)
+*$OMP+ PRIVATE(IPQD,IBM,SIGMA,V,ISFIX,IX,JX,IY,JY,IZ,JZ,AAA,BBB,CCC,DDD)
+*$OMP+ PRIVATE(IIX,IIY,IIZ,P,IEL)
+*$OMP+ PRIVATE(LHEX,IHEX_XY,IIHEX,DCOORD,ILOZLOOP,ILOZ,IL,I2,JL,J2)
+*$OMP+ PRIVATE(MUHTEMP,MUH,ETAHTEMP,ETAH,XIH,I3,I_FETCH,III,JJJ,IIM,JIM)
+*$OMP+ PRIVATE(TMPXNI,TMPXNJ,TMPXND,TMPXNK)
+*$OMP+ FIRSTPRIVATE(WX,WY,WZ,WX0,WY0,WZ0) SHARED(FUNKNO,IND)
+*$OMP+ REDUCTION(+:FLUX_G) COLLAPSE(2)      
+      ! LOOP FOR GROUPS TO EXECUTE IN PARALLEL
+      DO IG=1,NGEFF
+
+      ! LOOP OVER ALL DIRECTIONS
+      DO IPQD=1,NPQD(IND)
+      IF(.NOT.INCONV(IG)) CYCLE
+      M=INDANG(IPQD,IND)
+      IF(W(M).EQ.0.0) CYCLE
+
+      ! GET AND PRINT THREAD NUMBER 
+#if defined(_OPENMP)
+      ITID=omp_get_thread_num()
+#else
+      ITID=0
+#endif
+      IF(IMPX.GT.5) WRITE(IUNOUT,500) ITID,NGIND(IG),IPQD
+
+      ! INITIALIZE FLUX AND Z-BOUNDARY FLUXES
+      FLUX(:NM,:NSCT,:3*ISPLH**2,:NHEX,:LZ)=0.0D0
+      TMPXNK(:NMZ,:ISPLH,:ISPLH,:3,:NHEX)=0.0D0
+
+      ! PICK UP BOUNDARY ELEMENTS
+      IF((NCODE(5).NE.1).or.(NCODE(6).NE.1))THEN
+        IOF=(M-1)*(L5) + 1
+        TMPXNK(:NMZ,:ISPLH,:ISPLH,:3,:NHEX)=
+     >    RESHAPE(FUNKNO(LFLX+IOF:LFLX+IOF+L5,IG),
+     >    (/NMZ,ISPLH,ISPLH,3,NHEX/))
+      ENDIF
+      ! ACCOUNT FOR ALBEDO IN BOUNDARY ELEMENTS
+      IF(IND.LT.7) THEN
+        TMPXNK(:NMZ,:ISPLH,:ISPLH,:3,:NHEX)=
+     >    TMPXNK(:NMZ,:ISPLH,:ISPLH,:3,:NHEX)*ZCODE(5)
+      ELSE
+        TMPXNK(:NMZ,:ISPLH,:ISPLH,:3,:NHEX)=
+     >    TMPXNK(:NMZ,:ISPLH,:ISPLH,:3,:NHEX)*ZCODE(6)
+      ENDIF
+
+*----
+*  LOOP OVER Z-AXIS PLANES
+*----
+
+      DO K0=1,LZ
+      K=K0
+      IF(IND.GE.7) K=LZ+1-K0
+
+      ! INITIALIZE I,J,D-BOUNDARY FLUXES
+      TMPXNI(:NMX,:ISPLH,:NCOL)=0.0D0
+      TMPXNJ(:NMY,:ISPLH,:NCOL)=0.0D0
+      TMPXND(:NMX,:ISPLH,:NCOL)=0.0D0
+
+*----
+*  LOOP OVER CARTESIAN MAP OF HEXAGONAL DOMAIN
+*----
+
+      DO JJJ=1,NCOL
+      JIM=JJJ
+      ! Account for different sweep direction depending on angle 
+      IF((IND_XY.EQ.1).OR.(IND_XY.EQ.2).OR.(IND_XY.EQ.3)) JIM=NCOL+1-JIM
+        
+      DO III=1,NCOL
+      IIM=III
+      ! Account for different sweep direction depending on angle 
+      IF((IND_XY.EQ.2).OR.(IND_XY.EQ.3).OR.(IND_XY.EQ.4)) IIM=NCOL+1-IIM
+
+      ! For IND_XY 3 or 6, Cartesian axial coordinate map is swept 
+      ! vertically instead of horizontally. IM suffix is for 'IMmutable'
+      I=IIM
+      J=JIM
+      IF((IND_XY.EQ.3).OR.(IND_XY.EQ.6))THEN
+        I=JIM
+        J=IIM
+      ENDIF
+
+      ! If within corners of Cartesian axial coordinate map (where
+      ! there are no hexagons), skip loop
+      IF(TMPMAT(1,1,1,I,J,K).EQ.-1) CYCLE
+
+      ! Find in X-Y plane DRAGON geometry hexagonal index using I and J
+      LHEX=(COORDMAP(1,:).EQ.I .AND. COORDMAP(2,:).EQ.J)
+      IHEX_XY=0
+      DO IIHEX=1,NHEX
+        IF(LHEX(IIHEX)) THEN
+          IHEX_XY=IIHEX
+          EXIT
+        ENDIF
+      ENDDO
+      IF(IHEX_XY.EQ.0) CALL XABORT('SNFBH3: IHEX_XY FAILURE.')
+      ! Find D coordinate
+      DCOORD = ABS(COORDMAP(3,IHEX_XY))-NRINGS
+
+*----
+*  LOOP OVER LOZENGES
+*----
+
+      DO ILOZLOOP=1,3
+      ILOZ=LOZSWP(ILOZLOOP,IND_XY)
+
+      ! Get Jacobian elements values
+      AAA = JAC(1,1,ILOZ)
+      BBB = JAC(1,2,ILOZ)
+      CCC = JAC(2,1,ILOZ)
+      DDD = JAC(2,2,ILOZ)
+
+      ! CALL XABORT('SNFBH3: testing 19 ')
+*----
+*  LOOP OVER SUBMESH WITHIN EACH LOZENGE
+*----
+      DO IL=1,ISPLH
+      I2=IL
+      ! Account for different sweep direction depending on angle 
+      IF((ILOZ.EQ.1).OR.(ILOZ.EQ.2))THEN
+        IF((IND_XY.EQ.2).OR.(IND_XY.EQ.3).OR.(IND_XY.EQ.4))I2=ISPLH+1-I2
+      ELSEIF(ILOZ.EQ.3)THEN
+        IF((IND_XY.EQ.3).OR.(IND_XY.EQ.4).OR.(IND_XY.EQ.5))I2=ISPLH+1-I2
+      ENDIF
+
+      DO JL=1,ISPLH
+      J2=JL
+      ! Account for different sweep direction depending on angle 
+      IF((ILOZ.EQ.2).OR.(ILOZ.EQ.3))THEN
+        IF((IND_XY.EQ.4).OR.(IND_XY.EQ.5).OR.(IND_XY.EQ.6))J2=ISPLH+1-J2
+      ELSEIF(ILOZ.EQ.1)THEN
+        IF((IND_XY.EQ.3).OR.(IND_XY.EQ.4).OR.(IND_XY.EQ.5))J2=ISPLH+1-J2
+      ENDIF
+      ! READ IN XNI AND XNJ DEPENDING ON LOZENGE
+      I_FETCH=0
+      IF((ILOZ.EQ.1))THEN
+        ! Read boundary fluxes in reverse for lozenge A since affine
+        ! transformation of lozenges causes the D and I directions
+        ! of lozenges C and A respectively to be reversed
+        I_FETCH=ISPLH+1-I2
+        XNI(:) = TMPXNI(:,J2,J)
+        XNJ(:) = TMPXND(:,I_FETCH,DCOORD)
+      ELSEIF((ILOZ.EQ.2))THEN
+        XNI(:) = TMPXNI(:,J2,J)
+        XNJ(:) = TMPXNJ(:,I2,I)
+      ELSEIF((ILOZ.EQ.3))THEN
+        XNI(:) = TMPXND(:,J2,DCOORD)
+        XNJ(:) = TMPXNJ(:,I2,I)
+      ENDIF
+      XNK(:) = TMPXNK(:,I2,J2,ILOZ,IHEX_XY)
+
+      ! DATA
+      IBM=MAT(I2,J2,ILOZ,IHEX_XY,K)
+      ! Skip loop if virtual element
+      IF(IBM.EQ.0) CYCLE
+      SIGMA=TOTAL(IBM,IG)
+      V=VOL(I2,J2,ILOZ,IHEX_XY,K)
+
+      ! COMPUTE ADJUSTED DIRECTION COSINES
+      MUHTEMP  =  DA(1,K,M)
+      ETAHTEMP =  DB(1,K,M)
+      MUH = (MUHTEMP*DDD) - (ETAHTEMP*BBB)
+      ETAH = (-MUHTEMP*CCC) + (ETAHTEMP*AAA)
+      XIH  = DC(1,1,M)
+
+      ! IF(IND.EQ.12) CALL XABORT('SNFBH3: testing 60 ')
+    !   WRITE(*,*) (((((((K-1)*NHEX+(IHEX_XY-1))*3+(ILOZ-1))*ISPLH+
+    !  >  (J2-1))*ISPLH+(I2-1))*NSCT+(2-1))*NM)+1
+    !   WRITE(*,*) K, NHEX, IHEX_XY, ILOZ, ISPLH, J2, I2, NSCT, NM
+    !   CALL XABORT('SNFBH3: testing 2')
+      ! SOURCE DENSITY TERM
+      DO IEL=1,NM
+      Q(IEL)=0.0D0
+      DO P=1,NSCT
+      IOF=(((((((K-1)*NHEX+(IHEX_XY-1))*3+(ILOZ-1))*ISPLH+(J2-1))*
+     >  ISPLH+(I2-1))*NSCT+(P-1))*NM)+IEL
+      Q(IEL)=Q(IEL)+QEXT(IOF,IG)*MN(M,P)
+      ENDDO
+      ENDDO
+
+      ! CALL XABORT('SNFBH3: testing 17 ')
+      ISFIX=.FALSE.
+      DO WHILE (.NOT.ALL(ISFIX)) ! LOOP FOR ADAPTIVE CALCULATION
+
+      ! FLUX MOMENT COEFFICIENTS MATRIX
+      Q2(:NM,:NM+1)=0.0D0
+      DO IZ=1,IELEM
+      DO JZ=1,IELEM
+      DO IY=1,IELEM
+      DO JY=1,IELEM
+      DO IX=1,IELEM
+      DO JX=1,IELEM
+
+        II=IELEM**2*(IZ-1)+IELEM*(IY-1)+IX
+        JJ=IELEM**2*(JZ-1)+IELEM*(JY-1)+JX
+
+        ! IF(IPQD.EQ.3) CALL XABORT('SNFBH3: testing 69 ')
+      ! CALL XABORT('SNFBH3: testing 17 ')
+        ! DIAGONAL TERMS
+        IF(II.EQ.JJ) THEN
+          Q2(II,JJ)=SIGMA*V
+     1              +CST(IX)**2*WX(JX+1)*ABS(MUH)
+     2              +CST(IY)**2*WY(JY+1)*ABS(ETAH)
+     3              +CST(IZ)**2*WZ(JZ+1)*ABS(XIH)
+
+        ! IF(IND.EQ.12) CALL XABORT('SNFBH3: testing 70 ')
+      ! CALL XABORT('SNFBH3: testing 191 ')
+        ! UPPER DIAGONAL TERMS
+        ELSEIF(II.LT.JJ) THEN
+          ! CALL XABORT('SNFBH3: testing 1919 ')
+          IF(IZ.EQ.JZ) THEN
+            IF(IY.EQ.JY) THEN
+              ! X-SPACE TERMS
+              IF(MOD(IX+JX,2).EQ.1) THEN
+                Q2(II,JJ)=CST(IX)*CST(JX)*WX(JX+1)*MUH
+              ELSE
+                Q2(II,JJ)=CST(IX)*CST(JX)*WX(JX+1)*ABS(MUH)
+              ENDIF
+            ELSEIF(IX.EQ.JX) THEN
+              ! Y-SPACE TERMS
+              IF(MOD(IY+JY,2).EQ.1) THEN
+                Q2(II,JJ)=CST(IY)*CST(JY)*WY(JY+1)*ETAH
+              ELSE
+                Q2(II,JJ)=CST(IY)*CST(JY)*WY(JY+1)*ABS(ETAH)
+              ENDIF
+            ENDIF
+          ELSEIF(IY.EQ.JY.AND.IX.EQ.JX) THEN
+            ! Z-SPACE TERMS
+            IF(MOD(IZ+JZ,2).EQ.1) THEN
+              Q2(II,JJ)=CST(IZ)*CST(JZ)*WZ(JZ+1)*XIH
+            ELSE
+              Q2(II,JJ)=CST(IZ)*CST(JZ)*WZ(JZ+1)*ABS(XIH)
+            ENDIF
+          ENDIF
+        ! IF(IND.EQ.12) CALL XABORT('SNFBH3: testing 75 ')
+      ! CALL XABORT('SNFBH3: testing 19 ')
+
+         ! UNDER DIAGONAL TERMS
+        ELSE
+          ! CALL XABORT('SNFBH3: testing 1920 ')
+          IF(IZ.EQ.JZ) THEN
+            IF(IY.EQ.JY) THEN
+              ! X-SPACE TERMS
+              IF(MOD(IX+JX,2).EQ.1) THEN
+                Q2(II,JJ)=CST(IX)*CST(JX)*(WX(JX+1)-2)*MUH
+              ELSE
+                Q2(II,JJ)=CST(IX)*CST(JX)*WX(JX+1)*ABS(MUH)
+              ENDIF
+            ELSEIF(IX.EQ.JX) THEN
+              ! Y-SPACE TERMS
+              IF(MOD(IY+JY,2).EQ.1) THEN
+                Q2(II,JJ)=CST(IY)*CST(JY)*(WY(JY+1)-2)*ETAH
+              ELSE
+                Q2(II,JJ)=CST(IY)*CST(JY)*WY(JY+1)*ABS(ETAH)
+              ENDIF
+            ENDIF
+          ELSEIF(IY.EQ.JY.AND.IX.EQ.JX) THEN
+            ! Z-SPACE TERMS
+            IF(MOD(IZ+JZ,2).EQ.1) THEN
+              Q2(II,JJ)=CST(IZ)*CST(JZ)*(WZ(JZ+1)-2)*XIH
+            ELSE
+              Q2(II,JJ)=CST(IZ)*CST(JZ)*WZ(JZ+1)*ABS(XIH)
+            ENDIF
+          ENDIF
+        ENDIF
+      ENDDO
+      ENDDO
+      ENDDO
+      ENDDO
+      ENDDO
+      ENDDO
+
+      ! FLUX SOURCE VECTOR
+      DO IZ=1,IELEM
+      DO IY=1,IELEM
+      DO IX=1,IELEM
+        II=IELEM**2*(IZ-1)+IELEM*(IY-1)+IX
+        IIX=IELEM*(IZ-1)+IY
+        IIY=IELEM*(IZ-1)+IX
+        IIZ=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(IIX)*ABS(MUH)
+        ELSE
+          Q2(II,NM+1)=Q2(II,NM+1)-CST(IX)*(1+WX(1))
+     1                *XNI(IIX)*MUH
+        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(ETAH)
+        ELSE
+          Q2(II,NM+1)=Q2(II,NM+1)-CST(IY)*(1+WY(1))
+     1                *XNJ(IIY)*ETAH
+        ENDIF
+        ! Z-SPACE TERMS
+        IF(MOD(IZ,2).EQ.1) THEN
+          Q2(II,NM+1)=Q2(II,NM+1)+CST(IZ)*(1-WZ(1))
+     1                *XNK(IIZ)*ABS(XIH)
+        ELSE
+          Q2(II,NM+1)=Q2(II,NM+1)-CST(IZ)*(1+WZ(1))
+     1                *XNK(IIZ)*XIH
+        ENDIF
+      ENDDO
+      ENDDO
+      ENDDO
+
+      CALL ALSBD(NM,1,Q2,IER,NM)
+      IF(IER.NE.0) CALL XABORT('SNFBH3: SINGULAR MATRIX.')
+
+      ! ADAPTIVE CORRECTION OF WEIGHTING PARAMETERS
+      IF(ANY(ISADPT)) THEN
+        IF(ISADPT(1)) THEN
+          CALL SNADPT(IELEM,NM,IELEM**2,Q2(1:IELEM:1,NM+1),
+     1    XNI(:NMX),1.0,WX,ISFIX(1))
+        ELSE
+          ISFIX(1)=.TRUE.
+        ENDIF
+        IF(ISADPT(2)) THEN
+          CALL SNADPT(IELEM,NM,IELEM**2,Q2(1:IELEM**2:IELEM,NM+1),
+     1    XNJ(:NMY),1.0,WY,ISFIX(2))
+        ELSE
+          ISFIX(2)=.TRUE.
+        ENDIF
+        IF(ISADPT(3)) THEN
+          CALL SNADPT(IELEM,NM,IELEM**2,Q2(1:NM:IELEM**2,NM+1),
+     1    XNK(:NMZ),1.0,WZ,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
+      ! Read XNI/XNI/XNK into TMPXNI/J/D/K
+      IF((ILOZ.EQ.1).OR.(ILOZ.EQ.2))THEN
+        TMPXNI(:NMX,J2,J)=WX(1)*XNI(:NMX)
+      ELSE
+        TMPXND(:NMX,J2,DCOORD)=WX(1)*XNI(:NMX)
+      ENDIF
+      IF((ILOZ.EQ.2).OR.(ILOZ.EQ.3))THEN
+        TMPXNJ(:NMY,I2,I)=WY(1)*XNJ(:NMY)
+      ELSE
+        I3=I_FETCH
+        TMPXND(:NMY,I3,DCOORD)=WY(1)*XNJ(:NMY)
+      ENDIF
+      TMPXNK(:NMZ,I2,J2,ILOZ,IHEX_XY)=WZ(1)*XNK(:NMZ)
+      DO IZ=1,IELEM
+      DO IY=1,IELEM
+      DO IX=1,IELEM
+        II=IELEM**2*(IZ-1)+IELEM*(IY-1)+IX
+        IIX=IELEM*(IZ-1)+IY
+        IIY=IELEM*(IZ-1)+IX
+        IIZ=IELEM*(IY-1)+IX
+        ! X-SPACE
+        ! Assign I-boundary fluxes if lozenges A or B
+        IF((ILOZ.EQ.1).OR.(ILOZ.EQ.2))THEN
+        IF(MOD(IX,2).EQ.1) THEN
+          TMPXNI(IIX,J2,J)=TMPXNI(IIX,J2,J)+CST(IX)*WX(IX+1)
+     1                 *Q2(II,NM+1)
+        ELSE
+          TMPXNI(IIX,J2,J)=TMPXNI(IIX,J2,J)+CST(IX)*WX(IX+1)
+     1                 *Q2(II,NM+1)*SIGN(1.0,MUH)
+        ENDIF
+        ENDIF
+        ! Y-SPACE
+        ! Assign J-boundary fluxes if lozenges B or C
+        IF((ILOZ.EQ.2).OR.(ILOZ.EQ.3))THEN
+        IF(MOD(IY,2).EQ.1) THEN
+          TMPXNJ(IIY,I2,I)=TMPXNJ(IIY,I2,I)+CST(IY)*WY(IY+1)
+     1                 *Q2(II,NM+1)
+        ELSE
+          TMPXNJ(IIY,I2,I)=TMPXNJ(IIY,I2,I)+CST(IY)*WY(IY+1)
+     1                 *Q2(II,NM+1)*SIGN(1.0,ETAH)
+        ENDIF
+        ENDIF
+        ! D-SPACE
+        ! Assign D-boundary fluxes if lozenge A using XNJ
+        IF((ILOZ.EQ.1))THEN
+        I3=I_FETCH
+        IF(MOD(IY,2).EQ.1) THEN
+          TMPXND(IIY,I3,DCOORD)=TMPXND(IIY,I3,DCOORD)+CST(IY)*WY(IY+1)
+     1                 *Q2(II,NM+1)
+        ELSE
+          TMPXND(IIY,I3,DCOORD)=TMPXND(IIY,I3,DCOORD)+CST(IY)*WY(IY+1)
+     1                 *Q2(II,NM+1)*SIGN(1.0,ETAH)
+        ENDIF
+        ENDIF
+        ! Assign D-boundary fluxes if lozenge C using XNI
+        IF((ILOZ.EQ.3))THEN
+        IF(MOD(IX,2).EQ.1) THEN
+          TMPXND(IIX,J2,DCOORD)=TMPXND(IIX,J2,DCOORD)+CST(IX)*WX(IX+1)
+     1                 *Q2(II,NM+1)
+        ELSE
+          TMPXND(IIX,J2,DCOORD)=TMPXND(IIX,J2,DCOORD)+CST(IX)*WX(IX+1)
+     1                 *Q2(II,NM+1)*SIGN(1.0,MUH)
+        ENDIF
+        ENDIF
+        ! Z-SPACE
+        IF(MOD(IZ,2).EQ.1) THEN
+          TMPXNK(IIZ,I2,J2,ILOZ,IHEX_XY)=TMPXNK(IIZ,I2,J2,ILOZ,IHEX_XY)
+     1                 +CST(IZ)*WZ(IZ+1)*Q2(II,NM+1)
+        ELSE
+          TMPXNK(IIZ,I2,J2,ILOZ,IHEX_XY)=TMPXNK(IIZ,I2,J2,ILOZ,IHEX_XY)
+     1                 +CST(IZ)*WZ(IZ+1)*Q2(II,NM+1)*SIGN(1.0,XIH)
+        ENDIF
+      ENDDO
+      ENDDO
+      ENDDO
+      ! FLIP GRADIENTS IF NECESSARY
+      DO IZ=1,IELEM
+      DO IY=1,IELEM
+        IIX=IELEM*(IZ-1)+IY
+        IF((MOD(IY,2).EQ.0).AND.(ILOZ.EQ.3).AND.(IL.EQ.ISPLH))
+     1    TMPXND(IIX,J2,DCOORD)=TMPXND(IIX,J2,DCOORD)*(-1)
+      ENDDO
+      ENDDO
+      I3=I_FETCH
+      DO IZ=1,IELEM
+      DO IX=1,IELEM
+        IIY=IELEM*(IZ-1)+IX
+        IF((MOD(IX,2).EQ.0).AND.(ILOZ.EQ.1).AND.(JL.EQ.ISPLH))
+     1    TMPXND(IIY,I3,DCOORD)=TMPXND(IIY,I3,DCOORD)*(-1)
+      ENDDO
+      ENDDO
+      ! LFIXUP
+      IF(IELEM.EQ.1.AND.LFIXUP)THEN
+        IF((ILOZ.EQ.1).OR.(ILOZ.EQ.2))THEN
+          IF(TMPXNI(1,J2,J).LE.RLOG) TMPXNI(1,J2,J)=0.0
+        ELSE
+          IF(TMPXND(1,J2,DCOORD).LE.RLOG) TMPXND(1,J2,DCOORD)=0.0
+        ENDIF
+        IF((ILOZ.EQ.2).OR.(ILOZ.EQ.3))THEN
+          IF(TMPXNJ(1,I2,I).LE.RLOG) TMPXNJ(1,I2,I)=0.0
+        ELSE
+          I3=I_FETCH
+          IF(TMPXND(1,I3,DCOORD).LE.RLOG) TMPXND(1,I3,DCOORD)=0.0
+        ENDIF
+        IF(TMPXNK(1,I2,J2,ILOZ,IHEX_XY).LE.RLOG) 
+     1   TMPXNK(1,I2,J2,ILOZ,IHEX_XY)=0.0
+      ENDIF
+      WX=WX0
+      WY=WY0
+      WZ=WZ0
+
+      ! SAVE LEGENDRE MOMENT OF THE FLUX
+      IOF=((ILOZ-1)*ISPLH+(J2-1))*ISPLH+I2
+      DO P=1,NSCT
+      DO IEL=1,NM
+      FLUX(IEL,P,IOF,IHEX_XY,K) = FLUX(IEL,P,IOF,IHEX_XY,K)
+     1   +Q2(IEL,NM+1)*DN(P,M)
+      ENDDO
+      ENDDO
+
+      ENDDO ! END OF WITHIN LOZENGE J-LOOP
+      ENDDO ! END OF WITHIN LOZENGE I-LOOP
+
+      ENDDO ! END OF LOZENGE LOOP
+
+      ENDDO ! END OF I COLUMNS OF CARTESIAN MAP LOOP
+      ENDDO ! END OF J COLUMNS OF CARTESIAN MAP LOOP
+      ENDDO ! END OF Z-LOOP
+
+      ! SAVE FLUX INFORMATION
+      FLUX_G(:,:,:,:,:,IG)=FLUX_G(:,:,:,:,:,IG)+FLUX(:,:,:,:,:)
+
+      ! SAVE K-BOUNDARY CONDITIONS IF NOT VOID B.C.
+      IF((NCODE(5).NE.1).or.(NCODE(6).NE.1))THEN
+        IOF=(M-1)*(L5)
+        FUNKNO(LFLX+IOF+1:LFLX+IOF+L5,IG)=
+     >  RESHAPE(REAL(TMPXNK(:NMZ,:ISPLH,:ISPLH,:3,:NHEX)),
+     >  (/L5/))
+      ENDIF
+
+      ENDDO ! END OF DIRECTION LOOP
+      ENDDO ! END OF ENERGY LOOP
+*$OMP END PARALLEL DO
+      ENDDO ! END OF OCTANT LOOP
+
+      ! SAVE FLUX INFORMATION
+      DO IG=1,NGEFF
+        IF(.NOT.INCONV(IG)) CYCLE
+        FUNKNO(:LFLX,IG)=
+     1  RESHAPE(REAL(FLUX_G(:NM,:NSCT,:3*ISPLH**2,:NHEX,:LZ,IG)),
+     2  (/LFLX/))
+      ENDDO
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(FLUX_G,FLUX,INDANG,TMPXNI,TMPXNJ,TMPXND,TMPXNK,TMPMAT)
+      RETURN
+  500 FORMAT(16H SNFBH3: thread=,I8,12H --->(group=,I4,7H angle=,I4,1H))
+      END