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
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
19
[705]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',
[1918]42c    $file=trim(datadir)//'/EUVDAT/coln.dat',iostat=ierr)
43     $file=trim(datadir)//'/EUVDAT/param_v6/coln.dat',iostat=ierr)
[705]44
45      IF (ierr.NE.0) THEN
46       write(*,*)'cant find directory EUVDAT containing param_v6 subdir'
47       write(*,*)'(in aeronomars/param_read.F)'
[1918]48       write(*,*)'It should be in :', trim(datadir),'/'
[705]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:'
[1381]53       write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
[705]54       STOP
55      ENDIF
56 
57      !Tabulated photoabsorption coefficients
[1918]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')
[705]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
[1918]236      open(620,file=trim(datadir)//'/EUVDAT/param_v6/fit_js_e107.dat')
[705]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"
[1278]264      else if(solvaryear.eq.31) then
265         filename="e107_MY31.dat"
[1500]266      else if(solvaryear.eq.32) then
267         filename="e107_MY32.dat"
[1717]268      else if(solvaryear.eq.33) then
269         filename="e107_MY33.dat"
[705]270      else
271         write(*,*)"param_read_e107: "
272         write(*,*)"bad value for solvaryear in callphys.def"
[1717]273         write(*,*)"solvaryear must be between 24 and 33"
[705]274         stop
275      endif
276     
[1918]277      open(640,file=trim(datadir)//'/EUVDAT/param_v6/'//filename)
[705]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.
[1888]303         efiono2(inter,1)=0.
304         efiono2(inter,2)=0.
[705]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
[1918]319!      open(120,file=trim(datadir)//'/EUVDAT/param_v5/efdis_inter.dat')
[1888]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)
[705]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
[1918]378      open(130,file=trim(datadir)//'/EUVDAT'//
[1888]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)
[705]389      enddo
[1888]390      close(130)
391      do inter=17,36
392         efdisco2(inter)=1.
[705]393      enddo
[1888]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
[1918]411      open(131,file=trim(datadir)//'/EUVDAT'//
[1888]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)
[705]419      enddo
[1888]420      close(131)
421      do inter=24,36
422         efdiso2(inter)=1.
423      enddo
[705]424
[1888]425
[705]426      !For O(3p) total ionization under 91.1 nm
427      do inter=1,16
428         efiono3p(inter)=1.d0
429      enddo
430
[1888]431
[705]432      !2 channels for N2 ionization
[1918]433      open(132,file=trim(datadir)//'/EUVDAT'//
[1888]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)
[705]441      enddo
[1888]442      close(132)
443      do inter=16,36
444         efdisn2(inter)=1.
[705]445      enddo
[1888]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     
[705]455      !3 channels for CO ionization
[1918]456       open(133,file=trim(datadir)//'/EUVDAT'//
[1888]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)
[705]466      enddo
[1888]467      close(133)
468      do inter=17,36
469         efdisco(inter)=1.
[705]470      enddo
471
[1888]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
[705]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.