Context Navigation

photolysis_mod.F90 @ 3803

Last change on this file since 3803 was 3755, checked in by flefevre, 7 months ago

Implementation of Cl2SO2 chemistry

Cl2SO2 absorption cross-section file must be downloaded from:

https://owncloud.latmos.ipsl.fr/index.php/s/bkTjmcm6EdoxMnT

and placed in INPUT/cross_sections directory.

File size: 109.9 KB

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

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

Download in other formats:

Original Format