From 56054d90c3920afcf8ebd3ab7135031ea6982869 Mon Sep 17 00:00:00 2001
From: David Minton <daminton@purdue.edu>
Date: Sun, 11 Feb 2024 14:11:47 -0500
Subject: [PATCH] Fixed bugs that occur when ejecta thickness is near 0

---
 src/ejecta/ejecta_emplace.f90     |  1 +
 src/ejecta/ejecta_interpolate.f90 | 38 +++++++++++++++++++++++++++----
 src/util/util_sort_layer.f90      |  2 +-
 3 files changed, 35 insertions(+), 6 deletions(-)

diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90
index b6ec3381..fda8c4b4 100644
--- a/src/ejecta/ejecta_emplace.f90
+++ b/src/ejecta/ejecta_emplace.f90
@@ -258,6 +258,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ
             else
                call ejecta_interpolate(crater,domain,distance,ejb,ejtble,ebh,vsq=vsq,theta=ejtheta,erad=erad)
             end if
+            if (ebh < VSMALL) exit
             if (n > 1) then
                if ((abs(ebh0 - ebh) / ebh0) < domain%small) exit
             endif
diff --git a/src/ejecta/ejecta_interpolate.f90 b/src/ejecta/ejecta_interpolate.f90
index 3463c08e..1b181ee2 100644
--- a/src/ejecta/ejecta_interpolate.f90
+++ b/src/ejecta/ejecta_interpolate.f90
@@ -69,7 +69,11 @@ subroutine ejecta_interpolate(crater,domain,lrad,ejb,ejtble,ebh,vsq,theta,erad,m
       if (present(melt)) melt = ejb(k)%meltfrac - ((ejb(k)%meltfrac - ejb(k+1)%meltfrac) * frac) 
       if (present(erad)) erad = ejb(k)%erad - ((ejb(k)%erad - ejb(k+1)%erad) * frac) 
    end if
-   ebh = exp(ebh) 
+   if (ebh < LOGVSMALL) then 
+      ebh = 0.0_DP
+   else
+      ebh = exp(ebh) 
+   end if
    if (lrad > crater%ejdis) then
       ebh = 0._DP
       if (present(vsq)) vsq = 0._DP
@@ -77,10 +81,34 @@ subroutine ejecta_interpolate(crater,domain,lrad,ejb,ejtble,ebh,vsq,theta,erad,m
       if (present(melt)) melt = 0._DP
       if (present(erad)) erad = 0._DP
    else
-      if (present(vsq)) vsq = exp(vsq)
-      if (present(theta)) theta = exp(theta)
-      if (present(melt)) melt = exp(melt)
-      if (present(erad)) erad = exp(erad)
+      if (present(vsq)) then 
+         if (vsq < LOGVSMALL) then
+            vsq = 0.0_DP
+         else 
+            vsq = exp(vsq)
+         end if
+      end if
+      if (present(theta)) then
+         if (theta < LOGVSMALL) then
+            theta = 0.0_DP
+         else
+            theta = exp(theta)
+         end if
+      end if
+      if (present(melt)) then
+         if (melt < LOGVSMALL) then
+            melt = 0.0_DP
+         else
+            melt = exp(melt)
+         end if
+      end if
+      if (present(erad)) then
+         if (erad < LOGVSMALL) then
+            erad = 0.0_DP
+         else
+            erad = exp(erad)
+         end if
+      end if
    end if
    
    return
diff --git a/src/util/util_sort_layer.f90 b/src/util/util_sort_layer.f90
index c3d055b1..65cd70d0 100644
--- a/src/util/util_sort_layer.f90
+++ b/src/util/util_sort_layer.f90
@@ -57,7 +57,7 @@ subroutine util_sort_layer(user,surf,crater)
             temptime = surf(mx,my)%timestamp(1:user%numlayers)
 
             ! Sort the layers by crater diameter
-            call util_mrgrnk(surf(mx,my)%diam(1:user%numlayers),ind)
+            call util_mrgrnk(tempdiam,ind)
 
             do k=1,user%numlayers
                surf(mx,my)%diam(k)=tempdiam(ind(k))