source: trunk/LMDZ.MARS/libf/aeronomars/jthermcalc.F @ 784

Last change on this file since 784 was 658, checked in by emillour, 13 years ago

Mars GCM:

  • updated high atmosphere photochemistry (jthermcalc.F, param_v4.h, iono.h, paramfoto_compact.F, param_read.F , thermosphere.F).
  • minor change in calchim.F90 (to not use maxloc(zycol, dim = 2) function which seems to be a problem for g95) .
  • minor bug fix in perosat.F; set tendency on pdqscloud for h2o2 to zero if none is computed.
  • in "moldiff.F", changed "tridag" to "tridag_sp", "LUBKSB" to "LUBKSB_SP" and "LUDCMP" to "LUDCMP_SP" to avoid possible conflicts with same routines defined in "moldiff_red.F". Added use of automatic-sized array in "tridag" and "tridag_sp" local array "gam" (to avoid using an a priori oversized local array).

FGG+JYC+EM
File size: 59.0 KB
Line 
1c**********************************************************************
2
3      subroutine jthermcalc(ig,chemthermod,rm,nesptherm,tx,iz,zenit)
4
5
6c     feb 2002        fgg           first version
7c     nov 2002        fgg           second version
8c
9c modified from paramhr.F
10c MAC July 2003
11c**********************************************************************
12
13      implicit none
14
15c     common variables and constants
16      include "dimensions.h"
17      include "dimphys.h"
18      include 'param.h'
19      include 'param_v4.h'
20
21c     input and output variables
22
23      integer    ig
24      integer    chemthermod
25      integer    nesptherm                      !Number of species considered
26      real       rm(nlayermx,nesptherm)         !Densities (cm-3)
27      real       tx(nlayermx)                   !temperature
28      real       zenit                          !SZA
29      real       iz(nlayermx)                   !Local altitude
30
31
32c    local parameters and variables
33
34      real       co2colx(nlayermx)              !column density of CO2 (cm^-2)
35      real       o2colx(nlayermx)               !column density of O2(cm^-2)
36      real       o3pcolx(nlayermx)              !column density of O(3P)(cm^-2)
37      real       h2colx(nlayermx)               !H2 column density (cm-2)
38      real       h2ocolx(nlayermx)              !H2O column density (cm-2)
39      real       h2o2colx(nlayermx)             !column density of H2O2(cm^-2)
40      real       o3colx(nlayermx)               !O3 column density (cm-2)
41      real       n2colx(nlayermx)               !N2 column density (cm-2)
42      real       ncolx(nlayermx)                !N column density (cm-2)
43      real       nocolx(nlayermx)               !NO column density (cm-2)
44      real       cocolx(nlayermx)               !CO column density (cm-2)
45      real       hcolx(nlayermx)                !H column density (cm-2)
46      real       no2colx(nlayermx)              !NO2 column density (cm-2)
47      real       t2(nlayermx)
48      real       coltemp(nlayermx)
49      real       sigma(ninter,nlayermx)
50      real       alfa(ninter,nlayermx)
51     
52      integer    i,j,k,indexint                 !indexes
53      character  dn
54
55
56
57c     variables used in interpolation
58
59      real*8      auxcoltab(nz2)
60      real*8      auxjco2(nz2)
61      real*8      auxjo2(nz2)
62      real*8      auxjo3p(nz2)
63      real*8      auxjh2o(nz2)
64      real*8      auxjh2(nz2)
65      real*8      auxjh2o2(nz2)
66      real*8      auxjo3(nz2)
67      real*8      auxjn2(nz2)
68      real*8      auxjn(nz2)
69      real*8      auxjno(nz2)
70      real*8      auxjco(nz2)
71      real*8      auxjh(nz2)
72      real*8      auxjno2(nz2)
73      real*8      wp(nlayermx),wm(nlayermx)
74      real*8      auxcolinp(nlayermx)
75      integer     auxind(nlayermx)
76      integer     auxi
77      integer     ind
78      real*8      cortemp(nlayermx)
79
80      real*8      limdown                      !limits for interpolation
81      real*8      limup                        !        ""
82
83      !!!ATTENTION. Here i_co2 has to have the same value than in chemthermos.F90
84      !!!If the value is changed there, if has to be changed also here !!!!
85      integer,parameter :: i_co2=1
86
87
88c*************************PROGRAM STARTS*******************************
89     
90      if(zenit.gt.140.) then
91         dn='n'
92         else
93         dn='d'
94      end if
95      if(dn.eq.'n') then
96        return
97      endif
98     
99      !Initializing the photoabsorption coefficients
100      jfotsout(:,:,:)=0.
101
102      !Auxiliar temperature to take into account the temperature dependence
103      !of CO2 cross section
104      do i=1,nlayermx
105         t2(i)=tx(i)
106         if(t2(i).lt.195.0) t2(i)=195.0
107         if(t2(i).gt.295.0) t2(i)=295.0
108      end do
109
110      !Calculation of column amounts
111      call column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
112     $     co2colx,o2colx,o3pcolx,h2colx,h2ocolx,
113     $     h2o2colx,o3colx,n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
114
115      !Auxiliar column to include the temperature dependence
116      !of CO2 cross section
117      coltemp(nlayermx)=co2colx(nlayermx)*abs(t2(nlayermx)-t0(nlayermx))
118      do i=nlayermx-1,1,-1
119        coltemp(i)=!coltemp(i+1)+     PQ SE ELIMINA? REVISAR
120     $         ( rm(i,i_co2) + rm(i+1,i_co2) ) * 0.5
121     $         * 1e5 * (iz(i+1)-iz(i)) * abs(t2(i)-t0(i))
122      end do
123     
124      !Calculation of CO2 cross section at temperature t0(i)
125      do i=1,nlayermx
126         do indexint=24,32
127           sigma(indexint,i)=co2crsc195(indexint-23)
128           alfa(indexint,i)=((co2crsc295(indexint-23)
129     $          /sigma(indexint,i))-1.)/(295.-t0(i))
130        end do
131      end do
132
133! Interpolation to the tabulated photoabsorption coefficients for each species
134! in each spectral interval
135
136
137c     auxcolinp-> Actual atmospheric column
138c     auxj*-> Tabulated photoabsorption coefficients
139c     auxcoltab-> Tabulated atmospheric columns
140
141ccccccccccccccccccccccccccccccc
142c     0.1,5.0 (int 1)
143c
144c     Absorption by:
145c     CO2, O2, O, H2, N
146ccccccccccccccccccccccccccccccc
147
148c     Input atmospheric column
149      indexint=1
150      do i=1,nlayermx
151         auxcolinp(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint) +
152     $        o2colx(i)*crscabsi2(2,indexint) +
153     $        o3pcolx(i)*crscabsi2(3,indexint) +
154     $        h2colx(i)*crscabsi2(5,indexint) +
155     $        ncolx(i)*crscabsi2(9,indexint)
156      end do
157      limdown=1.e-20
158      limup=1.e26
159
160
161c     Interpolations
162
163      do i=1,nz2
164         auxi = nz2-i+1
165         !CO2 tabulated coefficient
166         auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
167         !O2 tabulated coefficient
168         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
169         !O3p tabulated coefficient
170         auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
171         !H2 tabulated coefficient
172         auxjh2(i) = jabsifotsintpar(auxi,5,indexint)
173         !Tabulated column
174         auxcoltab(i) = c1_16(auxi,indexint)
175      enddo
176      !Only if chemthermod.ge.2
177      !N tabulated coefficient
178      if(chemthermod.ge.2) then
179         do i=1,nz2
180            auxjn(i) = jabsifotsintpar(nz2-i+1,9,indexint)
181         enddo
182      endif
183
184      call interfast
185     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
186      do i=1,nlayermx
187         ind=auxind(i)
188         auxi=nlayermx-i+1
189         !CO2 interpolated coefficient
190         jfotsout(indexint,1,auxi) = wm(i)*auxjco2(ind+1) +
191     $        wp(i)*auxjco2(ind)
192         !O2 interpolated coefficient
193         jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
194     $        wp(i)*auxjo2(ind)
195         !O3p interpolated coefficient
196          jfotsout(indexint,3,auxi) = wm(i)*auxjo3p(ind+1) +
197     $        wp(i)*auxjo3p(ind)
198          !H2 interpolated coefficient
199          jfotsout(indexint,5,auxi) = wm(i)*auxjh2(ind+1) +
200     $         wp(i)*auxjh2(ind)
201      enddo
202      !Only if chemthermod.ge.2
203      !N interpolated coefficient
204      if(chemthermod.ge.2) then
205         do i=1,nlayermx
206            ind=auxind(i)
207            jfotsout(indexint,9,nlayermx-i+1) =  wm(i)*auxjn(ind+1) +
208     $         wp(i)*auxjn(ind)
209         enddo
210      endif
211         
212
213c     End interval 1
214
215
216ccccccccccccccccccccccccccccccc
217c     5-80.5nm (int 2-15)
218c
219c     Absorption by:
220c     CO2, O2, O, H2, N2, N,
221c     NO, CO, H, NO2
222ccccccccccccccccccccccccccccccc
223
224c     Input atmospheric column
225      do indexint=2,15
226         do i=1,nlayermx
227            auxcolinp(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
228     $           o2colx(i)*crscabsi2(2,indexint)+
229     $           o3pcolx(i)*crscabsi2(3,indexint)+
230     $           h2colx(i)*crscabsi2(5,indexint)+
231     $           n2colx(i)*crscabsi2(8,indexint)+
232     $           ncolx(i)*crscabsi2(9,indexint)+
233     $           nocolx(i)*crscabsi2(10,indexint)+
234     $           cocolx(i)*crscabsi2(11,indexint)+
235     $           hcolx(i)*crscabsi2(12,indexint)+
236     $           no2colx(i)*crscabsi2(13,indexint)
237         end do
238
239c     Interpolations
240
241         do i=1,nz2
242            auxi = nz2-i+1
243            !O2 tabulated coefficient
244            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
245            !O3p tabulated coefficient
246            auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
247            !CO2 tabulated coefficient
248            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
249            !H2 tabulated coefficient
250            auxjh2(i) = jabsifotsintpar(auxi,5,indexint)
251            !N2 tabulated coefficient
252            auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
253            !CO tabulated coefficient
254            auxjco(i) = jabsifotsintpar(auxi,11,indexint)
255            !H tabulated coefficient
256            auxjh(i) = jabsifotsintpar(auxi,12,indexint)
257            !tabulated column
258            auxcoltab(i) = c1_16(auxi,indexint)
259         enddo
260         !Only if chemthermod.ge.2
261         if(chemthermod.ge.2) then
262            do i=1,nz2
263               auxi = nz2-i+1
264               !N tabulated coefficient
265               auxjn(i) = jabsifotsintpar(auxi,9,indexint)
266               !NO tabulated coefficient
267               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
268               !NO2 tabulated coefficient
269               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
270            enddo
271         endif
272
273          call interfast(wm,wp,auxind,auxcolinp,nlayermx,
274     $        auxcoltab,nz2,limdown,limup)
275          do i=1,nlayermx
276             ind=auxind(i)
277             auxi = nlayermx-i+1
278             !O2 interpolated coefficient
279             jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
280     $            wp(i)*auxjo2(ind)
281             !O3p interpolated coefficient
282             jfotsout(indexint,3,auxi) = wm(i)*auxjo3p(ind+1) +
283     $            wp(i)*auxjo3p(ind)
284             !CO2 interpolated coefficient
285             jfotsout(indexint,1,auxi) = wm(i)*auxjco2(ind+1) +
286     $            wp(i)*auxjco2(ind)
287             !H2 interpolated coefficient
288             jfotsout(indexint,5,auxi) = wm(i)*auxjh2(ind+1) +
289     $            wp(i)*auxjh2(ind)
290             !N2 interpolated coefficient
291             jfotsout(indexint,8,auxi) = wm(i)*auxjn2(ind+1) +
292     $            wp(i)*auxjn2(ind)             
293             !CO interpolated coefficient
294             jfotsout(indexint,11,auxi) = wm(i)*auxjco(ind+1) +
295     $            wp(i)*auxjco(ind)
296             !H interpolated coefficient
297             jfotsout(indexint,12,auxi) = wm(i)*auxjh(ind+1) +
298     $            wp(i)*auxjh(ind)             
299          enddo
300          !Only if chemthermod.ge.2
301          if(chemthermod.ge.2) then
302             do i=1,nlayermx
303                ind=auxind(i)
304                auxi = nlayermx-i+1
305                !N interpolated coefficient
306                jfotsout(indexint,9,auxi) = wm(i)*auxjn(ind+1) +
307     $               wp(i)*auxjn(ind)
308                !NO interpolated coefficient
309                jfotsout(indexint,10,auxi)=wm(i)*auxjno(ind+1) +
310     $               wp(i)*auxjno(ind)
311                !NO2 interpolated coefficient
312                jfotsout(indexint,13,auxi)=wm(i)*auxjno2(ind+1)+
313     $               wp(i)*auxjno2(ind)
314             enddo
315          endif   
316      end do
317
318c     End intervals 2-15
319
320
321ccccccccccccccccccccccccccccccc
322c     80.6-90.8nm (int16)
323c
324c     Absorption by:
325c     CO2, O2, O, N2, N, NO,
326c     CO, H, NO2
327ccccccccccccccccccccccccccccccc
328
329c     Input atmospheric column
330      indexint=16
331      do i=1,nlayermx
332         auxcolinp(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
333     $        o2colx(i)*crscabsi2(2,indexint)+
334     $        o3pcolx(i)*crscabsi2(3,indexint)+
335     $        n2colx(i)*crscabsi2(8,indexint)+
336     $        ncolx(i)*crscabsi2(9,indexint)+
337     $        nocolx(i)*crscabsi2(10,indexint)+
338     $        cocolx(i)*crscabsi2(11,indexint)+
339     $        hcolx(i)*crscabsi2(12,indexint)+
340     $        no2colx(i)*crscabsi2(13,indexint)
341      end do
342
343c     Interpolations
344
345      do i=1,nz2
346         auxi = nz2-i+1
347         !O2 tabulated coefficient
348         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
349         !CO2 tabulated coefficient
350         auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
351         !O3p tabulated coefficient
352         auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
353         !N2 tabulated coefficient
354         auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
355         !CO tabulated coefficient
356         auxjco(i) = jabsifotsintpar(auxi,11,indexint)
357         !H tabulated coefficient
358         auxjh(i) = jabsifotsintpar(auxi,12,indexint)
359         !NO2 tabulated coefficient
360         auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
361         !Tabulated column
362         auxcoltab(i) = c1_16(auxi,indexint)
363      enddo
364      !Only if chemthermod.ge.2
365      if(chemthermod.ge.2) then
366         do i=1,nz2
367            auxi = nz2-i+1
368            !N tabulated coefficient
369            auxjn(i) = jabsifotsintpar(auxi,9,indexint)
370            !NO tabulated coefficient
371            auxjno(i) = jabsifotsintpar(auxi,10,indexint)
372            !NO2 tabulated coefficient
373            auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
374         enddo
375      endif
376
377      call interfast
378     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
379      do i=1,nlayermx
380         ind=auxind(i)
381         auxi = nlayermx-i+1
382         !O2 interpolated coefficient
383         jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
384     $            wp(i)*auxjo2(ind)
385         !CO2 interpolated coefficient
386         jfotsout(indexint,1,auxi) = wm(i)*auxjco2(ind+1) +
387     $            wp(i)*auxjco2(ind)
388         !O3p interpolated coefficient
389         jfotsout(indexint,3,auxi) = wm(i)*auxjo3p(ind+1) +
390     $            wp(i)*auxjo3p(ind)
391         !N2 interpolated coefficient
392         jfotsout(indexint,8,auxi) = wm(i)*auxjn2(ind+1) +
393     $            wp(i)*auxjn2(ind)
394         !CO interpolated coefficient
395         jfotsout(indexint,11,auxi) = wm(i)*auxjco(ind+1) +
396     $            wp(i)*auxjco(ind) 
397         !H interpolated coefficient
398         jfotsout(indexint,12,auxi) = wm(i)*auxjh(ind+1) +
399     $            wp(i)*auxjh(ind)         
400      enddo
401      !Only if chemthermod.ge.2
402      if(chemthermod.ge.2) then
403         do i=1,nlayermx
404            ind=auxind(i)
405            auxi = nlayermx-i+1
406            !N interpolated coefficient
407            jfotsout(indexint,9,auxi) = wm(i)*auxjn(ind+1) +
408     $           wp(i)*auxjn(ind)
409            !NO interpolated coefficient
410            jfotsout(indexint,10,auxi) = wm(i)*auxjno(ind+1) +
411     $           wp(i)*auxjno(ind)
412            !NO2 interpolated coefficient
413            jfotsout(indexint,13,auxi) = wm(i)*auxjno2(ind+1) +
414     $           wp(i)*auxjno2(ind)
415         enddo
416      endif
417c     End interval 16
418
419
420ccccccccccccccccccccccccccccccc
421c     90.9-119.5nm (int 17-24)
422c
423c     Absorption by:
424c     CO2, O2, N2, NO, CO, NO2
425ccccccccccccccccccccccccccccccc
426
427c     Input column
428
429      do i=1,nlayermx
430         auxcolinp(nlayermx-i+1) = co2colx(i) + o2colx(i) + n2colx(i) +
431     $        nocolx(i) + cocolx(i) + no2colx(i)
432      end do
433
434      do indexint=17,24
435
436c     Interpolations
437
438         do i=1,nz2
439            auxi = nz2-i+1
440            !CO2 tabulated coefficient
441            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
442            !O2 tabulated coefficient
443            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
444            !N2 tabulated coefficient
445            auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
446            !CO tabulated coefficient
447            auxjco(i) = jabsifotsintpar(auxi,11,indexint)           
448            !Tabulated column
449            auxcoltab(i) = c17_24(auxi)
450         enddo
451         !Only if chemthermod.ge.2
452         if(chemthermod.ge.2) then
453            do i=1,nz2
454               auxi = nz2-i+1
455               !NO tabulated coefficient
456               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
457               !NO2 tabulated coefficient
458               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
459            enddo
460         endif
461
462         call interfast
463     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
464         !Correction to include T variation of CO2 cross section
465         if(indexint.eq.24) then
466            do i=1,nlayermx
467               auxi = nlayermx-i+1
468               if(sigma(indexint,auxi)*
469     $              alfa(indexint,auxi)*coltemp(auxi)
470     $              .lt.60.) then
471                  cortemp(i)=exp(-sigma(indexint,auxi)*
472     $                alfa(indexint,auxi)*coltemp(auxi))
473               else
474                  cortemp(i)=0.
475               end if
476            enddo
477         else
478            do i=1,nlayermx
479               cortemp(i)=1.
480            enddo
481         end if
482         do i=1,nlayermx           
483            ind=auxind(i)
484            auxi = nlayermx-i+1
485            !O2 interpolated coefficient
486            jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
487     $           wp(i)*auxjo2(ind)) * cortemp(i)
488            !CO2 interpolated coefficient
489            jfotsout(indexint,1,auxi) = (wm(i)*auxjco2(ind+1) +
490     $           wp(i)*auxjco2(ind)) * cortemp(i)
491            if(indexint.eq.24) jfotsout(indexint,1,auxi)=
492     $           jfotsout(indexint,1,auxi)*
493     $           (1+alfa(indexint,auxi)*
494     $           (t2(auxi)-t0(auxi)))
495            !N2 interpolated coefficient
496            jfotsout(indexint,8,auxi) = (wm(i)*auxjn2(ind+1) +
497     $            wp(i)*auxjn2(ind)) * cortemp(i)           
498            !CO interpolated coefficient
499            jfotsout(indexint,11,auxi) = (wm(i)*auxjco(ind+1) +
500     $            wp(i)*auxjco(ind)) * cortemp(i)           
501         enddo
502         !Only if chemthermod.ge.2
503         if(chemthermod.ge.2) then
504            do i=1,nlayermx
505               ind=auxind(i)
506               auxi = nlayermx-i+1
507               !NO interpolated coefficient
508               jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
509     $              wp(i)*auxjno(ind)) * cortemp(i)
510               !NO2 interpolated coefficient
511               jfotsout(indexint,13,auxi)=(wm(i)*auxjno2(ind+1)+
512     $              wp(i)*auxjno2(ind)) * cortemp(i)
513            enddo
514         endif               
515      end do
516c     End intervals 17-24
517
518
519ccccccccccccccccccccccccccccccc
520c     119.6-167.0nm (int 25-29)
521c
522c     Absorption by:
523c     CO2, O2, H2O, H2O2, NO,
524c     CO, NO2
525ccccccccccccccccccccccccccccccc
526
527c     Input atmospheric column
528
529      do i=1,nlayermx
530         auxcolinp(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
531     $        h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
532      end do
533
534      do indexint=25,29
535
536c     Interpolations
537
538         do i=1,nz2
539            auxi = nz2-i+1
540            !CO2 tabulated coefficient
541            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
542            !O2 tabulated coefficient
543            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
544            !H2O tabulated coefficient
545            auxjh2o(i) = jabsifotsintpar(auxi,4,indexint)
546            !H2O2 tabulated coefficient
547            auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)           
548            !CO tabulated coefficient
549            auxjco(i) = jabsifotsintpar(auxi,11,indexint)           
550            !Tabulated column
551            auxcoltab(i) = c25_29(auxi)
552         enddo
553         !Only if chemthermod.ge.2
554         if(chemthermod.ge.2) then
555            do i=1,nz2
556               auxi = nz2-i+1
557               !NO tabulated coefficient
558               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
559               !NO2 tabulated coefficient
560               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
561            enddo
562         endif
563         call interfast
564     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
565         do i=1,nlayermx
566            ind=auxind(i)
567            auxi = nlayermx-i+1
568            !Correction to include T variation of CO2 cross section
569            if(sigma(indexint,auxi)*alfa(indexint,auxi)*
570     $           coltemp(auxi).lt.60.) then
571               cortemp(i)=exp(-sigma(indexint,auxi)*
572     $              alfa(indexint,auxi)*coltemp(auxi))
573            else
574               cortemp(i)=0.
575            end if
576            !CO2 interpolated coefficient
577            jfotsout(indexint,1,auxi) = (wm(i)*auxjco2(ind+1) +
578     $           wp(i)*auxjco2(ind)) * cortemp(i) *
579     $           (1+alfa(indexint,auxi)*
580     $           (t2(auxi)-t0(auxi)))
581            !O2 interpolated coefficient
582            jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
583     $           wp(i)*auxjo2(ind)) * cortemp(i)
584            !H2O interpolated coefficient
585            jfotsout(indexint,4,auxi) = (wm(i)*auxjh2o(ind+1) +
586     $           wp(i)*auxjh2o(ind)) * cortemp(i)
587            !H2O2 interpolated coefficient
588            jfotsout(indexint,6,auxi) = (wm(i)*auxjh2o2(ind+1) +
589     $           wp(i)*auxjh2o2(ind)) * cortemp(i)           
590            !CO interpolated coefficient
591            jfotsout(indexint,11,auxi) = (wm(i)*auxjco(ind+1) +
592     $           wp(i)*auxjco(ind)) * cortemp(i)
593         enddo
594         !Only if chemthermod.ge.2
595         if(chemthermod.ge.2) then
596            do i=1,nlayermx
597               ind=auxind(i)
598               auxi = nlayermx-i+1
599               !NO interpolated coefficient
600               jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
601     $              wp(i)*auxjno(ind)) * cortemp(i)
602               !NO2 interpolated coefficient
603               jfotsout(indexint,13,auxi)=(wm(i)*auxjno2(ind+1)+
604     $              wp(i)*auxjno2(ind)) * cortemp(i)
605            enddo
606         endif
607
608      end do
609
610c     End intervals 25-29
611
612
613cccccccccccccccccccccccccccccccc
614c     167.1-202.5nm (int 30-31)
615c   
616c     Absorption by:
617c     CO2, O2, H2O, H2O2, NO,
618c     NO2
619cccccccccccccccccccccccccccccccc
620
621c     Input atmospheric column
622
623      do i=1,nlayermx
624         auxcolinp(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
625     $        h2o2colx(i) + nocolx(i) + no2colx(i)
626      end do
627
628c     Interpolation
629
630      do indexint=30,31
631
632         do i=1,nz2
633            auxi = nz2-i+1
634            !CO2 tabulated coefficient
635            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
636            !O2 tabulated coefficient
637            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
638            !H2O tabulated coefficient
639            auxjh2o(i) = jabsifotsintpar(auxi,4,indexint)
640            !H2O2 tabulated coefficient
641            auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)           
642            !Tabulated column
643            auxcoltab(i) = c30_31(auxi)
644         enddo
645         !Only if chemthermod.ge.2
646         if(chemthermod.ge.2) then
647            do i=1,nz2
648               auxi = nz2-i+1
649               !NO tabulated coefficient
650               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
651               !NO2 tabulated coefficient
652               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
653            enddo
654         endif
655
656         call interfast
657     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
658         do i=1,nlayermx
659            ind=auxind(i)
660            auxi = nlayermx-i+1
661            !Correction to include T variation of CO2 cross section
662            if(sigma(indexint,auxi)*alfa(indexint,auxi)*
663     $           coltemp(auxi).lt.60.) then
664               cortemp(i)=exp(-sigma(indexint,auxi)*
665     $              alfa(indexint,auxi)*coltemp(auxi))
666            else
667               cortemp(i)=0.
668            end if
669            !CO2 interpolated coefficient
670            jfotsout(indexint,1,auxi) = (wm(i)*auxjco2(ind+1) +
671     $           wp(i)*auxjco2(ind)) * cortemp(i) *
672     $           (1+alfa(indexint,auxi)*
673     $           (t2(auxi)-t0(auxi)))
674            !O2 interpolated coefficient
675            jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
676     $            wp(i)*auxjo2(ind)) * cortemp(i)
677            !H2O interpolated coefficient
678            jfotsout(indexint,4,auxi) = (wm(i)*auxjh2o(ind+1) +
679     $            wp(i)*auxjh2o(ind)) * cortemp(i)
680            !H2O2 interpolated coefficient
681            jfotsout(indexint,6,auxi) = (wm(i)*auxjh2o2(ind+1) +
682     $            wp(i)*auxjh2o2(ind)) * cortemp(i)           
683         enddo
684         !Only if chemthermod.ge.2
685         if(chemthermod.ge.2) then
686            do i=1,nlayermx
687               ind=auxind(i)
688               auxi = nlayermx-i+1
689               !NO interpolated coefficient
690               jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
691     $              wp(i)*auxjno(ind)) * cortemp(i)
692               !NO2 interpolated coefficient
693               jfotsout(indexint,13,auxi)=(wm(i)*auxjno2(ind+1)+
694     $              wp(i)*auxjno2(ind)) * cortemp(i)
695            enddo
696         endif
697
698      end do
699
700c     End intervals 30-31
701
702
703ccccccccccccccccccccccccccccccc
704c     202.6-210.0nm (int 32)
705c
706c     Absorption by:
707c     CO2, O2, H2O2, NO, NO2
708ccccccccccccccccccccccccccccccc
709
710c     Input atmospheric column
711
712      indexint=32
713      do i=1,nlayermx
714         auxcolinp(nlayermx-i+1) =co2colx(i) + o2colx(i) + h2o2colx(i) +
715     $        nocolx(i) + no2colx(i)
716      end do
717
718c     Interpolation
719
720      do i=1,nz2
721         auxi = nz2-i+1
722         !CO2 tabulated coefficient
723         auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
724         !O2 tabulated coefficient
725         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
726         !H2O2 tabulated coefficient
727         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)         
728         !Tabulated column
729         auxcoltab(i) = c32(auxi)
730      enddo
731      !Only if chemthermod.ge.2
732      if(chemthermod.ge.2) then
733         do i=1,nz2
734            auxi = nz2-i+1
735            !NO tabulated coefficient
736            auxjno(i) = jabsifotsintpar(auxi,10,indexint)
737            !NO2 tabulated coefficient
738            auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
739         enddo
740      endif
741      call interfast
742     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
743      do i=1,nlayermx
744         ind=auxind(i)
745         auxi = nlayermx-i+1
746         !Correction to include T variation of CO2 cross section
747         if(sigma(indexint,nlayermx-i+1)*alfa(indexint,auxi)*
748     $        coltemp(auxi).lt.60.) then
749            cortemp(i)=exp(-sigma(indexint,auxi)*
750     $           alfa(indexint,auxi)*coltemp(auxi))
751         else
752            cortemp(i)=0.
753         end if
754         !CO2 interpolated coefficient
755         jfotsout(indexint,1,auxi) = (wm(i)*auxjco2(ind+1) +
756     $        wp(i)*auxjco2(ind)) * cortemp(i) *
757     $        (1+alfa(indexint,auxi)*
758     $        (t2(auxi)-t0(auxi)))
759         !O2 interpolated coefficient
760         jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
761     $        wp(i)*auxjo2(ind)) * cortemp(i)
762         !H2O2 interpolated coefficient
763         jfotsout(indexint,6,auxi) = (wm(i)*auxjh2o2(ind+1) +
764     $        wp(i)*auxjh2o2(ind)) * cortemp(i)         
765      enddo
766      !Only if chemthermod.ge.2
767      if(chemthermod.ge.2) then
768         do i=1,nlayermx
769            auxi = nlayermx-i+1
770            ind=auxind(i)
771            !NO interpolated coefficient
772            jfotsout(indexint,10,auxi) = (wm(i)*auxjno(ind+1) +
773     $           wp(i)*auxjno(ind)) * cortemp(i)
774           !NO2 interpolated coefficient
775            jfotsout(indexint,13,auxi) = (wm(i)*auxjno2(ind+1) +
776     $           wp(i)*auxjno2(ind)) * cortemp(i)
777         enddo
778      endif
779
780c     End of interval 32
781
782
783ccccccccccccccccccccccccccccccc
784c     210.1-231.0nm (int 33)
785c     
786c     Absorption by:
787c     O2, H2O2, NO2
788ccccccccccccccccccccccccccccccc
789
790c     Input atmospheric column
791
792      indexint=33
793      do i=1,nlayermx
794         auxcolinp(nlayermx-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
795      end do
796
797c     Interpolation
798
799      do i=1,nz2
800         auxi = nz2-i+1
801         !O2 tabulated coefficient
802         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
803         !H2O2 tabulated coefficient
804         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
805         !Tabulated column
806         auxcoltab(i) = c33(auxi)
807      enddo
808      !Only if chemthermod.ge.2
809      if(chemthermod.ge.2) then
810         do i=1,nz2
811            !NO2 tabulated coefficient
812            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
813         enddo
814      endif
815      call interfast
816     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
817      do i=1,nlayermx
818         ind=auxind(i)
819         auxi = nlayermx-i+1
820         !O2 interpolated coefficient
821         jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
822     $        wp(i)*auxjo2(ind)
823         !H2O2 interpolated coefficient
824         jfotsout(indexint,6,auxi) = wm(i)*auxjh2o2(ind+1) +
825     $        wp(i)*auxjh2o2(ind)         
826      enddo
827      !Only if chemthermod.ge.2
828      if(chemthermod.ge.2) then
829         do i=1,nlayermx
830            ind=auxind(i)
831            !NO2 interpolated coefficient
832            jfotsout(indexint,13,nlayermx-i+1) = wm(i)*auxjno2(ind+1) +
833     $           wp(i)*auxjno2(ind)
834         enddo
835      endif
836
837c     End of interval 33
838
839
840ccccccccccccccccccccccccccccccc
841c     231.1-240.0nm (int 34)
842c
843c     Absorption by:
844c     O2, H2O2, O3, NO2
845ccccccccccccccccccccccccccccccc
846
847c     Input atmospheric column
848
849      indexint=34
850      do i=1,nlayermx
851         auxcolinp(nlayermx-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) +
852     $        no2colx(i)
853      end do
854
855c     Interpolation
856
857      do i=1,nz2
858         auxi = nz2-i+1
859         !O2 tabulated coefficient
860         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
861         !H2O2 tabulated coefficient
862         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
863         !O3 tabulated coefficient
864         auxjo3(i) = jabsifotsintpar(auxi,7,indexint)         
865         !Tabulated column
866         auxcoltab(i) = c34(nz2-i+1)
867      enddo
868      !Only if chemthermod.ge.2
869      if(chemthermod.ge.2) then
870         do i=1,nz2
871            !NO2 tabulated coefficient
872            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
873         enddo
874      endif
875      call interfast
876     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
877      do i=1,nlayermx
878         ind=auxind(i)
879         auxi = nlayermx-i+1
880         !O2 interpolated coefficient
881         jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
882     $        wp(i)*auxjo2(ind)
883         !H2O2 interpolated coefficient
884         jfotsout(indexint,6,auxi) = wm(i)*auxjh2o2(ind+1) +
885     $        wp(i)*auxjh2o2(ind)
886         !O3 interpolated coefficient
887         jfotsout(indexint,7,auxi) = wm(i)*auxjo3(ind+1) +
888     $        wp(i)*auxjo3(ind)         
889      enddo
890      !Only if chemthermod.ge.2
891      if(chemthermod.ge.2) then
892         do i=1,nlayermx
893            ind=auxind(i)
894            !NO2 interpolated coefficient
895            jfotsout(indexint,13,nlayermx-i+1) = wm(i)*auxjno2(ind+1) +
896     $           wp(i)*auxjno2(ind)
897         enddo
898      endif
899
900c     End of interval 34     
901
902
903ccccccccccccccccccccccccccccccc
904c     240.1-337.7nm (int 35)
905c
906c     Absorption by:
907c     H2O2, O3, NO2
908ccccccccccccccccccccccccccccccc
909
910c     Input atmospheric column
911
912      indexint=35
913      do i=1,nlayermx
914         auxcolinp(nlayermx-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
915      end do
916
917c     Interpolation
918
919      do i=1,nz2
920         auxi = nz2-i+1
921         !H2O2 tabulated coefficient
922         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
923         !O3 tabulated coefficient
924         auxjo3(i) = jabsifotsintpar(auxi,7,indexint)
925         !Tabulated column
926         auxcoltab(i) = c35(auxi)
927      enddo
928      !Only if chemthermod.ge.2
929      if(chemthermod.ge.2) then
930         do i=1,nz2
931            !NO2 tabulated coefficient
932            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
933         enddo
934      endif
935      call interfast
936     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
937      do i=1,nlayermx
938         ind=auxind(i)
939         auxi = nlayermx-i+1
940         !H2O2 interpolated coefficient
941         jfotsout(indexint,6,auxi) = wm(i)*auxjh2o2(ind+1) +
942     $        wp(i)*auxjh2o2(ind)
943         !O3 interpolated coefficient
944         jfotsout(indexint,7,auxi) = wm(i)*auxjo3(ind+1) +
945     $        wp(i)*auxjo3(ind)         
946      enddo
947      if(chemthermod.ge.2) then
948         do i=1,nlayermx
949            ind=auxind(i)
950            !NO2 interpolated coefficient
951            jfotsout(indexint,13,nlayermx-i+1) = wm(i)*auxjno2(ind+1) +
952     $           wp(i)*auxjno2(ind)
953         enddo
954      endif
955
956c     End of interval 35
957
958ccccccccccccccccccccccccccccccc
959c     337.8-800.0 nm (int 36)
960c     
961c     Absorption by:
962c     O3, NO2
963ccccccccccccccccccccccccccccccc
964
965c     Input atmospheric column
966
967      indexint=36
968      do i=1,nlayermx
969         auxcolinp(nlayermx-i+1) = o3colx(i) + no2colx(i)
970      end do
971
972c     Interpolation
973
974      do i=1,nz2
975         auxi = nz2-i+1
976         !O3 tabulated coefficient
977         auxjo3(i) = jabsifotsintpar(auxi,7,indexint)         
978         !Tabulated column
979         auxcoltab(i) = c36(auxi)
980      enddo
981      !Only if chemthermod.ge.2
982      if(chemthermod.ge.2) then
983         do i=1,nz2
984            !NO2 tabulated coefficient
985            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
986         enddo
987      endif
988      call interfast
989     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
990      do i=1,nlayermx
991         ind=auxind(i)
992         !O3 interpolated coefficient
993         jfotsout(indexint,7,nlayermx-i+1) = wm(i)*auxjo3(ind+1) +
994     $        wp(i)*auxjo3(ind)         
995      enddo
996      !Only if chemthermod.ge.2
997      if(chemthermod.ge.2) then
998         do i=1,nlayermx
999            ind=auxind(i)
1000            !NO2 interpolated coefficient
1001            jfotsout(indexint,13,nlayermx-i+1) = wm(i)*auxjno2(ind+1) +
1002     $           wp(i)*auxjno2(ind)
1003         enddo
1004      endif
1005
1006c     End of interval 36
1007
1008c     End of interpolation to obtain photoabsorption rates
1009
1010
1011      return
1012
1013      end
1014
1015
1016
1017c**********************************************************************
1018c**********************************************************************
1019
1020      subroutine column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
1021     $     co2colx,o2colx,o3pcolx,h2colx,h2ocolx,h2o2colx,o3colx,
1022     $     n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
1023
1024c     nov 2002        fgg           first version
1025
1026c**********************************************************************
1027
1028      implicit none
1029
1030
1031c     common variables and constants
1032      include "dimensions.h"
1033      include "dimphys.h"
1034      include "tracer.h"
1035      include 'param.h'
1036      include 'param_v4.h'
1037      include 'callkeys.h'
1038
1039
1040
1041c    local parameters and variables
1042
1043
1044
1045c     input and output variables
1046
1047      integer    ig
1048      integer    chemthermod
1049      integer    nesptherm                      !# of species undergoing chemistry, input
1050      real       rm(nlayermx,nesptherm)         !densities (cm-3), input
1051      real       tx(nlayermx)                   !temperature profile, input
1052      real       iz(nlayermx+1)                 !height profile, input
1053      real       zenit                          !SZA, input
1054      real       co2colx(nlayermx)              !column density of CO2 (cm^-2), output
1055      real       o2colx(nlayermx)               !column density of O2(cm^-2), output
1056      real       o3pcolx(nlayermx)              !column density of O(3P)(cm^-2), output
1057      real       h2colx(nlayermx)               !H2 column density (cm-2), output
1058      real       h2ocolx(nlayermx)              !H2O column density (cm-2), output
1059      real       h2o2colx(nlayermx)             !column density of H2O2(cm^-2), output
1060      real       o3colx(nlayermx)               !O3 column density (cm-2), output
1061      real       n2colx(nlayermx)               !N2 column density (cm-2), output
1062      real       ncolx(nlayermx)                !N column density (cm-2), output
1063      real       nocolx(nlayermx)               !NO column density (cm-2), output
1064      real       cocolx(nlayermx)               !CO column density (cm-2), output
1065      real       hcolx(nlayermx)                !H column density (cm-2), output
1066      real       no2colx(nlayermx)              !NO2 column density (cm-2), output
1067     
1068
1069c     local variables
1070
1071      real       xx
1072      real       grav(nlayermx)
1073      real       Hco2,Ho3p,Ho2,Hh2,Hh2o,Hh2o2
1074      real       Ho3,Hn2,Hn,Hno,Hco,Hh,Hno2
1075
1076      real       co2x(nlayermx)
1077      real       o2x(nlayermx)
1078      real       o3px(nlayermx)
1079      real       cox(nlayermx)
1080      real       hx(nlayermx)
1081      real       h2x(nlayermx)
1082      real       h2ox(nlayermx)
1083      real       h2o2x(nlayermx)
1084      real       o3x(nlayermx)
1085      real       n2x(nlayermx)
1086      real       nx(nlayermx)
1087      real       nox(nlayermx)
1088      real       no2x(nlayermx)
1089
1090      integer    i,j,k,icol,indexint          !indexes
1091
1092c     variables for optical path calculation
1093
1094      integer    nz3
1095!      parameter  (nz3=nz*2)
1096
1097      integer    jj
1098      real*8      esp(nlayermx*2)
1099      real*8      ilayesp(nlayermx*2)
1100      real*8      szalayesp(nlayermx*2)
1101      integer     nlayesp
1102      real*8      zmini
1103      real*8      depth
1104      real*8      espco2, espo2, espo3p, esph2, esph2o, esph2o2,espo3
1105      real*8      espn2,espn,espno,espco,esph,espno2
1106      real*8      rcmnz, rcmmini
1107      real*8      szadeg
1108
1109
1110      ! Tracer indexes in the thermospheric chemistry:
1111      !!! ATTENTION. These values have to be identical to those in chemthermos.F90
1112      !!! If the values are changed there, the same has to be done here  !!!
1113      integer,parameter :: i_co2=1
1114      integer,parameter :: i_o2=2
1115      integer,parameter :: i_o=3
1116      integer,parameter :: i_co=4
1117      integer,parameter :: i_h=5
1118      integer,parameter :: i_h2=8
1119      integer,parameter :: i_h2o=9
1120      integer,parameter :: i_h2o2=10
1121      integer,parameter :: i_o3=12
1122      integer,parameter :: i_n2=13
1123      integer,parameter :: i_n=14
1124      integer,parameter :: i_no=15
1125      integer,parameter :: i_no2=17
1126
1127
1128c*************************PROGRAM STARTS*******************************
1129
1130      nz3 = nlayermx*2
1131      do i=1,nlayermx
1132         xx = ( radio + iz(i) ) * 1.e5
1133         grav(i) = gg * masa /(xx**2)
1134      end do
1135
1136      !Scale heights
1137      xx = kboltzman * tx(nlayermx) * n_avog / grav(nlayermx)
1138      Ho3p  = xx / mmol(igcm_o)
1139      Hco2  = xx / mmol(igcm_co2)
1140      Ho2   = xx / mmol(igcm_o2)
1141      Hh2   = xx / mmol(igcm_h2)
1142      Hh2o  = xx / mmol(igcm_h2o_vap)
1143      Hh2o2 = xx / mmol(igcm_h2o2)
1144      Hco   = xx / mmol(igcm_co)
1145      Hh    = xx / mmol(igcm_h)
1146      !Only if O3 chem. required
1147      if(chemthermod.ge.1)
1148     $     Ho3   = xx / mmol(igcm_o3)
1149      !Only if N or ion chem.
1150      if(chemthermod.ge.2) then
1151         Hn2   = xx / mmol(igcm_n2)
1152         Hn    = xx / mmol(igcm_n)
1153         Hno   = xx / mmol(igcm_no)
1154         Hno2  = xx / mmol(igcm_no2)
1155      endif
1156      ! first loop in altitude : initialisation
1157      do i=nlayermx,1,-1
1158         !Column initialisation
1159         co2colx(i)  = 0.
1160         o2colx(i)   = 0.
1161         o3pcolx(i)  = 0.
1162         h2colx(i)   = 0.
1163         h2ocolx(i)  = 0.
1164         h2o2colx(i) = 0.
1165         o3colx(i)   = 0.
1166         n2colx(i)   = 0.
1167         ncolx(i)    = 0.
1168         nocolx(i)   = 0.
1169         cocolx(i)   = 0.
1170         hcolx(i)    = 0.
1171         no2colx(i)  = 0.
1172         !Densities
1173         co2x(i)  = rm(i,i_co2)
1174         o2x(i)   = rm(i,i_o2)
1175         o3px(i)  = rm(i,i_o)
1176         h2x(i)   = rm(i,i_h2)
1177         h2ox(i)  = rm(i,i_h2o)
1178         h2o2x(i) = rm(i,i_h2o2)
1179         cox(i)   = rm(i,i_co)
1180         hx(i)    = rm(i,i_h)
1181         !Only if O3 chem. required
1182         if(chemthermod.ge.1)
1183     $        o3x(i)   = rm(i,i_o3)
1184         !Only if Nitrogen of ion chemistry requested
1185         if(chemthermod.ge.2) then
1186            n2x(i)   = rm(i,i_n2)
1187            nx(i)    = rm(i,i_n)
1188            nox(i)   = rm(i,i_no)
1189            no2x(i)  = rm(i,i_no2)
1190         endif
1191      enddo
1192      ! second loop in altitude : column calculations
1193      do i=nlayermx,1,-1
1194         !Routine to calculate the geometrical length of each layer
1195         call espesor_optico_A(ig,i,zenit,iz(i),nz3,iz,esp,ilayesp,
1196     $         szalayesp,nlayesp, zmini)
1197         if(ilayesp(nlayesp).eq.-1) then
1198            co2colx(i)=1.e25
1199            o2colx(i)=1.e25
1200            o3pcolx(i)=1.e25
1201            h2colx(i)=1.e25
1202            h2ocolx(i)=1.e25
1203            h2o2colx(i)=1.e25
1204            o3colx(i)=1.e25
1205            n2colx(i)=1.e25
1206            ncolx(i)=1.e25
1207            nocolx(i)=1.e25
1208            cocolx(i)=1.e25
1209            hcolx(i)=1.e25
1210            no2colx(i)=1.e25
1211         else
1212            rcmnz = ( radio + iz(nlayermx) ) * 1.e5
1213            rcmmini = ( radio + zmini ) * 1.e5
1214            !Column calculation taking into account the geometrical depth
1215            !calculated before
1216            do j=1,nlayesp
1217               jj=ilayesp(j)
1218               !Top layer
1219               if(jj.eq.nlayermx) then
1220                  if(zenit.le.60.) then
1221                     o3pcolx(i)=o3pcolx(i)+o3px(nlayermx)*Ho3p*esp(j)
1222     $                    *1.e-5
1223                     co2colx(i)=co2colx(i)+co2x(nlayermx)*Hco2*esp(j)
1224     $                    *1.e-5
1225                     h2o2colx(i)=h2o2colx(i)+
1226     $                    h2o2x(nlayermx)*Hh2o2*esp(j)*1.e-5
1227                     o2colx(i)=o2colx(i)+o2x(nlayermx)*Ho2*esp(j)
1228     $                    *1.e-5
1229                     h2colx(i)=h2colx(i)+h2x(nlayermx)*Hh2*esp(j)
1230     $                    *1.e-5
1231                     h2ocolx(i)=h2ocolx(i)+h2ox(nlayermx)*Hh2o*esp(j)
1232     $                    *1.e-5                     
1233                     cocolx(i)=cocolx(i)+cox(nlayermx)*Hco*esp(j)
1234     $                    *1.e-5
1235                     hcolx(i)=hcolx(i)+hx(nlayermx)*Hh*esp(j)
1236     $                    *1.e-5
1237                     !Only if O3 chemistry required
1238                     if(chemthermod.ge.1) o3colx(i)=
1239     $                    o3colx(i)+o3x(nlayermx)*Ho3*esp(j)
1240     $                    *1.e-5
1241                     !Only if N or ion chemistry requested
1242                     if(chemthermod.ge.2) then
1243                        n2colx(i)=n2colx(i)+n2x(nlayermx)*Hn2*esp(j)
1244     $                    *1.e-5
1245                        ncolx(i)=ncolx(i)+nx(nlayermx)*Hn*esp(j)
1246     $                       *1.e-5
1247                        nocolx(i)=nocolx(i)+nox(nlayermx)*Hno*esp(j)
1248     $                       *1.e-5
1249                        no2colx(i)=no2colx(i)+no2x(nlayermx)*Hno2*esp(j)
1250     $                       *1.e-5
1251                     endif
1252                  else if(zenit.gt.60.) then
1253                     espco2 =sqrt((rcmnz+Hco2)**2 -rcmmini**2) - esp(j)
1254                     espo2  = sqrt((rcmnz+Ho2)**2 -rcmmini**2) - esp(j)
1255                     espo3p = sqrt((rcmnz+Ho3p)**2 -rcmmini**2)- esp(j)
1256                     esph2  = sqrt((rcmnz+Hh2)**2 -rcmmini**2) - esp(j)
1257                     esph2o = sqrt((rcmnz+Hh2o)**2 -rcmmini**2)- esp(j)
1258                     esph2o2= sqrt((rcmnz+Hh2o2)**2-rcmmini**2)- esp(j)
1259                     espco  = sqrt((rcmnz+Hco)**2 -rcmmini**2) - esp(j)
1260                     esph   = sqrt((rcmnz+Hh)**2 -rcmmini**2)  - esp(j)
1261                     !Only if O3 chemistry required
1262                     if(chemthermod.ge.1)                     
1263     $                   espo3=sqrt((rcmnz+Ho3)**2-rcmmini**2)-esp(j)
1264                     !Only if N or ion chemistry requested
1265                     if(chemthermod.ge.2) then
1266                        espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j)
1267                        espn  =sqrt((rcmnz+Hn)**2-rcmmini**2)  - esp(j)
1268                        espno =sqrt((rcmnz+Hno)**2-rcmmini**2) - esp(j)
1269                        espno2=sqrt((rcmnz+Hno2)**2-rcmmini**2)- esp(j)
1270                     endif
1271                     co2colx(i) = co2colx(i) + espco2*co2x(nlayermx)
1272                     o2colx(i)  = o2colx(i)  + espo2*o2x(nlayermx)
1273                     o3pcolx(i) = o3pcolx(i) + espo3p*o3px(nlayermx)
1274                     h2colx(i)  = h2colx(i)  + esph2*h2x(nlayermx)
1275                     h2ocolx(i) = h2ocolx(i) + esph2o*h2ox(nlayermx)
1276                     h2o2colx(i)= h2o2colx(i)+ esph2o2*h2o2x(nlayermx)
1277                     cocolx(i)  = cocolx(i)  + espco*cox(nlayermx)
1278                     hcolx(i)   = hcolx(i)   + esph*hx(nlayermx)
1279                     !Only if O3 chemistry required
1280                     if(chemthermod.ge.1)                     
1281     $                  o3colx(i) = o3colx(i)  + espo3*o3x(nlayermx)
1282                     !Only if N or ion chemistry requested
1283                     if(chemthermod.ge.2) then
1284                        n2colx(i)  = n2colx(i)  + espn2*n2x(nlayermx)
1285                        ncolx(i)   = ncolx(i)   + espn*nx(nlayermx)
1286                        nocolx(i)  = nocolx(i)  + espno*nox(nlayermx)
1287                        no2colx(i) = no2colx(i) + espno2*no2x(nlayermx)
1288                     endif
1289                  endif !Of if zenit.lt.60
1290               !Other layers
1291               else
1292                  co2colx(i)  = co2colx(i) +
1293     $                 esp(j) * (co2x(jj)+co2x(jj+1)) / 2.
1294                  o2colx(i)   = o2colx(i) +
1295     $                 esp(j) * (o2x(jj)+o2x(jj+1)) / 2.
1296                  o3pcolx(i)  = o3pcolx(i) +
1297     $                 esp(j) * (o3px(jj)+o3px(jj+1)) / 2.
1298                  h2colx(i)   = h2colx(i) +
1299     $                 esp(j) * (h2x(jj)+h2x(jj+1)) / 2.
1300                  h2ocolx(i)  = h2ocolx(i) +
1301     $                 esp(j) * (h2ox(jj)+h2ox(jj+1)) / 2.
1302                  h2o2colx(i) = h2o2colx(i) +
1303     $                 esp(j) * (h2o2x(jj)+h2o2x(jj+1)) / 2.
1304                  cocolx(i)   = cocolx(i) +
1305     $                 esp(j) * (cox(jj)+cox(jj+1)) / 2.
1306                  hcolx(i)    = hcolx(i) +
1307     $                 esp(j) * (hx(jj)+hx(jj+1)) / 2.
1308                  !Only if O3 chemistry required
1309                  if(chemthermod.ge.1)
1310     $                 o3colx(i) = o3colx(i) +
1311     $                 esp(j) * (o3x(jj)+o3x(jj+1)) / 2.
1312                  !Only if N or ion chemistry requested
1313                  if(chemthermod.ge.2) then
1314                     n2colx(i)   = n2colx(i) +
1315     $                 esp(j) * (n2x(jj)+n2x(jj+1)) / 2.
1316                     ncolx(i)    = ncolx(i) +
1317     $                    esp(j) * (nx(jj)+nx(jj+1)) / 2.
1318                     nocolx(i)   = nocolx(i) +
1319     $                    esp(j) * (nox(jj)+nox(jj+1)) / 2.
1320                     no2colx(i)  = no2colx(i) +
1321     $                    esp(j) * (no2x(jj)+no2x(jj+1)) / 2.
1322                  endif
1323               endif  !Of if jj.eq.nlayermx
1324            end do    !Of do j=1,nlayesp
1325         endif        !Of ilayesp(nlayesp).eq.-1
1326      enddo           !Of do i=nlayermx,1,-1
1327      return
1328
1329
1330      end
1331
1332
1333c**********************************************************************
1334c**********************************************************************
1335
1336      subroutine interfast(wm,wp,nm,p,nlayer,pin,nl,limdown,limup)
1337C
1338C subroutine to perform linear interpolation in pressure from 1D profile
1339C escin(nl) sampled on pressure grid pin(nl) to profile
1340C escout(nlayer) on pressure grid p(nlayer).
1341C
1342      real*8 wm(nlayer),wp(nlayer),p(nlayer)
1343      integer nm(nlayer)
1344      real*8 pin(nl)
1345      real*8 limup,limdown
1346      integer nl,nlayer,n1,n,np,nini
1347      nini=1
1348      do n1=1,nlayer
1349         if(p(n1) .gt. limup .or. p(n1) .lt. limdown) then
1350            wm(n1) = 0.d0
1351            wp(n1) = 0.d0
1352         else
1353            do n = nini,nl-1
1354               if (p(n1).ge.pin(n).and.p(n1).le.pin(n+1)) then
1355                  nm(n1)=n
1356                  np=n+1
1357                  wm(n1)=abs(pin(n)-p(n1))/(pin(np)-pin(n))
1358                  wp(n1)=1.d0 - wm(n1)
1359                  nini = n
1360                  exit
1361               endif
1362            enddo
1363         endif
1364      enddo
1365      return
1366      end
1367
1368
1369c**********************************************************************
1370c**********************************************************************
1371
1372      subroutine espesor_optico_A (ig,capa, szadeg,z,
1373     @                   nz3,iz,esp,ilayesp,szalayesp,nlayesp, zmini)
1374
1375c     fgg              nov 03      Adaptation to Martian model
1376c     malv             jul 03      Corrected z grid. Split in alt & frec codes
1377c     fgg              feb 03      first version
1378*************************************************************************
1379
1380      implicit none
1381
1382
1383c     common variables and constants
1384
1385      include "dimensions.h"
1386      include "dimphys.h"
1387      include 'param.h'
1388      include 'param_v4.h'
1389
1390c     arguments
1391
1392      real        szadeg                ! I. SZA [rad]
1393      real        z                     ! I. altitude of interest [km]
1394      integer     nz3,ig                   ! I. dimension of esp, ylayesp, etc...
1395                                        !  (=2*nlayermx= max# of layers in ray path)
1396      real     iz(nlayermx+1)              ! I. Altitude of each layer
1397      real*8        esp(nz3)            ! O. layer widths after geometrically
1398                                        !    amplified; in [cm] except at TOA
1399                                        !    where an auxiliary value is used
1400      real*8        ilayesp(nz3)        ! O. Indexes of layers along ray path
1401      real*8        szalayesp(nz3)      ! O. SZA [deg]    "     "       "
1402      integer       nlayesp
1403!      real*8        nlayesp             ! O. # layers along ray path at this z
1404      real*8        zmini               ! O. Minimum altitud of ray path [km]
1405
1406
1407c     local variables and constants
1408
1409        integer     j,i,capa
1410        integer     jmin                  ! index of min.altitude along ray path
1411        real*8      szarad                ! SZA [deg]
1412        real*8      zz
1413        real*8      diz(nlayermx+1)
1414        real*8      rkmnz                 ! distance TOA to center of Planet [km]
1415        real*8      rkmmini               ! distance zmini to center of P [km]
1416        real*8      rkmj                  ! intermediate distance to C of P [km]
1417
1418c external function
1419        external  grid_R8          ! Returns index of layer containing the altitude
1420                                ! of interest, z; for example, if
1421                                ! zkm(i)=z or zkm(i)<z<zkm(i+1) => grid(z)=i
1422        integer   grid_R8
1423
1424*************************************************************************     
1425        szarad = dble(szadeg)*3.141592d0/180.d0
1426        zz=dble(z)
1427        do i=1,nlayermx
1428           diz(i)=dble(iz(i))
1429        enddo
1430        do j=1,nz3
1431          esp(j) = 0.d0
1432          szalayesp(j) = 777.d0
1433          ilayesp(j) = 0
1434        enddo
1435        nlayesp = 0
1436
1437        ! First case: szadeg<60
1438        ! The optical thickness will be given by  1/cos(sza)
1439        ! We deal with 2 different regions:
1440        !   1: First, all layers between z and ztop ("upper part of ray")
1441        !   2: Second, the layer at ztop
1442        if(szadeg.lt.60.d0) then
1443
1444           zmini = zz
1445           if(abs(zz-diz(nlayermx)).lt.1.d-3) goto 1357
1446           ! 1st Zone: Upper part of ray
1447           !
1448           do j=grid_R8(zz,diz,nlayermx),nlayermx-1
1449             nlayesp = nlayesp + 1
1450             ilayesp(nlayesp) = j
1451             esp(nlayesp) = (diz(j+1)-diz(j)) / cos(szarad)        ! [km]
1452             esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
1453             szalayesp(nlayesp) = szadeg
1454           end do
1455
1456           !
1457           ! 2nd Zone: Top layer
1458 1357      continue
1459           nlayesp = nlayesp + 1
1460           ilayesp(nlayesp) = nlayermx
1461           esp(nlayesp) = 1.d0 / cos(szarad)         ! aux. non-dimens. factor
1462           szalayesp(nlayesp) = szadeg
1463
1464
1465        ! Second case:  60 < szadeg < 90
1466        ! The optical thickness is evaluated.
1467        !    (the magnitude of the effect of not using cos(sza) is 3.e-5
1468        !     for z=60km & sza=30, and 5e-4 for z=60km & sza=60, approximately)
1469        ! We deal with 2 different regions:
1470        !   1: First, all layers between z and ztop ("upper part of ray")
1471        !   2: Second, the layer at ztop ("uppermost layer")
1472        else if(szadeg.le.90.d0.and.szadeg.ge.60.d0) then
1473
1474           zmini=(radio+zz)*sin(szarad)-radio
1475           rkmmini = radio + zmini
1476
1477           if(abs(zz-diz(nlayermx)).lt.1.d-4) goto 1470
1478
1479           ! 1st Zone: Upper part of ray
1480           !
1481           do j=grid_R8(zz,diz,nlayermx),nlayermx-1
1482              nlayesp = nlayesp + 1
1483              ilayesp(nlayesp) = j
1484              esp(nlayesp) =
1485     #             sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) -
1486     #             sqrt( (radio+diz(j))**2 - rkmmini**2 )           ! [km]
1487              esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
1488              rkmj = radio+diz(j)
1489              szalayesp(nlayesp) = asin( rkmmini/rkmj )             ! [rad]
1490              szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592 ! [deg]
1491           end do
1492 1470      continue
1493           ! 2nd Zone:  Uppermost layer of ray.
1494           !
1495           nlayesp = nlayesp + 1
1496           ilayesp(nlayesp) = nlayermx
1497           rkmnz = radio+diz(nlayermx)
1498           esp(nlayesp)  =  sqrt( rkmnz**2 - rkmmini**2 )       ! aux.factor[km]
1499           esp(nlayesp)  =  esp(nlayesp) * 1.d5                 ! aux.f. [cm]
1500           szalayesp(nlayesp) = asin( rkmmini/rkmnz )           ! [rad]
1501           szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592! [deg]
1502
1503
1504        ! Third case:  szadeg > 90
1505        ! The optical thickness is evaluated.
1506        ! We deal with 5 different regions:
1507        !   1: all layers between z and ztop ("upper part of ray")
1508        !   2: the layer at ztop ("uppermost layer")
1509        !   3: the lowest layer, at zmini
1510        !   4: the layers increasing from zmini to z (here SZA<90)
1511        !   5: the layers decreasing from z to zmini (here SZA>90)
1512        else if(szadeg.gt.90.d0) then
1513
1514           zmini=(radio+zz)*sin(szarad)-radio
1515           rkmmini = radio + zmini
1516
1517           if(zmini.lt.diz(1)) then         ! Can see the sun?  No => esp(j)=inft
1518             nlayesp = nlayesp + 1
1519             ilayesp(nlayesp) = - 1     ! Value to mark "no sun on view"
1520!             esp(nlayesp) = 1.e30
1521
1522           else
1523              jmin=grid_R8(zmini,diz,nlayermx)+1
1524             
1525
1526              if(abs(zz-diz(nlayermx)).lt.1.d-4) goto 9876
1527
1528              ! 1st Zone: Upper part of ray
1529              !
1530              do j=grid_R8(zz,diz,nlayermx),nlayermx-1
1531                nlayesp = nlayesp + 1
1532                ilayesp(nlayesp) = j
1533                esp(nlayesp) =
1534     $                sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) -
1535     $                sqrt( (radio+diz(j))**2 - rkmmini**2 )          ! [km]
1536                esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
1537                rkmj = radio+diz(j)
1538                szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
1539                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592      ! [deg]
1540              end do
1541
1542 9876         continue
1543              ! 2nd Zone:  Uppermost layer of ray.
1544              !
1545              nlayesp = nlayesp + 1
1546              ilayesp(nlayesp) = nlayermx
1547              rkmnz = radio+diz(nlayermx)
1548              esp(nlayesp) =  sqrt( rkmnz**2 - rkmmini**2 )      !aux.factor[km]
1549              esp(nlayesp) = esp(nlayesp) * 1.d5                 !aux.f.[cm]
1550              szalayesp(nlayesp) = asin( rkmmini/rkmnz )           ! [rad]
1551              szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
1552
1553              ! 3er Zone: Lowestmost layer of ray
1554              !
1555              if ( jmin .ge. 2 ) then      ! If above the planet's surface
1556                j=jmin-1
1557                nlayesp = nlayesp + 1
1558                ilayesp(nlayesp) = j
1559                esp(nlayesp) = 2. *
1560     $                 sqrt( (radio+diz(j+1))**2 -rkmmini**2 )       ! [km]
1561                esp(nlayesp) = esp(nlayesp) * 1.d5                   ! [cm]
1562                rkmj = radio+diz(j+1)
1563                szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad]
1564                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
1565              endif
1566
1567              ! 4th zone: Lower part of ray, increasing from zmin to z
1568              !    ( layers with SZA < 90 deg )
1569              do j=jmin,grid_R8(zz,diz,nlayermx)-1
1570                nlayesp = nlayesp + 1
1571                ilayesp(nlayesp) = j
1572                esp(nlayesp) =
1573     $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
1574     $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )       ! [km]
1575                esp(nlayesp) = esp(nlayesp) * 1.d5                     ! [cm]
1576                rkmj = radio+diz(j)
1577                szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
1578                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
1579              end do
1580
1581              ! 5th zone: Lower part of ray, decreasing from z to zmin
1582              !    ( layers with SZA > 90 deg )
1583              do j=grid_R8(zz,diz,nlayermx)-1, jmin, -1
1584                nlayesp = nlayesp + 1
1585                ilayesp(nlayesp) = j
1586                esp(nlayesp) =
1587     $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
1588     $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )        ! [km]
1589                esp(nlayesp) = esp(nlayesp) * 1.d5                      ! [cm]
1590                rkmj = radio+diz(j)
1591                szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj )          ! [rad]
1592                szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592 ! [deg]
1593              end do
1594
1595           end if
1596
1597        end if
1598
1599        return
1600
1601        end
1602
1603
1604
1605c**********************************************************************
1606c***********************************************************************
1607
1608        function grid_R8 (z, zgrid, nz)
1609
1610c Returns the index where z is located within vector zgrid
1611c The vector zgrid must be monotonously increasing, otherwise program stops.
1612c If z is outside zgrid limits, or zgrid dimension is nz<2, the program stops.
1613c
1614c FGG     Aug-2004     Correct z.lt.zgrid(i) to .le.
1615c MALV    Jul-2003
1616c***********************************************************************
1617
1618        implicit none
1619
1620c Arguments
1621        integer   nz
1622        real*8      z
1623        real*8      zgrid(nz)
1624        integer   grid_R8
1625
1626c Local 
1627        integer   i, nz1, nznew
1628
1629c*** CODE START
1630
1631        if ( z .lt. zgrid(1) .or. z.gt.zgrid(nz) ) then
1632           write (*,*) ' GRID/ z outside bounds of zgrid '
1633           write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz)
1634           stop ' Serious error in GRID.F '
1635        endif
1636        if ( nz .lt. 2 ) then
1637           write (*,*) ' GRID/ zgrid needs 2 points at least ! '
1638           stop ' Serious error in GRID.F '
1639        endif
1640        if ( zgrid(1) .ge. zgrid(nz) ) then
1641           write (*,*) ' GRID/ zgrid must increase with index'
1642           stop ' Serious error in GRID.F '
1643        endif
1644
1645        nz1 = 1
1646        nznew = nz/2
1647        if ( z .gt. zgrid(nznew) ) then
1648           nz1 = nznew
1649           nznew = nz
1650        endif
1651        do i=nz1+1,nznew
1652           if ( z. eq. zgrid(i) ) then
1653              grid_R8=i
1654              return
1655              elseif ( z .le. zgrid(i) ) then
1656              grid_R8 = i-1
1657              return
1658           endif
1659        enddo
1660        grid_R8 = nz
1661        return
1662
1663
1664
1665        end
1666
1667
1668
1669!c***************************************************
1670!c***************************************************
1671
1672      subroutine flujo(date)
1673
1674
1675!c     fgg           nov 2002     first version
1676!c***************************************************
1677
1678      implicit none
1679
1680
1681!     common variables and constants
1682      include "dimensions.h"
1683      include "dimphys.h"
1684      include "comsaison.h"
1685      include 'param.h'
1686      include 'param_v4.h'
1687
1688
1689!     Arguments
1690
1691      real date
1692
1693
1694!     Local variable and constants
1695
1696      integer i
1697      integer inter
1698      real    nada
1699
1700!c*************************************************
1701
1702      if(date.lt.1985.) date=1985.
1703      if(date.gt.2001.) date=2001.
1704     
1705      do i=1,ninter
1706         fluxtop(i)=1.
1707         !Variation of solar flux with 11 years solar cycle
1708         !For more details, see Gonzalez-Galindo et al. 2005
1709         !To be improved in next versions
1710        if(i.le.24) then
1711          fluxtop(i)=(((ct1(i)+p1(i)*date)/2.)                 
1712     $        *sin(2.*3.1416/11.*(date-1985.-3.1416))         
1713     $        +(ct2(i)+p2(i)*date)+1.)*fluxtop(i)
1714        end if
1715        fluxtop(i)=fluxtop(i)*(1.52/dist_sol)**2
1716      end do
1717!      fluxtop(1)=fluxtop(1)*10.
1718!      fluxtop(2)=fluxtop(2)*5.
1719     
1720      return
1721      end
Note: See TracBrowser for help on using the repository browser.