Attribute VB_Name = "Thermo_Funcs"
Option Explicit
'Thermodynamic chart software package.
'Developed by G. S. Stipanuk, Atmospheric Sciences Laboratory, White Sands
'Missile Range, New Mexico, 88002. Presented in "Algorithms for Generating
'a SKEW-T, log P diagram and computing selected meteorological quantities",
'U.S. Government publication ECOM-5515, published October, 1973, and available
'on microfiche from NTIS.
'Typed in (in FORTRAN) by Harold Reynolds, March 15, 1991.
'Translated into Visual Basic by Harold Reynolds, June 28 - July 2, 2009.
'The following subroutines approximate a thermodynamic chart.
'T is the temperature in Kelvin. Scalar in all functions except Z and CCL.
'TD is the dew point temperature. Ditto.
'P is the pressure in millibars. Ditto.
'TDS, TS, and PS are TD, T, and P at the surface.
'WBAR is the mean mixing ratio.
'O is really a theta.
'Soundings must be ordered by decreasing pressure.
'Harold's note: function names are capitalized, variables are not.
'LEGAL STUFF: To the best of my knowledge, these routines produce the correct
'results. However, YOU MUST VERIFY THEIR CORRECTNESS BEFORE USING THESE
'FUNCTIONS FOR ANY PROJECT! IF SOMETHING IS WRONG, I WOULD VERY MUCH LIKE
'TO CORRECT IT! YOU HAVE BEEN WARNED!
'-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Function ALCL(tds As Double, ts As Double, ps As Double) As Double
'Computes the pressure at the lifting condensation level.
'tds, ts in K, ps, alcl in mb. ABS = absolute value.
Dim aw As Double, ao As Double, pi As Double, x As Double
Dim i As Integer
aw = W(tds, ps)
ao = O(ts, ps)
pi = ps
For i = 1 To 10
x = 0.02 * (TMR(aw, pi) - TDA(ao, pi))
If Abs(x) < 0.01 Then Exit For
pi = pi * 2 ^ x
Next i
ALCL = pi
End Function
'-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Function CCL(pm As Double, p() As Double, t() As Double, td() As Double, _
wbar As Double, N As Integer) As Double
'Computes pressure at convective condensation level.
'N is the number of levels in the sounding. K is the last level below pm.
'PM is pressure at top of mixing layer.
'CCL and p in mb, t in Kelvin, war in g vapour/kg dry air.
Dim t() As Double, td() As Double, p() As Double, tq As Double, x As Double
Dim del As Double, pc As Double, a As Double
Dim k As Integer, j As Integer, i As Integer, L As Integer
Dim FoundIt As Boolean
wbar = 0
k = 1
Do While p(k) < pm
k = k + 1
Loop
k = k - 1
j = k - 1
If j >= 1 Then
'Compute the average mixing ratio. Log is natural logarithm
For i = 1 To j
L = i + 1
wbar = (W(td(i), p(i)) + W(td(L), p(L))) * Log(p(i) / p(L)) + wbar
Next i
End If
L = k + 1
tq = td(k) + (td(L) - td(k)) * Log(pm / p(k)) / Log(p(L) / p(k))
wbar = wbar + (W(td(k), p(k)) + W(tq, pm)) * Log(p(k) / pm)
wbar = wbar / (2 * Log(p(1) / pm))
'Find the level at which tmr - ts changes sign. TS is sounding temp.
FoundIt = False
For i = 1 To N
If TMR(wbar, p(i)) + 273.16 >= 0 Then
FoundIt = True
Exit For
End If
Next i
'Not found, exit with CCL = 0
CCL = 0
Exit Function
'Set up bisection routine
L = i - 1
del = p(L) - p(i)
pc = p(i) + 0.5 * del
a = (t(i) - t(L)) / Log(p(L) / pc) + 273.16
For i = 1 To 10
del = del / 2
x = TMR(wbar, pc) - t(L) - a * Log(p(L) / pc) + 273.16
'The SIGN(x,y) function is a FORTRAN intrinsic that replaces the sign of x
'with that of y. I
'had to make a separate function for it here.
pc = pc + SIGN(del, x)
Next i
CCL = pc
End Function
'-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Function CT(wbar As Double, pc As Double, ps As Double) As Double
'Computes the convective temperature.
'Wbar in g/kg, pc, ps in mb.
Dim tc As Double, ao As Double
tc = TMR(wbar, pc) + 273.16
ao = O(tc, pc)
CT = TDA(ao, ps)
End Function
'-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Function ESAT(t As Double) As Double
'Computes the saturation vapour pressure over water at temperature t.
'ESAT in mb, t in K. Log to base 10 is needed for this function.
Dim a0 As Double, a1 As Double, a2 As Double
a0 = 23.832241 - 5.02808 * Log10(t)
a1 = 0.00000013816 * 10 ^ (11.344 - 0.0303998 * t)
a2 = 0.0081328 * 10 ^ (3.49149 - 1302.8844 / t)
ESAT = 10 ^ (a0 - a1 + a2 - 2949.076 / t)
End Function
'-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Function FR(t As Double, td As Double) As Double
'Computes relative humidity. FR in percent, t, td in Kelvin.
FR = ESAT(td) / ESAT(t) * 100
End Function
'-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Function O(t As Double, p As Double) As Double
'Computes the dry adiabat through (t,p)
'O and T in K, p in mb.
O = t * (1000 / p) ^ 0.288
End Function
'-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Function OE(tds As Double, ts As Double, ps As Double) As Double
'Computes the potential equivalent / pseudo-equivalent temperature.
'tds, ts, OE in K, ps in mb.
'NOTE: The commented formula is the one from the paper and gives wrong results!
'I have used instead the formula from Holton, p. 331, with q as the mixing
'ratio of the parcel and T as the saturation temperature (temperature at the LCL).
Dim alift As Double, olift As Double, tlift As Double
alift = ALCL(tds, ts, ps)
olift = O(ts, ps)
tlift = TDA(olift, alift) + 273.16
OE = O(ts, ps) * Exp(2.6518986 * W(tds, ps) / tlift) - 273.16
'This is the code which gives the SATURATED OE
' atw = TW(tds,ts,ps) + 273.16
' OE = OS(atw,1000) -273.16
End Function
'-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Function OS(t As Double, p As Double) As Double
'Computes saturation adiabat through (t,p)
'OS, t in K, p in millibars (mb)
OS = t * (1000 / p) ^ 0.288 / Exp(-2.6518986 * W(t, p) / t)
End Function
'-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Function OW(tds As Double, ts As Double, ps As Double) As Double
'Computes potential wet bulb temperature.
'tds, ts and OW in K, ps in mb.
Dim atw As Double, aos As Double
atw = TW(tds, ts, ps) + 273.16
aos = OS(atw, ps)
OW = TSA(aos, 1000)
End Function
'-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Function TDA(O As Double, p As Double) As Double
'Computes temperature on a dry adiabat o (theta) at pressure p
'o, TDA in K, p in mb.
TDA = O * (p / 1000) ^ 0.288 - 273.16
End Function
'-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Function TE(tds As Double, ts As Double, ps As Double) As Double
'Computes equivalent temperature
'tds, td, TE in K, ps in mb.
Dim aoe As Double
aoe = OE(tds, ts, ps) + 273.16
TE = TDA(aoe, ps)
End Function
'-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Function TMR(W As Double, p As Double) As Double
'Computes temperature on mixing ratio w at pressure p.
'TMR in C, w in g/kg dry air, p in millibars.
Dim x As Double
x = Log10(W * p / (622 + W))
TMR = 10 ^ (0.0498646455 * x + 2.4082965) - 280.23475 + 38.9114 * _
((10 ^ (0.0915 * x) - 1.2035) ^ 2)
End Function
'-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Function TSA(OS As Double, p As Double) As Double
'Computes temperature on saturated adiabat os at pressure p.
'SIGN(a,b) replaces the algebraic sign of a with that of b.
Dim a As Double, tq As Double, d As Double, x As Double
Dim i As Integer
a = OS
tq = 253.16
d = 120
For i = 1 To 12
d = d / 2
'If the temperature difference x is small, exit the loop
x = a * Exp(-2.6518986 * W(tq, p) / tq) - tq * (1000 / p) ^ 0.288
If Abs(x) <= 0.0000001 Then Exit For
tq = tq + SIGN(d, x)
Next i
TSA = tq - 273.16
End Function
'-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Function TW(tds As Double, ts As Double, ps As Double) As Double
'Computes wet bulb temperature.
'tds, ts and TW in K, ps in mb.
Dim i As Integer
Dim aw As Double, ao As Double, pi As Double, x As Double
Dim ti As Double, aos As Double
aw = W(tds, ps)
ao = O(ts, ps)
pi = ps
For i = 1 To 10
x = 0.02 * (TMR(aw, pi) - TDA(ao, pi))
If Abs(x) <= 0.01 Then Exit For
pi = pi * 2 ^ x
Next i
ti = TDA(ao, pi) + 273.16
'The intersection has been found, now find a saturated adiabat through it.
aos = OS(ti, pi)
TW = TSA(aos, ps)
End Function
'-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Function W(t As Double, p As Double) As Double
'Computes the mixing ratio line through (t,p).
'T is in K, p in mb, W in g water /kg dry air.
Dim x As Double
x = ESAT(t)
W = 622 * x / (p - x)
End Function
'-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Function Z(pt As Double, p() As Double, t() As Double, td() As Double, _
N As Integer) As Double
'Computes thickness in metres from p(1) to pt.
Dim i As Integer, j As Integer
Dim a1 As Double, a2 As Double, Z1 As Double
Z1 = 0
For i = 1 To N
j = i + 1
If pt >= p(j) Then Exit For
a1 = t(j) * (1 + 0.0006078 * W(td(j), p(j)))
a2 = t(i) * (1 + 0.0006078 * W(td(i), p(i)))
Z1 = Z1 + 14.64285 * (a1 + a2) * Log(p(i) / p(j))
Next i
a1 = t(j) * (1 + 0.0006078 * W(td(j), p(j)))
a2 = t(i) * (1 + 0.0006078 * W(td(i), p(i)))
Z = Z1 + 14.64285 * (a1 + a2) * Log(p(i) / pt)
End Function
'-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Function Log10(ByVal x As Double) As Double
'Computes logarithm to base 10
Log10 = Log(x) / Log(10)
End Function
'-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Function SIGN(x As Double, y As Double) As Double
'Replaces the sign of x with that of y. This is used to mimic the intrinsic
'function SIGN in Fortran.
If y < 0 Then
SIGN = -Abs(x)
Else
SIGN = Abs(x)
End If
End Function