Improved complex crater slope collapse routine to add in perceptible noise

src/Makefile.am

@@ -82,7 +82,6 @@ crater/crater_soften.f90\
 crater/crater_soften_accumulate.f90\
 crater/crater_subpixel_diffusion.f90\
 crater/crater_make_list.f90\
-crater/crater_critical_slope.f90\
 crater/crater_degradation_function.f90\
 crater/crater_visibility.f90\
 crater/crater_get_degradation_state.f90\

src/crater/crater_form_interior.f90

@@ -60,53 +60,16 @@ subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative ,newele
    ! Executable code
    r = sqrt(x_relative**2+y_relative**2) / crater%frad
 
-   !rh = crater%rimheight 
-   !fld = -crater%floordepth 
-   !flrad = 0.5_DP * crater%floordiam / crater%frad 
-
-
-   ! Use polynomial crater profile similar to that of Fassett et al. (2014), but the parameters are set by the crater dimensions
-   !c1 = (fld - rh) / (flrad + flrad**2 / 3._DP - flrad**3 / 6._DP - 7._DP / 6._DP)
-   !c0 = rh - (7._DP / 6._DP) * c1
-   !c2 = c1 / 3._DP
-   !c3 = -c2 / 2._DP
-!
-!   if (r < flrad) then
-!      cform = fld 
-!   else
-!      cform = c0 + c1 * r + c2 * r**2 + c3 * r**3 
-!   end if
-
-
-   !if (r < r_floor) then
-   !   cform = -simple_depth_diam  * crater%fcrat
-   !else if (r < r_rim) then
-   !   cform =  (inner_c0 + inner_c1 * r + inner_c2 * r**2 + inner_c3 * r**3) * crater%fcrat 
-   !else 
-   !   cform =  (outer_c0 + outer_c1 * r + outer_c2 * r**2 + outer_c3 * r**3) * crater%fcrat
-   !end if      
-
    cform = crater_profile(user,crater,r)
    newdem = newelev + cform 
 
-   !if (crater%fcrat > crater%cxtran * 2) then 
-   if (crater%morphtype == "COMPLEX") then
-      ! Make this crater complex
-
-      !if (r < r_rim) then
-         if (newdem < (crater%melev - crater%floordepth)) then ! Flatten out the bottom of the crater
-            do layer = 1,user%numlayers ! Remove all pre-existing craters from this current pixel
-               call util_remove_from_layer(surfi,layer)
-            end do
-            newdem = crater%melev - crater%floordepth
-
-
-         end if
-      !end if
-
+   if (newdem < (crater%melev - crater%floordepth)) then 
+      newdem = crater%melev - crater%floordepth ! Flatten out the bottom of the crater regardless of the local slope
+      do layer = 1,user%numlayers ! Remove all pre-existing craters from the flat floor
+         call util_remove_from_layer(surfi,layer)
+      end do
    end if
 
-
    newdem = min(newdem,surfi%dem) ! Only allow excavation, no deposition
    elchange  = newdem - surfi%dem
    deltaMi = elchange

src/crater/crater_realistic_topography.f90

@@ -44,6 +44,7 @@ subroutine crater_realistic_topography(user,surf,crater,domain,ejecta_dem)
    ! Internal variables
    real(DP) :: deltaMtot
    integer(I4B) :: inc
+   real(DP),parameter :: complex_collapse_slope = 0.33_DP  ! Complex craters and basins undergo an extra strong slope collapse
 
    interface
       subroutine complex_floor(user,surf,crater,deltaMtot)
@@ -108,18 +109,26 @@ end subroutine crater_realistic_slope_texture
    end interface
 
    ! Executable code
-   call complex_terrace(user,surf,crater,deltaMtot)
-   call complex_wall_texture(user,surf,crater,domain,deltaMtot)
-   call complex_floor(user,surf,crater,deltaMtot)
-   call complex_peak(user,surf,crater,deltaMtot)
+
+   if (crater%morphtype .eq. 'COMPLEX') then
+      call complex_terrace(user,surf,crater,deltaMtot)
+      call complex_wall_texture(user,surf,crater,domain,deltaMtot)
+      call complex_floor(user,surf,crater,deltaMtot)
+      call complex_peak(user,surf,crater,deltaMtot)
+
+   else if (crater%morphtype .eq. 'SIMPLE') then
+      call complex_floor(user,surf,crater,deltaMtot)
+   end if
 
    ! Retrieve the size of the ejecta dem and correct for indexing
    inc = (size(ejecta_dem,1) - 1) / 2
    call ejecta_texture(user,surf,crater,deltaMtot,inc,ejecta_dem)
 
+   if (crater%morphtype .eq. 'COMPLEX') then
+      ! Do a final pass of the slope collapse with a shallower slope than normal to smooth out all of the sharp edges
+      call crater_slope_collapse(user,surf,crater,domain,(complex_collapse_slope * user%pix)**2,deltaMtot)
+   end if
 
-   ! Do a final pass of the slope collapse with a shallower slope than normal to smooth out all of the sharp edges
-   call crater_slope_collapse(user,surf,crater,domain,(0.33_DP * user%pix)**2,deltaMtot)
 
    return
 end subroutine crater_realistic_topography
@@ -633,7 +642,7 @@ subroutine ejecta_texture(user,surf,crater,deltaMtot,inc,ejecta_dem)
             if (insplat) noise = noise + (1.0_DP + splatnoise - hprof) * ejecta_dem(i,j) * splatmag
          end do
 
-         noise_dem(i,j) = noise_dem(i,j) + noise * areafrac
+         noise_dem(i,j) = noise_dem(i,j) + noise * areafrac 
 
       end do
    end do
@@ -683,10 +692,10 @@ subroutine crater_realistic_slope_texture(user,critical,inc,critarray)
    real(DP), dimension(2) :: rn
 
    ! Topographic noise parameters
-   integer(I4B), parameter :: num_octaves  = 8   ! Number of Perlin noise octaves
+   integer(I4B), parameter :: num_octaves  = 4   ! Number of Perlin noise octaves
    integer(I4B), parameter :: offset = 4000 ! Scales the random xy-offset so that each crater's random noise is unique 
-   real(DP), parameter :: xy_noise_fac = 0.05_DP  ! Spatial "size" of noise features at the first octave
-   real(DP), parameter :: noise_height = 0.3e0_DP  ! Magnitude of noise features at the first octave
+   real(DP), parameter :: xy_noise_fac = 0.125_DP  ! Spatial "size" of noise features at the first octave
+   real(DP), parameter :: noise_height = 0.4e0_DP  ! Magnitude of noise features at the first octave
    real(DP), parameter :: freq = 2.0_DP     ! Spatial size scale factor multiplier at each octave level
    real(DP), parameter :: pers = 1.00_DP  ! The relative size scaling at each octave level
 

src/crater/module_crater.f90

@@ -266,16 +266,6 @@ subroutine crater_make_list(domain,prod,ntotcrat,production_list)
       end subroutine crater_make_list
    end interface
 
-   interface
-      function crater_critical_slope(user,crater,lrad) result(critical)
-      use module_globals
-      type(usertype),intent(in) :: user
-      type(cratertype),intent(in) :: crater
-      real(DP),intent(in) :: lrad
-      real(DP) :: critical
-      end function crater_critical_slope
-   end interface
-
    interface
       function crater_degradation_function(user,r) result(Kd)
       use module_globals