Context Navigation

photolysis_mod.F90 @ 2082

Last change on this file since 2082 was 2044, checked in by flefevre, 7 years ago
Photolyse on-line : desormais par defaut
File size: 57.2 KB

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

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

Download in other formats:

Original Format