From 29de2ffe03617a70246ce6edf84d50f7dbff89a8 Mon Sep 17 00:00:00 2001 From: "Dr. Thomas Meissner" <75815736+thomasmeiss@users.noreply.github.com> Date: Fri, 12 Aug 2022 16:55:59 -0700 Subject: [PATCH] Add files via upload --- MAKE_L_BAND_ATM.f90 | 105 +++++++++++++ column.f90 | 42 +++++ dielectric.f90 | 131 ++++++++++++++++ fdabscoeff_v2.f90 | 318 ++++++++++++++++++++++++++++++++++++++ fdcldabs.f90 | 23 +++ findncep_025deg.f90 | 366 ++++++++++++++++++++++++++++++++++++++++++++ get_atm.f90 | 139 +++++++++++++++++ goff_gratch_vap.f90 | 32 ++++ 8 files changed, 1156 insertions(+) create mode 100644 MAKE_L_BAND_ATM.f90 create mode 100644 column.f90 create mode 100644 dielectric.f90 create mode 100644 fdabscoeff_v2.f90 create mode 100644 fdcldabs.f90 create mode 100644 findncep_025deg.f90 create mode 100644 get_atm.f90 create mode 100644 goff_gratch_vap.f90 diff --git a/MAKE_L_BAND_ATM.f90 b/MAKE_L_BAND_ATM.f90 new file mode 100644 index 0000000..9f4aaca --- /dev/null +++ b/MAKE_L_BAND_ATM.f90 @@ -0,0 +1,105 @@ +! Thomas Meissner +! Remote Sensing Systems +! RSS L-band ATMOSPHERIC Absorption +! 08/11/2022 + +! MAKE_L_BAND_ATM +! Create ancillary fileds for AO, AV, AL, TRAN, TBUP, TBDW +! at SMAP frequency and relevant Earth Incidence (EIA) range +! ingest NCEP GFS/GDAS 0.25 deg surface and atmospheric profiles +! for temperature, height, relative humidity, cloud water mixing ratio. +! write out as raw binary without header +! no space or time interpolation. +! +! RSS-Compile: Menu -> Tools -> compile64 + +! internal version history +! Change ilat and ilon to go from 1-721 and 1-1440, respectively - AManaster 11/2021 +! v2: change back from Wentz-Meissner 2016 to Liebe O2 absorption + + +include 'fdabscoeff_v2.f90' +include 'dielectric.f90' +include 'fdcldabs.f90' +include 'goff_gratch_vap.f90' +include 'column.f90' +include 'get_atm.f90' +include 'findncep_025deg.f90' + + + + program MAKE_L_BAND_ATM + implicit none + + real(4), parameter :: freq=1.413 ! center frequency + integer(4), parameter :: icld=1 ! 1 = include coukld liquid water absorption + + character(len=200), parameter :: outpath = 'sample_data\' ! output directory. specified by user. + character(len=200) :: outfile, str + + real(4) :: xtbup(39:41), xtbdw(39:41), xtran(39:41) + real(4) :: xao, xav, xal + real(4) :: tbup(1:1440,1:721,39:41), tbdw(1:1440,1:721,39:41), tran(1:1440,1:721,39:41) + real(4) :: ao(1:1440,1:721), av(1:1440,1:721), al(1:1440,1:721) + + integer(4) :: lyear, imon, iday, ihour + integer(4) :: ilat, ilon + + + write(*,*) ' L-BAND ATM for SMAP' + write(*,*) ' create ATM ABSORPTION, TRAN, TBUP, TBDW form 0.25 DEG NCEP ATMOSPHERIC PROFILES and SURFACE.' + write(*,*) ' output on regular 0.25 DEG NCEP GRID.' + write(*,*) + + + write(*,*)' enter year' + read(*,*) lyear + + write(*,*)' enter month' + read(*,*) imon + + write(*,*)' enter day of month' + read(*,*) iday + + write(*,*)'enter hour of day of analysis time: 00, 06, 12, 18' + read(*,*) ihour + + write(*,*) + write(*,*) ' start processing' + + write(*,*) ' time stamp ',lyear,imon,iday,ihour + + do ilat=1,721 + write(*,*) ' lat ',ilat + + do ilon=1,1440 + + call get_atm(icld,lyear,imon,iday,ihour,ilat,ilon,freq, xtbup,xtbdw,xtran,xao,xav,xal) + + ao(ilon,ilat)=xao + av(ilon,ilat)=xav + al(ilon,ilat)=xal + + tran(ilon,ilat,:)=xtran + tbup(ilon,ilat,:)=xtbup + tbdw(ilon,ilat,:)=xtbdw + + enddo !ilon + enddo !ilat + + + write(str,9001) lyear,imon,iday,ihour + 9001 format('atmos_ncep_',i4.4,'_',i2.2,'_',i2.2,'_',i2.2,'z.dat') + outfile = trim(outpath)//trim(str) + + open(unit=4,form='binary',file=outfile,action='write') + write(4) ao + write(4) av + write(4) al + write(4) tran + write(4) tbup + write(4) tbdw + close(4) + + stop ' norm end pgm MAKE_L_BAND_ATM' + end program MAKE_L_BAND_ATM \ No newline at end of file diff --git a/column.f90 b/column.f90 new file mode 100644 index 0000000..0ccda61 --- /dev/null +++ b/column.f90 @@ -0,0 +1,42 @@ +! columnar integral +! input: +! nlevel: number of profiles (nlevel = 0 -> sfc) +! z: altitude levels [in m] +! rho: density [in unit/m] (unit arbitrary) +! ip: =1 : linear varying profile -> arithmetic mean for integration +! =2 : approximately exponentially varying profile -> +! use average of arithmetic and geometric mean for integration +! this has smaller error in integration than either +! arithmetic (a) or geometric (g) mean +! (if profile is varying exactly exponentially with z, then the integration error +! is minimal for the combination: 2/3 g + 1/3 a) +! output: +! col [in unit] + + subroutine column(nlevel,z,rho,ip, col) + implicit none + + integer(4) :: ip, nlevel, i + real(4) :: col, dz , avg + real(4) :: z(0:nlevel),rho(0:nlevel) + + col = 0. + do i=1,nlevel + + if (z(i).le.z(i-1)) stop 'error in column, pgm stopped' + + dz = z(i) - z(i-1) + + if (ip.eq.1) then + avg = 0.5* (rho(i) + rho(i-1) ) + else + avg = 0.25* ( rho(i) + rho(i-1) + 2.*sqrt(rho(i-1)*rho(i)) ) + endif + + col = col + avg*dz + + enddo + + return + end subroutine column + diff --git a/dielectric.f90 b/dielectric.f90 new file mode 100644 index 0000000..953964a --- /dev/null +++ b/dielectric.f90 @@ -0,0 +1,131 @@ +! complex dielectric constant: eps +! t. meissner, february 2002 +! updated oct 2004 + +! input: +! name parameter unit range +! freq frequency [ghz] 1 to 400 +! t sst [k] 248.16 k (-25 c) to 313.16 k (40 c) for pure water +! 271.16 k (-2 c) to 307.16 k (34 c) for saline water +! s salinity [ppt] 0 to 40 +! +! output: +! eps complex dielectric constant +! negative imaginary part to be consistent with wentz1 conventionc +! + subroutine meissner(freq,t,s, eps) + implicit none + + real(4), parameter :: f0=17.97510 + real(4) freq,t,sst,s + real(4) e0s,e1s,e2s,n1s,n2s,sig + complex(4) :: j = (0.0,1.0), eps + + sst = t - 273.15 ! [celsius] + call dielectric_meissner_wentz(sst,s, e0s,e1s,e2s,n1s,n2s,sig) + +! debye law (2 relaxation wavelengths) + eps = (e0s - e1s)/(1.0 - j*(freq/n1s)) + & + (e1s - e2s)/(1.0 - j*(freq/n2s)) + e2s + & + j*sig*f0/freq + + + eps = conjg(eps) + + return + end subroutine meissner + +! this routine is same as the one in o:\rtm\geomod9c.f dated June 22 2009 + +! this routine was carefully editted by fjw to separate out the sst,s dependence +! i also permultiplied sst**2, sst**3 etc. +! i numerically verified that the editted version gave exactly the same results are thomas orginal + +! complex dielectric constant: eps +! t. meissner + f. wentz, ieee tgars, 42(9), 2004, 1836-1849 +! +! input: +! name parameter unit range +! sst sst [c] -25 c to 40 c for pure water +! -2 c to 34 c for saline water +! s salinity [ppt] 0 to 40 +! +! output: +! eps complex dielectric constant +! negative imaginary part to be consistent with wentz1 convention + + subroutine dielectric_meissner_wentz(sst_in,s, e0s,e1s,e2s,n1s,n2s,sig) + implicit none + real(4), intent(in) :: sst_in,s + real(4), intent(out) :: e0s,e1s,e2s,n1s,n2s,sig + + real(4), dimension(11), parameter :: & + x=(/ 5.7230e+00, 2.2379e-02, -7.1237e-04, 5.0478e+00, -7.0315e-02, 6.0059e-04, 3.6143e+00, & + 2.8841e-02, 1.3652e-01, 1.4825e-03, 2.4166e-04 /) + + + real(4), dimension(13), parameter :: & + z=(/ -3.56417e-03, 4.74868e-06, 1.15574e-05, 2.39357e-03, -3.13530e-05, & + 2.52477e-07, -6.28908e-03, 1.76032e-04, -9.22144e-05, -1.99723e-02, & + 1.81176e-04, -2.04265e-03, 1.57883e-04 /) ! 2004 + + real(4), dimension(3), parameter :: a0coef=(/ -0.33330E-02, 4.74868e-06, 0.0e+00/) + real(4), dimension(5), parameter :: b1coef=(/0.23232E-02, -0.79208E-04, 0.36764E-05, -0.35594E-06, 0.89795E-08/) + + real(4) :: e0,e1,e2,n1,n2 + real(4) :: a0,a1,a2,b1,b2 + real(4) :: sig35,r15,rtr15,alpha0,alpha1 + + real(4) sst,sst2,sst3,sst4,s2 + + sst=sst_in + ! if(sst.lt.-25) sst=-25 !pretects against n1 and n2 going zero for very cold water + if(sst.lt.-30.16) sst=-30.16 !pretects against n1 and n2 going zero for very cold water + + sst2=sst*sst + sst3=sst2*sst + sst4=sst3*sst + + s2=s*s + + ! pure water + e0 = (3.70886e4 - 8.2168e1*sst)/(4.21854e2 + sst) ! stogryn et al. + e1 = x(1) + x(2)*sst + x(3)*sst2 + n1 = (45.00 + sst)/(x(4) + x(5)*sst + x(6)*sst2) + e2 = x(7) + x(8)*sst + n2 = (45.00 + sst)/(x(9) + x(10)*sst + x(11)*sst2) + + ! saline water + ! conductivity [s/m] taken from stogryn et al. + sig35 = 2.903602 + 8.60700e-2*sst + 4.738817e-4*sst2 - 2.9910e-6*sst3 + 4.3047e-9*sst4 + r15 = s*(37.5109+5.45216*s+1.4409e-2*s2)/(1004.75+182.283*s+s2) + + alpha0 = (6.9431+3.2841*s-9.9486e-2*s2)/(84.850+69.024*s+s2) + alpha1 = 49.843 - 0.2276*s + 0.198e-2*s2 + rtr15 = 1.0 + (sst-15.0)*alpha0/(alpha1+sst) + + sig = sig35*r15*rtr15 + + ! permittivity + a0 = exp(a0coef(1)*s + a0coef(2)*s2 + a0coef(3)*s*sst) + e0s = a0*e0 + + if(sst.le.30) then + b1 = 1.0 + s*(b1coef(1) + b1coef(2)*sst + b1coef(3)*sst2 + b1coef(4)*sst3 + b1coef(5)*sst4) + else + b1 = 1.0 + s*(9.1873715e-04 + 1.5012396e-04*(sst-30)) + endif + n1s = n1*b1 + + a1 = exp(z(7)*s + z(8)*s2 + z(9)*s*sst) + e1s = e1*a1 + + ! b2 = 1.0 + s*(z(10) + z(11)*sst) + b2 = 1.0 + s*(z(10) + 0.5*z(11)*(sst + 30)) + n2s = n2*b2 + + a2 = 1.0 + s*(z(12) + z(13)*sst) + e2s = e2*a2 + + return + end subroutine dielectric_meissner_wentz diff --git a/fdabscoeff_v2.f90 b/fdabscoeff_v2.f90 new file mode 100644 index 0000000..fbf70e4 --- /dev/null +++ b/fdabscoeff_v2.f90 @@ -0,0 +1,318 @@ +! TM 09/19/2016 +! changed O2 absorption back to Liebe 1992 +! +! same as 'O:\skytemp3\atms_abs_routines.f' dated July 17 2015 except +! 1. routine fdabscoeff is added +! 2. fdaray is removed +! 3. ga=5.6e-3*(pdry+1.1*pwet)*1.097687*tht**0.8 !1.102393=(300./262.6)**0.7, 262.6 is avg teff from zonal_adjustment.m + + +! input: +! oxygen absorption from rosenkranz +! freq frequency [in ghz] +! p pressure [in h pa] +! t temperature [in k] +! pv water vapor pressure [in hpa] +! +! output: +! av water vapor absorption coefficients [neper/km] +! ao oxygen absortption coefficient [neper/km] + + subroutine fdabscoeff(freq,p,t,pv, av,ao) + implicit none + + real(4), parameter :: xnaper=0.2302585094 !convert db/km to naper/km + + real(4) freq,p,t,pv + real(4) av,ao + real(4) gamoxy,gamh2o + + call fdabsoxy_1992_modified(p,t,pv,freq, gamoxy) !gamoxy is db/km + call abh2o_rk_modified( p,t,pv,freq, gamh2o) !gamh2o is db/km + + ao=xnaper*gamoxy + av=xnaper*gamh2o + + return + end subroutine fdabscoeff + + +! june 11 2015 changed july 17 2015. oxyopc adjustment above 37 ghz changed slightly. see 'O:\gmi\abs_cal\memo20.txt' + +! this module contains the oxygen absoprtion routine, the vapor absorption routine and fdaray +! there were used for the april-june 2009 skytemp update. see 'memo8.txt' + + +! ================================================================================================================ +! ========================== modified version of Liebe 1992 oxygen model ========================================= +! ================================================================================================================ + +! This is from Atmospheric 60-GHz Oxygen Spectrum:.. Liebe, Rosenkranz, Hufford, 1992 +! coded: June 2009 1992 by f.wentz +! inputs: t, temperature (k) +! p, total pressure (mb) +! pv, water vapor pressure (mb) +! freq, frequency (ghz) +! output: gamoxy, oxygen absorption coefficient (db/km) + + +! It is the same as fdabsoxy_1989 except for the a5 and a6 coefs for finding delta have different values. +! Also this 1992 versions says delta is proprotional to total pressure p rather than pdry. +! Also note in this routine, the 1.e-3 scaling is done in the startup block. + +! compared to abo2_rk (the code Rosenkranz sent us), this routine gives very similar results if you set the apterm to 0. +! for my freqs 6-85 ghz, the largest dif was 0.0003 at the coldest vapor at 85.5 GHz. +! the apterm adds 0.003 at 85 ghz +! Apart from the apterm, you can essentially say fdabsoxy_1992 and abo2_rk are the same for my purposes. + +! this routine has been modified in the following ways (see 'memo8.txt') +! 1. non-resonance continuum temperature coef changed from 0.8 to 1.5 +! 2. a p*p continuum was added +! these modifications were done June 22 2009 + + subroutine fdabsoxy_1992_modified(p,t,pv,freq, gamoxy) + implicit none + + integer(4), parameter :: nlines=44 + integer(4) istart,i + real(4) p,t,pv,freq,gamoxy + real(4) tht,pwet,pdry,ga,gasq,delta,rnuneg,rnupos,ff,zterm,apterm,sftot,xterm + real(4) h(6,nlines),f0(nlines),a1(nlines),a2(nlines),a3(nlines),a4(nlines),a5(nlines),a6(nlines) + + real(8) sum + + data istart/1/ + data a4/38*0., 6*0.6/ + +! freq a1 a2 a3 a5 a6 + data h/ & + 50.474238, 0.94e-6, 9.694, 8.60e-3, 0.210, 0.685, & + 50.987749, 2.46e-6, 8.694, 8.70e-3, 0.190, 0.680, & + 51.503350, 6.08e-6, 7.744, 8.90e-3, 0.171, 0.673, & + 52.021410, 14.14e-6, 6.844, 9.20e-3, 0.144, 0.664, & + 52.542394, 31.02e-6, 6.004, 9.40e-3, 0.118, 0.653, & + 53.066907, 64.10e-6, 5.224, 9.70e-3, 0.114, 0.621, & + 53.595749, 124.70e-6, 4.484, 10.00e-3, 0.200, 0.508, & + 54.130000, 228.00e-6, 3.814, 10.20e-3, 0.291, 0.375, & + 54.671159, 391.80e-6, 3.194, 10.50e-3, 0.325, 0.265, & + 55.221367, 631.60e-6, 2.624, 10.79e-3, 0.224, 0.295, & + 55.783802, 953.50e-6, 2.119, 11.10e-3, -0.144, 0.613, & + 56.264775, 548.90e-6, 0.015, 16.46e-3, 0.339, -0.098, & + 56.363389, 1344.00e-6, 1.660, 11.44e-3, -0.258, 0.655, & + 56.968206, 1763.00e-6, 1.260, 11.81e-3, -0.362, 0.645, & + 57.612484, 2141.00e-6, 0.915, 12.21e-3, -0.533, 0.606, & + 58.323877, 2386.00e-6, 0.626, 12.66e-3, -0.178, 0.044, & + 58.446590, 1457.00e-6, 0.084, 14.49e-3, 0.650, -0.127, & + 59.164207, 2404.00e-6, 0.391, 13.19e-3, -0.628, 0.231, & + 59.590983, 2112.00e-6, 0.212, 13.60e-3, 0.665, -0.078, & + 60.306061, 2124.00e-6, 0.212, 13.82e-3, -0.613, 0.070, & + 60.434776, 2461.00e-6, 0.391, 12.97e-3, 0.606, -0.282, & + 61.150560, 2504.00e-6, 0.626, 12.48e-3, 0.090, -0.058, & + 61.800154, 2298.00e-6, 0.915, 12.07e-3, 0.496, -0.662, & + 62.411215, 1933.00e-6, 1.260, 11.71e-3, 0.313, -0.676, & + 62.486260, 1517.00e-6, 0.083, 14.68e-3, -0.433, 0.084, & + 62.997977, 1503.00e-6, 1.665, 11.39e-3, 0.208, -0.668, & + 63.568518, 1087.00e-6, 2.115, 11.08e-3, 0.094, -0.614, & + 64.127767, 733.50e-6, 2.620, 10.78e-3, -0.270, -0.289, & + 64.678903, 463.50e-6, 3.195, 10.50e-3, -0.366, -0.259, & + 65.224071, 274.80e-6, 3.815, 10.20e-3, -0.326, -0.368, & + 65.764772, 153.00e-6, 4.485, 10.00e-3, -0.232, -0.500, & + 66.302091, 80.09e-6, 5.225, 9.70e-3, -0.146, -0.609, & + 66.836830, 39.46e-6, 6.005, 9.40e-3, -0.147, -0.639, & + 67.369598, 18.32e-6, 6.845, 9.20e-3, -0.174, -0.647, & + 67.900867, 8.01e-6, 7.745, 8.90e-3, -0.198, -0.655, & + 68.431005, 3.30e-6, 8.695, 8.70e-3, -0.210, -0.660, & + 68.960311, 1.28e-6, 9.695, 8.60e-3, -0.220, -0.665, & + 118.750343, 945.00e-6, 0.009, 16.30e-3, -0.031, 0.008, & + 368.498350, 67.90e-6, 0.049, 19.20e-3, 0.0, 0.0, & + 424.763124, 638.00e-6, 0.044, 19.16e-3, 0.0, 0.0, & + 487.249370, 235.00e-6, 0.049, 19.20e-3, 0.0, 0.0, & + 715.393150, 99.60e-6, 0.145, 18.10e-3, 0.0, 0.0, & + 773.839675, 671.00e-6, 0.130, 18.10e-3, 0.0, 0.0, & + 834.145330, 180.00e-6, 0.147, 18.10e-3, 0.0, 0.0/ + + if(istart.eq.1) then + istart=0 + f0(:)=h(1,:) + a1(:)=h(2,:)/h(1,:) + a2(:)=h(3,:) + a3(:)=h(4,:) + a5(:)=0.001*h(5,:) + a6(:)=0.001*h(6,:) + endif + + tht = 300/t + pwet=0.1*pv + pdry=0.1*p-pwet + xterm=1-tht + + sum = 0. + do i=1,nlines + ga = a3(i)*(pdry*tht**(0.8-a4(i)) + 1.1*tht*pwet) + gasq=ga*ga + delta=(a5(i) + a6(i)*tht)*p*tht**0.8 + rnuneg = f0(i)-freq + rnupos = f0(i)+freq + ff = (ga-rnuneg*delta)/(gasq+rnuneg**2) + (ga-rnupos*delta)/(gasq+rnupos**2) + sum = sum + ff*a1(i)*exp(a2(i)*xterm) + enddo + if(sum.lt.0) sum=0 + +! add nonresonant contribution + +! ga=5.6e-3*(pdry+1.1*pwet)*tht**0.8 +!x ga=5.6e-3*(pdry+1.1*pwet)*tht**1.5 !modification 1 +! changed O2 absorption back to Liebe + ga=5.6e-3*(pdry+1.1*pwet)*1.097687*tht**0.8 + !1.102393=(300./262.6)**0.7, 262.6 is avg teff from zonal_adjustment.m + + zterm=ga*(1.+(freq/ga)**2) + apterm=1.4e-10*(1-1.2e-5*freq**1.5)*pdry*tht**1.5 + if(apterm.lt.0) apterm=0 + sftot=pdry*freq*tht**2 * (tht*sum + 6.14e-4/zterm + apterm) + + gamoxy=0.1820*freq*sftot +!x if(freq.gt.37) gamoxy=gamoxy + 0.1820*43.e-10 *pdry**2*tht**3*(freq-37.)**1.7 !prior to 7/17/2015 + if(freq.gt.37) gamoxy=gamoxy + 0.1820*26.e-10 *pdry**2*tht**3*(freq-37.)**1.8 !implemented 7/17/2015. + + return + end subroutine fdabsoxy_1992_modified + + +! ================================================================================================================ +! ========================== modified version of Rosenkranz water vapor model ==================================== +! ================================================================================================================ + + +! purpose- compute absorption coef in atmosphere due to water vapor +! +! calling sequence parameters- +! specifications +! name units i/o descripton valid range +! +! t kelvin i temperature +! p millibar i pressure .1 to 1000 +! f ghz i frequency 0 to 800 +! gamh2o db/km o absorption coefficient +! +! references- +! p.w. rosenkranz, radio science v.33, pp.919-928 (1998). +! +! 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, 95 pwr -add 1 line. +! july 7, 97 pwr -separate coeff. for self-broadening, +! revised continuum. +! dec. 11, 98 pwr - added comments + +! the routine is a modified version of abh2o_rk_reformat. +! this routine has been modified in the following three ways (see 'memo8.txt') +! 1. b1(1)=1.01*b1(1) :22 ghz line strength increase slightly +! 2. 22 ghz line shape below 22 ghz has been modified +! 3. foreign and self broadening continuum has been adjusted +! these modification were done June 22 2009 + + + + + + subroutine abh2o_rk_modified(p,t,pv,freq, gamh2o) + implicit none + + integer(4), parameter :: nlines=15 + + integer(4) istart,i + real(4) t,p,freq + real(4) b1(nlines),b2(nlines),b3(nlines),f0(nlines),b4(nlines),b5(nlines),b6(nlines) + real(4) pv,s,base,gamh2o + real(4) tht,pwet,pdry,ga,gasq,sftot,xterm,rnuneg,rnupos + real(8) sum + + real(4) chi,chisq,freqsq,f0sq,u + + data istart/1/ + +! line frequencies: + data f0/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: + data b1/ .1310e-13, .2273e-11, .8036e-13, .2694e-11, .2438e-10, & + .2179e-11, .4624e-12, .2562e-10, .8369e-12, .3263e-11, .6659e-12, & + .1531e-08, .1707e-10, .1011e-08, .4227e-10/ +! t coeff. of intensities: + data 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: + data b3/.0281, .0281, .023, .0278, .0287, .021, .0186, .0263, .0215, .0236, .026, .0321, .0244, .0306, .0267/ +! self-broadened width parameters at 300k: + data b5/.1349, .1491, .108, .135, .1541, .090, .0788, .1275, .0983, .1095, .1313, .1320, .1140, .1253, .1275/ +! t-exponent of air-broadening: + data b4 /.69, .64, .67, .68, .54, .63, .60, .66, .66, .65, .69, .69, .71, .68, .70/ +! t-exponent of self-broadening: + data b6/.61, .85, .54, .74, .89, .52, .50, .67, .65, .64, .72, 1.0, .68, .84, .78/ + + if(istart.eq.1) then + istart=0 + b1=1.8281089E+14*b1/f0**2 + b5=b5/b3 !convert b5 to Liebe notation + b1(1)=1.01*b1(1) !modification 1 + endif + + if(pv.le.0.) then + gamh2o=0 + return + endif + + + pwet=0.1*pv + pdry=0.1*p-pwet + tht = 300./t + xterm=1-tht + freqsq=freq*freq + + sum = 0. + do i=1,nlines + f0sq=f0(i)*f0(i) + ga=b3(i)*(pdry*tht**b4(i) + b5(i)*pwet*tht**b6(i)) + gasq = ga*ga + s = b1(i)*exp(b2(i)*xterm) + rnuneg = f0(i)-freq + rnupos = f0(i)+freq + base = ga/(562500. + gasq) !use clough's definition of local line contribution + + if(i.ne.1) then + if(abs(rnuneg).lt.750) sum = sum + s*(ga/(gasq + rnuneg**2) - base) + if(abs(rnupos).lt.750) sum = sum + s*(ga/(gasq + rnupos**2) - base) + + else + chi=0 + + if(freq.lt.19) then + u=abs(freq-19.)/16.5 + if(u.lt.0) u=0 + if(u.gt.1) u=1 + chi=ga*u*u*(3-2*u) !modification 2 + endif + + chisq=chi*chi + sum=sum + s*2*((ga-chi)*freqsq + (ga+chi)*(f0sq+gasq-chisq))/((freqsq-f0sq-gasq+chisq)**2 + 4*freqsq*gasq) + endif + + enddo + if(sum.lt.0) sum=0 + +!x sftot=pwet*freq*tht**3.5*(sum + 1.2957246e-6*pdry/tht**0.5 + 4.2952193e-5*pwet*tht**4) + sftot=pwet*freq*tht**3.5*(sum + 1.1*1.2957246e-6*pdry/tht**0.5 + 0.425*(freq**0.10)*4.2952193e-5*pwet*tht**4) !modification 3 + + gamh2o=0.1820*freq*sftot + + return + end subroutine abh2o_rk_modified + diff --git a/fdcldabs.f90 b/fdcldabs.f90 new file mode 100644 index 0000000..03c98b9 --- /dev/null +++ b/fdcldabs.f90 @@ -0,0 +1,23 @@ +! liquid cloud water absorption +! rayleigh +! freq: frequency [ghz] +! t: temperature [k] +! rhol: liquid cloud water density [g/m**3] +! output: +! al: cloud water absorption coefficient [neper/km] + + subroutine fdcldabs(freq,t,rhol, al) + implicit none + + real(4) :: freq,t,rhol,rhol0,al,wavlen,c=29.979,pi=3.14159 + complex(4) :: permit + + rhol0 = 1.0e-6*rhol ![g/cm**3] + call meissner(freq,t,0.0, permit) + wavlen = c/freq + al = (6.0*pi*rhol0/wavlen)*aimag((1.0-permit)/(2.0+permit)) ! nepers/cm + al = 1.0e5*al ! nepers/km + + return + end subroutine fdcldabs + diff --git a/findncep_025deg.f90 b/findncep_025deg.f90 new file mode 100644 index 0000000..df222c0 --- /dev/null +++ b/findncep_025deg.f90 @@ -0,0 +1,366 @@ +! 11/2021 (AManaster) +! Changed to read in 0.25 degree NCEP files. Changed 'ilon' and 'ilat' +! to indicies between 1-1440 and 1-721, respectively instead of integers +! representative of the actual lon and lat values (cannot for 0.25 degrees). + + +! 03/04/2013 +! changed path from Nasserver5 to ops1p +! +! 01/05/2006 +! was: if t(ipr) < 100 k : ierr=-10, bail out +! changed: if t(ipr) < 100 k, set rhocwat=0.0 +! +! 01/28/2003 +! rhol = liquid cloud water density +! +! ordered ncep profiles, surtep, wind, pwat ,colvap +! +! input: +! nmax: maximum number of levels +! lyear: long year (e.g. 1999) +! imon: month +! iday: day of the month +! ihour: uct hour +! ilat: latitude [between -90 and 90] +! ilon: longitude [between 0 and 360] + + +! output: +! colvap: ncep water vapor integrated [mm] +! colwat: ncep columnar liquid cloud water integrated [mm] +! pwat: ncep value for precipitable water [mm] +! cwat: ncep value for columnar cloud water (liquid + ice) [mm] +! p: air pressure profile [mb] (0:nmax) ordered +! t: air temperature profile [k] (0:nmax) +! z: elevation profile [m] +! pv: water vapor pressure profile [mb] (0:nmax) +! rhov: water vapor density [g /m**3] +! rhol: liquid water density [g /m**3] +! ibegin: index of surface level + + + subroutine findncep_025deg(lyear,imon,iday,ihour,ilat,ilon, colvap,colwat,pwat,cwat,p,t,pv,rhov,rhol,z,ibegin) + implicit none + + integer(4), parameter :: nmax =26 + real(4), parameter :: r_e = 6371000, rd=287.05, epsilon = 0.622 + + integer(4), intent(in) :: lyear,imon,iday,ihour,ilat,ilon + integer(4), intent(out) :: ibegin + + real(4), intent(out) :: colvap,colwat,pwat,cwat + real(4), dimension(0:nmax), intent(out) :: p,t,pv,rhov,rhol,z + + integer(4) :: ipr + real(4), dimension(0:nmax) :: rhocwat,rhoi,clwmr,hgt,rh + + real(4) :: p_sfc + real(4) :: p0,rhoair + + + p(0:nmax) = (/ 0., & + 1000.,975.,950.,925.,900.,850.,800.,750.,700.,650.,600., & + 550.,500.,450.,400.,350.,300.,250.,200.,150.,100., & + 70., 50., 30., 20. ,10. /) + + + call fnd_ncep_025deg(lyear,imon,iday,ihour,ilat,ilon, t,hgt,rh,clwmr,p_sfc,pwat,cwat) + + p0 = 0.01*p_sfc ! hpa + + ibegin=-1 + do ipr=1,nmax + if(p(ipr).le.p0) then !le is used rather than lt to be consistent with old sorting method + ibegin=ipr-1 + exit + endif + enddo + if(ibegin.eq.-1) stop 'ibegin.eq.-1 in findncep, pgm stopped' + + p( ibegin) =p0 + t( ibegin) =t(0) + hgt( ibegin) =hgt(0) + rh( ibegin) =rh(0) + clwmr(ibegin) =clwmr(0) + + + z = hgt * r_e / (r_e - hgt) + if(z(ibegin) .ge. z(ibegin+1)) z(ibegin) = z(ibegin+1) - 0.1 + +! transform rh -> water vapor pressure and density + call goff_gratch_vap(nmax,t,rh,p, pv,rhov) + +! convert ncep clwmr to rhocwat + clwmr(ibegin) = clwmr(ibegin+1) ! clwmr at surface + + do ipr = 0,nmax + if (t(ipr) >= 100) then + rhoair = p(ipr)/(rd*t(ipr)) ! density of dry air + rhoair = rhoair*(1.0 - (1.0-epsilon)*pv(ipr)/p(ipr)) + ! density of moist air (without cloud), unusual ncep definition for mixing ratio + + rhocwat(ipr) = 1.0e5*clwmr(ipr)*rhoair ! [g/m**3] + call cldwatice(t(ipr),rhocwat(ipr), rhol(ipr),rhoi(ipr)) + else ! unphysical temperature + rhocwat(ipr)=0.0 + rhol(ipr) = 0.0 + rhoi(ipr) = 0.0 + endif + + enddo + + where (rhov < 0) rhov=0.0 + where (rhol < 0) rhol=0.0 + + call column(nmax-ibegin,z(ibegin:nmax),rhov(ibegin:nmax),2, colvap) + call column(nmax-ibegin,z(ibegin:nmax),rhol(ibegin:nmax),1, colwat) + + colvap = colvap *1.e-3 ! g/m**2 -> mm = kg/m**2 + colwat = colwat *1.e-3 ! g/m**2 -> mm = kg/m**2 + + return + end subroutine findncep_025deg + + + + + subroutine cldwatice(t,rhoc, rhol,rhoi) +! t: temperature [k] +! rhoc: cloud water density [arbitrary unit] +! +! rhol : liquid part [same unit as rhoc] +! rhoi : ice part [same unit as rhoc] + implicit none + real, parameter :: twat = 273.15, tice = 253.15 + real, intent(in) :: t,rhoc + real(4), intent(out) :: rhol,rhoi + +! water or ice + if (t >= twat) then + rhol = rhoc ! all liquid water + rhoi = 0.0 + else if (t <= tice) then + rhol = 0.0 ! all ice + rhoi = rhoc + else + rhol = rhoc*(t-tice)/(twat-tice) + rhoi = rhoc - rhol + ! linear temperature interpolation in between + endif + + return + end subroutine cldwatice + +! 01/23/2006 +! kyle changed subroutine yread +! read fileheader between each profile level + +! +! ncep variables +! 1/4 deg resolution +! 4 times daily: 00z 06z 12z 18z +! +! thomas meissner : apr 1999 +! +! needs to be linked with time_routines.f +! +! +! varaible level (subfolder) +! sst sfc (0) +! t sfc (0) 2 m above ground +! rh sfc (0) +! hgt sfc (0) +! p sfc (0) +! pwat col +! cwat col +! t 1:nmax (prf) +! rh 1:nmax (prf) +! hgt 1:nmax (prf) +! clwmr 1:nmax (prf) +! +! named as year_month_day_xxz.dat : +! binary file: +! year month day hour (integer(4)) +! real*4 array +! grid: lat from +90 to -90 by -0.25 deg +! lon from 0 to 360 by 0.25 deg +! total size 4*4 + 360*181*4 = 260656 +! +! input: lyear (year, long form, e.g. 2005) integer +! imon (month, 1-12) integer +! iday (day of month, 1-31) integer +! ihour (hour of day) integer +! xlat (latitude) real(4) +! xlon (longitude) real(4) +! +! output: var (interpolated variable) real(4) +! +! + + subroutine fnd_ncep_025deg(lyear,imon,iday,ihour,ilat,ilon, t,hgt,rh,clwmr,p_sfc,pwat,cwat) + implicit none + + integer, parameter :: nmax =26 + + integer(4), intent(in) :: lyear,imon,iday,ihour,ilat,ilon + + real(4), intent(out) :: t(0:nmax),hgt(0:nmax),rh(0:nmax),clwmr(0:nmax),p_sfc,pwat,cwat + + integer(4), save :: lyearsv=-999,imonsv=-999,idaysv=-999,ihoursv=-999 + + integer(4) :: j,k + + real(4), save :: t_map(1440,721,0:nmax), hgt_map(1440,721,0:nmax), rh_map(1440,721,0:nmax),clwmr_map(1440,721,0:nmax) + real(4), save :: p_sfc_map(1440,721),pwat_map(1440,721),cwat_map(1440,721) + + if(lyear.ne.lyearsv .or. imon.ne.imonsv .or. iday.ne.idaysv .or. ihour.ne.ihoursv ) then + lyearsv=lyear + imonsv=imon + idaysv=iday + ihoursv=ihour + call read_ncep_025deg(lyear,imon,iday,ihour, t_map,hgt_map,rh_map,clwmr_map,p_sfc_map,pwat_map,cwat_map) + endif + + j= 722 - ilat + ! flips indicies around correctly since NCEP goes from +90 to -90 e.g., ilat = 1 (-90 on RSS grid) is NCEP idx 721 (-90 on NCEP grid) + k= ilon + + t= t_map(k,j,:) + hgt= hgt_map(k,j,:) + rh = rh_map(k,j,:) + clwmr=clwmr_map(k,j,:) + + p_sfc=p_sfc_map(k,j) + pwat = pwat_map(k,j) + cwat = cwat_map(k,j) + + return + end subroutine fnd_ncep_025deg + + + + subroutine read_ncep_025deg(lyear,imon,iday,ihour, t,hgt,rh,clwmr,p_sfc,pwat,cwat) + implicit none + + integer(4), parameter :: nmax = 26, nrh=21 + + integer(4), intent(in) :: lyear,imon,iday,ihour + real(4), intent(out) :: t(1440,721,0:nmax),hgt(1440,721,0:nmax),rh(1440,721,0:nmax),clwmr(1440,721,0:nmax) + real(4), intent(out) :: p_sfc(1440,721),pwat(1440,721),cwat(1440,721) + + character(len=200) :: filename + logical(4) :: lexist + + integer(4) :: kyear,kmon,kday,khour + integer(4) :: ilevel,icase,isleep + + ! location and name of NCEP files + ! specified by user + 9001 format('sample_data\NCEP_PROF_TMP_',i4.4,'_',i2.2,'_',i2.2,'_',i2.2,'z.dat') + 9002 format('sample_data\NCEP_SURF_TMP_2m_',i4.4,'_',i2.2,'_',i2.2,'_',i2.2,'z.dat') + 9003 format('sample_data\NCEP_PROF_RH_',i4.4,'_',i2.2,'_',i2.2,'_',i2.2,'z.dat') + 9004 format('sample_data\NCEP_SURF_RH_',i4.4,'_',i2.2,'_',i2.2,'_',i2.2,'z.dat') + 9005 format('sample_data\NCEP_PROF_HGT_',i4.4,'_',i2.2,'_',i2.2,'_',i2.2,'z.dat') + 9006 format('sample_data\NCEP_SURF_HGT_',i4.4,'_',i2.2,'_',i2.2,'_',i2.2,'z.dat') + 9007 format('sample_data\NCEP_SURF_PRES_',i4.4,'_',i2.2,'_',i2.2,'_',i2.2,'z.dat') + 9008 format('sample_data\NCEP_SURF_PWAT_',i4.4,'_',i2.2,'_',i2.2,'_',i2.2,'z.dat') + 9009 format('sample_data\NCEP_PROF_CLWMR_',i4.4,'_',i2.2,'_',i2.2,'_',i2.2,'z.dat') + 9010 format('sample_data\NCEP_SURF_CWAT_',i4.4,'_',i2.2,'_',i2.2,'_',i2.2,'z.dat') + + do icase=1,10 + if(icase.eq. 1) write(filename,9001) lyear,imon,iday,ihour + if(icase.eq. 2) write(filename,9002) lyear,imon,iday,ihour + if(icase.eq. 3) write(filename,9003) lyear,imon,iday,ihour + if(icase.eq. 4) write(filename,9004) lyear,imon,iday,ihour + if(icase.eq. 5) write(filename,9005) lyear,imon,iday,ihour + if(icase.eq. 6) write(filename,9006) lyear,imon,iday,ihour + if(icase.eq. 7) write(filename,9007) lyear,imon,iday,ihour + if(icase.eq. 8) write(filename,9008) lyear,imon,iday,ihour + if(icase.eq. 9) write(filename,9009) lyear,imon,iday,ihour + if(icase.eq.10) write(filename,9010) lyear,imon,iday,ihour + + isleep=1 + + 10 continue + inquire(file=filename,exist=lexist) + + if(.not.lexist) then + if(isleep.eq.1) write(*,1001) trim(filename) + 1001 format(' waiting for ',a100) + isleep=0 + call sleepqq(60000) !wait one minute + goto 10 + endif + + call sleepqq(10000) !wait 10 seconds for possible file download to complete + + open(unit=3,file=filename,status='old',form='binary',action='read') + + read(3) kyear, kmon, kday, khour + if(kyear.ne.lyear .or. kmon.ne.imon .or. kday.ne.iday .or. khour.ne.ihour) then + write(*,*) filename + write(*,*) kyear,kmon,kday,khour + write(*,*) lyear,imon,iday,ihour + write(*,*) ' header error in read_ncep_025deg, pgm stopped' + stop + endif + + + if(icase.eq.1) then + do ilevel=1,nmax + read(3) t(:,:,ilevel) + enddo + endif + + if(icase.eq.2) then + read(3) t(:,:,0) + endif + + if(icase.eq.3) then + do ilevel=1,nrh + read(3) rh(:,:,ilevel) + enddo + rh(:,:,nrh+1:nmax) = 0.0 + endif + + if(icase.eq.4) then + read(3) rh(:,:,0) + endif + + if(icase.eq.5) then + do ilevel=1,nmax + read(3) hgt(:,:,ilevel) + enddo + endif + + if(icase.eq.6) then + read(3) hgt(:,:,0) + endif + + if(icase.eq.7) then + read(3) p_sfc + endif + + if(icase.eq.8) then + read(3) pwat + endif + + if(icase.eq.9) then + do ilevel=1,nrh + read(3) clwmr(:,:,ilevel) + enddo + clwmr(:,:,nrh+1:nmax) = 0.0 + clwmr(:,:,0) = 0.0 + endif + + if(icase.eq.10) then + read(3) cwat + endif + + close(3) + + enddo !icase + + return + end subroutine read_ncep_025deg diff --git a/get_atm.f90 b/get_atm.f90 new file mode 100644 index 0000000..dcf8a06 --- /dev/null +++ b/get_atm.f90 @@ -0,0 +1,139 @@ +! adapted for SMAP +! EIA between 39 - 41 deg + + +! ncep 0.25 deg profiles +! +! rss vapor absorption +! liebe o2 absorption +! meissner - wentz dielectric model for cloud water +! +! icld = 0: no clouds +! icld = 1: ncep cloud density, only liquid part +! + + subroutine get_atm(icld,lyear,imonth,iday,ihour,ilat,ilon,freq, tbup,tbdw,tran,ao,av,al) + implicit none + + integer(4), parameter :: nmax = 26 + integer(4), intent(in) :: icld,lyear,imonth,iday,ihour,ilat,ilon + real(4), intent(in) :: freq + + real(4), intent(out) :: tbup(39:41),tbdw(39:41),tran(39:41) + real(4), intent(out) :: ao,av,al + + integer(4) :: ipr,ibegin,itht + real(4) :: tht,vap,cld,pwat,cwat + real(4) :: t(0:nmax),p(0:nmax),pv(0:nmax),z(0:nmax),rhov(0:nmax),rhol(0:nmax),rhol0(0:nmax) + real(4) :: abh2o(0:nmax),abo2(0:nmax),abcld(0:nmax),tabs(0:nmax) + + ! ncep parameters + call findncep_025deg(lyear,imonth,iday,ihour,ilat,ilon, vap,cld,pwat,cwat,p,t,pv,rhov,rhol0,z,ibegin) + + if (icld == 0) then ! no cloud + rhol= 0.0 + else + rhol= rhol0 + endif + + ! atmospheric absorption + ! o2 from rosenkranz + ! h2o vapor from amsr atbd + ! cloud with meissner dielectric constant + + + do ipr = ibegin,nmax + + call fdabscoeff(freq,p(ipr),t(ipr),pv(ipr), abh2o(ipr),abo2(ipr)) ! neper/km + + if (rhol(ipr) > 1.0e-7) then + call fdcldabs(freq,t(ipr),rhol(ipr), abcld(ipr)) ! nepers/km + else + abcld(ipr) = 0.0 + endif + + abh2o(ipr)=abh2o(ipr)/1000.0 ! neper/m + abo2(ipr) = abo2(ipr)/1000.0 ! neper/m + abcld(ipr)=abcld(ipr)/1000.0 ! neper/m + + enddo ! ipr + + ! vertical integrals + + call column(nmax-ibegin,z(ibegin:nmax),abh2o(ibegin:nmax),2, av) + call column(nmax-ibegin,z(ibegin:nmax), abo2(ibegin:nmax),2, ao) + call column(nmax-ibegin,z(ibegin:nmax),abcld(ibegin:nmax),1, al) + + ! total absorption + tabs(ibegin:nmax) = abh2o(ibegin:nmax) + abo2(ibegin:nmax) + abcld(ibegin:nmax) + + + do itht=39,41 + tht=float(itht) + call atm_tran(nmax-ibegin,tht,t(ibegin:nmax),z(ibegin:nmax),tabs(ibegin:nmax), tran(itht),tbdw(itht),tbup(itht)) + enddo + + return + end subroutine get_atm + + + ! compute atmospheric downwelling and upwelling brightness temperatures + ! and upward transmittance at each pressure level (altitude) + + ! input: + ! nlev number of atmosphere levels + ! tht earth incidence angle [in deg] + ! tabs(0:nlev) atmosphric absorptrion coefficients [nepers/m] + ! t(0:nlev) temperature profile[in k] + ! z(0:nlev) elevation (m) + + + ! output: + ! tran total atmospheric transmission + ! tbdw downwelling brightness temperature t_bd [in k] + ! tbup upwelling brightness temperature t_bu [in k] + + + subroutine atm_tran(nlev,tht,t,z,tabs, tran,tbdw,tbup) + implicit none + + real(4), parameter :: re=6378.135, delta=0.00035 + + integer(4), intent(in) :: nlev + real(4), intent(in) :: t(0:nlev),z(0:nlev),tabs(0:nlev) + real(4), intent(out) :: tran, tbdw, tbup + + real(4) :: opacty(nlev),tavg(nlev),ems(nlev) + real(4) :: sumop, sumdw, sumup, tbavg, dsdh, tht + + integer(4) :: i + + dsdh = (1.0+delta)/sqrt(cosd(tht)**2 + delta*(2+delta)) + + do i=1,nlev + opacty(i)=-dsdh*0.5*(tabs(i-1)+tabs(i))*(z(i)-z(i-1)) + tavg(i) =0.5*(t(i-1)+t(i)) + ems(i) =1.-exp(opacty(i)) + enddo + + sumop=0 + sumdw=0 + do i=1,nlev + sumdw=sumdw+(tavg(i)-t(1))*ems(i)*exp(sumop) + sumop=sumop+opacty(i) + enddo + + sumop=0 + sumup=0. + do i=nlev,1,-1 + sumup=sumup+(tavg(i)-t(1))*ems(i)*exp(sumop) + sumop=sumop+opacty(i) + enddo + + tran=exp(sumop) + tbavg=(1.-tran)*t(1) + tbdw=tbavg+sumdw + tbup=tbavg+sumup + + return + end subroutine atm_tran \ No newline at end of file diff --git a/goff_gratch_vap.f90 b/goff_gratch_vap.f90 new file mode 100644 index 0000000..9a76f47 --- /dev/null +++ b/goff_gratch_vap.f90 @@ -0,0 +1,32 @@ + subroutine goff_gratch_vap(nprofiles,T,RH,P, P_V,rho_V) + + +! calculates water vapor pressure P_V and density rho_V +! from relative humidity RH +! S. Cruz Pol, C. Ruf and S. Keihm, Radio Science 33 (5),1319 (1998) + +! water vapor pressure: P_v = P_s * RH +! water vapor density rho_v = F_w *(P_v * eps) / (R_d * T) +! + real, dimension(0:nprofiles) :: T,RH,P,P_v,rho_v,P_s + real, dimension(0:nprofiles) :: F_w,xi1,xi2,xi3 + real :: p_standard = 1013.246 ! hPa + real :: TS = 373.14 ! K (water boiling) + real :: T0 = 273.14 ! K + real :: eps = 1./1.607795 ! M(H2O)/M(air) + real :: R_d = 287.05 ! J /(K*kg) gas constant for dry air + + F_w = 1.0 + 1.E-4 * (5.92854 + 3.740346e-2 * P + & + 1.971198E-4 * (T-T0) * (800-P) + & + 6.045511E-6 * P * (T-T0)**2 ) ! deviation from ideal gas + + xi1 = -7.90298*(Ts/T - 1) + 5.02808*Log10(TS/T) + xi2 = -1.3816E-7 * 10**(11.344*(1.-T/TS) -1 ) + xi3 = 8.1328E-3 * (10**(-3.49149*(TS/T - 1)) - 1) + + P_s = p_standard * 10**(xi1 + xi2 + xi3) + P_v = P_s * RH* 0.01 ! mbar + rho_v = ((F_w * P_v * eps) / (R_d * T)) * 1.E5 ! [g/m**3] + + return + end subroutine goff_gratch_vap \ No newline at end of file