;+
; NAME:
;  rt_szss96_Cocool
; PURPOSE: (one line)
;  calculate CO cooling
; DESCRIPTION:
;  calculate CO cooling
; CATEGORY:
;  RT
; CALLING SEQUENCE:
;  q77 = rt_szss96_Cocool(t, p, rho_co)
; INPUTS:
;  t - temperature (K)
;  p - pressure (microbar)
;  rho_co - fraction (by weight) of CO (unitless)
; 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
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
; 
; EXAMPLE:
;
;  atm = rt_szss96_atmread('pluto.restart', 3.d)
;  t = atm.t  
;  p = atm.p
;  n_co = atm.n[*,1]
;  l77 = rt_szss96_Cocool(t, p, n_co)
; MODIFICATION HISTORY:
;  2008 Nov 9 - fixed documentation
;-

function rt_szss96_cocool, t, p, rho_co

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

    case codetype of
        'FORTRAN': begin

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

            szss96_get_gam212, x212, gam212
            if x212[0] eq 0. then szss96_init_gam212


            szss96_get_freng, vmco, strco, eneco
            if vmco[0] eq 0. then szss96_init_freng

            ; convert from standard atm structure to 
            ; fortran's expected format
            TEM = double(t)
            PRE = double(p / 10.d)
            RHO2 = double(rho_co)
;            BC = 1.3804D-23
;            DN1 = n_co * 1d6
;            DN = PRE/(BC * TEM)
;            RHO = double( ( DN1/ DN ) * (16/28.) )
            KM = n_elements(PRE)
            ICONT = 2L
            IYELLE = 2L

            szss96_qtcoz,TEM,PRE,RHO2,QT2,KM,ICONT

            ;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)

            QT2CGS = QT2 * FACTX1

            qt = QT2CGS

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

end