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

Last change on this file since 1771 was 1717, checked in by emillour, 8 years ago

Mars GCM:

  • add possibility to read in MY33 EUV forcing.

FGG

File size: 11.2 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 if(solvaryear.eq.32) then
266         filename="e107_MY32.dat"
267      else if(solvaryear.eq.33) then
268         filename="e107_MY33.dat"
269      else
270         write(*,*)"param_read_e107: "
271         write(*,*)"bad value for solvaryear in callphys.def"
272         write(*,*)"solvaryear must be between 24 and 33"
273         stop
274      endif
275     
276      open(640,file=trim(datafile)//'/EUVDAT/param_v6/'//filename)
277      read(640,*)
278      do i=1,669
279         read(640,*)date_e107(i),e107_tab(i)
280         write(*,*)'param_read_e107/292',date_e107(i),e107_tab(i),i
281      enddo
282      close(640)
283         
284
285c     dissociation and ionization efficiencies
286
287      do inter=1,ninter
288         efdisco2(inter)=0.
289         efdiso2(inter)=0.
290         efdish2(inter)=0.
291         efdish2o(inter)=0.
292         efdish2o2(inter)=0.
293         efdiso3(inter)=0.
294         efdisco(inter)=0.
295         efdisn2(inter)=0.
296         efdisno(inter)=0.
297         efdisno2(inter)=0.
298         efionco2(inter,1)=0.
299         efionco2(inter,2)=0.
300         efionco2(inter,3)=0.
301         efionco2(inter,4)=0.
302         efiono3p(inter)=0.
303         efionn2(inter,1)=0.
304         efionn2(inter,2)=0.
305         efionco(inter,1)=0.
306         efionco(inter,2)=0.
307         efionco(inter,3)=0.
308         efionn(inter)=0.
309         efionh(inter)=0.
310         efionno(inter)=0.
311      enddo
312
313
314c     CO2, O2, NO
315
316      open(120,file=trim(datafile)//'/EUVDAT/param_v5/efdis_inter.dat')
317      read(120,*)
318!      do i=1,21
319!         read(120,*)inter,efdisco2(inter),efdiso2(inter),efdisno(inter)
320      do inter=8,28
321         read(120,*)i,efdisco2(inter),efdiso2(inter),efdisno(inter)
322      enddo
323      do inter=29,ninter
324         efdisco2(inter)=1.
325         efdiso2(inter)=1.
326         efdisno(inter)=1.
327      enddo
328
329
330c     N2
331
332      efdisn2(15)=0.1
333      do inter=16,ninter
334         efdisn2(inter)=1.
335      enddo
336
337
338c     CO
339
340      efdisco(16)=0.5
341      do inter=17,ninter
342         efdisco(inter)=1.
343      enddo
344
345     
346c     O, N, H
347
348      do inter=1,ninter
349         efdiso(inter)=0.
350         efdisn(inter)=0.
351         efdish(inter)=0.
352      enddo
353
354
355c     H2O, H2O2, O3, NO2
356
357      do inter=25,31
358         efdish2o(inter)=1.
359      enddo
360      do inter=25,35
361         efdish2o2(inter)=1.
362      enddo
363      do inter=34,36
364         efdiso3(inter)=1.
365      enddo
366      do inter=27,36
367         efdisno2(inter)=1.
368      enddo
369      do inter=1,15
370         efdish2(inter)=1.
371      enddo
372         
373      !4 possible channels for CO2 ionization
374      do inter=14,16
375         efionco2(inter,1)=1.-efdisco2(inter)
376      enddo
377      efionco2(13,1)=0.805*(1.-efdisco2(13))
378      efionco2(13,2)=0.195*(1.-efdisco2(13))
379      do inter=11,12
380         efionco2(inter,3)=1.-efdisco2(inter)
381      enddo
382      efionco2(10,3)=0.9*(1.-efdisco2(10))
383      efionco2(10,4)=0.1*(1.-efdisco2(10))
384      do inter=2,9
385         efionco2(inter,4)=1.-efdisco2(inter)
386      enddo
387
388      !For O(3p) total ionization under 91.1 nm
389      do inter=1,16
390         efiono3p(inter)=1.d0
391      enddo
392
393      !2 channels for N2 ionization
394      do inter=9,15
395         efionn2(inter,1)=1.-efdisn2(inter)
396      enddo
397      do inter=2,8
398         efionn2(inter,2)=1.-efdisn2(inter)
399      enddo
400     
401      !3 channels for CO ionization
402      do inter=11,16
403         efionco(inter,1)=1.-efdisco(inter)
404      enddo
405      efionco(10,1)=0.87*(1.-efdisco(10))
406      efionco(10,2)=0.13*(1.-efdisco(10))
407      do inter=8,9
408         efionco(inter,2)=1.-efdisco(inter)
409      enddo
410      efionco(7,2)=0.1*(1.-efdisco(7))
411      efionco(7,3)=0.9*(1.-efdisco(7))
412      do inter=2,6
413         efionco(inter,3)=1.-efdisco(inter)
414      enddo
415
416      !Total ionization under 85 nm for N
417      do inter=1,16
418         efionn(inter)=1.
419      enddo
420
421      !NO
422      do inter=2,28
423         efionno(inter)=1.-efdisno(inter)
424      enddo
425
426      !Total ionization under 90 nm for H
427      do inter=3,16
428         efionh(inter)=1.
429      enddo
430
431
432      return
433
434
435      end
436
Note: See TracBrowser for help on using the repository browser.