;*************************************************
;plot annual sea surface temperature 
;Qian Zhang @ 11/11/2011
;************************************************
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/contributed.ncl"

begin

; READ land areas
  fm = addfile("csiro_landsea.nc","r")                 ; read the mask
  f0 = addfile("csiro_agcm_area.nc","r")               ; read grid area
  oro  = fm->mask
  area = f0->area                                      ; unit in m^2
  area = mask (area,oro,0)                             ; mask out the land
  asum = sum(area)                                      

; READ out_cable*.nc 

  startyr = 1
  endyr   = 220
  T    = new(endyr-startyr+1,"float")

  do i = startyr, endyr

    f = addfile("/short/p66/yxw599/coal_1/run1/archive/annual_"+sprinti("%0.5d",i)+".nc","r")
    Tmn = f->obgc01(0,0,:,:)
    Tmn = mask (Tmn,oro,0)                           ; mask out land grid
;     printVarSummary (Tmn)

    T(i-startyr) = dim_sum(Tmn*area)/asum                   

  end do

  y = ispan(startyr,endyr,1)
  
;---------------------------------------------------------------
; graph plotting

  wks      = gsn_open_wks("x11","sst_coal_20111212")

  res                  = True
  res@gsnMaximize      = True
  res@tiMainString     = "Surface Air Temperature"  ; add title
  res@vpWidthF         = 0.8                  ; set width of plot
  res@vpHeightF        = 0.6                  ; set height of plot
  res@xyLineThicknessF = 3
;  res@xyDashPatterns   = (/0,0,0,1,1,1/)
;  res@xyLineColors     = (/"black","darkgreen","red","black","darkgreen","red"/) 
  res@txFontHeightF    = 0.15               ; change title font heights
  res@trXMinF          = startyr                 ; y axis min
  res@trXMaxF          = endyr                 ; y axis max
  res@tiXAxisString    = "year"
  res@tiYAxisString    = "K"

; Manually create and attach legend
  res@pmLegendDisplayMode    = "Always"

  res@pmLegendSide           = "Top"               ; Change location of
  res@pmLegendParallelPosF   = .2                  ; move units right
  res@pmLegendOrthogonalPosF = -0.25                ; move units down

  res@pmLegendWidthF         = 0.12                 ; Change width and
  res@pmLegendHeightF        = 0.14                 ; height of legend.
  res@lgLabelFontHeightF     = .02                 ; change font height
;  res@lgOrientation         = "horizontal"
  res@lgPerimOn              = False               ; no box around
  res@xyExplicitLegendLabels = (/" sst"/)

  plot = gsn_csm_xy(wks,y,T,res)

end