Changeset 635 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Apr 26, 2012, 3:22:19 PM (13 years ago)
Author:
emillour
Message:

Mars GCM: Update of the chemistry package, including:

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

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

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

FGG+FL+EM

Location:
trunk/LMDZ.MARS/libf
Files:
5 added
10 deleted
13 edited

Legend:

Unmodified
Added
Removed

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

    r618 r635  
    9898      integer :: ig,l,i,iq,iqmax
    9999      integer :: foundswitch, lswitch
     100      integer,save :: chemthermod
    100101
    101102      integer,save :: i_co2  = 0
    … …  
    113114      integer,save :: i_n2   = 0
    114115      integer,save :: i_h2o  = 0
     116      integer,save :: i_n    = 0
     117      integer,save :: i_no   = 0
     118      integer,save :: i_no2  = 0
     119      integer,save :: i_n2d  = 0
     120      integer,save :: i_co2plus=0
     121      integer,save :: i_oplus=0
     122      integer,save :: i_o2plus=0
     123      integer,save :: i_coplus=0
     124      integer,save :: i_cplus=0
     125      integer,save :: i_nplus=0
     126      integer,save :: i_noplus=0
     127      integer,save :: i_n2plus=0
     128      integer,save :: i_hplus=0
     129      integer,save :: i_hco2plus=0
     130      integer,save :: i_elec=0
    115131
    116132      integer :: ig_vl1
    … …  
    157173
    158174         ! find index of chemical tracers to use
     175         ! Listed here are all tracers that can go into photochemistry
    159176         nbq = 0 ! to count number of tracers
     177         ! Species ALWAYS present if photochem=.T. or thermochem=.T.
    160178         i_co2 = igcm_co2
    161179         if (i_co2 == 0) then
    … …  
    248266         i_ch4 = igcm_ch4
    249267         if (i_ch4 == 0) then
    250             write(*,*) "calchim: no CH4 tracer !!!"
     268            write(*,*) "calchim: Error; no CH4 tracer !!!"
    251269            write(*,*) "CH4 will be ignored in the chemistry"
    252270         else
    … …  
    270288            niq(nbq) = i_h2o
    271289         end if
     290         !Check tracers needed for thermospheric chemistry
     291         if(thermochem) then
     292            chemthermod=0  !Default: C/O/H chemistry
     293            !Nitrogen chemistry
     294            !NO is used to determine if N chemistry is wanted
     295            !chemthermod=2 -> N chemistry
     296            i_no = igcm_no
     297            if (i_no == 0) then
     298               write(*,*) "calchim: no NO tracer"
     299               write(*,*) "C/O/H themosp chemistry only "
     300            else
     301               nbq = nbq + 1
     302               niq(nbq) = i_no
     303               chemthermod=2
     304               write(*,*) "calchim: NO in traceur.def"
     305               write(*,*) "Nitrogen chemistry included"
     306            end if
     307            ! N
     308            i_n = igcm_n
     309            if(chemthermod == 2) then
     310               if (i_n == 0) then
     311                  write(*,*) "calchim: Error; no N tracer !!!"
     312                  write(*,*) "N is needed if NO is in traceur.def"
     313                  stop
     314               else
     315                  nbq = nbq + 1
     316                  niq(nbq) = i_n
     317               end if
     318            else
     319               if (i_n /= 0) then
     320                  write(*,*) "calchim: Error: N is present, but NO is not!!!"
     321                  write(*,*) "Both must be in traceur.def if N chemistry is wanted"
     322                  stop
     323               endif
     324            endif    !Of if(chemthermod == 2)
     325            ! NO2
     326            i_no2 = igcm_no2
     327            if(chemthermod == 2) then
     328               if (i_no2 == 0) then
     329                  write(*,*) "calchim: Error; no NO2 tracer !!!"
     330                  write(*,*) "NO2 is needed if NO is in traceur.def"
     331                  stop
     332               else
     333                  nbq = nbq + 1
     334                  niq(nbq) = i_no2
     335               end if
     336            else
     337               if (i_no2 /= 0) then
     338                  write(*,*) "calchim: Error: N is present, but NO is not!!!"
     339                  write(*,*) "Both must be in traceur.def if N chemistry is wanted"
     340                  stop
     341               endif
     342            endif     !Of if(chemthermod == 2)
     343            ! N(2D)
     344            if(chemthermod == 2) then
     345               i_n2d = igcm_n2d
     346               if (i_n2d == 0) then
     347                  write(*,*) "calchim: Error; no N2D !!!"
     348                  write(*,*) "N2D is needed if NO is in traceur.def"
     349                  stop
     350               else
     351                  nbq = nbq + 1
     352                  niq(nbq) = i_n2d
     353               end if
     354            else
     355               if (i_n2d /= 0) then
     356                  write(*,*) "calchim: Error: N2D is present, but NO is not!!!"
     357                  write(*,*) "Both must be in traceur.def if N chemistry wanted"
     358                  stop
     359               endif
     360            endif    !Of if(chemthermod == 2)
     361            ! Ions
     362            ! O2+ is used to determine if ion chemistry is needed
     363            ! chemthermod=3 -> ion chemistry
     364            i_o2plus = igcm_o2plus
     365            if(chemthermod == 2) then
     366               if (i_o2plus == 0) then
     367                  write(*,*) "calchim: no O2+ tracer; no ion chemistry"
     368               else
     369                  nbq = nbq + 1
     370                  niq(nbq) = i_o2plus
     371                  chemthermod = 3
     372                  write(*,*) "calchim: O2+ in traceur.def"
     373                  write(*,*) "Ion chemistry included"
     374               end if
     375            else
     376               if (i_o2plus /= 0) then
     377                  write(*,*) "calchim: O2+ is present, but NO is not!!!"
     378                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     379                  stop
     380               endif
     381            endif   !Of if(chemthermod == 2)
     382            ! CO2+
     383            i_co2plus = igcm_co2plus
     384            if(chemthermod == 3) then
     385               if (i_co2plus == 0) then
     386                  write(*,*) "calchim: Error; no CO2+ tracer !!!"
     387                  write(*,*) "CO2+ is needed if O2+ is in traceur.def"
     388                  stop
     389               else
     390                  nbq = nbq + 1
     391                  niq(nbq) = i_co2plus
     392               end if
     393            else
     394               if (i_co2plus /= 0) then
     395                  write(*,*) "calchim: Error: CO2+ is present, but O2+ is not!!!"
     396                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     397                  stop
     398               endif
     399            endif    !Of if(chemthermod == 3)
     400            ! O+
     401            i_oplus = igcm_oplus
     402            if(chemthermod == 3) then
     403               if (i_oplus == 0) then
     404                  write(*,*) "calchim: Error; no O+ tracer !!!"
     405                  write(*,*) "O+ is needed if O2+ is in traceur.def"
     406                  stop
     407               else
     408                  nbq = nbq + 1
     409                  niq(nbq) = i_oplus
     410               end if
     411            else
     412               if (i_oplus /= 0) then
     413                  write(*,*) "calchim: Error: O+ is present, but O2+ is not!!!"
     414                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     415                  stop
     416               endif
     417            endif   !Of if (chemthermod == 3)
     418            ! CO+
     419            i_coplus = igcm_coplus
     420            if(chemthermod == 3) then
     421               if (i_coplus == 0) then
     422                  write(*,*) "calchim: Error; no CO+ tracer !!!"
     423                  write(*,*) "CO+ is needed if O2+ is in traceur.def"
     424                  stop
     425               else
     426                  nbq = nbq + 1
     427                  niq(nbq) = i_coplus
     428               end if
     429            else
     430               if (i_coplus /= 0) then
     431                  write(*,*) "calchim: Error: CO+ is present, but O2+ is not!!!"
     432                  write(*,*) " Both must be in traceur.def if ionosphere wanted"
     433                  stop
     434               endif
     435            endif   ! Of if (chemthermod == 3)
     436            ! C+
     437            i_cplus = igcm_cplus
     438            if(chemthermod == 3) then
     439               if (i_cplus == 0) then
     440                  write(*,*) "calchim: Error; no C+ tracer !!!"
     441                  write(*,*) "C+ is needed if O2+ is in traceur.def"
     442                  stop
     443               else
     444                  nbq = nbq + 1
     445                  niq(nbq) = i_cplus
     446               end if
     447            else
     448               if (i_cplus /= 0) then
     449                  write(*,*) "calchim: Error; C+ is present, but O2+ is not!!!"
     450                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     451                  stop
     452               endif
     453            endif   ! Of if (chemthermod == 3)
     454            ! N+
     455            i_nplus = igcm_nplus
     456            if(chemthermod == 3) then
     457               if (i_nplus == 0) then
     458                  write(*,*) "calchim: Error; no N+ tracer !!!"
     459                  write(*,*) "N+ is needed if O2+ is in traceur.def"
     460                  stop
     461               else
     462                  nbq = nbq + 1
     463                  niq(nbq) = i_nplus
     464               end if
     465            else
     466               if (i_nplus /= 0) then
     467                  write(*,*) "calchim: Error: N+ is present, but O2+ is not!!!"
     468                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     469                  stop
     470               endif
     471            endif   !Of if (chemthermod == 3)
     472            ! NO+
     473            i_noplus = igcm_noplus
     474            if(chemthermod == 3) then
     475               if (i_noplus == 0) then
     476                  write(*,*) "calchim: Error; no NO+ tracer !!!"
     477                  write(*,*) "NO+ is needed if O2+ is in traceur.def"
     478                  stop
     479               else
     480                  nbq = nbq + 1
     481                  niq(nbq) = i_noplus
     482               end if
     483            else
     484               if (i_noplus /= 0) then
     485                  write(*,*) "calchim: Error: NO+ is present, but O2+ is not!!!"
     486                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     487                  stop
     488               endif
     489            endif   !Of if (chemthermod == 3)
     490            ! N2+
     491            i_n2plus = igcm_n2plus
     492            if (chemthermod == 3) then
     493               if (i_n2plus == 0) then
     494                  write(*,*) "calchim: Error; no N2+ tracer !!!"
     495                  write(*,*) "N2+ is needed if O2+ is in traceur.def"
     496                  stop
     497               else
     498                  nbq = nbq + 1
     499                  niq(nbq) = i_n2plus
     500               end if
     501            else
     502               if (i_n2plus /= 0) then
     503                  write(*,*) "calchim: Error: N2+ is present, but O2+ is not!!!"
     504                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     505                  stop
     506               endif
     507            endif   !Of if (chemthermod == 3)
     508            !H+
     509            i_hplus = igcm_hplus
     510            if (chemthermod == 3) then
     511               if (i_hplus == 0) then
     512                  write(*,*) "calchim: Error; no H+ tracer !!!"
     513                  write(*,*) "H+ is needed if O2+ is in traceur.def"
     514                  stop
     515               else
     516                  nbq = nbq + 1
     517                  niq(nbq) = i_hplus
     518               end if
     519            else
     520               if (i_hplus /= 0) then
     521                  write(*,*) "calchim: Error: H+ is present, but O2+ is not!!!"
     522                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     523                  stop
     524               endif
     525            endif   !Of if (chemthermod == 3)
     526            ! HCO2+
     527            i_hco2plus = igcm_hco2plus
     528            if(chemthermod == 3) then
     529               if (i_hco2plus == 0) then
     530                  write(*,*) "calchim: Error; no HCO2+ tracer !!!"
     531                  write(*,*) "HCO2+ is needed if O2+ is in traceur.def"
     532                  stop
     533               else
     534                  nbq = nbq + 1
     535                  niq(nbq) = i_hco2plus
     536               end if
     537            else
     538               if (i_hco2plus /= 0) then
     539                  write(*,*) "calchim: Error: HCO2+ is present, but O2+ is not!!!"
     540                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     541                  stop
     542               endif
     543            endif    !Of if(chemthermod == 3)
     544            !e-
     545            i_elec = igcm_elec
     546            if(chemthermod == 3) then
     547               if (i_elec == 0) then
     548                  write(*,*) "calchim: Error; no e- tracer !!!"
     549                  write(*,*) "e- is needed if O2+ is in traceur.def"
     550                  stop
     551               else
     552                  nbq = nbq + 1
     553                  niq(nbq) = i_elec
     554               end if
     555            else
     556               if(i_elec /= 0) then
     557                  write(*,*) "calchim: Error: e- is present, but O2+ is not!!!"
     558                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     559                  stop
     560               endif
     561            endif   !Of if (chemthermod == 3)
     562         endif      !Of thermochem
    272563
    273564         write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
    274          
     565              
    275566         firstcall = .false.
    276567      end if ! if (firstcall)
    … …  
    279570
    280571      zycol(:,:)    = 0.
    281       dqchim(:,:,:) = 0
    282       dqschim(:,:)  = 0
     572      dqchim(:,:,:) = 0.
     573      dqschim(:,:)  = 0.
    283574
    284575!     latvl1= 22.27
    … …  
    292583
    293584      do ig = 1,ngridmx
    294 
     585         
    295586         foundswitch = 0
    296587         do l = 1,nlayermx
    … …  
    316607
    317608            if (photochem .and. thermochem) then
    318                if (foundswitch == 0 .and. pplay(ig,l) < 1.e-3) then
     609               if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then
    319610                  lswitch = l
    320611                  foundswitch = 1
    … …  
    357648
    358649!        chemistry in upper atmosphere
    359 
     650       
    360651         if (thermochem) then
    361             call chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress,  &
    362                              zlocal,szacol,ptimestep,zday)
     652            call chemthermos(ig,lswitch,chemthermod,zycol,ztemp,zdens,  &
     653                             zpress,zlocal,szacol,ptimestep,zday)
    363654         end if
    364655
    … …  
    392683            end do
    393684         end do ! of do l = 1,nlayermx
     685!         if(ig.eq.800)write(*,*)'calchim/686',dqchim(ig,27,23)
    394686
    395687!=======================================================================
    … …  
    419711      return
    420712      end
     713

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

    r618 r635  
    132132            cpi(nbq) = 1.870e3
    133133         end if
     134         if (igcm_n /= 0) then
     135            nbq = nbq + 1
     136            niq(nbq) = igcm_n
     137            aki(nbq) = 0.0
     138            cpi(nbq) = 0.0
     139         endif
     140         if(igcm_no /= 0) then
     141            nbq = nbq + 1
     142            niq(nbq) = igcm_no
     143            aki(nbq) = 0.0
     144            cpi(nbq) = 0.0
     145         endif
     146         if(igcm_no2 /= 0) then
     147            nbq = nbq + 1
     148            niq(nbq) = igcm_no2
     149            aki(nbq) = 0.0
     150            cpi(nbq) = 0.0
     151         endif
     152         if(igcm_n2d /= 0) then
     153            nbq = nbq + 1
     154            niq(nbq) = igcm_n2d
     155            aki(nbq) = 0.0
     156            cpi(nbq) = 0.0
     157         endif
     158         if(igcm_co2plus /= 0) then
     159            nbq = nbq + 1
     160            niq(nbq) = igcm_co2plus
     161            aki(nbq) = 0.0
     162            cpi(nbq) = 0.0
     163         endif
     164         if(igcm_oplus /= 0) then
     165            nbq = nbq + 1
     166            niq(nbq) = igcm_oplus
     167            aki(nbq) = 0.0
     168            cpi(nbq) = 0.0
     169         endif
     170         if(igcm_o2plus /= 0) then
     171            nbq = nbq + 1
     172            niq(nbq) = igcm_o2plus
     173            aki(nbq) = 0.0
     174            cpi(nbq) = 0.0
     175         endif
     176         if(igcm_coplus /= 0) then
     177            nbq = nbq + 1
     178            niq(nbq) = igcm_coplus
     179            aki(nbq) = 0.0
     180            cpi(nbq) = 0.0
     181         endif
     182         if(igcm_cplus /= 0) then
     183            nbq = nbq + 1
     184            niq(nbq) = igcm_cplus
     185            aki(nbq) = 0.0
     186            cpi(nbq) = 0.0
     187         endif
     188         if(igcm_nplus /= 0) then
     189            nbq = nbq + 1
     190            niq(nbq) = igcm_nplus
     191            aki(nbq) = 0.0
     192            cpi(nbq) = 0.0
     193         endif
     194         if(igcm_noplus /= 0) then
     195            nbq = nbq + 1
     196            niq(nbq) = igcm_noplus
     197            aki(nbq) = 0.0
     198            cpi(nbq) = 0.0
     199         endif
     200         if(igcm_n2plus /= 0) then
     201            nbq = nbq + 1
     202            niq(nbq) = igcm_n2plus
     203            aki(nbq) = 0.0
     204            cpi(nbq) = 0.0
     205         endif
     206         if(igcm_hplus /= 0) then
     207            nbq = nbq + 1
     208            niq(nbq) = igcm_hplus
     209            aki(nbq) = 0.0
     210            cpi(nbq) = 0.0
     211         endif
     212         if(igcm_hco2plus /= 0) then
     213            nbq = nbq + 1
     214            niq(nbq) = igcm_hco2plus
     215            aki(nbq) = 0.0
     216            cpi(nbq) = 0.0
     217         endif
     218         
     219           
    134220
    135221         firstcall = .false.

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

    r38 r635  
    11c**********************************************************************
    22
    3       subroutine hrtherm
    4      $ (co2x,o2x,o3px,h2x,h2ox,h2o2x,aux1,aux2,tx,nz,iz,date,zenit,jtot)
     3      subroutine hrtherm(ig,euvmod,rm,nespeuv,tx,iz,zenit,jtot)
    54
    65
    … …  
    1413c     common variables and constants
    1514
     15
     16      include 'dimensions.h'
     17      include 'dimphys.h'
    1618      include 'param.h'
    17       include 'param_v3.h'
    18       include 'callkeys.h'
     19      include 'param_v4.h'
    1920
    2021
    2122c    local parameters and variables
    2223
    23       real       xabsi(nabs,nzmax)                      !densities
    24       real       nada
    25       real       jergs(ninter,nabs,nzmax)
     24      real       xabsi(nabs,nlayermx)                   !densities
     25      real       jergs(ninter,nabs,nlayermx)
    2626     
    2727      integer    i,j,k,indexint          !indexes
    … …  
    3131c     input and output variables
    3232
    33       integer    nz                          !number of layers
    34       real       co2x(nz)                    !density of CO2(cm^-3)
    35       real       o2x(nz)                     !density of O2(cm^-3)
    36       real       o3px(nz)                    !density of O(3P)(cm^-3)
    37       real       h2x(nz)                     !density of H2(cm^-3)
    38       real       h2ox(nz)                    !density of H2O(cm^-3)
    39       real       h2o2x(nz)                   !density of H2O2(cm^-3)
    40       real       aux1(nz)                            !auxiliar variable
    41       real       aux2(nz)                    !auxiliar variable
    42       real       jtot(nz)                    !output: heating rate(erg/s)
    43       real       tx(nz)                      !temperature
    44       real       date
     33      integer    ig  ,euvmod
     34      integer    nespeuv
     35      real       rm(nlayermx,nespeuv)              !density matrix (cm^-3)
     36      real       jtot(nlayermx)                    !output: heating rate(erg/s)
     37      real       tx(nlayermx)                      !temperature
    4538      real       zenit
    46       real       iz(nz)
     39      real       iz(nlayermx)
    4740
    48       logical firstcall
    49       save firstcall
    50       data firstcall /.true./
     41      ! tracer indexes for the EUV heating:
     42!!! ATTENTION. These values have to be identical to those in chemthermos.F90
     43!!! If the values are changed there, the same has to be done here  !!!
     44      integer,parameter :: i_co2=1
     45      integer,parameter :: i_o2=2
     46      integer,parameter :: i_o=3
     47      integer,parameter :: i_co=4
     48      integer,parameter :: i_h=5
     49      integer,parameter :: i_h2=8
     50      integer,parameter :: i_h2o=9
     51      integer,parameter :: i_h2o2=10
     52      integer,parameter :: i_o3=12
     53      integer,parameter :: i_n2=13
     54      integer,parameter :: i_n=14
     55      integer,parameter :: i_no=15
     56      integer,parameter :: i_no2=17
    5157
    5258c*************************PROGRAM STARTS*******************************
    5359
    54 c      if (firstcall) then
    55 c        if(.not. thermochem) call param_read
    56 c        firstcall= .false.
    57 c      endif
    58 
    59       if(zenit.gt.100.) then
     60      !If nighttime, photoabsorption coefficient set to 0
     61      if(zenit.gt.140.) then
    6062         dn='n'
    6163         else
    … …  
    6365      end if
    6466      if(dn.eq.'n') then
    65         do i=1,nz                                   
     67        do i=1,nlayermx                                   
    6668              jtot(i)=0.
    6769        enddo       
    6870        return
    6971      endif
    70      
    71       do i=1,nz
    72          xabsi(1,i) = co2x(i)
    73          xabsi(2,i) = o2x(i)
    74          xabsi(3,i) = o3px(i)
    75          xabsi(4,i) = h2ox(i)
    76          xabsi(5,i) = h2x(i)
    77          xabsi(6,i) = h2o2x(i)
    78          jtot(i) = 0.
     72
     73      !initializations
     74      jergs(:,:,:)=0.
     75      xabsi(:,:)=0.
     76      jtot(:)=0.
     77      !All number densities to a single array, xabsi(species,layer)
     78      do i=1,nlayermx
     79         xabsi(1,i)  = rm(i,i_co2)
     80         xabsi(2,i)  = rm(i,i_o2)
     81         xabsi(3,i)  = rm(i,i_o)
     82         xabsi(4,i)  = rm(i,i_h2o)
     83         xabsi(5,i)  = rm(i,i_h2)
     84         xabsi(6,i)  = rm(i,i_h2o2)
     85         !Only if O3, N or ion chemistry requested
     86         if(euvmod.ge.1) then
     87            xabsi(7,i)  = rm(i,i_o3)
     88         endif
     89         !Only if N or ion chemistry requested
     90         if(euvmod.ge.2) then
     91            xabsi(8,i)  = rm(i,i_n2)
     92            xabsi(9,i)  = rm(i,i_n)
     93            xabsi(10,i) = rm(i,i_no)
     94            xabsi(13,i) = rm(i,i_no2)
     95         endif
     96         xabsi(11,i) = rm(i,i_co)
     97         xabsi(12,i) = rm(i,i_h)
    7998      end do
    8099
    81       if(.not. thermochem) then
    82         call jthermcalc
    83      $ (co2x,o2x,o3px,h2x,h2ox,h2o2x,aux1,aux2,tx,nz,iz,date,zenit)
    84       endif
     100      !Calculation of photoabsortion coefficient
     101      call jthermcalc(ig,euvmod,rm,nespeuv,tx,iz,zenit)
    85102
    86       do i=1,nz
     103      !Total photoabsorption coefficient
     104      do i=1,nlayermx
     105         jtot(i)=0.
    87106        do j=1,nabs
    88           do indexint=1,33
     107          do indexint=1,ninter
    89108            jergs(indexint,j,i) = jfotsout(indexint,j,i)
    90109     $              * xabsi (j,i) * fluxtop(indexint) 

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

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

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

    r38 r635  
    1 c***********************************************
     1!***********************************************
    22
    3 c       param.par
     3!       param.par
    44
    5 c       Parameters for paramhr.f
    6 c***********************************************
     5!       Parameters for paramhr.f
     6!***********************************************
    77
    88        integer    ninter
    9         parameter  (ninter=33)
     9        parameter  (ninter=36)
    1010
    1111        integer    nabs
    12         parameter  (nabs=6)
     12        parameter  (nabs=13)
    1313
    1414        integer    nz2
    15         parameter  (nz2=203)
     15        parameter  (nz2=253)
    1616
    1717        integer    ninter2
    … …  
    3333        parameter  (radio = 3390.)           !radio de Marte(km)
    3434
    35         real       pmco2
    36         parameter  (pmco2 = 44.)
    37 
    38         real       pmo3p
    39         parameter  (pmo3p = 16.)
    40 
    41         real       pmo2
    42         parameter  (pmo2 = 32.)
    43 
    44         real       pmh2o2
    45         parameter  (pmh2o2 = 36.)
    46 
    47         real       cpco2
    48         parameter  (cpco2=20.0e6)
    49 
    50         real       cpo2
    51         parameter  (cpo2=0.9194e7)
    52 
    53         real       cpo3p
    54         parameter  (cpo3p=1.30e7)
    55 
    5635        integer    nreact
    57         parameter  (nreact=22)
    58 
    59         integer    nzmax
    60         parameter  (nzmax=210)
     36        parameter  (nreact=93)
    6137
    6238
     39
     40

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

    r38 r635  
    55
    66c     common variables and constants
    7 
    8 #include "param.h"
    9 #include "param_v3.h"
    10 #include "datafile.h"
    11 
    12 
     7      include "dimensions.h"
     8      include "dimphys.h"
     9      include 'param.h'
     10      include 'param_v4.h'
     11      include 'datafile.h'
     12 
     13 
    1314c     local variables
    1415
    1516      integer    i,j,k,inter                          !indexes
    16       integer ierr
     17      integer ierr,lnblnk
    1718      real       nada
    18       logical,save :: firstcall=.true.
    19 
    20      
     19 
     20     
    2121c*************************+PROGRAM STARTS**************************
    2222
    2323
    24 c     data
    25       if (firstcall) then
    26       anchint(:)=(/4.9,5.3,6.1,12.8,2.0,6.9,11.8,0.4,1.9,4.9,0.9,4.0,
    27      $6.9,4.9,4.9,10.2,5.6,1.9,0.9,0.9,0.9,0.9,3.7
    28      $,13.1,3.9,6.0,4.8,23.9,8.4,8.5,26.8,7.5,127.6/)
    29 
    30       crscabsi2(1,1:ninter2)=(/0.0,4.42e-18,7.51e-18
    31      $,1.498e-17,2.344e-17,2.388e-17,2.927E-17,3.161e-17,3.161e-17,
    32      $3.4e-17,3.4e-17,2.588e-17,2.596e-17,2.248e-17,3.183e-17,
    33      $1.284e-17/)
    34 
    35       crscabsi2(2,1:ninter2)=(/2.736e-19,1.055e-18,4.789e-18
    36      $,1.063e-17,1.559e-17,1.698e-17,2.164e-17,2.241e-17,2.427e-17,
    37      $2.444e-17,2.502e-17,2.516e-17,2.907e-17,3.649e-17,3.066e-17,
    38      $2.093e-17/)
    39 
    40       crscabsi2(3,1:ninter2)=(/2.603e-19,6.985e-19,3.97e-18
    41      $,6.621e-18,8.433e-18,9.011e-18,9.862e-18,1.033e-17,1.012e-17,
    42      $1.033e-17,1.033e-17,1.034e-17,9.0e-18,3.564e-18,3.572e-18
    43      $,3.47E-18/)
    44 
    45       freccen(:)=(/3.4,7.5,14.5,23.0,30.3,34.1,49.6,50.5,52.5,56.0,
     24c     data for the UV heating tabulation
     25
     26      data (crscabsi2(1,j),j=1,16) /5.61031E-19,1.59677E-18,4.7072E-18,
     27     $     1.48254e-17,2.07445e-17,2.573e-17,2.901e-17,3.083e-17,
     28     $     3.217e-17,3.539e-17,3.658e-17,3.63e-17,3.41239e-17,
     29     $     2.71019e-17,4.93677e-17,1.64e-17/
     30
     31      data (crscabsi2(2,j),j=1,16) /0.27250E-18,0.11650E-17,0.39250E-17,
     32     $     0.10630E-16,0.15590E-16,0.17180E-16,0.19270E-16,0.22860E-16,
     33     $     0.24270E-16,0.24440E-16,0.25020E-16,0.26600E-16,0.25400E-16,
     34     $     0.35800E-16,0.25590E-16,0.16740E-16/
     35
     36      data (crscabsi2(3,j),j=1,16) /0.2776E-18,0.9792E-18,0.3313E-17,
     37     $     0.6621E-17,0.8481E-17,0.9146E-17,0.9414E-17,0.1039E-16,
     38     $     0.1012E-16,0.1033E-16,0.1033E-16,0.1033E-16,0.8268E-17,
     39     $     0.6563E-17,0.3506E-17,0.3470E-17/
     40
     41      data (crscabsi2(5,j),j=1,16) /.5E-20,.1077607E-19,.5670491E-19,
     42     $     .3322716E-18,.1054509E-17,.1700005E-17,.3171188E-17,
     43     $     .4734241E-17,.5108741E-17,.6022236E-17,.6741537E-17,
     44     $     .7277079E-17,.9070787E-17,.9708916E-17,.4026281E-17,0.0/
     45
     46      data (crscabsi2(8,j),j=1,16) /0.0, 7.44175e-19, 2.23167e-18,
     47     $    8.46200e-18,1.18275e-17,1.54900e-17,2.32475e-17,2.41373e-17,
     48     $     2.55482e-17,2.38431e-17,2.28600e-17,2.35067e-17,2.56000e-17,
     49     $     2.64636e-17,2.86260e-17,3.26561e-17/
     50
     51      data(crscabsi2(9,j),j=1,16) /3.48182e-20,3.37038e-19,1.03077e-18,
     52     $     4.01364e-18,6.45e-18,7.8e-18,1.0e-17,1.13500e-17,1.15500e-17,
     53     $     1.18000e-17,1.17500e-17,1.16000e-17,1.28667e-17,1.18500e-17,
     54     $     1.11000e-17,9.50000e-18/
     55
     56      data(crscabsi2(10,j),j=1,16) /0.0,9.39833e-19,2.87714e-18,
     57     $     9.66900e-18,1.37063e-17,1.61820e-17,2.30450e-17,2.63373e-17,
     58     $     2.63773e-17,2.67677e-17,2.64100e-17,2.53000e-17,2.18100e-17,
     59     $     2.04941e-17,2.28160e-17,2.93550e-17/
     60
     61      data(crscabsi2(11,j),j=1,16) /0.0,9.58555e-19,2.52767e-18,
     62     $     8.29700e-18,1.21850e-17,1.40500e-17,1.97025e-17,2.12018e-17,
     63     $     2.14673e-17,2.20331e-17,2.21500e-17,2.21600e-17,2.33200e-17,
     64     $     2.67800e-17,2.56400e-17,3.58561e-17/
     65
     66      data(crscabsi2(12,j),j=1,16) /0.0,1.0e-20,2.5e-20,1.30400e-19,
     67     $     2.93800e-19,4.36000e-19,8.49400e-19,1.29400e-18,1.40500e-18,
     68     $     1.67600e-18,1.93400e-18,2.12200e-18,2.75800e-18,3.48400e-18,
     69     $     4.17200e-18,5.26000e-18/
     70
     71      data(crscabsi2(13,j),j=1,16) /0.0,1.60e-18,4.99111e-18,1.48496e-17
     72     $     ,2.17395e-17,2.55857e-17,2.87754e-17,3.65571e-17,3.85691e-17,
     73     $     4.16286e-17,4.15117e-17,4.05901e-17,3.64000e-17,2.99670e-17,
     74     $     2.46796e-17,2.51789e-17/
     75
     76      data freccen /3.4,7.5,14.5,23.0,30.3,34.1,49.6,50.5,52.5,56.0,
    4677     $59.0,61.5,68.7,73.1,78.4,83.1,92.4,97.5,99.3,100.1,100.7,102.1,
    4778     $104.5,116.8,121.3,127.0,130.6,153.7,162.8,171.4
    48      $,195.6,206.3,273.5/)
    49 
    50       co2crsc195(:)=(/3.691e-19,4.44216e-20,3.86945e-19,5.94208e-19,
    51      $2.93217e-19,7.58769e-20,8.60192e-21,4.20007e-24,2.29996e-26/)
    52 
    53       co2crsc295(:)=(/3.691e-19,5.21572e-20,4.23488e-19,6.54728e-19,
    54      $3.30227e-19,1.03183e-19,1.55722e-20,1.72317e-23,7.0e-25/)
    55      
    56       firstcall=.false.
    57       endif ! of if (firstcall)
     79     $,195.6,206.3,222.0,236.0,289.0,600./
     80
     81      data co2crsc195/2.05864e-17,5.90557e-20,3.1027e-19,6.70653e-19,
     82     $4.55132e-19,8.87122e-20,1.32138e-20,7.22244e-23,2.88002e-26/
     83
     84      data co2crsc295/2.05897e-17,6.71104e-20,3.45509e-19,7.45711e-19,
     85     $4.82752e-19,1.11594e-19,1.98308e-20,1.3853e-22,2.1414e-25/
    5886
    5987c     Reads tabulated functions
    6088
    61       open(210,status='old',file=trim(datafile)//'/EUVDAT/coln.dat',
    62      & iostat=ierr)
    63 
    64       IF (ierr.NE.0) THEN
     89      !Tabulated column amount
     90      open(210, status = 'old',
     91c    $file=datafile(1:lnblnk(datafile))//'/EUVDAT/coln.dat',iostat=ierr)
     92     $file=datafile
     93     $   (1:lnblnk(datafile))//'/EUVDAT/param_v5/coln.dat',iostat=ierr)
     94
     95      IF (ierr.NE.0) THEN
    6596       write(*,*)'cant find directory EUVDAT and content coln.dat'
    66        write(*,*)'(in aeronomars/param_read.F)'
    67        write(*,*)'It should be in :',trim(datafile),'/'
     97       write(*,*)'(in aeronomars/param_read_iono.F)'
     98       write(*,*)'It should be in :', datafile,'/'
    6899       write(*,*)'1) You can change this directory address in '
    69100       write(*,*)'   file phymars/datafile.h'
    … …  
    73104       STOP
    74105      ENDIF
    75 
    76       open(220,file=trim(datafile)//'/EUVDAT/j2_an.dat')
    77       open(230,file=trim(datafile)//'/EUVDAT/j3_an.dat')
    78       open(240,file=trim(datafile)//'/EUVDAT/j1_an.dat')
    79       open(250,file=trim(datafile)//'/EUVDAT/j2_bn.dat')
    80       open(260,file=trim(datafile)//'/EUVDAT/j2_cn.dat')
    81       open(270,file=trim(datafile)//'/EUVDAT/j1_bn.dat')
    82       open(280,file=trim(datafile)//'/EUVDAT/j1_cn.dat')
    83       open(290,file=trim(datafile)//'/EUVDAT/j1_dn.dat')
    84       open(150,file=trim(datafile)//'/EUVDAT/j4n.dat')
    85       open(160,file=trim(datafile)//'/EUVDAT/j5n.dat')
    86       open(170,file=trim(datafile)//'/EUVDAT/j6n.dat')
    87 
    88       do i=210,290,10
    89          read(i,*)
    90          read(i,*)
    91       end do
    92 
    93       do i=150,170,10
    94          read(i,*)
    95          read(i,*)
    96       end do
    97 
    98       do i=nz2,1,-1
    99         read(210,*) c23(i),(c123(j,i),j=2,ninter2),c12(i),c1(i),ch2o2(i)
    100       end do
    101 
    102       do i=nz2,1,-1
    103          read(220,*) (jabsifotsint(j,2,i),j=1,ninter2)
    104       end do
    105      
    106       do i=nz2,1,-1
    107          read(230,*) (jabsifotsint(j,3,i),j=1,ninter2)
    108       end do
    109 
    110       do i=nz2,1,-1
    111          read(240,*) (jabsifotsint(j,1,i),j=2,ninter2)
    112       end do
    113 
    114       do i=nz2,1,-1
    115          read(250,*) (jabsifotsint(j,2,i),j=17,24)
    116       end do
    117 
    118       do i=nz2,1,-1
    119          read(260,*) (jabsifotsint(j,2,i),j=25,31)
    120       end do
    121 
    122       do i=nz2,1,-1
    123          read(270,*) (jabsifotsint(j,1,i),j=17,24)
    124       end do
    125 
    126       do i=nz2,1,-1
    127          read(280,*) (jabsifotsint(j,1,i),j=25,31)
    128       end do
    129 
    130       do i=nz2,1,-1
    131          read(290,*) jabsifotsint(32,1,i)
    132       end do
    133 
    134       do i=nz2,1,-1
    135          read(160,*) (jabsifotsint(j,5,i),j=1,15)
    136       end do
    137 
    138       do i=nz2,1,-1
    139          read(150,*) (jabsifotsint(j,4,i),j=25,31)
    140       end do
    141 
    142       do i=nz2,1,-1
    143          read(170,*) (jabsifotsint(j,6,i),j=25,33)
    144       end do
    145 
    146       do i=210,290,10
     106 
     107      !Tabulated photoabsorption coefficients
     108      open(220,file=datafile
     109     $     (1:lnblnk(datafile))//'/EUVDAT/param_v5/j2_an.dat')
     110      open(230,file=datafile
     111     $     (1:lnblnk(datafile))//'/EUVDAT/param_v5/j3_an.dat')
     112      open(240,file=datafile
     113     $     (1:lnblnk(datafile))//'/EUVDAT/param_v5/j1_an.dat')
     114      open(250,file=datafile
     115     $     (1:lnblnk(datafile))//'/EUVDAT/param_v5/j2_bn.dat')
     116      open(260,file=datafile
     117     $     (1:lnblnk(datafile))//'/EUVDAT/param_v5/j2_cn.dat')
     118      open(300,file=datafile
     119     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j2_dn.dat')
     120      open(270,file=datafile
     121     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j1_bn.dat')
     122      open(280,file=datafile
     123     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j1_cn.dat')
     124      open(290,file=datafile
     125     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j1_dn.dat')
     126      open(150,file=datafile
     127     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j4n.dat')
     128      open(160,file=datafile
     129     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j5n.dat')
     130      open(170,file=datafile
     131     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j6n.dat')
     132      open(180,file=datafile
     133     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j7n.dat')
     134      open(390,file=datafile
     135     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j8_an.dat')
     136      open(400,file=datafile
     137     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j8_bn.dat')
     138      open(410,file=datafile
     139     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j9n.dat')
     140      open(420,file=datafile
     141     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j10_an.dat')
     142      open(430,file=datafile
     143     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j10_bn.dat')
     144      open(440,file=datafile
     145     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j10_cn.dat')
     146      open(450,file=datafile
     147     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j11_an.dat')
     148      open(460,file=datafile
     149     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j11_bn.dat')
     150      open(470,file=datafile
     151     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j11_cn.dat')
     152      open(480,file=datafile
     153     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j12n.dat')
     154      open(490,file=datafile
     155     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j13_an.dat')
     156      open(500,file=datafile
     157     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j13_bn.dat')
     158      open(510,file=datafile
     159     $     (1:lnblnk(datafile))//'/EUVDAT//param_v5/j13_cn.dat')
     160
     161     
     162      do i=210,300,10
     163         read(i,*)
     164         read(i,*)
     165      end do
     166
     167      do i=150,180,10
     168         read(i,*)
     169         read(i,*)
     170      end do
     171
     172      do i=390,510,10
     173         read(i,*)
     174         read(i,*)
     175      enddo
     176
     177      do i=nz2,1,-1
     178         read(210,*) (c1_16(i,j),j=1,16),c17_24(i),c25_29(i),c30_31(i),
     179     $        c32(i),c33(i),c34(i),c35(i),c36(i)
     180      end do
     181
     182      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)
     285      enddo
     286
     287      do i=210,300,10
    147288         close(i)
    148289      end do
    149290
    150       do i=150,170,10
     291      do i=150,180,10
    151292         close(i)
    152293      end do
    153294
    154 
    155 c     reads t0
    156 
    157       open(120,file=trim(datafile)//'/EUVDAT/t0.dat')
    158       do i=1,201
    159          read(120,*)t0(i)
    160       end do
    161       close(120)
    162 
    163       open(100,file=trim(datafile)//'/EUVDAT/flujo.dat')
     295      do i=390,510,10
     296         close(i)
     297      enddo
     298
     299
     300c     set t0
     301
     302      do i=1,nz2
     303         t0(i)=195.
     304      end do
     305
     306
    164307      do i=1,ninter
    165          read(100,*) inter,fluxtophr(i)
     308         fluxtop(i)=1.
     309      end do
     310
     311      !Parameters for the variation of the solar flux with 11 years cycle
     312      open(100,file=datafile
     313     $     (1:lnblnk(datafile))//'/EUVDAT/param_v5/varflujo.dat')
     314      read(100,*)
     315      do i=1,24
     316         read(100,*) inter,ct1(i),p1(i),ct2(i),p2(i),nada
    166317      end do
    167318      close(100)
    168319
    169 
    170       open(99,file=trim(datafile)//'/EUVDAT/varflujo.dat')
    171       read(99,*)
    172       do i=1,24
    173          read(99,*) inter,ct1(i),p1(i),ct2(i),p2(i),nada
    174       end do
    175       close(99)
    176 
    177 c     eficiencias de disociacion (de Torr et al, 1979)
    178 
    179       do inter=1,11
    180          efdisco2(inter) = 0.
    181          efdiso2(inter) = 0.
    182       end do
    183 
     320c     dissociation and ionization efficiencies
     321
     322      do inter=1,ninter
     323         efdisco2(inter)=0.
     324         efdiso2(inter)=0.
     325         efdish2(inter)=0.
     326         efdish2o(inter)=0.
     327         efdish2o2(inter)=0.
     328         efdiso3(inter)=0.
     329         efdisco(inter)=0.
     330         efdisn2(inter)=0.
     331         efdisno(inter)=0.
     332         efdisno2(inter)=0.
     333         efionco2(inter,1)=0.
     334         efionco2(inter,2)=0.
     335         efionco2(inter,3)=0.
     336         efionco2(inter,4)=0.
     337         efiono3p(inter)=0.
     338         efionn2(inter,1)=0.
     339         efionn2(inter,2)=0.
     340         efionco(inter,1)=0.
     341         efionco(inter,2)=0.
     342         efionco(inter,3)=0.
     343         efionn(inter)=0.
     344         efionh(inter)=0.
     345         efionno(inter)=0.
     346      enddo
     347
     348
     349c     CO2, O2, NO
     350
     351      open(120,file=datafile
     352     $     (1:lnblnk(datafile))//'/EUVDAT/param_v5/efdis_inter.dat')
     353      read(120,*)
     354!      do i=1,21
     355!         read(120,*)inter,efdisco2(inter),efdiso2(inter),efdisno(inter)
     356      do inter=8,28
     357         read(120,*)i,efdisco2(inter),efdiso2(inter),efdisno(inter)
     358      enddo
     359      do inter=29,ninter
     360         efdisco2(inter)=1.
     361         efdiso2(inter)=1.
     362         efdisno(inter)=1.
     363      enddo
     364
     365
     366c     N2
     367
     368      efdisn2(15)=0.1
     369      do inter=16,ninter
     370         efdisn2(inter)=1.
     371      enddo
     372
     373
     374c     CO
     375
     376      efdisco(16)=0.5
     377      do inter=17,ninter
     378         efdisco(inter)=1.
     379      enddo
     380
     381     
     382c     O, N, H
     383
     384      do inter=1,ninter
     385         efdiso(inter)=0.
     386         efdisn(inter)=0.
     387         efdish(inter)=0.
     388      enddo
     389
     390
     391c     H2O, H2O2, O3, NO2
     392
     393      do inter=25,31
     394         efdish2o(inter)=1.
     395      enddo
     396      do inter=25,35
     397         efdish2o2(inter)=1.
     398      enddo
     399      do inter=34,36
     400         efdiso3(inter)=1.
     401      enddo
     402      do inter=27,36
     403         efdisno2(inter)=1.
     404      enddo
    184405      do inter=1,15
    185          efdish2(inter) = 1.
    186       end do
    187 
    188       efdisco2(12) = 0.183
    189       efdiso2(12) = 0.003
    190 
    191       efdisco2(13) = 0.163
    192       efdiso2(13) = 0.170
    193 
    194       efdisco2(14) = 0.243
    195       efdiso2(14) =0.180
    196 
    197       efdisco2(15) = 0.323
    198       efdiso2(15) =0.653
    199 
    200       efdisco2(16) = 0.235
    201       efdiso2(16) =0.616
    202 
    203       do inter=17,32
    204          efdisco2(inter) = 1.0
    205       end do
    206 
    207       efdiso2(17) = 0.399
    208       do inter=18,20
    209          efdiso2(inter) = 0.261
    210       end do
    211 
    212       do inter=21,23
    213          efdiso2(inter) = 0.755
    214       end do
    215 
    216 
    217       do inter=24,31
    218          efdiso2(inter) = 1.
    219       end do
    220 
    221       do inter=25,31
    222          efdish2o(inter) = 1.
    223       end do
    224 
    225       do inter=25,33
    226          efdish2o2(inter) = 1.
    227       end do
     406         efdish2(inter)=1.
     407      enddo
     408         
     409      !4 possible channels for CO2 ionization
     410      do inter=14,16
     411         efionco2(inter,1)=1.-efdisco2(inter)
     412      enddo
     413      efionco2(13,1)=0.805*(1.-efdisco2(13))
     414      efionco2(13,2)=0.195*(1.-efdisco2(13))
     415      do inter=11,12
     416         efionco2(inter,3)=1.-efdisco2(inter)
     417      enddo
     418      efionco2(10,3)=0.9*(1.-efdisco2(10))
     419      efionco2(10,4)=0.1*(1.-efdisco2(10))
     420      do inter=2,9
     421         efionco2(inter,4)=1.-efdisco2(inter)
     422      enddo
     423
     424      !For O(3p) total ionization under 91.1 nm
     425      do inter=1,16
     426         efiono3p(inter)=1.d0
     427      enddo
     428
     429      !2 channels for N2 ionization
     430      do inter=9,15
     431         efionn2(inter,1)=1.-efdisn2(inter)
     432      enddo
     433      do inter=2,8
     434         efionn2(inter,2)=1.-efdisn2(inter)
     435      enddo
     436     
     437      !3 channels for CO ionization
     438      do inter=11,16
     439         efionco(inter,1)=1.-efdisco(inter)
     440      enddo
     441      efionco(10,1)=0.87*(1.-efdisco(10))
     442      efionco(10,2)=0.13*(1.-efdisco(10))
     443      do inter=8,9
     444         efionco(inter,2)=1.-efdisco(inter)
     445      enddo
     446      efionco(7,2)=0.1*(1.-efdisco(7))
     447      efionco(7,3)=0.9*(1.-efdisco(7))
     448      do inter=2,6
     449         efionco(inter,3)=1.-efdisco(inter)
     450      enddo
     451
     452      !Total ionization under 85 nm for N
     453      do inter=1,16
     454         efionn(inter)=1.
     455      enddo
     456
     457      !NO
     458      do inter=2,28
     459         efionno(inter)=1.-efdisno(inter)
     460      enddo
     461
     462      !Total ionization under 90 nm for H
     463      do inter=3,16
     464         efionh(inter)=1.
     465      enddo
    228466
    229467
    … …  
    232470
    233471      end
     472

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

    r476 r635  
    1212!
    1313!     Franck Lefevre
    14 !     version 1.1 november 2011
     14!     version 1.2 april 2012
    1515!==========================================================================
    1616
    … …  
    4646! local
    4747
    48       integer l, ig
    49       real rho                     ! density (kg/m3)
    50       real dustnd                  ! uodated dust number density (kg/kg)
    51       real icend                   ! uodated ice number density (kg/kg)
    52       real ccntyp                  ! typical dust number density (#/kg)
    53                                    ! (microphys = false)
    54       real rdusttyp                ! typical dust radius (m)
    55                                    ! (microphys = false)
     48      integer    :: l, ig
     49      real       :: rho                     ! density (kg/m3)
     50      real       :: dustnd, icend           ! uodated dust and ice number densities (kg/kg)
     51      real, save :: factor_ice, factor_dust ! multiplying factor to compute total surface area
     52                                            ! from the mass-mean radius
     53      real       :: sigma_ice, sigma_dust   ! variance of the ice and dust distributions
     54      real       :: ccntyp                  ! typical dust number density (#/kg)
     55                                            ! (microphys = false)
     56      real       :: rdusttyp                ! typical dust radius (m)
     57                                            ! (microphys = false)
     58
     59      logical, save :: firstcall = .true.
    5660
    5761!==========================================================================
     62
     63      if (firstcall) then ! compute the multiplying factors
     64         sigma_dust  = varian
     65         sigma_ice   = sqrt(log(nuice_sed + 1.))
     66         factor_dust = exp(0.5*(log(sigma_dust))**2)
     67         factor_ice  = exp(0.5*(log(sigma_ice))**2)
     68         write(*,*) 'surfacearea : factor_dust = ', factor_dust
     69         write(*,*) 'surfacearea : factor_ice  = ', factor_ice
     70         firstcall = .false.
     71      end if
    5872
    5973      if (microphys) then ! improvedclouds
    … …  
    6983     $                + pdq(ig,l,igcm_ccn_number)*ptimestep
    7084!              dust surface area
    71                surfdust(ig,l) = dustnd*rho*tauscaling(ig)
     85               surfdust(ig,l) = factor_dust*dustnd*rho*tauscaling(ig)
    7286     $                          *4.*pi*rdust(ig,l)**2
    7387!              ice surface area
    74                surfice(ig,l)  = icend*rho*tauscaling(ig)
     88               surfice(ig,l)  = factor_ice*icend*rho*tauscaling(ig)
    7589     $                          *4.*pi*rice(ig,l)**2
    7690            end do
    … …  
    88102               ccntyp = ccntyp/ccn_factor
    89103               if (rice(ig,l) .gt. rdust(ig,l)) then
    90                   surfdust(ig,l) = ccntyp*(ccn_factor - 1.)*rho
    91      $                             *4.*pi*rdusttyp**2
    92                   surfice(ig,l)  = ccntyp*4.*pi*rice(ig,l)**2
     104                  surfdust(ig,l) = factor_dust*ccntyp*(ccn_factor - 1.)
     105     $                             *rho*4.*pi*rdusttyp**2
     106                  surfice(ig,l)  = factor_ice*ccntyp*4.*pi*rice(ig,l)**2
    93107               else
    94                   surfdust(ig,l) = ccntyp*ccn_factor*rho
    95      $                             *4.*pi*rdusttyp**2
     108                  surfdust(ig,l) = factor_dust*ccntyp*ccn_factor
     109     $                             *rho*4.*pi*rdusttyp**2
    96110                  surfice(ig,l)  = 0.
    97111               end if
    … …  
    105119        call wstats(ngrid,"surfdust", "Dust surface area",
    106120     $            "micron2 cm-3",3,surfdust*1.e6)
    107       endif
    108       call writediagfi(ngrid,"surfdust", "Dust cloud surface area",
    109      $            "micron2 cm-3",3,surfdust*1.e6)
    110       if (callstats) then
    111121        call wstats(ngrid,"surfice", "Ice cloud surface area",
    112122     $            "micron2 cm-3",3,surfice*1.e6)
    113123      endif
     124      call writediagfi(ngrid,"surfdust", "Dust surface area",
     125     $            "micron2 cm-3",3,surfdust*1.e6)
    114126      call writediagfi(ngrid,"surfice", "Ice cloud surface area",
    115127     $            "micron2 cm-3",3,surfice*1.e6)

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

    r517 r635  
    1212#include "comdiurn.h"
    1313#include "param.h"
    14 #include "param_v3.h"
     14#include "param_v4.h"
    1515#include "chimiedata.h"
    1616#include "conc.h"
    … …  
    5858      if (calleuv) then
    5959        call zerophys(ngridmx*nlayermx,zdteuv)
    60         call euvheat(pt,pdt,pplev,pplay,zzlay,dist_sol,
     60        call euvheat(pt,pdt,pplev,pplay,zzlay,
    6161     $               mu0,ptimestep,ptime,zday,pq,pdq,zdteuv)
    6262      endif
    … …  
    7979      if (callmoldiff) THEN
    8080        call zerophys(ngridmx*nlayermx*nqmx,zdqmoldiff)
    81         call moldiff_red(pplay,pplev,pt,pdt,pq,pdq,ptimestep,
     81        call moldiff(pplay,pplev,pt,pdt,pq,pdq,ptimestep,
    8282     &                   zzlay,zdteuv,zdtconduc,zdqmoldiff)
    8383      endif

  • trunk/LMDZ.MARS/libf/phymars/callkeys.h

    r552 r635  
    1616     
    1717      COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche    &
    18      &   ,dustbin,nqchem_min,nltemodel,nircorr
     18     &   ,dustbin,nltemodel,nircorr
    1919     
    2020      COMMON/callkeys_r/topdustref,solarcondate,semi,alphan,euveff,     &
    … …  
    5454      logical sedimentation,activice,water,microphys,caps
    5555      logical photochem
    56       integer nqchem_min
    5756      integer nltemodel
    5857      integer nircorr

  • trunk/LMDZ.MARS/libf/phymars/callsedim.F

    r626 r635  
    434434c     Update the dust particle size "rdust"
    435435c     -------------------------------------
    436       DO l = 1, nlay
     436      if (doubleq) then
     437       DO l = 1, nlay
    437438        DO ig=1,ngrid
    438439          rdust(ig,l)=
    … …  
    441442          rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
    442443        ENDDO
    443       ENDDO
     444       ENDDO
     445      endif ! of if (doubleq)
    444446     
    445447c     Update the ice particle size "rice"
    446448c     -------------------------------------
    447       IF(scavenging) THEN
     449      if (water) then
     450       IF(scavenging) THEN
    448451        DO l = 1, nlay
    449452          DO ig=1,ngrid
    … …  
    464467          ENDDO
    465468        ENDDO
    466       ELSE
     469       ELSE
    467470        DO l = 1, nlay
    468471          DO ig=1,ngrid
    … …  
    475478          ENDDO
    476479        ENDDO
    477       ENDIF
     480       ENDIF ! of IF(scavenging)
     481      endif ! of if (water)
    478482
    479483      RETURN

  • trunk/LMDZ.MARS/libf/phymars/initracer.F

    r459 r635  
    2424c                  q(iq=2) is the dust number mixing ratio
    2525
    26 c   If (photochem.or.thermochem) there is "ncomp" chemical species (ncomp
    27 c   is set in aeronomars/chimiedata.h) using the ncomp iq values starting at
    28 c      iq=nqchem_min = dustbin+1   (nqchem_min is defined in inifis.F)
    29 c   
    3026c
    3127c   author: F.Forget
    … …  
    4339#include "advtrac.h"
    4440#include "comgeomfi.h"
    45 #include "chimiedata.h"
    4641
    4742#include "surfdat.h"
    … …  
    5550c       et al. (2002), a value of 25 has been deduced;
    5651      real, parameter :: popratio = 25.
    57       logical :: oldnames ! =.true. if old tracer naming convention (q01,...)
    5852      character(len=20) :: txt ! to store some text
    5953
    … …  
    7165c-----------------------------------------------------------------------
    7266
    73       integer :: nqchem_start
    74 
    7567! Initialization: get tracer names from the dynamics and check if we are
    7668!                 using 'old' tracer convention ('q01',q02',...)
    7769!                 or new convention (full tracer names)
    7870      ! check if tracers have 'old' names
    79       oldnames=.false.
    8071      count=0
    8172      do iq=1,nqmx
    … …  
    9081        write(*,*) "initracer: tracers seem to follow old naming ",
    9182     &             "convention (q01,q02,...)"
    92         write(*,*) "   => will work for now ... "
    93         write(*,*) "      but you should run newstart to rename them"
    94         oldnames=.true.
     83        write(*,*) "you should run newstart to rename them"
     84        stop
    9585      endif
    9686
    97       ! copy/set tracer names
    98       if (oldnames) then ! old names (derived from iq & option values)
    99         ! 1. dust:
    100         if (dustbin.ne.0) then ! transported dust
    101           do iq=1,dustbin
    102             txt=" "
    103             write(txt,'(a4,i2.2)') 'dust',iq
    104             noms(iq)=txt
    105             mmol(iq)=100.
    106           enddo
    107           ! special case if doubleq
    108           if (doubleq) then
    109             if (dustbin.ne.2) then
    110               write(*,*) 'initracer: error expected dustbin=2'
    111             else
    112 !              noms(1)='dust_mass'   ! dust mass mixing ratio
    113 !              noms(2)='dust_number' ! dust number mixing ratio
    114 ! same as above, but avoid explict possible out-of-bounds on array
    115                noms(1)='dust_mass'   ! dust mass mixing ratio
    116                do iq=2,2
    117                noms(iq)='dust_number' ! dust number mixing ratio
    118                enddo
    119             endif
    120           endif
    121         endif
    122         ! 2. water & ice
    123         if (water) then
    124           noms(nqmx)='h2o_vap'
    125           mmol(nqmx)=18.
    126 !            noms(nqmx-1)='h2o_ice'
    127 !            mmol(nqmx-1)=18.
    128           !"loop" to avoid potential out-of-bounds in array
    129           do iq=nqmx-1,nqmx-1
    130             noms(iq)='h2o_ice'
    131             mmol(iq)=18.
    132           enddo
    133         endif
    134         ! 3. Chemistry
    135         if (photochem .or. callthermos) then
    136           if (doubleq) then
    137             nqchem_start=3
    138           else
    139             nqchem_start=dustbin+1
    140           end if
    141           do iq=nqchem_start,nqchem_start+ncomp-2
    142             noms(iq)=nomchem(iq-nqchem_start+1)
    143             mmol(iq)=mmolchem(iq-nqchem_start+1)
    144           enddo
    145         endif ! of if (photochem .or. callthermos)
    146         ! 4. Other tracers ????
    147         if ((dustbin.eq.0).and.(.not.water)) then
    148           noms(1)='co2'
    149           mmol(1)=44
    150           if (nqmx.eq.2) then
    151             noms(nqmx)='Ar_N2'
    152             mmol(nqmx)=30
    153           endif
    154         endif
    155       else
    156         ! copy tracer names from dynamics
    157         do iq=1,nqmx
    158           noms(iq)=tnom(iq)
    159         enddo
    160         ! mmol(:) array is initialized later (see below)
    161       endif ! of if (oldnames)
    162 
    163       ! special modification when using 'old' tracers:
    164       ! qsurf(nqmx) was h2o ice whereas q(nqmx) was water vapour
    165       ! and (if iceparty) q(nqmx-1) was null whereas q(nqmx-1) was water ice
    166       if (oldnames.and.water) then
    167         write(*,*)'initracer: moving surface water ice to index ',nqmx-1
    168         ! "loop" to avoid potential out-of-bounds on arrays
    169         do iq=nqmx-1,nqmx-1
    170           qsurf(1:ngridmx,iq)=qsurf(1:ngridmx,iq+1)
    171         enddo
    172         qsurf(1:ngridmx,nqmx)=0
    173       endif
    174      
    175 c------------------------------------------------------------
    176 c     Test Dimensions tracers
    177 c------------------------------------------------------------
    178 
    179 ! Ehouarn: testing number of tracers is obsolete...
    180 !      if(photochem.or.thermochem) then
    181 !          if (water) then
    182 !              if ((dustbin+ncomp+2).ne.nqmx) then
    183 !                 print*,'initracer: tracer dimension problem:'
    184 !                 print*,"(dustbin+ncomp+2).ne.nqmx"
    185 !                 print*,"ncomp: ",ncomp
    186 !                 print*,"dustbin: ",dustbin
    187 !                 print*,"nqmx: ",nqmx
    188 !                 print*,'Change ncomp in chimiedata.h'
    189 !               endif
    190 !          else
    191 !              if ((dustbin+ncomp+1).ne.nqmx) then
    192 !                 print*,'initracer: tracer dimension problem:'
    193 !                 print*,"(dustbin+ncomp+1).ne.nqmx"
    194 !                 print*,"ncomp: ",ncomp
    195 !                 print*,"dustbin: ",dustbin
    196 !                 print*,"nqmx: ",nqmx
    197 !                 print*,'Change ncomp in chimiedata.h'
    198 !                 STOP
    199 !               endif
    200 !            endif
    201 !      endif
     87      ! copy tracer names from dynamics
     88      do iq=1,nqmx
     89        noms(iq)=tnom(iq)
     90      enddo
    20291
    20392c------------------------------------------------------------
    20493c         NAME and molar mass of the tracer
    20594c------------------------------------------------------------
    206 
    20795   
    20896! Identify tracers by their names: (and set corresponding values of mmol)
    … …  
    246134      igcm_n2plus=0
    247135      igcm_hplus=0
     136      igcm_hco2plus=0
    248137      igcm_elec=0
    249 
    250138
    251139      ! 1. find dust tracers
    … …  
    432320          count=count+1
    433321        endif
     322        if (noms(iq).eq."hco2plus") then
     323          igcm_hco2plus=iq
     324          mmol(igcm_hco2plus)=45.
     325          count=count+1
     326        endif
    434327        if (noms(iq).eq."elec") then
    435328          igcm_elec=iq
    … …  
    454347        endif
    455348
    456 
    457349      enddo ! of do iq=1,nqmx
    458 !      count=count+nbqchem
    459350     
    460351      ! check that we identified all tracers:

  • trunk/LMDZ.MARS/libf/phymars/physiq.F

    r633 r635  
    131131#include "chimiedata.h"
    132132#include "param.h"
    133 #include "param_v3.h"
     133#include "param_v4.h"
    134134#include "conc.h"
    135135

  • trunk/LMDZ.MARS/libf/phymars/tracer.h

    r420 r635  
    6666      integer :: igcm_n2plus
    6767      integer :: igcm_hplus
     68      integer :: igcm_hco2plus
    6869      integer :: igcm_elec
    6970      ! other tracers
    … …  
    8586     & igcm_co2plus,igcm_oplus,igcm_o2plus,igcm_coplus,igcm_cplus,      &
    8687     & igcm_nplus,igcm_noplus,igcm_n2plus,igcm_hplus,igcm_elec,         &
    87      & igcm_ar_n2!,nbqchem,niqchem
     88     & igcm_hco2plus,igcm_ar_n2!,nbqchem,niqchem
    8889      COMMON/tracer3/noms
    8990!-----------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.