subroutine conduct(rplnt,tbot,tflux,zt,tmp,ntau,htcn,p)

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

        dimension zt(nd1),tmp(nd1),htcn(nd1),p(nd1,nd1)

           jt = 1
           cnb = fconduct(tmp(jt))
           cnc = fconduct(tmp(jt+1))
           cnp = half*(cnb+cnc)
           dcb = fdconduct(tmp(jt))
           dcc = fdconduct(tmp(jt+1))
           rb = rplnt + zt(jt)
           rc = rplnt + zt(jt+1)
           rp = half*(rb+rc)
           ra = two*rb - rc
           rm = half*(ra+rc)
           drp = rc - rb
           drm = rb - ra
           dr = half*(drp+drm)
           tc = (cnp/dr/drp)*(rp/rb)**2
           htcn(jt)=(tflux/dr)*(rm/rb)**2
     $      + tc*(tmp(jt+1)-tmp(jt))

           sb = tc*(one+half*(cnb/cnp)*dcb)
           sc = tc*(one+half*(cnc/cnp)*dcc)

           p(jt,jt)= - sb
           p(jt,jt+1)= sc
           
           do jt = 2, ntau-1
              cna = cnb
              cnb = cnc
              cnc = fconduct(tmp(jt+1))
              cnm = cnp
              cnp = half*(cnb+cnc)
              dca = dcb
              dcb = dcc
              dcc = fdconduct(tmp(jt+1))
              ra = rb
              rb = rc
              rc = rplnt + zt(jt+1)
              rm = rp
              rp = half*(rb + rc)
              drm = drp
              drp = rc - rb
              dr = half*(drm + drp)
              ta = (cnm/dr/drm)*(rm/rb)**2
              tc = (cnp/dr/drp)*(rp/rb)**2
              htcn(jt) = ta*(tmp(jt-1)-tmp(jt))+tc*(tmp(jt+1)-tmp(jt))

              sa = ta*(one+half*(cna/cnm)*dca)
              sb = ta*(one+half*(cnb/cnm)*dcb)
     $           + tc*(one+half*(cnb/cnp)*dcb)
              sc = tc*(one+half*(cnc/cnp)*dcc)

              p(jt,jt-1)= sa
              p(jt,jt)= - sb
              p(jt,jt+1)= sc

           enddo

           jt = ntau

           ra = rb
           rb = rc
           rc = two*rb - ra
           rm = rp
           rp = half*(rb + rc)
           drm = drp
           drp = rc - rb
           dr = half*(drm + drp)

           cna = cnb
           cnb = cnc
           cnc = fconduct(tbot)
           cnm = cnp
           cnp = half*(cnb+cnc)
           dca = dcb
           dcb = dcc
           dcc = fdconduct(tbot)
           ta = (cnm/dr/drm)*(rm/rb)**2
           tc = (cnp/dr/drp)*(rp/rb)**2
           htcn(jt) = ta*(tmp(jt-1)-tmp(jt))+tc*(tbot-tmp(jt))

           sa = ta*(one+half*(cna/cnm)*dca)
           sb = ta*(one+half*(cnb/cnm)*dcb)
     $        + tc*(one+half*(cnb/cnp)*dcb)
         
           p(jt,jt-1) =  sa
           p(jt,jt) = - sb

c           htcn(jt) = zero
c           p(jt,jt-1) = zero
c           p(jt,jt) = one

        return
        end

        function fconduct(temp)
        implicit double precision (a-h,o-z)
        parameter (a1=42.7d0,a2=2.61d0,a3=355.6d0,s=0.69d0)
        fconduct = a1*temp**s + a2*temp - a3   
        return
        end

        function fdconduct(temp)
        implicit double precision (a-h,o-z)
        parameter (a1=42.7d0,a2=2.61d0,s=0.69d0)
        fdconduct = (a1*s*temp**s + a2*temp)/fconduct(temp) 
        return
        end