pro convert_saber_netcdf_to_ascii, input_file_name, output_directory
;pro convert_saber_netcdf_to_ascii,input_file_name
;
; this routine converts a SABER netcdf file to an ASCII
; file that can be read by Maya's comparison code
;
read_netcdf,input_file_name,data,attributes,status
if status ne 0 then begin
  if status eq -1 then print, "BAD_PARAMS"
  if status eq -2 then print, "BAD_FILE"
  if status eq -3 then print, "BAD_DATA_FILE"
  if status eq -4 then print, "FILE_ALREADY_OPENED"
  stop
endif
this_event=0
print,data.event
read, 'Select Event from Event list> ', this_event
;
;             get file type, level2,level2a, level2b
;
file_type=strmid(attributes(1),18,7)
;
;             get the year and doy, then remove blanks
;
; first find the event that matches the requested event
; and assign the index number for that event
;
this_event_index=where(data.event eq this_event)
if this_event_index lt 0 then begin
  print,'*** ', this_event, ' not a valid event   ***'
  stop
endif
year_doy=string(data(this_event_index).date)
year_doy=strmid(year_doy,5,7)
;
;             now get orbit number and software version
;
orbit_number=strmid(attributes(15),18,5)
software_version=strmid(attributes(6),24,4)
;
;              now combine these into a string that will be the
;              first comment line in the ASCII file
;
first_ascii_line="* SABER atm profile from "+file_type+" "+orbit_number+" "+year_doy+" "+software_version
;
;              now the 2nd comment line
;              get the current date
;
current_date=systime()
second_ascii_line="* created  : "+current_date
;
;              and now the 3rd comment line
;              create event number as string 
;
ascii_event=string(this_event)
third_ascii_line="* Date: "+year_doy+"; Event:  "+ascii_event
;
;               grab data
;
altitudes=data(this_event_index).tpaltitude
ktemp=data(this_event_index).ktemp
pressure=data(this_event_index).pressure
o3_96=data(this_event_index).o3_96
h2o=data(this_event_index).h2o

;
;               We need a single latitude, longitude, and sza value for the event 
;               to put in the file header. We will pick values near 90km. 
;               Determine the data point number near 90km
;               to use in our selection of latitude, longitude and sza
;
size_it=size(altitudes)-1
for i=0,size_it(1) do begin
  if altitudes(i) gt 90.0 then counter=i
endfor
;
;                now the 4th comment line
;                we also need to grab the day/night flag 
;                in addition to the values at 90km
;
day_night=data(this_event_index).tpDN
tplatitude_90km=data(this_event_index).tplatitude(counter)
tplongitude_90km=data(this_event_index).tplongitude(counter)
tpSolarZen_90km=data(this_event_index).tpSolarZen(counter)
time_since_midnight_90km=data(this_event_index).time(counter)
int_tplongitude_90km=fix(tplongitude_90km)
int_tplatitude_90km=fix(tplatitude_90km)
if int_tplatitude_90km gt 0 then latitude_string=strtrim(string(int_tplatitude_90km),1)+"N"
if int_tplatitude_90km lt 0 then latitude_string=strtrim(string(int_tplatitude_90km),1)+"S"
int_tpSolarZen_90km=fix(tpSolarZen_90km)
fourth_ascii_line="* Latitude: " +latitude_string+$
                  "; Longitude: "+strtrim(string(int_tplongitude_90km),1)+ $
                  "; SZA  "      +strtrim(string(int_tpSolarZen_90km),1)+$
                  "(DN: "        +strtrim(string(day_night),1)+")"
;
;                 remove missing (-999) data
;
;screened_counter=where(altitudes ne -999.00)
screened_counter=where(ktemp ne -999.00)
altitudes_screened=altitudes(screened_counter)
ktemp_screened=ktemp(screened_counter)
pressure_screened=pressure(screened_counter)
o3_96_screened=o3_96(screened_counter)
h2o_screened=h2o(screened_counter)

size_it=size(altitudes_screened)
size_of_screened=size_it(1)
;
;                   create file name for output
;                   should be of the form saber_PT_2002135_02347_e033.dat_v1.08
;                   under a directory named of the form 2002135
;
zero='0'
event_as_string=strtrim(string(this_event),2)
if this_event le 99 then event_as_string=zero+strtrim(string(this_event),2)
if this_event le 9 then event_as_string=zero+zero+strtrim(string(this_event),2)
output_file_name=output_directory+"/"+year_doy+"/saber_PT_"+ year_doy+"_"+orbit_number+"_e"+event_as_string+".dat_v"+software_version
;output_file_name="saber_ascii_data/"+year_doy+"/saber_PT_"+ year_doy+"_"+orbit_number+"_e0"+event_as_string+".dat_v"+software_version
openw,1,output_file_name
;
;                   write out first 4 lines which are comments
;
printf,1,first_ascii_line
printf,1,second_ascii_line
printf,1,third_ascii_line
printf,1,fourth_ascii_line
;
;                   the values in some of the next lines may very well be read!
;
printf,1,size_of_screened,"  -> number of lines"
printf,1,format='(F7.2,"  -> latitude")',  tplatitude_90km
printf,1,format='(F7.2,"  -> SZA")', tpSolarZen_90km 
printf,1,format='(F7.2,"  -> longitude")', tplongitude_90km
printf,1,format='(6x,I6,  "  -> seconds since midnight")',   time_since_midnight_90km/1000
printf,1,"* (SABER data without processing)"
printf,1,"*   Alt         T         Pressure        O3          H2O"
printf,1,"*   (km)       (K)          (mb)        (ppmv)       (ppmv)"
;
;                   now write out the data starting at low altitudes and going up
;                   note: mixing ratios converted to ppmv
;
doloopcounter=size_of_screened-1
for i=doloopcounter,0, -1 do begin
   printf,1,format='(2x,F6.2,E13.5,E13.5,E13.5,E13.5)',$
   altitudes_screened(i),ktemp_screened(i),pressure_screened(i),$
   o3_96_screened(i)*1000000.,h2o_screened(i)*1000000.
endfor
close,1
end