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

Last change on this file since 3094 was 2836, checked in by abierjon, 2 years ago

VENUS GCM:

INCLUDING THE IONOSPHERE CODE IN VENUS GCM

ATTENTION: INCREASED MEMORY DEMAND

NEEDS AT LEAST 135 GB of allocated memory

==============================================================
===> LIST OF APPENDED FILES AND INTERNAL ADDITIONS <==========
==============================================================

NEW MODEL SPECIES

  • deftank/traceur-chemistry-IONOSPHERE.def

Neutrals: 4 Species: N, NO, NO2, N(2D)

Ions: 15 species:

CO2+, CO+, O+, O2+, H+,

N2+, H2O+, OH+, C+, HCO+,

H3O+, HCO2+, N+, NO+, elec

FROM 36 to 55 chemical species

NEW KEYWORD OF PHYSIQ.DEF

  • deftank/physiq-96x96x90-chemistry-IONOSPHERE.def
  • ok_ionchem: keyword supposed to activate ion chemistry.

(be careful that n, no, no2, n2d and the ion species are in the deftank/traceur-chemistry-IONOSPHERE.def)

  • ok_jonline: keyword supposed to activate the online photochemistry

(be careful that n, no, no2 and n2d are in the traceur-chemistry-IONOSPHERE.def)

==============================================================
===> LIST OF MODIFIED PROGRAMS <=========================
==============================================================

nonoro_gwd_ran_mod.F90

  • Change EPFLUXMAX value from 5.E-3 to 1.E-3

photochemistry_venus.F90 (krates; photolysis_online; indices;)

  • Import of keywords: ok_ionchem, tuneupperatm & ok_jonline
  • addition of ion species in order
  • Forcing electroneutrality
  • Update of the reactions a001 and a002 with taking into account the other species
  • Change of formula for a002 with the formulation of Baulch et al., 1976 (confirmed by Smith and Robertson, 2008)

photolysis_mod.F90

  • Modification of the rdsolarflux subroutine to include interpolation with Atlas1 and Atlas3 in connection with E10.7

concentration2.F90

  • Added chemical species n2d, no, no2 and n in the conductivity calculation.

The ions have been excluded because their sum is 105 times less dense than the neutrals and
their thermal conductivity is unknown

iono_h.F90

  • Addition of the phdisrate routine
  • replace the electronic temperature of Mars by that of

origin = 1: Theis et al. 1980 (Venus) with bilinear interpolation altitude/cos(SZA)
origin = 2: Theis et al. 1984 (Venus) with the formula of the electronic temperature at high solar activity

  • addition of an ion temperature model based on VIRA

cleshphys.h

  • added ok_ionchem & ok_jonline in COMMON/clesphys_l/

conf_phys.f90

  • add ok_ionchem & ok_jonline as parameters to read from physiq.def file set to .false. by default

chemparam_mod.F90

  • Add species in M_tr and corresponding i_X. Set all i_X to zero before reading traceur-chemistry-IONOSPHERE.def
  • Added Type_tr table to differentiate species: 1 == neutral, 2 == ION, 3 == liquid, 10 == others


euvheat.F90; hrtherm.F; jthermcalc_e107.F; param_read_e107.F

  • Normalization with Mars

A.M
File size: 13.7 KB
Line 
1      subroutine param_read_e107
2
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
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
32      character (len=100) :: datafile="HIGHATM"
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 datafile=/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         efiono2(inter,1)=0.
297         efiono2(inter,2)=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!      close(120)
325
326
327c     N2
328
329      efdisn2(15)=0.1
330      do inter=16,ninter
331         efdisn2(inter)=1.
332      enddo
333
334
335c     CO
336
337      efdisco(16)=0.5
338      do inter=17,ninter
339         efdisco(inter)=1.
340      enddo
341
342     
343c     O, N, H
344
345      do inter=1,ninter
346         efdiso(inter)=0.
347         efdisn(inter)=0.
348         efdish(inter)=0.
349      enddo
350
351
352c     H2O, H2O2, O3, NO2
353
354      do inter=25,31
355         efdish2o(inter)=1.
356      enddo
357      do inter=25,35
358         efdish2o2(inter)=1.
359      enddo
360      do inter=34,36
361         efdiso3(inter)=1.
362      enddo
363      do inter=27,36
364         efdisno2(inter)=1.
365      enddo
366      do inter=1,15
367         efdish2(inter)=1.
368      enddo
369         
370      !4 possible channels for CO2 ionization
371      open(130,file=trim(datafile)//'/EUVDAT'//
372     $     '/co2ion_branchingratio_schunkandnagy2000_param.dat')
373      do inter=1,16
374         read(130,*)i,nada,efionco2(inter,1),efionco2(inter,2),
375     $        efionco2(inter,3),efionco2(inter,4)
376         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
377         efdisco2(inter)=1.-nada
378         efionco2(inter,1)=(1.-efdisco2(inter))*efionco2(inter,1)
379         efionco2(inter,2)=(1.-efdisco2(inter))*efionco2(inter,2)
380         efionco2(inter,3)=(1.-efdisco2(inter))*efionco2(inter,3)
381         efionco2(inter,4)=(1.-efdisco2(inter))*efionco2(inter,4)
382      enddo
383      close(130)
384
385      do inter=17,36
386         efdisco2(inter)=1.
387      enddo
388
389!      do inter=14,16
390!         efionco2(inter,1)=1.-efdisco2(inter)
391!      enddo
392!      efionco2(13,1)=0.805*(1.-efdisco2(13))
393!      efionco2(13,2)=0.195*(1.-efdisco2(13))
394!      do inter=11,12
395!         efionco2(inter,3)=1.-efdisco2(inter)
396!      enddo
397!      efionco2(10,3)=0.9*(1.-efdisco2(10))
398!      efionco2(10,4)=0.1*(1.-efdisco2(10))
399!      do inter=2,9
400!         efionco2(inter,4)=1.-efdisco2(inter)
401!      enddo
402
403      !2 possible channels for O2 ionization
404      open(131,file=trim(datafile)//'/EUVDAT'//
405     $     '/o2ion_branchingratio_schunkandnagy2000_param.dat')
406      do inter=1,23
407         read(131,*)i,nada,efiono2(inter,1),efiono2(inter,2)
408         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
409         efdiso2(inter)=1.-nada
410         efiono2(inter,1)=(1.-efdiso2(inter))*efiono2(inter,1)
411         efiono2(inter,2)=(1.-efdiso2(inter))*efiono2(inter,2)
412      enddo
413      close(131)
414
415      do inter=24,36
416         efdiso2(inter)=1.
417      enddo
418
419
420      !For O(3p) total ionization under 91.1 nm
421      do inter=1,16
422         efiono3p(inter)=1.d0
423      enddo
424
425      !2 channels for N2 ionization
426      open(132,file=trim(datafile)//'/EUVDAT'//
427     $     '/n2ion_branchingratio_schunkandnagy2000_param.dat')
428      do inter=1,15
429         read(132,*)i,nada,efionn2(inter,1),efionn2(inter,2)
430         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
431         efdisn2(inter)=1.-nada
432         efionn2(inter,1)=(1.-efdisn2(inter))*efionn2(inter,1)
433         efionn2(inter,2)=(1.-efdisn2(inter))*efionn2(inter,2)
434      enddo
435      close(132)
436     
437      do inter=16,36
438         efdisn2(inter)=1.
439      enddo
440
441!      do inter=9,15
442!         efionn2(inter,1)=1.-efdisn2(inter)
443!      enddo
444!      do inter=2,8
445!         efionn2(inter,2)=1.-efdisn2(inter)
446!      enddo
447     
448      !3 channels for CO ionization
449       open(133,file=trim(datafile)//'/EUVDAT'//
450     $     '/coion_branchingratio_schunkandnagy2000_param.dat')
451      do inter=1,16
452         read(133,*)i,nada,efionco(inter,1),efionco(inter,2),
453     $        efionco(inter,3)
454         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
455         efdisco(inter)=1.-nada
456         efionco(inter,1)=(1.-efdisco(inter))*efionco(inter,1)
457         efionco(inter,2)=(1.-efdisco(inter))*efionco(inter,2)
458         efionco(inter,3)=(1.-efdisco(inter))*efionco(inter,3)
459      enddo
460      close(133)
461     
462     
463     
464      do inter=17,36
465         efdisco(inter)=1.
466      enddo
467     
468!      do inter=11,16
469!         efionco(inter,1)=1.-efdisco(inter)
470!      enddo
471!      efionco(10,1)=0.87*(1.-efdisco(10))
472!      efionco(10,2)=0.13*(1.-efdisco(10))
473!      do inter=8,9
474!         efionco(inter,2)=1.-efdisco(inter)
475!      enddo
476!      efionco(7,2)=0.1*(1.-efdisco(7))
477!      efionco(7,3)=0.9*(1.-efdisco(7))
478!      do inter=2,6
479!         efionco(inter,3)=1.-efdisco(inter)
480!      enddo
481
482
483      !Total ionization under 85 nm for N
484      do inter=1,16
485         efionn(inter)=1.
486      enddo
487
488      !NO
489      do inter=2,28
490         efionno(inter)=1.-efdisno(inter)
491      enddo
492
493      !Total ionization under 90 nm for H
494      do inter=3,16
495         efionh(inter)=1.
496      enddo
497
498
499      return
500
501
502      end
503
Note: See TracBrowser for help on using the repository browser.