Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_ra_gsfcsw.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 109.0 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format