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

Last change on this file since 2042 was 2042, checked in by emillour, 6 years ago

Mars GCM:
Modifications to use the parametrized photoabsorbtion coefficients;
a first step towards implementing ionospheric chemistry in the new
chemical solver:

  • change in species indexes in chemthermos.F90, paramfoto_compact.F, hrtherm.F and euvheat.F90
  • calchim.F90: added a variable in call to photochemistry
  • photochemistry.F90: added calls to jthermcalc_e107 and phdisrate, with an additionlal flag, jparam (.false. by default). The computed photodissociation coefficents are sent to v_phot, which is used in the chemistry. Thus concentrations computed in chimtogcm are now done over all atmospheric layers.

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