Update logic for when to compute o2_production

There is a possibility that photoC > 0 but NO3_V = NH4_V = 0 due to underflow
in the computation of those two values.

src/marbl_interior_tendency_mod.F90

@@ -3452,6 +3452,7 @@ subroutine compute_local_tendencies(km, marbl_tracer_indices, autotroph_derived_
  !-----------------------------------------------------------------------
  integer :: k, auto_ind, zoo_ind, n
  real(r8) :: auto_sum
+ real(r8) :: o2_production_denom
  !-----------------------------------------------------------------------
 
  associate( &
@@ -3691,15 +3692,17 @@ subroutine compute_local_tendencies(km, marbl_tracer_indices, autotroph_derived_
  o2_production(k) = c0
  do auto_ind = 1, autotroph_cnt
  if (.not. autotroph_settings(auto_ind)%Nfixer) then
- if (photoC(auto_ind,k) > c0) then
+ o2_production_denom = NO3_V(auto_ind,k) + NH4_V(auto_ind,k)
+ if (o2_production_denom > c0) then
  o2_production(k) = o2_production(k) + photoC(auto_ind,k) &
  * ((NO3_V(auto_ind,k) / (NO3_V(auto_ind,k) + NH4_V(auto_ind,k))) &
  / parm_Red_D_C_O2 &
  + (NH4_V(auto_ind,k) / (NO3_V(auto_ind,k) + NH4_V(auto_ind,k))) &
  / parm_Remin_D_C_O2)
  end if
  else
- if (photoC(auto_ind,k) > c0) then
+ o2_production_denom = NO3_V(auto_ind,k) + NH4_V(auto_ind,k) + Nfix(auto_ind,k)
+ if (o2_production_denom > c0) then
  o2_production(k) = o2_production(k) + photoC(auto_ind,k) &
  * ((NO3_V(auto_ind,k) / (NO3_V(auto_ind,k) + NH4_V(auto_ind,k) + Nfix(auto_ind,k))) &
  / parm_Red_D_C_O2 &