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

Last change on this file since 3807 was 3726, checked in by emillour, 2 months ago

Mars PCM:
Turn "callkeys.h" into module "callkeys_mod.F90"
EM

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