Algorithms, Comparisons and Source References by Schlatter and Baker
From schlatter@profsc.fsl.noaa.gov Wed Jun 12 16:25:46 1991
General information:
--------------------
This is an index of the thermodynamic subprograms located in THERMO.OLB. They
are listed according to the general type of calculation that is desired,
divided into the following categories:
(a) moisture parameters
(b) latent heat
(c) pressure
(d) temperature
(e) thickness
These algorithms were collected, edited, commented, and tested by Thomas W.
Schlatter and Donald V. Baker from August to October 1981 in the PROFS Program
Office, NOAA Environmental Research Laboratories, Boulder, Colorado. Where
possible, credit has been given to the original author of the algorithm and a
reference provided.
The input/output units are as follows:
Temperature Celsius
Pressure millibars
Relative humidity percent
Saturation specific humidity g vapor/kg moist air
Mixing ratio g vapor/kg dry air
Thickness meters
Precipitable water centimeters
Latent heat joules/kg
The following symbols are used in subprogram calls:
EW water vapor pressure
EQPT equivalent potential temperature
P pressure
PC pressure at the convective condensation level
PS surface pressure
RH relative humidity
T temperature
TC temperature at the lifting condensation level
TD dew point
TH potential temperature (a dry adiabat)
THW wet bulb potential temperature (a moist adiabat)
W mixing ratio
WBAR mean mixing ratio from surface to pressure pm
Note: all routines are function subprograms except for "PTLCL", which is a
subroutine.
Index: routines available for various calculations
---------------------------------------------------
(A) Moisture Parameters
To Calculate: Available Routine(s):
1. Saturation mixing ratio 1. a) WMR(P,T)
b) W(T,P)
2. Precipitable water 2. PRECPW(TD,P,N)
(Td and P are dimensioned by N)
3. Relative humidity 3. HUM(T,TD)
4. Saturation specific humidity 4. SSH(P,T)
(B) Latent Heat
To Calculate Latent Heat Of: Available Routines:
1. Evaporization/Condensation 1. HEATL(1,T)
2. Melting/Freezing 2. HEATL(2,T)
3. Sublimation/Deposition 3. HEATL(3,T)
(C) Pressure
To Calculate: Available Routine(s):
1. Pressure at Lifting 1. a) ALCL(T,TD,P)
Condensation Level b) PCON(P,T,TC)
c) (subroutine) PTLCL(P,T,TD,PC,TC)
2. Saturation vapor pressure 2. a) ESW(T)
over liquid water b) ES(T)
c) ESAT(T)
d) ESGG(T)
e) ESRW(T)
f) ESLO(T)
3. Saturation vapor pressure 3. a) ESICE(T)
over ice b) ESILO(T)
4. Pressure at Convective 4. PCCL(PM,P,T,TD,WBAR,N)
Condensation Level (P, T, and
Td are dimensioned by N)
(D) Temperature
To Calculate: Available Routine(s):
1. Difference in Wet Bulb Poten- 1. WOBF(T)
tial Temperatures for two
parcels of air at the same
temperature - one saturated
and the other completely dry
2. Equivalent temperature 2. TE(T,TD,P)
3. Equivalent potential 3. a) OE(T,TD,P)
temperature b) EPT(T,TD,P)
c) OS(T,P)
4. Dew point 4. a) DPT(EW)
b) DEWPT(EW)
c) DWPT(T,RH)
5. Wet bulb potential 5. a) POWT(T,P,TD)
temperature b) OW(T,TD,P)
c) THM(W,P)
6. Potential temperature 6. O(T,P)
7. Convective temperature 7. CT(WBAR,PC,PS)
8. Virtual temperature 8. TV(T,TD,P)
9. Wet bulb temperature 9. TW(T,TD,P)
10. Temperature at Lifting 10. a) TCON(T,TD)
Condensation Level b) TLCL(T,TD)
c) (subroutine) PTLCL(P,T,TD,PC,TC)
d) TLCL1(T,TD)
11. Temperature along moist 11. SATLFT(THW,P)
adiabat at given pressure
12. Temperature on a dry adiabat 12. TDA(TH,P)
at given pressure
13. Temperature on a moist adiabat 13. a) TSA(EQPT,P)
corresponding to an Equivalent b) TMLAPS(EQPT,P)
Potential Temperature "EQPT",
at given pressure
14. Temperature on a mixing ratio 14. TMR(W,P)
line at given pressure
(E) Thickness
To Calculate: Available Routine:
1. Thickness between the surface 1. Z(PT,P,T,TD,N)
and any other pressure level
PT in the sounding. (P,T, and
Td are dimensioned by N.)
Notes prepared by D. Baker, 24 Dec 86.
###############################################################################
General information:
--------------------
The following tables give results of the tests performed on
routines in [BAKER.THERMO]THERM.FOR. These routines calculate thermo-
dynamic parameters that are of interest when analyzing sound-
ings. the "Smithsonian Meterological Tables" were utilized to
test the accuracy of each routine where possible...otherwise,
a Skew T-Log p chart was used. In cases where more than one
routine calculates a given parameter, an efficiency test is also
included. For the efficiency tests, a specified number of calcula-
tions was done by each routine using identical input. After 10
trials, the average time needed by each routine was determined
and is listed, thereby indicating which routine is most efficient.
for further information involving any of the routines to follow,
refer to documentation in [BAKER.THERMO]THERM.FOR.
Computer used:
VAX 11/780 @ NOAA/ERL/PROFS...Boulder, Colorado.
Input/output units:
Temperature............celsius
Pressure...............millibars
Mixing ratio...........gm vapor/kg dry air
Relative humidity......percent
Sat specific humidity..gm vapor/kg moist air
Note: Smith = Smithsonian table value. Where these values do not
appear, the routines have been verified using a Skew T-
Log p chart due to lack of an applicable table. The number
in parentheses next to "Smith" indicates the page number
(6th revised edition) of the table used for the test.
###############################################################################
TEST OF FUNCTIONS ESW,ES,ESAT,ESGG,ESRW,ESLO
PURPOSE: CALCULATE SATURATION VAPOR PRESSURE OVER LIQUID WATER
GIVEN TEMPERATURE.
TEMP ESW ES ESAT ESGG ESRW ESLO SMITH.(340)
---- ----- ----- ------ ------ ------ ------ ------
-50 .06356 .06357 .06356 .06356 .06362 .06337 .06356
-30 .50878 .51036 .50878 .50880 .50843 .50877 .5088
-10 2.8627 2.8677 2.8626 2.8627 2.8625 2.8635 2.8627
0 6.1080 6.1121 6.1076 6.1078 6.1084 6.1078 6.1078
10 12.272 12.272 12.272 12.272 12.274 12.271 12.272
20 23.372 23.370 23.372 23.373 23.375 23.371 23.373
30 42.430 42.457 42.429 42.430 42.431 42.429 42.430
EFFICIENCY TEST: 5000 CALCULATIONS AT T=20C.
FUNC. T(SEC)
---- ------
ESW 1.12
ES 0.88
ESAT 4.50
ESGG 4.63
ESRW 1.08
ESLO 0.64
###############################################################################
TEST OF FUNCTIONS WMR,W
PURPOSE: CALCULATE MIXING RATIO GIVEN PRESSURE AND TEMPERATURE.
PRES TEMP WMR W SMITH.(302)
---- ---- ------ ------ ------
1000 30 27.699 27.560 27.69
1000 0 3.840 3.822 3.839
850 10 9.147 9.112 9.146
700 10 11.135 11.099 11.13
500 -20 1.568 1.564 1.568
400 -30 0.794 0.792 0.794
EFFICIENCY TEST: 2000 CALCULATIONS AT T=-10C, P=850MB.
FUNC. T(SEC)
----- ------
WMR 0.82
W 1.86
###############################################################################
TEST OF FUNCTIONS TCON,TLCL,TLCL1
PURPOSE: CALCULATE TEMPERATURE AT THE LCL GIVEN TEMPERATURE AND
DEW POINT.
TEMP DWPT TCON TLCL TLCL1 SMITH.(329)
---- ---- ------ ------ ------ ------
40 20 15.473 15.449 15.450 15.48
-20 -45 -48.703 -48.768 -48.689 -48.79
10 -20 -25.237 -25.295 -29.401 -25.29
0 -25 -29.277 -29.331 -29.401 -29.31
30 10 5.708 5.682 5.679 5.72
40 35 33.710 33.727 33.636 33.74
EFFICIENCY TEST: 5000 CALCULATIONS AT T=20C, TD= 10C.
FUNC. T(SEC)
----- ------
TCON 0.65
TLCL 0.96
TLCL1 1.70
###############################################################################
TEST OF FUNCTION PCON
PURPOSE: CALCULATE LCL PRESSURE GIVEN INITIAL TEMPERATURE, DEW
POINT, AND LCL TEMPERATURE (TEML).
PRES TEMP TEML PCON
---- ---- ---- ------
1000 35 20 839.59
975 30 20 866.89
950 15 -10 691.24
900 10 -25 566.87
500 -20 -20 500.00
500 -20 -50 321.40
###############################################################################
TEST OF SUBROUTINE PTLCL, AND FUNCTION ALCL
PURPOSE: GIVEN TEMPERATURE, DEW POINT, AND PRESSURE... PLTLCL
ESTIMATES LCL TEMPERATURE AND PRESSURE; ALCL CALC-
ULATES LCL PRESSURE ONLY.
PRES TEMP DWPT PTLCL:P PTLCL:T ALCL
---- ---- ---- ------- ------- ------
1000 35 20 811.46 16.63 806.71
900 25 20 837.04 18.83 840.12
1000 30 28 971.57 27.51 973.34
850 0 -30 544.92 -34.66 528.01
900 16 16 900.00 16.00 900.00
700 -10 -60 333.50 -65.69 300.97
###############################################################################
TEST OF FUNCTION WOBF
PURPOSE: CALCULATE THE DIFFERENCE WBPTS-WBPTD (SEE NOTES IN
[SCHLATTER]THERMO.FOR FOR EXPLANATION) GIVEN THE TEMPER-
ATURE. IN THIS CASE, EQUIVALENT POTENTIAL TEMPERATURE
WAS INPUT SO RESULTS COULD BE CHECKED USING THE
SMITHSONIAN TABLES.
EQPT WOBF(=EQPT-WBPT) SMITH.(319)
---- --------------- ------
205.24 165.2104 165.24
161.04 125.0629 125.04
113.04 83.1266 83.04
70.34 48.8294 48.34
54.84 37.4419 36.84
42.14 28.6543 28.14
6.74 8.8176 8.74
###############################################################################
TEST OF FUNCTIONS POWT,OW
PURPOSE: CALCULATE THE WET-BULB POTENTIAL TEMPERATURE GIVEN THE
TEMPERATURE, DEW POINT, AND PRESSURE.
TEMP DWPT PRES POWT OW
---- ---- ---- ----- -----
35 25 1000 27.55 27.61
10 -15 1000 2.13 1.83
-10 -60 850 -4.64 -4.68
15 10 900 16.29 16.42
37 30 1050 29.98 30.13
-5 -5 500 22.05 22.22
EFFICIENCY TEST: 500 CALCULATIONS AT T=10C, TD=-10C, P=850MB.
FUNC. T(SEC)
----- ------
POWT 0.53
OW 17.31
###############################################################################
TEST OF FUNCTIONS DPT,DEWPT
PURPOSE: CALCULATE DEW POINT GIVEN SATURATION VAPOR PRESSURE.
SAT.VP. DPT DEWPT SMITH.(351)
------- ------ ------ ------
0.1891 -40.002 -40.025 -40
1.2540 -20.000 -20.032 -20
2.8627 -10.000 -10.022 -10
6.1078 0.000 -0.010 0
12.272 10.000 10.000 10
23.373 20.000 20.003 20
EFFICIENCY TEST: 3000 CALCULATIONS AT SVP=.1891MB.
FUNC. T(SEC)
----- ------
DPT 1.92
DEWPT 0.59
###############################################################################
TEST OF FUNCTION SATLFT
PURPOSE: CALCULATE TEMPERATURE WHERE A MOIST ADIABAT CROSSES A GIVEN
PRESSURE GIVEN THE VALUE OF THE MOIST ADIABAT AND THE
PRESSURE.
PRES THETAW SATLFT SMITH.
----- ------ ------- ------
204.3 34 -19.791 -20
200 20 -60.356 -62
1071.4 38 40.151 40
501 32 9.941 10
422.7 0 -51.876 -52
###############################################################################
TEST OF FUNCTION HUM
PURPOSE: CALCULATE RELATIVE HUMIDITY GIVEN TEMPERATURE AND
DEW POINT.
TEMP DWPT HUM
---- ---- -----
35 30 75.46
25 10 38.77
0 -15 31.18
20 -10 12.22
30 28 89.09
TEST OF FUNCTION TW
PURPOSE: CALCULATE WET-BULB TEMPERATURE GIVEN TEMPERATURE, DEW
POINT, AND PRESSURE.
TEMP DWPT PRES TW
---- ---- ---- -----
30 15 1000 19.99
0 -20 700 -6.05
10 0 850 5.11
-15 -25 500 -17.39
25 20 900 21.57
###############################################################################
TEST OF FUNCTIONS OE,EPT
PURPOSE: CALCULATE EQUIVALENT POTENTIAL TEMPERATURE GIVEN TEMP-
ERATURE, DEW POINT, AND PRESSURE.
TEMP DWPT PRES OE EPT
---- ---- ---- ----- -----
30 15 1000 62.24 62.49
0 -20 700 33.10 33.01
10 0 850 36.97 37.01
-15 -25 500 45.04 45.02
25 20 900 85.02 84.60
EFFICIENCY TEST: 1000 CALCULATIONS AT T=20C, TD=10C, P=1000MB.
FUNC. T(SEC)
----- ------
OE 18.50
EPT 0.80
###############################################################################
TEST OF FUNCTION TE
PURPOSE: CALCULATE EQUIVALENT TEMPERATURE GIVEN TEMPERATURE, DEW
POINT, AND PRESSURE.
TEMP DWPT PRES TE
---- ---- ---- -----
30 15 1000 62.24
0 -20 700 3.31
10 0 850 22.88
-15 -25 500 -12.18
25 20 900 74.38
###############################################################################
TEST OF FUNCTION TDA
PURPOSE: CALCULATE THE TEMPERATURE ON A DRY ADIABAT GIVEN THE VALUE
OF THE DRY ADIABAT AND THE PRESSURE.
THETA PRES TDA SMITH.(308)
----- ---- ------ ------
132.74 250 -0.12 0
35.44 500 -20.05 -20
-14.66 1050 -11.03 -11
28.64 850 14.93 -15
18.24 700 -10.02 -10
###############################################################################
TEST OF FUNCTIONS TSA,TMLAPS
PURPOSE: CALCULATE TEMPERATURE ON A MOIST ADIABAT GIVEN
EQUIVALENT POTENTIAL TEMPERATURE AND PRESSURE.
EQPT PRES TSA TMLAPS SMITH.(319)
------ ------ ------ ------- ------
181.64 870.9 34.170 33.940 34
113.04 577.0 12.139 12.080 12
70.34 896.0 18.057 17.950 18
10.24 496.0 -41.766 -41.669 -42
54.84 309.5 -39.658 -39.600 -40
EFFICIENCY TEST: 500 CALCULATIONS AT EQPT=113.04, P=577MB
FUNC. T(SEC)
----- ------
TSA 7.16
TMLAPS 6.81
###############################################################################
TEST OF FUNCTION O
PURPOSE: CALCULATE POTENTIAL TEMPERATURE GIVEN TEMPERATURE AND
PRESSURE.
TEMP PRES O SMITH.(309)
---- ---- ----- ------
-35 850 -23.66 -23.66
-10 500 47.74 47.74
-20 950 -16.26 -16.26
10 1000 10.04 10.00
###############################################################################
TEST OF FUNCTION OS
PURPOSE: CALCULATE EQUIVALENT POTENTIAL TEMPERATURE FOR AIR SAT-
URATED AT GIVEN TEMPERATURE AND PRESSURE.
TEMP PRES OS SMITH.(319)
---- ----- ------ -------
0 275.3 179.88 181.64
-10 321.4 112.03 113.04
10 501.0 126.27 127.34
10 816.0 54.89 54.84
-50 282.7 47.63 48.14
###############################################################################
TEST OF FUNCTION TMR
PURPOSE: CALCULATE THE TEMPERATURE ALONG A GIVEN MIXING RATIO
LINE AT A GIVEN PRESSURE.
W PRES TMR SMITH.(302)
----- ---- ----- ------
9.146 850 10.07 10
1.120 700 -19.96 -20
7.710 500 0.03 0
1.960 400 -19.97 -20
14.95 1000 20.10 20
###############################################################################
TEST OF FUNCTION THM
PURPOSE: CALCULATE WET-BULB POTENTIAL TEMPERATURE OF A PARCEL
SATURATED AT A GIVEN TEMPERATURE AND PRESSURE.
TEMP PRES THM SMITH.(319)
---- ------ ------ ------
12 369.6 46.277 40
-32 188.2 41.376 32
38 1069.8 35.751 36
10 908. 14.205 14
-14 648. 7.841 8
-50 301.7 13.961 14
###############################################################################
TEST OF FUNCTION DWPT
PURPOSE: CALCULATE DEW POINT GIVEN TEMPERATURE AND RELATIVE HUM-
IDITY.
TEMP REL.HUM. TD TD*
---- -------- ------ -----
35 75.46 30.14 30
25 38.77 9.93 10
0 31.18 -15.19 -15
20 12.22 -10.16 -10
30 89.09 28.01 28
* INPUT RELATIVE HUMIDITY FROM OUTPUT OF FUNCTION HUM.
EXPECT DEW POINT OUTPUT TO BE APPROX. DEW POINT INPUT
USED IN "HUM". SEE TEST OF FUNCTION HUM FOR COMPARISON.
###############################################################################
TEST OF FUNCTION SSH
PURPOSE: CALCULATE SATURATION SPECIFIC HUMIDITY GIVEN THE
PRESSURE AND TEMPERATURE.
PRES TEMP SSH !CHECKED ON HAND CALCULATOR.
---- ---- -----
1000 10 7.70
1000 0 3.82
850 -10 2.11
700 0 5.46
500 -20 1.57
300 -30 1.06
###############################################################################
TEST OF FUNCTION TV
PURPOSE: CALCULATE VIRTUAL TEMPERATURE GIVEN TEMPERATURE, DEW POINT,
AND PRESSURE.
TEMP DWPT PRES TV !CHECKED ON HAND CALCULATOR.
---- ---- ---- -----
30 25 1000 33.69
20 0 1000 20.68
10 -10 850 10.36
0 -10 700 0.42
-20 -30 500 -19.90
-40 -60 300 -39.99
###############################################################################
TEST OF FUNCTIONS ESICE,ESILO
PURPOSE: CALCULATE SATURATION VAPOR PRESSURE OVER ICE GIVEN TEMP-
ERATURE.
TEMP ESICE ESILO SMITH.(360)
---- ------- ------- ------
-50 .03935 .03963 .03935
-40 .1283 .1283 .1283
-30 .3798 .3796 .3798
-20 1.032 1.032 1.032
-10 2.597 2.596 2.597
EFFICIENCY TEST: 2000 CALCULATIONS AT T=-10C.
FUNC. T(SEC)
----- ------
ESICE 1.10
ESILO 0.39
###############################################################################
TEST OF FUNCTION HEATL
PURPOSE: GIVEN TEMPERATURE, TO CALCULATE.....
1) LATENT HEAT OF EVAPORATION/CONDENSATION (KEY=1)
2) LATENT HEAT OF FREEZING/MELTING (KEY=2)
3) LATENT HEAT OF SUBLIMATION/DEPOSITION (KEY=3)
(RESULTS VERIFIED USING P.343, SMITH.)
TEMP KEY=1 SMITH. KEY=2 SMITH. KEY=3 SMITH.
---- -------- ------- -------- ------ --------- ------
-100 674.2192 674.4
-80 676.1257 676.3
-60 677.3396 677.5
-50 628.1667 629.3 48.7133 48.6
-40 56.1209 56.3 677.8610 678.0
-30 615.5021 615.0 62.9107 63.0
-20 69.0828 69.0 677.6898 677.9
-10 603.2438 603.0 74.6372 74.5 677.3444 677.5
0 597.2670 597.3 79.5740 79.7 676.8259 677.0
10 591.3918 591.7
20 585.6181 586.0
40 574.3755 574.7
50 568.9066 569.0
NOTE: HEATL IS DESIGNED TO RETURN THE LATENT HEATS IN UNITS OF
JOULES/KG. TO MAKE RESULTS COMPATIBLE WITH VALUES IN THE
SMITHSONIAN TABLES (IT-CAL/GRAM), ANSWERS WERE MULTIPLIED
BY THE CORRECTION FACTOR 2.3884E-04 TO OBTAIN THOSE SHOWN
ABOVE.
###############################################################################
USER_DEV:[GULIB.THERMOSRC]ALCL_TG.FOR;1
FUNCTION ALCL(T,TD,P)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C G.S. Stipanuk 1973 Original version.
C Reference Stipanuk paper entitled:
C "ALGORITHMS FOR GENERATING A SKEW-T, LOG P
C DIAGRAM AND COMPUTING SELECTED METEOROLOGICAL
C QUANTITIES."
C ATMOSPHERIC SCIENCES LABORATORY
C U.S. ARMY ELECTRONICS COMMAND
C WHITE SANDS MISSILE RANGE, NEW MEXICO 88002
C 33 PAGES
C Baker, Schlatter 17-MAY-1982
C THIS FUNCTION RETURNS THE PRESSURE ALCL (MB) OF THE LIFTING CONDEN-
C SATION LEVEL (LCL) FOR A PARCEL INITIALLY AT TEMPERATURE T (CELSIUS)
C DEW POINT TD (CELSIUS) AND PRESSURE P (MILLIBARS). ALCL IS COMPUTED
C BY AN ITERATIVE PROCEDURE DESCRIBED BY EQS. 8-12 IN STIPANUK (1973),
C PP.13-14.
C DETERMINE THE MIXING RATIO LINE THROUGH TD AND P.
AW = W(TD,P)
C DETERMINE THE DRY ADIABAT THROUGH T AND P.
AO = O(T,P)
C ITERATE TO LOCATE PRESSURE PI AT THE INTERSECTION OF THE TWO
C CURVES. PI HAS BEEN SET TO P FOR THE INITIAL GUESS.
3 CONTINUE
PI = P
DO 4 I= 1,10
X= .02*(TMR(AW,PI)-TDA(AO,PI))
IF (ABS(X).LT.0.01) GO TO 5
4 PI= PI*(2.**(X))
5 ALCL= PI
RETURN
ENTRY ALCLM(T,TD,P)
C FOR ENTRY ALCLM ONLY, T IS THE MEAN POTENTIAL TEMPERATURE (CELSIUS)
C AND TD IS THE MEAN MIXING RATIO (G/KG) OF THE LAYER CONTAINING THE
C PARCEL.
AW = TD
AO = T
GO TO 3
END
USER_DEV:[GULIB.THERMOSRC]CT_TG.FOR;1
FUNCTION CT(WBAR,PC,PS)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C G.S. Stipanuk 1973 Original version.
C Reference Stipanuk paper entitled:
C "ALGORITHMS FOR GENERATING A SKEW-T, LOG P
C DIAGRAM AND COMPUTING SELECTED METEOROLOGICAL
C QUANTITIES."
C ATMOSPHERIC SCIENCES LABORATORY
C U.S. ARMY ELECTRONICS COMMAND
C WHITE SANDS MISSILE RANGE, NEW MEXICO 88002
C 33 PAGES
C Baker, Schlatter 17-MAY-1982
C THIS FUNCTION RETURNS THE CONVECTIVE TEMPERATURE CT (CELSIUS)
C GIVEN THE MEAN MIXING RATIO WBAR (G/KG) IN THE SURFACE LAYER,
C THE PRESSURE PC (MB) AT THE CONVECTIVE CONDENSATION LEVEL (CCL)
C AND THE SURFACE PRESSURE PS (MB).
C COMPUTE THE TEMPERATURE (CELSIUS) AT THE CCL.
TC= TMR(WBAR,PC)
C COMPUTE THE POTENTIAL TEMPERATURE (CELSIUS), I.E., THE DRY
C ADIABAT AO THROUGH THE CCL.
AO= O(TC,PC)
C COMPUTE THE SURFACE TEMPERATURE ON THE SAME DRY ADIABAT AO.
CT= TDA(AO,PS)
RETURN
END
USER_DEV:[GULIB.THERMOSRC]DEWPT_TG.FOR;1
FUNCTION DEWPT(EW)
C THIS FUNCTION YIELDS THE DEW POINT DEWPT (CELSIUS), GIVEN THE
C WATER VAPOR PRESSURE EW (MILLIBARS).
C THE EMPIRICAL FORMULA APPEARS IN BOLTON, DAVID, 1980:
C "THE COMPUTATION OF EQUIVALENT POTENTIAL TEMPERATURE,"
C MONTHLY WEATHER REVIEW, VOL. 108, NO. 7 (JULY), P. 1047, EQ.(11).
C THE QUOTED ACCURACY IS 0.03C OR LESS FOR -35 < DEWPT < 35C.
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
ENL = ALOG(EW)
DEWPT = (243.5*ENL-440.8)/(19.48-ENL)
RETURN
END
USER_DEV:[GULIB.THERMOSRC]DPT_TG.FOR;1
FUNCTION DPT(EW)
C THIS FUNCTION RETURNS THE DEW POINT DPT (CELSIUS), GIVEN THE
C WATER VAPOR PRESSURE EW (MILLIBARS).
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
DATA ES0/6.1078/
C ES0 = SATURATION VAPOR PRESSURE (MB) OVER WATER AT 0C
C RETURN A FLAG VALUE IF THE VAPOR PRESSURE IS OUT OF RANGE.
IF (EW.GT..06.AND.EW.LT.1013.) GO TO 5
DPT = 9999.
RETURN
5 CONTINUE
C APPROXIMATE DEW POINT BY MEANS OF TETEN'S FORMULA.
C THE FORMULA APPEARS AS EQ.(8) IN BOLTON, DAVID, 1980:
C "THE COMPUTATION OF EQUIVALENT POTENTIAL TEMPERATURE,"
C MONTHLY WEATHER REVIEW, VOL 108, NO. 7 (JULY), P.1047.
C THE FORMULA IS EW(T) = ES0*10**(7.5*T/(T+237.3))
C OR EW(T) = ES0*EXP(17.269388*T/(T+237.3))
C THE INVERSE FORMULA IS USED BELOW.
X = ALOG(EW/ES0)
DNM = 17.269388-X
T = 237.3*X/DNM
FAC = 1./(EW*DNM)
C LOOP FOR ITERATIVE IMPROVEMENT OF THE ESTIMATE OF DEW POINT
10 CONTINUE
C GET THE PRECISE VAPOR PRESSURE CORRESPONDING TO T.
EDP = ESW(T)
C ESTIMATE THE CHANGE IN TEMPERATURE CORRESPONDING TO (EW-EDP)
C ASSUME THAT THE DERIVATIVE OF TEMPERATURE WITH RESPECT TO
C VAPOR PRESSURE (DTDEW) IS GIVEN BY THE DERIVATIVE OF THE
C INVERSE TETEN FORMULA.
DTDEW = (T+237.3)*FAC
DT = DTDEW*(EW-EDP)
T = T+DT
IF (ABS(DT).GT.1.E-04) GO TO 10
DPT = T
RETURN
END
USER_DEV:[GULIB.THERMOSRC]DWPT_TG.FOR;1
FUNCTION DWPT(T,RH)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C THIS FUNCTION RETURNS THE DEW POINT (CELSIUS) GIVEN THE TEMPERATURE
C (CELSIUS) AND RELATIVE HUMIDITY (%). THE FORMULA IS USED IN THE
C PROCESSING OF U.S. RAWINSONDE DATA AND IS REFERENCED IN PARRY, H.
C DEAN, 1969: "THE SEMIAUTOMATIC COMPUTATION OF RAWINSONDES,"
C TECHNICAL MEMORANDUM WBTM EDL 10, U.S. DEPARTMENT OF COMMERCE,
C ENVIRONMENTAL SCIENCE SERVICES ADMINISTRATION, WEATHER BUREAU,
C OFFICE OF SYSTEMS DEVELOPMENT, EQUIPMENT DEVELOPMENT LABORATORY,
C SILVER SPRING, MD (OCTOBER), PAGE 9 AND PAGE II-4, LINE 460.
X = 1.-0.01*RH
C COMPUTE DEW POINT DEPRESSION.
DPD =(14.55+0.114*T)*X+((2.5+0.007*T)*X)**3+(15.9+0.117*T)*X**14
DWPT = T-DPD
RETURN
END
USER_DEV:[GULIB.THERMOSRC]EPT_TG.FOR;1
FUNCTION EPT(T,TD,P)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C THIS FUNCTION RETURNS THE EQUIVALENT POTENTIAL TEMPERATURE EPT
C (CELSIUS) FOR A PARCEL OF AIR INITIALLY AT TEMPERATURE T (CELSIUS),
C DEW POINT TD (CELSIUS) AND PRESSURE P (MILLIBARS). THE FORMULA USED
C IS EQ.(43) IN BOLTON, DAVID, 1980: "THE COMPUTATION OF EQUIVALENT
C POTENTIAL TEMPERATURE," MONTHLY WEATHER REVIEW, VOL. 108, NO. 7
C (JULY), PP. 1046-1053. THE MAXIMUM ERROR IN EPT IN 0.3C. IN MOST
C CASES THE ERROR IS LESS THAN 0.1C.
C
C COMPUTE THE MIXING RATIO (GRAMS OF WATER VAPOR PER KILOGRAM OF
C DRY AIR).
W = WMR(P,TD)
C COMPUTE THE TEMPERATURE (CELSIUS) AT THE LIFTING CONDENSATION LEVEL.
TLCL = TCON(T,TD)
TK = T+273.15
TL = TLCL+273.15
PT = TK*(1000./P)**(0.2854*(1.-0.00028*W))
EPTK = PT*EXP((3.376/TL-0.00254)*W*(1.+0.00081*W))
EPT= EPTK-273.15
RETURN
END
USER_DEV:[GULIB.THERMOSRC]ESAT_TG.FOR;1
FUNCTION ESAT(T)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C G.S. Stipanuk 1973 Original version.
C Reference Stipanuk paper entitled:
C "ALGORITHMS FOR GENERATING A SKEW-T, LOG P
C DIAGRAM AND COMPUTING SELECTED METEOROLOGICAL
C QUANTITIES."
C ATMOSPHERIC SCIENCES LABORATORY
C U.S. ARMY ELECTRONICS COMMAND
C WHITE SANDS MISSILE RANGE, NEW MEXICO 88002
C 33 PAGES
C Baker, Schlatter 17-MAY-1982
C THIS FUNCTION RETURNS THE SATURATION VAPOR PRESSURE OVER
C WATER (MB) GIVEN THE TEMPERATURE (CELSIUS).
C THE ALGORITHM IS DUE TO NORDQUIST, W.S.,1973: "NUMERICAL APPROXIMA-
C TIONS OF SELECTED METEORLOLGICAL PARAMETERS FOR CLOUD PHYSICS PROB-
C LEMS," ECOM-5475, ATMOSPHERIC SCIENCES LABORATORY, U.S. ARMY
C ELECTRONICS COMMAND, WHITE SANDS MISSILE RANGE, NEW MEXICO 88002.
TK = T+273.15
P1 = 11.344-0.0303998*TK
P2 = 3.49149-1302.8844/TK
C1 = 23.832241-5.02808*ALOG10(TK)
ESAT = 10.**(C1-1.3816E-7*10.**P1+8.1328E-3*10.**P2-2949.076/TK)
RETURN
END
USER_DEV:[GULIB.THERMOSRC]ESGG_TG.FOR;1
FUNCTION ESGG(T)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C THIS FUNCTION RETURNS THE SATURATION VAPOR PRESSURE OVER LIQUID
C WATER ESGG (MILLIBARS) GIVEN THE TEMPERATURE T (CELSIUS). THE
C FORMULA USED, DUE TO GOFF AND GRATCH, APPEARS ON P. 350 OF THE
C SMITHSONIAN METEOROLOGICAL TABLES, SIXTH REVISED EDITION, 1963,
C BY ROLAND LIST.
DATA CTA,EWS,TS/273.15,1013.246,373.15/
C CTA = DIFFERENCE BETWEEN KELVIN AND CELSIUS TEMPERATURES
C EWS = SATURATION VAPOR PRESSURE (MB) OVER LIQUID WATER AT 100C
C TS = BOILING POINT OF WATER (K)
DATA C1, C2, C3, C4, C5, C6
1 / 7.90298, 5.02808, 1.3816E-7, 11.344, 8.1328E-3, 3.49149 /
TK = T+CTA
C GOFF-GRATCH FORMULA
RHS = -C1*(TS/TK-1.)+C2*ALOG10(TS/TK)-C3*(10.**(C4*(1.-TK/TS))
1 -1.)+C5*(10.**(-C6*(TS/TK-1.))-1.)+ALOG10(EWS)
ESW = 10.**RHS
IF (ESW.LT.0.) ESW = 0.
ESGG = ESW
RETURN
END
USER_DEV:[GULIB.THERMOSRC]ESICE_TG.FOR;2
FUNCTION ESICE(T)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C THIS FUNCTION RETURNS THE SATURATION VAPOR PRESSURE WITH RESPECT TO
C ICE ESICE (MILLIBARS) GIVEN THE TEMPERATURE T (CELSIUS).
C THE FORMULA USED IS BASED UPON THE INTEGRATION OF THE CLAUSIUS-
C CLAPEYRON EQUATION BY GOFF AND GRATCH. THE FORMULA APPEARS ON P.350
C OF THE SMITHSONIAN METEOROLOGICAL TABLES, SIXTH REVISED EDITION,
C 1963.
DATA CTA,EIS/273.15,6.1071/
C CTA = DIFFERENCE BETWEEN KELVIN AND CELSIUS TEMPERATURE
C EIS = SATURATION VAPOR PRESSURE (MB) OVER A WATER-ICE MIXTURE AT 0C
DATA C1,C2,C3/9.09718,3.56654,0.876793/
C C1,C2,C3 = EMPIRICAL COEFFICIENTS IN THE GOFF-GRATCH FORMULA
IF (T.LE.0.) GO TO 5
ESICE = 99999.
WRITE(6,3)ESICE
UNLOCK (6)
3 FORMAT(' SATURATION VAPOR PRESSURE FOR ICE CANNOT BE COMPUTED',
1 /' FOR TEMPERATURE > 0C. ESICE =',F7.0)
RETURN
5 CONTINUE
C FREEZING POINT OF WATER (K)
TF = CTA
TK = T+CTA
C GOFF-GRATCH FORMULA
RHS = -C1*(TF/TK-1.)-C2*ALOG10(TF/TK)+C3*(1.-TK/TF)+ALOG10(EIS)
ESI = 10.**RHS
IF (ESI.LT.0.) ESI = 0.
ESICE = ESI
RETURN
END
USER_DEV:[GULIB.THERMOSRC]ESILO_TG.FOR;2
FUNCTION ESILO(T)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C THIS FUNCTION RETURNS THE SATURATION VAPOR PRESSURE OVER ICE
C ESILO (MILLIBARS) GIVEN THE TEMPERATURE T (CELSIUS). THE FORMULA
C IS DUE TO LOWE, PAUL R., 1977: AN APPROXIMATING POLYNOMIAL FOR
C THE COMPUTATION OF SATURATION VAPOR PRESSURE, JOURNAL OF APPLIED
C METEOROLOGY, VOL. 16, NO. 1 (JANUARY), PP. 100-103.
C THE POLYNOMIAL COEFFICIENTS ARE A0 THROUGH A6.
DATA A0,A1,A2,A3,A4,A5,A6
1 /6.109177956, 5.034698970E-01, 1.886013408E-02,
2 4.176223716E-04, 5.824720280E-06, 4.838803174E-08,
3 1.838826904E-10/
IF (T.LE.0.) GO TO 5
ESILO = 9999.
WRITE(6,3)ESILO
UNLOCK (6)
3 FORMAT(' SATURATION VAPOR PRESSURE OVER ICE IS UNDEFINED FOR',
1 /' TEMPERATURE > 0C. ESILO =',F6.0)
RETURN
5 CONTINUE
ESILO = A0+T*(A1+T*(A2+T*(A3+T*(A4+T*(A5+A6*T)))))
RETURN
END
USER_DEV:[GULIB.THERMOSRC]ESLO_TG.FOR;1
FUNCTION ESLO(T)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C THIS FUNCTION RETURNS THE SATURATION VAPOR PRESSURE OVER LIQUID
C WATER ESLO (MILLIBARS) GIVEN THE TEMPERATURE T (CELSIUS). THE
C FORMULA IS DUE TO LOWE, PAUL R.,1977: AN APPROXIMATING POLYNOMIAL
C FOR THE COMPUTATION OF SATURATION VAPOR PRESSURE, JOURNAL OF APPLIED
C METEOROLOGY, VOL 16, NO. 1 (JANUARY), PP. 100-103.
C THE POLYNOMIAL COEFFICIENTS ARE A0 THROUGH A6.
DATA A0,A1,A2,A3,A4,A5,A6
1 /6.107799961, 4.436518521E-01, 1.428945805E-02,
2 2.650648471E-04, 3.031240396E-06, 2.034080948E-08,
3 6.136820929E-11/
ES = A0+T*(A1+T*(A2+T*(A3+T*(A4+T*(A5+A6*T)))))
IF (ES.LT.0.) ES = 0.
ESLO = ES
RETURN
END
USER_DEV:[GULIB.THERMOSRC]ESRW_TG.FOR;1
FUNCTION ESRW(T)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C THIS FUNCTION RETURNS THE SATURATION VAPOR PRESSURE OVER LIQUID
C WATER ESRW (MILLIBARS) GIVEN THE TEMPERATURE T (CELSIUS). THE
C FORMULA USED IS DUE TO RICHARDS, J.M., 1971: SIMPLE EXPRESSION
C FOR THE SATURATION VAPOUR PRESSURE OF WATER IN THE RANGE -50 TO
C 140C, BRITISH JOURNAL OF APPLIED PHYSICS, VOL. 4, PP.L15-L18.
C THE FORMULA WAS QUOTED MORE RECENTLY BY WIGLEY, T.M.L.,1974:
C COMMENTS ON 'A SIMPLE BUT ACCURATE FORMULA FOR THE SATURATION
C VAPOR PRESSURE OVER LIQUID WATER,' JOURNAL OF APPLIED METEOROLOGY,
C VOL. 13, NO. 5 (AUGUST) P.606.
DATA CTA,TS,EWS/273.15,373.15,1013.25/
C CTA = DIFFERENCE BETWEEN KELVIN AND CELSIUS TEMPERATURE
C TS = TEMPERATURE OF THE BOILING POINT OF WATER (K)
C EWS = SATURATION VAPOR PRESSURE OVER LIQUID WATER AT 100C
DATA C1, C2, C3, C4
1 / 13.3185,-1.9760,-0.6445,-0.1299 /
TK = T+CTA
X = 1.-TS/TK
PX = X*(C1+X*(C2+X*(C3+C4*X)))
VP = EWS*EXP(PX)
IF (VP.LT.0) VP = 0.
ESRW = VP
RETURN
END
USER_DEV:[GULIB.THERMOSRC]ESW_TG.FOR;1
FUNCTION ESW(T)
C THIS FUNCTION RETURNS THE SATURATION VAPOR PRESSURE ESW (MILLIBARS)
C OVER LIQUID WATER GIVEN THE TEMPERATURE T (CELSIUS). THE POLYNOMIAL
C APPROXIMATION BELOW IS DUE TO HERMAN WOBUS, A MATHEMATICIAN WHO
C WORKED AT THE NAVY WEATHER RESEARCH FACILITY, NORFOLK, VIRGINIA,
C BUT WHO IS NOW RETIRED. THE COEFFICIENTS OF THE POLYNOMIAL WERE
C CHOSEN TO FIT THE VALUES IN TABLE 94 ON PP. 351-353 OF THE SMITH-
C SONIAN METEOROLOGICAL TABLES BY ROLAND LIST (6TH EDITION). THE
C APPROXIMATION IS VALID FOR -50 < T < 100C.
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C
C ES0 = SATURATION VAPOR RESSURE OVER LIQUID WATER AT 0C
DATA ES0/6.1078/
POL = 0.99999683 + T*(-0.90826951E-02 +
1 T*(0.78736169E-04 + T*(-0.61117958E-06 +
2 T*(0.43884187E-08 + T*(-0.29883885E-10 +
3 T*(0.21874425E-12 + T*(-0.17892321E-14 +
4 T*(0.11112018E-16 + T*(-0.30994571E-19)))))))))
ESW = ES0/POL**8
RETURN
END
USER_DEV:[GULIB.THERMOSRC]ES_TG.FOR;1
FUNCTION ES(T)
C THIS FUNCTION RETURNS THE SATURATION VAPOR PRESSURE ES (MB) OVER
C LIQUID WATER GIVEN THE TEMPERATURE T (CELSIUS). THE FORMULA APPEARS
C IN BOLTON, DAVID, 1980: "THE COMPUTATION OF EQUIVALENT POTENTIAL
C TEMPERATURE," MONTHLY WEATHER REVIEW, VOL. 108, NO. 7 (JULY),
C P. 1047, EQ.(10). THE QUOTED ACCURACY IS 0.3% OR BETTER FOR
C -35 < T < 35C.
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C ES0 = SATURATION VAPOR PRESSURE OVER LIQUID WATER AT 0C
DATA ES0/6.1121/
ES = ES0*EXP(17.67*T/(T+243.5))
RETURN
END
USER_DEV:[GULIB.THERMOSRC]HEATL_TG.FOR;1
FUNCTION HEATL(KEY,T)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C THIS FUNCTION RETURNS THE LATENT HEAT OF
C EVAPORATION/CONDENSATION FOR KEY=1
C MELTING/FREEZING FOR KEY=2
C SUBLIMATION/DEPOSITION FOR KEY=3
C FOR WATER. THE LATENT HEAT HEATL (JOULES PER KILOGRAM) IS A
C FUNCTION OF TEMPERATURE T (CELSIUS). THE FORMULAS ARE POLYNOMIAL
C APPROXIMATIONS TO THE VALUES IN TABLE 92, P. 343 OF THE SMITHSONIAN
C METEOROLOGICAL TABLES, SIXTH REVISED EDITION, 1963 BY ROLAND LIST.
C THE APPROXIMATIONS WERE DEVELOPED BY ERIC SMITH AT COLORADO STATE
C UNIVERSITY.
C POLYNOMIAL COEFFICIENTS
DATA A0,A1,A2/ 3337118.5,-3642.8583, 2.1263947/
DATA B0,B1,B2/-1161004.0, 9002.2648,-12.931292/
DATA C0,C1,C2/ 2632536.8, 1726.9659,-3.6248111/
HLTNT = 0.
TK = T+273.15
IF (KEY.EQ.1) HLTNT=A0+A1*TK+A2*TK*TK
IF (KEY.EQ.2) HLTNT=B0+B1*TK+B2*TK*TK
IF (KEY.EQ.3) HLTNT=C0+C1*TK+C2*TK*TK
HEATL = HLTNT
RETURN
END
USER_DEV:[GULIB.THERMOSRC]HUM_TG.FOR;1
FUNCTION HUM(T,TD)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C G.S. Stipanuk 1973 Original version.
C Reference Stipanuk paper entitled:
C "ALGORITHMS FOR GENERATING A SKEW-T, LOG P
C DIAGRAM AND COMPUTING SELECTED METEOROLOGICAL
C QUANTITIES."
C ATMOSPHERIC SCIENCES LABORATORY
C U.S. ARMY ELECTRONICS COMMAND
C WHITE SANDS MISSILE RANGE, NEW MEXICO 88002
C 33 PAGES
C Baker, Schlatter 17-MAY-1982
C THIS FUNCTION RETURNS RELATIVE HUMIDITY (%) GIVEN THE
C TEMPERATURE T AND DEW POINT TD (CELSIUS). AS CALCULATED HERE,
C RELATIVE HUMIDITY IS THE RATIO OF THE ACTUAL VAPOR PRESSURE TO
C THE SATURATION VAPOR PRESSURE.
HUM= 100.*(ESAT(TD)/ESAT(T))
RETURN
END
USER_DEV:[GULIB.THERMOSRC]OE_TG.FOR;1
FUNCTION OE(T,TD,P)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C G.S. Stipanuk 1973 Original version.
C Reference Stipanuk paper entitled:
C "ALGORITHMS FOR GENERATING A SKEW-T, LOG P
C DIAGRAM AND COMPUTING SELECTED METEOROLOGICAL
C QUANTITIES."
C ATMOSPHERIC SCIENCES LABORATORY
C U.S. ARMY ELECTRONICS COMMAND
C WHITE SANDS MISSILE RANGE, NEW MEXICO 88002
C 33 PAGES
C Baker, Schlatter 17-MAY-1982
C THIS FUNCTION RETURNS EQUIVALENT POTENTIAL TEMPERATURE OE (CELSIUS)
C OF A PARCEL OF AIR GIVEN ITS TEMPERATURE T (CELSIUS), DEW POINT
C TD (CELSIUS) AND PRESSURE P (MILLIBARS).
C FIND THE WET BULB TEMPERATURE OF THE PARCEL.
ATW = TW(T,TD,P)
C FIND THE EQUIVALENT POTENTIAL TEMPERATURE.
OE = OS(ATW,P)
RETURN
END
USER_DEV:[GULIB.THERMOSRC]OS_TG.FOR;1
FUNCTION OS(T,P)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C G.S. Stipanuk 1973 Original version.
C Reference Stipanuk paper entitled:
C "ALGORITHMS FOR GENERATING A SKEW-T, LOG P
C DIAGRAM AND COMPUTING SELECTED METEOROLOGICAL
C QUANTITIES."
C ATMOSPHERIC SCIENCES LABORATORY
C U.S. ARMY ELECTRONICS COMMAND
C WHITE SANDS MISSILE RANGE, NEW MEXICO 88002
C 33 PAGES
C Baker, Schlatter 17-MAY-1982
C THIS FUNCTION RETURNS THE EQUIVALENT POTENTIAL TEMPERATURE OS
C (CELSIUS) FOR A PARCEL OF AIR SATURATED AT TEMPERATURE T (CELSIUS)
C AND PRESSURE P (MILLIBARS).
DATA B/2.6518986/
C B IS AN EMPIRICAL CONSTANT APPROXIMATELY EQUAL TO THE LATENT HEAT
C OF VAPORIZATION FOR WATER DIVIDED BY THE SPECIFIC HEAT AT CONSTANT
C PRESSURE FOR DRY AIR.
TK = T+273.15
OSK= TK*((1000./P)**.286)*(EXP(B*W(T,P)/TK))
OS= OSK-273.15
RETURN
END
USER_DEV:[GULIB.THERMOSRC]OW_TG.FOR;1
FUNCTION OW(T,TD,P)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C G.S. Stipanuk 1973 Original version.
C Reference Stipanuk paper entitled:
C "ALGORITHMS FOR GENERATING A SKEW-T, LOG P
C DIAGRAM AND COMPUTING SELECTED METEOROLOGICAL
C QUANTITIES."
C ATMOSPHERIC SCIENCES LABORATORY
C U.S. ARMY ELECTRONICS COMMAND
C WHITE SANDS MISSILE RANGE, NEW MEXICO 88002
C 33 PAGES
C Baker, Schlatter 17-MAY-1982
C THIS FUNCTION RETURNS THE WET-BULB POTENTIAL TEMPERATURE OW
C (CELSIUS) GIVEN THE TEMPERATURE T (CELSIUS), DEW POINT TD
C (CELSIUS), AND PRESSURE P (MILLIBARS). THE CALCULATION FOR OW IS
C VERY SIMILAR TO THAT FOR WET BULB TEMPERATURE. SEE P.13 STIPANUK (1973).
C FIND THE WET BULB TEMPERATURE OF THE PARCEL
ATW = TW(T,TD,P)
C FIND THE EQUIVALENT POTENTIAL TEMPERATURE OF THE PARCEL.
AOS= OS(ATW,P)
C FIND THE WET-BULB POTENTIAL TEMPERATURE.
OW= TSA(AOS,1000.)
RETURN
END
USER_DEV:[GULIB.THERMOSRC]O_TG.FOR;1
FUNCTION O(T,P)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C G.S. Stipanuk 1973 Original version.
C Reference Stipanuk paper entitled:
C "ALGORITHMS FOR GENERATING A SKEW-T, LOG P
C DIAGRAM AND COMPUTING SELECTED METEOROLOGICAL
C QUANTITIES."
C ATMOSPHERIC SCIENCES LABORATORY
C U.S. ARMY ELECTRONICS COMMAND
C WHITE SANDS MISSILE RANGE, NEW MEXICO 88002
C 33 PAGES
C Baker, Schlatter 17-MAY-1982
C THIS FUNCTION RETURNS POTENTIAL TEMPERATURE (CELSIUS) GIVEN
C TEMPERATURE T (CELSIUS) AND PRESSURE P (MB) BY SOLVING THE POISSON
C EQUATION.
TK= T+273.15
OK= TK*((1000./P)**.286)
O= OK-273.15
RETURN
END
USER_DEV:[GULIB.THERMOSRC]PCCL_TG.FOR;1
FUNCTION PCCL(PM,P,T,TD,WBAR,N)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C G.S. Stipanuk 1973 Original version.
C Reference Stipanuk paper entitled:
C "ALGORITHMS FOR GENERATING A SKEW-T, LOG P
C DIAGRAM AND COMPUTING SELECTED METEOROLOGICAL
C QUANTITIES."
C ATMOSPHERIC SCIENCES LABORATORY
C U.S. ARMY ELECTRONICS COMMAND
C WHITE SANDS MISSILE RANGE, NEW MEXICO 88002
C 33 PAGES
C Baker, Schlatter 17-MAY-1982
C THIS FUNCTION RETURNS THE PRESSURE AT THE CONVECTIVE CONDENSATION
C LEVEL GIVEN THE APPROPRIATE SOUNDING DATA.
C ON INPUT:
C P = PRESSURE (MILLIBARS). NOTE THAT P(I).GT.P(I+1).
C T = TEMPERATURE (CELSIUS)
C TD = DEW POINT (CELSIUS)
C N = NUMBER OF LEVELS IN THE SOUNDING AND THE DIMENSION OF
C P, T AND TD
C PM = PRESSURE (MILLIBARS) AT UPPER BOUNDARY OF THE LAYER FOR
C COMPUTING THE MEAN MIXING RATIO. P(1) IS THE LOWER
C BOUNDARY.
C ON OUTPUT:
C PCCL = PRESSURE (MILLIBARS) AT THE CONVECTIVE CONDENSATION LEVEL
C WBAR = MEAN MIXING RATIO (G/KG) IN THE LAYER BOUNDED BY
C PRESSURES P(1) AT THE BOTTOM AND PM AT THE TOP
C THE ALGORITHM IS DECRIBED ON P.17 OF STIPANUK, G.S.,1973:
C "ALGORITHMS FOR GENERATING A SKEW-T LOG P DIAGRAM AND COMPUTING
C SELECTED METEOROLOGICAL QUANTITIES," ATMOSPHERIC SCIENCES LABORA-
C TORY, U.S. ARMY ELECTRONICS COMMAND, WHITE SANDS MISSILE RANGE, NEW
C MEXICO 88002.
DIMENSION T(1),P(1),TD(1)
IF (PM.NE.P(1)) GO TO 5
WBAR= W(TD(1),P(1))
PC= PM
IF (ABS(T(1)-TD(1)).LT.0.05) GO TO 45
GO TO 25
5 CONTINUE
WBAR= 0.
K= 0
10 CONTINUE
K = K+1
IF (P(K).GT.PM) GO TO 10
K= K-1
J= K-1
IF(J.LT.1) GO TO 20
C COMPUTE THE AVERAGE MIXING RATIO....ALOG = NATURAL LOG
DO 15 I= 1,J
L= I+1
15 WBAR= (W(TD(I),P(I))+W(TD(L),P(L)))*ALOG(P(I)/P(L))
* +WBAR
20 CONTINUE
L= K+1
C ESTIMATE THE DEW POINT AT PRESSURE PM.
TQ= TD(K)+(TD(L)-TD(K))*(ALOG(PM/P(K)))/(ALOG(P(L)/P(K)))
WBAR= WBAR+(W(TD(K),P(K))+W(TQ,PM))*ALOG(P(K)/PM)
WBAR= WBAR/(2.*ALOG(P(1)/PM))
C FIND LEVEL AT WHICH THE MIXING RATIO LINE WBAR CROSSES THE
C ENVIRONMENTAL TEMPERATURE PROFILE.
25 CONTINUE
DO 30 J= 1,N
I= N-J+1
IF (P(I).LT.300.) GO TO 30
C TMR = TEMPERATURE (CELSIUS) AT PRESSURE P (MB) ALONG A MIXING
C RATIO LINE GIVEN BY WBAR (G/KG)
X= TMR(WBAR,P(I))-T(I)
IF (X.LE.0.) GO TO 35
30 CONTINUE
PCCL= 0.0
RETURN
C SET UP BISECTION ROUTINE
35 L = I
I= I+1
DEL= P(L)-P(I)
PC= P(I)+.5*DEL
A= (T(I)-T(L))/ALOG(P(L)/P(I))
DO 40 J= 1,10
DEL= DEL/2.
X= TMR(WBAR,PC)-T(L)-A*(ALOG(P(L)/PC))
C THE SIGN FUNCTION REPLACES THE SIGN OF THE FIRST ARGUMENT
C WITH THAT OF THE SECOND.
40 PC= PC+SIGN(DEL,X)
45 PCCL = PC
RETURN
END
USER_DEV:[GULIB.THERMOSRC]PCON_TG.FOR;1
FUNCTION PCON(P,T,TC)
C THIS FUNCTION RETURNS THE PRESSURE PCON (MB) AT THE LIFTED CONDENSA-
C TION LEVEL, GIVEN THE INITIAL PRESSURE P (MB) AND TEMPERATURE T
C (CELSIUS) OF THE PARCEL AND THE TEMPERATURE TC (CELSIUS) AT THE
C LCL. THE ALGORITHM IS EXACT. IT MAKES USE OF THE FORMULA FOR THE
C POTENTIAL TEMPERATURES CORRESPONDING TO T AT P AND TC AT PCON.
C THESE TWO POTENTIAL TEMPERATURES ARE EQUAL.
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C
DATA AKAPI/3.5037/
C AKAPI = (SPECIFIC HEAT AT CONSTANT PRESSURE FOR DRY AIR) /
C (GAS CONSTANT FOR DRY AIR)
C CONVERT T AND TC TO KELVIN TEMPERATURES.
TK = T+273.15
TCK = TC+273.15
PCON = P*(TCK/TK)**AKAPI
RETURN
END
USER_DEV:[GULIB.THERMOSRC]POWT_TG.FOR;1
FUNCTION POWT(T,P,TD)
C THIS FUNCTION YIELDS WET-BULB POTENTIAL TEMPERATURE POWT
C (CELSIUS), GIVEN THE FOLLOWING INPUT:
C T = TEMPERATURE (CELSIUS)
C P = PRESSURE (MILLIBARS)
C TD = DEW POINT (CELSIUS)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
DATA CTA,AKAP/273.15,0.28541/
C CTA = DIFFERENCE BETWEEN KELVIN AND CELSIUS TEMPERATURES
C AKAP = (GAS CONSTANT FOR DRY AIR) / (SPECIFIC HEAT AT
C CONSTANT PRESSURE FOR DRY AIR)
C COMPUTE THE POTENTIAL TEMPERATURE (CELSIUS)
PT = (T+CTA)*(1000./P)**AKAP-CTA
C COMPUTE THE LIFTING CONDENSATION LEVEL (LCL).
TC = TCON(T,TD)
C FOR THE ORIGIN OF THE FOLLOWING APPROXIMATION, SEE THE DOCUMEN-
C TATION FOR THE WOBUS FUNCTION.
POWT = PT-WOBF(PT)+WOBF(TC)
RETURN
END
USER_DEV:[GULIB.THERMOSRC]PRECPW_TG.FOR;1
FUNCTION PRECPW(TD,P,N)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C THIS FUNCTION COMPUTES TOTAL PRECIPITABLE WATER PRECPW (CM) IN A
C VERTICAL COLUMN OF AIR BASED UPON SOUNDING DATA AT N LEVELS:
C TD = DEW POINT (CELSIUS)
C P = PRESSURE (MILLIBARS)
C CALCULATIONS ARE DONE IN CGS UNITS.
DIMENSION TD(N),P(N)
C G = ACCELERATION DUE TO THE EARTH'S GRAVITY (CM/S**2)
DATA G/980.616/
C INITIALIZE VALUE OF PRECIPITABLE WATER
PW = 0.
NL = N-1
C CALCULATE THE MIXING RATIO AT THE LOWEST LEVEL.
WBOT = WMR(P(1),TD(1))
DO 5 I=1,NL
WTOP = WMR(P(I+1),TD(I+1))
C CALCULATE THE LAYER-MEAN MIXING RATIO (G/KG).
W = 0.5*(WTOP+WBOT)
C MAKE THE MIXING RATIO DIMENSIONLESS.
WL = .001*W
C CALCULATE THE SPECIFIC HUMIDITY.
QL = WL/(WL+1.)
C THE FACTOR OF 1000. BELOW CONVERTS FROM MILLIBARS TO DYNES/CM**2.
DP = 1000.*(P(I)-P(I+1))
PW = PW+(QL/G)*DP
WBOT = WTOP
5 CONTINUE
PRECPW = PW
RETURN
END
USER_DEV:[GULIB.THERMOSRC]PTLCL_TG.FOR;1
SUBROUTINE PTLCL(P,T,TD,PC,TC)
C THIS SUBROUTINE ESTIMATES THE PRESSURE PC (MB) AND THE TEMPERATURE
C TC (CELSIUS) AT THE LIFTED CONDENSATION LEVEL (LCL), GIVEN THE
C INITIAL PRESSURE P (MB), TEMPERATURE T (CELSIUS) AND DEW POINT
C (CELSIUS) OF THE PARCEL. THE APPROXIMATION IS THAT LINES OF
C CONSTANT POTENTIAL TEMPERATURE AND CONSTANT MIXING RATIO ARE
C STRAIGHT ON THE SKEW T/LOG P CHART.
C TETEN'S FORMULA FOR SATURATION VAPOR PRESSURE AS A FUNCTION OF
C PRESSURE WAS USED IN THE DERIVATION OF THE FORMULA BELOW. FOR
C ADDITIONAL DETAILS, SEE MATH NOTES BY T. SCHLATTER DATED 8 SEP 81.
C T. SCHLATTER, NOAA/ERL/PROFS PROGRAM OFFICE, BOULDER, COLORADO,
C WROTE THIS SUBROUTINE.
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C AKAP = (GAS CONSTANT FOR DRY AIR) / (SPECIFIC HEAT AT CONSTANT
C PRESSURE FOR DRY AIR)
C CTA = DIFFERENCE BETWEEN KELVIN AND CELSIUS TEMPERATURES
DATA AKAP,CTA/0.28541,273.15/
C1 = 4098.026/(TD+237.3)**2
C2 = 1./(AKAP*(T+CTA))
PC = P*EXP(C1*C2*(T-TD)/(C2-C1))
TC = T+C1*(T-TD)/(C2-C1)
RETURN
END
USER_DEV:[GULIB.THERMOSRC]SATLFT_TG.FOR;1
FUNCTION SATLFT(THW,P)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C INPUT: THW = WET-BULB POTENTIAL TEMPERATURE (CELSIUS).
C THW DEFINES A MOIST ADIABAT.
C P = PRESSURE (MILLIBARS)
C OUTPUT: SATLFT = TEMPERATURE (CELSIUS) WHERE THE MOIST ADIABAT
C CROSSES P
DATA CTA,AKAP/273.15,0.28541/
C CTA = DIFFERENCE BETWEEN KELVIN AND CELSIUS TEMPERATURES
C AKAP = (GAS CONSTANT FOR DRY AIR) / (SPECIFIC HEAT AT CONSTANT
C PRESSURE FOR DRY AIR)
C THE ALGORITHM BELOW CAN BEST BE UNDERSTOOD BY REFERRING TO A
C SKEW-T/LOG P CHART. IT WAS DEVISED BY HERMAN WOBUS, A MATHEMATI-
C CIAN FORMERLY AT THE NAVY WEATHER RESEARCH FACILITY BUT NOW RETIRED.
C THE VALUE RETURNED BY SATLFT CAN BE CHECKED BY REFERRING TO TABLE
C 78, PP.319-322, SMITHSONIAN METEOROLOGICAL TABLES, BY ROLAND LIST
C (6TH REVISED EDITION).
C
IF (P.NE.1000.) GO TO 5
SATLFT = THW
RETURN
5 CONTINUE
C COMPUTE TONE, THE TEMPERATURE WHERE THE DRY ADIABAT WITH VALUE THW
C (CELSIUS) CROSSES P.
PWRP = (P/1000.)**AKAP
TONE = (THW+CTA)*PWRP-CTA
C CONSIDER THE MOIST ADIABAT EW1 THROUGH TONE AT P. USING THE DEFINI-
C TION OF THE WOBUS FUNCTION (SEE DOCUMENTATION ON WOBF), IT CAN BE
C SHOWN THAT EONE = EW1-THW.
EONE = WOBF(TONE)-WOBF(THW)
RATE = 1.
GO TO 15
C IN THE LOOP BELOW, THE ESTIMATE OF SATLFT IS ITERATIVELY IMPROVED.
10 CONTINUE
C RATE IS THE RATIO OF A CHANGE IN T TO THE CORRESPONDING CHANGE IN
C E. ITS INITIAL VALUE WAS SET TO 1 ABOVE.
RATE = (TTWO-TONE)/(ETWO-EONE)
TONE = TTWO
EONE = ETWO
15 CONTINUE
C TTWO IS AN IMPROVED ESTIMATE OF SATLFT.
TTWO = TONE-EONE*RATE
C PT IS THE POTENTIAL TEMPERATURE (CELSIUS) CORRESPONDING TO TTWO AT P
PT = (TTWO+CTA)/PWRP-CTA
C CONSIDER THE MOIST ADIABAT EW2 THROUGH TTWO AT P. USING THE DEFINI-
C TION OF THE WOBUS FUNCTION, IT CAN BE SHOWN THAT ETWO = EW2-THW.
ETWO = PT+WOBF(TTWO)-WOBF(PT)-THW
C DLT IS THE CORRECTION TO BE SUBTRACTED FROM TTWO.
DLT = ETWO*RATE
IF (ABS(DLT).GT.0.1) GO TO 10
SATLFT = TTWO-DLT
RETURN
END
USER_DEV:[GULIB.THERMOSRC]SSH_TG.FOR;1
FUNCTION SSH(P,T)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C THIS FUNCTION RETURNS SATURATION SPECIFIC HUMIDITY SSH (GRAMS OF
C WATER VAPOR PER KILOGRAM OF MOIST AIR) GIVEN THE PRESSURE P
C (MILLIBARS) AND THE TEMPERATURE T (CELSIUS). THE EQUATION IS GIVEN
C IN STANDARD METEOROLOGICAL TEXTS. IF T IS DEW POINT (CELSIUS), THEN
C SSH RETURNS THE ACTUAL SPECIFIC HUMIDITY.
C COMPUTE THE DIMENSIONLESS MIXING RATIO.
W = .001*WMR(P,T)
C COMPUTE THE DIMENSIONLESS SATURATION SPECIFIC HUMIDITY.
Q = W/(1.+W)
SSH = 1000.*Q
RETURN
END
USER_DEV:[GULIB.THERMOSRC]TCON_TG.FOR;1
FUNCTION TCON(T,D)
C THIS FUNCTION RETURNS THE TEMPERATURE TCON (CELSIUS) AT THE LIFTING
C CONDENSATION LEVEL, GIVEN THE TEMPERATURE T (CELSIUS) AND THE
C DEW POINT D (CELSIUS).
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C
C COMPUTE THE DEW POINT DEPRESSION S.
S = T-D
C THE APPROXIMATION BELOW, A THIRD ORDER POLYNOMIAL IN S AND T,
C IS DUE TO HERMAN WOBUS. THE SOURCE OF DATA FOR FITTING THE
C POLYNOMIAL IS UNKNOWN.
DLT = S*(1.2185+1.278E-03*T+
1 S*(-2.19E-03+1.173E-05*S-5.2E-06*T))
TCON = T-DLT
RETURN
END
USER_DEV:[GULIB.THERMOSRC]TDA_TG.FOR;1
FUNCTION TDA(O,P)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C G.S. Stipanuk 1973 Original version.
C Reference Stipanuk paper entitled:
C "ALGORITHMS FOR GENERATING A SKEW-T, LOG P
C DIAGRAM AND COMPUTING SELECTED METEOROLOGICAL
C QUANTITIES."
C ATMOSPHERIC SCIENCES LABORATORY
C U.S. ARMY ELECTRONICS COMMAND
C WHITE SANDS MISSILE RANGE, NEW MEXICO 88002
C 33 PAGES
C Baker, Schlatter 17-MAY-1982
C THIS FUNCTION RETURNS THE TEMPERATURE TDA (CELSIUS) ON A DRY ADIABAT
C AT PRESSURE P (MILLIBARS). THE DRY ADIABAT IS GIVEN BY
C POTENTIAL TEMPERATURE O (CELSIUS). THE COMPUTATION IS BASED ON
C POISSON'S EQUATION.
OK= O+273.15
TDAK= OK*((P*.001)**.286)
TDA= TDAK-273.15
RETURN
END
USER_DEV:[GULIB.THERMOSRC]TE_TG.FOR;1
FUNCTION TE(T,TD,P)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C G.S. Stipanuk 1973 Original version.
C Reference Stipanuk paper entitled:
C "ALGORITHMS FOR GENERATING A SKEW-T, LOG P
C DIAGRAM AND COMPUTING SELECTED METEOROLOGICAL
C QUANTITIES."
C ATMOSPHERIC SCIENCES LABORATORY
C U.S. ARMY ELECTRONICS COMMAND
C WHITE SANDS MISSILE RANGE, NEW MEXICO 88002
C 33 PAGES
C Baker, Schlatter 17-MAY-1982
C THIS FUNCTION RETURNS THE EQUIVALENT TEMPERATURE TE (CELSIUS) OF A
C PARCEL OF AIR GIVEN ITS TEMPERATURE T (CELSIUS), DEW POINT (CELSIUS)
C AND PRESSURE P (MILLIBARS).
C CALCULATE EQUIVALENT POTENTIAL TEMPERATURE.
AOE = OE(T,TD,P)
C USE POISSONS'S EQUATION TO CALCULATE EQUIVALENT TEMPERATURE.
TE= TDA(AOE,P)
RETURN
END
USER_DEV:[GULIB.THERMOSRC]THM_TG.FOR;1
FUNCTION THM(T,P)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C THIS FUNCTION RETURNS THE WET-BULB POTENTIAL TEMPERATURE THM
C (CELSIUS) CORRESPONDING TO A PARCEL OF AIR SATURATED AT
C TEMPERATURE T (CELSIUS) AND PRESSURE P (MILLIBARS).
F(X) = 1.8199427E+01+X*( 2.1640800E-01+X*( 3.0716310E-04+X*
1 (-3.8953660E-06+X*( 1.9618200E-08+X*( 5.2935570E-11+X*
2 ( 7.3995950E-14+X*(-4.1983500E-17)))))))
THM = T
IF (P.EQ.1000.) RETURN
C COMPUTE THE POTENTIAL TEMPERATURE (CELSIUS).
THD = (T+273.15)*(1000./P)**.286-273.15
THM = THD+6.071*(EXP(T/F(T))-EXP(THD/F(THD)))
RETURN
END
USER_DEV:[GULIB.THERMOSRC]TLCL1_TG.FOR;1
FUNCTION TLCL1(T,TD)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C THIS FUNCTION RETURNS THE TEMPERATURE TLCL1 (CELSIUS) OF THE LIFTING
C CONDENSATION LEVEL (LCL) GIVEN THE INITIAL TEMPERATURE T (CELSIUS)
C AND DEW POINT TD (CELSIUS) OF A PARCEL OF AIR.
C ERIC SMITH AT COLORADO STATE UNIVERSITY HAS USED THE FORMULA
C BELOW, BUT ITS ORIGIN IS UNKNOWN.
DATA CTA/273.15/
C CTA = DIFFERENCE BETWEEN KELVIN AND CELSIUS TEMPERATURE
TK = T+CTA
C COMPUTE THE PARCEL VAPOR PRESSURE (MB).
ES = ESLO(TD)
TLCL = 2840./(3.5*ALOG(TK)-ALOG(ES)-4.805)+55.
TLCL1 = TLCL-CTA
RETURN
END
USER_DEV:[GULIB.THERMOSRC]TLCL_TG.FOR;1
FUNCTION TLCL(T,TD)
C THIS FUNCTION YIELDS THE TEMPERATURE TLCL (CELSIUS) OF THE LIFTING
C CONDENSATION LEVEL, GIVEN THE TEMPERATURE T (CELSIUS) AND THE
C DEW POINT TD (CELSIUS). THE FORMULA USED IS IN BOLTON, DAVID,
C 1980: "THE COMPUTATION OF EQUIVALENT POTENTIAL TEMPERATURE,"
C MONTHLY WEATHER REVIEW, VOL. 108, NO. 7 (JULY), P. 1048, EQ.(15).
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C CONVERT FROM CELSIUS TO KELVIN DEGREES.
TK = T+273.15
TDK = TD+273.15
A = 1./(TDK-56.)
B = ALOG(TK/TDK)/800.
TC = 1./(A+B)+56.
TLCL = TC-273.15
RETURN
END
USER_DEV:[GULIB.THERMOSRC]TMLAPS_TG.FOR;1
FUNCTION TMLAPS(THETAE,P)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C THIS FUNCTION RETURNS THE TEMPERATURE TMLAPS (CELSIUS) AT PRESSURE
C P (MILLIBARS) ALONG THE MOIST ADIABAT CORRESPONDING TO AN EQUIVALENT
C POTENTIAL TEMPERATURE THETAE (CELSIUS).
C THE ALGORITHM WAS WRITTEN BY ERIC SMITH AT COLORADO STATE
C UNIVERSITY.
DATA CRIT/0.1/
C CTA = DIFFERENCE BETWEEN KELVIN AND CELSIUS TEMPERATURES.
C CRIT = CONVERGENCE CRITERION (DEGREES KELVIN)
EQ0 = THETAE
C INITIAL GUESS FOR SOLUTION
TLEV = 25.
C COMPUTE THE SATURATION EQUIVALENT POTENTIAL TEMPERATURE CORRESPON-
C DING TO TEMPERATURE TLEV AND PRESSURE P.
EQ1 = EPT(TLEV,TLEV,P)
DIF = ABS(EQ1-EQ0)
IF (DIF.LT.CRIT) GO TO 3
IF (EQ1.GT.EQ0) GO TO 1
C DT IS THE INITIAL STEPPING INCREMENT.
DT = 10.
I = -1
GO TO 2
1 DT = -10.
I = 1
2 TLEV = TLEV+DT
EQ1 = EPT(TLEV,TLEV,P)
DIF = ABS(EQ1-EQ0)
IF (DIF.LT.CRIT) GO TO 3
J = -1
IF (EQ1.GT.EQ0) J=1
IF (I.EQ.J) GO TO 2
C THE SOLUTION HAS BEEN PASSED. REVERSE THE DIRECTION OF SEARCH
C AND DECREASE THE STEPPING INCREMENT.
TLEV = TLEV-DT
DT = DT/10.
GO TO 2
3 TMLAPS = TLEV
RETURN
END
USER_DEV:[GULIB.THERMOSRC]TMR_TG.FOR;1
FUNCTION TMR(W,P)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C G.S. Stipanuk 1973 Original version.
C Reference Stipanuk paper entitled:
C "ALGORITHMS FOR GENERATING A SKEW-T, LOG P
C DIAGRAM AND COMPUTING SELECTED METEOROLOGICAL
C QUANTITIES."
C ATMOSPHERIC SCIENCES LABORATORY
C U.S. ARMY ELECTRONICS COMMAND
C WHITE SANDS MISSILE RANGE, NEW MEXICO 88002
C 33 PAGES
C Baker, Schlatter 17-MAY-1982
C THIS FUNCTION RETURNS THE TEMPERATURE (CELSIUS) ON A MIXING
C RATIO LINE W (G/KG) AT PRESSURE P (MB). THE FORMULA IS GIVEN IN
C TABLE 1 ON PAGE 7 OF STIPANUK (1973).
C
C INITIALIZE CONSTANTS
DATA C1/.0498646455/,C2/2.4082965/,C3/7.07475/
DATA C4/38.9114/,C5/.0915/,C6/1.2035/
X= ALOG10(W*P/(622.+W))
TMRK= 10.**(C1*X+C2)-C3+C4*((10.**(C5*X)-C6)**2.)
TMR= TMRK-273.15
RETURN
END
USER_DEV:[GULIB.THERMOSRC]TSA_TG.FOR;1
FUNCTION TSA(OS,P)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C G.S. Stipanuk 1973 Original version.
C Reference Stipanuk paper entitled:
C "ALGORITHMS FOR GENERATING A SKEW-T, LOG P
C DIAGRAM AND COMPUTING SELECTED METEOROLOGICAL
C QUANTITIES."
C ATMOSPHERIC SCIENCES LABORATORY
C U.S. ARMY ELECTRONICS COMMAND
C WHITE SANDS MISSILE RANGE, NEW MEXICO 88002
C 33 PAGES
C Baker, Schlatter 17-MAY-1982
C THIS FUNCTION RETURNS THE TEMPERATURE TSA (CELSIUS) ON A SATURATION
C ADIABAT AT PRESSURE P (MILLIBARS). OS IS THE EQUIVALENT POTENTIAL
C TEMPERATURE OF THE PARCEL (CELSIUS). SIGN(A,B) REPLACES THE
C ALGEBRAIC SIGN OF A WITH THAT OF B.
C B IS AN EMPIRICAL CONSTANT APPROXIMATELY EQUAL TO 0.001 OF THE LATENT
C HEAT OF VAPORIZATION FOR WATER DIVIDED BY THE SPECIFIC HEAT AT CONSTANT
C PRESSURE FOR DRY AIR.
DATA B/2.6518986/
A= OS+273.15
C TQ IS THE FIRST GUESS FOR TSA.
TQ= 253.15
C D IS AN INITIAL VALUE USED IN THE ITERATION BELOW.
D= 120.
C ITERATE TO OBTAIN SUFFICIENT ACCURACY....SEE TABLE 1, P.8
C OF STIPANUK (1973) FOR EQUATION USED IN ITERATION.
DO 1 I= 1,12
TQK= TQ-273.15
D= D/2.
X= A*EXP(-B*W(TQK,P)/TQ)-TQ*((1000./P)**.286)
IF (ABS(X).LT.1E-7) GO TO 2
TQ= TQ+SIGN(D,X)
1 CONTINUE
2 TSA= TQ-273.15
RETURN
END
USER_DEV:[GULIB.THERMOSRC]TV_TG.FOR;1
FUNCTION TV(T,TD,P)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C THIS FUNCTION RETURNS THE VIRTUAL TEMPERATURE TV (CELSIUS) OF
C A PARCEL OF AIR AT TEMPERATURE T (CELSIUS), DEW POINT TD
C (CELSIUS), AND PRESSURE P (MILLIBARS). THE EQUATION APPEARS
C IN MOST STANDARD METEOROLOGICAL TEXTS.
DATA CTA,EPS/273.15,0.62197/
C CTA = DIFFERENCE BETWEEN KELVIN AND CELSIUS TEMPERATURES.
C EPS = RATIO OF THE MEAN MOLECULAR WEIGHT OF WATER (18.016 G/MOLE)
C TO THAT OF DRY AIR (28.966 G/MOLE)
TK = T+CTA
C CALCULATE THE DIMENSIONLESS MIXING RATIO.
W = .001*WMR(P,TD)
TV = TK*(1.+W/EPS)/(1.+W)-CTA
RETURN
END
USER_DEV:[GULIB.THERMOSRC]TW_TG.FOR;1
FUNCTION TW(T,TD,P)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C G.S. Stipanuk 1973 Original version.
C Reference Stipanuk paper entitled:
C "ALGORITHMS FOR GENERATING A SKEW-T, LOG P
C DIAGRAM AND COMPUTING SELECTED METEOROLOGICAL
C QUANTITIES."
C ATMOSPHERIC SCIENCES LABORATORY
C U.S. ARMY ELECTRONICS COMMAND
C WHITE SANDS MISSILE RANGE, NEW MEXICO 88002
C 33 PAGES
C Baker, Schlatter 17-MAY-1982
C THIS FUNCTION RETURNS THE WET-BULB TEMPERATURE TW (CELSIUS)
C GIVEN THE TEMPERATURE T (CELSIUS), DEW POINT TD (CELSIUS)
C AND PRESSURE P (MB). SEE P.13 IN STIPANUK (1973), REFERENCED
C ABOVE, FOR A DESCRIPTION OF THE TECHNIQUE.
C
C
C DETERMINE THE MIXING RATIO LINE THRU TD AND P.
AW = W(TD,P)
C
C DETERMINE THE DRY ADIABAT THRU T AND P.
AO = O(T,P)
PI = P
C ITERATE TO LOCATE PRESSURE PI AT THE INTERSECTION OF THE TWO
C CURVES . PI HAS BEEN SET TO P FOR THE INITIAL GUESS.
DO 4 I= 1,10
X= .02*(TMR(AW,PI)-TDA(AO,PI))
IF (ABS(X).LT.0.01) GO TO 5
4 PI= PI*(2.**(X))
C FIND THE TEMPERATURE ON THE DRY ADIABAT AO AT PRESSURE PI.
5 TI= TDA(AO,PI)
C THE INTERSECTION HAS BEEN LOCATED...NOW, FIND A SATURATION
C ADIABAT THRU THIS POINT. FUNCTION OS RETURNS THE EQUIVALENT
C POTENTIAL TEMPERATURE (C) OF A PARCEL SATURATED AT TEMPERATURE
C TI AND PRESSURE PI.
AOS= OS(TI,PI)
C FUNCTION TSA RETURNS THE WET-BULB TEMPERATURE (C) OF A PARCEL AT
C PRESSURE P WHOSE EQUIVALENT POTENTIAL TEMPERATURE IS AOS.
TW = TSA(AOS,P)
RETURN
END
USER_DEV:[GULIB.THERMOSRC]WMR_TG.FOR;1
FUNCTION WMR(P,T)
C THIS FUNCTION APPROXIMATES THE MIXING RATIO WMR (GRAMS OF WATER
C VAPOR PER KILOGRAM OF DRY AIR) GIVEN THE PRESSURE P (MB) AND THE
C TEMPERATURE T (CELSIUS). THE FORMULA USED IS GIVEN ON P. 302 OF THE
C SMITHSONIAN METEOROLOGICAL TABLES BY ROLAND LIST (6TH EDITION).
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C EPS = RATIO OF THE MEAN MOLECULAR WEIGHT OF WATER (18.016 G/MOLE)
C TO THAT OF DRY AIR (28.966 G/MOLE)
DATA EPS/0.62197/
C THE NEXT TWO LINES CONTAIN A FORMULA BY HERMAN WOBUS FOR THE
C CORRECTION FACTOR WFW FOR THE DEPARTURE OF THE MIXTURE OF AIR
C AND WATER VAPOR FROM THE IDEAL GAS LAW. THE FORMULA FITS VALUES
C IN TABLE 89, P. 340 OF THE SMITHSONIAN METEOROLOGICAL TABLES,
C BUT ONLY FOR TEMPERATURES AND PRESSURES NORMALLY ENCOUNTERED IN
C IN THE ATMOSPHERE.
X = 0.02*(T-12.5+7500./P)
WFW = 1.+4.5E-06*P+1.4E-03*X*X
FWESW = WFW*ESW(T)
R = EPS*FWESW/(P-FWESW)
C CONVERT R FROM A DIMENSIONLESS RATIO TO GRAMS/KILOGRAM.
WMR = 1000.*R
RETURN
END
USER_DEV:[GULIB.THERMOSRC]WOBF_TG.FOR;1
FUNCTION WOBF(T)
C THIS FUNCTION CALCULATES THE DIFFERENCE OF THE WET-BULB POTENTIAL
C TEMPERATURES FOR SATURATED AND DRY AIR GIVEN THE TEMPERATURE.
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C Baker, Schlatter 17-MAY-1982 Original version.
C LET WBPTS = WET-BULB POTENTIAL TEMPERATURE FOR SATURATED
C AIR AT TEMPERATURE T (CELSIUS). LET WBPTD = WET-BULB POTENTIAL
C TEMPERATURE FOR COMPLETELY DRY AIR AT THE SAME TEMPERATURE T.
C THE WOBUS FUNCTION WOBF (IN DEGREES CELSIUS) IS DEFINED BY
C WOBF(T) = WBPTS-WBPTD.
C ALTHOUGH WBPTS AND WBPTD ARE FUNCTIONS OF BOTH PRESSURE AND
C TEMPERATURE, THEIR DIFFERENCE IS A FUNCTION OF TEMPERATURE ONLY.
C TO UNDERSTAND WHY, CONSIDER A PARCEL OF DRY AIR AT TEMPERA-
C TURE T AND PRESSURE P. THE THERMODYNAMIC STATE OF THE PARCEL IS
C REPRESENTED BY A POINT ON A PSEUDOADIABATIC CHART. THE WET-BULB
C POTENTIAL TEMPERATURE CURVE (MOIST ADIABAT) PASSING THROUGH THIS
C POINT IS WBPTS. NOW T IS THE EQUIVALENT TEMPERATURE FOR ANOTHER
C PARCEL SATURATED AT SOME LOWER TEMPERATURE TW, BUT AT THE SAME
C PRESSURE P. TO FIND TW, ASCEND ALONG THE DRY ADIABAT THROUGH
C (T,P). AT A GREAT HEIGHT, THE DRY ADIABAT AND SOME MOIST
C ADIABAT WILL NEARLY COINCIDE. DESCEND ALONG THIS MOIST ADIABAT
C BACK TO P. THE PARCEL TEMPERATURE IS NOW TW. THE WET-BULB
C POTENTIAL TEMPERATURE CURVE (MOIST ADIABAT) THROUGH (TW,P) IS WBPTD.
C THE DIFFERENCE (WBPTS-WBPTD) IS PROPORTIONAL TO THE HEAT IMPARTED
C TO A PARCEL SATURATED AT TEMPERATURE TW IF ALL ITS WATER VAPOR
C WERE CONDENSED. SINCE THE AMOUNT OF WATER VAPOR A PARCEL CAN
C HOLD DEPENDS UPON TEMPERATURE ALONE, (WBPTD-WBPTS) MUST DEPEND
C ON TEMPERATURE ALONE.
C THE WOBUS FUNCTION IS USEFUL FOR EVALUATING SEVERAL THERMO-
C DYNAMIC QUANTITIES. BY DEFINITION:
C WOBF(T) = WBPTS-WBPTD. (1)
C IF T IS AT 1000 MB, THEN T IS A POTENTIAL TEMPERATURE PT AND
C WBPTS = PT. THUS
C WOBF(PT) = PT-WBPTD. (2)
C IF T IS AT THE CONDENSATION LEVEL, THEN T IS THE CONDENSATION
C TEMPERATURE TC AND WBPTS IS THE WET-BULB POTENTIAL TEMPERATURE
C WBPT. THUS
C WOBF(TC) = WBPT-WBPTD. (3)
C IF WBPTD IS ELIMINATED FROM (2) AND (3), THERE RESULTS
C WBPT = PT-WOBF(PT)+WOBF(TC).
C IF WBPTD IS ELIMINATED FROM (1) AND (2), THERE RESULTS
C WBPTS = PT-WOBF(PT)+WOBF(T).
C IF T IS AN EQUIVALENT POTENTIAL TEMPERATURE EPT (IMPLYING
C THAT THE AIR AT 1000 MB IS COMPLETELY DRY), THEN WBPTS = EPT
C AND WBPTD = WBPT. THUS
C WOBF(EPT) = EPT-WBPT.
C THIS FORM IS THE BASIS FOR A POLYNOMIAL APPROXIMATION TO WOBF.
C IN TABLE 78 ON PP.319-322 OF THE SMITHSONIAN METEOROLOGICAL
C TABLES BY ROLAND LIST (6TH REVISED EDITION), ONE FINDS WET-BULB
C POTENTIAL TEMPERATURES AND THE CORRESPONDING EQUIVALENT POTENTIAL
C TEMPERATURES LISTED TOGETHER. HERMAN WOBUS, A MATHEMATICIAN FOR-
C MERLY AT THE NAVY WEATHER RESEARCH FACILITY, NORFOLK, VIRGINIA,
C AND NOW RETIRED, COMPUTED THE COEFFICIENTS FOR THE POLYNOMIAL
C APPROXIMATION FROM NUMBERS IN THIS TABLE.
C
C NOTES BY T.W. SCHLATTER
C NOAA/ERL/PROFS PROGRAM OFFICE
C AUGUST 1981
X = T-20.
IF (X.GT.0.) GO TO 10
POL = 1. +X*(-8.8416605E-03
1 +X*( 1.4714143E-04 +X*(-9.6719890E-07
2 +X*(-3.2607217E-08 +X*(-3.8598073E-10)))))
WOBF = 15.130/POL**4
RETURN
10 CONTINUE
POL = 1. +X*( 3.6182989E-03
1 +X*(-1.3603273E-05 +X*( 4.9618922E-07
2 +X*(-6.1059365E-09 +X*( 3.9401551E-11
3 +X*(-1.2588129E-13 +X*( 1.6688280E-16)))))))
WOBF = 29.930/POL**4+0.96*X-14.8
RETURN
END
USER_DEV:[GULIB.THERMOSRC]W_TG.FOR;1
FUNCTION W(T,P)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C G.S. Stipanuk 1973 Original version.
C Reference Stipanuk paper entitled:
C "ALGORITHMS FOR GENERATING A SKEW-T, LOG P
C DIAGRAM AND COMPUTING SELECTED METEOROLOGICAL
C QUANTITIES."
C ATMOSPHERIC SCIENCES LABORATORY
C U.S. ARMY ELECTRONICS COMMAND
C WHITE SANDS MISSILE RANGE, NEW MEXICO 88002
C 33 PAGES
C Baker, Schlatter 17-MAY-1982
C THIS FUNCTION RETURNS THE MIXING RATIO (GRAMS OF WATER VAPOR PER
C KILOGRAM OF DRY AIR) GIVEN THE DEW POINT (CELSIUS) AND PRESSURE
C (MILLIBARS). IF THE TEMPERTURE IS INPUT INSTEAD OF THE
C DEW POINT, THEN SATURATION MIXING RATIO (SAME UNITS) IS RETURNED.
C THE FORMULA IS FOUND IN MOST METEOROLOGICAL TEXTS.
X= ESAT(T)
W= 622.*X/(P-X)
RETURN
END
USER_DEV:[GULIB.THERMOSRC]Z_TG.FOR;1
FUNCTION Z(PT,P,T,TD,N)
C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST'
C G.S. Stipanuk 1973 Original version.
C Reference Stipanuk paper entitled:
C "ALGORITHMS FOR GENERATING A SKEW-T, LOG P
C DIAGRAM AND COMPUTING SELECTED METEOROLOGICAL
C QUANTITIES."
C ATMOSPHERIC SCIENCES LABORATORY
C U.S. ARMY ELECTRONICS COMMAND
C WHITE SANDS MISSILE RANGE, NEW MEXICO 88002
C 33 PAGES
C Baker, Schlatter 17-MAY-1982
C THIS FUNCTION RETURNS THE THICKNESS OF A LAYER BOUNDED BY PRESSURE
C P(1) AT THE BOTTOM AND PRESSURE PT AT THE TOP.
C ON INPUT:
C P = PRESSURE (MB). NOTE THAT P(I).GT.P(I+1).
C T = TEMPERATURE (CELSIUS)
C TD = DEW POINT (CELSIUS)
C N = NUMBER OF LEVELS IN THE SOUNDING AND THE DIMENSION OF
C P, T AND TD
C ON OUTPUT:
C Z = GEOMETRIC THICKNESS OF THE LAYER (M)
C THE ALGORITHM INVOLVES NUMERICAL INTEGRATION OF THE HYDROSTATIC
C EQUATION FROM P(1) TO PT. IT IS DESCRIBED ON P.15 OF STIPANUK
C (1973).
DIMENSION T(1),P(1),TD(1),TK(100)
C C1 = .001*(1./EPS-1.) WHERE EPS = .62197 IS THE RATIO OF THE
C MOLECULAR WEIGHT OF WATER TO THAT OF
C DRY AIR. THE FACTOR 1000. CONVERTS THE
C MIXING RATIO W FROM G/KG TO A DIMENSION-
C LESS RATIO.
C C2 = R/(2.*G) WHERE R IS THE GAS CONSTANT FOR DRY AIR
C (287 KG/JOULE/DEG K) AND G IS THE ACCELERATION
C DUE TO THE EARTH'S GRAVITY (9.8 M/S**2). THE
C FACTOR OF 2 IS USED IN AVERAGING TWO VIRTUAL
C TEMPERATURES.
DATA C1/.0006078/,C2/14.64285/
DO 5 I= 1,N
TK(I)= T(I)+273.15
5 CONTINUE
Z= 0.0
IF (PT.LT.P(N)) GO TO 20
I= 0
10 I= I+1
J= I+1
IF (PT.GE.P(J)) GO TO 15
A1= TK(J)*(1.+C1*W(TD(J),P(J)))
A2= TK(I)*(1.+C1*W(TD(I),P(I)))
Z= Z+C2*(A1+A2)*(ALOG(P(I)/P(J)))
GO TO 10
15 CONTINUE
A1= TK(J)*(1.+C1*W(TD(J),P(J)))
A2= TK(I)*(1.+C1*W(TD(I),P(I)))
Z= Z+C2*(A1+A2)*(ALOG(P(I)/PT))
RETURN
20 Z= -1.0
RETURN
END