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/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
 a = addfile("/ptmp/notaro/run/ha.foam1-5lpj.ctl-pi-test.anni6351034.nc","r") ; open data file
 lat  = a->lat  ; retrieve lat, lon, time
 lon  = a->lon
 time = a->time

 temp = a->T(:,17,:,:)               ; retrieve sfc air T
 sfct = new( (/40,48,400/),float)
 sfct = temp(lat|:,lon|:,time|:)     ; reorder dims
 sfct!0 = "lat"
 sfct!1 = "lon"
 sfct!2 = "time"
 sfct&lat = lat
 sfct&lon = lon
 sfct&time = time

 precip = (a->PRECL)+(a->PRECC)    ; retrieve precip
 copy_VarCoords(temp,precip)       ; copy attributes from temp to precip
 pcp = new( (/40,48,400/),float)
 pcp = precip(lat|:,lon|:,time|:)  ; reorder dims
 copy_VarCoords(sfct,pcp)

 corr=new((/40,48/),float)
 corr=escorc(sfct,pcp)             ; correlation
 corr!0 = "lat"
 corr!1 = "lon"
 corr&lat = lat
 corr&lon = lon
 corr@_FillValue=1e+35
 print(max(corr))
 print(min(corr))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  wks = gsn_open_wks("ps","correl")          ; open a ps file
  plot=new(1,graphic)
  res                   = True      ; plot mods desired
  res@gsnDraw           = False      ; don't draw
  res@gsnFrame          = False     ; don't advance frame
  res@cnFillOn         = True               ; color Fill
  res@cnRasterModeOn   = False               ; Raster Mode
  res@cnLinesOn        =  False             ; Turn off contour lines
  res@cnLevelSelectionMode = "ExplicitLevels" ; set manual contour levels
  res@cnLevels = (/-0.65,-0.5,-0.35,-0.2,0.,0.2,0.35,0.5,0.65/)
  res@tiMainFontHeightF =  0.012
  res@tiMainString         = "Correl Annual Sfc Air T vs. Pcp (IntVeg 400Yrs)"
  res@mpOceanFillColor     = -1
  colors = \
   (/"white","black","midnight blue","light sky blue","green","pale green","white",\
     "white","yellow","orange","pink","red"/)
   gsn_define_colormap(wks,colors)
  res@vpXF            = 0.06      
  res@vpWidthF        = 0.78       
  res@vpYF            = 0.93
  res@vpHeightF       = 0.39
  res@pmLabelBarDisplayMode = False     ; Turn on a label bar.
  res@lbOrientation        = "vertical"
  res@tmXTOn = False
  res@tmYROn = False
  res@mpMinLatF            = -60.  
  res@mpMaxLatF            =  75.
  res@mpPerimOn          = True
  res@mpOutlineOn      = True           ; turn on map outline
  plot = gsn_csm_contour_map_ce(wks,corr,res) ; create plot

  pres             = True                 ; mod panel plot
  pres@gsnDraw=True
  pres@gsnFrame=True
  pres@gsnMaximize = True                        ; blow up plot
  pres@gsnPaperWidth=8.5
  pres@gsnPaperHeight=11.0
  pres@gsnPaperMargin=0.07
  pres@gsnPanelLeft = 0.04
  pres@gsnPanelRight = 0.96
  pres@gsnPanelTop = 0.96
  pres@gsnPanelBottom = 0.04
  pres@gsnPanelYWhiteSpacePercent = 3
  pres@gsnPanelRowSpec=True
  gsn_panel(wks,plot,(/1/),pres)        ; create panel plot

end