| blob | 6daff9b655a8814a30e80eef975565773765185c |
1 subroutine landem(theta,freq,mv,veg_frac,veg_tp,soil_tp, &
2 t_soil,t_skin,snow_depth,esh,esv)
4 !$$$ subprogram documentation block
5 ! . . . .
6 ! subprogram: landem noaa/nesdis emissivity model over land/ice
7 !
8 ! prgmmr: Fuzhong Weng org: nesdis date: 2000-11-28
9 ! Banghua Yan
10 !
11 ! abstract: noaa/nesdis emissivity model to compute microwave emissivity over
12 ! various land surface conditions (land/snow)
13 !
14 ! reference: weng, f, b. yan, and n. grody, 2001:
15 ! "A microwave land emissivity model", J. Geophys. Res., 106,
16 ! 20, 115-20, 123
17 !
18 ! version: beta
19 !
20 ! program history log:
21 ! 1999-07-01 weng
22 ! 2000-01-01 yan - add canopy scattering
23 ! 2000-11-28 weng,yan - parameterize model input and incorporate into NCEP ssi
24 ! 2004-08-04 treadon - add only on use declarations; add intent in/out
25 ! 2004-09-22 todling - replaced zsqrt by general interface sqrt
26 !
27 ! input argument list:
28 !
29 ! theta - local zenith angle (0 - 60.0)
30 ! freq - frequency in ghz ( 0 - 90.0)
31 ! mv - volumetric moisture content in soil (0.0 - 1.0) (gdas)
32 ! veg_frac - vegetation fraction (0 - 1.0) (gdas)
33 ! veg_tp - vegetation type (gdas, not used)
34 ! 1: broadleave evergreen trees
35 ! 2: broadleave deciduous trees
36 ! 3: broad & needle mixed forest
37 ! 4: needleleave evergreen trees
38 ! 5: needleleave deciduous trees
39 ! 6: broadleave tree with groundcover (savana)
40 ! 7: groundcover only (perenial groundcover)
41 ! 8: broadleave shrubs with perenial groundcover
42 ! 9: broadleave shrubs with bare soil
43 ! 10: dwarf trees & shrubs with bare soil
44 ! 11: bare soil'
45 ! 12: cultivations (use paramater 7)
46 ! 13: glacial
47 !
48 ! soil_tp - soil type (gdas, not used)
49 ! 1: loamy sand (coarse)
50 ! 2: silty clayloam (medium)
51 ! 3: light clay (fine)
52 ! 4: sand loam (coarse-medium)
53 ! 5: sandy clay (coarse-fine)
54 ! 6: clay loam (medium-fine)
55 ! 7: sandy clay loam (coarse-med-fine)
56 ! 8: loam (organic)
57 ! 9: ice (use loamy sand property)
58 !
59 ! t_soil - soil temperature (k) (gdas)
60 ! t_skin - scattering layer temperature (k) (gdas)
61 ! snow_depth - scatter medium depth (mm) (gdas)
62 !
63 !
64 ! output argument list:
65 !
66 ! esh - emissivity for horizontal polarization
67 ! esv - emissivity for vertical polarization
68 !
69 ! important internal variables:
70 !
71 ! rhob - bulk volume density of the soil (1.18-1.12)
72 ! rhos - density of the solids (2.65 g.cm^3 for solid soil material)
73 ! sand - sand fraction (sand + clay = 1.0)
74 ! clay - clay fraction
75 ! lai - leaf area index (eg. lai = 4.0 for corn leaves)
76 ! sigma - surface roughness formed between medium 1 and 2,
77 ! expressed as he standard deviation of roughtness height (mm)
78 ! leaf_thick -- leaf thickness (mm)
79 ! rad - radius of dense medium scatterers (mm)
80 ! va - fraction volume of dense medium scatterers(0.0 - 1.0)
81 ! ep - dielectric constant of ice or sand particles, complex value
82 ! (e.g, 3.0+i0.0)
83 !
84 !
85 ! remarks:
86 !
87 ! attributes:
88 ! language: f90
89 ! machine: ibm rs/6000 sp
90 !
91 !$$$
93 ! use kinds, only: r_kind
94 ! use constants, only: zero,one_tenth,half,one,three
95 implicit none
97 ! Declare passed variables
98 real(r_kind),intent(in):: theta,freq,mv,veg_frac,veg_tp,soil_tp,&
99 t_soil,t_skin,snow_depth
100 real(r_kind),intent(out):: esh,esv
102 ! Declare local parameters
103 real(r_kind),parameter:: rhob=1.18_r_kind
104 real(r_kind),parameter:: rhos = 2.65_r_kind
105 real(r_kind),parameter:: sand = 0.8_r_kind
106 real(r_kind),parameter:: clay = 0.2_r_kind
108 ! Declare local variables
109 real(r_kind) b,theta_i,theta_t,mu,r12_h,r12_v,r21_h,r21_v,r23_h,r23_v, &
110 t21_v,t21_h,t12_v,t12_h,gv,gh,ssalb_h,ssalb_v,tau_h,tau_v,mge,&
111 lai,leaf_thick,rad,sigma,va,ep_real,ep_imag
112 complex(r_kind) esoil, eveg, esnow, eair
114 eair = dcmplx(one,-zero)
115 ! if (snow_depth .gt.zero) then
116 if (snow_depth .gt.one_tenth) then
117 ! ice dielectric constant
118 ep_real = 3.2_r_kind
119 ep_imag = -0.0005_r_kind
120 sigma = one
121 va = 0.4_r_kind + 0.0004_r_kind*snow_depth
122 rad = one + 0.005_r_kind*snow_depth
123 call snow_diel(freq, ep_real, ep_imag, rad, va, esnow)
124 call soil_diel(freq, t_soil, mv, rhob, rhos, sand, clay, esoil)
125 ! theta_t = dasin(dreal(sin(theta)*sqrt(eair)/sqrt(esnow)))
126 ! call reflectance(eair, esnow, theta, theta_t, r12_v, r12_h)
127 ! call transmitance(eair, esnow, theta, theta_t, t12_v, t12_h)
128 theta_i = asin(real(sin(theta)*sqrt(eair)/sqrt(esnow)))
129 call reflectance(esnow, eair, theta_i, theta, r21_v, r21_h)
130 call transmitance(esnow, eair, theta_i, theta, t21_v, t21_h)
131 mu = cos(theta_i)
132 theta_t = asin(real(sin(theta_i)*sqrt(esnow)/sqrt(esoil)))
133 call reflectance(esnow, esoil, theta_i, theta_t, r23_v, r23_h)
134 call rough_reflectance(freq, theta_i, sigma, r23_v, r23_h)
135 call snow_optic(freq,rad,snow_depth,va,ep_real, ep_imag,gv,gh,&
136 ssalb_v,ssalb_h,tau_v,tau_h)
137 call two_stream_solution(t_skin,mu,gv,gh,ssalb_h,ssalb_v,tau_h,tau_v,r12_h, &
138 r12_v,r21_h,r21_v,r23_h,r23_v,t21_v,t21_h,t12_v,t12_h,esv,esh)
139 else
140 sigma = half
141 lai = three*veg_frac + half
142 mge = half*veg_frac
143 leaf_thick = 0.07_r_kind
144 mu = cos(theta)
145 ! r12_h = zero
146 ! r12_v = zero
147 r21_h = zero
148 r21_v = zero
149 t21_h = one
150 t21_v = one
151 ! t12_v = one
152 ! t12_h = one
153 call soil_diel(freq, t_soil, mv, rhob, rhos, sand, clay, esoil)
154 theta_t = asin(real(sin(theta)*sqrt(eair)/sqrt(esoil)))
155 call reflectance(eair, esoil, theta, theta_t, r23_v, r23_h)
156 call rough_reflectance(freq, theta, sigma, r23_v, r23_h)
157 call canopy_diel(freq, mge, eveg)
158 call canopy_optic(lai,freq,theta,eveg,leaf_thick,gv,gh,ssalb_v,ssalb_h,tau_v,tau_h)
159 call two_stream_solution(t_skin,mu,gv,gh,ssalb_h,ssalb_v,tau_h,tau_v,&
160 r12_h,r12_v,r21_h,r21_v,r23_h,r23_v,t21_v,t21_h,t12_v,t12_h,esv,esh)
161 endif
162 return
163 end subroutine landem
165 subroutine canopy_optic(lai,frequency,theta,esv,d,gv,gh,&
166 ssalb_v,ssalb_h,tau_v, tau_h)
168 !$$$ subprogram documentation block
169 ! . . . .
170 ! subprogram: canopy_optic compute optic parameters for canopy
171 !
172 ! prgmmr: org: nesdis date: 2000-11-28
173 !
174 ! abstract: compute optic parameters for canopy
175 !
176 ! program history log:
177 ! 2004-09-22 todling - replaced zsqrt by general interface sqrt; same for zexp/exp
178 !
179 ! input argument list:
180 !
181 ! lai - leaf area index
182 ! frequency - frequency (ghz)
183 ! theta - incident angle
184 ! esv - leaf dielectric constant
185 ! d - leaf thickness (mm)
186 !
187 ! output argument list:
188 !
189 ! gv - asymmetry factor for v pol
190 ! gh - asymmetry factor for h pol
191 ! ssalb_v - single scattering albedo at v. polarization
192 ! ssalb_h - single scattering albedo at h. polarization
193 ! tau_v - optical depth at v. polarization
194 ! tau_h - optical depth at h. polarization
195 !
196 ! remarks:
197 !
198 ! attributes:
199 ! language: f90
200 ! machine: ibm rs/6000 sp
201 !
202 !$$$
204 ! use kinds, only: r_kind
205 ! use constants, only: zero, half, one, two, four, pi
206 implicit none
207 real(r_kind) threshold
208 real(r_kind) frequency,theta,d,lai,ssalb_v,ssalb_h,tau_v,tau_h,gv, gh, mu
209 complex(r_kind) ix,k0,kz0,kz1,rhc,rvc,esv,expval1,factt,factrvc,factrhc
210 real(r_kind) rh,rvert,th,tv,av,ah,astar,rstar
211 threshold=0.999_r_kind
212 mu = cos(theta)
213 ix = dcmplx(zero,one)
214 k0 = dcmplx(two*pi*frequency/300.0_r_kind, zero) ! 1/mm
215 kz0 = k0*mu
216 kz1 = k0*sqrt(esv - sin(theta)**2)
217 rhc = (kz0 - kz1)/(kz0 + kz1)
218 rvc = (esv*kz0 - kz1)/(esv*kz0 + kz1)
219 expval1=exp(-two*ix*kz1*d)
220 factrvc=one-rvc**2*expval1
221 factrhc=one-rhc**2*expval1
222 factt=four*kz0*kz1*exp(ix*(kz0-kz1)*d)
223 rvert = abs(rvc*(one - expval1)/factrvc)**2
224 rh = abs(rhc*(one - expval1)/factrhc)**2
225 th = abs(factt/((kz1+kz0)**2*factrhc))**2
226 tv = abs(esv*factt/((kz1+esv*kz0)**2*factrvc))**2
227 ! av = one - rvert - tv
228 ! ah = one - rh - th
229 ! astar = av + ah
230 ! rstar = rvert + rh
231 ! randomly oriented leaves
232 gv = half
233 gh = half
234 ! tau_v = half*lai*(astar + rstar)
235 tau_v = half*lai*(two-tv-th)
236 tau_h = tau_v
237 ! tau_h = half*lai*(astar + rstar)
238 ! ssalb_v = rstar/ (astar + rstar)
239 ssalb_v = min((rvert+rh)/ (two-tv-th),threshold)
240 ssalb_h = ssalb_v
241 ! ssalb_h = rstar/ (astar + rstar)
242 return
243 end subroutine canopy_optic
245 subroutine snow_optic(frequency,a,h,f,ep_real,ep_imag,gv,gh,&
246 ssalb_v,ssalb_h,tau_v,tau_h)
248 !$$$ subprogram documentation block
249 ! . . . .
250 ! subprogram: landem comput optic parameters for snow
251 !
252 ! prgmmr: org: nesdis date: 2000-11-28
253 !
254 ! abstract: compute optic parameters for snow
255 !
256 ! program history log:
257 !
258 ! input argument list:
259 !
260 ! theta - local zenith angle (degree)
261 ! frequency - frequency (ghz)
262 ! ep_real - real part of dielectric constant of particles
263 ! ep_imag - imaginary part of dielectric constant of particles
264 ! a - particle radiu (mm)
265 ! h - snow depth(mm)
266 ! f - fraction volume of snow (0.0 - 1.0)
267 !
268 ! output argument list:
269 !
270 ! ssalb - single scattering albedo
271 ! tau - optical depth
272 ! g - asymmetry factor
273 !
274 ! important internal variables:
275 !
276 ! ks - scattering coeffcient (/mm)
277 ! ka - absorption coeffient (/mm)
278 ! kp - eigenvalue of two-stream approximation
279 ! y - = yr+iyi
280 !
281 ! remarks:
282 !
283 ! attributes:
284 ! language: f90
285 ! machine: ibm rs/6000 sp
286 !
287 !$$$
289 ! use kinds, only: r_kind
290 ! use constants, only: half, one, two, three, pi
291 implicit none
292 real(r_kind) yr,yi,ep_real,ep_imag
293 real(r_kind) frequency,a,h,f,ssalb_v,ssalb_h,tau_v,tau_h,gv,gh,k
294 real(r_kind) ks1,ks2,ks3,ks,kr1,kr2,kr3,kr,ki1,ki2,ki3,ki,ka
295 real(r_kind) fact1,fact2,fact3,fact4,fact5
296 k = two*pi/(300._r_kind/frequency)
297 yr = (ep_real - one)/(ep_real + two)
298 yi = - ep_imag/(ep_real + two)
299 fact1 = (one+two*f)**2
300 fact2 = one-f*yr
301 fact3 = (one-f)**4
302 fact4 = f*(k*a)**3
303 fact5 = one+two*f*yr
304 ks1 = k*sqrt(fact2/fact5)
305 ks2 = fact4*fact3/fact1
306 ks3 = (yr/fact2)**2
307 ks = ks1*ks2*ks3
308 kr1 = fact5/fact2
309 ! kr2 = two*fact4*fact3/fact1
310 kr2 = two*ks2
311 kr3 = two*yi*yr/(fact2**3)
312 kr = k*sqrt(kr1+kr2*kr3)
313 ki1 = three*f*yi/fact2**2
314 ! ki2 = two*fact4*fact3/fact1
315 ki2 = kr2
316 ! ki3 = (yr/fact2)**2
317 ki3 = ks3
318 ki = k**2/(two*kr)*(ki1+ki2*ki3)
319 ! ka = ki - ks
320 gv = half
321 gh = half
322 ssalb_v = ks / ki
323 if(ssalb_v .gt. 0.999_r_kind) ssalb_v = 0.999_r_kind
324 ssalb_h = ssalb_v
325 tau_v = two*ki*h
326 tau_h = tau_v
327 return
328 end subroutine snow_optic
329 subroutine soil_diel(freq,t_soil,vmc,rhob,rhos,sand,clay,esm)
331 !$$$ subprogram documentation block
332 ! . . . .
333 ! subprogram: soil_diel calculate the dielectric properties of soil
334 !
335 ! prgmmr: org: nesdis date: 2000-11-28
336 !
337 ! abstract: compute the dilectric constant of the bare soil
338 !
339 ! program history log:
340 !
341 ! input argument list:
342 !
343 ! theta - local zenith angle (degree)
344 ! frequency - frequency (ghz)
345 ! t_soil - soil temperature
346 ! vmc - volumetric moisture content (demensionless)
347 ! rhob - bulk volume density of the soil (1.18-1.12)
348 ! rhos - density of the solids (2.65 g.cm^3 for
349 ! solid soil material)
350 ! sand - sand fraction (sand + clay = 1.0)
351 ! clay - clay fraction
352 !
353 ! output argument list:
354 !
355 ! esm - dielectric constant for bare soil
356 !
357 ! important internal variables:
358 !
359 ! esof - the permittivity of free space
360 ! eswo - static dieletric constant
361 ! tauw - relaxation time of water
362 ! s - salinity
363 !
364 ! remarks:
365 !
366 ! attributes:
367 ! language: f90
368 ! machine: ibm rs/6000 sp
369 !
370 !$$$
372 ! use kinds, only: r_kind
373 ! use constants, only: zero, one, two, pi
374 implicit none
375 real(r_kind) esof
376 real(r_kind) f,a,b,tauw,freq,t_soil,vmc,rhob,rhos,sand,clay
377 real(r_kind) s,alpha,beta,ess,rhoef,t,eswi,eswo
378 complex(r_kind) ix,esm,esw,es1,es2
379 !
380 !
381 ! ix = dcmplx(zero, one)
382 ! s = zero
383 alpha = 0.65_r_kind
384 beta = 1.09_r_kind - 0.11_r_kind*sand + 0.18_r_kind*clay
385 ess = (1.01_r_kind + 0.44_r_kind*rhos)**2 - 0.062_r_kind !a2
386 rhoef = -1.645_r_kind + 1.939_r_kind*rhob - 0.020213_r_kind*sand + 0.01594_r_kind*clay !a4
387 t = t_soil - t_landem
388 f = freq*1.0e9_r_kind
389 ! the permittivity at the high frequency limit
390 eswi = 5.5_r_kind
391 ! the permittivity of free space (esof)
392 esof = 8.854e-12_r_kind
393 ! static dieletric constant (eswo)
394 eswo = 87.134_r_kind+(-1.949e-1_r_kind+(-1.276e-2_r_kind+2.491e-4_r_kind*t)*t)*t
395 ! a = one+(1.613e-5_r_kind*t-3.656e-3_r_kind+(3.210e-5_r_kind-4.232e-7_r_kind*s)*s)*s
396 ! eswo = eswo*a
397 ! relaxation time of water (tauw)
398 tauw = 1.1109e-10_r_kind+(-3.824e-12_r_kind+(6.938e-14_r_kind-5.096e-16_r_kind*t)*t)*t
399 ! b = one+(2.282e-5_r_kind*t-7.638e-4_r_kind+(-7.760e-6_r_kind+1.105e-8_r_kind*s)*s)*s
400 ! tauw = tauw*b
401 if (vmc .gt. zero) then
402 es1 = dcmplx(eswi, - rhoef *(rhos - rhob)/(two*pi*f*esof*rhos*vmc))
403 else
404 es1 = dcmplx(eswi, zero)
405 endif
406 es2 = dcmplx(eswo - eswi,zero)/dcmplx(one, f*tauw)
407 esw = es1 + es2
408 esm = one + (ess**alpha - one)*rhob/rhos + vmc**beta*esw**alpha - vmc
409 esm = esm**(one/alpha)
410 if(imag(esm) .ge.zero) esm = dcmplx(real(esm), -0.0001_r_kind)
411 return
412 end subroutine soil_diel
413 !
414 subroutine snow_diel(frequency,ep_real,ep_imag,rad,frac,ep_eff)
416 !$$$ subprogram documentation block
417 ! . . . .
418 ! subprogram: snow_diel compute dielectric constant of snow
419 !
420 ! prgmmr: org: nesdis date: 2000-11-28
421 !
422 ! abstract: compute dielectric constant of snow
423 !
424 !
425 ! program history log:
426 !
427 ! input argument list:
428 !
429 ! frequency - frequency (ghz)
430 ! ep_real - real part of dielectric constant of particle
431 ! ep_imag - imaginary part of dielectric constant of particle
432 ! rad - particle radiu (mm)
433 ! frac - fraction volume of snow (0.0 - 1.0)
434 !
435 ! output argument list:
436 !
437 ! ep_eff - dielectric constant of the dense medium
438 !
439 ! remarks:
440 !
441 ! attributes:
442 ! language: f90
443 ! machine: ibm rs/6000 sp
444 !
445 !$$$
447 ! use kinds, only: r_kind
448 ! use constants, only: zero, one, two, pi
449 implicit none
450 real(r_kind) ep_imag,ep_real
451 real(r_kind) frequency,rad,frac,k0,yr,yi
452 complex(r_kind) y,ep_r,ep_i,ep_eff,fracy
453 k0 = two*pi/(300.0_r_kind/frequency)
454 yr = (ep_real - one)/(ep_real + two)
455 yi = ep_imag/(ep_real + two)
456 y = dcmplx(yr, yi)
457 fracy=frac*y
458 ep_r = (one + two*fracy)/(one - fracy)
459 ep_i = two*fracy*y*(k0*rad)**3*(one-frac)**4/((one-fracy)**2*(one+two*frac)**2)
460 ep_eff = ep_r - dcmplx(zero,one)*ep_i
461 if (imag(ep_eff).ge.zero) ep_eff = dcmplx(real(ep_eff), -0.0001_r_kind)
462 return
463 end subroutine snow_diel
465 subroutine canopy_diel(frequency,mg,esv)
467 !
468 !$$$ subprogram documentation block
469 ! . . . .
470 ! subprogram: canopy_diel compute the dielectric constant of the vegetation canopy
471 !
472 ! prgmmr: org: nesdis date: 2000-11-28
473 !
474 ! abstract: compute the dielectric constant of the vegetation canopy
475 ! geomatrical optics approximation for vegetation canopy
476 ! work for horizontal leaves
477 !
478 ! program history log:
479 ! 2004-09-22 todling - replaced zsqrt by general interface sqrt
480 !
481 ! input argument list:
482 !
483 ! frequency - frequency (ghz)
484 ! mg - gravimetric water content
485 !
486 ! output argument list:
487 !
488 ! esv - dielectric constant of leaves
489 !
490 ! remarks:
491 !
492 ! references:
493 !
494 ! ulaby and el-rayer, 1987: microwave dielectric spectrum of vegetation part ii,
495 ! dual-dispersion model, ieee trans geosci. remote sensing, 25, 550-557
496 !
497 ! attributes:
498 ! language: f90
499 ! machine: ibm rs/6000 sp
500 !
501 !$$$
503 ! use kinds, only: r_kind
504 ! use constants, only: zero, one
505 implicit none
506 real(r_kind) frequency, mg, en, vf, vb
507 complex(r_kind) esv, xx
508 en = 1.7_r_kind - (0.74_r_kind - 6.16_r_kind*mg)*mg
509 vf = mg*(0.55_r_kind*mg - 0.076_r_kind)
510 vb = 4.64_r_kind*mg*mg/( one + 7.36_r_kind*mg*mg)
511 xx = dcmplx(zero,one)
512 esv = en + vf*(4.9_r_kind + 75.0_r_kind/(one + xx*frequency/18.0_r_kind)-xx*(18.0_r_kind/frequency)) + &
513 vb*(2.9_r_kind + 55.0_r_kind/(one + sqrt(xx*frequency/0.18_r_kind)))
514 if (imag(esv).ge.zero) esv = dcmplx(real(esv), -0.0001_r_kind)
515 return
516 end subroutine canopy_diel
517 subroutine reflectance(em1, em2, theta_i, theta_t, rvert, rh)
519 !$$$ subprogram documentation block
520 ! . . . .
521 ! subprogram: reflectance compute the surface reflectivity
522 !
523 ! prgmmr: org: nesdis date: 2000-11-28
524 !
525 ! abstract: compute the surface reflectivety using fresnel equations
526 ! for a rough surface having a standard deviation of height of sigma
527 !
528 ! program history log:
529 ! 2004-09-22 todling - replaced zsqrt by general interface sqrt; same for zabs/abs
530 !
531 ! input argument list:
532 ! theta_i - incident angle (degree)
533 ! theta_t - transmitted angle (degree)
534 ! em1 - dielectric constant of the medium 1
535 ! em2 - dielectric constant of the medium 2
536 !
537 ! output argument list:
538 !
539 ! rvert - reflectivity at vertical polarization
540 ! rh - reflectivity at horizontal polarization
541 !
542 ! remarks:
543 !
544 ! attributes:
545 ! language: f90
546 ! machine: ibm rs/6000 sp
547 !
548 !$$$
550 ! use kinds, only: r_kind
551 ! use constants, only: zero
552 implicit none
553 real(r_kind) theta_i, theta_t
554 real(r_kind) rh, rvert,cos_i,cos_t
555 complex(r_kind) em1, em2, m1, m2, angle_i, angle_t
556 ! compute the refractive index ratio between medium 2 and 1
557 ! using dielectric constant (n = sqrt(e))
558 cos_i=cos(theta_i)
559 cos_t=cos(theta_t)
560 angle_i = dcmplx(cos_i, zero)
561 angle_t = dcmplx(cos_t, zero)
562 m1 = sqrt(em1)
563 m2 = sqrt(em2)
564 rvert = (abs((m1*angle_t-m2*angle_i)/(m1*angle_t+m2*angle_i)))**2
565 rh = (abs((m1*angle_i-m2*angle_t)/(m1*angle_i+m2*angle_t)))**2
566 return
567 end subroutine reflectance
569 subroutine transmitance(em1,em2,theta_i,theta_t,tv,th)
571 !
572 !$$$ subprogram documentation block
573 ! . . . .
574 ! subprogram: transmitance calculate transmitance
575 !
576 ! prgmmr: org: nesdis date: 2000-11-28
577 !
578 ! abstract: compute transmitance
579 !
580 ! program history log:
581 ! 2004-09-22 todling - replaced zsqrt by general interface sqrt; same for zabs/abs
582 !
583 ! input argument list:
584 !
585 ! theta - local zenith angle (degree)
586 ! theta_i - incident angle (degree)
587 ! theta_t - transmitted angle (degree)
588 ! em1 - dielectric constant of the medium 1
589 ! em2 - dielectric constant of the medium 2
590 !
591 ! output argument list:
592 !
593 ! tv - transmisivity at vertical polarization
594 ! th - transmisivity at horizontal polarization
595 !
596 ! remarks:
597 !
598 ! attributes:
599 ! language: f90
600 ! machine: ibm rs/6000 sp
601 !
602 !$$$
604 ! use kinds, only: r_kind
605 ! use constants, only: zero, two
606 implicit none
607 real(r_kind) theta_i, theta_t
608 !
609 real(r_kind) th, tv, rr, cos_i,cos_t
610 !
611 complex(r_kind) em1, em2, m1, m2, angle_i, angle_t
612 ! compute the refractive index ratio between medium 2 and 1
613 ! using dielectric constant (n = sqrt(e))
614 cos_i=cos(theta_i)
615 cos_t=cos(theta_t)
616 angle_i = dcmplx(cos_i, zero)
617 angle_t = dcmplx(cos_t, zero)
618 m1 = sqrt(em1)
619 m2 = sqrt(em2)
620 rr = abs(m2/m1)*cos_t/cos_i
621 tv = rr*(abs(two*m1*angle_i/(m1*angle_t + m2*angle_i)))**2
622 th = rr*(abs(two*m1*angle_i/(m1*angle_i + m2*angle_t)))**2
623 return
624 end subroutine transmitance
626 subroutine rough_reflectance(frequency,theta,sigma,rvert,rh)
628 !$$$ subprogram documentation block
629 ! . . . .
630 ! subprogram: rought_reflectance calculate surface relectivity
631 !
632 ! prgmmr: org: nesdis date: 2000-11-28
633 !
634 ! abstract: compute the surface reflectivety for a rough surface
635 ! having a standard devoation of height of sigma
636 !
637 !
638 ! program history log:
639 !
640 ! input argument list:
641 !
642 ! frequency - frequency (ghz)
643 ! theta - local zenith angle (degree)
644 ! sigma - standard deviation of rough surface height
645 ! smooth surface:0.38, medium: 1.10, rough:2.15 cm
646 !
647 ! internal variables
648 !
649 !
650 ! output argument list:
651 !
652 ! rvert - reflectivity at vertical polarization
653 ! rh - reflectivity at horizontal polarization
654 !
655 !
656 ! important internal variables:
657 !
658 ! k0 - a propagation constant or wavenumber in a free space
659 !
660 ! remarks:
661 !
662 ! references:
663 !
664 ! wang, j. and b. j. choudhury, 1992: passive microwave radiation from soil: examples...
665 ! passive microwave remote sensing of .. ed. b. j. choudhury, etal vsp.
666 ! also wang and choudhury (1982)
667 !
668 ! attributes:
669 ! language: f90
670 ! machine: ibm rs/6000 sp
671 !
672 !$$$
674 ! use kinds, only: r_kind
675 ! use constants, only: one, two
676 implicit none
677 real(r_kind) theta, frequency
678 ! real(r_kind) p, q, sigma, rh, rvert, rh_s, rv_s
679 real(r_kind) p, q, rh, rvert, rh_s, rv_s, sigma
680 ! rh_s = rh
681 ! rv_s = rvert
682 rh_s = 0.3_r_kind*rh
683 rv_s = 0.3_r_kind*rvert
684 ! p = 0.3_r_kind
685 q = 0.35_r_kind*(one - exp(-0.60_r_kind*frequency*sigma**two))
686 ! rh = (q*rv_s + (one - q)*rh_s)*p
687 ! rv = (q*rh_s + (one - q)*rv_s)*p
688 rh = rh_s + q*(rv_s-rh_s)
689 rvert = rv_s + q*(rh_s-rv_s)
690 return
691 end subroutine rough_reflectance
693 subroutine two_stream_solution(b,mu,gv,gh,ssalb_h,ssalb_v,tau_h,tau_v,r12_h, &
694 r12_v,r21_h,r21_v,r23_h,r23_v,t21_v,t21_h,t12_v,t12_h,esv,esh)
696 !$$$ subprogram documentation block
697 ! . . . .
698 ! subprogram: two_stream_solution
699 !
700 ! prgmmr: org: nesdis date: 2000-11-28
701 !
702 ! abstract: two stream solution
703 !
704 ! version: beta
705 !
706 ! program history log:
707 !
708 ! input argument list:
709 !
710 ! b - scattering layer temperature (k) (gdas)
711 ! mu - cos(theta)
712 ! gv - asymmetry factor for v pol
713 ! gh - asymmetry factor for h pol
714 ! ssalb_v - single scattering albedo at v. polarization
715 ! ssalb_h - single scattering albedo at h. polarization
716 ! tau_v - optical depth at v. polarization
717 ! tau_h - optical depth at h. polarization
718 ! r12_v - reflectivity at vertical polarization
719 ! r12_h - reflectivity at horizontal polarization
720 ! r21_v - reflectivity at vertical polarization
721 ! r21_h - reflectivity at horizontal polarization
722 ! r23_v - reflectivity at vertical polarization
723 ! r23_h - reflectivity at horizontal polarization
724 ! t21_v - transmisivity at vertical polarization
725 ! t21_h - transmisivity at horizontal polarization
726 ! t12_v - transmisivity at vertical polarization
727 ! t12_h - transmisivity at horizontal polarization
728 !
729 ! output argument list:
730 !
731 ! esh - emissivity for horizontal polarization
732 ! esv - emissivity for vertical polarization
733 !
734 ! remarks:
735 !
736 ! attributes:
737 ! language: f90
738 ! machine: ibm rs/6000 sp
739 !
740 !$$$
742 ! use kinds, only: r_kind
743 ! use constants, only: one, two
744 implicit none
745 real(r_kind) b, mu, gv, gh, ssalb_h, ssalb_v, tau_h,tau_v,r12_h, &
746 r12_v, r21_h, r21_v, r23_h, r23_v, t21_v, t21_h, t12_v, t12_h, esv, esh
747 ! real(r_kind) esh0, esh1, esh2, esv0, esv1, esv2
748 ! real(r_kind) esh1, esv1
749 real(r_kind) alfa_v, alfa_h, kk_h, kk_v, gamma_h, gamma_v, beta_v, beta_h
750 real(r_kind) fact1,fact2
751 alfa_h = sqrt((one - ssalb_h)/(one - gh*ssalb_h))
752 kk_h = sqrt ((one - ssalb_h)*(one - gh*ssalb_h))/mu
753 beta_h = (one - alfa_h)/(one + alfa_h)
754 gamma_h = (beta_h -r23_h)/(one-beta_h*r23_h)
755 alfa_v = sqrt((one-ssalb_v)/(one - gv*ssalb_v))
756 kk_v = sqrt ((one-ssalb_v)*(one - gv*ssalb_v))/mu
757 beta_v = (one - alfa_v)/(one + alfa_v)
758 gamma_v = (beta_v -r23_v)/(one-beta_v*r23_v)
759 ! esh0 = i0/b*r12_h
760 ! esh0 = zero
761 ! esh1 = (one - beta_h)*(one + gamma_h*exp(-two*kk_h*tau_h))
762 ! esh1 = esh1/((one-beta_h*r21_h)-(beta_h-r21_h)*gamma_h*exp(-two*kk_h*tau_h))
763 ! esh2 = i0/b*t12_h*(beta_h-gamma_h*exp(-two*kk_h*tau_h))
764 ! esh2 = esh2 /((one-beta_h*r21_h)-(beta_h-r21_h)*gamma_h*exp(-two*kk_h*tau_h))
765 ! esh2 = zero
766 ! esv0 = i0/b*r12_v
767 ! esv0 = zero
768 ! esv1 = (one-beta_v)*(one+gamma_v*exp(-two*kk_v*tau_v))
769 ! esv1 = esv1/((one-beta_v*r21_v)-(beta_v-r21_v)*gamma_v*exp(-two*kk_v*tau_v))
770 ! esv2 = i0/b*t12_v*(beta_v - gamma_v*exp(-two*kk_v*tau_v))
771 ! esv2 = esv2 /((one-beta_v*r21_v)-(beta_v-r21_v)*beta_v*exp(-two*kk_v*tau_v))
772 ! esv2 = zero
773 ! esh = esh0 + t21_h*(esh1 + esh2)
774 ! esv = esv0 + t21_v*(esv1 + esv2)
775 fact1=gamma_h*exp(-two*kk_h*tau_h)
776 fact2=gamma_v*exp(-two*kk_v*tau_v)
777 esh = t21_h*(one - beta_h)*(one + fact1) &
778 /(one-beta_h*r21_h-(beta_h-r21_h)*fact1)
779 esv = t21_v*(one - beta_v)*(one + fact2) &
780 /(one-beta_v*r21_v-(beta_v-r21_v)*fact2)
781 return
782 end subroutine two_stream_solution