source: trunk/LMDZ.VENUS/libf/phyvenus/param_read_e107.F @ 2487

Last change on this file since 2487 was 2464, checked in by slebonnois, 4 years ago

SL: major update related to extension to 250 km. New EUV heating as used in Martian GCM, debug of a few routines, adding He, clean up, updating mmean regularly, modification of the order of processes, add the possibility to use Hedin composition in rad transfer
File size: 11.2 KB
RevLine 
[1310]1      subroutine param_read_e107
2
[2464]3      use dimphy
4      use param_v4_h, only: jfotsout,crscabsi2,
5     .    c1_16,c17_24,c25_29,c30_31,c32,c33,c34,c35,c36,
6     .    co2crsc195,co2crsc295,t0,
7     .    jabsifotsintpar,ninter,nz2,
8     .    nabs,e107,date_e107,e107_tab,
9     .    coefit0,coefit1,coefit2,coefit3,coefit4,
10     .    efdisco2,efdiso2,efdish2o,
11     .    efdish2o2,efdish2,efdiso3,
12     .    efdiso,efdisn,efdish,
13     .    efdisno,efdisn2,efdisno2,
14     .    efdisco,efionco2,efiono2,efionn2,
15     .    efionco,efiono3p,efionn,
16     .    efionno,efionh,
17     .    fluxtop,ct1,ct2,p1,p2
18
[1310]19      implicit none
20
21
22c     common variables and constants
23      include "clesphys.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
[2464]32      character (len=100) :: datafile="HIGHATM"
[1310]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/~forget/datagcm/datafile'
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
247c      if(solvaryear.eq.23) then
248c         filename="e107_MY23.dat"
249c      else if(solvaryear.eq.24) then
250c         filename="e107_MY24.dat"
251c      else if(solvaryear.eq.25) then
252c         filename="e107_MY25.dat"
253c      else if (solvaryear.eq.26) then
254c         filename="e107_MY26.dat"
255c      else if (solvaryear.eq.27) then
256c         filename="e107_MY27.dat"
257c     else if (solvaryear.eq.28) then
258c        filename="e107_MY28.dat"
259c      else if (solvaryear.eq.29) then
260c         filename="e107_MY29.dat"
261c      else if(solvaryear.eq.30) then
262c         filename="e107_MY30.dat"
263c      else
264c         write(*,*)"param_read_e107: "
265c         write(*,*)"bad value for solvaryear in callphys.def"
266c         write(*,*)"solvaryear must be between 24 and 30"
267c         stop
268c      endif
269     
270c      open(640,file=trim(datafile)//'/EUVDAT/param_v6/'//filename)
271c      read(640,*)
272c      do i=1,669
273c         read(640,*)date_e107(i),e107_tab(i)
274c         write(*,*)'param_read_e107/292',date_e107(i),e107_tab(i),i
275c      enddo
276c      close(640)
277         
278
279c     dissociation and ionization efficiencies
280
281!      do inter=1,ninter
282!         efdisco2(inter)=0.
283!         efdiso2(inter)=0.
284!         efdish2(inter)=0.
285!         efdish2o(inter)=0.
286!         efdish2o2(inter)=0.
287!         efdiso3(inter)=0.
288!         efdisco(inter)=0.
289!         efdisn2(inter)=0.
290!         efdisno(inter)=0.
291!         efdisno2(inter)=0.
292!         efionco2(inter,1)=0.
293!         efionco2(inter,2)=0.
294!         efionco2(inter,3)=0.
295!         efionco2(inter,4)=0.
296!         efiono3p(inter)=0.
297!         efionn2(inter,1)=0.
298!         efionn2(inter,2)=0.
299!         efionco(inter,1)=0.
300!         efionco(inter,2)=0.
301!         efionco(inter,3)=0.
302!         efionn(inter)=0.
303!         efionh(inter)=0.
304!         efionno(inter)=0.
305!      enddo
306
307
308c     CO2, O2, NO
309
310!      open(120,file=trim(datafile)//'/EUVDAT/param_v5/efdis_inter.dat')
311!      read(120,*)
312!!      do i=1,21
313!!         read(120,*)inter,efdisco2(inter),efdiso2(inter),efdisno(inter)
314!      do inter=8,28
315!         read(120,*)i,efdisco2(inter),efdiso2(inter),efdisno(inter)
316!      enddo
317!      do inter=29,ninter
318!         efdisco2(inter)=1.
319!         efdiso2(inter)=1.
320!         efdisno(inter)=1.
321!      enddo
322
323
324c     N2
325
326!      efdisn2(15)=0.1
327!      do inter=16,ninter
328!         efdisn2(inter)=1.
329!      enddo
330
331
332c     CO
333
334!      efdisco(16)=0.5
335!      do inter=17,ninter
336!         efdisco(inter)=1.
337!      enddo
338
339     
340c     O, N, H
341
342!      do inter=1,ninter
343!         efdiso(inter)=0.
344!         efdisn(inter)=0.
345!         efdish(inter)=0.
346!      enddo
347
348
349c     H2O, H2O2, O3, NO2
350
351!      do inter=25,31
352!         efdish2o(inter)=1.
353!      enddo
354!      do inter=25,35
355!         efdish2o2(inter)=1.
356!      enddo
357!      do inter=34,36
358!         efdiso3(inter)=1.
359!      enddo
360!      do inter=27,36
361!         efdisno2(inter)=1.
362!      enddo
363!      do inter=1,15
364!         efdish2(inter)=1.
365!      enddo
366         
367      !4 possible channels for CO2 ionization
368!      do inter=14,16
369!         efionco2(inter,1)=1.-efdisco2(inter)
370!      enddo
371!      efionco2(13,1)=0.805*(1.-efdisco2(13))
372!      efionco2(13,2)=0.195*(1.-efdisco2(13))
373!      do inter=11,12
374!         efionco2(inter,3)=1.-efdisco2(inter)
375!      enddo
376!      efionco2(10,3)=0.9*(1.-efdisco2(10))
377!      efionco2(10,4)=0.1*(1.-efdisco2(10))
378!      do inter=2,9
379!         efionco2(inter,4)=1.-efdisco2(inter)
380!      enddo
381
382      !For O(3p) total ionization under 91.1 nm
383!      do inter=1,16
384!         efiono3p(inter)=1.d0
385!      enddo
386
387      !2 channels for N2 ionization
388!      do inter=9,15
389!         efionn2(inter,1)=1.-efdisn2(inter)
390!      enddo
391!      do inter=2,8
392!         efionn2(inter,2)=1.-efdisn2(inter)
393!      enddo
394     
395      !3 channels for CO ionization
396!      do inter=11,16
397!         efionco(inter,1)=1.-efdisco(inter)
398!      enddo
399!      efionco(10,1)=0.87*(1.-efdisco(10))
400!      efionco(10,2)=0.13*(1.-efdisco(10))
401!      do inter=8,9
402!         efionco(inter,2)=1.-efdisco(inter)
403!      enddo
404!      efionco(7,2)=0.1*(1.-efdisco(7))
405!      efionco(7,3)=0.9*(1.-efdisco(7))
406!      do inter=2,6
407!         efionco(inter,3)=1.-efdisco(inter)
408!      enddo
409
410      !Total ionization under 85 nm for N
411!      do inter=1,16
412!         efionn(inter)=1.
413!      enddo
414
415      !NO
416!      do inter=2,28
417!         efionno(inter)=1.-efdisno(inter)
418!      enddo
419
420      !Total ionization under 90 nm for H
421!      do inter=3,16
422!         efionh(inter)=1.
423!      enddo
424
425
426      return
427
428
429      end
430
Note: See TracBrowser for help on using the repository browser.