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/crater/crater_realistic_topography.f90 b/src/crater/crater_realistic_topography.f90 index 6a185a26..ed1b29c4 100644 --- a/src/crater/crater_realistic_topography.f90 +++ b/src/crater/crater_realistic_topography.f90 @@ -642,11 +642,11 @@ subroutine ejecta_texture(user,surf,crater,deltaMtot,inc,ejecta_dem) splat_stretch = 16.0_DP splatmag = 0.10_DP - open(unit=12,file='params.txt',status='old') - read(12,*) num_octaves - read(12,*) xy_noise_fac - read(12,*) noise_height - close(12) + ! open(unit=12,file='params.txt',status='old') + ! read(12,*) num_octaves + ! read(12,*) xy_noise_fac + ! read(12,*) noise_height + ! close(12) ! Get the ejecta mass ejbmass = sum(ejecta_dem) diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index b9c40d63..04dee3f2 100644 --- a/src/ejecta/ejecta_emplace.f90 +++ b/src/ejecta/ejecta_emplace.f90 @@ -128,8 +128,8 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r if (user%doregotrack) call regolith_melt_zone(user,crater,crater%imp,crater%impvel,rm,dm) - crater%vdepth = crater%ejrim + crater%floordepth - crater%vrim = crater%ejrim + crater%rimheight + crater%vdepth = crater%rimheight + crater%floordepth + crater%vrim = crater%rimheight if (crater%ejdis <= crater%ejrad) return @@ -261,7 +261,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r areafrac = (1.0_DP - util_area_intersection(crater%ejrad,xbar,ybar,user%pix)) ebh = areafrac * ejdistribution(idistorted,jdistorted) * ebh - cumulative_elchange(i,j) = areafrac * cumulative_elchange(i,j) + ebh + cumulative_elchange(i,j) = areafrac * cumulative_elchange(i,j) + ebh + crater_profile(user, crater, lrad) if (user%dosoftening) then ! Do extra diffusive degradation over ejecta region diff --git a/src/ejecta/ejecta_table_define.f90 b/src/ejecta/ejecta_table_define.f90 index 1a91af58..01ab7ce8 100644 --- a/src/ejecta/ejecta_table_define.f90 +++ b/src/ejecta/ejecta_table_define.f90 @@ -75,7 +75,7 @@ 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) end if 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)