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"
begin
;  =================================
;   read lon-lat of 160 stations
;  =================================
;  lonlat=asciiread("/home/linx/work/160st/160LIST2.TXT",(/160,2/),"float")
;   lon   =lonlat(:,0)
;   lat   =lonlat(:,1)
;  =================================
;   read 160stations' obs. rainfall
;  =================================
   pt    = asciiread("/home/linx/work/160st/r16004.02",(/160,639/),"float")  ;from 1951.1 to 2004.2  53 yrs
   p     = pt(:,1:)                                    ;unit= mm,the first column is number of stations

;print(lon)  ; check

   xj    = new((/4,53/),"float")
   xj(0,:) = month_to_annual(p(150,0:635),0)   ; kashi
   xj(1,:) = month_to_annual(p(151,0:635),0)   ; hetian
   xj(2,:) = month_to_annual(p(152,0:635),0)   ; qiemo
   xj(3,:) = (xj(0,:)+xj(1,:))/2.0

   time = ispan(1951,2003,1)
   xj!1 = "time"
   xj&time = time


;  fi = addfile("/home/linx/work/160st/p160.nc","r")
   fi = addfile("/home/linx/CRU/CRUpre1901-2002.nc","r")
   printVarSummary(fi) 

   cru = fi->pre(time|600:1223,lat|:,lon|:)     ; load the CRU PRCP data from 1951 to 2002
   printVarSummary(cru) 
  
   xp160 = month_to_annual(cru,0)   ; time should be the leftmost dim when using month_to_annual
   xp160!0 = "time"
   printVarSummary(xp160)
   delete(cru)

; ************************************************
; reorder to get time as right most dimension
;***********************************************

   ann_cru = xp160(lat|:,lon|:,time|:)
   printVarSummary(ann_cru)
   delete(xp160)

   ccr1 = escorc(xj(3,0:51),ann_cru)   ; one-point correlation pattern
   prob1 = rtest(ccr1,52,0)

   ccr2 = escorc(xj(3,0:50),ann_cru(:,:,1:51))
   prob2 = rtest(ccr2,51,0)

   ccr3 = escorc(xj(3,1:51),ann_cru(:,:,0:50)) 
   prob3 = rtest(ccr3,51,0)

 ; printVarSummary(ccr)

;************************************************
; calculate cross correlations
;************************************************
;  maxlag = 2                             ; set lag
; note, the max lag should not be more than N/4

;  ccr = esccr(ann_cru,xj(3,0:51),maxlag)   ; calc the cross correlation

; copy meta data and coordinate variables using contributed functions
  copy_VarAtts(ann_cru, ccr1)                    
  copy_VarCoords_1(ann_cru,ccr1)
  copy_VarAtts(ann_cru, ccr2)
  copy_VarCoords_1(ann_cru,ccr2)
  copy_VarAtts(ann_cru, ccr3)
  copy_VarCoords_1(ann_cru,ccr3)


;************************************************
; plot the correlations
;************************************************
  res = True
  res@gsnDraw	= False
  res@gsnFrame	= False
  res@tfPolyDrawOrder    = "Predraw"  
  
  res@vpWidthF	= 2.0
  res@vpHeightF	= 0.6
  res@gsnPaperOrientation	= "Portrait"
  res@gsnMaximize	= True
  res@gsnPaperWidth	= 10
  res@gsnPaperHeight	= 40 

plot = new(3,graphic)
plotv = new(3,graphic)

wks = gsn_open_wks("eps","XinJiang.ann.PRCP.vs.CRU.corr")

gsn_define_colormap(wks,"BlWhRe")              ; choose colormap

res                      = True                ; make plot mods

  res@cnFillOn             = True                ; turn on color
  res@gsnSpreadColors      = True                ; use full colormap
  res@cnLinesOn            = False               ; turn off contour lines
  res@lbLabelAutoStride    = True                ; automatic lb label stride

  res@cnLevelSelectionMode = "ManualLevels"      ; manually set cn levels
  res@cnMinLevelValF       = -0.8                 ; min level
  res@cnMaxLevelValF       =  0.8                 ; max level
  res@cnLevelSpacingF      = .1                  ; contour level spacing

  res@gsnLeftString        = False
  res@gsnRightString       = False

  res@tiMainString         = "ANN PRCP Correlations: South XinJiang vs Global(CRU) "
  res@gsnCenterString      = "South Xinjiang lead 1yr"
  plot(0) = gsn_csm_contour_map_ce(wks,ccr2(:,:),res)  ; lead 1  
  res2    = True
  res2@cnLevelSelectionMode = "ManualLevels"  
  res2@cnMinLevelValF       = 0.05                  ; min level
  res2@cnMaxLevelValF       = 0.1                   ; max level
  res2@cnLevelSpacingF      = 0.05                  ; contour level spacing
  plotv(0)= gsn_csm_contour(wks, prob2, res2)
  overlay(plot(0),plotv(0)) 

  res@tiMainString         = False
  res@gsnCenterString      = "Same Year"
  plot(1) = gsn_csm_contour_map_ce(wks,ccr1(:,:),res)  ; lead 0 
  plotv(1)= gsn_csm_contour(wks, prob1, res2)
  overlay(plot(1),plotv(1))

  res@tiMainString         = False
  res@gsnCenterString      = "South Xinjiang lag 1yr"
  plot(2) = gsn_csm_contour_map_ce(wks,ccr3(:,:),res)  ; lag 1 
  plotv(2)= gsn_csm_contour(wks, prob3, res2)
  overlay(plot(2),plotv(2))
;********************************************************************************
  resp	= True
  resp@gsnFrame	= False
  gsn_panel(wks,plot,(/3,1/),resp)
  frame(wks)
 


end