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

Last change on this file since 1775 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

  • Remove all reference/use of nlayermx and dimphys.h
  • Made use of automatic arrays whenever arrays are needed with dimension nlayer
  • Remove lots of obsolete reference to dimensions.h
  • Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

  • Deleted param.h and put contents into module param_v4_h
  • Adapted testphys1d, newstart, etc...
  • Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

  • Optimized computations in paramfoto_compact (twice less dlog10 calculations)
  • Checked consistency before/after modification in debug mode
  • Checked performance is not impacted (same as before)
File size: 60.6 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_o2=2
1132      integer,parameter :: i_o=3
1133      integer,parameter :: i_co=4
1134      integer,parameter :: i_h=5
1135      integer,parameter :: i_h2=8
1136      integer,parameter :: i_h2o=9
1137      integer,parameter :: i_h2o2=10
1138      integer,parameter :: i_o3=12
1139      integer,parameter :: i_n2=13
1140      integer,parameter :: i_n=14
1141      integer,parameter :: i_no=15
1142      integer,parameter :: i_no2=17
1143
1144
1145c*************************PROGRAM STARTS*******************************
1146
1147      nz3 = nlayer*2
1148      do i=1,nlayer
1149         xx = ( radio + iz(i) ) * 1.e5
1150         grav(i) = gg * masa /(xx**2)
1151      end do
1152
1153      !Scale heights
1154      xx = kboltzman * tx(nlayer) * n_avog / grav(nlayer)
1155      Ho3p  = xx / mmol(igcm_o)
1156      Hco2  = xx / mmol(igcm_co2)
1157      Ho2   = xx / mmol(igcm_o2)
1158      Hh2   = xx / mmol(igcm_h2)
1159      Hh2o  = xx / mmol(igcm_h2o_vap)
1160      Hh2o2 = xx / mmol(igcm_h2o2)
1161      Hco   = xx / mmol(igcm_co)
1162      Hh    = xx / mmol(igcm_h)
1163      !Only if O3 chem. required
1164      if(chemthermod.ge.1)
1165     $     Ho3   = xx / mmol(igcm_o3)
1166      !Only if N or ion chem.
1167      if(chemthermod.ge.2) then
1168         Hn2   = xx / mmol(igcm_n2)
1169         Hn    = xx / mmol(igcm_n)
1170         Hno   = xx / mmol(igcm_no)
1171         Hno2  = xx / mmol(igcm_no2)
1172      endif
1173      ! first loop in altitude : initialisation
1174      do i=nlayer,1,-1
1175         !Column initialisation
1176         co2colx(i)  = 0.
1177         o2colx(i)   = 0.
1178         o3pcolx(i)  = 0.
1179         h2colx(i)   = 0.
1180         h2ocolx(i)  = 0.
1181         h2o2colx(i) = 0.
1182         o3colx(i)   = 0.
1183         n2colx(i)   = 0.
1184         ncolx(i)    = 0.
1185         nocolx(i)   = 0.
1186         cocolx(i)   = 0.
1187         hcolx(i)    = 0.
1188         no2colx(i)  = 0.
1189         !Densities
1190         co2x(i)  = rm(i,i_co2)
1191         o2x(i)   = rm(i,i_o2)
1192         o3px(i)  = rm(i,i_o)
1193         h2x(i)   = rm(i,i_h2)
1194         h2ox(i)  = rm(i,i_h2o)
1195         h2o2x(i) = rm(i,i_h2o2)
1196         cox(i)   = rm(i,i_co)
1197         hx(i)    = rm(i,i_h)
1198         !Only if O3 chem. required
1199         if(chemthermod.ge.1)
1200     $        o3x(i)   = rm(i,i_o3)
1201         !Only if Nitrogen of ion chemistry requested
1202         if(chemthermod.ge.2) then
1203            n2x(i)   = rm(i,i_n2)
1204            nx(i)    = rm(i,i_n)
1205            nox(i)   = rm(i,i_no)
1206            no2x(i)  = rm(i,i_no2)
1207         endif
1208      enddo
1209      ! second loop in altitude : column calculations
1210      do i=nlayer,1,-1
1211         !Routine to calculate the geometrical length of each layer
1212         call espesor_optico_A(ig,i,nlayer,zenit,iz(i),nz3,iz,esp,
1213     $         ilayesp,szalayesp,nlayesp, zmini)
1214         if(ilayesp(nlayesp).eq.-1) then
1215            co2colx(i)=1.e25
1216            o2colx(i)=1.e25
1217            o3pcolx(i)=1.e25
1218            h2colx(i)=1.e25
1219            h2ocolx(i)=1.e25
1220            h2o2colx(i)=1.e25
1221            o3colx(i)=1.e25
1222            n2colx(i)=1.e25
1223            ncolx(i)=1.e25
1224            nocolx(i)=1.e25
1225            cocolx(i)=1.e25
1226            hcolx(i)=1.e25
1227            no2colx(i)=1.e25
1228         else
1229            rcmnz = ( radio + iz(nlayer) ) * 1.e5
1230            rcmmini = ( radio + zmini ) * 1.e5
1231            !Column calculation taking into account the geometrical depth
1232            !calculated before
1233            do j=1,nlayesp
1234               jj=ilayesp(j)
1235               !Top layer
1236               if(jj.eq.nlayer) then
1237                  if(zenit.le.60.) then
1238                     o3pcolx(i)=o3pcolx(i)+o3px(nlayer)*Ho3p*esp(j)
1239     $                    *1.e-5
1240                     co2colx(i)=co2colx(i)+co2x(nlayer)*Hco2*esp(j)
1241     $                    *1.e-5
1242                     h2o2colx(i)=h2o2colx(i)+
1243     $                    h2o2x(nlayer)*Hh2o2*esp(j)*1.e-5
1244                     o2colx(i)=o2colx(i)+o2x(nlayer)*Ho2*esp(j)
1245     $                    *1.e-5
1246                     h2colx(i)=h2colx(i)+h2x(nlayer)*Hh2*esp(j)
1247     $                    *1.e-5
1248                     h2ocolx(i)=h2ocolx(i)+h2ox(nlayer)*Hh2o*esp(j)
1249     $                    *1.e-5                     
1250                     cocolx(i)=cocolx(i)+cox(nlayer)*Hco*esp(j)
1251     $                    *1.e-5
1252                     hcolx(i)=hcolx(i)+hx(nlayer)*Hh*esp(j)
1253     $                    *1.e-5
1254                     !Only if O3 chemistry required
1255                     if(chemthermod.ge.1) o3colx(i)=
1256     $                    o3colx(i)+o3x(nlayer)*Ho3*esp(j)
1257     $                    *1.e-5
1258                     !Only if N or ion chemistry requested
1259                     if(chemthermod.ge.2) then
1260                        n2colx(i)=n2colx(i)+n2x(nlayer)*Hn2*esp(j)
1261     $                    *1.e-5
1262                        ncolx(i)=ncolx(i)+nx(nlayer)*Hn*esp(j)
1263     $                       *1.e-5
1264                        nocolx(i)=nocolx(i)+nox(nlayer)*Hno*esp(j)
1265     $                       *1.e-5
1266                        no2colx(i)=no2colx(i)+no2x(nlayer)*Hno2*esp(j)
1267     $                       *1.e-5
1268                     endif
1269                  else if(zenit.gt.60.) then
1270                     espco2 =sqrt((rcmnz+Hco2)**2 -rcmmini**2) - esp(j)
1271                     espo2  = sqrt((rcmnz+Ho2)**2 -rcmmini**2) - esp(j)
1272                     espo3p = sqrt((rcmnz+Ho3p)**2 -rcmmini**2)- esp(j)
1273                     esph2  = sqrt((rcmnz+Hh2)**2 -rcmmini**2) - esp(j)
1274                     esph2o = sqrt((rcmnz+Hh2o)**2 -rcmmini**2)- esp(j)
1275                     esph2o2= sqrt((rcmnz+Hh2o2)**2-rcmmini**2)- esp(j)
1276                     espco  = sqrt((rcmnz+Hco)**2 -rcmmini**2) - esp(j)
1277                     esph   = sqrt((rcmnz+Hh)**2 -rcmmini**2)  - esp(j)
1278                     !Only if O3 chemistry required
1279                     if(chemthermod.ge.1)                     
1280     $                   espo3=sqrt((rcmnz+Ho3)**2-rcmmini**2)-esp(j)
1281                     !Only if N or ion chemistry requested
1282                     if(chemthermod.ge.2) then
1283                        espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j)
1284                        espn  =sqrt((rcmnz+Hn)**2-rcmmini**2)  - esp(j)
1285                        espno =sqrt((rcmnz+Hno)**2-rcmmini**2) - esp(j)
1286                        espno2=sqrt((rcmnz+Hno2)**2-rcmmini**2)- esp(j)
1287                     endif
1288                     co2colx(i) = co2colx(i) + espco2*co2x(nlayer)
1289                     o2colx(i)  = o2colx(i)  + espo2*o2x(nlayer)
1290                     o3pcolx(i) = o3pcolx(i) + espo3p*o3px(nlayer)
1291                     h2colx(i)  = h2colx(i)  + esph2*h2x(nlayer)
1292                     h2ocolx(i) = h2ocolx(i) + esph2o*h2ox(nlayer)
1293                     h2o2colx(i)= h2o2colx(i)+ esph2o2*h2o2x(nlayer)
1294                     cocolx(i)  = cocolx(i)  + espco*cox(nlayer)
1295                     hcolx(i)   = hcolx(i)   + esph*hx(nlayer)
1296                     !Only if O3 chemistry required
1297                     if(chemthermod.ge.1)                     
1298     $                  o3colx(i) = o3colx(i)  + espo3*o3x(nlayer)
1299                     !Only if N or ion chemistry requested
1300                     if(chemthermod.ge.2) then
1301                        n2colx(i)  = n2colx(i)  + espn2*n2x(nlayer)
1302                        ncolx(i)   = ncolx(i)   + espn*nx(nlayer)
1303                        nocolx(i)  = nocolx(i)  + espno*nox(nlayer)
1304                        no2colx(i) = no2colx(i) + espno2*no2x(nlayer)
1305                     endif
1306                  endif !Of if zenit.lt.60
1307               !Other layers
1308               else
1309                  co2colx(i)  = co2colx(i) +
1310     $                 esp(j) * (co2x(jj)+co2x(jj+1)) / 2.
1311                  o2colx(i)   = o2colx(i) +
1312     $                 esp(j) * (o2x(jj)+o2x(jj+1)) / 2.
1313                  o3pcolx(i)  = o3pcolx(i) +
1314     $                 esp(j) * (o3px(jj)+o3px(jj+1)) / 2.
1315                  h2colx(i)   = h2colx(i) +
1316     $                 esp(j) * (h2x(jj)+h2x(jj+1)) / 2.
1317                  h2ocolx(i)  = h2ocolx(i) +
1318     $                 esp(j) * (h2ox(jj)+h2ox(jj+1)) / 2.
1319                  h2o2colx(i) = h2o2colx(i) +
1320     $                 esp(j) * (h2o2x(jj)+h2o2x(jj+1)) / 2.
1321                  cocolx(i)   = cocolx(i) +
1322     $                 esp(j) * (cox(jj)+cox(jj+1)) / 2.
1323                  hcolx(i)    = hcolx(i) +
1324     $                 esp(j) * (hx(jj)+hx(jj+1)) / 2.
1325                  !Only if O3 chemistry required
1326                  if(chemthermod.ge.1)
1327     $                 o3colx(i) = o3colx(i) +
1328     $                 esp(j) * (o3x(jj)+o3x(jj+1)) / 2.
1329                  !Only if N or ion chemistry requested
1330                  if(chemthermod.ge.2) then
1331                     n2colx(i)   = n2colx(i) +
1332     $                 esp(j) * (n2x(jj)+n2x(jj+1)) / 2.
1333                     ncolx(i)    = ncolx(i) +
1334     $                    esp(j) * (nx(jj)+nx(jj+1)) / 2.
1335                     nocolx(i)   = nocolx(i) +
1336     $                    esp(j) * (nox(jj)+nox(jj+1)) / 2.
1337                     no2colx(i)  = no2colx(i) +
1338     $                    esp(j) * (no2x(jj)+no2x(jj+1)) / 2.
1339                  endif
1340               endif  !Of if jj.eq.nlayer
1341            end do    !Of do j=1,nlayesp
1342         endif        !Of ilayesp(nlayesp).eq.-1
1343      enddo           !Of do i=nlayer,1,-1
1344
1345
1346      return
1347
1348
1349      end
1350
1351
1352c**********************************************************************
1353c**********************************************************************
1354
1355      subroutine interfast(wm,wp,nm,p,nlayer,pin,nl,limdown,limup)
1356C
1357C subroutine to perform linear interpolation in pressure from 1D profile
1358C escin(nl) sampled on pressure grid pin(nl) to profile
1359C escout(nlayer) on pressure grid p(nlayer).
1360C
1361      real*8 wm(nlayer),wp(nlayer),p(nlayer)
1362      integer nm(nlayer)
1363      real*8 pin(nl)
1364      real*8 limup,limdown
1365      integer nl,nlayer,n1,n,np,nini
1366      nini=1
1367      do n1=1,nlayer
1368         if(p(n1) .gt. limup .or. p(n1) .lt. limdown) then
1369            wm(n1) = 0.d0
1370            wp(n1) = 0.d0
1371         else
1372            do n = nini,nl-1
1373               if (p(n1).ge.pin(n).and.p(n1).le.pin(n+1)) then
1374                  nm(n1)=n
1375                  np=n+1
1376                  wm(n1)=abs(pin(n)-p(n1))/(pin(np)-pin(n))
1377                  wp(n1)=1.d0 - wm(n1)
1378                  nini = n
1379                  exit
1380               endif
1381            enddo
1382         endif
1383      enddo
1384      return
1385      end
1386
1387
1388c**********************************************************************
1389c**********************************************************************
1390
1391      subroutine espesor_optico_A (ig,capa,nlayer, szadeg,z,
1392     @                   nz3,iz,esp,ilayesp,szalayesp,nlayesp, zmini)
1393
1394c     fgg              nov 03      Adaptation to Martian model
1395c     malv             jul 03      Corrected z grid. Split in alt & frec codes
1396c     fgg              feb 03      first version
1397*************************************************************************
1398
1399      use param_v4_h, only: radio
1400      implicit none
1401
1402c     arguments
1403
1404      real        szadeg                ! I. SZA [rad]
1405      real        z                     ! I. altitude of interest [km]
1406      integer     nz3,ig,nlayer         ! I. dimension of esp, ylayesp, etc...
1407                                        !  (=2*nlayer= max# of layers in ray path)
1408      real     iz(nlayer+1)              ! I. Altitude of each layer
1409      real*8        esp(nz3)            ! O. layer widths after geometrically
1410                                        !    amplified; in [cm] except at TOA
1411                                        !    where an auxiliary value is used
1412      real*8        ilayesp(nz3)        ! O. Indexes of layers along ray path
1413      real*8        szalayesp(nz3)      ! O. SZA [deg]    "     "       "
1414      integer       nlayesp
1415!      real*8        nlayesp             ! O. # layers along ray path at this z
1416      real*8        zmini               ! O. Minimum altitud of ray path [km]
1417
1418
1419c     local variables and constants
1420
1421        integer     j,i,capa
1422        integer     jmin                  ! index of min.altitude along ray path
1423        real*8      szarad                ! SZA [deg]
1424        real*8      zz
1425        real*8      diz(nlayer+1)
1426        real*8      rkmnz                 ! distance TOA to center of Planet [km]
1427        real*8      rkmmini               ! distance zmini to center of P [km]
1428        real*8      rkmj                  ! intermediate distance to C of P [km]
1429c external function
1430        external  grid_R8          ! Returns index of layer containing the altitude
1431                                ! of interest, z; for example, if
1432                                ! zkm(i)=z or zkm(i)<z<zkm(i+1) => grid(z)=i
1433        integer   grid_R8
1434
1435*************************************************************************     
1436        szarad = dble(szadeg)*3.141592d0/180.d0
1437        zz=dble(z)
1438        do i=1,nlayer
1439           diz(i)=dble(iz(i))
1440        enddo
1441        do j=1,nz3
1442          esp(j) = 0.d0
1443          szalayesp(j) = 777.d0
1444          ilayesp(j) = 0
1445        enddo
1446        nlayesp = 0
1447
1448        ! First case: szadeg<60
1449        ! The optical thickness will be given by  1/cos(sza)
1450        ! We deal with 2 different regions:
1451        !   1: First, all layers between z and ztop ("upper part of ray")
1452        !   2: Second, the layer at ztop
1453        if(szadeg.lt.60.d0) then
1454
1455           zmini = zz
1456           if(abs(zz-diz(nlayer)).lt.1.d-3) goto 1357
1457           ! 1st Zone: Upper part of ray
1458           !
1459           do j=grid_R8(zz,diz,nlayer),nlayer-1
1460             nlayesp = nlayesp + 1
1461             ilayesp(nlayesp) = j
1462             esp(nlayesp) = (diz(j+1)-diz(j)) / cos(szarad)        ! [km]
1463             esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
1464             szalayesp(nlayesp) = szadeg
1465           end do
1466
1467           !
1468           ! 2nd Zone: Top layer
1469 1357      continue
1470           nlayesp = nlayesp + 1
1471           ilayesp(nlayesp) = nlayer
1472           esp(nlayesp) = 1.d0 / cos(szarad)         ! aux. non-dimens. factor
1473           szalayesp(nlayesp) = szadeg
1474
1475
1476        ! Second case:  60 < szadeg < 90
1477        ! The optical thickness is evaluated.
1478        !    (the magnitude of the effect of not using cos(sza) is 3.e-5
1479        !     for z=60km & sza=30, and 5e-4 for z=60km & sza=60, approximately)
1480        ! We deal with 2 different regions:
1481        !   1: First, all layers between z and ztop ("upper part of ray")
1482        !   2: Second, the layer at ztop ("uppermost layer")
1483        else if(szadeg.le.90.d0.and.szadeg.ge.60.d0) then
1484
1485           zmini=(radio+zz)*sin(szarad)-radio
1486           rkmmini = radio + zmini
1487
1488           if(abs(zz-diz(nlayer)).lt.1.d-4) goto 1470
1489
1490           ! 1st Zone: Upper part of ray
1491           !
1492           do j=grid_R8(zz,diz,nlayer),nlayer-1
1493              nlayesp = nlayesp + 1
1494              ilayesp(nlayesp) = j
1495              esp(nlayesp) =
1496     #             sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) -
1497     #             sqrt( (radio+diz(j))**2 - rkmmini**2 )           ! [km]
1498              esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
1499              rkmj = radio+diz(j)
1500              szalayesp(nlayesp) = asin( rkmmini/rkmj )             ! [rad]
1501              szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592 ! [deg]
1502           end do
1503 1470      continue
1504           ! 2nd Zone:  Uppermost layer of ray.
1505           !
1506           nlayesp = nlayesp + 1
1507           ilayesp(nlayesp) = nlayer
1508           rkmnz = radio+diz(nlayer)
1509           esp(nlayesp)  =  sqrt( rkmnz**2 - rkmmini**2 )       ! aux.factor[km]
1510           esp(nlayesp)  =  esp(nlayesp) * 1.d5                 ! aux.f. [cm]
1511           szalayesp(nlayesp) = asin( rkmmini/rkmnz )           ! [rad]
1512           szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592! [deg]
1513
1514
1515        ! Third case:  szadeg > 90
1516        ! The optical thickness is evaluated.
1517        ! We deal with 5 different regions:
1518        !   1: all layers between z and ztop ("upper part of ray")
1519        !   2: the layer at ztop ("uppermost layer")
1520        !   3: the lowest layer, at zmini
1521        !   4: the layers increasing from zmini to z (here SZA<90)
1522        !   5: the layers decreasing from z to zmini (here SZA>90)
1523        else if(szadeg.gt.90.d0) then
1524
1525           zmini=(radio+zz)*sin(szarad)-radio
1526           !zmini should be lower than zz, as SZA<90. However, in situations
1527           !where SZA is very close to 90, rounding errors can make zmini
1528           !slightly higher than zz, causing problems in the determination
1529           !of the jmin index. A correction is implemented in the determination
1530           !of jmin, some lines below
1531           rkmmini = radio + zmini
1532
1533           if(zmini.lt.diz(1)) then         ! Can see the sun?  No => esp(j)=inft
1534             nlayesp = nlayesp + 1
1535             ilayesp(nlayesp) = - 1     ! Value to mark "no sun on view"
1536!             esp(nlayesp) = 1.e30
1537
1538           else
1539              jmin=grid_R8(zmini,diz,nlayer)+1
1540              !Correction for possible rounding errors when SZA very close
1541              !to 90 degrees
1542              if(jmin.gt.grid_R8(zz,diz,nlayer)) then
1543                 write(*,*)'jthermcalc warning: possible rounding error'
1544                 write(*,*)'point,sza,layer:',ig,szadeg,capa
1545                 jmin=grid_R8(zz,diz,nlayer)
1546              endif
1547
1548              if(abs(zz-diz(nlayer)).lt.1.d-4) goto 9876
1549
1550              ! 1st Zone: Upper part of ray
1551              !
1552              do j=grid_R8(zz,diz,nlayer),nlayer-1
1553                nlayesp = nlayesp + 1
1554                ilayesp(nlayesp) = j
1555                esp(nlayesp) =
1556     $                sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) -
1557     $                sqrt( (radio+diz(j))**2 - rkmmini**2 )          ! [km]
1558                esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
1559                rkmj = radio+diz(j)
1560                szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
1561                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592      ! [deg]
1562              end do
1563
1564 9876         continue
1565              ! 2nd Zone:  Uppermost layer of ray.
1566              !
1567              nlayesp = nlayesp + 1
1568              ilayesp(nlayesp) = nlayer
1569              rkmnz = radio+diz(nlayer)
1570              esp(nlayesp) =  sqrt( rkmnz**2 - rkmmini**2 )      !aux.factor[km]
1571              esp(nlayesp) = esp(nlayesp) * 1.d5                 !aux.f.[cm]
1572              szalayesp(nlayesp) = asin( rkmmini/rkmnz )           ! [rad]
1573              szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
1574
1575              ! 3er Zone: Lowestmost layer of ray
1576              !
1577              if ( jmin .ge. 2 ) then      ! If above the planet's surface
1578                j=jmin-1
1579                nlayesp = nlayesp + 1
1580                ilayesp(nlayesp) = j
1581                esp(nlayesp) = 2. *
1582     $                 sqrt( (radio+diz(j+1))**2 -rkmmini**2 )       ! [km]
1583                esp(nlayesp) = esp(nlayesp) * 1.d5                   ! [cm]
1584                rkmj = radio+diz(j+1)
1585                szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad]
1586                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
1587              endif
1588
1589              ! 4th zone: Lower part of ray, increasing from zmin to z
1590              !    ( layers with SZA < 90 deg )
1591              do j=jmin,grid_R8(zz,diz,nlayer)-1
1592                nlayesp = nlayesp + 1
1593                ilayesp(nlayesp) = j
1594                esp(nlayesp) =
1595     $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
1596     $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )       ! [km]
1597                esp(nlayesp) = esp(nlayesp) * 1.d5                     ! [cm]
1598                rkmj = radio+diz(j)
1599                szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
1600                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
1601              end do
1602
1603              ! 5th zone: Lower part of ray, decreasing from z to zmin
1604              !    ( layers with SZA > 90 deg )
1605              do j=grid_R8(zz,diz,nlayer)-1, jmin, -1
1606                nlayesp = nlayesp + 1
1607                ilayesp(nlayesp) = j
1608                esp(nlayesp) =
1609     $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
1610     $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )        ! [km]
1611                esp(nlayesp) = esp(nlayesp) * 1.d5                      ! [cm]
1612                rkmj = radio+diz(j)
1613                szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj )          ! [rad]
1614                szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592 ! [deg]
1615              end do
1616
1617           end if
1618
1619        end if
1620
1621
1622        return
1623
1624        end
1625
1626
1627
1628c**********************************************************************
1629c***********************************************************************
1630
1631        function grid_R8 (z, zgrid, nz)
1632
1633c Returns the index where z is located within vector zgrid
1634c The vector zgrid must be monotonously increasing, otherwise program stops.
1635c If z is outside zgrid limits, or zgrid dimension is nz<2, the program stops.
1636c
1637c FGG     Aug-2004     Correct z.lt.zgrid(i) to .le.
1638c MALV    Jul-2003
1639c***********************************************************************
1640
1641        implicit none
1642
1643c Arguments
1644        integer   nz
1645        real*8      z
1646        real*8      zgrid(nz)
1647        integer   grid_R8
1648
1649c Local 
1650        integer   i, nz1, nznew
1651
1652c*** CODE START
1653
1654        if ( z .lt. zgrid(1)  ) then
1655           write (*,*) ' GRID/ z outside bounds of zgrid '
1656           write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz)
1657           z = zgrid(1)
1658           write(*,*) 'WARNING: error in grid_r8 (jthermcalc.F)'
1659           write(*,*) 'Please check values of z and zgrid above'
1660        endif
1661        if (z .gt. zgrid(nz) ) then
1662           write (*,*) ' GRID/ z outside bounds of zgrid '
1663           write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz)
1664           z = zgrid(nz)
1665           write(*,*) 'WARNING: error in grid_r8 (jthermcalc.F)'
1666           write(*,*) 'Please check values of z and zgrid above'
1667        endif
1668        if ( nz .lt. 2 ) then
1669           write (*,*) ' GRID/ zgrid needs 2 points at least ! '
1670           stop ' Serious error in GRID.F '
1671        endif
1672        if ( zgrid(1) .ge. zgrid(nz) ) then
1673           write (*,*) ' GRID/ zgrid must increase with index'
1674           stop ' Serious error in GRID.F '
1675        endif
1676
1677        nz1 = 1
1678        nznew = nz/2
1679        if ( z .gt. zgrid(nznew) ) then
1680           nz1 = nznew
1681           nznew = nz
1682        endif
1683        do i=nz1+1,nznew
1684           if ( z. eq. zgrid(i) ) then
1685              grid_R8=i
1686              return
1687              elseif ( z .le. zgrid(i) ) then
1688              grid_R8 = i-1
1689              return
1690           endif
1691        enddo
1692        grid_R8 = nz
1693        return
1694
1695
1696
1697        end
1698
1699
1700
1701!c***************************************************
1702!c***************************************************
1703
1704      subroutine flujo(date)
1705
1706
1707!c     fgg           nov 2002     first version
1708!c***************************************************
1709
1710      use comsaison_h, only: dist_sol
1711      use param_v4_h, only: ninter,
1712     .                      fluxtop, ct1, ct2, p1, p2
1713      implicit none
1714
1715
1716!     common variables and constants
1717      include "callkeys.h"
1718
1719
1720!     Arguments
1721
1722      real date
1723
1724
1725!     Local variable and constants
1726
1727      integer i
1728      integer inter
1729      real    nada
1730
1731!c*************************************************
1732
1733      if(date.lt.1985.) date=1985.
1734      if(date.gt.2001.) date=2001.
1735     
1736      do i=1,ninter
1737         fluxtop(i)=1.
1738         !Variation of solar flux with 11 years solar cycle
1739         !For more details, see Gonzalez-Galindo et al. 2005
1740         !To be improved in next versions
1741        if(i.le.24.and.solvarmod.eq.0) then
1742            fluxtop(i)=(((ct1(i)+p1(i)*date)/2.)                 
1743     $           *sin(2.*3.1416/11.*(date-1985.-3.1416))         
1744     $           +(ct2(i)+p2(i)*date)+1.)*fluxtop(i)
1745         end if
1746         fluxtop(i)=fluxtop(i)*(1.52/dist_sol)**2
1747      end do
1748     
1749      return
1750      end
Note: See TracBrowser for help on using the repository browser.