	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
	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
