subroutine solve (a,b,x,n,iortho)
      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), b(maxcr6), x(maxcr6), sum,
     1 whammy, dabs
      nplus1=n+1
      nminus=n-1
c  find updated b vector
      do 20 j=1,nminus
      jmax=a(nplus1,j)
      whammy=b(jmax)
      b(jmax)=b(j)
      b(j)=whammy
      jplus1=j+1
      do 10 i=jplus1,n
   10 b(i)=b(i)+a(i,j)*b(j)
   20 continue
c  backsolve to find the solution
      i=n
      if (dabs(a(n,n)).gt.1.0d-7) x(n)=b(n)/a(n,n)
   30 do 80 j=2,n
      if (j.eq.2 .and. dabs(a(n,n)).le.1.0d-7) go to 50
      i=n-j+1
c  do 12 i=n-1,1
      sum=0.0
      iplus1=i+1
      do 40 l=iplus1,n
      if (dabs(a(i,l)).le.1.0d-8) go to 40
      sum=sum+a(i,l)*x(l)
   40 continue
      if (dabs(a(i,i)).gt.1.0d-7) go to 70
c  if multiple solutions exist, set x(i)=1 and continue
   50 x(i)=1.0
      if (i.eq.n) go to 30
      if (iortho.eq.0) go to 80
      ip1=i+1
      do 60 k=ip1,n
      if (dabs(a(i,k)).gt.1.0d-7) x(k)=0.0d0
   60 continue
      go to 80
   70 x(i)=(b(i)-sum)/a(i,i)
   80 continue
      return
      end