subroutine cascade(nlortz,wlortz,efc,rmass,tref,rh
     &   ,wav0,wdel,xu,wu,nu,z,tmp,prs,clm,den,ep1,ep2,ntau
     &   ,wav,str,eng,wid,nl,nlines,ht,th,ha)

        implicit double precision (a-h,o-z)
        include "parameters.dat"

c       ----------------------------------------------------
c                         Input Arrays
c       ----------------------------------------------------
        dimension z(1),tmp(1),prs(1),clm(1),den(1),ep1(1),ep2(1)
     &  ,xu(1),wu(1),nl(1),nlines(1),wav(1),str(1),eng(1),wid(1)

c       ----------------------------------------------------
c                         Output Arrays
c       ----------------------------------------------------
        dimension ht(1), th(1), ha(1)

c       ----------------------------------------------------
c                       Internal Arrays
c       ----------------------------------------------------
        dimension h(nd1),temp(nd1),pres(nd1),col(nd1),a(nd1,nd6)
     &  , x1(nd1), x2(nd1)

        data tmp0 / 296.d0 /
        data tmax / 7.d0 /

        do n = 1, ntau
           ha(n) = zero
           ht(n) = zero
           th(n) = zero
        end do

        do n = 1, ntau
           x2(n) = ep2(n)/(one+ep2(n))
           x1(n) = ep1(n)/(one+ep1(n))
        end do

        do n = 1, ntau-1
           col(n) = dsqrt(clm(n)*clm(n+1))
           temp(n) = half*(tmp(n)+tmp(n+1))
           pres(n) = dsqrt(prs(n)*prs(n+1))
        end do
        col(ntau) = col(ntau-1)**2/col(ntau-2)
        temp(ntau) = temp(ntau-1)
        pres(ntau) = pres(ntau-1)**2/pres(ntau-2)

c       -------------------------------------------------------
c                      Loop through near IR bands
c       -------------------------------------------------------
        do jb = 1, 3
        wend = wav(nl(jb)+nlines(jb)-1)
        wlast = wav(nl(jb))

c        do n = nl(jb), nl(jb)+nlines(jb)-1
c           write(unit=*,fmt='(i4,1p2e12.4)')
c     &      n,wav(n),str(n)
c        end do 
c        read (unit=*,fmt='(a)') iret

        do while (wlast .lt. wend)

           wnext = dmin1(wlast+wdel,wend)
           wmid = half*(wlast+wnext)
           delw = wnext - wlast
           solflux = solirflx(wmid)/rh**2

           lw1 = min(locate(wav(nl(jb)),nlines(jb),wlast)+1
     &      ,nlines(jb))+nl(jb)-1
           lw2 = max(locate(wav(nl(jb)),nlines(jb),wnext),1)+nl(jb)-1
           nlw = lw2 - lw1 + 1

           do i = 1, nu
              call random(xu(i),pres,temp,col,ntau,str(lw1)
     &         ,wid(lw1),eng(lw1),nlw,efc,rmass,wmid,delw,a(1,i))
           end do
           
           do n = 1, ntau
              sum = zero
              do i = 1, nu
                if(n .gt. 1) then
                   sum=sum+xu(i)*wu(i)*(a(n,i)-a(n-1,i))
     &              /(col(n)-col(n-1))
                else
                   sum=sum+xu(i)*wu(i)*a(1,i)/col(1)
                endif
              end do
              h(n) = half*solflux*delw*sum
           end do

c           write(unit=*,fmt=*) ' IR heating: ',wmid,a(1,1)

           do n = 1, ntau
              ha(n) = ha(n) + h(n)*den(n)
           end do

c          -----------------------------------------
c              Cascade heating rate (ergs/cm3/s)
c          -----------------------------------------
           fb = float(jb)
           dwr = (wmid-(fb+one)*wav0)/wmid
           rw = wav0/wmid
           fbr = fb*rw
           do n = 1, ntau
              ht(n) = ht(n)+x2(n)*(dwr+fbr*x1(n))*h(n)*den(n)
           end do

c          -----------------------------------------
c              nu-4 source function 
c          -----------------------------------------
           do n = 1, ntau
              th(n) = th(n) + rw*x2(n)*h(n)
           end do

           wlast = wnext

	end do

        end do


        return
        end


        function solirflx(ww)
        implicit double precision (a-h,o-z)
        parameter (one = 1.d0, c1 = 8.37d-10, c2 = 2.4945d-4)
        solirflx = c1*(ww**3)/(dexp(c2*ww) - one)
        return
        end