;*************************************************
;plot annual surface temperature on land.
;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,1)                             ; mask out the ocean
  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/out_cable_"+sprinti("%0.5d",i)+".nc","r")
    Tmn = f->Tair
    frac  = f->patchfrac

    Tann = dim_avg_n(Tmn,0)
    Tmap = dim_sum_n(Tann*frac,0)
    Tsum = dim_sum(Tmap*area)
;     printVarSummary (Tmap)

    T(i-startyr) = dim_sum(Tsum)/asum

  end do

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

  wks      = gsn_open_wks("pdf","tair_land_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.4                ; move units down

  res@pmLegendWidthF         = 0.12                 ; Change width and
  res@pmLegendHeightF        = 0.14                 ; height of legend.
  res@lgLabelFontHeightF     = .015                 ; change font height
;  res@lgOrientation         = "horizontal"
  res@lgPerimOn              = False               ; no box around
;  res@xyExplicitLegendLabels = (/" NOFERT_C"," NOFERT_CN"," NOFERT_CNP"," FERT_C"," FERT_CN"," FERT_CNP"/)

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

end