Context Navigation

photolysis_online.F @ 2325

Last change on this file since 2325 was 2170, checked in by flefevre, 5 years ago
la chimie des ions est désormais optionnelle (FGG) regroupement des options de la photochimie dans calchim_mod (FL) mise à jour cinétique CO + OH, inactive pour le moment (FL)
File size: 51.0 KB

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

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

Download in other formats:

Original Format