subroutine irheat(tmin,stot,opp0,wpf,opf,dopf,bpf,dbpf,bf1
     &   ,bf2,nopf,tbase,dtdc,xu,wu,nu,tmp,prs,den,clm,ep,dep,ts,ntau
     &   ,oppsum,bpmean,bj,s,heat,p)

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

c       -------------------------------------------------------------
c                        Input arrays
c       -------------------------------------------------------------
        dimension xu(nd6),wu(nd6),tmp(nd1),prs(nd1),clm(nd1),ep(nd1)
     &   ,ts(nd1),den(nd1),dep(nd1),wpf(nd2),opf(nd1,nd2),dopf(nd1,nd2)
     &   ,bpf(nd1,nd2),dbpf(nd1,nd2),bf1(nd2),bf2(nd2)

c       -------------------------------------------------------------
c                          Output arrays
c       -------------------------------------------------------------
        dimension bpmean(nd1),bj(nd1),s(nd1),heat(nd1),p(nd1,nd1)
     &  , oppsum(nd1)

c       -------------------------------------------------------------
c                         Internal arrays
c       -------------------------------------------------------------
        dimension uf(nd1,nd2), sf(nd1,nd2), dtau(nd1), ds(nd1,nd2) 
     &  , dbpmean(nd1), sb(nd1), dj(nd1,nd1), bb1(nd1), bb2(nd1)
     &  , bpi(nd1), epa(nd1), epb(nd1), tsa(nd1)

        do n = 1, ntau
           do m = 1, ntau
              p(m,n) = zero
           end do
        end do
        
        call means(wpf,opf,dopf,bpf,dbpf,nopf,ntau,oppsum,bpmean
     &  ,dbpmean)

        do n = 1, ntau
           bpi(n) = one/bpmean(n)
           epa(n) = one/(one+ep(n))
           epb(n) = ep(n)/(one+ep(n))
           tsa(n) = ts(n)*bpi(n)*epa(n)
        end do

        do k = 1, nopf
           do n = 1, ntau
              uf(n,k) = bpf(n,k)*bpi(n)*epa(n)
              sf(n,k) = (tsa(n)+epb(n))*bpf(n,k)
           end do
        end do 
       
        do n = 1, ntau - 1
           dtau(n) = clm(n + 1) - clm(n)
        end do

        call rybicki(tmin, dtau, uf, sf, ntau, xu, wu, nu, opf, wpf 
     &   , oppsum, nopf, bf1, bf2, bj)

c       -----------------------------------------------------------
c                           calculate S/B 
c       -----------------------------------------------------------
        do n = 1, ntau
           s(n) =  (bj(n)+ts(n))*bpi(n)*epa(n) + epb(n)  
           sb(n) = s(n) * bpmean(n)
        end do

c       -----------------------------------------------------------
c                      calculate heating rate
c       -----------------------------------------------------------
        fourpi = four * pi
        do n = 1, ntau
           heat(n)=fourpi*oppsum(n)*den(n)*(bj(n)-sb(n))
        end do

c       ------------------------------------------------------------
c              calculate derivative of source function
c       ------------------------------------------------------------
        do n = 1, ntau
           bb1(n) = (one-bj(n)/bpmean(n))*dep(n)/(one+ep(n))**2
     &         - bj(n)*dbpmean(n)/bpmean(n)**2/(one+ep(n))
           bb2(n) = (ep(n)+bj(n)/bpmean(n))/(one+ep(n))
        end do

        do k = 1, nopf
           do n = 1, ntau
              ds(n,k) = bb1(n)*bpf(n,k)+bb2(n)*dbpf(n,k)
           end do
        end do

c       ------------------------------------------------------------
c          drybicki calculates the temp derivative of the mean 
c          intensity
c       ------------------------------------------------------------
        call drybicki(tmin,dtau,uf,sb,ds,ntau,xu,wu,nu,opf,dopf,wpf
     &   ,nopf,bf1,oppsum,dj)

c       ------------------------------------------------------------
c               temperature derivative of heating rate
c       ------------------------------------------------------------
        do n = 1, ntau
           epn = ep(n)
           xx = one/(one+epn)
           fpx = fourpi*den(n)*xx
           p(n,n) = fpx*oppsum(n)*(-epn*dbpmean(n)+xx*(bj(n)
     &      -bpmean(n))*dep(n))
           fpp = fpx*epn
           do m = 1, ntau
              p(n,m) = p(n,m)+fpp*dj(n,m)
           end do
        end do

        return
        end