Initial commit from Polytechnique Montreal

1 files changed, 1129 insertions, 0 deletions

diff --git a/Dragon/src/SNFLUX.f b/Dragon/src/SNFLUX.f
new file mode 100644
index 0000000..09a832d
--- /dev/null
+++ b/Dragon/src/SNFLUX.f
@@ -0,0 +1,1129 @@
+*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