;
PRO waccmx_solar_f107_ap_monthly_yearly_1944, f107AllMonthly, apAllMonthly, f107YearlyMean, apYearlyMean
;
;Calculate WACCM-X 1944-1970 f10.7 and ap monthly and yearly means from daily values
;
;
;  Open f10.7 WACCM-X input file with daily values
;
;fid = ncdf_open("/glade/p/cesmdata/inputdata/atm/waccm/solar/wasolar_1948-2140_daily_c100805.nc")
fid = ncdf_open("/hao/aim/joemci/WACCMX/inputdata/proxy_solar_Solomon_Richmond_1845-2008_daily_c100423.nc")
;
;  Get daily f10.7 values
;
ncdf_varget, fid, 'f107', f107In
;
;  Get daily ap values
;
ncdf_varget, fid, 'ap', apIn
;
;Loop over years and months and calculate monthly f10.7 values
;
nYears = 27
;nYears = 51
nMonths = 12
iMonthAll = 0
iDayStartIn = 36342 ; This comes from looking at f107/ap file and getting index of Jan. 1, 1944
iLeap = 0
;iLeap = 2

f107AllMonthly = FLTARR(nYears*nMonths)
apAllMonthly = FLTARR(nYears*nMonths)

f107Monthly = FLTARR(nMonths)
apMonthly = FLTARR(nMonths)

f107YearlyMean = FLTARR(nYears)
apYearlyMean = FLTARR(nYears)

FOR iYear=0,nYears-1 DO BEGIN
  
  iLeap = iLeap + 1
  IF iLeap EQ 5 THEN iLeap = 1

  IF iYear EQ 0 THEN iDayIn = iDayStartIn - 1

  f107Monthly[*] = 0.
  apMonthly[*] = 0.

  FOR iMonth=0,nMonths-1 DO BEGIN

    nDaysOfMonth = 31

    IF iMonth EQ 1 THEN nDaysOfMonth = 28
    IF (iMonth EQ 1 AND iLeap EQ 4) THEN nDaysOfMonth = 29
    IF (iMonth EQ 3 OR iMonth EQ 5 OR iMonth EQ 8 OR iMonth EQ 10) THEN nDaysOfMonth = 30
    
    f107AllMonthly[iMonthAll] = 0.
    apAllMonthly[iMonthAll] = 0.

    FOR iDay = 0, nDaysOfMonth-1 DO BEGIN

      iDayIn = IDayIn + 1
      
      f107AllMonthly[iMonthAll] = f107AllMonthly[iMonthAll] + f107In[iDayIn]     
      apAllMonthly[iMonthAll] = apAllMonthly[iMonthAll] + apIn[iDayIn] 
      
      f107Monthly[iMonth] =  f107Monthly[iMonth] + f107In[iDayIn]   
      apMonthly[iMonth] =  apMonthly[iMonth] + apIn[iDayIn]   

    ENDFOR
    
    f107AllMonthly[iMonthAll] = f107AllMonthly[iMonthAll] / nDaysOfMonth
    apAllMonthly[iMonthAll] = apAllMonthly[iMonthAll] / nDaysOfMonth

    f107Monthly[iMonth] =  f107Monthly[iMonth] / nDaysOfMonth   
    apMonthly[iMonth] =  apMonthly[iMonth] / nDaysOfMonth   
     
    iMonthAll = iMonthAll + 1

  ENDFOR

  f107YearlyMean[iYear] = MEAN(f107Monthly[*])
  apYearlyMean[iYear] = MEAN(apMonthly[*])

ENDFOR

SET_PLOT, 'PS'

psfile = 'monthly_yearly_f10.7_ap_WACCMX_1944.ps'

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

MonthsAll = INDGEN(iMonthAll)

yTitle = 'f107'
xTitle = 'Months Since January, 1944'

Title = 'f10.7 1944-1970'
plot, MonthsAll, f107AllMonthly[0:iMonthAll-1],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

yTitle = 'ap'
Title = 'ap 1944-1970'
plot, MonthsAll, apAllMonthly[0:iMonthAll-1],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

Years = INDGEN(nYears)

yTitle = 'f107'
xTitle = 'Years Since 1944'

Title = 'f10.7 1944-1970'
plot, Years, f107YearlyMean[0:nYears-1],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

yTitle = 'ap'
Title = 'ap 1944-1970'
plot, Years, apYearlyMean[0:nYears-1],title=Title,xtitle=xTitle,ytitle=yTitle,xstyle=1,ystyle=1

device,/close

RETURN

END