ecearth3/sources/tm5mp/proj/cb05/ebischeme.F90

#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') rname, __FILE__, __LINE__ ; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      EBISCHEME
!
! !DESCRIPTION: Eulerian Backward Iteration (EBI) is a chemistry solver for
!               the modified CB05 scheme.
!\\
!\\
! !INTERFACE: 
!
module EBISCHEME
  !
  ! !USES:
  !
  use GO,              only : gol, goPr, goErr
  use dims,            only : nregions
  use tm5_distgrid,    only : dgrid, Get_DistGrid
  use chem_param     
#ifdef with_budgets
  use budget_global,   only : nbudg, nbud_vg
  use reaction_data,   only : nreac, nreacw
  use photolysis_data, only : nj
#endif

  implicit none
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: EBI
  !
  ! !PUBLIC DATA MEMBERS:
  !
#ifdef with_budgets
  real, dimension(nregions), public                ::  sum_wetchem
  real, dimension(nregions), public                ::  sum_chemistry
  real, dimension(nregions), public                ::  sum_deposition

  type, public :: buddep_data
     real, dimension(:,:,:), pointer               :: dry ! i,j,ntrace
  end type buddep_data
  
  type(buddep_data), dimension(nregions),  public  :: buddep_dat
                                                   
  real, dimension(nbudg, nbud_vg, nreac ), public  :: budrrg
  real, dimension(nbudg, nbud_vg, nj    ), public  :: budrjg
  real, dimension(nbudg, nbud_vg, nreacw), public  :: budrwg
  real, dimension(nbudg, nbud_vg, nmark ), public  :: budmarkg
  
  !allow saving of production of aqueous, gas-phase sulfate, elvoc, svoc

  type,public :: production_data
     real, dimension(:,:,:,:),pointer                   :: prod !i,j,lev,nprod
  end type production_data

  integer,parameter,public :: nprod_AC_o3=2
  integer,parameter,public :: nprod=12
  type(production_data), dimension(nregions),public:: diag_prod
  !The chemical production container for AerChemMIP output, needs own container, since the variable is zeroed when 
  ! writing out and could possibly conflict with the general output routine. 
   integer,parameter,public :: nprod_AC=3 ! N=3: AQSO4, GASSO4, SOA
  type(production_data), dimension(nregions),public:: AC_diag_prod
  type(production_data), dimension(nregions),public:: AC_O3_lp
  real,dimension(:,:),pointer::temp_prod_so4
  real,dimension(:,:),pointer::temp_prod_vocs
#endif
  logical, public::isoprene_on
  character(len=*), parameter  ::  mname = 'ebischeme'
  !
  ! !REVISION HISTORY: 
  ! 
  !      Feb 2007 - Elisabetta Vignati - changed for inclusion of cloud chemistry on modes
  !   22 Mar 2012 - Philippe Le Sager  - adapted for lon-lat MPI domain decomposition
  !   21 mar 2014 - Jason Williams - port the modified CB05 chemistry scheme
  !
  ! !REMARKS:
  !  (1) initializations are made in chemistry/chemistry_init
  !
  ! !TODO : FIXME check sum_chemistry and sum_deposition
  ! 
  !EOP
  !------------------------------------------------------------------------

contains

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EBI
  !
  ! !DESCRIPTION: perform Eulerian backwards chemistry at one model layer
  !               level in one region 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine EBI( region, level, is, js, rj, rr, y, ye, status) 
    !
    ! !USES:
    !
    use dims,           only : im, jm, nchem, tref, okdebug
    use Dims,           only : CheckShape
    use binas,          only : avog
    use chem_param,     only : xmso4
    use global_data,    only : region_dat
#ifndef without_dry_deposition
    use dry_deposition, only : vd  
#endif
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)        :: region, level, is, js
    real,    intent(in)        :: rr(is:,js:,:)    ! array of reaction rates
    real,    intent(in)        :: rj(is:,js:,:)    ! array of photolysis rates
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    real,    intent(inout)     :: y(is:,js:,:)     ! array of concentration
    real,    intent(inout)     :: ye(is:,js:,:)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)       :: status
    !
    ! !REVISION HISTORY: 
    !   22 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter          ::  rname = mname//'/ebi'
    
#ifdef with_budgets
    real, dimension(:,:,:), allocatable  :: cr2, cr3, cr4   ! reaction budgets
#endif
    real, dimension(:,:,:), allocatable  :: y0

    real    :: dtime, dt, dt2, fnoy
    integer :: iterebi, i, j, ib, maxit
    integer :: io, sfstart, sfend

    !SOA parameters
    real,parameter    :: total_soa_yield_isoprene=0.01
    real,parameter    :: total_soa_yield_terpenes=0.15
    real,parameter    :: y_oh_isop_elvoc=0.0003
    real,parameter    :: y_o3_isop_elvoc=0.0001
    real,parameter    :: y_oh_terp_elvoc=0.01
    real,parameter    :: y_o3_terp_elvoc=0.05

    real    :: y_oh_isop_svoc ! =0.01-0.0003
    real    :: y_o3_isop_svoc ! =0.01-0.0001
    real    :: y_oh_terp_svoc ! =0.15-0.01
    real    :: y_o3_terp_svoc ! =0.15-0.05
 
    ! For vectorization/blocking ....
    ! npts can be varied to optimize cache memory management.
    integer, parameter                 :: npts=15 
    integer, dimension(npts)           :: ipts, jpts
    real,    dimension(npts, nreac )   :: rrv
    real,    dimension(npts, nj    )   :: rjv
    real,    dimension(npts, ntrace)   :: vdv  ! deposition velocities
    real,    dimension(npts)           :: emino
    real,    dimension(npts, maxtrace) :: yv, y0v
    integer                       :: iv, itrace, ivt, n
    integer                       :: i1, j1, i2, j2

    ! --- BEGIN --------------------------------

    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    
    ! check arguments ...
    call CheckShape( (/i2-i1+1, j2-j1+1, nreac   /), shape(rr), status )
    IF_NOTOK_RETURN(status=1)
    call CheckShape( (/i2-i1+1, j2-j1+1, nj      /), shape(rj), status )
    IF_NOTOK_RETURN(status=1)
    call CheckShape( (/i2-i1+1, j2-j1+1, maxtrace/), shape(y ), status )
    IF_NOTOK_RETURN(status=1)
    call CheckShape( (/i2-i1+1, j2-j1+1, n_extra /), shape(ye), status )
    IF_NOTOK_RETURN(status=1)

    dtime=nchem/(2*tref(region))  
    !CMK  iterebi=max(1,nint(dtime/2400))  !needed if nchem <2400
    iterebi=max(1,nint(dtime/1350))  !needed if nchem <2400
    dt=dtime/iterebi

    ! define svoc yield as remainder total soa yield after elvoc yields
    ! total yields at present (to check see variable defs):
    ! 1%  for isoprene
    ! 15% for terpenes

    y_oh_isop_svoc=total_soa_yield_isoprene - y_oh_isop_elvoc
    y_o3_isop_svoc=total_soa_yield_isoprene - y_o3_isop_elvoc
    y_oh_terp_svoc=total_soa_yield_terpenes - y_oh_terp_elvoc
    y_o3_terp_svoc=total_soa_yield_terpenes - y_o3_terp_elvoc


    allocate( y0 (i1:i2, j1:j2, maxtrace))
#ifdef with_budgets
    allocate( cr2(i1:i2, j1:j2, nj      ))
    allocate( cr3(i1:i2, j1:j2, nreac   ))
    allocate( cr4(i1:i2, j1:j2, nreacw  ))
    allocate( temp_prod_so4(npts,2))
    allocate( temp_prod_vocs(npts,8))
#endif

    !*** SCALING OF NOx, which has changed due to transport/deposition
    ! This does not yet work. TODO: Make this working. (FIXME : ask VH what's this is about???)
    do j = j1, j2
    do i = i1, i2
          y(i,j,ino)  =max(1e-3,y(i,j,ino))
          y(i,j,ino2) =max(1e-3,y(i,j,ino2))
          y(i,j,ino3) =max(1e-3,y(i,j,ino3))
          y(i,j,in2o5)=max(1e-3,y(i,j,in2o5))
          y(i,j,ihno4)=max(1e-3,y(i,j,ihno4))
          fnoy=y(i,j,ino)+y(i,j,ino2)+y(i,j,ino3)+2.*y(i,j,in2o5)+y(i,j,ihno4)
          fnoy=y(i,j,inox)/fnoy
          y(i,j,ino)  =fnoy*y(i,j,ino)
          y(i,j,ino2) =fnoy*y(i,j,ino2)
          y(i,j,ino3) =fnoy*y(i,j,ino3)
          y(i,j,in2o5)=fnoy*y(i,j,in2o5)
          y(i,j,ihno4)=fnoy*y(i,j,ihno4)
    end do
    end do
    !
    ! set budget accumulators to zero
    !
#ifdef with_budgets
    cr2 = 0.0
    cr3 = 0.0
    cr4 = 0.0
    temp_prod_so4 = 0.0
    temp_prod_vocs = 0.0
#endif

    !===========================================================
    ! ** Start iterating over CHEMISTRY
    !===========================================================
    do ib=1,iterebi  
       maxit=8   !CMKTEMP
       if(level<=3) maxit = maxit*2   ! lowest layers more iterations

       y0 = y

       !-------------------------------
       ! wet sulphur/ammonia chemistry
       !------------------------------
#ifdef with_budgets
       call wetS(region, level, i1, j1, y0, dt, y, ye, cr4, status)
#else
       call wetS(region, level, i1, j1, y0, dt, y, ye, status)
#endif
       !-------------------------------------
       ! gasphase chemistry using EBI solver
       !-------------------------------------

       y0 = y ! reset initial concentrations

       ! ______do EBI solver_______

       dt2 = dt*dt
       
       ! block the input for EBI for efficiency:
       ! copy values with faster running index in inside loop
       iv = 0
       do j = j1, j2
       do i = i1, i2
             iv = iv+1
             ipts(iv) = i
             jpts(iv) = j

             if(iv==npts) then
                ! copy reaction rates...
                do itrace=1,nreac
                   do ivt=1,npts
                      rrv(ivt,itrace) = rr(ipts(ivt),jpts(ivt),itrace)
                   end do
                end do
                ! copy photolysis rates....
                do itrace=1,nj
                   do ivt=1,npts
                      rjv(ivt,itrace) = rj(ipts(ivt),jpts(ivt),itrace)
                   end do
                end do
                ! copy tracer concentrations ....
                do itrace=1,maxtrace
                   do ivt=1,npts
                      yv(ivt,itrace) = y(ipts(ivt),jpts(ivt),itrace)
                      y0v(ivt,itrace) = y0(ipts(ivt),jpts(ivt),itrace)
                   end do
                end do
                ! deposition ....
#ifndef without_dry_deposition
                if(level == 1) then 
                   do itrace=1,ntrace
                      do ivt=1,npts
                         vdv(ivt,itrace) = &
                              vd(region,itrace)%surf(ipts(ivt),jpts(ivt)) &
                              / ye(ipts(ivt),jpts(ivt),idz)   !1/s 
                      end do
                   end do
                else
                   vdv(:,:) = 0.0
                end if
#endif
                ! copy nox emissions....
                do ivt=1,npts
                   emino(ivt) = ye(ipts(ivt),jpts(ivt),ieno)
                end do

                ! do ebi solver....
                call DO_EBI(iv, npts, maxit, dt, dt2, rrv, rjv, vdv, emino, yv, y0v)

                do itrace=1,maxtrace
                   do ivt=1,npts
                      y(ipts(ivt),jpts(ivt),itrace)=yv(ivt,itrace)
                   end do
                end do
#ifdef with_budgets
#ifdef with_m7
                do ivt=1,npts
                   ! add change in gas-phase so4 production to diagnostic  and change to molec->mass kg
                   diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,1)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,1)+temp_prod_so4(ivt,1)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmso4
                   AC_diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,1)=AC_diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,1)+temp_prod_so4(ivt,1)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmso4

                   AC_diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,3)=AC_diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,3)&
                        +temp_prod_vocs(ivt,1)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc&
                        +temp_prod_vocs(ivt,2)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc&
                        +temp_prod_vocs(ivt,3)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc&
                        +temp_prod_vocs(ivt,4)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc&
                        +temp_prod_vocs(ivt,5)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc&
                        +temp_prod_vocs(ivt,6)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc&
                        +temp_prod_vocs(ivt,7)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc&
                        +temp_prod_vocs(ivt,8)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc




                   diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,5)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,5)+temp_prod_vocs(ivt,1)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc
                   diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,6)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,6)+temp_prod_vocs(ivt,2)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc
                   diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,7)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,7)+temp_prod_vocs(ivt,3)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc
                   diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,8)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,8)+temp_prod_vocs(ivt,4)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc


                   diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,9)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,9)+temp_prod_vocs(ivt,5)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc
                   diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,10)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,10)+temp_prod_vocs(ivt,6)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc
                   diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,11)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,11)+temp_prod_vocs(ivt,7)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc
                   diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,12)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,12)+temp_prod_vocs(ivt,8)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc


                end do
#endif
#endif
                iv=0
             end if
       end do
       end do
    
       ! do the 'remaining' points...
       if(iv > 0) then
          do itrace=1,nreac
             do ivt=1,iv
                rrv(ivt,itrace) = rr(ipts(ivt),jpts(ivt),itrace)
             end do
          end do
          do itrace=1,nj
             do ivt=1,iv
                rjv(ivt,itrace) = rj(ipts(ivt),jpts(ivt),itrace)
             end do
          end do
          do itrace=1,maxtrace
             do ivt=1,iv
                yv(ivt,itrace) = y(ipts(ivt),jpts(ivt),itrace)
                y0v(ivt,itrace) = y0(ipts(ivt),jpts(ivt),itrace)
             end do
          end do
          ! deposition ....
#ifndef without_dry_deposition
          if(level == 1) then 
             do itrace=1,ntrace
                do ivt=1,iv
                   vdv(ivt,itrace) = &
                        vd(region,itrace)%surf(ipts(ivt),jpts(ivt)) &
                        / ye(ipts(ivt),jpts(ivt),idz)   !1/s 
                end do
             end do
          else
             vdv(:,:) = 0.0
          end if
#endif
          do ivt=1,iv
             emino(ivt) = ye(ipts(ivt),jpts(ivt),ieno)
          end do

          !the actual EBI solver on remaining cells 
          call DO_EBI(iv, npts, maxit, dt, dt2, rrv, rjv, vdv, emino, yv, y0v)

          do itrace=1,maxtrace
             do ivt=1,iv
                y(ipts(ivt),jpts(ivt),itrace)=yv(ivt,itrace)
             end do
          end do
#ifdef with_budgets
          do ivt=1,iv
             ! add change in so4 production to diagnostic and change to molec->mass kg
             ! 
             diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,1)= diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,1)+temp_prod_so4(ivt,1)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmso4
             AC_diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,1)=AC_diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,1)+temp_prod_so4(ivt,1)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmso4
             !TOTAL SOA for AERCHEMMIP
             AC_diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,3)=AC_diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,3)&
                  +temp_prod_vocs(ivt,1)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc&
                  +temp_prod_vocs(ivt,2)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc&
                  +temp_prod_vocs(ivt,3)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc&
                  +temp_prod_vocs(ivt,4)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc&
                  +temp_prod_vocs(ivt,5)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc&
                  +temp_prod_vocs(ivt,6)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc&
                  +temp_prod_vocs(ivt,7)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc&
                  +temp_prod_vocs(ivt,8)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc
             
             diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,5)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,5)+temp_prod_vocs(ivt,1)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmso4
             diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,6)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,6)+temp_prod_vocs(ivt,2)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmso4
             diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,7)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,7)+temp_prod_vocs(ivt,3)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmso4
             diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,8)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,8)+temp_prod_vocs(ivt,4)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmso4
             
             
             diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,9)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,9)+temp_prod_vocs(ivt,5)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmso4
             diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,10)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,10)+temp_prod_vocs(ivt,6)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmso4
             diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,11)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,11)+temp_prod_vocs(ivt,7)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmso4
             diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,12)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,12)+temp_prod_vocs(ivt,8)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmso4
             

!*1e6*1e-3/avog*xmso4!yv(ivt,iso4)-y0v(ivt,iso4)
             ! add change in liquid so4 production to diagnostic
             !diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,2)=temp_prod_so4(ivt,2)!yv(ivt,iso4)-y0v(ivt,iso4)
          end do
#endif

       end if

       call NOYmass

       !-------------------------------------
       ! marked tracers
       ! apply after correction of nitrogen components
       !----------------------------------------------

       call MARK_TRAC(region, level, is, js, y, rr, rj, dt, ye)

       !------------------------------------------------------------
       ! increase budget accumulators cr2 and cr3 (cr4 is done in wetS)
       !------------------------------------------------------------

#ifdef with_budgets
       call incc2c3
#endif
       !===========================================================
       ! ** END iterating over CHEMISTRY
       !===========================================================
    end do   !iterebi

#ifdef with_budgets
    call REACBUD
#endif

    deallocate(y0)
#ifdef with_budgets
    deallocate(cr2)
    deallocate(cr3)
    deallocate(cr4)
    deallocate(temp_prod_so4)
    deallocate(temp_prod_vocs)
#endif

    ! ok
    status = 0

  contains 

    subroutine DO_EBI(lvec, npts, maxit, dt, dt2, rrv, rjv, vdv, emino, yv, y0v)

      ! INPUT/OUTPUT
      integer,intent(in) :: lvec, npts, maxit
      real,   intent(in) :: dt, dt2, rrv(npts,nreac), rjv(npts,nj), &
                            vdv(npts,ntrace), emino(npts),          &
                            y0v(npts,maxtrace)
      real,  intent(out) :: yv(npts,maxtrace)

      ! Local
      integer           :: ivec,iter
      real             :: r57,     r89, p1,  r12, r21,   xl1, p2,  xl2, p3, p32,  &
                          xl3,     x1,  x2,  x3,  c1,    c2,  c3,  y2,  xjt,r21t, &
                          r12t,    acub,  bcub,ccub,cubdet,dno2,r56, r65, r75,p5, &
                          xl5,     r66,   x5,  p6,  xl6,   x6,  c6,  x17, y1, c7, &
                          r98,     p8,    xl8, x4,  c5,    xl9, r101,r102,xl7,    &
                          c8 ,     x101,  r1920, r1919, p19, xl19,r2019, xl20,    &
                          xlhno3,  ph2o2, xlh2o2, pch2o, pco, phno3, xlch2o,      &
                          pch3o2,  xlch3o2, pch3o2h, xlch3o2h, pald2,             &
                          xlald2,  pmgly, xlmgly, pole, xleth, xlole,             &
                          xlisop,  prxpar, xlrxpar, ppar, xlpar, pror,            &
                          xlror,   pxo2,xlxo2,pxo2n, xlxo2n, prooh,xlterp,        &
			  xlethooh,xlethp,phcooh,pmcooh,xlc3h6, xlrooh,           &
                          porgntr, xlorgntr, xlco, qdms, pso2, qso2, qso2d,       &
                          qnh3,    qnh2o2, ppnh3, dnh3, pnh2, qnh2, qdms1, qdms2, &
                          pmsa,    pnh3, pispd,xlispd,xlacet,paco2,xlaco2,        &
			  pch3oh,  pc3h7o2,phypro2,xlc3h7o2,xlhypro2,pacet,       &
                          pelvoc,psvoc!RM

      do iter=1,maxit
         do ivec=1,lvec
            ! --- Short living compounds & groups
            ! --- First group: NO NO2 O3
            P1=rjv(ivec,jbno3)*yv(ivec,ino3)+rjv(ivec,jn2o5b)*yv(ivec,in2o5)&
               +rjv(ivec,jhono)*yv(ivec,ihono)+rrv(ivec,knh2o2)*yv(ivec,inh2)+emino(ivec)
            R12=0.
            R21=rrv(ivec,kho2no)*yv(ivec,iho2)+rrv(ivec,kmo2no)*yv(ivec,ich3o2)&
                 +rrv(ivec,kc79)*yv(ivec,ixo2)+rrv(ivec,kc46)*yv(ivec,ic2o3) &
                 +rrv(ivec,kaco2no)*yv(ivec,iaco2)+rrv(ivec,kc3h7o2no)*yv(ivec,ic3h7o2)&
                 +rrv(ivec,khypo2no)*yv(ivec,ihypro2)+rrv(ivec,knh2o2no)*yv(ivec,inh2o2)
            XL1=rrv(ivec,knono3)*yv(ivec,ino3)+rrv(ivec,kc81)*yv(ivec,ixo2n)&
                +rrv(ivec,knh2no)*yv(ivec,inh2)+rrv(ivec,khono)*yv(ivec,ioh)
            XL1 = XL1 + vdv(ivec,ino)
            P2=rjv(ivec,jhno3)*yv(ivec,ihno3)+rjv(ivec,jn2o5a)*yv(ivec,in2o5)&
                 +rrv(ivec,kn2o5)*yv(ivec,in2o5)+rjv(ivec,jano3)*yv(ivec,ino3)&
                 +yv(ivec,ihno4)*(rjv(ivec,jhno4)+rrv(ivec,khno4m)+rrv(ivec,khno4oh)*yv(ivec,ioh))&
                 +rrv(ivec,kohhono)*yv(ivec,ioh)*yv(ivec,ihono)&
                 +2.*rrv(ivec,knono3)*yv(ivec,ino3)*yv(ivec,ino)&
                 +rrv(ivec,kc48)*yv(ivec,ipan)+rrv(ivec,kc59)*yv(ivec,iole)*yv(ivec,ino3)&
                 +rjv(ivec,jorgn)*yv(ivec,iorgntr)&
                 +rjv(ivec,jpana)*yv(ivec,ipan)&
                 +0.2*rrv(ivec,kc78) *yv(ivec,iisop)*yv(ivec,ino3)&
                 +rrv(ivec,kno3mo2)*yv(ivec,ich3o2)*yv(ivec,ino3)&
                 +rrv(ivec,kno3c2o3)*yv(ivec,ic2o3)*yv(ivec,ino3)&
                 +rrv(ivec,kno3xo2)*yv(ivec,ixo2)*yv(ivec,ino3)&
                 +0.47*rrv(ivec,kno3terp)*yv(ivec,ino3)*yv(ivec,iterp)&
                 +rjv(ivec,jmo2no2a)*yv(ivec,ich3o2no2)&
                 +rrv(ivec,kch3o2no2m)*yv(ivec,ich3o2no2)                
            XL2=rrv(ivec,kno2oh)*yv(ivec,ioh)+rrv(ivec,kno2no3)*yv(ivec,ino3)&
                 +rrv(ivec,kno2ho2)*yv(ivec,iho2)+rrv(ivec,kno2o3)*yv(ivec,io3)&
                 +rrv(ivec,kc47)*yv(ivec,ic2o3)+rrv(ivec,knh2no2)*yv(ivec,inh2)&
                 +rrv(ivec,kmo2no2)*yv(ivec,ich3o2)
            XL2 = XL2 + vdv(ivec,ino2)
            P3=rjv(ivec,jano3)*yv(ivec,ino3)+rjv(ivec,jo2) ! JEW : O3P + O2 = O3
            XL3=rrv(ivec,ko3ho2)*yv(ivec,iho2)+rrv(ivec,ko3oh)*yv(ivec,ioh)&
                 +rrv(ivec,kno2o3)*yv(ivec,ino2)+rjv(ivec,jo3d)+rrv(ivec,ko3po3)&
                 +rrv(ivec,kc58)*yv(ivec,iole)&
                 +rrv(ivec,kc62)*yv(ivec,ieth)&
                 +rrv(ivec,kc77)*yv(ivec,iisop)&
                 +rrv(ivec,ko3c3h6)*yv(ivec,ic3h6)&
                 +rrv(ivec,ko3terp)*yv(ivec,iterp)&
                 +rrv(ivec,ko3ispd)*yv(ivec,iispd)&
                 +rrv(ivec,knh2o3)*yv(ivec,inh2)&
                 +rrv(ivec,knh2o2o3)*yv(ivec,inh2o2)
            XL3 = XL3 + vdv(ivec,io3)

            X1=y0v(ivec,ino)+P1*DT
            X2=y0v(ivec,ino2)+P2*DT
            X3=y0v(ivec,io3)+P3*DT
            C1=1.+XL1*DT
            C2=1.+XL2*DT
            C3=1.+XL3*DT
            Y1=rrv(ivec,knoo3)*DT
            R21T=R21*DT
            R12T=R12*DT
            XJT=rjv(ivec,jno2)*DT
            P32=0.4*rrv(ivec,kc50)*yv(ivec,ic2o3)*yv(ivec,iho2)
            ! --- Solve unknown: x
            ACUB=-2.*Y1*(C2+R12T+C2*R21T/C1)
            BCUB=2.*C1*C2*C3+2.*C1*C3*(R12T+XJT)+2.*C2*C3*R21T+&
                 2.*Y1*(R12T*(X1-X2)+2.*C2*R21T*X1/C1+C2*(X1+X3))
            CCUB=2.*C1*C3*X2*(R12T+XJT)-2.*C2*C3*X1*R21T+2.*Y1*X1*&
                 (X2*R12T-C2*X3-C2*R21T*X1/C1)
            CUBDET=BCUB*BCUB-4.*ACUB*CCUB
            DNO2=(-1.*BCUB+SQRT(CUBDET))/(2.*ACUB)
            dno2=min(x1,dno2)
            yv(ivec,ino2)=(X2+DNO2)/C2
            yv(ivec,ino)=(X1-DNO2)/C1
            yv(ivec,io3)=(X3+(P32*dt)+XJT*yv(ivec,ino2))/(C3+Y1*yv(ivec,ino))
            
            ! --- Second group: yv(ivec,iho2) yv(ivec,ioh) yv(ivec,ihno4) yv(ivec,ihono)
            R57=rjv(ivec,jhno4)+rrv(ivec,khno4m)
            R56=rrv(ivec,kcooh)*yv(ivec,ico)+rrv(ivec,ko3oh)*yv(ivec,io3)+rrv(ivec,khpoh)&
                 *yv(ivec,ih2o2)+rrv(ivec,kfrmoh)*yv(ivec,ich2o)+rrv(ivec,kh2oh)&
                 +rrv(ivec,kso2oh)*yv(ivec,iso2)

            p5=2.*rjv(ivec,jbch2o)*yv(ivec,ich2o)&
                 +rrv(ivec,kmo2no)*yv(ivec,ich3o2)*yv(ivec,ino)&
                 +rjv(ivec,jmepe)*yv(ivec,ich3o2h)&
                 +2.0*rjv(ivec,jald2)*yv(ivec,iald2)&
                 +rjv(ivec,jmgly)*yv(ivec,imgly)+0.11*rrv(ivec,kc52)*yv(ivec,ipar)*yv(ivec,ioh)&
                 +0.94*rrv(ivec,kc53)*yv(ivec,iror)+rrv(ivec,kc54)*yv(ivec,iror)&
                 +1.57*rrv(ivec,kc57)*yv(ivec,iole)*yv(ivec,ioh)&
                 +0.76*rrv(ivec,kc58)*yv(ivec,io3)*yv(ivec,iole)&
                 +0.56*rrv(ivec,kc59)*yv(ivec,ino3)*yv(ivec,iole)&
                 +rrv(ivec,kc61)*yv(ivec,ieth)*yv(ivec,ioh)&
                 +0.22*rrv(ivec,kc62)*yv(ivec,ieth)*yv(ivec,io3)&
                 +0.066*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)&
                 +rrv(ivec,kc41)*yv(ivec,ich2o)*yv(ivec,ino3)&
                 +0.9*rrv(ivec,kc84)*yv(ivec,ioh)*yv(ivec,iorgntr)&
                 +0.9*rjv(ivec,jorgn)*yv(ivec,iorgntr)&
                 +0.8*rrv(ivec,kc78)*yv(ivec,iisop)*yv(ivec,ino3)&
                 +0.74*rrv(ivec,kmo2mo2)*yv(ivec,ich3o2)*yv(ivec,ich3o2)&
                 +0.9*rjv(ivec,jrooh)*yv(ivec,irooh)&
                 +rrv(ivec,kno3mo2)*yv(ivec,ich3o2)*yv(ivec,ino3)&
                 +0.19*rrv(ivec,ko3c3h6)*yv(ivec,io3)*yv(ivec,ic3h6)&
                 +0.28*rrv(ivec,ko3terp)*yv(ivec,io3)*yv(ivec,iterp)&
                 +0.75*rrv(ivec,kno3terp)*yv(ivec,ino3)*yv(ivec,iterp)&
                 +0.154*rrv(ivec,ko3ispd)*yv(ivec,io3)*yv(ivec,iispd)&
                 +0.925*rrv(ivec,kno3ispd)*yv(ivec,ino3)*yv(ivec,iispd)&
                 +1.033*rjv(ivec,jispd)*yv(ivec,iispd)&
                 +rrv(ivec,kohch3oh)*yv(ivec,ich3oh)&
                 +rrv(ivec,kohhcooh)*yv(ivec,ihcooh)&
                 +rrv(ivec,kohethoh)*yv(ivec,iethoh)&
                 +1.22*rrv(ivec,kohterp)*yv(ivec,iterp)&
                 +rrv(ivec,kohc2h6)*yv(ivec,ic2h6)&
                 +0.5*rrv(ivec,kc76)*yv(ivec,ioh)*yv(ivec,iisop)& ! reduced according to archibald et al, AE, 2011
                 +0.503*rrv(ivec,kohispd)*yv(ivec,ioh)*yv(ivec,iispd)&
                 +0.5*rrv(ivec,kaco2mo2)*yv(ivec,iaco2)*yv(ivec,ich3o2)&
                 +rrv(ivec,kaco2no)*yv(ivec,iaco2)*yv(ivec,ino)&
                 +rrv(ivec,kc3h7o2no)*yv(ivec,ic3h7o2)*yv(ivec,ino)&
                 +rrv(ivec,kc3h7o2ho2)*yv(ivec,ihypro2)*yv(ivec,ino)&
                 +rjv(ivec,jmo2no2b)*yv(ivec,ich3o2no2)      
            XL5=rrv(ivec,kho2no)    *yv(ivec,ino)&
                 +rrv(ivec,kno2ho2) *yv(ivec,ino2)&
                 +rrv(ivec,ko3ho2)  *yv(ivec,io3)&
                 +rrv(ivec,kmo2ho2a) *yv(ivec,ich3o2)&
                 +rrv(ivec,kmo2ho2b) *yv(ivec,ich3o2)&
                 +rrv(ivec,kc50)    *yv(ivec,ic2o3)&
                 +rrv(ivec,kho2oh)  *yv(ivec,ioh)&
                 +rrv(ivec,kc82)    *yv(ivec,ixo2)&
                 +rrv(ivec,kc85)    *yv(ivec,ixo2n)&
                 +rrv(ivec,kno3ho2) *yv(ivec,ino3)&
                 +rrv(ivec,kaco2ho2)*yv(ivec,iaco2)&
                 +rrv(ivec,knh2ho2)*yv(ivec,inh2)&
                 +rrv(ivec,knh2o2ho2)*yv(ivec,inh2o2)&
                 +rrv(ivec,kho2_aer) &
                 +rrv(ivec,kho2_l)               
            R66=2.*rrv(ivec,kho2ho2)
            X5=y0v(ivec,iho2)+P5*DT
            R65=rrv(ivec,kho2no)*yv(ivec,ino)+rrv(ivec,ko3ho2)*yv(ivec,io3)

            P6=rjv(ivec,jhno3)*yv(ivec,ihno3)&
               +2.*rjv(ivec,jo3d)*yv(ivec,io3)&
               +2.*rjv(ivec,jh2o2)*yv(ivec,ih2o2)&
               +rjv(ivec,jmepe)*yv(ivec,ich3o2h)&
               +0.1*rrv(ivec,kc58)*yv(ivec,io3)*yv(ivec,iole)&
               +0.266*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)&
               +rjv(ivec,jrooh)*yv(ivec,irooh)&
               +0.12*rrv(ivec,kc62)*yv(ivec,ieth)*yv(ivec,io3)&
               +0.33*rrv(ivec,ko3c3h6)*yv(ivec,io3)*yv(ivec,ic3h6)&
               +0.57*rrv(ivec,ko3terp)*yv(ivec,iterp)*yv(ivec,io3)&
               +0.268*rrv(ivec,ko3ispd)*yv(ivec,io3)*yv(ivec,iispd)
            XL6=rrv(ivec,khno4oh)*yv(ivec,ihno4)&
               +rrv(ivec,kohhono)*yv(ivec,ihono)& 
               +rrv(ivec,kho2oh)*yv(ivec,iho2)&
               +rrv(ivec,kno2oh)*yv(ivec,ino2)&
               +rrv(ivec,khono)*yv(ivec,ino)&
               +rrv(ivec,kohhno3)*yv(ivec,ihno3)&
               +rrv(ivec,kcooh)*yv(ivec,ico)&
               +rrv(ivec,ko3oh)*yv(ivec,io3)&
               +rrv(ivec,khpoh)*yv(ivec,ih2o2)&
               +rrv(ivec,kfrmoh)*yv(ivec,ich2o)&
               +rrv(ivec,kch4oh)*yv(ivec,ich4)&
               +0.7*rrv(ivec,kohmper)*yv(ivec,ich3o2h)&
               +rrv(ivec,kc43)*yv(ivec,iald2)&
               +rrv(ivec,kc73)*yv(ivec,imgly)&
               +rrv(ivec,kc52)*yv(ivec,ipar)&
               +rrv(ivec,kc57)*yv(ivec,iole)&
               +rrv(ivec,kc61)*yv(ivec,ieth)&
               +rrv(ivec,kc76)*yv(ivec,iisop)&
               +0.77*rrv(ivec,kohrooh)*yv(ivec,irooh)& ! note: change from '0.7' to '0.77'
               +rrv(ivec,kc84)*yv(ivec,iorgntr)&
               +rrv(ivec,kh2oh)&
               +rrv(ivec,kso2oh)*yv(ivec,iso2) &  ! bug found by Jason 01/2008
               +(rrv(ivec,kdmsoha)+rrv(ivec,kdmsohb)) *yv(ivec,idms) &
               +rrv(ivec,knh3oh)*yv(ivec,inh3)&
               +rrv(ivec,knh2oh)*yv(ivec,inh2)&
               +rrv(ivec,kohch3oh)*yv(ivec,ich3oh)&
               +rrv(ivec,kohhcooh)*yv(ivec,ihcooh)&
               +rrv(ivec,kohethoh)*yv(ivec,iethoh)&
               +rrv(ivec,kohterp)*yv(ivec,iterp)&
               +rrv(ivec,kohispd)*yv(ivec,iispd)&
               +rrv(ivec,kohmcooh)*yv(ivec,imcooh)&
               +rrv(ivec,kohc2h6)*yv(ivec,ic2h6)&
               +rrv(ivec,kohc3h8)*yv(ivec,ic3h8)&
               +rrv(ivec,kohc3h6)*yv(ivec,ic3h6)&
               +rrv(ivec,kohacet)*yv(ivec,iacet)
            R101=rjv(ivec,jhono)               
            X6=y0v(ivec,ioh)+(P6*DT)+(R101*y0v(ivec,ihono)*DT)
            C6=1.+XL6*DT
            R75=rrv(ivec,kno2ho2)*yv(ivec,ino2)     
            R102=rrv(ivec,khono)*yv(ivec,ino)
            X101=rjv(ivec,jhono)+rrv(ivec,kohhono)*yv(ivec,ioh)
            X101= X101 + vdv(ivec,ihono)
            C8=1.+X101*DT
            XL7=rjv(ivec,jhno4)+rrv(ivec,khno4oh)*yv(ivec,ioh)+rrv(ivec,khno4m)
            XL7 = XL7 + vdv(ivec,ihno4)
            C7=1.+XL7*DT
            Y1=R57/C7
            Y2=R56/C6
            ACUB=R66*DT
            BCUB=1.+XL5*DT-DT2*(Y1*R75+Y2*R65)
            CCUB=-1.*X5-DT*(Y1*y0v(ivec,ihno4)+Y2*X6)
            CUBDET=BCUB*BCUB-4.*ACUB*CCUB
            CUBDET=MAX(CUBDET,1.E-20)
            yv(ivec,iho2)=max(0.1,(-1.*BCUB+SQRT(CUBDET))/(2.*ACUB))
            yv(ivec,ioh)=(X6+R65*yv(ivec,iho2)*DT)/C6
            yv(ivec,ihono)=(y0v(ivec,ihono)+R102*DT*yv(ivec,ioh))/C8
            yv(ivec,ihno4)=(y0v(ivec,ihno4)+R75*DT*yv(ivec,iho2))/C7
            !
            ! --- Third group: NO3 N2O5
            !
            R89=rjv(ivec,jn2o5a)+rjv(ivec,jn2o5b)+rrv(ivec,kn2o5)
            P8=rrv(ivec,kohhno3)*yv(ivec,ihno3)*yv(ivec,ioh)+rrv(ivec,kno2o3)*yv(ivec,ino2)*yv(ivec,io3)&
               +rjv(ivec,jmo2no2b)*yv(ivec,ich3o2no2)+rjv(ivec,jpanb)*yv(ivec,ipan) 
            XL8=rjv(ivec,jbno3)+rjv(ivec,jano3)&
                 +rrv(ivec,knono3)*yv(ivec,ino)&
                 +rrv(ivec,kno2no3)*yv(ivec,ino2)&
                 +rrv(ivec,kc44)*yv(ivec,iald2)&
                 +rrv(ivec,kc59)*yv(ivec,iole)&
                 +rrv(ivec,kc78)*yv(ivec,iisop)&
                 +rrv(ivec,kc41)*yv(ivec,ich2o)&
                 +rrv(ivec,kdmsno3)*yv(ivec,idms)&
                 +rrv(ivec,kno3ho2)*yv(ivec,iho2)&
                 +rrv(ivec,kno3mo2)*yv(ivec,ich3o2)&
                 +rrv(ivec,kno3c2o3)*yv(ivec,ic2o3)&
                 +rrv(ivec,kno3xo2)*yv(ivec,ixo2)&
                 +rrv(ivec,kno3c3h6)*yv(ivec,ic3h6)&
                 +rrv(ivec,kno3terp)*yv(ivec,iterp)&
                 +rrv(ivec,kno3ispd)*yv(ivec,iispd)
            XL8 = XL8 + vdv(ivec,ino3) 
            X4=y0v(ivec,ino3)+P8*DT
            C5=1.+XL8*DT
            R98=rrv(ivec,kno2no3)*yv(ivec,ino2)
            XL9=rjv(ivec,jn2o5a)+rjv(ivec,jn2o5b)+rrv(ivec,kn2o5)+rrv(ivec,kn2o5_aer)+rrv(ivec,kn2o5l)   !cmk rates now idependent from y
            XL9 = XL9 + vdv(ivec,in2o5)
            C6=1.+XL9*DT
            C7=(C5*C6-R89*R98*DT2)
            yv(ivec,in2o5)=(C5*y0v(ivec,in2o5)+R98*DT*X4)/C7
            yv(ivec,ino3)=(C6*X4+R89*DT*y0v(ivec,in2o5))/C7
            !
            ! --- Fourth group: C2O3 PAN
            !
            R1920=rrv(ivec,kc48)+rjv(ivec,jpana)
            R1919=rrv(ivec,kc49)
            P19=rrv(ivec,kc43)*yv(ivec,iald2)*yv(ivec,ioh)&
                 +rrv(ivec,kc44)*yv(ivec,iald2)*yv(ivec,ino3)&
                 +rrv(ivec,kc73)*yv(ivec,imgly)*yv(ivec,ioh) &
                 +rjv(ivec,jmgly)*yv(ivec,imgly)&
                 +0.2*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)&
                 +0.74*rjv(ivec,jrooh)*yv(ivec,irooh)&
                 +0.74*rrv(ivec,kc84)*yv(ivec,ioh)*yv(ivec,iorgntr)&
                 +0.74*rjv(ivec,jorgn)*yv(ivec,iorgntr)&
                 +0.39*rrv(ivec,ko3terp)*yv(ivec,io3)*yv(ivec,iterp)&
                 +0.498*rrv(ivec,kohispd)*yv(ivec,ioh)*yv(ivec,iispd)&
                 +0.114*rrv(ivec,ko3ispd)*yv(ivec,io3)*yv(ivec,iispd)&
                 +0.075*rrv(ivec,kno3ispd)*yv(ivec,ino3)*yv(ivec,iispd)&
                 +0.967*rjv(ivec,jispd)*yv(ivec,iispd)&
                 +0.3*rrv(ivec,kaco2mo2)*yv(ivec,iaco2)*yv(ivec,ich3o2)&
                 +rrv(ivec,kaco2no)*yv(ivec,iaco2)*yv(ivec,ino)&
                 +rjv(ivec,jb_acet)*yv(ivec,iacet)&
                 +rjv(ivec,jmo2no2b)*yv(ivec,ich3o2no2)      
            XL19=rrv(ivec,kc46)*yv(ivec,ino)+rrv(ivec,kc50)*yv(ivec,iho2)&
                 +rrv(ivec,kc47)*yv(ivec,ino2)&
                 +rrv(ivec,kno3c2o3)*yv(ivec,ino3)
            XL19 = XL19 + vdv(ivec,ic2o3)
            R2019=rrv(ivec,kc47)*yv(ivec,ino2)
            XL20=rrv(ivec,kc48)+rjv(ivec,jpana)+rjv(ivec,jpanb)
            XL20 = XL20 + vdv(ivec,ipan)
            ACUB=2*R1919*DT*(1+XL20*DT)
            BCUB=(1+XL20*DT)*(1+XL19*DT)-R1920*DT*R2019*DT
            CCUB=(1+XL20*DT)*(y0v(ivec,ic2o3)+P19*DT)+R1920*DT*y0v(ivec,ipan)
            CUBDET=BCUB*BCUB+4.*ACUB*CCUB
            yv(ivec,ic2o3)=max(1e-8,(-1.*BCUB+SQRT(CUBDET))/(2.*ACUB))    !cmk  put max here....
            yv(ivec,ipan)=(y0v(ivec,ipan)+R2019*yv(ivec,ic2o3)*DT)/(1.+XL20*DT)

            !
            ! ---- Fifth group : CH3O2 CH3O2NO2
            R1920=rrv(ivec,kch3o2no2m)+rjv(ivec,jmo2no2a)+rjv(ivec,jmo2no2b)
            R1919=rrv(ivec,kmo2mo2)
            P19=rrv(ivec,kch4oh)*yv(ivec,ich4)*yv(ivec,ioh) &
                 +0.7*rrv(ivec,kohmper)*yv(ivec,ioh)*yv(ivec,ich3o2h)&
                 +rrv(ivec,kno3c2o3)*yv(ivec,ino3)*yv(ivec,ic2o3)&
                 +rrv(ivec,kc46)*yv(ivec,ino)*yv(ivec,ic2o3)&
                 +2*rrv(ivec,kc49)*yv(ivec,ic2o3)*yv(ivec,ic2o3)&
                 +rjv(ivec,jpanb)*yv(ivec,ipan)&
                 +rjv(ivec,jald2)*yv(ivec,iald2)&
                 +0.74*rjv(ivec,jrooh)*yv(ivec,irooh)&
                 +0.74*rrv(ivec,kc84)*yv(ivec,iorgntr)&
                 +0.74*rjv(ivec,jorgn)*yv(ivec,iorgntr)&
                 +rrv(ivec,kohmcooh)*yv(ivec,ioh)*yv(ivec,imcooh)&
                 +0.31*rrv(ivec,ko3c3h6)*yv(ivec,io3)*yv(ivec,ic3h6)&
                 +0.39*rrv(ivec,ko3terp)*yv(ivec,iterp)*yv(ivec,io3)&
                 +2.0*rjv(ivec,ja_acet)*yv(ivec,iacet)+rjv(ivec,jb_acet)*yv(ivec,iacet)         
            XL19=rrv(ivec,kmo2no)*yv(ivec,ino)+rrv(ivec,kmo2ho2a)*yv(ivec,iho2)&
                 +rrv(ivec,kmo2ho2b)*yv(ivec,iho2)+rrv(ivec,kno3mo2)*yv(ivec,ino3)&
                 +rrv(ivec,kaco2mo2)*yv(ivec,iaco2)+rrv(ivec,kmo2no2)*yv(ivec,ino2)
            R2019=rrv(ivec,kmo2no2)*yv(ivec,ino2)
            XL20=rrv(ivec,kch3o2no2m)+rjv(ivec,jmo2no2a)+rjv(ivec,jmo2no2b)
            XL20 = XL20 + vdv(ivec,ich3o2no2)
            ACUB=2*R1919*DT*(1+XL20*DT)
            BCUB=(1+XL20*DT)*(1+XL19*DT)-R1920*DT*R2019*DT
            CCUB=(1+XL20*DT)*(y0v(ivec,ich3o2)+P19*DT)+R1920*DT*y0v(ivec,ich3o2no2)
            CUBDET=BCUB*BCUB+4.*ACUB*CCUB
            yv(ivec,ich3o2)=max(1e-8,(-1.*BCUB+SQRT(CUBDET))/(2.*ACUB))    !cmk  put max here....
            yv(ivec,ich3o2no2)=(y0v(ivec,ich3o2no2)+R2019*yv(ivec,ich3o2)*DT)/(1.+XL20*DT)      

            !
            ! -------- ISPD chemistry 
            !
            pispd=0.912*rrv(ivec,kc76)*yv(ivec,ioh)*yv(ivec,iisop)&
                +0.65*rrv(ivec,kc77)*yv(ivec,io3)*yv(ivec,iisop)&
                +0.2*rrv(ivec,kc78)*yv(ivec,ino3)*yv(ivec,iisop)
            xlispd=rrv(ivec,kohispd)*yv(ivec,ioh)+rrv(ivec,ko3ispd)*yv(ivec,io3)&
               +rrv(ivec,kno3ispd)*yv(ivec,ino3)+rjv(ivec,jispd)+vdv(ivec,iispd)

            yv(ivec,iispd)=(y0v(ivec,iispd)+pispd*DT)/(1.+xlispd*DT)
            !
            ! -------- ACO2 chemistry 
            !
            paco2=rrv(ivec,kohacet)*yv(ivec,iacet)*yv(ivec,ioh)
            
            xlaco2=rrv(ivec,kaco2ho2)*yv(ivec,iho2)+rrv(ivec,kaco2mo2)*yv(ivec,ich3o2)&
                   +rrv(ivec,kaco2no)*yv(ivec,ino)+rrv(ivec,kaco2xo2)*yv(ivec,ixo2)
            
            yv(ivec,iaco2)=(y0v(ivec,iaco2)+paco2*DT)/(1.+xlaco2*DT)
            
            !
            ! -------- C3H7O2 chemistry 
            !
            pc3h7o2=rrv(ivec,kohc3h8)*yv(ivec,ic3h8)*yv(ivec,ioh)
            
            xlc3h7o2=rrv(ivec,kc3h7o2no)*yv(ivec,ino)+rrv(ivec,kc3h7o2ho2)*yv(ivec,iho2)
            
            yv(ivec,ic3h7o2)=(y0v(ivec,ic3h7o2)+pc3h7o2*DT)/(1.+xlc3h7o2*DT)        
            
            !
            ! -------- HYPRO2 chemistry 
            !       
            
            phypro2=rrv(ivec,kohc3h6)*yv(ivec,ic3h6)*yv(ivec,ioh)
            
            xlhypro2=rrv(ivec,khypo2no)*yv(ivec,ino)+rrv(ivec,khypo2ho2)*yv(ivec,iho2)
            
            yv(ivec,ihypro2)=(y0v(ivec,ihypro2)+phypro2*DT)/(1.+xlhypro2*DT)                
            
            !       
            !
            ! --- CBM4 chem.(short living compounds & operators)
            !
            PRXPAR=0.11*rrv(ivec,kc52)*yv(ivec,ioh)*yv(ivec,ipar)&
               +2.1*rrv(ivec,kc53)*yv(ivec,iror)&
               +0.7*rrv(ivec,kc57)*yv(ivec,iole)*yv(ivec,ioh)&
               +rrv(ivec,kc58)*yv(ivec,io3)*yv(ivec,iole)&
               +rrv(ivec,kc59)*yv(ivec,iole)*yv(ivec,ino3)&
               +rrv(ivec,kohrooh)*yv(ivec,ioh)*yv(ivec,irooh)&
               +1.98*rjv(ivec,jrooh)*yv(ivec,irooh)&
               +1.98*rrv(ivec,kc84)*yv(ivec,ioh)*yv(ivec,iorgntr)&
               +1.98*rjv(ivec,jorgn)*yv(ivec,iorgntr)     


            XLRXPAR=rrv(ivec,kc83)*yv(ivec,ipar)
            yv(ivec,irxpar)=(y0v(ivec,irxpar)+PRXPAR*DT)/(1.+XLRXPAR*DT)
            XLISOP=rrv(ivec,kc76)*yv(ivec,ioh)+rrv(ivec,kc77)*yv(ivec,io3)+rrv(ivec,kc78)*yv(ivec,ino3)
            yv(ivec,iisop)=y0v(ivec,iisop)/(1.+XLISOP*DT)
            PROR=0.76*rrv(ivec,kc52)*yv(ivec,ipar)*yv(ivec,ioh)+0.02*rrv(ivec,kc53)*yv(ivec,iror)
            XLROR=rrv(ivec,kc53)+rrv(ivec,kc54)
            yv(ivec,iror)=(y0v(ivec,iror)+PROR*DT)/(1.+XLROR*DT)

            xlterp=rrv(ivec,kohterp)*yv(ivec,ioh)+rrv(ivec,ko3terp)*yv(ivec,io3)&
              +rrv(ivec,kno3terp)*yv(ivec,ino3)     
            yv(ivec,iterp)=y0v(ivec,iterp)/(1.+xlterp*dt)

#ifdef with_m7
            !ELVOC
            ! contribution from monoterpene
            pelvoc=y_oh_terp_elvoc*rrv(ivec,kohterp)*yv(ivec,ioh)*y0v(ivec,iterp)+y_o3_terp_elvoc*rrv(ivec,ko3terp)*yv(ivec,io3)*y0v(ivec,iterp)

#ifdef with_budgets

            temp_prod_vocs(ivec,1)=y_oh_terp_elvoc * rrv(ivec,kohterp)*yv(ivec,ioh)*y0v(ivec,iterp)*DT
            temp_prod_vocs(ivec,2)=y_o3_terp_elvoc * rrv(ivec,ko3terp)*yv(ivec,io3)*y0v(ivec,iterp)*DT

#endif

            !Contribution from isoprene

            if (isoprene_on)then
               pelvoc=pelvoc+ y_oh_isop_elvoc*rrv(ivec,kc76)*yv(ivec,ioh)*y0v(ivec,iisop)+y_o3_isop_elvoc*rrv(ivec,kc77)*yv(ivec,io3)*y0v(ivec,iisop)
#ifdef with_budgets
               temp_prod_vocs(ivec,3) = y_oh_isop_elvoc*rrv(ivec,kc76)*yv(ivec,ioh)*y0v(ivec,iisop)*DT
               temp_prod_vocs(ivec,4) = y_o3_isop_elvoc*rrv(ivec,kc77)*yv(ivec,io3)*y0v(ivec,iisop)*DT

#endif           
            endif


            yv(ivec,ielvoc)=y0v(ivec,ielvoc)+pelvoc*dt
            !SVOC
            ! contribution from monoterpene
            psvoc=y_oh_terp_svoc*rrv(ivec,kohterp)*yv(ivec,ioh)*y0v(ivec,iterp)+0.10*rrv(ivec,ko3terp)*yv(ivec,io3)*y0v(ivec,iterp)

#ifdef with_budgets

            temp_prod_vocs(ivec,5) = y_oh_terp_svoc*rrv(ivec,kohterp)*yv(ivec,ioh)*y0v(ivec,iterp)*DT
            temp_prod_vocs(ivec,6) = y_o3_terp_svoc*rrv(ivec,ko3terp)*yv(ivec,io3)*y0v(ivec,iterp)*DT
#endif
            !Contribution from isoprene
            if (isoprene_on)then
               psvoc=psvoc + y_oh_isop_svoc*rrv(ivec,kc76)*yv(ivec,ioh)*y0v(ivec,iisop)+y_o3_isop_svoc*rrv(ivec,kc77)*yv(ivec,io3)*y0v(ivec,iisop)

#ifdef with_budgets

               temp_prod_vocs(ivec,7)= y_oh_isop_svoc*rrv(ivec,kc76)*yv(ivec,ioh)*y0v(ivec,iisop)*DT
               temp_prod_vocs(ivec,8)= y_o3_isop_svoc*rrv(ivec,kc77)*yv(ivec,io3)*y0v(ivec,iisop)*DT
#endif
            end if
            yv(ivec,isvoc)=y0v(ivec,isvoc)+psvoc*dt

#endif

            pxo2=rrv(ivec,kc73)*yv(ivec,imgly)*yv(ivec,ioh)&
                 +0.87*rrv(ivec,kc52)*yv(ivec,ipar)*yv(ivec,ioh)&
                 +0.96*rrv(ivec,kc53)*yv(ivec,iror)&
                 +0.8*rrv(ivec,kc57)*yv(ivec,iole)*yv(ivec,ioh)&
                 +0.22*rrv(ivec,kc58)*yv(ivec,io3)*yv(ivec,iole)&
                 +0.91*rrv(ivec,kc59)*yv(ivec,iole)*yv(ivec,ino3)&
                 +rrv(ivec,kc61)*yv(ivec,ieth)*yv(ivec,ioh)&
                 +0.991*rrv(ivec,kc76)*yv(ivec,iisop)*yv(ivec,ioh)&
                 +rrv(ivec,kc78)*yv(ivec,iisop)*yv(ivec,ino3)&
                 +0.77*rrv(ivec,kohrooh)*yv(ivec,irooh)*yv(ivec,ioh)&
                 +0.5*rjv(ivec,jrooh)*yv(ivec,irooh)&
                 +0.51*rrv(ivec,kc84)*yv(ivec,ioh)*yv(ivec,iorgntr)&         
                 +0.51*rjv(ivec,jorgn)*yv(ivec,iorgntr)&
                 +0.2*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)&
                 +0.1*rrv(ivec,kohethoh)*yv(ivec,ioh)*yv(ivec,iethoh)&
                 +0.991*rrv(ivec,kohc2h6)*yv(ivec,ioh)*yv(ivec,ic2h6)&
                 +1.25*rrv(ivec,kohterp)*yv(ivec,ioh)*yv(ivec,iterp)&
                 +0.76*rrv(ivec,ko3terp)*yv(ivec,io3)*yv(ivec,iterp)&
                 +1.03*rrv(ivec,kno3terp)*yv(ivec,ino3)*yv(ivec,iterp)&
                 +0.713*rrv(ivec,kohispd)*yv(ivec,ioh)*yv(ivec,iispd)&
                 +0.064*rrv(ivec,ko3ispd)*yv(ivec,io3)*yv(ivec,iispd)&
                 +0.075*rrv(ivec,kno3ispd)*yv(ivec,ino3)*yv(ivec,iispd)&
                 +0.7*rjv(ivec,jispd)*yv(ivec,iispd)              
                             
            xlxo2=rrv(ivec,kc79)*yv(ivec,ino)+2.*rrv(ivec,kc80)*yv(ivec,ixo2)&
                 +rrv(ivec,kc82)*yv(ivec,iho2)+rrv(ivec,kno3xo2)*yv(ivec,ino3)&
                 +rrv(ivec,kaco2xo2)*yv(ivec,iaco2)+rrv(ivec,kxo2xo2n)*yv(ivec,ixo2n)
            yv(ivec,ixo2)=(y0v(ivec,ixo2)+pxo2*dt)/(1.+xlxo2*dt)

            pxo2n=0.13*rrv(ivec,kc52)*yv(ivec,ipar)*yv(ivec,ioh)&
                 +0.04*rrv(ivec,kc53)*yv(ivec,iror)&
                 +0.09*rrv(ivec,kc59)*yv(ivec,iole)*yv(ivec,ino3)&
                 +0.009*rrv(ivec,kohc2h6)*yv(ivec,ioh)*yv(ivec,ic2h6)&
                 +0.088*rrv(ivec,kc76)*yv(ivec,iisop)*yv(ivec,ioh)&
                 +0.25*rrv(ivec,kohterp)*yv(ivec,ioh)*yv(ivec,iterp)&
                 +0.18*rrv(ivec,ko3terp)*yv(ivec,io3)*yv(ivec,iterp)&
                 +0.25*rrv(ivec,kno3terp)*yv(ivec,ino3)*yv(ivec,iterp)
            xlxo2n=rrv(ivec,kc81)*yv(ivec,ino)+rrv(ivec,kc85)*yv(ivec,iho2)&
                   +rrv(ivec,kxo2xo2n)*yv(ivec,ixo2n)+rrv(ivec,kxo2n)*yv(ivec,ixo2n)
            yv(ivec,ixo2n)=(y0v(ivec,ixo2n)+pxo2n*dt)/(1.+xlxo2n*dt)
          
         end do !ivec
         if ( mod(iter,2) == 0 ) then

            do ivec=1,lvec
               ! --- Species with intermediate lifetimes
               ! --- Inorganic compounds (HNO3 H2O2)
               ! 
               PHNO3=rrv(ivec,kno2oh)*yv(ivec,ino2)*yv(ivec,ioh)&
                    +2.*(rrv(ivec,kn2o5_aer)+rrv(ivec,kn2o5l))*yv(ivec,in2o5)&
                    +rrv(ivec,kc44)*yv(ivec,iald2)*yv(ivec,ino3)&
                    +rrv(ivec,kc41)*yv(ivec,ich2o)*yv(ivec,ino3)&
                    +rrv(ivec,kc84)*yv(ivec,ioh)*yv(ivec,iorgntr)&
                    +rrv(ivec,kno3ho2)*yv(ivec,ino3)*yv(ivec,iho2)&
                    +0.15*rrv(ivec,kno3ispd)*yv(ivec,ino3)*yv(ivec,iispd)&
                    +rrv(ivec,kdmsno3)*yv(ivec,ino3)*yv(ivec,idms)&
                    +rrv(ivec,kno3_aer)*yv(ivec,ino3)
               XLHNO3=rjv(ivec,jhno3)+rrv(ivec,kohhno3)*yv(ivec,ioh)
               XLHNO3=XLHNO3 + vdv(ivec,ihno3)
               yv(ivec,ihno3)=(y0v(ivec,ihno3)+PHNO3*DT)/(1.+XLHNO3*DT)
               PH2O2=rrv(ivec,kho2ho2)*yv(ivec,iho2)*yv(ivec,iho2)&
                     +0.5*rrv(ivec,kho2_aer)*yv(ivec,iho2)
               XLH2O2=rjv(ivec,jh2o2)+rrv(ivec,khpoh)*yv(ivec,ioh)
               XLH2O2=XLH2O2 + vdv(ivec,ih2o2)
               yv(ivec,ih2o2)=(y0v(ivec,ih2o2)+PH2O2*DT)/(1.+XLH2O2*DT)
               ! --- CH4-chemistry (methyl peroxide formaldehyde)
               PCH3O2H=rrv(ivec,kmo2ho2a)*yv(ivec,ich3o2)*yv(ivec,iho2)
               XLCH3O2H=rrv(ivec,kohmper)*yv(ivec,ioh)+rjv(ivec,jmepe)
               XLCH3O2H=XLCH3O2H + vdv(ivec,ich3o2h)
               yv(ivec,ich3o2h)=(y0v(ivec,ich3o2h)+PCH3O2H*DT)/(1.+XLCH3O2H*DT)
               pch2o=0.3*rrv(ivec,kohmper)*yv(ivec,ich3o2h)*yv(ivec,ioh)&
                  +rrv(ivec,kmo2no)*yv(ivec,ich3o2)*yv(ivec,ino)&
                  +rrv(ivec,kmo2ho2b)*yv(ivec,ich3o2)*yv(ivec,iho2)&
                  +1.37*rrv(ivec,kmo2mo2)*yv(ivec,ich3o2)*yv(ivec,ich3o2)&
                  +rjv(ivec,jmepe)*yv(ivec,ich3o2h)&
                  +0.8*rrv(ivec,kc57)*yv(ivec,iole)*yv(ivec,ioh)&
                  +0.74*rrv(ivec,kc58)*yv(ivec,iole)*yv(ivec,io3)&
                  +rrv(ivec,kc59)*yv(ivec,iole)*yv(ivec,ino3)&
                  +1.56*rrv(ivec,kc61)*yv(ivec,ieth)*yv(ivec,ioh)&
                  +rrv(ivec,kc62)*yv(ivec,ieth)*yv(ivec,io3)&
                  +0.629*rrv(ivec,kc76)*yv(ivec,iisop)*yv(ivec,ioh)&
                  +0.6*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)&
                  +rrv(ivec,kohch3oh)*yv(ivec,ioh)*yv(ivec,ich3oh)&
                  +rrv(ivec,kno3mo2)*yv(ivec,ino3)*yv(ivec,ich3o2)&
                  +0.1*rrv(ivec,kohethoh)*yv(ivec,ioh)*yv(ivec,iethoh)&
                  +0.54*rrv(ivec,ko3c3h6)*yv(ivec,io3)*yv(ivec,ic3h6)&
                  +1.22*rrv(ivec,kohterp)*yv(ivec,ioh)*yv(ivec,iterp)&
                  +1.8*rrv(ivec,ko3terp)*yv(ivec,io3)*yv(ivec,iterp)&
                  +0.167*rrv(ivec,kohispd)*yv(ivec,ioh)*yv(ivec,iispd)&
                  +0.15*rrv(ivec,ko3ispd)*yv(ivec,io3)*yv(ivec,iispd)&
                  +0.282*rrv(ivec,kno3ispd)*yv(ivec,ino3)*yv(ivec,iispd)&
                  +0.9*rjv(ivec,jispd)*yv(ivec,iispd)&
                  +rrv(ivec,kaco2no)*yv(ivec,iaco2)*yv(ivec,ino)&
                  +rrv(ivec,khypo2no)*yv(ivec,ino)*yv(ivec,ihypro2)&
                  +rjv(ivec,jmo2no2b)*yv(ivec,ich3o2no2)           

               XLCH2O=rjv(ivec,jach2o)+rjv(ivec,jbch2o)+yv(ivec,ioh)*rrv(ivec,kfrmoh)&
                      +rrv(ivec,kc41)*yv(ivec,ino3)
               XLCH2O=XLCH2O + vdv(ivec,ich2o)
               yv(ivec,ich2o)=(y0v(ivec,ich2o)+PCH2O*DT)/(1.+XLCH2O*DT)
               ! --- CBIV-elements for higher HC-chemistry: ALD2 MGLY
               ! --- ETH OLE ISOP ROOH ORGNTR
               PALD2=0.11*rrv(ivec,kc52)*yv(ivec,ipar)*yv(ivec,ioh)&
                +1.1*rrv(ivec,kc53)*yv(ivec,iror)&
                +0.95*rrv(ivec,kc57)*yv(ivec,iole)*yv(ivec,ioh)&
                +0.5*rrv(ivec,kc58)*yv(ivec,iole)*yv(ivec,io3)&
                +0.91*rrv(ivec,kc59)*yv(ivec,iole)*yv(ivec,ino3)&
                +0.22*rrv(ivec,kc61)*yv(ivec,ieth)*yv(ivec,ioh)&
                +0.8*rrv(ivec,kc78)*yv(ivec,iisop)*yv(ivec,ino3)&
                +0.04*rrv(ivec,kohrooh)*yv(ivec,ioh)*yv(ivec,irooh)&
                +0.991*rrv(ivec,kohc2h6)*yv(ivec,ioh)*yv(ivec,ic2h6)&
                +0.3*rjv(ivec,jrooh)*yv(ivec,irooh)&
                +0.3*rrv(ivec,kc84)*yv(ivec,ioh)*yv(ivec,iorgntr)&
                +0.3*rjv(ivec,jorgn)*yv(ivec,iorgntr)&
                +rrv(ivec,kohethoh)*yv(ivec,ioh)*yv(ivec,iethoh)&
                +0.5*rrv(ivec,ko3c3h6)*yv(ivec,io3)*yv(ivec,ic3h6)&
                +0.47*rrv(ivec,kohterp)*yv(ivec,ioh)*yv(ivec,iterp)&
                +0.21*rrv(ivec,ko3terp)*yv(ivec,io3)*yv(ivec,iterp)&
                +0.47*rrv(ivec,kno3terp)*yv(ivec,ino3)*yv(ivec,iterp)&
                +0.273*rrv(ivec,kohispd)*yv(ivec,ioh)*yv(ivec,iispd)&
                +0.02*rrv(ivec,ko3ispd)*yv(ivec,io3)*yv(ivec,iispd)&
                +0.357*rrv(ivec,kno3ispd)*yv(ivec,ino3)*yv(ivec,iispd)&
                +0.067*rjv(ivec,jispd)*yv(ivec,iispd)&
                +0.7*rrv(ivec,kaco2mo2)*yv(ivec,iaco2)*yv(ivec,ich3o2)&
                +0.27*rrv(ivec,kc3h7o2no)*yv(ivec,ic3h7o2)*yv(ivec,ino)&
                +rrv(ivec,khypo2no)*yv(ivec,ihypro2)*yv(ivec,ino)               
                        
               XLALD2=rrv(ivec,kc43)*yv(ivec,ioh)+rrv(ivec,kc44)*yv(ivec,ino3)&
                    +rjv(ivec,jald2)
               XLALD2=XLALD2 + vdv(ivec,iald2)
               yv(ivec,iald2)=(y0v(ivec,iald2)+PALD2*DT)/(1.+XLALD2*DT)
               PMGLY=0.19*rrv(ivec,kohrooh)*yv(ivec,ioh)*yv(ivec,irooh)&
                 +0.168*rrv(ivec,kohispd)*yv(ivec,ioh)*yv(ivec,iispd)&
                 +0.85*rrv(ivec,ko3ispd)*yv(ivec,io3)*yv(ivec,iispd)&
                 +0.5*rrv(ivec,kaco2mo2)*yv(ivec,iaco2)*yv(ivec,ich3o2)
               XLMGLY=rrv(ivec,kc73)*yv(ivec,ioh)+rjv(ivec,jmgly)
               yv(ivec,imgly)=(y0v(ivec,imgly)+PMGLY*DT)/(1.+XLMGLY*DT)
               XLETH=rrv(ivec,kc61)*yv(ivec,ioh)+rrv(ivec,kc62)*yv(ivec,io3)
               yv(ivec,ieth)=y0v(ivec,ieth)/(1.+XLETH*DT)
               
               POLE=0.
               XLOLE=rrv(ivec,kc57)*yv(ivec,ioh)+rrv(ivec,kc58)*yv(ivec,io3)+rrv(ivec,kc59)*yv(ivec,ino3)
               yv(ivec,iole)=(y0v(ivec,iole)+POLE*DT)/(1.+XLOLE*DT)
               
               PROOH=rrv(ivec,kc82)*yv(ivec,ixo2)*yv(ivec,iho2)&
                    +rrv(ivec,kc85)*yv(ivec,iho2)*yv(ivec,ixo2n)&
                    +rrv(ivec,kaco2ho2)*yv(ivec,iaco2)*yv(ivec,iho2)&
                    +rrv(ivec,kc3h7o2ho2)*yv(ivec,ic3h7o2)*yv(ivec,iho2)&
                    +rrv(ivec,khypo2ho2)*yv(ivec,ihypro2)*yv(ivec,iho2)
               XLROOH=rjv(ivec,jrooh)+rrv(ivec,kohrooh)*yv(ivec,ioh)
               XLROOH = XLROOH + vdv(ivec,irooh)
               yv(ivec,irooh)=(y0v(ivec,irooh)+PROOH*DT)/(1.+XLROOH*DT)

              PORGNTR=rrv(ivec,kc81)*yv(ivec,ino)*yv(ivec,ixo2n)&
                +0.8*rrv(ivec,kc78)*yv(ivec,iisop)*yv(ivec,ino3)&
                +rrv(ivec,kno3c3h6)*yv(ivec,ic3h6)*yv(ivec,ino3)&
                +0.53*rrv(ivec,kno3terp)*yv(ivec,ino3)*yv(ivec,iterp)&
                +0.85*rrv(ivec,kno3ispd)*yv(ivec,ino3)*yv(ivec,iispd)
               XLORGNTR=rrv(ivec,kc84)*yv(ivec,ioh)+rjv(ivec,jorgn)
               XLORGNTR=XLORGNTR+vdv(ivec,iorgntr)

               yv(ivec,iorgntr)=(y0v(ivec,iorgntr)+PORGNTR*DT)/(1.+XLORGNTR*DT)
                       
               PACET=0.82*rrv(ivec,kc3h7o2no)*yv(ivec,ino)*yv(ivec,ic3h7o2)
               
               XLACET=rjv(ivec,ja_acet)+rjv(ivec,jb_acet)+rrv(ivec,kohacet)*yv(ivec,ioh)
               XLACET=XLACET+vdv(ivec,iacet)
                      
                yv(ivec,iacet)=(y0v(ivec,iacet)+PACET*DT)/(1.+XLACET*DT)       

               ! gas phase sulfur   & ammonia

               qdms1=rrv(ivec,kdmsoha)*yv(ivec,ioh)+rrv(ivec,kdmsno3)*yv(ivec,ino3)
               qdms2=rrv(ivec,kdmsohb)*yv(ivec,ioh)
               qdms=qdms1+qdms2
               yv(ivec,idms)=y0v(ivec,idms)/(1.+qdms*DT)
               pso2=yv(ivec,idms)*(qdms1+0.75*qdms2)
               pmsa=yv(ivec,idms)*0.25*qdms2
               qso2=rrv(ivec,kso2oh)*yv(ivec,ioh)
               qso2d=qso2 + vdv(ivec,iso2)
               yv(ivec,iso2)=(y0v(ivec,iso2)+pso2*DT) /(1.+qso2d*DT)  !qso2d includes deposition
               yv(ivec,imsa)=(y0v(ivec,imsa)+pmsa*DT) /(1.+vdv(ivec,imsa)*DT)  
#ifdef with_m7
               !VH: Do not apply dry deposition to H2SO4 : This deposition velocity represents aerosol deposition
               yv(ivec,iso4)=(y0v(ivec,iso4)+qso2*yv(ivec,iso2)*DT)   
#ifdef with_budgets
               ! leave out dt to get s-1 values
               temp_prod_so4(ivec,1)=qso2*yv(ivec,iso2)*DT

#endif
#else
               !VH: Do apply dry deposition to SO4_A : This deposition velocity does represent aerosol deposition
               !VH: Use the same aerosol deposition velocity for NO3_A deposition.
               yv(ivec,iso4)=(y0v(ivec,iso4)+qso2*yv(ivec,iso2)*DT) /(1. + vdv(ivec,iso4)*DT)   !corrected CMK qso2/qso2d
               yv(ivec,ino3_a)=y0v(ivec,ino3_a) /(1.+vdv(ivec,ino3_a)*DT)  ! VH/FD include dry dep for no3_a? (should be copy of so4)
#endif
               yv(ivec,inh4)=y0v(ivec,inh4)/(1.+vdv(ivec,inh4)*DT) ! This is just deposition
               pnh2=yv(ivec,inh2o2)*yv(ivec,ino)*rrv(ivec,knh2o2no)&
                    +yv(ivec,inh2o2)*yv(ivec,io3)*rrv(ivec,knh2o2o3)&
                    +yv(ivec,inh2o2)*yv(ivec,iho2)*rrv(ivec,knh2o2ho2)
               pnh3=rrv(ivec,Knh2ho2)*yv(ivec,inh2)*yv(ivec,iho2) 
               dnh3=yv(ivec,ioh)*rrv(ivec,knh3oh) + vdv(ivec,inh3)
               yv(ivec,inh3)=(y0v(ivec,inh3)+pnh3*DT)/(1.+dnh3*DT)
               ppnh3=yv(ivec,ioh)*yv(ivec,inh3)*rrv(ivec,knh3oh)
               qnh2= rrv(ivec,knh2oh)*yv(ivec,ioh)+rrv(ivec,knh2no)*yv(ivec,ino)&
                    +rrv(ivec,knh2no2)*yv(ivec,ino2)+rrv(ivec,knh2ho2)*yv(ivec,iho2)&
                    +rrv(ivec,knh2o3)*yv(ivec,io3)&
                    +rrv(ivec,knh2o2)
               yv(ivec,inh2)=(y0v(ivec,inh2)+ppnh3*dt+pnh2*dt)/(1.+qnh2*dt)
	       !
	       ! Now nh2o2 radical
	       !
	       qnh2o2=rrv(ivec,knh2o2no)*yv(ivec,ino)+rrv(ivec,knh2o2o3)*yv(ivec,io3)+rrv(ivec,knh2o2ho2)*yv(ivec,iho2)
	       yv(ivec,inh2o2)=(y0v(ivec,inh2o2)+(rrv(ivec,knh2o3)*yv(ivec,io3)*yv(ivec,inh2)*DT))/(1.+qnh2o2*DT)	       

             end do  !ivec

         end if
         
         if ( mod(iter,maxit) == 0 ) then

            ! --- Long living compounds
            do ivec=1,lvec
               yv(ivec,ich4)=y0v(ivec,ich4)/(1.+rrv(ivec,kch4oh)*yv(ivec,ioh)*DT)
               ! ch4loss=1.0-1.0/((1.+rrv(ivec,kch4oh)*yv(ivec,ioh)*DT))
               !methane loss?
               PCO=yv(ivec,ich2o)*(rjv(ivec,jach2o)+rjv(ivec,jbch2o)&
                +yv(ivec,ioh)*rrv(ivec,kfrmoh))&
                +rjv(ivec,jald2)*yv(ivec,iald2)&
                +rjv(ivec,jmgly)*yv(ivec,imgly)&
                +0.62*rrv(ivec,kc57)*yv(ivec,iole)*yv(ivec,ioh)&
                +0.65*rrv(ivec,kc58)*yv(ivec,iole)*yv(ivec,io3)&
                +0.56*rrv(ivec,kc59)*yv(ivec,iole)*yv(ivec,ino3)&
                +0.24*rrv(ivec,kc62)*yv(ivec,ieth)*yv(ivec,io3)&
                +0.066*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)&
                +rrv(ivec,kc41)*yv(ivec,ich2o)*yv(ivec,ino3)&
                +0.56*rrv(ivec,ko3c3h6)*yv(ivec,io3)*yv(ivec,ic3h6)&
                +0.47*rrv(ivec,kohterp)*yv(ivec,ioh)*yv(ivec,iterp)&
                +0.211*rrv(ivec,ko3terp)*yv(ivec,io3)*yv(ivec,iterp)&
                +0.47*rrv(ivec,kno3terp)*yv(ivec,ino3)*yv(ivec,iterp)&
                +0.334*rrv(ivec,kohispd)*yv(ivec,ioh)*yv(ivec,iispd)&
                +0.225*rrv(ivec,ko3ispd)*yv(ivec,io3)*yv(ivec,iispd)&
                +0.643*rrv(ivec,kno3ispd)*yv(ivec,ino3)*yv(ivec,iispd)&
                +0.333*rjv(ivec,jispd)*yv(ivec,iispd)&
                +rjv(ivec,ja_acet)*yv(ivec,iacet)       
               XLCO = rrv(ivec,kcooh)*yv(ivec,ioh)
               XLCO = XLCO + vdv(ivec,ico)
               yv(ivec,ico)=(y0v(ivec,ico)+PCO*DT)/(1.+XLCO*DT)
               !carbon monoxide loss?
               pch3oh=(0.63*rrv(ivec,kmo2mo2)*yv(ivec,ich3o2)*yv(ivec,ich3o2))+&
                      (0.5*rrv(ivec,kaco2mo2)*yv(ivec,ich3o2)*yv(ivec,iaco2))

               yv(ivec,ich3oh)=(y0v(ivec,ich3oh)+pch3oh*dt)/&
                     (1.+(vdv(ivec,ich3oh)+rrv(ivec,kohch3oh)*yv(ivec,ioh))*dt)
                     
               phcooh=(0.52*rrv(ivec,kc62)*yv(ivec,io3)*yv(ivec,ieth)) +&
                      (0.25*rrv(ivec,ko3c3h6)*yv(ivec,io3)*yv(ivec,ic3h6))           
                     
               yv(ivec,ihcooh)=(y0v(ivec,ihcooh)+(phcooh*dt))/&
                     (1.+rrv(ivec,kohhcooh)*yv(ivec,ioh)*dt)
                     
               pmcooh=0.4*rrv(ivec,kc50)*yv(ivec,ic2o3)*yv(ivec,iho2)
             
               yv(ivec,imcooh)=(y0v(ivec,imcooh)+(pmcooh*dt))/&
                     (1.+(vdv(ivec,imcooh)+rrv(ivec,kohmcooh)*yv(ivec,ioh))*dt)                      
                     
               yv(ivec,ic2h6)=y0v(ivec,ic2h6)/(1.+rrv(ivec,kohc2h6)*yv(ivec,ioh)*dt)
           
               yv(ivec,iethoh)=(y0v(ivec,iethoh)/(1.+(vdv(ivec,iethoh)+rrv(ivec,kohethoh)*yv(ivec,ioh))*dt))
           
               yv(ivec,ic3h8)=(y0v(ivec,ic3h8)/(1.+rrv(ivec,kohc3h8)*yv(ivec,ioh)*dt))
          
               xlc3h6=rrv(ivec,kohc3h6)*yv(ivec,ioh)&
                  +rrv(ivec,ko3c3h6)*yv(ivec,io3)&
                  +rrv(ivec,kno3c3h6)*yv(ivec,ino3)
               yv(ivec,ic3h6)=(y0v(ivec,ic3h6)/(1.+xlc3h6*dt))
                   
               ppar=0.35*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)&
                +2.4*rrv(ivec,kc78)*yv(ivec,iisop)*yv(ivec,ino3)&
                +5.0*rrv(ivec,kohterp)*yv(ivec,ioh)*yv(ivec,iterp)&
                +6.0*rrv(ivec,ko3terp)*yv(ivec,io3)*yv(ivec,iterp)&
                +6.0*rrv(ivec,kno3terp)*yv(ivec,ino3)*yv(ivec,iterp)&   
                +1.565*rrv(ivec,kohispd)*yv(ivec,ioh)*yv(ivec,iispd)&
                +0.36*rrv(ivec,ko3ispd)*yv(ivec,io3)*yv(ivec,iispd)&
                +1.282*rrv(ivec,kno3ispd)*yv(ivec,ino3)*yv(ivec,iispd)&
                +0.832*rjv(ivec,jispd)*yv(ivec,iispd)&
                +0.6*rrv(ivec,kaco2mo2)*yv(ivec,ich3o2)*yv(ivec,iaco2)  
               xlpar=rrv(ivec,kc52)*yv(ivec,ioh)+rrv(ivec,kc83)*yv(ivec,irxpar)
               yv(ivec,ipar)=(y0v(ivec,ipar)+ppar*dt)/(1.+xlpar*dt)
               
               !cmk ____added rn222 chemistry in EBI language
               yv(ivec,irn222) = y0v(ivec,irn222)/(1.+rrv(ivec,krn222)*dt)
               yv(ivec,ipb210) = y0v(ivec,ipb210)+y0v(ivec,irn222)-yv(ivec,irn222)

            end do   !ivec

         end if  
                 
      end do !ITER

    end subroutine DO_EBI



    subroutine NOYmass

      integer i,j,imax
      real :: ncormax,ncorav,totn,totn0,fnoy,fnoy1
      real :: ncorr,ncorr1,ncorr2,ncorr3,ncorr4,ncorr5, totdep
      logical :: nerr

      ncormax=0.
      ncorav=0.
      nerr=.false.
      imax = 0
      
      do j=j1,j2
      do i=i1,i2
            imax = imax + 1 
            !
            !** Guarantee exact mass conservation of NOY 
            !   (this may matter a few percent)
            ! 
            fnoy=y(i,j,ino)+y(i,j,ino2)+y(i,j,ino3)+2.*y(i,j,in2o5)+y(i,j,ihno4)
            if (level == 1) then
#ifndef without_dry_deposition
               totdep = (y(i,j,ino) *vd(region,ino )%surf(i,j)  + &
                         y(i,j,ino2)*vd(region,ino2)%surf(i,j)  + &
                         y(i,j,ino3)*vd(region,ino3)%surf(i,j)  + &
                         y(i,j,ihno3)*vd(region,ihno3)%surf(i,j)  + &
                         y(i,j,ipan)*vd(region,ipan)%surf(i,j)  + &
                         y(i,j,iorgntr)*vd(region,iorgntr)%surf(i,j)  + &
                         y(i,j,ich3o2no2)*vd(region,ich3o2no2)%surf(i,j)  + &
                         2*y(i,j,in2o5)*vd(region,in2o5)%surf(i,j)  + &
                         y(i,j,ihno4)*vd(region,ihno4)%surf(i,j) + &
                         y(i,j,ihono)*vd(region,ihono)%surf(i,j))*dt/ye(i,j,idz)  
#else
               totdep = 0.0
#endif
            else
               totdep  = 0.0 
            end if
            totn0=y0(i,j,inox)+y0(i,j,ihno3)+y0(i,j,ipan)+y0(i,j,ihono)+ &
                 y0(i,j,iorgntr) + y0(i,j,ich3o2no2) + ye(i,j,ieno)*dt -  totdep  
            ! note that emino is added here and the total deposition is subtracted
            !
            ! totn0 contains all nitrogen at beginning of timestep + nox emissions
            !
            !
            ! totn contains all nitrogen at end of timestep
            !
            totn=fnoy+y(i,j,ihno3)+y(i,j,ipan)+y(i,j,iorgntr)+y(i,j,ihono)+ y(i,j,ich3o2no2)
            ! correction factor for all nitrogen compounds
            ncorr=totn-totn0

            if ( totn < tiny(totn) ) cycle

            if ( (abs(ncorr)/totn) > 0.05 ) then    !CMK changed from 0.1 to 0.05

               nerr=.true.

               !AJS>>>
               !               print *,'NOYmass: N-error....',region,level,i,j,totn0,totn
               !               print *,'NOYmass: emino    ',ye(i,j,ieno)*dt/y0(i,j,iair)*1e9
               !               print *,'NOYmass: NO(0)    ', &
               !                    y0(i,j,ino)/y0(i,j,iair)*1e9,y(i,j,ino)/y(i,j,iair)*1e9
               !               print *,'NOYmass: NO2(0)   ', &
               !                    y0(i,j,ino2)/y0(i,j,iair)*1e9,y(i,j,ino2)/y(i,j,iair)*1e9
               !               print *,'NOYmass: O3(0)    ', &
               !                    y0(i,j,io3)/y0(i,j,iair)*1e9,y(i,j,io3)/y(i,j,iair)*1e9
               !               print *,'NOYmass: NO3(0)   ', &
               !                    y0(i,j,ino3)/y0(i,j,iair)*1e9,y(i,j,ino3)/y(i,j,iair)*1e9
               !               print *,'NOYmass: N2O5(0)  ', &
               !                    y0(i,j,in2o5)/y0(i,j,iair)*1e9,y(i,j,in2o5)/y(i,j,iair)*1e9
               !               print *,'NOYmass: HNO4(0)  ', &
               !                    y0(i,j,ihno4)/y0(i,j,iair)*1e9,y(i,j,ihno4)/y(i,j,iair)*1e9
               !               print *,'NOYmass: HNO3(0)  ', &
               !                    y0(i,j,ihno3)/y0(i,j,iair)*1e9,y(i,j,ihno3)/y(i,j,iair)*1e9
               !               print *,'NOYmass: PAN(0)   ', &
               !                    y0(i,j,ipan)/y0(i,j,iair)*1e9,y(i,j,ipan)/y(i,j,iair)*1e9
               !               print *,'NOYmass: ORGNT(0) ', &
               !                    y0(i,j,iorgntr)/y0(i,j,iair)*1e9,y(i,j,iorgntr)/y(i,j,iair)*1e9
               !               print *,'NOYmass: NOx(0)   ', &
               !                    y0(i,j,inox)/y0(i,j,iair)*1e9,y(i,j,inox)/y(i,j,iair)*1e9
               !               print *,'NOYmass: ',rj(i,j,jhno3),rr(i,j,kohhno3)*y(i,j,ioh), &
               !                    y(i,j,ioh)/y(i,j,iair)*1e9
               !<<<AJS

            end if
            ! maximum and average correction factor in this loop
            ncormax=max(abs(ncormax),abs(ncorr/totn)) 
            ncorav=ncorav+abs(ncorr/totn)
            !
            ! first correct hno3, pan and organic nitrates (added hono)
            ! (as a group of reservoir tracers)
            ! 
            totn=y(i,j,ihno3)+y(i,j,ipan)+y(i,j,iorgntr)+y(i,j,ihono)+y(i,j,ich3o2no2)
            if ( totn < tiny(totn) ) cycle
            ncorr1=y(i,j,ihno3)  *(1.-ncorr/totn)
            ncorr2=y(i,j,ipan)   *(1.-ncorr/totn)
            ncorr3=y(i,j,iorgntr)*(1.-ncorr/totn)
            ncorr4=y(i,j,ihono)*(1.-ncorr/totn)
            ncorr5=y(i,j,ich3o2no2)*(1.-ncorr/totn)
            y(i,j,ihno3)    =max(0.,ncorr1)
            y(i,j,ipan)     =max(0.,ncorr2)
            y(i,j,iorgntr)  =max(0.,ncorr3)
            y(i,j,ihono)    =max(0.,ncorr4)
            y(i,j,ich3o2no2)=max(0.,ncorr5)
            
            ncorr=ncorr1+ncorr2+ncorr3+ncorr4+ncorr5-y(i,j,ihno3)-y(i,j,ipan)-y(i,j,iorgntr)-y(i,j,ihono)-y(i,j,ich3o2no2)
            !
            ! the remainder is used to scale the noy components
            !
            fnoy1=(fnoy+ncorr)/fnoy
            y(i,j,ino)  =fnoy1*y(i,j,ino)
            y(i,j,ino2) =fnoy1*y(i,j,ino2)
            y(i,j,ino3) =fnoy1*y(i,j,ino3)
            y(i,j,in2o5)=fnoy1*y(i,j,in2o5)
            y(i,j,ihno4)=fnoy1*y(i,j,ihno4)
            y(i,j,inox) =y(i,j,ino)+y(i,j,ino2)+y(i,j,ino3)+2.*y(i,j,in2o5)+y(i,j,ihno4)

         end do

      end do

      if ( nerr ) then
         write (gol,'("NOYmass: N-mass balance error, ncorr>5% ")'); call goPr
         write (gol,'("  Maximum correction                     : ",f8.2)') ncormax; call goPr
         write (gol,'("  Average correction in this loop (imax) : ",f8.2," (",i6,")")') ncorav/imax, imax; call goPr
      end if

    end subroutine NOYmass

#ifdef with_budgets
    !--------------------------------------------------------------------------
    !                    TM5                                                  !
    !--------------------------------------------------------------------------
    !BOP
    !
    ! !IROUTINE:  
    !
    ! !DESCRIPTION: increase reaction budgets for each reaction
    !         arrays nrr and nrj determine which species are
    !         involved in a reaction
    !\\
    !\\
    ! !INTERFACE:
    !
    subroutine INCC2C3
      !
#ifdef with_tendencies  
      use TRACER_DATA,   only : PLC_AddValue, plc_ipr_lddep, plc_kg_from_tm
#endif

      !      integer, intent(out)    ::  status

      integer :: i01,n1,n2,jl,i,j
      ! nrj and nrr used for reaction budget calculations
      integer,dimension(nj),parameter    ::  nrj=(/io3,ino2,ih2o2,ihno3,ihno4,in2o5,in2o5,ich2o,ich2o, &
           ich3o2h,ino3,ino3,ipan,ipan,iorgntr,iald2,imgly,irooh,io2,iispd,iacet,iacet,  &
           ihono,ich3o2no2,ich3o2no2/)
      integer,dimension(nreac,2),parameter  ::  nrr = reshape((/ &
           ino    ,iho2  ,ich3o2,ino2   ,ioh    ,ino2   ,ino   ,ino2  ,in2o5, ihno4 , &
           ino2   ,ihno4 ,iair  ,ih2o   ,io3    ,ico    ,io3   ,ih2o2 ,ich2o, ich4  , &
           ioh    ,ioh   ,ich3o2,ich3o2 ,ich3o2 ,iho2   ,iho2  ,in2o5 ,ioh  , ich2o , &
           iald2  ,iald2 ,ic2o3 ,ic2o3  ,ipan   ,ic2o3  ,ic2o3 ,ipar  ,iror , iror  , &
           ioh    ,io3   ,ino3  ,ioh    ,io3    ,ioh    ,ioh   ,io3   ,ino3 , ixo2  , &
           ixo2   ,ixo2n ,ixo2  ,irxpar ,iorgntr,ixo2n  ,idms  ,idms  ,idms , iso2  , &
           inh3   ,inh3  ,inh2  ,inh2  ,inh2   ,inh2   ,inh2   ,inh2  ,ioh   ,ioh  , ino3  , &
           ino3   ,ino3  ,ino3  ,ioh    ,ioh    ,ioh    ,ioh   ,ioh   ,io3  , ino3  , &
           ioh    ,io3   ,ino3  ,ioh    ,io3    ,ino3   ,irn222,io3   ,iair , iacet , &
           iaco2  ,iaco2 ,iaco2 ,iaco2  ,ixo2   ,ixo2n  ,ino   ,iho2  ,ino  , iho2  , &
           in2o5  ,ino3  ,iho2  ,iho2   ,ino    ,io3    ,iho2, ioh   ,ioh   ,ich3o2,ich3o2no2, &
           !second reaction partner (if monmolecular = 0)
           io3    ,ino   ,ino   ,ioh    ,ihno3  ,io3    ,ino3   ,ino3  ,0      ,ioh    , &
           iho2   ,0     ,0     ,0      ,iho2   ,ioh    ,ioh    ,ioh   ,ioh    ,ioh    , &
           ich3o2h,irooh ,iho2  ,iho2   ,ich3o2 ,ioh    ,iho2   ,0     ,0      ,ino3   , &
           ioh    ,ino3  ,ino   ,ino2   ,0      ,ic2o3 ,iho2   ,ioh    ,0      ,0      , &
           iole   ,iole  ,iole  ,ieth   ,ieth   ,imgly ,iisop  ,iisop  ,iisop  ,ino    , &
           ixo2   ,ino   ,iho2  ,ipar   ,ioh    ,iho2  ,ioh    ,ioh    ,ino3   ,ioh    , &
           iacid  ,ioh   ,ioh   ,ino   ,ino2   ,iho2   ,0     ,io3    ,ich3oh ,ihcooh ,iho2   , &
           ich3o2 ,ic2o3 ,ixo2  ,imcooh ,ic2h6  ,iethoh,ic3h8  ,ic3h6  ,ic3h6  ,ic3h6  , &
           iterp  ,iterp ,iterp ,iispd  ,iispd  ,iispd ,0      ,iair   ,0      ,ioh    , &
           iho2   ,ich3o2,ino   ,ixo2   ,ixo2n  ,ixo2n ,ic3h7o2,ic3h7o2,ihypro2,ihypro2, &
           0      ,0     ,0     ,0      ,inh2o2 ,inh2o2,inh2o2, ino    ,ihono  ,ino2   ,0/),(/nreac,2/))

      real :: c1,xdep

      c1=dt*1000./xmair   !conversion to moles...

      ! Reaction budgets

      do i01=1,nj !photolysis rates
         n1=nrj(i01)
         do j=j1,j2
            do i=i1,i2
               if(n1 > 0) cr2(i,j,i01)=cr2(i,j,i01)+rj(i,j,i01)*y(i,j,n1)
            end do
         end do
      end do!i01=1,nj
      !
      do i01=1,nreac !reactions
         n1=nrr(i01,1) !make sure n1 > 0
         n2=nrr(i01,2)
         if (n2 > 0.) then
            do j=j1,j2
               do i=i1,i2
                  cr3(i,j,i01)= cr3(i,j,i01)+y(i,j,n1)*y(i,j,n2)*rr(i,j,i01)
               end do
            end do
         else
            do j=j1,j2
               do i=i1,i2
                  cr3(i,j,i01)= cr3(i,j,i01)+y(i,j,n1)*rr(i,j,i01)
               end do
            end do
         end if
      end do  !i01=1,nreac

      if ( level == 1 ) then   ! deposition budget

         do i01=1,ntrace
            if (fscale(i01) > 0) then 
               do j=j1,j2
                  do i=i1,i2
#ifndef without_dry_deposition
                     if (i01 .ne. inox) then
                        xdep = y(i,j,i01)*vd(region,i01)%surf(i,j)/ye(i,j,idz)* &
                             c1*ye(i,j,iairm)/y(i,j,iair)  !from updated concentrations
                     else ! compute nox deposition from other contributors!
                        xdep = (y(i,j,ino)  *vd(region,ino)%surf(i,j) + &
                             y(i,j,ino2 )*vd(region,ino2)%surf(i,j)+ &
                             y(i,j,ino3) *vd(region,ino3)%surf(i,j)+ &
                             y(i,j,ihno3)*vd(region,ihno3)%surf(i,j)+ &
                             2*y(i,j,in2o5)*vd(region,in2o5)%surf(i,j) + &
                             y(i,j,ihono)*vd(region,ihono)%surf(i,j)) &
                             /ye(i,j,idz)* &
                             c1*ye(i,j,iairm)/y(i,j,iair)  !from updated concentrations
                     endif
#else
                     xdep = 0.0
#endif
                     buddep_dat(region)%dry(i,j,i01) = &
                          buddep_dat(region)%dry(i,j,i01) + xdep
                     if ( i01 == 1 ) then   !seperate deposition from 'other' chemistry
#ifndef without_dry_deposition
                        sum_deposition(region) = sum_deposition(region) - &
                             xdep*ra(1)*1e-3   ! in kg
#endif
                        sum_chemistry(region) = sum_chemistry(region) +  &
                             (y(i,j,1)-y0(i,j,1))/y(i,j,iair)* &
                             ye(i,j,iairm)/xmair*ra(1) + xdep*ra(1)*1e-3  
                     end if

                     ! FIXME TENDENCIES
#ifdef with_tendencies  
                     ! Add deposition budget in kg to tendencies;
                     !   (mole tm tracer) * (g/mole) * (kg/g) = kg tm tracer
                     call PLC_AddValue( region, plc_ipr_lddep, i, j, 1, i01, &
                          xdep * ra(i01) * 1e-3 * plc_kg_from_tm(i01), &   ! kg plc tracer
                          status )
!                     IF_NOTOK_RETURN(status=1)
#endif

                  end do
               end do
            endif
         end do   !i01

      else   ! other layers

         do j=j1,j2
            do i=i1,i2
               sum_chemistry(region) = sum_chemistry(region) + &
                    (y(i,j,1)-y0(i,j,1))/y(i,j,iair)*ye(i,j,iairm)/xmair*ra(1)
            end do
         end do
      end if   !level ==1

    end subroutine INCC2C3

    !--------------------------------------------------------------------------
    !                    TM5                                                  !
    !--------------------------------------------------------------------------
    !BOP
    !
    ! !IROUTINE:  REACBUD
    !
    ! !DESCRIPTION: accumulate budgets, o3 P/L
    !\\
    !\\
    ! !INTERFACE:
    !
    SUBROUTINE REACBUD
      !
      ! !USES:
      !
      USE BUDGET_GLOBAL, ONLY : budg_dat, nzon_vg
      !
      !EOP
      !------------------------------------------------------------------------
      !BOC

      integer :: i01, i, j, nzone, nzone_v
      real    :: c1

      c1=dt*1000./xmair   !conversion to moles

      do j=j1,j2
         do i=i1,i2
            nzone=budg_dat(region)%nzong(i,j)    !global budget
            nzone_v=nzon_vg(level)    !level is passed to ebi...
            do i01=1,nj
               budrjg(nzone,nzone_v,i01)=budrjg(nzone,nzone_v,i01)+ &
                    cr2(i,j,i01)*c1*ye(i,j,iairm)/y(i,j,iair)   !units mole
            end do !nj
            ! ozone production on a 3D grid:
            ! defined as: HO2 + NO, CH3O2 + NO, XO2 + NO, C2O3 + NO
            o3p(region)%d3(i,j,level) = o3p(region)%d3(i,j,level) + &
                 (cr3(i,j,2) + cr3(i,j,3) + cr3(i,j,33) + cr3(i,j,50))   &
                 *c1*ye(i,j,iairm)/y(i,j,iair)   ! acc. moles/gridbox

            !AERCHEMMIP output
            AC_O3_lp(region)%prod(i,j,level,1)=AC_O3_lp(region)%prod(i,j,level,1)+&
                 (cr3(i,j,2)+ cr3(i,j,3) + cr3(i,j,33) + cr3(i,j,50))   &
                 *c1*ye(i,j,iairm)/y(i,j,iair) 

            o3l(region)%d3(i,j,level) = o3l(region)%d3(i,j,level) + &
                 (cr3(i,j,4)-2*cr3(i,j,7) + 2*cr3(i,j,6) + cr3(i,j,8) - cr3(i,j,9) + &
                 cr3(i,j,15) + cr3(i,j,17)  + cr3(i,j,42) - cr3(i,j,43)  + &
                 cr3(i,j,45) + cr3(i,j,48)  - 0.1*cr3(i,j,49) - cr3(i,j,55) + &
                 cr2(i,j,1) - cr2(i,j,4) - cr2(i,j,6) - cr2(i,j,13) - 2.*cr2(i,j,10)) &
                 *c1*ye(i,j,iairm)/y(i,j,iair)   ! acc. moles/gridbox

            !AERCHEMMIP output
            AC_O3_lp(region)%prod(i,j,level,2)=AC_O3_lp(region)%prod(i,j,level,2)+&
                  (cr3(i,j,4)-2*cr3(i,j,7) + 2*cr3(i,j,6) + cr3(i,j,8) - cr3(i,j,9) + &
                 cr3(i,j,15) + cr3(i,j,17)  + cr3(i,j,42) - cr3(i,j,43)  + &
                 cr3(i,j,45) + cr3(i,j,48)  - 0.1*cr3(i,j,49) - cr3(i,j,55) + &
                 cr2(i,j,1) - cr2(i,j,4) - cr2(i,j,6) - cr2(i,j,13) - 2.*cr2(i,j,10)) &
                 *c1*ye(i,j,iairm)/y(i,j,iair)   ! acc. moles/gridbox


            o3l(region)%d3(i,j,level) = o3l(region)%d3(i,j,level) - &
                 cr4(i,j,1)*(1000./xmair)*ye(i,j,iairm)/y(i,j,iair)   !O3 + SO2 (note negative)
            
            AC_O3_lp(region)%prod(i,j,level,2)=AC_O3_lp(region)%prod(i,j,level,2)+&
                 cr4(i,j,1)*(1000./xmair)*ye(i,j,iairm)/y(i,j,iair)   !O3 + SO2 (note negative)


            !write(4321,*) AC_O3_lp(region)%prod(i,j,level,2)
            !write(4322,*)  (cr3(i,j,2) + cr3(i,j,3) + cr3(i,j,33) + cr3(i,j,50))   &
            !     *c1*ye(i,j,iairm)/y(i,j,iair) 

            !write(4323,*) (cr3(i,j,4)-2*cr3(i,j,7) + 2*cr3(i,j,6) + cr3(i,j,8) - cr3(i,j,9) + &
            !     cr3(i,j,15) + cr3(i,j,17)  + cr3(i,j,42) - cr3(i,j,43)  + &
            !     cr3(i,j,45) + cr3(i,j,48)  - 0.1*cr3(i,j,49) - cr3(i,j,55) + &
            !     cr2(i,j,1) - cr2(i,j,4) - cr2(i,j,6) - cr2(i,j,13) - 2.*cr2(i,j,10)) &
            !     *c1*ye(i,j,iairm)/y(i,j,iair)   ! acc. moles/gridbox


            !write(4324,*)   cr4(i,j,1)*(1000./xmair)*ye(i,j,iairm)/y(i,j,iair)   !O3 + SO2 (note negative)



            !PLS            ! ch4+oh
            !PLS            ch4oh(region)%d3(i,j,level) = ch4oh(region)%d3(i,j,level) + &
            !PLS                 cr3(i,j,20)*c1*ye(i,j,iairm)/y(i,j,iair)   ! acc. moles/gridbox
            !PLS            ! S gas phase production (OH + SO2---> SO4, OH + DMS = 0.25 MSA)
            !PLS            so4pg(region)%d3(i,j,level) = so4pg(region)%d3(i,j,level) + &
            !PLS                         (0.25*cr3(i,j,58) + cr3(i,j,60))*c1*ye(i,j,iairm)/y(i,j,iair)   ! acc. moles/gridbox OH + SO2
            !PLS            ! SO4 wet production
            !PLS            so4pa(region)%d3(i,j,level) = so4pa(region)%d3(i,j,level) - &
            !PLS                         (cr4(i,j,1)+cr4(i,j,2))*(1000./xmair)*ye(i,j,iairm)/y(i,j,iair)   ! acc. moles/gridbox note neg.
            do i01=1,nreac
               budrrg(nzone,nzone_v,i01)=budrrg(nzone,nzone_v,i01)+ &
                    cr3(i,j,i01)*c1*ye(i,j,iairm)/y(i,j,iair)   !units mole
            end do
            do i01=1,nreacw
               budrwg(nzone,nzone_v,i01)=budrwg(nzone,nzone_v,i01)- &
                    cr4(i,j,i01)*(1000./xmair)*ye(i,j,iairm)/y(i,j,iair)   ! mole
               !note: changed sign to get 'positive' budget, just a
               !      matter of definition, !CMK
            end do

         end do
      end do

      !sum up global budgets (done in chemistry/chemistry_done)

    end subroutine REACBUD
    !EOC
#endif
    
  end subroutine EBI
  !EOC


  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    WETS
  !
  ! !DESCRIPTION: aqueous phase chemistry of sulfur (and other): oxidation of
  !               SO2 and uptake of other gases in the aqueous phase 
  !               Method : implicit solution of oxidation of SO2
  !\\
  !\\
  ! !INTERFACE:
  !
#ifdef with_budgets
  subroutine wetS(region, level, is, js, y0, dt, y, ye, c4, status)
#else 
  subroutine wetS(region, level, is, js, y0, dt, y, ye,     status)
#endif
    !
    ! !USES:
    !
    use Dims,           only: CheckShape, idatei
    use global_data,    only: region_dat
    use reaction_data,  only: nreacw, ntlow, kso2hp, kso2o3
    use chem_param
    use dims,           only: im, jm
    use Binas,          only: Avog
    use boundary,       only: LCMIP6_CO2, co2_glob
    !
    ! !INPUT PARAMETERS:
    !
    integer,intent(in)  :: region        !region  region under operation (provides im,jm,lm via chemistry.mod)
    integer,intent(in)  :: level, is, js ! vertical level, tile start indices
    real,   intent(in)  :: dt            ! chemistry timestep
    real,   intent(in)  :: y0(is:,js:,:) ! initial concentration
    !
    ! !OUTPUT PARAMETERS:
    !
    real,   intent(out) :: y(is:,js:,:) ! concentrations at time is t
    integer,intent(out) :: status
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    real, intent(inout) :: ye(is:,js:,:)   ! extra fields (temp, cc, pH)
#ifdef with_budgets
    real, intent(inout) :: c4(is:,js:,:)   ! budget accumulatior
#endif
    !
    ! !REVISION HISTORY: 
    !          ???? - Ad Jeuken (KNMI) and Frank Dentener (IMAU) - v0
    !      Jan 2002 - Maarten Krol (IMAU)      - adapted for TM5
    !      Feb 2007 - Elisabetta Vignati (JRC) - adapted for TM5/M7
    !   22 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/Wets'
    
    integer n,i,j,l,itemp,iter
    real x1,x2,x3,b1,b2,so2x,dh2o2,dso2,disc,dnh3,dn2o5,xso2o3a,xso2o3b,co2
    real,parameter :: rg=0.08314
    real,dimension(:,:),allocatable :: hkso2  ! Henry's constant for sulfur dioxide
    real,dimension(:,:),allocatable :: hkh2o2 ! Henry's constant for hydroperoxide
    real,dimension(:,:),allocatable :: hko3       ! Henry's constant for ozone
    real,dimension(:,:),allocatable :: dkso2      ! Dissociation constant for SO2
    real,dimension(:,:),allocatable :: dkhso3     ! Dissociation constant for HSO3-
    real,dimension(:,:),allocatable :: dkh2o      ! dissociation constant water
    real,dimension(:,:),allocatable :: dknh3      ! dissociation constant ammonia
    real,dimension(:,:),allocatable :: hknh3      ! Henry's constant ammonia
    real,dimension(:,:),allocatable :: hkco2      ! Henry's constant CO2
    real,dimension(:,:),allocatable :: dkco2      ! Dissociation constant CO2
    real phs4                       ! effective dissolvation of S(IV)
    real phso2                      ! effective dissolvation of SO2
    real phh2o2                     ! effective dissolvation of H2O2
    real phozone                    ! effective dissolvation of O3
    real,dimension(:,:),allocatable :: hplus      !concentration h+
    real,dimension(:,:),allocatable :: sulph  !accumul+coarse mode sulphate
    real a1,a2,a,b,c,z,ft           ! help variables
    real xcov,xliq,xl,temp,rt,ztr,h2o,air,press ! meteo
    real,dimension(:,:,:),allocatable   :: rw ! reaction rates
    logical,dimension(:,:),allocatable ::  cloudy

    integer            :: i1, i2, j1, j2
    real               :: l_sum_wet

    ! --- begin --------------------------------

    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    
    ! check arguments ...
    call CheckShape( (/i2-i1+1, j2-j1+1, maxtrace/), shape(y0), status )
    IF_NOTOK_RETURN(status=1)
    call CheckShape( (/i2-i1+1, j2-j1+1, maxtrace/), shape(y ), status )
    IF_NOTOK_RETURN(status=1)
    call CheckShape( (/i2-i1+1, j2-j1+1, n_extra /), shape(ye), status )
    IF_NOTOK_RETURN(status=1)
#ifdef with_budgets
   call CheckShape( (/i2-i1+1, j2-j1+1, nreacw  /), shape(c4), status )
   IF_NOTOK_RETURN(status=1)
#endif

    allocate(hkso2      (i1:i2, j1:j2))
    allocate(hkh2o2     (i1:i2, j1:j2))
    allocate(hko3       (i1:i2, j1:j2))
    allocate(dkso2      (i1:i2, j1:j2))
    allocate(dkhso3     (i1:i2, j1:j2))
    allocate(dkh2o      (i1:i2, j1:j2))
    allocate(dknh3      (i1:i2, j1:j2))
    allocate(hknh3      (i1:i2, j1:j2))
    allocate(hkco2      (i1:i2, j1:j2))
    allocate(dkco2      (i1:i2, j1:j2))
    allocate(hplus      (i1:i2, j1:j2))
    allocate(rw         (i1:i2, j1:j2, nreacw))
    allocate(cloudy     (i1:i2, j1:j2))
    allocate(sulph      (i1:i2, j1:j2)) 

    !-----------------------------
    ! wet phase reactions
    !-----------------------------
    rw   =0.0 
    hplus=0.0

    if (LCMIP6_CO2) then
       co2=co2_glob
    else
!    
! JEW: now scaled to 2000 to account for annual growth of ~2ppbv/yr-1
!    
       co2=3.69e-4  ! was parameter co2=3.75e-4,
    endif

#if defined (with_budgets)
    l_sum_wet = 0.0
#endif

    do j = j1, j2
    do i = i1, i2
          cloudy(i,j)=.false.

          ! lwc is dimensionless 
          if ((ye(i,j,ilwc) > 1e-10).and.(ye(i,j,icc) > 0.02)) then 
             cloudy(i,j)=.true. 
             TEMP=ye(i,j,i_temp)
             ZTR=(1./TEMP-1./298)
             RT=TEMP*rg
             ITEMP=nint(TEMP-float(ntlow))
             !
             !CEV sulph is the initial total sulphate content in accumulation+
             !coarse mode, the incloud production is calculated on bulk
             !characteristics, and redistributed on the modes depending on their
             !particle numbers
#ifdef with_m7
             !Stelios: small contributions from nucleation and Aitken mode
             !         as well as gas-phase should be added for pH calculation
             sulph(i,j)=y0(i,j,iso4acs)+y0(i,j,iso4cos)+&
                        y0(i,j,iso4nus)+y0(i,j,iso4ais)+y0(i,j,iso4)
#else
             sulph(i,j)=y0(i,j,iso4)
#endif
             !
             ! Henry and dissociation equilibria
             !
             dkh2o(i,j) =1.01e-14*exp(-6706.0 *ztr)   !h2o<=>hplus+so3--
             !bug hkco2(i,j) =3.4e-2*(2420.*ztr)           ! is already dimensionless
             hkco2(i,j) =3.4e-2*exp(2420.*ztr)           ! is already dimensionless
             dkco2(i,j) =4.5E-7*exp(-1000.*ztr)       !co2aq<=>hco3- + hplus
             hkso2(i,j) =henry(iso2,itemp)*rt         !dimensionless
             dknh3(i,j) =1.8e-5*exp(-450.*ztr)        !nh3<=>nh4+ + OH-
             hknh3(i,j) =henry(inh3,itemp)*rt         !dimensionless
             hkh2o2(i,j)=henry(ih2o2,itemp)*rt        !dimensionless
             hko3(i,j)  =henry(io3,itemp)*rt          !dimensionless
             dkso2(i,j) =1.7e-2*exp(2090.*ztr)        !so2<=>hso3m+hplus
             dkhso3(i,j)=6.6e-8*exp(1510.*ztr)        !hso3m<=>so3-- + hplus
             !
             ! calculate H+ from initial sulfate, ammonium, hno3, and nh3
             ! if solution is strongly acidic no further calculations are performed
             !

             xl=ye(i,j,ilwc)*Avog*1e-3/ye(i,j,icc)
             !x1 is initial strong acidity from SO4 and NO3
             !
             !acidity from strong acids alone
             !  
             !CMK hplus(i,j)=(2.*y0(i,j,iso4)+y0(i,j,imsa)-y0(i,j,inh4)+ &
             !CMK      y0(i,j,ihno3)+y0(i,j,ino3_a))/xl
             hplus(i,j)=(2.*sulph(i,j) + &
                  y0(i,j,imsa)-y0(i,j,inh4)+ &
                  y0(i,j,ihno3)+y0(i,j,ino3_a))/xl
          end if
    end do
    end do

    do iter=1,10
       do j=j1, j2
       do i=i1, i2
             ! only if solution pH>4.5
             if ( cloudy(i,j) .and. hplus(i,j) < 3e-5 ) then
                xl=ye(i,j,ilwc)*Avog*1e-3/ye(i,j,icc)
                !CEV  y0(i,j,iso4)---> sulph(i,j)
                x1=(2.*sulph(i,j)+y0(i,j,imsa)+y0(i,j,ihno3)+ &
                     y0(i,j,ino3_a))/xl    
                !x2 is initial total NHx
                x2=(y0(i,j,inh3)+y0(i,j,inh4))/xl
                !x3 is combined dissolution and solubility const for CO2
                x3=dkco2(i,j)*hkco2(i,j)*co2 
                a1=dkh2o(i,j)/dknh3(i,j)*(1.+1./hknh3(i,j)) ! integration constant
                a2=y0(i,j,iso2)/xl      !initial SO2
                ! trap division by zero ...
                if ( hplus(i,j) == 0.0 ) then
                   z = 0.0
                else
                   z = a2/( hplus(i,j)/dkso2(i,j)*(1.0+1.0/hkso2(i,j)) + dkhso3(i,j)/hplus(i,j) + 1.0 )
                end if
                ! solve quadratic equation for new H+ concentration:
                a=1.+x2/(a1+hplus(i,j))
                b=-x1-z
                c=-x3-2.*dkhso3(i,j)*z
                z=max(0.,(b*b-4.*a*c))
                hplus(i,j)=max(1.e-10,(-b+sqrt(z))/(2.*a))
             end if
       end do !
       end do ! i,j loop
    end do   !iter

    do j=j1,j2
    do i=i1,i2
          if (cloudy(i,j)) then
             temp=ye(i,j,i_temp)
             ZTR=(1./TEMP-1./298)
             xliq=ye(i,j,ilwc)/ye(i,j,icc)
             xl=ye(i,j,ilwc)*Avog*1e-3/ye(i,j,icc)
             ye(i,j,iph)=-log10(hplus(i,j))     ! pH for diagnostics 

             ! phase factor ratio of aqueous phase to gas phase concentration

             phs4   =hkso2(i,j) *(1.+dkso2(i,j)/hplus(i,j)+ &
                  dkso2(i,j)*dkhso3(i,j)/hplus(i,j)/hplus(i,j))*xliq
             phso2  =hkso2(i,j) *xliq
             phh2o2 =hkh2o2(i,j)*xliq
             phozone=hko3(i,j)  *xliq

             ! the original rate equations could be partly in look-up table

             rw(i,j,KSO2HP) =8e4*exp(-3560.*ztr)/(0.1+hplus(i,j))
             XSO2O3A=4.39e11*exp(-4131./temp)+2.56e3*exp(-966./temp)  !S(IV)
             XSO2O3B=2.56e3*exp(-966./temp)/hplus(i,j)  
             !divide by [H+]!S(IV)

             !  make rate constants dimensionless by multiplying 
             !  by (1./xliq/avo=6e20)  
             !  multiply with the fractions of concentrations residing 
             !  in the aqueous phase

             rw(i,j,KSO2HP)=rw(i,j,KSO2HP)/xl*phso2/(1.+phs4)*phh2o2/(1.+phh2o2)
             rw(i,j,KSO2O3)=(XSO2O3A+XSO2O3B)/xl*phs4/(1.+phs4)*phozone/ &
                  (1.+phozone)
          end if !cloudy
    end do ! 
    end do ! I,J, LOOP

    !-------------  Start main loop

    do j=j1,j2
    do i=i1,i2
          !
          ! only cloud chemistry if substantial amount of clouds are present 
          !
          if (cloudy(i,j)) then
             !
             ! oxidation of S(IV) by O3
             !
             so2x=y0(i,j,iso2)
             xcov=ye(i,j,icc)
             x1=min(100.,rw(i,j,kso2o3)*y0(i,j,io3)*dt)
             dso2=y0(i,j,iso2)*xcov*(exp(-x1)-1.)
             !only applied to xcov part of cloud
             !CMK         print *, i,j, xcov, x1, y0(i,j,iso2), dso2
             dso2=max(-y0(i,j,io3)*xcov,dso2)! limit to o3 availability
             !CEV     y(i,j,iso2)=y0(i,j,iso2)+dso2 
             !NOTE CMK: parallel MPI should take care here!
#ifdef with_m7
             ft = y0(i,j,iacs_n) + y0(i,j,icos_n)
             if (ft > tiny(ft)) then 
                y(i,j,iso4acs)=y0(i,j,iso4acs)-dso2*(y0(i,j,iacs_n)/ft)
                y(i,j,iso4cos)=y0(i,j,iso4cos)-dso2*(y0(i,j,icos_n)/ft)
                y(i,j,iso2)=y0(i,j,iso2)+dso2
                y(i,j,io3)=y0(i,j,io3)+dso2
                !AERHEMMIP
                ! Production of liquid phase so4
                ! D_prod_liq(i,j)=-dso2
#ifdef with_budgets
                c4(i,j,1)=c4(i,j,1)+dso2
                !conversion 1e-3 g->kg 1e6 cm-3 ->1m-3

                diag_prod(region)%prod(i,j,level,2)=diag_prod(region)%prod(i,j,level,2)-dso2/y(i,j,iair)*ye(i,j,iairm)/xmair*xmso4
                AC_diag_prod(region)%prod(i,j,level,2)=AC_diag_prod(region)%prod(i,j,level,2)-dso2/y(i,j,iair)*ye(i,j,iairm)/xmair*xmso4
#endif
             else
#ifdef with_budgets
                diag_prod(region)%prod(i,j,level,2)=diag_prod(region)%prod(i,j,level,2)+0.0
                AC_diag_prod(region)%prod(i,j,level,2)=AC_diag_prod(region)%prod(i,j,level,2)+0.0

#endif

                !CEV            y(i,j,iso4)=y0(i,j,iso4)-dso2
                y(i,j,iso4acs)=y0(i,j,iso4acs)
                y(i,j,iso4cos)=y0(i,j,iso4cos)
                y(i,j,iso2)=y0(i,j,iso2)
                y(i,j,io3)=y0(i,j,io3)
             endif
             !CEV          y(i,j,io3)=y0(i,j,io3)+dso2
#else
             ! gas phase chemistry:    ft = 0.
             y(i,j,iso4)=y0(i,j,iso4)-dso2
             y(i,j,iso2)=y0(i,j,iso2)+dso2
             y(i,j,io3)=y0(i,j,io3)+dso2
#ifdef with_budgets
             c4(i,j,1)=c4(i,j,1)+dso2
#endif
#endif

#ifdef with_budgets
             if ( io3 == 1 ) l_sum_wet = l_sum_wet- &
                  dso2 *ye(i,j,iairm)/ (fscale(1)*y(i,j,iair))
             if ( iso2 == 1 ) l_sum_wet = l_sum_wet+ &
                  dso2   *ye(i,j,iairm)/ (fscale(1)*y(i,j,iair))
             if ( iso4 == 1 ) l_sum_wet = l_sum_wet- &
                  dso2  *ye(i,j,iairm)/ (fscale(1)*y(i,j,iair))
             !CEV           c4(i,j,1)=c4(i,j,1)+dso2
#endif
             xliq=ye(i,j,ilwc)/ye(i,j,icc)
             !
             ! oxidation of S(IV) by H2O2
             !
             !*** here we explicitly solve the dv: 
             !    y'=P-Q*y-R*y*y (P and Q are 0=>b3=0.)
             !
             so2x=y(i,j,iso2)
             if ( so2x > tiny(so2x) ) then
                b1=rw(i,j,kso2hp)
                b2=b1*(y0(i,j,ih2o2)-so2x)
                disc=min(100.,sqrt(b2*b2))                ! disc is b2 for b3=0.0
                x1=(b2-disc)/(-2.*b1)                     ! in this case x1 =0.
                x2=(b2+disc)/(-2.*b1)
                x3=(so2x-x1)/(so2x-x2)*exp(-disc*dt)
                so2x=(x1-x2*x3)/(1.-x3)
                dso2=(so2x-y(i,j,iso2))*xcov
                dso2=max(dso2,-y0(i,j,ih2o2)*xcov)
                !CEV         y(i,j,iso2) =y(i,j,iso2)+dso2       ! dso2 is loss of SO2 and H2O2
                ! divide produced SO4 over CO/ACC 
#ifdef with_m7
                ft = y0(i,j,iacs_n) + y0(i,j,icos_n)
                if (ft > tiny(ft)) then 
                   y(i,j,iso4acs)=y(i,j,iso4acs)-dso2*(y0(i,j,iacs_n)/ft)
                   y(i,j,iso4cos)=y(i,j,iso4cos)-dso2*(y0(i,j,icos_n)/ft)
                   y(i,j,iso2)=y(i,j,iso2)+dso2
                   y(i,j,ih2o2)=y0(i,j,ih2o2)+dso2
#ifdef with_budgets
                   c4(i,j,2)=c4(i,j,2)+dso2
                   ! add amount liquid so4 production to diagnostic and change to molec->mass kg
                   !
                   diag_prod(region)%prod(i,j,level,2)=diag_prod(region)%prod(i,j,level,2)-dso2/y(i,j,iair)*ye(i,j,iairm)/xmair*xmso4
                   AC_diag_prod(region)%prod(i,j,level,2)=AC_diag_prod(region)%prod(i,j,level,2)-dso2/y(i,j,iair)*ye(i,j,iairm)/xmair*xmso4
                   !*1e6*1e-3/avog*xmso4
#endif
                else
#ifdef with_budgets
                   diag_prod(region)%prod(i,j,level,2)=diag_prod(region)%prod(i,j,level,2)+0.0
                   AC_diag_prod(region)%prod(i,j,level,2)=AC_diag_prod(region)%prod(i,j,level,2)+0.0
                   !*1e6*1e-3/avog*xmso4
#endif
                   y(i,j,ih2o2) =y0(i,j,ih2o2) 
                endif
#else
                ! gas - phase chemistry
                y(i,j,iso4)=y(i,j,iso4)-dso2
                y(i,j,iso2)=y(i,j,iso2)+dso2
                y(i,j,ih2o2)=y0(i,j,ih2o2)+dso2
#ifdef with_budgets
                c4(i,j,2)=c4(i,j,2)+dso2
#endif
#endif
#ifdef with_budgets
                if ( ih2o2 == 1 ) l_sum_wet = l_sum_wet- &
                     dso2  *ye(i,j,iairm)/ (fscale(1)*y(i,j,iair))
                if ( iso2 == 1 ) l_sum_wet = l_sum_wet- &
                     dso2   *ye(i,j,iairm)/ (fscale(1)*y(i,j,iair))
                if ( iso4 == 1 ) l_sum_wet = l_sum_wet+ &
                     dso2  *ye(i,j,iairm)/ (fscale(1)*y(i,j,iair))
                !CEV             c4(i,j,2)=c4(i,j,2)+dso2
#endif
             end if

             !
             ! NH3 uptake in cloud droplets is limited by SO4 availability
             ! no HNO3 is considered at this point
             ! assume instantaneous uptake of NH3 incloud  only in cloudy part
             !
             ! EQSAM gives hell to any NH3/NH4 interchange. It first drops both in one heap
             ! and then redistributes. So this cloud chemistry reaction does not matter.
          end if     !cloudy
    end do ! i,j,loop
    end do ! 

#ifdef with_budgets    
    sum_wetchem(region) = sum_wetchem(region) + l_sum_wet
#endif

    !free memory
    deallocate(hkso2  )
    deallocate(hkh2o2 )
    deallocate(hko3   )
    deallocate(dkso2  )
    deallocate(dkhso3 )
    deallocate(dkh2o  )
    deallocate(dknh3  )
    deallocate(hknh3  )
    deallocate(hkco2  )
    deallocate(dkco2  )
    deallocate(hplus  )
    deallocate(rw     )
    deallocate(cloudy )
    deallocate(sulph  )
    status = 0

  end subroutine WETS
  !EOC  
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    MARK_TRAC
  !
  ! !DESCRIPTION: calculate nox/pan/orgn/hno3 analogous to ebi scheme
  !               ozone production from marked nox
  !               simple nhx chemistry, scaled to real nhx
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine MARK_TRAC( region, level, is, js, y, rr, rj, dt, ye)
    !
    ! !USES:
    !
    use chem_param
    use Dims,           only : CheckShape
    use global_data,    only : region_dat
    use dims,           only : at, bt, im, jm
#ifdef with_budgets
    use budget_global,  only : budg_dat, nzon_vg
#endif
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in) :: region, level, is, js
    real                :: dt ! time step
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    real, intent(inout):: y (is:,js:,:) ! concentrations
    real, intent(in)   :: rr(is:,js:,:) ! reaction rates
    real, intent(in)   :: rj(is:,js:,:) ! photolysis rates
    real, intent(in)   :: ye(is:,js:,:) ! help fields ( air masses )
    !
    ! !REVISION HISTORY: 
    !      Jan 1999 - fjd
    !      Jan 2002 - MK          - adapted for TM5
    !   22 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer :: status, i1, i2, j1, j2

    ! --- begin --------------------------------
    
    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    
    ! check arguments ...
    call CheckShape( (/i2-i1+1, j2-j1+1, maxtrace/), shape(y ), status )
    call CheckShape( (/i2-i1+1, j2-j1+1, nreac   /), shape(rr), status )
    call CheckShape( (/i2-i1+1, j2-j1+1, nj      /), shape(rj), status )
    call CheckShape( (/i2-i1+1, j2-j1+1, n_extra /), shape(ye), status )

    call MARK_O3S
    !
    ! more marked tracers possible here
    !

  contains


    subroutine mark_o3s
      !---------------------------------------------------
      ! marked tracer  O3S stratospheric ozone
      !---------------------------------------------------

#ifndef without_dry_deposition
      use dry_deposition, only: vd
#endif

      integer :: i, j, nzone, nzone_v
      real    :: p3, xl3, o3old

      do j = j1, j2
      do i = i1, i2
            if (at(level+1)+bt(level+1)*1e5<= 14000) then ! 
               ! well, you want to count all layers below 140 hPa 
               ! (given surface pressure of 1e5 Pa)
               ! in the current model setup (25 layers) this means  
               ! 12077 + 1e5*0.00181 = 12258 Pa and above...

               ! p3: production of o3 in stratosphere
               P3 = rj(i,j,jano3)*y(i,j,ino3)+ &
                    rj(i,j,jno2)*y(i,j,ino2)
               XL3= rr(i,j,ko3ho2)*y(i,j,iho2)+&
                    rr(i,j,ko3oh)*y(i,j,ioh)+ &
                    rr(i,j,kno2o3)*y(i,j,ino2)+&
                    rj(i,j,jo3d)+&
                    rr(i,j,knoo3)*y(i,j,ino)+&
                    rr(i,j,kc62)*y(i,j,ieth)+&
                    rr(i,j,kc58)*y(i,j,iole)+&
                    rr(i,j,kc77)*y(i,j,iisop)+&
                    rr(i,j,ko3c3h6)*y(i,j,ic3h6)+&
                    rr(i,j,ko3terp)*y(i,j,iterp)+&
                    rr(i,j,ko3ispd)*y(i,j,iispd)
            else
               !
               ! these are only the net destruction reactions
               !
               P3 = 0.
               XL3= rr(i,j,ko3ho2)*y(i,j,iho2)+&
                    rr(i,j,ko3oh)*y(i,j,ioh)+&
                    rj(i,j,jo3d)+&
                    rr(i,j,kc62)*y(i,j,ieth)+&
                    rr(i,j,kc58)*y(i,j,iole)+&
                    rr(i,j,kc77)*y(i,j,iisop)+&
                    rr(i,j,ko3c3h6)*y(i,j,ic3h6)+&
                    rr(i,j,ko3terp)*y(i,j,iterp)+&
                    rr(i,j,ko3ispd)*y(i,j,iispd)+&
                    rr(i,j,knh2o3)*y(i,j,inh2)+&
                    rr(i,j,knh2o2o3)*y(i,j,inh2o2)
               ! add up deposition....
#ifndef without_dry_deposition
               if ( level == 1 ) &
                    XL3 = XL3 + vd(region,io3)%surf(i,j)/ye(i,j,idz)
#endif
            end if
            o3old=y(i,j,io3s)
            y(i,j,io3s)=(o3old+p3*dt)/(1.+xl3*dt)
#ifdef with_budgets
            nzone=budg_dat(region)%nzong(i,j)    ! global budget
            nzone_v=nzon_vg(level) 
            budmarkg(nzone,nzone_v,1)=budmarkg(nzone,nzone_v,1)+ &
                 (y(i,j,io3s)-o3old)*ye(i,j,iairm)*1000./xmair/y(i,j,iair)
#endif
      end do
      end do !i,j

    end subroutine MARK_O3S

  end subroutine MARK_TRAC
  !EOC

end module EBISCHEME