added files

.gitignore

@@ -1,3 +1,28 @@
+/*.tar.gz
+/tmp
+/dist
+/dist-extras
+/julia
+/usr
+/usr-staging
+/Make.user
+/julia-*
+/source-dist.tmp
+/source-dist.tmp1
+
+*.exe
+*.dll
+*.dwo
+*.do
+*.o
+*.obj
+*.so
+*.dylib
+*.dSYM
 *.jl.cov
 *.jl.*.cov
 *.jl.mem
+*.ji
+deps/deps.jl
+
+.DS_Store

REQUIRE

@@ -1 +1,2 @@
 julia 0.5
+Unitful

src/Thermo.jl

@@ -1,5 +1,192 @@
 module Thermo
 
-# package code goes here
+using Unitful
+using Unitful: K, °C, Pa, hPa
+using Unitful: NoUnits, Temperature, Pressure, Amount
+#using Unitful: Length, Mass
+#using Unitful: Units
+#=
+# extend Unitful.unit function to be able to find the units of arrays
+Unitful.unit(a::AbstractArray)=Unitful.unit(a[1])
+Unitful.dimension(a::AbstractArray)=Unitful.dimension(a[1])
+=#
 
-end # module
+# Practical extensions to Unitful:
+# Define a nondimensional number to include nondimensional unitless builtin Number types but not dimensional Quantities.
+# Define a union of NondimensionalNumber and DimensionlessQuantity types.
+# Includes dimensionless ratios of Pressure, Mass, Temperature quantities.
+# It's easier to dispatch on Unions than typeintersections.
+typealias NondimensionalNumber Union{Real, Complex} # not including Quantities
+typealias NondimensionalQuantity Union{NondimensionalNumber, DimensionlessQuantity}
+NdN = NondimensionalNumber # for short
+NdQ = NondimensionalQuantity
+
+# Extend uconvert to facilitate consistent interfaces for functions of Unitful and/or NondimensionalNumber inputs
+# by converting to expected units, dimensional or nondimensional, for all inputs as they come in. E.g:
+# most specific method calculates:
+#   theta(T,::Temperature, p::Pressure, p0::Pressure, kappa::NondimensionalNumber) = K(T)/Exner(p, p0, kappa)
+# one catch-all method handles all unit conversions:
+#   theta(T, p, p0, kappa) = theta(uconvert(K,T), uconvert(hPa,p), uconvert(hPa,p0), kappa)
+"""
+ucon(args...)=Unitful.uconvert(args...)
+ucon(a::Units, x::NondimensionalNumber) = a*x
+extends and shadows Unitful.uconvert to add units to NondimensionalNumbers
+x::NondimensionalNumber assumed to be in a::Units.
+examples:
+ucon(K, x*°C)      # Unitful.uconvert standard method
+ucon(K, x+273.15)  # extended to NondimensionalNumbers
+"""
+ucon(args...)=Unitful.uconvert(args...)
+ucon(a::Unitful.Units, x::NondimensionalNumber) = a*x
+ucon{T<:NondimensionalNumber}(a::Unitful.Units, x::AbstractArray{T}) = ucon.(a,x)
+#ucon(a::Unitful.Units, x::NondimensionalNumber, sourceunits::Unitful.Units) = ucon(a, x*sourceunits) # superfluous
+
+"""
+Nondimensional unit convert
+undconvert(a::Unitful.Units, x::Unitful.Quantity) = NoUnits(ucon(a,x)/a)
+undconvert(a::Unitful.Units, x::NondimensionalNumber) = x
+outputs nondimensional number assumed in the requested units. Assumes nondimensional x already is in units of a.
+"""
+undconvert(a::Unitful.Units, x::Unitful.Quantity) = NoUnits(ucon(a,x)/a)
+undconvert(a::Unitful.Units, x::NondimensionalNumber) = x # == NoUnits(ucon(a,a*x)/a)
+
+export theta, ev,qv,es,qs, Tlcl, Lv, dqsdT, moistad
+
+include("ThermoConstants.jl")
+using .constants_unitful
+using .constants_unitful: pref, Rd, Rv, Cp, Cpv, Cw
+using .constants_unitful: EarthGravity
+
+"T/Θ = Exner(p, p0=1000hPa, kappa=Rd/Cp)  Exner function"
+Exner(p::NdN  , p0::NdN, kappa=NoUnits(Rd/Cp)) = (p/p0)^kappa # both unitless
+Exner(p::Pressure, p0::Pressure, kappa=NoUnits(Rd/Cp)) = NoUnits(p/p0)^kappa # both pressure, but may be of different units
+#Exner(p::Pressure, p0::Pressure=hPa(pref), kappa=NoUnits(Rd/Cp)) = NoUnits(p/p0)^kappa # define intersecting specific case to avoid method ambiguity
+# weird cases where one of p,p0 has units but the other doesn't probably should be errors, but here we convert them to have same units as the unitful one:
+Exner(p::NdN, p0::Pressure=hPa(pref), kappa=NoUnits(Rd/Cp)) = (p/NoUnits(p0/unit(p0)))^kappa
+Exner(p::Pressure, p0::NdN=NoUnits(pref/hPa), kappa=NoUnits(Rd/Cp)) = (NoUnits(p/unit(p))/p0)^kappa
+# extend to broadcasting arrays, iterables with .
+Exner(p::AbstractArray, p0, kappa)=Exner.(p, p0, kappa)
+
+"Exner exponent kappa(qv[kg/kg])=(Rd/Cp)*(1-0.28*qv) varies with the specific humidity."
+kappa(qv::NdQ)=NoUnits(Rd/Cp*(1-0.28*qv))
+
+"""
+theta(T,p, p0=1000hPa, kappa=Rd/Cp)
+calculates potential temperature from
+temperature T (default K), pressure p (default hPa).
+Quantities with other compatible units may be specified
+by the Unitful module.
+"""
+# unitless input assumes T in K and returns unitless T values corresponding to K
+theta(T::NdN,      p, p0=hPa(pref), kappa=NoUnits(Rd/Cp)) =   T /Exner(p, p0,kappa)
+theta(T::Temperature, p, p0=hPa(pref), kappa=NoUnits(Rd/Cp)) = K(T)/Exner(p, p0,kappa)
+# output unitful quantities with units K, automatically converts Celsius to K
+theta(T::AbstractArray,p,p0=hPa(pref),kappa=NoUnits(Rd/Cp))=theta.(T,p,p0,kappa) # extend to arrays
+# theta(T,p::AbstractArray,p0=hPa(pref),kappa=NoUnits(Rd/Cp))=theta.(T,p,p0,kappa)
+# multiple dispatch of p is handed by Exner
+
+"""
+qv(p/ev)
+qv(p::Pressure,ev::Pressure)
+qv(p::Amount,ev::Amount)
+specific humidity [kg/kg]
+"""
+#qv(p::Pressure,ev::Pressure) = NoUnits(Rd/(Rv*(p/ev + (Cp/Rv-1.))))
+#qv(p::Amount,ev::Amount) = NoUnits(Rd/(Rv*(p/ev + (Cp/Rv-1.))))
+qv(poev::Number) = ucon(u"kg/kg",Rd/(Rv*(poev + (Cp/Rv-1.))))
+qv(p::Number,ev::Number) = qv(p/ev)
+#qv(p::Union{AbstractArray,Number},ev::Union{AbstractArray,Number})=qv.(p,ev) # handle abstractarrays
+
+"""
+es(T,p) [hPa] computes saturation vapor pressure based on Wexler's formula,
+with enhancement factor for moist air rather than water vapor.
+The enhancement factor requires a pressure.
+T [degrees C], p [hPa] (note the reversed input order), es [hPa]
+From A. L. Buck 1981: JAM, 20, 1527-1532.
+SPdeS 7 July 2004
+"""
+# procedure: convert to expected units, nondimensionalize, calculate
+# methods "diagonally dispatch" nondimensionalizing one argument at a time
+es(T::NdN,p::NdN=NoUnits(pref/hPa)) = 6.1121hPa*(1.0007 + 3.46e-8*p).*exp((17.502*T)./(240.97 + T))
+es(T::Temperature,p::Pressure=pref) = es(NoUnits(°C(T)/°C),NoUnits(p/hPa)) # resolves ambiguity below
+# weird corner cases probably should be errors, but are allowed for now:
+es(T::Temperature,p::NdN) = es(NoUnits(°C(T)/°C),p)
+es(T::NdN,   p::Pressure) = es(T,NoUnits(p/hPa))
+# NoUnits form expects T in °C, p in hPa
+# uses Buck 1981 numerical constants
+#es(T::AbstractArray,p               )=es.(T,p) # broadcasts
+#es(T::AbstractArray,p::AbstractArray)=es.(T,p)
+#es(T               ,p::AbstractArray)=es.(T,p)
+
+"""
+qs(T[°C],p[hPa]) [kg/kg]
+Saturation specific humidity based on Wexler's formula for es
+with enhancement factor. See es(T,p).
+From A. L. Buck 1981: JAM, 20, 1527-1532.
+SPdeS 7 July 2004
+"""
+qs(T::Temperature,p::Pressure=pref)=qv(p,es(T,p))
+qs(T::NdN,p::NdN=NoUnits(pref/hPa))=qv(p,NoUnits(es(T,p)/hPa))
+#qs(T::AbstractArray,p::AbstractArray)=qs.(T,p)  # broadcast
+#qs(T::AbstractArray,p)=qs.(T,p)
+#qs(T,p::AbstractArray)=qs.(T,p)
+
+"ev(p,qv) vapor pressure in units of p"
+ev(p::Pressure,qv::NdQ) = p/NoUnits(Rd/(Rv*qv) +1)
+ev(p::NdN,qv::NdQ) = p/NoUnits(Rd/(Rv*qv) +1)
+#ev(p,qv)=ev.(p,qv)
+#  ev = p.*qv./(Rd/Rv + qv)
+#  Tl = 2840./(3.5*log(T) - log(.01*ev) - 4.805) + 55
+
+"""
+moistad(T[°C], p[hPa])
+Moist adaibatic lapse rate of _potential_ temperature [K/m]
+from Wood and Bretherton 2006
+(also Rogers and Yau, but with mixing ratio instead of specific humidity).
+(c) Simon de Szoeke, 2011
+"""
+moistad(T::Temperature,p=pref)=moistad(NoUnits(°C(T)/°C), p)
+moistad(T::NdN,p::Pressure=pref)=moistad(T, NoUnits(p/hPa))
+function moistad(T::NdN,p::NdN=NoUnits(pref/hPa))
+    L=Lv(T)
+    qsat=qs(T,p)
+    (EarthGravity/Cp) * (1 - (1+L*qsat/(Rd*(T+273.15))) / (1+((L/(T+273.15))^2 *qsat/(Cp*Rv))) )
+end
+
+"Lv(T) [J/kg] Latent heat of vaporization of water as a function of temperature [°C]."
+Lv(T::NdN=0.0) = (2.501e6 + (Cpv-Cw)*T/unit(Cpv))*u"J/kg"
+Lv(T::Temperature=0.0°C) = (2.501e6u"J/kg" + (Cpv-Cw)*°C(T)) |> u"J/kg" # more specific method called first
+Lv{Y<:Union{Range,AbstractArray}}(T::Y)=Lv.(T)
+# Lv(T::AbstractArray)=Lv.(T)
+
+"dqsdT(T[K],p[hPa]) = dqs/dT = qs(p,T)*Lv(T)/(Rv*T^2) [1/K] using Wexler and Clausius-Clapeyron formulas"
+dqsdT(T::Temperature,p::Pressure=pref) = qs(T,p)*Lv(T)/(Rv*(K(T)^2))
+dqsdT(T::Temperature,p::NdN) = dqsdT(T,p*hPa)
+dqsdT(T::NdN,p::Pressure=pref) = dqsdT(T*K,p)
+dqsdT(T::NdN,p::NdN) = dqsdT(T*K,p*hPa)
+dqsdT(T::AbstractArray,p=pref) = dqsdT.(T,p) # broadcasts
+dqsdT(T,p::AbstractArray) = dqsdT.(T,p)
+
+
+"""
+Tlcl(T[K], p[hPa], qv[kg/kg])
+Tlcl(T[K], p[hPa], ev[Pa])
+Temperature at the LCL [K]. From Bolton, 1980, MWR, 108, 1046-1053.
+"""
+Tlcl(T::Temperature,p,qv) = Tlcl(NoUnits(K(T)/K),p,vapor)
+Tlcl(T::NdN,p::Pressure,qv) = Tlcl(T,NoUnits(p/Pa),vapor)
+Tlcl(T::NdN,p::NdN,ev::Pressure) = (2840./(3.5*log(T) - log(NoUnits(ev/hPa)           ) - 4.805) + 55.)*K
+Tlcl(T::NdN,p::NdN,qv::NdQ) =   (2840./(3.5*log(T) - log(p/(Rd/(Rv*NoUnits(qv)) +1)) - 4.805) + 55.)*K
+Tlcl(T,p,q)=Tlcl.(T,p,q)
+
+"theta_e(T[K],p[hPa],qv[kg/kg]) from eqn. 43 of Bolton, 1980, MWR, 108, 1046-1053."
+theta_e(T::Temperature,p::Pressure,qv::NdQ) = T/Exner(p,p0,kappa(qv)) * exp((3376./Tlcl(T,p,qv) - 2.54)*qv*(1 + 0.81*qv))
+theta_e(T,p,qv) = theta_e(ucon(K,T),ucon(hPa,p),NoUnits(ucon(u"kg/kg",qv)))
+
+include("theta_w.jl")
+
+"Blackbody radiation bb(T)"
+bb(t::Temperature)=ucon(u"W/m^2", Unitful.σ*K(t)^4) # specific method calculates
+bb(t)=bb(ucon(K, t)) # wrapper converts arguments' units
+
+end

src/ThermoConstants.jl

@@ -0,0 +1,104 @@
+"Useful constants for atmospheric science on Earth. Uses Unitful constants."
+module constants_unitful
+
+using Unitful
+
+#export UniversalGasConst, AvogadroNumber, BoltzmannConstant, SpeedOfLight
+#export PermittivityOfVacuum, PlanckConstant
+#export AirDryMolecularWeight, AirDrySpecificHeatConstPressure, AirDrySpecificHeatConstVolume
+#export AirDryGasConst, AirDryThermalConductivity
+#export Cp, Cpv, Cw
+#export H2OMolecularWeight, H2OGasConst, LatentHeatVapor, LatentHeatFusion
+#export EarthGravity, EarthAngularVelocity, EarthSunSurfaceDistance
+#export EarthRadius, SolarFluxTOA
+#export p0, KelvinCelsius
+
+# physical and chemical constants
+UniversalGasConst=Unitful.R # 8314.3 # J/K/kmol
+AvogadroNumber=Unitful.Na # 6.022e26 # Avogadro's number
+BoltzmannConstant=Unitful.k # 1.381e-23 # J/K/molecule Boltzmann's constant
+SpeedOfLight=Unitful.c0 #2.998e8 # m/s
+PermittivityOfVacuum=Unitful.ϵ0 # 8.85e-12  # C^-2 N^-1 m^-2 permittivity of vacuum
+PlanckConstant=Unitful.h # 6.6262e-34 # J*s
+
+# properties of dry air
+AirDryMolecularWeight=28.97u"kg/kmol"
+AirDrySpecificHeatConstPressure=1004.0u"J/K/kg" # Wallace and Hobbs
+AirDrySpecificHeatConstVolume=717.0u"J/K/kg"
+#AirDryGasConst=287.0u"J/K/kg"
+AirDryThermalConductivity=2.40e-2u"J/m/s/K" # at 0C
+Rd=287.04u"J/K/kg" # J/K/kg Bolton
+Cp=1005.7u"J/K/kg" # J/K/kg Bolton
+
+# properties of water
+H2OMolecularWeight=18.016u"g/mol"
+H2OGasConst=461u"J/K/kg"
+LatentHeatVapor=2.5e6u"J/kg" # at 0 C
+LatentHeatFusion=3.34e5u"K/kg"
+Rv=461.5u"J/kg/K" # Bolton
+# specific heats of vapor and liquid
+Cpv=1870u"J/kg/K" # Bolton
+Cw=4190u"J/kg/K"  # Bolton
+
+# properties of Earth
+EarthGravity=9.8u"m/s/s"            # m/s/s on gravity on Earth surface
+EarthAngularVelocity=7.292e-5u"s^-1"
+EarthSunSurfaceDistance=1.50e11u"m"
+EarthRadius=6.37e6u"m"
+SolarFluxTOA=1.38e3u"W/m^2"
+
+# reference values
+#hPa=100u"Pa"
+pref=1.e5u"Pa"
+#KelvinCelsius=273.15 # K
+end # module constants_unitful
+
+module constants_unitless
+
+#export UniversalGasConst, AvogadroNumber, BoltzmannConstant
+#export PermittivityOfVacuum, PlanckConstant, SpeedOfLight
+#export AirDryMolecularWeight, AirDrySpecificHeatConstPressure, AirDrySpecificHeatConstVolume
+#export AirDryThermalConductivity
+#export Rd, Cp, Cpv, Cw
+#export H2OMolecularWeight, H2OGasConst, LatentHeatVapor, LatentHeatFusion
+#export EarthGravity, EarthAngularVelocity, EarthRadius
+#export EarthSunSurfaceDistance, SolarFluxTOA
+
+# physical and chemical constants
+UniversalGasConst=8314.31 #J/K/kmol
+AvogadroNumber=6.022e26 # Avogadro's number
+BoltzmannConstant=1.381e-23 # J/K/molecule Boltzmann's constant
+SpeedOfLight=2.998e8 # m/s
+PermittivityOfVacuum=8.85e-12  # C^-2 N^-1 m^-2 permittivity of vacuum
+PlanckConstant=6.6262e-34 # J/s
+
+# properties of dry air
+AirDryMolecularWeight=28.97 # kg/kmol
+AirDrySpecificHeatConstPressure=1004.  # J/K/kg, Wallace and Hobbs
+AirDrySpecificHeatConstVolume=717.   # J/K/kg
+AirDryGasConst=287.   # J/K/kg
+AirDryThermalConductivity=2.40e-2 # J/m/s/K at 0C
+Rd=287.04 # J/K/kg Bolton
+Cp=1005.7 # J/K/kg Bolton
+
+# properties of water
+H2OMolecularWeight=18.016 # kg/kmol
+H2OGasConst=461        # J/K/kg
+LatentHeatVapor=2.5e6  # J kg-1 at 0 C
+LatentHeatFusion=3.34e5 # K kg-1
+Rv=461.5 # Bolton
+# specific heats of vapor and liquid
+Cpv=1870 # Bolton
+Cw=4190  # Bolton
+
+# properties of Earth
+EarthGravity=9.8                    # m/s/s on gravity on Earth surface
+EarthAngularVelocity=7.292e-5       # s-1
+EarthSunSurfaceDistance=1.50e11     # m
+EarthRadius=6.37e6                  # m
+SolarFluxTOA=1.38e3                 # W/m2
+
+# reference values
+pref=1.e5   # hPa
+KelvinCelsius=273.15 # K
+end # constants unitless

src/theta_w.jl

@@ -0,0 +1,66 @@
+"""
+theta_w(p[Pa],Temp[K],qv[kg/kg]) from eqn 3.8 of Davies-Jones 2008 and
+eqn. 43 of Bolton, 1980, MWR, 108, 1046-1053.
+"""
+theta_w(p::Number,T::Number,qv::Number) = theta_w(uconvert(hPa,p),uconvert(K,T),NoUnits(uconvert(u"kg/kg",qv))) # wrapper handles diverse units
+theta_w(p,T,qv)=theta_w.(p,T,qv) # array broadcasting
+function theta_w(p::Pressure,T::Temperature,qv::NondimensionalQuantity) # specific function proper
+  # Davies-Jones rational function coefficients
+  a0=  7.101574
+  a1=-20.68208
+  a2= 16.11182
+  a3=  2.574631
+  a4= -5.205688
+  b1= -3.552497
+  b2=  3.781782
+  b3= -0.6899655
+  b4= -0.5929340
+
+  # compute Theta_e from Bolton as first guess approx of thw
+  thw = theta_e(p,T,qv);
+
+  # evaluates Davies-Jones 2008 rational function fit for thw>=173.15 [here in K]
+  if thw>=173.15
+    X=thw/273.15
+    X2=X*X
+    X3=X2*X
+    X4=X2*X2
+    thw=thw-exp((a0+a1*X+a2*X2+a3*X3+a4*X4)/(1.0+b1*X+b2*X2+b3*X3+b4*X4));
+  end
+
+end
+
+#=
+# constants
+global Rd Cp
+A = 2675;  %K
+C = 273.15; %K
+lam = Cp/Rd;
+ep=0.622
+es0=6.112 % hPa
+a=17.67;
+b=243.5;
+
+% Davies-Jones 2008 method resembles Bolton
+if the<=257
+   dlnesdlnthe=a*b/((Temp-C+b).*(Temp-C+b));
+   thw=the-C-A*rs/(1+A*rs.*dlnesdlnthe));
+elseif the
+elseif 377<=the & the<674
+end
+
+% doesn't work at all
+% initialize thw
+thw=max(min(-35,Twet(th-273.15,qv,1e5,5)),85)+273.15;
+thw=th;
+
+for iter=1:niter
+   % es(T) in hPa from es.m Buck
+   esat=6.1121*(1.0007 + 3.46e-8*1000).*exp((17.502*(thw-273.15))./(240.97 + thw-273.15));
+   %hPa                          %hPa
+   rs=622.*esat./(1000-esat); % saturation mixing ratio in hPa, Bolton eqn. 41
+   the=thw.*exp( (3.376./thw-0.00254).*rs.*(1+0.81e-3*rs) );
+   rath=thetae./the; % ratio of Bolton's thetae to guess thetae
+   thw=thw.*rath.^0.287;
+end
+=#