;*************************************************
;plot annual NPP for 6 runs
;Qian Zhang @ 02/06/2010
;************************************************
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
    
  case = (/"run1","run2","run3","run4","run5","run6"/)

; 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 for 6 simulations

  nlines = 1 
  nyr    = 180
  NPP    = new((/nlines,nyr/),"float")
    
  do k = 0,nlines-1

    fils = systemfunc("ls /short/p66/yxw599/coal_1/"+case(k)+"/archive/out_cable_*.nc")
    f = addfiles(fils,"r")

    do i = 0, nyr-1

      NPPmn = f[i]->NEE
      frac  = f[i]->patchfrac

      NPPann = dim_avg_n(NPPmn,0)*86400*365*12*1.0E-6 ;gC/m2/year
      NPPmap = dim_sum_n(NPPann*frac,0)
      NPPsum = dim_sum(NPPmap*area)*1.0E-15     ; scale from g to GT
;      printVarSummary (NPPmap)

      NPP(k,i) = dim_sum(NPPsum)
;      print(NPP(i))

    end do

  end do

  y = ispan(1,nyr,1)
  
;------------------------------------------------------------
; graph plotting

  wks      = gsn_open_wks("pdf","NEE-spinup-C-20111208")

  res                  = True
  res@gsnMaximize      = True
;  res@tiMainString     = ""  ; 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@xyDashPatterns   = (/0/)
  res@xyLineColors     = (/"red"/) 
;  res@txFontHeightF    = 0.15                 ; change title font heights
;  res@trXMinF          = 1870                 ; y axis min
;  res@trXMaxF          = 2010                 ; y axis max
  res@tiXAxisString    = "year"
  res@tiYAxisString    = "NEE Gt C a~S~-1~N~"

; attach legend
  res@pmLegendDisplayMode    = "Always"

  res@pmLegendSide           = "Top"               ; Change location of 
  res@pmLegendParallelPosF   = .25                 ; move units right
  res@pmLegendOrthogonalPosF = -0.3                ; 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"/)
  res@xyExplicitLegendLabels = (/" SPIN-C"/)

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

end