cam6_4_172: Aqueous chemistry fixes

src/chemistry/carma_aero/carma_aero_gasaerexch.F90

@@ -507,6 +507,9 @@ subroutine carma_aero_gasaerexch_sub(  state, &
                           num_bin,          t,          pmid,       &
                           wetr_n,                   uptkrate    )
 
+  ! initialize to zero before summing
+  uptkrate_all(:,:) = 0._r8
+
   do m = 1, nbins
 
      write(fieldname,'("NUMDENS_bin",I2.2)') m

src/chemistry/carma_aero/sox_cldaero_mod.F90

@@ -230,6 +230,7 @@ subroutine sox_cldaero_update(  &
          dso4dt_aqrxn, dso4dt_hprxn, &
          dso4dt_gasuptk, dmsadt_gasuptk_toso4, &
          dqdt_aq, dqdt_wr, dqdt
+    real(r8) :: delnh3
 
     real(r8) :: fwetrem, uptkrate
 
@@ -295,6 +296,10 @@ subroutine sox_cldaero_update(  &
 
                 delso4_o3rxn = xso4(i,k) - xso4_init(i,k)
 
+                if (id_nh3>0) then
+                   delnh3 = nh3g(i,k) - xnh3(i,k)
+                endif
+
                 ! the factors are proportional to the activated particle MR for each
                 ! bin, which is the MR of cloud drops "associated with" the mode
                 ! thus we are assuming the cloud drop size is independent of the
@@ -384,7 +389,6 @@ subroutine sox_cldaero_update(  &
                 dqdt_aq = -dso4dt_aqrxn*cldfrc(i,k)
                 dqdt = dqdt_aq + dqdt_wr
                 qin(i,k,id_so2) = qin(i,k,id_so2) + dqdt * dtime
-                qin(i,k,id_so2) =  MAX( qin(i,k,id_so2),    small_value )
 
                 ! h2o2 -- the first order loss rate for h2o2 is frh2o2_c*clwlrat(i,k)
                 ! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*frh2o2_c*clwlrat(i,k)))) )
@@ -394,7 +398,13 @@ subroutine sox_cldaero_update(  &
                 dqdt_aq = -dso4dt_hprxn*cldfrc(i,k)
                 dqdt = dqdt_aq + dqdt_wr
                 qin(i,k,id_h2o2) = qin(i,k,id_h2o2) + dqdt * dtime
-                qin(i,k,id_h2o2) =  MAX( qin(i,k,id_h2o2),    small_value )
+
+                ! NH3
+                if (id_nh3>0) then
+                   dqdt_aq = delnh3/dtime*cldfrc(i,k)
+                   dqdt = dqdt_aq
+                   qin(i,k,id_nh3) = qin(i,k,id_nh3) + dqdt * dtime
+                endif
 
                 ! for SO4 from H2O2/O3 budgets
                 dqdt_aqhprxn(i,k) = dso4dt_hprxn*cldfrc(i,k)
@@ -408,31 +418,52 @@ subroutine sox_cldaero_update(  &
     !==============================================================
     ! ... Update the mixing ratios
     !==============================================================
+    do k = 1,pver
+
+       do n = 1, nbins
+          do l = 1, nspec(n)
+             mm = bin_idx(n, l)
+             call rad_aer_get_bin_props_by_idx(0, n, l,spectype=spectype)
+             if (trim(spectype) == 'sulfate') then
+                qcw(:,k,mm) = max(qcw(:,k,mm), small_value )
+             end if
+          end do
+       end do
+
+       qin(:ncol,k,id_so2)   = max( qin(:ncol,k,id_so2),   small_value )
+       qin(:ncol,k,id_h2o2)  = max( qin(:ncol,k,id_h2o2),  small_value )
+       qin(:ncol,k,id_h2so4) = max( qin(:ncol,k,id_h2so4), small_value )
+       if ( id_nh3 > 0 ) qin(:ncol,k,id_nh3) = max( qin(:ncol,k,id_nh3), small_value )
+
+    end do
 
     ! diagnostics
 
     specmw_so4_amode = 96.0_r8
-      do n = 1, nbins
-        ! while looking through all species, only dqdt_aqso4 from sulfates  is gt zero
-        do l = 1, nspec(n)
-           mm = bin_idx(n, l)
-           aqso4(:,n)=0._r8
-            do k=1,pver
-               do i=1,ncol
-                  aqso4(i,n)=aqso4(i,n)+dqdt_aqso4(i,k,mm)*specmw_so4_amode/mbar(i,k) &
-                       *pdel(i,k)/gravit ! kg/m2/s
-               enddo
-            enddo
-
-            aqh2so4(:,n)=0._r8
-            do k=1,pver
-               do i=1,ncol
-                  aqh2so4(i,n)=aqh2so4(i,n)+dqdt_aqh2so4(i,k,mm)*specmw_so4_amode/mbar(i,k) &
-                       *pdel(i,k)/gravit ! kg/m2/s
-               enddo
-            enddo
-         end do
-      end do
+    do n = 1, nbins
+       ! while looking through all species, only dqdt_aqso4 from sulfates  is gt zero
+       do l = 1, nspec(n)
+          call rad_aer_get_bin_props_by_idx(0, n, l,spectype=spectype)
+          if (trim(spectype) == 'sulfate') then
+             mm = bin_idx(n, l)
+             aqso4(:,n)=0._r8
+             do k=1,pver
+                do i=1,ncol
+                   aqso4(i,n)=aqso4(i,n)+dqdt_aqso4(i,k,mm)*specmw_so4_amode/mbar(i,k) &
+                        *pdel(i,k)/gravit ! kg/m2/s
+                enddo
+             enddo
+
+             aqh2so4(:,n)=0._r8
+             do k=1,pver
+                do i=1,ncol
+                   aqh2so4(i,n)=aqh2so4(i,n)+dqdt_aqh2so4(i,k,mm)*specmw_so4_amode/mbar(i,k) &
+                        *pdel(i,k)/gravit ! kg/m2/s
+                enddo
+             enddo
+          end if
+       end do
+    end do
 
     aqso4_h2o2(:) = 0._r8
     do k=1,pver

src/chemistry/modal_aero/sox_cldaero_mod.F90

@@ -77,8 +77,10 @@ function sox_cldaero_create_obj(cldfrc, qcw, lwc, cfact, ncol, loffset) result(
     integer :: i,k
 
     logical :: mode7
+    logical :: mode5
 
     mode7 = ntot_amode == 7
+    mode5 = ntot_amode == 5
 
     conc_obj => cldaero_allocate()
 
@@ -127,6 +129,27 @@ function sox_cldaero_create_obj(cldfrc, qcw, lwc, cfact, ncol, loffset) result(
             + qcw(:ncol,:,id_nh4_4a) &
             + qcw(:ncol,:,id_nh4_5a) &
             + qcw(:ncol,:,id_nh4_6a)
+    else if (mode5) then
+#if ( defined MODAL_AERO_5MODE )
+       id_so4_1a = lptr_so4_cw_amode(1) - loffset
+       id_so4_2a = lptr_so4_cw_amode(2) - loffset
+       id_so4_3a = lptr_so4_cw_amode(3) - loffset
+       id_so4_5a = lptr_so4_cw_amode(5) - loffset
+#endif
+       conc_obj%so4c(:ncol,:) &
+            = qcw(:,:,id_so4_1a) &
+            + qcw(:,:,id_so4_2a) &
+            + qcw(:,:,id_so4_3a) &
+            + qcw(:,:,id_so4_5a)
+
+        ! for 3-mode, so4 is assumed to be nh4hso4
+        ! the partial neutralization of so4 is handled by using a
+        !    -1 charge (instead of -2) in the electro-neutrality equation
+       conc_obj%nh4c(:ncol,:) = 0._r8
+
+       ! with 3-mode, assume so4 is nh4hso4, and so half-neutralized
+       conc_obj%so4_fact = 1._r8
+
     else
        id_so4_1a = lptr_so4_cw_amode(1) - loffset
        id_so4_2a = lptr_so4_cw_amode(2) - loffset