      subroutine aqoperator1(tbeg,deltat,gas,aerosol,lwc,
     &  rh,temp,p,iaq)
c
c
c
c.....INPUTS
c     TBEG : BEGINNING OF THE INTEGRATION (in min)
c     DELTAT : INTEGRATION TIME (in min)
c     GAS(NGAS) : GAS-PHASE CONCENTRATION VECTOR (in ppm)
c     AEROSOL(NSECT,NAERS) : AEROSOL SPECIES CONCENTRATION (in ug/m3) 
c     LWC : LIQUID WATER CONTENT (in g/m3)
c     RH : RELATIVE HUMIDITY  (in 0-1 scale)
c     TEMP : TEMPERATURE (in K)
c     P : PRESSURE (in atm)
c     IAQ : FLAG: 1 AT FIRST CALL, 0 THERE AFTER
c
c.....OUTPUTS
c     GAS(NGAS) : UPDATED GAS-PHASE CONCENTRATIONS (in ppm)
c     AEROSOL(NSECT,NAERS) : UPDATED AEROSOL SPECIES CONCENTRATION (in ug/m3) 
c     IAQ : FLAG SET TO ZERO
c
c.....VARIABLES
c     GCON(NGAS) : GAS-PHASE CONCENTRATIONS (in ppm) (THROUGH COMMON TO RATES)
c
c
c.....DIFFERENTIAL EQUATIONS ARE SOLVED FOR THE FOLLOWING SPECIES      
c
c       TOTAL                        GAS             AQUEOUS   
c-----------------------------------------------------------------------
c     1. FORMALDEHYDE TOTAL        1. SO2       1.  S(IV)
c     2. FORMIC ACID TOTAL         2. H2O2      2.  H2O2(aq)
c                                  3. NH3       3.  NITRATE
c                                               4.  CHLORIDE
c                                               5.  AMMONIUM
c                                               6.  SULFATE
c                                               7.  HSO5-
c                                               8.  HMSA
c
c
c.....Y MATRIX STRUCTURE  EVERYTHING IS (ug/m3)
c     ONLY 11 ODEs ARE SOLVED FOR THE WHOLE DISTRIBUTION
c
c  YAQ(1) = TOTAL FORMALDEHYDE
c  YAQ(2) = TOTAL FORMIC ACID
c  YAQ(3) = SO2 (g)
c  YAQ(4) = H2O2 (g)
c  YAQ(5) = NH3 (g)
c  YAQ(6) = SO2 (aq)
c  YAQ(7) = H2O2 (aq)
c  YAQ(8) = NITRATE (aq)
c  YAQ(9) = CHLORIDE (aq)          
c  YAQ(10) = AMMONIUM (aq)
c  YAQ(11) = SULFATE (aq)
c  YAQ(12) = HSO5- (aq)
c  YAQ(13) = HMSA (aq)
c     
      include 'aerpar.inc'
      include 'droppar.inc'
      include 'dropcom.inc'
      real*4 gas(ngas_aq), aerosol(nsect,naers)
      real*4 sodium(nsect)
      real*4 deltat, lwc, rh, temp
      real*4 gascon(ngas_aq)
      real*4 yaq(meqn1)
      real*4 watcont, watvap
      real*4 check(3,naers)
c
      common / gascon / gascon
      common / epcomy / ymin, hmax
      common / integration / wat1,relh1,therm,dt,dtemp
      common / aqrate / sodium, pres
      common / count/ icount, iprint
      common/aqchem1/watcont,watvap,crustal,salt
c
c     INITIALIZATION (ONLY ON FIRST CALL)
c     (CALCULATION OF SECTIONAL DIAMETERS, LOADING OF MOLEC. WEIGHTS)

      if (iaq .eq. 1) then
      call dropinit
      iaq = 0
      endif
      pres = p                                   ! pressure in atm
      call constants(temp)
c
c     ZERO THE YAQ MATRIX
c
      do i=1, meqn1
      yaq(i)=0.0
      enddo
c      
c     TRANSFER VIA COMMON VARIABLES TO AQRATE
c
      wat1 = lwc
      watcont = lwc
      relh1 = rh
      therm = temp
      dt = deltat
      icount=0
      iprint=20
c
c     TRANSFER OF ALL GAS-PHASE CONCENTRATIONS TO GASCON AND THEN
c     THROUGH COMMON BLOCK TO RATES
c
      do i=1, ngas_aq
      gascon(i) = gas(i)
      enddo
c
c     UNIT CONVERSION FACTORS     
c      
      fso2=1000.*pres*wso2/(8.314e-2*temp)     !SO2   conver. factor from ppm to ug/m3
      fh2o2=1000.*pres*wh2o2/(8.314e-2*temp)   !H2O2  conver. factor from ppm to ug/m3
      fhcho=1000.*pres*whcho/(8.314e-2*temp)   !HCHO  conver. factor from ppm to ug/m3
      fhcooh=1000.*pres*whcooh/(8.314e-2*temp) !HCOOH conver. factor from ppm to ug/m3
      fammon=1000.*pres*wnh3/(8.314e-2*temp)   !NH3   conver. factor from ppm to ug/m3
c    
c      
c   *** BEGINNING OF AQUEOUS-PHASE CALCULATION
c   LOADING OF THE CONCENTRATIONS IN THE YAQ VECTOR
c
c   THE TOTAL FORMALDEHYDE/FORMIC ACID ARE TRANSFERRED AS GAS-PHASE SPECIES
c    (THEY ARE NOT INCLUDED IN THE AQUEOUS OR AEROSOL MATRICES)
      yaq(1) = gas(nghcho)*fhcho           ! total HCHO (ug/m3)
      yaq(2) = gas(nghcooh)*fhcooh         ! total HCOOH (ug/m3)
c
c   GAS-PHASE SPECIES
c
      yaq(3) = gas(ngso2)*fso2             ! total SO2 (ug/m3)
      yaq(4) = gas(ngh2o2)*fh2o2           ! total H2O2 (ug/m3)                    
      yaq(5) = gas(nga) *fammon            ! NH3(g) in ug/m3
c
c   AQUEOUS-PHASE SPECIES
c
      do isect=1, nsect
c
c A FIXED SECTIONAL DIAMETER IS USED FOR THE ACTIVATION
c CALCULATION. THIS IS COMPATIBLE WITH THE ASSUMPTION THAT
c A CLOUD OR FOG EXISTS IF THE RH EXCEEDS A THRESHOLD VALUE
c (FOR EXAMPLE RH>85%).
c
c
c     ** THE DROPLETS ARE ASSUMED TO BY WITHOUT SO2 OR
c        H2O2 IN THE BEGINNING OF A TIMESTEP **
c      (THIS IS DONE TO AVOID CARRYING THE S(IV)/H2O2 CONCENTRATIONS
c      IN THE 3D SIMULATION **
      yaq(6) = 1.e-10     ! SO2 (aq) in ug/m3
      yaq(7) = 1.e-10     ! H2O2 (aq) in ug/m3
c
      if (daer(isect) .ge. dactiv) then
      yaq(8) = yaq(8) + aerosol(isect,nan)       ! NITRATE (aq) in ug/m3 
      yaq(9) = yaq(9) + aerosol(isect,nac)       ! CHLORIDE (aq) in ug/m3
      yaq(10) = yaq(10) + aerosol(isect,naa)     ! AMMONIUM (aq) in ug/m3   
      yaq(11) = yaq(11) + aerosol(isect,na4)     ! SULFATE (aq) in ug/m3
      yaq(12) = yaq(12) + aerosol(isect,nahso5)  ! HSO5- in ug/m3
      yaq(13) = yaq(13) + aerosol(isect,nahmsa)  ! HMSA in ug/m3
      endif
      enddo
      
c
c     CALCULATION OF THE DISSOLUTION OF THE AVAILABLE HNO3 AND HCl
c     TO THE DROPLETS/PARTICLES. WE ASSUME THAT THE TIMESCALE OF DISSOLUTION
c     IS SMALL AND THAT ALL HNO3 AND HCl ARE DISSOLVED (ZERO VAPOR PRESSURE)
c
      fhno3=1000.*whno3*pres/(8.314e-2*temp) !HNO3 conver. factor from ppm to ug/m3
      fhcl=1000.*whcl*pres/(8.314e-2*temp)   !HCl  conver. factor from ppm to ug/m3
c
      hno3 = gas(ngcn)*fhno3   ! HNO3(g) (in ug/m3)
      hcl =  gas(ngcc)*fhcl    ! HCl(g) (in ug/m3)
c      
c     SAVE THE OLD AQUEOUS-PHASE CONCENTRATIONS BEFORE THE INTEGRATION
c
      tnitold=yaq(8)
      chlorold=yaq(9)
      ammonold=yaq(10)
      sulfateold=yaq(11)
      hso5old=yaq(12)
      hmsaold=yaq(13)
c
c      ** ADD THE GAS-PHASE CONCENTRATIONS TO THE TOTAL FOR HCl AND HNO3 **
      gas(ngcn) = 0.0          ! all HNO3 is dissolved
      gas(ngcc) = 0.0          ! all HCl is dissolved
      yaq(8) = yaq(8)+hno3     ! HNO3 increase
      yaq(9) = yaq(9)+hcl      ! HCl increase
c
c
c.....THE YAQ MATRIX HAS BEEN INITIALIZED
c
c     CALCULATION OF SODIUM VECTOR (NEEDED FOR THE INTEGRATION ROUTINE)
c       (TRANSFERRED VIA COMMON STATEMENT AQRATESA)
c
      do isect=1, nsect
      sodium(isect) = aerosol(isect,nas)
      enddo
c
c     VARIABLES FOR INTEGRATION
c
      stime = tbeg*60.                        ! in seconds
      stout = (tbeg+deltat) * 60.              ! in seconds
c
c     CALCULATE CRUSTAL SPECIES CONCENTRATION INSIDE DROPLETS
c     (IT IS USED IN THE AQUEOUS-PHASE CHEMISTRY CALCULATION TO
c      ESTIMATE Fe and Mn CONCENTRATIONS
c
      crustal = 0.
      salt =0.
      do isect=1,nsect
c
      if (dd(isect) .ge. dactiv) then
      crustal=crustal+aerosol(isect,nar)
      salt=salt+aerosol(isect,nas)
      endif
c
      enddo   
c
c     INTEGRATE
c 
      sbef=0.
      do i=1,nsect
       sbef=sbef+aerosol(i,na4)*32./96.+aerosol(i,nahso5)*32./113.+
     &           aerosol(i,nahmsa)*32./111.
      enddo
      sbef=sbef+yaq(3)*32./64.

      call aqintegr1(meqn1, yaq, stime, stout, temp, p)

c
c     CALCULATE THE DIFFERENCES
c
      dnit = yaq(8) - tnitold
      dchlor = yaq(9) - chlorold
      dammon = yaq(10) - ammonold
      dsulf = yaq(11) - sulfateold
      dhso5 = yaq(12) - hso5old
      dhmsa = yaq(13) - hmsaold
c
c     ADD THE MASS CHANGE TO THE CORRESPONDING SECTION 
c
      do isect=1, nsect
      aerosol(isect,nan)=aerosol(isect,nan)+dnit*fdist(isect)    ! NITRATE (aq) in ug/m3
      aerosol(isect,nac)=aerosol(isect,nac)+dchlor*fdist(isect)  ! CHLORIDE (aq) in ug/m3
      aerosol(isect,naa)=aerosol(isect,naa)+dammon*fdist(isect)  ! AMMONIUM (aq) in ug/m3   
      aerosol(isect,na4)=aerosol(isect,na4)+dsulf*fdist(isect)   ! SULFATE (aq) in ug/m3
      aerosol(isect,nahso5)=aerosol(isect,nahso5)
     & +dhso5*fdist(isect)                                       ! hso5 (aq) in ug/m3
      aerosol(isect,nahmsa)=aerosol(isect,nahmsa)
     & +dhmsa*fdist(isect)                                       ! hmsa (aq) in ug/m3
      enddo
      
c
      saft=0.
      do i=1,nsect
       saft=saft+aerosol(i,na4)*32./96.+aerosol(i,nahso5)*32./113.+
     &           aerosol(i,nahmsa)*32./111.
      enddo
      saft=saft+(yaq(3)+0.7901*yaq(6))*32./64.
c     CHECK IF THE ALGORITHM RESULTED IN NEGATIVE CONCENTRATIONS
c     REPORT THE CORRECTIONS AT THE WARNING FILE  
c
c      do isect=1, nsect
c      do isp=1, naers     
c      if (aerosol(isect,isp) .lt. 0.0) then
c      write(2,*)'NEGATIVE CONCENTRATION AT', stime
c      write(2,*)'SECTION=',isect, 'SPECIES=',isp, aerosol(isect,isp)
c     
c      write(6,*)'NEGATIVE CONCENTRATION AT', stime      !added by TMR
c      write(6,*)'SECTION=',isect, ' SPECIES=',isp, aerosol(isect,isp)
c
cckf      aerosol(isect,isp)=1.e-12 
c      endif 
c      enddo 
c      enddo 

c  check(1,isp) = total aerosol concentration for species isp     
c  check(2,isp) = concentration below zero for species isp
c  check(3,isp) = (-) total mass added to the system to get a positive conc.

      do j = 1,3
      do isp = 1,naers
      check(j,isp) = 0.0
      enddo
      enddo
      
      do isp = 1,naers
      do isect=1,nsect
      check(1,isp) = check(1,isp) + aerosol(isect,isp)
      enddo
      enddo
      
c If negative in a section, change to a small positive value
c Account for this mass change by subtracting from other 
c non-negative sections.      
      
      do isp = 1,naers
      do isect = 1,nsect
      if (aerosol(isect,isp) .lt. 0.0) then
      check(2,isp) = check(2,isp)+aerosol(isect,isp)
      check(3,isp) = check(3,isp)+aerosol(isect,isp)-1.e-12 ! bkoo (10/07/03)
      aerosol(isect,isp) = 1.e-12
      endif
      enddo
      enddo
      
      do isect = 1,nsect
      
      if (aerosol(isect,nan) .gt. 1.e-12) then
      aerosol(isect,nan)=aerosol(isect,nan)+(aerosol(isect,nan))*
     & check(3,nan)/(check(1,nan)-check(2,nan))
      endif
      if (aerosol(isect,nac) .gt. 1.e-12) then
      aerosol(isect,nac)=aerosol(isect,nac)+(aerosol(isect,nac))*
     & check(3,nac)/(check(1,nac)-check(2,nac))
      endif
      if (aerosol(isect,naa) .gt. 1.e-12) then
      aerosol(isect,naa)=aerosol(isect,naa)+(aerosol(isect,naa))*
     & check(3,naa)/(check(1,naa)-check(2,naa))
      endif                  

c For HSO5 and HMSA, if negative and there is not enough mass in the
c remaining aerosol sections to account for the added mass, subtract
c from sulfate.
     
      if (check(1,nahso5)-check(2,nahso5) .lt. 0.-check(3,nahso5)) then
      check(3,na4) = check(3,na4)+check(3,nahso5)*96./113.
      else if (aerosol(isect,nahso5) .gt. 1.e-12) then
      aerosol(isect,nahso5)=aerosol(isect,nahso5)+
     & (aerosol(isect,nahso5))*check(3,nahso5)/
     & (check(1,nahso5)-check(2,nahso5))
      endif
      
      if (check(1,nahmsa)-check(2,nahmsa) .lt. 0.-check(3,nahmsa)) then
      check(3,na4) = check(3,na4)+check(3,nahmsa)*96./111.
      else if (aerosol(isect,nahmsa) .gt. 1.e-12) then
      aerosol(isect,nahmsa)=aerosol(isect,nahmsa)+
     & (aerosol(isect,nahmsa))*check(3,nahmsa)/
     & (check(1,nahmsa)-check(2,nahmsa))
      endif  
                  
      if (aerosol(isect,na4) .gt. 1.e-12) then
      aerosol(isect,na4)=aerosol(isect,na4)+(aerosol(isect,na4))*
     & check(3,na4)/(check(1,na4)-check(2,na4))
      endif
     
      enddo 

      saft1=0.
      do i=1,nsect
       saft1=saft1+aerosol(i,na4)*32./96.+aerosol(i,nahso5)*32./113.+
     &           aerosol(i,nahmsa)*32./111.
      enddo
      saft1=saft1+(yaq(3)+0.7901*yaq(6))*32./64.

c
c     RETURN YAQ RESULTS TO THEIR ORIGINAL MATRICES (GAS,AEROSOL)
c
c     ** GAS-PHASE SPECIES **     
c     IMPORTANT : THE DISSOLVED S(IV) and H2O2 ARE TRANSFERRED
c                 BACK TO THE GAS PHASE
c     
      gas(nghcho) = yaq(1)/fhcho              ! total HCHO (ppm)
      gas(nghcooh)= yaq(2)/fhcooh             ! total HCOOH (ppm)
      gas(ngso2)= (yaq(3)+yaq(6)*0.7901)/fso2       ! SO2(g) (ppm)
      gas(ngh2o2)=(yaq(4)+yaq(7))/fh2o2             ! H2O2(g) (ppm)                    
      gas(nga)=yaq(5)/fammon               ! NH3(g) (ppm)
c
c THE CODE NEGLECTS THE CHANGE IN THE SIZE DISTRIBUTION SHAPE
c OF THE NONVOLATILE AEROSOL COMPONENTS BECAUSE OF THE SULFATE
c PRODUCTION. THE CHANGE IN THE SULFATE DISTRIBUTION IS CALCULATED
c USING THE DISTRIBUTION FUNCTIONS WHILE THE CHANGE IN THE
c DISTRIBUTION OF THE VOLATILE INORGANIC AEROSOL COMPONENTS WILL
c CALCULATED BY THE AEROSOL MODULE AFTER CLOUD EVAPORATION
c
c
c.....END OF CODE
      return
      end
