;waccmx_50yr_yearly_pres_2010.pro
;
;Calculate and plot global monthly and yearly mean WACCM-X 1955-2005 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  
;
;
;Open first file to get lat/lon/pressure grid
;
fid = ncdf_open("/glade/scratch/joemci/f_1955-2005_waccmx_cesm1_1_beta08_pleiades/atm/hist/f_1955-2005_waccmx_cesm1_1_beta08.cam.h0.1955-01.nc")

ncdf_varget, fid, 'lat', lat
nLat = n_elements(lat)

dLat = FLTARR(nLat)

areaWeight = FLTARR(nLat)

ncdf_varget, fid, 'lon', lon
nLon = n_elements(lon)

ncdf_varget, fid, 'lev', pres
nPres = n_elements(pres)

ncdf_close,fid
;
;  Get area weights using spherical cap method.  Weights used to calculate global means
;
earthRadius = 6371.220 

FOR iLat = 0,nLat-1 DO BEGIN

  latForArea = ABS(lat)

  if iLat eq 0 then begin

    latArea = 2. * !PI * earthRadius*earthRadius * (1. - SIN(latForArea(1) * !PI/ 180.))

  endif else if iLat eq nLat-1 then begin

    latArea = 2. * !PI * earthRadius*earthRadius * (1. - SIN(latForArea(nLat-2) * !PI/ 180.))

  endif else begin 

    lat1 = (latForArea(iLat-1) + latForArea(iLat)) / 2.
    lat2 = (latForArea(iLat) + latForArea(iLat+1)) / 2.  
    latArea = 2. * !PI * earthRadius*earthRadius * ABS( (1. - SIN(lat1 * !PI/ 180.)) - (1. - SIN(lat2 * !PI/ 180.)) )

    if lat(iLat) lt 0. and lat(iLat+1) gt 0. then latArea = 2. * !PI * earthRadius*earthRadius * (1. - SIN(lat1 * !PI/ 180.))
    if lat(iLat) gt 0. and lat(iLat+1) lt 0. then latArea = 2. * !PI * earthRadius*earthRadius * (1. - SIN(lat2 * !PI/ 180.))

  endelse

  areaWeight[iLat] = latArea / nLon
  
ENDFOR

sumAreaWeight = TOTAL(areaWeight) * nLon

nYears = 56
nMonths = 12
iMonthAll = 0

CO2MonthlyAllGlobalMean = FLTARR(nYears*nMonths,nPres)
TMonthlyAllGlobalMean = FLTARR(nYears*nMonths,nPres)
geomAltMonthlyAllGlobalMean = FLTARR(nYears*nMonths,nPres)

CO2MonthlyGlobalMean = FLTARR(nMonths,nPres)
TMonthlyGlobalMean = FLTARR(nMonths,nPres)
geomAltMonthlyGlobalMean = FLTARR(nMonths,nPres)

CO2YearlyGlobalMean = FLTARR(nYears,nPres)
TYearlyGlobalMean = FLTARR(nYears,nPres)
geomAltYearlyGlobalMean = FLTARR(nYears,nPres)

CO2AllMean = FLTARR(nYears*nMonths,4)
TAllMean = FLTARR(nYears*nMonths,4)
UZonalMean = FLTARR(nYears*nMonths,4)

CO2WeightedMean = FLTARR(nPres)
TWeightedMean = FLTARR(nPres)
geomAltWeightedMean = FLTARR(nPres)
;
;Loop over years and months and read CO2 and temperature to compile 50 year monthly and yearly global mean values
;
FOR iYear=0,nYears-1 DO BEGIN
    
  CO2MonthlyGlobalMean[*,*] = 0.
  TMonthlyGlobalMean[*,*] = 0.
  geomAltMonthlyGlobalMean[*,*] = 0.

  cYear = STRTRIM(iYear+1955,2)
  print, 'Doing year ', cYear
 
  FOR iMonth=0,nMonths-1 DO BEGIN

    CO2WeightedMean[*] = 0.
    TWeightedMean[*] = 0.
    geomAltWeightedMean[*] = 0.
 
    ;
    ;Open file, read monthly file CO2 and T, close file
    ;
    cMonth = STRTRIM(iMonth+1,2)
    cYear = STRTRIM(iYear+1955,2)
  
    if iMonth lt 9 then cMonth = '0'+cMonth

;    print, 'Opening/reading file for month/year ', cMonth, '/', cYear
  
    fid = ncdf_open("/glade/scratch/joemci/f_1955-2005_waccmx_cesm1_1_beta08_pleiades/atm/hist/f_1955-2005_waccmx_cesm1_1_beta08.cam.h0."+cYear+"-"+cMonth+".nc")

    ncdf_varget, fid, 'CO2', CO2In  ; [lon,lat,vert]
    ncdf_varget, fid, 'T', TIn  ; [lon,lat,vert]
    ncdf_varget, fid, 'Z3', geopAlt

    geomAlt = geopAlt * (1. + geopAlt/earthRadius/1000.) /1000.    
;
;  Apply area weighting for each latitude and longitude to get global mean at each pressure level
;
    FOR iPres = 0,nPres-1 DO BEGIN
    
      FOR iLat = 0,nLat-1 DO BEGIN

        FOR iLon = 0,nLon-1 DO BEGIN
	      
	  CO2WeightedMean[iPres] = CO2WeightedMean[iPres] + CO2In[ILon,iLat,iPres] * areaWeight[iLat]
	  TWeightedMean[iPres] = TWeightedMean[iPres] + TIn[ILon,iLat,iPres] * areaWeight[iLat]
	  geomAltWeightedMean[iPres] = geomAltWeightedMean[iPres] + geomAlt[ILon,iLat,iPres] * areaWeight[iLat]
	  
        ENDFOR ;Lon

      ENDFOR ;Lat
    
      CO2MonthlyGlobalMean[iMonth,iPres] = CO2WeightedMean[iPres] / sumAreaWeight
      TMonthlyGlobalMean[iMonth,iPres] = TWeightedMean[iPres] / sumAreaWeight
      geomAltMonthlyGlobalMean[iMonth,iPres] = geomAltWeightedMean[iPres] / sumAreaWeight
    
      CO2MonthlyAllGlobalMean[iMonthAll,iPres] = CO2WeightedMean[iPres] / sumAreaWeight
      TMonthlyAllGlobalMean[iMonthAll,iPres] = TWeightedMean[iPres]/ sumAreaWeight
      geomAltMonthlyAllGlobalMean[iMonthAll,iPres] = geomAltWeightedMean[iPres]/ sumAreaWeight

    ENDFOR ;Pres
    
    iMonthAll = iMonthAll + 1
;
;  Get lowest top geometric altitude for whole 56 years to use in choosing of vertical levels to do real global mean
;    
    geomAltMonthlyGlobalMeanTop = geomAltMonthlyGlobalMean[iMonth,1]    
    
;    print, 'Geometric altitude at top level ', iMonth, geomAltAllMeanTop
    if geomAltMonthlyGlobalMeanTop lt 320. then print, 'Monthly geometric altitude at top level below 320', iMonth, geomAltMonthlyGlobalMeanTop
    if geomAltMonthlyGlobalMeanTop lt 350. then print, 'Monthly geometric altitude at top level below 350', iMonth, geomAltMonthlyGlobalMeanTop
;    if geomAltMonthlyGlobalMeanTop lt 400. then print, 'Monthly geometric altitude at top level below 400', iMonth, geomAltMonthlyGlobalMeanTop
;    if geomAltMonthlyGlobalMeanTop lt 450. then print, 'Monthly geometric altitude at top level below 450', iMonth, geomAltMonthlyGlobalMeanTop

    ncdf_close,fid

  ENDFOR ;Month
;
;  Get yearly mean from monthly means for each pressure level
;
  FOR iPres = 0,nPres-1 DO BEGIN

    CO2YearlyGlobalMean[iYear,iPres] = MEAN(CO2MonthlyGlobalMean[*,iPres])
    TYearlyGlobalMean[iYear,iPres] = MEAN(TMonthlyGlobalMean[*,iPres])
    geomAltYearlyGlobalMean[iYear,iPres] = MEAN(geomAltMonthlyGlobalMean[*,iPres])

  ENDFOR ;Pres

;print, 'TYearlyGlobalMean Year ', iYear, TYearlyGlobalMean[iYear,*]

    geomAltYearlyGlobalMeanTop = geomAltYearlyGlobalMean[iYear,1]    
    print, 'Yearly geometric altitude at top level ', iYear, geomAltYearlyGlobalMeanTop

ENDFOR ;Year
;
;  Plot monthly and yearly global mean data
;
SET_PLOT, 'PS'

psfile = 'monthly_CO2_T_AreaGlobalMean_WACCMX_Pres_2010.ps'

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

MonthsAll = INDGEN(iMonthAll)

FOR iPres = 0,nPres-1 DO BEGIN

 yTitle = 'CO2 mol/mol'
 xTitle = 'Months Since January, 1955'
 
 cPres = STRTRIM(pres[iPres],2)

 Title = 'CO2 1955-2010 Area Weighted Monthly Global Mean  '+cPres+' km'
 plot, MonthsAll, CO2MonthlyAllGlobalMean[0:iMonthAll-1,iPres],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

 yTitle = 'T (K)'

 Title = 'T 1955-2010 Area Weighted Monthly Global Mean  '+cPres+' km'
 plot, MonthsAll, TMonthlyAllGlobalMean[0:iMonthAll-1,iPres],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

ENDFOR ; Pressures

device,/close

psfile = 'yearly_CO2_T_AreaGlobalMean_WACCMX_Pres_2010.ps'

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

YearsAll = INDGEN(nYears)

FOR iPres = 0,nPres-1 DO BEGIN

 yTitle = 'CO2 mol/mol'
 xTitle = 'Years Since 1955'
 
 cPres = STRTRIM(pres[iPres],2)

 Title = 'CO2 1955-2010 Area Weighted Yearly Global Mean  '+cPres+' km'
 plot, YearsAll, CO2YearlyGlobalMean[0:nYears-1,iPres],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

 yTitle = 'T (K)'
 
 Title = 'T 1955-2010 Area Weighted Yearly Global Mean '+cPres+' km'
 plot, YearsAll, TYearlyGlobalMean[0:nYears-1,iPres],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

ENDFOR ; Pressures

device,/close

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

psfile = 'year_pres_CO2_T_50yr_AreaGlobalMean_WACCMX_2010.ps'

DEVICE, /COLOR, FILENAME = psfile
loadct, 39

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

levels = 21
clevels = findgen(levels) * 2.E-05
  
Title = ' WACCM-X CO2 1955-2010 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 1955-2010 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 1955-2010 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 1955-2010 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 1955-2010 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 1955-2010 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 1955-2010 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 1955-2010 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_2010, f107AllMonthly, apAllMonthly, f107YearlyMean, apYearlyMean
;
;  Now do regression of yearly mean temperature with f10.7, ap, and f10.7*f10.7
;
f107AllMonthlySquared = f107AllMonthly * f107AllMonthly
f107YearlyMeanSquared = f107YearlyMean * f107YearlyMean

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

psfile = 'yearly_CO2_T_AreaGlobalMean_WACCMX_Pres_Trends_2010.ps'

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

FOR iPres = 0,nPres-1 DO BEGIN

  CO2YearlyGMean_1955_2005 = CO2YearlyGlobalMean[0:nYears-1,iPres]
  regressCoeffs_1955_2005 = REGRESS(f107apYearlyMean_1955_2005, CO2YearlyGMean_1955_2005, YFIT=CO2YearlyGMeanFit_1955_2005)
  residualsCO2YearlyGMean_1955_2005 = CO2YearlyGMean_1955_2005 - CO2YearlyGMeanFit_1955_2005

  residualsCO2YearlyGMean_1970_2000 = residualsCO2YearlyGMean_1955_2005[15:nYears-6]
  ;
  ;  Get fit to residuals
  ;
  Years_1955_2005 = INDGEN(nYears)
  fitCoeffs = POLY_FIT(Years_1955_2005,residualsCO2YearlyGMean_1955_2005,1,YFIT=fitCO2YearlyGMean_1955_2005)
 
  Years_1970_2000 = INDGEN(nYears-20)
  fitCoeffs_1970_2000 = POLY_FIT(Years_1970_2000,residualsCO2YearlyGMean_1970_2000,1,YFIT=fitCO2YearlyGMean_1970_2000) 

  yTitle = 'CO2 mol/mol'
  xTitle = 'Years Since 1955'
 
  cPres = STRTRIM(pres[iPres],2)

  Title = 'CO2 1955-2010 Area Weighted Yearly Global Mean Residuals & Fit '+cPres+' km'
  plot, Years_1955_2005,residualsCO2YearlyGMean_1955_2005,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
  oplot, fitCO2YearlyGMean_1955_2005

  xTitle = 'Years Since 1970'

  Title = 'CO2 1970-2000 Area Weighted Yearly Global Mean Residuals & Fit '+cPres+' km'
  plot, Years_1970_2000,residualsCO2YearlyGMean_1970_2000,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
  oplot, fitCO2YearlyGMean_1970_2000

  TYearlyGMean_1955_2005 = TYearlyGlobalMean[0:nYears-1,iPres]
  regressCoeffs_1955_2005 = REGRESS(f107apYearlyMean_1955_2005, TYearlyGMean_1955_2005, YFIT=TYearlyGMeanFit_1955_2005)
  residualsTYearlyGMean_1955_2005 = TYearlyGMean_1955_2005 - TYearlyGMeanFit_1955_2005

  residualsTYearlyGMean_1970_2000 = residualsTYearlyGMean_1955_2005[15:nYears-6]
;
;  Get fit to residuals
;
  Years_1955_2005 = INDGEN(nYears)
  fitCoeffs = POLY_FIT(Years_1955_2005,residualsTYearlyGMean_1955_2005,1,YFIT=fitTYearlyGMean_1955_2005) 

  Years_1970_2000 = INDGEN(nYears-20)
  fitCoeffs = POLY_FIT(Years_1970_2000,residualsTYearlyGMean_1970_2000,1,YFIT=fitTYearlyGMean_1970_2000) 

  yTitle = 'T (K)'

  xTitle = 'Years Since 1955'

  Title = 'T 1955-2010 Area Weighted Yearly Global Mean Residuals & Fit '+cPres+' km'
  plot, Years_1955_2005,residualsTYearlyGMean_1955_2005,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
  oplot, fitTYearlyGMean_1955_2005

  xTitle = 'Years Since 1970'

  Title = 'T 1970-2000 Area Weighted Yearly Global Mean Residuals & Fit '+cPres+' km'
  plot, Years_1970_2000,residualsTYearlyGMean_1970_2000,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
  oplot, fitTYearlyGMean_1970_2000

ENDFOR

device,/close
;
;  Now do regression of monthly mean temperature with f10.7 and ap
;
f107apMonthlyMean_1955_2005 = FLTARR(3,iMonthAll)
f107apMonthlyMean_1955_2005[0,*] = f107AllMonthly[0:iMonthAll-1]
f107apMonthlyMean_1955_2005[1,*] = apAllMonthly[0:iMonthAll-1]
f107apMonthlyMean_1955_2005[2,*] = f107AllMonthlySquared[0:iMonthAll-1]

f107apMonthlyMean_1970_2000 = FLTARR(3,iMonthAll-180)
f107apMonthlyMean_1970_2000[0,*] = f107AllMonthly[180:iMonthAll-1]
f107apMonthlyMean_1970_2000[1,*] = apAllMonthly[180:iMonthAll-1]
f107apMonthlyMean_1970_2000[2,*] = f107AllMonthlySquared[180:iMonthAll-1]

psfile = 'monthly_CO2_T_AreaGlobalMean_WACCMX_Pres_Trends_2010.ps'

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

FOR iPres = 0,nPres-1 DO BEGIN

  CO2MonthlyAllGMean_1955_2005 = CO2MonthlyAllGlobalMean[*,iPres]
  regressCoeffs = REGRESS(f107apMonthlyMean_1955_2005, CO2MonthlyAllGMean_1955_2005, YFIT=CO2MonthlyGMeanFit_1955_2005)
  residualsCO2MonthlyGMean_1955_2005 = CO2MonthlyAllGMean_1955_2005 - CO2MonthlyGMeanFit_1955_2005

  residualsCO2MonthlyGMean_1970_2000 = residualsCO2MonthlyGMean_1955_2005[180:iMonthAll-61]
  ;
  ;  Get fit to residuals
  ;
  Months_1955_2005 = INDGEN(iMonthAll)
  fitCoeffs_1955_2005 = POLY_FIT(Months_1955_2005,residualsCO2MonthlyGMean_1955_2005,1,YFIT=fitCO2MonthlyGMean_1955_2005) 

  Months_1970_2000 = INDGEN(iMonthAll-240)
  fitCoeffs_1970_2000 = POLY_FIT(Months_1970_2000,residualsCO2MonthlyGMean_1970_2000[0:iMonthAll-241],1,YFIT=fitCO2MonthlyGMean_1970_2000) 

  yTitle = 'CO2 mol/mol'
  xTitle = 'Months Since January, 1955'
 
  cPres = STRTRIM(pres[iPres],2)

  Title = 'CO2 1955-2010 Area Weighted Monthly Global Mean Residuals and Fit '+cPres+' km'
  plot, Months_1955_2005,residualsCO2MonthlyGMean_1955_2005,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
  oplot, fitCO2MonthlyGMean_1955_2005

  xTitle = 'Months Since January, 1970'

  Title = 'CO2 1970-2000 Area Weighted Monthly Global Mean Residuals and Fit '+cPres+' km'
  plot, Months_1970_2000,residualsCO2MonthlyGMean_1970_2000,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
  oplot, fitCO2MonthlyGMean_1970_2000

  TMonthlyAllGMean_1955_2005 = TMonthlyAllGlobalMean[*,iPres]
  regressCoeffs = REGRESS(f107apMonthlyMean_1955_2005, TMonthlyAllGMean_1955_2005, YFIT=TMonthlyGMeanFit_1955_2005)
  residualsTMonthlyGMean_1955_2005= TMonthlyAllGMean_1955_2005 - TMonthlyGMeanFit_1955_2005

  residualsTMonthlyGMean_1970_2000 = residualsTMonthlyGMean_1955_2005[180:iMonthAll-61]
;
;  Get fit to residuals
;
  Months_1955_2005 = INDGEN(iMonthAll)
  fitCoeffs_1955_2005 = POLY_FIT(Months_1955_2005,residualsTMonthlyGMean_1955_2005,1,YFIT=fitTMonthlyGMean_1955_2005) 

  Months_1970_2000 = INDGEN(iMonthAll-240)
  fitCoeffs_1970_2000 = POLY_FIT(Months_1970_2000,residualsTMonthlyGMean_1970_2000[0:iMonthAll-241],1,YFIT=fitTMonthlyGMean_1970_2000) 

  yTitle = 'T (K)'

  Title = 'T 1955-2010 Area Weighted Monthly Global Mean Residuals and Fit '+cPres+' km'
  plot, Months_1955_2005,residualsTMonthlyGMean_1955_2005,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
  oplot, fitTMonthlyGMean_1955_2005

  xTitle = 'Months Since January, 1970'

  Title = 'T 1970-2000 Area Weighted Monthly Global Mean Residuals and Fit '+cPres+' km'
  plot, Months_1970_2000,residualsTMonthlyGMean_1970_2000,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
  oplot, fitTMonthlyGMean_1970_2000

ENDFOR

device,/close

END