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/pimodTC1/ha.foam1-5lpj.pimodTC1.anni10351134.nc", "r")   ; Ensemble 1  
 yrs = fspan(1.,100.,100)
 lat   = a->lat                                  ; get lat
 lon   = a->lon                                  ; get lon
 ta1   = a->T(:,17,:,:)-273.15                         ; get data
 ta1!0 = "time"
 ta1!1 = "lat"
 ta1!2 = "lon"
 ta1&time=yrs
 ta1&lat=lat
 ta1&lon=lon
 sfcta1 = new( (/40,48,100/),float)
 sfcta1 = ta1(lat|:,lon|:,time|:)

 b     = addfile ("/ptmp/notaro/pimodTC2/ha.foam1-5lpj.pimodTC2.anni10171116.nc", "r")    ; Ensemble 2
 ta2   = b->T(:,17,:,:)-273.15                         ; get data
 copy_VarCoords(ta1,ta2)
 sfcta2 = new( (/40,48,100/),float)
 sfcta2 = ta2(lat|:,lon|:,time|:)

 c     = addfile ("/ptmp/notaro/pimodTC3/ha.foam1-5lpj.pimodTC3.anni09991098.nc", "r")    ; Ensemble 3
 time3  = c->time
 ta3   = c->T(:,17,:,:)-273.15                         ; get data
 copy_VarCoords(ta1,ta3)
 sfcta3 = new( (/40,48,100/),float)
 sfcta3 = ta3(lat|:,lon|:,time|:)

 tens=new((/40,48,100/),float)  ; Ensemble av
 tens=(sfcta1+sfcta2+sfcta3)/3.
 copy_VarCoords(sfcta1,tens)

 times = ispan(1,100,1)
 yDtrend = dtrend_msg(times,tens,False,True)     ; Calculate the trend
 slope2d = onedtond(yDtrend@slope, (/40,48/) )   ; units of C/yr
 slope2d = slope2d*10.                           ; units of C/decade

 slope2d!0 = "lat"
 slope2d!1 = "lon"
 slope2d&lat = lat
 slope2d&lon = lon
 print(max(slope2d))
 print(min(slope2d))

 wks   = gsn_open_wks ("ps", "trend_ensemble" )            ; open workstation
 gsn_define_colormap(wks,"BlAqGrYeOrRevi200")  ; choose color map
 res          = True                             ; plot mods desired
 res@tiMainString        = "FOAM-LPJ PIMODTC1-3 Trend in Sfc Air Temp (C/decade) (100Yrs)"  ; main title
 res@lbPerimOn          = False
 res@mpPerimOn          = True
  res@mpOutlineOn      = True           ; turn on map outline
  res@cnLineThicknessF = 1.5                      ; thicker lines
  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 = "ManualLevels"     ; set manual contour levels
  res@cnMinLevelValF       = -0.10               ; set min contour level
  res@cnMaxLevelValF       = 0.40                 ; set max contour level
  res@cnLevelSpacingF      = 0.05               ; set contour spacing
  res@cnLineLabelFontHeightF = 0.009
  res@gsnDraw          = False   ; don't draw yet
  res@gsnFrame         = False   ; don't advance frame yet
  res@tmXTOn = False
  res@tmYROn = False
  plot = gsn_csm_contour_map_ce(wks,slope2d,res)
  plot = ColorNegDashZeroPosContour(plot,"blue","blue","red")
  draw(wks)
  frame(wks)

end