diff --git a/src/crater/crater_dimensions.f90 b/src/crater/crater_dimensions.f90 index f1f97f96..9a45e04f 100644 --- a/src/crater/crater_dimensions.f90 +++ b/src/crater/crater_dimensions.f90 @@ -69,6 +69,7 @@ subroutine crater_dimensions(user,crater,domain) ! Calculate the radius where the inner wall meets the original pre-existing surface ! This is used to demark the location where excavation transitions to deposition crater%ejrad = max(crater_profile_find_r_inner_wall(user,crater) * crater%frad, crater%rad) + crater%ejrim = 0.14_DP * (crater%fcrat * 1e-3 * 0.5_DP)**(0.74_DP) * 1e3_DP ! McGetchin et al. (1973) Thickness of ejecta at rim !find rim for counting purposes crater%frim = RIMFAC * crater%frad diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index a0530cc0..3db1fa06 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -251,10 +251,10 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt if (crater%ejdis > domain%smallest_ejecta) then ! Estimated size is big enough, so proceed with precise calculation if (user%doregotrack) then call ejecta_table_define(user,crater,domain,ejb,ejtble,melt) - call ejecta_interpolate(crater,domain,crater%frad,ejb(1:ejtble),ejtble,crater%ejrim) + !call ejecta_interpolate(crater,domain,crater%frad,ejb(1:ejtble),ejtble,crater%ejrim) else call ejecta_table_define(user,crater,domain,ejb,ejtble) - call ejecta_interpolate(crater,domain,crater%frad,ejb(1:ejtble),ejtble,crater%ejrim) + !call ejecta_interpolate(crater,domain,crater%frad,ejb(1:ejtble),ejtble,crater%ejrim) end if call ejecta_emplace(user,surf,crater,domain,ejb(1:ejtble),ejtble,ejbmass,age,age_resolution,ejecta_dem) else @@ -346,6 +346,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt hmin = minval(surf(:,:)%dem) if (any(surf(:,:)%dem /= surf(:,:)%dem)) then write(*,*) 'Invalid surface elevation detected. Halting.' + write(*,*) crater%imp, crater%impvel, crater%xl, crater%yl, crater%sinimpang exit end if end do ! end crater production loop diff --git a/src/crater/crater_profile.f90 b/src/crater/crater_profile.f90 index 3723d081..d597d9f1 100644 --- a/src/crater/crater_profile.f90 +++ b/src/crater/crater_profile.f90 @@ -48,7 +48,7 @@ function crater_profile(user,crater,r) result(h) if (r < flrad) then h = -crater%floordepth else if (r >= 1.0_DP) then - h = crater%rimheight * (r**(-RIMDROP)) + h = (crater%rimheight - crater%ejrim) * (r**(-RIMDROP)) else h = c0 + c1 * r + c2 * r**2 + c3 * r**3 end if diff --git a/src/ejecta/ejecta_table_define.f90 b/src/ejecta/ejecta_table_define.f90 index 1a91af58..91eb7b14 100644 --- a/src/ejecta/ejecta_table_define.f90 +++ b/src/ejecta/ejecta_table_define.f90 @@ -75,9 +75,9 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt) ! This will be replaced r = lrad / crater%frad if (lrad >= crater%frad) then - thick = crater%rimheight * r**(-3.0_DP) + thick = crater_profile(user, crater, r) + crater%ejrim * r**(-EJPROFILE) else - thick = max(crater_profile(user,crater,r),VSMALL) + thick = max(crater_profile(user,crater,r) + crater%ejrim,VSMALL) end if ejb(k)%thick = log(thick) ejb(k)%vesq = vejsq diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index 2d7a70e5..6157dc9d 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -310,6 +310,7 @@ module module_globals real(DP),parameter :: CT = KT * (PI*THIRD)**(SIXTH) real(DP),parameter :: RIMDROP = 4.20_DP ! Power law index for rim profile +real(DP), parameter :: EJPROFILE = 3.0_DP ! Power law index for ejecta profile real(DP),parameter :: TRSIM = 1.25_DP ! ? real(DP),parameter :: EXFAC = 0.1_DP ! Excavation depth relative to transient crater diameter real(DP),parameter :: CXEXPS = 1._DP / 0.885_DP - 1.0_DP ! Complex crater scaling exponent (see Croft 1985)