Context Navigation

photolysis_mod.F90 @ 3608

Last change on this file since 3608 was 3597, checked in by flefevre, 13 days ago

Updated S2 ultraviolet cross-section from Frandsen (personal communication, 2024).
The cross-section file can be downloaded from

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

and must be put in the INPUT/cross_sections directory.

File size: 99.9 KB

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

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

Download in other formats:

Original Format