! module_cam_mam_aerchem_driver.F
! created by r.c.easter, june 2010
!
! 2010-07-03 notes:
!
!--------------------------------------------------------------

      module module_cam_mam_aerchem_driver

      private
      public :: cam_mam_aerchem_driver

      contains


!==============================================================
      subroutine cam_mam_aerchem_driver(                            &
         id, curr_secs, ktau, dtstep, ktauc, dtstepc, config_flags, &
         t_phy, rho_phy, p_phy, p8w, alt, z, z_at_w, pbl_h, cldfra, &
         cldfra_mp_all, moist, chem,                                &
         dgnum, dgnumwet, wetdens_ap, del_h2so4_gasprod,            &
         dvmrdt_sv13d, dvmrcwdt_sv13d,                              & !Output from cloud chemistry
         ids,ide, jds,jde, kds,kde,                                 &
         ims,ime, jms,jme, kms,kme,                                 &
         its,ite, jts,jte, kts,kte                                  )
!-----------------------------------------------------------------------
! DESCRIPTION
!
! cam_mam_aerchem_driver is the interface between wrf-chem and the
!   cam_mam aerosol-chemistry routine cat computes condensation/evaporation
!   of trace gases to/from aerosol particles (AP).  It currently treats
!   water vapor and the 4 inorganic trace gases (nh3, h2so4, hno3, and hcl).
!   The aerosol-chemistry routine can work with either a sectional
!   (multiple size bins) or modal (multiple modes) representation.  
!
!   In both cases, condensation/evaporation to/from each bins/mode is 
!   first computed.  For sectional representation, AP mass and number 
!   are then transferred between size bins as a result of AP 
!   positive/negative growth.  Either a moving-center or two-moment
!   algorithm can be used to compute this transfer.
!
! cam_mam_aerchem_driver is organized as follows
!   loop over j and i
!      call aerchemistry to do the aerosol chemistry calculations
!          for timestep = dtstepc
!
!-----------------------------------------------------------------------

      use module_configure, only:  &
            grid_config_rec_type,                         &
            p_qv

      use module_state_description, only:  num_moist, num_chem, &
                                           param_first_scalar, CAMMGMPSCHEME

      use shr_kind_mod, only: r8 => shr_kind_r8
      use module_cam_support, only: pcols, pver, pverp, &
                                    pcnst => pcnst_runtime, &
                                    pcnst_non_chem => pcnst_non_chem_modal_aero, &
                                    gas_pcnst => gas_pcnst_modal_aero, &
                                    gas_pcnst_pos => gas_pcnst_modal_aero_pos, &
                                    endrun
      use constituents, only: cnst_name, cnst_get_ind
      use physconst, only: mwdry, pi

      use modal_aero_calcsize,    only:  modal_aero_calcsize_sub
      use modal_aero_coag,        only:  modal_aero_coag_sub
      use modal_aero_gasaerexch , only:  modal_aero_gasaerexch_sub
      use modal_aero_newnuc,      only:  modal_aero_newnuc_sub
      use modal_aero_wateruptake, only:  modal_aero_wateruptake_sub

      use modal_aero_data, only:  &
            alnsg_amode, dgnum_amode, dgnumlo_amode, dgnumhi_amode, &
            lmassptr_amode, lspectype_amode, &
            lptr_nacl_a_amode, lptr_so4_a_amode, lptr_soa_a_amode, &
            ntot_amode, nspec_amode, ntot_amode, numptr_amode, &
            specdens_amode, spechygro

      use module_data_cam_mam_asect, only:  &
            ai_phase, &
            factconv_chem_to_q, factconv_chem_to_qqcw, &
            lptr_chem_to_q, lptr_chem_to_qqcw, &
            massptr_aer, mw_q_array, mw_q_mo_array, &
            nsize_aer, ntype_aer, numptr_aer, waterptr_aer

      implicit none

!   subr arguments
      type(grid_config_rec_type), intent(in) :: config_flags
!   config_flags - configuration and control parameters

      integer, intent(in) ::             &
         id, ktau, ktauc,                &
         ids, ide, jds, jde, kds, kde,   &
         ims, ime, jms, jme, kms, kme,   &
         its, ite, jts, jte, kts, kte
!   id - domain index
!   ktau - time step number
!   ktauc - gas and aerosol chemistry time step number

!   [ids:ide, kds:kde, jds:jde] - spatial (x,z,y) indices for "domain"
!   [ims:ime, kms:kme, jms:jme] - spatial (x,z,y) indices for "memory"
!      Most arrays that are arguments to chem_driver 
!      are dimensioned with these spatial indices.
!   [its:ite, kts:kte, jts:jte] - spatial (x,z,y) indices for "tile"
!      chem_driver and routines under it do calculations
!      over these spatial indices.

      real(kind=8), intent(in) :: curr_secs
      real, intent(in) :: dtstep, dtstepc
!   dtstep - main model time step (s)
!   dtstepc - time step for gas and aerosol chemistry(s)

      real, intent(in),   &
         dimension( ims:ime, kms:kme, jms:jme ) :: &
         t_phy, rho_phy, p_phy, p8w, alt, &
         z, z_at_w, cldfra, cldfra_mp_all
!   t_phy - temperature at layer center (K)
!   rho_phy - air density at layer center (kg/m^3)
!   p_phy - air pressure at layer center (Pa)
!   p8w - air pressure at layer interfaces (Pa)
!   alt - dry air specific volume at layer center (m^3/kg)
!   z - height at layer center (m msl)
!   z_at_w - height at layer interfaces (m msl)
!   cldfra - cloud fraction of layer (--)
!   cldfra_mp_all- cloud fraction from CAMMGMP

      real, intent(in),   &
         dimension( ims:ime, jms:jme ) :: &
         pbl_h
!   pbl_h - pbl height (m agl)

      real, intent(in),   &
         dimension( ims:ime, kms:kme, jms:jme, 1:num_moist ) :: &
         moist
!   moist - mixing ratios of moisture species (water vapor, 
!      cloud water, ...) (kg/kg for mass species, #/kg for number species)
 
      real, intent(inout),   &
         dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
         chem
!   chem - mixing ratios of trace gase (ppm) and aerosol species
!      (ug/kg for mass species, #/kg for number species)

      real, intent(inout),   &
         dimension( ims:ime, kms:kme, jms:jme, ntot_amode ) :: &
         dgnum, dgnumwet, wetdens_ap
!   dgnum - median dry-diameter of number distribution for mode (m)
!   dgnumwet - median wet-diameter of number distribution for mode (m)
!   wetdens_ap - wet density of mode (kg/m^3)
      real, intent(in),   &
         dimension( ims:ime, kms:kme, jms:jme ) :: &
         del_h2so4_gasprod
!   del_h2so4_gasprod - h2so4 change from gas-phase chemistry (ppmv)

!
!tendencies:dvmrdt_sv13d,dvmrcwdt_sv13d are the tendencies which are passsed on from the CAM-MAM cloud chemistry
!           to gasaerexch subroutine 
!
    real, intent(in), dimension( ims:ime, kms:kme, jms:jme, 1:gas_pcnst_pos )   :: dvmrdt_sv13d,dvmrcwdt_sv13d 

!-----------------------------------------------------------------------
!   local variables
      integer :: aerchem_onoffb, aercoag_onoffb, aernewnuc_onoffb
      integer, parameter :: debug_level=0
      integer, parameter :: icol = 1
      integer, parameter :: idiag_calcsize=0, idiag_wateruptake=0
      integer, parameter :: idiag_gasaerexch=0, idiag_newnuc=0, idiag_coag=0
      integer :: i, icalcaer_flag, istat, itmpa, iwaterup_flag
      integer :: idiagaa, idiagaa_idel, idiagaa_jdel
      integer :: imozart, imozart_m1
      integer :: j
      integer :: k, kcam
      integer :: l, l1, l2, l3, lchnk, loffset
      integer :: l_mo_h2so4, l_mo_soag
      integer :: latndx(pcols), lonndx(pcols)
      integer :: n, ncol, nstep
      integer :: p1st
      
      real(r8) :: cld8(pcols,pver)
      real(r8) :: deltat8
      real(r8) :: del_h2so4_gasprod8(pcols,pver)
      real(r8) :: del_h2so4_aeruptk8(pcols,pver)
      real(r8) :: dgnum8(pcols,pver,ntot_amode)
      real(r8) :: dgnumwet8(pcols,pver,ntot_amode)
      real(r8) :: dqdt8(pcols,pver,pcnst)
      real(r8) :: dvmrdt_sv18(pcols,pver,gas_pcnst_pos), &
                  dvmrcwdt_sv18(pcols,pver,gas_pcnst_pos)
      real(r8) :: fconv
      real(r8) :: pblh8(pcols)
      real(r8) :: pdel8(pcols,pver), pmid8(pcols,pver)
      real(r8) :: qaerwat8(pcols,pver,ntot_amode)
      real(r8) :: q8(pcols,pver,pcnst), qqcw8(pcols,pver,pcnst)
      real(r8) :: qv8(pcols,pver)
      real(r8) :: rh8(pcols,pver)
      real(r8) :: tmpa, tmpb
      real(r8) :: tmpveca(gas_pcnst), tmpvecb(gas_pcnst)
      real(r8) :: tmpvmra(pver,gas_pcnst)
      real(r8) :: tmp_dgnum(ntot_amode), tmp_dstar(ntot_amode)
      real(r8) :: tmp_hygro(pver,ntot_amode)
      real(r8) :: tmp_num(ntot_amode), tmp_vdry(ntot_amode)
      real(r8) :: t8(pcols,pver)
      real(r8) :: vmr8(pcols,pver,gas_pcnst), vmrcw8(pcols,pver,gas_pcnst)
      real(r8) :: wetdens_ap8(pcols,pver,ntot_amode)
      real(r8) :: zm8(pcols,pver)

      logical :: aero_mmr_flag
      logical :: dotend(pcnst)
      logical :: h2o_mmr_flag

      character*100 msg


      itmpa = config_flags%aerchem_onoff
      if (itmpa <= 0) return

      aerchem_onoffb = 1
      aercoag_onoffb = 1
      aernewnuc_onoffb = 1
      if (itmpa == 99) then
         aerchem_onoffb = 0
      else if (itmpa >= 100) then
         aerchem_onoffb   = mod(itmpa/100,10)  ! 100s digit
         aernewnuc_onoffb = mod(itmpa/10, 10)  !  10s digit
         aercoag_onoffb   = mod(itmpa,    10)  !   1s digit
      end if
      if (aerchem_onoffb <= 0) then
         aercoag_onoffb = 0
         aernewnuc_onoffb = 0
      end if

      write(*,'(/a,3i5)') 'cam_mam_aerchem_driver - id, ktau, ktauc =', id, ktau, ktauc
      write(*,'(a,3(4x,2i5))') 'ids/e, j..., k... ', ids,ide, jds,jde, kds,kde
      write(*,'(a,3(4x,2i5))') 'ims/e, j..., k... ', ims,ime, jms,jme, kms,kme
      write(*,'(a,3(4x,2i5))') 'its/e, j..., k... ', its,ite, jts,jte, kts,kte
      write(*,'(a,3i10     )') 'pver, pverp, pcols       ', pver, pverp, pcols
      write(*,'(a,3i10     )') 'aerchem/newnuc/coag_onoffb', &
         aerchem_onoffb, aernewnuc_onoffb, aercoag_onoffb


      p1st = param_first_scalar
      imozart_m1 = pcnst_non_chem
      imozart = imozart_m1 + 1

      call cnst_get_ind( 'h2so4', l_mo_h2so4, .false. )
      l_mo_h2so4 = l_mo_h2so4 - imozart_m1
      if ((l_mo_h2so4 < 1) .or. (l_mo_h2so4 > gas_pcnst)) &
         call endrun( &
            'cam_mam_aerchem_driver error -- no h2so4 species' )
      write(*,*) 'l_mo_h2so4 = ', l_mo_h2so4

      call cnst_get_ind( 'soag', l_mo_soag, .false. )
      l_mo_soag = l_mo_soag - imozart_m1
      if ((l_mo_soag < 1) .or. (l_mo_soag > gas_pcnst)) &
         call endrun( &
            'cam_mam_aerchem_driver error -- no soag species' )
      write(*,*) 'l_mo_soag = ', l_mo_soag



!--------------------------------------------------------------------
!
! Main loop over the points in the tile and treat them each as a CAM chunk
!
!--------------------------------------------------------------------
      do 2920 j = jts, jte
!     write( msg, '(a,3i5)' ) 'cam_mam_aerchem_driver 2920 loop - j =', j
!     call wrf_debug( 100, msg )
      do 2910 i = its, ite

      lchnk = (j-jts)*(ite-its+1) + (i-its+1) !1-D index location from the 2-D tile
      ncol = 1

      latndx(icol) = j
      lonndx(icol) = i

      ! initialize these arrays to zero
      q8(:,:,:) = 0.0_r8
      qqcw8(:,:,:) = 0.0_r8
      qaerwat8(:,:,:) = 0.0_r8
      
      dgnum8(:,:,:) = 0.0_r8
      dgnumwet8(:,:,:) = 0.0_r8
      wetdens_ap8(:,:,:) = 0.0_r8

      ! pblh8 is msl, pbl_h is agl
      pblh8(icol) = pbl_h(i,j) + z_at_w(i,kts,j)
      ! force it to be at least the top of lowest layer
      pblh8(icol) = max( pblh8(icol), 1.0_r8*z_at_w(i,kts+1,j) )

      ! load the cam layer-center variables
      ! (reverse the vertical index, and convert units as needed)
      do k = kts, kte
         kcam = kte-k+1

!        cld8(icol,kcam)  = cldfra(i,k,j)
!        qh8(icol,kcam,1) = max( qv(i,k,j)/(1.+qv(i,k,j)), 1e-30 ) !values of 0 cause a crash in entropy
!        if( present(qc) ) then
!           ql8(icol,kcam) = qc(i,k,j)/(1.+qv(i,k,j)) !Convert to moist mix ratio
!        else
!           ql8(icol,kcam) = 0.
!        end if
!        if( present(qi) ) then
!           qi8(icol,kcam) = qi(i,k,j)/(1.+qv(i,k,j)) !Used in convtran, ditto for conversion
!        else
!           qi8(icol,kcam) = 0.
!        end if
!        zm8(icol,kcam)   = z(i,k,j)

         t8(icol,kcam)    = t_phy(i,k,j)
         pdel8(icol,kcam) = p8w(i,k,j) - p8w(i,k+1,j)
         pmid8(icol,kcam) = p_phy(i,k,j)
         zm8(icol,kcam)   = z(i,k,j)
         if(config_flags%mp_physics == CAMMGMPSCHEME) then
            cld8(icol,kcam) = cldfra_mp_all(i,k,j)
         else
            cld8(icol,kcam) = cldfra(i,k,j)
         endif
            
!         cld8(icol,kcam) = 0.0   ! *** temporary

         qv8(icol,kcam) = moist(i,k,j,p_qv)/(1.+moist(i,k,j,p_qv)) ! convert from dry to moist mix ratio

         do n = 1, ntype_aer
            dgnum8(icol,kcam,n) = dgnum(i,k,j,n)
            dgnumwet8(icol,kcam,n) = dgnumwet(i,k,j,n)
            wetdens_ap8(icol,kcam,n) = wetdens_ap(i,k,j,n)

            l = waterptr_aer(1,n)
            if ((l >= p1st) .and. (l <= num_chem)) &
               qaerwat8(icol,kcam,n) = chem(i,k,j,l)*factconv_chem_to_q(l)
         end do ! n

         do l = p1st, num_chem
            l2 = lptr_chem_to_q(l)
            if ((l2 >= 1) .and. (l2 <= pcnst)) then
               q8(icol,kcam,l2) = chem(i,k,j,l)*factconv_chem_to_q(l)
            end if
            if (i*j==1 .and. k==kts .and. l2==l_mo_h2so4+pcnst_non_chem) &
               write(*,'(a,1p,2e11.3)') '1,1,1 h2so4 chem & q8', &
               chem(i,k,j,l), q8(icol,kcam,l2)
            if (i*j==1 .and. k==kts .and. l2==l_mo_soag+pcnst_non_chem) &
               write(*,'(a,1p,2e11.3)') '1,1,1 soag  chem & q8', &
               chem(i,k,j,l), q8(icol,kcam,l2)
!            cycle ! this is temporary

            l2 = lptr_chem_to_qqcw(l)
            if ((l2 >= 1) .and. (l2 <= pcnst)) then
               qqcw8(icol,kcam,l2) = chem(i,k,j,l)*factconv_chem_to_q(l)
            end if
         end do ! l

         del_h2so4_gasprod8(icol,kcam) = del_h2so4_gasprod(i,k,j)
         if(config_flags%cldchem_onoff == 1 .and. config_flags%chem_opt >= 501 .AND. config_flags%chem_opt <= 504) then
            do l = 1 , gas_pcnst
               dvmrdt_sv18(1,kcam,l)   = dvmrdt_sv13d(i,k,j,l)
               dvmrcwdt_sv18(1,kcam,l) = dvmrcwdt_sv13d(i,k,j,l)
            enddo
         else
            dvmrdt_sv18(1,kcam,:)   = 0.0_r8
            dvmrcwdt_sv18(1,kcam,:) = 0.0_r8 
         endif
      end do ! k

      ! load the cam layer-interface variables
      ! (reverse the vertical index, and convert units as needed)
      do k = kts, kte+1
         kcam = kte-k+2

!        pint8(icol,kcam) = p8w(i,k,j)
!        zi8(icol,kcam)   = z_at_w(i,k,j)
      end do



!--------------------------------------------------------------------
!
! do modal_aero_calcsize_sub
!
!--------------------------------------------------------------------
      idiagaa = 1 ; idiagaa_idel = 100 ; idiagaa_jdel = 100

      nstep = ktau
      icalcaer_flag = 1
      loffset = 0             ! using entire q array, not the mozart "chem subset"
      aero_mmr_flag = .true.  ! q values are mass mixing ratios (kg/kg)

      deltat8 = dtstepc

      dotend(:) = .false.
      dqdt8(:,:,:) = 0.0_r8

      ! diagnostics for initial testing
      if (idiag_calcsize > 0) then
      if ((i==its) .and. (j==jts)) then
!        deltat8 = 86400.0*0.1   ! for testing, use bigger step size to see bigger changes
         do kcam = kts, kte, (kte-kts)
            k = kte-kcam+1
            tmp_num(1:ntot_amode) = q8(icol,kcam,numptr_amode(1:ntot_amode))
            ! adjust number mixrat here for testing
!           q8(icol,k,numptr_amode(1)) = q8(icol,k,numptr_amode(1)) * 1.0_r8
!           q8(icol,k,numptr_amode(2)) = q8(icol,k,numptr_amode(2)) / 8.0_r8
!           q8(icol,k,numptr_amode(3)) = q8(icol,k,numptr_amode(3)) * 8.0_r8
            do n = 1, ntot_amode
               ! calc dry volume
               tmp_vdry(n) = 0.0
               tmp_hygro(kcam,n) = 0.0
               do l1 = 1, nspec_amode(n)
                  l2 = lmassptr_amode(l1,n)
                  l3 = lspectype_amode(l1,n)
                  tmpa = q8(icol,kcam,l2)/specdens_amode(l3)
                  tmp_vdry(n) = tmp_vdry(n) + tmpa
                  tmp_hygro(kcam,n) = tmp_hygro(kcam,n) + tmpa*spechygro(l3)
               end do
               tmp_hygro(kcam,n) = tmp_hygro(kcam,n)/max(tmp_vdry(n),1.0e-35_r8)
               ! calc volume mean and number median diameter
               tmpa = tmp_vdry(n)*6.0/( pi*q8(icol,kcam,numptr_amode(n)) )
               tmp_dstar(n) = tmpa**(1.0/3.0)
               tmpb = alnsg_amode(n)
               tmp_dgnum(n) = tmp_dstar(n) * exp(-1.5*tmpb*tmpb)
            end do

            write(*,'(/a,4i5,1p,e12.4)') &
               'calcsize diagnostics before, i, j, k, kcam, deltat =', &
               i, j, k, kcam, deltat8
            write(*,'(a,1p,7e12.4))') 'tmp_vdry  ', &
               ( tmp_vdry(n), n=1,ntot_amode )
            write(*,'(a,1p,7e12.4))') 'tmp_dstar ', &
               ( tmp_dstar(n), n=1,ntot_amode )
            write(*,'(a,1p,7e12.4))') 'tmp_dgnum ', &
               ( tmp_dgnum(n), n=1,ntot_amode )
            write(*,'(a,1p,7e12.4))') 'dgnumlo   ', &
               ( dgnumlo_amode(n), n=1,ntot_amode )
            write(*,'(a,1p,7e12.4))') 'dgnumhi   ', &
               ( dgnumhi_amode(n), n=1,ntot_amode )
            write(*,'(a,1p,7e12.4))') 'dgnum     ', &
               ( dgnum8(icol,kcam,n), n=1,ntot_amode )
            write(*,'(a,1p,7e12.4))') 'number    ', &
               ( q8(icol,kcam,numptr_amode(n)), n=1,ntot_amode )
            write(*,'(a,1p,7e12.4))') 'tmp_num   ', &
               ( tmp_num(n), n=1,ntot_amode )
            write(*,'(a,1p,7e12.4))') 'q         '
            write(*,'(5(2x,a,1p,e12.4))') ( cnst_name(l)(1:6), &
               q8(icol,kcam,l), l=lmassptr_amode(1,1), pcnst )
         end do
      end if
      end if !(idiag_calcsize > 0)

      if (idiagaa > 0) then
         if (mod(i-its,idiagaa_idel)==0 .and. mod(j-jts,idiagaa_jdel)==0)   &
            print 93010, 'calling modal_aero_calcsize_sub - i,j =', i, j
      end if

!
! this is the call done in cam
!   call modal_aero_calcsize_sub(                                &
!                      lchnk,   ncol,    nstep,                  &
!                      1,                0,                      &
!                      .true.,           dt,                     &
!                      state%t,          state%pmid,             &
!                      state%pdel,       state%q,                &
!                      ptend%q,          ptend%lq,               &
!                      qqcw,             dgnum_pbuf              )
!
! this is the subr in wrfchem
!   subr------ modal_aero_calcsize_sub(                          &
!                      lchnk,   ncol,    nstep,   icalcaer_flag, &
!                      loffset,                                  &
!                      aero_mmr_flag,    deltat,                 &
!                      t,       pmid,    pdel,    q,             &
!                      dqdt,    dotend,  qqcw,    dgncur_a       )
!
      call modal_aero_calcsize_sub(                  &
           lchnk,   ncol,    nstep,   icalcaer_flag, &
           loffset,                                  &
           aero_mmr_flag,    deltat8,                &
           t8,      pmid8,   pdel8,   q8,            &
           dqdt8,   dotend,  qqcw8,   dgnum8         )

      if (idiagaa > 0) then
         if (mod(i-its,idiagaa_idel)==0 .and. mod(j-jts,idiagaa_jdel)==0)   &
            print 93010, 'backfrm modal_aero_calcsize_sub - i,j =', i, j
      end if

      ! diagnostics for initial testing
      if (idiag_calcsize > 0) then
      if ((i==its) .and. (j==jts)) then
         do kcam = kts, kte, (kte-kts)
            k = kte-kcam+1
            write(*,'(/a,4i5)') &
               'calcsize diagnostics after , i, j, k, kcam =', i, j, k, kcam
            write(*,'(a,1p,7e12.4))') 'dgnum     ', &
               ( dgnum8(icol,kcam,n), n=1,ntot_amode )
            write(*,'(a,1p,7e12.4))') 'number    ', &
               ( q8(icol,kcam,numptr_amode(n)), n=1,ntot_amode )
            write(*,'(a,1p,7e12.4))') ' +dt*tend ', &
               ( q8(icol,kcam,numptr_amode(n))+ &
                 deltat8*dqdt8(icol,kcam,numptr_amode(n)), n=1,ntot_amode )
            write(*,'(a,1p,7e12.4))') 'q+dt*tend '
            write(*,'(5(2x,a,1p,e12.4))') ( cnst_name(l)(1:6), &
               q8(icol,kcam,l)+deltat8*dqdt8(icol,kcam,l), l=lmassptr_amode(1,1), pcnst )
         end do
         write(*,'(a)') 
      end if
      end if !(idiag_calcsize > 0)

! apply tendencies
      do l = 1, pcnst
         if ( .not. dotend(l) ) cycle
         do k = 1, pver
            q8(icol,k,l) = q8(icol,k,l) + deltat8*dqdt8(icol,k,l)
            q8(icol,k,l) = max( q8(icol,k,l), 0.0_r8 )
         end do
      end do



!--------------------------------------------------------------------
!
! do modal_aero_wateruptake_sub
!
!--------------------------------------------------------------------
#ifdef OBSRH
! the modal_aero_wateruptake_sub OBSRH option is not implemented in wrf-chem
      msg = 'cam_mam_aerchem_driver calling modal_aero_wateruptake_sub' // &
            ' -- cannot do OBSRH'
      call wrf_error_fatal( msg )
#endif

      iwaterup_flag = 1
      h2o_mmr_flag = .true.

      dotend(:) = .false.
      dqdt8(:,:,:) = 0.0_r8
      rh8(:,:) = 0.0_r8

      ! diagnostics for initial testing
      if (idiag_wateruptake > 0) then
      if ((i==its) .and. (j==jts)) then
         do kcam = kts, kte, (kte-kts)
!           qv8(icol,kcam) = qv8(icol,kcam)*2.0   ! for testing
            k = kte-kcam+1
            write(*,'(/a,4i5,1p,e12.4)') &
               'wateruptake diagnostics before, i, j, k, kcam, deltat =', &
               i, j, k, kcam, deltat8
            write(*,'(a,1p,7e12.4))') 'tmp_hygro ', &
               ( tmp_hygro(kcam,n), n=1,ntot_amode )
            write(*,'(a,1p,7e12.4))') 'qaerwat   ', &
               ( qaerwat8(icol,kcam,n), n=1,ntot_amode )
            write(*,'(a,1p,7e12.4))') 'dgnum     ', &
               ( dgnum8(icol,kcam,n), n=1,ntot_amode )
            write(*,'(a,1p,7e12.4))') 'dgnumwet  ', &
               ( dgnumwet8(icol,kcam,n), n=1,ntot_amode )
            write(*,'(a,1p,7e12.4))') 'wetdens_ap', &
               ( wetdens_ap8(icol,kcam,n), n=1,ntot_amode )
         end do
      end if
      end if ! (idiag_wateruptake > 0)

      if (idiagaa > 0) then
         if (mod(i-its,idiagaa_idel)==0 .and. mod(j-jts,idiagaa_jdel)==0)   &
            print 93010, 'calling modal_aero_wateruptake_sub - i,j =', i, j
      end if

!
! this is the call done in cam
!   call modal_aero_wateruptake_sub(                             &
!                      lchnk,   ncol,    nstep,                  &
!                      1,                0,                      &
!                      .true.,           .true.,                 &
!                      dt,               state%q(:,:,1),         &
!                      state%t,          state%pmid,             &
!                      state%pdel, cldn, state%q,                &
!                      ptend%q,          ptend%lq,               &
!                      qaerwat,                                  &
!                      dgnum_pbuf,       dgnumwet_pbuf,          &
!                      wetdens_pbuf                              )
!
! this is the subr in wrfchem
!     subr------ modal_aero_wateruptake_sub(                &
!                lchnk, ncol, nstep,                        &
!                iwaterup_flag, loffset,                    &
!                aero_mmr_flag, h2o_mmr_flag,               &
!                deltat, h2ommr, t, pmid, pdel, cldn,       &
!                raer, raertend, dotend, qaerwat,           &
!                dgncur_a, dgncur_awet, wetdens             &
!#ifdef OBSRH
!                , obsrh                                    &
!#endif
!                                                           )
!
      call modal_aero_wateruptake_sub(                &
           lchnk, ncol, nstep,                        &
           iwaterup_flag, loffset,                    &
           aero_mmr_flag, h2o_mmr_flag,               &
           deltat8, qv8, t8, pmid8, pdel8, cld8,      &
           q8, dqdt8, dotend, qaerwat8,               &
           dgnum8, dgnumwet8, wetdens_ap8             )

      if (idiagaa > 0) then
         if (mod(i-its,idiagaa_idel)==0 .and. mod(j-jts,idiagaa_jdel)==0)   &
            print 93010, 'backfrm modal_aero_wateruptake_sub - i,j =', i, j
      end if

      ! diagnostics for initial testing
      if (idiag_wateruptake > 0) then
      if ((i==its) .and. (j==jts)) then
         do kcam = kts, kte, (kte-kts)
            k = kte-kcam+1
            write(*,'(/a,4i5,1p,e12.4)') &
               'wateruptake diagnostics after , i, j, k, kcam, rh =', &
               i, j, k, kcam, rh8(icol,kcam)
            write(*,'(a,1p,7e12.4))') 'qaerwat   ', &
               ( qaerwat8(icol,kcam,n), n=1,ntot_amode )
            write(*,'(a,1p,7e12.4))') 'dgnum     ', &
               ( dgnum8(icol,kcam,n), n=1,ntot_amode )
            write(*,'(a,1p,7e12.4))') 'dgnumwet  ', &
               ( dgnumwet8(icol,kcam,n), n=1,ntot_amode )
            write(*,'(a,1p,7e12.4))') 'wetdens_ap', &
               ( wetdens_ap8(icol,kcam,n), n=1,ntot_amode )
         end do
         write(*,'(a/(10f8.4))') 'rh(pver:1)', &
               ( rh8(icol,kcam), kcam=pver,1,-1 )
      end if
      end if ! (idiag_wateruptake > 0)



!--------------------------------------------------------------------
!
! do modal_aero_gasaerexch_sub
!
!--------------------------------------------------------------------
      if (aerchem_onoffb > 0) then

! transfer mass mixing ratios in q8 and qqcw8 to molar mixing ratios in vmr8 and vmrcw8
         do l2 = pcnst_non_chem+1, pcnst
            l3 = l2 - pcnst_non_chem
            fconv = mwdry/mw_q_array(l2)
            vmr8(  icol,:,l3) = q8(   icol,:,l2)*fconv
            vmrcw8(icol,:,l3) = qqcw8(icol,:,l2)*fconv
            if (i*j==1 .and. k==kts .and. l3==l_mo_h2so4) &
               write(*,'(a,1p,2e11.3)') '1,1,1 h2so4 q8 & vmr8', &
               q8(icol,pver,l2), vmr8(icol,pver,l3)
            if (i*j==1 .and. k==kts .and. l3==l_mo_soag) &
               write(*,'(a,1p,2e11.3)') '1,1,1 soag  q8 & vmr8', &
               q8(icol,pver,l2), vmr8(icol,pver,l3)
         end do

! also convert this 
         fconv = mwdry/mw_q_mo_array(l_mo_h2so4)
         del_h2so4_gasprod8(icol,:) = del_h2so4_gasprod8(icol,:)*fconv

! save current h2so4 vmr
         del_h2so4_aeruptk8(icol,:) = vmr8(icol,:,l_mo_h2so4)
         tmpvmra(:,:) = vmr8(icol,:,:)

! following is the cam code in mo_gas_phase_chemdr.F90
!        if (l_mo_h2so4 > 0) then
!           del_h2so4_aeruptk(1:ncol,:) = vmr(1:ncol,:,l_mo_h2so4)
!        else
!           del_h2so4_aeruptk(:,:) = 0.0_r8
!        endif
!        call modal_aero_gasaerexch_sub(                         &
!                          lchnk,    ncol,     nstep,            &
!                          imozart-1,          delt,             &
!                          latndx,   lonndx,                     &
!                          tfld,     pmid,     pdel,             &
!                          vmr,                vmrcw,            &
!                          dvmrdt_sv1,         dvmrcwdt_sv1,     &
!                          dgnum_pbuf,         dgnumwet_pbuf     )
!        if (l_mo_h2so4 > 0) then
!           del_h2so4_aeruptk(1:ncol,:) = vmr(1:ncol,:,l_mo_h2so4) - del_h2so4_aeruptk(1:ncol,:)
!        endif

! make call
         if (idiagaa > 0) then
            if (mod(i-its,idiagaa_idel)==0 .and. mod(j-jts,idiagaa_jdel)==0)   &
               print 93010, 'calling modal_aero_gasaerexch_sub - i,j =', i, j
         end if

         call modal_aero_gasaerexch_sub(                            &
                           lchnk,     ncol,       nstep,            &
                           imozart_m1,            deltat8,          &
                           latndx,    lonndx,                       &
                           t8,        pmid8,      pdel8,            &
                           vmr8,                  vmrcw8,           &
                           dvmrdt_sv18,           dvmrcwdt_sv18,    &
                           dgnum8,                dgnumwet8,        &
                           gas_pcnst_pos                             )

         if (idiagaa > 0) then
            if (mod(i-its,idiagaa_idel)==0 .and. mod(j-jts,idiagaa_jdel)==0)   &
               print 93010, 'backfrm modal_aero_gasaerexch_sub - i,j =', i, j
         end if

! diagnostic output after gasaerexch
         if (idiag_gasaerexch > 0) then
         if ((i == its) .and. (j == jts)) then
         do kcam = 1, pver
            if ((kcam /= pver) .and. (kcam /= pver-10)) cycle

            tmpveca(:) = 0.0
            tmpvecb(:) = 0.0
            do n = 1, ntot_amode
               l = lptr_so4_a_amode(n) - imozart_m1
               if ((l > 0) .and. (l <= gas_pcnst)) then
                  tmpveca(n+2) = tmpvmra(kcam,l)
                  tmpvecb(n+2) = vmr8(icol,kcam,l)
               end if
            end do
            tmpveca(2) = tmpvmra(kcam,l_mo_h2so4)
            tmpvecb(2) = vmr8(icol,kcam,l_mo_h2so4)
            tmpveca(1) = sum( tmpveca(2:ntot_amode+2) )
            tmpvecb(1) = sum( tmpvecb(2:ntot_amode+2) )
            write(*,'(/a,i3,1p,e13.5,10e11.3)') 'kcam, old so4 ...', &
               kcam, tmpveca(1:ntot_amode+2)
            write(*,'( a,i3,1p,e13.5,10e11.3)') 'kcam, new so4 ...', &
               kcam, tmpvecb(1:ntot_amode+2)

            tmpveca(:) = 0.0
            tmpvecb(:) = 0.0
            do n = 1, ntot_amode
               l = lptr_soa_a_amode(n) - imozart_m1
               if ((l > 0) .and. (l <= gas_pcnst)) then
                  tmpveca(n+2) = tmpvmra(kcam,l)
                  tmpvecb(n+2) = vmr8(icol,kcam,l)
               end if
            end do
            tmpveca(2) = tmpvmra(kcam,l_mo_soag)
            tmpvecb(2) = vmr8(icol,kcam,l_mo_soag)
            tmpveca(1) = sum( tmpveca(2:ntot_amode+2) )
            tmpvecb(1) = sum( tmpvecb(2:ntot_amode+2) )
            write(*,'(a,i3,1p,e13.5,10e11.3)') 'kcam, old soa ...', &
               kcam, tmpveca(1:ntot_amode+2)
            write(*,'(a,i3,1p,e13.5,10e11.3)') 'kcam, new soa ...', &
               kcam, tmpvecb(1:ntot_amode+2)

            tmpveca(:) = 0.0
            tmpvecb(:) = 0.0
            do n = 1, ntot_amode
               l = lptr_nacl_a_amode(n) - imozart_m1
               if ((l > 0) .and. (l <= gas_pcnst)) then
                  tmpveca(n+2) = tmpvmra(kcam,l)
                  tmpvecb(n+2) = vmr8(icol,kcam,l)
               end if
            end do
            tmpveca(1) = sum( tmpveca(2:ntot_amode+2) )
            tmpvecb(1) = sum( tmpvecb(2:ntot_amode+2) )
            write(*,'(a,i3,1p,e13.5,10e11.3)') 'kcam, old ncl ...', &
               kcam, tmpveca(1:ntot_amode+2)
            write(*,'(a,i3,1p,e13.5,10e11.3)') 'kcam, new ncl ...', &
               kcam, tmpvecb(1:ntot_amode+2)

            tmpveca(:) = 0.0
            tmpvecb(:) = 0.0
            do n = 1, ntot_amode
               l = numptr_amode(n) - imozart_m1
               if ((l > 0) .and. (l <= gas_pcnst)) then
                  tmpveca(n+2) = tmpvmra(kcam,l)
                  tmpvecb(n+2) = vmr8(icol,kcam,l)
               end if
            end do
            tmpveca(1) = sum( tmpveca(2:ntot_amode+2) )
            tmpvecb(1) = sum( tmpvecb(2:ntot_amode+2) )
            write(*,'(a,i3,1p,e13.5,10e11.3)') 'kcam, old num ...', &
               kcam, tmpveca(1:ntot_amode+2)
            write(*,'(a,i3,1p,e13.5,10e11.3)') 'kcam, new num ...', &
               kcam, tmpvecb(1:ntot_amode+2)

         end do ! kcam
         end if ! ((i == its) .and. (j == jts))
         end if ! (idiag_gasaerexch > 0)

! calc change to h2so4 vmr
         del_h2so4_aeruptk8(icol,:) = &
            vmr8(icol,:,l_mo_h2so4) - del_h2so4_aeruptk8(icol,:)

      end if ! (aerchem_onoffb > 0)



!--------------------------------------------------------------------
!
! do modal_aero_newnuc_sub
!
!--------------------------------------------------------------------
      if (aernewnuc_onoffb > 0) then
         tmpvmra(:,:) = vmr8(icol,:,:)

         if (idiagaa > 0) then
            if (mod(i-its,idiagaa_idel)==0 .and. mod(j-jts,idiagaa_jdel)==0)   &
               print 93010, 'calling modal_aero_newnuc_sub - i,j =', i, j
         end if

!     subr------ modal_aero_newnuc_sub(                             &
!                          lchnk,    ncol,     nstep,               &
!                          pcnstxx,  loffset,  deltat,              &
!                          latndx,   lonndx,                        &
!                          t,        pmid,     pdel,                &
!                          zm,       pblh,                          &
!                          qv,       cld,                           &
!                          q,                                       &
!                          del_h2so4_gasprod,  del_h2so4_aeruptk    )
         call modal_aero_newnuc_sub(                                &
                           lchnk,     ncol,       nstep,            &
                           imozart_m1,            deltat8,          &
                           latndx,    lonndx,                       &
                           t8,        pmid8,      pdel8,            &
                           zm8,       pblh8,                        &
                           qv8,       cld8,                         &
                           vmr8,                                    &
                           del_h2so4_gasprod8,  del_h2so4_aeruptk8, &
                           gas_pcnst                                )

         if (idiagaa > 0) then
            if (mod(i-its,idiagaa_idel)==0 .and. mod(j-jts,idiagaa_jdel)==0)   &
               print 93010, 'backfrm modal_aero_newnuc_sub - i,j =', i, j
         end if

! diagnostic output after newnuc
         if (idiag_newnuc > 0) then
         if ((i == its) .and. (j == jts)) then
         do kcam = 1, pver
            if ((kcam /= pver) .and. (kcam /= pver-10)) cycle

            tmpveca(:) = 0.0
            tmpvecb(:) = 0.0
            do n = 1, ntot_amode
               l = lptr_so4_a_amode(n) - imozart_m1
               if ((l > 0) .and. (l <= gas_pcnst)) then
                  tmpveca(n+2) = tmpvmra(kcam,l)
                  tmpvecb(n+2) = vmr8(icol,kcam,l)
               end if
            end do
            tmpveca(2) = tmpvmra(kcam,l_mo_h2so4)
            tmpvecb(2) = vmr8(icol,kcam,l_mo_h2so4)
            tmpveca(1) = sum( tmpveca(2:ntot_amode+2) )
            tmpvecb(1) = sum( tmpvecb(2:ntot_amode+2) )
            write(*,'(/a,i3,1p,e13.5,10e11.3)') 'kcam, old so4 ...', &
               kcam, tmpveca(1:ntot_amode+2)
            write(*,'( a,i3,1p,e13.5,10e11.3)') 'kcam, new so4 ...', &
               kcam, tmpvecb(1:ntot_amode+2)

            tmpveca(:) = 0.0
            tmpvecb(:) = 0.0
            do n = 1, ntot_amode
               l = lptr_soa_a_amode(n) - imozart_m1
               if ((l > 0) .and. (l <= gas_pcnst)) then
                  tmpveca(n+2) = tmpvmra(kcam,l)
                  tmpvecb(n+2) = vmr8(icol,kcam,l)
               end if
            end do
            tmpveca(2) = tmpvmra(kcam,l_mo_soag)
            tmpvecb(2) = vmr8(icol,kcam,l_mo_soag)
            tmpveca(1) = sum( tmpveca(2:ntot_amode+2) )
            tmpvecb(1) = sum( tmpvecb(2:ntot_amode+2) )
            write(*,'(a,i3,1p,e13.5,10e11.3)') 'kcam, old soa ...', &
               kcam, tmpveca(1:ntot_amode+2)
            write(*,'(a,i3,1p,e13.5,10e11.3)') 'kcam, new soa ...', &
               kcam, tmpvecb(1:ntot_amode+2)

            tmpveca(:) = 0.0
            tmpvecb(:) = 0.0
            do n = 1, ntot_amode
               l = lptr_nacl_a_amode(n) - imozart_m1
               if ((l > 0) .and. (l <= gas_pcnst)) then
                  tmpveca(n+2) = tmpvmra(kcam,l)
                  tmpvecb(n+2) = vmr8(icol,kcam,l)
               end if
            end do
            tmpveca(1) = sum( tmpveca(2:ntot_amode+2) )
            tmpvecb(1) = sum( tmpvecb(2:ntot_amode+2) )
            write(*,'(a,i3,1p,e13.5,10e11.3)') 'kcam, old ncl ...', &
               kcam, tmpveca(1:ntot_amode+2)
            write(*,'(a,i3,1p,e13.5,10e11.3)') 'kcam, new ncl ...', &
               kcam, tmpvecb(1:ntot_amode+2)

            tmpveca(:) = 0.0
            tmpvecb(:) = 0.0
            do n = 1, ntot_amode
               l = numptr_amode(n) - imozart_m1
               if ((l > 0) .and. (l <= gas_pcnst)) then
                  tmpveca(n+2) = tmpvmra(kcam,l)
                  tmpvecb(n+2) = vmr8(icol,kcam,l)
               end if
            end do
            tmpveca(1) = sum( tmpveca(2:ntot_amode+2) )
            tmpvecb(1) = sum( tmpvecb(2:ntot_amode+2) )
            write(*,'(a,i3,1p,e13.5,10e11.3)') 'kcam, old num ...', &
               kcam, tmpveca(1:ntot_amode+2)
            write(*,'(a,i3,1p,e13.5,10e11.3)') 'kcam, new num ...', &
               kcam, tmpvecb(1:ntot_amode+2)

         end do ! kcam
         end if ! ((i == its) .and. (j == jts))
         end if ! (idiag_newnuc > 0)

      end if ! (aernewnuc_onoffb > 0)



!--------------------------------------------------------------------
!
! do modal_aero_coag_sub
!
!--------------------------------------------------------------------
      if (aercoag_onoffb > 0) then
         tmpvmra(:,:) = vmr8(icol,:,:)

         if (idiagaa > 0) then
            if (mod(i-its,idiagaa_idel)==0 .and. mod(j-jts,idiagaa_jdel)==0)   &
               print 93010, 'calling modal_aero_coag_sub - i,j =', i, j
         end if

!     subr------ modal_aero_coag_sub(                               &
!                          lchnk,    ncol,     nstep,               &
!                          pcnstxx,  loffset,  deltat_main,         &
!                          latndx,   lonndx,                        &
!                          t,        pmid,     pdel,                &
!                          q,                                       &
!                          dgncur_a,           dgncur_awet,         &
!                          wetdens_a                                )
         call modal_aero_coag_sub(                                  &
                           lchnk,     ncol,       nstep,            &
                           imozart_m1,            deltat8,          &
                           latndx,    lonndx,                       &
                           t8,        pmid8,      pdel8,            &
                           vmr8,                                    &
                           dgnum8,                dgnumwet8,        &
                           wetdens_ap8,           gas_pcnst         )

         if (idiagaa > 0) then
            if (mod(i-its,idiagaa_idel)==0 .and. mod(j-jts,idiagaa_jdel)==0)   &
               print 93010, 'backfrm modal_aero_coag_sub - i,j =', i, j
         end if

! diagnostic output after coag
         if (idiag_coag > 0) then
         if ((i == its) .and. (j == jts)) then
         do kcam = 1, pver
            if ((kcam /= pver) .and. (kcam /= pver-10)) cycle

            tmpveca(:) = 0.0
            tmpvecb(:) = 0.0
            do n = 1, ntot_amode
               l = lptr_so4_a_amode(n) - imozart_m1
               if ((l > 0) .and. (l <= gas_pcnst)) then
                  tmpveca(n+2) = tmpvmra(kcam,l)
                  tmpvecb(n+2) = vmr8(icol,kcam,l)
               end if
            end do
            tmpveca(2) = tmpvmra(kcam,l_mo_h2so4)
            tmpvecb(2) = vmr8(icol,kcam,l_mo_h2so4)
            tmpveca(1) = sum( tmpveca(2:ntot_amode+2) )
            tmpvecb(1) = sum( tmpvecb(2:ntot_amode+2) )
            write(*,'(/a,i3,1p,e13.5,10e11.3)') 'kcam, old so4 ...', &
               kcam, tmpveca(1:ntot_amode+2)
            write(*,'( a,i3,1p,e13.5,10e11.3)') 'kcam, new so4 ...', &
               kcam, tmpvecb(1:ntot_amode+2)

            tmpveca(:) = 0.0
            tmpvecb(:) = 0.0
            do n = 1, ntot_amode
               l = lptr_nacl_a_amode(n) - imozart_m1
               if ((l > 0) .and. (l <= gas_pcnst)) then
                  tmpveca(n+2) = tmpvmra(kcam,l)
                  tmpvecb(n+2) = vmr8(icol,kcam,l)
               end if
            end do
            tmpveca(1) = sum( tmpveca(2:ntot_amode+2) )
            tmpvecb(1) = sum( tmpvecb(2:ntot_amode+2) )
            write(*,'(a,i3,1p,e13.5,10e11.3)') 'kcam, old ncl ...', &
               kcam, tmpveca(1:ntot_amode+2)
            write(*,'(a,i3,1p,e13.5,10e11.3)') 'kcam, new ncl ...', &
               kcam, tmpvecb(1:ntot_amode+2)

            tmpveca(:) = 0.0
            tmpvecb(:) = 0.0
            do n = 1, ntot_amode
               l = numptr_amode(n) - imozart_m1
               if ((l > 0) .and. (l <= gas_pcnst)) then
                  tmpveca(n+2) = tmpvmra(kcam,l)
                  tmpvecb(n+2) = vmr8(icol,kcam,l)
               end if
            end do
            tmpveca(1) = sum( tmpveca(2:ntot_amode+2) )
            tmpvecb(1) = sum( tmpvecb(2:ntot_amode+2) )
            write(*,'(a,i3,1p,e13.5,10e11.3)') 'kcam, old num ...', &
               kcam, tmpveca(1:ntot_amode+2)
            write(*,'(a,i3,1p,e13.5,10e11.3)') 'kcam, new num ...', &
               kcam, tmpvecb(1:ntot_amode+2)

         end do ! kcam
         end if ! ((i == its) .and. (j == jts))
         end if ! (idiag_coag > 0)

      end if ! (aercoag_onoffb > 0)



!--------------------------------------------------------------------
!
! after doing modal_aero_gasaerexch_sub (and newnuc, coag, ...),
! transfer final mixing ratios from vmr to q arrays
!
!--------------------------------------------------------------------
      if (aerchem_onoffb > 0) then

! transfer molar mixing ratios in vmr8 and vmrcw8 back to mass mixing ratios in q8 and qqcw8
         do l2 = pcnst_non_chem+1, pcnst
            l3 = l2 - pcnst_non_chem
            fconv = mw_q_array(l2)/mwdry
            q8(   icol,:,l2) = vmr8(  icol,:,l3)*fconv
            qqcw8(icol,:,l2) = vmrcw8(icol,:,l3)*fconv
            if (i*j==1 .and. k==kts .and. l3==l_mo_h2so4) &
               write(*,'(a,1p,2e11.3)') '1,1,1 h2so4 q8 & vmr8', &
               q8(icol,pver,l2), vmr8(icol,pver,l3)
            if (i*j==1 .and. k==kts .and. l3==l_mo_soag) &
               write(*,'(a,1p,2e11.3)') '1,1,1 soag  q8 & vmr8', &
               q8(icol,pver,l2), vmr8(icol,pver,l3)
         end do

      end if ! (aerchem_onoffb > 0)



!--------------------------------------------------------------------
!
! transfer updated data from cam to wrf-chem variables
!
!--------------------------------------------------------------------
      ! unload the cam layer-center variables that have been modified
      ! (reverse the vertical index, and convert units as needed)
      do k = kts, kte
         kcam = kte-k+1

         do n = 1, ntype_aer
            dgnum(i,k,j,n) = dgnum8(icol,kcam,n) 
            dgnumwet(i,k,j,n) = dgnumwet8(icol,kcam,n) 
            wetdens_ap(i,k,j,n) = wetdens_ap8(icol,kcam,n) 

            l = waterptr_aer(1,n)
            if ((l >= p1st) .and. (l <= num_chem)) &
               chem(i,k,j,l) = qaerwat8(icol,kcam,n)/factconv_chem_to_q(l)
         end do ! n

         do l = p1st, num_chem
            l2 = lptr_chem_to_q(l)
            if ((l2 >= 1) .and. (l2 <= pcnst)) then
               chem(i,k,j,l) = q8(icol,kcam,l2)/factconv_chem_to_q(l)
            end if
            if (i*j==1 .and. k==kts .and. l2==l_mo_h2so4+pcnst_non_chem) &
               write(*,'(a,1p,2e11.3)') '1,1,1 h2so4 chem & q8', &
               chem(i,k,j,l), q8(icol,kcam,l2)
            if (i*j==1 .and. k==kts .and. l2==l_mo_soag+pcnst_non_chem) &
               write(*,'(a,1p,2e11.3)') '1,1,1 soag  chem & q8', &
               chem(i,k,j,l), q8(icol,kcam,l2)
!            cycle

            l2 = lptr_chem_to_qqcw(l)
            if ((l2 >= 1) .and. (l2 <= pcnst)) then
               chem(i,k,j,l) = qqcw8(icol,kcam,l2)/factconv_chem_to_q(l)
            end if
         end do ! l

      end do



!--------------------------------------------------------------------
!
! End of main loop over the points in the tile (which treats them each as a CAM chunk)
!
!--------------------------------------------------------------------
2910  continue
2920  continue


      print 93010, 'leaving cam_mam_aerchem_driver - ktau =', ktau
93010 format( a, 8(1x,i6) )

      return
      end subroutine cam_mam_aerchem_driver


!-----------------------------------------------------------------------
   subroutine sum_pm_cam_mam (                                         &
         alt, chem,                                                    &
         pm2_5_dry, pm2_5_water, pm2_5_dry_ec, pm10,                   &
         ids,ide, jds,jde, kds,kde,                                    &
         ims,ime, jms,jme, kms,kme,                                    &
         its,ite, jts,jte, kts,kte                                     )

   USE module_state_description, only: num_chem
   USE module_data_mosaic_asect
   IMPLICIT NONE

   INTEGER,      INTENT(IN   )    ::                                   &
                                      ids,ide, jds,jde, kds,kde,       &
                                      ims,ime, jms,jme, kms,kme,       &
                                      its,ite, jts,jte, kts,kte

   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                       &
         INTENT(IN) :: alt

   REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),             &
         INTENT(IN ) :: chem

   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                       &
         INTENT(OUT) :: pm2_5_dry,pm2_5_water,pm2_5_dry_ec,pm10

   REAL :: mass

   INTEGER :: i,imax,j,jmax,k,kmax,n,itype,iphase

   imax = min(ite,ide-1)
   jmax = min(jte,jde-1)
   kmax = kte
!
! Sum over bins with center diameter < 2.5e-4 cm for pm2_5_dry,
! pm2_5_dry_ec, and pm2_5_water. All bins go into pm10
!
   pm2_5_dry(its:imax,kts:kmax,jts:jmax)    = 0.
   pm2_5_dry_ec(its:imax,kts:kmax,jts:jmax) = 0.
   pm2_5_water(its:imax,kts:kmax,jts:jmax)  = 0.
   pm10(its:imax,kts:kmax,jts:jmax)         = 0.

   do iphase=1,nphase_aer
   do itype=1,ntype_aer
   do n = 1, nsize_aer(itype)
      if (dcen_sect(n,itype) .le. 2.5e-4) then
         do j=jts,jmax
            do k=kts,kmax
               do i=its,imax
                  mass = chem(i,k,j,lptr_so4_aer(n,itype,iphase)) &
                       + chem(i,k,j,lptr_no3_aer(n,itype,iphase)) &
                       + chem(i,k,j,lptr_cl_aer(n,itype,iphase))  &
                       + chem(i,k,j,lptr_nh4_aer(n,itype,iphase)) &
                       + chem(i,k,j,lptr_na_aer(n,itype,iphase))  &
                       + chem(i,k,j,lptr_oin_aer(n,itype,iphase)) &
                       + chem(i,k,j,lptr_oc_aer(n,itype,iphase))  &
                       + chem(i,k,j,lptr_bc_aer(n,itype,iphase))
 
                  pm2_5_dry(i,k,j) = pm2_5_dry(i,k,j) + mass

                  pm2_5_dry_ec(i,k,j) = pm2_5_dry_ec(i,k,j)            &
                                      + chem(i,k,j,lptr_bc_aer(n,itype,iphase))

                  pm2_5_water(i,k,j) = pm2_5_water(i,k,j)              &
                                     + chem(i,k,j,waterptr_aer(n,itype))

                  pm10(i,k,j) = pm10(i,k,j) + mass
               enddo
            enddo
         enddo
      else
         do j=jts,jmax
            do k=kts,kmax
               do i=its,imax
                  pm10(i,k,j) = pm10(i,k,j)                              &
                              + chem(i,k,j,lptr_so4_aer(n,itype,iphase)) &
                          + chem(i,k,j,lptr_no3_aer(n,itype,iphase)) &
                          + chem(i,k,j,lptr_cl_aer(n,itype,iphase))  &
                          + chem(i,k,j,lptr_nh4_aer(n,itype,iphase)) &
                          + chem(i,k,j,lptr_na_aer(n,itype,iphase))  &
                          + chem(i,k,j,lptr_oin_aer(n,itype,iphase)) &
                          + chem(i,k,j,lptr_oc_aer(n,itype,iphase))  &
                          + chem(i,k,j,lptr_bc_aer(n,itype,iphase))
               enddo
            enddo
         enddo
      endif
   enddo ! size
   enddo ! type
   enddo ! phase

   !Convert the units from mixing ratio to concentration (ug m^-3)
   pm2_5_dry(its:imax,kts:kmax,jts:jmax) = pm2_5_dry(its:imax,kts:kmax,jts:jmax) &
                                           / alt(its:imax,kts:kmax,jts:jmax)
   pm2_5_dry_ec(its:imax,kts:kmax,jts:jmax) = pm2_5_dry_ec(its:imax,kts:kmax,jts:jmax) &
                                              / alt(its:imax,kts:kmax,jts:jmax)
   pm2_5_water(its:imax,kts:kmax,jts:jmax) = pm2_5_water(its:imax,kts:kmax,jts:jmax) &
                                             / alt(its:imax,kts:kmax,jts:jmax)
   pm10(its:imax,kts:kmax,jts:jmax) = pm10(its:imax,kts:kmax,jts:jmax) &
                                      / alt(its:imax,kts:kmax,jts:jmax)

   end subroutine sum_pm_cam_mam




!==============================================================
      end module module_cam_mam_aerchem_driver

