c    @(#)set3d.f	1.1
c
c
c
c
      subroutine set3d(kz,ifdiff)
c
c --- subroutine set3d sets the stratification and the
c --- associated matrices and constants
c
      include 'qg_par.inc'
      common/threed/hz(k0),scsig,sigz(k01),eigval(k0),ainv(k0,k0),
     +      amat(k0,k0)
      common/bsolv/ab(mx0),bb(mx0),cb(mx0),workb(mnwork)
      common/unit/iin,iout,idff,idif,iorr,iowr,iplt,ifndd
      dimension hz1(k0),strat(k0,k0)
      dimension e1(k0),e2(k0,k0),iwork(2*k0)
c
c --- ifdiff=1/0 - finite diff/collocation
c
      if(ifdiff.eq.0)go to 50
c
c --- calculate the related depths hz1
c
      ierr=0
      do 10 k=2,kz
        hz1(k)=0.5*(hz(k-1)+hz(k))
   10 continue
c
c --- step 1 - set the stratification matrix strat
c
      do 15 i=1,kz
      do 15 j=1,kz
        strat(i,j)=0.
   15 continue
      kzm1=kz-1
      do 20 i=1,kzm1
        a1=scsig*sigz(i+1)/(hz(i)*hz1(i+1))
        strat(i,i)=-a1
        strat(i,i+1)=a1
   20 continue
      do 25 i=2,kz
        a1=scsig*sigz(i)/(hz(i)*hz1(i))
        strat(i,i-1)=a1
        strat(i,i)=strat(i,i)-a1
   25 continue
c
c --- step 2 - get eigenvalues (eigval) and eigenvectors (amat)
c
      call eigrg1(k0,kz,strat,kz,eigval,e1,amat,e2,workb,ierr)
c
c --- check the results
c
      if(ierr.ne.0) go to 99
      do 30 i=1,kz
        if(e1(i).ne.0) go to 99
   30 continue
      xkz=float(kz)
      xkz=sqrt(xkz)
c
c --- step 3 - normalize the eigenvectors
c
      do 45 i=1,kz
        a1=0.
        do 35 i1=1,kz
          a1=a1+amat(i1,i)**2
   35   continue
        do 40 j=1,kz
          amat(j,i)=xkz*amat(j,i)/sqrt(a1)
   40   continue
   45 continue
c
c --- come directly here if using the collocation method
c
   50 continue
c
      call sort(eigval,amat,kz)
c
c --- step 4 - obtain inverse of amat(ainv).  since amat will be
c --- used later, a dummy array will be used in the inversion
c --- to prevent amat from being destroyed
c
      do 55 i=1,kz
      do 55 j=1,kz
        e2(i,j)=amat(i,j)
   55 continue
c
      call invmtx(e2,k0,ainv,k0,kz,det,iwork,ierr)
c
      if(ierr.ne.0) go to 990
      if(ifdiff.ne.0) then 
           write(iout,"(' Finite Difference method   Eigenvalues:')")
           write(iout,60)(eigval(i),i=1,kz)
      endif
      if(ifdiff.eq.0) then 
           write(iout,"('Collocation in depth.  Eigenvalues:')")
           write(iout,60)(eigval(i),i=1,kz)
      endif
    
   60 format(5(1x,1pe11.4,/))


      do 70 i=1,kz
        write(iout,75) i
        write(iout,60)(amat(j,i),j=1,kz)
   70 continue
   75 format(/'Eigenvector(',i1,') :')
      return
c
c --- come here if there is a problem getting the eigenvalues
c --- or eigenvectors
c
   99 continue
      write(iout,100)
  100 format(/'Something wrong in EIGRG1--program halted')
      write(iout,101)ierr
  101 format(' Ierr = ',i4)
      do 105 iii=1,kz
        write(iout,110)(e2(iii,jjj),jjj=1,kz)
  105 continue
  110 format(4(2x,d13.6))
      write(iout,115)
  115 format(///' E1')
      write(iout,110)(e1(iii),iii=1,kz)
      stop
c
c --- come here if there is a problem in the matrix inverter
c
  990 continue
      write(iout,991)
  991 format(/'Something wrong in INVMTX--program halted')
      stop
      end
