;*************************************************

; corel_2.ncl

;

; Concepts illustrated:

;   - Calculating positive and negative lags in a cross correlation

;   - Illustrate using coordinate subscripting syntax [ {,,,} ] to select locations

;   - Illustrate reversing array order via the  ::-1 syntax

;

;************************************************

;

; These files are loaded by default in NCL V6.2.0 and newer

; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  

; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin

;************************************************

; variable and file handling

;************************************************

   in1  = addfile("pc1.nc","r")
   in2  = addfile("pc2.nc","r")

   pc1 = in1->precip
   pc2 = in2->precip


;************************************************

; calculate cross correlations 

;           Note: ccr1(0)=ccr2(0) 

;************************************************

   maxlag = 25                                     ; set lag

   ccr1   = esccr(pc1,pc2,maxlag)                  ; calc positive lag cross cor: ccr1(maxlag+1)

   ccr2   = esccr(pc2,pc1,maxlag)                  ; calc negative lag cross cor: ccr2(maxlag+1)

;************************************************

;  combine pos and neg into one time series      

;          do not repeat the ccr(0) index

;************************************************

   totlag =2*maxlag +1                           ; set total lag



   ccrtot = new( totlag, typeof(ccr1) )          ; allocate memory 

   ccrtot(0:maxlag )   = ccr2(::-1)              ; ::-1 means reverse the order 

   ccrtot(maxlag+1:)   = ccr1(1:)                ; ccr1(0:mxlagmaxlag-1)

 

   x = ispan(-maxlag,maxlag,1)                   ; define x axis



  ;print("---")

  ;print("x="+sprintf("%7.3f", x)+"  ccrtot="+sprintf("%7.3f", ccrtot))



;************************************************

; plot the correlations

;************************************************

   wks  = gsn_open_wks("png","corel")                  ; send graphics to PNG file

   res                 = True                         ; make plot mods

   res@tiMainString    = "("+sprintf("%5.1f",ts1@lat)+","+sprintf("%5.1f",ts1@lon)+") and " \ 
                       + "("+sprintf("%5.1f",ts2@lat)+","+sprintf("%5.1f",ts2@lon)+")"  

   res@tiXAxisString   = "LAG"                        ; x-axis label



   res@trXMinF  = -maxlag

   res@trXMaxF  =  maxlag



   plot = gsn_xy(wks,x,ccrtot,res)                       ; plot correlation

end