module module_mosaic_ext
contains

      ! scm_temp - subr dumpxx is temporary for scm testing
      subroutine dumpxx( dumpch, dtchem, t_k, p_atm, ah2o, &
         jaerosolstate, jphase, jhyst_leg, &
         aer, gas, num_a, water_a, water_a_hyst, dp_dry_a, &
         mosaic_vars_aa )

      use module_data_mosaic_kind,  only: r8

      use module_data_mosaic_aero,  only: &
         iso4_a,       ino3_a,       icl_a,        inh4_a,       &
         ina_a,        ioin_a,       ioc_a,        ibc_a,        &
         ipcg2_b_c_a,  ipcg2_b_o_a,  ipcg1_b_c_a,  ipcg1_b_o_a,  &
         iopcg1_b_c_a, iopcg1_b_o_a, ipcg2_f_c_a,  ipcg2_f_o_a,  &
         ipcg1_f_c_a,  ipcg1_f_o_a,  iopcg1_f_c_a, iopcg1_f_o_a, &
         iant1_c_a,    iant1_o_a,    ibiog1_c_a,   ibiog1_o_a,   &
         aer_name, gas_name, &
         jtotal, naer, nbin_a, nbin_a_max, ngas_aerchtot, ngas_volatile, &
         mosaic_vars_aa_type

      implicit none

      integer, intent(in), dimension(nbin_a_max) :: jaerosolstate, jphase, jhyst_leg

      real(r8), intent(in) :: dtchem, t_k, p_atm, ah2o
      real(r8), intent(in), dimension(naer,3,nbin_a_max) :: aer
      real(r8), intent(in), dimension(ngas_volatile) :: gas
      real(r8), intent(in), dimension(nbin_a_max) :: num_a, water_a, water_a_hyst, dp_dry_a

      character(len=*), intent(in) :: dumpch

      type (mosaic_vars_aa_type), intent(in) :: mosaic_vars_aa

      integer, parameter :: lunaa = 131
      integer, parameter :: naer_dump = 24
      integer :: ii, iv, ix, jj, kk, ktau
      integer :: iaer_dump(naer_dump) = 24
      character(len=80)  :: fmtaa
      character(len=2)   :: tmpch2
      character(len=256) :: tmpchaa
      character(len= 32) :: tmpchbb

!     return  ! scm_temp - for non-scm might want to activate this line

      if ( mosaic_vars_aa%idiagbb_host < 100 ) return

      ii = mosaic_vars_aa%hostgridinfo(2)
      jj = mosaic_vars_aa%hostgridinfo(3)
      kk = mosaic_vars_aa%hostgridinfo(4)
      if ( ii*jj*kk /= 1 ) return

      ktau = mosaic_vars_aa%it_host
      iaer_dump = &
         (/ iso4_a,       ino3_a,       icl_a,        inh4_a,      &
            ina_a,        ioin_a,       ioc_a,        ibc_a,       &
            ipcg2_b_c_a,  ipcg2_b_o_a,  ipcg1_b_c_a,  ipcg1_b_o_a, &
            iopcg1_b_c_a, iopcg1_b_o_a, ipcg2_f_c_a,  ipcg2_f_o_a, &
            ipcg1_f_c_a,  ipcg1_f_o_a,  iopcg1_f_c_a, iopcg1_f_o_a, &
            iant1_c_a,    iant1_o_a,    ibiog1_c_a,   ibiog1_o_a /)

! qakola
      tmpch2 = dumpch
      write(lunaa,'(/2a,i5,3i3)') tmpch2, 'dump', &
         ktau, ii, jj, kk
      write(lunaa,'(a,f11.2,f11.4,f11.4      )') 't p a ', t_k, p_atm, ah2o
!     if (tmpch2 == 'aa') &
!     write(lunaa,'(a,f11.2                  )') 'cairm3', cair_mol_m3

!     write(lunaa,'(a,1p,8e11.3)') 'gas1:4    ', gas(1:4)

      fmtaa = '(a,8i11)'
      if (nbin_a > 8) fmtaa = '(a,8i11/(10x,8i11))'
      write(lunaa,fmtaa) 'jstate    ', jaerosolstate(1:nbin_a)
      write(lunaa,fmtaa) 'jphase    ', jphase(1:nbin_a)
      write(lunaa,fmtaa) 'jhyst     ', jhyst_leg(1:nbin_a)

      fmtaa = '(a,1p,8e11.3)'
      if (nbin_a > 8) fmtaa = '(a,1p,8e11.3/(10x,1p,8e11.3))'
      write(lunaa,fmtaa) 'num       ', num_a(1:nbin_a)
      write(lunaa,fmtaa) 'dpdry     ', dp_dry_a(1:nbin_a)
      write(lunaa,fmtaa) 'water     ', water_a(1:nbin_a)
      write(lunaa,fmtaa) 'hyswtr    ', water_a_hyst(1:nbin_a)

!     write(lunaa,'(a,1p,8e11.3)') 'so4       ', aer(iso4_a,jtotal,1:nbin_a)
!     write(lunaa,'(a,1p,8e11.3)') 'no3       ', aer(ino3_a,jtotal,1:nbin_a)
!     write(lunaa,'(a,1p,8e11.3)') 'cl        ', aer(icl_a, jtotal,1:nbin_a)
!     write(lunaa,'(a,1p,8e11.3)') 'nh4       ', aer(inh4_a,jtotal,1:nbin_a)
!!!   write(lunaa,'(a,1p,8e11.3)') 'msa       ', aer(imsa_a,jtotal,1:nbin_a)
!!!   write(lunaa,'(a,1p,8e11.3)') 'co3       ', aer(ico3_a,jtotal,1:nbin_a)
!     write(lunaa,'(a,1p,8e11.3)') 'na        ', aer(ina_a, jtotal,1:nbin_a)
!!!   write(lunaa,'(a,1p,8e11.3)') 'ca        ', aer(ica_a, jtotal,1:nbin_a)
!     write(lunaa,'(a,1p,8e11.3)') 'oin       ', aer(ioin_a,jtotal,1:nbin_a)
!     write(lunaa,'(a,1p,8e11.3)') 'oc        ', aer(ioc_a, jtotal,1:nbin_a)
!     write(lunaa,'(a,1p,8e11.3)') 'bc        ', aer(ibc_a, jtotal,1:nbin_a)

!need to dump soa species
      do ix = 1, naer_dump
         iv = iaer_dump(ix)
         if (iv < 1 .or. iv > naer) cycle
         fmtaa = '(a,1p,8e11.3)'
         write(tmpchaa,fmtaa) aer_name(iv)(1:10), aer(iv,jtotal,1:min(nbin_a,8))
         if (iv <= ngas_aerchtot) then
            write(tmpchbb,fmtaa) gas_name(iv)(1:10), gas(iv)
            tmpchaa = trim(tmpchaa) // '  ' // trim(tmpchbb)
         end if
         write(lunaa,'(a)') trim(tmpchaa)
         fmtaa = '(10x,1p,8e11.3)'
         if (nbin_a > 8) write(lunaa,fmtaa) aer(iv,jtotal,9:nbin_a)
      end do ! ix

      return
      end subroutine dumpxx

  !***********************************************************************
  ! determines phase state of an aerosol bin. includes kelvin effect.
  !
  ! author: Rahul A. Zaveri
  ! update: Sep 2015
  !-----------------------------------------------------------------------
  subroutine aerosol_phase_state( ibin, jaerosolstate, jphase, aer,  &
       jhyst_leg, electrolyte, epercent, kel, activity, mc, num_a, mass_wet_a, mass_dry_a,    &
       mass_soluble_a, vol_dry_a, vol_wet_a, water_a, water_a_hyst, water_a_up, aH2O_a,     &
       aH2O, ma, gam, log_gamZ, zc, za, gam_ratio,     &
       xeq_a, na_Ma, nc_Mc, xeq_c,       mw_electrolyte, partial_molar_vol, sigma_soln, T_K, & ! RAZ deleted a_zsr
       RH_pc, mw_aer_mac, dens_aer_mac, sigma_water, Keq_ll, Keq_sl, MW_a, MW_c,             &
       growth_factor, MDRH, MDRH_T, molality0, rtol_mesa, jsalt_present, jsalt_index,       &
       jsulf_poor, jsulf_rich, phi_salt_old,                                     &
       kappa_nonelectro, mosaic_vars_aa )

    use module_data_mosaic_aero,  only: r8, nbin_a_max,                                   &
         ngas_aerchtot, ngas_volatile, nelectrolyte,                                      &!Parameters
         Ncation, naer, jtotal, all_solid, jhyst_up, all_liquid, Nanion, nrxn_aer_ll,     &
         nrxn_aer_sl, nsalt, MDRH_T_NUM, jsulf_poor_NUM, jsulf_rich_NUM,                  &!Parameters
         inh4_a, ina_a, ica_a, ico3_a, imsa_a, icl_a, ino3_a, iso4_a,                     & ! TBD
         a_zsr,  b_zsr,  aw_min,                                                          &! RAZ added a_zsr,  b_zsr,  aw_min
         mosaic_vars_aa_type

    implicit none
    !Intent -ins

    integer, intent(in):: ibin
    integer, intent(in), dimension(nsalt) :: jsalt_index
    integer, intent(in), dimension(jsulf_poor_NUM) :: jsulf_poor
    integer, intent(in), dimension(jsulf_rich_NUM) :: jsulf_rich

    real(r8), intent(in) :: aH2O,T_K,RH_pc,rtol_mesa
    real(r8), intent(in), dimension(naer) :: mw_aer_mac,dens_aer_mac
    real(r8), intent(in), dimension(Ncation) :: zc,MW_c
    real(r8), intent(in), dimension(Nanion)  :: za,MW_a
    real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
    real(r8), intent(in), dimension(ngas_aerchtot) :: partial_molar_vol
    real(r8), intent(in), dimension(naer) :: kappa_nonelectro

    !Intent - inout
    integer, intent(inout), dimension(nsalt) :: jsalt_present
    integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate,jphase,jhyst_leg

    real(r8), intent(inout) :: sigma_water
    real(r8), intent(inout), dimension(Ncation) :: nc_Mc,xeq_c
    real(r8), intent(inout), dimension(Nanion)  :: xeq_a,na_Ma
    real(r8), intent(inout), dimension(nbin_a_max) :: num_a,mass_wet_a,mass_dry_a
    real(r8), intent(inout), dimension(nbin_a_max) :: mass_soluble_a,gam_ratio
    real(r8), intent(inout), dimension(nbin_a_max) :: vol_dry_a,vol_wet_a,water_a
    real(r8), intent(inout), dimension(nbin_a_max) :: water_a_hyst,water_a_up,aH2O_a
    real(r8), intent(inout), dimension(nbin_a_max) :: sigma_soln,growth_factor,MDRH
    real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: molality0            ! RAZ 5/20/2014
    real(r8), intent(inout), dimension (nrxn_aer_ll) :: Keq_ll
    real(r8), intent(inout), dimension (nrxn_aer_sl) :: Keq_sl
    real(r8), intent(inout), dimension(MDRH_T_NUM) :: MDRH_T
    real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kel
    real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
    real(r8), intent(inout), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
    real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
    real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: epercent
    real(r8), intent(inout), dimension(nsalt) :: phi_salt_old

    type (mosaic_vars_aa_type), intent(inout) :: mosaic_vars_aa

    ! local variables
    character(len=256) :: errmsg
    integer, parameter :: aer_pha_sta_diagaa = -1 !BALLI- changed from 100 to -1
    integer, parameter :: iter_kelvin_method =  3
    ! iter_kelvin_method = 1 - use rahuls original iteration method
    ! iter_kelvin_method = 2 - use bisection
    ! iter_kelvin_method = 3 - start with rahuls original iteration method, but if it fails, switch to bisection
    integer, parameter :: iter_kelvin_meth1_max = 10
    integer, parameter :: iter_kelvin_meth2_max = 100
    integer :: iaer, iv, itmpa
    integer :: iter_kelvin, iter_kelvin_meth1, iter_kelvin_state
    integer :: js, je

    real(r8) :: aer_H
    real(r8):: aH2O_range_bisect_toler
    real(r8) :: aH2O_a_new, aH2O_a_old, aH2O_a_oldn, aH2O_a_oldp, aH2O_a_del_state3
    real(r8), dimension(nbin_a_max) :: DpmV
    real(r8), dimension(nbin_a_max) :: kelvin
    real(r8) :: kelvin_old, kelvin_oldn, kelvin_oldp
    real(r8) :: kelvin_toler
    real(r8) :: rel_err, rel_err_old, rel_err_old2, rel_err_oldn, rel_err_oldp
    real(r8) :: term, tmpa
    real(r8) :: water_a_old, water_a_oldn, water_a_oldp

    if (aer_pha_sta_diagaa >= 3) &
    write(*,'(/a,5i5,2f12.8,1p,2e11.3)') 'aer_pha_sta_a', ibin, jhyst_leg(ibin), jaerosolstate(ibin), -1, 0, aH2O, aH2O_a(ibin)
    !aH2O = RH_pc*0.01 !**BALLI, this is already done in init subr
    aH2O_a(ibin) = aH2O
    kelvin(ibin) = 1.0
    do iv = 1, ngas_aerchtot
       kel(iv,ibin) = 1.0
    enddo

!   if(RH_pc .le. 97.0)then     ! RAZ
!      kelvin_toler = 1.e-4
!   else
!      kelvin_toler = 1.e-10    ! RAZ
!   endif
! define error tolerances become stricter as aH2O approaches 1.0
    kelvin_toler = 1.e-6_r8 * max( 1.0_r8-aH2O, 1.0e-4_r8 )
    aH2O_range_bisect_toler = 1.e-6_r8 * max( 1.0_r8-aH2O, 1.0e-4_r8 )

    ! calculate dry mass and dry volume of a bin
    mass_dry_a(ibin) = 0.0          ! initialize to 0.0
    vol_dry_a(ibin)  = 0.0          ! initialize to 0.0

    aer_H = (2.*aer(iso4_a,jtotal,ibin) +   &
         aer(ino3_a,jtotal,ibin) +   &
         aer(icl_a,jtotal,ibin)  +   &
         aer(imsa_a,jtotal,ibin) +   &
         2.*aer(ico3_a,jtotal,ibin))-   &
         (2.*aer(ica_a,jtotal,ibin)  +   &
         aer(ina_a,jtotal,ibin)  +   &
         aer(inh4_a,jtotal,ibin))
    aer_H = max(aer_H, 0.0d0)

    do iaer = 1, naer
       mass_dry_a(ibin) = mass_dry_a(ibin) +   &
            aer(iaer,jtotal,ibin)*mw_aer_mac(iaer)  ! ng/m^3(air)
       vol_dry_a(ibin)  = vol_dry_a(ibin) +   &
            aer(iaer,jtotal,ibin)*mw_aer_mac(iaer)/dens_aer_mac(iaer)       ! ncc/m^3(air)
    enddo
    mass_dry_a(ibin) = mass_dry_a(ibin) + aer_H
    vol_dry_a(ibin) = vol_dry_a(ibin) + aer_H

    mass_dry_a(ibin) = mass_dry_a(ibin)*1.e-15                      ! g/cc(air)
    vol_dry_a(ibin)  = vol_dry_a(ibin)*1.e-15                               ! cc(aer)/cc(air) or m^3/m^3(air)

    ! wet mass and wet volume
    mass_wet_a(ibin) = mass_dry_a(ibin) + water_a(ibin)*1.e-3               ! g/cc(air)
    vol_wet_a(ibin)  = vol_dry_a(ibin) + water_a(ibin)*1.e-3                ! cc(aer)/cc(air) or m^3/m^3(air)

    water_a_up(ibin) = aerosol_water_up(ibin,electrolyte,aer,kappa_nonelectro,a_zsr)   ! for hysteresis curve determination

    iter_kelvin = 0
    iter_kelvin_meth1 = 0

    iter_kelvin_state = 0
    if (iter_kelvin_method == 2) iter_kelvin_state = 2

    aH2O_a_old = aH2O
    kelvin_old = 1.0_r8
    rel_err_old = 1.0e30_r8
    rel_err_old2 = 1.0e30_r8
    water_a_old = 0.0_r8

    aH2O_a_del_state3 = 1.0e-3_r8
    aH2O_a_oldn = aH2O
    aH2O_a_oldp = aH2O
    kelvin_oldp = 1.0_r8
    kelvin_oldn = 1.0_r8
    rel_err_oldn = 1.0e30_r8
    rel_err_oldp = 1.0e30_r8
    water_a_oldp = 0.0_r8
    water_a_oldn = 0.0_r8
    aH2O_a_new = aH2O

10  iter_kelvin = iter_kelvin + 1
    aH2O_a(ibin) = aH2O_a_new

! RAZ uncommented the next 3 lines
      do je = 1, nelectrolyte
        molality0(je,ibin) = bin_molality(je,ibin,aH2O_a,b_zsr,a_zsr,aw_min)  ! compute aH2O dependent binary molalities  EFFI
      enddo
    call MESA( ibin, jaerosolstate, jphase, aer, jhyst_leg,         &
         electrolyte, epercent, activity, mc, num_a, mass_wet_a, mass_dry_a,              &
         mass_soluble_a, vol_dry_a, vol_wet_a, water_a, water_a_hyst, water_a_up, aH2O_a, &
         aH2O, ma, gam, log_gamZ, zc, za, gam_ratio, &
         xeq_a, na_Ma, nc_Mc, xeq_c, mw_electrolyte, mw_aer_mac, dens_aer_mac, Keq_ll,     &
         Keq_sl, MW_c, MW_a, growth_factor, MDRH, MDRH_T, molality0, rtol_mesa,            &
         jsalt_present, jsalt_index, jsulf_poor, jsulf_rich, phi_salt_old,    &
         kappa_nonelectro, mosaic_vars_aa )

    if(jaerosolstate(ibin) .eq. all_solid)then
       if (aer_pha_sta_diagaa >= 2) &
       write(*,'(a,5i5,2f12.8,1p,2e11.3)') 'aer_pha_sta_b', ibin, jhyst_leg(ibin), jaerosolstate(ibin), &
          iter_kelvin_state, iter_kelvin, aH2O, aH2O_a(ibin)
       return
    endif
    ! new wet mass and wet volume
    mass_wet_a(ibin) = mass_dry_a(ibin) + water_a(ibin)*1.e-3               ! g/cc(air)
    vol_wet_a(ibin)  = vol_dry_a(ibin) + water_a(ibin)*1.e-3                ! cc(aer)/cc(air) or m^3/m^3(air)

    call calculate_kelvin(ibin,num_a,vol_wet_a,aH2O_a,DpmV,kelvin,sigma_soln,T_K,  &
         sigma_water)
    !      kelvin(ibin) = 1.0
    kelvin(ibin) = max( kelvin(ibin), 1.0_r8 )
    if (water_a(ibin) <= 0.0_r8) kelvin(ibin) = 1.0_r8

    aH2O_a_new = aH2O/kelvin(ibin)

!   if(RH_pc .le. 97.0)then
!     rel_err = abs( (aH2O_a_new - aH2O_a(ibin))/aH2O_a(ibin))
!   else
!     if(water_a(ibin) .gt. 0.0)then
!       rel_err = abs( (water_a(ibin) - water_a_old)/water_a(ibin))
!     else
!       rel_err = 0.0 ! no soluble material is present
!     endif
!   endif
! the above rel_err involve differences between current and previous
!    iteration values, and is not suitable for bisection
! this rel_err below uses error from the exact solution, and is suitable for bisection
    rel_err = (aH2O_a(ibin)*kelvin(ibin) - aH2O) / max( aH2O, 0.01_r8 )

    if (aer_pha_sta_diagaa >= 10) &
    write(*,'(a,2i5, 1p,e10.2, 0p,f14.10, 2x,2f14.10, 2x,1p,2e18.10)') &
       'iter_kelvin', iter_kelvin_state, iter_kelvin, rel_err, kelvin(ibin), &
       aH2O_a(ibin), aH2O_a_new, water_a_old, water_a(ibin)

    if (abs(rel_err) <= kelvin_toler) then
       iter_kelvin_state = iter_kelvin_state + 100
       goto 90
    end if

    if (iter_kelvin_state <= 0) then
       ! doing rahuls original iteration method
       itmpa = 0
       if (iter_kelvin >= iter_kelvin_meth1_max) then
          itmpa = 1
       else if (iter_kelvin >= iter_kelvin_meth1_max) then
          tmpa = min( rel_err_old, rel_err_old2 )
          if (tmpa < 0.0_r8 .and. rel_err <= tmpa) itmpa = 1
          tmpa = max( rel_err_old, rel_err_old2 )
          if (tmpa > 0.0_r8 .and. rel_err >= tmpa) itmpa = 1
       end if

       if (itmpa > 0) then
          if (iter_kelvin_method <= 1) then
             ! quit if number of iterations is too large OR
             ! rel_err is outside the range of the previous two rel_err values,
             !    and one previous rel_err is positive, and one previous rel_err is negative
             aH2O_a(ibin) = aH2O_a_new   ! do this to get same output as prev version
             if (aer_pha_sta_diagaa >= 1) &
             write(*,'(a,5i5,2f12.8,1p,3e11.3)') 'iter_kelv_err', ibin, jhyst_leg(ibin), jaerosolstate(ibin), &
                iter_kelvin_state, iter_kelvin, aH2O, aH2O_a(ibin), rel_err, kelvin_toler
             iter_kelvin_state = 100
             goto 90
          else
             ! switch to method 2 but do not iterate yet
             iter_kelvin_state = 1
             iter_kelvin_meth1 = iter_kelvin
          end if
       else
          ! save current values to old then do next iteration
          aH2O_a_old = aH2O_a(ibin)
          kelvin_old = kelvin(ibin)
          rel_err_old2 = rel_err_old
          rel_err_old = rel_err
          water_a_old  = water_a(ibin)
       !        aH2O = aH2O_a_new
       !        call MTEM_compute_log_gamZ  ! recompute activity coeffs (for surface tension and solid-liquid equilibria)
          goto 10
       end if
    endif

    if (iter_kelvin_state == 1) then
       ! rahuls original iteration method failed, so do some things before switching to bisection
       iter_kelvin_state = 2
       if (rel_err < 0.0_r8) then
          ! current aH2O_a has negative rel_err so must start at the beginning
          aH2O_a_new = aH2O
          goto 10
       else
          ! current aH2O_a has positive rel_err and can be used in bisection
          ! do not iterate yet
          continue
       end if
    end if

    if (iter_kelvin_state == 2) then
       ! this is first "setup" step of bisection, and the algorithm is expecting that
       !    the current aH2O_a has hel_err be > 0, and can be used as one of the 2 bisection points
       if (rel_err < 0.0_r8) then
          ! error should be positive, so this is a fatal error
          if (aer_pha_sta_diagaa >= 1) &
             write(*,'(a,5i5,2f12.8,1p,3e11.3)') 'iter_kelv_er2', ibin, jhyst_leg(ibin), jaerosolstate(ibin), &
                iter_kelvin_state, iter_kelvin, aH2O, aH2O_a(ibin), rel_err, kelvin_toler
          iter_kelvin_state = 100
          goto 90
       end if
       ! current aH2O_a will work as one of the two initial bisection points
       ! (the one with a positive error)
       aH2O_a_oldp = aH2O_a(ibin)
       kelvin_oldp = kelvin(ibin)
       rel_err_oldp = rel_err
       water_a_oldp  = water_a(ibin)
       aH2O_a_new = min( aH2O/kelvin(ibin), 0.999999_r8 )   ! is this needed, or should it be 1.0, or ???
       iter_kelvin_state = 3
       goto 10
    end if

    if (iter_kelvin_state == 3) then
       ! this is the second "setup" step of bisection, and the algorithm is looking for an aH2O_a
       ! that has rel_err < 0, so that the "root" will be bracketed and bisection can begin
       if (rel_err < 0.0_r8) then
          ! current aH2O_a will work as one of the two initial bisection points
          ! (the one with a negative error)
          aH2O_a_oldn = aH2O_a(ibin)
          kelvin_oldn = kelvin(ibin)
          rel_err_oldn = rel_err
          water_a_oldn  = water_a(ibin)
          aH2O_a_new = 0.5_r8*(aH2O_a_oldn + aH2O_a_oldp)
          iter_kelvin_state = 4
          goto 10
       else
          ! need to find a point with a negative error
          if ( (rel_err >= rel_err_oldp) .or. &
               (aH2O_a_del_state3 >= 0.999_r8) ) then
             ! cannot find such a point -- this is a fatal error
             if (aer_pha_sta_diagaa >= 1) &
                write(*,'(a,5i5,2f12.8,1p,3e11.3)') 'iter_kelv_er3', ibin, jhyst_leg(ibin), jaerosolstate(ibin), &
                   iter_kelvin_state, iter_kelvin, aH2O, aH2O_a(ibin), rel_err, kelvin_toler
             iter_kelvin_state = 200
             goto 90
          else
             ! save current aH2O_a as the initial bisection point with positive error
             ! then calc aH2O_a_new = aH2O_a(ibin) - aH2O_a_del_state3
             ! which will hopefully have a negative error
             aH2O_a_oldp = aH2O_a(ibin)
             kelvin_oldp = kelvin(ibin)
             rel_err_oldp = rel_err
             water_a_oldp  = water_a(ibin)
             aH2O_a_new = aH2O_a(ibin) - aH2O_a_del_state3
             aH2O_a_del_state3 = aH2O_a_del_state3*1.5_r8
             if (aH2O_a_new .le. 0.01_r8) then
                aH2O_a_new = 0.01_r8
                aH2O_a_del_state3 = 1.0_r8
             end if
             goto 10
          end if
       end if
    end if

    if (iter_kelvin_state == 4) then
       ! at this point, the algorithm is doing bisection
       if ( iter_kelvin >= iter_kelvin_meth2_max + iter_kelvin_meth1 ) then
          ! maximum iterations is exceeded
          if (aer_pha_sta_diagaa >= 1) &
             write(*,'(a,5i5,2f12.8,1p,3e11.3)') 'iter_kelv_er4', ibin, jhyst_leg(ibin), jaerosolstate(ibin), &
                iter_kelvin_state, iter_kelvin, aH2O, aH2O_a(ibin), rel_err, kelvin_toler
          iter_kelvin_state = 301
          goto 90
       else if ( abs(aH2O_a_oldp - aH2O_a_oldn) <= aH2O_range_bisect_toler ) then
          ! the aH2O_a_oldp to aH2O_a_oldn range is very small, which is treated as convergence
!         if (aer_pha_sta_diagaa >= 1) &
!            write(*,'(a,5i5,2f12.8,1p,4e11.3)') 'iter_kelv_er5', ibin, jhyst_leg(ibin), jaerosolstate(ibin), &
!               iter_kelvin_state, iter_kelvin, aH2O, aH2O_a(ibin), rel_err, kelvin_toler, &
!               aH2O_range_bisect_toler
          iter_kelvin_state = 302
          goto 90
       end if
       ! decide if the current aH2O_a should replace the old negative-error point
       !    or the old positive-error point
       if (rel_err >= 0.0_r8) then
          if (rel_err >= rel_err_oldp) then
             ! current aH2O_a has positive error, but the error is not smaller
             !    than the old positive-error point -- this is a fatal error
             if (aer_pha_sta_diagaa >= 1) &
                write(*,'(a,5i5,2f12.8,1p,3e11.3)') 'iter_kelv_er6', ibin, jhyst_leg(ibin), jaerosolstate(ibin), &
                   iter_kelvin_state, iter_kelvin, aH2O, aH2O_a(ibin), rel_err, kelvin_toler
             iter_kelvin_state = 303
             goto 90
          else
             ! current aH2O_a has positive error and replaces the the old positive-error point
             aH2O_a_oldp = aH2O_a(ibin)
             kelvin_oldp = kelvin(ibin)
             rel_err_oldp = rel_err
             water_a_oldp  = water_a(ibin)
          end if
       else
          if (rel_err <= rel_err_oldn) then
             ! current aH2O_a has negative error, but the error is not smaller
             !    than the old negative-error point -- this is a fatal error
             if (aer_pha_sta_diagaa >= 1) &
                write(*,'(a,5i5,2f12.8,1p,3e11.3)') 'iter_kelv_er7', ibin, jhyst_leg(ibin), jaerosolstate(ibin), &
                   iter_kelvin_state, iter_kelvin, aH2O, aH2O_a(ibin), rel_err, kelvin_toler
             iter_kelvin_state = 304
             goto 90
          else
             ! current aH2O_a has negative error and replaces the the old negative-error point
             aH2O_a_oldn = aH2O_a(ibin)
             kelvin_oldn = kelvin(ibin)
             rel_err_oldn = rel_err
             water_a_oldn  = water_a(ibin)
          end if
       end if
       aH2O_a_new = 0.5_r8*(aH2O_a_oldn + aH2O_a_oldp)
       goto 10
    end if

    write(errmsg,'(a,4i5)') 'iter_kelv fatal err 1', ibin, iter_kelvin, iter_kelvin_state
    call wrf_error_fatal(trim(adjustl(errmsg)))

    ! kelvin iterations completed
90  if (iter_kelvin_state == 200) then
       ! select aH2O_a(ibin) or aH2O_a_oldp, whichever has lowest error
       if (abs(rel_err_oldp) < abs(rel_err)) then
          aH2O_a(ibin) = aH2O_a_oldp
          rel_err = rel_err_oldp
       end if
    else if (iter_kelvin_state >= 300 .and. iter_kelvin_state <= 304) then
       ! select aH2O_a(ibin) or aH2O_a_oldp or aH2O_a_oldn, whichever has lowest error
       tmpa = min( abs(rel_err_oldn), abs(rel_err_oldp), abs(rel_err) )
       if (abs(rel_err_oldp) == tmpa) then
          aH2O_a(ibin) = aH2O_a_oldp
          rel_err = rel_err_oldp
       else if (abs(rel_err_oldn) == tmpa) then
          aH2O_a(ibin) = aH2O_a_oldn
          rel_err = rel_err_oldn
       end if
    end if

    if(jaerosolstate(ibin) .eq. all_liquid)jhyst_leg(ibin) = jhyst_up

    ! now compute kelvin effect terms for condensing species (nh3, hno3, and hcl)
    do iv = 1,  ngas_aerchtot
       term = 4.*sigma_soln(ibin)*partial_molar_vol(iv)/   &
            (8.3144e7*T_K*DpmV(ibin))
       kel(iv,ibin) = 1. + term*(1. + 0.5*term*(1. + term/3.))
    enddo

    if (aer_pha_sta_diagaa >= 2) &
    write(*,'(a,5i5,2f12.8,1p,e11.3,e14.5)') 'aer_pha_sta_c', ibin, jhyst_leg(ibin), jaerosolstate(ibin), &
       iter_kelvin_state, iter_kelvin, aH2O, aH2O_a(ibin), rel_err, water_a(ibin)
    return
  end subroutine aerosol_phase_state

  !**********************************`*************************************
  ! MESA: Multicomponent Equilibrium Solver for Aerosols.
  ! Computes equilibrum solid and liquid phases by integrating
  ! pseudo-transient dissolution and precipitation reactions
  !
  ! author: Rahul A. Zaveri
  ! update: sep 2015
  !
  ! 9/3/2015 RAZ: Bugfix - fixed phase state calculations for aerosols that dont contain any salts,
  !               but can still contain water due to presence of BC, OC, SOA, and OIN, which are now
  !               allowed to absorb some water.
  !-----------------------------------------------------------------------
  subroutine MESA( ibin, jaerosolstate, jphase, aer, jhyst_leg,      &
       electrolyte, epercent, activity, mc, num_a, mass_wet_a, mass_dry_a, mass_soluble_a,    &
       vol_dry_a, vol_wet_a, water_a, water_a_hyst, water_a_up, aH2O_a, aH2O,                &
       ma, gam, log_gamZ, zc, za, gam_ratio, xeq_a,     &
       na_Ma, nc_Mc, xeq_c, mw_electrolyte, mw_aer_mac, dens_aer_mac, Keq_ll, Keq_sl, MW_c,    &
       MW_a, growth_factor, MDRH, MDRH_T, molality0, rtol_mesa, jsalt_present, jsalt_index,   &
       jsulf_poor, jsulf_rich, phi_salt_old,                                       &
       kappa_nonelectro, mosaic_vars_aa )

    use module_data_mosaic_aero,  only: r8, nbin_a_max, nelectrolyte, Ncation, naer,        &!Parameters
         jtotal, all_solid, jsolid, all_liquid, jliquid, jhyst_lo, mhyst_uporlo_jhyst,       &!Parameters
         jhyst_up, mhyst_uporlo_waterhyst, nsoluble, nsalt, Nanion, nrxn_aer_sl,            &
         nrxn_aer_ll, MDRH_T_NUM, jsulf_poor_NUM, jsulf_rich_NUM,                         &!Parameters
         ptol_mol_astem,  mhyst_force_lo,  mhyst_force_up,                               &!Input
         jcacl2, jcano3, mhyst_method, ioin_a, ibc_a, jcaco3, jcaso4,                     & !TBD
         mosaic_vars_aa_type

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    integer, intent(inout), dimension(nbin_a_max)  :: jhyst_leg
    integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate,jphase
    integer, intent(in), dimension(nsalt) :: jsalt_index
    integer, intent(inout), dimension(nsalt) :: jsalt_present
    integer, intent(in), dimension(jsulf_poor_NUM) :: jsulf_poor
    integer, intent(in), dimension(jsulf_rich_NUM) :: jsulf_rich

    real(r8), intent(in) :: aH2O,rtol_mesa
    real(r8), intent(in), dimension(naer) :: mw_aer_mac,dens_aer_mac
    real(r8), intent(in), dimension(Ncation) :: zc,MW_c
    real(r8), intent(inout), dimension(Ncation) :: nc_Mc,xeq_c
    real(r8), intent(in), dimension(Nanion)  :: za,MW_a
    real(r8), intent(inout), dimension(Nanion)  :: xeq_a,na_Ma
    real(r8), intent(inout), dimension(nbin_a_max) :: num_a,mass_wet_a,mass_dry_a
    real(r8), intent(inout), dimension(nbin_a_max) :: mass_soluble_a,vol_dry_a
    real(r8), intent(inout), dimension(nbin_a_max) :: vol_wet_a,gam_ratio
    real(r8), intent(inout), dimension(nbin_a_max) :: water_a,water_a_hyst,water_a_up
    real(r8), intent(inout), dimension(nbin_a_max) :: aH2O_a,growth_factor,MDRH
    real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
    real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
    real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
    real(r8), intent(inout), dimension(nrxn_aer_sl) :: Keq_sl
    real(r8), intent(inout), dimension(MDRH_T_NUM) :: MDRH_T
    real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
    real(r8), intent(inout), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
    real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
    real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: epercent
    real(r8), intent(inout), dimension(nsalt) :: phi_salt_old
    real(r8), intent(in), dimension(naer) :: kappa_nonelectro

    type (mosaic_vars_aa_type), intent(inout) :: mosaic_vars_aa

    ! local variables
    integer :: idissolved, j_index, jsalt_dum, jdum, js, je ! 9/3/2015 RAZ: added jsalt_dum
    real(r8) :: CRH, solids, sum_soluble, sum_insoluble, XT !BALLI** XT, should it be subr arg?
    !real(r8) :: aerosol_water                               ! mosaic func
    !real(r8) :: drh_mutual                          ! mosaic func
    real(r8) :: H_ion, sum_dum

    !! EFFI
    !! calculate percent composition
    sum_dum = 0.0
    do je = 1, nelectrolyte
       sum_dum = sum_dum + electrolyte(je,jtotal,ibin)
    enddo

    if(sum_dum .eq. 0.)sum_dum = 1.0

    do je = 1, nelectrolyte
       epercent(je,jtotal,ibin) = 100.*electrolyte(je,jtotal,ibin)/sum_dum
    enddo

    call calculate_XT(ibin,jtotal,XT,aer)

!! begin new algorithm - 6/3/2015 RAZ
    jsalt_dum = 0 ! 9/3/2015 RAZ
    do js = 1, nsalt
       jsalt_present(js) = 0                        ! default value - salt absent

       if(epercent(js,jtotal,ibin) .gt. ptol_mol_astem)then
          jsalt_present(js) = 1                     ! salt present
          jsalt_dum = jsalt_dum + jsalt_index(js)   ! 9/3/2015 RAZ
       endif
    enddo

    if( (epercent(jcano3,jtotal,ibin) .gt. ptol_mol_astem) .or. &
        (epercent(jcacl2,jtotal,ibin) .gt. ptol_mol_astem) )then
      CRH = 0.0  ! no crystrallization or efflorescence point
    else
      CRH = 0.35 ! default value
    endif

  ! now diagnose MDRH
    if(jsalt_dum .eq. 0)then			! no salts or acids are present ! 9/3/2015 RAZ: updated algorithm for jsalt_dum = 0

       CRH = 0.0
       MDRH(ibin) = 0.0
       jaerosolstate(ibin) = all_solid
       jphase(ibin)    = jsolid
       jhyst_leg(ibin) = jhyst_lo
       call adjust_solid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent,water_a)
       water_a(ibin) = aerosol_water(jtotal,ibin,jaerosolstate,jphase,jhyst_leg,   &	! 9/3/2015 RAZ: water due to nonelectrolytes (OC, BC, SOA, OIN)
        electrolyte,aer,kappa_nonelectro,num_a,mass_dry_a,mass_soluble_a,aH2O,molality0)
       return

    elseif(XT .lt. 1. .and. XT .gt. 0.0)then  	! excess sulfate, always liquid, MDRH=0.0
       MDRH(ibin) = 0.0
    elseif(XT .ge. 2.0 .or. XT .lt. 0.0)then  	! sulfate poor
       j_index = jsulf_poor(jsalt_dum)          ! 9/3/2015 RAZ
       MDRH(ibin) = MDRH_T(j_index)
    else					! sulfate rich
       j_index = jsulf_rich(jsalt_dum)          ! 9/3/2015 RAZ
       MDRH(ibin) = MDRH_T(j_index)
    endif

    CRH = min(CRH, MDRH(ibin)/100.0)		! 6/3/2015 RAZ

!! end new algorithm - 6/3/2015 RAZ

    ! modified step 1: 9/3/2015 RAZ
    ! step 1: check if aH2O is below CRH (crystallization or efflorescence point)
    if( aH2O_a(ibin).lt.CRH )then
       jaerosolstate(ibin) = all_solid
       jphase(ibin)    = jsolid
       jhyst_leg(ibin) = jhyst_lo
       call adjust_solid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent,water_a)
       water_a(ibin) = aerosol_water(jtotal,ibin,jaerosolstate,jphase,jhyst_leg,   &	! 9/3/2015 RAZ: water due to nonelectrolytes (OC, BC, SOA, OIN)
        electrolyte,aer,kappa_nonelectro,num_a,mass_dry_a,mass_soluble_a,aH2O,molality0)
       return
    endif

    ! step 2: check mhyst_method for supersaturation/metastable state
    jdum = 0
    if (mhyst_method == mhyst_uporlo_jhyst) then         ! BOX method/logic
       if (jhyst_leg(ibin) == jhyst_up) jdum = 1
    elseif (mhyst_method == mhyst_uporlo_waterhyst) then ! 3-D method/logic
       if (water_a_hyst(ibin) > 0.5*water_a_up(ibin)) jdum = 1
       !BSINGH - 05/28/2013(RCE updates)
    elseif (mhyst_method == mhyst_force_lo) then
       jdum = 0
    elseif (mhyst_method == mhyst_force_up) then
       jdum = 1
       !BSINGH - 05/28/2013(RCE updates ENDS)
    else
       call wrf_error_fatal('*** MESA - bad mhyst_method')
    endif
    if (jdum == 1) then ! the aerosol is fully deliquesced in metastable or subsaturated state
       call do_full_deliquescence(ibin,aer,electrolyte)

       !        call ions_to_electrolytes(jliquid,ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma,nc_Mc,xeq_c) ! for Li and Lu surface tension
       !        call compute_activities(ibin,jphase,aer,jhyst_leg,electrolyte, &
       !activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a,aH2O,ma,gam,log_gamZ,gam_ratio)              ! for Li and Lu surface tension

! MODIFIED LOGIC IF SOA, POA, BC, OIN ARE ASSUMED TO BE SLIGHTLY HYGROSCOPIC  RAZ 4/16/2014
!       sum_soluble = 0.0
!       do js = 1, nsoluble
!          sum_soluble = sum_soluble + electrolyte(js,jtotal,ibin)
!       enddo
!
!       solids = electrolyte(jcaso4,jtotal,ibin) +   &
!                electrolyte(jcaco3,jtotal,ibin) +   &
!                aer(ioin_a,jtotal,ibin)         +   &
!                aer(ibc_a,jtotal,ibin)
!
!
!       if(sum_soluble .le. 0.0 .and. solids .gt. 0.0)then ! RAZ modified logic
!
!          jdum = 0
!          jaerosolstate(ibin) = all_solid ! no soluble material present, so go back to solid state
!          jphase(ibin) = jsolid
!          call adjust_solid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent,water_a)
!
!          ! new wet mass and wet volume
!          mass_wet_a(ibin) = mass_dry_a(ibin) + water_a(ibin)*1.e-3 ! g/cc(air)
!          vol_wet_a(ibin)  = vol_dry_a(ibin) + water_a(ibin)*1.e-3  ! cc(aer)/cc(air) or m^3/m^3(air)
!          growth_factor(ibin) = mass_wet_a(ibin)/mass_dry_a(ibin)   ! mass growth factor
!
!          return
!
!       elseif(sum_soluble .gt. 0.0)then  ! RAZ modified logic
!
          jaerosolstate(ibin) = all_liquid
          jhyst_leg(ibin) = jhyst_up
          jphase(ibin) = jliquid
          water_a(ibin) = aerosol_water(jtotal,ibin,jaerosolstate,jphase,jhyst_leg,   &
               electrolyte,aer,kappa_nonelectro,num_a,mass_dry_a,mass_soluble_a,aH2O,molality0)
          if(water_a(ibin) .le. 0.0)then     ! one last attempt to catch bad input
             jdum = 0
             jaerosolstate(ibin) = all_solid ! no soluble material present
             jphase(ibin)    = jsolid
             jhyst_leg(ibin) = jhyst_lo
             call adjust_solid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent,water_a)
          else
             call adjust_liquid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent)
             call compute_activities(ibin,jaerosolstate,jphase,aer,jhyst_leg,      &
                  electrolyte,activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a, &
                  aH2O,ma,gam,log_gamZ,gam_ratio,Keq_ll,molality0,kappa_nonelectro)
          endif

          ! new wet mass and wet volume
          mass_wet_a(ibin) = mass_dry_a(ibin) + water_a(ibin)*1.e-3 ! g/cc(air)
          vol_wet_a(ibin)  = vol_dry_a(ibin) + water_a(ibin)*1.e-3  ! cc(aer)/cc(air) or m^3/m^3(air)
          growth_factor(ibin) = mass_wet_a(ibin)/mass_dry_a(ibin)   ! mass growth factor

          return

!       endif

    endif ! jdum

    ! step 3: diagnose phase state based on MDRH
    if(aH2O_a(ibin)*100. .lt. MDRH(ibin)) then
       jaerosolstate(ibin) = all_solid
       jphase(ibin) = jsolid
       jhyst_leg(ibin) = jhyst_lo
       call adjust_solid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent,water_a)
       return
    endif

    ! step 4: none of the above means it must be sub-saturated or mixed-phase
10  call do_full_deliquescence(ibin,aer,electrolyte)
    call MESA_PTC( ibin, jaerosolstate, jphase, aer, jhyst_leg,                   &
         electrolyte, epercent, activity, mc, num_a, mass_dry_a, mass_wet_a,      &
         mass_soluble_a, vol_dry_a, vol_wet_a, water_a, aH2O,                     &
         ma, gam, log_gamZ, zc, za, gam_ratio, xeq_a, na_Ma, nc_Mc, xeq_c,        &
         mw_electrolyte, mw_aer_mac, dens_aer_mac, Keq_sl, MW_c, MW_a, Keq_ll,    &
         growth_factor, molality0, rtol_mesa, jsalt_present, phi_salt_old,        &
         kappa_nonelectro, mosaic_vars_aa                                    )     ! determines jaerosolstate(ibin)
    return
  end subroutine MESA

  !***********************************************************************
  ! computes kelvin effect term (kelvin => 1.0)
  !
  ! author: Rahul A. Zaveri
  ! update: jan 2005
  !-----------------------------------------------------------------------
  subroutine calculate_kelvin(ibin,num_a,vol_wet_a,aH2O_a,DpmV,kelvin,sigma_soln,  &
       T_K,sigma_water)
    use module_data_mosaic_constants, only:  pi
    use module_data_mosaic_aero, only: r8,nbin_a_max                                   !Parameters

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    real(r8), intent(in) :: T_K,sigma_water
    real(r8), intent(in), dimension(nbin_a_max) :: num_a
    real(r8), intent(inout), dimension(nbin_a_max) :: sigma_soln
    real(r8), intent(inout), dimension(nbin_a_max) ::vol_wet_a,aH2O_a,DpmV,kelvin
    ! local variables
    integer je
    real(r8) :: term, sum_dum
    real(r8), dimension(nbin_a_max) :: volume_a

    volume_a(ibin) = vol_wet_a(ibin)                                ! [cc/cc(air)]
    DpmV(ibin)=(6.*volume_a(ibin)/(num_a(ibin)*pi))**(1./3.)        ! [cm]

    ! Li and Lu (2001) surface tension model:
    !      sum_dum = 0.0
    !      do je = 1, nelectrolyte
    !        sum_dum = sum_dum + G_MX(je)*
    !     &                      alog(1./(1.+K_MX(je)*activity(je,ibin)))
    !      enddo
    !      sigma_soln(ibin) = sigma_water + 8.3144e7*T_K*sum_dum

    ! simpler correlation for solution surface tension:
    sigma_soln(ibin) = sigma_water + 49.0*(1. - aH2O_a(ibin))       ! [dyn/cm]

    term = 72.*sigma_soln(ibin)/(8.3144e7*T_K*DpmV(ibin))           ! [-]
!    kelvin(ibin) = exp(term)
    kelvin(ibin) = 1. + term*(1. + 0.5*term*(1. + term/3.))

    return
  end subroutine calculate_kelvin

  !***********************************************************************
  ! computes sulfate ratio
  !
  ! author: Rahul A. Zaveri
  ! update: dec 1999
  !-----------------------------------------------------------------------
  subroutine calculate_XT(ibin,jp,XT,aer)
    use module_data_mosaic_aero, only: r8,naer,nbin_a_max,                         &
         imsa_a,iso4_a,ica_a,ina_a,inh4_a

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin, jp
    real(r8), intent(inout) :: XT
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer

    if( (aer(iso4_a,jp,ibin)+aer(imsa_a,jp,ibin)) .gt.0.0)then
       XT   = ( aer(inh4_a,jp,ibin) +   &
            aer(ina_a,jp,ibin)  +   &
            2.*aer(ica_a,jp,ibin) )/   &
            (aer(iso4_a,jp,ibin)+0.5*aer(imsa_a,jp,ibin))
    else
       XT   = -1.0
    endif

    return
  end subroutine calculate_XT

  !***********************************************************************
  ! called when aerosol bin is completely solid.
  !
  ! author: Rahul A. Zaveri
  ! update: jan 2005
  !-----------------------------------------------------------------------
  subroutine adjust_solid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent,  &
       water_a)

    use module_data_mosaic_aero, only: r8,nbin_a_max,naer,nelectrolyte,jsolid,     &!Parameters
         jhyst_lo,jtotal,jliquid,                                                  &!Parameters
         inh4_a,ino3_a,icl_a                                                        !TBD

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    integer, intent(inout), dimension(nbin_a_max) :: jphase,jhyst_leg
    real(r8), intent(inout), dimension(nbin_a_max) :: water_a
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte,epercent
    ! local variables
    integer iaer, je

    jphase(ibin)    = jsolid
    jhyst_leg(ibin) = jhyst_lo   ! lower curve
    water_a(ibin)   = 0.0

    ! transfer aer(jtotal) to aer(jsolid)
    do iaer = 1, naer
       aer(iaer, jsolid, ibin) = aer(iaer,jtotal,ibin)
       aer(iaer, jliquid,ibin) = 0.0
    enddo

    ! transfer electrolyte(jtotal) to electrolyte(jsolid)
    do je = 1, nelectrolyte
       electrolyte(je,jliquid,ibin) = 0.0
       epercent(je,jliquid,ibin)    = 0.0
       electrolyte(je,jsolid,ibin)  = electrolyte(je,jtotal,ibin)
       epercent(je,jsolid,ibin)     = epercent(je,jtotal,ibin)
    enddo

    ! update aer(jtotal) that may have been affected above
    aer(inh4_a,jtotal,ibin) = aer(inh4_a,jsolid,ibin)
    aer(ino3_a,jtotal,ibin) = aer(ino3_a,jsolid,ibin)
    aer(icl_a,jtotal,ibin)  = aer(icl_a,jsolid,ibin)

    return
  end subroutine adjust_solid_aerosol

  !***********************************************************************
  ! called when aerosol bin is completely liquid.
  !
  ! author: Rahul A. Zaveri
  ! update: jan 2005
  !-----------------------------------------------------------------------
  subroutine adjust_liquid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent) ! TOUCH

    use module_data_mosaic_aero, only: r8,nbin_a_max,naer,nelectrolyte,jliquid,    &!Parameters
         jhyst_up,jsolid,jtotal,                                                   &!Parameters
         jcaco3,jcaso4,                                                            &
         inh4_a,ina_a,ica_a,imsa_a,icl_a,ino3_a,iso4_a

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    integer, intent(inout), dimension(nbin_a_max) :: jphase,jhyst_leg

    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: epercent
    ! local variables
    integer :: iaer, je

    jphase(ibin)    = jliquid
    jhyst_leg(ibin) = jhyst_up   ! upper curve

    ! partition all electrolytes into liquid phase
    do je = 1, nelectrolyte
       electrolyte(je,jsolid,ibin)  = 0.0
       epercent(je,jsolid,ibin)     = 0.0
       electrolyte(je,jliquid,ibin) = electrolyte(je,jtotal,ibin)
       epercent(je,jliquid,ibin)    = epercent(je,jtotal,ibin)
    enddo
    ! except these electrolytes, which always remain in the solid phase
    electrolyte(jcaco3,jsolid,ibin) = electrolyte(jcaco3,jtotal,ibin)
    electrolyte(jcaso4,jsolid,ibin) = electrolyte(jcaso4,jtotal,ibin)
    epercent(jcaco3,jsolid,ibin)    = epercent(jcaco3,jtotal,ibin)
    epercent(jcaso4,jsolid,ibin)    = epercent(jcaso4,jtotal,ibin)
    electrolyte(jcaco3,jliquid,ibin)= 0.0
    electrolyte(jcaso4,jliquid,ibin)= 0.0
    epercent(jcaco3,jliquid,ibin)   = 0.0
    epercent(jcaso4,jliquid,ibin)   = 0.0

    ! partition all the aer species into
    ! solid phase
    do iaer = 1, naer
    aer(iaer,jsolid,ibin)  = aer(iaer,jtotal,ibin)
    end do
    aer(iso4_a,jsolid,ibin) = electrolyte(jcaso4,jsolid,ibin)
    aer(ino3_a,jsolid,ibin) = 0.0
    aer(icl_a,jsolid,ibin)  = 0.0
    aer(inh4_a,jsolid,ibin) = 0.0
    aer(imsa_a,jsolid,ibin) = 0.0
    aer(ina_a,jsolid,ibin)  = 0.0
    aer(ica_a,jsolid,ibin)  = electrolyte(jcaco3,jsolid,ibin) +   &
                              electrolyte(jcaso4,jsolid,ibin)
!   aer(ico3_a,jsolid,ibin) = aer(ico3_a,jtotal,ibin)
!   aer(ioc_a,jsolid,ibin)  = aer(ioc_a,jtotal,ibin)
!   aer(ibc_a,  jsolid,ibin)= aer(ibc_a,  jtotal,ibin)
!   aer(ioin_a, jsolid,ibin)= aer(ioin_a, jtotal,ibin)
!   aer(iaro1_a,jsolid,ibin)= aer(iaro1_a,jtotal,ibin)
!   aer(iaro2_a,jsolid,ibin)= aer(iaro2_a,jtotal,ibin)
!   aer(ialk1_a,jsolid,ibin)= aer(ialk1_a,jtotal,ibin)
!   aer(iole1_a,jsolid,ibin)= aer(iole1_a,jtotal,ibin)
!   aer(iapi1_a,jsolid,ibin)= aer(iapi1_a,jtotal,ibin)
!   aer(iapi2_a,jsolid,ibin)= aer(iapi2_a,jtotal,ibin)
!   aer(ilim1_a,jsolid,ibin)= aer(ilim1_a,jtotal,ibin)
!   aer(ilim2_a,jsolid,ibin)= aer(ilim2_a,jtotal,ibin)

    ! liquid-phase
    do iaer = 1, naer
    aer(iaer,jliquid,ibin)  = 0.0
    end do
    aer(iso4_a,jliquid,ibin) = aer(iso4_a,jtotal,ibin) -   &
                               aer(iso4_a,jsolid,ibin)
    aer(iso4_a,jliquid,ibin) = max(0.d0, aer(iso4_a,jliquid,ibin)) ! RAZ 4/16/2014
    aer(ino3_a,jliquid,ibin) = aer(ino3_a,jtotal,ibin)
    aer(icl_a,jliquid,ibin)  = aer(icl_a,jtotal,ibin)
    aer(inh4_a,jliquid,ibin) = aer(inh4_a,jtotal,ibin)
    aer(imsa_a,jliquid,ibin) = aer(imsa_a,jtotal,ibin)
    aer(ina_a,jliquid,ibin)  = aer(ina_a,jtotal,ibin)
    aer(ica_a,jliquid,ibin)  = aer(ica_a,jtotal,ibin) -   &
                               aer(ica_a,jsolid,ibin)
    aer(ica_a,jliquid,ibin)  = max(0.d0, aer(ica_a,jliquid,ibin)) ! RAZ 4/16/2014
!   aer(ico3_a,jliquid,ibin) = 0.0
!   aer(ioc_a,jliquid,ibin)  = 0.0
!   aer(ibc_a,  jliquid,ibin)= 0.0
!   aer(ioin_a, jliquid,ibin)= 0.0
!   aer(iaro1_a,jliquid,ibin)= 0.0
!   aer(iaro2_a,jliquid,ibin)= 0.0
!   aer(ialk1_a,jliquid,ibin)= 0.0
!   aer(iole1_a,jliquid,ibin)= 0.0
!   aer(iapi1_a,jliquid,ibin)= 0.0
!   aer(iapi2_a,jliquid,ibin)= 0.0
!   aer(ilim1_a,jliquid,ibin)= 0.0
!   aer(ilim2_a,jliquid,ibin)= 0.0

    return
  end subroutine adjust_liquid_aerosol

  !***********************************************************************
  ! this subroutine completely deliquesces an aerosol and partitions
  ! all the soluble electrolytes into the liquid phase and insoluble
  ! ones into the solid phase. It also calculates the corresponding
  ! aer(js,jliquid,ibin) and aer(js,jsolid,ibin) generic species
  ! concentrations
  !
  ! author: Rahul A. Zaveri
  ! update: jan 2005
  !-----------------------------------------------------------------------
  subroutine do_full_deliquescence(ibin,aer,electrolyte)    ! TOUCH
    use module_data_mosaic_aero, only: r8,naer,nbin_a_max,nelectrolyte,jtotal,jsolid, &!Parameters
         jliquid,                                                                     &!Parameters
         jcacl2,jcano3,ioin_a,jcaco3,jcaso4,                                          &
         inh4_a,ina_a,ica_a,imsa_a,icl_a,ino3_a,iso4_a

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    ! local variables
    integer ::  iaer, js

    ! partition all electrolytes into liquid phase
    do js = 1, nelectrolyte
       electrolyte(js,jsolid,ibin)  = 0.0
       electrolyte(js,jliquid,ibin) = electrolyte(js,jtotal,ibin)
    enddo
    !
    ! except these electrolytes, which always remain in the solid phase
    electrolyte(jcaco3,jsolid,ibin) = electrolyte(jcaco3,jtotal,ibin)
    electrolyte(jcaso4,jsolid,ibin) = electrolyte(jcaso4,jtotal,ibin)
    electrolyte(jcaco3,jliquid,ibin)= 0.0
    electrolyte(jcaso4,jliquid,ibin)= 0.0

    ! partition all the generic aer species into solid and liquid phases
    ! solid phase
    do iaer = 1, naer
    aer(iaer,jsolid,ibin)= aer(iaer,jtotal,ibin)
    end do
    aer(iso4_a,jsolid,ibin) = electrolyte(jcaso4,jsolid,ibin)
    aer(ino3_a,jsolid,ibin) = 0.0
    aer(icl_a, jsolid,ibin) = 0.0
    aer(inh4_a,jsolid,ibin) = 0.0
    aer(imsa_a,jsolid,ibin) = 0.0
    aer(ina_a, jsolid,ibin) = 0.0
    aer(ica_a, jsolid,ibin) = electrolyte(jcaco3,jsolid,ibin) +   &
                              electrolyte(jcaso4,jsolid,ibin)
!   aer(ico3_a,jsolid,ibin) = aer(ico3_a,jtotal,ibin)
!   aer(ioc_a, jsolid,ibin) = aer(ioc_a,jtotal,ibin)
!   aer(ibc_a, jsolid,ibin) = aer(ibc_a,jtotal,ibin)
!   aer(ioin_a,jsolid,ibin) = aer(ioin_a,jtotal,ibin)
!   aer(iaro1_a,jsolid,ibin)= aer(iaro1_a,jtotal,ibin)
!   aer(iaro2_a,jsolid,ibin)= aer(iaro2_a,jtotal,ibin)
!   aer(ialk1_a,jsolid,ibin)= aer(ialk1_a,jtotal,ibin)
!   aer(iole1_a,jsolid,ibin)= aer(iole1_a,jtotal,ibin)
!   aer(iapi1_a,jsolid,ibin)= aer(iapi1_a,jtotal,ibin)
!   aer(iapi2_a,jsolid,ibin)= aer(iapi2_a,jtotal,ibin)
!   aer(ilim1_a,jsolid,ibin)= aer(ilim1_a,jtotal,ibin)
!   aer(ilim2_a,jsolid,ibin)= aer(ilim2_a,jtotal,ibin)

    ! liquid-phase
    do iaer = 1, naer
    aer(iaer,jliquid,ibin) = 0.0
    end do
    aer(iso4_a,jliquid,ibin) = max(0.0_r8, aer(iso4_a,jtotal,ibin) -   &
                               electrolyte(jcaso4,jsolid,ibin))      ! added max() RAZ 4/16/2014
    aer(ino3_a,jliquid,ibin) = aer(ino3_a,jtotal,ibin)
    aer(icl_a, jliquid,ibin) = aer(icl_a,jtotal,ibin)
    aer(inh4_a,jliquid,ibin) = aer(inh4_a,jtotal,ibin)
    aer(imsa_a,jliquid,ibin) = aer(imsa_a,jtotal,ibin)
    aer(ina_a, jliquid,ibin) = aer(ina_a,jtotal,ibin)
    aer(ica_a, jliquid,ibin) = electrolyte(jcano3,jtotal,ibin) +   &
                               electrolyte(jcacl2,jtotal,ibin)
!   aer(ioc_a, jliquid,ibin) = 0.0
!   aer(ico3_a,jliquid,ibin) = 0.0
!   aer(ibc_a, jliquid,ibin) = 0.0
!   aer(ioin_a,jliquid,ibin) = 0.0
!   aer(iaro1_a,jliquid,ibin)= 0.0
!   aer(iaro2_a,jliquid,ibin)= 0.0
!   aer(ialk1_a,jliquid,ibin)= 0.0
!   aer(iole1_a,jliquid,ibin)= 0.0
!   aer(iapi1_a,jliquid,ibin)= 0.0
!   aer(iapi2_a,jliquid,ibin)= 0.0
!   aer(ilim1_a,jliquid,ibin)= 0.0
!   aer(ilim2_a,jliquid,ibin)= 0.0

    return
  end subroutine do_full_deliquescence

  !***********************************************************************
  ! MESA: Multicomponent Equilibrium Solver for Aerosol-phase
  ! computes equilibrum solid and liquid phases by integrating
  ! pseudo-transient dissolution and precipitation reactions
  !
  ! author: Rahul A. Zaveri
  ! update: jan 2005
  ! Reference: Zaveri R.A., R.C. Easter, and L.K. Peters, JGR, 2005b
  !-----------------------------------------------------------------------
  subroutine MESA_PTC(ibin, jaerosolstate, jphase, aer, jhyst_leg,  &
       electrolyte, epercent, activity, mc, num_a, mass_dry_a, mass_wet_a, mass_soluble_a,    &
       vol_dry_a, vol_wet_a, water_a, aH2O, ma, gam,   &
       log_gamZ, zc, za, gam_ratio, xeq_a, na_Ma, nc_Mc, xeq_c, mw_electrolyte, mw_aer_mac,     &
       dens_aer_mac, Keq_sl, MW_c, MW_a, Keq_ll, growth_factor, molality0, rtol_mesa,         &
       jsalt_present, phi_salt_old,                                                &
       kappa_nonelectro, mosaic_vars_aa )                ! TOUCH

    use module_data_mosaic_aero,  only: r8, nbin_a_max, nelectrolyte, Ncation, naer, nsalt,  &!Parameters
         jhyst_lo, mixed, all_liquid, jsolid, jliquid, jtotal, mYES,                         &!Parameters
         all_solid, Nanion, nrxn_aer_sl, nrxn_aer_ll,                                     &!Parameters
         ino3_a, iso4_a, ioc_a, ilim1_a, ilim2_a, inh4_a, ina_a, ica_a, ico3_a, imsa_a, icl_a, &
         mosaic_vars_aa_type

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    integer, intent(inout), dimension(nsalt) :: jsalt_present
    integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate,jphase,jhyst_leg

    real(r8), intent(in) :: aH2O,rtol_mesa
    real(r8), intent(in), dimension(naer) :: mw_aer_mac,dens_aer_mac
    real(r8), intent(in), dimension(Ncation) :: zc,MW_c
    real(r8), intent(inout), dimension(Ncation) :: nc_Mc,xeq_c
    real(r8), intent(in), dimension(Nanion)  :: za,MW_a
    real(r8), intent(inout), dimension(Nanion)  :: xeq_a,na_Ma
    real(r8), intent(inout), dimension(nbin_a_max) :: num_a,mass_dry_a,mass_wet_a
    real(r8), intent(inout), dimension(nbin_a_max) :: mass_soluble_a,vol_dry_a
    real(r8), intent(inout), dimension(nbin_a_max) :: growth_factor
    real(r8), intent(inout), dimension(nbin_a_max) :: vol_wet_a,water_a,gam_ratio
    real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
    real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
    real(r8), intent(inout), dimension(nrxn_aer_sl) :: Keq_sl
    real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
    real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
    real(r8), intent(inout), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
    real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
    real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: epercent
    real(r8), intent(inout), dimension(nsalt) :: phi_salt_old
    real(r8), intent(in), dimension(naer) :: kappa_nonelectro

    type (mosaic_vars_aa_type), intent(inout) :: mosaic_vars_aa

    ! local variables
    integer iaer, iconverge, iconverge_flux, iconverge_mass,   &
         idissolved, itdum, js, je, jp

    real(r8) :: tau_p(nsalt), tau_d(nsalt)
    real(r8) :: frac_solid, sumflux, hsalt_min, alpha, XT, dumdum,   &
         H_ion
    real(r8) :: phi_prod, alpha_fac, sum_dum
    real(r8) :: aer_H,hsalt_max
    real(r8), dimension(nelectrolyte) :: eleliquid
    real(r8), dimension(nbin_a_max) :: mass_dry_salt
    real(r8), dimension(nsalt) :: phi_salt,flux_sl,phi_bar,alpha_salt
    real(r8), dimension(nsalt) :: sat_ratio,hsalt

    ! function
    !real(r8) :: aerosol_water

    ! initialize
    itdum = 0               ! initialize time
    hsalt_max = 1.e25

    do js = 1, nsalt
       hsalt(js)     = 0.0
       sat_ratio(js) = 0.0
       phi_salt(js)  = 0.0
       flux_sl(js)   = 0.0
    enddo

    !! EFFI calculate percent composition
    sum_dum = 0.0
    do je = 1, nelectrolyte
       sum_dum = sum_dum + electrolyte(je,jtotal,ibin)
    enddo

    if(sum_dum .eq. 0.)sum_dum = 1.0

    do je = 1, nelectrolyte
       epercent(je,jtotal,ibin) = 100.*electrolyte(je,jtotal,ibin)/sum_dum
    enddo
    !! EFFI

    do js = 1, nsalt
       jsalt_present(js) = 0                        ! default value - salt absent
       if(epercent(js,jtotal,ibin) .gt. 1.0)then
          jsalt_present(js) = 1                     ! salt present
       endif
    enddo

    mass_dry_a(ibin) = 0.0

    aer_H = (2.*aer(iso4_a,jtotal,ibin) +   &
         aer(ino3_a,jtotal,ibin) +   &
         aer(icl_a,jtotal,ibin)  +   &
         aer(imsa_a,jtotal,ibin) +   &
         2.*aer(ico3_a,jtotal,ibin))-   &
         (2.*aer(ica_a,jtotal,ibin)  +   &
         aer(ina_a,jtotal,ibin)  +   &
         aer(inh4_a,jtotal,ibin))
    aer_H = max(aer_H, 0.0d0)

    do iaer = 1, naer
       mass_dry_a(ibin) = mass_dry_a(ibin) +   &
            aer(iaer,jtotal,ibin)*mw_aer_mac(iaer)  ! [ng/m^3(air)]
       vol_dry_a(ibin)  = vol_dry_a(ibin) +   &
            aer(iaer,jtotal,ibin)*mw_aer_mac(iaer)/dens_aer_mac(iaer)       ! ncc/m^3(air)
    enddo
    mass_dry_a(ibin) = mass_dry_a(ibin) + aer_H
    vol_dry_a(ibin) = vol_dry_a(ibin) + aer_H

    mass_dry_a(ibin) = mass_dry_a(ibin)*1.e-15                      ! [g/cc(air)]
    vol_dry_a(ibin) = vol_dry_a(ibin)*1.e-15                                ! [cc(aer)/cc(air)]

    mass_dry_salt(ibin) = 0.0               ! soluble salts only
    do je = 1, nsalt
       mass_dry_salt(ibin) = mass_dry_salt(ibin) +   &
            electrolyte(je,jtotal,ibin)*mw_electrolyte(je)*1.e-15   ! g/cc(air)
    enddo

    mosaic_vars_aa%jMESA_call = mosaic_vars_aa%jMESA_call + 1

    !----begin pseudo time continuation loop-------------------------------

    do 500 itdum = 1, mosaic_vars_aa%Nmax_MESA

       ! compute new salt fluxes
       call MESA_flux_salt(ibin,jaerosolstate,jphase, aer,jhyst_leg,electrolyte, &
            epercent,activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a,aH2O,ma,&
            gam,log_gamZ,zc,za,gam_ratio,xeq_a,na_Ma,nc_Mc,xeq_c,mw_electrolyte, &
            Keq_sl,MW_c,MW_a,Keq_ll,eleliquid,flux_sl,phi_salt,sat_ratio,        &
            molality0,jsalt_present,kappa_nonelectro)

       ! check convergence
       call MESA_convergence_criterion(ibin,iconverge_mass,iconverge_flux,idissolved, &
            aer,electrolyte,mass_dry_a,mass_dry_salt,mw_electrolyte,mw_aer_mac, &
            flux_sl,phi_salt,rtol_mesa)

       if(iconverge_mass .eq. mYES)then
          mosaic_vars_aa%iter_MESA(ibin) = mosaic_vars_aa%iter_MESA(ibin) + itdum
          mosaic_vars_aa%niter_MESA = mosaic_vars_aa%niter_MESA + float(itdum)
          mosaic_vars_aa%niter_MESA_max = max( mosaic_vars_aa%niter_MESA_max, itdum)
          jaerosolstate(ibin) = all_solid
          call adjust_solid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent,   &
               water_a)
          jhyst_leg(ibin) = jhyst_lo
          growth_factor(ibin) = 1.0
          return
       elseif(iconverge_flux .eq. mYES)then
          mosaic_vars_aa%iter_MESA(ibin) = mosaic_vars_aa%iter_MESA(ibin) + itdum
          mosaic_vars_aa%niter_MESA = mosaic_vars_aa%niter_MESA + float(itdum)
          mosaic_vars_aa%niter_MESA_max = max( mosaic_vars_aa%niter_MESA_max, itdum)
          jaerosolstate(ibin) = mixed
          vol_wet_a(ibin)  = vol_dry_a(ibin) + water_a(ibin)*1.e-3          ! cc(aer)/cc(air) or m^3/m^3(air)
          growth_factor(ibin) = mass_wet_a(ibin)/mass_dry_a(ibin)           ! mass growth factor

          if(idissolved .eq. myes)then
             jaerosolstate(ibin) = all_liquid
             !          jhyst_leg(ibin) = jhyst_up  ! ! do this later (to avoid tripping kelvin iterations)
          else
             jaerosolstate(ibin) = mixed
             jhyst_leg(ibin) = jhyst_lo
          endif

          ! calculate epercent(jsolid) composition in mixed-phase aerosol EFFI
          !!        sum_dum = 0.0
          !!        jp = jsolid
          !!        do je = 1, nelectrolyte
          !!          electrolyte(je,jp,ibin) = max(0.d0,electrolyte(je,jp,ibin)) ! remove -ve
          !!          sum_dum = sum_dum + electrolyte(je,jp,ibin)
          !!        enddo
          !!        electrolyte_sum(jp,ibin) = sum_dum
          !!        if(sum_dum .eq. 0.)sum_dum = 1.0
          !!        do je = 1, nelectrolyte
          !!          epercent(je,jp,ibin) = 100.*electrolyte(je,jp,ibin)/sum_dum
          !!        enddo

          return
       endif

       ! calculate hsalt(js)        ! time step
       hsalt_min = 1.e25

       do js = 1, nsalt

          phi_prod = phi_salt(js) * phi_salt_old(js)

          if(itdum .gt. 1 .and. phi_prod .gt. 0.0)then
             phi_bar(js) = (abs(phi_salt(js))-abs(phi_salt_old(js)))/   &
                  alpha_salt(js)
          else
             phi_bar(js) = 0.0                      ! oscillating, or phi_salt and/or phi_salt_old may be zero
          endif

          if(phi_bar(js) .lt. 0.0)then              ! good. phi getting lower. maybe able to take bigger alphas
             phi_bar(js) = max(phi_bar(js), -10.0d0)
             alpha_fac = 3.0*exp(phi_bar(js))
             alpha_salt(js) = min(alpha_fac*abs(phi_salt(js)), 0.9d0)
          elseif(phi_bar(js) .gt. 0.0)then  ! bad - phi is getting bigger. so be conservative with alpha
             alpha_salt(js) = min(abs(phi_salt(js)), 0.5d0)
          else                                      ! very bad - phi is oscillating. be very conservative
             alpha_salt(js) = min(abs(phi_salt(js))/3.0d0, 0.5d0)
          endif

          !        alpha_salt(js) = max(alpha_salt(js), 0.01)

          phi_salt_old(js) = phi_salt(js)           ! update old array

          if(flux_sl(js) .gt. 0.)then

             tau_p(js) = eleliquid(js)/flux_sl(js)  ! precipitation time scale
             if(tau_p(js) .eq. 0.0)then
                hsalt(js) = 1.e25
                flux_sl(js) = 0.0
                phi_salt(js)= 0.0
             else
                hsalt(js) = alpha_salt(js)*tau_p(js)
             endif

          elseif(flux_sl(js) .lt. 0.)then

             tau_p(js) = -eleliquid(js)/flux_sl(js) ! precipitation time scale
             tau_d(js) = -electrolyte(js,jsolid,ibin)/flux_sl(js) ! dissolution time scale
             if(tau_p(js) .eq. 0.0)then
                hsalt(js) = alpha_salt(js)*tau_d(js)
             else
                hsalt(js) = alpha_salt(js)*min(tau_p(js),tau_d(js))
             endif

          else

             hsalt(js) = 1.e25

          endif

          hsalt_min = min(hsalt(js), hsalt_min)

       enddo

       !---------------------------------

       ! integrate electrolyte(solid)
       do js = 1, nsalt
          electrolyte(js,jsolid,ibin) = (   &
               (electrolyte(js,jsolid,ibin))  +   &
               (hsalt(js)) * (flux_sl(js)) )
       enddo

       ! compute aer(solid) from electrolyte(solid)
       call electrolytes_to_ions(jsolid,ibin,aer,electrolyte)

       ! compute new electrolyte(liquid) from mass balance
       do iaer = 1, naer
          aer(iaer,jliquid,ibin) = ( (aer(iaer,jtotal,ibin)) -   &
               (aer(iaer,jsolid,ibin)) )
       enddo

       !---------------------------------

500 continue     ! end time continuation loop
    !--------------------------------------------------------------------
    mosaic_vars_aa%jMESA_fail = mosaic_vars_aa%jMESA_fail + 1
    mosaic_vars_aa%iter_MESA(ibin) = mosaic_vars_aa%iter_MESA(ibin) + itdum
    mosaic_vars_aa%niter_MESA = mosaic_vars_aa%niter_MESA + float(itdum)
    jaerosolstate(ibin) = mixed
    jhyst_leg(ibin) = jhyst_lo
    mass_wet_a(ibin) = mass_dry_a(ibin) + water_a(ibin)*1.e-3    ! g/cc(air)
    vol_wet_a(ibin)  = vol_dry_a(ibin) + water_a(ibin)*1.e-3     ! cc(aer)/cc(air) or m^3/m^3(air)
    growth_factor(ibin) = mass_wet_a(ibin)/mass_dry_a(ibin)      ! mass growth factor

    return
  end subroutine MESA_PTC

  !***********************************************************************
  ! part of MESA: calculates solid-liquid fluxes of soluble salts
  !
  ! author: Rahul A. Zaveri
  ! update: jan 2005
  !-----------------------------------------------------------------------
  subroutine MESA_flux_salt(ibin, jaerosolstate,jphase,aer,jhyst_leg,electrolyte,  &
       epercent,activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a,aH2O,ma,gam,   &
       log_gamZ,zc,za,gam_ratio,xeq_a,na_Ma,nc_Mc,xeq_c,mw_electrolyte,Keq_sl,MW_c,&
       MW_a,Keq_ll,eleliquid,flux_sl,phi_salt,sat_ratio,molality0,jsalt_present,   &
       kappa_nonelectro                                                            )      ! TOUCH

    use module_data_mosaic_aero, only: r8,nbin_a_max,nelectrolyte,Ncation,naer,    &!Parameters
         jliquid,nsalt,jsolid,Nanion,nrxn_aer_sl,nrxn_aer_ll,nrxn_aer_sl,          &!Parameter
         jna3hso4,ica_a,jcano3,jcacl2                                               !TBD

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    integer, intent(inout), dimension(nsalt) :: jsalt_present
    integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate,jphase,jhyst_leg

    real(r8), intent(in) :: aH2O
    real(r8), intent(inout), dimension(nsalt) :: flux_sl,phi_salt,sat_ratio
    real(r8), intent(in), dimension(Ncation) :: zc,MW_c
    real(r8), intent(inout), dimension(Ncation) :: nc_Mc,xeq_c
    real(r8), intent(in), dimension(Nanion)  :: za,MW_a
    real(r8), intent(inout), dimension(Nanion)  :: xeq_a,na_Ma
    real(r8), intent(inout), dimension(nbin_a_max) :: num_a,mass_dry_a,gam_ratio
    real(r8), intent(inout), dimension(nbin_a_max) :: mass_soluble_a,water_a
    real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
    real(r8), intent(inout), dimension(nelectrolyte) :: eleliquid
    real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
    real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
    real(r8), intent(inout), dimension(nrxn_aer_sl) :: Keq_sl
    real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
    real(r8), intent(inout), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
    real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
    real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte,epercent
    real(r8), intent(in), dimension(naer) :: kappa_nonelectro

    ! local variables
    integer js, je
    real(r8) :: XT, calcium, sum_salt, sum_dum !**BALLI XT should it be subr arg??
    real(r8), dimension(nsalt) :: frac_salt_liq,frac_salt_solid

    ! compute activities and water content
    call ions_to_electrolytes(jliquid,ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma,   &
         nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a)
    call compute_activities(ibin,jaerosolstate,jphase,aer,jhyst_leg,electrolyte,   &
         activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a,aH2O,ma,gam,log_gamZ, &
         gam_ratio,Keq_ll,molality0,kappa_nonelectro)
    activity(jna3hso4,ibin)   = 0.0

    if(water_a(ibin) .le. 0.0)then
       do js = 1, nsalt
          flux_sl(js) = 0.0
       enddo
       return
    endif

    call MESA_estimate_eleliquid(ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma,nc_Mc,  &
         xeq_c,mw_electrolyte,MW_c,MW_a,eleliquid)

    calcium = aer(ica_a,jliquid,ibin)

    !! EFFI calculate percent composition
    sum_dum = 0.0
    do je = 1, nelectrolyte
       sum_dum = sum_dum + electrolyte(je,jliquid,ibin)
    enddo

    if(sum_dum .eq. 0.)sum_dum = 1.0

    do je = 1, nelectrolyte
       epercent(je,jliquid,ibin) = 100.*electrolyte(je,jliquid,ibin)/sum_dum
    enddo
    !! EFFI

    ! calculate % electrolyte composition in the solid and liquid phases
    sum_salt = 0.0
    do js = 1, nsalt
       sum_salt = sum_salt + electrolyte(js,jsolid,ibin)
    enddo

    if(sum_salt .eq. 0.0)sum_salt = 1.0
    do js = 1, nsalt
       frac_salt_solid(js) = electrolyte(js,jsolid,ibin)/sum_salt
       frac_salt_liq(js)   = epercent(js,jliquid,ibin)/100.
    enddo

    ! compute salt fluxes
    do js = 1, nsalt             ! soluble solid salts

       ! compute new saturation ratio
       sat_ratio(js) = activity(js,ibin)/Keq_sl(js)
       ! compute relative driving force
       phi_salt(js)  = (sat_ratio(js) - 1.0)/max(sat_ratio(js),1.0d0)

       ! check if too little solid-phase salt is trying to dissolve
       if(sat_ratio(js)       .lt. 1.00 .and.   &
            frac_salt_solid(js) .lt. 0.01 .and.   &
            frac_salt_solid(js) .gt. 0.0)then
          call MESA_dissolve_small_salt(ibin,js,aer,electrolyte)
          call MESA_estimate_eleliquid(ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma,  &
               nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a,eleliquid)
          sat_ratio(js) = activity(js,ibin)/Keq_sl(js)
       endif

       ! compute flux
       flux_sl(js) = sat_ratio(js) - 1.0

       ! apply Heaviside function
       if( (sat_ratio(js)               .lt. 1.0 .and.   &
            electrolyte(js,jsolid,ibin) .eq. 0.0) .or.   &
            (calcium .gt. 0.0 .and. frac_salt_liq(js).lt.0.01).or.   &
            (calcium .gt. 0.0 .and. jsalt_present(js).eq.0) )then
          flux_sl(js) = 0.0
          phi_salt(js)= 0.0
       endif

    enddo

    ! force cacl2 and cano3 fluxes to zero
    sat_ratio(jcano3) = 1.0
    phi_salt(jcano3)  = 0.0
    flux_sl(jcano3)   = 0.0

    sat_ratio(jcacl2) = 1.0
    phi_salt(jcacl2)  = 0.0
    flux_sl(jcacl2)   = 0.0

    return
  end subroutine MESA_flux_salt

 !***********************************************************************
  ! computes activities
  !
  ! author: Rahul A. Zaveri
  ! update: jan 2007
  !-----------------------------------------------------------------------
  subroutine compute_activities(ibin,jaerosolstate,jphase,aer,jhyst_leg,           &
       electrolyte,activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a,aH2O,ma,gam,&
       log_gamZ,gam_ratio,Keq_ll,molality0,kappa_nonelectro)

    use module_data_mosaic_aero, only: r8,nbin_a_max,nelectrolyte,Ncation,naer,    &
         jliquid,Nanion,nrxn_aer_ll,                                               &
         iso4_a,ja_so4,ja_hso4,ino3_a,ja_no3,icl_a,ja_cl,imsa_a,ja_msa,ica_a,jc_ca,&
         inh4_a,jc_nh4,ina_a,jc_na,jc_h,jhcl,jhno3,jcacl2,jcano3,jnacl,jnano3,     &
         jna2so4,jnh4so4,jnh4cl,jnh4no3,jlvcite,jnh4hso4,jnh4msa,jna3hso4,jnahso4, &
         jnamsa,jcamsa2,jh2so4,jhhso4,jmsa                                          !TBD

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate,jphase,jhyst_leg

    real(r8), intent(in) :: aH2O
    real(r8), intent(in), dimension(nbin_a_max) :: num_a
    real(r8), intent(inout), dimension(nbin_a_max) :: mass_dry_a,mass_soluble_a
    real(r8), intent(inout), dimension(nbin_a_max) :: water_a,gam_ratio
    real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
    real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
    real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
    real(r8), intent(inout), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
    real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
    real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(in), dimension(naer) :: kappa_nonelectro

    ! local variables
    real(r8), dimension(nelectrolyte) :: log_gam
    integer jp, jA
    real(r8) :: XT, xmol(Nelectrolyte), sum_elec, dumK, c_bal, a_c !BALLI** should xt be subr arg??
    real(r8) :: quad, aq, bq, cq, xq, dum, mSULF
    !real(r8) :: aerosol_water     ! mosaic function

    water_a(ibin) = aerosol_water(jliquid,ibin,jaerosolstate,jphase, &
         jhyst_leg,electrolyte,aer,kappa_nonelectro,num_a,mass_dry_a,mass_soluble_a,aH2O, &
         molality0)      ! Kg/m^3(air)
    if(water_a(ibin) .eq. 0.0)return

    call calculate_XT(ibin,jliquid,XT,aer)

    if(XT.ge.2.0 .or. XT.lt.0.)then   ! changed .gt. to .ge.   RAZ 4/16/2014
       ! SULFATE POOR: fully dissociated electrolytes

       ! anion molalities (mol/kg water)
       ma(ja_so4,ibin)  = 1.e-9*aer(iso4_a,jliquid,ibin)/water_a(ibin)
       ma(ja_hso4,ibin) = 0.0
       ma(ja_no3,ibin)  = 1.e-9*aer(ino3_a,jliquid,ibin)/water_a(ibin)
       ma(ja_cl,ibin)   = 1.e-9*aer(icl_a, jliquid,ibin)/water_a(ibin)
       ma(ja_msa,ibin)  = 1.e-9*aer(imsa_a,jliquid,ibin)/water_a(ibin)

       ! cation molalities (mol/kg water)
       mc(jc_ca,ibin)   = 1.e-9*aer(ica_a, jliquid,ibin)/water_a(ibin)
       mc(jc_nh4,ibin)  = 1.e-9*aer(inh4_a,jliquid,ibin)/water_a(ibin)
       mc(jc_na,ibin)   = 1.e-9*aer(ina_a, jliquid,ibin)/water_a(ibin)
       a_c              = (   &
            (2.*ma(ja_so4,ibin)+   &
            ma(ja_no3,ibin)+   &
            ma(ja_cl,ibin) +   &
            ma(ja_msa,ibin)) -   &
            (2.*mc(jc_ca,ibin) +   &
            mc(jc_nh4,ibin)+   &
            mc(jc_na,ibin)) )

       mc(jc_h,ibin) = 0.5*( (a_c) +   &
            (sqrt(a_c**2 + 4.*Keq_ll(3))) )

       if(mc(jc_h,ibin) .le. 0.0)then   ! changed .eq. to .le. RAZ 4/16/2014
          mc(jc_h,ibin) = 1.e-10
       endif

       jp = jliquid

       sum_elec = 2.*electrolyte(jnh4no3,jp,ibin) +   &
            2.*electrolyte(jnh4cl,jp,ibin)  +   &
            3.*electrolyte(jnh4so4,jp,ibin) +   &
            3.*electrolyte(jna2so4,jp,ibin) +   &
            2.*electrolyte(jnano3,jp,ibin)  +   &
            2.*electrolyte(jnacl,jp,ibin)   +   &
            3.*electrolyte(jcano3,jp,ibin)  +   &
            3.*electrolyte(jcacl2,jp,ibin)  +   &
            2.*electrolyte(jhno3,jp,ibin)   +   &
            2.*electrolyte(jhcl,jp,ibin)

       if(sum_elec .eq. 0.0)then
          do jA = 1, nelectrolyte
             gam(jA,ibin) = 1.0
          enddo
          goto 10
       endif

       ! ionic mole fractions
       xmol(jnh4no3) = 2.*electrolyte(jnh4no3,jp,ibin)/sum_elec
       xmol(jnh4cl)  = 2.*electrolyte(jnh4cl,jp,ibin) /sum_elec
       xmol(jnh4so4) = 3.*electrolyte(jnh4so4,jp,ibin)/sum_elec
       xmol(jna2so4) = 3.*electrolyte(jna2so4,jp,ibin)/sum_elec
       xmol(jnano3)  = 2.*electrolyte(jnano3,jp,ibin) /sum_elec
       xmol(jnacl)   = 2.*electrolyte(jnacl,jp,ibin)  /sum_elec
       xmol(jcano3)  = 3.*electrolyte(jcano3,jp,ibin) /sum_elec
       xmol(jcacl2)  = 3.*electrolyte(jcacl2,jp,ibin) /sum_elec
       xmol(jhno3)   = 2.*electrolyte(jhno3,jp,ibin)  /sum_elec
       xmol(jhcl)    = 2.*electrolyte(jhcl,jp,ibin)   /sum_elec

       jA = jnh4so4
       if(xmol(jA).gt.0.0)then
          log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) +   &
               xmol(jnh4cl) *log_gamZ(jA,jnh4cl)  +   &
               xmol(jnh4so4)*log_gamZ(jA,jnh4so4) +   &
               xmol(jna2so4)*log_gamZ(jA,jna2so4) +   &
               xmol(jnano3) *log_gamZ(jA,jnano3)  +   &
               xmol(jnacl)  *log_gamZ(jA,jnacl)   +   &
               xmol(jcano3) *log_gamZ(jA,jcano3)  +   &
               xmol(jcacl2) *log_gamZ(jA,jcacl2)  +   &
               xmol(jhno3)  *log_gamZ(jA,jhno3)   +   &
               xmol(jhcl)   *log_gamZ(jA,jhcl)
          gam(jA,ibin) = 10.**log_gam(jA)
          activity(jnh4so4,ibin) = mc(jc_nh4,ibin)**2 * ma(ja_so4,ibin) *   &
               gam(jnh4so4,ibin)**3
       endif

! RAZ 11/7/2014
! always calculate gam(jnh4no3), even if xmol(jnh4no3) = 0. this to calculate gam_ratio
       jA = jnh4no3
!      if(xmol(jA).gt.0.0)then
          log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) +   &
               xmol(jnh4cl) *log_gamZ(jA,jnh4cl)  +   &
               xmol(jnh4so4)*log_gamZ(jA,jnh4so4) +   &
               xmol(jna2so4)*log_gamZ(jA,jna2so4) +   &
               xmol(jnano3) *log_gamZ(jA,jnano3)  +   &
               xmol(jnacl)  *log_gamZ(jA,jnacl)   +   &
               xmol(jcano3) *log_gamZ(jA,jcano3)  +   &
               xmol(jcacl2) *log_gamZ(jA,jcacl2)  +   &
               xmol(jhno3)  *log_gamZ(jA,jhno3)   +   &
               xmol(jhcl)   *log_gamZ(jA,jhcl)
          gam(jA,ibin) = 10.**log_gam(jA)
          activity(jnh4no3,ibin) = mc(jc_nh4,ibin) * ma(ja_no3,ibin) *   &
               gam(jnh4no3,ibin)**2
!      endif

       jA = jnh4cl
       if(xmol(jA).gt.0.0)then
          log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) +   &
               xmol(jnh4cl) *log_gamZ(jA,jnh4cl)  +   &
               xmol(jnh4so4)*log_gamZ(jA,jnh4so4) +   &
               xmol(jna2so4)*log_gamZ(jA,jna2so4) +   &
               xmol(jnano3) *log_gamZ(jA,jnano3)  +   &
               xmol(jnacl)  *log_gamZ(jA,jnacl)   +   &
               xmol(jcano3) *log_gamZ(jA,jcano3)  +   &
               xmol(jcacl2) *log_gamZ(jA,jcacl2)  +   &
               xmol(jhno3)  *log_gamZ(jA,jhno3)   +   &
               xmol(jhcl)   *log_gamZ(jA,jhcl)
          gam(jA,ibin) = 10.**log_gam(jA)
          activity(jnh4cl,ibin)  = mc(jc_nh4,ibin) * ma(ja_cl,ibin) *   &
               gam(jnh4cl,ibin)**2
       endif

       jA = jna2so4
       if(xmol(jA).gt.0.0)then
          log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) +   &
               xmol(jnh4cl) *log_gamZ(jA,jnh4cl)  +   &
               xmol(jnh4so4)*log_gamZ(jA,jnh4so4) +   &
               xmol(jna2so4)*log_gamZ(jA,jna2so4) +   &
               xmol(jnano3) *log_gamZ(jA,jnano3)  +   &
               xmol(jnacl)  *log_gamZ(jA,jnacl)   +   &
               xmol(jcano3) *log_gamZ(jA,jcano3)  +   &
               xmol(jcacl2) *log_gamZ(jA,jcacl2)  +   &
               xmol(jhno3)  *log_gamZ(jA,jhno3)   +   &
               xmol(jhcl)   *log_gamZ(jA,jhcl)
          gam(jA,ibin) = 10.**log_gam(jA)
          activity(jna2so4,ibin) = mc(jc_na,ibin)**2 * ma(ja_so4,ibin) *   &
               gam(jna2so4,ibin)**3
       endif

       jA = jnano3
       if(xmol(jA).gt.0.0)then
          log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) +   &
               xmol(jnh4cl) *log_gamZ(jA,jnh4cl)  +   &
               xmol(jnh4so4)*log_gamZ(jA,jnh4so4) +   &
               xmol(jna2so4)*log_gamZ(jA,jna2so4) +   &
               xmol(jnano3) *log_gamZ(jA,jnano3)  +   &
               xmol(jnacl)  *log_gamZ(jA,jnacl)   +   &
               xmol(jcano3) *log_gamZ(jA,jcano3)  +   &
               xmol(jcacl2) *log_gamZ(jA,jcacl2)  +   &
               xmol(jhno3)  *log_gamZ(jA,jhno3)   +   &
               xmol(jhcl)   *log_gamZ(jA,jhcl)
          gam(jA,ibin) = 10.**log_gam(jA)
          activity(jnano3,ibin)  = mc(jc_na,ibin) * ma(ja_no3,ibin) *   &
               gam(jnano3,ibin)**2
       endif

       jA = jnacl
       if(xmol(jA).gt.0.0)then
          log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) +   &
               xmol(jnh4cl) *log_gamZ(jA,jnh4cl)  +   &
               xmol(jnh4so4)*log_gamZ(jA,jnh4so4) +   &
               xmol(jna2so4)*log_gamZ(jA,jna2so4) +   &
               xmol(jnano3) *log_gamZ(jA,jnano3)  +   &
               xmol(jnacl)  *log_gamZ(jA,jnacl)   +   &
               xmol(jcano3) *log_gamZ(jA,jcano3)  +   &
               xmol(jcacl2) *log_gamZ(jA,jcacl2)  +   &
               xmol(jhno3)  *log_gamZ(jA,jhno3)   +   &
               xmol(jhcl)   *log_gamZ(jA,jhcl)
          gam(jA,ibin) = 10.**log_gam(jA)
          activity(jnacl,ibin)   = mc(jc_na,ibin) * ma(ja_cl,ibin) *   &
               gam(jnacl,ibin)**2
       endif

       !c      jA = jcano3
       !c      if(xmol(jA).gt.0.0)then
       !c      gam(jA,ibin) = 1.0
       !c      activity(jcano3,ibin)  = 1.0
       !c      endif

       !c      jA = jcacl2
       !c      if(xmol(jA).gt.0.0)then
       !c      gam(jA,ibin) = 1.0
       !c      activity(jcacl2,ibin)  = 1.0
       !c      endif

       jA = jcano3
       if(xmol(jA).gt.0.0)then
          log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) +   &
               xmol(jnh4cl) *log_gamZ(jA,jnh4cl)  +   &
               xmol(jnh4so4)*log_gamZ(jA,jnh4so4) +   &
               xmol(jna2so4)*log_gamZ(jA,jna2so4) +   &
               xmol(jnano3) *log_gamZ(jA,jnano3)  +   &
               xmol(jnacl)  *log_gamZ(jA,jnacl)   +   &
               xmol(jcano3) *log_gamZ(jA,jcano3)  +   &
               xmol(jcacl2) *log_gamZ(jA,jcacl2)  +   &
               xmol(jhno3)  *log_gamZ(jA,jhno3)   +   &
               xmol(jhcl)   *log_gamZ(jA,jhcl)
          gam(jA,ibin) = 10.**log_gam(jA)
          activity(jcano3,ibin)  = mc(jc_ca,ibin) * ma(ja_no3,ibin)**2 *   &
               gam(jcano3,ibin)**3
       endif

       jA = jcacl2
       if(xmol(jA).gt.0.0)then
          log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) +   &
               xmol(jnh4cl) *log_gamZ(jA,jnh4cl)  +   &
               xmol(jnh4so4)*log_gamZ(jA,jnh4so4) +   &
               xmol(jna2so4)*log_gamZ(jA,jna2so4) +   &
               xmol(jnano3) *log_gamZ(jA,jnano3)  +   &
               xmol(jnacl)  *log_gamZ(jA,jnacl)   +   &
               xmol(jcano3) *log_gamZ(jA,jcano3)  +   &
               xmol(jcacl2) *log_gamZ(jA,jcacl2)  +   &
               xmol(jhno3)  *log_gamZ(jA,jhno3)   +   &
               xmol(jhcl)   *log_gamZ(jA,jhcl)
          gam(jA,ibin) = 10.**log_gam(jA)
          activity(jcacl2,ibin)  = mc(jc_ca,ibin) * ma(ja_cl,ibin)**2 *   &
               gam(jcacl2,ibin)**3
       endif

       jA = jhno3
       log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) +   &
            xmol(jnh4cl) *log_gamZ(jA,jnh4cl)  +   &
            xmol(jnh4so4)*log_gamZ(jA,jnh4so4) +   &
            xmol(jna2so4)*log_gamZ(jA,jna2so4) +   &
            xmol(jnano3) *log_gamZ(jA,jnano3)  +   &
            xmol(jnacl)  *log_gamZ(jA,jnacl)   +   &
            xmol(jcano3) *log_gamZ(jA,jcano3)  +   &
            xmol(jcacl2) *log_gamZ(jA,jcacl2)  +   &
            xmol(jhno3)  *log_gamZ(jA,jhno3)   +   &
            xmol(jhcl)   *log_gamZ(jA,jhcl)
       gam(jA,ibin) = 10.**log_gam(jA)
       activity(jhno3,ibin)   = mc(jc_h,ibin) * ma(ja_no3,ibin) *   &
            gam(jhno3,ibin)**2

       jA = jhcl
       log_gam(jA) = xmol(jnh4no3)*log_gamZ(jA,jnh4no3) +   &
            xmol(jnh4cl) *log_gamZ(jA,jnh4cl)  +   &
            xmol(jnh4so4)*log_gamZ(jA,jnh4so4) +   &
            xmol(jna2so4)*log_gamZ(jA,jna2so4) +   &
            xmol(jnano3) *log_gamZ(jA,jnano3)  +   &
            xmol(jnacl)  *log_gamZ(jA,jnacl)   +   &
            xmol(jcano3) *log_gamZ(jA,jcano3)  +   &
            xmol(jcacl2) *log_gamZ(jA,jcacl2)  +   &
            xmol(jhno3)  *log_gamZ(jA,jhno3)   +   &
            xmol(jhcl)   *log_gamZ(jA,jhcl)
       gam(jA,ibin) = 10.**log_gam(jA)
       activity(jhcl,ibin)    = mc(jc_h,ibin) * ma(ja_cl,ibin) *   &
            gam(jhcl,ibin)**2

       !----
10     gam(jlvcite,ibin) = 1.0

       gam(jnh4hso4,ibin)= 1.0

       gam(jnh4msa,ibin) = 1.0

       gam(jna3hso4,ibin) = 1.0

       gam(jnahso4,ibin) = 1.0

       gam(jnamsa,ibin)  = 1.0

       gam(jcamsa2,ibin) = 1.0

       activity(jlvcite,ibin) = 0.0

       activity(jnh4hso4,ibin)= 0.0

       activity(jnh4msa,ibin) = mc(jc_nh4,ibin) * ma(ja_msa,ibin) *   &
            gam(jnh4msa,ibin)**2

       activity(jna3hso4,ibin)= 0.0

       activity(jnahso4,ibin) = 0.0

       activity(jnamsa,ibin) = mc(jc_na,ibin) * ma(ja_msa,ibin) *   &
            gam(jnamsa,ibin)**2

       activity(jcamsa2,ibin) = mc(jc_ca,ibin) * ma(ja_msa,ibin)**2 *   &
            gam(jcamsa2,ibin)**3

       gam_ratio(ibin) = gam(jnh4no3,ibin)**2/gam(jhno3,ibin)**2

    else
       !  SULFATE-RICH: solve for SO4= and HSO4- ions

       jp = jliquid

       sum_elec = 3.*electrolyte(jh2so4,jp,ibin)    +   &
            2.*electrolyte(jnh4hso4,jp,ibin)  +   &
            5.*electrolyte(jlvcite,jp,ibin)   +   &
            3.*electrolyte(jnh4so4,jp,ibin)   +   &
            2.*electrolyte(jnahso4,jp,ibin)   +   &
            5.*electrolyte(jna3hso4,jp,ibin)  +   &
            3.*electrolyte(jna2so4,jp,ibin)   +   &
            2.*electrolyte(jhno3,jp,ibin)     +   &
            2.*electrolyte(jhcl,jp,ibin)

       if(sum_elec .eq. 0.0)then
          do jA = 1, nelectrolyte
             gam(jA,ibin) = 1.0
          enddo
          goto 20
       endif

       xmol(jh2so4)  = 3.*electrolyte(jh2so4,jp,ibin)/sum_elec
       xmol(jnh4hso4)= 2.*electrolyte(jnh4hso4,jp,ibin)/sum_elec
       xmol(jlvcite) = 5.*electrolyte(jlvcite,jp,ibin)/sum_elec
       xmol(jnh4so4) = 3.*electrolyte(jnh4so4,jp,ibin)/sum_elec
       xmol(jnahso4) = 2.*electrolyte(jnahso4,jp,ibin)/sum_elec
       xmol(jna3hso4)= 5.*electrolyte(jna3hso4,jp,ibin)/sum_elec
       xmol(jna2so4) = 3.*electrolyte(jna2so4,jp,ibin)/sum_elec
       xmol(jhno3)   = 2.*electrolyte(jhno3,jp,ibin)/sum_elec
       xmol(jhcl)    = 2.*electrolyte(jhcl,jp,ibin)/sum_elec

       ! 2H.SO4
       jA = jh2so4
       log_gam(jA) = xmol(jh2so4)  *log_gamZ(jA,jh2so4)  +   &
            xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+   &
            xmol(jlvcite) *log_gamZ(jA,jlvcite) +   &
            xmol(jnh4so4) *log_gamZ(jA,jnh4so4) +   &
            xmol(jnahso4) *log_gamZ(jA,jnahso4) +   &
            xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+   &
            xmol(jna2so4) *log_gamZ(jA,jna2so4) +   &
            xmol(jhno3)   *log_gamZ(jA,jhno3)   +   &
            xmol(jhcl)    *log_gamZ(jA,jhcl)
       gam(jA,ibin) = 10.**log_gam(jA)

       ! H.HSO4
       jA = jhhso4
       log_gam(jA) = xmol(jh2so4)  *log_gamZ(jA,jh2so4)  +   &
            xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+   &
            xmol(jlvcite) *log_gamZ(jA,jlvcite) +   &
            xmol(jnh4so4) *log_gamZ(jA,jnh4so4) +   &
            xmol(jnahso4) *log_gamZ(jA,jnahso4) +   &
            xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+   &
            xmol(jna2so4) *log_gamZ(jA,jna2so4) +   &
            xmol(jhno3)   *log_gamZ(jA,jhno3)   +   &
            xmol(jhcl)    *log_gamZ(jA,jhcl)
       gam(jA,ibin) = 10.**log_gam(jA)

       ! NH4HSO4
       jA = jnh4hso4
       log_gam(jA) = xmol(jh2so4)  *log_gamZ(jA,jh2so4)  +   &
            xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+   &
            xmol(jlvcite) *log_gamZ(jA,jlvcite) +   &
            xmol(jnh4so4) *log_gamZ(jA,jnh4so4) +   &
            xmol(jnahso4) *log_gamZ(jA,jnahso4) +   &
            xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+   &
            xmol(jna2so4) *log_gamZ(jA,jna2so4) +   &
            xmol(jhno3)   *log_gamZ(jA,jhno3)   +   &
            xmol(jhcl)    *log_gamZ(jA,jhcl)
       gam(jA,ibin) = 10.**log_gam(jA)

       ! LETOVICITE
       jA = jlvcite
       log_gam(jA) = xmol(jh2so4)  *log_gamZ(jA,jh2so4)  +   &
            xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+   &
            xmol(jlvcite) *log_gamZ(jA,jlvcite) +   &
            xmol(jnh4so4) *log_gamZ(jA,jnh4so4) +   &
            xmol(jnahso4) *log_gamZ(jA,jnahso4) +   &
            xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+   &
            xmol(jna2so4) *log_gamZ(jA,jna2so4) +   &
            xmol(jhno3)   *log_gamZ(jA,jhno3)   +   &
            xmol(jhcl)    *log_gamZ(jA,jhcl)
       gam(jA,ibin) = 10.**log_gam(jA)

       ! (NH4)2SO4
       jA = jnh4so4
       log_gam(jA) = xmol(jh2so4)  *log_gamZ(jA,jh2so4)  +   &
            xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+   &
            xmol(jlvcite) *log_gamZ(jA,jlvcite) +   &
            xmol(jnh4so4) *log_gamZ(jA,jnh4so4) +   &
            xmol(jnahso4) *log_gamZ(jA,jnahso4) +   &
            xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+   &
            xmol(jna2so4) *log_gamZ(jA,jna2so4) +   &
            xmol(jhno3)   *log_gamZ(jA,jhno3)   +   &
            xmol(jhcl)    *log_gamZ(jA,jhcl)
       gam(jA,ibin) = 10.**log_gam(jA)

       ! NaHSO4
       jA = jnahso4
       log_gam(jA) = xmol(jh2so4)  *log_gamZ(jA,jh2so4)  +   &
            xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+   &
            xmol(jlvcite) *log_gamZ(jA,jlvcite) +   &
            xmol(jnh4so4) *log_gamZ(jA,jnh4so4) +   &
            xmol(jnahso4) *log_gamZ(jA,jnahso4) +   &
            xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+   &
            xmol(jna2so4) *log_gamZ(jA,jna2so4) +   &
            xmol(jhno3)   *log_gamZ(jA,jhno3)   +   &
            xmol(jhcl)    *log_gamZ(jA,jhcl)
       gam(jA,ibin) = 10.**log_gam(jA)

       ! Na3H(SO4)2
       jA = jna3hso4
       !      log_gam(jA) = xmol(jh2so4)  *log_gamZ(jA,jh2so4)  +
       !     &              xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+
       !     &              xmol(jlvcite) *log_gamZ(jA,jlvcite) +
       !     &              xmol(jnh4so4) *log_gamZ(jA,jnh4so4) +
       !     &              xmol(jnahso4) *log_gamZ(jA,jnahso4) +
       !     &              xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+
       !     &              xmol(jna2so4) *log_gamZ(jA,jna2so4) +
       !     &              xmol(jhno3)   *log_gamZ(jA,jhno3)   +
       !     &              xmol(jhcl)    *log_gamZ(jA,jhcl)
       !      gam(jA,ibin) = 10.**log_gam(jA)
       gam(jA,ibin) = 1.0

       ! Na2SO4
       jA = jna2so4
       log_gam(jA) = xmol(jh2so4)  *log_gamZ(jA,jh2so4)  +   &
            xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+   &
            xmol(jlvcite) *log_gamZ(jA,jlvcite) +   &
            xmol(jnh4so4) *log_gamZ(jA,jnh4so4) +   &
            xmol(jnahso4) *log_gamZ(jA,jnahso4) +   &
            xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+   &
            xmol(jna2so4) *log_gamZ(jA,jna2so4) +   &
            xmol(jhno3)   *log_gamZ(jA,jhno3)   +   &
            xmol(jhcl)    *log_gamZ(jA,jhcl)
       gam(jA,ibin) = 10.**log_gam(jA)

       ! HNO3
       jA = jhno3
       log_gam(jA) = xmol(jh2so4)  *log_gamZ(jA,jh2so4)  +   &
            xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+   &
            xmol(jlvcite) *log_gamZ(jA,jlvcite) +   &
            xmol(jnh4so4) *log_gamZ(jA,jnh4so4) +   &
            xmol(jnahso4) *log_gamZ(jA,jnahso4) +   &
            xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+   &
            xmol(jna2so4) *log_gamZ(jA,jna2so4) +   &
            xmol(jhno3)   *log_gamZ(jA,jhno3)   +   &
            xmol(jhcl)    *log_gamZ(jA,jhcl)
       gam(jA,ibin) = 10.**log_gam(jA)

       ! HCl
       jA = jhcl
       log_gam(jA) = xmol(jh2so4)  *log_gamZ(jA,jh2so4)  +   &
            xmol(jnh4hso4)*log_gamZ(jA,jnh4hso4)+   &
            xmol(jlvcite) *log_gamZ(jA,jlvcite) +   &
            xmol(jnh4so4) *log_gamZ(jA,jnh4so4) +   &
            xmol(jnahso4) *log_gamZ(jA,jnahso4) +   &
            xmol(jna3hso4)*log_gamZ(jA,jna3hso4)+   &
            xmol(jna2so4) *log_gamZ(jA,jna2so4) +   &
            xmol(jhno3)   *log_gamZ(jA,jhno3)   +   &
            xmol(jhcl)    *log_gamZ(jA,jhcl)
       gam(jA,ibin) = 10.**log_gam(jA)

20     gam(jnh4no3,ibin) = 1.0
       gam(jnh4cl,ibin)  = 1.0
       gam(jnano3,ibin)  = 1.0
       gam(jnacl,ibin)   = 1.0
       gam(jcano3,ibin)  = 1.0
       gam(jcacl2,ibin)  = 1.0

       gam(jnh4msa,ibin) = 1.0
       gam(jnamsa,ibin)  = 1.0
       gam(jcamsa2,ibin) = 1.0

       ! compute equilibrium pH
       ! cation molalities (mol/kg water)
       mc(jc_ca,ibin)   = 1.e-9*aer(ica_a,jliquid,ibin)/water_a(ibin)
       mc(jc_nh4,ibin)  = 1.e-9*aer(inh4_a,jliquid,ibin)/water_a(ibin)
       mc(jc_na,ibin)   = 1.e-9*aer(ina_a, jliquid,ibin)/water_a(ibin)

       ! anion molalities (mol/kg water)
       mSULF            = 1.e-9*aer(iso4_a,jliquid,ibin)/water_a(ibin)
       ma(ja_hso4,ibin) = 0.0
       ma(ja_so4,ibin)  = 0.0
       ma(ja_no3,ibin)  = 1.e-9*aer(ino3_a,jliquid,ibin)/water_a(ibin)
       ma(ja_cl,ibin)   = 1.e-9*aer(icl_a, jliquid,ibin)/water_a(ibin)
       ma(ja_msa,ibin)  = 1.e-9*aer(imsa_a,jliquid,ibin)/water_a(ibin)

       gam_ratio(ibin)  = gam(jnh4hso4,ibin)**2/gam(jhhso4,ibin)**2
       dumK = Keq_ll(1)*gam(jhhso4,ibin)**2/gam(jh2so4,ibin)**3

       c_bal =  mc(jc_nh4,ibin) + mc(jc_na,ibin) + 2.*mc(jc_ca,ibin)   &
            - ma(ja_no3,ibin) - ma(ja_cl,ibin) - mSULF - ma(ja_msa,ibin)

       aq = 1.0
       bq = dumK + c_bal
       cq = dumK*(c_bal - mSULF)

       !--quadratic solution
       if(bq .ne. 0.0)then
          xq = 4.*(1./bq)*(cq/bq)
       else
          xq = 1.e+6
       endif

       if(abs(xq) .lt. 1.e-6)then
          dum = xq*(0.5 + xq*(0.125 + xq*0.0625))
          quad = (-0.5*bq/aq)*dum
          if(quad .lt. 0.)then
             quad = -bq/aq - quad
          endif
       else
          quad = 0.5*(-bq+sqrt(bq*bq - 4.*cq))
       endif
       !--end of quadratic solution

       mc(jc_h,ibin) = max(quad, 1.d-7)
       ma(ja_so4,ibin) = mSULF*dumK/(mc(jc_h,ibin) + dumK)
       ma(ja_hso4,ibin)= mSULF - ma(ja_so4,ibin)

       activity(jcamsa2,ibin) = mc(jc_ca,ibin) * ma(ja_msa,ibin)**2 *   &
            gam(jcamsa2,ibin)**3

       activity(jnh4so4,ibin) = mc(jc_nh4,ibin)**2 * ma(ja_so4,ibin) *   &
            gam(jnh4so4,ibin)**3

       activity(jlvcite,ibin) = mc(jc_nh4,ibin)**3 * ma(ja_hso4,ibin) *   &
            ma(ja_so4,ibin) * gam(jlvcite,ibin)**5

       activity(jnh4hso4,ibin)= mc(jc_nh4,ibin) * ma(ja_hso4,ibin) *   &
            gam(jnh4hso4,ibin)**2

       activity(jnh4msa,ibin) = mc(jc_nh4,ibin) * ma(ja_msa,ibin) *   &
            gam(jnh4msa,ibin)**2

       activity(jna2so4,ibin) = mc(jc_na,ibin)**2 * ma(ja_so4,ibin) *   &
            gam(jna2so4,ibin)**3

       activity(jnahso4,ibin) = mc(jc_na,ibin) * ma(ja_hso4,ibin) *   &
            gam(jnahso4,ibin)**2

       activity(jnamsa,ibin)  = mc(jc_na,ibin) * ma(ja_msa,ibin) *   &
            gam(jnamsa,ibin)**2

       !      activity(jna3hso4,ibin)= mc(jc_na,ibin)**3 * ma(ja_hso4,ibin) *
       !     &                         ma(ja_so4,ibin) * gam(jna3hso4,ibin)**5

       activity(jna3hso4,ibin)= 0.0

       activity(jhno3,ibin)   = mc(jc_h,ibin) * ma(ja_no3,ibin) *   &
            gam(jhno3,ibin)**2

       activity(jhcl,ibin)    = mc(jc_h,ibin) * ma(ja_cl,ibin) *   &
            gam(jhcl,ibin)**2

       activity(jmsa,ibin)    = mc(jc_h,ibin) * ma(ja_msa,ibin) *   &
            gam(jmsa,ibin)**2

       ! sulfate-poor species
       activity(jnh4no3,ibin) = 0.0

       activity(jnh4cl,ibin)  = 0.0

       activity(jnano3,ibin)  = 0.0

       activity(jnacl,ibin)   = 0.0

       activity(jcano3,ibin)  = 0.0

       activity(jcacl2,ibin)  = 0.0

    endif
    return
  end subroutine compute_activities

  !***********************************************************************
  ! part of MESA: checks MESA convergence
  !
  ! author: Rahul A. Zaveri
  ! update: jan 2005
  !-----------------------------------------------------------------------
  subroutine MESA_convergence_criterion(ibin,iconverge_mass,iconverge_flux,        &
       idissolved,aer,electrolyte,mass_dry_a,mass_dry_salt,                        &
       mw_electrolyte,mw_aer_mac,flux_sl,phi_salt,rtol_mesa)  ! TOUCH

    use module_data_mosaic_aero, only: r8,nbin_a_max,naer,nelectrolyte,nsalt,      &!Parameters
         jsolid, mYES, xhyst_up_crustal_thresh,                                    &!Parameters
         mno,ioin_a,jcaso4,jcaco3         !TBD

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    integer, intent(inout) :: iconverge_mass, iconverge_flux, idissolved
    real(r8), intent(in) :: rtol_mesa
    real(r8), intent(inout), dimension(nsalt) :: flux_sl,phi_salt
    real(r8), intent(in),    dimension(nbin_a_max) :: mass_dry_a
    real(r8), intent(inout), dimension(nbin_a_max) :: mass_dry_salt
    real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
    real(r8), intent(in), dimension(naer) :: mw_aer_mac
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    ! local variables
    integer je, js, iaer
    real(r8) :: mass_solid, mass_solid_salt,frac_solid, XT, H_ion,   &
         crustal_solids, sumflux

    idissolved = mno             ! default = not completely dissolved

    ! check mass convergence
    iconverge_mass = mNO ! default value = no convergence

    !      call electrolytes_to_ions(jsolid,ibin,aer,electrolyte)
    !      mass_solid = 0.0
    !      do iaer = 1, naer
    !        mass_solid = mass_solid +
    !     &               aer(iaer,jsolid,ibin)*mw_aer_mac(iaer)*1.e-15  ! g/cc(air)
    !      enddo

    mass_solid_salt = 0.0
    do je = 1, nsalt
       mass_solid_salt = mass_solid_salt +   &
            electrolyte(je,jsolid,ibin)*mw_electrolyte(je)*1.e-15        ! g/cc(air)
    enddo

    !      frac_solid = mass_solid/mass_dry_a(ibin)

    if(mass_dry_salt(ibin) .le. 0.0)then
      frac_solid = 0.0
    else
      frac_solid = mass_solid_salt/mass_dry_salt(ibin)
    endif

    if(frac_solid .ge. 0.98)then
       iconverge_mass = mYES
       return
    endif

    ! check relative driving force convergence
    iconverge_flux = mYES
    do js = 1, nsalt
       if(abs(phi_salt(js)).gt. rtol_mesa)then
          iconverge_flux = mNO
          return
       endif
    enddo

    ! check if all the fluxes are zero

    sumflux = 0.0
    do js = 1, nsalt
       sumflux = sumflux + abs(flux_sl(js))
    enddo

!   crustal_solids = electrolyte(jcaco3,jsolid,ibin) +   &
!        electrolyte(jcaso4,jsolid,ibin) +   &
!        aer(ioin_a,jsolid,ibin)
    crustal_solids = electrolyte(jcaco3,jsolid,ibin)*mw_electrolyte(jcaco3) +  &
                     electrolyte(jcaso4,jsolid,ibin)*mw_electrolyte(jcaso4) +  &
                     aer(ioin_a,jsolid,ibin)*mw_aer_mac(ioin_a)

!   if(sumflux .eq. 0.0 .and. crustal_solids .eq. 0.0)then
    if ( sumflux .eq. 0.0 .and. &
         crustal_solids .le. xhyst_up_crustal_thresh*(mass_dry_a(ibin)*1.0e15) ) then
       ! crustal_solids is ng/m^3, mass_dry_a is g/cm^3
       idissolved = myes
    endif

    return
  end subroutine MESA_convergence_criterion

  !***********************************************************************
  ! computes ions from electrolytes
  !
  ! author: Rahul A. Zaveri
  ! update: jan 2005
  !-----------------------------------------------------------------------
  subroutine electrolytes_to_ions(jp,ibin,aer,electrolyte)

    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &!Parameters
         jh2so4,jnh4hso4,jlvcite,jnh4so4,jnahso4,jna3hso4,jna2so4,jcaso4,iso4_a,   &!TBD
         jhno3,jnh4no3,jcano3,jnano3,ino3_a,jhcl,jnh4cl,jcacl2,jnacl,icl_a,jmsa,   &!TBD
         jcamsa2,jnamsa,jnh4msa,imsa_a,jcaco3,ico3_a,ica_a,ina_a,inh4_a             !TBD
    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    ! local variables
    real(r8) :: sum_dum

    aer(iso4_a,jp,ibin) = electrolyte(jcaso4,jp,ibin)  +   &
         electrolyte(jna2so4,jp,ibin) +   &
         2.*electrolyte(jna3hso4,jp,ibin)+   &
         electrolyte(jnahso4,jp,ibin) +   &
         electrolyte(jnh4so4,jp,ibin) +   &
         2.*electrolyte(jlvcite,jp,ibin) +   &
         electrolyte(jnh4hso4,jp,ibin)+   &
         electrolyte(jh2so4,jp,ibin)

    aer(ino3_a,jp,ibin) = electrolyte(jnano3,jp,ibin)  +   &
         2.*electrolyte(jcano3,jp,ibin)  +   &
         electrolyte(jnh4no3,jp,ibin) +   &
         electrolyte(jhno3,jp,ibin)

    aer(icl_a,jp,ibin)  = electrolyte(jnacl,jp,ibin)   +   &
         2.*electrolyte(jcacl2,jp,ibin)  +   &
         electrolyte(jnh4cl,jp,ibin)  +   &
         electrolyte(jhcl,jp,ibin)

    aer(imsa_a,jp,ibin) = electrolyte(jnh4msa,jp,ibin) +   &
         electrolyte(jnamsa,jp,ibin)  +   &
         2.*electrolyte(jcamsa2,jp,ibin) +   &
         electrolyte(jmsa,jp,ibin)

    aer(ico3_a,jp,ibin) = electrolyte(jcaco3,jp,ibin)

    aer(ica_a,jp,ibin)  = electrolyte(jcaso4,jp,ibin)  +   &
         electrolyte(jcano3,jp,ibin)  +   &
         electrolyte(jcacl2,jp,ibin)  +   &
         electrolyte(jcaco3,jp,ibin)  +   &
         electrolyte(jcamsa2,jp,ibin)

    aer(ina_a,jp,ibin)  = electrolyte(jnano3,jp,ibin)  +   &
         electrolyte(jnacl,jp,ibin)   +   &
         2.*electrolyte(jna2so4,jp,ibin) +   &
         3.*electrolyte(jna3hso4,jp,ibin)+   &
         electrolyte(jnahso4,jp,ibin) +   &
         electrolyte(jnamsa,jp,ibin)

    aer(inh4_a,jp,ibin) = electrolyte(jnh4no3,jp,ibin) +   &
         electrolyte(jnh4cl,jp,ibin)  +   &
         2.*electrolyte(jnh4so4,jp,ibin) +   &
         3.*electrolyte(jlvcite,jp,ibin) +   &
         electrolyte(jnh4hso4,jp,ibin)+   &
         electrolyte(jnh4msa,jp,ibin)

    return
  end subroutine electrolytes_to_ions

  !***********************************************************************
  ! combinatorial method for computing electrolytes from ions
  !
  ! notes:
  !  - to be used for liquid-phase or total-phase only
  !  - transfers caso4 and caco3 from liquid to solid phase
  !
  ! author: Rahul A. Zaveri (based on code provided by A.S. Wexler)
  ! update: apr 2005
  !-----------------------------------------------------------------------
  subroutine ions_to_electrolytes(jp,ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma,    &
       nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a)

    use module_data_mosaic_aero, only: r8,naer,nbin_a_max,nelectrolyte,ncation,    &!Parameters
         nanion,jliquid,jsolid,                                                    &!Parameters
         ica_a,iso4_a,jcaso4,imsa_a,ina_a,inh4_a,ja_hso4,ja_so4,ino3_a,ja_no3,     &
         icl_a,ja_cl,ja_msa,jc_ca,jc_na,jc_nh4,jc_h,jna2so4,jnahso4,jnamsa,jnano3, &
         jnacl,jnh4so4,jnh4hso4,jnh4msa,jnh4no3,jnh4cl,jcano3,jcacl2,jcamsa2,      &
         jh2so4,jhno3,jhcl,jmsa,jlvcite,jna3hso4                                    !TBD
    implicit none

    ! subr arguments
    integer, intent(in) :: ibin, jp
    real(r8), intent(inout) :: XT
    real(r8), intent(in), dimension(Ncation) :: zc,MW_c
    real(r8), intent(inout), dimension(Ncation) :: nc_Mc,xeq_c
    real(r8), intent(in), dimension(Nanion) :: za,MW_a
    real(r8), intent(inout), dimension(Nanion) :: xeq_a,na_Ma
    real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    ! local variables
    character(len=256) :: errmsg
    integer iaer, je, jc, ja, icase
    real(r8) :: store(naer), sum_dum, sum_naza, sum_nczc, sum_na_nh4,   &
         f_nh4, f_na, xh, xb, xl, xs, cat_net, rem_nh4, rem_na
    real(r8) :: nc(ncation), na(nanion)

    if(jp .ne. jliquid)then
       call wrf_message(' jp must be jliquid')
       call wrf_message(' in ions_to_electrolytes sub')
       write(errmsg,*)' wrong jp = ', jp
       call wrf_error_fatal(trim(adjustl(errmsg)))

    endif

    ! remove negative concentrations, if any
    !      do iaer = 1, naer
    !        aer(iaer,jp,ibin) = max(0.0d0, aer(iaer,jp,ibin))    ! EFFI
    !      enddo

    ! first transfer caso4 from liquid to solid phase (caco3 should not be present here)
    store(ica_a)  = aer(ica_a, jp,ibin)
    store(iso4_a) = aer(iso4_a,jp,ibin)

    call form_caso4(store,jp,ibin,electrolyte)

    if(jp .eq. jliquid)then ! transfer caso4 from liquid to solid phase
       aer(ica_a,jliquid,ibin) = aer(ica_a,jliquid,ibin) -   &
            electrolyte(jcaso4,jliquid,ibin)

       aer(iso4_a,jliquid,ibin)= aer(iso4_a,jliquid,ibin)-   &
            electrolyte(jcaso4,jliquid,ibin)

       aer(ica_a,jsolid,ibin)  = aer(ica_a,jsolid,ibin) +   &
            electrolyte(jcaso4,jliquid,ibin)

       aer(iso4_a,jsolid,ibin) = aer(iso4_a,jsolid,ibin) +   &
            electrolyte(jcaso4,jliquid,ibin)

       electrolyte(jcaso4,jsolid,ibin)=electrolyte(jcaso4,jsolid,ibin)   &
            +electrolyte(jcaso4,jliquid,ibin)
       electrolyte(jcaso4,jliquid,ibin)= 0.0
    endif

    ! calculate sulfate ratio
    !      call calculate_XT(ibin,jp,XT,aer)              ! EFFI

    if( (aer(iso4_a,jp,ibin)+aer(imsa_a,jp,ibin)) .gt.0.0)then
       XT   = ( aer(inh4_a,jp,ibin) +   &
            aer(ina_a,jp,ibin)  +   &
            2.*aer(ica_a,jp,ibin) )/   &
            (aer(iso4_a,jp,ibin)+0.5*aer(imsa_a,jp,ibin))
    else
       XT   = -1.0
    endif

!    if(XT.ge.1.9999 .or. XT.lt.0.)then ! commented out by RAZ 4/16/2014
    if(XT.ge.2.0 .or. XT.lt.0.)then     ! Slightly different logic, consistent with that in compute_activities subr. RAZ 4/16/2014
       icase = 1  ! sulfate poor: near neutral (acidity is caused by HCl and/or HNO3)
    else
       icase = 2  ! sulfate rich: acidic (acidity is caused by excess SO4)
    endif

    ! initialize to zero
    do je = 1, nelectrolyte
       electrolyte(je,jp,ibin) = 0.0
    enddo

    !
    !---------------------------------------------------------
    ! initialize moles of ions depending on the sulfate domain

    if(icase.eq.1)then ! XT >= 2 or XT < 0: SULFATE POOR (OR NO SULFATE) DOMAIN. RAZ 4/16/2014

       na(ja_hso4)= 0.0
       na(ja_so4) = aer(iso4_a,jp,ibin)
       na(ja_no3) = aer(ino3_a,jp,ibin)
       na(ja_cl)  = aer(icl_a, jp,ibin)
       na(ja_msa) = aer(imsa_a,jp,ibin)

       nc(jc_ca)  = aer(ica_a, jp,ibin)
       nc(jc_na)  = aer(ina_a, jp,ibin)
       nc(jc_nh4) = aer(inh4_a,jp,ibin)

       cat_net = (   &
            (2.*na(ja_so4)+na(ja_no3)+na(ja_cl)+na(ja_msa)) -   &
            (2.*nc(jc_ca) +nc(jc_nh4)+nc(jc_na)) )

       if(cat_net .lt. 0.0)then

          nc(jc_h) = 0.0

       else  ! cat_net must be 0.0 or positive

          nc(jc_h) = cat_net

       endif

       ! now compute equivalent fractions
       sum_naza = 0.0
       do ja = 1, nanion
          sum_naza = sum_naza + na(ja)*za(ja)
       enddo

       sum_nczc = 0.0
       do jc = 1, ncation
          sum_nczc = sum_nczc + nc(jc)*zc(jc)
       enddo

       if(sum_naza .eq. 0. .or. sum_nczc .eq. 0.)then ! it's ok. this may happen if the aerosol is assumed to be composed of hygroscopic SOA, POA, BC, OIN, but does not contain any inorganic electrolytes
!          write(6,*)'ionic concentrations are zero in ibin', ibin  ! commented out by RAZ 4/16/2014
!          write(6,*)'sum_naza = ', sum_naza          ! commented out by RAZ 4/16/2014
!          write(6,*)'sum_nczc = ', sum_nczc          ! commented out by RAZ 4/16/2014
          return
       endif

       do ja = 1, nanion
          xeq_a(ja) = na(ja)*za(ja)/sum_naza
       enddo

       do jc = 1, ncation
          xeq_c(jc) = nc(jc)*zc(jc)/sum_nczc
       enddo

       na_Ma(ja_so4) = na(ja_so4) *MW_a(ja_so4)
       na_Ma(ja_no3) = na(ja_no3) *MW_a(ja_no3)
       na_Ma(ja_cl)  = na(ja_cl)  *MW_a(ja_cl)
       na_Ma(ja_msa) = na(ja_msa) *MW_a(ja_msa)
       na_Ma(ja_hso4)= na(ja_hso4)*MW_a(ja_hso4)

       nc_Mc(jc_ca)  = nc(jc_ca) *MW_c(jc_ca)
       nc_Mc(jc_na)  = nc(jc_na) *MW_c(jc_na)
       nc_Mc(jc_nh4) = nc(jc_nh4)*MW_c(jc_nh4)
       nc_Mc(jc_h)   = nc(jc_h)  *MW_c(jc_h)

       ! now compute electrolyte moles
       if(xeq_c(jc_na) .gt. 0. .and. xeq_a(ja_so4) .gt. 0.)then
          electrolyte(jna2so4,jp,ibin) = (xeq_c(jc_na) *na_Ma(ja_so4) +   &
               xeq_a(ja_so4)*nc_Mc(jc_na))/   &
               mw_electrolyte(jna2so4)
       endif

       electrolyte(jnahso4,jp,ibin) = 0.0

       if(xeq_c(jc_na) .gt. 0. .and. xeq_a(ja_msa) .gt. 0.)then
          electrolyte(jnamsa,jp,ibin)  = (xeq_c(jc_na) *na_Ma(ja_msa) +   &
               xeq_a(ja_msa)*nc_Mc(jc_na))/   &
               mw_electrolyte(jnamsa)
       endif

       if(xeq_c(jc_na) .gt. 0. .and. xeq_a(ja_no3) .gt. 0.)then
          electrolyte(jnano3,jp,ibin)  = (xeq_c(jc_na) *na_Ma(ja_no3) +   &
               xeq_a(ja_no3)*nc_Mc(jc_na))/   &
               mw_electrolyte(jnano3)
       endif

       if(xeq_c(jc_na) .gt. 0. .and. xeq_a(ja_cl) .gt. 0.)then
          electrolyte(jnacl,jp,ibin)   = (xeq_c(jc_na) *na_Ma(ja_cl) +   &
               xeq_a(ja_cl) *nc_Mc(jc_na))/   &
               mw_electrolyte(jnacl)
       endif

       if(xeq_c(jc_nh4) .gt. 0. .and. xeq_a(ja_so4) .gt. 0.)then
          electrolyte(jnh4so4,jp,ibin) = (xeq_c(jc_nh4)*na_Ma(ja_so4) +   &
               xeq_a(ja_so4)*nc_Mc(jc_nh4))/   &
               mw_electrolyte(jnh4so4)
       endif

       electrolyte(jnh4hso4,jp,ibin)= 0.0

       if(xeq_c(jc_nh4) .gt. 0. .and. xeq_a(ja_msa) .gt. 0.)then
          electrolyte(jnh4msa,jp,ibin) = (xeq_c(jc_nh4)*na_Ma(ja_msa) +   &
               xeq_a(ja_msa)*nc_Mc(jc_nh4))/   &
               mw_electrolyte(jnh4msa)
       endif

       if(xeq_c(jc_nh4) .gt. 0. .and. xeq_a(ja_no3) .gt. 0.)then
          electrolyte(jnh4no3,jp,ibin) = (xeq_c(jc_nh4)*na_Ma(ja_no3) +   &
               xeq_a(ja_no3)*nc_Mc(jc_nh4))/   &
               mw_electrolyte(jnh4no3)
       endif

       if(xeq_c(jc_nh4) .gt. 0. .and. xeq_a(ja_cl) .gt. 0.)then
          electrolyte(jnh4cl,jp,ibin)  = (xeq_c(jc_nh4)*na_Ma(ja_cl) +   &
               xeq_a(ja_cl) *nc_Mc(jc_nh4))/   &
               mw_electrolyte(jnh4cl)
       endif

       if(xeq_c(jc_ca) .gt. 0. .and. xeq_a(ja_no3) .gt. 0.0)then
          electrolyte(jcano3, jp,ibin) = (xeq_c(jc_ca) *na_Ma(ja_no3) +   &
               xeq_a(ja_no3)*nc_Mc(jc_ca))/   &
               mw_electrolyte(jcano3)
       endif

       if(xeq_c(jc_ca) .gt. 0. .and. xeq_a(ja_cl) .gt. 0.)then
          electrolyte(jcacl2,jp,ibin)  = (xeq_c(jc_ca) *na_Ma(ja_cl) +   &
               xeq_a(ja_cl) *nc_Mc(jc_ca))/   &
               mw_electrolyte(jcacl2)
       endif

       if(xeq_c(jc_ca) .gt. 0. .and. xeq_a(ja_msa) .gt. 0.)then
          electrolyte(jcamsa2,jp,ibin) = (xeq_c(jc_ca) *na_Ma(ja_msa) +   &
               xeq_a(ja_msa) *nc_Mc(jc_ca))/   &
               mw_electrolyte(jcamsa2)
       endif

       electrolyte(jh2so4, jp,ibin) = 0.0

       if(xeq_c(jc_h) .gt. 0. .and. xeq_a(ja_no3) .gt. 0.)then
          electrolyte(jhno3,jp,ibin)     = (xeq_c(jc_h)  *na_Ma(ja_no3) +   &
               xeq_a(ja_no3)*nc_Mc(jc_h))/   &
               mw_electrolyte(jhno3)
       endif

       if(xeq_c(jc_h) .gt. 0. .and. xeq_a(ja_cl) .gt. 0.)then
          electrolyte(jhcl,jp,ibin)    = (xeq_c(jc_h) *na_Ma(ja_cl) +   &
               xeq_a(ja_cl)*nc_Mc(jc_h))/   &
               mw_electrolyte(jhcl)
       endif

       if(xeq_c(jc_h) .gt. 0. .and. xeq_a(ja_msa) .gt. 0.)then
          electrolyte(jmsa,jp,ibin)    = (xeq_c(jc_h) *na_Ma(ja_msa) +   &
               xeq_a(ja_msa)*nc_Mc(jc_h))/   &
               mw_electrolyte(jmsa)
       endif

       !--------------------------------------------------------------------

    elseif(icase.eq.2)then ! XT < 2 : SULFATE RICH DOMAIN

       store(imsa_a) = aer(imsa_a,jp,ibin)
       store(ica_a)  = aer(ica_a, jp,ibin)

       call form_camsa2(store,jp,ibin,electrolyte)

       sum_na_nh4 = aer(ina_a,jp,ibin) + aer(inh4_a,jp,ibin)

       if(sum_na_nh4 .gt. 0.0)then
          f_na  = aer(ina_a,jp,ibin)/sum_na_nh4
          f_nh4 = aer(inh4_a,jp,ibin)/sum_na_nh4
       else
          f_na  = 0.0
          f_nh4 = 0.0
       endif

       ! first form msa electrolytes
       if(sum_na_nh4 .gt. store(imsa_a))then
          electrolyte(jnamsa,jp,ibin)  = f_na *store(imsa_a)
          electrolyte(jnh4msa,jp,ibin) = f_nh4*store(imsa_a)
          rem_na = max(0.0_r8, aer(ina_a,jp,ibin) - electrolyte(jnamsa,jp,ibin))  ! remaining na  RAZ 4/16/2014
          rem_nh4= max(0.0_r8, aer(inh4_a,jp,ibin)- electrolyte(jnh4msa,jp,ibin)) ! remaining nh4 RAZ 4/16/2014
       else
          electrolyte(jnamsa,jp,ibin)  = aer(ina_a,jp,ibin)
          electrolyte(jnh4msa,jp,ibin) = aer(inh4_a,jp,ibin)
          electrolyte(jmsa,jp,ibin)    = max(0.0_r8, store(imsa_a) - sum_na_nh4) ! RAZ 4/16/2014
          rem_nh4 = 0.0  ! remaining nh4
          rem_na  = 0.0  ! remaining na
       endif

       ! recompute XT
       if(aer(iso4_a,jp,ibin).gt.0.0)then
          XT = (rem_nh4 + rem_na)/aer(iso4_a,jp,ibin)
       else
          goto 10
       endif

       if(XT .le. 1.0)then            ! h2so4 + bisulfate
          xh = max(0.0_r8, (1.0_r8 - XT))   ! RAZ 4/16/2014
          xb = XT
          electrolyte(jh2so4,jp,ibin)   = xh*aer(iso4_a,jp,ibin)
          electrolyte(jnh4hso4,jp,ibin) = xb*f_nh4*aer(iso4_a,jp,ibin)
          electrolyte(jnahso4,jp,ibin)  = xb*f_na *aer(iso4_a,jp,ibin)
       elseif(XT .le. 1.5)then    ! bisulfate + letovicite
          xb = max(0.0_r8, 3.0_r8 - 2.0_r8*XT) ! RAZ 4/16/2014
          xl = max(0.0_r8, XT - 1.0_r8)     ! RAZ 4/16/2014
          electrolyte(jnh4hso4,jp,ibin) = xb*f_nh4*aer(iso4_a,jp,ibin)
          electrolyte(jnahso4,jp,ibin)  = xb*f_na *aer(iso4_a,jp,ibin)
          electrolyte(jlvcite,jp,ibin)  = xl*f_nh4*aer(iso4_a,jp,ibin)
          electrolyte(jna3hso4,jp,ibin) = xl*f_na *aer(iso4_a,jp,ibin)
       else                       ! letovicite + sulfate
          xl = max(0.0_r8, 2.0_r8 - XT)     ! RAZ 4/16/2014
          xs = max(0.0_r8, 2.0_r8*XT - 3.0_r8) ! RAZ 4/16/2014
          electrolyte(jlvcite,jp,ibin)  = xl*f_nh4*aer(iso4_a,jp,ibin)
          electrolyte(jna3hso4,jp,ibin) = xl*f_na *aer(iso4_a,jp,ibin)
          electrolyte(jnh4so4,jp,ibin)  = xs*f_nh4*aer(iso4_a,jp,ibin)
          electrolyte(jna2so4,jp,ibin)  = xs*f_na *aer(iso4_a,jp,ibin)
       endif

       electrolyte(jhno3,jp,ibin) = aer(ino3_a,jp,ibin)
       electrolyte(jhcl,jp,ibin)  = aer(icl_a,jp,ibin)

    endif
    !---------------------------------------------------------
    !
    ! calculate % composition  EFFI
10  sum_dum = 0.0
    !!      do je = 1, nelectrolyte
    !!        sum_dum = sum_dum + electrolyte(je,jp,ibin)
    !!      enddo
    !!
    !!      if(sum_dum .eq. 0.)sum_dum = 1.0
    !!      electrolyte_sum(jp,ibin) = sum_dum
    !!
    !!      do je = 1, nelectrolyte
    !!        epercent(je,jp,ibin) = 100.*electrolyte(je,jp,ibin)/sum_dum
    !!      enddo
    !!

    return
  end subroutine ions_to_electrolytes

  !***********************************************************************
  ! part of MESA: calculates liquid electrolytes from ions
  !
  ! notes:
  !  - this subroutine is to be used for liquid-phase or total-phase only
  !  - this sub transfers caso4 and caco3 from liquid to solid phase
  !
  ! author: Rahul A. Zaveri
  ! update: jan 2005
  !-----------------------------------------------------------------------
  subroutine MESA_estimate_eleliquid(ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma,    &
       nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a,eleliquid)    ! TOUCH
    use module_data_mosaic_aero, only: r8,naer,nbin_a_max,nelectrolyte,ncation,    &!Parameters
         nanion,jliquid,                                                           &!Parameters
         jh2so4,jhno3,jhcl,jmsa,jlvcite,jnh4no3,jnh4cl,jcamsa2,jcano3,jcacl2,      &
         jnano3,jnacl,jnh4so4,jnh4hso4,jnh4msa,jna2so4,jnahso4,jnamsa,iso4_a,      &
         ja_so4,ja_no3,ja_cl,imsa_a,ja_msa,jc_ca,ina_a,jc_na,inh4_a,jc_nh4,jc_h,   &
         ica_a,ino3_a,icl_a,ja_hso4

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    real(r8), intent(inout) :: XT
    real(r8), intent(in), dimension(Ncation) :: zc,MW_c
    real(r8), intent(inout), dimension(Ncation) :: nc_Mc,xeq_c
    real(r8), intent(in), dimension(Nanion) :: za,MW_a
    real(r8), intent(inout), dimension(Nanion) :: xeq_a,na_Ma
    real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
    real(r8), intent(inout), dimension(nelectrolyte) :: eleliquid
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    ! local variables
    integer iaer, je, jc, ja, icase, jp
    real(r8) :: store(naer), sum_dum, sum_naza, sum_nczc, sum_na_nh4,   &
         f_nh4, f_na, xh, xb, xl, xs, XT_d, XNa_d, XNH4_d,   &
         xdum, dum, cat_net
    real(r8) :: nc(ncation), na(nanion)
    real(r8) :: dum_ca, dum_no3, dum_cl, cano3, cacl2

    !nc(:) = 0.0_r8!BSINGH - initialized to zero

    ! remove negative concentrations, if any
    do iaer =  1, naer
       aer(iaer,jliquid,ibin) = max(0.0d0, aer(iaer,jliquid,ibin))
    enddo

    ! calculate sulfate ratio
    call calculate_XT(ibin,jliquid,XT,aer)

    if(XT .ge. 2.0 .or. XT.lt.0.)then
       icase = 1 ! near neutral (acidity is caused by HCl and/or HNO3)
    else
       icase = 2 ! acidic (acidity is caused by excess SO4)
    endif

    ! initialize to zero
    do je = 1, nelectrolyte
       eleliquid(je) = 0.0
    enddo

    !
    !---------------------------------------------------------
    ! initialize moles of ions depending on the sulfate domain

    jp = jliquid

    if(icase.eq.1)then ! XT >= 2 : SULFATE POOR DOMAIN

       dum_ca  = aer(ica_a,jp,ibin)
       dum_no3 = aer(ino3_a,jp,ibin)
       dum_cl  = aer(icl_a,jp,ibin)

       cano3   = min(dum_ca, 0.5*dum_no3)
       dum_ca  = max(0.d0, dum_ca - cano3)
       dum_no3 = max(0.d0, dum_no3 - 2.*cano3)

       cacl2   = min(dum_ca, 0.5*dum_cl)
       dum_ca  = max(0.d0, dum_ca - cacl2)
       dum_cl  = max(0.d0, dum_cl - 2.*cacl2)

       na(ja_hso4)= 0.0
       na(ja_so4) = aer(iso4_a,jp,ibin)
       na(ja_no3) = aer(ino3_a,jp,ibin)
       na(ja_cl)  = aer(icl_a, jp,ibin)
       na(ja_msa) = aer(imsa_a,jp,ibin)

       nc(jc_ca)  = aer(ica_a, jp,ibin)
       nc(jc_na)  = aer(ina_a, jp,ibin)
       nc(jc_nh4) = aer(inh4_a,jp,ibin)

       cat_net = (   &
            (2.*na(ja_so4)+na(ja_no3)+na(ja_cl)+na(ja_msa)) -   &
            (2.*nc(jc_ca) +nc(jc_nh4)+nc(jc_na)) )   ! RAZ 11/11/2014: bug fix. remove nc(jc_h)

       if(cat_net .lt. 0.0)then

          nc(jc_h) = 0.0

       else  ! cat_net must be 0.0 or positive

          nc(jc_h) = cat_net

       endif

       ! now compute equivalent fractions
       sum_naza = 0.0
       do ja = 1, nanion
          sum_naza = sum_naza + na(ja)*za(ja)
       enddo

       sum_nczc = 0.0
       do jc = 1, ncation
          sum_nczc = sum_nczc + nc(jc)*zc(jc)
       enddo

       if(sum_naza .eq. 0. .or. sum_nczc .eq. 0.)then
          write(6,*)'ionic concentrations are zero in ibin', ibin
          write(6,*)'sum_naza = ', sum_naza
          write(6,*)'sum_nczc = ', sum_nczc
          return
       endif

       do ja = 1, nanion
          xeq_a(ja) = na(ja)*za(ja)/sum_naza
       enddo

       do jc = 1, ncation
          xeq_c(jc) = nc(jc)*zc(jc)/sum_nczc
       enddo

       na_Ma(ja_so4) = na(ja_so4) *MW_a(ja_so4)
       na_Ma(ja_no3) = na(ja_no3) *MW_a(ja_no3)
       na_Ma(ja_cl)  = na(ja_cl)  *MW_a(ja_cl)
       na_Ma(ja_hso4)= na(ja_hso4)*MW_a(ja_hso4)
       na_Ma(ja_msa) = na(ja_msa) *MW_a(ja_msa)

       nc_Mc(jc_ca)  = nc(jc_ca) *MW_c(jc_ca)
       nc_Mc(jc_na)  = nc(jc_na) *MW_c(jc_na)
       nc_Mc(jc_nh4) = nc(jc_nh4)*MW_c(jc_nh4)
       nc_Mc(jc_h)   = nc(jc_h)  *MW_c(jc_h)

       ! now compute electrolyte moles
       eleliquid(jna2so4) = (xeq_c(jc_na) *na_Ma(ja_so4) +   &
            xeq_a(ja_so4)*nc_Mc(jc_na))/   &
            mw_electrolyte(jna2so4)

       eleliquid(jnahso4) = (xeq_c(jc_na) *na_Ma(ja_hso4) +   &
            xeq_a(ja_hso4)*nc_Mc(jc_na))/   &
            mw_electrolyte(jnahso4)

       eleliquid(jnamsa)  = (xeq_c(jc_na) *na_Ma(ja_msa) +   &
            xeq_a(ja_msa)*nc_Mc(jc_na))/   &
            mw_electrolyte(jnamsa)

       eleliquid(jnano3)  = (xeq_c(jc_na) *na_Ma(ja_no3) +   &
            xeq_a(ja_no3)*nc_Mc(jc_na))/   &
            mw_electrolyte(jnano3)

       eleliquid(jnacl)   = (xeq_c(jc_na) *na_Ma(ja_cl) +   &
            xeq_a(ja_cl) *nc_Mc(jc_na))/   &
            mw_electrolyte(jnacl)

       eleliquid(jnh4so4) = (xeq_c(jc_nh4)*na_Ma(ja_so4) +   &
            xeq_a(ja_so4)*nc_Mc(jc_nh4))/   &
            mw_electrolyte(jnh4so4)

       eleliquid(jnh4hso4)= (xeq_c(jc_nh4)*na_Ma(ja_hso4) +   &
            xeq_a(ja_hso4)*nc_Mc(jc_nh4))/   &
            mw_electrolyte(jnh4hso4)

       eleliquid(jnh4msa) = (xeq_c(jc_nh4) *na_Ma(ja_msa) +   &
            xeq_a(ja_msa)*nc_Mc(jc_nh4))/   &
            mw_electrolyte(jnh4msa)

       eleliquid(jnh4no3) = (xeq_c(jc_nh4)*na_Ma(ja_no3) +   &
            xeq_a(ja_no3)*nc_Mc(jc_nh4))/   &
            mw_electrolyte(jnh4no3)

       eleliquid(jnh4cl)  = (xeq_c(jc_nh4)*na_Ma(ja_cl) +   &
            xeq_a(ja_cl) *nc_Mc(jc_nh4))/   &
            mw_electrolyte(jnh4cl)

       eleliquid(jcamsa2) = (xeq_c(jc_ca) *na_Ma(ja_msa) +   &
            xeq_a(ja_msa)*nc_Mc(jc_ca))/   &
            mw_electrolyte(jcamsa2)

       eleliquid(jcano3)  = (xeq_c(jc_ca) *na_Ma(ja_no3) +   &
            xeq_a(ja_no3)*nc_Mc(jc_ca))/   &
            mw_electrolyte(jcano3)

       eleliquid(jcacl2)  = (xeq_c(jc_ca) *na_Ma(ja_cl) +   &
            xeq_a(ja_cl) *nc_Mc(jc_ca))/   &
            mw_electrolyte(jcacl2)

       eleliquid(jh2so4)  = (xeq_c(jc_h)   *na_Ma(ja_hso4) +   &
            xeq_a(ja_hso4)*nc_Mc(jc_h))/   &
            mw_electrolyte(jh2so4)

       eleliquid(jhno3)   = (xeq_c(jc_h)  *na_Ma(ja_no3) +   &
            xeq_a(ja_no3)*nc_Mc(jc_h))/   &
            mw_electrolyte(jhno3)

       eleliquid(jhcl)    = (xeq_c(jc_h) *na_Ma(ja_cl) +   &
            xeq_a(ja_cl)*nc_Mc(jc_h))/   &
            mw_electrolyte(jhcl)

       eleliquid(jmsa)    = (xeq_c(jc_h)  *na_Ma(ja_msa) +   &
            xeq_a(ja_msa)*nc_Mc(jc_h))/   &
            mw_electrolyte(jmsa)

       !--------------------------------------------------------------------

    elseif(icase.eq.2)then ! XT < 2 : SULFATE RICH DOMAIN

       jp = jliquid

       store(iso4_a) = aer(iso4_a,jp,ibin)
       store(imsa_a) = aer(imsa_a,jp,ibin)
       store(inh4_a) = aer(inh4_a,jp,ibin)
       store(ina_a)  = aer(ina_a, jp,ibin)
       store(ica_a)  = aer(ica_a, jp,ibin)

       call form_camsa2(store,jp,ibin,electrolyte)

       sum_na_nh4 = store(ina_a) + store(inh4_a)
       if(sum_na_nh4 .gt. 0.0)then
          f_nh4 = store(inh4_a)/sum_na_nh4
          f_na  = store(ina_a)/sum_na_nh4
       else
          f_nh4 = 0.0
          f_na  = 0.0
       endif

       ! first form msa electrolytes
       if(sum_na_nh4 .gt. store(imsa_a))then
          eleliquid(jnh4msa) = f_nh4*store(imsa_a)
          eleliquid(jnamsa)  = f_na *store(imsa_a)
          store(inh4_a)= store(inh4_a)-eleliquid(jnh4msa) ! remaining nh4
          store(ina_a) = store(ina_a) -eleliquid(jnamsa)  ! remaining na
       else
          eleliquid(jnh4msa) = store(inh4_a)
          eleliquid(jnamsa)  = store(ina_a)
          eleliquid(jmsa)    = store(imsa_a) - sum_na_nh4
          store(inh4_a)= 0.0  ! remaining nh4
          store(ina_a) = 0.0  ! remaining na
       endif

       if(store(iso4_a).eq.0.0)goto 10

       XT_d  = XT
       XNa_d = 1. + 0.5*store(ina_a)/store(iso4_a)
       xdum = store(iso4_a) - store(inh4_a)

       dum = ( (2.*store(iso4_a)) -   &
            (store(ina_a)) )
       if(store(inh4_a) .gt. 0.0 .and. dum .gt. 0.0)then
          XNH4_d = 2.*store(inh4_a)/   &
               (2.*store(iso4_a) - store(ina_a))
       else
          XNH4_d = 0.0
       endif

       IF(store(inh4_a) .gt. 0.0)THEN
          if(XT_d .ge. XNa_d)then
             eleliquid(jna2so4) = 0.5*store(ina_a)

             if(XNH4_d .ge. 5./3.)then
                eleliquid(jnh4so4) = 1.5*store(ina_a)   &
                     - 3.*xdum - store(inh4_a)
                eleliquid(jlvcite) = 2.*xdum + store(inh4_a)   &
                     - store(ina_a)
             elseif(XNH4_d .ge. 1.5)then
                eleliquid(jnh4so4) = store(inh4_a)/5.
                eleliquid(jlvcite) = store(inh4_a)/5.
             elseif(XNH4_d .ge. 1.0)then
                eleliquid(jnh4so4) = store(inh4_a)/6.
                eleliquid(jlvcite) = store(inh4_a)/6.
                eleliquid(jnh4hso4)= store(inh4_a)/6.
             endif

          elseif(XT_d .gt. 1.0)then
             eleliquid(jnh4so4)  = store(inh4_a)/6.
             eleliquid(jlvcite)  = store(inh4_a)/6.
             eleliquid(jnh4hso4) = store(inh4_a)/6.
             eleliquid(jna2so4)  = store(ina_a)/3.
             eleliquid(jnahso4)  = store(ina_a)/3.
          elseif(XT_d .le. 1.0)then
             eleliquid(jna2so4)  = store(ina_a)/4.
             eleliquid(jnahso4)  = store(ina_a)/2.
             eleliquid(jlvcite)  = store(inh4_a)/6.
             eleliquid(jnh4hso4) = store(inh4_a)/2.
          endif

       ELSE

          if(XT_d .gt. 1.0)then
             eleliquid(jna2so4) = store(ina_a) - store(iso4_a)
             eleliquid(jnahso4) = 2.*store(iso4_a) -   &
                  store(ina_a)
          else
             eleliquid(jna2so4) = store(ina_a)/4.
             eleliquid(jnahso4) = store(ina_a)/2.
          endif

       ENDIF

    endif
    !---------------------------------------------------------

10  return
  end subroutine MESA_estimate_eleliquid

  !***********************************************************************
  ! part of MESA: completely dissolves small amounts of soluble salts
  !
  ! author: Rahul A. Zaveri
  ! update: jan 2005
  !-----------------------------------------------------------------------
  subroutine MESA_dissolve_small_salt(ibin,js,aer,electrolyte)

    use module_data_mosaic_aero, only:r8,naer,nbin_a_max,nelectrolyte,jsolid,      &!Parameters
         jliquid,                                                                  &!Parameters
         jh2so4,jhno3,jhcl,jlvcite,jnh4no3,jnh4cl,jcamsa2,jcano3,jcacl2,jnano3,    &!TBD
         jnacl,jnh4so4,jnh4hso4,jnh4msa,jna2so4,jnahso4,jnamsa,iso4_a,ina_a,       &!TBD
         inh4_a,jna3hso4,jcaso4,jcaco3,ica_a,ino3_a,icl_a

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin, js
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    !Local variables
    integer :: jp

    jp = jsolid

    if(js .eq. jnh4so4)then
       aer(inh4_a,jliquid,ibin) = aer(inh4_a,jliquid,ibin) +   &
            2.*electrolyte(js,jsolid,ibin)
       aer(iso4_a,jliquid,ibin) = aer(iso4_a,jliquid,ibin) +   &
            electrolyte(js,jsolid,ibin)

       electrolyte(js,jsolid,ibin) = 0.0

       aer(inh4_a,jp,ibin) = electrolyte(jnh4no3,jp,ibin) +   &
            electrolyte(jnh4cl,jp,ibin)  +   &
            2.*electrolyte(jnh4so4,jp,ibin) +   &
            3.*electrolyte(jlvcite,jp,ibin) +   &
            electrolyte(jnh4hso4,jp,ibin)+   &
            electrolyte(jnh4msa,jp,ibin)

       aer(iso4_a,jp,ibin) = electrolyte(jcaso4,jp,ibin)  +   &
            electrolyte(jna2so4,jp,ibin) +   &
            2.*electrolyte(jna3hso4,jp,ibin)+   &
            electrolyte(jnahso4,jp,ibin) +   &
            electrolyte(jnh4so4,jp,ibin) +   &
            2.*electrolyte(jlvcite,jp,ibin) +   &
            electrolyte(jnh4hso4,jp,ibin)+   &
            electrolyte(jh2so4,jp,ibin)
       return
    endif

    if(js .eq. jlvcite)then
       aer(inh4_a,jliquid,ibin) = aer(inh4_a,jliquid,ibin) +   &
            3.*electrolyte(js,jsolid,ibin)
       aer(iso4_a,jliquid,ibin) = aer(iso4_a,jliquid,ibin) +   &
            2.*electrolyte(js,jsolid,ibin)

       electrolyte(js,jsolid,ibin) = 0.0

       aer(inh4_a,jp,ibin) = electrolyte(jnh4no3,jp,ibin) +   &
            electrolyte(jnh4cl,jp,ibin)  +   &
            2.*electrolyte(jnh4so4,jp,ibin) +   &
            3.*electrolyte(jlvcite,jp,ibin) +   &
            electrolyte(jnh4hso4,jp,ibin)+   &
            electrolyte(jnh4msa,jp,ibin)

       aer(iso4_a,jp,ibin) = electrolyte(jcaso4,jp,ibin)  +   &
            electrolyte(jna2so4,jp,ibin) +   &
            2.*electrolyte(jna3hso4,jp,ibin)+   &
            electrolyte(jnahso4,jp,ibin) +   &
            electrolyte(jnh4so4,jp,ibin) +   &
            2.*electrolyte(jlvcite,jp,ibin) +   &
            electrolyte(jnh4hso4,jp,ibin)+   &
            electrolyte(jh2so4,jp,ibin)
       return
    endif

    if(js .eq. jnh4hso4)then
       aer(inh4_a,jliquid,ibin) = aer(inh4_a,jliquid,ibin) +   &
            electrolyte(js,jsolid,ibin)
       aer(iso4_a,jliquid,ibin) = aer(iso4_a,jliquid,ibin) +   &
            electrolyte(js,jsolid,ibin)

       electrolyte(js,jsolid,ibin) = 0.0

       aer(inh4_a,jp,ibin) = electrolyte(jnh4no3,jp,ibin) +   &
            electrolyte(jnh4cl,jp,ibin)  +   &
            2.*electrolyte(jnh4so4,jp,ibin) +   &
            3.*electrolyte(jlvcite,jp,ibin) +   &
            electrolyte(jnh4hso4,jp,ibin)+   &
            electrolyte(jnh4msa,jp,ibin)

       aer(iso4_a,jp,ibin) = electrolyte(jcaso4,jp,ibin)  +   &
            electrolyte(jna2so4,jp,ibin) +   &
            2.*electrolyte(jna3hso4,jp,ibin)+   &
            electrolyte(jnahso4,jp,ibin) +   &
            electrolyte(jnh4so4,jp,ibin) +   &
            2.*electrolyte(jlvcite,jp,ibin) +   &
            electrolyte(jnh4hso4,jp,ibin)+   &
            electrolyte(jh2so4,jp,ibin)
       return
    endif

    if(js .eq. jna2so4)then
       aer(ina_a,jliquid,ibin)  = aer(ina_a,jliquid,ibin) +   &
            2.*electrolyte(js,jsolid,ibin)
       aer(iso4_a,jliquid,ibin) = aer(iso4_a,jliquid,ibin) +   &
            electrolyte(js,jsolid,ibin)

       electrolyte(js,jsolid,ibin) = 0.0

       aer(ina_a,jp,ibin)  = electrolyte(jnano3,jp,ibin)  +   &
            electrolyte(jnacl,jp,ibin)   +   &
            2.*electrolyte(jna2so4,jp,ibin) +   &
            3.*electrolyte(jna3hso4,jp,ibin)+   &
            electrolyte(jnahso4,jp,ibin) +   &
            electrolyte(jnamsa,jp,ibin)

       aer(iso4_a,jp,ibin) = electrolyte(jcaso4,jp,ibin)  +   &
            electrolyte(jna2so4,jp,ibin) +   &
            2.*electrolyte(jna3hso4,jp,ibin)+   &
            electrolyte(jnahso4,jp,ibin) +   &
            electrolyte(jnh4so4,jp,ibin) +   &
            2.*electrolyte(jlvcite,jp,ibin) +   &
            electrolyte(jnh4hso4,jp,ibin)+   &
            electrolyte(jh2so4,jp,ibin)
       return
    endif

    if(js .eq. jna3hso4)then
       aer(ina_a,jliquid,ibin)  = aer(ina_a,jliquid,ibin) +   &
            3.*electrolyte(js,jsolid,ibin)
       aer(iso4_a,jliquid,ibin) = aer(iso4_a,jliquid,ibin) +   &
            2.*electrolyte(js,jsolid,ibin)

       electrolyte(js,jsolid,ibin) = 0.0

       aer(ina_a,jp,ibin)  = electrolyte(jnano3,jp,ibin)  +   &
            electrolyte(jnacl,jp,ibin)   +   &
            2.*electrolyte(jna2so4,jp,ibin) +   &
            3.*electrolyte(jna3hso4,jp,ibin)+   &
            electrolyte(jnahso4,jp,ibin) +   &
            electrolyte(jnamsa,jp,ibin)

       aer(iso4_a,jp,ibin) = electrolyte(jcaso4,jp,ibin)  +   &
            electrolyte(jna2so4,jp,ibin) +   &
            2.*electrolyte(jna3hso4,jp,ibin)+   &
            electrolyte(jnahso4,jp,ibin) +   &
            electrolyte(jnh4so4,jp,ibin) +   &
            2.*electrolyte(jlvcite,jp,ibin) +   &
            electrolyte(jnh4hso4,jp,ibin)+   &
            electrolyte(jh2so4,jp,ibin)
       return
    endif

    if(js .eq. jnahso4)then
       aer(ina_a,jliquid,ibin)  = aer(ina_a,jliquid,ibin) +   &
            electrolyte(js,jsolid,ibin)
       aer(iso4_a,jliquid,ibin) = aer(iso4_a,jliquid,ibin) +   &
            electrolyte(js,jsolid,ibin)

       electrolyte(js,jsolid,ibin) = 0.0

       aer(ina_a,jp,ibin)  = electrolyte(jnano3,jp,ibin)  +   &
            electrolyte(jnacl,jp,ibin)   +   &
            2.*electrolyte(jna2so4,jp,ibin) +   &
            3.*electrolyte(jna3hso4,jp,ibin)+   &
            electrolyte(jnahso4,jp,ibin) +   &
            electrolyte(jnamsa,jp,ibin)

       aer(iso4_a,jp,ibin) = electrolyte(jcaso4,jp,ibin)  +   &
            electrolyte(jna2so4,jp,ibin) +   &
            2.*electrolyte(jna3hso4,jp,ibin)+   &
            electrolyte(jnahso4,jp,ibin) +   &
            electrolyte(jnh4so4,jp,ibin) +   &
            2.*electrolyte(jlvcite,jp,ibin) +   &
            electrolyte(jnh4hso4,jp,ibin)+   &
            electrolyte(jh2so4,jp,ibin)
       return
    endif

    if(js .eq. jnh4no3)then
       aer(inh4_a,jliquid,ibin) = aer(inh4_a,jliquid,ibin) +   &
            electrolyte(js,jsolid,ibin)
       aer(ino3_a,jliquid,ibin) = aer(ino3_a,jliquid,ibin) +   &
            electrolyte(js,jsolid,ibin)

       electrolyte(js,jsolid,ibin) = 0.0

       aer(inh4_a,jp,ibin) = electrolyte(jnh4no3,jp,ibin) +   &
            electrolyte(jnh4cl,jp,ibin)  +   &
            2.*electrolyte(jnh4so4,jp,ibin) +   &
            3.*electrolyte(jlvcite,jp,ibin) +   &
            electrolyte(jnh4hso4,jp,ibin)+   &
            electrolyte(jnh4msa,jp,ibin)

       aer(ino3_a,jp,ibin) = electrolyte(jnano3,jp,ibin)  +   &
            2.*electrolyte(jcano3,jp,ibin)  +   &
            electrolyte(jnh4no3,jp,ibin) +   &
            electrolyte(jhno3,jp,ibin)
       return
    endif

    if(js .eq. jnh4cl)then
       aer(inh4_a,jliquid,ibin) = aer(inh4_a,jliquid,ibin) +   &
            electrolyte(js,jsolid,ibin)
       aer(icl_a,jliquid,ibin)  = aer(icl_a,jliquid,ibin) +   &
            electrolyte(js,jsolid,ibin)

       electrolyte(js,jsolid,ibin) = 0.0

       aer(inh4_a,jp,ibin) = electrolyte(jnh4no3,jp,ibin) +   &
            electrolyte(jnh4cl,jp,ibin)  +   &
            2.*electrolyte(jnh4so4,jp,ibin) +   &
            3.*electrolyte(jlvcite,jp,ibin) +   &
            electrolyte(jnh4hso4,jp,ibin)+   &
            electrolyte(jnh4msa,jp,ibin)

       aer(icl_a,jp,ibin)  = electrolyte(jnacl,jp,ibin)   +   &
            2.*electrolyte(jcacl2,jp,ibin)  +   &
            electrolyte(jnh4cl,jp,ibin)  +   &
            electrolyte(jhcl,jp,ibin)
       return
    endif

    if(js .eq. jnano3)then
       aer(ina_a,jliquid,ibin)  = aer(ina_a,jliquid,ibin) +   &
            electrolyte(js,jsolid,ibin)
       aer(ino3_a,jliquid,ibin) = aer(ino3_a,jliquid,ibin) +   &
            electrolyte(js,jsolid,ibin)

       electrolyte(js,jsolid,ibin) = 0.0

       aer(ina_a,jp,ibin)  = electrolyte(jnano3,jp,ibin)  +   &
            electrolyte(jnacl,jp,ibin)   +   &
            2.*electrolyte(jna2so4,jp,ibin) +   &
            3.*electrolyte(jna3hso4,jp,ibin)+   &
            electrolyte(jnahso4,jp,ibin) +   &
            electrolyte(jnamsa,jp,ibin)

       aer(ino3_a,jp,ibin) = electrolyte(jnano3,jp,ibin)  +   &
            2.*electrolyte(jcano3,jp,ibin)  +   &
            electrolyte(jnh4no3,jp,ibin) +   &
            electrolyte(jhno3,jp,ibin)
       return
    endif

    if(js .eq. jnacl)then
       aer(ina_a,jliquid,ibin)  = aer(ina_a,jliquid,ibin) +   &
            electrolyte(js,jsolid,ibin)
       aer(icl_a,jliquid,ibin)  = aer(icl_a,jliquid,ibin) +   &
            electrolyte(js,jsolid,ibin)

       electrolyte(js,jsolid,ibin) = 0.0

       aer(ina_a,jp,ibin)  = electrolyte(jnano3,jp,ibin)  +   &
            electrolyte(jnacl,jp,ibin)   +   &
            2.*electrolyte(jna2so4,jp,ibin) +   &
            3.*electrolyte(jna3hso4,jp,ibin)+   &
            electrolyte(jnahso4,jp,ibin) +   &
            electrolyte(jnamsa,jp,ibin)

       aer(icl_a,jp,ibin)  = electrolyte(jnacl,jp,ibin)   +   &
            2.*electrolyte(jcacl2,jp,ibin)  +   &
            electrolyte(jnh4cl,jp,ibin)  +   &
            electrolyte(jhcl,jp,ibin)
       return
    endif

    if(js .eq. jcano3)then
       aer(ica_a,jliquid,ibin)  = aer(ica_a,jliquid,ibin) +   &
            electrolyte(js,jsolid,ibin)
       aer(ino3_a,jliquid,ibin) = aer(ino3_a,jliquid,ibin) +   &
            2.*electrolyte(js,jsolid,ibin)

       electrolyte(js,jsolid,ibin) = 0.0

       aer(ica_a,jp,ibin)  = electrolyte(jcaso4,jp,ibin)  +   &
            electrolyte(jcano3,jp,ibin)  +   &
            electrolyte(jcacl2,jp,ibin)  +   &
            electrolyte(jcaco3,jp,ibin)  +   &
            electrolyte(jcamsa2,jp,ibin)

       aer(ino3_a,jp,ibin) = electrolyte(jnano3,jp,ibin)  +   &
            2.*electrolyte(jcano3,jp,ibin)  +   &
            electrolyte(jnh4no3,jp,ibin) +   &
            electrolyte(jhno3,jp,ibin)
       return
    endif

    if(js .eq. jcacl2)then
       aer(ica_a,jliquid,ibin) = aer(ica_a,jliquid,ibin) +   &
            electrolyte(js,jsolid,ibin)
       aer(icl_a,jliquid,ibin) = aer(icl_a,jliquid,ibin) +   &
            2.*electrolyte(js,jsolid,ibin)

       electrolyte(js,jsolid,ibin) = 0.0

       aer(ica_a,jp,ibin)  = electrolyte(jcaso4,jp,ibin)  +   &
            electrolyte(jcano3,jp,ibin)  +   &
            electrolyte(jcacl2,jp,ibin)  +   &
            electrolyte(jcaco3,jp,ibin)  +   &
            electrolyte(jcamsa2,jp,ibin)

       aer(icl_a,jp,ibin)  = electrolyte(jnacl,jp,ibin)   +   &
            2.*electrolyte(jcacl2,jp,ibin)  +   &
            electrolyte(jnh4cl,jp,ibin)  +   &
            electrolyte(jhcl,jp,ibin)
       return
    endif

    return
  end subroutine MESA_dissolve_small_salt

  !***********************************************************************
  ! electrolyte formation subroutines
  !
  ! author: Rahul A. Zaveri
  ! update: june 2000
  !-----------------------------------------------------------------------
  subroutine form_caso4(store,jp,ibin,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         iso4_a,ica_a,jcaso4

    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jcaso4,jp,ibin) = min(store(ica_a),store(iso4_a))
    store(ica_a)  = ( (store(ica_a)) -   &
         (electrolyte(jcaso4,jp,ibin)) )
    store(iso4_a) = ( (store(iso4_a)) -   &
         (electrolyte(jcaso4,jp,ibin)) )
    store(ica_a)  = max(0.d0, store(ica_a))
    store(iso4_a) = max(0.d0, store(iso4_a))

    return
  end subroutine form_caso4

  subroutine form_camsa2(store,jp,ibin,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         imsa_a,ica_a,jcamsa2
    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jcamsa2,jp,ibin) = min(store(ica_a),0.5*store(imsa_a))
    store(ica_a)  = ( (store(ica_a)) -   &
         (electrolyte(jcamsa2,jp,ibin)) )
    store(imsa_a) = ( (store(imsa_a)) -   &
         (2.*electrolyte(jcamsa2,jp,ibin)) )
    store(ica_a)  = max(0.d0, store(ica_a))
    store(imsa_a) = max(0.d0, store(imsa_a))

    return
  end subroutine form_camsa2

  !***********************************************************************
  ! computes mass transfer coefficients for each condensing species for
  ! all the aerosol bins
  !
  ! author: Rahul A. Zaveri
  ! update: jan 2005
  !-----------------------------------------------------------------------
  subroutine aerosolmtc( jaerosolstate, num_a, Dp_wet_a, sigmag_a, P_atm, T_K, kg )

    use module_data_mosaic_aero,  only: r8, nbin_a_max, nbin_a, naer, naercomp,             &!Parameters
         ngas_aerchtot, ngas_volatile, nelectrolyte, ngas_ioa,                              &
         mMODAL, no_aerosol, mUNSTRUCTURED, mSECTIONAL, mSIZE_FRAMEWORK,                    &!Input
         isoa_first, mw_gas, v_molar_gas,                                                   &!TBD
         i_gas2bin_uptk_flag, m_gas2bin_uptk_flag,                                          &
         use_cam5mam_accom_coefs, mosaic_vars_aa_type

    implicit none

    !Subroutine Arguments
    integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate

    real(r8), intent(in) :: P_atm,T_K
    real(r8), intent(in), dimension(nbin_a_max) :: num_a
    real(r8), intent(in), dimension(nbin_a_max) :: sigmag_a
    real(r8), intent(inout), dimension(nbin_a_max) :: Dp_wet_a
    real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg

    ! local variables
    integer nghq
    parameter (nghq = 2)         ! gauss-hermite quadrature order
    integer ibin, iq, iv
    real(r8) :: tworootpi, root2, beta
    parameter (tworootpi = 3.5449077, root2 = 1.4142135, beta = 2.0)
    real(r8) :: cdum, Dp, Dp_avg, Fkn, Kn, lnsg, lnDpgn, lnDp, speed,   &
                sumghq, tmpa
    real(r8) :: xghq(nghq), wghq(nghq)                           ! quadrature abscissae and weights
    real(r8) :: accom(ngas_aerchtot)                             ! keep local
    real(r8) :: freepath(ngas_aerchtot)                          ! keep local
    real(r8) :: Dg(ngas_aerchtot)                                ! keep local
    !real(r8) :: fuchs_sutugin                                   ! mosaic func
    !real(r8) :: gas_diffusivity                                 ! mosaic func
    !real(r8) :: mean_molecular_speed                            ! mosaic func

    ! mass accommodation coefficients
    tmpa = 0.1
    if ( use_cam5mam_accom_coefs > 0 ) tmpa = 0.65
    accom(1:ngas_aerchtot) = tmpa  ! default
!   accom(ih2so4_g)  = tmpa
!   accom(ihno3_g)   = tmpa
!   accom(ihcl_g)    = tmpa
!   accom(inh3_g)    = tmpa
!   accom(imsa_g)    = tmpa
!   accom(iaro1_g)   = tmpa
!   accom(iaro2_g)   = tmpa
!   accom(ialk1_g)   = tmpa
!   accom(iole1_g)   = tmpa
!   accom(iapi1_g)   = tmpa
!   accom(iapi2_g)   = tmpa
!   accom(ilim1_g)   = tmpa
!   accom(ilim2_g)   = tmpa

    ! quadrature weights
    xghq(1) =  0.70710678
    xghq(2) = -0.70710678
    wghq(1) =  0.88622693
    wghq(2) =  0.88622693

    ! calculate gas diffusivity and mean free path for condensing gases
    ! ioa
    do iv = 1, ngas_ioa
       speed  = mean_molecular_speed(T_K,mw_gas(iv))     ! cm/s
       Dg(iv) = gas_diffusivity(T_K,P_atm,mw_gas(iv),v_molar_gas(iv)) ! cm^2/s
       freepath(iv) = 3.*Dg(iv)/speed                    ! cm
    enddo

    ! soa
    do iv = isoa_first, ngas_volatile
       speed = mean_molecular_speed(T_K,mw_gas(iv))      ! cm/s
       Dg(iv) = 0.1                                      ! cm^2/s
       freepath(iv) = 3.*Dg(iv)/speed
    enddo

! het-rct gases
    do iv = (ngas_volatile+1), ngas_aerchtot
       speed = mean_molecular_speed(t_k,mw_gas(iv))    ! cm/s
       dg(iv) = gas_diffusivity(t_k,p_atm,mw_gas(iv),v_molar_gas(iv)) ! cm^2/s
       freepath(iv) = 3.*dg(iv)/speed                  ! cm
    enddo

    ! calc mass transfer coefficients for gases over various aerosol bins

    if (mSIZE_FRAMEWORK .eq. mMODAL) then

       ! for modal approach
       do 10 ibin = 1, nbin_a

          if(jaerosolstate(ibin) .eq. no_aerosol)goto 10

          lnsg   = log(sigmag_a(ibin))

          ! following 2 lines were incorrect as Dp_wet_a is wet "average" Dp
          !       Dpgn_a(ibin) = Dp_wet_a(ibin)  ! cm
          !       lnDpgn = log(Dpgn_a(ibin))
          ! do this instead which gives
          ! lnDpgn = ln( wet geometric-mean Dp of number distribution )
          lnDpgn = log(Dp_wet_a(ibin)) - 1.5*lnsg*lnsg

          cdum   = tworootpi*num_a(ibin)*   &
               exp(beta*lnDpgn + 0.5*(beta*lnsg)**2)

          do 20 iv = 1, ngas_aerchtot

             sumghq = 0.0_r8
             do 30 iq = 1, nghq  ! sum over gauss-hermite quadrature points
                lnDp = lnDpgn + beta*lnsg**2 + root2*lnsg*xghq(iq)
                Dp = exp(lnDp)
                Kn = 2.*freepath(iv)/Dp
                Fkn = fuchs_sutugin(Kn,accom(iv))
                sumghq = sumghq + wghq(iq)*Dp*Fkn/(Dp**beta)
30              continue

                kg(iv,ibin) = cdum*Dg(iv)*sumghq         ! 1/s

20              continue
10     continue

    elseif ((mSIZE_FRAMEWORK .eq. mSECTIONAL   ) .or. &
         (mSIZE_FRAMEWORK .eq. mUNSTRUCTURED)) then

       ! for sectional approach
       do 11 ibin = 1, nbin_a

          if(jaerosolstate(ibin) .eq. no_aerosol)goto 11

          cdum  = 6.283185*Dp_wet_a(ibin)*num_a(ibin)

          do 21 iv = 1, ngas_aerchtot
             Kn = 2.*freepath(iv)/Dp_wet_a(ibin)
             Fkn = fuchs_sutugin(Kn,accom(iv))
             kg(iv,ibin) = cdum*Dg(iv)*Fkn              ! 1/s
21           continue

11     continue

    else
       call wrf_message('Error in the choice of mSIZE_FRAMEWORK')
       call wrf_error_fatal('Stopping in subr. aerosolmtc')
    endif

    if (m_gas2bin_uptk_flag <= 0) then
       ! when m_gas2bin_uptk_flag <= 0, some gases cannot condense onto every bin
       do ibin = 1, nbin_a
          do iv = 1, ngas_aerchtot
             if (i_gas2bin_uptk_flag(iv,ibin) <= 0) kg(iv,ibin) = 0.0
          end do
       end do
    end if

    return
  end subroutine aerosolmtc

  !***********************************************************************
  ! calculates dry and wet aerosol properties: density, refractive indices
  !
  ! author: Rahul A. Zaveri
  ! update: jan 2005
  !-----------------------------------------------------------------------
  subroutine calc_dry_n_wet_aerosol_props(                          &
     ibin, jaerosolstate, aer, electrolyte, water_a, num_a,         &  ! input
     dens_comp_a, mw_comp_a, dens_aer_mac, mw_aer_mac, ref_index_a, &  ! input
     Dp_dry_a, Dp_wet_a, dp_core_a,                                 &  ! output
     area_dry_a, area_wet_a, mass_dry_a, mass_wet_a,                &  ! output
     vol_dry_a, vol_wet_a, dens_dry_a, dens_wet_a,                  &  ! output
     ri_shell_a, ri_core_a, ri_avg_a                                )  ! output
    !      include 'v33com9a'

    use module_data_mosaic_constants,  only: piover4,piover6,third
    use module_data_mosaic_aero,  only: r8,nbin_a_max,naer,nelectrolyte,naercomp,  &!Parameters
         ngas_soa,no_aerosol,msectional,                                           &!Parameters
         maeroptic_aero,msize_framework,                                           &!Input
         inh4_a,ina_a,ica_a,ico3_a,imsa_a,icl_a,ino3_a,jtotal,iso4_a,ioc_a,joc,    &!TBD
         ibc_a,jbc,ioin_a,join,jh2o, isoa_first, jsoa_first

    use module_data_mosaic_asecthp, only: dcen_sect,isize_of_ibin,itype_of_ibin

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    integer, intent(in), dimension(nbin_a_max) :: jaerosolstate

    real(r8), intent(in), dimension(nbin_a_max) :: num_a
    real(r8), intent(in), dimension(naer) :: mw_aer_mac,dens_aer_mac
    real(r8), intent(in), dimension(naercomp) :: dens_comp_a,mw_comp_a
    real(r8), intent(inout), dimension(nbin_a_max) :: Dp_dry_a,Dp_wet_a
    real(r8), intent(inout), dimension(nbin_a_max) :: dp_core_a,vol_dry_a
    real(r8), intent(inout), dimension(nbin_a_max) :: vol_wet_a,dens_wet_a,water_a
    real(r8), intent(inout), dimension(nbin_a_max) :: area_dry_a,area_wet_a
    real(r8), intent(inout), dimension(nbin_a_max) :: mass_dry_a,mass_wet_a
    real(r8), intent(inout), dimension(nbin_a_max) :: dens_dry_a
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    complex, intent(in), dimension(naercomp) :: ref_index_a
    complex, intent(inout), dimension(nbin_a_max) :: ri_avg_a,ri_core_a,ri_shell_a
    ! local variables
    integer i, iaer, isize, itype, j, jc, je, k
    real(r8) :: aer_H, duma, vol_core, vol_shell, vol_dum
    real(r8),dimension(naercomp) :: comp_a
    complex rixvol_tot, rixvol_core, rixvol_shell

    ! calculate dry mass and dry volume of a bin
    mass_dry_a(ibin) = 0.0                ! initialize to 0.0
    vol_dry_a(ibin)  = 0.0                ! initialize to 0.0
    area_dry_a(ibin) = 0.0                ! initialize to 0.0

    if(jaerosolstate(ibin) .ne. no_aerosol)then

       aer_H = (2.*aer(iso4_a,jtotal,ibin) +   &
            aer(ino3_a,jtotal,ibin) +   &
            aer(icl_a,jtotal,ibin)  +   &
            aer(imsa_a,jtotal,ibin) +   &
            2.*aer(ico3_a,jtotal,ibin))-   &
            (2.*aer(ica_a,jtotal,ibin)  +   &
            aer(ina_a,jtotal,ibin)  +   &
            aer(inh4_a,jtotal,ibin))
       aer_H = max(aer_H, 0.0d0)

       do iaer = 1, naer
          mass_dry_a(ibin) = mass_dry_a(ibin) +   &
               aer(iaer,jtotal,ibin)*mw_aer_mac(iaer)     ! ng/m^3(air)
          vol_dry_a(ibin) = vol_dry_a(ibin) +   &
               aer(iaer,jtotal,ibin)*mw_aer_mac(iaer)/dens_aer_mac(iaer)          ! ncc/m^3(air)
       enddo
       mass_dry_a(ibin) = mass_dry_a(ibin) + aer_H
       vol_dry_a(ibin) = vol_dry_a(ibin) + aer_H

       mass_dry_a(ibin) = mass_dry_a(ibin)*1.e-15                 ! g/cc(air)
       vol_dry_a(ibin) = vol_dry_a(ibin)*1.e-15                   ! cc(aer)/cc(air)

       ! wet mass and wet volume
       mass_wet_a(ibin) = mass_dry_a(ibin) + water_a(ibin)*1.e-3  ! g/cc(air)
       vol_wet_a(ibin)  = vol_dry_a(ibin) + water_a(ibin)*1.e-3   ! cc(aer)/cc(air)

       ! calculate mean dry and wet particle densities
       dens_dry_a(ibin) = mass_dry_a(ibin)/vol_dry_a(ibin)                ! g/cc(aerosol)
       dens_wet_a(ibin) = mass_wet_a(ibin)/vol_wet_a(ibin)                ! g/cc(aerosol)

       ! calculate mean dry and wet particle diameters
       Dp_dry_a(ibin)=(vol_dry_a(ibin)/(piover6*num_a(ibin)))**third      ! cm
       Dp_wet_a(ibin)=(vol_wet_a(ibin)/(piover6*num_a(ibin)))**third      ! cm

       ! calculate mean dry and wet particle surface areas
       area_dry_a(ibin)= piover4*num_a(ibin)*Dp_dry_a(ibin)**2    ! cm^2/cc(air)
       area_wet_a(ibin)= piover4*num_a(ibin)*Dp_wet_a(ibin)**2    ! cm^2/cc(air)

       ! calculate volume average refractive index
       !   load comp_a array with component mass concentrations

       ! rahul had turned this off, but it is needed
       !        if(1 == 1)go to 100               ! TEMP
       if (maeroptic_aero <= 0) goto 100

       do je = 1, nelectrolyte
          comp_a(je)=electrolyte(je,jtotal,ibin)*mw_comp_a(je)*1.e-15     ! g/cc(air)
       enddo
       comp_a(joc)  = aer(ioc_a,  jtotal,ibin)*mw_comp_a(joc  )*1.e-15    ! g/cc(air)
       comp_a(jbc)  = aer(ibc_a,  jtotal,ibin)*mw_comp_a(jbc  )*1.e-15    ! g/cc(air)
       comp_a(join) = aer(ioin_a, jtotal,ibin)*mw_comp_a(join )*1.e-15    ! g/cc(air)

!      comp_a(jaro1)= aer(iaro1_a,jtotal,ibin)*mw_comp_a(jaro1)*1.e-15    ! g/cc(air)
!      comp_a(jaro2)= aer(iaro2_a,jtotal,ibin)*mw_comp_a(jaro2)*1.e-15    ! g/cc(air)
!      comp_a(jalk1)= aer(ialk1_a,jtotal,ibin)*mw_comp_a(jalk1)*1.e-15    ! g/cc(air)
!      comp_a(jole1)= aer(iole1_a,jtotal,ibin)*mw_comp_a(jole1)*1.e-15    ! g/cc(air)
!      comp_a(japi1)= aer(iapi1_a,jtotal,ibin)*mw_comp_a(japi1)*1.e-15    ! g/cc(air)
!      comp_a(japi2)= aer(iapi2_a,jtotal,ibin)*mw_comp_a(japi2)*1.e-15    ! g/cc(air)
!      comp_a(jlim1)= aer(ilim1_a,jtotal,ibin)*mw_comp_a(jlim1)*1.e-15    ! g/cc(air)
!      comp_a(jlim2)= aer(ilim2_a,jtotal,ibin)*mw_comp_a(jlim2)*1.e-15    ! g/cc(air)
       do k = 1, ngas_soa
       j = jsoa_first + k - 1
       i = isoa_first + k - 1
       comp_a(j) = aer(i,jtotal,ibin)*mw_comp_a(j)*1.e-15    ! g/cc(air)
       end do

       comp_a(jh2o) = water_a(ibin)*1.e-3                         ! g/cc(air)

       rixvol_tot   = (0.0,0.0)
       do jc = 1, naercomp
          comp_a(jc) = max( 0.0d0, comp_a(jc) )
          rixvol_tot = rixvol_tot   &
               + ref_index_a(jc)*comp_a(jc)/dens_comp_a(jc)
       enddo
       ri_avg_a(ibin) = rixvol_tot/vol_wet_a(ibin)

       !
       ! shell/core calcs - first set values to default (corresponding to zero core)
       !
       ri_shell_a(ibin) = ri_avg_a(ibin)
       ri_core_a(ibin)  = (0.0,0.0)
       Dp_core_a(ibin)  = 0.0

       ! sum ri*vol and vol for core species (bc and optionally oin=dust)
       ! currently just bc in core, but what about insoluble oin and dust species ???
       jc = jbc
       rixvol_core  = ref_index_a(jc)*comp_a(jc)/dens_comp_a(jc)
       vol_core = comp_a(jc)/dens_comp_a(jc)
       vol_core = max( 0.0d0, min( vol_core, vol_wet_a(ibin) ) )

       ! neglect core if (core volume) < 1.0d-9*(total volume)
       !              or (core volume) < 1.0d-22 cm3 = (0.58 nm)**3
       ! neglect shell using similar criteria
       vol_dum = max( 1.0d-22, 1.0d-9*vol_wet_a(ibin) )
       vol_shell = vol_wet_a(ibin) - vol_core
       if (vol_core >= vol_dum) then
          if (vol_shell < vol_dum) then
             ri_shell_a(ibin)  = (0.0,0.0)
             ri_core_a(ibin) = ri_avg_a(ibin)
             Dp_core_a(ibin) = Dp_wet_a(ibin)
          else
             ri_core_a(ibin) = rixvol_core/vol_core
             Dp_core_a(ibin) = Dp_wet_a(ibin)   &
                  * (vol_core/vol_wet_a(ibin))**third

             if (vol_shell >= vol_dum) then
                rixvol_shell = rixvol_tot - rixvol_core
                ri_shell_a(ibin) = rixvol_shell/vol_shell
             else
                ri_shell_a(ibin) = (0.0,0.0)
             endif
          endif
       endif

    else
       ! use defaults when (jaerosolstate(ibin) .eq. no_aerosol)

       dens_dry_a(ibin) = 1.0      ! g/cc(aerosol)
       dens_wet_a(ibin) = 1.0      ! g/cc(aerosol)
       !        Dp_dry_a(ibin) = dcen_sect(ibin)  ! cm
       !        Dp_wet_a(ibin) = dcen_sect(ibin)  ! cm
       if (msize_framework == msectional) then
          isize = isize_of_ibin(ibin)
          itype = itype_of_ibin(ibin)
          Dp_dry_a(ibin) = dcen_sect(isize,itype)
          Dp_wet_a(ibin) = Dp_dry_a(ibin)
       end if

       ri_avg_a(ibin) = (1.5,0.0)
       ri_shell_a(ibin) = (1.5,0.0)
       ri_core_a(ibin)  = (0.0,0.0)
       Dp_core_a(ibin)  = 0.0

    endif   ! if(jaerosolstate(ibin) .ne. no_aerosol)then

100 continue

    return
  end subroutine calc_dry_n_wet_aerosol_props

  !***********************************************************************
  ! forms electrolytes from ions
  !
  ! author: Rahul A. Zaveri
  ! update: june 2000
  !-----------------------------------------------------------------------
  subroutine form_electrolytes(jp,ibin,XT,aer,gas,electrolyte,total_species,tot_cl_in)
    use module_data_mosaic_aero, only: r8,naer,nbin_a_max,           &
         ngas_aerchtot,ngas_volatile,nelectrolyte,jsolid,            &
         imsa_a,iso4_a,ica_a,ina_a,inh4_a,ino3_a,icl_a,ico3_a

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin, jp
    real(r8), intent(inout) :: XT
    real(r8), intent(inout) :: tot_cl_in
    real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
    real(r8), intent(inout), dimension(ngas_volatile) :: total_species
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    ! local variables
    integer i, iXT_case, j, je
    real(r8) :: sum_dum, XNa_prime, XNH4_prime, XT_prime
    real(r8) :: store(naer)

    ! remove negative concentrations, if any
    !      do i=1,naer
    !        aer(i,jp,ibin) = max(0.0d0, aer(i,jp,ibin))  ! EFFI
    !      enddo

    !      call calculate_XT(ibin,jp,XT,aer)      ! EFFI

    if( (aer(iso4_a,jp,ibin)+aer(imsa_a,jp,ibin)) .gt.0.0)then
       XT   = ( aer(inh4_a,jp,ibin) +   &
            aer(ina_a,jp,ibin)  +   &
            2.*aer(ica_a,jp,ibin) )/   &
            (aer(iso4_a,jp,ibin)+0.5*aer(imsa_a,jp,ibin))
    else
       XT   = -1.0
    endif

!   if(XT .ge. 1.9999 .or. XT.lt.0.)then
    if(XT .ge. 2.0 .or. XT.lt.0.)then         ! RAZ 11/10/2014
       iXT_case = 1       ! near neutral (acidity is caused by HCl and/or HNO3)
    else
       iXT_case = 2       ! acidic (acidity is caused by excess SO4)
    endif

    ! initialize
    !
    ! put total aer(*) into store(*)
    store(iso4_a) = aer(iso4_a,jp,ibin)
    store(ino3_a) = aer(ino3_a,jp,ibin)
    store(icl_a)  = aer(icl_a, jp,ibin)
    store(imsa_a) = aer(imsa_a,jp,ibin)
    store(ico3_a) = aer(ico3_a,jp,ibin)
    store(inh4_a) = aer(inh4_a,jp,ibin)
    store(ina_a)  = aer(ina_a, jp,ibin)
    store(ica_a)  = aer(ica_a, jp,ibin)

    do j=1,nelectrolyte
       electrolyte(j,jp,ibin) = 0.0
    enddo

    !
    !---------------------------------------------------------
    !
    if(iXT_case.eq.1)then

       ! XT >= 2   : sulfate deficient
       call form_caso4(store,jp,ibin,electrolyte)
       call form_camsa2(store,jp,ibin,electrolyte)
       call form_na2so4(store,jp,ibin,electrolyte)
       call form_namsa(store,jp,ibin,electrolyte)
       call form_cano3(store,jp,ibin,electrolyte)
       call form_nano3(store,jp,ibin,electrolyte)
       call form_nacl(store,jp,ibin,aer,gas,electrolyte,total_species,tot_cl_in)
       call form_cacl2(store,jp,ibin,electrolyte)
       call form_caco3(store,jp,ibin,aer,electrolyte)
       call form_nh4so4(store,jp,ibin,electrolyte)
       call form_nh4msa(store,jp,ibin,electrolyte)
       call form_nh4no3(store,jp,ibin,electrolyte)
       call form_nh4cl(store,jp,ibin,electrolyte)
       call form_msa(store,jp,ibin,electrolyte)

       if(jp .eq. jsolid)then
          call degas_hno3(store,jp,ibin,aer,gas,electrolyte)
          call degas_hcl(store,jp,ibin,aer,gas,electrolyte)
          call degas_nh3(store,jp,ibin,aer,gas)
       else
          call form_hno3(store,jp,ibin,electrolyte)
          call form_hcl(store,jp,ibin,electrolyte)
          call degas_nh3(store,jp,ibin,aer,gas)
       endif

    elseif(iXT_case.eq.2)then

       ! XT < 2   : sulfate enough or sulfate excess

       call form_caso4(store,jp,ibin,electrolyte)
       call form_camsa2(store,jp,ibin,electrolyte)
       call form_namsa(store,jp,ibin,electrolyte)
       call form_nh4msa(store,jp,ibin,electrolyte)
       call form_msa(store,jp,ibin,electrolyte)

       if(store(iso4_a).eq.0.0)goto 10

       XT_prime =(store(ina_a)+store(inh4_a))/   &
            store(iso4_a)
       XNa_prime=0.5*store(ina_a)/store(iso4_a) + 1.

       if(XT_prime.ge.XNa_prime)then
          call form_na2so4(store,jp,ibin,electrolyte)
          XNH4_prime = 0.0
          if(store(iso4_a).gt.1.e-15)then
             XNH4_prime = store(inh4_a)/store(iso4_a)
          endif

          if(XNH4_prime .ge. 1.5)then
             call form_nh4so4_lvcite(store,jp,ibin,electrolyte)
          else
             call form_lvcite_nh4hso4(store,jp,ibin,electrolyte)
          endif

       elseif(XT_prime.ge.1.)then
          call form_nh4hso4(store,jp,ibin,electrolyte)
          call form_na2so4_nahso4(store,jp,ibin,electrolyte)
       elseif(XT_prime.lt.1.)then
          call form_nahso4(store,jp,ibin,electrolyte)
          call form_nh4hso4(store,jp,ibin,electrolyte)
          call form_h2so4(store,jp,ibin,electrolyte)
       endif

10     if(jp .eq. jsolid)then
          call degas_hno3(store,jp,ibin,aer,gas,electrolyte)
          call degas_hcl(store,jp,ibin,aer,gas,electrolyte)
          call degas_nh3(store,jp,ibin,aer,gas)
       else
          call form_hno3(store,jp,ibin,electrolyte)
          call form_hcl(store,jp,ibin,electrolyte)
          call degas_nh3(store,jp,ibin,aer,gas)
       endif

    endif ! case 1, 2

    ! re-calculate ions to eliminate round-off errors
    call electrolytes_to_ions(jp, ibin,aer,electrolyte)
    !---------------------------------------------------------
    !
    ! calculate % composition EFFI
    !!      sum_dum = 0.0
    !!      do je = 1, nelectrolyte
    !!        electrolyte(je,jp,ibin) = max(0.d0,electrolyte(je,jp,ibin)) ! remove -ve  EFFI
    !!        sum_dum = sum_dum + electrolyte(je,jp,ibin)
    !!      enddo
    !!
    !!      if(sum_dum .eq. 0.)sum_dum = 1.0
    !!      electrolyte_sum(jp,ibin) = sum_dum
    !!
    !!      do je = 1, nelectrolyte
    !!        epercent(je,jp,ibin) = 100.*electrolyte(je,jp,ibin)/sum_dum
    !!      enddo

    return
  end subroutine form_electrolytes

   subroutine form_na2so4(store,jp,ibin,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         iso4_a,ina_a,jna2so4
    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store(naer)
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    electrolyte(jna2so4,jp,ibin) = min(.5*store(ina_a),   &
         store(iso4_a))
    store(ina_a) =( (store(ina_a)) -   &
         (2.*electrolyte(jna2so4,jp,ibin)) )
    store(iso4_a)=( (store(iso4_a)) -   &
         (electrolyte(jna2so4,jp,ibin)) )
    store(ina_a) =max(0.d0, store(ina_a))
    store(iso4_a)=max(0.d0, store(iso4_a))

    return
  end subroutine form_na2so4

  subroutine form_nahso4(store,jp,ibin,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         iso4_a,ina_a,jnahso4

    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jnahso4,jp,ibin) = min(store(ina_a),   &
         store(iso4_a))
    store(ina_a)  = ( (store(ina_a)) -   &
         (electrolyte(jnahso4,jp,ibin)) )
    store(iso4_a) = ( (store(iso4_a)) -   &
         (electrolyte(jnahso4,jp,ibin)) )
    store(ina_a)  = max(0.d0, store(ina_a))
    store(iso4_a) = max(0.d0, store(iso4_a))

    return
  end subroutine form_nahso4

  subroutine form_namsa(store,jp,ibin,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         imsa_a,ina_a,jnamsa
    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer)  :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jnamsa,jp,ibin) = min(store(ina_a),   &
         store(imsa_a))
    store(ina_a)  = ( (store(ina_a)) -   &
         (electrolyte(jnamsa,jp,ibin)) )
    store(imsa_a) = ( (store(imsa_a)) -   &
         (electrolyte(jnamsa,jp,ibin)) )
    store(ina_a)  = max(0.d0, store(ina_a))
    store(imsa_a) = max(0.d0, store(imsa_a))

    return
  end subroutine form_namsa

  subroutine form_nano3(store,jp,ibin,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ino3_a,ina_a,jnano3
    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jnano3,jp,ibin)=min(store(ina_a),store(ino3_a))
    store(ina_a)  = ( (store(ina_a)) -   &
         (electrolyte(jnano3,jp,ibin)) )
    store(ino3_a) = ( (store(ino3_a)) -   &
         (electrolyte(jnano3,jp,ibin)) )
    store(ina_a)  = max(0.d0, store(ina_a))
    store(ino3_a) = max(0.d0, store(ino3_a))

    return
  end subroutine form_nano3

  subroutine form_cano3(store,jp,ibin,electrolyte)        ! Ca(NO3)2
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ino3_a,ica_a,jcano3
    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jcano3,jp,ibin) = min(store(ica_a),0.5*store(ino3_a))

    store(ica_a)  = ( (store(ica_a)) -   &
         (electrolyte(jcano3,jp,ibin)) )
    store(ino3_a) = ( (store(ino3_a)) -   &
         (2.*electrolyte(jcano3,jp,ibin)) )
    store(ica_a)  = max(0.d0, store(ica_a))
    store(ino3_a) = max(0.d0, store(ino3_a))

    return
  end subroutine form_cano3

  subroutine form_cacl2(store,jp,ibin,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         icl_a,ica_a,jcacl2

    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jcacl2,jp,ibin) = min(store(ica_a),0.5*store(icl_a))

    store(ica_a)  = ( (store(ica_a)) -   &
         (electrolyte(jcacl2,jp,ibin)) )
    store(icl_a)  = ( (store(icl_a)) -   &
         (2.*electrolyte(jcacl2,jp,ibin)) )
    store(ica_a)  = max(0.d0, store(ica_a))
    store(icl_a)  = max(0.d0, store(icl_a))

    return
  end subroutine form_cacl2

  subroutine form_caco3(store,jp,ibin,aer,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,jsolid,     &
         jtotal,                                                                   &
         ica_a,jcaco3,ico3_a

    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    if(jp.eq.jtotal .or. jp.eq.jsolid)then
       electrolyte(jcaco3,jp,ibin) = store(ica_a)

       aer(ico3_a,jp,ibin)= electrolyte(jcaco3,jp,ibin)   ! force co3 = caco3

       store(ica_a) = 0.0
       store(ico3_a)= 0.0
    endif

    return
  end subroutine form_caco3

  subroutine form_nacl(store,jp,ibin,aer,gas,electrolyte,total_species,tot_cl_in)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ngas_aerchtot,ngas_volatile,jtotal,jsolid,jliquid,                        &
         ina_a,jnacl,icl_a,ihcl_g

    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
    real(r8), intent(inout), dimension(ngas_volatile) :: total_species
    real(r8), intent(inout) :: tot_cl_in
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    ! local

    electrolyte(jnacl,jp,ibin) = store(ina_a)

    store(ina_a) = 0.0
    store(icl_a) = ( (store(icl_a)) -   &
         (electrolyte(jnacl,jp,ibin)) )

    if(store(icl_a) .lt. 0.)then                          ! cl deficit in aerosol. take some from gas
       aer(icl_a,jp,ibin)= aer(icl_a,jp,ibin)- store(icl_a)       ! update aer(icl_a)

       if(jp .ne. jtotal)then
          aer(icl_a,jtotal,ibin)= aer(icl_a,jliquid,ibin)+ &      ! update for jtotal
               aer(icl_a,jsolid,ibin)
       endif

       gas(ihcl_g) = gas(ihcl_g) + store(icl_a)                   ! update gas(ihcl_g)

       if(gas(ihcl_g) .lt. 0.0)then
          total_species(ihcl_g) = total_species(ihcl_g) - gas(ihcl_g)     ! update total_species
          tot_cl_in = tot_cl_in - gas(ihcl_g)                             ! update tot_cl_in
       endif

       gas(ihcl_g) = max(0.d0, gas(ihcl_g))                               ! restrict gas(ihcl_g) to >= 0.
       store(icl_a) = 0.                                          ! force store(icl_a) to 0.

    endif

    store(icl_a) = max(0.d0, store(icl_a))

    return
  end subroutine form_nacl

  subroutine form_nh4so4(store,jp,ibin,electrolyte)       ! (nh4)2so4
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         iso4_a,inh4_a,jnh4so4
    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jnh4so4,jp,ibin)= min(.5*store(inh4_a),   &
         store(iso4_a))
    store(inh4_a)= ( (store(inh4_a)) -   &
         (2.*electrolyte(jnh4so4,jp,ibin)) )
    store(iso4_a)= ( (store(iso4_a)) -   &
         (electrolyte(jnh4so4,jp,ibin)) )
    store(inh4_a) = max(0.d0, store(inh4_a))
    store(iso4_a) = max(0.d0, store(iso4_a))

    return
  end subroutine form_nh4so4

  subroutine form_nh4hso4(store,jp,ibin,electrolyte)      ! nh4hso4
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         iso4_a,inh4_a,jnh4hso4
    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jnh4hso4,jp,ibin) = min(store(inh4_a),   &
         store(iso4_a))
    store(inh4_a)= ( (store(inh4_a)) -   &
         (electrolyte(jnh4hso4,jp,ibin)) )
    store(iso4_a)= ( (store(iso4_a)) -   &
         (electrolyte(jnh4hso4,jp,ibin)) )
    store(inh4_a) = max(0.d0, store(inh4_a))
    store(iso4_a) = max(0.d0, store(iso4_a))

    return
  end subroutine form_nh4hso4

  subroutine form_nh4msa(store,jp,ibin,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         imsa_a,inh4_a,jnh4msa
    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jnh4msa,jp,ibin) = min(store(inh4_a),   &
         store(imsa_a))
    store(inh4_a) = ( (store(inh4_a)) -   &
         (electrolyte(jnh4msa,jp,ibin)) )
    store(imsa_a) = ( (store(imsa_a)) -   &
         (electrolyte(jnh4msa,jp,ibin)) )
    store(inh4_a) = max(0.d0, store(inh4_a))
    store(imsa_a) = max(0.d0, store(imsa_a))

    return
  end subroutine form_nh4msa

  subroutine form_nh4cl(store,jp,ibin,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         icl_a,inh4_a,jnh4cl
    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jnh4cl,jp,ibin) = min(store(inh4_a),   &
         store(icl_a))
    store(inh4_a) = ( (store(inh4_a)) -   &
         (electrolyte(jnh4cl,jp,ibin)) )
    store(icl_a)  = ( (store(icl_a)) -   &
         (electrolyte(jnh4cl,jp,ibin)) )
    store(inh4_a) = max(0.d0, store(inh4_a))
    store(icl_a)  = max(0.d0, store(icl_a))

    return
  end subroutine form_nh4cl

  subroutine form_nh4no3(store,jp,ibin,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ino3_a,inh4_a,jnh4no3
    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jnh4no3,jp,ibin) = min(store(inh4_a),   &
         store(ino3_a))
    store(inh4_a) = ( (store(inh4_a)) -   &
         (electrolyte(jnh4no3,jp,ibin)) )
    store(ino3_a) = ( (store(ino3_a)) -   &
         (electrolyte(jnh4no3,jp,ibin)) )
    store(inh4_a) = max(0.d0, store(inh4_a))
    store(ino3_a) = max(0.d0, store(ino3_a))

    return
  end subroutine form_nh4no3

  subroutine form_nh4so4_lvcite(store,jp,ibin,electrolyte) ! (nh4)2so4 + (nh4)3h(so4)2
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         iso4_a,inh4_a,jnh4so4,jlvcite
    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jnh4so4,jp,ibin)= ( (2.*store(inh4_a)) -   &
         (3.*store(iso4_a)) )
    electrolyte(jlvcite,jp,ibin)= ( (2.*store(iso4_a)) -   &
         (store(inh4_a)) )
    electrolyte(jnh4so4,jp,ibin)= max(0.d0,   &
         electrolyte(jnh4so4,jp,ibin))
    electrolyte(jlvcite,jp,ibin)= max(0.d0,   &
         electrolyte(jlvcite,jp,ibin))
    store(inh4_a) = 0.
    store(iso4_a) = 0.

    return
  end subroutine form_nh4so4_lvcite

  subroutine form_lvcite_nh4hso4(store,jp,ibin,electrolyte) ! (nh4)3h(so4)2 + nh4hso4
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         iso4_a,inh4_a,jlvcite,jnh4hso4
    implicit none

    ! subr arguments
    integer, intent(in) ::  jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jlvcite,jp,ibin) = ( (store(inh4_a)) -   &
         (store(iso4_a)) )
    electrolyte(jnh4hso4,jp,ibin)= ( (3.*store(iso4_a)) -   &
         (2.*store(inh4_a)) )
    electrolyte(jlvcite,jp,ibin) = max(0.d0,   &
         electrolyte(jlvcite,jp,ibin))
    electrolyte(jnh4hso4,jp,ibin)= max(0.d0,   &
         electrolyte(jnh4hso4,jp,ibin))
    store(inh4_a) = 0.
    store(iso4_a) = 0.

    return
  end subroutine form_lvcite_nh4hso4

  subroutine form_na2so4_nahso4(store,jp,ibin,electrolyte) ! na2so4 + nahso4
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         iso4_a,ina_a,jna2so4,jnahso4
    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jna2so4,jp,ibin)= ( (store(ina_a)) -   &
         (store(iso4_a)) )
    electrolyte(jnahso4,jp,ibin)= ( (2.*store(iso4_a))-   &
         (store(ina_a)) )
    electrolyte(jna2so4,jp,ibin)= max(0.d0,   &
         electrolyte(jna2so4,jp,ibin))
    electrolyte(jnahso4,jp,ibin)= max(0.d0,   &
         electrolyte(jnahso4,jp,ibin))
    store(ina_a)  = 0.
    store(iso4_a) = 0.

    !     write(6,*)'na2so4 + nahso4'

    return
  end subroutine form_na2so4_nahso4

  subroutine form_h2so4(store,jp,ibin,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         iso4_a,jh2so4
    implicit none

    ! subr arguments
    integer, intent(in) ::  jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jh2so4,jp,ibin) = max(0.0d0, store(iso4_a))
    store(iso4_a) = 0.0

    return
  end subroutine form_h2so4

  subroutine form_msa(store,jp,ibin,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         imsa_a,jmsa
    implicit none

    ! subr arguments
    integer, intent(in) ::  jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jmsa,jp,ibin) = max(0.0d0, store(imsa_a))
    store(imsa_a) = 0.0

    return
  end subroutine form_msa

  subroutine form_hno3(store,jp,ibin,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ino3_a,jhno3
    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jhno3,jp,ibin) = max(0.0d0, store(ino3_a))
    store(ino3_a) = 0.0

    return
  end subroutine form_hno3

  subroutine form_hcl(store,jp,ibin,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         icl_a,jhcl
    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    electrolyte(jhcl,jp,ibin) = max(0.0d0, store(icl_a))
    store(icl_a) = 0.0

    return
  end subroutine form_hcl

 subroutine degas_hno3(store,jp,ibin,aer,gas,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ngas_aerchtot,jtotal,jliquid,jsolid,                                      &
         ino3_a,ihno3_g,jhno3

    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout), dimension(naer) :: store
    real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    store(ino3_a) = max(0.0d0, store(ino3_a))
    gas(ihno3_g) = gas(ihno3_g) + store(ino3_a)
    aer(ino3_a,jp,ibin) = ( (aer(ino3_a,jp,ibin)) -   &
         (store(ino3_a)) )
    aer(ino3_a,jp,ibin) = max(0.0d0,aer(ino3_a,jp,ibin))

    ! also do it for jtotal
    if(jp .ne. jtotal)then
       aer(ino3_a,jtotal,ibin) = aer(ino3_a,jsolid, ibin) +   &
            aer(ino3_a,jliquid,ibin)
    endif

    electrolyte(jhno3,jp,ibin) = 0.0
    store(ino3_a) = 0.0

    return
  end subroutine degas_hno3

  subroutine degas_hcl(store,jp,ibin,aer,gas,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ngas_aerchtot,jtotal,jliquid,jsolid,                                      &
         icl_a,ihcl_g,jhcl
    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout) :: store(naer)
    real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte

    store(icl_a) = max(0.0d0, store(icl_a))
    gas(ihcl_g) = gas(ihcl_g) + store(icl_a)
    aer(icl_a,jp,ibin) = ( (aer(icl_a,jp,ibin)) -   &
         (store(icl_a)) )
    aer(icl_a,jp,ibin) = max(0.0d0,aer(icl_a,jp,ibin))

    ! also do it for jtotal
    if(jp .ne. jtotal)then
       aer(icl_a,jtotal,ibin) = aer(icl_a,jsolid, ibin) +   &
            aer(icl_a,jliquid,ibin)
    endif

    electrolyte(jhcl,jp,ibin) = 0.0
    store(icl_a) = 0.0

    return
  end subroutine degas_hcl

  subroutine degas_nh3(store,jp,ibin,aer,gas)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ngas_aerchtot,jtotal,jliquid,jsolid,                                      &
         inh3_g,inh4_a

    implicit none

    ! subr arguments
    integer, intent(in) :: jp, ibin
    real(r8), intent(inout) :: store(naer)
    real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer

    store(inh4_a) = max(0.0d0, store(inh4_a))
    gas(inh3_g) = gas(inh3_g) + store(inh4_a)
    aer(inh4_a,jp,ibin) = ( (aer(inh4_a,jp,ibin)) -   &
         (store(inh4_a)) )
    aer(inh4_a,jp,ibin) = max(0.0d0,aer(inh4_a,jp,ibin))

    ! also do it for jtotal
    if(jp .ne. jtotal)then
       aer(inh4_a,jtotal,ibin)= aer(inh4_a,jsolid, ibin) +   &
            aer(inh4_a,jliquid,ibin)
    endif

    store(inh4_a) = 0.0

    return
  end subroutine degas_nh3

  !***********************************************************************
  ! subroutines to absorb and degas small amounts of volatile species
  !
  ! author: Rahul A. Zaveri
  ! update: jun 2002
  !-----------------------------------------------------------------------
  !
  ! nh4no3 (liquid)
  subroutine absorb_tiny_nh4no3(ibin,aer,gas,electrolyte,delta_nh3_max,            &
       delta_hno3_max,electrolyte_sum)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ngas_aerchtot,jtotal,jliquid,jsolid,                                      &
         inh4_a,ino3_a,inh3_g,ihno3_g
    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    real(r8), intent(inout), dimension(nbin_a_max) :: delta_nh3_max,delta_hno3_max
    real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
    real(r8), intent(inout), dimension(3,nbin_a_max) :: electrolyte_sum
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    ! local variables
    integer je
    real(r8) :: small_aer, small_gas, small_amt

    !! EFFI
    electrolyte_sum(jtotal,ibin) = 0.0
    do je = 1, nelectrolyte
       electrolyte_sum(jtotal,ibin) = electrolyte_sum(jtotal,ibin) + &
            electrolyte(je,jtotal,ibin)
    enddo
    !! EFFI

    small_gas = 0.01 * min(delta_nh3_max(ibin),delta_hno3_max(ibin))
    small_aer = 0.01 * electrolyte_sum(jtotal,ibin)
    if(small_aer .eq. 0.0)small_aer = small_gas

    small_amt = min(small_gas, small_aer)

    aer(inh4_a,jliquid,ibin) = aer(inh4_a,jliquid,ibin) + small_amt
    aer(ino3_a,jliquid,ibin) = aer(ino3_a,jliquid,ibin) + small_amt

    ! update jtotal
    aer(inh4_a,jtotal,ibin)  = aer(inh4_a,jsolid,ibin) +   &
         aer(inh4_a,jliquid,ibin)
    aer(ino3_a,jtotal,ibin)  = aer(ino3_a,jsolid,ibin) +   &
         aer(ino3_a,jliquid,ibin)

    ! update gas
    gas(inh3_g)  = ((gas(inh3_g)) - (small_amt))
    gas(ihno3_g) = ((gas(ihno3_g)) - (small_amt))

    return
  end subroutine absorb_tiny_nh4no3

  !--------------------------------------------------------------------
  ! nh4cl (liquid)
  subroutine absorb_tiny_nh4cl(ibin,aer,gas,electrolyte,delta_nh3_max,             &
       delta_hcl_max,electrolyte_sum)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ngas_aerchtot,jtotal,jliquid,jsolid,                                      &
         inh4_a,icl_a,inh3_g,ihcl_g
    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    real(r8), intent(inout), dimension(nbin_a_max) :: delta_nh3_max,delta_hcl_max
    real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
    real(r8), intent(inout), dimension(3,nbin_a_max) :: electrolyte_sum
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    ! local variables
    integer je
    real(r8) :: small_aer, small_gas, small_amt

    !! EFFI
    electrolyte_sum(jtotal,ibin) = 0.0
    do je = 1, nelectrolyte
       electrolyte_sum(jtotal,ibin) = electrolyte_sum(jtotal,ibin) + &
            electrolyte(je,jtotal,ibin)
    enddo
    !! EFFI

    small_gas = 0.01 * min(delta_nh3_max(ibin), delta_hcl_max(ibin))
    small_aer = 0.01 * electrolyte_sum(jtotal,ibin)
    if(small_aer .eq. 0.0)small_aer = small_gas

    small_amt = min(small_gas, small_aer)

    aer(inh4_a,jliquid,ibin) = aer(inh4_a,jliquid,ibin) + small_amt
    aer(icl_a,jliquid,ibin)  = aer(icl_a,jliquid,ibin)  + small_amt

    ! update jtotal
    aer(inh4_a,jtotal,ibin)  = aer(inh4_a,jsolid,ibin) +   &
         aer(inh4_a,jliquid,ibin)
    aer(icl_a,jtotal,ibin)   = aer(icl_a,jsolid,ibin)  +   &
         aer(icl_a,jliquid,ibin)

    ! update gas
    gas(inh3_g) = ((gas(inh3_g)) - (small_amt))
    gas(ihcl_g) = ((gas(ihcl_g)) - (small_amt))

    return
  end subroutine absorb_tiny_nh4cl

  !--------------------------------------------------------------------
  ! hno3 (liquid)
  subroutine absorb_tiny_hno3(ibin,aer,gas,delta_hno3_max)        ! and degas tiny hcl
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ngas_aerchtot,jliquid,jsolid,jtotal,                                      &
         icl_a,ino3_a,ihno3_g,ihcl_g
    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    real(r8), intent(inout), dimension(nbin_a_max) :: delta_hno3_max
    real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    ! local variables
    real(r8) :: small_aer, small_amt, small_gas

    small_gas = 0.01 * delta_hno3_max(ibin)
    small_aer = 0.01 * aer(icl_a,jliquid,ibin)

    small_amt = min(small_gas, small_aer)

    ! absorb tiny hno3
    aer(ino3_a,jliquid,ibin) = aer(ino3_a,jliquid,ibin) + small_amt
    aer(ino3_a,jtotal,ibin)  = aer(ino3_a,jsolid,ibin) +   &
         aer(ino3_a,jliquid,ibin)
    gas(ihno3_g) = ((gas(ihno3_g))-(small_amt))

    ! degas tiny hcl
    aer(icl_a,jliquid,ibin)  = ((aer(icl_a,jliquid,ibin))-   &
         (small_amt))
    aer(icl_a,jtotal,ibin)   = aer(icl_a,jsolid,ibin) +   &
         aer(icl_a,jliquid,ibin)

    ! update gas
    gas(ihcl_g) = gas(ihcl_g) + small_amt

    return
  end subroutine absorb_tiny_hno3

  !--------------------------------------------------------------
  ! hcl (liquid)
  subroutine absorb_tiny_hcl(ibin,aer,gas,delta_hcl_max)  ! and degas tiny hno3
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ngas_aerchtot,jliquid,jtotal,jsolid,                                      &
         ino3_a,icl_a,ihcl_g,ihno3_g
    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    real(r8), intent(inout), dimension(nbin_a_max) :: delta_hcl_max
    real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    ! local variables
    real(r8) :: small_aer, small_amt, small_gas

    small_gas = 0.01 * delta_hcl_max(ibin)
    small_aer = 0.01 * aer(ino3_a,jliquid,ibin)

    small_amt = min(small_gas, small_aer)

    ! absorb tiny hcl
    aer(icl_a,jliquid,ibin)= aer(icl_a,jliquid,ibin) + small_amt
    aer(icl_a,jtotal,ibin) = aer(icl_a,jsolid,ibin) +   &
         aer(icl_a,jliquid,ibin)
    gas(ihcl_g) = ((gas(ihcl_g))-(small_amt))

    ! degas tiny hno3
    aer(ino3_a,jliquid,ibin) = ((aer(ino3_a,jliquid,ibin))-   &
         (small_amt))
    aer(ino3_a,jtotal,ibin)  = aer(ino3_a,jsolid,ibin) +   &
         aer(ino3_a,jliquid,ibin)

    ! update gas
    gas(ihno3_g) = gas(ihno3_g) + small_amt

    return
  end subroutine absorb_tiny_hcl

  !--------------------------------------------------------------
  ! nh4no3 (liquid)
  subroutine degas_tiny_nh4no3(ibin,aer,gas,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ngas_aerchtot,jliquid,jsolid,jtotal,                                      &
         jnh4no3,inh4_a,ino3_a,inh3_g,ihno3_g
    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    ! local variables
    real(r8) :: small_amt

    small_amt = 0.01 * electrolyte(jnh4no3,jliquid,ibin)

    aer(inh4_a,jliquid,ibin) = ((aer(inh4_a,jliquid,ibin))-   &
         (small_amt))
    aer(ino3_a,jliquid,ibin) = ((aer(ino3_a,jliquid,ibin))-   &
         (small_amt))

    ! update jtotal
    aer(inh4_a,jtotal,ibin)  = aer(inh4_a,jsolid,ibin) +   &
         aer(inh4_a,jliquid,ibin)
    aer(ino3_a,jtotal,ibin)  = aer(ino3_a,jsolid,ibin) +   &
         aer(ino3_a,jliquid,ibin)

    ! update gas
    gas(inh3_g)  = gas(inh3_g)  + small_amt
    gas(ihno3_g) = gas(ihno3_g) + small_amt

    return
  end subroutine degas_tiny_nh4no3

  !--------------------------------------------------------------------
  ! nh4cl (liquid)
  subroutine degas_tiny_nh4cl(ibin,aer,gas,electrolyte)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ngas_aerchtot,jliquid,jsolid,jtotal,                                      &
         jnh4cl,inh4_a,icl_a,inh3_g,ihcl_g
    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    ! local variables
    real(r8) :: small_amt

    small_amt = 0.01 * electrolyte(jnh4cl,jliquid,ibin)

    aer(inh4_a,jliquid,ibin) = ((aer(inh4_a,jliquid,ibin))-   &
         (small_amt))
    aer(icl_a,jliquid,ibin)  = ((aer(icl_a,jliquid,ibin))-   &
         (small_amt))

    ! update jtotal
    aer(inh4_a,jtotal,ibin)  = aer(inh4_a,jsolid,ibin) +   &
         aer(inh4_a,jliquid,ibin)
    aer(icl_a,jtotal,ibin)   = aer(icl_a,jsolid,ibin)  +   &
         aer(icl_a,jliquid,ibin)

    ! update gas
    gas(inh3_g) = gas(inh3_g) + small_amt
    gas(ihcl_g) = gas(ihcl_g) + small_amt

    return
  end subroutine degas_tiny_nh4cl

  !***********************************************************************
  ! subroutines to equilibrate volatile acids
  !
  ! author: Rahul A. Zaveri
  ! update: may 2002
  !-----------------------------------------------------------------------
  subroutine equilibrate_acids(ibin,aer,gas,electrolyte,activity,mc,water_a,       &
       total_species,tot_cl_in,ma,gam,Keq_ll,Keq_gl)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ngas_aerchtot,ngas_volatile,Ncation,Nanion,nrxn_aer_gl,nrxn_aer_ll,       &
         ihno3_g,ihcl_g
    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    real(r8), intent(inout), dimension(nbin_a_max) :: water_a
    real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
    real(r8), intent(inout), dimension(ngas_volatile) :: total_species
    real(r8), intent(inout) :: tot_cl_in
    real(r8), intent(inout), dimension(nrxn_aer_gl) :: Keq_gl
    real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
    real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
    real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
    real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    ! local variables

    if(gas(ihcl_g)*gas(ihno3_g) .gt. 0.)then
       call equilibrate_hcl_and_hno3(ibin,aer,gas,electrolyte,activity,mc,water_a, &
            total_species,tot_cl_in,ma,gam,Keq_ll,Keq_gl)
    elseif(gas(ihcl_g) .gt. 0.)then
       call equilibrate_hcl(ibin,aer,gas,electrolyte,activity,mc,water_a,          &
            total_species,tot_cl_in,ma,gam,Keq_ll,Keq_gl)
    elseif(gas(ihno3_g) .gt. 0.)then
       call equilibrate_hno3(ibin,aer,gas,electrolyte,activity,mc,water_a,         &
            total_species,tot_cl_in,ma,gam,Keq_ll,Keq_gl)
    endif

    return
  end subroutine equilibrate_acids

  ! only hcl
  subroutine equilibrate_hcl(ibin,aer,gas,electrolyte,activity,mc,water_a,         &
       total_species,tot_cl_in,ma,gam,Keq_ll,Keq_gl)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ngas_aerchtot,ngas_volatile,Ncation,jliquid,jsolid,jtotal,Nanion,         &
         nrxn_aer_gl,nrxn_aer_ll,                                                  &
         ja_so4,ja_hso4,ihcl_g,icl_a,jhcl,ino3_a,ica_a,inh4_a,ina_a,jc_h,jc_ca,    &
         jc_nh4,jc_na,ja_cl,ja_no3,jhno3,jnh4cl

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    real(r8), intent(inout), dimension(nbin_a_max) :: water_a
    real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
    real(r8), intent(inout), dimension(ngas_volatile) :: total_species
    real(r8), intent(inout) :: tot_cl_in
    real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
    real(r8), intent(inout), dimension(nrxn_aer_gl) ::Keq_gl
    real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
    real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
    real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    ! local variables
    real(r8) :: a, aerH, aerHSO4, aerSO4, b, c, dum, Kdash_hcl, mH, Tcl,   &
         W, XT, Z
    !real(r8) :: quadratic                                 ! mosaic func

    aerSO4 = ma(ja_so4,ibin)*water_a(ibin)*1.e+9
    aerHSO4= ma(ja_hso4,ibin)*water_a(ibin)*1.e+9

    Tcl = aer(icl_a,jliquid,ibin) + gas(ihcl_g)   ! nmol/m^3(air)
    Kdash_hcl = Keq_gl(4)*1.e+18/gam(jhcl,ibin)**2        ! (nmol^2/kg^2)/(nmol/m^3(air))
    Z = (   aer(ina_a, jliquid,ibin) +               &  ! nmol/m^3(air)
         aer(inh4_a,jliquid,ibin) +   &
         2.*aer(ica_a, jliquid,ibin) ) -   &
         (2.*aerSO4  +   &
         aerHSO4 +   &
         aer(ino3_a,jliquid,ibin) )

    W     = water_a(ibin)                         ! kg/m^3(air)

    Kdash_hcl = Keq_gl(4)*1.e+18/gam(jhcl,ibin)**2        ! (nmol^2/kg^2)/(nmol/m^3(air))
    a = 1.0
    b = ((Kdash_hcl*W) + (Z/W))*1.e-9
    c = Kdash_hcl*(Z - Tcl)*1.e-18

    dum = ((b*b)-(4.*a*c))
    if (dum .lt. 0.) return               ! no real root

    if(c .lt. 0.)then
       mH = quadratic(a,b,c)      ! mol/kg(water)
       aerH = mH*W*1.e+9
       aer(icl_a,jliquid,ibin) = ((aerH) + (Z))
    else
       mH = sqrt(Keq_ll(3))
    endif

    call form_electrolytes(jliquid,ibin,XT,aer,gas,electrolyte,total_species,tot_cl_in)

    ! update gas phase concentration
    gas(ihcl_g) = ( (Tcl)  - (aer(icl_a,jliquid,ibin))  )

    ! update the following molalities
    ma(ja_so4,ibin)  = 1.e-9*aerSO4/water_a(ibin)
    ma(ja_hso4,ibin) = 1.e-9*aerHSO4/water_a(ibin)
    ma(ja_no3,ibin)  = 1.e-9*aer(ino3_a,jliquid,ibin)/water_a(ibin)
    ma(ja_cl,ibin)   = 1.e-9*aer(icl_a, jliquid,ibin)/water_a(ibin)

    mc(jc_h,ibin)    = mH
    mc(jc_ca,ibin)   = 1.e-9*aer(ica_a, jliquid,ibin)/water_a(ibin)
    mc(jc_nh4,ibin)  = 1.e-9*aer(inh4_a,jliquid,ibin)/water_a(ibin)
    mc(jc_na,ibin)   = 1.e-9*aer(ina_a, jliquid,ibin)/water_a(ibin)

    ! update the following activities
    activity(jhcl,ibin)    = mc(jc_h,ibin)  *ma(ja_cl,ibin)  *   &
         gam(jhcl,ibin)**2

    activity(jhno3,ibin)   = mc(jc_h,ibin)  *ma(ja_no3,ibin) *   &
         gam(jhno3,ibin)**2

    activity(jnh4cl,ibin)  = mc(jc_nh4,ibin)*ma(ja_cl,ibin) *   &
         gam(jnh4cl,ibin)**2

    ! also update xyz(jtotal)
    aer(icl_a,jtotal,ibin) = aer(icl_a,jliquid,ibin) +   &
         aer(icl_a,jsolid,ibin)

    electrolyte(jhcl,jtotal,ibin) = electrolyte(jhcl,jliquid,ibin)

    return
  end subroutine equilibrate_hcl

  ! only hno3
  subroutine equilibrate_hno3(ibin,aer,gas,electrolyte,activity,mc,water_a,        &
       total_species,tot_cl_in,ma,gam,Keq_ll,Keq_gl)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ngas_aerchtot,ngas_volatile,Ncation,jliquid,jsolid,jtotal,Nanion,         &
         nrxn_aer_gl,nrxn_aer_ll,                                                  &
         ja_so4,ja_hso4,ihno3_g,ino3_a,jhno3,icl_a,ica_a,inh4_a,ina_a,jc_h,jc_ca,  &
         jc_nh4,jc_na,ja_cl,jhcl,ja_no3,jnh4no3

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    real(r8), intent(inout), dimension(nbin_a_max) :: water_a
    real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
    real(r8), intent(inout), dimension(ngas_volatile) :: total_species
    real(r8), intent(inout) :: tot_cl_in
    real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
    real(r8), intent(inout), dimension(nrxn_aer_gl) :: Keq_gl
    real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
    real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
    real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    ! local variables
    real(r8) :: a, aerH, aerHSO4, aerSO4, b, c, dum, Kdash_hno3, mH,   &
         Tno3, W, XT, Z
    !real(r8) :: quadratic                                 ! mosaic func

    aerSO4 = ma(ja_so4,ibin)*water_a(ibin)*1.e+9
    aerHSO4= ma(ja_hso4,ibin)*water_a(ibin)*1.e+9

    Tno3 = aer(ino3_a,jliquid,ibin) + gas(ihno3_g)        ! nmol/m^3(air)
    Kdash_hno3 = Keq_gl(3)*1.e+18/gam(jhno3,ibin)**2      ! (nmol^2/kg^2)/(nmol/m^3(air))
    Z = (   aer(ina_a, jliquid,ibin) +               &  ! nmol/m^3(air)
         aer(inh4_a,jliquid,ibin) +   &
         2.*aer(ica_a, jliquid,ibin) ) -   &
         (2.*aerSO4  +   &
         aerHSO4 +   &
         aer(icl_a,jliquid,ibin) )

    W     = water_a(ibin)                         ! kg/m^3(air)

    Kdash_hno3 = Keq_gl(3)*1.e+18/gam(jhno3,ibin)**2      ! (nmol^2/kg^2)/(nmol/m^3(air))
    a = 1.0
    b = ((Kdash_hno3*W) + (Z/W))*1.e-9
    c = Kdash_hno3*(Z - Tno3)*1.e-18

    dum = ((b*b)-(4.*a*c))
    if (dum .lt. 0.) return               ! no real root

    if(c .lt. 0.)then
       mH = quadratic(a,b,c)      ! mol/kg(water)
       aerH = mH*W*1.e+9
       aer(ino3_a,jliquid,ibin) = ((aerH) + (Z))
    else
       mH = sqrt(Keq_ll(3))
    endif

    call form_electrolytes(jliquid,ibin,XT,aer,gas,electrolyte,total_species,tot_cl_in)

    ! update gas phase concentration
    gas(ihno3_g)= ( (Tno3) - (aer(ino3_a,jliquid,ibin)) )

    ! update the following molalities
    ma(ja_so4,ibin)  = 1.e-9*aerSO4/water_a(ibin)
    ma(ja_hso4,ibin) = 1.e-9*aerHSO4/water_a(ibin)
    ma(ja_no3,ibin)  = 1.e-9*aer(ino3_a,jliquid,ibin)/water_a(ibin)
    ma(ja_cl,ibin)   = 1.e-9*aer(icl_a, jliquid,ibin)/water_a(ibin)

    mc(jc_h,ibin)    = mH
    mc(jc_ca,ibin)   = 1.e-9*aer(ica_a, jliquid,ibin)/water_a(ibin)
    mc(jc_nh4,ibin)  = 1.e-9*aer(inh4_a,jliquid,ibin)/water_a(ibin)
    mc(jc_na,ibin)   = 1.e-9*aer(ina_a, jliquid,ibin)/water_a(ibin)

    ! update the following activities
    activity(jhcl,ibin)    = mc(jc_h,ibin)  *ma(ja_cl,ibin)  *   &
         gam(jhcl,ibin)**2

    activity(jhno3,ibin)   = mc(jc_h,ibin)  *ma(ja_no3,ibin) *   &
         gam(jhno3,ibin)**2

    activity(jnh4no3,ibin) = mc(jc_nh4,ibin)*ma(ja_no3,ibin) *   &
         gam(jnh4no3,ibin)**2

    ! also update xyz(jtotal)
    aer(ino3_a,jtotal,ibin) = aer(ino3_a,jliquid,ibin) +   &
         aer(ino3_a,jsolid,ibin)

    electrolyte(jhno3,jtotal,ibin) = electrolyte(jhno3,jliquid,ibin)

    return
  end subroutine equilibrate_hno3

  ! both hcl and hno3
  subroutine equilibrate_hcl_and_hno3(ibin,aer,gas,electrolyte,activity,mc,        &
       water_a,total_species,tot_cl_in,ma,gam,Keq_ll,Keq_gl)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ngas_aerchtot,ngas_volatile,Ncation,jliquid,jsolid,jtotal,Nanion,         &
         nrxn_aer_gl,nrxn_aer_ll,                                                  &
         ja_so4,ja_hso4,ihcl_g,icl_a,ihno3_g,ino3_a,jhcl,jhno3,             &
         ica_a,inh4_a,ina_a,jc_h,jc_ca,jc_nh4,jc_na,ja_cl,ja_no3,jnh4no3,   &
         jnh4cl

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    real(r8), intent(inout), dimension(nbin_a_max) :: water_a
    real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
    real(r8), intent(inout), dimension(ngas_volatile) :: total_species
    real(r8), intent(inout) :: tot_cl_in
    real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
    real(r8), intent(inout), dimension(nrxn_aer_gl) :: Keq_gl
    real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
    real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
    real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    ! local variables
    real(r8) :: aerH, aerHSO4, aerSO4, Kdash_hcl, Kdash_hno3,   &
         mH, p, q, r, Tcl, Tno3, W, XT, Z
    !real(r8) :: cubic                                     ! mosaic func

    aerSO4 = ma(ja_so4,ibin)*water_a(ibin)*1.e+9
    aerHSO4= ma(ja_hso4,ibin)*water_a(ibin)*1.e+9

    Tcl  = aer(icl_a,jliquid,ibin)  + gas(ihcl_g) ! nmol/m^3(air)
    Tno3 = aer(ino3_a,jliquid,ibin) + gas(ihno3_g)        ! nmol/m^3(air)

    Kdash_hcl  = Keq_gl(4)*1.e+18/gam(jhcl,ibin)**2       ! (nmol^2/kg^2)/(nmol/m^3(air))
    Kdash_hno3 = Keq_gl(3)*1.e+18/gam(jhno3,ibin)**2      ! (nmol^2/kg^2)/(nmol/m^3(air))

    Z = (   aer(ina_a, jliquid,ibin) +               &  ! nmol/m^3(air)
         aer(inh4_a,jliquid,ibin) +   &
         2.*aer(ica_a, jliquid,ibin) ) -   &
         (2.*aerSO4 + aerHSO4 )

    W = water_a(ibin)

    Kdash_hcl  = Keq_gl(4)*1.e+18/gam(jhcl,ibin)**2       ! (nmol^2/kg^2)/(nmol/m^3(air))
    Kdash_hno3 = Keq_gl(3)*1.e+18/gam(jhno3,ibin)**2      ! (nmol^2/kg^2)/(nmol/m^3(air))

    p = (Z/W + W*(Kdash_hcl + Kdash_hno3))*1.e-9

    q = 1.e-18*Kdash_hcl*Kdash_hno3*W**2  +   &
         1.e-18*Z*(Kdash_hcl + Kdash_hno3) -   &
         1.e-18*Kdash_hcl*Tcl -   &
         1.e-18*Kdash_hno3*Tno3

    r = 1.e-18*Kdash_hcl*Kdash_hno3*W*(Z - Tcl - Tno3)*1.e-9

    mH = cubic(p,q,r)

    if(mH .gt. 0.0)then
       aerH = mH*W*1.e+9
       aer(ino3_a,jliquid,ibin) = Kdash_hno3*W*W*Tno3/   &
            (aerH + Kdash_hno3*W*W)
       aer(icl_a, jliquid,ibin) = Kdash_hcl*W*W*Tcl/   &
            (aerH + Kdash_hcl*W*W)
    else
       mH = sqrt(Keq_ll(3))
    endif

    call form_electrolytes(jliquid,ibin,XT,aer,gas,electrolyte,total_species,tot_cl_in)

    ! update gas phase concentration
    gas(ihno3_g)= ( (Tno3) - (aer(ino3_a,jliquid,ibin)) )
    gas(ihcl_g) = ( (Tcl)  - (aer(icl_a,jliquid,ibin))  )

    ! update the following molalities
    ma(ja_so4,ibin)  = 1.e-9*aerSO4/water_a(ibin)
    ma(ja_hso4,ibin) = 1.e-9*aerHSO4/water_a(ibin)
    ma(ja_no3,ibin)  = 1.e-9*aer(ino3_a,jliquid,ibin)/water_a(ibin)
    ma(ja_cl,ibin)   = 1.e-9*aer(icl_a, jliquid,ibin)/water_a(ibin)

    mc(jc_h,ibin)    = mH
    mc(jc_ca,ibin)   = 1.e-9*aer(ica_a, jliquid,ibin)/water_a(ibin)
    mc(jc_nh4,ibin)  = 1.e-9*aer(inh4_a,jliquid,ibin)/water_a(ibin)
    mc(jc_na,ibin)   = 1.e-9*aer(ina_a, jliquid,ibin)/water_a(ibin)

    ! update the following activities
    activity(jhcl,ibin)    = mc(jc_h,ibin)*ma(ja_cl,ibin)   *   &
         gam(jhcl,ibin)**2

    activity(jhno3,ibin)   = mc(jc_h,ibin)*ma(ja_no3,ibin)  *   &
         gam(jhno3,ibin)**2

    activity(jnh4no3,ibin) = mc(jc_nh4,ibin)*ma(ja_no3,ibin)*   &
         gam(jnh4no3,ibin)**2

    activity(jnh4cl,ibin)  = mc(jc_nh4,ibin)*ma(ja_cl,ibin) *   &
         gam(jnh4cl,ibin)**2

    ! also update xyz(jtotal)
    aer(icl_a,jtotal,ibin)  = aer(icl_a,jliquid,ibin) +   &
         aer(icl_a,jsolid,ibin)

    aer(ino3_a,jtotal,ibin) = aer(ino3_a,jliquid,ibin) +   &
         aer(ino3_a,jsolid,ibin)

    electrolyte(jhno3,jtotal,ibin) = electrolyte(jhno3,jliquid,ibin)
    electrolyte(jhcl, jtotal,ibin) = electrolyte(jhcl, jliquid,ibin)

    return
  end subroutine equilibrate_hcl_and_hno3

  !***********************************************************************
  ! subroutines to evaporate solid volatile species
  !
  ! author: Rahul A. Zaveri
  ! update: sep 2004
  !-----------------------------------------------------------------------
  !
  ! nh4no3 (solid)
  subroutine degas_solid_nh4no3(ibin,aer,gas,electrolyte,Keq_sg)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ngas_aerchtot,jsolid,jliquid,jtotal,nrxn_aer_sg,                          &
         ihno3_g,inh3_g,jnh4no3,inh4_a,ino3_a

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
    real(r8), intent(inout), dimension(nrxn_aer_sg) :: Keq_sg
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    ! local variables
    integer jp
    real(r8) :: a, b, c, xgas, XT
    !real(r8) :: quadratic                                 ! mosaic func

    jp = jsolid

    a = 1.0
    b = gas(inh3_g) + gas(ihno3_g)
    c = gas(inh3_g)*gas(ihno3_g) - Keq_sg(1)
    xgas = quadratic(a,b,c)

    if(xgas .ge. electrolyte(jnh4no3,jp,ibin))then ! degas all nh4no3

       gas(inh3_g) = gas(inh3_g)  + electrolyte(jnh4no3,jp,ibin)
       gas(ihno3_g)= gas(ihno3_g) + electrolyte(jnh4no3,jp,ibin)
       aer(inh4_a,jp,ibin) = aer(inh4_a,jp,ibin) -   &
            electrolyte(jnh4no3,jp,ibin)
       aer(ino3_a,jp,ibin) = aer(ino3_a,jp,ibin) -   &
            electrolyte(jnh4no3,jp,ibin)

    else  ! degas only xgas amount of nh4no3

       gas(inh3_g) = gas(inh3_g)  + xgas
       gas(ihno3_g)= gas(ihno3_g) + xgas
       aer(inh4_a,jp,ibin) = aer(inh4_a,jp,ibin) - xgas
       aer(ino3_a,jp,ibin) = aer(ino3_a,jp,ibin) - xgas
    endif

    ! update jtotal
    aer(inh4_a,jtotal,ibin)  = aer(inh4_a,jsolid,ibin) +   &
         aer(inh4_a,jliquid,ibin)
    aer(ino3_a,jtotal,ibin)  = aer(ino3_a,jsolid,ibin) +   &
         aer(ino3_a,jliquid,ibin)

    return
  end subroutine degas_solid_nh4no3

  ! nh4cl (solid)
  subroutine degas_solid_nh4cl(ibin,aer,gas,electrolyte,Keq_sg)
    use module_data_mosaic_aero, only: r8,naer,nelectrolyte,nbin_a_max,            &
         ngas_aerchtot,jsolid,jliquid,jtotal,nrxn_aer_sg,                          &
         ihcl_g,inh3_g,jnh4cl,inh4_a,icl_a
    implicit none

    ! subr arguments
    integer, intent(in) :: ibin
    real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
    real(r8), intent(inout), dimension(nrxn_aer_sg) :: Keq_sg
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    ! local variables
    integer jp
    real(r8) :: a, b, c, xgas, XT
    !real(r8) :: quadratic                                 ! mosaic func

    jp = jsolid

    a = 1.0
    b = gas(inh3_g) + gas(ihcl_g)
    c = gas(inh3_g)*gas(ihcl_g) - Keq_sg(2)
    xgas = quadratic(a,b,c)

    if(xgas .ge. electrolyte(jnh4cl,jp,ibin))then ! degas all nh4cl

       gas(inh3_g) = gas(inh3_g) + electrolyte(jnh4cl,jp,ibin)
       gas(ihcl_g) = gas(ihcl_g) + electrolyte(jnh4cl,jp,ibin)
       aer(inh4_a,jp,ibin) = aer(inh4_a,jp,ibin) -   &
            electrolyte(jnh4cl,jp,ibin)
       aer(icl_a,jp,ibin)  = aer(icl_a,jp,ibin) -   &
            electrolyte(jnh4cl,jp,ibin)

    else  ! degas only xgas amount of nh4cl

       gas(inh3_g) = gas(inh3_g) + xgas
       gas(ihcl_g) = gas(ihcl_g) + xgas
       aer(inh4_a,jp,ibin) = aer(inh4_a,jp,ibin) - xgas
       aer(icl_a,jp,ibin)  = aer(icl_a,jp,ibin)  - xgas

    endif

    ! update jtotal
    aer(inh4_a,jtotal,ibin)  = aer(inh4_a,jsolid,ibin) +   &
         aer(inh4_a,jliquid,ibin)
    aer(icl_a,jtotal,ibin)   = aer(icl_a,jsolid,ibin)  +   &
         aer(icl_a,jliquid,ibin)

    return
  end subroutine degas_solid_nh4cl

  !***********************************************************************
  ! conforms aerosol generic species to a valid electrolyte composition
  !
  ! author: Rahul A. Zaveri
  ! update: june 2000
  !-----------------------------------------------------------------------
  subroutine conform_electrolytes(jp,ibin,XT,aer,gas,electrolyte,total_species,tot_cl_in)

    use module_data_mosaic_aero, only: r8,naer,nbin_a_max,           &
         ngas_aerchtot, ngas_volatile, nelectrolyte,                 &
         imsa_a,iso4_a,ica_a,ina_a,inh4_a,ino3_a,icl_a,ico3_a

    implicit none

    ! subr arguments
    integer, intent(in) :: ibin, jp
    real(r8), intent(inout) :: XT
    real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
    real(r8), intent(inout), dimension(ngas_volatile) :: total_species
    real(r8), intent(inout) :: tot_cl_in
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    ! local variables
    integer i, iXT_case, je
    real(r8) :: sum_dum, XNa_prime, XNH4_prime, XT_prime
    real(r8) :: store(naer)

    ! remove negative concentrations, if any
    !      do i=1,naer
    !      aer(i,jp,ibin) = max(0.0d0, aer(i,jp,ibin))    ! EFFI
    !      enddo

    !      call calculate_XT(ibin,jp,XT,aer)      ! EFFI

    if( (aer(iso4_a,jp,ibin)+aer(imsa_a,jp,ibin)) .gt.0.0)then
       XT   = ( aer(inh4_a,jp,ibin) +   &
            aer(ina_a,jp,ibin)  +   &
            2.*aer(ica_a,jp,ibin) )/   &
            (aer(iso4_a,jp,ibin)+0.5*aer(imsa_a,jp,ibin))
    else
       XT   = -1.0
    endif

!   if(XT .ge. 1.9999 .or. XT.lt.0.)then ! RAZ 11/10/2014
    if(XT .ge. 2.0 .or. XT.lt.0.)then ! RAZ 11/10/2014
       iXT_case = 1       ! near neutral (acidity is caused by HCl and/or HNO3)
    else
       iXT_case = 2       ! acidic (acidity is caused by excess SO4)
    endif

    ! initialize
    !
    ! put total aer(*) into store(*)
    store(iso4_a) = aer(iso4_a,jp,ibin)
    store(ino3_a) = aer(ino3_a,jp,ibin)
    store(icl_a)  = aer(icl_a, jp,ibin)
    store(imsa_a) = aer(imsa_a,jp,ibin)
    store(ico3_a) = aer(ico3_a,jp,ibin)
    store(inh4_a) = aer(inh4_a,jp,ibin)
    store(ina_a)  = aer(ina_a, jp,ibin)
    store(ica_a)  = aer(ica_a, jp,ibin)

    do je=1,nelectrolyte
       electrolyte(je,jp,ibin) = 0.0
    enddo

    !
    !---------------------------------------------------------
    !
    if(iXT_case.eq.1)then

       ! XT >= 2   : sulfate deficient

       call form_caso4(store,jp,ibin,electrolyte)
       call form_camsa2(store,jp,ibin,electrolyte)
       call form_na2so4(store,jp,ibin,electrolyte)
       call form_namsa(store,jp,ibin,electrolyte)
       call form_cano3(store,jp,ibin,electrolyte)
       call form_nano3(store,jp,ibin,electrolyte)
       call form_nacl(store,jp,ibin,aer,gas,electrolyte,total_species,tot_cl_in)
       call form_cacl2(store,jp,ibin,electrolyte)
       call form_caco3(store,jp,ibin,aer,electrolyte)
       call form_nh4so4(store,jp,ibin,electrolyte)
       call form_nh4msa(store,jp,ibin,electrolyte)
       call form_nh4no3(store,jp,ibin,electrolyte)
       call form_nh4cl(store,jp,ibin,electrolyte)
       call form_msa(store,jp,ibin,electrolyte)
       call degas_hno3(store,jp,ibin,aer,gas,electrolyte)
       call degas_hcl(store,jp,ibin,aer,gas,electrolyte)
       call degas_nh3(store,jp,ibin,aer,gas)

    elseif(iXT_case.eq.2)then

       ! XT < 2   : sulfate enough or sulfate excess

       call form_caso4(store,jp,ibin,electrolyte)
       call form_camsa2(store,jp,ibin,electrolyte)
       call form_namsa(store,jp,ibin,electrolyte)
       call form_nh4msa(store,jp,ibin,electrolyte)
       call form_msa(store,jp,ibin,electrolyte)

       if(store(iso4_a).eq.0.0)goto 10

       XT_prime =(store(ina_a)+store(inh4_a))/   &
            store(iso4_a)
       XNa_prime=0.5*store(ina_a)/store(iso4_a) + 1.

       if(XT_prime.ge.XNa_prime)then
          call form_na2so4(store,jp,ibin,electrolyte)
          XNH4_prime = 0.0
          if(store(iso4_a).gt.1.e-15)then
             XNH4_prime = store(inh4_a)/store(iso4_a)
          endif

          if(XNH4_prime .ge. 1.5)then
             call form_nh4so4_lvcite(store,jp,ibin,electrolyte)
          else
             call form_lvcite_nh4hso4(store,jp,ibin,electrolyte)
          endif

       elseif(XT_prime.ge.1.)then
          call form_nh4hso4(store,jp,ibin,electrolyte)
          call form_na2so4_nahso4(store,jp,ibin,electrolyte)
       elseif(XT_prime.lt.1.)then
          call form_nahso4(store,jp,ibin,electrolyte)
          call form_nh4hso4(store,jp,ibin,electrolyte)
          call form_h2so4(store,jp,ibin,electrolyte)
       endif

10     call degas_hno3(store,jp,ibin,aer,gas,electrolyte)
       call degas_hcl(store,jp,ibin,aer,gas,electrolyte)
       call degas_nh3(store,jp,ibin,aer,gas)

    endif ! case 1, 2

    ! re-calculate ions to eliminate round-off errors
    call electrolytes_to_ions(jp, ibin,aer,electrolyte)
    !---------------------------------------------------------
    !
    ! calculate % composition  EFFI
    !!      sum_dum = 0.0
    !!      do je = 1, nelectrolyte
    !!        electrolyte(je,jp,ibin) = max(0.d0,electrolyte(je,jp,ibin)) ! remove -ve
    !!        sum_dum = sum_dum + electrolyte(je,jp,ibin)
    !!      enddo
    !!
    !!      if(sum_dum .eq. 0.)sum_dum = 1.0
    !!      electrolyte_sum(jp,ibin) = sum_dum
    !!
    !!      do je = 1, nelectrolyte
    !!        epercent(je,jp,ibin) = 100.*electrolyte(je,jp,ibin)/sum_dum
    !!      enddo
    !!
    !!
    return
  end subroutine conform_electrolytes

  !----------------------------------------------------------
  ! solution to x^3 + px^2 + qx + r = 0
  !
  function cubic( psngl, qsngl, rsngl )
    use module_data_mosaic_kind, only:  r8
    implicit none
    real(r8) :: cubic
    ! subr arguments
    real(r8) :: psngl, qsngl, rsngl
    ! local variables
    real(r8) :: p, q, r, A, B, D, M, N, third, y
    real(r8) :: k, phi, thesign, x(3), duma
    integer icase, kk

    third = 1.d0/3.d0

    q = (qsngl)
    p = (psngl)
    r = (rsngl)

    A = (1.d0/3.d0)*((3.d0*q) - (p*p))
    B = (1.d0/27.d0)*((2.d0*p*p*p) - (9.d0*p*q) + (27.d0*r))

    D = ( ((A*A*A)/27.d0) + ((B*B)/4.d0) )

    if(D .gt. 0.)then     !       => 1 real and 2 complex roots
       icase = 1
    elseif(D .eq. 0.)then !       => 3 real roots, atleast 2 identical
       icase = 2
    else  ! D < 0         => 3 distinct real roots
       icase = 3
    endif

    goto (1,2,3), icase

    ! case 1: D > 0
1   thesign = 1.
    if(B .gt. 0.)then
       B = -B
       thesign = -1.
    endif

    M = thesign*((-B/2.d0) + (sqrt(D)))**(third)
    N = thesign*((-B/2.d0) - (sqrt(D)))**(third)

    cubic = ( (M) + (N) - (p/3.d0) )
    return

    ! case 2: D = 0
2   thesign = 1.
    if(B .gt. 0.)then
       B = -B
       thesign = -1.
    endif

    M = thesign*(-B/2.d0)**third
    N = M

    x(1) = ( (M) + (N) - (p/3.d0) )
    x(2) = ( (-M/2.d0) + (-N/2.d0) - (p/3.d0) )
    x(2) = ( (-M/2.d0) + (-N/2.d0) - (p/3.d0) )

    cubic = 0.
    do kk = 1, 3
       if(x(kk).gt.cubic) cubic = x(kk)
    enddo
    return

    ! case 3: D < 0
3   if(B.gt.0.)then
       thesign = -1.
    elseif(B.lt.0.)then
       thesign = 1.
    endif

    ! rce 18-nov-2004 -- make sure that acos argument is between +/-1.0
    !     phi = acos(thesign*sqrt( (B*B/4.d0)/(-A*A*A/27.d0) ))   ! radians
    duma = thesign*sqrt( (B*B/4.d0)/(-A*A*A/27.d0) )
    duma = min( duma, +1.0d0 )
    duma = max( duma, -1.0d0 )
    phi  = acos( duma )   ! radians

    cubic = 0.
    do kk = 1, 3
       k = kk-1
       y = 2.*Sqrt(-A/3.)*cos(phi + 120.*k*0.017453293)
       x(kk) = ((y) - (p/3.d0))
       if(x(kk).gt.cubic) cubic = x(kk)
    enddo
    return

  end function cubic
   !----------------------------------------------------------

  !----------------------------------------------------------
  function quadratic(a,b,c)
    use module_data_mosaic_kind, only:  r8
    implicit none
    real(r8) :: quadratic
    ! subr. arguments
    real(r8) :: a, b, c
    ! local variables
    real(r8) :: x, dum, quad1, quad2

    if(b .ne. 0.0)then
       x = 4.*(a/b)*(c/b)
    else
       x = 1.e+6
    endif

    if(abs(x) .lt. 1.e-6)then
       dum = ( (0.5*x) +   &
            (0.125*x**2) +   &
            (0.0625*x**3) )

       quadratic = (-0.5*b/a)*dum

       if(quadratic .lt. 0.)then
          quadratic = -b/a - quadratic
       endif

    else
       quad1 = ((-b)+sqrt((b*b)-(4.*a*c)))/   &
            (2.*a)
       quad2 = ((-b)-sqrt((b*b)-(4.*a*c)))/   &
            (2.*a)

       quadratic = max(quad1, quad2)
    endif

    return
  end function quadratic
  !----------------------------------------------------------

  !----------------------------------------------------------
  function mean_molecular_speed(T, MW)    ! in cm/s
    use module_data_mosaic_kind, only:  r8
    implicit none
    real(r8) :: mean_molecular_speed
    ! subr. arguments
    real(r8) :: T, MW     ! T(K)

    mean_molecular_speed = 1.455e4 * sqrt(T/MW)

    return
  end function mean_molecular_speed
  !----------------------------------------------------------

  !----------------------------------------------------------
  function gas_diffusivity(T, P, MW, Vm)  ! in cm^2/s
    use module_data_mosaic_kind, only:  r8
    use module_data_mosaic_constants, only:  third
    implicit none
    real(r8) :: gas_diffusivity
    ! subr. arguments
    real(r8) :: MW, Vm, T, P      ! T(K), P(atm)

    gas_diffusivity = (1.0e-3 * T**1.75 * sqrt(1./MW + 0.035))/   &
         (P * (Vm**third + 2.7189)**2)

    return
  end function gas_diffusivity
  !----------------------------------------------------------

  !----------------------------------------------------------
  function fuchs_sutugin(rkn,a)
    use module_data_mosaic_kind, only:  r8
    implicit none
    real(r8) :: fuchs_sutugin
    ! subr. arguments
    real(r8) :: rkn, a
    ! local variables
    real(r8) :: rnum, denom

    rnum  = 0.75*a*(1. + rkn)
    denom = rkn**2 + rkn + 0.283*rkn*a + 0.75*a
    fuchs_sutugin = rnum/denom

    return
  end function fuchs_sutugin
  !----------------------------------------------------------

  !----------------------------------------------------------
  ! ZSR method at 60% RH
  !
  function aerosol_water_up( ibin, electrolyte, aer, kappa_nonelectro, a_zsr ) ! kg (water)/m^3 (air)

    use module_data_mosaic_aero, only: r8, nelectrolyte, naer, nbin_a_max, jtotal, &
        nsalt, ioc_a, ibc_a, ilim2_a, ioin_a, dens_aer_mac, mw_aer_mac ! RAZ 4/16/2014

    use module_data_mosaic_asecthp, only: dens_water_aer

    implicit none

    real(r8) :: aerosol_water_up
    ! subr. arguments
    integer, intent(in) :: ibin
    real(r8), intent(in), dimension (6,nelectrolyte) :: a_zsr
    real(r8), intent(in), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    real(r8), intent(in), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(in), dimension(naer) :: kappa_nonelectro

    ! local variables
    integer :: iaer, jp, je
    real(r8) :: tmpa, tmpb, aH2O_60  ! RAZ 4/16/2014
    ! function
    !real(r8) :: bin_molality_60

    aH2O_60 = 0.6

    jp = jtotal
    tmpa = 0.0_r8

    do je = 1, (nsalt+4)  ! include hno3 and hcl in water calculation
       tmpa = tmpa + electrolyte(je,jp,ibin)/bin_molality_60(je,a_zsr)
    enddo

!   tmpa = tmpa + &
!        ( (aer(ilim2_a,jp,ibin)/dens_aer_mac(ilim2_a))*kappa_nonelectro(ilim2_a) + &                      ! RCE 5/20/2015
!          (aer(ioin_a, jp,ibin)/dens_aer_mac(ioin_a ))*kappa_nonelectro(ioin_a ) + &                      !  "     "
!          (aer(ioc_a,  jp,ibin)/dens_aer_mac(ioc_a  ))*kappa_nonelectro(ioc_a  ) + &                      !  "     "
!          (aer(ibc_a,  jp,ibin)/dens_aer_mac(ibc_a  ))*kappa_nonelectro(ibc_a  ) )*aH2O_60/(1.0-aH2O_60)  ! RCE 5/20/2015

    tmpb = 0.0_r8
    do iaer = 1, naer
       if (kappa_nonelectro(iaer) > 0.0_r8) then
          tmpb = tmpb + (aer(iaer,jp,ibin)*mw_aer_mac(iaer)/dens_aer_mac(iaer))*kappa_nonelectro(iaer)
       end if
    end do
    tmpa = tmpa + tmpb * dens_water_aer * 1.0e-3 * aH2O_60/(1.0-aH2O_60)   ! RCE 1/6/2016 added 1.0e-3 factor

    aerosol_water_up = tmpa*1.e-9  ! kg(water)/m^3(air)

    return
  end function aerosol_water_up
  !----------------------------------------------------------

  !----------------------------------------------------------
  function bin_molality_60(je,a_zsr)            ! TOUCH
    use module_data_mosaic_aero, only: r8,nelectrolyte

    implicit none

    real(r8) :: bin_molality_60
    ! subr. arguments
    integer, intent(in) ::  je
    real(r8), intent(in), dimension (6,nelectrolyte) :: a_zsr
    ! local variables
    real(r8) :: aw, xm

    aw = 0.6_r8

    xm =     a_zsr(1,je) +   &
         aw*(a_zsr(2,je) +   &
         aw*(a_zsr(3,je) +   &
         aw*(a_zsr(4,je) +   &
         aw*(a_zsr(5,je) +   &
         aw* a_zsr(6,je) ))))

    bin_molality_60 = 55.509_r8*xm/(1. - xm)

    return
  end function bin_molality_60
  !----------------------------------------------------------

  !----------------------------------------------------------
  ! ZSR method
  function aerosol_water( jp, ibin, jaerosolstate, jphase, jhyst_leg, electrolyte, aer,   &
           kappa_nonelectro, num_a, mass_dry_a, mass_soluble_a, aH2O, molality0 ) ! kg (water)/m^3 (air). RAZ added aer

    use module_data_mosaic_aero, only: r8, nbin_a_max, nelectrolyte, nsoluble, naer,   &
         all_solid, jsolid, jhyst_lo, ioc_a, ibc_a, ilim2_a, ioin_a,   &   ! RAZ 4/16/2014
         jtotal, ah2o_max, dens_aer_mac, mw_aer_mac, ename

    use module_data_mosaic_asecthp, only: dens_water_aer

    implicit none

    real(r8) :: aerosol_water
    ! subr. arguments
    integer, intent(in) :: jp, ibin
    integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate,jphase,jhyst_leg

    real(r8), intent(in) :: aH2O
    real(r8), intent(in), dimension(nbin_a_max) :: num_a,mass_dry_a,mass_soluble_a
    real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
    real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
    real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
    real(r8), intent(in), dimension(naer) :: kappa_nonelectro

    ! local variables
    integer :: iaer, iclm_aer, jclm_aer, je
    real(r8) :: tmpa, tmpb
    ! function
    real(r8) :: bin_molality

    tmpa = 0.0_r8
    if (jaerosolstate(ibin) .ne. all_solid) then                     ! RAZ 5/24/2016 - added this "if" check
! should this be "do je = 1, nsalt+4" as in aerosol_water_up ?
!             or "do je = 1, jxxxx" where jxxxx is the index of the last electrolyte used here ?
    do je = 1, 19   ! include hno3 and hcl in water calculation
       tmpa = tmpa + electrolyte(je,jp,ibin)/molality0(je,ibin)      ! RAZ 5/20/2014
    enddo
    endif

! note that this only considers the ilim2_a soa species

!   tmpa = tmpa + &
!       ( (aer(ioc_a,  jtotal,ibin)/dens_aer_mac(ioc_a  ))*kappa_nonelectro(ioc_a  ) +   &   ! RCE 5/20/2015
!         (aer(ilim2_a,jtotal,ibin)/dens_aer_mac(ilim2_a))*kappa_nonelectro(ilim2_a) +   &   ! RCE 5/20/2015
!         (aer(ioin_a, jtotal,ibin)/dens_aer_mac(ioin_a ))*kappa_nonelectro(ioin_a ) +   &   ! RCE 5/20/2015
!         (aer(ibc_a,  jtotal,ibin)/dens_aer_mac(ibc_a  ))*kappa_nonelectro(ibc_a  ) )   &   ! RCE 5/20/2015
!       * 1.0e-3 * aH2O/(1.0-min(ah2o,ah2o_max))                                             ! RCE 5/20/2015 - need 1.0e-3 factor

    tmpb = 0.0_r8
    do iaer = 1, naer
       if (kappa_nonelectro(iaer) > 0.0_r8) then
          tmpb = tmpb + (aer(iaer,jtotal,ibin)*mw_aer_mac(iaer)/dens_aer_mac(iaer))*kappa_nonelectro(iaer)
       end if
    end do
    tmpa = tmpa + tmpb * dens_water_aer * 1.0e-3 * ah2o/(1.0-min(ah2o,ah2o_max))  ! ug(water)/m^3(air)

    aerosol_water = tmpa*1.e-9  ! kg(water)/m^3(air)

    if (aerosol_water .le. 0.0) then !BALLI- Commented out to avoid slow runtime.
       iclm_aer = 0 !BSINGH- THIS IS WRONG!!!
       jclm_aer = 0 !BSINGH- THIS IS WRONG!!!
       !write(6,*)'iclm  jclm  ibin  jp = ',   &
       !     iclm_aer, jclm_aer, ibin, jp      !BSINGH- iclm_aer and jclm_aer are never set but they are used here.***

       !write(6,*)'aH2O, water = ', aH2O, aerosol_water
       !write(6,*)'dry mass = ', mass_dry_a(ibin)
       !write(6,*)'soluble mass = ', mass_soluble_a(ibin)
       !write(6,*)'number = ', num_a(ibin)
       !do je = 1, nsoluble
       !   write(6,44)ename(je), electrolyte(je,jp,ibin)
       !enddo
       !write(6,*)'Error in water calculation'
       !write(6,*)'ibin = ', ibin
       !write(6,*)'water content cannot be negative or zero'
       !write(6,*)'setting jaerosolstate to all_solid'

       !        call print_input

       jaerosolstate(ibin) = all_solid
       jphase(ibin)    = jsolid
       jhyst_leg(ibin) = jhyst_lo

    endif

44  format(a7, 2x, e11.3)

    return
  end function aerosol_water

  !----------------------------------------------------------
  function bin_molality(je,ibin,aH2O_a,b_zsr,a_zsr,aw_min)
    use module_data_mosaic_aero, only:r8,  nbin_a_max, nelectrolyte

    implicit none

    real(r8) :: bin_molality
    ! subr. arguments
    integer, intent(in) :: je, ibin
    real(r8), intent(in), dimension(nbin_a_max) :: aH2O_a
    real(r8), intent(in), dimension(nelectrolyte) :: b_zsr,aw_min
    real(r8), intent(in), dimension (6,nelectrolyte) :: a_zsr
    ! local variables
    real(r8) :: aw, xm

    aw = max(aH2O_a(ibin), aw_min(je))
    aw = min(aw, 0.999999_r8)

    if(aw .lt. 0.97_r8)then

       xm =     a_zsr(1,je) +   &
            aw*(a_zsr(2,je) +   &
            aw*(a_zsr(3,je) +   &
            aw*(a_zsr(4,je) +   &
            aw*(a_zsr(5,je) +   &
            aw* a_zsr(6,je) ))))

       bin_molality = 55.509_r8*xm/(1. - xm)

    else

       bin_molality = -b_zsr(je)*log(aw)

    endif

    return
  end function bin_molality
  !----------------------------------------------------------

end module module_mosaic_ext