;;################################################################################# ;; FUNCTION WVSatPressureLiquidWater, TK, punit ;+ ;NAME: ; WVSatPressureLiquidWater ; ;PURPOSE: ; The recommended formulas for the calculation of equilibrium water vapor pressure ; in the pure phase (no other gases present) over a plane surface of liquid water ; or ice are those of Sonntag, 1994. (References below.) The formulas are quoted ; also in the papers of Leiterer et al., 1997 and Helten et al., 1999. I have ; verified that they are exactly the same in all three papers. (Sonntag gives ; the formulas for partial pressure in Pa and hPa, Leiterer et al. for hPa, ; Helten et al. for Pa.) ; ; The equation and coefficients below are for the equilibrium water vapor pressure ; in Pa/hPa and the temperature in K. ; ; The temperature of 0C corresponds to 273.15K. (Not 273.16K, as stated in the ; Leiterer paper.) ; ; References: ; Sonntag, D., Advancements in the field of hygrometry, ; Meteorologische Zeitschrift, 3, 51-66, 1994. ; ; Helten, M. et al, In-flight comparison of MOZAIC and POLINAT water ; vapor measurements, JGR, 104, 26.087-26.096, 1999. ; ; Leiterer U. et al, Improvements in Radiosonde Humidity Profiles ; Using RS80/RS90 Radiosondes of Vaisala, ; Beitr. Phys. Atmosph., 70(4), 319-336, 1997. ; ; --------------------------------------------------------------------- ; INPUT : TK DOUBLE temperature [K] ; punit STRING pressure unit 'Pa' or 'hPa' ; ; OUTPUT : equilibrium water vapor pressure in the pure phase ; (no other gases present) over a plane surface of ; liquid water [Pa] ; ; HISTORY: alpha version, 2003-04-02, TKS ; --------------------------------------------------------------------- ;- ;###################################################################### a = DOUBLE(-6096.9385) if (STRUPCASE(punit) EQ 'PA') then begin b = DOUBLE(21.2409642) endif else begin b = DOUBLE(16.635794) endelse c = DOUBLE(-2.711193e-2) d = DOUBLE(1.673952e-5) e = DOUBLE(2.433502) ew = EXP( DOUBLE(a/TK) + $ b + $ DOUBLE(c*TK) + $ DOUBLE(d*TK*TK) + $ DOUBLE(e*ALOG(TK)) ) return, DOUBLE(eW) ;; saturation pressure in Pa or hPa end ;; ;; ################################################################################# ;; FUNCTION WVSatPressureIce, TK, punit ;+ ;NAME: ; WVSatPressureIce ; ;PURPOSE: ; The recommended formulas for the calculation of equilibrium water vapor pressure ; in the pure phase (no other gases present) over a plane surface of liquid water ; or ice are those of Sonntag, 1994. (References below.) The formulas are quoted ; also in the papers of Leiterer et al., 1997 and Helten et al., 1999. I have ; verified that they are exactly the same in all three papers. (Sonntag gives ; the formulas for partial pressure in Pa and hPa, Leiterer et al. for hPa, ; Helten et al. for Pa.) ; ; The equation and coefficients below are for the equilibrium water vapor pressure ; in Pa/hPa and the temperature in K. ; ; The temperature of 0C corresponds to 273.15K. (Not 273.16K, as stated in the ; Leiterer paper.) ; ; References: ; Sonntag, D., Advancements in the field of hygrometry, ; Meteorologische Zeitschrift, 3, 51-66, 1994. ; ; Helten, M. et al, In-flight comparison of MOZAIC and POLINAT water ; vapor measurements, JGR, 104, 26.087-26.096, 1999. ; ; Leiterer U. et al, Improvements in Radiosonde Humidity Profiles ; Using RS80/RS90 Radiosondes of Vaisala, ; Beitr. Phys. Atmosph., 70(4), 319-336, 1997. ; --------------------------------------------------------------------- ; INPUT : TK DOUBLE temperature [K] ; punit STRING pressure unit 'Pa' or 'hPa' ; ; OUTPUT : equilibrium water vapor pressure in the pure phase ; (no other gases present) over ice [Pa] ; ; HISTORY: alpha version, 2003-04-02, TKS ; --------------------------------------------------------------------- ;- ;###################################################################### a = DOUBLE(-6024.5282) if (STRUPCASE(punit) EQ 'PA') then begin b = DOUBLE(29.32707) endif else begin b = DOUBLE(24.7219) endelse c = DOUBLE(1.0613868e-2) d = DOUBLE(-1.3198825e-5) e = DOUBLE(-0.49382577) ei = EXP( DOUBLE(a/TK) + $ b + $ DOUBLE(c*TK) + $ DOUBLE(d*TK*TK) + $ DOUBLE(e*ALOG(TK)) ) return, DOUBLE(ei) ;; saturation pressure in Pa or hPa end ;; ;; ################################################################################# ;; FUNCTION AirCorrFunWater, P, TK ; ; ---------------------------------------------------------------- ; INPUT ; P DOUBLE total pressure [Pa] ; TK DOUBLE temperature [K] ; ; OUTPUT ; fw DOUBLE correction function for moist air over water ; ---------------------------------------------------------------- P = (P * 1.00000E-2) ;; pressure in hPa TC = DOUBLE(TK - 273.15) ;; temperature in Celsius ew = WVSatPressureLiquidWater( DOUBLE(TK), 'hPa' ) a = ew * 1.000E-4 / (273.000E0+TC) b = 38.000E0 + $ 173.000E0*EXP(-TC/43.000E0)*(1.000E0-(ew/p)) + $ (6.390E0+4.280E0*EXP(-TC/107.000E0))*((p/ew)-1.000E0) fw = DOUBLE(1.000E0 + (a * b)) return, fw END ;; ;;################################################################################# ;; FUNCTION AirCorrFunIce, P, TK ; ; -------------------------------------------------------------- ; INPUT ; P DOUBLE total pressure [Pa] ; TK DOUBLE temperature [K] ; ; OUTPUT ; fi DOUBLE correction function for moist air over ice ; -------------------------------------------------------------- P = (P * 1.00000E-2) ;; pressure in hPa TC = DOUBLE(TK - 273.15) ;; temperature in Celsius ei = WVSatPressureLiquidIce( DOUBLE(TK), 'hPa' ) a = ei * 1.000E-5 / (273.000E0+TC) b = (2100.000E0-65.000E0*TC)*(1.000E0-(ei/p)) + $ (109.000E0-0.350E0*TC+(TC*TC/338.000E0))*((p/ei)-1.000E0) fi = DOUBLE(1.000E0 + (a * b)) return, fi END ;; ;; ################################################################################# ;; FUNCTION WaterVaporSatPressure, TK, PTOT=PTOT, $ phase=phase, punit=punit, corr=corr ; ; --------------------------------------------------------------------- ;+ ;NAME: ; WaterVaporSatPressure ; ; INPUT: ; ESSENTIAL INPUT PARAMETERS: ; TK DBLE temperature [T] ; OPTIONAL INPUT PARAMETERS (KEYWORDS): ; PTOT DBLE total pressure [Pa] ; PTOT has to be given if the correction for ; moist air should be done, i.e. corr='yes' ; phase STRING water phase, possible phases are 'water' and ; 'ice', default phase is liquid water. ; punit STRING pressure unit of the saturation pressure ; possible units are 'Pa' and 'hPa' and 'mbar' ; corr STRING consider/not consider the correction term for ; moist air instead of pure water vapor over ; waterr/ice possible values are 'yes' and 'no'. ; The default is 'no', i.e. no correction calculation. ; ; OUTPUT: ; es DBLE water vapor saturtation pressure [Pa, hPa, or mbar]. ; If no unit is specified by the input ; parameter 'punit' es is given in units of [Pa] ; If the calculation had some error detected, ; the return value is es = -9.99999e9 ; ; COMMENT: ; formula taken from: ; Sonntag, D., Advancements in the field of hygrometry, ; Meteorologische Zeitschrift, 3, 51-66, 1994. ; ; HISTORY: ; 2003-04-02, TKS: alpha version, ; 2003-04-04, TKS: added PTOT as optional if corr='no' no ; PTOT is needed ;- ; --------------------------------------------------------------------- ;; default value (return value for erroneous calculation) es = -9.99999e9 ;; check the input variable 'phase' ;; default phase is liquid water if not KEYWORD_SET(phase) then phase='WATER' ;; check the input variable 'pressure unit' ;; default unit is Pascal if not KEYWORD_SET(punit) then punit='PA' if (STRUPCASE(punit) EQ 'MBAR') then punit='HPA' ;; check the input variable for 'correction for moist air' ;; default is no correction will be applied for moist air instead of ;; pure water vapor if not KEYWORD_SET(corr) then corr='no' if ( (STRUPCASE(phase) NE 'WATER') AND (STRUPCASE(phase) NE 'ICE') ) then begin print,'WaterVaporSatPressure> !!! ERROR!' print,'WaterVaporSatPressure> !!! wrong water phase unit selected,' print,'WaterVaporSatPressure> !!! possible units are "water" and "ice"' return, es endif if ( (STRUPCASE(punit) NE 'PA') AND (STRUPCASE(punit) NE 'HPA') ) then begin print,'WaterVaporSatPressure> !!! ERROR!' print,'WaterVaporSatPressure> !!! wrong pressure unit selected,' print,'WaterVaporSatPressure> !!! possible units are "Pa" and "hPa"' return, es endif ;; check the input variable temperature if ( STRUPCASE(phase) EQ 'WATER' ) then begin if ((TK LT 173.15) OR (TK GT 373.15)) then begin print,'WaterVaporSatPressure> !!! ERROR!' print,'WaterVaporSatPressure> !!! temperature is out of the validity range' print,'WaterVaporSatPressure> !!! for the saturation over water: T=',TK,'K' print,'WaterVaporSatPressure> !!! valid range: 173.15K < T < 373.15K' return, es endif endif if ( STRUPCASE(phase) EQ 'ICE' ) then begin if ((TK LT 173.15) OR (TK GT 273.16)) then begin print,'WaterVaporSatPressure> !!! ERROR!' print,'WaterVaporSatPressure> !!! temperature is out of the validity range' print,'WaterVaporSatPressure> !!! for the saturation over ice: T=',TK,'K' print,'WaterVaporSatPressure> !!! valid range: 173.15K < T < 273.16K' return, es endif endif ;; calculation for pure water vapor over liquid water [Pa or hPa] if ( (STRUPCASE(phase) EQ 'WATER') AND (STRUPCASE(corr) EQ 'NO') ) then begin es = WVSatPressureLiquidWater( DOUBLE(TK), punit ) endif if ( (STRUPCASE(phase) EQ 'WATER') AND (STRUPCASE(corr) EQ 'YES') ) then begin ;; apply the correction term for moist air if necessary if ( KEYWORD_SET(PTOT) ) then begin ;; if the total pressure is provided by the user then perform ;; the correction calculation es = WVSatPressureLiquidWater( DOUBLE(TK), punit ) * $ AirCorrFunWater(PTOT, TK) endif else begin print,'WaterVaporSatPressure> !!! ERROR!' print,'WaterVaporSatPressure> !!! moist air correction should be applied' print,'WaterVaporSatPressure> !!! but no total pressure is given by' print,'WaterVaporSatPressure> !!! the input variable PTOT' endelse return, es;; [Pa or hPa or mbar] endif ;; calculation for pure water vapor over ice [Pa or hPa] if ( (STRUPCASE(phase) EQ 'ICE') AND (STRUPCASE(corr) EQ 'NO') ) then begin es = WVSatPressureIce( DOUBLE(TK), punit) endif if ( (STRUPCASE(phase) EQ 'ICE') AND (STRUPCASE(corr) EQ 'YES') ) then begin ;; apply the correction term for moist air if necessary if ( KEYWORD_SET(PTOT) ) then begin ;; if the total pressure is provided by the user then perform ;; the correction calculation es = WVSatPressureIce( DOUBLE(TK), punit) * $ AirCorrFunIce(DOUBLE(PTOT), DOUBLE(TK)) endif else begin ;; if the total pressure is NOT provided by the user then ;; print an error message and return the default es value print,'WaterVaporSatPressure> !!! ERROR!' print,'WaterVaporSatPressure> !!! moist air correction should be applied' print,'WaterVaporSatPressure> !!! but no total pressure is given by' print,'WaterVaporSatPressure> !!! the input variable PTOT' endelse return, es;; [Pa or hPa or mbar] endif return, es; water vapor saturation pressure [Pa or hPa or mbar] end ;;#################################################################################