Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history 
Added template code to generate realistic topography.
  • Loading branch information
@daminton
daminton committed Mar 12, 2020
1 parent 02ee261 commit 3b7b0ad
Show file tree
Hide file tree
Showing 9 changed files with 265 additions and 51 deletions.
5 changes: 3 additions & 2 deletions src/Makefile.am
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,9 @@ bin_PROGRAMS = CTEM
#AM_FCFLAGS = -O3 -p -g -openmp -debug all -traceback -CB -assume byterecl -m64 -heap-arrays -FR

#gfortran optimized flags
#AM_FCFLAGS = -O3 -fopenmp -ffree-form -g -fbounds-check -fbacktrace
AM_FCFLAGS = -O3 -fopenmp -ffree-form -g -fbounds-check -fbacktrace
#gfortran debug flags
AM_FCFLAGS = -O0 -g -fopenmp -fbounds-check -Wall -Warray-bounds -Warray-temporaries -Wimplicit-interface -ffree-form
#AM_FCFLAGS = -O0 -g -fopenmp -fbounds-check -Wall -Warray-bounds -Warray-temporaries -Wimplicit-interface -ffree-form

CTEM_SOURCES = globals/module_globals.f90\
util/module_util.f90\
Expand Down Expand Up @@ -71,6 +71,7 @@ crater/crater_generate.f90\
crater/crater_find_visible.f90\
crater/crater_averages.f90\
crater/crater_emplace.f90\
crater/crater_realistic_topography.f90\
crater/crater_form_interior.f90\
crater/crater_form_exterior.f90\
crater/crater_record.f90\
Expand Down
1 change: 1 addition & 0 deletions src/crater/crater_emplace.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,7 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot)

! Executable code


! determine area to effect
! First make the interior of the crater
inc = max(min(crater%fradpx,PBCLIM*user%gridsize),1) + 1
Expand Down
58 changes: 39 additions & 19 deletions src/crater/crater_form_interior.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,40 +33,60 @@ subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative ,newele
real(DP),intent(out) :: deltaMi

! Internal variables
real(DP) :: cform,newdem,elchange,pikeD,r
real(DP) :: cform,newdem,elchange,pikeD,r, Hpeak
integer(I4B) :: layer

! A list for poped data
! A list for popped data
type(regolisttype),pointer :: poppedlist


! Empirical crater shape parameters from Fassett et al. (2014)
real(DP),parameter :: r_floor = 0.2_DP
real(DP),parameter :: simple_depth_diam = 0.181_DP
real(DP),parameter :: r_rim = 0.98_DP
real(DP),parameter :: inner_c0 = -0.229_DP
real(DP),parameter :: inner_c1 = 0.228_DP
real(DP),parameter :: inner_c2 = 0.083_DP
real(DP),parameter :: inner_c3 = -0.039_DP

real(DP),parameter :: outer_c0 = 0.188_DP
real(DP),parameter :: outer_c1 = -0.187_DP
real(DP),parameter :: outer_c2 = 0.018_DP
real(DP),parameter :: outer_c3 = 0.015_DP


! Executable code

!change digital elevation map
r = sqrt(x_relative**2+y_relative**2) / crater%frad
! Use empirical crater form from Fassett et al. 2014
if (r < 0.2_DP) then
cform = -0.181_DP * crater%fcrat
else if (r < 0.98_DP) then
cform = (-0.229_DP + 0.228_DP * r + 0.083_DP * r**2 - 0.039_DP * r**3) * crater%fcrat
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 = (0.188_DP - 0.187_DP * r + 0.018_DP * r**2 + 0.015_DP * r**3) * crater%fcrat
cform = (outer_c0 + outer_c1 * r + outer_c2 * r**2 + outer_c3 * r**3) * crater%fcrat
end if
newdem = newelev + cform

pikeD = 1.044e3_DP * (crater%fcrat * 1e-3_DP)**(0.301_DP) ! Pike (1977)
if ((crater%fcrat > crater%cxtran * 2) .and. newdem < (crater%melev - pikeD)) then
newdem = crater%melev - pikeD ! Flatten out the bottom of the crater
!write(*,*) x_relative,y_relative,util_perlin_noise(x_relative,y_relative)*1e3_DP
!newdem = newdem + max(util_perlin_noise(x_relative/crater%fcrat,y_relative/crater%fcrat)*1e3_DP,0.0_DP)
!newdem = newdem + abs(util_perlin_noise(x_relative,y_relative)*1e3_DP)
end if
if (newdem < (crater%melev - user%deplimit)) then
newdem = crater%melev - user%deplimit ! 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
pikeD = min(1.044e3_DP * (crater%fcrat * 1e-3_DP)**(0.301_DP), user%deplimit) ! Pike (1977) depth/diameter ratio of complex craters
if (crater%fcrat > crater%cxtran * 2) then
! Make this crater complex

!if (r < r_rim) then
if (newdem < (crater%melev - pikeD)) 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 - pikeD


end if
!end if

end if


newdem = min(newdem,surfi%dem) ! Only allow excavation, no deposition
elchange = newdem - surfi%dem
deltaMi = elchange
Expand Down
4 changes: 4 additions & 0 deletions src/crater/crater_populate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -195,6 +195,10 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt

! Place crater onto the surface
call crater_emplace(user,surf,crater,domain,ejbmass)


if (user%dorealistic) call crater_realistic_topography(user,surf,crater,domain,ejbmass)



call ejecta_distance_estimate(user,crater,domain,crater%ejdis) ! Fast but imprecise estimate of the total ejecta distance
Expand Down
166 changes: 166 additions & 0 deletions src/crater/crater_realistic_topography.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,166 @@
!***** crater/crater_realistic_topography
! Name
! crater_realistic_topography
! SYNOPSIS
! This uses
! * module_globals
! * module_util
! * module_crater
!
! call crater_realistic_topography(user,surf,crater,domain,deltaMtot)
!
! DESCRIPTION
! Creates realistic topography for each fresh crater using Perlin noise.
!
!
! ARGUMENTS
! Input
! * user -- User input parameters
! * surf -- Surface grid
! * crater -- Crater dimension container
! * domain -- Simulation domain variable container
! * deltaMtot -- Total displaced mass used to calculate mass-conserving ejecta blanket
!
! Output
! * surf -- Outputs the new ejecta blanket onto the grid
! * crater -- May affects the value of the maximum affected distance
!
! Notes
!
!**********************************************************************************************************************************
subroutine crater_realistic_topography(user,surf,crater,domain,deltaMtot)
use module_globals
use module_util
use module_crater, EXCEPT_THIS_ONE => crater_realistic_topography
implicit none

! Arguments
type(usertype),intent(in) :: user
type(surftype),dimension(:,:),intent(inout) :: surf
type(cratertype),intent(inout) :: crater
type(domaintype),intent(in) :: domain
real(DP),intent(out) :: deltaMtot

! Internal variables
real(DP) :: lradsq,newelev, x_relative, y_relative
integer(I4B) :: xpi,ypi,i,j,inc,incsq,iradsq
real(DP) :: xp,yp,fradsq,deltaMi,rimheight
logical :: lastloop
real(DP), dimension(2) :: rn


! Topographic noise parameters
real(DP) :: xynoise, znoise
integer(I4B) :: octave
real(DP) :: noise


! Complex crater floor
integer(I4B), parameter :: complex_floor_num_octaves = 8 ! Number of Perlin noise octaves
integer(I4B), parameter :: complex_floor_offset = 10000 ! Scales the random xy-offset so that each crater's random noise is unique
real(DP), parameter :: complex_floor_xy_noise_fac = 2.0_DP ! Spatial "size" of noise features at the first octave
real(DP), parameter :: complex_floor_noise_height = 1e-3_DP ! Vertical height of noise features as a function of crater radius at the first octave
real(DP), parameter :: complex_floor_freq = 2.0_DP ! Spatial size scale factor multiplier at each octave level
real(DP), parameter :: complex_floor_pers = 0.80_DP ! The relative size scaling at each octave level

! Complex crater peak
integer(I4B), parameter :: complex_peak_num_octaves = 8 ! Number of Perlin noise octaves
integer(I4B), parameter :: complex_peak_offset = 1000 ! Scales the random xy-offset so that each crater's random noise is unique
real(DP), parameter :: complex_peak_rad = 0.20_DP ! Radial size of central peak relative to frad
real(DP), parameter :: complex_peak_xy_noise_fac = 16.0_DP ! Spatial "size" of noise features at the first octave
real(DP), parameter :: complex_peak_freq = 2.0_DP ! Spatial size scale factor multiplier at each octave level
real(DP), parameter :: complex_peak_pers = 0.80_DP ! The relative size scaling at each octave level

! Complex crater wall
integer(I4B), parameter :: complex_wall_num_octaves = 3 ! Number of Perlin noise octaves
integer(I4B), parameter :: complex_wall_offset = 2000 ! Scales the random xy-offset so that each crater's random noise is unique
real(DP), parameter :: complex_wall_xy_noise_fac = 2.0_DP ! Spatial "size" of noise features at the first octave
real(DP), parameter :: complex_wall_freq = 2.0_DP ! Spatial size scale factor multiplier at each octave level
real(DP), parameter :: complex_wall_pers = 0.80_DP ! The relative size scaling at each octave level
real(DP), parameter :: complex_wall_noise_height = 0.10_DP ! Vertical height of noise features as a function of crater radius at the first octave
real(DP), parameter :: complex_wall_slope = 25._DP * DEG2RAD
integer(I4B), parameter :: complex_wall_terracefac = 16 ! Spatial size factor for terraces relative to crater radius
integer(I4B) :: complex_wall_terrace_num


! Executable code

call random_number(rn)

! determine area to effect
! First make the interior of the crater
inc = max(min(crater%fradpx,PBCLIM*user%gridsize),1) + 1
crater%maxinc = max(crater%maxinc,inc)
fradsq = crater%frad**2
deltaMtot = 0.0_DP
incsq = inc**2

do j=-inc,inc ! Do the loop in pixel space
do i=-inc,inc
iradsq = i**2 + j**2
! find distance from crater center

! find elevation and grid point
newelev = crater%melev + ((i * crater%xslp) + (j * crater%yslp)) * user%pix
xpi = crater%xlpx + i
ypi = crater%ylpx + j

xp = xpi * user%pix
yp = ypi * user%pix

! periodic boundary conditions
call util_periodic(xpi,ypi,user%gridsize)
x_relative = (crater%xl - xp)
y_relative = (crater%yl - yp)

lradsq = x_relative**2 + y_relative**2

if (lradsq > crater%frad**2) cycle




deltaMtot = deltaMtot + deltaMi

end do
end do


! ! Make the floor topographically "noisy"
! noise = 0.0_DP
! do octave = 1, complex_floor_num_octaves
! xynoise = complex_floor_xy_noise_fac * complex_floor_freq ** octave / crater%fcrat
! znoise = complex_floor_noise_height * complex_floor_pers ** (octave - 1) * crater%fcrat
! noise = noise + util_perlin_noise(xynoise * x_relative + complex_floor_offset * rn(1), &
! xynoise * y_relative + complex_floor_offset * rn(2)) * znoise
! end do
! newdem = newdem + max(noise,0.0_DP)
!
! ! Add in "noisy" central peak
! if (r < complex_peak_rad) then
! Hpeak = 0.032e3_DP * (crater%fcrat * 1e-3_DP)**(0.900_DP) ! Pike (1977)
! noise = 0.0_DP
! do octave = 1, complex_peak_num_octaves
! xynoise = complex_peak_xy_noise_fac * complex_peak_freq ** octave / crater%fcrat
! znoise = Hpeak * (1.0_DP - r / complex_peak_rad) * complex_peak_pers ** (octave - 1)
! noise = noise + util_perlin_noise(xynoise * x_relative + complex_peak_offset * rn(1), &
! xynoise * y_relative + complex_peak_offset * rn(2)) * znoise
! end do
! newdem = newdem + max(noise,0.0_DP)
! end if
! else ! Terraced walls
! complex_wall_terrace_num = int(complex_wall_terracefac * r) + 1
! noise = 0.0_DP
! do octave = 1, complex_wall_num_octaves
! xynoise = complex_wall_xy_noise_fac * complex_wall_freq ** octave / crater%fcrat
! znoise = complex_wall_noise_height * complex_wall_pers ** (octave - 1) * crater%fcrat
! noise = noise + util_perlin_noise(xynoise * x_relative + complex_wall_terrace_num * complex_wall_offset * rn(1), &
! xynoise * y_relative + complex_wall_terrace_num * complex_wall_offset * rn(2)) &
! * znoise
! end do
! newdem = max(newdem + noise,crater%melev - pikeD) !/ cos(complex_wall_slope)


return
end subroutine crater_realistic_topography

12 changes: 12 additions & 0 deletions src/crater/module_crater.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -107,6 +107,18 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot)
end subroutine crater_emplace
end interface

interface
subroutine crater_realistic_topography(user,surf,crater,domain,deltaMtot)
use module_globals
implicit none
type(usertype),intent(in) :: user
type(surftype),dimension(:,:),intent(inout) :: surf
type(cratertype),intent(inout) :: crater
type(domaintype),intent(in) :: domain
real(DP),intent(out) :: deltaMtot
end subroutine crater_realistic_topography
end interface

interface
subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative, newelev,deltaMi)
use module_globals
Expand Down
37 changes: 20 additions & 17 deletions src/globals/module_globals.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,10 +150,10 @@ module module_globals
integer(I4B) :: gridsize ! Resolution
integer(I4B) :: numlayers ! Number of perched layers
real(DP) :: pix ! Pixel size (m)
real(DP) :: mu_b ! Crater scaling exponential constant (ignored for basins)
real(DP) :: kv_b ! Crater scaling linear constant
real(DP) :: mu_r ! Crater scaling exponential constant (ignored for basins)
real(DP) :: kv_r ! Crater scaling linear constant
real(DP) :: mu_b ! Crater scaling exponential constant (ignored for basins)
real(DP) :: kv_b ! Crater scaling linear constant
real(DP) :: mu_r ! Crater scaling exponential constant (ignored for basins)
real(DP) :: kv_r ! Crater scaling linear constant
integer(I4B) :: seed ! Random number generator seed (only used in non-IDL driven mode)
real(DP) :: trho_b ! Target bedrock density
real(DP) :: trho_r ! Target surface regolith layer density
Expand All @@ -165,23 +165,24 @@ module module_globals
real(DP) :: prho ! Projectile density
character(STRMAX) :: sfdfile ! Name of size distribution file
character(STRMAX) :: velfile ! Name of velocity distribution file
real(DP) :: seisk,cohaccel ! seismic keff, cohesion breaking acceleration
real(DP) :: seisk,cohaccel ! seismic keff, cohesion breaking acceleration

! Optional input variables
logical :: docollapse ! Set T to use the slope collapse model (turning off speeds up the code for testing)
logical :: doangle ! Set to F to only do vertical impacts, otherwise do range of angles (default is T)
logical :: doporosity ! Porosity on/off flg. Set to F to turn the model off. Default F.
real(DP) :: basinimp ! Impactor size to switch to multiring basin
real(DP) :: maxcrat ! fraction that maximum crater can be relative to grid
real(DP) :: deplimit ! complex crater depth limit
logical :: docollapse ! Set T to use the slope collapse model (turning off speeds up the code for testing)
logical :: doangle ! Set to F to only do vertical impacts, otherwise do range of angles (default is T)
logical :: doporosity ! Porosity on/off flg. Set to F to turn the model off. Default F.
real(DP) :: basinimp ! Impactor size to switch to multiring basin
real(DP) :: maxcrat ! fraction that maximum crater can be relative to grid
real(DP) :: deplimit ! complex crater depth limit
logical :: dorealistic ! Set to T to enable realistic crater morphology. Default is F.

! Seismic input variables
logical :: doseismic ! Set to T if you want to do the seismic shaking model
real(DP) :: seisq ! Seismic energy attenuation quality factor (Q)
real(DP) :: neff ! impact seismic energy efficiency factor
real(DP) :: tvel ! target P-wave (body wave) speed (m/s)
real(DP) :: tfrac ! mean free path for seismic wave scattering in medium
real(DP) :: regcoh ! target surface regolith layer cohesion
logical :: doseismic ! Set to T if you want to do the seismic shaking model
real(DP) :: seisq ! Seismic energy attenuation quality factor (Q)
real(DP) :: neff ! impact seismic energy efficiency factor
real(DP) :: tvel ! target P-wave (body wave) speed (m/s)
real(DP) :: tfrac ! mean free path for seismic wave scattering in medium
real(DP) :: regcoh ! target surface regolith layer cohesion

! Crater diffusion input parameters
real(DP) :: Kd1 ! Degradation function coefficient (from Minton et al. (2019))
Expand Down Expand Up @@ -230,6 +231,8 @@ module module_globals
real(DP) :: shadedminh ! Minimum height for shaded relief map (m)
real(DP) :: shadedmaxh ! Maximum height for shaded relief map (m)
character(STRMAX) :: sfdcompare ! Type of run: 0 for normal, 1 for statistical (domain is reset between intervals)


end type usertype

! Derived data type for the ejecta blanket table elements
Expand Down
7 changes: 7 additions & 0 deletions src/io/io_input.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,7 @@ subroutine io_input(infile,user)
user%rbreak = 0.1e3_DP
user%fe = 5.0_DP
user%ejecta_truncation = 5.0_DP
user%dorealistic = .false.

open(unit=LUN,file=infile,status="old",iostat=ierr)
if (ierr /= 0) then
Expand Down Expand Up @@ -466,6 +467,12 @@ subroutine io_input(infile,user)
token = line(ifirst:ilast)
read(token, *) user%fe

case ("DOREALISTIC")
ifirst = ilast + 1
call io_get_token(line, ilength, ifirst, ilast, ierr)
token = line(ifirst:ilast)
read(token, *) user%dorealistic

!**************************************************************************
! The following is for backwards compatibility with older style input files
!**************************************************************************
Expand Down
Loading

0 comments on commit 3b7b0ad

Please sign in to comment.