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

Last change on this file since 3026 was 2921, checked in by emillour, 22 months ago

Mars PCM:
Add the possibility to use MY36 dust scenario and MY36 EUV inputs
EM
File size: 14.4 KB
RevLine 
[705]1      subroutine param_read_e107
2
[1266]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,
[1888]13     .    efdisco,efionco2,efiono2,efionn2,
[1266]14     .    efionco,efiono3p,efionn,
15     .    efionno,efionh,
16     .    fluxtop,ct1,ct2,p1,p2
17
[1918]18      use datafile_mod, only: datadir
[2610]19      USE mod_phys_lmdz_para, ONLY: is_master
20      USE mod_phys_lmdz_transfert_para, ONLY: bcast
[1918]21
[705]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
[2610]39      if(is_master) then
[705]40
41c     Reads tabulated functions
42
43      !Tabulated column amount
44      open(210, status = 'old',
[1918]45c    $file=trim(datadir)//'/EUVDAT/coln.dat',iostat=ierr)
46     $file=trim(datadir)//'/EUVDAT/param_v6/coln.dat',iostat=ierr)
[705]47
48      IF (ierr.NE.0) THEN
49       write(*,*)'cant find directory EUVDAT containing param_v6 subdir'
50       write(*,*)'(in aeronomars/param_read.F)'
[1918]51       write(*,*)'It should be in :', trim(datadir),'/'
[705]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:'
[1381]56       write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
[705]57       STOP
58      ENDIF
59 
60      !Tabulated photoabsorption coefficients
[1918]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')
[705]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
[2610]226      endif !is master
[705]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
[2610]239      if(is_master) then
240
[705]241      !Parameters for the variation of the solar flux with 11 years cycle
[1918]242      open(620,file=trim(datadir)//'/EUVDAT/param_v6/fit_js_e107.dat')
[705]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
[2921]254      if((solvaryear.ge.23).and.(solvaryear.le.36)) then
[2689]255        write(filename,fmt="(a,i2,a)")"e107_MY",solvaryear,".dat"
[705]256      else
257         write(*,*)"param_read_e107: "
258         write(*,*)"bad value for solvaryear in callphys.def"
[2921]259         write(*,*)"solvaryear must be between 24 and 36"
[705]260         stop
261      endif
262     
[1918]263      open(640,file=trim(datadir)//'/EUVDAT/param_v6/'//filename)
[705]264      read(640,*)
265      do i=1,669
266         read(640,*)date_e107(i),e107_tab(i)
267         write(*,*)'param_read_e107/292',date_e107(i),e107_tab(i),i
268      enddo
269      close(640)
[2610]270
271      endif !is master
272
273      call bcast(t0)
274      call bcast(fluxtop)
275
276      call bcast(c1_16)
277      call bcast(c17_24)
278      call bcast(c25_29)
279      call bcast(c30_31)
280      call bcast(c32)
281      call bcast(c33)
282      call bcast(c34)
283      call bcast(c35)
284      call bcast(c36)
285
286      call bcast(jabsifotsintpar)
287
288      call bcast(coefit0)
289      call bcast(coefit1)
290      call bcast(coefit2)
291      call bcast(coefit3)
292      call bcast(coefit4)
293      call bcast(date_e107)
294      call bcast(e107_tab)
[705]295         
296
297c     dissociation and ionization efficiencies
298
299      do inter=1,ninter
300         efdisco2(inter)=0.
301         efdiso2(inter)=0.
302         efdish2(inter)=0.
303         efdish2o(inter)=0.
304         efdish2o2(inter)=0.
305         efdiso3(inter)=0.
306         efdisco(inter)=0.
307         efdisn2(inter)=0.
308         efdisno(inter)=0.
309         efdisno2(inter)=0.
310         efionco2(inter,1)=0.
311         efionco2(inter,2)=0.
312         efionco2(inter,3)=0.
313         efionco2(inter,4)=0.
[1888]314         efiono2(inter,1)=0.
315         efiono2(inter,2)=0.
[705]316         efiono3p(inter)=0.
317         efionn2(inter,1)=0.
318         efionn2(inter,2)=0.
319         efionco(inter,1)=0.
320         efionco(inter,2)=0.
321         efionco(inter,3)=0.
322         efionn(inter)=0.
323         efionh(inter)=0.
324         efionno(inter)=0.
325      enddo
326
327
328c     CO2, O2, NO
329
[1918]330!      open(120,file=trim(datadir)//'/EUVDAT/param_v5/efdis_inter.dat')
[1888]331!      read(120,*)
332!!      do i=1,21
333!!         read(120,*)inter,efdisco2(inter),efdiso2(inter),efdisno(inter)
334!      do inter=8,28
335!         read(120,*)i,efdisco2(inter),efdiso2(inter),efdisno(inter)
336!      enddo
337!      do inter=29,ninter
338!         efdisco2(inter)=1.
339!         efdiso2(inter)=1.
340!         efdisno(inter)=1.
341!      enddo
342!      close(120)
[705]343
344
345c     N2
346
347      efdisn2(15)=0.1
348      do inter=16,ninter
349         efdisn2(inter)=1.
350      enddo
351
352
353c     CO
354
355      efdisco(16)=0.5
356      do inter=17,ninter
357         efdisco(inter)=1.
358      enddo
359
360     
361c     O, N, H
362
363      do inter=1,ninter
364         efdiso(inter)=0.
365         efdisn(inter)=0.
366         efdish(inter)=0.
367      enddo
368
369
370c     H2O, H2O2, O3, NO2
371
372      do inter=25,31
373         efdish2o(inter)=1.
374      enddo
375      do inter=25,35
376         efdish2o2(inter)=1.
377      enddo
378      do inter=34,36
379         efdiso3(inter)=1.
380      enddo
381      do inter=27,36
382         efdisno2(inter)=1.
383      enddo
384      do inter=1,15
385         efdish2(inter)=1.
386      enddo
[2610]387
388      if(is_master) then
[705]389         
390      !4 possible channels for CO2 ionization
[1918]391      open(130,file=trim(datadir)//'/EUVDAT'//
[1888]392     $     '/co2ion_branchingratio_schunkandnagy2000_param.dat')
393      do inter=1,16
394         read(130,*)i,nada,efionco2(inter,1),efionco2(inter,2),
395     $        efionco2(inter,3),efionco2(inter,4)
396         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
397         efdisco2(inter)=1.-nada
398         efionco2(inter,1)=(1.-efdisco2(inter))*efionco2(inter,1)
399         efionco2(inter,2)=(1.-efdisco2(inter))*efionco2(inter,2)
400         efionco2(inter,3)=(1.-efdisco2(inter))*efionco2(inter,3)
401         efionco2(inter,4)=(1.-efdisco2(inter))*efionco2(inter,4)
[705]402      enddo
[1888]403      close(130)
[2610]404
405      endif !is_master
406
407      call bcast(i)
408      call bcast(nada)
409      call bcast(efionco2)
410      call bcast(efdisco2)
411
[1888]412      do inter=17,36
413         efdisco2(inter)=1.
[705]414      enddo
[1888]415
416!      do inter=14,16
417!         efionco2(inter,1)=1.-efdisco2(inter)
418!      enddo
419!      efionco2(13,1)=0.805*(1.-efdisco2(13))
420!      efionco2(13,2)=0.195*(1.-efdisco2(13))
421!      do inter=11,12
422!         efionco2(inter,3)=1.-efdisco2(inter)
423!      enddo
424!      efionco2(10,3)=0.9*(1.-efdisco2(10))
425!      efionco2(10,4)=0.1*(1.-efdisco2(10))
426!      do inter=2,9
427!         efionco2(inter,4)=1.-efdisco2(inter)
428!      enddo
429
[2610]430      if(is_master) then
[1888]431      !2 possible channels for O2 ionization
[1918]432      open(131,file=trim(datadir)//'/EUVDAT'//
[1888]433     $     '/o2ion_branchingratio_schunkandnagy2000_param.dat')
434      do inter=1,23
435         read(131,*)i,nada,efiono2(inter,1),efiono2(inter,2)
436         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
437         efdiso2(inter)=1.-nada
438         efiono2(inter,1)=(1.-efdiso2(inter))*efiono2(inter,1)
439         efiono2(inter,2)=(1.-efdiso2(inter))*efiono2(inter,2)
[705]440      enddo
[1888]441      close(131)
[2610]442
443      endif !is_master
444
445      call bcast(i)
446      call bcast(nada)
447      call bcast(efiono2)
448      call bcast(efdiso2)
449
[1888]450      do inter=24,36
451         efdiso2(inter)=1.
452      enddo
[705]453
[1888]454
[705]455      !For O(3p) total ionization under 91.1 nm
456      do inter=1,16
457         efiono3p(inter)=1.d0
458      enddo
459
[2610]460      if(is_master) then
[1888]461
[705]462      !2 channels for N2 ionization
[1918]463      open(132,file=trim(datadir)//'/EUVDAT'//
[1888]464     $     '/n2ion_branchingratio_schunkandnagy2000_param.dat')
465      do inter=1,15
466         read(132,*)i,nada,efionn2(inter,1),efionn2(inter,2)
467         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
468         efdisn2(inter)=1.-nada
469         efionn2(inter,1)=(1.-efdisn2(inter))*efionn2(inter,1)
470         efionn2(inter,2)=(1.-efdisn2(inter))*efionn2(inter,2)
[705]471      enddo
[1888]472      close(132)
[2610]473
474      endif
475
476      call bcast(i)
477      call bcast(nada)
478      call bcast(efionn2)
479      call bcast(efdisn2)
480
[1888]481      do inter=16,36
482         efdisn2(inter)=1.
[705]483      enddo
[1888]484
485!      do inter=9,15
486!         efionn2(inter,1)=1.-efdisn2(inter)
487!      enddo
488!      do inter=2,8
489!         efionn2(inter,2)=1.-efdisn2(inter)
490!      enddo
491 
[2610]492      if(is_master) then     
[705]493      !3 channels for CO ionization
[1918]494       open(133,file=trim(datadir)//'/EUVDAT'//
[1888]495     $     '/coion_branchingratio_schunkandnagy2000_param.dat')
496      do inter=1,16
497         read(133,*)i,nada,efionco(inter,1),efionco(inter,2),
498     $        efionco(inter,3)
499         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
500         efdisco(inter)=1.-nada
501         efionco(inter,1)=(1.-efdisco(inter))*efionco(inter,1)
502         efionco(inter,2)=(1.-efdisco(inter))*efionco(inter,2)
503         efionco(inter,3)=(1.-efdisco(inter))*efionco(inter,3)
[705]504      enddo
[1888]505      close(133)
[2610]506
507      endif  !     if(is_master)
508
509      call bcast(i)
510      call bcast(nada)
511      call bcast(efionco)
512      call bcast(efdisco)
513
[1888]514      do inter=17,36
515         efdisco(inter)=1.
[705]516      enddo
517
[1888]518!      do inter=11,16
519!         efionco(inter,1)=1.-efdisco(inter)
520!      enddo
521!      efionco(10,1)=0.87*(1.-efdisco(10))
522!      efionco(10,2)=0.13*(1.-efdisco(10))
523!      do inter=8,9
524!         efionco(inter,2)=1.-efdisco(inter)
525!      enddo
526!      efionco(7,2)=0.1*(1.-efdisco(7))
527!      efionco(7,3)=0.9*(1.-efdisco(7))
528!      do inter=2,6
529!         efionco(inter,3)=1.-efdisco(inter)
530!      enddo
531
532
[705]533      !Total ionization under 85 nm for N
534      do inter=1,16
535         efionn(inter)=1.
536      enddo
537
538      !NO
539      do inter=2,28
540         efionno(inter)=1.-efdisno(inter)
541      enddo
542
543      !Total ionization under 90 nm for H
544      do inter=3,16
545         efionh(inter)=1.
546      enddo
547
548
549      return
550
551
552      end
553
Note: See TracBrowser for help on using the repository browser.