;waccmx_50yr_yearly_pres.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 regression
;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("/hao/aim3/joemci/WACCMX/wa3548_2x_refb1.2.cam2.h0/atm/hist/Subset/wa3548_2x_refb1.2.cam2.h0.1955-01_TCO2.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 where need to calculate area above one latitude and subtract from area of next lower 
;  latitude and get area for latitude band.  Then divide by number of longitudes to get grid area weight latitude/longitude point.
;  Area 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 = 51
nMonths = 12
iMonthAll = 0

;nNewGeomAlt = 9 ;50km
;newGeomAlt = REVERSE(FINDGEN(nNewGeomAlt)*50.)
;nNewGeomAlt = 17 ;25km
;newGeomAlt = REVERSE(FINDGEN(nNewGeomAlt)*50.) / 2.
;nNewGeomAlt = 81 ;10km
nNewGeomAlt = 28 ;5km
newGeomAlt = REVERSE(FINDGEN(nNewGeomAlt)*5.)
;newGeomAlt[nNewGeomAlt-1] = newGeomAlt[nNewGeomAlt-1] + .1  ; Don't make bottom point zero km altitude

CO2MonthlyAllGlobalMean = FLTARR(nYears*nMonths,nNewGeomAlt)
TMonthlyAllGlobalMean = FLTARR(nYears*nMonths,nNewGeomAlt)

CO2MonthlyGlobalMean = FLTARR(nMonths,nNewGeomAlt)
TMonthlyGlobalMean = FLTARR(nMonths,nNewGeomAlt)

CO2InterpAlt = FLTARR(nLon,nLat,nNewGeomAlt)
TInterpAlt = FLTARR(nLon,nLat,nNewGeomAlt)

CO2WeightedMean = FLTARR(nNewGeomAlt)
TWeightedMean = FLTARR(nNewGeomAlt)

CO2YearlyGlobalMean = FLTARR(nYears,nNewGeomAlt)
TYearlyGlobalMean = FLTARR(nYears,nNewGeomAlt)
;
;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.

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

    CO2WeightedMean[*] = 0.
    TWeightedMean[*] = 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("/hao/aim3/joemci/WACCMX/wa3548_2x_refb1.2.cam2.h0/atm/hist/Subset/wa3548_2x_refb1.2.cam2.h0."+cYear+"-"+cMonth+"_TCO2.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.    
;
;  Do interpolation to altitude setting missing points to surface value
;    
    FOR iLat = 0,nLat-1 DO BEGIN

      FOR iLon = 0,nLon-1 DO BEGIN

        vfnd, geomAlt[iLon,iLat,*], CO2In[iLon,iLat,*], newGeomAlt, outData, cut=CO2In[iLon,iLat,nNewGeomAlt-1], linear=1
	CO2InterpAlt[iLon,iLat,*] = outData
        vfnd, geomAlt[iLon,iLat,*], TIn[iLon,iLat,*], newGeomAlt, outData, cut=TIn[iLon,iLat,nNewGeomAlt-1], linear=1
	TInterpAlt[iLon,iLat,*] = outData
;
;  Check for when input maximum altitude for this point is below interpolated altitude and fill with top valid point
;
        FOR iAlt = nNewGeomAlt-1,0,-1 DO BEGIN

          if MAX(geomAlt[iLon,iLat,*]) lt newGeomAlt[iAlt] then CO2InterpAlt[iLon,iLat,iAlt] = CO2InterpAlt[iLon,iLat,iAlt+1]
  
        ENDFOR
        
      ENDFOR ;Lon

    ENDFOR ;Lat
;
;  Apply area weighting for each latitude and longitude to get global mean at each altitude level
;
    FOR iAlt = 0,nNewGeomAlt-1 DO BEGIN
    
      FOR iLat = 0,nLat-1 DO BEGIN

        FOR iLon = 0,nLon-1 DO BEGIN
	      
	  CO2WeightedMean[iAlt] = CO2WeightedMean[iAlt] + CO2InterpAlt[ILon,iLat,iAlt] * areaWeight[iLat]
	  TWeightedMean[iAlt] = TWeightedMean[iAlt] + TInterpAlt[ILon,iLat,iAlt] * areaWeight[iLat]
	  
        ENDFOR ;Lon

      ENDFOR ;Lat
    
      CO2MonthlyGlobalMean[iMonth,iAlt] = CO2WeightedMean[iAlt] / sumAreaWeight
      TMonthlyGlobalMean[iMonth,iAlt] = TWeightedMean[iAlt] / sumAreaWeight
    
      CO2MonthlyAllGlobalMean[iMonthAll,iAlt] = CO2WeightedMean[iAlt] / sumAreaWeight
      TMonthlyAllGlobalMean[iMonthAll,iAlt] = TWeightedMean[iAlt]/ sumAreaWeight

    ENDFOR ;Altitude
    
    iMonthAll = iMonthAll + 1
    
;    print, 'Geometric altitude at top level ', iMonth, geomAltAllMeanTop
;    if geomAltMonthlyGlobalMeanTop lt 350. then print, 'Monthly geometric altitude at top level below 350', iMonth, geomAltMonthlyGlobalMeanTop
;    if geomAltMonthlyGlobalMeanTop gt 350. then print, 'Monthly geometric altitude at top level above 350', iMonth, geomAltMonthlyGlobalMeanTop
;    if geomAltMonthlyGlobalMeanTop gt 400. then print, 'Monthly geometric altitude at top level above 400', iMonth, geomAltMonthlyGlobalMeanTop
;    if geomAltMonthlyGlobalMeanTop gt 450. then print, 'Monthly geometric altitude at top level above 450', iMonth, geomAltMonthlyGlobalMeanTop

    ncdf_close,fid

  ENDFOR ;Month
;
;  Get yearly mean from monthly means for each altitude level
;
  FOR iAlt = 0,nNewGeomAlt-1 DO BEGIN

    CO2YearlyGlobalMean[iYear,iAlt] = MEAN(CO2MonthlyGlobalMean[*,iAlt])
    TYearlyGlobalMean[iYear,iAlt] = MEAN(TMonthlyGlobalMean[*,iAlt])

  ENDFOR ;Altitude

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

ENDFOR ;Year

;tempout = TMonthlyAllGlobalMean[*,7]
;
;fname='temperature_280km_globalmonthly_1955-2005.txt'
;OPENW,1,fname
;PRINTF,1,tempout,FORMAT='(F8.3,6X,F8.3,6X,F8.3,6X,F8.3,6X,F8.3,6X,F8.3,6X,F8.3,6X,F8.3,6X,F8.3)'
;;PRINTF,TMonthlyAllGlobalMean[*,7],FORMAT='(F7.3,6X,F7.3,6X,F7.3,6X,F7.3,6X,F7.3,6X,F7.3,6X,F7.3,6X,F7.3,6X,F7.3)'
;;PRINTF,TMonthlyAllGlobalMean[*,7]
;CLOSE,1
;
;tempout = TYearlyGlobalMean[*,7]
;
;fname='temperature_280km_globalyearly_1955-2005.txt'
;OPENW,1,fname
;PRINTF,1,tempout,FORMAT='(F8.3,6X,F8.3,6X,F8.3,6X,F8.3,6X,F8.3,6X,F8.3,6X,F8.3,6X,F8.3,6X,F8.3)'
;;PRINTF,TYearlyGlobalMean[*,7],FORMAT='(F7.3,6X,F7.3,6X,F7.3,6X,F7.3,6X,F7.3,6X,F7.3,6X,F7.3,6X,F7.3,6X,F7.3)'
;;PRINTF,TYearlyGlobalMean[*,7]
;CLOSE,1

SET_PLOT, 'PS'

psfile = 'monthly_CO2_T_AreaGlobalMean_WACCM_Alt_refb1.ps'

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

MonthsAll = INDGEN(iMonthAll)

FOR iAlt = 0,nNewGeomAlt-1 DO BEGIN

 yTitle = 'CO2 mol/mol'
 xTitle = 'Months Since January, 1955'
 
 cAlt = STRTRIM(newGeomAlt[iAlt],2)

 Title = 'CO2 1955-2005 Refb1 Area Weighted Monthly Global Mean  '+cAlt+' km'
 plot, MonthsAll, CO2MonthlyAllGlobalMean[0:iMonthAll-1,iAlt],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'CO2 1955-2005 Area Weighted Monthly Global Mean ~350km'
;plot, MonthsAll, CO2MonthlyAllGlobalMean[0:iMonthAll-1,1],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'CO2 1955-2005 Area Weighted Monthly Global Mean ~300km'
;plot, MonthsAll, CO2MonthlyAllGlobalMean[0:iMonthAll-1,2],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'CO2 1955-2005 Area Weighted Monthly Global Mean ~250km'
;plot, MonthsAll, CO2MonthlyAllGlobalMean[0:iMonthAll-1,3],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

;yrange=[0.0003,0.00038],
;yrange=[0.0003,0.00038],
;yrange=[0.0003,0.00038],


 yTitle = 'T (K)'

 Title = 'T 1955-2005 Refb1 Area Weighted Monthly Global Mean  '+cAlt+' km'
 plot, MonthsAll, TMonthlyAllGlobalMean[0:iMonthAll-1,iAlt],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'T 1955-2005 Area Weighted Monthly Global Mean ~350km'
;plot, MonthsAll, TMonthlyAllGlobalMean[0:iMonthAll-1,1],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'T 1955-2005 Area Weighted Monthly Global Mean ~300km'
;plot, MonthsAll, TMonthlyAllGlobalMean[0:iMonthAll-1,2],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'T 1955-2005 Area Weighted Monthly Global Mean ~250km'
;plot, MonthsAll, TMonthlyAllGlobalMean[0:iMonthAll-1,3],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

;yrange=[190.0,220.0],
;yrange=[222.0,238.0],
;yrange=[285.0,300.0],

ENDFOR ; Altitudes

device,/close

psfile = 'yearly_CO2_T_AreaGlobalMean_WACCM_Alt_refb1.ps'

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

YearsAll = INDGEN(nYears)

FOR iAlt = 0,nNewGeomAlt-1 DO BEGIN

 yTitle = 'CO2 mol/mol'
 xTitle = 'Years Since 1955'
 
 cAlt = STRTRIM(newGeomAlt[iAlt],2)

 Title = 'CO2 1955-2005 Refb1 Area Weighted Yearly Global Mean  '+cAlt+' km'
 plot, YearsAll, CO2YearlyGlobalMean[0:nYears-1,iAlt],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'CO2 1955-2005 Area Weighted Yearly Global Mean ~350km'
;plot, YearsAll, CO2YearlyGlobalMean[0:nYears-1,1],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'CO2 1955-2005 Area Weighted Yearly Global Mean ~300km'
;plot, YearsAll, CO2YearlyGlobalMean[0:nYears-1,2],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'CO2 1955-2005 Area Weighted Yearly Global Mean 250km'
;plot, YearsAll, CO2YearlyGlobalMean[0:nYears-1,3],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

;yrange=[0.0003,0.00038],
;yrange=[0.0003,0.00038],
;yrange=[0.0003,0.00038],


 yTitle = 'T (K)'
 
 Title = 'T 1955-2005 Refb1 Area Weighted Yearly Global Mean '+cAlt+' km'
 plot, YearsAll, TYearlyGlobalMean[0:nYears-1,0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'T 1955-2005 Area Weighted Yearly Global Mean ~350km'
;plot, YearsAll, TYearlyGlobalMean[0:nYears-1,1],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'T 1955-2005 Area Weighted Yearly Global Mean ~300km'
;plot, YearsAll, TYearlyGlobalMean[0:nYears-1,2],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'T 1955-2005 Area Weighted Yearly Global Mean ~250km'
;plot, YearsAll, TYearlyGlobalMean[0:nYears-1,3],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

;yrange=[190.0,220.0],
;yrange=[222.0,238.0],
;yrange=[285.0,300.0],

ENDFOR ; Altitudes

device,/close
;
;  Get monthly and yearly mean f107 and ap for regression
;
waccmx_solar_f107_ap_monthly_yearly, 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]

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

psfile = 'yearly_CO2_T_AreaGlobalMean_WACCM_Alt_Trends_refb1.ps'

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

FOR iAlt = 0,nNewGeomAlt-1 DO BEGIN

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

;  CO2YearlyGMean_1970_2000 = CO2YearlyGlobalMean[15:nYears-6,iAlt]
;  regressCoeffs_1970_2000 = REGRESS(f107apYearlyMean_1970_2000, CO2YearlyGMean_1970_2000, YFIT=CO2YearlyGMeanFit_1970_2000)
;  residualsCO2YearlyGMean_1970_2000 = CO2YearlyGMean_1970_2000 - CO2YearlyGMeanFit_1970_2000
   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'
 
  cAlt = STRTRIM(newGeomAlt[iAlt],2)

  Title = 'CO2 1955-2005 Refb1 Area Weighted Yearly Global Mean Residuals & Fit '+cAlt+' 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 Refb1 Area Weighted Yearly Global Mean Residuals & Fit '+cAlt+' km'
  plot, Years_1970_2000,residualsCO2YearlyGMean_1970_2000,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
  oplot, fitCO2YearlyGMean_1970_2000

;yrange=[0.0003,0.00038],
;yrange=[0.0003,0.00038],
;yrange=[0.0003,0.00038],

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

;  TYearlyGMean_1970_2000 = TYearlyGlobalMean[15:nYears-6,iAlt]
;  regressCoeffs_1970_2000 = REGRESS(f107apYearlyMean_1970_2000, TYearlyGMean_1970_2000, YFIT=TYearlyGMeanFit_1970_2000)
;  residualsTYearlyGMean_1970_2000 = TYearlyGMean_1970_2000 - TYearlyGMeanFit_1970_2000
  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_1970_2000 = POLY_FIT(Years_1970_2000,residualsTYearlyGMean_1970_2000,1,YFIT=fitTYearlyGMean_1970_2000) 

  yTitle = 'T (K)'

  xTitle = 'Years Since 1955'

  Title = 'T 1955-2005 Refb1 Area Weighted Yearly Global Mean Residuals & Fit '+cAlt+' km'
  plot, Years_1955_2005,residualsTYearlyGMean_1955_2005,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
  oplot, fitTYearlyGMean_1955_2005
;Title = 'T 1955-2005 Area Weighted Yearly Global Mean Residuals and Fit ~350km'
;plot, Years_1955_2005, [0:iMonthAll-1,1],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'T 1955-2005 Area Weighted Yearly Global Mean Residuals and Fit ~300km'
;plot, Years_1955_2005, [0:iMonthAll-1,2],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'T 1955-2005 Area Weighted Yearly Global Mean Residuals and Fit ~250km'
;plot, Years_1955_2005, [0:iMonthAll-1,3],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

  xTitle = 'Years Since 1970'

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

ENDFOR

;yrange=[190.0,220.0],
;yrange=[222.0,238.0],
;yrange=[285.0,300.0],

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-240)
f107apMonthlyMean_1970_2000[0,*] = f107AllMonthly[180:iMonthAll-61]
f107apMonthlyMean_1970_2000[1,*] = apAllMonthly[180:iMonthAll-61]
f107apMonthlyMean_1970_2000[2,*] = f107AllMonthlySquared[180:iMonthAll-61]

psfile = 'monthly_CO2_T_AreaGlobalMean_WACCM_Alt_Trends_refb1.ps'

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

FOR iAlt = 0,nNewGeomAlt-1 DO BEGIN

  CO2MonthlyAllGMean_1955_2005 = CO2MonthlyAllGlobalMean[*,iAlt]
;  CO2MonthlyAllGMean_1970_2000 = CO2MonthlyAllGlobalMean[180:iMonthAll-61,iAlt]
  regressCoeffs = REGRESS(f107apMonthlyMean_1955_2005, CO2MonthlyAllGMean_1955_2005, YFIT=CO2MonthlyGMeanFit_1955_2005)
  residualsCO2MonthlyGMean_1955_2005 = CO2MonthlyAllGMean_1955_2005 - CO2MonthlyGMeanFit_1955_2005
;  regressCoeffs = REGRESS(f107apMonthlyMean_1970_2000, CO2MonthlyAllGMean_1970_2000, YFIT=CO2MonthlyGMeanFit_1970_2000)
;  residualsCO2MonthlyGMean_1970_2000 = CO2MonthlyAllGMean_1970_2000 - CO2MonthlyGMeanFit_1970_2000
  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'
 
  cAlt = STRTRIM(newGeomAlt[iAlt],2)

  Title = 'CO2 1955-2005 Refb1 Area Weighted Monthly Global Mean Residuals and Fit '+cAlt+' 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 Refb1 Area Weighted Monthly Global Mean Residuals and Fit '+cAlt+' km'
  plot, Months_1970_2000,residualsCO2MonthlyGMean_1970_2000,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
  oplot, fitCO2MonthlyGMean_1970_2000

;Title = 'CO2 1955-2005 Area Weighted Monthly Global Mean Residuals and Fit ~400km'
;plot, Months_1955_2005, ,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'CO2 1955-2005 Area Weighted Monthly Global Mean Residuals and Fit ~350km'
;plot, Months_1955_2005, [0:iMonthAll-1,1],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'CO2 1955-2005 Area Weighted Monthly Global Mean Residuals and Fit ~300km'
;plot, Months_1955_2005, [0:iMonthAll-1,2],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'CO2 1955-2005 Area Weighted Monthly Global Mean Residuals and Fit ~250km'
;plot, Months_1955_2005, [0:iMonthAll-1,3],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

;yrange=[0.0003,0.00038],
;yrange=[0.0003,0.00038],
;yrange=[0.0003,0.00038],

  TMonthlyAllGMean_1955_2005 = TMonthlyAllGlobalMean[*,iAlt]
;  TMonthlyAllGMean_1970_2000 = TMonthlyAllGlobalMean[180:iMonthAll-61,iAlt]
  regressCoeffs = REGRESS(f107apMonthlyMean_1955_2005, TMonthlyAllGMean_1955_2005, YFIT=TMonthlyGMeanFit_1955_2005)
  residualsTMonthlyGMean_1955_2005 = TMonthlyAllGMean_1955_2005 - TMonthlyGMeanFit_1955_2005
;  regressCoeffs = REGRESS(f107apMonthlyMean_1970_2000, TMonthlyAllGMean_1970_2000, YFIT=TMonthlyGMeanFit_1970_2000)
;  residualsTMonthlyGMean_1970_2000 = TMonthlyAllGMean_1970_2000 - TMonthlyGMeanFit_1970_2000
  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-2005 Refb1 Area Weighted Monthly Global Mean Residuals and Fit '+cAlt+' km'
  plot, Months_1955_2005,residualsTMonthlyGMean_1955_2005,title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
  oplot, fitTMonthlyGMean_1955_2005
;Title = 'T 1955-2005 Area Weighted Monthly Global Mean Residuals and Fit ~350km'
;plot, Months_1955_2005, [0:iMonthAll-1,1],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'T 1955-2005 Area Weighted Monthly Global Mean Residuals and Fit ~300km'
;plot, Months_1955_2005, [0:iMonthAll-1,2],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
;Title = 'T 1955-2005 Area Weighted Monthly Global Mean Residuals and Fit ~250km'
;plot, Months_1955_2005, [0:iMonthAll-1,3],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

  xTitle = 'Months Since January, 1970'

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

;yrange=[190.0,220.0],
;yrange=[222.0,238.0],
;yrange=[285.0,300.0],

ENDFOR

device,/close

;plot, residualsTMonthlyGMean
;oplot, fitTMonthlyGMean
;
;plot, residualsTYearlyGMean
;oplot, fitTYearlyGMean

END