Context Navigation

photolysis_mod.F90 @ 3461

Last change on this file since 3461 was 2925, checked in by streel, 20 months ago

VENUS PCM:

Update of S2 photolysis : now computed within the online photolysis
Update of the reaction N + O + CO2 -> NO + CO2 to take in account CO2 efficiency
New reaction : N(2d) + CO -> N + CO

File size: 97.2 KB

Line
1	module photolysis_mod
2
3	implicit none
4
5	! photolysis
6
7	integer, save :: nphot = 24 ! number of photolysis
8
9	!$OMP THREADPRIVATE(nphot)
10
11	integer, parameter :: nabs = 21 ! number of absorbing gases
12
13	! spectral grid
14
15	integer, parameter :: nw = 194 ! number of spectral intervals (low-res)
16	integer, save :: mopt ! high-res/low-res switch
17
18	!$OMP THREADPRIVATE(mopt)
19
20	real, dimension(nw), save :: wl, wc, wu ! lower, center, upper wavelength for each interval
21
22	!$OMP THREADPRIVATE(wl, wc, wu)
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 :: xshcl ! hcl absorption cross-section (cm2)
39	real, dimension(nw), save :: xscl2 ! cl2 absorption cross-section (cm2)
40	real, dimension(nw), save :: xshocl ! hocl absorption cross-section (cm2)
41	real, dimension(nw), save :: xsso2_200, xsso2_298, xsso2_360 ! so2 absorption cross-section (cm2)
42	real, dimension(nw), save :: xsso ! so absorption cross-section (cm2)
43	real, dimension(nw), save :: xsso3 ! so3 absorption cross-section (cm2)
44	real, dimension(nw), save :: xsclo ! clo absorption cross-section (cm2)
45	real, dimension(nw), save :: xsocs ! cos absorption cross-section (cm2)
46	real, dimension(nw), save :: xss2 ! s2 absorption cross-section (cm2)
47	real, dimension(nw), save :: xscocl2 ! cocl2 absorption cross-section (cm2)
48	real, dimension(nw), save :: xsh2so4 ! h2so4 absorption cross-section (cm2)
49	real, dimension(nw), save :: xsh2 ! h2 absorption cross-section (cm2)
50	real, dimension(nw), save :: yieldh2 ! h2 photodissociation yield
51	real, dimension(nw), save :: xsno2, xsno2_220, xsno2_294 ! no2 absorption cross-section at 220-294 k (cm2)
52	real, dimension(nw), save :: yldno2_248, yldno2_298 ! no2 quantum yield at 248-298 k
53	real, dimension(nw), save :: xsno ! no absorption cross-section (cm2)
54	real, dimension(nw), save :: yieldno ! no photodissociation yield
55	real, dimension(nw), save :: xsn2 ! n2 absorption cross-section (cm2)
56	real, dimension(nw), save :: yieldn2 ! n2 photodissociation yield
57	real, dimension(nw), save :: xshdo ! hdo absorption cross-section (cm2)
58	real, dimension(nw), save :: albedo ! surface albedo
59
60	!$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)
61	!$OMP THREADPRIVATE(xsno,yieldno,xsn2,yieldn2,xshdo,albedo)
62
63	contains
64
65	subroutine init_photolysis
66
67	! initialise on-line photolysis
68
69	! mopt = 1 high-resolution
70	! mopt = 2 martian low-resolution (recommended for gcm use)
71	! mopt = 3 venusian low-resolution (recommended for gcm use)
72
73	mopt = 3
74
75	! set wavelength grid
76
77	call gridw(nw,wl,wc,wu,mopt)
78
79	! read and grid solar flux data
80
81	call rdsolarflux(nw,wl,wc,f)
82
83	! read and grid o2 cross-sections
84
85	call rdxso2(nw,wl,xso2_150,xso2_200,xso2_250,xso2_300,yieldo2)
86
87	! read and grid co2 cross-sections
88
89	call rdxsco2(nw,wl,xsco2_195,xsco2_295,xsco2_370,yieldco2)
90
91	! read and grid o3 cross-sections
92
93	call rdxso3(nw,wl,xso3_218,xso3_298)
94
95	! read and grid h2o cross-sections
96
97	call rdxsh2o(nw,wl,xsh2o)
98
99	! read and grid h2o2 cross-sections
100
101	call rdxsh2o2(nw,wl,xsh2o2)
102
103	! read and grid ho2 cross-sections
104
105	call rdxsho2(nw,wl,xsho2)
106
107	! read and grid hcl cross-sections
108
109	call rdxshcl(nw,wl,xshcl)
110
111	! read and grid cl2 cross-sections
112
113	call rdxscl2(nw,wl,xscl2)
114
115	! read and grid hocl cross-sections
116
117	call rdxshocl(nw,wl,xshocl)
118
119	! read and grid so2 cross-sections
120
121	call rdxsso2(nw,wl,xsso2_200,xsso2_298,xsso2_360)
122
123	! read and grid so cross-sections
124
125	call rdxsso(nw,wl,xsso)
126
127	! read and grid so3 cross-sections
128
129	call rdxsso3(nw,wl,xsso3)
130
131	! read and grid clo cross-sections
132
133	call rdxsclo(nw,wl,xsclo)
134
135	! read and grid ocs cross-sections
136
137	call rdxsocs(nw,wl,xsocs)
138
139	! read and grid cocl2 cross-sections
140
141	call rdxscocl2(nw,wl,xscocl2)
142
143	! read and grid s2 cross-sections
144
145	call rdxss2(nw,wl,xss2)
146
147	! read and grid h2so4 cross-sections
148
149	call rdxsh2so4(nw,wl,xsh2so4)
150
151	! read and grid h2 cross-sections
152
153	call rdxsh2(nw,wl,wc,xsh2,yieldh2)
154
155	! read and grid no cross-sections
156
157	call rdxsno(nw,wl,xsno,yieldno)
158
159	! read and grid no2 cross-sections
160
161	call rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294,yldno2_248,yldno2_298)
162
163	! read and grid n2 cross-sections
164
165	call rdxsn2(nw,wl,xsn2,yieldn2)
166
167	! read and grid hdo cross-sections
168
169	call rdxshdo(nw,wl,xshdo)
170
171	! set surface albedo
172
173	call setalb(nw,wl,albedo)
174
175	end subroutine init_photolysis
176
177	!==============================================================================
178
179	subroutine gridw(nw,wl,wc,wu,mopt)
180
181	! Create the wavelength grid for all interpolations and radiative transfer
182	! calculations. Grid may be irregularly spaced. Wavelengths are in nm.
183	! No gaps are allowed within the wavelength grid.
184
185	implicit none
186
187	! input
188
189	integer :: nw ! number of wavelength grid points
190	integer :: mopt ! high-res/low-res switch
191
192	! output
193
194	real, dimension(nw) :: wl, wc, wu ! lower, center, upper wavelength for each interval
195
196	! local
197
198	real :: wincr ! wavelength increment
199	integer :: iw, kw
200
201	! mopt = 1 high-resolution mode (3789 intervals)
202	!
203	! 0-108 nm : 1.0 nm
204	! 108-124 nm : 0.1 nm
205	! 124-175 nm : 0.5 nm
206	! 175-205 nm : 0.01 nm
207	! 205-365 nm : 0.5 nm
208	! 365-850 nm : 5.0 nm
209	!
210	! mopt = 2 low-resolution mode (162)
211	!
212	! 0-60 nm : 6.0 nm
213	! 60-80 nm : 2.0 nm
214	! 80-85 nm : 5.0 nm
215	! 85-117 nm : 2.0 nm
216	! 117-120 nm : 5.0 nm
217	! 120-123 nm : 0.2 nm
218	! 123-163 nm : 5.0 nm
219	! 163-175 nm : 2.0 nm
220	! 175-205 nm : 0.5 nm
221	! 205-245 nm : 5.0 nm
222	! 245-415 nm : 10.0 nm
223	! 415-815 nm : 50.0 nm
224	!
225	! mopt = 3 venusian low-resolution mode (194)
226	! (205-245 nm -> 1nm of resolution)
227	!
228	! 0-60 nm : 6.0 nm
229	! 60-80 nm : 2.0 nm
230	! 80-85 nm : 5.0 nm
231	! 85-117 nm : 2.0 nm
232	! 117-120 nm : 5.0 nm
233	! 120-123 nm : 0.2 nm
234	! 123-163 nm : 5.0 nm
235	! 163-175 nm : 2.0 nm
236	! 175-205 nm : 0.5 nm
237	! 205-245 nm : 1.0 nm
238	! 245-415 nm : 10.0 nm
239	! 415-815 nm : 50.0 nm
240
241	if (mopt == 1) then ! high-res
242
243	! define wavelength intervals of width 1.0 nm from 0 to 108 nm:
244
245	kw = 0
246	wincr = 1.0
247	do iw = 0, 107
248	kw = kw + 1
249	wl(kw) = real(iw)
250	wu(kw) = wl(kw) + wincr
251	wc(kw) = (wl(kw) + wu(kw))/2.
252	end do
253
254	! define wavelength intervals of width 0.1 nm from 108 to 124 nm:
255
256	wincr = 0.1
257	do iw = 1080, 1239, 1
258	kw = kw + 1
259	wl(kw) = real(iw)/10.
260	wu(kw) = wl(kw) + wincr
261	wc(kw) = (wl(kw) + wu(kw))/2.
262	end do
263
264	! define wavelength intervals of width 0.5 nm from 124 to 175 nm:
265
266	wincr = 0.5
267	do iw = 1240, 1745, 5
268	kw = kw + 1
269	wl(kw) = real(iw)/10.
270	wu(kw) = wl(kw) + wincr
271	wc(kw) = (wl(kw) + wu(kw))/2.
272	end do
273
274	! define wavelength intervals of width 0.01 nm from 175 to 205 nm:
275
276	wincr = 0.01
277	do iw = 17500, 20499, 1
278	kw = kw + 1
279	wl(kw) = real(iw)/100.
280	wu(kw) = wl(kw) + wincr
281	wc(kw) = (wl(kw) + wu(kw))/2.
282	end do
283
284	! define wavelength intervals of width 0.5 nm from 205 to 365 nm:
285
286	wincr = 0.5
287	do iw = 2050, 3645, 5
288	kw = kw + 1
289	wl(kw) = real(iw)/10.
290	wu(kw) = wl(kw) + wincr
291	wc(kw) = (wl(kw) + wu(kw))/2.
292	end do
293
294	! define wavelength intervals of width 5.0 nm from 365 to 855 nm:
295
296	wincr = 5.0
297	do iw = 365, 850, 5
298	kw = kw + 1
299	wl(kw) = real(iw)
300	wu(kw) = wl(kw) + wincr
301	wc(kw) = (wl(kw) + wu(kw))/2.
302	end do
303	wl(kw+1) = wu(kw)
304
305	!============================================================
306
307	else if (mopt == 2) then ! low-res
308
309	! define wavelength intervals of width 6.0 nm from 0 to 60 nm:
310
311	kw = 0
312	wincr = 6.0
313	DO iw = 0, 54, 6
314	kw = kw + 1
315	wl(kw) = real(iw)
316	wu(kw) = wl(kw) + wincr
317	wc(kw) = (wl(kw) + wu(kw))/2.
318	END DO
319
320	! define wavelength intervals of width 2.0 nm from 60 to 80 nm:
321
322	wincr = 2.0
323	DO iw = 60, 78, 2
324	kw = kw + 1
325	wl(kw) = real(iw)
326	wu(kw) = wl(kw) + wincr
327	wc(kw) = (wl(kw) + wu(kw))/2.
328	END DO
329
330	! define wavelength intervals of width 5.0 nm from 80 to 85 nm:
331
332	wincr = 5.0
333	DO iw = 80, 80
334	kw = kw + 1
335	wl(kw) = real(iw)
336	wu(kw) = wl(kw) + wincr
337	wc(kw) = (wl(kw) + wu(kw))/2.
338	END DO
339
340	! define wavelength intervals of width 2.0 nm from 85 to 117 nm:
341
342	wincr = 2.0
343	DO iw = 85, 115, 2
344	kw = kw + 1
345	wl(kw) = real(iw)
346	wu(kw) = wl(kw) + wincr
347	wc(kw) = (wl(kw) + wu(kw))/2.
348	END DO
349
350	! define wavelength intervals of width 3.0 nm from 117 to 120 nm:
351
352	wincr = 3.0
353	DO iw = 117, 117
354	kw = kw + 1
355	wl(kw) = real(iw)
356	wu(kw) = wl(kw) + wincr
357	wc(kw) = (wl(kw) + wu(kw))/2.
358	END DO
359
360	! define wavelength intervals of width 0.2 nm from 120 to 123 nm:
361
362	wincr = 0.2
363	DO iw = 1200, 1228, 2
364	kw = kw + 1
365	wl(kw) = real(iw)/10.
366	wu(kw) = wl(kw) + wincr
367	wc(kw) = (wl(kw) + wu(kw))/2.
368	ENDDO
369
370	! define wavelength intervals of width 5.0 nm from 123 to 163 nm:
371
372	wincr = 5.0
373	DO iw = 123, 158, 5
374	kw = kw + 1
375	wl(kw) = real(iw)
376	wu(kw) = wl(kw) + wincr
377	wc(kw) = (wl(kw) + wu(kw))/2.
378	ENDDO
379
380	! define wavelength intervals of width 2.0 nm from 163 to 175 nm:
381
382	wincr = 2.0
383	DO iw = 163, 173, 2
384	kw = kw + 1
385	wl(kw) = real(iw)
386	wu(kw) = wl(kw) + wincr
387	wc(kw) = (wl(kw) + wu(kw))/2.
388	ENDDO
389
390	! define wavelength intervals of width 0.5 nm from 175 to 205 nm:
391
392	wincr = 0.5
393	DO iw = 1750, 2045, 5
394	kw = kw + 1
395	wl(kw) = real(iw)/10.
396	wu(kw) = wl(kw) + wincr
397	wc(kw) = (wl(kw) + wu(kw))/2.
398	ENDDO
399
400	! define wavelength intervals of width 5.0 nm from 205 to 245 nm:
401
402	wincr = 5.
403	DO iw = 205, 240, 5
404	kw = kw + 1
405	wl(kw) = real(iw)
406	wu(kw) = wl(kw) + wincr
407	wc(kw) = (wl(kw) + wu(kw))/2.
408	ENDDO
409
410	! define wavelength intervals of width 10.0 nm from 245 to 415 nm:
411
412	wincr = 10.0
413	DO iw = 245, 405, 10
414	kw = kw + 1
415	wl(kw) = real(iw)
416	wu(kw) = wl(kw) + wincr
417	wc(kw) = (wl(kw) + wu(kw))/2.
418	ENDDO
419
420	! define wavelength intervals of width 50.0 nm from 415 to 815 nm:
421
422	wincr = 50.0
423	DO iw = 415, 815, 50
424	kw = kw + 1
425	wl(kw) = real(iw)
426	wu(kw) = wl(kw) + wincr
427	wc(kw) = (wl(kw) + wu(kw))/2.
428	ENDDO
429
430	wl(kw+1) = wu(kw)
431
432	!============================================================
433
434	else if (mopt == 3) then ! low-res
435
436	! define wavelength intervals of width 6.0 nm from 0 to 60 nm:
437
438	kw = 0
439	wincr = 6.0
440	DO iw = 0, 54, 6
441	kw = kw + 1
442	wl(kw) = real(iw)
443	wu(kw) = wl(kw) + wincr
444	wc(kw) = (wl(kw) + wu(kw))/2.
445	END DO
446
447	! define wavelength intervals of width 2.0 nm from 60 to 80 nm:
448
449	wincr = 2.0
450	DO iw = 60, 78, 2
451	kw = kw + 1
452	wl(kw) = real(iw)
453	wu(kw) = wl(kw) + wincr
454	wc(kw) = (wl(kw) + wu(kw))/2.
455	END DO
456
457	! define wavelength intervals of width 5.0 nm from 80 to 85 nm:
458
459	wincr = 5.0
460	DO iw = 80, 80
461	kw = kw + 1
462	wl(kw) = real(iw)
463	wu(kw) = wl(kw) + wincr
464	wc(kw) = (wl(kw) + wu(kw))/2.
465	END DO
466
467	! define wavelength intervals of width 2.0 nm from 85 to 117 nm:
468
469	wincr = 2.0
470	DO iw = 85, 115, 2
471	kw = kw + 1
472	wl(kw) = real(iw)
473	wu(kw) = wl(kw) + wincr
474	wc(kw) = (wl(kw) + wu(kw))/2.
475	END DO
476
477	! define wavelength intervals of width 3.0 nm from 117 to 120 nm:
478
479	wincr = 3.0
480	DO iw = 117, 117
481	kw = kw + 1
482	wl(kw) = real(iw)
483	wu(kw) = wl(kw) + wincr
484	wc(kw) = (wl(kw) + wu(kw))/2.
485	END DO
486
487	! define wavelength intervals of width 0.2 nm from 120 to 123 nm:
488
489	wincr = 0.2
490	DO iw = 1200, 1228, 2
491	kw = kw + 1
492	wl(kw) = real(iw)/10.
493	wu(kw) = wl(kw) + wincr
494	wc(kw) = (wl(kw) + wu(kw))/2.
495	ENDDO
496
497	! define wavelength intervals of width 5.0 nm from 123 to 163 nm:
498
499	wincr = 5.0
500	DO iw = 123, 158, 5
501	kw = kw + 1
502	wl(kw) = real(iw)
503	wu(kw) = wl(kw) + wincr
504	wc(kw) = (wl(kw) + wu(kw))/2.
505	ENDDO
506
507	! define wavelength intervals of width 2.0 nm from 163 to 175 nm:
508
509	wincr = 2.0
510	DO iw = 163, 173, 2
511	kw = kw + 1
512	wl(kw) = real(iw)
513	wu(kw) = wl(kw) + wincr
514	wc(kw) = (wl(kw) + wu(kw))/2.
515	ENDDO
516
517	! define wavelength intervals of width 0.5 nm from 175 to 205 nm:
518
519	wincr = 0.5
520	DO iw = 1750, 2045, 5
521	kw = kw + 1
522	wl(kw) = real(iw)/10.
523	wu(kw) = wl(kw) + wincr
524	wc(kw) = (wl(kw) + wu(kw))/2.
525	ENDDO
526
527	! define wavelength intervals of width 5.0 nm from 205 to 245 nm:
528
529	wincr = 1
530	DO iw = 205, 244, 1
531	kw = kw + 1
532	wl(kw) = real(iw)
533	wu(kw) = wl(kw) + wincr
534	wc(kw) = (wl(kw) + wu(kw))/2.
535	ENDDO
536
537	! define wavelength intervals of width 10.0 nm from 245 to 415 nm:
538
539	wincr = 10.0
540	DO iw = 245, 405, 10
541	kw = kw + 1
542	wl(kw) = real(iw)
543	wu(kw) = wl(kw) + wincr
544	wc(kw) = (wl(kw) + wu(kw))/2.
545	ENDDO
546
547	! define wavelength intervals of width 50.0 nm from 415 to 815 nm:
548
549	wincr = 50.0
550	DO iw = 415, 815, 50
551	kw = kw + 1
552	wl(kw) = real(iw)
553	wu(kw) = wl(kw) + wincr
554	wc(kw) = (wl(kw) + wu(kw))/2.
555	ENDDO
556
557	wl(kw+1) = wu(kw)
558
559	print*, 'number of spectral intervals : ', kw+1
560
561	endif
562
563	end subroutine gridw
564
565	!==============================================================================
566
567	subroutine rdsolarflux(nw,wl,wc,f)
568
569	! Read and re-grid solar flux data.
570
571	USE mod_phys_lmdz_para, ONLY: is_master
572	USE mod_phys_lmdz_transfert_para, ONLY: bcast
573
574	implicit none
575
576	#include "clesphys.h" ! fixed_euv_value
577
578	! input
579
580	integer :: nw ! number of wavelength grid points
581	real, dimension(nw) :: wl, wc ! lower and central wavelength for each interval
582
583	! output
584
585	real, dimension(nw) :: f ! solar flux (w.m-2.nm-1)
586
587	! local
588
589	integer, parameter :: kdata = 20000 ! max dimension of input solar flux
590	integer :: msun ! choice of solar flux
591	integer :: iw, nhead, ihead, n, i, ierr, kin1, kin2
592
593	real, parameter :: deltax = 1.e-4
594	! atlas1_thuillier_tuv.txt
595	real, dimension(kdata) :: x1, y1 ! input solar flux - HIGH SOLAR ACTIVITY
596	real :: E107_max = 196
597	! atlas3_thuillier_tuv.txt
598	real, dimension(kdata) :: x2, y2 ! input solar flux - "LOW" SOLAR ACTIVITY
599	real :: E107_min = 97
600	! be careful, x2 need to be equal to x1
601	real, dimension(kdata) :: x0, y0 ! input solar flux - interpolated solar activity
602	real, dimension(nw) :: yg0 ! gridded solar flux
603	real :: factor_interp
604
605	character(len=100) :: fil1, fil2
606
607	kin1 = 10 ! input logical unit
608	kin2 = 11 ! input logical unit
609
610	! select desired extra-terrestrial solar irradiance, using msun:
611
612	! 18 = atlas1_thuillier_tuv.txt 0-900 nm March 29, 1992
613	! Article: F10.7 = 192 s.f.u \| <F10.7>81d = 171 s.f.u \| <Rz> = 121 (sunsport number)
614	! Thuillier et al., Adv. Space. Res., 34, 256-261, 2004
615	! Model SOLAR2000 version 2015/12: E10.7 = 195.8 s.f.u \| <E10.7>81d = 196.9 s.f.u
616	! Choix de la limite haute: E10.7 = 196 s.f.u
617
618	! 20 = atlas3_thuillier_tuv.txt 0-900 nm November 11, 1994
619	! Article: F10.7 = 77.5 s.f.u \| <F10.7>81d = 83.5 s.f.u \| <Rz> = 20 (sunsport number)
620	! Thuillier et al., Adv. Space. Res., 34, 256-261, 2004
621	! Model SOLAR2000 version 2015/12: E10.7 = 96.9 s.f.u \| <E10.7>81d = 100.0 s.f.u
622	! Choix de la limite basse: E10.7 = 97 s.f.u
623
624	if (fixed_euv_value .ge. 196) THEN
625	msun = 18
626	print*, 'Atlas1 spectrum Thuiller chosen'
627	else
628	msun = 19
629	print*, 'interpolated Spectrum case'
630	end if
631
632	IF (is_master) THEN
633
634	fil1 = 'solar_fluxes/atlas1_thuillier_tuv.txt'
635	print*, 'solar flux : ', fil1
636
637	open(kin1, file=fil1, status='old', iostat=ierr)
638
639	if (ierr /= 0) THEN
640	write(,)'cant find solar flux : ', fil1
641	write(,)'It should be in : INPUT/solar_fluxes'
642	write(,)'1) You can change this directory address in '
643	write(,)' callphys.def with datadir=/path/to/dir'
644	write(,)'2) If necessary, /solar fluxes (and other datafiles)'
645	write(,)' can be obtained online on:'
646	write(,)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
647	stop
648	end if
649
650	nhead = 9
651	n = 19193
652	DO ihead = 1, nhead
653	READ(kin1,*)
654	ENDDO
655	DO i = 1, n
656	READ(kin1,*) x1(i), y1(i)
657	y1(i) = y1(i)*1.e-3 ! mw -> w
658	ENDDO
659	CLOSE (kin1)
660
661	fil2 = 'solar_fluxes/atlas3_thuillier_tuv.txt'
662	print*, 'solar flux : ', fil2
663
664	open(kin2, file=fil2, status='old', iostat=ierr)
665
666	if (ierr /= 0) THEN
667	write(,)'cant find solar flux : ', fil2
668	write(,)'It should be in : INPUT/solar_fluxes'
669	write(,)'1) You can change this directory address in '
670	write(,)' callphys.def with datadir=/path/to/dir'
671	write(,)'2) If necessary, /solar fluxes (and other datafiles)'
672	write(,)' can be obtained online on:'
673	write(,)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
674	stop
675	end if
676
677	nhead = 9
678	n = 19193
679	DO ihead = 1, nhead
680	READ(kin2,*)
681	ENDDO
682	DO i = 1, n
683	READ(kin2,*) x2(i), y2(i)
684	y2(i) = y2(i)*1.e-3 ! mw -> w
685	ENDDO
686	CLOSE (kin2)
687
688	IF (msun .eq. 18) THEN
689	DO i = 1, n
690	x0(i) = x1(i)
691	y0(i) = y1(i)
692	ENDDO
693	ELSE
694	! INTERPOLATION BETWEEN E107_min to E107_max and extrapolation below E107_min
695	factor_interp=(fixed_euv_value-E107_min)/(E107_max-E107_min)
696	DO i = 1, n
697	x0(i) = x1(i)
698	y0(i) = y2(i) + (y1(i) - y2(i))* factor_interp
699	ENDDO
700	ENDIF ! msun .ne. 18
701
702	CALL addpnt(x0,y0,kdata,n,x0(1)*(1.-deltax),0.)
703	CALL addpnt(x0,y0,kdata,n, 0.,0.)
704	CALL addpnt(x0,y0,kdata,n,x0(n)*(1.+deltax),0.)
705	CALL addpnt(x0,y0,kdata,n, 1.e+38,0.)
706	CALL inter2(nw,wl,yg0,n,x0,y0,ierr)
707
708	IF (ierr .NE. 0) THEN
709	WRITE(,) ierr, fil1, fil2
710	STOP
711	ENDIF
712
713	! convert to photon.s-1.nm-1.cm-2
714	! 5.039e11 = 1.e-41e-9/(hc = 6.62e-342.998e8)
715
716	DO iw = 1, nw-1
717	f(iw) = yg0(iw)wc(iw)5.039e11
718	ENDDO
719
720	endif !is_master
721
722	call bcast(f)
723
724	end subroutine rdsolarflux
725
726	!==============================================================================
727
728	subroutine addpnt ( x, y, ld, n, xnew, ynew )
729
730	!-----------------------------------------------------------------------------*
731	!= PURPOSE: =*
732	!= Add a point <xnew,ynew> to a set of data pairs <x,y>. x must be in =*
733	!= ascending order =*
734	!-----------------------------------------------------------------------------*
735	!= PARAMETERS: =*
736	!= X - REAL vector of length LD, x-coordinates (IO)=*
737	!= Y - REAL vector of length LD, y-values (IO)=*
738	!= LD - INTEGER, dimension of X, Y exactly as declared in the calling (I)=*
739	!= program =*
740	!= N - INTEGER, number of elements in X, Y. On entry, it must be: (IO)=*
741	!= N < LD. On exit, N is incremented by 1. =*
742	!= XNEW - REAL, x-coordinate at which point is to be added (I)=*
743	!= YNEW - REAL, y-value of point to be added (I)=*
744	!-----------------------------------------------------------------------------*
745
746	IMPLICIT NONE
747
748	! calling parameters
749
750	INTEGER ld, n
751	REAL x(ld), y(ld)
752	REAL xnew, ynew
753	INTEGER ierr
754
755	! local variables
756
757	INTEGER insert
758	INTEGER i
759
760	!-----------------------------------------------------------------------
761
762	! initialize error flag
763
764	ierr = 0
765
766	! check n<ld to make sure x will hold another point
767
768	IF (n .GE. ld) THEN
769	WRITE(0,*) '>>> ERROR (ADDPNT) <<< Cannot expand array '
770	WRITE(0,*) ' All elements used.'
771	STOP
772	ENDIF
773
774	insert = 1
775	i = 2
776
777	! check, whether x is already sorted.
778	! also, use this loop to find the point at which xnew needs to be inserted
779	! into vector x, if x is sorted.
780
781	10 CONTINUE
782	IF (i .LT. n) THEN
783	IF (x(i) .LT. x(i-1)) THEN
784	print*, x(i-1), x(i)
785	WRITE(0,*) '>>> ERROR (ADDPNT) <<< x-data must be in ascending order!'
786	STOP
787	ELSE
788	IF (xnew .GT. x(i)) insert = i + 1
789	ENDIF
790	i = i+1
791	GOTO 10
792	ENDIF
793
794	! if <xnew,ynew> needs to be appended at the end, just do so,
795	! otherwise, insert <xnew,ynew> at position INSERT
796
797	IF ( xnew .GT. x(n) ) THEN
798
799	x(n+1) = xnew
800	y(n+1) = ynew
801
802	ELSE
803
804	! shift all existing points one index up
805
806	DO i = n, insert, -1
807	x(i+1) = x(i)
808	y(i+1) = y(i)
809	ENDDO
810
811	! insert new point
812
813	x(insert) = xnew
814	y(insert) = ynew
815
816	ENDIF
817
818	! increase total number of elements in x, y
819
820	n = n+1
821
822	end subroutine addpnt
823
824	!==============================================================================
825
826	subroutine inter2(ng,xg,yg,n,x,y,ierr)
827
828	!-----------------------------------------------------------------------------*
829	!= PURPOSE: =*
830	!= Map input data given on single, discrete points onto a set of target =*
831	!= bins. =*
832	!= The original input data are given on single, discrete points of an =*
833	!= arbitrary grid and are being linearly interpolated onto a specified set =*
834	!= of target bins. In general, this is the case for most of the weighting =*
835	!= functions (action spectra, molecular cross section, and quantum yield =*
836	!= data), which have to be matched onto the specified wavelength intervals. =*
837	!= The average value in each target bin is found by averaging the trapezoi- =*
838	!= dal area underneath the input data curve (constructed by linearly connec-=*
839	!= ting the discrete input values). =*
840	!= Some caution should be used near the endpoints of the grids. If the =*
841	!= input data set does not span the range of the target grid, an error =*
842	!= message is printed and the execution is stopped, as extrapolation of the =*
843	!= data is not permitted. =*
844	!= If the input data does not encompass the target grid, use ADDPNT to =*
845	!= expand the input array. =*
846	!-----------------------------------------------------------------------------*
847	!= PARAMETERS: =*
848	!= NG - INTEGER, number of bins + 1 in the target grid (I)=*
849	!= XG - REAL, target grid (e.g., wavelength grid); bin i is defined (I)=*
850	!= as [XG(i),XG(i+1)] (i = 1..NG-1) =*
851	!= YG - REAL, y-data re-gridded onto XG, YG(i) specifies the value for (O)=*
852	!= bin i (i = 1..NG-1) =*
853	!= N - INTEGER, number of points in input grid (I)=*
854	!= X - REAL, grid on which input data are defined (I)=*
855	!= Y - REAL, input y-data (I)=*
856	!-----------------------------------------------------------------------------*
857
858	IMPLICIT NONE
859
860	! input:
861	INTEGER ng, n
862	REAL x(n), y(n), xg(ng)
863
864	! output:
865	REAL yg(ng)
866
867	! local:
868	REAL area, xgl, xgu
869	REAL darea, slope
870	REAL a1, a2, b1, b2
871	INTEGER ngintv
872	INTEGER i, k, jstart
873	INTEGER ierr
874	!_______________________________________________________________________
875
876	ierr = 0
877
878	! test for correct ordering of data, by increasing value of x
879
880	DO 10, i = 2, n
881	IF (x(i) .LE. x(i-1)) THEN
882	ierr = 1
883	WRITE(,)'data not sorted'
884	WRITE(,) x(i), x(i-1)
885	RETURN
886	ENDIF
887	10 CONTINUE
888
889	DO i = 2, ng
890	IF (xg(i) .LE. xg(i-1)) THEN
891	ierr = 2
892	WRITE(0,*) '>>> ERROR (inter2) <<< xg-grid not sorted!'
893	RETURN
894	ENDIF
895	ENDDO
896
897	! check for xg-values outside the x-range
898
899	IF ( (x(1) .GT. xg(1)) .OR. (x(n) .LT. xg(ng)) ) THEN
900	WRITE(0,*) '>>> ERROR (inter2) <<< Data do not span grid. '
901	WRITE(0,*) ' Use ADDPNT to expand data and re-run.'
902	STOP
903	ENDIF
904
905	! find the integral of each grid interval and use this to
906	! calculate the average y value for the interval
907	! xgl and xgu are the lower and upper limits of the grid interval
908
909	jstart = 1
910	ngintv = ng - 1
911	DO 50, i = 1,ngintv
912
913	! initialize:
914
915	area = 0.0
916	xgl = xg(i)
917	xgu = xg(i+1)
918
919	! discard data before the first grid interval and after the
920	! last grid interval
921	! for internal grid intervals, start calculating area by interpolating
922	! between the last point which lies in the previous interval and the
923	! first point inside the current interval
924
925	k = jstart
926	IF (k .LE. n-1) THEN
927
928	! if both points are before the first grid, go to the next point
929	30 CONTINUE
930	IF (x(k+1) .LE. xgl) THEN
931	jstart = k - 1
932	k = k+1
933	IF (k .LE. n-1) GO TO 30
934	ENDIF
935
936
937	! if the last point is beyond the end of the grid, complete and go to the next
938	! grid
939	40 CONTINUE
940	IF ((k .LE. n-1) .AND. (x(k) .LT. xgu)) THEN
941
942	jstart = k-1
943
944	! compute x-coordinates of increment
945
946	a1 = MAX(x(k),xgl)
947	a2 = MIN(x(k+1),xgu)
948
949	! if points coincide, contribution is zero
950
951	IF (x(k+1).EQ.x(k)) THEN
952	darea = 0.e0
953	ELSE
954	slope = (y(k+1) - y(k))/(x(k+1) - x(k))
955	b1 = y(k) + slope*(a1 - x(k))
956	b2 = y(k) + slope*(a2 - x(k))
957	darea = (a2 - a1)*(b2 + b1)/2.
958	ENDIF
959
960	! find the area under the trapezoid from a1 to a2
961
962	area = area + darea
963
964	! go to next point
965
966	k = k+1
967	GO TO 40
968
969	ENDIF
970	ENDIF
971
972	! calculate the average y after summing the areas in the interval
973
974	yg(i) = area/(xgu - xgl)
975
976	50 CONTINUE
977
978	end subroutine inter2
979
980	!==============================================================================
981
982	subroutine rdxsco2(nw,wl,xsco2_195,xsco2_295,xsco2_370,yieldco2)
983
984	!-----------------------------------------------------------------------------*
985	!= PURPOSE: =*
986	!= Read and grid CO2 absorption cross-sections and photodissociation yield =*
987	!-----------------------------------------------------------------------------*
988	!= PARAMETERS: =*
989	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
990	!= wavelength grid =*
991	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
992	!= working wavelength grid =*
993	!= XSCO2 - REAL, molecular absoprtion cross section (cm^2) of CO2 at (O)=*
994	!= each specified wavelength =*
995	!-----------------------------------------------------------------------------*
996
997	USE mod_phys_lmdz_para, ONLY: is_master
998	USE mod_phys_lmdz_transfert_para, ONLY: bcast
999
1000	implicit none
1001
1002	! input
1003
1004	integer :: nw ! number of wavelength grid points
1005	real, dimension(nw) :: wl ! lower and central wavelength for each interval
1006
1007	! output
1008
1009	real, dimension(nw) :: xsco2_195, xsco2_295, xsco2_370 ! co2 cross-sections (cm2)
1010	real, dimension(nw) :: yieldco2 ! co2 photodissociation yield
1011
1012	! local
1013
1014	integer, parameter :: kdata = 42000
1015	real, parameter :: deltax = 1.e-4
1016	real, dimension(kdata) :: x1, y1, y2, y3, xion, ion
1017	real, dimension(nw) :: yg
1018	real :: xl, xu
1019	integer :: ierr, i, l, n, n1, n2, n3, n4
1020	CHARACTER*100 fil
1021
1022	integer :: kin, kout ! input/ouput logical units
1023
1024	kin = 10
1025	kout = 30
1026
1027	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1028	!
1029	! CO2 absorption cross-sections
1030	!
1031	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1032	!
1033	! 195K: huestis and berkowitz (2010) + Starck et al. (2006)
1034	! + Yoshino et al. (1996) + Parkinson et al. (2003) + extrapolation
1035	!
1036	! 295K: huestis and berkowitz (2010) + Starck et al. (2006)
1037	! + Yoshino et al. (1996) + Parkinson et al. (2003) + extrapolation
1038	!
1039	! 370K: huestis and berkowitz (2010) + Starck et al. (2006)
1040	! + Lewis and Carver (1983) + extrapolation
1041	!
1042	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1043
1044	n1 = 40769
1045	n2 = 41586
1046	n3 = 10110
1047
1048	! 195K:
1049
1050	fil = 'cross_sections/co2_euv_uv_2018_195k.txt'
1051	print*, 'section efficace CO2 195K: ', fil
1052
1053	if(is_master) then
1054
1055	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1056	DO i = 1,11
1057	read(kin,*)
1058	END DO
1059
1060	DO i = 1, n1
1061	READ(kin,*) x1(i), y1(i)
1062	END DO
1063	CLOSE (kin)
1064
1065	CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
1066	CALL addpnt(x1,y1,kdata,n1, 0.,0.)
1067	CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
1068	CALL addpnt(x1,y1,kdata,n1, 1.e+38,0.)
1069	CALL inter2(nw,wl,yg,n1,x1,y1,ierr)
1070	IF (ierr .NE. 0) THEN
1071	WRITE(,) ierr, fil
1072	STOP
1073	ENDIF
1074
1075	DO l = 1, nw-1
1076	xsco2_195(l) = yg(l)
1077	END DO
1078
1079	endif !is_master
1080
1081	call bcast(xsco2_195)
1082
1083	! 295K:
1084
1085	fil = 'cross_sections/co2_euv_uv_2018_295k.txt'
1086	print*, 'section efficace CO2 295K: ', fil
1087
1088	if(is_master) then
1089
1090	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1091	DO i = 1,11
1092	read(kin,*)
1093	END DO
1094
1095	DO i = 1, n2
1096	READ(kin,*) x1(i), y1(i)
1097	END DO
1098	CLOSE (kin)
1099
1100	CALL addpnt(x1,y1,kdata,n2,x1(1)*(1.-deltax),0.)
1101	CALL addpnt(x1,y1,kdata,n2, 0.,0.)
1102	CALL addpnt(x1,y1,kdata,n2,x1(n2)*(1.+deltax),0.)
1103	CALL addpnt(x1,y1,kdata,n2, 1.e+38,0.)
1104	CALL inter2(nw,wl,yg,n2,x1,y1,ierr)
1105	IF (ierr .NE. 0) THEN
1106	WRITE(,) ierr, fil
1107	STOP
1108	ENDIF
1109
1110	DO l = 1, nw-1
1111	xsco2_295(l) = yg(l)
1112	END DO
1113
1114	endif !is_master
1115
1116	call bcast(xsco2_295)
1117
1118	! 370K:
1119
1120	fil = 'cross_sections/co2_euv_uv_2018_370k.txt'
1121	print*, 'section efficace CO2 370K: ', fil
1122
1123	if(is_master) then
1124
1125	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1126	DO i = 1,11
1127	read(kin,*)
1128	END DO
1129
1130	DO i = 1, n3
1131	READ(kin,*) x1(i), y1(i)
1132	END DO
1133	CLOSE (kin)
1134
1135	CALL addpnt(x1,y1,kdata,n3,x1(1)*(1.-deltax),0.)
1136	CALL addpnt(x1,y1,kdata,n3, 0.,0.)
1137	CALL addpnt(x1,y1,kdata,n3,x1(n3)*(1.+deltax),0.)
1138	CALL addpnt(x1,y1,kdata,n3, 1.e+38,0.)
1139	CALL inter2(nw,wl,yg,n3,x1,y1,ierr)
1140	IF (ierr .NE. 0) THEN
1141	WRITE(,) ierr, fil
1142	STOP
1143	ENDIF
1144
1145	DO l = 1, nw-1
1146	xsco2_370(l) = yg(l)
1147	END DO
1148
1149	endif !is_master
1150
1151	call bcast(xsco2_370)
1152
1153	! photodissociation yield:
1154
1155	fil = 'cross_sections/efdis_co2-o2_schunkandnagy2000.txt'
1156	print*, 'photodissociation yield CO2: ', fil
1157
1158	if(is_master) then
1159
1160	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1161
1162	do i = 1,3
1163	read(kin,*)
1164	end do
1165
1166	n4 = 17
1167	do i = 1, n4
1168	read(kin,*) xl, xu, ion(i)
1169	xion(i) = (xl + xu)/2.
1170	ion(i) = max(ion(i), 0.)
1171	end do
1172	close(kin)
1173
1174	CALL addpnt(xion,ion,kdata,n4,xion(1)*(1.-deltax),0.)
1175	CALL addpnt(xion,ion,kdata,n4, 0.,0.)
1176	CALL addpnt(xion,ion,kdata,n4,xion(n4)*(1.+deltax),1.)
1177	CALL addpnt(xion,ion,kdata,n4, 1.e+38,1.)
1178	CALL inter2(nw,wl,yieldco2,n4,xion,ion,ierr)
1179	IF (ierr .NE. 0) THEN
1180	WRITE(,) ierr, fil
1181	STOP
1182	ENDIF
1183
1184	endif !is_master
1185
1186	call bcast(yieldco2)
1187
1188	! DO l = 1, nw-1
1189	! write(kout,*) wl(l), xsco2_195(l),
1190	! $ xsco2_295(l),
1191	! $ xsco2_370(l),
1192	! $ yieldco2(l)
1193	! END DO
1194
1195	end subroutine rdxsco2
1196
1197	!==============================================================================
1198
1199	subroutine rdxso2(nw,wl,xso2_150,xso2_200,xso2_250,xso2_300,yieldo2)
1200
1201	!-----------------------------------------------------------------------------*
1202	!= PURPOSE: =*
1203	!= Read and grid O2 cross-sections and photodissociation yield =*
1204	!-----------------------------------------------------------------------------*
1205	!= PARAMETERS: =*
1206	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
1207	!= wavelength grid =*
1208	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
1209	!= working wavelength grid =*
1210	!= XSO2 - REAL, molecular absorption cross section =*
1211	!-----------------------------------------------------------------------------*
1212
1213	USE mod_phys_lmdz_para, ONLY: is_master
1214	USE mod_phys_lmdz_transfert_para, ONLY: bcast
1215
1216	implicit none
1217
1218	! input
1219
1220	integer :: nw ! number of wavelength grid points
1221	real, dimension(nw) :: wl ! lower and central wavelength for each interval
1222
1223	! output
1224
1225	real, dimension(nw) :: xso2_150, xso2_200, xso2_250, xso2_300 ! o2 cross-sections (cm2)
1226	real, dimension(nw) :: yieldo2 ! o2 photodissociation yield
1227
1228	! local
1229
1230	integer, parameter :: kdata = 18000
1231	real, parameter :: deltax = 1.e-4
1232	real, dimension(kdata) :: x1, y1, x2, y2, x3, y3, x4, y4
1233	real, dimension(kdata) :: xion, ion
1234	real :: factor, xl, xu, dummy
1235	integer :: i, ierr, n, n1, n2, n3, n4, nhead
1236	integer :: kin, kout ! input/output logical units
1237	character*100 fil
1238
1239	kin = 10
1240	kout = 30
1241
1242	! read o2 cross section data
1243
1244	nhead = 22
1245	n = 17434
1246
1247	fil = 'cross_sections/o2_composite_2018_150K.txt'
1248	print*, 'section efficace O2 150K: ', fil
1249
1250	if(is_master) then
1251
1252	open(kin, file=fil, status='old', iostat=ierr)
1253
1254	if (ierr /= 0) THEN
1255	write(,)'cant find O2 cross-sections : ', fil
1256	write(,)'It should be in :INPUT/cross_sections'
1257	write(,)'1) You can have to setup the link to the dir '
1258	write(,)'2) If necessary, /cross_sections (and other datafiles)'
1259	write(,)' can be obtained online on:'
1260	write(,)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
1261	stop
1262	end if
1263
1264	DO i = 1,nhead
1265	read(kin,*)
1266	END DO
1267
1268	n1 = n
1269	DO i = 1, n1
1270	READ(kin,*) x1(i), y1(i)
1271	END DO
1272	CLOSE (kin)
1273
1274	CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
1275	CALL addpnt(x1,y1,kdata,n1, 0.,0.)
1276	CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
1277	CALL addpnt(x1,y1,kdata,n1, 1.e+38,0.)
1278	CALL inter2(nw,wl,xso2_150,n1,x1,y1,ierr)
1279	IF (ierr .NE. 0) THEN
1280	WRITE(,) ierr, fil
1281	STOP
1282	ENDIF
1283
1284	endif !is_master
1285
1286	call bcast(xso2_150)
1287
1288	fil = 'cross_sections/o2_composite_2018_200K.txt'
1289	print*, 'section efficace O2 200K: ', fil
1290
1291	if(is_master) then
1292
1293	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1294
1295	DO i = 1,nhead
1296	read(kin,*)
1297	END DO
1298
1299	n2 = n
1300	DO i = 1, n2
1301	READ(kin,*) x2(i), y2(i)
1302	END DO
1303	CLOSE (kin)
1304
1305	CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
1306	CALL addpnt(x2,y2,kdata,n2, 0.,0.)
1307	CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
1308	CALL addpnt(x2,y2,kdata,n2, 1.e+38,0.)
1309	CALL inter2(nw,wl,xso2_200,n2,x2,y2,ierr)
1310	IF (ierr .NE. 0) THEN
1311	WRITE(,) ierr, fil
1312	STOP
1313	ENDIF
1314
1315	endif !is_master
1316
1317	call bcast(xso2_200)
1318
1319	fil = 'cross_sections/o2_composite_2018_250K.txt'
1320	print*, 'section efficace O2 250K: ', fil
1321
1322	if(is_master) then
1323
1324	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1325
1326	DO i = 1,nhead
1327	read(kin,*)
1328	END DO
1329
1330	n3 = n
1331	DO i = 1, n3
1332	READ(kin,*) x3(i), y3(i)
1333	END DO
1334	CLOSE (kin)
1335
1336	CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax),0.)
1337	CALL addpnt(x3,y3,kdata,n3, 0.,0.)
1338	CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.)
1339	CALL addpnt(x3,y3,kdata,n3, 1.e+38,0.)
1340	CALL inter2(nw,wl,xso2_250,n3,x3,y3,ierr)
1341	IF (ierr .NE. 0) THEN
1342	WRITE(,) ierr, fil
1343	STOP
1344	ENDIF
1345
1346	endif !is_master
1347
1348	call bcast(xso2_250)
1349
1350	fil = 'cross_sections/o2_composite_2018_300K.txt'
1351	print*, 'section efficace O2 300K: ', fil
1352
1353	if(is_master) then
1354
1355	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1356
1357	DO i = 1,nhead
1358	read(kin,*)
1359	END DO
1360
1361	n4 = n
1362	DO i = 1, n4
1363	READ(kin,*) x4(i), y4(i)
1364	END DO
1365	CLOSE (kin)
1366
1367	CALL addpnt(x4,y4,kdata,n4,x4(1)*(1.-deltax),0.)
1368	CALL addpnt(x4,y4,kdata,n4, 0.,0.)
1369	CALL addpnt(x4,y4,kdata,n4,x4(n4)*(1.+deltax),0.)
1370	CALL addpnt(x4,y4,kdata,n4, 1.e+38,0.)
1371	CALL inter2(nw,wl,xso2_300,n4,x4,y4,ierr)
1372	IF (ierr .NE. 0) THEN
1373	WRITE(,) ierr, fil
1374	STOP
1375	ENDIF
1376
1377	endif !is_master
1378
1379	call bcast(xso2_300)
1380
1381	! photodissociation yield
1382
1383	fil = 'cross_sections/efdis_co2-o2_schunkandnagy2000.txt'
1384
1385	if(is_master) then
1386
1387	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1388
1389	do i = 1,11
1390	read(kin,*)
1391	end do
1392
1393	n = 9
1394	DO i = 1, n
1395	READ(kin,*) xl, xu, dummy, ion(i)
1396	xion(i) = (xl + xu)/2.
1397	ion(i) = max(ion(i), 0.)
1398	END DO
1399	CLOSE (kin)
1400
1401	CALL addpnt(xion,ion,kdata,n,xion(1)*(1.-deltax),0.)
1402	CALL addpnt(xion,ion,kdata,n, 0.,0.)
1403	CALL addpnt(xion,ion,kdata,n,xion(n)*(1.+deltax),1.)
1404	CALL addpnt(xion,ion,kdata,n, 1.e+38,1.)
1405	CALL inter2(nw,wl,yieldo2,n,xion,ion,ierr)
1406	IF (ierr .NE. 0) THEN
1407	WRITE(,) ierr, fil
1408	STOP
1409	ENDIF
1410
1411	endif !is_master
1412
1413	call bcast(yieldo2)
1414
1415	end subroutine rdxso2
1416
1417	!==============================================================================
1418
1419	subroutine rdxso3(nw,wl,xso3_218,xso3_298)
1420
1421	!-----------------------------------------------------------------------------*
1422	!= PURPOSE: =*
1423	!= Read ozone molecular absorption cross section. Re-grid data to match =*
1424	!= specified wavelength working grid. =*
1425	!-----------------------------------------------------------------------------*
1426	!= PARAMETERS: =*
1427	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
1428	!= wavelength grid =*
1429	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
1430	!= working wavelength grid =*
1431	!= XSO3_218 REAL, molecular absoprtion cross section (cm^2) of O3 at (O)=*
1432	!= each specified wavelength (JPL 2006) 218 K =*
1433	!= XSO3_298 REAL, molecular absoprtion cross section (cm^2) of O3 at (O)=*
1434	!= each specified wavelength (JPL 2006) 298 K =*
1435	!-----------------------------------------------------------------------------*
1436
1437	USE mod_phys_lmdz_para, ONLY: is_master
1438	USE mod_phys_lmdz_transfert_para, ONLY: bcast
1439
1440	implicit none
1441
1442	! input
1443
1444	integer :: nw ! number of wavelength grid points
1445	real, dimension(nw) :: wl ! lower and central wavelength for each interval
1446
1447	! output
1448
1449	real, dimension(nw) :: xso3_218, xso3_298 ! o3 cross-sections (cm2)
1450
1451	! local
1452
1453	integer, parameter :: kdata = 200
1454	real, parameter :: deltax = 1.e-4
1455	real, dimension(kdata) :: x1, x2, y1, y2
1456	real, dimension(nw) :: yg
1457	real :: a1, a2
1458
1459	integer :: i, ierr, iw, n, n1, n2
1460	integer :: kin, kout ! input/output logical units
1461
1462	character*100 fil
1463
1464	kin = 10
1465
1466	! JPL 2006 218 K
1467
1468	fil = 'cross_sections/o3_cross-sections_jpl_2006_218K.txt'
1469	print*, 'section efficace O3 218K: ', fil
1470
1471	if(is_master) then
1472
1473	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1474	n1 = 167
1475	DO i = 1, n1
1476	READ(kin,*) a1, a2, y1(i)
1477	x1(i) = (a1+a2)/2.
1478	END DO
1479	CLOSE (kin)
1480
1481	CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
1482	CALL addpnt(x1,y1,kdata,n1, 0.,0.)
1483	CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
1484	CALL addpnt(x1,y1,kdata,n1, 1.e+38,0.)
1485	CALL inter2(nw,wl,yg,n1,x1,y1,ierr)
1486	IF (ierr .NE. 0) THEN
1487	WRITE(,) ierr, fil
1488	STOP
1489	ENDIF
1490
1491	DO iw = 1, nw-1
1492	xso3_218(iw) = yg(iw)
1493	END DO
1494
1495	endif !is_master
1496
1497	call bcast(xso3_218)
1498
1499	! JPL 2006 298 K
1500
1501	fil = 'cross_sections/o3_cross-sections_jpl_2006_298K.txt'
1502	print*, 'section efficace O3 298K: ', fil
1503
1504	if(is_master) then
1505
1506	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1507	n2 = 167
1508	DO i = 1, n2
1509	READ(kin,*) a1, a2, y2(i)
1510	x2(i) = (a1+a2)/2.
1511	END DO
1512	CLOSE (kin)
1513
1514	CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
1515	CALL addpnt(x2,y2,kdata,n2, 0.,0.)
1516	CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
1517	CALL addpnt(x2,y2,kdata,n2, 1.e+38,0.)
1518	CALL inter2(nw,wl,yg,n2,x2,y2,ierr)
1519	IF (ierr .NE. 0) THEN
1520	WRITE(,) ierr, fil
1521	STOP
1522	ENDIF
1523
1524	DO iw = 1, nw-1
1525	xso3_298(iw) = yg(iw)
1526	END DO
1527
1528	endif !is_master
1529
1530	call bcast(xso3_298)
1531
1532	end subroutine rdxso3
1533
1534	!==============================================================================
1535
1536	subroutine rdxss2(nw, wl, yg)
1537
1538	!-----------------------------------------------------------------------------*
1539	!= PURPOSE: =*
1540	!= Read S2 molecular absorption cross section. Re-grid data to match =*
1541	!= specified wavelength working grid. =*
1542	!-----------------------------------------------------------------------------*
1543	!= PARAMETERS: =*
1544	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
1545	!= wavelength grid =*
1546	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
1547	!= working wavelength grid =*
1548	!= YG - REAL, molecular absoprtion cross section (cm^2) of H2O at (O)=*
1549	!= each specified wavelength =*
1550	!-----------------------------------------------------------------------------*
1551
1552	USE mod_phys_lmdz_para, ONLY: is_master
1553	USE mod_phys_lmdz_transfert_para, ONLY: bcast
1554
1555	IMPLICIT NONE
1556
1557	! input
1558
1559	integer :: nw ! number of wavelength grid points
1560	real, dimension(nw) :: wl ! lower and central wavelength for each interval
1561
1562	! output
1563
1564	real, dimension(nw) :: yg ! h2o cross-sections (cm2)
1565
1566	! local
1567
1568	integer, parameter :: kdata = 1500
1569	real, parameter :: deltax = 1.e-4
1570	REAL x1(kdata)
1571	REAL y1(kdata)
1572	INTEGER ierr, dummy
1573	INTEGER i, n
1574
1575	CHARACTER*100 fil
1576	integer :: kin, kout ! input/output logical units
1577
1578	kin = 10
1579
1580	fil = 'cross_sections/s2_millsBestEst_1560_5059.txt'
1581	print*, 'section efficace S2: ', fil
1582
1583	if(is_master) then
1584
1585	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1586
1587	n = 1203
1588	DO i = 1, n
1589	READ(kin,*) x1(i), y1(i)
1590	x1(i) = x1(i)/10
1591	END DO
1592	CLOSE (kin)
1593
1594	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1595	CALL addpnt(x1,y1,kdata,n, 0.,0.)
1596	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1597	CALL addpnt(x1,y1,kdata,n, 1.e+38,0.)
1598	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
1599	IF (ierr .NE. 0) THEN
1600	WRITE(,) ierr, fil
1601	STOP
1602	ENDIF
1603
1604	endif !is_master
1605
1606	call bcast(yg)
1607
1608	end subroutine rdxss2
1609
1610	!==============================================================================
1611
1612	subroutine rdxsh2o(nw, wl, yg)
1613
1614	!-----------------------------------------------------------------------------*
1615	!= PURPOSE: =*
1616	!= Read H2O molecular absorption cross section. Re-grid data to match =*
1617	!= specified wavelength working grid. =*
1618	!-----------------------------------------------------------------------------*
1619	!= PARAMETERS: =*
1620	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
1621	!= wavelength grid =*
1622	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
1623	!= working wavelength grid =*
1624	!= YG - REAL, molecular absoprtion cross section (cm^2) of H2O at (O)=*
1625	!= each specified wavelength =*
1626	!-----------------------------------------------------------------------------*
1627
1628	USE mod_phys_lmdz_para, ONLY: is_master
1629	USE mod_phys_lmdz_transfert_para, ONLY: bcast
1630
1631	IMPLICIT NONE
1632
1633	! input
1634
1635	integer :: nw ! number of wavelength grid points
1636	real, dimension(nw) :: wl ! lower and central wavelength for each interval
1637
1638	! output
1639
1640	real, dimension(nw) :: yg ! h2o cross-sections (cm2)
1641
1642	! local
1643
1644	integer, parameter :: kdata = 500
1645	real, parameter :: deltax = 1.e-4
1646	REAL x1(kdata)
1647	REAL y1(kdata)
1648	INTEGER ierr
1649	INTEGER i, n
1650	CHARACTER*100 fil
1651	integer :: kin, kout ! input/output logical units
1652
1653	kin = 10
1654
1655	fil = 'cross_sections/h2o_composite_250K.txt'
1656	print*, 'section efficace H2O: ', fil
1657
1658	if(is_master) then
1659
1660	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1661
1662	DO i = 1,26
1663	read(kin,*)
1664	END DO
1665
1666	n = 420
1667	DO i = 1, n
1668	READ(kin,*) x1(i), y1(i)
1669	END DO
1670	CLOSE (kin)
1671
1672	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1673	CALL addpnt(x1,y1,kdata,n, 0.,0.)
1674	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1675	CALL addpnt(x1,y1,kdata,n, 1.e+38,0.)
1676	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
1677	IF (ierr .NE. 0) THEN
1678	WRITE(,) ierr, fil
1679	STOP
1680	ENDIF
1681
1682	endif !is_master
1683
1684	call bcast(yg)
1685
1686	end subroutine rdxsh2o
1687
1688	!==============================================================================
1689
1690	subroutine rdxshdo(nw, wl, yg)
1691
1692	!-----------------------------------------------------------------------------*
1693	!= PURPOSE: =*
1694	!= Read HDO molecular absorption cross section. Re-grid data to match =*
1695	!= specified wavelength working grid. =*
1696	!-----------------------------------------------------------------------------*
1697	!= PARAMETERS: =*
1698	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
1699	!= wavelength grid =*
1700	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
1701	!= working wavelength grid =*
1702	!= YG - REAL, molecular absoprtion cross section (cm^2) of HDO at (O)=*
1703	!= each specified wavelength =*
1704	!-----------------------------------------------------------------------------*
1705
1706	USE mod_phys_lmdz_para, ONLY: is_master
1707	USE mod_phys_lmdz_transfert_para, ONLY: bcast
1708
1709	IMPLICIT NONE
1710
1711	! input
1712
1713	integer :: nw ! number of wavelength grid points
1714	real, dimension(nw) :: wl ! lower and central wavelength for each interval
1715
1716	! output
1717
1718	real, dimension(nw) :: yg ! hdo cross-sections (cm2)
1719
1720	! local
1721
1722	integer, parameter :: kdata = 900
1723	real, parameter :: deltax = 1.e-4
1724	REAL x1(kdata)
1725	REAL y1(kdata)
1726	INTEGER ierr
1727	INTEGER i, n
1728	CHARACTER*100 fil
1729	integer :: kin, kout ! input/output logical units
1730
1731	kin = 10
1732
1733	fil = 'cross_sections/hdo_composite_295K.txt'
1734	print*, 'section efficace HDO: ', fil
1735
1736	if(is_master) then
1737
1738	OPEN(UNIT=kin,FILE=fil,STATUS='old')
1739
1740	DO i = 1,17
1741	read(kin,*)
1742	END DO
1743
1744	n = 806
1745	DO i = 1, n
1746	READ(kin,*) x1(i), y1(i)
1747	END DO
1748	CLOSE (kin)
1749
1750	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1751	CALL addpnt(x1,y1,kdata,n, 0.,0.)
1752	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1753	CALL addpnt(x1,y1,kdata,n, 1.e+38,0.)
1754	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
1755	IF (ierr .NE. 0) THEN
1756	WRITE(,) ierr, fil
1757	STOP
1758	ENDIF
1759
1760	endif !is_master
1761
1762	call bcast(yg)
1763
1764	end subroutine rdxshdo
1765
1766	!==============================================================================
1767
1768	subroutine rdxsh2o2(nw, wl, xsh2o2)
1769
1770	!-----------------------------------------------------------------------------*
1771	!= PURPOSE: =*
1772	!= Read and grid H2O2 cross-sections
1773	!= H2O2 + hv -> 2 OH =*
1774	!= Cross section: Schuergers and Welge, Z. Naturforsch. 23a (1968) 1508 =*
1775	!= from 125 to 185 nm, then JPL97 from 190 to 350 nm. =*
1776	!-----------------------------------------------------------------------------*
1777	!= PARAMETERS: =*
1778	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
1779	!= wavelength grid =*
1780	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
1781	!= working wavelength grid =*
1782	!-----------------------------------------------------------------------------*
1783
1784	USE mod_phys_lmdz_para, ONLY: is_master
1785	USE mod_phys_lmdz_transfert_para, ONLY: bcast
1786
1787	implicit none
1788
1789	! input
1790
1791	integer :: nw ! number of wavelength grid points
1792	real, dimension(nw) :: wl ! lower wavelength for each interval
1793
1794	! output
1795
1796	real, dimension(nw) :: xsh2o2 ! h2o2 cross-sections (cm2)
1797
1798	! local
1799
1800	real, parameter :: deltax = 1.e-4
1801	integer, parameter :: kdata = 100
1802	real, dimension(kdata) :: x1, y1
1803	real, dimension(nw) :: yg
1804	integer :: i, ierr, iw, n, idum
1805	integer :: kin, kout ! input/output logical units
1806	character*100 fil
1807
1808	kin = 10
1809
1810	! read cross-sections
1811
1812	fil = 'cross_sections/h2o2_composite.txt'
1813	print*, 'section efficace H2O2: ', fil
1814
1815	if(is_master) then
1816
1817	OPEN(kin,FILE=fil,STATUS='OLD')
1818	READ(kin,*) idum,n
1819	DO i = 1, idum-2
1820	READ(kin,*)
1821	ENDDO
1822	DO i = 1, n
1823	READ(kin,*) x1(i), y1(i)
1824	ENDDO
1825	CLOSE (kin)
1826
1827	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1828	CALL addpnt(x1,y1,kdata,n, 0.,0.)
1829	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1830	CALL addpnt(x1,y1,kdata,n, 1.e+38,0.)
1831	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
1832	IF (ierr .NE. 0) THEN
1833	WRITE(,) ierr, fil
1834	STOP
1835	ENDIF
1836
1837	DO iw = 1, nw - 1
1838	xsh2o2(iw) = yg(iw)
1839	END DO
1840
1841	endif !is_master
1842
1843	call bcast(xsh2o2)
1844
1845	end subroutine rdxsh2o2
1846
1847	!==============================================================================
1848
1849	subroutine rdxsho2(nw, wl, yg)
1850
1851	!-----------------------------------------------------------------------------*
1852	!= PURPOSE: =*
1853	!= Read ho2 cross-sections =*
1854	!= JPL 2006 recommendation =*
1855	!-----------------------------------------------------------------------------*
1856	!= PARAMETERS: =*
1857	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
1858	!= wavelength grid =*
1859	!-----------------------------------------------------------------------------*
1860
1861	USE mod_phys_lmdz_para, ONLY: is_master
1862	USE mod_phys_lmdz_transfert_para, ONLY: bcast
1863
1864	IMPLICIT NONE
1865
1866	! input
1867
1868	integer :: nw ! number of wavelength grid points
1869	real, dimension(nw) :: wl ! lower wavelength for each interval
1870
1871	! output
1872
1873	real, dimension(nw) :: yg ! ho2 cross-sections (cm2)
1874
1875	! local
1876
1877	real, parameter :: deltax = 1.e-4
1878	integer, parameter :: kdata = 100
1879	real, dimension(kdata) :: x1, y1
1880	integer :: i, n, ierr
1881	character*100 fil
1882	integer :: kin, kout ! input/output logical units
1883
1884	kin = 10
1885
1886	!*** cross sections from Sander et al. [2003]
1887
1888	fil = 'cross_sections/ho2_jpl2003.txt'
1889	print*, 'section efficace HO2: ', fil
1890
1891	if(is_master) then
1892
1893	OPEN(kin,FILE=fil,STATUS='OLD')
1894	READ(kin,*) n
1895	DO i = 1, n
1896	READ(kin,*) x1(i), y1(i)
1897	ENDDO
1898	CLOSE(kin)
1899
1900	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1901	CALL addpnt(x1,y1,kdata,n, 0.,0.)
1902	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1903	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
1904
1905	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
1906
1907	IF (ierr .NE. 0) THEN
1908	WRITE(,) ierr, fil
1909	STOP
1910	ENDIF
1911
1912	endif !is_master
1913
1914	call bcast(yg)
1915
1916	end subroutine rdxsho2
1917
1918	!==============================================================================
1919
1920	subroutine rdxshcl(nw, wl, yg)
1921
1922	!-----------------------------------------------------------------------------*
1923	!= PURPOSE: =*
1924	!= Read hcl cross-sections =*
1925	!= JPL 2006 recommendation =*
1926	!-----------------------------------------------------------------------------*
1927	!= PARAMETERS: =*
1928	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
1929	!= wavelength grid =*
1930	!-----------------------------------------------------------------------------*
1931
1932	USE mod_phys_lmdz_para, ONLY: is_master
1933	USE mod_phys_lmdz_transfert_para, ONLY: bcast
1934
1935	IMPLICIT NONE
1936
1937	! input
1938
1939	integer :: nw ! number of wavelength grid points
1940	real, dimension(nw) :: wl ! lower wavelength for each interval
1941
1942	! output
1943
1944	real, dimension(nw) :: yg ! hcl cross-sections (cm2)
1945
1946	! local
1947
1948	real, parameter :: deltax = 1.e-4
1949	integer, parameter :: kdata = 100
1950	real, dimension(kdata) :: x1, y1
1951	integer :: i, n, ierr
1952	character*100 fil
1953	integer :: kin, kout ! input/output logical units
1954
1955	kin = 10
1956
1957	!*** cross sections from JPL [2006]
1958
1959	fil = 'cross_sections/hcl_jpl2006.txt'
1960	print*, 'section efficace HCl: ', fil
1961
1962	if(is_master) then
1963	n = 31
1964	OPEN(kin,FILE=fil,STATUS='OLD')
1965	DO i = 1,4
1966	READ(kin,*)
1967	ENDDO
1968	DO i = 1,n
1969	READ(kin,*) x1(i), y1(i)
1970	ENDDO
1971	CLOSE(kin)
1972
1973	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1974	CALL addpnt(x1,y1,kdata,n, 0.,0.)
1975	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1976	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
1977	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
1978
1979
1980	IF (ierr .NE. 0) THEN
1981	WRITE(,) ierr, fil
1982	STOP
1983	ENDIF
1984
1985	endif !is_master
1986
1987	call bcast(yg)
1988
1989	end subroutine rdxshcl
1990
1991	!==============================================================================
1992
1993	subroutine rdxscl2(nw, wl, yg)
1994
1995	!-----------------------------------------------------------------------------*
1996	!= PURPOSE: =*
1997	!= Read cl2 cross-sections =*
1998	!= JPL 2006 recommendation =*
1999	!-----------------------------------------------------------------------------*
2000	!= PARAMETERS: =*
2001	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
2002	!= wavelength grid =*
2003	!-----------------------------------------------------------------------------*
2004
2005	USE mod_phys_lmdz_para, ONLY: is_master
2006	USE mod_phys_lmdz_transfert_para, ONLY: bcast
2007
2008	IMPLICIT NONE
2009
2010	! input
2011
2012	integer :: nw ! number of wavelength grid points
2013	real, dimension(nw) :: wl ! lower wavelength for each interval
2014
2015	! output
2016
2017	real, dimension(nw) :: yg ! cl2 cross-sections (cm2)
2018
2019	! local
2020
2021	real, parameter :: deltax = 1.e-4
2022	integer, parameter :: kdata = 100
2023	real, dimension(kdata) :: x1, y1
2024	integer :: i, n, ierr
2025	character*100 fil
2026	integer :: kin, kout ! input/output logical units
2027
2028	kin = 10
2029
2030	!*** cross sections from JPL [2006]
2031
2032	fil = 'cross_sections/cl2_jpl2006.txt'
2033	print*, 'section efficace Cl2: ', fil
2034
2035	if(is_master) then
2036
2037	n = 30
2038	OPEN(kin,FILE=fil,STATUS='OLD')
2039	DO i = 1,4
2040	READ(kin,*)
2041	ENDDO
2042	DO i = 1,n
2043	READ(kin,*) x1(i), y1(i)
2044	ENDDO
2045	CLOSE(kin)
2046
2047	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2048	CALL addpnt(x1,y1,kdata,n, 0.,0.)
2049	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2050	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
2051
2052	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2053
2054	IF (ierr .NE. 0) THEN
2055	WRITE(,) ierr, fil
2056	STOP
2057	ENDIF
2058
2059	endif !is_master
2060
2061	call bcast(yg)
2062
2063	end subroutine rdxscl2
2064
2065	!==============================================================================
2066
2067	subroutine rdxshocl(nw, wl, yg)
2068
2069	!-----------------------------------------------------------------------------*
2070	!= PURPOSE: =*
2071	!= Read hocl cross-sections =*
2072	!= JPL 2000 recommendation =*
2073	!-----------------------------------------------------------------------------*
2074	!= PARAMETERS: =*
2075	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
2076	!= wavelength grid =*
2077	!-----------------------------------------------------------------------------*
2078
2079	USE mod_phys_lmdz_para, ONLY: is_master
2080	USE mod_phys_lmdz_transfert_para, ONLY: bcast
2081
2082	IMPLICIT NONE
2083
2084	! input
2085
2086	integer :: nw ! number of wavelength grid points
2087	real, dimension(nw) :: wl ! lower wavelength for each interval
2088
2089	! output
2090
2091	real, dimension(nw) :: yg ! hocl cross-sections (cm2)
2092
2093	! local
2094
2095	real, parameter :: deltax = 1.e-4
2096	integer, parameter :: kdata = 200
2097	real, dimension(kdata) :: x1, y1
2098	integer :: i, n, ierr
2099	character*100 fil
2100	integer :: kin, kout ! input/output logical units
2101
2102	kin = 10
2103
2104	!*** cross sections from JPL [2000]
2105
2106	fil = 'cross_sections/HOCl_jpl2000.txt'
2107	print*, 'section efficace HOCl: ', fil
2108
2109	if(is_master) then
2110
2111	n = 111
2112	OPEN(kin,FILE=fil,STATUS='OLD')
2113	DO i = 1,6
2114	READ(kin,*)
2115	ENDDO
2116	DO i = 1,111
2117	READ(kin,*) y1(i)
2118	x1(i) = 200 + real(i-1)*2.
2119	y1(i) = y1(i) * 1.E-20
2120	ENDDO
2121	CLOSE(kin)
2122
2123	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2124	CALL addpnt(x1,y1,kdata,n, 0.,0.)
2125	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2126	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
2127
2128	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2129
2130	IF (ierr .NE. 0) THEN
2131	WRITE(,) ierr, fil
2132	STOP
2133	ENDIF
2134
2135	endif !is_master
2136
2137	call bcast(yg)
2138
2139	end subroutine rdxshocl
2140
2141	!==============================================================================
2142
2143	subroutine rdxsso2(nw,wl,xsso2_200,xsso2_298,xsso2_360)
2144
2145	!-----------------------------------------------------------------------------*
2146	!= PURPOSE: =*
2147	!= Read and grid SO2 absorption cross-sections and photodissociation yield =*
2148	!-----------------------------------------------------------------------------*
2149	!= PARAMETERS: =*
2150	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
2151	!= wavelength grid =*
2152	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
2153	!= working wavelength grid =*
2154	!= XSSO2 - REAL, molecular absoprtion cross section (cm^2) of SO2 at (O)=*
2155	!= each specified wavelength =*
2156	!-----------------------------------------------------------------------------*
2157
2158	USE mod_phys_lmdz_para, ONLY: is_master
2159	USE mod_phys_lmdz_transfert_para, ONLY: bcast
2160
2161	implicit none
2162
2163	! input
2164
2165	integer :: nw ! number of wavelength grid points
2166	real, dimension(nw) :: wl ! lower and central wavelength for each interval
2167
2168	! output
2169
2170	real, dimension(nw) :: xsso2_200, xsso2_298, xsso2_360 ! so2 cross-sections (cm2)
2171
2172	! local
2173
2174	integer, parameter :: kdata = 55000
2175	real, parameter :: deltax = 1.e-4
2176	real, dimension(kdata) :: x1, y1, y2, y3, xion, ion
2177	real, dimension(nw) :: yg
2178	real :: xl, xu
2179	integer :: ierr, i, l, n, n1, n2, n3, n4
2180	CHARACTER*100 fil
2181
2182	integer :: kin, kout ! input/ouput logical units
2183
2184	kin = 10
2185	kout = 30
2186
2187	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2188	!
2189	! SO2 absorption cross-sections
2190	!
2191	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2192	!
2193	! iopt = 1
2194	!
2195	! 200K: Wu et al. (2000) + Vandaele et al. (2009) + Hermans et al. (2009)
2196	! 7 lignes d'en-tete et n1 = 50276
2197	! fichier: so2_composite_200K.txt
2198	!
2199	! 298K: Wu et al. (2000) + Vandaele et al. (2009) + Hermans et al. (2009)
2200	! 7 lignes d'en-tete et n2 = 49833
2201	! fichier: so2_composite_298K.txt
2202	!
2203	! 360K: Wu et al. (2000) + Vandaele et al. (2009) + Hermans et al. (2009)
2204	! 7 lignes d'en-tete et n3 = 49261
2205	! fichier: so2_composite_360K.txt
2206	!
2207	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2208
2209
2210	n1 = 50276
2211	n2 = 49833
2212	n3 = 46261
2213
2214	! 200K:
2215
2216	fil = 'cross_sections/so2_composite_200K.txt'
2217	print*, 'section efficace SO2 195K: ', fil
2218
2219	if(is_master) then
2220
2221	OPEN(UNIT=kin,FILE=fil,STATUS='old')
2222	DO i = 1,7
2223	read(kin,*)
2224	END DO
2225
2226	DO i = 1, n1
2227	READ(kin,*) x1(i), y1(i)
2228	END DO
2229	CLOSE (kin)
2230
2231	CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
2232	CALL addpnt(x1,y1,kdata,n1, 0.,0.)
2233	CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
2234	CALL addpnt(x1,y1,kdata,n1, 1.e+38,0.)
2235	CALL inter2(nw,wl,yg,n1,x1,y1,ierr)
2236	IF (ierr .NE. 0) THEN
2237	WRITE(,) ierr, fil
2238	STOP
2239	ENDIF
2240
2241	DO l = 1, nw-1
2242	xsso2_200(l) = yg(l)
2243	END DO
2244
2245	endif !is_master
2246
2247	call bcast(xsso2_200)
2248
2249	! 298K:
2250
2251	fil = 'cross_sections/so2_composite_298K.txt'
2252	print*, 'section efficace SO2 295K: ', fil
2253
2254	if(is_master) then
2255
2256	OPEN(UNIT=kin,FILE=fil,STATUS='old')
2257	DO i = 1,7
2258	read(kin,*)
2259	END DO
2260
2261	DO i = 1, n2
2262	READ(kin,*) x1(i), y1(i)
2263	END DO
2264	CLOSE (kin)
2265
2266	CALL addpnt(x1,y1,kdata,n2,x1(1)*(1.-deltax),0.)
2267	CALL addpnt(x1,y1,kdata,n2, 0.,0.)
2268	CALL addpnt(x1,y1,kdata,n2,x1(n2)*(1.+deltax),0.)
2269	CALL addpnt(x1,y1,kdata,n2, 1.e+38,0.)
2270	CALL inter2(nw,wl,yg,n2,x1,y1,ierr)
2271	IF (ierr .NE. 0) THEN
2272	WRITE(,) ierr, fil
2273	STOP
2274	ENDIF
2275
2276	DO l = 1, nw-1
2277	xsso2_298(l) = yg(l)
2278	END DO
2279
2280	endif !is_master
2281
2282	call bcast(xsso2_298)
2283
2284	! 360K:
2285
2286	fil = 'cross_sections/so2_composite_360K.txt'
2287	print*, 'section efficace SO2 370K: ', fil
2288
2289	if(is_master) then
2290
2291	OPEN(UNIT=kin,FILE=fil,STATUS='old')
2292	DO i = 1,7
2293	read(kin,*)
2294	END DO
2295
2296	DO i = 1, n3
2297	READ(kin,*) x1(i), y1(i)
2298	END DO
2299	CLOSE (kin)
2300
2301	CALL addpnt(x1,y1,kdata,n3,x1(1)*(1.-deltax),0.)
2302	CALL addpnt(x1,y1,kdata,n3, 0.,0.)
2303	CALL addpnt(x1,y1,kdata,n3,x1(n3)*(1.+deltax),0.)
2304	CALL addpnt(x1,y1,kdata,n3, 1.e+38,0.)
2305	CALL inter2(nw,wl,yg,n3,x1,y1,ierr)
2306	IF (ierr .NE. 0) THEN
2307	WRITE(,) ierr, fil
2308	STOP
2309	ENDIF
2310
2311	DO l = 1, nw-1
2312	xsso2_360(l) = yg(l)
2313	END DO
2314
2315	endif !is_master
2316
2317	call bcast(xsso2_360)
2318
2319	! DO l = 1, nw-1
2320	! write(kout,*) wl(l), xsso2_200(l),
2321	! $ xsso2_298(l),
2322	! $ xsso2_360(l),
2323	! END DO
2324
2325	end subroutine rdxsso2
2326
2327	!==============================================================================
2328
2329	subroutine rdxsso(nw, wl, yg)
2330
2331	!-----------------------------------------------------------------------------*
2332	!= PURPOSE: =*
2333	!= Read SO cross-sections =*
2334	!= JPL 2006 recommendation =*
2335	!-----------------------------------------------------------------------------*
2336	!= PARAMETERS: =*
2337	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
2338	!= wavelength grid =*
2339	!-----------------------------------------------------------------------------*
2340
2341	USE mod_phys_lmdz_para, ONLY: is_master
2342	USE mod_phys_lmdz_transfert_para, ONLY: bcast
2343
2344	IMPLICIT NONE
2345
2346	! input
2347	integer :: nw ! number of wavelength grid points
2348	real, dimension(nw) :: wl ! lower wavelength for each interval
2349
2350	! output
2351
2352	real, dimension(nw) :: yg ! SO cross-sections (cm2)
2353
2354	! local
2355	real :: dummy
2356	real, parameter :: deltax = 1.e-4
2357	integer, parameter :: kdata = 1000
2358	real, dimension(kdata) :: x1, y1
2359	integer :: i, n, ierr
2360	character*100 fil
2361	integer :: kin, kout ! input/output logical units
2362
2363	kin = 10
2364
2365	!*** cross sections from JPL [2006]
2366
2367	fil = 'cross_sections/so_marcq.txt'
2368	print*, 'section efficace SO: ', fil
2369
2370	if(is_master) then
2371
2372	n = 851
2373	OPEN(kin,FILE=fil,STATUS='OLD')
2374	DO i = 1,n
2375	READ(kin,*) x1(i), dummy, dummy, dummy, dummy, y1(i)
2376	y1(i) = y1(i) * 0.88
2377	ENDDO
2378	CLOSE(kin)
2379
2380	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2381	CALL addpnt(x1,y1,kdata,n, 0.,0.)
2382	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2383	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
2384
2385	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2386
2387	IF (ierr .NE. 0) THEN
2388	WRITE(,) ierr, fil
2389	STOP
2390	ENDIF
2391
2392	endif !is_master
2393
2394	call bcast(yg)
2395
2396	end subroutine rdxsso
2397
2398	!==============================================================================
2399
2400	subroutine rdxsso3(nw, wl, yg)
2401
2402	!-----------------------------------------------------------------------------*
2403	!= PURPOSE: =*
2404	!= Read SO3 cross-sections =*
2405	!= Composite section 140-294 nm -> Pintze al. [2003] =*
2406	!= 296-330 nm -> Burkholder al. [1997] =*
2407	!-----------------------------------------------------------------------------*
2408	!= PARAMETERS: =*
2409	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
2410	!= wavelength grid =*
2411	!-----------------------------------------------------------------------------*
2412
2413	USE mod_phys_lmdz_para, ONLY: is_master
2414	USE mod_phys_lmdz_transfert_para, ONLY: bcast
2415
2416	IMPLICIT NONE
2417
2418	! input
2419
2420	integer :: nw ! number of wavelength grid points
2421	real, dimension(nw) :: wl ! lower wavelength for each interval
2422
2423	! output
2424
2425	real, dimension(nw) :: yg ! SO3 cross-sections (cm2)
2426
2427	! local
2428
2429	real, parameter :: deltax = 1.e-4
2430	integer, parameter :: kdata = 200
2431	real, dimension(kdata) :: x1, y1
2432	integer :: i, n, ierr
2433	character*100 fil
2434	integer :: kin, kout ! input/output logical units
2435
2436	kin = 10
2437
2438	!*** cross sections
2439
2440	fil = 'cross_sections/so3_composite.txt'
2441	print*, 'section efficace SO3: ', fil
2442
2443	if(is_master) then
2444
2445	n = 96
2446	OPEN(kin,FILE=fil,STATUS='OLD')
2447	DO i = 1,16
2448	READ(kin,*)
2449	ENDDO
2450	DO i = 1,n
2451	READ(kin,*) x1(i), y1(i)
2452	ENDDO
2453	CLOSE(kin)
2454
2455	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2456	CALL addpnt(x1,y1,kdata,n, 0.,0.)
2457	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2458	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
2459
2460	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2461
2462	IF (ierr .NE. 0) THEN
2463	WRITE(,) ierr, fil
2464	STOP
2465	ENDIF
2466
2467	endif !is_master
2468
2469	call bcast(yg)
2470
2471	end subroutine rdxsso3
2472
2473	!==============================================================================
2474
2475	subroutine rdxsclo(nw, wl, yg)
2476
2477	!-----------------------------------------------------------------------------*
2478	!= PURPOSE: =*
2479	!= Read ClO cross-sections =*
2480	!= From Trollier and al. [1990] =*
2481	!-----------------------------------------------------------------------------*
2482	!= PARAMETERS: =*
2483	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
2484	!= wavelength grid =*
2485	!-----------------------------------------------------------------------------*
2486
2487	USE mod_phys_lmdz_para, ONLY: is_master
2488	USE mod_phys_lmdz_transfert_para, ONLY: bcast
2489
2490	IMPLICIT NONE
2491
2492	! input
2493
2494	integer :: nw ! number of wavelength grid points
2495	real, dimension(nw) :: wl ! lower wavelength for each interval
2496
2497	! output
2498
2499	real, dimension(nw) :: yg ! hcl cross-sections (cm2)
2500
2501	! local
2502
2503	real, parameter :: deltax = 1.e-4
2504	integer, parameter :: kdata = 100
2505	real, dimension(kdata) :: x1, y1
2506	integer :: i, n, ierr
2507	character*100 fil
2508	integer :: kin, kout ! input/output logical units
2509
2510	kin = 10
2511
2512	!*** cross sections
2513
2514	fil = 'cross_sections/clo_xs_trolier_1990.txt'
2515	print*, 'section efficace ClO: ', fil
2516
2517	if(is_master) then
2518
2519	n = 81
2520	OPEN(kin,FILE=fil,STATUS='OLD')
2521	DO i = 1,6
2522	READ(kin,*)
2523	ENDDO
2524	DO i = 1,n
2525	READ(kin,*) x1(i), y1(i)
2526	ENDDO
2527	CLOSE(kin)
2528
2529	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2530	CALL addpnt(x1,y1,kdata,n, 0.,0.)
2531	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2532	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
2533
2534	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2535
2536	IF (ierr .NE. 0) THEN
2537	WRITE(,) ierr, fil
2538	STOP
2539	ENDIF
2540
2541	endif !is_master
2542
2543	call bcast(yg)
2544
2545	end subroutine rdxsclo
2546
2547	!==============================================================================
2548
2549	subroutine rdxsocs(nw, wl, yg)
2550
2551	!-----------------------------------------------------------------------------*
2552	!= PURPOSE: =*
2553	!= Read OCS cross-sections =*
2554	!= JPL 2011 recommendation =*
2555	!-----------------------------------------------------------------------------*
2556	!= PARAMETERS: =*
2557	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
2558	!= wavelength grid =*
2559	!-----------------------------------------------------------------------------*
2560
2561	USE mod_phys_lmdz_para, ONLY: is_master
2562	USE mod_phys_lmdz_transfert_para, ONLY: bcast
2563
2564	IMPLICIT NONE
2565
2566	! input
2567
2568	integer :: nw ! number of wavelength grid points
2569	real, dimension(nw) :: wl ! lower wavelength for each interval
2570
2571	! output
2572
2573	real, dimension(nw) :: yg ! hcl cross-sections (cm2)
2574
2575	! local
2576
2577	real, parameter :: deltax = 1.e-4
2578	integer, parameter :: kdata = 100
2579	real, dimension(kdata) :: x1, y1
2580	integer :: i, n, ierr
2581	character*100 fil
2582	integer :: kin, kout ! input/output logical units
2583
2584	kin = 10
2585
2586	!*** cross sections from JPL [2006]
2587
2588	fil = 'cross_sections/ocs_cross_sections_jpl2011.txt'
2589	print*, 'section efficace OCS: ', fil
2590
2591	if(is_master) then
2592
2593	n = 40
2594	OPEN(kin,FILE=fil,STATUS='OLD')
2595	DO i = 1,4
2596	READ(kin,*)
2597	ENDDO
2598	DO i = 1,n
2599	READ(kin,*) x1(i), y1(i)
2600	ENDDO
2601	CLOSE(kin)
2602
2603	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2604	CALL addpnt(x1,y1,kdata,n, 0.,0.)
2605	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2606	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
2607
2608	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2609
2610	IF (ierr .NE. 0) THEN
2611	WRITE(,) ierr, fil
2612	STOP
2613	ENDIF
2614
2615	endif !is_master
2616
2617	call bcast(yg)
2618
2619	end subroutine rdxsocs
2620
2621	!==============================================================================
2622
2623	subroutine rdxscocl2(nw, wl, yg)
2624
2625	!-----------------------------------------------------------------------------*
2626	!= PURPOSE: =*
2627	!= Read COCl2 cross-sections =*
2628	!= JPL 2011 recommendation =*
2629	!-----------------------------------------------------------------------------*
2630	!= PARAMETERS: =*
2631	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
2632	!= wavelength grid =*
2633	!-----------------------------------------------------------------------------*
2634
2635	USE mod_phys_lmdz_para, ONLY: is_master
2636	USE mod_phys_lmdz_transfert_para, ONLY: bcast
2637
2638	IMPLICIT NONE
2639
2640	! input
2641
2642	integer :: nw ! number of wavelength grid points
2643	real, dimension(nw) :: wl ! lower wavelength for each interval
2644
2645	! output
2646
2647	real, dimension(nw) :: yg ! COCl2 cross-sections (cm2)
2648
2649	! local
2650
2651	real, parameter :: deltax = 1.e-4
2652	integer, parameter :: kdata = 100
2653	real, dimension(kdata) :: x1, y1
2654	integer :: i, n, ierr
2655	character*100 fil
2656	integer :: kin, kout ! input/output logical units
2657
2658	kin = 10
2659
2660	!*** cross sections from JPL [2011]
2661
2662	fil = 'cross_sections/cocl2_cross_sections_jpl2011.txt'
2663	print*, 'section efficace COCl2: ', fil
2664
2665	if(is_master) then
2666
2667	n = 53
2668	OPEN(kin,FILE=fil,STATUS='OLD')
2669	DO i = 1,4
2670	READ(kin,*)
2671	ENDDO
2672	DO i = 1,n
2673	READ(kin,*) x1(i), y1(i)
2674	ENDDO
2675	CLOSE(kin)
2676
2677	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2678	CALL addpnt(x1,y1,kdata,n, 0.,0.)
2679	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2680	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
2681
2682	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2683
2684	IF (ierr .NE. 0) THEN
2685	WRITE(,) ierr, fil
2686	STOP
2687	ENDIF
2688
2689	endif !is_master
2690
2691	call bcast(yg)
2692
2693	end subroutine rdxscocl2
2694
2695	!==============================================================================
2696
2697	subroutine rdxsh2so4(nw, wl, yg)
2698
2699	!-----------------------------------------------------------------------------*
2700	!= PURPOSE: =*
2701	!= Read H2SO4 cross-sections =*
2702	!= JPL 2006 recommendation =*
2703	!-----------------------------------------------------------------------------*
2704	!= PARAMETERS: =*
2705	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
2706	!= wavelength grid =*
2707	!-----------------------------------------------------------------------------*
2708
2709	USE mod_phys_lmdz_para, ONLY: is_master
2710	USE mod_phys_lmdz_transfert_para, ONLY: bcast
2711
2712	IMPLICIT NONE
2713
2714	! input
2715
2716	integer :: nw ! number of wavelength grid points
2717	real, dimension(nw) :: wl ! lower wavelength for each interval
2718
2719	! output
2720
2721	real, dimension(nw) :: yg ! h2so4 cross-sections (cm2)
2722
2723	! local
2724
2725	real, parameter :: deltax = 1.e-4
2726	integer, parameter :: kdata = 100
2727	real, dimension(kdata) :: x1, y1
2728	integer :: i, n, ierr
2729	character*100 fil
2730	integer :: kin, kout ! input/output logical units
2731
2732	kin = 10
2733
2734	!*** cross sections from JPL [2006]
2735
2736	fil = 'cross_sections/h2so4_cross_sections.txt'
2737	print*, 'section efficace h2so4: ', fil
2738
2739	if(is_master) then
2740
2741	n = 22
2742	OPEN(kin,FILE=fil,STATUS='OLD')
2743	DO i = 1,3
2744	READ(kin,*)
2745	ENDDO
2746	DO i = 1,n
2747	READ(kin,*) x1(i), y1(i)
2748	ENDDO
2749	CLOSE(kin)
2750
2751	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2752	CALL addpnt(x1,y1,kdata,n, 0.,0.)
2753	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2754	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
2755
2756	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2757
2758	IF (ierr .NE. 0) THEN
2759	WRITE(,) ierr, fil
2760	STOP
2761	ENDIF
2762
2763	endif !is_master
2764
2765	call bcast(yg)
2766
2767	end subroutine rdxsh2so4
2768
2769	!==============================================================================
2770
2771	subroutine rdxsh2(nw, wl, wc, yg, yieldh2)
2772
2773	!-----------------------------------------------------------------------------*
2774	!= PURPOSE: =*
2775	!= Read h2 cross-sections and photodissociation yield =*
2776	!-----------------------------------------------------------------------------*
2777	!= PARAMETERS: =*
2778	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
2779	!= wavelength grid =*
2780	!-----------------------------------------------------------------------------*
2781
2782	USE mod_phys_lmdz_para, ONLY: is_master
2783	USE mod_phys_lmdz_transfert_para, ONLY: bcast
2784
2785	implicit none
2786
2787	! input
2788
2789	integer :: nw ! number of wavelength grid points
2790	real, dimension(nw) :: wl, wc ! lower and central wavelength for each interval
2791
2792	! output
2793
2794	real, dimension(nw) :: yg ! h2 cross-sections (cm2)
2795	real, dimension(nw) :: yieldh2 ! photodissociation yield
2796
2797	! local
2798
2799	integer, parameter :: kdata = 1000
2800	real, parameter :: deltax = 1.e-4
2801	real, dimension(kdata) :: x1, y1, x2, y2
2802	real :: xl, xu
2803	integer :: i, iw, n, ierr
2804	integer :: kin, kout ! input/output logical units
2805	character*100 fil
2806
2807	kin = 10
2808
2809	! h2 cross sections
2810
2811	fil = 'cross_sections/h2secef.txt'
2812	print*, 'section efficace H2: ', fil
2813
2814	if(is_master) then
2815
2816	OPEN(kin,FILE=fil,STATUS='OLD')
2817
2818	n = 792
2819	read(kin,*) ! avoid first line with wavelength = 0.
2820	DO i = 1, n
2821	READ(kin,*) x1(i), y1(i)
2822	ENDDO
2823	CLOSE(kin)
2824
2825	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2826	CALL addpnt(x1,y1,kdata,n, 0.,0.)
2827	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2828	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
2829	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2830
2831	IF (ierr .NE. 0) THEN
2832	WRITE(,) ierr, fil
2833	STOP
2834	ENDIF
2835
2836	endif !is_master
2837
2838	call bcast(yg)
2839
2840	! photodissociation yield
2841
2842	fil = 'cross_sections/h2_ionef_schunknagy2000.txt'
2843
2844	if(is_master) then
2845
2846	OPEN(UNIT=kin,FILE=fil,STATUS='old')
2847
2848	n = 19
2849	read(kin,*)
2850	DO i = 1, n
2851	READ(kin,*) xl, xu, y2(i)
2852	x2(i) = (xl + xu)/2.
2853	y2(i) = max(1. - y2(i),0.)
2854	END DO
2855	CLOSE (kin)
2856
2857	CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
2858	CALL addpnt(x2,y2,kdata,n, 0.,0.)
2859	CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
2860	CALL addpnt(x2,y2,kdata,n, 1.e+38,1.)
2861	CALL inter2(nw,wl,yieldh2,n,x2,y2,ierr)
2862	IF (ierr .NE. 0) THEN
2863	WRITE(,) ierr, fil
2864	STOP
2865	ENDIF
2866
2867	endif !is_master
2868
2869	call bcast(yieldh2)
2870
2871	end subroutine rdxsh2
2872
2873	!==============================================================================
2874
2875	subroutine rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294,yldno2_248, yldno2_298)
2876
2877	!-----------------------------------------------------------------------------*
2878	!= PURPOSE: =*
2879	!= read and grid cross section + quantum yield for NO2 =*
2880	!= photolysis =*
2881	!= Jenouvrier et al., 1996 200-238 nm
2882	!= Vandaele et al., 1998 238-666 nm 220K and 294K
2883	!= quantum yield from jpl 2006
2884	!-----------------------------------------------------------------------------*
2885	!= PARAMETERS: =*
2886	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
2887	!= wavelength grid =*
2888	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
2889	!= working wavelength grid =*
2890	!= SQ - REAL, cross section x quantum yield (cm^2) for each (O)=*
2891	!= photolysis reaction defined, at each defined wavelength and =*
2892	!= at each defined altitude level =*
2893	!-----------------------------------------------------------------------------*
2894
2895	USE mod_phys_lmdz_para, ONLY: is_master
2896	USE mod_phys_lmdz_transfert_para, ONLY: bcast
2897
2898	implicit none
2899
2900	! input
2901
2902	integer :: nw ! number of wavelength grid points
2903	real, dimension(nw) :: wl ! lower and central wavelength for each interval
2904
2905	! output
2906
2907	real, dimension(nw) :: xsno2, xsno2_220, xsno2_294 ! no2 cross-sections (cm2)
2908	real, dimension(nw) :: yldno2_248, yldno2_298 ! quantum yields at 248-298 k
2909
2910	! local
2911
2912	integer, parameter :: kdata = 28000
2913	real, parameter :: deltax = 1.e-4
2914	real, dimension(kdata) :: x1, x2, x3, x4, x5, y1, y2, y3, y4, y5
2915	real, dimension(nw) :: yg1, yg2, yg3, yg4, yg5
2916	real :: dum, qy
2917	integer :: i, iw, n, n1, n2, n3, n4, n5, ierr
2918	character*100 fil
2919	integer :: kin, kout ! input/output logical units
2920
2921	kin = 10
2922
2923	!*************** NO2 photodissociation
2924
2925	! Jenouvrier 1996 + Vandaele 1998 (JPL 2006)
2926
2927	fil = 'cross_sections/no2_xs_jenouvrier.txt'
2928	print*, 'section efficace NO2: ', fil
2929
2930	if (is_master) then
2931
2932	OPEN(UNIT=kin,FILE=fil,status='old')
2933	DO i = 1, 3
2934	READ(kin,*)
2935	END DO
2936	n1 = 10001
2937	DO i = 1, n1
2938	READ(kin,*) x1(i), y1(i)
2939	end do
2940
2941	CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax), 0.)
2942	CALL addpnt(x1,y1,kdata,n1, 0., 0.)
2943	CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
2944	CALL addpnt(x1,y1,kdata,n1, 1.e+38, 0.)
2945	CALL inter2(nw,wl,yg1,n1,x1,y1,ierr)
2946
2947	end if !is_master
2948
2949	call bcast(yg1)
2950
2951	fil = 'cross_sections/no2_xs_vandaele_294K.txt'
2952	print*, 'section efficace NO2: ', fil
2953
2954	if (is_master) then
2955
2956	OPEN(UNIT=kin,FILE=fil,status='old')
2957	DO i = 1, 3
2958	READ(kin,*)
2959	END DO
2960	n2 = 27993
2961	DO i = 1, n2
2962	READ(kin,*) x2(i), y2(i)
2963	end do
2964
2965	CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax), 0.)
2966	CALL addpnt(x2,y2,kdata,n2, 0., 0.)
2967	CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
2968	CALL addpnt(x2,y2,kdata,n2, 1.e+38, 0.)
2969	CALL inter2(nw,wl,yg2,n2,x2,y2,ierr)
2970
2971	end if !is_master
2972
2973	call bcast(yg2)
2974
2975	fil = 'cross_sections/no2_xs_vandaele_220K.txt'
2976	print*, 'section efficace NO2: ', fil
2977
2978	if (is_master) then
2979
2980	OPEN(UNIT=kin,FILE=fil,status='old')
2981	DO i = 1, 3
2982	READ(kin,*)
2983	END DO
2984	n3 = 27993
2985	do i = 1, n3
2986	READ(kin,*) x3(i), y3(i)
2987	end do
2988
2989	CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax), 0.)
2990	CALL addpnt(x3,y3,kdata,n3, 0., 0.)
2991	CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.)
2992	CALL addpnt(x3,y3,kdata,n3, 1.e+38, 0.)
2993	CALL inter2(nw,wl,yg3,n3,x3,y3,ierr)
2994
2995	do iw = 1, nw - 1
2996	xsno2(iw) = yg1(iw)
2997	xsno2_294(iw) = yg2(iw)
2998	xsno2_220(iw) = yg3(iw)
2999	end do
3000
3001	end if !is_master
3002
3003	call bcast(yg3)
3004	call bcast(xsno2)
3005	call bcast(xsno2_294)
3006	call bcast(xsno2_220)
3007
3008	! photodissociation efficiency from jpl 2006
3009
3010	fil = 'cross_sections/no2_yield_jpl2006.txt'
3011	print*, 'quantum yield NO2: ', fil
3012
3013	if (is_master) then
3014
3015	OPEN(UNIT=kin,FILE=fil,STATUS='old')
3016	DO i = 1, 5
3017	READ(kin,*)
3018	END DO
3019	n = 25
3020	n4 = n
3021	n5 = n
3022	DO i = 1, n
3023	READ(kin,*) x4(i), y4(i), y5(i)
3024	x5(i) = x4(i)
3025	END DO
3026	CLOSE(kin)
3027
3028	CALL addpnt(x4,y4,kdata,n4,x4(1)*(1.-deltax),y4(1))
3029	CALL addpnt(x4,y4,kdata,n4, 0.,y4(1))
3030	CALL addpnt(x4,y4,kdata,n4,x4(n4)*(1.+deltax), 0.)
3031	CALL addpnt(x4,y4,kdata,n4, 1.e+38, 0.)
3032	CALL inter2(nw,wl,yg4,n4,x4,y4,ierr)
3033	IF (ierr .NE. 0) THEN
3034	WRITE(,) ierr, fil
3035	STOP
3036	END IF
3037
3038	end if !is_master
3039
3040	call bcast(yg4)
3041
3042	if (is_master) then
3043
3044	CALL addpnt(x5,y5,kdata,n5,x5(1)*(1.-deltax),y5(1))
3045	CALL addpnt(x5,y5,kdata,n5, 0.,y5(1))
3046	CALL addpnt(x5,y5,kdata,n5,x5(n5)*(1.+deltax), 0.)
3047	CALL addpnt(x5,y5,kdata,n5, 1.e+38, 0.)
3048	CALL inter2(nw,wl,yg5,n5,x5,y5,ierr)
3049	IF (ierr .NE. 0) THEN
3050	WRITE(,) ierr, fil
3051	STOP
3052	END IF
3053
3054	do iw = 1, nw - 1
3055	yldno2_298(iw) = yg4(iw)
3056	yldno2_248(iw) = yg5(iw)
3057	end do
3058
3059	end if !is_master
3060
3061	call bcast(yg5)
3062	call bcast(yldno2_298)
3063	call bcast(yldno2_248)
3064
3065	end subroutine rdxsno2
3066
3067	!==============================================================================
3068
3069	subroutine rdxsno(nw, wl, yg, yieldno)
3070
3071	!-----------------------------------------------------------------------------*
3072	!= PURPOSE: =*
3073	!= Read NO cross-sections and photodissociation efficiency =*
3074	!= Lida et al 1986 (provided by Francisco Gonzalez-Galindo) =*
3075	!-----------------------------------------------------------------------------*
3076	!= PARAMETERS: =*
3077	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
3078	!= wavelength grid =*
3079	!-----------------------------------------------------------------------------*
3080
3081	USE mod_phys_lmdz_para, ONLY: is_master
3082	USE mod_phys_lmdz_transfert_para, ONLY: bcast
3083
3084	implicit none
3085
3086	! input
3087
3088	integer :: nw ! number of wavelength grid points
3089	real, dimension(nw) :: wl ! lower wavelength for each interval
3090
3091	! output
3092
3093	real, dimension(nw) :: yg ! no cross-sections (cm2)
3094	real, dimension(nw) :: yieldno ! no photodissociation efficiency
3095
3096	! local
3097
3098	integer, parameter :: kdata = 110
3099	real, parameter :: deltax = 1.e-4
3100	real, dimension(kdata) :: x1, y1, x2, y2
3101	integer :: i, iw, n, ierr
3102	character*100 fil
3103	integer :: kin, kout ! input/output logical units
3104
3105	kin = 10
3106
3107	! no cross-sections
3108
3109	fil = 'cross_sections/no_xs_francisco.txt'
3110	print*, 'section efficace NO: ', fil
3111
3112	if (is_master) then
3113
3114	OPEN(kin,FILE=fil,STATUS='OLD')
3115
3116	n = 99
3117	DO i = 1, n
3118	READ(kin,*) x1(i), y1(i)
3119	END DO
3120	CLOSE(kin)
3121
3122	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
3123	CALL addpnt(x1,y1,kdata,n, 0.,0.)
3124	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
3125	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
3126
3127	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
3128	IF (ierr .NE. 0) THEN
3129	WRITE(,) ierr, fil
3130	STOP
3131	END IF
3132
3133	end if !is_master
3134
3135	call bcast(yg)
3136
3137	! photodissociation yield
3138
3139	fil = 'cross_sections/noefdis.txt'
3140
3141	if (is_master) then
3142
3143	OPEN(UNIT=kin,FILE=fil,STATUS='old')
3144
3145	n = 33
3146	DO i = 1, n
3147	READ(kin,*) x2(n-i+1), y2(n-i+1)
3148	END DO
3149	CLOSE (kin)
3150
3151	CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
3152	CALL addpnt(x2,y2,kdata,n, 0.,0.)
3153	CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
3154	CALL addpnt(x2,y2,kdata,n, 1.e+38,1.)
3155	CALL inter2(nw,wl,yieldno,n,x2,y2,ierr)
3156	IF (ierr .NE. 0) THEN
3157	WRITE(,) ierr, fil
3158	STOP
3159	END IF
3160
3161	end if !is_master
3162
3163	call bcast(yieldno)
3164
3165	end subroutine rdxsno
3166
3167	!==============================================================================
3168
3169	subroutine rdxsn2(nw, wl, yg, yieldn2)
3170
3171	!-----------------------------------------------------------------------------*
3172	!= PURPOSE: =*
3173	!= Read n2 cross-sections and photodissociation yield =*
3174	!-----------------------------------------------------------------------------*
3175	!= PARAMETERS: =*
3176	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
3177	!= wavelength grid =*
3178	!-----------------------------------------------------------------------------*
3179
3180	USE mod_phys_lmdz_para, ONLY: is_master
3181	USE mod_phys_lmdz_transfert_para, ONLY: bcast
3182
3183	implicit none
3184
3185	! input
3186
3187	integer :: nw ! number of wavelength grid points
3188	real, dimension(nw) :: wl ! lower wavelength for each interval
3189
3190	! output
3191
3192	real, dimension(nw) :: yg ! n2 cross-sections (cm2)
3193	real, dimension(nw) :: yieldn2 ! n2 photodissociation yield
3194
3195	! local
3196
3197	integer, parameter :: kdata = 1100
3198	real, parameter :: deltax = 1.e-4
3199	real, dimension(kdata) :: x1, y1, x2, y2
3200	real :: xl, xu
3201	integer :: i, iw, n, ierr
3202	integer :: kin, kout ! input/output logical units
3203	character*100 fil
3204
3205	kin = 10
3206
3207	! n2 cross sections
3208
3209	fil = 'cross_sections/n2secef_01nm.txt'
3210	print*, 'section efficace N2: ', fil
3211
3212	if (is_master) then
3213
3214	OPEN(kin,FILE=fil,STATUS='OLD')
3215
3216	n = 1020
3217	DO i = 1, n
3218	READ(kin,*) x1(i), y1(i)
3219	END DO
3220	CLOSE(kin)
3221
3222	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
3223	CALL addpnt(x1,y1,kdata,n, 0.,0.)
3224	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
3225	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
3226
3227	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
3228
3229	IF (ierr .NE. 0) THEN
3230	WRITE(,) ierr, fil
3231	STOP
3232	END IF
3233
3234	end if !is_master
3235
3236	call bcast(yg)
3237
3238	! photodissociation yield
3239
3240	fil = 'cross_sections/n2_ionef_schunknagy2000.txt'
3241
3242	if (is_master) then
3243
3244	OPEN(UNIT=kin,FILE=fil,STATUS='old')
3245
3246	n = 19
3247	read(kin,*)
3248	DO i = 1, n
3249	READ(kin,*) xl, xu, y2(i)
3250	x2(i) = (xl + xu)/2.
3251	y2(i) = 1. - y2(i)
3252	END DO
3253	CLOSE (kin)
3254
3255	CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
3256	CALL addpnt(x2,y2,kdata,n, 0.,0.)
3257	CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
3258	CALL addpnt(x2,y2,kdata,n, 1.e+38,1.)
3259	CALL inter2(nw,wl,yieldn2,n,x2,y2,ierr)
3260	IF (ierr .NE. 0) THEN
3261	WRITE(,) ierr, fil
3262	STOP
3263	END IF
3264
3265	end if !is_master
3266
3267	call bcast(yieldn2)
3268
3269	end subroutine rdxsn2
3270
3271	!==============================================================================
3272
3273	subroutine setalb(nw,wl,albedo)
3274
3275	!-----------------------------------------------------------------------------*
3276	!= PURPOSE: =*
3277	!= Set the albedo of the surface. The albedo is assumed to be Lambertian, =*
3278	!= i.e., the reflected light is isotropic, and idependt of the direction =*
3279	!= of incidence of light. Albedo can be chosen to be wavelength dependent. =*
3280	!-----------------------------------------------------------------------------*
3281	!= PARAMETERS: =*
3282	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
3283	!= wavelength grid =*
3284	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
3285	!= working wavelength grid =*
3286	!= ALBEDO - REAL, surface albedo at each specified wavelength (O)=*
3287	!-----------------------------------------------------------------------------*
3288
3289	implicit none
3290
3291	! input: (wavelength working grid data)
3292
3293	INTEGER nw
3294	REAL wl(nw)
3295
3296	! output:
3297
3298	REAL albedo(nw)
3299
3300	! local:
3301
3302	INTEGER iw
3303	REAL alb
3304
3305	! 0.015: mean value from clancy et al., icarus, 49-63, 1999.
3306
3307	alb = 0.015
3308
3309	do iw = 1, nw - 1
3310	albedo(iw) = alb
3311	end do
3312
3313	end subroutine setalb
3314
3315	end module photolysis_mod

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

Download in other formats:

Original Format