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

Last change on this file since 2212 was 2212, checked in by aslmd, 5 years ago

MESOSCALE: precompiling flags 1) change in physiq_mod to avoid new GCM subgrid-scale slope wind param 2) addition in jthermcalc because abort_gcm is no longer available in mesoscale setting and this is about some lines related to pure testing.
File size: 61.2 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_gcm('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 wm(nlayer),wp(nlayer),p(nlayer)
1380      integer nm(nlayer)
1381      real*8 pin(nl)
1382      real*8 limup,limdown
1383      integer nl,nlayer,n1,n,np,nini
1384      nini=1
1385      do n1=1,nlayer
1386         if(p(n1) .gt. limup .or. p(n1) .lt. limdown) then
1387            wm(n1) = 0.d0
1388            wp(n1) = 0.d0
1389         else
1390            do n = nini,nl-1
1391               if (p(n1).ge.pin(n).and.p(n1).le.pin(n+1)) then
1392                  nm(n1)=n
1393                  np=n+1
1394                  wm(n1)=abs(pin(n)-p(n1))/(pin(np)-pin(n))
1395                  wp(n1)=1.d0 - wm(n1)
1396                  nini = n
1397                  exit
1398               endif
1399            enddo
1400         endif
1401      enddo
1402      return
1403      end
1404
1405
1406c**********************************************************************
1407c**********************************************************************
1408
1409      subroutine espesor_optico_A (ig,capa,nlayer, szadeg,z,
1410     @                   nz3,iz,esp,ilayesp,szalayesp,nlayesp, zmini)
1411
1412c     fgg              nov 03      Adaptation to Martian model
1413c     malv             jul 03      Corrected z grid. Split in alt & frec codes
1414c     fgg              feb 03      first version
1415*************************************************************************
1416
1417      use param_v4_h, only: radio
1418      implicit none
1419
1420c     arguments
1421
1422      real        szadeg                ! I. SZA [rad]
1423      real        z                     ! I. altitude of interest [km]
1424      integer     nz3,ig,nlayer         ! I. dimension of esp, ylayesp, etc...
1425                                        !  (=2*nlayer= max# of layers in ray path)
1426      real     iz(nlayer+1)              ! I. Altitude of each layer
1427      real*8        esp(nz3)            ! O. layer widths after geometrically
1428                                        !    amplified; in [cm] except at TOA
1429                                        !    where an auxiliary value is used
1430      real*8        ilayesp(nz3)        ! O. Indexes of layers along ray path
1431      real*8        szalayesp(nz3)      ! O. SZA [deg]    "     "       "
1432      integer       nlayesp
1433!      real*8        nlayesp             ! O. # layers along ray path at this z
1434      real*8        zmini               ! O. Minimum altitud of ray path [km]
1435
1436
1437c     local variables and constants
1438
1439        integer     j,i,capa
1440        integer     jmin                  ! index of min.altitude along ray path
1441        real*8      szarad                ! SZA [deg]
1442        real*8      zz
1443        real*8      diz(nlayer+1)
1444        real*8      rkmnz                 ! distance TOA to center of Planet [km]
1445        real*8      rkmmini               ! distance zmini to center of P [km]
1446        real*8      rkmj                  ! intermediate distance to C of P [km]
1447c external function
1448        external  grid_R8          ! Returns index of layer containing the altitude
1449                                ! of interest, z; for example, if
1450                                ! zkm(i)=z or zkm(i)<z<zkm(i+1) => grid(z)=i
1451        integer   grid_R8
1452
1453*************************************************************************     
1454        szarad = dble(szadeg)*3.141592d0/180.d0
1455        zz=dble(z)
1456        do i=1,nlayer
1457           diz(i)=dble(iz(i))
1458        enddo
1459        do j=1,nz3
1460          esp(j) = 0.d0
1461          szalayesp(j) = 777.d0
1462          ilayesp(j) = 0
1463        enddo
1464        nlayesp = 0
1465
1466        ! First case: szadeg<60
1467        ! The optical thickness will be given by  1/cos(sza)
1468        ! We deal with 2 different regions:
1469        !   1: First, all layers between z and ztop ("upper part of ray")
1470        !   2: Second, the layer at ztop
1471        if(szadeg.lt.60.d0) then
1472
1473           zmini = zz
1474           if(abs(zz-diz(nlayer)).lt.1.d-3) goto 1357
1475           ! 1st Zone: Upper part of ray
1476           !
1477           do j=grid_R8(zz,diz,nlayer),nlayer-1
1478             nlayesp = nlayesp + 1
1479             ilayesp(nlayesp) = j
1480             esp(nlayesp) = (diz(j+1)-diz(j)) / cos(szarad)        ! [km]
1481             esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
1482             szalayesp(nlayesp) = szadeg
1483           end do
1484
1485           !
1486           ! 2nd Zone: Top layer
1487 1357      continue
1488           nlayesp = nlayesp + 1
1489           ilayesp(nlayesp) = nlayer
1490           esp(nlayesp) = 1.d0 / cos(szarad)         ! aux. non-dimens. factor
1491           szalayesp(nlayesp) = szadeg
1492
1493
1494        ! Second case:  60 < szadeg < 90
1495        ! The optical thickness is evaluated.
1496        !    (the magnitude of the effect of not using cos(sza) is 3.e-5
1497        !     for z=60km & sza=30, and 5e-4 for z=60km & sza=60, approximately)
1498        ! We deal with 2 different regions:
1499        !   1: First, all layers between z and ztop ("upper part of ray")
1500        !   2: Second, the layer at ztop ("uppermost layer")
1501        else if(szadeg.le.90.d0.and.szadeg.ge.60.d0) then
1502
1503           zmini=(radio+zz)*sin(szarad)-radio
1504           rkmmini = radio + zmini
1505
1506           if(abs(zz-diz(nlayer)).lt.1.d-4) goto 1470
1507
1508           ! 1st Zone: Upper part of ray
1509           !
1510           do j=grid_R8(zz,diz,nlayer),nlayer-1
1511              nlayesp = nlayesp + 1
1512              ilayesp(nlayesp) = j
1513              esp(nlayesp) =
1514     #             sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) -
1515     #             sqrt( (radio+diz(j))**2 - rkmmini**2 )           ! [km]
1516              esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
1517              rkmj = radio+diz(j)
1518              szalayesp(nlayesp) = asin( rkmmini/rkmj )             ! [rad]
1519              szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592 ! [deg]
1520           end do
1521 1470      continue
1522           ! 2nd Zone:  Uppermost layer of ray.
1523           !
1524           nlayesp = nlayesp + 1
1525           ilayesp(nlayesp) = nlayer
1526           rkmnz = radio+diz(nlayer)
1527           esp(nlayesp)  =  sqrt( rkmnz**2 - rkmmini**2 )       ! aux.factor[km]
1528           esp(nlayesp)  =  esp(nlayesp) * 1.d5                 ! aux.f. [cm]
1529           szalayesp(nlayesp) = asin( rkmmini/rkmnz )           ! [rad]
1530           szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592! [deg]
1531
1532
1533        ! Third case:  szadeg > 90
1534        ! The optical thickness is evaluated.
1535        ! We deal with 5 different regions:
1536        !   1: all layers between z and ztop ("upper part of ray")
1537        !   2: the layer at ztop ("uppermost layer")
1538        !   3: the lowest layer, at zmini
1539        !   4: the layers increasing from zmini to z (here SZA<90)
1540        !   5: the layers decreasing from z to zmini (here SZA>90)
1541        else if(szadeg.gt.90.d0) then
1542
1543           zmini=(radio+zz)*sin(szarad)-radio
1544           !zmini should be lower than zz, as SZA<90. However, in situations
1545           !where SZA is very close to 90, rounding errors can make zmini
1546           !slightly higher than zz, causing problems in the determination
1547           !of the jmin index. A correction is implemented in the determination
1548           !of jmin, some lines below
1549           rkmmini = radio + zmini
1550
1551           if(zmini.lt.diz(1)) then         ! Can see the sun?  No => esp(j)=inft
1552             nlayesp = nlayesp + 1
1553             ilayesp(nlayesp) = - 1     ! Value to mark "no sun on view"
1554!             esp(nlayesp) = 1.e30
1555
1556           else
1557              jmin=grid_R8(zmini,diz,nlayer)+1
1558              !Correction for possible rounding errors when SZA very close
1559              !to 90 degrees
1560              if(jmin.gt.grid_R8(zz,diz,nlayer)) then
1561                 write(*,*)'jthermcalc warning: possible rounding error'
1562                 write(*,*)'point,sza,layer:',ig,szadeg,capa
1563                 jmin=grid_R8(zz,diz,nlayer)
1564              endif
1565
1566              if(abs(zz-diz(nlayer)).lt.1.d-4) goto 9876
1567
1568              ! 1st Zone: Upper part of ray
1569              !
1570              do j=grid_R8(zz,diz,nlayer),nlayer-1
1571                nlayesp = nlayesp + 1
1572                ilayesp(nlayesp) = j
1573                esp(nlayesp) =
1574     $                sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) -
1575     $                sqrt( (radio+diz(j))**2 - rkmmini**2 )          ! [km]
1576                esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
1577                rkmj = radio+diz(j)
1578                szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
1579                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592      ! [deg]
1580              end do
1581
1582 9876         continue
1583              ! 2nd Zone:  Uppermost layer of ray.
1584              !
1585              nlayesp = nlayesp + 1
1586              ilayesp(nlayesp) = nlayer
1587              rkmnz = radio+diz(nlayer)
1588              esp(nlayesp) =  sqrt( rkmnz**2 - rkmmini**2 )      !aux.factor[km]
1589              esp(nlayesp) = esp(nlayesp) * 1.d5                 !aux.f.[cm]
1590              szalayesp(nlayesp) = asin( rkmmini/rkmnz )           ! [rad]
1591              szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
1592
1593              ! 3er Zone: Lowestmost layer of ray
1594              !
1595              if ( jmin .ge. 2 ) then      ! If above the planet's surface
1596                j=jmin-1
1597                nlayesp = nlayesp + 1
1598                ilayesp(nlayesp) = j
1599                esp(nlayesp) = 2. *
1600     $                 sqrt( (radio+diz(j+1))**2 -rkmmini**2 )       ! [km]
1601                esp(nlayesp) = esp(nlayesp) * 1.d5                   ! [cm]
1602                rkmj = radio+diz(j+1)
1603                szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad]
1604                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
1605              endif
1606
1607              ! 4th zone: Lower part of ray, increasing from zmin to z
1608              !    ( layers with SZA < 90 deg )
1609              do j=jmin,grid_R8(zz,diz,nlayer)-1
1610                nlayesp = nlayesp + 1
1611                ilayesp(nlayesp) = j
1612                esp(nlayesp) =
1613     $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
1614     $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )       ! [km]
1615                esp(nlayesp) = esp(nlayesp) * 1.d5                     ! [cm]
1616                rkmj = radio+diz(j)
1617                szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
1618                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
1619              end do
1620
1621              ! 5th zone: Lower part of ray, decreasing from z to zmin
1622              !    ( layers with SZA > 90 deg )
1623              do j=grid_R8(zz,diz,nlayer)-1, jmin, -1
1624                nlayesp = nlayesp + 1
1625                ilayesp(nlayesp) = j
1626                esp(nlayesp) =
1627     $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
1628     $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )        ! [km]
1629                esp(nlayesp) = esp(nlayesp) * 1.d5                      ! [cm]
1630                rkmj = radio+diz(j)
1631                szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj )          ! [rad]
1632                szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592 ! [deg]
1633              end do
1634
1635           end if
1636
1637        end if
1638
1639
1640        return
1641
1642        end
1643
1644
1645
1646c**********************************************************************
1647c***********************************************************************
1648
1649        function grid_R8 (z, zgrid, nz)
1650
1651c Returns the index where z is located within vector zgrid
1652c The vector zgrid must be monotonously increasing, otherwise program stops.
1653c If z is outside zgrid limits, or zgrid dimension is nz<2, the program stops.
1654c
1655c FGG     Aug-2004     Correct z.lt.zgrid(i) to .le.
1656c MALV    Jul-2003
1657c***********************************************************************
1658
1659        implicit none
1660
1661c Arguments
1662        integer   nz
1663        real*8      z
1664        real*8      zgrid(nz)
1665        integer   grid_R8
1666
1667c Local 
1668        integer   i, nz1, nznew
1669
1670c*** CODE START
1671
1672        if ( z .lt. zgrid(1)  ) then
1673           write (*,*) ' GRID/ z outside bounds of zgrid '
1674           write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz)
1675           z = zgrid(1)
1676           write(*,*) 'WARNING: error in grid_r8 (jthermcalc.F)'
1677           write(*,*) 'Please check values of z and zgrid above'
1678        endif
1679        if (z .gt. zgrid(nz) ) then
1680           write (*,*) ' GRID/ z outside bounds of zgrid '
1681           write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz)
1682           z = zgrid(nz)
1683           write(*,*) 'WARNING: error in grid_r8 (jthermcalc.F)'
1684           write(*,*) 'Please check values of z and zgrid above'
1685        endif
1686        if ( nz .lt. 2 ) then
1687           write (*,*) ' GRID/ zgrid needs 2 points at least ! '
1688           stop ' Serious error in GRID.F '
1689        endif
1690        if ( zgrid(1) .ge. zgrid(nz) ) then
1691           write (*,*) ' GRID/ zgrid must increase with index'
1692           stop ' Serious error in GRID.F '
1693        endif
1694
1695        nz1 = 1
1696        nznew = nz/2
1697        if ( z .gt. zgrid(nznew) ) then
1698           nz1 = nznew
1699           nznew = nz
1700        endif
1701        do i=nz1+1,nznew
1702           if ( z. eq. zgrid(i) ) then
1703              grid_R8=i
1704              return
1705              elseif ( z .le. zgrid(i) ) then
1706              grid_R8 = i-1
1707              return
1708           endif
1709        enddo
1710        grid_R8 = nz
1711        return
1712
1713
1714
1715        end
1716
1717
1718
1719!c***************************************************
1720!c***************************************************
1721
1722      subroutine flujo(date)
1723
1724
1725!c     fgg           nov 2002     first version
1726!c***************************************************
1727
1728      use comsaison_h, only: dist_sol
1729      use param_v4_h, only: ninter,
1730     .                      fluxtop, ct1, ct2, p1, p2
1731      implicit none
1732
1733
1734!     common variables and constants
1735      include "callkeys.h"
1736
1737
1738!     Arguments
1739
1740      real date
1741
1742
1743!     Local variable and constants
1744
1745      integer i
1746      integer inter
1747      real    nada
1748
1749!c*************************************************
1750
1751      if(date.lt.1985.) date=1985.
1752      if(date.gt.2001.) date=2001.
1753     
1754      do i=1,ninter
1755         fluxtop(i)=1.
1756         !Variation of solar flux with 11 years solar cycle
1757         !For more details, see Gonzalez-Galindo et al. 2005
1758         !To be improved in next versions
1759        if(i.le.24.and.solvarmod.eq.0) then
1760            fluxtop(i)=(((ct1(i)+p1(i)*date)/2.)                 
1761     $           *sin(2.*3.1416/11.*(date-1985.-3.1416))         
1762     $           +(ct2(i)+p2(i)*date)+1.)*fluxtop(i)
1763         end if
1764         fluxtop(i)=fluxtop(i)*(1.52/dist_sol)**2
1765      end do
1766     
1767      return
1768      end
Note: See TracBrowser for help on using the repository browser.