SUBROUTINE MATINV (cxxx,Mx,JL)
C
C     MATRIX INVERTER FOR MATRIX B(Mx,Mx)
C
c     work arrays of the size m are needed below, set ms= largest
c     value of m in any call to matinv
c       
c     Method: Gaussian Elimination
c         declare variables double precision for greater
c         accuracy
c
      parameter msx=94
	REAL*8 cxxx(94,94)
      real*8 Bx(msx,msx),Cx(msx,msx)
      real*8 B1x(msx),B2x(msx),C2x(msx)
	do 11 i=1,mx
	do 11 j=1,mx
	bx(i,j)=cxxx(i,j)
11	continue
       do  i=1,mx
	b1x(i)=0.0
	b2x(i)=0.0
	c2x(i)=0.0
       enddo
      M1=Mx-1
      DO 20 J=1,Mx
      DO 10 I=1,Mx
   10 Cx(I,J)=0.0
   20 Cx(J,J)=1.0D0
      DO 100 I=1,M1
c       if(b(i,i) .eq. 0.0)b(i,i)=1.e-30
      B1x(1)=1.0D0/Bx(I,I)
C	IF(I.GE.80 .AND. JL.EQ.37)THEN
c	type *,'b1x bx i',b1x(1),bx(i,i),i
C	ENDIF
      Bx(I,I)=1.0D0
c       type *,'after do 100',b1x(1),bx(i,i),i,mx
      DO 30 J=1,I
      Cx(I,J)=Cx(I,J)*B1x(1)
C	IF(J.LE.2)THEN
C	IF(I.GE.80 .AND. JL.EQ.37)THEN
c	type *,'cx b1x I J',cx(i,j),b1x(1),i,j
C	ENDIF
C	ENDIF
   30 CONTINUE
      IP1=I+1
      DO 40 J=IP1,Mx
      Bx(I,J)=Bx(I,J)*B1x(1)
   40 CONTINUE
      DO 50 I1=IP1,Mx
      B1x(I1)=Bx(I1,I)
c	type *,'in do 50 loop',b1x(i1)
   50 CONTINUE
c        type *,'after 50 continue',i,mx
      DO 60 I1=IP1,Mx
      Bx(I1,I)=0.0
c	type *,' do 60 bx(i1,i)',bx(i1,i)
   60 CONTINUE
c        type *,'after 60 continue',i,Mx
      DO 70 J=1,Mx
      B2x(J)=Bx(I,J)
      C2x(J)=Cx(I,J)
C	IF(J.LE.2)THEN
C	IF(I.GE.80 .AND. JL.EQ.37)THEN
c	type *,' do 70 b2x(j) c2x(j) J I',b2x(j),c2x(j),J,I
C	ENDIF
C	ENDIF
   70 CONTINUE
c        type *,'after 70 continue',i,Mx
      DO 80 I1=IP1,Mx
      DO 80 J=1,I
      Cx(I1,J)=Cx(I1,J)-B1x(I1)*C2x(J)
C	IF(I.GE.80 .AND. JL.EQ.37)THEN
C	type *,' do 80 cx(i1,j) I1 J',cx(i1,j),I1,J
C	ENDIF
   80 CONTINUE
c        type *,'after 80 continue  i ip1',i,ip1,Mx 
      DO 90 I1=IP1,Mx
      DO 90 J=IP1,Mx
      Bx(I1,J)=Bx(I1,J)-B1x(I1)*B2x(J)
C	if(i.GT.79)then
C	IF(JL.EQ.37)THEN
C	if(i1.eq.92)then
C	if(j.gt.90 .and. j.lE.92)then
c	type *,' do 90 bx(i1,j)  I i1 j MX',bx(i1,j),I,i1,j,Mx
C	endif
C	endif
C	ENDIF
C	endif
   90 CONTINUE
c       type *,'after 90 continue',i
  100 CONTINUE
c       if(b(m,m) .eq. 0.0)b(m,m)=1.e-30
      B1x(1)=1.0/Bx(Mx,Mx)
      Bx(Mx,Mx)=1.0
C	IF(JL.EQ.37)THEN
c	type *,'after 100 cont b1x bx',b1x(1),bx(mx,mx)
C	ENDIF
      DO 110 J=1,Mx
      Cx(Mx,J)=Cx(Mx,J)*B1x(1)
c	type *,'do 110 cx',cx(mx,j)
  110 CONTINUE
c        type *,'after 110 continue' 
      DO 150 I=2,Mx
      I1=Mx-I+2
      DO 120 I2=1,I1
      B1x(I2)=Bx(I2,I1)
c	type *,'do 120 b1x',b1x(i2)
  120 CONTINUE
c        type *,'after 120 continue'
      DO 130 J=1,Mx
      C2x(J)=Cx(I1,J)
c	type *,'do 130 c2x',c2x(j)
  130 CONTINUE
c        type *,'after 130 continue'
      IM1=I1-1
      DO 140 I2=1,IM1
      DO 140 J=1,Mx
      Cx(I2,J)=Cx(I2,J)-B1x(I2)*C2x(J)
c	type *,'do 140 cx',cx(i2,j)
  140 CONTINUE
c        type *,'after 140 continue' 
 150  CONTINUE
      DO 160 I=1,Mx
      DO 160 J=1,Mx
      Bx(I,J)=Cx(I,J)
c	type *,'do 160 bx',bx(i,j)
  160 CONTINUE
	do 111 i=1,mx
	do 111 j=1,mx
	cxxx(i,j)=bx(i,j)
111	continue
      RETURN
      END