source: trunk/LMDZ.MARS/libf/aeronomars/jthermcalc_e107.F @ 3026

Last change on this file since 3026 was 2808, checked in by emillour, 2 years ago

Mars GCM:
Bug fix in jthermcalc_e107 (wrong indexes used in some computations involving
NO2 and H2).
AM
File size: 37.3 KB
Line 
1c**********************************************************************
2
3      subroutine jthermcalc_e107
4     $     (ig,nlayer,chemthermod,rm,nesptherm,tx,iz,zenit,zday)
5
6
7c     feb 2002        fgg           first version
8c     nov 2002        fgg           second version
9c
10c modified from paramhr.F
11c MAC July 2003
12c**********************************************************************
13
14      use param_v4_h, only: jfotsout,crscabsi2,
15     .    c1_16,c17_24,c25_29,c30_31,c32,c33,c34,c35,c36,
16     .    co2crsc195,co2crsc295,t0,
17     .    jabsifotsintpar,ninter,nz2,
18     .    nabs,e107,date_e107,e107_tab,
19     .    coefit0,coefit1,coefit2,coefit3,coefit4
20      use comsaison_h, only: dist_sol
21
22      implicit none
23
24      include "callkeys.h"
25
26c     input and output variables
27
28      integer,intent(in) :: ig,nlayer
29      integer,intent(in) :: chemthermod
30      integer,intent(in) :: nesptherm          !Number of species considered
31      real,intent(in) :: rm(nlayer,nesptherm)  !Densities (cm-3)
32      real,intent(in) :: tx(nlayer)            !temperature
33      real,intent(in) :: zenit         !SZA
34      real,intent(in) :: iz(nlayer)    !Local altitude
35      real,intent(in) :: zday          !Martian day after Ls=0
36
37! NB: no output variable! (computed jfotsout() is in module param_v4_h)
38
39c    local parameters and variables
40
41      real       co2colx(nlayer)              !column density of CO2 (cm^-2)
42      real       o2colx(nlayer)               !column density of O2(cm^-2)
43      real       o3pcolx(nlayer)              !column density of O(3P)(cm^-2)
44      real       h2colx(nlayer)               !H2 column density (cm-2)
45      real       h2ocolx(nlayer)              !H2O column density (cm-2)
46      real       h2o2colx(nlayer)             !column density of H2O2(cm^-2)
47      real       o3colx(nlayer)               !O3 column density (cm-2)
48      real       n2colx(nlayer)               !N2 column density (cm-2)
49      real       ncolx(nlayer)                !N column density (cm-2)
50      real       nocolx(nlayer)               !NO column density (cm-2)
51      real       cocolx(nlayer)               !CO column density (cm-2)
52      real       hcolx(nlayer)                !H column density (cm-2)
53      real       no2colx(nlayer)              !NO2 column density (cm-2)
54      real       t2(nlayer)
55      real       coltemp(nlayer)
56      real       sigma(ninter,nlayer)
57      real       alfa(ninter,nlayer)
58      real       realday
59     
60      integer    i,j,k,indexint                 !indexes
61      character  dn
62      integer    tinf,tsup
63
64
65
66c     variables used in interpolation
67
68      real*8      auxcoltab(nz2)
69      real*8      auxjco2(nz2)
70      real*8      auxjo2(nz2)
71      real*8      auxjo3p(nz2)
72      real*8      auxjh2o(nz2)
73      real*8      auxjh2(nz2)
74      real*8      auxjh2o2(nz2)
75      real*8      auxjo3(nz2)
76      real*8      auxjn2(nz2)
77      real*8      auxjn(nz2)
78      real*8      auxjno(nz2)
79      real*8      auxjco(nz2)
80      real*8      auxjh(nz2)
81      real*8      auxjno2(nz2)
82      real*8      wp(nlayer),wm(nlayer)
83      real*8      auxcolinp(nlayer)
84      integer     auxind(nlayer)
85      integer     auxi
86      integer     ind
87      real*8      cortemp(nlayer)
88
89      real*8      limdown                      !limits for interpolation
90      real*8      limup                        !        ""
91
92      !!!ATTENTION. Here i_co2 has to have the same value than in chemthermos.F90
93      !!!If the value is changed there, if has to be changed also here !!!!
94      integer,parameter :: i_co2=1
95
96
97c*************************PROGRAM STARTS*******************************
98     
99      if(zenit.gt.140.) then
100         dn='n'
101         else
102         dn='d'
103      end if
104      if(dn.eq.'n') then
105        return
106      endif
107     
108      !Initializing the photoabsorption coefficients
109      jfotsout(:,:,:)=0.
110
111      !Auxiliar temperature to take into account the temperature dependence
112      !of CO2 cross section
113      do i=1,nlayer
114         t2(i)=tx(i)
115         if(t2(i).lt.195.0) t2(i)=195.0
116         if(t2(i).gt.295.0) t2(i)=295.0
117      end do
118
119      !Calculation of column amounts
120      call column(ig,nlayer,chemthermod,rm,nesptherm,tx,iz,zenit,
121     $     co2colx,o2colx,o3pcolx,h2colx,h2ocolx,
122     $     h2o2colx,o3colx,n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
123
124      !Auxiliar column to include the temperature dependence
125      !of CO2 cross section
126      coltemp(nlayer)=co2colx(nlayer)*abs(t2(nlayer)-t0(nlayer))
127      do i=nlayer-1,1,-1
128        coltemp(i)=!coltemp(i+1)+     PQ SE ELIMINA? REVISAR
129     $         ( rm(i,i_co2) + rm(i+1,i_co2) ) * 0.5
130     $         * 1e5 * (iz(i+1)-iz(i)) * abs(t2(i)-t0(i))
131      end do
132     
133      !Calculation of CO2 cross section at temperature t0(i)
134      do i=1,nlayer
135         do indexint=24,32
136           sigma(indexint,i)=co2crsc195(indexint-23)
137           alfa(indexint,i)=((co2crsc295(indexint-23)
138     $          /sigma(indexint,i))-1.)/(295.-t0(i))
139        end do
140      end do
141
142      if (solvarmod==0) then
143        e107=fixed_euv_value
144      else
145        !E10.7 for the day: linear interpolation to tabulated values
146        realday=mod(zday,669.)
147        if(realday.lt.date_e107(1)) then
148         e107=e107_tab(1)
149        else if(realday.ge.date_e107(669)) then
150         e107=e107_tab(669)   
151        else if(realday.ge.date_e107(1).and.
152     $        realday.lt.date_e107(669)) then
153         do i=1,668
154            if(realday.ge.date_e107(i).and.
155     $           realday.lt.date_e107(i+1)) then
156               tinf=i
157               tsup=i+1
158               e107=e107_tab(tinf)+(realday-date_e107(tinf))*
159     $              (e107_tab(tsup)-e107_tab(tinf))
160            endif
161         enddo
162        endif
163      endif ! of if (solvarmod==0)
164
165      !Photoabsorption coefficients at TOA as a function of E10.7
166      do j=1,nabs
167         do indexint=1,ninter
168            jfotsout(indexint,j,nlayer)=coefit0(indexint,j)+
169     $           coefit1(indexint,j)*e107+coefit2(indexint,j)*e107**2+
170     $           coefit3(indexint,j)*e107**3+coefit4(indexint,j)*e107**4
171         enddo
172      enddo
173! Interpolation to the tabulated photoabsorption coefficients for each species
174! in each spectral interval
175
176
177c     auxcolinp-> Actual atmospheric column
178c     auxj*-> Tabulated photoabsorption coefficients
179c     auxcoltab-> Tabulated atmospheric columns
180
181ccccccccccccccccccccccccccccccc
182c     0.1,5.0 (int 1)
183c
184c     Absorption by:
185c     CO2, O2, O, H2, N
186ccccccccccccccccccccccccccccccc
187
188c     Input atmospheric column
189      indexint=1
190      do i=1,nlayer
191         auxcolinp(nlayer-i+1) = co2colx(i)*crscabsi2(1,indexint) +
192     $        o2colx(i)*crscabsi2(2,indexint) +
193     $        o3pcolx(i)*crscabsi2(3,indexint) +
194     $        h2colx(i)*crscabsi2(5,indexint) +
195     $        ncolx(i)*crscabsi2(9,indexint)
196      end do
197      limdown=1.e-20
198      limup=1.e26
199
200
201c     Interpolations
202
203      do i=1,nz2
204         auxi = nz2-i+1
205         !CO2 tabulated coefficient
206         auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
207         !O2 tabulated coefficient
208         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
209         !O3p tabulated coefficient
210         auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
211         !H2 tabulated coefficient
212         auxjh2(i) = jabsifotsintpar(auxi,5,indexint)
213         !Tabulated column
214         auxcoltab(i) = c1_16(auxi,indexint)
215      enddo
216      !Only if chemthermod.ge.2
217      !N tabulated coefficient
218      if(chemthermod.ge.2) then
219         do i=1,nz2
220            auxjn(i) = jabsifotsintpar(nz2-i+1,9,indexint)
221         enddo
222      endif
223
224      call interfast
225     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
226     
227      do i=1,nlayer
228         ind=auxind(i)
229         auxi=nlayer-i+1
230         !CO2 interpolated coefficient
231         jfotsout(indexint,1,auxi) = jfotsout(indexint,1,nlayer) *
232     $        (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind))
233         !O2 interpolated coefficient
234         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayer) *
235     $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
236         !O3p interpolated coefficient
237         jfotsout(indexint,3,auxi) = jfotsout(indexint,3,nlayer) *
238     $        (wm(i)*auxjo3p(ind+1) + wp(i)*auxjo3p(ind))
239         !H2 interpolated coefficient
240         jfotsout(indexint,5,auxi) = jfotsout(indexint,5,nlayer) *
241     $        (wm(i)*auxjh2(ind+1) + wp(i)*auxjh2(ind))
242      enddo
243      !Only if chemthermod.ge.2
244      !N interpolated coefficient
245      if(chemthermod.ge.2) then
246         do i=1,nlayer
247            ind=auxind(i)
248            jfotsout(indexint,9,nlayer-i+1) = 
249     $           jfotsout(indexint,9,nlayer) *
250     $           (wm(i)*auxjn(ind+1) + wp(i)*auxjn(ind))
251         enddo
252      endif
253         
254
255c     End interval 1
256
257
258ccccccccccccccccccccccccccccccc
259c     5-80.5nm (int 2-15)
260c
261c     Absorption by:
262c     CO2, O2, O, H2, N2, N,
263c     NO, CO, H, NO2
264ccccccccccccccccccccccccccccccc
265
266c     Input atmospheric column
267      do indexint=2,15
268         do i=1,nlayer
269            auxcolinp(nlayer-i+1) = co2colx(i)*crscabsi2(1,indexint)+
270     $           o2colx(i)*crscabsi2(2,indexint)+
271     $           o3pcolx(i)*crscabsi2(3,indexint)+
272     $           h2colx(i)*crscabsi2(5,indexint)+
273     $           n2colx(i)*crscabsi2(8,indexint)+
274     $           ncolx(i)*crscabsi2(9,indexint)+
275     $           nocolx(i)*crscabsi2(10,indexint)+
276     $           cocolx(i)*crscabsi2(11,indexint)+
277     $           hcolx(i)*crscabsi2(12,indexint)+
278     $           no2colx(i)*crscabsi2(13,indexint)
279         end do
280
281c     Interpolations
282
283         do i=1,nz2
284            auxi = nz2-i+1
285            !O2 tabulated coefficient
286            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
287            !O3p tabulated coefficient
288            auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
289            !CO2 tabulated coefficient
290            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
291            !H2 tabulated coefficient
292            auxjh2(i) = jabsifotsintpar(auxi,5,indexint)
293            !N2 tabulated coefficient
294            auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
295            !CO tabulated coefficient
296            auxjco(i) = jabsifotsintpar(auxi,11,indexint)
297            !H tabulated coefficient
298            auxjh(i) = jabsifotsintpar(auxi,12,indexint)
299            !tabulated column
300            auxcoltab(i) = c1_16(auxi,indexint)
301         enddo
302         !Only if chemthermod.ge.2
303         if(chemthermod.ge.2) then
304            do i=1,nz2
305               auxi = nz2-i+1
306               !N tabulated coefficient
307               auxjn(i) = jabsifotsintpar(auxi,9,indexint)
308               !NO tabulated coefficient
309               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
310               !NO2 tabulated coefficient
311               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
312            enddo
313         endif
314
315          call interfast(wm,wp,auxind,auxcolinp,nlayer,
316     $        auxcoltab,nz2,limdown,limup)
317
318          do i=1,nlayer
319             ind=auxind(i)
320             auxi = nlayer-i+1
321             !O2 interpolated coefficient
322             jfotsout(indexint,2,auxi) =
323     $            jfotsout(indexint,2,nlayer) *
324     $            (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
325             !O3p interpolated coefficient
326             jfotsout(indexint,3,auxi) =
327     $            jfotsout(indexint,3,nlayer) *
328     $            (wm(i)*auxjo3p(ind+1) + wp(i)*auxjo3p(ind))
329             !CO2 interpolated coefficient
330             jfotsout(indexint,1,auxi) =
331     $            jfotsout(indexint,1,nlayer) *
332     $            (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind))
333             !H2 interpolated coefficient
334             jfotsout(indexint,5,auxi) =
335     $            jfotsout(indexint,5,nlayer) *
336     $            (wm(i)*auxjh2(ind+1) + wp(i)*auxjh2(ind))
337             !N2 interpolated coefficient
338             jfotsout(indexint,8,auxi) =
339     $            jfotsout(indexint,8,nlayer) *
340     $            (wm(i)*auxjn2(ind+1) + wp(i)*auxjn2(ind))
341             !CO interpolated coefficient
342             jfotsout(indexint,11,auxi) =
343     $            jfotsout(indexint,11,nlayer) *
344     $            (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind))
345             !H interpolated coefficient
346             jfotsout(indexint,12,auxi) =
347     $            jfotsout(indexint,12,nlayer) *
348     $            (wm(i)*auxjh(ind+1) + wp(i)*auxjh(ind))
349          enddo
350          !Only if chemthermod.ge.2
351          if(chemthermod.ge.2) then
352             do i=1,nlayer
353                ind=auxind(i)
354                auxi = nlayer-i+1
355                !N interpolated coefficient
356                jfotsout(indexint,9,auxi) =
357     $               jfotsout(indexint,9,nlayer) *
358     $               (wm(i)*auxjn(ind+1) + wp(i)*auxjn(ind))
359                !NO interpolated coefficient
360                jfotsout(indexint,10,auxi)=
361     $               jfotsout(indexint,10,nlayer) *
362     $               (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind))
363                !NO2 interpolated coefficient
364                jfotsout(indexint,13,auxi)=
365     $               jfotsout(indexint,13,nlayer) *
366     $               (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
367             enddo
368          endif   
369      end do
370
371c     End intervals 2-15
372
373
374ccccccccccccccccccccccccccccccc
375c     80.6-90.8nm (int16)
376c
377c     Absorption by:
378c     CO2, O2, O, N2, N, NO,
379c     CO, H, NO2
380ccccccccccccccccccccccccccccccc
381
382c     Input atmospheric column
383      indexint=16
384      do i=1,nlayer
385         auxcolinp(nlayer-i+1) = co2colx(i)*crscabsi2(1,indexint)+
386     $        o2colx(i)*crscabsi2(2,indexint)+
387     $        o3pcolx(i)*crscabsi2(3,indexint)+
388     $        n2colx(i)*crscabsi2(8,indexint)+
389     $        ncolx(i)*crscabsi2(9,indexint)+
390     $        nocolx(i)*crscabsi2(10,indexint)+
391     $        cocolx(i)*crscabsi2(11,indexint)+
392     $        hcolx(i)*crscabsi2(12,indexint)+
393     $        no2colx(i)*crscabsi2(13,indexint)
394      end do
395
396c     Interpolations
397
398      do i=1,nz2
399         auxi = nz2-i+1
400         !O2 tabulated coefficient
401         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
402         !CO2 tabulated coefficient
403         auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
404         !O3p tabulated coefficient
405         auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
406         !N2 tabulated coefficient
407         auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
408         !CO tabulated coefficient
409         auxjco(i) = jabsifotsintpar(auxi,11,indexint)
410         !H tabulated coefficient
411         auxjh(i) = jabsifotsintpar(auxi,12,indexint)
412         !NO2 tabulated coefficient
413         auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
414         !Tabulated column
415         auxcoltab(i) = c1_16(auxi,indexint)
416      enddo
417      !Only if chemthermod.ge.2
418      if(chemthermod.ge.2) then
419         do i=1,nz2
420            auxi = nz2-i+1
421            !N tabulated coefficient
422            auxjn(i) = jabsifotsintpar(auxi,9,indexint)
423            !NO tabulated coefficient
424            auxjno(i) = jabsifotsintpar(auxi,10,indexint)
425            !NO2 tabulated coefficient
426            auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
427         enddo
428      endif
429
430      call interfast
431     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
432
433      do i=1,nlayer
434         ind=auxind(i)
435         auxi = nlayer-i+1
436         !O2 interpolated coefficient
437         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayer) *
438     $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
439         !CO2 interpolated coefficient
440         jfotsout(indexint,1,auxi) = jfotsout(indexint,1,nlayer) *
441     $        (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind))
442         !O3p interpolated coefficient
443         jfotsout(indexint,3,auxi) = jfotsout(indexint,3,nlayer) *
444     $        (wm(i)*auxjo3p(ind+1) + wp(i)*auxjo3p(ind))
445         !N2 interpolated coefficient
446         jfotsout(indexint,8,auxi) = jfotsout(indexint,8,nlayer) *
447     $        (wm(i)*auxjn2(ind+1) + wp(i)*auxjn2(ind))
448         !CO interpolated coefficient
449         jfotsout(indexint,11,auxi) =
450     $        jfotsout(indexint,11,nlayer) *
451     $        (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind))
452         !H interpolated coefficient
453         jfotsout(indexint,12,auxi) =
454     $        jfotsout(indexint,12,nlayer) *
455     $        (wm(i)*auxjh(ind+1) + wp(i)*auxjh(ind))
456      enddo
457      !Only if chemthermod.ge.2
458      if(chemthermod.ge.2) then
459         do i=1,nlayer
460            ind=auxind(i)
461            auxi = nlayer-i+1
462            !N interpolated coefficient
463            jfotsout(indexint,9,auxi) =
464     $           jfotsout(indexint,9,nlayer) *
465     $           (wm(i)*auxjn(ind+1) + wp(i)*auxjn(ind))
466            !NO interpolated coefficient
467            jfotsout(indexint,10,auxi) =
468     $           jfotsout(indexint,10,nlayer) *
469     $           (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind))
470            !NO2 interpolated coefficient
471            jfotsout(indexint,13,auxi) =
472     $           jfotsout(indexint,13,nlayer) *
473     $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
474         enddo
475      endif
476c     End interval 16
477
478
479ccccccccccccccccccccccccccccccc
480c     90.9-119.5nm (int 17-24)
481c
482c     Absorption by:
483c     CO2, O2, N2, NO, CO, NO2
484ccccccccccccccccccccccccccccccc
485
486c     Input column
487
488      do i=1,nlayer
489         auxcolinp(nlayer-i+1) = co2colx(i) + o2colx(i) + n2colx(i) +
490     $        nocolx(i) + cocolx(i) + no2colx(i)
491      end do
492
493      do indexint=17,24
494
495c     Interpolations
496
497         do i=1,nz2
498            auxi = nz2-i+1
499            !CO2 tabulated coefficient
500            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
501            !O2 tabulated coefficient
502            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
503            !N2 tabulated coefficient
504            auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
505            !CO tabulated coefficient
506            auxjco(i) = jabsifotsintpar(auxi,11,indexint)           
507            !Tabulated column
508            auxcoltab(i) = c17_24(auxi)
509         enddo
510         !Only if chemthermod.ge.2
511         if(chemthermod.ge.2) then
512            do i=1,nz2
513               auxi = nz2-i+1
514               !NO tabulated coefficient
515               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
516               !NO2 tabulated coefficient
517               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
518            enddo
519         endif
520
521         call interfast
522     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
523
524         !Correction to include T variation of CO2 cross section
525         if(indexint.eq.24) then
526            do i=1,nlayer
527               auxi = nlayer-i+1
528               if(sigma(indexint,auxi)*
529     $              alfa(indexint,auxi)*coltemp(auxi)
530     $              .lt.60.) then
531                  cortemp(i)=exp(-sigma(indexint,auxi)*
532     $                alfa(indexint,auxi)*coltemp(auxi))
533               else
534                  cortemp(i)=0.
535               end if
536            enddo
537         else
538            do i=1,nlayer
539               cortemp(i)=1.
540            enddo
541         end if
542         do i=1,nlayer           
543            ind=auxind(i)
544            auxi = nlayer-i+1
545            !O2 interpolated coefficient
546            jfotsout(indexint,2,auxi) =
547     $           jfotsout(indexint,2,nlayer) *
548     $           (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind)) *
549     $           cortemp(i)
550            !CO2 interpolated coefficient
551            jfotsout(indexint,1,auxi) =
552     $           jfotsout(indexint,1,nlayer) *
553     $           (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind))
554     $           * cortemp(i)
555            if(indexint.eq.24) jfotsout(indexint,1,auxi)=
556     $           jfotsout(indexint,1,auxi)*
557     $           (1+alfa(indexint,auxi)*
558     $           (t2(auxi)-t0(auxi)))
559            !N2 interpolated coefficient
560            jfotsout(indexint,8,auxi) =
561     $           jfotsout(indexint,8,nlayer) *
562     $           (wm(i)*auxjn2(ind+1) + wp(i)*auxjn2(ind)) *
563     $           cortemp(i)           
564            !CO interpolated coefficient
565            jfotsout(indexint,11,auxi) =
566     $           jfotsout(indexint,11,nlayer) *
567     $           (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind)) *
568     $           cortemp(i)           
569         enddo
570         !Only if chemthermod.ge.2
571         if(chemthermod.ge.2) then
572            do i=1,nlayer
573               ind=auxind(i)
574               auxi = nlayer-i+1
575               !NO interpolated coefficient
576               jfotsout(indexint,10,auxi)=
577     $              jfotsout(indexint,10,nlayer) *
578     $              (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind)) *
579     $              cortemp(i)
580               !NO2 interpolated coefficient
581               jfotsout(indexint,13,auxi)=
582     $              jfotsout(indexint,13,nlayer) *
583     $              (wm(i)*auxjno2(ind+1)+ wp(i)*auxjno2(ind)) *
584     $              cortemp(i)
585            enddo
586         endif               
587      end do
588c     End intervals 17-24
589
590
591ccccccccccccccccccccccccccccccc
592c     119.6-167.0nm (int 25-29)
593c
594c     Absorption by:
595c     CO2, O2, H2O, H2O2, NO,
596c     CO, NO2
597ccccccccccccccccccccccccccccccc
598
599c     Input atmospheric column
600
601      do i=1,nlayer
602         auxcolinp(nlayer-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
603     $        h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
604      end do
605
606      do indexint=25,29
607
608c     Interpolations
609
610         do i=1,nz2
611            auxi = nz2-i+1
612            !CO2 tabulated coefficient
613            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
614            !O2 tabulated coefficient
615            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
616            !H2O tabulated coefficient
617            auxjh2o(i) = jabsifotsintpar(auxi,4,indexint)
618            !H2O2 tabulated coefficient
619            auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)           
620            !CO tabulated coefficient
621            auxjco(i) = jabsifotsintpar(auxi,11,indexint)           
622            !Tabulated column
623            auxcoltab(i) = c25_29(auxi)
624         enddo
625         !Only if chemthermod.ge.2
626         if(chemthermod.ge.2) then
627            do i=1,nz2
628               auxi = nz2-i+1
629               !NO tabulated coefficient
630               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
631               !NO2 tabulated coefficient
632               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
633            enddo
634         endif
635         call interfast
636     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
637
638         do i=1,nlayer
639            ind=auxind(i)
640            auxi = nlayer-i+1
641            !Correction to include T variation of CO2 cross section
642            if(sigma(indexint,auxi)*alfa(indexint,auxi)*
643     $           coltemp(auxi).lt.60.) then
644               cortemp(i)=exp(-sigma(indexint,auxi)*
645     $              alfa(indexint,auxi)*coltemp(auxi))
646            else
647               cortemp(i)=0.
648            end if
649            !CO2 interpolated coefficient
650            jfotsout(indexint,1,auxi) =
651     $           jfotsout(indexint,1,nlayer) *
652     $           (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind)) *
653     $           cortemp(i) *
654     $           (1+alfa(indexint,auxi)*
655     $           (t2(auxi)-t0(auxi)))
656            !O2 interpolated coefficient
657            jfotsout(indexint,2,auxi) =
658     $           jfotsout(indexint,2,nlayer) *
659     $           (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind)) *
660     $           cortemp(i)
661            !H2O interpolated coefficient
662            jfotsout(indexint,4,auxi) =
663     $           jfotsout(indexint,4,nlayer) *
664     $           (wm(i)*auxjh2o(ind+1) + wp(i)*auxjh2o(ind)) *
665     $           cortemp(i)
666            !H2O2 interpolated coefficient
667            jfotsout(indexint,6,auxi) =
668     $           jfotsout(indexint,6,nlayer) *
669     $           (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind)) *
670     $           cortemp(i)           
671            !CO interpolated coefficient
672            jfotsout(indexint,11,auxi) =
673     $           jfotsout(indexint,11,nlayer) *
674     $           (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind)) *
675     $           cortemp(i)
676         enddo
677         !Only if chemthermod.ge.2
678         if(chemthermod.ge.2) then
679            do i=1,nlayer
680               ind=auxind(i)
681               auxi = nlayer-i+1
682               !NO interpolated coefficient
683               jfotsout(indexint,10,auxi)=
684     $              jfotsout(indexint,10,nlayer) *
685     $              (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind)) *
686     $              cortemp(i)
687               !NO2 interpolated coefficient
688               jfotsout(indexint,13,auxi)=
689     $              jfotsout(indexint,13,nlayer) *
690     $              (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind)) *
691     $              cortemp(i)
692            enddo
693         endif
694
695      end do
696
697c     End intervals 25-29
698
699
700cccccccccccccccccccccccccccccccc
701c     167.1-202.5nm (int 30-31)
702c   
703c     Absorption by:
704c     CO2, O2, H2O, H2O2, NO,
705c     NO2
706cccccccccccccccccccccccccccccccc
707
708c     Input atmospheric column
709
710      do i=1,nlayer
711         auxcolinp(nlayer-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
712     $        h2o2colx(i) + nocolx(i) + no2colx(i)
713      end do
714
715c     Interpolation
716
717      do indexint=30,31
718
719         do i=1,nz2
720            auxi = nz2-i+1
721            !CO2 tabulated coefficient
722            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
723            !O2 tabulated coefficient
724            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
725            !H2O tabulated coefficient
726            auxjh2o(i) = jabsifotsintpar(auxi,4,indexint)
727            !H2O2 tabulated coefficient
728            auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)           
729            !Tabulated column
730            auxcoltab(i) = c30_31(auxi)
731         enddo
732         !Only if chemthermod.ge.2
733         if(chemthermod.ge.2) then
734            do i=1,nz2
735               auxi = nz2-i+1
736               !NO tabulated coefficient
737               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
738               !NO2 tabulated coefficient
739               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
740            enddo
741         endif
742
743         call interfast
744     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
745
746         do i=1,nlayer
747            ind=auxind(i)
748            auxi = nlayer-i+1
749            !Correction to include T variation of CO2 cross section
750            if(sigma(indexint,auxi)*alfa(indexint,auxi)*
751     $           coltemp(auxi).lt.60.) then
752               cortemp(i)=exp(-sigma(indexint,auxi)*
753     $              alfa(indexint,auxi)*coltemp(auxi))
754            else
755               cortemp(i)=0.
756            end if
757            !CO2 interpolated coefficient
758            jfotsout(indexint,1,auxi) =
759     $           jfotsout(indexint,1,nlayer) *
760     $           (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind)) *
761     $           cortemp(i) *
762     $           (1+alfa(indexint,auxi)*
763     $           (t2(auxi)-t0(auxi)))
764            !O2 interpolated coefficient
765            jfotsout(indexint,2,auxi) =
766     $           jfotsout(indexint,2,nlayer) *
767     $           (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind)) *
768     $           cortemp(i)
769            !H2O interpolated coefficient
770            jfotsout(indexint,4,auxi) =
771     $           jfotsout(indexint,4,nlayer) *
772     $           (wm(i)*auxjh2o(ind+1) + wp(i)*auxjh2o(ind)) *
773     $           cortemp(i)
774            !H2O2 interpolated coefficient
775            jfotsout(indexint,6,auxi) =
776     $           jfotsout(indexint,6,nlayer) *
777     $           (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind)) *
778     $           cortemp(i)           
779         enddo
780         !Only if chemthermod.ge.2
781         if(chemthermod.ge.2) then
782            do i=1,nlayer
783               ind=auxind(i)
784               auxi = nlayer-i+1
785               !NO interpolated coefficient
786               jfotsout(indexint,10,auxi)=
787     $              jfotsout(indexint,10,nlayer) *
788     $              (wm(i)*auxjno(ind+1) +wp(i)*auxjno(ind)) *
789     $              cortemp(i)
790               !NO2 interpolated coefficient
791               jfotsout(indexint,13,auxi)=
792     $              jfotsout(indexint,13,nlayer) *
793     $              (wm(i)*auxjno2(ind+1)+wp(i)*auxjno2(ind)) *
794     $              cortemp(i)
795            enddo
796         endif
797
798      end do
799
800c     End intervals 30-31
801
802
803ccccccccccccccccccccccccccccccc
804c     202.6-210.0nm (int 32)
805c
806c     Absorption by:
807c     CO2, O2, H2O2, NO, NO2
808ccccccccccccccccccccccccccccccc
809
810c     Input atmospheric column
811
812      indexint=32
813      do i=1,nlayer
814         auxcolinp(nlayer-i+1) =co2colx(i) + o2colx(i) + h2o2colx(i) +
815     $        nocolx(i) + no2colx(i)
816      end do
817
818c     Interpolation
819
820      do i=1,nz2
821         auxi = nz2-i+1
822         !CO2 tabulated coefficient
823         auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
824         !O2 tabulated coefficient
825         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
826         !H2O2 tabulated coefficient
827         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)         
828         !Tabulated column
829         auxcoltab(i) = c32(auxi)
830      enddo
831      !Only if chemthermod.ge.2
832      if(chemthermod.ge.2) then
833         do i=1,nz2
834            auxi = nz2-i+1
835            !NO tabulated coefficient
836            auxjno(i) = jabsifotsintpar(auxi,10,indexint)
837            !NO2 tabulated coefficient
838            auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
839         enddo
840      endif
841      call interfast
842     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
843
844      do i=1,nlayer
845         ind=auxind(i)
846         auxi = nlayer-i+1
847         !Correction to include T variation of CO2 cross section
848         if(sigma(indexint,nlayer-i+1)*alfa(indexint,auxi)*
849     $        coltemp(auxi).lt.60.) then
850            cortemp(i)=exp(-sigma(indexint,auxi)*
851     $           alfa(indexint,auxi)*coltemp(auxi))
852         else
853            cortemp(i)=0.
854         end if
855         !CO2 interpolated coefficient
856         jfotsout(indexint,1,auxi) =
857     $        jfotsout(indexint,1,nlayer) *
858     $        (wm(i)*auxjco2(ind+1)+wp(i)*auxjco2(ind)) *
859     $        cortemp(i) *
860     $        (1+alfa(indexint,auxi)*
861     $        (t2(auxi)-t0(auxi)))
862         !O2 interpolated coefficient
863         jfotsout(indexint,2,auxi) =
864     $        jfotsout(indexint,2,nlayer) *
865     $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind)) *
866     $        cortemp(i)
867         !H2O2 interpolated coefficient
868         jfotsout(indexint,6,auxi) =
869     $        jfotsout(indexint,6,nlayer) *
870     $        (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind)) *
871     $        cortemp(i)         
872      enddo
873      !Only if chemthermod.ge.2
874      if(chemthermod.ge.2) then
875         do i=1,nlayer
876            auxi = nlayer-i+1
877            ind=auxind(i)
878            !NO interpolated coefficient
879            jfotsout(indexint,10,auxi) =
880     $           jfotsout(indexint,10,nlayer) *
881     $           (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind)) *
882     $           cortemp(i)
883           !NO2 interpolated coefficient
884            jfotsout(indexint,13,auxi) =
885     $           jfotsout(indexint,13,nlayer) *
886     $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind)) *
887     $           cortemp(i)
888         enddo
889      endif
890
891c     End of interval 32
892
893
894ccccccccccccccccccccccccccccccc
895c     210.1-231.0nm (int 33)
896c     
897c     Absorption by:
898c     O2, H2O2, NO2
899ccccccccccccccccccccccccccccccc
900
901c     Input atmospheric column
902
903      indexint=33
904      do i=1,nlayer
905         auxcolinp(nlayer-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
906      end do
907
908c     Interpolation
909
910      do i=1,nz2
911         auxi = nz2-i+1
912         !O2 tabulated coefficient
913         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
914         !H2O2 tabulated coefficient
915         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
916         !Tabulated column
917         auxcoltab(i) = c33(auxi)
918      enddo
919      !Only if chemthermod.ge.2
920      if(chemthermod.ge.2) then
921         do i=1,nz2
922            !NO2 tabulated coefficient
923            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
924         enddo
925      endif
926      call interfast
927     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
928
929      do i=1,nlayer
930         ind=auxind(i)
931         auxi = nlayer-i+1
932         !O2 interpolated coefficient
933         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayer) *
934     $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
935         !H2O2 interpolated coefficient
936         jfotsout(indexint,6,auxi) = jfotsout(indexint,6,nlayer) *
937     $        (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind))
938      enddo
939      !Only if chemthermod.ge.2
940      if(chemthermod.ge.2) then
941         do i=1,nlayer
942            ind=auxind(i)
943            !NO2 interpolated coefficient
944            jfotsout(indexint,13,nlayer-i+1) =
945     $           jfotsout(indexint,13,nlayer) *
946     $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
947         enddo
948      endif
949
950c     End of interval 33
951
952
953ccccccccccccccccccccccccccccccc
954c     231.1-240.0nm (int 34)
955c
956c     Absorption by:
957c     O2, H2O2, O3, NO2
958ccccccccccccccccccccccccccccccc
959
960c     Input atmospheric column
961
962      indexint=34
963      do i=1,nlayer
964         auxcolinp(nlayer-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) +
965     $        no2colx(i)
966      end do
967
968c     Interpolation
969
970      do i=1,nz2
971         auxi = nz2-i+1
972         !O2 tabulated coefficient
973         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
974         !H2O2 tabulated coefficient
975         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
976         !O3 tabulated coefficient
977         auxjo3(i) = jabsifotsintpar(auxi,7,indexint)         
978         !Tabulated column
979         auxcoltab(i) = c34(nz2-i+1)
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,nlayer,auxcoltab,nz2,limdown,limup)
990
991      do i=1,nlayer
992         ind=auxind(i)
993         auxi = nlayer-i+1
994         !O2 interpolated coefficient
995         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayer) *
996     $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
997         !H2O2 interpolated coefficient
998         jfotsout(indexint,6,auxi) = jfotsout(indexint,6,nlayer) *
999     $        (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind))
1000         !O3 interpolated coefficient
1001         jfotsout(indexint,7,auxi) = jfotsout(indexint,7,nlayer) *
1002     $        (wm(i)*auxjo3(ind+1) + wp(i)*auxjo3(ind))
1003      enddo
1004      !Only if chemthermod.ge.2
1005      if(chemthermod.ge.2) then
1006         do i=1,nlayer
1007            ind=auxind(i)
1008            !NO2 interpolated coefficient
1009            jfotsout(indexint,13,nlayer-i+1) =
1010     $           jfotsout(indexint,13,nlayer) *
1011     $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
1012         enddo
1013      endif
1014
1015c     End of interval 34     
1016
1017
1018ccccccccccccccccccccccccccccccc
1019c     240.1-337.7nm (int 35)
1020c
1021c     Absorption by:
1022c     H2O2, O3, NO2
1023ccccccccccccccccccccccccccccccc
1024
1025c     Input atmospheric column
1026
1027      indexint=35
1028      do i=1,nlayer
1029         auxcolinp(nlayer-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
1030      end do
1031
1032c     Interpolation
1033
1034      do i=1,nz2
1035         auxi = nz2-i+1
1036         !H2O2 tabulated coefficient
1037         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
1038         !O3 tabulated coefficient
1039         auxjo3(i) = jabsifotsintpar(auxi,7,indexint)
1040         !Tabulated column
1041         auxcoltab(i) = c35(auxi)
1042      enddo
1043      !Only if chemthermod.ge.2
1044      if(chemthermod.ge.2) then
1045         do i=1,nz2
1046            !NO2 tabulated coefficient
1047            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
1048         enddo
1049      endif
1050      call interfast
1051     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
1052
1053      do i=1,nlayer
1054         ind=auxind(i)
1055         auxi = nlayer-i+1
1056         !H2O2 interpolated coefficient
1057         jfotsout(indexint,6,auxi) = jfotsout(indexint,6,nlayer) *
1058     $        (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind))
1059         !O3 interpolated coefficient
1060         jfotsout(indexint,7,auxi) = jfotsout(indexint,7,nlayer) *
1061     $        (wm(i)*auxjo3(ind+1) + wp(i)*auxjo3(ind))
1062      enddo
1063      if(chemthermod.ge.2) then
1064         do i=1,nlayer
1065            ind=auxind(i)
1066            !NO2 interpolated coefficient
1067            jfotsout(indexint,13,nlayer-i+1) =
1068     $           jfotsout(indexint,13,nlayer) *
1069     $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
1070         enddo
1071      endif
1072
1073c     End of interval 35
1074
1075ccccccccccccccccccccccccccccccc
1076c     337.8-800.0 nm (int 36)
1077c     
1078c     Absorption by:
1079c     O3, NO2
1080ccccccccccccccccccccccccccccccc
1081
1082c     Input atmospheric column
1083
1084      indexint=36
1085      do i=1,nlayer
1086         auxcolinp(nlayer-i+1) = o3colx(i) + no2colx(i)
1087      end do
1088
1089c     Interpolation
1090
1091      do i=1,nz2
1092         auxi = nz2-i+1
1093         !O3 tabulated coefficient
1094         auxjo3(i) = jabsifotsintpar(auxi,7,indexint)         
1095         !Tabulated column
1096         auxcoltab(i) = c36(auxi)
1097      enddo
1098      !Only if chemthermod.ge.2
1099      if(chemthermod.ge.2) then
1100         do i=1,nz2
1101            !NO2 tabulated coefficient
1102            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
1103         enddo
1104      endif
1105      call interfast
1106     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
1107
1108      do i=1,nlayer
1109         ind=auxind(i)
1110         !O3 interpolated coefficient
1111         jfotsout(indexint,7,nlayer-i+1) =
1112     $        jfotsout(indexint,7,nlayer) *
1113     $        (wm(i)*auxjo3(ind+1) + wp(i)*auxjo3(ind))
1114      enddo
1115      !Only if chemthermod.ge.2
1116      if(chemthermod.ge.2) then
1117         do i=1,nlayer
1118            ind=auxind(i)
1119            !NO2 interpolated coefficient
1120            jfotsout(indexint,13,nlayer-i+1) =
1121     $           jfotsout(indexint,13,nlayer) *
1122     $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
1123         enddo
1124      endif
1125
1126c     End of interval 36
1127
1128c     End of interpolation to obtain photoabsorption rates
1129
1130c     Correction to the actual Sun-Mars distance
1131
1132      jfotsout(:,:,:)=jfotsout(:,:,:)*(1.52/dist_sol)**2
1133
1134      end
1135
1136
1137
1138
1139
1140
Note: See TracBrowser for help on using the repository browser.