; Concepts illustrated:
;   - Calculate latitude/longitude arrays that define the location 
;     of a circles around multiple central locations .
;   - Extracting variables from a variable of type 'list'
;   - Using gc_inout to mask data inside a great circle
;----------------------------------------------------------------------

;---Multiple locations for cyclone track
;---Read lat/lon data from Track file
; March
 filename = "track_mar.txt"
 ncols = numAsciiCol(filename) ; Calculate the number of columns
 data = readAsciiTable(filename,ncols,"float",0) ; Given the # of columns, we can use readAsciiTable to read this file
 slat=data(:,0)  ; Latitude  ; lat:storm centers
 slon=data(:,1)  ; Longitude  ; lon:storm centers
        
  srad      = (/ 500 /)                   ; storm radii (km)
  srad_unit = 1                           ; ; 1=km, 0=degrees
  N         = 10                          ; # of points; more points nicer 'circle'
  opt       = False

  circle    = geolocation_circle(slat, slon, srad, srad_unit, N, opt) ; circle is type list
  ;printVarSummary(circle)                  ; list containing two variables
  
  nLoc      = dimsizes(slat)               ; # of locations
  nRad      = dimsizes(srad)               ; # of circles at each location
  
  circle_lat = circle[0]                   ; For clarity: explicitly extract list elements
  circle_lon = circle[1]
  
  ;printVarSummary(circle_lat)
  ;printVarSummary(circle_lon)
  ;printMinMax(circle_lon,2)
  ;exit
  
 ;---Read in TRMM precipitation data
 ;March
 f2  = addfile("all_mar.nc", "r")
 prc2 =f2->pre(:,{75:100},{0:30})
 prc2_reorder = prc2(time|:,latitude|:,longitude|:)
 value1 = 1d20
 prc2_reorder@missing_value = value1                  ; Force change to CF recognized value
 prc2_reorder@_FillValue    = value1
 
; March
 ntim2 = (/26,27,28,29/) 
 tcmar = prc2_reorder(ntim2,:,:)
 tcavgmar = dim_avg_n_Wrap(tcmar,0)
 
 PRC = where(tcavgmar.lt.1,tcavgmar@_FillValue,tcavgmar)
 copy_VarMeta(tcavgmar,PRC)
 
 ;printVarSummary(PRC)
 ;printMinMax(PRC,0)
 
  latmin = 0
  latmax = 30
  lonmin = 75
  lonmax = 100  
   
;---Start the graphics
   plot      = new ( nLoc, "graphic") 
   wks_type = "png"
   wks_type@wkWidth = 2500
   wks_type@wkHeight = 2500 
   wks = gsn_open_wks(wks_type,"tc500km_mar3")

  res                       = True
  res@gsnMaximize           = True
  
  res@gsnAddCyclic          = False
  res@mpLimitMode = "LatLon"        ; select subregion 
  res@mpMinLonF             = lonmin
  res@mpMaxLonF             = lonmax
  res@mpMinLatF             = latmin
  res@mpMaxLatF             = latmax
  ;res@mpCenterLonF          = (lonmin+lonmax)/2
  ;res@mpProjection = "Mercator"
  res@mpDataBaseVersion    = "HighRes"                ; use GMT coastline
  res@pmTickMarkDisplayMode = "Always"     ; nicer tickmarks for regional plots
 
  res@cnFillOn              = True
  res@cnLinesOn             = False
  res@cnInfoLabelOn        = False              ; turn off contour info label
  res@cnLinesOn            = False              ; turn off contour lines
  res@cnLineLabelsOn       = False              ; turn off line labels
  res@mpFillOn             = False              ; turn off gray continents
  res@lbLabelBarOn         = True              ; add single label bar
  res@lbOrientation = "horizontal"
  ;res@cnFillMode            = "RasterFill"
  res@cnFillPalette         = "precip3_16lev"
  ;res@trGridType           = "TriangularMesh"
  res@gsnMaximize         = True   
  res@gsnPaperOrientation = "auto"
  res@gsnLeftString            = " "
  res@gsnRightString            = " "
  res@mpShapeMode  = "FreeAspect"
  res@vpWidthF      = 0.9
  res@vpHeightF     = 0.7 
 
  res@tmYROn      = False     ; turn off right and top tickmarks
  res@tmXTOn      = False	
  
    
  ;; Cyclone Tracking
;March
 filename4 = "00_bob01.txt"
 ncols = numAsciiCol(filename4) 
 data4= readAsciiTable(filename4,ncols,"float",0) 
 y4=data4(:,0)  
 x4=data4(:,1) 
 
 ; polyline resources
 resp = True                                ; mods yes
 resp@gsLineThicknessF = 12.0                ; line thickness
 resp@gsLineColor      = "black" 
 
 ;Add markers to the trajectories.
   amres                = True         ; marker resources for best track
   amres@gsMarkerIndex  = 16           ; marker style (filled circle)
   amres@gsMarkerSizeF  = 6.0          ; marker size
   amres@gsMarkerColor  = "red"      ; maker color  
   
 ;creat panel resources
resP             = True
resP@gsnMaximize = True
resP@gsnPanelLabelBar =True  ; add common label bar
;resP@gsnPanelMainString = "Yearly TCs mean rainfall"      
   
 clat = new((/15,15,10/),float)   ; arrays to hold great circle
 clon= new((/15,15,10/),float)
 
 do i =0,nRad-1
    do j = 0,nLoc-1
       clat(j,i,:) = circle_lat(j,i,:)
       clon(j,i,:) = circle_lon(j,i,:)
	
    clat1 = ndtooned(clat)
    clon1 = ndtooned(clon)

    ;printVarSummary(clat)
    ;printVarSummary(clon)
    ;printVarSummary(clon1)

    min_lat = min(clat1)
    min_lon = min(clon1)
    max_lat = max(clat1)
    max_lon = max(clon1)
   ;printMinMax(min_lat,0)
   ;printMinMax(min_lon,0)
   ;printMinMax(max_lat,0)
   ;printMinMax(max_lon,0)

   ;printMinMax(clon,2)
   ;printMinMax(clat1,0)
   ;printMinMax(clon1,0)
 
 
;---Subset the desired rectagle of data 
    R2 := PRC({min_lat:max_lat},{min_lon:max_lon})
    
;---Set points that are outside of the circle of data to missing
    lat2d := conform(R2,R2&latitude,0)
    lon2d := conform(R2,R2&longitude,1)
    in_circle := gc_inout(lat2d,lon2d,clat1,clon1)
    R2  = where(in_circle,R2,R2@_FillValue)
	
    ;printVarSummary(R2)
    ;printMinMax(R2,1)
   
   
  ;---Create a plot 
    res@gsnLeftString = "e) Mean TCs rainfall (mm): March"
    res@cnLevelSelectionMode = "ExplicitLevels"   ; set explicit contour levels
    res@cnLevels             = (/1,5,10,20,30,40,50,60,70,80/) 
   
	if (j.eq.0) then
       plot(j) = gsn_csm_contour_map(wks,R2,res)
    else
	   plot(j) = gsn_csm_contour_map(wks,R2,res)
       overlay(plot(0),plot(j))
    end if
   end do	
	
    ;d1 = gsn_add_polyline(wks,plot(0),x4,y4,resp)
    ;d2 = gsn_add_polymarker(wks,plot(0),x4(0),y4(0),amres)
   
    gsn_panel(wks,plot(0),(/1,1/),resP)      ; alternatively can just call draw(plot(0)) followed by frame(wks)
end do