"""This script computes the adiabatic temperature profile
for a single-component condensible atmosphere.  The
atmosphere may be subsaturated at the ground, in which
case, the profile follows the dry adiabat until it reaches
saturation, whereafter it follows the saturation
temperature curve given by the Clausius-Clapeyron relation
"""
#Data on section of text which this script is associated with
Chapter = '2.**'
Figure = '**'


import math,phys
from ClimateUtilities import *



#Set needed thermodynamic constants
stuff = phys.CO2
L = stuff.L_sublimation # Latent heat
cp = stuff.cp #Specific heat
R = stuff.R #Gas constant for the atmosphere
Tr = stuff.TriplePointT # Reference temperature
psat_ref = stuff.TriplePointP  # Saturation vapor pressure at
                               #reference temperature

#Function that solves for temperature at which saturation
#occurs, for a given pressure.  This uses the simplified
#form of Clausius-Clapeyron based on the assumption that
#the latent heat is constant.
def Tsat(p):
    return(1./(1./Tr - (R/L)*math.log(p/psat_ref)))

#Function giving the temperature of the dry adiabat
#p0 and T0 are the surface pressure and temperature
def Tdry(p,p0,T0):
    return(T0*(p/p0)**(R/cp))

#
#Function returning a list consisting of pressure
#and adiabatic temperature, ready for plotting
#p0 and T0 are surface pressure and temperature.
#N is the number of (logarithmically spaced)
#points to plot.
#
def MoistAdiabat1comp(p0,T0,N):
    p1 = .001*p0
    pfact = (p1/p0)**(1./N)
    # Make lists to hold the results
    plist = []
    Tlist = []
    # Do the computation
    p = p0
    while p >= p1:
        plist.append(p)
        Tlist.append(max(Tsat(p),Tdry(p,p0,T0)))
        p = p*pfact
    # Build the results into a curve object and return
    c = Curve()
    c.addCurve(plist,'p','pressure')
    c.addCurve(Tlist,'T','Temperature')
    #
    #Set plot options
    c.switchXY = c.reverseY = c.YlogAxis = True
    c.PlotTitle = '%s Condensing adiabat, unsaturated sfc'%stuff.formula
    c.Xlabel = 'Pressure (Pa)'
    c.Ylabel = 'Temperature (K)'
    return c

# Main script: Make a plot
result = MoistAdiabat1comp(425.,250.,100)
#plot(result)
plot(result)
result.dump()