load "./grid2geocircle.ncl"
;================================================================================================
;                               MAIN
;================================================================================================

;---Open & Read a rectilinear grid

   diri  = "/Users/shea/Data/GRIB/"
   fili  = systemfunc("cd "+diri+" ; ls ockhi03b*.grb2")
   pthi  =  diri+fili
   setfileoption("grb","SingleElementDimensions","Forecast_time")
   vName = "VVEL_P0_L100_GLL0"

   f     =  addfiles(pthi,"r")
   x     = f[:]->$vName$           ; [forecast_time0 | 5] x [lv_ISBL0 | 46] x [lat_0 | 721] x [lon_0 | 881] 
   printVarSummary(x)

   dimx  = dimsizes(x)
   rankx = dimsizes(dimx)
   ntim  = dimx(0)

   if (rankx.eq.4) then
       nlev = dimx(1)
       x = x(:,:,::-1,:)         ; linint2_points [used by rgrid2geolocation] requires ascending latitude order]
   else
       x = x(:,::-1,:)
   end if

;---Specify the central lat/lon point(s). 

   clat  = 10.0   ;  could be more than on: (/  0.0, 10.0 /) 
   clon  = 70.5   ;                         (/ 45.0, 70.5 /) 
   clat!0= "center_lat"
   clon!0= "center_lon"
   nctr  = dimsizes(clat)

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

   nrad         = 17
   radius       = fspan(0,4,nrad)        ; radii (degrees)
   radius!0     = "radius"
   radius@units = "degrees"
   radius@long_name = "Radii"
  ;nrad         = dimsizes(radius)

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

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

   xrad = x_geocircle[0]
 ;;plat = x_geocircle[1]   ; not needed      
 ;;plon = x_geocircle[2]
   printVarSummary(xrad)      ; [forecast_time0 | 5] x [center | 2] x [lv_ISBL0 | 46] x [radius | 33] x [circle | 90]

;---Average all 'geo-circles' for each time, each center, each pressure level and each radius

   xrad_avg = dim_avg_n_Wrap(xrad, rankx)  
   printVarSummary(xrad_avg)  ; [forecast_time0 | 5] x [center | 2] x [lv_ISBL0 | 46] x [radius | 33]

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

   ymdh = cd_calendar(x&forecast_time0,-3)  ; used in plot totle 
   print(ymdh)
   print("-----")

   pltType = "png"
   pltName = "rgrid2geocircle_"+vName
   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= "BlAqGrYeOrReVi200"
 
   extra = 0.25                     ; If you want to have a small areal buffer
   res@mpMinLatF    = min(x&lat_0)-extra
   res@mpMaxLatF    = max(x&lat_0)+extra
   res@mpMinLonF    = min(x&lon_0)-extra
   res@mpMaxLonF    = max(x&lon_0)+extra
 
 ;---Draw map of the original variable  at specified level (plev)
 
   nt_plot = ntim-1    ; arbitrary time index; (ntim-1) is the last time
   plev    = 50000     ; Pa
   res@tiMainString = ymdh(nt_plot)+": plev="+plev 

 ;;optsd = True
 ;;optsd@PrintStat = True
 ;;statsd = stat_dispersion(x(nt_plot,{plev},:,:), optsd )   ; explore data being plotted
 
   if (vName.eq."VVEL_P0_L100_GLL0") then            ; want symmetric 
      ;res@cnFillPalette = "ViBlGrWhYeOrRe"
       res@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
       res@cnMinLevelValF       = -2.00              ; set min contour level
       res@cnMaxLevelValF       =  2.00              ; set max contour level
       res@cnLevelSpacingF      =  0.05              ; set contour spacing
       res@cnFillMode           = "RasterFill"
   end if

   plot = gsn_csm_contour_map(wks,x(nt_plot,{plev},:,:), res)
 
   plres  = True
   plres@gsMarkerIndex = 2           ; Plus sign ( + )
  ;plres@gsMarkerIndex = 15          ; Circle with x
  ;plres@gsMarkerIndex = 16          ; Filled Circle
   plres@gsMarkerSizeF = 0.015     
   plres@gsMarkerThicknessF = 0.025     
   plres@gsMarkerColor = "foreground"; "white"  
   gsn_polymarker(wks, plot, clon, clat, plres) 
   frame(wks)

   delete( [/ res@mpMinLatF, res@mpMaxLatF, res@mpMinLonF, res@mpMaxLonF /] )  ; not appropriate for X-section
   if (isatt(res,"cnFillOn")) then
       delete( res@cnFillMode )              
   end if

;---Draw cross-section
 
   res@gsnFrame            = True
   res@tiMainFontHeightF   = 0.020     ; default 0.025
   res@lbOrientation       = "vertical"
 
   res@trYReverse          = True
  ;res@gsnYAxisIrregular2Log = True    ; Convert Y axis to logarithmic
   res@tiYAxisString       = "Pressure (Pa)"                   
   res@tiXAxisString       = "Radius (degrees)"  
 
   if (vName.eq."VVEL_P0_L100_GLL0") then       ; want symmetric 
       res@cnFillPalette = "ViBlGrWhYeOrRe"

     ;;optsd = True
     ;;optsd@PrintStat = True
     ;;stat_xrad_avg = stat_dispersion(xrad_avg(nt_plot,:,:,:), optsd )
 
       symMinMaxPlt (xrad_avg(nt_plot,:,:,:),20,False,res)    
       res@cnMinLevelValF       = -5.00              ; set min contour level
       res@cnMaxLevelValF       =  5.00              ; set max contour level
       res@cnLevelSpacingF      =  0.25              ; set contour spacing
   end if
 
;---For each time, center ... draw a contour
 
   do nt=nt_plot,nt_plot  ; nt=0,ntim-1  ==> all
     do nc=nctr-1,nctr-1  ; nc=0,nctr-1  ==> all
        res@tiMainString = "rgrid2geocircle: "+ymdh(nt)+": ["+sprintf("%5.1f",clat(nc)) \
                                                       +","+sprintf("%5.1f",clon(nc))+"]"
        plot = gsn_csm_contour(wks,xrad_avg(nt,nc,:,:),res)       ; contour the variable
     end do
   end do