      SUBROUTINE FORMA(MXX,MCOL,X,Y,KX,NTYPE,NUMNP,NUMEL,NCOL,ETA,XI,W,
     &     CONST,ADOT,FRACT,BDROCK,PSURF,RHOI,FLOWA,SLDGB,
     &     T,KODE,NUMGBC,IBFLUX,BFLUX,
     &        MBAND,LM,AADOT,AFRACT,ABDRCK,AFLOWA,ASLDGB,D,B,A)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION AADOT(MXX),AFRACT(MXX),AFLOWA(MXX)
      DIMENSION ABDRCK(MXX),ASLDGB(MXX),T(MXX),KODE(MXX),IBFLUX(MXX,2)
      DIMENSION BFLUX(MXX)
      DIMENSION ADOT(MXX),FRACT(MXX),PSURF(MXX),FLOWA(MXX)
      DIMENSION BDROCK(MXX),SLDGB(MXX)
      DIMENSION KX(MXX,4),X(MXX),Y(MXX),NTYPE(MXX),D(MXX),B(MXX)
      DIMENSION A(MCOL,MXX),LM(MXX),CONST(MXX)
      DIMENSION P(5),S(5,5),DD(5)
      DIMENSION DPSIX(9),DPSIY(9),DXDS(2,2),DSDX(2,2)
      DIMENSION PSI(4),DPSI(4,2)
      DIMENSION XY(2,4),XI(2,9),ETA(2,9),W(2,9)
C **********************************************************************
C     FORM CONDUCTIVITY MATRIX FOR COMPLETE BODY
C **********************************************************************
      DO 110 I=1,NUMNP
      D(I)=0.0
      B(I)=0.0
      DO 110 J=1,NCOL
  110 A(J,I)=0.0
      MBAND=(NCOL-1)/2
C     WRITE(12,2003)
C
      DO 160 N=1,NUMEL
      IF(NTYPE(N).EQ.1) THEN
        NODEN=4
        NINT=9
      ELSE
        NODEN=3
        NINT=4
      ENDIF
      DO 125 I=1,4
      LM(I)=KX(N,I)
125   CONTINUE
C
C
C     2. FORM ELEMENT CONDUCTIVITY MATRIX
C
      DO 130 I=1,4
      DD(I)=0.0
      P(I)=0.0
      DO 130 J=1,4
  130 S(I,J)=0.0
C
      I=LM(1)
      J=LM(2)
      K=LM(3)
      L=LM(4)
      XY(1,1)=X(I)
      XY(1,2)=X(J)
      XY(1,3)=X(K)
C     IF(NTYPE(N).EQ.1) XY(1,4)=X(L)
      XY(2,1)=Y(I)
      XY(2,2)=Y(J)
      XY(2,3)=Y(K)
C     IF(NTYPE(N).EQ.1) XY(2,4)=Y(L)
      IF(NODEN.EQ.4) THEN
        XY(1,4)=X(L)
        XY(2,4)=Y(L)
      ENDIF
C
C FORM ELEMENT MATRIX AND VECTORS
C
C BEGIN INTEGRATION POINT LOOP
      DO 139 L=1,NINT
      CALL SHAPE(NTYPE(N),XI(NTYPE(N),L),ETA(NTYPE(N),L),PSI,DPSI)
C CALCULATE DXDS...EQUATION (5.3.6)
      DO 133 I=1,2
      DO 133 J=1,2
      DXDS(I,J)=0.0
      DO 133 K=1,4
 133   DXDS(I,J)=DXDS(I,J)+DPSI(K,J)*XY(I,K)
C CALCULATE DSDX...EQUATION (5.2.7)
      DETJ=(DXDS(1,1)*DXDS(2,2)-DXDS(1,2)*DXDS(2,1))
      IF(DETJ.LE.0.0) GOTO 141
      DENOM=1./DETJ
      DSDX(1,1)=DXDS(2,2)*DENOM
      DSDX(2,2)=DXDS(1,1)*DENOM
      DSDX(1,2)=-DXDS(1,2)*DENOM
      DSDX(2,1)=-DXDS(2,1)*DENOM
C CALCULATE D(PSI)/DX...EQUATION (5.3.5)
      DO 136 I=1,4
      DPSIX(I)=DPSI(I,1)*DSDX(1,1)+DPSI(I,2)*DSDX(2,1)
      DPSIY(I)=DPSI(I,1)*DSDX(1,2)+DPSI(I,2)*DSDX(2,2)
 136  CONTINUE
C ACCUMULATE INTEGRATION POINT VALUES OF INTEGRALS
      FAC=DETJ*W(NTYPE(N),L)
      DO 138 I=1,4
      DD(I)=DD(I)+PSI(I)*FAC
C THIS IS LUMPED CAPACITANCE MATRIX
      P(I)=P(I)+AADOT(N)*PSI(I)*FAC
      DO 138 J=1,4
C
      IF(CONST(N).GT.1.E-30) THEN
        TERM1=CONST(N)*(DPSIX(I)*DPSIX(J)+DPSIY(I)*DPSIY(J))
        S(I,J)=S(I,J)+TERM1*FAC
      ENDIF
 138  CONTINUE
 139  CONTINUE
      GOTO 142
141   WRITE(12,1100) DETJ,((XY(MM,NN),NN=1,4),MM=1,2)
1100  FORMAT(' BAD JACOBIAN at 161',G13.6,/,4G13.6,/,4G13.6)
      STOP
142   CONTINUE
C
C     3. ADD ELEMENT CONDUCTIVITY TO COMPLETE CONDUCTIVITY MATRIX
C
      DO 155 L=1,4
      I=LM(L)
      D(I)=D(I)+DD(L)
C THIS (D) IS LUMPED CAPACITANCE MATRIX
      B(I)=B(I)+P(L)
      DO 155 M=1,4
      J=LM(M)-I+1
C OPTIONAL CHECK OF BANDWIDTH, CAN BE REMOVED FOR MORE SPEED
C     IF(J.GT.NCOL) THEN
C     WRITE(12,3000)
C     STOP
C     ENDIF
3000  FORMAT('INCREASE BANDWIDTH AND DIMENSION OF A()')
      A(J+MBAND,I)=A(J+MBAND,I)+S(L,M)
      ZZZ=J+MBAND
  155 CONTINUE
  160 CONTINUE
C
C
C     BOUNDARY CONDITIONS
C
C
      IF (NUMGBC) 2102,2102,1901
1901  CONTINUE
      DO 2101 N=1,NUMGBC
      I=IBFLUX(N,1)
      J=IBFLUX(N,2)
      TEMP=BFLUX(N)
      XL=SQRT((X(J)-X(I))**2+(Y(J)-Y(I))**2)
      TEMP=XL*TEMP*.5
      B(I)=B(I)+TEMP
      B(J)=B(J)+TEMP
2101  CONTINUE
2102  CONTINUE
C REPEAT FROM ADJUST HERE
C BRACNCH 215 BACK HERE, NOW MOVED
C
C     2. TEMPERATURE BOUNDARY CONDITIONS
C
C  *** THIS IS THE LOOP THAT CAUSES PROBLEMS IN VECTOR COMPILE
      DO 300 N=1,NUMNP
      IF(KODE(N)) 225,300,225
225   CONTINUE
      DO 275 M=MAX0(N-MBAND,1),MIN0(N+MBAND,NUMNP)
      K=N+MBAND-M+1
      B(M)=B(M)-A(K,M)*T(N)
275   CONTINUE
  300 CONTINUE
      DO 1300 N=1,NUMNP
      IF(KODE(N)) 1225,1300,1225
1225  CONTINUE
      DO 1250 M=1,MBAND
      A(M,N)=0.0
1250  A(MBAND+1+M,N)=0.0
      DO 1275 M=MAX0(N-MBAND,1),MIN0(N+MBAND,NUMNP)
      K=N+MBAND-M+1
1275  A(K,M)=0.0
1300  CONTINUE
      DO 1400 N=1,NUMNP
      IF(KODE(N)) 1325,1400,1325
1325  A(MBAND+1,N)=1.0
      B(N)=T(N)
1400  CONTINUE
      RETURN
      END
