Context Navigation

photolysis_online.F @ 3556

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

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

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

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

File size: 66.9 KB

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

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

Download in other formats:

Original Format