source: trunk/LMDZ.MARS/libf/aeronomars/jthermcalc.F

Last change on this file was 3464, checked in by emillour, 5 weeks 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
Line 
1      module jthermcalc_mod
2     
3      implicit none
4     
5      contains
6
7c**********************************************************************
8
9      subroutine jthermcalc(ig,nlayer,chemthermod,
10     .      rm,nesptherm,tx,iz,zenit)
11
12
13c     feb 2002        fgg           first version
14c     nov 2002        fgg           second version
15c
16c modified from paramhr.F
17c MAC July 2003
18c**********************************************************************
19
20      use param_v4_h, only: jfotsout,crscabsi2,
21     .    c1_16,c17_24,c25_29,c30_31,c32,c33,c34,c35,c36,
22     .    co2crsc195,co2crsc295,t0,
23     .    jabsifotsintpar,ninter,nz2
24      use jthermcalc_util, only: column, interfast
25     
26      implicit none
27
28c     input and output variables
29
30      integer    ig,nlayer
31      integer    chemthermod
32      integer    nesptherm                      !Number of species considered
33      real       rm(nlayer,nesptherm)         !Densities (cm-3)
34      real       tx(nlayer)                   !temperature
35      real       zenit                          !SZA
36      real       iz(nlayer)                   !Local altitude
37
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     
59      integer    i,j,k,indexint                 !indexes
60      character  dn
61
62
63
64c     variables used in interpolation
65
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)
80      real*8      wp(nlayer),wm(nlayer)
81      real*8      auxcolinp(nlayer)
82      integer     auxind(nlayer)
83      integer     auxi
84      integer     ind
85      real*8      cortemp(nlayer)
86
87      real*8      limdown                      !limits for interpolation
88      real*8      limup                        !        ""
89
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
94
95c*************************PROGRAM STARTS*******************************
96     
97      if(zenit.gt.140.) then
98         dn='n'
99         else
100         dn='d'
101      end if
102      if(dn.eq.'n') then
103        return
104      endif
105     
106      !Initializing the photoabsorption coefficients
107      jfotsout(:,:,:)=0.
108
109      !Auxiliar temperature to take into account the temperature dependence
110      !of CO2 cross section
111      do i=1,nlayer
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
117      !Calculation of column amounts
118      call column(ig,nlayer,chemthermod,rm,nesptherm,tx,iz,zenit,
119     $     co2colx,o2colx,o3pcolx,h2colx,h2ocolx,
120     $     h2o2colx,o3colx,n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
121
122      !Auxiliar column to include the temperature dependence
123      !of CO2 cross section
124      coltemp(nlayer)=co2colx(nlayer)*abs(t2(nlayer)-t0(nlayer))
125      do i=nlayer-1,1,-1
126        coltemp(i)=!coltemp(i+1)+     PQ SE ELIMINA? REVISAR
127     $         ( rm(i,i_co2) + rm(i+1,i_co2) ) * 0.5
128     $         * 1e5 * (iz(i+1)-iz(i)) * abs(t2(i)-t0(i))
129      end do
130     
131      !Calculation of CO2 cross section at temperature t0(i)
132      do i=1,nlayer
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
139
140! Interpolation to the tabulated photoabsorption coefficients for each species
141! in each spectral interval
142
143
144c     auxcolinp-> Actual atmospheric column
145c     auxj*-> Tabulated photoabsorption coefficients
146c     auxcoltab-> Tabulated atmospheric columns
147
148ccccccccccccccccccccccccccccccc
149c     0.1,5.0 (int 1)
150c
151c     Absorption by:
152c     CO2, O2, O, H2, N
153ccccccccccccccccccccccccccccccc
154
155c     Input atmospheric column
156      indexint=1
157      do i=1,nlayer
158         auxcolinp(nlayer-i+1) = co2colx(i)*crscabsi2(1,indexint) +
159     $        o2colx(i)*crscabsi2(2,indexint) +
160     $        o3pcolx(i)*crscabsi2(3,indexint) +
161     $        h2colx(i)*crscabsi2(5,indexint) +
162     $        ncolx(i)*crscabsi2(9,indexint)
163         
164         
165      end do
166      limdown=1.e-20
167      limup=1.e26
168
169c     Interpolations
170
171      do i=1,nz2
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)
183      enddo
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
191
192      call interfast
193     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
194      do i=1,nlayer
195         ind=auxind(i)
196         auxi=nlayer-i+1
197         ! Ehouarn: test
198          if ((ind+1.gt.nz2).or.(ind.le.0)) then
199            write(*,*) "jthercalc error: ind=",ind,ig,zenit
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
207            write(*,*) " auxcoltab(1:nz2)=",auxcoltab
208            write(*,*) " limdown=",limdown
209            write(*,*) " limup=",limup
210#ifndef MESOSCALE
211            call abort_physic('jthermcalc','error',1)
212#endif
213          endif
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)
226      enddo
227      !Only if chemthermod.ge.2
228      !N interpolated coefficient
229      if(chemthermod.ge.2) then
230         do i=1,nlayer
231            ind=auxind(i)
232            jfotsout(indexint,9,nlayer-i+1) =  wm(i)*auxjn(ind+1) +
233     $         wp(i)*auxjn(ind)
234         enddo
235      endif
236         
237
238c     End interval 1
239
240
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
248
249c     Input atmospheric column
250      do indexint=2,15
251         do i=1,nlayer
252            auxcolinp(nlayer-i+1) = co2colx(i)*crscabsi2(1,indexint)+
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)
262         end do
263
264c     Interpolations
265
266         do i=1,nz2
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)
284         enddo
285         !Only if chemthermod.ge.2
286         if(chemthermod.ge.2) then
287            do i=1,nz2
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)
295            enddo
296         endif
297
298          call interfast(wm,wp,auxind,auxcolinp,nlayer,
299     $        auxcoltab,nz2,limdown,limup)
300          do i=1,nlayer
301             ind=auxind(i)
302             auxi = nlayer-i+1
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) +
323     $            wp(i)*auxjh(ind)
324          enddo
325          !Only if chemthermod.ge.2
326          if(chemthermod.ge.2) then
327             do i=1,nlayer
328                ind=auxind(i)
329                auxi = nlayer-i+1
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   
341      end do
342
343c     End intervals 2-15
344
345
346ccccccccccccccccccccccccccccccc
347c     80.6-90.8nm (int16)
348c
349c     Absorption by:
350c     CO2, O2, O, N2, N, NO,
351c     CO, H, NO2
352ccccccccccccccccccccccccccccccc
353
354c     Input atmospheric column
355      indexint=16
356      do i=1,nlayer
357         auxcolinp(nlayer-i+1) = co2colx(i)*crscabsi2(1,indexint)+
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
368c     Interpolations
369
370      do i=1,nz2
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)
388      enddo
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
401
402      call interfast
403     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
404      do i=1,nlayer
405         ind=auxind(i)
406         auxi = nlayer-i+1
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)         
425      enddo
426      !Only if chemthermod.ge.2
427      if(chemthermod.ge.2) then
428         do i=1,nlayer
429            ind=auxind(i)
430            auxi = nlayer-i+1
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)
440         enddo
441      endif
442c     End interval 16
443
444
445ccccccccccccccccccccccccccccccc
446c     90.9-119.5nm (int 17-24)
447c
448c     Absorption by:
449c     CO2, O2, N2, NO, CO, NO2
450ccccccccccccccccccccccccccccccc
451
452c     Input column
453
454      do i=1,nlayer
455         auxcolinp(nlayer-i+1) = co2colx(i) + o2colx(i) + n2colx(i) +
456     $        nocolx(i) + cocolx(i) + no2colx(i)
457      end do
458
459      do indexint=17,24
460
461c     Interpolations
462
463         do i=1,nz2
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)
475         enddo
476         !Only if chemthermod.ge.2
477         if(chemthermod.ge.2) then
478            do i=1,nz2
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)
484            enddo
485         endif
486
487         call interfast
488     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
489         !Correction to include T variation of CO2 cross section
490         if(indexint.eq.24) then
491            do i=1,nlayer
492               auxi = nlayer-i+1
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.
500               end if
501            enddo
502         else
503            do i=1,nlayer
504               cortemp(i)=1.
505            enddo
506         end if
507         do i=1,nlayer           
508            ind=auxind(i)
509            auxi = nlayer-i+1
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
529            do i=1,nlayer
530               ind=auxind(i)
531               auxi = nlayer-i+1
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)
538            enddo
539         endif               
540      end do
541c     End intervals 17-24
542
543
544ccccccccccccccccccccccccccccccc
545c     119.6-167.0nm (int 25-29)
546c
547c     Absorption by:
548c     CO2, O2, H2O, H2O2, NO,
549c     CO, NO2
550ccccccccccccccccccccccccccccccc
551
552c     Input atmospheric column
553
554      do i=1,nlayer
555         auxcolinp(nlayer-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
556     $        h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
557      end do
558
559      do indexint=25,29
560
561c     Interpolations
562
563         do i=1,nz2
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)
577         enddo
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
589     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
590         do i=1,nlayer
591            ind=auxind(i)
592            auxi = nlayer-i+1
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.
600            end if
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)
618         enddo
619         !Only if chemthermod.ge.2
620         if(chemthermod.ge.2) then
621            do i=1,nlayer
622               ind=auxind(i)
623               auxi = nlayer-i+1
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)
630            enddo
631         endif
632
633      end do
634
635c     End intervals 25-29
636
637
638cccccccccccccccccccccccccccccccc
639c     167.1-202.5nm (int 30-31)
640c   
641c     Absorption by:
642c     CO2, O2, H2O, H2O2, NO,
643c     NO2
644cccccccccccccccccccccccccccccccc
645
646c     Input atmospheric column
647
648      do i=1,nlayer
649         auxcolinp(nlayer-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
650     $        h2o2colx(i) + nocolx(i) + no2colx(i)
651      end do
652
653c     Interpolation
654
655      do indexint=30,31
656
657         do i=1,nz2
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)
669         enddo
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
680
681         call interfast
682     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
683         do i=1,nlayer
684            ind=auxind(i)
685            auxi = nlayer-i+1
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.
693            end if
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)           
708         enddo
709         !Only if chemthermod.ge.2
710         if(chemthermod.ge.2) then
711            do i=1,nlayer
712               ind=auxind(i)
713               auxi = nlayer-i+1
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)
720            enddo
721         endif
722
723      end do
724
725c     End intervals 30-31
726
727
728ccccccccccccccccccccccccccccccc
729c     202.6-210.0nm (int 32)
730c
731c     Absorption by:
732c     CO2, O2, H2O2, NO, NO2
733ccccccccccccccccccccccccccccccc
734
735c     Input atmospheric column
736
737      indexint=32
738      do i=1,nlayer
739         auxcolinp(nlayer-i+1) =co2colx(i) + o2colx(i) + h2o2colx(i) +
740     $        nocolx(i) + no2colx(i)
741      end do
742
743c     Interpolation
744
745      do i=1,nz2
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)
755      enddo
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
767     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
768      do i=1,nlayer
769         ind=auxind(i)
770         auxi = nlayer-i+1
771         !Correction to include T variation of CO2 cross section
772         if(sigma(indexint,nlayer-i+1)*alfa(indexint,auxi)*
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.
778         end if
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)         
790      enddo
791      !Only if chemthermod.ge.2
792      if(chemthermod.ge.2) then
793         do i=1,nlayer
794            auxi = nlayer-i+1
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)
802         enddo
803      endif
804
805c     End of interval 32
806
807
808ccccccccccccccccccccccccccccccc
809c     210.1-231.0nm (int 33)
810c     
811c     Absorption by:
812c     O2, H2O2, NO2
813ccccccccccccccccccccccccccccccc
814
815c     Input atmospheric column
816
817      indexint=33
818      do i=1,nlayer
819         auxcolinp(nlayer-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
820      end do
821
822c     Interpolation
823
824      do i=1,nz2
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)
832      enddo
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
841     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
842      do i=1,nlayer
843         ind=auxind(i)
844         auxi = nlayer-i+1
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)         
851      enddo
852      !Only if chemthermod.ge.2
853      if(chemthermod.ge.2) then
854         do i=1,nlayer
855            ind=auxind(i)
856            !NO2 interpolated coefficient
857            jfotsout(indexint,13,nlayer-i+1) = wm(i)*auxjno2(ind+1) +
858     $           wp(i)*auxjno2(ind)
859         enddo
860      endif
861
862c     End of interval 33
863
864
865ccccccccccccccccccccccccccccccc
866c     231.1-240.0nm (int 34)
867c
868c     Absorption by:
869c     O2, H2O2, O3, NO2
870ccccccccccccccccccccccccccccccc
871
872c     Input atmospheric column
873
874      indexint=34
875      do i=1,nlayer
876         auxcolinp(nlayer-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) +
877     $        no2colx(i)
878      end do
879
880c     Interpolation
881
882      do i=1,nz2
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)
892      enddo
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
901     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
902      do i=1,nlayer
903         ind=auxind(i)
904         auxi = nlayer-i+1
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)         
914      enddo
915      !Only if chemthermod.ge.2
916      if(chemthermod.ge.2) then
917         do i=1,nlayer
918            ind=auxind(i)
919            !NO2 interpolated coefficient
920            jfotsout(indexint,13,nlayer-i+1) = wm(i)*auxjno2(ind+1) +
921     $           wp(i)*auxjno2(ind)
922         enddo
923      endif
924
925c     End of interval 34     
926
927
928ccccccccccccccccccccccccccccccc
929c     240.1-337.7nm (int 35)
930c
931c     Absorption by:
932c     H2O2, O3, NO2
933ccccccccccccccccccccccccccccccc
934
935c     Input atmospheric column
936
937      indexint=35
938      do i=1,nlayer
939         auxcolinp(nlayer-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
940      end do
941
942c     Interpolation
943
944      do i=1,nz2
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)
952      enddo
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
961     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
962      do i=1,nlayer
963         ind=auxind(i)
964         auxi = nlayer-i+1
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)         
971      enddo
972      if(chemthermod.ge.2) then
973         do i=1,nlayer
974            ind=auxind(i)
975            !NO2 interpolated coefficient
976            jfotsout(indexint,13,nlayer-i+1) = wm(i)*auxjno2(ind+1) +
977     $           wp(i)*auxjno2(ind)
978         enddo
979      endif
980
981c     End of interval 35
982
983ccccccccccccccccccccccccccccccc
984c     337.8-800.0 nm (int 36)
985c     
986c     Absorption by:
987c     O3, NO2
988ccccccccccccccccccccccccccccccc
989
990c     Input atmospheric column
991
992      indexint=36
993      do i=1,nlayer
994         auxcolinp(nlayer-i+1) = o3colx(i) + no2colx(i)
995      end do
996
997c     Interpolation
998
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
1008         do i=1,nz2
1009            !NO2 tabulated coefficient
1010            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
1011         enddo
1012      endif
1013      call interfast
1014     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
1015      do i=1,nlayer
1016         ind=auxind(i)
1017         !O3 interpolated coefficient
1018         jfotsout(indexint,7,nlayer-i+1) = wm(i)*auxjo3(ind+1) +
1019     $        wp(i)*auxjo3(ind)         
1020      enddo
1021      !Only if chemthermod.ge.2
1022      if(chemthermod.ge.2) then
1023         do i=1,nlayer
1024            ind=auxind(i)
1025            !NO2 interpolated coefficient
1026            jfotsout(indexint,13,nlayer-i+1) = wm(i)*auxjno2(ind+1) +
1027     $           wp(i)*auxjno2(ind)
1028         enddo
1029      endif
1030
1031c     End of interval 36
1032
1033c     End of interpolation to obtain photoabsorption rates
1034
1035
1036      end subroutine jthermcalc
1037     
1038      end module jthermcalc_mod
1039
1040
1041
Note: See TracBrowser for help on using the repository browser.