subroutine snwem_amsu(theta,frequency,snow_depth,ts,tba,tbb,esh,esv)
!$$$ subprogram documentation block
! . . . .
! subprogram: noaa/nesdis emissivity model over snow/ice for AMSU-A/B
!
! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
!
! abstract: noaa/nesdis emissivity model to compute microwave emissivity over
! snow for AMSU-A/B. The processing varies according to input parameters
! Option 1 : AMSU-A & B window channels of brightness temperatures (Tb)
! and surface temperature (Ts) are available
! Option 2 : AMSU-A window channels of Tb and Ts are available
! Option 3 : AMSU-A & B window channels of Tb are available
! Option 4 : AMSU-A window channels of Tb are available
! Option 5 : AMSU-B window channels of Tb and Ts are available
! Option 6 : AMSU-B window channels of Tb are available
! Option 7 : snow depth and Ts are available
!
! references:
! Yan, B., F. Weng and K.Okamoto,2004:
! "A microwave snow emissivity model, submitted to TGRS
!
! version: 3
!
! program history log:
! beta : November 28, 2000
!
! version 2.0: June 18, 2003.
!
! Version 2.0 enhances the capability/performance of beta version of
! land emissivity model (LandEM) over snow conditions. Two new subroutines
! (i.e., snowem_tb and six_indices) are added as replacements of the
! previous snow emissivity. If AMSU measurements are not available, the
! results are the same as these in beta version. The new snow emissivity
! model is empirically derived from satellite retrievals and ground-based
! measurements.
!
! version 3.0: August 18, 2003.
!
! Version 3.0 is an extended version of LandEM 2.0 over snow conditions.
! It covers seven different options (see below for details) for LandEM
! inputs over snow conditions. When All or limited AMSU measurements are
! available, one of the subroutines sem_ABTs, sem_ATs, sem_AB, sem_amsua,
! sem_BTs and sem_amsub, which are empirically derived from satellite
! retrievals and ground-based measurements, are called to simuate snow
! emissivity; when no AMSU measurements are avalaiable, the subroutine
! ALandEM_Snow is called where the results over snow conditions in beta
! version are adjusted with a bias correction that is obtained using a
! statistical algorithm. Thus, LandEM 3.0 significantly enhances the
! flexibility/performance of LandEM 2.0 in smulating emissivity over snow
! conditions.
!
!
! 2004-07-26 okamoto - modified the version 3.0 for use in gsi
! 2004-12-13 kokron - the 'depth' variable was declared then used before being given
! a value. changed to use snow_depth in all calculations
!
!
! input argument list:
! theta - local zenith angle in radian
! frequency - frequency in GHz
! t_skin - scattering layer temperature (K) (gdas)
! snow_depth - scatter medium depth (mm) (gdas)
! tba[1] ~ tba[4] - Tb at four AMSU-A window channels
! tba[1] : 23.8 GHz
! tba[2] : 31.4 GHz
! tba[3] : 50.3 GHz
! tba[4] : 89 GHz
! tbb[1] ~ tbb[2] - Tb at two AMSU-B window channels:
! tbb[1] : 89 GHz
! tbb[2] : 150 GHz
! When tba[ ] or tbb[ ] = -999.9: a missing value (no available data)
!
! output argument list:
! esv - emissivity at vertical polarization
! esh - emissivity at horizontal polarization
! snow_type - snow type (not output here)
! 1 : Wet Snow
! 2 : Grass_after_Snow
! 3 : RS_Snow (A)
! 4 : Powder Snow
! 5 : RS_Snow (B)
! 6 : RS_Snow (C)
! 7 : RS_Snow (D)
! 8 : Thin Crust Snow
! 9 : RS_Snow (E)
! 10: Bottom Crust Snow (A)
! 11: Shallow Snow
! 12: Deep Snow
! 13: Crust Snow
! 14: Medium Snow
! 15: Bottom Crust Snow (B)
! 16: Thick Crust Snow
! 999: AMSU measurements are not available or over non-snow conditions
! important internal variables/parameters:
!
! INDATA[1] ~ INDATA[7] - seven options calling different subroutines
! INDATA(1) = ABTs : call sem_ABTs
! INDATA(2) = ATs : call sem_ATs
! INDATA(3) = AMSUAB : call sem_AB
! INDATA(4) = AMSUA : call sem_amsua
! INDATA(5) = BTs : call sem_BTs
! INDATA(6) = AMSUB : call sem_amsub
! INDATA(7) = MODL : call ALandEM_Snow
! input_type - specific option index
! tb[1] ~ tb[5] - Tb at five AMSU-A & B window channels:
! tb[1] = tba[1]
! tb[2] = tba[2]
! tb[3] = tba[3]
! tb[4] = tba[4]
! tb[5] = tbb[2]
!
! remarks:
!
! Questions/comments: Please send to Fuzhong.Weng@noaa.gov or Banghua.Yan@noaa.gov
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$
! use kinds, only: r_kind,i_kind
implicit none
integer(i_kind) :: nch,nwcha,nwchb,nwch,nalg
Parameter(nch = 10, nwcha = 4, nwchb = 2, nwch = 5,nalg = 7)
real(r_kind) :: theta,frequency,ts,snow_depth
real(r_kind) :: em_vector(2),esv,esh
real(r_kind) :: tb(nwch),tba(nwcha),tbb(nwchb)
logical :: INDATA(nalg),AMSUAB,AMSUA,AMSUB,ABTs,ATs,BTs,MODL
integer(i_kind) :: snow_type,input_type,i,ich,np,k
Equivalence(INDATA(1), ABTs)
Equivalence(INDATA(2), ATs)
Equivalence(INDATA(3), AMSUAB)
Equivalence(INDATA(4), AMSUA)
Equivalence(INDATA(5), BTs)
Equivalence(INDATA(6), AMSUB)
Equivalence(INDATA(7), MODL)
! Initialization
call em_initialization(frequency,em_vector)
snow_type = -999
input_type = -999
do k = 1, nalg
INDATA(k) = .TRUE.
end do
! Read AMSU & Ts data and set available option
! Get five AMSU-A/B window measurements
tb(1) = tba(1); tb(2) = tba(2); tb(3) = tba(3); tb(4) = tba(4); tb(5) = tbb(2)
! Check available data
if((ts <= 100.0_r_kind) .or. (ts >= 320.0_r_kind) ) then
ABTs = .false.; ATs = .false.; BTs = .false.; MODL = .false.
end if
do i=1,nwcha
if((tba(i) <= 100.0_r_kind) .or. (tba(i) >= 320.0_r_kind) ) then
ABTs = .false.; ATs = .false.; AMSUAB = .false.; AMSUA = .false.
exit
end if
end do
do i=1,nwchb
if((tbb(i) <= 100.0_r_kind) .or. (tbb(i) >= 320.0_r_kind) ) then
ABTs = .false.; AMSUAB = .false.; BTs = .false.; AMSUB = .false.
exit
end if
end do
if((snow_depth < 0.0_r_kind) .or. (snow_depth >= 3000.0_r_kind)) MODL = .false.
if((frequency >= 80._r_kind) .and. (BTs)) then
ATs = .false.; AMSUAB = .false.
end if
! Check input type and call a specific Option/subroutine
do np = 1, nalg
if (INDATA(np)) then
input_type = np
exit
end if
end do
! write(UNIT=stdout,FMT='(a,2f6.1,i5,a,4f7.1,a,2f7.1)') 'freq,theta,input_tyep=',frequency,theta,input_type, ' tba=',tba,' tbb=',tbb
GET_option: SELECT CASE (input_type)
CASE (1)
call sem_ABTs(theta,frequency,tb,ts,snow_type,em_vector)
CASE (2)
call sem_ATs(theta,frequency,tba,ts,snow_type,em_vector)
CASE (3)
call sem_AB(theta,frequency,tb,snow_type,em_vector)
CASE (4)
call sem_amsua(theta,frequency,tba,snow_type,em_vector)
CASE(5)
call sem_BTs(theta,frequency,tbb,ts,snow_type,em_vector)
CASE(6)
call sem_amsub(theta,frequency,tbb,snow_type,em_vector)
CASE(7)
call ALandEM_Snow(theta,frequency,snow_depth,ts,snow_type,em_vector)
END SELECT GET_option
esv = em_vector(1)
esh = em_vector(2)
return
end subroutine snwem_amsu
!---------------------------------------------------------------------!
subroutine em_initialization(frequency,em_vector)
!$$$ subprogram documentation block
!
! subprogram: AMSU-A/B snow emissivity initialization
!
! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
!
! abstract: AMSU-A/B snow emissivity initialization
!
! program history log:
!
! input argument list:
!
! frequency - frequency in GHz
!
! output argument list:
!
! em_vector[1] and [2] - initial emissivity at two polarizations.
!
! important internal variables:
!
! freq[1~10] - ten frequencies for sixteen snow types of emissivity
! em[1~16,*] - sixteen snow emissivity spectra
! snow_type - snow type
! where it is initialized to as the type 4,i.e, Powder Snow
!
! remarks:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$
! use kinds, only: r_kind,i_kind
! use constants, only: one
implicit none
integer(i_kind) :: nch,ncand
Parameter(nch = 10,ncand=16)
real(r_kind) :: frequency,em_vector(*),freq(nch)
real(r_kind) :: em(ncand,nch)
real(r_kind) :: kratio, bconst,emissivity
integer(i_kind) :: snow_type,ich,k
save em
! Sixteen candidate snow emissivity spectra
data (em(1,k),k=1,nch)/0.87_r_kind,0.89_r_kind,0.91_r_kind,0.93_r_kind, &
0.94_r_kind,0.94_r_kind,0.94_r_kind,0.93_r_kind,0.92_r_kind,0.90_r_kind/
data (em(2,k),k=1,nch)/0.91_r_kind,0.91_r_kind,0.92_r_kind,0.91_r_kind, &
0.90_r_kind,0.90_r_kind,0.91_r_kind,0.91_r_kind,0.91_r_kind,0.86_r_kind/
data (em(3,k),k=1,nch)/0.90_r_kind,0.89_r_kind,0.88_r_kind,0.87_r_kind, &
0.86_r_kind,0.86_r_kind,0.85_r_kind,0.85_r_kind,0.82_r_kind,0.82_r_kind/
data (em(4,k),k=1,nch)/0.91_r_kind,0.91_r_kind,0.93_r_kind,0.93_r_kind, &
0.93_r_kind,0.93_r_kind,0.89_r_kind,0.88_r_kind,0.79_r_kind,0.79_r_kind/
data (em(5,k),k=1,nch)/0.90_r_kind,0.89_r_kind,0.88_r_kind,0.85_r_kind, &
0.84_r_kind,0.83_r_kind,0.83_r_kind,0.82_r_kind,0.79_r_kind,0.73_r_kind/
data (em(6,k),k=1,nch)/0.90_r_kind,0.89_r_kind,0.86_r_kind,0.82_r_kind, &
0.80_r_kind,0.79_r_kind,0.78_r_kind,0.78_r_kind,0.77_r_kind,0.77_r_kind/
data (em(7,k),k=1,nch)/0.88_r_kind,0.86_r_kind,0.85_r_kind,0.80_r_kind, &
0.78_r_kind,0.77_r_kind,0.77_r_kind,0.76_r_kind,0.72_r_kind,0.72_r_kind/
data (em(8,k),k=1,nch)/0.93_r_kind,0.94_r_kind,0.96_r_kind,0.96_r_kind, &
0.95_r_kind,0.93_r_kind,0.87_r_kind,0.86_r_kind,0.74_r_kind,0.65_r_kind/
data (em(9,k),k=1,nch)/0.87_r_kind,0.86_r_kind,0.84_r_kind,0.80_r_kind, &
0.76_r_kind,0.76_r_kind,0.75_r_kind,0.75_r_kind,0.70_r_kind,0.69_r_kind/
data (em(10,k),k=1,nch)/0.87_r_kind,0.86_r_kind,0.83_r_kind,0.77_r_kind, &
0.73_r_kind,0.68_r_kind,0.66_r_kind,0.66_r_kind,0.68_r_kind,0.67_r_kind/
data (em(11,k),k=1,nch)/0.89_r_kind,0.89_r_kind,0.88_r_kind,0.87_r_kind, &
0.86_r_kind,0.82_r_kind,0.77_r_kind,0.76_r_kind,0.69_r_kind,0.64_r_kind/
data (em(12,k),k=1,nch)/0.88_r_kind,0.87_r_kind,0.86_r_kind,0.83_r_kind, &
0.81_r_kind,0.77_r_kind,0.74_r_kind,0.73_r_kind,0.69_r_kind,0.64_r_kind/
data (em(13,k),k=1,nch)/0.86_r_kind,0.86_r_kind,0.86_r_kind,0.85_r_kind, &
0.82_r_kind,0.78_r_kind,0.69_r_kind,0.68_r_kind,0.51_r_kind,0.47_r_kind/
data (em(14,k),k=1,nch)/0.89_r_kind,0.88_r_kind,0.87_r_kind,0.83_r_kind, &
0.80_r_kind,0.75_r_kind,0.70_r_kind,0.70_r_kind,0.64_r_kind,0.60_r_kind/
data (em(15,k),k=1,nch)/0.91_r_kind,0.92_r_kind,0.93_r_kind,0.88_r_kind, &
0.84_r_kind,0.76_r_kind,0.66_r_kind,0.64_r_kind,0.48_r_kind,0.44_r_kind/
data (em(16,k),k=1,nch)/0.94_r_kind,0.95_r_kind,0.97_r_kind,0.91_r_kind, &
0.86_r_kind,0.74_r_kind,0.63_r_kind,0.63_r_kind,0.50_r_kind,0.45_r_kind/
data freq/4.9_r_kind,6.93_r_kind,10.65_r_kind,18.7_r_kind,23.8_r_kind, &
31.4_r_kind, 50.3_r_kind,52.5_r_kind, 89.0_r_kind, 150._r_kind/
! Initialization for emissivity at certain frequency
! In case of no any inputs available for various options
! A constant snow type & snow emissivity spectrum is assumed
! (e.g., powder) snow_type = 4
! Specify snow emissivity at required frequency
do ich = 2, nch
if(frequency < freq(1)) exit
if(frequency >= freq(nch)) exit
if(frequency < freq(ich)) then
emissivity = em(4,ich-1) + (em(4,ich) - em(4,ich-1)) &
*(frequency - freq(ich-1))/(freq(ich) - freq(ich-1))
exit
end if
end do
! Extrapolate to lower frequencies than 4.9GHz
if (frequency <= freq(1)) then
kratio = (em(4,2) - em(4,1))/(freq(2) - freq(1))
bconst = em(4,1) - kratio*freq(1)
emissivity = kratio*frequency + bconst
if(emissivity > one) emissivity = one
if(emissivity <= 0.8_r_kind) emissivity = 0.8_r_kind
end if
! Assume emissivity = constant at frequencies >= 150 GHz
if (frequency >= freq(nch)) emissivity = em(4,nch)
em_vector(1) = emissivity
em_vector(2) = emissivity
return
end subroutine em_initialization
!---------------------------------------------------------------------!
subroutine em_interpolate(frequency,discriminator,emissivity,snow_type)
!$$$ subprogram documentation block
!
! subprogram: determine snow_type and calculate emissivity
!
! prgmmr:Banghua Yan org: nesdis date: 2003-08-18
!
! abstract: 1. Find one snow emissivity spectrum to mimic the emission
! property of the realistic snow condition using a set of
! discrminators
! 2. Interpolate/extrapolate emissivity at a required frequency
!
! program history log:
!
! input argument list:
!
! frequency - frequency in GHz
! discriminators - emissivity discriminators at five AMSU-A & B window
! channels
! discriminator[1] : emissivity discriminator at 23.8 GHz
! discriminator[2] : emissivity discriminator at 31.4 GHz
! discriminator[3] : emissivity discriminator at 50.3 GHz
! discriminator[4] : emissivity discriminator at 89 GHz
! discriminator[5] : emissivity discriminator at 150 GHz
!
! Note: discriminator(1) and discriminator(3) are missing value in
! 'AMSU-B & Ts','AMUS-B' and 'MODL' options., which are defined to as -999.9,
! output argument list:
!
! em_vector[1] and [2] - emissivity at two polarizations.
! snow_type - snow type
!
! important internal variables:
!
! freq[1 ~ 10] - ten frequencies for sixteen snow types of emissivity
! em[1~16,*] - sixteen snow emissivity spectra
!
! remarks:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$
! use kinds, only: r_kind,i_kind
! use constants, only: zero, one
implicit none
integer(i_kind),parameter:: ncand = 16,nch =10
integer(i_kind):: ich,ichmin,ichmax,i,j,k,s,snow_type
real(r_kind) :: dem,demmin0
real(r_kind) :: em(ncand,nch)
real(r_kind) :: frequency,freq(nch),emissivity,discriminator(*),emis(nch)
real(r_kind) :: cor_factor,adjust_check,kratio, bconst
data freq/4.9_r_kind, 6.93_r_kind, 10.65_r_kind, 18.7_r_kind,&
23.8_r_kind, 31.4_r_kind, 50.3_r_kind, 52.5_r_kind, &
89.0_r_kind, 150._r_kind/
! Sixteen candidate snow emissivity spectra
data (em(1,k),k=1,nch)/0.87_r_kind,0.89_r_kind,0.91_r_kind,0.93_r_kind,0.94_r_kind,&
0.94_r_kind,0.94_r_kind,0.93_r_kind,0.92_r_kind,0.90_r_kind/
data (em(2,k),k=1,nch)/0.91_r_kind,0.91_r_kind,0.92_r_kind,0.91_r_kind,0.90_r_kind,&
0.90_r_kind,0.91_r_kind,0.91_r_kind,0.91_r_kind,0.86_r_kind/
data (em(3,k),k=1,nch)/0.90_r_kind,0.89_r_kind,0.88_r_kind,0.87_r_kind,0.86_r_kind,&
0.86_r_kind,0.85_r_kind,0.85_r_kind,0.82_r_kind,0.82_r_kind/
data (em(4,k),k=1,nch)/0.91_r_kind,0.91_r_kind,0.93_r_kind,0.93_r_kind,0.93_r_kind,&
0.93_r_kind,0.89_r_kind,0.88_r_kind,0.79_r_kind,0.79_r_kind/
data (em(5,k),k=1,nch)/0.90_r_kind,0.89_r_kind,0.88_r_kind,0.85_r_kind,0.84_r_kind,&
0.83_r_kind,0.83_r_kind,0.82_r_kind,0.79_r_kind,0.73_r_kind/
data (em(6,k),k=1,nch)/0.90_r_kind,0.89_r_kind,0.86_r_kind,0.82_r_kind,0.80_r_kind,&
0.79_r_kind,0.78_r_kind,0.78_r_kind,0.77_r_kind,0.77_r_kind/
data (em(7,k),k=1,nch)/0.88_r_kind,0.86_r_kind,0.85_r_kind,0.80_r_kind,0.78_r_kind,&
0.77_r_kind,0.77_r_kind,0.76_r_kind,0.72_r_kind,0.72_r_kind/
data (em(8,k),k=1,nch)/0.93_r_kind,0.94_r_kind,0.96_r_kind,0.96_r_kind,0.95_r_kind,&
0.93_r_kind,0.87_r_kind,0.86_r_kind,0.74_r_kind,0.65_r_kind/
data (em(9,k),k=1,nch)/0.87_r_kind,0.86_r_kind,0.84_r_kind,0.80_r_kind,0.76_r_kind,&
0.76_r_kind,0.75_r_kind,0.75_r_kind,0.70_r_kind,0.69_r_kind/
data (em(10,k),k=1,nch)/0.87_r_kind,0.86_r_kind,0.83_r_kind,0.77_r_kind,0.73_r_kind,&
0.68_r_kind,0.66_r_kind,0.66_r_kind,0.68_r_kind,0.67_r_kind/
data (em(11,k),k=1,nch)/0.89_r_kind,0.89_r_kind,0.88_r_kind,0.87_r_kind,0.86_r_kind,&
0.82_r_kind,0.77_r_kind,0.76_r_kind,0.69_r_kind,0.64_r_kind/
data (em(12,k),k=1,nch)/0.88_r_kind,0.87_r_kind,0.86_r_kind,0.83_r_kind,0.81_r_kind,&
0.77_r_kind,0.74_r_kind,0.73_r_kind,0.69_r_kind,0.64_r_kind/
data (em(13,k),k=1,nch)/0.86_r_kind,0.86_r_kind,0.86_r_kind,0.85_r_kind,0.82_r_kind,&
0.78_r_kind,0.69_r_kind,0.68_r_kind,0.51_r_kind,0.47_r_kind/
data (em(14,k),k=1,nch)/0.89_r_kind,0.88_r_kind,0.87_r_kind,0.83_r_kind,0.80_r_kind,&
0.75_r_kind,0.70_r_kind,0.70_r_kind,0.64_r_kind,0.60_r_kind/
data (em(15,k),k=1,nch)/0.91_r_kind,0.92_r_kind,0.93_r_kind,0.88_r_kind,0.84_r_kind,&
0.76_r_kind,0.66_r_kind,0.64_r_kind,0.48_r_kind,0.44_r_kind/
data (em(16,k),k=1,nch)/0.94_r_kind,0.95_r_kind,0.97_r_kind,0.91_r_kind,0.86_r_kind,&
0.74_r_kind,0.63_r_kind,0.63_r_kind,0.50_r_kind,0.45_r_kind/
save em
! Adjust unreasonable discriminator
if (discriminator(4) > discriminator(2)) &
discriminator(4) = discriminator(2) +(discriminator(5) - discriminator(2))* &
(150.0_r_kind - 89.0_r_kind)/(150.0_r_kind - 31.4_r_kind)
if ( (discriminator(3) /= -999.9_r_kind) .and. &
( ((discriminator(3)-0.01_r_kind) > discriminator(2)) .or. &
((discriminator(3)-0.01_r_kind) < discriminator(4))) ) &
discriminator(3) = discriminator(2) + &
(discriminator(4) - discriminator(2))*(89.0_r_kind - 50.3_r_kind) &
/ (89.0_r_kind - 31.4_r_kind)
! Find a snow emissivity spectrum
if(snow_type .eq. -999) then
demmin0 = 10.0_r_kind
do k = 1, ncand
dem = zero
ichmin = 1
ichmax = 3
if(discriminator(1) == -999.9_r_kind) then
ichmin = 2
ichmax = 2
end if
do ich = ichmin,ichmax
dem = dem + abs(discriminator(ich) - em(k,ich+4))
end do
do ich = 4,5
dem = dem + abs(discriminator(ich) - em(k,ich+5))
end do
if (dem < demmin0) then
demmin0 = dem
snow_type = k
end if
end do
end if
! Shift snow emissivity according to discriminator at 31.4 GHz
cor_factor = discriminator(2) - em(snow_type,6)
do ich = 1, nch
emis(ich) = em(snow_type,ich) + cor_factor
if(emis(ich) .gt. one) emis(ich) = one
if(emis(ich) .lt. 0.3_r_kind) emis(ich) = 0.3_r_kind
end do
! Emisivity data quality control
adjust_check = zero
do ich = 5, 9
if (ich .le. 7) then
if (discriminator(ich - 4) .ne. -999.9_r_kind) &
adjust_check = adjust_check + abs(emis(ich) - discriminator(ich - 4))
else
if (discriminator(ich - 4) .ne. -999.9_r_kind) &
adjust_check = adjust_check + abs(emis(ich+1) - discriminator(ich - 4))
end if
end do
if (adjust_check >= 0.04_r_kind) then
if (discriminator(1) /= -999.9_r_kind) then
if (discriminator(1) < emis(4)) then
emis(5) = emis(4) + &
(31.4_r_kind - 23.8_r_kind) * &
(discriminator(2) - emis(4))/(31.4_r_kind - 18.7_r_kind)
else
emis(5) = discriminator(1)
end if
end if
emis(6) = discriminator(2)
if (discriminator(3) /= -999.9_r_kind) then
emis(7) = discriminator(3)
else
! In case of missing the emissivity discriminator at 50.3 GHz
emis(7) = emis(6) + (89.0_r_kind - 50.3_r_kind) * &
(discriminator(4) - emis(6))/(89.0_r_kind - 31.4_r_kind)
end if
emis(8) = emis(7)
emis(9) = discriminator(4)
emis(10) = discriminator(5)
end if
! Estimate snow emissivity at a required frequency
do i = 2, nch
if(frequency < freq(1)) exit
if(frequency >= freq(nch)) exit
if(frequency < freq(i)) then
emissivity = emis(i-1) + (emis(i) - emis(i-1))*(frequency - freq(i-1)) &
/(freq(i) - freq(i-1))
exit
end if
end do
! Extrapolate to lower frequencies than 4.9GHz
if (frequency <= freq(1)) then
kratio = (emis(2) - emis(1))/(freq(2) - freq(1))
bconst = emis(1) - kratio*freq(1)
emissivity = kratio*frequency + bconst
if(emissivity > one) emissivity = one
if(emissivity <= 0.8_r_kind) emissivity = 0.8_r_kind
end if
! Assume emissivity = constant at frequencies >= 150 GHz
if (frequency >= freq(nch)) emissivity = emis(nch)
return
end subroutine em_interpolate
!---------------------------------------------------------------------!
subroutine sem_ABTs(theta,frequency,tb,ts,snow_type,em_vector)
!$$$ subprogram documentation block
!
! subprogram:
!
! prgmmr:Banghua Yan org: nesdis date: 2003-08-18
!
! abstract:
! Calculate the emissivity discriminators and interpolate/extrapolate
! emissivity at a required frequency with respect to secenery ABTs
!
! program history log:
!
! input argument list:
!
! frequency - frequency in GHz
! theta - local zenith angle (not used here)
! tb[1] ~ tb[5] - brightness temperature at five AMSU window channels:
! tb[1] : 23.8 GHz
! tb[2] : 31.4 GHz
! tb[3] : 50.3 GHz
! tb[4] : 89.0 GHz
! tb[5] : 150 GHz
!
! output argument list:
!
! em_vector[1] and [2] - emissivity at two polarizations.
! set esv = esh here and will be updated
! snow_type - snow type
!
! important internal variables:
!
! nind - number of threshold in decision trees
! to identify each snow type ( = 6)
! em(1~16,*) - sixteen snow emissivity spectra
! DI_coe - coefficients to generate six discriminators to describe
! the overall emissivity variability within a wider frequency range
! threshold - thresholds in decision trees to identify snow types
! index_in - six indices to discriminate snow type
!
! remarks:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$
! use kinds, only: r_kind,i_kind
implicit none
integer(i_kind),parameter:: ncand = 16,nch =10,nthresh=38
integer(i_kind),parameter:: nind=6,ncoe=8,nLIcoe=6,nHIcoe=12
integer(i_kind):: ich,i,j,k,num,npass,snow_type,md0,md1,nmodel(ncand-1)
real(r_kind) :: theta,frequency,tb150,LI,HI,DS1,DS2,DS3
real(r_kind) :: em(ncand,nch), em_vector(*)
real(r_kind) :: tb(*),freq(nch),DTB(nind-1),DI(nind-1), &
DI_coe(nind-1,0:ncoe-1),threshold(nthresh,nind), &
index_in(nind),threshold0(nind)
real(r_kind) :: LI_coe(0:nLIcoe-1),HI_coe(0:nHIcoe-1)
real(r_kind) :: ts,emissivity
real(r_kind) :: discriminator(5)
logical:: pick_status,tindex(nind)
save em,threshold,DI_coe,LI_coe, HI_coe,nmodel,freq
data freq/4.9_r_kind,6.93_r_kind,10.65_r_kind,18.7_r_kind,23.8_r_kind, &
31.4_r_kind, 50.3_r_kind,52.5_r_kind, 89.0_r_kind, 150._r_kind/
data nmodel/5,10,13,16,18,24,30,31,32,33,34,35,36,37,38/
! Fitting coefficients for five discriminators
data (DI_coe(1,k),k=0,ncoe-1)/ &
3.285557e-002_r_kind, 2.677179e-005_r_kind, &
4.553101e-003_r_kind, 5.639352e-005_r_kind, &
-1.825188e-004_r_kind, 1.636145e-004_r_kind, &
1.680881e-005_r_kind, -1.708405e-004_r_kind/
data (DI_coe(2,k),k=0,ncoe-1)/ &
-4.275539e-002_r_kind, -2.541453e-005_r_kind, &
4.154796e-004_r_kind, 1.703443e-004_r_kind, &
4.350142e-003_r_kind, 2.452873e-004_r_kind, &
-4.748506e-003_r_kind, 2.293836e-004_r_kind/
data (DI_coe(3,k),k=0,ncoe-1)/ &
-1.870173e-001_r_kind, -1.061678e-004_r_kind, &
2.364055e-004_r_kind, -2.834876e-005_r_kind, &
4.899651e-003_r_kind, -3.418847e-004_r_kind, &
-2.312224e-004_r_kind, 9.498600e-004_r_kind/
data (DI_coe(4,k),k=0,ncoe-1)/ &
-2.076519e-001_r_kind, 8.475901e-004_r_kind, &
-2.072679e-003_r_kind, -2.064717e-003_r_kind, &
2.600452e-003_r_kind, 2.503923e-003_r_kind, &
5.179711e-004_r_kind, 4.667157e-005_r_kind/
data (DI_coe(5,k),k=0,ncoe-1)/ &
-1.442609e-001_r_kind, -8.075003e-005_r_kind, &
-1.790933e-004_r_kind, -1.986887e-004_r_kind, &
5.495115e-004_r_kind, -5.871732e-004_r_kind, &
4.517280e-003_r_kind, 7.204695e-004_r_kind/
! Fitting coefficients for emissivity index at 31.4 GHz
data LI_coe/ &
7.963632e-001_r_kind, 7.215580e-003_r_kind, &
-2.015921e-005_r_kind, -1.508286e-003_r_kind, &
1.731405e-005_r_kind, -4.105358e-003_r_kind/
! Fitting coefficients for emissivity index at 150 GHz
data HI_coe/ &
1.012160e+000_r_kind, 6.100397e-003_r_kind, &
-1.774347e-005_r_kind, -4.028211e-003_r_kind, &
1.224470e-005_r_kind, 2.345612e-003_r_kind, &
-5.376814e-006_r_kind, -2.795332e-003_r_kind, &
8.072756e-006_r_kind, 3.529615e-003_r_kind, &
1.955293e-006_r_kind, -4.942230e-003_r_kind/
! Six thresholds for sixteen candidate snow types
! Note: some snow type contains several possible
! selections for six thresholds
!1 Wet Snow
data (threshold(1,k),k=1,6)/0.88_r_kind,0.86_r_kind,-999.9_r_kind,&
0.01_r_kind,0.01_r_kind,200._r_kind/
data (threshold(2,k),k=1,6)/0.88_r_kind,0.85_r_kind,-999.9_r_kind,&
0.06_r_kind,0.10_r_kind,200._r_kind/
data (threshold(3,k),k=1,6)/0.88_r_kind,0.83_r_kind,-0.02_r_kind,&
0.12_r_kind,0.16_r_kind,204._r_kind/
data (threshold(4,k),k=1,6)/0.90_r_kind,0.89_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
data (threshold(5,k),k=1,6)/0.92_r_kind,0.85_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
!2 Grass_after_Snow
data (threshold(6,k),k=1,6)/0.84_r_kind,0.83_r_kind,-999.9_r_kind,&
0.08_r_kind,0.10_r_kind,195._r_kind/
data (threshold(7,k),k=1,6)/0.85_r_kind,0.85_r_kind,-999.9_r_kind,&
0.10_r_kind,-999.9_r_kind,190._r_kind/
data (threshold(8,k),k=1,6)/0.86_r_kind,0.81_r_kind,-999.9_r_kind,&
0.12_r_kind,-999.9_r_kind,200._r_kind/
data (threshold(9,k),k=1,6)/0.86_r_kind,0.81_r_kind,0.0_r_kind,&
0.12_r_kind,-999.9_r_kind,189._r_kind/
data (threshold(10,k),k=1,6)/0.90_r_kind,0.81_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,195._r_kind/
!3 RS_Snow (A)
data (threshold(11,k),k=1,6)/0.80_r_kind,0.76_r_kind,-999.9_r_kind,&
0.05_r_kind,-999.9_r_kind,185._r_kind/
data (threshold(12,k),k=1,6)/0.82_r_kind,0.78_r_kind,-999.9_r_kind,&
-999.9_r_kind,0.25_r_kind,180._r_kind/
data (threshold(13,k),k=1,6)/0.90_r_kind,0.76_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,180._r_kind/
!4 Powder Snow
data (threshold(14,k),k=1,6)/0.89_r_kind,0.73_r_kind,-999.9_r_kind,&
0.20_r_kind,-999.9_r_kind,-999.9_r_kind/
data (threshold(15,k),k=1,6)/0.89_r_kind,0.75_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
data (threshold(16,k),k=1,6)/0.93_r_kind,0.72_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
!5 RS_Snow (B)
data (threshold(17,k),k=1,6)/0.82_r_kind,0.70_r_kind,-999.9_r_kind,&
0.20_r_kind,-999.9_r_kind,160._r_kind/
data (threshold(18,k),k=1,6)/0.83_r_kind,0.70_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,160._r_kind/
!6 RS_Snow (C)
data (threshold(19,k),k=1,6)/0.75_r_kind,0.76_r_kind,-999.9_r_kind,&
0.08_r_kind,-999.9_r_kind,172._r_kind/
data (threshold(20,k),k=1,6)/0.77_r_kind,0.72_r_kind,-999.9_r_kind,&
0.12_r_kind,0.15_r_kind,175._r_kind/
data (threshold(21,k),k=1,6)/0.78_r_kind,0.74_r_kind,-999.9_r_kind,&
-999.9_r_kind,0.20_r_kind,172._r_kind/
data (threshold(22,k),k=1,6)/0.80_r_kind,0.77_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,170._r_kind/
data (threshold(23,k),k=1,6)/0.82_r_kind,-999.9_r_kind,-999.9_r_kind,&
0.15_r_kind,0.22_r_kind,170._r_kind/
data (threshold(24,k),k=1,6)/0.82_r_kind,0.73_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,170._r_kind/
!7 RS_Snow (D)
data (threshold(25,k),k=1,6)/0.75_r_kind,0.70_r_kind,-999.9_r_kind,&
0.15_r_kind,0.25_r_kind,167._r_kind/
data (threshold(26,k),k=1,6)/0.77_r_kind,0.76_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
data (threshold(27,k),k=1,6)/0.80_r_kind,0.72_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
data (threshold(28,k),k=1,6)/0.77_r_kind,0.73_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
data (threshold(29,k),k=1,6)/0.81_r_kind,0.71_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
data (threshold(30,k),k=1,6)/0.82_r_kind,0.69_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
!8 Thin Crust Snow
data (threshold(31,k),k=1,6)/0.88_r_kind,0.58_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
!9 RS_Snow (E)
data (threshold(32,k),k=1,6)/0.73_r_kind,0.67_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
!10 Bottom Crust Snow (A)
data (threshold(33,k),k=1,6)/0.83_r_kind,0.66_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
!11 Shallow Snow
data (threshold(34,k),k=1,6)/0.82_r_kind,0.60_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
!12 Deep Snow
data (threshold(35,k),k=1,6)/0.77_r_kind,0.60_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
!13 Crust Snow
data (threshold(36,k),k=1,6)/0.77_r_kind,0.7_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
!14 Medium Snow
data (threshold(37,k),k=1,6)/-999.9_r_kind,0.55_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
!15 Bottom Crust Snow(B)
data (threshold(38,k),k=1,6)/0.74_r_kind,-999.9_r_kind,-999.9_r_kind,&
-999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
!16 Thick Crust Snow
! lowest priority: No constraints
! Sixteen candidate snow emissivity spectra
data (em(1,k),k=1,nch)/0.87_r_kind,0.89_r_kind,0.91_r_kind,0.93_r_kind,&
0.94_r_kind,0.94_r_kind,0.94_r_kind,0.93_r_kind,0.92_r_kind,0.90_r_kind/
data (em(2,k),k=1,nch)/0.91_r_kind,0.91_r_kind,0.92_r_kind,0.91_r_kind,&
0.90_r_kind,0.90_r_kind,0.91_r_kind,0.91_r_kind,0.91_r_kind,0.86_r_kind/
data (em(3,k),k=1,nch)/0.90_r_kind,0.89_r_kind,0.88_r_kind,0.87_r_kind,&
0.86_r_kind,0.86_r_kind,0.85_r_kind,0.85_r_kind,0.82_r_kind,0.82_r_kind/
data (em(4,k),k=1,nch)/0.91_r_kind,0.91_r_kind,0.93_r_kind,0.93_r_kind,&
0.93_r_kind,0.93_r_kind,0.89_r_kind,0.88_r_kind,0.79_r_kind,0.79_r_kind/
data (em(5,k),k=1,nch)/0.90_r_kind,0.89_r_kind,0.88_r_kind,0.85_r_kind,&
0.84_r_kind,0.83_r_kind,0.83_r_kind,0.82_r_kind,0.79_r_kind,0.73_r_kind/
data (em(6,k),k=1,nch)/0.90_r_kind,0.89_r_kind,0.86_r_kind,0.82_r_kind,&
0.80_r_kind,0.79_r_kind,0.78_r_kind,0.78_r_kind,0.77_r_kind,0.77_r_kind/
data (em(7,k),k=1,nch)/0.88_r_kind,0.86_r_kind,0.85_r_kind,0.80_r_kind,&
0.78_r_kind,0.77_r_kind,0.77_r_kind,0.76_r_kind,0.72_r_kind,0.72_r_kind/
data (em(8,k),k=1,nch)/0.93_r_kind,0.94_r_kind,0.96_r_kind,0.96_r_kind,&
0.95_r_kind,0.93_r_kind,0.87_r_kind,0.86_r_kind,0.74_r_kind,0.65_r_kind/
data (em(9,k),k=1,nch)/0.87_r_kind,0.86_r_kind,0.84_r_kind,0.80_r_kind,&
0.76_r_kind,0.76_r_kind,0.75_r_kind,0.75_r_kind,0.70_r_kind,0.69_r_kind/
data (em(10,k),k=1,nch)/0.87_r_kind,0.86_r_kind,0.83_r_kind,0.77_r_kind,&
.73_r_kind,0.68_r_kind,0.66_r_kind,0.66_r_kind,0.68_r_kind,0.67_r_kind/
data (em(11,k),k=1,nch)/0.89_r_kind,0.89_r_kind,0.88_r_kind,0.87_r_kind,&
0.86_r_kind,0.82_r_kind,0.77_r_kind,0.76_r_kind,0.69_r_kind,0.64_r_kind/
data (em(12,k),k=1,nch)/0.88_r_kind,0.87_r_kind,0.86_r_kind,0.83_r_kind,&
0.81_r_kind,0.77_r_kind,0.74_r_kind,0.73_r_kind,0.69_r_kind,0.64_r_kind/
data (em(13,k),k=1,nch)/0.86_r_kind,0.86_r_kind,0.86_r_kind,0.85_r_kind,&
0.82_r_kind,0.78_r_kind,0.69_r_kind,0.68_r_kind,0.51_r_kind,0.47_r_kind/
data (em(14,k),k=1,nch)/0.89_r_kind,0.88_r_kind,0.87_r_kind,0.83_r_kind,&
0.80_r_kind,0.75_r_kind,0.70_r_kind,0.70_r_kind,0.64_r_kind,0.60_r_kind/
data (em(15,k),k=1,nch)/0.91_r_kind,0.92_r_kind,0.93_r_kind,0.88_r_kind,&
0.84_r_kind,0.76_r_kind,0.66_r_kind,0.64_r_kind,0.48_r_kind,0.44_r_kind/
data (em(16,k),k=1,nch)/0.94_r_kind,0.95_r_kind,0.97_r_kind,0.91_r_kind,&
0.86_r_kind,0.74_r_kind,0.63_r_kind,0.63_r_kind,0.50_r_kind,0.45_r_kind/
!*** DEFINE SIX DISCRIMINATORS
dtb(1) = tb(1) - tb(2)
dtb(2) = tb(2) - tb(4)
dtb(3) = tb(2) - tb(5)
dtb(4) = tb(3) - tb(5)
dtb(5) = tb(4) - tb(5)
tb150 = tb(5)
LI = LI_coe(0)
do i=0,1
LI = LI + LI_coe(2*i+1)*tb(i+1) + LI_coe(2*i+2)*tb(i+1)*tb(i+1)
end do
LI = LI + LI_coe(nLIcoe-1)*ts
HI = HI_coe(0)
do i=0,4
HI = HI + HI_coe(2*i+1)*tb(i+1) + HI_coe(2*i+2)*tb(i+1)*tb(i+1)
end do
HI = HI + HI_coe(nHIcoe-1)*ts
do num=1,nind-1
DI(num) = DI_coe(num,0) + DI_coe(num,1)*tb(2)
do i=1,5
DI(num) = DI(num) + DI_coe(num,1+i)*DTB(i)
end do
DI(num) = DI(num) + DI_coe(num,ncoe-1)*ts
end do
!*** DEFINE FIVE INDIES
!HI = DI(0) - DI(3)
DS1 = DI(1) + DI(2)
DS2 = DI(4) + DI(5)
DS3 = DS1 + DS2 + DI(3)
index_in(1) = LI
index_in(2) = HI
index_in(3) = DS1
index_in(4) = DS2
index_in(5) = DS3
index_in(6) = tb150
!*** IDENTIFY SNOW TYPE
! Initialization
md0 = 1
snow_type = ncand
pick_status = .false.
! Pick one snow type
! Check all possible selections for six thresholds for each snow type
do i = 1, ncand - 1
md1 = nmodel(i)
do j = md0, md1
npass = 0
do k = 1 , nind
threshold0(k) = threshold(j,k)
end do
CALL six_indices(nind,index_in,threshold0,tindex)
! Corrections
if((i == 5) .and. (index_in(2) > 0.75_r_kind)) tindex(2) = .false.
if((i == 5) .and. (index_in(4) > 0.20_r_kind) &
.and. (index_in(1) > 0.88_r_kind)) tindex(1) = .false.
if((i == 10) .and. (index_in(1) <= 0.83_r_kind)) tindex(1) = .true.
if((i == 13) .and. (index_in(2) < 0.52_r_kind)) tindex(2) = .true.
do k = 1, nind
if(.not.tindex(k)) exit
npass = npass + 1
end do
if(npass == nind) exit
end do
if(npass == nind) then
pick_status = .true.
snow_type = i
end if
if(pick_status) exit
md0 = md1 + 1
end do
discriminator(1) = LI + DI(1)
discriminator(2) = LI
discriminator(3) = DI(4) + HI
discriminator(4) = LI - DI(2)
discriminator(5) = HI
call em_interpolate(frequency,discriminator,emissivity,snow_type)
em_vector(1) = emissivity
em_vector(2) = emissivity
return
end subroutine sem_ABTs
!---------------------------------------------------------------------!
subroutine six_indices(nind,index_in,threshold,tindex)
!$$$ subprogram documentation block
!
! subprogram:
!
! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
!
! abstract:
!
! program history log:
!
! input argument list:
!
! nind - Number of threshold in decision trees
! to identify each snow type ( = 6)
! index_in - six indices to discriminate snow type
! threshold - Thresholds in decision trees to identify snow types
!
! output argument list:
!
! tindex - state vaiable to show surface snow emissivity feature
! tindex[ ] = .T.: snow emissivity feature matches the
! corresponding threshold for certain snow type
! tindex[ ] = .F.: snow emissivity feature doesn't match the
! corresponding threshold for certain snow type
!
! remarks:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$
! use kinds, only: r_kind,i_kind
implicit none
integer(i_kind) :: i,nind
real(r_kind) :: index_in(*),threshold(*)
logical :: tindex(*)
do i=1,nind
tindex(i) = .false.
if (threshold(i) .eq. -999.9_r_kind) then
tindex(i) = .true.
else
if ( (i .le. 2) .or. (i .gt. (nind-1)) ) then
if (index_in(i) .ge. threshold(i)) tindex(i) = .true.
else
if (index_in(i) .le. threshold(i)) tindex(i) = .true.
end if
end if
end do
return
end subroutine six_indices
!---------------------------------------------------------------------!
subroutine sem_AB(theta,frequency,tb,snow_type,em_vector)
!$$$ subprogram documentation block
!
! subprogram:
!
! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
!
! abstract:
! Calculate the emissivity discriminators and interpolate/extrapolate
! emissivity at required frequency with respect to option AMSUAB
!
! program history log:
! 2004-10-28 treadon - correct problem in declared dimensions of array coe
!
! input argument list:
!
! frequency - frequency in GHz
! theta - local zenith angle (not used here)
! tb[1]~tb[5] - brightness temperature at five AMSU-A & B window channels:
! tb[1] : 23.8 GHz
! tb[2] : 31.4 GHz
! tb[3] : 50.3 GHz
! tb[4] : 89 GHz
! tb[5] : 150 GHz
!
! output argument list:
!
! em_vector[1] and [2] - emissivity at two polarizations.
! set esv = esh here and will be updated
! snow_type - snow type (reference [2])
!
! important internal variables:
!
! coe23 - fitting coefficients to estimate discriminator at 23.8 GHz
! coe31 - fitting coefficients to estimate discriminator at 31.4 GHz
! coe50 - fitting coefficients to estimate discriminator at 50.3 GHz
! coe89 - fitting coefficients to estimate discriminator at 89 GHz
! coe150 - fitting coefficients to estimate discriminator at 150 GHz
!
! remarks:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$
! use kinds, only: r_kind,i_kind
implicit none
integer(i_kind),parameter:: nch =10,nwch = 5,ncoe = 10
real(r_kind) :: tb(*),theta,frequency
real(r_kind) :: em_vector(*),emissivity,discriminator(nwch)
integer(i_kind) :: i,snow_type,k,ich,nvalid_ch
real(r_kind) :: coe23(0:ncoe),coe31(0:ncoe),coe50(0:ncoe),coe89(0:ncoe),coe150(0:ncoe)
real(r_kind) :: coe(nch*(ncoe+1))
Equivalence (coe(1),coe23)
Equivalence (coe(12),coe31)
Equivalence (coe(23),coe50)
Equivalence (coe(34),coe89)
Equivalence (coe(45),coe150)
! Fitting Coefficients at 23.8 GHz: Using Tb1 ~ Tb3
data (coe23(k),k=0,6)/&
-1.326040e+000_r_kind, 2.475904e-002_r_kind, &
-5.741361e-005_r_kind, -1.889650e-002_r_kind, &
6.177911e-005_r_kind, 1.451121e-002_r_kind, &
-4.925512e-005_r_kind/
! Fitting Coefficients at 31.4 GHz: Using Tb1 ~ Tb3
data (coe31(k),k=0,6)/ &
-1.250541e+000_r_kind, 1.911161e-002_r_kind, &
-5.460238e-005_r_kind, -1.266388e-002_r_kind, &
5.745064e-005_r_kind, 1.313985e-002_r_kind, &
-4.574811e-005_r_kind/
! Fitting Coefficients at 50.3 GHz: Using Tb1 ~ Tb3
data (coe50(k),k=0,6)/ &
-1.246754e+000_r_kind, 2.368658e-002_r_kind, &
-8.061774e-005_r_kind, -3.206323e-002_r_kind, &
1.148107e-004_r_kind, 2.688353e-002_r_kind, &
-7.358356e-005_r_kind/
! Fitting Coefficients at 89 GHz: Using Tb1 ~ Tb4
data (coe89(k),k=0,8)/ &
-1.278780e+000_r_kind, 1.625141e-002_r_kind, &
-4.764536e-005_r_kind, -1.475181e-002_r_kind, &
5.107766e-005_r_kind, 1.083021e-002_r_kind, &
-4.154825e-005_r_kind, 7.703879e-003_r_kind, &
-6.351148e-006_r_kind/
! Fitting Coefficients at 150 GHz: Using Tb1 ~ Tb5
data coe150/&
-1.691077e+000_r_kind, 3.352403e-002_r_kind, &
-7.310338e-005_r_kind, -4.396138e-002_r_kind, &
1.028994e-004_r_kind, 2.301014e-002_r_kind, &
-7.070810e-005_r_kind, 1.270231e-002_r_kind, &
-2.139023e-005_r_kind, -2.257991e-003_r_kind, &
1.269419e-005_r_kind/
save coe23,coe31,coe50,coe89,coe150
! Calculate emissivity discriminators at five AMSU window channels
do ich = 1, nwch
discriminator(ich) = coe(1+(ich-1)*11)
if (ich .le. 3) nvalid_ch = 3
if (ich .eq. 4) nvalid_ch = 4
if (ich .eq. 5) nvalid_ch = 5
do i=1,nvalid_ch
discriminator(ich) = discriminator(ich) + coe((ich-1)*11 + 2*i)*tb(i) + &
coe((ich-1)*11 + 2*i+1)*tb(i)*tb(i)
end do
end do
! Identify one snow emissivity spectrum and interpolate/extrapolate emissivity
! at a required frequency
call em_interpolate(frequency,discriminator,emissivity,snow_type)
em_vector(1) = emissivity
em_vector(2) = emissivity
return
end subroutine sem_AB
!---------------------------------------------------------------------!
subroutine sem_ATs(theta,frequency,tba,ts,snow_type,em_vector)
!$$$ subprogram documentation block
!
! subprogram:
!
! prgmmr:Banghua Yan org: nesdis date: 2003-08-18
!
! abstract:
! Calculate the emissivity discriminators and interpolate/extrapolate
! emissivity at required frequency with respect to secenery AMSUAB
!
! program history log:
!
! input argument list:
!
! frequency - frequency in GHz
! theta - local zenith angle (not used here)
! ts - surface temperature
! tba[1] ~ tba[4] - brightness temperature at five AMSU-A window channels:
! tba[1] : 23.8 GHz
! tba[2] : 31.4 GHz
! tba[3] : 50.3 GHz
! tba[4] : 89 GHz
! output argument list:
!
! em_vector[1] and [2] - emissivity at two polarizations.
! set esv = esh here and will be updated
! snow_type - snow type (reference [2])
!
! important internal variables:
!
! coe23 - fitting coefficients to estimate discriminator at 23.8 GHz
! coe31 - fitting coefficients to estimate discriminator at 31.4 GHz
! coe50 - fitting coefficients to estimate discriminator at 50.3 GHz
! coe89 - fitting coefficients to estimate discriminator at 89 GHz
! coe150 - fitting coefficients to estimate discriminator at 150 GHz
!
! remarks:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$
! use kinds, only: r_kind,i_kind
implicit none
integer(i_kind),parameter:: nch =10,nwch = 5,ncoe = 9
real(r_kind) :: tba(*),theta
real(r_kind) :: em_vector(*),emissivity,ts,frequency,discriminator(nwch)
integer(i_kind) :: snow_type,i,k,ich,nvalid_ch
real(r_kind) :: coe23(0:ncoe),coe31(0:ncoe),coe50(0:ncoe),coe89(0:ncoe),coe150(0:ncoe)
real(r_kind) :: coe(nch*(ncoe+1))
Equivalence (coe(1),coe23)
Equivalence (coe(11),coe31)
Equivalence (coe(21),coe50)
Equivalence (coe(31),coe89)
Equivalence (coe(41),coe150)
! Fitting Coefficients at 23.8 GHz: Using Tb1, Tb2 and Ts
data (coe23(k),k=0,5)/ &
8.210105e-001_r_kind, 1.216432e-002_r_kind, &
-2.113875e-005_r_kind, -6.416648e-003_r_kind, &
1.809047e-005_r_kind, -4.206605e-003_r_kind/
! Fitting Coefficients at 31.4 GHz: Using Tb1, Tb2 and Ts
data (coe31(k),k=0,5)/ &
7.963632e-001_r_kind, 7.215580e-003_r_kind, &
-2.015921e-005_r_kind, -1.508286e-003_r_kind, &
1.731405e-005_r_kind, -4.105358e-003_r_kind/
! Fitting Coefficients at 50.3 GHz: Using Tb1, Tb2, Tb3 and Ts
data (coe50(k),k=0,7)/ &
1.724160e+000_r_kind, 5.556665e-003_r_kind, &
-2.915872e-005_r_kind, -1.146713e-002_r_kind, &
4.724243e-005_r_kind, 3.851791e-003_r_kind, &
-5.581535e-008_r_kind, -5.413451e-003_r_kind/
! Fitting Coefficients at 89 GHz: Using Tb1 ~ Tb4 and Ts
data coe89/ &
9.962065e-001_r_kind, 1.584161e-004_r_kind, &
-3.988934e-006_r_kind, 3.427638e-003_r_kind, &
-5.084836e-006_r_kind, -6.178904e-004_r_kind, &
1.115315e-006_r_kind, 9.440962e-004_r_kind, &
9.711384e-006_r_kind, -4.259102e-003_r_kind/
! Fitting Coefficients at 150 GHz: Using Tb1 ~ Tb4 and Ts
data coe150/ &
-5.244422e-002_r_kind, 2.025879e-002_r_kind, &
-3.739231e-005_r_kind, -2.922355e-002_r_kind, &
5.810726e-005_r_kind, 1.376275e-002_r_kind, &
-3.757061e-005_r_kind, 6.434187e-003_r_kind, &
6.190403e-007_r_kind, -2.944785e-003_r_kind/
save coe23,coe31,coe50,coe89,coe150
! Calculate emissivity discriminators at five AMSU window channels
DO ich = 1, nwch
discriminator(ich) = coe(1+(ich-1)*10)
if (ich .le. 2) nvalid_ch = 2
if (ich .eq. 3) nvalid_ch = 3
if (ich .ge. 4) nvalid_ch = 4
do i=1,nvalid_ch
discriminator(ich) = discriminator(ich) + coe((ich-1)*10 + 2*i)*tba(i) + &
coe((ich-1)*10 + 2*i+1)*tba(i)*tba(i)
end do
discriminator(ich) = discriminator(ich) + coe( (ich-1)*10 + (nvalid_ch+1)*2 )*ts
end do
call em_interpolate(frequency,discriminator,emissivity,snow_type)
em_vector(1) = emissivity
em_vector(2) = emissivity
return
end subroutine sem_ATs
!---------------------------------------------------------------------!
subroutine sem_amsua(theta,frequency,tba,snow_type,em_vector)
!$$$ subprogram documentation block
!
! subprogram:
!
! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
!
! abstract:
! Calculate the emissivity discriminators and interpolate/extrapolate
! emissivity at required frequency with respect to secenery AMSUA
!
! program history log:
! 2004-10-28 treadon - correct problem in declared dimensions of array coe
!
! input argument list:
!
! frequency - frequency in GHz
! theta - local zenith angle (not used here)
! tba[1]~tba[4] - brightness temperature at five AMSU-A window channels:
! tba[1] : 23.8 GHz
! tba[2] : 31.4 GHz
! tba[3] : 50.3 GHz
! tba[4] : 89 GHz
!
! output argument list:
!
! em_vector(1) and (2) - emissivity at two polarizations.
! set esv = esh here and will be updated
! snow_type - snow type
!
! important internal variables:
!
! coe23 - fitting coefficients to estimate discriminator at 23.8 GHz
! coe31 - fitting coefficients to estimate discriminator at 31.4 GHz
! coe50 - fitting coefficients to estimate discriminator at 50.3 GHz
! coe89 - fitting coefficients to estimate discriminator at 89 GHz
! coe150 - fitting coefficients to estimate discriminator at 150 GHz
!
! remarks:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$
! use kinds, only: r_kind,i_kind
implicit none
integer(i_kind),parameter:: nch =10,nwch = 5,ncoe = 8
real(r_kind) :: tba(*),theta
real(r_kind) :: em_vector(*),emissivity,frequency,discriminator(nwch)
integer(i_kind) :: snow_type,i,k,ich,nvalid_ch
real(r_kind) :: coe23(0:ncoe),coe31(0:ncoe),coe50(0:ncoe),coe89(0:ncoe),coe150(0:ncoe)
real(r_kind) :: coe(nch*(ncoe+1))
Equivalence (coe(1),coe23)
Equivalence (coe(11),coe31)
Equivalence (coe(21),coe50)
Equivalence (coe(31),coe89)
Equivalence (coe(41),coe150)
! Fitting Coefficients at 23.8 GHz: Using Tb1 ~ Tb3
data (coe23(k),k=0,6)/ &
-1.326040e+000_r_kind, 2.475904e-002_r_kind, -5.741361e-005_r_kind, &
-1.889650e-002_r_kind, 6.177911e-005_r_kind, 1.451121e-002_r_kind, &
-4.925512e-005_r_kind/
! Fitting Coefficients at 31.4 GHz: Using Tb1 ~ Tb3
data (coe31(k),k=0,6)/ &
-1.250541e+000_r_kind, 1.911161e-002_r_kind, -5.460238e-005_r_kind, &
-1.266388e-002_r_kind, 5.745064e-005_r_kind, 1.313985e-002_r_kind, &
-4.574811e-005_r_kind/
! Fitting Coefficients at 50.3 GHz: Using Tb1 ~ Tb3
data (coe50(k),k=0,6)/ &
-1.246754e+000_r_kind, 2.368658e-002_r_kind, -8.061774e-005_r_kind, &
-3.206323e-002_r_kind, 1.148107e-004_r_kind, 2.688353e-002_r_kind, &
-7.358356e-005_r_kind/
! Fitting Coefficients at 89 GHz: Using Tb1 ~ Tb4
data coe89/ &
-1.278780e+000_r_kind, 1.625141e-002_r_kind, -4.764536e-005_r_kind, &
-1.475181e-002_r_kind, 5.107766e-005_r_kind, 1.083021e-002_r_kind, &
-4.154825e-005_r_kind, 7.703879e-003_r_kind, -6.351148e-006_r_kind/
! Fitting Coefficients at 150 GHz: Using Tb1 ~ Tb4
data coe150/ &
-1.624857e+000_r_kind, 3.138243e-002_r_kind, -6.757028e-005_r_kind, &
-4.178496e-002_r_kind, 9.691893e-005_r_kind, 2.165964e-002_r_kind, &
-6.702349e-005_r_kind, 1.111658e-002_r_kind, -1.050708e-005_r_kind/
save coe23,coe31,coe50,coe150
! Calculate emissivity discriminators at five AMSU window channels
do ich = 1, nwch
discriminator(ich) = coe(1+(ich-1)*10)
if (ich .le. 2) nvalid_ch = 3
if (ich .ge. 3) nvalid_ch = 4
do i=1,nvalid_ch
discriminator(ich) = discriminator(ich) + coe((ich-1)*10 + 2*i)*tba(i) + &
coe((ich-1)*10 + 2*i+1)*tba(i)*tba(i)
end do
end do
! Quality Control
if(discriminator(4) .gt. discriminator(2)) &
discriminator(4) = discriminator(2) + (150.0_r_kind - 89.0_r_kind)* &
(discriminator(5) - discriminator(2))/ &
(150.0_r_kind - 31.4_r_kind)
! Quality control at 50.3 GHz
if((discriminator(3) .gt. discriminator(2)) .or. &
(discriminator(3) .lt. discriminator(4))) &
discriminator(3) = discriminator(2) + (89.0_r_kind - 50.3_r_kind)* &
(discriminator(4) - discriminator(2))/(89.0_r_kind - 31.4_r_kind)
call em_interpolate(frequency,discriminator,emissivity,snow_type)
em_vector(1) = emissivity
em_vector(2) = emissivity
return
end subroutine sem_amsua
!---------------------------------------------------------------------!
subroutine sem_BTs(theta,frequency,tbb,ts,snow_type,em_vector)
!$$$ subprogram documentation block
!
! subprogram:
!
! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
!
! abstract:
! Calculate the emissivity discriminators and interpolate/extrapolate
! emissivity at required frequency with respect to secenery BTs
!
! program history log:
!
! input argument list:
!
! frequency - frequency in GHz
! theta - local zenith angle (not used here)
! ts - surface temperature in degree
! tbb[1] ~ tbb[2] - brightness temperature at five AMSU-B window channels:
! tbb[1] : 89 GHz
! tbb[2] : 150 GHz
!
! output argument list:
!
! em_vector(1) and (2) - emissivity at two polarizations.
! set esv = esh here and will be updated
! snow_type - snow type
!
! important internal variables:
!
! coe31 - fitting coefficients to estimate discriminator at 31.4 GHz
! coe89 - fitting coefficients to estimate discriminator at 89 GHz
! coe150 - fitting coefficients to estimate discriminator at 150 GHz
!
! remarks:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$
! use kinds, only: r_kind,i_kind
implicit none
integer(i_kind),parameter:: nch =10,nwch = 3,ncoe = 5
real(r_kind) :: tbb(*),theta
real(r_kind) :: em_vector(*),emissivity,ts,frequency,ed0(nwch),discriminator(5)
integer(i_kind) :: snow_type,i,k,ich,nvalid_ch
real(r_kind) :: coe31(0:ncoe),coe89(0:ncoe),coe150(0:ncoe)
real(r_kind) :: coe(nch*(ncoe+1))
Equivalence (coe(1),coe31)
Equivalence (coe(11),coe89)
Equivalence (coe(21),coe150)
! Fitting Coefficients at 31.4 GHz: Using Tb4, Tb5 and Ts
data coe31/ 3.110967e-001_r_kind, 1.100175e-002_r_kind, -1.677626e-005_r_kind, &
-4.020427e-003_r_kind, 9.242240e-006_r_kind, -2.363207e-003_r_kind/
! Fitting Coefficients at 89 GHz: Using Tb4, Tb5 and Ts
data coe89/ 1.148098e+000_r_kind, 1.452926e-003_r_kind, 1.037081e-005_r_kind, &
1.340696e-003_r_kind, -5.185640e-006_r_kind, -4.546382e-003_r_kind /
! Fitting Coefficients at 150 GHz: Using Tb4, Tb5 and Ts
data coe150/ 1.165323e+000_r_kind, -1.030435e-003_r_kind, 4.828009e-006_r_kind, &
4.851731e-003_r_kind, -2.588049e-006_r_kind, -4.990193e-003_r_kind/
save coe31,coe89,coe150
! Calculate emissivity discriminators at five AMSU window channels
do ich = 1, nwch
ed0(ich) = coe(1+(ich-1)*10)
nvalid_ch = 2
do i=1,nvalid_ch
ed0(ich) = ed0(ich) + coe((ich-1)*10 + 2*i)*tbb(i) + &
coe((ich-1)*10 + 2*i+1)*tbb(i)*tbb(i)
end do
ed0(ich) = ed0(ich) + coe( (ich-1)*10 + (nvalid_ch+1)*2 )*ts
end do
! Quality control
if(ed0(2) .gt. ed0(1)) &
ed0(2) = ed0(1) + (150.0_r_kind - 89.0_r_kind)*(ed0(3) - ed0(1)) / &
(150.0_r_kind - 31.4_r_kind)
! Match the format of the input variable
! Missing value at 23.8 GHz
discriminator(1) = -999.9_r_kind; discriminator(2) = ed0(1)
! Missing value at 50.3 GHz
discriminator(3) = -999.9_r_kind; discriminator(4) = ed0(2); discriminator(5) = ed0(3)
call em_interpolate(frequency,discriminator,emissivity,snow_type)
em_vector(1) = emissivity
em_vector(2) = emissivity
return
end subroutine sem_BTs
!---------------------------------------------------------------------!
subroutine sem_amsub(theta,frequency,tbb,snow_type,em_vector)
!$$$ subprogram documentation block
!
! subprogram:
!
! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
!
! abstract:
! Calculate the emissivity discriminators and interpolate/extrapolate
! emissivity at required frequency with respect to secenery AMSUB
!
! program history log:
! 2004-10-28 treadon - correct problem in declared dimensions of array coe
!
! input argument list:
!
! frequency - frequency in GHz
! theta - local zenith angle (not used here)
! tbb[1] ~ tbb[2] - brightness temperature at five AMSU-B window channels:
! tbb[1] : 89 GHz
! tbb[2] : 150 GHz
!
! output argument list:
! em_vector(1) and (2) - emissivity at two polarizations.
! set esv = esh here and will be updated
! snow_type - snow type (reference [2])
!
! important internal variables:
!
! coe31 - fitting coefficients to estimate discriminator at 31.4 GHz
! coe89 - fitting coefficients to estimate discriminator at 89 GHz
! coe150 - fitting coefficients to estimate discriminator at 150 GHz
!
! remarks:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$
! use kinds, only: r_kind,i_kind
implicit none
integer(i_kind),parameter:: nch =10,nwch = 3,ncoe = 4
real(r_kind) :: tbb(*)
real(r_kind) :: em_vector(*),emissivity,frequency,ed0(nwch),discriminator(5)
integer(i_kind) :: snow_type,i,k,ich,nvalid_ch
real(r_kind) :: coe31(0:ncoe),coe89(0:ncoe),coe150(0:ncoe)
real(r_kind) :: coe(nch*(ncoe+1))
real(r_kind) :: theta,dem,demmin0
Equivalence (coe(1),coe31)
Equivalence (coe(11),coe89)
Equivalence (coe(21),coe150)
! Fitting Coefficients at 31.4 GHz: Using Tb4, Tb5
data coe31/-4.015636e-001_r_kind,9.297894e-003_r_kind, -1.305068e-005_r_kind, &
3.717131e-004_r_kind, -4.364877e-006_r_kind/
! Fitting Coefficients at 89 GHz: Using Tb4, Tb5
data coe89/-2.229547e-001_r_kind, -1.828402e-003_r_kind,1.754807e-005_r_kind, &
9.793681e-003_r_kind, -3.137189e-005_r_kind/
! Fitting Coefficients at 150 GHz: Using Tb4, Tb5
data coe150/-3.395416e-001_r_kind,-4.632656e-003_r_kind,1.270735e-005_r_kind, &
1.413038e-002_r_kind,-3.133239e-005_r_kind/
save coe31,coe89,coe150
! Calculate emissivity discriminators at five AMSU window channels
do ich = 1, nwch
ed0(ich) = coe(1+(ich-1)*10)
nvalid_ch = 2
do i=1,nvalid_ch
ed0(ich) = ed0(ich) + coe((ich-1)*10 + 2*i)*tbb(i) + &
coe((ich-1)*10 + 2*i+1)*tbb(i)*tbb(i)
end do
end do
! Quality Control
if(ed0(2) .gt. ed0(1)) &
ed0(2) = ed0(1) + (150.0_r_kind - 89.0_r_kind) * &
(ed0(3) - ed0(1))/(150.0_r_kind - 31.4_r_kind)
! Match the format of the input variable
! Missing value at 23.8 GHz
discriminator(1) = -999.9_r_kind; discriminator(2) = ed0(1)
! Missing value at 50.3 GHz
discriminator(3) = -999.9_r_kind; discriminator(4) = ed0(2); discriminator(5) = ed0(3)
call em_interpolate(frequency,discriminator,emissivity,snow_type)
em_vector(1) = emissivity
em_vector(2) = emissivity
return
end subroutine sem_amsub
!---------------------------------------------------------------------!
subroutine ALandEM_Snow(theta,frequency,snow_depth,t_skin,snow_type,em_vector)
!$$$ subprogram documentation block
!
! subprogram:
!
! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
!
! abstract:
! Calculate the emissivity at required frequency with respect to option MODL
! using the LandEM and a bias correction algorithm, where the original LandEM with a
! bias correction algorithm is referred to as value-added LandEM or AlandEM.
!
! program history log:
!
! input argument list:
!
! frequency - frequency in GHz
! theta - local zenith angle (not used here)
! snow_depth - snow depth in mm
! t_skin - surface temperature
!
! output argument list:
!
! em_vector(1) and (2) - emissivity at two polarizations.
! set esv = esh here and will be updated
! snow_type - snow type
!
! important internal variables:
!
! esv_3w and esh_3w - initial emissivity discriminator at two polarizations
! at three AMSU window channels computed using LandEM
! esv_3w[1] and esh_3w[1] : 31.4 GHz
! esv_3w[2] and esh_3w[2] : 89 GHz
! esv_3w[3] and esh_3w[3] : 150 GHz
!
! remarks:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$
! use kinds, only: r_kind,i_kind
! use constants, only: zero, one
implicit none
integer(i_kind) :: nw_ind
parameter(nw_ind=3)
real(r_kind) theta, frequency, freq,snow_depth, mv, t_soil, t_skin, em_vector(2)
real(r_kind) esv,esh,esh0,esv0,theta0,b
integer(i_kind) snow_type,ich
real(r_kind) freq_3w(nw_ind),esh_3w(nw_ind),esv_3w(nw_ind)
complex(r_kind) eair
data freq_3w/31.4_r_kind,89.0_r_kind,150.0_r_kind/
eair = dcmplx(one,-zero)
b = t_skin
snow_type = -999
call snowem_default(theta,frequency,snow_depth,t_skin,b,esv0,esh0)
theta0 = theta
do ich = 1, nw_ind
freq =freq_3w(ich)
theta = theta0
call snowem_default(theta,freq,snow_depth,t_skin,b,esv,esh)
esv_3w(ich) = esv
esh_3w(ich) = esh
end do
call ems_adjust(theta,frequency,snow_depth,t_skin,esv_3w,esh_3w,em_vector,snow_type)
return
end subroutine ALandEM_Snow
!---------------------------------------------------------------------!
subroutine snowem_default(theta,freq,snow_depth,t_soil,b,esv,esh)
!$$$ subprogram documentation block
!
! subprogram:
!
! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
!
! abstract:
! Initialize discriminator using LandEM
!
! program history log:
! 2004-09-22 todling - replaced zsqrt by general interface sqrt
!
! input argument list:
!
! frequency - frequency in GHz
! theta - local zenith angle in radian
! snow_depth - snow depth in mm
! t_skin - surface temperature
!
! output argument list:
!
! esv - initial discriminator at vertical polarization
! esh - at horizontal polarization
!
! remarks:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$
! use kinds, only: r_kind
! use constants, only: zero, one
implicit none
real(r_kind) rhob,rhos,sand,clay
Parameter(rhob = 1.18_r_kind, rhos = 2.65_r_kind, &
sand = 0.8_r_kind, clay = 0.2_r_kind)
real(r_kind) theta, freq, mv, t_soil, snow_depth,b
real(r_kind) theta_i,theta_t, mu, r12_h, r12_v, r21_h, r21_v, r23_h, r23_v, &
t21_v, t21_h, t12_v, t12_h, gv, gh, ssalb_h,ssalb_v,tau_h, &
tau_v, esh, esv,rad, sigma, va ,ep_real,ep_imag
complex(r_kind) esoil, esnow, eair
eair = dcmplx(one,-zero)
! ep = dcmplx(3.2_r_kind,-0.0005_r_kind)
sigma = one
theta_i = theta
mv = 0.1_r_kind
ep_real = 3.2_r_kind
ep_imag = -0.0005_r_kind
va = 0.4_r_kind + 0.0004_r_kind*snow_depth
rad = one + 0.005_r_kind*snow_depth
call snow_diel(freq, ep_real, ep_imag, rad, va, esnow)
! call snow_diel(freq, ep, rad, va, esnow)
call soil_diel(freq, t_soil, mv, rhob, rhos, sand, clay, esoil)
theta_t = asin(real(sin(theta_i)*sqrt(eair)/sqrt(esnow)))
call reflectance(eair, esnow, theta_i, theta_t, r12_v, r12_h)
call transmitance(eair, esnow, theta_i, theta_t, t12_v, t12_h)
theta_t = theta
theta_i = asin(real(sin(theta_t)*sqrt(eair)/sqrt(esnow)))
call reflectance(esnow, eair, theta_i, theta_t, r21_v, r21_h)
call transmitance(esnow, eair, theta_i, theta_t, t21_v, t21_h)
mu = cos(theta_i)
theta_t = asin(real(sin(theta_i)*sqrt(esnow)/sqrt(esoil)))
call reflectance(esnow, esoil, theta_i, theta_t, r23_v, r23_h)
call rough_reflectance(freq, theta_i, sigma, r23_v, r23_h)
! call snow_optic(freq, rad, snow_depth, va, ep, gv, gh, ssalb_v, ssalb_h, tau_v, tau_h)
call snow_optic(freq,rad,snow_depth,va,ep_real, ep_imag,gv,gh,&
ssalb_v,ssalb_h,tau_v,tau_h)
call two_stream_solution(b,mu,gv,gh,ssalb_h, ssalb_v, tau_h, tau_v, r12_h, &
r12_v, r21_h, r21_v, r23_h, r23_v, t21_v, t21_h, t12_v, t12_h, esv, esh)
return
end subroutine snowem_default
!---------------------------------------------------------------------!
subroutine ems_adjust(theta,frequency,depth,ts,esv_3w,esh_3w,em_vector,snow_type)
!$$$ subprogram documentation block
!
! subprogram:
!
! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
!
! abstract:
! Calculate the emissivity discriminators and interpolate/extrapolate
! emissivity at required frequency with respect to secenery MODL
!
! program history log:
! 2004-10-28 treadon - remove nch from parameter declaration below (not used)
!
! input argument list:
!
! frequency - frequency in GHz
! theta - local zenith angle (not used here)
! depth - snow depth in mm
! ts - surface temperature
!
! output argument list:
!
! em_vector(1) and (2) - emissivity at two polarizations.
! set esv = esh here and will be updated
! snow_type - snow type
!
! important internal variables:
!
! dem_coe - fitting coefficients to compute discriminator correction value
! dem_coe[1,*] : 31.4 GHz
! dem_coe[2,*] : 89 GHz
! dem_coe[3,*] : 150 GHz
!
! remarks:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$
! use kinds, only: r_kind,r_double,i_kind
! use constants, only: one,deg2rad
implicit none
integer(i_kind),parameter:: nw_3=3
integer(i_kind),parameter:: ncoe=6
real(r_kind),parameter :: earthrad = 6374._r_kind, satheight = 833.4_r_kind
integer(i_kind) :: snow_type,ich,j,k
real(r_kind) :: theta,frequency,depth,ts,esv_3w(*),esh_3w(*)
real(r_kind) :: discriminator(5),emmod(nw_3),dem(nw_3)
real(r_kind) :: emissivity,em_vector(2)
real(r_double) :: dem_coe(nw_3,0:ncoe-1),sinthetas,costhetas
save dem_coe
data (dem_coe(1,k),k=0,ncoe-1)/ 2.306844e+000_r_double, -7.287718e-003_r_double, &
-6.433248e-004_r_double, 1.664216e-005_r_double, &
4.766508e-007_r_double, -1.754184e+000_r_double/
data (dem_coe(2,k),k=0,ncoe-1)/ 3.152527e+000_r_double, -1.823670e-002_r_double, &
-9.535361e-004_r_double, 3.675516e-005_r_double, &
9.609477e-007_r_double, -1.113725e+000_r_double/
data (dem_coe(3,k),k=0,ncoe-1)/ 3.492495e+000_r_double, -2.184545e-002_r_double, &
6.536696e-005_r_double, 4.464352e-005_r_double, &
-6.305717e-008_r_double, -1.221087e+000_r_double/
sinthetas = sin(theta*deg2rad)* earthrad/(earthrad + satheight)
sinthetas = sinthetas*sinthetas
costhetas = one - sinthetas
do ich = 1, nw_3
emmod(ich) = costhetas*esv_3w(ich) + sinthetas*esh_3w(ich)
end do
do ich=1,nw_3
dem(ich) = dem_coe(ich,0) + dem_coe(ich,1)*ts + dem_coe(ich,2)*depth + &
dem_coe(ich,3)*ts*ts + dem_coe(ich,4)*depth*depth + &
dem_coe(ich,5)*emmod(ich)
end do
emmod(1) = emmod(1) + dem(1)
emmod(2) = emmod(2) + dem(2)
emmod(3) = emmod(3) + dem(3)
! Match the format of the input variable
! Missing value at 23.8 GHz
discriminator(1) = -999.9_r_kind; discriminator(2) = emmod(1)
! Missing value at 50.3 GHz
discriminator(3) = -999.9_r_kind; discriminator(4) = emmod(2); discriminator(5) = emmod(3)
call em_interpolate(frequency,discriminator,emissivity,snow_type)
em_vector(1) = emissivity
em_vector(2) = emissivity
return
end subroutine ems_adjust