#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

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

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

#Heat flux added to atmosphere by dynamic flux
Fh = 70. #W/m**2

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

#Air thermal inertia parameters


#Atmospheric flux and radiative parameters
swAbsCoeff = .2 #Shortwave solar absorption
aTurb = 20. #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

#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):
    if T<273.15:
        albedo = .6
        H = 1.
    else:
        albedo = .1
        H = 50.
    S = (1.-albedo)*L*solar(lat,2.*pi*t/year,obliquity)
    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 \
                + Fh
    return Numeric.array([surfHeat/(rho*H*Cp),atmosHeat/(Ma*Cpa)])

dt = year/100.
c = Curve()

lat = 80.*degrad
albedo = .1
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.addCurve(time,'t','Time (years)')
c.addCurve(Ts,'Ts','Surface Temperature (Kelvin)')
c.addCurve(Ta,'Ta','Air temperature (Kelvin)')
plot(c)