From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001
From: stainer_t <thomas.stainer@oecd-nea.org>
Date: Mon, 8 Sep 2025 13:48:49 +0200
Subject: Initial commit from Polytechnique Montreal

---
 Trivac/src/BIVPRH.f | 469 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 469 insertions(+)
 create mode 100755 Trivac/src/BIVPRH.f

(limited to 'Trivac/src/BIVPRH.f')

diff --git a/Trivac/src/BIVPRH.f b/Trivac/src/BIVPRH.f
new file mode 100755
index 0000000..7034adc
--- /dev/null
+++ b/Trivac/src/BIVPRH.f
@@ -0,0 +1,469 @@
+*DECK BIVPRH
+      SUBROUTINE BIVPRH (MAXEV,MAXKN,IMPX,ISPLH,LX,IHEX,NCODE,ICODE,
+     1 ZCODE,MAT,SIDE,LL4,NELEM,VOL,KN,QFR,IQFR,BFR,MUW)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Numbering corresponding to a mesh corner finite difference or linear
+* Lagrangian finite element discretization of a 2-D hexagonal 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
+* MAXEV   maximum number of unknowns.
+* MAXKN   dimension for arrays KN, QFR and BFR.
+* IMPX    print parameter.
+* ISPLH   type of hexagonal mesh-splitting:
+*         =1: hexagonal elements; >1: 6*(ISPLH-1)**2 triangular elements
+*         per hexagon.
+* LX      number of hexagons.
+* IHEX    type of hexagonal boundary condition.
+* NCODE   type of boundary condition applied on each side
+*         (i=1: X-  i=2: X+  i=3: Y-  i=4: Y+):
+*         NCODE(I)=1: VOID;   NCODE(I)=2: REFL;   NCODE(I)=5: SYME;
+*         NCODE(I)=7: ZERO.
+* ICODE   physical albedo index on each side of the domain.
+* 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 hexagon.
+* SIDE    side of the hexagon.
+*
+*Parameters: output
+* LL4     order of system matrices.
+* NELEM   number of finite elements (hexagons or triangles) excluding
+*         the virtual elements.
+* VOL     volume of each hexagon.
+* KN      element-ordered unknown list.
+* QFR     element-ordered boundary conditions.
+* IQFR    element-ordered physical albedo indices.
+* BFR     element-ordered surface fractions.
+* MUW     compressed storage mode indices.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER MAXEV,MAXKN,IMPX,ISPLH,LX,IHEX,NCODE(4),ICODE(4),MAT(LX),
+     1 LL4,NELEM,KN(MAXKN),IQFR(MAXKN),MUW(LL4)
+      REAL ZCODE(4),SIDE,VOL(LX),QFR(MAXKN),BFR(MAXKN)
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER ISR(6,2),JCR(6),KK(6),ISRH(6,2),ISRT(3,2),ISRT2(3,2)
+      CHARACTER HSMG*131
+      LOGICAL LOG
+      INTEGER, DIMENSION(:), ALLOCATABLE :: IGAR,KN1
+      ALB(X)=0.5*(1.0-X)/(1.0+X)
+      DATA ISRH/2,1,4,5,6,3,1,4,5,6,3,2/
+      DATA ISRT/1,2,3,2,3,1/
+      DATA ISRT2/3,1,2,2,3,1/
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(IGAR(MAXEV),KN1(MAXKN))
+*
+      IF(IMPX.GT.0) WRITE(6,500)
+      IF(ISPLH.EQ.1) THEN
+         NSURF=6
+         MAXNH=LX
+      ELSE
+         NSURF=3
+         MAXNH=(6*(ISPLH-1)**2)*LX
+      ENDIF
+      CALL BIVSBH (MAXNH,MAXKN,IMPX,ISPLH,LX,SIDE,NELEM,IHEX,NCODE,
+     1 MAT,VOL,KN1,QFR)
+      IF(NELEM*NSURF.GT.MAXKN) THEN
+         WRITE(HSMG,'(28HBIVPRH: INSUFFICIENT MAXKN (,I7,10H). SHOULD ,
+     1   15HBE INCREASED TO,I7,1H.)') MAXKN,NELEM*NSURF
+         CALL XABORT(HSMG)
+      ENDIF
+*----
+*  PRODUCE STANDARD MESH CORNER FINITE DIFFERENCE NUMBERING.
+*----
+      DO 10 I=1,NELEM*(NSURF+1)
+      KN(I)=-98
+   10 CONTINUE
+      DO 20 IC=1,NSURF
+      IF(ISPLH.EQ.1) THEN
+         JCR(IC)=IC+3-(IC/4)*6
+         ISR(IC,1)=ISRH(IC,1)
+         ISR(IC,2)=ISRH(IC,2)
+      ELSE
+         IF(IC.EQ.1) JCR(1)=1
+         IF(IC.EQ.2) JCR(2)=3
+         IF(IC.EQ.3) JCR(3)=2
+         ISR(IC,1)=ISRT(IC,1)
+         ISR(IC,2)=ISRT(IC,2)
+      ENDIF
+   20 CONTINUE
+*----
+*  SET ZERO BOUNDARY CONDITIONS
+*----
+      NUM1=0
+      DO 30 KX=1,NELEM
+      DO 26 IC=1,NSURF
+      KY=ABS(KN1(NUM1+IC))
+      DO 25 I1=1,2
+      IF((KY.GT.NELEM).AND.(NCODE(1).EQ.7)) KN(NUM1+ISR(IC,I1))=-99
+   25 CONTINUE
+   26 CONTINUE
+      NUM1=NUM1+NSURF+1
+   30 CONTINUE
+*
+      SURFTOT=0.0
+      LL4=0
+      NUM1=0
+      DO 50 KX=1,NELEM
+      DO 40 IC=1,NSURF
+      IF(NSURF.EQ.6) THEN
+         BFR(NUM1+IC)=QFR(NUM1+IC)*QFR(NUM1+7)/(1.5*SQRT(3.0)*SIDE)
+         IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(NUM1+IC)=ALB(ZCODE(1))*QFR(NUM1+IC)*QFR(NUM1+7)/
+     1      (1.5*SQRT(3.0)*SIDE)
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(NUM1+IC)=QFR(NUM1+IC)*QFR(NUM1+7)/(1.5*SQRT(3.0)*SIDE)
+         ELSE
+            QFR(NUM1+IC)=0.0
+         ENDIF
+      ELSE
+         AA=SIDE/REAL(ISPLH-1)
+         BFR(NUM1+IC)=QFR(NUM1+IC)*QFR(NUM1+4)/(0.25*SQRT(3.0)*AA)
+         IF((NCODE(1).EQ.1).AND.(ICODE(1).EQ.0)) THEN
+            QFR(NUM1+IC)=ALB(ZCODE(1))*QFR(NUM1+IC)*QFR(NUM1+4)/
+     1      (0.25*SQRT(3.0)*AA)
+         ELSE IF(NCODE(1).EQ.1) THEN
+            QFR(NUM1+IC)=QFR(NUM1+IC)*QFR(NUM1+4)/(0.25*SQRT(3.0)*AA)
+         ELSE
+            QFR(NUM1+IC)=0.0
+         ENDIF
+      ENDIF
+      IQFR(NUM1+IC)=ICODE(1)
+      SURFTOT=SURFTOT+BFR(NUM1+IC)
+      KY0=KN1(NUM1+IC)
+      IF((KY0.LT.0).AND.(IHEX.NE.5).AND.(IHEX.NE.6)) GO TO 40
+      KY=ABS(KY0)
+      DO 35 I1=1,2
+      IND=ISR(IC,I1)
+      IF((KY.GT.NELEM).OR.(KY0.LT.0)) THEN
+         IF(KN(NUM1+IND).EQ.-98) THEN
+            LL4=LL4+1
+            KN(NUM1+IND)=LL4
+         ENDIF
+      ELSE
+         JND=ISR(JCR(IC),I1+1-(I1/2)*2)
+         IOF2=(KY-1)*(NSURF+1)+JND
+         IF(IOF2.GT.MAXKN) CALL XABORT('BIVPRH: ALGORITHM FAILURE 2.')
+         LOG=.FALSE.
+         IF(KN(IOF2).EQ.-99) THEN
+            KN(NUM1+IND)=-99
+         ELSE IF(KN(NUM1+IND).EQ.-98) THEN
+            LL4=LL4+1
+            KN(NUM1+IND)=LL4
+            IF(KY.NE.KX) KN(IOF2)=LL4
+            IF((KY.NE.KX).AND.(ISPLH.GT.1)) LOG=.TRUE.
+         ELSE IF(KN(NUM1+IND).EQ.-99) THEN
+            GO TO 35
+         ELSE IF((KN(IOF2).EQ.-98).AND.(KY.NE.KX)) THEN
+            KN(IOF2)=KN(NUM1+IND)
+            IF((KY.NE.KX).AND.(ISPLH.GT.1)) LOG=.TRUE.
+         ELSE IF((KN(NUM1+IND).NE.KN(IOF2)).AND.(KY.NE.KX)) THEN
+            CALL XABORT('BIVPRH: ALGORITHM FAILURE 3.')
+         ELSE IF((KY.NE.KX).AND.(ISPLH.GT.1)) THEN
+            LOG=.TRUE.
+         ENDIF
+         IF(LOG) THEN
+            KND=0
+            IF(JND.EQ.1) KND=2
+            IF(JND.EQ.2) KND=1
+            IF(JND.EQ.3) KND=3
+            KZ=KN1((KY-1)*4+ISRT2(JCR(IC),I1+1-(I1/2)*2))
+            IF((KZ.GT.0).AND.(KZ.LE.NELEM).AND.(KZ.NE.KY)) THEN
+               IF(KN((KZ-1)*4+KND).EQ.-99) THEN
+                  KN(NUM1+IND)=-99
+               ELSE
+                  KN((KZ-1)*4+KND)=KN(NUM1+IND)
+               ENDIF
+            ENDIF
+         ENDIF
+      ENDIF
+   35 CONTINUE
+   40 CONTINUE
+      KN(NUM1+NSURF+1)=KN1(KX*(NSURF+1))
+      NUM1=NUM1+NSURF+1
+   50 CONTINUE
+*----
+*  COMPUTE THE SURFACE FRACTIONS.
+*----
+      IF(SURFTOT.GT.0.0) THEN
+         DO 55 I=1,NUM1
+         BFR(I)=BFR(I)/SURFTOT
+   55    CONTINUE
+      ENDIF
+*
+      NUM1=0
+      DO 150 KX=1,NELEM
+      DO 60 IC=1,NSURF
+      KK(IC)=KN1(NUM1+IC)
+   60 CONTINUE
+      IF(ISPLH.EQ.1) THEN
+         IF((KX.EQ.1).AND.((IHEX.EQ.1).OR.(IHEX.EQ.10))) THEN
+            DO 70 I=1,6
+            KN(I)=KN(2)
+   70       CONTINUE
+         ELSE IF((KX.EQ.1).AND.((IHEX.EQ.2).OR.(IHEX.EQ.11))) THEN
+            DO 80 I=1,6
+            KN(I)=KN(2)
+   80       CONTINUE
+         ELSE IF((KX.EQ.1).AND.(IHEX.EQ.3)) THEN
+            KN(3)=KN(1)
+            KN(4)=KN(2)
+            KN(5)=KN(1)
+            KN(6)=KN(2)
+         ELSE IF((KX.EQ.1).AND.(IHEX.EQ.4)) THEN
+            KN(3)=KN(1)
+            KN(4)=KN(1)
+            KN(5)=KN(2)
+            KN(6)=KN(1)
+         ELSE IF((KX.EQ.1).AND.(IHEX.EQ.5)) THEN
+            KN(3)=KN(1)
+            KN(4)=KN(2)
+            KN(5)=KN(1)
+            KN(6)=KN(2)
+         ELSE IF((KX.EQ.1).AND.(IHEX.EQ.6)) THEN
+            KN(4)=KN(3)
+            KN(5)=KN(2)
+            KN(6)=KN(1)
+         ELSE IF((KK(1).EQ.-KK(5)).AND.(KK(2).EQ.-KK(4)).AND.
+     1   (KK(3).EQ.-KK(5)).AND.(KK(6).EQ.-KK(4))) THEN
+            DO 90 I=1,6
+            KN(NUM1+I)=KN(NUM1+6)
+   90       CONTINUE
+         ELSE IF((KK(1).EQ.-KK(3)).AND.(KK(4).EQ.-KK(2)).AND.
+     1   (KK(5).EQ.-KK(3)).AND.(KK(6).EQ.-KK(2))) THEN
+            DO 100 I=1,6
+            KN(NUM1+I)=KN(NUM1+4)
+  100       CONTINUE
+         ELSE IF((KK(1).EQ.-KK(6)).AND.(KK(2).EQ.-KK(5)).AND.
+     1   (KK(3).EQ.-KK(4))) THEN
+            KN(NUM1+3)=KN(NUM1+1)
+            KN(NUM1+6)=KN(NUM1+4)
+         ELSE IF((KK(5).EQ.-KK(4)).AND.(KK(6).EQ.-KK(3)).AND.
+     1   (KK(1).EQ.-KK(2))) THEN
+            KN(NUM1+4)=KN(NUM1+2)
+            KN(NUM1+5)=KN(NUM1+3)
+         ELSE IF((KK(1).EQ.-KK(3)).AND.(KK(4).EQ.-KK(3)).AND.
+     1   (KK(5).EQ.-KK(2)).AND.(KK(6).EQ.-KK(3))) THEN
+            KN(NUM1+1)=KN(NUM1+4)
+            KN(NUM1+2)=KN(NUM1+5)
+            KN(NUM1+3)=KN(NUM1+4)
+            KN(NUM1+6)=KN(NUM1+4)
+         ELSE IF((KK(2).EQ.-KK(6)).AND.(KK(3).EQ.-KK(5)).AND.
+     1   (KK(2).LT.0)) THEN
+            KN(NUM1+1)=KN(NUM1+2)
+            KN(NUM1+4)=KN(NUM1+3)
+            KN(NUM1+5)=KN(NUM1+6)
+         ELSE IF((KK(1).EQ.-KK(3)).AND.(KK(6).EQ.-KK(4)).AND.
+     1   (KK(1).LT.0)) THEN
+            KN(NUM1+1)=KN(NUM1+4)
+            KN(NUM1+2)=KN(NUM1+5)
+            KN(NUM1+3)=KN(NUM1+6)
+         ELSE IF((KK(4).EQ.-KK(2)).AND.(KK(5).EQ.-KK(1)).AND.
+     1   (KK(4).LT.0)) THEN
+            KN(NUM1+3)=KN(NUM1+2)
+            KN(NUM1+5)=KN(NUM1+4)
+            KN(NUM1+6)=KN(NUM1+1)
+         ELSE IF((KK(3).EQ.-KK(1)).AND.(KK(4).EQ.-KK(6)).AND.
+     1   (KK(3).LT.0)) THEN
+            KN(NUM1+4)=KN(NUM1+1)
+            KN(NUM1+5)=KN(NUM1+2)
+            KN(NUM1+6)=KN(NUM1+3)
+         ENDIF
+         DO 120 IC=1,NSURF
+         IF((KK(IC).LT.0).AND.((IHEX.EQ.5).OR.(IHEX.EQ.6)).AND.
+     1   (KX.GT.1)) THEN
+            IS=0
+            DO 110 I=1,6
+            IF(-KN1((-KK(IC)-1)*7+I).EQ.KX) IS=I
+  110       CONTINUE
+            IF(IS.EQ.0) CALL XABORT('BIVPRH: ALGORITHM FAILURE 4.')
+            KN((-KK(IC)-1)*7+ISRH(IS,2))=KN(NUM1+ISRH(IC,1))
+            KN((-KK(IC)-1)*7+ISRH(IS,1))=KN(NUM1+ISRH(IC,2))
+         ENDIF
+  120    CONTINUE
+      ELSE
+         IF((IHEX.NE.5).AND.(IHEX.NE.6)) THEN
+            IF((KK(1).EQ.-KK(2)).AND.(KK(1).LT.0)) THEN
+               KN(NUM1+1)=KN(NUM1+3)
+            ELSE IF((KK(1).EQ.-KK(3)).AND.(KK(1).LT.0)) THEN
+               KN(NUM1+2)=KN(NUM1+3)
+            ELSE IF((KK(2).EQ.-KK(3)).AND.(KK(2).LT.0)) THEN
+               KN(NUM1+2)=KN(NUM1+1)
+            ELSE IF((KK(2).EQ.-KK(1)).AND.(KK(2).LT.0)) THEN
+               KN(NUM1+3)=KN(NUM1+1)
+            ELSE IF((KK(3).EQ.-KK(1)).AND.(KK(3).LT.0)) THEN
+               KN(NUM1+3)=KN(NUM1+2)
+            ELSE IF((KK(3).EQ.-KK(2)).AND.(KK(3).LT.0)) THEN
+               KN(NUM1+1)=KN(NUM1+2)
+            ENDIF
+         ENDIF
+         DO 140 IC=1,NSURF
+         IF((KK(IC).LT.0).AND.((IHEX.EQ.5).OR.(IHEX.EQ.6))) THEN
+            IS=0
+            DO 130 I=1,3
+            IF(-KN1((-KK(IC)-1)*4+I).EQ.KX) IS=I
+  130       CONTINUE
+            IF(IS.EQ.0) CALL XABORT('BIVPRH: ALGORITHM FAILURE 5.')
+            KY=-KK(IC)
+            DO 135 I1=1,2
+            J1=I1+1-(I1/2)*2
+            IND=ISRT(IC,I1)
+            JND=ISRT(IS,J1)
+            KN((-KK(IC)-1)*4+JND)=KN(NUM1+IND)
+            KND=0
+            IF(JND.EQ.1) KND=2
+            IF(JND.EQ.2) KND=1
+            IF(JND.EQ.3) KND=3
+            KZ=KN1((-KK(IC)-1)*4+ISRT2(IS,J1))
+            IF((KZ.GT.0).AND.(KZ.LE.NELEM).AND.(KZ.NE.KY)) THEN
+               KNZ=KN((KZ-1)*4+KND)
+               KN((KZ-1)*4+KND)=KN(NUM1+IND)
+               IF(IHEX.EQ.6) THEN
+                  DO 132 L=1,NELEM
+                  DO 131 LC=1,NSURF
+                  IF(KN((L-1)*4+LC).EQ.KNZ) KN((L-1)*4+LC)=KN(NUM1+IND)
+  131             CONTINUE
+  132             CONTINUE
+               ENDIF
+            ENDIF
+  135       CONTINUE
+         ENDIF
+  140    CONTINUE
+      ENDIF
+      NUM1=NUM1+NSURF+1
+  150 CONTINUE
+      LL5=0
+      DO 170 I=1,MAXEV
+      IGAR(I)=0
+  170 CONTINUE
+      NUM1=0
+      DO 190 I=1,NELEM
+      DO 180 IC=1,NSURF
+      IND=KN(NUM1+IC)
+      IF(IND.GT.MAXEV) THEN
+         WRITE(HSMG,'(28HBIVPRH: INSUFFICIENT MAXEV (,I7,10H). SHOULD ,
+     1   15HBE INCREASED TO,I7,1H.)') MAXEV,IND
+         CALL XABORT(HSMG)
+      ELSE IF(IND.EQ.-98) THEN
+         CALL XABORT('BIVPRH: ALGORITHM FAILURE 6.')
+      ELSE IF(IND.EQ.-99) THEN
+         KN(NUM1+IC)=0
+      ELSE IF(IGAR(IND).EQ.0) THEN
+         LL5=LL5+1
+         IGAR(IND)=LL5
+      ENDIF
+  180 CONTINUE
+      NUM1=NUM1+NSURF+1
+  190 CONTINUE
+      NUM1=0
+      DO 210 I=1,NELEM
+      DO 200 IC=1,NSURF
+      IF(KN(NUM1+IC).NE.0) THEN
+         IF(IGAR(KN(NUM1+IC)).EQ.0) CALL XABORT('BIVPRH: ALGORITHM FAI'
+     1   //'LURE 7.')
+         KN(NUM1+IC)=IGAR(KN(NUM1+IC))
+      ENDIF
+  200 CONTINUE
+      NUM1=NUM1+NSURF+1
+  210 CONTINUE
+      LL4=LL5
+      IF(IMPX.GT.0) WRITE(6,570) LL4
+      IF(LL4.GT.MAXEV) THEN
+         WRITE(HSMG,'(28HBIVPRH: INSUFFICIENT MAXEV (,I7,10H). SHOULD ,
+     1   15HBE INCREASED TO,I7,1H.)') MAXEV,LL4
+         CALL XABORT(HSMG)
+      ENDIF
+      IF(LL4.GT.MAXEV) CALL XABORT('BIVPRH: INSUFFICIENT MAXEV.')
+*
+      IF((IMPX.GT.1).AND.(NSURF.EQ.6)) THEN
+         WRITE(6,510)
+         NUM1=0
+         DO 220 I=1,NELEM
+         WRITE(6,520) I,KN(NUM1+7),(KN(NUM1+J),J=1,6),(QFR(NUM1+J),
+     1   J=1,7)
+         NUM1=NUM1+7
+  220    CONTINUE
+         NUM1=0
+         WRITE (6,580)
+         DO 225 I=1,NELEM
+         IF(MAT(I).LE.0) GO TO 225
+         WRITE (6,590) I,(BFR(NUM1+J),J=1,6)
+         NUM1=NUM1+7
+  225    CONTINUE
+      ELSE IF((IMPX.GT.1).AND.(NSURF.EQ.3)) THEN
+         WRITE(6,530)
+         NUM1=0
+         DO 230 I=1,NELEM
+         WRITE(6,540) I,KN(NUM1+4),(KN(NUM1+J),J=1,3),(QFR(NUM1+J),
+     1   J=1,4),(BFR(NUM1+J),J=1,3)
+         NUM1=NUM1+4
+  230    CONTINUE
+      ENDIF
+*----
+*  COMPUTE THE SYSTEM MATRIX BANDWIDTH.
+*----
+      DO 240 I=1,LL4
+      MUW(I)=1
+  240 CONTINUE
+      NUM1=0
+      DO 270 K=1,NELEM
+      DO 260 I=1,NSURF
+      IND1=KN(NUM1+I)
+      IF(IND1.EQ.0) GO TO 260
+      DO 250 J=1,NSURF
+      IND2=KN(NUM1+J)
+      IF(IND2.EQ.0) GO TO 250
+      MUW(IND1)=MAX0(MUW(IND1),IND1-IND2+1)
+  250 CONTINUE
+  260 CONTINUE
+      NUM1=NUM1+NSURF+1
+  270 CONTINUE
+      IIMAX=0
+      DO 280 I=1,LL4
+      IIMAX=IIMAX+MUW(I)
+      MUW(I)=IIMAX
+  280 CONTINUE
+      IF(IMPX.GT.6) WRITE(6,550) 'MUW :',(MUW(I),I=1,LL4)
+      IF(IMPX.GT.2) WRITE(6,560) IIMAX
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(KN1,IGAR)
+      RETURN
+*
+  500 FORMAT(//52H BIVPRH: NUMBERING FOR A MESH CORNER FINITE DIFFEREN,
+     1 60HCE OR LINEAR LAGRANGIAN FINITE ELEMENT DISCRETIZATION IN HEX,
+     2 16HAGONAL GEOMETRY.)
+  510 FORMAT(/31H BIVPRH: NUMBERING OF UNKNOWNS./1X,30(1H-)/9X,
+     1 7HHEXAGON,3X,8HUNKNOWNS,29X,23HVOID BOUNDARY CONDITION,45X,
+     2 6HVOLUME)
+  520 FORMAT (1X,2I6,2X,6I6,2X,1P,6E11.2,5X,E13.6)
+  530 FORMAT(/31H BIVPRH: NUMBERING OF UNKNOWNS./1X,30(1H-)/9X,
+     1 7HHEXAGON,3X,8HUNKNOWNS,11X,23HVOID BOUNDARY CONDITION,12X,
+     2 6HVOLUME,13X,16HSURFACE FRACTION)
+  540 FORMAT (1X,2I6,2X,3I6,2X,1P,3E11.2,5X,E13.6,5X,3E10.2)
+  550 FORMAT(/1X,A5/(1X,20I6))
+  560 FORMAT(/52H NUMBER OF TERMS IN THE COMPRESSED SYSTEM MATRICES =,
+     > I6)
+  570 FORMAT(/39H BIVPRH: NUMBER OF UNKNOWNS PER GROUP =,I6/)
+  580 FORMAT (//17H SURFACE FRACTION//8H HEXAGON,5X,3HBFR)
+  590 FORMAT (3X,I4,4X,1P,6E10.2)
+      END
-- 
cgit v1.2.3