Ignore:
Timestamp:
May 12, 2012, 8:10:08 AM (13 years ago)
Author:
emillour
Message:

Mars GCM:

  • updated high atmosphere photochemistry (jthermcalc.F, param_v4.h, iono.h, paramfoto_compact.F, param_read.F , thermosphere.F).
  • minor change in calchim.F90 (to not use maxloc(zycol, dim = 2) function which seems to be a problem for g95) .
  • minor bug fix in perosat.F; set tendency on pdqscloud for h2o2 to zero if none is computed.
  • in "moldiff.F", changed "tridag" to "tridag_sp", "LUBKSB" to "LUBKSB_SP" and "LUDCMP" to "LUDCMP_SP" to avoid possible conflicts with same routines defined in "moldiff_red.F". Added use of automatic-sized array in "tridag" and "tridag_sp" local array "gam" (to avoid using an a priori oversized local array).

FGG+JYC+EM

Location:
trunk/LMDZ.MARS/libf/aeronomars
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/aeronomars/calchim.F90

    r635 r658  
    9595      integer,save :: nbq                   ! number of tracers used in the chemistry
    9696      integer,save :: niq(nqmx)             ! array storing the indexes of the tracers
    97       integer :: major(nlayermx)            ! index of major species
     97      integer :: iloc(1)            ! index of major species
    9898      integer :: ig,l,i,iq,iqmax
    9999      integer :: foundswitch, lswitch
    … …  
    171171            call read_phototable
    172172         end if
    173 
    174173         ! find index of chemical tracers to use
    175174         ! Listed here are all tracers that can go into photochemistry
    … …  
    660659                            zu, zv, zt, zycol, ptimestep, co2ice)
    661660         end if
    662 
    663661!=======================================================================
    664662!     tendencies
    … …  
    667665!     index of the most abundant species at each level
    668666
    669          major(:) = maxloc(zycol, dim = 2)
     667!         major(:) = maxloc(zycol, dim = 2)
    670668
    671669!     tendency for the most abundant species = - sum of others
    672 
    673670         do l = 1,nlayermx
    674             iqmax = major(l)
     671            iloc=maxloc(zycol(l,:))
     672            iqmax=iloc(1)
    675673            do i = 1,nbq
    676674               iq = niq(i) ! get tracer index
    … …  
    683681            end do
    684682         end do ! of do l = 1,nlayermx
    685 !         if(ig.eq.800)write(*,*)'calchim/686',dqchim(ig,27,23)
    686683
    687684!=======================================================================

  • trunk/LMDZ.MARS/libf/aeronomars/iono.h

    r635 r658  
    1919     $ taunoplus,taun2plus,tauhplus,tauhco2plus
    2020
    21       real*8    tauco2(nlayermx,nreact)
    22       real*8    tauo2(nlayermx,nreact)
    23       real*8    tauo3p(nlayermx,nreact)
    24       real*8    tauco(nlayermx,nreact)
    25       real*8    tauh(nlayermx,nreact)
    26       real*8    tauoh(nlayermx,nreact)
    27       real*8    tauho2(nlayermx,nreact)
    28       real*8    tauh2(nlayermx,nreact)
    29       real*8    tauh2o(nlayermx,nreact)
    30       real*8    tauo1d(nlayermx,nreact)
    31       real*8    tauh2o2(nlayermx,nreact)
    32       real*8    tauo3(nlayermx,nreact)
    33       real*8    taun(nlayermx,nreact)
    34       real*8    tauno(nlayermx,nreact)
    35       real*8    taun2(nlayermx,nreact)
    36       real*8    taun2d(nlayermx,nreact)
    37       real*8    tauno2(nlayermx,nreact)
    38       real*8    tauco2plus(nlayermx,nreact)
    39       real*8    tauoplus(nlayermx,nreact)
    40       real*8    tauo2plus(nlayermx, nreact)
    41       real*8    taucoplus(nlayermx, nreact)
    42       real*8    taucplus(nlayermx, nreact)
    43       real*8    taunplus(nlayermx, nreact)
    44       real*8    taunoplus(nlayermx, nreact)
    45       real*8    taun2plus(nlayermx, nreact)
    46       real*8    tauhplus(nlayermx, nreact)
    47       real*8    tauhco2plus(nlayermx, nreact)
     21      real*8    tauco2(nreact,nlayermx)
     22      real*8    tauo2(nreact,nlayermx)
     23      real*8    tauo3p(nreact,nlayermx)
     24      real*8    tauco(nreact,nlayermx)
     25      real*8    tauh(nreact,nlayermx)
     26      real*8    tauoh(nreact,nlayermx)
     27      real*8    tauho2(nreact,nlayermx)
     28      real*8    tauh2(nreact,nlayermx)
     29      real*8    tauh2o(nreact,nlayermx)
     30      real*8    tauo1d(nreact,nlayermx)
     31      real*8    tauh2o2(nreact,nlayermx)
     32      real*8    tauo3(nreact,nlayermx)
     33      real*8    taun(nreact,nlayermx)
     34      real*8    tauno(nreact,nlayermx)
     35      real*8    taun2(nreact,nlayermx)
     36      real*8    taun2d(nreact,nlayermx)
     37      real*8    tauno2(nreact,nlayermx)
     38      real*8    tauco2plus(nreact,nlayermx)
     39      real*8    tauoplus(nreact,nlayermx)
     40      real*8    tauo2plus(nreact,nlayermx)
     41      real*8    taucoplus(nreact,nlayermx)
     42      real*8    taucplus(nreact,nlayermx)
     43      real*8    taunplus(nreact,nlayermx)
     44      real*8    taunoplus(nreact,nlayermx)
     45      real*8    taun2plus(nreact,nlayermx)
     46      real*8    tauhplus(nreact,nlayermx)
     47      real*8    tauhco2plus(nreact,nlayermx)

  • trunk/LMDZ.MARS/libf/aeronomars/jthermcalc.F

    r635 r658  
    3232c    local parameters and variables
    3333
    34       real       aux1(nlayermx)                 !auxiliar variable
    35       real       aux2(nlayermx)                 !auxiliar variable
    3634      real       co2colx(nlayermx)              !column density of CO2 (cm^-2)
    3735      real       o2colx(nlayermx)               !column density of O2(cm^-2)
    … …  
    5755
    5856
    59 
    6057c     variables used in interpolation
    6158
    62       real aux3(nz2), aux4(nz2)
    63                                               !auxiliar variables
    64       real       limdown                      !limits for interpolation
    65       real       limup                        !        ""
     59      real*8      auxcoltab(nz2)
     60      real*8      auxjco2(nz2)
     61      real*8      auxjo2(nz2)
     62      real*8      auxjo3p(nz2)
     63      real*8      auxjh2o(nz2)
     64      real*8      auxjh2(nz2)
     65      real*8      auxjh2o2(nz2)
     66      real*8      auxjo3(nz2)
     67      real*8      auxjn2(nz2)
     68      real*8      auxjn(nz2)
     69      real*8      auxjno(nz2)
     70      real*8      auxjco(nz2)
     71      real*8      auxjh(nz2)
     72      real*8      auxjno2(nz2)
     73      real*8      wp(nlayermx),wm(nlayermx)
     74      real*8      auxcolinp(nlayermx)
     75      integer     auxind(nlayermx)
     76      integer     auxi
     77      integer     ind
     78      real*8      cortemp(nlayermx)
     79
     80      real*8      limdown                      !limits for interpolation
     81      real*8      limup                        !        ""
    6682
    6783      !!!ATTENTION. Here i_co2 has to have the same value than in chemthermos.F90
    6884      !!!If the value is changed there, if has to be changed also here !!!!
    6985      integer,parameter :: i_co2=1
     86
    7087
    7188c*************************PROGRAM STARTS*******************************
    … …  
    117134! in each spectral interval
    118135
    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
    135 c************************************************
    136 c     co2,0.1,6.0
    137 c************************************************
    138 
     136
     137c     auxcolinp-> Actual atmospheric column
     138c     auxj*-> Tabulated photoabsorption coefficients
     139c     auxcoltab-> Tabulated atmospheric columns
     140
     141ccccccccccccccccccccccccccccccc
     142c     0.1,5.0 (int 1)
     143c
     144c     Absorption by:
     145c     CO2, O2, O, H2, N
     146ccccccccccccccccccccccccccccccc
     147
     148c     Input atmospheric column
    139149      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) +
     150      do i=1,nlayermx
     151         auxcolinp(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint) +
    145152     $        o2colx(i)*crscabsi2(2,indexint) +
    146153     $        o3pcolx(i)*crscabsi2(3,indexint) +
    … …  
    148155     $        ncolx(i)*crscabsi2(9,indexint)
    149156      end do
     157      limdown=1.e-20
     158      limup=1.e26
     159
     160
     161c     Interpolations
     162
    150163      do i=1,nz2
    151          aux3(i) = jabsifotsintpar(indexint,1,nz2-i+1)
    152          aux4(i) = c1_16(nz2-i+1,indexint)
     164         auxi = nz2-i+1
     165         !CO2 tabulated coefficient
     166         auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
     167         !O2 tabulated coefficient
     168         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
     169         !O3p tabulated coefficient
     170         auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
     171         !H2 tabulated coefficient
     172         auxjh2(i) = jabsifotsintpar(auxi,5,indexint)
     173         !Tabulated column
     174         auxcoltab(i) = c1_16(auxi,indexint)
    153175      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)
     176      !Only if chemthermod.ge.2
     177      !N tabulated coefficient
     178      if(chemthermod.ge.2) then
     179         do i=1,nz2
     180            auxjn(i) = jabsifotsintpar(nz2-i+1,9,indexint)
     181         enddo
     182      endif
     183
     184      call interfast
     185     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
     186      do i=1,nlayermx
     187         ind=auxind(i)
     188         auxi=nlayermx-i+1
     189         !CO2 interpolated coefficient
     190         jfotsout(indexint,1,auxi) = wm(i)*auxjco2(ind+1) +
     191     $        wp(i)*auxjco2(ind)
     192         !O2 interpolated coefficient
     193         jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
     194     $        wp(i)*auxjo2(ind)
     195         !O3p interpolated coefficient
     196          jfotsout(indexint,3,auxi) = wm(i)*auxjo3p(ind+1) +
     197     $        wp(i)*auxjo3p(ind)
     198          !H2 interpolated coefficient
     199          jfotsout(indexint,5,auxi) = wm(i)*auxjh2(ind+1) +
     200     $         wp(i)*auxjh2(ind)
    159201      enddo
    160 
    161 c************************************************
    162 c     o2,0.1,6.0
    163 c************************************************
    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 
    187 c**************************************************         
    188 c     o3p,0.1,6.0
    189 c**************************************************
    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      
    212 c**************************************************
    213 c     h2,0.1,6.0
    214 c**************************************************
    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
     202      !Only if chemthermod.ge.2
     203      !N interpolated coefficient
    237204      if(chemthermod.ge.2) then
    238 
    239 c**************************************************
    240 c     n,0.1,6.0
    241 c**************************************************
    242 
    243 
    244          indexint=1
    245          limdown=1e-20
    246          limup=1e26
    247205         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)
     206            ind=auxind(i)
     207            jfotsout(indexint,9,nlayermx-i+1) =  wm(i)*auxjn(ind+1) +
     208     $         wp(i)*auxjn(ind)
     209         enddo
     210      endif
     211         
     212
     213c     End interval 1
     214
     215
     216ccccccccccccccccccccccccccccccc
     217c     5-80.5nm (int 2-15)
     218c
     219c     Absorption by:
     220c     CO2, O2, O, H2, N2, N,
     221c     NO, CO, H, NO2
     222ccccccccccccccccccccccccccccccc
     223
     224c     Input atmospheric column
     225      do indexint=2,15
    260226         do i=1,nlayermx
    261             jfotsout(indexint,9,i) = aux1(nlayermx-i+1)
    262          enddo
    263 
    264       endif  !Of chemthermod >= 2
    265      
    266 
    267 c**************************************************
    268 c     o2, 5-80.5nm
    269 c**************************************************
    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)+
     227            auxcolinp(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
    278228     $           o2colx(i)*crscabsi2(2,indexint)+
    279229     $           o3pcolx(i)*crscabsi2(3,indexint)+
    … …  
    286236     $           no2colx(i)*crscabsi2(13,indexint)
    287237         end do
     238
     239c     Interpolations
     240
    288241         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)
     242            auxi = nz2-i+1
     243            !O2 tabulated coefficient
     244            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
     245            !O3p tabulated coefficient
     246            auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
     247            !CO2 tabulated coefficient
     248            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
     249            !H2 tabulated coefficient
     250            auxjh2(i) = jabsifotsintpar(auxi,5,indexint)
     251            !N2 tabulated coefficient
     252            auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
     253            !CO tabulated coefficient
     254            auxjco(i) = jabsifotsintpar(auxi,11,indexint)
     255            !H tabulated coefficient
     256            auxjh(i) = jabsifotsintpar(auxi,12,indexint)
     257            !tabulated column
     258            auxcoltab(i) = c1_16(auxi,indexint)
     259         enddo
     260         !Only if chemthermod.ge.2
     261         if(chemthermod.ge.2) then
     262            do i=1,nz2
     263               auxi = nz2-i+1
     264               !N tabulated coefficient
     265               auxjn(i) = jabsifotsintpar(auxi,9,indexint)
     266               !NO tabulated coefficient
     267               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
     268               !NO2 tabulated coefficient
     269               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
     270            enddo
     271         endif
     272
     273          call interfast(wm,wp,auxind,auxcolinp,nlayermx,
     274     $        auxcoltab,nz2,limdown,limup)
    294275          do i=1,nlayermx
    295             jfotsout(indexint,2,i) = aux1(nlayermx-i+1)
     276             ind=auxind(i)
     277             auxi = nlayermx-i+1
     278             !O2 interpolated coefficient
     279             jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
     280     $            wp(i)*auxjo2(ind)
     281             !O3p interpolated coefficient
     282             jfotsout(indexint,3,auxi) = wm(i)*auxjo3p(ind+1) +
     283     $            wp(i)*auxjo3p(ind)
     284             !CO2 interpolated coefficient
     285             jfotsout(indexint,1,auxi) = wm(i)*auxjco2(ind+1) +
     286     $            wp(i)*auxjco2(ind)
     287             !H2 interpolated coefficient
     288             jfotsout(indexint,5,auxi) = wm(i)*auxjh2(ind+1) +
     289     $            wp(i)*auxjh2(ind)
     290             !N2 interpolated coefficient
     291             jfotsout(indexint,8,auxi) = wm(i)*auxjn2(ind+1) +
     292     $            wp(i)*auxjn2(ind)             
     293             !CO interpolated coefficient
     294             jfotsout(indexint,11,auxi) = wm(i)*auxjco(ind+1) +
     295     $            wp(i)*auxjco(ind)
     296             !H interpolated coefficient
     297             jfotsout(indexint,12,auxi) = wm(i)*auxjh(ind+1) +
     298     $            wp(i)*auxjh(ind)             
    296299          enddo
    297       end do
    298      
    299 
    300 c**************************************************
    301 c     o3p, 5-80.5nm
    302 c**************************************************
    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 
    328 c**************************************************
    329 c     co2, 5-80.5nm
    330 c**************************************************
    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 
    357 c**************************************************
    358 c     h2, 5-80.5nm
    359 c**************************************************
    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 
    386 c**************************************************
    387 c     n2, 5-80.5nm
    388 c**************************************************
    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 
    415 c**************************************************
    416 c     co, 5-80.5nm
    417 c**************************************************
    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 
    444 c**************************************************
    445 c     h, 5-80.5nm
    446 c**************************************************
    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
    475 c**************************************************
    476 c     n, 5-80.5nm
    477 c**************************************************
    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 
    504 c**************************************************
    505 c     no, 5-80.5nm
    506 c**************************************************
    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 
    532 c**************************************************
    533 c     no2, 5-80.5nm
    534 c**************************************************
    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          
    563 c**************************************************
    564 c     o2, 80.6-90.8nm
    565 c**************************************************
    566 
     300          !Only if chemthermod.ge.2
     301          if(chemthermod.ge.2) then
     302             do i=1,nlayermx
     303                ind=auxind(i)
     304                auxi = nlayermx-i+1
     305                !N interpolated coefficient
     306                jfotsout(indexint,9,auxi) = wm(i)*auxjn(ind+1) +
     307     $               wp(i)*auxjn(ind)
     308                !NO interpolated coefficient
     309                jfotsout(indexint,10,auxi)=wm(i)*auxjno(ind+1) +
     310     $               wp(i)*auxjno(ind)
     311                !NO2 interpolated coefficient
     312                jfotsout(indexint,13,auxi)=wm(i)*auxjno2(ind+1)+
     313     $               wp(i)*auxjno2(ind)
     314             enddo
     315          endif   
     316      end do
     317
     318c     End intervals 2-15
     319
     320
     321ccccccccccccccccccccccccccccccc
     322c     80.6-90.8nm (int16)
     323c
     324c     Absorption by:
     325c     CO2, O2, O, N2, N, NO,
     326c     CO, H, NO2
     327ccccccccccccccccccccccccccccccc
     328
     329c     Input atmospheric column
    567330      indexint=16
    568331      do i=1,nlayermx
    569          aux2(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
     332         auxcolinp(nlayermx-i+1) = co2colx(i)*crscabsi2(1,indexint)+
    570333     $        o2colx(i)*crscabsi2(2,indexint)+
    571334     $        o3pcolx(i)*crscabsi2(3,indexint)+
    … …  
    577340     $        no2colx(i)*crscabsi2(13,indexint)
    578341      end do
     342
     343c     Interpolations
     344
    579345      do i=1,nz2
    580          aux3(i) = jabsifotsintpar(indexint,2,nz2-i+1)
    581          aux4(i) = c1_16(nz2-i+1,indexint)
     346         auxi = nz2-i+1
     347         !O2 tabulated coefficient
     348         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
     349         !CO2 tabulated coefficient
     350         auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
     351         !O3p tabulated coefficient
     352         auxjo3p(i) = jabsifotsintpar(auxi,3,indexint)
     353         !N2 tabulated coefficient
     354         auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
     355         !CO tabulated coefficient
     356         auxjco(i) = jabsifotsintpar(auxi,11,indexint)
     357         !H tabulated coefficient
     358         auxjh(i) = jabsifotsintpar(auxi,12,indexint)
     359         !NO2 tabulated coefficient
     360         auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
     361         !Tabulated column
     362         auxcoltab(i) = c1_16(auxi,indexint)
    582363      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)
     364      !Only if chemthermod.ge.2
     365      if(chemthermod.ge.2) then
     366         do i=1,nz2
     367            auxi = nz2-i+1
     368            !N tabulated coefficient
     369            auxjn(i) = jabsifotsintpar(auxi,9,indexint)
     370            !NO tabulated coefficient
     371            auxjno(i) = jabsifotsintpar(auxi,10,indexint)
     372            !NO2 tabulated coefficient
     373            auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
     374         enddo
     375      endif
     376
     377      call interfast
     378     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
     379      do i=1,nlayermx
     380         ind=auxind(i)
     381         auxi = nlayermx-i+1
     382         !O2 interpolated coefficient
     383         jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
     384     $            wp(i)*auxjo2(ind)
     385         !CO2 interpolated coefficient
     386         jfotsout(indexint,1,auxi) = wm(i)*auxjco2(ind+1) +
     387     $            wp(i)*auxjco2(ind)
     388         !O3p interpolated coefficient
     389         jfotsout(indexint,3,auxi) = wm(i)*auxjo3p(ind+1) +
     390     $            wp(i)*auxjo3p(ind)
     391         !N2 interpolated coefficient
     392         jfotsout(indexint,8,auxi) = wm(i)*auxjn2(ind+1) +
     393     $            wp(i)*auxjn2(ind)
     394         !CO interpolated coefficient
     395         jfotsout(indexint,11,auxi) = wm(i)*auxjco(ind+1) +
     396     $            wp(i)*auxjco(ind) 
     397         !H interpolated coefficient
     398         jfotsout(indexint,12,auxi) = wm(i)*auxjh(ind+1) +
     399     $            wp(i)*auxjh(ind)         
    587400      enddo
    588 
    589 
    590 c**************************************************
    591 c     co2, 80.6-90.8nm
    592 c**************************************************
    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 
    617 c**************************************************
    618 c     o3p, 80.6-90.8nm
    619 c**************************************************
    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 
    644 c**************************************************
    645 c     n2, 80.6-90.8nm
    646 c**************************************************
    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 
    671 c**************************************************
    672 c     co, 80.6-90.8nm
    673 c************************************************************
    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 
    698 c**************************************************
    699 c     h, 80.6-90.8nm
    700 c**************************************************
    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
     401      !Only if chemthermod.ge.2
    726402      if(chemthermod.ge.2) then
    727 
    728 c**************************************************
    729 c     n, 80.6-90.8nm
    730 c**************************************************
    731 
    732          indexint=16
    733403         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
     404            ind=auxind(i)
     405            auxi = nlayermx-i+1
     406            !N interpolated coefficient
     407            jfotsout(indexint,9,auxi) = wm(i)*auxjn(ind+1) +
     408     $           wp(i)*auxjn(ind)
     409            !NO interpolated coefficient
     410            jfotsout(indexint,10,auxi) = wm(i)*auxjno(ind+1) +
     411     $           wp(i)*auxjno(ind)
     412            !NO2 interpolated coefficient
     413            jfotsout(indexint,13,auxi) = wm(i)*auxjno2(ind+1) +
     414     $           wp(i)*auxjno2(ind)
     415         enddo
     416      endif
     417c     End interval 16
     418
     419
     420ccccccccccccccccccccccccccccccc
     421c     90.9-119.5nm (int 17-24)
     422c
     423c     Absorption by:
     424c     CO2, O2, N2, NO, CO, NO2
     425ccccccccccccccccccccccccccccccc
     426
     427c     Input column
     428
     429      do i=1,nlayermx
     430         auxcolinp(nlayermx-i+1) = co2colx(i) + o2colx(i) + n2colx(i) +
     431     $        nocolx(i) + cocolx(i) + no2colx(i)
     432      end do
     433
     434      do indexint=17,24
     435
     436c     Interpolations
     437
    744438         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 
    755 c**************************************************
    756 c     no, 80.6-90.8nm
    757 c**************************************************
    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 
    782 c***********************************************************
    783 c     no2, 80.6-90.8nm
    784 c**************************************************
    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 
    810 c**************************************************
    811 c      co2, 90.9-119.5nm
    812 c**************************************************
    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)))
     439            auxi = nz2-i+1
     440            !CO2 tabulated coefficient
     441            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
     442            !O2 tabulated coefficient
     443            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
     444            !N2 tabulated coefficient
     445            auxjn2(i) = jabsifotsintpar(auxi,8,indexint)
     446            !CO tabulated coefficient
     447            auxjco(i) = jabsifotsintpar(auxi,11,indexint)           
     448            !Tabulated column
     449            auxcoltab(i) = c17_24(auxi)
     450         enddo
     451         !Only if chemthermod.ge.2
     452         if(chemthermod.ge.2) then
     453            do i=1,nz2
     454               auxi = nz2-i+1
     455               !NO tabulated coefficient
     456               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
     457               !NO2 tabulated coefficient
     458               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
     459            enddo
     460         endif
     461
     462         call interfast
     463     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
     464         !Correction to include T variation of CO2 cross section
     465         if(indexint.eq.24) then
     466            do i=1,nlayermx
     467               auxi = nlayermx-i+1
     468               if(sigma(indexint,auxi)*
     469     $              alfa(indexint,auxi)*coltemp(auxi)
     470     $              .lt.60.) then
     471                  cortemp(i)=exp(-sigma(indexint,auxi)*
     472     $                alfa(indexint,auxi)*coltemp(auxi))
    834473               else
    835                   jfotsout(indexint,1,i)=0.
    836                end if
    837             end if
    838          enddo
    839       end do
    840 
    841 
    842 c**************************************************
    843 c     o2, 90.9-119.5nm
    844 c**************************************************
    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 
    871 c**************************************************
    872 c     n2, 90.9-119.5nm
    873 c**************************************************
    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 
    900 c**************************************************
    901 c     co, 90.9-119.5nm
    902 c**************************************************
    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 
    932 c**************************************************
    933 c     no, 90.9-119.5nm
    934 c**************************************************
    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
     474                  cortemp(i)=0.
    957475               end if
    958476            enddo
    959          end do
    960 
    961 
    962 c**************************************************
    963 c     no2, 90.9-119.5nm
    964 c**************************************************
    965      
    966          do indexint=17,24
     477         else
    967478            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
     479               cortemp(i)=1.
     480            enddo
     481         end if
     482         do i=1,nlayermx           
     483            ind=auxind(i)
     484            auxi = nlayermx-i+1
     485            !O2 interpolated coefficient
     486            jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
     487     $           wp(i)*auxjo2(ind)) * cortemp(i)
     488            !CO2 interpolated coefficient
     489            jfotsout(indexint,1,auxi) = (wm(i)*auxjco2(ind+1) +
     490     $           wp(i)*auxjco2(ind)) * cortemp(i)
     491            if(indexint.eq.24) jfotsout(indexint,1,auxi)=
     492     $           jfotsout(indexint,1,auxi)*
     493     $           (1+alfa(indexint,auxi)*
     494     $           (t2(auxi)-t0(auxi)))
     495            !N2 interpolated coefficient
     496            jfotsout(indexint,8,auxi) = (wm(i)*auxjn2(ind+1) +
     497     $            wp(i)*auxjn2(ind)) * cortemp(i)           
     498            !CO interpolated coefficient
     499            jfotsout(indexint,11,auxi) = (wm(i)*auxjco(ind+1) +
     500     $            wp(i)*auxjco(ind)) * cortemp(i)           
     501         enddo
     502         !Only if chemthermod.ge.2
     503         if(chemthermod.ge.2) then
     504            do i=1,nlayermx
     505               ind=auxind(i)
     506               auxi = nlayermx-i+1
     507               !NO interpolated coefficient
     508               jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
     509     $              wp(i)*auxjno(ind)) * cortemp(i)
     510               !NO2 interpolated coefficient
     511               jfotsout(indexint,13,auxi)=(wm(i)*auxjno2(ind+1)+
     512     $              wp(i)*auxjno2(ind)) * cortemp(i)
     513            enddo
     514         endif               
     515      end do
     516c     End intervals 17-24
     517
     518
     519ccccccccccccccccccccccccccccccc
     520c     119.6-167.0nm (int 25-29)
     521c
     522c     Absorption by:
     523c     CO2, O2, H2O, H2O2, NO,
     524c     CO, NO2
     525ccccccccccccccccccccccccccccccc
     526
     527c     Input atmospheric column
     528
     529      do i=1,nlayermx
     530         auxcolinp(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
     531     $        h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
     532      end do
     533
     534      do indexint=25,29
     535
     536c     Interpolations
     537
     538         do i=1,nz2
     539            auxi = nz2-i+1
     540            !CO2 tabulated coefficient
     541            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
     542            !O2 tabulated coefficient
     543            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
     544            !H2O tabulated coefficient
     545            auxjh2o(i) = jabsifotsintpar(auxi,4,indexint)
     546            !H2O2 tabulated coefficient
     547            auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)           
     548            !CO tabulated coefficient
     549            auxjco(i) = jabsifotsintpar(auxi,11,indexint)           
     550            !Tabulated column
     551            auxcoltab(i) = c25_29(auxi)
     552         enddo
     553         !Only if chemthermod.ge.2
     554         if(chemthermod.ge.2) then
    971555            do i=1,nz2
    972                aux3(i) = jabsifotsintpar(indexint,13,nz2-i+1)
    973                aux4(i) = c17_24(nz2-i+1)
     556               auxi = nz2-i+1
     557               !NO tabulated coefficient
     558               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
     559               !NO2 tabulated coefficient
     560               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
    974561            enddo
    975             call interpfast
    976      $           (aux1,aux2,nlayermx,aux3,aux4,nz2,limdown,limup)
     562         endif
     563         call interfast
     564     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
     565         do i=1,nlayermx
     566            ind=auxind(i)
     567            auxi = nlayermx-i+1
     568            !Correction to include T variation of CO2 cross section
     569            if(sigma(indexint,auxi)*alfa(indexint,auxi)*
     570     $           coltemp(auxi).lt.60.) then
     571               cortemp(i)=exp(-sigma(indexint,auxi)*
     572     $              alfa(indexint,auxi)*coltemp(auxi))
     573            else
     574               cortemp(i)=0.
     575            end if
     576            !CO2 interpolated coefficient
     577            jfotsout(indexint,1,auxi) = (wm(i)*auxjco2(ind+1) +
     578     $           wp(i)*auxjco2(ind)) * cortemp(i) *
     579     $           (1+alfa(indexint,auxi)*
     580     $           (t2(auxi)-t0(auxi)))
     581            !O2 interpolated coefficient
     582            jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
     583     $           wp(i)*auxjo2(ind)) * cortemp(i)
     584            !H2O interpolated coefficient
     585            jfotsout(indexint,4,auxi) = (wm(i)*auxjh2o(ind+1) +
     586     $           wp(i)*auxjh2o(ind)) * cortemp(i)
     587            !H2O2 interpolated coefficient
     588            jfotsout(indexint,6,auxi) = (wm(i)*auxjh2o2(ind+1) +
     589     $           wp(i)*auxjh2o2(ind)) * cortemp(i)           
     590            !CO interpolated coefficient
     591            jfotsout(indexint,11,auxi) = (wm(i)*auxjco(ind+1) +
     592     $           wp(i)*auxjco(ind)) * cortemp(i)
     593         enddo
     594         !Only if chemthermod.ge.2
     595         if(chemthermod.ge.2) then
    977596            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
     597               ind=auxind(i)
     598               auxi = nlayermx-i+1
     599               !NO interpolated coefficient
     600               jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
     601     $              wp(i)*auxjno(ind)) * cortemp(i)
     602               !NO2 interpolated coefficient
     603               jfotsout(indexint,13,auxi)=(wm(i)*auxjno2(ind+1)+
     604     $              wp(i)*auxjno2(ind)) * cortemp(i)
    988605            enddo
    989          end do
    990 
    991          endif !Of chemthermod >= 2
    992 
    993 
    994 c**************************************************
    995 c      co2, 119.6-167.0nm
    996 c**************************************************
    997 
    998       do indexint=25,29
     606         endif
     607
     608      end do
     609
     610c     End intervals 25-29
     611
     612
     613cccccccccccccccccccccccccccccccc
     614c     167.1-202.5nm (int 30-31)
     615c   
     616c     Absorption by:
     617c     CO2, O2, H2O, H2O2, NO,
     618c     NO2
     619cccccccccccccccccccccccccccccccc
     620
     621c     Input atmospheric column
     622
     623      do i=1,nlayermx
     624         auxcolinp(nlayermx-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
     625     $        h2o2colx(i) + nocolx(i) + no2colx(i)
     626      end do
     627
     628c     Interpolation
     629
     630      do indexint=30,31
     631
     632         do i=1,nz2
     633            auxi = nz2-i+1
     634            !CO2 tabulated coefficient
     635            auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
     636            !O2 tabulated coefficient
     637            auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
     638            !H2O tabulated coefficient
     639            auxjh2o(i) = jabsifotsintpar(auxi,4,indexint)
     640            !H2O2 tabulated coefficient
     641            auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)           
     642            !Tabulated column
     643            auxcoltab(i) = c30_31(auxi)
     644         enddo
     645         !Only if chemthermod.ge.2
     646         if(chemthermod.ge.2) then
     647            do i=1,nz2
     648               auxi = nz2-i+1
     649               !NO tabulated coefficient
     650               auxjno(i) = jabsifotsintpar(auxi,10,indexint)
     651               !NO2 tabulated coefficient
     652               auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
     653            enddo
     654         endif
     655
     656         call interfast
     657     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
    999658         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
     659            ind=auxind(i)
     660            auxi = nlayermx-i+1
     661            !Correction to include T variation of CO2 cross section
     662            if(sigma(indexint,auxi)*alfa(indexint,auxi)*
     663     $           coltemp(auxi).lt.60.) then
     664               cortemp(i)=exp(-sigma(indexint,auxi)*
     665     $              alfa(indexint,auxi)*coltemp(auxi))
     666            else
     667               cortemp(i)=0.
     668            end if
     669            !CO2 interpolated coefficient
     670            jfotsout(indexint,1,auxi) = (wm(i)*auxjco2(ind+1) +
     671     $           wp(i)*auxjco2(ind)) * cortemp(i) *
     672     $           (1+alfa(indexint,auxi)*
     673     $           (t2(auxi)-t0(auxi)))
     674            !O2 interpolated coefficient
     675            jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
     676     $            wp(i)*auxjo2(ind)) * cortemp(i)
     677            !H2O interpolated coefficient
     678            jfotsout(indexint,4,auxi) = (wm(i)*auxjh2o(ind+1) +
     679     $            wp(i)*auxjh2o(ind)) * cortemp(i)
     680            !H2O2 interpolated coefficient
     681            jfotsout(indexint,6,auxi) = (wm(i)*auxjh2o2(ind+1) +
     682     $            wp(i)*auxjh2o2(ind)) * cortemp(i)           
     683         enddo
     684         !Only if chemthermod.ge.2
     685         if(chemthermod.ge.2) then
     686            do i=1,nlayermx
     687               ind=auxind(i)
     688               auxi = nlayermx-i+1
     689               !NO interpolated coefficient
     690               jfotsout(indexint,10,auxi)=(wm(i)*auxjno(ind+1) +
     691     $              wp(i)*auxjno(ind)) * cortemp(i)
     692               !NO2 interpolated coefficient
     693               jfotsout(indexint,13,auxi)=(wm(i)*auxjno2(ind+1)+
     694     $              wp(i)*auxjno2(ind)) * cortemp(i)
     695            enddo
     696         endif
     697
     698      end do
     699
     700c     End intervals 30-31
     701
     702
     703ccccccccccccccccccccccccccccccc
     704c     202.6-210.0nm (int 32)
     705c
     706c     Absorption by:
     707c     CO2, O2, H2O2, NO, NO2
     708ccccccccccccccccccccccccccccccc
     709
     710c     Input atmospheric column
     711
     712      indexint=32
     713      do i=1,nlayermx
     714         auxcolinp(nlayermx-i+1) =co2colx(i) + o2colx(i) + h2o2colx(i) +
     715     $        nocolx(i) + no2colx(i)
     716      end do
     717
     718c     Interpolation
     719
     720      do i=1,nz2
     721         auxi = nz2-i+1
     722         !CO2 tabulated coefficient
     723         auxjco2(i) = jabsifotsintpar(auxi,1,indexint)
     724         !O2 tabulated coefficient
     725         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
     726         !H2O2 tabulated coefficient
     727         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)         
     728         !Tabulated column
     729         auxcoltab(i) = c32(auxi)
     730      enddo
     731      !Only if chemthermod.ge.2
     732      if(chemthermod.ge.2) then
    1003733         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)
     734            auxi = nz2-i+1
     735            !NO tabulated coefficient
     736            auxjno(i) = jabsifotsintpar(auxi,10,indexint)
     737            !NO2 tabulated coefficient
     738            auxjno2(i) = jabsifotsintpar(auxi,13,indexint)
     739         enddo
     740      endif
     741      call interfast
     742     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
     743      do i=1,nlayermx
     744         ind=auxind(i)
     745         auxi = nlayermx-i+1
     746         !Correction to include T variation of CO2 cross section
     747         if(sigma(indexint,nlayermx-i+1)*alfa(indexint,auxi)*
     748     $        coltemp(auxi).lt.60.) then
     749            cortemp(i)=exp(-sigma(indexint,auxi)*
     750     $           alfa(indexint,auxi)*coltemp(auxi))
     751         else
     752            cortemp(i)=0.
     753         end if
     754         !CO2 interpolated coefficient
     755         jfotsout(indexint,1,auxi) = (wm(i)*auxjco2(ind+1) +
     756     $        wp(i)*auxjco2(ind)) * cortemp(i) *
     757     $        (1+alfa(indexint,auxi)*
     758     $        (t2(auxi)-t0(auxi)))
     759         !O2 interpolated coefficient
     760         jfotsout(indexint,2,auxi) = (wm(i)*auxjo2(ind+1) +
     761     $        wp(i)*auxjo2(ind)) * cortemp(i)
     762         !H2O2 interpolated coefficient
     763         jfotsout(indexint,6,auxi) = (wm(i)*auxjh2o2(ind+1) +
     764     $        wp(i)*auxjh2o2(ind)) * cortemp(i)         
     765      enddo
     766      !Only if chemthermod.ge.2
     767      if(chemthermod.ge.2) then
    1009768         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 
    1021 c**************************************************
    1022 c      o2, 119.6-167.0nm
    1023 c**************************************************
    1024 
    1025       do indexint=25,29
     769            auxi = nlayermx-i+1
     770            ind=auxind(i)
     771            !NO interpolated coefficient
     772            jfotsout(indexint,10,auxi) = (wm(i)*auxjno(ind+1) +
     773     $           wp(i)*auxjno(ind)) * cortemp(i)
     774           !NO2 interpolated coefficient
     775            jfotsout(indexint,13,auxi) = (wm(i)*auxjno2(ind+1) +
     776     $           wp(i)*auxjno2(ind)) * cortemp(i)
     777         enddo
     778      endif
     779
     780c     End of interval 32
     781
     782
     783ccccccccccccccccccccccccccccccc
     784c     210.1-231.0nm (int 33)
     785c     
     786c     Absorption by:
     787c     O2, H2O2, NO2
     788ccccccccccccccccccccccccccccccc
     789
     790c     Input atmospheric column
     791
     792      indexint=33
     793      do i=1,nlayermx
     794         auxcolinp(nlayermx-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
     795      end do
     796
     797c     Interpolation
     798
     799      do i=1,nz2
     800         auxi = nz2-i+1
     801         !O2 tabulated coefficient
     802         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
     803         !H2O2 tabulated coefficient
     804         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
     805         !Tabulated column
     806         auxcoltab(i) = c33(auxi)
     807      enddo
     808      !Only if chemthermod.ge.2
     809      if(chemthermod.ge.2) then
     810         do i=1,nz2
     811            !NO2 tabulated coefficient
     812            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
     813         enddo
     814      endif
     815      call interfast
     816     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
     817      do i=1,nlayermx
     818         ind=auxind(i)
     819         auxi = nlayermx-i+1
     820         !O2 interpolated coefficient
     821         jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
     822     $        wp(i)*auxjo2(ind)
     823         !H2O2 interpolated coefficient
     824         jfotsout(indexint,6,auxi) = wm(i)*auxjh2o2(ind+1) +
     825     $        wp(i)*auxjh2o2(ind)         
     826      enddo
     827      !Only if chemthermod.ge.2
     828      if(chemthermod.ge.2) then
    1026829         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
     830            ind=auxind(i)
     831            !NO2 interpolated coefficient
     832            jfotsout(indexint,13,nlayermx-i+1) = wm(i)*auxjno2(ind+1) +
     833     $           wp(i)*auxjno2(ind)
     834         enddo
     835      endif
     836
     837c     End of interval 33
     838
     839
     840ccccccccccccccccccccccccccccccc
     841c     231.1-240.0nm (int 34)
     842c
     843c     Absorption by:
     844c     O2, H2O2, O3, NO2
     845ccccccccccccccccccccccccccccccc
     846
     847c     Input atmospheric column
     848
     849      indexint=34
     850      do i=1,nlayermx
     851         auxcolinp(nlayermx-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) +
     852     $        no2colx(i)
     853      end do
     854
     855c     Interpolation
     856
     857      do i=1,nz2
     858         auxi = nz2-i+1
     859         !O2 tabulated coefficient
     860         auxjo2(i) = jabsifotsintpar(auxi,2,indexint)
     861         !H2O2 tabulated coefficient
     862         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
     863         !O3 tabulated coefficient
     864         auxjo3(i) = jabsifotsintpar(auxi,7,indexint)         
     865         !Tabulated column
     866         auxcoltab(i) = c34(nz2-i+1)
     867      enddo
     868      !Only if chemthermod.ge.2
     869      if(chemthermod.ge.2) then
    1030870         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)
     871            !NO2 tabulated coefficient
     872            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
     873         enddo
     874      endif
     875      call interfast
     876     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
     877      do i=1,nlayermx
     878         ind=auxind(i)
     879         auxi = nlayermx-i+1
     880         !O2 interpolated coefficient
     881         jfotsout(indexint,2,auxi) = wm(i)*auxjo2(ind+1) +
     882     $        wp(i)*auxjo2(ind)
     883         !H2O2 interpolated coefficient
     884         jfotsout(indexint,6,auxi) = wm(i)*auxjh2o2(ind+1) +
     885     $        wp(i)*auxjh2o2(ind)
     886         !O3 interpolated coefficient
     887         jfotsout(indexint,7,auxi) = wm(i)*auxjo3(ind+1) +
     888     $        wp(i)*auxjo3(ind)         
     889      enddo
     890      !Only if chemthermod.ge.2
     891      if(chemthermod.ge.2) then
    1036892         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 
    1047 c**************************************************
    1048 c      h2o, 119.6-167.0nm
    1049 c**************************************************
    1050 
    1051       do indexint=25,29
     893            ind=auxind(i)
     894            !NO2 interpolated coefficient
     895            jfotsout(indexint,13,nlayermx-i+1) = wm(i)*auxjno2(ind+1) +
     896     $           wp(i)*auxjno2(ind)
     897         enddo
     898      endif
     899
     900c     End of interval 34     
     901
     902
     903ccccccccccccccccccccccccccccccc
     904c     240.1-337.7nm (int 35)
     905c
     906c     Absorption by:
     907c     H2O2, O3, NO2
     908ccccccccccccccccccccccccccccccc
     909
     910c     Input atmospheric column
     911
     912      indexint=35
     913      do i=1,nlayermx
     914         auxcolinp(nlayermx-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
     915      end do
     916
     917c     Interpolation
     918
     919      do i=1,nz2
     920         auxi = nz2-i+1
     921         !H2O2 tabulated coefficient
     922         auxjh2o2(i) = jabsifotsintpar(auxi,6,indexint)
     923         !O3 tabulated coefficient
     924         auxjo3(i) = jabsifotsintpar(auxi,7,indexint)
     925         !Tabulated column
     926         auxcoltab(i) = c35(auxi)
     927      enddo
     928      !Only if chemthermod.ge.2
     929      if(chemthermod.ge.2) then
     930         do i=1,nz2
     931            !NO2 tabulated coefficient
     932            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
     933         enddo
     934      endif
     935      call interfast
     936     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
     937      do i=1,nlayermx
     938         ind=auxind(i)
     939         auxi = nlayermx-i+1
     940         !H2O2 interpolated coefficient
     941         jfotsout(indexint,6,auxi) = wm(i)*auxjh2o2(ind+1) +
     942     $        wp(i)*auxjh2o2(ind)
     943         !O3 interpolated coefficient
     944         jfotsout(indexint,7,auxi) = wm(i)*auxjo3(ind+1) +
     945     $        wp(i)*auxjo3(ind)         
     946      enddo
     947      if(chemthermod.ge.2) then
    1052948         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
     949            ind=auxind(i)
     950            !NO2 interpolated coefficient
     951            jfotsout(indexint,13,nlayermx-i+1) = wm(i)*auxjno2(ind+1) +
     952     $           wp(i)*auxjno2(ind)
     953         enddo
     954      endif
     955
     956c     End of interval 35
     957
     958ccccccccccccccccccccccccccccccc
     959c     337.8-800.0 nm (int 36)
     960c     
     961c     Absorption by:
     962c     O3, NO2
     963ccccccccccccccccccccccccccccccc
     964
     965c     Input atmospheric column
     966
     967      indexint=36
     968      do i=1,nlayermx
     969         auxcolinp(nlayermx-i+1) = o3colx(i) + no2colx(i)
     970      end do
     971
     972c     Interpolation
     973
     974      do i=1,nz2
     975         auxi = nz2-i+1
     976         !O3 tabulated coefficient
     977         auxjo3(i) = jabsifotsintpar(auxi,7,indexint)         
     978         !Tabulated column
     979         auxcoltab(i) = c36(auxi)
     980      enddo
     981      !Only if chemthermod.ge.2
     982      if(chemthermod.ge.2) then
    1056983         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)
     984            !NO2 tabulated coefficient
     985            auxjno2(i) = jabsifotsintpar(nz2-i+1,13,indexint)
     986         enddo
     987      endif
     988      call interfast
     989     $     (wm,wp,auxind,auxcolinp,nlayermx,auxcoltab,nz2,limdown,limup)
     990      do i=1,nlayermx
     991         ind=auxind(i)
     992         !O3 interpolated coefficient
     993         jfotsout(indexint,7,nlayermx-i+1) = wm(i)*auxjo3(ind+1) +
     994     $        wp(i)*auxjo3(ind)         
     995      enddo
     996      !Only if chemthermod.ge.2
     997      if(chemthermod.ge.2) then
    1062998         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 
    1073 c**************************************************
    1074 c      h2o2,119.6-167.0nm
    1075 c**************************************************
    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 
    1099 c**************************************************
    1100 c      co, 119.6-167.0nm
    1101 c**************************************************
    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 
    1128 c**************************************************
    1129 c      no, 119.6-167.0nm
    1130 c**************************************************
    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 
    1155 c**************************************************
    1156 c      no2, 119.6-167.0nm
    1157 c**************************************************
    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 
    1184 c**************************************************
    1185 c      co2, 167.0-202.5nm
    1186 c**************************************************
    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 
    1211 c**************************************************
    1212 c      o2, 167.0-202.5nm
    1213 c**************************************************
    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 
    1237 c**************************************************
    1238 c      h2o, 167.0-202.5nm
    1239 c**************************************************
    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 
    1263 c**************************************************
    1264 c      h2o2, 167.0-202.5nm
    1265 c**************************************************
    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      
    1292 c**************************************************
    1293 c      no, 167.0-202.5nm
    1294 c**************************************************
    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 
    1319 c**************************************************
    1320 c      no2, 167.0-202.5nm
    1321 c**************************************************
    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 
    1347 c**************************************************
    1348 c     co2, 202.6-210.0nm
    1349 c**************************************************
    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 
    1373 c**************************************************
    1374 c     o2, 202.6-210.0nm
    1375 c**************************************************
    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 
    1398 c**************************************************
    1399 c     h2o2, 202.6-210.0nm
    1400 c**************************************************
    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 
    1426 c**************************************************
    1427 c     no, 202.6-210.0nm
    1428 c**************************************************
    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 
    1452 c**************************************************
    1453 c     no2, 202.6-210.0nm
    1454 c**************************************************
    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 
    1480 c**************************************************
    1481 c     o2, 210.0-231.0nm
    1482 c**************************************************
    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 
    1499 c**************************************************
    1500 c     h2o2, 210.1-231.0nm
    1501 c**************************************************
    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 
    1523 c**************************************************
    1524 c     no2, 210.1-231.0nm
    1525 c**************************************************
    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 
    1546 c**************************************************
    1547 c     o2, 231.-240.nm
    1548 c**************************************************
    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 
    1568 c**************************************************
    1569 c     h2o2, 231.0-240.nm
    1570 c**************************************************
    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 
    1593 c**************************************************
    1594 c     o3, 231.0-240.nm
    1595 c**************************************************
    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 
    1620 c**************************************************
    1621 c     no2, 231.0-240.nm
    1622 c**************************************************
    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 
    1644 c**************************************************
    1645 c     h2o2, 240.0-337.7nm
    1646 c**************************************************
    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 
    1668 c**************************************************
    1669 c     o3, 240.0-337.7nm
    1670 c**************************************************
    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 
    1694 c**************************************************
    1695 c     no2, 240.0-337.7nm
    1696 c**************************************************
    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 
    1720 c**************************************************
    1721 c     o3, 337.7-800.0nm
    1722 c**************************************************
    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 
     999            ind=auxind(i)
     1000            !NO2 interpolated coefficient
     1001            jfotsout(indexint,13,nlayermx-i+1) = wm(i)*auxjno2(ind+1) +
     1002     $           wp(i)*auxjno2(ind)
     1003         enddo
    17401004      endif
    17411005
    1742 
    1743       !Only if Nitrogen or ionospheric chemistry requested
    1744       if(chemthermod.ge.2) then
    1745 
    1746 c**************************************************
    1747 c     no2, 337.7-800.0nm
    1748 c**************************************************
    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
     1006c     End of interval 36
     1007
     1008c     End of interpolation to obtain photoabsorption rates
     1009
     1010
    17671011      return
    17681012
    … …  
    18851129
    18861130      nz3 = nlayermx*2
    1887       szadeg = zenit*180./3.141592
    18881131      do i=1,nlayermx
    18891132         xx = ( radio + iz(i) ) * 1.e5
    … …  
    19111154         Hno2  = xx / mmol(igcm_no2)
    19121155      endif
     1156      ! first loop in altitude : initialisation
    19131157      do i=nlayermx,1,-1
    19141158         !Column initialisation
    … …  
    19451189            no2x(i)  = rm(i,i_no2)
    19461190         endif
     1191      enddo
     1192      ! second loop in altitude : column calculations
     1193      do i=nlayermx,1,-1
    19471194         !Routine to calculate the geometrical length of each layer
    19481195         call espesor_optico_A(ig,i,zenit,iz(i),nz3,iz,esp,ilayesp,
    … …  
    19711218               !Top layer
    19721219               if(jj.eq.nlayermx) then
    1973                   if(szadeg.le.60.) then
     1220                  if(zenit.le.60.) then
    19741221                     o3pcolx(i)=o3pcolx(i)+o3px(nlayermx)*Ho3p*esp(j)
    19751222     $                    *1.e-5
    … …  
    20031250     $                       *1.e-5
    20041251                     endif
    2005                   else if(szadeg.gt.60.) then
     1252                  else if(zenit.gt.60.) then
    20061253                     espco2 =sqrt((rcmnz+Hco2)**2 -rcmmini**2) - esp(j)
    20071254                     espo2  = sqrt((rcmnz+Ho2)**2 -rcmmini**2) - esp(j)
    … …  
    20401287                        no2colx(i) = no2colx(i) + espno2*no2x(nlayermx)
    20411288                     endif
    2042                   endif !Of if szadeg.lt.60
     1289                  endif !Of if zenit.lt.60
    20431290               !Other layers
    20441291               else
    … …  
    20861333c**********************************************************************
    20871334c**********************************************************************
    2088       subroutine interpfast(escout,p,nlayer,escin,pin,nl,limdown,limup)
     1335
     1336      subroutine interfast(wm,wp,nm,p,nlayer,pin,nl,limdown,limup)
    20891337C
    20901338C subroutine to perform linear interpolation in pressure from 1D profile
    … …  
    20921340C escout(nlayer) on pressure grid p(nlayer).
    20931341C
    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
     1342      real*8 wm(nlayer),wp(nlayer),p(nlayer)
     1343      integer nm(nlayer)
     1344      real*8 pin(nl)
     1345      real*8 limup,limdown
     1346      integer nl,nlayer,n1,n,np,nini
     1347      nini=1
     1348      do n1=1,nlayer
    21001349         if(p(n1) .gt. limup .or. p(n1) .lt. limdown) then
    2101             escout(n1) = 0.0
     1350            wm(n1) = 0.d0
     1351            wp(n1) = 0.d0
    21021352         else
    2103             do n = nm,nl-1
     1353            do n = nini,nl-1
    21041354               if (p(n1).ge.pin(n).and.p(n1).le.pin(n+1)) then
    2105                   nm=n
     1355                  nm(n1)=n
    21061356                  np=n+1
    2107                   wm=abs(pin(nm)-p(n1))/(pin(np)-pin(nm))
    2108                   wp=1.0 - wm
    2109                   goto 33
     1357                  wm(n1)=abs(pin(n)-p(n1))/(pin(np)-pin(n))
     1358                  wp(n1)=1.d0 - wm(n1)
     1359                  nini = n
     1360                  exit
    21101361               endif
    21111362            enddo
    2112  33         escout(n1) = escin(np)*wm + escin(nm)*wp
    21131363         endif
    2114     5 continue
     1364      enddo
    21151365      return
    21161366      end
    2117 
    21181367
    21191368

  • trunk/LMDZ.MARS/libf/aeronomars/moldiff.F

    r414 r658  
    279279c Inverting the alfa matrix
    280280c
    281         call ludcmp(alf,ncompmoldiff-1,ncompmoldiff-1,indx,d,ierr)
     281        call ludcmp_sp(alf,ncompmoldiff-1,ncompmoldiff-1,indx,d,ierr)
    282282
    283283c       TEMPORAIRE *****************************
    284284        if (ierr.ne.0) then
    285             write(*,*) 'In moldiff: Problem in LUDCMP with matrix alf'
     285            write(*,*)'In moldiff: Problem in LUDCMP_SP with matrix alf'
    286286            write(*,*) 'Singular matrix ?'
    287287c           write(*,*) 'Matrix alf = ', alf
    … …  
    294294c       *******************************************
    295295        do n=1,ncompmoldiff-1
    296           call lubksb(alf,ncompmoldiff-1,ncompmoldiff-1,indx,y(1,n))
     296          call lubksb_sp(alf,ncompmoldiff-1,ncompmoldiff-1,indx,y(1,n))
    297297          do nn=1,ncompmoldiff-1
    298298            alfinv(l,nn,n)=y(nn,n)/hh
    … …  
    415415        btri(1)=btri(1)+atri(1)
    416416
    417         call tridag(atri,btri,ctri,rtri,qtri,nz-2)
     417        call tridag_sp(atri,btri,ctri,rtri,qtri,nz-2)
    418418
    419419        do l=2,nz-1
    … …  
    464464c    ********************************************************************
    465465 
    466       subroutine tridag(a,b,c,r,u,n)
    467       parameter (nmax=100)
     466      subroutine tridag_sp(a,b,c,r,u,n)
     467c      parameter (nmax=100)
    468468c      dimension gam(nmax),a(n),b(n),c(n),r(n),u(n)
    469       real gam(nmax),a(n),b(n),c(n),r(n),u(n)
     469      real gam(n),a(n),b(n),c(n),r(n),u(n)
    470470      if(b(1).eq.0.)then
    471         stop 'tridag: error: b(1)=0 !!! '
     471        stop 'tridag_sp: error: b(1)=0 !!! '
    472472      endif
    473473      bet=b(1)
    … …  
    477477        bet=b(j)-a(j)*gam(j)
    478478        if(bet.eq.0.) then
    479           stop 'tridag: error: bet=0 !!! '
     479          stop 'tridag_sp: error: bet=0 !!! '
    480480        endif
    481481        u(j)=(r(j)-a(j)*u(j-1))/bet
    … …  
    491491c    ********************************************************************
    492492
    493       SUBROUTINE LUBKSB(A,N,NP,INDX,B)
     493      SUBROUTINE LUBKSB_SP(A,N,NP,INDX,B)
    494494
    495495      implicit none
    … …  
    530530c    ********************************************************************
    531531
    532       SUBROUTINE LUDCMP(A,N,NP,INDX,D,ierr)
     532      SUBROUTINE LUDCMP_SP(A,N,NP,INDX,D,ierr)
    533533
    534534      implicit none
    … …  
    55055011      CONTINUE
    551551        IF (AAMAX.EQ.0.) then
    552             write(*,*) 'In moldiff: Problem in LUDCMP with matrix A'
     552            write(*,*) 'In moldiff: Problem in LUDCMP_SP with matrix A'
    553553            write(*,*) 'Singular matrix ?'
    554554c           write(*,*) 'Matrix A = ', A

  • trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90

    r645 r658  
    726726
    727727      subroutine tridag(a,b,c,r,u,n)
    728       parameter (nmax=4000)
     728!      parameter (nmax=4000)
    729729!      dimension gam(nmax),a(n),b(n),c(n),r(n),u(n)
    730       real*8 gam(nmax),a(n),b(n),c(n),r(n),u(n)
     730      real*8 gam(n),a(n),b(n),c(n),r(n),u(n)
    731731      if(b(1).eq.0.)then
    732732        stop 'tridag: error: b(1)=0 !!! '

  • trunk/LMDZ.MARS/libf/aeronomars/param_read.F

    r635 r658  
    9595      IF (ierr.NE.0) THEN
    9696       write(*,*)'cant find directory EUVDAT and content coln.dat'
    97        write(*,*)'(in aeronomars/param_read_iono.F)'
     97       write(*,*)'(in aeronomars/param_read.F)'
    9898       write(*,*)'It should be in :', datafile,'/'
    9999       write(*,*)'1) You can change this directory address in '
    … …  
    181181
    182182      do i=nz2,1,-1
    183          read(220,*) (jabsifotsintpar(j,2,i),j=1,16)
    184       end do
    185      
    186       do i=nz2,1,-1
    187          read(230,*) (jabsifotsintpar(j,3,i),j=1,16)
    188       end do
    189 
    190       do i=nz2,1,-1
    191          read(240,*) (jabsifotsintpar(j,1,i),j=1,16)
    192       end do
    193 
    194       do i=nz2,1,-1
    195          read(250,*) (jabsifotsintpar(j,2,i),j=17,24)
    196       end do
    197 
    198 
    199       do i=nz2,1,-1
    200          read(260,*) (jabsifotsintpar(j,2,i),j=25,31)
    201       end do
    202 
    203       do i=nz2,1,-1
    204          read(270,*) (jabsifotsintpar(j,1,i),j=17,24)
    205       end do
    206 
    207       do i=nz2,1,-1
    208          read(280,*) (jabsifotsintpar(j,1,i),j=25,31)
    209       end do
    210 
    211       do i=nz2,1,-1
    212          read(290,*) jabsifotsintpar(32,1,i)
    213       end do
    214 
    215       do i=nz2,1,-1
    216          read(300,*) (jabsifotsintpar(j,2,i),j=32,34)
    217       end do
    218 
    219       do i=nz2,1,-1
    220          read(160,*) (jabsifotsintpar(j,5,i),j=1,15)
    221       end do
    222 
    223       do i=nz2,1,-1
    224          read(150,*) (jabsifotsintpar(j,4,i),j=25,31)
    225       end do
    226 
    227       do i=nz2,1,-1
    228          read(170,*) (jabsifotsintpar(j,6,i),j=25,35)
    229       end do
    230 
    231       do i=nz2,1,-1
    232          read(180,*) (jabsifotsintpar(j,7,i),j=34,36)
    233       end do
    234 
    235       do i=nz2,1,-1
    236          read(390,*) (jabsifotsintpar(j,8,i),j=2,16)
    237       enddo
    238 
    239       do i=nz2,1,-1
    240          read(400,*) (jabsifotsintpar(j,8,i),j=17,24)
    241       enddo
    242 
    243       do i=nz2,1,-1
    244          read(410,*) (jabsifotsintpar(j,9,i),j=1,16)
    245       enddo
    246 
    247       do i=nz2,1,-1
    248          read(420,*) (jabsifotsintpar(j,10,i),j=2,16)
    249       enddo
    250 
    251       do i=nz2,1,-1
    252          read(430,*) (jabsifotsintpar(j,10,i),j=17,24)
    253       enddo
    254 
    255       do i=nz2,1,-1
    256          read(440,*) (jabsifotsintpar(j,10,i),j=25,32)
    257       enddo
    258 
    259       do i=nz2,1,-1
    260          read(450,*) (jabsifotsintpar(j,11,i),j=2,16)
    261       enddo
    262 
    263       do i=nz2,1,-1
    264          read(460,*) (jabsifotsintpar(j,11,i),j=17,24)
    265       enddo
    266 
    267       do i=nz2,1,-1
    268          read(470,*) (jabsifotsintpar(j,11,i),j=25,29)
    269       enddo
    270      
    271       do i=nz2,1,-1
    272          read(480,*) (jabsifotsintpar(j,12,i),j=2,16)
    273       enddo
    274 
    275       do i=nz2,1,-1
    276          read(490,*) (jabsifotsintpar(j,13,i),j=2,16)
    277       enddo
    278      
    279       do i=nz2,1,-1
    280          read(500,*) (jabsifotsintpar(j,13,i),j=17,24)
    281       enddo
    282      
    283       do i=nz2,1,-1
    284          read(510,*) (jabsifotsintpar(j,13,i),j=25,36)
     183         read(220,*) (jabsifotsintpar(i,2,j),j=1,16)
     184      end do
     185     
     186      do i=nz2,1,-1
     187         read(230,*) (jabsifotsintpar(i,3,j),j=1,16)
     188      end do
     189
     190      do i=nz2,1,-1
     191         read(240,*) (jabsifotsintpar(i,1,j),j=1,16)
     192      end do
     193
     194      do i=nz2,1,-1
     195         read(250,*) (jabsifotsintpar(i,2,j),j=17,24)
     196      end do
     197
     198
     199      do i=nz2,1,-1
     200         read(260,*) (jabsifotsintpar(i,2,j),j=25,31)
     201      end do
     202
     203      do i=nz2,1,-1
     204         read(270,*) (jabsifotsintpar(i,1,j),j=17,24)
     205      end do
     206
     207      do i=nz2,1,-1
     208         read(280,*) (jabsifotsintpar(i,1,j),j=25,31)
     209      end do
     210
     211      do i=nz2,1,-1
     212         read(290,*) jabsifotsintpar(i,1,32)
     213      end do
     214
     215      do i=nz2,1,-1
     216         read(300,*) (jabsifotsintpar(i,2,j),j=32,34)
     217      end do
     218
     219      do i=nz2,1,-1
     220         read(160,*) (jabsifotsintpar(i,5,j),j=1,15)
     221      end do
     222
     223      do i=nz2,1,-1
     224         read(150,*) (jabsifotsintpar(i,4,j),j=25,31)
     225      end do
     226
     227      do i=nz2,1,-1
     228         read(170,*) (jabsifotsintpar(i,6,j),j=25,35)
     229      end do
     230
     231      do i=nz2,1,-1
     232         read(180,*) (jabsifotsintpar(i,7,j),j=34,36)
     233      end do
     234
     235      do i=nz2,1,-1
     236         read(390,*) (jabsifotsintpar(i,8,j),j=2,16)
     237      enddo
     238
     239      do i=nz2,1,-1
     240         read(400,*) (jabsifotsintpar(i,8,j),j=17,24)
     241      enddo
     242
     243      do i=nz2,1,-1
     244         read(410,*) (jabsifotsintpar(i,9,j),j=1,16)
     245      enddo
     246
     247      do i=nz2,1,-1
     248         read(420,*) (jabsifotsintpar(i,10,j),j=2,16)
     249      enddo
     250
     251      do i=nz2,1,-1
     252         read(430,*) (jabsifotsintpar(i,10,j),j=17,24)
     253      enddo
     254
     255      do i=nz2,1,-1
     256         read(440,*) (jabsifotsintpar(i,10,j),j=25,32)
     257      enddo
     258
     259      do i=nz2,1,-1
     260         read(450,*) (jabsifotsintpar(i,11,j),j=2,16)
     261      enddo
     262
     263      do i=nz2,1,-1
     264         read(460,*) (jabsifotsintpar(i,11,j),j=17,24)
     265      enddo
     266
     267      do i=nz2,1,-1
     268         read(470,*) (jabsifotsintpar(i,11,j),j=25,29)
     269      enddo
     270     
     271      do i=nz2,1,-1
     272         read(480,*) (jabsifotsintpar(i,12,j),j=2,16)
     273      enddo
     274
     275      do i=nz2,1,-1
     276         read(490,*) (jabsifotsintpar(i,13,j),j=2,16)
     277      enddo
     278     
     279      do i=nz2,1,-1
     280         read(500,*) (jabsifotsintpar(i,13,j),j=17,24)
     281      enddo
     282     
     283      do i=nz2,1,-1
     284         read(510,*) (jabsifotsintpar(i,13,j),j=25,36)
    285285      enddo
    286286

  • trunk/LMDZ.MARS/libf/aeronomars/param_v4.h

    r635 r658  
    2929      real fluxtop(ninter)
    3030      real freccen(ninter)          !representative wavelenght
    31       real jabsifotsintpar(ninter,nabs,nz2)
     31      real jabsifotsintpar(nz2,nabs,ninter)
    3232
    3333

  • trunk/LMDZ.MARS/libf/aeronomars/paramfoto_compact.F

    r635 r658  
    1 
    21c**********************************************************************
    32
    … …  
    14491448           enddo
    14501449           incremento = ( tehanson(i2)-tehanson(i1) ) /
    1451         1       ( zhanson(i2)-zhanson(i1) )
     1450     $          ( zhanson(i2)-zhanson(i1) )
    14521451           t_elect = tehanson(i1) + (zkm-zhanson(i1)) * incremento
    14531452        endif
    … …  
    18141813         !Initialization of lifetimes
    18151814         do j = 1, nreact
    1816             tauco2(i,j)  = 1.d30
    1817             tauo2(i,j)   = 1.d30
    1818             tauo3p(i,j)  = 1.d30
    1819             tauco(i,j)   = 1.d30
    1820             tauh(i,j)    = 1.d30
    1821             tauoh(i,j)   = 1.d30
    1822             tauho2(i,j)  = 1.d30
    1823             tauh2(i,j)   = 1.d30
    1824             tauh2o(i,j)  = 1.d30
    1825             tauo1d(i,j)  = 1.d30
    1826             tauh2o2(i,j) = 1.d30
    1827             tauo3(i,j)   = 1.d30
    1828             taun2(i,j)   = 1.d30
    1829             taun(i,j)    = 1.d30
    1830             tauno(i,j)   = 1.d30
    1831             taun2d(i,j)  = 1.d30
    1832             tauno2(i,j)  = 1.d30
    1833             tauco2plus(i,j) = 1.d30
    1834             tauoplus(i,j)   = 1.d30
    1835             tauo2plus(i,j)  = 1.d30
    1836             taucoplus(i,j)  = 1.d30
    1837             taucplus(i,j)   = 1.d30
    1838             taunplus(i,j)   = 1.d30
    1839             taun2plus(i,j)  = 1.d30
    1840             taunoplus(i,j)  = 1.d30
    1841             tauhplus(i,j)   = 1.d30
    1842             tauhco2plus(i,j)= 1.d30
     1815            tauco2(j,i)  = 1.d30
     1816            tauo2(j,i)   = 1.d30
     1817            tauo3p(j,i)  = 1.d30
     1818            tauco(j,i)   = 1.d30
     1819            tauh(j,i)    = 1.d30
     1820            tauoh(j,i)   = 1.d30
     1821            tauho2(j,i)  = 1.d30
     1822            tauh2(j,i)   = 1.d30
     1823            tauh2o(j,i)  = 1.d30
     1824            tauo1d(j,i)  = 1.d30
     1825            tauh2o2(j,i) = 1.d30
     1826            tauo3(j,i)   = 1.d30
     1827            taun2(j,i)   = 1.d30
     1828            taun(j,i)    = 1.d30
     1829            tauno(j,i)   = 1.d30
     1830            taun2d(j,i)  = 1.d30
     1831            tauno2(j,i)  = 1.d30
     1832            tauco2plus(j,i) = 1.d30
     1833            tauoplus(j,i)   = 1.d30
     1834            tauo2plus(j,i)  = 1.d30
     1835            taucoplus(j,i)  = 1.d30
     1836            taucplus(j,i)   = 1.d30
     1837            taunplus(j,i)   = 1.d30
     1838            taun2plus(j,i)  = 1.d30
     1839            taunoplus(j,i)  = 1.d30
     1840            tauhplus(j,i)   = 1.d30
     1841            tauhco2plus(j,i)= 1.d30
    18431842         end do
    18441843
    18451844         !Lifetime of each species in each reaction
    1846          if(jdistot8(1,i).gt.1.d-30) tauco2(i,1) = 1.d0 / jdistot8(1,i)
     1845         if(jdistot8(1,i).gt.1.d-30) tauco2(1,i) = 1.d0 / jdistot8(1,i)
    18471846                 
    18481847         if(ch2*o2xini*co2xini.gt.1.d-30)
    1849      $      tauh(i,2) = 1.d0 / (ch2 * o2xini * co2xini) 
     1848     $      tauh(2,i) = 1.d0 / (ch2 * o2xini * co2xini) 
    18501849         if(ch2*hxini*co2xini.gt.1.d-30)
    1851      $      tauo2(i,2) = 1.d0 / (ch2 * hxini * co2xini)
     1850     $      tauo2(2,i) = 1.d0 / (ch2 * hxini * co2xini)
    18521851 
    1853          if(ch3*o3pxini.gt.1.d-30) tauho2(i,3) = 1.d0 /
     1852         if(ch3*o3pxini.gt.1.d-30) tauho2(3,i) = 1.d0 /
    18541853     $        (ch3 * o3pxini)
    1855          if(ch3*ho2xini.gt.1.d-30) tauo3p(i,3) = 1.d0 /
     1854         if(ch3*ho2xini.gt.1.d-30) tauo3p(3,i) = 1.d0 /
    18561855     $        (ch3 * ho2xini)
    18571856 
    1858          if(ch4*coxini.gt.1.d-30) tauoh(i,4) = 1.d0 /
     1857         if(ch4*coxini.gt.1.d-30) tauoh(4,i) = 1.d0 /
    18591858     $        (ch4 * coxini)
    1860          if(ch4*ohxini.gt.1.d-30) tauco(i,4) = 1.d0 /
     1859         if(ch4*ohxini.gt.1.d-30) tauco(4,i) = 1.d0 /
    18611860     $        (ch4 * ohxini)
    18621861 
    1863          if(ch5*ho2xini.gt.1.d-30)tauho2(i,5)=1.d0 /
     1862         if(ch5*ho2xini.gt.1.d-30)tauho2(5,i)=1.d0 /
    18641863     $        (2.d0*ch5*ho2xini)
    18651864
    18661865                 
    1867          if(jdistot8(6,i).gt.1.d-30) tauh2o2(i,6) = 1.d0 / jdistot8(6,i)
    1868 
    1869          if(ch7*ohxini.gt.1.d-30) tauho2(i,7) = 1.d0 /
     1866         if(jdistot8(6,i).gt.1.d-30) tauh2o2(6,i) = 1.d0 / jdistot8(6,i)
     1867
     1868         if(ch7*ohxini.gt.1.d-30) tauho2(7,i) = 1.d0 /
    18701869     $        (ch7 * ohxini)
    1871          if(ch7*ho2xini.gt.1.d-30) tauoh(i,7) = 1.d0 /
     1870         if(ch7*ho2xini.gt.1.d-30) tauoh(7,i) = 1.d0 /
    18721871     $        (ch7 * ho2xini)
    18731872 
    1874          if(jdistot8(4,i).gt.1.d-30) tauh2o(i,8) = 1.d0 / jdistot8(4,i)
    1875 
    1876          if(ch9*o1dxini.gt.1.d-30) tauh2o(i,9) = 1.d0 /
     1873         if(jdistot8(4,i).gt.1.d-30) tauh2o(8,i) = 1.d0 / jdistot8(4,i)
     1874
     1875         if(ch9*o1dxini.gt.1.d-30) tauh2o(9,i) = 1.d0 /
    18771876     $        (ch9 * o1dxini)
    1878          if(ch9*h2oxini.gt.1.d-30) tauo1d(i,9) = 1.d0 /
     1877         if(ch9*h2oxini.gt.1.d-30) tauo1d(9,i) = 1.d0 /
    18791878     $        (ch9 * h2oxini)
    18801879
    18811880         if(ch10*o3pxini*co2xini.gt.1.d-30)
    1882      $     tauo3p(i,10) = 1.d0 /
     1881     $     tauo3p(10,i) = 1.d0 /
    18831882     $        (2.d0 * ch10 * o3pxini * co2xini)
    18841883
    1885          if(ch11*o3pxini.gt.1.d-30) tauoh(i,11)=1.d0 /
     1884         if(ch11*o3pxini.gt.1.d-30) tauoh(11,i)=1.d0 /
    18861885     $        (ch11 * o3pxini)
    1887          if(ch11*ohxini.gt.1.d-30) tauo3p(i,11) = 1.d0 /
     1886         if(ch11*ohxini.gt.1.d-30) tauo3p(11,i) = 1.d0 /
    18881887     $        (ch11 * ohxini)
    18891888
    1890          if(jdistot8(2,i).gt.1.d-30) tauo2(i,12) = 1.d0 / jdistot8(2,i)
    1891 
    1892          if(ch13*hxini.gt.1.d-30) tauho2(i,13) = 1.d0 /
     1889         if(jdistot8(2,i).gt.1.d-30) tauo2(12,i) = 1.d0 / jdistot8(2,i)
     1890
     1891         if(ch13*hxini.gt.1.d-30) tauho2(13,i) = 1.d0 /
    18931892     $        (ch13 * hxini)
    1894          if(ch13*ho2xini.gt.1.d-30) tauh(i,13) = 1.d0 /
     1893         if(ch13*ho2xini.gt.1.d-30) tauh(13,i) = 1.d0 /
    18951894     $        (ch13 * ho2xini)
    18961895
    1897          if(ch14*o1dxini.gt.1.d-30) tauh2(i,14) = 1.d0 /
     1896         if(ch14*o1dxini.gt.1.d-30) tauh2(14,i) = 1.d0 /
    18981897     $        (ch14 * o1dxini)
    1899          if(ch14*h2xini.gt.1.d-30) tauo1d(i,14) = 1.d0 /
     1898         if(ch14*h2xini.gt.1.d-30) tauo1d(14,i) = 1.d0 /
    19001899     $        (ch14 * h2xini)
    19011900
    1902          if(ch15*ohxini.gt.1.d-30) tauh2(i,15) = 1.d0 /
     1901         if(ch15*ohxini.gt.1.d-30) tauh2(15,i) = 1.d0 /
    19031902     $        (ch15 * ohxini)
    1904          if(ch15*h2xini.gt.1.d-30) tauoh(i,15) = 1.d0 /
     1903         if(ch15*h2xini.gt.1.d-30) tauoh(15,i) = 1.d0 /
    19051904     $        (ch15 * h2xini)
    19061905
    1907          if(jdistot8_b(1,i).gt.1.d-30) tauco2(i,16)=1.d0/jdistot8_b(1,i)
    1908 
    1909          if(jdistot8_b(2,i).gt.1.d-30) tauo2(i,17)=1.d0/jdistot8_b(2,i)
    1910 
    1911          if(ch18*ohxini.gt.1.d-30) tauh2o2(i,18)=1.d0 /
     1906         if(jdistot8_b(1,i).gt.1.d-30) tauco2(16,i)=1.d0/jdistot8_b(1,i)
     1907
     1908         if(jdistot8_b(2,i).gt.1.d-30) tauo2(17,i)=1.d0/jdistot8_b(2,i)
     1909
     1910         if(ch18*ohxini.gt.1.d-30) tauh2o2(18,i)=1.d0 /
    19121911     $        (ch18 * ohxini)
    1913          if(ch18*h2o2xini.gt.1.d-30) tauoh(i,18)=1.d0 /
     1912         if(ch18*h2o2xini.gt.1.d-30) tauoh(18,i)=1.d0 /
    19141913     $        (ch18 * h2o2xini)
    19151914
    1916          if(ch19*co2xini.gt.1.d-30)tauo1d(i,19)=1.d0 /
     1915         if(ch19*co2xini.gt.1.d-30)tauo1d(19,i)=1.d0 /
    19171916     $        (ch19 * co2xini)
    19181917
    1919          if(ch20*o2xini.gt.1.d-30)tauo1d(i,20)= 1.d0 /
     1918         if(ch20*o2xini.gt.1.d-30)tauo1d(20,i)= 1.d0 /
    19201919     $        (ch20 * o2xini)
    19211920
    1922          if(ch21*o2xini*co2xini.gt.1.d-30)tauo3p(i,21)= 1.d0 /
     1921         if(ch21*o2xini*co2xini.gt.1.d-30)tauo3p(21,i)= 1.d0 /
    19231922     $        (ch21 * o2xini * co2xini)
    1924          if(ch21*o3pxini*co2xini.gt.1.d-30) tauo2(i,21) = 1.d0 /
     1923         if(ch21*o3pxini*co2xini.gt.1.d-30) tauo2(21,i) = 1.d0 /
    19251924     $        (ch21 * o3pxini * co2xini)
    19261925
    19271926         !Only if O3, N or ion chemistry requested
    19281927         if(chemthermod.ge.1) then
    1929             if(ch22*hxini.gt.1.d-30) tauo3(i,22) = 1.d0 /
     1928            if(ch22*hxini.gt.1.d-30) tauo3(22,i) = 1.d0 /
    19301929     $           (ch22 * hxini)
    1931             if(ch22*o3xini.gt.1.d-30) tauh(i,22) = 1.d0 /
     1930            if(ch22*o3xini.gt.1.d-30) tauh(22,i) = 1.d0 /
    19321931     $           (ch22 * o3xini)
    19331932           
    1934             if(ch23*ohxini.gt.1.d-30) tauo3(i,23) = 1.d0 /
     1933            if(ch23*ohxini.gt.1.d-30) tauo3(23,i) = 1.d0 /
    19351934     $           (ch23 * ohxini)
    1936             if(ch23*o3xini.gt.1.d-30) tauoh(i,23) = 1.d0 /
     1935            if(ch23*o3xini.gt.1.d-30) tauoh(23,i) = 1.d0 /
    19371936     $           (ch23 * o3xini)
    19381937           
    1939             if(ch24 * ho2xini.gt.1.d-30)tauo3(i,24)= 1.d0 /
     1938            if(ch24 * ho2xini.gt.1.d-30)tauo3(24,i)= 1.d0 /
    19401939     $           (ch24 * ho2xini)
    1941             if(ch24 * o3xini.gt.1.d-30)tauho2(i,24)= 1.d0 /
     1940            if(ch24 * o3xini.gt.1.d-30)tauho2(24,i)= 1.d0 /
    19421941     $           (ch24 * o3xini)
    19431942           
    1944             if(jdistot8(7,i).gt.1.d-30) tauo3(i,25)=1.d0 /
     1943            if(jdistot8(7,i).gt.1.d-30) tauo3(25,i)=1.d0 /
    19451944     $           jdistot8(7,i)
    19461945
    1947             if(jdistot8_b(7,i).gt.1.d-30) tauo3(i,26)= 1.d0 /
     1946            if(jdistot8_b(7,i).gt.1.d-30) tauo3(26,i)= 1.d0 /
    19481947     $           jdistot8_b(7,i)
    19491948
    19501949         endif    !Of chemthermod.ge.1
    19511950
    1952          if(jdistot8(5,i).gt.1.d-30) tauh2(i,27)= 1.d0/jdistot8(5,i)
     1951         if(jdistot8(5,i).gt.1.d-30) tauh2(27,i)= 1.d0/jdistot8(5,i)
    19531952
    19541953!         if(ch28 * h2oxini.gt.1.d-30) tauo3(i,28) = 1.d0/(ch28*h2oxini)
    … …  
    19571956         !Only if N or ion chemistry requested
    19581957         if(chemthermod.ge.2) then
    1959             if(jdistot8(8,i).gt.1.d-30) taun2(i,28) = 1.d0 /
     1958            if(jdistot8(8,i).gt.1.d-30) taun2(28,i) = 1.d0 /
    19601959     $           jdistot8(8,i)
    19611960
    1962             if(jdistot8(10,i).gt.1.d-30) tauno(i,29) = 1.d0 /
     1961            if(jdistot8(10,i).gt.1.d-30) tauno(29,i) = 1.d0 /
    19631962     $           jdistot8(10,i)
    19641963
    1965             if(ch30 * noxini.gt.1.d-30) taun(i,30) = 1.d0 /
     1964            if(ch30 * noxini.gt.1.d-30) taun(30,i) = 1.d0 /
    19661965     $           (ch30 * noxini)
    1967             if(ch30 * nxini.gt.1.d-30) tauno(i,30) = 1.d0 /
     1966            if(ch30 * nxini.gt.1.d-30) tauno(30,i) = 1.d0 /
    19681967     $           (ch30 * nxini)
    19691968
    1970             if(ch31 * o1dxini.gt.1.d-30) taun2(i,31) = 1.d0 /
     1969            if(ch31 * o1dxini.gt.1.d-30) taun2(31,i) = 1.d0 /
    19711970     $           (ch31 * o1dxini)
    1972             if(ch31 * n2xini.gt.1.d-30) tauo1d(i,31) = 1.d0 /
     1971            if(ch31 * n2xini.gt.1.d-30) tauo1d(31,i) = 1.d0 /
    19731972     $           (ch31 * n2xini)
    19741973
    1975             if(ch32 * o2xini.gt.1.d-30) taun(i,32) = 1.d0 /
     1974            if(ch32 * o2xini.gt.1.d-30) taun(32,i) = 1.d0 /
    19761975     $           (ch32 * o2xini)
    1977             if(ch32 * nxini.gt.1.d-30) tauo2(i,32) = 1.d0 /
     1976            if(ch32 * nxini.gt.1.d-30) tauo2(32,i) = 1.d0 /
    19781977     $           (ch32 * nxini)
    19791978           
    1980             if(ch33 * ohxini.gt.1.d-30) taun(i,33) = 1.d0 /
     1979            if(ch33 * ohxini.gt.1.d-30) taun(33,i) = 1.d0 /
    19811980     $           (ch33 * ohxini)
    1982             if(ch33 * nxini.gt.1.d-30) tauoh(i,33) = 1.d0 /
     1981            if(ch33 * nxini.gt.1.d-30) tauoh(33,i) = 1.d0 /
    19831982     $           (ch33 * nxini)
    19841983
    1985             if(ch34 * o3xini.gt.1.d-30) taun(i,34) = 1.d0 /
     1984            if(ch34 * o3xini.gt.1.d-30) taun(34,i) = 1.d0 /
    19861985     $           (ch34 * o3xini)
    1987             if(ch34 * nxini.gt.1.d-30) tauo3(i,34) = 1.d0 /
     1986            if(ch34 * nxini.gt.1.d-30) tauo3(34,i) = 1.d0 /
    19881987     $           (ch34 * nxini)
    19891988           
    1990             if(ch35 * ho2xini.gt.1.d-30) taun(i,35) = 1.d0 /
     1989            if(ch35 * ho2xini.gt.1.d-30) taun(35,i) = 1.d0 /
    19911990     $           (ch35 * ho2xini)
    1992             if(ch35 * nxini.gt.1.d-30) tauho2(i,35) = 1.d0 /
     1991            if(ch35 * nxini.gt.1.d-30) tauho2(35,i) = 1.d0 /
    19931992     $           (ch35 * nxini)
    19941993           
    1995             if(ch36 * o3pxini.gt.1.d-30) taun2d(i,36) = 1.d0 /
     1994            if(ch36 * o3pxini.gt.1.d-30) taun2d(36,i) = 1.d0 /
    19961995     $           (ch36 * o3pxini)
    1997             if(ch36 * n2dxini.gt.1.d-30) tauo3p(i,36) = 1.d0 /
     1996            if(ch36 * n2dxini.gt.1.d-30) tauo3p(36,i) = 1.d0 /
    19981997     $           (ch36 * n2dxini)
    19991998           
    2000             if(ch37 * n2xini.gt.1.d-30) taun2d(i,37) = 1.d0 /
     1999            if(ch37 * n2xini.gt.1.d-30) taun2d(37,i) = 1.d0 /
    20012000     $           (ch37 * n2xini)
    2002             if(ch37 * n2dxini.gt.1.d-30) taun2(i,37) = 1.d0 /
     2001            if(ch37 * n2dxini.gt.1.d-30) taun2(37,i) = 1.d0 /
    20032002     $           (ch37 * n2dxini)
    20042003           
    2005             if(ch38 * co2xini.gt.1.d-30) taun2d(i,38) = 1.d0 /
     2004            if(ch38 * co2xini.gt.1.d-30) taun2d(38,i) = 1.d0 /
    20062005     $           (ch38 * co2xini)
    2007             if(ch38 * n2dxini.gt.1.d-30) tauco2(i,38) = 1.d0 /
     2006            if(ch38 * n2dxini.gt.1.d-30) tauco2(38,i) = 1.d0 /
    20082007     $           (ch38 * n2dxini)
    20092008           
    2010             if(ch39 * ho2xini.gt.1.d-30) tauno(i,39) = 1.d0 /
     2009            if(ch39 * ho2xini.gt.1.d-30) tauno(39,i) = 1.d0 /
    20112010     $           (ch39 * ho2xini)
    2012             if(ch39 * noxini.gt.1.d-30) tauho2(i,39) = 1.d0 /
     2011            if(ch39 * noxini.gt.1.d-30) tauho2(39,i) = 1.d0 /
    20132012     $           (ch39 * noxini)
    20142013           
    2015             if(ch40 * noxini * co2xini.gt.1.d-30) tauo3p(i,40) = 1.d0 /
     2014            if(ch40 * noxini * co2xini.gt.1.d-30) tauo3p(40,i) = 1.d0 /
    20162015     $           (ch40 * noxini * co2xini)
    2017             if(ch40 * o3pxini * co2xini.gt.1.d-30) tauno(i,40) = 1.d0 /
     2016            if(ch40 * o3pxini * co2xini.gt.1.d-30) tauno(40,i) = 1.d0 /
    20182017     $           (ch40 * o3pxini * co2xini)
    20192018           
    2020             if(ch41 * no2xini.gt.1.d-30) tauo3p(i,41) = 1.d0 /
     2019            if(ch41 * no2xini.gt.1.d-30) tauo3p(41,i) = 1.d0 /
    20212020     $           (ch41 * no2xini)
    2022             if(ch41 * o3pxini.gt.1.d-30) tauno2(i,41) = 1.d0 /
     2021            if(ch41 * o3pxini.gt.1.d-30) tauno2(41,i) = 1.d0 /
    20232022     $           (ch41 * o3pxini)
    20242023           
    2025             if(ch42 * noxini.gt.1.d-30) tauo3(i,42) = 1.d0 /
     2024            if(ch42 * noxini.gt.1.d-30) tauo3(42,i) = 1.d0 /
    20262025     $           (ch42 * noxini)
    2027             if(ch42 * o3xini.gt.1.d-30) tauno(i,42) = 1.d0 /
     2026            if(ch42 * o3xini.gt.1.d-30) tauno(42,i) = 1.d0 /
    20282027     $           (ch42 * o3xini)
    20292028           
    2030             if(ch43 * no2xini.gt.1.d-30) tauh(i,43) = 1.d0 /
     2029            if(ch43 * no2xini.gt.1.d-30) tauh(43,i) = 1.d0 /
    20312030     $           (ch43 * no2xini)
    2032             if(ch43 * hxini.gt.1.d-30) tauno2(i,43) = 1.d0 /
     2031            if(ch43 * hxini.gt.1.d-30) tauno2(43,i) = 1.d0 /
    20332032     $           (ch43 * hxini)
    20342033
    2035             if(jdistot8(13,i).gt.1.d-30) tauno2(i,44) = 1.d0 /
     2034            if(jdistot8(13,i).gt.1.d-30) tauno2(44,i) = 1.d0 /
    20362035     $           jdistot8(13,i)
    20372036
    2038             if(ch45 * nxini.gt.1.d-30) tauo3p(i,45) = 1.d0 /
     2037            if(ch45 * nxini.gt.1.d-30) tauo3p(45,i) = 1.d0 /
    20392038     $           (ch45 * nxini)
    2040             if(ch45 * o3pxini.gt.1.d-30) taun(i,45) = 1.d0 /
     2039            if(ch45 * o3pxini.gt.1.d-30) taun(45,i) = 1.d0 /
    20412040     $           (ch45 * o3pxini)
    20422041
    … …  
    20472046         !Only if ion chemistry requested
    20482047         if(chemthermod.eq.3) then
    2049             if(ch46 * co2plusxini .gt.1.d-30) tauo2(i,46) =
     2048            if(ch46 * co2plusxini .gt.1.d-30) tauo2(46,i) =
    20502049     @           1.d0/(ch46*co2plusxini)
    2051             if(ch46 * o2xini .gt.1.d-30) tauco2plus(i,46) =
     2050            if(ch46 * o2xini .gt.1.d-30) tauco2plus(46,i) =
    20522051     @           1.d0/(ch46*o2xini)
    20532052
    2054             if ( ch47*o3pxini .gt. 1.d-30 ) tauco2plus(i,47) =
     2053            if ( ch47*o3pxini .gt. 1.d-30 ) tauco2plus(47,i) =
    20552054     @           1.d0/( ch47*o3pxini )
    2056             if ( ch47*co2plusxini .gt. 1.d-30 ) tauo3p(i,47) =
     2055            if ( ch47*co2plusxini .gt. 1.d-30 ) tauo3p(47,i) =
    20572056     @           1.d0/( ch47*co2plusxini )
    20582057           
    2059             if ( ch48*o3pxini .gt. 1.d-30 ) tauco2plus(i,48) =
     2058            if ( ch48*o3pxini .gt. 1.d-30 ) tauco2plus(48,i) =
    20602059     @           1.d0/(ch48*o3pxini)       
    2061             if ( ch48*co2plusxini .gt. 1.d-30 ) tauo3p(i,48) =
     2060            if ( ch48*co2plusxini .gt. 1.d-30 ) tauo3p(48,i) =
    20622061     @           1.d0/(ch48*co2plusxini)         
    20632062
    2064             if ( ch49*electxini .gt. 1.d-30 ) tauo2plus(i,49) =
     2063            if ( ch49*electxini .gt. 1.d-30 ) tauo2plus(49,i) =
    20652064     @           1.d0/(ch49*electxini)   
    20662065
    2067             if ( ch50*co2xini  .gt. 1.d-30 ) tauoplus(i,50) =
     2066            if ( ch50*co2xini  .gt. 1.d-30 ) tauoplus(50,i) =
    20682067     @           1.d0/(ch50*co2xini)
    2069             if ( ch50*oplusxini .gt. 1.d-30 ) tauco2(i,50) =
     2068            if ( ch50*oplusxini .gt. 1.d-30 ) tauco2(50,i) =
    20702069     @           1.d0/(ch50*oplusxini)
    20712070
    2072             if ( jion8(1,i,1).gt.1.d-30 ) tauco2(i,51) =
     2071            if ( jion8(1,i,1).gt.1.d-30 ) tauco2(51,i) =
    20732072     $           1.d0 / jion8(1,i,1)
    20742073
    2075             if ( jion8(1,i,2).gt.1.d-30 ) tauco2(i,52) =
     2074            if ( jion8(1,i,2).gt.1.d-30 ) tauco2(52,i) =
    20762075     $           1.d0 / jion8(1,i,2)
    20772076
    2078             if ( jion8(1,i,3).gt.1.d-30 ) tauco2(i,53) =
     2077            if ( jion8(1,i,3).gt.1.d-30 ) tauco2(53,i) =
    20792078     $           1.d0 / jion8(1,i,3)
    20802079
    2081             if ( jion8(1,i,4).gt.1.d-30 ) tauco2(i,54) =
     2080            if ( jion8(1,i,4).gt.1.d-30 ) tauco2(54,i) =
    20822081     $           1.d0 / jion8(1,i,4)
    20832082               
    2084             if ( ch55*electxini .gt. 1.d-30 ) tauco2plus(i,55) =
     2083            if ( ch55*electxini .gt. 1.d-30 ) tauco2plus(55,i) =
    20852084     @           1.d0/(ch55*electxini)   
    20862085
    2087             if ( ch56*oplusxini .gt. 1.d-30 ) tauco2(i,56) =
     2086            if ( ch56*oplusxini .gt. 1.d-30 ) tauco2(56,i) =
    20882087     @           1.d0/(ch56*oplusxini)
    2089             if ( ch56*co2xini .gt. 1.d-30 ) tauoplus(i,56) =
     2088            if ( ch56*co2xini .gt. 1.d-30 ) tauoplus(56,i) =
    20902089     @           1.d0/(ch56*co2xini)   
    20912090   
    2092             if ( ch57*coplusxini .gt. 1.d-30 ) tauco2(i,57) =
     2091            if ( ch57*coplusxini .gt. 1.d-30 ) tauco2(57,i) =
    20932092     @           1.d0/(ch57*coplusxini) 
    2094             if ( ch57*co2xini .gt. 1.d-30 ) taucoplus(i,57) =
     2093            if ( ch57*co2xini .gt. 1.d-30 ) taucoplus(57,i) =
    20952094     @           1.d0/(ch57*co2xini)
    20962095
    2097             if ( ch58*coplusxini .gt. 1.d-30 ) tauo3p(i,58) =
     2096            if ( ch58*coplusxini .gt. 1.d-30 ) tauo3p(58,i) =
    20982097     @           1.d0/(ch58*coplusxini)
    2099             if ( ch58*o3pxini .gt. 1.d-30 ) taucoplus(i,58) =
     2098            if ( ch58*o3pxini .gt. 1.d-30 ) taucoplus(58,i) =
    21002099     @           1.d0/(ch58*o3pxini)
    21012100
    2102             if ( ch59*cplusxini .gt. 1.d-30 ) tauco2(i,59) =
     2101            if ( ch59*cplusxini .gt. 1.d-30 ) tauco2(59,i) =
    21032102     @           1.d0/(ch59*cplusxini) 
    2104             if ( ch59*co2xini .gt. 1.d-30 ) taucplus(i,59) =
     2103            if ( ch59*co2xini .gt. 1.d-30 ) taucplus(59,i) =
    21052104     @           1.d0/(ch59*co2xini)                     
    21062105         
    2107             if ( jion8(2,i,1).gt.1.d-30 ) tauo2(i,60) =
     2106            if ( jion8(2,i,1).gt.1.d-30 ) tauo2(60,i) =
    21082107     $           1.d0 / jion8(2,i,1)
    21092108
    2110             if ( jion8(3,i,1).gt.1.d-30 ) tauo3p(i,61) =
     2109            if ( jion8(3,i,1).gt.1.d-30 ) tauo3p(61,i) =
    21112110     $           1.d0 / jion8(3,i,1) 
    21122111
    2113             if ( ch62*co2plusxini .gt. 1.d-30 ) tauno(i,62) =
     2112            if ( ch62*co2plusxini .gt. 1.d-30 ) tauno(62,i) =
    21142113     @           1.d0/(ch62*co2plusxini)       
    2115             if ( ch62*noxini .gt. 1.d-30 ) tauco2plus(i,62) =
     2114            if ( ch62*noxini .gt. 1.d-30 ) tauco2plus(62,i) =
    21162115     @           1.d0/(ch62*noxini)                     
    21172116 
    2118             if ( ch63*co2plusxini .gt. 1.d-30 ) taun(i,63) =
     2117            if ( ch63*co2plusxini .gt. 1.d-30 ) taun(63,i) =
    21192118     @           1.d0/(ch63*cplusxini) 
    2120             if ( ch63*nxini .gt. 1.d-30 ) tauco2plus(i,63) =
     2119            if ( ch63*nxini .gt. 1.d-30 ) tauco2plus(63,i) =
    21212120     @           1.d0/(ch63*nxini)                       
    21222121     
    2123             if ( ch64*o2plusxini .gt. 1.d-30 ) tauno(i,64) =
     2122            if ( ch64*o2plusxini .gt. 1.d-30 ) tauno(64,i) =
    21242123     @           1.d0/(ch64*o2plusxini)
    2125             if ( ch64*noxini .gt. 1.d-30 ) tauo2plus(i,64) =
     2124            if ( ch64*noxini .gt. 1.d-30 ) tauo2plus(64,i) =
    21262125     @           1.d0/(ch64*noxini)                     
    21272126   
    2128             if ( ch65*o2plusxini .gt. 1.d-30 ) taun2(i,65) =
     2127            if ( ch65*o2plusxini .gt. 1.d-30 ) taun2(65,i) =
    21292128     @           1.d0/(ch65*o2plusxini)
    2130             if ( ch65*n2xini .gt. 1.d-30 ) tauo2plus(i,65) =
     2129            if ( ch65*n2xini .gt. 1.d-30 ) tauo2plus(65,i) =
    21312130     @           1.d0/(ch65*n2xini)                     
    21322131   
    2133             if ( ch66*o2plusxini .gt. 1.d-30 ) taun(i,66) =
     2132            if ( ch66*o2plusxini .gt. 1.d-30 ) taun(66,i) =
    21342133     @           1.d0/(ch66*o2plusxini)
    2135             if ( ch66*nxini .gt. 1.d-30 ) tauo2plus(i,66) =
     2134            if ( ch66*nxini .gt. 1.d-30 ) tauo2plus(66,i) =
    21362135     @           1.d0/(ch66*nxini)                       
    21372136   
    2138             if ( ch67*oplusxini .gt. 1.d-30 ) taun2(i,67) =
     2137            if ( ch67*oplusxini .gt. 1.d-30 ) taun2(67,i) =
    21392138     @           1.d0/(ch67*oplusxini) 
    2140             if ( ch67*n2xini .gt. 1.d-30 ) tauoplus(i,67) =
     2139            if ( ch67*n2xini .gt. 1.d-30 ) tauoplus(67,i) =
    21412140     @           1.d0/(ch67*n2xini)                     
    21422141
    2143             if ( ch68*n2plusxini .gt. 1.d-30 ) tauco2(i,68) =
     2142            if ( ch68*n2plusxini .gt. 1.d-30 ) tauco2(68,i) =
    21442143     @           1.d0/(ch68*n2plusxini)
    2145             if ( ch68*co2xini .gt. 1.d-30 ) taun2plus(i,68) =
     2144            if ( ch68*co2xini .gt. 1.d-30 ) taun2plus(68,i) =
    21462145     @           1.d0/(ch68*co2xini)                     
    21472146   
    2148             if ( ch69*cplusxini .gt. 1.d-30 ) tauco2(i,69) =
     2147            if ( ch69*cplusxini .gt. 1.d-30 ) tauco2(69,i) =
    21492148     @           1.d0/(ch69*cplusxini) 
    2150             if ( ch69*co2xini .gt. 1.d-30 ) taucplus(i,69) =
     2149            if ( ch69*co2xini .gt. 1.d-30 ) taucplus(69,i) =
    21512150     @           1.d0/(ch69*co2xini)                     
    21522151           
    2153             if ( ch70*n2plusxini .gt. 1.d-30 ) tauco(i,70) =
     2152            if ( ch70*n2plusxini .gt. 1.d-30 ) tauco(70,i) =
    21542153     @           1.d0/(ch70*n2plusxini)
    2155             if ( ch70*coxini .gt. 1.d-30 ) taun2plus(i,70) =
     2154            if ( ch70*coxini .gt. 1.d-30 ) taun2plus(70,i) =
    21562155     @           1.d0/(ch70*coxini)                     
    21572156           
    2158             if ( ch71*electxini .gt. 1.d-30 ) taun2plus(i,71) =
     2157            if ( ch71*electxini .gt. 1.d-30 ) taun2plus(71,i) =
    21592158     @           1.d0/(ch71*electxini) 
    21602159           
    2161             if ( ch72*n2plusxini .gt. 1.d-30 ) tauo3p(i,72) =
     2160            if ( ch72*n2plusxini .gt. 1.d-30 ) tauo3p(72,i) =
    21622161     @           1.d0/(ch72*n2plusxini)
    2163             if ( ch72*o3pxini .gt. 1.d-30 ) taun2plus(i,72) =
     2162            if ( ch72*o3pxini .gt. 1.d-30 ) taun2plus(72,i) =
    21642163     @           1.d0/(ch72*o3pxini)                     
    21652164                 
    2166             if ( ch73*coplusxini .gt. 1.d-30 ) tauh(i,73) =
     2165            if ( ch73*coplusxini .gt. 1.d-30 ) tauh(73,i) =
    21672166     @           1.d0/(ch73*coplusxini)
    2168             if ( ch73*hxini .gt. 1.d-30 ) taucoplus(i,73) =
     2167            if ( ch73*hxini .gt. 1.d-30 ) taucoplus(73,i) =
    21692168     @           1.d0/(ch73*hxini)                       
    21702169         
    2171             if ( ch74*oplusxini .gt. 1.d-30 ) tauh(i,74) =
     2170            if ( ch74*oplusxini .gt. 1.d-30 ) tauh(74,i) =
    21722171     @           1.d0/(ch74*oplusxini) 
    2173             if ( ch74*hxini .gt. 1.d-30 ) tauoplus(i,74) =
     2172            if ( ch74*hxini .gt. 1.d-30 ) tauoplus(74,i) =
    21742173     @           1.d0/(ch74*hxini)     
    21752174 
    2176             if ( ch75*electxini .gt. 1.d-30 ) taunoplus(i,75) =
     2175            if ( ch75*electxini .gt. 1.d-30 ) taunoplus(75,i) =
    21772176     @           1.d0/(ch75*electxini)                   
    21782177           
    2179             if ( ch76*hplusxini .gt. 1.d-30 ) tauo3p(i,76) =
     2178            if ( ch76*hplusxini .gt. 1.d-30 ) tauo3p(76,i) =
    21802179     @           1.d0/(ch76*hplusxini) 
    2181             if ( ch76*o3pxini .gt. 1.d-30 ) tauhplus(i,76) =
     2180            if ( ch76*o3pxini .gt. 1.d-30 ) tauhplus(76,i) =
    21822181     @           1.d0/(ch76*o3pxini)   
    21832182         
    2184             if( jion8(11,i,1).gt.1.d-30 )  tauco(i,77) =
     2183            if( jion8(11,i,1).gt.1.d-30 )  tauco(77,i) =
    21852184     $           1.d0 / jion8(11,i,1)
    21862185
    2187             if( jion8(11,i,2).gt.1.d-30 )  tauco(i,78) =
     2186            if( jion8(11,i,2).gt.1.d-30 )  tauco(78,i) =
    21882187     $           1.d0 / jion8(11,i,2) 
    21892188
    2190          !if( jion8(11,i,3).gt.1.d-30 ) tauco(i,79) =
     2189         !if( jion8(11,i,3).gt.1.d-30 ) tauco(79,i) =
    21912190!     $        1.d0 / jion8(11,i,3)
    21922191
    2193             if( jion8(10,i,1).gt.1.d-30 )  tauno(i,80) =
     2192            if( jion8(10,i,1).gt.1.d-30 )  tauno(80,i) =
    21942193     $           1.d0 / jion8(10,i,1)
    21952194         
    2196             if( jion8(8,i,1).gt.1.d-30 )   taun2(i,81)  =
     2195            if( jion8(8,i,1).gt.1.d-30 )   taun2(81,i)  =
    21972196     $           1.d0 / jion8(8,i,1)
    21982197
    2199             if( jion8(8,i,2).gt.1.d-30 )   taun2(i,82)  =
     2198            if( jion8(8,i,2).gt.1.d-30 )   taun2(82,i)  =
    22002199     $           1.d0 / jion8(8,i,2)
    22012200                 
    2202             if( jion8(12,i,1).gt.1.d-30 )  tauh(i,83)  =
     2201            if( jion8(12,i,1).gt.1.d-30 )  tauh(83,i)  =
    22032202     $           1.d0 / jion8(12,i,1)
    22042203
    2205             if( jion8(9,i,1).gt.1.d-30 )   taun(i,84)  =
     2204            if( jion8(9,i,1).gt.1.d-30 )   taun(84,i)  =
    22062205     $           1.d0 / jion8(9,i,1)
    22072206
    2208             if ( ch85*nplusxini .gt. 1.d-30 ) tauco2(i,85) =
     2207            if ( ch85*nplusxini .gt. 1.d-30 ) tauco2(85,i) =
    22092208     @           1.d0/(ch85*nplusxini) 
    2210             if ( ch85*co2xini .gt. 1.d-30 ) taunplus(i,85) =
     2209            if ( ch85*co2xini .gt. 1.d-30 ) taunplus(85,i) =
    22112210     @           1.d0/(ch85*co2xini)
    22122211
    2213             if ( ch86*co2plusxini .gt. 1.d-30) tauh2(i,86) =
     2212            if ( ch86*co2plusxini .gt. 1.d-30) tauh2(86,i) =
    22142213     $           1.d0/(ch86*co2plusxini)
    2215             if ( ch86*h2xini .gt. 1.d-30) tauco2plus(i,86) =
     2214            if ( ch86*h2xini .gt. 1.d-30) tauco2plus(86,i) =
    22162215     $           1.d0/(ch86*h2xini)
    22172216         
    2218             if ( ch87*electxini .gt. 1.d-30) tauhco2plus(i,87) =
     2217            if ( ch87*electxini .gt. 1.d-30) tauhco2plus(87,i) =
    22192218     $                  1.d0/(ch87*electxini)
    22202219
    22212220            if ( jion8(9,i,1)*ionsec_nplus(zenit,zx(i)).gt.1.d-30)
    2222      $           taun(i,88) = 1.d0 /
     2221     $           taun(88,i) = 1.d0 /
    22232222     $           (jion8(9,i,1)*ionsec_nplus(zenit,zx(i)))
    22242223
    22252224            if ( jion8(8,i,1)*ionsec_n2plus(zenit,zx(i)).gt.1.d-30)
    2226      $           taun2(i,89) = 1.d0 /
     2225     $           taun2(89,i) = 1.d0 /
    22272226     $           (jion8(8,i,1)*ionsec_n2plus(zenit,zx(i)))
    22282227           
    22292228            if ( jion8(3,i,1)*ionsec_oplus(zenit,zx(i)).gt.1.d-30)
    2230      $           tauo3p(i,90) = 1.d0 /
     2229     $           tauo3p(90,i) = 1.d0 /
    22312230     $           (jion8(3,i,1)*ionsec_oplus(zenit,zx(i)))
    22322231           
    22332232            if (jion8(11,i,1)*ionsec_coplus(zenit,zx(i)).gt.1.d-30)
    2234      $           tauco(i,91) = 1.d0 /
     2233     $           tauco(91,i) = 1.d0 /
    22352234     $           (jion8(11,i,1)*ionsec_coplus(zenit,zx(i)))
    22362235
    22372236            if (jion8(1,i,1)*ionsec_co2plus(zenit,zx(i)).gt.1.d-30)
    2238      $           tauco2(i,92) = 1.d0 /
     2237     $           tauco2(92,i) = 1.d0 /
    22392238     $           (jion8(1,i,1)*ionsec_co2plus(zenit,zx(i)))
    22402239
    22412240            if ( jion8(2,i,1)*ionsec_o2plus(zenit,zx(i)).gt.1.d-30)
    2242      $           tauo2(i,93) = 1.d0 /
     2241     $           tauo2(93,i) = 1.d0 /
    22432242     $           (jion8(2,i,1)*ionsec_o2plus(zenit,zx(i)))
    22442243
    … …  
    22762275
    22772276         do j=1,nreact
    2278             tminco2(i)  = min(tminco2(i),tauco2(i,j))
    2279             tmino2(i)   = min(tmino2(i),tauo2(i,j))
    2280             tmino3p(i)  = min(tmino3p(i),tauo3p(i,j))
    2281             tminco(i)   = min(tminco(i),tauco(i,j))
    2282             tminh(i)    = min(tminh(i),tauh(i,j))
    2283             tminoh(i)   = min(tminoh(i),tauoh(i,j))
    2284             tminho2(i)  = min(tminho2(i),tauho2(i,j))
    2285             tminh2(i)   = min(tminh2(i),tauh2(i,j))
    2286             tminh2o(i)  = min(tminh2o(i),tauh2o(i,j))
    2287             tmino1d(i)  = min(tmino1d(i),tauo1d(i,j))
    2288             tminh2o2(i) = min(tminh2o2(i),tauh2o2(i,j))
    2289             tmino3(i)   = min(tmino3(i),tauo3(i,j))
    2290             tminn(i)    = min(tminn(i),taun(i,j))
    2291             tminno(i)   = min(tminno(i),tauno(i,j))
    2292             tminn2(i)   = min(tminn2(i),taun2(i,j))
    2293             tminn2d(i)  = min(tminn2d(i),taun2d(i,j))
    2294             tminno2(i)  = min(tminno2(i),tauno2(i,j))
     2277            tminco2(i)  = min(tminco2(i),tauco2(j,i))
     2278            tmino2(i)   = min(tmino2(i),tauo2(j,i))
     2279            tmino3p(i)  = min(tmino3p(i),tauo3p(j,i))
     2280            tminco(i)   = min(tminco(i),tauco(j,i))
     2281            tminh(i)    = min(tminh(i),tauh(j,i))
     2282            tminoh(i)   = min(tminoh(i),tauoh(j,i))
     2283            tminho2(i)  = min(tminho2(i),tauho2(j,i))
     2284            tminh2(i)   = min(tminh2(i),tauh2(j,i))
     2285            tminh2o(i)  = min(tminh2o(i),tauh2o(j,i))
     2286            tmino1d(i)  = min(tmino1d(i),tauo1d(j,i))
     2287            tminh2o2(i) = min(tminh2o2(i),tauh2o2(j,i))
     2288            tmino3(i)   = min(tmino3(i),tauo3(j,i))
     2289            tminn(i)    = min(tminn(i),taun(j,i))
     2290            tminno(i)   = min(tminno(i),tauno(j,i))
     2291            tminn2(i)   = min(tminn2(i),taun2(j,i))
     2292            tminn2d(i)  = min(tminn2d(i),taun2d(j,i))
     2293            tminno2(i)  = min(tminno2(i),tauno2(j,i))
    22952294            !
    2296             tminco2plus(i) = min(tminco2plus(i),tauco2plus(i,j))
    2297             tminoplus(i)   = min(tminoplus(i),tauoplus(i,j))
    2298             tmino2plus(i)  = min(tmino2plus(i),tauo2plus(i,j))
    2299             tmincoplus(i)  = min(tmincoplus(i),taucoplus(i,j))
    2300             tmincplus(i)   = min(tmincplus(i),taucplus(i,j))
    2301             tminnplus(i)   = min(tminnplus(i),taunplus(i,j))
    2302             tminn2plus(i)  = min(tminn2plus(i),taun2plus(i,j))
    2303             tminnoplus(i)  = min(tminnoplus(i),taunoplus(i,j))
    2304             tminhplus(i)   = min(tminhplus(i),tauhplus(i,j))
    2305             tminhco2plus(i)= min(tminhco2plus(i),tauhco2plus(i,j))
     2295            tminco2plus(i) = min(tminco2plus(i),tauco2plus(j,i))
     2296            tminoplus(i)   = min(tminoplus(i),tauoplus(j,i))
     2297            tmino2plus(i)  = min(tmino2plus(i),tauo2plus(j,i))
     2298            tmincoplus(i)  = min(tmincoplus(i),taucoplus(j,i))
     2299            tmincplus(i)   = min(tmincplus(i),taucplus(j,i))
     2300            tminnplus(i)   = min(tminnplus(i),taunplus(j,i))
     2301            tminn2plus(i)  = min(tminn2plus(i),taun2plus(j,i))
     2302            tminnoplus(i)  = min(tminnoplus(i),taunoplus(j,i))
     2303            tminhplus(i)   = min(tminhplus(i),tauhplus(j,i))
     2304            tminhco2plus(i)= min(tminhco2plus(i),tauhco2plus(j,i))
    23062305         end do
    23072306
    … …  
    24412440
    24422441      subroutine timemarching( ig,i,chemthermod,n_comp_en_EQ,
    2443      $     compmin,tmin,timefrac_sec, deltat )
     2442     $     compmin,tmin,timefrac_sec, deltat,fmargin1 )
    24442443
    24452444 
    … …  
    24782477ccccccccccccccc CODE STARTS
    24792478
    2480       fmargin1=1.
     2479!      fmargin1=1.
    24812480      fmargin2=5.
    24822481
    … …  
    36613660      Lo1dtot(i) = Lo1d(i,9) + Lo1d(i,14) + Lo1d(i,19) + Lo1d(i,20) +
    36623661     $     Lo1d(i,31)
    3663      
    36643662
    36653663c     H2O2

  • trunk/LMDZ.MARS/libf/aeronomars/perosat.F

    r38 r658  
    127127c
    128128           zynew(1)=zysat(1)
     129        else
     130          pdqscloud(ig,igcm_h2o2)=0
    129131        end if
    130132c

  • trunk/LMDZ.MARS/libf/aeronomars/thermosphere.F

    r635 r658  
    7979      if (callmoldiff) THEN
    8080        call zerophys(ngridmx*nlayermx*nqmx,zdqmoldiff)
    81         call moldiff(pplay,pplev,pt,pdt,pq,pdq,ptimestep,
     81        call moldiff_red(pplay,pplev,pt,pdt,pq,pdq,ptimestep,
    8282     &                   zzlay,zdteuv,zdtconduc,zdqmoldiff)
    8383      endif
Note: See TracChangeset for help on using the changeset viewer.