pro rt_cascadedriver

    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'


        if (iverbose ge 2) then print, 'Methane cascade heating'

        im = 0
        ib = 0
        ep1 = rt_epsln(atm.t, atm.n[*,0],$
                       band.arad[ib,im], band.pcol[ib,im],$
                       band.rdmol[im], band.rmcol[im])
        ib=3
        ep2 = rt_epsln(atm.t, atm.n[*,0],$
                       band.arad[ib,im], band.pcol[ib,im],$
                       band.rdmol[im], band.rmcol[im])
        if (iverbose ge 2) then begin
            print, 'ep1'
            print, ep1[0:2], form='(3e11.3)'
            print, '..', ep1[nz-2:nz-1], form='(A11,2e11.3)'
            print, 'ep2'
            print, ep2[0:2], form='(3e11.3)'
            print, '..', ep2[nz-2:nz-1], form='(A11,2e11.3)'
; FORTRAN OUTPUT
;  ep1
;   0.710E-07  0.981E-07  0.136E-06
;          ..  0.524E+00  0.723E+00
; ep2
;   0.867E-04  0.120E-03  0.167E-03
;          ..  0.640E+03  0.883E+03

        endif

        ts = dblarr(atm.nz, 4)
        hcsc = rt_cascade(band.qrot[im],!phys.m_u*band.rmass[im], $
                          dat.rh,1306.d0,dat.wdel,xu,wu,dat.nu,$
                          atm,ep1,ep2,ssolar,ts)

        ts[*,im] = ts[*,im]/(sfund.stot[im]*4.*!pi)
        if (iverbose ge 2) then begin
            print, 'hcsc'
            print, hcsc[0:2], form='(3e11.3)'
            print, '..', hcsc[nz-2:nz-1], form='(A11,2e11.3)'
        endif

;  FORTRAN OUTPUT
;  0.895E-15  0.146E-14  0.239E-14
;         ..  0.776E-07  0.992E-07


  
end