source: trunk/LMDZ.VENUS/libf/phyvenus/photolysis_mod.F90 @ 3556

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

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

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

​https://owncloud.latmos.ipsl.fr/index.php/s/3xLzwn7G587sTdB
File size: 99.2 KB
RevLine 
[2780]1module photolysis_mod
2
3  implicit none
4
5! photolysis
6
[2925]7  integer, save :: nphot = 24             ! number of photolysis
[2780]8
9!$OMP THREADPRIVATE(nphot)
10
[2925]11  integer, parameter :: nabs  = 21        ! number of absorbing gases
[2780]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)
[3530]42  real, dimension(nw), save :: xsso_150, xsso_250                     ! so absorption cross-section (cm2)
[2780]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)
[2925]46  real, dimension(nw), save :: xss2                                   ! s2 absorption cross-section (cm2)
[2780]47  real, dimension(nw), save :: xscocl2                                ! cocl2 absorption cross-section (cm2)
48  real, dimension(nw), save :: xsh2so4                                ! h2so4 absorption cross-section (cm2)
[2795]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)
[2780]58  real, dimension(nw), save :: albedo                                 ! surface albedo
59
[3530]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)
[2780]62
63contains
64
65subroutine 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
[3530]125  call rdxsso(nw,wl,xsso_150,xsso_250)
[2780]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
[2925]143! read and grid s2 cross-sections
144
145  call rdxss2(nw,wl,xss2)
146
[2780]147! read and grid h2so4 cross-sections
148
149  call rdxsh2so4(nw,wl,xsh2so4)
150
151! read and grid h2 cross-sections
152
[2795]153  call rdxsh2(nw,wl,wc,xsh2,yieldh2)
[2780]154
155! read and grid no cross-sections
156
[2795]157  call rdxsno(nw,wl,xsno,yieldno)
[2780]158
159! read and grid no2 cross-sections
160
[2795]161  call rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294,yldno2_248,yldno2_298)
[2780]162
163! read and grid n2 cross-sections
164
[2795]165  call rdxsn2(nw,wl,xsn2,yieldn2)
[2780]166
167! read and grid hdo cross-sections
168
[2795]169  call rdxshdo(nw,wl,xshdo)
[2780]170
171! set surface albedo
172
173  call setalb(nw,wl,albedo)
174
175end 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
[2836]576#include "clesphys.h" ! fixed_euv_value
577
[2780]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
[2836]591      integer :: iw, nhead, ihead, n, i, ierr, kin1, kin2
[2780]592
593      real, parameter :: deltax = 1.e-4
[2836]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
[2780]606
[2836]607      kin1 = 10    ! input logical unit
608      kin2 = 11    ! input logical unit
609     
[2780]610! select desired extra-terrestrial solar irradiance, using msun:
611
[2836]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)
[2780]614!      Thuillier et al., Adv. Space. Res., 34, 256-261, 2004
[2836]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
[2780]617
[2836]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
[2780]623
[2836]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
[2780]631
[2836]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
[2780]676
[2836]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)
[2780]687
[2836]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)
[2780]707
[2836]708        IF (ierr .NE. 0) THEN
709          WRITE(*,*) ierr, fil1, fil2
710          STOP
711        ENDIF
[2780]712
713!     convert to photon.s-1.nm-1.cm-2
714!     5.039e11 = 1.e-4*1e-9/(hc = 6.62e-34*2.998e8)
715
[2836]716        DO iw = 1, nw-1
717          f(iw) = yg0(iw)*wc(iw)*5.039e11
718        ENDDO
[2780]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
[2925]1536      subroutine rdxss2(nw, wl, yg)
1537
1538!-----------------------------------------------------------------------------*
1539!=  PURPOSE:                                                                 =*
1540!=  Read S2 molecular absorption cross section.  Re-grid data to match      =*
1541!=  specified wavelength working grid.                                       =*
1542!-----------------------------------------------------------------------------*
1543!=  PARAMETERS:                                                              =*
1544!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
1545!=           wavelength grid                                                 =*
1546!=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
1547!=           working wavelength grid                                         =*
1548!=  YG     - REAL, molecular absoprtion cross section (cm^2) of H2O at    (O)=*
1549!=           each specified wavelength                                       =*
1550!-----------------------------------------------------------------------------*
1551
1552      USE mod_phys_lmdz_para, ONLY: is_master
1553      USE mod_phys_lmdz_transfert_para, ONLY: bcast
1554
1555      IMPLICIT NONE
1556
1557!     input
1558
1559      integer :: nw               ! number of wavelength grid points
1560      real, dimension(nw) :: wl   ! lower and central wavelength for each interval
1561
1562!     output
1563
1564      real, dimension(nw) :: yg   ! h2o cross-sections (cm2)
1565
1566!     local
1567
1568      integer, parameter :: kdata = 1500
1569      real, parameter :: deltax = 1.e-4
1570      REAL x1(kdata)
1571      REAL y1(kdata)
1572      INTEGER ierr, dummy
1573      INTEGER i, n
1574
1575      CHARACTER*100 fil
1576      integer :: kin, kout ! input/output logical units
1577
1578      kin = 10
1579
1580      fil = 'cross_sections/s2_millsBestEst_1560_5059.txt'
1581      print*, 'section efficace S2: ', fil
1582
1583      if(is_master) then
1584
1585      OPEN(UNIT=kin,FILE=fil,STATUS='old')
1586
1587      n = 1203
1588      DO i = 1, n
1589        READ(kin,*) x1(i), y1(i)
1590        x1(i) = x1(i)/10
1591      END DO
1592      CLOSE (kin)
1593
1594      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1595      CALL addpnt(x1,y1,kdata,n,          0.,0.)
1596      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1597      CALL addpnt(x1,y1,kdata,n,      1.e+38,0.)
1598      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
1599      IF (ierr .NE. 0) THEN
1600         WRITE(*,*) ierr, fil
1601         STOP
1602      ENDIF
1603
1604      endif !is_master
1605
1606      call bcast(yg)
1607
1608      end subroutine rdxss2
1609
1610!==============================================================================
1611
[2780]1612      subroutine rdxsh2o(nw, wl, yg)
1613
1614!-----------------------------------------------------------------------------*
1615!=  PURPOSE:                                                                 =*
1616!=  Read H2O molecular absorption cross section.  Re-grid data to match      =*
1617!=  specified wavelength working grid.                                       =*
1618!-----------------------------------------------------------------------------*
1619!=  PARAMETERS:                                                              =*
1620!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
1621!=           wavelength grid                                                 =*
1622!=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
1623!=           working wavelength grid                                         =*
1624!=  YG     - REAL, molecular absoprtion cross section (cm^2) of H2O at    (O)=*
1625!=           each specified wavelength                                       =*
1626!-----------------------------------------------------------------------------*
1627
1628      USE mod_phys_lmdz_para, ONLY: is_master
1629      USE mod_phys_lmdz_transfert_para, ONLY: bcast
1630
1631      IMPLICIT NONE
1632
1633!     input
1634
1635      integer :: nw               ! number of wavelength grid points
1636      real, dimension(nw) :: wl   ! lower and central wavelength for each interval
1637
1638!     output
1639
1640      real, dimension(nw) :: yg   ! h2o cross-sections (cm2)
1641
1642!     local
1643
1644      integer, parameter :: kdata = 500
1645      real, parameter :: deltax = 1.e-4
1646      REAL x1(kdata)
1647      REAL y1(kdata)
1648      INTEGER ierr
1649      INTEGER i, n
1650      CHARACTER*100 fil
1651      integer :: kin, kout ! input/output logical units
1652
1653      kin = 10
1654
1655      fil = 'cross_sections/h2o_composite_250K.txt'
1656      print*, 'section efficace H2O: ', fil
1657
1658      if(is_master) then
1659
1660      OPEN(UNIT=kin,FILE=fil,STATUS='old')
1661
1662      DO i = 1,26
1663         read(kin,*)
1664      END DO
1665
1666      n = 420
1667      DO i = 1, n
1668         READ(kin,*) x1(i), y1(i)
1669      END DO
1670      CLOSE (kin)
1671
1672      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1673      CALL addpnt(x1,y1,kdata,n,          0.,0.)
1674      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1675      CALL addpnt(x1,y1,kdata,n,      1.e+38,0.)
1676      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
1677      IF (ierr .NE. 0) THEN
1678         WRITE(*,*) ierr, fil
1679         STOP
1680      ENDIF
1681
1682      endif !is_master
1683
1684      call bcast(yg)
1685     
1686      end subroutine rdxsh2o
1687
1688!==============================================================================
1689
[2795]1690       subroutine rdxshdo(nw, wl, yg)
[2780]1691
1692!-----------------------------------------------------------------------------*
1693!=  PURPOSE:                                                                 =*
1694!=  Read HDO molecular absorption cross section.  Re-grid data to match      =*
1695!=  specified wavelength working grid.                                       =*
1696!-----------------------------------------------------------------------------*
1697!=  PARAMETERS:                                                              =*
1698!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
1699!=           wavelength grid                                                 =*
1700!=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
1701!=           working wavelength grid                                         =*
1702!=  YG     - REAL, molecular absoprtion cross section (cm^2) of HDO at    (O)=*
1703!=           each specified wavelength                                       =*
1704!-----------------------------------------------------------------------------*
1705
[2795]1706      USE mod_phys_lmdz_para, ONLY: is_master
1707      USE mod_phys_lmdz_transfert_para, ONLY: bcast
[2780]1708
[2795]1709      IMPLICIT NONE
[2780]1710
1711!     input
1712
[2795]1713      integer :: nw               ! number of wavelength grid points
1714      real, dimension(nw) :: wl   ! lower and central wavelength for each interval
[2780]1715
1716!     output
1717
[2795]1718      real, dimension(nw) :: yg   ! hdo cross-sections (cm2)
[2780]1719
1720!     local
1721
[2795]1722      integer, parameter :: kdata = 900
1723      real, parameter :: deltax = 1.e-4
1724      REAL x1(kdata)
1725      REAL y1(kdata)
1726      INTEGER ierr
1727      INTEGER i, n
1728      CHARACTER*100 fil
1729      integer :: kin, kout ! input/output logical units
[2780]1730
[2795]1731      kin = 10
[2780]1732
[2799]1733      fil = 'cross_sections/hdo_composite_295K.txt'
[2795]1734      print*, 'section efficace HDO: ', fil
[2780]1735
[2795]1736      if(is_master) then
[2780]1737
[2795]1738      OPEN(UNIT=kin,FILE=fil,STATUS='old')
[2780]1739
[2795]1740      DO i = 1,17
1741         read(kin,*)
1742      END DO
[2780]1743
[2795]1744      n = 806
1745      DO i = 1, n
1746         READ(kin,*) x1(i), y1(i)
1747      END DO
1748      CLOSE (kin)
[2780]1749
[2795]1750      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1751      CALL addpnt(x1,y1,kdata,n,          0.,0.)
1752      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1753      CALL addpnt(x1,y1,kdata,n,      1.e+38,0.)
1754      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
1755      IF (ierr .NE. 0) THEN
1756         WRITE(*,*) ierr, fil
1757         STOP
1758      ENDIF
[2780]1759
[2795]1760      endif !is_master
[2780]1761
[2795]1762      call bcast(yg)
[2780]1763     
[2795]1764      end subroutine rdxshdo
[2780]1765
1766!==============================================================================
1767
1768      subroutine rdxsh2o2(nw, wl, xsh2o2)
1769
1770!-----------------------------------------------------------------------------*
1771!=  PURPOSE:                                                                 =*
1772!=  Read and grid H2O2 cross-sections
1773!=         H2O2 + hv -> 2 OH                                                 =*
1774!=  Cross section:  Schuergers and Welge, Z. Naturforsch. 23a (1968) 1508    =*
1775!=                  from 125 to 185 nm, then JPL97 from 190 to 350 nm.       =*
1776!-----------------------------------------------------------------------------*
1777!=  PARAMETERS:                                                              =*
1778!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
1779!=           wavelength grid                                                 =*
1780!=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
1781!=           working wavelength grid                                         =*
1782!-----------------------------------------------------------------------------*
1783
1784      USE mod_phys_lmdz_para, ONLY: is_master
1785      USE mod_phys_lmdz_transfert_para, ONLY: bcast
1786
1787      implicit none
1788
1789!     input
1790
1791      integer :: nw               ! number of wavelength grid points
1792      real, dimension(nw) :: wl   ! lower wavelength for each interval
1793
1794!     output
1795
1796      real, dimension(nw) :: xsh2o2   ! h2o2 cross-sections (cm2)
1797
1798!     local
1799
1800      real, parameter :: deltax = 1.e-4
1801      integer, parameter :: kdata = 100
1802      real, dimension(kdata) :: x1, y1
1803      real, dimension(nw)    :: yg
1804      integer :: i, ierr, iw, n, idum
1805      integer :: kin, kout ! input/output logical units
1806      character*100 fil
1807
1808      kin = 10
1809
1810!     read cross-sections
1811
1812      fil = 'cross_sections/h2o2_composite.txt'
1813      print*, 'section efficace H2O2: ', fil
1814
1815      if(is_master) then
1816
1817      OPEN(kin,FILE=fil,STATUS='OLD')
1818      READ(kin,*) idum,n
1819      DO i = 1, idum-2
1820         READ(kin,*)
1821      ENDDO
1822      DO i = 1, n
1823         READ(kin,*) x1(i), y1(i)
1824      ENDDO
1825      CLOSE (kin)
1826
1827      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1828      CALL addpnt(x1,y1,kdata,n,               0.,0.)
1829      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1830      CALL addpnt(x1,y1,kdata,n,           1.e+38,0.)
1831      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
1832      IF (ierr .NE. 0) THEN
1833         WRITE(*,*) ierr, fil
1834         STOP
1835      ENDIF
1836
1837      DO iw = 1, nw - 1
1838         xsh2o2(iw) = yg(iw)
1839      END DO
1840
1841      endif !is_master
1842
1843      call bcast(xsh2o2)
1844
1845      end subroutine rdxsh2o2
1846
1847!==============================================================================
1848
1849      subroutine rdxsho2(nw, wl, yg)
1850
1851!-----------------------------------------------------------------------------*
1852!=  PURPOSE:                                                                 =*
1853!=  Read ho2 cross-sections                                                  =*
1854!=  JPL 2006 recommendation                                                  =*
1855!-----------------------------------------------------------------------------*
1856!=  PARAMETERS:                                                              =*
1857!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
1858!=           wavelength grid                                                 =*
1859!-----------------------------------------------------------------------------*
1860
1861      USE mod_phys_lmdz_para, ONLY: is_master
1862      USE mod_phys_lmdz_transfert_para, ONLY: bcast
1863
1864      IMPLICIT NONE
1865
1866!     input
1867
1868      integer :: nw               ! number of wavelength grid points
1869      real, dimension(nw) :: wl   ! lower wavelength for each interval
1870
1871!     output
1872
1873      real, dimension(nw) :: yg   ! ho2 cross-sections (cm2)
1874
1875!     local
1876
1877      real, parameter :: deltax = 1.e-4
1878      integer, parameter :: kdata = 100
1879      real, dimension(kdata) :: x1, y1
1880      integer :: i, n, ierr
1881      character*100 fil
1882      integer :: kin, kout ! input/output logical units
1883
1884      kin = 10
1885
1886!*** cross sections from Sander et al. [2003]
1887
1888      fil = 'cross_sections/ho2_jpl2003.txt'
1889      print*, 'section efficace HO2: ', fil
1890
1891      if(is_master) then
1892
1893      OPEN(kin,FILE=fil,STATUS='OLD')
1894      READ(kin,*) n
1895      DO i = 1, n
1896        READ(kin,*) x1(i), y1(i)
1897      ENDDO
1898      CLOSE(kin)
1899
1900      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1901      CALL addpnt(x1,y1,kdata,n,          0.,0.)
1902      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1903      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
1904
1905      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
1906 
1907      IF (ierr .NE. 0) THEN
1908        WRITE(*,*) ierr, fil
1909        STOP
1910      ENDIF
1911
1912      endif !is_master
1913
1914      call bcast(yg)
1915 
1916      end subroutine rdxsho2
1917
1918!==============================================================================
1919
1920      subroutine rdxshcl(nw, wl, yg)
1921
1922!-----------------------------------------------------------------------------*
1923!=  PURPOSE:                                                                 =*
1924!=  Read hcl cross-sections                                                  =*
1925!=  JPL 2006 recommendation                                                  =*
1926!-----------------------------------------------------------------------------*
1927!=  PARAMETERS:                                                              =*
1928!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
1929!=           wavelength grid                                                 =*
1930!-----------------------------------------------------------------------------*
1931
1932      USE mod_phys_lmdz_para, ONLY: is_master
1933      USE mod_phys_lmdz_transfert_para, ONLY: bcast
1934
1935      IMPLICIT NONE
1936
1937!     input
1938
1939      integer :: nw               ! number of wavelength grid points
1940      real, dimension(nw) :: wl   ! lower wavelength for each interval
1941
1942!     output
1943
1944      real, dimension(nw) :: yg   ! hcl cross-sections (cm2)
1945
1946!     local
1947
1948      real, parameter :: deltax = 1.e-4
1949      integer, parameter :: kdata = 100
1950      real, dimension(kdata) :: x1, y1
1951      integer :: i, n, ierr
1952      character*100 fil
1953      integer :: kin, kout ! input/output logical units
1954
1955      kin = 10
1956
1957!*** cross sections from JPL [2006]
1958
1959      fil = 'cross_sections/hcl_jpl2006.txt'
1960      print*, 'section efficace HCl: ', fil
1961
1962      if(is_master) then
1963      n = 31
1964      OPEN(kin,FILE=fil,STATUS='OLD')
1965      DO i = 1,4
1966         READ(kin,*)
1967      ENDDO   
1968      DO i = 1,n
1969        READ(kin,*) x1(i), y1(i)
1970      ENDDO
1971      CLOSE(kin)
1972
1973      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1974      CALL addpnt(x1,y1,kdata,n,          0.,0.)
1975      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1976      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
1977      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
1978
1979 
1980      IF (ierr .NE. 0) THEN
1981        WRITE(*,*) ierr, fil
1982        STOP
1983      ENDIF
1984
1985      endif !is_master
1986
1987      call bcast(yg)
1988 
1989      end subroutine rdxshcl
1990
1991!==============================================================================
1992
1993      subroutine rdxscl2(nw, wl, yg)
1994
1995!-----------------------------------------------------------------------------*
1996!=  PURPOSE:                                                                 =*
1997!=  Read cl2 cross-sections                                                  =*
1998!=  JPL 2006 recommendation                                                  =*
1999!-----------------------------------------------------------------------------*
2000!=  PARAMETERS:                                                              =*
2001!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
2002!=           wavelength grid                                                 =*
2003!-----------------------------------------------------------------------------*
2004
2005      USE mod_phys_lmdz_para, ONLY: is_master
2006      USE mod_phys_lmdz_transfert_para, ONLY: bcast
2007
2008      IMPLICIT NONE
2009
2010!     input
2011
2012      integer :: nw               ! number of wavelength grid points
2013      real, dimension(nw) :: wl   ! lower wavelength for each interval
2014
2015!     output
2016
2017      real, dimension(nw) :: yg   ! cl2 cross-sections (cm2)
2018
2019!     local
2020
2021      real, parameter :: deltax = 1.e-4
2022      integer, parameter :: kdata = 100
2023      real, dimension(kdata) :: x1, y1
2024      integer :: i, n, ierr
2025      character*100 fil
2026      integer :: kin, kout ! input/output logical units
2027
2028      kin = 10
2029
2030!*** cross sections from JPL [2006]
2031
2032      fil = 'cross_sections/cl2_jpl2006.txt'
2033      print*, 'section efficace Cl2: ', fil
2034
2035      if(is_master) then
2036
2037      n = 30
2038      OPEN(kin,FILE=fil,STATUS='OLD')
2039      DO i = 1,4
2040         READ(kin,*)
2041      ENDDO   
2042      DO i = 1,n
2043        READ(kin,*) x1(i), y1(i)
2044      ENDDO
2045      CLOSE(kin)
2046
2047      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2048      CALL addpnt(x1,y1,kdata,n,          0.,0.)
2049      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2050      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
2051
2052      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2053 
2054      IF (ierr .NE. 0) THEN
2055        WRITE(*,*) ierr, fil
2056        STOP
2057      ENDIF
2058
2059      endif !is_master
2060
2061      call bcast(yg)
2062 
2063      end subroutine rdxscl2
2064
2065!==============================================================================
2066
2067      subroutine rdxshocl(nw, wl, yg)
2068
2069!-----------------------------------------------------------------------------*
2070!=  PURPOSE:                                                                 =*
2071!=  Read hocl cross-sections                                                  =*
2072!=  JPL 2000 recommendation                                                  =*
2073!-----------------------------------------------------------------------------*
2074!=  PARAMETERS:                                                              =*
2075!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
2076!=           wavelength grid                                                 =*
2077!-----------------------------------------------------------------------------*
2078
2079      USE mod_phys_lmdz_para, ONLY: is_master
2080      USE mod_phys_lmdz_transfert_para, ONLY: bcast
2081
2082      IMPLICIT NONE
2083
2084!     input
2085
2086      integer :: nw               ! number of wavelength grid points
2087      real, dimension(nw) :: wl   ! lower wavelength for each interval
2088
2089!     output
2090
2091      real, dimension(nw) :: yg   ! hocl cross-sections (cm2)
2092
2093!     local
2094
2095      real, parameter :: deltax = 1.e-4
2096      integer, parameter :: kdata = 200
2097      real, dimension(kdata) :: x1, y1
2098      integer :: i, n, ierr
2099      character*100 fil
2100      integer :: kin, kout ! input/output logical units
2101
2102      kin = 10
2103
2104!*** cross sections from JPL [2000]
2105
2106      fil = 'cross_sections/HOCl_jpl2000.txt'
2107      print*, 'section efficace HOCl: ', fil
2108
2109      if(is_master) then
2110
2111         n = 111
2112      OPEN(kin,FILE=fil,STATUS='OLD')
2113      DO i = 1,6
2114         READ(kin,*)
2115      ENDDO   
2116      DO i = 1,111
2117        READ(kin,*) y1(i)
2118        x1(i) = 200 + real(i-1)*2.
2119        y1(i) = y1(i) * 1.E-20
2120      ENDDO
2121      CLOSE(kin)
2122
2123      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2124      CALL addpnt(x1,y1,kdata,n,          0.,0.)
2125      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2126      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
2127
2128      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2129 
2130      IF (ierr .NE. 0) THEN
2131        WRITE(*,*) ierr, fil
2132        STOP
2133      ENDIF
2134
2135      endif !is_master
2136
2137      call bcast(yg)
2138 
2139      end subroutine rdxshocl
2140
2141!==============================================================================
2142
2143      subroutine rdxsso2(nw,wl,xsso2_200,xsso2_298,xsso2_360)
2144
2145!-----------------------------------------------------------------------------*
2146!=  PURPOSE:                                                                 =*
2147!=  Read and grid SO2 absorption cross-sections and photodissociation yield   =*
2148!-----------------------------------------------------------------------------*
2149!=  PARAMETERS:                                                              =*
2150!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
2151!=           wavelength grid                                                 =*
2152!=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
2153!=           working wavelength grid                                         =*
2154!=  XSSO2  - REAL, molecular absoprtion cross section (cm^2) of SO2 at    (O)=*
2155!=           each specified wavelength                                       =*
2156!-----------------------------------------------------------------------------*
2157
2158      USE mod_phys_lmdz_para, ONLY: is_master
2159      USE mod_phys_lmdz_transfert_para, ONLY: bcast
2160
2161      implicit none
2162
2163!     input
2164
2165      integer :: nw               ! number of wavelength grid points
2166      real, dimension(nw) :: wl   ! lower and central wavelength for each interval
2167
2168!     output
2169
2170      real, dimension(nw) :: xsso2_200, xsso2_298, xsso2_360 ! so2 cross-sections (cm2)
2171
2172!     local
2173
2174      integer, parameter :: kdata = 55000
2175      real, parameter :: deltax = 1.e-4
2176      real, dimension(kdata) :: x1, y1, y2, y3, xion, ion
2177      real, dimension(nw) :: yg
2178      real :: xl, xu
2179      integer :: ierr, i, l, n, n1, n2, n3, n4
2180      CHARACTER*100 fil
2181
2182      integer :: kin, kout ! input/ouput logical units
2183
2184      kin  = 10
2185      kout = 30
2186
2187!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2188!
2189!     SO2 absorption cross-sections
2190!
2191!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2192!
2193!     iopt = 1
2194!
2195!     200K: Wu et al. (2000) + Vandaele et al. (2009) + Hermans et al. (2009)
2196!           7 lignes d'en-tete et n1 = 50276
2197!           fichier: so2_composite_200K.txt
2198!
2199!     298K: Wu et al. (2000) + Vandaele et al. (2009) + Hermans et al. (2009)
2200!           7 lignes d'en-tete et n2 = 49833
2201!           fichier: so2_composite_298K.txt
2202!
2203!     360K: Wu et al. (2000) + Vandaele et al. (2009) + Hermans et al. (2009)
2204!           7 lignes d'en-tete et n3 = 49261
2205!           fichier: so2_composite_360K.txt
2206!
2207!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2208
2209
2210      n1 = 50276
2211      n2 = 49833
2212      n3 = 46261
2213
2214!     200K:
2215
2216      fil = 'cross_sections/so2_composite_200K.txt'
2217      print*, 'section efficace SO2 195K: ', fil
2218
2219      if(is_master) then
2220
2221      OPEN(UNIT=kin,FILE=fil,STATUS='old')
2222      DO i = 1,7
2223         read(kin,*)
2224      END DO
2225
2226      DO i = 1, n1
2227         READ(kin,*) x1(i), y1(i)
2228      END DO
2229      CLOSE (kin)
2230
2231      CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
2232      CALL addpnt(x1,y1,kdata,n1,          0.,0.)
2233      CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
2234      CALL addpnt(x1,y1,kdata,n1,      1.e+38,0.)
2235      CALL inter2(nw,wl,yg,n1,x1,y1,ierr)
2236      IF (ierr .NE. 0) THEN
2237         WRITE(*,*) ierr, fil
2238         STOP
2239      ENDIF
2240     
2241      DO l = 1, nw-1
2242         xsso2_200(l) = yg(l)
2243      END DO
2244
2245      endif !is_master
2246
2247      call bcast(xsso2_200)
2248
2249!     298K:
2250
2251      fil = 'cross_sections/so2_composite_298K.txt'
2252      print*, 'section efficace SO2 295K: ', fil
2253
2254      if(is_master) then
2255
2256      OPEN(UNIT=kin,FILE=fil,STATUS='old')
2257      DO i = 1,7
2258         read(kin,*)
2259      END DO
2260
2261      DO i = 1, n2
2262         READ(kin,*) x1(i), y1(i)
2263      END DO
2264      CLOSE (kin)
2265
2266      CALL addpnt(x1,y1,kdata,n2,x1(1)*(1.-deltax),0.)
2267      CALL addpnt(x1,y1,kdata,n2,          0.,0.)
2268      CALL addpnt(x1,y1,kdata,n2,x1(n2)*(1.+deltax),0.)
2269      CALL addpnt(x1,y1,kdata,n2,      1.e+38,0.)
2270      CALL inter2(nw,wl,yg,n2,x1,y1,ierr)
2271      IF (ierr .NE. 0) THEN
2272         WRITE(*,*) ierr, fil
2273         STOP
2274      ENDIF
2275
2276      DO l = 1, nw-1
2277         xsso2_298(l) = yg(l)
2278      END DO
2279
2280      endif !is_master
2281
2282      call bcast(xsso2_298)
2283
2284!     360K:
2285
2286      fil = 'cross_sections/so2_composite_360K.txt'
2287      print*, 'section efficace SO2 370K: ', fil
2288
2289      if(is_master) then
2290
2291      OPEN(UNIT=kin,FILE=fil,STATUS='old')
2292      DO i = 1,7
2293         read(kin,*)
2294      END DO
2295
2296      DO i = 1, n3
2297         READ(kin,*) x1(i), y1(i)
2298      END DO
2299      CLOSE (kin)
2300
2301      CALL addpnt(x1,y1,kdata,n3,x1(1)*(1.-deltax),0.)
2302      CALL addpnt(x1,y1,kdata,n3,          0.,0.)
2303      CALL addpnt(x1,y1,kdata,n3,x1(n3)*(1.+deltax),0.)
2304      CALL addpnt(x1,y1,kdata,n3,      1.e+38,0.)
2305      CALL inter2(nw,wl,yg,n3,x1,y1,ierr)
2306      IF (ierr .NE. 0) THEN
2307         WRITE(*,*) ierr, fil
2308         STOP
2309      ENDIF
2310
2311      DO l = 1, nw-1
2312         xsso2_360(l) = yg(l)
2313      END DO
2314
2315      endif !is_master
2316
2317      call bcast(xsso2_360)
2318
2319!     DO l = 1, nw-1
2320!        write(kout,*) wl(l), xsso2_200(l),
2321!    $                        xsso2_298(l),
2322!    $                        xsso2_360(l),
2323!     END DO
2324
2325      end subroutine rdxsso2
2326
2327!==============================================================================
2328
[3530]2329      subroutine rdxsso(nw,wl,xsso_150,xsso_250)
[2780]2330
2331!-----------------------------------------------------------------------------*
2332!=  PURPOSE:                                                                 =*
[3530]2333!=  Read SO cross-sections                                                   =*
2334!=  Phillips (1981) or Heays et al. (2023)                                   =*
[2780]2335!-----------------------------------------------------------------------------*
2336!=  PARAMETERS:                                                              =*
2337!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
2338!=           wavelength grid                                                 =*
2339!-----------------------------------------------------------------------------*
2340
2341      USE mod_phys_lmdz_para, ONLY: is_master
2342      USE mod_phys_lmdz_transfert_para, ONLY: bcast
2343
[3530]2344      implicit none
[2780]2345
2346!     input
[3530]2347
[2780]2348      integer :: nw               ! number of wavelength grid points
[3530]2349      real, dimension(nw) :: wl   ! lower wavelength for each interval
[2780]2350
2351!     output
2352
[3530]2353      real, dimension(nw) :: xsso_150, xsso_250 ! so cross-sections (cm2)
[2780]2354
2355!     local
[3530]2356
2357      integer, parameter :: kdata = 1000
2358      integer :: i, l, n, n1, n2, ierr, iopt
2359      integer :: kin, kout ! input/output logical units
2360      real, dimension(nw) :: yg, yg1, yg2
2361      real, dimension(kdata) :: x1, y1, x2, y2
[2780]2362      real :: dummy
2363      real, parameter :: deltax = 1.e-4
2364      character*100 fil
2365
2366      kin = 10
2367
[3530]2368!     iopt = 1 : Phillips (1981)                       
2369!     iopt = 2 : Heays et al., Molecular Physics, 2023
[2780]2370
[3530]2371      iopt = 2
[2780]2372
[3530]2373      if (iopt == 1) then
[2780]2374
[3530]2375!        cross sections from Phillips (1981)
2376
2377         fil = 'cross_sections/so_marcq.txt'
2378
2379         if (is_master) then
2380
2381         print*, 'section efficace SO: ', fil
2382
[2780]2383         n = 851
[3530]2384         OPEN(kin,FILE=fil,STATUS='OLD')
2385         DO i = 1,n
2386            READ(kin,*) x1(i), dummy, dummy, dummy, dummy, y1(i)
2387            y1(i) = y1(i)*0.88  ! from Mills, PhD, 1998
2388         ENDDO
2389         CLOSE(kin)
[2780]2390
[3530]2391         CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2392         CALL addpnt(x1,y1,kdata,n,          0.,0.)
2393         CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2394         CALL addpnt(x1,y1,kdata,n,        1E38,0.)
[2780]2395
[3530]2396         do i = 1,n
2397            print*, x1(i), y1(i)
2398         end do
2399
2400         CALL inter2(nw,wl,yg,n,x1,y1,ierr)
[2780]2401 
[3530]2402         IF (ierr .NE. 0) THEN
2403           WRITE(*,*) ierr, fil
2404           STOP
2405         ENDIF
[2780]2406
[3530]2407         DO l = 1, nw-1
2408            xsso_150(l) = yg(l)
2409            xsso_250(l) = yg(l)
2410         END DO
[2780]2411
[3530]2412         endif !is_master
2413
2414      else if (iopt == 2) then
2415
2416!        cross sections from Heays et al. (2023)
2417!        values given at 150 K and 250 K
2418
2419         fil = 'cross_sections/SO_heays_2023.txt'
2420
2421         if (is_master) then
2422
2423         print*, 'section efficace SO: ', fil
2424         
2425         n1 = 525
2426         n2 = n1
2427         OPEN(kin,FILE=fil,STATUS='OLD')
2428
2429         DO i = 1, 3
2430            READ(kin,*)
2431         ENDDO
2432
2433         DO i = 1,n1
2434            READ(kin,*) x1(i), dummy, dummy, dummy, dummy, dummy, &
2435                        dummy, dummy, dummy, y1(i), y2(i)
2436            x2(i) = x1(i)
2437         ENDDO
2438         CLOSE(kin)
2439
2440         CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
2441         CALL addpnt(x1,y1,kdata,n1,               0.,0.)
2442         CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
2443         CALL addpnt(x1,y1,kdata,n1,           1.e+38,0.)
2444
2445         CALL inter2(nw,wl,yg1,n1,x1,y1,ierr)
2446
2447         IF (ierr .NE. 0) THEN
2448            WRITE(*,*) ierr, fil
2449            STOP
2450         ENDIF
2451
2452         CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
2453         CALL addpnt(x2,y2,kdata,n2,               0.,0.)
2454         CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
2455         CALL addpnt(x2,y2,kdata,n2,           1.e+38,0.)
2456
2457         CALL inter2(nw,wl,yg2,n2,x2,y2,ierr)
2458
2459         IF (ierr .NE. 0) THEN
2460            WRITE(*,*) ierr, fil
2461            STOP
2462         ENDIF
2463
2464         DO l = 1, nw-1
2465            xsso_150(l) = yg1(l)
2466            xsso_250(l) = yg2(l)
2467         END DO
2468
2469         endif !is_master
2470
2471      end if ! iopt
2472
[2780]2473      call bcast(yg)
2474 
2475      end subroutine rdxsso
2476
2477!==============================================================================
2478
2479      subroutine rdxsso3(nw, wl, yg)
2480
2481!-----------------------------------------------------------------------------*
2482!=  PURPOSE:                                                                 =*
2483!=  Read SO3 cross-sections                                                  =*
2484!=  Composite section 140-294 nm -> Pintze al. [2003]                        =*
2485!=                    296-330 nm -> Burkholder al. [1997]                    =*
2486!-----------------------------------------------------------------------------*
2487!=  PARAMETERS:                                                              =*
2488!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
2489!=           wavelength grid                                                 =*
2490!-----------------------------------------------------------------------------*
2491
2492      USE mod_phys_lmdz_para, ONLY: is_master
2493      USE mod_phys_lmdz_transfert_para, ONLY: bcast
2494
2495      IMPLICIT NONE
2496
2497!     input
2498
2499      integer :: nw               ! number of wavelength grid points
2500      real, dimension(nw) :: wl   ! lower wavelength for each interval
2501
2502!     output
2503
2504      real, dimension(nw) :: yg   ! SO3 cross-sections (cm2)
2505
2506!     local
2507
2508      real, parameter :: deltax = 1.e-4
2509      integer, parameter :: kdata = 200
2510      real, dimension(kdata) :: x1, y1
2511      integer :: i, n, ierr
2512      character*100 fil
2513      integer :: kin, kout ! input/output logical units
2514
2515      kin = 10
2516
2517!*** cross sections
2518
2519      fil = 'cross_sections/so3_composite.txt'
2520      print*, 'section efficace SO3: ', fil
2521
2522      if(is_master) then
2523
2524         n = 96
2525      OPEN(kin,FILE=fil,STATUS='OLD')
2526      DO i = 1,16
2527         READ(kin,*)
2528      ENDDO   
2529      DO i = 1,n
2530        READ(kin,*) x1(i), y1(i)
2531      ENDDO
2532      CLOSE(kin)
2533
2534      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2535      CALL addpnt(x1,y1,kdata,n,          0.,0.)
2536      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2537      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
2538
2539      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2540 
2541      IF (ierr .NE. 0) THEN
2542        WRITE(*,*) ierr, fil
2543        STOP
2544      ENDIF
2545
2546      endif !is_master
2547
2548      call bcast(yg)
2549 
2550      end subroutine rdxsso3
2551
2552!==============================================================================
2553
2554      subroutine rdxsclo(nw, wl, yg)
2555
2556!-----------------------------------------------------------------------------*
2557!=  PURPOSE:                                                                 =*
2558!=  Read ClO cross-sections                                                  =*
2559!=  From Trollier and al. [1990]                                             =*
2560!-----------------------------------------------------------------------------*
2561!=  PARAMETERS:                                                              =*
2562!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
2563!=           wavelength grid                                                 =*
2564!-----------------------------------------------------------------------------*
2565
2566      USE mod_phys_lmdz_para, ONLY: is_master
2567      USE mod_phys_lmdz_transfert_para, ONLY: bcast
2568
2569      IMPLICIT NONE
2570
2571!     input
2572
2573      integer :: nw               ! number of wavelength grid points
2574      real, dimension(nw) :: wl   ! lower wavelength for each interval
2575
2576!     output
2577
2578      real, dimension(nw) :: yg   ! hcl cross-sections (cm2)
2579
2580!     local
2581
2582      real, parameter :: deltax = 1.e-4
2583      integer, parameter :: kdata = 100
2584      real, dimension(kdata) :: x1, y1
2585      integer :: i, n, ierr
2586      character*100 fil
2587      integer :: kin, kout ! input/output logical units
2588
2589      kin = 10
2590
2591!*** cross sections
2592
2593      fil = 'cross_sections/clo_xs_trolier_1990.txt'
2594      print*, 'section efficace ClO: ', fil
2595
2596      if(is_master) then
2597
2598         n = 81
2599      OPEN(kin,FILE=fil,STATUS='OLD')
2600      DO i = 1,6
2601         READ(kin,*)
2602      ENDDO   
2603      DO i = 1,n
2604        READ(kin,*) x1(i), y1(i)
2605      ENDDO
2606      CLOSE(kin)
2607
2608      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2609      CALL addpnt(x1,y1,kdata,n,          0.,0.)
2610      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2611      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
2612
2613      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2614 
2615      IF (ierr .NE. 0) THEN
2616        WRITE(*,*) ierr, fil
2617        STOP
2618      ENDIF
2619
2620      endif !is_master
2621
2622      call bcast(yg)
2623 
2624      end subroutine rdxsclo
2625
2626!==============================================================================
2627
2628      subroutine rdxsocs(nw, wl, yg)
2629
2630!-----------------------------------------------------------------------------*
2631!=  PURPOSE:                                                                 =*
2632!=  Read OCS cross-sections                                                  =*
2633!=  JPL 2011 recommendation                                                  =*
2634!-----------------------------------------------------------------------------*
2635!=  PARAMETERS:                                                              =*
2636!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
2637!=           wavelength grid                                                 =*
2638!-----------------------------------------------------------------------------*
2639
2640      USE mod_phys_lmdz_para, ONLY: is_master
2641      USE mod_phys_lmdz_transfert_para, ONLY: bcast
2642
2643      IMPLICIT NONE
2644
2645!     input
2646
2647      integer :: nw               ! number of wavelength grid points
2648      real, dimension(nw) :: wl   ! lower wavelength for each interval
2649
2650!     output
2651
2652      real, dimension(nw) :: yg   ! hcl cross-sections (cm2)
2653
2654!     local
2655
2656      real, parameter :: deltax = 1.e-4
2657      integer, parameter :: kdata = 100
2658      real, dimension(kdata) :: x1, y1
2659      integer :: i, n, ierr
2660      character*100 fil
2661      integer :: kin, kout ! input/output logical units
2662
2663      kin = 10
2664
2665!*** cross sections from JPL [2006]
2666
2667      fil = 'cross_sections/ocs_cross_sections_jpl2011.txt'
2668      print*, 'section efficace OCS: ', fil
2669
2670      if(is_master) then
2671
2672         n = 40
2673      OPEN(kin,FILE=fil,STATUS='OLD')
2674      DO i = 1,4
2675         READ(kin,*)
2676      ENDDO   
2677      DO i = 1,n
2678        READ(kin,*) x1(i), y1(i)
2679      ENDDO
2680      CLOSE(kin)
2681
2682      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2683      CALL addpnt(x1,y1,kdata,n,          0.,0.)
2684      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2685      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
2686
2687      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2688 
2689      IF (ierr .NE. 0) THEN
2690        WRITE(*,*) ierr, fil
2691        STOP
2692      ENDIF
2693
2694      endif !is_master
2695
2696      call bcast(yg)
2697 
2698      end subroutine rdxsocs
2699
2700!==============================================================================
2701
2702      subroutine rdxscocl2(nw, wl, yg)
2703
2704!-----------------------------------------------------------------------------*
2705!=  PURPOSE:                                                                 =*
2706!=  Read COCl2 cross-sections                                                =*
2707!=  JPL 2011 recommendation                                                  =*
2708!-----------------------------------------------------------------------------*
2709!=  PARAMETERS:                                                              =*
2710!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
2711!=           wavelength grid                                                 =*
2712!-----------------------------------------------------------------------------*
2713
2714      USE mod_phys_lmdz_para, ONLY: is_master
2715      USE mod_phys_lmdz_transfert_para, ONLY: bcast
2716
2717      IMPLICIT NONE
2718
2719!     input
2720
2721      integer :: nw               ! number of wavelength grid points
2722      real, dimension(nw) :: wl   ! lower wavelength for each interval
2723
2724!     output
2725
2726      real, dimension(nw) :: yg   ! COCl2 cross-sections (cm2)
2727
2728!     local
2729
2730      real, parameter :: deltax = 1.e-4
2731      integer, parameter :: kdata = 100
2732      real, dimension(kdata) :: x1, y1
2733      integer :: i, n, ierr
2734      character*100 fil
2735      integer :: kin, kout ! input/output logical units
2736
2737      kin = 10
2738
2739!*** cross sections from JPL [2011]
2740
2741      fil = 'cross_sections/cocl2_cross_sections_jpl2011.txt'
2742      print*, 'section efficace COCl2: ', fil
2743
2744      if(is_master) then
2745
2746         n = 53
2747      OPEN(kin,FILE=fil,STATUS='OLD')
2748      DO i = 1,4
2749         READ(kin,*)
2750      ENDDO   
2751      DO i = 1,n
2752        READ(kin,*) x1(i), y1(i)
2753      ENDDO
2754      CLOSE(kin)
2755
2756      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2757      CALL addpnt(x1,y1,kdata,n,          0.,0.)
2758      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2759      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
2760
2761      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2762 
2763      IF (ierr .NE. 0) THEN
2764        WRITE(*,*) ierr, fil
2765        STOP
2766      ENDIF
2767
2768      endif !is_master
2769
2770      call bcast(yg)
2771 
2772      end subroutine rdxscocl2
2773
2774!==============================================================================
2775
2776      subroutine rdxsh2so4(nw, wl, yg)
2777
2778!-----------------------------------------------------------------------------*
2779!=  PURPOSE:                                                                 =*
2780!=  Read H2SO4 cross-sections                                                  =*
2781!=  JPL 2006 recommendation                                                  =*
2782!-----------------------------------------------------------------------------*
2783!=  PARAMETERS:                                                              =*
2784!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
2785!=           wavelength grid                                                 =*
2786!-----------------------------------------------------------------------------*
2787
2788      USE mod_phys_lmdz_para, ONLY: is_master
2789      USE mod_phys_lmdz_transfert_para, ONLY: bcast
2790
2791      IMPLICIT NONE
2792
2793!     input
2794
2795      integer :: nw               ! number of wavelength grid points
2796      real, dimension(nw) :: wl   ! lower wavelength for each interval
2797
2798!     output
2799
2800      real, dimension(nw) :: yg   ! h2so4 cross-sections (cm2)
2801
2802!     local
2803
2804      real, parameter :: deltax = 1.e-4
2805      integer, parameter :: kdata = 100
2806      real, dimension(kdata) :: x1, y1
2807      integer :: i, n, ierr
2808      character*100 fil
2809      integer :: kin, kout ! input/output logical units
2810
2811      kin = 10
2812
2813!*** cross sections from JPL [2006]
2814
2815      fil = 'cross_sections/h2so4_cross_sections.txt'
2816      print*, 'section efficace h2so4: ', fil
2817
2818      if(is_master) then
2819
2820         n = 22
2821      OPEN(kin,FILE=fil,STATUS='OLD')
2822      DO i = 1,3
2823         READ(kin,*)
2824      ENDDO   
2825      DO i = 1,n
2826        READ(kin,*) x1(i), y1(i)
2827      ENDDO
2828      CLOSE(kin)
2829
2830      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2831      CALL addpnt(x1,y1,kdata,n,          0.,0.)
2832      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2833      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
2834
2835      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
2836 
2837      IF (ierr .NE. 0) THEN
2838        WRITE(*,*) ierr, fil
2839        STOP
2840      ENDIF
2841
2842      endif !is_master
2843
2844      call bcast(yg)
2845 
2846      end subroutine rdxsh2so4
2847
2848!==============================================================================
2849
[2795]2850       subroutine rdxsh2(nw, wl, wc, yg, yieldh2)
[2780]2851
2852!-----------------------------------------------------------------------------*
2853!=  PURPOSE:                                                                 =*
2854!=  Read h2 cross-sections and photodissociation yield                       =*
2855!-----------------------------------------------------------------------------*
2856!=  PARAMETERS:                                                              =*
2857!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
2858!=           wavelength grid                                                 =*
2859!-----------------------------------------------------------------------------*
2860
[2795]2861      USE mod_phys_lmdz_para, ONLY: is_master
2862      USE mod_phys_lmdz_transfert_para, ONLY: bcast
[2780]2863
[2795]2864      implicit none
[2780]2865
2866!     input
2867
[2795]2868      integer :: nw                   ! number of wavelength grid points
2869      real, dimension(nw) :: wl, wc   ! lower and central wavelength for each interval
[2780]2870
2871!     output
2872
[2795]2873      real, dimension(nw) :: yg        ! h2 cross-sections (cm2)
2874      real, dimension(nw) :: yieldh2   ! photodissociation yield
[2780]2875
2876!     local
2877
[2795]2878      integer, parameter :: kdata = 1000
2879      real, parameter :: deltax = 1.e-4
2880      real, dimension(kdata) :: x1, y1, x2, y2
2881      real :: xl, xu
2882      integer :: i, iw, n, ierr
2883      integer :: kin, kout ! input/output logical units
2884      character*100 fil
[2780]2885
[2795]2886      kin = 10
[2780]2887
2888!     h2 cross sections
2889
[2799]2890      fil = 'cross_sections/h2secef.txt'
[2795]2891      print*, 'section efficace H2: ', fil
[2780]2892
[2795]2893      if(is_master) then
[2780]2894
[2795]2895      OPEN(kin,FILE=fil,STATUS='OLD')
[2780]2896
[2795]2897      n = 792
2898      read(kin,*) ! avoid first line with wavelength = 0.
2899      DO i = 1, n
2900        READ(kin,*) x1(i), y1(i)
2901      ENDDO
2902      CLOSE(kin)
[2780]2903
[2795]2904      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
2905      CALL addpnt(x1,y1,kdata,n,          0.,0.)
2906      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
2907      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
2908      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
[2780]2909 
[2795]2910      IF (ierr .NE. 0) THEN
2911        WRITE(*,*) ierr, fil
2912        STOP
2913      ENDIF
[2780]2914
[2795]2915      endif !is_master
[2780]2916
[2795]2917      call bcast(yg)
[2780]2918 
2919!     photodissociation yield
2920
[2799]2921      fil = 'cross_sections/h2_ionef_schunknagy2000.txt'
[2780]2922
[2795]2923      if(is_master) then
[2780]2924
[2795]2925      OPEN(UNIT=kin,FILE=fil,STATUS='old')
[2780]2926
[2795]2927      n = 19
2928      read(kin,*)
2929      DO i = 1, n
2930         READ(kin,*) xl, xu, y2(i)
2931         x2(i) = (xl + xu)/2.
2932         y2(i) = max(1. - y2(i),0.)
2933      END DO
2934      CLOSE (kin)
[2780]2935
[2795]2936      CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
2937      CALL addpnt(x2,y2,kdata,n,          0.,0.)
2938      CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
2939      CALL addpnt(x2,y2,kdata,n,      1.e+38,1.)
2940      CALL inter2(nw,wl,yieldh2,n,x2,y2,ierr)
2941      IF (ierr .NE. 0) THEN
2942        WRITE(*,*) ierr, fil
2943        STOP
2944      ENDIF
[2780]2945
[2795]2946      endif !is_master
[2780]2947
[2795]2948      call bcast(yieldh2)
[2780]2949
[2795]2950      end subroutine rdxsh2
[2780]2951
2952!==============================================================================
[2795]2953
2954      subroutine rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294,yldno2_248, yldno2_298)
2955
[2780]2956!-----------------------------------------------------------------------------*
2957!=  PURPOSE:                                                                 =*
2958!=  read and grid cross section + quantum yield for NO2                      =*
2959!=  photolysis                                                               =*
2960!=  Jenouvrier et al., 1996  200-238 nm
2961!=  Vandaele et al., 1998    238-666 nm 220K and 294K
2962!=  quantum yield from jpl 2006
2963!-----------------------------------------------------------------------------*
2964!=  PARAMETERS:                                                              =*
2965!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
2966!=           wavelength grid                                                 =*
2967!=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
2968!=           working wavelength grid                                         =*
2969!=  SQ     - REAL, cross section x quantum yield (cm^2) for each          (O)=*
2970!=           photolysis reaction defined, at each defined wavelength and     =*
2971!=           at each defined altitude level                                  =*
2972!-----------------------------------------------------------------------------*
[2795]2973
2974      USE mod_phys_lmdz_para, ONLY: is_master
2975      USE mod_phys_lmdz_transfert_para, ONLY: bcast
2976
2977      implicit none
2978
[2780]2979!     input
[2795]2980
2981      integer :: nw               ! number of wavelength grid points
2982      real, dimension(nw) :: wl   ! lower and central wavelength for each interval
2983     
[2780]2984!     output
[2795]2985
2986      real, dimension(nw) :: xsno2, xsno2_220, xsno2_294 ! no2 cross-sections (cm2)
2987      real, dimension(nw) :: yldno2_248, yldno2_298      ! quantum yields at 248-298 k
2988
[2780]2989!     local
[2795]2990
2991      integer, parameter :: kdata = 28000
2992      real, parameter :: deltax = 1.e-4
2993      real, dimension(kdata) :: x1, x2, x3, x4, x5, y1, y2, y3, y4, y5
2994      real, dimension(nw) :: yg1, yg2, yg3, yg4, yg5
2995      real :: dum, qy
2996      integer :: i, iw, n, n1, n2, n3, n4, n5, ierr
2997      character*100 fil
2998      integer :: kin, kout ! input/output logical units
2999
3000      kin = 10
3001
[2780]3002!*************** NO2 photodissociation
[2795]3003
[2780]3004!  Jenouvrier 1996 + Vandaele 1998 (JPL 2006) 
[2795]3005
3006      fil = 'cross_sections/no2_xs_jenouvrier.txt'
3007      print*, 'section efficace NO2: ', fil
3008
3009      if (is_master) then
3010
3011      OPEN(UNIT=kin,FILE=fil,status='old')
3012      DO i = 1, 3
3013         READ(kin,*)
3014      END DO
3015      n1 = 10001
3016      DO i = 1, n1
3017         READ(kin,*) x1(i), y1(i)
3018      end do
3019
3020      CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax), 0.)
3021      CALL addpnt(x1,y1,kdata,n1,               0., 0.)
3022      CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
3023      CALL addpnt(x1,y1,kdata,n1,           1.e+38, 0.)
3024      CALL inter2(nw,wl,yg1,n1,x1,y1,ierr)
3025
3026      end if !is_master
3027
3028      call bcast(yg1)
3029
3030      fil = 'cross_sections/no2_xs_vandaele_294K.txt'
3031      print*, 'section efficace NO2: ', fil
3032
3033      if (is_master) then
3034
3035      OPEN(UNIT=kin,FILE=fil,status='old')
3036      DO i = 1, 3
3037        READ(kin,*)
3038      END DO
3039      n2 = 27993
3040      DO i = 1, n2
3041         READ(kin,*) x2(i), y2(i)
3042      end do
3043
3044      CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax), 0.)
3045      CALL addpnt(x2,y2,kdata,n2,               0., 0.)
3046      CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
3047      CALL addpnt(x2,y2,kdata,n2,           1.e+38, 0.)
3048      CALL inter2(nw,wl,yg2,n2,x2,y2,ierr)
3049
3050      end if !is_master
3051
3052      call bcast(yg2)
3053
3054      fil = 'cross_sections/no2_xs_vandaele_220K.txt'
3055      print*, 'section efficace NO2: ', fil
3056
3057      if (is_master) then
3058
3059      OPEN(UNIT=kin,FILE=fil,status='old')
3060      DO i = 1, 3
3061         READ(kin,*)
3062      END DO
3063      n3 = 27993
3064      do i = 1, n3
3065         READ(kin,*) x3(i), y3(i)
3066      end do
3067
3068      CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax), 0.)
3069      CALL addpnt(x3,y3,kdata,n3,               0., 0.)
3070      CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.)
3071      CALL addpnt(x3,y3,kdata,n3,           1.e+38, 0.)
3072      CALL inter2(nw,wl,yg3,n3,x3,y3,ierr)
3073
3074      do iw = 1, nw - 1
3075         xsno2(iw)     = yg1(iw)
3076         xsno2_294(iw) = yg2(iw)
3077         xsno2_220(iw) = yg3(iw)
3078      end do
3079
3080      end if !is_master
3081
3082      call bcast(yg3)
3083      call bcast(xsno2)
3084      call bcast(xsno2_294)
3085      call bcast(xsno2_220)
3086
[2780]3087!     photodissociation efficiency from jpl 2006
[2795]3088
3089      fil = 'cross_sections/no2_yield_jpl2006.txt'
3090      print*, 'quantum yield NO2: ', fil
3091
3092      if (is_master) then
3093
3094      OPEN(UNIT=kin,FILE=fil,STATUS='old')
3095      DO i = 1, 5
3096         READ(kin,*)
3097      END DO
3098      n = 25
3099      n4 = n
3100      n5 = n
3101      DO i = 1, n
3102         READ(kin,*) x4(i), y4(i), y5(i)
3103         x5(i) = x4(i)
3104      END DO
3105      CLOSE(kin)
3106
3107      CALL addpnt(x4,y4,kdata,n4,x4(1)*(1.-deltax),y4(1))
3108      CALL addpnt(x4,y4,kdata,n4,               0.,y4(1))
3109      CALL addpnt(x4,y4,kdata,n4,x4(n4)*(1.+deltax),  0.)
3110      CALL addpnt(x4,y4,kdata,n4,           1.e+38,   0.)
3111      CALL inter2(nw,wl,yg4,n4,x4,y4,ierr)
3112      IF (ierr .NE. 0) THEN
3113         WRITE(*,*) ierr, fil
3114         STOP
3115      END IF
3116
3117      end if !is_master
3118
3119      call bcast(yg4)
3120
3121      if (is_master) then
3122
3123      CALL addpnt(x5,y5,kdata,n5,x5(1)*(1.-deltax),y5(1))
3124      CALL addpnt(x5,y5,kdata,n5,               0.,y5(1))
3125      CALL addpnt(x5,y5,kdata,n5,x5(n5)*(1.+deltax),  0.)
3126      CALL addpnt(x5,y5,kdata,n5,           1.e+38,   0.)
3127      CALL inter2(nw,wl,yg5,n5,x5,y5,ierr)
3128      IF (ierr .NE. 0) THEN
3129         WRITE(*,*) ierr, fil
3130         STOP
3131      END IF
3132
3133      do iw = 1, nw - 1
3134         yldno2_298(iw) = yg4(iw)
3135         yldno2_248(iw) = yg5(iw)
3136      end do
3137
3138      end if !is_master
3139
3140      call bcast(yg5)
3141      call bcast(yldno2_298)
3142      call bcast(yldno2_248)
3143     
3144      end subroutine rdxsno2
3145
[2780]3146!==============================================================================
3147
[2795]3148      subroutine rdxsno(nw, wl, yg, yieldno)
[2780]3149
3150!-----------------------------------------------------------------------------*
3151!=  PURPOSE:                                                                 =*
3152!=  Read NO cross-sections  and photodissociation efficiency                 =*
3153!=  Lida et al 1986 (provided by Francisco Gonzalez-Galindo)                 =*
3154!-----------------------------------------------------------------------------*
3155!=  PARAMETERS:                                                              =*
3156!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
3157!=           wavelength grid                                                 =*
3158!-----------------------------------------------------------------------------*
3159
[2795]3160      USE mod_phys_lmdz_para, ONLY: is_master
3161      USE mod_phys_lmdz_transfert_para, ONLY: bcast
[2780]3162
[2795]3163      implicit none
[2780]3164
3165!     input
3166
[2795]3167      integer :: nw               ! number of wavelength grid points
3168      real, dimension(nw) :: wl   ! lower wavelength for each interval
[2780]3169
3170!     output
3171
[2795]3172      real, dimension(nw) :: yg        ! no cross-sections (cm2)
3173      real, dimension(nw) :: yieldno   ! no photodissociation efficiency
[2780]3174
3175!     local
3176
[2795]3177      integer, parameter :: kdata = 110
3178      real, parameter :: deltax = 1.e-4
3179      real, dimension(kdata) :: x1, y1, x2, y2
3180      integer :: i, iw, n, ierr
3181      character*100 fil
3182      integer :: kin, kout ! input/output logical units
[2780]3183
[2795]3184      kin = 10
[2780]3185
3186!     no cross-sections
3187
[2795]3188      fil = 'cross_sections/no_xs_francisco.txt'
3189      print*, 'section efficace NO: ', fil
[2780]3190
[2795]3191      if (is_master) then
[2780]3192
[2795]3193      OPEN(kin,FILE=fil,STATUS='OLD')
[2780]3194
[2795]3195      n = 99
3196      DO i = 1, n
3197        READ(kin,*) x1(i), y1(i)
3198      END DO
3199      CLOSE(kin)
[2780]3200
[2795]3201      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
3202      CALL addpnt(x1,y1,kdata,n,          0.,0.)
3203      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
3204      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
[2780]3205
[2795]3206      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
3207      IF (ierr .NE. 0) THEN
3208         WRITE(*,*) ierr, fil
3209         STOP
3210      END IF
[2780]3211
[2795]3212      end if !is_master
[2780]3213
[2795]3214      call bcast(yg)
[2780]3215 
3216!     photodissociation yield
3217
[2795]3218      fil = 'cross_sections/noefdis.txt'
[2780]3219
[2795]3220      if (is_master) then
[2780]3221
[2795]3222      OPEN(UNIT=kin,FILE=fil,STATUS='old')
[2780]3223
[2795]3224      n = 33
3225      DO i = 1, n
3226         READ(kin,*) x2(n-i+1), y2(n-i+1)
3227      END DO
3228      CLOSE (kin)
[2780]3229
[2795]3230      CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
3231      CALL addpnt(x2,y2,kdata,n,          0.,0.)
3232      CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
3233      CALL addpnt(x2,y2,kdata,n,      1.e+38,1.)
3234      CALL inter2(nw,wl,yieldno,n,x2,y2,ierr)
3235      IF (ierr .NE. 0) THEN
3236        WRITE(*,*) ierr, fil
3237        STOP
3238      END IF
[2780]3239
[2795]3240      end if !is_master
[2780]3241
[2795]3242      call bcast(yieldno)
[2780]3243
[2795]3244      end subroutine rdxsno
[2780]3245
3246!==============================================================================
3247
[2795]3248      subroutine rdxsn2(nw, wl, yg, yieldn2)
[2780]3249
3250!-----------------------------------------------------------------------------*
3251!=  PURPOSE:                                                                 =*
3252!=  Read n2 cross-sections and photodissociation yield                       =*
3253!-----------------------------------------------------------------------------*
3254!=  PARAMETERS:                                                              =*
3255!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
3256!=           wavelength grid                                                 =*
3257!-----------------------------------------------------------------------------*
3258
[2795]3259      USE mod_phys_lmdz_para, ONLY: is_master
3260      USE mod_phys_lmdz_transfert_para, ONLY: bcast
[2780]3261
[2795]3262      implicit none
[2780]3263
3264!     input
3265
[2795]3266      integer :: nw               ! number of wavelength grid points
3267      real, dimension(nw) :: wl   ! lower wavelength for each interval
[2780]3268
3269!     output
3270
[2795]3271      real, dimension(nw) :: yg        ! n2 cross-sections (cm2)
3272      real, dimension(nw) :: yieldn2   ! n2 photodissociation yield
[2780]3273
3274!     local
3275
[2795]3276      integer, parameter :: kdata = 1100
3277      real, parameter :: deltax = 1.e-4
3278      real, dimension(kdata) :: x1, y1, x2, y2
3279      real :: xl, xu
3280      integer :: i, iw, n, ierr
3281      integer :: kin, kout ! input/output logical units
3282      character*100 fil
[2780]3283
[2795]3284      kin = 10
[2780]3285
3286!     n2 cross sections
3287
[2795]3288      fil = 'cross_sections/n2secef_01nm.txt'
3289      print*, 'section efficace N2: ', fil
[2780]3290
[2795]3291      if (is_master) then
[2780]3292
[2795]3293      OPEN(kin,FILE=fil,STATUS='OLD')
[2780]3294
[2795]3295      n = 1020
3296      DO i = 1, n
3297        READ(kin,*) x1(i), y1(i)
3298      END DO
3299      CLOSE(kin)
[2780]3300
[2795]3301      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
3302      CALL addpnt(x1,y1,kdata,n,          0.,0.)
3303      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
3304      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
[2780]3305
[2795]3306      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
[2780]3307 
[2795]3308      IF (ierr .NE. 0) THEN
3309        WRITE(*,*) ierr, fil
3310        STOP
3311      END IF
[2780]3312
[2795]3313      end if !is_master
[2780]3314
[2795]3315      call bcast(yg)
[2780]3316 
3317!     photodissociation yield
3318
[2795]3319      fil = 'cross_sections/n2_ionef_schunknagy2000.txt'
[2780]3320
[2795]3321      if (is_master) then
[2780]3322
[2795]3323      OPEN(UNIT=kin,FILE=fil,STATUS='old')
[2780]3324
[2795]3325      n = 19
3326      read(kin,*)
3327      DO i = 1, n
3328         READ(kin,*) xl, xu, y2(i)
3329         x2(i) = (xl + xu)/2.
3330         y2(i) = 1. - y2(i)
3331      END DO
3332      CLOSE (kin)
[2780]3333
[2795]3334      CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
3335      CALL addpnt(x2,y2,kdata,n,          0.,0.)
3336      CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
3337      CALL addpnt(x2,y2,kdata,n,      1.e+38,1.)
3338      CALL inter2(nw,wl,yieldn2,n,x2,y2,ierr)
3339      IF (ierr .NE. 0) THEN
3340        WRITE(*,*) ierr, fil
3341        STOP
3342      END IF
[2780]3343
[2795]3344      end if !is_master
[2780]3345
[2795]3346      call bcast(yieldn2)
[2780]3347
[2795]3348      end subroutine rdxsn2
[2780]3349
3350!==============================================================================
3351 
3352      subroutine setalb(nw,wl,albedo)
3353 
3354!-----------------------------------------------------------------------------*
3355!=  PURPOSE:                                                                 =*
3356!=  Set the albedo of the surface.  The albedo is assumed to be Lambertian,  =*
3357!=  i.e., the reflected light is isotropic, and idependt of the direction    =*
3358!=  of incidence of light.  Albedo can be chosen to be wavelength dependent. =*
3359!-----------------------------------------------------------------------------*
3360!=  PARAMETERS:                                                              =*
3361!=  NW      - INTEGER, number of specified intervals + 1 in working       (I)=*
3362!=            wavelength grid                                                =*
3363!=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
3364!=           working wavelength grid                                         =*
3365!=  ALBEDO  - REAL, surface albedo at each specified wavelength           (O)=*
3366!-----------------------------------------------------------------------------*
3367
3368      implicit none
3369
3370! input: (wavelength working grid data)
3371
3372      INTEGER nw
3373      REAL wl(nw)
3374
3375! output:
3376
3377      REAL albedo(nw)
3378
3379! local:
3380
3381      INTEGER iw
3382      REAL alb
3383
3384!     0.015: mean value from clancy et al., icarus, 49-63, 1999.
3385
3386      alb = 0.015
3387
3388      do iw = 1, nw - 1
3389         albedo(iw) = alb
3390      end do
3391
3392      end subroutine setalb
3393 
3394end module photolysis_mod
Note: See TracBrowser for help on using the repository browser.