PRO read_atox_athy_profiles,pfile

;  sample read software for SABER updated nightime atomic oxygen and atomic hydrogene profile netCDF data files
;  provides examples of reading several parameters and creating a plot of atomic oxygen
;  data selected using several pieces of metadata as search criteria

;
;  INPUTS:
;  pfile - the full directory path and file name of the SABER atomic oxygen file to be read
;  File names are of the form
;  atox_athy_night_YYyyyy_Vv.nc
;  where y is the four-digit yaw year like 2004 and
;  where v is the product version number like 1.0 or 1.02
;  
;  Specify month, pressure level, and latitude range in the code below at ***
;
;  NOTE that the data for a given calendar day includes those orbits that begin on that day.
;  The last orbit that begins before midnight usually extends into the next day, but all of its
;  scans are included for the day in which the orbit originates.

;  *** edit data selection criteria here
;  NOTE:  data are provided on the following pressure level grid
; 0.0100000
; 0.00794328
; 0.00630957
; 0.00501187
; 0.00398107
; 0.00316227
; 0.00251188
; 0.00199526
; 0.00158489
; 0.00125892
; 0.00100000
; 0.000794328
; 0.000630957
; 0.000501187
; 0.000398107
; 0.000316227

mon = 3
pres = 1.e-3
min_lat = -45.0
max_lat = -40.0

min_lat_str = lat_to_latstr(min_lat,'(F5.1)')
max_lat_str = lat_to_latstr(max_lat,'(F5.1)')

months = ['January','February','March','April','May','June','July','August','September','October','November','December']
monthdays = [0,31,28,31,30,31,30,31,31,30,31,30]

;  make sure the file exists
IF FILE_TEST(pfile) NE 1 THEN BEGIN
  PRINT,'File not found: ',pfile
  RETURN
ENDIF

;  get the directory path for the file and the file name without suffix
fdir = FILE_DIRNAME(pfile)
fname = FILE_BASENAME(pfile,'.nc')

;   plot atomic oxygen for the month of March at a specific pressure
;   level within +/- 5 degrees of the equator (values hard-coded below)

;  get the year from the file name
fparts = STRSPLIT(fname,'_',/EXTRACT)
yr = FIX(STRMID(fparts[3],2,4))

; build the plot title from the metadata specified
title = 'Nighttime Atomic Oxygen at Pressure ' + STRING(pres,FORMAT='(F5.3)') +  ' hPa ' + $
  'in '+ months[mon-1] + ' ' + STRING(yr,FORMAT='(I4)') + '  Latitude range: ' + $
  min_lat_str + ' to ' + max_lat_str

;  open the NetCDF file
fid = NCDF_OPEN(pfile,/NOWRITE)

;  read the data for the necessary variables from the file
NCDF_VARGET,fid,'pressure',pressure
NCDF_VARGET,fid,'year',year
NCDF_VARGET,fid,'day',day
NCDF_VARGET,fid,'time',time
NCDF_VARGET,fid,'lat',lat
NCDF_VARGET,fid,'qatox',atox
numscans = N_ELEMENTS(time)

;  convert the year, day of year, and time to Julian dates
jdates = DBLARR(numscans)
FOR t=0,numscans-1 DO BEGIN
  doy_to_mdy,day[t],year[t],mm,dd
  ; convert the SABER time value of "msec since midnight" to floating point
  ; hours since midnight
  jdates[t] = JULDAY(mm,dd,year[t],time[t]*1.e-3/3600.)
ENDFOR

; get the start and end days-of-year for the requested month
mdy_to_doy,stday,yr,mon,1
monthend = monthdays[mon]
IF (((yr MOD 4) EQ 0) AND (mon GE 2)) THEN monthend = monthend + 1
mdy_to_doy,endday,yr,mon,monthend

; select the scans that match the search criteria
scanloc_in_month = WHERE(year EQ yr AND day GE stday AND day LE endday AND lat GE min_lat AND lat LE max_lat,mct)
; get the index in the pressure array that matches the specified pressure level
presloc = WHERE(pressure EQ pres,pct)

;  plot atox values that match the search criteria (y-axis) vs.time (x-axis)
j = jdates[scanloc_in_month]
a = REFORM(atox[presloc,MIN(scanloc_in_month):MAX(scanloc_in_month)])
good_atox_loc = WHERE(a GT 0,gct)
IF gct GT 0 THEN BEGIN
  good_dates = j[good_atox_loc]
  good_atox = a[good_atox_loc] 
  p1 = PLOT(good_dates,good_atox,DIMENSIONS=[2000,500],$
  TITLE=title,MIN_VALUE=0,$
  XRANGE=[MIN(good_dates)-1,MAX(good_dates)],XTITLE=months[mon-1],$
  XTICKUNITS='Day',XTICKINTERVAL=1,XMINOR=4,$
  YTITLE='Volume Mixing Ratio (log scale)',/YLOG,YMINOR=4,COLOR='blue',THICK=2,FONT_SIZE=18)
ENDIF ELSE BEGIN
  PRINT,'No good data values found for specified selection.'
ENDELSE
stop

END

PRO doy_to_mdy,doy,year,month,day
  ;  converts day of year number for the given year
  ;  to month and day
  ;  INPUTS:
  ;    doy - day of year number
  ;    year
  ;  OUTPUTS:
  ;    month - numeric
  ;    day - day of month
  ;
  ;  LAH - 9/30/2008

  mdays = [0,31,59,90,120,151,181,212,243,273,304,334]
  leap_mdays = [0,31,60,91,121,152,182,213,244,274,305,335,366]

  IF ((year MOD 4) EQ 0) AND ((year MOD 100) NE 0) THEN mdays = leap_mdays  ;  it's a leap year

  month = MAX(WHERE(doy GT mdays,mcount)) + 1
  day = doy - mdays[month-1]

END

PRO mdy_to_doy,doy,year,month,day
  ;  converts month day year to day of year number
  ;  INPUTS:
  ;    day - day of month
  ;    month - numeric, 1-12
  ;    year
  ;  OUTPUTS:
  ;    doy - numeric day of year
  ;
  ;  LAH - 10/21/2008

  mdays = [0,31,59,90,120,151,181,212,243,273,304,334]
  leap_mdays = [0,31,60,91,121,152,182,213,244,274,305,335,366]


  IF ((year MOD 4) EQ 0) AND ((year MOD 100) NE 0) THEN mdays = leap_mdays  ;  it's a leap year

  ; get number of days in preceding months using month number,
  ; then just add days into current month

  doy = mdays[month-1] + day

END

FUNCTION lat_to_latstr,lat,fmt

  ; converts numeric latitude +/- 90 to latitude string with N or S (or neither, if equator)
  ; lat is numeric lat to be converted
  ; fmt is an IDL format string to determine the display of the string latitude, e.g., '(F5.1)'
  ; returns the latitude string
  IF lat LT 0 THEN latstr = STRTRIM(STRING(ABS(lat),FORMAT=fmt) + ' S',2) $
  ELSE IF lat GT 0 THEN latstr = STRTRIM(STRING(ABS(lat),FORMAT=fmt) + ' N',2) $
  ELSE latstr = STRTRIM(STRING(lat,FORMAT=fmt),2)

  RETURN,latstr

END