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

Last change on this file since 3567 was 3464, checked in by emillour, 2 months ago

Mars PCM:
Some tidying in aeronomars:

  • make a jthermcalc_util.F to contain utility routines (used by jthermcal & jthermcalc_e107). Also add the possibility (turned off by default) in the interfast routine to do extra sanity checks.
  • turn chemthermos, euvheat, hrtherm, jthermcalc, jthermcalc_e107, paramphoto_compact and photochemistry into modules.

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