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

Last change on this file since 643 was 635, checked in by emillour, 13 years ago

Mars GCM: Update of the chemistry package, including:

  • 93 reactions are accounted for (instead of 22); tracking 28 species (instead of 11)
  • computation of photoabsorption using raytracing
  • improved time stepping in the photochemistry
  • updated parameters (cross-sections); with this new version input files

are in 'EUV/param_v5' of "datafile" directory.

  • transition between lower and upper atmosphere chemistry set to 0.1 Pa (calchim.F90)
  • Lots of code clean-up: removed obsolete files column.F, param_v3.h, flujo.F, phdisrate.F, ch.F, interpfast.F, paramfoto.F, getch.F Converted chemtermos.F -> chemthermos.F90 and euvheat.F -> euvheat.F90. Added paramfoto_compact.F , param_v4.h and iono.h
  • Upadted surfacearea.F
  • Cleaned initracer.F and callkeys.h (removed use of obsolete "nqchem" and "oldnames" case when initializing tracers).
  • Minor correction in "callsedim": compute "rdust" and/or "rice" only when it makes sense.

FGG+FL+EM
File size: 84.4 KB
Line 
1c**********************************************************************
2
3      subroutine jthermcalc(ig,chemthermod,rm,nesptherm,tx,iz,zenit)
4
5
6c     feb 2002        fgg           first version
7c     nov 2002        fgg           second version
8c
9c modified from paramhr.F
10c MAC July 2003
11c**********************************************************************
12
13      implicit none
14
15c     common variables and constants
16      include "dimensions.h"
17      include "dimphys.h"
18      include 'param.h'
19      include 'param_v4.h'
20
21c     input and output variables
22
23      integer    ig
24      integer    chemthermod
25      integer    nesptherm                      !Number of species considered
26      real       rm(nlayermx,nesptherm)         !Densities (cm-3)
27      real       tx(nlayermx)                   !temperature
28      real       zenit                          !SZA
29      real       iz(nlayermx)                   !Local altitude
30
31
32c    local parameters and variables
33
34      real       aux1(nlayermx)                 !auxiliar variable
35      real       aux2(nlayermx)                 !auxiliar variable
36      real       co2colx(nlayermx)              !column density of CO2 (cm^-2)
37      real       o2colx(nlayermx)               !column density of O2(cm^-2)
38      real       o3pcolx(nlayermx)              !column density of O(3P)(cm^-2)
39      real       h2colx(nlayermx)               !H2 column density (cm-2)
40      real       h2ocolx(nlayermx)              !H2O column density (cm-2)
41      real       h2o2colx(nlayermx)             !column density of H2O2(cm^-2)
42      real       o3colx(nlayermx)               !O3 column density (cm-2)
43      real       n2colx(nlayermx)               !N2 column density (cm-2)
44      real       ncolx(nlayermx)                !N column density (cm-2)
45      real       nocolx(nlayermx)               !NO column density (cm-2)
46      real       cocolx(nlayermx)               !CO column density (cm-2)
47      real       hcolx(nlayermx)                !H column density (cm-2)
48      real       no2colx(nlayermx)              !NO2 column density (cm-2)
49      real       t2(nlayermx)
50      real       coltemp(nlayermx)
51      real       sigma(ninter,nlayermx)
52      real       alfa(ninter,nlayermx)
53     
54      integer    i,j,k,indexint                 !indexes
55      character  dn
56
57
58
59
60c     variables used in interpolation
61
62      real aux3(nz2), aux4(nz2)
63                                              !auxiliar variables
64      real       limdown                      !limits for interpolation
65      real       limup                        !        ""
66
67      !!!ATTENTION. Here i_co2 has to have the same value than in chemthermos.F90
68      !!!If the value is changed there, if has to be changed also here !!!!
69      integer,parameter :: i_co2=1
70
71c*************************PROGRAM STARTS*******************************
72     
73      if(zenit.gt.140.) then
74         dn='n'
75         else
76         dn='d'
77      end if
78      if(dn.eq.'n') then
79        return
80      endif
81     
82      !Initializing the photoabsorption coefficients
83      jfotsout(:,:,:)=0.
84
85      !Auxiliar temperature to take into account the temperature dependence
86      !of CO2 cross section
87      do i=1,nlayermx
88         t2(i)=tx(i)
89         if(t2(i).lt.195.0) t2(i)=195.0
90         if(t2(i).gt.295.0) t2(i)=295.0
91      end do
92
93      !Calculation of column amounts
94      call column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
95     $     co2colx,o2colx,o3pcolx,h2colx,h2ocolx,
96     $     h2o2colx,o3colx,n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
97
98      !Auxiliar column to include the temperature dependence
99      !of CO2 cross section
100      coltemp(nlayermx)=co2colx(nlayermx)*abs(t2(nlayermx)-t0(nlayermx))
101      do i=nlayermx-1,1,-1
102        coltemp(i)=!coltemp(i+1)+     PQ SE ELIMINA? REVISAR
103     $         ( rm(i,i_co2) + rm(i+1,i_co2) ) * 0.5
104     $         * 1e5 * (iz(i+1)-iz(i)) * abs(t2(i)-t0(i))
105      end do
106     
107      !Calculation of CO2 cross section at temperature t0(i)
108      do i=1,nlayermx
109         do indexint=24,32
110           sigma(indexint,i)=co2crsc195(indexint-23)
111           alfa(indexint,i)=((co2crsc295(indexint-23)
112     $          /sigma(indexint,i))-1.)/(295.-t0(i))
113        end do
114      end do
115
116! Interpolation to the tabulated photoabsorption coefficients for each species
117! in each spectral interval
118
119!AQUI SE PODRIAN AGRUPAR CALCULOS PARA AHORRAR TIEMPO DE CPU
120!P.EJ. LOS CALCULOS DE AUX2 y AUX4 NO ES NECESARIO REPETIRLOS
121!PARA CADA ESPECIE EN UN INTERVALO.
122!REVISAR
123
124!TAMBIEN SE PODRIA PONER, EN LUGAR DE CONDICIONES SOBRE CHEMTHERMOD
125!PARA VER QUE ESPECIES SE CONSIDERAN, PONER CONDICION SOBRE LA EXISTENCIA DE
126!CADA ESPECIE EN TRACEUR.DEF. ASI SE CALCULARIA LA FOTOABSORCION DE TODAS LAS
127!ESPECIES INCLUIDAS, AUNQUE LUEGO NO SE TENGA EN CUENTA EN LA QUIMICA (P.EJ.,
128!SI HAY NO2 PERO NO NO, NO SE CALCULARIA QUIMICA DEL NITROGENO PERO SE PODRIA
129!TENER EN CUENTA PARA EL CALENTAMIENTO
130!ESTUDIAR EN EL FUTURO
131
132!OTRA POSIBILIDAD ES SUSTITUIR LA ESCRITURA EN DURO DE LOS INDICES
133!DE LAS ESPECIES EN JABSIFOTS, CRSCABSI, ETC. POR INDICES DEL TIPO I_*
134!ESTUDIAR EN EL FUTURO
135c************************************************
136c     co2,0.1,6.0
137c************************************************
138
139      indexint=1
140      limdown=1e-20
141      limup=1e26
142      do i=1,nlayermx
143         aux1(i)=0.
144         aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint) +
145     $        o2colx(i)*crscabsi2(2,indexint) +
146     $        o3pcolx(i)*crscabsi2(3,indexint) +
147     $        h2colx(i)*crscabsi2(5,indexint) +
148     $        ncolx(i)*crscabsi2(9,indexint)
149      end do
150      do i=1,nz2
151         aux3(i) = jabsifotsintpar(indexint,1,nz2-i+1)
152         aux4(i) = c1_16(nz2-i+1,indexint)
153      enddo
154
155      call interpfast
156     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
157      do i=1,nlayermx
158         jfotsout(indexint,1,i) = aux1(nlayermx-i+1)
159      enddo
160
161c************************************************
162c     o2,0.1,6.0
163c************************************************
164
165      indexint=1
166      limdown=1e-20
167      limup=1e26
168      do i=1,nlayermx
169         aux1(i)=0.
170         aux2(nlayermx-i+1) =  co2colx(i)*crscabsi2(1,indexint) +
171     $        o2colx(i)*crscabsi2(2,indexint) +
172     $        o3pcolx(i)*crscabsi2(3,indexint) +
173     $        h2colx(i)*crscabsi2(5,indexint) +
174     $        ncolx(i)*crscabsi2(9,indexint)
175      end do
176      do i=1,nz2
177         aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
178         aux4(i) = c1_16(nz2-i+1,indexint)
179      enddo
180
181      call interpfast
182     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
183      do i=1,nlayermx
184         jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
185      enddo
186
187c**************************************************         
188c     o3p,0.1,6.0
189c**************************************************
190
191      indexint=1
192      limdown=1e-20
193      limup=1e26
194      do i=1,nlayermx
195         aux2(nlayermx-i+1) =  co2colx(i)*crscabsi2(1,indexint) +
196     $        o2colx(i)*crscabsi2(2,indexint) +
197     $        o3pcolx(i)*crscabsi2(3,indexint) +
198     $        h2colx(i)*crscabsi2(5,indexint) +
199     $        ncolx(i)*crscabsi2(9,indexint)
200      end do
201      do i=1,nz2
202         aux3(i) = jabsifotsintpar(indexint,3,nz2-i+1)
203         aux4(i) = c1_16(nz2-i+1,indexint)
204      enddo
205      call interpfast
206     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
207      do i=1,nlayermx
208         jfotsout(indexint,3,i) = aux1(nlayermx-i+1)
209      enddo
210
211     
212c**************************************************
213c     h2,0.1,6.0
214c**************************************************
215
216      indexint=1
217      limdown=1e-20
218      limup=1e26
219      do i=1,nlayermx
220         aux2(nlayermx-i+1) =  co2colx(i)*crscabsi2(1,indexint) +
221     $        o2colx(i)*crscabsi2(2,indexint) +
222     $        o3pcolx(i)*crscabsi2(3,indexint) +
223     $        h2colx(i)*crscabsi2(5,indexint) +
224     $        ncolx(i)*crscabsi2(9,indexint)
225      end do
226      do i=1,nz2
227         aux3(i) = jabsifotsintpar(indexint,5,nz2-i+1)
228         aux4(i) = c1_16(nz2-i+1,indexint)
229      enddo
230      call interpfast
231     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
232      do i=1,nlayermx
233         jfotsout(indexint,5,i) = aux1(nlayermx-i+1)
234      enddo
235
236      !Only if Nitrogen or ionospheric chemistry requested
237      if(chemthermod.ge.2) then
238
239c**************************************************
240c     n,0.1,6.0
241c**************************************************
242
243
244         indexint=1
245         limdown=1e-20
246         limup=1e26
247         do i=1,nlayermx
248            aux2(nlayermx-i+1) =  co2colx(i)*crscabsi2(1,indexint) +
249     $           o2colx(i)*crscabsi2(2,indexint) +
250     $           o3pcolx(i)*crscabsi2(3,indexint) +
251     $           h2colx(i)*crscabsi2(5,indexint) +
252     $           ncolx(i)*crscabsi2(9,indexint)
253         end do
254         do i=1,nz2
255            aux3(i) = jabsifotsintpar(indexint,9,nz2-i+1)
256            aux4(i) = c1_16(nz2-i+1,indexint)
257         enddo
258         call interpfast
259     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
260         do i=1,nlayermx
261            jfotsout(indexint,9,i) = aux1(nlayermx-i+1)
262         enddo
263
264      endif  !Of chemthermod >= 2
265     
266
267c**************************************************
268c     o2, 5-80.5nm
269c**************************************************
270
271      limdown=1e-20
272      limup=1e26
273      do indexint=2,15
274         do i=nlayermx-1,1,-1
275         end do
276         do i=1,nlayermx
277            aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
278     $           o2colx(i)*crscabsi2(2,indexint)+
279     $           o3pcolx(i)*crscabsi2(3,indexint)+
280     $           h2colx(i)*crscabsi2(5,indexint)+
281     $           n2colx(i)*crscabsi2(8,indexint)+
282     $           ncolx(i)*crscabsi2(9,indexint)+
283     $           nocolx(i)*crscabsi2(10,indexint)+
284     $           cocolx(i)*crscabsi2(11,indexint)+
285     $           hcolx(i)*crscabsi2(12,indexint)+
286     $           no2colx(i)*crscabsi2(13,indexint)
287         end do
288         do i=1,nz2
289            aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
290            aux4(i) = c1_16(nz2-i+1,indexint)
291         enddo
292          call interpfast
293     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
294          do i=1,nlayermx
295            jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
296          enddo
297      end do
298     
299
300c**************************************************
301c     o3p, 5-80.5nm
302c**************************************************
303
304      do indexint=2,15
305         do i=1,nlayermx
306            aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
307     $           o2colx(i)*crscabsi2(2,indexint)+
308     $           o3pcolx(i)*crscabsi2(3,indexint)+
309     $           h2colx(i)*crscabsi2(5,indexint)+
310     $           n2colx(i)*crscabsi2(8,indexint)+
311     $           ncolx(i)*crscabsi2(9,indexint)+
312     $           nocolx(i)*crscabsi2(10,indexint)+
313     $           cocolx(i)*crscabsi2(11,indexint)+
314     $           hcolx(i)*crscabsi2(12,indexint)+
315     $           no2colx(i)*crscabsi2(13,indexint)
316         end do
317         do i=1,nz2
318            aux3(i) = jabsifotsintpar(indexint,3,nz2-i+1)
319            aux4(i) = c1_16(nz2-i+1,indexint)
320         enddo
321         call interpfast
322     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
323         do i=1,nlayermx
324            jfotsout(indexint,3,i) = aux1(nlayermx-i+1)
325         enddo
326      end do
327
328c**************************************************
329c     co2, 5-80.5nm
330c**************************************************
331
332      do indexint=2,15
333         do i=1,nlayermx
334            aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
335     $           o2colx(i)*crscabsi2(2,indexint)+
336     $           o3pcolx(i)*crscabsi2(3,indexint)+
337     $           h2colx(i)*crscabsi2(5,indexint)+
338     $           n2colx(i)*crscabsi2(8,indexint)+
339     $           ncolx(i)*crscabsi2(9,indexint)+
340     $           nocolx(i)*crscabsi2(10,indexint)+
341     $           cocolx(i)*crscabsi2(11,indexint)+
342     $           hcolx(i)*crscabsi2(12,indexint)+
343     $           no2colx(i)*crscabsi2(13,indexint)
344         end do
345         do i=1,nz2
346            aux3(i) = jabsifotsintpar(indexint,1,nz2-i+1)
347            aux4(i) = c1_16(nz2-i+1,indexint)
348         enddo
349         call interpfast
350     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
351         do i=1,nlayermx
352            jfotsout(indexint,1,i) = aux1(nlayermx-i+1)
353         enddo
354      end do
355
356
357c**************************************************
358c     h2, 5-80.5nm
359c**************************************************
360
361      do indexint=2,15
362         do i=1,nlayermx
363            aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
364     $           o2colx(i)*crscabsi2(2,indexint)+
365     $           o3pcolx(i)*crscabsi2(3,indexint)+
366     $           h2colx(i)*crscabsi2(5,indexint)+
367     $           n2colx(i)*crscabsi2(8,indexint)+
368     $           ncolx(i)*crscabsi2(9,indexint)+
369     $           nocolx(i)*crscabsi2(10,indexint)+
370     $           cocolx(i)*crscabsi2(11,indexint)+
371     $           hcolx(i)*crscabsi2(12,indexint)+
372     $           no2colx(i)*crscabsi2(13,indexint)
373         end do
374         do i=1,nz2
375            aux3(i) = jabsifotsintpar(indexint,5,nz2-i+1)
376            aux4(i) = c1_16(nz2-i+1,indexint)
377         enddo
378         call interpfast
379     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
380         do i=1,nlayermx
381            jfotsout(indexint,5,i) = aux1(nlayermx-i+1)
382         enddo
383      end do
384
385
386c**************************************************
387c     n2, 5-80.5nm
388c**************************************************
389
390      do indexint=2,15
391         do i=1,nlayermx
392            aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
393     $           o2colx(i)*crscabsi2(2,indexint)+
394     $           o3pcolx(i)*crscabsi2(3,indexint)+
395     $           h2colx(i)*crscabsi2(5,indexint)+
396     $           n2colx(i)*crscabsi2(8,indexint)+
397     $           ncolx(i)*crscabsi2(9,indexint)+
398     $           nocolx(i)*crscabsi2(10,indexint)+
399     $           cocolx(i)*crscabsi2(11,indexint)+
400     $           hcolx(i)*crscabsi2(12,indexint)+
401     $           no2colx(i)*crscabsi2(13,indexint)
402         end do
403         do i=1,nz2
404            aux3(i) = jabsifotsintpar(indexint,8,nz2-i+1)
405            aux4(i) = c1_16(nz2-i+1,indexint)
406         enddo
407         call interpfast
408     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
409         do i=1,nlayermx
410            jfotsout(indexint,8,i) = aux1(nlayermx-i+1)
411         enddo
412      end do
413
414
415c**************************************************
416c     co, 5-80.5nm
417c**************************************************
418
419      do indexint=2,15
420         do i=1,nlayermx
421            aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
422     $           o2colx(i)*crscabsi2(2,indexint)+
423     $           o3pcolx(i)*crscabsi2(3,indexint)+
424     $           h2colx(i)*crscabsi2(5,indexint)+
425     $           n2colx(i)*crscabsi2(8,indexint)+
426     $           ncolx(i)*crscabsi2(9,indexint)+
427     $           nocolx(i)*crscabsi2(10,indexint)+
428     $           cocolx(i)*crscabsi2(11,indexint)+
429     $           hcolx(i)*crscabsi2(12,indexint)+
430     $           no2colx(i)*crscabsi2(13,indexint)
431         end do
432         do i=1,nz2
433            aux3(i) = jabsifotsintpar(indexint,11,nz2-i+1)
434            aux4(i) = c1_16(nz2-i+1,indexint)
435         enddo
436         call interpfast
437     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
438         do i=1,nlayermx
439            jfotsout(indexint,11,i) = aux1(nlayermx-i+1)
440         enddo
441      end do
442
443
444c**************************************************
445c     h, 5-80.5nm
446c**************************************************
447
448      do indexint=2,15
449         do i=1,nlayermx
450            aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
451     $           o2colx(i)*crscabsi2(2,indexint)+
452     $           o3pcolx(i)*crscabsi2(3,indexint)+
453     $           h2colx(i)*crscabsi2(5,indexint)+
454     $           n2colx(i)*crscabsi2(8,indexint)+
455     $           ncolx(i)*crscabsi2(9,indexint)+
456     $           nocolx(i)*crscabsi2(10,indexint)+
457     $           cocolx(i)*crscabsi2(11,indexint)+
458     $           hcolx(i)*crscabsi2(12,indexint)+
459     $           no2colx(i)*crscabsi2(13,indexint)
460         end do
461         do i=1,nz2
462            aux3(i) = jabsifotsintpar(indexint,12,nz2-i+1)
463            aux4(i) = c1_16(nz2-i+1,indexint)
464         enddo
465         call interpfast
466     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
467         do i=1,nlayermx
468            jfotsout(indexint,12,i) = aux1(nlayermx-i+1)
469         enddo
470      end do
471
472
473      !Only if Nitrogen or ionospheric chemistry requested
474      if(chemthermod.ge.2) then
475c**************************************************
476c     n, 5-80.5nm
477c**************************************************
478
479         do indexint=2,15
480            do i=1,nlayermx
481               aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
482     $              o2colx(i)*crscabsi2(2,indexint)+
483     $              o3pcolx(i)*crscabsi2(3,indexint)+
484     $              h2colx(i)*crscabsi2(5,indexint)+
485     $              n2colx(i)*crscabsi2(8,indexint)+
486     $              ncolx(i)*crscabsi2(9,indexint)+
487     $              nocolx(i)*crscabsi2(10,indexint)+
488     $              cocolx(i)*crscabsi2(11,indexint)+
489     $              hcolx(i)*crscabsi2(12,indexint)+
490     $              no2colx(i)*crscabsi2(13,indexint)
491            end do
492            do i=1,nz2
493               aux3(i) = jabsifotsintpar(indexint,9,nz2-i+1)
494               aux4(i) = c1_16(nz2-i+1,indexint)
495            enddo
496            call interpfast
497     $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
498            do i=1,nlayermx
499               jfotsout(indexint,9,i) = aux1(nlayermx-i+1)
500            enddo
501         end do
502
503
504c**************************************************
505c     no, 5-80.5nm
506c**************************************************
507
508         do indexint=2,15
509            do i=1,nlayermx
510               aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
511     $              o2colx(i)*crscabsi2(2,indexint)+
512     $              o3pcolx(i)*crscabsi2(3,indexint)+
513     $              h2colx(i)*crscabsi2(5,indexint)+
514     $              n2colx(i)*crscabsi2(8,indexint)+
515     $              ncolx(i)*crscabsi2(9,indexint)+
516     $              nocolx(i)*crscabsi2(10,indexint)+
517     $              cocolx(i)*crscabsi2(11,indexint)+
518     $              hcolx(i)*crscabsi2(12,indexint)+
519     $              no2colx(i)*crscabsi2(13,indexint)
520            end do
521            do i=1,nz2
522               aux3(i) = jabsifotsintpar(indexint,10,nz2-i+1)
523               aux4(i) = c1_16(nz2-i+1,indexint)
524            enddo
525            call interpfast
526     $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
527            do i=1,nlayermx
528               jfotsout(indexint,10,i) = aux1(nlayermx-i+1)
529            enddo
530         end do
531
532c**************************************************
533c     no2, 5-80.5nm
534c**************************************************
535
536         do indexint=2,15
537            do i=1,nlayermx
538               aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
539     $              o2colx(i)*crscabsi2(2,indexint)+
540     $              o3pcolx(i)*crscabsi2(3,indexint)+
541     $              h2colx(i)*crscabsi2(5,indexint)+
542     $              n2colx(i)*crscabsi2(8,indexint)+
543     $              ncolx(i)*crscabsi2(9,indexint)+
544     $              nocolx(i)*crscabsi2(10,indexint)+
545     $              cocolx(i)*crscabsi2(11,indexint)+
546     $              hcolx(i)*crscabsi2(12,indexint)+
547     $              no2colx(i)*crscabsi2(13,indexint)
548            end do
549            do i=1,nz2
550               aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
551               aux4(i) = c1_16(nz2-i+1,indexint)
552            enddo
553            call interpfast
554     $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
555            do i=1,nlayermx
556               jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
557            enddo
558         end do
559     
560         endif  !Of chemthermod >= 2
561
562         
563c**************************************************
564c     o2, 80.6-90.8nm
565c**************************************************
566
567      indexint=16
568      do i=1,nlayermx
569         aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
570     $        o2colx(i)*crscabsi2(2,indexint)+
571     $        o3pcolx(i)*crscabsi2(3,indexint)+
572     $        n2colx(i)*crscabsi2(8,indexint)+
573     $        ncolx(i)*crscabsi2(9,indexint)+
574     $        nocolx(i)*crscabsi2(10,indexint)+
575     $        cocolx(i)*crscabsi2(11,indexint)+
576     $        hcolx(i)*crscabsi2(12,indexint)+
577     $        no2colx(i)*crscabsi2(13,indexint)
578      end do
579      do i=1,nz2
580         aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
581         aux4(i) = c1_16(nz2-i+1,indexint)
582      enddo
583      call interpfast
584     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
585      do i=1,nlayermx
586         jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
587      enddo
588
589
590c**************************************************
591c     co2, 80.6-90.8nm
592c**************************************************
593
594      indexint=16
595      do i=1,nlayermx
596         aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
597     $        o2colx(i)*crscabsi2(2,indexint)+
598     $        o3pcolx(i)*crscabsi2(3,indexint)+
599     $        n2colx(i)*crscabsi2(8,indexint)+
600     $        ncolx(i)*crscabsi2(9,indexint)+
601     $        nocolx(i)*crscabsi2(10,indexint)+
602     $        cocolx(i)*crscabsi2(11,indexint)+
603     $        hcolx(i)*crscabsi2(12,indexint)+
604     $        no2colx(i)*crscabsi2(13,indexint)
605      end do
606      do i=1,nz2
607         aux3(i) = jabsifotsintpar(indexint,1,nz2-i+1)
608         aux4(i) = c1_16(nz2-i+1,indexint)
609      enddo
610      call interpfast
611     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
612      do i=1,nlayermx
613         jfotsout(indexint,1,i) = aux1(nlayermx-i+1)
614      enddo
615
616
617c**************************************************
618c     o3p, 80.6-90.8nm
619c**************************************************
620
621      indexint=16
622      do i=1,nlayermx
623         aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
624     $        o2colx(i)*crscabsi2(2,indexint)+
625     $        o3pcolx(i)*crscabsi2(3,indexint)+
626     $        n2colx(i)*crscabsi2(8,indexint)+
627     $        ncolx(i)*crscabsi2(9,indexint)+
628     $        nocolx(i)*crscabsi2(10,indexint)+
629     $        cocolx(i)*crscabsi2(11,indexint)+
630     $        hcolx(i)*crscabsi2(12,indexint)+
631     $        no2colx(i)*crscabsi2(13,indexint)
632      end do
633      do i=1,nz2
634         aux3(i) = jabsifotsintpar(indexint,3,nz2-i+1)
635         aux4(i) = c1_16(nz2-i+1,indexint)
636      enddo
637      call interpfast
638     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
639      do i=1,nlayermx
640         jfotsout(indexint,3,i) = aux1(nlayermx-i+1)
641      enddo
642     
643
644c**************************************************
645c     n2, 80.6-90.8nm
646c**************************************************
647
648      indexint=16
649      do i=1,nlayermx
650         aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
651     $        o2colx(i)*crscabsi2(2,indexint)+
652     $        o3pcolx(i)*crscabsi2(3,indexint)+
653     $        n2colx(i)*crscabsi2(8,indexint)+
654     $        ncolx(i)*crscabsi2(9,indexint)+
655     $        nocolx(i)*crscabsi2(10,indexint)+
656     $        cocolx(i)*crscabsi2(11,indexint)+
657     $        hcolx(i)*crscabsi2(12,indexint)+
658     $        no2colx(i)*crscabsi2(13,indexint)
659      end do
660      do i=1,nz2
661         aux3(i) = jabsifotsintpar(indexint,8,nz2-i+1)
662         aux4(i) = c1_16(nz2-i+1,indexint)
663      enddo
664      call interpfast
665     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
666      do i=1,nlayermx
667         jfotsout(indexint,8,i) = aux1(nlayermx-i+1)
668      enddo
669
670
671c**************************************************
672c     co, 80.6-90.8nm
673c************************************************************
674
675      indexint=16
676      do i=1,nlayermx
677         aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
678     $        o2colx(i)*crscabsi2(2,indexint)+
679     $        o3pcolx(i)*crscabsi2(3,indexint)+
680     $        n2colx(i)*crscabsi2(8,indexint)+
681     $        ncolx(i)*crscabsi2(9,indexint)+
682     $        nocolx(i)*crscabsi2(10,indexint)+
683     $        cocolx(i)*crscabsi2(11,indexint)+
684     $        hcolx(i)*crscabsi2(12,indexint)+
685     $        no2colx(i)*crscabsi2(13,indexint)
686      end do
687      do i=1,nz2
688         aux3(i) = jabsifotsintpar(indexint,11,nz2-i+1)
689         aux4(i) = c1_16(nz2-i+1,indexint)
690      enddo
691      call interpfast
692     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
693      do i=1,nlayermx
694         jfotsout(indexint,11,i) = aux1(nlayermx-i+1)
695      enddo
696
697
698c**************************************************
699c     h, 80.6-90.8nm
700c**************************************************
701
702      indexint=16
703      do i=1,nlayermx
704         aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
705     $        o2colx(i)*crscabsi2(2,indexint)+
706     $        o3pcolx(i)*crscabsi2(3,indexint)+
707     $        n2colx(i)*crscabsi2(8,indexint)+
708     $        ncolx(i)*crscabsi2(9,indexint)+
709     $        nocolx(i)*crscabsi2(10,indexint)+
710     $        cocolx(i)*crscabsi2(11,indexint)+
711     $        hcolx(i)*crscabsi2(12,indexint)+
712     $        no2colx(i)*crscabsi2(13,indexint)
713      end do
714      do i=1,nz2
715         aux3(i) = jabsifotsintpar(indexint,12,nz2-i+1)
716         aux4(i) = c1_16(nz2-i+1,indexint)
717      enddo
718      call interpfast
719     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
720      do i=1,nlayermx
721         jfotsout(indexint,12,i) = aux1(nlayermx-i+1)
722      enddo
723
724
725      !Only if Nitrogen or ionospheric chemistry requested
726      if(chemthermod.ge.2) then
727
728c**************************************************
729c     n, 80.6-90.8nm
730c**************************************************
731
732         indexint=16
733         do i=1,nlayermx
734            aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
735     $           o2colx(i)*crscabsi2(2,indexint)+
736     $           o3pcolx(i)*crscabsi2(3,indexint)+
737     $           n2colx(i)*crscabsi2(8,indexint)+
738     $           ncolx(i)*crscabsi2(9,indexint)+
739     $           nocolx(i)*crscabsi2(10,indexint)+
740     $           cocolx(i)*crscabsi2(11,indexint)+
741     $           hcolx(i)*crscabsi2(12,indexint)+
742     $           no2colx(i)*crscabsi2(13,indexint)
743         end do
744         do i=1,nz2
745            aux3(i) = jabsifotsintpar(indexint,9,nz2-i+1)
746            aux4(i) = c1_16(nz2-i+1,indexint)
747         enddo
748         call interpfast
749     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
750         do i=1,nlayermx
751            jfotsout(indexint,9,i) = aux1(nlayermx-i+1)
752         enddo
753
754
755c**************************************************
756c     no, 80.6-90.8nm
757c**************************************************
758
759         indexint=16
760         do i=1,nlayermx
761            aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
762     $           o2colx(i)*crscabsi2(2,indexint)+
763     $           o3pcolx(i)*crscabsi2(3,indexint)+
764     $           n2colx(i)*crscabsi2(8,indexint)+
765     $           ncolx(i)*crscabsi2(9,indexint)+
766     $           nocolx(i)*crscabsi2(10,indexint)+
767     $           cocolx(i)*crscabsi2(11,indexint)+
768     $           hcolx(i)*crscabsi2(12,indexint)+
769     $           no2colx(i)*crscabsi2(13,indexint)
770         end do
771         do i=1,nz2
772            aux3(i) = jabsifotsintpar(indexint,10,nz2-i+1)
773            aux4(i) = c1_16(nz2-i+1,indexint)
774         enddo
775         call interpfast
776     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
777         do i=1,nlayermx
778            jfotsout(indexint,10,i) = aux1(nlayermx-i+1)
779         enddo
780
781
782c***********************************************************
783c     no2, 80.6-90.8nm
784c**************************************************
785
786         indexint=16
787         do i=1,nlayermx
788            aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
789     $           o2colx(i)*crscabsi2(2,indexint)+
790     $           o3pcolx(i)*crscabsi2(3,indexint)+
791     $           n2colx(i)*crscabsi2(8,indexint)+
792     $           ncolx(i)*crscabsi2(9,indexint)+
793     $           nocolx(i)*crscabsi2(10,indexint)+
794     $           cocolx(i)*crscabsi2(11,indexint)+
795     $           hcolx(i)*crscabsi2(12,indexint)+
796     $           no2colx(i)*crscabsi2(13,indexint)
797         end do
798         do i=1,nz2
799            aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
800            aux4(i) = c1_16(nz2-i+1,indexint)
801         enddo
802         call interpfast
803     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
804         do i=1,nlayermx
805            jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
806         enddo
807     
808         endif !Of chemthermod >= 2
809
810c**************************************************
811c      co2, 90.9-119.5nm
812c**************************************************
813
814      limdown=1e-20
815      limup=1e26
816      do indexint=17,24
817         do i=1,nlayermx
818            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + n2colx(i) +
819     $           nocolx(i) + cocolx(i) + no2colx(i)
820         end do
821         do i=1,nz2
822            aux3(i) = jabsifotsintpar(indexint,1,nz2-i+1)
823            aux4(i) = c17_24(nz2-i+1)
824         enddo
825         call interpfast
826     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
827         do i=1,nlayermx
828            jfotsout(indexint,1,i) = aux1(nlayermx-i+1)
829            if(indexint.eq.24) then
830          if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
831               jfotsout(indexint,1,i) = aux1(nlayermx-i+1)
832     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))*
833     $           (1+alfa(indexint,i)*(t2(i)-t0(i)))
834               else
835                  jfotsout(indexint,1,i)=0.
836               end if
837            end if
838         enddo
839      end do
840
841
842c**************************************************
843c     o2, 90.9-119.5nm
844c**************************************************
845     
846      do indexint=17,24
847         do i=1,nlayermx
848            aux2(nlayermx-i+1) = o2colx(i) + co2colx(i) + n2colx(i) +
849     $           nocolx(i) + cocolx(i) + no2colx(i)
850         end do
851         do i=1,nz2
852            aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
853            aux4(i) = c17_24(nz2-i+1)
854         enddo
855         call interpfast
856     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
857         do i=1,nlayermx
858            jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
859            if(indexint.eq.24) then
860          if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
861               jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
862     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
863               else
864                  jfotsout(indexint,2,i)=0.
865               end if
866            end if
867         enddo
868      end do
869
870
871c**************************************************
872c     n2, 90.9-119.5nm
873c**************************************************
874     
875      do indexint=17,24
876         do i=1,nlayermx
877            aux2(nlayermx-i+1) = o2colx(i) + co2colx(i) + n2colx(i) +
878     $           nocolx(i) + cocolx(i) + no2colx(i)
879         end do
880         do i=1,nz2
881            aux3(i) = jabsifotsintpar(indexint,8,nz2-i+1)
882            aux4(i) = c17_24(nz2-i+1)
883         enddo
884         call interpfast
885     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
886         do i=1,nlayermx
887            jfotsout(indexint,8,i) = aux1(nlayermx-i+1)
888            if(indexint.eq.24) then
889          if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
890               jfotsout(indexint,8,i) = aux1(nlayermx-i+1)
891     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
892               else
893                  jfotsout(indexint,8,i)=0.
894               end if
895            end if
896         enddo
897      end do
898
899
900c**************************************************
901c     co, 90.9-119.5nm
902c**************************************************
903     
904      do indexint=17,24
905         do i=1,nlayermx
906            aux2(nlayermx-i+1) = o2colx(i) + co2colx(i) + n2colx(i) +
907     $           nocolx(i) + cocolx(i) + no2colx(i)
908         end do
909         do i=1,nz2
910            aux3(i) = jabsifotsintpar(indexint,11,nz2-i+1)
911            aux4(i) = c17_24(nz2-i+1)
912         enddo
913         call interpfast
914     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
915         do i=1,nlayermx
916            jfotsout(indexint,11,i) = aux1(nlayermx-i+1)
917            if(indexint.eq.24) then
918          if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
919               jfotsout(indexint,11,i) = aux1(nlayermx-i+1)
920     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
921               else
922                  jfotsout(indexint,11,i)=0.
923               end if
924            end if
925         enddo
926      end do
927
928
929      !Only if Nitrogen or ionospheric chemistry requested
930      if(chemthermod.ge.2) then
931
932c**************************************************
933c     no, 90.9-119.5nm
934c**************************************************
935     
936         do indexint=17,24
937            do i=1,nlayermx
938               aux2(nlayermx-i+1) = o2colx(i) + co2colx(i) + n2colx(i) +
939     $              nocolx(i) + cocolx(i) + no2colx(i)
940            end do
941            do i=1,nz2
942               aux3(i) = jabsifotsintpar(indexint,10,nz2-i+1)
943               aux4(i) = c17_24(nz2-i+1)
944            enddo
945            call interpfast
946     $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
947            do i=1,nlayermx
948               jfotsout(indexint,10,i) = aux1(nlayermx-i+1)
949               if(indexint.eq.24) then
950                  if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i)
951     $                 .lt.60.) then
952                     jfotsout(indexint,10,i) = aux1(nlayermx-i+1)
953     $             *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
954                  else
955                     jfotsout(indexint,10,i)=0.
956                  end if
957               end if
958            enddo
959         end do
960
961
962c**************************************************
963c     no2, 90.9-119.5nm
964c**************************************************
965     
966         do indexint=17,24
967            do i=1,nlayermx
968               aux2(nlayermx-i+1) = o2colx(i) + co2colx(i) + n2colx(i) +
969     $              nocolx(i) + cocolx(i) + no2colx(i)
970            end do
971            do i=1,nz2
972               aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
973               aux4(i) = c17_24(nz2-i+1)
974            enddo
975            call interpfast
976     $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
977            do i=1,nlayermx
978               jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
979               if(indexint.eq.24) then
980                  if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i)
981     $                 .lt.60.) then
982                     jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
983     $              *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
984                  else
985                     jfotsout(indexint,13,i)=0.
986                  end if
987               end if
988            enddo
989         end do
990
991         endif !Of chemthermod >= 2
992
993
994c**************************************************
995c      co2, 119.6-167.0nm
996c**************************************************
997
998      do indexint=25,29
999         do i=1,nlayermx
1000            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
1001     $           h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
1002         end do
1003         do i=1,nz2
1004            aux3(i) = jabsifotsintpar(indexint,1,nz2-i+1)
1005            aux4(i) = c25_29(nz2-i+1)
1006         enddo
1007         call interpfast
1008     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1009         do i=1,nlayermx
1010          if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
1011               jfotsout(indexint,1,i) = aux1(nlayermx-i+1)
1012     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
1013     $           *(1+alfa(indexint,i)*(t2(i)-t0(i)))               
1014            else
1015               jfotsout(indexint,1,i) = 0.
1016            end if
1017         enddo
1018      end do
1019
1020
1021c**************************************************
1022c      o2, 119.6-167.0nm
1023c**************************************************
1024
1025      do indexint=25,29
1026         do i=1,nlayermx
1027            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
1028     $           h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
1029         end do
1030         do i=1,nz2
1031            aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
1032            aux4(i) = c25_29(nz2-i+1)
1033         enddo
1034         call interpfast
1035     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1036         do i=1,nlayermx
1037           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
1038            jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
1039     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
1040            else
1041               jfotsout(indexint,2,i)=0.
1042            end if
1043         enddo
1044      end do
1045
1046
1047c**************************************************
1048c      h2o, 119.6-167.0nm
1049c**************************************************
1050
1051      do indexint=25,29
1052         do i=1,nlayermx
1053            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
1054     $           h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
1055         end do
1056         do i=1,nz2
1057            aux3(i) = jabsifotsintpar(indexint,4,nz2-i+1)
1058            aux4(i) = c25_29(nz2-i+1)
1059         enddo
1060         call interpfast
1061     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1062         do i=1,nlayermx
1063           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
1064            jfotsout(indexint,4,i) = aux1(nlayermx-i+1)
1065     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
1066            else
1067               jfotsout(indexint,4,i) = 0.
1068            end if
1069         enddo
1070      end do
1071
1072
1073c**************************************************
1074c      h2o2,119.6-167.0nm
1075c**************************************************
1076
1077      do indexint=25,29
1078         do i=1,nlayermx
1079            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
1080     $           h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
1081         end do
1082         do i=1,nz2
1083            aux3(i) = jabsifotsintpar(indexint,6,nz2-i+1)
1084            aux4(i) = c25_29(nz2-i+1)
1085         enddo
1086         call interpfast
1087     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1088         do i=1,nlayermx
1089           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
1090            jfotsout(indexint,6,i) = aux1(nlayermx-i+1)
1091     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
1092            else
1093               jfotsout(indexint,6,i) = 0.
1094            end if
1095         enddo
1096      end do
1097
1098
1099c**************************************************
1100c      co, 119.6-167.0nm
1101c**************************************************
1102
1103      do indexint=25,29
1104         do i=1,nlayermx
1105            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
1106     $           h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
1107         end do
1108         do i=1,nz2
1109            aux3(i) = jabsifotsintpar(indexint,11,nz2-i+1)
1110            aux4(i) = c25_29(nz2-i+1)
1111         enddo
1112         call interpfast
1113     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1114         do i=1,nlayermx
1115           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
1116            jfotsout(indexint,11,i) = aux1(nlayermx-i+1)
1117     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
1118            else
1119               jfotsout(indexint,11,i) = 0.
1120            end if
1121         enddo
1122      end do
1123
1124
1125      !Only if Nitrogen or ionospheric chemistry requested
1126      if(chemthermod.ge.2) then
1127
1128c**************************************************
1129c      no, 119.6-167.0nm
1130c**************************************************
1131
1132         do indexint=25,29
1133            do i=1,nlayermx
1134               aux2(nlayermx-i+1)=co2colx(i) + o2colx(i) + h2ocolx(i) +
1135     $              h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
1136            end do
1137            do i=1,nz2
1138               aux3(i) = jabsifotsintpar(indexint,10,nz2-i+1)
1139               aux4(i) = c25_29(nz2-i+1)
1140            enddo
1141            call interpfast
1142     $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1143            do i=1,nlayermx
1144               if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i)
1145     $              .lt.60.) then
1146                  jfotsout(indexint,10,i) = aux1(nlayermx-i+1)
1147     $             *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
1148               else
1149                  jfotsout(indexint,10,i) = 0.
1150               end if
1151            enddo
1152         end do
1153
1154
1155c**************************************************
1156c      no2, 119.6-167.0nm
1157c**************************************************
1158
1159         do indexint=25,29
1160            do i=1,nlayermx
1161               aux2(nlayermx-i+1)=co2colx(i) + o2colx(i) + h2ocolx(i) +
1162     $              h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
1163            end do
1164            do i=1,nz2
1165               aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
1166               aux4(i) = c25_29(nz2-i+1)
1167            enddo
1168            call interpfast
1169     $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1170            do i=1,nlayermx
1171               if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i)
1172     $              .lt.60.) then
1173                  jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
1174     $             *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
1175               else
1176                  jfotsout(indexint,13,i) = 0.
1177               end if
1178            enddo
1179         end do
1180
1181         endif !Of chemthermod >= 2
1182
1183
1184c**************************************************
1185c      co2, 167.0-202.5nm
1186c**************************************************
1187
1188      do indexint=30,31
1189         do i=1,nlayermx
1190            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
1191     $           h2o2colx(i) + nocolx(i) + no2colx(i)
1192         end do
1193         do i=1,nz2
1194            aux3(i) = jabsifotsintpar(indexint,1,nz2-i+1)
1195            aux4(i) = c30_31(nz2-i+1)
1196         enddo
1197         call interpfast
1198     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1199         do i=1,nlayermx
1200          if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
1201               jfotsout(indexint,1,i) = aux1(nlayermx-i+1)
1202     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
1203     $           *(1+alfa(indexint,i)*(t2(i)-t0(i)))               
1204            else
1205               jfotsout(indexint,1,i) = 0.
1206            end if
1207         enddo
1208      end do
1209
1210
1211c**************************************************
1212c      o2, 167.0-202.5nm
1213c**************************************************
1214
1215      do indexint=30,31
1216         do i=1,nlayermx
1217            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
1218     $           h2o2colx(i) + nocolx(i) + no2colx(i)
1219         end do
1220         do i=1,nz2
1221            aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
1222            aux4(i) = c30_31(nz2-i+1)
1223         enddo
1224         call interpfast
1225     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1226         do i=1,nlayermx
1227           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
1228            jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
1229     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
1230            else
1231               jfotsout(indexint,2,i)=0.
1232            end if
1233         enddo
1234      end do
1235
1236
1237c**************************************************
1238c      h2o, 167.0-202.5nm
1239c**************************************************
1240
1241      do indexint=30,31
1242         do i=1,nlayermx
1243            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
1244     $           h2o2colx(i) + nocolx(i) + no2colx(i)
1245         end do
1246         do i=1,nz2
1247            aux3(i) = jabsifotsintpar(indexint,4,nz2-i+1)
1248            aux4(i) = c30_31(nz2-i+1)
1249         enddo
1250         call interpfast
1251     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1252         do i=1,nlayermx
1253           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
1254            jfotsout(indexint,4,i) = aux1(nlayermx-i+1)
1255     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
1256            else
1257               jfotsout(indexint,4,i) = 0.
1258            end if
1259         enddo
1260      end do
1261
1262
1263c**************************************************
1264c      h2o2, 167.0-202.5nm
1265c**************************************************
1266
1267      do indexint=30,31
1268         do i=1,nlayermx
1269            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
1270     $           h2o2colx(i) + nocolx(i) + no2colx(i)
1271         end do
1272         do i=1,nz2
1273            aux3(i) = jabsifotsintpar(indexint,6,nz2-i+1)
1274            aux4(i) = c30_31(nz2-i+1)
1275         enddo
1276         call interpfast
1277     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1278         do i=1,nlayermx
1279           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
1280            jfotsout(indexint,6,i) = aux1(nlayermx-i+1)
1281     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
1282            else
1283               jfotsout(indexint,6,i) = 0.
1284            end if
1285         enddo
1286      end do
1287
1288
1289      !Only if Nitrogen or ionospheric chemistry requested
1290      if(chemthermod.ge.2) then
1291     
1292c**************************************************
1293c      no, 167.0-202.5nm
1294c**************************************************
1295
1296         do indexint=30,31
1297            do i=1,nlayermx
1298               aux2(nlayermx-i+1)=co2colx(i) + o2colx(i) + h2ocolx(i) +
1299     $              h2o2colx(i) + nocolx(i) + no2colx(i)
1300            end do
1301            do i=1,nz2
1302               aux3(i) = jabsifotsintpar(indexint,10,nz2-i+1)
1303               aux4(i) = c30_31(nz2-i+1)
1304            enddo
1305            call interpfast
1306     $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1307            do i=1,nlayermx
1308               if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i)
1309     $              .lt.60.) then
1310                  jfotsout(indexint,10,i) = aux1(nlayermx-i+1)
1311     $              *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
1312               else
1313                  jfotsout(indexint,10,i) = 0.
1314               end if
1315            enddo
1316         end do
1317
1318
1319c**************************************************
1320c      no2, 167.0-202.5nm
1321c**************************************************
1322
1323         do indexint=30,31
1324            do i=1,nlayermx
1325               aux2(nlayermx-i+1)=co2colx(i) + o2colx(i) + h2ocolx(i) +
1326     $              h2o2colx(i) + nocolx(i) + no2colx(i)
1327            end do
1328            do i=1,nz2
1329               aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
1330               aux4(i) = c30_31(nz2-i+1)
1331            enddo
1332            call interpfast
1333     $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1334            do i=1,nlayermx
1335               if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i)
1336     $              .lt.60.) then
1337                  jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
1338     $              *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
1339               else
1340                  jfotsout(indexint,13,i) = 0.
1341               end if
1342            enddo
1343         end do
1344         
1345      endif !Of chemthermod >= 2
1346
1347c**************************************************
1348c     co2, 202.6-210.0nm
1349c**************************************************
1350
1351      indexint=32
1352      do i=1,nlayermx
1353         aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2o2colx(i) +
1354     $        nocolx(i) + no2colx(i)
1355      end do
1356      do i=1,nz2
1357         aux3(i) = jabsifotsintpar(indexint,1,nz2-i+1)
1358         aux4(i) = c32(nz2-i+1)
1359      enddo
1360      call interpfast
1361     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1362      do i=1,nlayermx
1363         if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
1364            jfotsout(indexint,1,i) = aux1(nlayermx-i+1)
1365     $           *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
1366     $           *(1+alfa(indexint,i)*(t2(i)-t0(i)))
1367         else
1368            jfotsout(indexint,1,i)=0.
1369         end if
1370      enddo
1371
1372
1373c**************************************************
1374c     o2, 202.6-210.0nm
1375c**************************************************
1376
1377      indexint=32
1378      do i=1,nlayermx
1379         aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2o2colx(i) +
1380     $        nocolx(i) + no2colx(i)
1381      end do
1382      do i=1,nz2
1383         aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
1384         aux4(i) = c32(nz2-i+1)
1385      enddo
1386      call interpfast
1387     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1388      do i=1,nlayermx
1389         if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
1390            jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
1391     $           *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
1392         else
1393            jfotsout(indexint,2,i)=0.
1394         end if
1395      enddo
1396
1397
1398c**************************************************
1399c     h2o2, 202.6-210.0nm
1400c**************************************************
1401
1402      indexint=32
1403      do i=1,nlayermx
1404         aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2o2colx(i) +
1405     $        nocolx(i) + no2colx(i)
1406      end do
1407      do i=1,nz2
1408         aux3(i) = jabsifotsintpar(indexint,6,nz2-i+1)
1409         aux4(i) = c32(nz2-i+1)
1410      enddo
1411      call interpfast
1412     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1413      do i=1,nlayermx
1414         if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
1415            jfotsout(indexint,6,i) = aux1(nlayermx-i+1)
1416     $           *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
1417         else
1418            jfotsout(indexint,6,i)=0.
1419         end if
1420      enddo
1421
1422
1423      !Only if Nitrogen or ionospheric chemistry requested
1424      if(chemthermod.eq.2) then
1425
1426c**************************************************
1427c     no, 202.6-210.0nm
1428c**************************************************
1429
1430         indexint=32
1431         do i=1,nlayermx
1432            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2o2colx(i) +
1433     $           nocolx(i) + no2colx(i)
1434         end do
1435         do i=1,nz2
1436            aux3(i) = jabsifotsintpar(indexint,10,nz2-i+1)
1437            aux4(i) = c32(nz2-i+1)
1438         enddo
1439         call interpfast
1440     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1441         do i=1,nlayermx
1442            if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i)
1443     $           .lt.60.) then
1444               jfotsout(indexint,10,i) = aux1(nlayermx-i+1)
1445     $              *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
1446            else
1447               jfotsout(indexint,10,i)=0.
1448            end if
1449         enddo
1450
1451
1452c**************************************************
1453c     no2, 202.6-210.0nm
1454c**************************************************
1455
1456         indexint=32
1457         do i=1,nlayermx
1458            aux2(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2o2colx(i) +
1459     $           nocolx(i) + no2colx(i)
1460         end do
1461         do i=1,nz2
1462            aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
1463            aux4(i) = c32(nz2-i+1)
1464         enddo
1465         call interpfast
1466     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1467         do i=1,nlayermx
1468            if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i)
1469     $           .lt.60.) then
1470               jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
1471     $              *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
1472            else
1473               jfotsout(indexint,13,i)=0.
1474            end if
1475         enddo
1476         
1477      endif !Of chemthermod >= 2
1478
1479
1480c**************************************************
1481c     o2, 210.0-231.0nm
1482c**************************************************
1483
1484      indexint=33
1485      do i=1,nlayermx
1486         aux2(nlayermx-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
1487      end do
1488      do i=1,nz2
1489         aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
1490         aux4(i) = c33(nz2-i+1)
1491      enddo
1492      call interpfast
1493     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1494      do i=1,nlayermx
1495         jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
1496      enddo
1497
1498
1499c**************************************************
1500c     h2o2, 210.1-231.0nm
1501c**************************************************
1502
1503      indexint=33
1504      limdown=1.e-20
1505      limup=1e25
1506      do i=1,nlayermx
1507         aux2(nlayermx-i+1) = h2o2colx(i) + o2colx(i) + no2colx(i)
1508      end do
1509      do i=1,nz2
1510         aux3(i) = jabsifotsintpar(indexint,6,nz2-i+1)
1511         aux4(i) = c33(nz2-i+1)
1512      enddo
1513      call interpfast
1514     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1515      do i=1,nlayermx
1516         jfotsout(indexint,6,i) = aux1(nlayermx-i+1)
1517      enddo
1518
1519
1520      !Only if Nitrogen or ionospheric chemistry requested
1521      if(chemthermod.ge.2) then
1522
1523c**************************************************
1524c     no2, 210.1-231.0nm
1525c**************************************************
1526
1527         indexint=33
1528         limdown=1.e-20
1529         limup=1e25
1530         do i=1,nlayermx
1531            aux2(nlayermx-i+1) = h2o2colx(i) + o2colx(i) + no2colx(i)
1532         end do
1533         do i=1,nz2
1534            aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
1535            aux4(i) = c33(nz2-i+1)
1536         enddo
1537         call interpfast
1538     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1539         do i=1,nlayermx
1540            jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
1541         enddo
1542
1543      endif !Of chemthermod >= 2
1544
1545
1546c**************************************************
1547c     o2, 231.-240.nm
1548c**************************************************
1549
1550      indexint=34
1551      limdown=1.e-20
1552      limup=1e25
1553      do i=1,nlayermx
1554         aux2(nlayermx-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) +
1555     $        no2colx(i)
1556      end do
1557      do i=1,nz2
1558         aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
1559         aux4(i) = c34(nz2-i+1)
1560      enddo
1561      call interpfast
1562     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1563      do i=1,nlayermx
1564         jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
1565      enddo
1566
1567
1568c**************************************************
1569c     h2o2, 231.0-240.nm
1570c**************************************************
1571
1572      indexint=34
1573      limdown=1.e-20
1574      limup=1e25
1575      do i=1,nlayermx
1576         aux2(nlayermx-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) +
1577     $        no2colx(i)
1578      end do
1579      do i=1,nz2
1580         aux3(i) = jabsifotsintpar(indexint,6,nz2-i+1)
1581         aux4(i) = c34(nz2-i+1)
1582      enddo
1583      call interpfast
1584     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1585      do i=1,nlayermx
1586         jfotsout(indexint,6,i) = aux1(nlayermx-i+1)
1587      enddo
1588
1589
1590      !Only if Ozone, Nitrogen or ionospheric chem. requested
1591      if(chemthermod.ge.1) then
1592
1593c**************************************************
1594c     o3, 231.0-240.nm
1595c**************************************************
1596
1597         indexint=34
1598         limdown=1.e-20
1599         limup=1e25
1600         do i=1,nlayermx
1601            aux2(nlayermx-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) +
1602     $           no2colx(i)
1603         end do
1604         do i=1,nz2
1605            aux3(i) = jabsifotsintpar(indexint,7,nz2-i+1)
1606            aux4(i) = c34(nz2-i+1)
1607         enddo
1608         call interpfast
1609     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1610         do i=1,nlayermx
1611            jfotsout(indexint,7,i) = aux1(nlayermx-i+1)
1612         enddo
1613
1614      endif !Of chemthermod.ge.1
1615
1616
1617      !Only if nitrogen or ionospheric chemistry requested
1618      if(chemthermod.ge.2) then
1619
1620c**************************************************
1621c     no2, 231.0-240.nm
1622c**************************************************
1623
1624         indexint=34
1625         limdown=1.e-20
1626         limup=1e25
1627         do i=1,nlayermx
1628            aux2(nlayermx-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) +
1629     $           no2colx(i)
1630         end do
1631         do i=1,nz2
1632            aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
1633            aux4(i) = c34(nz2-i+1)
1634         enddo
1635         call interpfast
1636     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1637         do i=1,nlayermx
1638            jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
1639         enddo
1640
1641      endif !Of chemthermod >= 2
1642
1643
1644c**************************************************
1645c     h2o2, 240.0-337.7nm
1646c**************************************************
1647
1648      indexint=35
1649      limdown=1.e-20
1650      limup=1e25
1651      do i=1,nlayermx
1652         aux2(nlayermx-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
1653      end do
1654      do i=1,nz2
1655         aux3(i) = jabsifotsintpar(indexint,6,nz2-i+1)
1656         aux4(i) = c35(nz2-i+1)
1657      enddo
1658      call interpfast
1659     $     (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1660      do i=1,nlayermx
1661         jfotsout(indexint,6,i) = aux1(nlayermx-i+1)
1662      enddo
1663     
1664
1665      !Only if Ozone, Nitrogen or ionospheric chem. requested
1666      if(chemthermod.ge.1) then
1667
1668c**************************************************
1669c     o3, 240.0-337.7nm
1670c**************************************************
1671
1672         indexint=35
1673         limdown=1.e-20
1674         limup=1e25
1675         do i=1,nlayermx
1676            aux2(nlayermx-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
1677         end do
1678         do i=1,nz2
1679            aux3(i) = jabsifotsintpar(indexint,7,nz2-i+1)
1680            aux4(i) = c35(nz2-i+1)
1681         enddo
1682         call interpfast
1683     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1684         do i=1,nlayermx
1685            jfotsout(indexint,7,i) = aux1(nlayermx-i+1)
1686         enddo
1687
1688      endif !Of chemthermod.ge.1
1689
1690
1691      !Only if Nitrogen or ionospheric chemistry requested
1692      if(chemthermod.ge.2) then
1693
1694c**************************************************
1695c     no2, 240.0-337.7nm
1696c**************************************************
1697
1698         indexint=35
1699         limdown=1.e-20
1700         limup=1e25
1701         do i=1,nlayermx
1702            aux2(nlayermx-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
1703         end do
1704         do i=1,nz2
1705            aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
1706            aux4(i) = c35(nz2-i+1)
1707         enddo
1708         call interpfast
1709     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1710         do i=1,nlayermx
1711            jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
1712         enddo
1713
1714      endif !Of chemthermod.ge.2
1715
1716
1717      !Only if Ozone, Nitrogen or ionospheric chem. requested
1718      if(chemthermod.ge.1) then
1719
1720c**************************************************
1721c     o3, 337.7-800.0nm
1722c**************************************************
1723
1724         indexint=36
1725         limdown=1.e-20
1726         limup=1e25
1727         do i=1,nlayermx
1728            aux2(nlayermx-i+1) = o3colx(i) + no2colx(i)
1729         end do
1730         do i=1,nz2
1731            aux3(i) = jabsifotsintpar(indexint,7,nz2-i+1)
1732            aux4(i) = c36(nz2-i+1)
1733         enddo
1734         call interpfast
1735     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1736         do i=1,nlayermx
1737            jfotsout(indexint,7,i) = aux1(nlayermx-i+1)
1738         enddo
1739
1740      endif
1741
1742
1743      !Only if Nitrogen or ionospheric chemistry requested
1744      if(chemthermod.ge.2) then
1745
1746c**************************************************
1747c     no2, 337.7-800.0nm
1748c**************************************************
1749
1750         indexint=36
1751         limdown=1.e-20
1752         limup=1e25
1753         do i=1,nlayermx
1754            aux2(nlayermx-i+1) = o3colx(i) + no2colx(i)
1755         end do
1756         do i=1,nz2
1757            aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
1758            aux4(i) = c36(nz2-i+1)
1759         enddo
1760         call interpfast
1761     $        (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
1762         do i=1,nlayermx
1763            jfotsout(indexint,13,i) = aux1(nlayermx-i+1)
1764         enddo
1765
1766      endif !Of chemthermod.ge.2
1767      return
1768
1769      end
1770
1771
1772
1773c**********************************************************************
1774c**********************************************************************
1775
1776      subroutine column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
1777     $     co2colx,o2colx,o3pcolx,h2colx,h2ocolx,h2o2colx,o3colx,
1778     $     n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
1779
1780c     nov 2002        fgg           first version
1781
1782c**********************************************************************
1783
1784      implicit none
1785
1786
1787c     common variables and constants
1788      include "dimensions.h"
1789      include "dimphys.h"
1790      include "tracer.h"
1791      include 'param.h'
1792      include 'param_v4.h'
1793      include 'callkeys.h'
1794
1795
1796
1797c    local parameters and variables
1798
1799
1800
1801c     input and output variables
1802
1803      integer    ig
1804      integer    chemthermod
1805      integer    nesptherm                      !# of species undergoing chemistry, input
1806      real       rm(nlayermx,nesptherm)         !densities (cm-3), input
1807      real       tx(nlayermx)                   !temperature profile, input
1808      real       iz(nlayermx+1)                 !height profile, input
1809      real       zenit                          !SZA, input
1810      real       co2colx(nlayermx)              !column density of CO2 (cm^-2), output
1811      real       o2colx(nlayermx)               !column density of O2(cm^-2), output
1812      real       o3pcolx(nlayermx)              !column density of O(3P)(cm^-2), output
1813      real       h2colx(nlayermx)               !H2 column density (cm-2), output
1814      real       h2ocolx(nlayermx)              !H2O column density (cm-2), output
1815      real       h2o2colx(nlayermx)             !column density of H2O2(cm^-2), output
1816      real       o3colx(nlayermx)               !O3 column density (cm-2), output
1817      real       n2colx(nlayermx)               !N2 column density (cm-2), output
1818      real       ncolx(nlayermx)                !N column density (cm-2), output
1819      real       nocolx(nlayermx)               !NO column density (cm-2), output
1820      real       cocolx(nlayermx)               !CO column density (cm-2), output
1821      real       hcolx(nlayermx)                !H column density (cm-2), output
1822      real       no2colx(nlayermx)              !NO2 column density (cm-2), output
1823     
1824
1825c     local variables
1826
1827      real       xx
1828      real       grav(nlayermx)
1829      real       Hco2,Ho3p,Ho2,Hh2,Hh2o,Hh2o2
1830      real       Ho3,Hn2,Hn,Hno,Hco,Hh,Hno2
1831
1832      real       co2x(nlayermx)
1833      real       o2x(nlayermx)
1834      real       o3px(nlayermx)
1835      real       cox(nlayermx)
1836      real       hx(nlayermx)
1837      real       h2x(nlayermx)
1838      real       h2ox(nlayermx)
1839      real       h2o2x(nlayermx)
1840      real       o3x(nlayermx)
1841      real       n2x(nlayermx)
1842      real       nx(nlayermx)
1843      real       nox(nlayermx)
1844      real       no2x(nlayermx)
1845
1846      integer    i,j,k,icol,indexint          !indexes
1847
1848c     variables for optical path calculation
1849
1850      integer    nz3
1851!      parameter  (nz3=nz*2)
1852
1853      integer    jj
1854      real*8      esp(nlayermx*2)
1855      real*8      ilayesp(nlayermx*2)
1856      real*8      szalayesp(nlayermx*2)
1857      integer     nlayesp
1858      real*8      zmini
1859      real*8      depth
1860      real*8      espco2, espo2, espo3p, esph2, esph2o, esph2o2,espo3
1861      real*8      espn2,espn,espno,espco,esph,espno2
1862      real*8      rcmnz, rcmmini
1863      real*8      szadeg
1864
1865
1866      ! Tracer indexes in the thermospheric chemistry:
1867      !!! ATTENTION. These values have to be identical to those in chemthermos.F90
1868      !!! If the values are changed there, the same has to be done here  !!!
1869      integer,parameter :: i_co2=1
1870      integer,parameter :: i_o2=2
1871      integer,parameter :: i_o=3
1872      integer,parameter :: i_co=4
1873      integer,parameter :: i_h=5
1874      integer,parameter :: i_h2=8
1875      integer,parameter :: i_h2o=9
1876      integer,parameter :: i_h2o2=10
1877      integer,parameter :: i_o3=12
1878      integer,parameter :: i_n2=13
1879      integer,parameter :: i_n=14
1880      integer,parameter :: i_no=15
1881      integer,parameter :: i_no2=17
1882
1883
1884c*************************PROGRAM STARTS*******************************
1885
1886      nz3 = nlayermx*2
1887      szadeg = zenit*180./3.141592
1888      do i=1,nlayermx
1889         xx = ( radio + iz(i) ) * 1.e5
1890         grav(i) = gg * masa /(xx**2)
1891      end do
1892
1893      !Scale heights
1894      xx = kboltzman * tx(nlayermx) * n_avog / grav(nlayermx)
1895      Ho3p  = xx / mmol(igcm_o)
1896      Hco2  = xx / mmol(igcm_co2)
1897      Ho2   = xx / mmol(igcm_o2)
1898      Hh2   = xx / mmol(igcm_h2)
1899      Hh2o  = xx / mmol(igcm_h2o_vap)
1900      Hh2o2 = xx / mmol(igcm_h2o2)
1901      Hco   = xx / mmol(igcm_co)
1902      Hh    = xx / mmol(igcm_h)
1903      !Only if O3 chem. required
1904      if(chemthermod.ge.1)
1905     $     Ho3   = xx / mmol(igcm_o3)
1906      !Only if N or ion chem.
1907      if(chemthermod.ge.2) then
1908         Hn2   = xx / mmol(igcm_n2)
1909         Hn    = xx / mmol(igcm_n)
1910         Hno   = xx / mmol(igcm_no)
1911         Hno2  = xx / mmol(igcm_no2)
1912      endif
1913      do i=nlayermx,1,-1
1914         !Column initialisation
1915         co2colx(i)  = 0.
1916         o2colx(i)   = 0.
1917         o3pcolx(i)  = 0.
1918         h2colx(i)   = 0.
1919         h2ocolx(i)  = 0.
1920         h2o2colx(i) = 0.
1921         o3colx(i)   = 0.
1922         n2colx(i)   = 0.
1923         ncolx(i)    = 0.
1924         nocolx(i)   = 0.
1925         cocolx(i)   = 0.
1926         hcolx(i)    = 0.
1927         no2colx(i)  = 0.
1928         !Densities
1929         co2x(i)  = rm(i,i_co2)
1930         o2x(i)   = rm(i,i_o2)
1931         o3px(i)  = rm(i,i_o)
1932         h2x(i)   = rm(i,i_h2)
1933         h2ox(i)  = rm(i,i_h2o)
1934         h2o2x(i) = rm(i,i_h2o2)
1935         cox(i)   = rm(i,i_co)
1936         hx(i)    = rm(i,i_h)
1937         !Only if O3 chem. required
1938         if(chemthermod.ge.1)
1939     $        o3x(i)   = rm(i,i_o3)
1940         !Only if Nitrogen of ion chemistry requested
1941         if(chemthermod.ge.2) then
1942            n2x(i)   = rm(i,i_n2)
1943            nx(i)    = rm(i,i_n)
1944            nox(i)   = rm(i,i_no)
1945            no2x(i)  = rm(i,i_no2)
1946         endif
1947         !Routine to calculate the geometrical length of each layer
1948         call espesor_optico_A(ig,i,zenit,iz(i),nz3,iz,esp,ilayesp,
1949     $         szalayesp,nlayesp, zmini)
1950         if(ilayesp(nlayesp).eq.-1) then
1951            co2colx(i)=1.e25
1952            o2colx(i)=1.e25
1953            o3pcolx(i)=1.e25
1954            h2colx(i)=1.e25
1955            h2ocolx(i)=1.e25
1956            h2o2colx(i)=1.e25
1957            o3colx(i)=1.e25
1958            n2colx(i)=1.e25
1959            ncolx(i)=1.e25
1960            nocolx(i)=1.e25
1961            cocolx(i)=1.e25
1962            hcolx(i)=1.e25
1963            no2colx(i)=1.e25
1964         else
1965            rcmnz = ( radio + iz(nlayermx) ) * 1.e5
1966            rcmmini = ( radio + zmini ) * 1.e5
1967            !Column calculation taking into account the geometrical depth
1968            !calculated before
1969            do j=1,nlayesp
1970               jj=ilayesp(j)
1971               !Top layer
1972               if(jj.eq.nlayermx) then
1973                  if(szadeg.le.60.) then
1974                     o3pcolx(i)=o3pcolx(i)+o3px(nlayermx)*Ho3p*esp(j)
1975     $                    *1.e-5
1976                     co2colx(i)=co2colx(i)+co2x(nlayermx)*Hco2*esp(j)
1977     $                    *1.e-5
1978                     h2o2colx(i)=h2o2colx(i)+
1979     $                    h2o2x(nlayermx)*Hh2o2*esp(j)*1.e-5
1980                     o2colx(i)=o2colx(i)+o2x(nlayermx)*Ho2*esp(j)
1981     $                    *1.e-5
1982                     h2colx(i)=h2colx(i)+h2x(nlayermx)*Hh2*esp(j)
1983     $                    *1.e-5
1984                     h2ocolx(i)=h2ocolx(i)+h2ox(nlayermx)*Hh2o*esp(j)
1985     $                    *1.e-5                     
1986                     cocolx(i)=cocolx(i)+cox(nlayermx)*Hco*esp(j)
1987     $                    *1.e-5
1988                     hcolx(i)=hcolx(i)+hx(nlayermx)*Hh*esp(j)
1989     $                    *1.e-5
1990                     !Only if O3 chemistry required
1991                     if(chemthermod.ge.1) o3colx(i)=
1992     $                    o3colx(i)+o3x(nlayermx)*Ho3*esp(j)
1993     $                    *1.e-5
1994                     !Only if N or ion chemistry requested
1995                     if(chemthermod.ge.2) then
1996                        n2colx(i)=n2colx(i)+n2x(nlayermx)*Hn2*esp(j)
1997     $                    *1.e-5
1998                        ncolx(i)=ncolx(i)+nx(nlayermx)*Hn*esp(j)
1999     $                       *1.e-5
2000                        nocolx(i)=nocolx(i)+nox(nlayermx)*Hno*esp(j)
2001     $                       *1.e-5
2002                        no2colx(i)=no2colx(i)+no2x(nlayermx)*Hno2*esp(j)
2003     $                       *1.e-5
2004                     endif
2005                  else if(szadeg.gt.60.) then
2006                     espco2 =sqrt((rcmnz+Hco2)**2 -rcmmini**2) - esp(j)
2007                     espo2  = sqrt((rcmnz+Ho2)**2 -rcmmini**2) - esp(j)
2008                     espo3p = sqrt((rcmnz+Ho3p)**2 -rcmmini**2)- esp(j)
2009                     esph2  = sqrt((rcmnz+Hh2)**2 -rcmmini**2) - esp(j)
2010                     esph2o = sqrt((rcmnz+Hh2o)**2 -rcmmini**2)- esp(j)
2011                     esph2o2= sqrt((rcmnz+Hh2o2)**2-rcmmini**2)- esp(j)
2012                     espco  = sqrt((rcmnz+Hco)**2 -rcmmini**2) - esp(j)
2013                     esph   = sqrt((rcmnz+Hh)**2 -rcmmini**2)  - esp(j)
2014                     !Only if O3 chemistry required
2015                     if(chemthermod.ge.1)                     
2016     $                   espo3=sqrt((rcmnz+Ho3)**2-rcmmini**2)-esp(j)
2017                     !Only if N or ion chemistry requested
2018                     if(chemthermod.ge.2) then
2019                        espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j)
2020                        espn  =sqrt((rcmnz+Hn)**2-rcmmini**2)  - esp(j)
2021                        espno =sqrt((rcmnz+Hno)**2-rcmmini**2) - esp(j)
2022                        espno2=sqrt((rcmnz+Hno2)**2-rcmmini**2)- esp(j)
2023                     endif
2024                     co2colx(i) = co2colx(i) + espco2*co2x(nlayermx)
2025                     o2colx(i)  = o2colx(i)  + espo2*o2x(nlayermx)
2026                     o3pcolx(i) = o3pcolx(i) + espo3p*o3px(nlayermx)
2027                     h2colx(i)  = h2colx(i)  + esph2*h2x(nlayermx)
2028                     h2ocolx(i) = h2ocolx(i) + esph2o*h2ox(nlayermx)
2029                     h2o2colx(i)= h2o2colx(i)+ esph2o2*h2o2x(nlayermx)
2030                     cocolx(i)  = cocolx(i)  + espco*cox(nlayermx)
2031                     hcolx(i)   = hcolx(i)   + esph*hx(nlayermx)
2032                     !Only if O3 chemistry required
2033                     if(chemthermod.ge.1)                     
2034     $                  o3colx(i) = o3colx(i)  + espo3*o3x(nlayermx)
2035                     !Only if N or ion chemistry requested
2036                     if(chemthermod.ge.2) then
2037                        n2colx(i)  = n2colx(i)  + espn2*n2x(nlayermx)
2038                        ncolx(i)   = ncolx(i)   + espn*nx(nlayermx)
2039                        nocolx(i)  = nocolx(i)  + espno*nox(nlayermx)
2040                        no2colx(i) = no2colx(i) + espno2*no2x(nlayermx)
2041                     endif
2042                  endif !Of if szadeg.lt.60
2043               !Other layers
2044               else
2045                  co2colx(i)  = co2colx(i) +
2046     $                 esp(j) * (co2x(jj)+co2x(jj+1)) / 2.
2047                  o2colx(i)   = o2colx(i) +
2048     $                 esp(j) * (o2x(jj)+o2x(jj+1)) / 2.
2049                  o3pcolx(i)  = o3pcolx(i) +
2050     $                 esp(j) * (o3px(jj)+o3px(jj+1)) / 2.
2051                  h2colx(i)   = h2colx(i) +
2052     $                 esp(j) * (h2x(jj)+h2x(jj+1)) / 2.
2053                  h2ocolx(i)  = h2ocolx(i) +
2054     $                 esp(j) * (h2ox(jj)+h2ox(jj+1)) / 2.
2055                  h2o2colx(i) = h2o2colx(i) +
2056     $                 esp(j) * (h2o2x(jj)+h2o2x(jj+1)) / 2.
2057                  cocolx(i)   = cocolx(i) +
2058     $                 esp(j) * (cox(jj)+cox(jj+1)) / 2.
2059                  hcolx(i)    = hcolx(i) +
2060     $                 esp(j) * (hx(jj)+hx(jj+1)) / 2.
2061                  !Only if O3 chemistry required
2062                  if(chemthermod.ge.1)
2063     $                 o3colx(i) = o3colx(i) +
2064     $                 esp(j) * (o3x(jj)+o3x(jj+1)) / 2.
2065                  !Only if N or ion chemistry requested
2066                  if(chemthermod.ge.2) then
2067                     n2colx(i)   = n2colx(i) +
2068     $                 esp(j) * (n2x(jj)+n2x(jj+1)) / 2.
2069                     ncolx(i)    = ncolx(i) +
2070     $                    esp(j) * (nx(jj)+nx(jj+1)) / 2.
2071                     nocolx(i)   = nocolx(i) +
2072     $                    esp(j) * (nox(jj)+nox(jj+1)) / 2.
2073                     no2colx(i)  = no2colx(i) +
2074     $                    esp(j) * (no2x(jj)+no2x(jj+1)) / 2.
2075                  endif
2076               endif  !Of if jj.eq.nlayermx
2077            end do    !Of do j=1,nlayesp
2078         endif        !Of ilayesp(nlayesp).eq.-1
2079      enddo           !Of do i=nlayermx,1,-1
2080      return
2081
2082
2083      end
2084
2085
2086c**********************************************************************
2087c**********************************************************************
2088      subroutine interpfast(escout,p,nlayer,escin,pin,nl,limdown,limup)
2089C
2090C subroutine to perform linear interpolation in pressure from 1D profile
2091C escin(nl) sampled on pressure grid pin(nl) to profile
2092C escout(nlayer) on pressure grid p(nlayer).
2093C
2094      real escout(nlayer),p(nlayer)
2095      real escin(nl),pin(nl),wm,wp
2096      real limup,limdown
2097      integer nl,nlayer,n1,n,nm,np
2098      nm=1
2099      do 5 n1=1,nlayer
2100         if(p(n1) .gt. limup .or. p(n1) .lt. limdown) then
2101            escout(n1) = 0.0
2102         else
2103            do n = nm,nl-1
2104               if (p(n1).ge.pin(n).and.p(n1).le.pin(n+1)) then
2105                  nm=n
2106                  np=n+1
2107                  wm=abs(pin(nm)-p(n1))/(pin(np)-pin(nm))
2108                  wp=1.0 - wm
2109                  goto 33
2110               endif
2111            enddo
2112 33         escout(n1) = escin(np)*wm + escin(nm)*wp
2113         endif
2114    5 continue
2115      return
2116      end
2117
2118
2119
2120c**********************************************************************
2121c**********************************************************************
2122
2123      subroutine espesor_optico_A (ig,capa, szadeg,z,
2124     @                   nz3,iz,esp,ilayesp,szalayesp,nlayesp, zmini)
2125
2126c     fgg              nov 03      Adaptation to Martian model
2127c     malv             jul 03      Corrected z grid. Split in alt & frec codes
2128c     fgg              feb 03      first version
2129*************************************************************************
2130
2131      implicit none
2132
2133
2134c     common variables and constants
2135
2136      include "dimensions.h"
2137      include "dimphys.h"
2138      include 'param.h'
2139      include 'param_v4.h'
2140
2141c     arguments
2142
2143      real        szadeg                ! I. SZA [rad]
2144      real        z                     ! I. altitude of interest [km]
2145      integer     nz3,ig                   ! I. dimension of esp, ylayesp, etc...
2146                                        !  (=2*nlayermx= max# of layers in ray path)
2147      real     iz(nlayermx+1)              ! I. Altitude of each layer
2148      real*8        esp(nz3)            ! O. layer widths after geometrically
2149                                        !    amplified; in [cm] except at TOA
2150                                        !    where an auxiliary value is used
2151      real*8        ilayesp(nz3)        ! O. Indexes of layers along ray path
2152      real*8        szalayesp(nz3)      ! O. SZA [deg]    "     "       "
2153      integer       nlayesp
2154!      real*8        nlayesp             ! O. # layers along ray path at this z
2155      real*8        zmini               ! O. Minimum altitud of ray path [km]
2156
2157
2158c     local variables and constants
2159
2160        integer     j,i,capa
2161        integer     jmin                  ! index of min.altitude along ray path
2162        real*8      szarad                ! SZA [deg]
2163        real*8      zz
2164        real*8      diz(nlayermx+1)
2165        real*8      rkmnz                 ! distance TOA to center of Planet [km]
2166        real*8      rkmmini               ! distance zmini to center of P [km]
2167        real*8      rkmj                  ! intermediate distance to C of P [km]
2168
2169c external function
2170        external  grid_R8          ! Returns index of layer containing the altitude
2171                                ! of interest, z; for example, if
2172                                ! zkm(i)=z or zkm(i)<z<zkm(i+1) => grid(z)=i
2173        integer   grid_R8
2174
2175*************************************************************************     
2176        szarad = dble(szadeg)*3.141592d0/180.d0
2177        zz=dble(z)
2178        do i=1,nlayermx
2179           diz(i)=dble(iz(i))
2180        enddo
2181        do j=1,nz3
2182          esp(j) = 0.d0
2183          szalayesp(j) = 777.d0
2184          ilayesp(j) = 0
2185        enddo
2186        nlayesp = 0
2187
2188        ! First case: szadeg<60
2189        ! The optical thickness will be given by  1/cos(sza)
2190        ! We deal with 2 different regions:
2191        !   1: First, all layers between z and ztop ("upper part of ray")
2192        !   2: Second, the layer at ztop
2193        if(szadeg.lt.60.d0) then
2194
2195           zmini = zz
2196           if(abs(zz-diz(nlayermx)).lt.1.d-3) goto 1357
2197           ! 1st Zone: Upper part of ray
2198           !
2199           do j=grid_R8(zz,diz,nlayermx),nlayermx-1
2200             nlayesp = nlayesp + 1
2201             ilayesp(nlayesp) = j
2202             esp(nlayesp) = (diz(j+1)-diz(j)) / cos(szarad)        ! [km]
2203             esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
2204             szalayesp(nlayesp) = szadeg
2205           end do
2206
2207           !
2208           ! 2nd Zone: Top layer
2209 1357      continue
2210           nlayesp = nlayesp + 1
2211           ilayesp(nlayesp) = nlayermx
2212           esp(nlayesp) = 1.d0 / cos(szarad)         ! aux. non-dimens. factor
2213           szalayesp(nlayesp) = szadeg
2214
2215
2216        ! Second case:  60 < szadeg < 90
2217        ! The optical thickness is evaluated.
2218        !    (the magnitude of the effect of not using cos(sza) is 3.e-5
2219        !     for z=60km & sza=30, and 5e-4 for z=60km & sza=60, approximately)
2220        ! We deal with 2 different regions:
2221        !   1: First, all layers between z and ztop ("upper part of ray")
2222        !   2: Second, the layer at ztop ("uppermost layer")
2223        else if(szadeg.le.90.d0.and.szadeg.ge.60.d0) then
2224
2225           zmini=(radio+zz)*sin(szarad)-radio
2226           rkmmini = radio + zmini
2227
2228           if(abs(zz-diz(nlayermx)).lt.1.d-4) goto 1470
2229
2230           ! 1st Zone: Upper part of ray
2231           !
2232           do j=grid_R8(zz,diz,nlayermx),nlayermx-1
2233              nlayesp = nlayesp + 1
2234              ilayesp(nlayesp) = j
2235              esp(nlayesp) =
2236     #             sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) -
2237     #             sqrt( (radio+diz(j))**2 - rkmmini**2 )           ! [km]
2238              esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
2239              rkmj = radio+diz(j)
2240              szalayesp(nlayesp) = asin( rkmmini/rkmj )             ! [rad]
2241              szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592 ! [deg]
2242           end do
2243 1470      continue
2244           ! 2nd Zone:  Uppermost layer of ray.
2245           !
2246           nlayesp = nlayesp + 1
2247           ilayesp(nlayesp) = nlayermx
2248           rkmnz = radio+diz(nlayermx)
2249           esp(nlayesp)  =  sqrt( rkmnz**2 - rkmmini**2 )       ! aux.factor[km]
2250           esp(nlayesp)  =  esp(nlayesp) * 1.d5                 ! aux.f. [cm]
2251           szalayesp(nlayesp) = asin( rkmmini/rkmnz )           ! [rad]
2252           szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592! [deg]
2253
2254
2255        ! Third case:  szadeg > 90
2256        ! The optical thickness is evaluated.
2257        ! We deal with 5 different regions:
2258        !   1: all layers between z and ztop ("upper part of ray")
2259        !   2: the layer at ztop ("uppermost layer")
2260        !   3: the lowest layer, at zmini
2261        !   4: the layers increasing from zmini to z (here SZA<90)
2262        !   5: the layers decreasing from z to zmini (here SZA>90)
2263        else if(szadeg.gt.90.d0) then
2264
2265           zmini=(radio+zz)*sin(szarad)-radio
2266           rkmmini = radio + zmini
2267
2268           if(zmini.lt.diz(1)) then         ! Can see the sun?  No => esp(j)=inft
2269             nlayesp = nlayesp + 1
2270             ilayesp(nlayesp) = - 1     ! Value to mark "no sun on view"
2271!             esp(nlayesp) = 1.e30
2272
2273           else
2274              jmin=grid_R8(zmini,diz,nlayermx)+1
2275             
2276
2277              if(abs(zz-diz(nlayermx)).lt.1.d-4) goto 9876
2278
2279              ! 1st Zone: Upper part of ray
2280              !
2281              do j=grid_R8(zz,diz,nlayermx),nlayermx-1
2282                nlayesp = nlayesp + 1
2283                ilayesp(nlayesp) = j
2284                esp(nlayesp) =
2285     $                sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) -
2286     $                sqrt( (radio+diz(j))**2 - rkmmini**2 )          ! [km]
2287                esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
2288                rkmj = radio+diz(j)
2289                szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
2290                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592      ! [deg]
2291              end do
2292
2293 9876         continue
2294              ! 2nd Zone:  Uppermost layer of ray.
2295              !
2296              nlayesp = nlayesp + 1
2297              ilayesp(nlayesp) = nlayermx
2298              rkmnz = radio+diz(nlayermx)
2299              esp(nlayesp) =  sqrt( rkmnz**2 - rkmmini**2 )      !aux.factor[km]
2300              esp(nlayesp) = esp(nlayesp) * 1.d5                 !aux.f.[cm]
2301              szalayesp(nlayesp) = asin( rkmmini/rkmnz )           ! [rad]
2302              szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
2303
2304              ! 3er Zone: Lowestmost layer of ray
2305              !
2306              if ( jmin .ge. 2 ) then      ! If above the planet's surface
2307                j=jmin-1
2308                nlayesp = nlayesp + 1
2309                ilayesp(nlayesp) = j
2310                esp(nlayesp) = 2. *
2311     $                 sqrt( (radio+diz(j+1))**2 -rkmmini**2 )       ! [km]
2312                esp(nlayesp) = esp(nlayesp) * 1.d5                   ! [cm]
2313                rkmj = radio+diz(j+1)
2314                szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad]
2315                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
2316              endif
2317
2318              ! 4th zone: Lower part of ray, increasing from zmin to z
2319              !    ( layers with SZA < 90 deg )
2320              do j=jmin,grid_R8(zz,diz,nlayermx)-1
2321                nlayesp = nlayesp + 1
2322                ilayesp(nlayesp) = j
2323                esp(nlayesp) =
2324     $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
2325     $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )       ! [km]
2326                esp(nlayesp) = esp(nlayesp) * 1.d5                     ! [cm]
2327                rkmj = radio+diz(j)
2328                szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
2329                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
2330              end do
2331
2332              ! 5th zone: Lower part of ray, decreasing from z to zmin
2333              !    ( layers with SZA > 90 deg )
2334              do j=grid_R8(zz,diz,nlayermx)-1, jmin, -1
2335                nlayesp = nlayesp + 1
2336                ilayesp(nlayesp) = j
2337                esp(nlayesp) =
2338     $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
2339     $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )        ! [km]
2340                esp(nlayesp) = esp(nlayesp) * 1.d5                      ! [cm]
2341                rkmj = radio+diz(j)
2342                szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj )          ! [rad]
2343                szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592 ! [deg]
2344              end do
2345
2346           end if
2347
2348        end if
2349
2350        return
2351
2352        end
2353
2354
2355
2356c**********************************************************************
2357c***********************************************************************
2358
2359        function grid_R8 (z, zgrid, nz)
2360
2361c Returns the index where z is located within vector zgrid
2362c The vector zgrid must be monotonously increasing, otherwise program stops.
2363c If z is outside zgrid limits, or zgrid dimension is nz<2, the program stops.
2364c
2365c FGG     Aug-2004     Correct z.lt.zgrid(i) to .le.
2366c MALV    Jul-2003
2367c***********************************************************************
2368
2369        implicit none
2370
2371c Arguments
2372        integer   nz
2373        real*8      z
2374        real*8      zgrid(nz)
2375        integer   grid_R8
2376
2377c Local 
2378        integer   i, nz1, nznew
2379
2380c*** CODE START
2381
2382        if ( z .lt. zgrid(1) .or. z.gt.zgrid(nz) ) then
2383           write (*,*) ' GRID/ z outside bounds of zgrid '
2384           write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz)
2385           stop ' Serious error in GRID.F '
2386        endif
2387        if ( nz .lt. 2 ) then
2388           write (*,*) ' GRID/ zgrid needs 2 points at least ! '
2389           stop ' Serious error in GRID.F '
2390        endif
2391        if ( zgrid(1) .ge. zgrid(nz) ) then
2392           write (*,*) ' GRID/ zgrid must increase with index'
2393           stop ' Serious error in GRID.F '
2394        endif
2395
2396        nz1 = 1
2397        nznew = nz/2
2398        if ( z .gt. zgrid(nznew) ) then
2399           nz1 = nznew
2400           nznew = nz
2401        endif
2402        do i=nz1+1,nznew
2403           if ( z. eq. zgrid(i) ) then
2404              grid_R8=i
2405              return
2406              elseif ( z .le. zgrid(i) ) then
2407              grid_R8 = i-1
2408              return
2409           endif
2410        enddo
2411        grid_R8 = nz
2412        return
2413
2414
2415
2416        end
2417
2418
2419
2420!c***************************************************
2421!c***************************************************
2422
2423      subroutine flujo(date)
2424
2425
2426!c     fgg           nov 2002     first version
2427!c***************************************************
2428
2429      implicit none
2430
2431
2432!     common variables and constants
2433      include "dimensions.h"
2434      include "dimphys.h"
2435      include "comsaison.h"
2436      include 'param.h'
2437      include 'param_v4.h'
2438
2439
2440!     Arguments
2441
2442      real date
2443
2444
2445!     Local variable and constants
2446
2447      integer i
2448      integer inter
2449      real    nada
2450
2451!c*************************************************
2452
2453      if(date.lt.1985.) date=1985.
2454      if(date.gt.2001.) date=2001.
2455     
2456      do i=1,ninter
2457         fluxtop(i)=1.
2458         !Variation of solar flux with 11 years solar cycle
2459         !For more details, see Gonzalez-Galindo et al. 2005
2460         !To be improved in next versions
2461        if(i.le.24) then
2462          fluxtop(i)=(((ct1(i)+p1(i)*date)/2.)                 
2463     $        *sin(2.*3.1416/11.*(date-1985.-3.1416))         
2464     $        +(ct2(i)+p2(i)*date)+1.)*fluxtop(i)
2465        end if
2466        fluxtop(i)=fluxtop(i)*(1.52/dist_sol)**2
2467      end do
2468!      fluxtop(1)=fluxtop(1)*10.
2469!      fluxtop(2)=fluxtop(2)*5.
2470     
2471      return
2472      end
Note: See TracBrowser for help on using the repository browser.