#This is a data conversion utility script which extracts
#monthly mean surface temperature data for the ERBE years
#from the CDC-derived netcdf monthly-mean version of the
#NCEP reanalysis data product. It writes out the lat-lon
#data into a series of netcdf files (one per year for
#the ERBE years) and a corresponding set of plain text
#zonal mean files.  This surface temperature data is
#intended for use with the ERBE radiation budget data,
#for example for exercises comparing the observed OLR
#with the upward infrared emitted by the surface.

#The CDC data was obtained from http://**.

#****UNDER CONSTRUCTION*****

#ToDo: Provide netCDF datasets of upper air temperature and humidity
#data as well, for more sophisticated radiation exercises in next
#chapter.

#**ToDo: Interpolate to the same lat-lon grid as the ERBE data

import cdms,MV,Numeric
from ClimateUtilities import *

#Location of the data
datapath = '/Users/rtp1/bigdata/ncep_mm_surface/'
surfTfilename = 'air.mon.mean.nc'

#-------Define dictionaries and constants--------------------
months = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sept',
          'Oct','Nov','Dec']
ERBEyears = [1985,1986,1987,1988,1989]

f = cdms.open(datapath+surfTfilename)
yearHours = 365.24*24. #ToDo: Replace this with exact value

T = f['air']
t = 1. + T.getTime()[:]/yearHours #Time in decimal years
lat = T.getLatitude()[:]

for year in ERBEyears:
    start = Numeric.searchsorted(t,year)
    print start,year
    c = Curve()
    #A temporary hack to interpolate to the ERBE grid
    latE = .5*(lat[0:-1]+lat[1:])
    c.addCurve(latE,'lat')
    T1 = T[start:start+12] + 273.15 #Extract and convert to K
    Tbar = MV.average(T1,axis = 2) #Zonal mean
    Tann = MV.average(Tbar) # Annual average
    #Temporary interpolation hack
    TbarE = .5*(Tbar[:,0:-1]+Tbar[:,1:])
    TannE = MV.average(TbarE)
    for month in range(12):
        c.addCurve(TbarE[month],months[month])
    outfile = 'TsBar%d.txt'%year
    c.dump(outfile)