Merge branch 'realistic'

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

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 

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

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)

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

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

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)