; FUNCTION gasabsr02, f, TK, rhowv, Pa, absair, abswv, bad=bad, volume=volume ;----------------------------------------------------------------- ; ROSENKRANZ (1998) model -- 2002 update ; ; Reference: ; "Water vapor microwave continuum absorption: a comparison of ; measurements and results" Radio Science Vol.33, PP.919-928 ; (1998); V.34, P.1025 (1999). ; ; Top level interface (GasabsRxx) by G. Petty ; ; This subroutine calculates the mass extinction coefficient of the ; dry and vapor components of air. ; ;---------------------------------------------------------------- ; This file contains subroutines for computing atmospheric absorption at ; microwave wavelengths, as supplied by P. Rosenkranz (2/98) via his ; anonymous ftp server (mesa.mit.edu; login anonymous; go to phil/lbl_rt). ; Consolidated into one file by G. Petty. ; Converted to Fortran 90 (9/20/02) by M. Walters. ; Converted to IDL (9/03) by C. O'Dell ;---------------------------------------------------------------- ; ; ; INPUT VARIABLES ; F = frequency (GHz), ; Tk = absolute temperature (K) ; Rhowv = water vapor density (kg/m**3). ; Pa = Total air pressure (Pascals). ; *note. F must be a scalar. Tk, Rhowv, and Pa must ; be arrays of the same size and shape. ; OUTPUT VARS ; ; Absair = extinction by dry air (meters squared per kg *moist air*) ; Abswv = extinction by water vapor (meters squared per kg water vapor) ; *note: output arrays will be the same size and shape as Tk, Rhowv, Pa. ; real :: Pmb, vapden ; real :: e, q, Tv, rhoair ; real :: abs_n2, abs_o2, abs_h2o ; KEYWORDS ; bad : set this to an unnamed variable will return the indices ; where the input data are out of bounds, and hence ; the calculations are meaningless. ; ; volume: set this keyword to report the output in m^-1 (volume extinction) ; as opposed to m^2 kg^-1 (mass extinction) ;-------------------------------------------------------------------------- ;---------------------------------------------------------------- ; Subroutines from P. W. Rosenkranz ; Copyright (c) 2002 Massachusetts Institute of Technology ;---------------------------------------------------------------- function abs_n2, T, P, f ; abs_n2 = absorption coefficient due to Nitrogen in air (Neper/km) ; T = temperature (K) ; P = pressure (mb) ; f = frequency (GHz) (valid 0-1000 GHz) ; ; 5/22/02 P.Rosenkranz ; ; Equations based on: ; Borysow, A, and L. Frommhold, ; Astrophysical Journal, v.311, pp.1043-1057 (1986) ; with modification of 1.29 to account for O2-O2 and O2-N2 ; collisions, as suggested by ; J.R. Pardo, E.Serabyn, J.Cernicharo, J. Quant. Spectros. ; Radiat. Trans. v.68, pp.419-433 (2001). ; ; Conversion to Fortran 90 (9/21/02) by M. Walters ; Converted to IDL (9/03) by C. O'Dell FDEPEN = .5 + .5 / (1. + (f / 450.)^2) alpha = 1.29 * 6.5E-14 * FDEPEN * P^2 * f^2 * (300./T)^3.6 return, alpha end ; function abs_n2 function abs_o2,TEMP, PRES, VAPDEN, FREQ ; PURPOSE: RETURNS ABSORPTION COEFFICIENT DUE TO OXYGEN IN AIR, ; IN NEPERS/KM ; ; 5/1/95 P. Rosenkranz ; 11/5/97 P. Rosenkranz - 1- line modification. ; 12/16/98 PWR - updated submm freq's and intensities from HITRAN96 ; 8/21/02 PWR - shift at 118 and revised width at 425 ; ; 9/21/02 M. Walters - coverted to Fortran 90 ; Converted to IDL (9/03) by C. O'Dell ; ARGUMENTS: ; ; NAME UNITS DESCRIPTION VALID RANGE ;---------------------------------------------------------------------------- ; TEMP Kelvin temperature UNCERTAIN, but believed to be ; valid for atmosphere ; PRES millibars pressure 3 TO 1000 ; VAPDEN g/m**3 water vapor density (enters linewidth calculation due to ; greater broadening efficiency of H2O) ; FREQ GHz frequency 0 TO 900 ; ; ; REFERENCES FOR EQUATIONS AND COEFFICIENTS: ; P.W. Rosenkranz, CHAP. 2 and appendix, in ATMOSPHERIC REMOTE SENSING ; BY MICROWAVE RADIOMETRY (M.A. Janssen, ed., 1993). ; HJ. Liebe et al, JQSRT V.48, pp.629-643 (1992). ; M.J. Schwartz, Ph.D. thesis, M.I.T. (1998). ; M.Yu. Tretyakov et al, J. Mol. Spect. v.208, pp.110-112 (2001). ; A.F. Krupnov et al, J. Mol. Spect. (2002). ; SUBMILLIMETER LINE INTENSITIES FROM HITRAN96. ; ; This version differs from Liebe's MPM92 in these significant respects: ; 1. It uses the modification of the 1- line width temperature dependence ; recommended by Schwartz: (1/T). ; 2. The 1- line has the shift measured by Tretyakov et al. Temperature ; dependence of shift is unknown. ; 3. It uses the same temperature dependence (X) for submillimeter line ; widths as in the 60 GHz band: (1/T)**0.8 ; 4. The 425 GHz line width is from Krupnov et al. ;---------------------------------------------------------------------------- ; Local variables: ; real :: TH, TH1, B, PRESWV, PRESDA, DEN, DENS, DFNR, DF ; real :: SUM, STR, Y, SF1, SF2, FCEN WB300 = 0.56 X = 0.8 ; real, save :: W300(40), Y300(40), V(40), BE(40) ; integer :: K ; Lines are arranged 1-,1+,3-,3+,etc. in spin-rotation spectrum F = [ $ 118.7503, 56.2648, 62.4863, 58.4466, 60.3061, 59.5910, $ 59.1642, 60.4348, 58.3239, 61.1506, 57.6125, 61.8002, $ 56.9682, 62.4112, 56.3634, 62.9980, 55.7838, 63.5685, $ 55.2214, 64.1278, 54.6712, 64.6789, 54.1300, 65.2241, $ 53.5957, 65.7648, 53.0669, 66.3021, 52.5424, 66.8368, $ 52.0214, 67.3696, 51.5034, 67.9009, 368.4984, 424.7632, $ 487.2494, 715.3931, 773.8397, 834.1458 ] S300 = [ $ .2936E-14, .8079E-15, .2480E-14, .2228E-14, $ .3351E-14, .3292E-14, .3721E-14, .3891E-14, $ .3640E-14, .4005E-14, .3227E-14, .3715E-14, $ .2627E-14, .3156E-14, .1982E-14, .2477E-14, $ .1391E-14, .1808E-14, .9124E-15, .1230E-14, $ .5603E-15, .7842E-15, .3228E-15, .4689E-15, $ .1748E-15, .2632E-15, .8898E-16, .1389E-15, $ .4264E-16, .6899E-16, .1924E-16, .3229E-16, $ .8191E-17, .1423E-16, .6494E-15, .7083E-14, $ .3025E-14, .1835E-14, .1158E-13, .3993E-14 ] BE = [ .009, .015, .083, .084,.212,.212,.391,.391, .626,.626, $ .915,.915, 1.260,1.260, 1.660, 1.665, 2.119, 2.115, 2.624, 2.625, $ 3.194,3.194, 3.814, 3.814, 4.484, 4.484, 5.224, 5.224, $ 6.004, 6.004, 6.844, 6.884, $ 7.744, 7.774, .048, .044, .049, .145, .141, .145 ] ; Widths in MHz/mb W300 = [ 1.63, 1.646, 1.468, 1.449, 1.382, 1.360, $ 1.319, 1.297, 1.266, 1.248, 1.221, 1.207, 1.181, 1.171, $ 1.144, 1.139, 1.110, 1.108, 1.079, 1.078, 1.05, 1.05, $ 1.02, 1.02, 1.00, 1.00, .97, 0.97, .94, .94, .92, .92, .89, .89, $ 1.64, 1.64, 1.64, 1.81, 1.81, 1.81 ] Y300 = [ -0.0233, 0.2408, -0.3486, 0.5227, $ -0.5430, 0.5877, -0.3970, 0.3237, -0.1348, 0.0311, $ 0.0725, -0.1663, 0.2832, -0.3629, 0.3970, -0.4599, $ 0.4695, -0.5199, 0.5187, -0.5597, 0.5903, -0.6246, $ 0.6656, -0.6942, 0.7086, -0.7325, 0.7348, -0.7546, $ 0.7702, -0.7864, 0.8083, -0.8210, 0.8439, -0.8529, $ 0., 0., 0., 0., 0., 0. ] V = [ 0.0079, -0.0978, 0.0844, -0.1273, $ 0.0699, -0.0776, 0.2309, -0.2825, 0.0436, -0.0584, $ 0.6056, -0.6619, 0.6451, -0.6759, 0.6547, -0.6675, $ 0.6135, -0.6139, 0.2952, -0.2895, 0.2654, -0.2590, $ 0.3750, -0.3680, 0.5085, -0.5002, 0.6206, -0.6091, $ 0.6526, -0.6393, 0.6640, -0.6475, 0.6729, -0.6545, $ 0., 0., 0., 0., 0., 0. ] TH = 300. / TEMP B = TH^X PRESWV = VAPDEN * TEMP / 217. PRESDA = PRES - PRESWV DEN = .001 * (PRESDA * B + 1.1 * PRESWV * TH) DENS = .001 * (PRESDA + 1.1 * PRESWV) * TH DFNR = WB300 * DEN SUM = 1.6E-17 * FREQ*FREQ*DFNR / (TH * (FREQ*FREQ + DFNR*DFNR)) for k = 0, n_elements(F)-1 do begin ; cycle through our 40 lines ; exception for 1- line if (k EQ 0) then begin DF = W300[0] * DENS FCEN = F[0] - .14 * DENS endif else begin DF = W300[k] * DEN FCEN = F[k] endelse Y = .001 * PRES * B * (Y300[k] + V[k] * (TH-1)) STR = S300[k] * exp( - BE[k] * (TH-1)) SF1 = (DF + (FREQ - FCEN) * Y) / ((FREQ - FCEN)^2 + DF^2) SF2 = (DF - (FREQ + FCEN) * Y) / ((FREQ + FCEN)^2 + DF^2) SUM = SUM + STR * (SF1 + SF2) * (FREQ / F[k] )^2 endfor alpha = 0.5034E12 * SUM * PRESDA * TH^3 / !pi return, alpha end ; function abs_o2 ;---------------------------------------------------------------- function abs_h2o, T, P, RHO, F ; Purpose: compute absorption coef in atmosphere due to water vapor ; ; Calling sequence parameters: ; ; NAME UNITS I/O DESCRIPTON VALID RANGE ;-------------------------------------------------------------- ; T Kelvin I temperature ; P millibar I pressure .1 to 1000 ; RHO g/m**3 I water vapor density ; F GHz I frequency 0 to 800 ; alpha nepers/km O absorption coefficient ; ; References: ; P.W. Rosenkranz, RADIO SCIENCE V.33, PP.919-928 (1998); V.34, P.1025 ; ; Line intensities selection threshold = ; half of continuum absorption at 1000 mb. ; Widths measured at 22,183,380 GHz, others calculated. ; A.BAUER ET AL.ASA WORKSHOP (SEPT. 1989) (380GHz). ; ; Revision history: ; DATE- Oct. 6, 1988 P. W. Rosenkranz - eqs as publ. in 1993. ; Oct. 4, 1995 PWR - use clough's definition of local line ; contribution, hitran intensities, add 7 lines. ; Oct. 24, 1995 PWR - add 1 line. ; July 7, 1997 PWR - separate coeff. for self-broadening, ; revised continuum. ; Aug. 28, 2002 PWR - corrected line intensities ; Sep. 21, 2002 M. Walters - convert to Fortran 90 ; Sep. 18, 2003 Conversion to IDL by C. O'Dell ; Local variables: ; integer :: I, J ; real :: PVAP, PDA, DEN, TI, TI2, SUM, WIDTH, WSQ, S, BASE, RES, CON ; Line frequencies: FL = [ $ 22.2351, 183.3101, 321.2256, 325.1529, 380.1974, 439.1508, $ 443.0183, 448.0011, 470.8890, 474.6891, 488.4911, 556.9360, $ 620.7008, 752.0332, 916.1712 ] ; Line intensities at 300K: S1 = [ $ .1314E-13, .2279E-11, .8058E-13, .2701E-11, .2444E-10, $ .2185E-11, .4637E-12, .2568E-10, .8392E-12, .3272E-11, $ .6676E-12, .1535E-08, .1711E-10, .1014E-08, .4238E-10 ] ; T coeff. of intensities: B2 = [ $ 2.144, .668, 6.179, 1.541, 1.048, 3.595, 5.048, 1.405, $ 3.597, 2.379, 2.852, .159, 2.391, .396, 1.441 ] ; air-broadened width parameters at 300K: W3 = [ $ .00281, .00281, .0023, .00278, .00287, $ .0021, .00186, .00263, .00215, .00236, $ .0026, .00321, .00244, .00306, .00267 ] ; T-exponent of air-broadening: X = [ .69, .64, .67, .68, .54, .63, .60, $ .66, .66, .65, .69, .69, .71, .68, .70 ] ; self-broadened width parameters at 300K: WS = [ $ .01349, .01491, .0108, .0135, .01541, $ .0090, .00788, .01275, .00983, .01095, $ .01313, .01320, .01140, .01253, .01275 ] ; T-exponent of self-broadening: XS = [ .61, .85, .54, .74, .89, .52, .50, $ .67, .65, .64, .72, 1.0, .68, .84, .78 ] ; determine output size and shape alpha = f*0 + t*0 + p*0 + rho*0 ; get the parts where there is no absorption (b/c no water vapor!) wle0 = where(RHO LE 0) wgt0 = where(RHO GT 0) if wgt0[0] ne -1 then begin PVAP = RHO * T / 217. ; const includes isotopic abundance TI = 300. / T ; continuum terms CON = (5.43E-10 * (P-PVAP) * TI^3 + 1.8E-8 * PVAP * TI^7.5) * PVAP*F^2 ; add resonances SUM = alpha*0. for i = 0, n_elements(fl)-1 do begin WIDTH = W3[i] * (P-PVAP) * TI^X[i] + WS[i] * PVAP * TI^XS[i] WSQ = WIDTH * WIDTH S = S1[i] * (TI^2.5) * exp(B2[i] * (1. - TI)) DF1 = F - FL[i] DF2 = F + FL[i] ; use clough's definition of local line contribution BASE = WIDTH / (562500. + WSQ) ; do for positive and negative resonance RES1 = alpha*0. RES2 = RES1 if n_elements(f) eq 1 then begin ; the one frequency case if abs(df1) LT 750. then RES1 = RES1 + WIDTH/(DF1^2 + WSQ) - BASE if abs(df2) LT 750. then RES2 = RES2 + WIDTH/(DF2^2 + WSQ) - BASE endif else begin ; the multifrequency case wdf1 = where(abs(df1) LT 750.) wdf2 = where(abs(df2) LT 750.) if wdf1[0] ne -1 then RES1[wdf1] = WIDTH/(DF1[wdf1]^2 + WSQ) - BASE if wdf2[0] ne -1 then RES2[wdf2] = WIDTH/(DF2[wdf2]^2 + WSQ) - BASE endelse SUM = SUM + S * (RES1+RES2) * (F/FL[i])^2 endfor alpha = .3183E-4 * (3.335E16*RHO) * SUM + CON endif if wle0[0] ne -1 then alpha[wle0] = 0. return, alpha end ; function abs_h2o PRO gasabsr02, f, TK, rhowv, Pa, absair, abswv, bad=bad, volume=volume ; INPUT VARIABLES ; F = frequency (GHz), ; Tk = absolute temperature (K) ; Rhowv = water vapor density (kg/m**3). ; Pa = Total air pressure (Pascals). ; *note. F must be a scalar. Tk, Rhowv, and Pa must ; be arrays of the same size and shape. ; OUTPUT VARS ; ; Absair = extinction by dry air (meters squared per kg *moist air*) ; Abswv = extinction by water vapor (meters squared per kg water vapor) ; *note: output arrays will be the same size and shape as Tk, Rhowv, Pa. ; real :: Pmb, vapden ; real :: e, q, Tv, rhoair ; real :: abs_n2, abs_o2, abs_h2o ; KEYWORDS ; bad : set this to an unnamed variable will return the indices ; where the input data are out of bounds, and hence ; the calculations are meaningless. ; ; volume: set this keyword to report the output in m^-1 (volume extinction) ; as opposed to m^2 kg^-1 (mass extinction) ;-------------------------------------------------------------------------- ; check for "reasonable" input values if ((f LE 0) OR (F GT 800)) then begin print, 'Frequency out of range in gasabs02' STOP endif bad = where((Pa LT 10) OR (Pa GT 1e5)) ; still try to do calc on all elements! ; convert pressure from Pa to Mb Pmb = Pa / 100.0 ; convert vapor density from kg/m**3 to g/m**3 vapden = rhowv * 1000.0 absair = abs_n2(TK, Pmb, f) + abs_o2(TK, Pmb, vapden, f) abswv = abs_h2o(TK, Pmb, vapden, f) ; convert above from Np/km to m**2/kg ; convert vapor density to vapor pressure e = rhowv * (TK * 461.5) ; calculate specific humidity q = 0.622 * e / Pa ; calculate virtual temperature Tv = (1. + 0.61*q) * TK ; moist air density rhoair = Pa / (Tv * 287.06) ; convert from km^-1 to m^-1 absair = absair * 0.001 abswv = abswv * 0.001 if not keyword_set(volume) then begin absair = absair / rhoair abswv = abswv / rhowv endif END