c=======================================================================


      subroutine thflux(lev,nz,taug,bplnk,fup,fdn),2
c
c compute upward and downward thermal flux at level lev
c both surface reflectance and atmospheric scattering are 
c assumed negligible.  The transmission is assumed to be given
c by exp(-taug).
c
      dimension taug(*),bplnk(*)
      parameter (mxly=50)
      parameter (pi=3.141592654)
      dimension e3(mxly),e4(mxly)
c
      recurs(x,f)=(exp(-x)-x*f)/3
c
      e3(lev)=1./2.
      e4(lev)=1./3.
c
c upward flux (contribution from below level LEV)
c
      do 10 i=lev+1,nz
        t=taug(i)-taug(lev)
        e3(i)=expint(3,t)
        e4(i)=recurs(t,e3(i))
 10   continue
c 
      fup=bplnk(nz)*e3(nz)

cdb1  write(*,'(3x,7a11)') 't','dt','e3','e4','slope','fup','pi*bplnk'
      do 20 i=lev,nz-1
        dt=taug(i+1)-taug(i)
        if(dt.eq.0.) goto 20
        slope=(bplnk(i+1)-bplnk(i))/dt
        fup=fup+e3(i)*bplnk(i)-e3(i+1)*bplnk(i+1)
     &         +slope*(e4(i)-e4(i+1))
cdb1    write(*,'(i3,1p7e11.3)') i,taug(i)-taug(1),dt,
cdb1 &                e3(i),e4(i),slope,2*pi*fup,pi*bplnk(i)
 20   continue        
cdb1  write(*,'(i3,1p7e11.3)') nz,0.,0.,e3(nz),e4(nz),
cdb1 &                0.,2*pi*fdn,pi*bplnk(nz)
      fup=2*pi*fup
c
c downward flux (contribution from above level LEV)
c
      do 30 i=1,lev-1
        t=taug(lev)-taug(i)
        e3(i)=expint(3,t)
        e4(i)=recurs(t,e3(i))
 30   continue

      fdn=0.
cdb2  write(*,'(3x,6a11)') 't','dt','e3','e4','slope','fdn','pi*bplnk'
      do 40 i=1,lev-1
        dt=taug(i+1)-taug(i)
        if(dt.eq.0.) goto 40
        slope=(bplnk(i+1)-bplnk(i))/dt
        fdn=fdn+e3(i+1)*bplnk(i+1)-e3(i)*bplnk(i)
     &         -slope*(e4(i+1)-e4(i))
cdb2  write(*,'(i3,1p7e11.3)') i,taug(nz)-taug(i),dt,
cdb2 &                e3(i),e4(i),slope,2*pi*fdn,pi*bplnk(i)
 40   continue
cdb2  write(*,'(i3,1p7e11.3)') nz,0.,0.,e3(nz),e4(nz),
cdb2 &                         0,2*pi*fdn,pi*bplnk(nz)
      fdn=2*pi*fdn
      return
      end