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

Last change on this file since 2325 was 1684, checked in by emillour, 8 years ago

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