diff --git a/src/crater/crater_emplace.f90 b/src/crater/crater_emplace.f90 index a05c5938..850321f6 100644 --- a/src/crater/crater_emplace.f90 +++ b/src/crater/crater_emplace.f90 @@ -101,7 +101,7 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot) lradsq = x_relative**2 + y_relative**2 if (lradsq > crater%frad**2) cycle - call crater_form_interior(user,surf(xpi,ypi),crater,x_relative, y_relative,newelev,deltaMi) + call crater_form_interior(user,surf(xpi,ypi),crater,x_relative,y_relative,newelev,deltaMi) deltaMtot = deltaMtot + deltaMi ! do porosity computation if (user%doporosity) diff --git a/src/crater/crater_form_exterior_func.f90 b/src/crater/crater_form_exterior_func.f90 index 46a94dbe..31bdc48c 100644 --- a/src/crater/crater_form_exterior_func.f90 +++ b/src/crater/crater_form_exterior_func.f90 @@ -39,7 +39,7 @@ function crater_form_exterior_func(user,surf,crater,domain,rd,deltaMtot,lastloop type(surftype) :: surfi ! Now make the exterior of the crater - inc = int(crater%frad/user%pix*(domain%small/crater%rheight)**(-1._DP/RIMDROP)) ! Maximum distance of crater form + inc = int(crater%frad/user%pix*(domain%small/crater%rimheight)**(-1._DP/RIMDROP)) ! Maximum distance of crater form inc = max(min(max(crater%rimdispx,inc),PBCLIM*user%gridsize),1) crater%maxinc = max(crater%maxinc,inc) diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index 969ecd0e..bf8ef50f 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -164,7 +164,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt ! Stamp the current time onto the crater crater%timestamp = real(curyear + real(icrater,kind=DP) / real(ntotcrat,kind=DP) * user%interval,kind=SP) ! Find the visible crater parameters - call crater_find_visible(user,crater,domain) + call crater_dimensions(user,crater,domain) ! Crater is big enough to keep, so record it into the true distribution ntrue = ntrue + 1 diff --git a/src/crater/crater_realistic_topography.f90 b/src/crater/crater_realistic_topography.f90 index 9fd12825..b24f55b0 100644 --- a/src/crater/crater_realistic_topography.f90 +++ b/src/crater/crater_realistic_topography.f90 @@ -39,52 +39,51 @@ subroutine crater_realistic_topography(user,surf,crater,domain,deltaMtot) type(surftype),dimension(:,:),intent(inout) :: surf type(cratertype),intent(inout) :: crater type(domaintype),intent(in) :: domain - real(DP),intent(out) :: deltaMtot + real(DP),intent(inout) :: deltaMtot ! Internal variables - real(DP) :: lradsq,newelev, x_relative, y_relative + real(DP) :: lradsq,x_relative, y_relative integer(I4B) :: xpi,ypi,i,j,inc,incsq,iradsq - real(DP) :: xp,yp,fradsq,deltaMi,rimheight + real(DP) :: r,xp,yp,fradsq,deltaMi,rimheight logical :: lastloop real(DP), dimension(2) :: rn + interface + subroutine complex_floor(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi) + use module_globals + implicit none + type(usertype),intent(in) :: user + type(surftype),intent(inout) :: surfi + type(cratertype),intent(in) :: crater + real(DP),intent(in) :: r,x_relative, y_relative + real(DP),dimension(:),intent(in) :: rn + real(DP),intent(out) :: deltaMi + end subroutine complex_floor - ! 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 + subroutine complex_peak(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi) + use module_globals + implicit none + type(usertype),intent(in) :: user + type(surftype),intent(inout) :: surfi + type(cratertype),intent(in) :: crater + real(DP),intent(in) :: r,x_relative, y_relative + real(DP),dimension(:),intent(in) :: rn + real(DP),intent(out) :: deltaMi + end subroutine complex_peak + subroutine complex_wall(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi) + use module_globals + implicit none + type(usertype),intent(in) :: user + type(surftype),intent(inout) :: surfi + type(cratertype),intent(in) :: crater + real(DP),intent(in) :: r,x_relative, y_relative + real(DP),dimension(:),intent(in) :: rn + real(DP),intent(out) :: deltaMi + end subroutine complex_wall + end interface ! Executable code - call random_number(rn) ! determine area to effect @@ -92,7 +91,6 @@ subroutine crater_realistic_topography(user,surf,crater,domain,deltaMtot) 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 @@ -101,7 +99,6 @@ subroutine crater_realistic_topography(user,surf,crater,domain,deltaMtot) ! 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 @@ -116,9 +113,14 @@ subroutine crater_realistic_topography(user,surf,crater,domain,deltaMtot) lradsq = x_relative**2 + y_relative**2 if (lradsq > crater%frad**2) cycle + + r = sqrt(x_relative**2+y_relative**2) / crater%frad - - + if (crater%fcrat > crater%cxtran * 2) then + call complex_floor(user,surf(xpi,ypi),crater,r,x_relative,y_relative,rn,deltaMi) + call complex_peak (user,surf(xpi,ypi),crater,r,x_relative,y_relative,rn,deltaMi) + call complex_wall (user,surf(xpi,ypi),crater,r,x_relative,y_relative,rn,deltaMi) + end if deltaMtot = deltaMtot + deltaMi @@ -126,41 +128,163 @@ subroutine crater_realistic_topography(user,surf,crater,domain,deltaMtot) 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 +subroutine complex_peak(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi) + 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),intent(inout) :: surfi + type(cratertype),intent(in) :: crater + real(DP),intent(in) :: r,x_relative, y_relative + real(DP),dimension(:),intent(in) :: rn + real(DP),intent(out) :: deltaMi + + ! Complex crater peak + real(DP), parameter :: complex_peak_rad = 0.20_DP ! Radial size of central peak relative to frad + 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_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 + + ! Internal variables + real(DP) :: newdem,elchange,pikeD,Hpeak + + ! Topographic noise parameters + real(DP) :: xynoise, znoise + integer(I4B) :: octave + real(DP) :: noise + + if (r > complex_peak_rad) return + ! Add in "noisy" central peak + newdem = surfi%dem + 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) + + elchange = newdem - surfi%dem + deltaMi = elchange + surfi%dem = newdem + +end subroutine complex_peak + +subroutine complex_floor(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi) + 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),intent(inout) :: surfi + type(cratertype),intent(in) :: crater + real(DP),intent(in) :: r,x_relative, y_relative + real(DP),dimension(:),intent(in) :: rn + real(DP),intent(out) :: deltaMi + + ! Complex crater floor parameters + 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 + + ! Internal variables + real(DP) :: newdem,elchange,pikeD,Hpeak + + ! Topographic noise parameters + real(DP) :: xynoise, znoise + integer(I4B) :: octave + real(DP) :: noise + + newdem = surfi%dem + ! 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) + + elchange = newdem - surfi%dem + deltaMi = elchange + surfi%dem = newdem + +end subroutine complex_floor + + +subroutine complex_wall(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi) + 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),intent(inout) :: surfi + type(cratertype),intent(in) :: crater + real(DP),intent(in) :: r,x_relative, y_relative + real(DP),dimension(:),intent(in) :: rn + real(DP),intent(out) :: deltaMi + + ! Internal variables + real(DP) :: newdem,elchange,pikeD,pikeFloor,Hpeak + + ! Topographic noise parameters + real(DP) :: xynoise, znoise + integer(I4B) :: octave + real(DP) :: noise + + ! 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 = 32 ! Spatial size factor for terraces relative to crater radius + integer(I4B) :: complex_wall_terrace_num + + + newdem = surfi%dem + pikeD = min(1.044e3_DP * (crater%fcrat * 1e-3_DP)**(0.301_DP), user%deplimit) ! Pike (1977) depth/diameter ratio of complex craters + !pikeFloor = + pikeFloor = 0.7_DP + if (r < pikeFloor) return + 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) + + + elchange = newdem - surfi%dem + deltaMi = elchange + surfi%dem = newdem + +end subroutine complex_wall + + + + diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index a9df7e6c..c4fa0167 100644 --- a/src/ejecta/ejecta_emplace.f90 +++ b/src/ejecta/ejecta_emplace.f90 @@ -94,7 +94,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) real(DP),intent(in) :: deltaMtot ! Internal variables - real(DP) :: lrad,lradsq,cdepth + real(DP) :: lrad,lradsq integer(I4B),parameter :: MAXLOOP = 100 ! Maximum number of times to loop the ejecta angle correction calculation integer(I4B) :: xpi,ypi,i,j,k,n,inc,incsq,iradsq,idistorted,jdistorted real(DP) :: xp,yp,fradsq,fradpxsq,radsq,ebh,ejdissq,ejbmass,fmasscons,areafrac,xbar,ybar,krad,kdiffmax @@ -121,9 +121,8 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) !call regolith_melt_zone(user,crater,crater%imp,crater%impvel,rm,dm) - cdepth = DDRATIO * crater%fcrat - crater%vdepth = crater%ejrim + cdepth - crater%vrim = crater%ejrim + crater%rheight + crater%vdepth = crater%ejrim + crater%floordepth + crater%vrim = crater%ejrim + crater%rimheight if (crater%ejdis <= crater%rad) return