;waccmx_50yr_monthly_levs.pro
;
;Plot WACCM-X 1955-2005 CO2 and Temperature at 4 levels near the surface, in the stratosphere, in the mesosphere and in the thermosphere. 
;Use h0 monthly data to plot values near the equator and 180 degrees longitude (Mauna Loa)
;
;
;Loop over years and months and read CO2 and temperature to complile 50 year monthly values
;
nYears = 51
nMonths = 12
iMonthAll = 0

CO2AllEq = FLTARR(nYears*nMonths,4)
TAllEq = FLTARR(nYears*nMonths,4)

CO2AllEqMean = FLTARR(nYears*nMonths,4)
TAllEqMean = FLTARR(nYears*nMonths,4)
UEqZonalMean = FLTARR(nYears*nMonths,4)

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
    
    CO2AllEq[iMonthAll,0] = CO2In[72,48,7]  ;Thermosphere near 1e-07hPa ~280km
    CO2AllEq[iMonthAll,1] = CO2In[72,48,31] ;Mesosphere near 0.018hPa ~75km
    CO2AllEq[iMonthAll,2] = CO2In[72,48,50] ;Stratosphere near 8.6hPa  ~32km
    CO2AllEq[iMonthAll,3] = CO2In[72,48,80] ;Surface
    
    CO2AllEqMean[iMonthAll,0] = MEAN(CO2In[*,*,7])  ;Thermosphere near 1e-07hPa ~280km
    CO2AllEqMean[iMonthAll,1] = MEAN(CO2In[*,*,31]) ;Mesosphere near 0.018hPa ~75km
    CO2AllEqMean[iMonthAll,2] = MEAN(CO2In[*,*,50]) ;Stratosphere near 8.6hPa  ~32km
    CO2AllEqMean[iMonthAll,3] = MEAN(CO2In[*,*,80]) ;Surface

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

    ncdf_varget, fid, 'U', UIn
    
    UEqZonalMean[iMonthAll,0] = MEAN(UIn[*,48,7])  ;Thermosphere near 1e-07hPa ~280km
    UEqZonalMean[iMonthAll,1] = MEAN(UIn[*,48,31]) ;Mesosphere near 0.018hPa ~75km
    UEqZonalMean[iMonthAll,2] = MEAN(UIn[*,48,50]) ;Stratosphere near 8.6hPa  ~32km
    UEqZonalMean[iMonthAll,3] = MEAN(UIn[*,48,80]) ;Surface
;
;  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/6370000.0) /1000.    

    geomAltAllMeanTop = MEAN(geomAlt[*,*,1])
    
    
;    print, 'Geometric altitude at top level ', geomAltAllMeanTop
    if geomAltAllMeanTop lt 340. then print, 'Geopotential altitude at top level below 330', 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_GlobalMean_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