Merge pull request #1548 from fvitt/aqueous_chemistry_fixes

fvitt · web-flow · commit 35d419d3689e · 2026-05-14T08:12:04.000-06:00
cam6_4_172: Aqueous chemistry fixes
diff --git a/doc/ChangeLog b/doc/ChangeLog
@@ -1,5 +1,93 @@
 ===============================================================
 
+Tag name: cam6_4_172
+Originator(s): fvitt
+Date: 14 May 2026
+One-line Summary: Aqueous chemistry bug fixes
+Github PR URL: https://github.com/ESCOMP/CAM/pull/1548
+
+Purpose of changes (include the issue number and title text for each relevant GitHub issue):
+
+Issue #1544
+Fix the following aqueous chemistry bugs:
+ - 5th mode is not included in the updates when MAM5 is used
+ - Failure to update NH3 when CARMA is used
+
+Describe any changes made to build system: N/A
+
+Describe any changes made to the namelist: N/A
+
+List any changes to the defaults for the boundary datasets: N/A
+
+Describe any substantial timing or memory changes: N/A
+
+Code reviewed by: cacraigucar tilmes
+
+List all files eliminated:
+
+List all files added and what they do:
+
+List all existing files that have been modified, and describe the changes:
+M       src/chemistry/carma_aero/carma_aero_gasaerexch.F90
+ - initialize diagnostic varible to zero
+
+M       src/chemistry/carma_aero/sox_cldaero_mod.F90
+ - include NH3 gas loss
+ - enforce lower limits
+
+M       src/chemistry/modal_aero/sox_cldaero_mod.F90
+ - include the 5th sulfate mode
+
+If there were any failures reported from running test_driver.sh on any test
+platform, and checkin with these failures has been OK'd by the gatekeeper,
+then copy the lines from the td.*.status files for the failed tests to the
+appropriate machine below.  All failed tests must be justified.
+
+derecho/intel/aux_cam:
+  DIFF ERP_Ld3.ne16pg3_ne16pg3_mg17.FHISTC_WAt1ma.derecho_intel.cam-reduced_hist1d
+  DIFF ERP_Ld3.ne30pg3_ne30pg3_mt232.FHISTC_MTt4s.derecho_intel.cam-outfrq1d_aoa
+  DIFF ERP_Ln9.f09_f09_mg17.FCSD_HCO.derecho_intel.cam-outfrq9s
+  DIFF ERP_Ln9.ne30pg3_ne30pg3_mg17.FCnudged.derecho_intel.cam-outfrq9s
+  DIFF ERP_Ln9.ne30pg3_ne30pg3_mg17.FHISTC_WAma.derecho_intel.cam-outfrq9s
+  DIFF ERS_Ln9.f09_f09_mg17.FX2000.derecho_intel.cam-outfrq9s
+  DIFF ERS_Ln9.f19_f19_mg17.FXSD.derecho_intel.cam-outfrq9s
+  DIFF ERS_Ln9.ne30pg3_ne30pg3_mg17.FHISTC_WXma.derecho_intel.cam-outfrq9s_ctem
+  DIFF SMS_C2_D_Ln9.ne16pg3_ne16pg3_mg17.FHISTC_WXma.derecho_intel.cam-outfrq9s
+  DIFF SMS_D_Ln9.f09_f09_mg17.FCts2nudged.derecho_intel.cam-outfrq9s_leapday
+  DIFF SMS_D_Ln9.f09_f09_mg17.FCvbsxHIST.derecho_intel.cam-outfrq9s
+  DIFF SMS_D_Ln9.f19_f19_mg17.FCARMA2000climo.derecho_intel.cam-outfrq9s
+  DIFF SMS_D_Ln9.f19_f19_mg17.FCHIST_SLH.derecho_intel.cam-outfrq9s
+  DIFF SMS_D_Ln9.f19_f19_mg17.FWma2000climo.derecho_intel.cam-outfrq9s
+  DIFF SMS_D_Ln9.f19_f19_mg17.QPC2000climo.derecho_intel.cam-outfrq3s_usecase
+  DIFF SMS_D_Ln9.ne30pg3_ne30pg3_mt232.1850C_CAM70%MT%CT4S2_CLM60%SP_CICE%PRES_DOCN%DOM_MOSART_SGLC_SWAV_SESP.derecho_intel.cam-outfrq9s
+  DIFF SMS_D_Ln9_P1280x1.ne30pg3_ne30pg3_mt232.FHISTC_MTt1s.derecho_intel.cam-outfrq9s_Leung_dust
+  DIFF SMS_D_Ln9_P1536x1.ne0CONUSne30x8_ne0CONUSne30x8_mt12.FCHIST.derecho_intel.cam-outfrq9s
+  DIFF SMS_Ld1.f09_f09_mg17.FW2000climo.derecho_intel.cam-outfrq1d
+  DIFF SMS_Ld1.ne30pg3_ne30pg3_mg17.FC2010climo.derecho_intel.cam-outfrq1d
+  DIFF SMS_Lh12.f09_f09_mg17.FCSD_HCO.derecho_intel.cam-outfrq3h
+  DIFF SMS_Ln9.f09_f09_mg17.FW1850.derecho_intel.cam-reduced_hist3s
+  DIFF SMS_Ln9.ne30pg3_ne30pg3_mg17.FW2000climo.derecho_intel.cam-outfrq9s_rrtmgp
+  - expected differences due to bug fixes for MAM5 and CARMA
+
+derecho/nvhpc/aux_cam: PASS
+
+izumi/nag/aux_cam:
+  FAIL ERC_D_Ln9.f10_f10_mt232.FHIST_C5.izumi_nag.cam-outfrq3s_subcol COMPARE_base_rest
+  - pre-existing failure -- see https://github.com/ESCOMP/CAM/issues/1514
+
+  DIFF ERC_D_Ln9.f10_f10_mg37.QPWmaC6.izumi_nag.cam-outfrq3s
+  - expected differences due to bug fix for MAM5
+
+izumi/gnu/aux_cam:
+  DIFF SMS_D_Ln9.f10_f10_mg37.FWmaHIST.izumi_gnu.cam-outfrq9s_mee_fluxes
+  - expected differences due to bug fix for MAM5
+
+Summarize any changes to answers:
+  larger than roundoff for MAM5 and CARMA configurations
+
+===============================================================
+===============================================================
+
 Tag name: cam6_4_171
 Originator(s): fvitt
 Date: 1 May 2026
diff --git a/src/chemistry/carma_aero/carma_aero_gasaerexch.F90 b/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
diff --git a/src/chemistry/carma_aero/sox_cldaero_mod.F90 b/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
diff --git a/src/chemistry/modal_aero/sox_cldaero_mod.F90 b/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
