subroutine rybicki_lte(dtau, bp, ntau, xu, wu, nu
     & , opp, wpp, npp, b1, q)

      implicit double precision (o-z, a-h)
      parameter (zero = 0.d0, half = 0.5d0, one = 1.d0, two = 2.d0
     & , nd1 = 101)

c     --------------------------------------------------------
c                        Input Arrays
c     --------------------------------------------------------
      dimension dtau(1), bp(nd1,1), xu(1), wu(1)
     & , opp(nd1, 1), wpp(1),  b1(1)

c     --------------------------------------------------------
c                        Output Arrays
c     --------------------------------------------------------
      dimension q(1)

c     --------------------------------------------------------
c                        Internal Arrays
c     --------------------------------------------------------
      dimension s(nd1), e(nd1,nd1), ta(nd1), tb(nd1), tc(nd1)
     & , w(nd1), tu(nd1,nd1), gf(nd1), dt(nd1)
      integer indx(nd1)

      data tmin / 1.d-10 /

c     --------------------------------------------------------
c                      Initialize Arrays
c     --------------------------------------------------------
      do k = 1, ntau
         q(k) = zero
         do m = 1, ntau
            e(m,k) = zero
         end do
         e(k,k) = one
      end do

c     ---------------------------------------------------------
c                  Loop Through Frequencies 
c     ---------------------------------------------------------
      do j = 1, npp

            bf1 = b1(j)
            wgf = wpp(j)
            do k = 1, ntau
               gf(k) = opp(k,j)
            end do

c           -----------------------------------------------
c                  Monochromatic optical depth
c           -----------------------------------------------
            do k = 1, ntau - 1
               dt(k) = half*(gf(k)+gf(k+1))*dtau(k)
            end do

c           -----------------------------------------------
c                 find optically thin level
c           -----------------------------------------------
            k = 1
            tot = dt(k) 
            do while((tot .lt. tmin).and.(k.lt.ntau))
               k = k + 1
               tot = tot + dt(k)
            end do
cd            ntmn = max(k-1,1)
            ntmn = max(k,1)

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

               do k = 1, ntau
                  q(k) = q(k) +half*wgf*gf(k)*bf1
               end do

            else

               do i = 1, nu

c                 -------------------------------------
c                    Calculate W matrix (Diagonal)
c                 -------------------------------------
                  do k = 1, ntau
                     w(k) = wu(i)*gf(k)*wgf
                  end do 

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

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

c                 -------------------------------------
c                     Calculate T(-1)*B store in s
c                 -------------------------------------
                  do k = 1, ntmn - 1
                     s(k) = zero
                  end do
                  do k = ntmn, ntau
                     s(k) = bp(k)
                  end do
                  s(ntau) = s(ntau) + two*xu(i)*bf1/dt(ntau-1)

                  call lubksbt(ta, tb, tc, s, ntau)

c                 -----------------------------------
c                    Calculate W*T(-1)*S add into Q
c                 ------------------------------------
                  do k = 1, ntau
                     q(k) = q(k) + w(k) * s(k)
                  end do

               end do 

            end if 

         end do


      return 
      end