;+
; NAME:
;  rt_szss96_Ch4heat33
; PURPOSE: (one line)
;  calculate CH4 cooling at 7.7 micorn
; DESCRIPTION:
;  calculate CH4 cooling at 7.7 micorn
; CATEGORY:
;  RT
; CALLING SEQUENCE:
;  h33 = rt_szss96_Ch4heat33(t, p, n_ch4, gamn4=gamn4)
; INPUTS:
;  t - temperature (K)
;  p - pressure (microbar)
;  n_ch4 - number density of CH4 (cm^-3)
; OPTIONAL INPUT PARAMETERS:
; KEYWORD INPUT PARAMETERS:
;  codetype
;    'FORTRAN' - call routines from plutot.f
;    'FORTEST' - native IDL routines that reproduce the fortran
;    'IDL' - idl code, without the restriction of needing to reproduce fortran
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;  returns cooling rate in erg/cm^3/s
;  gamn4 : cool-to-space Gamma
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
; 
; EXAMPLE:
;
;  atm = rt_szss96_atmread('pluto.restart', 3.d)
;  t = atm.t  
;  p = atm.p
;  n_ch4 = atm.n[*,1]
;  l77 = rt_szss96_Ch4heat33(t, p, n_ch4)
; MODIFICATION HISTORY:
;  2008 Nov 9 - fixed documentation
;-

function rt_szss96_Ch4heat33, t, p, n_ch4, ncol_ch4, gamn4=gamn4

    fcn = 'rt_szss96__Ch4heat33'
    
    if not keyword_set(codetype) then codetype = 'FORTRAN'

    case codetype of
        'FORTRAN': begin

            szss96_get_ch4all,GG,WW,TREF
            if tref[0] ne 40.d then begin
                szss96_input6
            endif

            szss96_get_ch4v3, sigx
            if sigx[0,0] lt 1d-51 then begin
                szss96_input4
            endif

            szss96_get_parair, R00,WTMOL
            if WTMOL lt 1. then begin
                R00 = 1150.0D3
                WTMOL = 28.0D
                szss96_set_parair,R00, WTMOL
            endif

            ; convert from standard atm structure to 
            ; fortran's expected format
            AN = double(n_ch4 * 1d6)
            BN = double(ncol_ch4 * 1d4)
            TEM = double(t)
            PRE = double(p / 10.d)
            KM = n_elements(PRE)
            ICONT = 2L
            IYELLE = 2L
            IZHU = 1L
            

            if not keyword_set(gamn4) then begin
                BC = 1.3804D-23
                DN1 = n_ch4 * 1d6
                DN = PRE/(BC * TEM)
                RHO = double( ( DN1/ DN ) * (16/28.) )
                szss96_qtch4z,TEM,PRE,RHO,QT,KM,ICONT,GAMN4,IYELLE
            endif

            szss96_htch4,HT,AN,BN,TEM,PRE,KM,ICONT,IYELLE,GAMN4,IZHU

            ;convert from K/s to erg/cm^3/s
            RAIR=8314.0D0/WTMOL  ; gas constant
            CP=RAIR*7.0D0/2.0D0  ; specific heat at the constant pressure
            RHO = PRE/(RAIR*TEM) ; kg m**(-3)
            FACTX1=CP*RHO*10.0D0 ; To convert K/sec to erg cm**(-3) s**(-1)

            HT1CGS = HT * FACTX1

            ht = HT1CGS

            return, ht
        end
        'FORTEST': begin
            print, fcn, ': code type ', codetype, ' not implemented'
            return, -1
        end
        else: begin
            print, fcn, ': code type ', codetype, ' not implemented'
            return, -1
        end
    endcase

end