Merge remote-tracking branch 'origin/release-v4.5.2'
[WRF.git] / var / da / da_ssmi / tbatmos.inc
blobab534ee5241bf744f7d0e477c943bdadfcb5ba36
1 !==============================================================
3 !  SYNOPSIS OF subroutineS TB, TBATMOS
5 !===============================================================
6 !  subroutine TBATMOS(ifREQ,THETA,P0,WV,HWV,TA,GAMMA,LW,ZCLD,TBUP,
7 !                    TBDN,TAUATM)
9 !  This routine calculates the upwelling and downwelling microwave 
10 !  atmospheric brightness temperatures and transmittance at SSM/I
11 !  frequencies.  It is a generalization of "TBCLEAR" and includes a 
12 !  cloud layer of known height and total liquid water content.
13 !------------------------------------------
14 !  Input :   ifREQ = (1,2,3, or 4 ) for (19.35, 22.235, 37, or 85.5 GHz)
15 !            THETA = incidence angle (deg.)
17 !                                                       ( approximate range )
18 !            P0    = surface pressure (mb)                      (940 -1030)
19 !            WV    = precipitable water (kg/m**2)                 (0 - 70)
20 !            HWV   = water vapor density scale height (km)      (0.5 - 3.0)
21 !            TA, GAMMA = "effective" surface air temperature 
22 !                               and lapse rate (K ; km)
23 !                     T(z) = Ta - gamma*z              (263 - 303 ; 4.0 - 6.5)
25 !            LW, ZCLD  = total cloud liquid water content (kg/m**2) and height (km)
27 !  Output :  TBUP, TBDN = atmospheric component of upwelling and downwelling
28 !                         brightness temperature (K)
29 !            TAUATM    = total atmospheric transmittance at specified incidence angle.
31 !  Subroutines called : EFFHT, SIGMA_V
34 !===================================================================
36 !   NOTES
38 !===================================================================
41 !    1) This model is described in detail by Petty (1990).  
42 !    Frequency dependent
43 !    coefficients of the model are based on radiative transfer
44 !    integrations of 16,893 radiosonde profiles from 29 ships, islands and
45 !    coastal stations, representing all major climatic regimes and all
46 !    seasons.  Absorption coefficients used in these integrations were
47 !    obtained from an adaptation of the Millimeter-wave Propagation Model
48 !    (MPM) of Liebe (1985) with updated line parameters (Liebe 1988,
49 !    personal communication).  Empirical corrections to total water
50 !    absorption values have been added, based on comparisons with
51 !    radiosondes.
52 !    
53 !    2) The frequency-dependent component of the model is only 
54 !    intended to give meaningful results for
55 !    combinations of input parameters which are meteorologically and
56 !    physically reasonable under coastal or maritime conditions near sea
57 !    level.  The user must ensure that input values are not only
58 !    individually reasonable, but also are mutually consistent.  See
59 !    Petty (1990) for details.
60 !    
61 !    3) Computed brightness temperatures from this model are generally
62 !    valid only for (gaseous optical depth)/cos(theta) of order unity or
63 !    less.  This is usually of concern only for theta > 75 deg. in a moist
64 !    atmosphere at 22 and 85 GHz.  An ad hoc correction has been added to
65 !    TBCLEAR which greatly improves the calculated downwelling brightness
66 !    temperature for near-grazing angles, so that its contribution to the
67 !    diffuse reflected component from the sea surface may be correctly
68 !    estimated.  No such correction is applied to upwelling brightness
69 !    temperatures, since there is no contribution to SSM/I observations
70 !    from angles other than at the nadir angle.  The correction is also
71 !    not yet available for downwelling brightness temperatures from
72 !    TBATMOS. 
73 !                  
75 !    The author welcomes feedback from other workers concerning the
76 !    accuracy and overall performance of the routines presented here.
77 !    
78 !    
79 !    References:
80 !    
81 !       Liebe, H.J., 1985 :  An updated model for millimeter
82 !         wave propagation in moist air.  Radio Science, 20, pp. 1069-1089
83 !    
84 !       Liebe, H.J., 1988, personal communication: updated absorption line
85 !         parameters. 
86 !    
87 !       Petty, G.W. 1990:  On the Response of the SSM/I to the 
88 !         Marine Environment --- Implications for Atmospheric Parameter
89 !         Retrievals.  Ph.D. Dissertation, Dept. of Atmospheric Sciences,
90 !         University of Washington
92 !       Petty, G.W., and K.B. Katsaros, 1992: The response
93 !        of the Special Sensor Microwave/Imager to the
94 !        marine environment. Part I: An analytic model for
95 !        the atmospheric component of observed brightness
96 !        temperatures.  J. Atmos. Ocean. Tech., 9, 746--761
97 !       
98 !       Petty, G.W., and K.B. Katsaros, 1993: The response
99 !        of the Special Sensor Microwave/Imager to the
100 !        marine environment. Part II: A parameterization of
101 !        roughness effects on sea surface emission and
102 !        reflection.  Submitted to J. Atmos. Ocean. Tech.
103 !       
104 !    
105 !===============================================================
107 !  fortran source code listings for TBATMOS,
108 !                       SIGMA_V, EFFHT
110 !===============================================================
111 !  subroutine TBATMOS(ifREQ,THETA,P0,WV,HWV,TA,GAMMA,LW,ZCLD,TBUP,
112 !                    TBDN,TAUATM)
114 !  This routine calculates the upwelling and downwelling microwave 
115 !  atmospheric brightness temperatures and transmittance at an SSM/I     
116 !  frequency.  It is a generalization of "TBCLEAR" and includes a 
117 !  cloud layer of known height and total liquid water content.
119 !  Input :   ifREQ = (1,2,3, or 4 ) for (19.35, 22.235, 37, or 85.5 GHz)
120 !            THETA = incidence angle (deg.)
122 !                                                       ( approximate range )
123 !            P0    = surface pressure (mb)                      (940 -1030)
124 !            WV    = precipitable water (kg/m**2)                 (0 - 70)
125 !            HWV   = water vapor density scale height (km)      (0.5 - 3.0)
126 !            TA, GAMMA = "effective" surface air temperature 
127 !                               and lapse rate (K ; km)
128 !                     T(z) = Ta - gamma*z              (263 - 303 ; 4.0 - 6.5)
130 !            LW, ZCLD  = total cloud liquid water content (kg/m**2) and height (km)
132 !  Output :  TBUP, TBDN = atmospheric component of upwelling and downwelling
133 !                         brightness temperature (K)
134 !            TAUATM    = total atmospheric transmittance at specified incidence angle.
136 !  Subroutines called : EFFHT, SIGMA_V
137 !------------------------------------------------------------
138       subroutine tbatmos(ifreq,theta,p0,wv,hwv,ta,gamma,lw,zcld, &
139           tbup,tbdn,tauatm)
140       integer ifreq
141       real theta,p0,wv,hwv,ta,gamma,lw,zcld,tbup,tbdn,tauatm
142       real mu,hdn,hup,hdninf,hupinf
144       real b1(4),b2(4),b3(4)
145       real c(4),d1(4),d2(4),d3(4),zeta(4),kw0(4),kw1(4),kw2(4),kw3(4)
146       real tau,tau1,tau2,taucld
147       real tcld,tc,em,em1
148       real sigv,sigo,sig,sig1,sigcld
149       real teff1dn,teff1up,teffdn,teffup
150       real tbcld,tbclrdn,tbclrup,tb1dn,tb1up,tb2dn,tb2up
151       real otbar,tc2,tc3,hv,ho,alph
153       data b1/-.46847e-1,-.57752e-1,-.18885,.10990/
154       data b2/.26640e-4,.31662e-4,.9832e-4,.60531e-4/
155       data b3/.87560e+1,.10961e+2,.36678e+2,-.37578e+2/
156       data c/ .9207,   1.208,     .8253,     .8203/
157       data zeta/4.2,4.2,4.2,2.9/
158       data d1/-.35908e+1,-.38921e+1,-.43072e+1,-.17020e+0/
159       data d2/ .29797e-1, .31054e-1, .32801e-1, .13610e-1/
160       data d3/-.23174e-1,-.23543e-1,-.24101e-1,-.15776e+0/
161       data kw0/ .786e-1, .103,    .267,    .988/
162       data kw1/-.230e-2,-.296e-2,-.673e-2,-.107e-1/
163       data kw2/ .448e-4, .557e-4, .975e-4,-.535e-4/
164       data kw3/-.464e-6,-.558e-6,-.724e-6, .115e-5/
165 ! mu = secant(theta)
166       mu = 1.0/cos(theta*0.0174533)
167 ! get water vapor optical depth
168 !=====
169       call cal_sigma_v(ifreq,p0,wv,hwv,ta,gamma,sigv)
170 ! otbar = one over "mean" temperature
171       otbar = 1.0/(ta - gamma*zeta(ifreq))
172 ! sigo = dry air optical depth
173       sigo = b1(ifreq) + b2(ifreq)*p0  + b3(ifreq)*otbar
174 ! cloud parameters
175       tcld   = ta - gamma*zcld
176       tc = tcld - t_kelvin
177       tc2 = tc*tc
178       tc3 = tc2*tc
179       sigcld = (kw0(ifreq) + tc*kw1(ifreq) + tc2*kw2(ifreq) +  &
180            tc3*kw3(ifreq))*lw
181       taucld = exp(-mu*sigcld)
182       tbcld = (1.0 - taucld)*tcld
183 ! hv, ho = effective absorber scale heights for vapor, dry air
184       hv = c(ifreq)*hwv
185       ho = d1(ifreq) + d2(ifreq)*ta + d3(ifreq)*gamma
186 ! get effective emission heights for layer 1 and total atmosphere
187       call effht(ho,hv,sigo,sigv,mu,zcld,hdn,hup, &
188          hdninf,hupinf)
189 ! atmospheric transmittances in layer one and two, and combined
190       sig = sigo + sigv
191       sig1 = sigo*(1.0-exp(-zcld/ho)) + sigv*(1.0-exp(-zcld/hv))
192       tau  = exp(-mu*sig)
193       tau1 = exp(-mu*sig1)
194       tau2 = tau/tau1
195 ! atmospheric "emissivity"
196       em1  = 1.0 - tau1
197       em   = 1.0 - tau
198 ! downwelling and upwelling brightness temperature for each layer
199       teff1dn = ta - gamma*hdn
200       teff1up = ta - gamma*hup
201       teffdn = ta - gamma*hdninf
202       teffup = ta - gamma*hupinf
203       tbclrdn = teffdn*em
204       tbclrup = teffup*em
206       tb1dn = em1*teff1dn
207       tb1up = em1*teff1up
208       tb2dn = (tbclrdn - tb1dn)/tau1
209       tb2up =  tbclrup - tau2*tb1up
210 ! total downwelling and upwelling brightness temperature and transmittance
211       tbdn  = tb1dn + tau1*(tbcld + taucld*tb2dn)
212       tbup  = tb2up + tau2*(tbcld + taucld*tb1up)
213       tauatm = tau*taucld
215 ! the following lines apply an ad hoc correction to improve fit 
216 ! at large angles and/or high gaseous opacities 
217 !  (downwelling brightness temperatures only)
219       alph = (0.636619*atan(mu*sig))**2
220       tbdn = (1.0-alph)*tbdn + em*alph*ta
222       end subroutine tbatmos