;***************************************************************
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

;************************************************
; open file and read in variable
;***********************************************
  fils1 = ("ts.cfsr.nc")
  f1    = addfile(fils1, "r") 
  x1    = f1->ts

x1@_FillValue = 1e20
x1=where(x1.ge.-3, x1, x1@_FillValue)
printVarSummary( x1 ) 
                      
  fils2 = ("ts.bcc.nc")
  f2    = addfile(fils2, "r")  
  t     = f2->ts
 
x2=t-273.15
copy_VarCoords(t, x2)
x2@_FillValue = 1e20
x2=where(x2.ge.-3, x2, x2@_FillValue)

printVarSummary( x2 )
printMinMax( x2, True )

  
; Compute climatologies and anomalies. This removes the annual cycle. 

  tmp1Clm = clmMonTLL(x1)
  tmp1Anom= calcMonAnomTLL(x1, tmp1Clm)
  printVarSummary(tmp1Anom)

  tmp2Clm = clmMonTLL(x2)
  tmp2Anom= calcMonAnomTLL(x2, tmp2Clm)
  printVarSummary(tmp2Anom)
  
; For convenience, reorder and store the data

  TMP1 = tmp1Anom(lat|:,lon|:,time|:)
  TMP2 = tmp2Anom(lat|:,lon|:,time|:)
 
  delete( [/tmp1Anom, tmp2Anom /] )   ; no longer needed

; Compute cross correlation coefficients at 0-lag

  r   = escorc( TMP1, TMP2 )   ; (nlat,mlon)

  copy_VarCoords(TMP1(:,:,0),r)
  r@long_name = "linear cross correlation"
  printVarSummary(r)
  printMinMax( r, True )

; Calculate degrees of freedom for each time series pair
; This estimates the dof in autocorrelated series. 
  
  siglvl = 0.05
  dof1   = equiv_sample_size (TMP1, siglvl,0)
  dof2   = equiv_sample_size (TMP2, siglvl,0)

  printMinMax(dof1, 0)
  printMinMax(dof2, 0)
  
 ; For each grid point time series pair, choose the lesser abs value
 
  dof  = where(abs(dof1).le.abs(dof2), dof1, dof2) 

  copy_VarCoords(r,dof)
  dof@long_name = "degrees of freedom"
  printVarSummary(dof)
  printMinMax(dof,0)
  
  rsig = rtest( r, round(dof,3), 0)

  rsig@long_name = "significance of r"
  copy_VarCoords(r,rsig)
  printVarSummary(rsig)
  printMinMax(rsig, 0)

; ========================= PLOT 1 ==============================

  wks   = gsn_open_wks ("png", "ts.cfsrxmiroc" )   ; open workstation 

  gsn_define_colormap(wks,"cmp_b2r")
 
  res                      = True                ; make plot mods


  res@cnFillOn             = True                ; turn on color

  ;res@gsnAddCyclic = False

  res@gsnSpreadColors      = True                ; use full colormap


  res@lbLabelAutoStride    = True                ; automatic lb label stride

  res@cnLinesOn            = False               ; turn off contour lines

  res@lbLabelFontHeightF   = 0.015           ;aumenta numeros label barra
  res@cnLevelSelectionMode = "ManualLevels"  ; set manual contour levels
  res@cnMinLevelValF       = -1.             ; set min contour level
  res@cnMaxLevelValF       =  1.             ; set max contour level
  res@cnLevelSpacingF      =  0.1            ; set contour spacing

  res@gsnCenterString      = "CFSR x MIROC5"                               
  res@gsnStringFontHeightF = 0.025
 
          
  res@mpLandFillColor      = "white"
  res@gsnRightString   = ""    ;
  res@gsnLeftString    = ""    ; 
  res@mpDataBaseVersion = "HighRes"
  plot = gsn_csm_contour_map_ce(wks,r, res)  

  plot = ZeroNegDashLineContour (plot)
; ========================= PLOT 2 ==============================
  res2 = True                            ; res2 probability plots

  res2@gsnDraw             = False       ; Do not draw plot
  res2@gsnFrame            = False       ; Do not advance frome

  ;res2@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
  ;res2@cnMinLevelValF      = 0.0          ; set min contour level
  ;res2@cnMaxLevelValF      = 1.0          ; set max contour level
  ;res2@cnLevelSpacingF     = 0.05        ; set contour spacing

  res2@cnInfoLabelOn       = False       ; turn off info label

  res2@cnLinesOn           = False       ; do not draw contour lines
  res2@cnLineLabelsOn      = False       ; do not draw contour labels

  plot2   = gsn_csm_contour(wks,gsn_add_cyclic_point(rsig), res2) 
  plot2   = ShadeGtContour(plot2, 0.95, 17)  ; shade all areas more than the
                                             ; 0.05 contour level
  overlay (plot, plot2)

  draw (plot)
  frame(wks)
end