;waccmx_1944_plot_yearly_pres.pro
;
;Plot global monthly and yearly mean WACCM-X 1944-1970 CO2 and Temperature at all pressure levels from the surface to the thermosphere using data from 
;h0 monthly files.  Global mean is calculated using an area weighted algorithm with the area calculated using the spherical cap method.  Then a
;analysis is performed as well as a fit to residuals to get the trend  
;
;  Restore the .sav files with monthly and yearly data
;
print, 'Restoring yearly global mean data '

;restore, file='/glade/scratch/joemci/TrendSavFiles/CO2_T_Yearly_GlobalMean_Pres_1944_1962.sav'
restore, file='TrendSavFiles/CO2_T_Yearly_GlobalMean_Pres_1944_1962.sav'

;  Plot monthly and yearly global mean data
;
!P.THICK = 3
!x.thick = 3
!y.thick = 3

SET_PLOT, 'PS'

psfile = 'yearly_1944-1970_CO2_T_AreaGlobalMean_WACCMX_Pres.ps'

DEVICE, /COLOR, FILENAME = psfile
loadct, 39, ncolors=10

YearsAll = INDGEN(nYears)

CYears = [' ','1949',' ','1959',' ','1969',' ','1990',' ','2000',' ','2010']

FOR iPres = 0,nPres-1 DO BEGIN

 yTitle = 'CO2 (mol/mol)'
 xTitle = 'Year'
 
 cPresTemp = STRING(pres[iPres],FORMAT='(E8.1)')
 cPres = STRTRIM(cPresTemp,2)

 Title = 'CO2 1944-1970 Yearly Global Mean '+cPres+' hPa'
 plot, YearsAll, CO2YearlyGlobalMean[0:nYears-1,iPres],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1,charsize=1.3,xmargin=[12,4],xticks=11,xtickname=CYears

 yTitle = 'T (K)'
 
 Title = 'T 1944-1970 Yearly Global Mean '+cPres+' hPa'
 plot, YearsAll, TYearlyGlobalMean[0:nYears-1,iPres],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1,charsize=1.3,xticks=11,xtickname=CYears

 yTitle = 'NO (mol/mol)'
 
 Title = 'NO 1944-1970 Yearly Global Mean '+cPres+' hPa'
 plot, YearsAll, NOYearlyGlobalMean[0:nYears-1,iPres],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1,charsize=1.3,xmargin=[12,4],xticks=11,xtickname=CYears

ENDFOR ; Pressures

yTitle = 'T (K)'
 
Title = 'T 1944-1970 Yearly Global Mean 2 Meter Reference Height '
plot, YearsAll, TREFHTYearlyGlobalMean[0:nYears-1],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1,charsize=1.3,xticks=11,xtickname=CYears

device,/close

CYears = ['1955','1960','1965','1970','1975','1980','1985','1990','1995','2000','2005','2010']

psfile = 'year_pres_CO2_T_1944-1970_AreaGlobalMean_WACCMX.ps'

DEVICE, /COLOR, FILENAME = psfile
loadct, 39

yTitle = 'Pressure (hPa)'
xTitle = 'Year'

levels = 21
clevels = findgen(levels) * 2.E-05
  
Title = ' WACCM-X CO2 1944-1970 Area Weighted Global Mean '
CONTOUR, CO2YearlyGlobalMean[0:nYears-1,18:80],YearsAll,pres[18:80],/fill,levels=clevels,/ylog,yrange=[1100.,5.E-05],ystyle=1,xticks=11,xtickname=CYears,Title=Title,xtitle=xTitle,ytitle=yTitle
;CONTOUR, CO2YearlyGlobalMean[0:nYears-1,11:80], YearsAll, pres[11:80], nlev=24,/fill,c_colors=BYTSCL(indgen(24)), /ylog, yrange=[1000.,10E-06]
;,xticks=51,xtickname=CYears
;,xrange=[0,11],yrange=[-90
;.0,90.0], xstyle=1,ystyle=1,Title=Title,xtitle=xTitle,ytitle=yTitle,ymargin=[7,3],xticks=11,xtickname=CMonths
;Title = ' WACCM-X CO2 1944-1970 Area Weighted Global Mean '
CONTOUR, CO2YearlyGlobalMean[0:nYears-1,18:80], YearsAll, pres[18:80], levels=clevels,/follow, /ylog, yrange=[1100.,5.E-05],/overplot,ystyle=1
;CONTOUR, CO2YearlyGlobalMean[0:nYears-1,11:80], YearsAll, pres[11:80], nlev=12,/follow, /ylog, yrange=[1000.,10E-06],/overplot
;,xticks=51,xtickn=Months,xtitle='Year'
;,xrange=[0,11], yrange=[-90.0,90.0],xstyle=1,ystyle=
;1,/overplot,c_charsize=1.2, xtickname = CMonths
;  axis, xaxis=1,/save,xticks=12,xtickn=Months,xtitle='Month'

clevels = [1.E-12,5.E-12,1.E-11,5.E-11,1.E-10,5.E-10,1.E-09,5.E-09,1.E-08,5.E-08,1.E-07,5.E-07,1.E-06,5.E-06,1.E-05,5.E-05,1.E-04]

;Title = ' WACCM-X CO2 1944-1970 Area Weighted Global Mean '
CONTOUR, CO2YearlyGlobalMean[0:nYears-1,0:20], YearsAll, pres[0:20],/fill,levels=clevels, /ylog, yrange=[5.E-05,3.E-09],ystyle=1,xticks=11,xtickname=CYears,Title=Title,xtitle=xTitle,ytitle=yTitle
;CONTOUR, CO2YearlyGlobalMean[0:nYears-1,0:10], YearsAll, pres[0:10], nlev=24,/fill,c_colors=BYTSCL(indgen(24)), /ylog, yrange=[10E-06,10E-10]
;,xticks=51,xtickname=CYears
;,xrange=[0,11],yrange=[-90
;.0,90.0], xstyle=1,ystyle=1,Title=Title,xtitle=xTitle,ytitle=yTitle,ymargin=[7,3],xticks=11,xtickname=CMonths
;Title = ' WACCM-X CO2 1944-1970 Area Weighted Global Mean '
CONTOUR, CO2YearlyGlobalMean[0:nYears-1,0:20], YearsAll, pres[0:20], levels=clevels,/follow, /ylog, yrange=[5.E-05,3.E-09],/overplot,ystyle=1
;CONTOUR, CO2YearlyGlobalMean[0:nYears-1,0:10], YearsAll, pres[0:10], nlev=12,/follow, /ylog, yrange=[10E-06,10E-10],/overplot
;,xticks=51,xtickn=Months,xtitle='Year'
;,xrange=[0,11], yrange=[-90.0,90.0],xstyle=1,ystyle=
;1,/overplot,c_charsize=1.2, xtickname = CMonths
;  axis, xaxis=1,/save,xticks=12,xtickn=Months,xtitle='Month'

levels = 21
clevels = findgen(levels) * 16.

Title = ' WACCM-X T 1944-1970 Area Weighted Global Mean '
CONTOUR, TYearlyGlobalMean[0:nYears-1,18:80], YearsAll, pres[18:80],/fill,levels=clevels, /ylog, yrange=[1100.,5.E-05],ystyle=1,xticks=11,xtickname=CYears,Title=Title,xtitle=xTitle,ytitle=yTitle
;CONTOUR, TYearlyGlobalMean[0:nYears-1,15:80], YearsAll, pres[15:80], nlev=24,/fill,c_colors=BYTSCL(indgen(24)), /ylog, yrange=[1100.,5.E-05]
;,xticks=51,xtickname=CYears
;,xrange=[0,11],yrange=[-90
;.0,90.0], xstyle=1,ystyle=1,Title=Title,xtitle=xTitle,ytitle=yTitle,ymargin=[7,3],xticks=11,xtickname=CMonths
;Title = ' WACCM-X T 1944-1970 Area Weighted Global Mean '
CONTOUR, TYearlyGlobalMean[0:nYears-1,18:80], YearsAll, pres[18:80], levels=clevels,/follow, /ylog, yrange=[1100.,5.E-05],/overplot,ystyle=1
;CONTOUR, TYearlyGlobalMean[0:nYears-1,15:80], YearsAll, pres[15:80], nlev=12,/follow, /ylog, yrange=[1100.,5.E-05],/overplot
;,xticks=51,xtickn=Months,xtitle='Year'
;,xrange=[0,11], yrange=[-90.0,90.0],xstyle=1,ystyle=1,/overplot,c_charsize=1.2, xtickname = CMonths
;  axis, xaxis=1,/save,xticks=12,xtickn=Months,xtitle='Month'

levels = 16
clevels = findgen(levels) * 75.

;Title = ' WACCM-X T 1944-1970 Area Weighted Global Mean '
CONTOUR, TYearlyGlobalMean[0:nYears-1,0:20], YearsAll, pres[0:20],/fill,levels=clevels, /ylog, yrange=[5.E-05,3.E-09],ystyle=1,xticks=11,xtickname=CYears,Title=Title,xtitle=xTitle,ytitle=yTitle
;CONTOUR, TYearlyGlobalMean[0:nYears-1,0:15], YearsAll, pres[0:15], nlev=24,/fill,c_colors=BYTSCL(indgen(24)), /ylog, yrange=[10E-04,10E-10]
;,xticks=51,xtickname=CYears
;,xrange=[0,11],yrange=[-90
;.0,90.0], xstyle=1,ystyle=1,Title=Title,xtitle=xTitle,ytitle=yTitle,ymargin=[7,3],xticks=11,xtickname=CMonths
;Title = ' WACCM-X T 1944-1970 Area Weighted Global Mean '
CONTOUR, TYearlyGlobalMean[0:nYears-1,0:20], YearsAll, pres[0:20], levels=clevels,/follow, /ylog, yrange=[5.E-05,3.E-09],/overplot,ystyle=1
;CONTOUR, TYearlyGlobalMean[0:nYears-1,0:15], YearsAll, pres[0:15], nlev=12,/follow, /ylog, yrange=[10E-04,10E-10],/overplot
;,xticks=51,xtickn=Months,xtitle='Year'
;,xrange=[0,11], yrange=[-90.0,90.0],xstyle=1,ystyle=
;1,/overplot,c_charsize=1.2, xtickname = CMonths
;  axis, xaxis=1,/save,xticks=12,xtickn=Months,xtitle='Month'

device,/close
;
;  Get monthly and yearly mean f107 and ap for regression
;
waccmx_solar_f107_ap_monthly_yearly_1944, f107AllMonthly, apAllMonthly, f107YearlyMean, apYearlyMean
;
;  Now do regression of yearly mean temperature with f10.7, ap, and f10.7*f10.7
;
f107YearlyMeanSquared = f107YearlyMean * f107YearlyMean

f107apYearlyMean_1944_1970 = FLTARR(3,nYears)
f107apYearlyMean_1944_1970[0,*] = f107YearlyMean[0:nYears-1]
f107apYearlyMean_1944_1970[1,*] = apYearlyMean[0:nYears-1]
f107apYearlyMean_1944_1970[2,*] = f107YearlyMeanSquared[0:nYears-1]

;f107apYearlyMean_1970_2000 = FLTARR(3,nYears-25)
;f107apYearlyMean_1970_2000[0,*] = f107YearlyMean[15:nYears-11]
;f107apYearlyMean_1970_2000[1,*] = apYearlyMean[15:nYears-11]
;f107apYearlyMean_1970_2000[2,*] = f107YearlyMeanSquared[15:nYears-11]

residualsCO2YearlyGMean_1944_1970 = FLTARR(nYears,nPres)
residualsTYearlyGMean_1944_1970 = FLTARR(nYears,nPres)
residualsDensYearlyGMean_1944_1970 = FLTARR(nYears,nPres)
residualsNOYearlyGMean_1944_1970 = FLTARR(nYears,nPres)
fitCO2YearlyGMean_1944_1970 = FLTARR(nYears,nPres)
fitTYearlyGMean_1944_1970 = FLTARR(nYears,nPres)
fitDensYearlyGMean_1944_1970 = FLTARR(nYears,nPres)
fitNOYearlyGMean_1944_1970 = FLTARR(nYears,nPres)

;residualsCO2YearlyGMean_1970_2000 = FLTARR(nYears-25,nPres)
;residualsTYearlyGMean_1970_2000 = FLTARR(nYears-25,nPres)
;residualsDensYearlyGMean_1970_2000 = FLTARR(nYears-25,nPres)
;residualsNOYearlyGMean_1970_2000 = FLTARR(nYears-25,nPres)
;fitCO2YearlyGMean_1970_2000 = FLTARR(nYears-25,nPres)
;fitTYearlyGMean_1970_2000 = FLTARR(nYears-25,nPres)
;fitDensYearlyGMean_1970_2000 = FLTARR(nYears-25,nPres)
;fitNOYearlyGMean_1970_2000 = FLTARR(nYears-25,nPres)

mean_CO2_1955_2010 = FLTARR(nPres)
;mean_CO2_1970_2000 = FLTARR(nPres)

psfile = 'yearly_1944-1970_CO2_T_AreaGlobalMean_WACCMX_Pres_Trends.ps'

CYears = [' ','1960',' ','1970',' ','1980',' ','1990',' ','2000',' ','2010']
;CYears1970 = ['1970','1975','1980','1985','1990','1995','2000']

DEVICE, /COLOR, FILENAME = psfile
loadct, 39, ncolors=10

FOR iPres = 0,nPres-1 DO BEGIN

  CO2YearlyGMean_1944_1970 = CO2YearlyGlobalMean[0:nYears-1,iPres]
  regressCoeffs_1944_1970 = REGRESS(f107apYearlyMean_1944_1970, CO2YearlyGMean_1944_1970, YFIT=CO2YearlyGMeanFit_1944_1970)
  residualsCO2YearlyGMean_1944_1970[*,iPres] = CO2YearlyGMean_1944_1970 - CO2YearlyGMeanFit_1944_1970
  
  mean_CO2_1955_2010[iPres] = MEAN(CO2YearlyGMean_1944_1970)

;  residualsCO2YearlyGMean_1970_2000[*,iPres] = residualsCO2YearlyGMean_1944_1970[15:nYears-11,iPres]

;  mean_CO2_1970_2000[iPres] = MEAN(CO2YearlyGMean_1944_1970[15:nYears-11])
  ;
  ;  Get fit to residuals
  ;
  Years_1944_1970 = INDGEN(nYears)
  fitCoeffs = POLY_FIT(Years_1944_1970,residualsCO2YearlyGMean_1944_1970[*,iPres],1,YFIT=fitCO2YGMean_1944_1970)
  fitCO2YearlyGMean_1944_1970[*,iPres] = fitCO2YGMean_1944_1970

  perCO2Increase = (fitCO2YGMean_1944_1970[nYears-1] - fitCO2YGMean_1944_1970[0]) / mean_CO2_1955_2010[iPres] * 100.
  cperCO2Increase = STRTRIM(perCO2Increase,2)
  perCO2IncreaseDec = perCO2Increase / 5.6
  cperCO2IncreaseDec = STRTRIM(perCO2IncreaseDec,2)
 
;  Years_1970_2000 = INDGEN(nYears-25)
;  fitCoeffs_1970_2000 = POLY_FIT(Years_1970_2000,residualsCO2YearlyGMean_1970_2000[*,iPres],1,YFIT=fitCO2YGMean_1970_2000) 
;  fitCO2YearlyGMean_1970_2000[*,iPres] = fitCO2YGMean_1970_2000

  yTitle = 'CO2 (mol/mol)'
  xTitle = 'Year'
 
  cPresTemp = STRING(pres[iPres],FORMAT='(E8.1)')
  cPres = STRTRIM(cPresTemp,2)

  Title = 'CO2 1944-1970 Yearly Global Mean '+cPres+' hPa'
  plot, Years_1944_1970,residualsCO2YearlyGMean_1944_1970[*,iPres],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1,charsize=1.3,xmargin=[12,4],xticks=11,xtickname=CYears
  oplot, fitCO2YearlyGMean_1944_1970[*,iPres],color=230, thick=6
  XYOUTS,0.21,0.88,'CO2 % Change = '+cperCO2Increase+'%',/normal,color=230,charsize=1.1, charthick=3
  XYOUTS,0.21,0.84, '('+cperCO2IncreaseDec+' %/decade)',/normal,color=230,charsize=1.1, charthick=3

;  perCO2Increase = (fitCO2YGMean_1970_2000[nYears-26] - fitCO2YGMean_1970_2000[0]) / mean_CO2_1970_2000[iPres] * 100.
;  cperCO2Increase = STRTRIM(perCO2Increase,2)
;  perCO2IncreaseDec = perCO2Increase / 3.1
;  cperCO2IncreaseDec = STRTRIM(perCO2IncreaseDec,2)
 
;  Title = 'CO2 1970-2000 Yearly Global Mean '+cPres+' hPa'
;  plot, Years_1970_2000,residualsCO2YearlyGMean_1970_2000[*,iPres],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1,charsize=1.3,xmargin=[12,4],xticks=6,xtickname=CYears1970
;  oplot, fitCO2YearlyGMean_1970_2000[*,iPres],color=230, thick=6
;  XYOUTS,0.21,0.88,'CO2 % Change = '+cperCO2Increase+'%',/normal,color=230,charsize=1.1, charthick=3
;  XYOUTS,0.21,0.84, '('+cperCO2IncreaseDec+' %/decade)',/normal,color=230,charsize=1.1, charthick=3
  ;
  ;  Do temperature
  ;
  TYearlyGMean_1944_1970 = TYearlyGlobalMean[0:nYears-1,iPres]
  regressCoeffs_1944_1970 = REGRESS(f107apYearlyMean_1944_1970, TYearlyGMean_1944_1970, YFIT=TYearlyGMeanFit_1944_1970)
  residualsTYearlyGMean_1944_1970[*,iPres] = TYearlyGMean_1944_1970 - TYearlyGMeanFit_1944_1970

;  residualsTYearlyGMean_1970_2000[*,iPres] = residualsTYearlyGMean_1944_1970[15:nYears-11,iPres]
;
;  Get fit to residuals
;
  fitCoeffs = POLY_FIT(Years_1944_1970,residualsTYearlyGMean_1944_1970[*,iPres],1,YFIT=fitTYGMean_1944_1970) 
  fitTYearlyGMean_1944_1970[*,iPres] = fitTYGMean_1944_1970

  trendFit = fitTYGMean_1944_1970[nYears-1] - fitTYGMean_1944_1970[0]
  ctrendFit = STRTRIM(trendFit,2)
  trendFitPerDecade = trendFit / 5.6
  ctrendFitPerDecade = STRTRIM(trendFitPerDecade,2)

;  fitCoeffs = POLY_FIT(Years_1970_2000,residualsTYearlyGMean_1970_2000[*,iPres],1,YFIT=fitTYGMean_1970_2000) 
;  fitTYearlyGMean_1970_2000[*,iPres] = fitTYGMean_1970_2000

  yTitle = 'T (K)'

  xTitle = 'Year'

  Title = 'T 1944-1970 Yearly Global Mean '+cPres+' hPa'
  plot, Years_1944_1970,residualsTYearlyGMean_1944_1970[*,iPres],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1,charsize=1.3,xticks=11,xtickname=CYears
  oplot, fitTYearlyGMean_1944_1970[*,iPres],color=230, thick=6
  XYOUTS,0.21,0.25,'T Change = '+ctrendFit+' K',/normal,color=230,charsize=1.1, charthick=3
  XYOUTS,0.21,0.21, '('+ctrendFitPerDecade+' K/decade)',/normal,color=230,charsize=1.1, charthick=3

;  trendFit = fitTYGMean_1970_2000[nYears-26] - fitTYGMean_1970_2000[0]
;  ctrendFit = STRTRIM(trendFit,2)
;  trendFitPerDecade = trendFit / 3.1
;  ctrendFitPerDecade = STRTRIM(trendFitPerDecade,2)

;  Title = 'T 1970-2000 Yearly Global Mean '+cPres+' hPa'
;  plot, Years_1970_2000,residualsTYearlyGMean_1970_2000[*,iPres],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1,charsize=1.3,xticks=6,xtickname=CYears1970
;  oplot, fitTYearlyGMean_1970_2000[*,iPres],color=230, thick=6
;  XYOUTS,0.21,0.25,'T Change = '+ctrendFit+' K',/normal,color=230,charsize=1.1, charthick=3
;  XYOUTS,0.21,0.21, '('+ctrendFitPerDecade+' K/decade)',/normal,color=230,charsize=1.1, charthick=3
  ;
  ; Do NO
  ;
  NOYearlyGMean_1944_1970 = NOYearlyGlobalMean[0:nYears-1,iPres]
  regressCoeffs_1944_1970 = REGRESS(f107apYearlyMean_1944_1970, NOYearlyGMean_1944_1970, YFIT=NOYearlyGMeanFit_1944_1970)
  residualsNOYearlyGMean_1944_1970[*,iPres] = NOYearlyGMean_1944_1970 - NOYearlyGMeanFit_1944_1970

;  residualsNOYearlyGMean_1970_2000[*,iPres] = residualsNOYearlyGMean_1944_1970[15:nYears-11,iPres]
  ;
  ;  Get fit to residuals
  ;
  Years_1944_1970 = INDGEN(nYears)
  fitCoeffs = POLY_FIT(Years_1944_1970,residualsNOYearlyGMean_1944_1970[*,iPres],1,YFIT=fitNOYGMean_1944_1970)
  fitNOYearlyGMean_1944_1970[*,iPres] = fitNOYGMean_1944_1970
 
;  Years_1970_2000 = INDGEN(nYears-25)
;  fitCoeffs_1970_2000 = POLY_FIT(Years_1970_2000,residualsNOYearlyGMean_1970_2000[*,iPres],1,YFIT=fitNOYGMean_1970_2000) 
;  fitNOYearlyGMean_1970_2000[*,iPres] = fitNOYGMean_1970_2000

  yTitle = 'NO (mol/mol)'
  xTitle = 'Year'
 
  cPresTemp = STRING(pres[iPres],FORMAT='(E8.1)')
  cPres = STRTRIM(cPresTemp,2)

  Title = 'NO 1944-1970 Yearly Global Mean Residuals '+cPres+' hPa'
  plot, Years_1944_1970,residualsNOYearlyGMean_1944_1970[*,iPres],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1,charsize=1.3,xmargin=[12,4],xticks=11,xtickname=CYears
  oplot, fitNOYearlyGMean_1944_1970[*,iPres],color=230, thick=6

;  xTitle = 'Year'
;
;  Title = 'NO 1970-2000 Yearly Global Mean Residuals '+cPres+' hPa'
;  plot, Years_1970_2000,residualsNOYearlyGMean_1970_2000[*,iPres],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1,charsize=1.3,xmargin=[12,4],xticks=6,xtickname=CYears1970
;  oplot, fitNOYearlyGMean_1970_2000[*,iPres],color=230, thick=6

ENDFOR

  fitCoeffs = POLY_FIT(Years_1944_1970,TREFHTYearlyGlobalMean[*],1,YFIT=fitTREFHTYGMean_1944_1970)

  yTitle = 'TREFHT 2 Meter Air Temperature (K)'
  xTitle = 'Year'

  Title = 'TREFHT 1944-1970 Yearly Global Mean 2 Meter'
  plot, Years_1944_1970,TREFHTYearlyGlobalMean[*],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1,charsize=1.3,xmargin=[12,4],xticks=11,xtickname=CYears
  oplot, fitTREFHTYGMean_1944_1970[*],color=230, thick=6

;  fitCoeffs_1970_2000 = POLY_FIT(Years_1970_2000,TREFHTYearlyGlobalMean[15:nYears-11],1,YFIT=fitTREFHTYGMean_1970_2000) 
;
;  Title = 'TREFHT 1970-2000 Yearly Global Mean 2 Meter'
;  plot, Years_1970_2000,TREFHTYearlyGlobalMean[15:nYears-11],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1,charsize=1.3,xmargin=[12,4],xticks=6,xtickname=CYears1970
;  oplot, fitTREFHTYGMean_1970_2000[*],color=230, thick=6

device,/close
;
;  Plot year/altitude trend cross sections
;
psfile = 'year_pres_CO2_T_1944-1970_AreaGlobalMean_WACCMX_Trends.ps'

DEVICE, /COLOR, FILENAME = psfile
loadct, 39

CYears = ['1955','1960','1965','1970','1975','1980','1985','1990','1995','2000','2005','2010']

yTitle = 'Pressure (hPa)'
xTitle = 'Year'

levels = 21
clevels = findgen(levels)  * 3.5E-06 - 3.E-05
;clevels = findgen(levels) * 2.E-05
  
Title = ' WACCM-X CO2 1944-1970 Global Mean'
CONTOUR, residualsCO2YearlyGMean_1944_1970[0:nYears-1,36:62],YearsAll,pres[36:62],/fill,levels=clevels,yrange=[1000.0,120.],ystyle=1,xticks=11,$
                                                                                        xtickname=CYears,Title=Title,xtitle=xTitle,ytitle=yTitle,/ylog
CONTOUR, residualsCO2YearlyGMean_1944_1970[0:nYears-1,36:62],YearsAll,pres[36:62],levels=clevels,/follow, yrange=[1000.,120.],/overplot,ystyle=1

;clevels = [1.E-12,5.E-12,1.E-11,5.E-11,1.E-10,5.E-10,1.E-09,5.E-09,1.E-08,5.E-08,1.E-07,5.E-07,1.E-06,5.E-06,1.E-05,5.E-05,1.E-04]
clevels = findgen(levels)  * 5.5E-07 - 6.E-06

Title = ' WACCM-X CO2 1944-1970 Global Mean'
CONTOUR, residualsCO2YearlyGMean_1944_1970[0:nYears-1,0:38],YearsAll,pres[0:38],/fill,levels=clevels, yrange=[1000.,310.],ystyle=1,xticks=11,$
                                                                                     xtickname=CYears,Title=Title,xtitle=xTitle,ytitle=yTitle,/ylog
CONTOUR, residualsCO2YearlyGMean_1944_1970[0:nYears-1,0:38],YearsAll,pres[0:38],levels=clevels,/follow,yrange=[120.,310.],/overplot,ystyle=1

levels = 10
clevels = findgen(levels) / 1.25 - 4.

Title = ' WACCM-X T 1944-1970 Global Mean'
CONTOUR, residualsTYearlyGMean_1944_1970[0:nYears-1,36:62],YearsAll,pres[36:62],/fill,levels=clevels, yrange=[1000.,120.],ystyle=1,xticks=11,$
                                                                                    xtickname=CYears,Title=Title,xtitle=xTitle,ytitle=yTitle,/ylog
CONTOUR, residualsTYearlyGMean_1944_1970[0:nYears-1,36:62],YearsAll,pres[36:62],levels=clevels,/follow, yrange=[1000.,120],/overplot,ystyle=1

levels = 10
clevels = findgen(levels) * 14. - 68.

Title = ' WACCM-X T 1944-1970 Global Mean'
CONTOUR, residualsTYearlyGMean_1944_1970[0:nYears-1,0:38],YearsAll,pres[0:38],/fill,levels=clevels, yrange=[1000.,310.],ystyle=1,xticks=11,$
                                                                                   xtickname=CYears,Title=Title,xtitle=xTitle,ytitle=yTitle,/ylog
CONTOUR, residualsTYearlyGMean_1944_1970[0:nYears-1,0:38],YearsAll,pres[0:38],levels=clevels,/follow, yrange=[1000.,310.],/overplot,ystyle=1

levels = 21
clevels = findgen(levels)  * 3.5E-06 - 3.E-05
;clevels = findgen(levels) * 2.E-05
  
Title = ' WACCM-X NO 1944-1970 Global Mean'
CONTOUR, residualsNOYearlyGMean_1944_1970[0:nYears-1,36:62],YearsAll,pres[36:62],/fill,levels=clevels,yrange=[1000.0,120.],ystyle=1,xticks=11,$
                                                                                        xtickname=CYears,Title=Title,xtitle=xTitle,ytitle=yTitle,/ylog
CONTOUR, residualsNOYearlyGMean_1944_1970[0:nYears-1,36:62],YearsAll,pres[36:62],levels=clevels,/follow, yrange=[1000.,120.],/overplot,ystyle=1

;clevels = [1.E-12,5.E-12,1.E-11,5.E-11,1.E-10,5.E-10,1.E-09,5.E-09,1.E-08,5.E-08,1.E-07,5.E-07,1.E-06,5.E-06,1.E-05,5.E-05,1.E-04]
clevels = findgen(levels)  * 5.5E-07 - 6.E-06

Title = ' WACCM-X NO 1944-1970 Global Mean'
CONTOUR, residualsNOYearlyGMean_1944_1970[0:nYears-1,0:38],YearsAll,pres[0:38],/fill,levels=clevels, yrange=[1000.,310.],ystyle=1,xticks=11,$
                                                                                     xtickname=CYears,Title=Title,xtitle=xTitle,ytitle=yTitle,/ylog
CONTOUR, residualsNOYearlyGMean_1944_1970[0:nYears-1,0:38],YearsAll,pres[0:38],levels=clevels,/follow,yrange=[120.,310.],/overplot,ystyle=1

device,/close
;
;  Plot trend altitude vertical profiles for global mean
;
psfile = 'CO2_T_pres_1944_AreaGlobalMean_WACCMX_Trends.ps'

DEVICE, /COLOR, FILENAME = psfile
loadct, 39

yTitle = 'Pressure (hPa)'

xTitle = 'CO2 Residual Trend (mol/mol)'
;
diffsfitCO2YearlyGMean_1944_1970 = fitCO2YearlyGMean_1944_1970[nYears-1,*] - fitCO2YearlyGMean_1944_1970[0,*]
;
perdiffsfitCO2YearlyGMean_1944_1970 = FLTARR(nPres)

geomAltMean = FLTARR(nPres)

FOR iPres=0,nPres-1 DO BEGIN

  perdiffsfitCO2YearlyGMean_1944_1970[iPres] = diffsfitCO2YearlyGMean_1944_1970[iPres] / $
                                         ((CO2YearlyGlobalMean[nYears-1,iPres] + CO2YearlyGlobalMean[0,iPres]) / 2.) * 100.
 
  geomAltMean[iPres] = MEAN(geomAltYearlyGlobalMean[0:nYears-1,iPres])

ENDFOR

xTitle = 'CO2 Trend (%)'

Title = 'CO2 1944-1970 Yearly Global Mean % Trend'
plot, REFORM(perdiffsfitCO2YearlyGMean_1944_1970),pres,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=8,xrange=[-30.,25.],/ylog, yrange=[pres[nPres-19],pres[3]],xmargin=[10,6]
AXIS, YAXIS=1, YSTYLE = 1, YTITLE = 'Mean Altitude (km)',YRANGE=[geomAltMean[nPres-19],geomAltMean[3]]

;diffsfitCO2YearlyGMean_1970_2000 = fitCO2YearlyGMean_1970_2000[nYears-26,*] - fitCO2YearlyGMean_1970_2000[0,*]
;
;perdiffsfitCO2YearlyGMean_1970_2000 = FLTARR(nPres)
;
;FOR iPres=0,nPres-1 DO BEGIN
;
;  perdiffsfitCO2YearlyGMean_1970_2000[iPres] = diffsfitCO2YearlyGMean_1970_2000[iPres] / $
;                                         ((CO2YearlyGlobalMean[nYears-11,iPres] + CO2YearlyGlobalMean[0,iPres]) / 2.) * 100.
;
;ENDFOR

;xTitle = 'CO2 Trend (%)'
;
;Title = 'CO2 1970-2000 Yearly Global Mean % Trend'
;plot, REFORM(perdiffsfitCO2YearlyGMean_1970_2000),pres,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=8,xrange=[12.,30.],/ylog, yrange=[pres[nPres-19],pres[3]],xmargin=[10,6]
;AXIS, YAXIS=1, YSTYLE = 1, YTITLE = 'Mean Altitude (km)',YRANGE=[geomAltMean[nPres-19],geomAltMean[3]]

xTitle = 'T Trend (K)'

diffsfitTYearlyGMean_1944_1970 = fitTYearlyGMean_1944_1970[nYears-1,*] - fitTYearlyGMean_1944_1970[0,*]
diffsdecfitTYearlyGMean_1944_1970 = diffsfitTYearlyGMean_1944_1970 / 5.5

Title = 'T 1944-1970 Yearly Global Mean Residual Trend '
plot, REFORM(diffsfitTYearlyGMean_1944_1970),pres,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=8,xrange=[-20.,30.],/ylog, yrange=[pres[nPres-19],pres[3]],xmargin=[10,6]
AXIS, YAXIS=1, YSTYLE = 1, YTITLE = 'Mean Altitude (km)',YRANGE=[geomAltMean[nPres-19],geomAltMean[3]]

xTitle = 'T Trend (K/decade)'

Title = 'T 1944-1970 Yearly Global Mean Residual Trend Per Decade'
plot, REFORM(diffsdecfitTYearlyGMean_1944_1970),pres,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=8,xrange=[-4.,5.],/ylog, yrange=[pres[nPres-19],pres[3]],xmargin=[10,6]
AXIS, YAXIS=1, YSTYLE = 1, YTITLE = 'Mean Altitude (km)',YRANGE=[geomAltMean[nPres-19],geomAltMean[3]]

;diffsfitTYearlyGMean_1970_2000 = fitTYearlyGMean_1970_2000[nYears-26,*] - fitTYearlyGMean_1970_2000[0,*]
;diffsdecfitTYearlyGMean_1970_2000 = diffsfitTYearlyGMean_1970_2000 / 3.0
;
;xTitle = 'T Trend (K)'
;
;Title = 'T 1970-2000 Yearly Global Mean Residual Trend '
;plot, REFORM(diffsfitTYearlyGMean_1970_2000),pres,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=8,xrange=[-15.,5.],/ylog, yrange=[pres[nPres-19],pres[3]],xmargin=[10,6]
;AXIS, YAXIS=1, YSTYLE = 1, YTITLE = 'Mean Altitude (km)',YRANGE=[geomAltMean[nPres-19],geomAltMean[3]]
;
;xTitle = 'T Trend (K/decade)'
;
;Title = 'T 1970-2000 Yearly Global Mean Residual Trend Per Decade'
;plot, REFORM(diffsdecfitTYearlyGMean_1970_2000),pres,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=8,xrange=[-5.,2.],/ylog, yrange=[pres[nPres-19],pres[3]],xmargin=[10,6]
;AXIS, YAXIS=1, YSTYLE = 1, YTITLE = 'Mean Altitude (km)',YRANGE=[geomAltMean[nPres-19],geomAltMean[3]]
;
; Do NO
;
diffsfitNOYearlyGMean_1944_1970 = fitNOYearlyGMean_1944_1970[nYears-1,*] - fitNOYearlyGMean_1944_1970[0,*]
;
perdiffsfitNOYearlyGMean_1944_1970 = FLTARR(nPres-1)

FOR iPres=0,nPres-2 DO BEGIN

  perdiffsfitNOYearlyGMean_1944_1970[iPres] = diffsfitNOYearlyGMean_1944_1970[iPres] / $
                                         ((NOYearlyGlobalMean[nYears-1,iPres] + NOYearlyGlobalMean[0,iPres]) / 2.) * 100.

ENDFOR

xTitle = 'NO Residual Trend (%)'

Title = 'NO 1944-1970 Yearly Global Mean Residual % Trend'
plot, REFORM(perdiffsfitNOYearlyGMean_1944_1970),pres,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=8,xrange=[-15.,50.],/ylog, yrange=[pres[nPres-19],pres[3]],xmargin=[10,6]
AXIS, YAXIS=1, YSTYLE = 1, YTITLE = 'Mean Altitude (km)',YRANGE=[geomAltMean[nPres-19],geomAltMean[3]]

;diffsfitNOYearlyGMean_1970_2000 = fitNOYearlyGMean_1970_2000[nYears-26,*] - fitNOYearlyGMean_1970_2000[0,*]
;
;perdiffsfitNOYearlyGMean_1970_2000 = FLTARR(nPres-1)
;
;FOR iPres=0,nPres-2 DO BEGIN
;
;  perdiffsfitNOYearlyGMean_1970_2000[iPres] = diffsfitNOYearlyGMean_1970_2000[iPres] / $
;                                         ((NOYearlyGlobalMean[nYears-11,iPres] + NOYearlyGlobalMean[0,iPres]) / 2.) * 100.
;
;ENDFOR
;
;xTitle = 'NO Residual Trend (%)'
;
;Title = 'NO 1970-2000 Yearly Global Mean Residual % Trend'
;plot, REFORM(perdiffsfitNOYearlyGMean_1970_2000),pres,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=8,xrange=[-25.,15.],/ylog, yrange=[pres[nPres-19],pres[3]],xmargin=[10,6]
;AXIS, YAXIS=1, YSTYLE = 1, YTITLE = 'Mean Altitude (km)',YRANGE=[geomAltMean[nPres-19],geomAltMean[3]],ylog=0

device,/close

END