Changeset 2542 for trunk/LMDZ.GENERIC/libf/aeronostd

Timestamp:

Jul 5, 2021, 2:44:34 PM (3 years ago)

Author:

aslmd

Message:

Generic GCM:

Large update of the chemical modules

• Read chemical network from input files
• Init chemistry with ModernTrac?
• Photolysis online calculation

YJ

Location:

trunk/LMDZ.GENERIC/libf/aeronostd

Files:

• calchim_asis.F90 (modified) (17 diffs)
• chimiedata.h (deleted)
• chimiedata_h.F90 (added)
• concentrations.F (modified) (5 diffs)
• deposition_source.F90 (added)
• pfunc.F90 (added)
• photochemistry_asis.F90 (modified) (33 diffs)
• photolysis_asis.F90 (modified) (14 diffs)
• photolysis_mod.F90 (added)
• photolysis_online.F (added)
• read_phototable.F90 (modified) (5 diffs)
• types_asis.F90 (modified) (2 diffs)

trunk/LMDZ.GENERIC/libf/aeronostd/calchim_asis.F90

@@ -3 +3 @@
                          zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,            &
                          tauref,co2ice,                                     &
-                         pu,pdu,pv,pdv,surfdust,surfice)
-
-      use tracer_h, only:   igcm_co2, igcm_co, igcm_o, igcm_o1d, igcm_o2, &
-                            igcm_o3, igcm_h, igcm_h2, igcm_oh, igcm_ho2,  &
-                            igcm_h2o2, igcm_ch4, igcm_n2, igcm_h2o_vap,   &
-                            igcm_n, igcm_no, igcm_no2, igcm_n2d,   &
-                            mmol
+                         pu,pdu,pv,pdv,surfdust,surfice,icount,zdtchim)
+
+      use tracer_h, only: mmol, noms, nesp, is_chim
 
       use conc_mod, only: mmean ! mean molecular mass of the atmosphere
       USE comcstfi_mod
       use callkeys_mod
+      use time_phylmdz_mod, only: ecritphy, iphysiq ! output rate set by ecritphy
+      use types_asis,       only: v_phot_3d, jlabel, nb_phot_hv_max
+      use chimiedata_h
+      use radcommon_h,      only: glat
 
       implicit none
@@ -35 +35 @@
 !    update 11/12/2013 optimization (Franck Lefevre)
 !    update 20/10/2017 adapted to LMDZ GENERIC+cosmetic changes (Benjamin Charnay)
+!    update 06/03/2021 generic tracer/network + photolysis online (Yassin Jaziri)
 !
 !   Arguments:
@@ -62 +63 @@
 !  Output:
 !
-!    dqchim(ngrid,nlayer,nq) ! tendencies on pq due to chemistry
-!    dqschim(ngrid,nq)         ! tendencies on qsurf 
-!
-!=======================================================================
-
-#include "chimiedata.h"
+!    dqchim(ngrid,nlayer,nesp)   ! tendencies on pq due to chemistry
+!    dqschim(ngrid,nesp)         ! tendencies on qsurf 
+!
+!=======================================================================
 
 !     input:
@@ -76 +75 @@
       real :: ptimestep
       real :: pplay(ngrid,nlayer)    ! pressure at the middle of the layers
-      real :: zzlay(ngrid,nlayer)    ! pressure at the middle of the layers
+      real :: zzlay(ngrid,nlayer)    ! altitude at the middle of the layers
       real :: pplev(ngrid,nlayer+1)  ! intermediate pressure levels
       real :: zzlev(ngrid,nlayer+1)  ! altitude at layer boundaries
@@ -97 +96 @@
 !     output:
 
-      real :: dqchim(ngrid,nlayer,nq)   ! tendencies on pq due to chemistry
-      real :: dqschim(ngrid,nq)         ! tendencies on qsurf 
+      real :: dqchim(ngrid,nlayer,nesp)   ! tendencies on pq due to chemistry
+      real :: dqschim(ngrid,nesp)         ! tendencies on qsurf 
 
 !     local variables:
 
-      integer,save :: nbq                   ! number of tracers used in the chemistry
-      integer,allocatable,save :: niq(:)    ! array storing the indexes of the tracers
-      integer :: iloc(1)                    ! index of major species
-      integer :: ig,l,i,iq,iqmax
-      integer :: foundswitch, lswitch
+      integer :: iloc(1)                  ! index of major species
+      integer :: ig,l,i,iq,iqmax,iesp
+      integer :: foundswitch, lswitch     ! to switch between photochem and thermochem ? (YJ)
       integer,save :: chemthermod
 
-      integer,save :: i_co2  = 0
-      integer,save :: i_co   = 0
-      integer,save :: i_o    = 0
-      integer,save :: i_o1d  = 0
-      integer,save :: i_o2   = 0
-      integer,save :: i_o3   = 0
-      integer,save :: i_h    = 0
-      integer,save :: i_h2   = 0
-      integer,save :: i_oh   = 0
-      integer,save :: i_ho2  = 0
-      integer,save :: i_h2o2 = 0
-      integer,save :: i_ch4  = 0
-      integer,save :: i_n2   = 0
-      integer,save :: i_h2o  = 0
-      integer,save :: i_n    = 0
-      integer,save :: i_no    = 0
-      integer,save :: i_no2    = 0
-      integer,save :: i_n2d    = 0
-
-      integer :: ig_vl1
-
-      real    :: latvl1, lonvl1
-      real    :: zq(ngrid,nlayer,nq)   ! pq+pdq*ptimestep before chemistry
+      real    :: zq(ngrid,nlayer,nesp) ! pq+pdq*ptimestep before chemistry
                                        ! new mole fraction after
       real    :: zt(ngrid,nlayer)      ! temperature
@@ -137 +112 @@
       real    :: zv(ngrid,nlayer)      ! v component of the wind
       real    :: taucol                ! optical depth at 7 hPa
-
+      real    :: xmmol(nlayer,nesp)    ! mmol/mmean but only for chemical species
+ 
       logical,save :: firstcall = .true.
-      logical,save :: depos = .false.  ! switch for dry deposition
 
 !     for each column of atmosphere:
@@ -147 +122 @@
       real :: ztemp(nlayer)        ! Temperature (K)
       real :: zlocal(nlayer)       ! Altitude (km)
-      real :: zycol(nlayer,nq)     ! Composition (mole fractions)
+      real :: zycol(nlayer,nesp)   ! Composition (mole fractions)
       real :: zmmean(nlayer)       ! Mean molecular mass (g.mole-1)
       real :: szacol               ! Solar zenith angle
@@ -154 +129 @@
       real :: surfice1d(nlayer)    ! Ice surface area (cm2/cm3)
       real :: surfdust1d(nlayer)   ! Dust surface area (cm2/cm3)
-      real :: jo3(nlayer)          ! Photodissociation rate O3->O1D (s-1)
 
       integer :: iter(nlayer)      !  Number of chemical iterations
                                    !  within one physical timestep 
-
+      integer :: icount
 !     for output:
 
-      logical :: output             ! to issue calls to writediagfi and stats
+      logical :: output              ! to issue calls to writediagfi and stats
       parameter (output = .true.)
-      real :: jo3_3d(ngrid,nlayer)  ! Photodissociation rate O3->O1D (s-1)
-      real :: iter_3d(ngrid,nlayer) ! Number of chemical iterations
-                                    !  within one physical timestep
+      real :: iter_3d(ngrid,nlayer)  ! Number of chemical iterations
+                                     !  within one physical timestep
+      integer           ::  ierr
+      real              ::  zdtchim(ngrid,nlayer)    ! temperature modification by chemistry
+      real              ::  dEzchim(ngrid,nlayer)    ! energy modification by chemistry
+      real              ::  zdtchim_output(ngrid)    ! flux modification by chemistry in W.m-2
 
 !=======================================================================
@@ -173 +150 @@
       if (firstcall) then
 
-         if (photochem) then
-            print*,'calchim: Read photolysis lookup table'
-            call read_phototable
-         end if
-         ! find index of chemical tracers to use
-         allocate(niq(nq))
-         ! Listed here are all tracers that can go into photochemistry
-         nbq = 0 ! to count number of tracers
-         ! Species ALWAYS present if photochem=.T. or thermochem=.T.
-         i_co2 = igcm_co2
-         if (i_co2 == 0) then
-            write(*,*) "calchim: Error; no CO2 tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_co2
-         end if
-         i_co = igcm_co
-         if (i_co == 0) then
-            write(*,*) "calchim: Error; no CO tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_co
-         end if
-         i_o = igcm_o
-         if (i_o == 0) then
-            write(*,*) "calchim: Error; no O tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_o
-         end if
-         i_o1d = igcm_o1d
-         if (i_o1d == 0) then
-            write(*,*) "calchim: Error; no O1D tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_o1d
-         end if
-         i_o2 = igcm_o2
-         if (i_o2 == 0) then
-            write(*,*) "calchim: Error; no O2 tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_o2
-         end if
-         i_o3 = igcm_o3
-         if (i_o3 == 0) then
-            write(*,*) "calchim: Error; no O3 tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_o3
-         end if
-         i_h = igcm_h
-         if (i_h == 0) then
-            write(*,*) "calchim: Error; no H tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_h
-         end if
-         i_h2 = igcm_h2
-         if (i_h2 == 0) then
-            write(*,*) "calchim: Error; no H2 tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_h2
-         end if
-         i_oh = igcm_oh
-         if (i_oh == 0) then
-            write(*,*) "calchim: Error; no OH tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_oh
-         end if
-         i_ho2 = igcm_ho2
-         if (i_ho2 == 0) then
-            write(*,*) "calchim: Error; no HO2 tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_ho2
-         end if
-         i_h2o2 = igcm_h2o2
-         if (i_h2o2 == 0) then
-            write(*,*) "calchim: Error; no H2O2 tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_h2o2
-         end if
-         i_ch4 = igcm_ch4
-         if (i_ch4 == 0) then
-            write(*,*) "calchim: Error; no CH4 tracer !!!"
-            write(*,*) "CH4 will be ignored in the chemistry"
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_ch4
-         end if
-         i_n2 = igcm_n2
-         if (i_n2 == 0) then
-            write(*,*) "calchim: Error; no N2 tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_n2
-         end if
-         i_n = igcm_n
-         if (i_n == 0) then
-            if (photochem) then
-               write(*,*) "calchim: Error; no N tracer !!!"
-               write(*,*) "N will be ignored in the chemistry"
+!! Moved to routine indice in photochemistry_asis
+!! because nb_phot_hv_max value needed in order
+!! to choose if we call read_phototable or not.
+!! A cleaner solution need to be find.
+!         if (photochem .and. .not. jonline) then
+!           print*,'calchim: Read photolysis lookup table'
+!           call read_phototable
+!         end if
+
+         if (.not.allocated(SF_mode))       allocate(SF_mode(nesp))
+         if (.not.allocated(SF_value))      allocate(SF_value(nesp))
+         if (.not.allocated(prod_rate))     allocate(prod_rate(nesp))
+         if (.not.allocated(surface_flux))  allocate(surface_flux(ngrid,nesp))
+         if (.not.allocated(surface_flux2)) allocate(surface_flux2(ngrid,nesp))
+         if (.not.allocated(escape))        allocate(escape(ngrid,nesp))
+         if (.not.allocated(chemnoms))      allocate(chemnoms(nesp))
+
+         surface_flux(:,:)  = 0.
+         surface_flux2(:,:) = 0.
+         escape(:,:)        = 0.
+         SF_mode(:)         = 2
+         SF_value(:)        = 0.
+         prod_rate(:)       = 0.
+         iter_3d(:,:)       = 0.
+         iter(:)            = 0.
+
+         call ini_tracchim
+
+         ! Sanity check mmol /= 0. in chemistry
+         do iq = 1,nq
+            if (is_chim(iq).eq.1 .and. mmol(iq).eq.0.) then
+               write(*,*) 'Error in calchim:'
+               write(*,*) 'Mmol cannot be equal to 0 for chemical species'
+               stop
             end if
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_n
-         end if
-         i_n2d = igcm_n2d
-         if (i_n2d == 0) then
-            if (photochem) then
-               write(*,*) "calchim: Error; no N2D tracer !!!"
-               write(*,*) "N2d will be ignored in the chemistry"
-            end if
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_n2d
-         end if
-         i_no = igcm_no
-         if (i_no == 0) then
-            if (photochem) then
-               write(*,*) "calchim: Error; no NO tracer !!!"
-               write(*,*) "NO will be ignored in the chemistry"
-            end if
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_no
-         end if
-         i_no2 = igcm_no2
-         if (i_no2 == 0) then
-            if (photochem) then
-               write(*,*) "calchim: Error; no NO2 tracer !!!"
-               write(*,*) "NO2 will be ignored in the chemistry"
-            end if
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_no2
-         end if
-         i_h2o = igcm_h2o_vap
-         if (i_h2o == 0) then
-            write(*,*) "calchim: Error; no water vapor tracer !!!"
-            stop
-         else
-            nbq = nbq + 1
-            niq(nbq) = i_h2o
-         end if
-         !Check tracers needed for thermospheric chemistry
-!         if(thermochem) then
-!            chemthermod=0  !Default: C/O/H chemistry
-!            !Nitrogen chemistry
-!            !NO is used to determine if N chemistry is wanted
-!            !chemthermod=2 -> N chemistry
-!            if (i_no == 0) then
-!               write(*,*) "calchim: no NO tracer"
-!               write(*,*) "C/O/H themosp chemistry only "
-!            else
-!               chemthermod=2
-!               write(*,*) "calchim: NO in traceur.def"
-!               write(*,*) "Nitrogen chemistry included"
-!            end if
-!            ! N
-!            if(chemthermod == 2) then
-!               if (i_n == 0) then
-!                  write(*,*) "calchim: Error; no N tracer !!!"
-!                  write(*,*) "N is needed if NO is in traceur.def"
-!                  stop
-!               end if
-!            ! NO2
-!               if (i_no2 == 0) then
-!                  write(*,*) "calchim: Error; no NO2 tracer !!!"
-!                  write(*,*) "NO2 is needed if NO is in traceur.def"
-!                  stop
-!               end if
-!            ! N(2D)
-!               if (i_n2d == 0) then
-!                  write(*,*) "calchim: Error; no N2D !!!"
-!                  write(*,*) "N2D is needed if NO is in traceur.def"
-!                  stop
-!               end if
-!            endif    !Of if(chemthermod == 2)
-!         endif      !Of thermochem
-
-         write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
-               
+         end do
+
          firstcall = .false.
       end if ! if (firstcall)
@@ -380 +195 @@
       dqchim(:,:,:) = 0.
       dqschim(:,:)  = 0.
-
-!     latvl1= 22.27
-!     lonvl1= -47.94
-!     ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.)  -2 )*iim +    &
-!             int(1.5+(lonvl1+180)*iim/360.)
 
 !=======================================================================
@@ -394 +204 @@
          foundswitch = 0
          do l = 1,nlayer
-            do i = 1,nbq
-               iq = niq(i) ! get tracer index
-               zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
-               zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
-            end do
+            iesp = 0
+            do iq = 1,nq
+               if (is_chim(iq).eq.1) then
+                  iesp = iesp + 1
+                  zq(ig,l,iesp) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
+                  xmmol(l,iesp) = mmol(iq)/mmean(ig,l)
+                  zycol(l,iesp) = zq(ig,l,iesp)/xmmol(l,iesp)
+               end if
+            end do
+
             zt(ig,l)  = pt(ig,l) + pdt(ig,l)*ptimestep
             zu(ig,l)  = pu(ig,l) + pdu(ig,l)*ptimestep
@@ -430 +245 @@
          end do ! of do l=1,nlayer
 
-         szacol = acos(mu0(ig))*180./pi
-         taucol = tauref(ig)*(700./610.)  ! provisoire en attente de nouveau jmars
-         fractcol=fract(ig)
+         szacol   = acos(mu0(ig))*180./pi
+         taucol   = tauref(ig)*(700./610.)  ! provisoire en attente de nouveau jmars
+         fractcol = fract(ig)
 
 !=======================================================================
@@ -440 +255 @@
 !        chemistry in lower atmosphere
 
-         if (photochem) then
-
-            call photochemistry_asis(nlayer,nq,ngrid,                          &
-                                ig,lswitch,zycol,szacol,fractcol,ptimestep,    &
-                                zpress,ztemp,zdens,zmmean,dist_sol,            &
-                                surfdust1d,surfice1d,jo3,taucol,iter)
-
-!        ozone photolysis, for output
-
+!         if (photochem) then
+
+            call photochemistry_asis(nlayer,ngrid,                         &
+                                ig,lswitch,zycol,szacol,fractcol,ptimestep,&
+                                zpress,zlocal,ztemp,zdens,zmmean,dist_sol, &
+                                surfdust1d,surfice1d,taucol,iter,zdtchim(ig,:))
+
+            ! diagnostic photochemical heating
+            zdtchim_output(ig) = 0.
             do l = 1,nlayer
-               jo3_3d(ig,l) = jo3(l)
+              dEzchim(ig,l) = zdtchim(ig,l)*cpp*(pplev(ig,l) - pplev(ig,l+1))/glat(ig)
+              zdtchim_output(ig) = zdtchim_output(ig) + zdtchim(ig,l)*cpp*(pplev(ig,l) - pplev(ig,l+1))/glat(ig)
+            end do
+
+            do l = 1,nlayer
                iter_3d(ig,l) = iter(l)
             end do
-
-!        condensation of h2o2
-
-!            call perosat(ngrid, nlayer, nq,                        &
-!                         ig,ptimestep,pplev,pplay,                 &
-!                         ztemp,zycol,dqcloud,dqscloud)
-         end if
 
 !        chemistry in upper atmosphere
@@ -470 +282 @@
 !        dry deposition
 
-!         if (depos) then
-!            call deposition(ngrid, nlayer, nq,                      &
-!                            ig, ig_vl1, pplay, pplev, zzlay, zzlev, & 
-!                            zu, zv, zt, zycol, ptimestep, co2ice)
-!         end if
+         if (depos) then
+            call deposition_source(ngrid, nlayer, nq,   &
+                        ig, zzlay, zzlev, zdens, zycol, ptimestep)
+            surface_flux2(ig,:) = surface_flux2(ig,:) + surface_flux(ig,:)
+            if (ngrid==1) then
+              if(mod(icount,ecritphy).eq.0) then
+                surface_flux2(ig,:) = surface_flux2(ig,:)/float(ecritphy)
+              endif
+            else
+              if(mod(icount*iphysiq,ecritphy).eq.0) then
+                surface_flux2(ig,:) = surface_flux2(ig,:)*iphysiq/float(ecritphy)
+              endif
+            endif
+         end if
 
 !=======================================================================
@@ -489 +310 @@
             iloc=maxloc(zycol(l,:))
             iqmax=iloc(1)
-            do i = 1,nbq
-               iq = niq(i) ! get tracer index
+            do iq = 1,nesp
                if (iq /= iqmax) then
-                  dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l) - zq(ig,l,iq))/ptimestep
+                  dqchim(ig,l,iq) = (zycol(l,iq)*xmmol(l,iq) - zq(ig,l,iq))/ptimestep
                   dqchim(ig,l,iq) = max(dqchim(ig,l,iq),-zq(ig,l,iq)/ptimestep)
                   dqchim(ig,l,iqmax) = dqchim(ig,l,iqmax) - dqchim(ig,l,iq) 
                end if
             end do
+
          end do ! of do l = 1,nlayer
 
@@ -508 +329 @@
 !=======================================================================
 
-      end do ! of do ig=1,ngrid
+      end do ! of do ig=1,ngridbidon(ig,:) = 1
 
 !=======================================================================
@@ -516 +337 @@
 ! value of parameter 'output' to trigger writting of outputs
 ! is set above at the declaration of the variable.
-
       if (photochem .and. output) then
-            call writediagfi(ngrid,'jo3','j o3->o1d',    &
-                             's-1',3,jo3_3d(1,1))
-            call writediagfi(ngrid,'iter','iterations',  &
+
+         if (callstats) then
+            ! photoloysis
+            do i=1,nb_phot_hv_max
+               call wstats(ngrid,trim(jlabel(i,1)),jlabel(i,1),    &
+                                's-1',3,v_phot_3d(1,1,i))
+            end do
+            ! iter
+            call wstats(ngrid,'iter','iterations',  &
                              ' ',3,iter_3d(1,1))
-            if (callstats) then
-               call wstats(ngrid,'jo3','j o3->o1d',       &
-                           's-1',3,jo3_3d(1,1))
-               call wstats(ngrid,'mmean','mean molecular mass',       &
+            ! phothcemical heating
+            call wstats(ngrid,'zdtchim','dT chim',    &
+                             'K.s-1',3,zdtchim)
+            call wstats(ngrid,'dEzchim','dE chim',    &
+                             'W m-2',3,dEzchim)
+            call wstats(ngrid,"Ezchim","integrated dT chim","W m-2",2,zdtchim_output)
+            ! Mean molecular mass
+            call wstats(ngrid,'mmean','mean molecular mass',       &
                            'g.mole-1',3,mmean(1,1))
+            ! Chemical tendencies
+            do iesp=1,nesp
+               call wstats(ngrid,trim(chemnoms(iesp))//'_dq',trim(chemnoms(iesp))//'_dq',    & 
+                               'kg/kg s^-1',3,dqchim(1,1,iesp) )
+            end do
+            if (depos) then
+               ! Surface fluxes
+               do iesp=1,nesp
+                  call wstats(ngrid,trim(chemnoms(iesp))//'_flux_mean',trim(chemnoms(iesp))//' mean flux',    & 
+                               'molecule.m-2.s-1',2,surface_flux2(1,indexchim(trim(chemnoms(iesp)))))
+                  call wstats(ngrid,trim(chemnoms(iesp))//'_flux',trim(chemnoms(iesp))//' flux',    & 
+                               'molecule.m-2.s-1',2,surface_flux(1,indexchim(trim(chemnoms(iesp)))))
+               end do
+            endif ! end depos
+         endif ! end callstats
+
+         ! photoloysis
+         do i=1,nb_phot_hv_max
+            call writediagfi(ngrid,trim(jlabel(i,1)),jlabel(i,1),    &
+                             's-1',3,v_phot_3d(1,1,i))
+         end do
+         ! iter
+         call writediagfi(ngrid,'iter','iterations',  &
+                          ' ',3,iter_3d(1,1))
+         ! phothcemical heating
+         call writediagfi(ngrid,'zdtchim','dT chim',    &
+                          'K.s-1',3,zdtchim)
+         call writediagfi(ngrid,'dEzchim','dE chim',    &
+                          'W m-2',3,dEzchim)
+         call writediagfi(ngrid,"Ezchim","integrated dT chim","W m-2",2,zdtchim_output)
+         ! Mean molecular mass
+         call writediagfi(ngrid,'mmean','mean molecular mass',       &
+                        'g.mole-1',3,mmean(1,1))
+         ! Chemical tendencies
+         do iesp=1,nesp
+            call writediagfi(ngrid,trim(chemnoms(iesp))//'_dq',trim(chemnoms(iesp))//'_dq',    & 
+                            'kg/kg s^-1',3,dqchim(1,1,iesp) )
+         end do
+         if (depos) then
+            ! Surface fluxes
+            do iesp=1,nesp
+               call writediagfi(ngrid,trim(chemnoms(iesp))//'_flux_mean',trim(chemnoms(iesp))//' mean flux',    & 
+                            'molecule.m-2.s-1',2,surface_flux2(1,indexchim(trim(chemnoms(iesp)))))
+               call writediagfi(ngrid,trim(chemnoms(iesp))//'_flux',trim(chemnoms(iesp))//' flux',    & 
+                            'molecule.m-2.s-1',2,surface_flux(1,indexchim(trim(chemnoms(iesp)))))
+            end do
+            ! Restart flux average
+            if (ngrid == 1) then
+               if(mod(icount,ecritphy).eq.0) then
+                  surface_flux2(:,:) = 0.0
+               endif
+            else
+               if(mod(icount*iphysiq,ecritphy).eq.0) then
+                  surface_flux2(:,:) = 0.0
+               endif
             endif
+         endif ! end depos
+
       end if ! of if (output)
 
       return
-      end
+
+
+      contains
+
+
+!======================================================================
+
+      subroutine ini_tracchim
+
+!======================================================================
+! YJ Modern tracer.def
+! Get in tracer.def chemical parameters
+!======================================================================
+
+         use chimiedata_h
+         use tracer_h, only:   is_chim
+         implicit none
+
+         ! local
+         character(len=500) :: tracline ! to store a line of text
+         integer            :: nq       ! number of tracers
+         integer            :: iesp, iq
+
+         ! Get nq
+         open(90,file='traceur.def',status='old',form='formatted',iostat=ierr)
+         if (ierr.eq.0) then
+           READ(90,'(A)') tracline
+           IF (trim(tracline).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
+            READ(tracline,*,iostat=ierr) nq ! Try standard traceur.def
+           ELSE
+              DO
+                READ(90,'(A)',iostat=ierr) tracline
+                IF (ierr.eq.0) THEN
+                  IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
+                    READ(tracline,*,iostat=ierr) nq
+                    EXIT
+                  ENDIF
+                ELSE ! If pb, or if reached EOF without having found number of tracer
+                   write(*,*) "calchim: error reading line of tracers"
+                   write(*,*) "   (first lines of traceur.def) "
+                   stop
+                ENDIF
+              ENDDO
+           ENDIF ! if modern or standard traceur.def
+         else
+            write(*,*) "calchim: error opening traceur.def"
+            stop
+         endif
+
+         ! Get data of tracers
+         iesp = 0
+         do iq=1,nq
+            read(90,'(A)') tracline
+            if (is_chim(iq).eq.1) then
+               iesp = iesp + 1
+               read(tracline,*) chemnoms(iesp)
+               write(*,*)"ini_tracchim: iq=",iq,"noms(iq)=",trim(chemnoms(iesp))
+               if (index(tracline,'SF_mode='    ) /= 0) then
+                  read(tracline(index(tracline,'SF_mode=')+len('SF_mode='):),*) SF_mode(iesp)
+                  write(*,*) ' Parameter value (traceur.def) : SF_mode=', SF_mode(iesp)
+               else
+                  write(*,*) ' Parameter value (default)     : SF_mode=', SF_mode(iesp)
+               end if
+               if (index(tracline,'SF_value='    ) /= 0) then
+                  read(tracline(index(tracline,'SF_value=')+len('SF_value='):),*) SF_value(iesp)
+                  write(*,*) ' Parameter value (traceur.def) : SF_value=', SF_value(iesp)
+               else
+                  write(*,*) ' Parameter value (default)     : SF_value=', SF_value(iesp)
+               end if
+               if (index(tracline,'prod_rate='    ) /= 0) then
+                  read(tracline(index(tracline,'prod_rate=')+len('prod_rate='):),*) prod_rate(iesp)
+                  write(*,*) ' Parameter value (traceur.def) : prod_rate=', prod_rate(iesp)
+               else
+                  write(*,*) ' Parameter value (default)     : prod_rate=', prod_rate(iesp)
+               end if
+            end if
+         end do
+
+         close(90)
+
+      end subroutine ini_tracchim
+
+      end subroutine calchim_asis
+

trunk/LMDZ.GENERIC/libf/aeronostd/concentrations.F

@@ -2 +2 @@
      &                          pplay,pt,pdt,pq,pdq,ptimestep)
                                              
-      use tracer_h, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,
-     &                      igcm_o2, igcm_o3, igcm_h, igcm_h2,
-     &                      igcm_oh, igcm_ho2, igcm_n2, igcm_ar,
-     &                      igcm_h2o_vap, igcm_n, igcm_no, igcm_no2,
-     &                      igcm_n2d, igcm_ch4, 
-     &                      igcm_ch3, igcm_ch, igcm_3ch2, igcm_1ch2,      
-     &                      igcm_cho, igcm_ch2o, igcm_ch3o,               
-     &                      igcm_c, igcm_c2, igcm_c2h, igcm_c2h2,         
-     &                      igcm_c2h3, igcm_c2h4, igcm_c2h6, igcm_ch2co,  
-     &                      igcm_ch3co, igcm_hcaer,
-     &                      igcm_h2o2, mmol
+      use tracer_h, only: mmol, noms, aki, cpi, nesp
 
-      use conc_mod, only: mmean, Akknew, rnew, cpnew
+      use conc_mod, only: mmean, akknew, rnew, cpnew
       USE comcstfi_mod                   
       use callkeys_mod
+      use chimiedata_h
       implicit none
 
@@ -28 +19 @@
 !
 ! version: April 2012 - Franck Lefevre
+! update 06/03/2021 cpi/aki input (Yassin Jaziri)
 !=======================================================================
 
-!     declarations
- 
-#include "chimiedata.h"
 !     input/output
 
-      integer,intent(in) :: ngrid ! number of atmospheric columns
+      integer,intent(in) :: ngrid  ! number of atmospheric columns
       integer,intent(in) :: nlayer ! number of atmospheric layers
-      integer,intent(in) :: nq ! number of tracers
-      real,intent(in) :: pplay(ngrid,nlayer)
-      real,intent(in) :: pt(ngrid,nlayer)
-      real,intent(in) :: pdt(ngrid,nlayer)
-      real,intent(in) :: pq(ngrid,nlayer,nq)
-      real,intent(in) :: pdq(ngrid,nlayer,nq)
-      real,intent(in) :: ptimestep
+      integer,intent(in) :: nq     ! number of tracers
+      real,   intent(in) :: pplay(ngrid,nlayer)
+      real,   intent(in) :: pt(ngrid,nlayer)
+      real,   intent(in) :: pdt(ngrid,nlayer)
+      real,   intent(in) :: pq(ngrid,nlayer,nq)
+      real,   intent(in) :: pdq(ngrid,nlayer,nq)
+      real,   intent(in) :: ptimestep
 
 !     local variables
 
-      integer       :: i, l, ig, iq
-      integer, save :: nbq
-      integer,allocatable,save :: niq(:)
+      integer       :: l, ig, iq
       real          :: ni(nq), ntot
       real          :: zq(ngrid, nlayer, nq)
       real          :: zt(ngrid, nlayer)
-      real,allocatable,save    :: aki(:)
-      real,allocatable,save    :: cpi(:)
-
-      logical, save :: firstcall = .true.
-
-
-      if (firstcall) then
-
-         ! allocate local saved arrays:
-         allocate(aki(nq))
-         allocate(cpi(nq))
-         allocate(niq(nq))
-!        find index of chemical tracers to use
-!        initialize thermal conductivity and specific heat coefficients
-!        !? values are estimated
-
-         nbq = 0 ! to count number of tracers used in this subroutine
-
-         if (igcm_co2 /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_co2
-            aki(nbq) = 3.072e-4
-            cpi(nbq) = 0.834e3
-         end if
-         if (igcm_co /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_co
-            aki(nbq) = 4.87e-4
-            cpi(nbq) = 1.034e3
-         end if
-         if (igcm_o /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_o
-            aki(nbq) = 7.59e-4
-            cpi(nbq) = 1.3e3
-         end if
-         if (igcm_o1d /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_o1d
-            aki(nbq) = 7.59e-4  !?
-            cpi(nbq) = 1.3e3    !?
-         end if
-         if (igcm_o2 /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_o2
-            aki(nbq) = 5.68e-4
-            cpi(nbq) = 0.9194e3
-         end if
-         if (igcm_o3 /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_o3
-            aki(nbq) = 3.00e-4  !?
-            cpi(nbq) = 0.800e3  !?
-         end if
-         if (igcm_h /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_h
-            aki(nbq) = 0.0
-            cpi(nbq) = 20.780e3
-         end if
-         if (igcm_h2 /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_h2
-            aki(nbq) = 36.314e-4
-            cpi(nbq) = 14.266e3
-         end if
-         if (igcm_oh /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_oh
-            aki(nbq)  = 7.00e-4 !?
-            cpi(nbq)  = 1.045e3
-         end if
-         if (igcm_ho2 /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_ho2
-            aki(nbq) = 0.0
-            cpi(nbq) = 1.065e3  !?
-         end if
-         if (igcm_h2o2 /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_h2o2
-            aki(nbq) = 0.0
-            cpi(nbq) = 1.065e3  !?
-         end if
-         if (igcm_h2o_vap /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_h2o_vap
-            aki(nbq) = 0.0
-            cpi(nbq) = 1.870e3
-         end if
-         if (igcm_n /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_n
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_n2d /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_n2d
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_no /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_no
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_no2 /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_no2
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if (igcm_n2 /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_n2
-            aki(nbq) = 5.6e-4
-            cpi(nbq) = 1.034e3
-         end if
-         if(igcm_ch4 /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_ch4
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_ch3 /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_ch3
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_ch /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_ch
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_1ch2 /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_1ch2
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_3ch2 /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_3ch2
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_cho /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_cho
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_ch2o /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_ch2o
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_ch3o /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_ch3o
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_c /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_c
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_c2 /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_c2
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_c2h /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_c2h
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_c2h2 /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_c2h2
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_c2h3 /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_c2h3
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_c2h4 /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_c2h4
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_c2h6 /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_c2h6
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_ch2co /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_ch2co
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_ch3co /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_ch3co
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if(igcm_hcaer /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_hcaer
-            aki(nbq) = 0.0
-            cpi(nbq) = 0.0
-         endif
-         if (igcm_ar /= 0) then
-            nbq = nbq + 1
-            niq(nbq) = igcm_ar
-            aki(nbq) = 0.0      !?
-            cpi(nbq) = 1.000e3  !?
-         end if
-
-
-         ! tell the world about it:
-         write(*,*) "concentrations: firstcall, nbq=",nbq
-         write(*,*) "  niq(1:nbq)=",niq(1:nbq)
-         write(*,*) "  aki(1:nbq)=",aki(1:nbq)
-         write(*,*) "  cpi(1:nbq)=",cpi(1:nbq)
-
-         firstcall = .false.
-
-      end if ! if (firstcall)
 
 !     update temperature
@@ -311 +53 @@
       do l = 1,nlayer
          do ig = 1,ngrid
-            do i = 1,nbq
-               iq = niq(i) 
+            do iq = 1,nq
                zq(ig,l,iq) = max(1.e-30, pq(ig,l,iq)
      $                                 + pdq(ig,l,iq)*ptimestep)
@@ -325 +66 @@
       do l = 1,nlayer
          do ig = 1,ngrid
-            do i = 1,nbq
-               iq = niq(i)
-               mmean(ig,l) = mmean(ig,l) + zq(ig,l,iq)/mmol(iq)
+            do iq = 1,nq
+               if (mmol(iq).ne.0.) then
+                  mmean(ig,l) = mmean(ig,l) + zq(ig,l,iq)/mmol(iq)
+               end if
             end do
             mmean(ig,l) = 1./mmean(ig,l)
@@ -343 +85 @@
          do ig = 1,ngrid
             ntot = pplay(ig,l)/(1.381e-23*zt(ig,l))*1.e-6  ! in #/cm3
-            do i = 1,nbq
-               iq = niq(i) 
+            do iq = 1,nq
                ni(iq) = ntot*zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
-               cpnew(ig,l) = cpnew(ig,l) + ni(iq)*cpi(i)
-               akknew(ig,l) = akknew(ig,l) + ni(iq)*aki(i)
+               cpnew(ig,l)  = cpnew(ig,l) + ni(iq)*cpi(iq)
+               akknew(ig,l) = akknew(ig,l) + ni(iq)*aki(iq)
             end do 
-            cpnew(ig,l) = cpnew(ig,l)/ntot
+            cpnew(ig,l)  = cpnew(ig,l)/ntot
             akknew(ig,l) = akknew(ig,l)/ntot
          end do

trunk/LMDZ.GENERIC/libf/aeronostd/photochemistry_asis.F90

@@ -4 +4 @@
 !
 !     Author: Franck Lefevre
+!             Benjamin Charnay
+!             Yassin Jaziri
 !     ------
 !
-!     Version: 10/11/2015
+!     Version: 27/05/2016
+!     update 06/03/2021 generic tracer/network + photolysis online (Yassin Jaziri)
 !
 !*****************************************************************
 
-subroutine photochemistry_asis(nlayer, nq, ngrid,                               &
+subroutine photochemistry_asis(nlayer, ngrid,                                   &
                           ig, lswitch, zycol, sza, fractcol, ptimestep, press,  &
-                          temp, dens, zmmean, dist_sol, surfdust1d,             &
-                          surfice1d, jo3, tau, iter)
+                          alt, temp, dens, zmmean, dist_sol, surfdust1d,        &
+                          surfice1d, tau, iter,zdtchim)
 
       use callkeys_mod
+      use comcstfi_mod, only:  r,cpp,avocado,pi
+      use tracer_h
+      use types_asis
+      use chimiedata_h
+      use photolysis_mod
+
 implicit none
 
-#include "chimiedata.h" 
-
 !===================================================================
 !     inputs:
@@ -25 +32 @@
 
 integer, intent(in) :: nlayer ! number of atmospheric layers
-integer, intent(in) :: nq     ! number of tracers
-integer,intent(in) :: ngrid   ! number of atmospheric columns
+integer, intent(in) :: ngrid  ! number of atmospheric columns
+
 integer :: ig                 ! grid point index
       
@@ -33 +40 @@
 real :: ptimestep             ! physics timestep (s)
 real :: press(nlayer)         ! pressure (hpa)
+real :: alt(nlayer)           ! altitude (km)
 real :: temp(nlayer)          ! temperature (k)
 real :: dens(nlayer)          ! density (cm-3)
@@ -45 +53 @@
 !===================================================================
       
-real :: zycol(nlayer,nq)       ! chemical species volume mixing ratio
+real :: zycol(nlayer,nesp)    ! chemical species volume mixing ratio
 
 !===================================================================
@@ -51 +59 @@
 !===================================================================
       
-integer :: iter(nlayer)        ! iteration counter
-real    :: jo3(nlayer)         ! photodissociation rate o3 -> o1d
+integer :: iter(nlayer)       ! iteration counter
+real    :: zdtchim(nlayer)    ! temperature modification by ozone
 
 !===================================================================
@@ -59 +67 @@
 
 integer :: phychemrat         ! (physical timestep)/(nominal chemical timestep)
-integer :: j_o3_o1d, ilev
-integer :: iesp, nesp
+integer :: ilev, iesp, iphot, iw
 integer :: lswitch
 
 logical, save :: firstcall = .true.
-
-parameter (nesp = 17)         ! number of species in the chemistry
-
-! tracer indexes in the chemistry:
-
-integer,parameter :: i_co2  =  1
-integer,parameter :: i_co   =  2
-integer,parameter :: i_o    =  3
-integer,parameter :: i_o1d  =  4
-integer,parameter :: i_o2   =  5
-integer,parameter :: i_o3   =  6
-integer,parameter :: i_h    =  7
-integer,parameter :: i_h2   =  8
-integer,parameter :: i_oh   =  9
-integer,parameter :: i_ho2  = 10
-integer,parameter :: i_h2o2 = 11
-integer,parameter :: i_h2o  = 12
-integer,parameter :: i_n    = 13
-integer,parameter :: i_n2d  = 14
-integer,parameter :: i_no   = 15
-integer,parameter :: i_no2  = 16
-integer,parameter :: i_n2   = 17
 
 real :: stimestep           ! standard timestep for the chemistry (s) 
@@ -95 +80 @@
 real :: j(nlayer,nd)        ! interpolated photolysis rates (s-1)
 real :: time                ! internal time (between 0 and ptimestep, in s)
+real :: rho(nlayer)         ! mass density (kg/m-3)
 
 real, dimension(nlayer,nesp)            :: rm   ! mixing ratios 
@@ -103 +89 @@
 ! reaction rates
  
-real (kind = 8), dimension(nlayer,      nb_phot_max) :: v_phot
-real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
-real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
-logical :: hetero_dust, hetero_ice
+real (kind = 8), allocatable, save     :: v_phot(:,:)
+real (kind = 8), allocatable, save     :: v_3(:,:)
+real (kind = 8), allocatable, save     :: v_4(:,:)
+real (kind = 8), allocatable, save     :: e_phot(:,:)    ! photolysis rates by energie (J.mol-1.s-1)
+
+integer, save                          :: nreact         ! number of reactions in reactions files
+integer, allocatable, save             :: rtype(:)       ! reaction rate type
+integer, allocatable, save             :: third_body(:)  ! if the reaction have a third body: index of the third body, else zero
+logical, allocatable, save             :: three_prod(:)  ! if the reaction have a three defferent proucts egal .true.
 
 ! matrix
 
-real (kind = 8), dimension(nesp,nesp) :: mat, mat1
-integer, dimension(nesp)              :: indx
-integer                               :: code
+real (kind = 8), dimension(nesp,nesp)  :: mat, mat1
+integer, dimension(nesp)               :: indx
+integer                                :: code
 
 ! production and loss terms (for first-guess solution only)
 
-real (kind = 8), dimension(nesp) :: prod, loss
-
-! curvatures
-
-real :: ratio, curv, e, e1, e2, e3
+real (kind = 8), dimension(nesp)       :: prod, loss
 
 !===================================================================
@@ -128 +115 @@
 if (firstcall) then
    print*,'photochemistry: initialize indexes'
-   call indice(i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
-               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
-               i_n, i_n2d, i_no, i_no2, i_n2)
+   call indice(nreact,rtype,third_body,three_prod)
+   allocate(v_phot(nlayer,nb_phot_max))
+   allocate(v_3(nlayer,nb_reaction_3_max))
+   allocate(v_4(nlayer,nb_reaction_4_max))
+   allocate(v_phot_3d(ngrid,nlayer,nb_phot_hv_max))
+   allocate(e_phot(nlayer,nb_phot_max))
+   v_phot(:,:)      = 0.
+   v_3(:,:)         = 0.
+   v_4(:,:)         = 0.
+   v_phot_3d(:,:,:) = 0.
+   e_phot(:,:)      = 0.
+
+   ! Need to be done after subroutine indice
+   if (jonline) then
+     print*,'calchim: Read UV absorption cross-sections'
+     ! Jonline cannot run without photolysis
+     if (nb_phot_hv_max == 0) then
+       print*,'Jonline cannot run without photolysis'
+       print*,'set jonline=.false. in callphys.def'
+       print*,'or set photolysis reactions in monoreact file'
+       stop
+     end if
+     call init_photolysis(nlayer,nreact,three_prod)     ! on-line photolysis
+     allocate(albedo_snow_chim(nw))
+     allocate(albedo_co2_ice_chim(nw))
+    ! Step 1 : Initialisation of the Spectral Albedos.
+     DO iw=1,nw
+      albedo_snow_chim(iw)=albedosnow
+      albedo_co2_ice_chim(iw)=albedoco2ice
+
+      ! Spectral Snow Albedo Calculation.
+      call albedo_snow_calc(wc(iw)*1.0e-3,        &
+                     albedo_snow_chim(iw),  & 
+                     albedosnow)
+      
+     ENDDO
+   end if
+
    firstcall = .false.
 end if
@@ -138 +160 @@
 !===================================================================
 
-call gcmtochim(nlayer, nq, zycol, lswitch, nesp,               &
-               i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
-               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
-               i_n, i_n2d, i_no, i_no2, i_n2, dens, rm, c) 
+call gcmtochim(nlayer, zycol, lswitch, nesp, dens, rm, c) 
 
 !===================================================================
@@ -147 +166 @@
 !===================================================================
 
-call photolysis_asis(nlayer, ngrid, lswitch, press, temp, sza, fractcol, tau, zmmean, dist_sol, & 
-                     rm(:,i_co2), rm(:,i_o3), v_phot)
-
-! save o3 photolysis for output
-
-j_o3_o1d = 5
-jo3(:) = v_phot(:,j_o3_o1d)
+if (jonline) then
+   if (sza <= 95.) then
+      call photolysis_online(nlayer, alt, press, temp, zmmean, rm,   &
+                             tau, sza, dist_sol, v_phot, e_phot, ig, ngrid, nreact, three_prod)
+      if (ngrid.eq.1) then
+        do iphot = 1,nb_phot_hv_max
+          v_phot(:,iphot) = v_phot(:,iphot)* 0.25 / cos(sza*pi/180.) ! globally averaged = divide by 4
+          e_phot(:,iphot) = e_phot(:,iphot)* 0.25 / cos(sza*pi/180.) ! globally averaged = divide by 4
+        end do
+      elseif(diurnal .eqv. .false.) then
+        do iphot = 1,nb_phot_hv_max
+          v_phot(:,iphot) = v_phot(:,iphot) * fractcol
+          e_phot(:,iphot) = e_phot(:,iphot) * fractcol
+        end do
+      endif
+   else ! night
+      v_phot(:,:) = 0.
+      e_phot(:,:) = 0.
+   end if
+   ! save photolysis for output
+   v_phot_3d(ig,:,:) = v_phot(:,1:nb_phot_hv_max)
+else if (nb_phot_hv_max /= 0) then
+   call photolysis_asis(nlayer, ngrid, lswitch, press, temp, sza, fractcol,tau, zmmean, dist_sol, & 
+                     rm(:,indexchim('co2')), rm(:,indexchim('o3')), rm(:,indexchim('ch4')), v_phot)
+   ! save photolysis for output
+   v_phot_3d(ig,:,:) = v_phot(:,1:nb_phot_hv_max)
+end if
 
 !===================================================================
 !     reaction rates                                     
 !===================================================================
-!     switches for heterogeneous chemistry
-!     hetero_ice  : reactions on ice clouds
-!     hetero_dust : reactions on dust    
-!===================================================================
-
-hetero_dust = .false.
-hetero_ice  = .false.
-
-call reactionrates(nlayer, lswitch, dens, c(:,i_co2), c(:,i_o2), &
-                   press, temp, hetero_dust, hetero_ice,         &
-                   surfdust1d, surfice1d, v_phot, v_3, v_4)
+
+call reactionrates(nlayer, nreact, rtype, third_body, three_prod, c, lswitch, dens, &
+                   press, temp, zmmean, sza, surfdust1d, surfice1d, v_phot, v_3, v_4)
+
+zdtchim(:) = 0.
+rho(:)     = (press(:)*1e2)/(r*temp(:))
+
+!===================================================================
+!     temperature modification by photolysis reaction                                    
+!===================================================================
+
+if (photoheat .and. jonline) then ! for now just with jonline
+
+  do iphot = 1,nb_phot_hv_max
+    zdtchim(:nlayer-1) = zdtchim(:nlayer-1) + e_phot(:nlayer-1,iphot)*c(:nlayer-1,indexchim(trim(jlabel(iphot,2))))/(cpp*(rho(:nlayer-1)*1e-6)*avocado)
+  end do
+
+else
+
+!! o + o2 -> o3
+!!dE = 24 ! kcal.mol-1
+!!zdtchim(:) = zdtchim(:) + 4.184*24e3*v_4(:,1)*c(:,indexchim('o'))*c(:,indexchim('o2'))*1e6/(cpp*rho*avocado)
+!
+!! o3 -> o2 + o1d
+!! E(250nm) = 480 kJ.mol-1
+!j_o3_o1d   = 3
+!zdtchim(:) = zdtchim(:) + 480e3*v_phot(:,j_o3_o1d)*c(:,indexchim('o3'))/(cpp*(rho*1e-6)*avocado)
+!
+!! o3 -> o2 + o
+!! E(350nm) = 343 kJ.mol-1
+!j_o3_o     = 4
+!zdtchim(:) = zdtchim(:) + 343e3*v_phot(:,j_o3_o)*c(:,indexchim('o3'))/(cpp*(rho*1e-6)*avocado)
+
+zdtchim(:) = 0.
+
+end if
 
 !===================================================================
@@ -182 +246 @@
 
 ctimestep = ptimestep/real(phychemrat)
-
-!print*, "stimestep  = ", stimestep
-!print*, "ptimestep  = ", ptimestep
-!print*, "phychemrat = ", phychemrat
-!print*, "ctimestep  = ", ctimestep
-!stop
 
 !===================================================================
@@ -207 +265 @@
 
    iter(ilev) = iter(ilev) + 1    ! iteration counter
-
+  
 !  first-guess: fill matrix
 
@@ -268 +326 @@
 !===================================================================
 
-call chimtogcm(nlayer, nq, zycol, lswitch, nesp,              &
-               i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,      &
-               i_h2, i_oh, i_ho2, i_h2o2, i_h2o,              &
-               i_n, i_n2d, i_no, i_no2, i_n2, dens, c) 
+call chimtogcm(nlayer, zycol, lswitch, nesp, dens, c) 
 
 contains
@@ -318 +373 @@
 real (kind = 8), parameter :: dtmin   = 10.      ! minimum time step (s)
 real (kind = 8), parameter :: vmrtol  = 1.e-11   ! absolute tolerance on vmr
-real (kind = 8), parameter :: rtol    = 1./0.05   ! 1/rtol recommended value : 0.1-0.02
+real (kind = 8), parameter :: rtol    = 1./0.05  ! 1/rtol recommended value : 0.1-0.02
 integer,         parameter :: niter   = 3        ! number of iterations
 real (kind = 8), parameter :: coefmax = 2.
@@ -382 +437 @@
 ! timestep correction
 
-coef = max(coefmin, min(coefmax,0.8/sqrt(e)))
+if (e <= 0.) then
+  coef = coefmax
+else
+  coef = max(coefmin, min(coefmax,0.8/sqrt(e)))
+end if
 
 dttest = max(dtmin,dttest*coef)
@@ -396 +455 @@
 
 
-
-
 !======================================================================
 
- subroutine reactionrates(nlayer,                               &
-                          lswitch, dens, co2, o2, press, t,     &
-                          hetero_dust, hetero_ice,              &
+ subroutine reactionrates(nlayer, nreact, rtype, third_body, three_prod, c, &
+                          lswitch, dens, press, t, zmmean, sza, &
                           surfdust1d, surfice1d,                &
                           v_phot, v_3, v_4)
@@ -418 +474 @@
 
 use comcstfi_mod
+use types_asis
+use pfunc
+use tracer_h
+use chimiedata_h
 
 implicit none
-
-#include "chimiedata.h"
 
 !----------------------------------------------------------------------
@@ -427 +485 @@
 !----------------------------------------------------------------------
 
-integer, intent(in)     :: nlayer            ! number of atmospheric layers
-integer                 :: lswitch           ! interface level between lower
-                                             ! atmosphere and thermosphere chemistries
-real, dimension(nlayer) :: dens              ! total number density (molecule.cm-3)
-real, dimension(nlayer) :: press             ! pressure (hPa)
-real, dimension(nlayer) :: t                 ! temperature (K)
-real, dimension(nlayer) :: surfdust1d        ! dust surface area (cm2.cm-3)
-real, dimension(nlayer) :: surfice1d         ! ice surface area (cm2.cm-3)
-real (kind = 8), dimension(nlayer) :: co2    ! co2 number density (molecule.cm-3)
-real (kind = 8), dimension(nlayer) :: o2     ! o2 number density (molecule.cm-3)
-logical :: hetero_dust, hetero_ice           ! switches for heterogeneous chemistry
+integer, intent(in)     :: nlayer             ! number of atmospheric layers
+integer, intent(in)     :: nreact             ! number of reactions in reactions files
+integer                 :: lswitch            ! interface level between lower
+                                              ! atmosphere and thermosphere chemistries
+real, dimension(nlayer) :: dens               ! total number density (molecule.cm-3)
+real, dimension(nlayer) :: press              ! pressure (hPa)
+real, dimension(nlayer) :: t                  ! temperature (K)
+real, dimension(nlayer) :: zmmean             ! mean molar mass (g/mole)
+real                    :: sza                ! solar zenith angle (deg)
+real, dimension(nlayer) :: surfdust1d         ! dust surface area (cm2.cm-3)
+real, dimension(nlayer) :: surfice1d          ! ice surface area (cm2.cm-3)
+real (kind = 8), dimension(nlayer,nesp) :: c  ! species number density (molecule.cm-3)
+
+integer, intent(in)     :: rtype(nreact)      ! reaction rate type
+integer, intent(in)     :: third_body(nreact) ! if the reaction have a third body: index of the third body, else zero
+logical, intent(in)     :: three_prod(nreact) ! if the reaction have three different products egal .true.
 
 !----------------------------------------------------------------------
@@ -451 +514 @@
 !----------------------------------------------------------------------
 
-integer          :: ilev
+logical,save     :: firstcall = .true.
+integer          :: ilev, ireact
 integer          :: nb_phot, nb_reaction_3, nb_reaction_4
-real :: ak0, ak1, xpo, rate
-real :: k1a0, k1b0, k1ainf, k1a, k1b, fc, fx, x, y, gam
-real, dimension(nlayer) :: deq
-real, dimension(nlayer) :: a001, a002, a003,                           &
-                             b001, b002, b003, b004, b005, b006, b007,   &
-                             b008, b009,                                 &
-                             c001, c002, c003, c004, c005, c006, c007,   &
-                             c008, c009, c010, c011, c012, c013, c014,   &
-                             c015, c016, c017, c018,                     &
-                             d001, d002, d003, d004, d005, d006, d007,   &
-                             d008, d009,                                 &
-                             e001, e002,                                 &
-                             h001, h002, h003, h004, h005
+integer          :: nb_hv, nb_pfunc1, nb_pfunc2, nb_pfunc3
+real, dimension(nlayer) :: ratec ! rate coefficient
 
 !----------------------------------------------------------------------
@@ -471 +524 @@
 !----------------------------------------------------------------------
 
-      nb_phot       = 11       ! jmars.20140930 reduit de 13 a 11
+      nb_phot       = nb_phot_hv_max
       nb_reaction_3 = 0
       nb_reaction_4 = 0
 
+      nb_hv         = 0
+      nb_pfunc1     = 0
+      nb_pfunc2     = 0
+      nb_pfunc3     = 0
+
 !----------------------------------------------------------------------
-!     oxygen reactions 
+!     reactions 
 !----------------------------------------------------------------------
 
-!---  a001: o + o2 + co2 -> o3 + co2
-
-!     jpl 2003
-!
-!     co2 efficiency as a third body (2.075)
-!     from sehested et al., j. geophys. res., 100, 1995.
-
-      a001(:) = 2.075*6.0e-34*(t(:)/300.)**(-2.4)*dens(:)
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = a001(:)
-
-!---  a002: o + o + co2 -> o2 + co2
-
-!     Tsang and Hampson, J. Chem. Phys. Ref. Data, 15, 1087, 1986
-
-!     a002(:) = 2.5*5.2e-35*exp(900./t(:))*dens(:)
-
-!     Campbell and Gray, Chem. Phys. Lett., 18, 607, 1973
-
-!     a002(:) = 1.2e-32*(300./t(:))**(2.0)*dens(:)  ! yung expression
-      a002(:) = 2.5*9.46e-34*exp(485./t(:))*dens(:) ! nist expression
-
-      nb_reaction_3 = nb_reaction_3 + 1
-      v_3(:,nb_reaction_3) = a002(:)
-
-!---  a003: o + o3 -> o2 + o2
-
-!     jpl 2003
-
-      a003(:) = 8.0e-12*exp(-2060./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = a003(:)
-
-!----------------------------------------------------------------------
-!     o(1d) reactions
-!----------------------------------------------------------------------
-
-!---  b001: o(1d) + co2  -> o + co2
-
-!     jpl 2006
-
-      b001(:) = 7.5e-11*exp(115./t(:))
-   
+ireact = 1
+
+! Reaction from monoreact file are implemented first
+do while(nb_phot<nb_phot_max-nphot_hard_coding)
+
+  if (rtype(ireact)/=0) then     ! skip photolysis reactions
+    ! get right coefficient rate function
+    if (rtype(ireact)==1) then
+      nb_pfunc1 = nb_pfunc1 + 1
+      if (third_body(ireact)==0) then            !! third_body check
+        ratec = pfunc1(nlayer,t,dens,pfunc1_param(nb_pfunc1))
+      else
+        ratec = pfunc1(nlayer,t,c(:,third_body(ireact)),pfunc1_param(nb_pfunc1))
+      end if
+    else if (rtype(ireact)==2) then
+      nb_pfunc2 = nb_pfunc2 + 1
+      if (third_body(ireact)==0) then            !! third_body check
+        ratec = pfunc2(nlayer,t,dens,pfunc2_param(nb_pfunc2))
+      else
+        ratec = pfunc2(nlayer,t,c(:,third_body(ireact)),pfunc2_param(nb_pfunc2))
+      end if
+    else if (rtype(ireact)==3) then
+      nb_pfunc3 = nb_pfunc3 + 1
+      if (third_body(ireact)==0) then            !! third_body check
+        ratec = pfunc3(nlayer,t,dens,pfunc3_param(nb_pfunc3))
+      else
+        ratec = pfunc3(nlayer,t,c(:,third_body(ireact)),pfunc3_param(nb_pfunc3))
+      end if
+    else
+      print*, 'Error in reactionrates: wrong index coefficient rate vphot'
+      print*, 'rtype(ireact) = ',rtype(ireact)
+      print*, 'It should be 1 or 2 '
+      print*, 'Please check the reaction ',ireact
+      if (ireact>nreact) print*, 'Please check the reaction count into the code beacause ireact > nreact is not possible'
+      stop
+    end if
+
+    ! fill the right type index
+    nb_phot = nb_phot + 1
+    v_phot(:,nb_phot) = ratec(:)
+    if (three_prod(ireact)) then             ! three_prod check
       nb_phot = nb_phot + 1
-      v_phot(:,nb_phot) = b001(:)*co2(:)
-
-!---  b002: o(1d) + h2o  -> oh + oh
-
-!     jpl 2006
- 
-      b002(:) = 1.63e-10*exp(60./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = b002(:)
-
-!---  b003: o(1d) + h2  -> oh + h
-
-!     jpl 2011
-
-      b003(:) = 1.2e-10
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = b003(:)
-
-!---  b004: o(1d) + o2  -> o + o2
-
-!     jpl 2006
-
-      b004(:) = 3.3e-11*exp(55./t(:))
-
-      nb_phot = nb_phot + 1
-      v_phot(:,nb_phot) = b004(:)*o2(:)
-    
-!---  b005: o(1d) + o3  -> o2 + o2
-
-!     jpl 2003
-
-      b005(:) = 1.2e-10
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = b005(:)
-    
-!---  b006: o(1d) + o3  -> o2 + o + o
-
-!     jpl 2003
-
-      b006(:) = 1.2e-10
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = b006(:)
-    
-!---  b007: o(1d) + ch4 -> ch3 + oh
-
-!     jpl 2003
-
-      b007(:) = 1.5e-10*0.75
-
-!---  b008: o(1d) + ch4 -> ch3o + h
-
-!     jpl 2003
-
-      b008(:) = 1.5e-10*0.20
-!
-!---  b009: o(1d) + ch4 -> ch2o + h2
-
-!     jpl 2003
-
-      b009(:) = 1.5e-10*0.05
-
-!----------------------------------------------------------------------
-!     hydrogen reactions
-!----------------------------------------------------------------------
-
-!---  c001: o + ho2 -> oh + o2
-
-!     jpl 2003
-
-      c001(:) = 3.0e-11*exp(200./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c001(:)
-
-!---  c002: o + oh -> o2 + h
-
-!     jpl 2011
-
-      c002(:) = 1.8e-11*exp(180./t(:))
-
-!     robertson and smith, j. chem. phys. a 110, 6673, 2006
-
-!     c002(:) = 11.2e-11*t(:)**(-0.32)*exp(177./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c002(:)
-
-!---  c003: h + o3 -> oh + o2
-
-!     jpl 2003
-
-      c003(:) = 1.4e-10*exp(-470./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c003(:)
-
-!---  c004: h + ho2 -> oh + oh
-
-!     jpl 2006
-
-      c004(:) = 7.2e-11
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c004(:)
-
-!---  c005: h + ho2 -> h2 + o2
-
-!     jpl 2006
-
-      c005(:) = 6.9e-12
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c005(:)
-
-!---  c006: h + ho2 -> h2o + o
-
-!     jpl 2006
-
-      c006(:) = 1.6e-12
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c006(:)
-
-!---  c007: oh + ho2 -> h2o + o2
-
-!     jpl 2003
-
-!     canty et al., grl, 2006 suggest to increase this rate
-!     by 20%. not done here.
-
-      c007(:) = 4.8e-11*exp(250./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c007(:)
-
-!---  c008: ho2 + ho2 -> h2o2 + o2
-
-!     jpl 2006
-
-!     c008(:) = 3.5e-13*exp(430./t(:))
-
-!     christensen et al., grl, 13, 2002
-
-      c008(:) = 1.5e-12*exp(19./t(:))
-
-      nb_reaction_3 = nb_reaction_3 + 1
-      v_3(:,nb_reaction_3) = c008(:)
-
-!---  c009: oh + h2o2 -> h2o + ho2
-
-!     jpl 2006
-
-      c009(:) = 1.8e-12
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c009(:)
-
-!---  c010: oh + h2 -> h2o + h
-
-!     jpl 2006
-
-      c010(:) = 2.8e-12*exp(-1800./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c010(:)
-
-!---  c011: h + o2 + co2 -> ho2 + co2
-
-!     jpl 2011
-
-      do ilev = 1,lswitch-1
-         ak0 = 2.5*4.4e-32*(t(ilev)/300.)**(-1.3)
-         ak1 = 7.5e-11*(t(ilev)/300.)**(0.2)
-
-         rate = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1)
-         xpo = 1./(1. + alog10((ak0*dens(ilev))/ak1)**2)
-         c011(ilev) = rate*0.6**xpo
-      end do
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c011(:)
-
-!---  c012: o + h2o2 -> oh + ho2
-
-!     jpl 2003
-
-      c012(:) = 1.4e-12*exp(-2000./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c012(:)
-
-!---  c013: oh + oh -> h2o + o
-
-!     jpl 2006
-
-      c013(:) = 1.8e-12
-
-      nb_reaction_3 = nb_reaction_3 + 1
-      v_3(:,nb_reaction_3) = c013(:)
-
-!---  c014: oh + o3 -> ho2 + o2
-
-!     jpl 2003
-
-      c014(:) = 1.7e-12*exp(-940./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c014(:)
-
-!---  c015: ho2 + o3 -> oh + o2 + o2
-
-!     jpl 2003
-
-      c015(:) = 1.0e-14*exp(-490./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = c015(:)
-
-!---  c016: ho2 + ho2 + co2 -> h2o2 + o2 + co2
-
-!     jpl 2011
-
-      c016(:) = 2.5*2.1e-33*exp(920./t(:))*dens(:)
-
-      nb_reaction_3 = nb_reaction_3 + 1
-      v_3(:,nb_reaction_3) = c016(:)
-
-!---  c017: oh + oh + co2 -> h2o2 + co2
-
-!     jpl 2003
-
-      do ilev = 1,lswitch-1
-         ak0 = 2.5*6.9e-31*(t(ilev)/300.)**(-1.0)
-         ak1 = 2.6e-11*(t(ilev)/300.)**(0.0)
-
-         rate = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1)
-         xpo = 1./(1. + alog10((ak0*dens(ilev))/ak1)**2)
-         c017(ilev) = rate*0.6**xpo
-      end do
-
-      nb_reaction_3 = nb_reaction_3 + 1
-      v_3(:,nb_reaction_3) = c017(:)
-
-!---  c018: h + h + co2 -> h2 + co2
-
-!     baulch et al., 2005
-
-      c018(:) = 2.5*1.8e-30*(t(:)**(-1.0))*dens(:)
-
-      nb_reaction_3 = nb_reaction_3 + 1
-      v_3(:,nb_reaction_3) = c018(:)
-
-!----------------------------------------------------------------------
-!     nitrogen reactions
-!----------------------------------------------------------------------
-
-!---  d001: no2 + o -> no + o2
-
-!     jpl 2006
-
-      d001(:) = 5.1e-12*exp(210./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = d001(:)
-
-!---  d002: no + o3 -> no2 + o2
-
-!     jpl 2006
-
-      d002(:) = 3.0e-12*exp(-1500./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = d002(:)
-
-!---  d003: no + ho2 -> no2 + oh
-
-!     jpl 2011
-
-      d003(:) = 3.3e-12*exp(270./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = d003(:)
-
-!---  d004: n + no -> n2 + o
-
-!     jpl 2011
-
-      d004(:) = 2.1e-11*exp(100./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = d004(:)
-
-!---  d005: n + o2 -> no + o
-
-!     jpl 2011
-
-      d005(:) = 1.5e-11*exp(-3600./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = d005(:)
-
-!---  d006: no2 + h -> no + oh
-
-!     jpl 2011
-
-      d006(:) = 4.0e-10*exp(-340./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = d006(:)
-
-!---  d007: n + o -> no
- 
-      d007(:) = 2.8e-17*(300./t(:))**0.5
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = d007(:)
-
-!---  d008: n + ho2 -> no + oh
-
-!     brune et al., j. chem. phys., 87, 1983 
-
-      d008(:) = 2.19e-11
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = d008(:)
-
-!---  d009: n + oh -> no + h
-
-!     atkinson et al., j. phys. chem. ref. data, 18, 881, 1989 
-
-      d009(:) = 3.8e-11*exp(85./t(:))
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = d009(:)
-
-!----------------------------------------------------------------------
-!     carbon reactions
-!----------------------------------------------------------------------
-
-!---  e001: oh + co -> co2 + h
-
-!     jpl 2003
-
-!     e001(:) = 1.5e-13*(1 + 0.6*press(:)/1013.)
-
-!     mccabe et al., grl, 28, 3135, 2001
-
-!     e001(:) = 1.57e-13 + 3.54e-33*dens(:)
-
-!     jpl 2006
-
-!     ak0 = 1.5e-13*(t(:)/300.)**(0.6)
-!     ak1 = 2.1e-9*(t(:)/300.)**(6.1)
-!     rate1 = ak0/(1. + ak0/(ak1/dens(:)))
-!     xpo1 = 1./(1. + alog10(ak0/(ak1/dens(:)))**2)
-
-!     ak0 = 5.9e-33*(t(:)/300.)**(-1.4)
-!     ak1 = 1.1e-12*(t(:)/300.)**(1.3)
-!     rate2 = (ak0*dens(:))/(1. + ak0*dens(:)/ak1)
-!     xpo2 = 1./(1. + alog10((ak0*dens(:))/ak1)**2)
-
-!     e001(:) = rate1*0.6**xpo1 + rate2*0.6**xpo2
-
-!     joshi et al., 2006
-
-      do ilev = 1,lswitch-1
-         k1a0 = 1.34*2.5*dens(ilev)                                  &
-               *1/(1/(3.62e-26*t(ilev)**(-2.739)*exp(-20./t(ilev)))  &
-               + 1/(6.48e-33*t(ilev)**(0.14)*exp(-57./t(ilev))))     ! typo in paper corrected
-         k1b0 = 1.17e-19*t(ilev)**(2.053)*exp(139./t(ilev))          &
-              + 9.56e-12*t(ilev)**(-0.664)*exp(-167./t(ilev))
-         k1ainf = 1.52e-17*t(ilev)**(1.858)*exp(28.8/t(ilev))        &
-                + 4.78e-8*t(ilev)**(-1.851)*exp(-318./t(ilev))
-         x = k1a0/(k1ainf - k1b0)
-         y = k1b0/(k1ainf - k1b0)
-         fc = 0.628*exp(-1223./t(ilev)) + (1. - 0.628)*exp(-39./t(ilev))  &
-            + exp(-t(ilev)/255.)
-         fx = fc**(1./(1. + (alog(x))**2))                           ! typo in paper corrected
-         k1a = k1a0*((1. + y)/(1. + x))*fx
-         k1b = k1b0*(1./(1.+x))*fx
-
-         e001(ilev) = k1a + k1b
-      end do
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = e001(:)
-
-!---  e002: o + co + m -> co2 + m
-
-!     tsang and hampson, 1986.
-
-      e002(:) = 2.5*6.5e-33*exp(-2184./t(:))*dens(:)
-
-      nb_reaction_4 = nb_reaction_4 + 1
-      v_4(:,nb_reaction_4) = e002(:)
-
-!----------------------------------------------------------------------
-!     heterogeneous chemistry 
-!----------------------------------------------------------------------
-
-      if (hetero_ice) then
-
-!        k = (surface*v*gamma)/4 (s-1)
-!        v = 100*sqrt(8rt/(pi*m))  (cm s-1)
- 
-!---     h001: ho2 + ice -> products
- 
-!        cooper and abbatt, 1996: gamma = 0.025
-      
-         gam = 0.025
-         h001(:) = surfice1d(:)       &
-                   *100.*sqrt(8.*8.31*t(:)/(33.e-3*pi))*gam/4.
- 
-!        h002: oh + ice -> products
- 
-!        cooper and abbatt, 1996: gamma = 0.03
- 
-         gam = 0.03
-         h002(:) = surfice1d(:)       &
-                   *100.*sqrt(8.*8.31*t(:)/(17.e-3*pi))*gam/4.
-
-!---     h003: h2o2 + ice -> products
- 
-!        gamma = 0.    test value
- 
-         gam = 0.
-         h003(:) = surfice1d(:)       &
-                   *100.*sqrt(8.*8.31*t(:)/(34.e-3*pi))*gam/4.
-      else
-         h001(:) = 0.
-         h002(:) = 0.
-         h003(:) = 0.
-      end if
-
-      nb_phot = nb_phot + 1
-      v_phot(:,nb_phot) = h001(:)
-
-      nb_phot = nb_phot + 1
-      v_phot(:,nb_phot) = h002(:)
-
-      nb_phot = nb_phot + 1
-      v_phot(:,nb_phot) = h003(:)
-
-      if (hetero_dust) then
- 
-!---     h004: ho2 + dust -> products
- 
-!        jacob, 2000: gamma = 0.2
-!        see dereus et al., atm. chem. phys., 2005
- 
-         gam = 0.2
-         h004(:) = surfdust1d(:)  &
-                   *100.*sqrt(8.*8.31*t(:)/(33.e-3*pi))*gam/4.
- 
-!---     h005: h2o2 + dust -> products
- 
-!        gamma = 5.e-4
-!        see dereus et al., atm. chem. phys., 2005
- 
-         gam = 5.e-4
-         h005(:) = surfdust1d(:)  &
-                   *100.*sqrt(8.*8.31*t(:)/(34.e-3*pi))*gam/4.
-      else
-         h004(:) = 0.
-         h005(:) = 0.
-      end if
- 
-      nb_phot = nb_phot + 1
-      v_phot(:,nb_phot) = h004(:)
-
-      nb_phot = nb_phot + 1
-      v_phot(:,nb_phot) = h005(:)
+      v_phot(:,nb_phot) = ratec(:)
+    end if
+  else
+    nb_hv = nb_hv + 1
+  end if
+
+  ireact = ireact + 1
+
+end do
+
+! Reaction from bimolreact file are implemented secondly
+do while(nb_reaction_4<nb_reaction_4_max-n4_hard_coding)
+
+  ! get right coefficient rate function
+  if (rtype(ireact)==1) then
+    nb_pfunc1 = nb_pfunc1 + 1
+    if (third_body(ireact)==0) then            !! third_body check
+      ratec = pfunc1(nlayer,t,dens,pfunc1_param(nb_pfunc1))
+    else
+      ratec = pfunc1(nlayer,t,c(:,third_body(ireact)),pfunc1_param(nb_pfunc1))
+    end if
+  else if (rtype(ireact)==2) then
+    nb_pfunc2 = nb_pfunc2 + 1
+    if (third_body(ireact)==0) then            !! third_body check
+      ratec = pfunc2(nlayer,t,dens,pfunc2_param(nb_pfunc2))
+    else
+      ratec = pfunc2(nlayer,t,c(:,third_body(ireact)),pfunc2_param(nb_pfunc2))
+    end if
+  else if (rtype(ireact)==3) then
+    nb_pfunc3 = nb_pfunc3 + 1
+    if (third_body(ireact)==0) then            !! third_body check
+      ratec = pfunc3(nlayer,t,dens,pfunc3_param(nb_pfunc3))
+    else
+      ratec = pfunc3(nlayer,t,c(:,third_body(ireact)),pfunc3_param(nb_pfunc3))
+    end if
+  else
+    print*, 'Error in reactionrates: wrong index coefficient rate v4'
+    print*, 'rtype(ireact) = ',rtype(ireact)
+    print*, 'It should be 1 or 2 '
+    print*, 'Please check the reaction ',ireact-nb_phot_max
+    if (ireact>nreact) print*, 'Please check the reaction count into the code beacause ireact > nreact is not possible'
+    stop
+  end if
+
+  ! fill the right type index
+  nb_reaction_4 = nb_reaction_4 + 1
+  v_4(:,nb_reaction_4) = ratec(:)
+  if (three_prod(ireact)) then             ! three_prod check
+    nb_reaction_4 = nb_reaction_4 + 1
+    v_4(:,nb_reaction_4) = ratec(:)
+  end if
+
+  ireact = ireact + 1
+
+end do
+
+! Reaction from quadrareact file are implemented thirdly
+do while(nb_reaction_3<nb_reaction_3_max-n3_hard_coding)
+
+  ! get right coefficient rate function
+  if (rtype(ireact)==1) then
+    nb_pfunc1 = nb_pfunc1 + 1
+    if (third_body(ireact)==0) then            !! third_body check
+      ratec = pfunc1(nlayer,t,dens,pfunc1_param(nb_pfunc1))
+    else
+      ratec = pfunc1(nlayer,t,c(:,third_body(ireact)),pfunc1_param(nb_pfunc1))
+    end if
+  else if (rtype(ireact)==2) then
+    nb_pfunc2 = nb_pfunc2 + 1
+    if (third_body(ireact)==0) then            !! third_body check
+      ratec = pfunc2(nlayer,t,dens,pfunc2_param(nb_pfunc2))
+    else
+      ratec = pfunc2(nlayer,t,c(:,third_body(ireact)),pfunc2_param(nb_pfunc2))
+    end if
+  else if (rtype(ireact)==3) then
+    nb_pfunc3 = nb_pfunc3 + 1
+    if (third_body(ireact)==0) then            !! third_body check
+      ratec = pfunc3(nlayer,t,dens,pfunc3_param(nb_pfunc3))
+    else
+      ratec = pfunc3(nlayer,t,c(:,third_body(ireact)),pfunc3_param(nb_pfunc3))
+    end if
+  else
+    print*, 'Error in reactionrates: wrong index coefficient rate v3'
+    print*, 'rtype(ireact) = ',rtype(ireact)
+    print*, 'It should be 1 or 2 '
+    print*, 'Please check the reaction ',ireact-nb_phot_max-nb_reaction_4_max
+    if (ireact>nreact) print*, 'Please check the reaction count into the code beacause ireact > nreact is not possible'
+    stop
+  end if
+
+  ! fill the right type index
+  nb_reaction_3 = nb_reaction_3 + 1
+  v_3(:,nb_reaction_3) = ratec(:)
+  if (three_prod(ireact)) then             ! three_prod check
+    nb_reaction_3 = nb_reaction_3 + 1
+    v_3(:,nb_reaction_3) = ratec(:)
+  end if
+
+  ireact = ireact + 1
+
+end do
+
+call reactionrates_HC(nlayer,nesp,dens,t,press,zmmean,sza,c,v_phot,v_4,v_3,nb_phot,nb_reaction_4,nb_reaction_3)
+
+!===========================================================
+!  check dimensions 
+!===========================================================
+
+if (firstcall) then
+ ireact = ireact-1
+ if (ireact /= nreact) print*, 'wrong dimensions in reactionrates : number of reaction should be ', nreact,' and not ', ireact
+ if ((nb_phot /= nb_phot_max)             .or.  &
+     (nb_reaction_3 /= nb_reaction_3_max) .or.  &
+     (nb_reaction_4 /= nb_reaction_4_max)) then
+    print*, 'nb_phot       = ', nb_phot
+    print*, 'nb_reaction_4 = ', nb_reaction_4
+    print*, 'nb_reaction_3 = ', nb_reaction_3
+    print*, 'wrong dimensions in reactionrates' 
+    print*, 'nb_phot_max       = ', nb_phot_max
+    print*, 'nb_reaction_4_max = ', nb_reaction_4_max
+    print*, 'nb_reaction_3_max = ', nb_reaction_3_max
+    print*, '------ hard coding reaction ------'
+    print*, 'nphot_hard_coding = ', nphot_hard_coding
+    print*, 'n4_hard_coding    = ', n4_hard_coding
+    print*, 'n3_hard_coding    = ', n3_hard_coding
+    stop
+ end if
+ if ((nb_hv /= nb_hv_max)         .or.  &
+     (nb_pfunc1 /= nb_pfunc1_max) .or.  &
+     (nb_pfunc2 /= nb_pfunc2_max) .or.  &
+     (nb_pfunc3 /= nb_pfunc3_max)) then
+    print*, 'nb_hv         = ', nb_hv
+    print*, 'nb_pfunc1     = ', nb_pfunc1
+    print*, 'nb_pfunc2     = ', nb_pfunc2
+    print*, 'nb_pfunc3     = ', nb_pfunc3
+    print*, 'wrong dimensions in reactionrates'
+    print*, 'nb_hv_max     = ', nb_hv_max
+    print*, 'nb_pfunc1_max = ', nb_pfunc1_max
+    print*, 'nb_pfunc2_max = ', nb_pfunc2_max
+    print*, 'nb_pfunc3_max = ', nb_pfunc3_max
+    stop
+ end if
+firstcall = .false.
+end if 
 
 end subroutine reactionrates
@@ -1013 +735 @@
 
 use types_asis
+use chimiedata_h
 
 implicit none
-
-#include "chimiedata.h"
 
 ! input
@@ -1122 +843 @@
 !================================================================
 
- subroutine indice(i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
-                   i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
-                   i_n, i_n2d, i_no, i_no2, i_n2)
+ subroutine indice(nreact,rtype,third_body,three_prod)
 
 !================================================================
@@ -1139 +858 @@
 
 use types_asis
+use datafile_mod
+use ioipsl_getin_p_mod, only: getin_p
+use tracer_h, only: nesp
+use chimiedata_h
+use callkeys_mod, only: jonline
 
 implicit none
 
-#include "chimiedata.h"
-
-! input
-
-integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
-           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
-           i_n, i_n2d, i_no, i_no2, i_n2
+! output
+
+integer, intent(out)              :: nreact          ! number of reactions in reactions files
+integer, allocatable, intent(out) :: rtype(:)        ! reaction rate type
+integer, allocatable, intent(out) :: third_body(:)   ! if the reaction have a third body: index of the third body, else zero
+logical, allocatable, intent(out) :: three_prod(:)   ! if the reaction have a three defferent proucts egal .true.
 
 ! local
 
 integer :: nb_phot, nb_reaction_3, nb_reaction_4
-integer :: i_dummy
-
-allocate (indice_phot(nb_phot_max))
-allocate (indice_3(nb_reaction_3_max))
-allocate (indice_4(nb_reaction_4_max))
-
-i_dummy = 1
-
-nb_phot       = 0
-nb_reaction_3 = 0
-nb_reaction_4 = 0
+integer :: nb_hv, nb_pfunc1, nb_pfunc2, nb_pfunc3
+integer :: iq, ireact
+
+character(len = 128)              :: reactfile       ! reactions table file name
+character(len = 128)              :: monoreact       ! photochemical reactions table file name
+character(len = 128)              :: bimolreact      ! bimolecular reactions table file name
+character(len = 128)              :: quadrareact     ! quadratic reactions table file name
+integer                           :: nbq             ! number of species in reactions file
+integer                           :: pnlines         ! number of lines in photochemical reactions file
+integer                           :: bnlines         ! number of lines in bimolecular reactions file
+integer                           :: qnlines         ! number of lines in quadratic reactions file
+integer                           :: pnreact         ! number of reaction in photochemical reactions file
+integer                           :: bnreact         ! number of reaction in bimolecular reactions file
+integer                           :: qnreact         ! number of reaction in quadratic reactions file
+integer                           :: specheck(nesp)  ! to count the species of traceur.def in reactions file
+integer                           :: specheckr(nesp) ! to count the reactant species of traceur.def in reactions file
+integer                           :: specheckp(nesp) ! to count the product species of traceur.def in reactions file
+
+nb_phot        = 0
+nb_reaction_3  = 0
+nb_reaction_4  = 0
+nb_phot_hv_max = 0
+
+nb_hv          = 0
+nb_pfunc1      = 0
+nb_pfunc2      = 0
+nb_pfunc3      = 0
 
 !===========================================================
-!      O2 + hv -> O + O
+!  Read chemical reactions
 !===========================================================
 
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_o2, 2.0, i_o, 0.0, i_dummy) 
+! Get number of reactions
+pnlines        = 0
+bnlines        = 0
+qnlines        = 0
+nreact         = 0
+pnreact        = 0
+bnreact        = 0
+qnreact        = 0
+
+write(*,*) "Read photochemical reaction"
+ monoreact = "monoreact" ! default
+call getin_p("monoreact",monoreact) ! default path
+write(*,*) " monoreact = ",trim(monoreact)
+
+write(*,*) "Read bimolecular reaction"
+ bimolreact = "bimolreact" ! default
+call getin_p("bimolreact",bimolreact) ! default path
+write(*,*) " bimolreact = ",trim(bimolreact)
+
+write(*,*) "Read quadratic reaction"
+ quadrareact = "quadrareact" ! default
+call getin_p("quadrareact",quadrareact) ! default path
+write(*,*) " quadrareact = ",trim(quadrareact)
+
+call count_react(monoreact,pnlines,pnreact,nb_phot,nb_hv,nb_pfunc1,nb_pfunc2,nb_pfunc3)
+call count_react(bimolreact,bnlines,bnreact,nb_reaction_4,nb_hv,nb_pfunc1,nb_pfunc2,nb_pfunc3)
+call count_react(quadrareact,qnlines,qnreact,nb_reaction_3,nb_hv,nb_pfunc1,nb_pfunc2,nb_pfunc3)
+nreact = pnreact + bnreact + qnreact
+
+!!!!!!!!!!! Hard coding reaction !!!!!!!!!!!!!!!!!!!!!!!!!!!
+nb_phot       = nb_phot       + nphot_hard_coding
+nb_reaction_4 = nb_reaction_4 + n4_hard_coding
+nb_reaction_3 = nb_reaction_3 + n3_hard_coding
+!!!!!!!!!!! END Hard coding reaction !!!!!!!!!!!!!!!!!!!!!!!
+
+write(*,*)'number of reaction in reaction files is = ',nreact
+print*, 'nb_phot       = ', nb_phot
+print*, 'nb_reaction_4 = ', nb_reaction_4
+print*, 'nb_reaction_3 = ', nb_reaction_3
+print*, '****************'
+print*, 'nb_hv         = ', nb_hv
+print*, 'nb_pfunc1     = ', nb_pfunc1
+print*, 'nb_pfunc2     = ', nb_pfunc2
+print*, 'nb_pfunc3     = ', nb_pfunc3
+nb_phot_max       = nb_phot
+nb_reaction_4_max = nb_reaction_4
+nb_reaction_3_max = nb_reaction_3
+nb_hv_max         = nb_hv
+nb_pfunc1_max     = nb_pfunc1
+nb_pfunc2_max     = nb_pfunc2
+nb_pfunc3_max     = nb_pfunc3
+
+! Allocate
+allocate(indice_phot(nb_phot_max))
+allocate(indice_4(nb_reaction_4_max))
+allocate(indice_3(nb_reaction_3_max))
+allocate(rtype(nreact))
+allocate(third_body(nreact))
+allocate(three_prod(nreact))
+allocate(jlabel(nb_phot_hv_max,2))
+allocate(jparamline(nb_phot_hv_max))
+allocate(pfunc1_param(nb_pfunc1_max))
+allocate(pfunc2_param(nb_pfunc2_max))
+allocate(pfunc3_param(nb_pfunc3_max))
+
+! Get reactants, products and number of species in reactfile
+! inititialize variables
+nbq             = 0
+specheck(:)     = 0
+specheckr(:)    = 0
+specheckp(:)    = 0
+rtype(:)        = 0
+third_body(:)   = 0
+pfunc1_param(:) = rtype1(0.,0.,0.,0.,0.)
+pfunc2_param(:) = rtype2(0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.)
+pfunc3_param(:) = rtype3(0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.)
+nb_pfunc1       = 0
+nb_pfunc2       = 0
+nb_pfunc3       = 0
+three_prod(:)   = .false.
+jlabel(:,:)     = ''
+jparamline(:)   = ''
+
+call get_react(monoreact,pnlines,pnreact,rtype(1:pnreact),third_body(1:pnreact),three_prod(1:pnreact),    &
+               nb_phot,specheck,specheckr,specheckp,'vphot',nbq,nb_pfunc1,nb_pfunc2,nb_pfunc3)
+call get_react(bimolreact,bnlines,bnreact,rtype(pnreact+1:pnreact+bnreact),third_body(pnreact+1:pnreact+bnreact), &
+               three_prod(pnreact+1:pnreact+bnreact),nb_reaction_4,specheck,specheckr,specheckp,'v4',nbq,nb_pfunc1,nb_pfunc2,nb_pfunc3)
+call get_react(quadrareact,qnlines,qnreact,rtype(pnreact+bnreact+1:nreact),third_body(pnreact+bnreact+1:nreact),  &
+               three_prod(pnreact+bnreact+1:nreact),nb_reaction_3,specheck,specheckr,specheckp,'v3',nbq,nb_pfunc1,nb_pfunc2,nb_pfunc3)
+
+! rewrite jlabel corretly for output
+do ireact=1,nb_phot_hv_max
+  if (three_prod(ireact)) then
+    jlabel(ireact+1:nb_phot_hv_max,2) = jlabel(ireact:nb_phot_hv_max-1,2)
+    jlabel(ireact+1:nb_phot_hv_max,1) = jlabel(ireact:nb_phot_hv_max-1,1)
+    jlabel(ireact,1) = trim(jlabel(ireact,1))//'_a'
+    jlabel(ireact+1,1) = trim(jlabel(ireact+1,1))//'_b'
+  end if
+end do
+
+!!!!!!!!!!! Hard coding reaction !!!!!!!!!!!!!!!!!!!!!!!!!!!
+call indice_HC(nb_phot,nb_reaction_4,nb_reaction_3)
+!!!!!!!!!!! END Hard coding reaction !!!!!!!!!!!!!!!!!!!!!!!
+
+write(*,*)'number of species in reaction files is = ',nbq
+write(*,*)'species in reaction files are:'
+
+do iq=1,nesp
+  if (specheck(iq)==1) print*, chemnoms(iq)
+end do
 
 !===========================================================
-!      O2 + hv -> O + O(1D)
+!  check species only destroy or product
 !===========================================================
 
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_o2, 1.0, i_o, 1.0, i_o1d) 
+do iq=1,nesp
+  if (specheckr(iq)/=specheckp(iq)) then
+    if (specheckr(iq)==0 .and. specheckp(iq)==1) then
+      print*, 'WARNING: ', chemnoms(iq),' is only product'
+    else if (specheckr(iq)==1 .and. specheckp(iq)==0) then
+      print*, 'WARNING: ', chemnoms(iq),' is only destroy'
+    else
+      print*, 'Error in indice: bad value in specheckr or specheckp'
+      stop
+    end if
+  end if
+end do 
 
 !===========================================================
-!      CO2 + hv -> CO + O
+!  check stochiometry
 !===========================================================
 
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co, 1.0, i_o) 
-
-!===========================================================
-!      CO2 + hv -> CO + O(1D)
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co, 1.0, i_o1d) 
-
-!===========================================================
-!      O3 + hv -> O2 + O(1D)
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_o3, 1.0, i_o2, 1.0, i_o1d) 
-
-!===========================================================
-!      O3 + hv -> O2 + O
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_o3, 1.0, i_o2, 1.0, i_o) 
-
-!===========================================================
-!      H2O + hv -> H + OH
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_h2o, 1.0, i_h, 1.0, i_oh) 
-
-!===========================================================
-!      H2O2 + hv -> OH + OH
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 2.0, i_oh, 0.0, i_dummy) 
-
-!===========================================================
-!      HO2 + hv -> OH + O
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_ho2, 1.0, i_oh, 1.0, i_o) 
-
-!===========================================================
-!      NO + hv -> N + O
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_n, 1.0, i_o) 
-
-!===========================================================
-!      NO2 + hv -> NO + O
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_no2, 1.0, i_no, 1.0, i_o) 
-
-!===========================================================
-!      a001 : O + O2 + CO2 -> O3 + CO2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o2, 1.0, i_o3, 0.0, i_dummy) 
-
-!===========================================================
-!      a002 : O + O + CO2 -> O2 + CO2
-!===========================================================
-
-nb_reaction_3 = nb_reaction_3 + 1
-
-indice_3(nb_reaction_3) = z3spec(2.0, i_o, 1.0, i_o2, 0.0, i_dummy) 
-
-!===========================================================
-!      a003 : O + O3 -> O2 + O2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy) 
-
-!===========================================================
-!      b001 : O(1D) + CO2 -> O + CO2
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy) 
-
-!===========================================================
-!      b002 : O(1D) + H2O -> OH + OH
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2o, 2.0, i_oh, 0.0, i_dummy) 
-
-!===========================================================
-!      b003 : O(1D) + H2 -> OH + H
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_h2, 1.0, i_oh, 1.0, i_h) 
-
-!===========================================================
-!      b004 : O(1D) + O2 -> O + O2
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_o1d, 1.0, i_o, 0.0, i_dummy) 
-
-!===========================================================
-!      b005 : O(1D) + O3 -> O2 + O2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 2.0, i_o2, 0.0, i_dummy) 
-
-!===========================================================
-!      b006 : O(1D) + O3 -> O2 + O + O
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o1d, 1.0, i_o3, 1.0, i_o2, 2.0, i_o) 
-
-!===========================================================
-!      c001 : O + HO2 -> OH + O2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_ho2, 1.0, i_oh, 1.0, i_o2) 
-
-!===========================================================
-!      c002 : O + OH -> O2 + H
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_oh, 1.0, i_o2, 1.0, i_h) 
-
-!===========================================================
-!      c003 : H + O3 -> OH + O2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o3, 1.0, i_oh, 1.0, i_o2) 
-
-!===========================================================
-!      c004 : H + HO2 -> OH + OH
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 2.0, i_oh, 0.0, i_dummy) 
-
-!===========================================================
-!      c005 : H + HO2 -> H2 + O2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2, 1.0, i_o2) 
-
-!===========================================================
-!      c006 : H + HO2 -> H2O + O
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o) 
-
-!===========================================================
-!      c007 : OH + HO2 -> H2O + O2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_ho2, 1.0, i_h2o, 1.0, i_o2) 
-
-!===========================================================
-!      c008 : HO2 + HO2 -> H2O2 + O2
-!===========================================================
-
-nb_reaction_3 = nb_reaction_3 + 1
-
-indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2) 
-
-!===========================================================
-!      c009 : OH + H2O2 -> H2O + HO2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2o2, 1.0, i_h2o, 1.0, i_ho2) 
-
-!===========================================================
-!      c010 : OH + H2 -> H2O + H
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_h2, 1.0, i_h2o, 1.0, i_h) 
-
-!===========================================================
-!      c011 : H + O2 + CO2 -> HO2 + CO2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_h, 1.0, i_o2, 1.0, i_ho2, 0.0, i_dummy) 
-
-!===========================================================
-!      c012 : O + H2O2 -> OH + HO2
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_o, 1.0, i_h2o2, 1.0, i_oh, 1.0, i_ho2) 
-
-!===========================================================
-!      c013 : OH + OH -> H2O + O
-!===========================================================
-
-nb_reaction_3 = nb_reaction_3 + 1
-
-indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o, 1.0, i_o) 
-
-!===========================================================
-!      c014 : OH + O3 -> HO2 + O2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_oh, 1.0, i_o3, 1.0, i_ho2, 1.0, i_o2) 
-
-!===========================================================
-!      c015 : HO2 + O3 -> OH + O2 + O2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_ho2, 1.0, i_o3, 1.0, i_oh, 2.0, i_o2) 
-
-!===========================================================
-!      c016 : HO2 + HO2 + CO2 -> H2O2 + O2 + CO2 
-!===========================================================
-
-nb_reaction_3 = nb_reaction_3 + 1
-
-indice_3(nb_reaction_3) = z3spec(2.0, i_ho2, 1.0, i_h2o2, 1.0, i_o2) 
-
-!===========================================================
-!      c017 : OH + OH + CO2 -> H2O2 + CO2 
-!===========================================================
-
-nb_reaction_3 = nb_reaction_3 + 1
-
-indice_3(nb_reaction_3) = z3spec(2.0, i_oh, 1.0, i_h2o2, 0.0, i_dummy) 
-
-!===========================================================
-!      c018 : H + H + CO2 -> H2 + CO2 
-!===========================================================
-
-nb_reaction_3 = nb_reaction_3 + 1
-
-indice_3(nb_reaction_3) = z3spec(2.0, i_h, 1.0, i_h2, 0.0, i_dummy) 
-
-!===========================================================
-!      d001 : NO2 + O -> NO + O2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_o, 1.0, i_no, 1.0, i_o2) 
-
-!===========================================================
-!      d002 : NO + O3 -> NO2 + O2 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_o3, 1.0, i_no2, 1.0, i_o2) 
-
-!===========================================================
-!      d003 : NO + HO2 -> NO2 + OH 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_ho2, 1.0, i_no2, 1.0, i_oh) 
-
-!===========================================================
-!      d004 : N + NO -> N2 + O 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_no, 1.0, i_n2, 1.0, i_o) 
-
-!===========================================================
-!      d005 : N + O2 -> NO + O 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o2, 1.0, i_no, 1.0, i_o) 
-
-!===========================================================
-!      d006 : NO2 + H -> NO + OH 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_h, 1.0, i_no, 1.0, i_oh) 
-
-!===========================================================
-!      d007 : N + O -> NO
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o, 1.0, i_no, 0.0, i_dummy) 
-
-!===========================================================
-!      d008 : N + HO2 -> NO + OH
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_ho2, 1.0, i_no, 1.0, i_oh) 
-
-!===========================================================
-!      d009 : N + OH -> NO + H
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_oh, 1.0, i_no, 1.0, i_h) 
-
-!===========================================================
-!      e001 : CO + OH -> CO2 + H 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_oh, 1.0, i_co2, 1.0, i_h) 
-
-!===========================================================
-!      e002 : CO + O + M -> CO2 + M 
-!===========================================================
-
-nb_reaction_4 = nb_reaction_4 + 1
-
-indice_4(nb_reaction_4) = z4spec(1.0, i_co, 1.0, i_o, 1.0, i_co2, 0.0, i_dummy) 
-
-!===========================================================
-!      h001: HO2 + ice -> products
-!            treated as
-!            HO2 -> 0.5 H2O + 0.75 O2
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_ho2, 0.5, i_h2o, 0.75, i_o2) 
-
-!===========================================================
-!      h002: OH + ice -> products
-!            treated as
-!            OH -> 0.5 H2O + 0.25 O2
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_oh, 0.5, i_h2o, 0.25, i_o2) 
-
-!===========================================================
-!      h003: H2O2 + ice -> products
-!            treated as
-!            H2O2 -> H2O + 0.5 O2
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 1.0, i_h2o, 0.5, i_o2) 
-
-!===========================================================
-!      h004: HO2 + dust -> products
-!            treated as
-!            HO2 -> 0.5 H2O + 0.75 O2
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_ho2, 0.5, i_h2o, 0.75, i_o2) 
-
-!===========================================================
-!      h005: H2O2 + dust -> products
-!            treated as
-!            H2O2 -> H2O + 0.5 O2
-!===========================================================
-
-nb_phot = nb_phot + 1
-
-indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 1.0, i_h2o, 0.5, i_o2) 
+! If you found a way
 
 !===========================================================
@@ -1611 +1046 @@
 !===========================================================
 
-print*, 'nb_phot       = ', nb_phot
-print*, 'nb_reaction_4 = ', nb_reaction_4
-print*, 'nb_reaction_3 = ', nb_reaction_3
+if (jonline) then
+  nd = nb_hv_max
+else if (nb_phot_hv_max /= 0) then
+  print*,'calchim: Read photolysis lookup table'
+  call read_phototable
+end if
 
 if ((nb_phot /= nb_phot_max)             .or.  &
     (nb_reaction_3 /= nb_reaction_3_max) .or.  &
-    (nb_reaction_4 /= nb_reaction_4_max)) then
+    (nb_reaction_4 /= nb_reaction_4_max) .or.  &
+    (nd /= nb_hv_max)) then
+   print*, 'nb_phot       = ', nb_phot
+   print*, 'nb_reaction_4 = ', nb_reaction_4
+   print*, 'nb_reaction_3 = ', nb_reaction_3
+   print*, 'nd            = ', nd
    print*, 'wrong dimensions in indice' 
+   print*, 'nb_phot_max       = ', nb_phot_max
+   print*, 'nb_reaction_4_max = ', nb_reaction_4_max
+   print*, 'nb_reaction_3_max = ', nb_reaction_3_max
+   print*, 'nb_phot_hv_max    = ', nb_phot_hv_max
+   print*, 'nb_hv_max         = ', nb_hv_max
    stop
 end if  
@@ -1626 +1074 @@
 !*****************************************************************
 
-      subroutine gcmtochim(nlayer, nq, zycol, lswitch, nesp,         &
-                           i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
-                           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,         &
-                           i_n, i_n2d, i_no, i_no2, i_n2,            &
-                           dens, rm, c) 
+      subroutine gcmtochim(nlayer, zycol, lswitch, nesp, dens, rm, c) 
 
 !*****************************************************************
-
-      use tracer_h, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,         &
-     &                      igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,  &
-     &                      igcm_ho2, igcm_h2o2, igcm_h2o_vap,           &
-     &                      igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2
 
       use callkeys_mod
@@ -1649 +1088 @@
       
       integer, intent(in) :: nlayer ! number of atmospheric layers
-      integer, intent(in) :: nq     ! number of tracers in the gcm
-      integer :: nesp               ! number of species in the chemistry
+      integer, intent(in) :: nesp   ! number of species in the chemistry
       integer :: lswitch            ! interface level between chemistries
 
-      integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,           &
-                 i_h2, i_oh, i_ho2, i_h2o2, i_h2o,                   &
-                 i_n, i_n2d, i_no, i_no2, i_n2
-
-      real :: zycol(nlayer,nq)      ! volume mixing ratios in the gcm
+      real :: zycol(nlayer,nesp)    ! volume mixing ratios in the gcm
       real :: dens(nlayer)          ! total number density (molecule.cm-3) 
 
@@ -1672 +1106 @@
 
       integer      :: l, iesp
-      logical,save :: firstcall = .true.
       
+      rm(:,:) = 0.
       
-!     first call initializations
-
-      if (firstcall) then
-
-!       identify the indexes of the tracers we need
-
-         if (igcm_co2 == 0) then
-            write(*,*) "gcmtochim: Error; no CO2 tracer !!!"
-            stop
-         endif
-         if (igcm_co == 0) then
-            write(*,*) "gcmtochim: Error; no CO tracer !!!"
-            stop
-         end if
-         if (igcm_o == 0) then
-            write(*,*) "gcmtochim: Error; no O tracer !!!"
-            stop
-         end if
-         if (igcm_o1d == 0) then
-            write(*,*) "gcmtochim: Error; no O1D tracer !!!"
-            stop
-         end if
-         if (igcm_o2 == 0) then
-            write(*,*) "gcmtochim: Error; no O2 tracer !!!"
-            stop
-         end if
-         if (igcm_o3 == 0) then
-            write(*,*) "gcmtochim: Error; no O3 tracer !!!"
-            stop
-         end if
-         if (igcm_h == 0) then
-            write(*,*) "gcmtochim: Error; no H tracer !!!"
-            stop
-         end if
-         if (igcm_h2 == 0) then
-            write(*,*) "gcmtochim: Error; no H2 tracer !!!"
-            stop
-         end if
-         if (igcm_oh == 0) then
-            write(*,*) "gcmtochim: Error; no OH tracer !!!"
-            stop
-         end if
-         if (igcm_ho2 == 0) then
-            write(*,*) "gcmtochim: Error; no HO2 tracer !!!"
-            stop
-         end if
-         if (igcm_h2o2 == 0) then
-            write(*,*) "gcmtochim: Error; no H2O2 tracer !!!"
-            stop
-         end if
-         if (igcm_n == 0) then
-            write(*,*) "gcmtochim: Error; no N tracer !!!"
-!            stop
-         end if
-         if (igcm_n2d == 0) then
-            write(*,*) "gcmtochim: Error; no N2D tracer !!!"
-!            stop
-         end if
-         if (igcm_no == 0) then
-            write(*,*) "gcmtochim: Error; no NO tracer !!!"
-!            stop
-         end if
-         if (igcm_no2 == 0) then
-            write(*,*) "gcmtochim: Error; no NO2 tracer !!!"
-!            stop
-         end if
-         if (igcm_n2 == 0) then
-            write(*,*) "gcmtochim: Error; no N2 tracer !!!"
-            stop
-         end if
-         if (igcm_h2o_vap == 0) then
-            write(*,*) "gcmtochim: Error; no water vapor tracer !!!"
-            stop
-         end if
-         firstcall = .false.
-      end if ! of if (firstcall)
-
 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 !     initialise mixing ratios
@@ -1757 +1114 @@
 
       do l = 1,lswitch-1
-         rm(l,i_co2)  = zycol(l, igcm_co2)
-         rm(l,i_co)   = zycol(l, igcm_co)
-         rm(l,i_o)    = zycol(l, igcm_o)
-         rm(l,i_o1d)  = zycol(l, igcm_o1d)
-         rm(l,i_o2)   = zycol(l, igcm_o2)
-         rm(l,i_o3)   = zycol(l, igcm_o3)
-         rm(l,i_h)    = zycol(l, igcm_h)
-         rm(l,i_h2)   = zycol(l, igcm_h2)
-         rm(l,i_oh)   = zycol(l, igcm_oh)
-         rm(l,i_ho2)  = zycol(l, igcm_ho2)
-         rm(l,i_h2o2) = zycol(l, igcm_h2o2)
-         rm(l,i_h2o)  = zycol(l, igcm_h2o_vap)
-         rm(l,i_n)    = zycol(l, igcm_n)
-         rm(l,i_n2d)  = zycol(l, igcm_n2d)
-         rm(l,i_no)   = zycol(l, igcm_no)
-         rm(l,i_no2)  = zycol(l, igcm_no2)
-         rm(l,i_n2)   = zycol(l, igcm_n2)
+         rm(l,:) = zycol(l,:)
       end do 
 
@@ -1794 +1135 @@
 !*****************************************************************
  
-      subroutine chimtogcm(nlayer, nq, zycol, lswitch, nesp,         &
-                           i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &
-                           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,         &
-                           i_n, i_n2d, i_no, i_no2, i_n2, dens, c) 
+      subroutine chimtogcm(nlayer, zycol, lswitch, nesp, dens, c) 
+
  
 !*****************************************************************
- 
-      use tracer_h, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,          &
-                            igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,   &
-                            igcm_ho2, igcm_h2o2, igcm_h2o_vap,            &
-                            igcm_n, igcm_n2d, igcm_no, igcm_no2, igcm_n2
 
       use callkeys_mod
@@ -1816 +1150 @@
  
       integer, intent(in) :: nlayer  ! number of atmospheric layers
-      integer, intent(in) :: nq      ! number of tracers in the gcm
-      integer :: nesp                ! number of species in the chemistry
+      integer, intent(in) :: nesp    ! number of species in the chemistry
       integer :: lswitch             ! interface level between chemistries
-      integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,       &
-                 i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               &
-                 i_n, i_n2d, i_no, i_no2, i_n2
 
       real :: dens(nlayer)     ! total number density (molecule.cm-3) 
@@ -1830 +1160 @@
 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        
-      real zycol(nlayer,nq)  ! volume mixing ratios in the gcm
+      real zycol(nlayer,nesp)  ! volume mixing ratios in the gcm
 
 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
@@ -1843 +1173 @@
 
       do l = 1,lswitch-1
-         zycol(l, igcm_co2)     = c(l,i_co2)/dens(l) 
-         zycol(l, igcm_co)      = c(l,i_co)/dens(l) 
-         zycol(l, igcm_o)       = c(l,i_o)/dens(l) 
-         zycol(l, igcm_o1d)     = c(l,i_o1d)/dens(l)
-         zycol(l, igcm_o2)      = c(l,i_o2)/dens(l) 
-         zycol(l, igcm_o3)      = c(l,i_o3)/dens(l) 
-         zycol(l, igcm_h)       = c(l,i_h)/dens(l)  
-         zycol(l, igcm_h2)      = c(l,i_h2)/dens(l) 
-         zycol(l, igcm_oh)      = c(l,i_oh)/dens(l) 
-         zycol(l, igcm_ho2)     = c(l,i_ho2)/dens(l) 
-         zycol(l, igcm_h2o2)    = c(l,i_h2o2)/dens(l)
-         zycol(l, igcm_h2o_vap) = c(l,i_h2o)/dens(l)
-         zycol(l, igcm_n)       = c(l,i_n)/dens(l)
-         zycol(l, igcm_n2d)     = c(l,i_n2d)/dens(l)
-         zycol(l, igcm_no)      = c(l,i_no)/dens(l)
-         zycol(l, igcm_no2)     = c(l,i_no2)/dens(l)
-         zycol(l, igcm_n2)      = c(l,i_n2)/dens(l)
+         zycol(l,:)  = c(l,:)/dens(l) 
       end do 
 
       end subroutine chimtogcm
 
+!*****************************************************************
+
+    subroutine split_str(line,words,n,nmax)
+
+!*****************************************************************
+
+        implicit none
+        character(*), intent(in)  :: line
+        integer,      intent(in)  :: nmax
+        character(*), intent(out) :: words(nmax)
+        integer,      intent(out) :: n          ! number of words at the end
+        integer :: ios
+        character(100) :: buf(100)  ! large buffer!
+    
+        n = 0
+        do
+            n = n + 1
+            read(line,*,iostat=ios) buf(1:n)  ! use list-directed input
+            if (ios==0) then
+                words(1:n) = buf(1:n)   ! if success, copy to the original array
+            else
+                n = n-1
+                exit       ! if all the words are obtained, finish
+            endif
+            if (n>nmax) then
+              print*,'Error in split_str: to much words'
+              print*,'limit = ',nmax
+              print*,'change it, if you want, with increasing nmax'
+              stop
+            end if
+        enddo
+    end subroutine split_str
+
+!*****************************************************************
+
+    subroutine count_react(reactfile,nlines,nreact,nrtype,nb_hv,nb_pfunc1,nb_pfunc2,nb_pfunc3)
+
+!*****************************************************************
+
+        use types_asis, only : nb_phot_hv_max
+
+        implicit none
+        character(*), intent(in)     :: reactfile ! name of the file to read 
+        integer,      intent(inout)  :: nlines    ! number of lines in the file
+        integer,      intent(out)    :: nreact    ! real number of reaction
+        integer,      intent(inout)  :: nrtype    ! number of reaction for calculation
+        integer,      intent(inout)  :: nb_hv, nb_pfunc1, nb_pfunc2, nb_pfunc3
+
+        ! local
+        character(len = 50) :: reactline          ! all reactants of one reaction
+        character(len = 50) :: prodline           ! all produts of one reaction
+        integer             :: typerate           ! get the type of the rate reaction coefficient (pfunc_i)
+        integer             :: nwp                ! number of products for a reaction
+        integer,parameter   :: nmax = 5           ! number max of reactants or products
+        character(len = 24) :: words(nmax)        ! to get words in reactants and products line
+        integer             :: ierr
+
+        nreact     = 0
+
+        open(402, form = 'formatted', status = 'old',                &
+              file ='chemnetwork/'//trim(reactfile),iostat=ierr)
+
+        if (ierr /= 0) THEN
+           write(*,*)'Error : cannot open chemical reaction file : chemnetwork/'//trim(reactfile)
+           write(*,*)'It should be in your simulation directory in chemnetwork directory'
+           write(*,*)'   You can change the input chemical reactions file name in'
+           write(*,*)'   callphys.def with:'
+           write(*,*)'   monoreact=filename or bimolreact=filename or quadrareact=filename'
+           write(*,*)'   depending on what chemical reaction type it is'
+           stop
+        end if
+    
+        do
+          read(402,'(A,A,I2)',iostat=ierr) reactline, prodline, typerate
+          if (ierr<0) exit
+          nlines = nlines + 1
+          if (reactline(1:1)/='!' .and. reactline(1:1)/='') then
+            ! count the number of reaction
+            nreact = nreact + 1
+            call split_str(prodline,words,nwp,nmax)
+            nrtype = nrtype + 1
+            ! check three product reaction
+            if (nwp>2 .and. trim(words(1))/=trim(words(2)) &
+                      .and. trim(words(1))/=trim(words(3)) &
+                      .and. trim(words(2))/=trim(words(3))) nrtype = nrtype + 1
+
+            ! count the number of each rate reaction coefficient type
+            if (typerate==0) then
+              nb_hv = nb_hv + 1
+              nb_phot_hv_max = nb_phot_hv_max + 1
+              if (nwp>2 .and. trim(words(1))/=trim(words(2)) &
+                        .and. trim(words(1))/=trim(words(3)) &
+                        .and. trim(words(2))/=trim(words(3))) nb_phot_hv_max = nb_phot_hv_max + 1
+            else if (typerate==1) then
+              nb_pfunc1 = nb_pfunc1 + 1
+            else if (typerate==2) then
+              nb_pfunc2 = nb_pfunc2 + 1
+            else if (typerate==3) then
+              nb_pfunc3 = nb_pfunc3 + 1
+            else
+              print*, 'Error in indice: wrong index coefficient rate line ',nlines
+              print*, 'in file : chemnetwork/'//trim(reactfile)
+              print*, 'It should be 0 for photolysis reations and 1 or 2 for the others'
+              print*, 'And not : ', typerate
+              stop
+            end if
+
+          end if
+
+        end do
+
+        close(402)
+
+    end subroutine count_react
+
+!*****************************************************************
+
+    subroutine get_react(reactfile,nlines,nreact,rtype,third_body,three_prod, &
+                         nrtype,specheck,specheckr,specheckp,typeindice,nbq,  &
+                         init_nb_pfunc1,init_nb_pfunc2,init_nb_pfunc3)
+
+!*****************************************************************
+
+        use types_asis
+        use tracer_h
+        use chimiedata_h, only: indexchim
+        use callkeys_mod, only: jonline
+
+        implicit none
+        character(*), intent(in)     :: reactfile          ! name of the file to read 
+        integer,      intent(in)     :: nlines             ! number of lines in the file
+        integer,      intent(in)     :: nreact             ! real number of reaction in the file
+        integer,      intent(inout)  :: rtype(nreact)      ! reaction rate type
+        integer,      intent(inout)  :: third_body(nreact) ! if the reaction have a third body: index of the third body, else zero
+        logical,      intent(inout)  :: three_prod(nreact) ! if the reaction have a three defferent proucts egal .true.
+        integer,      intent(out)    :: nrtype             ! number of calculation reactions
+        integer,      intent(inout)  :: specheck(nesp)     ! to count the species of traceur.def in reactions file
+        integer,      intent(inout)  :: specheckr(nesp)    ! to count the reactant species of traceur.def in reactions file
+        integer,      intent(inout)  :: specheckp(nesp)    ! to count the product species of traceur.def in reactions file
+        character(*), intent(in)     :: typeindice         ! reaction type (v3, v4 or vphot)
+        integer,      intent(inout)  :: nbq                ! number of species in reactions file
+        integer,      intent(inout)  :: init_nb_pfunc1     ! in : initial value of nb_pfunc1 - out : final value of nb_pfunc1
+        integer,      intent(inout)  :: init_nb_pfunc2     ! in : initial value of nb_pfunc2 - out : final value of nb_pfunc2
+        integer,      intent(inout)  :: init_nb_pfunc3     ! in : initial value of nb_pfunc3 - out : final value of nb_pfunc3
+
+        ! local
+        character(len = 50)  :: reactline                  ! all reactants of one reaction
+        character(len = 50)  :: prodline                   ! all produts of one reaction
+        integer              :: nwr                        ! number of reactants for a reaction
+        integer              :: nwp                        ! number of products for a reaction
+        integer,parameter    :: nmax = 5                   ! number max of reactants or products
+        character(len = 24)  :: words(nmax)                ! to get words in reactants and products line
+        real                 :: nindice(2*nmax)            ! stoichiometry of species (for indice variables)
+        integer              :: iindice(2*nmax)            ! indice of species (for indice variables)
+        character(len = 500) :: paramline                  ! line of reactions parameters
+        character(len = 50)  :: reactants(nreact,nmax)     ! reactions reactants
+        character(len = 50)  :: products(nreact,nmax)      ! reactions products
+        logical              :: spedouble                  ! check if a specie appears twice in reactants or products
+        integer              :: ierr, ilines, ireact, i_dummy, iw, iwhere, i
+        integer              :: nb_hv, nb_pfunc1, nb_pfunc2, nb_pfunc3
+
+        i_dummy   = 1
+        nrtype    = 0
+        ireact    = 0
+        nb_hv     = 0
+        nb_pfunc1 = init_nb_pfunc1
+        nb_pfunc2 = init_nb_pfunc2
+        nb_pfunc3 = init_nb_pfunc3
+
+        open(402, form = 'formatted', status = 'old',                &
+              file ='chemnetwork/'//trim(reactfile),iostat=ierr)
+
+        if (ierr /= 0) THEN
+           write(*,*)'Error : cannot open chemical reaction file : chemnetwork/'//trim(reactfile)
+           write(*,*)'It should be in your simulation directory in chemnetwork directory'
+           write(*,*)'   You can change the input chemical reactions file name in'
+           write(*,*)'   callphys.def with:'
+           write(*,*)'   monoreact=filename or bimolreact=filename or quadrareact=filename'
+           write(*,*)'   depending on what chemical reaction type it is'
+           stop
+        end if
+
+        do ilines=1,nlines
+          paramline = ''
+
+          read(402,'(A,A,A)') reactline, prodline, paramline
+
+          ! continue only if it's not a comment line
+          if (reactline(1:1)/='!' .and. reactline(1:1)/='') then
+
+            ! increment number of reactions and init
+            ireact = ireact + 1
+        !!!!!!!!!!!!!!!!!!!!!!!!!!! for fill indice part
+            nrtype = nrtype + 1
+            nindice(:) = 0.0
+            iindice(:) = i_dummy
+        !!!!!!!!!!!!!!!!!!!!!!!!!!! end
+            ! get indice, rate type and parameters
+            if (trim(paramline)=='') then
+              print*, 'Error in reactfile: where are the parameters - line',ilines
+              stop
+            else
+              read(paramline,*) rtype(ireact)
+              if (rtype(ireact)==1) then
+                nb_pfunc1 = nb_pfunc1 + 1
+                read(paramline,*) rtype(ireact), pfunc1_param(nb_pfunc1)
+              else if (rtype(ireact)==0) then
+                nb_hv = nb_hv + 1
+                if (jonline) then
+                  read(paramline,'(I5,A,A)') rtype(ireact), jlabel(nb_hv,1), jparamline(nb_hv)
+                else
+                  read(paramline,*) rtype(ireact), jlabel(nb_hv,1)
+                end if
+              else if (rtype(ireact)==2) then
+                nb_pfunc2 = nb_pfunc2 + 1
+                read(paramline,*) rtype(ireact), pfunc2_param(nb_pfunc2)
+              else if (rtype(ireact)==3) then
+                nb_pfunc3 = nb_pfunc3 + 1
+                read(paramline,*) rtype(ireact), pfunc3_param(nb_pfunc3)
+              end if
+            end if
+        
+            ! get reactants
+            call split_str(reactline,words,nwr,nmax)
+            if (rtype(ireact)==0) jlabel(nb_hv,2) = words(1)
+            ! loop on reactants
+            do iw=1,nwr
+              ! store reactants in variable 'reactants'
+              reactants(ireact,iw) = trim(words(iw))
+              ! check third body and exit reactants loop
+              if (reactants(ireact,iw)=='M') then
+                if (iw==nwr) then
+                  exit
+                else if (iw==nwr-1) then
+                  third_body(ireact) = indexchim(words(iw+1))
+                  exit
+                else
+                  print*, 'Error in reactfile: just only one specie can be after M corresponding to the third body - line',ilines
+                  stop
+                end if
+              end if
+              if (trim(words(iw))/='hv' .and. trim(words(iw))/='dummy') then
+                iwhere = indexchim(words(iw))
+                ! check if species are chemical tracers
+                if (iwhere>nesp) then
+                  print*, 'Error: in ', trim(reactfile)
+                  print*, 'check if the specie', trim(words(iw)),' is include into chemical tracers in traceur.def'
+                  stop
+                end if
+                ! to count the species used in 'reactfile'
+                if (specheck(iwhere)==0) then
+                  specheckr(iwhere) = 1
+                  specheck(iwhere) = 1
+                  nbq = nbq + 1
+                else if (specheckr(iwhere)==0) then
+                  specheckr(iwhere) = 1
+                end if
+        
+        !!!!!!!!!!!!!!!!!!! for fill indice part
+                ! fill stochiometry and indice of rection species depending of reaction type
+                if (trim(typeindice)=='v3') then
+                  nindice(1) = 2.0
+                  iindice(1) = indexchim(words(iw))
+                  if (nwr>3 .or. nwr<2) print*, 'Error in reactfile: wrong number of reactants for v3 reaction line',ilines
+                  if (nwr==2 .and. trim(words(1))/=trim(words(2))) print*, 'Error in reactfile: both reactants should be the same for v3 reaction line',ilines
+                else if (trim(typeindice)=='v4') then
+                  nindice(iw) = 1.0
+                  iindice(iw) = indexchim(words(iw))
+                else if (trim(typeindice)=='vphot') then
+                  nindice(1) = 1.0
+                  if (iw>2) then
+                    print*, 'Something weird in your photolysis reaction'
+                    print*, 'You should have 1 reactants and hv'
+                    print*, 'Reactants are: ',words
+                    stop
+                  end if
+                  iindice(1) = indexchim(words(iw))
+                end if
+        !!!!!!!!!!!!!!!!!!! end
+        
+              end if
+            end do
+
+            ! same as reactants but for the products
+            call split_str(prodline,words,nwp,nmax)
+            do iw=1,nwp
+              spedouble = .false.
+              products(ireact,iw) = trim(words(iw))
+              if (trim(words(iw))/='hv' .and. trim(words(iw))/='dummy' .and. trim(words(iw))/='M') then
+                iwhere = indexchim(words(iw))
+                if (iwhere>nesp) then
+                  print*, 'Error: in ', trim(reactfile)
+                  print*, 'check if the specie', trim(words(iw)),' is include into chemical tracers in traceur.def'
+                  stop
+                end if
+                if (specheck(iwhere)==0) then
+                  specheckp(iwhere) = 1
+                  specheck(iwhere) = 1
+                  nbq = nbq + 1
+                else if (specheckp(iwhere)==0) then
+                  specheckp(iwhere) = 1
+                end if
+        
+        !!!!!!!!!!!!!!!!!!!!!!!!!! for fill indice part
+                if (trim(typeindice)=='v3' .or. trim(typeindice)=='vphot') then
+                    iindice(1+iw) = indexchim(words(iw))
+                  do i=1,iw-1
+                    if (iindice(1+iw)==iindice(1+i)) then
+                      nindice(1+i)  = nindice(1+i) + 1.0
+                      iindice(1+iw) = i_dummy
+                      spedouble = .true.
+                    end if
+                  end do
+                  if (.not. spedouble) nindice(1+iw) = 1.0
+                else if (trim(typeindice)=='v4') then
+                  iindice(2+iw) = indexchim(words(iw))
+                  do i=1,iw-1
+                    if (iindice(2+iw)==iindice(2+i)) then
+                      nindice(2+i)  = nindice(2+i) + 1.0
+                      iindice(2+iw) = i_dummy
+                      spedouble = .true.
+                    end if
+                  end do
+                  if (.not. spedouble) nindice(2+iw) = 1.0
+                end if
+        !!!!!!!!!!!!!!!!!!!!!!!!!!! end
+              else
+                print*, 'Error: no hv, dummy or M in products'
+                stop
+              end if
+            end do
+
+            ! fill indice variables
+            if (trim(typeindice)=='v3') then
+              if (nindice(4)/=0.0) then ! reaction with 3 products
+                if (nindice(3)==0.0) then ! 2 are the same species
+                  indice_3(nrtype) = z3spec(nindice(1), iindice(1), nindice(2), iindice(2), nindice(4), iindice(4))
+                else ! reaction with 3 different products
+                  indice_3(nrtype) = z3spec(nindice(1)/2., iindice(1), nindice(2), iindice(2), 0.0, i_dummy)
+                  nrtype = nrtype + 1
+                  indice_3(nrtype) = z3spec(nindice(1)/2., iindice(1), nindice(3), iindice(3), nindice(4), iindice(4))
+                  three_prod(ireact) = .true.
+                end if
+              else ! reaction with 1 or 2 products
+                indice_3(nrtype) = z3spec(nindice(1), iindice(1), nindice(2), iindice(2), nindice(3), iindice(3))
+              end if
+            else if (trim(typeindice)=='v4') then
+              if (nindice(5)/=0.0) then ! reaction with 3 products
+                if (nindice(4)==0.0) then ! 2 are the same species
+                  indice_4(nrtype) = z4spec(nindice(1), iindice(1), nindice(2), iindice(2), nindice(3), iindice(3), nindice(5), iindice(5))
+                else ! reaction with 3 different products
+                  indice_4(nrtype) = z4spec(nindice(1)/2., iindice(1), nindice(2)/2., iindice(2), nindice(3), iindice(3), nindice(4)/2., iindice(4))
+                  nrtype = nrtype + 1
+                  indice_4(nrtype) = z4spec(nindice(1)/2., iindice(1), nindice(2)/2., iindice(2), nindice(5), iindice(5), nindice(4)/2., iindice(4))
+                  three_prod(ireact) = .true.
+                end if
+              else ! reaction with 1 or 2 products
+                indice_4(nrtype) = z4spec(nindice(1), iindice(1), nindice(2), iindice(2), nindice(3), iindice(3), nindice(4), iindice(4))
+              end if
+            else if (trim(typeindice)=='vphot') then
+              if (nindice(4)/=0.0) then ! reaction with 3 products
+                if (nindice(3)==0.0) then ! 2 are the same species
+                  indice_phot(nrtype) = z3spec(nindice(1), iindice(1), nindice(2), iindice(2), nindice(4), iindice(4))
+                else ! reaction with 3 different products
+                  indice_phot(nrtype) = z3spec(nindice(1)/2., iindice(1), nindice(2), iindice(2), 0.0, i_dummy)
+                  nrtype = nrtype + 1
+                  indice_phot(nrtype) = z3spec(nindice(1)/2., iindice(1), nindice(3), iindice(3), nindice(4), iindice(4))
+                  three_prod(ireact) = .true.
+                end if
+              else ! reaction with 1 or 2 products
+                indice_phot(nrtype) = z3spec(nindice(1), iindice(1), nindice(2), iindice(2), nindice(3), iindice(3))
+              end if
+            end if
+        
+          end if
+
+        end do
+
+        init_nb_pfunc1 = nb_pfunc1
+        init_nb_pfunc2 = nb_pfunc2
+        init_nb_pfunc3 = nb_pfunc3
+
+        close(402)
+
+    end subroutine get_react
+
 end subroutine photochemistry_asis
-

trunk/LMDZ.GENERIC/libf/aeronostd/photolysis_asis.F90

@@ -1 +1 @@
 !==========================================================================
 
-      subroutine photolysis_asis(nlayer, ngrid,                                 &
-                                 lswitch, press, temp, sza, fractcol, tauref,   &
-                                 zmmean, dist_sol, rmco2, rmo3, v_phot)
+      subroutine photolysis_asis(nlayer, ngrid,                    &
+                                 lswitch, press, temp, sza,fractcol, tauref,   &
+                                 zmmean, dist_sol, rmco2, rmo3, rmch4, v_phot, e_phot, nreact, three_prod)
 
 !==========================================================================
@@ -9 +9 @@
       use comcstfi_mod
       use callkeys_mod
+      use types_asis
+      use chimiedata_h
 
       implicit none
 
-#include "chimiedata.h" 
+!#include "chimiedata.h" 
 
 !==========================================================================
@@ -19 +21 @@
        
       integer, intent(in) :: nlayer ! number of atmospheric layers
-      integer,intent(in) :: ngrid   ! number of atmospheric columns
+      integer, intent(in) :: ngrid  ! number of atmospheric columns
       integer :: lswitch            ! interface level between chemistries
       real :: press(nlayer)         ! pressure (hPa)
@@ -26 +28 @@
       real :: fractcol              ! day fraction
       real :: tauref                ! optical depth at 7 hpa
-      real :: zmmean(nlayer)        ! mean molecular mass (g)
+      real :: zmmean(nlayer)        ! mean molecular mass (g/mol)
       real :: dist_sol              ! sun distance (AU) 
       real :: rmco2(nlayer)         ! co2 mixing ratio
       real :: rmo3(nlayer)          ! ozone mixing ratio
+      real :: rmch4(nlayer)         ! ch4 mixing ratio
+      integer, intent(in) :: nreact             ! number of reactions in reactions files
+      logical, intent(in) :: three_prod(nreact) ! if the reaction have a three defferent products egal .true.
 
 !==========================================================================
@@ -36 +41 @@
 
       real (kind = 8), dimension(nlayer,nb_phot_max) :: v_phot
+      real (kind = 8), dimension(nlayer,nb_phot_max) :: e_phot
 
 !==========================================================================
@@ -42 +48 @@
 
       integer :: icol, ij, indsza, indtau, indcol, indozo, indtemp,     &
-                 iozo, isza, itau, it, l
-
-      integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d,        &
-                 j_o3_o, j_h2o, j_hdo, j_h2o2, j_ho2, j_no, j_no2,      &
-                 j_hno3, j_hno4,                                        &
-                 j_ch4_ch3_h, j_ch4_1ch2_h2, j_ch4_3ch2_h_h,            &
-                 j_ch4_ch_h2_h, j_ch3o2h, j_ch2o_hco, j_ch2o_co,        &
-                 j_ch3oh, j_c2h6, j_hcl, j_hocl, j_clo, j_so2, j_so,    &
-                 j_h2s, j_so3
+                 iozo, isza, itau, it, ich4, indch4, l, nb_phot 
 
       real :: col(nlayer)                 ! overhead air column   (molecule cm-2)
       real :: colo3(nlayer)               ! overhead ozone column (molecule cm-2)
-      real :: poids(2,2,2,2,2)            ! 5D interpolation weights 
+      real :: colch4(nlayer)               ! overhead CH4 column (molecule cm-2)
+      real :: tauch4(nlayer)               ! estimation of optical depth by CH4
+      real :: ch4_equ(nlayer)               ! equivalent constant mixing ratio for the same column of CH4
+      real :: poids(2,2,2,2,2,2)          ! 6D interpolation weights 
       real :: tref                        ! temperature  at 1.9 hPa in the gcm (K)
       real :: table_temp(ntemp)           ! temperatures at 1.9 hPa in jmars   (K)
-      real :: cinf, csup, cicol, ciozo, cisza, citemp, citau
+      real :: ch4ref                      ! ch4 mixing ratio at top of the atmosphere
+      real :: cinf, csup, cicol, ciozo, cisza, citemp, citau, cich4
       real :: colo3min, dp, coef
       real :: ratio_o3(nlayer)
       real :: tau
       real :: j(nlayer,nd)
+      real :: e(nlayer,nd)
 
 !==========================================================================
@@ -73 +76 @@
 !==========================================================================
       
-      table_temp(1) = 226.2
-      table_temp(2) = 206.2
-      table_temp(3) = 186.2
-      table_temp(4) = 169.8
+!      table_temp(1) = 226.2
+!      table_temp(2) = 206.2
+!      table_temp(3) = 186.2
+!      table_temp(4) = 169.8
+
+!      table_temp(2) = 186.2
+      table_temp(1) = 176.2
 
 !==========================================================================
@@ -90 +96 @@
          end if
       end do
-      cisza = (sza - szatab(indsza))  &
+
+      if(nsza.eq.1) then
+        cisza = 0.
+        indsza=1
+      else
+        cisza = (sza - szatab(indsza))  &
              /(szatab(indsza + 1) - szatab(indsza))
+      endif
 
 !==========================================================================
@@ -108 +120 @@
          end if
       end do
-      citau = (tau - tautab(indtau))     &
+
+      if(ntau.eq.1) then
+        citau=0.
+        indtau=1
+      else
+        citau = (tau - tautab(indtau))     &
              /(tautab(indtau + 1) - tautab(indtau))
-
-!==========================================================================
-!     co2 and ozone columns
+      endif
+
+
+
+!==========================================================================
+!     air and ozone columns
 !==========================================================================
 
 !     co2 column at model top (molecule.cm-2)
 
-      col(lswitch-1) = 6.022e22*rmco2(lswitch-1)*press(lswitch-1)*100.  &
-                       /(zmmean(lswitch-1)*g)
+!      col(lswitch-1) = 6.022e22*rmco2(lswitch-1)*press(lswitch-1)*100.  &
+!                       /(zmmean(lswitch-1)*g)
+      col(lswitch-1) = 6.022e22*press(lswitch-1)*100./(zmmean(lswitch-1)*g)
+
 
 !     ozone column at model top
 
       colo3(lswitch-1) = 0.
-
 !     co2 and ozone columns for other levels (molecule.cm-2)
 
       do l = lswitch-2,1,-1
          dp = (press(l) - press(l+1))*100.
-         col(l) = col(l+1) + (rmco2(l+1) + rmco2(l))*0.5   &
-                             *6.022e22*dp/(zmmean(l)*g)
+!         col(l) = col(l+1) + (rmco2(l+1) + rmco2(l))*0.5   &
+!                             *6.022e22*dp/(zmmean(l)*g)
+         col(l) = col(l+1) +  6.022e22*dp/(zmmean(l)*g)
          col(l) = min(col(l), colairtab(1))
          colo3(l) = colo3(l+1) + (rmo3(l+1) + rmo3(l))*0.5 &
                                  *6.022e22*dp/(zmmean(l)*g)
-      end do
-
-!     ratio ozone column/minimal theoretical column (0.1 micron-atm) 
-
-!     ro3 = 7.171e-10 is the o3 mixing ratio for a uniform
-!     profile giving a column 0.1 micron-atmosphere at
-!     a surface pressure of 10 hpa.
+
+      end do
+
+!     ratio ozone column/minimal theoretical column (10 micron-atm) 
+
+!     ro3 = 1.227e-10 is the o3 mixing ratio for a uniform
+!     profile giving a column 10 micron-atmosphere at
+!     a surface pressure of 1 bar.
 
       do l = 1,lswitch-1
-         colo3min    = col(l)*7.171e-10
+!         colo3min    = col(l)*7.171e-10
+         colo3min    = col(l)*1.227e-10*(g/9.81)*(mugaz/28)
          ratio_o3(l) = colo3(l)/colo3min
-         ratio_o3(l) = min(ratio_o3(l), table_ozo(nozo)*10.)
-         ratio_o3(l) = max(ratio_o3(l), 1.)
-      end do
+         ratio_o3(l) = min(ratio_o3(l), table_ozo(nozo))
+         ratio_o3(l) = max(ratio_o3(l), 0.)
+      end do
+
+!      print*,'co3(1)=',colo3(1)
+!      print*,'col(1)=',col(1)
+!      print*,'ratio_o3(1)=',ratio_o3(1)
+!      print*,'maxval(ratio_o3)=',maxval(ratio_o3(:))
+!      print*,'maxval(ozo)=',table_ozo(nozo)/10.
+
+!==========================================================================
+!     ch4 dependence
+!==========================================================================
+
+!     1) search for temperature at 1.9 hPa (tref): vertical
+!     interpolation
+
+      ch4ref = rmch4(lswitch-2)
+      colch4(lswitch-1) = 0.
+      ch4_equ(lswitch-1) = 0.
+      do l = lswitch-2,1,-1
+         dp = (press(l) - press(l+1))*100.
+         colch4(l) = colch4(l+1) + (rmch4(l+1) + rmch4(l))*0.5 &
+                                 *6.022e22*dp/(zmmean(l)*g)
+         ch4_equ(l)=colch4(l)/col(l)
+!         tauch4(l)=1.8e-21*colch4(l)
+!         if(tauch4(l).ge.1.0) exit
+      end do
+!      ch4ref = (rmch4(l+1)*(tauch4(l)-1)+rmch4(l)*(1-tauch4(l+1))) &
+!                                 /(tauch4(l)-tauch4(l+1))
+
+!     2) interpolation in CH4
+
+!      ch4ref = min(ch4ref,table_ch4(nch4))
+!      ch4ref = max(ch4ref,table_ch4(1))
+
+!      indch4 = nch4 - 1
+!      do ich4 = 1,nch4
+!         if (table_ch4(ich4) >= ch4ref) then
+!            indch4 = ich4 - 1
+!            indch4 = max(indch4, 1)
+!            exit
+!         end if
+!      end do
+!      cich4 = (ch4ref - table_ch4(indch4))     &
+!             /(table_ch4(indch4 + 1) - table_ch4(indch4))
+
+
 
 !==========================================================================
@@ -169 +238 @@
       tref = max(tref,table_temp(ntemp))
 
-      do it = 2, ntemp
-         if (table_temp(it) <= tref) then
-            citemp = (log(tref) - log(table_temp(it)))              &
-                    /(log(table_temp(it-1)) - log(table_temp(it)))
-            indtemp = it - 1
-            exit
-         end if
-      end do
+      if(ntemp.eq.1) then
+        citemp = 1.
+        indtemp = 1
+      else
+        do it = 2, ntemp
+           if (table_temp(it) <= tref) then
+              citemp = (log(tref) - log(table_temp(it)))              &
+                      /(log(table_temp(it-1)) - log(table_temp(it)))
+              indtemp = it - 1
+              exit
+           end if
+        end do
+      endif
+
+
 
 !==========================================================================
@@ -200 +276 @@
          indozo = nozo - 1
          do iozo = 1,nozo
-            if (table_ozo(iozo)*10. >= ratio_o3(l)) then
+            if (table_ozo(iozo) >= ratio_o3(l)) then
                indozo = iozo - 1
                indozo = max(indozo, 1)
@@ -206 +282 @@
             end if
          end do
-         ciozo = (ratio_o3(l) - table_ozo(indozo)*10.)             &
-                /(table_ozo(indozo + 1)*10. - table_ozo(indozo)*10.)
-
-!     4-dimensional interpolation weights
-
-!     poids(temp,sza,co2,o3,tau)
-
-         poids(1,1,1,1,1) = citemp*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)
-         poids(1,1,1,2,1) = citemp*(1.-cisza)*cicol*ciozo*(1.-citau)
-         poids(1,1,2,1,1) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)
-         poids(1,1,2,2,1) = citemp*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)
-         poids(1,2,1,1,1) = citemp*cisza*cicol*(1.-ciozo)*(1.-citau)
-         poids(1,2,1,2,1) = citemp*cisza*cicol*ciozo*(1.-citau)
-         poids(1,2,2,1,1) = citemp*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)
-         poids(1,2,2,2,1) = citemp*cisza*(1.-cicol)*ciozo*(1.-citau)
-         poids(2,1,1,1,1) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)
-         poids(2,1,1,2,1) = (1.-citemp)*(1.-cisza)*cicol*ciozo*(1.-citau)
-         poids(2,1,2,1,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)
-         poids(2,1,2,2,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)
-         poids(2,2,1,1,1) = (1.-citemp)*cisza*cicol*(1.-ciozo)*(1.-citau)
-         poids(2,2,1,2,1) = (1.-citemp)*cisza*cicol*ciozo*(1.-citau)
-         poids(2,2,2,1,1) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)
-         poids(2,2,2,2,1) = (1.-citemp)*cisza*(1.-cicol)*ciozo*(1.-citau)
-!
-         poids(1,1,1,1,2) = citemp*(1.-cisza)*cicol*(1.-ciozo)*citau
-         poids(1,1,1,2,2) = citemp*(1.-cisza)*cicol*ciozo*citau
-         poids(1,1,2,1,2) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau
-         poids(1,1,2,2,2) = citemp*(1.-cisza)*(1.-cicol)*ciozo*citau
-         poids(1,2,1,1,2) = citemp*cisza*cicol*(1.-ciozo)*citau
-         poids(1,2,1,2,2) = citemp*cisza*cicol*ciozo*citau
-         poids(1,2,2,1,2) = citemp*cisza*(1.-cicol)*(1.-ciozo)*citau
-         poids(1,2,2,2,2) = citemp*cisza*(1.-cicol)*ciozo*citau
-         poids(2,1,1,1,2) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*citau
-         poids(2,1,1,2,2) = (1.-citemp)*(1.-cisza)*cicol*ciozo*citau
-         poids(2,1,2,1,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau
-         poids(2,1,2,2,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*citau
-         poids(2,2,1,1,2) = (1.-citemp)*cisza*cicol*(1.-ciozo)*citau
-         poids(2,2,1,2,2) = (1.-citemp)*cisza*cicol*ciozo*citau
-         poids(2,2,2,1,2) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*citau
-         poids(2,2,2,2,2) = (1.-citemp)*cisza*(1.-cicol)*ciozo*citau
+
+      if(nozo.eq.1) then
+         ciozo = 0.
+      else
+         ciozo = (ratio_o3(l) - table_ozo(indozo))             &
+                /(table_ozo(indozo + 1) - table_ozo(indozo))
+      endif
+
+!     2) interpolation in CH4
+
+      ch4ref = min(ch4_equ(l),table_ch4(nch4))
+      ch4ref = max(ch4ref,table_ch4(1))
+
+      indch4 = nch4 - 1
+      do ich4 = 1,nch4
+         if (table_ch4(ich4) >= ch4ref) then
+            indch4 = ich4 - 1
+            indch4 = max(indch4, 1)
+            exit
+         end if
+      end do
+      if(nch4.eq.1) then
+        cich4=0.
+        indch4=1
+      else
+        cich4 = (ch4ref - table_ch4(indch4))     &
+               /(table_ch4(indch4 + 1) - table_ch4(indch4))
+      endif
+
+!     5-dimensional interpolation weights
+
+!     poids(temp,sza,co2,o3,tau,ch4)
+
+         poids(1,1,1,1,1,1) = citemp*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)*(1.-cich4)
+         poids(1,1,1,2,1,1) = citemp*(1.-cisza)*cicol*ciozo*(1.-citau)*(1.-cich4)       
+         poids(1,1,2,1,1,1) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)*(1.-cich4)
+         poids(1,1,2,2,1,1) = citemp*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)*(1.-cich4) 
+         poids(1,2,1,1,1,1) = citemp*cisza*cicol*(1.-ciozo)*(1.-citau)*(1.-cich4)               
+         poids(1,2,1,2,1,1) = citemp*cisza*cicol*ciozo*(1.-citau)*(1.-cich4)                    
+         poids(1,2,2,1,1,1) = citemp*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)*(1.-cich4)          
+         poids(1,2,2,2,1,1) = citemp*cisza*(1.-cicol)*ciozo*(1.-citau)*(1.-cich4)               
+!         poids(2,1,1,1,1,1) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)*(1.-cich4)     
+!         poids(2,1,1,2,1,1) = (1.-citemp)*(1.-cisza)*cicol*ciozo*(1.-citau)*(1.-cich4)          
+!         poids(2,1,2,1,1,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)*(1.-cich4)
+!         poids(2,1,2,2,1,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)*(1.-cich4)     
+!         poids(2,2,1,1,1,1) = (1.-citemp)*cisza*cicol*(1.-ciozo)*(1.-citau)*(1.-cich4)          
+!         poids(2,2,1,2,1,1) = (1.-citemp)*cisza*cicol*ciozo*(1.-citau)*(1.-cich4)               
+!         poids(2,2,2,1,1,1) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)*(1.-cich4)     
+!         poids(2,2,2,2,1,1) = (1.-citemp)*cisza*(1.-cicol)*ciozo*(1.-citau)*(1.-cich4)          
+!!
+!         poids(1,1,1,1,2,1) = citemp*(1.-cisza)*cicol*(1.-ciozo)*citau*(1.-cich4)               
+!         poids(1,1,1,2,2,1) = citemp*(1.-cisza)*cicol*ciozo*citau*(1.-cich4)                    
+!         poids(1,1,2,1,2,1) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau*(1.-cich4)          
+!         poids(1,1,2,2,2,1) = citemp*(1.-cisza)*(1.-cicol)*ciozo*citau*(1.-cich4)               
+!         poids(1,2,1,1,2,1) = citemp*cisza*cicol*(1.-ciozo)*citau*(1.-cich4)                    
+!         poids(1,2,1,2,2,1) = citemp*cisza*cicol*ciozo*citau*(1.-cich4)                         
+!         poids(1,2,2,1,2,1) = citemp*cisza*(1.-cicol)*(1.-ciozo)*citau*(1.-cich4)               
+!         poids(1,2,2,2,2,1) = citemp*cisza*(1.-cicol)*ciozo*citau*(1.-cich4)                    
+!         poids(2,1,1,1,2,1) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*citau*(1.-cich4)          
+!         poids(2,1,1,2,2,1) = (1.-citemp)*(1.-cisza)*cicol*ciozo*citau*(1.-cich4)               
+!         poids(2,1,2,1,2,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau*(1.-cich4)     
+!         poids(2,1,2,2,2,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*citau*(1.-cich4)          
+!         poids(2,2,1,1,2,1) = (1.-citemp)*cisza*cicol*(1.-ciozo)*citau*(1.-cich4)               
+!         poids(2,2,1,2,2,1) = (1.-citemp)*cisza*cicol*ciozo*citau*(1.-cich4)                    
+!         poids(2,2,2,1,2,1) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*citau*(1.-cich4)          
+!         poids(2,2,2,2,2,1) = (1.-citemp)*cisza*(1.-cicol)*ciozo*citau*(1.-cich4)               
+!      
+         poids(1,1,1,1,1,2) = citemp*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)*cich4               
+         poids(1,1,1,2,1,2) = citemp*(1.-cisza)*cicol*ciozo*(1.-citau)*cich4                    
+         poids(1,1,2,1,1,2) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)*cich4          
+         poids(1,1,2,2,1,2) = citemp*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)*cich4               
+         poids(1,2,1,1,1,2) = citemp*cisza*cicol*(1.-ciozo)*(1.-citau)*cich4                    
+         poids(1,2,1,2,1,2) = citemp*cisza*cicol*ciozo*(1.-citau)*cich4                         
+         poids(1,2,2,1,1,2) = citemp*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)*cich4               
+         poids(1,2,2,2,1,2) = citemp*cisza*(1.-cicol)*ciozo*(1.-citau)*cich4                    
+!         poids(2,1,1,1,1,2) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)*cich4          
+!         poids(2,1,1,2,1,2) = (1.-citemp)*(1.-cisza)*cicol*ciozo*(1.-citau)*cich4               
+!         poids(2,1,2,1,1,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)*cich4     
+!         poids(2,1,2,2,1,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)*cich4          
+!         poids(2,2,1,1,1,2) = (1.-citemp)*cisza*cicol*(1.-ciozo)*(1.-citau)*cich4               
+!         poids(2,2,1,2,1,2) = (1.-citemp)*cisza*cicol*ciozo*(1.-citau)*cich4                    
+!         poids(2,2,2,1,1,2) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)*cich4          
+!         poids(2,2,2,2,1,2) = (1.-citemp)*cisza*(1.-cicol)*ciozo*(1.-citau)*cich4               
+!!
+!         poids(1,1,1,1,2,2) = citemp*(1.-cisza)*cicol*(1.-ciozo)*citau*cich4                    
+!         poids(1,1,1,2,2,2) = citemp*(1.-cisza)*cicol*ciozo*citau*cich4                         
+!         poids(1,1,2,1,2,2) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau*cich4               
+!         poids(1,1,2,2,2,2) = citemp*(1.-cisza)*(1.-cicol)*ciozo*citau*cich4                    
+!         poids(1,2,1,1,2,2) = citemp*cisza*cicol*(1.-ciozo)*citau*cich4                         
+!         poids(1,2,1,2,2,2) = citemp*cisza*cicol*ciozo*citau*cich4                              
+!         poids(1,2,2,1,2,2) = citemp*cisza*(1.-cicol)*(1.-ciozo)*citau*cich4                    
+!         poids(1,2,2,2,2,2) = citemp*cisza*(1.-cicol)*ciozo*citau*cich4                         
+!         poids(2,1,1,1,2,2) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*citau*cich4               
+!         poids(2,1,1,2,2,2) = (1.-citemp)*(1.-cisza)*cicol*ciozo*citau*cich4                    
+!         poids(2,1,2,1,2,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau*cich4          
+!         poids(2,1,2,2,2,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*citau*cich4               
+!         poids(2,2,1,1,2,2) = (1.-citemp)*cisza*cicol*(1.-ciozo)*citau*cich4                    
+!         poids(2,2,1,2,2,2) = (1.-citemp)*cisza*cicol*ciozo*citau*cich4                         
+!         poids(2,2,2,1,2,2) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*citau*cich4               
+!         poids(2,2,2,2,2,2) = (1.-citemp)*cisza*(1.-cicol)*ciozo*citau*cich4                    
+
+
+
+
+
+
+
 
 !     4-dimensional interpolation in the lookup table
 
          do ij = 1,nd
-            j(l,ij) =                                                                & 
-            poids(1,1,1,1,1)*jphot(indtemp,indsza,indcol,indozo,indtau,ij)           &
-          + poids(1,1,1,2,1)*jphot(indtemp,indsza,indcol,indozo+1,indtau,ij)         &
-          + poids(1,1,2,1,1)*jphot(indtemp,indsza,indcol+1,indozo,indtau,ij)         &
-          + poids(1,1,2,2,1)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau,ij)       &
-          + poids(1,2,1,1,1)*jphot(indtemp,indsza+1,indcol,indozo,indtau,ij)         &
-          + poids(1,2,1,2,1)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau,ij)       &
-          + poids(1,2,2,1,1)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau,ij)       &
-          + poids(1,2,2,2,1)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau,ij)     &
-          + poids(2,1,1,1,1)*jphot(indtemp+1,indsza,indcol,indozo,indtau,ij)         &
-          + poids(2,1,1,2,1)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau,ij)       &
-          + poids(2,1,2,1,1)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau,ij)       &
-          + poids(2,1,2,2,1)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau,ij)     &
-          + poids(2,2,1,1,1)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau,ij)       &
-          + poids(2,2,1,2,1)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau,ij)     &
-          + poids(2,2,2,1,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau,ij)     &
-          + poids(2,2,2,2,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau,ij)   &
-!
-          + poids(1,1,1,1,2)*jphot(indtemp,indsza,indcol,indozo,indtau+1,ij)         &
-          + poids(1,1,1,2,2)*jphot(indtemp,indsza,indcol,indozo+1,indtau+1,ij)       &
-          + poids(1,1,2,1,2)*jphot(indtemp,indsza,indcol+1,indozo,indtau+1,ij)       &
-          + poids(1,1,2,2,2)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau+1,ij)     &
-          + poids(1,2,1,1,2)*jphot(indtemp,indsza+1,indcol,indozo,indtau+1,ij)       &
-          + poids(1,2,1,2,2)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau+1,ij)     &
-          + poids(1,2,2,1,2)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau+1,ij)     &
-          + poids(1,2,2,2,2)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau+1,ij)   &
-          + poids(2,1,1,1,2)*jphot(indtemp+1,indsza,indcol,indozo,indtau+1,ij)       &
-          + poids(2,1,1,2,2)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau+1,ij)     &
-          + poids(2,1,2,1,2)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau+1,ij)     &
-          + poids(2,1,2,2,2)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau+1,ij)   &
-          + poids(2,2,1,1,2)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau+1,ij)     &
-          + poids(2,2,1,2,2)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau+1,ij)   &
-          + poids(2,2,2,1,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau+1,ij)   &
-          + poids(2,2,2,2,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau+1,ij)
+            j(l,ij) =                                                                         & 
+            poids(1,1,1,1,1,1)*jphot(indtemp,indsza,indcol,indozo,indtau,indch4,ij)           &
+          + poids(1,1,1,2,1,1)*jphot(indtemp,indsza,indcol,indozo+1,indtau,indch4,ij)         &
+          + poids(1,1,2,1,1,1)*jphot(indtemp,indsza,indcol+1,indozo,indtau,indch4,ij)         &
+          + poids(1,1,2,2,1,1)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau,indch4,ij)       &
+          + poids(1,2,1,1,1,1)*jphot(indtemp,indsza+1,indcol,indozo,indtau,indch4,ij)         &
+          + poids(1,2,1,2,1,1)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau,indch4,ij)       &
+          + poids(1,2,2,1,1,1)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau,indch4,ij)       &
+          + poids(1,2,2,2,1,1)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau,indch4,ij)     &
+!          + poids(2,1,1,1,1,1)*jphot(indtemp+1,indsza,indcol,indozo,indtau,indch4,ij)         &
+!          + poids(2,1,1,2,1,1)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau,indch4,ij)       &
+!          + poids(2,1,2,1,1,1)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau,indch4,ij)       &
+!          + poids(2,1,2,2,1,1)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau,indch4,ij)     &
+!          + poids(2,2,1,1,1,1)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau,indch4,ij)       &
+!          + poids(2,2,1,2,1,1)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau,indch4,ij)     &
+!          + poids(2,2,2,1,1,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau,indch4,ij)     &
+!          + poids(2,2,2,2,1,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau,indch4,ij)   &
+!!
+!          + poids(1,1,1,1,2,1)*jphot(indtemp,indsza,indcol,indozo,indtau+1,indch4,ij)         &
+!          + poids(1,1,1,2,2,1)*jphot(indtemp,indsza,indcol,indozo+1,indtau+1,indch4,ij)       &
+!          + poids(1,1,2,1,2,1)*jphot(indtemp,indsza,indcol+1,indozo,indtau+1,indch4,ij)       &
+!          + poids(1,1,2,2,2,1)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau+1,indch4,ij)     &
+!          + poids(1,2,1,1,2,1)*jphot(indtemp,indsza+1,indcol,indozo,indtau+1,indch4,ij)       &
+!          + poids(1,2,1,2,2,1)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau+1,indch4,ij)     &
+!          + poids(1,2,2,1,2,1)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau+1,indch4,ij)     &
+!          + poids(1,2,2,2,2,1)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau+1,indch4,ij)   &
+!          + poids(2,1,1,1,2,1)*jphot(indtemp+1,indsza,indcol,indozo,indtau+1,indch4,ij)       &
+!          + poids(2,1,1,2,2,1)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau+1,indch4,ij)     &
+!          + poids(2,1,2,1,2,1)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau+1,indch4,ij)     &
+!          + poids(2,1,2,2,2,1)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau+1,indch4,ij)   &
+!          + poids(2,2,1,1,2,1)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau+1,indch4,ij)     &
+!          + poids(2,2,1,2,2,1)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau+1,indch4,ij)   &
+!          + poids(2,2,2,1,2,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau+1,indch4,ij)   &
+!          + poids(2,2,2,2,2,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau+1,indch4,ij)
+! CH4
+          + poids(1,1,1,1,1,2)*jphot(indtemp,indsza,indcol,indozo,indtau,indch4+1,ij)         &
+          + poids(1,1,1,2,1,2)*jphot(indtemp,indsza,indcol,indozo+1,indtau,indch4+1,ij)       &
+          + poids(1,1,2,1,1,2)*jphot(indtemp,indsza,indcol+1,indozo,indtau,indch4+1,ij)       &
+          + poids(1,1,2,2,1,2)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau,indch4+1,ij)     &
+          + poids(1,2,1,1,1,2)*jphot(indtemp,indsza+1,indcol,indozo,indtau,indch4+1,ij)       &
+          + poids(1,2,1,2,1,2)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau,indch4+1,ij)     &
+          + poids(1,2,2,1,1,2)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau,indch4+1,ij)     &
+          + poids(1,2,2,2,1,2)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau,indch4+1,ij) 
+!          + poids(2,1,1,1,1,2)*jphot(indtemp+1,indsza,indcol,indozo,indtau,indch4+1,ij) &
+!          + poids(2,1,1,2,1,2)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau,indch4+1,ij) &
+!          + poids(2,1,2,1,1,2)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau,indch4+1,ij) &
+!          + poids(2,1,2,2,1,2)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau,indch4+1,ij) &
+!          + poids(2,2,1,1,1,2)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau,indch4+1,ij) &
+!          + poids(2,2,1,2,1,2)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau,indch4+1,ij) &
+!          + poids(2,2,2,1,1,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau,indch4+1,ij) &
+!          + poids(2,2,2,2,1,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau,indch4+1,ij) &
+!!
+!          + poids(1,1,1,1,2,2)*jphot(indtemp,indsza,indcol,indozo,indtau+1,indch4+1,ij) &
+!          + poids(1,1,1,2,2,2)*jphot(indtemp,indsza,indcol,indozo+1,indtau+1,indch4+1,ij) &
+!          + poids(1,1,2,1,2,2)*jphot(indtemp,indsza,indcol+1,indozo,indtau+1,indch4+1,ij) &
+!          + poids(1,1,2,2,2,2)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau+1,indch4+1,ij) &
+!          + poids(1,2,1,1,2,2)*jphot(indtemp,indsza+1,indcol,indozo,indtau+1,indch4+1,ij) &
+!          + poids(1,2,1,2,2,2)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau+1,indch4+1,ij) &
+!          + poids(1,2,2,1,2,2)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau+1,indch4+1,ij) &
+!          + poids(1,2,2,2,2,2)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau+1,indch4+1,ij) &
+!          + poids(2,1,1,1,2,2)*jphot(indtemp+1,indsza,indcol,indozo,indtau+1,indch4+1,ij) &
+!          + poids(2,1,1,2,2,2)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau+1,indch4+1,ij) &
+!          + poids(2,1,2,1,2,2)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau+1,indch4+1,ij) &
+!          + poids(2,1,2,2,2,2)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau+1,indch4+1,ij) &
+!          + poids(2,2,1,1,2,2)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau+1,indch4+1,ij) &
+!          + poids(2,2,1,2,2,2)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau+1,indch4+1,ij) &
+!          + poids(2,2,2,1,2,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau+1,indch4+1,ij) &
+!          + poids(2,2,2,2,2,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau+1,indch4+1,ij)
+
          end do
+
+      if (photoheat) then
+         do ij = 1,nd
+            e(l,ij) =                                                                         & 
+            poids(1,1,1,1,1,1)*ephot(indtemp,indsza,indcol,indozo,indtau,indch4,ij)           &
+          + poids(1,1,1,2,1,1)*ephot(indtemp,indsza,indcol,indozo+1,indtau,indch4,ij)         &
+          + poids(1,1,2,1,1,1)*ephot(indtemp,indsza,indcol+1,indozo,indtau,indch4,ij)         &
+          + poids(1,1,2,2,1,1)*ephot(indtemp,indsza,indcol+1,indozo+1,indtau,indch4,ij)       &
+          + poids(1,2,1,1,1,1)*ephot(indtemp,indsza+1,indcol,indozo,indtau,indch4,ij)         &
+          + poids(1,2,1,2,1,1)*ephot(indtemp,indsza+1,indcol,indozo+1,indtau,indch4,ij)       &
+          + poids(1,2,2,1,1,1)*ephot(indtemp,indsza+1,indcol+1,indozo,indtau,indch4,ij)       &
+          + poids(1,2,2,2,1,1)*ephot(indtemp,indsza+1,indcol+1,indozo+1,indtau,indch4,ij)     &
+!          + poids(2,1,1,1,1,1)*ephot(indtemp+1,indsza,indcol,indozo,indtau,indch4,ij)         &
+!          + poids(2,1,1,2,1,1)*ephot(indtemp+1,indsza,indcol,indozo+1,indtau,indch4,ij)       &
+!          + poids(2,1,2,1,1,1)*ephot(indtemp+1,indsza,indcol+1,indozo,indtau,indch4,ij)       &
+!          + poids(2,1,2,2,1,1)*ephot(indtemp+1,indsza,indcol+1,indozo+1,indtau,indch4,ij)     &
+!          + poids(2,2,1,1,1,1)*ephot(indtemp+1,indsza+1,indcol,indozo,indtau,indch4,ij)       &
+!          + poids(2,2,1,2,1,1)*ephot(indtemp+1,indsza+1,indcol,indozo+1,indtau,indch4,ij)     &
+!          + poids(2,2,2,1,1,1)*ephot(indtemp+1,indsza+1,indcol+1,indozo,indtau,indch4,ij)     &
+!          + poids(2,2,2,2,1,1)*ephot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau,indch4,ij)   &
+!!
+!          + poids(1,1,1,1,2,1)*ephot(indtemp,indsza,indcol,indozo,indtau+1,indch4,ij)         &
+!          + poids(1,1,1,2,2,1)*ephot(indtemp,indsza,indcol,indozo+1,indtau+1,indch4,ij)       &
+!          + poids(1,1,2,1,2,1)*ephot(indtemp,indsza,indcol+1,indozo,indtau+1,indch4,ij)       &
+!          + poids(1,1,2,2,2,1)*ephot(indtemp,indsza,indcol+1,indozo+1,indtau+1,indch4,ij)     &
+!          + poids(1,2,1,1,2,1)*ephot(indtemp,indsza+1,indcol,indozo,indtau+1,indch4,ij)       &
+!          + poids(1,2,1,2,2,1)*ephot(indtemp,indsza+1,indcol,indozo+1,indtau+1,indch4,ij)     &
+!          + poids(1,2,2,1,2,1)*ephot(indtemp,indsza+1,indcol+1,indozo,indtau+1,indch4,ij)     &
+!          + poids(1,2,2,2,2,1)*ephot(indtemp,indsza+1,indcol+1,indozo+1,indtau+1,indch4,ij)   &
+!          + poids(2,1,1,1,2,1)*ephot(indtemp+1,indsza,indcol,indozo,indtau+1,indch4,ij)       &
+!          + poids(2,1,1,2,2,1)*ephot(indtemp+1,indsza,indcol,indozo+1,indtau+1,indch4,ij)     &
+!          + poids(2,1,2,1,2,1)*ephot(indtemp+1,indsza,indcol+1,indozo,indtau+1,indch4,ij)     &
+!          + poids(2,1,2,2,2,1)*ephot(indtemp+1,indsza,indcol+1,indozo+1,indtau+1,indch4,ij)   &
+!          + poids(2,2,1,1,2,1)*ephot(indtemp+1,indsza+1,indcol,indozo,indtau+1,indch4,ij)     &
+!          + poids(2,2,1,2,2,1)*ephot(indtemp+1,indsza+1,indcol,indozo+1,indtau+1,indch4,ij)   &
+!          + poids(2,2,2,1,2,1)*ephot(indtemp+1,indsza+1,indcol+1,indozo,indtau+1,indch4,ij)   &
+!          + poids(2,2,2,2,2,1)*ephot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau+1,indch4,ij)
+! CH4
+          + poids(1,1,1,1,1,2)*ephot(indtemp,indsza,indcol,indozo,indtau,indch4+1,ij)         &
+          + poids(1,1,1,2,1,2)*ephot(indtemp,indsza,indcol,indozo+1,indtau,indch4+1,ij)       &
+          + poids(1,1,2,1,1,2)*ephot(indtemp,indsza,indcol+1,indozo,indtau,indch4+1,ij)       &
+          + poids(1,1,2,2,1,2)*ephot(indtemp,indsza,indcol+1,indozo+1,indtau,indch4+1,ij)     &
+          + poids(1,2,1,1,1,2)*ephot(indtemp,indsza+1,indcol,indozo,indtau,indch4+1,ij)       &
+          + poids(1,2,1,2,1,2)*ephot(indtemp,indsza+1,indcol,indozo+1,indtau,indch4+1,ij)     &
+          + poids(1,2,2,1,1,2)*ephot(indtemp,indsza+1,indcol+1,indozo,indtau,indch4+1,ij)     &
+          + poids(1,2,2,2,1,2)*ephot(indtemp,indsza+1,indcol+1,indozo+1,indtau,indch4+1,ij) 
+!          + poids(2,1,1,1,1,2)*ephot(indtemp+1,indsza,indcol,indozo,indtau,indch4+1,ij) &
+!          + poids(2,1,1,2,1,2)*ephot(indtemp+1,indsza,indcol,indozo+1,indtau,indch4+1,ij) &
+!          + poids(2,1,2,1,1,2)*ephot(indtemp+1,indsza,indcol+1,indozo,indtau,indch4+1,ij) &
+!          + poids(2,1,2,2,1,2)*ephot(indtemp+1,indsza,indcol+1,indozo+1,indtau,indch4+1,ij) &
+!          + poids(2,2,1,1,1,2)*ephot(indtemp+1,indsza+1,indcol,indozo,indtau,indch4+1,ij) &
+!          + poids(2,2,1,2,1,2)*ephot(indtemp+1,indsza+1,indcol,indozo+1,indtau,indch4+1,ij) &
+!          + poids(2,2,2,1,1,2)*ephot(indtemp+1,indsza+1,indcol+1,indozo,indtau,indch4+1,ij) &
+!          + poids(2,2,2,2,1,2)*ephot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau,indch4+1,ij) &
+!!
+!          + poids(1,1,1,1,2,2)*ephot(indtemp,indsza,indcol,indozo,indtau+1,indch4+1,ij) &
+!          + poids(1,1,1,2,2,2)*ephot(indtemp,indsza,indcol,indozo+1,indtau+1,indch4+1,ij) &
+!          + poids(1,1,2,1,2,2)*ephot(indtemp,indsza,indcol+1,indozo,indtau+1,indch4+1,ij) &
+!          + poids(1,1,2,2,2,2)*ephot(indtemp,indsza,indcol+1,indozo+1,indtau+1,indch4+1,ij) &
+!          + poids(1,2,1,1,2,2)*ephot(indtemp,indsza+1,indcol,indozo,indtau+1,indch4+1,ij) &
+!          + poids(1,2,1,2,2,2)*ephot(indtemp,indsza+1,indcol,indozo+1,indtau+1,indch4+1,ij) &
+!          + poids(1,2,2,1,2,2)*ephot(indtemp,indsza+1,indcol+1,indozo,indtau+1,indch4+1,ij) &
+!          + poids(1,2,2,2,2,2)*ephot(indtemp,indsza+1,indcol+1,indozo+1,indtau+1,indch4+1,ij) &
+!          + poids(2,1,1,1,2,2)*ephot(indtemp+1,indsza,indcol,indozo,indtau+1,indch4+1,ij) &
+!          + poids(2,1,1,2,2,2)*ephot(indtemp+1,indsza,indcol,indozo+1,indtau+1,indch4+1,ij) &
+!          + poids(2,1,2,1,2,2)*ephot(indtemp+1,indsza,indcol+1,indozo,indtau+1,indch4+1,ij) &
+!          + poids(2,1,2,2,2,2)*ephot(indtemp+1,indsza,indcol+1,indozo+1,indtau+1,indch4+1,ij) &
+!          + poids(2,2,1,1,2,2)*ephot(indtemp+1,indsza+1,indcol,indozo,indtau+1,indch4+1,ij) &
+!          + poids(2,2,1,2,2,2)*ephot(indtemp+1,indsza+1,indcol,indozo+1,indtau+1,indch4+1,ij) &
+!          + poids(2,2,2,1,2,2)*ephot(indtemp+1,indsza+1,indcol+1,indozo,indtau+1,indch4+1,ij) &
+!          + poids(2,2,2,2,2,2)*ephot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau+1,indch4+1,ij)
+
+         end do
+       end if
 
 !     correction for sun distance
@@ -291 +543 @@
 !            j(l,ij) = j(l,ij)*(1.52/dist_sol)**2.
             j(l,ij) = j(l,ij)*(1.0/dist_sol)**2.
-
+            if (photoheat) e(l,ij) = e(l,ij)*(1.0/dist_sol)**2.
+  
             ! Only during daylight.
             if((ngrid.eq.1))then
-                  j(l,ij)= j(l,ij)* 0.25 ! globally averaged = divide by 4
+                  j(l,ij)= j(l,ij)* 0.25 / cos(sza*pi/180.) ! globally averaged = divide by 4
+                  if (photoheat) e(l,ij)= e(l,ij)* 0.25 / cos(sza*pi/180.) ! globally averaged = divide by 4
             elseif(diurnal .eqv. .false.) then
                   j(l,ij)= j(l,ij) * fractcol
+                  if (photoheat) e(l,ij)= e(l,ij) * fractcol
             endif
+
+
          end do 
 
+
+
 !==========================================================================
 !     end of loop over vertical levels
@@ -313 +572 @@
 
          j(:,:) = 0.
+         if (photoheat) e(:,:) = 0.
 
       end if
 
-! photodissociation rates numbering in the lookup table
-
-! jmars.20140930
-
-
-      j_o2_o         =  1      ! o2 + hv     -> o + o
-      j_o2_o1d       =  2      ! o2 + hv     -> o + o(1d)
-      j_co2_o        =  3      ! co2 + hv    -> co + o
-      j_co2_o1d      =  4      ! co2 + hv    -> co + o(1d)
-      j_o3_o1d       =  5      ! o3 + hv     -> o2 + o(1d)
-      j_o3_o         =  6      ! o3 + hv     -> o2 + o
-      j_h2o          =  7      ! h2o + hv    -> h + oh
-      j_h2o2         =  8      ! h2o2 + hv   -> oh + oh
-      j_ho2          =  9      ! ho2 + hv    -> oh + o
-      j_no           =  10     ! no + hv     -> n + o
-      j_no2          =  11     ! no2 + hv    -> no + o
-      j_hno3         =  12     ! hno3 + hv   -> no2 + oh
-      j_hno4         =  13     ! hno4 + hv   -> no2 + ho2
-
-! jmars.20111014
-
-!     j_o2_o         =  1      ! o2 + hv     -> o + o
-!     j_o2_o1d       =  2      ! o2 + hv     -> o + o(1d)
-!     j_co2_o        =  3      ! co2 + hv    -> co + o
-!     j_co2_o1d      =  4      ! co2 + hv    -> co + o(1d)
-!     j_o3_o1d       =  5      ! o3 + hv     -> o2 + o(1d)
-!     j_o3_o         =  6      ! o3 + hv     -> o2 + o
-!     j_h2o          =  7      ! h2o + hv    -> h + oh
-!     j_hdo          =  8      ! hdo + hv    -> d + oh
-!     j_h2o2         =  9      ! h2o2 + hv   -> oh + oh
-!     j_ho2          =  10     ! ho2 + hv    -> oh + o
-!     j_no2          =  11     ! no2 + hv    -> no + o
-!     j_ch4_ch3_h    =  12     ! ch4 + hv    -> ch3 + h
-!     j_ch4_1ch2_h2  =  13     ! ch4 + hv    -> 1ch2 + h2
-!     j_ch4_3ch2_h_h =  14     ! ch4 + hv    -> 3ch2 + h + h
-!     j_ch4_ch_h2_h  =  15     ! ch4 + hv    -> ch + h2 + h
-!     j_ch3o2h       =  16     ! ch3o2h + hv -> ch3o + oh
-!     j_ch2o_hco     =  17     ! ch2o + hv   -> h + hco
-!     j_ch2o_co      =  18     ! ch2o + hv   -> h2 + co
-!     j_ch3oh        =  19     ! ch3oh + hv  -> ch3o + h
-!     j_c2h6         =  20     ! c2h6 + hv   -> products
-!     j_hcl          =  21     ! hcl + hv    -> h + cl
-!     j_hocl         =  22     ! hocl + hv   -> oh + cl
-!     j_clo          =  23     ! clo + hv    -> o + cl
-!     j_so2          =  24     ! so2 + hv    -> so + o
-!     j_so           =  25     ! so + hv     -> s + o
-!     j_h2s          =  26     ! h2s + hv    -> hs + s
-!     j_so3          =  27     ! so2 + hv    -> so2 + o
-!     j_hno3         =  28     ! hno3 + hv   -> oh + no2
-!     j_hno4         =  29     ! hno4 + hv   -> ho2 + no2
-
 ! fill v_phot array
 
       v_phot(:,:) = 0.
-
-      do l = 1,lswitch-1
-         v_phot(l, 1) = j(l,j_o2_o)
-         v_phot(l, 2) = j(l,j_o2_o1d)
-         v_phot(l, 3) = j(l,j_co2_o)
-         v_phot(l, 4) = j(l,j_co2_o1d)
-         v_phot(l, 5) = j(l,j_o3_o1d)
-         v_phot(l, 6) = j(l,j_o3_o)
-         v_phot(l, 7) = j(l,j_h2o)
-         v_phot(l, 8) = j(l,j_h2o2)
-         v_phot(l, 9) = j(l,j_ho2)
-         v_phot(l,10) = j(l,j_no)
-         v_phot(l,11) = j(l,j_no2)
+      e_phot(:,:) = 0.
+
+! Order of photolysis reaction has to be the same in monoreact file and the phototable file
+      ij      = 0
+      nb_phot = 0
+      do while(nb_phot<nb_phot_hv_max)
+        ij      = ij + 1
+        nb_phot = nb_phot + 1
+        do l = 1,lswitch-1
+          v_phot(l,nb_phot) = j(l,ij)
+          if (photoheat) e_phot(l,nb_phot) = e(l,ij)
+        end do
+        if (three_prod(ij)) then
+          nb_phot = nb_phot + 1
+          do l = 1,lswitch-1
+            v_phot(l,nb_phot) = j(l,ij)
+            if (photoheat) e_phot(l,nb_phot) = e(l,ij)
+          end do
+        end if
       end do
 

trunk/LMDZ.GENERIC/libf/aeronostd/read_phototable.F90

@@ -13 +13 @@
 !
 !   Author:   Franck Lefevre
+!   update 06/03/2021 dimension set in table + CH4 dimension + photoheat (Yassin Jaziri)
 !
 !   Arguments:
 !   ----------
 !
-!   The output variable is jphot and is put in common chimiedata.
+!   The output variable is jphot/ephot and is put in common chimiedata.
 !
 !***********************************************************************
@@ -24 +25 @@
       use ioipsl_getin_p_mod, only: getin_p
       use datafile_mod
+      use callkeys_mod
+      use chimiedata_h
       implicit none
 
-#include "chimiedata.h"
+!#include "chimiedata.h"
 
 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
@@ -32 +35 @@
 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
-      integer :: fic, ij, iozo, isza, itemp, iz, itau, ierr
-      real    :: xsza
+      integer :: fic, ij, iozo, isza, itemp, iz, itau, ich4, ierr
+      real    :: xsza, dummy
 
       character(len = 128) :: phototable ! photolysis table file name
+      character(len = 128) :: ephototable ! energie photolysis table file name
 
 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
@@ -73 +77 @@
       print*, 'read photolysis lookup table ',trim(phototable)
 
+      read(fic,*) nz, nsza, nozo, dummy, ntemp, ntau, nch4, nd
+
+      allocate(jphot(ntemp,nsza,nz,nozo,ntau,nch4,nd))
+      allocate(ephot(ntemp,nsza,nz,nozo,ntau,nch4,nd))
+      allocate(colairtab(nz))
+      allocate(szatab(nsza))
+      allocate(table_ozo(nozo))
+      allocate(tautab(ntau))
+      allocate(table_ch4(nch4))
+
       do itau = 1,ntau
          do itemp = 1,ntemp
             do iozo = 1,nozo
-               do isza = 1,nsza
-                  do iz = nz,1,-1
-                     read(fic,*) colairtab(iz), xsza, table_ozo(iozo)
-                     read(fic,'(7e11.4)') (jphot(itemp,isza,iz,iozo,itau,ij), ij= 1,nd)
-                     do ij = 1,nd
-                        if (jphot(itemp,isza,iz,iozo,itau,ij) == 1.e-30) then
-                           jphot(itemp,isza,iz,iozo,itau,ij) = 0.
-                        end if
+               do ich4 =1,nch4
+                  do isza = 1,nsza
+                     do iz = nz,1,-1
+                        read(fic,*) colairtab(iz), szatab(isza),table_ozo(iozo), dummy, dummy, dummy, table_ch4(ich4)
+                        read(fic,'(7e11.4)') (jphot(itemp,isza,iz,iozo,itau,ich4,ij), ij= 1,nd)
+                        do ij = 1,nd
+                           if (jphot(itemp,isza,iz,iozo,itau,ich4,ij) == 1.e-30) then
+                              jphot(itemp,isza,iz,iozo,itau,ich4,ij) = 0.
+                           end if
+                        end do
                      end do
                   end do
@@ -94 +110 @@
       close(fic)
 
+      if (photoheat) then 
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+! set energie photolysis table input file name
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+
+
+! look for a " ephototable= ..." option in def files
+     write(*,*) "Directory where external input files are:"
+      ephototable = "zdtearth_N+_CO2_0.0004_O2_0.21_G11" ! default
+     call getin_p("ephototable",ephototable) ! default path
+     write(*,*) " ephototable = ",trim(ephototable)
+
+
+      fic = 81
+
+      open(fic, form = 'formatted', status = 'old',                &
+           file =trim(datadir)//"/"//trim(ephototable),iostat=ierr)
+
+      if (ierr /= 0) THEN
+        write(*,*)'Error : cannot open energie photolysis lookup table ', trim(ephototable)
+        write(*,*)'It should be in :',trim(datadir),'/'
+        write(*,*)'1) You can change this directory in callphys.def'
+        write(*,*)'   with:'
+        write(*,*)'   datadir=/path/to/the/directory'
+        write(*,*)'2) You can change the input ephototable file name in'
+        write(*,*)'   callphys.def with:'
+        write(*,*)'   ephototable=filename'
+        stop
+      end if
+
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+! read energie photolys table
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      print*, 'read energie photolysis lookup table ',trim(ephototable)
+
+      do itau = 1,ntau
+         do itemp = 1,ntemp
+            do iozo = 1,nozo
+               do ich4 =1,nch4
+                  do isza = 1,nsza
+                     do iz = nz,1,-1
+                        read(fic,*) colairtab(iz), xsza,table_ozo(iozo), dummy, dummy, dummy, table_ch4(ich4)
+                        read(fic,'(7e11.4)') (ephot(itemp,isza,iz,iozo,itau,ich4,ij), ij= 1,nd)
+                        do ij = 1,nd
+                           if (ephot(itemp,isza,iz,iozo,itau,ich4,ij) == 1.e-30) then
+                              ephot(itemp,isza,iz,iozo,itau,ich4,ij) = 0.
+                           end if
+                        end do
+                     end do
+                  end do
+               end do
+            end do
+         end do
+      end do
+
+      print*, 'lookup etable...ok'
+      close(fic)
+      end if
+
       return
       end

trunk/LMDZ.GENERIC/libf/aeronostd/types_asis.F90

@@ -7 +7 @@
 
 type z3spec
-        real(kind=jprb)    :: z1
-        integer(kind=jpim) :: z2
-        real(kind=jprb)    :: z3
-        integer(kind=jpim) :: z4
-        real(kind=jprb)    :: z5
-        integer(kind=jpim) :: z6
+    real(kind=jprb)    :: z1
+    integer(kind=jpim) :: z2
+    real(kind=jprb)    :: z3
+    integer(kind=jpim) :: z4
+    real(kind=jprb)    :: z5
+    integer(kind=jpim) :: z6
 end type z3spec
 type z4spec
-        real(kind=jprb)    :: z1
-        integer(kind=jpim) :: z2
-        real(kind=jprb)    :: z3
-        integer(kind=jpim) :: z4
-        real(kind=jprb)    :: z5
-        integer(kind=jpim) :: z6
-        real(kind=jprb)    :: z7
-        integer(kind=jpim) :: z8
+    real(kind=jprb)    :: z1
+    integer(kind=jpim) :: z2
+    real(kind=jprb)    :: z3
+    integer(kind=jpim) :: z4
+    real(kind=jprb)    :: z5
+    integer(kind=jpim) :: z6
+    real(kind=jprb)    :: z7
+    integer(kind=jpim) :: z8
 end type z4spec
 
@@ -31 +31 @@
 type(z4spec), allocatable, save :: indice_4(:)
 
+! dimension of indexes variables
+
+integer, save :: nb_phot_max       ! dimension of phot reaction, including photolysis and quenching reaction
+integer, save :: nb_reaction_3_max
+integer, save :: nb_reaction_4_max
+integer, save :: nb_phot_hv_max    ! dimension of photolysis, including three product photolysis
+
+! photolysis reaction rate and label
+
+real (kind = 8),      allocatable, save :: v_phot_3d(:,:,:)
+character(len = 20),  allocatable, save :: jlabel(:,:)        ! photolysis label and species name
+character(len = 300), allocatable, save :: jparamline(:)      ! line of jonline parameters
+
+! pfunc type
+
+type rtype1
+    real(kind=jprb)    :: a
+    real(kind=jprb)    :: b
+    real(kind=jprb)    :: c
+    real(kind=jprb)    :: t0
+    real(kind=jprb)    :: d
+end type rtype1
+
+type rtype2
+    real(kind=jprb)    :: k0
+    real(kind=jprb)    :: n
+    real(kind=jprb)    :: a
+    real(kind=jprb)    :: kinf
+    real(kind=jprb)    :: m
+    real(kind=jprb)    :: b
+    real(kind=jprb)    :: t0
+    real(kind=jprb)    :: fc
+    real(kind=jprb)    :: g
+    real(kind=jprb)    :: h
+    real(kind=jprb)    :: dup
+    real(kind=jprb)    :: ddown
+end type rtype2
+
+type rtype3
+    real(kind=jprb)    :: k0
+    real(kind=jprb)    :: n
+    real(kind=jprb)    :: a
+    real(kind=jprb)    :: kinf
+    real(kind=jprb)    :: m
+    real(kind=jprb)    :: b
+    real(kind=jprb)    :: t0
+    real(kind=jprb)    :: atroe
+    real(kind=jprb)    :: btroe
+    real(kind=jprb)    :: ctroe
+    real(kind=jprb)    :: dtroe
+    real(kind=jprb)    :: dup
+    real(kind=jprb)    :: ddown
+end type rtype3
+
+! pfunc parameters for the reaction rates
+
+type(rtype1), allocatable, save :: pfunc1_param(:)
+type(rtype2), allocatable, save :: pfunc2_param(:)
+type(rtype3), allocatable, save :: pfunc3_param(:)
+
+! dimension of pfunc type variables
+
+integer, save :: nb_hv_max
+integer, save :: nb_pfunc1_max
+integer, save :: nb_pfunc2_max
+integer, save :: nb_pfunc3_max
+
 end module types_asis