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

Last change on this file since 1356 was 1310, checked in by slebonnois, 11 years ago

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