Three bugs in micro_pumas_v1.F90 / micro_pumas_utils.F90 (accre_sees_auto droplet number, vtrmi_factor, ice lambda lower bound)

Start by comparing micro_pumas_v1.F90 around lines 2442 and 3350–3370 with the issue’s proposed changes, then inspect micro_pumas_utils.F90 around line 379 and the earlier v1.17 form. The report describes tested climate runs but names no test files or commands. Done means correcting all three reported issues and checking their process-level effects, including the noted vtrmi_factor tuning caveat.

Hi all, some bugs were identified with Claude by Brian Dobbins, and I (and Claude) confirmed and investigated.

Found in pumas_cam-release_v1.41 (cesm3_0_alpha11c); the same lines are unchanged in main_cam and pumas_cam-release_v1.42. All three are active in the default CAM7 configuration (micro_mg_warm_rain='kk2000', micro_mg_accre_sees_auto=.true., micro_mg_vtrmi_factor=1.39).

1. accre_sees_auto subtracts the rain-number source instead of the droplet-number sink

micro_pumas_v1.F90 ~L2442:

ntmp(i,k) = ntmp(i,k) - nprc(i,k)*deltat    ! nprc = rain-number gain (prc/droplet_mass_25um)

The droplet-number loss from autoconversion is nprc1 (= prc*ncic/qcic). The rest of the code agrees: it is used in the nc conservation check and mapped from nctend_KK2000. With the bug, npra = pra*ntmp/ctmp is inflated by (1 − f·m̄/m25)/(1 − f), where f = prc·Δt/qcic. Fix: use nprc1. Present since the option was added (1e3466f/d69eff4, July 2021).

2. micro_mg_vtrmi_factor applied twice, and never to number

micro_pumas_v1.F90 ~L3350–3370: the factor multiplies vtrmi once before the small/large-ice blend and again after it. For 18 < r < 20 µm the large-ice term therefore carries 1.39² = 1.93, and the fall speed jumps at r = 20 µm. fni (the number-weighted speed) is never scaled, so mass and number sediment with inconsistent scalings. Fix: blend the unscaled speeds, then apply the factor once to both fi and fni. Same code in v1.17.

3. Ice slope lower bound is dcs/2 instead of 1/(2*dcs)

micro_pumas_utils.F90 L379:

ice_lambda_bounds(1) = 1._r8/2._r8*dcs      ! evaluates to dcs/2 = 2.25e-4 m^-1

The intended bound is 1/(2·dcs), a 900 µm ceiling on the mean ice diameter for dcs = 450 µm. As written the ceiling never engages, and the model carries ice with mean sizes of several cm and near-zero number. v1.17 had the correct form 1._r8/[2._r8*dcs, 1.e-6_r8]. This is a regression from 5434e2b (2022-05-12, OpenACC merge). Fix: 1._r8/(2._r8*dcs).

Tested impact (F2000dev, ne30pg3; 3 control + 3 all-fix members × 6 months)

• Process level: each fix works as intended.
  • Ice cells with mean diameter ≥ 1 mm go from ~6×10⁻⁵ of ice-bearing cells to 0.
  • Ice-number sedimentation shifts downward, and REI drops 3–10 µm at 400–800 hPa in the tropics.
  • The accretion number/mass ratio no longer grows with autoconversion fraction.
• Climate, fixed SST: small.
  • IWP +2.7% (p≈0.01–0.07).
  • LWP +1.4% (p≈0.04).
  • Net TOA −0.16 ± 0.24 W m⁻² (1 SE), not significant.
  • No field-significant map changes.
• Caveat: vtrmi_factor=1.39 was tuned with fni unscaled, so it may be worth revisiting.

Proposed fix

--- a/micro_pumas_v1.F90
+++ b/micro_pumas_v1.F90
@@ -2439,7 +2439,9 @@
            if (accre_sees_auto) then
               rtmp(i,k) = rtmp(i,k) + prc(i,k)*deltat
               ctmp(i,k) = ctmp(i,k) - prc(i,k)*deltat
-              ntmp(i,k) = ntmp(i,k) - nprc(i,k)*deltat
+              ! remove the droplet-number sink (nprc1), not the
+              ! rain-number source (nprc), from the provisional droplet number
+              ntmp(i,k) = ntmp(i,k) - nprc1(i,k)*deltat
            endif
         end do
      end do
@@ -3349,9 +3351,7 @@
         if (dumi(i,k).ge.qsmall) then
            proc_rates%vtrmi(i,k)=min(ain(i,k)*gamma_bi_plus4/(6._r8*lami(i,k)**bi), &
                 1.2_r8*rhof(i,k))
-           proc_rates%vtrmi(i,k)=proc_rates%vtrmi(i,k)*micro_mg_vtrmi_factor

-           fi(i,k) = g*rho(i,k)*proc_rates%vtrmi(i,k)
            fni(i,k) = g*rho(i,k)* &
                 min(ain(i,k)*gamma_bi_plus1/lami(i,k)**bi,1.2_r8*rhof(i,k))

@@ -3365,15 +3365,19 @@
                  (1._r8 - ifrac) * &
                  min(ajn(i,k)*gamma_bj_plus4/(6._r8*lami(i,k)**bj), &
                  1.2_r8*rhof(i,k))
-              proc_rates%vtrmi(i,k)=proc_rates%vtrmi(i,k)*micro_mg_vtrmi_factor

-              fi(i,k)  = g*rho(i,k)*proc_rates%vtrmi(i,k)
               fni(i,k) = ifrac * fni(i,k) + &
                  (1._r8 - ifrac) * &
                  g*rho(i,k)* &
                  min(ajn(i,k)*gamma_bj_plus1/lami(i,k)**bj,1.2_r8*rhof(i,k))
            end if

+           ! apply micro_mg_vtrmi_factor once, after the small/large
+           ! blend, to both the mass- and number-weighted ice fall speeds
+           proc_rates%vtrmi(i,k) = proc_rates%vtrmi(i,k)*micro_mg_vtrmi_factor
+           fi(i,k)  = g*rho(i,k)*proc_rates%vtrmi(i,k)
+           fni(i,k) = fni(i,k)*micro_mg_vtrmi_factor
+
            ! Fix ice fall speed following IFS microphysics
            if (ifs_sed) then
               fi(i,k)=g*rho(i,k)*0.1_r8
--- a/micro_pumas_utils.F90
+++ b/micro_pumas_utils.F90
@@ -376,7 +376,8 @@

   ! Mean ice diameter can not grow bigger than twice the autoconversion
   ! threshold for snow.
-  ice_lambda_bounds(1) = 1._r8/2._r8*dcs
+  ! 1/(2*dcs), as in pumas v1.17; 1._r8/2._r8*dcs evaluates to dcs/2
+  ice_lambda_bounds(1) = 1._r8/(2._r8*dcs)
   ice_lambda_bounds(2) = 1._r8/1.e-6_r8

   mg_ice_props = MGHydrometeorProps(rhoi, dsph, &

Notes:

• vtrmi fix: outside the 18–20 µm blend band, fi and the VTRMI diagnostic stay bit-for-bit identical. The behaviour change comes from scaling fni.
• The ifs_sed override still sets both speeds afterwards, as before.
• The accre fix only touches the accre_sees_auto branch.
• All three fixes are answer-changing for CAM7 defaults.