# Atmospheric thermodynamics.
# Functions first, then examples.
# This version uses a very approximate form of the potential temperature,
# ThetaE = Theta + 3000*w where w is the water vapor mixing ratio in kg/kg and Theta is potential temperature.  
# Put in values of Po, To, and Tdew in lines 141-143 for the LCL calculation. 

import math
# Global variables.
eo = 6.112      # partial pressure of water vapor at saturation at 0 C.
Ro = 8314.3     # Universal gas constant in Joules/K kmole
RD = Ro/28.97   # Dry air gas constant in Joules/K kg
Rv = Ro/18.02   # Water vapor gas constant in Joules/K kg
epsilon=RD/Rv   # Ratio of molecular mass for H20 and dry air.
cp = 7.0*RD/2.0 # Specific heat at constant pressure.

# Function for saturation vapor pressure over water in mb for T in C.
# From Bolton.
def esatWater(T):
    return eo*math.exp(17.67*T/(243.5+T))
    
# Function to calculate the latent heat of water from T in C,
# From https://en.wikipedia.org/wiki/Latent_heat.
def LatentHeatWater(Tc):
    Llv=(2500.8-2.36*Tc+0.0016*Tc*Tc-0.00006*Tc*Tc*Tc)*1000.0
    return Llv # J/kg.
    
# Function for RH from T and Tdew in C.
# Calls: esatWater function.  
def RelativeHumidity(Tdew,T):
    return 100.0*esatWater(Tdew)/esatWater(T)

# Water vapor mixing ratio in g/kg
# Inputs: water vapor pressure, e, in mb and total pressure, P, in mb.
def mixingRatio(e,P):
    return epsilon*e*1000./P

# Potential temperature calculator, returning in Celsius.
# T is entered in Celsius, P in mb.
def PotTemp(T,P):
    theta=(T+273.15)*((1000.0/P)**(2.0/7.0))
    return theta-273.15
    
# Pressure altitude in meters.
# Inputs:
# To, Po, starting temperature in Kelvin and mb.
# Pz pressure in mb where the pressure altitude is desired.
def PressureAltitude(To,Po,Pz):
    gamma=9.8/1000.0 # Lapse rate in Kelvin / meter
    Ro = 8314.3 # Universal gas constant in Joules/K kmole
    RD = Ro/28.97 # Dry air gas constant in Joules/K kg
    g  = 9.81 # Acceleration due to gravity in N/kg.
    expon=RD*gamma/g
    PressureRatio=Pz/Po
    z=(To/gamma)*(1-PressureRatio**expon)
    return z
    
# Get the dewpoint temperature from ThetaE.
# Less accurate than from the first law.
# INPUTS:
#     Wetbulb and air temperature in C, Tw and To
#     pressure in mb, Po.
# OUTPUTS:
#     Tdew in Celsius.
#     Relative humidity in %.
# Theory: 
#     Equivalent potential temperature can be obtained from either
#     Tw and Po, or from Tdew, To, Po.  Solve for Tdew.
def getTdewFromThetaE(Tw,To,Po):
    R=Po*((PotTemp(Tw,Po)-PotTemp(To,Po))/3000. + mixingRatio(esatWater(Tw),Po)/1000.0)/(epsilon*eo)
    R=math.log(R)
    Tdew=243.5*R/(17.67-R)
    return Tdew,RelativeHumidity(Tdew,To)
    
# Get the dewpoint temperature directly from first law of thermo.
# INPUTS:
#     Wetbulb and air temperature in Celsius, Tw and To
#     air pressure in mb, Po.
# OUTPUTS:
#     Tdew in Celsius.
#     Relative humidity in %.
def getTdew(Tw,To,Po):
    R=esatWater(Tw)/eo - Po*cp*(To-Tw)/(epsilon*eo*LatentHeatWater(To))
    R=math.log(R)
    Tdew=243.5*R/(17.67-R)
    return Tdew,RelativeHumidity(Tdew,To)


# Get LCL P and T given Tdew, To, and Po at the level o.
# Theory: mixing ratio and potential temperature are constant (conserved) below the LCL.
# Solve each of these equations for mixing ratio and potential temperature for pressure.
# Set these equations for pressure equal to each other and solve for temperature.
# Then solve for pressure.  
# INPUTS:
#   Po in mb
#   To in Celsius
#   Tdew in Celsius.
# OUTPUTS:
#   P in mb, pressure at the LCL.
#   Tc in Celsius, temperature at the LCL.
#   zo in meters, the height of the starting level o as a pressure altitude.
#   z in meters, the height of the LCL level pressure altitude.
#   thetaK in Kelvin, the potential temperature below the LCL (constant value).
#   wo in kg/kg, the water vapor mixing ratio below the LCL (constant value)
#   Tw in Celsius, wet bulb temperature at Po.
#   ThetaW in Kelvin, wet bulb temperature at 1000 mb.
def LCL(Po,To,Tdew):
    epsilon=0.622
    e0=esatWater(0.0) # Saturation vapor pressure at T=0 C.
    zo=PressureAltitude(285.0,1013.25,Po) # Pressure altitude at the initial point.
    esatTdew=esatWater(Tdew) # Saturation vapor pressure at the dew point temp.
    wo=mixingRatio(esatTdew,Po)/1000.0 # kg/kg value.
    beta=epsilon*e0/(wo*1000.0) # Intermediate value.
    thetaK=PotTemp(To,Po) + 273.15 # Potential temperature in Kelvin.
    thetaE_K=thetaK+3000.0*wo
    Tguess=-10.0 + 273.15 # Initial guess in Kelvin.
#print " Itterating to get the LCL temperature"
    for i in range(200):
        T=273.15+(3.5*math.log(Tguess/thetaK)-math.log(beta))*(243.5+Tguess-273.15)/17.67
#    print "Tguess={0:1.3f}".format(Tguess),"K"," T={0:1.3f}".format(T),"K"
        Tguess=T
    Tc=T-273.15 # temperature at the LCL in Celsius.
    P=1000.0*((T/thetaK)**3.5) # Pressure at the LCL in mb.
    z=PressureAltitude(To+273.15,Po,P)+zo # Pressure altitude at the LCL in meters.
# Get the wetbulb temperature and the wetbulb potential temperature.
    PressureValues=[Po,1000.0] # Pressure values for wet bulb and wetbulb potential temp are needed.
    nudge = 0.3 # Used for itterating towards a solution, use the new value with the old value.
    for pref in PressureValues:  # Loop over pressure values
        Tguess=(Tdew+To)/2.0 + 273.15 # Initial estimate for wetbulb
        for i in range(100):
            Tw=nudge*(thetaE_K-3.0*mixingRatio(esatWater(Tguess-273.15),pref))*((pref/1000.0)**(2./7.)) + (1.0-nudge)*Tguess
            Tguess=Tw
        if pref == Po:
            Twvalue=Tw-273.15
        elif pref == 1000.0:
            thetaW=Tw
    return P,Tc,zo,z,thetaK,wo,thetaE_K,Twvalue,thetaW


# Demonstrate the LCL calculator.
# Input values
Po=870.0 # Pressure at the measurement level.
To=20.0  # Temperature at the measurement level.
Tdew=5.0 # Dewpoint temperature at the measurement level.
# Call the function
P,Tc,zo,z,thetaK,wo,thetaE_K,Tw,thetaW = LCL(Po,To,Tdew)

# Summary the LCL calculation output. 
print "LCL CALCULATOR RESULTS"
print "INPUT VALUES:"
print "     Po={0:1.1f}".format(Po),"mb"
print "     To={0:1.2f}".format(To),"C"
print "     Tdew={0:1.2f}".format(Tdew),"C"
print "     Calculated Pressure Altitude={0:1.1f}".format(zo),"meters"
print "CONSTANT VALUES BELOW THE LCL:"
print "     Potential Temperature={0:1.1f}".format(thetaK),"K ={0:1.1f}".format(thetaK-273.15),"C"
print "     Water vapor mixing ratio={0:1.2e}".format(wo),"kg/kg ={0:1.2f}".format(wo*1000.0),"g/kg"
print "CONSTANT VALUE EVERYWHERE ON THE PARCEL TRAJECTORY:"
print "     Equivalent Potential Temperature={0:1.1f}".format(thetaE_K),"K ={0:1.1f}".format(thetaE_K-273.15),"C"
print "WET BULB RELATED VALUES:"
print "     Wetbulb Temperature={0:1.1f}".format(Tw+273.15),"K ={0:1.1f}".format(Tw),"C"
print "     Wetbulb Potential Temperature={0:1.1f}".format(thetaW),"K ={0:1.1f}".format(thetaW-273.15),"C"
print "LCL VALUES:"
print "     P={0:1.1f}".format(P),"mb"
print "     T={0:1.1f}".format(Tc+273.15),"K ={0:1.2f}".format(Tc),"C"
print "     Calculated Pressure Altitude={0:1.1f}".format(z),"meters\n"


# Demonstrate saturation water vapor partial pressure calculator.
print "SATURATON VAPOR PRESSURE CALCULATOR RESULTS"
esatTo=esatWater(To)
print "esatWater(To)={0:1.2f}".format(esatTo),"mb\n"

# Demonstrate relative humidity calculator.
print "RELATIVE HUMIDITY CALCULATOR RESULTS"
RH = RelativeHumidity(Tdew,To)
print "RH(Tdew={0:1.1f}".format(Tdew),"C, T={0:1.1f}".format(To),"C) = {0:1.2f}".format(RH),"%\n"

# Demonstrate water vapor mixing ratio calculator.
# Calculate the volume mixing ratio.  
print "SATURATION WATER VAPOR MIXING RATIO CALCULATOR RESULTS"
# P=860.0          # total pressure in mb.
# T=275.0-273.15   # temperature in Celsius.
P=Po
T=To
eps=0.622
Wsat=mixingRatio(esatWater(T),P)
print "INPUT VALUES:"
print "     P={0:1.1f}".format(P),"mb"
print "     T={0:1.2f}".format(T),"C"
print "Calculated Wsat={0:1.2f}".format(Wsat),"g/kg\n"
#print "Calculated VMRsat={0:1.5e}".format(0.001*Wsat/eps),"partsH20perPartsDryAir\n"

# Demonstrate water vapor mixing ratio calculator.
# Calculate the volume mixing ratio.  
print "WATER VAPOR MIXING RATIO AT THE WET BULB TEMPERATURE "
# P=860.0          # total pressure in mb.
# T=275.0-273.15   # temperature in Celsius.
P=Po
T=Tw
eps=0.622
Wsat=mixingRatio(esatWater(T),P)
print "INPUT VALUES:"
print "     P={0:1.1f}".format(P),"mb"
print "     T={0:1.2f}".format(T),"C"
print "Calculated Wsat={0:1.2f}".format(Wsat),"g/kg"
#print "Calculated VMRsat={0:1.5e}".format(0.001*Wsat/eps),"partsH20perPartsDryAir\n"

# Demonstrate the pressure altitude calculator.
#To=290.0 # Kelvin Assumed temperature value at the surface.
#Po=1013.25 # mb Taken to be sea level pressure. 
#Pz=400.0 # mb Pressure at the unknown altitude.
#z=PressureAltitude(To,Po,Pz)
#print "PRESSURE ALTITUDE CALCULATOR RESULTS"
#print "INPUT VALUES:"
#print "     Po={0:1.1f}".format(Po),"mb"
#print "     To={0:1.2f}".format(To),"C"
#print "     Pz={0:1.1f}".format(Pz),"mb"
#print "Calculated Pressure Altitude={0:1.1f}".format(z),"meters\n"

# Demonstrate the Tdew calculator from Tw, To, and Po,
# using ThetaE.
#To=21.3  # Celsius air temperature value in the class room.
#Po=871.0 # mb Taken to be pressure in the class room. 
#Tw=12.0  # Celsius air temperature value in the class room..
#Tdew,RH=getTdewFromThetaE(Tw,To,Po)
#print "TDEW CALCULATOR RESULTS FROM Tw, To, and Po"
#print "This version uses ThetaE"
#print "INPUT VALUES:"
#print "     Po={0:1.2f}".format(Po),"mb"
#print "     To={0:1.2f}".format(To),"C"
#print "     Tw={0:1.2f}".format(Tw),"C"
#print "Tdew={0:1.2f}".format(Tdew),"C"
#print "RH={0:1.2f}".format(RH),"%\n"

# Demonstrate the Tdew calculator from Tw, To, and Po,
# using the first law of thermo directly.
#To=21.3  # Celsius air temperature value in the class room.
#Po=871.0 # mb Taken to be pressure in the class room. 
#Tw=12.0  # Celsius air temperature value in the class room..
#Tdew,RH=getTdew(Tw,To,Po)
#print "TDEW CALCULATOR RESULTS FROM Tw, To, and Po"
#print "This version uses the first law of thermo"
#print "INPUT VALUES:"
#print "     Po={0:1.2f}".format(Po),"mb"
#print "     To={0:1.2f}".format(To),"C"
#print "     Tw={0:1.2f}".format(Tw),"C"
#print "Tdew={0:1.2f}".format(Tdew),"C"
#print "RH={0:1.2f}".format(RH),"%\n"