path: root/Trivac/src/TRIPKN.f

*DECK TRIPKN
      SUBROUTINE TRIPKN (IELEM,LX,LY,LZ,L4,CYLIND,XXX,YYY,ZZZ,XX,YY,ZZ,
     1 DD,KN,QFR,IQFR,VOL,MAT,NCODE,ICODE,ZCODE,IMPX)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Numbering corresponding to a primal formulation of the finite element
* discretization in a 3-D geometry.
*
*Copyright:
* Copyright (C) 2002 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.
* LX      number of elements along the X axis.
* LY      number of elements along the Y axis.
* LZ      number of elements along the Z axis.
* CYLIND  cylindrical geometry flag (set with CYLIND=.true.).
* IELEM   degree of the Lagrangian finite elements: =1 (linear);
*         =2 (parabolic); =3 (cubic); =4 (quartic).
* 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+;
*         NCODE(I)=1: VOID;  NCODE(I)=2: REFL;  NCODE(I)=4: TRAN;
*         NCODE(I)=5: SYME;  NCODE(I)=7: ZERO;  NCODE(I)=20: CYLI.
* 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).
* 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.
*
*Parameters: output
* L4      total number of unknown (variational coefficients) per
*         energy group (order of system matrices).
* VOL     volume of each element.
* XX      X-directed mesh spacings.
* YY      Y-directed mesh spacings.
* ZZ      Z-directed mesh spacings.
* DD      used with cylindrical geometry.
* KN      element-ordered unknown list.
* QFR     element-ordered boundary conditions.
* IQFR    element-ordered physical albedo indices.
*
*-----------------------------------------------------------------------
*
*----
*  SUBROUTINE ARGUMENTS
*----
      INTEGER IELEM,LX,LY,LZ,L4,KN(LX*LY*LZ*(IELEM+1)**3),MAT(LX*LY*LZ),
     1 IQFR(6*LX*LY*LZ),NCODE(6),ICODE(6),IMPX
      REAL XXX(LX+1),YYY(LY+1),ZZZ(LZ+1),XX(LX*LY*LZ),YY(LX*LY*LZ),
     1 ZZ(LX*LY*LZ),DD(LX*LY*LZ),QFR(6*LX*LY*LZ),VOL(LX*LY*LZ),ZCODE(6)
      LOGICAL CYLIND
*----
*  LOCAL VARIABLES
*----
      LOGICAL LL1,LL2
      INTEGER, DIMENSION(:), ALLOCATABLE :: IP,IWRK
*
      ALB(X)=0.5*(1.0-X)/(1.0+X)
*----
*  SCRATCH STORAGE ALLOCATION
*----
      ALLOCATE(IP((1+IELEM*LX)*(1+IELEM*LY)*(1+IELEM*LZ)),
     1 IWRK((1+IELEM*LX)*(1+IELEM*LY)*(1+IELEM*LZ)))
*
      IF(IMPX.GT.0) WRITE(6,500) LX,LY,LZ
      MAXIP=(1+IELEM*LX)*(1+IELEM*LY)*(1+IELEM*LZ)
      LC=1+IELEM
      LL=LC*LC*LC
*----
*  IDENTIFICATION OF THE GEOMETRY. MAIN LOOP OVER THE ELEMENTS
*----
      IX=LX*(LC-1)+1
      IXY=(LY*(LC-1)+1)*IX
      IXYZ=(LZ*(LC-1)+1)*IXY
      NUM1=0
      NUM2=0
      KEL=0
      DO 182 K0=1,LZ
      DO 181 K1=1,LY
      DO 180 K2=1,LX
      KEL=KEL+1
      XX(KEL)=0.0
      YY(KEL)=0.0
      ZZ(KEL)=0.0
      VOL(KEL)=0.0
      IF(MAT(KEL).LE.0) GO TO 180
      XX(KEL)=XXX(K2+1)-XXX(K2)
      YY(KEL)=YYY(K1+1)-YYY(K1)
      ZZ(KEL)=ZZZ(K0+1)-ZZZ(K0)
      IF(CYLIND) DD(KEL)=0.5*(XXX(K2)+XXX(K2+1))
      IND1=(LC-1)*((K0-1)*IXY+(K1-1)*IX+(K2-1))
      L=0
      DO 12 I=1,LC
      DO 11 J=1,LC
      DO 10 K=1,LC
      L=L+1
      KN(NUM1+L)=IND1+(I-1)*IXY+(J-1)*IX+K
   10 CONTINUE
   11 CONTINUE
   12 CONTINUE
      DO 20 IC=1,6
      QFR(NUM2+IC)=0.0
      IQFR(NUM2+IC)=0
   20 CONTINUE
      KK1=KEL-1
      KK2=KEL+1
      KK3=KEL-LX
      KK4=KEL+LX
      KK5=KEL-LX*LY
      KK6=KEL+LX*LY
      FRX=1.0
      FRY=1.0
      FRZ=1.0
*----
*  VOID, REFL OR ZERO BOUNDARY CONTITION
*----
      IF(K2.EQ.1) THEN
         LL1=.TRUE.
      ELSE
         LL1=(MAT(KK1).EQ.0)
      ENDIF
      IF(LL1) THEN
         IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
            QFR(NUM2+1)=ALB(ZCODE(1))
         ELSE IF(NCODE(1).EQ.1) THEN
            QFR(NUM2+1)=1.0
            IQFR(NUM2+1)=ICODE(1)
         ELSE IF(NCODE(1).EQ.7) THEN
            L=0
            DO 32 I=1,LC
            DO 31 J=1,LC
            DO 30 K=1,LC
            L=L+1
            IF(K.EQ.1) KN(NUM1+L)=0
   30       CONTINUE
   31       CONTINUE
   32       CONTINUE
         ENDIF
      ENDIF
*
      IF(K2.EQ.LX) THEN
         LL1=.TRUE.
      ELSE
         LL1=(MAT(KK2).EQ.0)
      ENDIF
      IF(LL1) THEN
         IF((NCODE(2).EQ.1).AND.(ICODE(2).EQ.0)) THEN
            QFR(NUM2+2)=ALB(ZCODE(2))
         ELSE IF(NCODE(2).EQ.1) THEN
            QFR(NUM2+2)=1.0
            IQFR(NUM2+2)=ICODE(2)
         ELSE IF(NCODE(2).EQ.7) THEN
            L=0
            DO 42 I=1,LC
            DO 41 J=1,LC
            DO 40 K=1,LC
            L=L+1
            IF(K.EQ.LC) KN(NUM1+L)=0
   40       CONTINUE
   41       CONTINUE
   42       CONTINUE
         ENDIF
      ENDIF
*
      IF(K1.EQ.1) THEN
         LL1=.TRUE.
      ELSE
         LL1=(MAT(KK3).EQ.0)
      ENDIF
      IF(LL1) THEN
         IF((NCODE(3).EQ.1).AND.(ICODE(3).EQ.0)) THEN
            QFR(NUM2+3)=ALB(ZCODE(3))
         ELSE IF(NCODE(3).EQ.1) THEN
            QFR(NUM2+3)=1.0
            IQFR(NUM2+3)=ICODE(3)
         ELSE IF(NCODE(3).EQ.7) THEN
            L=0
            DO 52 I=1,LC
            DO 51 J=1,LC
            DO 50 K=1,LC
            L=L+1
            IF(J.EQ.1) KN(NUM1+L)=0
   50       CONTINUE
   51       CONTINUE
   52       CONTINUE
         ENDIF
      ENDIF
*
      IF(K1.EQ.LY) THEN
         LL1=.TRUE.
      ELSE
         LL1=(MAT(KK4).EQ.0)
      ENDIF
      IF(LL1) THEN
         IF((NCODE(4).EQ.1).AND.(ICODE(4).EQ.0)) THEN
            QFR(NUM2+4)=ALB(ZCODE(4))
         ELSE IF(NCODE(4).EQ.1) THEN
            QFR(NUM2+4)=1.0
            IQFR(NUM2+4)=ICODE(4)
         ELSE IF(NCODE(4).EQ.7) THEN
            L=0
            DO 62 I=1,LC
            DO 61 J=1,LC
            DO 60 K=1,LC
            L=L+1
            IF(J.EQ.LC) KN(NUM1+L)=0
   60       CONTINUE
   61       CONTINUE
   62       CONTINUE
         ENDIF
      ENDIF
*
      IF(K0.EQ.1) THEN
         LL1=.TRUE.
      ELSE
         LL1=(MAT(KK5).EQ.0)
      ENDIF
      IF(LL1) THEN
         IF((NCODE(5).EQ.1).AND.(ICODE(5).EQ.0)) THEN
            QFR(NUM2+5)=ALB(ZCODE(5))
         ELSE IF(NCODE(5).EQ.1) THEN
            QFR(NUM2+5)=1.0
            IQFR(NUM2+5)=ICODE(5)
         ELSE IF(NCODE(5).EQ.7) THEN
            L=0
            DO 72 I=1,LC
            DO 71 J=1,LC
            DO 70 K=1,LC
            L=L+1
            IF(I.EQ.1) KN(NUM1+L)=0
   70       CONTINUE
   71       CONTINUE
   72       CONTINUE
         ENDIF
      ENDIF
*
      IF(K0.EQ.LZ) THEN
         LL1=.TRUE.
      ELSE
         LL1=(MAT(KK6).EQ.0)
      ENDIF
      IF(LL1) THEN
         IF((NCODE(6).EQ.1).AND.(ICODE(6).EQ.0)) THEN
            QFR(NUM2+6)=ALB(ZCODE(6))
         ELSE IF(NCODE(6).EQ.1) THEN
            QFR(NUM2+6)=1.0
            IQFR(NUM2+6)=ICODE(6)
         ELSE IF(NCODE(6).EQ.7) THEN
            L=0
            DO 82 I=1,LC
            DO 81 J=1,LC
            DO 80 K=1,LC
            L=L+1
            IF(I.EQ.LC) KN(NUM1+L)=0
   80       CONTINUE
   81       CONTINUE
   82       CONTINUE
         ENDIF
      ENDIF
*----
*  TRAN BOUNDARY CONDITION
*----
      IF((K2.EQ.LX).AND.(NCODE(2).EQ.4)) THEN
         DO 91 I=1,LC
         DO 90 J=1,LC
         M=(I-1)*LC*LC+(J-1)*LC+LC
         KN(NUM1+M)=KN(NUM1+M)-IX+1
   90    CONTINUE
   91    CONTINUE
      ENDIF
      IF((K1.EQ.LY).AND.(NCODE(4).EQ.4)) THEN
         DO 101 I=1,LC
         DO 100 K=1,LC
         M=(I-1)*LC*LC+(LC-1)*LC+K
         KN(NUM1+M)=KN(NUM1+M)-IXY+IX
  100    CONTINUE
  101    CONTINUE
      ENDIF
      IF((K0.EQ.LZ).AND.(NCODE(6).EQ.4)) THEN
         DO 111 J=1,LC
         DO 110 K=1,LC
         M=(LC-1)*LC*LC+(J-1)*LC+K
         KN(NUM1+M)=KN(NUM1+M)-IXYZ+IXY
  110    CONTINUE
  111    CONTINUE
      ENDIF
*----
*  SYME BOUNDARY CONDITION
*----
      IF((NCODE(1).EQ.5).AND.(K2.EQ.1)) THEN
         QFR(NUM2+1)=QFR(NUM2+2)
         IQFR(NUM2+1)=IQFR(NUM2+2)
         FRX=0.5
         DO 122 I=1,LC
         DO 121 J=1,LC
         DO 120 K=1,(LC+1)/2
         L=(I-1)*LC*LC+(J-1)*LC+K
         M=(I-1)*LC*LC+(J-1)*LC+(LC-K+1)
         KN(NUM1+L)=KN(NUM1+M)
  120    CONTINUE
  121    CONTINUE
  122    CONTINUE
      ELSE IF((NCODE(2).EQ.5).AND.(K2.EQ.LX)) THEN
         QFR(NUM2+2)=QFR(NUM2+1)
         IQFR(NUM2+2)=IQFR(NUM2+1)
         FRX=0.5
         DO 132 I=1,LC
         DO 131 J=1,LC
         DO 130 K=(LC+2)/2,LC
         L=(I-1)*LC*LC+(J-1)*LC+K
         M=(I-1)*LC*LC+(J-1)*LC+(LC-K+1)
         KN(NUM1+L)=KN(NUM1+M)
  130    CONTINUE
  131    CONTINUE
  132    CONTINUE
      ENDIF
      IF((NCODE(3).EQ.5).AND.(K1.EQ.1)) THEN
         QFR(NUM2+3)=QFR(NUM2+4)
         IQFR(NUM2+3)=IQFR(NUM2+4)
         FRY=0.5
         DO 142 I=1,LC
         DO 141 J=1,(LC+1)/2
         DO 140 K=1,LC
         L=(I-1)*LC*LC+(J-1)*LC+K
         M=(I-1)*LC*LC+(LC-J)*LC+K
         KN(NUM1+L)=KN(NUM1+M)
  140    CONTINUE
  141    CONTINUE
  142    CONTINUE
      ELSE IF((NCODE(4).EQ.5).AND.(K1.EQ.LY)) THEN
         QFR(NUM2+4)=QFR(NUM2+3)
         IQFR(NUM2+4)=IQFR(NUM2+3)
         FRY=0.5
         DO 152 I=1,LC
         DO 151 J=(LC+2)/2,LC
         DO 150 K=1,LC
         L=(I-1)*LC*LC+(J-1)*LC+K
         M=(I-1)*LC*LC+(LC-J)*LC+K
         KN(NUM1+L)=KN(NUM1+M)
  150    CONTINUE
  151    CONTINUE
  152    CONTINUE
      ENDIF
      IF((NCODE(5).EQ.5).AND.(K0.EQ.1)) THEN
         QFR(NUM2+5)=QFR(NUM2+6)
         IQFR(NUM2+5)=IQFR(NUM2+6)
         FRZ=0.5
         DO 162 I=1,(LC+1)/2
         DO 161 J=1,LC
         DO 160 K=1,LC
         L=(I-1)*LC*LC+(J-1)*LC+K
         M=(LC-I)*LC*LC+(J-1)*LC+K
         KN(NUM1+L)=KN(NUM1+M)
  160    CONTINUE
  161    CONTINUE
  162    CONTINUE
      ELSE IF((NCODE(6).EQ.5).AND.(K0.EQ.LZ)) THEN
         QFR(NUM2+6)=QFR(NUM2+5)
         IQFR(NUM2+6)=IQFR(NUM2+5)
         FRZ=0.5
         DO 172 I=(LC+2)/2,LC
         DO 171 J=1,LC
         DO 170 K=1,LC
         L=(I-1)*LC*LC+(J-1)*LC+K
         M=(LC-I)*LC*LC+(J-1)*LC+K
         KN(NUM1+L)=KN(NUM1+M)
  170    CONTINUE
  171    CONTINUE
  172    CONTINUE
      ENDIF
*
      VOL0=XX(KEL)*YY(KEL)*ZZ(KEL)*FRX*FRY*FRZ
      IF(CYLIND) VOL0=6.2831853072*DD(KEL)*VOL0
      VOL(KEL)=VOL0
      QFR(NUM2+1)=QFR(NUM2+1)*VOL0/XX(KEL)
      QFR(NUM2+2)=QFR(NUM2+2)*VOL0/XX(KEL)
      QFR(NUM2+3)=QFR(NUM2+3)*VOL0/YY(KEL)
      QFR(NUM2+4)=QFR(NUM2+4)*VOL0/YY(KEL)
      QFR(NUM2+5)=QFR(NUM2+5)*VOL0/ZZ(KEL)
      QFR(NUM2+6)=QFR(NUM2+6)*VOL0/ZZ(KEL)
      NUM1=NUM1+LL
      NUM2=NUM2+6
  180 CONTINUE
  181 CONTINUE
  182 CONTINUE
* END OF THE MAIN LOOP OVER ELEMENTS.
*
*----
*  PROCESSING OF 1-D AND 2-D CASES
*----
      LL1=(LX.EQ.1).AND.(NCODE(1).EQ.2).AND.(NCODE(2).EQ.5)
     1 .AND.(IELEM.GT.1)
      LL2=(LX.EQ.1).AND.(NCODE(1).EQ.5).AND.(NCODE(2).EQ.2)
     1 .AND.(IELEM.GT.1)
      IF(LL1.OR.LL2) THEN
         NUM1=0
         DO 200 KEL=1,LX*LY*LZ
         IF(MAT(KEL).EQ.0) GO TO 200
         DO 192 I=1,LC
         DO 191 J=1,LC
         DO 190 K=2,LC
         L=(I-1)*LC*LC+(J-1)*LC+K
         M=(I-1)*LC*LC+(J-1)*LC+1
         KN(NUM1+L)=KN(NUM1+M)
  190    CONTINUE
  191    CONTINUE
  192    CONTINUE
         NUM1=NUM1+LL
  200    CONTINUE
      ENDIF
      LL1=(LY.EQ.1).AND.(NCODE(3).EQ.2).AND.(NCODE(4).EQ.5)
     1 .AND.(IELEM.GT.1)
      LL2=(LY.EQ.1).AND.(NCODE(3).EQ.5).AND.(NCODE(4).EQ.2)
     1 .AND.(IELEM.GT.1)
      IF(LL1.OR.LL2) THEN
         NUM1=0
         DO 220 KEL=1,LX*LY*LZ
         IF(MAT(KEL).EQ.0) GO TO 220
         DO 212 I=1,LC
         DO 211 J=2,LC
         DO 210 K=1,LC
         L=(I-1)*LC*LC+(J-1)*LC+K
         M=(I-1)*LC*LC+K
         KN(NUM1+L)=KN(NUM1+M)
  210    CONTINUE
  211    CONTINUE
  212    CONTINUE
         NUM1=NUM1+LL
  220    CONTINUE
      ENDIF
      LL1=(LZ.EQ.1).AND.(NCODE(5).EQ.2).AND.(NCODE(6).EQ.5)
     1 .AND.(IELEM.GT.1)
      LL2=(LZ.EQ.1).AND.(NCODE(5).EQ.5).AND.(NCODE(6).EQ.2)
     1 .AND.(IELEM.GT.1)
      IF(LL1.OR.LL2) THEN
         NUM1=0
         DO 240 KEL=1,LX*LY*LZ
         IF(MAT(KEL).EQ.0) GO TO 240
         DO 232 I=2,LC
         DO 231 J=1,LC
         DO 230 K=1,LC
         L=(I-1)*LC*LC+(J-1)*LC+K
         M=(J-1)*LC+K
         KN(NUM1+L)=KN(NUM1+M)
  230    CONTINUE
  231    CONTINUE
  232    CONTINUE
         NUM1=NUM1+LL
  240    CONTINUE
      ENDIF
*----
*  JUXTAPOSITION OF A CHECKERBOARD OVER THE REACTOR DOMAIN
*----
      LZTOT=LZ*(LC-1)+1
      LYTOT=LY*(LC-1)+1
      LXTOT=LX*(LC-1)+1
      DO 250 I=1,LXTOT*LYTOT*LZTOT
      IWRK(I)=-1
  250 CONTINUE
      NUM1=0
      KEL=0
      DO 272 K0=1,LZ
      LK0=(K0-1)*(LC-1)
      DO 271 K1=1,LY
      LK1=(K1-1)*(LC-1)
      DO 270 K2=1,LX
      KEL=KEL+1
      IF(MAT(KEL).EQ.0) GO TO 270
      LK2=(K2-1)*(LC-1)
      L=0
      DO 262 IK0=LK0+1,LK0+LC
      I0=(IK0-1)*LXTOT*LYTOT
      DO 261 IK1=LK1+1,LK1+LC
      I1=I0+(IK1-1)*LXTOT
      DO 260 IK2=LK2+1,LK2+LC
      I2=I1+IK2
      L=L+1
      IND1=KN(NUM1+L)
      IF(IND1.EQ.0) THEN
         IWRK(I2)=0
         GO TO 260
      ENDIF
      IF(IWRK(I2).EQ.-1) THEN
         IWRK(I2)=IND1
      ELSE IF(IWRK(I2).EQ.0) THEN
         KN(NUM1+L)=0
      ELSE IF(IWRK(I2).NE.IND1) THEN
         CALL XABORT('TRIPKN: FAILURE OF THE RENUMBERING ALGORITHM(1).')
      ENDIF
  260 CONTINUE
  261 CONTINUE
  262 CONTINUE
      NUM1=NUM1+LL
  270 CONTINUE
  271 CONTINUE
  272 CONTINUE
*----
*  CALCULATION OF PERMUTATION VECTOR IP AND RENUMBERING OF UNKNOWNS
*----
      DO 280 I=1,MAXIP
      IP(I)=0
  280 CONTINUE
      L4=0
      IF(NCODE(1).EQ.5) THEN
         K2MIN=1+LC/2
      ELSE
         K2MIN=1
      ENDIF
      DO 292 K0=1,LZTOT
      IK0=(K0-1)*LXTOT*LYTOT
      DO 291 K1=1,LYTOT
      IK1=IK0+(K1-1)*LXTOT
      DO 290 K2=K2MIN,LXTOT
      I=IWRK(IK1+K2)
      IF(I.LE.0) GO TO 290
      IF(I.GT.MAXIP) THEN
         CALL XABORT('TRIPKN: FAILURE OF THE RENUMBERING ALGORITHM(2).')
      ENDIF
      IF(IP(I).EQ.0) THEN
         L4=L4+1
         IP(I)=L4
      ENDIF
  290 CONTINUE
  291 CONTINUE
  292 CONTINUE
      DO 300 K=1,NUM1
      KNK=KN(K)
      IF(KNK.NE.0) KN(K)=IP(KNK)
  300 CONTINUE
      IF(IMPX.GT.0) WRITE (6,510) L4
      IF(IMPX.GT.2) WRITE (6,520) (VOL(I),I=1,LX*LY*LZ)
      IF(L4.EQ.0) THEN
         CALL XABORT('TRIPKN: FAILURE OF THE RENUMBERING ALGORITHM(3).')
      ENDIF
*
      IF(IMPX.LT.2) RETURN
      IF(IELEM.EQ.1) THEN
         WRITE (6,530)
         NUM1=0
         NUM2=0
         DO 310 KEL=1,LX*LY*LZ
         IF(MAT(KEL).LE.0) GO TO 310
         WRITE (6,540) KEL,(KN(NUM1+I),I=1,LL),(QFR(NUM2+I),I=1,6)
         NUM1=NUM1+LL
         NUM2=NUM2+6
  310    CONTINUE
      ELSE
         WRITE (6,590)
         NUM1=0
         DO 320 KEL=1,LX*LY*LZ
         IF(MAT(KEL).LE.0) GO TO 320
         WRITE (6,600) KEL,(KN(NUM1+I),I=1,LL)
         NUM1=NUM1+LL
  320    CONTINUE
         WRITE (6,610)
         NUM2=0
         DO 330 KEL=1,LX*LY*LZ
         IF(MAT(KEL).LE.0) GO TO 330
         WRITE (6,620) KEL,(QFR(NUM2+I),I=1,6)
         NUM2=NUM2+6
  330    CONTINUE
      ENDIF
*----
*  SCRATCH STORAGE DEALLOCATION
*----
      DEALLOCATE(IP,IWRK)
      RETURN
*
  500 FORMAT(/38H TRIPKN: PRIMAL FINITE ELEMENT METHOD.//7H NUMBER,
     1 27H OF ELEMENTS ALONG X AXIS =,I3/20X,14HALONG Y AXIS =,I3/
     2 20X,14HALONG Z AXIS =,I3)
  510 FORMAT(31H NUMBER OF UNKNOWNS PER GROUP =,I8)
  520 FORMAT(/20H VOLUMES PER ELEMENT/(1X,1P,10E13.4))
  530 FORMAT(/22H NUMBERING OF UNKNOWNS//8H ELEMENT,5X,7HNUMBERS,
     1 41X,23HVOID BOUNDARY CONDITION)
  540 FORMAT(1X,I6,2X,8I6,2X,1P,6E11.2)
  590 FORMAT(/22H NUMBERING OF UNKNOWNS//5H ELE-/5H MENT,3X,
     1 7HNUMBERS)
  600 FORMAT(1X,I6,2X,20I6/(9X,20I6))
  610 FORMAT(///24H VOID BOUNDARY CONDITION//8H ELEMENT,5X,3HQFR)
  620 FORMAT(1X,I6,4X,1P,6E11.2)
      END