source: trunk/LMDZ.VENUS/libf/phyvenus/jthermcalc_e107.F @ 2732

Last change on this file since 2732 was 2464, checked in by slebonnois, 4 years ago

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