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

Last change on this file since 2042 was 1918, checked in by emillour, 7 years ago

Mars GCM:
Code cleanup:

  • remove "comorbit.h" since it is no longer used.
  • turn "datafile.h" into module datafile_mod.F90 (and rename variable "datafile" as "datadir" since it stores the path to the datafile directory).

EM

File size: 13.8 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
271         write(*,*)"param_read_e107: "
272         write(*,*)"bad value for solvaryear in callphys.def"
273         write(*,*)"solvaryear must be between 24 and 33"
274         stop
275      endif
276     
277      open(640,file=trim(datadir)//'/EUVDAT/param_v6/'//filename)
278      read(640,*)
279      do i=1,669
280         read(640,*)date_e107(i),e107_tab(i)
281         write(*,*)'param_read_e107/292',date_e107(i),e107_tab(i),i
282      enddo
283      close(640)
284         
285
286c     dissociation and ionization efficiencies
287
288      do inter=1,ninter
289         efdisco2(inter)=0.
290         efdiso2(inter)=0.
291         efdish2(inter)=0.
292         efdish2o(inter)=0.
293         efdish2o2(inter)=0.
294         efdiso3(inter)=0.
295         efdisco(inter)=0.
296         efdisn2(inter)=0.
297         efdisno(inter)=0.
298         efdisno2(inter)=0.
299         efionco2(inter,1)=0.
300         efionco2(inter,2)=0.
301         efionco2(inter,3)=0.
302         efionco2(inter,4)=0.
303         efiono2(inter,1)=0.
304         efiono2(inter,2)=0.
305         efiono3p(inter)=0.
306         efionn2(inter,1)=0.
307         efionn2(inter,2)=0.
308         efionco(inter,1)=0.
309         efionco(inter,2)=0.
310         efionco(inter,3)=0.
311         efionn(inter)=0.
312         efionh(inter)=0.
313         efionno(inter)=0.
314      enddo
315
316
317c     CO2, O2, NO
318
319!      open(120,file=trim(datadir)//'/EUVDAT/param_v5/efdis_inter.dat')
320!      read(120,*)
321!!      do i=1,21
322!!         read(120,*)inter,efdisco2(inter),efdiso2(inter),efdisno(inter)
323!      do inter=8,28
324!         read(120,*)i,efdisco2(inter),efdiso2(inter),efdisno(inter)
325!      enddo
326!      do inter=29,ninter
327!         efdisco2(inter)=1.
328!         efdiso2(inter)=1.
329!         efdisno(inter)=1.
330!      enddo
331!      close(120)
332
333
334c     N2
335
336      efdisn2(15)=0.1
337      do inter=16,ninter
338         efdisn2(inter)=1.
339      enddo
340
341
342c     CO
343
344      efdisco(16)=0.5
345      do inter=17,ninter
346         efdisco(inter)=1.
347      enddo
348
349     
350c     O, N, H
351
352      do inter=1,ninter
353         efdiso(inter)=0.
354         efdisn(inter)=0.
355         efdish(inter)=0.
356      enddo
357
358
359c     H2O, H2O2, O3, NO2
360
361      do inter=25,31
362         efdish2o(inter)=1.
363      enddo
364      do inter=25,35
365         efdish2o2(inter)=1.
366      enddo
367      do inter=34,36
368         efdiso3(inter)=1.
369      enddo
370      do inter=27,36
371         efdisno2(inter)=1.
372      enddo
373      do inter=1,15
374         efdish2(inter)=1.
375      enddo
376         
377      !4 possible channels for CO2 ionization
378      open(130,file=trim(datadir)//'/EUVDAT'//
379     $     '/co2ion_branchingratio_schunkandnagy2000_param.dat')
380      do inter=1,16
381         read(130,*)i,nada,efionco2(inter,1),efionco2(inter,2),
382     $        efionco2(inter,3),efionco2(inter,4)
383         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
384         efdisco2(inter)=1.-nada
385         efionco2(inter,1)=(1.-efdisco2(inter))*efionco2(inter,1)
386         efionco2(inter,2)=(1.-efdisco2(inter))*efionco2(inter,2)
387         efionco2(inter,3)=(1.-efdisco2(inter))*efionco2(inter,3)
388         efionco2(inter,4)=(1.-efdisco2(inter))*efionco2(inter,4)
389      enddo
390      close(130)
391      do inter=17,36
392         efdisco2(inter)=1.
393      enddo
394
395!      do inter=14,16
396!         efionco2(inter,1)=1.-efdisco2(inter)
397!      enddo
398!      efionco2(13,1)=0.805*(1.-efdisco2(13))
399!      efionco2(13,2)=0.195*(1.-efdisco2(13))
400!      do inter=11,12
401!         efionco2(inter,3)=1.-efdisco2(inter)
402!      enddo
403!      efionco2(10,3)=0.9*(1.-efdisco2(10))
404!      efionco2(10,4)=0.1*(1.-efdisco2(10))
405!      do inter=2,9
406!         efionco2(inter,4)=1.-efdisco2(inter)
407!      enddo
408
409
410      !2 possible channels for O2 ionization
411      open(131,file=trim(datadir)//'/EUVDAT'//
412     $     '/o2ion_branchingratio_schunkandnagy2000_param.dat')
413      do inter=1,23
414         read(131,*)i,nada,efiono2(inter,1),efiono2(inter,2)
415         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
416         efdiso2(inter)=1.-nada
417         efiono2(inter,1)=(1.-efdiso2(inter))*efiono2(inter,1)
418         efiono2(inter,2)=(1.-efdiso2(inter))*efiono2(inter,2)
419      enddo
420      close(131)
421      do inter=24,36
422         efdiso2(inter)=1.
423      enddo
424
425
426      !For O(3p) total ionization under 91.1 nm
427      do inter=1,16
428         efiono3p(inter)=1.d0
429      enddo
430
431
432      !2 channels for N2 ionization
433      open(132,file=trim(datadir)//'/EUVDAT'//
434     $     '/n2ion_branchingratio_schunkandnagy2000_param.dat')
435      do inter=1,15
436         read(132,*)i,nada,efionn2(inter,1),efionn2(inter,2)
437         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
438         efdisn2(inter)=1.-nada
439         efionn2(inter,1)=(1.-efdisn2(inter))*efionn2(inter,1)
440         efionn2(inter,2)=(1.-efdisn2(inter))*efionn2(inter,2)
441      enddo
442      close(132)
443      do inter=16,36
444         efdisn2(inter)=1.
445      enddo
446
447!      do inter=9,15
448!         efionn2(inter,1)=1.-efdisn2(inter)
449!      enddo
450!      do inter=2,8
451!         efionn2(inter,2)=1.-efdisn2(inter)
452!      enddo
453 
454     
455      !3 channels for CO ionization
456       open(133,file=trim(datadir)//'/EUVDAT'//
457     $     '/coion_branchingratio_schunkandnagy2000_param.dat')
458      do inter=1,16
459         read(133,*)i,nada,efionco(inter,1),efionco(inter,2),
460     $        efionco(inter,3)
461         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
462         efdisco(inter)=1.-nada
463         efionco(inter,1)=(1.-efdisco(inter))*efionco(inter,1)
464         efionco(inter,2)=(1.-efdisco(inter))*efionco(inter,2)
465         efionco(inter,3)=(1.-efdisco(inter))*efionco(inter,3)
466      enddo
467      close(133)
468      do inter=17,36
469         efdisco(inter)=1.
470      enddo
471
472!      do inter=11,16
473!         efionco(inter,1)=1.-efdisco(inter)
474!      enddo
475!      efionco(10,1)=0.87*(1.-efdisco(10))
476!      efionco(10,2)=0.13*(1.-efdisco(10))
477!      do inter=8,9
478!         efionco(inter,2)=1.-efdisco(inter)
479!      enddo
480!      efionco(7,2)=0.1*(1.-efdisco(7))
481!      efionco(7,3)=0.9*(1.-efdisco(7))
482!      do inter=2,6
483!         efionco(inter,3)=1.-efdisco(inter)
484!      enddo
485
486
487      !Total ionization under 85 nm for N
488      do inter=1,16
489         efionn(inter)=1.
490      enddo
491
492      !NO
493      do inter=2,28
494         efionno(inter)=1.-efdisno(inter)
495      enddo
496
497      !Total ionization under 90 nm for H
498      do inter=3,16
499         efionh(inter)=1.
500      enddo
501
502
503      return
504
505
506      end
507
Note: See TracBrowser for help on using the repository browser.