Changeset 2553 for trunk/LMDZ.GENERIC


Ignore:
Timestamp:
Jul 21, 2021, 10:44:46 AM (3 years ago)
Author:
yjaziri
Message:

Generic GCM:

Chemistry: one input file instead of three for the chemical network
Automatic detection of mono/bi-molecular or quadratic reactions

YJ

Location:
trunk/LMDZ.GENERIC
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/README

    r2550 r2553  
    16791679Some OpenMP fixes in routines initracer.F, nonoro_gwd_ran_mod.F90,
    16801680phys_state_var_mod.F90 and sugas_corrk.F90
     1681
     1682== 21/07/2021 == YJ
     1683Chemistry: one input file instead of three for the chemical network
     1684Automatic detection of mono/bi-molecular or quadratic reactions

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

    r2542 r2553  
    9595
    9696integer, save                          :: nreact         ! number of reactions in reactions files
    97 integer, allocatable, save             :: rtype(:)       ! reaction rate type
    98 integer, allocatable, save             :: third_body(:)  ! if the reaction have a third body: index of the third body, else zero
    99 logical, allocatable, save             :: three_prod(:)  ! if the reaction have a three defferent proucts egal .true.
    10097
    10198! matrix
    … …  
    115112if (firstcall) then
    116113   print*,'photochemistry: initialize indexes'
    117    call indice(nreact,rtype,third_body,three_prod)
     114   call indice(nreact)
    118115   allocate(v_phot(nlayer,nb_phot_max))
    119116   allocate(v_3(nlayer,nb_reaction_3_max))
    … …  
    135132       print*,'set jonline=.false. in callphys.def'
    136133       print*,'or set photolysis reactions in monoreact file'
    137        stop
     134       call abort
    138135     end if
    139      call init_photolysis(nlayer,nreact,three_prod)     ! on-line photolysis
     136     call init_photolysis(nlayer,nreact)     ! on-line photolysis
    140137     allocate(albedo_snow_chim(nw))
    141138     allocate(albedo_co2_ice_chim(nw))
    … …  
    169166   if (sza <= 95.) then
    170167      call photolysis_online(nlayer, alt, press, temp, zmmean, rm,   &
    171                              tau, sza, dist_sol, v_phot, e_phot, ig, ngrid, nreact, three_prod)
     168                             tau, sza, dist_sol, v_phot, e_phot, ig, ngrid, nreact)
    172169      if (ngrid.eq.1) then
    173170        do iphot = 1,nb_phot_hv_max
    … …  
    198195!===================================================================
    199196
    200 call reactionrates(nlayer, nreact, rtype, third_body, three_prod, c, lswitch, dens, &
     197call reactionrates(nlayer, nreact, c, lswitch, dens, &
    201198                   press, temp, zmmean, sza, surfdust1d, surfice1d, v_phot, v_3, v_4)
    202199
    … …  
    296293#else
    297294   write(*,*) "photochemistry_asis error, missing LAPACK routine dgesv"
    298    stop
     295   call abort
    299296#endif
    300297
    … …  
    412409#else
    413410   write(*,*) "photochemistry_asis error, missing LAPACK routine dgesv"
    414    stop
     411   call abort
    415412#endif
    416413
    … …  
    457454!======================================================================
    458455
    459  subroutine reactionrates(nlayer, nreact, rtype, third_body, three_prod, c, &
     456 subroutine reactionrates(nlayer, nreact, c,                    &
    460457                          lswitch, dens, press, t, zmmean, sza, &
    461458                          surfdust1d, surfice1d,                &
    … …  
    498495real (kind = 8), dimension(nlayer,nesp) :: c  ! species number density (molecule.cm-3)
    499496
    500 integer, intent(in)     :: rtype(nreact)      ! reaction rate type
    501 integer, intent(in)     :: third_body(nreact) ! if the reaction have a third body: index of the third body, else zero
    502 logical, intent(in)     :: three_prod(nreact) ! if the reaction have three different products egal .true.
    503 
    504497!----------------------------------------------------------------------
    505498!     output
    … …  
    537530!----------------------------------------------------------------------
    538531
    539 ireact = 1
    540 
    541 ! Reaction from monoreact file are implemented first
    542 do while(nb_phot<nb_phot_max-nphot_hard_coding)
    543 
    544   if (rtype(ireact)/=0) then     ! skip photolysis reactions
    545     ! get right coefficient rate function
    546     if (rtype(ireact)==1) then
    547       nb_pfunc1 = nb_pfunc1 + 1
    548       if (third_body(ireact)==0) then            !! third_body check
    549         ratec = pfunc1(nlayer,t,dens,pfunc1_param(nb_pfunc1))
    550       else
    551         ratec = pfunc1(nlayer,t,c(:,third_body(ireact)),pfunc1_param(nb_pfunc1))
    552       end if
    553     else if (rtype(ireact)==2) then
    554       nb_pfunc2 = nb_pfunc2 + 1
    555       if (third_body(ireact)==0) then            !! third_body check
    556         ratec = pfunc2(nlayer,t,dens,pfunc2_param(nb_pfunc2))
    557       else
    558         ratec = pfunc2(nlayer,t,c(:,third_body(ireact)),pfunc2_param(nb_pfunc2))
    559       end if
    560     else if (rtype(ireact)==3) then
    561       nb_pfunc3 = nb_pfunc3 + 1
    562       if (third_body(ireact)==0) then            !! third_body check
    563         ratec = pfunc3(nlayer,t,dens,pfunc3_param(nb_pfunc3))
    564       else
    565         ratec = pfunc3(nlayer,t,c(:,third_body(ireact)),pfunc3_param(nb_pfunc3))
    566       end if
     532do ireact = 1,nreact
     533  if (reactions(ireact)%rtype==0) then
     534    nb_hv = nb_hv + 1
     535    cycle
     536  end if
     537
     538  ! get right coefficient rate function
     539  if (reactions(ireact)%rtype==1) then
     540    nb_pfunc1 = nb_pfunc1 + 1
     541    if (reactions(ireact)%third_body==0) then            !! third_body check
     542      ratec = pfunc1(nlayer,t,dens,pfunc1_param(nb_pfunc1))
    567543    else
    568       print*, 'Error in reactionrates: wrong index coefficient rate vphot'
    569       print*, 'rtype(ireact) = ',rtype(ireact)
    570       print*, 'It should be 1 or 2 '
    571       print*, 'Please check the reaction ',ireact
    572       if (ireact>nreact) print*, 'Please check the reaction count into the code beacause ireact > nreact is not possible'
    573       stop
     544      ratec = pfunc1(nlayer,t,c(:,reactions(ireact)%third_body),pfunc1_param(nb_pfunc1))
    574545    end if
    575 
    576     ! fill the right type index
     546  else if (reactions(ireact)%rtype==2) then
     547    nb_pfunc2 = nb_pfunc2 + 1
     548    if (reactions(ireact)%third_body==0) then            !! third_body check
     549      ratec = pfunc2(nlayer,t,dens,pfunc2_param(nb_pfunc2))
     550    else
     551      ratec = pfunc2(nlayer,t,c(:,reactions(ireact)%third_body),pfunc2_param(nb_pfunc2))
     552    end if
     553  else if (reactions(ireact)%rtype==3) then
     554    nb_pfunc3 = nb_pfunc3 + 1
     555    if (reactions(ireact)%third_body==0) then            !! third_body check
     556      ratec = pfunc3(nlayer,t,dens,pfunc3_param(nb_pfunc3))
     557    else
     558      ratec = pfunc3(nlayer,t,c(:,reactions(ireact)%third_body),pfunc3_param(nb_pfunc3))
     559    end if
     560  else
     561    print*, 'Error in reactionrates: wrong index coefficient rate vphot'
     562    print*, 'rtype(ireact) = ',reactions(ireact)%rtype
     563    print*, 'It should be 1 or 2 '
     564    print*, 'Please check the reaction ',ireact
     565    if (ireact>nreact) print*, 'Please check the reaction count into the code beacause ireact > nreact is not possible'
     566    call abort
     567  end if
     568
     569  ! fill the right type index
     570  if (reactions(ireact)%vtype=="v4") then
     571    nb_reaction_4 = nb_reaction_4 + 1
     572    v_4(:,nb_reaction_4) = ratec(:)
     573    if (reactions(ireact)%three_prod) then             ! three_prod check
     574      nb_reaction_4 = nb_reaction_4 + 1
     575      v_4(:,nb_reaction_4) = ratec(:)
     576    end if
     577  else if (reactions(ireact)%vtype=="v3") then
     578    nb_reaction_3 = nb_reaction_3 + 1
     579    v_3(:,nb_reaction_3) = ratec(:)
     580    if (reactions(ireact)%three_prod) then             ! three_prod check
     581      nb_reaction_3 = nb_reaction_3 + 1
     582      v_3(:,nb_reaction_3) = ratec(:)
     583    end if
     584  else if (reactions(ireact)%vtype=="vphot") then
    577585    nb_phot = nb_phot + 1
    578586    v_phot(:,nb_phot) = ratec(:)
    579     if (three_prod(ireact)) then             ! three_prod check
     587    if (reactions(ireact)%three_prod) then             ! three_prod check
    580588      nb_phot = nb_phot + 1
    581589      v_phot(:,nb_phot) = ratec(:)
    582590    end if
    583591  else
    584     nb_hv = nb_hv + 1
     592    print*, 'Error in reactionrates: wrong type coefficient rate'
     593    print*, 'vtype(ireact) = ',reactions(ireact)%vtype
     594    print*, 'It should be v4, v3 or vphot '
     595    print*, 'Please check the reaction ',ireact
     596    if (ireact>nreact) print*, 'Please check the reaction count into the code beacause ireact > nreact is not possible'
     597    call abort
    585598  end if
    586 
    587   ireact = ireact + 1
    588 
    589 end do
    590 
    591 ! Reaction from bimolreact file are implemented secondly
    592 do while(nb_reaction_4<nb_reaction_4_max-n4_hard_coding)
    593 
    594   ! get right coefficient rate function
    595   if (rtype(ireact)==1) then
    596     nb_pfunc1 = nb_pfunc1 + 1
    597     if (third_body(ireact)==0) then            !! third_body check
    598       ratec = pfunc1(nlayer,t,dens,pfunc1_param(nb_pfunc1))
    599     else
    600       ratec = pfunc1(nlayer,t,c(:,third_body(ireact)),pfunc1_param(nb_pfunc1))
    601     end if
    602   else if (rtype(ireact)==2) then
    603     nb_pfunc2 = nb_pfunc2 + 1
    604     if (third_body(ireact)==0) then            !! third_body check
    605       ratec = pfunc2(nlayer,t,dens,pfunc2_param(nb_pfunc2))
    606     else
    607       ratec = pfunc2(nlayer,t,c(:,third_body(ireact)),pfunc2_param(nb_pfunc2))
    608     end if
    609   else if (rtype(ireact)==3) then
    610     nb_pfunc3 = nb_pfunc3 + 1
    611     if (third_body(ireact)==0) then            !! third_body check
    612       ratec = pfunc3(nlayer,t,dens,pfunc3_param(nb_pfunc3))
    613     else
    614       ratec = pfunc3(nlayer,t,c(:,third_body(ireact)),pfunc3_param(nb_pfunc3))
    615     end if
    616   else
    617     print*, 'Error in reactionrates: wrong index coefficient rate v4'
    618     print*, 'rtype(ireact) = ',rtype(ireact)
    619     print*, 'It should be 1 or 2 '
    620     print*, 'Please check the reaction ',ireact-nb_phot_max
    621     if (ireact>nreact) print*, 'Please check the reaction count into the code beacause ireact > nreact is not possible'
    622     stop
    623   end if
    624 
    625   ! fill the right type index
    626   nb_reaction_4 = nb_reaction_4 + 1
    627   v_4(:,nb_reaction_4) = ratec(:)
    628   if (three_prod(ireact)) then             ! three_prod check
    629     nb_reaction_4 = nb_reaction_4 + 1
    630     v_4(:,nb_reaction_4) = ratec(:)
    631   end if
    632 
    633   ireact = ireact + 1
    634 
    635 end do
    636 
    637 ! Reaction from quadrareact file are implemented thirdly
    638 do while(nb_reaction_3<nb_reaction_3_max-n3_hard_coding)
    639 
    640   ! get right coefficient rate function
    641   if (rtype(ireact)==1) then
    642     nb_pfunc1 = nb_pfunc1 + 1
    643     if (third_body(ireact)==0) then            !! third_body check
    644       ratec = pfunc1(nlayer,t,dens,pfunc1_param(nb_pfunc1))
    645     else
    646       ratec = pfunc1(nlayer,t,c(:,third_body(ireact)),pfunc1_param(nb_pfunc1))
    647     end if
    648   else if (rtype(ireact)==2) then
    649     nb_pfunc2 = nb_pfunc2 + 1
    650     if (third_body(ireact)==0) then            !! third_body check
    651       ratec = pfunc2(nlayer,t,dens,pfunc2_param(nb_pfunc2))
    652     else
    653       ratec = pfunc2(nlayer,t,c(:,third_body(ireact)),pfunc2_param(nb_pfunc2))
    654     end if
    655   else if (rtype(ireact)==3) then
    656     nb_pfunc3 = nb_pfunc3 + 1
    657     if (third_body(ireact)==0) then            !! third_body check
    658       ratec = pfunc3(nlayer,t,dens,pfunc3_param(nb_pfunc3))
    659     else
    660       ratec = pfunc3(nlayer,t,c(:,third_body(ireact)),pfunc3_param(nb_pfunc3))
    661     end if
    662   else
    663     print*, 'Error in reactionrates: wrong index coefficient rate v3'
    664     print*, 'rtype(ireact) = ',rtype(ireact)
    665     print*, 'It should be 1 or 2 '
    666     print*, 'Please check the reaction ',ireact-nb_phot_max-nb_reaction_4_max
    667     if (ireact>nreact) print*, 'Please check the reaction count into the code beacause ireact > nreact is not possible'
    668     stop
    669   end if
    670 
    671   ! fill the right type index
    672   nb_reaction_3 = nb_reaction_3 + 1
    673   v_3(:,nb_reaction_3) = ratec(:)
    674   if (three_prod(ireact)) then             ! three_prod check
    675     nb_reaction_3 = nb_reaction_3 + 1
    676     v_3(:,nb_reaction_3) = ratec(:)
    677   end if
    678 
    679   ireact = ireact + 1
    680599
    681600end do
    … …  
    688607
    689608if (firstcall) then
    690  ireact = ireact-1
    691  if (ireact /= nreact) print*, 'wrong dimensions in reactionrates : number of reaction should be ', nreact,' and not ', ireact
    692609 if ((nb_phot /= nb_phot_max)             .or.  &
    693610     (nb_reaction_3 /= nb_reaction_3_max) .or.  &
    … …  
    704621    print*, 'n4_hard_coding    = ', n4_hard_coding
    705622    print*, 'n3_hard_coding    = ', n3_hard_coding
    706     stop
     623    call abort
    707624 end if
    708625 if ((nb_hv /= nb_hv_max)         .or.  &
    … …  
    719636    print*, 'nb_pfunc2_max = ', nb_pfunc2_max
    720637    print*, 'nb_pfunc3_max = ', nb_pfunc3_max
    721     stop
     638    call abort
    722639 end if
    723640firstcall = .false.
    … …  
    843760!================================================================
    844761
    845  subroutine indice(nreact,rtype,third_body,three_prod)
     762 subroutine indice(nreact)
    846763
    847764!================================================================
    … …  
    869786
    870787integer, intent(out)              :: nreact          ! number of reactions in reactions files
    871 integer, allocatable, intent(out) :: rtype(:)        ! reaction rate type
    872 integer, allocatable, intent(out) :: third_body(:)   ! if the reaction have a third body: index of the third body, else zero
    873 logical, allocatable, intent(out) :: three_prod(:)   ! if the reaction have a three defferent proucts egal .true.
    874788
    875789! local
    … …  
    880794
    881795character(len = 128)              :: reactfile       ! reactions table file name
    882 character(len = 128)              :: monoreact       ! photochemical reactions table file name
    883 character(len = 128)              :: bimolreact      ! bimolecular reactions table file name
    884 character(len = 128)              :: quadrareact     ! quadratic reactions table file name
    885796integer                           :: nbq             ! number of species in reactions file
    886 integer                           :: pnlines         ! number of lines in photochemical reactions file
    887 integer                           :: bnlines         ! number of lines in bimolecular reactions file
    888 integer                           :: qnlines         ! number of lines in quadratic reactions file
     797integer                           :: nlines          ! number of lines in reactions file
    889798integer                           :: pnreact         ! number of reaction in photochemical reactions file
    890799integer                           :: bnreact         ! number of reaction in bimolecular reactions file
    … …  
    909818
    910819! Get number of reactions
    911 pnlines        = 0
    912 bnlines        = 0
    913 qnlines        = 0
     820nlines         = 0
    914821nreact         = 0
    915822pnreact        = 0
    … …  
    917824qnreact        = 0
    918825
    919 write(*,*) "Read photochemical reaction"
    920  monoreact = "monoreact" ! default
    921 call getin_p("monoreact",monoreact) ! default path
    922 write(*,*) " monoreact = ",trim(monoreact)
    923 
    924 write(*,*) "Read bimolecular reaction"
    925  bimolreact = "bimolreact" ! default
    926 call getin_p("bimolreact",bimolreact) ! default path
    927 write(*,*) " bimolreact = ",trim(bimolreact)
    928 
    929 write(*,*) "Read quadratic reaction"
    930  quadrareact = "quadrareact" ! default
    931 call getin_p("quadrareact",quadrareact) ! default path
    932 write(*,*) " quadrareact = ",trim(quadrareact)
    933 
    934 call count_react(monoreact,pnlines,pnreact,nb_phot,nb_hv,nb_pfunc1,nb_pfunc2,nb_pfunc3)
    935 call count_react(bimolreact,bnlines,bnreact,nb_reaction_4,nb_hv,nb_pfunc1,nb_pfunc2,nb_pfunc3)
    936 call count_react(quadrareact,qnlines,qnreact,nb_reaction_3,nb_hv,nb_pfunc1,nb_pfunc2,nb_pfunc3)
    937 nreact = pnreact + bnreact + qnreact
     826write(*,*) "Read reaction file"
     827reactfile = "reactfile" ! default
     828call getin_p("reactfile",reactfile) ! default path
     829write(*,*) " reactfile = ",trim(reactfile)
     830
     831call count_react(reactfile,nlines,nreact,nb_phot,nb_reaction_4,nb_reaction_3,nb_hv,nb_pfunc1,nb_pfunc2,nb_pfunc3)
    938832
    939833!!!!!!!!!!! Hard coding reaction !!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    964858allocate(indice_4(nb_reaction_4_max))
    965859allocate(indice_3(nb_reaction_3_max))
    966 allocate(rtype(nreact))
    967 allocate(third_body(nreact))
    968 allocate(three_prod(nreact))
     860allocate(reactions(nreact))
    969861allocate(jlabel(nb_phot_hv_max,2))
    970 allocate(jparamline(nb_phot_hv_max))
     862allocate(jparamline(nb_hv_max))
    971863allocate(pfunc1_param(nb_pfunc1_max))
    972864allocate(pfunc2_param(nb_pfunc2_max))
    … …  
    979871specheckr(:)    = 0
    980872specheckp(:)    = 0
    981 rtype(:)        = 0
    982 third_body(:)   = 0
     873reactions(:)    = reaction('',-1,0,.false.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.)
    983874pfunc1_param(:) = rtype1(0.,0.,0.,0.,0.)
    984875pfunc2_param(:) = rtype2(0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.)
    985876pfunc3_param(:) = rtype3(0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.)
    986 nb_pfunc1       = 0
    987 nb_pfunc2       = 0
    988 nb_pfunc3       = 0
    989 three_prod(:)   = .false.
     877nb_phot         = 0
     878nb_reaction_4   = 0
     879nb_reaction_3   = 0
    990880jlabel(:,:)     = ''
    991881jparamline(:)   = ''
    992882
    993 call get_react(monoreact,pnlines,pnreact,rtype(1:pnreact),third_body(1:pnreact),three_prod(1:pnreact),    &
    994                nb_phot,specheck,specheckr,specheckp,'vphot',nbq,nb_pfunc1,nb_pfunc2,nb_pfunc3)
    995 call get_react(bimolreact,bnlines,bnreact,rtype(pnreact+1:pnreact+bnreact),third_body(pnreact+1:pnreact+bnreact), &
    996                three_prod(pnreact+1:pnreact+bnreact),nb_reaction_4,specheck,specheckr,specheckp,'v4',nbq,nb_pfunc1,nb_pfunc2,nb_pfunc3)
    997 call get_react(quadrareact,qnlines,qnreact,rtype(pnreact+bnreact+1:nreact),third_body(pnreact+bnreact+1:nreact),  &
    998                three_prod(pnreact+bnreact+1:nreact),nb_reaction_3,specheck,specheckr,specheckp,'v3',nbq,nb_pfunc1,nb_pfunc2,nb_pfunc3)
     883call get_react(reactfile,nlines,nreact,specheck,specheckr,specheckp,nbq,nb_phot,nb_reaction_4,nb_reaction_3)
    999884
    1000885! rewrite jlabel corretly for output
    1001886do ireact=1,nb_phot_hv_max
    1002   if (three_prod(ireact)) then
     887  if (reactions(ireact)%three_prod) then
    1003888    jlabel(ireact+1:nb_phot_hv_max,2) = jlabel(ireact:nb_phot_hv_max-1,2)
    1004889    jlabel(ireact+1:nb_phot_hv_max,1) = jlabel(ireact:nb_phot_hv_max-1,1)
    … …  
    1031916    else
    1032917      print*, 'Error in indice: bad value in specheckr or specheckp'
    1033       stop
     918      call abort
    1034919    end if
    1035920  end if
    … …  
    1067952   print*, 'nb_phot_hv_max    = ', nb_phot_hv_max
    1068953   print*, 'nb_hv_max         = ', nb_hv_max
    1069    stop
     954   call abort
    1070955end if 
    1071956
    … …  
    11951080        do
    11961081            n = n + 1
     1082            if (n>nmax) then
     1083              print*,'Error in split_str: to much words'
     1084              print*,'limit = ',nmax
     1085              print*,'change it, if you want, with increasing nmax'
     1086              print*,'line is:',line
     1087              call abort
     1088            end if
    11971089            read(line,*,iostat=ios) buf(1:n)  ! use list-directed input
    11981090            if (ios==0) then
    … …  
    12021094                exit       ! if all the words are obtained, finish
    12031095            endif
    1204             if (n>nmax) then
    1205               print*,'Error in split_str: to much words'
    1206               print*,'limit = ',nmax
    1207               print*,'change it, if you want, with increasing nmax'
    1208               stop
    1209             end if
    12101096        enddo
    12111097    end subroutine split_str
    … …  
    12131099!*****************************************************************
    12141100
    1215     subroutine count_react(reactfile,nlines,nreact,nrtype,nb_hv,nb_pfunc1,nb_pfunc2,nb_pfunc3)
     1101    subroutine find_vtype(reactline,vtype)
    12161102
    12171103!*****************************************************************
    12181104
     1105        use tracer_h,     only: nesp
     1106        use chimiedata_h, only: indexchim
     1107
     1108        implicit none
     1109
     1110        character(len = 50), intent(in)    :: reactline ! all reactants of one reaction
     1111        character(*),        intent(inout) :: vtype     ! "v3", "v4" or "vphot"
     1112
     1113        integer             :: nwr                      ! number of reactant for a reaction
     1114        integer,parameter   :: nmax = 5                 ! number max of reactants or products
     1115        character(len = 24) :: words(nmax)              ! to get words in reactants and products line
     1116        integer             :: stochio(nesp)            ! stoichiometry of reaction
     1117        integer             :: iword
     1118
     1119        ! init
     1120        stochio(:) = 0
     1121
     1122        ! split reactant
     1123        call split_str(reactline,words,nwr,nmax)
     1124
     1125        ! set stochio
     1126        do iword = 1,nwr
     1127            if (trim(words(iword))=="M") exit
     1128            if (trim(words(iword))/="hv"         &
     1129                .and. trim(words(iword))/="dummy") stochio(indexchim(words(iword))) = stochio(indexchim(words(iword))) + 1
     1130        end do
     1131
     1132        ! set vtype
     1133        if (sum(stochio)==1) then
     1134            vtype = "vphot"
     1135        else if (sum(stochio)==2) then
     1136            if (any(stochio==2)) then
     1137                vtype = "v3"
     1138            else
     1139                vtype = "v4"
     1140            end if
     1141        else if (sum(stochio)==3) then
     1142            if (any(stochio==2)) then
     1143                vtype = "v4"
     1144            else
     1145                print*,'Error in photochemistry_asis (find_vtype):'
     1146                print*,'3 different reactants not implemented'
     1147                call abort
     1148            end if
     1149        else
     1150            print*,'Error in photochemistry_asis (find_vtype):'
     1151            print*,'more than 2 different reactants not implemented'
     1152            call abort
     1153        end if
     1154
     1155    end subroutine find_vtype
     1156
     1157!*****************************************************************
     1158
     1159    subroutine count_react(reactfile,nlines,nreact,nb_phot,nb_reaction_4,nb_reaction_3,nb_hv,nb_pfunc1,nb_pfunc2,nb_pfunc3)
     1160
     1161!*****************************************************************
     1162
    12191163        use types_asis, only : nb_phot_hv_max
    12201164
    12211165        implicit none
    1222         character(*), intent(in)     :: reactfile ! name of the file to read
    1223         integer,      intent(inout)  :: nlines    ! number of lines in the file
    1224         integer,      intent(out)    :: nreact    ! real number of reaction
    1225         integer,      intent(inout)  :: nrtype    ! number of reaction for calculation
    1226         integer,      intent(inout)  :: nb_hv, nb_pfunc1, nb_pfunc2, nb_pfunc3
     1166        character(*), intent(in)     :: reactfile     ! name of the file to read
     1167        integer,      intent(inout)  :: nlines        ! number of lines in the file
     1168        integer,      intent(out)    :: nreact        ! real number of reaction
     1169        integer,      intent(inout)  :: nb_phot       ! dimension of "vphot" reaction
     1170        integer,      intent(inout)  :: nb_reaction_4 ! dimension of "v4" reaction
     1171        integer,      intent(inout)  :: nb_reaction_3 ! dimension of "v3" reaction
     1172        integer,      intent(inout)  :: nb_hv         ! number of typerate 0 reaction
     1173        integer,      intent(inout)  :: nb_pfunc1     ! number of typerate 1 reaction
     1174        integer,      intent(inout)  :: nb_pfunc2     ! number of typerate 2 reaction
     1175        integer,      intent(inout)  :: nb_pfunc3     ! number of typerate 3 reaction
    12271176
    12281177        ! local
    1229         character(len = 50) :: reactline          ! all reactants of one reaction
    1230         character(len = 50) :: prodline           ! all produts of one reaction
    1231         integer             :: typerate           ! get the type of the rate reaction coefficient (pfunc_i)
    1232         integer             :: nwp                ! number of products for a reaction
    1233         integer,parameter   :: nmax = 5           ! number max of reactants or products
    1234         character(len = 24) :: words(nmax)        ! to get words in reactants and products line
     1178        character(len = 50) :: reactline              ! all reactants of one reaction
     1179        character(len = 50) :: prodline               ! all produts of one reaction
     1180        integer             :: typerate               ! get the type of the rate reaction coefficient (pfunc_i)
     1181        character(5)        :: vtype                  ! "v3", "v4" or "vphot"
     1182        integer             :: nwp                    ! number of products for a reaction
     1183        integer,parameter   :: nmax = 5               ! number max of reactants or products
     1184        character(len = 24) :: words(nmax)            ! to get words in reactants and products line
    12351185        integer             :: ierr
    12361186
    1237         nreact     = 0
     1187        nreact = 0
     1188        vtype  = ''
    12381189
    12391190        open(402, form = 'formatted', status = 'old',                &
    … …  
    12471198           write(*,*)'   monoreact=filename or bimolreact=filename or quadrareact=filename'
    12481199           write(*,*)'   depending on what chemical reaction type it is'
    1249            stop
     1200           call abort
    12501201        end if
    12511202   
    … …  
    12531204          read(402,'(A,A,I2)',iostat=ierr) reactline, prodline, typerate
    12541205          if (ierr<0) exit
     1206          ! count the number of lines
    12551207          nlines = nlines + 1
    12561208          if (reactline(1:1)/='!' .and. reactline(1:1)/='') then
    12571209            ! count the number of reaction
    12581210            nreact = nreact + 1
     1211
     1212            call find_vtype(reactline,vtype)
    12591213            call split_str(prodline,words,nwp,nmax)
    1260             nrtype = nrtype + 1
    1261             ! check three product reaction
    1262             if (nwp>2 .and. trim(words(1))/=trim(words(2)) &
    1263                       .and. trim(words(1))/=trim(words(3)) &
    1264                       .and. trim(words(2))/=trim(words(3))) nrtype = nrtype + 1
     1214
     1215            ! count the dimension of each rate reaction coefficient type
     1216            if (trim(vtype)=="vphot") then
     1217                nb_phot = nb_phot + 1
     1218                ! check three product reaction
     1219                if (nwp>2 .and. trim(words(1))/=trim(words(2)) &
     1220                          .and. trim(words(1))/=trim(words(3)) &
     1221                          .and. trim(words(2))/=trim(words(3))) nb_phot = nb_phot + 1
     1222            else if (trim(vtype)=="v4") then
     1223                nb_reaction_4 = nb_reaction_4 + 1
     1224                ! check three product reaction
     1225                if (nwp>2 .and. trim(words(1))/=trim(words(2)) &
     1226                          .and. trim(words(1))/=trim(words(3)) &
     1227                          .and. trim(words(2))/=trim(words(3))) nb_reaction_4 = nb_reaction_4 + 1
     1228            else if (trim(vtype)=="v3") then
     1229                nb_reaction_3 = nb_reaction_3 + 1
     1230                ! check three product reaction
     1231                if (nwp>2 .and. trim(words(1))/=trim(words(2)) &
     1232                          .and. trim(words(1))/=trim(words(3)) &
     1233                          .and. trim(words(2))/=trim(words(3))) nb_reaction_3 = nb_reaction_3 + 1
     1234            else
     1235                print*,'Error in photochemistry_asis (count_react):'
     1236                print*,'vtype not found'
     1237                call abort
     1238            end if
    12651239
    12661240            ! count the number of each rate reaction coefficient type
    … …  
    12821256              print*, 'It should be 0 for photolysis reations and 1 or 2 for the others'
    12831257              print*, 'And not : ', typerate
    1284               stop
     1258              call abort
    12851259            end if
    12861260
    … …  
    12951269!*****************************************************************
    12961270
    1297     subroutine get_react(reactfile,nlines,nreact,rtype,third_body,three_prod, &
    1298                          nrtype,specheck,specheckr,specheckp,typeindice,nbq,  &
    1299                          init_nb_pfunc1,init_nb_pfunc2,init_nb_pfunc3)
     1271    subroutine get_react(reactfile,nlines,nreact,specheck,specheckr,specheckp, &
     1272                         nbq,nb_phot,nb_reaction_4,nb_reaction_3)
    13001273
    13011274!*****************************************************************
    … …  
    13101283        integer,      intent(in)     :: nlines             ! number of lines in the file
    13111284        integer,      intent(in)     :: nreact             ! real number of reaction in the file
    1312         integer,      intent(inout)  :: rtype(nreact)      ! reaction rate type
    1313         integer,      intent(inout)  :: third_body(nreact) ! if the reaction have a third body: index of the third body, else zero
    1314         logical,      intent(inout)  :: three_prod(nreact) ! if the reaction have a three defferent proucts egal .true.
    1315         integer,      intent(out)    :: nrtype             ! number of calculation reactions
     1285        integer,      intent(out)    :: nb_phot            ! number of "vphot" calculation reactions
     1286        integer,      intent(out)    :: nb_reaction_4      ! number of "vphot" calculation reactions
     1287        integer,      intent(out)    :: nb_reaction_3      ! number of "vphot" calculation reactions
    13161288        integer,      intent(inout)  :: specheck(nesp)     ! to count the species of traceur.def in reactions file
    13171289        integer,      intent(inout)  :: specheckr(nesp)    ! to count the reactant species of traceur.def in reactions file
    13181290        integer,      intent(inout)  :: specheckp(nesp)    ! to count the product species of traceur.def in reactions file
    1319         character(*), intent(in)     :: typeindice         ! reaction type (v3, v4 or vphot)
    13201291        integer,      intent(inout)  :: nbq                ! number of species in reactions file
    1321         integer,      intent(inout)  :: init_nb_pfunc1     ! in : initial value of nb_pfunc1 - out : final value of nb_pfunc1
    1322         integer,      intent(inout)  :: init_nb_pfunc2     ! in : initial value of nb_pfunc2 - out : final value of nb_pfunc2
    1323         integer,      intent(inout)  :: init_nb_pfunc3     ! in : initial value of nb_pfunc3 - out : final value of nb_pfunc3
    13241292
    13251293        ! local
    … …  
    13331301        integer              :: iindice(2*nmax)            ! indice of species (for indice variables)
    13341302        character(len = 500) :: paramline                  ! line of reactions parameters
    1335         character(len = 50)  :: reactants(nreact,nmax)     ! reactions reactants
    1336         character(len = 50)  :: products(nreact,nmax)      ! reactions products
    1337         logical              :: spedouble                  ! check if a specie appears twice in reactants or products
    1338         integer              :: ierr, ilines, ireact, i_dummy, iw, iwhere, i
    1339         integer              :: nb_hv, nb_pfunc1, nb_pfunc2, nb_pfunc3
    1340 
    1341         i_dummy   = 1
    1342         nrtype    = 0
    1343         ireact    = 0
    1344         nb_hv     = 0
    1345         nb_pfunc1 = init_nb_pfunc1
    1346         nb_pfunc2 = init_nb_pfunc2
    1347         nb_pfunc3 = init_nb_pfunc3
     1303        integer              :: stochio(nesp)              ! stoichiometry of reaction
     1304        integer              :: ierr, ilines, ireact, i_dummy, iw, ispe, iloc
     1305        integer              :: nb_phot_hv, nb_hv, nb_pfunc1, nb_pfunc2, nb_pfunc3
     1306
     1307        i_dummy       = 1
     1308        ireact        = 0
     1309        nb_phot       = 0
     1310        nb_phot_hv    = 0
     1311        nb_hv         = 0
     1312        nb_pfunc1     = 0
     1313        nb_pfunc2     = 0
     1314        nb_pfunc3     = 0
    13481315
    13491316        open(402, form = 'formatted', status = 'old',                &
    … …  
    13571324           write(*,*)'   monoreact=filename or bimolreact=filename or quadrareact=filename'
    13581325           write(*,*)'   depending on what chemical reaction type it is'
    1359            stop
     1326           call abort
    13601327        end if
    13611328
    … …  
    13681335          if (reactline(1:1)/='!' .and. reactline(1:1)/='') then
    13691336
    1370             ! increment number of reactions and init
     1337            ! increment number of reactions
    13711338            ireact = ireact + 1
    1372         !!!!!!!!!!!!!!!!!!!!!!!!!!! for fill indice part
    1373             nrtype = nrtype + 1
    1374             nindice(:) = 0.0
    1375             iindice(:) = i_dummy
    1376         !!!!!!!!!!!!!!!!!!!!!!!!!!! end
    1377             ! get indice, rate type and parameters
     1339
     1340            ! fill reaction vtype
     1341            call find_vtype(reactline,reactions(ireact)%vtype)
     1342            if (trim(reactions(ireact)%vtype)=="v4") then
     1343                nb_reaction_4 = nb_reaction_4 + 1
     1344            else if (trim(reactions(ireact)%vtype)=="v3") then
     1345                nb_reaction_3 = nb_reaction_3 + 1
     1346            else if (trim(reactions(ireact)%vtype)=="vphot") then
     1347                nb_phot = nb_phot + 1
     1348            else
     1349                print*,'Error in photochemistry_asis (get_react):'
     1350                print*,'vtype not found'
     1351                call abort
     1352            end if
     1353
     1354            ! fill reaction rate type and parameters
    13781355            if (trim(paramline)=='') then
    1379               print*, 'Error in reactfile: where are the parameters - line',ilines
    1380               stop
     1356                print*, 'Error in reactfile: where are the parameters - line',ilines
     1357                call abort
    13811358            else
    1382               read(paramline,*) rtype(ireact)
    1383               if (rtype(ireact)==1) then
    1384                 nb_pfunc1 = nb_pfunc1 + 1
    1385                 read(paramline,*) rtype(ireact), pfunc1_param(nb_pfunc1)
    1386               else if (rtype(ireact)==0) then
    1387                 nb_hv = nb_hv + 1
    1388                 if (jonline) then
    1389                   read(paramline,'(I5,A,A)') rtype(ireact), jlabel(nb_hv,1), jparamline(nb_hv)
     1359                read(paramline,*) reactions(ireact)%rtype
     1360                if (reactions(ireact)%rtype==1) then
     1361                  nb_pfunc1 = nb_pfunc1 + 1
     1362                  read(paramline,*) reactions(ireact)%rtype, pfunc1_param(nb_pfunc1)
     1363                else if (reactions(ireact)%rtype==0) then
     1364                  nb_hv = nb_hv + 1
     1365                  nb_phot_hv = nb_phot_hv + 1
     1366                  if (jonline) then
     1367                    read(paramline,'(I5,A,A)') reactions(ireact)%rtype, jlabel(nb_hv,1), jparamline(nb_hv)
     1368                  else
     1369                    read(paramline,*) reactions(ireact)%rtype, jlabel(nb_hv,1)
     1370                  end if
     1371                else if (reactions(ireact)%rtype==2) then
     1372                  nb_pfunc2 = nb_pfunc2 + 1
     1373                  read(paramline,*) reactions(ireact)%rtype, pfunc2_param(nb_pfunc2)
     1374                else if (reactions(ireact)%rtype==3) then
     1375                  nb_pfunc3 = nb_pfunc3 + 1
     1376                  read(paramline,*) reactions(ireact)%rtype, pfunc3_param(nb_pfunc3)
     1377                end if
     1378            end if
     1379
     1380            ! WARNING: the implementation is limited to 3 different reactants (for now only 2) and 3 different products
     1381            nindice(:) = 0.0     ! 3 first indice for reactants, 3 following for products
     1382            iindice(:) = i_dummy ! 3 first indice for reactants, 3 following for products
     1383
     1384            !-----------------!
     1385            !--- reactants ---!
     1386            !-----------------!
     1387
     1388            ! init for reactant
     1389            stochio(:) = 0
     1390            ! split reactant
     1391            call split_str(reactline,words,nwr,nmax)
     1392            ! set species that are photolysed
     1393            if (reactions(ireact)%rtype==0) jlabel(nb_hv,2) = words(1)
     1394            ! find reactant stochio
     1395            do iw = 1,nwr
     1396                ! fill third body index
     1397                if (trim(words(iw))=="M") then
     1398                    if (iw==nwr) then
     1399                        exit
     1400                    else if (iw==nwr-1) then
     1401                        reactions(ireact)%third_body = indexchim(words(iw+1))
     1402                        exit
     1403                    else
     1404                        print*, 'Error in reactfile: just only one specie can be after M corresponding to the third body - line',ilines
     1405                        call abort
     1406                    end if
     1407                end if
     1408                ! count stochio
     1409                if (trim(words(iw))/="hv"         &
     1410                    .and. trim(words(iw))/="dummy") stochio(indexchim(words(iw))) = stochio(indexchim(words(iw))) + 1
     1411            end do
     1412            ! get reaction reactant stochio and indice
     1413            iloc = 0
     1414            do ispe = 1,nesp
     1415                if (stochio(ispe)==0) cycle
     1416                iloc = iloc + 1
     1417                nindice(iloc) = stochio(ispe)
     1418                iindice(iloc) = ispe
     1419                ! check up: to count the species used in 'reactfile'
     1420                if (specheck(ispe)==0) then
     1421                    specheckr(ispe) = 1
     1422                    specheck(ispe) = 1
     1423                    nbq = nbq + 1
     1424                else if (specheckr(ispe)==0) then
     1425                    specheckr(ispe) = 1
     1426                end if
     1427            end do
     1428            ! fill reaction reactant stochio and indice
     1429            reactions(ireact)%ir1 = nindice(1)
     1430            reactions(ireact)%r1  = iindice(1)
     1431            reactions(ireact)%ir2 = nindice(2)
     1432            reactions(ireact)%r2  = iindice(2)
     1433            reactions(ireact)%ir3 = nindice(3)
     1434            reactions(ireact)%r3  = iindice(3)
     1435
     1436            !----------------!
     1437            !--- products ---!
     1438            !----------------!
     1439
     1440            ! init for products
     1441            stochio(:) = 0
     1442            ! split products
     1443            call split_str(prodline,words,nwp,nmax)
     1444            ! find products stochio
     1445            do iw = 1,nwp
     1446                stochio(indexchim(words(iw))) = stochio(indexchim(words(iw))) + 1
     1447            end do
     1448            ! get reaction product stochio and indice
     1449            iloc = 3
     1450            do ispe = 1,nesp
     1451                if (stochio(ispe)==0) cycle
     1452                iloc = iloc + 1
     1453                nindice(iloc) = stochio(ispe)
     1454                iindice(iloc) = ispe
     1455                ! check up: to count the species used in 'reactfile'
     1456                if (specheck(ispe)==0) then
     1457                    specheckp(ispe) = 1
     1458                    specheck(ispe) = 1
     1459                    nbq = nbq + 1
     1460                else if (specheckp(ispe)==0) then
     1461                    specheckp(ispe) = 1
     1462                end if
     1463            end do
     1464            ! fill reaction product stochio and indice
     1465            reactions(ireact)%ip1 = nindice(4)
     1466            reactions(ireact)%p1  = iindice(4)
     1467            reactions(ireact)%ip2 = nindice(5)
     1468            reactions(ireact)%p2  = iindice(5)
     1469            reactions(ireact)%ip3 = nindice(6)
     1470            reactions(ireact)%p3  = iindice(6)
     1471            ! set reaction three prod to true with checking in position 6 if there is three different products
     1472            if (nindice(6)/=0) reactions(ireact)%three_prod = .true.
     1473
     1474
     1475            !-------------------------!
     1476            !--- fill indice array ---!
     1477            !-------------------------!
     1478
     1479            if (trim(reactions(ireact)%vtype)=="v4") then
     1480                if (nindice(6)==0) then ! reaction with 1 or 2 products
     1481                    ! z4spec(ir1,r1,ir2,r2,ip1,p1,ip2,p2)
     1482                    indice_4(nb_reaction_4) = z4spec(nindice(1), iindice(1), nindice(2), iindice(2), nindice(4), iindice(4), nindice(5), iindice(5))
     1483                else ! reaction with 3 different products
     1484                    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))
     1485                    nb_reaction_4 = nb_reaction_4 + 1
     1486                    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))
     1487                end if
     1488            else if (trim(reactions(ireact)%vtype)=="v3") then
     1489                if (nindice(6)==0) then ! reaction with 1 or 2 products
     1490                    ! z3spec(ir1,r1,ip1,p1,ip2,p2)
     1491                    indice_3(nb_reaction_3) = z3spec(nindice(1), iindice(1), nindice(4), iindice(4), nindice(5), iindice(5))
     1492                else ! reaction with 3 different products
     1493                    indice_3(nb_reaction_3) = z3spec(nindice(1)/2., iindice(1), nindice(4), iindice(4), 0.0, i_dummy)
     1494                    nb_reaction_3 = nb_reaction_3 + 1
     1495                    indice_3(nb_reaction_3) = z3spec(nindice(1)/2., iindice(1), nindice(5), iindice(5), nindice(6), iindice(6))
     1496                end if
     1497            else if (trim(reactions(ireact)%vtype)=="vphot") then
     1498                if (reactions(ireact)%rtype==0) then
     1499                    if (nindice(6)==0) then ! reaction with 1 or 2 products
     1500                        ! z3spec(ir1,r1,ip1,p1,ip2,p2)
     1501                        indice_phot(nb_phot_hv) = z3spec(nindice(1), iindice(1), nindice(4), iindice(4), nindice(5), iindice(5))
     1502                    else ! reaction with 3 different products
     1503                        indice_phot(nb_phot_hv) = z3spec(nindice(1)/2., iindice(1), nindice(4), iindice(4), 0.0, i_dummy)
     1504                        nb_phot = nb_phot + 1
     1505                        nb_phot_hv = nb_phot_hv + 1
     1506                        indice_phot(nb_phot_hv) = z3spec(nindice(1)/2., iindice(1), nindice(5), iindice(5), nindice(6), iindice(6))
     1507                    end if
    13901508                else
    1391                   read(paramline,*) rtype(ireact), jlabel(nb_hv,1)
     1509                    if (nindice(6)==0) then ! reaction with 1 or 2 products
     1510                        ! z3spec(ir1,r1,ip1,p1,ip2,p2)
     1511                        indice_phot(nb_phot_hv_max+nb_phot-nb_phot_hv) = z3spec(nindice(1), iindice(1), nindice(4), iindice(4), nindice(5), iindice(5))
     1512                    else ! reaction with 3 different products
     1513                        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)
     1514                        nb_phot = nb_phot + 1
     1515                        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))
     1516                    end if
    13921517                end if
    1393               else if (rtype(ireact)==2) then
    1394                 nb_pfunc2 = nb_pfunc2 + 1
    1395                 read(paramline,*) rtype(ireact), pfunc2_param(nb_pfunc2)
    1396               else if (rtype(ireact)==3) then
    1397                 nb_pfunc3 = nb_pfunc3 + 1
    1398                 read(paramline,*) rtype(ireact), pfunc3_param(nb_pfunc3)
    1399               end if
     1518            else
     1519                print*,'Error in photochemistry_asis (get_react):'
     1520                print*,'vtype',reactions(ireact)%vtype,' not implemented'
     1521                call abort
    14001522            end if
    14011523       
    1402             ! get reactants
    1403             call split_str(reactline,words,nwr,nmax)
    1404             if (rtype(ireact)==0) jlabel(nb_hv,2) = words(1)
    1405             ! loop on reactants
    1406             do iw=1,nwr
    1407               ! store reactants in variable 'reactants'
    1408               reactants(ireact,iw) = trim(words(iw))
    1409               ! check third body and exit reactants loop
    1410               if (reactants(ireact,iw)=='M') then
    1411                 if (iw==nwr) then
    1412                   exit
    1413                 else if (iw==nwr-1) then
    1414                   third_body(ireact) = indexchim(words(iw+1))
    1415                   exit
    1416                 else
    1417                   print*, 'Error in reactfile: just only one specie can be after M corresponding to the third body - line',ilines
    1418                   stop
    1419                 end if
    1420               end if
    1421               if (trim(words(iw))/='hv' .and. trim(words(iw))/='dummy') then
    1422                 iwhere = indexchim(words(iw))
    1423                 ! check if species are chemical tracers
    1424                 if (iwhere>nesp) then
    1425                   print*, 'Error: in ', trim(reactfile)
    1426                   print*, 'check if the specie', trim(words(iw)),' is include into chemical tracers in traceur.def'
    1427                   stop
    1428                 end if
    1429                 ! to count the species used in 'reactfile'
    1430                 if (specheck(iwhere)==0) then
    1431                   specheckr(iwhere) = 1
    1432                   specheck(iwhere) = 1
    1433                   nbq = nbq + 1
    1434                 else if (specheckr(iwhere)==0) then
    1435                   specheckr(iwhere) = 1
    1436                 end if
    1437        
    1438         !!!!!!!!!!!!!!!!!!! for fill indice part
    1439                 ! fill stochiometry and indice of rection species depending of reaction type
    1440                 if (trim(typeindice)=='v3') then
    1441                   nindice(1) = 2.0
    1442                   iindice(1) = indexchim(words(iw))
    1443                   if (nwr>3 .or. nwr<2) print*, 'Error in reactfile: wrong number of reactants for v3 reaction line',ilines
    1444                   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
    1445                 else if (trim(typeindice)=='v4') then
    1446                   nindice(iw) = 1.0
    1447                   iindice(iw) = indexchim(words(iw))
    1448                 else if (trim(typeindice)=='vphot') then
    1449                   nindice(1) = 1.0
    1450                   if (iw>2) then
    1451                     print*, 'Something weird in your photolysis reaction'
    1452                     print*, 'You should have 1 reactants and hv'
    1453                     print*, 'Reactants are: ',words
    1454                     stop
    1455                   end if
    1456                   iindice(1) = indexchim(words(iw))
    1457                 end if
    1458         !!!!!!!!!!!!!!!!!!! end
    1459        
    1460               end if
    1461             end do
    1462 
    1463             ! same as reactants but for the products
    1464             call split_str(prodline,words,nwp,nmax)
    1465             do iw=1,nwp
    1466               spedouble = .false.
    1467               products(ireact,iw) = trim(words(iw))
    1468               if (trim(words(iw))/='hv' .and. trim(words(iw))/='dummy' .and. trim(words(iw))/='M') then
    1469                 iwhere = indexchim(words(iw))
    1470                 if (iwhere>nesp) then
    1471                   print*, 'Error: in ', trim(reactfile)
    1472                   print*, 'check if the specie', trim(words(iw)),' is include into chemical tracers in traceur.def'
    1473                   stop
    1474                 end if
    1475                 if (specheck(iwhere)==0) then
    1476                   specheckp(iwhere) = 1
    1477                   specheck(iwhere) = 1
    1478                   nbq = nbq + 1
    1479                 else if (specheckp(iwhere)==0) then
    1480                   specheckp(iwhere) = 1
    1481                 end if
    1482        
    1483         !!!!!!!!!!!!!!!!!!!!!!!!!! for fill indice part
    1484                 if (trim(typeindice)=='v3' .or. trim(typeindice)=='vphot') then
    1485                     iindice(1+iw) = indexchim(words(iw))
    1486                   do i=1,iw-1
    1487                     if (iindice(1+iw)==iindice(1+i)) then
    1488                       nindice(1+i)  = nindice(1+i) + 1.0
    1489                       iindice(1+iw) = i_dummy
    1490                       spedouble = .true.
    1491                     end if
    1492                   end do
    1493                   if (.not. spedouble) nindice(1+iw) = 1.0
    1494                 else if (trim(typeindice)=='v4') then
    1495                   iindice(2+iw) = indexchim(words(iw))
    1496                   do i=1,iw-1
    1497                     if (iindice(2+iw)==iindice(2+i)) then
    1498                       nindice(2+i)  = nindice(2+i) + 1.0
    1499                       iindice(2+iw) = i_dummy
    1500                       spedouble = .true.
    1501                     end if
    1502                   end do
    1503                   if (.not. spedouble) nindice(2+iw) = 1.0
    1504                 end if
    1505         !!!!!!!!!!!!!!!!!!!!!!!!!!! end
    1506               else
    1507                 print*, 'Error: no hv, dummy or M in products'
    1508                 stop
    1509               end if
    1510             end do
    1511 
    1512             ! fill indice variables
    1513             if (trim(typeindice)=='v3') then
    1514               if (nindice(4)/=0.0) then ! reaction with 3 products
    1515                 if (nindice(3)==0.0) then ! 2 are the same species
    1516                   indice_3(nrtype) = z3spec(nindice(1), iindice(1), nindice(2), iindice(2), nindice(4), iindice(4))
    1517                 else ! reaction with 3 different products
    1518                   indice_3(nrtype) = z3spec(nindice(1)/2., iindice(1), nindice(2), iindice(2), 0.0, i_dummy)
    1519                   nrtype = nrtype + 1
    1520                   indice_3(nrtype) = z3spec(nindice(1)/2., iindice(1), nindice(3), iindice(3), nindice(4), iindice(4))
    1521                   three_prod(ireact) = .true.
    1522                 end if
    1523               else ! reaction with 1 or 2 products
    1524                 indice_3(nrtype) = z3spec(nindice(1), iindice(1), nindice(2), iindice(2), nindice(3), iindice(3))
    1525               end if
    1526             else if (trim(typeindice)=='v4') then
    1527               if (nindice(5)/=0.0) then ! reaction with 3 products
    1528                 if (nindice(4)==0.0) then ! 2 are the same species
    1529                   indice_4(nrtype) = z4spec(nindice(1), iindice(1), nindice(2), iindice(2), nindice(3), iindice(3), nindice(5), iindice(5))
    1530                 else ! reaction with 3 different products
    1531                   indice_4(nrtype) = z4spec(nindice(1)/2., iindice(1), nindice(2)/2., iindice(2), nindice(3), iindice(3), nindice(4)/2., iindice(4))
    1532                   nrtype = nrtype + 1
    1533                   indice_4(nrtype) = z4spec(nindice(1)/2., iindice(1), nindice(2)/2., iindice(2), nindice(5), iindice(5), nindice(4)/2., iindice(4))
    1534                   three_prod(ireact) = .true.
    1535                 end if
    1536               else ! reaction with 1 or 2 products
    1537                 indice_4(nrtype) = z4spec(nindice(1), iindice(1), nindice(2), iindice(2), nindice(3), iindice(3), nindice(4), iindice(4))
    1538               end if
    1539             else if (trim(typeindice)=='vphot') then
    1540               if (nindice(4)/=0.0) then ! reaction with 3 products
    1541                 if (nindice(3)==0.0) then ! 2 are the same species
    1542                   indice_phot(nrtype) = z3spec(nindice(1), iindice(1), nindice(2), iindice(2), nindice(4), iindice(4))
    1543                 else ! reaction with 3 different products
    1544                   indice_phot(nrtype) = z3spec(nindice(1)/2., iindice(1), nindice(2), iindice(2), 0.0, i_dummy)
    1545                   nrtype = nrtype + 1
    1546                   indice_phot(nrtype) = z3spec(nindice(1)/2., iindice(1), nindice(3), iindice(3), nindice(4), iindice(4))
    1547                   three_prod(ireact) = .true.
    1548                 end if
    1549               else ! reaction with 1 or 2 products
    1550                 indice_phot(nrtype) = z3spec(nindice(1), iindice(1), nindice(2), iindice(2), nindice(3), iindice(3))
    1551               end if
    1552             end if
    1553        
    1554           end if
    1555 
    1556         end do
    1557 
    1558         init_nb_pfunc1 = nb_pfunc1
    1559         init_nb_pfunc2 = nb_pfunc2
    1560         init_nb_pfunc3 = nb_pfunc3
     1524          end if ! end if comment line
     1525
     1526        end do ! end loop on lines
    15611527
    15621528        close(402)

  • trunk/LMDZ.GENERIC/libf/aeronostd/photolysis_mod.F90

    r2542 r2553  
    4545    contains
    4646 
    47     subroutine init_photolysis(nlayer,nreact,three_prod)
    48  
     47    subroutine init_photolysis(nlayer,nreact)
     48 
     49    use types_asis, only: reactions
    4950    use tracer_h
    5051    use chimiedata_h, only: indexchim
    … …  
    5354    integer, intent(in) :: nlayer              ! number of atmospheric layers
    5455    integer, intent(in) :: nreact              ! number of reactions in reactions files
    55     logical, intent(in) :: three_prod(nreact)  ! if the reaction have a three defferent products egal .true.
    56     integer :: iphot,tdim_max,i3prod,nd,ij,iw,ilay
     56    integer :: iphot,tdim_max,i3prod,ij,iw,ilay,ireact
    5757    integer :: specheck(nesp)
    5858 
    … …  
    9494      endif
    9595    end do
    96     nd = nb_hv_max
    9796 
    9897    ! Get temperature dimension and allocate
    99     allocate(tdim(nd))
     98    allocate(tdim(nb_hv_max))
    10099    tdim(:)  = 1
    101100    tdim_max = 1
    102     do iphot=1,nd
     101    do iphot=1,nb_hv_max
    103102      read(jparamline(iphot),*) tdim(iphot)
    104103      if (tdim(iphot) > tdim_max) tdim_max = tdim(iphot)
    … …  
    106105 
    107106    ! allocation
    108     allocate(qy(nw,nd))
    109     allocate(xs(tdim_max,nw,nd))
     107    allocate(qy(nw,nb_hv_max))
     108    allocate(xs(tdim_max,nw,nb_hv_max))
    110109    allocate(sj(nlayer,nw,nb_phot_hv_max))
    111     allocate(xs_temp(tdim_max,nd))
     110    allocate(xs_temp(tdim_max,nb_hv_max))
    112111    xs(:,:,:)    = 0.
    113112    xs_temp(:,:) = 0.
    … …  
    119118 
    120119    i3prod = 0
    121  
    122     ! read and grid cross-sections (photolysis should be first reactions in monoreact file)
    123     do iphot=1,nd
     120    ireact = 0
     121 
     122    ! read and grid cross-sections
     123    do iphot=1,nb_hv_max
     124      ireact = ireact + 1
    124125      call rdxsi(nw,wl,xs(:tdim(iphot),:,iphot),jparamline(iphot),xs_temp(:tdim(iphot),iphot),qy(:,iphot),tdim(iphot),jlabel(iphot+i3prod,2))
    125       if (three_prod(iphot)) i3prod = i3prod + 1
     126      do while(reactions(ireact)%rtype/=0)
     127        ireact = ireact + 1
     128      end do
     129      if (reactions(ireact)%three_prod) i3prod = i3prod + 1
    126130    end do
    127131 
    128132    ! init sj for no temperature dependent cross sections (tdim.eq.1)
    129     iphot = 0
    130     ij    = 0
     133    iphot  = 0
     134    ij     = 0
     135    ireact = 0
    131136    do while(iphot<nb_phot_hv_max)
    132        ij = ij + 1
    133        iphot = iphot + 1
     137       ij     = ij     + 1
     138       iphot  = iphot  + 1
     139       ireact = ireact + 1
    134140       if (tdim(ij).eq.1) then
    135141         do iw = 1,nw-1
    … …  
    139145         end do
    140146       endif
    141        if (three_prod(ij)) then
     147       do while(reactions(ireact)%rtype/=0)
     148         ireact = ireact + 1
     149       end do
     150       if (reactions(ireact)%three_prod) then
    142151         iphot = iphot + 1
    143152         if (tdim(ij).eq.1) then
    … …  
    882891        print*, 'ERROR: temperature cross section file'
    883892        print*, '       has to be sorted from low to high values'
    884         print*, '       Check monoreact file'
     893        print*, '       Check reactfile file'
    885894        stop
    886895      end if
    … …  
    903912         write(*,*)'2) You can check if the file is in datadir/cross_sections/'
    904913         write(*,*)'3) You can change the input cross_sections file name in'
    905          write(*,*)'   chemestry/monoreact file for the specie:',trim(species)
     914         write(*,*)'   chemestry/reactfile file for the specie:',trim(species)
    906915         stop
    907916      end if
    … …  
    960969       write(*,*)'2) You can check if the file is in datadir/cross_sections/'
    961970       write(*,*)'3) You can change the input photodissociation yield file name in'
    962        write(*,*)'   chemestry/monoreact file for the specie:',trim(species)
     971       write(*,*)'   chemestry/reactfile file for the specie:',trim(species)
    963972       stop
    964973    end if

  • trunk/LMDZ.GENERIC/libf/aeronostd/photolysis_online.F

    r2542 r2553  
    44     $           zmmean, rm,
    55     $           tau, sza, dist_sol, v_phot, e_phot, ig, ngrid,
    6      $           nreact, three_prod)
     6     $           nreact)
    77
    88!***********************************************************************
    … …  
    2121      use tracer_h
    2222      use chimiedata_h, only: indexchim
    23       use types_asis,   only: nb_phot_hv_max, nb_phot_max, jlabel
     23      use types_asis,   only: nb_phot_hv_max, nb_phot_max, jlabel,
     24     $                        reactions
    2425
    2526      implicit none
    … …  
    3839      integer, intent(in) :: ig                                     ! grid point index
    3940      integer, intent(in) :: nreact                                 ! number of reactions in reactions files
    40       logical, intent(in) :: three_prod(nreact)                     ! if the reaction have a three defferent products egal .true.
    4141
    4242!     output
    … …  
    7171      real,    dimension(0:nlayer,nlayer) :: dsdh
    7272
    73       integer :: i, ilay, iw, ialt, iphot, ispe, ij, igas
     73      integer :: i, ilay, iw, ialt, iphot, ispe, ij, igas, ireact
    7474      real    :: deltaj
    7575
    … …  
    8585!==== set temperature-dependent cross-sections and optical depths
    8686
    87       iphot = 0
    88       ij    = 0
    89       igas  = 0
     87      iphot  = 0
     88      ij     = 0
     89      igas   = 0
     90      ireact = 0
    9091
    9192      do while(iphot<nb_phot_hv_max)
    92         ij = ij + 1
    93         iphot = iphot + 1
     93        ij     = ij     + 1
     94        iphot  = iphot  + 1
     95        ireact = ireact + 1
    9496
    9597        if (tdim(ij).eq.1) then
    … …  
    106108            igas = igas + 1
    107109          end if
    108           if (three_prod(ij)) iphot = iphot + 1
     110          do while(reactions(ireact)%rtype/=0)
     111            ireact = ireact + 1
     112          end do
     113          if (reactions(ireact)%three_prod) iphot = iphot + 1
    109114        else
    110115          call setsj(nb_phot_hv_max,nlayer,nw,temp,tdim(ij),
    … …  
    131136          end do
    132137
    133           if (three_prod(ij)) then
     138          do while(reactions(ireact)%rtype/=0)
     139            ireact = ireact + 1
     140          end do
     141          if (reactions(ireact)%three_prod) then
    134142            iphot = iphot + 1
    135143            do iw = 1,nw-1

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

    r2542 r2553  
    8585end type rtype3
    8686
     87! reaction type
     88
     89type reaction
     90    character(5)    :: vtype
     91    integer         :: rtype
     92    integer         :: third_body
     93    logical         :: three_prod
     94    integer         :: ir1
     95    integer         :: r1
     96    integer         :: ir2
     97    integer         :: r2
     98    integer         :: ir3
     99    integer         :: r3
     100    integer         :: ip1
     101    integer         :: p1
     102    integer         :: ip2
     103    integer         :: p2
     104    integer         :: ip3
     105    integer         :: p3
     106end type reaction
     107
    87108! pfunc parameters for the reaction rates
    88109
    … …  
    90111type(rtype2), allocatable, save :: pfunc2_param(:)
    91112type(rtype3), allocatable, save :: pfunc3_param(:)
     113
     114! reactions attributs
     115
     116type(reaction), allocatable, save :: reactions(:)
    92117
    93118! dimension of pfunc type variables
Note: See TracChangeset for help on using the changeset viewer.