pro traj_read,file=file,theta=theta,lats=latsave,lons=lonsave,sza=szasave,temp=temsave,$
   time=time,pv=qsave,press=psave,ntsteps=nsteps,nparcels=nparcels,back=back,$
   hemisphere=qm,year=year,basicfacts=basicfacts,revdimensions=revdimensions,$
   scalartag=scalartag,err=err,parcel_no=parcel_no,prec_times=prec_times, $
   contour_ind=contour_ind,nosource=nosource
;+
; NAME: traj_read
;
; PURPOSE: Read in trajectory xdr type files
;
; CATEGORY: trajectory 
;
; CALLING SEQUENCE:
;   traj_read,file=file,theta=theta,lats=lats,lons=lons,sza=sza,temp=temp,$
;   time=time,pv=pv,press=press,ntsteps=ntsteps,nparcels=nparcels,back=back,$
;   hemisphere=qm,year=year,basicfacts=basicfacts,contour_ind=contour_ind,$
;   scalartag=scalartag,prec_times=prec_times,source=source,nosource=nosource
;
;
; INPUTS:
;    file name without .xdr appendix 
;
; OPTIONAL INPUT PARAMETERS:
;
; KEYWORD PARAMETERS:
;    basicfacts = keyword switch which has program read in only the
;       beginning information such as the theta, nsteps, nparcels etc.
;       in the file - normal reads may be much longer
;    revdimesion= kewword switch which reverses the dimensions as
;       given below 
;    scalartag  an array defined by the user with nparcel elements.
;      to get scalartag you must give a value in the call list
;       e.g.
;       scalartag=1
;       traj_read,.....,scalartag=scalartag
;       (scalartag will now be read and returned as an array)
;     err is nonzero if an error occurred 
;     parcel_no= parcel_no , if parcel_no is specified then only the 
;        data for a singel parcel is returned
;     prec_times = additional array used if times for trajectories are offset
;     contour_ind = contour index for each parcel - used if multiple contours
;     source = source character string for the data source
;     nosource = keyword for old style files
;        In the older files the source variable was written instead of 
;         the scalar tag varaible. If you get the error message:
;              % Unable to allocate memory: ensuring string length.
;         Then set the keyword /nosource and rerun the program
; KEYWORD OUTPUTS:
;            
;
;            ntsteps = number of time steps
;            nparcels = number of parcels
;            hemisphere = +1 for NH = -1 for SH
;            year = year number eg 89
;            theta = potential temperature
;            back = 1 for foreward trajectory, -1 for back trajectories
;            time= fltarr(ntsteps) time in days (using daynumbers)
;            lats= fltarr(ntsteps,nparcels), latitudes
;            lons= "", longitudes
;            temp= "", temperatures (K)
;            sza=  "", solar zenith angles
;            pv=   "" , potential vorticity
;            press="", pressure
;            scalartag = an array defined by the user with nparcels elements. This
;                user defined array was created when the trajectory model was
;                called
;            prec_times = times for individual parcels if they are
;                not sequenced on the same time track.
;            contour_ind = contour index for parcels
;            source = string (3chr) for source 
;    
;
; OPTIONAL OUTPUT PARAMETERS:
;
;
; COMMON BLOCKS:
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
;
; REQUIRED ROUTINES:
;
; MODIFICATION HISTORY: 
; MRS 9/1990
;-

get_lun,unit
err=0
openr,unit,file+'.xdr',/xdr,error=err
if err ne 0 then begin
   print,'From traj_read: File '+file+'.xdr not found'
   free_lun,unit
   return
endif
tlen=1
npts=long(1)
qm=1
back=1
theta=0.
yr=0
      on_ioerror,nosta
forrd,unit,npts,tlen,theta,qm,back,yr
goto,rxx
nosta:
   npts=fix(1)
   forrd,unit,npts,tlen,theta,qm,back,yr
rxx:
   year=yr
   nparcels=npts
   nsteps=tlen
   if keyword_set(basicfacts) then  begin
     free_lun,unit
     return
endif
latsave=fltarr(npts,tlen) & lonsave=latsave
qsave=latsave & szasave=latsave
temsave=latsave & psave=latsave & prec_times=latsave
time=fltarr(tlen)
scalartag=findgen(npts)
source='NMC'
contour_ind=fix(latsave)

forrd,unit,latsave,lonsave,qsave,temsave,time,szasave,psave
if not keyword_set(nosource)then begin 
      on_ioerror,nost
      forrd,unit,source,scalartag,contour_ind,prec_times
   endif else begin
      on_ioerror,nost
      forrd,unit,scalartag
      source='NMC'
   endelse
nost:

free_lun,unit
if n_elements(parcel_no) ne 0 then begin
   if parcel_no gt npts-1 then begin
     print,'Requested parcel_no '+str(parcel_no)+' ; maximum number = '+str(npts)
     err=1 
     return
     endif
 latsave = reform(latsave(parcel_no,*))
 lonsave = reform(lonsave(parcel_no,*))
 qsave = reform(qsave(parcel_no,*))
 temsave = reform(temsave(parcel_no,*))
 szasave = reform(szasave(parcel_no,*))
 psave = reform(psave(parcel_no,*))
 scalartag=scalartag(parcel_no)
 contour_ind=reform(contour_ind(parcel_no,*))
 prec_times=reform(prec_times(parcel_no,*))
 return
  endif
if keyword_set(revdimensions) then return

latsave=(transpose(latsave))
lonsave=(transpose(lonsave))
szasave=(transpose(szasave))
temsave=(transpose(temsave))
psave=(transpose(psave))
qsave=(transpose(qsave))
contour_ind=(transpose(contour_ind))
prec_times=(transpose(prec_times))

return
end