;  This example code illustrates how to access and visualize CDPC CloudSat 
; 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).

; Tested under: NCL 6.0.0
; Last updated: 2011-11-16

load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"


  latS     = -20.0
  latN     =  20.0

  file_name = "2010128055614_21420_CS_2B-GEOPROF_GRANULE_P_R04_E03.hdf"

  eos_file=addfile(file_name+ ".he2","r")    ; .he2 handle as a HDF4-EOS2 file 
  print(eos_file)                            ; overview of EOS file as NCL 'sees' it

  data_raw = eos_file->Radar_Reflectivity_2B_GEOPROF ; (nray_2B_GEOPROF, nbin_2B_GEOPROF)
                                                     ; (37082, 125)

  lat = eos_file->Latitude_2B_GEOPROF        ; (nray_2B_GEOPROF) ==> (37082)
  lon = eos_file->Longitude_2B_GEOPROF       ; (nray_2B_GEOPROF)
  lev = eos_file->Height_2B_GEOPROF          ; (nray_2B_GEOPROF, nbin_2B_GEOPROF)

  lev@long_name  = lev@hdfeos_name           ; long_name is COARDS & CF conformant 
  lev@units      = eos_file@Height_units_2B_GEOPROF

  time           = eos_file->Profile_time_2B_GEOPROF          ; (nray_2B_GEOPROF)
  time@long_name = eos_file@Profile_time_long_name_2B_GEOPROF
  time@units     = eos_file@Profile_time_units_2B_GEOPROF

  hdf4_file=addfile(file_name,"r")           ; open as standard hdf file 
  print(hdf4_file)                           ; file overview 

  data_hdf4 = hdf4_file->Radar_Reflectivity
  printVarSummary(data_hdf4)                 ; variable overview
  print("data_hdf4: min="+min(data_hdf4)+"  max="+max(data_hdf4))

; Process valid_range. Fill value and missing value will be handled by this
; since they are outside of range values.
  data  = where(data_raw.gt.data_hdf4@valid_range(0) .and. data_raw.lt.data_hdf4@valid_range(1) \
                ,data_raw, data_hdf4@_FillValue)

; Apply scale factor according to the data spec [1].
  data@_FillValue = data_hdf4@_FillValue
  dataf = tofloat(data)
  dataf = dataf / data_hdf4@factor

  dataf!0 = data_raw!0                        ; assign original named dimensions
  dataf!1 = data_raw!1
  printVarSummary(dataf)

; Although 2D height values are all slightly different at each profile 
; time, the difference is not significant. Pick the first one only since
; NCL doesn't allow 2D array for axis variable.
;;dataf&nbin_2B_GEOPROF = lev(0,:)
;;dataf&nray_2B_GEOPROF = time

; The desired contour plot wants the x-axis to be time. Reorder
; so time is the right dimension. The plot template expects this.

  data_plt = dataf(nbin_2B_GEOPROF|:,nray_2B_GEOPROF|:)
  lat2d          = conform(data_plt, lat, 1)
  data_plt       = where(lat2d.lt.latS .or. lat2d.gt.latN \
                        ,data_plt@_FillValue, data_plt)

                 ; arbitrary filter to reduce plot clutter
  data_plt       = where(data_plt.lt.-20, data_plt@_FillValue, data_plt) 

  if (latS.lt.0) then
      latLabelS = abs(latS)+"S"
  else
      latLabelS = latS+"N"
  end if

  if (latN.lt.0) then
      latLabelN = abs(latN)+"S"
  else
      latLabelN = latN+"N"
  end if

  latLabel = latLabelS+"-"+latLabelN

  pltTitle = str_get_cols(file_name,0,18)+"_"+latLabel
;;xwks = gsn_open_wks("pdf", file_name); open workstation
  xwks = gsn_open_wks("pdf", pltTitle) ; open workstation

  gsn_define_colormap(xwks,"BlAqGrYeOrReVi200") ; define colormap

  res = True                           ; plot mods desired
  res@cnFillOn = True                  ; enable contour fill
  res@gsnMaximize = True               ; make plot large
  res@gsnPaperOrientation = "portrait" ; force portrait orientation
  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
;;res@tiMainString = file_name
  res@tiMainString = str_get_cols(file_name,0,18)+": "+latLabel
  res@tiXAxisString = time@long_name+" ("+time@units+")"
  res@tiYAxisString = lev@long_name+" ("+lev@units+")" 
  res@gsnLeftString=data_hdf4@long_name 
  res@gsnRightString=data_hdf4@units 
  
  res@sfXArray   = conform(data_plt,time, 1)
  res@sfYArray   = lev(nbin_2B_GEOPROF|:,nray_2B_GEOPROF|:)
  res@trGridType = "Triangularmesh"


  plot=gsn_csm_contour(xwks, data_plt, res)

;----------------------------------------------------------------
; References
;
; [1] http://www.cloudsat.cira.colostate.edu/dataSpecs.php