Context Navigation

photolysis_mod.F90 @ 2041

Last change on this file since 2041 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: 56.1 KB

Line
1	module photolysis_mod
2
3	implicit none
4
5	integer, parameter :: nphot = 13 ! number of photolysis
6	integer, parameter :: nabs = 10 ! number of absorbing gases
7
8	! number of reactions in chemical solver
9
10	integer, parameter :: nb_phot_max = nphot + 9 ! photolysis + quenching/heterogeneous
11	integer, parameter :: nb_reaction_3_max = 6 ! quadratic
12	integer, parameter :: nb_reaction_4_max = 31 ! bimolecular
13
14	! spectral grid
15
16	integer, parameter :: nw = 162 ! number of spectral intervals (low-res)
17	integer, save :: mopt ! high-res/low-res switch
18
19	real, dimension(nw), save :: wl, wc, wu ! lower, center, upper wavelength for each interval
20
21	! solar flux
22
23	real, dimension(nw), save :: f ! solar flux (w.m-2.nm-1) at 1 au
24
25	! cross-sections and quantum yields
26
27	real, dimension(nw), save :: xsco2_195, xsco2_295, xsco2_370 ! co2 absorption cross-section at 195-295-370 k (cm2)
28	real, dimension(nw), save :: yieldco2 ! co2 photodissociation yield
29	real, dimension(nw), save :: xso2_150, xso2_200, xso2_250, xso2_300 ! o2 absorption cross-section at 150-200-250-300 k (cm2)
30	real, dimension(nw), save :: yieldo2 ! o2 photodissociation yield
31	real, dimension(nw), save :: xso3_218, xso3_298 ! o3 absorption cross-section at 218-298 k (cm2)
32	real, dimension(nw), save :: xsh2o ! h2o absorption cross-section (cm2)
33	real, dimension(nw), save :: xsh2o2 ! h2o2 absorption cross-section (cm2)
34	real, dimension(nw), save :: xsho2 ! ho2 absorption cross-section (cm2)
35	real, dimension(nw), save :: xsh2 ! h2 absorption cross-section (cm2)
36	real, dimension(nw), save :: yieldh2 ! h2 photodissociation yield
37	real, dimension(nw), save :: xsno2, xsno2_220, xsno2_294 ! no2 absorption cross-section at 220-294 k (cm2)
38	real, dimension(nw), save :: yldno2_248, yldno2_298 ! no2 quantum yield at 248-298 k
39	real, dimension(nw), save :: xsno ! no absorption cross-section (cm2)
40	real, dimension(nw), save :: yieldno ! no photodissociation yield
41	real, dimension(nw), save :: xsn2 ! n2 absorption cross-section (cm2)
42	real, dimension(nw), save :: yieldn2 ! n2 photodissociation yield
43	real, dimension(nw), save :: albedo ! surface albedo
44
45	contains
46
47	subroutine init_photolysis
48
49	! initialise on-line photolysis
50
51	! mopt = 1 high-resolution
52	! mopt = 2 low-resolution (recommended for gcm use)
53
54	mopt = 2
55
56	! set wavelength grid
57
58	call gridw(nw,wl,wc,wu,mopt)
59
60	! read and grid solar flux data
61
62	call rdsolarflux(nw,wl,wc,f)
63
64	! read and grid o2 cross-sections
65
66	call rdxso2(nw,wl,xso2_150,xso2_200,xso2_250,xso2_300,yieldo2)
67
68	! read and grid co2 cross-sections
69
70	call rdxsco2(nw,wl,xsco2_195,xsco2_295,xsco2_370,yieldco2)
71
72	! read and grid o3 cross-sections
73
74	call rdxso3(nw,wl,xso3_218,xso3_298)
75
76	! read and grid h2o cross-sections
77
78	call rdxsh2o(nw,wl,xsh2o)
79
80	! read and grid h2o2 cross-sections
81
82	call rdxsh2o2(nw,wl,xsh2o2)
83
84	! read and grid ho2 cross-sections
85
86	call rdxsho2(nw,wl,xsho2)
87
88	! read and grid h2 cross-sections
89
90	call rdxsh2(nw,wl,wc,xsh2,yieldh2)
91
92	! read and grid no cross-sections
93
94	call rdxsno(nw,wl,xsno,yieldno)
95
96	! read and grid no2 cross-sections
97
98	call rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294,yldno2_248,yldno2_298)
99
100	! read and grid n2 cross-sections
101
102	call rdxsn2(nw,wl,xsn2,yieldn2)
103
104	! set surface albedo
105
106	call setalb(nw,wl,albedo)
107
108	end subroutine init_photolysis
109
110	!==============================================================================
111
112	subroutine gridw(nw,wl,wc,wu,mopt)
113
114	! Create the wavelength grid for all interpolations and radiative transfer
115	! calculations. Grid may be irregularly spaced. Wavelengths are in nm.
116	! No gaps are allowed within the wavelength grid.
117
118	implicit none
119
120	! input
121
122	integer :: nw ! number of wavelength grid points
123	integer :: mopt ! high-res/low-res switch
124
125	! output
126
127	real, dimension(nw) :: wl, wc, wu ! lower, center, upper wavelength for each interval
128
129	! local
130
131	real :: wincr ! wavelength increment
132	integer :: iw, kw
133
134	! mopt = 1 high-resolution mode (3789 intervals)
135	!
136	! 0-108 nm : 1.0 nm
137	! 108-124 nm : 0.1 nm
138	! 124-175 nm : 0.5 nm
139	! 175-205 nm : 0.01 nm
140	! 205-365 nm : 0.5 nm
141	! 365-850 nm : 5.0 nm
142	!
143	! mopt = 2 low-resolution mode
144	!
145	! 0-60 nm : 6.0 nm
146	! 60-80 nm : 2.0 nm
147	! 80-85 nm : 5.0 nm
148	! 85-117 nm : 2.0 nm
149	! 117-120 nm : 5.0 nm
150	! 120-123 nm : 0.2 nm
151	! 123-163 nm : 5.0 nm
152	! 163-175 nm : 2.0 nm
153	! 175-205 nm : 0.5 nm
154	! 205-245 nm : 5.0 nm
155	! 245-415 nm : 10.0 nm
156	! 415-815 nm : 50.0 nm
157
158	if (mopt == 1) then ! high-res
159
160	! define wavelength intervals of width 1.0 nm from 0 to 108 nm:
161
162	kw = 0
163	wincr = 1.0
164	do iw = 0, 107
165	kw = kw + 1
166	wl(kw) = real(iw)
167	wu(kw) = wl(kw) + wincr
168	wc(kw) = (wl(kw) + wu(kw))/2.
169	end do
170
171	! define wavelength intervals of width 0.1 nm from 108 to 124 nm:
172
173	wincr = 0.1
174	do iw = 1080, 1239, 1
175	kw = kw + 1
176	wl(kw) = real(iw)/10.
177	wu(kw) = wl(kw) + wincr
178	wc(kw) = (wl(kw) + wu(kw))/2.
179	end do
180
181	! define wavelength intervals of width 0.5 nm from 124 to 175 nm:
182
183	wincr = 0.5
184	do iw = 1240, 1745, 5
185	kw = kw + 1
186	wl(kw) = real(iw)/10.
187	wu(kw) = wl(kw) + wincr
188	wc(kw) = (wl(kw) + wu(kw))/2.
189	end do
190
191	! define wavelength intervals of width 0.01 nm from 175 to 205 nm:
192
193	wincr = 0.01
194	do iw = 17500, 20499, 1
195	kw = kw + 1
196	wl(kw) = real(iw)/100.
197	wu(kw) = wl(kw) + wincr
198	wc(kw) = (wl(kw) + wu(kw))/2.
199	end do
200
201	! define wavelength intervals of width 0.5 nm from 205 to 365 nm:
202
203	wincr = 0.5
204	do iw = 2050, 3645, 5
205	kw = kw + 1
206	wl(kw) = real(iw)/10.
207	wu(kw) = wl(kw) + wincr
208	wc(kw) = (wl(kw) + wu(kw))/2.
209	end do
210
211	! define wavelength intervals of width 5.0 nm from 365 to 855 nm:
212
213	wincr = 5.0
214	do iw = 365, 850, 5
215	kw = kw + 1
216	wl(kw) = real(iw)
217	wu(kw) = wl(kw) + wincr
218	wc(kw) = (wl(kw) + wu(kw))/2.
219	end do
220	wl(kw+1) = wu(kw)
221
222	!============================================================
223
224	else if (mopt == 2) then ! low-res
225
226	! define wavelength intervals of width 6.0 nm from 0 to 60 nm:
227
228	kw = 0
229	wincr = 6.0
230	DO iw = 0, 54, 6
231	kw = kw + 1
232	wl(kw) = real(iw)
233	wu(kw) = wl(kw) + wincr
234	wc(kw) = (wl(kw) + wu(kw))/2.
235	END DO
236
237	! define wavelength intervals of width 2.0 nm from 60 to 80 nm:
238
239	wincr = 2.0
240	DO iw = 60, 78, 2
241	kw = kw + 1
242	wl(kw) = real(iw)
243	wu(kw) = wl(kw) + wincr
244	wc(kw) = (wl(kw) + wu(kw))/2.
245	END DO
246
247	! define wavelength intervals of width 5.0 nm from 80 to 85 nm:
248
249	wincr = 5.0
250	DO iw = 80, 80
251	kw = kw + 1
252	wl(kw) = real(iw)
253	wu(kw) = wl(kw) + wincr
254	wc(kw) = (wl(kw) + wu(kw))/2.
255	END DO
256
257	! define wavelength intervals of width 2.0 nm from 85 to 117 nm:
258
259	wincr = 2.0
260	DO iw = 85, 115, 2
261	kw = kw + 1
262	wl(kw) = real(iw)
263	wu(kw) = wl(kw) + wincr
264	wc(kw) = (wl(kw) + wu(kw))/2.
265	END DO
266
267	! define wavelength intervals of width 3.0 nm from 117 to 120 nm:
268
269	wincr = 3.0
270	DO iw = 117, 117
271	kw = kw + 1
272	wl(kw) = real(iw)
273	wu(kw) = wl(kw) + wincr
274	wc(kw) = (wl(kw) + wu(kw))/2.
275	END DO
276
277	! define wavelength intervals of width 0.2 nm from 120 to 123 nm:
278
279	wincr = 0.2
280	DO iw = 1200, 1228, 2
281	kw = kw + 1
282	wl(kw) = real(iw)/10.
283	wu(kw) = wl(kw) + wincr
284	wc(kw) = (wl(kw) + wu(kw))/2.
285	ENDDO
286
287	! define wavelength intervals of width 5.0 nm from 123 to 163 nm:
288
289	wincr = 5.0
290	DO iw = 123, 158, 5
291	kw = kw + 1
292	wl(kw) = real(iw)
293	wu(kw) = wl(kw) + wincr
294	wc(kw) = (wl(kw) + wu(kw))/2.
295	ENDDO
296
297	! define wavelength intervals of width 2.0 nm from 163 to 175 nm:
298
299	wincr = 2.0
300	DO iw = 163, 173, 2
301	kw = kw + 1
302	wl(kw) = real(iw)
303	wu(kw) = wl(kw) + wincr
304	wc(kw) = (wl(kw) + wu(kw))/2.
305	ENDDO
306
307	! define wavelength intervals of width 0.5 nm from 175 to 205 nm:
308
309	wincr = 0.5
310	DO iw = 1750, 2045, 5
311	kw = kw + 1
312	wl(kw) = real(iw)/10.
313	wu(kw) = wl(kw) + wincr
314	wc(kw) = (wl(kw) + wu(kw))/2.
315	ENDDO
316
317	! define wavelength intervals of width 5.0 nm from 205 to 245 nm:
318
319	wincr = 5.
320	DO iw = 205, 240, 5
321	kw = kw + 1
322	wl(kw) = real(iw)
323	wu(kw) = wl(kw) + wincr
324	wc(kw) = (wl(kw) + wu(kw))/2.
325	ENDDO
326
327	! define wavelength intervals of width 10.0 nm from 245 to 415 nm:
328
329	wincr = 10.0
330	DO iw = 245, 405, 10
331	kw = kw + 1
332	wl(kw) = real(iw)
333	wu(kw) = wl(kw) + wincr
334	wc(kw) = (wl(kw) + wu(kw))/2.
335	ENDDO
336
337	! define wavelength intervals of width 50.0 nm from 415 to 815 nm:
338
339	wincr = 50.0
340	DO iw = 415, 815, 50
341	kw = kw + 1
342	wl(kw) = real(iw)
343	wu(kw) = wl(kw) + wincr
344	wc(kw) = (wl(kw) + wu(kw))/2.
345	ENDDO
346
347	wl(kw+1) = wu(kw)
348
349	end if ! mopt
350
351	print*, 'number of spectral intervals : ', kw+1
352
353	end subroutine gridw
354
355	!==============================================================================
356
357	subroutine rdsolarflux(nw,wl,wc,f)
358
359	! Read and re-grid solar flux data.
360
361	use datafile_mod, only: datadir
362
363	implicit none
364
365	! input
366
367	integer :: nw ! number of wavelength grid points
368	real, dimension(nw) :: wl, wc ! lower and central wavelength for each interval
369
370	! output
371
372	real, dimension(nw) :: f ! solar flux (w.m-2.nm-1)
373
374	! local
375
376	integer, parameter :: kdata = 20000 ! max dimension of input solar flux
377	integer :: msun ! choice of solar flux
378	integer :: iw, nhead, ihead, n, i, ierr, kin
379
380	real, parameter :: deltax = 1.e-4
381	real, dimension(kdata) :: x1, y1 ! input solar flux
382	real, dimension(nw) :: yg1 ! gridded solar flux
383
384	character(len=100) :: fil
385
386	kin = 10 ! input logical unit
387
388	! select desired extra-terrestrial solar irradiance, using msun:
389
390	! 18 = atlas3_thuillier_tuv.txt 0-900 nm November 1994
391	! Thuillier et al., Adv. Space. Res., 34, 256-261, 2004
392
393	msun = 18
394
395	if (msun == 18) THEN
396
397	fil = trim(datadir)//'/solar_fluxes/atlas3_thuillier_tuv.txt'
398	print*, 'solar flux : ', fil
399	OPEN(UNIT=kin,FILE=fil,STATUS='old')
400	nhead = 9
401	n = 19193
402	DO ihead = 1, nhead
403	READ(kin,*)
404	ENDDO
405	DO i = 1, n
406	READ(kin,*) x1(i), y1(i)
407	y1(i) = y1(i)*1.e-3 ! mw -> w
408	ENDDO
409	CLOSE (kin)
410	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
411	CALL addpnt(x1,y1,kdata,n, 0.,0.)
412	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
413	CALL addpnt(x1,y1,kdata,n, 1.e+38,0.)
414	CALL inter2(nw,wl,yg1,n,x1,y1,ierr)
415
416	IF (ierr .NE. 0) THEN
417	WRITE(,) ierr, fil
418	STOP
419	ENDIF
420
421	! convert to photon.s-1.nm-1.cm-2
422	! 5.039e11 = 1.e-41e-9/(hc = 6.62e-342.998e8)
423
424	DO iw = 1, nw-1
425	f(iw) = yg1(iw)wc(iw)5.039e11
426	! write(25,*) iw, wc(iw), f(iw)
427	ENDDO
428
429	end if
430
431	end subroutine rdsolarflux
432
433	!==============================================================================
434
435	subroutine addpnt ( x, y, ld, n, xnew, ynew )
436
437	!-----------------------------------------------------------------------------*
438	!= PURPOSE: =*
439	!= Add a point <xnew,ynew> to a set of data pairs <x,y>. x must be in =*
440	!= ascending order =*
441	!-----------------------------------------------------------------------------*
442	!= PARAMETERS: =*
443	!= X - REAL vector of length LD, x-coordinates (IO)=*
444	!= Y - REAL vector of length LD, y-values (IO)=*
445	!= LD - INTEGER, dimension of X, Y exactly as declared in the calling (I)=*
446	!= program =*
447	!= N - INTEGER, number of elements in X, Y. On entry, it must be: (IO)=*
448	!= N < LD. On exit, N is incremented by 1. =*
449	!= XNEW - REAL, x-coordinate at which point is to be added (I)=*
450	!= YNEW - REAL, y-value of point to be added (I)=*
451	!-----------------------------------------------------------------------------*
452
453	IMPLICIT NONE
454
455	! calling parameters
456
457	INTEGER ld, n
458	REAL x(ld), y(ld)
459	REAL xnew, ynew
460	INTEGER ierr
461
462	! local variables
463
464	INTEGER insert
465	INTEGER i
466
467	!-----------------------------------------------------------------------
468
469	! initialize error flag
470
471	ierr = 0
472
473	! check n<ld to make sure x will hold another point
474
475	IF (n .GE. ld) THEN
476	WRITE(0,*) '>>> ERROR (ADDPNT) <<< Cannot expand array '
477	WRITE(0,*) ' All elements used.'
478	STOP
479	ENDIF
480
481	insert = 1
482	i = 2
483
484	! check, whether x is already sorted.
485	! also, use this loop to find the point at which xnew needs to be inserted
486	! into vector x, if x is sorted.
487
488	10 CONTINUE
489	IF (i .LT. n) THEN
490	IF (x(i) .LT. x(i-1)) THEN
491	print*, x(i-1), x(i)
492	WRITE(0,*) '>>> ERROR (ADDPNT) <<< x-data must be in ascending order!'
493	STOP
494	ELSE
495	IF (xnew .GT. x(i)) insert = i + 1
496	ENDIF
497	i = i+1
498	GOTO 10
499	ENDIF
500
501	! if <xnew,ynew> needs to be appended at the end, just do so,
502	! otherwise, insert <xnew,ynew> at position INSERT
503
504	IF ( xnew .GT. x(n) ) THEN
505
506	x(n+1) = xnew
507	y(n+1) = ynew
508
509	ELSE
510
511	! shift all existing points one index up
512
513	DO i = n, insert, -1
514	x(i+1) = x(i)
515	y(i+1) = y(i)
516	ENDDO
517
518	! insert new point
519
520	x(insert) = xnew
521	y(insert) = ynew
522
523	ENDIF
524
525	! increase total number of elements in x, y
526
527	n = n+1
528
529	end subroutine addpnt
530
531	!==============================================================================
532
533	subroutine inter2(ng,xg,yg,n,x,y,ierr)
534
535	!-----------------------------------------------------------------------------*
536	!= PURPOSE: =*
537	!= Map input data given on single, discrete points onto a set of target =*
538	!= bins. =*
539	!= The original input data are given on single, discrete points of an =*
540	!= arbitrary grid and are being linearly interpolated onto a specified set =*
541	!= of target bins. In general, this is the case for most of the weighting =*
542	!= functions (action spectra, molecular cross section, and quantum yield =*
543	!= data), which have to be matched onto the specified wavelength intervals. =*
544	!= The average value in each target bin is found by averaging the trapezoi- =*
545	!= dal area underneath the input data curve (constructed by linearly connec-=*
546	!= ting the discrete input values). =*
547	!= Some caution should be used near the endpoints of the grids. If the =*
548	!= input data set does not span the range of the target grid, an error =*
549	!= message is printed and the execution is stopped, as extrapolation of the =*
550	!= data is not permitted. =*
551	!= If the input data does not encompass the target grid, use ADDPNT to =*
552	!= expand the input array. =*
553	!-----------------------------------------------------------------------------*
554	!= PARAMETERS: =*
555	!= NG - INTEGER, number of bins + 1 in the target grid (I)=*
556	!= XG - REAL, target grid (e.g., wavelength grid); bin i is defined (I)=*
557	!= as [XG(i),XG(i+1)] (i = 1..NG-1) =*
558	!= YG - REAL, y-data re-gridded onto XG, YG(i) specifies the value for (O)=*
559	!= bin i (i = 1..NG-1) =*
560	!= N - INTEGER, number of points in input grid (I)=*
561	!= X - REAL, grid on which input data are defined (I)=*
562	!= Y - REAL, input y-data (I)=*
563	!-----------------------------------------------------------------------------*
564
565	IMPLICIT NONE
566
567	! input:
568	INTEGER ng, n
569	REAL x(n), y(n), xg(ng)
570
571	! output:
572	REAL yg(ng)
573
574	! local:
575	REAL area, xgl, xgu
576	REAL darea, slope
577	REAL a1, a2, b1, b2
578	INTEGER ngintv
579	INTEGER i, k, jstart
580	INTEGER ierr
581	!_______________________________________________________________________
582
583	ierr = 0
584
585	! test for correct ordering of data, by increasing value of x
586
587	DO 10, i = 2, n
588	IF (x(i) .LE. x(i-1)) THEN
589	ierr = 1
590	WRITE(,)'data not sorted'
591	WRITE(,) x(i), x(i-1)
592	RETURN
593	ENDIF
594	10 CONTINUE
595
596	DO i = 2, ng
597	IF (xg(i) .LE. xg(i-1)) THEN
598	ierr = 2
599	WRITE(0,*) '>>> ERROR (inter2) <<< xg-grid not sorted!'
600	RETURN
601	ENDIF
602	ENDDO
603
604	! check for xg-values outside the x-range
605
606	IF ( (x(1) .GT. xg(1)) .OR. (x(n) .LT. xg(ng)) ) THEN
607	WRITE(0,*) '>>> ERROR (inter2) <<< Data do not span grid. '
608	WRITE(0,*) ' Use ADDPNT to expand data and re-run.'
609	STOP
610	ENDIF
611
612	! find the integral of each grid interval and use this to
613	! calculate the average y value for the interval
614	! xgl and xgu are the lower and upper limits of the grid interval
615
616	jstart = 1
617	ngintv = ng - 1
618	DO 50, i = 1,ngintv
619
620	! initialize:
621
622	area = 0.0
623	xgl = xg(i)
624	xgu = xg(i+1)
625
626	! discard data before the first grid interval and after the
627	! last grid interval
628	! for internal grid intervals, start calculating area by interpolating
629	! between the last point which lies in the previous interval and the
630	! first point inside the current interval
631
632	k = jstart
633	IF (k .LE. n-1) THEN
634
635	! if both points are before the first grid, go to the next point
636	30 CONTINUE
637	IF (x(k+1) .LE. xgl) THEN
638	jstart = k - 1
639	k = k+1
640	IF (k .LE. n-1) GO TO 30
641	ENDIF
642
643
644	! if the last point is beyond the end of the grid, complete and go to the next
645	! grid
646	40 CONTINUE
647	IF ((k .LE. n-1) .AND. (x(k) .LT. xgu)) THEN
648
649	jstart = k-1
650
651	! compute x-coordinates of increment
652
653	a1 = MAX(x(k),xgl)
654	a2 = MIN(x(k+1),xgu)
655
656	! if points coincide, contribution is zero
657
658	IF (x(k+1).EQ.x(k)) THEN
659	darea = 0.e0
660	ELSE
661	slope = (y(k+1) - y(k))/(x(k+1) - x(k))
662	b1 = y(k) + slope*(a1 - x(k))
663	b2 = y(k) + slope*(a2 - x(k))
664	darea = (a2 - a1)*(b2 + b1)/2.
665	ENDIF
666
667	! find the area under the trapezoid from a1 to a2
668
669	area = area + darea
670
671	! go to next point
672
673	k = k+1
674	GO TO 40
675
676	ENDIF
677	ENDIF
678
679	! calculate the average y after summing the areas in the interval
680
681	yg(i) = area/(xgu - xgl)
682
683	50 CONTINUE
684
685	end subroutine inter2
686
687	!==============================================================================
688
689	subroutine rdxsco2(nw,wl,xsco2_195,xsco2_295,xsco2_370,yieldco2)
690
691	!-----------------------------------------------------------------------------*
692	!= PURPOSE: =*
693	!= Read and grid CO2 absorption cross-sections and photodissociation yield =*
694	!-----------------------------------------------------------------------------*
695	!= PARAMETERS: =*
696	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
697	!= wavelength grid =*
698	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
699	!= working wavelength grid =*
700	!= XSCO2 - REAL, molecular absoprtion cross section (cm^2) of CO2 at (O)=*
701	!= each specified wavelength =*
702	!-----------------------------------------------------------------------------*
703
704	use datafile_mod, only: datadir
705
706	implicit none
707
708	! input
709
710	integer :: nw ! number of wavelength grid points
711	real, dimension(nw) :: wl ! lower and central wavelength for each interval
712
713	! output
714
715	real, dimension(nw) :: xsco2_195, xsco2_295, xsco2_370 ! co2 cross-sections (cm2)
716	real, dimension(nw) :: yieldco2 ! co2 photodissociation yield
717
718	! local
719
720	integer, parameter :: kdata = 42000
721	real, parameter :: deltax = 1.e-4
722	real, dimension(kdata) :: x1, y1, y2, y3, xion, ion
723	real, dimension(nw) :: yg
724	real :: xl, xu
725	integer :: ierr, i, l, n, n1, n2, n3, n4
726	CHARACTER*100 fil
727
728	integer :: kin, kout ! input/ouput logical units
729
730	kin = 10
731	kout = 30
732
733	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
734	!
735	! CO2 absorption cross-sections
736	!
737	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
738	!
739	! 195K: huestis and berkowitz (2010) + Starck et al. (2006)
740	! + Yoshino et al. (1996) + Parkinson et al. (2003) + extrapolation
741	!
742	! 295K: huestis and berkowitz (2010) + Starck et al. (2006)
743	! + Yoshino et al. (1996) + Parkinson et al. (2003) + extrapolation
744	!
745	! 370K: huestis and berkowitz (2010) + Starck et al. (2006)
746	! + Lewis and Carver (1983) + extrapolation
747	!
748	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
749
750	n1 = 40769
751	n2 = 41586
752	n3 = 10110
753
754	! 195K:
755
756	fil = trim(datadir)//'/cross_sections/co2_euv_uv_2018_195k.txt'
757	print*, 'section efficace CO2 195K: ', fil
758
759	OPEN(UNIT=kin,FILE=fil,STATUS='old')
760	DO i = 1,11
761	read(kin,*)
762	END DO
763
764	DO i = 1, n1
765	READ(kin,*) x1(i), y1(i)
766	END DO
767	CLOSE (kin)
768
769	CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
770	CALL addpnt(x1,y1,kdata,n1, 0.,0.)
771	CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
772	CALL addpnt(x1,y1,kdata,n1, 1.e+38,0.)
773	CALL inter2(nw,wl,yg,n1,x1,y1,ierr)
774	IF (ierr .NE. 0) THEN
775	WRITE(,) ierr, fil
776	STOP
777	ENDIF
778
779	DO l = 1, nw-1
780	xsco2_195(l) = yg(l)
781	END DO
782
783	! 295K:
784
785	fil = trim(datadir)//'/cross_sections/co2_euv_uv_2018_295k.txt'
786	print*, 'section efficace CO2 295K: ', fil
787
788	OPEN(UNIT=kin,FILE=fil,STATUS='old')
789	DO i = 1,11
790	read(kin,*)
791	END DO
792
793	DO i = 1, n2
794	READ(kin,*) x1(i), y1(i)
795	END DO
796	CLOSE (kin)
797
798	CALL addpnt(x1,y1,kdata,n2,x1(1)*(1.-deltax),0.)
799	CALL addpnt(x1,y1,kdata,n2, 0.,0.)
800	CALL addpnt(x1,y1,kdata,n2,x1(n2)*(1.+deltax),0.)
801	CALL addpnt(x1,y1,kdata,n2, 1.e+38,0.)
802	CALL inter2(nw,wl,yg,n2,x1,y1,ierr)
803	IF (ierr .NE. 0) THEN
804	WRITE(,) ierr, fil
805	STOP
806	ENDIF
807
808	DO l = 1, nw-1
809	xsco2_295(l) = yg(l)
810	END DO
811
812	! 370K:
813
814	fil = trim(datadir)//'/cross_sections/co2_euv_uv_2018_370k.txt'
815	print*, 'section efficace CO2 370K: ', fil
816
817	OPEN(UNIT=kin,FILE=fil,STATUS='old')
818	DO i = 1,11
819	read(kin,*)
820	END DO
821
822	DO i = 1, n3
823	READ(kin,*) x1(i), y1(i)
824	END DO
825	CLOSE (kin)
826
827	CALL addpnt(x1,y1,kdata,n3,x1(1)*(1.-deltax),0.)
828	CALL addpnt(x1,y1,kdata,n3, 0.,0.)
829	CALL addpnt(x1,y1,kdata,n3,x1(n3)*(1.+deltax),0.)
830	CALL addpnt(x1,y1,kdata,n3, 1.e+38,0.)
831	CALL inter2(nw,wl,yg,n3,x1,y1,ierr)
832	IF (ierr .NE. 0) THEN
833	WRITE(,) ierr, fil
834	STOP
835	ENDIF
836
837	DO l = 1, nw-1
838	xsco2_370(l) = yg(l)
839	END DO
840
841	! photodissociation yield:
842
843	fil = trim(datadir)//'/cross_sections/efdis_co2-o2_schunkandnagy2000.txt'
844	print*, 'photodissociation yield CO2: ', fil
845
846	OPEN(UNIT=kin,FILE=fil,STATUS='old')
847
848	do i = 1,3
849	read(kin,*)
850	end do
851
852	n4 = 17
853	do i = 1, n4
854	read(kin,*) xl, xu, ion(i)
855	xion(i) = (xl + xu)/2.
856	ion(i) = max(ion(i), 0.)
857	end do
858	close(kin)
859
860	CALL addpnt(xion,ion,kdata,n4,xion(1)*(1.-deltax),0.)
861	CALL addpnt(xion,ion,kdata,n4, 0.,0.)
862	CALL addpnt(xion,ion,kdata,n4,xion(n4)*(1.+deltax),1.)
863	CALL addpnt(xion,ion,kdata,n4, 1.e+38,1.)
864	CALL inter2(nw,wl,yieldco2,n4,xion,ion,ierr)
865	IF (ierr .NE. 0) THEN
866	WRITE(,) ierr, fil
867	STOP
868	ENDIF
869
870	! DO l = 1, nw-1
871	! write(kout,*) wl(l), xsco2_195(l),
872	! $ xsco2_295(l),
873	! $ xsco2_370(l),
874	! $ yieldco2(l)
875	! END DO
876
877	end subroutine rdxsco2
878
879	!==============================================================================
880
881	subroutine rdxso2(nw,wl,xso2_150,xso2_200,xso2_250,xso2_300,yieldo2)
882
883	!-----------------------------------------------------------------------------*
884	!= PURPOSE: =*
885	!= Read and grid O2 cross-sections and photodissociation yield =*
886	!-----------------------------------------------------------------------------*
887	!= PARAMETERS: =*
888	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
889	!= wavelength grid =*
890	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
891	!= working wavelength grid =*
892	!= XSO2 - REAL, molecular absorption cross section =*
893	!-----------------------------------------------------------------------------*
894
895	use datafile_mod, only: datadir
896
897	implicit none
898
899	! input
900
901	integer :: nw ! number of wavelength grid points
902	real, dimension(nw) :: wl ! lower and central wavelength for each interval
903
904	! output
905
906	real, dimension(nw) :: xso2_150, xso2_200, xso2_250, xso2_300 ! o2 cross-sections (cm2)
907	real, dimension(nw) :: yieldo2 ! o2 photodissociation yield
908
909	! local
910
911	integer, parameter :: kdata = 18000
912	real, parameter :: deltax = 1.e-4
913	real, dimension(kdata) :: x1, y1, x2, y2, x3, y3, x4, y4
914	real, dimension(kdata) :: xion, ion
915	real :: factor, xl, xu, dummy
916	integer :: i, ierr, n, n1, n2, n3, n4, nhead
917	integer :: kin, kout ! input/output logical units
918	character*100 fil
919
920	kin = 10
921	kout = 30
922
923	! read o2 cross section data
924
925	nhead = 22
926	n = 17434
927
928	fil = trim(datadir)//'/cross_sections/o2_composite_2018_150K.txt'
929	print*, 'section efficace O2 150K: ', fil
930	OPEN(UNIT=kin,FILE=fil,STATUS='old')
931
932	DO i = 1,nhead
933	read(kin,*)
934	END DO
935
936	n1 = n
937	DO i = 1, n1
938	READ(kin,*) x1(i), y1(i)
939	END DO
940	CLOSE (kin)
941
942	CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
943	CALL addpnt(x1,y1,kdata,n1, 0.,0.)
944	CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
945	CALL addpnt(x1,y1,kdata,n1, 1.e+38,0.)
946	CALL inter2(nw,wl,xso2_150,n1,x1,y1,ierr)
947	IF (ierr .NE. 0) THEN
948	WRITE(,) ierr, fil
949	STOP
950	ENDIF
951
952	fil = trim(datadir)//'/cross_sections/o2_composite_2018_200K.txt'
953	print*, 'section efficace O2 200K: ', fil
954	OPEN(UNIT=kin,FILE=fil,STATUS='old')
955
956	DO i = 1,nhead
957	read(kin,*)
958	END DO
959
960	n2 = n
961	DO i = 1, n2
962	READ(kin,*) x2(i), y2(i)
963	END DO
964	CLOSE (kin)
965
966	CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
967	CALL addpnt(x2,y2,kdata,n2, 0.,0.)
968	CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
969	CALL addpnt(x2,y2,kdata,n2, 1.e+38,0.)
970	CALL inter2(nw,wl,xso2_200,n2,x2,y2,ierr)
971	IF (ierr .NE. 0) THEN
972	WRITE(,) ierr, fil
973	STOP
974	ENDIF
975
976	fil = trim(datadir)//'/cross_sections/o2_composite_2018_250K.txt'
977	print*, 'section efficace O2 250K: ', fil
978	OPEN(UNIT=kin,FILE=fil,STATUS='old')
979
980	DO i = 1,nhead
981	read(kin,*)
982	END DO
983
984	n3 = n
985	DO i = 1, n3
986	READ(kin,*) x3(i), y3(i)
987	END DO
988	CLOSE (kin)
989
990	CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax),0.)
991	CALL addpnt(x3,y3,kdata,n3, 0.,0.)
992	CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.)
993	CALL addpnt(x3,y3,kdata,n3, 1.e+38,0.)
994	CALL inter2(nw,wl,xso2_250,n3,x3,y3,ierr)
995	IF (ierr .NE. 0) THEN
996	WRITE(,) ierr, fil
997	STOP
998	ENDIF
999
1000	fil = trim(datadir)//'/cross_sections/o2_composite_2018_300K.txt'
1001	print*, 'section efficace O2 300K: ', fil
1002	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1003
1004	DO i = 1,nhead
1005	read(kin,*)
1006	END DO
1007
1008	n4 = n
1009	DO i = 1, n4
1010	READ(kin,*) x4(i), y4(i)
1011	END DO
1012	CLOSE (kin)
1013
1014	CALL addpnt(x4,y4,kdata,n4,x4(1)*(1.-deltax),0.)
1015	CALL addpnt(x4,y4,kdata,n4, 0.,0.)
1016	CALL addpnt(x4,y4,kdata,n4,x4(n4)*(1.+deltax),0.)
1017	CALL addpnt(x4,y4,kdata,n4, 1.e+38,0.)
1018	CALL inter2(nw,wl,xso2_300,n4,x4,y4,ierr)
1019	IF (ierr .NE. 0) THEN
1020	WRITE(,) ierr, fil
1021	STOP
1022	ENDIF
1023
1024	! photodissociation yield
1025
1026	fil = trim(datadir)//'/cross_sections/efdis_co2-o2_schunkandnagy2000.txt'
1027	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1028
1029	do i = 1,11
1030	read(kin,*)
1031	end do
1032
1033	n = 9
1034	DO i = 1, n
1035	READ(kin,*) xl, xu, dummy, ion(i)
1036	xion(i) = (xl + xu)/2.
1037	ion(i) = max(ion(i), 0.)
1038	END DO
1039	CLOSE (kin)
1040
1041	CALL addpnt(xion,ion,kdata,n,xion(1)*(1.-deltax),0.)
1042	CALL addpnt(xion,ion,kdata,n, 0.,0.)
1043	CALL addpnt(xion,ion,kdata,n,xion(n)*(1.+deltax),1.)
1044	CALL addpnt(xion,ion,kdata,n, 1.e+38,1.)
1045	CALL inter2(nw,wl,yieldo2,n,xion,ion,ierr)
1046	IF (ierr .NE. 0) THEN
1047	WRITE(,) ierr, fil
1048	STOP
1049	ENDIF
1050
1051	end subroutine rdxso2
1052
1053	!==============================================================================
1054
1055	subroutine rdxso3(nw,wl,xso3_218,xso3_298)
1056
1057	!-----------------------------------------------------------------------------*
1058	!= PURPOSE: =*
1059	!= Read ozone molecular absorption cross section. Re-grid data to match =*
1060	!= specified wavelength working grid. =*
1061	!-----------------------------------------------------------------------------*
1062	!= PARAMETERS: =*
1063	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
1064	!= wavelength grid =*
1065	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
1066	!= working wavelength grid =*
1067	!= XSO3_218 REAL, molecular absoprtion cross section (cm^2) of O3 at (O)=*
1068	!= each specified wavelength (JPL 2006) 218 K =*
1069	!= XSO3_298 REAL, molecular absoprtion cross section (cm^2) of O3 at (O)=*
1070	!= each specified wavelength (JPL 2006) 298 K =*
1071	!-----------------------------------------------------------------------------*
1072
1073	use datafile_mod, only: datadir
1074
1075	implicit none
1076
1077	! input
1078
1079	integer :: nw ! number of wavelength grid points
1080	real, dimension(nw) :: wl ! lower and central wavelength for each interval
1081
1082	! output
1083
1084	real, dimension(nw) :: xso3_218, xso3_298 ! o3 cross-sections (cm2)
1085
1086	! local
1087
1088	integer, parameter :: kdata = 200
1089	real, parameter :: deltax = 1.e-4
1090	real, dimension(kdata) :: x1, x2, y1, y2
1091	real, dimension(nw) :: yg
1092	real :: a1, a2
1093
1094	integer :: i, ierr, iw, n, n1, n2
1095	integer :: kin, kout ! input/output logical units
1096
1097	character*100 fil
1098
1099	! JPL 2006 218 K
1100
1101	fil = trim(datadir)//'/cross_sections/o3_cross-sections_jpl_2006_218K.txt'
1102	print*, 'section efficace O3 218K: ', fil
1103
1104	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1105	n1 = 167
1106	DO i = 1, n1
1107	READ(kin,*) a1, a2, y1(i)
1108	x1(i) = (a1+a2)/2.
1109	END DO
1110	CLOSE (kin)
1111
1112	CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
1113	CALL addpnt(x1,y1,kdata,n1, 0.,0.)
1114	CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
1115	CALL addpnt(x1,y1,kdata,n1, 1.e+38,0.)
1116	CALL inter2(nw,wl,yg,n1,x1,y1,ierr)
1117	IF (ierr .NE. 0) THEN
1118	WRITE(,) ierr, fil
1119	STOP
1120	ENDIF
1121
1122	DO iw = 1, nw-1
1123	xso3_218(iw) = yg(iw)
1124	END DO
1125
1126	! JPL 2006 298 K
1127
1128	fil = trim(datadir)//'/cross_sections/o3_cross-sections_jpl_2006_298K.txt'
1129	print*, 'section efficace O3 298K: ', fil
1130
1131	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1132	n2 = 167
1133	DO i = 1, n2
1134	READ(kin,*) a1, a2, y2(i)
1135	x2(i) = (a1+a2)/2.
1136	END DO
1137	CLOSE (kin)
1138
1139	CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
1140	CALL addpnt(x2,y2,kdata,n2, 0.,0.)
1141	CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
1142	CALL addpnt(x2,y2,kdata,n2, 1.e+38,0.)
1143	CALL inter2(nw,wl,yg,n2,x2,y2,ierr)
1144	IF (ierr .NE. 0) THEN
1145	WRITE(,) ierr, fil
1146	STOP
1147	ENDIF
1148
1149	DO iw = 1, nw-1
1150	xso3_298(iw) = yg(iw)
1151	END DO
1152
1153	end subroutine rdxso3
1154
1155	!==============================================================================
1156
1157	subroutine rdxsh2o(nw, wl, yg)
1158
1159	!-----------------------------------------------------------------------------*
1160	!= PURPOSE: =*
1161	!= Read H2O molecular absorption cross section. Re-grid data to match =*
1162	!= specified wavelength working grid. =*
1163	!-----------------------------------------------------------------------------*
1164	!= PARAMETERS: =*
1165	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
1166	!= wavelength grid =*
1167	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
1168	!= working wavelength grid =*
1169	!= YG - REAL, molecular absoprtion cross section (cm^2) of H2O at (O)=*
1170	!= each specified wavelength =*
1171	!-----------------------------------------------------------------------------*
1172
1173	use datafile_mod, only: datadir
1174
1175	IMPLICIT NONE
1176
1177	! input
1178
1179	integer :: nw ! number of wavelength grid points
1180	real, dimension(nw) :: wl ! lower and central wavelength for each interval
1181
1182	! output
1183
1184	real, dimension(nw) :: yg ! h2o cross-sections (cm2)
1185
1186	! local
1187
1188	integer, parameter :: kdata = 500
1189	real, parameter :: deltax = 1.e-4
1190	REAL x1(kdata)
1191	REAL y1(kdata)
1192	INTEGER ierr
1193	INTEGER i, n
1194	CHARACTER*100 fil
1195	integer :: kin, kout ! input/output logical units
1196
1197	kin = 10
1198
1199	fil = trim(datadir)//'/cross_sections/h2o_composite_250K.txt'
1200	print*, 'section efficace H2O: ', fil
1201
1202	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1203
1204	DO i = 1,26
1205	read(kin,*)
1206	END DO
1207
1208	n = 420
1209	DO i = 1, n
1210	READ(kin,*) x1(i), y1(i)
1211	END DO
1212	CLOSE (kin)
1213
1214	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1215	CALL addpnt(x1,y1,kdata,n, 0.,0.)
1216	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1217	CALL addpnt(x1,y1,kdata,n, 1.e+38,0.)
1218	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
1219	IF (ierr .NE. 0) THEN
1220	WRITE(,) ierr, fil
1221	STOP
1222	ENDIF
1223
1224	end subroutine rdxsh2o
1225
1226	!==============================================================================
1227
1228	subroutine rdxsh2o2(nw, wl, xsh2o2)
1229
1230	!-----------------------------------------------------------------------------*
1231	!= PURPOSE: =*
1232	!= Read and grid H2O2 cross-sections
1233	!= H2O2 + hv -> 2 OH =*
1234	!= Cross section: Schuergers and Welge, Z. Naturforsch. 23a (1968) 1508 =*
1235	!= from 125 to 185 nm, then JPL97 from 190 to 350 nm. =*
1236	!-----------------------------------------------------------------------------*
1237	!= PARAMETERS: =*
1238	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
1239	!= wavelength grid =*
1240	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
1241	!= working wavelength grid =*
1242	!-----------------------------------------------------------------------------*
1243
1244	use datafile_mod, only: datadir
1245
1246	implicit none
1247
1248	! input
1249
1250	integer :: nw ! number of wavelength grid points
1251	real, dimension(nw) :: wl ! lower wavelength for each interval
1252
1253	! output
1254
1255	real, dimension(nw) :: xsh2o2 ! h2o2 cross-sections (cm2)
1256
1257	! local
1258
1259	real, parameter :: deltax = 1.e-4
1260	integer, parameter :: kdata = 100
1261	real, dimension(kdata) :: x1, y1
1262	real, dimension(nw) :: yg
1263	integer :: i, ierr, iw, n, idum
1264	integer :: kin, kout ! input/output logical units
1265	character*100 fil
1266
1267	kin = 10
1268
1269	! read cross-sections
1270
1271	fil = trim(datadir)//'/cross_sections/h2o2_composite.txt'
1272	print*, 'section efficace H2O2: ', fil
1273
1274	OPEN(kin,FILE=fil,STATUS='OLD')
1275	READ(kin,*) idum,n
1276	DO i = 1, idum-2
1277	READ(kin,*)
1278	ENDDO
1279	DO i = 1, n
1280	READ(kin,*) x1(i), y1(i)
1281	ENDDO
1282	CLOSE (kin)
1283
1284	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1285	CALL addpnt(x1,y1,kdata,n, 0.,0.)
1286	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1287	CALL addpnt(x1,y1,kdata,n, 1.e+38,0.)
1288	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
1289	IF (ierr .NE. 0) THEN
1290	WRITE(,) ierr, fil
1291	STOP
1292	ENDIF
1293
1294	DO iw = 1, nw - 1
1295	xsh2o2(iw) = yg(iw)
1296	END DO
1297
1298	end subroutine rdxsh2o2
1299
1300	!==============================================================================
1301
1302	subroutine rdxsho2(nw, wl, yg)
1303
1304	!-----------------------------------------------------------------------------*
1305	!= PURPOSE: =*
1306	!= Read ho2 cross-sections =*
1307	!= JPL 2006 recommendation =*
1308	!-----------------------------------------------------------------------------*
1309	!= PARAMETERS: =*
1310	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
1311	!= wavelength grid =*
1312	!-----------------------------------------------------------------------------*
1313
1314	use datafile_mod, only: datadir
1315
1316	IMPLICIT NONE
1317
1318	! input
1319
1320	integer :: nw ! number of wavelength grid points
1321	real, dimension(nw) :: wl ! lower wavelength for each interval
1322
1323	! output
1324
1325	real, dimension(nw) :: yg ! ho2 cross-sections (cm2)
1326
1327	! local
1328
1329	real, parameter :: deltax = 1.e-4
1330	integer, parameter :: kdata = 100
1331	real, dimension(kdata) :: x1, y1
1332	integer :: i, n, ierr
1333	character*100 fil
1334	integer :: kin, kout ! input/output logical units
1335
1336	kin = 10
1337
1338	!*** cross sections from Sander et al. [2003]
1339
1340	fil = trim(datadir)//'/cross_sections/ho2_jpl2003.txt'
1341	print*, 'section efficace HO2: ', fil
1342
1343	OPEN(kin,FILE=fil,STATUS='OLD')
1344	READ(kin,*) n
1345	DO i = 1, n
1346	READ(kin,*) x1(i), y1(i)
1347	ENDDO
1348	CLOSE(kin)
1349
1350	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1351	CALL addpnt(x1,y1,kdata,n, 0.,0.)
1352	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1353	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
1354
1355	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
1356
1357	IF (ierr .NE. 0) THEN
1358	WRITE(,) ierr, fil
1359	STOP
1360	ENDIF
1361
1362	end subroutine rdxsho2
1363
1364	!==============================================================================
1365
1366	subroutine rdxsh2(nw, wl, wc, yg, yieldh2)
1367
1368	!-----------------------------------------------------------------------------*
1369	!= PURPOSE: =*
1370	!= Read h2 cross-sections and photodissociation yield =*
1371	!-----------------------------------------------------------------------------*
1372	!= PARAMETERS: =*
1373	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
1374	!= wavelength grid =*
1375	!-----------------------------------------------------------------------------*
1376
1377	use datafile_mod, only: datadir
1378
1379	implicit none
1380
1381	! input
1382
1383	integer :: nw ! number of wavelength grid points
1384	real, dimension(nw) :: wl, wc ! lower and central wavelength for each interval
1385
1386	! output
1387
1388	real, dimension(nw) :: yg ! h2 cross-sections (cm2)
1389	real, dimension(nw) :: yieldh2 ! photodissociation yield
1390
1391	! local
1392
1393	integer, parameter :: kdata = 1000
1394	real, parameter :: deltax = 1.e-4
1395	real, dimension(kdata) :: x1, y1, x2, y2
1396	real :: xl, xu
1397	integer :: i, iw, n, ierr
1398	integer :: kin, kout ! input/output logical units
1399	character*100 fil
1400
1401	kin = 10
1402
1403	! h2 cross sections
1404
1405	fil = trim(datadir)//'/cross_sections/h2secef.txt'
1406	print*, 'section efficace H2: ', fil
1407
1408	OPEN(kin,FILE=fil,STATUS='OLD')
1409
1410	n = 792
1411	read(kin,*) ! avoid first line with wavelength = 0.
1412	DO i = 1, n
1413	READ(kin,*) x1(i), y1(i)
1414	ENDDO
1415	CLOSE(kin)
1416
1417	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1418	CALL addpnt(x1,y1,kdata,n, 0.,0.)
1419	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1420	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
1421	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
1422
1423	IF (ierr .NE. 0) THEN
1424	WRITE(,) ierr, fil
1425	STOP
1426	ENDIF
1427
1428	! photodissociation yield
1429
1430	fil = trim(datadir)//'/cross_sections/h2_ionef_schunknagy2000.txt'
1431	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1432
1433	n = 19
1434	read(kin,*)
1435	DO i = 1, n
1436	READ(kin,*) xl, xu, y2(i)
1437	x2(i) = (xl + xu)/2.
1438	y2(i) = max(1. - y2(i),0.)
1439	END DO
1440	CLOSE (kin)
1441
1442	CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
1443	CALL addpnt(x2,y2,kdata,n, 0.,0.)
1444	CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
1445	CALL addpnt(x2,y2,kdata,n, 1.e+38,1.)
1446	CALL inter2(nw,wl,yieldh2,n,x2,y2,ierr)
1447	IF (ierr .NE. 0) THEN
1448	WRITE(,) ierr, fil
1449	STOP
1450	ENDIF
1451
1452	end subroutine rdxsh2
1453
1454	!==============================================================================
1455
1456	subroutine rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294,yldno2_248, yldno2_298)
1457
1458	!-----------------------------------------------------------------------------*
1459	!= PURPOSE: =*
1460	!= read and grid cross section + quantum yield for NO2 =*
1461	!= photolysis =*
1462	!= Jenouvrier et al., 1996 200-238 nm
1463	!= Vandaele et al., 1998 238-666 nm 220K and 294K
1464	!= quantum yield from jpl 2006
1465	!-----------------------------------------------------------------------------*
1466	!= PARAMETERS: =*
1467	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
1468	!= wavelength grid =*
1469	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
1470	!= working wavelength grid =*
1471	!= SQ - REAL, cross section x quantum yield (cm^2) for each (O)=*
1472	!= photolysis reaction defined, at each defined wavelength and =*
1473	!= at each defined altitude level =*
1474	!-----------------------------------------------------------------------------*
1475
1476	use datafile_mod, only: datadir
1477
1478	implicit none
1479
1480	! input
1481
1482	integer :: nw ! number of wavelength grid points
1483	real, dimension(nw) :: wl ! lower and central wavelength for each interval
1484
1485	! output
1486
1487	real, dimension(nw) :: xsno2, xsno2_220, xsno2_294 ! no2 cross-sections (cm2)
1488	real, dimension(nw) :: yldno2_248, yldno2_298 ! quantum yields at 248-298 k
1489
1490	! local
1491
1492	integer, parameter :: kdata = 28000
1493	real, parameter :: deltax = 1.e-4
1494	real, dimension(kdata) :: x1, x2, x3, x4, x5, y1, y2, y3, y4, y5
1495	real, dimension(nw) :: yg1, yg2, yg3, yg4, yg5
1496	real :: dum, qy
1497	integer :: i, iw, n, n1, n2, n3, n4, n5, ierr
1498	character*100 fil
1499	integer :: kin, kout ! input/output logical units
1500
1501	kin = 10
1502
1503	!*************** NO2 photodissociation
1504
1505	! Jenouvrier 1996 + Vandaele 1998 (JPL 2006)
1506
1507	fil = trim(datadir)//'/cross_sections/no2_xs_jenouvrier.txt'
1508	print*, 'section efficace NO2: ', fil
1509
1510	OPEN(UNIT=kin,FILE=fil,status='old')
1511	DO i = 1, 3
1512	READ(kin,*)
1513	ENDDO
1514	n1 = 10001
1515	DO i = 1, n1
1516	READ(kin,*) x1(i), y1(i)
1517	end do
1518
1519	CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax), 0.)
1520	CALL addpnt(x1,y1,kdata,n1, 0., 0.)
1521	CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
1522	CALL addpnt(x1,y1,kdata,n1, 1.e+38, 0.)
1523	CALL inter2(nw,wl,yg1,n1,x1,y1,ierr)
1524
1525	fil = trim(datadir)//'/cross_sections/no2_xs_vandaele_294K.txt'
1526	print*, 'section efficace NO2: ', fil
1527
1528	OPEN(UNIT=kin,FILE=fil,status='old')
1529	DO i = 1, 3
1530	READ(kin,*)
1531	ENDDO
1532	n2 = 27993
1533	DO i = 1, n2
1534	READ(kin,*) x2(i), y2(i)
1535	end do
1536
1537	CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax), 0.)
1538	CALL addpnt(x2,y2,kdata,n2, 0., 0.)
1539	CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
1540	CALL addpnt(x2,y2,kdata,n2, 1.e+38, 0.)
1541	CALL inter2(nw,wl,yg2,n2,x2,y2,ierr)
1542
1543	fil = trim(datadir)//'/cross_sections/no2_xs_vandaele_220K.txt'
1544	print*, 'section efficace NO2: ', fil
1545
1546	OPEN(UNIT=kin,FILE=fil,status='old')
1547	DO i = 1, 3
1548	READ(kin,*)
1549	ENDDO
1550	n3 = 27993
1551	DO i = 1, n3
1552	READ(kin,*) x3(i), y3(i)
1553	end do
1554
1555	CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax), 0.)
1556	CALL addpnt(x3,y3,kdata,n3, 0., 0.)
1557	CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.)
1558	CALL addpnt(x3,y3,kdata,n3, 1.e+38, 0.)
1559	CALL inter2(nw,wl,yg3,n3,x3,y3,ierr)
1560
1561	do iw = 1, nw - 1
1562	xsno2(iw) = yg1(iw)
1563	xsno2_294(iw) = yg2(iw)
1564	xsno2_220(iw) = yg3(iw)
1565	end do
1566
1567	! photodissociation efficiency from jpl 2006
1568
1569	fil = trim(datadir)//'/cross_sections/no2_yield_jpl2006.txt'
1570	print*, 'quantum yield NO2: ', fil
1571
1572	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1573	DO i = 1, 5
1574	READ(kin,*)
1575	ENDDO
1576	n = 25
1577	n4 = n
1578	n5 = n
1579	DO i = 1, n
1580	READ(kin,*) x4(i), y4(i), y5(i)
1581	x5(i) = x4(i)
1582	ENDDO
1583	CLOSE(kin)
1584
1585	CALL addpnt(x4,y4,kdata,n4,x4(1)*(1.-deltax),y4(1))
1586	CALL addpnt(x4,y4,kdata,n4, 0.,y4(1))
1587	CALL addpnt(x4,y4,kdata,n4,x4(n4)*(1.+deltax), 0.)
1588	CALL addpnt(x4,y4,kdata,n4, 1.e+38, 0.)
1589	CALL inter2(nw,wl,yg4,n4,x4,y4,ierr)
1590	IF (ierr .NE. 0) THEN
1591	WRITE(,) ierr, fil
1592	STOP
1593	ENDIF
1594
1595	CALL addpnt(x5,y5,kdata,n5,x5(1)*(1.-deltax),y5(1))
1596	CALL addpnt(x5,y5,kdata,n5, 0.,y5(1))
1597	CALL addpnt(x5,y5,kdata,n5,x5(n5)*(1.+deltax), 0.)
1598	CALL addpnt(x5,y5,kdata,n5, 1.e+38, 0.)
1599	CALL inter2(nw,wl,yg5,n5,x5,y5,ierr)
1600	IF (ierr .NE. 0) THEN
1601	WRITE(,) ierr, fil
1602	STOP
1603	ENDIF
1604
1605	do iw = 1, nw - 1
1606	yldno2_298(iw) = yg4(iw)
1607	yldno2_248(iw) = yg5(iw)
1608	end do
1609
1610	end subroutine rdxsno2
1611
1612	!==============================================================================
1613
1614	subroutine rdxsno(nw, wl, yg, yieldno)
1615
1616	!-----------------------------------------------------------------------------*
1617	!= PURPOSE: =*
1618	!= Read NO cross-sections and photodissociation efficiency =*
1619	!= Lida et al 1986 (provided by Francisco Gonzalez-Galindo) =*
1620	!-----------------------------------------------------------------------------*
1621	!= PARAMETERS: =*
1622	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
1623	!= wavelength grid =*
1624	!-----------------------------------------------------------------------------*
1625
1626	use datafile_mod, only: datadir
1627
1628	implicit none
1629
1630	! input
1631
1632	integer :: nw ! number of wavelength grid points
1633	real, dimension(nw) :: wl ! lower wavelength for each interval
1634
1635	! output
1636
1637	real, dimension(nw) :: yg ! no cross-sections (cm2)
1638	real, dimension(nw) :: yieldno ! no photodissociation efficiency
1639
1640	! local
1641
1642	integer, parameter :: kdata = 110
1643	real, parameter :: deltax = 1.e-4
1644	real, dimension(kdata) :: x1, y1, x2, y2
1645	integer :: i, iw, n, ierr
1646	character*100 fil
1647	integer :: kin, kout ! input/output logical units
1648
1649	kin = 10
1650
1651	! no cross-sections
1652
1653	fil = trim(datadir)//'/cross_sections/no_xs_francisco.txt'
1654	print*, 'section efficace NO: ', fil
1655	OPEN(kin,FILE=fil,STATUS='OLD')
1656
1657	n = 99
1658	DO i = 1, n
1659	READ(kin,*) x1(i), y1(i)
1660	ENDDO
1661	CLOSE(kin)
1662
1663	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1664	CALL addpnt(x1,y1,kdata,n, 0.,0.)
1665	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1666	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
1667
1668	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
1669	IF (ierr .NE. 0) THEN
1670	WRITE(,) ierr, fil
1671	STOP
1672	ENDIF
1673
1674	! photodissociation yield
1675
1676	fil = trim(datadir)//'/cross_sections/noefdis.txt'
1677	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1678
1679	n = 33
1680	DO i = 1, n
1681	READ(kin,*) x2(n-i+1), y2(n-i+1)
1682	END DO
1683	CLOSE (kin)
1684
1685	CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
1686	CALL addpnt(x2,y2,kdata,n, 0.,0.)
1687	CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
1688	CALL addpnt(x2,y2,kdata,n, 1.e+38,1.)
1689	CALL inter2(nw,wl,yieldno,n,x2,y2,ierr)
1690	IF (ierr .NE. 0) THEN
1691	WRITE(,) ierr, fil
1692	STOP
1693	ENDIF
1694
1695	end subroutine rdxsno
1696
1697	!==============================================================================
1698
1699	subroutine rdxsn2(nw, wl, yg, yieldn2)
1700
1701	!-----------------------------------------------------------------------------*
1702	!= PURPOSE: =*
1703	!= Read n2 cross-sections and photodissociation yield =*
1704	!-----------------------------------------------------------------------------*
1705	!= PARAMETERS: =*
1706	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
1707	!= wavelength grid =*
1708	!-----------------------------------------------------------------------------*
1709
1710	use datafile_mod, only: datadir
1711
1712	implicit none
1713
1714	! input
1715
1716	integer :: nw ! number of wavelength grid points
1717	real, dimension(nw) :: wl ! lower wavelength for each interval
1718
1719	! output
1720
1721	real, dimension(nw) :: yg ! n2 cross-sections (cm2)
1722	real, dimension(nw) :: yieldn2 ! n2 photodissociation yield
1723
1724	! local
1725
1726	integer, parameter :: kdata = 1100
1727	real, parameter :: deltax = 1.e-4
1728	real, dimension(kdata) :: x1, y1, x2, y2
1729	real :: xl, xu
1730	integer :: i, iw, n, ierr
1731	integer :: kin, kout ! input/output logical units
1732	character*100 fil
1733
1734	kin = 10
1735
1736	! n2 cross sections
1737
1738	fil = trim(datadir)//'/cross_sections/n2secef_01nm.txt'
1739	print*, 'section efficace N2: ', fil
1740	OPEN(kin,FILE=fil,STATUS='OLD')
1741
1742	n = 1020
1743	DO i = 1, n
1744	READ(kin,*) x1(i), y1(i)
1745	ENDDO
1746	CLOSE(kin)
1747
1748	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1749	CALL addpnt(x1,y1,kdata,n, 0.,0.)
1750	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1751	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
1752
1753	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
1754
1755	IF (ierr .NE. 0) THEN
1756	WRITE(,) ierr, fil
1757	STOP
1758	ENDIF
1759
1760	! photodissociation yield
1761
1762	fil = trim(datadir)//'/cross_sections/n2_ionef_schunknagy2000.txt'
1763	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1764
1765	n = 19
1766	read(kin,*)
1767	DO i = 1, n
1768	READ(kin,*) xl, xu, y2(i)
1769	x2(i) = (xl + xu)/2.
1770	y2(i) = 1. - y2(i)
1771	END DO
1772	CLOSE (kin)
1773
1774	CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
1775	CALL addpnt(x2,y2,kdata,n, 0.,0.)
1776	CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
1777	CALL addpnt(x2,y2,kdata,n, 1.e+38,1.)
1778	CALL inter2(nw,wl,yieldn2,n,x2,y2,ierr)
1779	IF (ierr .NE. 0) THEN
1780	WRITE(,) ierr, fil
1781	STOP
1782	ENDIF
1783
1784	end subroutine rdxsn2
1785
1786	!==============================================================================
1787
1788	subroutine setalb(nw,wl,albedo)
1789
1790	!-----------------------------------------------------------------------------*
1791	!= PURPOSE: =*
1792	!= Set the albedo of the surface. The albedo is assumed to be Lambertian, =*
1793	!= i.e., the reflected light is isotropic, and idependt of the direction =*
1794	!= of incidence of light. Albedo can be chosen to be wavelength dependent. =*
1795	!-----------------------------------------------------------------------------*
1796	!= PARAMETERS: =*
1797	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
1798	!= wavelength grid =*
1799	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
1800	!= working wavelength grid =*
1801	!= ALBEDO - REAL, surface albedo at each specified wavelength (O)=*
1802	!-----------------------------------------------------------------------------*
1803
1804	implicit none
1805
1806	! input: (wavelength working grid data)
1807
1808	INTEGER nw
1809	REAL wl(nw)
1810
1811	! output:
1812
1813	REAL albedo(nw)
1814
1815	! local:
1816
1817	INTEGER iw
1818	REAL alb
1819
1820	! 0.015: mean value from clancy et al., icarus, 49-63, 1999.
1821
1822	alb = 0.015
1823
1824	do iw = 1, nw - 1
1825	albedo(iw) = alb
1826	end do
1827
1828	end subroutine setalb
1829
1830	end module photolysis_mod

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

Download in other formats:

Original Format