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

Last change on this file since 1443 was 1442, checked in by slebonnois, 10 years ago

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