Context Navigation

photolysis_online.F @ 3803

Last change on this file since 3803 was 3755, checked in by flefevre, 7 months ago

Implementation of Cl2SO2 chemistry

Cl2SO2 absorption cross-section file must be downloaded from:

https://owncloud.latmos.ipsl.fr/index.php/s/bkTjmcm6EdoxMnT

and placed in INPUT/cross_sections directory.

File size: 68.7 KB

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

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

Download in other formats:

Original Format