From 584cc82cd38f248d842052c086aed82141565c9e Mon Sep 17 00:00:00 2001
From: daminton <daminton@purdue.edu>
Date: Mon, 5 Dec 2016 21:29:26 +0000
Subject: [PATCH] Ejecta modified to be more self-consistent with empirical
 data on ejecta thickness (calibrated for strengthless soft rock)

---
 src/ejecta/ejecta_table_define.f90 | 76 +++++++++++++-----------------
 1 file changed, 32 insertions(+), 44 deletions(-)

diff --git a/src/ejecta/ejecta_table_define.f90 b/src/ejecta/ejecta_table_define.f90
index 73fccd2d..b6efb0a1 100644
--- a/src/ejecta/ejecta_table_define.f90
+++ b/src/ejecta/ejecta_table_define.f90
@@ -42,10 +42,14 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt)
    ! Executable code
 
    ! Get estimate of size of ejb table
-   crater%ejdis = DISEJB * 2.3_DP * crater%frad**(1.006_DP)  ! Continuous ejecta distance From Melosh (1989) eq. 6.3.1
+   if (.not.user%discontinuous) then
+      crater%ejdis =  2 * 2.3_DP * crater%frad**(1.006_DP)  ! Continuous ejecta distance From Melosh (1989) eq. 6.3.1
+   else
+      crater%ejdis = DISEJB * 2.3_DP * crater%frad**(1.006_DP)  ! Continuous ejecta distance From Melosh (1989) eq. 6.3.1
+   end if
                                                         ! We go out a factor of 3 to get the discontinuous ejecta thickness 
-   domain%ejbres = (crater%ejdis - crater%frad) / EJBTABSIZE
-   lrad = crater%frad
+   domain%ejbres = (crater%ejdis - crater%rad) / EJBTABSIZE
+   lrad = crater%frad 
    erad = crater%rad
    ejtble = EJBTABSIZE
    firstrun = .true.
@@ -60,53 +64,37 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt)
          dimp = crater%imp
          vimp = crater%impvel
       end if 
+   
 
       call regolith_melt_zone(user,crater,dimp,vimp,rmelt,depthb)
+   end if
  
-      !write(*,*) '     lrad/Df                   vej                       ebh                       melt & 
-      !           fraction              melt thickness'
-      do k = 0,EJBTABSIZE
-         call ejecta_rootfind(user,crater,domain,erad,lrad,vejsq,ejang,firstrun)
-         if (k >= 1) then
-            call ejecta_thickness(user,crater,eradold,erad,lrad - domain%ejbres,lrad,thick)
-            ejb(k)%lrad = log(lrad - 0.5_DP * domain%ejbres)
-            ejb(k)%thick = log(thick) 
-            ejb(k)%vesq = vejsq
-            ejb(k)%angle = ejang
-            ejb(k)%erad = erad
+   !write(*,*) '     lrad/Df                   vej                       ebh                       melt & 
+   !           fraction              melt thickness'
+   do k = 0,EJBTABSIZE
+      call ejecta_rootfind(user,crater,domain,erad,lrad,vejsq,ejang,firstrun)
+      if (k >= 1) then
+         call ejecta_thickness(user,crater,eradold,erad,lrad - domain%ejbres,lrad,thick)
+         ejb(k)%lrad = log(lrad - 0.5_DP * domain%ejbres)
+         ejb(k)%thick = log(thick) 
+         ejb(k)%vesq = vejsq
+         ejb(k)%angle = ejang
+         ejb(k)%erad = erad
+         if (present(melt)) then
             call regolith_melt_fraction(dimp,depthb,erad,eradold,rmelt,melt)
             ejb(k)%meltfrac = melt
-            !write(*,*) lrad/crater%frad,sqrt(vejsq),thick,melt,thick*melt
-            if ((thick <= VSMALL) .or. (abs(eradold - erad) < VSMALL)) then
-               ejtble = k
-               crater%ejdis = lrad
-               exit
-            end if
          end if
-         lrad = lrad + domain%ejbres
-         eradold = erad
-      end do
-      !write(*,*) 'A MELT ZONE of ',crater%frad,' meter-sized crater: ',rmelt,'at a rim',ejb(1)%meltfrac
-   else 
-     do k = 0,EJBTABSIZE
-        call ejecta_rootfind(user,crater,domain,erad,lrad,vejsq,ejang,firstrun)
-        if (k >= 1) then
-           call ejecta_thickness(user,crater,eradold,erad,lrad - domain%ejbres,lrad,thick)
-           ejb(k)%lrad = log(lrad - 0.5_DP * domain%ejbres)
-           ejb(k)%thick = log(thick) 
-           ejb(k)%vesq = vejsq
-           ejb(k)%angle = ejang
-           ejb(k)%erad = erad
-           if ((thick <= VSMALL) .or. (abs(eradold - erad) < VSMALL)) then
-              ejtble = k
-              crater%ejdis = lrad
-              exit
-           end if
-        end if
-        lrad = lrad + domain%ejbres
-        eradold = erad
-     end do
-   end if
+         !write(*,*) lrad/crater%rad,erad/crater%rad,sqrt(vejsq),thick !,melt,thick*melt
+         if ((thick <= VSMALL) .or. (abs(eradold - erad) < VSMALL)) then
+            ejtble = k
+            crater%ejdis = lrad
+            exit
+         end if
+      end if
+      lrad = lrad + domain%ejbres
+      eradold = erad
+   end do
+   !write(*,*) 'A MELT ZONE of ',crater%frad,' meter-sized crater: ',rmelt,'at a rim',ejb(1)%meltfrac
    ! Get pixel space distance
    crater%ejdispx = nint(crater%ejdis / user%pix)