From 61e86da5453e109500528478d70cdf720aab1fc5 Mon Sep 17 00:00:00 2001
From: Paul Romano <romano7@gmail.com>
Date: Mon, 13 Feb 2012 16:57:59 -0500
Subject: [PATCH] Added checks for divide-by-zero on histogram interpolation.

---
 src/physics.F90 | 40 ++++++++++++++++++++++++++++++----------
 1 file changed, 30 insertions(+), 10 deletions(-)

diff --git a/src/physics.F90 b/src/physics.F90
index cce1abda14..a898ec849a 100644
--- a/src/physics.F90
+++ b/src/physics.F90
@@ -1145,7 +1145,7 @@ contains
        xi = prn()
        lc = lc + 2
        c_k = rxn % adist % data(lc + 2*NP + 1)
-       do k = 1, NP-1
+       do k = 1, NP - 1
           c_k1 = rxn % adist % data(lc + 2*NP + k+1)
           if (xi < c_k1) exit
           c_k = c_k1
@@ -1155,7 +1155,11 @@ contains
        mu0 = rxn % adist % data(lc + k)
        if (interp == HISTOGRAM) then
           ! Histogram interpolation
-          mu = mu0 + (xi - c_k)/p0
+          if (p0 > ZERO) then
+             mu = mu0 + (xi - c_k)/p0
+          else
+             mu = mu0
+          end if
 
        elseif (interp == LINEAR_LINEAR) then
           ! Linear-linear interpolation
@@ -1462,7 +1466,7 @@ contains
        r1 = prn()
        lc = lc + 2 ! start of EOUT
        c_k = edist % data(lc + 2*NP + 1)
-       do k = 1, NP-1
+       do k = 1, NP - 1
           c_k1 = edist % data(lc + 2*NP + k+1)
           if (r1 < c_k1) exit
           c_k = c_k1
@@ -1472,7 +1476,11 @@ contains
        p_l_k = edist % data(lc+NP+k)
        if (INTT == HISTOGRAM) then
           ! Histogram interpolation
-          E_out = E_l_k + (r1 - c_k)/p_l_k
+          if (p_l_k > ZERO) then
+             E_out = E_l_k + (r1 - c_k)/p_l_k
+          else
+             E_out = E_l_k
+          end if
 
        elseif (INTT == LINEAR_LINEAR) then
           ! Linear-linear interpolation -- not sure how you come about the
@@ -1689,7 +1697,7 @@ contains
        r1 = prn()
        lc = lc + 2 ! start of EOUT
        c_k = edist % data(lc + 2*NP + 1)
-       do k = 1, NP-1
+       do k = 1, NP - 1
           c_k1 = edist % data(lc + 2*NP + k+1)
           if (r1 < c_k1) exit
           c_k = c_k1
@@ -1699,7 +1707,11 @@ contains
        p_l_k = edist % data(lc+NP+k)
        if (INTT == HISTOGRAM) then
           ! Histogram interpolation
-          E_out = E_l_k + (r1 - c_k)/p_l_k
+          if (p_l_k > ZERO) then
+             E_out = E_l_k + (r1 - c_k)/p_l_k
+          else
+             E_out = E_l_k
+          end if
 
           ! Determine Kalbach-Mann parameters
           KM_R = edist % data(lc + 3*NP + k)
@@ -1832,7 +1844,7 @@ contains
        r1 = prn()
        lc = lc + 2 ! start of EOUT
        c_k = edist % data(lc + 2*NP + 1)
-       do k = 1, NP-1
+       do k = 1, NP - 1
           c_k1 = edist % data(lc + 2*NP + k+1)
           if (r1 < c_k1) exit
           c_k = c_k1
@@ -1842,7 +1854,11 @@ contains
        p_l_k = edist % data(lc+NP+k)
        if (INTT == HISTOGRAM) then
           ! Histogram interpolation
-          E_out = E_l_k + (r1 - c_k)/p_l_k
+          if (p_l_k > ZERO) then
+             E_out = E_l_k + (r1 - c_k)/p_l_k
+          else
+             E_out = E_l_k
+          end if
 
        elseif (INTT == LINEAR_LINEAR) then
           ! Linear-linear interpolation -- not sure how you come about the
@@ -1887,7 +1903,7 @@ contains
        r3 = prn()
        lc = lc + 2
        c_k = edist % data(lc + 2*NP + 1)
-       do k = 1, NP-1
+       do k = 1, NP - 1
           c_k1 = edist % data(lc + 2*NP + k+1)
           if (r3 < c_k1) exit
           c_k = c_k1
@@ -1897,7 +1913,11 @@ contains
        mu_k = edist % data(lc + k)
        if (JJ == HISTOGRAM) then
           ! Histogram interpolation
-          mu_out = mu_k + (r3 - c_k)/p_k
+          if (p_k > ZERO) then
+             mu_out = mu_k + (r3 - c_k)/p_k
+          else
+             mu_out = mu_k
+          end if
 
        elseif (JJ == LINEAR_LINEAR) then
           ! Linear-linear interpolation -- not sure how you come about the