Initial commit from Polytechnique Montreal

1 files changed, 959 insertions, 0 deletions

diff --git a/Trivac/src/BIVSFH.f b/Trivac/src/BIVSFH.f
new file mode 100755
index 0000000..d439f4a
--- /dev/null
+++ b/Trivac/src/BIVSFH.f
@@ -0,0 +1,959 @@
+*DECK BIVSFH
+      SUBROUTINE BIVSFH (MAXEV,NBLOS,IMPX,ISPLH,IELEM,LXH,MAT,SIDE,
+     1 NCODE,ICODE,ZCODE,LL4,VOL,IDL,IPERT,KN,QFR,IQFR,BFR,MU)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Numbering corresponding to a Thomas-Raviart-Schneider finite element
+* discretization of a 2-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
+* MAXEV   allocated storage for vector MU.
+* NBLOS   number of lozenges per direction, taking into account
+*         mesh-splitting.
+* IMPX    print parameter.
+* 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.
+* MAT     mixture index assigned to each lozenge.
+* SIDE    side of a lozenge.
+* 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.
+* 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.
+* BFR     element-ordered surface fractions.
+* MU      compressed storage mode indices.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER MAXEV,NBLOS,IMPX,ISPLH,IELEM,LXH,MAT(3,ISPLH**2,LXH),
+     1 NCODE(4),ICODE(4),LL4,IDL(3,NBLOS),IPERT(NBLOS),
+     2 KN(NBLOS,4+6*IELEM*(IELEM+1)),IQFR(NBLOS,6),MU(MAXEV)
+      REAL SIDE,ZCODE(4),VOL(3,NBLOS),QFR(NBLOS,6),BFR(NBLOS,6)
+*----
+*  LOCAL VARIABLES
+*----
+      LOGICAL COND,LL1,LL2
+      INTEGER, DIMENSION(:),ALLOCATABLE :: IP,I1,I3,I4,I5
+      INTEGER, DIMENSION(:,:),ALLOCATABLE :: IZGLOB
+      INTEGER, DIMENSION(:,:,:),ALLOCATABLE :: IJP
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(IJP(LXH,ISPLH,ISPLH),IP(MAXEV),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('BIVSFH: 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('BIVSFH: INVALID VALUE OF LXH(2)')
+      ENDIF
+      ALLOCATE(I1(3*LXH),I3(2*LXH),I4(NBLOS),I5(NBLOS))
+      DO 30 I=1,LXH
+      I3(I)=I
+      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)
+      DO 42 I=1,LXH
+      IOFW=I1(I)
+      IOFX=I1(LXH+I)
+      IOFY=I1(2*LXH+I)
+      DO 41 IM=1,ISPLH
+      DO 40 JM=1,ISPLH
+      IZGLOB(IJP(IOFW,IM,JM),1)=I4(I)
+      IZGLOB(IJP(IOFX,IM,JM),2)=I4(I)
+      IZGLOB(IJP(IOFY,IM,JM),3)=I4(I)
+   40 CONTINUE
+   41 CONTINUE
+   42 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,:4+6*IELEM*(IELEM+1))=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
+         DO 90 I=1,LXH
+         IF(IZGLOB(I,1).EQ.0) GO TO 90
+         KN(I4(I),2)=I5(I3(I))+LT4
+   90    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 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
+         II2=ISPLH*(JKEEP*ISPLH+(ISPLH-JM)*KMAX+K0-1)+IM
+         IF(IZGLOB(II1,1).EQ.0) GO TO 120
+         KN(I4(II1),2)=I5(II2)+LT4
+  110    CONTINUE
+  115    CONTINUE
+  120    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
+         DO 140 I=1,LXH
+         IF(IZGLOB(I,1).EQ.0) GO TO 140
+         KN(I4(I),3)=I5(I3(LXH+I))+2*LT4
+  140    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 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
+         II2=ISPLH*(JKEEP*ISPLH+(ISPLH-IM)*KMAX+K0-1)+(ISPLH-JM+1)
+         IF(IZGLOB(II1,1).EQ.0) GO TO 170
+         KN(I4(II1),3)=I5(II2)+2*LT4
+  160    CONTINUE
+  165    CONTINUE
+  170    CONTINUE
+      ENDIF
+*----
+*  COMPUTE THE FLUX PERMUTATION PART OF MATRIX KN (Y <--> W)
+*----
+      IF(ISPLH.EQ.1) THEN
+         DO 180 I=1,LXH
+         IF(IZGLOB(I,1).EQ.0) GO TO 180
+         KN(I4(I),4)=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 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
+         II2=ISPLH*(IKEEP*ISPLH+(JM-1)*JMAX+J00-1)+(ISPLH-IM+1)
+         IF(IZGLOB(II1,1).EQ.0) GO TO 210
+         KN(I4(II1),4)=I4(II2)
+  200    CONTINUE
+  205    CONTINUE
+  210    CONTINUE
+      ENDIF
+      DEALLOCATE(I5,I4,I3,I1)
+*----
+*  SET THE CURRENT NUMBERING PART OF MATRIX KN AND MATRIX QFR (W-AXIS)
+*----
+      LL4W0=(2*NBLOS*IELEM+(2*NBC-1)*ISPLH)*IELEM
+      LL4F=3*LT4*IELEM*IELEM
+      QFR(:NBLOS,:6)=0.0
+      IQFR(:NBLOS,:6)=0
+      BFR(:NBLOS,:6)=0.0
+      ALBEDO=0.5*(1.0-ZCODE(1))/(1.0+ZCODE(1))
+      NELEM=IELEM*(IELEM+1)
+      NB1=(2*NBC*IELEM*ISPLH+1)*IELEM*ISPLH
+      KEL=0
+      NDDIR=LL4F
+      NUM=0
+      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
+      KN(NUM,1)=NUM
+      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
+      DO 250 I=1,IELEM+1
+      LCOUR=LCOUR+1
+      ITEMP = NDDIR
+     >      + (JEL-1)*(2*(NBC+JSTAGE-1)*IELEM*ISPLH+1)*IELEM
+     >      + (IRANG-1)*(2*IELEM*ISPLH)
+     >      + (IEL-1)*IELEM
+     >      + (J-1)*(2*(NBC+JSTAGE-1)*IELEM*ISPLH+1) + I
+      IF(LCOUR.GT.NELEM) CALL XABORT('BIVSFH: bug1')
+      IF(KEL.GT.NBLOS) CALL XABORT('BIVSFH: bug2')
+      KN(NUM,4+LCOUR)=ITEMP
+      KN(NUM,4+NELEM+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
+            KN(NUM,4+(I-1)*(IELEM+1)+1)=0
+  260       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(NUM,1)=SIDE/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(NUM,1)=SIDE
+            IQFR(NUM,1)=ICODE(1)
+         ENDIF
+         IF((NCODE(1).EQ.1).OR.(NCODE(1).EQ.7)) BFR(NUM,1)=SIDE
+      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
+            KN(NUM,4+NELEM+I*(IELEM+1))=0
+  270       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(NUM,2)=SIDE/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(NUM,2)=SIDE
+            IQFR(NUM,2)=ICODE(1)
+         ENDIF
+         IF((NCODE(1).EQ.1).OR.(NCODE(1).EQ.7)) BFR(NUM,2)=SIDE
+      ENDIF
+  280 CONTINUE
+  281 CONTINUE
+  282 CONTINUE
+      NDDIR=NDDIR+NB1+(2*(JSTAGE-1)*IELEM*ISPLH)*IELEM*ISPLH
+  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
+      KN(NUM,1)=NUM
+      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
+      DO 300 I=1,IELEM+1
+      LCOUR=LCOUR+1
+      ITEMP = NDDIR
+     >      + (JEL-1)*(2*(2*NBC-1+NBC-JSTAGE)*IELEM*ISPLH+1)*IELEM
+     >      + (IRANG-1)*(2*IELEM*ISPLH)
+     >      + (IEL-1)*IELEM
+     >      + (J-1)*(2*(2*NBC-1+NBC-JSTAGE)*IELEM*ISPLH+1) + I
+      IF(LCOUR.GT.NELEM) CALL XABORT('BIVSFH: bug3')
+      IF(KEL.GT.NBLOS) CALL XABORT('BIVSFH: bug4')
+      KN(NUM,4+LCOUR)=ITEMP
+      KN(NUM,4+NELEM+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
+            KN(NUM,4+(I-1)*(IELEM+1)+1)=0
+  310       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(NUM,1)=SIDE/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(NUM,1)=SIDE
+            IQFR(NUM,1)=ICODE(1)
+         ENDIF
+         IF((NCODE(1).EQ.1).OR.(NCODE(1).EQ.7)) BFR(NUM,1)=SIDE
+      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
+            KN(NUM,4+NELEM+I*(IELEM+1))=0
+  320       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(NUM,2)=SIDE/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(NUM,2)=SIDE
+            IQFR(NUM,2)=ICODE(1)
+         ENDIF
+         IF((NCODE(1).EQ.1).OR.(NCODE(1).EQ.7)) BFR(NUM,2)=SIDE
+      ENDIF
+  330 CONTINUE
+  331 CONTINUE
+  332 CONTINUE
+      NDDIR=NDDIR+(2*(2*NBC-1)*IELEM*ISPLH+1)*IELEM*ISPLH
+     >           -(2*(JSTAGE-NBC)*IELEM*ISPLH)*IELEM*ISPLH
+  340 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,2)-LT4)=NUM
+  350 CONTINUE
+      KEL=0
+      NUM=0
+      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
+      DO 360 I=1,IELEM+1
+      LCOUR=LCOUR+1
+      ITEMP = NDDIR
+     >      + (JEL-1)*(2*(NBC+JSTAGE-1)*IELEM*ISPLH+1)*IELEM
+     >      + (IRANG-1)*(2*IELEM*ISPLH)
+     >      + (IEL-1)*IELEM
+     >      + (J-1)*(2*(NBC+JSTAGE-1)*IELEM*ISPLH+1) + I
+      IF(LCOUR.GT.NELEM) CALL XABORT('BIVSFH: bug5')
+      IF(KEL.GT.NBLOS) CALL XABORT('BIVSFH: bug6')
+      KN(IP(NUM),4+2*NELEM+LCOUR)=ITEMP
+      KN(IP(NUM),4+3*NELEM+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
+            KN(IP(NUM),4+2*NELEM+(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/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(IP(NUM),3)=SIDE
+            IQFR(NUM,3)=ICODE(1)
+         ENDIF
+         IF((NCODE(1).EQ.1).OR.(NCODE(1).EQ.7)) BFR(NUM,3)=SIDE
+      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
+            KN(IP(NUM),4+3*NELEM+I*(IELEM+1))=0
+  380       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(IP(NUM),4)=SIDE/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(IP(NUM),4)=SIDE
+            IQFR(NUM,4)=ICODE(1)
+         ENDIF
+         IF((NCODE(1).EQ.1).OR.(NCODE(1).EQ.7)) BFR(NUM,4)=SIDE
+      ENDIF
+  390 CONTINUE
+  391 CONTINUE
+  392 CONTINUE
+      NDDIR=NDDIR+NB1+(2*(JSTAGE-1)*IELEM*ISPLH)*IELEM*ISPLH
+  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
+      DO 410 I=1,IELEM+1
+      LCOUR=LCOUR+1
+      ITEMP = NDDIR
+     >      + (JEL-1)*(2*(2*NBC-1+NBC-JSTAGE)*IELEM*ISPLH+1)*IELEM
+     >      + (IRANG-1)*(2*IELEM*ISPLH)
+     >      + (IEL-1)*IELEM
+     >      + (J-1)*(2*(2*NBC-1+NBC-JSTAGE)*IELEM*ISPLH+1) + I
+      IF(LCOUR.GT.NELEM) CALL XABORT('BIVSFH: bug7')
+      IF(KEL.GT.NBLOS) CALL XABORT('BIVSFH: bug8')
+      KN(IP(NUM),4+2*NELEM+LCOUR)=ITEMP
+      KN(IP(NUM),4+3*NELEM+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
+            KN(IP(NUM),4+2*NELEM+(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/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(IP(NUM),3)=SIDE
+            IQFR(NUM,3)=ICODE(1)
+         ENDIF
+         IF((NCODE(1).EQ.1).OR.(NCODE(1).EQ.7)) BFR(NUM,3)=SIDE
+      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
+            KN(IP(NUM),4+3*NELEM+I*(IELEM+1))=0
+  430       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(IP(NUM),4)=SIDE/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(IP(NUM),4)=SIDE
+            IQFR(NUM,4)=ICODE(1)
+         ENDIF
+         IF((NCODE(1).EQ.1).OR.(NCODE(1).EQ.7)) BFR(NUM,4)=SIDE
+      ENDIF
+  440 CONTINUE
+  441 CONTINUE
+  442 CONTINUE
+      NDDIR=NDDIR+(2*(2*NBC-1)*IELEM*ISPLH+1)*IELEM*ISPLH
+     >           -(2*(JSTAGE-NBC)*IELEM*ISPLH)*IELEM*ISPLH
+  450 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,3)-2*LT4)=NUM
+  460 CONTINUE
+      KEL=0
+      NUM=0
+      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
+      DO 470 I=1,IELEM+1
+      LCOUR=LCOUR+1
+      ITEMP = NDDIR
+     >      + (JEL-1)*(2*(NBC+JSTAGE-1)*IELEM*ISPLH+1)*IELEM
+     >      + (IRANG-1)*(2*IELEM*ISPLH)
+     >      + (IEL-1)*IELEM
+     >      + (J-1)*(2*(NBC+JSTAGE-1)*IELEM*ISPLH+1) + I
+      IF(LCOUR.GT.NELEM) CALL XABORT('BIVSFH: bug9')
+      IF(KEL.GT.NBLOS) CALL XABORT('BIVSFH: bug10')
+      KN(IP(NUM),4+4*NELEM+LCOUR)=ITEMP
+      KN(IP(NUM),4+5*NELEM+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
+            KN(IP(NUM),4+4*NELEM+(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/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(IP(NUM),5)=SIDE
+            IQFR(NUM,5)=ICODE(1)
+         ENDIF
+         IF((NCODE(1).EQ.1).OR.(NCODE(1).EQ.7)) BFR(NUM,5)=SIDE
+      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
+            KN(IP(NUM),4+5*NELEM+I*(IELEM+1))=0
+  490       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(IP(NUM),6)=SIDE/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(IP(NUM),6)=SIDE
+            IQFR(NUM,6)=ICODE(1)
+         ENDIF
+         IF((NCODE(1).EQ.1).OR.(NCODE(1).EQ.7)) BFR(NUM,6)=SIDE
+      ENDIF
+  500 CONTINUE
+  501 CONTINUE
+  502 CONTINUE
+      NDDIR=NDDIR+NB1+(2*(JSTAGE-1)*IELEM*ISPLH)*IELEM*ISPLH
+  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
+      DO 520 I=1,IELEM+1
+      LCOUR=LCOUR+1
+      ITEMP = NDDIR
+     >      + (JEL-1)*(2*(2*NBC-1+NBC-JSTAGE)*IELEM*ISPLH+1)*IELEM
+     >      + (IRANG-1)*(2*IELEM*ISPLH)
+     >      + (IEL-1)*IELEM
+     >      + (J-1)*(2*(2*NBC-1+NBC-JSTAGE)*IELEM*ISPLH+1) + I
+      IF(LCOUR.GT.NELEM) CALL XABORT('BIVSFH: bug11')
+      IF(KEL.GT.NBLOS) CALL XABORT('BIVSFH: bug12')
+      KN(IP(NUM),4+4*NELEM+LCOUR)=ITEMP
+      KN(IP(NUM),4+5*NELEM+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
+            KN(IP(NUM),4+4*NELEM+(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/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(IP(NUM),5)=SIDE
+            IQFR(NUM,5)=ICODE(1)
+         ENDIF
+         IF((NCODE(1).EQ.1).OR.(NCODE(1).EQ.7)) BFR(NUM,5)=SIDE
+      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
+            KN(IP(NUM),4+5*NELEM+I*(IELEM+1))=0
+  540       CONTINUE
+         ELSE IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(IP(NUM),6)=SIDE/ALBEDO
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(IP(NUM),6)=SIDE
+            IQFR(NUM,6)=ICODE(1)
+         ENDIF
+         IF((NCODE(1).EQ.1).OR.(NCODE(1).EQ.7)) BFR(NUM,6)=SIDE
+      ENDIF
+  550 CONTINUE
+  551 CONTINUE
+  552 CONTINUE
+      NDDIR=NDDIR+(2*(2*NBC-1)*IELEM*ISPLH+1)*IELEM*ISPLH
+     >           -(2*(JSTAGE-NBC)*IELEM*ISPLH)*IELEM*ISPLH
+  560 CONTINUE
+*----
+*  COMPUTE THE SURFACE FRACTIONS
+*----
+      SURFTOT=0.0
+      DO 566 I=1,NBLOS
+      DO 565 J=1,6
+      SURFTOT=SURFTOT+BFR(I,J)
+  565 CONTINUE
+  566 CONTINUE
+      IF(SURFTOT.GT.0.0) THEN
+         DO 575 I=1,NBLOS
+         DO 570 J=1,6
+         BFR(I,J)=BFR(I,J)/SURFTOT
+  570    CONTINUE
+  575    CONTINUE
+      ENDIF
+*----
+*  REORDER THE UNKNOWNS AND REMOVE THE UNUSED UNKNOWNS INDICES FROM KN
+*----
+      IP(:LL4F+3*LL4W0)=0
+      LL4=0
+      DO 591 KEL=1,LT4
+      DO 582 IFLUX=1,4
+      NUM=KN(KEL,IFLUX)
+      DO 581 K2=1,IELEM
+      DO 580 K1=1,IELEM
+      JND1=(NUM-1)*IELEM**2+(K2-1)*IELEM+K1
+      IF(JND1.GT.MAXEV) CALL XABORT('BIVSFH: MAXEV OVERFLOW(1).')
+      IF(IP(JND1).EQ.0) THEN
+         LL4=LL4+1
+         IP(JND1)=LL4
+      ENDIF
+  580 CONTINUE
+  581 CONTINUE
+  582 CONTINUE
+      DO 590 ICOUR=1,6*NELEM
+      IND=ABS(KN(KEL,4+ICOUR))
+      IF(IND.GT.MAXEV) CALL XABORT('BIVSFH: MAXEV OVERFLOW(2).')
+      IF(IND.NE.0) THEN
+         IF(IP(IND).EQ.0) THEN
+            LL4=LL4+1
+            IP(IND)=LL4
+         ENDIF
+      ENDIF
+  590 CONTINUE
+  591 CONTINUE
+      DO 605 KEL=1,LT4
+      DO 595 IFLUX=1,4
+      NUM=KN(KEL,IFLUX)
+      KN(KEL,IFLUX)=IP((NUM-1)*IELEM**2+1)
+  595 CONTINUE
+      DO 600 ICOUR=1,6*NELEM
+      IF(KN(KEL,4+ICOUR).NE.0) THEN
+         IND=KN(KEL,4+ICOUR)
+         KN(KEL,4+ICOUR)=SIGN(IP(ABS(IND)),IND)
+      ENDIF
+  600 CONTINUE
+  605 CONTINUE
+*----
+*  PRINT A FEW GEOMETRY CHARACTERISTICS
+*----
+      IF(IMPX.GT.0) THEN
+         write(6,*) ' '
+         write(6,*) 'ISPLH =',ISPLH
+         write(6,*) 'IELEM =',IELEM
+         write(6,*) 'NELEM =',NELEM
+         write(6,*) 'NBLOS =',NBLOS
+         write(6,*) 'LL4F  =',LL4F
+         write(6,*) 'LL4   =',LL4
+         write(6,*) 'NBC   =',NBC
+      ENDIF
+*----
+*  SET IPERT
+*----
+      KEL=0
+      DO 613 JSTAGE=1,NBC
+      DO 612 JEL=1,ISPLH
+      DO 611 IRANG=1,NBC+JSTAGE-1
+      DO 610 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
+  610 CONTINUE
+  611 CONTINUE
+  612 CONTINUE
+  613 CONTINUE
+      DO 623 JSTAGE=NBC+1,2*NBC-1
+      DO 622 JEL=1,ISPLH
+      DO 621 IRANG=1,(2*NBC-2)-(JSTAGE-NBC-1)
+      DO 620 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
+  620 CONTINUE
+  621 CONTINUE
+  622 CONTINUE
+  623 CONTINUE
+      IF(KEL.NE.NBLOS) CALL XABORT('BIVSFH: IPERT FAILURE.')
+*----
+*  SET IDL AND VOL
+*----
+      NUM=0
+      IDL(:3,:NBLOS)=0
+      VOL(:3,:NBLOS)=0.0
+      DO 630 KEL=1,NBLOS
+      KEL2=IPERT(KEL)
+      IF(KEL2.EQ.0) GO TO 630
+      NUM=NUM+1
+      IDL(:3,KEL2)=KN(NUM,:3)
+      VOL(:3,KEL2)=2.59807587*SIDE*SIDE/REAL(3)
+  630 CONTINUE
+      IF(IMPX.GT.2) THEN
+         WRITE(6,800) 'MAT',(((MAT(I,J,K),I=1,3),J=1,ISPLH**2),K=1,LXH)
+         WRITE(6,800) 'IDL',((IDL(I,J),I=1,3),J=1,NBLOS)
+         WRITE(6,810) 'VOL',((VOL(I,J),I=1,3),J=1,NBLOS)
+      ENDIF
+*----
+*  COMPUTE THE SYSTEM MATRIX BANDWIDTH.
+*----
+      MU(:LL4)=1
+      NUM=0
+      DO 690 KEL=1,NBLOS
+      IF(IZGLOB(KEL,1).EQ.0) GO TO 690
+      NUM=NUM+1
+      DO 663 K4=0,1
+      DO 662 K3=0,IELEM-1
+      DO 661 K2=1,IELEM+1
+      INW1=ABS(KN(NUM,4+K4*NELEM+K3*(IELEM+1)+K2))
+      INX1=ABS(KN(NUM,4+(K4+2)*NELEM+K3*(IELEM+1)+K2))
+      INY1=ABS(KN(NUM,4+(K4+4)*NELEM+K3*(IELEM+1)+K2))
+      DO 650 K1=1,IELEM+1
+      INW2=ABS(KN(NUM,4+K4*NELEM+K3*(IELEM+1)+K1))
+      INX2=ABS(KN(NUM,4+(K4+2)*NELEM+K3*(IELEM+1)+K1))
+      INY2=ABS(KN(NUM,4+(K4+4)*NELEM+K3*(IELEM+1)+K1))
+      IF((INW2.NE.0).AND.(INW1.NE.0)) THEN
+         MU(INW1)=MAX(MU(INW1),INW1-INW2+1)
+         MU(INW2)=MAX(MU(INW2),INW2-INW1+1)
+      ENDIF
+      IF((INX2.NE.0).AND.(INX1.NE.0)) THEN
+         MU(INX1)=MAX(MU(INX1),INX1-INX2+1)
+         MU(INX2)=MAX(MU(INX2),INX2-INX1+1)
+      ENDIF
+      IF((INY2.NE.0).AND.(INY1.NE.0)) THEN
+         MU(INY1)=MAX(MU(INY1),INY1-INY2+1)
+         MU(INY2)=MAX(MU(INY2),INY2-INY1+1)
+      ENDIF
+  650 CONTINUE
+      DO 660 K1=0,IELEM-1
+      IF(K4.EQ.0) THEN
+         JND1=KN(NUM,1)+K3*IELEM+K1
+         JND2=KN(NUM,2)+K3*IELEM+K1
+         JND3=KN(NUM,3)+K3*IELEM+K1
+      ELSE
+         JND1=KN(NUM,2)+K1*IELEM+K3
+         JND2=KN(NUM,3)+K1*IELEM+K3
+         JND3=KN(NUM,4)+K1*IELEM+K3
+      ENDIF
+      IF(INW1.NE.0) THEN
+         MU(JND1)=MAX(MU(JND1),JND1-INW1+1)
+         MU(INW1)=MAX(MU(INW1),INW1-JND1+1)
+      ENDIF
+      IF(INX1.NE.0) THEN
+         MU(JND2)=MAX(MU(JND2),JND2-INX1+1)
+         MU(INX1)=MAX(MU(INX1),INX1-JND2+1)
+      ENDIF
+      IF(INY1.NE.0) THEN
+         MU(JND3)=MAX(MU(JND3),JND3-INY1+1)
+         MU(INY1)=MAX(MU(INY1),INY1-JND3+1)
+      ENDIF
+  660 CONTINUE
+  661 CONTINUE
+  662 CONTINUE
+  663 CONTINUE
+      ITRS=0
+      DO I=1,LT4
+         IF(KN(I,1).EQ.KN(NUM,4)) THEN
+            ITRS=I
+            GO TO 670
+         ENDIF
+      ENDDO
+      CALL XABORT('BIVSFH: ITRS FAILURE.')
+  670 DO 685 I=1,NELEM
+      INW1=ABS(KN(ITRS,4+I))
+      INX1=ABS(KN(NUM,4+2*NELEM+I))
+      INY1=ABS(KN(NUM,4+4*NELEM+I))
+      DO 680 J=1,NELEM
+      INW2=ABS(KN(NUM,4+NELEM+J))
+      INX2=ABS(KN(NUM,4+3*NELEM+J))
+      INY2=ABS(KN(NUM,4+5*NELEM+J))
+      IF((INY2.NE.0).AND.(INW1.NE.0)) THEN
+         MU(INW1)=MAX(MU(INW1),INW1-INY2+1)
+         MU(INY2)=MAX(MU(INY2),INY2-INW1+1)
+      ENDIF
+      IF((INW2.NE.0).AND.(INX1.NE.0)) THEN
+         MU(INX1)=MAX(MU(INX1),INX1-INW2+1)
+         MU(INW2)=MAX(MU(INW2),INW2-INX1+1)
+      ENDIF
+      IF((INX2.NE.0).AND.(INY1.NE.0)) THEN
+         MU(INY1)=MAX(MU(INY1),INY1-INX2+1)
+         MU(INX2)=MAX(MU(INX2),INX2-INY1+1)
+      ENDIF
+  680 CONTINUE
+  685 CONTINUE
+  690 CONTINUE
+      MUMAX=0
+      IIMAX=0
+      DO 700 I=1,LL4
+      MUMAX=MAX(MUMAX,MU(I))
+      IIMAX=IIMAX+MU(I)
+      MU(I)=IIMAX
+  700 CONTINUE
+*
+      IF(IMPX.GT.0) WRITE(6,820) LL4
+      IF(IMPX.GT.2) THEN
+         WRITE (6,830) MUMAX,IIMAX
+         WRITE (6,840)
+         DO 710 K=1,LXH*ISPLH**2
+         WRITE (6,850) K,(IZGLOB(K,I),I=1,3)
+  710    CONTINUE
+         WRITE (6,860)
+         DO 720 K=1,LT4
+         WRITE (6,870) K,(KN(K,I),I=1,4+2*NELEM)
+         WRITE (6,880) 'X',(KN(K,I),I=4+2*NELEM+1,4+4*NELEM)
+         WRITE (6,880) 'Y',(KN(K,I),I=4+4*NELEM+1,4+6*NELEM)
+  720    CONTINUE
+         WRITE (6,890)
+         DO 730 K=1,LXH*ISPLH**2
+         WRITE (6,900) K,(QFR(K,I),I=1,6)
+  730    CONTINUE
+      ENDIF
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(IZGLOB,IJP,IP)
+      RETURN
+*
+  800 FORMAT(1X,A3/14(2X,I6))
+  810 FORMAT(1X,A3/7(2X,E12.5))
+  820 FORMAT(31H NUMBER OF UNKNOWNS PER GROUP =,I6)
+  830 FORMAT(/41H BIVSFH: MAXIMUM BANDWIDTH FOR MATRICES =,I6/9X,
+     1 51HNUMBER OF TERMS IN THE COMPRESSED SYSTEM MATRICES =,I10)
+  840 FORMAT(/22H NUMBERING OF HEXAGONS/1X,21(1H-)//8H ELEMENT,4X,
+     1 24H W ----- X ----- Y -----)
+  850 FORMAT(1X,I6,5X,3I8)
+  860 FORMAT(/22H NUMBERING OF UNKNOWNS/1X,21(1H-)//8H ELEMENT,5X,
+     1 27H---> W ---> X ---> Y ---> W,4X,8HCURRENTS,89(1H.))
+  870 FORMAT(1X,I6,5X,4I7,4X,1HW,12I8:/(45X,12I8))
+  880 FORMAT(44X,A1,12I8:/(45X,12I8))
+  890 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-))
+  900 FORMAT(1X,I6,5X,1P,10E10.1/(12X,1P,10E10.1))
+      END