Context Navigation

photolysis_online.F @ 3461

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

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

Download in other formats:

Original Format