load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/diagnostics_cam.ncl"
; import data
f=addfile("/mnt/geog/k/kg/kg312/Exercise/onset_x.nc","r")
x=f->precip

; SJ onset
; Perform a 10-day running average
x_sm_SJ=runave_n_Wrap(x,10,0,0) ; this implies that 4 points will be lost at the leftsmost side and 5 at the rightmost side
; zonal average
x_zonal_SJ=dim_avg_n_Wrap(x_sm_SJ,2)
;Interpolating data to pass from a (0.5,1.5,2.5,...) grid to a (0,1,2,3,...) grid
xi_SJ=x_zonal_SJ&lat ; get the latitude coordinates of xi_SJ(input coordinates)
xo_SJ=ispan(0,21,1)  ; define the new latitude coordinates (output coordinates)
copy_VarMeta(xi_SJ,xo_SJ)
xo_SJ=dim_pqsort( xo_SJ,-1); Sort in the same order as xi_SJ
x_regrid_SJ=linint1_n_Wrap(xi_SJ,x_zonal_SJ,False,xo_SJ,0,1); False is to say not to process cyclically, 0 is set by default; 1 means the dimension of the latitude in x_zonal_SJ

; Change the type of lat cordinate of x_regrid_SJ to double
Lat=int2dble(x_regrid_SJ&lat)
delete(x_regrid_SJ&lat)
x_regrid_SJ&lat=Lat

; Changing time to yyyymmdd format
time=f->time
TIME=cd_calendar(time,-2)
;time in fraction of year
yrfrac=yyyymmdd_to_yyyyfrac(TIME,0)
x_regrid_SJ&time=yrfrac
;Extract the year
TIM=cd_calendar(time,0)
year=floattointeger(TIM(0,0))
;Plotting
wks=gsn_open_wks("png", "time_lat"+tostring(year) )
res                  = True                ; plot mods desired
res@tiMainString     = "Sultan_Janicot_"+tostring(year)     ; title
res@tmXBLabelStride  = 2                   ; tick mark label stride
res@tiYAxisString    = "time"              ; y axis title

res@cnFillOn         = True                ; color on
res@lbLabelStride    = 1                   ; every other label
res@lbOrientation    = "Vertical"          ; vertical label bar
res@cnLinesOn        = True               ; turn off contour lines
;res@cnFillPalette    = "ViBlGrWhYeOrRe"    ; set color map
res@cnLevelSpacingF  = 2                   ; contour spacing
res@sfXArray =x_regrid_SJ&lat
res@sfYArray =x_regrid_SJ&time

plot=gsn_csm_lat_time(wks, x_regrid_SJ, res)