```
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
```