;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
;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 = 51
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 51 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 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 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
;
;  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 = FLTARR(3,nYears-20)
f107apYearlyMean[0,*] = f107YearlyMean[15:nYears-6]
f107apYearlyMean[1,*] = apYearlyMean[15:nYears-6]
f107apYearlyMean[2,*] = f107YearlyMeanSquared[15:nYears-6]

TYearlyGMean = TYearlyGlobalMean[15:nYears-6,7]
regressCoeffs = REGRESS(f107apYearlyMean, TYearlyGMean, YFIT=TYearlyGMeanFit)
residualsTYearlyGMean = TYearlyGMean - TYearlyGMeanFit
;
;  Get fit to residuals
;
Years = INDGEN(nYears-20)
fitCoeffs = POLY_FIT(Years,residualsTYearlyGMean,1,YFIT=fitTYearlyGMean) 
;
;  Now do regression of monthly mean temperature with f10.7 and ap
;
;f107apMonthlyMean = FLTARR(3,iMonthAll)
;f107apMonthlyMean[0,*] = f107AllMonthly[0:iMonthAll-1]
;f107apMonthlyMean[1,*] = apAllMonthly[0:iMonthAll-1]
;f107apMonthlyMean[2,*] = f107AllMonthlySquared[0:iMonthAll-1]
f107apMonthlyMean = FLTARR(3,iMonthAll-180)
f107apMonthlyMean[0,*] = f107AllMonthly[180:iMonthAll-1]
f107apMonthlyMean[1,*] = apAllMonthly[180:iMonthAll-1]
f107apMonthlyMean[2,*] = f107AllMonthlySquared[180:iMonthAll-1]

;TMonthlyAllGMean = TMonthlyAllGlobalMean[*,7]
TMonthlyAllGMean = TMonthlyAllGlobalMean[180:iMonthAll-1,7]
regressCoeffs = REGRESS(f107apMonthlyMean, TMonthlyAllGMean, YFIT=TMonthlyGMeanFit)
residualsTMonthlyGMean = TMonthlyAllGMean - TMonthlyGMeanFit
;
;  Get fit to residuals
;
;Months = INDGEN(iMonthAll)
;fitCoeffs = POLY_FIT(Months,residualsTMonthlyGMean,1,YFIT=fitTMonthlyGMean) 

Months = INDGEN(iMonthAll-180)
fitCoeffs = POLY_FIT(Months,residualsTMonthlyGMean[0:iMonthAll-181],1,YFIT=fitTMonthlyGMean) 

plot, residualsTMonthlyGMean
oplot, fitTMonthlyGMean

plot, residualsTYearlyGMean
oplot, fitTYearlyGMean

;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_WACCMX.ps'

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

MonthsAll = INDGEN(iMonthAll)

yTitle = 'CO2 mol/mol'
xTitle = 'Months Since January, 1955'

Title = 'CO2 1955-2005 Area Weighted Global Mean ~1e-07hPa ~280km'
plot, MonthsAll, CO2MonthlyAllGlobalMean[0:iMonthAll-1,7],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'CO2 1955-2005 Area Weighted Global Mean ~0.018hPa ~75km'
plot, MonthsAll, CO2MonthlyAllGlobalMean[0:iMonthAll-1,31],yrange=[0.0003,0.00038],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'CO2 1955-2005 Area Weighted Global Mean ~8.6hPa  ~32km'
plot, MonthsAll, CO2MonthlyAllGlobalMean[0:iMonthAll-1,50],yrange=[0.0003,0.00038],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'CO2 1955-2005 Area Weighted Global Mean Surface'
plot, MonthsAll, CO2MonthlyAllGlobalMean[0:iMonthAll-1,80],yrange=[0.0003,0.00038],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

yTitle = 'T (K)'

Title = 'T 1955-2005 Area Weighted Global Mean ~1e-07hPa ~280km'
plot, MonthsAll, TMonthlyAllGlobalMean[0:iMonthAll-1,7],yrange=[500.0,1300.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'T 1955-2005 Area Weighted Global Mean ~0.018hPa ~75km'
plot, MonthsAll, TMonthlyAllGlobalMean[0:iMonthAll-1,31],yrange=[190.0,220.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'T 1955-2005 Area Weighted Global Mean ~8.6hPa  ~32km'
plot, MonthsAll, TMonthlyAllGlobalMean[0:iMonthAll-1,50],yrange=[222.0,238.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'T 1955-2005 Area Weighted Global Mean Surface'
plot, MonthsAll, TMonthlyAllGlobalMean[0:iMonthAll-1,80],yrange=[285.0,300.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

device,/close

psfile = 'yearly_CO2_T_AreaGlobalMean_WACCMX.ps'

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

YearsAll = INDGEN(nYears)

yTitle = 'CO2 mol/mol'
xTitle = 'Years Since 1955'

Title = 'CO2 1955-2005 Area Weighted Global Mean ~1e-07hPa ~280km'
plot, YearsAll, CO2YearlyGlobalMean[0:nYears-1,7],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'CO2 1955-2005 Area Weighted Global Mean ~0.018hPa ~75km'
plot, YearsAll, CO2YearlyGlobalMean[0:nYears-1,31],yrange=[0.0003,0.00038],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'CO2 1955-2005 Area Weighted Global Mean ~8.6hPa  ~32km'
plot, YearsAll, CO2YearlyGlobalMean[0:nYears-1,50],yrange=[0.0003,0.00038],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'CO2 1955-2005 Area Weighted Global Mean Surface'
plot, YearsAll, CO2YearlyGlobalMean[0:nYears-1,80],yrange=[0.0003,0.00038],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

yTitle = 'T (K)'

Title = 'T 1955-2005 Area Weighted Global Mean ~1e-07hPa ~280km'
plot, YearsAll, TYearlyGlobalMean[0:nYears-1,7],yrange=[500.0,1300.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'T 1955-2005 Area Weighted Global Mean ~0.018hPa ~75km'
plot, YearsAll, TYearlyGlobalMean[0:nYears-1,31],yrange=[190.0,220.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'T 1955-2005 Area Weighted Global Mean ~8.6hPa  ~32km'
plot, YearsAll, TYearlyGlobalMean[0:nYears-1,50],yrange=[222.0,238.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'T 1955-2005 Area Weighted Global Mean Surface'
plot, YearsAll, TYearlyGlobalMean[0:nYears-1,80],yrange=[285.0,300.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

device,/close

END



  
rad    = 4.0*atan(1.0)/180.0 ;pi/180.
earthRadius	= 6371220.0 
rr     = re*rad
                                                    
dLat(0)  = abs( lat(2) - lat(1) ) * rr
dLat(1:nlat-2)  = abs( lat(2:nLat-1) - lat(0:nLat-3) ) * rr * 0.5   
dLat(nLat-1)	= abs( lat(nLat-1) - lat(nLat-2 )) * rr

;  areaWeight   = wtLat * wtLon                                ; cell area function of latitude only

ncdf_varget, fid, 'lev', pres

FOR iYear=0,nYears-1 DO BEGIN

    cYear = STRTRIM(iYear+1955,2)
  print, 'Doing year ', cYear
 
  FOR iMonth=0,nMonths-1 DO BEGIN
    ;
    ;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]
    
    CO2All[iMonthAll,0] = CO2In[72,48,7]  ;Thermosphere near 1e-07hPa ~280km
    CO2All[iMonthAll,1] = CO2In[72,48,31] ;Mesosphere near 0.018hPa ~75km
    CO2All[iMonthAll,2] = CO2In[72,48,50] ;Stratosphere near 8.6hPa  ~32km
    CO2All[iMonthAll,3] = CO2In[72,48,80] ;Surface

 q_ave = SUM[q*wgty*wgtx]/SUM[wgty*wgtx]
 
    FOR iLat=0,nLat-1 DO BEGIN
 
      FOR iLon=0,nLon-1 DO BEGIN
 
      CO2WeightedZonalMean[iLat] = MEAN(CO2In[*,iLat,7]) * dLat[iLat]
 
      ENDFOR
   
    ENDFOR

      CO2GlobalMean[iMonth,7] = SUM(CO2WeightedZonalMean) / SUM(dLat)
    
    CO2AllMean[iMonthAll,0] = SUM(CO2ZonalMean[*] * dLat[*])  ;Thermosphere near 1e-07hPa ~280km

    CO2ZonalMean[iMonthAll,0] = MEAN(CO2In[*,iLat,7])  ;Thermosphere near 1e-07hPa ~280km
    CO2ZonalMean[iMonthAll,1] = MEAN(CO2In[*,iLat,31]) ;Mesosphere near 0.018hPa ~75km
    CO2AllMean[iMonthAll,2] = MEAN(CO2In[*,iLat,50]) ;Stratosphere near 8.6hPa  ~32km
    CO2AllMean[iMonthAll,3] = MEAN(CO2In[*,iLat,80]) ;Surface

    ncdf_varget, fid, 'T', TIn
    
    TAll[iMonthAll,0] = TIn[72,48,7]  ;Thermosphere near 1e-07hPa ~280km
    TAll[iMonthAll,1] = TIn[72,48,31] ;Mesosphere near 0.018hPa ~75km
    TAll[iMonthAll,2] = TIn[72,48,50] ;Stratosphere near 8.6hPa  ~32km
    TAll[iMonthAll,3] = TIn[72,48,80] ;Surface
    
    TAllMean[iMonthAll,0] = MEAN(TIn[*,*,7])  ;Thermosphere near 1e-07hPa ~280km
    TAllMean[iMonthAll,1] = MEAN(TIn[*,*,31]) ;Mesosphere near 0.018hPa ~75km
    TAllMean[iMonthAll,2] = MEAN(TIn[*,*,50]) ;Stratosphere near 8.6hPa  ~32km
    TAllMean[iMonthAll,3] = MEAN(TIn[*,*,80]) ;Surface

    ncdf_varget, fid, 'U', UIn
    
    UZonalMean[iMonthAll,0] = MEAN(UIn[*,48,7])  ;Thermosphere near 1e-07hPa ~280km
    UZonalMean[iMonthAll,1] = MEAN(UIn[*,48,31]) ;Mesosphere near 0.018hPa ~75km
    UZonalMean[iMonthAll,2] = MEAN(UIn[*,48,50]) ;Stratosphere near 8.6hPa  ~32km
    UZonalMean[iMonthAll,3] = MEAN(UIn[*,48,80]) ;Surface
    
    rhoAllMean[iMonthAll,0] = pres / R / TAll
;
;  Get lowest top geometric altitude for whole 51 years to use in choosing of vertical levels to do real global mean
;    
    ncdf_varget, fid, 'Z3', geopAlt

    geomAlt = geopAlt * (1. + geopAlt/earthRadius) /1000.    

    geomAltAllMeanTop = MEAN(geomAlt[*,*,1])    
    
;    print, 'Geometric altitude at top level ', geomAltAllMeanTop
    if geomAltAllMeanTop lt 340. then print, 'Geopotential altitude at top level below 340', geomAltAllMeanTop
    
    iMonthAll = iMonthAll + 1

  ENDFOR

ENDFOR

SET_PLOT, 'PS'

psfile = 'monthly_CO2_T_Eq_180Lon_WACCMX.ps'

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

MonthsAll = INDGEN(iMonthAll)

yTitle = 'CO2 mol/mol'
xTitle = 'Months Since January, 1955'

Title = 'CO2 1955-2005 Equator 180 Degrees Longitude ~1e-07hPa ~280km'
plot, MonthsAll, CO2AllEq[0:iMonthAll-1,0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'CO2 1955-2005 Equator 180 Degrees Longitude ~0.018hPa ~75km'
plot, MonthsAll, CO2AllEq[0:iMonthAll-1,1],yrange=[0.0003,0.00038],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'CO2 1955-2005 Equator 180 Degrees Longitude ~8.6hPa  ~32km'
plot, MonthsAll, CO2AllEq[0:iMonthAll-1,2],yrange=[0.0003,0.00038],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'CO2 1955-2005 Equator 180 Degrees Longitude Surface'
plot, MonthsAll, CO2AllEq[0:iMonthAll-1,3],yrange=[0.0003,0.00038],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

yTitle = 'T (K)'

Title = 'T 1955-2005 Equator 180 Degrees Longitude ~1e-07hPa ~280km'
plot, MonthsAll, TAllEq[0:iMonthAll-1,0],yrange=[500.0,1300.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'T 1955-2005 Equator 180 Degrees Longitude ~0.018hPa ~75km'
plot, MonthsAll, TAllEq[0:iMonthAll-1,1],yrange=[180.0,220.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'T 1955-2005 Equator 180 Degrees Longitude ~8.6hPa  ~32km'
plot, MonthsAll, TAllEq[0:iMonthAll-1,2],yrange=[226.0,238.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'T 1955-2005 Equator 180 Degrees Longitude Surface'
plot, MonthsAll, TAllEq[0:iMonthAll-1,3],yrange=[296.0,302.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

device,/close

psfile = 'monthly_CO2_T_AreaGlobalMean_WACCMX.ps'

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

MonthsAll = INDGEN(iMonthAll)

yTitle = 'CO2 mol/mol'
xTitle = 'Months Since January, 1955'

Title = 'CO2 1955-2005 Global Mean ~1e-07hPa ~280km'
plot, MonthsAll, CO2AllEqMean[0:iMonthAll-1,0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'CO2 1955-2005 Global Mean ~0.018hPa ~75km'
plot, MonthsAll, CO2AllEqMean[0:iMonthAll-1,1],yrange=[0.0003,0.00038],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'CO2 1955-2005 Global Mean ~8.6hPa  ~32km'
plot, MonthsAll, CO2AllEqMean[0:iMonthAll-1,2],yrange=[0.0003,0.00038],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'CO2 1955-2005 Global Mean Surface'
plot, MonthsAll, CO2AllEqMean[0:iMonthAll-1,3],yrange=[0.0003,0.00038],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

yTitle = 'T (K)'

Title = 'T 1955-2005 Global Mean ~1e-07hPa ~280km'
plot, MonthsAll, TAllEqMean[0:iMonthAll-1,0],yrange=[500.0,1300.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'T 1955-2005 Global Mean ~0.018hPa ~75km'
plot, MonthsAll, TAllEqMean[0:iMonthAll-1,1],yrange=[190.0,220.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'T 1955-2005 Global Mean ~8.6hPa  ~32km'
plot, MonthsAll, TAllEqMean[0:iMonthAll-1,2],yrange=[222.0,238.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'T 1955-2005 Global Mean Surface'
plot, MonthsAll, TAllEqMean[0:iMonthAll-1,3],yrange=[270.0,290.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

device,/close

psfile = 'monthly_U_GlobalMean_WACCMX.ps'

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

MonthsAll = INDGEN(iMonthAll)

yTitle = 'U (m/s)'
xTitle = 'Months Since January, 2001'

Title = 'U 2001-2005 Zonal Mean At Equator ~1e-07hPa ~280km'
plot, MonthsAll, UEqZonalMean[552:iMonthAll-1,0],xrange=[0,60],yrange=[20.0,60.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'U 2001-2005 Zonal Mean At Equator ~0.018hPa ~75km'
plot, MonthsAll, UEqZonalMean[552:iMonthAll-1,1],xrange=[0,60],yrange=[-20.0,50.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'U 2001-2005 Zonal Mean At Equator ~8.6hPa  ~32km'
plot, MonthsAll, UEqZonalMean[552:iMonthAll-1,2],xrange=[0,60],yrange=[-30.0,25.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'U 2001-2005 Zonal Mean At Equator Surface'
plot, MonthsAll, UEqZonalMean[552:iMonthAll-1,3],xrange=[0,60],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

device,/close

psfile = 'monthly_U_GlobalMean_WACCMX_1955-2005.ps'

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

MonthsAll = INDGEN(iMonthAll)

yTitle = 'U (m/s)'
xTitle = 'Months Since January, 1955'

Title = 'U 1955-2005 Zonal Mean At Equator ~1e-07hPa ~280km'
plot, MonthsAll, UEqZonalMean[0:iMonthAll-1,0],yrange=[20.0,60.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'U 1955-2005 Zonal Mean At Equator ~0.018hPa ~75km'
plot, MonthsAll, UEqZonalMean[0:iMonthAll-1,1],yrange=[-20.0,50.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'U 1955-2005 Zonal Mean At Equator ~8.6hPa  ~32km'
plot, MonthsAll, UEqZonalMean[0:iMonthAll-1,2],yrange=[-30.0,25.0],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1
Title = 'U 1955-2005 Zonal Mean At Equator Surface'
plot, MonthsAll, UEqZonalMean[0:iMonthAll-1,3],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

device,/close

;plot,months,REFORM(t12spc_mean[*,73,34,ilev]),xrange=[-1,12],ytitle=yTitle,yrange=[0,1.5],xticks=13,xtickname = CMonths,thick=4,_EXTRA=gang_plot_pos(5,1,3) 

END