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

Last change on this file since 2325 was 2137, checked in by emillour, 6 years ago

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