Context Navigation

photolysis_online.F @ 2157

Last change on this file since 2157 was 2041, checked in by flefevre, 7 years ago
Photolyse on-line : Mise a jour de la section efficace Rayleigh (Itiaksov et al., 2008)
File size: 51.0 KB

Line
1	!==============================================================================
2
3	subroutine photolysis_online(nlayer, alt, press, temp,
4	$ zmmean, i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,
5	$ i_h2, i_oh, i_ho2, i_h2o2, i_h2o,
6	$ i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm,
7	$ tau, sza, dist_sol, v_phot)
8
9	use photolysis_mod
10
11	implicit none
12
13	! input
14
15	integer, intent(in) :: nesp ! total number of chemical species
16	integer, intent(in) :: nlayer
17	integer, intent(in) :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,
18	$ i_h2, i_oh, i_ho2, i_h2o2, i_h2o,
19	$ i_n, i_n2d, i_no, i_no2, i_n2
20
21	real, dimension(nlayer), intent(in) :: press, temp, zmmean ! pressure (hpa)/temperature (k)/mean molecular mass (g.mol-1)
22	real, dimension(nlayer), intent(in) :: alt ! altitude (km)
23	real, dimension(nlayer,nesp), intent(in) :: rm ! mixing ratios
24	real, intent(in) :: tau ! integrated aerosol optical depth at the surface
25	real, intent(in) :: sza ! solar zenith angle (degrees)
26	real, intent(in) :: dist_sol ! solar distance (au)
27
28	! output
29
30	real (kind = 8), dimension(nlayer,nb_phot_max) :: v_phot ! photolysis rates (s-1)
31
32	! solar flux at mars
33
34	real, dimension(nw) :: fmars ! solar flux (w.m-2.nm-1)
35	real :: factor
36
37	! cross-sections
38
39	real, dimension(nlayer,nw,nphot) :: sj ! general cross-section array (cm2)
40
41	! atmosphere
42
43	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
44	real, dimension(nlayer,nw) :: dtrl ! rayleigh optical depth
45	real, dimension(nlayer,nw) :: dtaer ! aerosol optical depth
46	real, dimension(nlayer,nw) :: omaer ! aerosol single scattering albedo
47	real, dimension(nlayer,nw) :: gaer ! aerosol asymmetry parameter
48	real, dimension(nlayer,nw) :: dtcld ! cloud optical depth
49	real, dimension(nlayer,nw) :: omcld ! cloud single scattering albedo
50	real, dimension(nlayer,nw) :: gcld ! cloud asymmetry parameter
51	real, dimension(nlayer,nw,nabs) :: dtgas ! optical depth for each gas
52	real, dimension(nlayer,nw) :: dagas ! total gas optical depth
53	real, dimension(nlayer) :: edir, edn, eup ! normalised irradiances
54	real, dimension(nlayer) :: fdir, fdn, fup ! normalised actinic fluxes
55	real, dimension(nlayer) :: saflux ! total actinic flux
56
57	integer, dimension(0:nlayer) :: nid
58	real, dimension(0:nlayer,nlayer) :: dsdh
59
60	integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d, j_o3_o,
61	$ j_h2o, j_h2o2, j_ho2, j_h2, j_no, j_no2, j_n2
62
63	integer :: a_o2, a_co2, a_o3, a_h2o, a_h2o2, a_ho2, a_h2, a_no,
64	$ a_no2, a_n2
65	integer :: i, ilay, iw, ialt
66	real :: deltaj
67
68	! absorbing gas numbering
69
70	a_o2 = 1 ! o2
71	a_co2 = 2 ! co2
72	a_o3 = 3 ! o3
73	a_h2o = 4 ! h2o
74	a_h2o2 = 5 ! h2o2
75	a_ho2 = 6 ! ho2
76	a_h2 = 7 ! h2
77	a_no = 8 ! no
78	a_no2 = 9 ! no2
79	a_n2 = 10 ! no2
80
81	! photodissociation rates numbering.
82	! photodissociations must be ordered the same way in subroutine "indice"
83
84	j_o2_o = 1 ! o2 + hv -> o + o
85	j_o2_o1d = 2 ! o2 + hv -> o + o(1d)
86	j_co2_o = 3 ! co2 + hv -> co + o
87	j_co2_o1d = 4 ! co2 + hv -> co + o(1d)
88	j_o3_o1d = 5 ! o3 + hv -> o2 + o(1d)
89	j_o3_o = 6 ! o3 + hv -> o2 + o
90	j_h2o = 7 ! h2o + hv -> h + oh
91	j_h2o2 = 8 ! h2o2 + hv -> oh + oh
92	j_ho2 = 9 ! ho2 + hv -> oh + o
93	j_h2 = 10 ! h2 + hv -> h + h
94	j_no = 11 ! no + hv -> n + o
95	j_no2 = 12 ! no2 + hv -> no + o
96	j_n2 = 13 ! n2 + hv -> n + n
97
98	!==== air column increments and rayleigh optical depth
99
100	call setair(nlayer, nw, wl, wc, press, temp, zmmean, colinc, dtrl)
101
102	!==== set temperature-dependent cross-sections and optical depths
103
104	! o2:
105
106	call seto2(nphot, nlayer, nw, wc, mopt, temp, xso2_150, xso2_200,
107	$ xso2_250, xso2_300, yieldo2, j_o2_o, j_o2_o1d,
108	$ colinc, rm(:,i_o2), dtgas(:,:,a_o2), sj)
109
110	! co2:
111
112	call setco2(nphot, nlayer, nw, wc, temp, xsco2_195, xsco2_295,
113	$ xsco2_370, yieldco2, j_co2_o, j_co2_o1d,
114	$ colinc, rm(:,i_co2), dtgas(:,:,a_co2), sj)
115
116	! o3:
117
118	call seto3(nphot, nlayer, nw, wc, temp, xso3_218, xso3_298,
119	$ j_o3_o, j_o3_o1d,
120	$ colinc, rm(:,i_o3), dtgas(:,:,a_o3), sj)
121
122	! h2o2:
123
124	call seth2o2(nphot, nlayer, nw, wc, temp, xsh2o2, j_h2o2,
125	$ colinc, rm(:,i_h2o2), dtgas(:,:,a_h2o2), sj)
126
127	! no2:
128
129	call setno2(nphot, nlayer, nw, wc, temp, xsno2, xsno2_220,
130	$ xsno2_294, yldno2_248, yldno2_298, j_no2,
131	$ colinc, rm(:,i_no2), dtgas(:,:,a_no2), sj)
132
133	!==== temperature independent optical depths and cross-sections
134
135	! optical depths
136
137	do ilay = 1,nlayer
138	do iw = 1,nw-1
139	dtgas(ilay,iw,a_h2o) = colinc(ilay)rm(ilay,i_h2o)xsh2o(iw)
140	dtgas(ilay,iw,a_ho2) = colinc(ilay)rm(ilay,i_ho2)xsho2(iw)
141	dtgas(ilay,iw,a_h2) = colinc(ilay)rm(ilay,i_h2)xsh2(iw)
142	dtgas(ilay,iw,a_no) = colinc(ilay)rm(ilay,i_no)xsno(iw)
143	dtgas(ilay,iw,a_n2) = colinc(ilay)rm(ilay,i_n2)xsn2(iw)
144	end do
145	end do
146
147	! total gas optical depth
148
149	dagas(:,:) = 0.
150
151	do ilay = 1,nlayer
152	do iw = 1,nw-1
153	do i = 1,nabs
154	dagas(ilay,iw) = dagas(ilay,iw) + dtgas(ilay,iw,i)
155	end do
156	end do
157	end do
158
159	! cross-sections
160
161	do ilay = 1,nlayer
162	do iw = 1,nw-1
163	sj(ilay,iw,j_h2o) = xsh2o(iw) ! h2o
164	sj(ilay,iw,j_ho2) = xsho2(iw) ! ho2
165	sj(ilay,iw,j_h2) = xsh2(iw)*yieldh2(iw) ! h2
166	sj(ilay,iw,j_no) = xsno(iw)*yieldno(iw) ! no
167	sj(ilay,iw,j_n2) = xsn2(iw)*yieldn2(iw) ! n2
168	end do
169	end do
170
171	!==== set aerosol properties and optical depth
172
173	call setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer)
174
175	!==== set cloud properties and optical depth
176
177	call setcld(nlayer,nw,dtcld,omcld,gcld)
178
179	!==== slant path lengths in spherical geometry
180
181	call sphers(nlayer,alt,sza,dsdh,nid)
182
183	!==== solar flux at mars
184
185	factor = (1./dist_sol)**2.
186	do iw = 1,nw-1
187	fmars(iw) = f(iw)*factor
188	end do
189
190	!==== initialise photolysis rates
191
192	v_phot(:,1:nphot) = 0.
193
194	!==== start of wavelength lopp
195
196	do iw = 1,nw-1
197
198	! monochromatic radiative transfer. outputs are:
199	! normalized irradiances edir(nlayer), edn(nlayer), eup(nlayer)
200	! normalized actinic fluxes fdir(nlayer), fdn(nlayer), fup(nlayer)
201	! where
202	! dir = direct beam, dn = down-welling diffuse, up = up-welling diffuse
203
204	call rtlink(nlayer, nw, iw, albedo(iw), sza, dsdh, nid, dtrl,
205	$ dagas, dtcld, omcld, gcld, dtaer, omaer, gaer,
206	$ edir, edn, eup, fdir, fdn, fup)
207
208
209	! spherical actinic flux
210
211	do ilay = 1,nlayer
212	saflux(ilay) = fmars(iw)*(fdir(ilay) + fdn(ilay) + fup(ilay))
213	end do
214
215	! photolysis rate integration
216
217	do i = 1,nphot
218	do ilay = 1,nlayer
219	deltaj = saflux(ilay)*sj(ilay,iw,i)
220	v_phot(ilay,i) = v_phot(ilay,i) + deltaj*(wu(iw)-wl(iw))
221	end do
222	end do
223
224	! eliminate small values
225
226	where (v_phot(:,1:nphot) < 1.e-30)
227	v_phot(:,1:nphot) = 0.
228	end where
229
230	end do ! iw
231
232	contains
233
234	!==============================================================================
235
236	subroutine setair(nlev, nw, wl, wc, press, temp, zmmean,
237	$ colinc, dtrl)
238
239	-----------------------------------------------------------------------------
240	= PURPOSE: =
241	= computes air column increments and rayleigh optical depth =
242	-----------------------------------------------------------------------------
243
244	implicit none
245
246	! input:
247
248	integer :: nlev, nw
249
250	real, dimension(nw) :: wl, wc ! lower and central wavelength grid (nm)
251	real, dimension(nlev) :: press, temp, zmmean ! pressure (hpa), temperature (k), molecular mass (g.mol-1)
252
253	! output:
254
255	real, dimension(nlev) :: colinc ! air column increments (molecule.cm-2)
256	real, dimension(nlev,nw) :: dtrl ! rayleigh optical depth
257
258	! local:
259
260	real, parameter :: avo = 6.022e23
261	real, parameter :: g = 3.72
262	real :: dp, nu
263	real, dimension(nw) :: srayl
264	integer :: ilev, iw
265
266	! compute column increments
267
268	do ilev = 1, nlev-1
269	dp = (press(ilev) - press(ilev+1))*100.
270	colinc(ilev) = avo0.1dp/(zmmean(ilev)*g)
271	end do
272	colinc(nlev) = avo0.1press(nlev)100./(zmmean(nlev)g)
273
274	do iw = 1, nw - 1
275
276	! co2 rayleigh cross-section
277	! ityaksov et al., chem. phys. lett., 462, 31-34, 2008
278
279	nu = 1./(wc(iw)*1.e-7)
280	srayl(iw) = 1.78e-26nu*(4. + 0.625)
281	srayl(iw) = srayl(iw)*1.e-20 ! cm2
282
283	do ilev = 1, nlev
284	dtrl(ilev,iw) = colinc(ilev)*srayl(iw) ! cm2
285	end do
286	end do
287
288	end subroutine setair
289
290	!==============================================================================
291
292	subroutine setco2(nd, nlayer, nw, wc, tlay, xsco2_195, xsco2_295,
293	$ xsco2_370, yieldco2, j_co2_o, j_co2_o1d,
294	$ colinc, rm, dt, sj)
295
296	!-----------------------------------------------------------------------------*
297	!= PURPOSE: =*
298	!= Set up the CO2 temperature-dependent cross-sections and optical depth =*
299	!-----------------------------------------------------------------------------*
300
301	implicit none
302
303	! input:
304
305	integer :: nd ! number of photolysis rates
306	integer :: nlayer ! number of vertical layers
307	integer :: nw ! number of wavelength grid points
308	integer :: j_co2_o, j_co2_o1d ! photolysis indexes
309	real, dimension(nw) :: wc ! central wavelength for each interval
310	real, dimension(nw) :: xsco2_195, xsco2_295, xsco2_370 ! co2 cross-sections (cm2)
311	real, dimension(nw) :: yieldco2 ! co2 photodissociation yield
312	real, dimension(nlayer) :: tlay ! temperature (k)
313	real, dimension(nlayer) :: rm ! co2 mixing ratio
314	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
315
316	! output:
317
318	real, dimension(nlayer,nw) :: dt ! optical depth
319	real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2)
320
321	! local:
322
323	integer :: extrapol
324	integer :: i, l
325	real :: temp, sco2
326
327	! extrapol = 0 no extrapolation below 195 k
328	! extrapol = 1 extrapolation below 195 k
329
330	extrapol = 0
331
332	do i = 1, nlayer
333	if (extrapol == 1) then
334	temp = tlay(i)
335	else
336	temp = max(tlay(i), 195.)
337	end if
338	temp = min(temp, 370.)
339	do l = 1, nw-1
340	if (temp <= 295.) then
341	if (xsco2_195(l) /= 0. .and. xsco2_295(l) /= 0.) then
342	sco2 = alog(xsco2_195(l))
343	$ + (alog(xsco2_295(l)) - alog(xsco2_195(l)))
344	$ /(295. - 195.)*(temp - 195.)
345	sco2 = exp(sco2)
346	else
347	sco2 = 0.
348	end if
349	else
350	if (xsco2_295(l) /= 0. .and. xsco2_370(l) /= 0.) then
351	sco2 = alog(xsco2_295(l))
352	$ + (alog(xsco2_370(l)) - alog(xsco2_295(l)))
353	$ /(370. - 295.)*(temp - 295.)
354	sco2 = exp(sco2)
355	else
356	sco2 = 0.
357	end if
358	end if
359
360	! optical depth
361
362	dt(i,l) = colinc(i)rm(i)sco2
363
364	! production of o(1d) for wavelengths shorter than 167 nm
365
366	if (wc(l) >= 167.) then
367	sj(i,l,j_co2_o) = sco2*yieldco2(l)
368	sj(i,l,j_co2_o1d) = 0.
369	else
370	sj(i,l,j_co2_o) = 0.
371	sj(i,l,j_co2_o1d) = sco2*yieldco2(l)
372	end if
373	end do
374	end do
375
376	end subroutine setco2
377
378	!==============================================================================
379
380	subroutine seto2(nd, nlayer, nw, wc, mopt, tlay, xso2_150,
381	$ xso2_200, xso2_250, xso2_300, yieldo2, j_o2_o,
382	$ j_o2_o1d, colinc, rm, dt, sj)
383
384	!-----------------------------------------------------------------------------*
385	!= PURPOSE: =*
386	!= Set up the O2 temperature-dependent cross-sections and optical depth =*
387	!-----------------------------------------------------------------------------*
388
389	implicit none
390
391	! input:
392
393	integer :: nd ! number of photolysis rates
394	integer :: nlayer ! number of vertical layers
395	integer :: nw ! number of wavelength grid points
396	integer :: mopt ! high-res/low-res switch
397	integer :: j_o2_o, j_o2_o1d ! photolysis indexes
398	real, dimension(nw) :: wc ! central wavelength for each interval
399	real, dimension(nw) :: xso2_150, xso2_200, xso2_250, ! o2 cross-sections (cm2)
400	$ xso2_300
401	real, dimension(nw) :: yieldo2 ! o2 photodissociation yield
402	real, dimension(nlayer) :: tlay ! temperature (k)
403	real, dimension(nlayer) :: rm ! o2 mixing ratio
404	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
405
406	! output:
407
408	real, dimension(nlayer,nw) :: dt ! optical depth
409	real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2)
410
411	! local:
412
413	integer :: ilev, iw
414	real :: temp
415	real :: xso2, factor
416
417	! correction by factor if low-resolution in schumann-runge bands
418
419	if (mopt == 1) then
420	factor = 1.
421	else if (mopt == 2) then
422	factor = 0.8
423	end if
424
425	! calculate temperature dependance
426
427	do ilev = 1,nlayer
428	temp = max(tlay(ilev),150.)
429	temp = min(temp, 300.)
430	do iw = 1, nw-1
431	if (tlay(ilev) > 250.) then
432	xso2 = xso2_250(iw) + (xso2_300(iw) - xso2_250(iw))
433	$ /(300. - 250.)*(temp - 250.)
434	else if (tlay(ilev) > 200.) then
435	xso2 = xso2_200(iw) + (xso2_250(iw) - xso2_200(iw))
436	$ /(250. - 200.)*(temp - 200.)
437	else
438	xso2 = xso2_150(iw) + (xso2_200(iw) - xso2_150(iw))
439	$ /(200. - 150.)*(temp - 150.)
440	end if
441
442	if (wc(iw) > 180. .and. wc(iw) < 200.) then
443	xso2 = xso2*factor
444	end if
445
446	! optical depth
447
448	dt(ilev,iw) = colinc(ilev)rm(ilev)xso2
449
450	! production of o(1d) for wavelengths shorter than 175 nm
451
452	if (wc(iw) >= 175.) then
453	sj(ilev,iw,j_o2_o) = xso2*yieldo2(iw)
454	sj(ilev,iw,j_o2_o1d) = 0.
455	else
456	sj(ilev,iw,j_o2_o) = 0.
457	sj(ilev,iw,j_o2_o1d) = xso2*yieldo2(iw)
458	end if
459
460	end do
461	end do
462
463	end subroutine seto2
464
465	!==============================================================================
466
467	subroutine seto3(nd, nlayer, nw, wc, tlay, xso3_218, xso3_298,
468	$ j_o3_o, j_o3_o1d,
469	$ colinc, rm, dt, sj)
470
471	!-----------------------------------------------------------------------------*
472	!= PURPOSE: =*
473	!= Set up the O3 temperature dependent cross-sections and optical depth =*
474	!-----------------------------------------------------------------------------*
475
476	implicit none
477
478	! input:
479
480	integer :: nd ! number of photolysis rates
481	integer :: nlayer ! number of vertical layers
482	integer :: nw ! number of wavelength grid points
483	integer :: j_o3_o, j_o3_o1d ! photolysis indexes
484	real, dimension(nw) :: wc ! central wavelength for each interval
485	real, dimension(nw) :: xso3_218, xso3_298 ! o3 cross-sections (cm2)
486	real, dimension(nlayer) :: tlay ! temperature (k)
487	real, dimension(nlayer) :: rm ! o3 mixing ratio
488	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
489
490	! output:
491
492	real, dimension(nlayer,nw) :: dt ! optical depth
493	real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2)
494
495	! local:
496	!
497	integer :: ilev, iw
498	real :: temp
499	real, dimension(nw) :: xso3(nw)
500	real, dimension(nw) :: qy1d ! quantum yield for o(1d) production
501	real :: q1, q2, a1, a2, a3
502
503	do ilev = 1, nlayer
504	temp = max(tlay(ilev), 218.)
505	temp = min(temp,298.)
506	do iw = 1, nw-1
507	xso3(iw) = xso3_218(iw) + (xso3_298(iw) - xso3_218(iw))
508	$ /(298. - 218.) *(temp - 218.)
509
510	! optical depth
511
512	dt(ilev,iw) = colinc(ilev)rm(ilev)xso3(iw)
513
514	end do
515
516	! calculate quantum yield for o(1d) production (jpl 2006)
517
518	temp = max(tlay(ilev),200.)
519	temp = min(temp,320.)
520	do iw = 1, nw-1
521	if (wc(iw) <= 306.) then
522	qy1d(iw) = 0.90
523	else if (wc(iw) > 306. .and. wc(iw) < 328.) then
524	q1 = 1.
525	q2 = exp(-825.518/(0.695*temp))
526	a1 = (304.225 - wc(iw))/5.576
527	a2 = (314.957 - wc(iw))/6.601
528	a3 = (310.737 - wc(iw))/2.187
529	qy1d(iw) = (q1/(q1 + q2))0.8036exp(-(a1a1a1*a1))
530	$ + (q2/(q1 + q2))8.9061(temp/300.)**2.
531	$ exp(-(a2a2))
532	$ + 0.1192(temp/300.)1.5exp(-(a3*a3))
533	$ + 0.0765
534	else if (wc(iw) >= 328. .and. wc(iw) <= 340.) then
535	qy1d(iw) = 0.08
536	else
537	qy1d(iw) = 0.
538	endif
539	end do
540	do iw = 1, nw-1
541	sj(ilev,iw,j_o3_o) = xso3(iw)*(1. - qy1d(iw))
542	sj(ilev,iw,j_o3_o1d) = xso3(iw)*qy1d(iw)
543	end do
544	end do
545
546	end subroutine seto3
547
548	!==============================================================================
549
550	subroutine seth2o2(nd, nlayer, nw, wc, tlay, xsh2o2, j_h2o2,
551	$ colinc, rm, dt, sj)
552
553	!-----------------------------------------------------------------------------*
554	!= PURPOSE: =*
555	!= Set up the h2o2 temperature dependent cross-sections and optical depth =*
556	!-----------------------------------------------------------------------------*
557
558	implicit none
559
560	! input:
561
562	integer :: nd ! number of photolysis rates
563	integer :: nlayer ! number of vertical layers
564	integer :: nw ! number of wavelength grid points
565	integer :: j_h2o2 ! photolysis index
566	real, dimension(nw) :: wc ! central wavelength for each interval
567	real, dimension(nw) :: xsh2o2 ! h2o2 cross-sections (cm2)
568	real, dimension(nlayer) :: tlay ! temperature (k)
569	real, dimension(nlayer) :: rm ! h2o2 mixing ratio
570	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
571
572	! output:
573
574	real, dimension(nlayer,nw) :: dt ! optical depth
575	real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2)
576
577	! local:
578
579	integer :: ilev, iw
580	real :: a0, a1, a2, a3, a4, a5, a6, a7
581	real :: b0, b1, b2, b3, b4
582	real :: lambda, suma, sumb, chi, temp, xs
583
584	A0 = 6.4761E+04
585	A1 = -9.2170972E+02
586	A2 = 4.535649
587	A3 = -4.4589016E-03
588	A4 = -4.035101E-05
589	A5 = 1.6878206E-07
590	A6 = -2.652014E-10
591	A7 = 1.5534675E-13
592
593	B0 = 6.8123E+03
594	B1 = -5.1351E+01
595	B2 = 1.1522E-01
596	B3 = -3.0493E-05
597	B4 = -1.0924E-07
598
599	! temperature dependance: jpl 2006
600
601	do ilev = 1,nlayer
602	temp = min(max(tlay(ilev),200.),400.)
603	chi = 1./(1. + exp(-1265./temp))
604	do iw = 1, nw-1
605	if ((wc(iw) >= 260.) .and. (wc(iw) < 350.)) then
606	lambda = wc(iw)
607	sumA = ((((((A7lambda + A6)lambda + A5)*lambda +
608	$ A4)lambda +A3)lambda + A2)*lambda +
609	$ A1)*lambda + A0
610	sumB = (((B4lambda + B3)lambda + B2)*lambda +
611	$ B1)*lambda + B0
612	xs = (chisumA + (1. - chi)sumB)*1.e-21
613	sj(ilev,iw,j_h2o2) = xs
614	else
615	sj(ilev,iw,j_h2o2) = xsh2o2(iw)
616	end if
617
618	! optical depth
619
620	dt(ilev,iw) = colinc(ilev)rm(ilev)sj(ilev,iw,j_h2o2)
621	end do
622	end do
623
624	end subroutine seth2o2
625
626	!==============================================================================
627
628	subroutine setno2(nd, nlayer, nw, wc, tlay, xsno2, xsno2_220,
629	$ xsno2_294, yldno2_248, yldno2_298, j_no2,
630	$ colinc, rm, dt, sj)
631
632	!-----------------------------------------------------------------------------*
633	!= PURPOSE: =*
634	!= Set up the no2 temperature-dependent cross-sections and optical depth =*
635	!-----------------------------------------------------------------------------*
636
637	implicit none
638
639	! input:
640
641	integer :: nd ! number of photolysis rates
642	integer :: nlayer ! number of vertical layers
643	integer :: nw ! number of wavelength grid points
644	integer :: j_no2 ! photolysis index
645	real, dimension(nw) :: wc ! central wavelength for each interval
646	real, dimension(nw) :: xsno2, xsno2_220, xsno2_294 ! no2 absorption cross-section at 220-294 k (cm2)
647	real, dimension(nw) :: yldno2_248, yldno2_298 ! no2 quantum yield at 248-298 k
648	real, dimension(nlayer) :: tlay ! temperature (k)
649	real, dimension(nlayer) :: rm ! no2 mixing ratio
650	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
651
652	! output:
653
654	real, dimension(nlayer,nw) :: dt ! optical depth
655	real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2)
656
657	! local:
658
659	integer :: ilev, iw
660	real :: temp, qy
661
662	! temperature dependance: jpl 2006
663
664	do ilev = 1,nlayer
665	temp = max(220.,min(tlay(ilev),294.))
666	do iw = 1, nw - 1
667	if (wc(iw) < 238.) then
668	sj(ilev,iw,j_no2) = xsno2(iw)
669	else
670	sj(ilev,iw,j_no2) = xsno2_220(iw)
671	$ + (xsno2_294(iw) - xsno2_220(iw))
672	$ /(294. - 220.)*(temp - 220.)
673	end if
674
675	! optical depth
676
677	dt(ilev,iw) = colinc(ilev)rm(ilev)sj(ilev,iw,j_no2)
678	end do
679	end do
680
681	! quantum yield: jpl 2006
682
683	do ilev = 1,nlayer
684	temp = max(248.,min(tlay(ilev),298.))
685	do iw = 1, nw - 1
686	qy = yldno2_248(iw) + (yldno2_298(iw) - yldno2_248(iw))
687	$ /(298. - 248.)*(temp - 248.)
688	sj(ilev,iw,j_no2) = sj(ilev,iw,j_no2)*qy
689	end do
690	end do
691
692	end subroutine setno2
693
694	!==============================================================================
695
696	subroutine setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer)
697
698	!-----------------------------------------------------------------------------*
699	!= PURPOSE: =*
700	!= Set aerosol properties for each specified altitude layer. Properties =*
701	!= may be wavelength dependent. =*
702	!-----------------------------------------------------------------------------*
703
704	implicit none
705
706	! input
707
708	integer :: nlayer ! number of vertical layers
709	integer :: nw ! number of wavelength grid points
710	real, dimension(nlayer) :: alt ! altitude (km)
711	real :: tau ! integrated aerosol optical depth at the surface
712
713	! output
714
715	real, dimension(nlayer,nw) :: dtaer ! aerosol optical depth
716	real, dimension(nlayer,nw) :: omaer ! aerosol single scattering albedo
717	real, dimension(nlayer,nw) :: gaer ! aerosol asymmetry parameter
718
719	! local
720
721	integer :: ilay, iw
722	real, dimension(nlayer) :: aer ! dust extinction
723	real :: omega, g, scaleh, gamma
724	real :: dz, tautot, q0
725
726	omega = 0.622 ! single scattering albedo : wolff et al.(2010) at 258 nm
727	g = 0.88 ! asymmetry factor : mateshvili et al. (2007) at 210 nm
728	scaleh = 10. ! scale height (km)
729	gamma = 0.03 ! conrath parameter
730
731	dtaer(:,:) = 0.
732	omaer(:,:) = 0.
733	gaer(:,:) = 0.
734
735	! optical depth profile:
736
737	tautot = 0.
738	do ilay = 1, nlayer-1
739	dz = alt(ilay+1) - alt(ilay)
740	tautot = tautot + exp(gamma(1. - exp(alt(ilay)/scaleh)))dz
741	end do
742
743	q0 = tau/tautot
744	do ilay = 1, nlayer-1
745	dz = alt(ilay+1) - alt(ilay)
746	dtaer(ilay,:) = q0exp(gamma(1. - exp(alt(ilay)/scaleh)))*dz
747	omaer(ilay,:) = omega
748	gaer(ilay,:) = g
749	end do
750
751	end subroutine setaer
752
753	!==============================================================================
754
755	subroutine setcld(nlayer,nw,dtcld,omcld,gcld)
756
757	!-----------------------------------------------------------------------------*
758	!= PURPOSE: =*
759	!= Set cloud properties for each specified altitude layer. Properties =*
760	!= may be wavelength dependent. =*
761	!-----------------------------------------------------------------------------*
762
763	implicit none
764
765	! input
766
767	integer :: nlayer ! number of vertical layers
768	integer :: nw ! number of wavelength grid points
769
770	! output
771
772	real, dimension(nlayer,nw) :: dtcld ! cloud optical depth
773	real, dimension(nlayer,nw) :: omcld ! cloud single scattering albedo
774	real, dimension(nlayer,nw) :: gcld ! cloud asymmetry parameter
775
776	! local
777
778	integer :: ilay, iw
779
780	! dtcld : optical depth
781	! omcld : single scattering albedo
782	! gcld : asymmetry factor
783
784	do ilay = 1, nlayer - 1
785	do iw = 1, nw - 1
786	dtcld(ilay,iw) = 0. ! no clouds for the moment
787	omcld(ilay,iw) = 0.99
788	gcld(ilay,iw) = 0.85
789	end do
790	end do
791
792	end subroutine setcld
793
794	!==============================================================================
795
796	subroutine sphers(nlev, z, zen, dsdh, nid)
797
798	!-----------------------------------------------------------------------------*
799	!= PURPOSE: =*
800	!= Calculate slant path over vertical depth ds/dh in spherical geometry. =*
801	!= Calculation is based on: A.Dahlback, and K.Stamnes, A new spheric model =*
802	!= for computing the radiation field available for photolysis and heating =*
803	!= at twilight, Planet.Space Sci., v39, n5, pp. 671-683, 1991 (Appendix B) =*
804	!-----------------------------------------------------------------------------*
805	!= PARAMETERS: =*
806	!= NZ - INTEGER, number of specified altitude levels in the working (I)=*
807	!= grid =*
808	!= Z - REAL, specified altitude working grid (km) (I)=*
809	!= ZEN - REAL, solar zenith angle (degrees) (I)=*
810	!= DSDH - REAL, slant path of direct beam through each layer crossed (O)=*
811	!= when travelling from the top of the atmosphere to layer i; =*
812	!= DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 =*
813	!= NID - INTEGER, number of layers crossed by the direct beam when (O)=*
814	!= travelling from the top of the atmosphere to layer i; =*
815	!= NID(i), i = 0..NZ-1 =*
816	!-----------------------------------------------------------------------------*
817
818	implicit none
819
820	! input
821
822	integer, intent(in) :: nlev
823	real, dimension(nlev), intent(in) :: z
824	real, intent(in) :: zen
825
826	! output
827
828	INTEGER nid(0:nlev)
829	REAL dsdh(0:nlev,nlev)
830
831	! more program constants
832
833	REAL re, ze(nlev)
834	REAL dr
835	real radius
836	parameter (radius = 3393.)
837
838	! local
839
840	real :: pi, zenrad, rpsinz, rj, rjp1, dsj, dhj, ga, gb, sm
841	integer :: i, j, k, id, nlay
842
843	REAL zd(0:nlev-1)
844
845	!-----------------------------------------------------------------------------
846
847	pi = acos(-1.0)
848	dr = pi/180.
849	zenrad = zen*dr
850
851	! number of layers:
852
853	nlay = nlev - 1
854
855	! include the elevation above sea level to the radius of Mars:
856
857	re = radius + z(1)
858
859	! correspondingly z changed to the elevation above Mars surface:
860
861	DO k = 1, nlev
862	ze(k) = z(k) - z(1)
863	END DO
864
865	! inverse coordinate of z
866
867	zd(0) = ze(nlev)
868	DO k = 1, nlay
869	zd(k) = ze(nlev - k)
870	END DO
871
872	! initialise dsdh(i,j), nid(i)
873
874	nid(:) = 0.
875	dsdh(:,:) = 0.
876
877	! calculate ds/dh of every layer
878
879	do i = 0,nlay
880	rpsinz = (re + zd(i))*sin(zenrad)
881
882	IF ( (zen .GT. 90.0) .AND. (rpsinz .LT. re) ) THEN
883	nid(i) = -1
884	ELSE
885
886	! Find index of layer in which the screening height lies
887
888	id = i
889	if (zen > 90.) then
890	do j = 1,nlay
891	IF( (rpsinz .LT. ( zd(j-1) + re ) ) .AND.
892	$ (rpsinz .GE. ( zd(j) + re )) ) id = j
893	end do
894	end if
895
896	do j = 1,id
897	sm = 1.0
898	IF (j .EQ. id .AND. id .EQ. i .AND. zen .GT. 90.0)
899	$ sm = -1.0
900
901	rj = re + zd(j-1)
902	rjp1 = re + zd(j)
903
904	dhj = zd(j-1) - zd(j)
905
906	ga = rjrj - rpsinzrpsinz
907	gb = rjp1rjp1 - rpsinzrpsinz
908
909	ga = max(ga, 0.)
910	gb = max(gb, 0.)
911
912	IF (id.GT.i .AND. j.EQ.id) THEN
913	dsj = sqrt(ga)
914	ELSE
915	dsj = sqrt(ga) - sm*sqrt(gb)
916	END IF
917	dsdh(i,j) = dsj/dhj
918	end do
919	nid(i) = id
920	end if
921	end do ! i = 0,nlay
922
923	end subroutine sphers
924
925	!==============================================================================
926
927	SUBROUTINE rtlink(nlev, nw, iw, ag, zen, dsdh, nid, dtrl,
928	$ dagas, dtcld, omcld, gcld, dtaer, omaer, gaer,
929	$ edir, edn, eup, fdir, fdn, fup)
930
931	implicit none
932
933	! input
934
935	integer, intent(in) :: nlev, nw, iw ! number of wavelength grid points
936	REAL ag
937	REAL zen
938	REAL dsdh(0:nlev,nlev)
939	INTEGER nid(0:nlev)
940
941	REAL dtrl(nlev,nw)
942	REAL dagas(nlev,nw)
943	REAL dtcld(nlev,nw), omcld(nlev,nw), gcld(nlev,nw)
944	REAL dtaer(nlev,nw), omaer(nlev,nw), gaer(nlev,nw)
945
946	! output
947
948	REAL edir(nlev), edn(nlev), eup(nlev)
949	REAL fdir(nlev), fdn(nlev), fup(nlev)
950
951	! local:
952
953	REAL dt(nlev), om(nlev), g(nlev)
954	REAL dtabs,dtsct,dscld,dsaer,dacld,daaer
955	INTEGER i, ii
956	real, parameter :: largest = 1.e+36
957
958	! specific two ps2str
959
960	REAL ediri(nlev), edni(nlev), eupi(nlev)
961	REAL fdiri(nlev), fdni(nlev), fupi(nlev)
962
963	logical, save :: delta = .true.
964
965	!_______________________________________________________________________
966
967	! initialize:
968
969	do i = 1, nlev
970	fdir(i) = 0.
971	fup(i) = 0.
972	fdn(i) = 0.
973	edir(i) = 0.
974	eup(i) = 0.
975	edn(i) = 0.
976	end do
977
978	do i = 1, nlev - 1
979	dscld = dtcld(i,iw)*omcld(i,iw)
980	dacld = dtcld(i,iw)*(1.-omcld(i,iw))
981
982	dsaer = dtaer(i,iw)*omaer(i,iw)
983	daaer = dtaer(i,iw)*(1.-omaer(i,iw))
984
985	dtsct = dtrl(i,iw) + dscld + dsaer
986	dtabs = dagas(i,iw) + dacld + daaer
987
988	dtabs = amax1(dtabs,1./largest)
989	dtsct = amax1(dtsct,1./largest)
990
991	! invert z-coordinate:
992
993	ii = nlev - i
994	dt(ii) = dtsct + dtabs
995	om(ii) = dtsct/(dtsct + dtabs)
996	IF(dtsct .EQ. 1./largest) om(ii) = 1./largest
997	g(ii) = (gcld(i,iw)*dscld +
998	$ gaer(i,iw)*dsaer)/dtsct
999	end do
1000
1001	! call rt routine:
1002
1003	call ps2str(nlev, zen, ag, dt, om, g,
1004	$ dsdh, nid, delta,
1005	$ fdiri, fupi, fdni, ediri, eupi, edni)
1006
1007	! output (invert z-coordinate)
1008
1009	do i = 1, nlev
1010	ii = nlev - i + 1
1011	fdir(i) = fdiri(ii)
1012	fup(i) = fupi(ii)
1013	fdn(i) = fdni(ii)
1014	edir(i) = ediri(ii)
1015	eup(i) = eupi(ii)
1016	edn(i) = edni(ii)
1017	end do
1018
1019	end subroutine rtlink
1020
1021	=============================================================================
1022
1023	subroutine ps2str(nlev,zen,rsfc,tauu,omu,gu,
1024	$ dsdh, nid, delta,
1025	$ fdr, fup, fdn, edr, eup, edn)
1026
1027	!-----------------------------------------------------------------------------*
1028	!= PURPOSE: =*
1029	!= Solve two-stream equations for multiple layers. The subroutine is based =*
1030	!= on equations from: Toon et al., J.Geophys.Res., v94 (D13), Nov 20, 1989.=*
1031	!= It contains 9 two-stream methods to choose from. A pseudo-spherical =*
1032	!= correction has also been added. =*
1033	!-----------------------------------------------------------------------------*
1034	!= PARAMETERS: =*
1035	!= NLEVEL - INTEGER, number of specified altitude levels in the working (I)=*
1036	!= grid =*
1037	!= ZEN - REAL, solar zenith angle (degrees) (I)=*
1038	!= RSFC - REAL, surface albedo at current wavelength (I)=*
1039	!= TAUU - REAL, unscaled optical depth of each layer (I)=*
1040	!= OMU - REAL, unscaled single scattering albedo of each layer (I)=*
1041	!= GU - REAL, unscaled asymmetry parameter of each layer (I)=*
1042	!= DSDH - REAL, slant path of direct beam through each layer crossed (I)=*
1043	!= when travelling from the top of the atmosphere to layer i; =*
1044	!= DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 =*
1045	!= NID - INTEGER, number of layers crossed by the direct beam when (I)=*
1046	!= travelling from the top of the atmosphere to layer i; =*
1047	!= NID(i), i = 0..NZ-1 =*
1048	!= DELTA - LOGICAL, switch to use delta-scaling (I)=*
1049	!= .TRUE. -> apply delta-scaling =*
1050	!= .FALSE.-> do not apply delta-scaling =*
1051	!= FDR - REAL, contribution of the direct component to the total (O)=*
1052	!= actinic flux at each altitude level =*
1053	!= FUP - REAL, contribution of the diffuse upwelling component to (O)=*
1054	!= the total actinic flux at each altitude level =*
1055	!= FDN - REAL, contribution of the diffuse downwelling component to (O)=*
1056	!= the total actinic flux at each altitude level =*
1057	!= EDR - REAL, contribution of the direct component to the total (O)=*
1058	!= spectral irradiance at each altitude level =*
1059	!= EUP - REAL, contribution of the diffuse upwelling component to (O)=*
1060	!= the total spectral irradiance at each altitude level =*
1061	!= EDN - REAL, contribution of the diffuse downwelling component to (O)=*
1062	= the total spectral irradiance at each altitude level =
1063	!-----------------------------------------------------------------------------*
1064
1065	implicit none
1066
1067	! input:
1068
1069	INTEGER nlev
1070	REAL zen, rsfc
1071	REAL tauu(nlev), omu(nlev), gu(nlev)
1072	REAL dsdh(0:nlev,nlev)
1073	INTEGER nid(0:nlev)
1074	LOGICAL delta
1075
1076	! output:
1077
1078	REAL fup(nlev),fdn(nlev),fdr(nlev)
1079	REAL eup(nlev),edn(nlev),edr(nlev)
1080
1081	! local:
1082
1083	REAL tausla(0:nlev), tauc(0:nlev)
1084	REAL mu2(0:nlev), mu, sum
1085
1086	! internal coefficients and matrix
1087
1088	REAL lam(nlev),taun(nlev),bgam(nlev)
1089	REAL e1(nlev),e2(nlev),e3(nlev),e4(nlev)
1090	REAL cup(nlev),cdn(nlev),cuptn(nlev),cdntn(nlev)
1091	REAL mu1(nlev)
1092	INTEGER row
1093	REAL a(2nlev),b(2nlev),d(2nlev),e(2nlev),y(2*nlev)
1094
1095	! other:
1096
1097	REAL pifs, fdn0
1098	REAL gi(nlev), omi(nlev), tempg
1099	REAL f, g, om
1100	REAL gam1, gam2, gam3, gam4
1101	real, parameter :: largest = 1.e+36
1102	real, parameter :: precis = 1.e-7
1103
1104	! For calculations of Associated Legendre Polynomials for GAMA1,2,3,4
1105	! in delta-function, modified quadrature, hemispheric constant,
1106	! Hybrid modified Eddington-delta function metods, p633,Table1.
1107	! W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
1108	! W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440,
1109	! uncomment the following two lines and the appropriate statements further
1110	! down.
1111	! REAL YLM0, YLM2, YLM4, YLM6, YLM8, YLM10, YLM12, YLMS, BETA0,
1112	! > BETA1, BETAn, amu1, subd
1113
1114	REAL expon, expon0, expon1, divisr, temp, up, dn
1115	REAL ssfc
1116	INTEGER nlayer, mrows, lev
1117
1118	INTEGER i, j
1119
1120	! Some additional program constants:
1121
1122	real pi, dr
1123	REAL eps
1124	PARAMETER (eps = 1.E-3)
1125	!_______________________________________________________________________
1126
1127	! MU = cosine of solar zenith angle
1128	! RSFC = surface albedo
1129	! TAUU = unscaled optical depth of each layer
1130	! OMU = unscaled single scattering albedo
1131	! GU = unscaled asymmetry factor
1132	! KLEV = max dimension of number of layers in atmosphere
1133	! NLAYER = number of layers in the atmosphere
1134	! NLEVEL = nlayer + 1 = number of levels
1135
1136	! initial conditions: pi*solar flux = 1; diffuse incidence = 0
1137
1138	pifs = 1.
1139	fdn0 = 0.
1140
1141	nlayer = nlev - 1
1142
1143	pi = acos(-1.)
1144	dr = pi/180.
1145	mu = COS(zen*dr)
1146
1147	!************* compute coefficients for each layer:
1148	! GAM1 - GAM4 = 2-stream coefficients, different for different approximations
1149	! EXPON0 = calculation of e when TAU is zero
1150	! EXPON1 = calculation of e when TAU is TAUN
1151	! CUP and CDN = calculation when TAU is zero
1152	! CUPTN and CDNTN = calc. when TAU is TAUN
1153	! DIVISR = prevents division by zero
1154
1155	do j = 0, nlev
1156	tauc(j) = 0.
1157	tausla(j) = 0.
1158	mu2(j) = 1./SQRT(largest)
1159	end do
1160
1161	IF (.NOT. delta) THEN
1162	DO i = 1, nlayer
1163	gi(i) = gu(i)
1164	omi(i) = omu(i)
1165	taun(i) = tauu(i)
1166	END DO
1167	ELSE
1168
1169	! delta-scaling. Have to be done for delta-Eddington approximation,
1170	! delta discrete ordinate, Practical Improved Flux Method, delta function,
1171	! and Hybrid modified Eddington-delta function methods approximations
1172
1173	DO i = 1, nlayer
1174	f = gu(i)*gu(i)
1175	gi(i) = (gu(i) - f)/(1 - f)
1176	omi(i) = (1 - f)omu(i)/(1 - omu(i)f)
1177	taun(i) = (1 - omu(i)f)tauu(i)
1178	END DO
1179	END IF
1180
1181	! calculate slant optical depth at the top of the atmosphere when zen>90.
1182	! in this case, higher altitude of the top layer is recommended.
1183
1184	IF (zen .GT. 90.0) THEN
1185	IF (nid(0) .LT. 0) THEN
1186	tausla(0) = largest
1187	ELSE
1188	sum = 0.0
1189	DO j = 1, nid(0)
1190	sum = sum + 2.taun(j)dsdh(0,j)
1191	END DO
1192	tausla(0) = sum
1193	END IF
1194	END IF
1195
1196	DO 11, i = 1, nlayer
1197	g = gi(i)
1198	om = omi(i)
1199	tauc(i) = tauc(i-1) + taun(i)
1200
1201	! stay away from 1 by precision. For g, also stay away from -1
1202
1203	tempg = AMIN1(abs(g),1. - precis)
1204	g = SIGN(tempg,g)
1205	om = AMIN1(om,1.-precis)
1206
1207	! calculate slant optical depth
1208
1209	IF (nid(i) .LT. 0) THEN
1210	tausla(i) = largest
1211	ELSE
1212	sum = 0.0
1213	DO j = 1, MIN(nid(i),i)
1214	sum = sum + taun(j)*dsdh(i,j)
1215	END DO
1216	DO j = MIN(nid(i),i)+1,nid(i)
1217	sum = sum + 2.taun(j)dsdh(i,j)
1218	END DO
1219	tausla(i) = sum
1220	IF (tausla(i) .EQ. tausla(i-1)) THEN
1221	mu2(i) = SQRT(largest)
1222	ELSE
1223	mu2(i) = (tauc(i)-tauc(i-1))/(tausla(i)-tausla(i-1))
1224	mu2(i) = SIGN( AMAX1(ABS(mu2(i)),1./SQRT(largest)),
1225	$ mu2(i) )
1226	END IF
1227	END IF
1228
1229	!** the following gamma equations are from pg 16,289, Table 1
1230	!** save mu1 for each approx. for use in converting irradiance to actinic flux
1231
1232	! Eddington approximation(Joseph et al., 1976, JAS, 33, 2452):
1233
1234	c gam1 = (7. - om(4. + 3.g))/4.
1235	c gam2 = -(1. - om(4. - 3.g))/4.
1236	c gam3 = (2. - 3.gmu)/4.
1237	c gam4 = 1. - gam3
1238	c mu1(i) = 0.5
1239
1240	* quadrature (Liou, 1973, JAS, 30, 1303-1326; 1974, JAS, 31, 1473-1475):
1241
1242	c gam1 = 1.7320508(2. - om(1. + g))/2.
1243	c gam2 = 1.7320508om(1. - g)/2.
1244	c gam3 = (1. - 1.7320508gmu)/2.
1245	c gam4 = 1. - gam3
1246	c mu1(i) = 1./sqrt(3.)
1247
1248	* hemispheric mean (Toon et al., 1089, JGR, 94, 16287):
1249
1250	gam1 = 2. - om*(1. + g)
1251	gam2 = om*(1. - g)
1252	gam3 = (2. - g*mu)/4.
1253	gam4 = 1. - gam3
1254	mu1(i) = 0.5
1255
1256	* PIFM (Zdunkovski et al.,1980, Conrib.Atmos.Phys., 53, 147-166):
1257	c GAM1 = 0.25(8. - OM(5. + 3.*G))
1258	c GAM2 = 0.75OM(1.-G)
1259	c GAM3 = 0.25(2.-3.G*MU)
1260	c GAM4 = 1. - GAM3
1261	c mu1(i) = 0.5
1262
1263	* delta discrete ordinates (Schaller, 1979, Contrib.Atmos.Phys, 52, 17-26):
1264	c GAM1 = 0.51.7320508(2. - OM*(1. + G))
1265	c GAM2 = 0.51.7320508OM*(1.-G)
1266	c GAM3 = 0.5(1.-1.7320508G*MU)
1267	c GAM4 = 1. - GAM3
1268	c mu1(i) = 1./sqrt(3.)
1269
1270	* Calculations of Associated Legendre Polynomials for GAMA1,2,3,4
1271	* in delta-function, modified quadrature, hemispheric constant,
1272	* Hybrid modified Eddington-delta function metods, p633,Table1.
1273	* W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
1274	* W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440
1275	c YLM0 = 2.
1276	c YLM2 = -3.GMU
1277	c YLM4 = 0.875G3MU(5.MU**2-3.)
1278	c YLM6=-0.171875G5MU(15.-70.MU*2+63.MU**4)
1279	c YLM8=+0.073242G7MU(-35.+315.MU*2-693.MU**4
1280	c +429.MU**6)
1281	c YLM10=-0.008118G9MU(315.-4620.MU*2+18018.MU**4
1282	c -25740.MU*6+12155.MU**8)
1283	c YLM12=0.003685G11MU(-693.+15015.MU*2-90090.MU**4
1284	c +218790.MU*6-230945.MU*8+88179.MU**10)
1285	c YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
1286	c YLMS=0.25*YLMS
1287	c BETA0 = YLMS
1288	c
1289	c amu1=1./1.7320508
1290	c YLM0 = 2.
1291	c YLM2 = -3.Gamu1
1292	c YLM4 = 0.875G3amu1(5.amu1**2-3.)
1293	c YLM6=-0.171875G5amu1(15.-70.amu1*2+63.amu1**4)
1294	c YLM8=+0.073242G7amu1(-35.+315.amu1*2-693.amu1**4
1295	c +429.amu1**6)
1296	c YLM10=-0.008118G9amu1(315.-4620.amu1*2+18018.amu1**4
1297	c -25740.amu1*6+12155.amu1**8)
1298	c YLM12=0.003685G11amu1(-693.+15015.amu1*2-90090.amu1**4
1299	c +218790.amu1*6-230945.amu1*8+88179.amu1**10)
1300	c YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
1301	c YLMS=0.25*YLMS
1302	c BETA1 = YLMS
1303	c
1304	c BETAn = 0.25(2. - 1.5G-0.21875G3-0.085938G**5
1305	c -0.045776G**7)
1306
1307
1308	* Hybrid modified Eddington-delta function(Meador and Weaver,1980,JAS,37,630):
1309	c subd=4.(1.-GG*(1.-MU))
1310	c GAM1 = (7.-3.GG-OM(4.+3.G)+OMGG(4.BETA0+3.*G))/subd
1311	c GAM2 =-(1.-GG-OM(4.-3.G)-OMGG(4.BETA0+3.G-4.))/subd
1312	c GAM3 = BETA0
1313	c GAM4 = 1. - GAM3
1314	c mu1(i) = (1. - gg(1.- mu) )/(2. - g*g)
1315
1316	*****
1317	* delta function (Meador, and Weaver, 1980, JAS, 37, 630):
1318	c GAM1 = (1. - OM*(1. - beta0))/MU
1319	c GAM2 = OM*BETA0/MU
1320	c GAM3 = BETA0
1321	c GAM4 = 1. - GAM3
1322	c mu1(i) = mu
1323	*****
1324	* modified quadrature (Meador, and Weaver, 1980, JAS, 37, 630):
1325	c GAM1 = 1.7320508(1. - OM(1. - beta1))
1326	c GAM2 = 1.7320508OMbeta1
1327	c GAM3 = BETA0
1328	c GAM4 = 1. - GAM3
1329	c mu1(i) = 1./sqrt(3.)
1330
1331	* hemispheric constant (Toon et al., 1989, JGR, 94, 16287):
1332	c GAM1 = 2.(1. - OM(1. - betan))
1333	c GAM2 = 2.OMBETAn
1334	c GAM3 = BETA0
1335	c GAM4 = 1. - GAM3
1336	c mu1(i) = 0.5
1337
1338	*****
1339
1340	* lambda = pg 16,290 equation 21
1341	* big gamma = pg 16,290 equation 22
1342	* if gam2 = 0., then bgam = 0.
1343
1344	lam(i) = sqrt(gam1gam1 - gam2gam2)
1345
1346	IF (gam2 .NE. 0.) THEN
1347	bgam(i) = (gam1 - lam(i))/gam2
1348	ELSE
1349	bgam(i) = 0.
1350	END IF
1351
1352	expon = EXP(-lam(i)*taun(i))
1353
1354	* e1 - e4 = pg 16,292 equation 44
1355
1356	e1(i) = 1. + bgam(i)*expon
1357	e2(i) = 1. - bgam(i)*expon
1358	e3(i) = bgam(i) + expon
1359	e4(i) = bgam(i) - expon
1360
1361	* the following sets up for the C equations 23, and 24
1362	* found on page 16,290
1363	* prevent division by zero (if LAMBDA=1/MU, shift 1/MU^2 by EPS = 1.E-3
1364	* which is approx equiv to shifting MU by 0.5EPS (MU)**3
1365
1366	expon0 = EXP(-tausla(i-1))
1367	expon1 = EXP(-tausla(i))
1368
1369	divisr = lam(i)lam(i) - 1./(mu2(i)mu2(i))
1370	temp = AMAX1(eps,abs(divisr))
1371	divisr = SIGN(temp,divisr)
1372
1373	up = ompifs((gam1 - 1./mu2(i))gam3 + gam4gam2)/divisr
1374	dn = ompifs((gam1 + 1./mu2(i))gam4 + gam2gam3)/divisr
1375
1376	* cup and cdn are when tau is equal to zero
1377	* cuptn and cdntn are when tau is equal to taun
1378
1379	cup(i) = up*expon0
1380	cdn(i) = dn*expon0
1381	cuptn(i) = up*expon1
1382	cdntn(i) = dn*expon1
1383
1384	11 CONTINUE
1385
1386	*************** set up matrix ****
1387	* ssfc = pg 16,292 equation 37 where pi Fs is one (unity).
1388
1389	ssfc = rsfcmuEXP(-tausla(nlayer))*pifs
1390
1391	* MROWS = the number of rows in the matrix
1392
1393	mrows = 2*nlayer
1394
1395	* the following are from pg 16,292 equations 39 - 43.
1396	* set up first row of matrix:
1397
1398	i = 1
1399	a(1) = 0.
1400	b(1) = e1(i)
1401	d(1) = -e2(i)
1402	e(1) = fdn0 - cdn(i)
1403
1404	row=1
1405
1406	* set up odd rows 3 thru (MROWS - 1):
1407
1408	i = 0
1409	DO 20, row = 3, mrows - 1, 2
1410	i = i + 1
1411	a(row) = e2(i)e3(i) - e4(i)e1(i)
1412	b(row) = e1(i)e1(i + 1) - e3(i)e3(i + 1)
1413	d(row) = e3(i)e4(i + 1) - e1(i)e2(i + 1)
1414	e(row) = e3(i)*(cup(i + 1) - cuptn(i)) +
1415	$ e1(i)*(cdntn(i) - cdn(i + 1))
1416	20 CONTINUE
1417
1418	* set up even rows 2 thru (MROWS - 2):
1419
1420	i = 0
1421	DO 30, row = 2, mrows - 2, 2
1422	i = i + 1
1423	a(row) = e2(i + 1)e1(i) - e3(i)e4(i + 1)
1424	b(row) = e2(i)e2(i + 1) - e4(i)e4(i + 1)
1425	d(row) = e1(i + 1)e4(i + 1) - e2(i + 1)e3(i + 1)
1426	e(row) = (cup(i + 1) - cuptn(i))*e2(i + 1) -
1427	$ (cdn(i + 1) - cdntn(i))*e4(i + 1)
1428	30 CONTINUE
1429
1430	* set up last row of matrix at MROWS:
1431
1432	row = mrows
1433	i = nlayer
1434
1435	a(row) = e1(i) - rsfc*e3(i)
1436	b(row) = e2(i) - rsfc*e4(i)
1437	d(row) = 0.
1438	e(row) = ssfc - cuptn(i) + rsfc*cdntn(i)
1439
1440	* solve tri-diagonal matrix:
1441
1442	CALL tridiag(a, b, d, e, y, mrows)
1443
1444	**** unfold solution of matrix, compute output fluxes:
1445
1446	row = 1
1447	lev = 1
1448	j = 1
1449
1450	* the following equations are from pg 16,291 equations 31 & 32
1451
1452	fdr(lev) = EXP( -tausla(0) )
1453	edr(lev) = mu * fdr(lev)
1454	edn(lev) = fdn0
1455	eup(lev) = y(row)e3(j) - y(row + 1)e4(j) + cup(j)
1456	fdn(lev) = edn(lev)/mu1(lev)
1457	fup(lev) = eup(lev)/mu1(lev)
1458
1459	DO 60, lev = 2, nlayer + 1
1460	fdr(lev) = EXP(-tausla(lev-1))
1461	edr(lev) = mu *fdr(lev)
1462	edn(lev) = y(row)e3(j) + y(row + 1)e4(j) + cdntn(j)
1463	eup(lev) = y(row)e1(j) + y(row + 1)e2(j) + cuptn(j)
1464	fdn(lev) = edn(lev)/mu1(j)
1465	fup(lev) = eup(lev)/mu1(j)
1466
1467	row = row + 2
1468	j = j + 1
1469	60 CONTINUE
1470
1471	end subroutine ps2str
1472
1473	=============================================================================
1474
1475	subroutine tridiag(a,b,c,r,u,n)
1476
1477	!_______________________________________________________________________
1478	! solves tridiagonal system. From Numerical Recipies, p. 40
1479	!_______________________________________________________________________
1480
1481	IMPLICIT NONE
1482
1483	! input:
1484
1485	INTEGER n
1486	REAL a, b, c, r
1487	DIMENSION a(n),b(n),c(n),r(n)
1488
1489	! output:
1490
1491	REAL u
1492	DIMENSION u(n)
1493
1494	! local:
1495
1496	INTEGER j
1497
1498	REAL bet, gam
1499	DIMENSION gam(n)
1500	!_______________________________________________________________________
1501
1502	IF (b(1) .EQ. 0.) STOP 1001
1503	bet = b(1)
1504	u(1) = r(1)/bet
1505	DO 11, j = 2, n
1506	gam(j) = c(j - 1)/bet
1507	bet = b(j) - a(j)*gam(j)
1508	IF (bet .EQ. 0.) STOP 2002
1509	u(j) = (r(j) - a(j)*u(j - 1))/bet
1510	11 CONTINUE
1511	DO 12, j = n - 1, 1, -1
1512	u(j) = u(j) - gam(j + 1)*u(j + 1)
1513	12 CONTINUE
1514	!_______________________________________________________________________
1515
1516	end subroutine tridiag
1517
1518	end subroutine photolysis_online

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

Download in other formats:

Original Format