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