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/kutz/ha.foam1-5lpj.ctl-pi-test.SLP.ann.6351034.nc","r") ; open data file 1
 lat = a->lat     ; retrieve lat, lon, time from file a
 lon = a->lon
 time1 = a->time

  psl1   = a->psl  ; retrieve slp from file a
  slp1 = new( (/40,48,400/),float)
  slp1 = psl1(lat|:,lon|:,time|:)   ; reorder slp dimensions
 slp1!0 = "lat"
 slp1!1 = "lon"
 slp1!2 = "time"
 slp1&lat = lat
 slp1&lon = lon
 slp1@time = time1
 slp1@_FillValue = 1e+35

 var1 = new((/40,48/),float)
 var1 = dim_variance(slp1)  ; calculate variance
 var1!0 = "lat"
 var1!1 = "lon"
 var1&lat = lat
 var1&lon = lon
 var1@_FillValue = 1e+35

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

  b = addfile("/ptmp/notaro/0K0KFIXVEG/ha.foam1-5lpj.0k0kfixveg.anni6351034.slp.nc","r") ; open data file 2
  time2 = b->time   ; retrieve time from file b

  psl2   = b->psl    ; retrieve slp from file b
  slp2 = new( (/40,48,400/),float)   
  slp2 = psl2(lat|:,lon|:,time|:)   ; reorder slp dimensions
 slp2!0 = "lat"
 slp2!1 = "lon"
 slp2!2 = "time"
 slp2&lat = lat
 slp2&lon = lon
 slp2@time = time2
 slp2@_FillValue = 1e+35

 var2 = new((/40,48/),float) 
 var2 = dim_variance(slp2)    ; calculate variance
 copy_VarCoords(var1,var2)
 var2@_FillValue = 1e+35

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

 perdiff=new((/40,48/),float)
 perdiff=(var1-var2)/var2*100.   ; % difference in variance
 copy_VarCoords(var1,perdiff)
 perdiff=smth9(perdiff,0.5,0.25,True)

 diff=new((/40,48/),float)
 diff=(var1-var2)    ; diff in variance
 copy_VarCoords(var1,diff)

 xEqv = equiv_sample_size (slp1, 0.1,0)  ; equivalent sample size
 yEqv = equiv_sample_size (slp2, 0.1,0) 
 gwt = latGauWgt (40, "lat", "gaussian weights", "dimension_less ")
 xN   = wgt_areaave (xEqv, gwt, 1., 0)   
 yN   = wgt_areaave (yEqv, gwt, 1., 0) 

 alpha=new((/40,48/),float)
 alpha = 100.*(1. - ftest(var1,xN,var2,yN,0)) ; compute significance
 copy_VarCoords(var1,alpha)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  wks = gsn_open_wks("ps","variance_diff")          ; open a ps file
  plot=new(5,graphic)

  res          = True                             ; plot mods desired
  res@gsnDraw=False
  res@gsnFrame=False
  res@tiMainString        = "Var of Annual SLP - Run B 400Yrs 280ppmv"
  res@tiMainFontHeightF = 0.011
  res@lbPerimOn          = False
  res@mpPerimOn          = True
  res@mpOutlineOn      = True           ; turn on map outline
  res@cnLineThicknessF = 1.5                      ; thicker lines
  res@tmXTOn=False
  res@tmYROn=False
  gsn_define_colormap(wks,"wxpEnIR")  ; choose color map
  res@lbLabelAutoStride    = True
  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            = (/1.,2.,3.,4.,5.,6.,7.,8./)
  res@cnLineLabelFontHeightF = 0.009
  res@cnMonoLineDashPattern = False
  res@cnLineDashPatterns = (/0,0,0,0,0,0,0,0/)
  res@cnFillOn=True
  res@cnFillColors=(/-1,-1,-1,-1,-1,-1,-1,-1,-1/)
  res@lbLabelBarOn        = False
  res@pmLabelBarDisplayMode =  False
  res@mpLandFillColor = -1
  plot(0) = gsn_csm_contour_map_ce(wks,var2,res)

  res@tiMainString        = "Var of Annual SLP - Run A 400Yrs 280ppmv"
  plot(1) = gsn_csm_contour_map_ce(wks,var1,res)

  delete(res@cnLevels)
  delete(res@cnLineDashPatterns)
  delete(res@cnFillColors)
  res@tiMainString        = "% Diff in Var of Annual SLP (Run A-Run B)"
  res@cnLevels = (/-38.,-36.,-32.,-28.,-24.,-20.,-16.,-12.,-8.,-4.,4.,8.,12.,16.,20.,24.,28.,32.,36.,38./)
  res@cnLineDashPatterns = (/2,2,2,2,2,2,2,2,2,2,0,0,0,0,0,0,0,0,0,0/)
  res@cnFillColors=(/-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,51,51,51,51,70,70,70,70,70,70/)
  plot(2) = gsn_csm_contour_map_ce(wks,perdiff,res)

  res@tiMainString        = "Diff in Var of Annual SLP (Run A-Run B)"
  res@cnLevels = (/-.90,-.85,-.75,-.65,-.55,-.45,-.35,-.25,-.15,-.05,.05,.15,.25,.35,.45,.55,.65,.75,.85,.90/)
  res@cnLineDashPatterns = (/2,2,2,2,2,2,2,2,2,2,0,0,0,0,0,0,0,0,0,0/)
  res@cnFillColors=(/-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,51,51,51,51,70,70,70,70,70,70/)
  plot(3) = gsn_csm_contour_map_ce(wks,diff,res)

  res@cnRasterModeOn=True
  res@tiMainString        = "F-Test Alpha (Shading 90%+)"
  res@cnInfoLabelOn         = False
  res@cnLineLabelsOn       = False              ; turn off contour line labels
  res@cnLinesOn = False
  delete(res@cnLevels)
  delete(res@cnLineDashPatterns)
  delete(res@cnFillColors)
  res@cnLevels = (/90./)
  res@cnMonoLineDashPattern = False
  res@cnLineDashPatterns = (/0/)
  res@cnFillColors=(/-1,85/)
  plot(4) = gsn_csm_contour_map_ce(wks,alpha,res)

  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.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,(/2,2,1/),pres)        ; create panel plot

end