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

Last change on this file since 706 was 705, checked in by emillour, 13 years ago

Mars GCM:

  • Added possibility to run with a varying EUV cycle following real one. The flag solvarmod=1 triggers this behaviour, with companion flag solvaryear=## , where ## is the Mars Year (from 23 to 30). Setting solvarmod=0 reverts to 'old' behaviour, where there is a constant EUV forcing throughout the run, set by the "solarcondate" flag.
  • Needs corresponding input data files ("param_v6" subdirectory of "EUV" subdirectory in "datadir").
  • Added files jthermcalc_e107.F and param_read_e107.F in "aeronomars", modified files euvheat.F90, hrtherm.F, chemthermos.F90, param_v4.h and param_read.F in "aeronomars" and files inifis.F, physiq.F and callkeys.h in "phymars".

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