;
;  This Script makes spatial plot of nearSurface Rain observed from radar 
;;;;   The Dimesions of the above files is 9246X 49                           
;
;
    load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
    load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin

;  1. Read file and Print Meta Info

 FILES       = systemfunc("ls ./DATA_FNL/2A25*.*HDF")
                   print(FILES)
                   print(FILES(0))
 F           = addfiles(FILES,"r")
                   printVarSummary(F)
;               print(F)
;               print(F[1])
 nf          = dimsizes(FILES)
                   print(nf)
 Rows        = nf/2
 Cols        = nf/2

;
  Wks_Rain      =   gsn_open_wks("png","NearSurfRainfnl")
                    gsn_define_colormap(Wks_Rain,"default")

  PLOT          = new(nf,graphic)
  IND_ID        = new(nf,graphic)

;;;;;;;   SHAPE FILES ADDITION ;;;;;;;;;;;;;;;;;;;;

    IND_Shape_Fname         =  "/home/aditya/geeta/ncl/TS_rain-U/IND_adm1.shp"
    Shpres                  =  True
    Shpres@gsLineColor      = "Blue" 
    Shpres@gsLineThicknessF = 0.5          
 

     res                      =   True
     res@mpMinLatF            =   5.00
     res@mpMaxLatF            =   20.0
     res@mpMinLonF            =   72.0
     res@mpMaxLonF            =   85.0
     res@mpFillOn             =   False

     res@gsnDraw              =   False
     res@gsnFrame             =   False
     res@gsnMaximize          =   True

     res@cnFillOn             =   True 
     res@cnLevelSelectionMode =   "ExplicitLevels"
     res@cnLevels             =   (/1,10,30/)            
;     res@cnLevels             =   ispan(2,20,4)    
     NLevels                  =   dimsizes(res@cnLevels)
;;  
     res@tmXTOn               =   False
     res@tmYROn               =   False

     res@lbLabelBarOn         =   False
;
     do   ii     =   0, nf-1
     For_Date    =   str_split(FILES(ii),".")
                     print(For_Date)
;     print(F[ii])
     longitude   =   F[ii]->Longitude
     latitude    =   F[ii]->Latitude
     Surf_Rain   =   F[ii]->nearSurfRain
     Surf_Rain@_FillValue   =   -99                     
     Surf_Rain   = where(Surf_Rain.le.0,Surf_Rain@_FillValue,Surf_Rain)

;          print(";;;; of Values< "+ res@cnLevels(0) + " = " +Surf_Rain.lt.res@cnLevels(0))

     Surf_Rain!0         =  "longitude"
     Surf_Rain!1         =  "latitude"
     Surf_Rain@lon2d     =  longitude              
     Surf_Rain@lat2d     =  latitude                       
     longitude@units     = "degrees_east"
     latitude@units      = "degrees_north"
;                            printVarSummary(Surf_Rain)
     res@gsnCenterString = For_Date(1)+"."+For_Date(2)
     res@gsnLeftString   = "     "                           
     res@gsnRightString  = "     "                           
     PLOT(ii)            =  gsn_csm_contour_map(Wks_Rain,Surf_Rain,res) ; keep Long in X axis
     IND_ID(ii)          =  gsn_add_shapefile_polylines(Wks_Rain,PLOT(ii),IND_Shape_Fname,Shpres) 

     delete(Surf_Rain)
     delete(latitude)
     delete(longitude)
     delete(For_Date)
     end do
;  
;  Maximise the Plots in the Frame
    Maximise_Plot          =  True
                             maximize_output(Wks_Rain,Maximise_Plot)
     
;   PANEL PLOTS
    Panel_Res                     = True
    Panel_Res@gsnPanelLabelBar    = True     
    Panel_Res@txString            = "Near Surface Rain (mm/hr)"
    Panel_Res@gsnPanelWhiteSpacePercent = 0.0 

                   gsn_panel(Wks_Rain,PLOT,(/Rows+1,Cols+1/),Panel_Res) ; keep Long in X axis
end