path: root/Dragon/src/SNFKH3.F

*DECK SNFKH3
      SUBROUTINE SNFKH3(NKBA,NUN,NGEFF,IMPX,INCONV,NGIND,NHEX,LZ,ISPLH,
     1 SIDE,IELEM,NM,NMX,NMY,NMZ,NMAT,NPQ,NSCT,MAT,VOL,TOTAL,NCODE,
     2 ZCODE,QEXT,LFIXUP,DU,DE,DZ,W,MRMZ,DC,DB,DA,MN,DN,WX,WY,WZ,CST,
     3 ISADPT,LOZSWP,COORDMAP,FUNKNO)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Perform one inner iteration for solving SN equations in 3D Cartesian
* geometry for the HODD method. KBA-like 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
* NKBA     number of macrocells along the z-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.
* 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),ZCODE(6),TOTAL(0:NMAT,NGEFF),
     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(:) :: III,JJJ,KKK
      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

      ! CALL XABORT('SNFKH3: testing 1')
*----
*  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 IZ=1,LZ
        DO IHEX_XY=1,NHEX
          TMPMAT(:,:,:,COORDMAP(1,IHEX_XY),COORDMAP(2,IHEX_XY),IZ) = 
     >         MAT(:,:,:,IHEX_XY,IZ)
        ENDDO
      ENDDO
*----
*  SCRATCH STORAGE ALLOCATION
*----
      ALLOCATE(INDANG(NPQ,12))
      ALLOCATE(FLUX(NM,NSCT))
      ALLOCATE(FLUX_G(NM,NSCT,3*ISPLH**2,NHEX,LZ,NGEFF))
      ALLOCATE(TMPXNI(NMX,ISPLH,NCOL,LZ,NPQ,NGEFF))
      ALLOCATE(TMPXNJ(NMY,ISPLH,NCOL,LZ,NPQ,NGEFF))
      ALLOCATE(TMPXND(NMX,ISPLH,NCOL,LZ,NPQ,NGEFF))
      ALLOCATE(TMPXNK(NMZ,ISPLH,ISPLH,3,NHEX,NPQ,NGEFF))
*----
*  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
*----
*  NUMBER OF MACROCELLS (MACRO*)
*  NUMBER OF LZ LAYERS IN EACH MACROCELL (NCELL*)
*----
      MACROZ=NKBA
      NCELLZ = (LZ + MACROZ - 1) / MACROZ
*----
*  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

      TMPXNI(:NMX,:ISPLH,:NCOL,:LZ,:NPQ,:NGEFF)=0.0D0
      TMPXNJ(:NMY,:ISPLH,:NCOL,:LZ,:NPQ,:NGEFF)=0.0D0
      TMPXND(:NMX,:ISPLH,:NCOL,:LZ,:NPQ,:NGEFF)=0.0D0
      TMPXNK(:NMZ,:ISPLH,:ISPLH,:3,:NHEX,:NPQ,:NGEFF)=0.0D0

*----
*  LOOP OVER WAVEFRONTS 
*----
      DO IDI=1,NCOL+NCOL+MACROZ-2
*----
*  SET SWEEP INDICES
*----
      NMAX=MIN(NCOL,IDI)*MIN(NCOL,IDI)
      ALLOCATE(III(NMAX),JJJ(NMAX),KKK(NMAX))
      NCEL=0

      J_STT=MAX(1,IDI-NCOL-MACROZ+2)
      J_END=MIN(NCOL,IDI)

      DO J_MC=J_STT,J_END
      JIM=J_MC
      ! 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

      IF((IND_XY.EQ.1).OR.(IND_XY.EQ.4)) THEN
        I_STT=MAX(1,IDI-J_MC-MACROZ+2)
        I_END=MIN(NCOL,IDI-J_MC+1)
      ELSE
        I_STT=MAX(1,NRINGS-J_MC+1)
        I_END=MIN(NCOL,IDI-J_MC+(NRINGS+1-J_MC))
      ENDIF

      DO I_MC=I_STT,I_END
      IIM=I_MC
      ! 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,1).EQ.-1) CYCLE

      ! Find I coordinate of macrocell
      IF((IND_XY.EQ.1).OR.(IND_XY.EQ.4)) THEN
        K_MC=IDI-I_MC-J_MC+2
      ELSE
        K_MC=IDI-I_MC+NRINGS-((J_MC-1)*2)
      ENDIF
      IF((K_MC.GT.MACROZ)) CYCLE
      K=K_MC

      NCEL=NCEL+1
      IF(NCEL.GT.NMAX) CALL XABORT('SNFDH3: NMAX OVERFLOW.')
      III(NCEL)=I
      JJJ(NCEL)=J
      KKK(NCEL)=K
      ENDDO ! I_MC

      ENDDO ! J_MC
*
*----
*  MAIN SWAPPING LOOPS FOR SN FLUX CALCULATION
*  LOOP OVER MACROCELLS IN WAVEFRONT AND ENERGY
*----
*$OMP  PARALLEL DO 
*$OMP+ PRIVATE(ITID,FLUX,M,IG,XNI,XNJ,XNK,Q,Q2,IOF,IER,II,JJ,I,J,K,K_MC)
*$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)
*$OMP+ SHARED(IND,TMPXNI,TMPXNJ,TMPXND,TMPXNK,III,JJJ,KKK,FUNKNO)
*$OMP+ FIRSTPRIVATE(WX,WY,WZ,WX0,WY0,WZ0)
*$OMP+ COLLAPSE(2)

      ! LOOP FOR MACROCELLS TO EXECUTE IN PARALLEL
      DO ICEL=1,NCEL

      ! 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

      I=III(ICEL)
      J=JJJ(ICEL)
      K_MC=KKK(ICEL)

      ! 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('SNFDH3: IHEX_XY FAILURE.')
      ! Find D coordinate
      DCOORD = ABS(COORDMAP(3,IHEX_XY))-NRINGS

      IF(IDI.EQ.1)THEN
      ! 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,IPQD,IG)=
     >    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,IPQD,IG)=
     >    TMPXNK(:NMZ,:ISPLH,:ISPLH,:3,:NHEX,IPQD,IG)*ZCODE(5)
      ELSE
        TMPXNK(:NMZ,:ISPLH,:ISPLH,:3,:NHEX,IPQD,IG)=
     >    TMPXNK(:NMZ,:ISPLH,:ISPLH,:3,:NHEX,IPQD,IG)*ZCODE(6)
      ENDIF
      ENDIF

*----
*  LOOP OVER Z-AXIS PLANES IN MACROCELL
*----

      DO IMZ=1,MIN(NCELLZ,LZ-(K_MC-1)*NCELLZ)
      K=(K_MC-1)*NCELLZ+IMZ
      IF(IND.GE.7) K=LZ+1-K

*----
*  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('SNFKH3: 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

      FLUX(:IELEM**3,:NSCT)=0.0D0

      ! 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,K,IPQD,IG)
        XNJ(:) = TMPXND(:,I_FETCH,DCOORD,K,IPQD,IG)
      ELSEIF((ILOZ.EQ.2))THEN
        XNI(:) = TMPXNI(:,J2,J,K,IPQD,IG)
        XNJ(:) = TMPXNJ(:,I2,I,K,IPQD,IG)
      ELSEIF((ILOZ.EQ.3))THEN
        XNI(:) = TMPXND(:,J2,DCOORD,K,IPQD,IG)
        XNJ(:) = TMPXNJ(:,I2,I,K,IPQD,IG)
      ENDIF
      XNK(:) = TMPXNK(:,I2,J2,ILOZ,IHEX_XY,IPQD,IG)

      ! 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('SNFKH3: 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('SNFKH3: 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('SNFKH3: 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('SNFKH3: testing 69 ')
      ! CALL XABORT('SNFKH3: 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('SNFKH3: testing 70 ')
      ! CALL XABORT('SNFKH3: testing 191 ')
        ! UPPER DIAGONAL TERMS
        ELSEIF(II.LT.JJ) THEN
          ! CALL XABORT('SNFKH3: 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('SNFKH3: testing 75 ')
      ! CALL XABORT('SNFKH3: testing 19 ')

         ! UNDER DIAGONAL TERMS
        ELSE
          ! CALL XABORT('SNFKH3: 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('SNFKH3: 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,K,IPQD,IG)=WX(1)*XNI(:NMX)
      ELSE
        TMPXND(:NMX,J2,DCOORD,K,IPQD,IG)=WX(1)*XNI(:NMX)
      ENDIF
      IF((ILOZ.EQ.2).OR.(ILOZ.EQ.3))THEN
        TMPXNJ(:NMY,I2,I,K,IPQD,IG)=WY(1)*XNJ(:NMY)
      ELSE
        I3=I_FETCH
        TMPXND(:NMY,I3,DCOORD,K,IPQD,IG)=WY(1)*XNJ(:NMY)
      ENDIF
      TMPXNK(:NMZ,I2,J2,ILOZ,IHEX_XY,IPQD,IG)=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,K,IPQD,IG)=TMPXNI(IIX,J2,J,K,IPQD,IG)+
     1      CST(IX)*WX(IX+1)*Q2(II,NM+1)
        ELSE
          TMPXNI(IIX,J2,J,K,IPQD,IG)=TMPXNI(IIX,J2,J,K,IPQD,IG)+
     1      CST(IX)*WX(IX+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,K,IPQD,IG)=TMPXNJ(IIY,I2,I,K,IPQD,IG)+
     1      CST(IY)*WY(IY+1)*Q2(II,NM+1)
        ELSE
          TMPXNJ(IIY,I2,I,K,IPQD,IG)=TMPXNJ(IIY,I2,I,K,IPQD,IG)+
     1      CST(IY)*WY(IY+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,K,IPQD,IG)=
     1      TMPXND(IIY,I3,DCOORD,K,IPQD,IG)+CST(IY)*WY(IY+1)*Q2(II,NM+1)
        ELSE
          TMPXND(IIY,I3,DCOORD,K,IPQD,IG)=
     1      TMPXND(IIY,I3,DCOORD,K,IPQD,IG)+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,K,IPQD,IG)=
     1      TMPXND(IIX,J2,DCOORD,K,IPQD,IG)+CST(IX)*WX(IX+1)*Q2(II,NM+1)
        ELSE
          TMPXND(IIX,J2,DCOORD,K,IPQD,IG)=
     1      TMPXND(IIX,J2,DCOORD,K,IPQD,IG)+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,IPQD,IG)=
     1      TMPXNK(IIZ,I2,J2,ILOZ,IHEX_XY,IPQD,IG)+
     1      CST(IZ)*WZ(IZ+1)*Q2(II,NM+1)
        ELSE
          TMPXNK(IIZ,I2,J2,ILOZ,IHEX_XY,IPQD,IG)=
     1      TMPXNK(IIZ,I2,J2,ILOZ,IHEX_XY,IPQD,IG)+
     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,K,IPQD,IG)=
     1    TMPXND(IIX,J2,DCOORD,K,IPQD,IG)*(-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,K,IPQD,IG)=
     1    TMPXND(IIY,I3,DCOORD,K,IPQD,IG)*(-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,K,IPQD,IG).LE.RLOG) 
     1      TMPXNI(1,J2,J,K,IPQD,IG)=0.0
        ELSE
          IF(TMPXND(1,J2,DCOORD,K,IPQD,IG).LE.RLOG) 
     1      TMPXND(1,J2,DCOORD,K,IPQD,IG)=0.0
        ENDIF
        IF((ILOZ.EQ.2).OR.(ILOZ.EQ.3))THEN
          IF(TMPXNJ(1,I2,I,K,IPQD,IG).LE.RLOG) 
     1      TMPXNJ(1,I2,I,K,IPQD,IG)=0.0
        ELSE
          I3=I_FETCH
          IF(TMPXND(1,I3,DCOORD,K,IPQD,IG).LE.RLOG) 
     1      TMPXND(1,I3,DCOORD,K,IPQD,IG)=0.0
        ENDIF
        IF(TMPXNK(1,I2,J2,ILOZ,IHEX_XY,IPQD,IG).LE.RLOG) 
     1   TMPXNK(1,I2,J2,ILOZ,IHEX_XY,IPQD,IG)=0.0
      ENDIF
      WX=WX0
      WY=WY0
      WZ=WZ0

      ! 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

      IOF=((ILOZ-1)*ISPLH+(J2-1))*ISPLH+I2
      FLUX_G(:,:,IOF,IHEX_XY,K,IG)=FLUX_G(:,:,IOF,IHEX_XY,K,IG)+
     1  FLUX(:,:)

      ENDDO ! END OF WITHIN LOZENGE J-LOOP
      ENDDO ! END OF WITHIN LOZENGE I-LOOP

      ENDDO ! END OF LOZENGE LOOP
      ENDDO ! END OF Z-LOOP

      ! SAVE K-BOUNDARY CONDITIONS IF NOT VOID B.C.
      IF(IDI.EQ.NCOL+NCOL+MACROZ-2)THEN
      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,IPQD,IG)),
     >  (/L5/))
      ENDIF
      ENDIF

      ENDDO ! END OF DIRECTION LOOP
      ENDDO ! END OF ENERGY LOOP
      ENDDO ! END OF MACROCELL (WITHIN ONE WAVEFRONT) LOOP
*$OMP END PARALLEL DO
      DEALLOCATE(JJJ,III,KKK)
      ENDDO ! END OF WAVEFRONT LOOP
      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 SNFKH3: thread=,I8,12H --->(group=,I4,7H angle=,I4,1H))
      END