; HDF4 format reader for NASA QuikScat daily marine winds data sets
; - also dumps a selected subregion to a CSV file for import
;   into other programs
;
; David E. Atkinson, IARC/University of Alaska Fairbanks
; January 2006
;
; load the libraries just in case 

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

begin           
  fil   =  "QS_XWGRD3_2003120.20031231631.hdf"
  latS  =  20.         ; min lat
  latN  =  60.         ; max lat
  lonL  =  115.        ; min lon
  lonR  =  175.        ; max lon

  filc  = stringtochar(fil)     
 ;print(filc)
  yr1   = stringtointeger((/filc(10:13)/)) 
  

  file1   = addfile(fil,"r")
  print (file1)  

  v  = short2flt( file1->des_avg_wind_speed )
  v@_FillValue = 0.000                          ; original missing value

  v@_FillValue = 1.e20                          ; change 1.e20

  nlat= 720
  mlon= 1440
  lat = latGlobeFo (nlat, "lat", "latitude", "degrees_north")
  lon = lonGlobeFo (mlon, "lon", "longitude", "degrees_east")

  v!0 = "lat"
  v!1 = "lon"

  v&lat  = lat
  v&lon =  lon 

  printVarSummary(v)
  printMinMax(v,True)

;grab a subset of data to test
;Dimensions and sizes:   [lat | 151] x [lon | 241]
;Coordinates: 
;           lat: [22.625..60.125]
;           lon: [115.125..175.125]
;sub_lat_min_pos = 450
;sub_lat_max_pos = 600
;sub_lon_min_pos = 460
;sub_lon_max_pos = 700
;vp = v(sub_lat_min_pos:sub_lat_max_pos,sub_lon_min_pos:sub_lon_max_pos)
;printVarSummary(vp)
;printMinMax(vp,True)

  jlat = ind(lat.ge.latS .and. lat.le.latN)
  ilon = ind(lon.ge.lonL .and. lon.le.lonR)
  printVarSummary(v(jlat,ilon))
  printMinMax(v(jlat,ilon),True)
 ;print(v(jlat,ilon))                   ; activate if you want to see values



  wks = gsn_open_wks("X11", "HDF QuikScat test")
 ;gsn_define_colormap(wks,"testcmap")      ; choose colormap
  gsn_define_colormap(wks,"BlGrYeOrReVi200")      ; choose colormap [202 colors]
  i   = NhlNewColor(wks,0.8,0.8,0.8)      ; add gray to map

  res                       = True         ; plot mods desired
  res@cnFillOn              = True         ; turn on color
  res@cnLinesOn             = False        ; turn off contour lines

  res@cnFillMode            = "RasterFill"       ; Raster Mode

 ;res@cnLevelSelectionMode  = "ManualLevels"
 ;res@cnMinLevelValF        = 0 
 ;res@cnMaxLevelValF        = 50 
 ;res@cnLevelSpacingF       = 10 
 ;res@cnLevelSpacingF       = 2          ; contour spacing

  res@lbLabelAutoStride     = True
  res@lbBoxLinesOn          = False        ; turn off box between colors

  res@gsnSpreadColors       = True         ; use full colormap
  res@gsnSpreadColorStart   = 10           ; start at color 10
  res@gsnSpreadColorEnd     = 185          ; end at color 185
  res@gsnAddCyclic          = False

  res@mpProjection          = "Satellite"  ; choose map projection
  res@mpCenterLonF          = 145.         ; choose center lon
  res@mpCenterLatF          =  40.         ; choose center lat
  res@mpLimitMode           = "LatLon"     ; required
  res@mpMinLatF             =  latS         ; min lat
  res@mpMaxLatF             =  latN        ; max lat
  res@mpMinLonF             =  lonL        ; min lon
  res@mpMaxLonF             =  lonR        ; max lon
  res@mpGridAndLimbOn       = True         ; turn on lat/lon lines
 ;res@mpGridMaskMode        = "MaskLand"   ; Mask grid over land.
  res@mpFillDrawOrder       = "PostDraw"
  res@mpPerimDrawOrder      = "PostDraw"
  res@gsnMaximize           = True         ; enlarge plot if ps or eps

  res@tiMainString          = fil

  plot = gsn_csm_contour_map(wks,v(jlat,ilon) ,res)         ; plot subset

; Dump the extracted subset to a CSV file for import into 
; a GIS system or something. Use "grid2triple"

  gisData = grid2triple(lon(ilon),lat(jlat),v(jlat,ilon))   ; make x, y, value
 ;print(gisData(0,:)+"  "+gisData(1,:)+"  "+gisData(2,:) )

  gisData!0   = "row"    ; name dimensions so we can reorder for GIS_DATA
  gisData!1   = "col" 

  opt         = True
  opt@fout    = "GIS_DATA_"+filc(0:28)
  system("/bin/rm -f "+opt@fout)
  write_matrix (gisData(col|:,row|:), "3f10.3", opt)
end