Initial commit from Polytechnique Montreal

1 files changed, 1022 insertions, 0 deletions

diff --git a/Trivac/src/TRISFH.f b/Trivac/src/TRISFH.f
new file mode 100755
index 0000000..f9cfd62
--- /dev/null
+++ b/Trivac/src/TRISFH.f
@@ -0,0 +1,1022 @@
+*DECK TRISFH
+      SUBROUTINE TRISFH (IMPX,MAXKN,MAXIP,NBLOS,ISPLH,IELEM,LXH,LZ,MAT,
+     1 SIDE,ZZZ,NCODE,ICODE,ZCODE,LL4,LL4F,LL4W,LL4X,LL4Y,LL4Z,VOL,
+     2 IDL,IPERT,ZZ,FRZ,KN,QFR,IQFR)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Numbering corresponding to a Thomas-Raviart-Schneider finite element
+* discretization of a 3-D hexagonal geometry.
+*
+*Copyright:
+* Copyright (C) 2006 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
+* IMPX    print parameter.
+* MAXKN   number of components in KN.
+* MAXIP   maximum number of currents
+* NBLOS   number of lozenges per direction in 3D with mesh-splitting.
+* ISPLH   mesh-splitting in 3*ISPLH**2 lozenges per hexagon.
+* IELEM   degree of the Lagrangian finite elements: =1 (linear);
+*         =2 (parabolic); =3 (cubic).
+* LXH     number of hexagons in a plane.
+* LZ      number of axial planes.
+* MAT     mixture index assigned to each lozenge.
+* SIDE    side of a lozenge.
+* ZZZ     Z-coordinates of the axial planes.
+* NCODE   type of boundary condition applied on each side (I=1: hbc):
+*         NCODE(I)=1: VOID;          =2: REFL;       =6: ALBE;
+*                 =5: SYME;          =7: ZERO.
+* ICODE   physical albedo index on each side of the domain.
+* ZCODE   albedo corresponding to boundary condition 'VOID' on each
+*         side (ZCODE(I)=0.0 by default).
+*
+*Parameters: output
+* LL4     order of the system matrices.
+* LL4F    number of flux unknowns.
+* LL4W    number of W-directed currents
+* LL4X    number of X-directed currents
+* LL4Y    number of Y-directed currents
+* LL4Z    number of Z-directed currents
+* ZZ      Z-sides of each hexagon.
+* FRZ     volume fractions for the axial SYME boundary condition.
+* VOL     volume of each lozenge.
+* IDL     position of the average flux component associated with each
+*         lozenge.
+* IPERT   mixture permutation index.
+* KN      ADI permutation indices for the volumes and currents.
+* QFR     element-ordered boundary conditions.
+* IQFR    element-ordered physical albedo indices.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER IMPX,MAXKN,MAXIP,NBLOS,ISPLH,IELEM,LXH,LZ,
+     1 MAT(3,ISPLH**2,LXH*LZ),NCODE(6),ICODE(6),LL4,LL4F,LL4W,LL4X,
+     2 LL4Y,LL4Z,IDL(3,NBLOS),IPERT(NBLOS),KN(NBLOS,MAXKN/NBLOS),
+     3 IQFR(NBLOS,8)
+      REAL SIDE,ZZZ(LZ+1),ZCODE(6),VOL(3,NBLOS),ZZ(3,NBLOS),
+     1 FRZ(NBLOS),QFR(NBLOS,8)
+*----
+*  LOCAL VARIABLES
+*----
+      LOGICAL COND,LL1,LL2
+      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: IJP
+      INTEGER, DIMENSION(:,:), ALLOCATABLE :: IZGLOB
+      INTEGER, DIMENSION(:), ALLOCATABLE :: IP,I1,I3,I4,I5
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(IJP(LXH,ISPLH,ISPLH),IP(MAXIP),IZGLOB(NBLOS,3))
+*----
+*  THOMAS-RAVIART-SCHNEIDER SPECIFIC NUMEROTATION
+*----
+      NBC=INT((SQRT(REAL((4*LXH-1)/3))+1.)/2.)
+      IF(LXH.NE.1+3*NBC*(NBC-1)) CALL XABORT('TRISFH: INVALID VALUE OF '
+     1 //'LXH(1).')
+      IF(ISPLH.EQ.1) THEN
+         DO 10 I=1,LXH
+         IJP(I,1,1)=I
+   10    CONTINUE
+      ELSE
+         I=0
+         DO 23 I0=1,2*NBC-1
+         JMAX=NBC+I0-1
+         IF(I0.GE.NBC) JMAX=3*NBC-I0-1
+         IKEEP=I
+         DO 22 J0=1,JMAX
+         I=I+1
+         DO 21 IM=1,ISPLH
+         DO 20 JM=1,ISPLH
+         IJP(I,IM,JM)=ISPLH*(IKEEP*ISPLH+(IM-1)*JMAX+J0-1)+JM
+   20    CONTINUE
+   21    CONTINUE
+   22    CONTINUE
+   23    CONTINUE
+         IF(I.NE.LXH) CALL XABORT('TRISFH: INVALID VALUE OF LXH(2)')
+      ENDIF
+      ALLOCATE(I1(3*LXH),I3(2*LXH),I4(NBLOS),I5(NBLOS))
+      DO 25 I=1,LXH
+      I3(I)=I
+   25 CONTINUE
+      DO 30 I=1,LXH*LZ
+      I4(I)=0
+      IF(MAT(1,1,I).GT.0) I4(I)=I
+   30 CONTINUE
+      IZGLOB(:NBLOS,:3)=0
+      J1=2+3*(NBC-1)*(NBC-2)
+      IF(NBC.EQ.1) J1=1
+      J3=J1+2*NBC-2
+      J5=J3+2*NBC-2
+      CALL BIVPER(J1,1,LXH,LXH,I1(1),I3)
+      CALL BIVPER(J3,3,LXH,LXH,I1(LXH+1),I3)
+      CALL BIVPER(J5,5,LXH,LXH,I1(2*LXH+1),I3)
+      I=0
+      DO 43 IZ=1,LZ
+      DO 42 IX=1,LXH
+      I=I+1
+      IOFW=I1(IX)
+      IOFX=I1(LXH+IX)
+      IOFY=I1(2*LXH+IX)
+      DO 41 IM=1,ISPLH
+      DO 40 JM=1,ISPLH
+      IZGLOB((IZ-1)*LXH*ISPLH**2+IJP(IOFW,IM,JM),1)=I4(I)
+      IZGLOB((IZ-1)*LXH*ISPLH**2+IJP(IOFX,IM,JM),2)=I4(I)
+      IZGLOB((IZ-1)*LXH*ISPLH**2+IJP(IOFY,IM,JM),3)=I4(I)
+   40 CONTINUE
+   41 CONTINUE
+   42 CONTINUE
+   43 CONTINUE
+      DO 50 I=1,LXH
+      II1=I1(I)
+      II2=I1(LXH+I)
+      II3=I1(2*LXH+I)
+      I3(II1)=II2
+      I3(LXH+II1)=II3
+   50 CONTINUE
+*----
+*  COMPUTE THE FLUX PERMUTATION PART OF MATRIX KN (W <--> X)
+*----
+      KN(:NBLOS,:3+6*IELEM*IELEM*(IELEM+2))=0
+      LT4=0
+      DO 70 II2=1,NBLOS
+      I=IZGLOB(II2,1)
+      I4(II2)=0
+      IF(I.NE.0) THEN
+         LT4=LT4+1
+         I4(II2)=LT4
+      ENDIF
+   70 CONTINUE
+      LT4=0
+      DO 80 II2=1,NBLOS
+      I=IZGLOB(II2,2)
+      I5(II2)=0
+      IF(I.NE.0) THEN
+         LT4=LT4+1
+         I5(II2)=LT4
+      ENDIF
+   80 CONTINUE
+      IF(ISPLH.EQ.1) THEN
+         I=0
+         DO 95 IZ=1,LZ
+         DO 90 IX=1,LXH
+         I=I+1
+         IF(IZGLOB(I,1).EQ.0) GO TO 90
+         IOF=(IZ-1)*LXH+I3(IX)
+         KN(I4(I),1)=I5(IOF)+LT4
+   90    CONTINUE
+   95    CONTINUE
+      ELSE
+         I=0
+         DO 105 I0=1,2*NBC-1
+         JMAX=NBC+I0-1
+         IF(I0.GE.NBC) JMAX=3*NBC-I0-1
+         IKEEP=I
+         DO 100 J0=1,JMAX
+         I=I+1
+         I1(I)=JMAX
+         I1(LXH+I)=IKEEP
+         I1(2*LXH+I)=J0
+  100    CONTINUE
+  105    CONTINUE
+         DO 125 IZ=1,LZ
+         DO 120 I=1,LXH
+         JMAX=I1(I)
+         IKEEP=I1(LXH+I)
+         J00=I1(2*LXH+I)
+         KMAX=I1(I3(I))
+         JKEEP=I1(LXH+I3(I))
+         K0=I1(2*LXH+I3(I))
+         DO 115 IM=1,ISPLH
+         DO 110 JM=1,ISPLH
+         II1=ISPLH*(IKEEP*ISPLH+(IM-1)*JMAX+J00-1)+JM
+         IOF1=(IZ-1)*LXH*ISPLH**2+II1
+         IF(IZGLOB(IOF1,1).EQ.0) GO TO 120
+         II2=ISPLH*(JKEEP*ISPLH+(ISPLH-JM)*KMAX+K0-1)+IM
+         IOF2=(IZ-1)*LXH*ISPLH**2+II2
+         KN(I4(IOF1),1)=I5(IOF2)+LT4
+  110    CONTINUE
+  115    CONTINUE
+  120    CONTINUE
+  125    CONTINUE
+      ENDIF
+*----
+*  COMPUTE THE FLUX PERMUTATION PART OF MATRIX KN (X <--> Y)
+*----
+      LT4=0
+      DO 130 II2=1,NBLOS
+      I=IZGLOB(II2,3)
+      I5(II2)=0
+      IF(I.NE.0) THEN
+         LT4=LT4+1
+         I5(II2)=LT4
+      ENDIF
+  130 CONTINUE
+      IF(ISPLH.EQ.1) THEN
+         I=0
+         DO 145 IZ=1,LZ
+         DO 140 IX=1,LXH
+         I=I+1
+         IF(IZGLOB(I,1).EQ.0) GO TO 140
+         IOF=(IZ-1)*LXH+I3(LXH+IX)
+         KN(I4(I),2)=I5(IOF)+2*LT4
+  140    CONTINUE
+  145    CONTINUE
+      ELSE
+         I=0
+         DO 155 I0=1,2*NBC-1
+         JMAX=NBC+I0-1
+         IF(I0.GE.NBC) JMAX=3*NBC-I0-1
+         IKEEP=I
+         DO 150 J0=1,JMAX
+         I=I+1
+         I1(I)=JMAX
+         I1(LXH+I)=IKEEP
+         I1(2*LXH+I)=J0
+  150    CONTINUE
+  155    CONTINUE
+         DO 175 IZ=1,LZ
+         DO 170 I=1,LXH
+         JMAX=I1(I)
+         IKEEP=I1(LXH+I)
+         J00=I1(2*LXH+I)
+         KMAX=I1(I3(LXH+I))
+         JKEEP=I1(LXH+I3(LXH+I))
+         K0=I1(2*LXH+I3(LXH+I))
+         DO 165 IM=1,ISPLH
+         DO 160 JM=1,ISPLH
+         II1=ISPLH*(IKEEP*ISPLH+(IM-1)*JMAX+J00-1)+JM
+         IOF1=(IZ-1)*LXH*ISPLH**2+II1
+         IF(IZGLOB(IOF1,1).EQ.0) GO TO 170
+         II2=ISPLH*(JKEEP*ISPLH+(ISPLH-IM)*KMAX+K0-1)+(ISPLH-JM+1)
+         IOF2=(IZ-1)*LXH*ISPLH**2+II2
+         KN(I4(IOF1),2)=I5(IOF2)+2*LT4
+  160    CONTINUE
+  165    CONTINUE
+  170    CONTINUE
+  175    CONTINUE
+      ENDIF
+*----
+*  COMPUTE THE FLUX PERMUTATION PART OF MATRIX KN (Y <--> W)
+*----
+      IF(ISPLH.EQ.1) THEN
+         DO 180 I=1,LXH*LZ
+         IF(IZGLOB(I,1).EQ.0) GO TO 180
+         KN(I4(I),3)=I4(I)
+  180    CONTINUE
+      ELSE
+         I=0
+         DO 195 I0=1,2*NBC-1
+         JMAX=NBC+I0-1
+         IF(I0.GE.NBC) JMAX=3*NBC-I0-1
+         IKEEP=I
+         DO 190 J0=1,JMAX
+         I=I+1
+         I1(I)=JMAX
+         I1(LXH+I)=IKEEP
+         I1(2*LXH+I)=J0
+  190    CONTINUE
+  195    CONTINUE
+         DO 215 IZ=1,LZ
+         DO 210 I=1,LXH
+         JMAX=I1(I)
+         IKEEP=I1(LXH+I)
+         J00=I1(2*LXH+I)
+         DO 205 IM=1,ISPLH
+         DO 200 JM=1,ISPLH
+         II1=ISPLH*(IKEEP*ISPLH+(IM-1)*JMAX+J00-1)+JM
+         IOF1=(IZ-1)*LXH*ISPLH**2+II1
+         IF(IZGLOB(IOF1,1).EQ.0) GO TO 210
+         II2=ISPLH*(IKEEP*ISPLH+(JM-1)*JMAX+J00-1)+(ISPLH-IM+1)
+         IOF2=(IZ-1)*LXH*ISPLH**2+II2
+         KN(I4(IOF1),3)=I4(IOF2)
+  200    CONTINUE
+  205    CONTINUE
+  210    CONTINUE
+  215    CONTINUE
+      ENDIF
+      DEALLOCATE(I5,I4,I3,I1)
+*----
+*  SET THE CURRENT NUMBERING PART OF MATRIX KN AND MATRIX QFR (W-AXIS)
+*----
+      LL4W0=(2*LXH*ISPLH*IELEM+2*NBC-1)*ISPLH*LZ*IELEM**2
+      LL4Z0=3*LXH*(LZ+1)*(ISPLH**2)*IELEM**2
+      LL4F=3*LT4*IELEM**3
+      QFR(:NBLOS,:8)=0.0
+      IQFR(:NBLOS,:8)=0
+      ALBEDO=0.5*(1.0-ZCODE(1))/(1.0+ZCODE(1))
+      NELEH=(IELEM+1)*IELEM**2
+      NELEZ=6*IELEM**2
+      NB1=2*NBC*ISPLH*IELEM+1
+      NB2=2*(2*NBC-1)*ISPLH*IELEM+1
+      KEL=0
+      NDDIR=0
+      NUM=0
+      DO 345 IZ=1,LZ
+      FRACT=1.0
+      IF((NCODE(5).EQ.5).AND.(IZ.EQ.1)) FRACT=0.5
+      IF((NCODE(6).EQ.5).AND.(IZ.EQ.LZ)) FRACT=0.5
+      DZZ=ZZZ(IZ+1)-ZZZ(IZ)
+      DO 290 JSTAGE=1,NBC
+      DO 282 JEL=1,ISPLH
+      DO 281 IRANG=1,NBC+JSTAGE-1
+      DO 280 IEL=1,ISPLH
+      KEL=KEL+1
+      IF(IZGLOB(KEL,1).EQ.0) GO TO 280
+      NUM=NUM+1
+      IF((IRANG.EQ.1).AND.(IEL.EQ.1)) THEN
+         LL1=.TRUE.
+      ELSE
+         LL1=(IZGLOB(KEL-1,1).EQ.0)
+      ENDIF
+      IF((IRANG.EQ.NBC+JSTAGE-1).AND.(IEL.EQ.ISPLH)) THEN
+         LL2=.TRUE.
+      ELSE
+         LL2=(IZGLOB(KEL+1,1).EQ.0)
+      ENDIF
+      LCOUR=0
+      DO 255 J=1,IELEM**2
+      DO 250 I=1,IELEM+1
+      LCOUR=LCOUR+1
+      ITEMP = NDDIR
+     >      + (JEL-1)*(2*(NBC+JSTAGE-1)*IELEM*ISPLH+1)*IELEM**2
+     >      + (IRANG-1)*(2*IELEM*ISPLH)
+     >      + (IEL-1)*IELEM
+     >      + (J-1)*(2*(NBC+JSTAGE-1)*IELEM*ISPLH+1) + I
+      IF(LCOUR.GT.NELEH) CALL XABORT('TRISFH: BUG1')
+      IF(KEL.GT.NBLOS) CALL XABORT('TRISFH: BUG2')
+      KN(NUM,3+LCOUR)=ITEMP
+      KN(NUM,3+NELEH+LCOUR)=ITEMP+IELEM*ISPLH
+  250 CONTINUE
+  255 CONTINUE
+      IF(LL1) THEN
+         COND=(NCODE(1).EQ.2).OR.((NCODE(1).EQ.1).AND.(ZCODE(1).EQ.1.0))
+         IF(COND) THEN
+            DO 260 I=1,IELEM**2
+            KN(NUM,3+(I-1)*(IELEM+1)+1)=0
+  260       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(NUM,1)=SIDE*DZZ*FRACT/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(NUM,1)=SIDE*DZZ*FRACT
+            IQFR(NUM,1)=ICODE(1)
+         ENDIF
+      ENDIF
+      IF(LL2) THEN
+         COND=(NCODE(1).EQ.2).OR.((NCODE(1).EQ.1).AND.(ZCODE(1).EQ.1.0))
+         IF(COND) THEN
+            DO 270 I=1,IELEM**2
+            KN(NUM,3+NELEH+I*(IELEM+1))=0
+  270       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(NUM,2)=SIDE*DZZ*FRACT/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(NUM,2)=SIDE*DZZ*FRACT
+            IQFR(NUM,1)=ICODE(1)
+         ENDIF
+      ENDIF
+  280 CONTINUE
+  281 CONTINUE
+  282 CONTINUE
+      NDDIR=NDDIR+(NB1+2*(JSTAGE-1)*ISPLH*IELEM)*ISPLH*IELEM**2
+  290 CONTINUE
+*
+      DO 340 JSTAGE=NBC+1,2*NBC-1
+      DO 332 JEL=1,ISPLH
+      DO 331 IRANG=1,(2*NBC-2)-(JSTAGE-NBC-1)
+      DO 330 IEL=1,ISPLH
+      KEL=KEL+1
+      IF(IZGLOB(KEL,1).EQ.0) GO TO 330
+      NUM=NUM+1
+      IF((IRANG.EQ.1).AND.(IEL.EQ.1)) THEN
+         LL1=.TRUE.
+      ELSE
+         LL1=(IZGLOB(KEL-1,1).EQ.0)
+      ENDIF
+      IF((IRANG.EQ.(2*NBC-2)-(JSTAGE-NBC-1)).AND.(IEL.EQ.ISPLH)) THEN
+         LL2=.TRUE.
+      ELSE
+         LL2=(IZGLOB(KEL+1,1).EQ.0)
+      ENDIF
+      LCOUR=0
+      DO 305 J=1,IELEM**2
+      DO 300 I=1,IELEM+1
+      LCOUR=LCOUR+1
+      ITEMP = NDDIR
+     >      + (JEL-1)*(2*(2*NBC-1+NBC-JSTAGE)*IELEM*ISPLH+1)*IELEM**2
+     >      + (IRANG-1)*(2*IELEM*ISPLH)
+     >      + (IEL-1)*IELEM
+     >      + (J-1)*(2*(2*NBC-1+NBC-JSTAGE)*IELEM*ISPLH+1) + I
+      IF(LCOUR.GT.NELEH) CALL XABORT('TRISFH: BUG3')
+      IF(KEL.GT.NBLOS) CALL XABORT('TRISFH: BUG4')
+      KN(NUM,3+LCOUR)=ITEMP
+      KN(NUM,3+NELEH+LCOUR)=ITEMP+IELEM*ISPLH
+  300 CONTINUE
+  305 CONTINUE
+      IF(LL1) THEN
+         COND=(NCODE(1).EQ.2).OR.((NCODE(1).EQ.1).AND.(ZCODE(1).EQ.1.0))
+         IF(COND) THEN
+            DO 310 I=1,IELEM**2
+            KN(NUM,3+(I-1)*(IELEM+1)+1)=0
+  310       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(NUM,1)=SIDE*DZZ*FRACT/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(NUM,1)=SIDE*DZZ*FRACT
+            IQFR(NUM,1)=ICODE(1)
+         ENDIF
+      ENDIF
+      IF(LL2) THEN
+         COND=(NCODE(1).EQ.2).OR.((NCODE(1).EQ.1).AND.(ZCODE(1).EQ.1.0))
+         IF(COND) THEN
+            DO 320 I=1,IELEM**2
+            KN(NUM,3+NELEH+I*(IELEM+1))=0
+  320       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(NUM,2)=SIDE*DZZ*FRACT/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(NUM,2)=SIDE*DZZ*FRACT
+            IQFR(NUM,2)=ICODE(1)
+         ENDIF
+      ENDIF
+  330 CONTINUE
+  331 CONTINUE
+  332 CONTINUE
+      NDDIR=NDDIR+(NB2-2*(JSTAGE-NBC)*ISPLH*IELEM)*ISPLH*IELEM**2
+  340 CONTINUE
+  345 CONTINUE
+*----
+*  SET THE CURRENT NUMBERING PART OF MATRIX KN AND MATRIX QFR (X-AXIS)
+*----
+      IP(:NBLOS)=0
+      DO 350 NUM=1,LT4
+      IP(KN(NUM,1)-LT4)=NUM
+  350 CONTINUE
+      KEL=0
+      NUM=0
+      DO 455 IZ=1,LZ
+      FRACT=1.0
+      IF((NCODE(5).EQ.5).AND.(IZ.EQ.1)) FRACT=0.5
+      IF((NCODE(6).EQ.5).AND.(IZ.EQ.LZ)) FRACT=0.5
+      DZZ=ZZZ(IZ+1)-ZZZ(IZ)
+      DO 400 JSTAGE=1,NBC
+      DO 392 JEL=1,ISPLH
+      DO 391 IRANG=1,NBC+JSTAGE-1
+      DO 390 IEL=1,ISPLH
+      KEL=KEL+1
+      IF(IZGLOB(KEL,2).EQ.0) GO TO 390
+      NUM=NUM+1
+      IF((IRANG.EQ.1).AND.(IEL.EQ.1)) THEN
+         LL1=.TRUE.
+      ELSE
+         LL1=(IZGLOB(KEL-1,2).EQ.0)
+      ENDIF
+      IF((IRANG.EQ.NBC+JSTAGE-1).AND.(IEL.EQ.ISPLH)) THEN
+         LL2=.TRUE.
+      ELSE
+         LL2=(IZGLOB(KEL+1,2).EQ.0)
+      ENDIF
+      LCOUR=0
+      DO 365 J=1,IELEM**2
+      DO 360 I=1,IELEM+1
+      LCOUR=LCOUR+1
+      ITEMP = NDDIR
+     >      + (JEL-1)*(2*(NBC+JSTAGE-1)*IELEM*ISPLH+1)*IELEM**2
+     >      + (IRANG-1)*(2*IELEM*ISPLH)
+     >      + (IEL-1)*IELEM
+     >      + (J-1)*(2*(NBC+JSTAGE-1)*IELEM*ISPLH+1) + I
+      IF(LCOUR.GT.NELEH) CALL XABORT('TRISFH: BUG5')
+      IF(KEL.GT.NBLOS) CALL XABORT('TRISFH: BUG6')
+      KN(IP(NUM),3+2*NELEH+LCOUR)=ITEMP
+      KN(IP(NUM),3+3*NELEH+LCOUR)=ITEMP+IELEM*ISPLH
+  360 CONTINUE
+  365 CONTINUE
+      IF(LL1) THEN
+         COND=(NCODE(1).EQ.2).OR.((NCODE(1).EQ.1).AND.(ZCODE(1).EQ.1.0))
+         IF(COND) THEN
+            DO 370 I=1,IELEM**2
+            KN(IP(NUM),3+2*NELEH+(I-1)*(IELEM+1)+1)=0
+  370       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(IP(NUM),3)=SIDE*DZZ*FRACT/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(IP(NUM),3)=SIDE*DZZ*FRACT
+            IQFR(IP(NUM),3)=ICODE(1)
+         ENDIF
+      ENDIF
+      IF(LL2) THEN
+         COND=(NCODE(1).EQ.2).OR.((NCODE(1).EQ.1).AND.(ZCODE(1).EQ.1.0))
+         IF(COND) THEN
+            DO 380 I=1,IELEM**2
+            KN(IP(NUM),3+3*NELEH+I*(IELEM+1))=0
+  380       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(IP(NUM),4)=SIDE*DZZ*FRACT/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(IP(NUM),4)=SIDE*DZZ*FRACT
+            IQFR(IP(NUM),4)=ICODE(1)
+         ENDIF
+      ENDIF
+  390 CONTINUE
+  391 CONTINUE
+  392 CONTINUE
+      NDDIR=NDDIR+(NB1+2*(JSTAGE-1)*ISPLH*IELEM)*ISPLH*IELEM**2
+  400 CONTINUE
+*
+      DO 450 JSTAGE=NBC+1,2*NBC-1
+      DO 442 JEL=1,ISPLH
+      DO 441 IRANG=1,(2*NBC-2)-(JSTAGE-NBC-1)
+      DO 440 IEL=1,ISPLH
+      KEL=KEL+1
+      IF(IZGLOB(KEL,2).EQ.0) GO TO 440
+      NUM=NUM+1
+      IF((IRANG.EQ.1).AND.(IEL.EQ.1)) THEN
+         LL1=.TRUE.
+      ELSE
+         LL1=(IZGLOB(KEL-1,2).EQ.0)
+      ENDIF
+      IF((IRANG.EQ.(2*NBC-2)-(JSTAGE-NBC-1)).AND.(IEL.EQ.ISPLH)) THEN
+         LL2=.TRUE.
+      ELSE
+         LL2=(IZGLOB(KEL+1,2).EQ.0)
+      ENDIF
+      LCOUR=0
+      DO 415 J=1,IELEM**2
+      DO 410 I=1,IELEM+1
+      LCOUR=LCOUR+1
+      ITEMP = NDDIR
+     >      + (JEL-1)*(2*(2*NBC-1+NBC-JSTAGE)*IELEM*ISPLH+1)*IELEM**2
+     >      + (IRANG-1)*(2*IELEM*ISPLH)
+     >      + (IEL-1)*IELEM
+     >      + (J-1)*(2*(2*NBC-1+NBC-JSTAGE)*IELEM*ISPLH+1) + I
+      IF(LCOUR.GT.NELEH) CALL XABORT('TRISFH: BUG7')
+      IF(KEL.GT.NBLOS) CALL XABORT('TRISFH: BUG8')
+      KN(IP(NUM),3+2*NELEH+LCOUR)=ITEMP
+      KN(IP(NUM),3+3*NELEH+LCOUR)=ITEMP+IELEM*ISPLH
+  410 CONTINUE
+  415 CONTINUE
+      IF(LL1) THEN
+         COND=(NCODE(1).EQ.2).OR.((NCODE(1).EQ.1).AND.(ZCODE(1).EQ.1.0))
+         IF(COND) THEN
+            DO 420 I=1,IELEM**2
+            KN(IP(NUM),3+2*NELEH+(I-1)*(IELEM+1)+1)=0
+  420       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(IP(NUM),3)=SIDE*DZZ*FRACT/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(IP(NUM),3)=SIDE*DZZ*FRACT
+            IQFR(IP(NUM),3)=ICODE(1)
+         ENDIF
+      ENDIF
+      IF(LL2) THEN
+         COND=(NCODE(1).EQ.2).OR.((NCODE(1).EQ.1).AND.(ZCODE(1).EQ.1.0))
+         IF(COND) THEN
+            DO 430 I=1,IELEM**2
+            KN(IP(NUM),3+3*NELEH+I*(IELEM+1))=0
+  430       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(IP(NUM),4)=SIDE*DZZ*FRACT/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(IP(NUM),4)=SIDE*DZZ*FRACT
+            IQFR(IP(NUM),4)=ICODE(1)
+         ENDIF
+      ENDIF
+  440 CONTINUE
+  441 CONTINUE
+  442 CONTINUE
+      NDDIR=NDDIR+(NB2-2*(JSTAGE-NBC)*ISPLH*IELEM)*ISPLH*IELEM**2
+  450 CONTINUE
+  455 CONTINUE
+*----
+*  SET THE CURRENT NUMBERING PART OF MATRIX KN AND MATRIX QFR (Y-AXIS)
+*----
+      IP(:NBLOS)=0
+      DO 460 NUM=1,LT4
+      IP(KN(NUM,2)-2*LT4)=NUM
+  460 CONTINUE
+      KEL=0
+      NUM=0
+      DO 565 IZ=1,LZ
+      FRACT=1.0
+      IF((NCODE(5).EQ.5).AND.(IZ.EQ.1)) FRACT=0.5
+      IF((NCODE(6).EQ.5).AND.(IZ.EQ.LZ)) FRACT=0.5
+      DZZ=ZZZ(IZ+1)-ZZZ(IZ)
+      DO 510 JSTAGE=1,NBC
+      DO 502 JEL=1,ISPLH
+      DO 501 IRANG=1,NBC+JSTAGE-1
+      DO 500 IEL=1,ISPLH
+      KEL=KEL+1
+      IF(IZGLOB(KEL,3).EQ.0) GO TO 500
+      NUM=NUM+1
+      IF((IRANG.EQ.1).AND.(IEL.EQ.1)) THEN
+         LL1=.TRUE.
+      ELSE
+         LL1=(IZGLOB(KEL-1,3).EQ.0)
+      ENDIF
+      IF((IRANG.EQ.NBC+JSTAGE-1).AND.(IEL.EQ.ISPLH)) THEN
+         LL2=.TRUE.
+      ELSE
+         LL2=(IZGLOB(KEL+1,3).EQ.0)
+      ENDIF
+      LCOUR=0
+      DO 475 J=1,IELEM**2
+      DO 470 I=1,IELEM+1
+      LCOUR=LCOUR+1
+      ITEMP = NDDIR
+     >      + (JEL-1)*(2*(NBC+JSTAGE-1)*IELEM*ISPLH+1)*IELEM**2
+     >      + (IRANG-1)*(2*IELEM*ISPLH)
+     >      + (IEL-1)*IELEM
+     >      + (J-1)*(2*(NBC+JSTAGE-1)*IELEM*ISPLH+1) + I
+      IF(LCOUR.GT.NELEH) CALL XABORT('TRISFH: BUG9')
+      IF(KEL.GT.NBLOS) CALL XABORT('TRISFH: BUG10')
+      KN(IP(NUM),3+4*NELEH+LCOUR)=ITEMP
+      KN(IP(NUM),3+5*NELEH+LCOUR)=ITEMP+IELEM*ISPLH
+  470 CONTINUE
+  475 CONTINUE
+      IF(LL1) THEN
+         COND=(NCODE(1).EQ.2).OR.((NCODE(1).EQ.1).AND.(ZCODE(1).EQ.1.0))
+         IF(COND) THEN
+            DO 480 I=1,IELEM**2
+            KN(IP(NUM),3+4*NELEH+(I-1)*(IELEM+1)+1)=0
+  480       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(IP(NUM),5)=SIDE*DZZ*FRACT/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(IP(NUM),5)=SIDE*DZZ*FRACT
+            IQFR(IP(NUM),5)=ICODE(1)
+         ENDIF
+      ENDIF
+      IF(LL2) THEN
+         COND=(NCODE(1).EQ.2).OR.((NCODE(1).EQ.1).AND.(ZCODE(1).EQ.1.0))
+         IF(COND) THEN
+            DO 490 I=1,IELEM**2
+            KN(IP(NUM),3+5*NELEH+I*(IELEM+1))=0
+  490       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(IP(NUM),6)=SIDE*DZZ*FRACT/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(IP(NUM),6)=SIDE*DZZ*FRACT
+            IQFR(IP(NUM),6)=ICODE(1)
+         ENDIF
+      ENDIF
+  500 CONTINUE
+  501 CONTINUE
+  502 CONTINUE
+      NDDIR=NDDIR+(NB1+2*(JSTAGE-1)*ISPLH*IELEM)*ISPLH*IELEM**2
+  510 CONTINUE
+*
+      DO 560 JSTAGE=NBC+1,2*NBC-1
+      DO 552 JEL=1,ISPLH
+      DO 551 IRANG=1,(2*NBC-2)-(JSTAGE-NBC-1)
+      DO 550 IEL=1,ISPLH
+      KEL=KEL+1
+      IF(IZGLOB(KEL,3).EQ.0) GO TO 550
+      NUM=NUM+1
+      IF((IRANG.EQ.1).AND.(IEL.EQ.1)) THEN
+         LL1=.TRUE.
+      ELSE
+         LL1=(IZGLOB(KEL-1,3).EQ.0)
+      ENDIF
+      IF((IRANG.EQ.(2*NBC-2)-(JSTAGE-NBC-1)).AND.(IEL.EQ.ISPLH)) THEN
+         LL2=.TRUE.
+      ELSE
+         LL2=(IZGLOB(KEL+1,3).EQ.0)
+      ENDIF
+      LCOUR=0
+      DO 525 J=1,IELEM**2
+      DO 520 I=1,IELEM+1
+      LCOUR=LCOUR+1
+      ITEMP = NDDIR
+     >      + (JEL-1)*(2*(2*NBC-1+NBC-JSTAGE)*IELEM*ISPLH+1)*IELEM**2
+     >      + (IRANG-1)*(2*IELEM*ISPLH)
+     >      + (IEL-1)*IELEM
+     >      + (J-1)*(2*(2*NBC-1+NBC-JSTAGE)*IELEM*ISPLH+1) + I
+      IF(LCOUR.GT.NELEH) CALL XABORT('TRISFH: BUG11')
+      IF(KEL.GT.NBLOS) CALL XABORT('TRISFH: BUG12')
+      KN(IP(NUM),3+4*NELEH+LCOUR)=ITEMP
+      KN(IP(NUM),3+5*NELEH+LCOUR)=ITEMP+IELEM*ISPLH
+  520 CONTINUE
+  525 CONTINUE
+      IF(LL1) THEN
+         COND=(NCODE(1).EQ.2).OR.((NCODE(1).EQ.1).AND.(ZCODE(1).EQ.1.0))
+         IF(COND) THEN
+            DO 530 I=1,IELEM**2
+            KN(IP(NUM),3+4*NELEH+(I-1)*(IELEM+1)+1)=0
+  530       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(IP(NUM),5)=SIDE*DZZ*FRACT/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(IP(NUM),5)=SIDE*DZZ*FRACT
+            IQFR(IP(NUM),5)=ICODE(1)
+         ENDIF
+      ENDIF
+      IF(LL2) THEN
+         COND=(NCODE(1).EQ.2).OR.((NCODE(1).EQ.1).AND.(ZCODE(1).EQ.1.0))
+         IF(COND) THEN
+            DO 540 I=1,IELEM**2
+            KN(IP(NUM),3+5*NELEH+I*(IELEM+1))=0
+  540       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(IP(NUM),6)=SIDE*DZZ*FRACT/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(IP(NUM),6)=SIDE*DZZ*FRACT
+            IQFR(IP(NUM),6)=ICODE(1)
+         ENDIF
+      ENDIF
+  550 CONTINUE
+  551 CONTINUE
+  552 CONTINUE
+      NDDIR=NDDIR+(NB2-2*(JSTAGE-NBC)*ISPLH*IELEM)*ISPLH*IELEM**2
+  560 CONTINUE
+  565 CONTINUE
+*----
+*  SET THE CURRENT NUMBERING PART OF MATRIX KN AND MATRIX QFR (Z-AXIS)
+*----
+      KEL=0
+      NUM=0
+      DO 635 IZ=1,LZ
+      DO 630 IX=1,LXH*ISPLH**2
+      KEL=KEL+1
+      IF(IZGLOB(KEL,1).EQ.0) GO TO 630
+      NUM=NUM+1
+      IF(IZ.EQ.1) THEN
+         LL1=.TRUE.
+      ELSE
+         LL1=(IZGLOB(KEL-LXH*ISPLH**2,1).EQ.0)
+      ENDIF
+      IF(IZ.EQ.LZ) THEN
+         LL2=.TRUE.
+      ELSE
+         LL2=(IZGLOB(KEL+LXH*ISPLH**2,1).EQ.0)
+      ENDIF
+      DO 572 K=0,2 ! THREE LOZENGES PER HEXAGON
+      DO 571 I=0,1 ! FACE ZINF/ZSUP
+      DO 570 J=1,IELEM**2
+      LCOUR=(2*K+I)*IELEM**2+J
+      IF(LCOUR.GT.NELEZ) CALL XABORT('TRISFH: BUG11')
+      IF(KEL.GT.NBLOS) CALL XABORT('TRISFH: BUG12')
+      ITEMP = NDDIR
+     >      + 3*(IX-1)*(LZ+1)*IELEM**2
+     >      + K*(LZ+1)*IELEM**2
+     >      + (J-1)*(LZ+1) + IZ + I
+      KN(NUM,3+6*NELEH+LCOUR)=ITEMP
+  570 CONTINUE
+  571 CONTINUE
+  572 CONTINUE
+*
+*     REFL OR ALBE BOUNDARY CONDITION
+      IF(LL1) THEN
+         COND=(NCODE(5).EQ.2).OR.((NCODE(5).EQ.1).AND.(ZCODE(5).EQ.1.0))
+         IF(COND) THEN
+            DO 585 K=0,2
+            DO 580 J=1,IELEM**2
+            LCINF=2*K*IELEM**2+J
+            KN(NUM,3+6*NELEH+LCINF)=0
+  580       CONTINUE
+  585       CONTINUE
+         ELSE IF((NCODE(5).EQ.1).AND.(ICODE(5).EQ.0)) THEN
+            ALBEDO=0.5*(1.0-ZCODE(5))/(1.0+ZCODE(5))
+            QFR(NUM,7)=0.8660254038*SIDE*SIDE/ALBEDO
+         ELSE IF(NCODE(5).EQ.1) THEN
+            QFR(NUM,7)=0.8660254038*SIDE*SIDE
+            IQFR(NUM,7)=ICODE(5)
+         ENDIF
+      ENDIF
+      IF(LL2) THEN
+         COND=(NCODE(6).EQ.2).OR.((NCODE(6).EQ.1).AND.(ZCODE(6).EQ.1.0))
+         IF(COND) THEN
+            DO 595 K=0,2
+            DO 590 J=1,IELEM**2
+            LCSUP=(2*K+1)*IELEM**2+J
+            KN(NUM,3+6*NELEH+LCSUP)=0
+  590       CONTINUE
+  595       CONTINUE
+         ELSE IF((NCODE(6).EQ.1).AND.(ICODE(6).EQ.0)) THEN
+            ALBEDO=0.5*(1.0-ZCODE(6))/(1.0+ZCODE(6))
+            QFR(NUM,8)=0.8660254038*SIDE*SIDE/ALBEDO
+         ELSE IF(NCODE(6).EQ.1) THEN
+            QFR(NUM,8)=0.8660254038*SIDE*SIDE
+            IQFR(NUM,8)=ICODE(6)
+         ENDIF
+      ENDIF
+*     TRAN BOUNDARY CONDITION
+      IF((IZ.EQ.LZ).AND.(NCODE(6).EQ.4)) THEN
+         DO 605 K=0,2
+         DO 600 J=1,IELEM**2
+         LCSUP=(2*K+1)*IELEM**2+J
+         KN(NUM,3+6*NELEH+LCSUP)=KN(NUM,3+6*NELEH+LCSUP)-LZ
+  600    CONTINUE
+  605    CONTINUE
+      ENDIF
+*     SYME BOUNDARY CONDITION
+      IF((NCODE(5).EQ.5).AND.(IZ.EQ.1)) THEN
+         QFR(NUM,7)=QFR(NUM,8)
+         IQFR(NUM,7)=IQFR(NUM,8)
+         DO 615 K=0,2
+         DO 610 J=1,IELEM**2
+         LCINF=2*K*IELEM**2+J
+         LCSUP=(2*K+1)*IELEM**2+J
+         KN(NUM,3+6*NELEH+LCINF)=-KN(NUM,3+6*NELEH+LCSUP)
+  610    CONTINUE
+  615    CONTINUE
+      ELSE IF((NCODE(6).EQ.5).AND.(IZ.EQ.LZ)) THEN
+         QFR(NUM,8)=QFR(NUM,7)
+         IQFR(NUM,8)=IQFR(NUM,7)
+         DO 625 K=0,2
+         DO 620 J=1,IELEM**2
+         LCINF=2*K*IELEM**2+J
+         LCSUP=(2*K+1)*IELEM**2+J
+         KN(NUM,3+6*NELEH+LCSUP)=-KN(NUM,3+6*NELEH+LCINF)
+  620    CONTINUE
+  625    CONTINUE
+      ENDIF
+  630 CONTINUE
+  635 CONTINUE
+*----
+*  REMOVING THE UNUSED UNKNOWNS INDICES FROM KN
+*----
+      IP(:3*LL4W0+LL4Z0)=0
+      DO 645 KEL=1,LT4
+      DO 640 ICOUR=1,6*NELEH+NELEZ
+      IND=ABS(KN(KEL,3+ICOUR))
+      IF(IND.GT.MAXIP) CALL XABORT('TRISFH: MAXIP OVERFLOW.')
+      IF(IND.NE.0) IP(IND)=1
+  640 CONTINUE
+  645 CONTINUE
+      LL4W=0
+      DO 650 IND=1,LL4W0
+      IF(IP(IND).EQ.1) THEN
+         LL4W=LL4W+1
+         IP(IND)=LL4W
+      ENDIF
+  650 CONTINUE
+      LL4X=0
+      DO 660 IND=1,LL4W0
+      IF(IP(LL4W0+IND).EQ.1) THEN
+         LL4X=LL4X+1
+         IP(LL4W0+IND)=LL4W+LL4X
+      ENDIF
+  660 CONTINUE
+      LL4Y=0
+      DO 670 IND=1,LL4W0
+      IF(IP(2*LL4W0+IND).EQ.1) THEN
+         LL4Y=LL4Y+1
+         IP(2*LL4W0+IND)=LL4W+LL4X+LL4Y
+      ENDIF
+  670 CONTINUE
+      LL4Z=0
+      DO 680 IND=1,LL4Z0
+      IF(IP(3*LL4W0+IND).EQ.1) THEN
+         LL4Z=LL4Z+1
+         IP(3*LL4W0+IND)=LL4W+LL4X+LL4Y+LL4Z
+      ENDIF
+  680 CONTINUE
+      DO 695 KEL=1,LT4
+      DO 690 ICOUR=1,6*NELEH+NELEZ
+      IF(KN(KEL,3+ICOUR).NE.0) THEN
+         IND=KN(KEL,3+ICOUR)
+         KN(KEL,3+ICOUR)=SIGN(IP(ABS(IND)),IND)
+      ENDIF
+  690 CONTINUE
+  695 CONTINUE
+      LL4=LL4F+LL4W+LL4X+LL4Y+LL4Z
+*----
+*  PRINT A FEW GEOMETRY CHARACTERISTICS
+*----
+      IF(IMPX.GT.0) THEN
+         write(6,*) ' '
+         write(6,*) 'ISPLH =',ISPLH
+         write(6,*) 'IELEM =',IELEM
+         write(6,*) 'NELEH =',NELEH
+         write(6,*) 'NELEZ =',NELEZ
+         write(6,*) 'NBLOS =',NBLOS
+         write(6,*) 'LL4F  =',LL4F
+         write(6,*) 'LL4W  =',LL4W
+         write(6,*) 'LL4X  =',LL4X
+         write(6,*) 'LL4Y  =',LL4Y
+         write(6,*) 'LL4Z  =',LL4Z
+         write(6,*) 'NBC   =',NBC
+      ENDIF
+*----
+*  SET IPERT
+*----
+      KEL=0
+      DO 714 IZ=1,LZ
+      DO 703 JSTAGE=1,NBC
+      DO 702 JEL=1,ISPLH
+      DO 701 IRANG=1,NBC+JSTAGE-1
+      DO 700 IEL=1,ISPLH
+      KEL=KEL+1
+      IHEX=IZGLOB(KEL,1)
+      IF(IHEX.EQ.0) THEN
+        IPERT(KEL)=0
+      ELSE
+        IPERT(KEL)=(IHEX-1)*ISPLH**2+(IEL-1)*ISPLH+JEL
+      ENDIF
+      IF(IPERT(KEL).GT.NBLOS) call XABORT('TRISFH: NBLOS OVERFLOW(1)')
+  700 CONTINUE
+  701 CONTINUE
+  702 CONTINUE
+  703 CONTINUE
+      DO 713 JSTAGE=NBC+1,2*NBC-1
+      DO 712 JEL=1,ISPLH
+      DO 711 IRANG=1,(2*NBC-2)-(JSTAGE-NBC-1)
+      DO 710 IEL=1,ISPLH
+      KEL=KEL+1
+      IHEX=IZGLOB(KEL,1)
+      IF(IHEX.EQ.0) THEN
+        IPERT(KEL)=0
+      ELSE
+        IPERT(KEL)=(IHEX-1)*ISPLH**2+(IEL-1)*ISPLH+JEL
+      ENDIF
+      IF(IPERT(KEL).GT.NBLOS) call XABORT('TRISFH: NBLOS OVERFLOW(2)')
+  710 CONTINUE
+  711 CONTINUE
+  712 CONTINUE
+  713 CONTINUE
+  714 CONTINUE
+      IF(KEL.NE.NBLOS) CALL XABORT('TRISFH: IPERT FAILURE.')
+*----
+*  SET IDL, VOL, FRZ AND ZZ
+*----
+      NUM=0
+      IDL(:3,:NBLOS)=0
+      VOL(:3,:NBLOS)=0.0
+      FRZ(:NBLOS)=0.0
+      ZZ(:3,:NBLOS)=0.0
+      DO 725 IZ=1,LZ
+      FRACT=1.0
+      IF((NCODE(5).EQ.5).AND.(IZ.EQ.1)) FRACT=0.5
+      IF((NCODE(6).EQ.5).AND.(IZ.EQ.LZ)) FRACT=0.5
+      DZ=ZZZ(IZ+1)-ZZZ(IZ)
+      DO 720 J=1,LXH*ISPLH**2
+      KEL=(IZ-1)*LXH*ISPLH**2+J
+      KEL2=IPERT(KEL)
+      IF(KEL2.EQ.0) GO TO 720
+      NUM=NUM+1
+      IDL(1,KEL2)=(NUM-1)*IELEM**3+1
+      IDL(2,KEL2)=(KN(NUM,1)-1)*IELEM**3+1
+      IDL(3,KEL2)=(KN(NUM,2)-1)*IELEM**3+1
+      VOL(:3,KEL2)=2.59807587*SIDE*SIDE*DZ*FRACT/REAL(3)
+      FRZ(KEL)=FRACT
+      ZZ(:3,KEL2)=DZ
+  720 CONTINUE
+  725 CONTINUE
+      IF(IMPX.GT.2) THEN
+         WRITE(6,790) 'MAT',(((MAT(I,J,K),I=1,3),J=1,ISPLH**2),
+     1   K=1,LXH*LZ)
+         WRITE(6,790) 'IDL',((IDL(I,J),I=1,3),J=1,NBLOS)
+         WRITE(6,800) 'ZZ ',((ZZ(I,J),I=1,3),J=1,NBLOS)
+         WRITE(6,800) 'VOL',((VOL(I,J),I=1,3),J=1,NBLOS)
+      ENDIF
+*
+      IF(IMPX.GT.0) WRITE(6,810) LL4
+      IF(IMPX.GT.2) THEN
+         WRITE (6,830)
+         DO 730 K=1,NBLOS
+         WRITE (6,840) K,(IZGLOB(K,I),I=1,3)
+  730    CONTINUE
+         WRITE (6,850)
+         DO 740 K=1,LT4
+         WRITE (6,860) K,(KN(K,I),I=1,3+2*NELEH)
+         WRITE (6,870) 'X',(KN(K,I),I=3+2*NELEH+1,3+4*NELEH)
+         WRITE (6,870) 'Y',(KN(K,I),I=3+4*NELEH+1,3+6*NELEH)
+         IF(LL4Z.GT.0) THEN
+            WRITE (6,870) 'Z',(KN(K,I),I=3+6*NELEH+1,3+6*NELEH+NELEZ)
+         ENDIF
+  740    CONTINUE
+         WRITE (6,880)
+         DO 750 K=1,LT4
+         WRITE (6,890) K,(QFR(K,I),I=1,8)
+  750    CONTINUE
+      ENDIF
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(IZGLOB,IP,IJP)
+      RETURN
+*
+  790 FORMAT(1X,A3/14(2X,I6))
+  800 FORMAT(1X,A3/7(2X,E12.5))
+  810 FORMAT(31H NUMBER OF UNKNOWNS PER GROUP =,I8)
+  830 FORMAT(/22H NUMBERING OF HEXAGONS/1X,21(1H-)//8H ELEMENT,4X,
+     1 24H W ----- X ----- Y -----)
+  840 FORMAT(1X,I6,5X,3I8)
+  850 FORMAT(/22H NUMBERING OF UNKNOWNS/1X,21(1H-)//8H ELEMENT,5X,
+     1 20H---> X ---> Y ---> W,4X,8HCURRENTS,89(1H.))
+  860 FORMAT(1X,I6,5X,3I7,4X,1HW,12I8:/(38X,12I8))
+  870 FORMAT(37X,A1,12I8:/(38X,12I8))
+  880 FORMAT(/8H ELEMENT,3X,23HVOID BOUNDARY CONDITION/15X,7(1H-),
+     1 3H W ,7(1H-),3X,7(1H-),3H X ,7(1H-),3X,7(1H-),3H Y ,7(1H-),
+     2 3X,7(1H-),3H Z ,7(1H-))
+  890 FORMAT(1X,I6,5X,1P,10E10.1/(12X,1P,10E10.1))
+      END