source: trunk/LMDZ.VENUS/libf/phyvenus/jthermcalc.F @ 1543

Last change on this file since 1543 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: 58.9 KB
Line 
1c**********************************************************************
2
3      subroutine jthermcalc(ig,chemthermod,rm,nesptherm,tx,iz,zenit)
4
5
6c     feb 2002        fgg           first version
7c     nov 2002        fgg           second version
8c     
9c     mar 2014        gg            update for Venus GCM
10c
11c modified from paramhr.F
12c MAC July 2003
13c**********************************************************************
14      use dimphy
15      use conc
16c      use chemparam_mod
17      implicit none
18
19c     common variables and constants
20      include "param.h"
21      include "param_v4.h"
22
23c     input and output variables
24
25      integer    ig
26      integer    chemthermod
27      integer    nesptherm                  !Number of species considered
28      real       rm(klev,nesptherm)         !Densities (cm-3)
29      real       tx(klev)                   !temperature
30      real       zenit                      !SZA
31      real       iz(klev)                   !Local altitude
32
33
34c    local parameters and variables
35
36      real       co2colx(klev)              !column density of CO2 (cm^-2)
37      real       o3pcolx(klev)              !column density of O(3P)(cm^-2)
38      real       n2colx(klev)               !N2 column density (cm-2)
39      real       cocolx(klev)               !CO column density (cm-2)
40c      real       o2colx(klev)               !column density of O2(cm^-2)
41c      real       h2colx(klev)               !H2 column density (cm-2)
42c      real       h2ocolx(klev)              !H2O column density (cm-2)
43c      real       h2o2colx(klev)             !column density of H2O2(cm^-2)
44c      real       o3colx(klev)               !O3 column density (cm-2)
45c     real       hcolx(klev)                !H column density (cm-2)
46c     real       no2colx(klev)              !NO2 column density (cm-2)
47c     real       nocolx(klev)               !NO column density (cm-2)
48      real       t2(klev)
49      real       coltemp(klev)
50      real       sigma(ninter,klev)
51      real       alfa(ninter,klev)
52   
53      integer    i,j,k,indexint                 !indexes
54      character  dn
55
56
57
58c     variables used in interpolation
59
60      real*8      auxcoltab(nz2)
61      real*8      auxjco2(nz2)
62c     real*8      auxjo2(nz2)
63      real*8      auxjo3p(nz2)
64c     real*8      auxjh2o(nz2)
65c     real*8      auxjh2(nz2)
66c     real*8      auxjh2o2(nz2)
67c     real*8      auxjo3(nz2)
68      real*8      auxjn2(nz2)
69c      real*8      auxjn(nz2)
70c     real*8      auxjno(nz2)
71      real*8      auxjco(nz2)
72c     real*8      auxjh(nz2)
73c     real*8      auxjno2(nz2)
74      real*8      wp(klev),wm(klev)
75      real*8      auxcolinp(klev)
76      integer     auxind(klev)
77      integer     auxi
78      integer     ind
79      real*8      cortemp(klev)
80
81      real*8      limdown                      !limits for interpolation
82      real*8      limup                        !        ""
83
84 
85! Tracer indexes in the thermospheric chemistry:
86!!! ATTENTION. These values have to be identical to those in euvheat.F90
87!!! If the values are changed there, the same has to be done here  !!!
88
89      integer,parameter :: ix_co2=1
90      integer,parameter :: ix_n2=13
91c      integer,parameter :: i_n=14
92      integer,parameter :: ix_o=3
93      integer,parameter :: ix_co=4
94
95
96c*************************PROGRAM STARTS*******************************
97     
98      if(zenit.gt.140.) then
99         dn='n'
100         else
101         dn='d'
102      end if
103      if(dn.eq.'n') then
104        return
105      endif
106     
107      !Initializing the photoabsorption coefficients
108      jfotsout(:,:,:)=0.
109
110      !Auxiliar temperature to take into account the temperature dependence
111      !of CO2 cross section
112      do i=1,klev
113         t2(i)=tx(i)
114         if(t2(i).lt.195.0) t2(i)=195.0
115         if(t2(i).gt.295.0) t2(i)=295.0
116      end do
117
118      !Calculation of column amounts
119      call column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
120     $     co2colx,o3pcolx, n2colx,cocolx)
121
122      !Auxiliar column to include the temperature dependence
123      !of CO2 cross section
124      coltemp(klev)=co2colx(klev)*abs(t2(klev)-t0(klev))
125      do i=klev-1,1,-1
126        coltemp(i)=!coltemp(i+1)+     PQ SE ELIMINA? REVISAR
127     $         ( rm(i,ix_co2) + rm(i+1,ix_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,klev
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,klev
158         auxcolinp(klev-i+1) = co2colx(i)*crscabsi2(1,indexint) +
159c     $        o2colx(i)*crscabsi2(2,indexint) +
160     $        o3pcolx(i)*crscabsi2(3,indexint)
161c     $        h2colx(i)*crscabsi2(5,indexint) +
162      end do
163      limdown=1.e-20
164      limup=1.e26
165
166
167c     Interpolations
168
169      do i=1,nz2
170         auxi = nz2-i+1
171         !CO2 tabulated coefficient
172         auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
173         !O2 tabulated coefficient
174c         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
175         !O3p tabulated coefficient
176         auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
177         !H2 tabulated coefficient
178c         auxjh2(i) = jabsifotsintpar(auxi,5,indexint)
179         !N tabulated coefficient
180c         auxjn(i) = jabsifotsintpar(auxi,9,indexint)
181         !Tabulated column
182         auxcoltab(i) = c1_16(auxi,indexint)
183      enddo
184
185      !Only if chemthermod.ge.2
186      !N tabulated coefficient
187c      if(chemthermod.ge.2) then
188c         do i=1,nz2
189c            auxjn(i) = jabsifotsintpar(nz2-i+1,9,indexint)
190c         enddo
191c      endif
192
193      call interfast
194     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
195      do i=1,klev
196         ind=auxind(i)
197         auxi=klev-i+1
198         !CO2 interpolated coefficient
199         jfotsout(indexint,1,auxi) = wm(i)*auxjco2(ind+1) +
200     $        wp(i)*auxjco2(ind)
201         !O2 interpolated coefficient
202c         jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
203c     $        wp(i)*auxjo2(ind)
204         !O3p interpolated coefficient
205          jfotsout(indexint,3,auxi) = wm(i)*auxjo3p(ind+1) +
206     $        wp(i)*auxjo3p(ind)
207          !H2 interpolated coefficient
208c          jfotsout(indexint,5,auxi) = wm(i)*auxjh2(ind+1) +
209c     $         wp(i)*auxjh2(ind)
210c      !N interpolated coefficient
211c            jfotsout(indexint,9,auxi) =  wm(i)*auxjn(ind+1) +
212c     $         wp(i)*auxjn(ind)
213
214C      print*, '--- L214 jthermcal.F ---'
215C      print*, jfotsout(indexint,1,auxi)
216c      STOP         
217
218      enddo
219
220      !Only if chemthermod.ge.2
221      !N interpolated coefficient
222c      if(chemthermod.ge.2) then
223c         do i=1,klev
224c            ind=auxind(i)
225c            jfotsout(indexint,9,klev-i+1) =  wm(i)*auxjn(ind+1) +
226c     $         wp(i)*auxjn(ind)
227c         enddo
228c      endif
229
230       
231c     End interval 1
232
233
234ccccccccccccccccccccccccccccccc
235c     5-80.5nm (int 2-15)
236c
237c     Absorption by:
238c     CO2, O2, O, H2, N2, N,
239c     NO, CO, H, NO2
240ccccccccccccccccccccccccccccccc
241
242c     Input atmospheric column
243      do indexint=2,15
244         do i=1,klev
245            auxcolinp(klev-i+1) = co2colx(i)*crscabsi2(1,indexint)+
246     $           o3pcolx(i)*crscabsi2(3,indexint)+
247     $           n2colx(i)*crscabsi2(8,indexint)+
248     $           cocolx(i)*crscabsi2(11,indexint)
249
250c     $           o2colx(i)*crscabsi2(2,indexint)+
251c     $           h2colx(i)*crscabsi2(5,indexint)+
252c     $           nocolx(i)*crscabsi2(10,indexint)+
253c     $           hcolx(i)*crscabsi2(12,indexint)+
254c     $           no2colx(i)*crscabsi2(13,indexint)
255         end do
256
257c     Interpolations
258
259         do i=1,nz2
260            auxi = nz2-i+1
261            !O2 tabulated coefficient
262c            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
263            !O3p tabulated coefficient
264            auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
265            !CO2 tabulated coefficient
266            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
267            !H2 tabulated coefficient
268c            auxjh2(i) = jabsifotsintpar(auxi,5,indexint)
269            !N2 tabulated coefficient
270            auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
271            !N tabulated coefficient
272c            auxjn(i) = jabsifotsintpar(auxi,9,indexint)
273            !CO tabulated coefficient
274            auxjco(i) = jabsifotsintpar(auxi,11,indexint)
275            !H tabulated coefficient
276c            auxjh(i) = jabsifotsintpar(auxi,12,indexint)
277            !tabulated column
278            auxcoltab(i) = c1_16(auxi,indexint)
279         enddo
280c         !Only if chemthermod.ge.2
281c         if(chemthermod.ge.2) then
282c            do i=1,nz2
283c               auxi = nz2-i+1
284c               !N tabulated coefficient
285c               auxjn(i) = jabsifotsintpar(auxi,9,indexint)
286c               !NO tabulated coefficient
287c               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
288c               !NO2 tabulated coefficient
289c               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
290c            enddo
291c         endif
292
293          call interfast(wm,wp,auxind,auxcolinp,klev,
294     $        auxcoltab,nz2,limdown,limup)
295          do i=1,klev
296             ind=auxind(i)
297             auxi = klev-i+1
298             !O2 interpolated coefficient
299c             jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
300c     $            wp(i)*auxjo2(ind)
301             !O3p interpolated coefficient
302             jfotsout(indexint,3,auxi) = wm(i)*auxjo3p(ind+1) +
303     $            wp(i)*auxjo3p(ind)
304             !CO2 interpolated coefficient
305             jfotsout(indexint,1,auxi) = wm(i)*auxjco2(ind+1) +
306     $            wp(i)*auxjco2(ind)
307             !H2 interpolated coefficient
308c             jfotsout(indexint,5,auxi) = wm(i)*auxjh2(ind+1) +
309c     $            wp(i)*auxjh2(ind)
310             !N2 interpolated coefficient
311             jfotsout(indexint,8,auxi) = wm(i)*auxjn2(ind+1) +
312     $            wp(i)*auxjn2(ind)             
313             !N interpolated coefficient
314c             jfotsout(indexint,9,auxi) = wm(i)*auxjn(ind+1) +
315c     $            wp(i)*auxjn(ind)
316             !CO interpolated coefficient
317             jfotsout(indexint,11,auxi) = wm(i)*auxjco(ind+1) +
318     $            wp(i)*auxjco(ind)
319             !H interpolated coefficient
320c             jfotsout(indexint,12,auxi) = wm(i)*auxjh(ind+1) +
321c     $            wp(i)*auxjh(ind)             
322          enddo
323          !Only if chemthermod.ge.2
324c          if(chemthermod.ge.2) then
325c             do i=1,klev
326c                ind=auxind(i)
327c                auxi = klev-i+1
328c                !N interpolated coefficient
329c                jfotsout(indexint,9,auxi) = wm(i)*auxjn(ind+1) +
330c     $               wp(i)*auxjn(ind)
331c                !NO interpolated coefficient
332c                jfotsout(indexint,10,auxi)=wm(i)*auxjno(ind+1) +
333c     $               wp(i)*auxjno(ind)
334c                !NO2 interpolated coefficient
335c                jfotsout(indexint,13,auxi)=wm(i)*auxjno2(ind+1)+
336c     $               wp(i)*auxjno2(ind)
337c             enddo
338c          endif   
339      end do
340
341c     End intervals 2-15
342
343       
344ccccccccccccccccccccccccccccccc
345c     80.6-90.8nm (int16)
346c
347c     Absorption by:
348c     CO2, O2, O, N2, N, NO,
349c     CO, H, NO2
350ccccccccccccccccccccccccccccccc
351
352c     Input atmospheric column
353      indexint=16
354      do i=1,klev
355         auxcolinp(klev-i+1) = co2colx(i)*crscabsi2(1,indexint)+
356c     $        o2colx(i)*crscabsi2(2,indexint)+
357     $        o3pcolx(i)*crscabsi2(3,indexint)+
358     $        n2colx(i)*crscabsi2(8,indexint)+
359     $        cocolx(i)*crscabsi2(11,indexint)
360c     $        hcolx(i)*crscabsi2(12,indexint)+
361c     $        no2colx(i)*crscabsi2(13,indexint)
362      end do
363
364c     Interpolations
365
366      do i=1,nz2
367         auxi = nz2-i+1
368         !O2 tabulated coefficient
369c         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
370         !CO2 tabulated coefficient
371         auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
372         !O3p tabulated coefficient
373         auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
374         !N2 tabulated coefficient
375         auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
376         !CO tabulated coefficient
377         auxjco(i) = jabsifotsintpar(auxi,11,indexint)
378c        !N tabulated coefficient
379c         auxjn(i) = jabsifotsintpar(auxi,9,indexint)
380
381         !NO tabulated coefficient
382c         auxjno(i) = jabsifotsintpar(auxi,10,indexint)
383         !H tabulated coefficient
384c         auxjh(i) = jabsifotsintpar(auxi,12,indexint)
385         !NO2 tabulated coefficient
386c         auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
387
388         !Tabulated column
389         auxcoltab(i) = c1_16(auxi,indexint)
390      enddo
391      !Only if chemthermod.ge.2
392c      if(chemthermod.ge.2) then
393c         do i=1,nz2
394c            auxi = nz2-i+1
395c            !N tabulated coefficient
396c            auxjn(i) = jabsifotsintpar(auxi,9,indexint)
397c            !NO tabulated coefficient
398c            auxjno(i) = jabsifotsintpar(auxi,10,indexint)
399c            !NO2 tabulated coefficient
400c            auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
401c         enddo
402c      endif
403
404      call interfast
405     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
406      do i=1,klev
407         ind=auxind(i)
408         auxi = klev-i+1
409         !O2 interpolated coefficient
410c         jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
411c     $            wp(i)*auxjo2(ind)
412         !CO2 interpolated coefficient
413         jfotsout(indexint,1,auxi) = wm(i)*auxjco2(ind+1) +
414     $            wp(i)*auxjco2(ind)
415         !O3p interpolated coefficient
416         jfotsout(indexint,3,auxi) = wm(i)*auxjo3p(ind+1) +
417     $            wp(i)*auxjo3p(ind)
418         !N2 interpolated coefficient
419         jfotsout(indexint,8,auxi) = wm(i)*auxjn2(ind+1) +
420     $            wp(i)*auxjn2(ind)
421         !CO interpolated coefficient
422         jfotsout(indexint,11,auxi) = wm(i)*auxjco(ind+1) +
423     $            wp(i)*auxjco(ind) 
424         !N interpolated coefficient
425c          jfotsout(indexint,9,auxi) = wm(i)*auxjn(ind+1) +
426c     $           wp(i)*auxjn(ind)
427
428c         !H interpolated coefficient
429c         jfotsout(indexint,12,auxi) = wm(i)*auxjh(ind+1) +
430c     $            wp(i)*auxjh(ind)         
431      enddo
432      !Only if chemthermod.ge.2
433c      if(chemthermod.ge.2) then
434c         do i=1,klev
435c            ind=auxind(i)
436c            auxi = klev-i+1
437c            !N interpolated coefficient
438c            jfotsout(indexint,9,auxi) = wm(i)*auxjn(ind+1) +
439c     $           wp(i)*auxjn(ind)
440c            !NO interpolated coefficient
441c            jfotsout(indexint,10,auxi) = wm(i)*auxjno(ind+1) +
442c     $           wp(i)*auxjno(ind)
443c            !NO2 interpolated coefficient
444c            jfotsout(indexint,13,auxi) = wm(i)*auxjno2(ind+1) +
445c     $           wp(i)*auxjno2(ind)
446c         enddo
447c      endif
448c     End interval 16
449
450       
451ccccccccccccccccccccccccccccccc
452c     90.9-119.5nm (int 17-24)
453c
454c     Absorption by:
455c     CO2, O2, N2, NO, CO, NO2
456ccccccccccccccccccccccccccccccc
457
458c     Input column
459
460      do i=1,klev
461         auxcolinp(klev-i+1) = co2colx(i) +  n2colx(i) +
462     $         + cocolx(i)
463      end do
464
465      do indexint=17,24
466
467c     Interpolations
468
469         do i=1,nz2
470            auxi = nz2-i+1
471            !CO2 tabulated coefficient
472            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
473            !O2 tabulated coefficient
474c            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
475            !N2 tabulated coefficient
476            auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
477            !CO tabulated coefficient
478            auxjco(i) = jabsifotsintpar(auxi,11,indexint)           
479            !Tabulated column
480            auxcoltab(i) = c17_24(auxi)
481         enddo
482         !Only if chemthermod.ge.2
483c         if(chemthermod.ge.2) then
484c            do i=1,nz2
485c               auxi = nz2-i+1
486c               !NO tabulated coefficient
487c               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
488c               !NO2 tabulated coefficient
489c               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
490c            enddo
491c         endif
492
493         call interfast
494     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
495         !Correction to include T variation of CO2 cross section
496         if(indexint.eq.24) then
497            do i=1,klev
498               auxi = klev-i+1
499               if(sigma(indexint,auxi)*
500     $              alfa(indexint,auxi)*coltemp(auxi)
501     $              .lt.60.) then
502                  cortemp(i)=exp(-sigma(indexint,auxi)*
503     $                alfa(indexint,auxi)*coltemp(auxi))
504               else
505                  cortemp(i)=0.
506               end if
507            enddo
508         else
509            do i=1,klev
510               cortemp(i)=1.
511            enddo
512         end if
513         do i=1,klev           
514            ind=auxind(i)
515            auxi = klev-i+1
516            !O2 interpolated coefficient
517c            jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
518c     $           wp(i)*auxjo2(ind)) * cortemp(i)
519            !CO2 interpolated coefficient
520            jfotsout(indexint,1,auxi) = (wm(i)*auxjco2(ind+1) +
521     $           wp(i)*auxjco2(ind)) * cortemp(i)
522            if(indexint.eq.24) jfotsout(indexint,1,auxi)=
523     $           jfotsout(indexint,1,auxi)*
524     $           (1+alfa(indexint,auxi)*
525     $           (t2(auxi)-t0(auxi)))
526            !N2 interpolated coefficient
527            jfotsout(indexint,8,auxi) = (wm(i)*auxjn2(ind+1) +
528     $            wp(i)*auxjn2(ind)) * cortemp(i)           
529            !CO interpolated coefficient
530            jfotsout(indexint,11,auxi) = (wm(i)*auxjco(ind+1) +
531     $            wp(i)*auxjco(ind)) * cortemp(i)           
532         enddo
533         !Only if chemthermod.ge.2
534c         if(chemthermod.ge.2) then
535c            do i=1,klev
536c               ind=auxind(i)
537c               auxi = klev-i+1
538c               !NO interpolated coefficient
539c               jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
540c     $              wp(i)*auxjno(ind)) * cortemp(i)
541c               !NO2 interpolated coefficient
542c               jfotsout(indexint,13,auxi)=(wm(i)*auxjno2(ind+1)+
543c     $              wp(i)*auxjno2(ind)) * cortemp(i)
544c            enddo
545c         endif               
546      end do
547c     End intervals 17-24
548
549       
550ccccccccccccccccccccccccccccccc
551c     119.6-167.0nm (int 25-29)
552c
553c     Absorption by:
554c     CO2, O2, H2O, H2O2, NO,
555c     CO, NO2
556ccccccccccccccccccccccccccccccc
557
558c     Input atmospheric column
559
560      do i=1,klev
561c         auxcolinp(klev-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
562c     $        h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
563         auxcolinp(klev-i+1) = co2colx(i)  + cocolx(i)
564
565      end do
566
567      do indexint=25,29
568
569c     Interpolations
570
571         do i=1,nz2
572            auxi = nz2-i+1
573            !CO2 tabulated coefficient
574            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
575            !O2 tabulated coefficient
576c            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
577            !H2O tabulated coefficient
578c            auxjh2o(i) = jabsifotsintpar(auxi,4,indexint)
579            !H2O2 tabulated coefficient
580c            auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)           
581            !CO tabulated coefficient
582            auxjco(i) = jabsifotsintpar(auxi,11,indexint)           
583            !Tabulated column
584            auxcoltab(i) = c25_29(auxi)
585         enddo
586         !Only if chemthermod.ge.2
587c         if(chemthermod.ge.2) then
588c            do i=1,nz2
589c               auxi = nz2-i+1
590c               !NO tabulated coefficient
591c               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
592c               !NO2 tabulated coefficient
593c               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
594c            enddo
595c         endif
596         call interfast
597     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
598         do i=1,klev
599            ind=auxind(i)
600            auxi = klev-i+1
601            !Correction to include T variation of CO2 cross section
602            if(sigma(indexint,auxi)*alfa(indexint,auxi)*
603     $           coltemp(auxi).lt.60.) then
604               cortemp(i)=exp(-sigma(indexint,auxi)*
605     $              alfa(indexint,auxi)*coltemp(auxi))
606            else
607               cortemp(i)=0.
608            end if
609            !CO2 interpolated coefficient
610            jfotsout(indexint,1,auxi) = (wm(i)*auxjco2(ind+1) +
611     $           wp(i)*auxjco2(ind)) * cortemp(i) *
612     $           (1+alfa(indexint,auxi)*
613     $           (t2(auxi)-t0(auxi)))
614            !O2 interpolated coefficient
615c            jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
616c     $           wp(i)*auxjo2(ind)) * cortemp(i)
617            !H2O interpolated coefficient
618c            jfotsout(indexint,4,auxi) = (wm(i)*auxjh2o(ind+1) +
619c     $           wp(i)*auxjh2o(ind)) * cortemp(i)
620            !H2O2 interpolated coefficient
621c            jfotsout(indexint,6,auxi) = (wm(i)*auxjh2o2(ind+1) +
622c     $           wp(i)*auxjh2o2(ind)) * cortemp(i)           
623            !CO interpolated coefficient
624            jfotsout(indexint,11,auxi) = (wm(i)*auxjco(ind+1) +
625     $           wp(i)*auxjco(ind)) * cortemp(i)
626         enddo
627         !Only if chemthermod.ge.2
628c         if(chemthermod.ge.2) then
629c            do i=1,klev
630c               ind=auxind(i)
631c               auxi = klev-i+1
632c               !NO interpolated coefficient
633c               jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
634c     $              wp(i)*auxjno(ind)) * cortemp(i)
635c               !NO2 interpolated coefficient
636c               jfotsout(indexint,13,auxi)=(wm(i)*auxjno2(ind+1)+
637c     $              wp(i)*auxjno2(ind)) * cortemp(i)
638c            enddo
639c         endif
640
641      end do
642
643c     End intervals 25-29
644
645
646cccccccccccccccccccccccccccccccc
647c     167.1-202.5nm (int 30-31)
648c   
649c     Absorption by:
650c     CO2, O2, H2O, H2O2, NO,
651c     NO2
652cccccccccccccccccccccccccccccccc
653
654c     Input atmospheric column
655
656      do i=1,klev
657         auxcolinp(klev-i+1) = co2colx(i)
658      end do
659
660c     Interpolation
661
662      do indexint=30,31
663
664         do i=1,nz2
665            auxi = nz2-i+1
666            !CO2 tabulated coefficient
667            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
668            !O2 tabulated coefficient
669c            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
670            !H2O tabulated coefficient
671c            auxjh2o(i) = jabsifotsintpar(auxi,4,indexint)
672            !H2O2 tabulated coefficient
673c            auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)           
674            !Tabulated column
675c            auxcoltab(i) = c30_31(auxi)
676         enddo
677         !Only if chemthermod.ge.2
678c         if(chemthermod.ge.2) then
679c            do i=1,nz2
680c               auxi = nz2-i+1
681c               !NO tabulated coefficient
682c               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
683c               !NO2 tabulated coefficient
684c               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
685c            enddo
686c         endif
687
688         call interfast
689     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
690         do i=1,klev
691            ind=auxind(i)
692            auxi = klev-i+1
693            !Correction to include T variation of CO2 cross section
694            if(sigma(indexint,auxi)*alfa(indexint,auxi)*
695     $           coltemp(auxi).lt.60.) then
696               cortemp(i)=exp(-sigma(indexint,auxi)*
697     $              alfa(indexint,auxi)*coltemp(auxi))
698            else
699               cortemp(i)=0.
700            end if
701            !CO2 interpolated coefficient
702            jfotsout(indexint,1,auxi) = (wm(i)*auxjco2(ind+1) +
703     $           wp(i)*auxjco2(ind)) * cortemp(i) *
704     $           (1+alfa(indexint,auxi)*
705     $           (t2(auxi)-t0(auxi)))
706            !O2 interpolated coefficient
707c            jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
708c     $            wp(i)*auxjo2(ind)) * cortemp(i)
709            !H2O interpolated coefficient
710c            jfotsout(indexint,4,auxi) = (wm(i)*auxjh2o(ind+1) +
711c     $            wp(i)*auxjh2o(ind)) * cortemp(i)
712            !H2O2 interpolated coefficient
713c            jfotsout(indexint,6,auxi) = (wm(i)*auxjh2o2(ind+1) +
714c     $            wp(i)*auxjh2o2(ind)) * cortemp(i)           
715         enddo
716c         !Only if chemthermod.ge.2
717c         if(chemthermod.ge.2) then
718c            do i=1,klev
719c               ind=auxind(i)
720c               auxi = klev-i+1
721c               !NO interpolated coefficient
722c               jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
723c     $              wp(i)*auxjno(ind)) * cortemp(i)
724c               !NO2 interpolated coefficient
725c               jfotsout(indexint,13,auxi)=(wm(i)*auxjno2(ind+1)+
726c     $              wp(i)*auxjno2(ind)) * cortemp(i)
727c            enddo
728c         endif
729
730      end do
731
732c     End intervals 30-31
733
734
735ccccccccccccccccccccccccccccccc
736c     202.6-210.0nm (int 32)
737c
738c     Absorption by:
739c     CO2, O2, H2O2, NO, NO2
740ccccccccccccccccccccccccccccccc
741
742c     Input atmospheric column
743
744      indexint=32
745      do i=1,klev
746         auxcolinp(klev-i+1) =co2colx(i)
747      end do
748
749c     Interpolation
750
751      do i=1,nz2
752         auxi = nz2-i+1
753         !CO2 tabulated coefficient
754         auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
755         !O2 tabulated coefficient
756c         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
757         !H2O2 tabulated coefficient
758c         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)         
759         !Tabulated column
760         auxcoltab(i) = c32(auxi)
761      enddo
762      !Only if chemthermod.ge.2
763c      if(chemthermod.ge.2) then
764c         do i=1,nz2
765c            auxi = nz2-i+1
766c            !NO tabulated coefficient
767c            auxjno(i) = jabsifotsintpar(auxi,10,indexint)
768c            !NO2 tabulated coefficient
769c            auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
770c         enddo
771c      endif
772      call interfast
773     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
774      do i=1,klev
775         ind=auxind(i)
776         auxi = klev-i+1
777         !Correction to include T variation of CO2 cross section
778         if(sigma(indexint,klev-i+1)*alfa(indexint,auxi)*
779     $        coltemp(auxi).lt.60.) then
780            cortemp(i)=exp(-sigma(indexint,auxi)*
781     $           alfa(indexint,auxi)*coltemp(auxi))
782         else
783            cortemp(i)=0.
784         end if
785         !CO2 interpolated coefficient
786         jfotsout(indexint,1,auxi) = (wm(i)*auxjco2(ind+1) +
787     $        wp(i)*auxjco2(ind)) * cortemp(i) *
788     $        (1+alfa(indexint,auxi)*
789     $        (t2(auxi)-t0(auxi)))
790         !O2 interpolated coefficient
791c         jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
792c     $        wp(i)*auxjo2(ind)) * cortemp(i)
793         !H2O2 interpolated coefficient
794c         jfotsout(indexint,6,auxi) = (wm(i)*auxjh2o2(ind+1) +
795c     $        wp(i)*auxjh2o2(ind)) * cortemp(i)         
796      enddo
797      !Only if chemthermod.ge.2
798c      if(chemthermod.ge.2) then
799c         do i=1,klev
800c            auxi = klev-i+1
801c            ind=auxind(i)
802c            !NO interpolated coefficient
803c            jfotsout(indexint,10,auxi) = (wm(i)*auxjno(ind+1) +
804c     $           wp(i)*auxjno(ind)) * cortemp(i)
805           !NO2 interpolated coefficient
806c            jfotsout(indexint,13,auxi) = (wm(i)*auxjno2(ind+1) +
807c     $           wp(i)*auxjno2(ind)) * cortemp(i)
808c         enddo
809c      endif
810
811c     End of interval 32
812
813
814ccccccccccccccccccccccccccccccc
815c     210.1-231.0nm (int 33)
816c     
817c     Absorption by:
818c     O2, H2O2, NO2
819ccccccccccccccccccccccccccccccc
820
821c     Input atmospheric column
822
823c      indexint=33
824c      do i=1,klev
825c         auxcolinp(klev-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
826c      end do
827
828c     Interpolation
829
830c      do i=1,nz2
831c         auxi = nz2-i+1
832         !O2 tabulated coefficient
833c         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
834         !H2O2 tabulated coefficient
835c         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
836         !Tabulated column
837c        auxcoltab(i) = c33(auxi)
838c     enddo
839      !Only if chemthermod.ge.2
840c     if(chemthermod.ge.2) then
841c        do i=1,nz2
842            !NO2 tabulated coefficient
843c           auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
844c        enddo
845c     endif
846c     call interfast
847c    $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
848c     do i=1,klev
849c        ind=auxind(i)
850c        auxi = klev-i+1
851         !O2 interpolated coefficient
852c        jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
853c    $        wp(i)*auxjo2(ind)
854c        !H2O2 interpolated coefficient
855c        jfotsout(indexint,6,auxi) = wm(i)*auxjh2o2(ind+1) +
856c    $        wp(i)*auxjh2o2(ind)         
857c     enddo
858      !Only if chemthermod.ge.2
859c     if(chemthermod.ge.2) then
860c        do i=1,klev
861c           ind=auxind(i)
862c           !NO2 interpolated coefficient
863c           jfotsout(indexint,13,klev-i+1) = wm(i)*auxjno2(ind+1) +
864c    $           wp(i)*auxjno2(ind)
865c        enddo
866c     endif
867
868c     End of interval 33
869
870
871ccccccccccccccccccccccccccccccc
872c     231.1-240.0nm (int 34)
873c
874c     Absorption by:
875c     O2, H2O2, O3, NO2
876ccccccccccccccccccccccccccccccc
877
878c     Input atmospheric column
879
880c     indexint=34
881c     do i=1,klev
882c        auxcolinp(klev-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) +
883c    $        no2colx(i)
884c     end do
885
886c     Interpolation
887
888c     do i=1,nz2
889c        auxi = nz2-i+1
890         !O2 tabulated coefficient
891c        auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
892         !H2O2 tabulated coefficient
893c        auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
894         !O3 tabulated coefficient
895c        auxjo3(i) = jabsifotsintpar(auxi,7,indexint)         
896         !Tabulated column
897c        auxcoltab(i) = c34(nz2-i+1)
898c     enddo
899      !Only if chemthermod.ge.2
900c     if(chemthermod.ge.2) then
901c        do i=1,nz2
902            !NO2 tabulated coefficient
903c           auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
904c        enddo
905c     endif
906c     call interfast
907c    $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
908c     do i=1,klev
909c        ind=auxind(i)
910c        auxi = klev-i+1
911         !O2 interpolated coefficient
912c        jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
913c    $        wp(i)*auxjo2(ind)
914         !H2O2 interpolated coefficient
915c        jfotsout(indexint,6,auxi) = wm(i)*auxjh2o2(ind+1) +
916c    $        wp(i)*auxjh2o2(ind)
917         !O3 interpolated coefficient
918c        jfotsout(indexint,7,auxi) = wm(i)*auxjo3(ind+1) +
919c    $        wp(i)*auxjo3(ind)         
920c     enddo
921      !Only if chemthermod.ge.2
922c     if(chemthermod.ge.2) then
923c        do i=1,klev
924c           ind=auxind(i)
925            !NO2 interpolated coefficient
926c           jfotsout(indexint,13,klev-i+1) = wm(i)*auxjno2(ind+1) +
927c    $           wp(i)*auxjno2(ind)
928c        enddo
929c     endif
930
931c     End of interval 34     
932
933
934ccccccccccccccccccccccccccccccc
935c     240.1-337.7nm (int 35)
936c
937c     Absorption by:
938c     H2O2, O3, NO2
939ccccccccccccccccccccccccccccccc
940
941c     Input atmospheric column
942
943      indexint=35
944c      do i=1,klev
945c         auxcolinp(klev-i+1) = o3colx(i)
946c      end do
947c
948c     Interpolation
949
950c      do i=1,nz2
951c         auxi = nz2-i+1
952         !H2O2 tabulated coefficient
953c         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
954         !O3 tabulated coefficient
955c         auxjo3(i) = jabsifotsintpar(auxi,7,indexint)
956         !Tabulated column
957c         auxcoltab(i) = c35(auxi)
958c      enddo
959      !Only if chemthermod.ge.2
960c      if(chemthermod.ge.2) then
961c         do i=1,nz2
962c            !NO2 tabulated coefficient
963c            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
964c         enddo
965c      endif
966c      call interfast
967c     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
968c      do i=1,klev
969c         ind=auxind(i)
970c         auxi = klev-i+1
971c         !H2O2 interpolated coefficient
972c         jfotsout(indexint,6,auxi) = wm(i)*auxjh2o2(ind+1) +
973c     $        wp(i)*auxjh2o2(ind)
974c         !O3 interpolated coefficient
975c         jfotsout(indexint,7,auxi) = wm(i)*auxjo3(ind+1) +
976c     $        wp(i)*auxjo3(ind)         
977c      enddo
978c      if(chemthermod.ge.2) then
979c         do i=1,klev
980c            ind=auxind(i)
981c            !NO2 interpolated coefficient
982c            jfotsout(indexint,13,klev-i+1) = wm(i)*auxjno2(ind+1) +
983c     $           wp(i)*auxjno2(ind)
984c         enddo
985c      endif
986
987c     End of interval 35
988
989ccccccccccccccccccccccccccccccc
990c     337.8-800.0 nm (int 36)
991c     
992c     Absorption by:
993c     O3, NO2
994ccccccccccccccccccccccccccccccc
995
996c     Input atmospheric column
997
998      indexint=36
999c      do i=1,klev
1000c         auxcolinp(klev-i+1) = o3colx(i)
1001c      end do
1002
1003c     Interpolation
1004
1005c      do i=1,nz2
1006c         auxi = nz2-i+1
1007         !O3 tabulated coefficient
1008c         auxjo3(i) = jabsifotsintpar(auxi,7,indexint)         
1009         !Tabulated column
1010c         auxcoltab(i) = c36(auxi)
1011c      enddo
1012      !Only if chemthermod.ge.2
1013c      if(chemthermod.ge.2) then
1014c         do i=1,nz2
1015c            !NO2 tabulated coefficient
1016c            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
1017c         enddo
1018c      endif
1019c      call interfast
1020c     $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
1021c      do i=1,klev
1022c         ind=auxind(i)
1023c         !O3 interpolated coefficient
1024c         jfotsout(indexint,7,klev-i+1) = wm(i)*auxjo3(ind+1) +
1025c     $        wp(i)*auxjo3(ind)         
1026c      enddo
1027c      !Only if chemthermod.ge.2
1028c      if(chemthermod.ge.2) then
1029c         do i=1,klev
1030c            ind=auxind(i)
1031c           !NO2 interpolated coefficient
1032c            jfotsout(indexint,13,klev-i+1) = wm(i)*auxjno2(ind+1) +
1033c     $           wp(i)*auxjno2(ind)
1034c         enddo
1035c      endif
1036
1037c     End of interval 36
1038
1039c     End of interpolation to obtain photoabsorption rates
1040
1041       
1042      return
1043      end
1044
1045
1046
1047c**********************************************************************
1048c**********************************************************************
1049
1050      subroutine column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
1051     $     co2colx,o3pcolx, n2colx, cocolx)
1052
1053c     mar 2014        gg            adapted to Venus GCM
1054c     nov 2002        fgg           first version
1055
1056c**********************************************************************
1057      use dimphy
1058      use conc
1059      implicit none
1060
1061
1062c     common variables and constants
1063c#include "tracer.h"
1064#include "param.h"
1065#include "param_v4.h"
1066#include "clesphys.h"
1067#include "mmol.h"
1068
1069
1070c    local parameters and variables
1071
1072c     input and output variables
1073
1074      integer    ig
1075      integer    chemthermod
1076      integer    nesptherm                   !# of species undergoing chemistry, input
1077      real       rm(klev,nesptherm)         !densities (cm-3), input
1078      real       tx(klev)                   !temperature profile, input
1079      real       iz(klev+1)                 !height profile, input
1080      real       zenit                      !SZA, input
1081      real       co2colx(klev)              !column density of CO2 (cm^-2), output
1082      real       o3pcolx(klev)              !column density of O(3P)(cm^-2), output
1083      real       n2colx(klev)               !N2 column density (cm-2), output
1084      real       cocolx(klev)               !CO column density (cm-2), output
1085
1086c      real       o2colx(klev)               !column density of O2(cm^-2), output
1087c      real       h2colx(klev)               !H2 column density (cm-2), output
1088c      real       h2ocolx(klev)              !H2O column density (cm-2), output
1089c      real       h2o2colx(klev)             !column density of H2O2(cm^-2), output
1090c      real       o3colx(klev)               !O3 column density (cm-2), output
1091c      real       nocolx(klev)               !NO column density (cm-2), output
1092c      real       hcolx(klev)                !H column density (cm-2), output
1093c      real       no2colx(klev)              !NO2 column density (cm-2), output
1094     
1095
1096c     local variables
1097
1098      real       xx
1099      real       grav(klev)
1100      real       Hco2,Ho3p,Ho2,Hh2,Hh2o,Hh2o2
1101      real       Ho3,Hn2,Hn,Hno,Hco,Hh,Hno2
1102
1103      real       co2x(klev)
1104      real       o3px(klev)
1105      real       cox(klev)
1106      real       n2x(klev)
1107      real       nx(klev)
1108
1109c      real       o2x(klev)
1110c      real       o3x(klev)
1111c      real       hx(klev)
1112c      real       h2x(klev)
1113c      real       h2ox(klev)
1114c      real       h2o2x(klev)
1115c      real       nox(klev)
1116c      real       no2x(klev)
1117
1118      integer    i,j,k,icol,indexint          !indexes
1119
1120c     variables for optical path calculation
1121
1122      integer    nz3
1123!      parameter  (nz3=nz*2)
1124
1125      integer    jj
1126      real*8      esp(klev*2)
1127      real*8      ilayesp(klev*2)
1128      real*8      szalayesp(klev*2)
1129      integer     nlayesp
1130      real*8      zmini
1131      real*8      depth
1132      real*8      espco2, espo2, espo3p, esph2, esph2o, esph2o2,espo3
1133      real*8      espn2,espn,espno,espco,esph,espno2
1134      real*8      rcmnz, rcmmini
1135      real*8      szadeg
1136
1137! Tracer indexes in the thermospheric chemistry:
1138!!! ATTENTION. These values have to be identical to those in euvheat.F90
1139!!! If the values are changed there, the same has to be done here  !!!
1140
1141      integer,parameter :: ix_co2=1
1142      integer,parameter :: ix_n2=13
1143      integer,parameter :: ix_o=3
1144      integer,parameter :: ix_co=4
1145
1146c*************************PROGRAM STARTS*******************************
1147
1148      nz3 = klev*2
1149      do i=1,klev
1150        xx = ( radio + iz(i) ) * 1.e5    ! conversion [km] ---> [cm]
1151         grav(i) = gg * masa /(xx**2)    ! [cm/s2]
1152      end do
1153
1154      !Scale heights  H = kT /Mg  -->  [cm]
1155      xx = kboltzman * tx(klev) * n_avog / grav(klev)   ! g cm mol-1
1156
1157      Ho3p  = xx / mmolo
1158      Hco2  = xx / mmolco2
1159      Hco   = xx / mmolco
1160      Hn2   = xx / mmoln2
1161      Hn    = xx / mmoln
1162
1163      !Only if O3 chem. required
1164c      if(chemthermod.ge.1)
1165!     $     Ho3   = xx / mmol(igcm_o3)
1166c     $     Ho3   = xx / mmolo3
1167c      !Only if N or ion chem.
1168c      if(chemthermod.ge.2) then
1169c         Hn2   = xx / mmoln2
1170c         Hn    = xx / mmoln
1171c         Hno   = xx / mmolno
1172c         Hno2  = xx / mmolno2
1173c      endif
1174      ! first loop in altitude : initialisation
1175      do i=klev,1,-1
1176         !Column initialisation
1177         co2colx(i)  = 0.
1178         o3pcolx(i)  = 0.
1179         n2colx(i)   = 0.
1180         cocolx(i)   = 0.
1181
1182         !--Densities [cm-3]
1183         co2x(i)  = rm(i,ix_co2)
1184         o3px(i)  = rm(i,ix_o)
1185         cox(i)   = rm(i,ix_co)
1186         n2x(i)   = rm(i,ix_n2)
1187
1188c         write(*,*), '--jthermcalc--', co2x(i)
1189
1190         !Only if O3 chem. required
1191c         if(chemthermod.ge.1)
1192c     $        o3x(i)   = rm(i,i_o3)
1193c         !Only if Nitrogen of ion chemistry requested
1194c         if(chemthermod.ge.2) then
1195c            n2x(i)   = rm(i,i_n2)
1196c            nx(i)    = rm(i,i_n)
1197c            nox(i)   = rm(i,i_no)
1198c            no2x(i)  = rm(i,i_no2)
1199c         endif
1200      enddo     ! end first loop
1201      ! second loop in altitude : column calculations
1202      do i=klev,1,-1
1203         !Routine to calculate the geometrical length of each layer
1204         call espesor_optico_A(ig,i,zenit,iz(i),nz3,iz,esp,ilayesp,
1205     $         szalayesp,nlayesp, zmini)
1206         if(ilayesp(nlayesp).eq.-1) then
1207            co2colx(i)=1.e25
1208            o3pcolx(i)=1.e25
1209            n2colx(i)=1.e25
1210            cocolx(i)=1.e25
1211
1212c            o2colx(i)=1.e25
1213c            o3pcolx(i)=1.e25
1214c            h2colx(i)=1.e25
1215c            h2ocolx(i)=1.e25
1216c            h2o2colx(i)=1.e25
1217c            o3colx(i)=1.e25
1218c            ncolx(i)=1.e25
1219c            nocolx(i)=1.e25
1220c            cocolx(i)=1.e25
1221c            hcolx(i)=1.e25
1222c            no2colx(i)=1.e25
1223         else
1224            rcmnz = ( radio + iz(klev) ) * 1.e5    ! km --> cm
1225            rcmmini = ( radio + zmini ) * 1.e5
1226            !Column calculation taking into account the geometrical depth
1227            !calculated before
1228            do j=1,nlayesp
1229               jj=ilayesp(j)
1230               !Top layer
1231               if(jj.eq.klev) then
1232                  if(zenit.le.60.) then
1233                     o3pcolx(i)=o3pcolx(i)+o3px(klev)*Ho3p*esp(j)
1234     $                    *1.e-5
1235                     co2colx(i)=co2colx(i)+co2x(klev)*Hco2*esp(j)
1236     $                    *1.e-5
1237                     cocolx(i)=cocolx(i)+cox(klev)*Hco*esp(j)
1238     $                    *1.e-5
1239                     n2colx(i)=n2colx(i)+n2x(klev)*Hn2*esp(j)
1240     $                    *1.e-5
1241
1242c                     h2o2colx(i)=h2o2colx(i)+
1243c     $                    h2o2x(klev)*Hh2o2*esp(j)*1.e-5
1244c                     o2colx(i)=o2colx(i)+o2x(klev)*Ho2*esp(j)
1245c     $                    *1.e-5
1246c                     h2colx(i)=h2colx(i)+h2x(klev)*Hh2*esp(j)
1247c     $                    *1.e-5
1248c                     h2ocolx(i)=h2ocolx(i)+h2ox(klev)*Hh2o*esp(j)
1249c     $                    *1.e-5                     
1250c                     cocolx(i)=cocolx(i)+cox(klev)*Hco*esp(j)
1251c     $                    *1.e-5
1252c                     hcolx(i)=hcolx(i)+hx(klev)*Hh*esp(j)
1253c     $                    *1.e-5
1254                     !Only if O3 chemistry required
1255c                     if(chemthermod.ge.1) o3colx(i)=
1256c     $                    o3colx(i)+o3x(klev)*Ho3*esp(j)
1257c     $                    *1.e-5
1258                     !Only if N or ion chemistry requested
1259c                     if(chemthermod.ge.2) then
1260c                        n2colx(i)=n2colx(i)+n2x(klev)*Hn2*esp(j)
1261c     $                    *1.e-5
1262
1263c                     endif
1264                  else if(zenit.gt.60.) then
1265                     espco2 =sqrt((rcmnz+Hco2)**2 -rcmmini**2) - esp(j)
1266                     espo3p = sqrt((rcmnz+Ho3p)**2 -rcmmini**2)- esp(j)
1267                     espco  = sqrt((rcmnz+Hco)**2 -rcmmini**2) - esp(j)
1268                     espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j)
1269                     espn  =sqrt((rcmnz+Hn)**2-rcmmini**2)  - esp(j)
1270
1271c                     espo2  = sqrt((rcmnz+Ho2)**2 -rcmmini**2) - esp(j)
1272c                     esph2  = sqrt((rcmnz+Hh2)**2 -rcmmini**2) - esp(j)
1273c                     esph2o = sqrt((rcmnz+Hh2o)**2 -rcmmini**2)- esp(j)
1274c                     esph2o2= sqrt((rcmnz+Hh2o2)**2-rcmmini**2)- esp(j)
1275c                     esph   = sqrt((rcmnz+Hh)**2 -rcmmini**2)  - esp(j)
1276                     !Only if O3 chemistry required
1277c                     if(chemthermod.ge.1)                     
1278c     $                   espo3=sqrt((rcmnz+Ho3)**2-rcmmini**2)-esp(j)
1279c                     !Only if N or ion chemistry requested
1280c                     if(chemthermod.ge.2) then
1281c                        espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j)
1282c                        espn  =sqrt((rcmnz+Hn)**2-rcmmini**2)  - esp(j)
1283c                        espno =sqrt((rcmnz+Hno)**2-rcmmini**2) - esp(j)
1284c                        espno2=sqrt((rcmnz+Hno2)**2-rcmmini**2)- esp(j)
1285c                     endif
1286
1287                     co2colx(i) = co2colx(i) + espco2*co2x(klev)
1288                     o3pcolx(i) = o3pcolx(i) + espo3p*o3px(klev)
1289                     cocolx(i)  = cocolx(i)  + espco*cox(klev)
1290                     n2colx(i)  = n2colx(i)  + espn2*n2x(klev)
1291
1292c                     o2colx(i)  = o2colx(i)  + espo2*o2x(klev)
1293c                     h2colx(i)  = h2colx(i)  + esph2*h2x(klev)
1294c                     h2ocolx(i) = h2ocolx(i) + esph2o*h2ox(klev)
1295c                     h2o2colx(i)= h2o2colx(i)+ esph2o2*h2o2x(klev)
1296c                     cocolx(i)  = cocolx(i)  + espco*cox(klev)
1297c                     hcolx(i)   = hcolx(i)   + esph*hx(klev)
1298                     !Only if O3 chemistry required
1299c                     if(chemthermod.ge.1)                     
1300c     $                  o3colx(i) = o3colx(i)  + espo3*o3x(klev)
1301c                     !Only if N or ion chemistry requested
1302c                     if(chemthermod.ge.2) then
1303c                        n2colx(i)  = n2colx(i)  + espn2*n2x(klev)
1304c                        ncolx(i)   = ncolx(i)   + espn*nx(klev)
1305c                        nocolx(i)  = nocolx(i)  + espno*nox(klev)
1306c                        no2colx(i) = no2colx(i) + espno2*no2x(klev)
1307c                     endif
1308                  endif !Of if zenit.lt.60
1309               !Other layers
1310               else
1311                  co2colx(i)  = co2colx(i) +
1312     $                 esp(j) * (co2x(jj)+co2x(jj+1)) / 2.
1313                  o3pcolx(i)  = o3pcolx(i) +
1314     $                 esp(j) * (o3px(jj)+o3px(jj+1)) / 2.
1315                  cocolx(i)   = cocolx(i) +
1316     $                 esp(j) * (cox(jj)+cox(jj+1)) / 2.
1317                  n2colx(i)   = n2colx(i) +
1318     $                 esp(j) * (n2x(jj)+n2x(jj+1)) / 2.
1319
1320c
1321c                  o2colx(i)   = o2colx(i) +
1322c     $                 esp(j) * (o2x(jj)+o2x(jj+1)) / 2.
1323c                  h2colx(i)   = h2colx(i) +
1324c     $                 esp(j) * (h2x(jj)+h2x(jj+1)) / 2.
1325c                  h2ocolx(i)  = h2ocolx(i) +
1326c     $                 esp(j) * (h2ox(jj)+h2ox(jj+1)) / 2.
1327c                  h2o2colx(i) = h2o2colx(i) +
1328c     $                 esp(j) * (h2o2x(jj)+h2o2x(jj+1)) / 2.
1329c                  hcolx(i)    = hcolx(i) +
1330c     $                 esp(j) * (hx(jj)+hx(jj+1)) / 2.
1331                  !Only if O3 chemistry required
1332c                  if(chemthermod.ge.1)
1333c     $                 o3colx(i) = o3colx(i) +
1334c     $                 esp(j) * (o3x(jj)+o3x(jj+1)) / 2.
1335c                  !Only if N or ion chemistry requested
1336c                  if(chemthermod.ge.2) then
1337c                     n2colx(i)   = n2colx(i) +
1338c     $                 esp(j) * (n2x(jj)+n2x(jj+1)) / 2.
1339c                     ncolx(i)    = ncolx(i) +
1340c     $                    esp(j) * (nx(jj)+nx(jj+1)) / 2.
1341c                     nocolx(i)   = nocolx(i) +
1342c     $                    esp(j) * (nox(jj)+nox(jj+1)) / 2.
1343c                     no2colx(i)  = no2colx(i) +
1344c     $                    esp(j) * (no2x(jj)+no2x(jj+1)) / 2.
1345c                  endif
1346
1347               endif  !Of if jj.eq.klev
1348            end do    !Of do j=1,nlayesp
1349         endif        !Of ilayesp(nlayesp).eq.-1
1350      enddo           !Of do i=klev,1,-1
1351      return
1352
1353
1354      end
1355
1356
1357c**********************************************************************
1358c**********************************************************************
1359
1360      subroutine interfast(wm,wp,nm,p,nlayer,pin,nl,limdown,limup)
1361C
1362C subroutine to perform linear interpolation in pressure from 1D profile
1363C escin(nl) sampled on pressure grid pin(nl) to profile
1364C escout(nlayer) on pressure grid p(nlayer).
1365C
1366      real*8 wm(nlayer),wp(nlayer),p(nlayer)
1367      integer nm(nlayer)
1368      real*8 pin(nl)
1369      real*8 limup,limdown
1370      integer nl,nlayer,n1,n,np,nini
1371      nini=1
1372      do n1=1,nlayer
1373         if(p(n1) .gt. limup .or. p(n1) .lt. limdown) then
1374            wm(n1) = 0.d0
1375            wp(n1) = 0.d0
1376         else
1377            do n = nini,nl-1
1378               if (p(n1).ge.pin(n).and.p(n1).le.pin(n+1)) then
1379                  nm(n1)=n
1380                  np=n+1
1381                  wm(n1)=abs(pin(n)-p(n1))/(pin(np)-pin(n))
1382                  wp(n1)=1.d0 - wm(n1)
1383                  nini = n
1384                  exit
1385               endif
1386            enddo
1387         endif
1388      enddo
1389      return
1390      end
1391
1392
1393c**********************************************************************
1394c**********************************************************************
1395
1396      subroutine espesor_optico_A (ig,capa, szadeg,z,
1397     @                   nz3,iz,esp,ilayesp,szalayesp,nlayesp, zmini)
1398
1399c     fgg              nov 03      Adaptation to Martian model
1400c     malv             jul 03      Corrected z grid. Split in alt & frec codes
1401c     fgg              feb 03      first version
1402*************************************************************************
1403      use dimphy
1404      implicit none
1405
1406
1407c     common variables and constants
1408#include "param.h"
1409#include "param_v4.h"
1410
1411c     arguments
1412
1413      real        szadeg                ! I. SZA [rad]
1414      real        z                     ! I. altitude of interest [km]
1415      integer     nz3,ig                   ! I. dimension of esp, ylayesp, etc...
1416                                        !  (=2*klev= max# of layers in ray path)
1417      real     iz(klev+1)              ! I. Altitude of each layer
1418      real*8        esp(nz3)            ! O. layer widths after geometrically
1419                                        !    amplified; in [cm] except at TOA
1420                                        !    where an auxiliary value is used
1421      real*8        ilayesp(nz3)        ! O. Indexes of layers along ray path
1422      real*8        szalayesp(nz3)      ! O. SZA [deg]    "     "       "
1423      integer       nlayesp
1424!      real*8        nlayesp             ! O. # layers along ray path at this z
1425      real*8        zmini               ! O. Minimum altitud of ray path [km]
1426
1427
1428c     local variables and constants
1429
1430        integer     j,i,capa
1431        integer     jmin                  ! index of min.altitude along ray path
1432        real*8      szarad                ! SZA [deg]
1433        real*8      zz
1434        real*8      diz(klev+1)
1435        real*8      rkmnz                 ! distance TOA to center of Planet [km]
1436        real*8      rkmmini               ! distance zmini to center of P [km]
1437        real*8      rkmj                  ! intermediate distance to C of P [km]
1438
1439c external function
1440        external  grid_R8          ! Returns index of layer containing the altitude
1441                                ! of interest, z; for example, if
1442                                ! zkm(i)=z or zkm(i)<z<zkm(i+1) => grid(z)=i
1443        integer   grid_R8
1444
1445*************************************************************************     
1446        szarad = dble(szadeg)*3.141592d0/180.d0
1447        zz=dble(z)
1448        do i=1,klev
1449           diz(i)=dble(iz(i))
1450        enddo
1451        do j=1,nz3
1452          esp(j) = 0.d0
1453          szalayesp(j) = 777.d0
1454          ilayesp(j) = 0
1455        enddo
1456        nlayesp = 0
1457
1458        ! First case: szadeg<60
1459        ! The optical thickness will be given by  1/cos(sza)
1460        ! We deal with 2 different regions:
1461        !   1: First, all layers between z and ztop ("upper part of ray")
1462        !   2: Second, the layer at ztop
1463        if(szadeg.lt.60.d0) then
1464
1465           zmini = zz
1466           if(abs(zz-diz(klev)).lt.1.d-3) goto 1357
1467           ! 1st Zone: Upper part of ray
1468           !
1469           do j=grid_R8(zz,diz,klev),klev-1
1470             nlayesp = nlayesp + 1
1471             ilayesp(nlayesp) = j
1472             esp(nlayesp) = (diz(j+1)-diz(j)) / cos(szarad)        ! [km]
1473             esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
1474             szalayesp(nlayesp) = szadeg
1475           end do
1476
1477           !
1478           ! 2nd Zone: Top layer
1479 1357      continue
1480           nlayesp = nlayesp + 1
1481           ilayesp(nlayesp) = klev
1482           esp(nlayesp) = 1.d0 / cos(szarad)         ! aux. non-dimens. factor
1483           szalayesp(nlayesp) = szadeg
1484
1485
1486        ! Second case:  60 < szadeg < 90
1487        ! The optical thickness is evaluated.
1488        !    (the magnitude of the effect of not using cos(sza) is 3.e-5
1489        !     for z=60km & sza=30, and 5e-4 for z=60km & sza=60, approximately)
1490        ! We deal with 2 different regions:
1491        !   1: First, all layers between z and ztop ("upper part of ray")
1492        !   2: Second, the layer at ztop ("uppermost layer")
1493        else if(szadeg.le.90.d0.and.szadeg.ge.60.d0) then
1494
1495           zmini=(radio+zz)*sin(szarad)-radio
1496           rkmmini = radio + zmini
1497
1498           if(abs(zz-diz(klev)).lt.1.d-4) goto 1470
1499
1500           ! 1st Zone: Upper part of ray
1501           !
1502           do j=grid_R8(zz,diz,klev),klev-1
1503              nlayesp = nlayesp + 1
1504              ilayesp(nlayesp) = j
1505              esp(nlayesp) =
1506     #             sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) -
1507     #             sqrt( (radio+diz(j))**2 - rkmmini**2 )           ! [km]
1508              esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
1509              rkmj = radio+diz(j)
1510              szalayesp(nlayesp) = asin( rkmmini/rkmj )             ! [rad]
1511              szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592 ! [deg]
1512           end do
1513 1470      continue
1514           ! 2nd Zone:  Uppermost layer of ray.
1515           !
1516           nlayesp = nlayesp + 1
1517           ilayesp(nlayesp) = klev
1518           rkmnz = radio+diz(klev)
1519           esp(nlayesp)  =  sqrt( rkmnz**2 - rkmmini**2 )       ! aux.factor[km]
1520           esp(nlayesp)  =  esp(nlayesp) * 1.d5                 ! aux.f. [cm]
1521           szalayesp(nlayesp) = asin( rkmmini/rkmnz )           ! [rad]
1522           szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592! [deg]
1523
1524
1525        ! Third case:  szadeg > 90
1526        ! The optical thickness is evaluated.
1527        ! We deal with 5 different regions:
1528        !   1: all layers between z and ztop ("upper part of ray")
1529        !   2: the layer at ztop ("uppermost layer")
1530        !   3: the lowest layer, at zmini
1531        !   4: the layers increasing from zmini to z (here SZA<90)
1532        !   5: the layers decreasing from z to zmini (here SZA>90)
1533        else if(szadeg.gt.90.d0) then
1534
1535           zmini=(radio+zz)*sin(szarad)-radio
1536           rkmmini = radio + zmini
1537
1538           if(zmini.lt.diz(1)) then         ! Can see the sun?  No => esp(j)=inft
1539             nlayesp = nlayesp + 1
1540             ilayesp(nlayesp) = - 1     ! Value to mark "no sun on view"
1541!             esp(nlayesp) = 1.e30
1542
1543           else
1544              jmin=grid_R8(zmini,diz,klev)+1
1545             
1546
1547              if(abs(zz-diz(klev)).lt.1.d-4) goto 9876
1548
1549              ! 1st Zone: Upper part of ray
1550              !
1551              do j=grid_R8(zz,diz,klev),klev-1
1552                nlayesp = nlayesp + 1
1553                ilayesp(nlayesp) = j
1554                esp(nlayesp) =
1555     $                sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) -
1556     $                sqrt( (radio+diz(j))**2 - rkmmini**2 )          ! [km]
1557                esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
1558                rkmj = radio+diz(j)
1559                szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
1560                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592      ! [deg]
1561              end do
1562
1563 9876         continue
1564              ! 2nd Zone:  Uppermost layer of ray.
1565              !
1566              nlayesp = nlayesp + 1
1567              ilayesp(nlayesp) = klev
1568              rkmnz = radio+diz(klev)
1569              esp(nlayesp) =  sqrt( rkmnz**2 - rkmmini**2 )      !aux.factor[km]
1570              esp(nlayesp) = esp(nlayesp) * 1.d5                 !aux.f.[cm]
1571              szalayesp(nlayesp) = asin( rkmmini/rkmnz )           ! [rad]
1572              szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
1573
1574              ! 3er Zone: Lowestmost layer of ray
1575              !
1576              if ( jmin .ge. 2 ) then      ! If above the planet's surface
1577                j=jmin-1
1578                nlayesp = nlayesp + 1
1579                ilayesp(nlayesp) = j
1580                esp(nlayesp) = 2. *
1581     $                 sqrt( (radio+diz(j+1))**2 -rkmmini**2 )       ! [km]
1582                esp(nlayesp) = esp(nlayesp) * 1.d5                   ! [cm]
1583                rkmj = radio+diz(j+1)
1584                szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad]
1585                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
1586              endif
1587
1588              ! 4th zone: Lower part of ray, increasing from zmin to z
1589              !    ( layers with SZA < 90 deg )
1590              do j=jmin,grid_R8(zz,diz,klev)-1
1591                nlayesp = nlayesp + 1
1592                ilayesp(nlayesp) = j
1593                esp(nlayesp) =
1594     $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
1595     $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )       !  [km]
1596                esp(nlayesp) = esp(nlayesp) * 1.d5                     ! [cm]
1597                rkmj = radio+diz(j)
1598                szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
1599                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
1600              end do
1601
1602              ! 5th zone: Lower part of ray, decreasing from z to zmin
1603              !    ( layers with SZA > 90 deg )
1604              do j=grid_R8(zz,diz,klev)-1, jmin, -1
1605                nlayesp = nlayesp + 1
1606                ilayesp(nlayesp) = j
1607                esp(nlayesp) =
1608     $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
1609     $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )        ! [km]
1610                esp(nlayesp) = esp(nlayesp) * 1.d5                      ! [cm]
1611                rkmj = radio+diz(j)
1612                szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj )          ! [rad]
1613                szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592 ! [deg]
1614              end do
1615
1616           end if
1617
1618        end if
1619
1620        return
1621
1622        end
1623
1624
1625c**********************************************************************
1626c***********************************************************************
1627
1628        function grid_R8 (z, zgrid, nz)
1629
1630c Returns the index where z is located within vector zgrid
1631c The vector zgrid must be monotonously increasing, otherwise program stops.
1632c If z is outside zgrid limits, or zgrid dimension is nz<2, the program stops.
1633c
1634c FGG     Aug-2004     Correct z.lt.zgrid(i) to .le.
1635c MALV    Jul-2003
1636c***********************************************************************
1637
1638        implicit none
1639
1640c Arguments
1641        integer   nz
1642        real*8      z
1643        real*8      zgrid(nz)
1644        integer   grid_R8
1645
1646c Local 
1647        integer   i, nz1, nznew
1648
1649c*** CODE START
1650
1651        if ( z .lt. zgrid(1) .or. z.gt.zgrid(nz) ) then
1652           write (*,*) ' GRID/ z outside bounds of zgrid '
1653           write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz)
1654           stop ' Serious error in GRID.F '
1655        endif
1656        if ( nz .lt. 2 ) then
1657           write (*,*) ' GRID/ zgrid needs 2 points at least ! '
1658           stop ' Serious error in GRID.F '
1659        endif
1660        if ( zgrid(1) .ge. zgrid(nz) ) then
1661           write (*,*) ' GRID/ zgrid must increase with index'
1662           stop ' Serious error in GRID.F '
1663        endif
1664
1665        nz1 = 1
1666        nznew = nz/2
1667        if ( z .gt. zgrid(nznew) ) then
1668           nz1 = nznew
1669           nznew = nz
1670        endif
1671        do i=nz1+1,nznew
1672           if ( z. eq. zgrid(i) ) then
1673              grid_R8=i
1674              return
1675              elseif ( z .le. zgrid(i) ) then
1676              grid_R8 = i-1
1677              return
1678           endif
1679        enddo
1680        grid_R8 = nz
1681        return
1682
1683        end
1684
1685
1686
1687!c***************************************************
1688!c***************************************************
1689
1690      subroutine flujo(date)
1691
1692
1693!c     fgg           nov 2002     first version
1694!c***************************************************
1695      use dimphy
1696      use conc
1697      implicit none
1698
1699
1700!     common variables and constants
1701      include "param.h"
1702      include 'param_v4.h'
1703      include "clesphys.h"
1704
1705!     Arguments
1706
1707       real date
1708c       integer, parameter :: dateyr = 2006
1709
1710!     Local variable and constants
1711!    dist_sol : distance venus - soleil
1712
1713      real, parameter :: dist_sol=0.72333
1714      integer i
1715      integer inter
1716      real    nada
1717
1718!c*************************************************
1719
1720      if(date.lt.1985.) date=1985.
1721      if(date.gt.2001.) date=2001.
1722     
1723      do i=1,ninter
1724         fluxtop(i)=1.
1725         !Variation of solar flux with 11 years solar cycle
1726         !For more details, see Gonzalez-Galindo et al. 2005
1727         !To be improved in next versions
1728        if(i.le.24 .and.solvarmod.eq.0) then
1729          fluxtop(i)=(((ct1(i)+p1(i)*date)/2.)                 
1730     $        *sin(2.*3.1416/11.*(date-1985.-3.1416))         
1731     $        +(ct2(i)+p2(i)*date)+1.)*fluxtop(i)
1732         
1733        end if
1734 !     The solar flux calculated 
1735 !      is  corrected for
1736 !     the actual Venus-Sun dist
1737        fluxtop(i)=fluxtop(i)*(1./dist_sol)**2
1738
1739
1740       end do
1741     
1742      return
1743      end
Note: See TracBrowser for help on using the repository browser.