;load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/conDIVibuted.ncl"
load "$NCARG_ROOT/lib/lingjian.ncl"

begin

	HRES = 2.5 ; horizontal resolution
	VAR = "pr"
   	fname = "PR.Filtered.2.5.nc"
   	f = addfile(fname, "r")
   	PR = f->$str_lower(VAR)$(:,{-40:250})
   	dim = dimsizes(PR)
   	; smooth in time by 5 days that will give the track a width for 5 days
  	PR = runave_n_Wrap(PR,5,0,0)
	; get the standard deviation in time for each individual longitude points
	stdPR = dim_stddev_n_Wrap(PR(:,{-40:240}),0)	
		
	NTIME = dim(0)
	delete(dim)
		
	VAR = "RAY"
    name = "MJO.Track.nc"
    print(name)
    f= addfile(name,"r")       
    tmp = f->$VAR$(:,:,{-40:250})
 	tmp = runave_n_Wrap(tmp,5,0,0)	; do the five day running mean the same as above
	TRACK = tmp(:,:,{-40:240})
	lon = tofloat(f->lon({-40:240}))
	speed = TRACK&speed
	time = TRACK&time
	delete(tmp)

;================================================
;set criterion
;================================================
	dim = dimsizes(TRACK)	
	printVarSummary(TRACK)
	NTIME = dim(0)
	NSPED = dim(1)
	NLONS = dim(2)
	
	OTRACK = TRACK ; OTRACK is used to store orignal value for all tracks
	; normalization in longitude 
	do i = 0 , NLONS - 1
		TRACK(:,:,i) = TRACK(:,:,i) / stdPR (i)
	end do		

	attr = TRACK(:,:,0:4)
	attr = attr@_FillValue
	attr!2 = "attrs"
	; 	0 	=>	amp
	;	1	=>	west boundary
	;	2	=>	east boundary
	;	3	=>	distance
	;	4	=>	total amp	
; do integration for each ray
; 1) the integration only applied to those points value > 1std
; if the value is smaller than 1 std, then set them to _FillVlaue

	Threshold = 1.
 	TRACK = where(TRACK .ge. Threshold, OTRACK, TRACK@_FillValue)
 	
;2) find the longest range for the area with points > than 1std
	do i = 90 , NTIME - 91
		do j = 0 , NSPED - 1
			if (all (ismissing(TRACK(i,j,:)))) then
				; if all the points along the track is _FillValue, then the amp is set to 0.
				attr(i,j,:) = 0
			else
				; set up a new array with two ends are _FillValue for easy handling
				tmpray = new(NLONS+2,float,TRACK@_FillValue)
				tmpray(1:NLONS) = TRACK(i,j,:)
				tmporay = tmpray
				tmporay(1:NLONS) = OTRACK(i,j,:)
				; used to store several pair of west and east boundaries
				pair = new((/100,2/),integer,-999)
				range = new(100,integer,-999)
				ix = 0
				iy = 0
			
				do k = 1 , NLONS
					; find the index of the first points that is not _FillValue
					if (.not. ismissing(tmpray(k)) .and. ismissing(tmpray(k-1))) then
						pair(ix,0) = k - 1
						ix = ix + 1  
					end if
					; find the index of the last points that is not _FillValue
					if (.not. ismissing(tmpray(k)) .and. ismissing(tmpray(k+1))) then
						pair(iy,1) = k - 1
						iy = iy + 1  
					end if
				end do
				
				range = pair(:,1) - pair(:,0)
				
				len = ix ; number of pairs
				
				; if only one pair is found
				if (len .eq. 1) then
					attr(i,j,0) = (/sum(TRACK(i,j,pair(0,0):pair(0,1)))/) * HRES
					attr(i,j,1) = (/lon(pair(0,0))/)
					attr(i,j,2) = (/lon(pair(0,1))/)
					attr(i,j,3) = (/lon(pair(0,1)) - lon(pair(0,0))/)
					attr(i,j,4) = (/sum(OTRACK(i,j,pair(0,0):pair(0,1)))/) * HRES					
				else 
					; more than one pair is found, find the longest one				
					; combined the pair that distance is smaller than 10 degree in longitude
					do k = 1, len - 1
						ik = len - k
						if (pair(ik,0) .le. pair(ik-1,1) + 4 ) then ; the gap can not more than 10 degree,
							if ( avg(tmporay(pair(ik-1,1):pair(ik,0))) .gt. 0) then ; the averaged value for the gap must larger than 0
								pair(ik-1,1) = pair(ik,1)
								pair(ik:98,:) = pair(ik+1:99,:) 
;								range(ik-1) = range(ik-1) + range(ik)
;								range(ik:98) = range(ik+1:99)
							end if
						end if	
					end do					  
					
					dist = pair(:,1) - pair(:,0)
					idx = maxind(dist)
					attr(i,j,0) = (/sum(TRACK(i,j,pair(idx,0):pair(idx,1)))/) * HRES ; store the total amplitude integrated along the track with the threshold (large > threhold)
					attr(i,j,1) = (/lon(pair(idx,0))/) 	; store the west boundary of the track
					attr(i,j,2) = (/lon(pair(idx,1))/)	; store the east boundary of the track
					attr(i,j,3) = (/lon(pair(idx,1)) - lon(pair(idx,0))/)		; store the propagation range for the track
					attr(i,j,4) = (/sum(OTRACK(i,j,pair(idx,0):pair(idx,1)))/) * HRES ; store the total amplitude integrated along the track without any threshold
				end if
			end if
		end do
		print(i + " done!!")
	end do
			
	fname = "Track.ATTR.nc"
	system("rm -rf " + fname)
	f = addfile(fname,"c")
	f->attr = attr 

end