subroutine rtefe2(dtau,s,rup,rdown,ri)
C     ======================================
C
C     formal solver of the radiative transfer equation 
C     for one frequency, angle, and for completely known source function;
c     original Feautrier (second-order) scheme
c
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      parameter (one=1.d0)
      dimension dtau(mdepth),s(mdepth),ri(mdepth)
      dimension a(mdepth),b(mdepth),c(mdepth),d(mdepth),
     *          f(mdepth),v(mdepth),z(mdepth)
c
c     set up the global tridiagonal matrix
c               
c     upper boundary condition 
c                 
      id=1
      cc=two/dtau(id)
      c(id)=cc/dtau(id)
      b(id)=one+cc+c(id)
      a(id)=0.
      v(id)=s(id)+cc*rup
c               
c     normal depth points 
c              
      do id=2,nd-1
         dtinv=two/(dtau(id-1)+dtau(id))  
         a(id)=dtinv/dtau(id-1)
         c(id)=dtinv/dtau(id)  
         b(id)=one+a(id)+c(id)
         v(id)=s(id)
      enddo
c               
c     lower boundary condition 
c              
      id=nd
      aa=two/dtau(id-1)
      a(id)=aa/dtau(id-1)
      b(id)=one+aa+a(id)
      if(rdown.eq.0.) b(id)=one+a(id)
      c(id)=0.
      v(id)=s(id)+aa*rdown
c               
c     ---------------------------------------------------
c     solution by elimination
c     1. forward sweep
c               
      f(1)=(b(1)-c(1))/c(1)
      d(1)=one/(one+f(1))
      z(1)=v(1)/b(1)
c 
c     ii) normal depth points
c
      do id=2,nd-1
         f(id)=(b(id)-a(id)-c(id)+a(id)*f(id-1)*d(id-1))/c(id)
         d(id)=one/(one+f(id))
         z(id)=(v(id)+a(id)*z(id-1))*d(id)/c(id)
      enddo
c
c     iii) upper boundary
c
      id=nd
      z(id)=(v(id)+a(id)*z(id-1))/(b(id)-a(id)*d(id-1))
c
c     2. backward elimination 
c   
      ri(nd)=z(nd)
      do id=nd-1,1,-1
         ri(id)=ri(id+1)*d(id)+z(id)
      enddo
c
      return
      end