C     ----------------------------------------------
      SUBROUTINE WSHAPE1(KXX,KYX,KYY,KZX,KZY,KZZ,SS,
     +    XE,YE,ZE,AE,BE,VOL,L2)
C     ----------------------------------------------
C
C     SUBROUTINE AUTHOR:
C     TIM DURBIN
C     TIMOTHY J. DURBIN, INC.
C     5330 PRIMROSE DRIVE, SUITE 228
C     FAIR OAKS, CA 95628
C     (916) 536-2314
C
C     PURPOSE:
C     THIS SUBROUTINE PERFORMS NECESSARY INTEGRATIONS THAT DETERMINE THE
C     VALUES OF THE COEFFICIENTS IN MATRICES A AND B AND VECTOR F.
C     THE FIRST ENTRY POINT IS CALLED ONCE BY THE MAIN PROGRAM. THE
C     SECOND ENTRY POINT IS CALLED ONCE FOR EACH TIMESTEP ITERATION. THE
C     THIRD ENTRY POINT IS CALLED ONCE FOR EACH TIMESTEP.
C
C
      INCLUDE 'DIMEN.INC'
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
C
C     INTERNAL SCALARS AND ARRAYS.
C     KXX         COMPONENT OF HYDRAULIC CONDUCTIVITY TENSOR.
C     KXY         COMPONENT OF HYDRAULIC CONDUCTIVITY TENSOR.
C     KXZ         COMPONENT OF HYDRAULIC CONDUCTIVITY TENSOR.
C     KYX         COMPONENT OF HYDRAULIC CONDUCTIVITY TENSOR.
C     KYY         COMPONENT OF HYDRAULIC CONDUCTIVITY TENSOR.
C     KYZ         COMPONENT OF HYDRAULIC CONDUCTIVITY TENSOR.
C     KZX         COMPONENT OF HYDRAULIC CONDUCTIVITY TENSOR.
C     KZY         COMPONENT OF HYDRAULIC CONDUCTIVITY TENSOR.
C     KZZ         COMPONENT OF HYDRAULIC CONDUCTIVITY TENSOR.
C     LEAK        FAULT LEAKANCE.
C     XE()        ELEMENT COORDINATES.
C     YE()        ELEMENT COORDINATES.
C     ZE()        ELEMENT COORDINATES.
C     AE()        ELEMENT CONDUCTANCE MATRIX.
C     BE()        ELEMENT STORAGE MATRIX.
C     HE()        ELEMENT HEADS.
C     VOL         ELEMENT VOLUME.
C     AREA        ELEMENT AREA.
C     B()         BASIS FUNCTION COEFICIENT.
C     C()         BASIS FUNCTION COEFICIENT.
C     D()         BASIS FUNCTION COEFICIENT.
C
C     INTERNAL ARRAYS.
      REAL*8 KXX,KXY,KXZ,KYX,KYY,KYZ,KZX,KZY,KZZ,LEAK
      REAL*8 XE(4),YE(4),ZE(4),AE(4,4),BE(4,4),HE(4)
      REAL*8 VOL,AREA
      REAL*8 B(4),C(4),D(4)
      SAVE
C
C     TETRAHEDRAL ELEMENTS.
C     REASSIGN NODAL COORDINATES.
      X1=XE(1)
      X2=XE(2)
      X3=XE(3)
      X4=XE(4)
      Y1=YE(1)
      Y2=YE(2)
      Y3=YE(3)
      Y4=YE(4)
      Z1=ZE(1)
      Z2=ZE(2)
      Z3=ZE(3)
      Z4=ZE(4)
C
C     COFACTORS OF THE VOLUME DETERMINATE.
      A1=X2*(Y3*Z4-Y4*Z3)-X3*(Y2*Z4-Y4*Z2)+X4*(Y2*Z3-Y3*Z2)
      B(1)=-Y3*Z4+Y4*Z3+Y2*Z4-Y4*Z2-Y2*Z3+Y3*Z2
      B(2)=+Y3*Z4-Y4*Z3-Y1*Z4+Y4*Z1+Y1*Z3-Y3*Z1
      B(3)=-Y2*Z4+Y4*Z2+Y1*Z4-Y4*Z1-Y1*Z2+Y2*Z1
      B(4)=+Y2*Z3-Y3*Z2-Y1*Z3+Y3*Z1+Y1*Z2-Y2*Z1
      C(1)=X3*Z4-X4*Z3-X2*Z4+X4*Z2+X2*Z3-X3*Z2
      C(2)=-X3*Z4+X4*Z3+X1*Z4-X4*Z1-X1*Z3+X3*Z1
      C(3)=X2*Z4-X4*Z2-X1*Z4+X4*Z1+X1*Z2-X2*Z1
      C(4)=-X2*Z3+X3*Z2+X1*Z3-X3*Z1-X1*Z2+X2*Z1
      D(1)=-X3*Y4+X4*Y3+X2*Y4-X4*Y2-X2*Y3+X3*Y2
      D(2)=+X3*Y4-X4*Y3-X1*Y4+X4*Y1+X1*Y3-X3*Y1
      D(3)=-X2*Y4+X4*Y2+X1*Y4-X4*Y1-X1*Y2+X2*Y1
      D(4)=+X2*Y3-X3*Y2-X1*Y3+X3*Y1+X1*Y2-X2*Y1
C
C     VOLUME FOR TETRAHEDRON.
      VOL=(A1+B(1)*X1+C(1)*Y1+D(1)*Z1)/6.0
      IF (VOL.LE.0.0) THEN
          WRITE(*,902) L2
          WRITE(*,904)
          STOP
      ENDIF
C
C     ELEMENT STIFFNESS MATRIX (AE) AND DYNAMIC MATRIX (BE).
      KXY=KYX
      KXZ=KZX
      KYZ=KZY
      DO I=1,4
          DO J=1,4
              AE(I,J)=(KXX*B(I)*B(J)+KXY*B(I)*C(J)+KXZ*B(I)*D(J)+
     +            KYX*C(I)*B(J)+KYY*C(I)*C(J)+KYZ*C(I)*D(J)+
     +            KZX*D(I)*B(J)+KZY*D(I)*C(J)+KZZ*D(I)*D(J))/
     +            (36.0*VOL)
              IF (I.NE.J) THEN
                  BE(I,J)=SS*VOL/20.0
              ELSE
                  BE(I,J)=SS*VOL/10.0
              ENDIF
          END DO
      END DO
      RETURN
C
C     FORMATS:
902   FORMAT(1X,'*ERROR* NEGATIVE OR ZERO VOLUME FOR TETRAHEDRON',I10)
904   FORMAT(1X,'*ERROR* SIMULATION TERMINATED FOR ERROR IN WSHAPE')
C
C     ---------------------------------------------
      ENTRY WSHAPE2(TRAN,STOR,XE,YE,AE,BE,AREA,IFE)
C     ---------------------------------------------
C
C     TRIANGULAR ELEMENTS.
C     REASSIGN NODAL COORDINATES.
      X1=XE(1)
      X2=XE(2)
      X3=XE(3)
      Y1=YE(1)
      Y2=YE(2)
      Y3=YE(3)
C
C     AREA OF PRISMATIC ELEMENT
      B(1)=Y2-Y3
      B(2)=Y3-Y1
      B(3)=Y1-Y2
      C(1)=X3-X2
      C(2)=X1-X3
      C(3)=X2-X1
      AREA=(X2*Y3-X3*Y2-X1*Y3+X3*Y1+X1*Y2-X2*Y1)/2.0
      IF (AREA.LE.0) THEN
          WRITE(*,906) L2
          WRITE(*,904)
          STOP
      ENDIF
C
C     ELEMENT STIFFNESS MATRIX (AE) AND DYNAMIC MATRIX (BE).
      DO I=1,3
          DO J=1,3
              AE(I,J)=TRAN*(B(I)*B(J)+C(I)*C(J))/
     +                (4.0*AREA)
              BE(I,J)=STOR*AREA/12.0
              IF (I.EQ.J) BE(I,J)=BE(I,J)*2.0
          END DO
      END DO
      RETURN
C
C     FORMATS:
906   FORMAT(1X,'*ERROR* NEGATIVE OR ZERO AREA FOR ELEMENT',I10)
C
C     ----------------------
      ENTRY WSHAPE3(LEAK,AE)
C     ----------------------
C
C     LINK ELEMENTS.
      DO I=1,2
          DO J=1,3
              IF (I.NE.J) THEN
                  AE(I,J)=-LEAK
              ELSE
                  AE(I,J)=LEAK
              ENDIF
          END DO
      END DO
      RETURN
C
C     -----------------------------------------------------------
      ENTRY WSHAPE4(XE,YE,ZE,HE,KXX,KYX,KYY,KZX,KZY,KZZ,QX,QY,QZ)
C     -----------------------------------------------------------
C
C     GROUNDWATER DARCY VELOCITIES.
C     REASSIGN NODAL COORDINATES.
      X1=XE(1)
      X2=XE(2)
      X3=XE(3)
      X4=XE(4)
      Y1=YE(1)
      Y2=YE(2)
      Y3=YE(3)
      Y4=YE(4)
      Z1=ZE(1)
      Z2=ZE(2)
      Z3=ZE(3)
      Z4=ZE(4)
C
C     COFACTORS OF THE VOLUME DETERMINATE.
      A1=X2*(Y3*Z4-Y4*Z3)-X3*(Y2*Z4-Y4*Z2)+X4*(Y2*Z3-Y3*Z2)
      B(1)=-Y3*Z4+Y4*Z3+Y2*Z4-Y4*Z2-Y2*Z3+Y3*Z2
      B(2)=+Y3*Z4-Y4*Z3-Y1*Z4+Y4*Z1+Y1*Z3-Y3*Z1
      B(3)=-Y2*Z4+Y4*Z2+Y1*Z4-Y4*Z1-Y1*Z2+Y2*Z1
      B(4)=+Y2*Z3-Y3*Z2-Y1*Z3+Y3*Z1+Y1*Z2-Y2*Z1
      C(1)=X3*Z4-X4*Z3-X2*Z4+X4*Z2+X2*Z3-X3*Z2
      C(2)=-X3*Z4+X4*Z3+X1*Z4-X4*Z1-X1*Z3+X3*Z1
      C(3)=X2*Z4-X4*Z2-X1*Z4+X4*Z1+X1*Z2-X2*Z1
      C(4)=-X2*Z3+X3*Z2+X1*Z3-X3*Z1-X1*Z2+X2*Z1
      D(1)=-X3*Y4+X4*Y3+X2*Y4-X4*Y2-X2*Y3+X3*Y2
      D(2)=+X3*Y4-X4*Y3-X1*Y4+X4*Y1+X1*Y3-X3*Y1
      D(3)=-X2*Y4+X4*Y2+X1*Y4-X4*Y1-X1*Y2+X2*Y1
      D(4)=+X2*Y3-X3*Y2-X1*Y3+X3*Y1+X1*Y2-X2*Y1
C
C     VOLUME FOR TETRAHEDRON.
      VOL=(A1+B(1)*X1+C(1)*Y1+D(1)*Z1)/6.0

C
C     HYDRAULIC GRADIENTS.
      DHDX=0.0
      DHDY=0.0
      DHDZ=0.0
      DO I=1,4
          DHDX=DHDX+HE(I)*B(I)
          DHDY=DHDY+HE(I)*C(I)
          DHDZ=DHDZ+HE(I)*D(I)
      END DO
      DHDX=DHDX/(6.0*VOL)
      DHDY=DHDY/(6.0*VOL)
      DHDZ=DHDZ/(6.0*VOL)
C
C     VELOCITY COMPONENTS
      KXY=KYX
      KXZ=KZX
      KYZ=KZY
      QX=-KXX*DHDX-KXY*DHDY-KXZ*DHDZ
      QY=-KYX*DHDX-KYY*DHDY-KYZ*DHDZ
      QZ=-KZX*DHDX-KZY*DHDY-KZZ*DHDZ
      RETURN
      END SUBROUTINE
