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

Last change on this file since 1198 was 705, checked in by emillour, 12 years ago

Mars GCM:

  • Added possibility to run with a varying EUV cycle following real one. The flag solvarmod=1 triggers this behaviour, with companion flag solvaryear=## , where ## is the Mars Year (from 23 to 30). Setting solvarmod=0 reverts to 'old' behaviour, where there is a constant EUV forcing throughout the run, set by the "solarcondate" flag.
  • Needs corresponding input data files ("param_v6" subdirectory of "EUV" subdirectory in "datadir").
  • Added files jthermcalc_e107.F and param_read_e107.F in "aeronomars", modified files euvheat.F90, hrtherm.F, chemthermos.F90, param_v4.h and param_read.F in "aeronomars" and files inifis.F, physiq.F and callkeys.h in "phymars".

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