| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104 |
- # Misc :
- import sys
- import numpy as nmp
- import math
- rt0 = 273.16
- grav = 9.8 # gravity
- Rgas = 287.04
- Patm = 101000.
- ctv = 0.608 # for virtual temperature
- R_dry = 287.05 # Specific gas constant for dry air [J/K/kg]
- R_vap = 461.495 # Specific gas constant for water vapor [J/K/kg]
- reps0 = R_dry/R_vap # ratio of gas constant for dry air and water vapor => ~ 0.622
- cte = 0.622
- kappa = 0.4 # Von Karman's constant
- Cp = 1000.5
- Pi = 3.141592654
- eps_w = 0.987 # emissivity of water
- sigma = 5.67E-8 # Stefan Boltzman constamt
- alfa = 0.066 # Surface albedo over ocean
- rtt0 = 273.16 # triple point of temperature [K]
- sensit = 0.1
- def Lvap(zsst):
- #
- # INPUT : zsst => water temperature in [K]
- # OUTPUT : Lvap => Latent Heat of Vaporization [J/Kg]
- return ( 2.501 - 0.00237*(zsst - rt0) )*1.E6
- def e_air(q_air, zslp):
- #
- #--------------------------------------------------------------------
- # **** Function e_air ****
- #
- # Gives vapour pressure of air from pressure and specific humidity
- #
- #--------------------------------------------------------------------
- #
- diff = 1.E8
- e_old = q_air*zslp/reps0
- while diff > 1.:
- #print "Again... diff = ", diff
- ee = q_air/reps0*(zslp - (1. - reps0)*e_old)
- diff = nmp.sum(abs( ee - e_old ))
- e_old = ee
- return ee
- ### Update: June 2019, LB:
- def e_sat(rT):
- #
- #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- # rT: air temperature [K]
- # e_sat: water vapor at saturation [Pa]
- #
- # Recommended by WMO
- #
- # Goff, J. A., 1957: Saturation pressure of water on the new kelvin
- # temperature scale. Transactions of the American society of heating
- # and ventilating engineers, 347.
- #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- #
- ztmp1 = nmp.zeros(nmp.shape(rT))
- ztmp2 = nmp.zeros(nmp.shape(rT))
- ztmp2 = rT/rtt0
- ztmp1 = 1./ztmp2
- #
- #e_sat = 100.*( 10.^(10.79574*(1. - ztmp) - 5.028*LOG10(rT/rtt0) \
- # + 1.50475*10.^(-4)*(1. - 10.^(-8.2969*(rT/rtt0 - 1.)) ) \
- # + 0.42873*10.^(-3)*(10.^(4.76955*(1. - ztmp)) - 1.) + 0.78614) )
- #
- e_sat = 100.*( nmp.power(10.,(10.79574*(1. - ztmp1) - 5.028*nmp.log10(ztmp2) \
- + 1.50475*0.0001*(1. - nmp.power(10.,(-8.2969*(ztmp2 - 1.))) ) \
- + 0.42873*0.001 *(nmp.power(10.,(4.76955*(1. - ztmp1))) - 1.) + 0.78614 ) ) )
- #
- del ztmp1, ztmp2
- #
- return e_sat
- def q_air_dp(da, slp):
- #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- # Air specific humidity from dew point temperature
- # da !: dew-point temperature [K]
- # slp !: atmospheric pressure [Pa]
- #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- #
- es = e_sat(da)
- q_air_dp = es*reps0/(slp - (1. - reps0)*es)
- return q_air_dp
|
