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

Last change on this file was 3466, checked in by emillour, 2 months ago

Mars PCM:
More tidying in aeronomars:

  • remove unused "inv.F" and remove "dtridgl.F" which is not used here and is a duplicate of the "dtridgl" routine in phymars/swr_toon.F
  • turn aeronomars routines to modules, for those which aren't in modules yet.

EM
File size: 14.5 KB
Line 
1      module param_read_e107_mod
2     
3      implicit none
4     
5      contains
6     
7      subroutine param_read_e107
8
9      use param_v4_h, only: jfotsout,crscabsi2,
10     .    c1_16,c17_24,c25_29,c30_31,c32,c33,c34,c35,c36,
11     .    co2crsc195,co2crsc295,t0,
12     .    jabsifotsintpar,ninter,nz2,
13     .    nabs,e107,date_e107,e107_tab,
14     .    coefit0,coefit1,coefit2,coefit3,coefit4,
15     .    efdisco2,efdiso2,efdish2o,
16     .    efdish2o2,efdish2,efdiso3,
17     .    efdiso,efdisn,efdish,
18     .    efdisno,efdisn2,efdisno2,
19     .    efdisco,efionco2,efiono2,efionn2,
20     .    efionco,efiono3p,efionn,
21     .    efionno,efionh,
22     .    fluxtop,ct1,ct2,p1,p2
23
24      use datafile_mod, only: datadir
25      USE mod_phys_lmdz_para, ONLY: is_master
26      USE mod_phys_lmdz_transfert_para, ONLY: bcast
27
28      implicit none
29
30
31c     common variables and constants
32      include "callkeys.h"
33 
34 
35c     local variables
36
37      integer    i,j,k,inter                          !indexes
38      integer ierr,lnblnk
39      real       nada
40      character*13  filename
41 
42     
43c*************************+PROGRAM STARTS**************************
44
45      if(is_master) then
46
47c     Reads tabulated functions
48
49      !Tabulated column amount
50      open(210, status = 'old',
51c    $file=trim(datadir)//'/EUVDAT/coln.dat',iostat=ierr)
52     $file=trim(datadir)//'/EUVDAT/param_v6/coln.dat',iostat=ierr)
53
54      IF (ierr.NE.0) THEN
55       write(*,*)'cant find directory EUVDAT containing param_v6 subdir'
56       write(*,*)'(in aeronomars/param_read.F)'
57       write(*,*)'It should be in :', trim(datadir),'/'
58       write(*,*)'1) You can change this directory address in '
59       write(*,*)'   callphys.def with datadir=/path/to/dir'
60       write(*,*)'2) If necessary, EUVDAT (and other datafiles)'
61       write(*,*)'   can be obtained online on:'
62       write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
63       STOP
64      ENDIF
65 
66      !Tabulated photoabsorption coefficients
67      open(220,file=trim(datadir)//'/EUVDAT/param_v6/trans2_an.dat')
68      open(230,file=trim(datadir)//'/EUVDAT/param_v6/trans3_an.dat')
69      open(240,file=trim(datadir)//'/EUVDAT/param_v6/trans1_an.dat')
70      open(250,file=trim(datadir)//'/EUVDAT/param_v6/trans2_bn.dat')
71      open(260,file=trim(datadir)//'/EUVDAT/param_v6/trans2_cn.dat')
72      open(300,file=trim(datadir)//'/EUVDAT/param_v6/trans2_dn.dat')
73      open(270,file=trim(datadir)//'/EUVDAT/param_v6/trans1_bn.dat')
74      open(280,file=trim(datadir)//'/EUVDAT/param_v6/trans1_cn.dat')
75      open(290,file=trim(datadir)//'/EUVDAT/param_v6/trans1_dn.dat')
76      open(150,file=trim(datadir)//'/EUVDAT/param_v6/trans4n.dat')
77      open(160,file=trim(datadir)//'/EUVDAT/param_v6/trans5n.dat')
78      open(170,file=trim(datadir)//'/EUVDAT/param_v6/trans6n.dat')
79      open(180,file=trim(datadir)//'/EUVDAT/param_v6/trans7n.dat')
80      open(390,file=trim(datadir)//'/EUVDAT/param_v6/trans8_an.dat')
81      open(400,file=trim(datadir)//'/EUVDAT/param_v6/trans8_bn.dat')
82      open(410,file=trim(datadir)//'/EUVDAT/param_v6/trans9n.dat')
83      open(420,file=trim(datadir)//'/EUVDAT/param_v6/trans10_an.dat')
84      open(430,file=trim(datadir)//'/EUVDAT/param_v6/trans10_bn.dat')
85      open(440,file=trim(datadir)//'/EUVDAT/param_v6/trans10_cn.dat')
86      open(450,file=trim(datadir)//'/EUVDAT/param_v6/trans11_an.dat')
87      open(460,file=trim(datadir)//'/EUVDAT/param_v6/trans11_bn.dat')
88      open(470,file=trim(datadir)//'/EUVDAT/param_v6/trans11_cn.dat')
89      open(480,file=trim(datadir)//'/EUVDAT/param_v6/trans12n.dat')
90      open(490,file=trim(datadir)//'/EUVDAT/param_v6/trans13_an.dat')
91      open(500,file=trim(datadir)//'/EUVDAT/param_v6/trans13_bn.dat')
92      open(510,file=trim(datadir)//'/EUVDAT/param_v6/trans13_cn.dat')
93
94     
95      do i=210,300,10
96         read(i,*)
97         read(i,*)
98      end do
99
100      do i=150,180,10
101         read(i,*)
102         read(i,*)
103      end do
104
105      do i=390,510,10
106         read(i,*)
107         read(i,*)
108      enddo
109
110      do i=nz2,1,-1
111         read(210,*) (c1_16(i,j),j=1,16),c17_24(i),c25_29(i),c30_31(i),
112     $        c32(i),c33(i),c34(i),c35(i),c36(i)
113      end do
114
115      do i=nz2,1,-1
116         read(220,*) (jabsifotsintpar(i,2,j),j=1,16)
117      end do
118     
119      do i=nz2,1,-1
120         read(230,*) (jabsifotsintpar(i,3,j),j=1,16)
121      end do
122
123      do i=nz2,1,-1
124         read(240,*) (jabsifotsintpar(i,1,j),j=1,16)
125      end do
126
127      do i=nz2,1,-1
128         read(250,*) (jabsifotsintpar(i,2,j),j=17,24)
129      end do
130
131
132      do i=nz2,1,-1
133         read(260,*) (jabsifotsintpar(i,2,j),j=25,31)
134      end do
135
136      do i=nz2,1,-1
137         read(270,*) (jabsifotsintpar(i,1,j),j=17,24)
138      end do
139
140      do i=nz2,1,-1
141         read(280,*) (jabsifotsintpar(i,1,j),j=25,31)
142      end do
143
144      do i=nz2,1,-1
145         read(290,*) jabsifotsintpar(i,1,32)
146      end do
147
148      do i=nz2,1,-1
149         read(300,*) (jabsifotsintpar(i,2,j),j=32,34)
150      end do
151
152      do i=nz2,1,-1
153         read(160,*) (jabsifotsintpar(i,5,j),j=1,15)
154      end do
155
156      do i=nz2,1,-1
157         read(150,*) (jabsifotsintpar(i,4,j),j=25,31)
158      end do
159
160      do i=nz2,1,-1
161         read(170,*) (jabsifotsintpar(i,6,j),j=25,35)
162      end do
163
164      do i=nz2,1,-1
165         read(180,*) (jabsifotsintpar(i,7,j),j=34,36)
166      end do
167
168      do i=nz2,1,-1
169         read(390,*) (jabsifotsintpar(i,8,j),j=2,16)
170      enddo
171
172      do i=nz2,1,-1
173         read(400,*) (jabsifotsintpar(i,8,j),j=17,24)
174      enddo
175
176      do i=nz2,1,-1
177         read(410,*) (jabsifotsintpar(i,9,j),j=1,16)
178      enddo
179
180      do i=nz2,1,-1
181         read(420,*) (jabsifotsintpar(i,10,j),j=2,16)
182      enddo
183
184      do i=nz2,1,-1
185         read(430,*) (jabsifotsintpar(i,10,j),j=17,24)
186      enddo
187
188      do i=nz2,1,-1
189         read(440,*) (jabsifotsintpar(i,10,j),j=25,32)
190      enddo
191
192      do i=nz2,1,-1
193         read(450,*) (jabsifotsintpar(i,11,j),j=2,16)
194      enddo
195
196      do i=nz2,1,-1
197         read(460,*) (jabsifotsintpar(i,11,j),j=17,24)
198      enddo
199
200      do i=nz2,1,-1
201         read(470,*) (jabsifotsintpar(i,11,j),j=25,29)
202      enddo
203     
204      do i=nz2,1,-1
205         read(480,*) (jabsifotsintpar(i,12,j),j=2,16)
206      enddo
207
208      do i=nz2,1,-1
209         read(490,*) (jabsifotsintpar(i,13,j),j=2,16)
210      enddo
211     
212      do i=nz2,1,-1
213         read(500,*) (jabsifotsintpar(i,13,j),j=17,24)
214      enddo
215     
216      do i=nz2,1,-1
217         read(510,*) (jabsifotsintpar(i,13,j),j=25,36)
218      enddo
219
220      do i=210,300,10
221         close(i)
222      end do
223
224      do i=150,180,10
225         close(i)
226      end do
227
228      do i=390,510,10
229         close(i)
230      enddo
231
232      endif !is master
233
234c     set t0
235
236      do i=1,nz2
237         t0(i)=195.
238      end do
239
240
241      do i=1,ninter
242         fluxtop(i)=1.
243      end do
244
245      if(is_master) then
246
247      !Parameters for the variation of the solar flux with 11 years cycle
248      open(620,file=trim(datadir)//'/EUVDAT/param_v6/fit_js_e107.dat')
249      do i=1,ninter
250         read(620,*)
251         do j=1,nabs
252            read(620,*)coefit0(i,j),coefit1(i,j),coefit2(i,j),
253     $           coefit3(i,j),coefit4(i,j)
254         enddo
255      enddo
256      close(620)
257
258
259      !Tabulated values for E10.7
260      if((solvaryear.ge.23).and.(solvaryear.le.36)) then
261        write(filename,fmt="(a,i2,a)")"e107_MY",solvaryear,".dat"
262      else
263         write(*,*)"param_read_e107: "
264         write(*,*)"bad value for solvaryear in callphys.def"
265         write(*,*)"solvaryear must be between 24 and 36"
266         stop
267      endif
268     
269      open(640,file=trim(datadir)//'/EUVDAT/param_v6/'//filename)
270      read(640,*)
271      do i=1,669
272         read(640,*)date_e107(i),e107_tab(i)
273         write(*,*)'param_read_e107/292',date_e107(i),e107_tab(i),i
274      enddo
275      close(640)
276
277      endif !is master
278
279      call bcast(t0)
280      call bcast(fluxtop)
281
282      call bcast(c1_16)
283      call bcast(c17_24)
284      call bcast(c25_29)
285      call bcast(c30_31)
286      call bcast(c32)
287      call bcast(c33)
288      call bcast(c34)
289      call bcast(c35)
290      call bcast(c36)
291
292      call bcast(jabsifotsintpar)
293
294      call bcast(coefit0)
295      call bcast(coefit1)
296      call bcast(coefit2)
297      call bcast(coefit3)
298      call bcast(coefit4)
299      call bcast(date_e107)
300      call bcast(e107_tab)
301         
302
303c     dissociation and ionization efficiencies
304
305      do inter=1,ninter
306         efdisco2(inter)=0.
307         efdiso2(inter)=0.
308         efdish2(inter)=0.
309         efdish2o(inter)=0.
310         efdish2o2(inter)=0.
311         efdiso3(inter)=0.
312         efdisco(inter)=0.
313         efdisn2(inter)=0.
314         efdisno(inter)=0.
315         efdisno2(inter)=0.
316         efionco2(inter,1)=0.
317         efionco2(inter,2)=0.
318         efionco2(inter,3)=0.
319         efionco2(inter,4)=0.
320         efiono2(inter,1)=0.
321         efiono2(inter,2)=0.
322         efiono3p(inter)=0.
323         efionn2(inter,1)=0.
324         efionn2(inter,2)=0.
325         efionco(inter,1)=0.
326         efionco(inter,2)=0.
327         efionco(inter,3)=0.
328         efionn(inter)=0.
329         efionh(inter)=0.
330         efionno(inter)=0.
331      enddo
332
333
334c     CO2, O2, NO
335
336!      open(120,file=trim(datadir)//'/EUVDAT/param_v5/efdis_inter.dat')
337!      read(120,*)
338!!      do i=1,21
339!!         read(120,*)inter,efdisco2(inter),efdiso2(inter),efdisno(inter)
340!      do inter=8,28
341!         read(120,*)i,efdisco2(inter),efdiso2(inter),efdisno(inter)
342!      enddo
343!      do inter=29,ninter
344!         efdisco2(inter)=1.
345!         efdiso2(inter)=1.
346!         efdisno(inter)=1.
347!      enddo
348!      close(120)
349
350
351c     N2
352
353      efdisn2(15)=0.1
354      do inter=16,ninter
355         efdisn2(inter)=1.
356      enddo
357
358
359c     CO
360
361      efdisco(16)=0.5
362      do inter=17,ninter
363         efdisco(inter)=1.
364      enddo
365
366     
367c     O, N, H
368
369      do inter=1,ninter
370         efdiso(inter)=0.
371         efdisn(inter)=0.
372         efdish(inter)=0.
373      enddo
374
375
376c     H2O, H2O2, O3, NO2
377
378      do inter=25,31
379         efdish2o(inter)=1.
380      enddo
381      do inter=25,35
382         efdish2o2(inter)=1.
383      enddo
384      do inter=34,36
385         efdiso3(inter)=1.
386      enddo
387      do inter=27,36
388         efdisno2(inter)=1.
389      enddo
390      do inter=1,15
391         efdish2(inter)=1.
392      enddo
393
394      if(is_master) then
395         
396      !4 possible channels for CO2 ionization
397      open(130,file=trim(datadir)//'/EUVDAT'//
398     $     '/co2ion_branchingratio_schunkandnagy2000_param.dat')
399      do inter=1,16
400         read(130,*)i,nada,efionco2(inter,1),efionco2(inter,2),
401     $        efionco2(inter,3),efionco2(inter,4)
402         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
403         efdisco2(inter)=1.-nada
404         efionco2(inter,1)=(1.-efdisco2(inter))*efionco2(inter,1)
405         efionco2(inter,2)=(1.-efdisco2(inter))*efionco2(inter,2)
406         efionco2(inter,3)=(1.-efdisco2(inter))*efionco2(inter,3)
407         efionco2(inter,4)=(1.-efdisco2(inter))*efionco2(inter,4)
408      enddo
409      close(130)
410
411      endif !is_master
412
413      call bcast(i)
414      call bcast(nada)
415      call bcast(efionco2)
416      call bcast(efdisco2)
417
418      do inter=17,36
419         efdisco2(inter)=1.
420      enddo
421
422!      do inter=14,16
423!         efionco2(inter,1)=1.-efdisco2(inter)
424!      enddo
425!      efionco2(13,1)=0.805*(1.-efdisco2(13))
426!      efionco2(13,2)=0.195*(1.-efdisco2(13))
427!      do inter=11,12
428!         efionco2(inter,3)=1.-efdisco2(inter)
429!      enddo
430!      efionco2(10,3)=0.9*(1.-efdisco2(10))
431!      efionco2(10,4)=0.1*(1.-efdisco2(10))
432!      do inter=2,9
433!         efionco2(inter,4)=1.-efdisco2(inter)
434!      enddo
435
436      if(is_master) then
437      !2 possible channels for O2 ionization
438      open(131,file=trim(datadir)//'/EUVDAT'//
439     $     '/o2ion_branchingratio_schunkandnagy2000_param.dat')
440      do inter=1,23
441         read(131,*)i,nada,efiono2(inter,1),efiono2(inter,2)
442         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
443         efdiso2(inter)=1.-nada
444         efiono2(inter,1)=(1.-efdiso2(inter))*efiono2(inter,1)
445         efiono2(inter,2)=(1.-efdiso2(inter))*efiono2(inter,2)
446      enddo
447      close(131)
448
449      endif !is_master
450
451      call bcast(i)
452      call bcast(nada)
453      call bcast(efiono2)
454      call bcast(efdiso2)
455
456      do inter=24,36
457         efdiso2(inter)=1.
458      enddo
459
460
461      !For O(3p) total ionization under 91.1 nm
462      do inter=1,16
463         efiono3p(inter)=1.d0
464      enddo
465
466      if(is_master) then
467
468      !2 channels for N2 ionization
469      open(132,file=trim(datadir)//'/EUVDAT'//
470     $     '/n2ion_branchingratio_schunkandnagy2000_param.dat')
471      do inter=1,15
472         read(132,*)i,nada,efionn2(inter,1),efionn2(inter,2)
473         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
474         efdisn2(inter)=1.-nada
475         efionn2(inter,1)=(1.-efdisn2(inter))*efionn2(inter,1)
476         efionn2(inter,2)=(1.-efdisn2(inter))*efionn2(inter,2)
477      enddo
478      close(132)
479
480      endif
481
482      call bcast(i)
483      call bcast(nada)
484      call bcast(efionn2)
485      call bcast(efdisn2)
486
487      do inter=16,36
488         efdisn2(inter)=1.
489      enddo
490
491!      do inter=9,15
492!         efionn2(inter,1)=1.-efdisn2(inter)
493!      enddo
494!      do inter=2,8
495!         efionn2(inter,2)=1.-efdisn2(inter)
496!      enddo
497 
498      if(is_master) then     
499      !3 channels for CO ionization
500       open(133,file=trim(datadir)//'/EUVDAT'//
501     $     '/coion_branchingratio_schunkandnagy2000_param.dat')
502      do inter=1,16
503         read(133,*)i,nada,efionco(inter,1),efionco(inter,2),
504     $        efionco(inter,3)
505         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
506         efdisco(inter)=1.-nada
507         efionco(inter,1)=(1.-efdisco(inter))*efionco(inter,1)
508         efionco(inter,2)=(1.-efdisco(inter))*efionco(inter,2)
509         efionco(inter,3)=(1.-efdisco(inter))*efionco(inter,3)
510      enddo
511      close(133)
512
513      endif  !     if(is_master)
514
515      call bcast(i)
516      call bcast(nada)
517      call bcast(efionco)
518      call bcast(efdisco)
519
520      do inter=17,36
521         efdisco(inter)=1.
522      enddo
523
524!      do inter=11,16
525!         efionco(inter,1)=1.-efdisco(inter)
526!      enddo
527!      efionco(10,1)=0.87*(1.-efdisco(10))
528!      efionco(10,2)=0.13*(1.-efdisco(10))
529!      do inter=8,9
530!         efionco(inter,2)=1.-efdisco(inter)
531!      enddo
532!      efionco(7,2)=0.1*(1.-efdisco(7))
533!      efionco(7,3)=0.9*(1.-efdisco(7))
534!      do inter=2,6
535!         efionco(inter,3)=1.-efdisco(inter)
536!      enddo
537
538
539      !Total ionization under 85 nm for N
540      do inter=1,16
541         efionn(inter)=1.
542      enddo
543
544      !NO
545      do inter=2,28
546         efionno(inter)=1.-efdisno(inter)
547      enddo
548
549      !Total ionization under 90 nm for H
550      do inter=3,16
551         efionh(inter)=1.
552      enddo
553
554
555      end subroutine param_read_e107
556
557      end module param_read_e107_mod
Note: See TracBrowser for help on using the repository browser.