program decomp

c **********************************************************************
c *                                                                    *
c *                                                                    *
c *      Program decomp calcualtes the eigenvector expansion of an     *
c *      arbitrary temperature profile.                                *  
c *                                                                    *
c *                                                                    *
c *                                                    Roger Yelle     *
c *                                                    27 JUNE 1991    *
c *                                                                    *
c **********************************************************************

        implicit double precision (o-z, a-h)
        include "parameters.dat"
        character fname*64, header*64

c       ---------------------------------------------------------
c            Arrays : Model atmosphere
c       ---------------------------------------------------------
        dimension tmp1(nd1),tmp2(nd1),dtmp(nd1),coef(nd1)
     &  , evalue(nd1), evector(nd1,nd1), a(nd1,nd1), indx(nd1)

c       ----------------------------------------------------------
c                        Read reference atmosphere
c       ----------------------------------------------------------
        ntau = 51
        open (unit=30, file='heat_a.out') 
           do n = 1, ntau
              read(unit=30, fmt=906) zn,tmp1(n)
           end do
        close(unit=30)
  906   format(1x,0pf7.1,1x,0pf7.2)
c       ----------------------------------------------------------
c                        Read perturbed atmosphere
c       ----------------------------------------------------------
        open (unit=30, file='heat_hot_b.out') 
           do n = 1, ntau
              read(unit=30, fmt=906) zn,tmp2(n)
           end do
        close(unit=30)
c       ----------------------------------------------------------
c                    Read eigenvectors and eigenvalues
c       ----------------------------------------------------------
        open (unit=30, file='eigenvectors.dat') 
           do n = 1, ntau
              read(unit=30,fmt='(a)') header 
              read(unit=30,fmt=907) zn,(evector(n,j),j=1,ntau)
           end do
        close(unit=30)
  907   format(1x,1p10e11.3)

        open (unit=30, file='eigenvalues.dat') 
           do n = 1, ntau
              read(unit=30,fmt='(1x,1pe15.7)') evalue(n)
           end do
        close(unit=30)

c       ----------------------------------------------------------
c              Calculate eigenvector decomposition
c       ----------------------------------------------------------
        do n = 1, ntau
           dtmp(n) = tmp2(n) - tmp1(n)
           coef(n) = dtmp(n)
           do j = 1, ntau
              a(j,n) = evector(j,n)
           end do
        end do

        call ludcmp(a, ntau, nd1, indx, dinx, ierr)
        call lubksb(a, ntau, nd1, indx, coef)

c         check
        do n = 1, ntau
           sum = zero
           do j = 1, ntau
              sum = sum + coef(j)*evector(n,j)
           end do
           write(unit=*,fmt=*) n,sum,dtmp(n)
        end do 

c       -----------------------------------------------------------
c                         write results to disc
c       -----------------------------------------------------------
        write(unit=*,fmt=920)
        read (unit=*,fmt='(a)') fname
        open(unit=60,file=fname)
           do j = 1, ntau
              write(unit=60,fmt='(1x,i3,1p3e13.5)') 
     &        j,coef(j),evalue(j),one/evalue(j)
           end do
        close(unit=60)
 920    format(' Enter output filename > ',$)

        stop
        end

      subroutine ludcmp(a, n, np, indx, d, ierr)
      parameter (nmax = 1106, tiny = 1.0d-30, zero = 0.d0)
      implicit double precision (o-z, a-h)
      dimension a(np, np), indx(n), vv(nmax)

      ierr = 0
      d = 1.d0
      do 12 i = 1, n
         aamax = 0.d0
         do 11 j = 1, n
            if (dabs(a(i,j)) .gt. aamax) aamax = dabs(a(i,j))
   11    continue
         if (aamax .eq. 0.d0) then
            write(unit=*,fmt=*) 'Singular matrix.'
            ierr = 1
            return
         end if
         vv(i) = 1.d0 / aamax
   12 continue

      do 19 j = 1, n
         if (j .gt. 1) then
            do 14 i = 1, j - 1
               sum = a(i,j)
               if (i .gt. 1) then
                  do 13 k = 1, i - 1
                     sum = sum - (a(i,k) * a(k,j))
   13             continue
                  a(i,j) = sum
               end if
   14       continue
         end if
         aamax = 0.d0
         do 16 i = j, n
            sum = a(i,j)
            if (j .gt. 1) then
               do 15 k = 1, j - 1
                  sum = sum - (a(i,k) * a(k,j))
   15          continue
               a(i,j) = sum
            end if
            dum = vv(i) * dabs(sum)
            if (dum .ge. aamax) then
               imax = i
               aamax = dum
            end if
   16    continue
         if (j .ne. imax) then
            do 17 k = 1, n
               dum = a(imax,k)
               a(imax,k) = a(j,k)
               a(j,k) = dum
   17       continue
            d = - d
            vv(imax) = vv(j)
         end if
         indx(j) = imax
         if (j .ne. n) then
            if (a(j,j) .eq. 0.d0) a(j,j) = tiny
               dum = 1.d0 / a(j,j)
               do 18 i = j + 1, n
                  a(i,j) = a(i,j) * dum
   18          continue
         end if
   19 continue
      if (a(n,n) .eq. 0.d0) a(n,n) = tiny
      return 
      end


      subroutine lubksb(a, n, np, indx, b)
      implicit double precision (o-z, a-h)
      dimension a(np, np), indx(n), b(n)
      ii = 0
      do 12 i = 1, n
         ll = indx(i)
         sum = b(ll)
         b(ll) = b(i)
         if (ii .ne. 0) then
            do 11 j = ii, i - 1
               sum = sum - (a(i,j) * b(j))
   11       continue
         else if (sum .ne. 0.d0) then
            ii = i
         end if
         b(i) = sum
   12 continue

      do 14 i = n, 1, -1
         sum = b(i)
         if (i .lt. n) then
            do 13 j = i + 1, n
               sum = sum - (a(i,j) * b(j))
   13       continue
         end if
         b(i) = sum / a(i,i)
   14 continue

      return 
      end