;; sstsnowleadlagcc.ncl Dr. Sujata Mandke  16-11-2018
;;--------------
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"        
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
begin 
;;;-----------------------------------
;; To calculate Lag-lead monthly correlation 
;;;;between monthly observed  OISST v2  sst 
;;;averaged over East Pacific(10S-10N;110W-80W) 
;;;and snow monthly averaged over Eurasia(50N-70N;20E-140E)
;;;
;;-------------------------------
   siglevl=0.05
;********************************
; Read observed monthly Snow   data
;********************************
   
   ymdStrt = 19820101   ; start yyyymmdd
   ymdLast = 20061201   ; end yyyymmdd
   
   diri    = "/moes/home/amin/cmip5/amitapap/data/" ; input directory   
   fili    = "swe_mon_jan1979tomay2007smmrssmi.nc" ; input snow file 
;;;********************************

   in1  = addfile(diri+fili ,"r")
   
   tm    = in1->time              ; all times on file
   ymd     = cd_calendar(tm, -2)  ; yyyymmdd
   
;;;;
   iStrt   = ind(ymd.eq.ymdStrt) ; index start
   iLast   = ind(ymd.eq.ymdLast) ; index last 
   
   delete(tm)
   delete(ymd)

;***********************************************************    
   sn= in1->swe(iStrt:iLast,{50:70},{20:140})  ;only specific region averaged over total Eurasia
   sneur=wgt_areaave(sn,1.0,1.0,0);;;create area averaged snow over Eurasia
;; 
    nr=dimsizes(sneur)
;************************************************
; Read monthly OISST v2   data 
;********************************
   ymdSt = 19820101   ; start yyyymmdd 
   ymdLa = 20061201   ; End yyyymmdd 

   vars     = "sst"   ; name of variable in file  
   dirs    = "/iitm3/erpas-res/skm/data/sst/oisst/mon/" ; input directory   
   fils    = "sst.mnmean.dec1981-oct2018.nc" ; input file 
;***********************************************************
; Read user specified time and create required yyyyddd                    
;***********************************************************
   in2      = addfile (dirs+fils , "r")                          

   tms    = in2->time             ; all times on file
   ymds   = cd_calendar(tms, -2)  ; yyyymmdd
;;;;
   iSts   = ind(ymds.eq.ymdSt)   ; index start
   iLas   = ind(ymds.eq.ymdLa)   ; index last 
      
   delete(tms)
   delete(ymds)

;***********************************************************   
   st= in2->sst(iSts:iLas,{-10:10},{250:280});only specific region East Pacific (110W-80W)ocean
   stepac=wgt_areaave(st,1.0,1.0,0)   ;;;create area averaged SST over East Pacific ocean
;************************************************
;************************************************
; calculate cross correlations for observed data 
;           Note: ccr1(0)=ccr2(0) 
;************************************************
   maxlag = 12                          ; set lag
   ccr1   = esccr(sneur,stepac,maxlag)  ; sneur lead stepac;calc positive lag cross cor: ccr1(maxlag+1)
   ccr2   = esccr(stepac,sneur,maxlag)  ; stepac lead sneur;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)  ; stepac lead sneur;negative lag;:::-1 means reverse the order 
   ccrtot(maxlag+1:)   = ccr1(1:)    ; sneur lead stepac;positive lag; ccr1(0:mxlagmaxlag-1)
 
   x = ispan(-maxlag,maxlag,1)       ; define x axis
;;---Follows testing significance of correlation---
;;
prob  = rtest(ccrtot, nr, 0)
;**************************
; plot the correlations and significance
;***************************
  wks = gsn_open_wks("eps","sstepacsnowleadlagcorsig")  ; send graphics to eps file
  
  res          = True
  res@tiMainString      = "Lead/lag correlation between Eurasian snow and East Pacific SST"  
  res@tmXBLabelFontHeightF = 0.027
  res@tmYLLabelFontHeightF = 0.027
  res@tmXBLabelFontAspectF = 0.9
  res@tmYLLabelFontAspectF = 0.9
;;-----------------------------
   res@tiMainString    = "Eurasian Snow-East Pacific SST lead-lag correlation"   
   res@tmXTOn  = False    ; turn off top xaxis tick marks
   res@tmYROn  = False    ; turn off right yaxis tick marks   
   res@tiXAxisString   = "LAG (Months)"  ; x-axis label
   res@tiYAxisString = "CORRELATION"     ; y-axis label
   res@trXMinF  = -maxlag  ; set minimum x-axis value
   res@trXMaxF  =  maxlag  ; set maximum x-axis value
;;;;
   plot1 = gsn_xy(wks,x,ccrtot,res)  ; plot correlation
;************************************************
end