source: trunk/LMDZ.VENUS/libf/phyvenus/jthermcalc_e107.F @ 2203

Last change on this file since 2203 was 1530, checked in by emillour, 9 years ago

Venus and Titan GCMs:
Updates in the physics to keep up with updates in LMDZ5 (up to
LMDZ5 trunk, rev 2350) concerning dynamics/physics separation:

  • Adapted makelmdz and makelmdz_fcm script to stop if trying to compile 1d model or newstart or start2archive in parallel.
  • got rid of references to "dimensions.h" in physics. Within physics packages, use nbp_lon (=iim), nbp_lat (=jjmp1) and nbp_lev (=llm) from module mod_grid_phy_lmdz (in phy_common) instead. Only partially done for Titan, because of many hard-coded commons; a necessary first step will be to clean these up (using modules).

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