;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"
;This program is used to calculate the MJO tracking

load "$NCARG_ROOT/lib/lingjian.ncl"

begin
;************************************************
;define the parameter in this script
;************************************************
    FPATH = "./"
;******************************************************
; define the MaDIVix
;******************************************************
	VAR = "pr"
    name = FPATH+ VAR + ".Filtered.2.5.nc"
    print(name)
    f= addfile(name,"r")       
    PR = f->$str_lower(VAR)$
    dim = dimsizes(PR)
   	NYEAR = dim(0) / 365
	NTIME = dim(0)

    lon = f->lon
         
    time = f->time
    NLONS = dimsizes(lon)
	
	; define the eastward propagation speed, from 2 - 16 m/s with interval of 1 m/s
	speed = ispan(20,160,10) / 10.
	; calculate the slope of the line
	K = 20000*1000/speed/24/3600 /180
	
	;fixed the starting longitude to 30
			
	NK = dimsizes(speed)		
	RAY = new((/NTIME,NK,NLONS/),float, -999)
	RAY!0 = "time"
	RAY!1 = "speed"
	RAY!2 = "lon"
	RAY&time = ispan(1,NTIME,1)
	RAY&speed = speed
	RAY&lon = lon
		
	do ik = 0 , NK - 1
		dt  = K(ik) * doubletofloat(lon - lon({90}))
		L1 = ceil(dt) - dt
		L2 = 1 - L1
		
		do i = 0 , NLONS - 1
			xt = floattointeger(ceil(dt(i)))
			PR1 = PR({90+xt:NTIME-120-1+xt},i)
			PR2 = PR({90+xt-1:NTIME-120-1+xt-1},i)
			
			XPR = PR1 * L2(i) + PR2 * L1(i)
			RAY({90:NTIME-120-1},ik,i) = (/XPR/)	
			delete(PR1)
			delete(PR2)
			delete(XPR)		
		end do
		print(speed(ik) + " done!")
	end do
	
	
	fname = "MJO.Track.nc"
	system("rm -rf " + fname)
	f = addfile(fname, "c")
	f->RAY = RAY
;	
	
	
end