| blob | 86af2696849aca7442b04287b8f480e828909591 |
1 subroutine snwem_amsu(theta,frequency,snow_depth,ts,tba,tbb,esh,esv)
3 !$$$ subprogram documentation block
4 ! . . . .
5 ! subprogram: noaa/nesdis emissivity model over snow/ice for AMSU-A/B
6 !
7 ! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
8 !
9 ! abstract: noaa/nesdis emissivity model to compute microwave emissivity over
10 ! snow for AMSU-A/B. The processing varies according to input parameters
11 ! Option 1 : AMSU-A & B window channels of brightness temperatures (Tb)
12 ! and surface temperature (Ts) are available
13 ! Option 2 : AMSU-A window channels of Tb and Ts are available
14 ! Option 3 : AMSU-A & B window channels of Tb are available
15 ! Option 4 : AMSU-A window channels of Tb are available
16 ! Option 5 : AMSU-B window channels of Tb and Ts are available
17 ! Option 6 : AMSU-B window channels of Tb are available
18 ! Option 7 : snow depth and Ts are available
19 !
20 ! references:
21 ! Yan, B., F. Weng and K.Okamoto,2004:
22 ! "A microwave snow emissivity model, submitted to TGRS
23 !
24 ! version: 3
25 !
26 ! program history log:
27 ! beta : November 28, 2000
28 !
29 ! version 2.0: June 18, 2003.
30 !
31 ! Version 2.0 enhances the capability/performance of beta version of
32 ! land emissivity model (LandEM) over snow conditions. Two new subroutines
33 ! (i.e., snowem_tb and six_indices) are added as replacements of the
34 ! previous snow emissivity. If AMSU measurements are not available, the
35 ! results are the same as these in beta version. The new snow emissivity
36 ! model is empirically derived from satellite retrievals and ground-based
37 ! measurements.
38 !
39 ! version 3.0: August 18, 2003.
40 !
41 ! Version 3.0 is an extended version of LandEM 2.0 over snow conditions.
42 ! It covers seven different options (see below for details) for LandEM
43 ! inputs over snow conditions. When All or limited AMSU measurements are
44 ! available, one of the subroutines sem_ABTs, sem_ATs, sem_AB, sem_amsua,
45 ! sem_BTs and sem_amsub, which are empirically derived from satellite
46 ! retrievals and ground-based measurements, are called to simuate snow
47 ! emissivity; when no AMSU measurements are avalaiable, the subroutine
48 ! ALandEM_Snow is called where the results over snow conditions in beta
49 ! version are adjusted with a bias correction that is obtained using a
50 ! statistical algorithm. Thus, LandEM 3.0 significantly enhances the
51 ! flexibility/performance of LandEM 2.0 in smulating emissivity over snow
52 ! conditions.
53 !
54 !
55 ! 2004-07-26 okamoto - modified the version 3.0 for use in gsi
56 ! 2004-12-13 kokron - the 'depth' variable was declared then used before being given
57 ! a value. changed to use snow_depth in all calculations
58 !
59 !
60 ! input argument list:
61 ! theta - local zenith angle in radian
62 ! frequency - frequency in GHz
63 ! t_skin - scattering layer temperature (K) (gdas)
64 ! snow_depth - scatter medium depth (mm) (gdas)
65 ! tba[1] ~ tba[4] - Tb at four AMSU-A window channels
66 ! tba[1] : 23.8 GHz
67 ! tba[2] : 31.4 GHz
68 ! tba[3] : 50.3 GHz
69 ! tba[4] : 89 GHz
70 ! tbb[1] ~ tbb[2] - Tb at two AMSU-B window channels:
71 ! tbb[1] : 89 GHz
72 ! tbb[2] : 150 GHz
73 ! When tba[ ] or tbb[ ] = -999.9: a missing value (no available data)
74 !
75 ! output argument list:
76 ! esv - emissivity at vertical polarization
77 ! esh - emissivity at horizontal polarization
78 ! snow_type - snow type (not output here)
79 ! 1 : Wet Snow
80 ! 2 : Grass_after_Snow
81 ! 3 : RS_Snow (A)
82 ! 4 : Powder Snow
83 ! 5 : RS_Snow (B)
84 ! 6 : RS_Snow (C)
85 ! 7 : RS_Snow (D)
86 ! 8 : Thin Crust Snow
87 ! 9 : RS_Snow (E)
88 ! 10: Bottom Crust Snow (A)
89 ! 11: Shallow Snow
90 ! 12: Deep Snow
91 ! 13: Crust Snow
92 ! 14: Medium Snow
93 ! 15: Bottom Crust Snow (B)
94 ! 16: Thick Crust Snow
95 ! 999: AMSU measurements are not available or over non-snow conditions
96 ! important internal variables/parameters:
97 !
98 ! INDATA[1] ~ INDATA[7] - seven options calling different subroutines
99 ! INDATA(1) = ABTs : call sem_ABTs
100 ! INDATA(2) = ATs : call sem_ATs
101 ! INDATA(3) = AMSUAB : call sem_AB
102 ! INDATA(4) = AMSUA : call sem_amsua
103 ! INDATA(5) = BTs : call sem_BTs
104 ! INDATA(6) = AMSUB : call sem_amsub
105 ! INDATA(7) = MODL : call ALandEM_Snow
106 ! input_type - specific option index
107 ! tb[1] ~ tb[5] - Tb at five AMSU-A & B window channels:
108 ! tb[1] = tba[1]
109 ! tb[2] = tba[2]
110 ! tb[3] = tba[3]
111 ! tb[4] = tba[4]
112 ! tb[5] = tbb[2]
113 !
114 ! remarks:
115 !
116 ! Questions/comments: Please send to Fuzhong.Weng@noaa.gov or Banghua.Yan@noaa.gov
117 !
118 ! attributes:
119 ! language: f90
120 ! machine: ibm rs/6000 sp
121 !
122 !$$$
124 ! use kinds, only: r_kind,i_kind
125 implicit none
127 integer(i_kind) :: nch,nwcha,nwchb,nwch,nalg
128 Parameter(nch = 10, nwcha = 4, nwchb = 2, nwch = 5,nalg = 7)
129 real(r_kind) :: theta,frequency,ts,snow_depth
130 real(r_kind) :: em_vector(2),esv,esh
131 real(r_kind) :: tb(nwch),tba(nwcha),tbb(nwchb)
132 logical :: INDATA(nalg),AMSUAB,AMSUA,AMSUB,ABTs,ATs,BTs,MODL
133 integer(i_kind) :: snow_type,input_type,i,ich,np,k
135 Equivalence(INDATA(1), ABTs)
136 Equivalence(INDATA(2), ATs)
137 Equivalence(INDATA(3), AMSUAB)
138 Equivalence(INDATA(4), AMSUA)
139 Equivalence(INDATA(5), BTs)
140 Equivalence(INDATA(6), AMSUB)
141 Equivalence(INDATA(7), MODL)
143 ! Initialization
144 call em_initialization(frequency,em_vector)
145 snow_type = -999
146 input_type = -999
147 do k = 1, nalg
148 INDATA(k) = .TRUE.
149 end do
151 ! Read AMSU & Ts data and set available option
152 ! Get five AMSU-A/B window measurements
153 tb(1) = tba(1); tb(2) = tba(2); tb(3) = tba(3); tb(4) = tba(4); tb(5) = tbb(2)
155 ! Check available data
156 if((ts <= 100.0_r_kind) .or. (ts >= 320.0_r_kind) ) then
157 ABTs = .false.; ATs = .false.; BTs = .false.; MODL = .false.
158 end if
159 do i=1,nwcha
160 if((tba(i) <= 100.0_r_kind) .or. (tba(i) >= 320.0_r_kind) ) then
161 ABTs = .false.; ATs = .false.; AMSUAB = .false.; AMSUA = .false.
162 exit
163 end if
164 end do
165 do i=1,nwchb
166 if((tbb(i) <= 100.0_r_kind) .or. (tbb(i) >= 320.0_r_kind) ) then
167 ABTs = .false.; AMSUAB = .false.; BTs = .false.; AMSUB = .false.
168 exit
169 end if
170 end do
171 if((snow_depth < 0.0_r_kind) .or. (snow_depth >= 3000.0_r_kind)) MODL = .false.
172 if((frequency >= 80._r_kind) .and. (BTs)) then
173 ATs = .false.; AMSUAB = .false.
174 end if
176 ! Check input type and call a specific Option/subroutine
177 do np = 1, nalg
178 if (INDATA(np)) then
179 input_type = np
180 exit
181 end if
182 end do
183 ! 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
185 GET_option: SELECT CASE (input_type)
186 CASE (1)
187 call sem_ABTs(theta,frequency,tb,ts,snow_type,em_vector)
188 CASE (2)
189 call sem_ATs(theta,frequency,tba,ts,snow_type,em_vector)
190 CASE (3)
191 call sem_AB(theta,frequency,tb,snow_type,em_vector)
192 CASE (4)
193 call sem_amsua(theta,frequency,tba,snow_type,em_vector)
194 CASE(5)
195 call sem_BTs(theta,frequency,tbb,ts,snow_type,em_vector)
196 CASE(6)
197 call sem_amsub(theta,frequency,tbb,snow_type,em_vector)
198 CASE(7)
199 call ALandEM_Snow(theta,frequency,snow_depth,ts,snow_type,em_vector)
200 END SELECT GET_option
202 esv = em_vector(1)
203 esh = em_vector(2)
205 return
206 end subroutine snwem_amsu
209 !---------------------------------------------------------------------!
210 subroutine em_initialization(frequency,em_vector)
212 !$$$ subprogram documentation block
213 !
214 ! subprogram: AMSU-A/B snow emissivity initialization
215 !
216 ! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
217 !
218 ! abstract: AMSU-A/B snow emissivity initialization
219 !
220 ! program history log:
221 !
222 ! input argument list:
223 !
224 ! frequency - frequency in GHz
225 !
226 ! output argument list:
227 !
228 ! em_vector[1] and [2] - initial emissivity at two polarizations.
229 !
230 ! important internal variables:
231 !
232 ! freq[1~10] - ten frequencies for sixteen snow types of emissivity
233 ! em[1~16,*] - sixteen snow emissivity spectra
234 ! snow_type - snow type
235 ! where it is initialized to as the type 4,i.e, Powder Snow
236 !
237 ! remarks:
238 !
239 ! attributes:
240 ! language: f90
241 ! machine: ibm rs/6000 sp
242 !
243 !$$$
245 ! use kinds, only: r_kind,i_kind
246 ! use constants, only: one
247 implicit none
249 integer(i_kind) :: nch,ncand
250 Parameter(nch = 10,ncand=16)
251 real(r_kind) :: frequency,em_vector(*),freq(nch)
252 real(r_kind) :: em(ncand,nch)
253 real(r_kind) :: kratio, bconst,emissivity
254 integer(i_kind) :: snow_type,ich,k
255 save em
257 ! Sixteen candidate snow emissivity spectra
258 data (em(1,k),k=1,nch)/0.87_r_kind,0.89_r_kind,0.91_r_kind,0.93_r_kind, &
259 0.94_r_kind,0.94_r_kind,0.94_r_kind,0.93_r_kind,0.92_r_kind,0.90_r_kind/
260 data (em(2,k),k=1,nch)/0.91_r_kind,0.91_r_kind,0.92_r_kind,0.91_r_kind, &
261 0.90_r_kind,0.90_r_kind,0.91_r_kind,0.91_r_kind,0.91_r_kind,0.86_r_kind/
262 data (em(3,k),k=1,nch)/0.90_r_kind,0.89_r_kind,0.88_r_kind,0.87_r_kind, &
263 0.86_r_kind,0.86_r_kind,0.85_r_kind,0.85_r_kind,0.82_r_kind,0.82_r_kind/
264 data (em(4,k),k=1,nch)/0.91_r_kind,0.91_r_kind,0.93_r_kind,0.93_r_kind, &
265 0.93_r_kind,0.93_r_kind,0.89_r_kind,0.88_r_kind,0.79_r_kind,0.79_r_kind/
266 data (em(5,k),k=1,nch)/0.90_r_kind,0.89_r_kind,0.88_r_kind,0.85_r_kind, &
267 0.84_r_kind,0.83_r_kind,0.83_r_kind,0.82_r_kind,0.79_r_kind,0.73_r_kind/
268 data (em(6,k),k=1,nch)/0.90_r_kind,0.89_r_kind,0.86_r_kind,0.82_r_kind, &
269 0.80_r_kind,0.79_r_kind,0.78_r_kind,0.78_r_kind,0.77_r_kind,0.77_r_kind/
270 data (em(7,k),k=1,nch)/0.88_r_kind,0.86_r_kind,0.85_r_kind,0.80_r_kind, &
271 0.78_r_kind,0.77_r_kind,0.77_r_kind,0.76_r_kind,0.72_r_kind,0.72_r_kind/
272 data (em(8,k),k=1,nch)/0.93_r_kind,0.94_r_kind,0.96_r_kind,0.96_r_kind, &
273 0.95_r_kind,0.93_r_kind,0.87_r_kind,0.86_r_kind,0.74_r_kind,0.65_r_kind/
274 data (em(9,k),k=1,nch)/0.87_r_kind,0.86_r_kind,0.84_r_kind,0.80_r_kind, &
275 0.76_r_kind,0.76_r_kind,0.75_r_kind,0.75_r_kind,0.70_r_kind,0.69_r_kind/
276 data (em(10,k),k=1,nch)/0.87_r_kind,0.86_r_kind,0.83_r_kind,0.77_r_kind, &
277 0.73_r_kind,0.68_r_kind,0.66_r_kind,0.66_r_kind,0.68_r_kind,0.67_r_kind/
278 data (em(11,k),k=1,nch)/0.89_r_kind,0.89_r_kind,0.88_r_kind,0.87_r_kind, &
279 0.86_r_kind,0.82_r_kind,0.77_r_kind,0.76_r_kind,0.69_r_kind,0.64_r_kind/
280 data (em(12,k),k=1,nch)/0.88_r_kind,0.87_r_kind,0.86_r_kind,0.83_r_kind, &
281 0.81_r_kind,0.77_r_kind,0.74_r_kind,0.73_r_kind,0.69_r_kind,0.64_r_kind/
282 data (em(13,k),k=1,nch)/0.86_r_kind,0.86_r_kind,0.86_r_kind,0.85_r_kind, &
283 0.82_r_kind,0.78_r_kind,0.69_r_kind,0.68_r_kind,0.51_r_kind,0.47_r_kind/
284 data (em(14,k),k=1,nch)/0.89_r_kind,0.88_r_kind,0.87_r_kind,0.83_r_kind, &
285 0.80_r_kind,0.75_r_kind,0.70_r_kind,0.70_r_kind,0.64_r_kind,0.60_r_kind/
286 data (em(15,k),k=1,nch)/0.91_r_kind,0.92_r_kind,0.93_r_kind,0.88_r_kind, &
287 0.84_r_kind,0.76_r_kind,0.66_r_kind,0.64_r_kind,0.48_r_kind,0.44_r_kind/
288 data (em(16,k),k=1,nch)/0.94_r_kind,0.95_r_kind,0.97_r_kind,0.91_r_kind, &
289 0.86_r_kind,0.74_r_kind,0.63_r_kind,0.63_r_kind,0.50_r_kind,0.45_r_kind/
290 data freq/4.9_r_kind,6.93_r_kind,10.65_r_kind,18.7_r_kind,23.8_r_kind, &
291 31.4_r_kind, 50.3_r_kind,52.5_r_kind, 89.0_r_kind, 150._r_kind/
293 ! Initialization for emissivity at certain frequency
294 ! In case of no any inputs available for various options
295 ! A constant snow type & snow emissivity spectrum is assumed
296 ! (e.g., powder) snow_type = 4
298 ! Specify snow emissivity at required frequency
299 do ich = 2, nch
300 if(frequency < freq(1)) exit
301 if(frequency >= freq(nch)) exit
302 if(frequency < freq(ich)) then
303 emissivity = em(4,ich-1) + (em(4,ich) - em(4,ich-1)) &
304 *(frequency - freq(ich-1))/(freq(ich) - freq(ich-1))
305 exit
306 end if
307 end do
309 ! Extrapolate to lower frequencies than 4.9GHz
310 if (frequency <= freq(1)) then
311 kratio = (em(4,2) - em(4,1))/(freq(2) - freq(1))
312 bconst = em(4,1) - kratio*freq(1)
313 emissivity = kratio*frequency + bconst
314 if(emissivity > one) emissivity = one
315 if(emissivity <= 0.8_r_kind) emissivity = 0.8_r_kind
316 end if
319 ! Assume emissivity = constant at frequencies >= 150 GHz
320 if (frequency >= freq(nch)) emissivity = em(4,nch)
321 em_vector(1) = emissivity
322 em_vector(2) = emissivity
324 return
325 end subroutine em_initialization
329 !---------------------------------------------------------------------!
330 subroutine em_interpolate(frequency,discriminator,emissivity,snow_type)
332 !$$$ subprogram documentation block
333 !
334 ! subprogram: determine snow_type and calculate emissivity
335 !
336 ! prgmmr:Banghua Yan org: nesdis date: 2003-08-18
337 !
338 ! abstract: 1. Find one snow emissivity spectrum to mimic the emission
339 ! property of the realistic snow condition using a set of
340 ! discrminators
341 ! 2. Interpolate/extrapolate emissivity at a required frequency
342 !
343 ! program history log:
344 !
345 ! input argument list:
346 !
347 ! frequency - frequency in GHz
348 ! discriminators - emissivity discriminators at five AMSU-A & B window
349 ! channels
350 ! discriminator[1] : emissivity discriminator at 23.8 GHz
351 ! discriminator[2] : emissivity discriminator at 31.4 GHz
352 ! discriminator[3] : emissivity discriminator at 50.3 GHz
353 ! discriminator[4] : emissivity discriminator at 89 GHz
354 ! discriminator[5] : emissivity discriminator at 150 GHz
355 !
356 ! Note: discriminator(1) and discriminator(3) are missing value in
357 ! 'AMSU-B & Ts','AMUS-B' and 'MODL' options., which are defined to as -999.9,
358 ! output argument list:
359 !
360 ! em_vector[1] and [2] - emissivity at two polarizations.
361 ! snow_type - snow type
362 !
363 ! important internal variables:
364 !
365 ! freq[1 ~ 10] - ten frequencies for sixteen snow types of emissivity
366 ! em[1~16,*] - sixteen snow emissivity spectra
367 !
368 ! remarks:
369 !
370 ! attributes:
371 ! language: f90
372 ! machine: ibm rs/6000 sp
373 !
374 !$$$
376 ! use kinds, only: r_kind,i_kind
377 ! use constants, only: zero, one
378 implicit none
380 integer(i_kind),parameter:: ncand = 16,nch =10
381 integer(i_kind):: ich,ichmin,ichmax,i,j,k,s,snow_type
382 real(r_kind) :: dem,demmin0
383 real(r_kind) :: em(ncand,nch)
384 real(r_kind) :: frequency,freq(nch),emissivity,discriminator(*),emis(nch)
385 real(r_kind) :: cor_factor,adjust_check,kratio, bconst
386 data freq/4.9_r_kind, 6.93_r_kind, 10.65_r_kind, 18.7_r_kind,&
387 23.8_r_kind, 31.4_r_kind, 50.3_r_kind, 52.5_r_kind, &
388 89.0_r_kind, 150._r_kind/
390 ! Sixteen candidate snow emissivity spectra
391 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,&
392 0.94_r_kind,0.94_r_kind,0.93_r_kind,0.92_r_kind,0.90_r_kind/
393 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,&
394 0.90_r_kind,0.91_r_kind,0.91_r_kind,0.91_r_kind,0.86_r_kind/
395 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,&
396 0.86_r_kind,0.85_r_kind,0.85_r_kind,0.82_r_kind,0.82_r_kind/
397 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,&
398 0.93_r_kind,0.89_r_kind,0.88_r_kind,0.79_r_kind,0.79_r_kind/
399 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,&
400 0.83_r_kind,0.83_r_kind,0.82_r_kind,0.79_r_kind,0.73_r_kind/
401 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,&
402 0.79_r_kind,0.78_r_kind,0.78_r_kind,0.77_r_kind,0.77_r_kind/
403 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,&
404 0.77_r_kind,0.77_r_kind,0.76_r_kind,0.72_r_kind,0.72_r_kind/
405 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,&
406 0.93_r_kind,0.87_r_kind,0.86_r_kind,0.74_r_kind,0.65_r_kind/
407 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,&
408 0.76_r_kind,0.75_r_kind,0.75_r_kind,0.70_r_kind,0.69_r_kind/
409 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,&
410 0.68_r_kind,0.66_r_kind,0.66_r_kind,0.68_r_kind,0.67_r_kind/
411 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,&
412 0.82_r_kind,0.77_r_kind,0.76_r_kind,0.69_r_kind,0.64_r_kind/
413 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,&
414 0.77_r_kind,0.74_r_kind,0.73_r_kind,0.69_r_kind,0.64_r_kind/
415 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,&
416 0.78_r_kind,0.69_r_kind,0.68_r_kind,0.51_r_kind,0.47_r_kind/
417 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,&
418 0.75_r_kind,0.70_r_kind,0.70_r_kind,0.64_r_kind,0.60_r_kind/
419 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,&
420 0.76_r_kind,0.66_r_kind,0.64_r_kind,0.48_r_kind,0.44_r_kind/
421 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,&
422 0.74_r_kind,0.63_r_kind,0.63_r_kind,0.50_r_kind,0.45_r_kind/
423 save em
425 ! Adjust unreasonable discriminator
426 if (discriminator(4) > discriminator(2)) &
427 discriminator(4) = discriminator(2) +(discriminator(5) - discriminator(2))* &
428 (150.0_r_kind - 89.0_r_kind)/(150.0_r_kind - 31.4_r_kind)
429 if ( (discriminator(3) /= -999.9_r_kind) .and. &
430 ( ((discriminator(3)-0.01_r_kind) > discriminator(2)) .or. &
431 ((discriminator(3)-0.01_r_kind) < discriminator(4))) ) &
432 discriminator(3) = discriminator(2) + &
433 (discriminator(4) - discriminator(2))*(89.0_r_kind - 50.3_r_kind) &
434 / (89.0_r_kind - 31.4_r_kind)
436 ! Find a snow emissivity spectrum
437 if(snow_type .eq. -999) then
438 demmin0 = 10.0_r_kind
439 do k = 1, ncand
440 dem = zero
441 ichmin = 1
442 ichmax = 3
443 if(discriminator(1) == -999.9_r_kind) then
444 ichmin = 2
445 ichmax = 2
446 end if
447 do ich = ichmin,ichmax
448 dem = dem + abs(discriminator(ich) - em(k,ich+4))
449 end do
450 do ich = 4,5
451 dem = dem + abs(discriminator(ich) - em(k,ich+5))
452 end do
453 if (dem < demmin0) then
454 demmin0 = dem
455 snow_type = k
456 end if
457 end do
458 end if
460 ! Shift snow emissivity according to discriminator at 31.4 GHz
461 cor_factor = discriminator(2) - em(snow_type,6)
462 do ich = 1, nch
463 emis(ich) = em(snow_type,ich) + cor_factor
464 if(emis(ich) .gt. one) emis(ich) = one
465 if(emis(ich) .lt. 0.3_r_kind) emis(ich) = 0.3_r_kind
466 end do
468 ! Emisivity data quality control
469 adjust_check = zero
470 do ich = 5, 9
471 if (ich .le. 7) then
472 if (discriminator(ich - 4) .ne. -999.9_r_kind) &
473 adjust_check = adjust_check + abs(emis(ich) - discriminator(ich - 4))
474 else
475 if (discriminator(ich - 4) .ne. -999.9_r_kind) &
476 adjust_check = adjust_check + abs(emis(ich+1) - discriminator(ich - 4))
477 end if
478 end do
480 if (adjust_check >= 0.04_r_kind) then
481 if (discriminator(1) /= -999.9_r_kind) then
482 if (discriminator(1) < emis(4)) then
483 emis(5) = emis(4) + &
484 (31.4_r_kind - 23.8_r_kind) * &
485 (discriminator(2) - emis(4))/(31.4_r_kind - 18.7_r_kind)
486 else
487 emis(5) = discriminator(1)
488 end if
489 end if
490 emis(6) = discriminator(2)
491 if (discriminator(3) /= -999.9_r_kind) then
492 emis(7) = discriminator(3)
493 else
494 ! In case of missing the emissivity discriminator at 50.3 GHz
495 emis(7) = emis(6) + (89.0_r_kind - 50.3_r_kind) * &
496 (discriminator(4) - emis(6))/(89.0_r_kind - 31.4_r_kind)
497 end if
498 emis(8) = emis(7)
499 emis(9) = discriminator(4)
500 emis(10) = discriminator(5)
501 end if
503 ! Estimate snow emissivity at a required frequency
504 do i = 2, nch
505 if(frequency < freq(1)) exit
506 if(frequency >= freq(nch)) exit
507 if(frequency < freq(i)) then
508 emissivity = emis(i-1) + (emis(i) - emis(i-1))*(frequency - freq(i-1)) &
509 /(freq(i) - freq(i-1))
510 exit
511 end if
512 end do
514 ! Extrapolate to lower frequencies than 4.9GHz
515 if (frequency <= freq(1)) then
516 kratio = (emis(2) - emis(1))/(freq(2) - freq(1))
517 bconst = emis(1) - kratio*freq(1)
518 emissivity = kratio*frequency + bconst
519 if(emissivity > one) emissivity = one
520 if(emissivity <= 0.8_r_kind) emissivity = 0.8_r_kind
521 end if
523 ! Assume emissivity = constant at frequencies >= 150 GHz
524 if (frequency >= freq(nch)) emissivity = emis(nch)
526 return
527 end subroutine em_interpolate
530 !---------------------------------------------------------------------!
531 subroutine sem_ABTs(theta,frequency,tb,ts,snow_type,em_vector)
533 !$$$ subprogram documentation block
534 !
535 ! subprogram:
536 !
537 ! prgmmr:Banghua Yan org: nesdis date: 2003-08-18
538 !
539 ! abstract:
540 ! Calculate the emissivity discriminators and interpolate/extrapolate
541 ! emissivity at a required frequency with respect to secenery ABTs
542 !
543 ! program history log:
544 !
545 ! input argument list:
546 !
547 ! frequency - frequency in GHz
548 ! theta - local zenith angle (not used here)
549 ! tb[1] ~ tb[5] - brightness temperature at five AMSU window channels:
550 ! tb[1] : 23.8 GHz
551 ! tb[2] : 31.4 GHz
552 ! tb[3] : 50.3 GHz
553 ! tb[4] : 89.0 GHz
554 ! tb[5] : 150 GHz
555 !
556 ! output argument list:
557 !
558 ! em_vector[1] and [2] - emissivity at two polarizations.
559 ! set esv = esh here and will be updated
560 ! snow_type - snow type
561 !
562 ! important internal variables:
563 !
564 ! nind - number of threshold in decision trees
565 ! to identify each snow type ( = 6)
566 ! em(1~16,*) - sixteen snow emissivity spectra
567 ! DI_coe - coefficients to generate six discriminators to describe
568 ! the overall emissivity variability within a wider frequency range
569 ! threshold - thresholds in decision trees to identify snow types
570 ! index_in - six indices to discriminate snow type
571 !
572 ! remarks:
573 !
574 ! attributes:
575 ! language: f90
576 ! machine: ibm rs/6000 sp
577 !
578 !$$$
580 ! use kinds, only: r_kind,i_kind
581 implicit none
583 integer(i_kind),parameter:: ncand = 16,nch =10,nthresh=38
584 integer(i_kind),parameter:: nind=6,ncoe=8,nLIcoe=6,nHIcoe=12
585 integer(i_kind):: ich,i,j,k,num,npass,snow_type,md0,md1,nmodel(ncand-1)
586 real(r_kind) :: theta,frequency,tb150,LI,HI,DS1,DS2,DS3
587 real(r_kind) :: em(ncand,nch), em_vector(*)
588 real(r_kind) :: tb(*),freq(nch),DTB(nind-1),DI(nind-1), &
589 DI_coe(nind-1,0:ncoe-1),threshold(nthresh,nind), &
590 index_in(nind),threshold0(nind)
591 real(r_kind) :: LI_coe(0:nLIcoe-1),HI_coe(0:nHIcoe-1)
592 real(r_kind) :: ts,emissivity
593 real(r_kind) :: discriminator(5)
594 logical:: pick_status,tindex(nind)
595 save em,threshold,DI_coe,LI_coe, HI_coe,nmodel,freq
597 data freq/4.9_r_kind,6.93_r_kind,10.65_r_kind,18.7_r_kind,23.8_r_kind, &
598 31.4_r_kind, 50.3_r_kind,52.5_r_kind, 89.0_r_kind, 150._r_kind/
599 data nmodel/5,10,13,16,18,24,30,31,32,33,34,35,36,37,38/
601 ! Fitting coefficients for five discriminators
602 data (DI_coe(1,k),k=0,ncoe-1)/ &
603 3.285557e-002_r_kind, 2.677179e-005_r_kind, &
604 4.553101e-003_r_kind, 5.639352e-005_r_kind, &
605 -1.825188e-004_r_kind, 1.636145e-004_r_kind, &
606 1.680881e-005_r_kind, -1.708405e-004_r_kind/
607 data (DI_coe(2,k),k=0,ncoe-1)/ &
608 -4.275539e-002_r_kind, -2.541453e-005_r_kind, &
609 4.154796e-004_r_kind, 1.703443e-004_r_kind, &
610 4.350142e-003_r_kind, 2.452873e-004_r_kind, &
611 -4.748506e-003_r_kind, 2.293836e-004_r_kind/
612 data (DI_coe(3,k),k=0,ncoe-1)/ &
613 -1.870173e-001_r_kind, -1.061678e-004_r_kind, &
614 2.364055e-004_r_kind, -2.834876e-005_r_kind, &
615 4.899651e-003_r_kind, -3.418847e-004_r_kind, &
616 -2.312224e-004_r_kind, 9.498600e-004_r_kind/
617 data (DI_coe(4,k),k=0,ncoe-1)/ &
618 -2.076519e-001_r_kind, 8.475901e-004_r_kind, &
619 -2.072679e-003_r_kind, -2.064717e-003_r_kind, &
620 2.600452e-003_r_kind, 2.503923e-003_r_kind, &
621 5.179711e-004_r_kind, 4.667157e-005_r_kind/
622 data (DI_coe(5,k),k=0,ncoe-1)/ &
623 -1.442609e-001_r_kind, -8.075003e-005_r_kind, &
624 -1.790933e-004_r_kind, -1.986887e-004_r_kind, &
625 5.495115e-004_r_kind, -5.871732e-004_r_kind, &
626 4.517280e-003_r_kind, 7.204695e-004_r_kind/
628 ! Fitting coefficients for emissivity index at 31.4 GHz
629 data LI_coe/ &
630 7.963632e-001_r_kind, 7.215580e-003_r_kind, &
631 -2.015921e-005_r_kind, -1.508286e-003_r_kind, &
632 1.731405e-005_r_kind, -4.105358e-003_r_kind/
634 ! Fitting coefficients for emissivity index at 150 GHz
635 data HI_coe/ &
636 1.012160e+000_r_kind, 6.100397e-003_r_kind, &
637 -1.774347e-005_r_kind, -4.028211e-003_r_kind, &
638 1.224470e-005_r_kind, 2.345612e-003_r_kind, &
639 -5.376814e-006_r_kind, -2.795332e-003_r_kind, &
640 8.072756e-006_r_kind, 3.529615e-003_r_kind, &
641 1.955293e-006_r_kind, -4.942230e-003_r_kind/
643 ! Six thresholds for sixteen candidate snow types
644 ! Note: some snow type contains several possible
645 ! selections for six thresholds
647 !1 Wet Snow
648 data (threshold(1,k),k=1,6)/0.88_r_kind,0.86_r_kind,-999.9_r_kind,&
649 0.01_r_kind,0.01_r_kind,200._r_kind/
650 data (threshold(2,k),k=1,6)/0.88_r_kind,0.85_r_kind,-999.9_r_kind,&
651 0.06_r_kind,0.10_r_kind,200._r_kind/
652 data (threshold(3,k),k=1,6)/0.88_r_kind,0.83_r_kind,-0.02_r_kind,&
653 0.12_r_kind,0.16_r_kind,204._r_kind/
654 data (threshold(4,k),k=1,6)/0.90_r_kind,0.89_r_kind,-999.9_r_kind,&
655 -999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
656 data (threshold(5,k),k=1,6)/0.92_r_kind,0.85_r_kind,-999.9_r_kind,&
657 -999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
659 !2 Grass_after_Snow
660 data (threshold(6,k),k=1,6)/0.84_r_kind,0.83_r_kind,-999.9_r_kind,&
661 0.08_r_kind,0.10_r_kind,195._r_kind/
662 data (threshold(7,k),k=1,6)/0.85_r_kind,0.85_r_kind,-999.9_r_kind,&
663 0.10_r_kind,-999.9_r_kind,190._r_kind/
664 data (threshold(8,k),k=1,6)/0.86_r_kind,0.81_r_kind,-999.9_r_kind,&
665 0.12_r_kind,-999.9_r_kind,200._r_kind/
666 data (threshold(9,k),k=1,6)/0.86_r_kind,0.81_r_kind,0.0_r_kind,&
667 0.12_r_kind,-999.9_r_kind,189._r_kind/
668 data (threshold(10,k),k=1,6)/0.90_r_kind,0.81_r_kind,-999.9_r_kind,&
669 -999.9_r_kind,-999.9_r_kind,195._r_kind/
671 !3 RS_Snow (A)
672 data (threshold(11,k),k=1,6)/0.80_r_kind,0.76_r_kind,-999.9_r_kind,&
673 0.05_r_kind,-999.9_r_kind,185._r_kind/
674 data (threshold(12,k),k=1,6)/0.82_r_kind,0.78_r_kind,-999.9_r_kind,&
675 -999.9_r_kind,0.25_r_kind,180._r_kind/
676 data (threshold(13,k),k=1,6)/0.90_r_kind,0.76_r_kind,-999.9_r_kind,&
677 -999.9_r_kind,-999.9_r_kind,180._r_kind/
679 !4 Powder Snow
680 data (threshold(14,k),k=1,6)/0.89_r_kind,0.73_r_kind,-999.9_r_kind,&
681 0.20_r_kind,-999.9_r_kind,-999.9_r_kind/
682 data (threshold(15,k),k=1,6)/0.89_r_kind,0.75_r_kind,-999.9_r_kind,&
683 -999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
684 data (threshold(16,k),k=1,6)/0.93_r_kind,0.72_r_kind,-999.9_r_kind,&
685 -999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
687 !5 RS_Snow (B)
688 data (threshold(17,k),k=1,6)/0.82_r_kind,0.70_r_kind,-999.9_r_kind,&
689 0.20_r_kind,-999.9_r_kind,160._r_kind/
690 data (threshold(18,k),k=1,6)/0.83_r_kind,0.70_r_kind,-999.9_r_kind,&
691 -999.9_r_kind,-999.9_r_kind,160._r_kind/
693 !6 RS_Snow (C)
694 data (threshold(19,k),k=1,6)/0.75_r_kind,0.76_r_kind,-999.9_r_kind,&
695 0.08_r_kind,-999.9_r_kind,172._r_kind/
696 data (threshold(20,k),k=1,6)/0.77_r_kind,0.72_r_kind,-999.9_r_kind,&
697 0.12_r_kind,0.15_r_kind,175._r_kind/
698 data (threshold(21,k),k=1,6)/0.78_r_kind,0.74_r_kind,-999.9_r_kind,&
699 -999.9_r_kind,0.20_r_kind,172._r_kind/
700 data (threshold(22,k),k=1,6)/0.80_r_kind,0.77_r_kind,-999.9_r_kind,&
701 -999.9_r_kind,-999.9_r_kind,170._r_kind/
702 data (threshold(23,k),k=1,6)/0.82_r_kind,-999.9_r_kind,-999.9_r_kind,&
703 0.15_r_kind,0.22_r_kind,170._r_kind/
704 data (threshold(24,k),k=1,6)/0.82_r_kind,0.73_r_kind,-999.9_r_kind,&
705 -999.9_r_kind,-999.9_r_kind,170._r_kind/
707 !7 RS_Snow (D)
708 data (threshold(25,k),k=1,6)/0.75_r_kind,0.70_r_kind,-999.9_r_kind,&
709 0.15_r_kind,0.25_r_kind,167._r_kind/
710 data (threshold(26,k),k=1,6)/0.77_r_kind,0.76_r_kind,-999.9_r_kind,&
711 -999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
712 data (threshold(27,k),k=1,6)/0.80_r_kind,0.72_r_kind,-999.9_r_kind,&
713 -999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
714 data (threshold(28,k),k=1,6)/0.77_r_kind,0.73_r_kind,-999.9_r_kind,&
715 -999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
717 data (threshold(29,k),k=1,6)/0.81_r_kind,0.71_r_kind,-999.9_r_kind,&
718 -999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
719 data (threshold(30,k),k=1,6)/0.82_r_kind,0.69_r_kind,-999.9_r_kind,&
720 -999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
722 !8 Thin Crust Snow
723 data (threshold(31,k),k=1,6)/0.88_r_kind,0.58_r_kind,-999.9_r_kind,&
724 -999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
726 !9 RS_Snow (E)
727 data (threshold(32,k),k=1,6)/0.73_r_kind,0.67_r_kind,-999.9_r_kind,&
728 -999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
730 !10 Bottom Crust Snow (A)
731 data (threshold(33,k),k=1,6)/0.83_r_kind,0.66_r_kind,-999.9_r_kind,&
732 -999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
734 !11 Shallow Snow
735 data (threshold(34,k),k=1,6)/0.82_r_kind,0.60_r_kind,-999.9_r_kind,&
736 -999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
738 !12 Deep Snow
739 data (threshold(35,k),k=1,6)/0.77_r_kind,0.60_r_kind,-999.9_r_kind,&
740 -999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
742 !13 Crust Snow
743 data (threshold(36,k),k=1,6)/0.77_r_kind,0.7_r_kind,-999.9_r_kind,&
744 -999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
746 !14 Medium Snow
747 data (threshold(37,k),k=1,6)/-999.9_r_kind,0.55_r_kind,-999.9_r_kind,&
748 -999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
750 !15 Bottom Crust Snow(B)
751 data (threshold(38,k),k=1,6)/0.74_r_kind,-999.9_r_kind,-999.9_r_kind,&
752 -999.9_r_kind,-999.9_r_kind,-999.9_r_kind/
754 !16 Thick Crust Snow
755 ! lowest priority: No constraints
757 ! Sixteen candidate snow emissivity spectra
758 data (em(1,k),k=1,nch)/0.87_r_kind,0.89_r_kind,0.91_r_kind,0.93_r_kind,&
759 0.94_r_kind,0.94_r_kind,0.94_r_kind,0.93_r_kind,0.92_r_kind,0.90_r_kind/
760 data (em(2,k),k=1,nch)/0.91_r_kind,0.91_r_kind,0.92_r_kind,0.91_r_kind,&
761 0.90_r_kind,0.90_r_kind,0.91_r_kind,0.91_r_kind,0.91_r_kind,0.86_r_kind/
762 data (em(3,k),k=1,nch)/0.90_r_kind,0.89_r_kind,0.88_r_kind,0.87_r_kind,&
763 0.86_r_kind,0.86_r_kind,0.85_r_kind,0.85_r_kind,0.82_r_kind,0.82_r_kind/
764 data (em(4,k),k=1,nch)/0.91_r_kind,0.91_r_kind,0.93_r_kind,0.93_r_kind,&
765 0.93_r_kind,0.93_r_kind,0.89_r_kind,0.88_r_kind,0.79_r_kind,0.79_r_kind/
766 data (em(5,k),k=1,nch)/0.90_r_kind,0.89_r_kind,0.88_r_kind,0.85_r_kind,&
767 0.84_r_kind,0.83_r_kind,0.83_r_kind,0.82_r_kind,0.79_r_kind,0.73_r_kind/
768 data (em(6,k),k=1,nch)/0.90_r_kind,0.89_r_kind,0.86_r_kind,0.82_r_kind,&
769 0.80_r_kind,0.79_r_kind,0.78_r_kind,0.78_r_kind,0.77_r_kind,0.77_r_kind/
770 data (em(7,k),k=1,nch)/0.88_r_kind,0.86_r_kind,0.85_r_kind,0.80_r_kind,&
771 0.78_r_kind,0.77_r_kind,0.77_r_kind,0.76_r_kind,0.72_r_kind,0.72_r_kind/
772 data (em(8,k),k=1,nch)/0.93_r_kind,0.94_r_kind,0.96_r_kind,0.96_r_kind,&
773 0.95_r_kind,0.93_r_kind,0.87_r_kind,0.86_r_kind,0.74_r_kind,0.65_r_kind/
774 data (em(9,k),k=1,nch)/0.87_r_kind,0.86_r_kind,0.84_r_kind,0.80_r_kind,&
775 0.76_r_kind,0.76_r_kind,0.75_r_kind,0.75_r_kind,0.70_r_kind,0.69_r_kind/
776 data (em(10,k),k=1,nch)/0.87_r_kind,0.86_r_kind,0.83_r_kind,0.77_r_kind,&
777 .73_r_kind,0.68_r_kind,0.66_r_kind,0.66_r_kind,0.68_r_kind,0.67_r_kind/
778 data (em(11,k),k=1,nch)/0.89_r_kind,0.89_r_kind,0.88_r_kind,0.87_r_kind,&
779 0.86_r_kind,0.82_r_kind,0.77_r_kind,0.76_r_kind,0.69_r_kind,0.64_r_kind/
780 data (em(12,k),k=1,nch)/0.88_r_kind,0.87_r_kind,0.86_r_kind,0.83_r_kind,&
781 0.81_r_kind,0.77_r_kind,0.74_r_kind,0.73_r_kind,0.69_r_kind,0.64_r_kind/
782 data (em(13,k),k=1,nch)/0.86_r_kind,0.86_r_kind,0.86_r_kind,0.85_r_kind,&
783 0.82_r_kind,0.78_r_kind,0.69_r_kind,0.68_r_kind,0.51_r_kind,0.47_r_kind/
784 data (em(14,k),k=1,nch)/0.89_r_kind,0.88_r_kind,0.87_r_kind,0.83_r_kind,&
785 0.80_r_kind,0.75_r_kind,0.70_r_kind,0.70_r_kind,0.64_r_kind,0.60_r_kind/
786 data (em(15,k),k=1,nch)/0.91_r_kind,0.92_r_kind,0.93_r_kind,0.88_r_kind,&
787 0.84_r_kind,0.76_r_kind,0.66_r_kind,0.64_r_kind,0.48_r_kind,0.44_r_kind/
788 data (em(16,k),k=1,nch)/0.94_r_kind,0.95_r_kind,0.97_r_kind,0.91_r_kind,&
789 0.86_r_kind,0.74_r_kind,0.63_r_kind,0.63_r_kind,0.50_r_kind,0.45_r_kind/
791 !*** DEFINE SIX DISCRIMINATORS
793 dtb(1) = tb(1) - tb(2)
794 dtb(2) = tb(2) - tb(4)
795 dtb(3) = tb(2) - tb(5)
796 dtb(4) = tb(3) - tb(5)
797 dtb(5) = tb(4) - tb(5)
798 tb150 = tb(5)
800 LI = LI_coe(0)
801 do i=0,1
802 LI = LI + LI_coe(2*i+1)*tb(i+1) + LI_coe(2*i+2)*tb(i+1)*tb(i+1)
803 end do
804 LI = LI + LI_coe(nLIcoe-1)*ts
806 HI = HI_coe(0)
807 do i=0,4
808 HI = HI + HI_coe(2*i+1)*tb(i+1) + HI_coe(2*i+2)*tb(i+1)*tb(i+1)
809 end do
810 HI = HI + HI_coe(nHIcoe-1)*ts
812 do num=1,nind-1
813 DI(num) = DI_coe(num,0) + DI_coe(num,1)*tb(2)
814 do i=1,5
815 DI(num) = DI(num) + DI_coe(num,1+i)*DTB(i)
816 end do
817 DI(num) = DI(num) + DI_coe(num,ncoe-1)*ts
818 end do
820 !*** DEFINE FIVE INDIES
821 !HI = DI(0) - DI(3)
822 DS1 = DI(1) + DI(2)
823 DS2 = DI(4) + DI(5)
824 DS3 = DS1 + DS2 + DI(3)
826 index_in(1) = LI
827 index_in(2) = HI
828 index_in(3) = DS1
829 index_in(4) = DS2
830 index_in(5) = DS3
831 index_in(6) = tb150
833 !*** IDENTIFY SNOW TYPE
836 ! Initialization
837 md0 = 1
838 snow_type = ncand
839 pick_status = .false.
841 ! Pick one snow type
842 ! Check all possible selections for six thresholds for each snow type
843 do i = 1, ncand - 1
844 md1 = nmodel(i)
845 do j = md0, md1
846 npass = 0
847 do k = 1 , nind
848 threshold0(k) = threshold(j,k)
849 end do
850 CALL six_indices(nind,index_in,threshold0,tindex)
852 ! Corrections
853 if((i == 5) .and. (index_in(2) > 0.75_r_kind)) tindex(2) = .false.
854 if((i == 5) .and. (index_in(4) > 0.20_r_kind) &
855 .and. (index_in(1) > 0.88_r_kind)) tindex(1) = .false.
856 if((i == 10) .and. (index_in(1) <= 0.83_r_kind)) tindex(1) = .true.
857 if((i == 13) .and. (index_in(2) < 0.52_r_kind)) tindex(2) = .true.
858 do k = 1, nind
859 if(.not.tindex(k)) exit
860 npass = npass + 1
861 end do
862 if(npass == nind) exit
863 end do
865 if(npass == nind) then
866 pick_status = .true.
867 snow_type = i
868 end if
869 if(pick_status) exit
870 md0 = md1 + 1
871 end do
873 discriminator(1) = LI + DI(1)
874 discriminator(2) = LI
875 discriminator(3) = DI(4) + HI
876 discriminator(4) = LI - DI(2)
877 discriminator(5) = HI
879 call em_interpolate(frequency,discriminator,emissivity,snow_type)
881 em_vector(1) = emissivity
882 em_vector(2) = emissivity
884 return
885 end subroutine sem_ABTs
888 !---------------------------------------------------------------------!
889 subroutine six_indices(nind,index_in,threshold,tindex)
891 !$$$ subprogram documentation block
892 !
893 ! subprogram:
894 !
895 ! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
896 !
897 ! abstract:
898 !
899 ! program history log:
900 !
901 ! input argument list:
902 !
903 ! nind - Number of threshold in decision trees
904 ! to identify each snow type ( = 6)
905 ! index_in - six indices to discriminate snow type
906 ! threshold - Thresholds in decision trees to identify snow types
907 !
908 ! output argument list:
909 !
910 ! tindex - state vaiable to show surface snow emissivity feature
911 ! tindex[ ] = .T.: snow emissivity feature matches the
912 ! corresponding threshold for certain snow type
913 ! tindex[ ] = .F.: snow emissivity feature doesn't match the
914 ! corresponding threshold for certain snow type
915 !
916 ! remarks:
917 !
918 ! attributes:
919 ! language: f90
920 ! machine: ibm rs/6000 sp
921 !
922 !$$$
924 ! use kinds, only: r_kind,i_kind
925 implicit none
927 integer(i_kind) :: i,nind
928 real(r_kind) :: index_in(*),threshold(*)
929 logical :: tindex(*)
931 do i=1,nind
932 tindex(i) = .false.
933 if (threshold(i) .eq. -999.9_r_kind) then
934 tindex(i) = .true.
935 else
936 if ( (i .le. 2) .or. (i .gt. (nind-1)) ) then
937 if (index_in(i) .ge. threshold(i)) tindex(i) = .true.
938 else
939 if (index_in(i) .le. threshold(i)) tindex(i) = .true.
940 end if
941 end if
942 end do
943 return
945 end subroutine six_indices
948 !---------------------------------------------------------------------!
949 subroutine sem_AB(theta,frequency,tb,snow_type,em_vector)
951 !$$$ subprogram documentation block
952 !
953 ! subprogram:
954 !
955 ! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
956 !
957 ! abstract:
958 ! Calculate the emissivity discriminators and interpolate/extrapolate
959 ! emissivity at required frequency with respect to option AMSUAB
960 !
961 ! program history log:
962 ! 2004-10-28 treadon - correct problem in declared dimensions of array coe
963 !
964 ! input argument list:
965 !
966 ! frequency - frequency in GHz
967 ! theta - local zenith angle (not used here)
968 ! tb[1]~tb[5] - brightness temperature at five AMSU-A & B window channels:
969 ! tb[1] : 23.8 GHz
970 ! tb[2] : 31.4 GHz
971 ! tb[3] : 50.3 GHz
972 ! tb[4] : 89 GHz
973 ! tb[5] : 150 GHz
974 !
975 ! output argument list:
976 !
977 ! em_vector[1] and [2] - emissivity at two polarizations.
978 ! set esv = esh here and will be updated
979 ! snow_type - snow type (reference [2])
980 !
981 ! important internal variables:
982 !
983 ! coe23 - fitting coefficients to estimate discriminator at 23.8 GHz
984 ! coe31 - fitting coefficients to estimate discriminator at 31.4 GHz
985 ! coe50 - fitting coefficients to estimate discriminator at 50.3 GHz
986 ! coe89 - fitting coefficients to estimate discriminator at 89 GHz
987 ! coe150 - fitting coefficients to estimate discriminator at 150 GHz
988 !
989 ! remarks:
990 !
991 ! attributes:
992 ! language: f90
993 ! machine: ibm rs/6000 sp
994 !
995 !$$$
996 ! use kinds, only: r_kind,i_kind
997 implicit none
999 integer(i_kind),parameter:: nch =10,nwch = 5,ncoe = 10
1000 real(r_kind) :: tb(*),theta,frequency
1001 real(r_kind) :: em_vector(*),emissivity,discriminator(nwch)
1002 integer(i_kind) :: i,snow_type,k,ich,nvalid_ch
1003 real(r_kind) :: coe23(0:ncoe),coe31(0:ncoe),coe50(0:ncoe),coe89(0:ncoe),coe150(0:ncoe)
1004 real(r_kind) :: coe(nch*(ncoe+1))
1006 Equivalence (coe(1),coe23)
1007 Equivalence (coe(12),coe31)
1008 Equivalence (coe(23),coe50)
1009 Equivalence (coe(34),coe89)
1010 Equivalence (coe(45),coe150)
1012 ! Fitting Coefficients at 23.8 GHz: Using Tb1 ~ Tb3
1013 data (coe23(k),k=0,6)/&
1014 -1.326040e+000_r_kind, 2.475904e-002_r_kind, &
1015 -5.741361e-005_r_kind, -1.889650e-002_r_kind, &
1016 6.177911e-005_r_kind, 1.451121e-002_r_kind, &
1017 -4.925512e-005_r_kind/
1019 ! Fitting Coefficients at 31.4 GHz: Using Tb1 ~ Tb3
1020 data (coe31(k),k=0,6)/ &
1021 -1.250541e+000_r_kind, 1.911161e-002_r_kind, &
1022 -5.460238e-005_r_kind, -1.266388e-002_r_kind, &
1023 5.745064e-005_r_kind, 1.313985e-002_r_kind, &
1024 -4.574811e-005_r_kind/
1026 ! Fitting Coefficients at 50.3 GHz: Using Tb1 ~ Tb3
1027 data (coe50(k),k=0,6)/ &
1028 -1.246754e+000_r_kind, 2.368658e-002_r_kind, &
1029 -8.061774e-005_r_kind, -3.206323e-002_r_kind, &
1030 1.148107e-004_r_kind, 2.688353e-002_r_kind, &
1031 -7.358356e-005_r_kind/
1033 ! Fitting Coefficients at 89 GHz: Using Tb1 ~ Tb4
1034 data (coe89(k),k=0,8)/ &
1035 -1.278780e+000_r_kind, 1.625141e-002_r_kind, &
1036 -4.764536e-005_r_kind, -1.475181e-002_r_kind, &
1037 5.107766e-005_r_kind, 1.083021e-002_r_kind, &
1038 -4.154825e-005_r_kind, 7.703879e-003_r_kind, &
1039 -6.351148e-006_r_kind/
1041 ! Fitting Coefficients at 150 GHz: Using Tb1 ~ Tb5
1042 data coe150/&
1043 -1.691077e+000_r_kind, 3.352403e-002_r_kind, &
1044 -7.310338e-005_r_kind, -4.396138e-002_r_kind, &
1045 1.028994e-004_r_kind, 2.301014e-002_r_kind, &
1046 -7.070810e-005_r_kind, 1.270231e-002_r_kind, &
1047 -2.139023e-005_r_kind, -2.257991e-003_r_kind, &
1048 1.269419e-005_r_kind/
1050 save coe23,coe31,coe50,coe89,coe150
1052 ! Calculate emissivity discriminators at five AMSU window channels
1053 do ich = 1, nwch
1054 discriminator(ich) = coe(1+(ich-1)*11)
1055 if (ich .le. 3) nvalid_ch = 3
1056 if (ich .eq. 4) nvalid_ch = 4
1057 if (ich .eq. 5) nvalid_ch = 5
1058 do i=1,nvalid_ch
1059 discriminator(ich) = discriminator(ich) + coe((ich-1)*11 + 2*i)*tb(i) + &
1060 coe((ich-1)*11 + 2*i+1)*tb(i)*tb(i)
1061 end do
1062 end do
1063 ! Identify one snow emissivity spectrum and interpolate/extrapolate emissivity
1064 ! at a required frequency
1065 call em_interpolate(frequency,discriminator,emissivity,snow_type)
1067 em_vector(1) = emissivity
1068 em_vector(2) = emissivity
1070 return
1071 end subroutine sem_AB
1074 !---------------------------------------------------------------------!
1075 subroutine sem_ATs(theta,frequency,tba,ts,snow_type,em_vector)
1077 !$$$ subprogram documentation block
1078 !
1079 ! subprogram:
1080 !
1081 ! prgmmr:Banghua Yan org: nesdis date: 2003-08-18
1082 !
1083 ! abstract:
1084 ! Calculate the emissivity discriminators and interpolate/extrapolate
1085 ! emissivity at required frequency with respect to secenery AMSUAB
1086 !
1087 ! program history log:
1088 !
1089 ! input argument list:
1090 !
1091 ! frequency - frequency in GHz
1092 ! theta - local zenith angle (not used here)
1093 ! ts - surface temperature
1094 ! tba[1] ~ tba[4] - brightness temperature at five AMSU-A window channels:
1095 ! tba[1] : 23.8 GHz
1096 ! tba[2] : 31.4 GHz
1097 ! tba[3] : 50.3 GHz
1098 ! tba[4] : 89 GHz
1099 ! output argument list:
1100 !
1101 ! em_vector[1] and [2] - emissivity at two polarizations.
1102 ! set esv = esh here and will be updated
1103 ! snow_type - snow type (reference [2])
1104 !
1105 ! important internal variables:
1106 !
1107 ! coe23 - fitting coefficients to estimate discriminator at 23.8 GHz
1108 ! coe31 - fitting coefficients to estimate discriminator at 31.4 GHz
1109 ! coe50 - fitting coefficients to estimate discriminator at 50.3 GHz
1110 ! coe89 - fitting coefficients to estimate discriminator at 89 GHz
1111 ! coe150 - fitting coefficients to estimate discriminator at 150 GHz
1112 !
1113 ! remarks:
1114 !
1115 ! attributes:
1116 ! language: f90
1117 ! machine: ibm rs/6000 sp
1118 !
1119 !$$$
1121 ! use kinds, only: r_kind,i_kind
1122 implicit none
1124 integer(i_kind),parameter:: nch =10,nwch = 5,ncoe = 9
1125 real(r_kind) :: tba(*),theta
1126 real(r_kind) :: em_vector(*),emissivity,ts,frequency,discriminator(nwch)
1127 integer(i_kind) :: snow_type,i,k,ich,nvalid_ch
1128 real(r_kind) :: coe23(0:ncoe),coe31(0:ncoe),coe50(0:ncoe),coe89(0:ncoe),coe150(0:ncoe)
1129 real(r_kind) :: coe(nch*(ncoe+1))
1131 Equivalence (coe(1),coe23)
1132 Equivalence (coe(11),coe31)
1133 Equivalence (coe(21),coe50)
1134 Equivalence (coe(31),coe89)
1135 Equivalence (coe(41),coe150)
1137 ! Fitting Coefficients at 23.8 GHz: Using Tb1, Tb2 and Ts
1138 data (coe23(k),k=0,5)/ &
1139 8.210105e-001_r_kind, 1.216432e-002_r_kind, &
1140 -2.113875e-005_r_kind, -6.416648e-003_r_kind, &
1141 1.809047e-005_r_kind, -4.206605e-003_r_kind/
1143 ! Fitting Coefficients at 31.4 GHz: Using Tb1, Tb2 and Ts
1144 data (coe31(k),k=0,5)/ &
1145 7.963632e-001_r_kind, 7.215580e-003_r_kind, &
1146 -2.015921e-005_r_kind, -1.508286e-003_r_kind, &
1147 1.731405e-005_r_kind, -4.105358e-003_r_kind/
1149 ! Fitting Coefficients at 50.3 GHz: Using Tb1, Tb2, Tb3 and Ts
1150 data (coe50(k),k=0,7)/ &
1151 1.724160e+000_r_kind, 5.556665e-003_r_kind, &
1152 -2.915872e-005_r_kind, -1.146713e-002_r_kind, &
1153 4.724243e-005_r_kind, 3.851791e-003_r_kind, &
1154 -5.581535e-008_r_kind, -5.413451e-003_r_kind/
1156 ! Fitting Coefficients at 89 GHz: Using Tb1 ~ Tb4 and Ts
1157 data coe89/ &
1158 9.962065e-001_r_kind, 1.584161e-004_r_kind, &
1159 -3.988934e-006_r_kind, 3.427638e-003_r_kind, &
1160 -5.084836e-006_r_kind, -6.178904e-004_r_kind, &
1161 1.115315e-006_r_kind, 9.440962e-004_r_kind, &
1162 9.711384e-006_r_kind, -4.259102e-003_r_kind/
1164 ! Fitting Coefficients at 150 GHz: Using Tb1 ~ Tb4 and Ts
1165 data coe150/ &
1166 -5.244422e-002_r_kind, 2.025879e-002_r_kind, &
1167 -3.739231e-005_r_kind, -2.922355e-002_r_kind, &
1168 5.810726e-005_r_kind, 1.376275e-002_r_kind, &
1169 -3.757061e-005_r_kind, 6.434187e-003_r_kind, &
1170 6.190403e-007_r_kind, -2.944785e-003_r_kind/
1172 save coe23,coe31,coe50,coe89,coe150
1174 ! Calculate emissivity discriminators at five AMSU window channels
1175 DO ich = 1, nwch
1176 discriminator(ich) = coe(1+(ich-1)*10)
1177 if (ich .le. 2) nvalid_ch = 2
1178 if (ich .eq. 3) nvalid_ch = 3
1179 if (ich .ge. 4) nvalid_ch = 4
1180 do i=1,nvalid_ch
1181 discriminator(ich) = discriminator(ich) + coe((ich-1)*10 + 2*i)*tba(i) + &
1182 coe((ich-1)*10 + 2*i+1)*tba(i)*tba(i)
1183 end do
1184 discriminator(ich) = discriminator(ich) + coe( (ich-1)*10 + (nvalid_ch+1)*2 )*ts
1185 end do
1187 call em_interpolate(frequency,discriminator,emissivity,snow_type)
1189 em_vector(1) = emissivity
1190 em_vector(2) = emissivity
1192 return
1193 end subroutine sem_ATs
1195 !---------------------------------------------------------------------!
1196 subroutine sem_amsua(theta,frequency,tba,snow_type,em_vector)
1198 !$$$ subprogram documentation block
1199 !
1200 ! subprogram:
1201 !
1202 ! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
1203 !
1204 ! abstract:
1205 ! Calculate the emissivity discriminators and interpolate/extrapolate
1206 ! emissivity at required frequency with respect to secenery AMSUA
1207 !
1208 ! program history log:
1209 ! 2004-10-28 treadon - correct problem in declared dimensions of array coe
1210 !
1211 ! input argument list:
1212 !
1213 ! frequency - frequency in GHz
1214 ! theta - local zenith angle (not used here)
1215 ! tba[1]~tba[4] - brightness temperature at five AMSU-A window channels:
1216 ! tba[1] : 23.8 GHz
1217 ! tba[2] : 31.4 GHz
1218 ! tba[3] : 50.3 GHz
1219 ! tba[4] : 89 GHz
1220 !
1221 ! output argument list:
1222 !
1223 ! em_vector(1) and (2) - emissivity at two polarizations.
1224 ! set esv = esh here and will be updated
1225 ! snow_type - snow type
1226 !
1227 ! important internal variables:
1228 !
1229 ! coe23 - fitting coefficients to estimate discriminator at 23.8 GHz
1230 ! coe31 - fitting coefficients to estimate discriminator at 31.4 GHz
1231 ! coe50 - fitting coefficients to estimate discriminator at 50.3 GHz
1232 ! coe89 - fitting coefficients to estimate discriminator at 89 GHz
1233 ! coe150 - fitting coefficients to estimate discriminator at 150 GHz
1234 !
1235 ! remarks:
1236 !
1237 ! attributes:
1238 ! language: f90
1239 ! machine: ibm rs/6000 sp
1240 !
1241 !$$$
1243 ! use kinds, only: r_kind,i_kind
1244 implicit none
1246 integer(i_kind),parameter:: nch =10,nwch = 5,ncoe = 8
1247 real(r_kind) :: tba(*),theta
1248 real(r_kind) :: em_vector(*),emissivity,frequency,discriminator(nwch)
1249 integer(i_kind) :: snow_type,i,k,ich,nvalid_ch
1250 real(r_kind) :: coe23(0:ncoe),coe31(0:ncoe),coe50(0:ncoe),coe89(0:ncoe),coe150(0:ncoe)
1251 real(r_kind) :: coe(nch*(ncoe+1))
1253 Equivalence (coe(1),coe23)
1254 Equivalence (coe(11),coe31)
1255 Equivalence (coe(21),coe50)
1256 Equivalence (coe(31),coe89)
1257 Equivalence (coe(41),coe150)
1259 ! Fitting Coefficients at 23.8 GHz: Using Tb1 ~ Tb3
1260 data (coe23(k),k=0,6)/ &
1261 -1.326040e+000_r_kind, 2.475904e-002_r_kind, -5.741361e-005_r_kind, &
1262 -1.889650e-002_r_kind, 6.177911e-005_r_kind, 1.451121e-002_r_kind, &
1263 -4.925512e-005_r_kind/
1265 ! Fitting Coefficients at 31.4 GHz: Using Tb1 ~ Tb3
1266 data (coe31(k),k=0,6)/ &
1267 -1.250541e+000_r_kind, 1.911161e-002_r_kind, -5.460238e-005_r_kind, &
1268 -1.266388e-002_r_kind, 5.745064e-005_r_kind, 1.313985e-002_r_kind, &
1269 -4.574811e-005_r_kind/
1271 ! Fitting Coefficients at 50.3 GHz: Using Tb1 ~ Tb3
1272 data (coe50(k),k=0,6)/ &
1273 -1.246754e+000_r_kind, 2.368658e-002_r_kind, -8.061774e-005_r_kind, &
1274 -3.206323e-002_r_kind, 1.148107e-004_r_kind, 2.688353e-002_r_kind, &
1275 -7.358356e-005_r_kind/
1277 ! Fitting Coefficients at 89 GHz: Using Tb1 ~ Tb4
1278 data coe89/ &
1279 -1.278780e+000_r_kind, 1.625141e-002_r_kind, -4.764536e-005_r_kind, &
1280 -1.475181e-002_r_kind, 5.107766e-005_r_kind, 1.083021e-002_r_kind, &
1281 -4.154825e-005_r_kind, 7.703879e-003_r_kind, -6.351148e-006_r_kind/
1283 ! Fitting Coefficients at 150 GHz: Using Tb1 ~ Tb4
1284 data coe150/ &
1285 -1.624857e+000_r_kind, 3.138243e-002_r_kind, -6.757028e-005_r_kind, &
1286 -4.178496e-002_r_kind, 9.691893e-005_r_kind, 2.165964e-002_r_kind, &
1287 -6.702349e-005_r_kind, 1.111658e-002_r_kind, -1.050708e-005_r_kind/
1289 save coe23,coe31,coe50,coe150
1292 ! Calculate emissivity discriminators at five AMSU window channels
1293 do ich = 1, nwch
1294 discriminator(ich) = coe(1+(ich-1)*10)
1295 if (ich .le. 2) nvalid_ch = 3
1296 if (ich .ge. 3) nvalid_ch = 4
1297 do i=1,nvalid_ch
1298 discriminator(ich) = discriminator(ich) + coe((ich-1)*10 + 2*i)*tba(i) + &
1299 coe((ich-1)*10 + 2*i+1)*tba(i)*tba(i)
1300 end do
1301 end do
1303 ! Quality Control
1304 if(discriminator(4) .gt. discriminator(2)) &
1305 discriminator(4) = discriminator(2) + (150.0_r_kind - 89.0_r_kind)* &
1306 (discriminator(5) - discriminator(2))/ &
1307 (150.0_r_kind - 31.4_r_kind)
1309 ! Quality control at 50.3 GHz
1310 if((discriminator(3) .gt. discriminator(2)) .or. &
1311 (discriminator(3) .lt. discriminator(4))) &
1312 discriminator(3) = discriminator(2) + (89.0_r_kind - 50.3_r_kind)* &
1313 (discriminator(4) - discriminator(2))/(89.0_r_kind - 31.4_r_kind)
1315 call em_interpolate(frequency,discriminator,emissivity,snow_type)
1317 em_vector(1) = emissivity
1318 em_vector(2) = emissivity
1320 return
1321 end subroutine sem_amsua
1324 !---------------------------------------------------------------------!
1325 subroutine sem_BTs(theta,frequency,tbb,ts,snow_type,em_vector)
1327 !$$$ subprogram documentation block
1328 !
1329 ! subprogram:
1330 !
1331 ! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
1332 !
1333 ! abstract:
1334 ! Calculate the emissivity discriminators and interpolate/extrapolate
1335 ! emissivity at required frequency with respect to secenery BTs
1336 !
1337 ! program history log:
1338 !
1339 ! input argument list:
1340 !
1341 ! frequency - frequency in GHz
1342 ! theta - local zenith angle (not used here)
1343 ! ts - surface temperature in degree
1344 ! tbb[1] ~ tbb[2] - brightness temperature at five AMSU-B window channels:
1345 ! tbb[1] : 89 GHz
1346 ! tbb[2] : 150 GHz
1347 !
1348 ! output argument list:
1349 !
1350 ! em_vector(1) and (2) - emissivity at two polarizations.
1351 ! set esv = esh here and will be updated
1352 ! snow_type - snow type
1353 !
1354 ! important internal variables:
1355 !
1356 ! coe31 - fitting coefficients to estimate discriminator at 31.4 GHz
1357 ! coe89 - fitting coefficients to estimate discriminator at 89 GHz
1358 ! coe150 - fitting coefficients to estimate discriminator at 150 GHz
1359 !
1360 ! remarks:
1361 !
1362 ! attributes:
1363 ! language: f90
1364 ! machine: ibm rs/6000 sp
1365 !
1366 !$$$
1367 ! use kinds, only: r_kind,i_kind
1368 implicit none
1370 integer(i_kind),parameter:: nch =10,nwch = 3,ncoe = 5
1371 real(r_kind) :: tbb(*),theta
1372 real(r_kind) :: em_vector(*),emissivity,ts,frequency,ed0(nwch),discriminator(5)
1373 integer(i_kind) :: snow_type,i,k,ich,nvalid_ch
1374 real(r_kind) :: coe31(0:ncoe),coe89(0:ncoe),coe150(0:ncoe)
1375 real(r_kind) :: coe(nch*(ncoe+1))
1377 Equivalence (coe(1),coe31)
1378 Equivalence (coe(11),coe89)
1379 Equivalence (coe(21),coe150)
1381 ! Fitting Coefficients at 31.4 GHz: Using Tb4, Tb5 and Ts
1382 data coe31/ 3.110967e-001_r_kind, 1.100175e-002_r_kind, -1.677626e-005_r_kind, &
1383 -4.020427e-003_r_kind, 9.242240e-006_r_kind, -2.363207e-003_r_kind/
1384 ! Fitting Coefficients at 89 GHz: Using Tb4, Tb5 and Ts
1385 data coe89/ 1.148098e+000_r_kind, 1.452926e-003_r_kind, 1.037081e-005_r_kind, &
1386 1.340696e-003_r_kind, -5.185640e-006_r_kind, -4.546382e-003_r_kind /
1387 ! Fitting Coefficients at 150 GHz: Using Tb4, Tb5 and Ts
1388 data coe150/ 1.165323e+000_r_kind, -1.030435e-003_r_kind, 4.828009e-006_r_kind, &
1389 4.851731e-003_r_kind, -2.588049e-006_r_kind, -4.990193e-003_r_kind/
1390 save coe31,coe89,coe150
1392 ! Calculate emissivity discriminators at five AMSU window channels
1393 do ich = 1, nwch
1394 ed0(ich) = coe(1+(ich-1)*10)
1395 nvalid_ch = 2
1396 do i=1,nvalid_ch
1397 ed0(ich) = ed0(ich) + coe((ich-1)*10 + 2*i)*tbb(i) + &
1398 coe((ich-1)*10 + 2*i+1)*tbb(i)*tbb(i)
1399 end do
1400 ed0(ich) = ed0(ich) + coe( (ich-1)*10 + (nvalid_ch+1)*2 )*ts
1401 end do
1403 ! Quality control
1404 if(ed0(2) .gt. ed0(1)) &
1405 ed0(2) = ed0(1) + (150.0_r_kind - 89.0_r_kind)*(ed0(3) - ed0(1)) / &
1406 (150.0_r_kind - 31.4_r_kind)
1408 ! Match the format of the input variable
1409 ! Missing value at 23.8 GHz
1410 discriminator(1) = -999.9_r_kind; discriminator(2) = ed0(1)
1411 ! Missing value at 50.3 GHz
1412 discriminator(3) = -999.9_r_kind; discriminator(4) = ed0(2); discriminator(5) = ed0(3)
1414 call em_interpolate(frequency,discriminator,emissivity,snow_type)
1416 em_vector(1) = emissivity
1417 em_vector(2) = emissivity
1419 return
1420 end subroutine sem_BTs
1423 !---------------------------------------------------------------------!
1424 subroutine sem_amsub(theta,frequency,tbb,snow_type,em_vector)
1427 !$$$ subprogram documentation block
1428 !
1429 ! subprogram:
1430 !
1431 ! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
1432 !
1433 ! abstract:
1434 ! Calculate the emissivity discriminators and interpolate/extrapolate
1435 ! emissivity at required frequency with respect to secenery AMSUB
1436 !
1437 ! program history log:
1438 ! 2004-10-28 treadon - correct problem in declared dimensions of array coe
1439 !
1440 ! input argument list:
1441 !
1442 ! frequency - frequency in GHz
1443 ! theta - local zenith angle (not used here)
1444 ! tbb[1] ~ tbb[2] - brightness temperature at five AMSU-B window channels:
1445 ! tbb[1] : 89 GHz
1446 ! tbb[2] : 150 GHz
1447 !
1448 ! output argument list:
1449 ! em_vector(1) and (2) - emissivity at two polarizations.
1450 ! set esv = esh here and will be updated
1451 ! snow_type - snow type (reference [2])
1452 !
1453 ! important internal variables:
1454 !
1455 ! coe31 - fitting coefficients to estimate discriminator at 31.4 GHz
1456 ! coe89 - fitting coefficients to estimate discriminator at 89 GHz
1457 ! coe150 - fitting coefficients to estimate discriminator at 150 GHz
1458 !
1459 ! remarks:
1460 !
1461 ! attributes:
1462 ! language: f90
1463 ! machine: ibm rs/6000 sp
1464 !
1465 !$$$
1466 ! use kinds, only: r_kind,i_kind
1467 implicit none
1469 integer(i_kind),parameter:: nch =10,nwch = 3,ncoe = 4
1470 real(r_kind) :: tbb(*)
1471 real(r_kind) :: em_vector(*),emissivity,frequency,ed0(nwch),discriminator(5)
1472 integer(i_kind) :: snow_type,i,k,ich,nvalid_ch
1473 real(r_kind) :: coe31(0:ncoe),coe89(0:ncoe),coe150(0:ncoe)
1474 real(r_kind) :: coe(nch*(ncoe+1))
1475 real(r_kind) :: theta,dem,demmin0
1477 Equivalence (coe(1),coe31)
1478 Equivalence (coe(11),coe89)
1479 Equivalence (coe(21),coe150)
1481 ! Fitting Coefficients at 31.4 GHz: Using Tb4, Tb5
1482 data coe31/-4.015636e-001_r_kind,9.297894e-003_r_kind, -1.305068e-005_r_kind, &
1483 3.717131e-004_r_kind, -4.364877e-006_r_kind/
1484 ! Fitting Coefficients at 89 GHz: Using Tb4, Tb5
1485 data coe89/-2.229547e-001_r_kind, -1.828402e-003_r_kind,1.754807e-005_r_kind, &
1486 9.793681e-003_r_kind, -3.137189e-005_r_kind/
1487 ! Fitting Coefficients at 150 GHz: Using Tb4, Tb5
1488 data coe150/-3.395416e-001_r_kind,-4.632656e-003_r_kind,1.270735e-005_r_kind, &
1489 1.413038e-002_r_kind,-3.133239e-005_r_kind/
1490 save coe31,coe89,coe150
1492 ! Calculate emissivity discriminators at five AMSU window channels
1493 do ich = 1, nwch
1494 ed0(ich) = coe(1+(ich-1)*10)
1495 nvalid_ch = 2
1496 do i=1,nvalid_ch
1497 ed0(ich) = ed0(ich) + coe((ich-1)*10 + 2*i)*tbb(i) + &
1498 coe((ich-1)*10 + 2*i+1)*tbb(i)*tbb(i)
1499 end do
1500 end do
1502 ! Quality Control
1503 if(ed0(2) .gt. ed0(1)) &
1504 ed0(2) = ed0(1) + (150.0_r_kind - 89.0_r_kind) * &
1505 (ed0(3) - ed0(1))/(150.0_r_kind - 31.4_r_kind)
1507 ! Match the format of the input variable
1508 ! Missing value at 23.8 GHz
1509 discriminator(1) = -999.9_r_kind; discriminator(2) = ed0(1)
1510 ! Missing value at 50.3 GHz
1511 discriminator(3) = -999.9_r_kind; discriminator(4) = ed0(2); discriminator(5) = ed0(3)
1513 call em_interpolate(frequency,discriminator,emissivity,snow_type)
1515 em_vector(1) = emissivity
1516 em_vector(2) = emissivity
1518 return
1519 end subroutine sem_amsub
1522 !---------------------------------------------------------------------!
1523 subroutine ALandEM_Snow(theta,frequency,snow_depth,t_skin,snow_type,em_vector)
1526 !$$$ subprogram documentation block
1527 !
1528 ! subprogram:
1529 !
1530 ! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
1531 !
1532 ! abstract:
1533 ! Calculate the emissivity at required frequency with respect to option MODL
1534 ! using the LandEM and a bias correction algorithm, where the original LandEM with a
1535 ! bias correction algorithm is referred to as value-added LandEM or AlandEM.
1536 !
1537 ! program history log:
1538 !
1539 ! input argument list:
1540 !
1541 ! frequency - frequency in GHz
1542 ! theta - local zenith angle (not used here)
1543 ! snow_depth - snow depth in mm
1544 ! t_skin - surface temperature
1545 !
1546 ! output argument list:
1547 !
1548 ! em_vector(1) and (2) - emissivity at two polarizations.
1549 ! set esv = esh here and will be updated
1550 ! snow_type - snow type
1551 !
1552 ! important internal variables:
1553 !
1554 ! esv_3w and esh_3w - initial emissivity discriminator at two polarizations
1555 ! at three AMSU window channels computed using LandEM
1556 ! esv_3w[1] and esh_3w[1] : 31.4 GHz
1557 ! esv_3w[2] and esh_3w[2] : 89 GHz
1558 ! esv_3w[3] and esh_3w[3] : 150 GHz
1559 !
1560 ! remarks:
1561 !
1562 ! attributes:
1563 ! language: f90
1564 ! machine: ibm rs/6000 sp
1565 !
1566 !$$$
1568 ! use kinds, only: r_kind,i_kind
1569 ! use constants, only: zero, one
1570 implicit none
1572 integer(i_kind) :: nw_ind
1573 parameter(nw_ind=3)
1574 real(r_kind) theta, frequency, freq,snow_depth, mv, t_soil, t_skin, em_vector(2)
1575 real(r_kind) esv,esh,esh0,esv0,theta0,b
1576 integer(i_kind) snow_type,ich
1577 real(r_kind) freq_3w(nw_ind),esh_3w(nw_ind),esv_3w(nw_ind)
1578 complex(r_kind) eair
1579 data freq_3w/31.4_r_kind,89.0_r_kind,150.0_r_kind/
1581 eair = dcmplx(one,-zero)
1582 b = t_skin
1583 snow_type = -999
1585 call snowem_default(theta,frequency,snow_depth,t_skin,b,esv0,esh0)
1587 theta0 = theta
1588 do ich = 1, nw_ind
1589 freq =freq_3w(ich)
1590 theta = theta0
1591 call snowem_default(theta,freq,snow_depth,t_skin,b,esv,esh)
1592 esv_3w(ich) = esv
1593 esh_3w(ich) = esh
1594 end do
1596 call ems_adjust(theta,frequency,snow_depth,t_skin,esv_3w,esh_3w,em_vector,snow_type)
1598 return
1599 end subroutine ALandEM_Snow
1602 !---------------------------------------------------------------------!
1603 subroutine snowem_default(theta,freq,snow_depth,t_soil,b,esv,esh)
1605 !$$$ subprogram documentation block
1606 !
1607 ! subprogram:
1608 !
1609 ! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
1610 !
1611 ! abstract:
1612 ! Initialize discriminator using LandEM
1613 !
1614 ! program history log:
1615 ! 2004-09-22 todling - replaced zsqrt by general interface sqrt
1616 !
1617 ! input argument list:
1618 !
1619 ! frequency - frequency in GHz
1620 ! theta - local zenith angle in radian
1621 ! snow_depth - snow depth in mm
1622 ! t_skin - surface temperature
1623 !
1624 ! output argument list:
1625 !
1626 ! esv - initial discriminator at vertical polarization
1627 ! esh - at horizontal polarization
1628 !
1629 ! remarks:
1630 !
1631 ! attributes:
1632 ! language: f90
1633 ! machine: ibm rs/6000 sp
1634 !
1635 !$$$
1636 ! use kinds, only: r_kind
1637 ! use constants, only: zero, one
1638 implicit none
1640 real(r_kind) rhob,rhos,sand,clay
1641 Parameter(rhob = 1.18_r_kind, rhos = 2.65_r_kind, &
1642 sand = 0.8_r_kind, clay = 0.2_r_kind)
1643 real(r_kind) theta, freq, mv, t_soil, snow_depth,b
1644 real(r_kind) theta_i,theta_t, mu, r12_h, r12_v, r21_h, r21_v, r23_h, r23_v, &
1645 t21_v, t21_h, t12_v, t12_h, gv, gh, ssalb_h,ssalb_v,tau_h, &
1646 tau_v, esh, esv,rad, sigma, va ,ep_real,ep_imag
1647 complex(r_kind) esoil, esnow, eair
1649 eair = dcmplx(one,-zero)
1650 ! ep = dcmplx(3.2_r_kind,-0.0005_r_kind)
1651 sigma = one
1652 theta_i = theta
1653 mv = 0.1_r_kind
1654 ep_real = 3.2_r_kind
1655 ep_imag = -0.0005_r_kind
1656 va = 0.4_r_kind + 0.0004_r_kind*snow_depth
1657 rad = one + 0.005_r_kind*snow_depth
1659 call snow_diel(freq, ep_real, ep_imag, rad, va, esnow)
1660 ! call snow_diel(freq, ep, rad, va, esnow)
1661 call soil_diel(freq, t_soil, mv, rhob, rhos, sand, clay, esoil)
1662 theta_t = asin(real(sin(theta_i)*sqrt(eair)/sqrt(esnow)))
1663 call reflectance(eair, esnow, theta_i, theta_t, r12_v, r12_h)
1664 call transmitance(eair, esnow, theta_i, theta_t, t12_v, t12_h)
1666 theta_t = theta
1667 theta_i = asin(real(sin(theta_t)*sqrt(eair)/sqrt(esnow)))
1668 call reflectance(esnow, eair, theta_i, theta_t, r21_v, r21_h)
1669 call transmitance(esnow, eair, theta_i, theta_t, t21_v, t21_h)
1671 mu = cos(theta_i)
1672 theta_t = asin(real(sin(theta_i)*sqrt(esnow)/sqrt(esoil)))
1673 call reflectance(esnow, esoil, theta_i, theta_t, r23_v, r23_h)
1674 call rough_reflectance(freq, theta_i, sigma, r23_v, r23_h)
1676 ! call snow_optic(freq, rad, snow_depth, va, ep, gv, gh, ssalb_v, ssalb_h, tau_v, tau_h)
1677 call snow_optic(freq,rad,snow_depth,va,ep_real, ep_imag,gv,gh,&
1678 ssalb_v,ssalb_h,tau_v,tau_h)
1680 call two_stream_solution(b,mu,gv,gh,ssalb_h, ssalb_v, tau_h, tau_v, r12_h, &
1681 r12_v, r21_h, r21_v, r23_h, r23_v, t21_v, t21_h, t12_v, t12_h, esv, esh)
1682 return
1683 end subroutine snowem_default
1686 !---------------------------------------------------------------------!
1687 subroutine ems_adjust(theta,frequency,depth,ts,esv_3w,esh_3w,em_vector,snow_type)
1690 !$$$ subprogram documentation block
1691 !
1692 ! subprogram:
1693 !
1694 ! prgmmr: Banghua Yan org: nesdis date: 2003-08-18
1695 !
1696 ! abstract:
1697 ! Calculate the emissivity discriminators and interpolate/extrapolate
1698 ! emissivity at required frequency with respect to secenery MODL
1699 !
1700 ! program history log:
1701 ! 2004-10-28 treadon - remove nch from parameter declaration below (not used)
1702 !
1703 ! input argument list:
1704 !
1705 ! frequency - frequency in GHz
1706 ! theta - local zenith angle (not used here)
1707 ! depth - snow depth in mm
1708 ! ts - surface temperature
1709 !
1710 ! output argument list:
1711 !
1712 ! em_vector(1) and (2) - emissivity at two polarizations.
1713 ! set esv = esh here and will be updated
1714 ! snow_type - snow type
1715 !
1716 ! important internal variables:
1717 !
1718 ! dem_coe - fitting coefficients to compute discriminator correction value
1719 ! dem_coe[1,*] : 31.4 GHz
1720 ! dem_coe[2,*] : 89 GHz
1721 ! dem_coe[3,*] : 150 GHz
1722 !
1723 ! remarks:
1724 !
1725 ! attributes:
1726 ! language: f90
1727 ! machine: ibm rs/6000 sp
1728 !
1729 !$$$
1731 ! use kinds, only: r_kind,r_double,i_kind
1732 ! use constants, only: one,deg2rad
1733 implicit none
1735 integer(i_kind),parameter:: nw_3=3
1736 integer(i_kind),parameter:: ncoe=6
1737 real(r_kind),parameter :: earthrad = 6374._r_kind, satheight = 833.4_r_kind
1738 integer(i_kind) :: snow_type,ich,j,k
1739 real(r_kind) :: theta,frequency,depth,ts,esv_3w(*),esh_3w(*)
1740 real(r_kind) :: discriminator(5),emmod(nw_3),dem(nw_3)
1741 real(r_kind) :: emissivity,em_vector(2)
1742 real(r_double) :: dem_coe(nw_3,0:ncoe-1),sinthetas,costhetas
1744 save dem_coe
1746 data (dem_coe(1,k),k=0,ncoe-1)/ 2.306844e+000_r_double, -7.287718e-003_r_double, &
1747 -6.433248e-004_r_double, 1.664216e-005_r_double, &
1748 4.766508e-007_r_double, -1.754184e+000_r_double/
1749 data (dem_coe(2,k),k=0,ncoe-1)/ 3.152527e+000_r_double, -1.823670e-002_r_double, &
1750 -9.535361e-004_r_double, 3.675516e-005_r_double, &
1751 9.609477e-007_r_double, -1.113725e+000_r_double/
1752 data (dem_coe(3,k),k=0,ncoe-1)/ 3.492495e+000_r_double, -2.184545e-002_r_double, &
1753 6.536696e-005_r_double, 4.464352e-005_r_double, &
1754 -6.305717e-008_r_double, -1.221087e+000_r_double/
1756 sinthetas = sin(theta*deg2rad)* earthrad/(earthrad + satheight)
1757 sinthetas = sinthetas*sinthetas
1758 costhetas = one - sinthetas
1759 do ich = 1, nw_3
1760 emmod(ich) = costhetas*esv_3w(ich) + sinthetas*esh_3w(ich)
1761 end do
1762 do ich=1,nw_3
1763 dem(ich) = dem_coe(ich,0) + dem_coe(ich,1)*ts + dem_coe(ich,2)*depth + &
1764 dem_coe(ich,3)*ts*ts + dem_coe(ich,4)*depth*depth + &
1765 dem_coe(ich,5)*emmod(ich)
1766 end do
1767 emmod(1) = emmod(1) + dem(1)
1768 emmod(2) = emmod(2) + dem(2)
1769 emmod(3) = emmod(3) + dem(3)
1771 ! Match the format of the input variable
1773 ! Missing value at 23.8 GHz
1774 discriminator(1) = -999.9_r_kind; discriminator(2) = emmod(1)
1776 ! Missing value at 50.3 GHz
1777 discriminator(3) = -999.9_r_kind; discriminator(4) = emmod(2); discriminator(5) = emmod(3)
1779 call em_interpolate(frequency,discriminator,emissivity,snow_type)
1781 em_vector(1) = emissivity
1782 em_vector(2) = emissivity
1784 return
1785 end subroutine ems_adjust