;----------------------------------------------------------------------
; This script plots both ascending and descending wind speed on
; separate map plots. This is so you can see that they are on
; different parts of the globe.
;----------------------------------------------------------------------
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"

;----------------------------------------------------------------------
; Function to read a variable off a file and unscale it.
;----------------------------------------------------------------------
function unscale_data(a,varname)
local data,data_unscaled
begin
  data = a->$varname$
  printVarSummary(data)
  printMinMax(data,0)
  data_unscaled=tofloat(data)
  data_actual=data_unscaled*0.01
  copy_VarMeta(data,data_actual)
  return(data_actual)
end

;----------------------------------------------------------------------
; Procedure to attached lat/lon arrays as coordinate arrays to an
; existing variale.
;----------------------------------------------------------------------
procedure add_latlon(data,lat,lon)
begin
  data!0 = "lat"
  data!1 = "lon"
  data&lat = lat
  data&lon = lon
end


;----------------------------------------------------------------------
; Main code
;----------------------------------------------------------------------
begin
;---Open HDF5 file
  dir = "/home/user/jayanti/scatsat-1/SWATH DATA/2017/Jan/"
  dir = "./"
  f   = addfile(dir + "S1L3WW2017001_25km.h5","r")

;---Create lat/lon coordinate arrays
  dims   = getfilevardimsizes(f,"Ascending_wind_speed")
  numlat = dims(0)  ; 720
  numlon = dims(1)  ;  1440

  lon = (360./numlon)*(ispan(0,numlon-1,1)+0.5)
  lat = (180./numlat)*(ispan(0,numlat-1,1)+0.5)-90
  lat!0   = "lat"
  lon!0   = "lon"
  lat@units="degrees_north"
  lon@units="degrees_east"

;---Read and unscale ascending/descending wind speed
  asc_wind_speed = unscale_data(f,"Ascending_wind_speed")
  des_wind_speed = unscale_data(f,"Descending_wind_speed")

;---Set 0 values to missing
  asc_wind_speed@_FillValue = -999.
  des_wind_speed@_FillValue = -999.
  asc_wind_speed = where(asc_wind_speed.eq.0,asc_wind_speed@_FillValue,asc_wind_speed)
  des_wind_speed = where(des_wind_speed.eq.0,des_wind_speed@_FillValue,des_wind_speed)

;---Attach lat/lon coordinate arrays
  add_latlon(asc_wind_speed,lat,lon)
  add_latlon(des_wind_speed,lat,lon)

;---LOOK AT YOUR DATA!
  printVarSummary(asc_wind_speed)
  printVarSummary(des_wind_speed)
  printMinMax(asc_wind_speed,0)
  printMinMax(des_wind_speed,0)

;----------------------------------------------------------------------
; Start the graphics
;----------------------------------------------------------------------
  xwks = gsn_open_wks("png","asc_des_wspd_maps")
  cmap = read_colormap_file("BlAqGrYeOrRe")

;---Set resources for contour/map
  res              = True
  res@gsnMaximize     = True

  res@cnFillOn        = True
  res@cnFillMode      = "RasterFill"  ; can be faster
  res@cnFillPalette   = cmap(32:,:)  ; start at green, colors 0-31 are bluish

  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels             = (/1,1.5,2,2.25,2.5,2.75,3,3.5,4,4.5,5.2,5.6,6,6.4,6.8,7.2,7.6,8.5,\
                               11.8,12.4,14.0,14.5,15,30/)


  res@cnInfoLabelOn   = False
  res@cnLinesOn       = False
  res@cnLineLabelsOn  = False

  res@mpFillOn              = False    ; zoom in on map area of interest
; res@mpMinLatF             = -60
; res@mpMaxLatF             =  60
; res@mpMinLonF             =  30
; res@mpMaxLonF             = 120

  res@tiMainString = "Ascending wind speed"
  plot_asc = gsn_csm_contour_map(xwks,asc_wind_speed,res)

  res@tiMainString = "Descending wind speed"
  plot_des = gsn_csm_contour_map(xwks,des_wind_speed,res)
end