;+
; NAME:
;  rt_szss96_pre
; PURPOSE: (one line)
;  Return the pressure
; DESCRIPTION:
;  Return the pressure
; CATEGORY:
;  RT
; CALLING SEQUENCE:
;  p = rt_szss96_pre(
; INPUTS:
;  r00 - Pluto's radius (cm)
;  t - Temperature (K) [array]
;  rch4 - volume mixing ratio of CH4 (n_ch4/n_tot)  [array or scalar]
; OPTIONAL INPUT PARAMETERS:
; KEYWORD INPUT PARAMETERS:
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;    t   Temperature [nz]  (K)
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
; 
; EXAMPLE:
;
;  t = rt_szss96_pre('pluto.restart', 106)
; 
; MODIFICATION HISTORY:
;  Written 2007 Aug 31, Leslie Young SwRI
;  based on plutot.f, written by Xun Zhu
;-

function rt_szss96_p, z, t, codetype = codetype

    fcn = 'rt_szss96_p'

    if not keyword_set(codetype) then codetype = 'FORTEST'
    
    case codetype of
        'FORTEST': begin
            ZM = Z/1.0D2       ; Z in m
            TEM = t
            P00 = 0.3d
            R00=1150.0D3        ; Pluto's radius in m
            KM = n_elements(z)
            HSCAL = dblarr(KM)
            PRE = dblarr(KM)
;            BC=1.3804D-23       ; Boltzmann constant: J K
            GG0=0.59D0          ; Surface gravity: m s**(-2)
            RCH4X0=0.030D0 ; CH4 volume mixing ratio at surface      ************
            WTMOL=28.0D0*(1.0D0-RCH4X0)+16.0D0*RCH4X0 ; mean molecular weight
            RGAS=8314.0D0/WTMOL ; gas constant for N2: J K**(-1) kg**(-1)
            GRV=GG0*(R00/(R00+ZM))^2 ; gravity
            HSCAL[1-1]=RGAS*TEM[1-1]/GG0
            PRE[1-1]=P00
            for K=2,KM do begin
                HSCAL[K-1]=RGAS*(TEM[K-1]+TEM[K-2])/(GRV[K-1]+GRV[K-2])
                FAC2=(ZM[K-1]-ZM[K-2])/HSCAL[K-1]
                PRE[K-1]=PRE[K-2]*EXP(-double(FAC2))
            endfor
;            for K=1,KM do begin
;                RHO[K-1]=PRE[K-1]/(RGAS*TEM[K-1]) ; mass density: kg m**(-3)
;                DN[K-1]=PRE[K-1]/(BC*TEM[K-1]) ; number density: m**(-3)
;            endfor
;            BN[KM-1]=DN[KM-1]*HSCAL[KM-1] ; column number density: m**(-2)
;            for K=2,KM do begin
;                IK=KM-K+1
;                FACC1=(ZM[IK]-ZM[IK-1])*(DN[IK]+DN[IK-1])/2.0D0
;                BN[IK-1]=BN[IK-1]+FACC1
;            endfor
            p = PRE * 10.d
        end
    endcase

    return, p

end