subroutine radioheat(tmin,wav,str,eng,wid,nlines,tbase,tmp
     &   ,den,clm,prs,ntau,xu,wu,nu,ht,p)

        implicit double precision (a-h,o-z)
        include "parameters.dat"
        parameter (tmp0 = 296.d0,rmass=4.48d-23,tref=150.d0
     &   ,nlortz=300,wlortz=1.2d-7,strmin=1.d-22)
        dimension wav(1),str(1),eng(1),wid(1),tmp(1),den(1)
     &   ,clm(1),prs(1),xu(1),wu(1),ht(1),p(nd1,1)
        dimension f(nd1,nd2),df(nd1,nd2),bp(nd1),dbp(nd1)
     &   ,w(nd2),x(nd2),strn(nd1),dj(nd1,nd1),dstrn(nd1)
     &   ,hckt(nd1),uc(nd1),qrot(nd1),dcol(nd1),h(nd1)

        uprm = dsqrt(two*rkb/rmass)/cs
        wdopref0 = uprm*dsqrt(tref)
        nmax = 1000

        do n = 1, ntau
           do m = 1, ntau
              p(m,n) = zero
           end do
           ht(n) = zero
        end do

        do n = 1, ntau
           hckt(n) = hck*(one/tmp(n)-one/tmp0)
           qrot(n) = tmp0/tmp(n)
           uc(n) = uprm*dsqrt(tmp(n))
        end do

        do n = 1, ntau-1
           dcol(n) = clm(n+1)-clm(n) 
        end do

        tmpmax = 0.d0
        tmpmin = 10000.d0
        imin = 1
        imax = 1
        do n = 1, ntau
           if(tmp(n) .lt. tmpmin) then
              tmpmin = tmp(n)
              imin = n
           end if
           if(tmp(n) .gt. tmpmax) then
              tmpmax = tmp(n)
              imax = n
           endif
        end do

        do lw = 1, nlines

c           write(unit=*,fmt='(i5,f10.6,1x,e10.3,1x,f11.6)') 
c     &     lw,wav(lw),str(lw),eng(lw)

           str1 = qrot(imin)*str(lw)*dexp(-hckt(imin)*eng(lw))
           str2 = qrot(imax)*str(lw)*dexp(-hckt(imax)*eng(lw))
           if((str1 .gt. strmin).or.(str2 .gt. strmin)) then

              wavlw = wav(lw)
              wdopref = wdopref0*wavlw
              wstep = wdopref
              fmin = 1.d-2*strmin/wdopref
              bb = bplanck(wavlw,tbase)
              hckw = hck*wavlw
              hckwt0 = hckw/tmp0
              strn0 = str(lw)/(one-dexp(-hckwt0))

              do n = ntau, 1, -1
                 bp(n) = bplanck(wavlw,tmp(n))
                 dbp(n) = bp(n)*dplanck(wavlw,tmp(n))/tmp(n)
                 wdop = wavlw*uc(n)
                 hckwt = hckw/tmp(n)
                 strn(n) = qrot(n)*strn0*dexp(-hckt(n)*eng(lw))
     &            *(one-dexp(-hckwt))
                 alor=1.d-6*prs(n)*wid(lw)*dsqrt(tmp(n)/tmp0)/wdop
                 dstrn(n)=hck*eng(lw)/tmp(n)-hckwt*dexp(-hckwt)
     &            /(one-dexp(-hckwt))-one
                 fac = strn(n)/sqpi/wdop

                 if(n .eq. ntau) then
                    wlast = zero 
                    i = 0
 10                 continue
                       wnext = wlast+wstep
                       i = i + 1
                       x(i) = half*(wnext+wlast)
                       w(i) = two*(wnext-wlast)
                       call voigt(x(i)/wdop,alor,vgt)
                       f(n,i) = fac*vgt
                       df(n,i) = (f(n,i)/tmp(n))*dstrn(n)
                       wlast = wnext
                       if(i .gt. 1) then
                          if(f(n,i)/f(n,i-1) .gt. 0.8d0) 
     &                     wstep = 1.2d0*wstep
                       end if
c                       write(unit=*,fmt=810) i,x(i),w(i),f(n,i)
                    if((i .lt. nmax).and.(f(n,i).gt.fmin)) goto 10
                    iwx = i 
c                    read(unit=*,fmt='(a)') icr
                 else
                    do i = 1, iwx
                       call voigt(x(i)/wdop,alor,vgt)
                       f(n,i) = fac*vgt
                       df(n,i) = (f(n,i)/tmp(n))*dstrn(n)
                    end do
                 end if
              end do

              call drybicki_lte(tmin,dcol,bp,dbp,ntau,xu,wu,nu,f
     &         ,df,w,iwx,bb,h,dj)

              do n = 1, ntau
                 fac=four*pi*den(n)
                 ht(n)=ht(n)+fac*h(n)
                 do m = 1, ntau
                    p(n,m)=p(n,m)+fac*dj(n,m)
                 end do
              end do

           end if

        end do

c        do n = 1, ntau
c           write(unit=*,fmt='(i5,7e11.3)')
c     &      n,ht(n),dpt(n),(dj(m,n),m=1,ntau)
c        end do
c        read (unit=*,fmt='(a)') icr

 810    format(i6,' x = ',e11.3,' w = ',e11.3,' f = ',e11.3)
        return
        end