module_ra_gsfcsw.F
References to this file elsewhere.
1 !Comment the following out to turn off aerosol-radiation
2 !feedback between MOSAIC and GSFCSW. wig, 21-Feb-2005
3
4 MODULE module_ra_gsfcsw
5
6 REAL, PARAMETER, PRIVATE :: thresh=1.e-9
7 REAL, SAVE :: center_lat
8
9 ! Assign co2 and trace gases amount (units are parts/part by volumn)
10
11 REAL, PARAMETER, PRIVATE :: co2 = 300.e-6
12
13 CONTAINS
14
15 !------------------------------------------------------------------
16 ! urban related variable are added to arguments of gsfcswrad
17 !------------------------------------------------------------------
18 SUBROUTINE GSFCSWRAD(rthraten,gsw,xlat,xlong &
19 ,dz8w,rho_phy &
20 ,alb,t3d,qv3d,qc3d,qr3d &
21 ,qi3d,qs3d,qg3d,qndrop3d &
22 ,p3d,p8w3d,pi3d,cldfra3d,rswtoa &
23 ,gmt,cp,g,julday,xtime,declin,solcon &
24 ,radfrq,degrad,taucldi,taucldc,warm_rain &
25 ,tauaer300,tauaer400,tauaer600,tauaer999 & ! jcb
26 ,gaer300,gaer400,gaer600,gaer999 & ! jcb
27 ,waer300,waer400,waer600,waer999 & ! jcb
28 ,aer_ra_feedback &
29 ,f_qv,f_qc,f_qr,f_qi,f_qs,f_qg,f_qndrop &
30 ,ids,ide, jds,jde, kds,kde &
31 ,ims,ime, jms,jme, kms,kme &
32 ,its,ite, jts,jte, kts,kte &
33 ,cosz_urb2d,omg_urb2d ) !Optional urban
34 !------------------------------------------------------------------
35 IMPLICIT NONE
36 !------------------------------------------------------------------
37 INTEGER, PARAMETER :: np = 75
38
39 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
40 ims,ime, jms,jme, kms,kme, &
41 its,ite, jts,jte, kts,kte
42 LOGICAL, INTENT(IN ) :: warm_rain
43
44 INTEGER, INTENT(IN ) :: JULDAY
45
46
47 REAL, INTENT(IN ) :: RADFRQ,DEGRAD, &
48 XTIME,DECLIN,SOLCON
49 !
50 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
51 INTENT(IN ) :: P3D, &
52 P8W3D, &
53 pi3D, &
54 T3D, &
55 dz8w, &
56 rho_phy, &
57 CLDFRA3D
58
59
60 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
61 INTENT(INOUT) :: RTHRATEN
62 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
63 OPTIONAL, &
64 INTENT(INOUT) :: taucldi, &
65 taucldc
66 !
67 REAL, DIMENSION( ims:ime, jms:jme ), &
68 INTENT(IN ) :: XLAT, &
69 XLONG, &
70 ALB
71 !
72 REAL, DIMENSION( ims:ime, jms:jme ), &
73 INTENT(INOUT) :: GSW, &
74 RSWTOA
75 !
76 REAL, INTENT(IN ) :: GMT,CP,G
77 !
78
79 !
80 ! Optional
81 !
82 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
83 INTENT(IN ) :: tauaer300,tauaer400,tauaer600,tauaer999, & ! jcb
84 gaer300,gaer400,gaer600,gaer999, & ! jcb
85 waer300,waer400,waer600,waer999 ! jcb
86
87 INTEGER, INTENT(IN ), OPTIONAL :: aer_ra_feedback
88
89 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
90 OPTIONAL, &
91 INTENT(IN ) :: &
92 QV3D, &
93 QC3D, &
94 QR3D, &
95 QI3D, &
96 QS3D, &
97 QG3D, &
98 QNDROP3D
99
100 LOGICAL, OPTIONAL, INTENT(IN ) :: &
101 F_QV,F_QC,F_QR,F_QI,F_QS,F_QG, &
102 F_QNDROP
103
104 ! LOCAL VARS
105
106 REAL, DIMENSION( its:ite ) :: &
107 ts, &
108 cosz, &
109 fp, &
110 rsuvbm, &
111 rsuvdf, &
112 rsirbm, &
113 rsirdf, &
114 p400, &
115 p700
116
117 INTEGER, DIMENSION( its:ite ) :: &
118 ict, &
119 icb
120
121 REAL, DIMENSION( its:ite, kts-1:kte, 2 ) :: taucld
122
123 REAL, DIMENSION( its:ite, kts-1:kte+1 ) :: flx, &
124 flxd
125 !
126 REAL, DIMENSION( its:ite, kts-1:kte ) :: O3
127 !
128 REAL, DIMENSION( its:ite, kts-1:kte, 11 ) :: &
129 taual, &
130 ssaal, &
131 asyal
132
133 REAL, DIMENSION( its:ite, kts-1:kte, 2 ) :: &
134 reff, &
135 cwc
136 REAL, DIMENSION( its: ite, kts-1:kte+1 ) :: &
137 P8W2D
138 REAL, DIMENSION( its: ite, kts-1:kte ) :: &
139 TTEN2D, &
140 qndrop2d, &
141 SH2D, &
142 P2D, &
143 T2D, &
144 fcld2D
145 real, DIMENSION( its:ite , kts:kte+1 ) :: phyd
146 real, DIMENSION( its:ite , kts:kte ) :: phydmid
147
148 REAL, DIMENSION( np, 5 ) :: pres, &
149 ozone
150 REAL, DIMENSION( np ) :: p
151
152 LOGICAL :: cldwater,overcast, predicate
153 !
154 INTEGER :: i,j,K,NK,ib,kk,mix,mkx
155
156 ! iprof = 1 : mid-latitude summer profile
157 ! = 2 : mid-latitude winter profile
158 ! = 3 : sub-arctic summer profile
159 ! = 4 : sub-arctic winter profile
160 ! = 5 : tropical profile
161 !
162
163 INTEGER :: iprof, &
164 is_summer, &
165 ie_summer, &
166 lattmp
167
168
169 !
170 REAL :: XLAT0,XLONG0
171 REAL :: fac,latrmp
172 REAL :: xt24,tloctm,hrang,xxlat
173
174 !URBAN
175 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: COSZ_URB2D !urban
176 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: OMG_URB2D !urban
177
178 real, dimension(11) :: midbands ! jcb
179 data midbands/.2,.235,.27,.2875,.3025,.305,.3625,.55,1.92,1.745,6.135/ ! jcb
180 real :: ang,slope ! jcb
181 character(len=200) :: msg !wig
182 real pi, third, relconst, lwpmin, rhoh2o
183 !
184 !--------------------------------------------------------------------------------
185 ! data set 1
186 ! mid-latitude summer (75 levels) : p(mb) o3(g/g)
187 ! surface temp = 294.0
188 !
189 data (pres(i,1),i=1,np)/ &
190 0.0006244, 0.0008759, 0.0012286, 0.0017234, 0.0024174, &
191 0.0033909, 0.0047565, 0.0066720, 0.0093589, 0.0131278, &
192 0.0184145, 0.0258302, 0.0362323, 0.0508234, 0.0712906, &
193 0.1000000, 0.1402710, 0.1967600, 0.2759970, 0.3871430, &
194 0.5430, 0.7617, 1.0685, 1.4988, 2.1024, 2.9490, &
195 4.1366, 5.8025, 8.1392, 11.4170, 16.0147, 22.4640, &
196 31.5105, 44.2001, 62.0000, 85.7750, 109.5500, 133.3250, &
197 157.1000, 180.8750, 204.6500, 228.4250, 252.2000, 275.9750, &
198 299.7500, 323.5250, 347.3000, 371.0750, 394.8500, 418.6250, &
199 442.4000, 466.1750, 489.9500, 513.7250, 537.5000, 561.2750, &
200 585.0500, 608.8250, 632.6000, 656.3750, 680.1500, 703.9250, &
201 727.7000, 751.4750, 775.2500, 799.0250, 822.8000, 846.5750, &
202 870.3500, 894.1250, 917.9000, 941.6750, 965.4500, 989.2250, &
203 1013.0000/
204 !
205 data (ozone(i,1),i=1,np)/ &
206 0.1793E-06, 0.2228E-06, 0.2665E-06, 0.3104E-06, 0.3545E-06, &
207 0.3989E-06, 0.4435E-06, 0.4883E-06, 0.5333E-06, 0.5786E-06, &
208 0.6241E-06, 0.6698E-06, 0.7157E-06, 0.7622E-06, 0.8557E-06, &
209 0.1150E-05, 0.1462E-05, 0.1793E-05, 0.2143E-05, 0.2512E-05, &
210 0.2902E-05, 0.3313E-05, 0.4016E-05, 0.5193E-05, 0.6698E-05, &
211 0.8483E-05, 0.9378E-05, 0.9792E-05, 0.1002E-04, 0.1014E-04, &
212 0.9312E-05, 0.7834E-05, 0.6448E-05, 0.5159E-05, 0.3390E-05, &
213 0.1937E-05, 0.1205E-05, 0.8778E-06, 0.6935E-06, 0.5112E-06, &
214 0.3877E-06, 0.3262E-06, 0.2770E-06, 0.2266E-06, 0.2020E-06, &
215 0.1845E-06, 0.1679E-06, 0.1519E-06, 0.1415E-06, 0.1317E-06, &
216 0.1225E-06, 0.1137E-06, 0.1055E-06, 0.1001E-06, 0.9487E-07, &
217 0.9016E-07, 0.8641E-07, 0.8276E-07, 0.7930E-07, 0.7635E-07, &
218 0.7347E-07, 0.7065E-07, 0.6821E-07, 0.6593E-07, 0.6368E-07, &
219 0.6148E-07, 0.5998E-07, 0.5859E-07, 0.5720E-07, 0.5582E-07, &
220 0.5457E-07, 0.5339E-07, 0.5224E-07, 0.5110E-07, 0.4999E-07/
221
222 !--------------------------------------------------------------------------------
223 ! data set 2
224 ! mid-latitude winter (75 levels) : p(mb) o3(g/g)
225 ! surface temp = 272.2
226 !
227 data (pres(i,2),i=1,np)/ &
228 0.0006244, 0.0008759, 0.0012286, 0.0017234, 0.0024174, &
229 0.0033909, 0.0047565, 0.0066720, 0.0093589, 0.0131278, &
230 0.0184145, 0.0258302, 0.0362323, 0.0508234, 0.0712906, &
231 0.1000000, 0.1402710, 0.1967600, 0.2759970, 0.3871430, &
232 0.5430, 0.7617, 1.0685, 1.4988, 2.1024, 2.9490, &
233 4.1366, 5.8025, 8.1392, 11.4170, 16.0147, 22.4640, &
234 31.5105, 44.2001, 62.0000, 85.9000, 109.8000, 133.7000, &
235 157.6000, 181.5000, 205.4000, 229.3000, 253.2000, 277.1000, &
236 301.0000, 324.9000, 348.8000, 372.7000, 396.6000, 420.5000, &
237 444.4000, 468.3000, 492.2000, 516.1000, 540.0000, 563.9000, &
238 587.8000, 611.7000, 635.6000, 659.5000, 683.4000, 707.3000, &
239 731.2000, 755.1000, 779.0000, 802.9000, 826.8000, 850.7000, &
240 874.6000, 898.5000, 922.4000, 946.3000, 970.2000, 994.1000, &
241 1018.0000/
242 !
243 data (ozone(i,2),i=1,np)/ &
244 0.2353E-06, 0.3054E-06, 0.3771E-06, 0.4498E-06, 0.5236E-06, &
245 0.5984E-06, 0.6742E-06, 0.7511E-06, 0.8290E-06, 0.9080E-06, &
246 0.9881E-06, 0.1069E-05, 0.1152E-05, 0.1319E-05, 0.1725E-05, &
247 0.2145E-05, 0.2581E-05, 0.3031E-05, 0.3497E-05, 0.3980E-05, &
248 0.4478E-05, 0.5300E-05, 0.6725E-05, 0.8415E-05, 0.1035E-04, &
249 0.1141E-04, 0.1155E-04, 0.1143E-04, 0.1093E-04, 0.1060E-04, &
250 0.9720E-05, 0.8849E-05, 0.7424E-05, 0.6023E-05, 0.4310E-05, &
251 0.2820E-05, 0.1990E-05, 0.1518E-05, 0.1206E-05, 0.9370E-06, &
252 0.7177E-06, 0.5450E-06, 0.4131E-06, 0.3277E-06, 0.2563E-06, &
253 0.2120E-06, 0.1711E-06, 0.1524E-06, 0.1344E-06, 0.1199E-06, &
254 0.1066E-06, 0.9516E-07, 0.8858E-07, 0.8219E-07, 0.7598E-07, &
255 0.6992E-07, 0.6403E-07, 0.5887E-07, 0.5712E-07, 0.5540E-07, &
256 0.5370E-07, 0.5214E-07, 0.5069E-07, 0.4926E-07, 0.4785E-07, &
257 0.4713E-07, 0.4694E-07, 0.4676E-07, 0.4658E-07, 0.4641E-07, &
258 0.4634E-07, 0.4627E-07, 0.4619E-07, 0.4612E-07, 0.4605E-07/
259
260
261 !--------------------------------------------------------------------------------
262 ! data set 3
263 ! sub-arctic summer (75 levels) : p(mb) o3(g/g)
264 ! surface temp = 287.0
265 !
266 data (pres(i,3),i=1,np)/ &
267 0.0006244, 0.0008759, 0.0012286, 0.0017234, 0.0024174, &
268 0.0033909, 0.0047565, 0.0066720, 0.0093589, 0.0131278, &
269 0.0184145, 0.0258302, 0.0362323, 0.0508234, 0.0712906, &
270 0.1000000, 0.1402710, 0.1967600, 0.2759970, 0.3871430, &
271 0.5430, 0.7617, 1.0685, 1.4988, 2.1024, 2.9490, &
272 4.1366, 5.8025, 8.1392, 11.4170, 16.0147, 22.4640, &
273 31.5105, 44.2001, 62.0000, 85.7000, 109.4000, 133.1000, &
274 156.8000, 180.5000, 204.2000, 227.9000, 251.6000, 275.3000, &
275 299.0000, 322.7000, 346.4000, 370.1000, 393.8000, 417.5000, &
276 441.2000, 464.9000, 488.6000, 512.3000, 536.0000, 559.7000, &
277 583.4000, 607.1000, 630.8000, 654.5000, 678.2000, 701.9000, &
278 725.6000, 749.3000, 773.0000, 796.7000, 820.4000, 844.1000, &
279 867.8000, 891.5000, 915.2000, 938.9000, 962.6000, 986.3000, &
280 1010.0000/
281 !
282 data (ozone(i,3),i=1,np)/ &
283 0.1728E-06, 0.2131E-06, 0.2537E-06, 0.2944E-06, 0.3353E-06, &
284 0.3764E-06, 0.4176E-06, 0.4590E-06, 0.5006E-06, 0.5423E-06, &
285 0.5842E-06, 0.6263E-06, 0.6685E-06, 0.7112E-06, 0.7631E-06, &
286 0.1040E-05, 0.1340E-05, 0.1660E-05, 0.2001E-05, 0.2362E-05, &
287 0.2746E-05, 0.3153E-05, 0.3762E-05, 0.4988E-05, 0.6518E-05, &
288 0.8352E-05, 0.9328E-05, 0.9731E-05, 0.8985E-05, 0.7632E-05, &
289 0.6814E-05, 0.6384E-05, 0.5718E-05, 0.4728E-05, 0.4136E-05, &
290 0.3033E-05, 0.2000E-05, 0.1486E-05, 0.1121E-05, 0.8680E-06, &
291 0.6474E-06, 0.5164E-06, 0.3921E-06, 0.2996E-06, 0.2562E-06, &
292 0.2139E-06, 0.1723E-06, 0.1460E-06, 0.1360E-06, 0.1267E-06, &
293 0.1189E-06, 0.1114E-06, 0.1040E-06, 0.9678E-07, 0.8969E-07, &
294 0.8468E-07, 0.8025E-07, 0.7590E-07, 0.7250E-07, 0.6969E-07, &
295 0.6694E-07, 0.6429E-07, 0.6208E-07, 0.5991E-07, 0.5778E-07, &
296 0.5575E-07, 0.5403E-07, 0.5233E-07, 0.5067E-07, 0.4904E-07, &
297 0.4721E-07, 0.4535E-07, 0.4353E-07, 0.4173E-07, 0.3997E-07/
298
299
300 !--------------------------------------------------------------------------------
301 ! data set 3
302 ! sub-arctic winter (75 levels) : p(mb) o3(g/g)
303 ! surface temp = 257.1
304 !
305 data (pres(i,4),i=1,np)/ &
306 0.0006244, 0.0008759, 0.0012286, 0.0017234, 0.0024174, &
307 0.0033909, 0.0047565, 0.0066720, 0.0093589, 0.0131278, &
308 0.0184145, 0.0258302, 0.0362323, 0.0508234, 0.0712906, &
309 0.1000000, 0.1402710, 0.1967600, 0.2759970, 0.3871430, &
310 0.5430, 0.7617, 1.0685, 1.4988, 2.1024, 2.9490, &
311 4.1366, 5.8025, 8.1392, 11.4170, 16.0147, 22.4640, &
312 31.5105, 44.2001, 62.0000, 85.7750, 109.5500, 133.3250, &
313 157.1000, 180.8750, 204.6500, 228.4250, 252.2000, 275.9750, &
314 299.7500, 323.5250, 347.3000, 371.0750, 394.8500, 418.6250, &
315 442.4000, 466.1750, 489.9500, 513.7250, 537.5000, 561.2750, &
316 585.0500, 608.8250, 632.6000, 656.3750, 680.1500, 703.9250, &
317 727.7000, 751.4750, 775.2500, 799.0250, 822.8000, 846.5750, &
318 870.3500, 894.1250, 917.9000, 941.6750, 965.4500, 989.2250, &
319 1013.0000/
320 !
321 data (ozone(i,4),i=1,np)/ &
322 0.2683E-06, 0.3562E-06, 0.4464E-06, 0.5387E-06, 0.6333E-06, &
323 0.7301E-06, 0.8291E-06, 0.9306E-06, 0.1034E-05, 0.1140E-05, &
324 0.1249E-05, 0.1360E-05, 0.1474E-05, 0.1855E-05, 0.2357E-05, &
325 0.2866E-05, 0.3383E-05, 0.3906E-05, 0.4437E-05, 0.4975E-05, &
326 0.5513E-05, 0.6815E-05, 0.8157E-05, 0.1008E-04, 0.1200E-04, &
327 0.1242E-04, 0.1250E-04, 0.1157E-04, 0.1010E-04, 0.9063E-05, &
328 0.8836E-05, 0.8632E-05, 0.8391E-05, 0.7224E-05, 0.6054E-05, &
329 0.4503E-05, 0.3204E-05, 0.2278E-05, 0.1833E-05, 0.1433E-05, &
330 0.9996E-06, 0.7440E-06, 0.5471E-06, 0.3944E-06, 0.2852E-06, &
331 0.1977E-06, 0.1559E-06, 0.1333E-06, 0.1126E-06, 0.9441E-07, &
332 0.7678E-07, 0.7054E-07, 0.6684E-07, 0.6323E-07, 0.6028E-07, &
333 0.5746E-07, 0.5468E-07, 0.5227E-07, 0.5006E-07, 0.4789E-07, &
334 0.4576E-07, 0.4402E-07, 0.4230E-07, 0.4062E-07, 0.3897E-07, &
335 0.3793E-07, 0.3697E-07, 0.3602E-07, 0.3506E-07, 0.3413E-07, &
336 0.3326E-07, 0.3239E-07, 0.3153E-07, 0.3069E-07, 0.2987E-07/
337
338 !--------------------------------------------------------------------------------
339 ! data set 4
340 ! tropical (75 levels) : p(mb) o3(g/g)
341 ! surface temp = 300.0
342 !
343 data (pres(i,5),i=1,np)/ &
344 0.0006244, 0.0008759, 0.0012286, 0.0017234, 0.0024174, &
345 0.0033909, 0.0047565, 0.0066720, 0.0093589, 0.0131278, &
346 0.0184145, 0.0258302, 0.0362323, 0.0508234, 0.0712906, &
347 0.1000000, 0.1402710, 0.1967600, 0.2759970, 0.3871430, &
348 0.5430, 0.7617, 1.0685, 1.4988, 2.1024, 2.9490, &
349 4.1366, 5.8025, 8.1392, 11.4170, 16.0147, 22.4640, &
350 31.5105, 44.2001, 62.0000, 85.7750, 109.5500, 133.3250, &
351 157.1000, 180.8750, 204.6500, 228.4250, 252.2000, 275.9750, &
352 299.7500, 323.5250, 347.3000, 371.0750, 394.8500, 418.6250, &
353 442.4000, 466.1750, 489.9500, 513.7250, 537.5000, 561.2750, &
354 585.0500, 608.8250, 632.6000, 656.3750, 680.1500, 703.9250, &
355 727.7000, 751.4750, 775.2500, 799.0250, 822.8000, 846.5750, &
356 870.3500, 894.1250, 917.9000, 941.6750, 965.4500, 989.2250, &
357 1013.0000/
358 !
359 data (ozone(i,5),i=1,np)/ &
360 0.1993E-06, 0.2521E-06, 0.3051E-06, 0.3585E-06, 0.4121E-06, &
361 0.4661E-06, 0.5203E-06, 0.5748E-06, 0.6296E-06, 0.6847E-06, &
362 0.7402E-06, 0.7959E-06, 0.8519E-06, 0.9096E-06, 0.1125E-05, &
363 0.1450E-05, 0.1794E-05, 0.2156E-05, 0.2538E-05, 0.2939E-05, &
364 0.3362E-05, 0.3785E-05, 0.4753E-05, 0.6005E-05, 0.7804E-05, &
365 0.9635E-05, 0.1023E-04, 0.1067E-04, 0.1177E-04, 0.1290E-04, &
366 0.1134E-04, 0.9223E-05, 0.6667E-05, 0.3644E-05, 0.1545E-05, &
367 0.5355E-06, 0.2523E-06, 0.2062E-06, 0.1734E-06, 0.1548E-06, &
368 0.1360E-06, 0.1204E-06, 0.1074E-06, 0.9707E-07, 0.8960E-07, &
369 0.8419E-07, 0.7962E-07, 0.7542E-07, 0.7290E-07, 0.7109E-07, &
370 0.6940E-07, 0.6786E-07, 0.6635E-07, 0.6500E-07, 0.6370E-07, &
371 0.6244E-07, 0.6132E-07, 0.6022E-07, 0.5914E-07, 0.5884E-07, &
372 0.5855E-07, 0.5823E-07, 0.5772E-07, 0.5703E-07, 0.5635E-07, &
373 0.5570E-07, 0.5492E-07, 0.5412E-07, 0.5335E-07, 0.5260E-07, &
374 0.5167E-07, 0.5063E-07, 0.4961E-07, 0.4860E-07, 0.4761E-07/
375
376 !--------------------------------------------------------------------------------
377
378 #ifdef WRF_CHEM
379 IF ( aer_ra_feedback == 1) then
380 IF ( .NOT. &
381 ( PRESENT(tauaer300) .AND. &
382 PRESENT(tauaer400) .AND. &
383 PRESENT(tauaer600) .AND. &
384 PRESENT(tauaer999) .AND. &
385 PRESENT(gaer300) .AND. &
386 PRESENT(gaer400) .AND. &
387 PRESENT(gaer600) .AND. &
388 PRESENT(gaer999) .AND. &
389 PRESENT(waer300) .AND. &
390 PRESENT(waer400) .AND. &
391 PRESENT(waer600) .AND. &
392 PRESENT(waer999) ) ) THEN
393 CALL wrf_error_fatal ( 'Warning: missing fields required for aerosol radiation' )
394 ENDIF
395 ENDIF
396 #endif
397 cldwater = .true.
398 overcast = .false.
399
400 mix=ite-its+1
401 mkx=kte-kts+1
402
403 is_summer=80
404 ie_summer=265
405
406 ! testing, need to change iprof, which is function of lat and julian day
407 ! iprof = 1 : mid-latitude summer profile
408 ! = 2 : mid-latitude winter profile
409 ! = 3 : sub-arctic summer profile
410 ! = 4 : sub-arctic winter profile
411 ! = 5 : tropical profile
412
413 IF (abs(center_lat) .le. 30. ) THEN ! tropic
414 iprof = 5
415 ELSE
416 IF (center_lat .gt. 0.) THEN
417 IF (center_lat .gt. 60. ) THEN ! arctic
418 IF (JULDAY .gt. is_summer .and. JULDAY .lt. ie_summer ) THEN
419 ! arctic summer
420 iprof = 3
421 ELSE
422 ! arctic winter
423 iprof = 4
424 ENDIF
425 ELSE ! midlatitude
426 IF (JULDAY .gt. is_summer .and. JULDAY .lt. ie_summer ) THEN
427 ! north midlatitude summer
428 iprof = 1
429 ELSE
430 ! north midlatitude winter
431 iprof = 2
432 ENDIF
433 ENDIF
434
435 ELSE
436 IF (center_lat .lt. -60. ) THEN ! antarctic
437 IF (JULDAY .lt. is_summer .or. JULDAY .gt. ie_summer ) THEN
438 ! antarctic summer
439 iprof = 3
440 ELSE
441 ! antarctic winter
442 iprof = 4
443 ENDIF
444 ELSE ! midlatitude
445 IF (JULDAY .lt. is_summer .or. JULDAY .gt. ie_summer ) THEN
446 ! south midlatitude summer
447 iprof = 1
448 ELSE
449 ! south midlatitude winter
450 iprof = 2
451 ENDIF
452 ENDIF
453
454 ENDIF
455 ENDIF
456
457
458 j_loop: DO J=jts,jte
459
460 DO K=kts,kte
461 DO I=its,ite
462 cwc(i,k,1) = 0.
463 cwc(i,k,2) = 0.
464 ENDDO
465 ENDDO
466
467 DO K=1,np
468 p(k)=pres(k,iprof)
469 ENDDO
470
471 do k = kts,kte+1
472 do i = its,ite
473 if(k.eq.kts)then
474 phyd(i,k)=p8w3d(i,kts,j)
475 else
476 phyd(i,k)=phyd(i,k-1) - g*rho_phy(i,k-1,j)*dz8w(i,k-1,j)
477 phydmid(i,k-1)=0.5*(phyd(i,k-1)+phyd(i,k))
478 endif
479 enddo
480 enddo
481 ! normalize full pressure range
482 do k = kts+1,kte+1
483 do i = its,ite
484 if(k.eq.kts+1)fp(i) = (p8w3d(i,kts,j)-p8w3d(i,kte+1,j))/(phyd(i,kts)-phyd(i,kte+1))
485 phyd(i,k)=phyd(i,k-1) - g*rho_phy(i,k-1,j)*dz8w(i,k-1,j)*fp(i)
486 phydmid(i,k-1)=0.5*(phyd(i,k-1)+phyd(i,k))
487 enddo
488 enddo
489
490 ! reverse vars
491 !
492 DO K=kts,kte+1
493 DO I=its,ite
494 NK=kme-K+kms
495 P8W2D(I,K)=phyd(I,NK)*0.01 ! P8w2D is in mb
496 ENDDO
497 ENDDO
498
499 DO I=its,ite
500 P8W2D(I,0)=.0
501 ENDDO
502 !
503 DO K=kts,kte
504 DO I=its,ite
505 NK=kme-1-K+kms
506 TTEN2D(I,K)=0.
507 T2D(I,K)=T3D(I,NK,J)
508
509 ! SH2D specific humidity
510 SH2D(I,K)=QV3D(I,NK,J)/(1.+QV3D(I,NK,J))
511 SH2D(I,K)=max(0.,SH2D(I,K))
512 cwc(I,K,2)=QC3D(I,NK,J)
513 cwc(I,K,2)=max(0.,cwc(I,K,2))
514
515 P2D(I,K)=phydmid(I,NK)*0.01 ! P2D is in mb
516 fcld2D(I,K)=CLDFRA3D(I,NK,J)
517 ENDDO
518 ENDDO
519
520 ! This logic is tortured because cannot test F_QI unless
521 ! it is present, and order of evaluation of expressions
522 ! is not specified in Fortran
523
524 IF ( PRESENT ( F_QI ) ) THEN
525 predicate = F_QI
526 ELSE
527 predicate = .FALSE.
528 ENDIF
529
530 IF (.NOT. warm_rain .AND. .NOT. predicate ) THEN
531 DO K=kts,kte
532 DO I=its,ite
533 IF (T2D(I,K) .lt. 273.15) THEN
534 cwc(I,K,1)=cwc(I,K,2)
535 cwc(I,K,2)=0.
536 ENDIF
537 ENDDO
538 ENDDO
539 ENDIF
540
541 IF ( PRESENT( F_QNDROP ) ) THEN
542 IF ( F_QNDROP ) THEN
543 DO K=kts,kte
544 DO I=its,ite
545 NK=kme-1-K+kms
546 qndrop2d(I,K)=qndrop3d(I,NK,j)
547 ENDDO
548 ENDDO
549 qndrop2d(:,kts-1)=0.
550 END IF
551 END IF
552
553 DO I=its,ite
554 TTEN2D(I,0)=0.
555 T2D(I,0)=T2D(I,1)
556 ! SH2D specific humidity
557 SH2D(I,0)=0.5*SH2D(i,1)
558 cwc(I,0,2)=0.
559 cwc(I,0,1)=0.
560 P2D(I,0)=0.5*(P8W2D(I,0)+P8W2D(I,1))
561 fcld2D(I,0)=0.
562 ENDDO
563 !
564 IF ( PRESENT( F_QI ) .AND. PRESENT( qi3d) ) THEN
565 IF ( (F_QI) ) THEN
566 DO K=kts,kte
567 DO I=its,ite
568 NK=kme-1-K+kms
569 cwc(I,K,1)=QI3D(I,NK,J)
570 cwc(I,K,1)=max(0.,cwc(I,K,1))
571 ENDDO
572 ENDDO
573 ENDIF
574 ENDIF
575 !
576 ! ... Vertical profiles for ozone
577 !
578 call o3prof (np, p, ozone(1,iprof), its, ite, kts-1, kte, P2D, O3)
579
580 ! ... Vertical profiles for effective particle size
581 !
582 pi = 4.*atan(1.0)
583 third=1./3.
584 rhoh2o=1.e3
585 relconst=3/(4.*pi*rhoh2o)
586 ! minimun liquid water path to calculate rel
587 ! corresponds to optical depth of 1.e-3 for radius 4 microns.
588 lwpmin=3.e-5
589 do k = kts-1, kte
590 do i = its, ite
591 reff(i,k,2) = 10.
592 if( PRESENT( F_QNDROP ) ) then
593 if( F_QNDROP ) then
594 if ( cwc(i,k,2)*(P8W2D(I,K+1)-P8W2D(I,K)).gt.lwpmin.and. &
595 qndrop2d(i,k).gt.1000. ) then
596 reff(i,k,2)=(relconst*cwc(i,k,2)/qndrop2d(i,k))**third ! effective radius in m
597 ! apply scaling from Martin et al., JAS 51, 1830.
598 reff(i,k,2)=1.1*reff(i,k,2)
599 reff(i,k,2)=reff(i,k,2)*1.e6 ! convert from m to microns
600 reff(i,k,2)=max(reff(i,k,2),4.)
601 reff(i,k,2)=min(reff(i,k,2),20.)
602 end if
603 end if
604 end if
605 reff(i,k,1) = 80.
606 end do
607 end do
608 !
609 ! ... Level indices separating high, middle and low clouds
610 !
611 do i = its, ite
612 p400(i) = 1.e5
613 p700(i) = 1.e5
614 enddo
615
616 do k = kts-1,kte+1
617 do i = its, ite
618 if (abs(P8W2D(i,k) - 400.) .lt. p400(i)) then
619 p400(i) = abs(P8W2D(i,k) - 400.)
620 ict(i) = k
621 endif
622 if (abs(P8W2D(i,k) - 700.) .lt. p700(i)) then
623 p700(i) = abs(P8W2D(i,k) - 700.)
624 icb(i) = k
625 endif
626 end do
627 end do
628
629 !wig beg
630 ! ... Aerosol effects. Added aerosol feedbacks with MOSAIC, Dec. 2005.
631 !
632 do ib = 1, 11
633 do k = kts-1,kte
634 do i = its,ite
635 taual(i,k,ib) = 0.
636 ssaal(i,k,ib) = 0.
637 asyal(i,k,ib) = 0.
638 end do
639 end do
640 end do
641
642 #ifdef WRF_CHEM
643 IF ( AER_RA_FEEDBACK == 1) then
644 !wig end
645 do ib = 1, 11
646 do k = kts-1,kte-1 !wig
647 do i = its,ite
648
649 ! taual(i,kte-k,ib) = 0.
650 ! ssaal(i,kte-k,ib) = 0.
651 ! asyal(i,kte-k,ib) = 0.
652
653 !jcb beg
654 ! convert optical properties at 300,400,600, and 999 to conform to the band wavelengths
655 ! these are: 200,235,270,287.5,302.5,305,362.5,550,1920,1745,6135; why the emphasis on the UV?
656 ! taual - use angstrom exponent
657 if(tauaer300(i,k+1,j).gt.thresh .and. tauaer999(i,k+1,j).gt.thresh) then
658 ang=log(tauaer300(i,k+1,j)/tauaer999(i,k+1,j))/log(999./300.)
659 ! write(6,*)i,k,ang,tauaer300(i,k+1,j),tauaer999(i,k+1,j)
660 taual(i,kte-k,ib)=tauaer400(i,k+1,j)*(0.4/midbands(ib))**ang ! notice reserved variable
661 ! write(6,10001)i,k,ang,tauaer300(i,k+1,j),tauaer999(i,k+1,j),midbands(ib),taual(i,k,ib)
662 !10001 format(i3,i3,5f12.6)
663
664 ! ssa - linear interpolation; extrapolation
665 slope=(waer600(i,k+1,j)-waer400(i,k+1,j))/.2
666 ssaal(i,kte-k,ib) = slope*(midbands(ib)-.6)+waer600(i,k+1,j) ! notice reversed variables
667 if(ssaal(i,kte-k,ib).lt.0.4) ssaal(i,kte-k,ib)=0.4
668 if(ssaal(i,kte-k,ib).ge.1.0) ssaal(i,kte-k,ib)=1.0
669
670 ! g - linear interpolation;extrapolation
671 slope=(gaer600(i,k+1,j)-gaer400(i,k+1,j))/.2
672 asyal(i,kte-k,ib) = slope*(midbands(ib)-.6)+gaer600(i,k+1,j) ! notice reversed varaibles
673 if(asyal(i,kte-k,ib).lt.0.5) asyal(i,kte-k,ib)=0.5
674 if(asyal(i,kte-k,ib).ge.1.0) asyal(i,kte-k,ib)=1.0
675 endif
676 !jcb end
677 end do
678 end do
679 end do
680
681 !wig beg
682 do ib = 1, 11
683 do i = its,ite
684 slope = 0. !use slope as a sum holder
685 do k = kts-1,kte
686 slope = slope + taual(i,k,ib)
687 end do
688 if( slope < 0. ) then
689 write(msg,'("ERROR: Negative total optical depth of ",f8.2," at point i,j,ib=",3i5)') slope,i,j,ib
690 call wrf_error_fatal(msg)
691 else if( slope > 5. ) then
692 call wrf_message("-------------------------")
693 write(msg,'("WARNING: Large total optical depth of ",f8.2," at point i,j,ib=",3i5)') slope,i,j,ib
694 call wrf_message(msg)
695
696 call wrf_message("Diagnostics 1: k, tauaer300, tauaer400, tauaer600, tauaer999")
697 do k=kts,kte
698 write(msg,'(i4,4f8.2)') k, tauaer300(i,k,j), tauaer400(i,k,j), &
699 tauaer600(i,k,j), tauaer999(i,k,j)
700 call wrf_message(msg)
701 end do
702
703 call wrf_message("Diagnostics 2: k, gaer300, gaer400, gaer600, gaer999")
704 do k=kts,kte
705 write(msg,'(i4,4f8.2)') k, gaer300(i,k,j), gaer400(i,k,j), &
706 gaer600(i,k,j), gaer999(i,k,j)
707 call wrf_message(msg)
708 end do
709
710 call wrf_message("Diagnostics 3: k, waer300, waer400, waer600, waer999")
711 do k=kts,kte
712 write(msg,'(i4,4f8.2)') k, waer300(i,k,j), waer400(i,k,j), &
713 waer600(i,k,j), waer999(i,k,j)
714 call wrf_message(msg)
715 end do
716
717 call wrf_message("Diagnostics 4: k, ssaal, asyal, taual")
718 do k=kts-1,kte
719 write(msg,'(i4,3f8.2)') k, ssaal(i,k,ib), asyal(i,k,ib), taual(i,k,ib)
720 call wrf_message(msg)
721 end do
722 call wrf_message("-------------------------")
723 end if
724 end do
725 end do
726 !wig end
727 endif
728 #endif
729 !
730 ! ... Initialize output arrays
731 !
732 do ib = 1, 2
733 do k = kts-1, kte
734 do i = its, ite
735 taucld(i,k,ib) = 0.
736 end do
737 end do
738 end do
739 !
740 do k = kts-1,kte+1
741 do i = its,ite
742 flx(i,k) = 0.
743 flxd(i,k) = 0.
744 end do
745 end do
746 !
747 ! ... Solar zenith angle
748 !
749 do i = its,ite
750 xt24 = mod(xtime + radfrq * 0.5, 1440.)
751 tloctm = GMT + xt24 / 60. + XLONG(i,j) / 15.
752 hrang = 15. * (tloctm - 12.) * degrad
753 xxlat = XLAT(i,j) * degrad
754 cosz(i) = sin(xxlat) * sin(declin) + &
755 cos(xxlat) * cos(declin) * cos(hrang)
756 !urban
757 if(present(COSZ_URB2D)) COSZ_URB2D(i,j)=cosz(i) !urban
758 if(present(OMG_URB2D)) OMG_URB2D(i,j)=hrang !urban
759 rsuvbm(i) = ALB(i,j)
760 rsuvdf(i) = ALB(i,j)
761 rsirbm(i) = ALB(i,j)
762 rsirdf(i) = ALB(i,j)
763 end do
764
765 call sorad (mix,1,1,mkx+1,p8w2D,t2D,sh2D,o3, &
766 overcast,cldwater,cwc,taucld,reff,fcld2D,ict,icb,&
767 taual,ssaal,asyal, &
768 cosz,rsuvbm,rsuvdf,rsirbm,rsirdf, &
769 flx,flxd)
770 !
771 ! ... Convert the units of flx and flc from fraction to w/m^2
772 !
773 do k = kts, kte
774 do i = its, ite
775 nk=kme-1-k+kms
776 if(present(taucldc)) taucldc(i,nk,j)=taucld(i,k,2)
777 if(present(taucldi)) taucldi(i,nk,j)=taucld(i,k,1)
778 enddo
779 enddo
780
781 do k = kts, kte+1
782 do i = its, ite
783 if (cosz(i) .lt. thresh) then
784 flx(i,k) = 0.
785 else
786 flx(i,k) = flx(i,k) * SOLCON * cosz(i)
787 endif
788 end do
789 end do
790 !
791 ! ... Calculate heating rate (deg/sec)
792 !
793 fac = .01 * g / Cp
794 do k = kts, kte
795 do i = its, ite
796 if (cosz(i) .gt. thresh) then
797 TTEN2D(i,k) = - fac * (flx(i,k) - flx(i,k+1))/ &
798 (p8w2d(i,k)-p8w2d(i,k+1))
799 endif
800 end do
801 end do
802
803 ! upward top of atmosphere
804 do i = its, ite
805 if (cosz(i) .le. thresh) then
806 RSWTOA(i,j) = 0.
807 else
808 RSWTOA(i,j) = flx(i,kts) - flxd(i,kts) * SOLCON * cosz(i)
809 ! print *,'cosz,rswtoa=',cosz(i),rswtoa(i,j)
810 endif
811 end do
812 !
813 ! ... Absorbed part in surface energy budget
814 !
815 do i = its, ite
816 if (cosz(i) .le. thresh) then
817 GSW(i,j) = 0.
818 else
819 GSW(i,j) = (1. - rsuvbm(i)) * flxd(i,kte+1) * SOLCON * cosz(i)
820 endif
821 end do
822
823 DO K=kts,kte
824 NK=kme-1-K+kms
825 DO I=its,ite
826 RTHRATEN(I,K,J)=RTHRATEN(I,K,J)+TTEN2D(I,NK)/pi3D(I,K,J)
827 ENDDO
828 ENDDO
829 !
830 ENDDO j_loop
831
832 END SUBROUTINE GSFCSWRAD
833
834 !********************* Version Solar-6 (May 8, 1997) *****************
835
836 subroutine sorad (m,n,ndim,np,pl,ta,wa,oa, &
837 overcast,cldwater,cwc,taucld,reff,fcld,ict,icb, &
838 taual,ssaal,asyal, &
839 cosz,rsuvbm,rsuvdf,rsirbm,rsirdf, &
840 flx,flxd)
841
842 !************************************************************************
843 !
844 ! Version Solar-6 (May 8, 1997)
845 !
846 ! New feature of this version is:
847 ! (1) An option is added for scaling the cloud optical thickness. If
848 ! the fractional cloud cover, fcld, in an atmospheric model is alway
849 ! either 1 or 0 (i.e. partly cloudy sky is not allowed), it does
850 ! not require the scaling of cloud optical thickness, and the
851 ! option "overcast" can be set to .true. Computation is faster
852 ! with this option than with overcast=.false.
853 !
854 !**********************************************************************
855 !
856 ! Version Solar-5 (April 1997)
857 !
858 ! New features of this version are:
859 ! (1) Cloud optical properties can be computed from cloud water/ice
860 ! amount and the effective particle size.
861 ! (2) Aerosol optical properties are functions of height and band.
862 ! (3) A maximum-random cloud overlapping approximation is applied.
863 !
864 !*********************************************************************
865 !
866 ! This routine computes solar fluxes due to the absoption by water
867 ! vapor, ozone, co2, o2, clouds, and aerosols and due to the
868 ! scattering by clouds, aerosols, and gases.
869 !
870 ! The solar spectrum is divided into one UV+visible band and three IR
871 ! bands separated by the wavelength 0.7 micron. The UV+visible band
872 ! is further divided into eight sub-bands.
873 !
874 ! This is a vectorized code. It computes fluxes simultaneously for
875 ! (m x n) soundings, which is a subset of (m x ndim) soundings.
876 ! In a global climate model, m and ndim correspond to the numbers of
877 ! grid boxes in the zonal and meridional directions, respectively.
878 !
879 ! Ice and liquid cloud particles are allowed to co-exist in a layer.
880 !
881 ! There is an option of providing either cloud ice/water mixing ratio
882 ! (cwc) or thickness (taucld). If the former is provided, set
883 ! cldwater=.true., and taucld will be computed from cwc and reff as a
884 ! function of spectra band. Otherwise, set cldwater=.false., and
885 ! specify taucld, independent of spectral band.
886 !
887 ! If no information is available for reff, a default value of
888 ! 10 micron for liquid water and 75 micron for ice can be used.
889 ! For a clear layer, reff can be set to any values except zero.
890 !
891 ! The maximum-random assumption is applied for treating cloud
892 ! overlapping.
893
894 ! Clouds are grouped into high, middle, and low clouds separated by
895 ! the level indices ict and icb. For detail, see subroutine cldscale.
896 !
897 ! In a high spatial-resolution atmospheric model, fractional cloud cover
898 ! might be computed to be either 0 or 1. In such a case, scaling of the
899 ! cloud optical thickness is not necessary, and the computation can be
900 ! made faster by setting overcast=.true. The option overcast=.false.
901 ! can be applied to any values of the fractional cloud cover, but the
902 ! computation is slower.
903 !
904 ! Aerosol optical thickness, single-scattering albaedo, and asymmtry
905 ! factor can be specified as functions of height and spectral band.
906 !
907 !----- Input parameters:
908 ! units size
909 ! number of soundings in zonal direction (m) n/d 1
910 ! number of soundings in meridional direction (n) n/d 1
911 ! maximum number of soundings in n/d 1
912 ! meridional direction (ndim>=n)
913 ! number of atmospheric layers (np) n/d 1
914 ! level pressure (pl) mb m*ndim*(np+1)
915 ! layer temperature (ta) k m*ndim*np
916 ! layer specific humidity (wa) gm/gm m*ndim*np
917 ! layer ozone concentration (oa) gm/gm m*ndim*np
918 ! co2 mixing ratio by volumn (co2) pppv 1
919 ! option for scaling cloud optical thickness n/d 1
920 ! overcast="true" if scaling is NOT required
921 ! overcast="fasle" if scaling is required
922 ! option for cloud optical thickness n/d 1
923 ! cldwater="true" if cwc is provided
924 ! cldwater="false" if taucld is provided
925 ! cloud water mixing ratio (cwc) gm/gm m*ndim*np*2
926 ! index 1 for ice particles
927 ! index 2 for liquid drops
928 ! cloud optical thickness (taucld) n/d m*ndim*np*2
929 ! index 1 for ice particles
930 ! index 2 for liquid drops
931 ! effective cloud-particle size (reff) micrometer m*ndim*np*2
932 ! index 1 for ice particles
933 ! index 2 for liquid drops
934 ! cloud amount (fcld) fraction m*ndim*np
935 ! level index separating high and middle n/d 1
936 ! clouds (ict)
937 ! level index separating middle and low n/d 1
938 ! clouds (icb)
939 ! aerosol optical thickness (taual) n/d m*ndim*np*11
940 ! aerosol single-scattering albedo (ssaal) n/d m*ndim*np*11
941 ! aerosol asymmetry factor (asyal) n/d m*ndim*np*11
942 ! in the uv region :
943 ! index 1 for the 0.175-0.225 micron band
944 ! index 2 for the 0.225-0.245; 0.260-0.280 micron band
945 ! index 3 for the 0.245-0.260 micron band
946 ! index 4 for the 0.280-0.295 micron band
947 ! index 5 for the 0.295-0.310 micron band
948 ! index 6 for the 0.310-0.320 micron band
949 ! index 7 for the 0.325-0.400 micron band
950 ! in the par region :
951 ! index 8 for the 0.400-0.700 micron band
952 ! in the infrared region :
953 ! index 9 for the 0.700-1.220 micron band
954 ! index 10 for the 1.220-2.270 micron band
955 ! index 11 for the 2.270-10.00 micron band
956 ! cosine of solar zenith angle (cosz) n/d m*ndim
957 ! uv+visible sfc albedo for beam radiation
958 ! for wavelengths<0.7 micron (rsuvbm) fraction m*ndim
959 ! uv+visible sfc albedo for diffuse radiation
960 ! for wavelengths<0.7 micron (rsuvdf) fraction m*ndim
961 ! ir sfc albedo for beam radiation
962 ! for wavelengths>0.7 micron (rsirbm) fraction m*ndim
963 ! ir sfc albedo for diffuse radiation (rsirdf) fraction m*ndim
964 !
965 !----- Output parameters
966 !
967 ! all-sky flux (downward minus upward) (flx) fraction m*ndim*(np+1)
968 ! clear-sky flux (downward minus upward) (flc) fraction m*ndim*(np+1)
969 ! all-sky direct downward uv (0.175-0.4 micron)
970 ! flux at the surface (fdiruv) fraction m*ndim
971 ! all-sky diffuse downward uv flux at
972 ! the surface (fdifuv) fraction m*ndim
973 ! all-sky direct downward par (0.4-0.7 micron)
974 ! flux at the surface (fdirpar) fraction m*ndim
975 ! all-sky diffuse downward par flux at
976 ! the surface (fdifpar) fraction m*ndim
977 ! all-sky direct downward ir (0.7-10 micron)
978 ! flux at the surface (fdirir) fraction m*ndim
979 ! all-sky diffuse downward ir flux at
980 ! the surface (fdifir) fraction m*ndim
981 !
982 !----- Notes:
983 !
984 ! (1) The unit of "flux" is fraction of the incoming solar radiation
985 ! at the top of the atmosphere. Therefore, fluxes should
986 ! be equal to "flux" multiplied by the extra-terrestrial solar
987 ! flux and the cosine of solar zenith angle.
988 ! (2) pl(i,j,1) is the pressure at the top of the model, and
989 ! pl(i,j,np+1) is the surface pressure.
990 ! (3) the pressure levels ict and icb correspond approximately
991 ! to 400 and 700 mb.
992 ! (4) if overcast='true', the clear-sky flux, flc, is not computed.
993 !
994 !**************************************************************************
995 implicit none
996 !**************************************************************************
997
998 !-----input parameters
999
1000 integer m,n,ndim,np
1001 integer ict(m,ndim),icb(m,ndim)
1002 real pl(m,ndim,np+1),ta(m,ndim,np),wa(m,ndim,np),oa(m,ndim,np)
1003 real cwc(m,ndim,np,2),taucld(m,ndim,np,2),reff(m,ndim,np,2), &
1004 fcld(m,ndim,np)
1005 real taual(m,ndim,np,11),ssaal(m,ndim,np,11),asyal(m,ndim,np,11)
1006 real cosz(m,ndim),rsuvbm(m,ndim),rsuvdf(m,ndim), &
1007 rsirbm(m,ndim),rsirdf(m,ndim)
1008 logical overcast,cldwater
1009
1010 !-----output parameters
1011
1012 real flx(m,ndim,np+1),flc(m,ndim,np+1)
1013 real flxu(m,ndim,np+1),flxd(m,ndim,np+1)
1014 real fdiruv (m,ndim),fdifuv (m,ndim)
1015 real fdirpar(m,ndim),fdifpar(m,ndim)
1016 real fdirir (m,ndim),fdifir (m,ndim)
1017
1018 !-----temporary array
1019
1020 integer i,j,k
1021 real cwp(m,n,np,2)
1022 real dp(m,n,np),wh(m,n,np),oh(m,n,np),scal(m,n,np)
1023 real swh(m,n,np+1),so2(m,n,np+1),df(m,n,np+1)
1024 real sdf(m,n),sclr(m,n),csm(m,n),x
1025
1026 do j= 1, n
1027 do i= 1, m
1028 if (pl(i,j,1) .eq. 0.0) then
1029 pl(i,j,1)=1.0e-4
1030 endif
1031 enddo
1032 enddo
1033
1034 do j= 1, n
1035 do i= 1, m
1036
1037 swh(i,j,1)=0.
1038 so2(i,j,1)=0.
1039
1040 !-----csm is the effective secant of the solar zenith angle
1041 ! see equation (12) of Lacis and Hansen (1974, JAS)
1042
1043 csm(i,j)=35./sqrt(1224.*cosz(i,j)*cosz(i,j)+1.)
1044
1045 enddo
1046 enddo
1047
1048 do k= 1, np
1049 do j= 1, n
1050 do i= 1, m
1051
1052 !-----compute layer thickness and pressure-scaling function.
1053 ! indices for the surface level and surface layer
1054 ! are np+1 and np, respectively.
1055
1056 dp(i,j,k)=pl(i,j,k+1)-pl(i,j,k)
1057 scal(i,j,k)=dp(i,j,k)*(.5*(pl(i,j,k)+pl(i,j,k+1))/300.)**.8
1058
1059 !-----compute scaled water vapor amount, unit is g/cm**2
1060 ! note: the sign prior to the constant 0.00135 was incorrectly
1061 ! set to negative in the previous version
1062
1063 wh(i,j,k)=1.02*wa(i,j,k)*scal(i,j,k)* &
1064 (1.+0.00135*(ta(i,j,k)-240.)) +1.e-11
1065 swh(i,j,k+1)=swh(i,j,k)+wh(i,j,k)
1066
1067 !-----compute ozone amount, unit is (cm-atm)stp
1068 ! the number 466.7 is a conversion factor from g/cm**2 to (cm-atm)stp
1069
1070 oh(i,j,k)=1.02*oa(i,j,k)*dp(i,j,k)*466.7 +1.e-11
1071
1072 !-----compute layer cloud water amount (gm/m**2)
1073 ! the index is 1 for ice crystals and 2 for liquid drops
1074
1075 cwp(i,j,k,1)=1.02*10000.*cwc(i,j,k,1)*dp(i,j,k)
1076 cwp(i,j,k,2)=1.02*10000.*cwc(i,j,k,2)*dp(i,j,k)
1077
1078 enddo
1079 enddo
1080 enddo
1081
1082 !-----initialize fluxes for all-sky (flx), clear-sky (flc), and
1083 ! flux reduction (df)
1084
1085 do k=1, np+1
1086 do j=1, n
1087 do i=1, m
1088 flx(i,j,k)=0.
1089 flc(i,j,k)=0.
1090 flxu(i,j,k)=0.
1091 flxd(i,j,k)=0.
1092 df(i,j,k)=0.
1093 enddo
1094 enddo
1095 enddo
1096
1097 !-----compute solar uv and par fluxes
1098
1099 call soluv (m,n,ndim,np,oh,dp,overcast,cldwater, &
1100 cwp,taucld,reff,ict,icb,fcld,cosz, &
1101 taual,ssaal,asyal,csm,rsuvbm,rsuvdf, &
1102 flx,flc,flxu,flxd,fdiruv,fdifuv,fdirpar,fdifpar)
1103
1104 !-----compute and update solar ir fluxes
1105
1106 call solir (m,n,ndim,np,wh,overcast,cldwater, &
1107 cwp,taucld,reff,ict,icb,fcld,cosz, &
1108 taual,ssaal,asyal,csm,rsirbm,rsirdf, &
1109 flx,flc,flxu,flxd,fdirir,fdifir)
1110
1111 !-----compute scaled o2 amount, unit is (cm-atm)stp.
1112
1113 do k= 1, np
1114 do j= 1, n
1115 do i= 1, m
1116 so2(i,j,k+1)=so2(i,j,k)+165.22*scal(i,j,k)
1117 enddo
1118 enddo
1119 enddo
1120
1121 !-----compute flux reduction due to oxygen following
1122 ! chou (J. climate, 1990). The fraction 0.0287 is the
1123 ! extraterrestrial solar flux in the o2 bands.
1124
1125 do k= 2, np+1
1126 do j= 1, n
1127 do i= 1, m
1128 x=so2(i,j,k)*csm(i,j)
1129 df(i,j,k)=df(i,j,k)+0.0287*(1.-exp(-0.00027*sqrt(x)))
1130 enddo
1131 enddo
1132 enddo
1133
1134 !-----compute scaled co2 amounts. unit is (cm-atm)stp.
1135
1136 do k= 1, np
1137 do j= 1, n
1138 do i= 1, m
1139 so2(i,j,k+1)=so2(i,j,k)+co2*789.*scal(i,j,k)+1.e-11
1140 enddo
1141 enddo
1142 enddo
1143
1144 !-----compute and update flux reduction due to co2 following
1145 ! chou (J. Climate, 1990)
1146
1147 call flxco2(m,n,np,so2,swh,csm,df)
1148
1149 !-----adjust for the effect of o2 cnd co2 on clear-sky fluxes.
1150
1151 do k= 2, np+1
1152 do j= 1, n
1153 do i= 1, m
1154 flc(i,j,k)=flc(i,j,k)-df(i,j,k)
1155 enddo
1156 enddo
1157 enddo
1158
1159 !-----adjust for the all-sky fluxes due to o2 and co2. It is
1160 ! assumed that o2 and co2 have no effects on solar radiation
1161 ! below clouds.
1162
1163 do j=1,n
1164 do i=1,m
1165 sdf(i,j)=0.0
1166 sclr(i,j)=1.0
1167 enddo
1168 enddo
1169
1170 do k=1,np
1171 do j=1,n
1172 do i=1,m
1173
1174 !-----sclr is the fraction of clear sky.
1175 ! sdf is the flux reduction below clouds.
1176
1177 if(fcld(i,j,k).gt.0.01) then
1178 sdf(i,j)=sdf(i,j)+df(i,j,k)*sclr(i,j)*fcld(i,j,k)
1179 sclr(i,j)=sclr(i,j)*(1.-fcld(i,j,k))
1180 endif
1181 flx(i,j,k+1)=flx(i,j,k+1)-sdf(i,j)-df(i,j,k+1)*sclr(i,j)
1182 flxu(i,j,k+1)=flxu(i,j,k+1)-sdf(i,j)-df(i,j,k+1)*sclr(i,j)
1183 flxd(i,j,k+1)=flxd(i,j,k+1)-sdf(i,j)-df(i,j,k+1)*sclr(i,j) ! SG: same as flux????
1184
1185 enddo
1186 enddo
1187 enddo
1188
1189 !-----adjustment for the direct downward ir flux.
1190
1191 do j= 1, n
1192 do i= 1, m
1193 flc(i,j,np+1)=flc(i,j,np+1)+df(i,j,np+1)*rsirbm(i,j)
1194 flx(i,j,np+1)=flx(i,j,np+1)+(sdf(i,j)+ &
1195 df(i,j,np+1)*sclr(i,j))*rsirbm(i,j)
1196 flxu(i,j,np+1)=flxu(i,j,np+1)+(sdf(i,j)+ &
1197 df(i,j,np+1)*sclr(i,j))*rsirbm(i,j)
1198 flxd(i,j,np+1)=flxd(i,j,np+1)+(sdf(i,j)+ &
1199 df(i,j,np+1)*sclr(i,j))*rsirbm(i,j)
1200 fdirir(i,j)=fdirir(i,j)-(sdf(i,j)+df(i,j,np+1)*sclr(i,j))
1201 enddo
1202 enddo
1203
1204 end subroutine sorad
1205
1206 !************************************************************************
1207
1208 subroutine soluv (m,n,ndim,np,oh,dp,overcast,cldwater, &
1209 cwp,taucld,reff,ict,icb,fcld,cosz, &
1210 taual,ssaal,asyal,csm,rsuvbm,rsuvdf, &
1211 flx,flc,flxu,flxd,fdiruv,fdifuv,fdirpar,fdifpar)
1212
1213 !************************************************************************
1214 ! compute solar fluxes in the uv+par region. the spectrum is
1215 ! grouped into 8 bands:
1216 !
1217 ! Band Micrometer
1218 !
1219 ! UV-C 1. .175 - .225
1220 ! 2. .225 - .245
1221 ! .260 - .280
1222 ! 3. .245 - .260
1223 !
1224 ! UV-B 4. .280 - .295
1225 ! 5. .295 - .310
1226 ! 6. .310 - .320
1227 !
1228 ! UV-A 7. .320 - .400
1229 !
1230 ! PAR 8. .400 - .700
1231 !
1232 !----- Input parameters: units size
1233 !
1234 ! number of soundings in zonal direction (m) n/d 1
1235 ! number of soundings in meridional direction (n) n/d 1
1236 ! maximum number of soundings in n/d 1
1237 ! meridional direction (ndim)
1238 ! number of atmospheric layers (np) n/d 1
1239 ! layer ozone content (oh) (cm-atm)stp m*n*np
1240 ! layer pressure thickness (dp) mb m*n*np
1241 ! option for scaling cloud optical thickness n/d 1
1242 ! overcast="true" if scaling is NOT required
1243 ! overcast="fasle" if scaling is required
1244 ! input option for cloud optical thickness n/d 1
1245 ! cldwater="true" if taucld is provided
1246 ! cldwater="false" if cwp is provided
1247 ! cloud water amount (cwp) gm/m**2 m*n*np*2
1248 ! index 1 for ice particles
1249 ! index 2 for liquid drops
1250 ! cloud optical thickness (taucld) n/d m*ndim*np*2
1251 ! index 1 for ice paticles
1252 ! index 2 for liquid particles
1253 ! effective cloud-particle size (reff) micrometer m*ndim*np*2
1254 ! index 1 for ice paticles
1255 ! index 2 for liquid particles
1256 ! level indiex separating high and n/d m*n
1257 ! middle clouds (ict)
1258 ! level indiex separating middle and n/d m*n
1259 ! low clouds (icb)
1260 ! cloud amount (fcld) fraction m*ndim*np
1261 ! cosine of solar zenith angle (cosz) n/d m*ndim
1262 ! aerosol optical thickness (taual) n/d m*ndim*np*11
1263 ! aerosol single-scattering albedo (ssaal) n/d m*ndim*np*11
1264 ! aerosol asymmetry factor (asyal) n/d m*ndim*np*11
1265 ! cosecant of the solar zenith angle (csm) n/d m*n
1266 ! uv+par surface albedo for beam fraction m*ndim
1267 ! radiation (rsuvbm)
1268 ! uv+par surface albedo for diffuse fraction m*ndim
1269 ! radiation (rsuvdf)
1270 !
1271 !---- temporary array
1272 !
1273 ! scaled cloud optical thickness n/d m*n*np
1274 ! for beam radiation (tauclb)
1275 ! scaled cloud optical thickness n/d m*n*np
1276 ! for diffuse radiation (tauclf)
1277 !
1278 !----- output (updated) parameters:
1279 !
1280 ! all-sky net downward flux (flx) fraction m*ndim*(np+1)
1281 ! clear-sky net downward flux (flc) fraction m*ndim*(np+1)
1282 ! all-sky direct downward uv flux at
1283 ! the surface (fdiruv) fraction m*ndim
1284 ! all-sky diffuse downward uv flux at
1285 ! the surface (fdifuv) fraction m*ndim
1286 ! all-sky direct downward par flux at
1287 ! the surface (fdirpar) fraction m*ndim
1288 ! all-sky diffuse downward par flux at
1289 ! the surface (fdifpar) fraction m*ndim
1290 !
1291 !***********************************************************************
1292 implicit none
1293 !***********************************************************************
1294
1295 !-----input parameters
1296
1297 integer m,n,ndim,np
1298 integer ict(m,ndim),icb(m,ndim)
1299 real taucld(m,ndim,np,2),reff(m,ndim,np,2),fcld(m,ndim,np)
1300 real cc(m,n,3),cosz(m,ndim)
1301 real cwp(m,n,np,2),oh(m,n,np),dp(m,n,np)
1302 real taual(m,ndim,np,11),ssaal(m,ndim,np,11),asyal(m,ndim,np,11)
1303 real rsuvbm(m,ndim),rsuvdf(m,ndim),csm(m,n)
1304 logical overcast,cldwater
1305
1306 !-----output (updated) parameter
1307
1308 real flx(m,ndim,np+1),flc(m,ndim,np+1)
1309 real flxu(m,ndim,np+1),flxd(m,ndim,np+1)
1310 real fdiruv (m,ndim),fdifuv (m,ndim)
1311 real fdirpar(m,ndim),fdifpar(m,ndim)
1312
1313 !-----static parameters
1314
1315 integer nband
1316 parameter (nband=8)
1317 real hk(nband),xk(nband),ry(nband)
1318 real aig(3),awg(3)
1319
1320 !-----temporary array
1321
1322 integer i,j,k,ib
1323 real tauclb(m,n,np),tauclf(m,n,np),asycl(m,n,np)
1324 real taurs,tauoz,tausto,ssatau,asysto,tauto,ssato,asyto
1325 real taux,reff1,reff2,g1,g2
1326 real td(m,n,np+1,2),rr(m,n,np+1,2),tt(m,n,np+1,2), &
1327 rs(m,n,np+1,2),ts(m,n,np+1,2)
1328 real fall(m,n,np+1),fclr(m,n,np+1),fsdir(m,n),fsdif(m,n)
1329 real fallu(m,n,np+1),falld(m,n,np+1)
1330 real asyclt(m,n)
1331 real rr1t(m,n),tt1t(m,n),td1t(m,n),rs1t(m,n),ts1t(m,n)
1332 real rr2t(m,n),tt2t(m,n),td2t(m,n),rs2t(m,n),ts2t(m,n)
1333
1334 !-----hk is the fractional extra-terrestrial solar flux in each
1335 ! of the 8 bands. the sum of hk is 0.47074.
1336
1337 data hk/.00057, .00367, .00083, .00417, &
1338 .00600, .00556, .05913, .39081/
1339
1340 !-----xk is the ozone absorption coefficient. unit: /(cm-atm)stp
1341
1342 data xk /30.47, 187.2, 301.9, 42.83, &
1343 7.09, 1.25, 0.0345, 0.0539/
1344
1345 !-----ry is the extinction coefficient for Rayleigh scattering.
1346 ! unit: /mb.
1347
1348 data ry /.00604, .00170, .00222, .00132, &
1349 .00107, .00091, .00055, .00012/
1350
1351 !-----coefficients for computing the asymmetry factor of ice clouds
1352 ! from asycl=aig(*,1)+aig(*,2)*reff+aig(*,3)*reff**2, independent
1353 ! of spectral band.
1354
1355 data aig/.74625000,.00105410,-.00000264/
1356
1357 !-----coefficients for computing the asymmetry factor of liquid
1358 ! clouds from asycl=awg(*,1)+awg(*,2)*reff+awg(*,3)*reff**2,
1359 ! independent of spectral band.
1360
1361 data awg/.82562000,.00529000,-.00014866/
1362
1363 !-----initialize fdiruv, fdifuv, surface reflectances and transmittances.
1364 ! cc is the maximum cloud cover in each of the three cloud groups.
1365
1366 do j= 1, n
1367 do i= 1, m
1368 fdiruv(i,j)=0.0
1369 fdifuv(i,j)=0.0
1370 rr(i,j,np+1,1)=rsuvbm(i,j)
1371 rr(i,j,np+1,2)=rsuvbm(i,j)
1372 rs(i,j,np+1,1)=rsuvdf(i,j)
1373 rs(i,j,np+1,2)=rsuvdf(i,j)
1374 td(i,j,np+1,1)=0.0
1375 td(i,j,np+1,2)=0.0
1376 tt(i,j,np+1,1)=0.0
1377 tt(i,j,np+1,2)=0.0
1378 ts(i,j,np+1,1)=0.0
1379 ts(i,j,np+1,2)=0.0
1380 cc(i,j,1)=0.0
1381 cc(i,j,2)=0.0
1382 cc(i,j,3)=0.0
1383 enddo
1384 enddo
1385
1386
1387 !-----compute cloud optical thickness
1388
1389 if (cldwater) then
1390
1391 do k= 1, np
1392 do j= 1, n
1393 do i= 1, m
1394 taucld(i,j,k,1)=cwp(i,j,k,1)*( 3.33e-4+2.52/reff(i,j,k,1))
1395 taucld(i,j,k,2)=cwp(i,j,k,2)*(-6.59e-3+1.65/reff(i,j,k,2))
1396 enddo
1397 enddo
1398 enddo
1399
1400 endif
1401
1402 !-----options for scaling cloud optical thickness
1403
1404 if (overcast) then
1405
1406 do k= 1, np
1407 do j= 1, n
1408 do i= 1, m
1409 tauclb(i,j,k)=taucld(i,j,k,1)+taucld(i,j,k,2)
1410 tauclf(i,j,k)=tauclb(i,j,k)
1411 enddo
1412 enddo
1413 enddo
1414
1415 do k= 1, 3
1416 do j= 1, n
1417 do i= 1, m
1418 cc(i,j,k)=1.0
1419 enddo
1420 enddo
1421 enddo
1422
1423 else
1424
1425 !-----scale cloud optical thickness in each layer from taucld (with
1426 ! cloud amount fcld) to tauclb and tauclf (with cloud amount cc).
1427 ! tauclb is the scaled optical thickness for beam radiation and
1428 ! tauclf is for diffuse radiation.
1429
1430 call cldscale(m,n,ndim,np,cosz,fcld,taucld,ict,icb, &
1431 cc,tauclb,tauclf)
1432
1433 endif
1434
1435 !-----compute cloud asymmetry factor for a mixture of
1436 ! liquid and ice particles. unit of reff is micrometers.
1437
1438 do k= 1, np
1439
1440 do j= 1, n
1441 do i= 1, m
1442
1443 asyclt(i,j)=1.0
1444
1445 taux=taucld(i,j,k,1)+taucld(i,j,k,2)
1446 if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then
1447
1448 reff1=min(reff(i,j,k,1),130.)
1449 reff2=min(reff(i,j,k,2),20.0)
1450
1451 g1=(aig(1)+(aig(2)+aig(3)*reff1)*reff1)*taucld(i,j,k,1)
1452 g2=(awg(1)+(awg(2)+awg(3)*reff2)*reff2)*taucld(i,j,k,2)
1453 asyclt(i,j)=(g1+g2)/taux
1454
1455 endif
1456
1457 enddo
1458 enddo
1459
1460 do j=1,n
1461 do i=1,m
1462 asycl(i,j,k)=asyclt(i,j)
1463 enddo
1464 enddo
1465
1466 enddo
1467
1468 !-----integration over spectral bands
1469
1470 do 100 ib=1,nband
1471
1472 do 300 k= 1, np
1473
1474 do j= 1, n
1475 do i= 1, m
1476
1477 !-----compute ozone and rayleigh optical thicknesses
1478
1479 taurs=ry(ib)*dp(i,j,k)
1480 tauoz=xk(ib)*oh(i,j,k)
1481
1482 !-----compute clear-sky optical thickness, single scattering albedo,
1483 ! and asymmetry factor
1484
1485 tausto=taurs+tauoz+taual(i,j,k,ib)+1.0e-8
1486 ssatau=ssaal(i,j,k,ib)*taual(i,j,k,ib)+taurs
1487 asysto=asyal(i,j,k,ib)*ssaal(i,j,k,ib)*taual(i,j,k,ib)
1488
1489 tauto=tausto
1490 ssato=ssatau/tauto+1.0e-8
1491 ssato=min(ssato,0.999999)
1492 asyto=asysto/(ssato*tauto)
1493
1494 !-----compute reflectance and transmittance for cloudless layers
1495
1496 !- for direct incident radiation
1497
1498 call deledd (tauto,ssato,asyto,csm(i,j), &
1499 rr1t(i,j),tt1t(i,j),td1t(i,j))
1500
1501 !- for diffuse incident radiation
1502
1503 call sagpol (tauto,ssato,asyto,rs1t(i,j),ts1t(i,j))
1504
1505 !-----compute reflectance and transmittance for cloud layers
1506
1507 if (tauclb(i,j,k).lt.0.01 .or. fcld(i,j,k).lt.0.01) then
1508
1509 rr2t(i,j)=rr1t(i,j)
1510 tt2t(i,j)=tt1t(i,j)
1511 td2t(i,j)=td1t(i,j)
1512 rs2t(i,j)=rs1t(i,j)
1513 ts2t(i,j)=ts1t(i,j)
1514
1515 else
1516
1517 !-- for direct incident radiation
1518
1519 tauto=tausto+tauclb(i,j,k)
1520 ssato=(ssatau+tauclb(i,j,k))/tauto+1.0e-8
1521 ssato=min(ssato,0.999999)
1522 asyto=(asysto+asycl(i,j,k)*tauclb(i,j,k))/(ssato*tauto)
1523
1524 call deledd (tauto,ssato,asyto,csm(i,j), &
1525 rr2t(i,j),tt2t(i,j),td2t(i,j))
1526
1527 !-- for diffuse incident radiation
1528
1529 tauto=tausto+tauclf(i,j,k)
1530 ssato=(ssatau+tauclf(i,j,k))/tauto+1.0e-8
1531 ssato=min(ssato,0.999999)
1532 asyto=(asysto+asycl(i,j,k)*tauclf(i,j,k))/(ssato*tauto)
1533
1534 call sagpol (tauto,ssato,asyto,rs2t(i,j),ts2t(i,j))
1535
1536 endif
1537
1538 enddo
1539 enddo
1540
1541 do j=1,n
1542 do i=1,m
1543 rr(i,j,k,1)=rr1t(i,j)
1544 enddo
1545 enddo
1546 do j=1,n
1547 do i=1,m
1548 tt(i,j,k,1)=tt1t(i,j)
1549 enddo
1550 enddo
1551 do j=1,n
1552 do i=1,m
1553 td(i,j,k,1)=td1t(i,j)
1554 enddo
1555 enddo
1556 do j=1,n
1557 do i=1,m
1558 rs(i,j,k,1)=rs1t(i,j)
1559 enddo
1560 enddo
1561 do j=1,n
1562 do i=1,m
1563 ts(i,j,k,1)=ts1t(i,j)
1564 enddo
1565 enddo
1566
1567 do j=1,n
1568 do i=1,m
1569 rr(i,j,k,2)=rr2t(i,j)
1570 enddo
1571 enddo
1572 do j=1,n
1573 do i=1,m
1574 tt(i,j,k,2)=tt2t(i,j)
1575 enddo
1576 enddo
1577 do j=1,n
1578 do i=1,m
1579 td(i,j,k,2)=td2t(i,j)
1580 enddo
1581 enddo
1582 do j=1,n
1583 do i=1,m
1584 rs(i,j,k,2)=rs2t(i,j)
1585 enddo
1586 enddo
1587 do j=1,n
1588 do i=1,m
1589 ts(i,j,k,2)=ts2t(i,j)
1590 enddo
1591 enddo
1592
1593 300 continue
1594
1595 !-----flux calculations
1596
1597 call cldflx (m,n,np,ict,icb,overcast,cc,rr,tt,td,rs,ts, &
1598 fclr,fall,fallu,falld,fsdir,fsdif)
1599
1600 do k= 1, np+1
1601 do j= 1, n
1602 do i= 1, m
1603 flx(i,j,k)=flx(i,j,k)+fall(i,j,k)*hk(ib)
1604 flxu(i,j,k)=flxu(i,j,k)+fallu(i,j,k)*hk(ib)
1605 flxd(i,j,k)=flxd(i,j,k)+falld(i,j,k)*hk(ib)
1606 enddo
1607 enddo
1608 do j= 1, n
1609 do i= 1, m
1610 flc(i,j,k)=flc(i,j,k)+fclr(i,j,k)*hk(ib)
1611 enddo
1612 enddo
1613 enddo
1614
1615 !-----compute downward surface fluxes in the UV and par regions
1616
1617 if(ib.lt.8) then
1618 do j=1,n
1619 do i=1,m
1620 fdiruv(i,j)=fdiruv(i,j)+fsdir(i,j)*hk(ib)
1621 fdifuv(i,j)=fdifuv(i,j)+fsdif(i,j)*hk(ib)
1622 enddo
1623 enddo
1624 else
1625 do j=1,n
1626 do i=1,m
1627 fdirpar(i,j)=fsdir(i,j)*hk(ib)
1628 fdifpar(i,j)=fsdif(i,j)*hk(ib)
1629 enddo
1630 enddo
1631 endif
1632
1633 100 continue
1634
1635 end subroutine soluv
1636
1637 !************************************************************************
1638
1639 subroutine solir (m,n,ndim,np,wh,overcast,cldwater, &
1640 cwp,taucld,reff,ict,icb,fcld,cosz, &
1641 taual,ssaal,asyal,csm,rsirbm,rsirdf, &
1642 flx,flc,flxu,flxd,fdirir,fdifir)
1643
1644 !************************************************************************
1645 ! compute solar flux in the infrared region. The spectrum is divided
1646 ! into three bands:
1647 !
1648 ! band wavenumber(/cm) wavelength (micron)
1649 ! 1( 9) 14300-8200 0.70-1.22
1650 ! 2(10) 8200-4400 1.22-2.27
1651 ! 3(11) 4400-1000 2.27-10.0
1652 !
1653 !----- Input parameters: units size
1654 !
1655 ! number of soundings in zonal direction (m) n/d 1
1656 ! number of soundings in meridional direction (n) n/d 1
1657 ! maximum number of soundings in n/d 1
1658 ! meridional direction (ndim)
1659 ! number of atmospheric layers (np) n/d 1
1660 ! layer scaled-water vapor content (wh) gm/cm^2 m*n*np
1661 ! option for scaling cloud optical thickness n/d 1
1662 ! overcast="true" if scaling is NOT required
1663 ! overcast="fasle" if scaling is required
1664 ! input option for cloud optical thickness n/d 1
1665 ! cldwater="true" if taucld is provided
1666 ! cldwater="false" if cwp is provided
1667 ! cloud water concentration (cwp) gm/m**2 m*n*np*2
1668 ! index 1 for ice particles
1669 ! index 2 for liquid drops
1670 ! cloud optical thickness (taucld) n/d m*ndim*np*2
1671 ! index 1 for ice paticles
1672 ! effective cloud-particle size (reff) micrometer m*ndim*np*2
1673 ! index 1 for ice paticles
1674 ! index 2 for liquid particles
1675 ! level index separating high and n/d m*n
1676 ! middle clouds (ict)
1677 ! level index separating middle and n/d m*n
1678 ! low clouds (icb)
1679 ! cloud amount (fcld) fraction m*ndim*np
1680 ! aerosol optical thickness (taual) n/d m*ndim*np*11
1681 ! aerosol single-scattering albedo (ssaal) n/d m*ndim*np*11
1682 ! aerosol asymmetry factor (asyal) n/d m*ndim*np*11
1683 ! cosecant of the solar zenith angle (csm) n/d m*n
1684 ! near ir surface albedo for beam fraction m*ndim
1685 ! radiation (rsirbm)
1686 ! near ir surface albedo for diffuse fraction m*ndim
1687 ! radiation (rsirdf)
1688 !
1689 !---- temporary array
1690 !
1691 ! scaled cloud optical thickness n/d m*n*np
1692 ! for beam radiation (tauclb)
1693 ! scaled cloud optical thickness n/d m*n*np
1694 ! for diffuse radiation (tauclf)
1695 !
1696 !----- output (updated) parameters:
1697 !
1698 ! all-sky flux (downward-upward) (flx) fraction m*ndim*(np+1)
1699 ! clear-sky flux (downward-upward) (flc) fraction m*ndim*(np+1)
1700 ! all-sky direct downward ir flux at
1701 ! the surface (fdirir) fraction m*ndim
1702 ! all-sky diffuse downward ir flux at
1703 ! the surface (fdifir) fraction m*ndim
1704 !
1705 !**********************************************************************
1706 implicit none
1707 !**********************************************************************
1708
1709 !-----input parameters
1710
1711 integer m,n,ndim,np
1712 integer ict(m,ndim),icb(m,ndim)
1713 real cwp(m,n,np,2),taucld(m,ndim,np,2),reff(m,ndim,np,2)
1714 real fcld(m,ndim,np),cc(m,n,3),cosz(m,ndim)
1715 real rsirbm(m,ndim),rsirdf(m,ndim)
1716 real taual(m,ndim,np,11),ssaal(m,ndim,np,11),asyal(m,ndim,np,11)
1717 real wh(m,n,np),csm(m,n)
1718 logical overcast,cldwater
1719
1720 !-----output (updated) parameters
1721
1722 real flx(m,ndim,np+1),flc(m,ndim,np+1)
1723 real flxu(m,ndim,np+1),flxd(m,ndim,np+1)
1724 real fdirir(m,ndim),fdifir(m,ndim)
1725
1726 !-----static parameters
1727
1728 integer nk,nband
1729 parameter (nk=10,nband=3)
1730 real xk(nk),hk(nband,nk),aib(nband,2),awb(nband,2)
1731 real aia(nband,3),awa(nband,3),aig(nband,3),awg(nband,3)
1732
1733 !-----temporary array
1734
1735 integer ib,iv,ik,i,j,k
1736 real tauclb(m,n,np),tauclf(m,n,np)
1737 real ssacl(m,n,np),asycl(m,n,np)
1738 real rr(m,n,np+1,2),tt(m,n,np+1,2),td(m,n,np+1,2), &
1739 rs(m,n,np+1,2),ts(m,n,np+1,2)
1740 real fall(m,n,np+1),fclr(m,n,np+1)
1741 real fallu(m,n,np+1),falld(m,n,np+1)
1742 real fsdir(m,n),fsdif(m,n)
1743
1744 real tauwv,tausto,ssatau,asysto,tauto,ssato,asyto
1745 real taux,reff1,reff2,w1,w2,g1,g2
1746 real ssaclt(m,n),asyclt(m,n)
1747 real rr1t(m,n),tt1t(m,n),td1t(m,n),rs1t(m,n),ts1t(m,n)
1748 real rr2t(m,n),tt2t(m,n),td2t(m,n),rs2t(m,n),ts2t(m,n)
1749
1750 !-----water vapor absorption coefficient for 10 k-intervals.
1751 ! unit: cm^2/gm
1752
1753 data xk/ &
1754 0.0010, 0.0133, 0.0422, 0.1334, 0.4217, &
1755 1.334, 5.623, 31.62, 177.8, 1000.0/
1756
1757 !-----water vapor k-distribution function,
1758 ! the sum of hk is 0.52926. unit: fraction
1759
1760 data hk/ &
1761 .20673,.08236,.01074, .03497,.01157,.00360, &
1762 .03011,.01133,.00411, .02260,.01143,.00421, &
1763 .01336,.01240,.00389, .00696,.01258,.00326, &
1764 .00441,.01381,.00499, .00115,.00650,.00465, &
1765 .00026,.00244,.00245, .00000,.00094,.00145/
1766
1767 !-----coefficients for computing the extinction coefficient of
1768 ! ice clouds from b=aib(*,1)+aib(*,2)/reff
1769
1770 data aib/ &
1771 .000333, .000333, .000333, &
1772 2.52, 2.52, 2.52/
1773
1774 !-----coefficients for computing the extinction coefficient of
1775 ! water clouds from b=awb(*,1)+awb(*,2)/reff
1776
1777 data awb/ &
1778 -0.0101, -0.0166, -0.0339, &
1779 1.72, 1.85, 2.16/
1780
1781
1782 !-----coefficients for computing the single scattering albedo of
1783 ! ice clouds from ssa=1-(aia(*,1)+aia(*,2)*reff+aia(*,3)*reff**2)
1784
1785 data aia/ &
1786 -.00000260, .00215346, .08938331, &
1787 .00000746, .00073709, .00299387, &
1788 .00000000,-.00000134,-.00001038/
1789
1790 !-----coefficients for computing the single scattering albedo of
1791 ! liquid clouds from ssa=1-(awa(*,1)+awa(*,2)*reff+awa(*,3)*reff**2)
1792
1793 data awa/ &
1794 .00000007,-.00019934, .01209318, &
1795 .00000845, .00088757, .01784739, &
1796 -.00000004,-.00000650,-.00036910/
1797
1798 !-----coefficients for computing the asymmetry factor of ice clouds
1799 ! from asycl=aig(*,1)+aig(*,2)*reff+aig(*,3)*reff**2
1800
1801 data aig/ &
1802 .74935228, .76098937, .84090400, &
1803 .00119715, .00141864, .00126222, &
1804 -.00000367,-.00000396,-.00000385/
1805
1806 !-----coefficients for computing the asymmetry factor of liquid clouds
1807 ! from asycl=awg(*,1)+awg(*,2)*reff+awg(*,3)*reff**2
1808
1809 data awg/ &
1810 .79375035, .74513197, .83530748, &
1811 .00832441, .01370071, .00257181, &
1812 -.00023263,-.00038203, .00005519/
1813
1814 !-----initialize surface fluxes, reflectances, and transmittances.
1815 ! cc is the maximum cloud cover in each of the three cloud groups.
1816
1817 do j= 1, n
1818 do i= 1, m
1819 fdirir(i,j)=0.0
1820 fdifir(i,j)=0.0
1821 rr(i,j,np+1,1)=rsirbm(i,j)
1822 rr(i,j,np+1,2)=rsirbm(i,j)
1823 rs(i,j,np+1,1)=rsirdf(i,j)
1824 rs(i,j,np+1,2)=rsirdf(i,j)
1825 td(i,j,np+1,1)=0.0
1826 td(i,j,np+1,2)=0.0
1827 tt(i,j,np+1,1)=0.0
1828 tt(i,j,np+1,2)=0.0
1829 ts(i,j,np+1,1)=0.0
1830 ts(i,j,np+1,2)=0.0
1831 cc(i,j,1)=0.0
1832 cc(i,j,2)=0.0
1833 cc(i,j,3)=0.0
1834 enddo
1835 enddo
1836
1837 !-----integration over spectral bands
1838
1839 do 100 ib=1,nband
1840
1841 iv=ib+8
1842
1843 !-----compute cloud optical thickness
1844
1845 if (cldwater) then
1846
1847 do k= 1, np
1848 do j= 1, n
1849 do i= 1, m
1850 taucld(i,j,k,1)=cwp(i,j,k,1)*(aib(ib,1) &
1851 +aib(ib,2)/reff(i,j,k,1))
1852 taucld(i,j,k,2)=cwp(i,j,k,2)*(awb(ib,1) &
1853 +awb(ib,2)/reff(i,j,k,2))
1854 enddo
1855 enddo
1856 enddo
1857
1858 endif
1859
1860 !-----options for scaling cloud optical thickness
1861
1862 if (overcast) then
1863
1864 do k= 1, np
1865 do j= 1, n
1866 do i= 1, m
1867 tauclb(i,j,k)=taucld(i,j,k,1)+taucld(i,j,k,2)
1868 tauclf(i,j,k)=tauclb(i,j,k)
1869 enddo
1870 enddo
1871 enddo
1872
1873 do k= 1, 3
1874 do j= 1, n
1875 do i= 1, m
1876 cc(i,j,k)=1.0
1877 enddo
1878 enddo
1879 enddo
1880
1881 else
1882
1883 !-----scale cloud optical thickness in each layer from taucld (with
1884 ! cloud amount fcld) to tauclb and tauclf (with cloud amount cc).
1885 ! tauclb is the scaled optical thickness for beam radiation and
1886 ! tauclf is for diffuse radiation.
1887
1888 call cldscale(m,n,ndim,np,cosz,fcld,taucld,ict,icb, &
1889 cc,tauclb,tauclf)
1890
1891 endif
1892
1893 !-----compute cloud single scattering albedo and asymmetry factor
1894 ! for a mixture of ice and liquid particles.
1895
1896 do k= 1, np
1897
1898 do j= 1, n
1899 do i= 1, m
1900
1901 ssaclt(i,j)=1.0
1902 asyclt(i,j)=1.0
1903
1904 taux=taucld(i,j,k,1)+taucld(i,j,k,2)
1905 if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then
1906
1907 reff1=min(reff(i,j,k,1),130.)
1908 reff2=min(reff(i,j,k,2),20.0)
1909
1910 w1=(1.-(aia(ib,1)+(aia(ib,2)+ &
1911 aia(ib,3)*reff1)*reff1))*taucld(i,j,k,1)
1912 w2=(1.-(awa(ib,1)+(awa(ib,2)+ &
1913 awa(ib,3)*reff2)*reff2))*taucld(i,j,k,2)
1914 ssaclt(i,j)=(w1+w2)/taux
1915
1916 g1=(aig(ib,1)+(aig(ib,2)+aig(ib,3)*reff1)*reff1)*w1
1917 g2=(awg(ib,1)+(awg(ib,2)+awg(ib,3)*reff2)*reff2)*w2
1918 asyclt(i,j)=(g1+g2)/(w1+w2)
1919
1920 endif
1921
1922 enddo
1923 enddo
1924
1925 do j=1,n
1926 do i=1,m
1927 ssacl(i,j,k)=ssaclt(i,j)
1928 enddo
1929 enddo
1930 do j=1,n
1931 do i=1,m
1932 asycl(i,j,k)=asyclt(i,j)
1933 enddo
1934 enddo
1935
1936 enddo
1937
1938 !-----integration over the k-distribution function
1939
1940 do 200 ik=1,nk
1941
1942 do 300 k= 1, np
1943
1944 do j= 1, n
1945 do i= 1, m
1946
1947 tauwv=xk(ik)*wh(i,j,k)
1948
1949 !-----compute clear-sky optical thickness, single scattering albedo,
1950 ! and asymmetry factor.
1951
1952 tausto=tauwv+taual(i,j,k,iv)+1.0e-8
1953 ssatau=ssaal(i,j,k,iv)*taual(i,j,k,iv)
1954 asysto=asyal(i,j,k,iv)*ssaal(i,j,k,iv)*taual(i,j,k,iv)
1955
1956 !-----compute reflectance and transmittance for cloudless layers
1957
1958 tauto=tausto
1959 ssato=ssatau/tauto+1.0e-8
1960
1961 if (ssato .gt. 0.001) then
1962
1963 ssato=min(ssato,0.999999)
1964 asyto=asysto/(ssato*tauto)
1965
1966 !- for direct incident radiation
1967
1968 call deledd (tauto,ssato,asyto,csm(i,j), &
1969 rr1t(i,j),tt1t(i,j),td1t(i,j))
1970
1971 !- for diffuse incident radiation
1972
1973 call sagpol (tauto,ssato,asyto,rs1t(i,j),ts1t(i,j))
1974
1975 else
1976
1977 td1t(i,j)=exp(-tauto*csm(i,j))
1978 ts1t(i,j)=exp(-1.66*tauto)
1979 tt1t(i,j)=0.0
1980 rr1t(i,j)=0.0
1981 rs1t(i,j)=0.0
1982
1983 endif
1984
1985 !-----compute reflectance and transmittance for cloud layers
1986
1987 if (tauclb(i,j,k).lt.0.01 .or. fcld(i,j,k).lt.0.01) then
1988
1989 rr2t(i,j)=rr1t(i,j)
1990 tt2t(i,j)=tt1t(i,j)
1991 td2t(i,j)=td1t(i,j)
1992 rs2t(i,j)=rs1t(i,j)
1993 ts2t(i,j)=ts1t(i,j)
1994
1995 else
1996
1997 !- for direct incident radiation
1998
1999 tauto=tausto+tauclb(i,j,k)
2000 ssato=(ssatau+ssacl(i,j,k)*tauclb(i,j,k))/tauto+1.0e-8
2001 ssato=min(ssato,0.999999)
2002 asyto=(asysto+asycl(i,j,k)*ssacl(i,j,k)*tauclb(i,j,k))/ &
2003 (ssato*tauto)
2004
2005 call deledd (tauto,ssato,asyto,csm(i,j), &
2006 rr2t(i,j),tt2t(i,j),td2t(i,j))
2007
2008 !- for diffuse incident radiation
2009
2010 tauto=tausto+tauclf(i,j,k)
2011 ssato=(ssatau+ssacl(i,j,k)*tauclf(i,j,k))/tauto+1.0e-8
2012 ssato=min(ssato,0.999999)
2013 asyto=(asysto+asycl(i,j,k)*ssacl(i,j,k)*tauclf(i,j,k))/ &
2014 (ssato*tauto)
2015
2016 call sagpol (tauto,ssato,asyto,rs2t(i,j),ts2t(i,j))
2017
2018 endif
2019
2020 enddo
2021 enddo
2022
2023 do j=1,n
2024 do i=1,m
2025 rr(i,j,k,1)=rr1t(i,j)
2026 enddo
2027 enddo
2028 do j=1,n
2029 do i=1,m
2030 tt(i,j,k,1)=tt1t(i,j)
2031 enddo
2032 enddo
2033 do j=1,n
2034 do i=1,m
2035 td(i,j,k,1)=td1t(i,j)
2036 enddo
2037 enddo
2038 do j=1,n
2039 do i=1,m
2040 rs(i,j,k,1)=rs1t(i,j)
2041 enddo
2042 enddo
2043 do j=1,n
2044 do i=1,m
2045 ts(i,j,k,1)=ts1t(i,j)
2046 enddo
2047 enddo
2048
2049 do j=1,n
2050 do i=1,m
2051 rr(i,j,k,2)=rr2t(i,j)
2052 enddo
2053 enddo
2054 do j=1,n
2055 do i=1,m
2056 tt(i,j,k,2)=tt2t(i,j)
2057 enddo
2058 enddo
2059 do j=1,n
2060 do i=1,m
2061 td(i,j,k,2)=td2t(i,j)
2062 enddo
2063 enddo
2064 do j=1,n
2065 do i=1,m
2066 rs(i,j,k,2)=rs2t(i,j)
2067 enddo
2068 enddo
2069 do j=1,n
2070 do i=1,m
2071 ts(i,j,k,2)=ts2t(i,j)
2072 enddo
2073 enddo
2074
2075 300 continue
2076
2077 !-----flux calculations
2078
2079 call cldflx (m,n,np,ict,icb,overcast,cc,rr,tt,td,rs,ts, &
2080 fclr,fall,fallu,falld,fsdir,fsdif)
2081
2082 do k= 1, np+1
2083 do j= 1, n
2084 do i= 1, m
2085 flx(i,j,k) = flx(i,j,k)+fall(i,j,k)*hk(ib,ik)
2086 flxu(i,j,k) = flxu(i,j,k)+fallu(i,j,k)*hk(ib,ik)
2087 flxd(i,j,k) = flxd(i,j,k)+falld(i,j,k)*hk(ib,ik)
2088 enddo
2089 enddo
2090 do j= 1, n
2091 do i= 1, m
2092 flc(i,j,k) = flc(i,j,k)+fclr(i,j,k)*hk(ib,ik)
2093 enddo
2094 enddo
2095 enddo
2096
2097 !-----compute downward surface fluxes in the ir region
2098
2099 do j= 1, n
2100 do i= 1, m
2101 fdirir(i,j) = fdirir(i,j)+fsdir(i,j)*hk(ib,ik)
2102 fdifir(i,j) = fdifir(i,j)+fsdif(i,j)*hk(ib,ik)
2103 enddo
2104 enddo
2105
2106 200 continue
2107 100 continue
2108
2109 end subroutine solir
2110
2111 !********************************************************************
2112
2113 subroutine cldscale (m,n,ndim,np,cosz,fcld,taucld,ict,icb, &
2114 cc,tauclb,tauclf)
2115
2116 !********************************************************************
2117 !
2118 ! This subroutine computes the high, middle, and
2119 ! low cloud amounts and scales the cloud optical thickness.
2120 !
2121 ! To simplify calculations in a cloudy atmosphere, clouds are
2122 ! grouped into high, middle and low clouds separated by the levels
2123 ! ict and icb (level 1 is the top of the model atmosphere).
2124 !
2125 ! Within each of the three groups, clouds are assumed maximally
2126 ! overlapped, and the cloud cover (cc) of a group is the maximum
2127 ! cloud cover of all the layers in the group. The optical thickness
2128 ! (taucld) of a given layer is then scaled to new values (tauclb and
2129 ! tauclf) so that the layer reflectance corresponding to the cloud
2130 ! cover cc is the same as the original reflectance with optical
2131 ! thickness taucld and cloud cover fcld.
2132 !
2133 !---input parameters
2134 !
2135 ! number of grid intervals in zonal direction (m)
2136 ! number of grid intervals in meridional direction (n)
2137 ! maximum number of grid intervals in meridional direction (ndim)
2138 ! number of atmospheric layers (np)
2139 ! cosine of the solar zenith angle (cosz)
2140 ! fractional cloud cover (fcld)
2141 ! cloud optical thickness (taucld)
2142 ! index separating high and middle clouds (ict)
2143 ! index separating middle and low clouds (icb)
2144 !
2145 !---output parameters
2146 !
2147 ! fractional cover of high, middle, and low clouds (cc)
2148 ! scaled cloud optical thickness for beam radiation (tauclb)
2149 ! scaled cloud optical thickness for diffuse radiation (tauclf)
2150 !
2151 !********************************************************************
2152 implicit none
2153 !********************************************************************
2154
2155 !-----input parameters
2156
2157 integer m,n,ndim,np
2158 integer ict(m,ndim),icb(m,ndim)
2159 real cosz(m,ndim),fcld(m,ndim,np),taucld(m,ndim,np,2)
2160
2161 !-----output parameters
2162
2163 real cc(m,n,3),tauclb(m,n,np),tauclf(m,n,np)
2164
2165 !-----temporary variables
2166
2167 integer i,j,k,im,it,ia,kk
2168 real fm,ft,fa,xai,taux
2169
2170 !-----pre-computed table
2171
2172 integer nm,nt,na
2173 parameter (nm=11,nt=9,na=11)
2174 real dm,dt,da,t1,caib(nm,nt,na),caif(nt,na)
2175 parameter (dm=0.1,dt=0.30103,da=0.1,t1=-0.9031)
2176
2177 !-----include the pre-computed table of mcai for scaling the cloud optical
2178 ! thickness under the assumption that clouds are maximally overlapped
2179 !
2180 ! caib is for scaling the cloud optical thickness for direct radiation
2181 ! caif is for scaling the cloud optical thickness for diffuse radiation
2182
2183
2184 data ((caib(1,i,j),j=1,11),i=1,9)/ &
2185 .000,0.068,0.140,0.216,0.298,0.385,0.481,0.586,0.705,0.840,1.000, &
2186 .000,0.052,0.106,0.166,0.230,0.302,0.383,0.478,0.595,0.752,1.000, &
2187 .000,0.038,0.078,0.120,0.166,0.218,0.276,0.346,0.438,0.582,1.000, &
2188 .000,0.030,0.060,0.092,0.126,0.164,0.206,0.255,0.322,0.442,1.000, &
2189 .000,0.025,0.051,0.078,0.106,0.136,0.170,0.209,0.266,0.462,1.000, &
2190 .000,0.023,0.046,0.070,0.095,0.122,0.150,0.187,0.278,0.577,1.000, &
2191 .000,0.022,0.043,0.066,0.089,0.114,0.141,0.187,0.354,0.603,1.000, &
2192 .000,0.021,0.042,0.063,0.086,0.108,0.135,0.214,0.349,0.565,1.000, &
2193 .000,0.021,0.041,0.062,0.083,0.105,0.134,0.202,0.302,0.479,1.000/
2194 data ((caib(2,i,j),j=1,11),i=1,9)/ &
2195 .000,0.088,0.179,0.272,0.367,0.465,0.566,0.669,0.776,0.886,1.000, &
2196 .000,0.079,0.161,0.247,0.337,0.431,0.531,0.637,0.749,0.870,1.000, &
2197 .000,0.065,0.134,0.207,0.286,0.372,0.466,0.572,0.692,0.831,1.000, &
2198 .000,0.049,0.102,0.158,0.221,0.290,0.370,0.465,0.583,0.745,1.000, &
2199 .000,0.037,0.076,0.118,0.165,0.217,0.278,0.354,0.459,0.638,1.000, &
2200 .000,0.030,0.061,0.094,0.130,0.171,0.221,0.286,0.398,0.631,1.000, &
2201 .000,0.026,0.052,0.081,0.111,0.146,0.189,0.259,0.407,0.643,1.000, &
2202 .000,0.023,0.047,0.072,0.098,0.129,0.170,0.250,0.387,0.598,1.000, &
2203 .000,0.022,0.044,0.066,0.090,0.118,0.156,0.224,0.328,0.508,1.000/
2204 data ((caib(3,i,j),j=1,11),i=1,9)/ &
2205 .000,0.094,0.189,0.285,0.383,0.482,0.582,0.685,0.788,0.894,1.000, &
2206 .000,0.088,0.178,0.271,0.366,0.465,0.565,0.669,0.776,0.886,1.000, &
2207 .000,0.079,0.161,0.247,0.337,0.431,0.531,0.637,0.750,0.870,1.000, &
2208 .000,0.066,0.134,0.209,0.289,0.375,0.470,0.577,0.697,0.835,1.000, &
2209 .000,0.050,0.104,0.163,0.227,0.300,0.383,0.483,0.606,0.770,1.000, &
2210 .000,0.038,0.080,0.125,0.175,0.233,0.302,0.391,0.518,0.710,1.000, &
2211 .000,0.031,0.064,0.100,0.141,0.188,0.249,0.336,0.476,0.689,1.000, &
2212 .000,0.026,0.054,0.084,0.118,0.158,0.213,0.298,0.433,0.638,1.000, &
2213 .000,0.023,0.048,0.074,0.102,0.136,0.182,0.254,0.360,0.542,1.000/
2214 data ((caib(4,i,j),j=1,11),i=1,9)/ &
2215 .000,0.096,0.193,0.290,0.389,0.488,0.589,0.690,0.792,0.896,1.000, &
2216 .000,0.092,0.186,0.281,0.378,0.477,0.578,0.680,0.785,0.891,1.000, &
2217 .000,0.086,0.174,0.264,0.358,0.455,0.556,0.660,0.769,0.882,1.000, &
2218 .000,0.074,0.153,0.235,0.323,0.416,0.514,0.622,0.737,0.862,1.000, &
2219 .000,0.061,0.126,0.195,0.271,0.355,0.449,0.555,0.678,0.823,1.000, &
2220 .000,0.047,0.098,0.153,0.215,0.286,0.370,0.471,0.600,0.770,1.000, &
2221 .000,0.037,0.077,0.120,0.170,0.230,0.303,0.401,0.537,0.729,1.000, &
2222 .000,0.030,0.062,0.098,0.138,0.187,0.252,0.343,0.476,0.673,1.000, &
2223 .000,0.026,0.053,0.082,0.114,0.154,0.207,0.282,0.391,0.574,1.000/
2224 data ((caib(5,i,j),j=1,11),i=1,9)/ &
2225 .000,0.097,0.194,0.293,0.392,0.492,0.592,0.693,0.794,0.897,1.000, &
2226 .000,0.094,0.190,0.286,0.384,0.483,0.584,0.686,0.789,0.894,1.000, &
2227 .000,0.090,0.181,0.274,0.370,0.468,0.569,0.672,0.778,0.887,1.000, &
2228 .000,0.081,0.165,0.252,0.343,0.439,0.539,0.645,0.757,0.874,1.000, &
2229 .000,0.069,0.142,0.218,0.302,0.392,0.490,0.598,0.717,0.850,1.000, &
2230 .000,0.054,0.114,0.178,0.250,0.330,0.422,0.529,0.656,0.810,1.000, &
2231 .000,0.042,0.090,0.141,0.200,0.269,0.351,0.455,0.589,0.764,1.000, &
2232 .000,0.034,0.070,0.112,0.159,0.217,0.289,0.384,0.515,0.703,1.000, &
2233 .000,0.028,0.058,0.090,0.128,0.174,0.231,0.309,0.420,0.602,1.000/
2234 data ((caib(6,i,j),j=1,11),i=1,9)/ &
2235 .000,0.098,0.196,0.295,0.394,0.494,0.594,0.695,0.796,0.898,1.000, &
2236 .000,0.096,0.193,0.290,0.389,0.488,0.588,0.690,0.792,0.895,1.000, &
2237 .000,0.092,0.186,0.281,0.378,0.477,0.577,0.680,0.784,0.891,1.000, &
2238 .000,0.086,0.174,0.264,0.358,0.455,0.556,0.661,0.769,0.882,1.000, &
2239 .000,0.075,0.154,0.237,0.325,0.419,0.518,0.626,0.741,0.865,1.000, &
2240 .000,0.062,0.129,0.201,0.279,0.366,0.462,0.571,0.694,0.836,1.000, &
2241 .000,0.049,0.102,0.162,0.229,0.305,0.394,0.501,0.631,0.793,1.000, &
2242 .000,0.038,0.080,0.127,0.182,0.245,0.323,0.422,0.550,0.730,1.000, &
2243 .000,0.030,0.064,0.100,0.142,0.192,0.254,0.334,0.448,0.627,1.000/
2244 data ((caib(7,i,j),j=1,11),i=1,9)/ &
2245 .000,0.098,0.198,0.296,0.396,0.496,0.596,0.696,0.797,0.898,1.000, &
2246 .000,0.097,0.194,0.293,0.392,0.491,0.591,0.693,0.794,0.897,1.000, &
2247 .000,0.094,0.190,0.286,0.384,0.483,0.583,0.686,0.789,0.894,1.000, &
2248 .000,0.089,0.180,0.274,0.369,0.467,0.568,0.672,0.778,0.887,1.000, &
2249 .000,0.081,0.165,0.252,0.344,0.440,0.541,0.646,0.758,0.875,1.000, &
2250 .000,0.069,0.142,0.221,0.306,0.397,0.496,0.604,0.722,0.854,1.000, &
2251 .000,0.056,0.116,0.182,0.256,0.338,0.432,0.540,0.666,0.816,1.000, &
2252 .000,0.043,0.090,0.143,0.203,0.273,0.355,0.455,0.583,0.754,1.000, &
2253 .000,0.034,0.070,0.111,0.157,0.210,0.276,0.359,0.474,0.650,1.000/
2254 data ((caib(8,i,j),j=1,11),i=1,9)/ &
2255 .000,0.099,0.198,0.298,0.398,0.497,0.598,0.698,0.798,0.899,1.000, &
2256 .000,0.098,0.196,0.295,0.394,0.494,0.594,0.695,0.796,0.898,1.000, &
2257 .000,0.096,0.193,0.290,0.390,0.489,0.589,0.690,0.793,0.896,1.000, &
2258 .000,0.093,0.186,0.282,0.379,0.478,0.578,0.681,0.786,0.892,1.000, &
2259 .000,0.086,0.175,0.266,0.361,0.458,0.558,0.663,0.771,0.883,1.000, &
2260 .000,0.076,0.156,0.240,0.330,0.423,0.523,0.630,0.744,0.867,1.000, &
2261 .000,0.063,0.130,0.203,0.282,0.369,0.465,0.572,0.694,0.834,1.000, &
2262 .000,0.049,0.102,0.161,0.226,0.299,0.385,0.486,0.611,0.774,1.000, &
2263 .000,0.038,0.078,0.122,0.172,0.229,0.297,0.382,0.498,0.672,1.000/
2264 data ((caib(9,i,j),j=1,11),i=1,9)/ &
2265 .000,0.099,0.199,0.298,0.398,0.498,0.598,0.699,0.799,0.899,1.000, &
2266 .000,0.099,0.198,0.298,0.398,0.497,0.598,0.698,0.798,0.899,1.000, &
2267 .000,0.098,0.196,0.295,0.394,0.494,0.594,0.695,0.796,0.898,1.000, &
2268 .000,0.096,0.193,0.290,0.389,0.488,0.588,0.690,0.792,0.895,1.000, &
2269 .000,0.092,0.185,0.280,0.376,0.474,0.575,0.678,0.782,0.890,1.000, &
2270 .000,0.084,0.170,0.259,0.351,0.447,0.547,0.652,0.762,0.878,1.000, &
2271 .000,0.071,0.146,0.224,0.308,0.398,0.494,0.601,0.718,0.850,1.000, &
2272 .000,0.056,0.114,0.178,0.248,0.325,0.412,0.514,0.638,0.793,1.000, &
2273 .000,0.042,0.086,0.134,0.186,0.246,0.318,0.405,0.521,0.691,1.000/
2274 data ((caib(10,i,j),j=1,11),i=1,9)/ &
2275 .000,0.100,0.200,0.300,0.400,0.500,0.600,0.700,0.800,0.900,1.000, &
2276 .000,0.100,0.200,0.300,0.400,0.500,0.600,0.700,0.800,0.900,1.000, &
2277 .000,0.100,0.200,0.300,0.400,0.500,0.600,0.700,0.800,0.900,1.000, &
2278 .000,0.100,0.199,0.298,0.398,0.498,0.598,0.698,0.798,0.899,1.000, &
2279 .000,0.098,0.196,0.294,0.392,0.491,0.590,0.691,0.793,0.896,1.000, &
2280 .000,0.092,0.185,0.278,0.374,0.470,0.570,0.671,0.777,0.886,1.000, &
2281 .000,0.081,0.162,0.246,0.333,0.424,0.521,0.625,0.738,0.862,1.000, &
2282 .000,0.063,0.128,0.196,0.270,0.349,0.438,0.540,0.661,0.809,1.000, &
2283 .000,0.046,0.094,0.146,0.202,0.264,0.337,0.426,0.542,0.710,1.000/
2284 data ((caib(11,i,j),j=1,11),i=1,9)/ &
2285 .000,0.101,0.202,0.302,0.402,0.502,0.602,0.702,0.802,0.901,1.000, &
2286 .000,0.102,0.202,0.303,0.404,0.504,0.604,0.703,0.802,0.902,1.000, &
2287 .000,0.102,0.205,0.306,0.406,0.506,0.606,0.706,0.804,0.902,1.000, &
2288 .000,0.104,0.207,0.309,0.410,0.510,0.609,0.707,0.805,0.902,1.000, &
2289 .000,0.106,0.208,0.309,0.409,0.508,0.606,0.705,0.803,0.902,1.000, &
2290 .000,0.102,0.202,0.298,0.395,0.493,0.590,0.690,0.790,0.894,1.000, &
2291 .000,0.091,0.179,0.267,0.357,0.449,0.545,0.647,0.755,0.872,1.000, &
2292 .000,0.073,0.142,0.214,0.290,0.372,0.462,0.563,0.681,0.822,1.000, &
2293 .000,0.053,0.104,0.158,0.217,0.281,0.356,0.446,0.562,0.726,1.000/
2294 data ((caif(i,j),j=1,11),i=1,9)/ &
2295 .000,0.099,0.198,0.297,0.397,0.496,0.597,0.697,0.798,0.899,1.000, &
2296 .000,0.098,0.196,0.294,0.394,0.494,0.594,0.694,0.796,0.898,1.000, &
2297 .000,0.096,0.192,0.290,0.388,0.487,0.587,0.689,0.792,0.895,1.000, &
2298 .000,0.092,0.185,0.280,0.376,0.476,0.576,0.678,0.783,0.890,1.000, &
2299 .000,0.085,0.173,0.263,0.357,0.454,0.555,0.659,0.768,0.881,1.000, &
2300 .000,0.076,0.154,0.237,0.324,0.418,0.517,0.624,0.738,0.864,1.000, &
2301 .000,0.063,0.131,0.203,0.281,0.366,0.461,0.567,0.688,0.830,1.000, &
2302 .000,0.052,0.107,0.166,0.232,0.305,0.389,0.488,0.610,0.770,1.000, &
2303 .000,0.043,0.088,0.136,0.189,0.248,0.317,0.400,0.510,0.675,1.000/
2304
2305 !-----clouds within each of the high, middle, and low clouds are assumed
2306 ! to be maximally overlapped, and the cloud cover (cc) for a group
2307 ! (high, middle, or low) is the maximum cloud cover of all the layers
2308 ! within a group
2309
2310 do j=1,n
2311 do i=1,m
2312 cc(i,j,1)=0.0
2313 cc(i,j,2)=0.0
2314 cc(i,j,3)=0.0
2315 enddo
2316 enddo
2317 do j=1,n
2318 do i=1,m
2319 do k=1,ict(i,j)-1
2320 cc(i,j,1)=max(cc(i,j,1),fcld(i,j,k))
2321 enddo
2322 enddo
2323 enddo
2324
2325 do j=1,n
2326 do i=1,m
2327 do k=ict(i,j),icb(i,j)-1
2328 cc(i,j,2)=max(cc(i,j,2),fcld(i,j,k))
2329 enddo
2330 enddo
2331 enddo
2332
2333 do j=1,n
2334 do i=1,m
2335 do k=icb(i,j),np
2336 cc(i,j,3)=max(cc(i,j,3),fcld(i,j,k))
2337 enddo
2338 enddo
2339 enddo
2340
2341 !-----scale the cloud optical thickness.
2342 ! taucld(i,j,k,1) is the optical thickness for ice particles, and
2343 ! taucld(i,j,k,2) is the optical thickness for liquid particles.
2344
2345 do j=1,n
2346 do i=1,m
2347
2348 do k=1,np
2349
2350 if(k.lt.ict(i,j)) then
2351 kk=1
2352 elseif(k.ge.ict(i,j) .and. k.lt.icb(i,j)) then
2353 kk=2
2354 else
2355 kk=3
2356 endif
2357
2358 tauclb(i,j,k) = 0.0
2359 tauclf(i,j,k) = 0.0
2360
2361 taux=taucld(i,j,k,1)+taucld(i,j,k,2)
2362 if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then
2363
2364 !-----normalize cloud cover
2365
2366 fa=fcld(i,j,k)/cc(i,j,kk)
2367
2368 !-----table look-up
2369
2370 taux=min(taux,32.)
2371
2372 fm=cosz(i,j)/dm
2373 ft=(log10(taux)-t1)/dt
2374 fa=fa/da
2375
2376 im=int(fm+1.5)
2377 it=int(ft+1.5)
2378 ia=int(fa+1.5)
2379
2380 im=max(im,2)
2381 it=max(it,2)
2382 ia=max(ia,2)
2383
2384 im=min(im,nm-1)
2385 it=min(it,nt-1)
2386 ia=min(ia,na-1)
2387
2388 fm=fm-float(im-1)
2389 ft=ft-float(it-1)
2390 fa=fa-float(ia-1)
2391
2392 !-----scale cloud optical thickness for beam radiation.
2393 ! the scaling factor, xai, is a function of the solar zenith
2394 ! angle, optical thickness, and cloud cover.
2395
2396 xai= (-caib(im-1,it,ia)*(1.-fm)+ &
2397 caib(im+1,it,ia)*(1.+fm))*fm*.5+caib(im,it,ia)*(1.-fm*fm)
2398
2399 xai=xai+(-caib(im,it-1,ia)*(1.-ft)+ &
2400 caib(im,it+1,ia)*(1.+ft))*ft*.5+caib(im,it,ia)*(1.-ft*ft)
2401
2402 xai=xai+(-caib(im,it,ia-1)*(1.-fa)+ &
2403 caib(im,it,ia+1)*(1.+fa))*fa*.5+caib(im,it,ia)*(1.-fa*fa)
2404
2405 xai= xai-2.*caib(im,it,ia)
2406 xai=max(xai,0.0)
2407
2408 tauclb(i,j,k) = taux*xai
2409
2410 !-----scale cloud optical thickness for diffuse radiation.
2411 ! the scaling factor, xai, is a function of the cloud optical
2412 ! thickness and cover but not the solar zenith angle.
2413
2414 xai= (-caif(it-1,ia)*(1.-ft)+ &
2415 caif(it+1,ia)*(1.+ft))*ft*.5+caif(it,ia)*(1.-ft*ft)
2416
2417 xai=xai+(-caif(it,ia-1)*(1.-fa)+ &
2418 caif(it,ia+1)*(1.+fa))*fa*.5+caif(it,ia)*(1.-fa*fa)
2419
2420 xai= xai-caif(it,ia)
2421 xai=max(xai,0.0)
2422
2423 tauclf(i,j,k) = taux*xai
2424
2425 endif
2426
2427 enddo
2428 enddo
2429 enddo
2430
2431 end subroutine cldscale
2432
2433 !*********************************************************************
2434
2435 subroutine deledd(tau,ssc,g0,csm,rr,tt,td)
2436
2437 !*********************************************************************
2438 !
2439 !-----uses the delta-eddington approximation to compute the
2440 ! bulk scattering properties of a single layer
2441 ! coded following King and Harshvardhan (JAS, 1986)
2442 !
2443 ! inputs:
2444 !
2445 ! tau: the effective optical thickness
2446 ! ssc: the effective single scattering albedo
2447 ! g0: the effective asymmetry factor
2448 ! csm: the effective secant of the zenith angle
2449 !
2450 ! outputs:
2451 !
2452 ! rr: the layer reflection of the direct beam
2453 ! tt: the layer diffuse transmission of the direct beam
2454 ! td: the layer direct transmission of the direct beam
2455 !
2456 !*********************************************************************
2457 implicit none
2458 !*********************************************************************
2459
2460 real zero,one,two,three,four,fourth,seven,thresh
2461 parameter (one =1., three=3.)
2462 parameter (two =2., seven=7.)
2463 parameter (four=4., fourth=.25)
2464 parameter (zero=0., thresh=1.e-8)
2465
2466 !-----input parameters
2467 real tau,ssc,g0,csm
2468
2469 !-----output parameters
2470 real rr,tt,td
2471
2472 !-----temporary parameters
2473
2474 real zth,ff,xx,taup,sscp,gp,gm1,gm2,gm3,akk,alf1,alf2, &
2475 all,bll,st7,st8,cll,dll,fll,ell,st1,st2,st3,st4
2476
2477 !---------------------------------------------------------------------
2478
2479 zth = one / csm
2480
2481 ! delta-eddington scaling of single scattering albedo,
2482 ! optical thickness, and asymmetry factor,
2483 ! K & H eqs(27-29)
2484
2485 ff = g0*g0
2486 xx = one-ff*ssc
2487 taup= tau*xx
2488 sscp= ssc*(one-ff)/xx
2489 gp = g0/(one+g0)
2490
2491 ! gamma1, gamma2, and gamma3. see table 2 and eq(26) K & H
2492 ! ssc and gp are the d-s single scattering
2493 ! albedo and asymmetry factor.
2494
2495 xx = three*gp
2496 gm1 = (seven - sscp*(four+xx))*fourth
2497 gm2 = -(one - sscp*(four-xx))*fourth
2498
2499 ! akk is k as defined in eq(25) of K & H
2500
2501 akk = sqrt((gm1+gm2)*(gm1-gm2))
2502
2503 xx = akk * zth
2504 st7 = one - xx
2505 st8 = one + xx
2506 st3 = st7 * st8
2507
2508 if (abs(st3) .lt. thresh) then
2509 zth = zth + 0.001
2510 xx = akk * zth
2511 st7 = one - xx
2512 st8 = one + xx
2513 st3 = st7 * st8
2514 endif
2515
2516 ! extinction of the direct beam transmission
2517
2518 td = exp(-taup/zth)
2519
2520 ! alf1 and alf2 are alpha1 and alpha2 from eqs (23) & (24) of K & H
2521
2522 gm3 = (two - zth*three*gp)*fourth
2523 xx = gm1 - gm2
2524 alf1 = gm1 - gm3 * xx
2525 alf2 = gm2 + gm3 * xx
2526
2527 ! all is last term in eq(21) of K & H
2528 ! bll is last term in eq(22) of K & H
2529
2530 xx = akk * two
2531 all = (gm3 - alf2 * zth )*xx*td
2532 bll = (one - gm3 + alf1*zth)*xx
2533
2534 xx = akk * gm3
2535 cll = (alf2 + xx) * st7
2536 dll = (alf2 - xx) * st8
2537
2538 xx = akk * (one-gm3)
2539 fll = (alf1 + xx) * st8
2540 ell = (alf1 - xx) * st7
2541
2542 st2 = exp(-akk*taup)
2543 st4 = st2 * st2
2544
2545 st1 = sscp / ((akk+gm1 + (akk-gm1)*st4) * st3)
2546
2547 ! rr is r-hat of eq(21) of K & H
2548 ! tt is diffuse part of t-hat of eq(22) of K & H
2549
2550 rr = ( cll-dll*st4 -all*st2)*st1
2551 tt = - ((fll-ell*st4)*td-bll*st2)*st1
2552
2553 rr = max(rr,zero)
2554 tt = max(tt,zero)
2555
2556 end subroutine deledd
2557
2558 !*********************************************************************
2559
2560 subroutine sagpol(tau,ssc,g0,rll,tll)
2561
2562 !*********************************************************************
2563 !-----transmittance (tll) and reflectance (rll) of diffuse radiation
2564 ! follows Sagan and Pollock (JGR, 1967).
2565 ! also, eq.(31) of Lacis and Hansen (JAS, 1974).
2566 !
2567 !-----input parameters:
2568 !
2569 ! tau: the effective optical thickness
2570 ! ssc: the effective single scattering albedo
2571 ! g0: the effective asymmetry factor
2572 !
2573 !-----output parameters:
2574 !
2575 ! rll: the layer reflection of diffuse radiation
2576 ! tll: the layer transmission of diffuse radiation
2577 !
2578 !*********************************************************************
2579 implicit none
2580 !*********************************************************************
2581
2582 real one,three,four
2583 parameter (one=1., three=3., four=4.)
2584
2585 !-----output parameters:
2586
2587 real tau,ssc,g0
2588
2589 !-----output parameters:
2590
2591 real rll,tll
2592
2593 !-----temporary arrays
2594
2595 real xx,uuu,ttt,emt,up1,um1,st1
2596
2597 xx = one-ssc*g0
2598 uuu = sqrt( xx/(one-ssc))
2599 ttt = sqrt( xx*(one-ssc)*three )*tau
2600 emt = exp(-ttt)
2601 up1 = uuu + one
2602 um1 = uuu - one
2603 xx = um1*emt
2604 st1 = one / ((up1+xx) * (up1-xx))
2605 rll = up1*um1*(one-emt*emt)*st1
2606 tll = uuu*four*emt *st1
2607
2608 end subroutine sagpol
2609
2610 !*******************************************************************
2611
2612 subroutine cldflx (m,n,np,ict,icb,overcast,cc,rr,tt,td,rs,ts,&
2613 fclr,fall,fallu,falld,fsdir,fsdif)
2614
2615 !*******************************************************************
2616 ! compute upward and downward fluxes using a two-stream adding method
2617 ! following equations (3)-(5) of Chou (1992, JAS).
2618 !
2619 ! clouds are grouped into high, middle, and low clouds which are
2620 ! assumed randomly overlapped. It involves eight sets of calculations.
2621 ! In each set of calculations, each atmospheric layer is homogeneous,
2622 ! either totally filled with clouds or without clouds.
2623
2624 ! input parameters:
2625 !
2626 ! m: number of soundings in zonal direction
2627 ! n: number of soundings in meridional direction
2628 ! np: number of atmospheric layers
2629 ! ict: the level separating high and middle clouds
2630 ! icb: the level separating middle and low clouds
2631 ! cc: effective cloud covers for high, middle and low clouds
2632 ! tt: diffuse transmission of a layer illuminated by beam radiation
2633 ! td: direct beam tranmssion
2634 ! ts: transmission of a layer illuminated by diffuse radiation
2635 ! rr: reflection of a layer illuminated by beam radiation
2636 ! rs: reflection of a layer illuminated by diffuse radiation
2637 !
2638 ! output parameters:
2639 !
2640 ! fclr: clear-sky flux (downward minus upward)
2641 ! fall: all-sky flux (downward minus upward)
2642 ! fsdir: surface direct downward flux
2643 ! fsdif: surface diffuse downward flux
2644 !
2645 !*********************************************************************c
2646 implicit none
2647 !*********************************************************************c
2648
2649 !-----input parameters
2650
2651 integer m,n,np
2652 integer ict(m,n),icb(m,n)
2653
2654 real rr(m,n,np+1,2),tt(m,n,np+1,2),td(m,n,np+1,2)
2655 real rs(m,n,np+1,2),ts(m,n,np+1,2)
2656 real cc(m,n,3)
2657 logical overcast
2658
2659 !-----temporary array
2660
2661 integer i,j,k,ih,im,is,itm
2662 real rra(m,n,np+1,2,2),tta(m,n,np+1,2,2),tda(m,n,np+1,2,2)
2663 real rsa(m,n,np+1,2,2),rxa(m,n,np+1,2,2)
2664 real ch(m,n),cm(m,n),ct(m,n),flxdn(m,n,np+1)
2665 real flxdnu(m,n,np+1),flxdnd(m,n,np+1)
2666 real fdndir(m,n),fdndif(m,n),fupdif
2667 real denm,xx
2668
2669 !-----output parameters
2670
2671 real fclr(m,n,np+1),fall(m,n,np+1)
2672 real fallu(m,n,np+1),falld(m,n,np+1)
2673 real fsdir(m,n),fsdif(m,n)
2674
2675 !-----initialize all-sky flux (fall) and surface downward fluxes
2676
2677 do k=1,np+1
2678 do j=1,n
2679 do i=1,m
2680 fclr(i,j,k)=0.0
2681 fall(i,j,k)=0.0
2682 fallu(i,j,k)=0.0
2683 falld(i,j,k)=0.0
2684 enddo
2685 enddo
2686 enddo
2687
2688 do j=1,n
2689 do i=1,m
2690 fsdir(i,j)=0.0
2691 fsdif(i,j)=0.0
2692 enddo
2693 enddo
2694
2695 !-----compute transmittances and reflectances for a composite of
2696 ! layers. layers are added one at a time, going down from the top.
2697 ! tda is the composite transmittance illuminated by beam radiation
2698 ! tta is the composite diffuse transmittance illuminated by
2699 ! beam radiation
2700 ! rsa is the composite reflectance illuminated from below
2701 ! by diffuse radiation
2702 ! tta and rsa are computed from eqs. (4b) and (3b) of Chou
2703
2704 itm=1
2705
2706 !-----if overcas.=.true., set itm=2, and only one set of fluxes is computed
2707
2708 if (overcast) itm=2
2709
2710 !-----for high clouds. indices 1 and 2 denote clear and cloudy
2711 ! situations, respectively.
2712
2713 do 10 ih=itm,2
2714
2715 do j= 1, n
2716 do i= 1, m
2717 tda(i,j,1,ih,1)=td(i,j,1,ih)
2718 tta(i,j,1,ih,1)=tt(i,j,1,ih)
2719 rsa(i,j,1,ih,1)=rs(i,j,1,ih)
2720 tda(i,j,1,ih,2)=td(i,j,1,ih)
2721 tta(i,j,1,ih,2)=tt(i,j,1,ih)
2722 rsa(i,j,1,ih,2)=rs(i,j,1,ih)
2723 enddo
2724 enddo
2725
2726 do j= 1, n
2727 do i= 1, m
2728 do k= 2, ict(i,j)-1
2729 denm = ts(i,j,k,ih)/( 1.-rsa(i,j,k-1,ih,1)*rs(i,j,k,ih))
2730 tda(i,j,k,ih,1)= tda(i,j,k-1,ih,1)*td(i,j,k,ih)
2731 tta(i,j,k,ih,1)= tda(i,j,k-1,ih,1)*tt(i,j,k,ih) &
2732 +(tda(i,j,k-1,ih,1)*rr(i,j,k,ih) &
2733 *rsa(i,j,k-1,ih,1)+tta(i,j,k-1,ih,1))*denm
2734 rsa(i,j,k,ih,1)= rs(i,j,k,ih)+ts(i,j,k,ih) &
2735 *rsa(i,j,k-1,ih,1)*denm
2736 tda(i,j,k,ih,2)= tda(i,j,k,ih,1)
2737 tta(i,j,k,ih,2)= tta(i,j,k,ih,1)
2738 rsa(i,j,k,ih,2)= rsa(i,j,k,ih,1)
2739 enddo
2740 enddo
2741 enddo
2742
2743 !-----for middle clouds
2744
2745 do 10 im=itm,2
2746
2747 do j= 1, n
2748 do i= 1, m
2749 do k= ict(i,j), icb(i,j)-1
2750 denm = ts(i,j,k,im)/( 1.-rsa(i,j,k-1,ih,im)*rs(i,j,k,im))
2751 tda(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*td(i,j,k,im)
2752 tta(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*tt(i,j,k,im) &
2753 +(tda(i,j,k-1,ih,im)*rr(i,j,k,im) &
2754 *rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm
2755 rsa(i,j,k,ih,im)= rs(i,j,k,im)+ts(i,j,k,im) &
2756 *rsa(i,j,k-1,ih,im)*denm
2757 enddo
2758 enddo
2759 enddo
2760
2761 10 continue
2762
2763 !-----layers are added one at a time, going up from the surface.
2764 ! rra is the composite reflectance illuminated by beam radiation
2765 ! rxa is the composite reflectance illuminated from above
2766 ! by diffuse radiation
2767 ! rra and rxa are computed from eqs. (4a) and (3a) of Chou
2768
2769 !-----for the low clouds
2770
2771 do 20 is=itm,2
2772
2773 do j= 1, n
2774 do i= 1, m
2775 rra(i,j,np+1,1,is)=rr(i,j,np+1,is)
2776 rxa(i,j,np+1,1,is)=rs(i,j,np+1,is)
2777 rra(i,j,np+1,2,is)=rr(i,j,np+1,is)
2778 rxa(i,j,np+1,2,is)=rs(i,j,np+1,is)
2779 enddo
2780 enddo
2781
2782 do j= 1, n
2783 do i= 1, m
2784 do k=np,icb(i,j),-1
2785 denm=ts(i,j,k,is)/( 1.-rs(i,j,k,is)*rxa(i,j,k+1,1,is) )
2786 rra(i,j,k,1,is)=rr(i,j,k,is)+(td(i,j,k,is) &
2787 *rra(i,j,k+1,1,is)+tt(i,j,k,is)*rxa(i,j,k+1,1,is))*denm
2788 rxa(i,j,k,1,is)= rs(i,j,k,is)+ts(i,j,k,is) &
2789 *rxa(i,j,k+1,1,is)*denm
2790 rra(i,j,k,2,is)=rra(i,j,k,1,is)
2791 rxa(i,j,k,2,is)=rxa(i,j,k,1,is)
2792 enddo
2793 enddo
2794 enddo
2795
2796 !-----for middle clouds
2797
2798 do 20 im=itm,2
2799
2800 do j= 1, n
2801 do i= 1, m
2802 do k= icb(i,j)-1,ict(i,j),-1
2803 denm=ts(i,j,k,im)/( 1.-rs(i,j,k,im)*rxa(i,j,k+1,im,is) )
2804 rra(i,j,k,im,is)= rr(i,j,k,im)+(td(i,j,k,im) &
2805 *rra(i,j,k+1,im,is)+tt(i,j,k,im)*rxa(i,j,k+1,im,is))*denm
2806 rxa(i,j,k,im,is)= rs(i,j,k,im)+ts(i,j,k,im) &
2807 *rxa(i,j,k+1,im,is)*denm
2808 enddo
2809 enddo
2810 enddo
2811
2812 20 continue
2813
2814 !-----integration over eight sky situations.
2815 ! ih, im, is denotes high, middle and low cloud groups.
2816
2817 do 100 ih=itm,2
2818
2819 !-----clear portion
2820
2821 if(ih.eq.1) then
2822 do j=1,n
2823 do i=1,m
2824 ch(i,j)=1.0-cc(i,j,1)
2825 enddo
2826 enddo
2827
2828 else
2829
2830 !-----cloudy portion
2831
2832 do j=1,n
2833 do i=1,m
2834 ch(i,j)=cc(i,j,1)
2835 enddo
2836 enddo
2837
2838 endif
2839
2840 do 100 im=itm,2
2841
2842 !-----clear portion
2843
2844 if(im.eq.1) then
2845
2846 do j=1,n
2847 do i=1,m
2848 cm(i,j)=ch(i,j)*(1.0-cc(i,j,2))
2849 enddo
2850 enddo
2851
2852 else
2853
2854 !-----cloudy portion
2855
2856 do j=1,n
2857 do i=1,m
2858 cm(i,j)=ch(i,j)*cc(i,j,2)
2859 enddo
2860 enddo
2861
2862 endif
2863
2864 do 100 is=itm,2
2865
2866 !-----clear portion
2867
2868 if(is.eq.1) then
2869
2870 do j=1,n
2871 do i=1,m
2872 ct(i,j)=cm(i,j)*(1.0-cc(i,j,3))
2873 enddo
2874 enddo
2875
2876 else
2877
2878 !-----cloudy portion
2879
2880 do j=1,n
2881 do i=1,m
2882 ct(i,j)=cm(i,j)*cc(i,j,3)
2883 enddo
2884 enddo
2885
2886 endif
2887
2888 !-----add one layer at a time, going down.
2889
2890 do j= 1, n
2891 do i= 1, m
2892 do k= icb(i,j), np
2893 denm = ts(i,j,k,is)/( 1.-rsa(i,j,k-1,ih,im)*rs(i,j,k,is) )
2894 tda(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*td(i,j,k,is)
2895 tta(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*tt(i,j,k,is) &
2896 +(tda(i,j,k-1,ih,im)*rr(i,j,k,is) &
2897 *rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm
2898 rsa(i,j,k,ih,im)= rs(i,j,k,is)+ts(i,j,k,is) &
2899 *rsa(i,j,k-1,ih,im)*denm
2900 enddo
2901 enddo
2902 enddo
2903
2904 !-----add one layer at a time, going up.
2905
2906 do j= 1, n
2907 do i= 1, m
2908 do k= ict(i,j)-1,1,-1
2909 denm =ts(i,j,k,ih)/(1.-rs(i,j,k,ih)*rxa(i,j,k+1,im,is))
2910 rra(i,j,k,im,is)= rr(i,j,k,ih)+(td(i,j,k,ih) &
2911 *rra(i,j,k+1,im,is)+tt(i,j,k,ih)*rxa(i,j,k+1,im,is))*denm
2912 rxa(i,j,k,im,is)= rs(i,j,k,ih)+ts(i,j,k,ih) &
2913 *rxa(i,j,k+1,im,is)*denm
2914 enddo
2915 enddo
2916 enddo
2917
2918 !-----compute fluxes following eq (5) of Chou (1992)
2919
2920 ! fdndir is the direct downward flux
2921 ! fdndif is the diffuse downward flux
2922 ! fupdif is the diffuse upward flux
2923
2924 do k=2,np+1
2925 do j=1, n
2926 do i=1, m
2927 denm= 1./(1.- rxa(i,j,k,im,is)*rsa(i,j,k-1,ih,im))
2928 fdndir(i,j)= tda(i,j,k-1,ih,im)
2929 xx = tda(i,j,k-1,ih,im)*rra(i,j,k,im,is)
2930 fdndif(i,j)= (xx*rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm
2931 fupdif= (xx+tta(i,j,k-1,ih,im)*rxa(i,j,k,im,is))*denm
2932 flxdn(i,j,k)=fdndir(i,j)+fdndif(i,j)-fupdif
2933 flxdnu(i,j,k)=-fupdif
2934 flxdnd(i,j,k)=fdndir(i,j)+fdndif(i,j)
2935 enddo
2936 enddo
2937 enddo
2938
2939 do j=1, n
2940 do i=1, m
2941 flxdn(i,j,1)=1.0-rra(i,j,1,im,is)
2942 flxdnu(i,j,1)=-rra(i,j,1,im,is)
2943 flxdnd(i,j,1)=1.0
2944 enddo
2945 enddo
2946
2947 !-----summation of fluxes over all (eight) sky situations.
2948
2949 do k=1,np+1
2950 do j=1,n
2951 do i=1,m
2952 if(ih.eq.1 .and. im.eq.1 .and. is.eq.1) then
2953 fclr(i,j,k)=flxdn(i,j,k)
2954 endif
2955 fall(i,j,k)=fall(i,j,k)+flxdn(i,j,k)*ct(i,j)
2956 fallu(i,j,k)=fallu(i,j,k)+flxdnu(i,j,k)*ct(i,j)
2957 falld(i,j,k)=falld(i,j,k)+flxdnd(i,j,k)*ct(i,j)
2958 enddo
2959 enddo
2960 enddo
2961
2962 do j=1,n
2963 do i=1,m
2964 fsdir(i,j)=fsdir(i,j)+fdndir(i,j)*ct(i,j)
2965 fsdif(i,j)=fsdif(i,j)+fdndif(i,j)*ct(i,j)
2966 enddo
2967 enddo
2968
2969 100 continue
2970
2971 end subroutine cldflx
2972
2973 !*****************************************************************
2974
2975 subroutine flxco2(m,n,np,swc,swh,csm,df)
2976
2977 !*****************************************************************
2978
2979 !-----compute the reduction of clear-sky downward solar flux
2980 ! due to co2 absorption.
2981
2982 implicit none
2983
2984 !-----input parameters
2985
2986 integer m,n,np
2987 real csm(m,n),swc(m,n,np+1),swh(m,n,np+1),cah(22,19)
2988
2989 !-----output (undated) parameter
2990
2991 real df(m,n,np+1)
2992
2993 !-----temporary array
2994
2995 integer i,j,k,ic,iw
2996 real xx,clog,wlog,dc,dw,x1,x2,y2
2997
2998 !********************************************************************
2999 !-----include co2 look-up table
3000
3001 data ((cah(i,j),i=1,22),j= 1, 5)/ &
3002 0.9923, 0.9922, 0.9921, 0.9920, 0.9916, 0.9910, 0.9899, 0.9882, &
3003 0.9856, 0.9818, 0.9761, 0.9678, 0.9558, 0.9395, 0.9188, 0.8945, &
3004 0.8675, 0.8376, 0.8029, 0.7621, 0.7154, 0.6647, 0.9876, 0.9876, &
3005 0.9875, 0.9873, 0.9870, 0.9864, 0.9854, 0.9837, 0.9811, 0.9773, &
3006 0.9718, 0.9636, 0.9518, 0.9358, 0.9153, 0.8913, 0.8647, 0.8350, &
3007 0.8005, 0.7599, 0.7133, 0.6627, 0.9808, 0.9807, 0.9806, 0.9805, &
3008 0.9802, 0.9796, 0.9786, 0.9769, 0.9744, 0.9707, 0.9653, 0.9573, &
3009 0.9459, 0.9302, 0.9102, 0.8866, 0.8604, 0.8311, 0.7969, 0.7565, &
3010 0.7101, 0.6596, 0.9708, 0.9708, 0.9707, 0.9705, 0.9702, 0.9697, &
3011 0.9687, 0.9671, 0.9647, 0.9612, 0.9560, 0.9483, 0.9372, 0.9221, &
3012 0.9027, 0.8798, 0.8542, 0.8253, 0.7916, 0.7515, 0.7054, 0.6551, &
3013 0.9568, 0.9568, 0.9567, 0.9565, 0.9562, 0.9557, 0.9548, 0.9533, &
3014 0.9510, 0.9477, 0.9428, 0.9355, 0.9250, 0.9106, 0.8921, 0.8700, &
3015 0.8452, 0.8171, 0.7839, 0.7443, 0.6986, 0.6486/
3016
3017 data ((cah(i,j),i=1,22),j= 6,10)/ &
3018 0.9377, 0.9377, 0.9376, 0.9375, 0.9372, 0.9367, 0.9359, 0.9345, &
3019 0.9324, 0.9294, 0.9248, 0.9181, 0.9083, 0.8948, 0.8774, 0.8565, &
3020 0.8328, 0.8055, 0.7731, 0.7342, 0.6890, 0.6395, 0.9126, 0.9126, &
3021 0.9125, 0.9124, 0.9121, 0.9117, 0.9110, 0.9098, 0.9079, 0.9052, &
3022 0.9012, 0.8951, 0.8862, 0.8739, 0.8579, 0.8385, 0.8161, 0.7900, &
3023 0.7585, 0.7205, 0.6760, 0.6270, 0.8809, 0.8809, 0.8808, 0.8807, &
3024 0.8805, 0.8802, 0.8796, 0.8786, 0.8770, 0.8747, 0.8712, 0.8659, &
3025 0.8582, 0.8473, 0.8329, 0.8153, 0.7945, 0.7697, 0.7394, 0.7024, &
3026 0.6588, 0.6105, 0.8427, 0.8427, 0.8427, 0.8426, 0.8424, 0.8422, &
3027 0.8417, 0.8409, 0.8397, 0.8378, 0.8350, 0.8306, 0.8241, 0.8148, &
3028 0.8023, 0.7866, 0.7676, 0.7444, 0.7154, 0.6796, 0.6370, 0.5897, &
3029 0.7990, 0.7990, 0.7990, 0.7989, 0.7988, 0.7987, 0.7983, 0.7978, &
3030 0.7969, 0.7955, 0.7933, 0.7899, 0.7846, 0.7769, 0.7664, 0.7528, &
3031 0.7357, 0.7141, 0.6866, 0.6520, 0.6108, 0.5646/
3032
3033 data ((cah(i,j),i=1,22),j=11,15)/ &
3034 0.7515, 0.7515, 0.7515, 0.7515, 0.7514, 0.7513, 0.7511, 0.7507, &
3035 0.7501, 0.7491, 0.7476, 0.7450, 0.7409, 0.7347, 0.7261, 0.7144, &
3036 0.6992, 0.6793, 0.6533, 0.6203, 0.5805, 0.5357, 0.7020, 0.7020, &
3037 0.7020, 0.7019, 0.7019, 0.7018, 0.7017, 0.7015, 0.7011, 0.7005, &
3038 0.6993, 0.6974, 0.6943, 0.6894, 0.6823, 0.6723, 0.6588, 0.6406, &
3039 0.6161, 0.5847, 0.5466, 0.5034, 0.6518, 0.6518, 0.6518, 0.6518, &
3040 0.6518, 0.6517, 0.6517, 0.6515, 0.6513, 0.6508, 0.6500, 0.6485, &
3041 0.6459, 0.6419, 0.6359, 0.6273, 0.6151, 0.5983, 0.5755, 0.5458, &
3042 0.5095, 0.4681, 0.6017, 0.6017, 0.6017, 0.6017, 0.6016, 0.6016, &
3043 0.6016, 0.6015, 0.6013, 0.6009, 0.6002, 0.5989, 0.5967, 0.5932, &
3044 0.5879, 0.5801, 0.5691, 0.5535, 0.5322, 0.5043, 0.4700, 0.4308, &
3045 0.5518, 0.5518, 0.5518, 0.5518, 0.5518, 0.5518, 0.5517, 0.5516, &
3046 0.5514, 0.5511, 0.5505, 0.5493, 0.5473, 0.5441, 0.5393, 0.5322, &
3047 0.5220, 0.5076, 0.4878, 0.4617, 0.4297, 0.3929/
3048
3049 data ((cah(i,j),i=1,22),j=16,19)/ &
3050 0.5031, 0.5031, 0.5031, 0.5031, 0.5031, 0.5030, 0.5030, 0.5029, &
3051 0.5028, 0.5025, 0.5019, 0.5008, 0.4990, 0.4960, 0.4916, 0.4850, &
3052 0.4757, 0.4624, 0.4441, 0.4201, 0.3904, 0.3564, 0.4565, 0.4565, &
3053 0.4565, 0.4564, 0.4564, 0.4564, 0.4564, 0.4563, 0.4562, 0.4559, &
3054 0.4553, 0.4544, 0.4527, 0.4500, 0.4460, 0.4400, 0.4315, 0.4194, &
3055 0.4028, 0.3809, 0.3538, 0.3227, 0.4122, 0.4122, 0.4122, 0.4122, &
3056 0.4122, 0.4122, 0.4122, 0.4121, 0.4120, 0.4117, 0.4112, 0.4104, &
3057 0.4089, 0.4065, 0.4029, 0.3976, 0.3900, 0.3792, 0.3643, 0.3447, &
3058 0.3203, 0.2923, 0.3696, 0.3696, 0.3696, 0.3696, 0.3696, 0.3696, &
3059 0.3695, 0.3695, 0.3694, 0.3691, 0.3687, 0.3680, 0.3667, 0.3647, &
3060 0.3615, 0.3570, 0.3504, 0.3409, 0.3279, 0.3106, 0.2892, 0.2642/
3061
3062 !********************************************************************
3063 !-----table look-up for the reduction of clear-sky solar
3064 ! radiation due to co2. The fraction 0.0343 is the
3065 ! extraterrestrial solar flux in the co2 bands.
3066
3067 do k= 2, np+1
3068 do j= 1, n
3069 do i= 1, m
3070 xx=1./.3
3071 clog=log10(swc(i,j,k)*csm(i,j))
3072 wlog=log10(swh(i,j,k)*csm(i,j))
3073 ic=int( (clog+3.15)*xx+1.)
3074 iw=int( (wlog+4.15)*xx+1.)
3075 if(ic.lt.2)ic=2
3076 if(iw.lt.2)iw=2
3077 if(ic.gt.22)ic=22
3078 if(iw.gt.19)iw=19
3079 dc=clog-float(ic-2)*.3+3.
3080 dw=wlog-float(iw-2)*.3+4.
3081 x1=cah(1,iw-1)+(cah(1,iw)-cah(1,iw-1))*xx*dw
3082 x2=cah(ic-1,iw-1)+(cah(ic-1,iw)-cah(ic-1,iw-1))*xx*dw
3083 y2=x2+(cah(ic,iw-1)-cah(ic-1,iw-1))*xx*dc
3084 if (x1.lt.y2) x1=y2
3085 df(i,j,k)=df(i,j,k)+0.0343*(x1-y2)
3086 enddo
3087 enddo
3088 enddo
3089
3090 end subroutine flxco2
3091
3092 !*****************************************************************
3093
3094 subroutine o3prof (np, pres, ozone, its, ite, kts, kte, p, o3)
3095
3096 !*****************************************************************
3097 implicit none
3098 !*****************************************************************
3099 !
3100 integer iprof,m,np,its,ite,kts,kte
3101 integer i,k,ko,kk
3102 real pres(np),ozone(np)
3103 real p(its:ite,kts:kte),o3(its:ite,kts:kte)
3104
3105 ! Statement function
3106
3107 real Linear, x1, y1, x2, y2, x
3108 Linear(x1, y1, x2, y2, x) = &
3109 (y1 * (x2 - x) + y2 * (x - x1)) / (x2 - x1)
3110 !
3111 do k = 1,np
3112 pres(k) = alog(pres(k))
3113 enddo
3114 do k = kts,kte
3115 do i = its, ite
3116 p(i,k) = alog(p(i,k))
3117 end do
3118 end do
3119
3120 ! assume the pressure at model top is greater than pres(1)
3121 ! if it is not, this part needs to change
3122
3123 do i = its, ite
3124 ko = 1
3125 do k = kts+1, kte
3126 do while (ko .lt. np .and. p(i,k) .gt. pres(ko))
3127 ko = ko + 1
3128 end do
3129 o3(i,k) = Linear (pres(ko), ozone(ko), &
3130 pres(ko-1), ozone(ko-1), &
3131 p(i,k))
3132 ko = ko - 1
3133 end do
3134 end do
3135
3136 ! calculate top lay O3
3137
3138 do i = its, ite
3139 ko = 1
3140 k = kts
3141 do while (ko .le. np .and. p(i,k) .gt. pres(ko))
3142 ko = ko + 1
3143 end do
3144 IF (ko-1 .le. 1) then
3145 O3(i,k)=ozone(k)
3146 ELSE
3147 O3(i,k)=0.
3148 do kk=ko-2,1,-1
3149 O3(i,k)=O3(i,k)+ozone(kk)*(pres(kk+1)-pres(kk))
3150 enddo
3151 O3(i,k)=O3(i,k)/(pres(ko-1)-pres(1))
3152 ENDIF
3153 ! print*,'O3=',i,k,ko,O3(i,k),p(i,k),ko,pres(ko),pres(ko-1)
3154 end do
3155
3156 end subroutine o3prof
3157
3158 !-----------------------------------------
3159 SUBROUTINE gsfc_swinit(cen_lat, allowed_to_read)
3160
3161 REAL, INTENT(IN ) :: cen_lat
3162 LOGICAL, INTENT(IN ) :: allowed_to_read
3163
3164 center_lat=cen_lat
3165
3166 END SUBROUTINE gsfc_swinit
3167
3168
3169 END MODULE module_ra_gsfcsw