; 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           

file1 = addfile("/winraid/Proposals/Clara_Falk/QS_XWGRD3_2003120.20031231631.hdf","r")
print (file1)  

v0 = file1->des_avg_wind_speed
v1  = short2flt(v0)
v  = v1;*v0@scale_factor

;v@_FillValue=0.000

printVarSummary(v)
la = ispan(-89875,89875,250)
lat = la/1000.
lat!0 = "lat"
lat@units="degrees_north"
lat&lat = lat

lo = ispan(125,359875,250)
lon = lo/1000.
lon!0 = "lon"
lon@units="degrees_east"
lon&lon = lon

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

v&lat  = lat
v&lon =  lon 

;grab a subset of data to test
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)

  wks = gsn_open_wks("X11", "HDF QuikScat test")

  gsn_define_colormap(wks,"testcmap")      ; choose colormap
  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@cnLevelSpacingF       = 5          ; contour spacing
;  res@cnFillDrawOrder       = "PreDraw"    ; fill first

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

  res@lbLabelStride         = 6            ; skip every other label
  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@mpProjection          = "Satellite"  ; choose map projection
  res@mpCenterLonF          = 145.         ; choose center lon
  res@mpCenterLatF          =  40.         ; choose center lat
  res@mpLimitMode           = "LatLon"     ; required
  res@mpMinLatF             =  20.         ; min lat
  res@mpMaxLatF             =  60.         ; max lat
  res@mpMinLonF             =  100.        ; min lon
  res@mpMaxLonF             =  190.        ; max lon
;  res@mpOutlineBoundarySets = "NoBoundaries"
;  res@mpFillBoundarySets = "NoBoundaries"

;  res@pmTickMarkDisplayMode = "Always"     ; turn on automatic tickmarks
  res@mpGridAndLimbOn       = True         ; turn on lat/lon lines
  res@mpGridMaskMode        = "MaskLand"   ; Mask grid over land.
  res@gsnMaximize           = True         ; enlarge plot 
  res@tiMainString          = "Example QuikScat HDF plot"
  res@gsnAddCyclic          = False

      plot = gsn_csm_contour_map(wks,vp,res)



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

vp = v(sub_lat_min_pos:sub_lat_max_pos,sub_lon_min_pos:sub_lon_max_pos)
lat_sub=lat(sub_lat_min_pos:sub_lat_max_pos)
lon_sub=lon(sub_lon_min_pos:sub_lon_max_pos)

len0=dimsizes(vp)
len = len0(0)*len0(1)
print(len0(0))

table = new(len+1,"string")
table(0) ="lat,lon,val" 
curloc = 1

do lac= 0, len0(0)-1
  do loc = 0, len0(1)-1 
    cval = vp(lac,loc)
    latv = lat_sub(lac)
    lonv = lon_sub(loc)
        if (lonv.ge.180) then 
           lonv=lonv-360
	end if
      table(curloc) = sprintf("%6.3f",latv) + ","+ sprintf("%7.3f",lonv) +","+ sprintf("%6.3f",cval)
      curloc = curloc +1

    end do
end do
asciiwrite("~/nclstuff/HDFtestDump.csv",table)


end