Modified ejecta emplacement so that it is more correctly handles asym…

…metric ejecta on slopes and rough surfaces

src/ejecta/ejecta_blanket.f90

@@ -93,7 +93,7 @@ subroutine ejecta_blanket(user,crater,domain,erad,lrad,vejsq,ejang,firstrun)
    ! lrad = erad + Psi * user%trad
 
    ! flat plane landing distance calculation
-   lrad = erad + vejsq*sin(2*ejang)/user%gaccel
+   lrad = erad + vejsq * sin(2 * ejang) / user%gaccel
 
    return
 end subroutine ejecta_blanket

src/ejecta/ejecta_emplace.f90

@@ -94,7 +94,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,ejbmass)
    real(DP),intent(out) :: ejbmass
 
    ! Internal variables
-   real(DP) :: lrad,lrad0,lradsq,cdepth,distance,erad,craterslope,baseline,inneredge,outeredge
+   real(DP) :: lrad,lrad0,lradsq,cdepth,distance,erad,craterslope,landslope,baseline,inneredge,outeredge,lradtrue
    integer(I4B),parameter :: MAXLOOP = 10 ! Maximum number of times to loop the ejecta angle correction calculation
    integer(I4B) :: xpi,ypi,i,j,k,n,inc,incsq,iradsq
    real(DP) :: xp,yp,fradsq,radsq,ebh,ejdissq,continuous
@@ -148,7 +148,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,ejbmass)
    if (crater%ejdis <= crater%rad) return
 
    ! determine area to effect
-   inc = max(min(crater%ejdispx + 1,PBCLIM*user%gridsize),1)
+   inc = max(min(nint(1.5_DP * crater%ejdis / user%pix) + 1,PBCLIM*user%gridsize),1)
    crater%maxinc = max(crater%maxinc,inc)
    radsq = crater%rad**2
    incsq = inc**2
@@ -286,14 +286,24 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,ejbmass)
 
          ! Find angle-corrected landing distance
          distance = sqrt((crater%xl - xp)**2 + (crater%yl - yp)**2)
-         baseline = crater%melev + ((i * crater%xslp) + (j * crater%yslp)) * user%pix - surf(xpi,ypi)%dem
+
+         baseline = crater%melev + ((i * crater%xslp) + (j * crater%yslp)) * user%pix 
          craterslope = atan(baseline / distance)
+
+
          lrad = distance
          lrad0 = lrad
          do n = 1,MAXLOOP
-            k = max(min(1 + int((lrad - inneredge) / (outeredge - inneredge) * (EJBTABSIZE - 1.0_DP)),ejtble),1)
+            !k = max(min(1 + int((lrad - inneredge) / (outeredge - inneredge) * (EJBTABSIZE - 1.0_DP)),ejtble),1)
             call ejecta_interpolate(crater,domain,lrad,ejb,ejtble,ebh,vsq=vsq,theta=ejtheta,erad=erad)
-            lrad =  distance * (erad + sqrt(vsq) * sin(2*ejtheta)) / (erad + sqrt(vsq) * (sin(2 * (ejtheta + craterslope))))
+            erad = exp(erad)
+            baseline = surf(xpi,ypi)%dem - erad * sin(craterslope)
+            landslope = atan(baseline / (distance - erad * cos(craterslope)))
+
+            lradtrue = erad * cos(craterslope) + vsq / (user%gaccel * cos(landslope)) * &
+                                             (sin(2 * ejtheta + 2 * craterslope - landslope) - sin(landslope))
+            lrad = lrad * distance / lradtrue
+            if (lrad > outeredge) exit
             if (abs(lrad - lrad0) < domain%small) exit
             lrad0 = lrad
          end do

src/ejecta/ejecta_table_define.f90

@@ -45,7 +45,7 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt)
    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
+      crater%ejdis = 2 * 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%rad) / EJBTABSIZE
@@ -79,7 +79,7 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt)
          ejb(k)%thick = log(thick) 
          ejb(k)%vesq = vejsq
          ejb(k)%angle = ejang
-         ejb(k)%erad = erad
+         ejb(k)%erad = log(erad)
          if (present(melt)) then
             call regolith_melt_fraction(dimp,depthb,erad,eradold,rmelt,melt)
             ejb(k)%meltfrac = melt
@@ -96,6 +96,7 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt)
    end do
    !write(*,*) 'A MELT ZONE of ',crater%frad,' meter-sized crater: ',rmelt,'at a rim',ejb(1)%meltfrac
    ! Get pixel space distance
+   crater%ejdis = crater%ejdis / 2.0_DP
    crater%ejdispx = nint(crater%ejdis / user%pix)
 
    return