Context Navigation

photolysis_mod.F90 @ 2833

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

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

Download in other formats:

Original Format