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

Last change on this file since 1571 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

  • Remove all reference/use of nlayermx and dimphys.h
  • Made use of automatic arrays whenever arrays are needed with dimension nlayer
  • Remove lots of obsolete reference to dimensions.h
  • Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

  • Deleted param.h and put contents into module param_v4_h
  • Adapted testphys1d, newstart, etc...
  • Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

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