c=========================================================================
c  03/09/03: bkoo
c            - modified to eliminate organics (eqparto now deals with organics)
c            - moved call to wdiameter
c            - added call to eqparto
c  03/07/03: bkoo
c            - commented out last call to step (redundant)
c            - brought initial call to step out of IF statement
c              so that step can be called when HYBR is selected
c  12/05/02: tmg
c            - eliminated call to subroutine eqneut
c  10/22/02: tmg
c            - added parameter to subroutine step because not all sections
c              are always in use
c  01/25/02: tmg
c            - use dry and wet basis for diameter as appropriate
c  09/20/00: bkoo
c            - combine eqpart & eqparth
c  06/08/00: bkoo
c            - add coagulation & nucleation
c  Apr 2000: bkoo
c            - put a counter to avoid infinity loop
c            - fix negative concentrations
c            - need codes for organics evaporation
c=========================================================================
c
cgy supress error messages to units 39 and 90 (6/12/00)
c
c  EQILIBRIUM PARTITIONING OF CONDENSING OR EVAPORATING SPECIES
C  ACROSS A SIZE AND COMPOSITION RESOLVED AEROSOL DISTRIBUTION
C
C  THIS SUBROUTINE ASSUMES THAT GAS-AEROSOL EQUILIBRIUM IS ACHIEVED 
C  IN FOR THE SUM OF THE SIZE SECTIONAL COMPOSITIONS.
C
C  WRITTEN BY KEVIN P. CAPALDO
c  DEPARTMENT OF CHEMICAL ENGINEERING
C  CARNEGIE MELLON UNIVERSITY
C  JUNE 1998
C
c   SUBROUTINES REQUIRED
C     ISRPIA        EQUILIBRIUM THERMODYNAMIC SUBROUTINE ISORROPIA
C     ORGANICS      GIVES THE EQUILIBRIUM VAPOR PRESSURE OF ORGANIC SPECIES
C     STEP          CALLS ISRPIA IN REVERSE MODE FOR EACH SECTION
C     
C   VARIABLES
c     T             TIME IN SECONDS
C     Q(NTOTALx2)   ARRAY OF AEROSOL SPECIES FOR EACH SECTION (UG/M3) 
C                   THEN GASES (PPM)
C     DQDT(NTOTAL)  THE AEROSOL AND GAS DERIVATIVES (PER SECOND)
C     DQ(NSP)       THE TOTAL MASS OF EACH SPECIES TRANSFERED 
C                   FROM THE GAS TO THE AEROSOL PHASE (uMOLES/m3)
C          THE FOLLOWING TWO DIMENSIONAL ARRAYS ARE NECESARY FOR THE 
C          LINKING OF THE WEIGHTING FACTORS FOR NH4, NO3 AND CL.
C     FRQ(NSEC,nsp)   FRACTION OF DQ TRANSFERED TO EACH AEROSOL SECTION
C
      SUBROUTINE EQPART(t,q)

      include 'dynamic.inc'
      INCLUDE 'equaer.inc'                  ! ISORROPIA declarations

cbk      REAL*8 DQ(NSP),FRQ(nsec,nsp),WI(5),q0(3),q(ntotal)
cbk      real*8 qsav(ntotal), DQsav(nsp), accom(nsp)
      REAL*8 DQ(nexti),FRQ(nsec,nexti),WI(5),q0(3),q(ntotal) ! bkoo (03/09/03)
      real*8 qsav(nexti*nsec), DQsav(nexti), accom(nexti)    ! bkoo (03/09/03)
      real*8 qq, frq0(nsec,nexti) ! bkoo (10/07/03)
      logical done

      if(aerm.eq.'EQUI') then	
c
c     STEP 1/3: CALCULATE NUCLEATION RATE FOR THE WHOLE STEP
c
       call nucl(q)
c
c     STEP 2/3: CALCULATE COAGULATION RATE FOR THE WHOLE STEP
c
       call coagul(q)

cbk       call step(nsec,q) ! tmg (10/22/02)
      endif
      call step(nsecx2,q) ! bkoo (03/07/03)
      call wdiameter(q) ! bkoo (03/09/03)
      call eqparto(t,q) ! bkoo (03/09/03)
c
C     STEP 1:  DETERMINE BULK EQUILIBRIUM
c
C  INPUT:
c     IPROB = 1=Reverse prob, 0=Foreward prob
C  1. [WI] 
C     DOUBLE PRECISION array of length [NCOMP=5]. [NCOMP] is defined in 
C     include file 'isrpia.inc'.
C     Total aerosol concentrations, expressed in moles/m3.  
C     WI(1) - sodium,   expressed as [Na]
C     WI(2) - sulfate,  expressed as [H2SO4]
C     WI(3) - ammonium, expressed as [NH3]
C     WI(4) - nitrate,  expressed as [HNO3]
C     WI(5) - chloride, expressed as [HCl]
C
C  2. [RHI] 
C     DOUBLE PRECISION variable.  
C     Ambient relative humidity expressed on a (0,1) scale.
C
C  3. [TEMPI]
C     DOUBLE PRECISION variable. 
C     Ambient temperature expressed in Kelvins. 
C
C *** CONVERT INPUT CONCENTRATIONS TO moles/m3 **************************
C
      ng = nsp*nsecx2                       ! gases
      prod=rgas*temp/(1.01325D5*pres)       ! conversion from ppm to umoles/m3
                                            ! pres (bkoo, 06/09/00)
      WI(1) = 0.0D0
      WI(2) = Q(ng+ih2SO4) / PROD*1.0d-6
      WI(3) = q(ng+iNH3)   / PROD*1.0d-6
      WI(4) = q(ng+ihNO3)  / PROD*1.0d-6
      WI(5) = q(ng+ihCL)   / PROD*1.0d-6
      do k=1,nsecx2                         ! aerosols
         nn=(k-1)*nsp
         WI(1) = wi(1)+ q(nn+KNa) /emw(KNa) *1.D-6
         WI(2) = wi(2)+ q(nn+KSO4)/emw(KSO4)*1.D-6
         WI(3) = wi(3)+ q(nn+KNH4)/emw(KNH4)*1.D-6
         WI(4) = wi(4)+ q(nn+KNO3)/emw(KNO3)*1.D-6
         WI(5) = wi(5)+ q(nn+KCL) /emw(KCL) *1.D-6
      enddo
      RHI=rh
      TEMPI=temp
      IPROB = 0
C
C *** CALL ISORROPIA+ ***************************************************
C
      CALL ISRPIA ( WI, RHI, TEMPI, IPROB )

c     initialize DQ & ACCOM arrays - bkoo (05/24/01)
cbk      do i=1,nsp
      do i=1,nexti ! bkoo (03/09/03)
         dq(i)    = 0.0d0
         accom(i) = 1.0d0
      enddo
C
C   STEP 2:  DETERMINE THE TOTAL TRANSFER BETWEEN GAS AND AEROSOL
C
c      DQ(KH2O)=0.0D0                        ! WATER DONE LATER
c      DQ(KNA) =0.0D0                        ! SODIUM IS IN AEROSOL PHASE ONLY
      DQ(KSO4)=Q(ng+ih2so4)/prod            ! all H2SO4 condenses
      DQ(KNH4)=Q(NG+INH3)  /PROD - dmax1(GNH3 *1.0d6, 0.d0) ! bkoo (10/07/03)
      DQ(KNO3)=Q(NG+IHNO3) /PROD - dmax1(GHNO3*1.0d6, 0.d0) ! bkoo (10/07/03)
      DQ(KCL) =Q(NG+IHCL)  /PROD - dmax1(GHCL *1.0d6, 0.d0) ! bkoo (10/07/03)
C    DO ORGANICS  **ASSUME PS HAS NO SIZE OR COMPOSITION DEPENDANCE **
cbk  removed - bkoo (03/09/03)
cbk      DO IOG = 1,NORG-1
cbk         DQ(NEXTI+IOG)=(Q(NG+IHCL+IOG)-PS(1,IHCL+IOG)/
cbk     &                 (1.01325d-1*pres))/prod	! pres (bkoo, 06/09/00)
cbk      ENDDO
c
c   STEP 2B: PARTITION DQ's NH4, NO3, AND CL INTO NH4NO3, NH4CL, AND
C            EXCESS NH4 (which is associated with SO4)
C            DQ(KNO3) => now represents the transfer of NH4NO3
C            DQ(KCl) => now represents the transfer of NH4Cl
C            DQ(KNH4) => now represents the transfer of NH4 with SO4
c
      DQ(KNH4)=DQ(KNH4)-DQ(KNO3)-DQ(KCL)
c
c     Save initial aerosol concentrations
cbk      do i=1,ntotalx2
cbk        qsav(i)=q(i)
      do i=1,nsecx2                          ! bkoo (03/09/03)
       do ii=1,nexti                         ! bkoo (03/09/03)
        qsav((i-1)*nexti+ii)=q((i-1)*nsp+ii) ! bkoo (03/09/03)
       enddo                                 ! bkoo (03/09/03)
      enddo
c
c     Its possible to not be able to reach equilibrium for NH4Cl and NH4NO3
c     due to diferences between aerosol sectional vs. bulk composition 
      do i=1,3
        q0(i)=0.0d0
      enddo
      do isec=1,nsecx2
        q0(1)=q0(1)+q((isec-1)*nsp+kno3)    ! initial aerosol NO3
        q0(2)=q0(2)+q((isec-1)*nsp+knh4)    ! initial aerosol NH4
        q0(3)=q0(3)+q((isec-1)*nsp+kcl)     ! initial aerosol Cl
      enddo
C
c   STEP 3: DETERMINE THE RELATIVE RATES OF MASS TRANSFER FOR EACH SECTION
c       if we assume composition changes between size sections are not 
c       important, then the rate of transfer is proportional to the mass
c       mass transfer rate dependance on particle size

cbk   removed - bkoo (03/09/03)
cbk      accom(KH2O)=1.0
cbk      accom(KNa)=1.0
      accom(KSO4)=delta(ih2so4)
      accom(kno3)=delta(ihno3)
      accom(knh4)=delta(inh3)
      accom(kcl) =delta(ihcl)
cbk   removed - bkoo (03/09/03)
cbk      do isp=1,norg
cbk        accom(kcl+isp)=delta(ihcl+isp)
cbk      enddo
cbk      do isp=1,ninert
cbk        accom(kcl+norg+isp)=0.1
cbk      enddo

cbk      call wdiameter(q) ! tmg (01/25/02)

      rlambda=0.065d0
c
c calculate factors
c
cbk      do isp=1,nsp
      do isp=1,nexti ! bkoo (03/09/03)
       frqtot = 0.0
       do isec = 1,nsecx2
        frq(isec,isp) = qn(isec)*
     &      dsec(isec)/(1.0+rlambda/(accom(isp)*dsec(isec)))
        frqtot = frqtot + frq(isec,isp)
       enddo
c
c normalize
c
       do isec = 1,nsecx2
        frq(isec,isp) = frq(isec,isp)/frqtot
        frq0(isec,isp) = frq(isec,isp) ! save frq - bkoo (10/07/03)
       enddo
      enddo

      iter=0 ! counter for escaping out of an infinity loop (bkoo: Apr, 2000)
 80   iter=iter+1 ! moved - bkoo (10/07/03)
c     save dq's
cbk      do i=1,nsp
      do i=1,nexti ! bkoo (03/09/03)
        dqsav(i)=dq(i)
      enddo
c
c   FIRST condense all condensing species
c
cbk      do isp = 1,nsp
      do isp = 1,nexti ! bkoo (03/09/03)
       if(dq(isp).gt.0.0d0) then
        do isec=1,nsecx2
         INDX=(ISEC-1)*NSP
         q(indx+isp)=q(indx+isp)+frq(isec,isp)*dq(isp)*emw(isp)
c     special cases for remaping of dq
         if(isp.eq.kNO3.or.isp.eq.kCl) then
          q(indx+kNH4)= q(indx+knh4)+frq(isec,isp)*dq(isp)*emw(knh4)
         endif
        enddo
        dq(isp)=0.0d0
       endif
      enddo
c
c   SECOND evaporate all evaporating species
c   only NH4NO3 and NH4Cl can evaporate
c
C     CHECK FOR COMPLETE EVAPORATION
 100  frt=1.0d0
      frtcl = 0.0 ! bkoo (10/07/03)
      isp=kCl ! Cl first b/c NH4Cl forms before NH4NO3 when Na exists
      do m=1,2
       if(dq(isp).lt.0.0d0) then            ! evaporating
        do isec=1,nsecx2
         INDX=(ISEC-1)*NSP
         dqfx=DQ(ISP)*FRQ(ISEC,isp)
         IF(-dqfx.gt.tinys) then            ! evaporating significantly
          if(Q(INDX+ISP).LT.-dqfx*emw(isp)) then ! not enough NO3 or Cl
           frtq=-q(indx+isp)/dqfx/emw(isp)
           if(frtq.lt.frt) then
            frt=frtq
            ispsav=isp 	                    ! species that evaporates first
            isecsav=isec                    ! section that evaporates first
            if(q(indx+isp).lt.tinys) goto 150
           endif
          endif
          qq = q(indx+knh4) ! bkoo (10/07/03)
          if(isp.eq.KNO3)   ! bkoo (10/07/03)
     &              qq = qq + frtcl*frq(isec,KCL)*dq(KCL)*emw(knh4)
cbk          if(q(indx+knh4).lt.-dqfx*emw(kNH4)) then ! not enough NH4+
cbk           frtq=-q(indx+knh4)/dqfx/emw(kNH4)
          if(qq.lt.-dqfx*emw(kNH4)) then    ! not enough NH4+
           frtq=-qq/dqfx/emw(kNH4)
           if(frtq.lt.frt) then
            frt=frtq
            ispsav=isp                      ! species that evaporates first
            isecsav=isec                    ! section that evaporates first
            if(q(indx+knh4).lt.tinys) goto 150
           endif
          endif
         endif
        enddo
        if(isp.eq.KCL) frtcl = frt ! bkoo (10/07/03)
       endif
       isp=kNO3
      enddo
c partition species up to evaporation point
      done=.true.
      do m=1,2
       do isec=1,nsecx2
        INDX=(ISEC-1)*NSP
        q(indx+isp)=q(indx+isp)+frt*frq(isec,isp)*dq(isp)*emw(isp)
        q(indx+kNH4)=q(indx+knh4)+frt*frq(isec,isp)*dq(isp)*emw(knh4)
       enddo
       dq(isp)=(1-frt)*dq(isp)
       if(abs(dq(isp)).gt.1e-5) done=.false. ! tolerance for solution
       isp=kCl
      enddo
c check for complete solution
      if(done) then
       goto 200
      endif
c
c adjust frq to account for the totally evaporated species
c
 150  frqtot=0.0 ! bkoo (10/07/03) ...
      do isec=1,nsecx2
        if(isec.eq.isecsav) frq(isec,ispsav)=0.0d0
        frqtot = frqtot + frq(isec,ispsav)
      enddo
      if(frqtot.gt.tinys) then ! renormalize
        do isec=1,nsecx2
          frq(isec,ispsav)=frq(isec,ispsav)/frqtot
        enddo
      else
        if(iter.gt.itmaxeq) then ! moved - bkoo (10/07/03)
c         if(min(dq(KNO3),dq(KCL)).lt.-0.3) then                               ! bkoo_dbg
c         call get_param(igrdchm,ichm,jchm,kchm,iout,idiag)                    ! bkoo_dbg
c         write(*,'(A8,2E15.5,4I4)')'EQUI-F: ',dq(KNO3),dq(KCL),ichm,jchm,kchm ! bkoo_dbg
c         endif                                                                ! bkoo_dbg
         goto 200
        endif
        goto 180 ! can't evaporate given species from any section
      endif

      goto 100
c
c   step 4: when sectional compositions prevent achievement of equilibrium
c           gas phase concentrations we need to adjust the non limiting
c           species so that that species does achieve equilibrium.
c           For example:  if NH4Cl cannot evaporate because a section 
c           runs out of NH4 and this is the only particle with any Cl
c           then KNH4Cl will not equal [NH3]final [HCL]final becuase 
c           this equilibrium could not be achieved.  However the 
c           equilibrium transfer flux of NH4NO3 was calculated assuming 
c           that all of the NH4Cl would evaporate.  Since it didn't 
c           KNH4NO3 will not equal [NH3]final [HNO3]final.  To correct
c           in this case we will evaporate enough NH4NO3 to establish
c           equilibrium. 
c
c           only evaporating species can encounter this problem
c           only if NH4 becomes zero in a section can this problem occur
c
 180  if(ispsav.eq.kno3) then
         isp2=kcl                     ! correcting species
         g2=ghcl*1.0d6                ! original equilibrium for HCL (umol/m3)
      else
         isp2=kno3                    ! correcting species
         g2=ghno3*1.0d6               ! original equilibrium for HNO3 (umol/m3)
      endif
      bb=gnh3*1.0d6+dq(ispsav)+g2
      dq(isp2)=(-bb+sqrt(bb**2-4*dq(ispsav)*g2))/2.0d0

      dq(isp2) = dmax1(dq(isp2), dq(knh4)+dqsav(ispsav)-dq(ispsav) ! bkoo (10/07/03)
     &                          +dqsav(isp2)-q(ng+inh3)/prod)

      if(iter.gt.itmaxeq+1) dq(isp2)=0.0d0 ! bkoo (11/14/01)
c reset aerosol concentrations and dq's for NO3, NH4 and Cl
      do isec=1,nsecx2
         indx=(isec-1)*nsp
cbk         q(indx+kno3)=qsav(indx+kno3) ! initial aerosol no3
cbk         q(indx+knh4)=qsav(indx+knh4) ! initial aerosol nh4
cbk         q(indx+kcl)=qsav(indx+kcl)   ! initial aerosol Cl
         indx2=(isec-1)*nexti          ! bkoo (03/09/03)
         q(indx+kno3)=qsav(indx2+kno3) ! bkoo (03/09/03)
         q(indx+knh4)=qsav(indx2+knh4) ! bkoo (03/09/03)
         q(indx+kcl)=qsav(indx2+kcl)   ! bkoo (03/09/03)
      enddo
c restore frq - bkoo (10/07/03)
      do isp=1,nexti
        do isec=1,nsecx2
          frq(isec,isp) = frq0(isec,isp)
        enddo
      enddo
c adjust DQ
cbk      do i=1,nsp
      do i=1,nexti ! bkoo (03/09/03)
         if(i.eq.ispsav) then
            dq(i)=dqsav(i)-dq(i)
         elseif(i.eq.isp2) then
            dq(i)=dqsav(i)-dq(i)
         elseif(i.eq.knh4) then
            dq(i)=dqsav(i)
         else
            dq(i)=0.0d0
         endif
      enddo
      goto 80
c      endif
C      
C   STEP 4: ASSIGN EQUILIBRIUM VALUES TO GASES
c   Its possible to not be able to reach equilibrium for NH4Cl and NH4NO3
c   do to diferences between aerosol sectional vs. bulk composition 
C   So we treat NH3, HNO3, and HCl conservatively
C
 200  Q(NG+IH2SO4)=0.0D0              ! all sulfuric acid condenses
      Q(NG+IHNO3)= Q(NG+IHNO3)+q0(1)*PROD/GMW(IHNO3)
      Q(NG+INH3) = Q(NG+INH3) +q0(2)*PROD/GMW(INH3)
      Q(NG+IHCL) = Q(NG+IHCL) +q0(3)*PROD/GMW(IHCL)
      do isec=1,nsecx2
       Q(NG+IHNO3)= Q(NG+IHNO3)-Q((ISEC-1)*NSP+KNO3)*PROD/GMW(IHNO3)
       Q(NG+INH3) = Q(NG+INH3) -Q((ISEC-1)*NSP+KNH4)*PROD/GMW(INH3)
       Q(NG+IHCL) = Q(NG+IHCL) -Q((ISEC-1)*NSP+KCL) *PROD/GMW(IHCL)
      enddo
c     correct negative NH3 - bkoo (10/07/03)
      dq(KNH4) = q(ng+inh3) / prod ! umol/m3
      if(dq(KNH4).lt.-tinys) then
        iter = 0
 300    frt = 1.d0
        do isec = 1,nsecx2
          indx=(isec-1)*nsp       
          if(frq0(isec,KNH4).gt.0.d0) frt = dmax1( dmin1(q(indx+KNH4)
     &             /(-dq(KNH4)*frq0(isec,KNH4)*emw(KNH4)),frt), 0.d0)
        enddo
        frqtot = 0.d0
        do isec = 1,nsecx2
          indx=(isec-1)*nsp
          q(indx+KNH4) = dmax1(q(indx+KNH4)
     &                  +frt*dq(KNH4)*frq0(isec,KNH4)*emw(KNH4),0.d0)
          if(q(indx+KNH4).lt.tinys) frq0(isec,KNH4) = 0.d0
          frqtot = frqtot + frq0(isec,KNH4)
        enddo
        q(ng+inh3) = q(ng+inh3) - frt*dq(KNH4)*prod
        ! check if we should evaporate more
        dq(KNH4) = (1.d0 - frt) * dq(KNH4)
        if(dq(KNH4).lt.-tinys) then
          if(frqtot.gt.tinys) then ! we have sections which are not empty    
            if(iter.le.itmaxeq) then ! check infinite loop
              iter = iter + 1
              do isec = 1,nsecx2
                frq0(isec,KNH4) = frq0(isec,KNH4) / frqtot
              enddo
              goto 300
            endif
          endif
          ! we need to evaporate more to achieve equilibrium
          ! but we completely evaporate the species in all sections
          ! or exceed itermax
c          write(*,*)'EQUI-F2: dq(KNH4)=',dq(KNH4),' iter=',iter ! bkoo_dbg
        endif
      endif
C   DO ORGANICS  **ASSUME PS HAS NO COMPOSITION DEPENDANCE **
cbk removed - bkoo (03/09/03)
cbk      DO IOG = 1,NORG-1
cbk       Q(NG+IHCL+IOG) =PS(1,IHCL+NORG)/(1.01325d-1*pres) ! pres (bkoo, 06/09/00)
cbk      ENDDO
  
      call negchk(t,q,nsecx2)

      call eqneut(0, q) ! tmg (12/05/02)

      if(aerm.eq.'EQUI') then ! (tmg,01/31/02)
C
C   STEP 5: CALL EQUILIBRIUM CODE FOR EACH SECTION TO DETERMINE WATER
C
cbk        CALL STEP(nsec,Q)  ! tmg (10/22/02) bkoo (03/07/03)
C
C   STEP 6: Determine new diameter
C
        call ddiameter(q)
      endif
C
      RETURN
      END     

c ----------------------------------------------------------------------------
c neutralize the acidic or basicness of the aerosol in each section
c
c is   : 0 (false) for equilibrium section
c        1 (true)  for dynamic     section
c q    : aerosol [ug/m3] and gas [ppm] species array
c        dimension - ntotal  for equilibrium section of EQUI &
c                                dynamic     section of MADM
c                    ntotale for equilibrium section of HYBR
c                    ntotald for dynamic     section of HYBR
c
c CHTOT: Total charge balance [umole/m3]
c
      subroutine eqneut(is, q)

      include 'dynamic.inc'

      real*8      q(ntotal),chtot,chmove,chno3,chcl,chnh4,prod
      integer     nsx,ng,na,chspc(nexti)
      integer     jani,jcat ! bkoo (10/07/03)
      integer     kerr ! bkoo (01/06/04)

c charges of H2O, Na, SO4, NO3, NH4, Cl
      data chspc / 0, 1, -2, -1, 1, -1 /

      nsx  = nsecx * is + nsecx2 * (1 - is)
      ng   = nsp * nsx
      prod = rgas*temp/(1.01325d5*pres) ! PPM / prod = UMOLE/M3

c      if(is) call negchk(0.d0,q,nsecx)

      do i = 1, nsx
        na = (i-1)*nsp
        chtot = 0.0d0
        jani = 0 ! bkoo (10/07/03)
        jcat = 0 ! bkoo (10/07/03)
        do j = 2, nexti ! the first element is H2O
          chtot = chtot + DBLE(chspc(j)) * q(na+j) / emw(j)
          if(chspc(j).lt.0 .and. q(na+j).gt.tinys) jani = 1 ! bkoo (10/07/03)
          if(chspc(j).gt.0 .and. q(na+j).gt.tinys) jcat = 1 ! bkoo (10/07/03)
        enddo
c BASIC:
        if(chtot .gt. 0.0d0 .and. jani .eq. 0) then ! bkoo (10/07/03)
          chno3 = (q(ng+IHNO3)/prod)
          chcl  = (q(ng+IHCL) /prod)
          chnh4 = (q(na+KNH4) /emw(KNH4))
c first move HCl in
          chmove     = max(min(chtot, chcl),0.0d0)
          q(na+KCL)  = q(na+KCL)  + chmove * emw(KCL)
          q(ng+IHCL) = max(q(ng+IHCL) - chmove * prod,0.0d0)
          chtot      = chtot - chmove
c then move HNO3 in
          chmove     = max(min(chtot, chno3),0.0d0)
          q(na+KNO3) = q(na+KNO3) + chmove * emw(KNO3)
          q(ng+IHNO3)= max(q(ng+IHNO3)- chmove * prod,0.0d0)
          chtot      = chtot - chmove
c then move NH4 out
          chmove     = max(min(chtot, chnh4),0.0d0)
          q(na+KNH4) = max(q(na+KNH4) - chmove * emw(KNH4),0.0d0)
          q(ng+INH3) = q(ng+INH3) + chmove * prod
          chtot      = chtot - chmove
        endif                                       ! bkoo (10/07/03)
c ACIDIC:
        if(chtot .lt. 0.0d0 .and. jcat .eq. 0) then ! bkoo (10/07/03)
          chtot = -1.0d0 * chtot
          chno3 = (q(na+KNO3)/emw(KNO3))
          chcl  = (q(na+KCL) /emw(KCL))
          chnh4 = (q(ng+INH3)/prod)
c first move NH3 in
          chmove     = max(min(chtot, chnh4),0.0d0)
          q(na+KNH4) = q(na+KNH4) + chmove * emw(KNH4)
          q(ng+INH3) = max(q(ng+INH3) - chmove * prod,0.0d0)
          chtot      = chtot - chmove
c then move HCl out
          chmove     = max(min(chtot, chcl),0.0d0)
          q(na+KCL)  = max(q(na+KCL)  - chmove * emw(KCL),0.0d0)
          q(ng+IHCL) = q(ng+IHCL) + chmove * prod
          chtot      = chtot - chmove
c then move HNO3 out
          chmove     = max(min(chtot, chno3),0.0d0)
          q(na+KNO3) = max(q(na+KNO3) - chmove * emw(KNO3),0.0d0)
          q(ng+IHNO3)= q(ng+IHNO3)+ chmove * prod
          chtot      = chtot - chmove
        endif
      enddo
c
c     check the charge balance in each section - bkoo (01/06/04)
c
c bkoo_dbg =>
      do i = 1, nsx
        na = (i-1)*nsp
        kerr = 0
        if (q(na+KNA).gt.tinys .or. q(na+KNH4).gt.tinys) then
          if (q(na+KSO4).le.tinys .and. q(na+KNO3).le.tinys .and.
     &        q(na+KCL).le.tinys) kerr = 1
cbk          if ( (2.*q(na+KSO4)/emw(KSO4)+q(na+KNO3)/emw(KNO3)+q(na+KCL)
cbk     &            /emw(KCL))/(q(na+KNA)/emw(KNA)+q(na+KNH4)/emw(KNH4))
cbk     &            .gt.0.2e10) kerr = 1
        else
          if (q(na+KNO3).gt.tinys .or. q(na+KCL).gt.tinys) kerr = 1
        endif
c        if (kerr.ne.0) then
c          call get_param(igrdchm,ichm,jchm,kchm,iout,idiag)
c          write(iout,'(//,A)')'CHECK in EQNEUT'
c          write(iout,*)' q(',i,'): ',(q(na+j),j=2,nexti)
c          write(iout,*)' igrd,i,j,k: ',igrdchm,ichm,jchm,kchm
c        endif
      enddo
c bkoo_dbg <=

      return
      end

