ecearth3/sources/tm5mp/base/ebischeme__dummy.F90

!#################################################################
!
! Eulerian backward Iteration
! Chemistry solver for the CBM4 scheme 
!
!### macro's #####################################################
!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') 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"
!
!#################################################################

module ebischeme
!
!  use GO, only : gol, goPr, goErr
!  
!  implicit none
!
!  private
!
!  public :: ebi
!
!
!contains
!
!
!  subroutine ebi(region,level,rj,rr,y,ye)
!    !
!    ! perform Eulerian backwards chemistry at one model layer level in region 
!    ! input: 
!    ! region : region 
!    ! level:   model layer...
!    ! rj    : array of photolysis rates
!    ! rr    : array of reaction rates
!    ! y     : vector of concentrations (to be returned...)
!    !
!#ifdef MPI
!    use mpi_const,only : lmar,lmloc,myid
!#endif
!    use chem_param     
!    use dims, only: isr, ier, jsr,jer,im ,jm, nchem, tref, okdebug
!    use global_data, only: region_dat
!    use toolbox, only: dumpfield
!    use dry_deposition, only: vd  
!    !type(emis_data),dimension(nregions,ntrace) :: vd
!
!    implicit none
!
!    ! input/output
!    integer, intent(in)                                   :: region,level
!    real,dimension(im(region),jm(region),nreac),intent(in):: rr
!    real,dimension(im(region),jm(region),nj   ),intent(in):: rj
!    real,dimension(im(region),jm(region),maxtrace)        :: y
!    real,dimension(im(region),jm(region),n_extra)         :: ye
!
!    ! local
!    integer,dimension(:,:),pointer        :: zoomed
!    real,dimension(:,:,:),allocatable     :: cr2,cr3,cr4   !reaction budgets...
!    real,dimension(:,:,:),allocatable     :: y0
!
!    real    :: dtime
!    integer :: iterebi,i,j,ib,maxit,iter
!    integer :: offsetl
!
!    real :: R57,R89,P1,R12,R21,XL1,P2,XL2,P3,XL3,X1,X2,X3
!    real :: C1,C2,C3,Y2,XJT,R21T,R12T,R12TC,R21TC,XJTC,ACUB,BCUB,CCUB
!    real :: CUBDET,DNO2,R56,R65,R75,P5,XL5,R66,X5,P6,XL6,X6,C6,XL7
!    real :: Y1,C7,R98,P8,XL8,X4,C5,XL9,R1920,R1919,P19,XL19,R2019,XL20
!    real :: XLHNO3,PH2O2,XLH2O2,PCH2O,PCO,PHNO3,XLCH2O,PCH3O2,XLCH3O2
!    real :: PCH3O2H,XLCH3O2H,PALD2,XLALD2,PMGLY,XLMGLY,POLE,XLETH
!    real :: XLOLE,XLISOP,PRXPAR,XLRXPAR,PPAR,XLPAR,PROR,XLROR,PXO2
!    real :: XLXO2,PXO2N,XLXO2N,PROOH,XLROOH,PORGNTR,XLORGNTR,XLCO
!
!    real :: dt,dt2
!    real :: qdms,pso2,qso2,qso2d,qnh3,pnh2,qnh2,qdms1,qdms2,pmsa,fnoy
!
!    integer :: io,sfstart,sfend
!
!    ! 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
!
!    ! start
!
!    zoomed => region_dat(region)%zoomed
!
!    offsetl=0
!#ifdef MPI
!    if(myid>0) offsetl=sum(lmar(0:myid-1))
!#endif
!    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
!
!    if ( okdebug ) print *, 'EBI: called with:',region,offsetl+level, &
!         'with dt:',dt,iterebi
!
!    !cmkif(level==1) then
!    !cmkio = sfstart('ebi.hdf',dfacc_create)
!    !cmkcall io_write3d_32(io,im(region),'im',jm(region),'jm', &
!    !  nreac,'nreac',rr,'rr')
!    !cmkcall io_write3d_32(io,im(region),'im',jm(region),'jm', &
!    !  nj   ,'nj   ',rj,'rj')
!    !cmkcall io_write3d_32(io,im(region),'im',jm(region),'jm', &
!    !  maxtrace,'maxtrace',y,'y')
!    !cmkcall io_write3d_32(io,im(region),'im',jm(region),'jm', &
!    !  n_extra,'n_extra',ye,'ye')
!    !cmkio = sfend(io)
!    !cmkend if
!
!    allocate(y0(im(region),jm(region),maxtrace))
!    allocate(cr2(im(region),jm(region),nj   ))
!    allocate(cr3(im(region),jm(region),nreac))
!    allocate(cr4(im(region),jm(region),nreacw))
!
!    !*** SCALING OF NOx, which has changed due to transport/deposition
!    do j=jsr(region),jer(region)
!       do i=isr(region),ier(region)
!          if(zoomed(i,j)/=region) cycle
!          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
!    !
!    cr2=0.
!    cr3=0.
!    cr4=0.
!    !===========================================================
!    ! ** Start iterating over CHEMISTRY
!    !===========================================================
!    do ib=1,iterebi  
!       maxit=8   !CMKTEMP
!       if(offsetl+level<=3) maxit = maxit*2   ! lowest layers more iterations
!
!       y0 = y
!       !-------------------------------
!       ! wet sulphur/ammonia chemistry
!       !------------------------------
!       call wetS(region,level,y0,dt,y,ye,cr4)
!       !-------------------------------------
!       ! gasphase chemistry using EBI solver
!       !-------------------------------------
!
!       !cmk NOTE this statement was  there before but someone (me?) removed it.
!       y0 = y   
!
!       !cmk ______do EBI solver_______
!
!       !if (level  ==  19) then
!       !call dumpfield(0,'dumpb.hdf',rr,'rr')
!       !call dumpfield(1,'dumpb.hdf',rj,'rj')
!       !call dumpfield(1,'dumpb.hdf',y,'y')
!       !call dumpfield(1,'dumpb.hdf',ye,'ye')
!       !end if
!
!       dt2 = dt*dt
!       ! block the input for EBI for efficiency
!       ! copy values with faster running index in inside loop
!       iv = 0
!       do j=jsr(region),jer(region)
!          do i=isr(region),ier(region)
!             if(zoomed(i,j)/=region) cycle
!             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 ....
!                if(offsetl + 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
!                ! copy nox emissions....
!                do ivt=1,npts
!                   emino(ivt) = ye(ipts(ivt),jpts(ivt),ieno)
!                end do
!                ! do ebi solver....
!
!                call do_ebi(npts)    !the actual EBI solver on vectorized arrays  
!
!                do itrace=1,maxtrace
!                   do ivt=1,npts
!                      y(ipts(ivt),jpts(ivt),itrace)=yv(ivt,itrace)
!                   end do
!                end do
!                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 ....
!          if(offsetl + 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
!          do ivt=1,iv
!             emino(ivt) = ye(ipts(ivt),jpts(ivt),ieno)
!          end do
!
!          call do_ebi(iv)    !the actual EBI solver on remaining cells  
!
!          do itrace=1,maxtrace
!             do ivt=1,iv
!                y(ipts(ivt),jpts(ivt),itrace)=yv(ivt,itrace)
!             end do
!          end do
!       end if
!
!       call NOymass
!
!       !-------------------------------------
!       ! marked tracers
!       ! apply after correction of nitrogen components
!       !----------------------------------------------
!
!       call mark_trac(region,level,y,rr,rj,dt,ye)
!
!       !------------------------------------------------------------
!       ! increase budget accumulators cr2 and cr3 (cr4 is done in wetS)
!       !------------------------------------------------------------
!
!       call incc2c3
!       !===========================================================
!       ! ** END iterating over CHEMISTRY
!       !===========================================================
!    end do   !iterebi
!
!    call reacbud   !add budgets for this timestep
!
!    deallocate(y0)
!    deallocate(cr2)
!    deallocate(cr3)
!    deallocate(cr4)
!
!    nullify(zoomed)
!
!  contains 
!
!    subroutine do_ebi(lvec)
!
!      integer, intent(in) :: lvec
!      integer             :: ivec
!
!      do ITER=1,MAXIT
!
!
!         do ivec=1,lvec
!            ! --- Short living compounds & groups
!            ! --- First group: NO NO2 O3
!            P1=rjv(ivec,jbno3)*yv(ivec,ino3)+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)
!            XL1=rrv(ivec,knono3)*yv(ivec,ino3)+rrv(ivec,kc81)*yv(ivec,ixo2n)
!            XL1 = XL1 + vdv(ivec,ino)
!            P2=rjv(ivec,jhno3)*yv(ivec,ihno3)+rjv(ivec,jn2o5)*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))+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)&
!                 +rrv(ivec,kc84)*yv(ivec,iorgntr)*yv(ivec,ioh)+rjv(ivec,jorgn)&
!                 *yv(ivec,iorgntr)+rjv(ivec,jpan)*yv(ivec,ipan)+0.1*rrv(ivec,kc78) *yv(ivec,iisop)*yv(ivec,ino3)
!            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)
!            XL2 = XL2 + vdv(ivec,ino2)
!            P3=rjv(ivec,jano3)*yv(ivec,ino3)
!            XL3=rrv(ivec,ko3ho2)*yv(ivec,iho2)+rrv(ivec,ko3oh)*yv(ivec,ioh)+rrv(ivec,kno2o3)&
!                 *yv(ivec,ino2)+rjv(ivec,jo3d)+rrv(ivec,kc58)*yv(ivec,iole)&
!                 +rrv(ivec,kc62)*yv(ivec,ieth)+rrv(ivec,kc77)*yv(ivec,iisop)
!            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
!            ! --- Oplossen onbekende 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+XJT*yv(ivec,ino2))/(C3+Y1*yv(ivec,ino))
!            ! --- Second group: yv(ivec,iho2) yv(ivec,ioh) yv(ivec,ihno4)
!            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)
!            P5=2.*rjv(ivec,jbch2o)*yv(ivec,ich2o)+rrv(ivec,kc46)*yv(ivec,ic2o3)*yv(ivec,ino)&
!                 +rrv(ivec,kmo2no)*yv(ivec,ich3o2)*yv(ivec,ino)+0.66*rrv(ivec,kmo2mo2)&
!                 *yv(ivec,ich3o2)*yv(ivec,ich3o2)+rjv(ivec,jmepe)*yv(ivec,ich3o2h)&
!                 +2.*rrv(ivec,kc49)*yv(ivec,ic2o3)*yv(ivec,ic2o3)+2.*rjv(ivec,j45)&
!                 *yv(ivec,iald2)+rjv(ivec,j74)*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)&
!                 +rrv(ivec,kc57)*yv(ivec,iole)*yv(ivec,ioh)+0.25*rrv(ivec,kc58)*yv(ivec,io3)&
!                 *yv(ivec,iole)+rrv(ivec,kc61)*yv(ivec,ieth)*yv(ivec,ioh)+0.26*rrv(ivec,kc62)&
!                 *yv(ivec,ieth)*yv(ivec,io3)+0.85*rrv(ivec,kc76)*yv(ivec,iisop)*yv(ivec,ioh)&
!                 +0.3*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)+rrv(ivec,kc41)*yv(ivec,ich2o)&
!                 *yv(ivec,ino3)+rjv(ivec,jorgn)*yv(ivec,iorgntr)+0.9*rrv(ivec,kc78)*yv(ivec,iisop)*yv(ivec,ino3)
!            XL5=rrv(ivec,kho2no)*yv(ivec,ino)+rrv(ivec,kno2ho2)*yv(ivec,ino2)&
!                 +rrv(ivec,ko3ho2)*yv(ivec,io3)+rrv(ivec,kmo2ho2)*yv(ivec,ich3o2)&
!                 +rrv(ivec,kho2oh)*yv(ivec,ioh)+rrv(ivec,kc82)*yv(ivec,ixo2)+rrv(ivec,kc85)*yv(ivec,ixo2n)
!            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.79*rrv(ivec,kc50)*yv(ivec,ic2o3)*yv(ivec,iho2)+0.4*rrv(ivec,kc58)&
!                 *yv(ivec,io3)*yv(ivec,iole)+0.28*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)
!            XL6=rrv(ivec,khno4oh)*yv(ivec,ihno4)+rrv(ivec,kho2oh)*yv(ivec,iho2)&
!                 +rrv(ivec,kno2oh)*yv(ivec,ino2)+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.7*rrv(ivec,kohrooh)*yv(ivec,irooh)&
!                 +rrv(ivec,kc84)*yv(ivec,iorgntr)+rrv(ivec,kh2oh)&
!                 +(rrv(ivec,kdmsoha)+rrv(ivec,kdmsohb)) *yv(ivec,idms) +rrv(ivec,knh3oh)*yv(ivec,inh3)  !sulfur
!
!            X6=y0v(ivec,ioh)+P6*DT
!            C6=1.+XL6*DT
!            R75=rrv(ivec,kno2ho2)*yv(ivec,ino2)
!            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,ihno4)=(y0v(ivec,ihno4)+R75*DT*yv(ivec,iho2))/C7
!            ! --- Third group: NO3 N2O5
!            R89=rjv(ivec,jn2o5)+rrv(ivec,kn2o5)
!            P8=rrv(ivec,kohhno3)*yv(ivec,ihno3)*yv(ivec,ioh)+rrv(ivec,kno2o3)*yv(ivec,ino2)*yv(ivec,io3)
!            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)
!            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,jn2o5)+rrv(ivec,kn2o5)+rrv(ivec,kn2o5aq)+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,jpan)
!            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,j74)&
!                 *yv(ivec,imgly)+0.15*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)
!            XL19=rrv(ivec,kc46)*yv(ivec,ino)+rrv(ivec,kc50)*yv(ivec,iho2)+rrv(ivec,kc47)*yv(ivec,ino2)
!            XL19 = XL19 + vdv(ivec,ic2o3)
!            R2019=rrv(ivec,kc47)*yv(ivec,ino2)
!            XL20=rrv(ivec,kc48)+rjv(ivec,jpan)
!            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)
!            ! --- CH4 chemistry (short living radicals)
!            PCH3O2=rrv(ivec,kch4oh)*yv(ivec,ich4)*yv(ivec,ioh)+0.7*rrv(ivec,kohmper)*yv(ivec,ioh)*yv(ivec,ich3o2h)
!            XLCH3O2=rrv(ivec,kmo2no)*yv(ivec,ino)+rrv(ivec,kmo2ho2)*yv(ivec,iho2)&
!                 +2*rrv(ivec,kmo2mo2)*yv(ivec,ich3o2)
!            yv(ivec,ich3o2)=(y0v(ivec,ich3o2)+PCH3O2*DT)/(1.+XLCH3O2*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)+rrv(ivec,kc57)*yv(ivec,iole)*yv(ivec,ioh)+0.9*rrv(ivec,kc58)&
!                 *yv(ivec,io3)*yv(ivec,iole)+rrv(ivec,kc59)*yv(ivec,iole)*yv(ivec,ino3)
!            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)
!            PXO2=rrv(ivec,kc46)*yv(ivec,ic2o3)*yv(ivec,ino)+2.*rrv(ivec,kc49)&
!                 *yv(ivec,ic2o3)*yv(ivec,ic2o3)+rrv(ivec,kc50)*yv(ivec,ic2o3)&
!                 *yv(ivec,iho2)+rjv(ivec,j45)*yv(ivec,iald2)+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)+rrv(ivec,kc57)*yv(ivec,iole)*yv(ivec,ioh)&
!                 +0.29*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.85*rrv(ivec,kc76)*yv(ivec,iisop)*yv(ivec,ioh)+0.7*rrv(ivec,kohrooh)&
!                 *yv(ivec,irooh)*yv(ivec,ioh)+rrv(ivec,kc84)*yv(ivec,ioh)*yv(ivec,iorgntr)&
!                 +0.18*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)
!            XLXO2=rrv(ivec,kc79)*yv(ivec,ino)+2.*rrv(ivec,kc80)*yv(ivec,ixo2)+rrv(ivec,kc82)*yv(ivec,iho2)
!            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.15*rrv(ivec,kc76)*yv(ivec,iisop)*yv(ivec,ioh)
!            XLXO2N=rrv(ivec,kc81)*yv(ivec,ino)+rrv(ivec,kc85)*yv(ivec,iho2)
!            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,kn2o5aq)+rrv(ivec,kn2o5l))&
!                    *yv(ivec,in2o5)+rrv(ivec,kc44)*yv(ivec,iald2)*yv(ivec,ino3)+rrv(ivec,kc41)*yv(ivec,ich2o)*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)
!               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,kmo2ho2)*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=rrv(ivec,kc46)*yv(ivec,ic2o3)*yv(ivec,ino)+0.3*rrv(ivec,kohmper)&
!                    *yv(ivec,ich3o2h)*yv(ivec,ioh)+rrv(ivec,kmo2no)*yv(ivec,ich3o2)*yv(ivec,ino)&
!                    +1.33*rrv(ivec,kmo2mo2)*yv(ivec,ich3o2)*yv(ivec,ich3o2)+rjv(ivec,jmepe)&
!                    *yv(ivec,ich3o2h)+2.*rrv(ivec,kc49)*yv(ivec,ic2o3)*yv(ivec,ic2o3)&
!                    +rrv(ivec,kc50)*yv(ivec,ic2o3)*yv(ivec,iho2)+rjv(ivec,j45)*yv(ivec,iald2)&
!                    +rrv(ivec,kc57)*yv(ivec,iole)*yv(ivec,ioh)+0.64*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.61*rrv(ivec,kc76)*yv(ivec,iisop)*yv(ivec,ioh)+0.9*rrv(ivec,kc77)&
!                    *yv(ivec,iisop)*yv(ivec,io3)+0.03*rrv(ivec,kc78)*yv(ivec,iisop)*yv(ivec,ino3)
!               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)+rrv(ivec,kc57)*yv(ivec,iole)*yv(ivec,ioh)+0.44*rrv(ivec,kc58)&
!                    *yv(ivec,iole)*yv(ivec,io3)+rrv(ivec,kc59)*yv(ivec,iole)*yv(ivec,ino3)&
!                    +0.22*rrv(ivec,kc61)*yv(ivec,ieth)*yv(ivec,ioh)+0.12*rrv(ivec,kc78)*yv(ivec,iisop)*yv(ivec,ino3)
!               XLALD2=rrv(ivec,kc43)*yv(ivec,ioh)+rrv(ivec,kc44)*yv(ivec,ino3)+rjv(ivec,j45)
!               XLALD2=XLALD2 + vdv(ivec,iald2)
!               yv(ivec,iald2)=(y0v(ivec,iald2)+PALD2*DT)/(1.+XLALD2*DT)
!               PMGLY=0.03*rrv(ivec,kc76)*yv(ivec,iisop)*yv(ivec,ioh)+0.03*rrv(ivec,kc77)&
!                    *yv(ivec,iisop)*yv(ivec,io3)+0.08*rrv(ivec,kc78)*yv(ivec,iisop)*yv(ivec,ino3)
!               XLMGLY=rrv(ivec,kc73)*yv(ivec,ioh)+rjv(ivec,j74)
!               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.58*rrv(ivec,kc76)*yv(ivec,iisop)*yv(ivec,ioh)+0.55*rrv(ivec,kc77)&
!                    *yv(ivec,iisop)*yv(ivec,io3)+0.45*rrv(ivec,kc78)*yv(ivec,iisop) *yv(ivec,ino3)
!               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)+0.21*rrv(ivec,kc50)&
!                    *yv(ivec,ic2o3)*yv(ivec,iho2)+rrv(ivec,kc85)*yv(ivec,iho2)*yv(ivec,ixo2n)
!               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.9*rrv(ivec,kc78)*yv(ivec,iisop)*yv(ivec,ino3)
!               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)
!
!               ! 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)  
!               yv(ivec,iso4)=(y0v(ivec,iso4)+qso2*yv(ivec,iso2)*DT) /(1. + vdv(ivec,iso4)*DT)   !corredted CMK qso2/qso2d
!               yv(ivec,iacid)=(y0v(ivec,iacid)+(pmsa+2.*qso2*yv(ivec,iso2))*DT)/&
!                    (1.+rrv(ivec,knh3SO4)*yv(ivec,inh3)*DT)
!               yv(ivec,inh4)=(y0v(ivec,inh4)+yv(ivec,iacid)*rrv(ivec,knh3SO4)*yv(ivec,inh3)*DT)/(1.+vdv(ivec,inh4)*DT)
!               pnh2=yv(ivec,ioh)*rrv(ivec,knh3oh)
!               qnh3=yv(ivec,iacid)*rrv(ivec,knh3SO4)+pnh2
!               qnh3 = qnh3 + vdv(ivec,inh3)
!               yv(ivec,inh3)=y0v(ivec,inh3)/(1.+qnh3*DT)
!               qnh2= rrv(ivec,knh2no)*yv(ivec,ino)+rrv(ivec,knh2no2)*yv(ivec,ino2)&
!                    +rrv(ivec,knh2ho2)*yv(ivec,iho2) +rrv(ivec,knh2o2)+rrv(ivec,knh2o3)*yv(ivec,io3)
!               yv(ivec,inh2)=(y0v(ivec,inh2)+yv(ivec,inh3)*pnh2*dt)/(1.+qnh2*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)
!               PCO=yv(ivec,ich2o)*(rjv(ivec,jach2o)+rjv(ivec,jbch2o)+yv(ivec,ioh)&
!                    *rrv(ivec,kfrmoh))+rjv(ivec,j45)*yv(ivec,iald2)+rjv(ivec,j74)*yv(ivec,imgly)&
!                    +0.37*rrv(ivec,kc58)*yv(ivec,iole)*yv(ivec,io3)+0.43*rrv(ivec,kc62)&
!                    *yv(ivec,ieth)*yv(ivec,io3)+0.36*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)&
!                    +rrv(ivec,kc41)*yv(ivec,ich2o)*yv(ivec,ino3)
!               XLCO = rrv(ivec,kcooh)*yv(ivec,ioh)
!               XLCO = XLCO + vdv(ivec,ico)
!               yv(ivec,ico)=(y0v(ivec,ico)+PCO*DT)/(1.+XLCO*DT)
!               PPAR=0.63*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)+0.63*rrv(ivec,kc76)*yv(ivec,iisop)*yv(ivec,ioh)
!               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)
!               !if(yv(ivec,ipb210) < 0.0 .or. yv(ivec,irn222) < 0.0) then
!               !   print *, 'Negatives .....rn222, pb210', y0v(ivec,irn222), yv(ivec,irn222) , y0v(ivec,ipb210),  yv(ivec,ipb210)
!               !end if
!            end do   !ivec
!
!         end if
!
!      end do !ITER
!
!    end subroutine do_ebi
!
!
!    subroutine NOYmass
!#ifdef MPI    
!      use mpi_comm,only: stopmpi
!#endif
!      implicit none
!      integer i,j,imax, offsetl
!      real :: ncormax,ncorav,totn,totn0,fnoy,fnoy1
!      real :: ncorr,ncorr1,ncorr2,ncorr3, totdep
!      logical :: nerr
!
!      offsetl=0
!#ifdef MPI
!      if(myid>0) offsetl=sum(lmar(0:myid-1))
!#endif
!
!      ncormax=0.
!      ncorav=0.
!      nerr=.false.
!      imax = 0
!      do j=jsr(region),jer(region)
!         do i=isr(region),ier(region)
!            if(zoomed(i,j)/=region) cycle
!            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+offsetl == 1) then
!               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)  + &
!                    2*y(i,j,in2o5)*vd(region,in2o5)%surf(i,j)  + &
!                    y(i,j,ihno4)*vd(region,ihno4)%surf(i,j) )*dt/ye(i,j,idz)  
!            else
!               totdep  = 0.0 
!            end if
!            totn0=y0(i,j,inox)+y0(i,j,ihno3)+y0(i,j,ipan)+ &
!                 y0(i,j,iorgntr) + 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)
!            ! 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.
!               print *,'NOYmass: N-error....',region,offsetl+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
!
!            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 
!            ! (as a group of reservoir tracers)
!            ! 
!            totn=y(i,j,ihno3)+y(i,j,ipan)+y(i,j,iorgntr)
!            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)
!            y(i,j,ihno3)  =max(0.,ncorr1)
!            y(i,j,ipan)   =max(0.,ncorr2)
!            y(i,j,iorgntr)=max(0.,ncorr3)
!            ncorr=ncorr1+ncorr2+ncorr3-y(i,j,ihno3)-y(i,j,ipan)- y(i,j,iorgntr)
!            !
!            ! 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 ) print*,'NOYmass: N-mass balance error, ncorr>5% ',&
!           'Maximum correction:',  ncormax,& 
!           'Average correction in this loop:',  ncorav/imax, imax, '(imax)'
!
!    end subroutine NOYmass
!
!
!
!    subroutine incc2c3
!      !
!      use budget_global,only : buddep_dat, sum_deposition, sum_chemistry
!      use budget_fine,only: depdry
!      implicit none
!      integer :: i1,n1,n2,jl,i,j,offsetl
!      ! nrj and nrr used for reaction budget calculations
!      integer,dimension(nj   ),parameter    ::  nrj=(/io3,ino2,ih2o2,ihno3,ihno4,in2o5,ich2o,ich2o, &
!           ich3o2h,ino3,ino3,ipan,iorgntr,iald2,imgly,irooh/)
!      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, iho2,iho2,in2o5,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,irn222, &
!           !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, ich3o2, ioh,iho2,0,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,ino,ino2,iho2,0,io3,0/),(/nreac,2/))
!
!      real :: c1,xdep
!
!      c1=dt*1000./xmair   !conversion to moles...
!
!      offsetl=0
!#ifdef MPI
!      if(myid>0) offsetl=sum(lmar(0:myid-1))
!#endif
!
!      ! reaction budgets
!
!      do i1=1,nj    !photolysis rates
!         n1=nrj(i1)
!         do j=jsr(region),jer(region)
!            do i=isr(region),ier(region)
!               if(zoomed(i,j)/=region) cycle
!               if(n1 > 0) cr2(i,j,i1)=cr2(i,j,i1)+rj(i,j,i1)*y(i,j,n1)
!            end do
!         end do
!      end do!i1=1,nj   
!      !
!      do i1=1,nreac !reactions
!         n1=nrr(i1,1) !make sure n1 > 0
!         n2=nrr(i1,2)
!         if (n2 > 0.) then
!            do j=jsr(region),jer(region)
!               do i=isr(region),ier(region)
!                  if(zoomed(i,j)/=region) cycle
!                  cr3(i,j,i1)= cr3(i,j,i1)+y(i,j,n1)*y(i,j,n2)*rr(i,j,i1)
!               end do
!            end do
!         else
!            do j=jsr(region),jer(region)
!               do i=isr(region),ier(region)
!                  if(zoomed(i,j)/=region) cycle
!                  cr3(i,j,i1)= cr3(i,j,i1)+y(i,j,n1)*rr(i,j,i1)
!               end do
!            end do
!         end if
!      end do  !i1=1,nreac
!
!      if ( offsetl + level == 1 ) then   ! deposition budget
!
!         do i1=1,ntrace
!            do j=jsr(region),jer(region)
!               do i=isr(region),ier(region)
!                  if(zoomed(i,j)/=region) cycle
!                  xdep = y(i,j,i1)*vd(region,i1)%surf(i,j)/ye(i,j,idz)* &
!                       c1*ye(i,j,iairm)/y(i,j,iair)  !from updated concentrations
!                  buddep_dat(region)%dry(i,j,i1) = &
!                       buddep_dat(region)%dry(i,j,i1) + xdep
!                  if ( region == nregions ) then
!                     depdry(i,j,i1) = depdry(i,j,i1) + xdep   !in mole
!                  end if
!                  if ( i1 == 1 ) then   !seperate deposition from 'other' chemistry
!                     sum_deposition(region) = sum_deposition(region) - &
!                          xdep*ra(1)*1e-3   ! in kg
!                     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
!               end do
!            end do
!         end do   !i1
!      else   ! other layers
!         do j=jsr(region),jer(region)
!            do i=isr(region),ier(region)
!               if(zoomed(i,j)/=region) cycle
!               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
!
!
!
!    subroutine reacbud
!      !------------------------------------------------------------------------
!      !
!      ! REACBUD increase reaction budgets for each reaction
!      !         arrays nrr and nrj determine which species are
!      !         involved in a reaction
!      !
!      !
!      !------------------------------------------------------------------------
!      use budget_global,only : budrjg,budrrg,budrwg,budg_dat,nzon_vg
!      use budget_fine,only: budrj,budrr,budrw,nzon,nzon_v
!      implicit none
!      integer :: i1,i,j,nzone,nzone_v
!      real    :: c1
!      !
!      ! attribute to regions
!      !
!      c1=dt*1000./xmair   !conversion to moles...
!      do j=jsr(region),jer(region)
!         do i=isr(region),ier(region)
!            if(zoomed(i,j)/=region) cycle
!            if(region==nregions) then    !finest region budget....
!               nzone=nzon(i,j)
!               nzone_v=nzon_v(offsetl+level)    !level is passed to ebi...
!               do i1=1,nj   
!                  budrj(nzone,nzone_v,i1)=budrj(nzone,nzone_v,i1)+ &
!                       cr2(i,j,i1)*c1*ye(i,j,iairm)/y(i,j,iair)   !units mole
!               end do !nj   
!               do i1=1,nreac
!                  budrr(nzone,nzone_v,i1)=budrr(nzone,nzone_v,i1)+ &
!                       cr3(i,j,i1)*c1*ye(i,j,iairm)/y(i,j,iair)   !units mole
!               end do
!               do i1=1,nreacw
!                  budrw(nzone,nzone_v,i1)=budrw(nzone,nzone_v,i1)- &
!                       cr4(i,j,i1)*(1000./xmair)*ye(i,j,iairm)/y(i,j,iair)
!                  !note: changed sign to get 'positive' budget, just a 
!                  !      matter of definition, !CMK
!               end do
!            end if
!            nzone=budg_dat(region)%nzong(i,j)    !global budget
!            nzone_v=nzon_vg(offsetl+level)    !level is passed to ebi...
!            do i1=1,nj   
!               budrjg(nzone,nzone_v,i1)=budrjg(nzone,nzone_v,i1)+ &
!                    cr2(i,j,i1)*c1*ye(i,j,iairm)/y(i,j,iair)   !units mole
!            end do !nj   
!            do i1=1,nreac
!               budrrg(nzone,nzone_v,i1)=budrrg(nzone,nzone_v,i1)+ &
!                    cr3(i,j,i1)*c1*ye(i,j,iairm)/y(i,j,iair)   !units mole
!            end do
!            do i1=1,nreacw
!               budrwg(nzone,nzone_v,i1)=budrwg(nzone,nzone_v,i1)- &
!                    cr4(i,j,i1)*(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
!
!    end subroutine REACBUD
!
!
!  end subroutine ebi
!
!
!
!  subroutine wetS(region,level,y0,dt,y,ye,c4)
!    !**********************************************************************
!    !     
!    !wetS - aqueous phase chemistry of sulfur  (and other)
!    !programmed by Ad Jeuken (KNMI) and Frank Dentener (IMAU)
!    !adapted for TM5 by Maarten Krol (IMAU) 1-2002
!    !
!    !purpose
!    !-------
!    !oxidation of SO2 and uptake of other gases in the aqueous phase
!    !
!    !interface
!    !---------
!    !call wetS(region,level,zoomed,y0,dt,y,ye,c4)
!    !region  region under operation (provides im,jm,lm via chemistry.mod)
!    !level   vertical level 
!    !zoomed       which cells to skip (chemistrty is done by child)
!    !y0    initial concentration
!    !dt    chemistry timestep
!    !yt    concentrations at time is t
!    !ye    extra fields (temperature, clouds, ....)
!    !c4    budget accumulatior
!    !   
!    !method
!    !------
!    !implicit solution of oxidation of SO2
!    !
!    !external
!    !--------
!    !none
!    !
!    !reference
!    !---------
!    !-
!    !**********************************************************************
!
!    use global_data,  only: region_dat
!    use reaction_data,only: nreacw,ntlow,kso2hp,kso2o3
!    use chem_param
!    use budget_global, only: sum_wetchem
!    use dims, only: isr, jsr, ier,jer, im, jm
!    use Binas, only: Avog
!
!    !    use toolbox, only: dumpfield, escape_tm
!    !    use dims, only: lm
!    implicit none
!
!    ! input/output
!
!    integer,intent(in)                                        :: region
!    integer,intent(in)                                        :: level
!    real,intent(in),dimension(im(region),jm(region),maxtrace) :: y0
!    real,intent(in)                                           :: dt 
!    real,intent(out),dimension(im(region),jm(region),maxtrace):: y 
!    real,dimension(im(region),jm(region),n_extra)             :: ye   !extra fields (temp, cc, pH)
!    real,dimension(im(region),jm(region),nreacw),intent(inout):: c4
!
!    ! local
!
!    integer,dimension(:,:),pointer                            :: zoomed
!    integer n,i,j,l,itemp,iter
!    real x1,x2,x3,b1,b2,so2x,dh2o2,dso2,disc,dnh3,dn2o5,xso2o3a,xso2o3b
!    real,parameter :: co2=3.20e-4,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 a1,a2,a,b,c,z              ! help variables
!    real xcov,xliq,xl,temp,rt,ztr,h2o,air,press ! meteo
!    real,dimension(:,:,:),allocatable   :: rw ! reaction rates
!    logical,dimension(:,:),allocatable ::  cloudy
!    !    character(len=2) :: levelc
!
!    ! start
!    
!    zoomed => region_dat(region)%zoomed
!
!    allocate(hkso2      (im(region),jm(region)))
!    allocate(hkh2o2     (im(region),jm(region)))
!    allocate(hko3       (im(region),jm(region)))
!    allocate(dkso2      (im(region),jm(region)))
!    allocate(dkhso3     (im(region),jm(region)))
!    allocate(dkh2o      (im(region),jm(region)))
!    allocate(dknh3      (im(region),jm(region)))
!    allocate(hknh3      (im(region),jm(region)))
!    allocate(hkco2      (im(region),jm(region)))
!    allocate(dkco2      (im(region),jm(region)))
!    allocate(hplus      (im(region),jm(region)))
!    allocate(rw (im(region),jm(region),nreacw))
!    allocate(cloudy (im(region),jm(region)))
!
!    !-----------------------------
!    ! wet phase reactions
!    !-----------------------------
!    rw   =0.0 
!    hplus=0.0
!
!    do j=jsr(region),jer(region)
!       do i=isr(region),ier(region)
!          if(zoomed(i,j)/=region) cycle
!          cloudy(i,j)=.false.
!
!          ! lwc is dimensionless 
!          if ((ye(i,j,ilwc) > 1e-10).and.(ye(i,j,icc) > 0.01)) then 
!             cloudy(i,j)=.true. 
!             TEMP=ye(i,j,i_temp)
!             ZTR=(1./TEMP-1./298)
!             RT=TEMP*rg
!             ITEMP=nint(TEMP-float(ntlow))
!             !
!             ! Henry and dissociation equilibria
!             !
!             dkh2o(i,j) =1.01e-14*exp(-6706.0 *ztr)   !h2o<=>hplus+so3--
!             hkco2(i,j) =3.4e-2*(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
!             !  
!             hplus(i,j)=(2.*y0(i,j,iso4)+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=jsr(region),jer(region)
!          do i=isr(region),ier(region)
!             if ( zoomed(i,j) /= region ) cycle
!             ! 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)
!                x1=(2.*y0(i,j,iso4)+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
!                z=a2/(hplus(i,j)/dkso2(i,j)*(1.+1./hkso2(i,j))+ &
!                     dkhso3(i,j)/hplus(i,j)+1.)
!                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=jsr(region),jer(region)
!       do i=isr(region),ier(region)
!          if(zoomed(i,j)/=region) cycle
!          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)+ &
!                  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
!!    write(levelc,'(i2.2)') level
!!    if(level == 1) then 
!!       call dumpfield(0,'rw.hdf',im(region),jm(region),nreacw,1,rw,'rw'//levelc)
!!       call dumpfield(1,'rw.hdf',im(region),jm(region),n_extra,1,ye,'ye'//levelc)
!!    else
!!       call dumpfield(1,'rw.hdf',im(region),jm(region),nreacw,1,rw,'rw'//levelc)
!!       call dumpfield(1,'rw.hdf',im(region),jm(region),n_extra,1,ye,'ye'//levelc)
!!    end if
!!    if(level == lm(1) ) call escape_tm(' forced stop ')
!
!    ! Start main loop
!    do j=jsr(region),jer(region)
!       do i=isr(region),ier(region)
!          if(zoomed(i,j)/=region) cycle
!          !
!          ! 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
!             y(i,j,iso2)=y0(i,j,iso2)+dso2 
!             !NOTE CMK: paralel MPI should take care here!
!             y(i,j,iso4)=y0(i,j,iso4)-dso2
!             y(i,j,io3)=y0(i,j,io3)+dso2
!             if ( io3 == 1 ) sum_wetchem(region) = sum_wetchem(region)- &
!                  dso2 *ye(i,j,iairm)/ (fscale(1)*y(i,j,iair))
!             if ( iso2 == 1 ) sum_wetchem(region) = sum_wetchem(region)+ &
!                  dso2   *ye(i,j,iairm)/ (fscale(1)*y(i,j,iair))
!             if ( iso4 == 1 ) sum_wetchem(region) = sum_wetchem(region)- &
!                  dso2  *ye(i,j,iairm)/ (fscale(1)*y(i,j,iair))
!             c4(i,j,1)=c4(i,j,1)+dso2
!             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)
!               y(i,j,iso2) =y(i,j,iso2)+dso2       ! dso2 is loss of SO2 and H2O2
!               y(i,j,iso4)=y(i,j,iso4)-dso2
!               y(i,j,ih2o2) =y0(i,j,ih2o2)+dso2 
!               if ( ih2o2 == 1 ) sum_wetchem(region) = sum_wetchem(region)- &
!                    dso2  *ye(i,j,iairm)/ (fscale(1)*y(i,j,iair))
!               if ( iso2 == 1 ) sum_wetchem(region) = sum_wetchem(region)- &
!                    dso2   *ye(i,j,iairm)/ (fscale(1)*y(i,j,iair))
!               if ( iso4 == 1 ) sum_wetchem(region) = sum_wetchem(region)+ &
!                    dso2  *ye(i,j,iairm)/ (fscale(1)*y(i,j,iair))
!               c4(i,j,2)=c4(i,j,2)+dso2
!             end if
!
!             !
!             ! NH3 uptake in cloud droplets is limited by H2SO4 availability
!             ! no HNO3 is considered at this point
!             ! assume instantaneous uptake of NH3 incloud  only in cloudy part
!             !
!             dnh3=max((2.*y(i,j,iso4)+y0(i,j,imsa)-y0(i,j,inh4))*xcov,0.)
!             dnh3=max(-y0(i,j,inh3)*xcov,-dnh3)
!             y(i,j,inh3)=y0(i,j,inh3)+dnh3                 ! dnh3 is loss of NH3  
!             y(i,j,inh4)=y0(i,j,inh4)-dnh3
!             if ( inh3 == 1 ) sum_wetchem(region) = sum_wetchem(region) - &
!                  dnh3*ye(i,j,iairm)/ (fscale(1)*y(i,j,iair))
!             if ( inh4 == 1 ) sum_wetchem(region) = sum_wetchem(region) + &
!                  dnh3*ye(i,j,iairm)/ (fscale(1)*y(i,j,iair))
!             c4(i,j,3)=c4(i,j,3)+dnh3
!          end if     !cloudy
!       end do ! i,j,loop
!    end do ! 
!
!    !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 )
!
!    nullify(zoomed)
!
!  !  write(levelc, '(i2.2)') level
!  !  call dumpfield(1,'wetS.hdf', im(region), jm(region), maxtrace, 1, y0, 'y0'//levelc)
!  !  call dumpfield(1,'wetS.hdf', im(region), jm(region), n_extra , 1, ye, 'ye'//levelc)
!  !  call dumpfield(1,'wetS.hdf', ntracet, ntemp, 1 , 1, henry, 'henry'//levelc)
!  !  call dumpfield(1,'wetS.hdf', im(region), jm(region), maxtrace, 1, y, 'y'//levelc)
!
!  end subroutine wets
!
!
!
!  subroutine mark_trac(region,level,y,rr,rj,dt,ye)
!    !      ---------
!    !   call subroutine mark_trac(region,level,zoomed,y,rr,rj,dt,ye)
!    !      region,level               :: where and which cells
!    !      zoomed                     :: which cells are done by child?
!    !      y    :: concentrations in layer
!    !      rr   :: reaction rates
!    !      rj   :: photolysis rates
!    !      dt   :: time step
!    !      ye   :: help fields ( air masses )
!    !         
!    !      method
!    !      ------
!    !      calculate nox/pan/orgn/hno3 analogous to ebi scheme
!    !      ozone production from marked nox
!    !      simple nhx chemistry, scaled to real nhx
!    !
!    !      fjd Mon Aug 10 16:55:35 MET 1998/Fri Jan  1 17:08:37 MET 1999
!    !      mk   adapted for TM5 jan/2002
!    !-------------------------------------------------------------
!    use global_data,   only : region_dat
!    use budget_global, only : budmarkg,budg_dat,nzon_vg
!    use budget_fine,   only : budmark,nzon,nzon_v
!    use chem_param
!    use dims, only: isr, ier, jsr, jer, at, bt, im, jm
!
!    implicit none
!
!    ! input/output
!    integer, intent(in)                            :: region,level
!    real,dimension(im(region),jm(region),maxtrace) :: y
!    real,dimension(im(region),jm(region),nreac),intent(in):: rr
!    real,dimension(im(region),jm(region),nj   ),intent(in):: rj
!    real                                           :: dt
!    real, dimension(im(region),jm(region),n_extra) :: ye
!
!    ! local
!    integer, dimension(:,:),pointer                :: zoomed
!
!    ! start
!
!    zoomed => region_dat(region)%zoomed
!
!    call mark_o3s
!    !
!    ! more marked tracers possible here
!    !
!    nullify(zoomed)
!
!  contains
!
!
!    subroutine mark_o3s
!      !---------------------------------------------------
!      ! marked tracer  O3S stratospheric ozone
!      !---------------------------------------------------
!
!#ifdef MPI
!      use mpi_const,only: myid,lmar
!#endif
!      use dry_deposition, only: vd
!      implicit none
!      integer :: i,j,nzone,nzone_v
!      real                           :: p3,xl3,o3old
!      integer :: offsetl
!
!      offsetl=0
!#ifdef MPI
!      if(myid>0) offsetl=sum(lmar(0:myid-1))
!#endif
!
!      do j=jsr(region),jer(region)
!         do i=isr(region),ier(region)
!            if(zoomed(i,j)/=region) cycle
!            if (at(offsetl+level+1)+bt(offsetl+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)
!            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)
!               ! add up deposition....
!               if ( offsetl + level == 1 ) &
!                    XL3 = XL3 + vd(region,io3)%surf(i,j)/ye(i,j,idz)
!            end if
!            o3old=y(i,j,io3s)
!            y(i,j,io3s)=(o3old+p3*dt)/(1.+xl3*dt)
!            if ( region == nregions ) then
!               nzone=nzon(i,j)
!               nzone_v=nzon_v(level+offsetl) 
!               ! budget in mole
!               budmark(nzone,nzone_v,1)=budmark(nzone,nzone_v,1)+ &
!                    (y(i,j,io3s)-o3old)*ye(i,j,iairm)*1000./xmair/y(i,j,iair)
!            end if
!            nzone=budg_dat(region)%nzong(i,j)    ! global budget
!            nzone_v=nzon_vg(level+offsetl) 
!            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)
!         end do
!
!      end do !i,j, l
!
!    end subroutine mark_o3s
!
!  end subroutine mark_trac
!
!
!
end module ebischeme