Update production calculations

P_CaCO3%prod(k) and P_SiO2%prod(k) were not being computed correctly in cases
where there were multiple calcifiers and silicifiers, respectively. POP%prod(k)
was missing the zoo_graze_poc term, which is non-zero when zooplankton grazes
on non-zooplankton.

The default CESM configuration only has one calcifier and one silicifier, and
the single zooplankton does not graze on itself, so this should not be answer
changing when PFT_defaults = "CESM2".

src/marbl_mod.F90

@@ -3996,7 +3996,7 @@ subroutine marbl_compute_large_detritus_prod(k, auto_cnt, zoo_cnt,          &
     !  large detritus P
     !-----------------------------------------------------------------------
 
-    POP%prod(k) = Qp_zoo * sum(zoo_loss_poc(:)) + sum(remaining_P_pop(:))
+    POP%prod(k) = Qp_zoo * (sum(zoo_loss_poc(:)) + sum(zoo_graze_poc(:))) + sum(remaining_P_pop(:))
 
     if (POC%prod(k) <= c0) POP%prod(k) = c0
 
@@ -4005,22 +4005,25 @@ subroutine marbl_compute_large_detritus_prod(k, auto_cnt, zoo_cnt,          &
     !  33% of CaCO3 is remin when phyto are grazed
     !-----------------------------------------------------------------------
 
+    P_CaCO3%prod(k) = c0
     do auto_ind = 1, auto_cnt
        if (marbl_tracer_indices%auto_inds(auto_ind)%CaCO3_ind > 0) then
-          P_CaCO3%prod(k) = ((c1 - f_graze_CaCO3_remin) * auto_graze(auto_ind) + &
-               auto_loss(auto_ind) + auto_agg(auto_ind)) * QCaCO3(auto_ind)
-          P_CaCO3_ALT_CO2%prod(k) = P_CaCO3%prod(k)
+          P_CaCO3%prod(k) = P_CaCO3%prod(k) + ((c1 - f_graze_CaCO3_remin) * auto_graze(auto_ind) &
+               + auto_loss(auto_ind) + auto_agg(auto_ind)) * QCaCO3(auto_ind)
        endif
     end do
+    P_CaCO3_ALT_CO2%prod(k) = P_CaCO3%prod(k)
 
     !-----------------------------------------------------------------------
     !  large detritus SiO2
     !  grazed diatom SiO2, 60% is remineralized
     !-----------------------------------------------------------------------
 
+    P_SiO2%prod(k) = c0
     do auto_ind = 1, auto_cnt
        if (marbl_tracer_indices%auto_inds(auto_ind)%Si_ind > 0) then
-          P_SiO2%prod(k) = Qsi(auto_ind) * ((c1 - f_graze_si_remin) * auto_graze(auto_ind) + auto_agg(auto_ind) &
+          P_SiO2%prod(k) = P_SiO2%prod(k) + Qsi(auto_ind) &
+               * ((c1 - f_graze_si_remin) * auto_graze(auto_ind) + auto_agg(auto_ind) &
                + autotrophs(auto_ind)%loss_poc * auto_loss(auto_ind))
        endif
     end do