Context Navigation

photolysis_online.F @ 3556

Last change on this file since 3556 was 3530, checked in by flefevre, 3 weeks ago

Implementation of temperature-dependent SO absorption cross-sections from Heays et al., 2023.

The SO cross-section file must be downloaded from the link below and put in the /cross_sections directory:

https://owncloud.latmos.ipsl.fr/index.php/s/3xLzwn7G587sTdB

File size: 66.9 KB

Rev	Line
[2780]	1	!==============================================================================
	2
	3	subroutine photolysis_online(nlayer, nb_phot_max,
[2974]	4	$ alt, press, temp, mmean,
[2851]	5	$ i_co2, i_co, i_o, i_o1d, i_o2, i_o3,i_h2,
[2780]	6	$ i_oh, i_ho2, i_h2o2, i_h2o,i_h,i_hcl,
[2925]	7	$ i_cl2, i_hocl, i_so2, i_so, i_so3, i_s2,
[2780]	8	$ i_clo, i_ocs, i_cocl2, i_h2so4, i_cl,
[2795]	9	$ i_no2, i_no, i_n2, i_n2d, nesp, rm,
[2780]	10	$ sza, dist_sol, v_phot)
	11
	12	use photolysis_mod
	13
	14	implicit none
[2836]	15
[2780]	16	! input
	17
	18	integer, intent(in) :: nesp ! total number of chemical species
	19	integer, intent(in) :: nlayer
	20	integer, intent(in) :: nb_phot_max
	21	integer, intent(in) :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3,
	22	$ i_oh, i_ho2, i_h2o2, i_h2o, i_h, i_hcl,
[2925]	23	$ i_cl2, i_hocl, i_so2, i_so, i_so3,i_s2,
	24	$ i_clo,i_ocs, i_cl, i_cocl2, i_h2so4,
	25	$ i_no2,i_no, i_n2, i_n2d, i_h2
[2780]	26
[2974]	27	real, dimension(nlayer), intent(in) :: press, temp, mmean ! pressure (hpa)/temperature (k)/mean molecular mass (g.mol-1)
[2780]	28	real, dimension(nlayer), intent(in) :: alt ! altitude (km)
	29	real, dimension(nlayer,nesp), intent(in) :: rm ! mixing ratios
	30	real :: tau ! integrated aerosol optical depth at the surface
	31	real, intent(in) :: sza ! solar zenith angle (degrees)
	32	real, intent(in) :: dist_sol ! solar distance (au)
	33
	34	! output
	35
	36	real (kind = 8), dimension(nlayer,nb_phot_max) :: v_phot ! photolysis rates (s-1)
	37
	38	! solar flux at venus
	39
	40	real, dimension(nw) :: fvenus ! solar flux (w.m-2.nm-1)
	41	real :: factor
	42
	43	! cross-sections
	44
	45	real, dimension(nlayer,nw,nphot) :: sj ! general cross-section array (cm2)
	46
	47	! atmosphere
	48
[2974]	49	real, dimension(nlayer+1) :: zpress, zalt, ztemp, zmmean ! pressure (hpa)/altitude (km)/temperature (k)/mean molecular mass (g.mol-1)
[2780]	50
[2974]	51	real, dimension(nlayer+1) :: colinc ! air column increment (molecule.cm-2)
	52	real, dimension(nlayer+1,nw) :: dtrl ! rayleigh optical depth
	53	real, dimension(nlayer+1,nw) :: dtaer ! aerosol optical depth
	54	real, dimension(nlayer+1,nw) :: omaer ! aerosol single scattering albedo
	55	real, dimension(nlayer+1,nw) :: gaer ! aerosol asymmetry parameter
	56	real, dimension(nlayer+1,nw) :: dtcld ! cloud optical depth
	57	real, dimension(nlayer+1,nw) :: omcld ! cloud single scattering albedo
	58	real, dimension(nlayer+1,nw) :: gcld ! cloud asymmetry parameter
	59	real, dimension(nlayer+1,nw,nabs) :: dtgas ! optical depth for each gas
	60	real, dimension(nlayer+1,nw) :: dagas ! total gas optical depth
	61	real, dimension(nlayer+1) :: edir, edn, eup ! normalised irradiances
	62	real, dimension(nlayer+1) :: fdir, fdn, fup ! normalised actinic fluxes
	63	real, dimension(nlayer+1) :: saflux ! total actinic flux
[2780]	64
[2974]	65	integer, dimension(0:nlayer+1) :: nid
	66	real, dimension(0:nlayer+1,nlayer+1) :: dsdh
	67
[2780]	68	integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d, j_o3_o,
	69	$ j_h2o, j_h2o2, j_ho2, j_h, j_hcl, j_cl2, j_hocl, j_so2,
[2925]	70	$ j_so, j_so3,j_s2, j_clo, j_ocs, j_cocl2, j_h2so4,
	71	$ j_no2, j_no, j_n2, j_h2
[2780]	72
	73	integer :: a_o2, a_co2, a_o3, a_h2o, a_h2o2, a_ho2, a_hcl, a_cl2,
[2925]	74	$ a_hocl, a_so2, a_so, a_so3, a_s2, a_clo, a_ocs,
	75	$ a_cocl2, a_h2so4, a_no2, a_no, a_n2, a_h2
[2974]	76	integer :: nlev, i, ilay, ilev, iw, ialt
[2780]	77	real :: deltaj
[2795]	78	logical :: deutchem
[2780]	79
	80	! absorbing gas numbering
	81
	82	a_o2 = 1 ! o2
	83	a_co2 = 2 ! co2
	84	a_o3 = 3 ! o3
[2851]	85	a_h2 = 4 ! h2
	86	a_h2o = 5 ! h2o
	87	a_h2o2 = 6 ! h2o2
	88	a_ho2 = 7 ! ho2
	89	a_hcl = 8 ! hcl
	90	a_cl2 = 9 ! cl2
	91	a_hocl = 10 ! hocl
	92	a_so2 = 11 ! so2
	93	a_so = 12 ! so
	94	a_so3 = 13 ! so3
[2925]	95	a_s2 = 14 ! s2
	96	a_clo = 15 ! clo
	97	a_ocs = 16 ! ocs
	98	a_cocl2 = 17 ! cocl2
	99	a_h2so4 = 18 ! h2so4
	100	a_no2 = 19 ! no2
	101	a_no = 20 ! no
	102	a_n2 = 21 ! n2
[2780]	103
	104	! photodissociation rates numbering.
	105	! photodissociations must be ordered the same way in subroutine "indice"
	106
	107	j_o2_o = 1 ! o2 + hv -> o + o
	108	j_o2_o1d = 2 ! o2 + hv -> o + o(1d)
	109	j_co2_o = 3 ! co2 + hv -> co + o
	110	j_co2_o1d = 4 ! co2 + hv -> co + o(1d)
	111	j_o3_o1d = 5 ! o3 + hv -> o2 + o(1d)
	112	j_o3_o = 6 ! o3 + hv -> o2 + o
[2851]	113	j_h2 = 7 ! h2 + hv -> h + h
	114	j_h2o = 8 ! h2o + hv -> h + oh
	115	j_ho2 = 9 ! ho2 + hv -> oh + o
	116	j_h2o2 = 10 ! h2o2 + hv -> oh + oh
	117	j_hcl = 11 ! hcl + hv -> h + cl
	118	j_cl2 = 12 ! cl2 + hv -> cl + cl
	119	j_hocl = 13 ! hocl + hv -> oh + cl
	120	j_so2 = 14 ! so2 + hv -> so + o
	121	j_so = 15 ! so + hv -> s + o
	122	j_so3 = 16 ! so3 + hv -> so2 + o
[2925]	123	j_s2 = 17 ! s2 + hv -> s + s
	124	j_clo = 18 ! clo + hv -> cl + o
	125	j_ocs = 19 ! ocs + hv -> co + s
	126	j_cocl2 = 20 ! cocl2 + hv -> 2cl + co
	127	j_h2so4 = 21 ! h2so4 + hv -> so3 + h2o
	128	j_no2 = 22 ! no2 + hv -> no + o
	129	j_no = 23 ! no + hv -> n + o
	130	j_n2 = 24 ! n2 + hv -> n(2d) + n
[2780]	131
[2795]	132	! j_hdo_od = ! hdo + hv -> od + h
	133	! j_hdo_d = ! hdo + hv -> d + oh
	134
[2974]	135	!==== define working vertical grid for the uv radiative code
	136
	137	nlev = nlayer + 1
	138
	139	do ilev = 1,nlev-1
	140	zpress(ilev) = press(ilev)
	141	zalt(ilev) = alt(ilev)
	142	ztemp(ilev) = temp(ilev)
	143	zmmean(ilev) = mmean(ilev)
	144	end do
	145
	146	zpress(nlev) = 0. ! top of the atmosphere
	147	zalt(nlev) = zalt(nlev-1) + (zalt(nlev-1) - zalt(nlev-2))
	148	ztemp(nlev) = ztemp(nlev-1)
	149	zmmean(nlev) = zmmean(nlev-1)
	150
[2780]	151	!==== air column increments and rayleigh optical depth
	152
[2974]	153	call setair(nlev, nw, wl, wc, zpress, ztemp, zmmean, colinc, dtrl)
[2780]	154
	155	!==== set temperature-dependent cross-sections and optical depths
	156
[2974]	157	dtgas(:,:,:) = 0.
	158
[2780]	159	! o2:
	160
	161	call seto2(nphot, nlayer, nw, wc, mopt, temp, xso2_150, xso2_200,
	162	$ xso2_250, xso2_300, yieldo2, j_o2_o, j_o2_o1d,
[2974]	163	$ colinc(1:nlayer), rm(:,i_o2),
	164	$ dtgas(1:nlayer,:,a_o2), sj)
[2780]	165
	166	! co2:
	167
	168	call setco2(nphot, nlayer, nw, wc, temp, xsco2_195, xsco2_295,
	169	$ xsco2_370, yieldco2, j_co2_o, j_co2_o1d,
[2974]	170	$ colinc(1:nlayer), rm(:,i_co2),
	171	$ dtgas(1:nlayer,:,a_co2), sj)
[2780]	172
	173	! o3:
	174
	175	call seto3(nphot, nlayer, nw, wc, temp, xso3_218, xso3_298,
[2974]	176	$ j_o3_o, j_o3_o1d, colinc(1:nlayer), rm(:,i_o3),
	177	$ dtgas(1:nlayer,:,a_o3), sj)
[2780]	178
	179	! h2o2:
	180
	181	call seth2o2(nphot, nlayer, nw, wc, temp, xsh2o2, j_h2o2,
[2974]	182	$ colinc(1:nlayer), rm(:,i_h2o2),
	183	$ dtgas(1:nlayer,:,a_h2o2), sj)
[2780]	184
[3530]	185	! so:
	186
	187	call setso(nphot, nlayer, nw, wc, temp, xsso_150, xsso_250,
	188	$ j_so, colinc(1:nlayer), rm(:,i_so),
	189	$ dtgas(1:nlayer,:,a_so), sj)
	190
[2780]	191	! so2:
	192
[3530]	193	call setso2(nphot, nlayer, nw, wc, temp, xsso2_200, xsso2_298,
	194	$ xsso2_360, j_so2, colinc(1:nlayer), rm(:,i_so2),
[2974]	195	$ dtgas(1:nlayer,:,a_so2), sj)
[2780]	196
	197	! no2:
	198
[2795]	199	call setno2(nphot, nlayer, nw, wc, temp, xsno2, xsno2_220,
	200	$ xsno2_294, yldno2_248, yldno2_298, j_no2,
[2974]	201	$ colinc(1:nlayer), rm(:,i_no2),
	202	$ dtgas(1:nlayer,:,a_no2), sj)
[2780]	203
	204	!==== temperature independent optical depths and cross-sections
	205
	206	! optical depths
	207
	208	do ilay = 1,nlayer
	209	do iw = 1,nw-1
[2851]	210	dtgas(ilay,iw,a_h2) = colinc(ilay)rm(ilay,i_h2)xsh2(iw)
[2780]	211	dtgas(ilay,iw,a_h2o) = colinc(ilay)rm(ilay,i_h2o)xsh2o(iw)
	212	dtgas(ilay,iw,a_ho2) = colinc(ilay)rm(ilay,i_ho2)xsho2(iw)
	213	dtgas(ilay,iw,a_hcl) = colinc(ilay)rm(ilay,i_hcl)xshcl(iw)
	214	dtgas(ilay,iw,a_cl2) = colinc(ilay)rm(ilay,i_cl2)xscl2(iw)
	215	dtgas(ilay,iw,a_hocl) = colinc(ilay)*rm(ilay,i_hocl)
	216	$ *xshocl(iw)
	217	dtgas(ilay,iw,a_so3) = colinc(ilay)rm(ilay,i_so3)xsso3(iw)
[2925]	218	dtgas(ilay,iw,a_s2) = colinc(ilay)rm(ilay,i_s2)xss2(iw)
[2780]	219	dtgas(ilay,iw,a_clo) = colinc(ilay)rm(ilay,i_clo)xsclo(iw)
	220	dtgas(ilay,iw,a_ocs) = colinc(ilay)rm(ilay,i_ocs)xsocs(iw)
	221	dtgas(ilay,iw,a_cocl2) = colinc(ilay)*rm(ilay,i_cocl2)
	222	$ *xscocl2(iw)
	223	dtgas(ilay,iw,a_h2so4) = colinc(ilay)*rm(ilay,i_h2so4)
	224	$ *xsh2so4(iw)
[2795]	225	dtgas(ilay,iw,a_no) = colinc(ilay)rm(ilay,i_no)xsno(iw)
	226	dtgas(ilay,iw,a_n2) = colinc(ilay)rm(ilay,i_n2)xsn2(iw)
[2780]	227	end do
	228	end do
	229
	230	! total gas optical depth
	231
	232	dagas(:,:) = 0.
	233
	234	do ilay = 1,nlayer
	235	do iw = 1,nw-1
	236	do i = 1,nabs
	237	dagas(ilay,iw) = dagas(ilay,iw) + dtgas(ilay,iw,i)
	238	end do
	239	end do
	240	end do
	241
	242	! cross-sections
	243
	244	do ilay = 1,nlayer
	245	do iw = 1,nw-1
[2851]	246	sj(ilay,iw,j_h2) = xsh2(iw) ! h2
[2780]	247	sj(ilay,iw,j_h2o) = xsh2o(iw) ! h2o
	248	sj(ilay,iw,j_ho2) = xsho2(iw) ! ho2
	249	sj(ilay,iw,j_hcl) = xshcl(iw) ! hcl
	250	sj(ilay,iw,j_cl2) = xscl2(iw) ! cl2
	251	sj(ilay,iw,j_hocl) = xshocl(iw) ! hocl
[2925]	252	sj(ilay,iw,j_s2) = xss2(iw) ! s2
[2780]	253	sj(ilay,iw,j_so3) = xsso3(iw) ! so3
	254	sj(ilay,iw,j_clo) = xsclo(iw) ! clo
	255	sj(ilay,iw,j_ocs) = xsocs(iw) ! ocs
	256	sj(ilay,iw,j_cocl2) = xscocl2(iw) ! cocl2
	257	sj(ilay,iw,j_h2so4) = xsh2so4(iw) ! h2so4
[2795]	258	sj(ilay,iw,j_no) = xsno(iw)*yieldno(iw) ! no
	259	sj(ilay,iw,j_n2) = xsn2(iw)*yieldn2(iw) ! n2
[2780]	260	end do
	261	end do
	262
[2795]	263	! hdo cross section
	264
[2780]	265	! if (deutchem) then
[2974]	266	! do ilay = 1,nlayer
	267	! do iw = 1,nw-1
	268	! !Two chanels for hdo, with same cross section
[2780]	269	! sj(ilay,iw,j_hdo_od) = 0.5*xshdo(iw) ! hdo -> od + h
	270	! sj(ilay,iw,j_hdo_d) = 0.5*xshdo(iw) ! hdo -> d + oh
[2974]	271	! end do
	272	! end do
	273	! end if
[2780]	274
[2795]	275	!==== set aerosol properties and optical depth
[2780]	276
[2795]	277	tau = 0. ! no solid aerosols for the present time
[2780]	278
[2974]	279	call setaer(nlev,zalt,tau,nw,dtaer,omaer,gaer)
[2780]	280
	281	!==== set cloud properties and optical depth
	282
[2974]	283	call setcld(nlev,zalt,nw,wl,dtcld,omcld,gcld)
[2780]	284
	285	!==== slant path lengths in spherical geometry
	286
[2974]	287	call sphers(nlev,zalt,sza,dsdh,nid)
[2780]	288
	289	!==== solar flux at venus
	290
	291	factor = (1./dist_sol)**2.
	292	do iw = 1,nw-1
	293	fvenus(iw) = f(iw)*factor
	294	end do
	295
	296	!==== initialise photolysis rates
	297
	298	v_phot(:,1:nphot) = 0.
	299
	300	!==== start of wavelength lopp
	301
	302	do iw = 1,nw-1
	303
	304	! monochromatic radiative transfer. outputs are:
	305	! normalized irradiances edir(nlayer), edn(nlayer), eup(nlayer)
	306	! normalized actinic fluxes fdir(nlayer), fdn(nlayer), fup(nlayer)
	307	! where
	308	! dir = direct beam, dn = down-welling diffuse, up = up-welling diffuse
	309
[2974]	310	call rtlink(nlev, nw, iw, albedo(iw), sza, dsdh, nid, dtrl,
[2780]	311	$ dagas, dtcld, omcld, gcld, dtaer, omaer, gaer,
	312	$ edir, edn, eup, fdir, fdn, fup)
	313
	314
	315	! spherical actinic flux
	316
	317	do ilay = 1,nlayer
	318	saflux(ilay) = fvenus(iw)*(fdir(ilay) + fdn(ilay) + fup(ilay))
	319	end do
	320
	321	! photolysis rate integration
	322
	323	do i = 1,nphot
	324	do ilay = 1,nlayer
	325	deltaj = saflux(ilay)*sj(ilay,iw,i)
	326	v_phot(ilay,i) = v_phot(ilay,i) + deltaj*(wu(iw)-wl(iw))
	327	end do
	328	end do
	329
	330	! eliminate small values
	331
	332	where (v_phot(:,1:nphot) < 1.e-30)
	333	v_phot(:,1:nphot) = 0.
	334	end where
[2795]	335
[2780]	336	end do ! iw
	337
	338	contains
	339
	340	!==============================================================================
	341
	342	subroutine setair(nlev, nw, wl, wc, press, temp, zmmean,
	343	$ colinc, dtrl)
	344
	345	-----------------------------------------------------------------------------
	346	= PURPOSE: =
	347	= computes air column increments and rayleigh optical depth =
	348	-----------------------------------------------------------------------------
	349
	350	implicit none
	351
	352	! input:
	353
	354	integer :: nlev, nw
	355
	356	real, dimension(nw) :: wl, wc ! lower and central wavelength grid (nm)
	357	real, dimension(nlev) :: press, temp, zmmean ! pressure (hpa), temperature (k), molecular mass (g.mol-1)
	358
	359	! output:
	360
	361	real, dimension(nlev) :: colinc ! air column increments (molecule.cm-2)
	362	real, dimension(nlev,nw) :: dtrl ! rayleigh optical depth
	363
	364	! local:
	365
	366	real, parameter :: avo = 6.022e23
	367	real, parameter :: g = 8.87
	368	real :: dp, nu
	369	real, dimension(nw) :: srayl
	370	integer :: ilev, iw
	371
	372	! compute column increments
	373
	374	do ilev = 1, nlev-1
	375	dp = (press(ilev) - press(ilev+1))*100.
	376	colinc(ilev) = avo0.1dp/(zmmean(ilev)*g)
	377	end do
[2974]	378	colinc(nlev) = 0.
[2780]	379
	380	do iw = 1, nw - 1
	381
	382	! co2 rayleigh cross-section
	383	! ityaksov et al., chem. phys. lett., 462, 31-34, 2008
	384
	385	nu = 1./(wc(iw)*1.e-7)
	386	srayl(iw) = 1.78e-26nu*(4. + 0.625)
	387	srayl(iw) = srayl(iw)*1.e-20 ! cm2
	388
	389	do ilev = 1, nlev
	390	dtrl(ilev,iw) = colinc(ilev)*srayl(iw) ! cm2
	391	end do
	392	end do
	393
	394	end subroutine setair
	395
	396	!==============================================================================
	397
	398	subroutine setco2(nd, nlayer, nw, wc, tlay, xsco2_195, xsco2_295,
	399	$ xsco2_370, yieldco2, j_co2_o, j_co2_o1d,
	400	$ colinc, rm, dt, sj)
	401
	402	!-----------------------------------------------------------------------------*
	403	!= PURPOSE: =*
	404	!= Set up the CO2 temperature-dependent cross-sections and optical depth =*
	405	!-----------------------------------------------------------------------------*
	406
	407	implicit none
	408
	409	! input:
	410
	411	integer :: nd ! number of photolysis rates
	412	integer :: nlayer ! number of vertical layers
	413	integer :: nw ! number of wavelength grid points
	414	integer :: j_co2_o, j_co2_o1d ! photolysis indexes
	415	real, dimension(nw) :: wc ! central wavelength for each interval
	416	real, dimension(nw) :: xsco2_195, xsco2_295, xsco2_370 ! co2 cross-sections (cm2)
	417	real, dimension(nw) :: yieldco2 ! co2 photodissociation yield
	418	real, dimension(nlayer) :: tlay ! temperature (k)
	419	real, dimension(nlayer) :: rm ! co2 mixing ratio
	420	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
	421
	422	! output:
	423
	424	real, dimension(nlayer,nw) :: dt ! optical depth
	425	real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2)
	426
	427	! local:
	428
	429	integer :: extrapol
	430	integer :: i, l
	431	real :: temp, sco2
	432
	433	! extrapol = 0 no extrapolation below 195 k
	434	! extrapol = 1 extrapolation below 195 k
	435
	436	extrapol = 0
	437
	438	do i = 1, nlayer
	439	if (extrapol == 1) then
	440	temp = tlay(i)
	441	else
	442	temp = max(tlay(i), 195.)
	443	end if
	444	temp = min(temp, 370.)
	445	do l = 1, nw-1
	446	if (temp <= 295.) then
	447	if (xsco2_195(l) /= 0. .and. xsco2_295(l) /= 0.) then
	448	sco2 = alog(xsco2_195(l))
	449	$ + (alog(xsco2_295(l)) - alog(xsco2_195(l)))
	450	$ /(295. - 195.)*(temp - 195.)
	451	sco2 = exp(sco2)
	452	else
	453	sco2 = 0.
	454	end if
	455	else
	456	if (xsco2_295(l) /= 0. .and. xsco2_370(l) /= 0.) then
	457	sco2 = alog(xsco2_295(l))
	458	$ + (alog(xsco2_370(l)) - alog(xsco2_295(l)))
	459	$ /(370. - 295.)*(temp - 295.)
	460	sco2 = exp(sco2)
	461	else
	462	sco2 = 0.
	463	end if
	464	end if
	465
	466	! optical depth
	467
	468	dt(i,l) = colinc(i)rm(i)sco2
	469
	470	! production of o(1d) for wavelengths shorter than 167 nm
	471
	472	if (wc(l) >= 167.) then
	473	sj(i,l,j_co2_o) = sco2*yieldco2(l)
	474	sj(i,l,j_co2_o1d) = 0.
	475	else
	476	sj(i,l,j_co2_o) = 0.
	477	sj(i,l,j_co2_o1d) = sco2*yieldco2(l)
	478	end if
	479	end do
	480	end do
	481
	482	end subroutine setco2
	483
	484	!==============================================================================
	485
	486	subroutine seto2(nd, nlayer, nw, wc, mopt, tlay, xso2_150,
	487	$ xso2_200, xso2_250, xso2_300, yieldo2, j_o2_o,
	488	$ j_o2_o1d, colinc, rm, dt, sj)
	489
	490	!-----------------------------------------------------------------------------*
	491	!= PURPOSE: =*
	492	!= Set up the O2 temperature-dependent cross-sections and optical depth =*
	493	!-----------------------------------------------------------------------------*
	494
	495	implicit none
	496
	497	! input:
	498
	499	integer :: nd ! number of photolysis rates
	500	integer :: nlayer ! number of vertical layers
	501	integer :: nw ! number of wavelength grid points
	502	integer :: mopt ! high-res/low-res switch
	503	integer :: j_o2_o, j_o2_o1d ! photolysis indexes
	504	real, dimension(nw) :: wc ! central wavelength for each interval
	505	real, dimension(nw) :: xso2_150, xso2_200, xso2_250, ! o2 cross-sections (cm2)
	506	$ xso2_300
	507	real, dimension(nw) :: yieldo2 ! o2 photodissociation yield
	508	real, dimension(nlayer) :: tlay ! temperature (k)
	509	real, dimension(nlayer) :: rm ! o2 mixing ratio
	510	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
	511
	512	! output:
	513
	514	real, dimension(nlayer,nw) :: dt ! optical depth
	515	real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2)
	516
	517	! local:
	518
	519	integer :: ilev, iw
	520	real :: temp
	521	real :: xso2, factor
	522
	523	! correction by factor if low-resolution in schumann-runge bands
	524
	525	if (mopt == 1) then
	526	factor = 1.
	527	else if (mopt == 2 .or. mopt == 3) then
	528	factor = 0.8
	529	end if
	530
	531	! calculate temperature dependance
	532
	533	do ilev = 1,nlayer
	534	temp = max(tlay(ilev),150.)
	535	temp = min(temp, 300.)
	536	do iw = 1, nw-1
	537	if (tlay(ilev) > 250.) then
	538	xso2 = xso2_250(iw) + (xso2_300(iw) - xso2_250(iw))
	539	$ /(300. - 250.)*(temp - 250.)
	540	else if (tlay(ilev) > 200.) then
	541	xso2 = xso2_200(iw) + (xso2_250(iw) - xso2_200(iw))
	542	$ /(250. - 200.)*(temp - 200.)
	543	else
	544	xso2 = xso2_150(iw) + (xso2_200(iw) - xso2_150(iw))
	545	$ /(200. - 150.)*(temp - 150.)
	546	end if
	547
	548	if (wc(iw) > 180. .and. wc(iw) < 200.) then
	549	xso2 = xso2*factor
	550	end if
	551
	552	! optical depth
	553
	554	dt(ilev,iw) = colinc(ilev)rm(ilev)xso2
	555
	556	! production of o(1d) for wavelengths shorter than 175 nm
	557
	558	if (wc(iw) >= 175.) then
	559	sj(ilev,iw,j_o2_o) = xso2*yieldo2(iw)
	560	sj(ilev,iw,j_o2_o1d) = 0.
	561	else
	562	sj(ilev,iw,j_o2_o) = 0.
	563	sj(ilev,iw,j_o2_o1d) = xso2*yieldo2(iw)
	564	end if
	565
	566	end do
	567	end do
	568
	569	end subroutine seto2
	570
	571	!==============================================================================
	572
	573	subroutine seto3(nd, nlayer, nw, wc, tlay, xso3_218, xso3_298,
	574	$ j_o3_o, j_o3_o1d,
	575	$ colinc, rm, dt, sj)
	576
	577	!-----------------------------------------------------------------------------*
	578	!= PURPOSE: =*
	579	!= Set up the O3 temperature dependent cross-sections and optical depth =*
	580	!-----------------------------------------------------------------------------*
	581
	582	implicit none
	583
	584	! input:
	585
	586	integer :: nd ! number of photolysis rates
	587	integer :: nlayer ! number of vertical layers
	588	integer :: nw ! number of wavelength grid points
	589	integer :: j_o3_o, j_o3_o1d ! photolysis indexes
	590	real, dimension(nw) :: wc ! central wavelength for each interval
	591	real, dimension(nw) :: xso3_218, xso3_298 ! o3 cross-sections (cm2)
	592	real, dimension(nlayer) :: tlay ! temperature (k)
	593	real, dimension(nlayer) :: rm ! o3 mixing ratio
	594	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
	595
	596	! output:
	597
	598	real, dimension(nlayer,nw) :: dt ! optical depth
	599	real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2)
	600
	601	! local:
	602	!
	603	integer :: ilev, iw
	604	real :: temp
	605	real, dimension(nw) :: xso3(nw)
	606	real, dimension(nw) :: qy1d ! quantum yield for o(1d) production
	607	real :: q1, q2, a1, a2, a3
	608
	609	do ilev = 1, nlayer
	610	temp = max(tlay(ilev), 218.)
	611	temp = min(temp,298.)
	612	do iw = 1, nw-1
	613	xso3(iw) = xso3_218(iw) + (xso3_298(iw) - xso3_218(iw))
	614	$ /(298. - 218.) *(temp - 218.)
	615
	616	! optical depth
	617
	618	dt(ilev,iw) = colinc(ilev)rm(ilev)xso3(iw)
	619
	620	end do
	621
	622	! calculate quantum yield for o(1d) production (jpl 2006)
	623
	624	temp = max(tlay(ilev),200.)
	625	temp = min(temp,320.)
	626	do iw = 1, nw-1
	627	if (wc(iw) <= 306.) then
	628	qy1d(iw) = 0.90
	629	else if (wc(iw) > 306. .and. wc(iw) < 328.) then
	630	q1 = 1.
	631	q2 = exp(-825.518/(0.695*temp))
	632	a1 = (304.225 - wc(iw))/5.576
	633	a2 = (314.957 - wc(iw))/6.601
	634	a3 = (310.737 - wc(iw))/2.187
	635	qy1d(iw) = (q1/(q1 + q2))0.8036exp(-(a1a1a1*a1))
	636	$ + (q2/(q1 + q2))8.9061(temp/300.)**2.
	637	$ exp(-(a2a2))
	638	$ + 0.1192(temp/300.)1.5exp(-(a3*a3))
	639	$ + 0.0765
	640	else if (wc(iw) >= 328. .and. wc(iw) <= 340.) then
	641	qy1d(iw) = 0.08
	642	else
	643	qy1d(iw) = 0.
	644	endif
	645	end do
	646	do iw = 1, nw-1
	647	sj(ilev,iw,j_o3_o) = xso3(iw)*(1. - qy1d(iw))
	648	sj(ilev,iw,j_o3_o1d) = xso3(iw)*qy1d(iw)
	649	end do
	650	end do
	651
	652	end subroutine seto3
	653	!==============================================================================
	654
[3530]	655	subroutine setso(nd, nlayer, nw, wc, tlay, xsso_150, xsso_250,
	656	$ j_so, colinc, rm, dt, sj)
	657
	658	!-----------------------------------------------------------------------------*
	659	!= PURPOSE: =*
	660	!= Set up the so temperature dependent cross-sections and optical depth =*
	661	!-----------------------------------------------------------------------------*
	662
	663	implicit none
	664
	665	! input:
	666
	667	integer :: nd ! number of photolysis rates
	668	integer :: nlayer ! number of vertical layers
	669	integer :: nw ! number of wavelength grid points
	670	integer :: j_so ! photolysis indexe
	671	real, dimension(nw) :: wc ! central wavelength for each interval
	672	real, dimension(nw) :: xsso_150, xsso_250 ! so cross-sections (cm2)
	673	real, dimension(nlayer) :: tlay ! temperature (k)
	674	real, dimension(nlayer) :: rm ! so mixing ratio
	675	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
	676
	677	! output:
	678
	679	real, dimension(nlayer,nw) :: dt ! optical depth
	680	real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2)
	681
	682	! local:
	683	!
	684	integer :: ilev, iw
	685	real :: temp, xsso
	686
	687	! calculate temperature dependance
	688
	689	do ilev = 1,nlayer
	690	temp = max(tlay(ilev),150.)
	691	temp = min(temp, 250.)
	692	do iw = 1, nw-1
	693	if (xsso_150(iw) /= 0. .and. xsso_250(iw) /= 0.) then
	694	xsso = log(xsso_150(iw))
	695	$ + (log(xsso_250(iw)) - log(xsso_150(iw)))
	696	$ /(250. - 150.)*(temp - 150.)
	697
	698	sj(ilev,iw,j_so) = exp(xsso)
	699	else
	700	sj(ilev,iw,j_so) = 0.
	701	end if
	702
	703	! optical depth
	704
	705	dt(ilev,iw) = colinc(ilev)rm(ilev)sj(ilev,iw,j_so)
	706
	707	end do
	708	end do
	709
	710	end subroutine setso
	711
	712	!==============================================================================
	713
[2780]	714	subroutine setso2(nd, nlayer, nw, wc, tlay, xsso2_200, xsso2_298,
	715	$ xsso2_360, j_so2, colinc, rm, dt, sj)
	716
	717	!-----------------------------------------------------------------------------*
	718	!= PURPOSE: =*
	719	!= Set up the so2 temperature dependent cross-sections and optical depth =*
	720	!-----------------------------------------------------------------------------*
	721
	722	implicit none
	723
	724	! input:
	725
	726	integer :: nd ! number of photolysis rates
	727	integer :: nlayer ! number of vertical layers
	728	integer :: nw ! number of wavelength grid points
	729	integer :: j_so2 ! photolysis indexe
	730	real, dimension(nw) :: wc ! central wavelength for each interval
	731	real, dimension(nw) :: xsso2_200, xsso2_298, xsso2_360 ! so2 cross-sections (cm2)
	732	real, dimension(nlayer) :: tlay ! temperature (k)
	733	real, dimension(nlayer) :: rm ! so2 mixing ratio
	734	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
	735
	736	! output:
	737
	738	real, dimension(nlayer,nw) :: dt ! optical depth
	739	real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2)
	740
	741	! local:
	742	!
	743	integer :: ilev, iw
	744	real :: temp , xsso2
	745
	746
	747	! calculate temperature dependance
	748	do ilev = 1,nlayer
	749	temp = max(tlay(ilev),200.)
	750	temp = min(temp, 360.)
	751	do iw = 1, nw-1
	752	if (tlay(ilev) < 298.) then
	753	xsso2 = xsso2_200(iw) + (xsso2_298(iw) - xsso2_200(iw))
	754	$ /(298. - 200.)*(temp - 200.)
	755	else
	756	xsso2 = xsso2_298(iw) + (xsso2_360(iw) - xsso2_298(iw))
	757	$ /(360. - 298.)*(temp - 298.)
	758	end if
	759	! 219 nm photolysis treshold
	760	if (wc(iw) <= 219.) then
	761	sj(ilev,iw,j_so2) = xsso2
	762	else
	763	sj(ilev,iw,j_so2) = 0.
	764	end if
	765
	766	! optical depth
	767
	768	dt(ilev,iw) = colinc(ilev)rm(ilev)xsso2
	769
	770	end do
	771	end do
	772
	773	end subroutine setso2
	774
	775	!==============================================================================
	776
	777	subroutine seth2o2(nd, nlayer, nw, wc, tlay, xsh2o2, j_h2o2,
	778	$ colinc, rm, dt, sj)
	779
	780	!-----------------------------------------------------------------------------*
	781	!= PURPOSE: =*
	782	!= Set up the h2o2 temperature dependent cross-sections and optical depth =*
	783	!-----------------------------------------------------------------------------*
	784
	785	implicit none
	786
	787	! input:
	788
	789	integer :: nd ! number of photolysis rates
	790	integer :: nlayer ! number of vertical layers
	791	integer :: nw ! number of wavelength grid points
	792	integer :: j_h2o2 ! photolysis index
	793	real, dimension(nw) :: wc ! central wavelength for each interval
	794	real, dimension(nw) :: xsh2o2 ! h2o2 cross-sections (cm2)
	795	real, dimension(nlayer) :: tlay ! temperature (k)
	796	real, dimension(nlayer) :: rm ! h2o2 mixing ratio
	797	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
	798
	799	! output:
	800
	801	real, dimension(nlayer,nw) :: dt ! optical depth
	802	real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2)
	803
	804	! local:
	805
	806	integer :: ilev, iw
	807	real :: a0, a1, a2, a3, a4, a5, a6, a7
	808	real :: b0, b1, b2, b3, b4
	809	real :: lambda, suma, sumb, chi, temp, xs
	810
	811	A0 = 6.4761E+04
	812	A1 = -9.2170972E+02
	813	A2 = 4.535649
	814	A3 = -4.4589016E-03
	815	A4 = -4.035101E-05
	816	A5 = 1.6878206E-07
	817	A6 = -2.652014E-10
	818	A7 = 1.5534675E-13
	819
	820	B0 = 6.8123E+03
	821	B1 = -5.1351E+01
	822	B2 = 1.1522E-01
	823	B3 = -3.0493E-05
	824	B4 = -1.0924E-07
	825
	826	! temperature dependance: jpl 2006
	827
	828	do ilev = 1,nlayer
	829	temp = min(max(tlay(ilev),200.),400.)
	830	chi = 1./(1. + exp(-1265./temp))
	831	do iw = 1, nw-1
	832	if ((wc(iw) >= 260.) .and. (wc(iw) < 350.)) then
	833	lambda = wc(iw)
	834	sumA = ((((((A7lambda + A6)lambda + A5)*lambda +
	835	$ A4)lambda +A3)lambda + A2)*lambda +
	836	$ A1)*lambda + A0
	837	sumB = (((B4lambda + B3)lambda + B2)*lambda +
	838	$ B1)*lambda + B0
	839	xs = (chisumA + (1. - chi)sumB)*1.e-21
	840	sj(ilev,iw,j_h2o2) = xs
	841	else
	842	sj(ilev,iw,j_h2o2) = xsh2o2(iw)
	843	end if
	844
	845	! optical depth
	846
	847	dt(ilev,iw) = colinc(ilev)rm(ilev)sj(ilev,iw,j_h2o2)
	848	end do
	849	end do
	850
	851	end subroutine seth2o2
	852
	853	!==============================================================================
[2795]	854
	855	subroutine setno2(nd, nlayer, nw, wc, tlay, xsno2, xsno2_220,
	856	$ xsno2_294, yldno2_248, yldno2_298, j_no2,
	857	$ colinc, rm, dt, sj)
	858
[2780]	859	!-----------------------------------------------------------------------------*
	860	!= PURPOSE: =*
	861	!= Set up the no2 temperature-dependent cross-sections and optical depth =*
	862	!-----------------------------------------------------------------------------*
[2795]	863
	864	implicit none
	865
[2780]	866	! input:
	867
[2795]	868	integer :: nd ! number of photolysis rates
	869	integer :: nlayer ! number of vertical layers
	870	integer :: nw ! number of wavelength grid points
	871	integer :: j_no2 ! photolysis index
	872	real, dimension(nw) :: wc ! central wavelength for each interval
	873	real, dimension(nw) :: xsno2, xsno2_220, xsno2_294 ! no2 absorption cross-section at 220-294 k (cm2)
	874	real, dimension(nw) :: yldno2_248, yldno2_298 ! no2 quantum yield at 248-298 k
	875	real, dimension(nlayer) :: tlay ! temperature (k)
	876	real, dimension(nlayer) :: rm ! no2 mixing ratio
	877	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
	878
[2780]	879	! output:
	880
[2795]	881	real, dimension(nlayer,nw) :: dt ! optical depth
	882	real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2)
[2780]	883
	884	! local:
	885
[2795]	886	integer :: ilev, iw
	887	real :: temp, qy
[2780]	888
	889	! temperature dependance: jpl 2006
	890
[2795]	891	do ilev = 1,nlayer
	892	temp = max(220.,min(tlay(ilev),294.))
	893	do iw = 1, nw - 1
	894	if (wc(iw) < 238.) then
	895	sj(ilev,iw,j_no2) = xsno2(iw)
	896	else
	897	sj(ilev,iw,j_no2) = xsno2_220(iw)
	898	$ + (xsno2_294(iw) - xsno2_220(iw))
	899	$ /(294. - 220.)*(temp - 220.)
	900	end if
[2780]	901
	902	! optical depth
	903
[2795]	904	dt(ilev,iw) = colinc(ilev)rm(ilev)sj(ilev,iw,j_no2)
	905	end do
	906	end do
[2780]	907
	908	! quantum yield: jpl 2006
	909
[2795]	910	do ilev = 1,nlayer
	911	temp = max(248.,min(tlay(ilev),298.))
	912	do iw = 1, nw - 1
	913	qy = yldno2_248(iw) + (yldno2_298(iw) - yldno2_248(iw))
	914	$ /(298. - 248.)*(temp - 248.)
	915	sj(ilev,iw,j_no2) = sj(ilev,iw,j_no2)*qy
	916	end do
	917	end do
[2780]	918
[2795]	919	end subroutine setno2
	920
[2780]	921	!==============================================================================
	922
[2974]	923	subroutine setaer(nlev,zalt,tau,nw,dtaer,omaer,gaer)
[2780]	924
	925	!-----------------------------------------------------------------------------*
	926	!= PURPOSE: =*
	927	!= Set aerosol properties for each specified altitude layer. Properties =*
	928	!= may be wavelength dependent. =*
	929	!-----------------------------------------------------------------------------*
	930
	931	implicit none
	932
	933	! input
	934
[2974]	935	integer :: nlev ! number of vertical layers
	936	integer :: nw ! number of wavelength grid points
	937	real, dimension(nlev) :: zalt ! altitude (km)
	938	real :: tau ! integrated aerosol optical depth at the surface
[2780]	939
	940	! output
	941
[2974]	942	real, dimension(nlev,nw) :: dtaer ! aerosol optical depth
	943	real, dimension(nlev,nw) :: omaer ! aerosol single scattering albedo
	944	real, dimension(nlev,nw) :: gaer ! aerosol asymmetry parameter
[2780]	945
	946	! local
	947
[2974]	948	integer :: ilev, iw
	949	real, dimension(nlev) :: aer ! dust extinction
[2780]	950	real :: omega, g, scaleh, gamma
	951	real :: dz, tautot, q0
	952
	953	omega = 0.622 ! single scattering albedo : wolff et al.(2010) at 258 nm
	954	g = 0.88 ! asymmetry factor : mateshvili et al. (2007) at 210 nm
	955	scaleh = 10. ! scale height (km)
	956	gamma = 0.03 ! conrath parameter
	957
[2795]	958	dtaer(:,:) = 0.
	959	omaer(:,:) = 0.
	960	gaer(:,:) = 0.
	961
	962	omaer(:,:) = omega
	963	gaer(:,:) =g
[2780]	964
[2795]	965	! optical depth profile:
	966
	967	tautot = 0.
[2974]	968	do ilev = 1, nlev-1
	969	dz = zalt(ilev+1) - zalt(ilev)
	970	tautot = tautot + exp(gamma(1. - exp(zalt(ilev)/scaleh)))dz
[2795]	971	end do
	972
	973	q0 = tau/tautot
[2974]	974	do ilev = 1, nlev-1
	975	dz = zalt(ilev+1) - zalt(ilev)
	976	dtaer(ilev,:) = q0exp(gamma(1. - exp(zalt(ilev)/scaleh)))*dz
	977	omaer(ilev,:) = omega
	978	gaer(ilev,:) = g
[2795]	979	end do
	980
[2780]	981	end subroutine setaer
	982
	983	!==============================================================================
	984	SUBROUTINE setcld(nz,z,nw,wl,dtcld,omcld,gcld)
	985
	986	-----------------------------------------------------------------------------
	987	= PURPOSE: =
	988	= Set cloud properties for each specified altitude layer. Properties =
	989	= may be wavelength dependent. =
	990	-----------------------------------------------------------------------------
	991	= PARAMETERS: =
	992	= NZ - INTEGER, number of specified altitude levels in the working (I)=
	993	= grid =
	994	= Z - REAL, specified altitude working grid (km) (I)=
	995	= NW - INTEGER, number of specified intervals + 1 in working (I)=
	996	= wavelength grid =
	997	= WL - REAL, vector of lower limits of wavelength intervals in (I)=
	998	= working wavelength grid =
	999	= DTCLD - REAL, optical depth due to absorption by clouds at each (O)=
	1000	= altitude and wavelength =
	1001	= OMCLD - REAL, single scattering albedo due to clouds at each (O)=
	1002	= defined altitude and wavelength =
	1003	= GCLD - REAL, cloud asymmetry factor at each defined altitude and (O)=
	1004	= wavelength =
	1005	-----------------------------------------------------------------------------
	1006	= EDIT HISTORY: =
	1007	= 12/94 Bug fix =
	1008	-----------------------------------------------------------------------------
	1009	= This program is free software; you can redistribute it and/or modify =
	1010	= it under the terms of the GNU General Public License as published by the =
	1011	= Free Software Foundation; either version 2 of the license, or (at your =
	1012	= option) any later version. =
	1013	= The TUV package is distributed in the hope that it will be useful, but =
	1014	= WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTIBI- =
	1015	= LITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public =
	1016	= License for more details. =
	1017	= To obtain a copy of the GNU General Public License, write to: =
	1018	= Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. =
	1019	-----------------------------------------------------------------------------
	1020	= To contact the authors, please mail to: =
	1021	= Sasha Madronich, NCAR/ACD, P.O.Box 3000, Boulder, CO, 80307-3000, USA or =
	1022	= send email to: sasha@ucar.edu =
	1023	-----------------------------------------------------------------------------
	1024	= Copyright (C) 1994,95,96 University Corporation for Atmospheric Research =
	1025	-----------------------------------------------------------------------------
	1026
	1027	IMPLICIT NONE
	1028
	1029	INTEGER kdata
	1030	INTEGER kout
	1031	PARAMETER(kdata=12)
	1032	PARAMETER(kout=53)
	1033
	1034	* input: (grids)
	1035	REAL wl(nw)
	1036	REAL z(nz)
	1037	INTEGER nz
	1038	INTEGER nw
	1039
	1040	* Output:
	1041	REAL dtcld(nz,nw), omcld(nz,nw), gcld(nz,nw)
	1042
	1043	* local:
	1044
	1045	logical clouds
	1046
	1047	* specified data:
	1048	REAL zd(kdata), cd(kdata), omd(kdata), gd(kdata)
	1049	REAL womd(kdata), wgd(kdata)
	1050
	1051	* other:
	1052	REAL cz(nz)
	1053	REAL omz(nz)
	1054	REAL gz(nz)
	1055	INTEGER i, iw, n
	1056
	1057	*_______________________________________________________________________
	1058	c
	1059	c initialize
	1060	c
	1061	do iw = 1, nw-1
	1062	do i = 1, nz-1
	1063	dtcld(i,iw) = 0.
	1064	omcld(i,iw) = 0.
	1065	gcld(i,iw) = 0.
[2795]	1066	end do
[2780]	1067	end do
[2795]	1068	c
[2780]	1069	c if you do not want any clouds, set clouds = .false.
[2795]	1070	c
[2780]	1071	clouds = .true.
[2795]	1072	c
[2780]	1073	if (.not. clouds) then
	1074	write(kout,*) 'no clouds'
	1075	return
	1076	end if
[2795]	1077	c
[2780]	1078	* cloud properties are set for each layer (not each level)
	1079
	1080	* Set as many clouds as want here:
	1081	* First choose a cloud grid, zd(n), in km above sea level
	1082	* Can allow altitude variation of omega, g:
	1083
	1084	n = 12
	1085
	1086	zd(1) = 0.
	1087	cd(1) = 0.
	1088	zd(2) = 30.
	1089	cd(2) = 0.25
	1090	zd(3) = 48.
	1091	cd(3) = 5.84
	1092	zd(4) = 50.
	1093	cd(4) = 5.48
	1094	zd(5) = 54.
	1095	cd(5) = 3.79
	1096	zd(6) = 57.
	1097	cd(6) = 2.1
	1098	zd(7) = 60.
	1099	cd(7) = 3.44
	1100	zd(8) = 62.
	1101	cd(8) = 5.0
	1102	zd(9) = 65.
	1103	cd(9) = 3.48
	1104	zd(10) = 70.
	1105	cd(10) = 0.8
	1106	zd(11) = 78.
	1107	cd(11) = 0.2
	1108	zd(12) = 80.
	1109	cd(12) = 0.
	1110
	1111	do i = 1,n
	1112	omd(i) = 0.999999 ! zhang et al., icarus, 2011
	1113	gd(i) = 0.74 ! zhang et al., icarus, 2011
	1114	end do
	1115
	1116	******************
	1117
	1118	* compute integrals and averages over grid layers:
	1119	* for g and omega, use averages weigthed by optical depth
	1120
[2795]	1121	! DO 11, i = 1, n !*** CHANGED!!See header!!***
[2780]	1122	DO 11, i = 1, n-1
	1123	womd(i) = omd(i) * cd(i)
	1124	wgd(i) = gd(i) * cd(i)
	1125	11 CONTINUE
	1126	CALL inter3(nz,z,cz, n, zd,cd,0)
	1127	CALL inter3(nz,z,omz, n, zd,womd,0)
	1128	CALL inter3(nz,z,gz , n, zd,wgd,0)
	1129
[2795]	1130
[2780]	1131	! Imprimer Cz et imprimer cd pour comparer
	1132
[2795]	1133
[2780]	1134	DO 15, i = 1, nz-1
	1135	IF (cz(i) .GT. 0.) THEN
	1136	omz(i) = omz(i)/cz(i)
	1137	gz(i) = gz(i) /cz(i)
	1138	ELSE
	1139	omz(i) = 1.
	1140	gz(i) = 0.
	1141	ENDIF
	1142	15 CONTINUE
	1143
	1144	* assign at all wavelengths
	1145	* (can move wavelength loop outside if want to vary with wavelength)
	1146
	1147	DO 17, iw = 1, nw-1
	1148	DO 16, i = 1, nz-1
	1149	dtcld(i,iw) = cz(i)
	1150	omcld(i,iw) = omz(i)
	1151	gcld (i,iw) = gz(i)
	1152	16 CONTINUE
	1153	17 CONTINUE
	1154
[2795]	1155	*_______________________________________________________________________
	1156
[2780]	1157	RETURN
[2795]	1158	END
[2780]	1159
	1160	!==============================================================================
	1161
	1162	subroutine sphers(nlev, z, zen, dsdh, nid)
	1163
	1164	!-----------------------------------------------------------------------------*
	1165	!= PURPOSE: =*
	1166	!= Calculate slant path over vertical depth ds/dh in spherical geometry. =*
	1167	!= Calculation is based on: A.Dahlback, and K.Stamnes, A new spheric model =*
	1168	!= for computing the radiation field available for photolysis and heating =*
	1169	!= at twilight, Planet.Space Sci., v39, n5, pp. 671-683, 1991 (Appendix B) =*
	1170	!-----------------------------------------------------------------------------*
	1171	!= PARAMETERS: =*
	1172	!= NZ - INTEGER, number of specified altitude levels in the working (I)=*
	1173	!= grid =*
	1174	!= Z - REAL, specified altitude working grid (km) (I)=*
	1175	!= ZEN - REAL, solar zenith angle (degrees) (I)=*
	1176	!= DSDH - REAL, slant path of direct beam through each layer crossed (O)=*
	1177	!= when travelling from the top of the atmosphere to layer i; =*
	1178	!= DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 =*
	1179	!= NID - INTEGER, number of layers crossed by the direct beam when (O)=*
	1180	!= travelling from the top of the atmosphere to layer i; =*
	1181	!= NID(i), i = 0..NZ-1 =*
	1182	!-----------------------------------------------------------------------------*
	1183
	1184	implicit none
	1185
	1186	! input
	1187
	1188	integer, intent(in) :: nlev
	1189	real, dimension(nlev), intent(in) :: z
	1190	real, intent(in) :: zen
	1191
	1192	! output
	1193
	1194	INTEGER nid(0:nlev)
	1195	REAL dsdh(0:nlev,nlev)
	1196
	1197	! more program constants
	1198
	1199	REAL re, ze(nlev)
	1200	REAL dr
	1201	real radius
	1202	parameter (radius = 6052.)
	1203
	1204	! local
	1205
	1206	real :: pi, zenrad, rpsinz, rj, rjp1, dsj, dhj, ga, gb, sm
	1207	integer :: i, j, k, id, nlay
	1208
	1209	REAL zd(0:nlev-1)
	1210
	1211	!-----------------------------------------------------------------------------
	1212
	1213	pi = acos(-1.0)
	1214	dr = pi/180.
	1215	zenrad = zen*dr
	1216
	1217	! number of layers:
	1218
	1219	nlay = nlev - 1
	1220
	1221	! include the elevation above sea level to the radius of Venus:
	1222
	1223	re = radius + z(1)
	1224
	1225	! correspondingly z changed to the elevation above Venus surface:
	1226
	1227	DO k = 1, nlev
	1228	ze(k) = z(k) - z(1)
	1229	END DO
	1230
	1231	! inverse coordinate of z
	1232
	1233	zd(0) = ze(nlev)
	1234	DO k = 1, nlay
	1235	zd(k) = ze(nlev - k)
	1236	END DO
	1237
	1238	! initialise dsdh(i,j), nid(i)
	1239
	1240	nid(:) = 0.
	1241	dsdh(:,:) = 0.
	1242
	1243	! calculate ds/dh of every layer
	1244
	1245	do i = 0,nlay
	1246	rpsinz = (re + zd(i))*sin(zenrad)
	1247
	1248	IF ( (zen .GT. 90.0) .AND. (rpsinz .LT. re) ) THEN
	1249	nid(i) = -1
	1250	ELSE
	1251
	1252	! Find index of layer in which the screening height lies
	1253
	1254	id = i
	1255	if (zen > 90.) then
	1256	do j = 1,nlay
	1257	IF( (rpsinz .LT. ( zd(j-1) + re ) ) .AND.
	1258	$ (rpsinz .GE. ( zd(j) + re )) ) id = j
	1259	end do
	1260	end if
	1261
	1262	do j = 1,id
	1263	sm = 1.0
	1264	IF (j .EQ. id .AND. id .EQ. i .AND. zen .GT. 90.0)
	1265	$ sm = -1.0
	1266
	1267	rj = re + zd(j-1)
	1268	rjp1 = re + zd(j)
	1269
	1270	dhj = zd(j-1) - zd(j)
	1271
	1272	ga = rjrj - rpsinzrpsinz
	1273	gb = rjp1rjp1 - rpsinzrpsinz
	1274
	1275	ga = max(ga, 0.)
	1276	gb = max(gb, 0.)
	1277
	1278	IF (id.GT.i .AND. j.EQ.id) THEN
	1279	dsj = sqrt(ga)
	1280	ELSE
	1281	dsj = sqrt(ga) - sm*sqrt(gb)
	1282	END IF
	1283	dsdh(i,j) = dsj/dhj
	1284	end do
	1285	nid(i) = id
	1286	end if
	1287	end do ! i = 0,nlay
	1288
	1289	end subroutine sphers
	1290
	1291	!==============================================================================
	1292
	1293	SUBROUTINE rtlink(nlev, nw, iw, ag, zen, dsdh, nid, dtrl,
	1294	$ dagas, dtcld, omcld, gcld, dtaer, omaer, gaer,
	1295	$ edir, edn, eup, fdir, fdn, fup)
	1296
	1297	implicit none
	1298
	1299	! input
	1300
	1301	integer, intent(in) :: nlev, nw, iw ! number of wavelength grid points
	1302	REAL ag
	1303	REAL zen
	1304	REAL dsdh(0:nlev,nlev)
	1305	INTEGER nid(0:nlev)
	1306
	1307	REAL dtrl(nlev,nw)
	1308	REAL dagas(nlev,nw)
	1309	REAL dtcld(nlev,nw), omcld(nlev,nw), gcld(nlev,nw)
	1310	REAL dtaer(nlev,nw), omaer(nlev,nw), gaer(nlev,nw)
	1311
	1312	! output
	1313
	1314	REAL edir(nlev), edn(nlev), eup(nlev)
	1315	REAL fdir(nlev), fdn(nlev), fup(nlev)
	1316
	1317	! local:
	1318
	1319	REAL dt(nlev), om(nlev), g(nlev)
	1320	REAL dtabs,dtsct,dscld,dsaer,dacld,daaer
	1321	INTEGER i, ii
	1322	real, parameter :: largest = 1.e+36
	1323
	1324	! specific two ps2str
	1325
	1326	REAL ediri(nlev), edni(nlev), eupi(nlev)
	1327	REAL fdiri(nlev), fdni(nlev), fupi(nlev)
	1328
	1329	logical, save :: delta = .true.
	1330
	1331	!$OMP THREADPRIVATE(delta)
	1332
	1333	!_______________________________________________________________________
	1334
	1335	! initialize:
	1336
	1337	do i = 1, nlev
	1338	fdir(i) = 0.
	1339	fup(i) = 0.
	1340	fdn(i) = 0.
	1341	edir(i) = 0.
	1342	eup(i) = 0.
	1343	edn(i) = 0.
	1344	end do
	1345
	1346	do i = 1, nlev - 1
	1347	dscld = dtcld(i,iw)*omcld(i,iw)
	1348	dacld = dtcld(i,iw)*(1.-omcld(i,iw))
	1349
	1350	dsaer = dtaer(i,iw)*omaer(i,iw)
	1351	daaer = dtaer(i,iw)*(1.-omaer(i,iw))
	1352
	1353	dtsct = dtrl(i,iw) + dscld + dsaer
	1354	dtabs = dagas(i,iw) + dacld + daaer
	1355
	1356	dtabs = amax1(dtabs,1./largest)
	1357	dtsct = amax1(dtsct,1./largest)
	1358
	1359	! invert z-coordinate:
	1360
	1361	ii = nlev - i
	1362	dt(ii) = dtsct + dtabs
	1363	om(ii) = dtsct/(dtsct + dtabs)
	1364	IF(dtsct .EQ. 1./largest) om(ii) = 1./largest
	1365	g(ii) = (gcld(i,iw)*dscld +
	1366	$ gaer(i,iw)*dsaer)/dtsct
	1367	end do
	1368
	1369	! call rt routine:
	1370
	1371	call ps2str(nlev, zen, ag, dt, om, g,
	1372	$ dsdh, nid, delta,
	1373	$ fdiri, fupi, fdni, ediri, eupi, edni)
	1374
	1375	! output (invert z-coordinate)
	1376
	1377	do i = 1, nlev
	1378	ii = nlev - i + 1
	1379	fdir(i) = fdiri(ii)
	1380	fup(i) = fupi(ii)
	1381	fdn(i) = fdni(ii)
	1382	edir(i) = ediri(ii)
	1383	eup(i) = eupi(ii)
	1384	edn(i) = edni(ii)
	1385	end do
	1386
	1387	end subroutine rtlink
	1388
	1389	=============================================================================
	1390
	1391	subroutine ps2str(nlev,zen,rsfc,tauu,omu,gu,
	1392	$ dsdh, nid, delta,
	1393	$ fdr, fup, fdn, edr, eup, edn)
	1394
	1395	!-----------------------------------------------------------------------------*
	1396	!= PURPOSE: =*
	1397	!= Solve two-stream equations for multiple layers. The subroutine is based =*
	1398	!= on equations from: Toon et al., J.Geophys.Res., v94 (D13), Nov 20, 1989.=*
	1399	!= It contains 9 two-stream methods to choose from. A pseudo-spherical =*
	1400	!= correction has also been added. =*
	1401	!-----------------------------------------------------------------------------*
	1402	!= PARAMETERS: =*
	1403	!= NLEVEL - INTEGER, number of specified altitude levels in the working (I)=*
	1404	!= grid =*
	1405	!= ZEN - REAL, solar zenith angle (degrees) (I)=*
	1406	!= RSFC - REAL, surface albedo at current wavelength (I)=*
	1407	!= TAUU - REAL, unscaled optical depth of each layer (I)=*
	1408	!= OMU - REAL, unscaled single scattering albedo of each layer (I)=*
	1409	!= GU - REAL, unscaled asymmetry parameter of each layer (I)=*
	1410	!= DSDH - REAL, slant path of direct beam through each layer crossed (I)=*
	1411	!= when travelling from the top of the atmosphere to layer i; =*
	1412	!= DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 =*
	1413	!= NID - INTEGER, number of layers crossed by the direct beam when (I)=*
	1414	!= travelling from the top of the atmosphere to layer i; =*
	1415	!= NID(i), i = 0..NZ-1 =*
	1416	!= DELTA - LOGICAL, switch to use delta-scaling (I)=*
	1417	!= .TRUE. -> apply delta-scaling =*
	1418	!= .FALSE.-> do not apply delta-scaling =*
	1419	!= FDR - REAL, contribution of the direct component to the total (O)=*
	1420	!= actinic flux at each altitude level =*
	1421	!= FUP - REAL, contribution of the diffuse upwelling component to (O)=*
	1422	!= the total actinic flux at each altitude level =*
	1423	!= FDN - REAL, contribution of the diffuse downwelling component to (O)=*
	1424	!= the total actinic flux at each altitude level =*
	1425	!= EDR - REAL, contribution of the direct component to the total (O)=*
	1426	!= spectral irradiance at each altitude level =*
	1427	!= EUP - REAL, contribution of the diffuse upwelling component to (O)=*
	1428	!= the total spectral irradiance at each altitude level =*
	1429	!= EDN - REAL, contribution of the diffuse downwelling component to (O)=*
	1430	= the total spectral irradiance at each altitude level =
	1431	!-----------------------------------------------------------------------------*
	1432
	1433	implicit none
	1434
	1435	! input:
	1436
	1437	INTEGER nlev
	1438	REAL zen, rsfc
	1439	REAL tauu(nlev), omu(nlev), gu(nlev)
	1440	REAL dsdh(0:nlev,nlev)
	1441	INTEGER nid(0:nlev)
	1442	LOGICAL delta
	1443
	1444	! output:
	1445
	1446	REAL fup(nlev),fdn(nlev),fdr(nlev)
	1447	REAL eup(nlev),edn(nlev),edr(nlev)
	1448
	1449	! local:
	1450
	1451	REAL tausla(0:nlev), tauc(0:nlev)
	1452	REAL mu2(0:nlev), mu, sum
	1453
	1454	! internal coefficients and matrix
	1455
	1456	REAL lam(nlev),taun(nlev),bgam(nlev)
	1457	REAL e1(nlev),e2(nlev),e3(nlev),e4(nlev)
	1458	REAL cup(nlev),cdn(nlev),cuptn(nlev),cdntn(nlev)
	1459	REAL mu1(nlev)
	1460	INTEGER row
	1461	REAL a(2nlev),b(2nlev),d(2nlev),e(2nlev),y(2*nlev)
	1462
	1463	! other:
	1464
	1465	REAL pifs, fdn0
	1466	REAL gi(nlev), omi(nlev), tempg
	1467	REAL f, g, om
	1468	REAL gam1, gam2, gam3, gam4
	1469	real, parameter :: largest = 1.e+36
	1470	real, parameter :: precis = 1.e-7
	1471
	1472	! For calculations of Associated Legendre Polynomials for GAMA1,2,3,4
	1473	! in delta-function, modified quadrature, hemispheric constant,
	1474	! Hybrid modified Eddington-delta function metods, p633,Table1.
	1475	! W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
	1476	! W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440,
	1477	! uncomment the following two lines and the appropriate statements further
	1478	! down.
	1479	! REAL YLM0, YLM2, YLM4, YLM6, YLM8, YLM10, YLM12, YLMS, BETA0,
	1480	! > BETA1, BETAn, amu1, subd
	1481
	1482	REAL expon, expon0, expon1, divisr, temp, up, dn
	1483	REAL ssfc
	1484	INTEGER nlayer, mrows, lev
	1485
	1486	INTEGER i, j
	1487
	1488	! Some additional program constants:
	1489
	1490	real pi, dr
	1491	REAL eps
	1492	PARAMETER (eps = 1.E-3)
	1493	!_______________________________________________________________________
	1494
	1495	! MU = cosine of solar zenith angle
	1496	! RSFC = surface albedo
	1497	! TAUU = unscaled optical depth of each layer
	1498	! OMU = unscaled single scattering albedo
	1499	! GU = unscaled asymmetry factor
	1500	! KLEV = max dimension of number of layers in atmosphere
	1501	! NLAYER = number of layers in the atmosphere
	1502	! NLEVEL = nlayer + 1 = number of levels
	1503
	1504	! initial conditions: pi*solar flux = 1; diffuse incidence = 0
	1505
	1506	pifs = 1.
	1507	fdn0 = 0.
	1508
	1509	nlayer = nlev - 1
	1510
	1511	pi = acos(-1.)
	1512	dr = pi/180.
	1513	mu = COS(zen*dr)
	1514
	1515	!************* compute coefficients for each layer:
	1516	! GAM1 - GAM4 = 2-stream coefficients, different for different approximations
	1517	! EXPON0 = calculation of e when TAU is zero
	1518	! EXPON1 = calculation of e when TAU is TAUN
	1519	! CUP and CDN = calculation when TAU is zero
	1520	! CUPTN and CDNTN = calc. when TAU is TAUN
	1521	! DIVISR = prevents division by zero
	1522	do j = 0, nlev
	1523	tauc(j) = 0.
	1524	tausla(j) = 0.
	1525	mu2(j) = 1./SQRT(largest)
	1526	end do
	1527
	1528	IF (.NOT. delta) THEN
	1529	DO i = 1, nlayer
	1530	gi(i) = gu(i)
	1531	omi(i) = omu(i)
	1532	taun(i) = tauu(i)
	1533	END DO
	1534	ELSE
	1535
	1536	! delta-scaling. Have to be done for delta-Eddington approximation,
	1537	! delta discrete ordinate, Practical Improved Flux Method, delta function,
	1538	! and Hybrid modified Eddington-delta function methods approximations
	1539
	1540	DO i = 1, nlayer
	1541	f = gu(i)*gu(i)
	1542	gi(i) = (gu(i) - f)/(1 - f)
	1543	omi(i) = (1 - f)omu(i)/(1 - omu(i)f)
	1544	taun(i) = (1 - omu(i)f)tauu(i)
	1545	END DO
	1546	END IF
	1547
	1548	! calculate slant optical depth at the top of the atmosphere when zen>90.
	1549	! in this case, higher altitude of the top layer is recommended.
	1550
	1551	IF (zen .GT. 90.0) THEN
	1552	IF (nid(0) .LT. 0) THEN
	1553	tausla(0) = largest
	1554	ELSE
	1555	sum = 0.0
	1556	DO j = 1, nid(0)
	1557	sum = sum + 2.taun(j)dsdh(0,j)
	1558	END DO
	1559	tausla(0) = sum
	1560	END IF
	1561	END IF
	1562
	1563	DO 11, i = 1, nlayer
	1564	g = gi(i)
	1565	om = omi(i)
	1566	tauc(i) = tauc(i-1) + taun(i)
	1567
	1568	! stay away from 1 by precision. For g, also stay away from -1
	1569
	1570	tempg = AMIN1(abs(g),1. - precis)
	1571	g = SIGN(tempg,g)
	1572	om = AMIN1(om,1.-precis)
	1573
	1574	! calculate slant optical depth
	1575
	1576	IF (nid(i) .LT. 0) THEN
	1577	tausla(i) = largest
	1578	ELSE
	1579	sum = 0.0
	1580	DO j = 1, MIN(nid(i),i)
	1581	sum = sum + taun(j)*dsdh(i,j)
	1582	END DO
	1583	DO j = MIN(nid(i),i)+1,nid(i)
	1584	sum = sum + 2.taun(j)dsdh(i,j)
	1585	END DO
	1586	tausla(i) = sum
	1587	IF (tausla(i) .EQ. tausla(i-1)) THEN
	1588	mu2(i) = SQRT(largest)
	1589	ELSE
	1590	mu2(i) = (tauc(i)-tauc(i-1))/(tausla(i)-tausla(i-1))
	1591	mu2(i) = SIGN( AMAX1(ABS(mu2(i)),1./SQRT(largest)),
	1592	$ mu2(i) )
	1593	END IF
	1594	END IF
	1595
	1596	!** the following gamma equations are from pg 16,289, Table 1
	1597	!** save mu1 for each approx. for use in converting irradiance to actinic flux
	1598
	1599	! Eddington approximation(Joseph et al., 1976, JAS, 33, 2452):
	1600
	1601	c gam1 = (7. - om(4. + 3.g))/4.
	1602	c gam2 = -(1. - om(4. - 3.g))/4.
	1603	c gam3 = (2. - 3.gmu)/4.
	1604	c gam4 = 1. - gam3
	1605	c mu1(i) = 0.5
	1606
	1607	* quadrature (Liou, 1973, JAS, 30, 1303-1326; 1974, JAS, 31, 1473-1475):
	1608
	1609	c gam1 = 1.7320508(2. - om(1. + g))/2.
	1610	c gam2 = 1.7320508om(1. - g)/2.
	1611	c gam3 = (1. - 1.7320508gmu)/2.
	1612	c gam4 = 1. - gam3
	1613	c mu1(i) = 1./sqrt(3.)
	1614
	1615	* hemispheric mean (Toon et al., 1089, JGR, 94, 16287):
	1616
	1617	gam1 = 2. - om*(1. + g)
	1618	gam2 = om*(1. - g)
	1619	gam3 = (2. - g*mu)/4.
	1620	gam4 = 1. - gam3
	1621	mu1(i) = 0.5
	1622
	1623	* PIFM (Zdunkovski et al.,1980, Conrib.Atmos.Phys., 53, 147-166):
	1624	c GAM1 = 0.25(8. - OM(5. + 3.*G))
	1625	c GAM2 = 0.75OM(1.-G)
	1626	c GAM3 = 0.25(2.-3.G*MU)
	1627	c GAM4 = 1. - GAM3
	1628	c mu1(i) = 0.5
	1629
	1630	* delta discrete ordinates (Schaller, 1979, Contrib.Atmos.Phys, 52, 17-26):
	1631	c GAM1 = 0.51.7320508(2. - OM*(1. + G))
	1632	c GAM2 = 0.51.7320508OM*(1.-G)
	1633	c GAM3 = 0.5(1.-1.7320508G*MU)
	1634	c GAM4 = 1. - GAM3
	1635	c mu1(i) = 1./sqrt(3.)
	1636
	1637	* Calculations of Associated Legendre Polynomials for GAMA1,2,3,4
	1638	* in delta-function, modified quadrature, hemispheric constant,
	1639	* Hybrid modified Eddington-delta function metods, p633,Table1.
	1640	* W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
	1641	* W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440
	1642	c YLM0 = 2.
	1643	c YLM2 = -3.GMU
	1644	c YLM4 = 0.875G3MU(5.MU**2-3.)
	1645	c YLM6=-0.171875G5MU(15.-70.MU*2+63.MU**4)
	1646	c YLM8=+0.073242G7MU(-35.+315.MU*2-693.MU**4
	1647	c +429.MU**6)
	1648	c YLM10=-0.008118G9MU(315.-4620.MU*2+18018.MU**4
	1649	c -25740.MU*6+12155.MU**8)
	1650	c YLM12=0.003685G11MU(-693.+15015.MU*2-90090.MU**4
	1651	c +218790.MU*6-230945.MU*8+88179.MU**10)
	1652	c YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
	1653	c YLMS=0.25*YLMS
	1654	c BETA0 = YLMS
	1655	c
	1656	c amu1=1./1.7320508
	1657	c YLM0 = 2.
	1658	c YLM2 = -3.Gamu1
	1659	c YLM4 = 0.875G3amu1(5.amu1**2-3.)
	1660	c YLM6=-0.171875G5amu1(15.-70.amu1*2+63.amu1**4)
	1661	c YLM8=+0.073242G7amu1(-35.+315.amu1*2-693.amu1**4
	1662	c +429.amu1**6)
	1663	c YLM10=-0.008118G9amu1(315.-4620.amu1*2+18018.amu1**4
	1664	c -25740.amu1*6+12155.amu1**8)
	1665	c YLM12=0.003685G11amu1(-693.+15015.amu1*2-90090.amu1**4
	1666	c +218790.amu1*6-230945.amu1*8+88179.amu1**10)
	1667	c YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
	1668	c YLMS=0.25*YLMS
	1669	c BETA1 = YLMS
	1670	c
	1671	c BETAn = 0.25(2. - 1.5G-0.21875G3-0.085938G**5
	1672	c -0.045776G**7)
	1673
	1674
	1675	* Hybrid modified Eddington-delta function(Meador and Weaver,1980,JAS,37,630):
	1676	c subd=4.(1.-GG*(1.-MU))
	1677	c GAM1 = (7.-3.GG-OM(4.+3.G)+OMGG(4.BETA0+3.*G))/subd
	1678	c GAM2 =-(1.-GG-OM(4.-3.G)-OMGG(4.BETA0+3.G-4.))/subd
	1679	c GAM3 = BETA0
	1680	c GAM4 = 1. - GAM3
	1681	c mu1(i) = (1. - gg(1.- mu) )/(2. - g*g)
	1682
	1683	*****
	1684	* delta function (Meador, and Weaver, 1980, JAS, 37, 630):
	1685	c GAM1 = (1. - OM*(1. - beta0))/MU
	1686	c GAM2 = OM*BETA0/MU
	1687	c GAM3 = BETA0
	1688	c GAM4 = 1. - GAM3
	1689	c mu1(i) = mu
	1690	*****
	1691	* modified quadrature (Meador, and Weaver, 1980, JAS, 37, 630):
	1692	c GAM1 = 1.7320508(1. - OM(1. - beta1))
	1693	c GAM2 = 1.7320508OMbeta1
	1694	c GAM3 = BETA0
	1695	c GAM4 = 1. - GAM3
	1696	c mu1(i) = 1./sqrt(3.)
	1697
	1698	* hemispheric constant (Toon et al., 1989, JGR, 94, 16287):
	1699	c GAM1 = 2.(1. - OM(1. - betan))
	1700	c GAM2 = 2.OMBETAn
	1701	c GAM3 = BETA0
	1702	c GAM4 = 1. - GAM3
	1703	c mu1(i) = 0.5
	1704
	1705	*****
	1706
	1707	* lambda = pg 16,290 equation 21
	1708	* big gamma = pg 16,290 equation 22
	1709	* if gam2 = 0., then bgam = 0.
	1710
	1711	lam(i) = sqrt(gam1gam1 - gam2gam2)
	1712
	1713	IF (gam2 .NE. 0.) THEN
	1714	bgam(i) = (gam1 - lam(i))/gam2
	1715	ELSE
	1716	bgam(i) = 0.
	1717	END IF
	1718
	1719	expon = EXP(-lam(i)*taun(i))
	1720
	1721	* e1 - e4 = pg 16,292 equation 44
	1722
	1723	e1(i) = 1. + bgam(i)*expon
	1724	e2(i) = 1. - bgam(i)*expon
	1725	e3(i) = bgam(i) + expon
	1726	e4(i) = bgam(i) - expon
	1727
	1728	* the following sets up for the C equations 23, and 24
	1729	* found on page 16,290
	1730	* prevent division by zero (if LAMBDA=1/MU, shift 1/MU^2 by EPS = 1.E-3
	1731	* which is approx equiv to shifting MU by 0.5EPS (MU)**3
	1732
	1733	expon0 = EXP(-tausla(i-1))
	1734	expon1 = EXP(-tausla(i))
	1735
	1736	divisr = lam(i)lam(i) - 1./(mu2(i)mu2(i))
	1737	temp = AMAX1(eps,abs(divisr))
	1738	divisr = SIGN(temp,divisr)
	1739
	1740	up = ompifs((gam1 - 1./mu2(i))gam3 + gam4gam2)/divisr
	1741	dn = ompifs((gam1 + 1./mu2(i))gam4 + gam2gam3)/divisr
	1742
	1743	* cup and cdn are when tau is equal to zero
	1744	* cuptn and cdntn are when tau is equal to taun
	1745
	1746	cup(i) = up*expon0
	1747	cdn(i) = dn*expon0
	1748	cuptn(i) = up*expon1
	1749	cdntn(i) = dn*expon1
	1750
	1751	11 CONTINUE
	1752
	1753	*************** set up matrix ****
	1754	* ssfc = pg 16,292 equation 37 where pi Fs is one (unity).
	1755
	1756	ssfc = rsfcmuEXP(-tausla(nlayer))*pifs
	1757
	1758	* MROWS = the number of rows in the matrix
	1759
	1760	mrows = 2*nlayer
	1761
	1762	* the following are from pg 16,292 equations 39 - 43.
	1763	* set up first row of matrix:
	1764
	1765	i = 1
	1766	a(1) = 0.
	1767	b(1) = e1(i)
	1768	d(1) = -e2(i)
	1769	e(1) = fdn0 - cdn(i)
	1770
	1771	row=1
	1772
	1773	* set up odd rows 3 thru (MROWS - 1):
	1774
	1775	i = 0
	1776	DO 20, row = 3, mrows - 1, 2
	1777	i = i + 1
	1778	a(row) = e2(i)e3(i) - e4(i)e1(i)
	1779	b(row) = e1(i)e1(i + 1) - e3(i)e3(i + 1)
	1780	d(row) = e3(i)e4(i + 1) - e1(i)e2(i + 1)
	1781	e(row) = e3(i)*(cup(i + 1) - cuptn(i)) +
	1782	$ e1(i)*(cdntn(i) - cdn(i + 1))
	1783	20 CONTINUE
	1784
	1785	* set up even rows 2 thru (MROWS - 2):
	1786
	1787	i = 0
	1788	DO 30, row = 2, mrows - 2, 2
	1789	i = i + 1
	1790	a(row) = e2(i + 1)e1(i) - e3(i)e4(i + 1)
	1791	b(row) = e2(i)e2(i + 1) - e4(i)e4(i + 1)
	1792	d(row) = e1(i + 1)e4(i + 1) - e2(i + 1)e3(i + 1)
	1793	e(row) = (cup(i + 1) - cuptn(i))*e2(i + 1) -
	1794	$ (cdn(i + 1) - cdntn(i))*e4(i + 1)
	1795	30 CONTINUE
	1796
	1797	* set up last row of matrix at MROWS:
	1798
	1799	row = mrows
	1800	i = nlayer
	1801
	1802	a(row) = e1(i) - rsfc*e3(i)
	1803	b(row) = e2(i) - rsfc*e4(i)
	1804	d(row) = 0.
	1805	e(row) = ssfc - cuptn(i) + rsfc*cdntn(i)
	1806
	1807	* solve tri-diagonal matrix:
	1808
	1809	CALL tridiag(a, b, d, e, y, mrows)
	1810
	1811	**** unfold solution of matrix, compute output fluxes:
	1812
	1813	row = 1
	1814	lev = 1
	1815	j = 1
	1816
	1817	* the following equations are from pg 16,291 equations 31 & 32
	1818
	1819	fdr(lev) = EXP( -tausla(0) )
	1820	edr(lev) = mu * fdr(lev)
	1821	edn(lev) = fdn0
	1822	eup(lev) = y(row)e3(j) - y(row + 1)e4(j) + cup(j)
	1823	fdn(lev) = edn(lev)/mu1(lev)
	1824	fup(lev) = eup(lev)/mu1(lev)
	1825	DO 60, lev = 2, nlayer + 1
	1826	fdr(lev) = EXP(-tausla(lev-1))
	1827	edr(lev) = mu *fdr(lev)
	1828	edn(lev) = y(row)e3(j) + y(row + 1)e4(j) + cdntn(j)
	1829	eup(lev) = y(row)e1(j) + y(row + 1)e2(j) + cuptn(j)
	1830	fdn(lev) = edn(lev)/mu1(j)
	1831	fup(lev) = eup(lev)/mu1(j)
	1832
	1833	row = row + 2
	1834	j = j + 1
	1835	60 CONTINUE
	1836
	1837	end subroutine ps2str
	1838
	1839	=============================================================================
	1840
	1841	subroutine tridiag(a,b,c,r,u,n)
	1842
	1843	!_______________________________________________________________________
	1844	! solves tridiagonal system. From Numerical Recipies, p. 40
	1845	!_______________________________________________________________________
	1846
	1847	IMPLICIT NONE
	1848
	1849	! input:
	1850
	1851	INTEGER n
	1852	REAL a, b, c, r
	1853	DIMENSION a(n),b(n),c(n),r(n)
	1854
	1855	! output:
	1856
	1857	REAL u
	1858	DIMENSION u(n)
	1859
	1860	! local:
	1861
	1862	INTEGER j
	1863
	1864	REAL bet, gam
	1865	DIMENSION gam(n)
	1866	!_______________________________________________________________________
	1867
	1868	IF (b(1) .EQ. 0.) STOP 1001
	1869	bet = b(1)
	1870	u(1) = r(1)/bet
	1871	DO 11, j = 2, n
	1872	gam(j) = c(j - 1)/bet
	1873	bet = b(j) - a(j)*gam(j)
	1874	IF (bet .EQ. 0.) STOP 2002
	1875	u(j) = (r(j) - a(j)*u(j - 1))/bet
	1876	11 CONTINUE
	1877	DO 12, j = n - 1, 1, -1
	1878	u(j) = u(j) - gam(j + 1)*u(j + 1)
	1879	12 CONTINUE
	1880	!_______________________________________________________________________
	1881
	1882	end subroutine tridiag
	1883
	1884	=============================================================================
	1885
	1886
	1887	SUBROUTINE inter3(ng,xg,yg, n,x,y, FoldIn)
	1888	IMPLICIT NONE
	1889
	1890	* input:
	1891	INTEGER n, ng
	1892	REAL xg(ng)
	1893	REAL x(n), y(n)
	1894
	1895	INTEGER FoldIn
	1896
	1897	* output:
	1898	REAL yg(ng)
	1899
	1900	* local:
	1901	REAL a1, a2, sum
	1902	REAL tail
	1903	INTEGER jstart, i, j, k
	1904	*_______________________________________________________________________
	1905
	1906	* check whether flag given is legal
	1907	IF ((FoldIn .NE. 0) .AND. (FoldIn .NE. 1)) THEN
	1908	WRITE(0,*) '>>> ERROR (inter3) <<< Value for FOLDIN invalid. '
	1909	WRITE(0,*) ' Must be 0 or 1'
	1910	STOP
	1911	ENDIF
	1912
	1913	* do interpolation
	1914
	1915	jstart = 1
	1916
	1917	DO 30, i = 1, ng - 1
	1918
	1919	yg(i) = 0.
	1920	sum = 0.
	1921	j = jstart
	1922
	1923	IF (j .LE. n-1) THEN
	1924
	1925	20 CONTINUE
	1926
	1927	IF (x(j+1) .LT. xg(i)) THEN
	1928	jstart = j
	1929	j = j+1
	1930	IF (j .LE. n-1) GO TO 20
	1931	ENDIF
	1932	25 CONTINUE
	1933
	1934	IF ((x(j) .LE. xg(i+1)) .AND. (j .LE. n-1)) THEN
	1935	a1 = AMAX1(x(j),xg(i))
	1936	a2 = AMIN1(x(j+1),xg(i+1))
	1937	sum = sum + y(j) * (a2-a1)/(x(j+1)-x(j))
	1938	j = j+1
	1939	GO TO 25
	1940
	1941	ENDIF
	1942	yg(i) = sum
	1943
	1944	ENDIF
	1945
	1946	30 CONTINUE
	1947
	1948
	1949	* if wanted, integrate data "overhang" and fold back into last bin
	1950
	1951	IF (FoldIn .EQ. 1) THEN
	1952
	1953	j = j-1
	1954	a1 = xg(ng) ! upper limit of last interpolated bin
	1955	a2 = x(j+1) ! upper limit of last input bin considered
	1956
	1957	* do folding only if grids don't match up and there is more input
	1958	IF ((a2 .GT. a1) .OR. (j+1 .LT. n)) THEN
	1959	tail = y(j) * (a2-a1)/(x(j+1)-x(j))
	1960	DO k = j+1, n-1
	1961	tail = tail + y(k) * (x(k+1)-x(k))
	1962	ENDDO
	1963	yg(ng-1) = yg(ng-1) + tail
	1964	ENDIF
	1965
	1966	ENDIF
	1967	*_______________________________________________________________________
	1968
	1969	RETURN
	1970	end subroutine inter3
	1971
	1972	end subroutine photolysis_online

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

Download in other formats:

Original Format