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

Last change on this file since 1574 was 1500, checked in by emillour, 10 years ago

Mars GCM:

  • set things in param_read_e107.F and read_dust_scenario.F90 to read MY32 EUV and dust forcings.

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