path: root/Dragon/src/SNTT3D.f

*DECK SNTT3D
      SUBROUTINE SNTT3D(IGE,IMPX,LX,LY,LZ,SIDE,IELEM,NLF,NPQ,NSCT,IQUAD,
     1 NCODE,ZCODE,MAT,XXX,YYY,ZZZ,VOL,IDL,DU,DE,DZ,W,MRMX,MRMY,MRMZ,DC,
     2 DB,DA,PL,LL4,NUN,EELEM,WX,WE,CST,IBFP,ISCHM,ESCHM,IGLK,MN,DN,IL,
     3 IM,ISCAT)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Numbering corresponding to a 3-D Cartesian with discrete ordinates
* approximation of the flux.
*
*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): N. Martin and C. Bienvenue
*
*Parameters: input
* IGE     geometry type with (=0 for Cartesian, =2 for hexagonal).
* IMPX    print parameter.
* LX      number of elements along the X axis;
*         OR, number of elements in X-Y plane in hexagonal geometry,
*         including lozenges and further submeshing of lozenges
*         e.g. domain of 7 hex. by 5 levels with a submeshing of 2
*         will have LX=7*3*2*2=84 
* LY      number of elements along the Y axis.
* LZ      number of elements along the Z axis.
* SIDE    side of hexagon.
* IELEM   measure of order of the spatial approximation polynomial:
*         =1 constant - only for HODD, classical diamond scheme 
*         (default for HODD);
*         =2 linear - default for DG;
*         =3 parabolic;
*         =4 cubic - only for DG.
* NLF     SN order for the flux (even number).
* NPQ     number of SN directions in eight octants.
* NSCT    maximum number of spherical harmonics moments of the flux.
* IQUAD   type of SN quadrature (1 Level symmetric, type IQUAD;
*         4 Legendre-Chebyshev; 5 symmetric Legendre-Chebyshev;
*         6 quadruple range).
* NCODE   type of boundary condition applied on each side
*         (i=1 X-;  i=2 X+;  i=3 Y-;  i=4 Y+;  i=5 Z-;  i=6 Z+):
*         =1 VOID;   =2 REFL;  =4 TRAN.
* ZCODE   ZCODE(I) is the albedo corresponding to boundary condition
*         'VOID' on each side (ZCODE(I)=0.0 by default).
* MAT     mixture index assigned to each element.
* XXX     Cartesian coordinates along the X axis.
* YYY     Cartesian coordinates along the Y axis.
* ZZZ     Cartesian coordinates along the Z axis.
* EELEM   measure of order of the energy approximation polynomial:
*         =1 constant - default for HODD;
*         =2 linear - default for DG;
*         >3 higher orders.
* IBFP    type of energy proparation relation:
*         =0 no Fokker-Planck term;
*         =1 Galerkin type;
*         =2 heuristic Przybylski and Ligou type.
* ISCHM   method of spatial discretisation:
*         =1 High-Order Diamond Differencing (HODD) - default;
*         =2 Discontinuous Galerkin finite element method (DG);
*         =3 Adaptive weighted method (AWD).
* ESCHM   method of energy discretisation:
*         =1 High-Order Diamond Differencing (HODD) - default;
*         =2 Discontinuous Galerkin finite element method (DG);
*         =3 Adaptive weighted method (AWD).
* IGLK    angular interpolation type:
*         =0 classical SN method.
*         =1 Galerkin quadrature method (M = inv(D))
*         =2 Galerkin quadrature method (D = inv(M))
* ISCAT   maximum number of spherical harmonics moments of the flux.
*
*Parameters: output
* VOL     volume of each element.
* IDL     isotropic flux indices.
* DU      first direction cosines ($\\mu$).
* DE      second direction cosines ($\\eta$).
* DZ      third direction cosines ($\\xi$).
* W       weights.
* MRMX    quadrature index.
* MRMY    quadrature index.
* MRMZ    quadrature index.
* DC      diamond-scheme parameter.
* DB      diamond-scheme parameter.
* DA      diamond-scheme parameter.
* PL      discrete values of the spherical harmonics corresponding
*         to the 3D SN quadrature.
* LL4     number of unknowns being solved for, over the domain. This 
*         includes the various moments of the isotropic (and if present,
*         anisotropic) flux. 
* NUN     total number of unknowns stored in the FLUX vector per group.
*         This includes LL4 (see above) as well as any surface boundary
*         fluxes, if present.
* WX      spatial closure relation weighting factors.
* WE      energy closure relation weighting factors.
* CST     constants for the polynomial approximations.
* MN      moment-to-discrete matrix.
* DN      discrete-to-moment matrix.
* IL      indexes (l) of each spherical harmonics in the
*         interpolation basis.
* IM      indexes (m) of each spherical harmonics in the
*         interpolation basis.
*
*-----------------------------------------------------------------------
*
*----
*  SUBROUTINE ARGUMENTS
*----
      INTEGER IMPX,LX,LY,LZ,IELEM,NLF,NPQ,NSCT,IQUAD,NCODE(6),
     1 MAT(LX,LY,LZ),IDL(LX*LY*LZ),MRMX(NPQ),MRMY(NPQ),MRMZ(NPQ),EELEM,
     2 IBFP,ISCHM,ESCHM,IL(NSCT),IM(NSCT),ISCAT,IGLK
      REAL ZCODE(6),VOL(LX,LY,LZ),XXX(LX+1),YYY(LY+1),ZZZ(LZ+1),
     1 DU(NPQ),DE(NPQ),DZ(NPQ),W(NPQ),DC(LX,LY,NPQ),DB(LX,LZ,NPQ),
     2 DA(LY,LZ,NPQ),PL(NSCT,NPQ),WX(IELEM+1),
     3 WE(EELEM+1),CST(MAX(IELEM,EELEM)),MN(NPQ,NSCT),DN(NSCT,NPQ)
*----
*  LOCAL VARIABLES
*----
      CHARACTER HSMG*131
      LOGICAL L1,L2,L3,L4,L5,L6
      PARAMETER(RLOG=1.0E-8,PI=3.141592654)
      REAL PX,PE
      DOUBLE PRECISION MND(NPQ,NPQ),NORM,IPROD
      INTEGER, ALLOCATABLE, DIMENSION(:) :: JOP
      REAL, ALLOCATABLE, DIMENSION(:) :: XX,YY,ZZ,UU,WW,TPQ,UPQ,VPQ,WPQ
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: V,V2
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: U
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: RLM
*----
*  SCRATCH STORAGE ALLOCATION
*----
      ALLOCATE(XX(LX),YY(LY),ZZ(LZ))
*----
*  UNFOLD HEIGHT-OCTANT QUADRATURES.
*----
      IF(MOD(NLF,2).EQ.1) CALL XABORT('SNTT3D: EVEN NLF EXPECTED.')
      IF(IQUAD.EQ.10) THEN
         NPQ0=NLF**2/4
      ELSE
         NPQ0=NLF*(NLF/2+1)/4
      ENDIF
      ALLOCATE(JOP(NLF/2),UU(NLF/2),WW(NLF/2),TPQ(NPQ0),UPQ(NPQ0),
     1 VPQ(NPQ0),WPQ(NPQ0))
      IF(IQUAD.EQ.1) THEN
!Level-symmetric quadrature of type 1
         CALL SNQU01(NLF,JOP,UU,WW,TPQ,UPQ,VPQ,WPQ)
      ELSE IF(IQUAD.EQ.2) THEN
!Level-symmetric quadrature of type 2
         CALL SNQU02(NLF,JOP,UU,WW,TPQ,UPQ,VPQ,WPQ)
      ELSE IF(IQUAD.EQ.3) THEN
!Snow Level-symmetric type quadrature 
         CALL SNQU03(NLF,JOP,UU,WW,TPQ,UPQ,VPQ,WPQ)
      ELSE IF(IQUAD.EQ.4) THEN
!Legendre-Chebyshev quadrature 
         CALL SNQU04(NLF,JOP,UU,WW,TPQ,UPQ,VPQ,WPQ)
      ELSE IF(IQUAD.EQ.5) THEN
!Symmetric Legendre-Chebyshev quadrature
         CALL SNQU05(NLF,TPQ,UPQ,VPQ,WPQ)
      ELSE IF(IQUAD.EQ.6) THEN
!Quadruple Range quadrature
         CALL SNQU06(NLF,TPQ,UPQ,VPQ,WPQ)
      ELSE IF(IQUAD.EQ.10) THEN
         CALL SNQU10(NLF,JOP,UU,WW,TPQ,UPQ,VPQ,WPQ)
      ELSE
         CALL XABORT('SNTT3D: UNKNOWN QUADRATURE TYPE.')
      ENDIF
      N=0
      IOF=0
      DO 320 I=1,NLF/2
         JOF=IOF+NLF-2*I+2
         DO 330 J=0,NLF/2-I
            IOF=IOF+1
            KOF=IOF+(NLF+2)*NLF/4
            LOF=IOF+(NLF+2)*NLF/2
            MRMX(IOF)=JOF
            MRMY(IOF)=KOF
            MRMZ(IOF)=LOF
            DU(IOF)=-UPQ(N+J+1)
            DE(IOF)=-VPQ(N+J+1)
            DZ(IOF)=-TPQ(N+J+1)
            W(IOF)=WPQ(N+J+1)
            JOF=JOF-1
 330     CONTINUE
         DO 340 J=NLF/2-I,0,-1
            IOF=IOF+1
            KOF=IOF+(NLF+2)*NLF/4
            LOF=IOF+(NLF+2)*NLF/2
            MRMX(IOF)=JOF
            MRMY(IOF)=KOF
            MRMZ(IOF)=LOF
            DU(IOF)=UPQ(N+J+1)
            DE(IOF)=-VPQ(N+J+1)
            DZ(IOF)=-TPQ(N+J+1)
            W(IOF)=WPQ(N+J+1)
            JOF=JOF-1
 340     CONTINUE
         N=N+NLF/2-I+1
 320  CONTINUE
      N=0
      DO 350 I=1,NLF/2
         JOF=IOF+NLF-2*I+2
         DO 360 J=0,NLF/2-I
            IOF=IOF+1
            KOF=IOF-(NLF+2)*NLF/4
            LOF=IOF+(NLF+2)*NLF/2
            MRMX(IOF)=JOF
            MRMY(IOF)=KOF
            MRMZ(IOF)=LOF
            DU(IOF)=-UPQ(N+J+1)
            DE(IOF)=VPQ(N+J+1)
            DZ(IOF)=-TPQ(N+J+1)
            W(IOF)=WPQ(N+J+1)
            JOF=JOF-1
 360     CONTINUE
         DO 370 J=NLF/2-I,0,-1
            IOF=IOF+1
            KOF=IOF-(NLF+2)*NLF/4
            LOF=IOF+(NLF+2)*NLF/2
            MRMX(IOF)=JOF
            MRMY(IOF)=KOF
            MRMZ(IOF)=LOF
            DU(IOF)=UPQ(N+J+1)
            DE(IOF)=VPQ(N+J+1)
            DZ(IOF)=-TPQ(N+J+1)
            W(IOF)=WPQ(N+J+1)
            JOF=JOF-1
 370     CONTINUE
         N=N+NLF/2-I+1
 350  CONTINUE
      N=0
      DO 380 I=1,NLF/2
         JOF=IOF+NLF-2*I+2
         DO 390 J=0,NLF/2-I
            IOF=IOF+1
            KOF=IOF+(NLF+2)*NLF/4
            LOF=IOF-(NLF+2)*NLF/2
            MRMX(IOF)=JOF
            MRMY(IOF)=KOF
            MRMZ(IOF)=LOF
            DU(IOF)=-UPQ(N+J+1)
            DE(IOF)=-VPQ(N+J+1)
            DZ(IOF)=TPQ(N+J+1)
            W(IOF)=WPQ(N+J+1)
            JOF=JOF-1
 390     CONTINUE
         DO 400 J=NLF/2-I,0,-1
            IOF=IOF+1
            KOF=IOF+(NLF+2)*NLF/4
            LOF=IOF-(NLF+2)*NLF/2
            MRMX(IOF)=JOF
            MRMY(IOF)=KOF
            MRMZ(IOF)=LOF
            DU(IOF)=UPQ(N+J+1)
            DE(IOF)=-VPQ(N+J+1)
            DZ(IOF)=TPQ(N+J+1)
            W(IOF)=WPQ(N+J+1)
            JOF=JOF-1
 400     CONTINUE
         N=N+NLF/2-I+1
 380  CONTINUE
      N=0
      DO 410 I=1,NLF/2
         JOF=IOF+NLF-2*I+2
         DO 420 J=0,NLF/2-I
            IOF=IOF+1
            KOF=IOF-(NLF+2)*NLF/4
            LOF=IOF-(NLF+2)*NLF/2
            MRMX(IOF)=JOF
            MRMY(IOF)=KOF
            MRMZ(IOF)=LOF
            DU(IOF)=-UPQ(N+J+1)
            DE(IOF)=VPQ(N+J+1)
            DZ(IOF)=TPQ(N+J+1)
            W(IOF)=WPQ(N+J+1)
            JOF=JOF-1
 420     CONTINUE
         DO 430 J=NLF/2-I,0,-1
            IOF=IOF+1
            KOF=IOF-(NLF+2)*NLF/4
            LOF=IOF-(NLF+2)*NLF/2
            MRMX(IOF)=JOF
            MRMY(IOF)=KOF
            MRMZ(IOF)=LOF
            DU(IOF)=UPQ(N+J+1)
            DE(IOF)=VPQ(N+J+1)
            DZ(IOF)=TPQ(N+J+1)
            W(IOF)=WPQ(N+J+1)
            JOF=JOF-1
 430     CONTINUE
         N=N+NLF/2-I+1
 410  CONTINUE
      DEALLOCATE(WPQ,VPQ,UPQ,TPQ,WW,UU,JOP)
      IF(IMPX.GE.4) THEN
         WRITE(6,'(/41H SNTT3D:HEIGHT-OCTANT ANGULAR QUADRATURES:/26X,
     1   2HMU,9X,3HETA,10X,2HXI,6X,6HWEIGHT)')
         SUM=0.0
         DO 70 N=1,NPQ
         SUM=SUM+W(N)
         WRITE(6,'(1X,4I5,1P,4E12.4)') N,MRMX(N),MRMY(N),MRMZ(N),DU(N),
     1   DE(N),DZ(N),W(N)
   70    CONTINUE
         WRITE(6,'(54X,10(1H-)/52X,1P,E12.4)') SUM
      ENDIF
*----
*  IDENTIFICATION OF THE GEOMETRY.
*----
      IF(IGE.EQ.0) THEN
* ----------
*        3D CARTESIAN
* ----------
         DO 83 N=1,NPQ
         VU=DU(N)
         VE=DE(N)
         VZ=DZ(N)
         DO 82 I=1,LX
         DO 81 J=1,LY
         DO 80 K=1,LZ
         XX(I)=XXX(I+1)-XXX(I)
         YY(J)=YYY(J+1)-YYY(J)
         ZZ(K)=ZZZ(K+1)-ZZZ(K)
         DA(J,K,N)=VU*YY(J)*ZZ(K)
         DB(I,K,N)=VE*XX(I)*ZZ(K)
         DC(I,J,N)=VZ*XX(I)*YY(J)
         VOL(I,J,K)=XX(I)*YY(J)*ZZ(K)
   80    CONTINUE
   81    CONTINUE
   82    CONTINUE
   83    CONTINUE
      ELSEIF(IGE.EQ.2) THEN
* ----------
*        3D HEXAGONAL
* ----------
         DET = SQRT(3.0)*(SIDE**2)/2.0 
         DO 93 N=1,NPQ
         VU=DU(N)
         VE=DE(N)
         VZ=DZ(N)
         DO 92 K=1,LZ
         DO 91 J=1,LY
         DO 90 I=1,LX
         ZZ(K)=ZZZ(K+1)-ZZZ(K)
         DA(J,K,N)=VU*ZZ(K)
         DB(I,K,N)=VE*ZZ(K)
         DC(I,J,N)=VZ*DET
         VOL(I,J,K)=DET*ZZ(K)
   90    CONTINUE
   91    CONTINUE
   92    CONTINUE
   93    CONTINUE
      ENDIF
*----
*  GENERATE SPHERICAL HARMONICS FOR SCATTERING SOURCE.
*----
      IOF=0
      DO 211 L=0,ISCAT-1
      DO 210 M=-L,L   
      IOF=IOF+1
      IF(IOF.GT.ISCAT) GO TO 211
      DO 200 N=1,NPQ
      PL(IOF,N)=PNSH(L,M,DU(N),DE(N),DZ(N))
  200 CONTINUE
  210 CONTINUE
  211 CONTINUE
*----
*  GENERATE MAPPING MATRIX FOR GALERKIN QUADRATURE METHOD
*----
      MN(:NPQ,:NSCT)=0.0
      DN(:NSCT,:NPQ)=0.0
      IL(:NSCT)=0
      IM(:NSCT)=0
      IF(IGLK.NE.0) THEN
         ALLOCATE(U(NPQ,NPQ),RLM(NPQ,ISCAT,2*ISCAT-1),V(NPQ),V2(NPQ))
         RLM(:NPQ,:ISCAT-1,:2*ISCAT-1)=0.0
         DO L=0,ISCAT-1
         DO M=-L,L
         DO N=1,NPQ
            RLM(N,L+1,M+L+1)=PNSH(L,M,DU(N),DE(N),DZ(N))
         ENDDO
         ENDDO
         ENDDO
         ! GRAM-SCHMIDT PROCEDURE TO FIND INDEPENDANT SET
         ! OF SPHERICAL HARMONICS WITH ANY QUADRATURE
         U(:NPQ,:NPQ)=0.0D0
         NORM=0.0D0
         DO N=1,NPQ
            NORM=NORM+RLM(N,1,1)**2
         ENDDO
         NORM=SQRT(NORM)
         DO N=1,NPQ 
            IF(IGLK.EQ.1) THEN
               MND(1,N)=W(N)*RLM(N,1,1)
            ELSEIF(IGLK.EQ.2) THEN
               MND(N,1)=(2.0*L+1.0)/(4.0*PI)*RLM(N,1,1)
            ELSE
               CALL XABORT('UNKNOWN GALERKIN QUADRATURE METHOD.')
            ENDIF
            U(N,1)=RLM(N,1,1)/NORM
         ENDDO
         IND=1
         ! ITERATE OVER THE SPHERICAL HARMONICS
         DO 212 L=1,ISCAT-1
         DO 213 M=-L,L
         V2(:NPQ)=0.0D0
         DO N=1,IND
         IPROD=0.0D0
         DO N2=1,NPQ
         IPROD=IPROD+U(N2,N)*RLM(N2,L+1,M+L+1)
         ENDDO
         DO N2=1,NPQ
         V2(N2)=V2(N2)+IPROD*U(N2,N)
         ENDDO
         ENDDO
         V(:NPQ)=0.0D0
         DO N=1,NPQ
         V(N)=RLM(N,L+1,M+L+1)-V2(N)
         ENDDO
         NORM=0.0D0
         DO N=1,NPQ
            NORM=NORM+V(N)**2
         ENDDO
         NORM=SQRT(NORM)
         ! KEEP THE SPHERICAL HARMONICS IF IT IS INDEPENDANT
         IF(NORM.GE.1.0E-5) THEN
            IND=IND+1
            DO N=1,NPQ
               U(N,IND)=V(N)/NORM
               IF(IGLK.EQ.1) THEN
                  MND(IND,N)=W(N)*RLM(N,L+1,M+L+1)
               ELSEIF(IGLK.EQ.2) THEN
                  MND(N,IND)=(2.0*L+1.0)/(4.0*PI)*RLM(N,L+1,M+L+1)
               ELSE
                  CALL XABORT('UNKNOWN GALERKIN QUADRATURE METHOD.')
               ENDIF
            ENDDO
            IL(IND)=L
            IM(IND)=M
         ENDIF
         IF(IND.EQ.NPQ) GOTO 217
  213    ENDDO
  212    ENDDO
         CALL XABORT('SNTT3D: THE'// 
     1   ' GRAM-SCHMIDTH PROCEDURE TO FIND A SUITABLE INTERPOLATION'//
     2   ' BASIS REQUIRE HIGHER LEGENDRE ORDER.')
         ! FIND INVERSE MATRIX    
  217    IF(IGLK.EQ.1) THEN    
            DN=REAL(MND)
            CALL ALINVD(NPQ,MND,NPQ,IER)
            IF(IER.NE.0) CALL XABORT('SNTT3D: SINGULAR MATRIX.')
            MN=REAL(MND)
         ELSEIF(IGLK.EQ.2) THEN
            MN=REAL(MND)
            CALL ALINVD(NPQ,MND,NPQ,IER)
            IF(IER.NE.0) CALL XABORT('SNTT3D: SINGULAR MATRIX.')
            DN=REAL(MND)
         ELSE
            CALL XABORT('UNKNOWN GALERKIN QUADRATURE METHOD.')
         ENDIF
         DEALLOCATE(U,RLM,V,V2)
      ELSE
         IND=1
         DO L=0,ISCAT-1
         DO 218 M=-L,L
         IL(IND)=L
         IM(IND)=M
         DO N=1,NPQ
         DN(IND,N)=W(N)*PNSH(L,M,DU(N),DE(N),DZ(N))
         MN(N,IND)=(2.0*L+1.0)/(4.0*PI)
     1   *PNSH(L,M,DU(N),DE(N),DZ(N))
         ENDDO
         IND=IND+1
  218    ENDDO
         ENDDO
      ENDIF
*----
* GENERATE THE WEIGHTING PARAMETERS OF THE CLOSURE RELATION.
*----
      PX=1
      PE=1
      IF(ISCHM.EQ.1.OR.ISCHM.EQ.3) THEN
        PX=1
      ELSEIF(ISCHM.EQ.2) THEN
        PX=0
      ELSE
        CALL XABORT('SNTT3D: UNKNOWN TYPE OF SPATIAL CLOSURE RELATION.')
      ENDIF
      IF(MOD(IELEM,2).EQ.1) THEN
        WX(1)=-PX
        WX(2:IELEM+1:2)=1+PX
        IF(IELEM.GE.2) WX(3:IELEM+1:2)=1-PX
      ELSE
        WX(1)=PX
        WX(2:IELEM+1:2)=1-PX
        IF(IELEM.GE.2) WX(3:IELEM+1:2)=1+PX
      ENDIF  
      IF(IBFP.NE.0) THEN
      IF(ESCHM.EQ.1.OR.ESCHM.EQ.3) THEN
        PE=1
      ELSEIF(ESCHM.EQ.2) THEN
        PE=0
      ELSE
        CALL XABORT('SNTT3D: UNKNOWN TYPE OF ENERGY CLOSURE RELATION.')
      ENDIF
      IF(MOD(EELEM,2).EQ.1) THEN
        WE(1)=-PE
        WE(2:EELEM+1:2)=1+PE
        IF(EELEM.GE.2) WE(3:EELEM+1:2)=1-PE
      ELSE
        WE(1)=PE
        WE(2:EELEM+1:2)=1-PE
        IF(EELEM.GE.2) WE(3:EELEM+1:2)=1+PE
      ENDIF
      ENDIF
      ! NORMALIZED LEGENDRE POLYNOMIAL CONSTANTS
      DO IEL=1,MAX(IELEM,EELEM)
        CST(IEL)=SQRT(2.0*IEL-1.0)
      ENDDO
*----
*  COMPUTE ISOTROPIC FLUX INDICES.
*----
      NM=IELEM**3*EELEM
      NMX=IELEM**2*EELEM
      NMY=NMX
      NMZ=NMX
      NME=IELEM**3
      LL4=LX*LY*LZ*NSCT*NM
      IF(IGE.EQ.0)THEN
        NUN=LL4+(NMX*LY*LZ+NMY*LX*LZ+NMZ*LX*LY)*NPQ
      ELSEIF(IGE.EQ.2)THEN
        IF((NCODE(5)==1).and.(NCODE(5)==NCODE(6)))THEN
           NUN=LL4
        ELSE
           NUN=LL4+(LX*LY*LZ)*NMZ*NPQ
        ENDIF
      ELSE
         CALL XABORT('SNTT3D: CHECK SPATIAL SCHEME DISCRETISATION '//
     1       'PARAMETER.')
      ENDIF
      DO I=1,LX*LY*LZ
         IDL(I)=(I-1)*NSCT*NM+1
      ENDDO
*----
*  SET BOUNDARY CONDITIONS.
*----
      DO 240 I=1,6
      IF(NCODE(I).NE.1) ZCODE(I)=1.0
      IF(NCODE(I).EQ.5) CALL XABORT('SNTT3D: SYME BC NOT ALLOWED.')
      IF(NCODE(I).EQ.7) CALL XABORT('SNTT3D: ZERO FLUX BC NOT ALLOWED.')
  240 CONTINUE
*----
*  CHECK FOR INVALID VIRTUAL ELEMENTS.
*----
      DO 292 I=2,LX-1
      DO 291 J=2,LY-1
      DO 290 K=2,LZ-1
      IF(MAT(I,J,K).EQ.0) THEN
         L1=(NCODE(1).NE.1)
         DO 251 J1=1,J-1
         DO 250 K1=1,K-1
         L1=L1.OR.(MAT(I,J1,K1).NE.0)
  250    CONTINUE
  251    CONTINUE
         L2=(NCODE(2).NE.1)
         DO 261 J1=J+1,LY
         DO 260 K1=K+1,LZ
         L2=L2.OR.(MAT(I,J1,K1).NE.0)
  260    CONTINUE
  261    CONTINUE
         L3=(NCODE(3).NE.1)
         DO 271 I1=1,I-1
         DO 270 K1=1,K-1
         L3=L3.OR.(MAT(I1,J,K1).NE.0)
  270    CONTINUE
  271    CONTINUE
         L4=(NCODE(4).NE.1)
         DO 281 I1=I+1,LX
         DO 280 K1=K+1,LZ
         L4=L4.OR.(MAT(I1,J,K1).NE.0)
  280    CONTINUE
  281    CONTINUE
         L5=(NCODE(5).NE.1)
         DO 301 I1=1,I-1
         DO 300 J1=1,J-1
         L5=L5.OR.(MAT(I1,J1,K).NE.0)
  300    CONTINUE
  301    CONTINUE
         L6=(NCODE(6).NE.0)
         DO 311 I1=I+1,LX
         DO 310 J1=I+1,LY
         L6=L6.OR.(MAT(I1,J1,K).NE.0)
  310    CONTINUE
  311    CONTINUE
         IF(L1.AND.L2.AND.L3.AND.L4.AND.L5.AND.L6) THEN
            WRITE(HSMG,'(17HSNTT3D: ELEMENT (,I3,1H,,I3,11H) CANNOT BE,
     1      9H VIRTUAL.)') I,J,K
            CALL XABORT(HSMG)
         ENDIF
      ENDIF
  290 CONTINUE
  291 CONTINUE
  292 CONTINUE
*----
*  SCRATCH STORAGE DEALLOCATION
*----
      DEALLOCATE(ZZ,YY,XX)
      RETURN
      END