;waccmx_1944_yearly_pres_savout.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.e13b02.F_1955_2005_WACCMX_CN.f19_f19.003_py/f.e13b02.F_1955_2005_WACCMX_CN.f19_f19.003.cam.h0.1944-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 = 19
yearOne = 1944
nMonths = 12
iMonthAll = 0

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

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

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

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

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

    CO2WeightedMean[*] = 0.
    TWeightedMean[*] = 0.
    geomAltWeightedMean[*] = 0.
    NOWeightedMean[*] = 0.
    TREFHTWeightedMean = 0.
 
    ;
    ;Open file, read monthly file CO2 and T, close file
    ;
    cMonth = STRTRIM(iMonth+1,2)
    cYear = STRTRIM(iYear+1944,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.e13b02.F_1955_2005_WACCMX_CN.f19_f19.003_py/f.e13b02.F_1955_2005_WACCMX_CN.f19_f19.003.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
    ncdf_varget, fid, 'NO', NOIn
    ncdf_varget, fid, 'TREFHT', TREFHTIn

    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]
	  NOWeightedMean[iPres] = NOWeightedMean[iPres] + NOIn[ILon,iLat,iPres] * areaWeight[iLat]
	  IF iPres EQ 0 THEN TREFHTWeightedMean = TREFHTWeightedMean + TREFHTIn[ILon,iLat] * areaWeight[iLat]
	  
        ENDFOR ;Lon

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

    ENDFOR ;Pres
  
    TREFHTMonthlyGlobalMean[iMonth] = TREFHTWeightedMean / sumAreaWeight
    
    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]    
    
    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])
    NOYearlyGlobalMean[iYear,iPres] = MEAN(NOMonthlyGlobalMean[*,iPres])

  ENDFOR ;Pres
 
  TREFHTYearlyGlobalMean[iYear] = MEAN(TREFHTMonthlyGlobalMean[*])

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

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

ENDFOR ;Year

save,iMonthAll,nPres,pres,CO2MonthlyAllGlobalMean,TMonthlyAllGlobalMean,geomAltMonthlyAllGlobalMean,NOMonthlyAllGlobalMean,TREFHTMonthlyAllGlobalMean, file ='/glade/scratch/joemci/TrendSavFiles/CO2_T_Monthly_GlobalMean_Pres_1944_1947.sav'
save, yearOne,nYears,nPres,pres,CO2YearlyGlobalMean,TYearlyGlobalMean,geomAltYearlyGlobalMean,NOYearlyGlobalMean,TREFHTYearlyGlobalMean, file ='/glade/scratch/joemci/TrendSavFiles/CO2_T_Yearly_GlobalMean_Pres_1944_1947.sav'

END