;------ cross correlation of zonal precipitation 1951-2000

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

;start and end dates for 10-year windows - CanESM2, CCSM4, CNRM-CM3 and GISS-E2-H
;R script used to get this info: gcm_sim_historical_average_timeseries_over_region.R

;   startDate startDate_index    endDate endDate_index
; 1951-01-01            1213  1960-12-01          1332
; 1961-01-01            1333  1970-12-01          1452
; 1971-01-01            1453  1980-12-01          1572
; 1981-01-01            1573  1990-12-01          1692
; 1991-01-01            1693  2000-12-01          1812

;start and end dates for 10-year windows - HadCM3
;   startDate startDate_index    endDate endDate_index
; 1951-01-01            1094  1960-12-01          1213
; 1961-01-01            1214  1970-12-01          1333
; 1971-01-01            1334  1980-12-01          1453
; 1981-01-01            1454  1990-12-01          1573
; 1991-01-01            1574  2000-12-01          1693

; start and end date indices for models CanESM2, CCSM4, CNRM-CM3, GISS-E2-H
stDate1 = (/195101,196101,197101,198101,199101/)
enDate1 = (/196012,197012,198012,199012,200012/)

;----- data locations
  diri = (/ "/home/fas/saiers/nra23/scratch/climate/global/forcings/1x1deg/monthly/prec/",\
            "/home/fas/saiers/nra23/scratch/climate/global/gcm/IPCC_5/historical/monthly/NorESM1-M/"/)

;------ file names
  fname=(/"prec_monthly_1948-2008.nc",\
          "pr_Amon_NorESM1-M_historical_r1i1p1_185001-200512.nc"/)

  title=(/"Observed","r1i1p1"/)

  out_fname = (/"pr_crosscorr_zonal_NorESM1-1948-2005"/)

;-------- define lat-lon box
  lat1 = -16.
  lat2 = 16.
  lon1 = 12.
  lon2 = 40.

  wks_type = "ps"
  wks_type@wkPaperSize = "letter"
  wks_type@wkOrientation = "landscape"
  wks   = gsn_open_wks(wks_type ,"test-crosscorr")

  gsn_define_colormap(wks, "precip3_16lev")
  plot = new(1,graphic)

;---- read the observed data file
  f      = addfile (diri(0) + fname(0), "r")
  lat    = f->lat
  lon    = f->lon
  time   = f->time
  time1  = cd_calendar(time, -1)
;---- observed data name is different from others
     P   = f->data
     P   = P * 30.		; units are in mm/day
     P@units  = "mm/month"

;------- find indices for the lat-lon box
  lat1_id = ind(lat1.gt.lat)
  lat2_id = ind(lat2.lt.lat)
  lon1_id = ind(lon1.gt.lon)
  lon2_id = ind(lon2.lt.lon)
  ilt1    = lat1_id(dimsizes(lat1_id)-1)
  iln1    = lon1_id(dimsizes(lon1_id)-1)
  ilt2    = lat2_id(0) + 1
  iln2    = lon2_id(0) + 1
;------- find time indices
  t1      = ind(stDate1(0).eq.time1)
  t2      = ind(enDate1(0).eq.time1)
;------- use functions in contributed.ncl, time must be multiple of 12
  Pclm = clmMonTLL(P(t1:t2,ilt1:ilt2,iln1:iln2))
  Pzone1 = dim_avg_Wrap(Pclm(lat|:,month|:,lon|:))
;------ delete variables for reuse for the next dataset
delete([/lat,lon,time,time1,f,P/])
delete([/lat1_id,lat2_id,lon1_id,lon2_id,ilt1,ilt2,\
         iln1,iln2,t1,t2,Pclm/])

;---- read gcm output
  f    = addfile (diri(1) + fname(1), "r")
  lat = f->lat
  lon = f->lon
  time = f->time
  time1 = cd_calendar(time, -1)
; observed data name is different from others. Correct for this
   P    = f->pr
   P    = P * 86400. * 30.	; units are in kg/m2/sec
   P@units  = "mm/month"
;------- find indices for the lat-lon box
  lat1_id = ind(lat1.gt.lat)
  lat2_id = ind(lat2.lt.lat)
  lon1_id = ind(lon1.gt.lon)
  lon2_id = ind(lon2.lt.lon)
  ilt1    = lat1_id(dimsizes(lat1_id)-1)
  iln1    = lon1_id(dimsizes(lon1_id)-1)
  ilt2    = lat2_id(0) + 1
  iln2    = lon2_id(0) + 1
;------- find time indices
  t1      = ind(stDate1(0).eq.time1)
  t2      = ind(enDate1(0).eq.time1)

;------- use functions in contributed.ncl, time must be multiple of 12
  Pclm = clmMonTLL(P(t1:t2,ilt1:ilt2,iln1:iln2))
  Pzone2 = dim_avg_Wrap(Pclm(lat|:,month|:,lon|:))

;------- cross correlation
  ccr  = escorc(Pzone1,Pzone2)

;------ delete variables 
delete([/lat,lon,time,time1,f,P/])
delete([/lat1_id,lat2_id,lon1_id,lon2_id,ilt1,ilt2,\
         iln1,iln2,t1,t2,Pclm/])

;------ plot resources
  res                      = True

  res@gsnDraw              = False
  res@gsnFrame             = False
  res@gsnMaximize          = True
  res@cnFillOn             = True            ; turn on color
  res@gsnSpreadColors      = True            ; use full color table
  res@cnLinesOn            = False           ; no contour lines
  res@cnInfoLabelOn        = False
  res@pmTickMarkDisplayMode  = "Always"      ; turn on fancy tickmarks
  																																																																																																																																																																			
  res@tmXBMode             = "Explicit"      ; label independently
  res@tmXBValues           = ispan(0,11,1)
  res@tmXBLabels           = (/"J","F","M","A","M","J","J",\
                               "A","S","O","N","D"/)
  
  res@cnLevelSelectionMode = "ManualLevels"
  res@cnLevelSpacingF      = .1
  res@cnMinLevelValF       = -1.
  res@cnMaxLevelValF       = 1.
  res@lbLabelBarOn         = False

  res@tiMainString         = title(0)
  res@tiMainFontHeightF    = 0.02
  res@gsnLeftString = ""
  res@gsnRightString = ""
 
  plot(0) = gsn_csm_lat_time(wks, ccr, res) 

;---- panel plot create
  pres                    = True         ; Set panel resources.
  pres@gsnPanelLabelBar   = True         ; Turn on panel labelbar.
  pres@lbLabelStride      = 2            ; print every label marker
  pres@lbLabelFontHeightF = 0.010        ; label font height
  pres@gsnFrame           =False         ; do not advance the frame
  gsn_panel(wks,plot,(/1,2/),pres)

;  drawNDCGrid(wks)

  txres1 = True
  txres1@txFontHeightF = 0.012
  gsn_text_ndc(wks,"NorESM1-M",0.15,1.1,txres1)

  txres2 = True
  txres2@txFontHeightF = 0.012
  gsn_text_ndc(wks,"Precipitation (mm/month)",0.5,0.001,txres2)

  frame(wks)

; remove plot and wks objects
  destroy(plot)
  destroy(wks)

end