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

Last change on this file since 2883 was 2674, checked in by emillour, 3 years ago

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