; Concepts illustrated:
;   - Specifying the central lat and lon
;   - Specifying the radial distances of points to be interpolated
;   - Use local function 'rgrid2geocircle' located in 'grid2geolocation'
;     to get the radial locations via 'nggcog'; use bilinear
;     interpolation (linint2_points) to get the values at various radii.
;   - Plotting the results     
;     Contour radial values         

;-------------------------------------------------
load "./grid2geocircle.ncl"        ; radial interpolation library
;================================================================================================
;                               MAIN
;================================================================================================

;---Read in TRMM precipitation data

 f   = addfile("aila.nc", "r")
 prc = f->pcp(:,{0:30},{75:100})    ; prc[time | 32] x [latitude | 120] x [longitude | 100]
 
 ;printVarSummary(prc)
   
  dimprc  = dimsizes(prc)
  rankprc = dimsizes(dimprc) 
  ntim  = dimprc(0)
   
;---Specify the central lat/lon point(s). 

   clat  = 16.5   
   clon  = 88.0   
   clat!0= "center_lat"
   clon!0= "center_lon"
   nctr  = dimsizes(clat)

;---Specify radii: can be unequally spaced but must be monotonically increasing

   nrad             = 40
   radius           = fspan(0,300,nrad)  ; radii (kilometers)
   radius!0         = "radius"
   radius@units     = "kilometers"
   radius@long_name = "Radii"
   nrad             = dimsizes(radius) ; 40

;---Calculate the interpolated value and lat/lon locations of each circle at each radius
    
   N = 180
   prc_geocircle = rgrid2geocircle(prc, clat, clon, radius, N, False)
   ;printVarSummary(prc_geocircle)
   

;---Clarity: Explicitly extract the returned variables from the returned list variable

   prcrad = prc_geocircle[0]
   plat   = prc_geocircle[1]          ; radial [azimuthal] latitudes 
   plon   = prc_geocircle[2]          ;   "        :       longutudes
  ;printVarSummary(prcrad)            ; [time | 32] x [center | 1] x [radius | 17] x [circle | 180]
  ;printVarSummary(plat)              ; [center | 1] x [radius | 17] x [circle | 180]          


  time_avg = dim_avg_n_Wrap(prcrad,0)
  ;printVarSummary(time_avg)
  

;================================================================================================
;  PLOT
;================================================================================================

   pltType = "png"
   pltName = "center_1"
   pltDir  = "./"
   pltPath = pltDir+pltName

   wks  = gsn_open_wks(pltType,pltPath)


   res                 = True          ; plot mods desired
   res@gsnMaximize     = True          ; maximize plot size
   res@gsnFrame        = False
   res@gsnAddCyclic    = False
 
   res@cnFillOn        = True          ; turn on color
   res@cnLinesOn       = False         ; turn off contour lines
   res@cnLineLabelsOn  = False      ; turn off contour line labels
   res@cnFillPalette   = "matlab_jet"
 
   extra = 0.25                                ; If you want to have a small areal buffer
   res@mpMinLatF    = min(prc&latitude)-extra
   res@mpMaxLatF    = max(prc&latitude)+extra
   res@mpMinLonF    = min(prc&longitude)-extra
   res@mpMaxLonF    = max(prc&longitude)+extra
                                            
; add marker for central location

   plres                    = True
   plres@gsMarkerIndex      = 16          ; Filled Circle
   plres@gsMarkerSizeF      = 0.015     
   plres@gsMarkerThicknessF = 0.10      
   plres@gsMarkerColor      = "black"    
   
   
; Azimuthal style plot

   nc_plot = 0                             ; currently only 1 center [index 0]
   radData = time_avg(nc_plot,:,:)
   ;printVarSummary(radData)                ; [radius | 17] x [circle | 180]


   res@sfXArray  = plon(nc_plot,:,:)       ;---Set coordinates          
   res@sfYArray  = plat(nc_plot,:,:)      
   ;printVarSummary(plon(nc_plot,:,:))      ; radius | 17] x [circle | 180]      
   ;printMinMax(res@sfXArray,0)
   ;printMinMax(res@sfYArray,0)

   extra = 0.10                            ; If you want to have a small areal buffer
   res@mpMinLatF = min(plat)-extra 
   res@mpMaxLatF = max(plat)+extra
   res@mpMinLonF = min(plon)-extra
   res@mpMaxLonF = max(plon)+extra
    
 
   plot = gsn_csm_contour_map(wks,radData, res)
   gsn_polymarker(wks, plot, clon, clat, plres) 
   frame(wks)
                                                                               ; X-section is not a map
   delete( [/ res@mpMinLatF, res@mpMaxLatF, res@mpMinLonF, res@mpMaxLonF /] )  ; not appropriate for X-section
   delete( [/ res@sfXArray, res@sfYArray /] )                                  ; not appropriate for X-section