pro rt_irheatdriver

    iverbose = 2
    physconstants

    restore, 'band.sav'
    restore,'atm.sav'
    nz = atm.nz
    restore,'dat.sav'
    gauleg, 0.d, 1.d, xu, wu, dat.nu
    restore,'sfund.sav'
;    restore,'ssolar.sav'


    imol = 0

; -----  snip 

    ib = 0
    ep = rt_epsln(atm.t, atm.n[*,0],$
                  band.arad[ib,imol], band.pcol[ib,imol],$
                  band.rdmol[imol], band.rmcol[imol])
    dep = rt_depsln(atm.t, atm.n[*,0],$
                    band.arad[ib,imol], band.pcol[ib,imol],$
                    band.rdmol[imol], band.rmcol[imol]) * (ep/atm.t)
    if (iverbose ge 2) then begin
        print, 'ep'
        print, ep[0:2], form='(3e11.3)'
        print, '..', ep[nz-2:nz-1], form='(A11,2e11.3)'
        print, 'dep'
        print, dep[0:2], form='(3e11.3)'
        print, '..', dep[nz-2:nz-1], form='(A11,2e11.3)'
    endif
; FORTRAN
; ep
;  0.710E-07  0.981E-07  0.136E-06
;         ..  0.524E+00  0.723E+00
; dep
;  0.195E-09  0.271E-09  0.380E-09
;         ..  0.149E-02  0.206E-02

    if (fmix[imol] gt 0.d) then do begin

;           call absorbk(imol,nlortz,wlortz(imol),qrot(imol)
;     &     ,wav(nl(imol)),str(nl(imol)),eng(nl(imol)),wid(nl(imol))
;     &     ,nlines(imol),otol,amu*rmass(imol),wstart(imol)
;     &     ,wend(imol),wdel,wstep(imol),wdopref(imol),opp0(imol)
;     &     ,tbase,tmp,prs,clm(1,imol),ntau,wpf,opf,dopf,nopf,bpf
;     &     ,dbpf,bf1,bf2)

;           call abscheck(wpf,opf,dopf,bpf,dbpf,bf1,bf2,nopf,ntau)

;           call irheat(tmin,stot(imol),opp0(imol),wpf,opf,dopf,bpf
;     &      ,dbpf,bf1,bf2,nopf,tbase,dtdc(imol),xu,wu,nu,tmp,prs
;     &      ,den(1,imol),clm(1,imol),ep,dep,ts(1,imol),ntau
;     &      ,oppsum(1,imol),bpmean(1,imol),bj(1,imol),s(1,imol)
;     &      ,hir(1,imol),pir(1,1,imol))

    endif                       ; if fmix
    
; -----  snip 
    


  
end