;
;  This example code illustrates how to access and visualize GES DISC MLS
;  O3 HDF-EOS5 Swath file in NCL. 
;
;  If you have any questions, suggestions, comments on this example, please 
; use the HDF-EOS Forum (http://hdfeos.org/forums). 
;
;  If you would like to see an  example of any other NASA HDF/HDF-EOS data 
; product that is not listed in the HDF-EOS Comprehensive Examples page 
; (http://hdfeos.org/zoo), feel free to contact us at eoshelp@hdfgroup.org or 
; post it at the HDF-EOS Forum (http://hdfeos.org/forums).

; Usage:save this script and run 
; 
; $ncl MLS-Aura_L2GP-O3_v04-23-c02_2019d001.he5.ncl
;
; Tested under: NCL 6.6.2
; Last updated: 2020-11-10

begin
  file_name="MLS-Aura_L2GP-O3_v04-23-c02_2019d001.he5"
  eos_file=addfile(file_name, "r") ; Read file. 

;  print(eos_file)

; Read geo-location/time variables.
  pressure=eos_file->Pressure_O3
  pressure@units = eos_file->Pressure_O3@Units
  pressure@long_name = eos_file->Pressure_O3@Title

  lat=eos_file->Latitude_O3

; The file has "degress_north".
  lat@units = eos_file->Latitude_O3@Units
  lat@long_name = eos_file->Latitude_O3@Title

  time = eos_file->Time_O3
  time@units = "seconds since 1993-1-1 00:00:0.0"

; Create date time string.
; cd_calendar() is supported since NCL 6.0.0.
; Use ut_calendar() for older version.
  utc_date = cd_calendar(time(0), 0)
  year   = tointeger(utc_date(:,0))    
  month  = tointeger(utc_date(:,1))    
  day    = tointeger(utc_date(:,2))
  hour   = tointeger(utc_date(:,3))
  minute = tointeger(utc_date(:,4))
  second = utc_date(:,5)
  date_str = sprinti("%0.4i", year)  + "-" + sprinti("%0.2i", month) + "-" + \
             sprinti("%0.2i ", day) +  sprinti("%0.2i", hour) + ":" + \ 
             sprinti("%0.2i", minute) + ":" + sprintf("%0.2f", second)
             
  utc_date = cd_calendar(time(70), 0)
  year   = tointeger(utc_date(:,0))    
  month  = tointeger(utc_date(:,1))    
  day    = tointeger(utc_date(:,2))
  hour   = tointeger(utc_date(:,3))
  minute = tointeger(utc_date(:,4))
  second = utc_date(:,5)
  date_str2 = sprinti("%0.4i", year)  + "-" + sprinti("%0.2i", month) + "-" + \
              sprinti("%0.2i ", day) +  sprinti("%0.2i", hour) + ":" + \ 
              sprinti("%0.2i", minute) + ":" + sprintf("%0.2f", second) 

; Read the dataset.
  data2D=eos_file->L2gpValue_O3
  data = data2D
  data@long_name = data2D@Title   + " from " + date_str + "~C~ to " + date_str2
  data@units = data2D@Units
  data@_FillValue = data2D@_FillValue

; Associate data with coordinate variables.
  data&nLevels_O3 = pressure;
  data&nTimes_O3 = lat;

  xwks = gsn_open_wks ("png", file_name + ".ncl") 
  gsn_define_colormap(xwks,"BlAqGrYeOrReVi200")

  res = True ; plot mods desired
  res@cnFillOn = True ; enable contour fill
  res@gsnMaximize = True; make plot large
  res@cnLinesOn = False ; turn off contour lines
  res@cnLineLabelsOn = False; turn off contour line labels
  res@gsnSpreadColors = True ; use the entire color spectrum
  res@cnFillMode = "RasterFill" ; faster
  res@lbOrientation = "vertical" ; vertical labels
  res@cnMissingValFillPattern = 0 ; missing value pattern is set to "SolidFill"
  res@cnMissingValFillColor = 0; white color for missing values
  res@lbLabelAutoStride = True ; ensure no label overlap
  
; Put high pressure at the bottom of Y-axis.
  res@trYReverse=True
  
  res@tmYLMode   = "Explicit"
  res@tmYLValues = (/1000,100,10,1,0.1,0.0/)
  res@tmYLLabels = "" + res@tmYLValues
  gsn_define_colormap(xwks,"BlAqGrYeOrReVi200") ; define colormap
  res@tiMainString = file_name
  res@tiYAxisString = pressure@long_name + " (" + pressure@units + ")"
; Select a region that latitude values are monotoic. 
; Otherwise, NCL can't plot latitude label with a warning message:
;
;  "ScalarFieldSetValues: irregular coordinate array \
;     sfXArray non-monotonic: defaulting sfXArray"
;
; You can use HDFView to see the latitude values.
  datas = data(0:70,:)
  plot=gsn_csm_contour(xwks, datas(nLevels_O3|:,nTimes_O3|:), res) 
end

; References
;
; [1] http://hdfeos.org/software/ncl.php