Context Navigation

photolysis_mod.F90 @ 3567

Last change on this file since 3567 was 3530, checked in by flefevre, 3 weeks ago

Implementation of temperature-dependent SO absorption cross-sections from Heays et al., 2023.

The SO cross-section file must be downloaded from the link below and put in the /cross_sections directory:

https://owncloud.latmos.ipsl.fr/index.php/s/3xLzwn7G587sTdB

File size: 99.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_150, xsso_250 ! 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_150,xsso_250)
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,xsso_150,xsso_250)
2330
2331	!-----------------------------------------------------------------------------*
2332	!= PURPOSE: =*
2333	!= Read SO cross-sections =*
2334	!= Phillips (1981) or Heays et al. (2023) =*
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
2348	integer :: nw ! number of wavelength grid points
2349	real, dimension(nw) :: wl ! lower wavelength for each interval
2350
2351	! output
2352
2353	real, dimension(nw) :: xsso_150, xsso_250 ! so cross-sections (cm2)
2354
2355	! local
2356
2357	integer, parameter :: kdata = 1000
2358	integer :: i, l, n, n1, n2, ierr, iopt
2359	integer :: kin, kout ! input/output logical units
2360	real, dimension(nw) :: yg, yg1, yg2
2361	real, dimension(kdata) :: x1, y1, x2, y2
2362	real :: dummy
2363	real, parameter :: deltax = 1.e-4
2364	character*100 fil
2365
2366	kin = 10
2367
2368	! iopt = 1 : Phillips (1981)
2369	! iopt = 2 : Heays et al., Molecular Physics, 2023
2370
2371	iopt = 2
2372
2373	if (iopt == 1) then
2374
2375	! cross sections from Phillips (1981)
2376
2377	fil = 'cross_sections/so_marcq.txt'
2378
2379	if (is_master) then
2380
2381	print*, 'section efficace SO: ', fil
2382
2383	n = 851
2384	OPEN(kin,FILE=fil,STATUS='OLD')
2385	DO i = 1,n
2386	READ(kin,*) x1(i), dummy, dummy, dummy, dummy, y1(i)
2387	y1(i) = y1(i)*0.88 ! from Mills, PhD, 1998
2388	ENDDO
2389	CLOSE(kin)
2390
2391	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2392	CALL addpnt(x1,y1,kdata,n, 0.,0.)
2393	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2394	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
2395
2396	do i = 1,n
2397	print*, x1(i), y1(i)
2398	end do
2399
2400	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2401
2402	IF (ierr .NE. 0) THEN
2403	WRITE(,) ierr, fil
2404	STOP
2405	ENDIF
2406
2407	DO l = 1, nw-1
2408	xsso_150(l) = yg(l)
2409	xsso_250(l) = yg(l)
2410	END DO
2411
2412	endif !is_master
2413
2414	else if (iopt == 2) then
2415
2416	! cross sections from Heays et al. (2023)
2417	! values given at 150 K and 250 K
2418
2419	fil = 'cross_sections/SO_heays_2023.txt'
2420
2421	if (is_master) then
2422
2423	print*, 'section efficace SO: ', fil
2424
2425	n1 = 525
2426	n2 = n1
2427	OPEN(kin,FILE=fil,STATUS='OLD')
2428
2429	DO i = 1, 3
2430	READ(kin,*)
2431	ENDDO
2432
2433	DO i = 1,n1
2434	READ(kin,*) x1(i), dummy, dummy, dummy, dummy, dummy, &
2435	dummy, dummy, dummy, y1(i), y2(i)
2436	x2(i) = x1(i)
2437	ENDDO
2438	CLOSE(kin)
2439
2440	CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
2441	CALL addpnt(x1,y1,kdata,n1, 0.,0.)
2442	CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
2443	CALL addpnt(x1,y1,kdata,n1, 1.e+38,0.)
2444
2445	CALL inter2(nw,wl,yg1,n1,x1,y1,ierr)
2446
2447	IF (ierr .NE. 0) THEN
2448	WRITE(,) ierr, fil
2449	STOP
2450	ENDIF
2451
2452	CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
2453	CALL addpnt(x2,y2,kdata,n2, 0.,0.)
2454	CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
2455	CALL addpnt(x2,y2,kdata,n2, 1.e+38,0.)
2456
2457	CALL inter2(nw,wl,yg2,n2,x2,y2,ierr)
2458
2459	IF (ierr .NE. 0) THEN
2460	WRITE(,) ierr, fil
2461	STOP
2462	ENDIF
2463
2464	DO l = 1, nw-1
2465	xsso_150(l) = yg1(l)
2466	xsso_250(l) = yg2(l)
2467	END DO
2468
2469	endif !is_master
2470
2471	end if ! iopt
2472
2473	call bcast(yg)
2474
2475	end subroutine rdxsso
2476
2477	!==============================================================================
2478
2479	subroutine rdxsso3(nw, wl, yg)
2480
2481	!-----------------------------------------------------------------------------*
2482	!= PURPOSE: =*
2483	!= Read SO3 cross-sections =*
2484	!= Composite section 140-294 nm -> Pintze al. [2003] =*
2485	!= 296-330 nm -> Burkholder al. [1997] =*
2486	!-----------------------------------------------------------------------------*
2487	!= PARAMETERS: =*
2488	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
2489	!= wavelength grid =*
2490	!-----------------------------------------------------------------------------*
2491
2492	USE mod_phys_lmdz_para, ONLY: is_master
2493	USE mod_phys_lmdz_transfert_para, ONLY: bcast
2494
2495	IMPLICIT NONE
2496
2497	! input
2498
2499	integer :: nw ! number of wavelength grid points
2500	real, dimension(nw) :: wl ! lower wavelength for each interval
2501
2502	! output
2503
2504	real, dimension(nw) :: yg ! SO3 cross-sections (cm2)
2505
2506	! local
2507
2508	real, parameter :: deltax = 1.e-4
2509	integer, parameter :: kdata = 200
2510	real, dimension(kdata) :: x1, y1
2511	integer :: i, n, ierr
2512	character*100 fil
2513	integer :: kin, kout ! input/output logical units
2514
2515	kin = 10
2516
2517	!*** cross sections
2518
2519	fil = 'cross_sections/so3_composite.txt'
2520	print*, 'section efficace SO3: ', fil
2521
2522	if(is_master) then
2523
2524	n = 96
2525	OPEN(kin,FILE=fil,STATUS='OLD')
2526	DO i = 1,16
2527	READ(kin,*)
2528	ENDDO
2529	DO i = 1,n
2530	READ(kin,*) x1(i), y1(i)
2531	ENDDO
2532	CLOSE(kin)
2533
2534	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2535	CALL addpnt(x1,y1,kdata,n, 0.,0.)
2536	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2537	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
2538
2539	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2540
2541	IF (ierr .NE. 0) THEN
2542	WRITE(,) ierr, fil
2543	STOP
2544	ENDIF
2545
2546	endif !is_master
2547
2548	call bcast(yg)
2549
2550	end subroutine rdxsso3
2551
2552	!==============================================================================
2553
2554	subroutine rdxsclo(nw, wl, yg)
2555
2556	!-----------------------------------------------------------------------------*
2557	!= PURPOSE: =*
2558	!= Read ClO cross-sections =*
2559	!= From Trollier and al. [1990] =*
2560	!-----------------------------------------------------------------------------*
2561	!= PARAMETERS: =*
2562	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
2563	!= wavelength grid =*
2564	!-----------------------------------------------------------------------------*
2565
2566	USE mod_phys_lmdz_para, ONLY: is_master
2567	USE mod_phys_lmdz_transfert_para, ONLY: bcast
2568
2569	IMPLICIT NONE
2570
2571	! input
2572
2573	integer :: nw ! number of wavelength grid points
2574	real, dimension(nw) :: wl ! lower wavelength for each interval
2575
2576	! output
2577
2578	real, dimension(nw) :: yg ! hcl cross-sections (cm2)
2579
2580	! local
2581
2582	real, parameter :: deltax = 1.e-4
2583	integer, parameter :: kdata = 100
2584	real, dimension(kdata) :: x1, y1
2585	integer :: i, n, ierr
2586	character*100 fil
2587	integer :: kin, kout ! input/output logical units
2588
2589	kin = 10
2590
2591	!*** cross sections
2592
2593	fil = 'cross_sections/clo_xs_trolier_1990.txt'
2594	print*, 'section efficace ClO: ', fil
2595
2596	if(is_master) then
2597
2598	n = 81
2599	OPEN(kin,FILE=fil,STATUS='OLD')
2600	DO i = 1,6
2601	READ(kin,*)
2602	ENDDO
2603	DO i = 1,n
2604	READ(kin,*) x1(i), y1(i)
2605	ENDDO
2606	CLOSE(kin)
2607
2608	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2609	CALL addpnt(x1,y1,kdata,n, 0.,0.)
2610	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2611	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
2612
2613	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2614
2615	IF (ierr .NE. 0) THEN
2616	WRITE(,) ierr, fil
2617	STOP
2618	ENDIF
2619
2620	endif !is_master
2621
2622	call bcast(yg)
2623
2624	end subroutine rdxsclo
2625
2626	!==============================================================================
2627
2628	subroutine rdxsocs(nw, wl, yg)
2629
2630	!-----------------------------------------------------------------------------*
2631	!= PURPOSE: =*
2632	!= Read OCS cross-sections =*
2633	!= JPL 2011 recommendation =*
2634	!-----------------------------------------------------------------------------*
2635	!= PARAMETERS: =*
2636	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
2637	!= wavelength grid =*
2638	!-----------------------------------------------------------------------------*
2639
2640	USE mod_phys_lmdz_para, ONLY: is_master
2641	USE mod_phys_lmdz_transfert_para, ONLY: bcast
2642
2643	IMPLICIT NONE
2644
2645	! input
2646
2647	integer :: nw ! number of wavelength grid points
2648	real, dimension(nw) :: wl ! lower wavelength for each interval
2649
2650	! output
2651
2652	real, dimension(nw) :: yg ! hcl cross-sections (cm2)
2653
2654	! local
2655
2656	real, parameter :: deltax = 1.e-4
2657	integer, parameter :: kdata = 100
2658	real, dimension(kdata) :: x1, y1
2659	integer :: i, n, ierr
2660	character*100 fil
2661	integer :: kin, kout ! input/output logical units
2662
2663	kin = 10
2664
2665	!*** cross sections from JPL [2006]
2666
2667	fil = 'cross_sections/ocs_cross_sections_jpl2011.txt'
2668	print*, 'section efficace OCS: ', fil
2669
2670	if(is_master) then
2671
2672	n = 40
2673	OPEN(kin,FILE=fil,STATUS='OLD')
2674	DO i = 1,4
2675	READ(kin,*)
2676	ENDDO
2677	DO i = 1,n
2678	READ(kin,*) x1(i), y1(i)
2679	ENDDO
2680	CLOSE(kin)
2681
2682	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2683	CALL addpnt(x1,y1,kdata,n, 0.,0.)
2684	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2685	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
2686
2687	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2688
2689	IF (ierr .NE. 0) THEN
2690	WRITE(,) ierr, fil
2691	STOP
2692	ENDIF
2693
2694	endif !is_master
2695
2696	call bcast(yg)
2697
2698	end subroutine rdxsocs
2699
2700	!==============================================================================
2701
2702	subroutine rdxscocl2(nw, wl, yg)
2703
2704	!-----------------------------------------------------------------------------*
2705	!= PURPOSE: =*
2706	!= Read COCl2 cross-sections =*
2707	!= JPL 2011 recommendation =*
2708	!-----------------------------------------------------------------------------*
2709	!= PARAMETERS: =*
2710	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
2711	!= wavelength grid =*
2712	!-----------------------------------------------------------------------------*
2713
2714	USE mod_phys_lmdz_para, ONLY: is_master
2715	USE mod_phys_lmdz_transfert_para, ONLY: bcast
2716
2717	IMPLICIT NONE
2718
2719	! input
2720
2721	integer :: nw ! number of wavelength grid points
2722	real, dimension(nw) :: wl ! lower wavelength for each interval
2723
2724	! output
2725
2726	real, dimension(nw) :: yg ! COCl2 cross-sections (cm2)
2727
2728	! local
2729
2730	real, parameter :: deltax = 1.e-4
2731	integer, parameter :: kdata = 100
2732	real, dimension(kdata) :: x1, y1
2733	integer :: i, n, ierr
2734	character*100 fil
2735	integer :: kin, kout ! input/output logical units
2736
2737	kin = 10
2738
2739	!*** cross sections from JPL [2011]
2740
2741	fil = 'cross_sections/cocl2_cross_sections_jpl2011.txt'
2742	print*, 'section efficace COCl2: ', fil
2743
2744	if(is_master) then
2745
2746	n = 53
2747	OPEN(kin,FILE=fil,STATUS='OLD')
2748	DO i = 1,4
2749	READ(kin,*)
2750	ENDDO
2751	DO i = 1,n
2752	READ(kin,*) x1(i), y1(i)
2753	ENDDO
2754	CLOSE(kin)
2755
2756	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2757	CALL addpnt(x1,y1,kdata,n, 0.,0.)
2758	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2759	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
2760
2761	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2762
2763	IF (ierr .NE. 0) THEN
2764	WRITE(,) ierr, fil
2765	STOP
2766	ENDIF
2767
2768	endif !is_master
2769
2770	call bcast(yg)
2771
2772	end subroutine rdxscocl2
2773
2774	!==============================================================================
2775
2776	subroutine rdxsh2so4(nw, wl, yg)
2777
2778	!-----------------------------------------------------------------------------*
2779	!= PURPOSE: =*
2780	!= Read H2SO4 cross-sections =*
2781	!= JPL 2006 recommendation =*
2782	!-----------------------------------------------------------------------------*
2783	!= PARAMETERS: =*
2784	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
2785	!= wavelength grid =*
2786	!-----------------------------------------------------------------------------*
2787
2788	USE mod_phys_lmdz_para, ONLY: is_master
2789	USE mod_phys_lmdz_transfert_para, ONLY: bcast
2790
2791	IMPLICIT NONE
2792
2793	! input
2794
2795	integer :: nw ! number of wavelength grid points
2796	real, dimension(nw) :: wl ! lower wavelength for each interval
2797
2798	! output
2799
2800	real, dimension(nw) :: yg ! h2so4 cross-sections (cm2)
2801
2802	! local
2803
2804	real, parameter :: deltax = 1.e-4
2805	integer, parameter :: kdata = 100
2806	real, dimension(kdata) :: x1, y1
2807	integer :: i, n, ierr
2808	character*100 fil
2809	integer :: kin, kout ! input/output logical units
2810
2811	kin = 10
2812
2813	!*** cross sections from JPL [2006]
2814
2815	fil = 'cross_sections/h2so4_cross_sections.txt'
2816	print*, 'section efficace h2so4: ', fil
2817
2818	if(is_master) then
2819
2820	n = 22
2821	OPEN(kin,FILE=fil,STATUS='OLD')
2822	DO i = 1,3
2823	READ(kin,*)
2824	ENDDO
2825	DO i = 1,n
2826	READ(kin,*) x1(i), y1(i)
2827	ENDDO
2828	CLOSE(kin)
2829
2830	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2831	CALL addpnt(x1,y1,kdata,n, 0.,0.)
2832	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2833	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
2834
2835	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2836
2837	IF (ierr .NE. 0) THEN
2838	WRITE(,) ierr, fil
2839	STOP
2840	ENDIF
2841
2842	endif !is_master
2843
2844	call bcast(yg)
2845
2846	end subroutine rdxsh2so4
2847
2848	!==============================================================================
2849
2850	subroutine rdxsh2(nw, wl, wc, yg, yieldh2)
2851
2852	!-----------------------------------------------------------------------------*
2853	!= PURPOSE: =*
2854	!= Read h2 cross-sections and photodissociation yield =*
2855	!-----------------------------------------------------------------------------*
2856	!= PARAMETERS: =*
2857	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
2858	!= wavelength grid =*
2859	!-----------------------------------------------------------------------------*
2860
2861	USE mod_phys_lmdz_para, ONLY: is_master
2862	USE mod_phys_lmdz_transfert_para, ONLY: bcast
2863
2864	implicit none
2865
2866	! input
2867
2868	integer :: nw ! number of wavelength grid points
2869	real, dimension(nw) :: wl, wc ! lower and central wavelength for each interval
2870
2871	! output
2872
2873	real, dimension(nw) :: yg ! h2 cross-sections (cm2)
2874	real, dimension(nw) :: yieldh2 ! photodissociation yield
2875
2876	! local
2877
2878	integer, parameter :: kdata = 1000
2879	real, parameter :: deltax = 1.e-4
2880	real, dimension(kdata) :: x1, y1, x2, y2
2881	real :: xl, xu
2882	integer :: i, iw, n, ierr
2883	integer :: kin, kout ! input/output logical units
2884	character*100 fil
2885
2886	kin = 10
2887
2888	! h2 cross sections
2889
2890	fil = 'cross_sections/h2secef.txt'
2891	print*, 'section efficace H2: ', fil
2892
2893	if(is_master) then
2894
2895	OPEN(kin,FILE=fil,STATUS='OLD')
2896
2897	n = 792
2898	read(kin,*) ! avoid first line with wavelength = 0.
2899	DO i = 1, n
2900	READ(kin,*) x1(i), y1(i)
2901	ENDDO
2902	CLOSE(kin)
2903
2904	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2905	CALL addpnt(x1,y1,kdata,n, 0.,0.)
2906	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2907	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
2908	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2909
2910	IF (ierr .NE. 0) THEN
2911	WRITE(,) ierr, fil
2912	STOP
2913	ENDIF
2914
2915	endif !is_master
2916
2917	call bcast(yg)
2918
2919	! photodissociation yield
2920
2921	fil = 'cross_sections/h2_ionef_schunknagy2000.txt'
2922
2923	if(is_master) then
2924
2925	OPEN(UNIT=kin,FILE=fil,STATUS='old')
2926
2927	n = 19
2928	read(kin,*)
2929	DO i = 1, n
2930	READ(kin,*) xl, xu, y2(i)
2931	x2(i) = (xl + xu)/2.
2932	y2(i) = max(1. - y2(i),0.)
2933	END DO
2934	CLOSE (kin)
2935
2936	CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
2937	CALL addpnt(x2,y2,kdata,n, 0.,0.)
2938	CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
2939	CALL addpnt(x2,y2,kdata,n, 1.e+38,1.)
2940	CALL inter2(nw,wl,yieldh2,n,x2,y2,ierr)
2941	IF (ierr .NE. 0) THEN
2942	WRITE(,) ierr, fil
2943	STOP
2944	ENDIF
2945
2946	endif !is_master
2947
2948	call bcast(yieldh2)
2949
2950	end subroutine rdxsh2
2951
2952	!==============================================================================
2953
2954	subroutine rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294,yldno2_248, yldno2_298)
2955
2956	!-----------------------------------------------------------------------------*
2957	!= PURPOSE: =*
2958	!= read and grid cross section + quantum yield for NO2 =*
2959	!= photolysis =*
2960	!= Jenouvrier et al., 1996 200-238 nm
2961	!= Vandaele et al., 1998 238-666 nm 220K and 294K
2962	!= quantum yield from jpl 2006
2963	!-----------------------------------------------------------------------------*
2964	!= PARAMETERS: =*
2965	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
2966	!= wavelength grid =*
2967	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
2968	!= working wavelength grid =*
2969	!= SQ - REAL, cross section x quantum yield (cm^2) for each (O)=*
2970	!= photolysis reaction defined, at each defined wavelength and =*
2971	!= at each defined altitude level =*
2972	!-----------------------------------------------------------------------------*
2973
2974	USE mod_phys_lmdz_para, ONLY: is_master
2975	USE mod_phys_lmdz_transfert_para, ONLY: bcast
2976
2977	implicit none
2978
2979	! input
2980
2981	integer :: nw ! number of wavelength grid points
2982	real, dimension(nw) :: wl ! lower and central wavelength for each interval
2983
2984	! output
2985
2986	real, dimension(nw) :: xsno2, xsno2_220, xsno2_294 ! no2 cross-sections (cm2)
2987	real, dimension(nw) :: yldno2_248, yldno2_298 ! quantum yields at 248-298 k
2988
2989	! local
2990
2991	integer, parameter :: kdata = 28000
2992	real, parameter :: deltax = 1.e-4
2993	real, dimension(kdata) :: x1, x2, x3, x4, x5, y1, y2, y3, y4, y5
2994	real, dimension(nw) :: yg1, yg2, yg3, yg4, yg5
2995	real :: dum, qy
2996	integer :: i, iw, n, n1, n2, n3, n4, n5, ierr
2997	character*100 fil
2998	integer :: kin, kout ! input/output logical units
2999
3000	kin = 10
3001
3002	!*************** NO2 photodissociation
3003
3004	! Jenouvrier 1996 + Vandaele 1998 (JPL 2006)
3005
3006	fil = 'cross_sections/no2_xs_jenouvrier.txt'
3007	print*, 'section efficace NO2: ', fil
3008
3009	if (is_master) then
3010
3011	OPEN(UNIT=kin,FILE=fil,status='old')
3012	DO i = 1, 3
3013	READ(kin,*)
3014	END DO
3015	n1 = 10001
3016	DO i = 1, n1
3017	READ(kin,*) x1(i), y1(i)
3018	end do
3019
3020	CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax), 0.)
3021	CALL addpnt(x1,y1,kdata,n1, 0., 0.)
3022	CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
3023	CALL addpnt(x1,y1,kdata,n1, 1.e+38, 0.)
3024	CALL inter2(nw,wl,yg1,n1,x1,y1,ierr)
3025
3026	end if !is_master
3027
3028	call bcast(yg1)
3029
3030	fil = 'cross_sections/no2_xs_vandaele_294K.txt'
3031	print*, 'section efficace NO2: ', fil
3032
3033	if (is_master) then
3034
3035	OPEN(UNIT=kin,FILE=fil,status='old')
3036	DO i = 1, 3
3037	READ(kin,*)
3038	END DO
3039	n2 = 27993
3040	DO i = 1, n2
3041	READ(kin,*) x2(i), y2(i)
3042	end do
3043
3044	CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax), 0.)
3045	CALL addpnt(x2,y2,kdata,n2, 0., 0.)
3046	CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
3047	CALL addpnt(x2,y2,kdata,n2, 1.e+38, 0.)
3048	CALL inter2(nw,wl,yg2,n2,x2,y2,ierr)
3049
3050	end if !is_master
3051
3052	call bcast(yg2)
3053
3054	fil = 'cross_sections/no2_xs_vandaele_220K.txt'
3055	print*, 'section efficace NO2: ', fil
3056
3057	if (is_master) then
3058
3059	OPEN(UNIT=kin,FILE=fil,status='old')
3060	DO i = 1, 3
3061	READ(kin,*)
3062	END DO
3063	n3 = 27993
3064	do i = 1, n3
3065	READ(kin,*) x3(i), y3(i)
3066	end do
3067
3068	CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax), 0.)
3069	CALL addpnt(x3,y3,kdata,n3, 0., 0.)
3070	CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.)
3071	CALL addpnt(x3,y3,kdata,n3, 1.e+38, 0.)
3072	CALL inter2(nw,wl,yg3,n3,x3,y3,ierr)
3073
3074	do iw = 1, nw - 1
3075	xsno2(iw) = yg1(iw)
3076	xsno2_294(iw) = yg2(iw)
3077	xsno2_220(iw) = yg3(iw)
3078	end do
3079
3080	end if !is_master
3081
3082	call bcast(yg3)
3083	call bcast(xsno2)
3084	call bcast(xsno2_294)
3085	call bcast(xsno2_220)
3086
3087	! photodissociation efficiency from jpl 2006
3088
3089	fil = 'cross_sections/no2_yield_jpl2006.txt'
3090	print*, 'quantum yield NO2: ', fil
3091
3092	if (is_master) then
3093
3094	OPEN(UNIT=kin,FILE=fil,STATUS='old')
3095	DO i = 1, 5
3096	READ(kin,*)
3097	END DO
3098	n = 25
3099	n4 = n
3100	n5 = n
3101	DO i = 1, n
3102	READ(kin,*) x4(i), y4(i), y5(i)
3103	x5(i) = x4(i)
3104	END DO
3105	CLOSE(kin)
3106
3107	CALL addpnt(x4,y4,kdata,n4,x4(1)*(1.-deltax),y4(1))
3108	CALL addpnt(x4,y4,kdata,n4, 0.,y4(1))
3109	CALL addpnt(x4,y4,kdata,n4,x4(n4)*(1.+deltax), 0.)
3110	CALL addpnt(x4,y4,kdata,n4, 1.e+38, 0.)
3111	CALL inter2(nw,wl,yg4,n4,x4,y4,ierr)
3112	IF (ierr .NE. 0) THEN
3113	WRITE(,) ierr, fil
3114	STOP
3115	END IF
3116
3117	end if !is_master
3118
3119	call bcast(yg4)
3120
3121	if (is_master) then
3122
3123	CALL addpnt(x5,y5,kdata,n5,x5(1)*(1.-deltax),y5(1))
3124	CALL addpnt(x5,y5,kdata,n5, 0.,y5(1))
3125	CALL addpnt(x5,y5,kdata,n5,x5(n5)*(1.+deltax), 0.)
3126	CALL addpnt(x5,y5,kdata,n5, 1.e+38, 0.)
3127	CALL inter2(nw,wl,yg5,n5,x5,y5,ierr)
3128	IF (ierr .NE. 0) THEN
3129	WRITE(,) ierr, fil
3130	STOP
3131	END IF
3132
3133	do iw = 1, nw - 1
3134	yldno2_298(iw) = yg4(iw)
3135	yldno2_248(iw) = yg5(iw)
3136	end do
3137
3138	end if !is_master
3139
3140	call bcast(yg5)
3141	call bcast(yldno2_298)
3142	call bcast(yldno2_248)
3143
3144	end subroutine rdxsno2
3145
3146	!==============================================================================
3147
3148	subroutine rdxsno(nw, wl, yg, yieldno)
3149
3150	!-----------------------------------------------------------------------------*
3151	!= PURPOSE: =*
3152	!= Read NO cross-sections and photodissociation efficiency =*
3153	!= Lida et al 1986 (provided by Francisco Gonzalez-Galindo) =*
3154	!-----------------------------------------------------------------------------*
3155	!= PARAMETERS: =*
3156	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
3157	!= wavelength grid =*
3158	!-----------------------------------------------------------------------------*
3159
3160	USE mod_phys_lmdz_para, ONLY: is_master
3161	USE mod_phys_lmdz_transfert_para, ONLY: bcast
3162
3163	implicit none
3164
3165	! input
3166
3167	integer :: nw ! number of wavelength grid points
3168	real, dimension(nw) :: wl ! lower wavelength for each interval
3169
3170	! output
3171
3172	real, dimension(nw) :: yg ! no cross-sections (cm2)
3173	real, dimension(nw) :: yieldno ! no photodissociation efficiency
3174
3175	! local
3176
3177	integer, parameter :: kdata = 110
3178	real, parameter :: deltax = 1.e-4
3179	real, dimension(kdata) :: x1, y1, x2, y2
3180	integer :: i, iw, n, ierr
3181	character*100 fil
3182	integer :: kin, kout ! input/output logical units
3183
3184	kin = 10
3185
3186	! no cross-sections
3187
3188	fil = 'cross_sections/no_xs_francisco.txt'
3189	print*, 'section efficace NO: ', fil
3190
3191	if (is_master) then
3192
3193	OPEN(kin,FILE=fil,STATUS='OLD')
3194
3195	n = 99
3196	DO i = 1, n
3197	READ(kin,*) x1(i), y1(i)
3198	END DO
3199	CLOSE(kin)
3200
3201	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
3202	CALL addpnt(x1,y1,kdata,n, 0.,0.)
3203	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
3204	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
3205
3206	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
3207	IF (ierr .NE. 0) THEN
3208	WRITE(,) ierr, fil
3209	STOP
3210	END IF
3211
3212	end if !is_master
3213
3214	call bcast(yg)
3215
3216	! photodissociation yield
3217
3218	fil = 'cross_sections/noefdis.txt'
3219
3220	if (is_master) then
3221
3222	OPEN(UNIT=kin,FILE=fil,STATUS='old')
3223
3224	n = 33
3225	DO i = 1, n
3226	READ(kin,*) x2(n-i+1), y2(n-i+1)
3227	END DO
3228	CLOSE (kin)
3229
3230	CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
3231	CALL addpnt(x2,y2,kdata,n, 0.,0.)
3232	CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
3233	CALL addpnt(x2,y2,kdata,n, 1.e+38,1.)
3234	CALL inter2(nw,wl,yieldno,n,x2,y2,ierr)
3235	IF (ierr .NE. 0) THEN
3236	WRITE(,) ierr, fil
3237	STOP
3238	END IF
3239
3240	end if !is_master
3241
3242	call bcast(yieldno)
3243
3244	end subroutine rdxsno
3245
3246	!==============================================================================
3247
3248	subroutine rdxsn2(nw, wl, yg, yieldn2)
3249
3250	!-----------------------------------------------------------------------------*
3251	!= PURPOSE: =*
3252	!= Read n2 cross-sections and photodissociation yield =*
3253	!-----------------------------------------------------------------------------*
3254	!= PARAMETERS: =*
3255	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
3256	!= wavelength grid =*
3257	!-----------------------------------------------------------------------------*
3258
3259	USE mod_phys_lmdz_para, ONLY: is_master
3260	USE mod_phys_lmdz_transfert_para, ONLY: bcast
3261
3262	implicit none
3263
3264	! input
3265
3266	integer :: nw ! number of wavelength grid points
3267	real, dimension(nw) :: wl ! lower wavelength for each interval
3268
3269	! output
3270
3271	real, dimension(nw) :: yg ! n2 cross-sections (cm2)
3272	real, dimension(nw) :: yieldn2 ! n2 photodissociation yield
3273
3274	! local
3275
3276	integer, parameter :: kdata = 1100
3277	real, parameter :: deltax = 1.e-4
3278	real, dimension(kdata) :: x1, y1, x2, y2
3279	real :: xl, xu
3280	integer :: i, iw, n, ierr
3281	integer :: kin, kout ! input/output logical units
3282	character*100 fil
3283
3284	kin = 10
3285
3286	! n2 cross sections
3287
3288	fil = 'cross_sections/n2secef_01nm.txt'
3289	print*, 'section efficace N2: ', fil
3290
3291	if (is_master) then
3292
3293	OPEN(kin,FILE=fil,STATUS='OLD')
3294
3295	n = 1020
3296	DO i = 1, n
3297	READ(kin,*) x1(i), y1(i)
3298	END DO
3299	CLOSE(kin)
3300
3301	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
3302	CALL addpnt(x1,y1,kdata,n, 0.,0.)
3303	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
3304	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
3305
3306	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
3307
3308	IF (ierr .NE. 0) THEN
3309	WRITE(,) ierr, fil
3310	STOP
3311	END IF
3312
3313	end if !is_master
3314
3315	call bcast(yg)
3316
3317	! photodissociation yield
3318
3319	fil = 'cross_sections/n2_ionef_schunknagy2000.txt'
3320
3321	if (is_master) then
3322
3323	OPEN(UNIT=kin,FILE=fil,STATUS='old')
3324
3325	n = 19
3326	read(kin,*)
3327	DO i = 1, n
3328	READ(kin,*) xl, xu, y2(i)
3329	x2(i) = (xl + xu)/2.
3330	y2(i) = 1. - y2(i)
3331	END DO
3332	CLOSE (kin)
3333
3334	CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
3335	CALL addpnt(x2,y2,kdata,n, 0.,0.)
3336	CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
3337	CALL addpnt(x2,y2,kdata,n, 1.e+38,1.)
3338	CALL inter2(nw,wl,yieldn2,n,x2,y2,ierr)
3339	IF (ierr .NE. 0) THEN
3340	WRITE(,) ierr, fil
3341	STOP
3342	END IF
3343
3344	end if !is_master
3345
3346	call bcast(yieldn2)
3347
3348	end subroutine rdxsn2
3349
3350	!==============================================================================
3351
3352	subroutine setalb(nw,wl,albedo)
3353
3354	!-----------------------------------------------------------------------------*
3355	!= PURPOSE: =*
3356	!= Set the albedo of the surface. The albedo is assumed to be Lambertian, =*
3357	!= i.e., the reflected light is isotropic, and idependt of the direction =*
3358	!= of incidence of light. Albedo can be chosen to be wavelength dependent. =*
3359	!-----------------------------------------------------------------------------*
3360	!= PARAMETERS: =*
3361	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
3362	!= wavelength grid =*
3363	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
3364	!= working wavelength grid =*
3365	!= ALBEDO - REAL, surface albedo at each specified wavelength (O)=*
3366	!-----------------------------------------------------------------------------*
3367
3368	implicit none
3369
3370	! input: (wavelength working grid data)
3371
3372	INTEGER nw
3373	REAL wl(nw)
3374
3375	! output:
3376
3377	REAL albedo(nw)
3378
3379	! local:
3380
3381	INTEGER iw
3382	REAL alb
3383
3384	! 0.015: mean value from clancy et al., icarus, 49-63, 1999.
3385
3386	alb = 0.015
3387
3388	do iw = 1, nw - 1
3389	albedo(iw) = alb
3390	end do
3391
3392	end subroutine setalb
3393
3394	end module photolysis_mod

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

Download in other formats:

Original Format