From cc2b0060283de414241b81898b009337b853cb54 Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 10 Apr 2020 15:09:42 -0400 Subject: [PATCH] Improved complex crater slope collapse routine to add in perceptible noise --- src/Makefile.am | 1 - src/crater/crater_critical_slope.f90 | 44 -------------------- src/crater/crater_form_interior.f90 | 47 +++------------------- src/crater/crater_realistic_topography.f90 | 29 ++++++++----- src/crater/module_crater.f90 | 10 ----- 5 files changed, 24 insertions(+), 107 deletions(-) delete mode 100644 src/crater/crater_critical_slope.f90 diff --git a/src/Makefile.am b/src/Makefile.am index 1e0e80b9..edd252e2 100755 --- a/src/Makefile.am +++ b/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\ diff --git a/src/crater/crater_critical_slope.f90 b/src/crater/crater_critical_slope.f90 deleted file mode 100644 index 4348617b..00000000 --- a/src/crater/crater_critical_slope.f90 +++ /dev/null @@ -1,44 +0,0 @@ -! Project : CTEM -! Language : Fortran 2003 -! -! Description : Forms the maximum slope for the initial crater morphology -! -! -! Input -! Arguments : -! -! Output -! Arguments : -! -! -! Notes : -! -!********************************************************************************************************************************** - - -function crater_critical_slope(user,crater,lrad) result(critical) - use module_globals - use module_util - use module_crater, EXCEPT_THIS_ONE => crater_critical_slope - implicit none - - ! Arguments - type(usertype),intent(in) :: user - type(cratertype),intent(in) :: crater - real(DP),intent(in) :: lrad - real(DP) :: critical - - ! Internal variables - real(DP) :: r - - r = lrad / crater%frad - if (r < 0.2_DP) then - critical = 0.0_DP - else if (r < 0.98_DP) then - critical = max(abs(0.228_DP + 2 * 0.083_DP * r - 3 * 0.039_DP * r**2),CRITSLP) - else - critical = max(abs(0.187_DP - 2 * 0.018_DP * r - 3 * 0.015_DP * r**2),CRITSLP) - end if - -end function crater_critical_slope - diff --git a/src/crater/crater_form_interior.f90 b/src/crater/crater_form_interior.f90 index 369ba7a9..955055e6 100644 --- a/src/crater/crater_form_interior.f90 +++ b/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 diff --git a/src/crater/crater_realistic_topography.f90 b/src/crater/crater_realistic_topography.f90 index 3919ea75..a3ee3261 100644 --- a/src/crater/crater_realistic_topography.f90 +++ b/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 diff --git a/src/crater/module_crater.f90 b/src/crater/module_crater.f90 index 466d42f5..822d5f5c 100644 --- a/src/crater/module_crater.f90 +++ b/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