Fixed bugs that occur when ejecta thickness is near 0

MintonGroup · Feb 11, 2024 · 56054d9 · 56054d9
1 parent ab4551f
commit 56054d9
Show file tree

Hide file tree

Showing 3 changed files with 35 additions and 6 deletions.
diff --git 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
@@ -69,18 +69,46 @@ 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
       if (present(theta)) theta = 0._DP
       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
@@ -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))