;*************************************************
; 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
;************************************************
   in  = addfile("b003_TS_200-299.nc","r")
   ts = in->TS

;************************************************
; extract 2 arbitrary time series from 3d data 
; Two method: (i) standard indexing; (ii) coordinate subscripting
;************************************************

  ;ts1=ts(:,45,64)                               ; (i) extract time series via indexing
  ;ts2=ts(:,23,117)

   lat1 = 37.1
   lon1 = 180  
   lat2 = -23.5
   lon2 = 329.1
   ts1  = ts(:,{lat1},{lon1})                    ; extract time series via
   ts2  = ts(:,{lat2},{lon2})                    ; (ii) coordinate subscripting {...}

                                                 ; print variable overview
  ;printVarSummary(ts1)                          ; ts1@lat , ts1@lon are the selected grid pts
  ;printVarSummary(ts2)                          ; ts2@lat , ts2@lon

;************************************************
; calculate cross correlations 
;           Note: ccr1(0)=ccr2(0) 
;************************************************
   maxlag = 25                                     ; set lag
   ccr1   = esccr(ts1,ts2,maxlag)                  ; calc positive lag cross cor: ccr1(maxlag+1)
   ccr2   = esccr(ts2,ts1,maxlag)                  ; calc negative lag cross cor: ccr2(maxlag+1)
  ;print("---")                 
  ;print("ccr1="+sprintf("%7.3f", ccr1)+"  ccr2="+sprintf("%7.3f", ccr2) \
  ;             +"  ccr2(::-1)="+sprintf("%7.3f", ccr2(::-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