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

Last change on this file since 2530 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
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"
[2137]270      else if(solvaryear.eq.34) then
271         filename="e107_MY34.dat"
[705]272      else
273         write(*,*)"param_read_e107: "
274         write(*,*)"bad value for solvaryear in callphys.def"
[1717]275         write(*,*)"solvaryear must be between 24 and 33"
[705]276         stop
277      endif
278     
[1918]279      open(640,file=trim(datadir)//'/EUVDAT/param_v6/'//filename)
[705]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.
[1888]305         efiono2(inter,1)=0.
306         efiono2(inter,2)=0.
[705]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
[1918]321!      open(120,file=trim(datadir)//'/EUVDAT/param_v5/efdis_inter.dat')
[1888]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)
[705]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
[1918]380      open(130,file=trim(datadir)//'/EUVDAT'//
[1888]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)
[705]391      enddo
[1888]392      close(130)
393      do inter=17,36
394         efdisco2(inter)=1.
[705]395      enddo
[1888]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
[1918]413      open(131,file=trim(datadir)//'/EUVDAT'//
[1888]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)
[705]421      enddo
[1888]422      close(131)
423      do inter=24,36
424         efdiso2(inter)=1.
425      enddo
[705]426
[1888]427
[705]428      !For O(3p) total ionization under 91.1 nm
429      do inter=1,16
430         efiono3p(inter)=1.d0
431      enddo
432
[1888]433
[705]434      !2 channels for N2 ionization
[1918]435      open(132,file=trim(datadir)//'/EUVDAT'//
[1888]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)
[705]443      enddo
[1888]444      close(132)
445      do inter=16,36
446         efdisn2(inter)=1.
[705]447      enddo
[1888]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     
[705]457      !3 channels for CO ionization
[1918]458       open(133,file=trim(datadir)//'/EUVDAT'//
[1888]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)
[705]468      enddo
[1888]469      close(133)
470      do inter=17,36
471         efdisco(inter)=1.
[705]472      enddo
473
[1888]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
[705]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.