C     -------------------
      SUBROUTINE WPARAMS1
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 INPUTS THE HYDRAULIC PARAMETERS FOR EACH ELEMENT,
C     CONSTRUCTS THE CONDUCTANCE AND CAPACITANCE MATRICIES, AND
C     EVALUATES THE GROUNDWATER-STORAGE CHANGE WITHIN A TIMESTEP. THE
C     FIRST ENTRY POINT IS CALLED ONCE BY THE MAIN PROGRAM. THE SECOND
C     ENTRY POINT IS CALLED ONCE FOR EACH TIMESTEP ITERATION. THE THIRD
C     ENTRY POINT IS CALL ONCE FOR EACH TIMESTEP.
C
      INCLUDE 'DIMEN.INC'
      INCLUDE 'WUNITS.INC'
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
C
C     INTERNAL SCALARS AND ARRAYS.
C     INFILE      NAME OF INPUT FILE.
C     OUTFILE     OUTPUT FILE NAME.
C     OUTEXT      EXTENSION FOR OUTPUT FILE.
C     XCNAME      COMPARTMENT NAME.
C     XKX()       PRINCIPAL COMPONENT OF HYDRAULIC CONDUCTIVITY TENSOR
C                 FOR A WEDGE ELEMENT FOR A LOCAL COORDINATE SYSTEM
C                 ALIGNED WITH THE TENSOR AXES.
C     XKY()       PRINCIPAL COMPONENT.
C     XKZ()       PRINCIPAL COMPONENT.
C     SS()        SPECIFIC STORAGE FOR A WEDGE ELEMENT.
C     SY()        SPECIFIC YIELD FOR A WEDGE ELEMENT.
C     KXX()       DIAGONAL COMPONENT OF HYDRAULIC CONDUCTIVITY TENSOR
C                 FOR A WEDGE ELEMENT FOR THE GLOBAL MODEL COORDINATE
C                 SYSTEM.
C     KXY         SYMMETRICAL OFF-DIAGONAL TERM.
C     KXZ         SYMMETRICAL OFF-DIAGONAL TERM.
C     KYX()       OFF-DIAGONAL TERM.
C     KYY()       DIAGONAL TERM.
C     KYZ         SYMMETRICAL OFF-DIAGONAL TERM.
C     KZX()       OFF-DIAGONAL TERM.
C     KZY()       OFF-DIAGONAL TERM.
C     KZZ()       DIAGONAL TERM.
C     ALPHA1()    COUNTERCLOCKWISE ROTATION WITHIN AN ELEMENT OF THE
C                 X-AXIS OF HYDRAULIC CONDUCTIVITY TENSOR IN THE X-Y
C                 PLANE OF GLOBAL MODEL COORDINATE SYSTEM.
C     ALPHA2()    COUNTERCLOCKWISE ROTATION WITHIN AN ELEMENT OF THE
C                 X-AXIS OF HYDRAULIC CONDUCTIVITY TENSOR FROM THE X-Y
C                 PLANE OF GLOBAL MODEL COORDINATE SYSTEM.
C     ALPHA3()    COUNTERCLOCKWISE ROTATION WITHIN AN ELEMENT OF THE
C                 Y-AXIS OF HYDRAULIC CONDUCTIVITY TENSOR FROM THE
C                 X-Y PLANE OF GLOBAL MODEL COORDINATE SYSTEM.
C     XITOP()     FLAG INDICATION THE ELEMENT IS IN TOP OF MESH.
C     XE()        COORDINATES OF NODES IN TETRAHEDRON.
C     YE()        COORDINATES OF NODES IN TETRAHEDRON.
C     ZE()        COORDINATES OF NODES IN TETRAHEDRON.
C     HE()        HEADS AT NODES IN TETRAHEDRON.
C     AE()        ELEMENTAL CONDUCTANCE MATRIX.
C     BE()        ELEMENTAL CAPACITANCE MATRIX.
C     H()         HEADS AT NODES FOR CURRENT TIMESTEP.
C     HL()        HEADS AT NODES FOR LAST TIMESTEP.
C     NGROUP      NUMBER OF PARAMETER GROUPS FOR A COMPARTMENT.
C     FACKX       SCALING FACTOR FOR HYDRAULIC CONDUCTIVITY,
C                 IN THE X-DIRECTION.
C     FACKY       SCALING FACTOR FOR HYDRAULIC CONDUCTIVITY,
C                 IN THE Y-DIRECTION.
C     FACKZ       SCALING FACTOR FOR HYDRAULIC CONDUCTIVITY,
C                 IN THE Z-DIRECTION.
C     FACSS       SCALING FACTOR FOR SPECIFIC STORAGE.
C     FACSY       SCALING FACTOR FOR SPECIFIC YIELD.
C     IECHO1      SWITCH FOR ECHO OF HYDRAULIC PARAMETERS.
C     LL1         FIRST ELEMENT NUMBER WITHIN RANGE FOR ELEMENT GROUP.
C     LL2         LAST ELEMENT NUMBE WITHIN RANGE FOR ELEMENT GROUP.
C     NGROUP      NUMBER OF ROTATIONAL GROUPS FOR COMPARTMENT.
C     XALPHA1     COUNTER-CLOCKWISE ANGLE IN X-Y PLANE FROM X-COORDINATE
C                 AXIS TO X-PRINCIPAL AXIS OF HYDRAULIC-CONDUCTIVITY
C                 TENSOR (DEGREES).
C     XALPHA2     COUNTER-CLOCKWISE ANGEL IN VERTICAL PLANE FROM X-Y
C                 PLANE TO X-PRINCIPAL AXIS OF HYDRAULIC-CONDUCTIVITY
C                 TENSOR (DEGREES).
C     XALPHA3     COUNTER-CLOCKWISE ANGEL IN PLANE PERPENDICULAR TO
C                 X-PRINCIPAL AXIS OF HYDRAULIC-CONDUCTIVITY TENSOR FROM
C                 INTERSECTION WITH X-Y PLANE TO Y-PRINCIPAL AXIS OF
C                 HYDRAULIC-CONDUCTIVITY TENSOR (DEGREES).
C
      CHARACTER*80 INFILE,OUTFILE,OUTEXT
      CHARACTER*20 XCNAME
      REAL*8 KX(MAXNE),KY(MAXNE),KZ(MAXNE),SS(MAXNE),SY(MAXNE)
      REAL*8 KXX(MAXNE2),KYX(MAXNE2),KYY(MAXNE2),KZX(MAXNE2),
     +    KZY(MAXNE2),KZZ(MAXNE2),VOL(MAXNE2),AREA(MAXNE2)
      REAL*8 KXY,KXZ,KYZ
      REAL*8 ALPHA1(MAXNE),ALPHA2(MAXNE),ALPHA3(MAXNE)
      REAL*8 TEMP(3,3),XK(3,3),XH(3,3),XHT(3,3),XXK(3,3)
      REAL*8 H(MAXNN),HL(MAXNN)
      REAL*8 XE(4),YE(4),ZE(4),AE(4,4),BE(4,4),HE(4)
      INTEGER*4 ITOP(MAXNE),XITOP
C
C     COMMON BLOCKS.
C     FNAME()     NAMES OF INPUT FILES.
C     RNAME1      SIMULATION NAME FOR GROUNDWATER-FLOW FILES.
C     RNAME2      SIMULATION NAME FOR SOLUTE-TRANSPORT FILES.
C     NCELL       NUMBER OF COMPARTMENTS.
C     CNAME()     NAMES FOR COMPARTMENTS.
C     NNC()       NUMBER OF NODES IN A COMPARTMENT.
C     NEC()       NUMBER OF WEDGE ELEMENTS WITHIN A COMPARTMENT.
C     X()         X-COORDINATE OF A NODE.
C     Y()         Y-COORDINATE OF A NODE.
C     Z()         Z-COORDINATE OF A NODE.
C     IN()        ELEMENT INCIDENCE FOR A WEDGE ELEMENT.
C     IN2()       ELEMENT INCIDENCES FOR A TETRAHEDRAL ELEMENT,
C                     WELL LINK, FAULT ELEMENT, OR FAULT LINK.
C     IEL()       PARENT WEDGE ELEMENT FOR TETRAHEDRAL ELEMENT.
C     NE          GLOBAL NUMBER OF WEDGE ELEMENTS.
C     NE2         GLOBAL NUMBER OF TETRAHEDRAL ELEMENTS.
C     NE3         GLOBAL NUMBER OF TETRAHEDRAL ELEMENS AND WELL LINKS.
C     NE4         GLOBAL NUMBER OF TETHEDRALAL ELEMENTS, WELL LINKS,
C                     AND FAULT ELEMENTS.
C     NE5         GLOBAL NUMBER OF TETHEDRALAL ELEMENTS, WELL LINKS,
C                     FAULT ELEMENTS, AND FAULT LINKS.
C     NN          GLOBAL NUMBER OF MESH NODES.
C     NN3         GLOBAL NUMBER OF MESH NODES AND WELL NODES.
C     NN5         GLOBAL NUMBER OF MESH NODES, WELL NODES, AND
C                     FAULT NODES.
C     ISTEP       TIME STEP NUMBER WITHIN FLOW SIMULATION.
C     ISTEP1      STARTING STEP NUMBER FOR SIMULATION.
C     ISTEP2      ENDING STEP FOR SIMULATION.
C     DELT        LENGTH OF TIME STEP.
C     STIME       INITIAL ELAPSED TIME FOR SIMULATION.
C     IPRN        FLAG TO INDICATED PRINT OUTPUT FOR TIME STEP.
C     IPLT        FLAG TO INDICATE PLOT OUTPUT FOR TIME STEP.
C     NTSTEP      NUMBER OF TRANSPORT STEPS WITHIN FLOW STEP.
C     ITSTEP      TIME STEP NUMBER WITHIN TRANSPORT SIMULATION.
C     A()         VECTOR OF NON-ZERO COEFFICIENTS FOR CONDUCTANCE
C                     MATRIX.
C     B()         VECTOR OF NON-ZERO COEFFICIENTS FOR CAPACITANCE
C                     MATRIX.
C     IA()        VECTOR INDICATING POSITIONS OF ROW ENDS IN
C                     VECTOR A() OR B().
C     JA()        VECTOR INDICATION COLUMN NUMBER ASSOCIATED WITH
C                     ENTRIES IN VECTORS A() OR B().
C     NIA         NUMBER OF ROWS REPRESENTED IN VECTORS A() OR B().
C     NJA         NUMBER OF NON-ZERO COEFFICIENTS IN VECTORS A() OR B().
C     F()         RIGHTHAND VECTOR FOR SYSTEM OF EQUATIONS.
C     QX()        DARCY GROUNDWATER VELOCITY IN THE X-DIRECTION.
C     QY()        DARCY GROUNDWATER VELOCITY IN THE Y-DIRECTION.
C     QZ()        DARCY GROUNDWATER VELOCITY IN THE Z-DIRECTION.
C     ITRAN       SWITCH FOR TRANSPORT SIMULATION.
C     IRPAR       INTERNAL READ FOR PAR UNIT NUMBER.
C     IWPAR       INTERNAL WRITE FOR PAR UNIT NUMBER.
C
      CHARACTER*80 FNAME(30)
      CHARACTER*12 RNAME1,RNAME2
      CHARACTER*20 CNAME(MAXCEL)
      INTEGER*4 NNC(MAXCEL),NEC(MAXCEL)
      INTEGER*4 IN(MAXNE,6),IN2(MAXNE2,4),IEL(MAXNE2)
      REAL*8 X(MAXNN),Y(MAXNN),Z(MAXNN)
      INTEGER*4 IA(MAXNN+1),JA(MAXNON)
      REAL*8 A(MAXNON),B(MAXNON),F(MAXNN)
      REAL*8 QX(MAXNE2),QY(MAXNE2),QZ(MAXNE2)
      COMMON /NAMES/ FNAME,RNAME1,RNAME2
      COMMON /CELLS/ NCELL,CNAME,NNC,NEC
      COMMON /MESH/ X,Y,Z,IN,IN2,IEL,NE,NE2,NE3,NE4,NE5,NN,NN3,NN5
      COMMON /STEPS/ ISTEP,ISTEP1,ISTEP2,DELT,STIME,IPRN,IPLT,NTSTEP,
     +    ITSTEP
      COMMON /MATRIX/ A,B,IA,JA,NIA,NJA,F
      COMMON /SHARE/ QX,QY,QZ,ITRAN
      SAVE
C
C     OPEN INPUT AND OUTPUT FILES.
      WRITE(*,902)
      OUTEXT='.PAR.OUT'
      INFILE=FNAME(4)
      CALL OPENFILES(RNAME1,INFILE,OUTFILE,OUTEXT,IRPAR,IWPAR)
      WRITE(IWPAR,904)
      WRITE(IWPAR,906) INFILE,OUTFILE
C
C     INITIALIZATIONS.
      DO L=1,NE
          KX(L)=0.0
          KY(L)=0.0
          KZ(L)=0.0
          ALPHA1(L)=0.0
          ALPHA2(L)=0.0
          ALPHA3(L)=0.0
      END DO
      PI=3.14159
C
C     LOOP OVER COMPARTMENTS.
      ILAST=0
      DO ICELL=1,NCELL
C
C         INPUT CONDUCTIVITY, SPECIFIC STORAGE, AND SPECIFIC YIELD.
          READ(IRPAR,908) XCNAME
          READ(IRPAR,910) NGROUP,FACKX,FACKY,FACKZ,FACSS,FACSY,
     +        IECHO1
          WRITE(IWPAR,912) NCELL,ICELL,XCNAME,NGROUP,
     +        FACKX,FACKY,FACKZ,FACSS,FACSY,IECHO1
C
C         CHECK COMPARTMENT NAME.
          IF (XCNAME.NE.CNAME(ICELL)) THEN
              WRITE(*,914)
              WRITE(*,916)
              STOP
          ENDIF
C
C         LOOP OVER PARAMETER GROUPS.
          IF (IECHO1.NE.0) WRITE(IWPAR,917)
          DO IGROUP=1,NGROUP
C
C             INPUT PARAMETER VALUES FOR GROUP.
              READ(IRPAR,918) LL1,LL2,XKX,XKY,XKZ,XSS,XSY,XITOP
              L1=LL1+ILAST
              L2=LL2+ILAST
              IF (IECHO1.NE.0) WRITE(IWPAR,920) LL1,LL2,L1,L2,
     +            XKX,XKY,XKZ,XSS,XSY,XITOP
C
C             SCALE PARAMETER VALUES.
              XKX=XKX*FACKX
              XKY=XKY*FACKY
              XKZ=XKZ*FACKZ
              XSS=XSS*FACSS
              XSY=XSY*FACSY
C
C             ASSIGN PARAMETER VALUES TO GLOBAL ARRAY.
              DO L=L1,L2
                  KX(L)=XKX
                  KY(L)=XKY
                  KZ(L)=XKZ
                  SS(L)=XSS
                  SY(L)=XSY
                  ITOP(L)=XITOP
              END DO
C
C             CHECK FOR ZERO HYDRAULIC CONDUCTIVITY.
              TEST=XKX*XKY*XKZ
              IF (TEST.EQ.0.0) THEN
                  WRITE(*,922)
                  WRITE(*,916)
                  STOP
              ENDIF
C
C             CHECK DIMENSIONING.
              IF ((LL1.GT.NEC(ICELL)).OR.(LL2.GT.NEC(ICELL))) THEN
                  WRITE(*,924)
                  WRITE(*,916)
                  STOP
              ENDIF
          END DO
C
C         CHECK THAT ALL ELEMENTS ASSIGNED PARAMETERS
          DO LL=1,NEC(ICELL)
              L=LL+ILAST
              IF (KX(L).LE.0) THEN
                  WRITE(*,928)
                  WRITE(*,916)
                  STOP
              ENDIF
          END DO
C
C         INPUT AXIS ROTATIONS.
          READ(IRPAR,942) NGROUP,FAC1,FAC2,FAC3,IECHO2
          IF (IECHO2.NE.0) WRITE(IWPAR,944) NCELL,ICELL,XCNAME,NGROUP,
     +        FAC1,FAC2,FAC3,IECHO2
          IF (IECHO2.NE.0) WRITE(IWPAR,946)
          DO IGROUP=1,NGROUP
C
C             INPUT PARAMETER VALUES FOR GROUP.
              READ(IRPAR,948) LL1,LL2,XALPHA1,XALPHA2,XALPHA3
              L1=LL1+ILAST
              L2=LL2+ILAST
              IF (IECHO2.NE.0) WRITE(IWPAR,950) LL1,LL2,L1,L2,
     +            XALPHA1,XALPHA2,XALPHA3
C
C             CHECK DIMENSIONING.
              IF ((LL1.GT.NEC(ICELL)).OR.(LL2.GT.NEC(ICELL))) THEN
                  WRITE(*,924)
                  WRITE(*,916)
                  STOP
              ENDIF
C
C             SCALE PARAMETER VALUES AND CONVERT TO RADIANS.
              XALPHA1=XALPHA1*FAC1*PI/180.0
              XALPHA2=XALPHA2*FAC2*PI/180.0
              XALPHA3=XALPHA3*FAC3*PI/180.0
C
C             ASSIGN PARAMETERS TO GLOBAL ARRAY.
              DO L=L1,L2
                  ALPHA1(L)=XALPHA1
                  ALPHA2(L)=XALPHA2
                  ALPHA3(L)=XALPHA3
              END DO
          END DO
          ILAST=ILAST+NEC(ICELL)
      END DO
C
C     CLOSE FILES AND EXIT WPARAMS.
      WRITE(IWPAR,952)
      WRITE(*,954)
      CLOSE(IRPAR)
      CLOSE(IWPAR)

      RETURN
C
C     FORMATS.
902   FORMAT(1X,'*START READING IN WPARAM')
904   FORMAT(1X,'**START READING IN WPARAM**')
906   FORMAT(/1X,'FILE NAMES:'/
     +    1X,'INPUT FILE..............................',A80/
     +    1X,'OUTPUT FILE.............................',A80)
908   FORMAT(A20)
910   FORMAT(I10,5F10.0,I10)
912   FORMAT(/1X,'AQUIFER CHARACTERISTICS:'/
     +    1X,'NUMBER OF COMPARTMENTS..................',I10/
     +    1X,'COMPARTMENT NUMBER......................',I10/
     +    1X,'COMPARTMENT NAME........................',A20/
     +    1X,'NUMBER OF GROUPS........................',I10/
     +    1X,'FACTOR FOR KX...........................',1PE10.2/
     +    1X,'FACTOR FOR KY...........................',1PE10.2/
     +    1X,'FACTOR FOR KZ...........................',1PE10.2/
     +    1X,'FACTOR FOR SS...........................',1PE10.2/
     +    1X,'FACTOR FOR SY...........................',1PE10.2/
     +    1X,'SWITCH FOR ECHO OF PARAMETERS...........',I10)
914   FORMAT(1X,'*ERROR* COMPARTMENT NAME NOT IN LIST')
915   FORMAT(1X,'*ERROR* NUMBER OF TABLES EXCEEDS DIMENSION')
916   FORMAT(1X,'*ERROR* SIMULATION TERMINATED FOR ERROR IN PARAM')
917   FORMAT(/1X,'HYDRAULIC PARAMETERS:'/
     +    1X,5X,'LOCAL',5X,'LOCAL',4X,'GLOBAL',4X,'GLOBAL'/
     +    1X,6X,'ELEM',6X,'ELEM',6X,'ELEM',6X,'ELEM',
     +    8X,'KX',8X,'KY',8X,'KZ',8X,'SS',8X,'SY',7X,'TOP')
918   FORMAT(2I10,5F10.0,I10)
920   FORMAT(1X,4I10,5(1PE10.2),I10)
922   FORMAT(1X,'*ERROR* HYDRAULIC CONDUCTIVITY EQUALS ZERO')
924   FORMAT(1X,'*ERROR* ELEMENT NUMBER EXCEEDS NUMBER OF ELEMENTS')
926   FORMAT(1X,'*ERROR* TABLE NUMBER EXCEEDS NUMBER OF TABLES')
927   FORMAT(1X,'*ERROR* NO TABLE ASSIGNED')
928   FORMAT(1X,'*ERROR* PARAMETERS NOT ASSIGNED TO ALL ELEMENTS')
930   FORMAT(/1X,'UNSATURATED PARAMETER TABLES:'/
     +    1X,'NUMBER OF TABLES........................',I10)
934   FORMAT(I10,F10.0)
936   FORMAT(/1X,'TABLE PARAMETERS:'/
     +    1X,'LOCAL TABLE NUMBER......................',I10/
     +    1X,'GLOBAL TABLE NUMBER.....................',I10/
     +    1X,'NUMBER OF POINTS IN TABLE...............',I10/
     +    1X,'MAXIMUM INTEGRATION SPACING.............',F10.1//
     +    5X,'POINT',3X,'PRESSURE',7X,'KRH',7X,'KRV',9X,'C',8X,'SW')
937   FORMAT(1X,'*ERROR* NUMBER OF POINTS EXCEEDS DIMENSIONS')
938   FORMAT(5F10.0)
940   FORMAT(1X,I10,5(1PE10.2))
942   FORMAT(I10,3F10.0,I10)
944   FORMAT(/1X,'CONDUCTIVITY ORIENTATION:'/
     +    1X,'NUMBER OF COMPARTMENTS..................',I10/
     +    1X,'COMPARTMENT NUMBER......................',I10/
     +    1X,'COMPARTMENT NAME........................',A20/
     +    1X,'NUMBER OF GROUPS........................',I10/
     +    1X,'FACTOR FOR FIRST ANGLE..................',1PE10.2/
     +    1X,'FACTOR FOR SECOND ANGLE.................',1PE10.2/
     +    1X,'FACTOR FOR THIRD........................',1PE10.2/
     +    1X,'SWITCH FOR ECOHO OF PARAMETERS..........',I10)
946   FORMAT(/1X,'DIRECTION PARAMETERS:'/
     +    1X,5X,'LOCAL',5X,'LOCAL',4X,'GLOBAL',4X,'GLOBAL'/
     +    1X,4(6X,'ELEM'),3(5X,'ANGLE'))
948   FORMAT(2I10,3F10.0)
950   FORMAT(1X,4I10,3F10.1)
952   FORMAT(/1X,'**END READING IN PARAM**')
954   FORMAT(1X,'*END   READING IN WPARAMS')
C
C     ---------------------
      ENTRY WPARAMS2(H,ISS)
C     ---------------------
C
C     LOOP OVER ELEMENTS.
      DO L2=1,NE2
          L=IEL(L2)
          XKXX=KX(L)
          XKYY=KY(L)
          XKZZ=KZ(L)
C
C         ROTATE TENSOR INTO GLOBAL COORDINATES
          A1=ALPHA1(L)
          A2=ALPHA2(L)
          A3=ALPHA3(L)
C
C         CONSTRUCT ROTATION MATRIX.
          XH(1,1)=COS(A1)*COS(A2)
          XH(1,2)=-SIN(A1)*COS(A3)-COS(A1)*SIN(A2)*SIN(A3)
          XH(1,3)=SIN(A1)*SIN(A3)-COS(A1)*SIN(A2)*COS(A3)
          XH(2,1)=SIN(A1)*COS(A2)
          XH(2,2)=COS(A1)*COS(A3)-SIN(A1)*SIN(A2)*SIN(A3)
          XH(2,3)=-COS(A1)*SIN(A3)-SIN(A1)*SIN(A2)*COS(A3)
          XH(3,1)=SIN(A2)
          XH(3,2)=COS(A2)*SIN(A3)
          XH(3,3)=COS(A2)*COS(A3)
C
C         CONSTRUCT TRANSPOSE OF RORATION MATRIX.
          XHT(1,1)=XH(1,1)
          XHT(1,2)=XH(2,1)
          XHT(1,3)=XH(3,1)
          XHT(2,1)=XH(1,2)
          XHT(2,2)=XH(2,2)
          XHT(2,3)=XH(3,2)
          XHT(3,1)=XH(1,3)
          XHT(3,2)=XH(2,3)
          XHT(3,3)=XH(3,3)
C
C         CONSTRUCT UNROTATED CONDUCTIVITY MATRIX.
          XK(1,1)=XKXX
          XK(1,2)=0
          XK(1,3)=0
          XK(2,1)=0
          XK(2,2)=XKYY
          XK(2,3)=0
          XK(3,1)=0
          XK(3,2)=0
          XK(3,3)=XKZZ
C
C         MULTIPLY MATRICIES XH*XK*XHT.
          DO II=1,3
              DO JJ=1,3
                  TEMP(II,JJ)=0.0
                  DO J=1,3
                      TEMP(II,JJ)=TEMP(II,JJ)+XH(II,J)*XK(J,JJ)
                  END DO
              END DO
          END DO
          DO II=1,3
              DO JJ=1,II
                  XXK(II,JJ)=0.0
                  DO J=1,3
                      XXK(II,JJ)=XXK(II,JJ)+TEMP(II,J)*XHT(J,JJ)
                  END DO
              END DO
          END DO
C
C         FILL CONDUCTIVITY TENSOR.
          KXX(L2)=XXK(1,1)
          KYX(L2)=XXK(2,1)
          KYY(L2)=XXK(2,2)
          KZX(L2)=XXK(3,1)
          KZY(L2)=XXK(3,2)
          KZZ(L2)=XXK(3,3)
C
C         VOLUME INTEGHRATIONS FOR CONDUCTIVITY AND SPECIFIC STORAGE.
          DO II=1,4
              I=IN2(L2,II)
              XE(II)=X(I)
              YE(II)=Y(I)
              ZE(II)=Z(I)
          END DO
          CALL WSHAPE1(KXX(L2),KYX(L2),KYY(L2),KZX(L2),KZY(L2),KZZ(L2),
     +        SS(L),XE,YE,ZE,AE,BE,VOL(L2),L2)
          VOL(L2)=VOL(L2)*SS(L)
          DO II=1,4
              IROW=IN2(L2,II)
              DO JJ=1,4
                  ICOL=IN2(L2,JJ)
                  CALL PACK2(IROW,ICOL,AE(II,JJ),1)
                  IF (.NOT.(ISS.EQ.1 .AND. ISTEP.EQ.ISTEP1))
     +                 CALL PACK2(IROW,ICOL,BE(II,JJ),2)
              END DO
          END DO
      END DO
C
C     ADD FREE-SURFACE BOUNDARY CONDITION.
      IF (ISS.EQ.1) RETURN
      DO L=1,NE
          AREA(L)=0.0
          IF (ITOP(L).NE.1) CYCLE
          DO II=1,3
              I=IN(L,II)
              XE(II)=X(I)
              YE(II)=Y(I)
          END DO
          CALL WSHAPE2(0.0D+00,SY(L),XE,YE,AE,BE,AREA(L),L)
          AREA(L)=AREA(L)*SY(L)
          DO II=1,3
              IROW=IN(L,II)
              DO JJ=1,3
                  ICOL=IN(L,JJ)
                  CALL PACK2(IROW,ICOL,BE(II,JJ),2)
              END DO
          END DO
      END DO
      RETURN
C
C     --------------------------
      ENTRY WPARAMS3(H,HL,SUMPS)
C     --------------------------
C
C     COMPUTE STORAGE CHANGE WITHIN MESH DOMAIN.
      SUMPS=0.0
      IF (ISS.EQ.1 .AND. ISTEP.EQ.ISTEP1) RETURN
C
C     STORAGE CHANGE DUE TO SPECIFIC STORAGE.
      DO L2=1,NE2
          DELH=0.0
          DO J=1,4
              II=IN2(L2,J)
              DELH=DELH+(H(II)-HL(II))/4.0
          END DO
          SUMPS=SUMPS+DELH*VOL(L2)/DELT
      END DO
C
C     STORAGE CHANGE DUE TO SPECIFIC YIELD.
      DO L=1,NE
          IF (ITOP(L).NE.1) CYCLE
          DELH=0.0
          DO J=1,3
              II=IN(L,J)
              DELH=DELH+(H(II)-HL(II))/3.0
          END DO
          SUMPS=SUMPS+DELH*AREA(L)/DELT
      END DO
      RETURN
C
C     -----------------
      ENTRY WPARAMS4(H)
C     -----------------
C
C     COMPUTE GROUNDWATER DARCY VELOCITIES.
      DO L2=1,NE2
          DO I=1,4
              II=IN2(L2,I)
              XE(I)=X(II)
              YE(I)=Y(II)
              ZE(I)=Z(II)
              HE(I)=H(II)
          END DO
          CALL WSHAPE4(XE,YE,ZE,HE,KXX(L2),KYX(L2),KYY(L2),KZX(L2),
     +        KZY(L2),KZZ(L2),QX(L2),QY(L2),QZ(L2))
      END DO
      RETURN
      END SUBROUTINE
