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