!
! NESDIS_AMSU_SnowEM_Module
!
! Module containing the AMSU microwave snow emissivity model
!
! References:
! Yan,B., F.Weng, and K.Okamoto, 2004, A microwave snow emissivity model,
! 8th Specialist Meeting on Microwave Radiometry and Remote Sensing Applications,
! 24-27 February, 2004, Rome, Italy.
!
!
! CREATION HISTORY:
! Written by: Banghua Yan, 03-Jun-2005, banghua.yan@noaa.gov
! Fuzhong Weng, fuzhong.weng@noaa.gov
!
MODULE NESDIS_AMSU_SnowEM_Module
! -----------------
! Environment setup
! -----------------
! Module use
USE Type_Kinds, ONLY: fp, Double
USE NESDIS_LandEM_Module
USE NESDIS_SnowEM_Parameters
! Disable implicit typing
IMPLICIT NONE
! ------------
! Visibilities
! ------------
PRIVATE
PUBLIC :: NESDIS_AMSU_SNOWEM
! -----------------
! Module parameters
! -----------------
! Version Id for the module
CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
'$Id: NESDIS_AMSU_SnowEM_Module.f90 99117 2017-11-27 18:37:14Z tong.zhu@noaa.gov $'
CONTAINS
!################################################################################
!################################################################################
!## ##
!## ## PUBLIC MODULE ROUTINES ## ##
!## ##
!################################################################################
!################################################################################
!-------------------------------------------------------------------------------------------------------------
!
! NAME:
! NESDIS_AMSU_SNOWEM
!
! PURPOSE:
! Subroutine to simulate microwave emissivity over snow conditions from AMSU measurements at window
! channels.
!
! REFERENCES:
! Yan, B., F. Weng and K.Okamoto,2004: "A microwave snow emissivity model, 8th Specialist Meeting on
! Microwave Radiometry and Remote Sension Applications,24-27 February, 2004, Rome, Italy.
!
! CATEGORY:
! CRTM : Surface : MW SNOWEM
!
! LANGUAGE:
! Fortran-95
!
! CALLING SEQUENCE:
! CALL NESDIS_AMSU_SNOWEM
!
! INPUT ARGUMENTS:
!
! Frequency Frequency User defines
! This is the "I" dimension
! UNITS: GHz
! TYPE: REAL( fp )
! DIMENSION: Scalar
!
!
! Satellite_Angle The local zenith angle in degree for AMSU measurements.
! ** NOTE: THIS IS A MANDATORY MEMBER **
! ** OF THIS STRUCTURE **
! UNITS: Degrees
! TYPE: REAL( fp )
! DIMENSION: Rank-1, (I)
!
! User_Angle The local angle value in degree user defines.
! ** NOTE: THIS IS A MANDATORY MEMBER **
! ** OF THIS STRUCTURE **
! UNITS: Degrees
! TYPE: REAL( fp )
! DIMENSION: Rank-1, (I)
!
!
! Tba BRIGHTNESS TEMPERATURES AT FOUR AMSU-A WINDOW CHANNELS
! UNITS: Kelvin, K
! TYPE: REAL( fp )
! DIMENSION 4*1 SCALAR
!
! WHICH ARE
! tba[1] = TB at 23.8 GHz
! tba[2] = TB at: 31.4 GHz
! tba[3] = TB at 50.3 GHz
! tba[4] = TB at 89 GHz
!
! Tbb BRIGHTNESS TEMPERATURES AT TWO AMSU-B WINDOW CHANNELS
! UNITS: Kelvin, K
! TYPE: REAL( fp )
! DIMENSION 2*1 SCALAR
!
! WHICH ARE
!
! tbb[1] = TB at 89 GHz
! tbb[2] = TB at 150 GHz
!
!
! Ts = Land_Temperature: The land surface temperature.
! UNITS: Kelvin, K
! TYPE: REAL( fp )
! DIMENSION: Scalar
!
!
! Snow_Depth: The snow depth.
! UNITS: mm
! TYPE: REAL( fp )
! DIMENSION: Scalar
!
! **** IMPORTANT NOTES ****
!
! When one variable among Tba[], Tbb[] and Ts are not available, set -999.0
!
!
!
!
! OUTPUT ARGUMENTS:
!
! Emissivity_H: The surface emissivity at a horizontal polarization.
! ** NOTE: THIS IS A MANDATORY MEMBER **
! ** OF THIS STRUCTURE **
! UNITS: N/A
! TYPE: REAL( fp )
! DIMENSION: Scalar
!
! Emissivity_V: The surface emissivity at a vertical polarization.
! ** NOTE: THIS IS A MANDATORY MEMBER **
! ** OF THIS STRUCTURE **
! UNITS: N/A
! TYPE: REAL( fp )
! DIMENSION: Scalar
!
!
! IINTERNAL ARGUMENTS:
!
! input_type (specific option index = 1 ~ 7):
!
! input_type = 1 : AMSU-A & B window channels of Tb and Ts are available (call AMSU_ABTs)
!
! input_type = 2 : AMSU-A window channels of Tb and Ts are available (call AMSU_ATs)
!
! input_type = 3 : AMSU-A & B window channels of Tb are available (call AMSU_AB)
!
! input_type = 4 : AMSU-A window channels of Tb are available (call AMSU_amsua)
!
! input_type = 5 : AMSU-B window channels of Tb and Ts are available (call AMSU_BTs)
!
! input_type = 6 : AMSU-B window channels of Tb are available (call AMSU_amsub)
!
! input_type = 7 : snow depth and Ts are available (call ALandEM_Snow)
!
!
! OPTIONAL OUTPUT ARGUMENTS:
!
! 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
!
! CALLS:
! AMSU_ABTs : Subroutine to calculate the microwave snow emissivity from AMSU-A/B TB and Ts
!
! AMSU_ATs : Subroutine to calculate the microwave snow emissivity from AMSU-A TB and Ts
!
! AMSU_AB : Subroutine to calculate the microwave snow emissivity from AMSU-A/B TB
!
! AMSU_amsua : Subroutine to calculate the microwave snow emissivity from AMSU-A TB
!
! AMSU_BTs : Subroutine to calculate the microwave snow emissivity from AMSU-B TB and Ts
!
! AMSU_amsub : Subroutine to calculate the microwave snow emissivity from AMSU-B TB
!
! AMSU_ALandEM_Snow : Subroutine to calculate the microwave snow emissivity from Ts and Snow Depth
!
! em_initialization : Subroutine to initialization snow emissivity
!
!
! SIDE EFFECTS:
! None.
!
! RESTRICTIONS:
! None.
!
! COMMENTS:
! Note the INTENT on the output SensorData argument is IN OUT rather than
! just OUT. This is necessary because the argument may be defined upon
! input. To prevent memory leaks, the IN OUT INTENT is a must.
!
! CREATION HISTORY:
! Written by: Banghua Yan, QSS Group Inc., Banghua.Yan@noaa.gov (03-June-2005)
!
!
! and Fuzhong Weng, NOAA/NESDIS/ORA, Fuzhong.Weng@noaa.gov
!
! Copyright (C) 2005 Fuzhong Weng and Banghua Yan
!
! This program is free software; you can redistribute it and/or modify it under the terms of the GNU
! General Public License as published by the Free Software Foundation; either version 2 of the License,
! or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even
! the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
! License for more details.
!
! You should have received a copy of the GNU General Public License along with this program; if not, write
! to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
!
!------------------------------------------------------------------------------------------------------------
subroutine NESDIS_AMSU_SNOWEM(Satellite_Angle, & ! INPUT
User_Angle, & ! INPUT
frequency, & ! INPUT
Snow_Depth, & ! INPUT
Ts, & ! INPUT
tba, & ! INPUT
tbb, & ! INPUT
Emissivity_H, & ! OUPUT
Emissivity_V) ! OUTPUT
integer,PARAMETER :: AMSU_ABTs_ALG = 1
integer,PARAMETER :: AMSU_ATs_ALG = 2
integer,PARAMETER :: AMSU_AB_ALG = 3
integer,PARAMETER :: AMSU_amsua_ALG = 4
integer,PARAMETER :: AMSU_BTs_ALG = 5
integer,PARAMETER :: AMSU_amsub_ALG = 6
integer,PARAMETER :: AMSU_ALandEM_Snow_ALG = 7
integer, parameter :: nch = 10, nwcha = 4, nwchb = 2, nwch = 5,nalg = 7
integer :: snow_type,input_type,i,np,k
real(fp) :: Satellite_Angle,User_Angle,frequency,Ts,Snow_Depth
real(fp) :: em_vector(2),esh1,esv1,esh2,esv2,desh,desv,dem
real(fp) :: tb(nwch),tba(nwcha),tbb(nwchb)
real(fp), intent(out) :: Emissivity_H,Emissivity_V
logical :: INDATA(nalg)
! Initialization
call em_initialization(frequency,em_vector)
snow_type = INVALID_SNOW_TYPE
input_type = INVALID_SNOW_TYPE
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 <= 150.0_fp) .or. (Ts >= 290.0_fp) ) then
INDATA(1:2) = .false.; INDATA(5) = .false.; INDATA(7) = .false.
end if
do i=1,nwcha
if((tba(i) <= 100.0_fp) .or. (tba(i) >= 290.0_fp) ) then
INDATA(1:4) = .false.
exit
end if
end do
do i=1,nwchb
if((tbb(i) <= 100.0_fp) .or. (tbb(i) >= 290.0_fp) ) then
INDATA(1) = .false.; INDATA(3) = .false.; INDATA(5:6) = .false.
exit
end if
end do
if((Snow_Depth < 0.0_fp) .or. (Snow_Depth >= 3000.0_fp)) INDATA(7) = .false.
if((frequency >= 80._fp) .and. (INDATA(5))) then
INDATA(2:3) = .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
! GET EMISSIVITY AT SATELLITE'S ZENITH ANGLE
GET_option: SELECT CASE (input_type)
CASE (AMSU_ABTs_ALG)
call AMSU_ABTs(frequency,tb,Ts,snow_type,em_vector)
CASE (AMSU_ATs_ALG)
call AMSU_ATs(frequency,tba,Ts,snow_type,em_vector)
CASE (AMSU_AB_ALG)
call AMSU_AB(frequency,tb,snow_type,em_vector)
CASE (AMSU_amsua_ALG)
call AMSU_amsua(frequency,tba,snow_type,em_vector)
CASE(AMSU_BTs_ALG)
call AMSU_BTs(frequency,tbb,Ts,snow_type,em_vector)
CASE(AMSU_amsub_ALG)
call AMSU_amsub(frequency,tbb,snow_type,em_vector)
CASE(AMSU_ALandEM_Snow_ALG)
call AMSU_ALandEM_Snow(Satellite_Angle,frequency,Snow_Depth,Ts,snow_type,em_vector)
END SELECT GET_option
! Get the emissivity angle dependence
call NESDIS_LandEM(Satellite_Angle,frequency,0.0_fp,0.0_fp,Ts,Ts,0.0_fp,9,13,2.0_fp,esh1,esv1)
call NESDIS_LandEM(User_Angle,frequency,0.0_fp,0.0_fp,Ts,Ts,0.0_fp,9,13,2.0_fp,esh2,esv2)
desh = esh1 - esh2
desv = esv1 - esv2
dem = ( desh + desv ) * 0.5_fp
! Emissivity at User's Angle
Emissivity_H = em_vector(1) - dem; Emissivity_V = em_vector(2)- dem
if (Emissivity_H > one) Emissivity_H = one
if (Emissivity_V > one) Emissivity_V = one
if (Emissivity_H < 0.3_fp) Emissivity_H = 0.3_fp
if (Emissivity_V < 0.3_fp) Emissivity_V = 0.3_fp
return
end subroutine NESDIS_AMSU_SNOWEM
!##################################################################################
!##################################################################################
!## ##
!## ## PRIVATE MODULE ROUTINES ## ##
!## ##
!##################################################################################
!##################################################################################
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
!
!----------------------------------------------------------------------------------------------------------!
integer :: nch,ncand
Parameter(nch = 10,ncand=16)
real(fp) :: frequency,em_vector(*),freq(nch)
real(fp) :: em(ncand,nch)
real(fp) :: kratio, bconst,emissivity
integer :: ich
! Sixteen candidate snow emissivity spectra
em(1, 1: N_FREQUENCY) = WET_SNOW_EMISS(1:N_FREQUENCY)
em(2, 1: N_FREQUENCY) = GRASS_AFTER_SNOW_EMISS(1:N_FREQUENCY)
em(3, 1: N_FREQUENCY) = RS_SNOW_A_EMISS(1:N_FREQUENCY)
em(4, 1: N_FREQUENCY) = POWDER_SNOW_EMISS(1:N_FREQUENCY)
em(5, 1: N_FREQUENCY) = RS_SNOW_B_EMISS(1:N_FREQUENCY)
em(6, 1: N_FREQUENCY) = RS_SNOW_C_EMISS(1:N_FREQUENCY)
em(7, 1: N_FREQUENCY) = RS_SNOW_D_EMISS(1:N_FREQUENCY)
em(8, 1: N_FREQUENCY) = THIN_CRUST_SNOW_EMISS(1:N_FREQUENCY)
em(9, 1: N_FREQUENCY) = RS_SNOW_E_EMISS(1:N_FREQUENCY)
em(10, 1: N_FREQUENCY) = BOTTOM_CRUST_SNOW_A_EMISS(1:N_FREQUENCY)
em(11, 1: N_FREQUENCY) = SHALLOW_SNOW_EMISS(1:N_FREQUENCY)
em(12, 1: N_FREQUENCY) = DEEP_SNOW_EMISS(1:N_FREQUENCY)
em(13, 1: N_FREQUENCY) = CRUST_SNOW_EMISS(1:N_FREQUENCY)
em(14, 1: N_FREQUENCY) = MEDIUM_SNOW_EMISS(1:N_FREQUENCY)
em(15, 1: N_FREQUENCY) = BOTTOM_CRUST_SNOW_B_EMISS(1:N_FREQUENCY)
em(16, 1: N_FREQUENCY) = THICK_CRUST_SNOW_EMISS(1:N_FREQUENCY)
freq = FREQUENCY_DEFAULT
! 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_fp) emissivity = 0.8_fp
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:
!
! emissivity - weighted emissivity from both V- and H- Pols.
! 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
!
!----------------------------------------------------------------------------------------------------------!
integer,parameter:: ncand = 16,nch =10
integer:: ich,ichmin,ichmax,i,k,snow_type
real(fp) :: dem,demmin0
real(fp) :: em(ncand,nch)
real(fp) :: frequency,freq(nch),emissivity,discriminator(*),emis(nch)
real(fp) :: cor_factor,adjust_check,kratio, bconst
! Sixteen candidate snow emissivity spectra
em(1, 1: N_FREQUENCY) = WET_SNOW_EMISS(1:N_FREQUENCY)
em(2, 1: N_FREQUENCY) = GRASS_AFTER_SNOW_EMISS(1:N_FREQUENCY)
em(3, 1: N_FREQUENCY) = RS_SNOW_A_EMISS(1:N_FREQUENCY)
em(4, 1: N_FREQUENCY) = POWDER_SNOW_EMISS(1:N_FREQUENCY)
em(5, 1: N_FREQUENCY) = RS_SNOW_B_EMISS(1:N_FREQUENCY)
em(6, 1: N_FREQUENCY) = RS_SNOW_C_EMISS(1:N_FREQUENCY)
em(7, 1: N_FREQUENCY) = RS_SNOW_D_EMISS(1:N_FREQUENCY)
em(8, 1: N_FREQUENCY) = THIN_CRUST_SNOW_EMISS(1:N_FREQUENCY)
em(9, 1: N_FREQUENCY) = RS_SNOW_E_EMISS(1:N_FREQUENCY)
em(10, 1: N_FREQUENCY) = BOTTOM_CRUST_SNOW_A_EMISS(1:N_FREQUENCY)
em(11, 1: N_FREQUENCY) = SHALLOW_SNOW_EMISS(1:N_FREQUENCY)
em(12, 1: N_FREQUENCY) = DEEP_SNOW_EMISS(1:N_FREQUENCY)
em(13, 1: N_FREQUENCY) = CRUST_SNOW_EMISS(1:N_FREQUENCY)
em(14, 1: N_FREQUENCY) = MEDIUM_SNOW_EMISS(1:N_FREQUENCY)
em(15, 1: N_FREQUENCY) = BOTTOM_CRUST_SNOW_B_EMISS(1:N_FREQUENCY)
em(16, 1: N_FREQUENCY) = THICK_CRUST_SNOW_EMISS(1:N_FREQUENCY)
freq = FREQUENCY_DEFAULT
! Adjust unreasonable discriminator
if (discriminator(4) > discriminator(2)) &
discriminator(4) = discriminator(2) +(discriminator(5) - discriminator(2))* &
(150.0_fp - 89.0_fp)/(150.0_fp - 31.4_fp)
if ( (discriminator(3) /= -999.9_fp) .and. &
( ((discriminator(3)-0.01_fp) > discriminator(2)) .or. &
((discriminator(3)-0.01_fp) < discriminator(4))) ) &
discriminator(3) = discriminator(2) + &
(discriminator(4) - discriminator(2))*(89.0_fp - 50.3_fp) &
/ (89.0_fp - 31.4_fp)
! Find a snow emissivity spectrum
if(snow_type .eq. -999) then
demmin0 = 10.0_fp
do k = 1, ncand
dem = zero
ichmin = 1
ichmax = 3
if(discriminator(1) == -999.9_fp) 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_fp) emis(ich) = 0.3_fp
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_fp) &
adjust_check = adjust_check + abs(emis(ich) - discriminator(ich - 4))
else
if (discriminator(ich - 4) .ne. -999.9_fp) &
adjust_check = adjust_check + abs(emis(ich+1) - discriminator(ich - 4))
end if
end do
if (adjust_check >= 0.04_fp) then
if (discriminator(1) /= -999.9_fp) then
if (discriminator(1) < emis(4)) then
emis(5) = emis(4) + &
(31.4_fp - 23.8_fp) * &
(discriminator(2) - emis(4))/(31.4_fp - 18.7_fp)
else
emis(5) = discriminator(1)
end if
end if
emis(6) = discriminator(2)
if (discriminator(3) /= -999.9_fp) then
emis(7) = discriminator(3)
else
! In case of missing the emissivity discriminator at 50.3 GHz
emis(7) = emis(6) + (89.0_fp - 50.3_fp) * &
(discriminator(4) - emis(6))/(89.0_fp - 31.4_fp)
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_fp) emissivity = 0.8_fp
end if
! Assume emissivity = constant at frequencies >= 150 GHz
if (frequency >= freq(nch)) emissivity = emis(nch)
return
end subroutine em_interpolate
subroutine AMSU_ABTs(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 (currently, 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
!
!----------------------------------------------------------------------------------------------------------!
integer,parameter:: ncand = 16,nch =10,nthresh=38
integer,parameter:: nind=6,ncoe=8,nLIcoe=6,nHIcoe=12
integer:: i,j,k,num,npass,snow_type,md0,md1,nmodel(ncand-1)
real(fp) :: frequency,tb150,LI,HI,DS1,DS2,DS3
real(fp) :: em(ncand,nch), em_vector(*)
real(fp) :: 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(fp) :: LI_coe(0:nLIcoe-1),HI_coe(0:nHIcoe-1)
real(fp) :: ts,emissivity
real(fp) :: discriminator(5)
logical:: pick_status,tindex(nind)
save threshold,DI_coe,LI_coe, HI_coe,nmodel
data nmodel/5,10,13,16,18,24,30,31,32,33,34,35,36,37,38/
! Fitting coefficients for five discriminators
DI_coe(1,0:ncoe-1) = (/ &
3.285557e-002_fp, 2.677179e-005_fp, &
4.553101e-003_fp, 5.639352e-005_fp, &
-1.825188e-004_fp, 1.636145e-004_fp, &
1.680881e-005_fp, -1.708405e-004_fp/)
DI_coe(2,0:ncoe-1) = (/ &
-4.275539e-002_fp, -2.541453e-005_fp, &
4.154796e-004_fp, 1.703443e-004_fp, &
4.350142e-003_fp, 2.452873e-004_fp, &
-4.748506e-003_fp, 2.293836e-004_fp/)
DI_coe(3,0:ncoe-1) = (/ &
-1.870173e-001_fp, -1.061678e-004_fp, &
2.364055e-004_fp, -2.834876e-005_fp, &
4.899651e-003_fp, -3.418847e-004_fp, &
-2.312224e-004_fp, 9.498600e-004_fp/)
DI_coe(4,0:ncoe-1) = (/ &
-2.076519e-001_fp, 8.475901e-004_fp, &
-2.072679e-003_fp, -2.064717e-003_fp, &
2.600452e-003_fp, 2.503923e-003_fp, &
5.179711e-004_fp, 4.667157e-005_fp/)
DI_coe(5,0:ncoe-1) = (/ &
-1.442609e-001_fp, -8.075003e-005_fp, &
-1.790933e-004_fp, -1.986887e-004_fp, &
5.495115e-004_fp, -5.871732e-004_fp, &
4.517280e-003_fp, 7.204695e-004_fp/)
! Fitting coefficients for emissivity index at 31.4 GHz
LI_coe = (/ &
7.963632e-001_fp, 7.215580e-003_fp, &
-2.015921e-005_fp, -1.508286e-003_fp, &
1.731405e-005_fp, -4.105358e-003_fp/)
! Fitting coefficients for emissivity index at 150 GHz
HI_coe = (/ &
1.012160e+000_fp, 6.100397e-003_fp, &
-1.774347e-005_fp, -4.028211e-003_fp, &
1.224470e-005_fp, 2.345612e-003_fp, &
-5.376814e-006_fp, -2.795332e-003_fp, &
8.072756e-006_fp, 3.529615e-003_fp, &
1.955293e-006_fp, -4.942230e-003_fp/)
! Six thresholds for sixteen candidate snow types
! Note: some snow type contains several possible
! selections for six thresholds
!1 Wet Snow
threshold(1,1:6) = (/0.88_fp,0.86_fp,-999.9_fp,&
0.01_fp,0.01_fp,200._fp/)
threshold(2,1:6) = (/0.88_fp,0.85_fp,-999.9_fp,&
0.06_fp,0.10_fp,200._fp/)
threshold(3,1:6) = (/0.88_fp,0.83_fp,-0.02_fp,&
0.12_fp,0.16_fp,204._fp/)
threshold(4,1:6) = (/0.90_fp,0.89_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,-999.9_fp/)
threshold(5,1:6) = (/0.92_fp,0.85_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,-999.9_fp/)
!2 Grass_after_Snow
threshold(6,1:6) = (/0.84_fp,0.83_fp,-999.9_fp,&
0.08_fp,0.10_fp,195._fp/)
threshold(7,1:6) = (/0.85_fp,0.85_fp,-999.9_fp,&
0.10_fp,-999.9_fp,190._fp/)
threshold(8,1:6) = (/0.86_fp,0.81_fp,-999.9_fp,&
0.12_fp,-999.9_fp,200._fp/)
threshold(9,1:6) = (/0.86_fp,0.81_fp,0.0_fp,&
0.12_fp,-999.9_fp,189._fp/)
threshold(10,1:6) = (/0.90_fp,0.81_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,195._fp/)
!3 RS_Snow (A)
threshold(11,1:6) = (/0.80_fp,0.76_fp,-999.9_fp,&
0.05_fp,-999.9_fp,185._fp/)
threshold(12,1:6) = (/0.82_fp,0.78_fp,-999.9_fp,&
-999.9_fp,0.25_fp,180._fp/)
threshold(13,1:6) = (/0.90_fp,0.76_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,180._fp/)
!4 Powder Snow
threshold(14,1:6) = (/0.89_fp,0.73_fp,-999.9_fp,&
0.20_fp,-999.9_fp,-999.9_fp/)
threshold(15,1:6) = (/0.89_fp,0.75_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,-999.9_fp/)
threshold(16,1:6) = (/0.93_fp,0.72_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,-999.9_fp/)
!5 RS_Snow (B)
threshold(17,1:6) = (/0.82_fp,0.70_fp,-999.9_fp,&
0.20_fp,-999.9_fp,160._fp/)
threshold(18,1:6) = (/0.83_fp,0.70_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,160._fp/)
!6 RS_Snow (C)
threshold(19,1:6) = (/0.75_fp,0.76_fp,-999.9_fp,&
0.08_fp,-999.9_fp,172._fp/)
threshold(20,1:6) = (/0.77_fp,0.72_fp,-999.9_fp,&
0.12_fp,0.15_fp,175._fp/)
threshold(21,1:6) = (/0.78_fp,0.74_fp,-999.9_fp,&
-999.9_fp,0.20_fp,172._fp/)
threshold(22,1:6) = (/0.80_fp,0.77_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,170._fp/)
threshold(23,1:6) = (/0.82_fp,-999.9_fp,-999.9_fp,&
0.15_fp,0.22_fp,170._fp/)
threshold(24,1:6) = (/0.82_fp,0.73_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,170._fp/)
!7 RS_Snow (D)
threshold(25,1:6) = (/0.75_fp,0.70_fp,-999.9_fp,&
0.15_fp,0.25_fp,167._fp/)
threshold(26,1:6) = (/0.77_fp,0.76_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,-999.9_fp/)
threshold(27,1:6) = (/0.80_fp,0.72_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,-999.9_fp/)
threshold(28,1:6) = (/0.77_fp,0.73_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,-999.9_fp/)
threshold(29,1:6) = (/0.81_fp,0.71_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,-999.9_fp/)
threshold(30,1:6) = (/0.82_fp,0.69_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,-999.9_fp/)
!8 Thin Crust Snow
threshold(31,1:6) = (/0.88_fp,0.58_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,-999.9_fp/)
!9 RS_Snow (E)
threshold(32,1:6) = (/0.73_fp,0.67_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,-999.9_fp/)
!10 Bottom Crust Snow (A)
threshold(33,1:6) = (/0.83_fp,0.66_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,-999.9_fp/)
!11 Shallow Snow
threshold(34,1:6) = (/0.82_fp,0.60_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,-999.9_fp/)
!12 Deep Snow
threshold(35,1:6) = (/0.77_fp,0.60_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,-999.9_fp/)
!13 Crust Snow
threshold(36,1:6) = (/0.77_fp,0.7_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,-999.9_fp/)
!14 Medium Snow
threshold(37,1:6) = (/-999.9_fp,0.55_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,-999.9_fp/)
!15 Bottom Crust Snow(B)
threshold(38,1:6) = (/0.74_fp,-999.9_fp,-999.9_fp,&
-999.9_fp,-999.9_fp,-999.9_fp/)
!16 Thick Crust Snow
! lowest priority: No constraints
! Sixteen candidate snow emissivity spectra
em(1, 1: N_FREQUENCY) = WET_SNOW_EMISS(1:N_FREQUENCY)
em(2, 1: N_FREQUENCY) = GRASS_AFTER_SNOW_EMISS(1:N_FREQUENCY)
em(3, 1: N_FREQUENCY) = RS_SNOW_A_EMISS(1:N_FREQUENCY)
em(4, 1: N_FREQUENCY) = POWDER_SNOW_EMISS(1:N_FREQUENCY)
em(5, 1: N_FREQUENCY) = RS_SNOW_B_EMISS(1:N_FREQUENCY)
em(6, 1: N_FREQUENCY) = RS_SNOW_C_EMISS(1:N_FREQUENCY)
em(7, 1: N_FREQUENCY) = RS_SNOW_D_EMISS(1:N_FREQUENCY)
em(8, 1: N_FREQUENCY) = THIN_CRUST_SNOW_EMISS(1:N_FREQUENCY)
em(9, 1: N_FREQUENCY) = RS_SNOW_E_EMISS(1:N_FREQUENCY)
em(10, 1: N_FREQUENCY) = BOTTOM_CRUST_SNOW_A_EMISS(1:N_FREQUENCY)
em(11, 1: N_FREQUENCY) = SHALLOW_SNOW_EMISS(1:N_FREQUENCY)
em(12, 1: N_FREQUENCY) = DEEP_SNOW_EMISS(1:N_FREQUENCY)
em(13, 1: N_FREQUENCY) = CRUST_SNOW_EMISS(1:N_FREQUENCY)
em(14, 1: N_FREQUENCY) = MEDIUM_SNOW_EMISS(1:N_FREQUENCY)
em(15, 1: N_FREQUENCY) = BOTTOM_CRUST_SNOW_B_EMISS(1:N_FREQUENCY)
em(16, 1: N_FREQUENCY) = THICK_CRUST_SNOW_EMISS(1:N_FREQUENCY)
freq = FREQUENCY_DEFAULT
!*** 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_fp)) tindex(2) = .false.
if((i == 5) .and. (index_in(4) > 0.20_fp) &
.and. (index_in(1) > 0.88_fp)) tindex(1) = .false.
if((i == 10) .and. (index_in(1) <= 0.83_fp)) tindex(1) = .true.
if((i == 13) .and. (index_in(2) < 0.52_fp)) 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 AMSU_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
!
!----------------------------------------------------------------------------------------------------------!
integer :: i,nind
real(fp) :: index_in(*),threshold(*)
logical :: tindex(*)
do i=1,nind
tindex(i) = .false.
if (threshold(i) .eq. -999.9_fp) 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 AMSU_AB(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:
!
! 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:
!
! coe - fitting coefficients to estimate discriminator at 23.8 ~ 150 GHz
!
! remarks:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!----------------------------------------------------------------------------------------------------------!
integer,parameter:: nch =10,nwch = 5,ncoe = 10
real(fp) :: tb(*),frequency
real(fp) :: em_vector(*),emissivity,discriminator(nwch)
integer :: i,snow_type,k,ich,nvalid_ch
real(fp) :: coe(nwch*(ncoe+1))
save coe
! Fitting Coefficients at 23.8 GHz: Using Tb1 ~ Tb3
coe(1:7) = (/&
-1.326040e+000_fp, 2.475904e-002_fp, &
-5.741361e-005_fp, -1.889650e-002_fp, &
6.177911e-005_fp, 1.451121e-002_fp, &
-4.925512e-005_fp/)
! Fitting Coefficients at 31.4 GHz: Using Tb1 ~ Tb3
coe(12:18) = (/ &
-1.250541e+000_fp, 1.911161e-002_fp, &
-5.460238e-005_fp, -1.266388e-002_fp, &
5.745064e-005_fp, 1.313985e-002_fp, &
-4.574811e-005_fp/)
! Fitting Coefficients at 50.3 GHz: Using Tb1 ~ Tb3
coe(23:29) = (/ &
-1.246754e+000_fp, 2.368658e-002_fp, &
-8.061774e-005_fp, -3.206323e-002_fp, &
1.148107e-004_fp, 2.688353e-002_fp, &
-7.358356e-005_fp/)
! Fitting Coefficients at 89 GHz: Using Tb1 ~ Tb4
coe(34:42) = (/ &
-1.278780e+000_fp, 1.625141e-002_fp, &
-4.764536e-005_fp, -1.475181e-002_fp, &
5.107766e-005_fp, 1.083021e-002_fp, &
-4.154825e-005_fp, 7.703879e-003_fp, &
-6.351148e-006_fp/)
! Fitting Coefficients at 150 GHz: Using Tb1 ~ Tb5
coe(45:55) = (/&
-1.691077e+000_fp, 3.352403e-002_fp, &
-7.310338e-005_fp, -4.396138e-002_fp, &
1.028994e-004_fp, 2.301014e-002_fp, &
-7.070810e-005_fp, 1.270231e-002_fp, &
-2.139023e-005_fp, -2.257991e-003_fp, &
1.269419e-005_fp/)
! save coe
! 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 AMSU_AB
subroutine AMSU_ATs(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:
!
! coe - fitting coefficients to estimate discriminator at 23.8 ~ 150 GHz
!
! remarks:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!-----------------------------------------------------------------------------------------------------------!
integer,parameter:: nch =10,nwch = 5,ncoe = 9
real(fp) :: tba(*)
real(fp) :: em_vector(*),emissivity,ts,frequency,discriminator(nwch)
integer :: snow_type,i,k,ich,nvalid_ch
real(fp) :: coe(nch*(ncoe+1))
save coe
! Fitting Coefficients at 23.8 GHz: Using Tb1, Tb2 and Ts
coe(1:6) = (/ &
8.210105e-001_fp, 1.216432e-002_fp, &
-2.113875e-005_fp, -6.416648e-003_fp, &
1.809047e-005_fp, -4.206605e-003_fp /)
! Fitting Coefficients at 31.4 GHz: Using Tb1, Tb2 and Ts
coe(11:16) = (/ &
7.963632e-001_fp, 7.215580e-003_fp, &
-2.015921e-005_fp, -1.508286e-003_fp, &
1.731405e-005_fp, -4.105358e-003_fp /)
! Fitting Coefficients at 50.3 GHz: Using Tb1, Tb2, Tb3 and Ts
coe(21:28) = (/ &
1.724160e+000_fp, 5.556665e-003_fp, &
-2.915872e-005_fp, -1.146713e-002_fp, &
4.724243e-005_fp, 3.851791e-003_fp, &
-5.581535e-008_fp, -5.413451e-003_fp /)
! Fitting Coefficients at 89 GHz: Using Tb1 ~ Tb4 and Ts
coe(31:40) = (/ &
9.962065e-001_fp, 1.584161e-004_fp, &
-3.988934e-006_fp, 3.427638e-003_fp, &
-5.084836e-006_fp, -6.178904e-004_fp, &
1.115315e-006_fp, 9.440962e-004_fp, &
9.711384e-006_fp, -4.259102e-003_fp /)
! Fitting Coefficients at 150 GHz: Using Tb1 ~ Tb4 and Ts
coe(41:50) = (/ &
-5.244422e-002_fp, 2.025879e-002_fp, &
-3.739231e-005_fp, -2.922355e-002_fp, &
5.810726e-005_fp, 1.376275e-002_fp, &
-3.757061e-005_fp, 6.434187e-003_fp, &
6.190403e-007_fp, -2.944785e-003_fp/)
! save coe
! 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 AMSU_ATs
subroutine AMSU_amsua(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:
!
! 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:
!
! coe - fitting coefficients to estimate discriminator at 23.8 ~ 150 GHz
!
! remarks:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!--------------------------------------------------------------------------------------------------------!
integer,parameter:: nch =10,nwch = 5,ncoe = 8
real(fp) :: tba(*)
real(fp) :: em_vector(*),emissivity,frequency,discriminator(nwch)
integer :: snow_type,i,k,ich,nvalid_ch
real(fp) :: coe(50)
save coe
! Fitting Coefficients at 23.8 GHz: Using Tb1 ~ Tb3
coe(1:7) = (/ &
-1.326040e+000_fp, 2.475904e-002_fp, -5.741361e-005_fp, &
-1.889650e-002_fp, 6.177911e-005_fp, 1.451121e-002_fp, &
-4.925512e-005_fp/)
! Fitting Coefficients at 31.4 GHz: Using Tb1 ~ Tb3
coe(11:17) = (/ &
-1.250541e+000_fp, 1.911161e-002_fp, -5.460238e-005_fp, &
-1.266388e-002_fp, 5.745064e-005_fp, 1.313985e-002_fp, &
-4.574811e-005_fp/)
! Fitting Coefficients at 50.3 GHz: Using Tb1 ~ Tb3
coe(21:27) = (/ &
-1.246754e+000_fp, 2.368658e-002_fp, -8.061774e-005_fp, &
-3.206323e-002_fp, 1.148107e-004_fp, 2.688353e-002_fp, &
-7.358356e-005_fp/)
! Fitting Coefficients at 89 GHz: Using Tb1 ~ Tb4
coe(31:39) = (/ &
-1.278780e+000_fp, 1.625141e-002_fp, -4.764536e-005_fp, &
-1.475181e-002_fp, 5.107766e-005_fp, 1.083021e-002_fp, &
-4.154825e-005_fp, 7.703879e-003_fp, -6.351148e-006_fp/)
! Fitting Coefficients at 150 GHz: Using Tb1 ~ Tb4
coe(41:49) = (/ &
-1.624857e+000_fp, 3.138243e-002_fp, -6.757028e-005_fp, &
-4.178496e-002_fp, 9.691893e-005_fp, 2.165964e-002_fp, &
-6.702349e-005_fp, 1.111658e-002_fp, -1.050708e-005_fp/)
! save coe
! 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_fp - 89.0_fp)* &
(discriminator(5) - discriminator(2))/ &
(150.0_fp - 31.4_fp)
! Quality control at 50.3 GHz
if((discriminator(3) .gt. discriminator(2)) .or. &
(discriminator(3) .lt. discriminator(4))) &
discriminator(3) = discriminator(2) + (89.0_fp - 50.3_fp)* &
(discriminator(4) - discriminator(2))/(89.0_fp - 31.4_fp)
call em_interpolate(frequency,discriminator,emissivity,snow_type)
em_vector(1) = emissivity
em_vector(2) = emissivity
return
end subroutine AMSU_amsua
subroutine AMSU_BTs(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:
!
! coe - fitting coefficients to estimate discriminator at 31.4 ~ 150 GHz
!
! remarks:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!-------------------------------------------------------------------------------------------------------!
integer,parameter:: nch =10,nwch = 3,ncoe = 5
real(fp) :: tbb(*)
real(fp) :: em_vector(*),emissivity,ts,frequency,ed0(nwch),discriminator(5)
integer :: snow_type,i,k,ich,nvalid_ch
real(fp) :: coe(nch*(ncoe+1))
save coe
! Fitting Coefficients at 31.4 GHz: Using Tb4, Tb5 and Ts
coe(1:6) = (/ 3.110967e-001_fp, 1.100175e-002_fp, -1.677626e-005_fp, &
-4.020427e-003_fp, 9.242240e-006_fp, -2.363207e-003_fp/)
! Fitting Coefficients at 89 GHz: Using Tb4, Tb5 and Ts
coe(11:16) = (/ 1.148098e+000_fp, 1.452926e-003_fp, 1.037081e-005_fp, &
1.340696e-003_fp, -5.185640e-006_fp, -4.546382e-003_fp /)
! Fitting Coefficients at 150 GHz: Using Tb4, Tb5 and Ts
coe(21:26) = (/ 1.165323e+000_fp, -1.030435e-003_fp, 4.828009e-006_fp, &
4.851731e-003_fp, -2.588049e-006_fp, -4.990193e-003_fp/)
! save coe
! 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_fp - 89.0_fp)*(ed0(3) - ed0(1)) / &
(150.0_fp - 31.4_fp)
! Match the format of the input variable
! Missing value at 23.8 GHz
discriminator(1) = -999.9_fp; discriminator(2) = ed0(1)
! Missing value at 50.3 GHz
discriminator(3) = -999.9_fp; 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 AMSU_BTs
subroutine AMSU_amsub(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:
!
! 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:
!
! coe - fitting coefficients to estimate discriminator at 31.4 ~ 150 GHz
!
! remarks:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!-------------------------------------------------------------------------------------------------------!
integer,parameter:: nch =10,nwch = 3,ncoe = 4
real(fp) :: tbb(*)
real(fp) :: em_vector(*),emissivity,frequency,ed0(nwch),discriminator(5)
integer :: snow_type,i,k,ich,nvalid_ch
real(fp) :: coe(50)
save coe
! Fitting Coefficients at 31.4 GHz: Using Tb4, Tb5
coe(1:5) = (/-4.015636e-001_fp,9.297894e-003_fp, -1.305068e-005_fp, &
3.717131e-004_fp, -4.364877e-006_fp/)
! Fitting Coefficients at 89 GHz: Using Tb4, Tb5
coe(11:15) = (/-2.229547e-001_fp, -1.828402e-003_fp,1.754807e-005_fp, &
9.793681e-003_fp, -3.137189e-005_fp/)
! Fitting Coefficients at 150 GHz: Using Tb4, Tb5
coe(21:25) = (/-3.395416e-001_fp,-4.632656e-003_fp,1.270735e-005_fp, &
1.413038e-002_fp,-3.133239e-005_fp/)
! save coe
! 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_fp - 89.0_fp) * &
(ed0(3) - ed0(1))/(150.0_fp - 31.4_fp)
! Match the format of the input variable
! Missing value at 23.8 GHz
discriminator(1) = -999.9_fp; discriminator(2) = ed0(1)
! Missing value at 50.3 GHz
discriminator(3) = -999.9_fp; 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 AMSU_amsub
subroutine AMSU_ALandEM_Snow(theta,frequency,snow_depth,ts,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 NESDIS_LandEM and a bias correction algorithm, where the original NESDIS_LandEM with a
! bias correction algorithm is referred to as value-added NESDIS_LandEM or AlandEM.
!
! program history log:
!
! input argument list:
!
! frequency - frequency in GHz
! theta - local zenith angle in degree
! snow_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:
!
! esv_3w and esh_3w - initial emissivity discriminator at two polarizations
! at three AMSU window channels computed using NESDIS_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
!
!------------------------------------------------------------------------------------------------------------
integer :: nw_ind
parameter(nw_ind=3)
real(fp) theta, frequency, freq,snow_depth, ts, em_vector(2)
real(fp) esv,esh,esh0,esv0,theta0
integer snow_type,ich
real(fp) freq_3w(nw_ind),esh_3w(nw_ind),esv_3w(nw_ind)
complex(fp) eair
save freq_3w
freq_3w = (/31.4_fp,89.0_fp,150.0_fp/)
eair = cmplx(one,-zero,fp)
snow_type = -999
call NESDIS_LandEM(theta, frequency,0.0_fp,0.0_fp,ts,ts,0.0_fp,9,13,snow_depth,esh0,esv0)
theta0 = theta
do ich = 1, nw_ind
freq =freq_3w(ich)
theta = theta0
call NESDIS_LandEM(theta, freq,0.0_fp,0.0_fp,ts,ts,0.0_fp,9,13,snow_depth,esh,esv)
esv_3w(ich) = esv
esh_3w(ich) = esh
end do
call ems_adjust(theta,frequency,snow_depth,ts,esv_3w,esh_3w,em_vector,snow_type)
return
end subroutine AMSU_ALandEM_Snow
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:
!
! input argument list:
!
! frequency - frequency in GHz
! theta - local zenith angle in degree
! 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
!
!------------------------------------------------------------------------------------------------------------
integer,parameter:: nch=10,nw_3=3
integer,parameter:: ncoe=6
real(fp),parameter :: earthrad = 6374._fp, satheight = 833.4_fp
integer :: snow_type,ich,k
real(fp) :: theta,frequency,depth,ts,esv_3w(*),esh_3w(*)
real(fp) :: discriminator(5),emmod(nw_3),dem(nw_3)
real(fp) :: emissivity,em_vector(2)
real(Double) :: dem_coe(nw_3,0:ncoe-1),sinthetas,costhetas,deg2rad
save dem_coe
dem_coe(1,0:ncoe-1) = (/ 2.306844e+000_Double, -7.287718e-003_Double, &
-6.433248e-004_Double, 1.664216e-005_Double, &
4.766508e-007_Double, -1.754184e+000_Double/)
dem_coe(2,0:ncoe-1) = (/ 3.152527e+000_Double, -1.823670e-002_Double, &
-9.535361e-004_Double, 3.675516e-005_Double, &
9.609477e-007_Double, -1.113725e+000_Double/)
dem_coe(3,0:ncoe-1) = (/ 3.492495e+000_Double, -2.184545e-002_Double, &
6.536696e-005_Double, 4.464352e-005_Double, &
-6.305717e-008_Double, -1.221087e+000_Double/)
!
deg2rad = 3.14159_fp*pi/180.0_fp
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_fp
discriminator(2) = emmod(1)
! Missing value at 50.3 GHz
discriminator(3) = -999.9_fp
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
END MODULE NESDIS_AMSU_SnowEM_Module