source: trunk/LMDZ.MARS/libf/aeronomars/param_read_e107.F @ 2673

Last change on this file since 2673 was 2610, checked in by romain.vande, 3 years ago

LMDZ_MARS RV : Open_MP; Reading files in parallel for thermosphere parametrisation
(callthermos = .true.)

File size: 15.1 KB
Line 
1      subroutine param_read_e107
2
3      use param_v4_h, only: jfotsout,crscabsi2,
4     .    c1_16,c17_24,c25_29,c30_31,c32,c33,c34,c35,c36,
5     .    co2crsc195,co2crsc295,t0,
6     .    jabsifotsintpar,ninter,nz2,
7     .    nabs,e107,date_e107,e107_tab,
8     .    coefit0,coefit1,coefit2,coefit3,coefit4,
9     .    efdisco2,efdiso2,efdish2o,
10     .    efdish2o2,efdish2,efdiso3,
11     .    efdiso,efdisn,efdish,
12     .    efdisno,efdisn2,efdisno2,
13     .    efdisco,efionco2,efiono2,efionn2,
14     .    efionco,efiono3p,efionn,
15     .    efionno,efionh,
16     .    fluxtop,ct1,ct2,p1,p2
17
18      use datafile_mod, only: datadir
19      USE mod_phys_lmdz_para, ONLY: is_master
20      USE mod_phys_lmdz_transfert_para, ONLY: bcast
21
22      implicit none
23
24
25c     common variables and constants
26      include "callkeys.h"
27 
28 
29c     local variables
30
31      integer    i,j,k,inter                          !indexes
32      integer ierr,lnblnk
33      real       nada
34      character*13  filename
35 
36     
37c*************************+PROGRAM STARTS**************************
38
39      if(is_master) then
40
41c     Reads tabulated functions
42
43      !Tabulated column amount
44      open(210, status = 'old',
45c    $file=trim(datadir)//'/EUVDAT/coln.dat',iostat=ierr)
46     $file=trim(datadir)//'/EUVDAT/param_v6/coln.dat',iostat=ierr)
47
48      IF (ierr.NE.0) THEN
49       write(*,*)'cant find directory EUVDAT containing param_v6 subdir'
50       write(*,*)'(in aeronomars/param_read.F)'
51       write(*,*)'It should be in :', trim(datadir),'/'
52       write(*,*)'1) You can change this directory address in '
53       write(*,*)'   callphys.def with datadir=/path/to/dir'
54       write(*,*)'2) If necessary, EUVDAT (and other datafiles)'
55       write(*,*)'   can be obtained online on:'
56       write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
57       STOP
58      ENDIF
59 
60      !Tabulated photoabsorption coefficients
61      open(220,file=trim(datadir)//'/EUVDAT/param_v6/trans2_an.dat')
62      open(230,file=trim(datadir)//'/EUVDAT/param_v6/trans3_an.dat')
63      open(240,file=trim(datadir)//'/EUVDAT/param_v6/trans1_an.dat')
64      open(250,file=trim(datadir)//'/EUVDAT/param_v6/trans2_bn.dat')
65      open(260,file=trim(datadir)//'/EUVDAT/param_v6/trans2_cn.dat')
66      open(300,file=trim(datadir)//'/EUVDAT/param_v6/trans2_dn.dat')
67      open(270,file=trim(datadir)//'/EUVDAT/param_v6/trans1_bn.dat')
68      open(280,file=trim(datadir)//'/EUVDAT/param_v6/trans1_cn.dat')
69      open(290,file=trim(datadir)//'/EUVDAT/param_v6/trans1_dn.dat')
70      open(150,file=trim(datadir)//'/EUVDAT/param_v6/trans4n.dat')
71      open(160,file=trim(datadir)//'/EUVDAT/param_v6/trans5n.dat')
72      open(170,file=trim(datadir)//'/EUVDAT/param_v6/trans6n.dat')
73      open(180,file=trim(datadir)//'/EUVDAT/param_v6/trans7n.dat')
74      open(390,file=trim(datadir)//'/EUVDAT/param_v6/trans8_an.dat')
75      open(400,file=trim(datadir)//'/EUVDAT/param_v6/trans8_bn.dat')
76      open(410,file=trim(datadir)//'/EUVDAT/param_v6/trans9n.dat')
77      open(420,file=trim(datadir)//'/EUVDAT/param_v6/trans10_an.dat')
78      open(430,file=trim(datadir)//'/EUVDAT/param_v6/trans10_bn.dat')
79      open(440,file=trim(datadir)//'/EUVDAT/param_v6/trans10_cn.dat')
80      open(450,file=trim(datadir)//'/EUVDAT/param_v6/trans11_an.dat')
81      open(460,file=trim(datadir)//'/EUVDAT/param_v6/trans11_bn.dat')
82      open(470,file=trim(datadir)//'/EUVDAT/param_v6/trans11_cn.dat')
83      open(480,file=trim(datadir)//'/EUVDAT/param_v6/trans12n.dat')
84      open(490,file=trim(datadir)//'/EUVDAT/param_v6/trans13_an.dat')
85      open(500,file=trim(datadir)//'/EUVDAT/param_v6/trans13_bn.dat')
86      open(510,file=trim(datadir)//'/EUVDAT/param_v6/trans13_cn.dat')
87
88     
89      do i=210,300,10
90         read(i,*)
91         read(i,*)
92      end do
93
94      do i=150,180,10
95         read(i,*)
96         read(i,*)
97      end do
98
99      do i=390,510,10
100         read(i,*)
101         read(i,*)
102      enddo
103
104      do i=nz2,1,-1
105         read(210,*) (c1_16(i,j),j=1,16),c17_24(i),c25_29(i),c30_31(i),
106     $        c32(i),c33(i),c34(i),c35(i),c36(i)
107      end do
108
109      do i=nz2,1,-1
110         read(220,*) (jabsifotsintpar(i,2,j),j=1,16)
111      end do
112     
113      do i=nz2,1,-1
114         read(230,*) (jabsifotsintpar(i,3,j),j=1,16)
115      end do
116
117      do i=nz2,1,-1
118         read(240,*) (jabsifotsintpar(i,1,j),j=1,16)
119      end do
120
121      do i=nz2,1,-1
122         read(250,*) (jabsifotsintpar(i,2,j),j=17,24)
123      end do
124
125
126      do i=nz2,1,-1
127         read(260,*) (jabsifotsintpar(i,2,j),j=25,31)
128      end do
129
130      do i=nz2,1,-1
131         read(270,*) (jabsifotsintpar(i,1,j),j=17,24)
132      end do
133
134      do i=nz2,1,-1
135         read(280,*) (jabsifotsintpar(i,1,j),j=25,31)
136      end do
137
138      do i=nz2,1,-1
139         read(290,*) jabsifotsintpar(i,1,32)
140      end do
141
142      do i=nz2,1,-1
143         read(300,*) (jabsifotsintpar(i,2,j),j=32,34)
144      end do
145
146      do i=nz2,1,-1
147         read(160,*) (jabsifotsintpar(i,5,j),j=1,15)
148      end do
149
150      do i=nz2,1,-1
151         read(150,*) (jabsifotsintpar(i,4,j),j=25,31)
152      end do
153
154      do i=nz2,1,-1
155         read(170,*) (jabsifotsintpar(i,6,j),j=25,35)
156      end do
157
158      do i=nz2,1,-1
159         read(180,*) (jabsifotsintpar(i,7,j),j=34,36)
160      end do
161
162      do i=nz2,1,-1
163         read(390,*) (jabsifotsintpar(i,8,j),j=2,16)
164      enddo
165
166      do i=nz2,1,-1
167         read(400,*) (jabsifotsintpar(i,8,j),j=17,24)
168      enddo
169
170      do i=nz2,1,-1
171         read(410,*) (jabsifotsintpar(i,9,j),j=1,16)
172      enddo
173
174      do i=nz2,1,-1
175         read(420,*) (jabsifotsintpar(i,10,j),j=2,16)
176      enddo
177
178      do i=nz2,1,-1
179         read(430,*) (jabsifotsintpar(i,10,j),j=17,24)
180      enddo
181
182      do i=nz2,1,-1
183         read(440,*) (jabsifotsintpar(i,10,j),j=25,32)
184      enddo
185
186      do i=nz2,1,-1
187         read(450,*) (jabsifotsintpar(i,11,j),j=2,16)
188      enddo
189
190      do i=nz2,1,-1
191         read(460,*) (jabsifotsintpar(i,11,j),j=17,24)
192      enddo
193
194      do i=nz2,1,-1
195         read(470,*) (jabsifotsintpar(i,11,j),j=25,29)
196      enddo
197     
198      do i=nz2,1,-1
199         read(480,*) (jabsifotsintpar(i,12,j),j=2,16)
200      enddo
201
202      do i=nz2,1,-1
203         read(490,*) (jabsifotsintpar(i,13,j),j=2,16)
204      enddo
205     
206      do i=nz2,1,-1
207         read(500,*) (jabsifotsintpar(i,13,j),j=17,24)
208      enddo
209     
210      do i=nz2,1,-1
211         read(510,*) (jabsifotsintpar(i,13,j),j=25,36)
212      enddo
213
214      do i=210,300,10
215         close(i)
216      end do
217
218      do i=150,180,10
219         close(i)
220      end do
221
222      do i=390,510,10
223         close(i)
224      enddo
225
226      endif !is master
227
228c     set t0
229
230      do i=1,nz2
231         t0(i)=195.
232      end do
233
234
235      do i=1,ninter
236         fluxtop(i)=1.
237      end do
238
239      if(is_master) then
240
241      !Parameters for the variation of the solar flux with 11 years cycle
242      open(620,file=trim(datadir)//'/EUVDAT/param_v6/fit_js_e107.dat')
243      do i=1,ninter
244         read(620,*)
245         do j=1,nabs
246            read(620,*)coefit0(i,j),coefit1(i,j),coefit2(i,j),
247     $           coefit3(i,j),coefit4(i,j)
248         enddo
249      enddo
250      close(620)
251
252
253      !Tabulated values for E10.7
254      if(solvaryear.eq.23) then
255         filename="e107_MY23.dat"
256      else if(solvaryear.eq.24) then
257         filename="e107_MY24.dat"
258      else if(solvaryear.eq.25) then
259         filename="e107_MY25.dat"
260      else if (solvaryear.eq.26) then
261         filename="e107_MY26.dat"
262      else if (solvaryear.eq.27) then
263         filename="e107_MY27.dat"
264      else if (solvaryear.eq.28) then
265         filename="e107_MY28.dat"
266      else if (solvaryear.eq.29) then
267         filename="e107_MY29.dat"
268      else if(solvaryear.eq.30) then
269         filename="e107_MY30.dat"
270      else if(solvaryear.eq.31) then
271         filename="e107_MY31.dat"
272      else if(solvaryear.eq.32) then
273         filename="e107_MY32.dat"
274      else if(solvaryear.eq.33) then
275         filename="e107_MY33.dat"
276      else if(solvaryear.eq.34) then
277         filename="e107_MY34.dat"
278      else
279         write(*,*)"param_read_e107: "
280         write(*,*)"bad value for solvaryear in callphys.def"
281         write(*,*)"solvaryear must be between 24 and 33"
282         stop
283      endif
284     
285      open(640,file=trim(datadir)//'/EUVDAT/param_v6/'//filename)
286      read(640,*)
287      do i=1,669
288         read(640,*)date_e107(i),e107_tab(i)
289         write(*,*)'param_read_e107/292',date_e107(i),e107_tab(i),i
290      enddo
291      close(640)
292
293      endif !is master
294
295      call bcast(t0)
296      call bcast(fluxtop)
297
298      call bcast(c1_16)
299      call bcast(c17_24)
300      call bcast(c25_29)
301      call bcast(c30_31)
302      call bcast(c32)
303      call bcast(c33)
304      call bcast(c34)
305      call bcast(c35)
306      call bcast(c36)
307
308      call bcast(jabsifotsintpar)
309
310      call bcast(coefit0)
311      call bcast(coefit1)
312      call bcast(coefit2)
313      call bcast(coefit3)
314      call bcast(coefit4)
315      call bcast(date_e107)
316      call bcast(e107_tab)
317         
318
319c     dissociation and ionization efficiencies
320
321      do inter=1,ninter
322         efdisco2(inter)=0.
323         efdiso2(inter)=0.
324         efdish2(inter)=0.
325         efdish2o(inter)=0.
326         efdish2o2(inter)=0.
327         efdiso3(inter)=0.
328         efdisco(inter)=0.
329         efdisn2(inter)=0.
330         efdisno(inter)=0.
331         efdisno2(inter)=0.
332         efionco2(inter,1)=0.
333         efionco2(inter,2)=0.
334         efionco2(inter,3)=0.
335         efionco2(inter,4)=0.
336         efiono2(inter,1)=0.
337         efiono2(inter,2)=0.
338         efiono3p(inter)=0.
339         efionn2(inter,1)=0.
340         efionn2(inter,2)=0.
341         efionco(inter,1)=0.
342         efionco(inter,2)=0.
343         efionco(inter,3)=0.
344         efionn(inter)=0.
345         efionh(inter)=0.
346         efionno(inter)=0.
347      enddo
348
349
350c     CO2, O2, NO
351
352!      open(120,file=trim(datadir)//'/EUVDAT/param_v5/efdis_inter.dat')
353!      read(120,*)
354!!      do i=1,21
355!!         read(120,*)inter,efdisco2(inter),efdiso2(inter),efdisno(inter)
356!      do inter=8,28
357!         read(120,*)i,efdisco2(inter),efdiso2(inter),efdisno(inter)
358!      enddo
359!      do inter=29,ninter
360!         efdisco2(inter)=1.
361!         efdiso2(inter)=1.
362!         efdisno(inter)=1.
363!      enddo
364!      close(120)
365
366
367c     N2
368
369      efdisn2(15)=0.1
370      do inter=16,ninter
371         efdisn2(inter)=1.
372      enddo
373
374
375c     CO
376
377      efdisco(16)=0.5
378      do inter=17,ninter
379         efdisco(inter)=1.
380      enddo
381
382     
383c     O, N, H
384
385      do inter=1,ninter
386         efdiso(inter)=0.
387         efdisn(inter)=0.
388         efdish(inter)=0.
389      enddo
390
391
392c     H2O, H2O2, O3, NO2
393
394      do inter=25,31
395         efdish2o(inter)=1.
396      enddo
397      do inter=25,35
398         efdish2o2(inter)=1.
399      enddo
400      do inter=34,36
401         efdiso3(inter)=1.
402      enddo
403      do inter=27,36
404         efdisno2(inter)=1.
405      enddo
406      do inter=1,15
407         efdish2(inter)=1.
408      enddo
409
410      if(is_master) then
411         
412      !4 possible channels for CO2 ionization
413      open(130,file=trim(datadir)//'/EUVDAT'//
414     $     '/co2ion_branchingratio_schunkandnagy2000_param.dat')
415      do inter=1,16
416         read(130,*)i,nada,efionco2(inter,1),efionco2(inter,2),
417     $        efionco2(inter,3),efionco2(inter,4)
418         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
419         efdisco2(inter)=1.-nada
420         efionco2(inter,1)=(1.-efdisco2(inter))*efionco2(inter,1)
421         efionco2(inter,2)=(1.-efdisco2(inter))*efionco2(inter,2)
422         efionco2(inter,3)=(1.-efdisco2(inter))*efionco2(inter,3)
423         efionco2(inter,4)=(1.-efdisco2(inter))*efionco2(inter,4)
424      enddo
425      close(130)
426
427      endif !is_master
428
429      call bcast(i)
430      call bcast(nada)
431      call bcast(efionco2)
432      call bcast(efdisco2)
433
434      do inter=17,36
435         efdisco2(inter)=1.
436      enddo
437
438!      do inter=14,16
439!         efionco2(inter,1)=1.-efdisco2(inter)
440!      enddo
441!      efionco2(13,1)=0.805*(1.-efdisco2(13))
442!      efionco2(13,2)=0.195*(1.-efdisco2(13))
443!      do inter=11,12
444!         efionco2(inter,3)=1.-efdisco2(inter)
445!      enddo
446!      efionco2(10,3)=0.9*(1.-efdisco2(10))
447!      efionco2(10,4)=0.1*(1.-efdisco2(10))
448!      do inter=2,9
449!         efionco2(inter,4)=1.-efdisco2(inter)
450!      enddo
451
452      if(is_master) then
453      !2 possible channels for O2 ionization
454      open(131,file=trim(datadir)//'/EUVDAT'//
455     $     '/o2ion_branchingratio_schunkandnagy2000_param.dat')
456      do inter=1,23
457         read(131,*)i,nada,efiono2(inter,1),efiono2(inter,2)
458         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
459         efdiso2(inter)=1.-nada
460         efiono2(inter,1)=(1.-efdiso2(inter))*efiono2(inter,1)
461         efiono2(inter,2)=(1.-efdiso2(inter))*efiono2(inter,2)
462      enddo
463      close(131)
464
465      endif !is_master
466
467      call bcast(i)
468      call bcast(nada)
469      call bcast(efiono2)
470      call bcast(efdiso2)
471
472      do inter=24,36
473         efdiso2(inter)=1.
474      enddo
475
476
477      !For O(3p) total ionization under 91.1 nm
478      do inter=1,16
479         efiono3p(inter)=1.d0
480      enddo
481
482      if(is_master) then
483
484      !2 channels for N2 ionization
485      open(132,file=trim(datadir)//'/EUVDAT'//
486     $     '/n2ion_branchingratio_schunkandnagy2000_param.dat')
487      do inter=1,15
488         read(132,*)i,nada,efionn2(inter,1),efionn2(inter,2)
489         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
490         efdisn2(inter)=1.-nada
491         efionn2(inter,1)=(1.-efdisn2(inter))*efionn2(inter,1)
492         efionn2(inter,2)=(1.-efdisn2(inter))*efionn2(inter,2)
493      enddo
494      close(132)
495
496      endif
497
498      call bcast(i)
499      call bcast(nada)
500      call bcast(efionn2)
501      call bcast(efdisn2)
502
503      do inter=16,36
504         efdisn2(inter)=1.
505      enddo
506
507!      do inter=9,15
508!         efionn2(inter,1)=1.-efdisn2(inter)
509!      enddo
510!      do inter=2,8
511!         efionn2(inter,2)=1.-efdisn2(inter)
512!      enddo
513 
514      if(is_master) then     
515      !3 channels for CO ionization
516       open(133,file=trim(datadir)//'/EUVDAT'//
517     $     '/coion_branchingratio_schunkandnagy2000_param.dat')
518      do inter=1,16
519         read(133,*)i,nada,efionco(inter,1),efionco(inter,2),
520     $        efionco(inter,3)
521         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
522         efdisco(inter)=1.-nada
523         efionco(inter,1)=(1.-efdisco(inter))*efionco(inter,1)
524         efionco(inter,2)=(1.-efdisco(inter))*efionco(inter,2)
525         efionco(inter,3)=(1.-efdisco(inter))*efionco(inter,3)
526      enddo
527      close(133)
528
529      endif  !     if(is_master)
530
531      call bcast(i)
532      call bcast(nada)
533      call bcast(efionco)
534      call bcast(efdisco)
535
536      do inter=17,36
537         efdisco(inter)=1.
538      enddo
539
540!      do inter=11,16
541!         efionco(inter,1)=1.-efdisco(inter)
542!      enddo
543!      efionco(10,1)=0.87*(1.-efdisco(10))
544!      efionco(10,2)=0.13*(1.-efdisco(10))
545!      do inter=8,9
546!         efionco(inter,2)=1.-efdisco(inter)
547!      enddo
548!      efionco(7,2)=0.1*(1.-efdisco(7))
549!      efionco(7,3)=0.9*(1.-efdisco(7))
550!      do inter=2,6
551!         efionco(inter,3)=1.-efdisco(inter)
552!      enddo
553
554
555      !Total ionization under 85 nm for N
556      do inter=1,16
557         efionn(inter)=1.
558      enddo
559
560      !NO
561      do inter=2,28
562         efionno(inter)=1.-efdisno(inter)
563      enddo
564
565      !Total ionization under 90 nm for H
566      do inter=3,16
567         efionh(inter)=1.
568      enddo
569
570
571      return
572
573
574      end
575
Note: See TracBrowser for help on using the repository browser.