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

Last change on this file since 1403 was 1381, checked in by emillour, 10 years ago

Mars GCM:

EM

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