Context Navigation

photolysis_mod.F90 @ 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: 56.9 KB

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

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

Download in other formats:

Original Format