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

Last change on this file since 1242 was 762, checked in by emillour, 12 years ago

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