Index: trunk/LMDZ.GENERIC/libf/aeronostd/photochemistry_asis.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/aeronostd/photochemistry_asis.F90	(revision 2542)
+++ trunk/LMDZ.GENERIC/libf/aeronostd/photochemistry_asis.F90	(revision 2553)
@@ -95,7 +95,4 @@
 
 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
@@ -115,5 +112,5 @@
 if (firstcall) then
    print*,'photochemistry: initialize indexes'
-   call indice(nreact,rtype,third_body,three_prod)
+   call indice(nreact)
    allocate(v_phot(nlayer,nb_phot_max))
    allocate(v_3(nlayer,nb_reaction_3_max))
@@ -135,7 +132,7 @@
        print*,'set jonline=.false. in callphys.def'
        print*,'or set photolysis reactions in monoreact file'
-       stop
+       call abort
      end if
-     call init_photolysis(nlayer,nreact,three_prod)     ! on-line photolysis
+     call init_photolysis(nlayer,nreact)     ! on-line photolysis
      allocate(albedo_snow_chim(nw))
      allocate(albedo_co2_ice_chim(nw))
@@ -169,5 +166,5 @@
    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)
+                             tau, sza, dist_sol, v_phot, e_phot, ig, ngrid, nreact)
       if (ngrid.eq.1) then
         do iphot = 1,nb_phot_hv_max
@@ -198,5 +195,5 @@
 !===================================================================
 
-call reactionrates(nlayer, nreact, rtype, third_body, three_prod, c, lswitch, dens, &
+call reactionrates(nlayer, nreact, c, lswitch, dens, &
                    press, temp, zmmean, sza, surfdust1d, surfice1d, v_phot, v_3, v_4)
 
@@ -296,5 +293,5 @@
 #else
    write(*,*) "photochemistry_asis error, missing LAPACK routine dgesv"
-   stop
+   call abort
 #endif
 
@@ -412,5 +409,5 @@
 #else
    write(*,*) "photochemistry_asis error, missing LAPACK routine dgesv"
-   stop
+   call abort
 #endif
 
@@ -457,5 +454,5 @@
 !======================================================================
 
- subroutine reactionrates(nlayer, nreact, rtype, third_body, three_prod, c, &
+ subroutine reactionrates(nlayer, nreact, c,                    &
                           lswitch, dens, press, t, zmmean, sza, &
                           surfdust1d, surfice1d,                &
@@ -498,8 +495,4 @@
 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.
-
 !----------------------------------------------------------------------
 !     output
@@ -537,145 +530,71 @@
 !----------------------------------------------------------------------
 
-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
+do ireact = 1,nreact
+  if (reactions(ireact)%rtype==0) then
+    nb_hv = nb_hv + 1
+    cycle
+  end if
+
+  ! get right coefficient rate function
+  if (reactions(ireact)%rtype==1) then
+    nb_pfunc1 = nb_pfunc1 + 1
+    if (reactions(ireact)%third_body==0) then            !! third_body check
+      ratec = pfunc1(nlayer,t,dens,pfunc1_param(nb_pfunc1))
     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
+      ratec = pfunc1(nlayer,t,c(:,reactions(ireact)%third_body),pfunc1_param(nb_pfunc1))
     end if
-
-    ! fill the right type index
+  else if (reactions(ireact)%rtype==2) then
+    nb_pfunc2 = nb_pfunc2 + 1
+    if (reactions(ireact)%third_body==0) then            !! third_body check
+      ratec = pfunc2(nlayer,t,dens,pfunc2_param(nb_pfunc2))
+    else
+      ratec = pfunc2(nlayer,t,c(:,reactions(ireact)%third_body),pfunc2_param(nb_pfunc2))
+    end if
+  else if (reactions(ireact)%rtype==3) then
+    nb_pfunc3 = nb_pfunc3 + 1
+    if (reactions(ireact)%third_body==0) then            !! third_body check
+      ratec = pfunc3(nlayer,t,dens,pfunc3_param(nb_pfunc3))
+    else
+      ratec = pfunc3(nlayer,t,c(:,reactions(ireact)%third_body),pfunc3_param(nb_pfunc3))
+    end if
+  else
+    print*, 'Error in reactionrates: wrong index coefficient rate vphot'
+    print*, 'rtype(ireact) = ',reactions(ireact)%rtype
+    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'
+    call abort
+  end if
+
+  ! fill the right type index
+  if (reactions(ireact)%vtype=="v4") then
+    nb_reaction_4 = nb_reaction_4 + 1
+    v_4(:,nb_reaction_4) = ratec(:)
+    if (reactions(ireact)%three_prod) then             ! three_prod check
+      nb_reaction_4 = nb_reaction_4 + 1
+      v_4(:,nb_reaction_4) = ratec(:)
+    end if
+  else if (reactions(ireact)%vtype=="v3") then
+    nb_reaction_3 = nb_reaction_3 + 1
+    v_3(:,nb_reaction_3) = ratec(:)
+    if (reactions(ireact)%three_prod) then             ! three_prod check
+      nb_reaction_3 = nb_reaction_3 + 1
+      v_3(:,nb_reaction_3) = ratec(:)
+    end if
+  else if (reactions(ireact)%vtype=="vphot") then
     nb_phot = nb_phot + 1
     v_phot(:,nb_phot) = ratec(:)
-    if (three_prod(ireact)) then             ! three_prod check
+    if (reactions(ireact)%three_prod) then             ! three_prod check
       nb_phot = nb_phot + 1
       v_phot(:,nb_phot) = ratec(:)
     end if
   else
-    nb_hv = nb_hv + 1
+    print*, 'Error in reactionrates: wrong type coefficient rate'
+    print*, 'vtype(ireact) = ',reactions(ireact)%vtype
+    print*, 'It should be v4, v3 or vphot '
+    print*, 'Please check the reaction ',ireact
+    if (ireact>nreact) print*, 'Please check the reaction count into the code beacause ireact > nreact is not possible'
+    call abort
   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
@@ -688,6 +607,4 @@
 
 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.  &
@@ -704,5 +621,5 @@
     print*, 'n4_hard_coding    = ', n4_hard_coding
     print*, 'n3_hard_coding    = ', n3_hard_coding
-    stop
+    call abort
  end if
  if ((nb_hv /= nb_hv_max)         .or.  &
@@ -719,5 +636,5 @@
     print*, 'nb_pfunc2_max = ', nb_pfunc2_max
     print*, 'nb_pfunc3_max = ', nb_pfunc3_max
-    stop
+    call abort
  end if
 firstcall = .false.
@@ -843,5 +760,5 @@
 !================================================================
 
- subroutine indice(nreact,rtype,third_body,three_prod)
+ subroutine indice(nreact)
 
 !================================================================
@@ -869,7 +786,4 @@
 
 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
@@ -880,11 +794,6 @@
 
 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                           :: nlines          ! number of lines in reactions file
 integer                           :: pnreact         ! number of reaction in photochemical reactions file
 integer                           :: bnreact         ! number of reaction in bimolecular reactions file
@@ -909,7 +818,5 @@
 
 ! Get number of reactions
-pnlines        = 0
-bnlines        = 0
-qnlines        = 0
+nlines         = 0
 nreact         = 0
 pnreact        = 0
@@ -917,23 +824,10 @@
 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
+write(*,*) "Read reaction file"
+reactfile = "reactfile" ! default
+call getin_p("reactfile",reactfile) ! default path
+write(*,*) " reactfile = ",trim(reactfile)
+
+call count_react(reactfile,nlines,nreact,nb_phot,nb_reaction_4,nb_reaction_3,nb_hv,nb_pfunc1,nb_pfunc2,nb_pfunc3)
 
 !!!!!!!!!!! Hard coding reaction !!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -964,9 +858,7 @@
 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(reactions(nreact))
 allocate(jlabel(nb_phot_hv_max,2))
-allocate(jparamline(nb_phot_hv_max))
+allocate(jparamline(nb_hv_max))
 allocate(pfunc1_param(nb_pfunc1_max))
 allocate(pfunc2_param(nb_pfunc2_max))
@@ -979,26 +871,19 @@
 specheckr(:)    = 0
 specheckp(:)    = 0
-rtype(:)        = 0
-third_body(:)   = 0
+reactions(:)    = reaction('',-1,0,.false.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,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.
+nb_phot         = 0
+nb_reaction_4   = 0
+nb_reaction_3   = 0
 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)
+call get_react(reactfile,nlines,nreact,specheck,specheckr,specheckp,nbq,nb_phot,nb_reaction_4,nb_reaction_3)
 
 ! rewrite jlabel corretly for output
 do ireact=1,nb_phot_hv_max
-  if (three_prod(ireact)) then
+  if (reactions(ireact)%three_prod) 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)
@@ -1031,5 +916,5 @@
     else
       print*, 'Error in indice: bad value in specheckr or specheckp'
-      stop
+      call abort
     end if
   end if
@@ -1067,5 +952,5 @@
    print*, 'nb_phot_hv_max    = ', nb_phot_hv_max
    print*, 'nb_hv_max         = ', nb_hv_max
-   stop
+   call abort
 end if  
 
@@ -1195,4 +1080,11 @@
         do
             n = n + 1
+            if (n>nmax) then
+              print*,'Error in split_str: to much words'
+              print*,'limit = ',nmax
+              print*,'change it, if you want, with increasing nmax'
+              print*,'line is:',line
+              call abort
+            end if
             read(line,*,iostat=ios) buf(1:n)  ! use list-directed input
             if (ios==0) then
@@ -1202,10 +1094,4 @@
                 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
@@ -1213,27 +1099,92 @@
 !*****************************************************************
 
-    subroutine count_react(reactfile,nlines,nreact,nrtype,nb_hv,nb_pfunc1,nb_pfunc2,nb_pfunc3)
+    subroutine find_vtype(reactline,vtype)
 
 !*****************************************************************
 
+        use tracer_h,     only: nesp
+        use chimiedata_h, only: indexchim
+
+        implicit none
+
+        character(len = 50), intent(in)    :: reactline ! all reactants of one reaction
+        character(*),        intent(inout) :: vtype     ! "v3", "v4" or "vphot"
+
+        integer             :: nwr                      ! number of reactant 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             :: stochio(nesp)            ! stoichiometry of reaction
+        integer             :: iword
+
+        ! init
+        stochio(:) = 0
+
+        ! split reactant
+        call split_str(reactline,words,nwr,nmax)
+
+        ! set stochio
+        do iword = 1,nwr
+            if (trim(words(iword))=="M") exit
+            if (trim(words(iword))/="hv"         &
+                .and. trim(words(iword))/="dummy") stochio(indexchim(words(iword))) = stochio(indexchim(words(iword))) + 1
+        end do
+
+        ! set vtype
+        if (sum(stochio)==1) then
+            vtype = "vphot"
+        else if (sum(stochio)==2) then
+            if (any(stochio==2)) then
+                vtype = "v3"
+            else
+                vtype = "v4"
+            end if
+        else if (sum(stochio)==3) then
+            if (any(stochio==2)) then
+                vtype = "v4"
+            else
+                print*,'Error in photochemistry_asis (find_vtype):'
+                print*,'3 different reactants not implemented'
+                call abort
+            end if
+        else
+            print*,'Error in photochemistry_asis (find_vtype):'
+            print*,'more than 2 different reactants not implemented'
+            call abort
+        end if
+
+    end subroutine find_vtype
+
+!*****************************************************************
+
+    subroutine count_react(reactfile,nlines,nreact,nb_phot,nb_reaction_4,nb_reaction_3,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
+        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)  :: nb_phot       ! dimension of "vphot" reaction
+        integer,      intent(inout)  :: nb_reaction_4 ! dimension of "v4" reaction
+        integer,      intent(inout)  :: nb_reaction_3 ! dimension of "v3" reaction
+        integer,      intent(inout)  :: nb_hv         ! number of typerate 0 reaction
+        integer,      intent(inout)  :: nb_pfunc1     ! number of typerate 1 reaction
+        integer,      intent(inout)  :: nb_pfunc2     ! number of typerate 2 reaction
+        integer,      intent(inout)  :: nb_pfunc3     ! number of typerate 3 reaction
 
         ! 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
+        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)
+        character(5)        :: vtype                  ! "v3", "v4" or "vphot"
+        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
+        nreact = 0
+        vtype  = ''
 
         open(402, form = 'formatted', status = 'old',                &
@@ -1247,5 +1198,5 @@
            write(*,*)'   monoreact=filename or bimolreact=filename or quadrareact=filename'
            write(*,*)'   depending on what chemical reaction type it is'
-           stop
+           call abort
         end if
     
@@ -1253,14 +1204,37 @@
           read(402,'(A,A,I2)',iostat=ierr) reactline, prodline, typerate
           if (ierr<0) exit
+          ! count the number of lines
           nlines = nlines + 1
           if (reactline(1:1)/='!' .and. reactline(1:1)/='') then
             ! count the number of reaction
             nreact = nreact + 1
+
+            call find_vtype(reactline,vtype)
             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 dimension of each rate reaction coefficient type
+            if (trim(vtype)=="vphot") then
+                nb_phot = nb_phot + 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))) nb_phot = nb_phot + 1
+            else if (trim(vtype)=="v4") then
+                nb_reaction_4 = nb_reaction_4 + 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))) nb_reaction_4 = nb_reaction_4 + 1
+            else if (trim(vtype)=="v3") then
+                nb_reaction_3 = nb_reaction_3 + 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))) nb_reaction_3 = nb_reaction_3 + 1
+            else
+                print*,'Error in photochemistry_asis (count_react):'
+                print*,'vtype not found'
+                call abort
+            end if
 
             ! count the number of each rate reaction coefficient type
@@ -1282,5 +1256,5 @@
               print*, 'It should be 0 for photolysis reations and 1 or 2 for the others'
               print*, 'And not : ', typerate
-              stop
+              call abort
             end if
 
@@ -1295,7 +1269,6 @@
 !*****************************************************************
 
-    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)
+    subroutine get_react(reactfile,nlines,nreact,specheck,specheckr,specheckp, &
+                         nbq,nb_phot,nb_reaction_4,nb_reaction_3)
 
 !*****************************************************************
@@ -1310,16 +1283,11 @@
         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(out)    :: nb_phot            ! number of "vphot" calculation reactions
+        integer,      intent(out)    :: nb_reaction_4      ! number of "vphot" calculation reactions
+        integer,      intent(out)    :: nb_reaction_3      ! number of "vphot" 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
@@ -1333,17 +1301,16 @@
         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
+        integer              :: stochio(nesp)              ! stoichiometry of reaction
+        integer              :: ierr, ilines, ireact, i_dummy, iw, ispe, iloc
+        integer              :: nb_phot_hv, nb_hv, nb_pfunc1, nb_pfunc2, nb_pfunc3
+
+        i_dummy       = 1
+        ireact        = 0
+        nb_phot       = 0
+        nb_phot_hv    = 0
+        nb_hv         = 0
+        nb_pfunc1     = 0
+        nb_pfunc2     = 0
+        nb_pfunc3     = 0
 
         open(402, form = 'formatted', status = 'old',                &
@@ -1357,5 +1324,5 @@
            write(*,*)'   monoreact=filename or bimolreact=filename or quadrareact=filename'
            write(*,*)'   depending on what chemical reaction type it is'
-           stop
+           call abort
         end if
 
@@ -1368,195 +1335,194 @@
           if (reactline(1:1)/='!' .and. reactline(1:1)/='') then
 
-            ! increment number of reactions and init
+            ! increment number of reactions
             ireact = ireact + 1
-        !!!!!!!!!!!!!!!!!!!!!!!!!!! for fill indice part
-            nrtype = nrtype + 1
-            nindice(:) = 0.0
-            iindice(:) = i_dummy
-        !!!!!!!!!!!!!!!!!!!!!!!!!!! end
-            ! get indice, rate type and parameters
+
+            ! fill reaction vtype
+            call find_vtype(reactline,reactions(ireact)%vtype)
+            if (trim(reactions(ireact)%vtype)=="v4") then
+                nb_reaction_4 = nb_reaction_4 + 1
+            else if (trim(reactions(ireact)%vtype)=="v3") then
+                nb_reaction_3 = nb_reaction_3 + 1
+            else if (trim(reactions(ireact)%vtype)=="vphot") then
+                nb_phot = nb_phot + 1
+            else
+                print*,'Error in photochemistry_asis (get_react):'
+                print*,'vtype not found'
+                call abort
+            end if
+
+            ! fill reaction rate type and parameters
             if (trim(paramline)=='') then
-              print*, 'Error in reactfile: where are the parameters - line',ilines
-              stop
+                print*, 'Error in reactfile: where are the parameters - line',ilines
+                call abort
             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)
+                read(paramline,*) reactions(ireact)%rtype
+                if (reactions(ireact)%rtype==1) then
+                  nb_pfunc1 = nb_pfunc1 + 1
+                  read(paramline,*) reactions(ireact)%rtype, pfunc1_param(nb_pfunc1)
+                else if (reactions(ireact)%rtype==0) then
+                  nb_hv = nb_hv + 1
+                  nb_phot_hv = nb_phot_hv + 1
+                  if (jonline) then
+                    read(paramline,'(I5,A,A)') reactions(ireact)%rtype, jlabel(nb_hv,1), jparamline(nb_hv)
+                  else
+                    read(paramline,*) reactions(ireact)%rtype, jlabel(nb_hv,1)
+                  end if
+                else if (reactions(ireact)%rtype==2) then
+                  nb_pfunc2 = nb_pfunc2 + 1
+                  read(paramline,*) reactions(ireact)%rtype, pfunc2_param(nb_pfunc2)
+                else if (reactions(ireact)%rtype==3) then
+                  nb_pfunc3 = nb_pfunc3 + 1
+                  read(paramline,*) reactions(ireact)%rtype, pfunc3_param(nb_pfunc3)
+                end if
+            end if
+
+            ! WARNING: the implementation is limited to 3 different reactants (for now only 2) and 3 different products
+            nindice(:) = 0.0     ! 3 first indice for reactants, 3 following for products
+            iindice(:) = i_dummy ! 3 first indice for reactants, 3 following for products
+
+            !-----------------!
+            !--- reactants ---!
+            !-----------------!
+
+            ! init for reactant
+            stochio(:) = 0
+            ! split reactant
+            call split_str(reactline,words,nwr,nmax)
+            ! set species that are photolysed
+            if (reactions(ireact)%rtype==0) jlabel(nb_hv,2) = words(1)
+            ! find reactant stochio
+            do iw = 1,nwr
+                ! fill third body index
+                if (trim(words(iw))=="M") then
+                    if (iw==nwr) then
+                        exit
+                    else if (iw==nwr-1) then
+                        reactions(ireact)%third_body = 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
+                        call abort
+                    end if
+                end if
+                ! count stochio
+                if (trim(words(iw))/="hv"         &
+                    .and. trim(words(iw))/="dummy") stochio(indexchim(words(iw))) = stochio(indexchim(words(iw))) + 1
+            end do
+            ! get reaction reactant stochio and indice
+            iloc = 0
+            do ispe = 1,nesp
+                if (stochio(ispe)==0) cycle
+                iloc = iloc + 1
+                nindice(iloc) = stochio(ispe)
+                iindice(iloc) = ispe
+                ! check up: to count the species used in 'reactfile'
+                if (specheck(ispe)==0) then
+                    specheckr(ispe) = 1
+                    specheck(ispe) = 1
+                    nbq = nbq + 1
+                else if (specheckr(ispe)==0) then
+                    specheckr(ispe) = 1
+                end if
+            end do
+            ! fill reaction reactant stochio and indice
+            reactions(ireact)%ir1 = nindice(1)
+            reactions(ireact)%r1  = iindice(1)
+            reactions(ireact)%ir2 = nindice(2)
+            reactions(ireact)%r2  = iindice(2)
+            reactions(ireact)%ir3 = nindice(3)
+            reactions(ireact)%r3  = iindice(3)
+
+            !----------------!
+            !--- products ---!
+            !----------------!
+
+            ! init for products
+            stochio(:) = 0
+            ! split products
+            call split_str(prodline,words,nwp,nmax)
+            ! find products stochio
+            do iw = 1,nwp
+                stochio(indexchim(words(iw))) = stochio(indexchim(words(iw))) + 1
+            end do
+            ! get reaction product stochio and indice
+            iloc = 3
+            do ispe = 1,nesp
+                if (stochio(ispe)==0) cycle
+                iloc = iloc + 1
+                nindice(iloc) = stochio(ispe)
+                iindice(iloc) = ispe
+                ! check up: to count the species used in 'reactfile'
+                if (specheck(ispe)==0) then
+                    specheckp(ispe) = 1
+                    specheck(ispe) = 1
+                    nbq = nbq + 1
+                else if (specheckp(ispe)==0) then
+                    specheckp(ispe) = 1
+                end if
+            end do
+            ! fill reaction product stochio and indice
+            reactions(ireact)%ip1 = nindice(4)
+            reactions(ireact)%p1  = iindice(4)
+            reactions(ireact)%ip2 = nindice(5)
+            reactions(ireact)%p2  = iindice(5)
+            reactions(ireact)%ip3 = nindice(6)
+            reactions(ireact)%p3  = iindice(6)
+            ! set reaction three prod to true with checking in position 6 if there is three different products
+            if (nindice(6)/=0) reactions(ireact)%three_prod = .true.
+
+
+            !-------------------------!
+            !--- fill indice array ---!
+            !-------------------------!
+
+            if (trim(reactions(ireact)%vtype)=="v4") then
+                if (nindice(6)==0) then ! reaction with 1 or 2 products
+                    ! z4spec(ir1,r1,ir2,r2,ip1,p1,ip2,p2)
+                    indice_4(nb_reaction_4) = z4spec(nindice(1), iindice(1), nindice(2), iindice(2), nindice(4), iindice(4), nindice(5), iindice(5))
+                else ! reaction with 3 different products
+                    indice_4(nb_reaction_4) = z4spec(nindice(1)/2., iindice(1), nindice(2)/2., iindice(2), nindice(4), iindice(4), nindice(5)/2., iindice(5))
+                    nb_reaction_4 = nb_reaction_4 + 1
+                    indice_4(nb_reaction_4) = z4spec(nindice(1)/2., iindice(1), nindice(2)/2., iindice(2), nindice(6), iindice(6), nindice(5)/2., iindice(5))
+                end if
+            else if (trim(reactions(ireact)%vtype)=="v3") then
+                if (nindice(6)==0) then ! reaction with 1 or 2 products
+                    ! z3spec(ir1,r1,ip1,p1,ip2,p2)
+                    indice_3(nb_reaction_3) = z3spec(nindice(1), iindice(1), nindice(4), iindice(4), nindice(5), iindice(5))
+                else ! reaction with 3 different products
+                    indice_3(nb_reaction_3) = z3spec(nindice(1)/2., iindice(1), nindice(4), iindice(4), 0.0, i_dummy)
+                    nb_reaction_3 = nb_reaction_3 + 1
+                    indice_3(nb_reaction_3) = z3spec(nindice(1)/2., iindice(1), nindice(5), iindice(5), nindice(6), iindice(6))
+                end if
+            else if (trim(reactions(ireact)%vtype)=="vphot") then
+                if (reactions(ireact)%rtype==0) then
+                    if (nindice(6)==0) then ! reaction with 1 or 2 products
+                        ! z3spec(ir1,r1,ip1,p1,ip2,p2)
+                        indice_phot(nb_phot_hv) = z3spec(nindice(1), iindice(1), nindice(4), iindice(4), nindice(5), iindice(5))
+                    else ! reaction with 3 different products
+                        indice_phot(nb_phot_hv) = z3spec(nindice(1)/2., iindice(1), nindice(4), iindice(4), 0.0, i_dummy)
+                        nb_phot = nb_phot + 1
+                        nb_phot_hv = nb_phot_hv + 1
+                        indice_phot(nb_phot_hv) = z3spec(nindice(1)/2., iindice(1), nindice(5), iindice(5), nindice(6), iindice(6))
+                    end if
                 else
-                  read(paramline,*) rtype(ireact), jlabel(nb_hv,1)
+                    if (nindice(6)==0) then ! reaction with 1 or 2 products
+                        ! z3spec(ir1,r1,ip1,p1,ip2,p2)
+                        indice_phot(nb_phot_hv_max+nb_phot-nb_phot_hv) = z3spec(nindice(1), iindice(1), nindice(4), iindice(4), nindice(5), iindice(5))
+                    else ! reaction with 3 different products
+                        indice_phot(nb_phot_hv_max+nb_phot-nb_phot_hv) = z3spec(nindice(1)/2., iindice(1), nindice(4), iindice(4), 0.0, i_dummy)
+                        nb_phot = nb_phot + 1
+                        indice_phot(nb_phot_hv_max+nb_phot-nb_phot_hv) = z3spec(nindice(1)/2., iindice(1), nindice(5), iindice(5), nindice(6), iindice(6))
+                    end if
                 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
+            else
+                print*,'Error in photochemistry_asis (get_react):'
+                print*,'vtype',reactions(ireact)%vtype,' not implemented'
+                call abort
             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
+          end if ! end if comment line
+
+        end do ! end loop on lines
 
         close(402)
Index: trunk/LMDZ.GENERIC/libf/aeronostd/photolysis_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/aeronostd/photolysis_mod.F90	(revision 2542)
+++ trunk/LMDZ.GENERIC/libf/aeronostd/photolysis_mod.F90	(revision 2553)
@@ -45,6 +45,7 @@
     contains
   
-    subroutine init_photolysis(nlayer,nreact,three_prod)
-  
+    subroutine init_photolysis(nlayer,nreact)
+  
+    use types_asis, only: reactions
     use tracer_h
     use chimiedata_h, only: indexchim
@@ -53,6 +54,5 @@
     integer, intent(in) :: nlayer              ! number of atmospheric layers
     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.
-    integer :: iphot,tdim_max,i3prod,nd,ij,iw,ilay
+    integer :: iphot,tdim_max,i3prod,ij,iw,ilay,ireact
     integer :: specheck(nesp)
   
@@ -94,11 +94,10 @@
       endif
     end do
-    nd = nb_hv_max
   
     ! Get temperature dimension and allocate
-    allocate(tdim(nd))
+    allocate(tdim(nb_hv_max))
     tdim(:)  = 1
     tdim_max = 1
-    do iphot=1,nd
+    do iphot=1,nb_hv_max
       read(jparamline(iphot),*) tdim(iphot)
       if (tdim(iphot) > tdim_max) tdim_max = tdim(iphot)
@@ -106,8 +105,8 @@
   
     ! allocation
-    allocate(qy(nw,nd))
-    allocate(xs(tdim_max,nw,nd))
+    allocate(qy(nw,nb_hv_max))
+    allocate(xs(tdim_max,nw,nb_hv_max))
     allocate(sj(nlayer,nw,nb_phot_hv_max))
-    allocate(xs_temp(tdim_max,nd))
+    allocate(xs_temp(tdim_max,nb_hv_max))
     xs(:,:,:)    = 0.
     xs_temp(:,:) = 0.
@@ -119,17 +118,24 @@
   
     i3prod = 0
-  
-    ! read and grid cross-sections (photolysis should be first reactions in monoreact file)
-    do iphot=1,nd
+    ireact = 0
+  
+    ! read and grid cross-sections
+    do iphot=1,nb_hv_max
+      ireact = ireact + 1
       call rdxsi(nw,wl,xs(:tdim(iphot),:,iphot),jparamline(iphot),xs_temp(:tdim(iphot),iphot),qy(:,iphot),tdim(iphot),jlabel(iphot+i3prod,2))
-      if (three_prod(iphot)) i3prod = i3prod + 1
+      do while(reactions(ireact)%rtype/=0)
+        ireact = ireact + 1
+      end do
+      if (reactions(ireact)%three_prod) i3prod = i3prod + 1
     end do
   
     ! init sj for no temperature dependent cross sections (tdim.eq.1)
-    iphot = 0
-    ij    = 0
+    iphot  = 0
+    ij     = 0
+    ireact = 0
     do while(iphot<nb_phot_hv_max)
-       ij = ij + 1
-       iphot = iphot + 1
+       ij     = ij     + 1
+       iphot  = iphot  + 1
+       ireact = ireact + 1
        if (tdim(ij).eq.1) then
          do iw = 1,nw-1
@@ -139,5 +145,8 @@
          end do
        endif
-       if (three_prod(ij)) then
+       do while(reactions(ireact)%rtype/=0)
+         ireact = ireact + 1
+       end do
+       if (reactions(ireact)%three_prod) then
          iphot = iphot + 1
          if (tdim(ij).eq.1) then
@@ -882,5 +891,5 @@
         print*, 'ERROR: temperature cross section file'
         print*, '       has to be sorted from low to high values'
-        print*, '       Check monoreact file'
+        print*, '       Check reactfile file'
         stop
       end if
@@ -903,5 +912,5 @@
          write(*,*)'2) You can check if the file is in datadir/cross_sections/'
          write(*,*)'3) You can change the input cross_sections file name in'
-         write(*,*)'   chemestry/monoreact file for the specie:',trim(species)
+         write(*,*)'   chemestry/reactfile file for the specie:',trim(species)
          stop
       end if
@@ -960,5 +969,5 @@
        write(*,*)'2) You can check if the file is in datadir/cross_sections/'
        write(*,*)'3) You can change the input photodissociation yield file name in'
-       write(*,*)'   chemestry/monoreact file for the specie:',trim(species)
+       write(*,*)'   chemestry/reactfile file for the specie:',trim(species)
        stop
     end if
Index: trunk/LMDZ.GENERIC/libf/aeronostd/photolysis_online.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/aeronostd/photolysis_online.F	(revision 2542)
+++ trunk/LMDZ.GENERIC/libf/aeronostd/photolysis_online.F	(revision 2553)
@@ -4,5 +4,5 @@
      $           zmmean, rm,
      $           tau, sza, dist_sol, v_phot, e_phot, ig, ngrid,
-     $           nreact, three_prod)
+     $           nreact)
 
 !***********************************************************************
@@ -21,5 +21,6 @@
       use tracer_h
       use chimiedata_h, only: indexchim
-      use types_asis,   only: nb_phot_hv_max, nb_phot_max, jlabel
+      use types_asis,   only: nb_phot_hv_max, nb_phot_max, jlabel,
+     $                        reactions
 
       implicit none
@@ -38,5 +39,4 @@
       integer, intent(in) :: ig                                     ! grid point index
       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.
 
 !     output
@@ -71,5 +71,5 @@
       real,    dimension(0:nlayer,nlayer) :: dsdh
 
-      integer :: i, ilay, iw, ialt, iphot, ispe, ij, igas
+      integer :: i, ilay, iw, ialt, iphot, ispe, ij, igas, ireact
       real    :: deltaj
 
@@ -85,11 +85,13 @@
 !==== set temperature-dependent cross-sections and optical depths
 
-      iphot = 0
-      ij    = 0
-      igas  = 0
+      iphot  = 0
+      ij     = 0
+      igas   = 0
+      ireact = 0
 
       do while(iphot<nb_phot_hv_max)
-        ij = ij + 1
-        iphot = iphot + 1
+        ij     = ij     + 1
+        iphot  = iphot  + 1
+        ireact = ireact + 1
 
         if (tdim(ij).eq.1) then
@@ -106,5 +108,8 @@
             igas = igas + 1
           end if
-          if (three_prod(ij)) iphot = iphot + 1
+          do while(reactions(ireact)%rtype/=0)
+            ireact = ireact + 1
+          end do
+          if (reactions(ireact)%three_prod) iphot = iphot + 1
         else
           call setsj(nb_phot_hv_max,nlayer,nw,temp,tdim(ij),
@@ -131,5 +136,8 @@
           end do
 
-          if (three_prod(ij)) then
+          do while(reactions(ireact)%rtype/=0)
+            ireact = ireact + 1
+          end do
+          if (reactions(ireact)%three_prod) then
             iphot = iphot + 1
             do iw = 1,nw-1
Index: trunk/LMDZ.GENERIC/libf/aeronostd/types_asis.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/aeronostd/types_asis.F90	(revision 2542)
+++ trunk/LMDZ.GENERIC/libf/aeronostd/types_asis.F90	(revision 2553)
@@ -85,4 +85,25 @@
 end type rtype3
 
+! reaction type
+
+type reaction
+    character(5)    :: vtype 
+    integer         :: rtype
+    integer         :: third_body
+    logical         :: three_prod
+    integer         :: ir1
+    integer         :: r1
+    integer         :: ir2
+    integer         :: r2
+    integer         :: ir3
+    integer         :: r3
+    integer         :: ip1
+    integer         :: p1
+    integer         :: ip2
+    integer         :: p2
+    integer         :: ip3
+    integer         :: p3
+end type reaction
+
 ! pfunc parameters for the reaction rates
 
@@ -90,4 +111,8 @@
 type(rtype2), allocatable, save :: pfunc2_param(:)
 type(rtype3), allocatable, save :: pfunc3_param(:)
+
+! reactions attributs
+
+type(reaction), allocatable, save :: reactions(:)
 
 ! dimension of pfunc type variables
