#---------------------Description--------------------------------------
#
#Computes baseline evaporation, which is the evaporation when Tg = Tsa,
#and also the latent heat, sensible heat, and radiative coupling coefficients.
#Used in making various tables in the Surface Energy Budget chapter
#
#**ToDo: Clean up setting of parameters, esp. es. It's very
#confusing.
#-----------------------------------------------------------------------

#Data on section of text which this script is associated with
Chapter = '6.**'
Figure = ''
Table = 'table:TypicalLatentFlux,table:CouplingCoeffs'
#

import phys

#Function definitions

#Computes the boundary layer density,
#taking into account mixture of condensible
#and noncondensible
def rho(Tbar,noncondensible, condensible):
    rhoa = pa/(noncondensible.R * Tbar)
    rhow = hsa*es(Tbar)/(condensible.R *Tbar)
    return rhoa+rhow
    

#Computes specific heat, taking into account
#mix of condensible and noncondensible
def cp(Tbar,noncondensible,condensible):   
    rhoa = pa/(noncondensible.R * Tbar)
    rhow = hsa*es(Tbar)/(condensible.R *Tbar)
    q = rhow/(rhoa+rhow)
    cpw = condensible.cp
    cpa = noncondensible.cp
    return q*cpw + (1.-q)*cpa

#Computes Richardson number
def Ri(Tg,Ta,rh,dry =False,z=10.):
    eta = rh*es(Ta)/(pa + es(Ta))
    etag = es(Tg)/(pa + es(Tg))
    if dry:
        eta = etag = 0.
    eps = condensible.MolecularWeight/noncon.MolecularWeight
    beta = g*(1.-(Ta/Tg)*(1.+(eps-1.)*etag)/(1.+(eps-1.)*eta))
    return -z*beta/(U*U)


#Computes CD as function of Ri in stable case
def CD(Ri,z=10.,z0=.001):
    Ric = .2
    if Ri > Ric:
        return 0.
    if Ri < 0.:
        Ri = 0. #Use neutral case if unstable
    return ((.4)**2/(math.log(z/z0)**2)) * (1.- Ri/Ric)**2

def evap(gas,T):
       es = phys.satvps_function(gas)
       if T < gas.TriplePointT:
         L = gas.L_sublimation
       else:
         L = gas.L_vaporization
       A = (L/(gas.R*T))
       Fstar = CD0*U*es(T)
       Eo = A*(1.-hsa)*Fstar
       b = A*A*Fstar/T
       return Eo,b

CD0 = .0015 #Drag coefficient (Doesn't take buoyancy effect into account)
U = 10    #Wind speed (m/s)
hsa = .7  #Boundary layer relative humidity
g = 9.8 #Used for Richardson number only

#-------Print out data for the latent heat flux table-------

print 'water',230.,evap(phys.water,230)
print 'water',273.,evap(phys.water,273.)
print 'water',280.,evap(phys.water,280.)
print 'water',300.,evap(phys.water,300.)
print 'water',320.,evap(phys.water,320.)

print 'CO2',150.,evap(phys.CO2,150.)
print 'CO2',160.,evap(phys.CO2,160.)
print 'CH4',80.,evap(phys.CH4,80.)
print 'CH4',95.,evap(phys.CH4, 95.)


#------Now print out data for the coupling constant table---------
def CoeffTable(case,Tlist):
    print 'case  T     bir     bsens     blat'
    for T in Tlist:
        bir = 4.*phys.sigma*T**3
        #Need to compute rho and Cp here, for mix of gases
        bsens = rho(T,noncon, condensible)*cp(T,noncon,condensible)*CD0*U
        blat = evap(condensible,T)[1]
        print case,T,bir,bsens,blat

#---Case 1: Water in air 
noncon = phys.air
condensible = phys.water
pa = 1.e5
es = phys.satvps_function(condensible)
Tlist = [250.,280.,300.,320.]
#
CoeffTable('water+air',Tlist)

#---Case 2: CH4 in N2 (Titan) 
noncon = phys.N2
condensible = phys.CH4
pa = 1.5e5
es = phys.satvps_function(condensible)
Tlist = [85.,90.,95.]
#
CoeffTable('CH4+N2',Tlist)