subroutine gauss (a,n)
      integer maxcor, maxcr2, maxcr4, maxcr5, maxcr6, maxcr7, maxcr8
c     parameter (maxcor=2, maxcr2=maxcor*maxcor, maxcr4=maxcr2*maxcr2,
c    1 maxcr5=maxcr2*(maxcr2+1)/2, maxcr6=maxcor*2, maxcr7=maxcr6+1,
c    2 maxcr8=maxcr6*maxcr7+1)
      parameter (maxcor=2, maxcr2=4, maxcr4=16,
     1 maxcr5=10, maxcr6=4, maxcr7=5,
     2 maxcr8=21)
      double precision a(maxcr7,maxcr6), max, dummy(maxcr6), factor,
     1 dabs
      nminus=n-1
      do 80 ij=1,nminus
      ijpl1=ij+1
c  find maximum non-zero element in ij column from ij to n rows
      max=a(ij,ij)
      jmax=ij
      do 20 i=ijpl1,n
      if (dabs(max).le.1.0d-9) go to 10
      if (dabs(a(i,ij)).le.1.0d-9) go to 20
      if (dabs(a(i,ij)).le.dabs(max)) go to 20
   10 max=a(i,ij)
      jmax=i
   20 continue
      if (max.eq.0.0) go to 60
      a(n+1,ij)=dble(jmax)
c swap places with ij row & jmax row
      if (ij.eq.i) go to 40
      do 30 i=ij,n
      dummy(i)=a(ij,i)
      a(ij,i)=a(jmax,i)
   30 a(jmax,i)=dummy(i)
c  multiply all rows after the ij row by multiplier
   40 if (a(ij,ij).eq.0.0) go to 60
      do 70 i=ijpl1,n
      factor=-(a(i,ij)/a(ij,ij))
      do 50 j=ij,n
   50 a(i,j)=a(i,j)+factor*a(ij,j)
      go to 70
c if the column is all zeros, set a(i,ij)=factor=0
   60 factor=0.0
      a(n+1,ij)=dble(ij)
   70 a(i,ij)=factor
   80 continue
      a(n+1,n)=dble(n)
      return
      end