fix: Trap floating point overflows

galacticusorg · Sep 18, 2023 · 6f74f91 · 6f74f91
1 parent 4260c63
commit 6f74f91
Showing 1 changed file with 9 additions and 7 deletions.
diff --git a/source/nodes.operators.physics.CGM.chemistry.F90 b/source/nodes.operators.physics.CGM.chemistry.F90
@@ -461,13 +461,15 @@ subroutine cgmChemistryComputeState(self,node,temperature,chemicalDensities,mass
     ! Truncate masses to zero to avoid unphysical behavior.
     call chemicalMasses%enforcePositive()
     massChemicals=chemicalMasses%sumOver()
-    if     (                                                     &
-         &   massChemicals > 0.0d0                               &
-         &  .and.                                                &
-         &   massChemicals > hotHalo%mass()                      &
-         & ) call chemicalMasses%scale(                          &
-         &                             +hotHalo%mass          () &
-         &                             /        massChemicals    &
+    if     (                                                                        &
+         &             massChemicals >         0.0d0                                &
+         &  .and.                                                                   &
+         &             massChemicals >         hotHalo%mass()                       &
+         &  .and.                                                                   &
+         &   -exponent(massChemicals)+exponent(hotHalo%mass()) < maxExponent(0.0d0) &
+         & ) call chemicalMasses%scale(                                             &
+         &                             +hotHalo%mass          ()                    &
+         &                             /        massChemicals                       &
          &                            )
     ! Scale all chemical masses by their mass in atomic mass units to get a number density.
     call chemicalMasses%massToNumber(chemicalDensities)