begin  
;  f      = addfile ("../SST/sst.mnmean.nc", "r")
;
;  TIME   = f->time
;  YYYY   = cd_calendar(TIME,-1)/100      ;cd_calendar (time, option) ;option = -1 or 1. The values returned will be in the format YYYYMM and will be type double for option = 1, and type integer for option = -1.        
  ;print(YYYY)
  
  ; define starting and ending year
  yrStrt = 1980
  yrLast = 2021
  
;  iYYYY  = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)  ;Returns the indices where the input is True. A one-dimensional logical array.
  ;print(iYYYY)
 
;  sst    = f->sst(iYYYY,:,:)
  ;printVarSummary(sst)  

  b = addfile("/project/mojave/observations/OBS-SST/hadisst.187001-202212.nc","r")
  sstf = b->sst
  sst = sstf({yrStrt*100+1:yrLast*100+12},:,:)
  
  sstSON = month_to_season (sst, "SON")
  ;printVarSummary(sstSON)
  
  ano_sstSON = dim_rmvmean_n_Wrap(sstSON,0)
  ;printVarSummary(ano_sstSON)


  lines     = asciiread("AMM_DMI_NINO_PDO_ELI_SepNovAvg_1980_2021.csv",-1,"string")
  data      = lines(1:)           ; remove header line
 
  nino      = tofloat(str_get_field(data,3,","))
  ;print(nino)
 
 
  ;-----------------------------------
  ; lag-0 correlation      
  ;-----------------------------------
  
  ano_nino = nino - avg(nino)
  copy_VarCoords(nino, ano_nino)
  
  r_sst  = escorc(ano_nino,ano_sstSON(lat|:,lon|:,time|:))        ; ccr(lat,lon)
  copy_VarCoords(ano_sstSON(0,:,:), r_sst)
  r_sst@long_name = "Correlation: r_sst"

  p = rtest(r_sst,42,0)
  p@_FillValue = 1.e20
  copy_VarCoords(r_sst, p)
  p = where(ismissing(r_sst),p@_FillValue,p)

  psig = 0.05                       ; test significance level                     
  
  if (any(p.le.psig)) then
      ;print(p)
      print("r is significant at the 95% level")
  else
      print("r is NOT significant at the 95% level")
  end if
  

  pltType = "png"            
  pltDir  = "./"
  pltName = "sst_SepNov_cor"
  pltPath = pltDir + pltName 
  
  pltType@wkWidth  = 2500
  pltType@wkHeight = 2500 

  wks  = gsn_open_wks(pltType,pltPath)


  res                      = True         
  res@gsnDraw              = False        ; don't draw yet
  res@gsnFrame             = False        ; don't advance frame yet
  res@gsnAddCyclic         = False        ; plotted dataa are not cyclic
  
  res@mpFillOn             = False        ; turn off map fill
  res@mpFillDrawOrder      = "PostDraw"          ; color of land
  res@mpLandFillColor      = "grey"              ; color of land
  res@mpGeophysicalLineThicknessF = 4
  res@mpGeophysicalLineColor      = "blue"
  res@mpNationalLineColor         = "blue"
  res@mpOutlineBoundarySets       = "AllBoundaries"
  res@pmTickMarkDisplayMode       = "Always"
  
  res@cnFillOn             = True         ; turn on color fill
  res@cnLinesOn            = False        ; True is default
  res@cnLineLabelsOn       = False            ; turn off line labels
  res@cnInfoLabelOn        = False            ; turn off info label
  
  res@cnFillMode           = "RasterFill"
  res@cnRasterSmoothingOn  = True
  res@cnFillPalette        = "BlWhRe"     ; set color map
  res@lbLabelBarOn         = True        ; turn off individual lb's
  
  
  ;regional plot
  ;res@mpLimitMode          = "LatLon"
  ;res@mpMinLatF            = 0         
  ;res@mpMaxLatF            = 30
  ;res@mpMinLonF            = 100
  ;res@mpMaxLonF            = 180
  res@mpCenterLonF         = 180.
  
  
  res2 = True                            ; res2 probability plots
  res2@gsnDraw             = False       ; Do not draw plot
  res2@gsnFrame            = False       ; Do not advance frome
  res2@cnFillOn            = True
  

  res2@cnLevelSelectionMode = "ExplicitLevels" ; set manual contour levels
  res2@cnLevels = (/0.05/)
  res2@cnInfoLabelOn       = False
  res2@cnLinesOn           = False      ; do not draw contour lines
  res2@cnLineLabelsOn      = False       ; do not draw contour labels
  res2@lbLabelBarOn        = False        
  
  res2@gsnLeftString       = ""
  res2@gsnRightString      = ""
  res2@cnMonoFillPattern = False
  res2@cnFillPatterns = (/17,0/)
  res2@cnFillColors = (/"green","transparent"/)   
  res2@cnFillScaleF       =  0.6  ;5
  res2@cnFillDotSizeF = .0020     ; .0001
  res2@cnFillDrawOrder = "Draw"
  
  
  opt                      = True
  opt@gsnShadeFillType     = "pattern"
  opt@gsnShadeLow          = 17          ; stipple 
  ;opt@gsnShadeHigh         = 4         
  opt@gsnShadeFillScaleF   = 0.25         ; add extra density
  opt@gsnShadeFillDotSizeF = 0.0015
 
  plot = gsn_csm_contour_map(wks,r_sst, res)  ; create plot 
  
  plot2 = gsn_csm_contour(wks,p, res2)  
;  plot2 = gsn_contour_shade(plot2,0.05,999.,opt)\
  
  overlay (plot, plot2)
 
 draw(plot)
 frame(wks)
end