;************************************************
;plot annual surface temperature stsu
;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"

function readT(path,startyr,endyr,area,asum)
local f,Tmn,yr,Tann,T
begin

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

 do i = 1,5
  f = addfile(path+"_"+sprinti("%0.1d",i)+"/stsu_testc.nc","r")
  Tmn = f->tsu

  Tann = dim_avg_n(Tmn,1)
  T(i-1,:) = dim_sum_n(Tann*conform(Tann,area,(/1,2/)),(/1,2/))/asum
 end do

  T@units = Tmn@units

  T@_FillValue = -9999.

  delete(f)

  return T

end

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
  asum = sum(area)                                      

  startyr = 1
  endyr = 251
  startyrx = 157
  endyrx = 251

; READ stsc*.nc 
  T1 = readT("/data/noflush/zha11a/coal_1/coal_5runs/run1",startyr,endyr,area,asum)
  T2 = readT("/data/noflush/zha11a/coal_1/coal_5runs/run2",startyr,endyr,area,asum)
  T3 = readT("/data/noflush/zha11a/coal_1/coal_5runs/run3",startyr,endyr,area,asum)

  T1x = readT("/data/noflush/zha11a/coal_1/coal_5runs/run1x",startyrx,endyrx,area,asum)
  T2x = readT("/data/noflush/zha11a/coal_1/coal_5runs/run2x",startyrx,endyrx,area,asum)
  T3x = readT("/data/noflush/zha11a/coal_1/coal_5runs/run3x",startyrx,endyrx,area,asum)

  iyear = (endyr-startyr+1) - (endyrx-startyrx+1)

  T1y = T1
  T2y = T2
  T3y = T3

;  T1y = T1x@_FillValue
;  T2y = T2x@_FillValue
;  T3y = T3x@_FillValue

  T1y(:,iyear:) = T1x
  T2y(:,iyear:) = T2x
  T3y(:,iyear:) = T3x

  startyr = startyr + 1850
  endyr = endyr + 1850

  startyrx = startyrx + 2005
  endyrx = endyrx + 2005

  y = ispan(startyr,endyr,1)

; average the 5 ensembles.
  L1_avg = dim_avg_n(T1,0)
  L2_avg = dim_avg_n(T2,0)
  L3_avg = dim_avg_n(T3,0)
  L1x_avg = dim_avg_n(T1y,0)
  L2x_avg = dim_avg_n(T2y,0)
  L3x_avg = dim_avg_n(T3y,0)

; write as NETCDF
  fo = addfile("tsu_global_r5.nc","c")
  fo->L1 = T1
  fo->L2 = T2
  fo->L3 = T3
  fo->L1x = T1y
  fo->L2x = T2y
  fo->L3x = T3y
  fo->L1_avg = L1_avg
  fo->L2_avg = L2_avg
  fo->L3_avg = L3_avg
  fo->L1x_avg = L1x_avg
  fo->L2x_avg = L2x_avg
  fo->L3x_avg = L3x_avg

; add mov. average
;  T1 = runave(T1,10,0) - 273.16
;  T2 = runave(T2,10,0) - 273.16
;  T3 = runave(T3,10,0) - 273.16
;  T1y = runave(T1y,10,0) - 273.16
;  T2y = runave(T2y,10,0) - 273.16
;  T3y = runave(T3y,10,0) - 273.16

; convert from K to celcius degree
  L1_avg = L1_avg - 273.16
  L2_avg = L2_avg - 273.16
  L3_avg = L3_avg - 273.16
  L1x_avg = L1x_avg - 273.16
  L2x_avg = L2x_avg - 273.16
  L3x_avg = L3x_avg - 273.16

;---------------------------------------------------------------
; graph plotting

  wks = gsn_open_wks("pdf","tsu_global_r5")

  gsn_define_colormap(wks,"BkBlAqGrYeOrReViWh200")

  res                  = True
  res@gsnMaximize      = True
  res@tiMainString     = "Global Surface Temperature (~S~o~N~C)"  ; add title
  res@vpWidthF         = 0.8                   ; set width of plot
  res@vpHeightF        = 0.6                   ; set height of plot
  res@xyLineThicknessF = 2.5
  res@xyDashPatterns   = (/0,0,0,2,2,2/)
  res@xyMonoLineColor  = False
; res@xyLineColors     = (/"black","darkgreen","red","black","darkgreen","red"/)
  res@xyLineColors     = (/2,75,154,2,75,154/)
  res@txFontHeightF    = 0.15                  ; change title font heights
  res@trXMinF          = startyr               ; y axis min
  res@trXMaxF          = endyr                 ; y axis max
;  res@trYMinF          = 284.8                 ; y axis min
;  res@trYMaxF          = 288.4                 ; y axis max
;  res@trYMinF          = 11                    ; y axis min
;  res@trYMaxF          = 16                    ; 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   = 0.2                 ; move units right
  res@pmLegendOrthogonalPosF = -0.5                ; move units down

  res@pmLegendWidthF         = 0.12                ; Change width and
  res@pmLegendHeightF        = 0.20                ; height of legend.
  res@lgLabelFontHeightF     = 0.02                ; change font height
; res@lgOrientation          = "horizontal"
  res@lgPerimOn              = False               ; no box around
  res@xyExplicitLegendLabels = (/" C-RCP8.5", " CN-RCP8.5", " CNP-RCP8.5", " C-RCP2.6", " CN-RCP2.6", " CNP-RCP2.6"/)

  plot = gsn_csm_xy(wks,y,(/L1_avg,L2_avg,L3_avg,L1x_avg,L2x_avg,L3x_avg/),res)


end