;+
; NAME:
;  rt_szss96_atmread
; PURPOSE: (one line)
;  read a model atmosphere
; DESCRIPTION:
;  read a model atmosphere, return atmospheric structure
; CATEGORY:
;  RT
; CALLING SEQUENCE:
;  atm = rt_szss96_atmread(fname, pmax)
; INPUTS:
;  fname = file name of atmospheric structure
;  pmax = maximum pressure (ubar) 
; 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 an atmosphere structure, with the following fields
;    nz  number of atmospheric levels
;    ns  number of species
;    z   altitudes [nz]  (cm)
;    p   pressures [nz] (cm)
;    t   temperatures [nz] (K)
;    n   number density [nz,ns] (cm^-3)
;    ncol   number column density [nz,ns] (cm^-2)
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
;  Modifies atm
; PROCEDURE:
; 
; EXAMPLE:
;
;  atm = rt_szss96_atmread('pluto.restart', 3.d)
; 
; MODIFICATION HISTORY:
;  Written 2004 May 8, Leslie Young SwRI
;  based on atmread2.f, written by Roger Yelle, 1990
;-

function rt_szss96_atmread, fname, pmax, codetype = codetype

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

    case codetype of
        'FORTRAN': begin
            mode = 2L
            init_atm,ZM,TEM,PRE,RHO,RHO1,RHO2,DN,DN1,BNS,BNS1,KM,MODE
            nz = KM
            z = ZM * 100.d
            p = PRE * 10.d
            t = TEM
            rch4 = rho1 / (16.0D0/28.0D0) ; volume mixing ratio
            rco  = rho2 / (28.0D0/28.0D0) ; volume mixing ratio
            ns = 3
            species = [22, 6, 5] ; species codes from HITRAN
            ntot = dn * 1d-6
            n = dblarr(nz,ns)
            n[*,1] = dn1 * 1d-6
            n[*,2] = ntot * (rco)
            n[*,0] = ntot - (n[*,1] + n[*,0])
            ncoltot = bns * 1d-4
            ncol = dblarr(nz,ns)
            ncol[*,1] = bns1 * 1d-4
            ncol[*,2] = ncoltot * (rco)
            ncol[*,0] = ncoltot - (ncol[*,1] + ncol[*,2])
            return, {nz:nz, $
                     ns:ns, $
                     species:species, $
                     z:z , $
                     p:p, $
                     t:t, $
                     ntot:ntot, $
                     n:n,  $
                     ncoltot:ncoltot, $
                     ncol: ncol }
        end
        'FORTEST': begin
            DELZ0_cm =1.0D5     ; delta z of the lowest layer
            FAC1 = 1.03D0
            nz = 106
            
            z = rt_szss96_zgrid(nz, delz0_cm, fac1)
            t = rt_szss96_told(fn, nz)
            p = rt_szss96_p(z, t)

            return, {nz:nz, $
                     z:z }
        end
        else: begin
            print, fcn, ': code type ', codetype, ' not implemented'
            return, -1
        end
    endcase

end