wrf svn trunk commit r4103
[wrffire.git] / wrfv2_fire / var / da / da_ssmi / effht.inc
blobc69b29e3c53747435d8386c57b1903abf039dfda
1 !******************************************
3 ! subroutine EFFHT(HO,HV,SIGO,SIGV,MU,ZCLD,HDN,HUP,HDNinF,HUPinF)
5 ! This subroutine returns "effective emitting heights" for an atmosphere
6 ! with gaseous absorption approximated by two exponentially-decaying
7 ! profiles with distinct scale heights.  For applications at
8 ! SSM/I frequencies (i.e., 19, 22, 37, 86 GHz), these absorption profiles 
9 ! correspond to water vapor and dry air absorption, respectively.  
11 ! Input parameters: HO, HV  --  scale heights of absorption corresponding to 
12 !                               dry air and water vapor (km):
14 !    SIGO, SIGV --  total optical depth due to each constituent
15 !                      (the present model is valid only for (sigo+sigv < 1.0) )
16
17 !         MU    --  secant(theta), where theta is the path angle from vertical
19 !        ZCLD   --  upper limit (km) of the atmospheric layer for which
20 !                   hdn, hup are to be calculated.  Layer extends down to the
21 !                   surface.  These parameters permit the separate calculation
22 !                   brightness temperature contributions from a lower and
23 !                   upper atmospheric layer, separated by a cloud layer at
24 !                   height zcld.
26 ! Output parameters:  
28 !  HDN, HUP     --  "effective emitting height" of an atmospheric layer
29 !            bounded below by the surface and above at height zcld.
30 !            "Effective emitting temperature" for upward and downward
31 !            microwave radiation emitted by that layer is then given by
32 !            (Ta - gamma*hup) and (Ta - gamma*hdn), respectively, where
33 !            Ta and gamma are effective surface temperature (deg K)
34 !            and lapse rate (K / km), respectively.  Brightness temperatures
35 !            due to emission from this layer are then (1 - tau)(Ta - gamma*hup)
36 !            and (1 - tau)(Ta - gamma*hdn), respectively, where tau is
37 !            the slant path transmittance of the layer at angle theta
39 ! HDinF,HUPinF  -- same as hdn and hup, but for the entire
40 !                  atmosphere  (i.e., zcld set to infinity)
43 !      This version of EFFHT executes 83 floating point multiplications,
44 !      11 floating point divisions, and 2 calls to exp() .
46 !----------------------------------------------
47       subroutine effht(ho,hv,sigo,sigv,mu,zcld,hdn,hup,hdninf,hupinf)
49       real ho,hv,sigo,sigv,mu,zcld,gint,zgint,hint,zhint, &
50            ginf,zginf,hinf,zhinf,hdn,hup,hdninf,hupinf
52       real   hoinv,hvinv,chio,chiv,ezho,ezhv,alpha,alph2,alph3
53       real   beta,beta2,beta3,mu2,mualph,cplus,cmin,dplus,dmin
54       real   chiov,chivv,chioo,chioov,chiovv,chiooo,chivvv
55       real   h11,h21,h12
56       real   sigoo,sigov,sigvv,sigooo,sigoov,sigovv,sigvvv
57       real   ezhoo,ezhov,ezhvv,ezhooo,ezhoov,ezhovv,ezhvvv
58       real   s,sprim,t,tprim,u,uprim,term1,term2,term3
59       real   halfmu,halfmu2,sixthmu2,etnthmu2,quartmu
62       hoinv = 1.0d0/ho
63       hvinv = 1.0d0/hv
64       chio = zcld*hoinv
65       chiv = zcld*hvinv
66       ezho = sigo*exp(-chio)
67       ezhv = sigv*exp(-chiv)
68       alpha = sigo + sigv
69       alph2 = alpha*alpha
70       alph3 = alpha*alph2
71       beta =  ezho + ezhv
72       beta2 = beta*beta
73       beta3 = beta*beta2
75       mu2 = mu*mu
76       halfmu = 0.5d0*mu
77       sixthmu2 = mu2/6.0d0
78       etnthmu2 = mu2/18.0d0
79       quartmu = 0.25d0*mu
80       halfmu2 = 0.5d0*mu2
81       mualph = mu*alpha
82       cplus = 1.0d0 + mualph
83       cmin  = 1.0d0 - mualph
84       dplus = halfmu2*alph2
85       dmin  = dplus
86       dplus = cplus + dplus
87       dmin  = cmin  + dmin
89       h11 = hoinv + hvinv
90       h21 = 1.0d0/(h11 + hvinv)
91       h12 = 1.0d0/(h11 + hoinv)
92       h11 = 1.0d0/h11
93       chiov = 1.0d0 + chio + chiv
94       chioo = 1.0d0 + chio + chio
95       chivv = 1.0d0 + chiv + chiv
96       chioov = chioo + chiv
97       chiovv = chio + chivv
98       chiooo = chioo + chio
99       chivvv = chivv + chiv
100       chio = 1.0d0 + chio
101       chiv = 1.0d0 + chiv
102       sigov = sigo*sigv
103       sigoo = sigo*sigo
104       sigvv = sigv*sigv
105       sigooo = sigoo*sigo
106       sigoov = sigoo*sigv
107       sigovv = sigo*sigvv
108       sigvvv = sigvv*sigv
109       ezhoo = ezho*ezho
110       ezhov = ezho*ezhv
111       ezhvv = ezhv*ezhv
112       ezhovv = ezho*ezhvv
113       ezhoov = ezhoo*ezhv
114       ezhooo = ezhoo*ezho
115       ezhvvv = ezhvv*ezhv
116       s     = sigo*ho   +    sigv*hv
117       sprim = ezho*ho*chio + ezhv*hv*chiv
118       t     = sigoo*ho    +    4.0d0*sigov*h11     +   sigvv*hv
119       tprim = ezhoo*ho*chioo + 4.0d0*ezhov*h11*chiov + ezhvv*hv*chivv
120       u     = sigooo*ho+9.0d0*(sigovv*h21+sigoov*h12)+sigvvv*hv
121       uprim = ezhvvv*hv*chivvv +  &
122            9.0d0*(ezhovv*h21*chiovv + ezhoov*h12*chioov) + &
123            ezhooo*ho*chiooo
125       term1 = s - sprim
126       term2 = quartmu*(t - tprim)
127       term3 = etnthmu2*(u - uprim)
128       zgint = dmin*term1 +  cmin*term2 + term3
129       zhint = -dplus*term1 + cplus*term2 - term3
130       term2 = quartmu*t
131       term3 = etnthmu2*u
132       zginf = dmin*s +  cmin*term2 + term3
133       zhinf = -dplus*s + cplus*term2 - term3
134       term1 = alpha - beta
135       term2 = halfmu*(alph2 - beta2)
136       term3 = sixthmu2*(alph3 - beta3)
137       gint  = dmin*term1 +  cmin*term2 + term3
138       hint  = -dplus*term1 + cplus*term2 - term3
139       term2 = halfmu*alph2
140       term3 = sixthmu2*alph3
141       ginf  = dmin*alpha +  cmin*term2 + term3
142       hinf  = -dplus*alpha + cplus*term2 - term3
143       hdn   = zgint/gint
144       hup   = zhint/hint
145       hdninf = zginf/ginf
146       hupinf = zhinf/hinf
148       end subroutine effht