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/0K0KFIXVEG/ha.foam1-5lpj.0k0kfixveg.anni6351034.nc","r") ; open file 1
 gw = a->gw
 lat = a->lat
 lon = a->lon
 time1 = a->time

  sfct1   = (a->T(:,17,:,:))-273.15     ; retrieve sfc air T from file a
  sfct1!0="time"
  sfct1!1="lat"
  sfct1!2="lon"
 sfct1&lat=lat
 sfct1&lon=lon
 temp1 = new( (/40,48,400/),float)
 temp1 = sfct1(lat|:,lon|:,time|:)     ; reorder dims
 temp1!0 = "lat"
 temp1!1 = "lon"
 temp1!2 = "time"
 temp1&lat = lat
 temp1&lon = lon
 temp1@time = time1
 temp1@_FillValue = 1e+35

 avtemp1=dim_avg(temp1)    ; av temp
 avtemp1!0="lat"
 avtemp1!1="lon"
 avtemp1&lat=lat
 avtemp1&lon=lon

 var1=dim_variance(temp1)  ; var of temp
 var1!0="lat"
 var1!1="lon"
 var1&lat=lat
 var1&lon=lon

;;;;;;;;;;;;;;;;;;;;;

  b = addfile("/ptmp/notaro/run/ha.foam1-5lpj.ctl-pi-test.anni6351034.nc","r")   ; open file 2
 time2 = b->time

  sfct2   = (b->T(:,17,:,:))-273.15  ; retrieve sfc air T from file b
  sfct2!0="time"
  sfct2!1="lat"
  sfct2!2="lon"
 sfct2&lat=lat
 sfct2&lon=lon
 temp2 = new( (/40,48,400/),float)
 temp2 = sfct2(lat|:,lon|:,time|:)  ; reorder dims
 temp2!0 = "lat"
 temp2!1 = "lon"
 temp2!2 = "time"
 temp2&lat = lat
 temp2&lon = lon
 temp2@time = time2
 temp2@_FillValue = 1e+35

 avtemp2=dim_avg(temp2)           ; average temp
 copy_VarCoords(avtemp1,avtemp2)

 var2=dim_variance(temp2)   ; variance of temp
 copy_VarCoords(var1,var2)

;;;;;;;;;;;;;;;

 avdiff=new((/40,48/),float)  ; av diff in temp
 avdiff=avtemp2-avtemp1
 copy_VarCoords(var1,avdiff)

 print(max(avdiff))
 print(min(avdiff))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  Eqv1 = equiv_sample_size (temp1, 0.1,0) ; critical sig lvl of 0.1 (equiv sample size)
  Eqv2 = equiv_sample_size (temp2, 0.1,0)
  N1   = wgt_areaave (Eqv1, gw, 1., 0)    
  N2   = wgt_areaave (Eqv2, gw, 1., 0) 
  prob = ttest(avtemp1,var1,N1,avtemp2,var2,N2, True, True)  ; t-test
  prob!0="var"
  prob!1="lat"
  prob!2="lon"
  prob&lat=lat
  prob&lon=lon
;;;;;;;;;;;;;;;;;;;;;;

 wks   = gsn_open_wks ("ps", "mean_diff" )            ; open workstation
 plot = new(1,"graphic")

 res          = True                             ; plot mods desired
 gsn_define_colormap(wks,"wxpEnIR")  ; choose color map
 res@gsnDraw=False
 res@gsnFrame=False
 res@gsnAddCyclic = True
 res@tiMainString        = "Diff in Av Sfc Air T(C) FOAM-LPJ Run A(400Yrs)-Run B(400Yrs)"
 res@tiMainFontHeightF = 0.011
 res@lbPerimOn          = False
 res@mpPerimOn          = True
  res@mpOutlineOn      = True           ; turn on map outline
  res@cnLineThicknessF = 1.5                      ; thicker lines
  res@lbLabelAutoStride    = True
  res@tmYROn=False
  res@tmXTOn=False
  res@cnInfoLabelOn         = False
  res@cnLineLabelsOn       = True              ; turn off contour line labels
  res@cnLineLabelPlacementMode = "computed"
  res@cnLineLabelInterval = 1
  res@cnLineLabelPerimOn         = False           ; line label box on
  res@cnLevelSelectionMode = "ExplicitLevels"   ; set explicit contour levels
  res@cnLevels            = (/-0.20,-.16,-.12,-.08,-.04,.04,.08,.12,.16,.20/)
  res@cnMonoLineDashPattern = False
  res@cnLineDashPatterns = (/2,2,2,2,2,0,0,0,0,0/)
  res@cnLineLabelFontHeightF = 0.009
  res@cnFillOn=False
  res@cnLineThicknessF = 4.0
  res@lbLabelBarOn        = False
  res@pmLabelBarDisplayMode =  False
  res@mpLandFillColor = -1
  res@mpOceanFillColor = -1
  res@mpCenterLonF = 0.
  plota1 = gsn_csm_contour_map_ce(wks,avdiff,res) 

  resa = True                            ; res2 probability plots
  resa@gsnDraw             = False       ; Do not draw plot
  resa@gsnFrame            = False       ; Do not advance frome
  resa@gsnAddCyclic = True
  resa@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
  resa@cnMinLevelValF      = 0.00        ; set min contour level
  resa@cnMaxLevelValF      = 0.15        ; set max contour level
  resa@cnLevelSpacingF     = 0.05        ; set contour spacing
  resa@cnInfoLabelOn       = False
  resa@cnLinesOn           = False       ; do not draw contour lines
  resa@cnLineLabelsOn      = False       ; do not draw contour labels
  resa@cnFillScaleF        = 0.9         ; add extra density
  plota2   = gsn_csm_contour(wks,prob(0,:,:), resa) 
  plota2   = ShadeLtContour(plota2, 0.1, 16)   ; shade signif <0.1
  overlay (plota1, plota2)
  plot(0)=plota1

  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@gsnPaperOrientation="portrait"
  pres@gsnPaperMargin=0.1
  pres@gsnPanelLeft = 0.02
  pres@gsnPanelRight = 0.98
  pres@gsnPanelBottom = 0.02
  pres@gsnPanelTop = 0.98
  pres@gsnPanelRowSpec=True
  gsn_panel(wks,plot,(/1/),pres)        ; create panel plot

end