Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_ra_gsfcsw.F @ 3568

Last change on this file since 3568 was 2759, checked in by aslmd, 3 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 110.3 KB

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

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

Download in other formats:

Original Format