diff --git a/src/ejecta/ejecta_blanket.f90 b/src/ejecta/ejecta_blanket.f90 index 232e1ebc..8362c0bb 100644 --- a/src/ejecta/ejecta_blanket.f90 +++ b/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 diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index adab74b6..b1a0422c 100644 --- a/src/ejecta/ejecta_emplace.f90 +++ b/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 diff --git a/src/ejecta/ejecta_table_define.f90 b/src/ejecta/ejecta_table_define.f90 index b6efb0a1..920cbb5a 100644 --- a/src/ejecta/ejecta_table_define.f90 +++ b/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