! needs post process ntape
!                             File begin.f90
!
      PROGRAM MAIN
!-->                         ***BEGIN.FEA***
!              A BEGINNER'S PROGRAM FOR FINITE ELEMENT ANALYSIS
!                            COPYRIGHT 1985, IN
!    "Finite Element Analysis for Undergraduates", Academic Press
!                             DR. J.E. AKIN
!   ************************************************************
!     First draft conversion from f77 to f90, Nov. 2004.
!     Needs dynamic array allocations, MATMUL, etc.
      COMMON R (1000), I (500)
      CHARACTER(70) TITLE
      DATA LIMITR, LIMITI / 4000, 1000 /
!     --- READ AND ECHO TITLE AND CONTROL ITEMS ---
      READ (5, * ) TITLE
      WRITE (6, * ) ' TITLE: ', TITLE
      CALL ZEROI (I, LIMITI)
      CALL ZEROA (R, LIMITR)
      READ  (5, * ) M, NE, NG, N, NSPACE, NPROP, NMAT, NFORCE
      WRITE (6, 30) M, NE, NG, N, NSPACE, NPROP, NMAT, NFORCE
   30 FORMAT (                                                  &
     & ' NUMBER OF NODES IN SYSTEM...................=',I4,/,   &
     & ' NUMBER OF ELEMENTS IN SYSTEM................=',I4,/,   &
     & ' NUMBER OF DEGREES OF FREEDOM PER NODE.......=',I4,/,   &
     & ' NUMBER OF NODES PER ELEMENT.................=',I4,/,   &
     & ' DIMENSION OF SPACE..........................=',I4,/,   &
     & ' NUMBER OF PROPERTIES PER MATERIAL...........=',I4,/,   &
     & ' NUMBER OF DIFFERENT MATERIALS...............=',I4,/,   &
     & ' NUMBER OF NODAL FORCE COMPONENT INPUTS......=' ,I4)
!     ----- REAL ARRAY STORAGE POINTERS, K1 TO K10 -------
!     K1-X, K2-COORD, K3-BC, K4-C, K5-CC, K6-D, K7-DD, K8-S,
!     K9-PROP, K10-SS
      K1 = 1
      K2 = K1 + M * NSPACE
      K3 = K2 + N * NSPACE
      K4 = K3 + M * NG
      K5 = K4 + N * NG
      K6 = K5 + M * NG
      K7 = K6 + N * NG
      K8 = K7 + M * NG
      K9 = K8 + (N * NG) **2
      K10 = K9 + NPROP * NMAT
      LEFTR = LIMITR - K10
      LIMITB = LEFTR / (M * NG)
!     ---- INTEGER ARRAY STORAGE POINTERS, L1 TO L5 -----
!     L1-IBC, L2-NODES, L3-MAT, L4-LNODES, L5-INDEX, L6=KODES
      L1 = 1
      L2 = L1 + M
      L3 = L2 + NE * N
      L4 = L3 + NE
      L5 = L4 + N
      L6 = L5 + N * NG
      LEFTI = LIMITI - NG
!     ----- CHECK STORAGE -----
      IF (LEFTR.GT.0.AND.LEFTI.GT.0.AND.LIMITB.GT.1) GOTO 40
      WRITE (6, * ) ' STORAGE EXCEEDED, STOP', LEFTR, LEFTI, LIMITB
      STOP ' STORAGE EXCEEDED, INCREASE COMMON'
!     -----  CALL REAL MAIN PROGRAM  -----
!     CALL  BEGIN  (M, NE, NG, N, NSPACE, NPROP, NMAT, LIMITB,
!    1            X, COORD, BC, C, CC, D, DD, S, PROP, SS, IBC,
!    2            NODES, MAT, LNODES, INDEX, KODES, NFORCE )
   40 CALL BEGIN (M, NE, NG, N, NSPACE, NPROP, NMAT, LIMITB, R (K1),  &
      R (K2), R (K3), R (K4), R (K5), R (K6), R (K7), R (K8), R (K9), &
      R (K10), I (L1), I (L2), I (L3), I (L4), I (L5), I (L6), NFORCE)
      WRITE (6, * ) ' NORMAL ENDING OF BEGIN.FEM'
      STOP
      END PROGRAM MAIN
!
      SUBROUTINE AT (N)
!  ******************************************************************
!      A DEBUGGING AID
!  ******************************************************************
      WRITE (6, 5) N
    5 FORMAT (' ======>>> AT ',I9)
      END SUBROUTINE AT
!
      SUBROUTINE BEGIN (M, NE, NG, N, NSPACE, NPROP, NMAT, MAXBAN, X, &
      COORD, BC, C, CC, D, DD, S, PROP, SS, IBC, NODES, MAT, LNODES,  &
      INDEX, KODES, NFORCE)
!   ******************************************************************
!           THE REAL MAIN PROGRAM
!   ******************************************************************
!     M = NUMBER OF NODES OF SYSTEM
!     N = NUMBER OF NODES PER ELEMENT
!     NELFRE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
!     NSPACE = DIMENSION OF SPACE
!     COORD = SPATIAL COORDINATES OF ELEMENT'S NODES
!     S  = ELEMENT SQUARE (STIFFNESS) MATRIX
!     C  = ELEMENT COLUMN (FORCING) MATRIX
      DIMENSION X (M, NSPACE), COORD (N, NSPACE), BC (M, NG), C (N*NG)&
      , CC (M * NG), D (N * NG), DD (M * NG), S (N * NG, N * NG),     &
      PROP (NPROP, NMAT), SS (M * NG, MAXBAN), IBC (M), NODES (NE, N),&
      MAT (NE), LNODES (N), INDEX (N * NG), KODES (NG)
      NTAPE = 8 !bb
      NGPLUS = NG + 1
      NDFREE = NG * M
      LEMFRE = NG * N
!     ---  READ  NODAL POINT DATA ---
!     NODAL POINT DATA:
!        J--------NODE NUMBER
!        IBC------BOUNDARY CONDITION FLAG, ONE DIGIT PER DOF
!                =0   --IF NO B.C. AT THAT DOF
!        X(J,K)---K SPATIAL COORDINATES
!        BC(J,L)--B. C. VALUES TO BE USED IF IBC DIGIT L NOT ZERO
      DO I = 1, M
        READ (5, 900) J, IBC (J), (X (J, K), K = 1, NSPACE), (BC (J, L), L = 1, NG)
!       PRINT *,  J, IBC (J), (X (J, K), K = 1, NSPACE), (BC (J, L), L = 1, NG)
      END DO
  900 FORMAT ( 2 I5, (7 F10.0) )
      WRITE (6, 901)
  901 FORMAT ( 15X, '* ECHO OF NODAL DATA *' ,/, ' B. C. CODES' )
      WRITE (6, 902) IBC
  902 FORMAT ( 5 I12, (/, 5 I12 ) )
      WRITE (6, 903)
  903 FORMAT ( ' X - COORDINATES' )
      WRITE (6, 904) (X (J, 1), J = 1, M)
  904 FORMAT ( 5( 1PE12.4 ), (/, 5 ( 1PE12.5) ) )
      IF (NSPACE.LT.2) GOTO 950
      WRITE (6, 905)
  905 FORMAT ( ' Y - COORDINATES' )
      WRITE (6, 904) (X (J, 2), J = 1, M)
      IF (NSPACE.LT.3) GOTO 950
      WRITE (6, 906)
  906 FORMAT ( ' Z - COORDINATES ')
      WRITE (6, 904) (X (J, 3), J = 1, M)
  950 WRITE (6, 907)
  907 FORMAT ( ' B. C. VALUES' )
      DO IG = 1, NG
        WRITE (6, 908) IG
        908 FORMAT ( ' FOR D.O.F. NUMBER', I2 )
        WRITE (6, 904) (BC (J, IG), J = 1, M)
      END DO
!     --- READ  ELEMENT DATA ---
!     ELEMENT DATA:
!        J------ELEMENT NUMBER
!        MAT----MATERIAL NUMBER OF ELEMENT
!        NODES--ELEMENT CONNECTIVITY
      WRITE (6, 910)
  910 FORMAT ( 15X, '* ECHO OF ELEMENTS *', /,                 &
     & ' ELEM  MAT  NODAL LIST' )
      DO I = 1, NE
        READ  (5, 1002) J, MAT (J), (NODES (J, K), K = 1, N)
        WRITE (6, 1002) J, MAT (J), (NODES (J, K), K = 1, N)
        1002 FORMAT ( 2 I5, (10 I5) )
      END DO
!     --- READ  PROPERTY DATA ---
!     PROPERTY DATA:
!        J----------MATERIAL NUMBER
!        PROP(I,J)--I TH PROPERTY FOR THE MATERIAL,1<I<NPROP
      WRITE (6, 911)
  911 FORMAT ( 15X, '* ECHO OF PROPERTY DATA *' , /,           &
     & ' MAT  PROPERTY LISTS' )
      DO I = 1, NMAT
      READ (5, 1003) J, (PROP (K, J), K = 1, NPROP)
 1003 FORMAT ( I5, ( 8 F10.0 ) )
      WRITE (6, 912) J, (PROP (K, J), K = 1, NPROP)
  912 FORMAT ( I5,  5 ( 1PE12.4 ), (/, 5X, 5 (1PE12.4) ) )
      END DO
!     READ FORCE DATA
      IF (NFORCE.LT.1) GOTO 10
      WRITE (6, 913)
  913 FORMAT ( 15X, '* ECHO NODAL FORCE DATA *', /,            &
     & ' NODE  DIR  VALUE' )
      DO I = 1, NFORCE
        READ (5, 1006) NODE, IDIR, VALUE
        1006 FORMAT ( 2 I5, F10.0 )
        WRITE (6, 914) NODE, IDIR, VALUE
        914 FORMAT ( 2 I5, 1PE12.4 )
        IDOF = NG * (NODE-1) + IDIR
        CC (IDOF) = VALUE
      END DO
   10 CONTINUE
!     DETERMINE SYSTEM HALF-BANDWIDTH
      CALL SYSBAN (NE, N, NG, IBW, NODES, LNODES)
      IF (IBW.GT.MAXBAN) WRITE (6, 1004)
 1004 FORMAT ('BANDWIDTH EXCEEDED')
!     CALCULATE AND ASSEMBLE ELEMENT MATRICES
      PRINT *, 'INPUT COMPLETED, BEGINNING ELEMENT ASSEMBLY'
      DO IE = 1, NE
        DO IG = 1, N
          LNODES (IG) = NODES (IE, IG)
          DO K = 1, NSPACE
            COORD (IG, K) = X (LNODES (IG), K)
          END DO
        END DO
        LMAT = MAT (IE)
        PRINT *, 'ELEMENT = ', IE, ', MATERIAL = ', LMAT
        CALL LSQCOL (N, NSPACE, NPROP, LEMFRE, COORD, PROP (1, LMAT), C, S)
        CALL INDXEL (N, NG, LEMFRE, LNODES, INDEX)
        CALL STORSQ (INDEX, S, SS, NDFREE, LEMFRE, IBW)
        CALL STORCL (INDEX, C, CC, NDFREE, LEMFRE)
        !a print *,' ELEM SQ : '
        !a call rprint (S, LEMFRE, LEMFRE, 1)
        !a print *,' ELEM COL: '
        !a call rprint (COL, 1, LEMFRE, 1)
      END DO
!     ASSEMBLY COMPLETED
      PRINT *,'ASSEMBLY COMPLETED, APPLYING BOUNDARY CONDITIONS'
      !a print *,'ASSEMBLED BAND '
      !a call rprint (SS, NDFREE, IBW, 1)
      !a print *,'ASSEMBLED COL  '
      !a call rprint (CC, 1, NDFREE, 1)
!     APPLY BOUNDARY CONDITIONS ON NODAL POINT VALUES
      DO J = 1, M
        IF (IBC (J) .EQ.0) GOTO 4
        CALL PTCODE (J, NG, IBC (J), KODES)
        CALL INDXPT (J, NG, INDEX)
        DO IG = 1, NG
          IF (KODES (IG) .EQ.0) GOTO 5
          INDX = INDEX (IG)
          CALL MODIFY (SS, NDFREE, IBW, INDX, BC (J, IG), CC)
        5 CONTINUE
        END DO
        4 CONTINUE
      END DO
      !a print *,'MODIFIED BAND '
      !a call rprint (SS, NDFREE, IBW, 1)
      !a print *,'MODIFIED COL  '
      !a call rprint (CC, 1, NDFREE, 1)

!     SOLVE FOR UNKNOWNS,   DD
      PRINT *, 'BC APPLIED, SOLVING SYSTEM EQUATIONS'
      CALL FACTOR (NDFREE, IBW, SS)
      CALL SOLVE  (NDFREE, IBW, SS, CC, DD)
!     ------- SOLUTION COMPLETED --------
!     PRINT RESULTS
      WRITE (6, 1005) NSPACE, NG
 1005 FORMAT ( '  ** COMPUTED RESULTS **', /,                  &
     & ' NODE,', I2, ' COORDINATES,', I2, ' DEGREES OF FREEDOM')
      CALL WRTPT (M, NG, NDFREE, NSPACE, X, DD, INDEX)
!     --------- POST PROCESSING --------------
      IF ( NTAPE .GT. 0 ) CALL POST(NE,N,LEMFRE,NDFREE,NODES,LNODE,INDEX,DD,D,NG)
      WRITE (6, * ) ' NORMAL ENDING OF BEGIN.FEA'
      STOP ' NORMAL ENDING OF BEGIN.FEA'
      END SUBROUTINE BEGIN
!
      SUBROUTINE DER8Q (S, T, DH)
!   ******************************************************************
!      LOCAL DERIVATIVES OF Q8 ELEMENT
!   ******************************************************************
      DIMENSION DH (2, 8)
      SP = 1. + S
      SM = 1. - S
      TP = 1. + T
      TM = 1. - T
      DH (1, 1) = - 0.25 * TM * (SM + SM + TM - 3.)
      DH (2, 1) = - 0.25 * SM * (TM + SM + TM - 3.)
      DH (1, 2) = 0.25 * TM * (SP + SP + TM - 3.)
      DH (2, 2) = - 0.25 * SP * (TM + SP + TM - 3.)
      DH (1, 3) = 0.25 * TP * (SP + SP + TP - 3.)
      DH (2, 3) = 0.25 * SP * (TP + SP + TP - 3.)
      DH (1, 4) = - 0.25 * TP * (SM + SM + TP - 3.)
      DH (2, 4) = 0.25 * SM * (TP + SM + TP - 3.)
      DH (1, 5) = - S * TM
      DH (2, 5) = - 0.5 * (1. - S * S)
      DH (1, 6) = 0.5 * (1. - T * T)
      DH (2, 6) = - T * SP
      DH (1, 7) = - S * TP
      DH (2, 7) = 0.5 * (1. - S * S)
      DH (1, 8) = - 0.5 * (1. - T * T)
      DH (2, 8) = - T * SM
      END SUBROUTINE DER8Q
!
      SUBROUTINE ELBAND (NUMNOD, NDF, IBW, LNODES)
!   ******************************************************************
!     FIND ELEMENT BANDWIDTH
!   ******************************************************************
      DIMENSION LNODES (NUMNOD)
      NLESS = NUMNOD-1
      DO I = 1, NLESS
        II = I + 1
        DO J = II, NUMNOD
          NEW = NDF * (IABS (LNODES (J) - LNODES (I) ) + 1)
          IF (NEW.GT.IBW) IBW = NEW
        END DO
      END DO
      END SUBROUTINE ELBAND
!
      SUBROUTINE FACTOR (N, IBW, S)
!   ******************************************************************
!     FACTORING OF A BANDED SYMMETRIC MATRIX
!   ******************************************************************
      DIMENSION S (N, IBW)
      TEMP = S (1, 1)
      DO J = 2, IBW
        S (1, J) = S (1, J) / TEMP
      END DO
      DO I = 2, N
        LL = I - 1
        NN = N - LL
       20 IF (NN.GT.IBW) NN = IBW
        DO J = 1, NN
          L = IBW - J
          SUM = 0.
          IF (L.EQ.0) GOTO 75
            IF (LL.LT.L) L = LL
         30 DO K = 1, L
            K1 = I - K
            K2 = 1 + K
            K3 = J + K
            SUM = SUM + S (K1, K2) * S (K1, K3) * S (K1, 1)
          END DO
       75 S (I, J) = S (I, J) - SUM
          IF (J.EQ.1) GOTO 100
          S (I, J) = S (I, J) / S (I, 1)
          100 CONTINUE
        END DO
      END DO
      END SUBROUTINE FACTOR
!
      SUBROUTINE  GAUSCO (N,A,W)
!     * * * * * * * * * * * * * * * * * * * * * * * * * *
!     GAUSSIAN QUADRATURE ABSCISSAE AND WEIGHT COEFFS
!     * * * * * * * * * * * * * * * * * * * * * * * * * *
      DIMENSION  A(N), W(N)
      DATA  NMAX/4/
!     N = NO. OF GAUSS POINTS IN ONE DIMENSION
!     A = ABSCISSAE, W = WEIGHTS  OF GAUSS POINTS
!     NMAX = MAX. NO. OF POINTS TABULATED HEREIN
      NLES1 = N - 1
      GO TO (10,20,30,40), NLES1
!--->>>  WARNING, ONE POINT RULE OMITTED FROM OPTIONS
!     N = 2
 10   A(1) = 0.5773502691896260
      A(2) = -A(1)
      W(1) = 1.0
      W(2) = 1.0
      RETURN
!     N = 3
 20   A(1) = 0.7745966692414830
      A(2) = 0.0
      A(3) = -A(1)
      W(1) = 0.5555555555555560
      W(2) = 0.8888888888888890
      W(3) = W(1)
      RETURN
!     N = 4
 30   A(1) = 0.8611363115940530
      A(2) = 0.3399810435848560
      A(3) = -A(2)
      A(4) = -A(1)
      W(1) = 0.3478548451374540
      W(2) = 0.6521451548625460
      W(3) = W(2)
      W(4) = W(1)
      RETURN
 40   WRITE (6,5000) NMAX
 5000 FORMAT ('N >',I3,' COEFFICIENTS NOT TABULATED')
      END SUBROUTINE GAUSCO

      SUBROUTINE  GAUS2D (NQP,GPT,GWT,NIP,PT,WT)
!     * * * * * * * * * * * * * * * * * * * * * * * * * *
!     USE 1-D GAUSSIAN DATA TO GENERATE
!     QUADRATURE DATA FOR A QUADRILATERAL
!     * * * * * * * * * * * * * * * * * * * * * * * * * *
  !b  DIMENSION GPT(NQP), GWT(NQP), PT(2,1), WT(1)
      DIMENSION GPT(NQP), GWT(NQP), PT(2,*), WT(*) !b
      !a DIMENSION GPT(NQP), GWT(NQP), PT(2,16), WT(16) !b
!     DIM PT(2,NIP), WT(NIP)
!     NQP = NO. TABULATED 1-D POINTS, NIP = NQP*NQP
!     GPT, GWT = TABULATED 1-D QUADRATURE RULES
!     PT, WT = CALCULATED COORDS AND WEIGHTS IN QUAD
      CALL  GAUSCO (NQP,GPT,GWT)
      NIP = NQP*NQP
      K = 0
      DO 20  I = 1,NQP
        DO 10  J = 1,NQP
          K = K + 1
          WT(K) = GWT(I)*GWT(J)
          PT(1,K) = GPT(J)
          PT(2,K) = GPT(I)
   10   CONTINUE
   20 CONTINUE
      END SUBROUTINE GAUS2D

      SUBROUTINE GDERIV (NSPACE, N, AJINV, DELTA, GLOBAL)
!   ******************************************************************
!      CONVERT LOCAL DERIVATIVES TO PHYSICAL DERIVATIVES
!   ******************************************************************
      DIMENSION AJINV (NSPACE, NSPACE), DELTA (NSPACE, N), &
                GLOBAL (NSPACE, N)
      DO I = 1, NSPACE
        DO J = 1, N
          SUM = 0.0
          DO K = 1, NSPACE
            SUM = SUM + AJINV (I, K) * DELTA (K, J)
          END DO
          GLOBAL (I, J) = SUM
        END DO
      END DO
      END SUBROUTINE GDERIV
!
      SUBROUTINE I2BY2 (A, AINV, DET)
!   ******************************************************************
!      INVERT 2 BY 2 MATRIX
!   ******************************************************************
      DIMENSION A (2, 2), AINV (2, 2)
      DET = A (1, 1) * A (2, 2) - A (1, 2) * A (2, 1)
      AINV (1, 1) = A (2, 2) / DET
      AINV (1, 2) = - A (1, 2) / DET
      AINV (2, 1) = - A (2, 1) / DET
      AINV (2, 2) = A (1, 1) / DET
      END SUBROUTINE I2BY2
!
      SUBROUTINE I3BY3 (A, AINV, DET)
!   ******************************************************************
!      INVERT 3 BY 3 MATRIX
!   ******************************************************************
      DIMENSION A (3, 3), AINV (3, 3)
      AINV (1, 1) = A (2, 2) * A (3, 3) - A (3, 2) * A (2, 3)
      AINV (2, 1) = - A (2, 1) * A (3, 3) + A (3, 1) * A (2, 3)
      AINV (3, 1) = A (2, 1) * A (3, 2) - A (3, 1) * A (2, 2)
      AINV (1, 2) = - A (1, 2) * A (3, 3) + A (3, 2) * A (1, 3)
      AINV (2, 2) = A (1, 1) * A (3, 3) - A (3, 1) * A (1, 3)
      AINV (3, 2) = - A (1, 1) * A (3, 2) + A (3, 1) * A (1, 2)
      AINV (1, 3) = A (1, 2) * A (2, 3) - A (2, 2) * A (1, 3)
      AINV (2, 3) = - A (1, 1) * A (2, 3) + A (2, 1) * A (1, 3)
      AINV (3, 3) = A (1, 1) * A (2, 2) - A (2, 1) * A (1, 2)
      DET = A (1, 1) * AINV (1, 1) + A (1, 2) * AINV (2, 1) &
          + A (1, 3) * AINV (3, 1)
      DO J = 1, 3
        DO I = 1, 3
          AINV (I, J) = AINV (I, J) / DET
        END DO
      END DO
      END SUBROUTINE I3BY3
!
      SUBROUTINE INDXEL (N, NG, LEMFRE, LNODES, INDEX)
!   ******************************************************************
!       COMPUTE DEGREES OF FREEDOM ON AN ELEMENT
!   ******************************************************************
      DIMENSION INDEX (LEMFRE), LNODES (N)
      DO K = 1, N
        DO IG = 1, NG
          IELM = NG * (K - 1) + IG
          INDEX (IELM) = NG * (LNODES (K) - 1) + IG
        END DO
      END DO
      END SUBROUTINE INDXEL
!
      SUBROUTINE INDXPT (IPT, NG, INDEX)
!   ******************************************************************
!        COMPUTE DEGREES OF FREEDOM NUMBERS AT A POINT
!   ******************************************************************
      DIMENSION INDEX (NG)
        DO J = 1, NG
        INDEX (J) = NG * (IPT - 1) + J
      END DO
      END SUBROUTINE INDXPT
!
      SUBROUTINE INVDET (AJ, AJINV, DET, NSPACE)
!   ******************************************************************
!      INVERT 2 BY 2 OR 3 BY 3 MATRIX AND RETURN DETERMINANT
!   ******************************************************************
      DIMENSION AJ (NSPACE, NSPACE), AJINV (NSPACE, NSPACE)
      GOTO (100, 200, 300, 400), NSPACE
  100 DET = AJ (1, 1)
      AJINV (1, 1) = 1.0 / DET
      RETURN
  200 CALL I2BY2 (AJ, AJINV, DET)
      RETURN
  300 CALL I3BY3 (AJ, AJINV, DET)
      RETURN
  400 WRITE (6, 1000)
 1000 FORMAT (' INVALID DATA SUBR INVDET')
      END SUBROUTINE INVDET
!
      SUBROUTINE JACOB (N, NSPACE, DELTA, COORD, AJ)
!   ******************************************************************
!      COMPUTE GEOMETRIC JACOBIAN
!   ******************************************************************
      DIMENSION DELTA (NSPACE, N), COORD (N, NSPACE), &
                AJ (NSPACE, NSPACE)
      DO I = 1, NSPACE
        DO J = 1, NSPACE
          SUM = 0.0
          DO K = 1, N
            SUM = SUM + DELTA (I, K) * COORD (K, J)
          END DO
          AJ (I, J) = SUM
        END DO
      END DO
      END SUBROUTINE JACOB
!
!      SUBROUTINE  LSQCOL (N,NSPACE,NPROP,LEMFRE,COORD,PROP,C,S)
!!     **********************************************************
!!                GENERATE ELEMENT SQUARE MATRIX
!!     **********************************************************
!      DIMENSION  COORD(N, NSPACE), C(LEMFRE), PROP(NPROP),  &
!                 S(LEMFRE, LEMFRE)
!!     .....................................................
!!     *** PROBLEM DEPENDENT STATEMENTS FOLLOW ***
!!     .....................................................
!!!==>>  HEAT CONDUCTION BAR WITH SOURCE
!!!        K*U,XX + Q = 0,          NPROP = 3
!!!     K - CONDUCTIVITY, A - AREA, Q - SOURCE PER UNIT VOLUME
!!!     USE A = 1 FOR DATA PER UNIT LENGTH
!!      DIST = COORD(2,1) - COORD(1,1)
!!      A = PROP(2)
!!      IF ( A .LE. 0.0 )  A = 1.0
!!      EK = PROP(1)*A/DIST
!!      S(1,1) =  EK
!!      S(2,1) = -EK
!!      S(1,2) = -EK
!!      S(2,2) =  EK
!!!-->  SOURCE TERM
!!      Q = PROP(3)
!!      EC = 0.5*Q*DIST*A
!!      C(1) = EC
!!      C(2) = EC
!
!! Linear Triangle, K_xx U,xx + 2 K_xy U,xy + K_yy U,yy + Q = 0
!  REAL :: X_I, X_J, X_K, Y_I, Y_J, Y_K     ! Global coordinates
!  REAL :: A_I, A_J, A_K, B_I, B_J, B_K     ! Standard geometry
!  REAL :: C_I, C_J, C_K, X_CG, Y_CG, TWO_A ! Standard geometry
!  REAL :: B (NSPACE, LEMFRE), E (NSPACE, NSPACE)
!
!!   DEFINE NODAL COORDINATES, CCW: I, J, K
!    X_I = COORD (1,1)
!    X_J = COORD (2,1)
!    X_K = COORD (3,1)
!    Y_I = COORD (1,2)
!    Y_J = COORD (2,2)
!    Y_K = COORD (3,2)
!
!!   DEFINE CENTROID COORDINATES (QUADRATURE POINT)
!   X_CG = (X_I + X_J + X_K)/3.d0
!   Y_CG = (Y_I + Y_J + Y_K)/3.d0
!
!!   GEOMETRIC PARAMETERS: H_I (X,Y) = (A_I + B_I*X + C_I*Y)/TWO_A
!   A_I = X_J * Y_K - X_K * Y_J
!   B_I = Y_J - Y_K
!   C_I = X_K - X_J
!   A_J = X_K * Y_I - X_I * Y_K
!   B_J = Y_K - Y_I
!   C_J = X_I - X_K
!   A_K = X_I * Y_J - X_J * Y_I
!   B_K = Y_I - Y_J
!   C_K = X_J - X_I
!
!!   CALCULATE TWICE ELEMENT AREA
!   TWO_A = A_I + A_J + A_K              ! = B_J*C_K - B_K*C_J also
!
!!   DEFINE 2 BY 3 GRADIENT MATRIX, B (= DGH)
!   B (1, 1:3) = (/ B_I, B_J, B_K /) / TWO_A  ! DH/DX, row 1
!   B (2, 1:3) = (/ C_I, C_J, C_K /) / TWO_A  ! DH/DY, row 2
!
!!   DEFINE PROPERTIES: 1-K_xx, 2-K_yy, 3-K_xy, 4-Source, 5-thick
!   E (1, 1) = PROP (1)
!   E (1, 2) = 0.d0
!   E (2, 1) = 0.d0
!   E (2, 2) = PROP (1)
!
!!   CONDUCTION MATRIX, WITH CONSTANT JACOBIAN
!   S = MATMUL ( TRANSPOSE (B), MATMUL (E, B) ) * TWO_A * 0.5d0
!
!!   SOURCE VECTOR: C(1:3) = SOURCE_PER_UNIT_AREA * AREA / 3
!   C = PROP (2) * TWO_A / 6.d0
!
!!  End of application dependent code
!      END SUBROUTINE  LSQCOL

      SUBROUTINE MADD (A, B, C, M, N)
!   ******************************************************************
!     ADD TWO MATRICES:  C = A + B
!   ******************************************************************
      DIMENSION A (1), B (1), C (1)
      MN = M * N
      DO I = 1, MN
        C (I) = A (I) + B (I)
      END DO
      END SUBROUTINE MADD
!
      SUBROUTINE  MMULT (A,B,C,L,M,N) !b
!   *****************************************************************
!      MULTIPLY TWO MATRICES:  C = A * B       !b
!   *****************************************************************
      DIMENSION A(L,M), B(M,N), C(L,N) !b
      ZERO = 0.0
      DO 30  I = 1,L
        DO 20  J = 1,N
        SUM = ZERO
          DO 10  K = 1,M
            BKJ = B(K,J)
            IF ( BKJ.EQ.ZERO )  GO TO 10
              AIK = A(I,K)
              IF ( AIK.EQ.ZERO )  GO TO 10
                SUM = SUM + AIK*BKJ
 10       END DO
          C(I,J) = SUM
 20     END DO
 30   END DO
      END SUBROUTINE MMULT  !b

      !b SUBROUTINE MMULT (A, B, C, N)
      SUBROUTINE SQMULT (A, B, C, N)  !b
!   *****************************************************************
!      MULTIPLY TWO SQUARE MATRICES:  C = A * B
!   *****************************************************************
      DIMENSION A (N, N), B (N, N), C (N, N)
      DO I = 1, N
        DO J = 1, N
          C (I, J) = 0.0
          DO K = 1, N
            C (I, J) = C (I, J) + A (I, K) * B (K, J)
          END DO
        END DO
      END DO
      !b END SUBROUTINE MMULT
      END SUBROUTINE SQMULT                          !b
!
      SUBROUTINE MODIFY (SS, NE, IBW, I, U, F)
!   ******************************************************************
!     MODIFY COMPACTED SQ MATRIX FOR ESSENTIAL BOUNDARY CONDITIONS
!   ******************************************************************
      DIMENSION SS (NE, IBW), F (NE)
!        AND SYSTEM COLUMN MATRIX
      IL = I - IBW + 1
      IF (IL.LE.0) IL = 1
      DO J = IL, I
        JJ = I - J + 1
        F (J) = F (J) - SS (J, JJ) * U
        SS (J, JJ) = 0.0
      END DO
      DO JJ = 2, IBW
        J = I + JJ - 1
        !b F (J) = F (J) - SS (I, JJ) * U
        IF ( J .LE. NE ) F (J) = F (J) - SS (I, JJ) * U !b author
        SS (I, JJ) = 0.0
      END DO
      SS (I, 1) = 1.0
      F (I) = U
      END SUBROUTINE MODIFY
!
      SUBROUTINE MTMULT (A, B, C, N)
!   ******************************************************************
!     (N*N) TRANPOSE MULTIPLICATION.  C=AT*B
!             WHERE AT = A TRANPOSE
!   ******************************************************************
      DIMENSION A (N, N), B (N, N), C (N, N)
      DO I = 1, N
        DO J = 1, N
          C (I, J) = 0.0
          DO K = 1, N
            C (I, J) = C (I, J) + A (K, I) * B (K, J)
          END DO
        END DO
      END DO
      END SUBROUTINE MTMULT
!
      SUBROUTINE PTCODE (JPT, NG, KODE, KODES)
!   ******************************************************************
!     EXTRACT LIST OF BOUNDARY CONDITION CODES AT A POINT
!   ******************************************************************
      DIMENSION KODES (NG)
      NGPLUS = NG + 1
      IOLD = KODE
      ISUM = 0
      DO I = 1, NG
        II = NGPLUS - I
        INEW = IOLD / 10
        IK = IOLD-INEW * 10
        ISUM = ISUM + IK * 10** (I - 1)
        IOLD = INEW
        KODES (II) = IK
      END DO
      IF (KODE.GT.ISUM) WRITE (6, 5000) JPT
 5000 FORMAT ('BC NOT RIGHT JUSTIFIED AT NODE',I5)
      END SUBROUTINE PTCODE
!
      SUBROUTINE SHP8Q (S, T, H)
!   ******************************************************************
!      LOCAL INTERPLOATIONS OF Q8 ELEMENT
!   ******************************************************************
      DIMENSION H (8)
      SP = 1. + S
      SM = 1. - S
      TP = 1. + T
      TM = 1. - T
      H (1) = 0.25 * SM * TM * (SM + TM - 3.)
      H (2) = 0.25 * SP * TM * (SP + TM - 3.)
      H (3) = 0.25 * SP * TP * (SP + TP - 3.)
      H (4) = 0.25 * SM * TP * (SM + TP - 3.)
      H (5) = 0.5 * TM * (1. - S * S)
      H (6) = 0.5 * SP * (1. - T * T)
      H (7) = 0.5 * TP * (1. - S * S)
      H (8) = 0.5 * SM * (1. - T * T)
      END SUBROUTINE SHP8Q
!
      SUBROUTINE SOLVE (N, IBW, S, P, D)
!   ******************************************************************
!     PART TWO OF CHOLESKY-GAUSSIAN SOLUTION
!     FOLLOWS FACTOR OF SYMMETRIC EQUATIONS
!   ******************************************************************
      DIMENSION S (N, IBW), P (N), D (N)
      D (1) = P (1) / S (1, 1)
      DO I = 2, N
        II = I + 1
        J = II - IBW
        IF (II.LE.IBW) J = 1
        IK = I - 1
        SUM = 0.
        DO K = J, IK
          KK = II - K
          SUM = SUM + S (K, KK) * S (K, 1) * D (K)
        END DO
        D (I) = (P (I) - SUM) / S (I, 1)
      END DO
      DO NN = 2, N
        I = N + 1 - NN
        LL = I - 1
        J = LL + IBW
        IF (J.GT.N) J = N
          L = I + 1
          SUM = 0.
          DO K = L, J
            KK = K - LL
            SUM = SUM + S (I, KK) * D (K)
          END DO
          D (I) = D (I) - SUM
      END DO
      END SUBROUTINE SOLVE
!
      SUBROUTINE STORCL (N, C, CC, NE, NDF)
!   ******************************************************************
!     STORE ELEMENT COLUMN MATRIX IN SYSTEM COLUMN MATRIX
!   ******************************************************************
      DIMENSION C (NDF), CC (NE), N (NDF)
      DO I = 1, NDF
        J = N (I)
        CC (J) = C (I) + CC (J)
      END DO
      END SUBROUTINE STORCL
!
      SUBROUTINE STORSQ (INDEX, S, SS, NE, NDF, IBW)
!   ******************************************************************
!     STORE ELEM SQ MATRIX IN UPPER HALF BAND OF COMPACT SYS SQ MATRIX
!   ******************************************************************
      DIMENSION INDEX (NDF), S (NDF, NDF), SS (NE, IBW)
!     I AND J ARE THE ROW AND COLUMN POSITIONS IN THE
!        UNPACKED SYSTEM MATRIX
!     JJ IS THE COLUMN POSITION IN THE PACKED SYSTEM MATRIX
      DO L = 1, NDF
        I = INDEX (L)
        DO K = 1, NDF
          J = INDEX (K)
          IF (I.GT.J) GOTO 1
            JJ = J - I + 1
            SS (I, JJ) = SS (I, JJ) + S (L, K)
        1 CONTINUE
        END DO
      END DO
      END SUBROUTINE STORSQ
!
      SUBROUTINE SYSBAN (NE, N, NG, IBW, NODES, LNODES)
!   ******************************************************************
!     FIND SYSTEM HALF BANDWIDTH FROM ELEMENT DATA
!   ******************************************************************
      DIMENSION NODES (NE, N), LNODES (N)
      IBW = 1
      DO I = 1, NE
        DO J = 1, N
          LNODES (J) = NODES (I, J)
        END DO
        CALL ELBAND (N, NG, IBW, LNODES)
      END DO
      END SUBROUTINE SYSBAN
!
      SUBROUTINE WRTPT (M, NG, NDF, NSPACE, X, DD, INDEX)
!   ******************************************************************
!     PRINT SOLUTION D AT A SPECIFIED POINT
!   ******************************************************************
      DIMENSION X (M, NSPACE), DD (NDF), INDEX (NG)
      DO I = 1, M
        CALL INDXPT (I, NG, INDEX)
        WRITE (6, 1004) I, (X (I, L), L = 1, NSPACE), (DD (INDEX (K)),&
        K = 1, NG)
        1004 FORMAT ( I5, (7 ( 1X, 1PE12.5 )) )
      END DO
      END SUBROUTINE WRTPT
!
      SUBROUTINE ZEROA (N, A)
!   ******************************************************************
!     ZERO A REAL ARRAY ! not needed, simply set A = 0
!   ******************************************************************
      DIMENSION A (N)
      DO J = 1, N
        A (J) = 0.0
      END DO
      END SUBROUTINE ZEROA
!
      SUBROUTINE ZEROI (N, I)
!   ******************************************************************
!     ZERO AN INTEGER ARRAY ! not needed, simply set I = 0
!   ******************************************************************
      DIMENSION I (N)
      DO J = 1, N
        I (J) = 0
      END DO
      END SUBROUTINE ZEROI

      SUBROUTINE LSQCOL (N,NSPACE,NPROP,LEMFRE,COORD,PROP,C,S) 
! * * * * * * * * * * * * * * * * * * * ** * * * * * * * * * * * * * 
!   GENERATE ELEMENT SQUARE MATRIX 
! * * * * * * * * * * * * * * * * * * * ** * * * * * * * * * * * * * 
!DP IMPLICIT REAL*8(A-H, O-Z) 
      DIMENSION  COORD (N, NSPACE),D(LEMFRE),PROP(NPROP), &
                 S(LEMFRE,LEMFRE) 
!
! *** ELSQ PROBLEM DEPENDENT STATEMENTS FOLLOW *** 
!
! SOLUTION OF LAPLACE EQUATION IN TWO DIMENSIONS 
! USING 8 NODE ISOPARAMETRIC QUADRILATERIAL ELEMENT 
! NSPACE = 2, N = 8, NG = 1, NELFRE = 8 
      COMMON / ELARG/ GPT (4), GWT (4), PT (2,16),WT (16), NIP 
      !b COMMON /ELARGE2/ XY(2), H(8), GRAD(2), AJ (2,2), AJINV (2,2), &
      COMMON /ELARG2/ XY(2), H(8), GRAD(2), AJ (2,2), AJINV (2,2), & !b
                       DELTA(2,8),GLOBAL (2,8),COL(8) 
      !b COMMON /ELARGE3/ XK, YK 
      COMMON /ELARG3/ XK, YK, DPDZ !b spelling, add source dp/dz
      EXTERNAL LAPL8Q !b
      DATA KALL, NQPOLD/1, 0/ 
      IF (KALL.EQ.0)  GO TO 10 
! DEFINE PROPERTIES 
  10  XK = PROP (1) 
      YK = PROP (2) 
      !b NQP =PROP (3) 
      DPDZ = PROP (3) !b  ! pressure gradient
      NQP = 2 !b
! DEFAULT NQP = 2 SINCE LPROP(1) = 0    IF NLPFIX = 0 
      !b IF (NQP.LT.2)   NQP = 2 
      !b IF (NQP.GT.4)   NQP = 4 
      IF (NQPOLD.EQ.NQP)   GO TO 20 
! CALCULATE QUADRATURE DATA (  IF REQUIRED  ) 
      NQPOLD = NQP 
      CALL GAUS2D (NQP, GPT, GWT, NIP, PT, WT) 
! STORE NUMBER OF POINTS FOR GRADIENT CALCULATIONS 
  20  CONTINUE
      !a WRITE (6,*)  'NUMBER OF QUADRATURE POINTS = ', NIP 
! NUMERICAL INTEGRATION LOOP 
      CALL ISOPAR (N, NSPACE, LEMFRE, NIP, S, C, PT, WT, H, &
                   DELTA, GLOBAL, COORD, XY, AJ, AJINV, LAPL8Q) 
      END  SUBROUTINE LSQCOL 

      SUBROUTINE LAPL8Q (WT,DET,H,DGH,XPT,N, NSPACE, &
                         LEMFRE, COL, SQ) 
! * * * * * * * * * * * * * * * * * * * ** * * * * * * * * * * * * * 
! PROBLEM DEPENDENT INTEGRAND EVALUATION IN 
! AN ISOPARAMETRIC ELEMENT 
! * * * * * * * * * * * * * * * * * * * ** * * * * * * * * * * * * * 
!DP IMPLICIT REAL*8 (A-H, O-Z) 
      DIMENSION COL(LEMFRE),SQ(LEMFRE,LEMFRE),H(N), &
                DGH(NSPACE,N), XPT(NSPACE) 
!
! *** LAPL8Q PROBLEM DEPENDENT STATEMENTS FOLLOW *** 
!
      !b COMMON /ELARG3/ XK, YK
      COMMON /ELARG3/ XK, YK, DPDZ !b omitted press grad
! STORE COORDINATE AND DERIVATIVE MATRIX AT THE POINT 
      !b WRITE (6,*)    XPT, DGH  ! ntape1
      IF ( NTAPE > 0 ) WRITE (NTAPE,*)    XPT, DGH  !bb ntape1
! EVALUATE PRODUCTS 
      !b DET WT = DET*WT  ! spell
      DETWT = DET*WT !b
      DO 40 J = 1,N 
        DO 30 I =1, J 
          SQ (I,J) = SQ (I,J) + DETWT*(XK*DGH(1,I)*DGH(1,J) &
                              +YK*DGH(2,I)*DGH(2,J)) 
          SQ (J,I) = SQ(I,J) 
  30    END DO
        COL (J) = COL (J) + DETWT*DPDZ*H(J) !b pressure grad source
  40  END DO
      END SUBROUTINE LAPL8Q
  
      SUBROUTINE ISOPAR (N, NSPACE, LEMFRE, NIP, SQ, C, QPT, &
                         QWT, H, DLH, DGH, COORD, XPT, AJ,   & !b coord
                         AJINV,LAPL8Q) 
                     !b  QWT, H, DLH, DGH, COCRD, XPT, AJ,   & ! spell error
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
! NUMERICAL INTEGRATION IN AN ISOPARAMETRIC ELEMENT 
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
!DP IMPLICIT REAL*8 (A-H, O-Z) 
      external  lapl8q !b
      !b DIMENSION COL(LEMFRE), SQ(LEMFRE,LEMFRE), QWT(NIP), &
      DIMENSION C(LEMFRE), SQ(LEMFRE,LEMFRE), QWT(NIP), &    !b
                QPT(NSPACE,NIP), H(N), DLH(NSPACE,N),     &
                DGH(NSPACE,N), COORD(N,NSPACE),XPT(NSPACE), &
                AJ(NSPACE,NSPACE), AJINV(NSPACE,NSPACE) 
! ZERO INTEGRANDS 
      CALL ZEROA (LEMFRE, C) 
      LSQ = LEMFRE*LEMFRE 
      CALL ZEROA (LSQ, SQ) 
! BEGIN INTEGRATION 
      DO 100    IP = 1,NIP 
        !a print *, 'ISOPAR: IP = ', IP
! EVALUATE INTERPOLATION FUNCTIONS 
        S = QPT (1, IP) !b fatal error
        T = QPT (2, IP) !b fatal error
        CALL SHP8Q (S,T,H) !b S, T undefined
! FIND GLOBAL COORD, XPT = H*COORD 
        !b CALL MMULT (H, COORD, XPT, I, N, NSPACE)  ! fatal error
        CALL MMULT (H, COORD, XPT, 1, N, NSPACE)  !b fatal error
! FIND LOCAL DERIVATIVES 
        !b CALL DER8Q(S,T,DH)  ! typo
        CALL DER8Q(S,T,DLH) !b
! FIND JACOBIAN AT THE PT 
        CALL JACOB (N, NSPACE, DLH, COORD, AJ) 
! FORM INVERSE AND DETERMINANT OF JACOBIAN 
        CALL INVDET (AJ, AJINV, DET, NSPACE) 
! EVALUATE GLOBAL DERIVATIVES 
        CALL GDERIV (NSPACE, N, AJINV, DLH, DGH) 
!  UPDATE ELEMENT MATRICES AT INTEGRATION PT !b omitted
        CALL LAPL8Q (QWT(IP),DET,H,DGH,XPT,N, NSPACE,LEMFRE,C,SQ)  !b omitted
  100 END DO
      END SUBROUTINE ISOPAR 

      SUBROUTINE POST(NE,N,LEMFRE,NDFREE,NODES,LNODE,INDEX,DD,D,NG) 
! * * * * * * * * * * * * * * * * * * * ** * * * * * * * * * * * * * 
!  POST PROCESS AN ELEMENT
! * * * * * * * * * * * * * * * * * * * ** * * * * * * * * * * * * * 
      DIMENSION DD(NDFREE), D(LEMFRE),NODES(NE,N),LNODE(N),INDEX(LEMFRE) 
! LOOP OVER ELEMENTS 
      DO 10 IE=1,NE 
! EXTRACT ELEMENT'S NODE 
        DO 20 IG=1,N 
          LNODE(I)=NODES(IE,IG) 
  20   END DO   
!     EXTRACT NODAL PARAMETERS OF THE ELEMENT 
        CALL INDXEL(N,NG,LEMFRE,LNODE,INDEX)     
!     EXTRACT ELEMENT DOF FROM SYSTEM DOF 
        DO 30 I=1,LEMFRE 
          D(I)=0.D0 
          IF(INDEX(I).GT.0)  D(I)=DD(INDEX(I)) 
  30    END DO        
!     PROBLEM DEPENDENT CALCULATIONS AND OUTPUT 
        CALL POSTEL(LEMFRE,D,IE,IT,NITER,NE,M) 
  10  END DO   
      END SUBROUTINE POST 

      SUBROUTINE POSTEL (LEMFRE,D,IE,IT,NITER,NE,M) 
! * * * * * * * * * * * * * * * * * * * ** * * * * * * * * * * * * * 
!   ELEMENT LEVEL POST SOLUTION CALCULATIONS 
! * * * * * * * * * * * * * * * * * * * ** * * * * * * * * * * * * * 
      DIMENSION D(LEMFRE) 
!
! *** POSTEL PROBLEM DEPENDENT STATEMENTS FOLLOW *** 
!
!  GRADIENT CALCULATION FOR 
! SOLUTION OF LAPLACE EQUATION IN TWO DIMENSIONS 
! USING 8 NODE ISOPARAMETRIC QUADRILATERIAL ELEMENT 
      COMMON/ELARG2/ XY(2),H(8),GRAD(2),AJ(2,2),AJINV(2,2), &
                     DELTA(2,8),GLOBAL(2,8),COL(8) 
! NE = TOTAL NUMBER OF ELEMENTS 
! XY = GLOBAL COORDINATES OF QUADRATURE POINT 
! GRAD = GLOBAL COMPONENTS OF GRADIENT VECTOR 
! H(I) = INTERPOLATION FUNCTION FOR NODE I   1<=I<=N 
! GLOBAL (J, I) = J TH GLOBAL DERIVATIVE OF H(I) 
      DATA KALL, SIZE, DERMAX /1, 0.07,0/ 
 
      IF (  KALL.EQ.0 )     GO TO 10 
      KALL = 0 
      WRITE (6,5000) 
 5000 FORMAT (' *** GLOBAL DERIVATIVES AT INTEGRATION POINTS ***', &
           /, 'POINT       X        Y      DP/DX        DP/DY')
  10  CONTINUE
! BEING ELEMENT POST SOLUTION ANALYSIS 
      WRITE (6,5010)    IE 
 5010 FORMAT (1H0,20X, 'ELEMENT NUMBER ',I3) 
!     READ NUMBER OF POINTS TO BE CALCULATED 
      READ(5,*)  NIP 
      DO 20   J = 1, NIP 
! READ COORDS. AND DERIVATIVE MATRIX 
        WRITE(6,*)  XY, GLOBAL 
! CALCULATE DERIVATIVES, GRAD = GLOBAL* D 
        !b CALL MMULT1(GLOBAL, D, GRAD, NSPACE, LEMFRE) ! fatal error
        CALL MMULT (GLOBAL, D, GRAD, NSPACE, LEMFRE, 1)
! PRINT COORDINATES AND GRADIENT AT THE POINT 
        WRITE (6,5020)      J, XY, GRAD 
 5020   FORMAT (I7,2E12.4,2E15.5) 
      GO TO 20  !b ?? no following statement
! ??? something missing here ? !b
  20 CONTINUE 
! ARE GRADIENT CALCULATIONS COMPLETE FOR ALL ELEMENTS 
      IF ( IE .LT.NE ) RETURN 
      END SUBROUTINE POSTEL 

      SUBROUTINE  RPRINT (A,NR,NC,IOPT)
!     * * * * * * * * * * * * * * * * * * * * * * * * * *
!              PRINTING OF REAL MATRIX A(NR,NC)
!     * * * * * * * * * * * * * * * * * * * * * * * * * *
!DP   IMPLICIT REAL*8 (A)
      DIMENSION  A(*), NCOL(8)
!     IF IOPT=0 USE F FORMAT, OTHERWISE USE E FORMAT
      DO 50  J = 1,NC,8
      JL1 = J - 1
      MAXCOL = 1
      K = NC - JL1
      MAXCOL = MIN0 (K,8)
      MXCLL1 = MAXCOL - 1
      DO 10  L = 1,MAXCOL
      NCOL(L) = L + JL1
   10 END DO
      WRITE (6,5000)  ( NCOL(N),N=1,MAXCOL )
 5000 FORMAT ('  ROW/COL',I10,7I14)
      DO 40  N = 1,NR
      NL = N + (J-1)*NR
      NH = NL + MXCLL1*NR
      IF ( IOPT )  30,20,30
   20 WRITE (6,5010)  N,( A(I),I=NL,NH,NR )
 5010 FORMAT (I6,8F14.4)
      GO TO 40
   30 WRITE (6,5020)  N,( A(I),I=NL,NH,NR )
 5020 FORMAT (I6,8(1PE14.5))
   40 END DO   
   50 END DO   
      END SUBROUTINE RPRINT