da_tbatmos_adj.inc
References to this file elsewhere.
1 subroutine da_tbatmos_adj(ifreq,theta,p0,wv,hwv,ta,gamma,lw,zcld, &
2 tbup,tbdn,tauatm, &
3 ADJ_theta,ADJ_p0,ADJ_wv,ADJ_hwv,ADJ_ta,ADJ_gamma, &
4 ADJ_lw,ADJ_zcld,ADJ_tbup,ADJ_tbdn, &
5 ADJ_tauatm )
6
7 implicit none
8
9 !-----------------------------------------------------------------
10 ! Purpose: TBD
11 ! Output : ADJ_p0,ADJ_wv,ADJ_hwv,ADJ_ta,ADJ_gamma,ADJ_lw,ADJ_zcld
12 ! ADJ_theta (somtime theta is a variable)
13 ! Input : ADJ_tbup,ADJ_tbdn,ADJ_tauatm
14 ! Output mean fields : tbup,tbdn,tauatm
15 !-----------------------------------------------------------------
16
17 integer,intent(in ) :: ifreq
18 real ,intent(in ) :: theta,p0,wv,hwv,ta,gamma,lw,zcld
19 real ,intent(inout) :: ADJ_p0,ADJ_wv,ADJ_hwv,ADJ_ta, &
20 ADJ_gamma,ADJ_lw,ADJ_zcld,ADJ_theta
21 real ,intent(inout) :: ADJ_tbup,ADJ_tbdn,ADJ_tauatm
22 real, intent(out) :: tbup,tbdn,tauatm
23
24 real :: tbdn_save
25
26 real :: mu,hdn,hup,hdninf,hupinf,ADJ_mu
27
28 real :: b1(4),b2(4),b3(4)
29 real :: c(4),d1(4),d2(4),d3(4),zeta(4),kw0(4),kw1(4),kw2(4),kw3(4)
30 real :: tau,tau1,tau2,taucld
31 real :: tcld,tc,em,em1
32 real :: sigv,sigo,sig,sig1,sigcld
33 real :: teff1dn,teff1up,teffdn,teffup
34 real :: tbcld,tbclrdn,tbclrup,tb1dn,tb1up,tb2dn,tb2up
35 real :: otbar,tc2,tc3,hv,ho,alph
36 real :: ADJ_sigv,ADJ_otbar,ADJ_sigo,ADJ_tcld,ADJ_tc,ADJ_tc2,ADJ_tc3
37 real :: ADJ_sigcld,ADJ_taucld,ADJ_tbcld,ADJ_hv,ADJ_ho
38 real :: ADJ_hdn,ADJ_hup,ADJ_hdninf,ADJ_sig,ADJ_sig1,ADJ_tau,ADJ_tau1
39 real :: ADJ_tau2,ADJ_em1,ADJ_teff1dn,ADJ_hupinf,ADJ_em,ADJ_teff1up
40 real :: ADJ_teffdn,ADJ_teffup,ADJ_tbclrdn,ADJ_tbclrup,ADJ_tb1dn,ADJ_tb1up
41 real :: ADJ_tb2dn,ADJ_tb2up,ADJ_alph
42
43
44 data b1/-.46847e-1,-.57752e-1,-.18885,.10990/
45 data b2/.26640e-4,.31662e-4,.9832e-4,.60531e-4/
46 data b3/.87560e+1,.10961e+2,.36678e+2,-.37578e+2/
47 data c/ .9207, 1.208, .8253, .8203/
48 data zeta/4.2,4.2,4.2,2.9/
49 data d1/-.35908e+1,-.38921e+1,-.43072e+1,-.17020e+0/
50 data d2/ .29797e-1, .31054e-1, .32801e-1, .13610e-1/
51 data d3/-.23174e-1,-.23543e-1,-.24101e-1,-.15776e+0/
52 data kw0/ .786e-1, .103, .267, .988/
53 data kw1/-.230e-2,-.296e-2,-.673e-2,-.107e-1/
54 data kw2/ .448e-4, .557e-4, .975e-4,-.535e-4/
55 data kw3/-.464e-6,-.558e-6,-.724e-6, .115e-5/
56
57 mu=0.;hdn=0.;hup=0.;hdninf=0.;hupinf=0.;ADJ_mu=0.
58 tcld=0.;tc=0.;em=0.;em1=0.
59 sigv=0.;sigo=0.;sig=0.;sig1=0.;sigcld=0.
60 teff1dn=0.;teff1up=0.;teffdn=0.;teffup=0.
61 tbcld=0.;tbclrdn=0.;tbclrup=0.;tb1dn=0.;tb1up=0.;tb2dn=0.;tb2up=0.
62 otbar=0.;tc2=0.;tc3=0.;hv=0.;ho=0.;alph=0.
63 ADJ_sigv=0.;ADJ_otbar=0.;ADJ_sigo=0.;ADJ_tcld=0.;
64 ADJ_tc=0.;ADJ_tc2=0.;ADJ_tc3=0.
65 ADJ_sigcld=0.;ADJ_taucld=0.;ADJ_tbcld=0.;ADJ_hv=0.;ADJ_ho=0.
66 ADJ_hdn=0.;ADJ_hup=0.;ADJ_hdninf=0.;ADJ_sig=0.;ADJ_sig1=0.
67 ADJ_tau=0.;ADJ_tau1=0.
68 ADJ_tau2=0.;ADJ_em1=0.;ADJ_teff1dn=0.;ADJ_hupinf=0.;ADJ_em=0.
69 ADJ_teff1up=0.;ADJ_teffdn=0;ADJ_teffup=0.;ADJ_tbclrdn=0.
70 ADJ_tbclrup=0.;ADJ_tb1dn=0.;ADJ_tb1up=0.
71 ADJ_tb2dn=0.;ADJ_tb2up=0.;ADJ_alph=0.
72 tau=0.;tau1=0.;tau2=0.;taucld=0.
73 tcld=0.;tc=0.;em=0.;em1=0.
74 sigv=0.;sigo=0.;sig=0.;sig1=0.;sigcld=0.
75 teff1dn=0.;teff1up=0.;teffdn=0.;teffup=0.
76
77
78 ! mu = secant(theta)
79 ! somtime theta is a variable
80
81 mu = 1.0/cos(theta*0.0174533)
82
83 ! get water vapor optical depth
84
85 call cal_sigma_v(ifreq,p0,wv,hwv,ta,gamma,sigv)
86
87 ! otbar = one over "mean" temperature
88
89 otbar = 1.0/(ta - gamma*zeta(ifreq))
90
91 ! sigo = dry air optical depth
92
93 sigo = b1(ifreq) + b2(ifreq)* p0 + b3(ifreq)* otbar
94
95 ! cloud parameters
96
97 tcld = ta - gamma*zcld
98 tc = tcld - t_kelvin
99 tc2 = tc*tc
100 tc3 = tc2*tc
101 sigcld = ( kw0(ifreq) + tc*kw1(ifreq) + tc2*kw2(ifreq) + &
102 tc3*kw3(ifreq) )*lw
103 taucld = exp(-mu*sigcld)
104 tbcld = (1.0 - taucld)*tcld
105
106 ! hv, ho = effective absorber scale heights for vapor, dry air
107
108 hv = c(ifreq)* hwv
109 ho = d1(ifreq) + d2(ifreq)* ta + d3(ifreq)* gamma
110
111 ! get effective emission heights for layer 1 and total atmosphere
112
113
114 call effht(ho,hv,sigo,sigv,mu,zcld,hdn,hup, hdninf,hupinf)
115
116 ! atmospheric transmittances in layer one and two, and combined
117
118 sig = sigo + sigv
119 sig1 = sigo*(1.0-exp(-zcld/ho)) + sigv*(1.0-exp(-zcld/hv))
120 tau = exp(-mu*sig)
121 tau1 = exp(-mu*sig1)
122 tau2 = tau/tau1
123
124 ! atmospheric "emissivity"
125
126 em1 = 1.0 - tau1
127 em = 1.0 - tau
128
129 ! downwelling and upwelling brightness temperature for each layer
130
131 teff1dn = ta - gamma*hdn
132 teff1up = ta - gamma*hup
133 teffdn = ta - gamma*hdninf
134 teffup = ta - gamma*hupinf
135 tbclrdn = teffdn*em
136 tbclrup = teffup*em
137
138 tb1dn = em1*teff1dn
139 tb1up = em1*teff1up
140 tb2dn = (tbclrdn - tb1dn)/tau1
141 tb2up = tbclrup - tau2*tb1up
142
143 ! total downwelling and upwelling brightness temperature and transmittance
144
145 tbdn = tb1dn + tau1*(tbcld + taucld*tb2dn)
146 tbup = tb2up + tau2*(tbcld + taucld*tb1up)
147 tauatm = tau*taucld
148
149 ! the following lines apply an ad hoc correction to improve fit
150 ! at large angles and/or high gaseous opacities
151 ! (downwelling brightness temperatures only)
152
153 alph = (0.636619*atan(mu*sig))**2
154 tbdn_save = tbdn
155 tbdn = (1.0-alph)*tbdn + em*alph*ta
156
157
158 ! start
159
160 tbdn = tbdn_save
161
162 ADJ_alph = - ADJ_tbdn*tbdn
163 ADJ_em = ADJ_tbdn*alph*ta
164 ADJ_alph = em*ADJ_tbdn*ta + ADJ_alph
165 ADJ_ta = em*alph*ADJ_tbdn + ADJ_ta
166 ADJ_tbdn = (1.0-alph)*ADJ_tbdn
167
168 if (abs(sig) .gt. 0.) then
169 ADJ_mu = 2.*0.636619*0.636619*ADJ_alph*sig*atan(mu*sig)/(1.+mu*mu*sig*sig)
170 ADJ_sig = 2.*0.636619*0.636619*mu*ADJ_alph*atan(mu*sig)/(1.+mu*mu*sig*sig)
171 else
172 ADJ_mu = 0.
173 ADJ_sig = 0.
174 end if
175
176 ADJ_tau = ADJ_tauatm*taucld
177 ADJ_taucld = tau*ADJ_tauatm
178 ADJ_tb2up = ADJ_tbup
179 ADJ_tau2 = ADJ_tbup*(tbcld + taucld*tb1up)
180 ADJ_tbcld = tau2*ADJ_tbup
181 ADJ_taucld = tau2*ADJ_tbup*tb1up + ADJ_taucld
182 ADJ_tb1up = tau2*taucld*ADJ_tbup
183 ADJ_tb1dn = ADJ_tbdn
184 ADJ_tau1 = ADJ_tbdn*(tbcld + taucld*tb2dn)
185 ADJ_tbcld = tau1*ADJ_tbdn + ADJ_tbcld
186 ADJ_taucld = tau1*ADJ_tbdn*tb2dn + ADJ_taucld
187 ADJ_tb2dn = tau1*taucld*ADJ_tbdn
188
189 ADJ_tbclrup = ADJ_tb2up
190 ADJ_tau2 = - ADJ_tb2up*tb1up + ADJ_tau2
191 ADJ_tb1up = - tau2*ADJ_tb2up + ADJ_tb1up
192
193 ADJ_tbclrdn = ADJ_tb2dn/tau1
194 ADJ_tb1dn = - ADJ_tb2dn/tau1 + ADJ_tb1dn
195 ADJ_tau1 = - tb2dn*ADJ_tb2dn/tau1 + ADJ_tau1
196
197 ADJ_em1 = ADJ_tb1up*teff1up
198 ADJ_teff1up = em1*ADJ_tb1up
199
200 ADJ_em1 = ADJ_tb1dn*teff1dn + ADJ_em1
201 ADJ_teff1dn = em1*ADJ_tb1dn
202
203 ADJ_teffup = ADJ_tbclrup*em
204 ADJ_em = teffup*ADJ_tbclrup + ADJ_em
205
206 ADJ_teffdn = ADJ_tbclrdn*em
207 ADJ_em = teffdn*ADJ_tbclrdn + ADJ_em
208
209 ADJ_ta = ADJ_teffup + ADJ_ta
210 ADJ_gamma = - ADJ_teffup*hupinf + ADJ_gamma
211 ADJ_hupinf = - gamma*ADJ_teffup
212
213 ADJ_ta = ADJ_teffdn + ADJ_ta
214 ADJ_gamma = - ADJ_teffdn*hdninf + ADJ_gamma
215 ADJ_hdninf = - gamma*ADJ_teffdn
216
217 ADJ_ta = ADJ_teff1up + ADJ_ta
218 ADJ_gamma = - ADJ_teff1up*hup + ADJ_gamma
219 ADJ_hup = - gamma*ADJ_teff1up
220
221 ADJ_ta = ADJ_teff1dn + ADJ_ta
222 ADJ_gamma = - ADJ_teff1dn*hdn + ADJ_gamma
223 ADJ_hdn = - gamma*ADJ_teff1dn
224
225 ADJ_tau = - ADJ_em + ADJ_tau
226
227 ADJ_tau1 = - ADJ_em1 + ADJ_tau1
228
229 ADJ_tau = ADJ_tau2/tau1 + ADJ_tau
230 ADJ_tau1 = - tau2*ADJ_tau2/tau1 + ADJ_tau1
231
232 ADJ_sig1 = - mu*ADJ_tau1*tau1
233 ADJ_mu = - ADJ_tau1*sig1*tau1 + ADJ_mu
234
235 ADJ_mu = - ADJ_tau*sig*tau + ADJ_mu
236 ADJ_sig = - mu*ADJ_tau*tau + ADJ_sig
237
238 ADJ_sigo = ADJ_sig1*(1.0-exp(-zcld/ho))
239 ADJ_sigv = ADJ_sig1*(1.0-exp(-zcld/hv))
240 ADJ_zcld = sigo*ADJ_sig1/ho*exp(-zcld/ho) + ADJ_zcld
241 ADJ_ho = - sigo*zcld*ADJ_sig1/(ho*ho)*exp(-zcld/ho)
242 ADJ_zcld = sigv*ADJ_sig1/hv*exp(-zcld/hv) + ADJ_zcld
243 ADJ_hv = - sigv*zcld*ADJ_sig1/(hv*hv)*exp(-zcld/hv)
244
245 ADJ_sigo = ADJ_sig + ADJ_sigo
246 ADJ_sigv = ADJ_sig + ADJ_sigv
247
248 call da_effht_adj(ho,hv,sigo,sigv,mu,zcld,hdn,hup, &
249 hdninf,hupinf, &
250 ADJ_ho,ADJ_hv,ADJ_sigo,ADJ_sigv,ADJ_mu, &
251 ADJ_zcld,ADJ_hdn,ADJ_hup,ADJ_hdninf, &
252 ADJ_hupinf )
253
254 ADJ_ta = d2(ifreq)*ADJ_ho + ADJ_ta
255 ADJ_gamma = d3(ifreq)*ADJ_ho + ADJ_gamma
256 ADJ_hwv = c(ifreq)*ADJ_hv + ADJ_hwv
257
258 ADJ_taucld = - ADJ_tbcld*tcld + ADJ_taucld
259 ADJ_tcld = (1.0 - taucld)*ADJ_tbcld
260
261 ADJ_mu = - ADJ_taucld*sigcld*taucld + ADJ_mu
262 ADJ_sigcld = - mu*ADJ_taucld*taucld
263
264 ADJ_tc = ADJ_sigcld*kw1(ifreq)*lw
265 ADJ_tc2 = ADJ_sigcld*kw2(ifreq)*lw
266 ADJ_tc3 = ADJ_sigcld*kw3(ifreq)*lw
267 ADJ_lw = (kw0(ifreq)+tc*kw1(ifreq)+tc2*kw2(ifreq)+tc3*kw3(ifreq)) &
268 *ADJ_sigcld + ADJ_lw
269
270 ADJ_tc2 = ADJ_tc3*tc + ADJ_tc2
271 ADJ_tc = tc2*ADJ_tc3 + ADJ_tc
272 ADJ_tc = 2.*tc*ADJ_tc2 + ADJ_tc
273
274 ADJ_tcld = ADJ_tc + ADJ_tcld
275
276 ADJ_ta = ADJ_tcld + ADJ_ta
277 ADJ_gamma = - ADJ_tcld*zcld + ADJ_gamma
278 ADJ_zcld = - gamma*ADJ_tcld + ADJ_zcld
279
280 ADJ_p0 = b2(ifreq)*ADJ_sigo + ADJ_p0
281 ADJ_otbar = b3(ifreq)*ADJ_sigo
282
283 ADJ_ta = - otbar*otbar*ADJ_otbar + ADJ_ta
284 ADJ_gamma = otbar*otbar*ADJ_otbar*zeta(ifreq) + ADJ_gamma
285
286 call da_sigma_v_adj(ifreq,p0,wv,hwv,ta,gamma,sigv, &
287 ADJ_p0,ADJ_wv,ADJ_hwv,ADJ_ta, &
288 ADJ_gamma,ADJ_sigv )
289
290 ADJ_theta = mu*mu*0.0174533*ADJ_mu*sin(theta*0.0174533) + ADJ_theta
291
292
293 end subroutine da_tbatmos_adj
294
295