#Two layer version of the seasonal cycle model,
#tracking surface temperature separately from air temperature.
#It is especially useful for illustrating the seasonal
#cycle of the Southern hemisphere
#
#Modified version to compute seasonal cycle variation over
#Milankovic cycles, and over CO2 cycles. (Move this to NZ glacier
#research directory once the solar module is put in a publically
#accessible place)

#Milankovic data
datapath ='/Users/rtp1/Havsornen/Teaching/GeoSci232/WorkbookDatasets/'
dataset = 'Chapter7Data/Milankovic/orbit91'




from math import *
from solar import * #This module is currently in the local chapter directory
from ClimateUtilities import *
import phys

#Get the orbital parameters data
orbitParams = readTable(datapath+dataset)

L = 1370.
year = 365.*24.*3600.
degrad = pi/180.

#Ocean thermal inertia parameters
rho = 1000.
H = 100.
Cp = 4218. #Specific heat for liquid water

#Air thermal inertia parameters

#Atmospheric flux and radiative parameters
swAbsCoeff = .2 #Shortwave solar absorption
aTurb = 10. #Turbulent exchange coeff, W/m**2/K
aRad = 4.*phys.sigma*(280.)**3 # Radiative exchange coeff
E0 = 30.  #Base state evaporation
Ma = 1.e4 #Atmospheric mass in kg
Cpa = phys.air.cp
#
lat = -45.
albedo = .23
#
#Function to look up and return Earth's orbital parameters
#The argument is time before present, in thousands of years.
#The function returns a triple consisting of eccentricity,
#obliquity (degrees) and the angle between the position of
#the Northern Hemisphere Summer Solstice and the position of
#the perihelion.
#
def getOrbit(time):
    e = orbitParams['ECC'][time]
    obl = orbitParams['OBL'][time]
    prec = orbitParams['OMEGA'][time] + 90. #relative to NH solstice
    if prec > 360.:
        prec -= 360.
    return e,obl,prec

#ToDo: Make OLR a parameter, and allow for
#use of polynomial fit.  Provide a module containing
#several OLR fits as part of radiation chapter.
#def OLR(T):
#    return phys.sigma*T**4
def OLR(T):
    a,b = (172.22409548794681, 2.232345632274972)
    return a + b*(T-250.)

def estar(T):
    return .3 #Replace with polynomial fit later

#T is a vector of temperatures. T[0] is the surface
#temperature and T[1] is the low level atmosphere temperature
#ToDo: Put in effect of atmospheric transparency on OLR. The
#OLR here is the correct OLR in the opaque limit
def f(t,T,params):
    #Find solar flux for the current time of year
    #**ToDo: For extreme eccentricities, need to compute
    #length of year changes properly.
    tt = orb['t'][-1]*t/year
    while tt > orb['t'][-1]:
        tt -= orb['t'][-1]
    orbitTime = Numeric.searchsorted(orb['t'],tt)
    fluxFactor = flux[orbitTime]
    #
    S = (1.-albedo)*L*fluxFactor
    surfHeat = (1.-swAbsCoeff)*S - estar(T[1])*phys.sigma*T[1]**4 \
               - (aTurb+aRad)*(T[0]-T[1]) - E0
    atmosHeat = swAbsCoeff*S + (aTurb+aRad)*(T[0]-T[1]) \
                + estar(T[1])*phys.sigma*T[1]**4 - OLR(T[1]) + E0
    return Numeric.array([surfHeat/(rho*H*Cp),atmosHeat/(Ma*Cpa)])

#Epoch is time B.P. in millennia
def Tcomp(epoch):
    global orb,flux
    #Get Milankovic fluxes
    e,obl,precess = getOrbit(epoch)
    orb  = orbit(e,precess*degrad)#Start at the solstice
    flux = orbitFlux(lat*degrad,obl*degrad,precess*degrad,orb)

    dt = year/100.
    fi = integrator(f,0.,Numeric.array([300.,300.]),dt) # Initial temperature is 300K
    time = []
    Tlist = []
    #First do a 50 year integration to reach equilibrium
    n = int(50*year/dt)

    for i in range(n):
        fi.next()
        
    Ts = []
    Ta = []
    #Now continue for 1 year and save the results
    for i in range(int(year/dt)):
        result = fi.next()
        time.append(result[0]/year - 50.)
        Ts.append(result[1][0])
        Ta.append(result[1][1])

    c = Curve()
    c.addCurve(time,'t','Time (years)')
    c.addCurve(Ts,'Ts','Surface Temperature (Kelvin)')
    c.addCurve(Ta,'Ta','Air temperature (Kelvin)')
    return c

c = Curve()
for epoch in range(0,400,5):
    c1 = Tcomp(epoch)
    if epoch == 0:
        c.addCurve(c1['t'],'time')
    c.addCurve(c1['Ta'],'Ta%d'%epoch)
plot(c)


#Compute degree days
names = c.listVariables()[1:]
TddList = []
c3 = Curve()
c3.addCurve(range(0,400,5),'TimeBP')
T0 = 282.
for name in names:
	T = c[name]
	n = sum(Numeric.where(T>T0,T-T0,0.))
	TddList.append(n)

	
c3.addCurve(TddList)
c3['v1'] = c3['v1']*365./100.
plot(c3)