subroutine drybicki_lte(tmin,dtau,bp,dbp,ntau,xu,wu,nu,opp,dopp
     & ,wpp,npp,bb,h,f)

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

c     ---------------------------------------------------------------
c                   Input arrays and variables
c     ---------------------------------------------------------------
      dimension dtau(nd1),bp(nd1),dbp(nd1),xu(nd6),wu(nd6)
     & , opp(nd1,nd2), dopp(nd1,nd2), wpp(nd2)

c     ---------------------------------------------------------------
c                         Output arrays
c     ---------------------------------------------------------------
      dimension h(nd1), f(nd1,nd1)

c     ---------------------------------------------------------------
c                   Internal arrays and variables
c     ---------------------------------------------------------------
      dimension gf(nd1), dgf(nd1), dt(nd1), w(nd1), ta(nd1)
     & , tb(nd1), tc(nd1), fj(nd1), xj(nd1)
     & , va(nd1), vb(nd1), vc(nd1), v(nd1,nd1)

c     ---------------------------------------------------------------
c                        Initialize arrays
c     ---------------------------------------------------------------
      do j1 = 1, ntau
         do j2 = 1, ntau
            f(j2,j1) = zero
         enddo
         h(j1) =zero 
      end do
         
c     ---------------------------------------------------------------
c                          Loop through frequencies 
c     ---------------------------------------------------------------
      do jf = 1, npp  

         wgf = wpp(jf)

         do jt = 1, ntau 
            gf(jt) = opp(jt,jf)
            dgf(jt) = dopp(jt,jf)
         end do

c        ----------------------------------------------
c            Calculate monochromatic optical depth
c        ----------------------------------------------
         do jt = 1, ntau - 1
            dt(jt) = half*(gf(jt)+gf(jt+1))*dtau(jt)
c            dt(jt) = dmax1(half*(gf(jt)+gf(jt+1))*dtau(jt),tmin)
         end do

c        ----------------------------------------------
c                Find optically thin level
c        ----------------------------------------------
         tmax = zero 
         do jt = 1, ntau -1
            tmax = tmax + dt(jt)
         end do

         k = 1
         tot = dt(k) 
         do while((tot .lt. tmin).and.(k.lt.ntau))
            k = k + 1
            tot = tot + dt(k)
         end do
         ntmn = max(k,1)

c         if (iverbose .ge. 2) write (*,*) "drybicki_lte: ntmn = ",ntmn

         if ((tot .lt. tmin).or.(ntmn .eq. ntau)) then

            do jt = 1, ntau
               h(jt)=h(jt)+wgf*gf(jt)*(half*bb-bp(jt))
            end do

            do jt = 1, ntau
               f(jt,jt)=f(jt,jt)+wgf*dgf(jt)*(half*bb-bp(jt))
     &          -wgf*gf(jt)*dbp(jt)
 	    end do

         else

c         ------------------------------------------------
c                    Loop through zenith angles
c         ------------------------------------------------
          do ju = 1, nu

c            -------------------------------------
c                Calculate W matrix (Diagonal)
c            -------------------------------------
             wght=wu(ju)*wgf
             do jt = 1, ntau
               w(jt)=wght*gf(jt)
             end do

c            -------------------------------------
c               Calculate T matrix (Tridiagonal)
c            -------------------------------------
             do jt = 1, ntmn-1
                ta(jt) = zero
c                tb(jt) = one
c                tc(jt) = -one
                tb(jt) = one/tmin/tmin
                tc(jt) = -one/tmin/tmin
             end do
             call tmatrix(xu(ju), dt(ntmn), ntau-ntmn+1
     &        , ta(ntmn), tb(ntmn), tc(ntmn))

c           -------------------------------------
c              Calculate LU decomposition of T
c           -------------------------------------
c            call ludcmpt(ta,tb,tc,ntau)

C           -------------------------------------
C           Calculate Feautrier Intensities, fj
C           -------------------------------------
            do jt = 1, ntmn - 1
               fj(jt) = zero
            end do
            do jt = ntmn, ntau
               fj(jt) = bp(jt)
            end do
            fj(ntau) = fj(ntau) + (two*xu(ju)/dt(ntau-1))*bb

c           call lubksbt(ta,tb,tc,fj,ntau)
            call tridag(ta,tb,tc,fj,ntau,ierr)
               if(ierr .eq. 1) then
                   write(unit=*,fmt=*) ' Error in drybicki_lte'
                   do k = 1, ntau
                      write(unit=*,fmt=*) ta(k),tb(k),tc(k)
                   end do
                   stop
               end if

c           -------------------------------------
c                 Add into heating rate
c           -------------------------------------
            do jt = 1, ntau
               h(jt) = h(jt)+w(jt)*(fj(jt)-bp(jt))
            end do

c           -------------------------------------
c             Add derivative of W terms into F
c           -------------------------------------
            do jt = 1, ntau
               f(jt,jt)=f(jt,jt)+wght*dgf(jt)*(fj(jt)-bp(jt))
     &          -w(jt)*dbp(jt)
 	    end do

c           -------------------------------------
c             Calculate V matrix (Tridiagonal)
c           -------------------------------------
            fup = bb
            fdn = zero
            call vmatrix(xu(ju),dtau(ntmn),gf(ntmn),dgf(ntmn)
     &       ,fj(ntmn),ntau-ntmn+1,fup,fdn,va(ntmn),vb(ntmn)
     &       ,vc(ntmn))

            do j1 = 1, ntau
               do j2 = 1, ntau
                  v(j2,j1) = zero
               end do
            end do

            v(ntmn,ntmn) = dbp(ntmn)-vb(ntmn)
            v(ntmn,ntmn+1) = -vc(ntmn)
            do j1 = ntmn+1, ntau-1
               v(j1,j1-1) = -va(j1)
               v(j1,j1) = dbp(j1)-vb(j1)
               v(j1,j1+1) = -vc(j1)
            end do 
            v(ntau,ntau-1) = -va(ntau)
            v(ntau,ntau) = dbp(ntau)-vb(ntau) 

c           -------------------------------------
c               Calculate T(-1)*V store in V
c           -------------------------------------
c            do jt = 1, ntau
c               call lubksbt(ta,tb,tc,v(1,jt),ntau)
c            end do
            do jt = 1, ntau
               call tridag(ta,tb,tc,v(1,jt),ntau,ierr)
               if(ierr .eq. 1) then
                   write(unit=*,fmt=*) ' Error in drybicki_lte'
                   do k = 1, ntau
                      write(unit=*,fmt=*) ta(k),tb(k),tc(k)
                   end do
                   stop
               end if
            end do

c           -------------------------------------
c              Calculate W*T(-1)*V add into F
c           -------------------------------------
  	    do j1 = 1, ntau
  	       do j2 = 1, ntau
  	          f(j2,j1)=f(j2,j1)+w(j2)*v(j2,j1)
  	       end do
            end do

         enddo
      end if
      enddo

      return
      end