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

Last change on this file since 2512 was 2429, checked in by emillour, 5 years ago

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