Initial commit from Polytechnique Montreal

1 files changed, 621 insertions, 0 deletions

diff --git a/Trivac/src/BIVCOL.f b/Trivac/src/BIVCOL.f
new file mode 100755
index 0000000..733855e
--- /dev/null
+++ b/Trivac/src/BIVCOL.f
@@ -0,0 +1,621 @@
+*DECK BIVCOL
+      SUBROUTINE BIVCOL (IPTRK,IMPX,IELEM,ICOL)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Selection of the unit matrices (mass, stiffness, etc.) for a finite
+* element approximation.
+*
+*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
+* IPTRK   L_TRACK pointer to the tracking information.
+* IMPX    print parameter.
+* IELEM   degree of the finite elements: =1: (linear polynomials);
+*         =2: (parabolic polynomials);   =3: (cubic polynomials);
+*         =4: (quartic polynomials).
+* ICOL    type of quadrature used to integrate the mass matrices:
+*         =1: (analytic integration);  =2: (Gauss-Lobatto quadrature)
+*         =3: (Gauss-Legendre quadrature).
+*         IELEM=1 with ICOL=2 is equivalent to finite differences.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPTRK
+      INTEGER IMPX,IELEM,ICOL
+*----
+*  LOCAL VARIABLES
+*----
+      CHARACTER*40 HTYPE
+      DOUBLE PRECISION DSUM
+      REAL EL(2,2),TL(2),TSL(2),RL(2,2),RSL(2,2),QSL(2,2),TSL1(2),
+     1 TSL2(2),RL2(2,2),RSL2(2,2)
+      REAL RHA6(6,6),QHA6(6,6),RHL6(6,6),QHL6(6,6),RTA(3,3),QTA(3,3),
+     1 RTL(3,3),QTL(3,3)
+      REAL EP(3,3),TP(3),TSP(3),RP(3,3),VP(3,2),HP(3,2),RSP(3,3),
+     1 QSP(3,3),EP1(3,3),TP1(3),TSP1(3),VP1(3,2),HP1(3,2),QSP1(3,3),
+     2 EP2(3,3),TP2(3),TSP2(3),RP2(3,3),VP2(3,2),HP2(3,2),RSP2(3,3),
+     3 QSP2(3,3)
+      REAL EC(4,4),TC(4),TSC(4),RC(4,4),VC(4,3),HC(4,3),RSC(4,4),
+     1 QSC(4,4),EC1(4,4),TC1(4),TSC1(4),VC1(4,3),HC1(4,3),QSC1(4,4),
+     2 EC2(4,4),TC2(4),TSC2(4),RC2(4,4),VC2(4,3),HC2(4,3),RSC2(4,4),
+     3 QSC2(4,4)
+      REAL EQ(5,5),VQ(5,4),HQ(5,4),TQ(5),TSQ(5),QSQ(5,5)
+      REAL RLQ(2,2),RLR(2,2),RL1Q(2,2),RL1R(2,2),RL2Q(2,2),RL2R(2,2)
+*-----------------------------------------------------------------------
+* THE BIVAC REFERENCE ELEMENT IS DEFINED BETWEEN POINTS -1/2 AND +1/2.
+* THE COLLOCATION POLYNOMIALS CORRESPONDING TO APPROXIMATIONS LL2$, PL3,
+* PL3$, PL3#, CL4, CL4$, CL4# AND QL5$ ARE PARTIALLY OR COMPLETELY
+* ORTHONORMALIZED IN ORDER TO PRODUCE A SPARSE MASS MATRIX.
+*-----------------------------------------------------------------------
+*
+******************* FINITE ELEMENT BASIC MATRICES **********************
+*                                                                      *
+      REAL T(5),TS(5),R(25),RS(25),Q(25),QS(25),V(20),H(20),E(25),
+     1 RH(6,6),QH(6,6),RT(3,3),QT(3,3),R1DQ(4),R1DR(4)
+*                                                                      *
+* LC      : NUMBER OF POLYNOMIALS IN A COMPLETE 1-D BASIS.             *
+* T       : CARTESIAN LINEAR PRODUCT VECTOR.                           *
+* TS      : CYLINDRICAL LINEAR PRODUCT VECTOR.                         *
+* R       : CARTESIAN MASS MATRIX.                                     *
+* RS      : CYLINDRICAL MASS MATRIX.                                   *
+* Q       : CARTESIAN STIFFNESS MATRIX.                                *
+* QS      : CYLINDRICAL STIFFNESS MATRIX.                              *
+* V       : NODAL COUPLING MATRIX.                                     *
+* H       : PIOLAT (HEXAGONAL) COUPLING MATRIX.                        *
+* E       : POLYNOMIAL COEFFICIENTS.                                   *
+* RH      : HEXAGONAL MASS MATRIX.                                     *
+* QH      : HEXAGONAL STIFFNESS MATRIX.                                *
+* RT      : TRIANGULAR MASS MATRIX.                                    *
+* QT      : TRIANGULAR STIFFNESS MATRIX.                               *
+* R1DQ    : SPHERICAL MASS MATRIX.                                     *
+* R1DR    : SPHERICAL MASS MATRIX.                                     *
+*                                                                      *
+************************************************************************
+*
+*----
+*  LINEAR LAGRANGIAN POLYNOMIALS.
+*----
+      DATA EL/0.5,-1.0,0.5,1.0/
+      DATA TL/0.5,0.5/
+      DATA RL/
+     $ 0.333333333333, 0.166666666667, 0.166666666667, 0.333333333333/
+*     CYLINDRICAL OPTION MATRICES:
+      DATA TSL/-0.083333333333, 0.083333333333/
+      DATA RSL/
+     $-0.083333333333, 0.000000000000, 0.000000000000, 0.083333333333/
+      DATA QSL/
+     $ 0.000000000000, 0.000000000000, 0.000000000000, 0.000000000000/
+*     SPHERICAL OPTION MATRICES (ANALYTIC INTEGRATION):
+      DATA RLQ/-.083333333333, 0.0, 0.0, 0.083333333333/
+      DATA RLR/
+     $ 0.033333333333, 0.008333333333, 0.008333333333, 0.033333333333/
+*     GAUSS-LOBATTO (FINITE DIFFERENCES) MATRICES:
+      DATA TSL1/-0.25,0.25/
+*     SPHERICAL OPTION MATRICES (GAUSS-LOBATTO):
+      DATA RL1Q/-0.166666666667, 0.0, 0.0, 0.166666666667/
+      DATA RL1R/0.041666666667, 0.0, 0.0, 0.041666666667/
+*     GAUSS-LEGENDRE (SUPERCONVERGENT) MATRICES:
+      DATA TSL2/0.0,0.0/
+      DATA RL2/0.25,0.25,0.25,0.25/
+      DATA RSL2/0.0,0.0,0.0,0.0/
+*     SPHERICAL OPTION MATRICES (SUPERCONVERGENT):
+      DATA RL2Q/
+     $ -0.03472222222,-0.0069444444444,-0.0069444444444, 0.048611111111/
+      DATA RL2R/
+     $  0.00925925926, 0.0185185185185, 0.0185185185185, 0.037037037037/
+*
+*     ANALYTIC INTEGRATION FOR HEXAGON AND TRIANGLE.
+      DATA RHA6/
+     > 0.158470 , 0.086580 , 0.036760 , 0.027870 , 0.036760 , 0.086580,
+     > 0.086580 , 0.158470 , 0.086580 , 0.036760 , 0.027870 , 0.036760,
+     > 0.036760 , 0.086580 , 0.158470 , 0.086580 , 0.036760 , 0.027870,
+     > 0.027870 , 0.036760 , 0.086580 , 0.158470 , 0.086580 , 0.036760,
+     > 0.036760 , 0.027870 , 0.036760 , 0.086580 , 0.158470 , 0.086580,
+     > 0.086580 , 0.036760 , 0.027870 , 0.036760 , 0.086580 , 0.158470/
+      DATA QHA6/
+     > 0.760640 ,-0.161980 ,-0.169310 ,-0.098060 ,-0.169310 ,-0.161980,
+     >-0.161980 , 0.760640 ,-0.161980 ,-0.169310 ,-0.098060 ,-0.169310,
+     >-0.169310 ,-0.161980 , 0.760640 ,-0.161980 ,-0.169310 ,-0.098060,
+     >-0.098060 ,-0.169310 ,-0.161980 , 0.760640 ,-0.161980 ,-0.169310,
+     >-0.169310 ,-0.098060 ,-0.169310 ,-0.161980 , 0.760640 ,-0.161980,
+     >-0.161980 ,-0.169310 ,-0.098060 ,-0.169310 ,-0.161980 , 0.760640/
+      DATA RTA/
+     > 1.0, 0.5, 0.5, 0.5, 1.0, 0.5, 0.5, 0.5, 1.0/
+      DATA QTA/
+     > 1.0,-0.5,-0.5,-0.5, 1.0,-0.5,-0.5,-0.5, 1.0/
+*
+*     GAUSS-LOBATTO INTEGRATION FOR HEXAGON AND TRIANGLE.
+      DATA RHL6/
+     > 1.000000 , 0.000000 , 0.000000 , 0.000000 , 0.000000 , 0.000000,
+     > 0.000000 , 1.000000 , 0.000000 , 0.000000 , 0.000000 , 0.000000,
+     > 0.000000 , 0.000000 , 1.000000 , 0.000000 , 0.000000 , 0.000000,
+     > 0.000000 , 0.000000 , 0.000000 , 1.000000 , 0.000000 , 0.000000,
+     > 0.000000 , 0.000000 , 0.000000 , 0.000000 , 1.000000 , 0.000000,
+     > 0.000000 , 0.000000 , 0.000000 , 0.000000 , 0.000000 , 1.000000/
+      DATA QHL6/
+     > 1.666667 ,-1.000000 , 0.166667 ,-1.000000 , 0.166667 , 0.000000,
+     >-1.000000 , 1.666667 ,-1.000000 , 0.166667 , 0.000000 , 0.166667,
+     > 0.166667 ,-1.000000 , 1.666667 , 0.000000 , 0.166667 ,-1.000000,
+     >-1.000000 , 0.166667 , 0.000000 , 1.666667 ,-1.000000 , 0.166667,
+     > 0.166667 , 0.000000 , 0.166667 ,-1.000000 , 1.666667 ,-1.000000,
+     > 0.000000 , 0.166667 ,-1.000000 , 0.166667 ,-1.000000 , 1.666667/
+      DATA RTL/
+     > 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0/
+      DATA QTL/
+     > 1.0,-0.5,-0.5,-0.5, 1.0,-0.5,-0.5,-0.5, 1.0/
+*----
+*  PARABOLIC LAGRANGIAN POLYNOMIALS.
+*----
+      DATA EP/
+     $-0.125000000000,-1.000000000000, 2.500000000000,
+     $ 1.250000000000, 0.000000000000,-5.000000000000,
+     $-0.125000000000, 1.000000000000, 2.500000000000/
+      DATA TP/
+     $ 0.083333333333, 0.833333333333, 0.083333333333/
+      DATA RP/
+     $ 0.125000000000, 0.000000000000,-0.041666666667,
+     $ 0.000000000000, 0.833333333333, 0.000000000000,
+     $-0.041666666667, 0.000000000000, 0.125000000000/
+      DATA VP/
+     $-1.000000000000, 0.000000000000, 1.000000000000,
+     $ 1.443375672974,-2.886751345948, 1.443375672974/
+      DATA HP/
+     $ 0.083333333333, 0.833333333333, 0.083333333333,
+     $-0.288675134595, 0.000000000000, 0.288675134595/
+*     CYLINDRICAL OPTION MATRICES:
+      DATA TSP/
+     $-0.083333333333, 0.000000000000, 0.083333333333/
+      DATA RSP/
+     $-0.041666666667,-0.041666666667, 0.000000000000,
+     $-0.041666666667, 0.000000000000, 0.041666666667,
+     $ 0.000000000000, 0.041666666667, 0.041666666667/
+      DATA QSP/
+     $-0.833333333333, 0.833333333333, 0.000000000000,
+     $ 0.833333333333, 0.000000000000,-0.833333333333,
+     $ 0.000000000000,-0.833333333333, 0.833333333333/
+*     GAUSS-LOBATTO (VARIATIONAL COLLOCATION METHOD) MATRICES:
+      DATA EP1/0.0,-1.0,2.0,1.0,0.0,-4.0,0.0,1.0,2.0/
+      DATA TP1/
+     $ 0.166666666667, 0.666666666667, 0.166666666667/
+      DATA VP1/
+     $-1.000000000000, 0.0000000000000, 1.000000000000,
+     $ 1.154700538379,-2.3094010767585, 1.154700538379/
+      DATA HP1/
+     $ 0.166666666667, 0.666666666667, 0.166666666667,
+     $-0.288675134595, 0.000000000000, 0.288675134595/
+*     CYLINDRICAL GAUSS-LOBATTO (VARIATIONAL COLLOCATION METHOD)
+*     MATRICES.
+      DATA TSP1/
+     $-0.083333333333, 0.000000000000, 0.083333333333/
+      DATA QSP1/
+     $-0.666666666667, 0.666666666667, 0.000000000000,
+     $ 0.666666666667, 0.000000000000,-0.666666666667,
+     $ 0.000000000000,-0.666666666667, 0.666666666667/
+*     GAUSS-LEGENDRE (SUPERCONVERGENT) MATRICES:
+      DATA EP2/
+     $-0.250000000000,-1.000000000000, 3.000000000000,
+     $ 1.500000000000, 0.000000000000,-6.000000000000,
+     $-0.250000000000, 1.000000000000, 3.000000000000/
+      DATA TP2/
+     $ 0.000000000000, 1.000000000000, 0.000000000000/
+      DATA RP2/
+     $ 0.083333333333, 0.000000000000,-0.083333333333,
+     $ 0.000000000000, 1.000000000000, 0.000000000000,
+     $-0.083333333333, 0.000000000000, 0.083333333333/
+      DATA VP2/
+     $-1.000000000000, 0.000000000000, 1.000000000000,
+     $ 1.732050807569,-3.464101615138, 1.732050807569/
+      DATA HP2/
+     $ 0.000000000000, 1.000000000000, 0.000000000000,
+     $-0.288675134595, 0.000000000000, 0.288675134595/
+*     CYLINDRICAL GAUSS-LEGENDRE (SUPERCONVERGENT) MATRICES:
+      DATA TSP2/
+     $-0.083333333333, 0.000000000000, 0.083333333333/
+      DATA RSP2/
+     $ 0.000000000000,-0.083333333333, 0.000000000000,
+     $-0.083333333333, 0.000000000000, 0.083333333333,
+     $ 0.000000000000, 0.083333333333, 0.000000000000/
+      DATA QSP2/
+     $-1.000000000000, 1.000000000000, 0.000000000000,
+     $ 1.000000000000, 0.000000000000,-1.000000000000,
+     $ 0.000000000000,-1.000000000000, 1.000000000000/
+*----
+*  CUBIC LAGRANGIAN POLYNOMIALS.
+*----
+      DATA EC/
+     $-0.125000000000, 0.750000000000, 2.500000000000, -7.000000000000,
+     $ 0.625000000000,-3.307189138831,-2.500000000000, 13.228756555323,
+     $ 0.625000000000, 3.307189138831,-2.500000000000,-13.228756555323,
+     $-0.125000000000,-0.750000000000, 2.500000000000,  7.000000000000/
+      DATA TC/
+     $ 0.083333333333, 0.416666666667, 0.416666666667, 0.083333333333/
+      DATA RC/
+     $ 0.066666666667, 0.000000000000, 0.000000000000, 0.016666666667,
+     $ 0.000000000000, 0.416666666667, 0.000000000000, 0.000000000000,
+     $ 0.000000000000, 0.000000000000, 0.416666666667, 0.000000000000,
+     $ 0.016666666667, 0.000000000000, 0.000000000000, 0.066666666667/
+      DATA VC/
+     $-1.000000000000, 0.000000000000, 0.000000000000,1.000000000000,
+     $ 1.443375672974,-1.443375672974,-1.443375672974,1.443375672974,
+     $-1.565247584250, 2.958039891550,-2.958039891550,1.565247584250/
+      DATA HC/
+     $ 0.083333333333, 0.416666666667, 0.416666666667, 0.083333333333,
+     $-0.086602540378,-0.381881307913, 0.381881307913, 0.086602540378,
+     $ 0.186338998125,-0.186338998125,-0.186338998125, 0.186338998125/
+*     CYLINDRICAL OPTION MATRICES:
+      DATA TSC /
+     $-0.025000000000,-0.110239637961, 0.110239637961, 0.025000000000/
+      DATA RSC/
+     $-0.025000000000,-0.015748519709, 0.015748519709, 0.000000000000,
+     $-0.015748519709,-0.078742598544, 0.000000000000,-0.015748519709,
+     $ 0.015748519709, 0.000000000000, 0.078742598544, 0.015748519709,
+     $ 0.000000000000,-0.015748519709, 0.015748519709, 0.025000000000/
+      DATA QSC/
+     $-2.000000000000, 2.102396379610,-0.102396379610, 0.000000000000,
+     $ 2.102396379610,-2.204792759220, 0.000000000000, 0.102396379610,
+     $-0.102396379610, 0.000000000000, 2.204792759220,-2.102396379610,
+     $ 0.000000000000, 0.102396379610,-2.102396379610, 2.000000000000/
+*     GAUSS-LOBATTO (VARIATIONAL COLLOCATION METHOD) MATRICES:
+      DATA EC1/
+     $-0.125000000000, 0.250000000000, 2.500000000000, -5.000000000000,
+     $ 0.625000000000,-2.795084971875,-2.500000000000, 11.180339887499,
+     $ 0.625000000000, 2.795084971875,-2.500000000000,-11.180339887499,
+     $-0.125000000000,-0.250000000000, 2.500000000000,  5.000000000000/
+      DATA TC1/
+     $ 0.083333333333, 0.416666666667, 0.416666666667, 0.083333333333/
+      DATA VC1/
+     $-1.000000000000, 0.000000000000, 0.000000000000,1.000000000000,
+     $ 1.443375672974,-1.443375672974,-1.443375672974,1.443375672974,
+     $-1.118033988750, 2.500000000000,-2.500000000000,1.118033988750/
+      DATA HC1/
+     $ 0.083333333333, 0.416666666667, 0.416666666667, 0.083333333333,
+     $-0.144337567297,-0.322748612184, 0.322748612184, 0.144337567297,
+     $ 0.186338998125,-0.186338998125,-0.186338998125, 0.186338998125/
+*     CYLINDRICAL GAUSS-LOBATTO (VARIATIONAL COLLOCATION METHOD)
+*     MATRICES:
+      DATA TSC1/
+     $-0.041666666667,-0.093169499062, 0.093169499062, 0.041666666667/
+      DATA QSC1/
+     $-1.666666666667, 1.765028323958,-0.098361657292, 0.000000000000,
+     $ 1.765028323958,-1.863389981250, 0.000000000000, 0.098361657292,
+     $-0.098361657292, 0.000000000000, 1.863389981250,-1.765028323958,
+     $ 0.000000000000, 0.098361657292,-1.765028323958, 1.666666666667/
+*     GAUSS-LEGENDRE (SUPERCONVERGENT) MATRICES:
+      DATA EC2/
+     $-0.125000000000, 1.500000000000, 2.500000000000,-10.000000000000,
+     $ 0.625000000000,-3.952847075210,-2.500000000000, 15.811388300842,
+     $ 0.625000000000, 3.952847075210,-2.500000000000,-15.811388300842,
+     $-0.125000000000,-1.500000000000, 2.500000000000, 10.000000000000/
+      DATA TC2/
+     $ 0.083333333333, 0.416666666667, 0.416666666667, 0.083333333333/
+      DATA RC2/
+     $ 0.041666666667, 0.000000000000, 0.000000000000, 0.041666666667,
+     $ 0.000000000000, 0.416666666667, 0.000000000000, 0.000000000000,
+     $ 0.000000000000, 0.000000000000, 0.416666666667, 0.000000000000,
+     $ 0.041666666667, 0.000000000000, 0.000000000000, 0.041666666667/
+      DATA VC2/
+     $-1.000000000000, 0.000000000000, 0.000000000000,1.000000000000,
+     $ 1.443375672974,-1.443375672974,-1.443375672974,1.443375672974,
+     $-2.236067977500, 3.535533905933,-3.535533905933,2.236067977500/
+      DATA HC2/
+     $ 0.083333333333, 0.416666666667, 0.416666666667, 0.083333333333,
+     $ 0.000000000000,-0.456435464588, 0.456435464588, 0.000000000000,
+     $ 0.186338998125,-0.186338998125,-0.186338998125, 0.186338998125/
+*     CYLINDRICAL GAUSS-LEGENDRE (SUPERCONVERGENT) MATRICES:
+      DATA TSC2/
+     $ 0.000000000000,-0.131761569174, 0.131761569174, 0.000000000000/
+      DATA RSC2/
+     $ 0.000000000000,-0.032940392293, 0.032940392293, 0.000000000000,
+     $-0.032940392293,-0.065880784587, 0.000000000000,-0.032940392293,
+     $ 0.032940392293, 0.000000000000, 0.065880784587, 0.032940392293,
+     $ 0.000000000000,-0.032940392293, 0.032940392293, 0.000000000000/
+      DATA QSC2/
+     $-2.500000000000, 2.567615691737,-0.067615691737, 0.000000000000,
+     $ 2.567615691737,-2.635231383474, 0.000000000000, 0.067615691737,
+     $-0.067615691737, 0.000000000000, 2.635231383474,-2.567615691737,
+     $ 0.000000000000, 0.067615691737,-2.567615691737, 2.500000000000/
+*----
+*  QUARTIC LAGRANGIAN POLYNOMIALS.
+*----
+      DATA EQ/
+     $  0.000000000000, 0.750000000000,-1.500000000000,-7.000000000000,
+     $ 14.000000000000, 0.000000000000,-2.673169155391, 8.166666666667,
+     $ 10.692676621564,-32.666666666667, 1.000000000000, 0.000000000000,
+     $-13.333333333333, 0.000000000000,37.333333333333, 0.000000000000,
+     $  2.673169155391,8.166666666667,-10.692676621564,-32.666666666667,
+     $  0.000000000000,-0.750000000000,-1.500000000000, 7.000000000000,
+     $ 14.000000000000/
+      DATA TQ/
+     $ 0.050000000000, 0.272222222222, 0.355555555556, 0.272222222222,
+     $ 0.050000000000/
+      DATA VQ/
+     $-1.000000000000, 0.00000000000, 0.000000000000, 0.00000000000,
+     $ 1.000000000000, 1.55884572681,-0.943005439677,-1.23168057427,
+     $-0.943005439677, 1.55884572681,-1.565247584250, 2.39095517873,
+     $ 0.000000000000,-2.39095517873, 1.565247584250, 1.058300524426,
+     $-2.46936789032 , 2.8221347318 ,-2.46936789032 , 1.058300524426/
+      DATA HQ/
+     $ 0.050000000000, 0.272222222222, 0.355555555556, 0.272222222222,
+     $ 0.050000000000,-0.086602540378,-0.308670986291, 0.000000000000,
+     $ 0.308670986291, 0.086602540378, 0.111803398875, 0.086958199125,
+     $-0.397523196000, 0.086958199125, 0.111803398875,-0.132287565553,
+     $ 0.202072594216, 0.000000000000,-0.202072594216, 0.132287565553/
+*     CYLINDRICAL GAUSS-LOBATTO (VARIATIONAL COLLOCATION METHOD)
+*     MATRICES:
+      DATA TSQ/
+     $-0.025000000000,-0.089105638513, 0.000000000000, 0.089105638513,
+     $ 0.025000000000/
+      DATA QSQ/
+     $-3.000000000000, 3.237234826568,-0.266666666667, 0.029431840099,
+     $ 0.000000000000, 3.237234826568,-4.158263130608, 0.950460144139,
+     $ 0.000000000000,-0.029431840099,-0.266666666667, 0.950460144139,
+     $ 0.000000000000,-0.950460144139, 0.266666666667, 0.029431840099,
+     $ 0.000000000000,-0.950460144139, 4.158263130608,-3.237234826568,
+     $ 0.000000000000,-0.029431840099, 0.266666666667,-3.237234826568,
+     $ 3.000000000000/
+*
+      LC=IELEM+1
+      IF((IELEM.EQ.1).AND.(ICOL.EQ.1)) THEN
+*        LL2
+         HTYPE='LINEAR LAGRANGIAN POLYNOMIALS'
+         DO 20 I=1,LC
+         T(I)=TL(I)
+         TS(I)=TSL(I)
+         DO 10 J=1,LC
+         R((J-1)*LC+I)=RL(I,J)
+         RS((J-1)*LC+I)=RSL(I,J)
+         QS((J-1)*LC+I)=QSL(I,J)
+         R1DQ((J-1)*LC+I)=RLQ(I,J)
+         R1DR((J-1)*LC+I)=RLR(I,J)
+         E((J-1)*LC+I)=EL(I,J)
+   10    CONTINUE
+   20    CONTINUE
+         V(1)=-1.0
+         V(2)=1.0
+         H(1)=0.5
+         H(2)=0.5
+         DO 40 I=1,6
+         DO 30 J=1,6
+         RH(I,J)=RHA6(I,J)
+         QH(I,J)=QHA6(I,J)
+   30    CONTINUE
+   40    CONTINUE
+         DO 60 I=1,3
+         DO 50 J=1,3
+         RT(I,J)=RTA(I,J)*SQRT(3.0)/24.0
+         QT(I,J)=QTA(I,J)/SQRT(3.0)
+   50    CONTINUE
+   60    CONTINUE
+      ELSE IF((IELEM.EQ.1).AND.(ICOL.EQ.2)) THEN
+*        LL2$
+         HTYPE='FINITE DIFFERENCES'
+         DO 80 I=1,LC
+         T(I)=TL(I)
+         TS(I)=TSL1(I)
+         DO 70 J=1,LC
+         R((J-1)*LC+I)=0.0
+         RS((J-1)*LC+I)=0.0
+         QS((J-1)*LC+I)=QSL(I,J)
+         R1DQ((J-1)*LC+I)=RL1Q(I,J)
+         R1DR((J-1)*LC+I)=RL1R(I,J)
+         E((J-1)*LC+I)=EL(I,J)
+   70    CONTINUE
+         R((I-1)*LC+I)=TL(I)
+         RS((I-1)*LC+I)=TSL1(I)
+   80    CONTINUE
+         V(1)=-1.0
+         V(2)=1.0
+         H(1)=0.5
+         H(2)=0.5
+         DO 100 I=1,6
+         DO 90 J=1,6
+         RH(I,J)=RHL6(I,J)*SQRT(3.0)/4.0
+         QH(I,J)=QHL6(I,J)*SQRT(3.0)
+   90    CONTINUE
+  100    CONTINUE
+         DO 120 I=1,3
+         DO 110 J=1,3
+         RT(I,J)=RTL(I,J)*SQRT(3.0)/12.0
+         QT(I,J)=QTL(I,J)/SQRT(3.0)
+  110    CONTINUE
+  120    CONTINUE
+      ELSE IF((IELEM.EQ.1).AND.(ICOL.EQ.3)) THEN
+*        LL2#
+         HTYPE='SUPERCONVERGENT LINEAR POLYNOMIALS'
+         DO 140 I=1,LC
+         T(I)=TL(I)
+         TS(I)=TSL2(I)
+         DO 130 J=1,LC
+         R((J-1)*LC+I)=RL2(I,J)
+         RS((J-1)*LC+I)=RSL2(I,J)
+         QS((J-1)*LC+I)=QSL(I,J)
+         R1DQ((J-1)*LC+I)=RL2Q(I,J)
+         R1DR((J-1)*LC+I)=RL2R(I,J)
+         E((J-1)*LC+I)=EL(I,J)
+  130    CONTINUE
+  140    CONTINUE
+         V(1)=-1.0
+         V(2)=1.0
+         H(1)=0.5
+         H(2)=0.5
+      ELSE IF((IELEM.EQ.2).AND.(ICOL.EQ.1)) THEN
+*        PL3
+         HTYPE='PARABOLIC LAGRANGIAN POLYNOMIALS'
+         DO 170 I=1,LC
+         T(I)=TP(I)
+         TS(I)=TSP(I)
+         DO 150 J=1,LC-1
+         V((J-1)*LC+I)=VP(I,J)
+         H((J-1)*LC+I)=HP(I,J)
+  150    CONTINUE
+         DO 160 J=1,LC
+         R((J-1)*LC+I)=RP(I,J)
+         RS((J-1)*LC+I)=RSP(I,J)
+         QS((J-1)*LC+I)=QSP(I,J)
+         E((J-1)*LC+I)=EP(I,J)
+  160    CONTINUE
+  170    CONTINUE
+      ELSE IF((IELEM.EQ.2).AND.(ICOL.EQ.2)) THEN
+*        PL3$
+         HTYPE='PARABOLIC COLLOCATION METHOD'
+         DO 200 I=1,LC
+         T(I)=TP1(I)
+         TS(I)=TSP1(I)
+         DO 180 J=1,LC-1
+         V((J-1)*LC+I)=VP1(I,J)
+         H((J-1)*LC+I)=HP1(I,J)
+  180    CONTINUE
+         DO 190 J=1,LC
+         R((J-1)*LC+I)=0.0
+         RS((J-1)*LC+I)=0.0
+         QS((J-1)*LC+I)=QSP1(I,J)
+         E((J-1)*LC+I)=EP1(I,J)
+  190    CONTINUE
+         R((I-1)*LC+I)=TP1(I)
+         RS((I-1)*LC+I)=TSP1(I)
+  200    CONTINUE
+      ELSE IF((IELEM.EQ.2).AND.(ICOL.EQ.3)) THEN
+*        PL3#
+         HTYPE='PARABOLIC SUPERCONVERGENT POLYNOMIALS'
+         DO 230 I=1,LC
+         T(I)=TP2(I)
+         TS(I)=TSP2(I)
+         DO 210 J=1,LC-1
+         V((J-1)*LC+I)=VP2(I,J)
+         H((J-1)*LC+I)=HP2(I,J)
+  210    CONTINUE
+         DO 220 J=1,LC
+         R((J-1)*LC+I)=RP2(I,J)
+         RS((J-1)*LC+I)=RSP2(I,J)
+         QS((J-1)*LC+I)=QSP2(I,J)
+         E((J-1)*LC+I)=EP2(I,J)
+  220    CONTINUE
+  230    CONTINUE
+      ELSE IF((IELEM.EQ.3).AND.(ICOL.EQ.1)) THEN
+*        CL4
+         HTYPE='CUBIC LAGRANGIAN POLYNOMIALS'
+         DO 260 I=1,LC
+         T(I)=TC(I)
+         TS(I)=TSC(I)
+         DO 240 J=1,LC-1
+         V((J-1)*LC+I)=VC(I,J)
+         H((J-1)*LC+I)=HC(I,J)
+  240    CONTINUE
+         DO 250 J=1,LC
+         R((J-1)*LC+I)=RC(I,J)
+         RS((J-1)*LC+I)=RSC(I,J)
+         QS((J-1)*LC+I)=QSC(I,J)
+         E((J-1)*LC+I)=EC(I,J)
+  250    CONTINUE
+  260    CONTINUE
+      ELSE IF((IELEM.EQ.3).AND.(ICOL.EQ.2)) THEN
+*        CL4$
+         HTYPE='CUBIC COLLOCATION METHOD'
+         DO 290 I=1,LC
+         T(I)=TC1(I)
+         TS(I)=TSC1(I)
+         DO 270 J=1,LC-1
+         V((J-1)*LC+I)=VC1(I,J)
+         H((J-1)*LC+I)=HC1(I,J)
+  270    CONTINUE
+         DO 280 J=1,LC
+         R((J-1)*LC+I)=0.0
+         RS((J-1)*LC+I)=0.0
+         QS((J-1)*LC+I)=QSC1(I,J)
+         E((J-1)*LC+I)=EC1(I,J)
+  280    CONTINUE
+         R((I-1)*LC+I)=TC1(I)
+         RS((I-1)*LC+I)=TSC1(I)
+  290    CONTINUE
+      ELSE IF((IELEM.EQ.3).AND.(ICOL.EQ.3)) THEN
+*        CL4#
+         HTYPE='SUPERCONVERGENT CUBIC POLYNOMIALS'
+         DO 320 I=1,LC
+         T(I)=TC2(I)
+         TS(I)=TSC2(I)
+         DO 300 J=1,LC-1
+         V((J-1)*LC+I)=VC2(I,J)
+         H((J-1)*LC+I)=HC2(I,J)
+  300    CONTINUE
+         DO 310 J=1,LC
+         R((J-1)*LC+I)=RC2(I,J)
+         RS((J-1)*LC+I)=RSC2(I,J)
+         QS((J-1)*LC+I)=QSC2(I,J)
+         E((J-1)*LC+I)=EC2(I,J)
+  310    CONTINUE
+  320    CONTINUE
+      ELSE IF((IELEM.EQ.4).AND.(ICOL.EQ.2)) THEN
+*        QL5$
+         HTYPE='QUARTIC COLLOCATION METHOD'
+         DO 350 I=1,LC
+         T(I)=TQ(I)
+         TS(I)=TSQ(I)
+         DO 330 J=1,LC-1
+         V((J-1)*LC+I)=VQ(I,J)
+         H((J-1)*LC+I)=HQ(I,J)
+  330    CONTINUE
+         DO 340 J=1,LC
+         R((J-1)*LC+I)=0.0
+         RS((J-1)*LC+I)=0.0
+         QS((J-1)*LC+I)=QSQ(I,J)
+         E((J-1)*LC+I)=EQ(I,J)
+  340    CONTINUE
+         R((I-1)*LC+I)=TQ(I)
+         RS((I-1)*LC+I)=TSQ(I)
+  350    CONTINUE
+      ELSE
+         CALL XABORT('BIVCOL: TYPE OF FINITE ELEMENT NOT AVAILABLE.')
+      ENDIF
+*----
+*  COMPUTE THE CARTESIAN STIFFNESS MATRIX FROM THE TENSORIAL PRODUCT OF
+*  TWO NODAL COUPLING MATRICES.
+*----
+      DO 380 I=1,LC
+      DO 370 J=1,LC
+      DSUM=0.0D0
+      DO 360 K=1,LC-1
+      DSUM=DSUM+V((K-1)*LC+I)*V((K-1)*LC+J)
+  360 CONTINUE
+      Q((J-1)*LC+I)=REAL(DSUM)
+  370 CONTINUE
+  380 CONTINUE
+      IF(IMPX.GT.0) WRITE (6,'(/9H BIVCOL: ,A40)') HTYPE
+*----
+*  SAVE THE UNIT MATRICES ON LCM.
+*----
+      CALL LCMSIX(IPTRK,'BIVCOL',1)
+      CALL LCMPUT(IPTRK,'T',LC,2,T)
+      CALL LCMPUT(IPTRK,'TS',LC,2,TS)
+      CALL LCMPUT(IPTRK,'R',LC*LC,2,R)
+      CALL LCMPUT(IPTRK,'RS',LC*LC,2,RS)
+      IF(IELEM.EQ.1) THEN
+         CALL LCMPUT(IPTRK,'RSH1',LC*LC,2,R1DQ)
+         CALL LCMPUT(IPTRK,'RSH2',LC*LC,2,R1DR)
+      ENDIF
+      CALL LCMPUT(IPTRK,'Q',LC*LC,2,Q)
+      CALL LCMPUT(IPTRK,'QS',LC*LC,2,QS)
+      CALL LCMPUT(IPTRK,'V',LC*(LC-1),2,V)
+      CALL LCMPUT(IPTRK,'H',LC*(LC-1),2,H)
+      CALL LCMPUT(IPTRK,'E',LC*LC,2,E)
+      IF((IELEM.EQ.1).AND.(ICOL.LE.2)) THEN
+         CALL LCMPUT(IPTRK,'RH',36,2,RH)
+         CALL LCMPUT(IPTRK,'QH',36,2,QH)
+         CALL LCMPUT(IPTRK,'RT',9,2,RT)
+         CALL LCMPUT(IPTRK,'QT',9,2,QT)
+      ENDIF
+      CALL LCMSIX(IPTRK,' ',2)
+      RETURN
+      END