Updated geometry to better separate the raised rim model and the ejec…

…ta blanket model to get better consistency
MintonGroup · Feb 9, 2022 · a99bd2a · a99bd2a
1 parent 584649c
commit a99bd2a
Show file tree

Hide file tree

Showing 5 changed files with 8 additions and 5 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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)