;----------------------------------------------------------------------
; This script plots both ascending and descending wind speed over 
; a single map plot. Note that these fields span a large area of the
; globe, so if you zoom in on India, you will not see the full fields.
;----------------------------------------------------------------------
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_overlay")
  cmap = read_colormap_file("BlAqGrYeOrRe")

;---Set resources specific to contour plots
  cnres                 = True
  cnres@gsnDraw         = False    ; we will draw and frame later
  cnres@gsnFrame        = False    ; after doing an overlay

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

  cnres@cnLevelSelectionMode = "ExplicitLevels"
  cnres@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/)

  cnres@cnInfoLabelOn   = False
  cnres@cnLinesOn       = False
  cnres@cnLineLabelsOn  = False

;---Set resources specific to contour/map plot
  mpres                       = cnres    ; copy all the contour resources first
  mpres@mpFillOn              = False    ; zoom in on map area of interest

;---We can't zoom in because the two wind speeds are in different parts of the globe.
;  mpres@mpMinLatF             = -60
;  mpres@mpMaxLatF             =  60
;  mpres@mpMinLonF             =  30
;  mpres@mpMaxLonF             = 120

  mpres@tiMainString  = "Ascending  / descending wind speed"
  plot_asc = gsn_csm_contour_map(xwks,asc_wind_speed,mpres)

;
; For descending wind speed, create only a contour plot.
; We can then overlay this on the contour/map plot
; we just created for ascending wind speed.
;
  cnres@lbLabelBarOn  = False
  plot_des = gsn_csm_contour(xwks,des_wind_speed,cnres)

  overlay(plot_asc,plot_des)

  draw(plot_asc)    ; This will draw both the ascending and descending wind speed.
  frame(xwks)

end