Add files via upload

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   

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
+

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