Context Navigation

photolysis_mod.F90 @ 2599

Last change on this file since 2599 was 2489, checked in by flefevre, 4 years ago
Photolysis : minor cosmetic changes
File size: 59.5 KB

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

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

Download in other formats:

Original Format