diff --git a/src/crater/crater_emplace.f90 b/src/crater/crater_emplace.f90 index 850321f6..b67c0983 100644 --- a/src/crater/crater_emplace.f90 +++ b/src/crater/crater_emplace.f90 @@ -79,6 +79,9 @@ subroutine crater_emplace(user,surf,crater,domain,deltaMtot) fradsq = crater%frad**2 deltaMtot = 0.0_DP !ejbmass incsq = inc**2 + + ! TEMP UNTIL BETTER FORMULA + crater%ejrad = crater%frad ! This loop may not be parallelizable because of the linked list operation inside crater_form_interior do j=-inc,inc ! Do the loop in pixel space do i=-inc,inc diff --git a/src/crater/crater_form_interior.f90 b/src/crater/crater_form_interior.f90 index 3f5b76bc..a1ea3572 100644 --- a/src/crater/crater_form_interior.f90 +++ b/src/crater/crater_form_interior.f90 @@ -27,7 +27,7 @@ subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative ,newele ! Arguments type(usertype),intent(in) :: user type(surftype),intent(inout) :: surfi - type(cratertype),intent(in) :: crater + type(cratertype),intent(inout) :: crater real(DP),intent(in) :: x_relative, y_relative real(DP),intent(in) :: newelev real(DP),intent(out) :: deltaMi @@ -66,7 +66,6 @@ subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative ,newele ! 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 @@ -77,6 +76,8 @@ subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative ,newele else cform = c0 + c1 * r + c2 * r**2 + c3 * r**3 end if + ! TEMP UNTIL I WRITE THE SOLUTION PROPERLY + if (cform > 0.0_DP .and. r * crater%frad < crater%ejrad) crater%ejrad = r * crater%frad !if (r < r_floor) then ! cform = -simple_depth_diam * crater%fcrat diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index d1e2f03b..fca76d77 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -196,11 +196,6 @@ 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 ! For very steep size distributions, only a fraction of the ! craters are retained. The full ejecta_table_define function @@ -220,6 +215,8 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt else ejtble = 0 end if + + if (user%dorealistic) call crater_realistic_topography(user,surf,crater,domain,ejbmass) ! Collapse any remaining unstable slopes if (user%docollapse) call crater_slope_collapse(user,surf,crater,domain,ejbmass) diff --git a/src/crater/crater_realistic_topography.f90 b/src/crater/crater_realistic_topography.f90 index 5c4ef4e6..b42f3fdd 100644 --- a/src/crater/crater_realistic_topography.f90 +++ b/src/crater/crater_realistic_topography.f90 @@ -42,96 +42,210 @@ subroutine crater_realistic_topography(user,surf,crater,domain,deltaMtot) real(DP),intent(inout) :: deltaMtot ! Internal variables - real(DP) :: lradsq,x_relative, y_relative + real(DP) :: lradsq,xbar, ybar integer(I4B) :: xpi,ypi,i,j,inc,incsq,iradsq - real(DP) :: r,xp,yp,fradsq,deltaMi,rimheight + real(DP) :: r,xp,yp,deltaMi,rimheight,rad logical :: lastloop real(DP), dimension(2) :: rn interface - subroutine complex_floor(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi) + subroutine complex_floor(user,surf,crater,rn,deltaMtot) 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 + type(surftype),dimension(:,:),intent(inout) :: surf + type(cratertype),intent(inout) :: crater real(DP),dimension(:),intent(in) :: rn - real(DP),intent(out) :: deltaMi + real(DP),intent(inout) :: deltaMtot end subroutine complex_floor - subroutine complex_peak(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi) + subroutine complex_peak(user,surf,crater,rn,deltaMtot) 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 + type(surftype),dimension(:,:),intent(inout) :: surf + type(cratertype),intent(inout) :: crater real(DP),dimension(:),intent(in) :: rn - real(DP),intent(out) :: deltaMi + real(DP),intent(inout) :: deltaMtot end subroutine complex_peak - subroutine complex_wall(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi) + subroutine complex_wall_texture(user,surf,crater,rn,deltaMtot) 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 + type(surftype),dimension(:,:),intent(inout) :: surf + type(cratertype),intent(inout) :: crater real(DP),dimension(:),intent(in) :: rn - real(DP),intent(out) :: deltaMi - end subroutine complex_wall + real(DP),intent(inout) :: deltaMtot + + end subroutine complex_wall_texture + + subroutine complex_terrace(user,surf,crater,rn,deltaMtot) + use module_globals + implicit none + type(usertype),intent(in) :: user + type(surftype),dimension(:,:),intent(inout) :: surf + type(cratertype),intent(inout) :: crater + real(DP),dimension(:),intent(in) :: rn + real(DP),intent(inout) :: deltaMtot + end subroutine complex_terrace + + subroutine complex_rim(user,surf,crater,rn,deltaMtot) + use module_globals + implicit none + type(usertype),intent(in) :: user + type(surftype),dimension(:,:),intent(inout) :: surf + type(cratertype),intent(inout) :: crater + real(DP),dimension(:),intent(in) :: rn + real(DP),intent(inout) :: deltaMtot + end subroutine complex_rim end interface ! Executable code call random_number(rn) + + call complex_rim(user,surf,crater,rn,deltaMtot) + call complex_terrace(user,surf,crater,rn,deltaMtot) + call complex_wall_texture(user,surf,crater,rn,deltaMtot) + call complex_floor(user,surf,crater,rn,deltaMtot) + call complex_peak(user,surf,crater,rn,deltaMtot) - ! determine area to effect - ! First make the interior of the crater - inc = max(min(crater%fradpx,PBCLIM*user%gridsize),1) + 1 + return +end subroutine crater_realistic_topography + +subroutine complex_peak(user,surf,crater,rn,deltaMtot) + use module_globals + use module_util + use module_crater, EXCEPT_THIS_ONE => crater_realistic_topography + implicit none + type(usertype),intent(in) :: user + type(surftype),dimension(:,:),intent(inout) :: surf + type(cratertype),intent(inout) :: crater + real(DP),dimension(:),intent(in) :: rn + real(DP),intent(inout) :: deltaMtot + + ! Internal variables + real(DP) :: newdem,elchange,xbar,ybar,areafrac,r + integer(I4B) :: i,j,inc,xpi,ypi + + ! Complex crater peak + real(DP), parameter :: rad = 0.20_DP ! Radial size of central peak relative to frad + integer(I4B), parameter :: num_octaves = 8 ! Number of Perlin noise octaves + integer(I4B), parameter :: offset = 1000 ! Scales the random xy-offset so that each crater's random noise is unique + real(DP), parameter :: xy_noise_fac = 32.0_DP ! Spatial "size" 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 = 0.80_DP ! The relative size scaling at each octave level + + ! Topographic noise parameters + real(DP) :: xynoise, znoise + integer(I4B) :: octave + real(DP) :: noise + + + inc = max(min(nint(rad * crater%frad / user%pix),PBCLIM*user%gridsize),1) + 1 crater%maxinc = max(crater%maxinc,inc) - fradsq = crater%frad**2 - 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 + + do j = -inc,inc + do i = -inc,inc 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 - - 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 + newdem = surf(xpi,ypi)%dem + + xbar = xpi * user%pix - crater%xl + ybar = ypi * user%pix - crater%yl + + areafrac = util_area_intersection(rad * crater%frad,xbar,ybar,user%pix) + + r = sqrt(xbar**2 + ybar**2) / crater%frad - deltaMtot = deltaMtot + deltaMi + noise = 0.0_DP + do octave = 1, num_octaves + xynoise = xy_noise_fac * freq ** (octave - 1) / crater%fcrat + znoise = crater%peakheight * (1.0_DP - r / rad) * pers ** (octave - 1) + noise = noise + util_perlin_noise(xynoise * xbar + offset * rn(1), & + xynoise * ybar + offset * rn(2)) * znoise + end do + newdem = newdem + max(noise,0.0_DP) * areafrac + elchange = newdem - surf(xpi,ypi)%dem + deltaMtot = deltaMtot + elchange + surf(xpi,ypi)%dem = newdem end do end do + return + +end subroutine complex_peak + +subroutine complex_floor(user,surf,crater,rn,deltaMtot) + use module_globals + use module_util + use module_crater, EXCEPT_THIS_ONE => crater_realistic_topography + implicit none + type(usertype),intent(in) :: user + type(surftype),dimension(:,:),intent(inout) :: surf + type(cratertype),intent(inout) :: crater + real(DP),dimension(:),intent(in) :: rn + real(DP),intent(inout) :: deltaMtot + + ! Internal variables + real(DP) :: newdem,elchange,xbar,ybar,areafrac + integer(I4B) :: i,j,inc,xpi,ypi + + ! Complex crater floor parameters + integer(I4B), parameter :: num_octaves = 8 ! Number of Perlin noise octaves + integer(I4B), parameter :: offset = 10000 ! Scales the random xy-offset so that each crater's random noise is unique + real(DP), parameter :: xy_noise_fac = 4.0_DP ! Spatial "size" of noise features at the first octave + real(DP), parameter :: noise_height = 1e-3_DP ! Vertical height of noise features as a function of crater radius at the first octave + real(DP), parameter :: freq = 2.0_DP ! Spatial size scale factor multiplier at each octave level + real(DP), parameter :: pers = 0.80_DP ! The relative size scaling at each octave level + + + ! Topographic noise parameters + real(DP) :: xynoise, znoise + integer(I4B) :: octave + real(DP) :: noise + inc = max(min(nint(0.5_DP * crater%floordiam / user%pix),PBCLIM*user%gridsize),1) + 1 + crater%maxinc = max(crater%maxinc,inc) + + do j = -inc,inc + do i = -inc,inc + xpi = crater%xlpx + i + ypi = crater%ylpx + j + + ! periodic boundary conditions + call util_periodic(xpi,ypi,user%gridsize) + newdem = surf(xpi,ypi)%dem + + xbar = xpi * user%pix - crater%xl + ybar = ypi * user%pix - crater%yl + + areafrac = util_area_intersection(0.5_DP * crater%floordiam,xbar,ybar,user%pix) + + ! Make the floor topographically "noisy" + noise = 0.0_DP + do octave = 1, num_octaves + xynoise = xy_noise_fac * freq ** (octave - 1) / crater%fcrat + znoise = noise_height * pers ** (octave - 1) * crater%fcrat + noise = noise + util_perlin_noise(xynoise * xbar + offset * rn(1), & + xynoise * ybar + offset * rn(2)) * znoise + end do + newdem = newdem + max(noise,0.0_DP) * areafrac + + elchange = newdem - surf(xpi,ypi)%dem + deltaMtot = deltaMtot + elchange + surf(xpi,ypi)%dem = newdem + end do + end do return -end subroutine crater_realistic_topography +end subroutine complex_floor + -subroutine complex_peak(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi) +subroutine complex_wall_texture(user,surf,crater,rn,deltaMtot) use module_globals use module_util use module_crater, EXCEPT_THIS_ONE => crater_realistic_topography @@ -139,47 +253,73 @@ subroutine complex_peak(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi) ! 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 + type(surftype),dimension(:,:),intent(inout) :: surf + type(cratertype),intent(inout) :: crater 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 + real(DP),intent(inout) :: deltaMtot ! Internal variables - real(DP) :: newdem,elchange + real(DP) :: newdem,elchange,xbar,ybar,areafrac,flr,r + integer(I4B) :: i,j,inc,xpi,ypi ! Topographic noise parameters real(DP) :: xynoise, znoise integer(I4B) :: octave - real(DP) :: noise + real(DP) :: noise,dnoise - if (r > complex_peak_rad) return - ! Add in "noisy" central peak - newdem = surfi%dem - noise = 0.0_DP - do octave = 1, complex_peak_num_octaves - xynoise = complex_peak_xy_noise_fac * complex_peak_freq ** octave / crater%fcrat - znoise = crater%peakheight * (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 + ! Complex crater all texture parameters + real(DP), parameter :: rad = 0.20_DP ! Radial size of central peak relative to frad + integer(I4B), parameter :: num_octaves = 8 ! 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 = 32.0_DP ! Spatial "size" of noise features at the first octave + real(DP), parameter :: noise_height = 1.0e-3_DP ! Spatial "size" 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 = 0.80_DP ! The relative size scaling at each octave level -end subroutine complex_peak -subroutine complex_floor(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi) + inc = max(min(nint(1.5_DP * crater%frad / user%pix),PBCLIM*user%gridsize),1) + 1 + crater%maxinc = max(crater%maxinc,inc) + + flr = crater%floordiam / crater%fcrat + + do j = -inc,inc + do i = -inc,inc + xpi = crater%xlpx + i + ypi = crater%ylpx + j + + ! periodic boundary conditions + call util_periodic(xpi,ypi,user%gridsize) + newdem = surf(xpi,ypi)%dem + + xbar = xpi * user%pix - crater%xl + ybar = ypi * user%pix - crater%yl + + r = sqrt(xbar**2 + ybar**2) / crater%frad + areafrac = 1.0 - util_area_intersection(0.4_DP * crater%floordiam,xbar,ybar,user%pix) + areafrac = areafrac * util_area_intersection(1.5_DP * crater%frad,xbar,ybar,user%pix) + areafrac = areafrac / max(r**2,0.1_DP) + + ! Add some roughness to the walls + noise = 0.0_DP + do octave = 1, num_octaves + xynoise = xy_noise_fac * freq ** (octave - 1) / crater%fcrat + znoise = noise_height * pers ** (octave - 1) * crater%fcrat + noise = noise + util_perlin_noise(xynoise * xbar + offset * rn(1), & + xynoise * ybar + offset * rn(2))* znoise + end do + newdem = newdem + noise * areafrac + + elchange = newdem - surf(xpi,ypi)%dem + deltaMtot = deltaMtot + elchange + surf(xpi,ypi)%dem = newdem + end do + end do + + return +end subroutine complex_wall_texture + +subroutine complex_rim(user,surf,crater,rn,deltaMtot) + !Makes terraced walls by applying noisy topographic diffusion in discrete radial zones along the wall use module_globals use module_util use module_crater, EXCEPT_THIS_ONE => crater_realistic_topography @@ -187,47 +327,74 @@ subroutine complex_floor(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi) ! 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 + type(surftype),dimension(:,:),intent(inout) :: surf + type(cratertype),intent(inout) :: crater 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 + real(DP),intent(inout) :: deltaMtot ! Internal variables - real(DP) :: newdem,elchange + real(DP) :: newdem,elchange,rad,xbar,ybar + integer(I4B) :: xpi,ypi,i,j,inc,maxhits ! Topographic noise parameters real(DP) :: xynoise, znoise - integer(I4B) :: octave - real(DP) :: noise + real(DP) :: noise,dnoise + + ! Complex crater wall + integer(I4B) :: num_octaves ! Number of Perlin noise octaves + integer(I4B) :: offset ! Scales the random xy-offset so that each crater's random noise is unique + real(DP) :: xy_noise_fac ! Spatial "size" of noise features at the first octave + real(DP) :: freq ! Spatial size scale factor multiplier at each octave level + real(DP) :: pers ! The relative size scaling at each octave level + real(DP) :: noise_height ! Vertical height of noise features as a function of crater radius at the first octave + + ! determine area to effect + ! First make the interior of the crater + rad = 1.25_DP * crater%frad + + ! Create diffusion noise for terraces + num_octaves = 4 + offset = 3000 + xy_noise_fac = 5.00_DP + noise_height = 1.00_DP + + inc = max(min(nint(rad / user%pix),PBCLIM*user%gridsize),1) + 1 + crater%maxinc = max(crater%maxinc,inc) + + do j = -inc,inc + do i = -inc,inc + xpi = crater%xlpx + i + ypi = crater%ylpx + j - 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 + ! periodic boundary conditions + call util_periodic(xpi,ypi,user%gridsize) + newdem = surf(xpi,ypi)%dem + + xbar = xpi * user%pix - crater%xl + ybar = ypi * user%pix - crater%yl + + ! Make scalloped rim + znoise = noise_height * crater%fcrat + xynoise = xy_noise_fac / crater%fcrat + dnoise = util_perlin_noise(xynoise * xbar + offset * rn(1), & + xynoise * ybar + offset * rn(2)) * znoise + noise = sqrt(sqrt(dnoise**2)) + !newdem = min(crater%melev + crater%ejrim - noise , newdem) + if (crater%melev + crater%ejrim - noise < newdem) newdem = newdem - 1._DP * crater%ejrim + + elchange = newdem - surf(xpi,ypi)%dem + deltaMtot = deltaMtot + elchange + surf(xpi,ypi)%dem = newdem + end do end do - newdem = newdem + max(noise,0.0_DP) - elchange = newdem - surfi%dem - deltaMi = elchange - surfi%dem = newdem + return -end subroutine complex_floor +end subroutine complex_rim -subroutine complex_wall(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi) +subroutine complex_terrace(user,surf,crater,rn,deltaMtot) + !Makes terraced walls by applying noisy topographic diffusion in discrete radial zones along the wall use module_globals use module_util use module_crater, EXCEPT_THIS_ONE => crater_realistic_topography @@ -235,51 +402,138 @@ subroutine complex_wall(user,surfi,crater,r,x_relative,y_relative,rn,deltaMi) ! 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 + type(surftype),dimension(:,:),intent(inout) :: surf + type(cratertype),intent(inout) :: crater real(DP),dimension(:),intent(in) :: rn - real(DP),intent(out) :: deltaMi + real(DP),intent(inout) :: deltaMtot ! Internal variables - real(DP) :: newdem,elchange + real(DP) :: newdem,elchange,rad,xbar,ybar,r,flr + integer(I4B) :: xpi,ypi,i,j,inc,maxhits ! Topographic noise parameters real(DP) :: xynoise, znoise integer(I4B) :: octave - real(DP) :: noise + real(DP) :: noise,dnoise ! Complex crater wall - integer(I4B), parameter :: complex_wall_num_octaves = 4 ! 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 = 1.5_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.70_DP ! The relative size scaling at each octave level - real(DP), parameter :: complex_wall_noise_height = 0.02_DP ! Vertical height of noise features as a function of crater radius at the first octave - integer(I4B), parameter :: complex_wall_terracefac = 8 ! Spatial size factor for terraces relative to crater radius - integer(I4B) :: complex_wall_terrace_num - - - newdem = surfi%dem - if (r < crater%floordiam / crater%fcrat) 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 + integer(I4B) :: num_octaves ! Number of Perlin noise octaves + integer(I4B) :: offset ! Scales the random xy-offset so that each crater's random noise is unique + real(DP) :: xy_noise_fac ! Spatial "size" of noise features at the first octave + real(DP) :: freq ! Spatial size scale factor multiplier at each octave level + real(DP) :: pers ! The relative size scaling at each octave level + real(DP) :: noise_height ! Vertical height of noise features as a function of crater radius at the first octave + real(DP) :: scal_height ! Vertical height of noise features as a function of crater radius at the first octave + real(DP) :: diff_height ! Magnitude of diffusion noise + integer(I4B) :: terracefac ! Spatial size factor for terraces relative to crater radius + integer(I4B) :: terrace_num ! The size of each terrace will be crater radius / terrace_num + + ! Internal variables + real(DP),dimension(:,:),allocatable :: kdiff,cumulative_elchange,areafrac + integer(I4B),dimension(:,:,:),allocatable :: indarray + + + ! determine area to effect + ! First make the interior of the crater + rad = 1.25_DP * crater%frad + + + ! Create diffusion noise for terraces + num_octaves = 4 + xy_noise_fac = 4.00_DP + diff_height = 8e-5_DP * (crater%fcrat)**2 + noise_height = 0.02_DP + freq = 2.0_DP + pers = 0.80_DP + terracefac = 8 + flr = crater%floordiam / crater%fcrat + + inc = max(min(nint(rad / user%pix),PBCLIM*user%gridsize),1) + 1 + crater%maxinc = max(crater%maxinc,inc) + + allocate(cumulative_elchange(-inc:inc,-inc:inc)) + allocate(kdiff(-inc:inc,-inc:inc)) + allocate(areafrac(-inc:inc,-inc:inc)) + allocate(indarray(2,-inc:inc,-inc:inc)) + + cumulative_elchange = 0.0_DP + kdiff = 0.0_DP + indarray = inc - 1 + areafrac = 0.0_DP + + do j = -inc,inc + do i = -inc,inc + xpi = crater%xlpx + i + ypi = crater%ylpx + j + + ! periodic boundary conditions + call util_periodic(xpi,ypi,user%gridsize) + newdem = surf(xpi,ypi)%dem + + indarray(1,i,j) = xpi + indarray(2,i,j) = ypi + + ! Cut a hole out from the floor and also outside the rim of the crater + xbar = xpi * user%pix - crater%xl + ybar = ypi * user%pix - crater%yl + + r = sqrt(xbar**2 + ybar**2) / crater%frad + + areafrac(i,j) = (1.0_DP - util_area_intersection(0.95_DP * crater%floordiam / 2,xbar,ybar,user%pix)) + areafrac(i,j) = areafrac(i,j) * util_area_intersection(2.0_DP*crater%frad,xbar,ybar,user%pix) + + areafrac(i,j) = areafrac(i,j) / max(r**4,1._DP) + + ! Make the terraces + ! This shifts the Perlin noise pattern at discrete radial distances, simulating terraces + terrace_num = int(terracefac * ((r - 1._DP) / (flr - 1._DP))**(1.5_DP)) + noise = 0.0_DP + if ((r <= 1.0_DP).and.(r >= flr)) then + offset = 2000 + xynoise = xy_noise_fac / crater%fcrat + znoise = noise_height * crater%fcrat + noise = noise + util_perlin_noise(xynoise * xbar + terrace_num * offset * rn(1), & + xynoise * ybar + terrace_num * offset * rn(2))* znoise + newdem = max(newdem + noise,crater%melev - crater%floordepth) + end if + elchange = newdem - surf(xpi,ypi)%dem + deltaMtot = deltaMtot + elchange + surf(xpi,ypi)%dem = newdem + + ! Now smooth out the sharp edges of the terraces with noisy diffusion + noise = 0.0_DP + if (r <= 1.0_DP) then + do octave = 1, num_octaves + xynoise = xy_noise_fac * freq ** (octave - 1) / crater%fcrat + znoise = diff_height * pers ** (octave - 1) + noise = noise + util_perlin_noise(xynoise * xbar + offset * rn(1), & + xynoise * ybar + offset * rn(2)) * znoise + end do + end if + kdiff(i,j) = noise + + end do end do - newdem = max(newdem + noise,crater%melev - crater%floordepth) + kdiff(:,:) = kdiff(:,:) - minval(kdiff(:,:)) + ! Cut out the holes + kdiff(:,:) = kdiff(:,:) * areafrac(:,:) !* 0.0_DP - elchange = newdem - surfi%dem - deltaMi = elchange - surfi%dem = newdem + maxhits = 1 + call util_diffusion_solver(user,surf,2 * inc + 1,indarray,kdiff,cumulative_elchange,maxhits) -end subroutine complex_wall + do j = -inc,inc + do i = -inc,inc + xpi = indarray(1,i,j) + ypi = indarray(2,i,j) + surf(xpi,ypi)%dem = surf(xpi,ypi)%dem + cumulative_elchange(i,j) + surf(xpi,ypi)%ejcov = max(surf(xpi,ypi)%ejcov + cumulative_elchange(i,j),0.0_DP) + end do + end do + deallocate(cumulative_elchange,kdiff,indarray,areafrac) + return +end subroutine complex_terrace diff --git a/src/crater/module_crater.f90 b/src/crater/module_crater.f90 index 7c3adc3f..a573bf97 100644 --- a/src/crater/module_crater.f90 +++ b/src/crater/module_crater.f90 @@ -125,7 +125,7 @@ subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative, newele implicit none type(usertype),intent(in) :: user type(surftype),intent(inout) :: surfi - type(cratertype),intent(in) :: crater + type(cratertype),intent(inout) :: crater real(DP),intent(in) :: x_relative, y_relative real(DP),intent(in) :: newelev real(DP),intent(out) :: deltaMi diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index c4fa0167..078f2ac5 100644 --- a/src/ejecta/ejecta_emplace.f90 +++ b/src/ejecta/ejecta_emplace.f90 @@ -90,7 +90,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) type(cratertype),intent(inout) :: crater type(domaintype),intent(in) :: domain integer(I4B),intent(in) :: ejtble - type(ejbtype),dimension(ejtble),intent(in) :: ejb + type(ejbtype),dimension(ejtble),intent(inout) :: ejb real(DP),intent(in) :: deltaMtot ! Internal variables @@ -124,7 +124,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) crater%vdepth = crater%ejrim + crater%floordepth crater%vrim = crater%ejrim + crater%rimheight - if (crater%ejdis <= crater%rad) return + if (crater%ejdis <= crater%ejrad) return ! determine area to effect inc = max(min(nint(min(PI * user%trad / user%pix, min(crater%ejdis, crater%frad * user%ejecta_truncation / user%pix))) + 1, & @@ -138,13 +138,10 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) inc = max(inc,dradsq) dradsq = dradsq**2 - if (user%dosoftening) then - !kdiffmax = user%Kd1 * crater%frad**(user%psi) - kdiffmax = crater_degradation_function(user,crater%frad) - end if + if (user%dosoftening) kdiffmax = crater_degradation_function(user,crater%frad) crater%maxinc = max(crater%maxinc,inc) - radsq = crater%rad**2 + radsq = crater%ejrad**2 fradsq = crater%frad**2 fradpxsq = crater%fradpx**2 ejdissq = crater%ejdis**2 @@ -206,7 +203,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) ! This must be done iteratively because the ejection distance and ejection angle vary distance = lrad maxslp = -huge(maxslp) - klo = int((log(lrad) - log(crater%rad)) / domain%ejbres) + klo = int((log(lrad) - log(crater%ejrad)) / domain%ejbres) do n = 1,MAXLOOP call ejecta_interpolate(crater,domain,distance,ejb,ejtble,ebh,vsq=vsq,theta=ejtheta,erad=erad,melt=melt) if ((n > 1).and.((abs(ebh0 - ebh) / ebh0) < domain%small)) exit @@ -245,13 +242,13 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) if (abs(jdistorted) > inc) cycle iradsq = idistorted**2 + jdistorted**2 - if ((iradsq > incsq).or.(distance <= crater%rad)) cycle + if ((iradsq > incsq).or.(distance <= crater%ejrad)) cycle ! we need to cut a hole out from the inside of the crater xbar = xpi * user%pix - crater%xl ybar = ypi * user%pix - crater%yl - areafrac = (1.0_DP - util_area_intersection(crater%rad,xbar,ybar,user%pix)) + 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 @@ -288,6 +285,8 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) ! Do mass conservation by adjusting ejecta thickness fmasscons = (-deltaMtot)/ ejbmass cumulative_elchange = cumulative_elchange * fmasscons + crater%ejrim = crater%ejrim * fmasscons + ejb(:)%thick = ejb(:)%thick * fmasscons maxhits = 1 ! Create box for soften calculation (will be no bigger than the grid itself) if (2 * inc + 1 < user%gridsize) then diff --git a/src/ejecta/ejecta_interpolate.f90 b/src/ejecta/ejecta_interpolate.f90 index a8ba7931..c1086650 100644 --- a/src/ejecta/ejecta_interpolate.f90 +++ b/src/ejecta/ejecta_interpolate.f90 @@ -39,8 +39,8 @@ subroutine ejecta_interpolate(crater,domain,lrad,ejb,ejtble,ebh,vsq,theta,erad,m ! Executable code ! Locate ourselves in the table - inneredge = crater%rad - outeredge = crater%rad * exp(domain%ejbres * EJBTABSIZE) + inneredge = crater%ejrad + outeredge = crater%ejrad * exp(domain%ejbres * EJBTABSIZE) k = max(min(1 + int((log(lrad) - log(inneredge)) / (log(outeredge) - log(inneredge)) * (EJBTABSIZE - 1.0_DP)),ejtble),1) loglrad = log(lrad) logtablerad = ejb(k)%lrad diff --git a/src/ejecta/ejecta_rootfind.f90 b/src/ejecta/ejecta_rootfind.f90 index 716da6cf..2bd48e45 100644 --- a/src/ejecta/ejecta_rootfind.f90 +++ b/src/ejecta/ejecta_rootfind.f90 @@ -78,7 +78,7 @@ subroutine ejecta_rootfind(user,crater,domain,erad,lrad,vejsq,ejang,firstrun) f2=ejecta_blanket_func(user,crater,domain,x2,lrad,vejsq,ejang,firstrun) end if end do - if ((erad<=crater%rad).and.(erad>0._DP)) exit + if ((erad<=crater%ejrad).and.(erad>0._DP)) exit factor=0.5_DP*(factor+1._DP) end do bracket @@ -153,7 +153,7 @@ subroutine ejecta_rootfind(user,crater,domain,erad,lrad,vejsq,ejang,firstrun) end if niter = i erad = resultat - if ((erad <= crater%rad).and.(lrad >= crater%rad).and.(erad >= 0._DP)) exit + if ((erad <= crater%ejrad).and.(lrad >= crater%ejrad).and.(erad >= 0._DP)) exit factor = 0.5_DP * (factor + 1.0_DP) ! Failed. Try again with a new factor end do everything diff --git a/src/ejecta/ejecta_table_define.f90 b/src/ejecta/ejecta_table_define.f90 index bbb732e3..5d5dd1ca 100644 --- a/src/ejecta/ejecta_table_define.f90 +++ b/src/ejecta/ejecta_table_define.f90 @@ -39,15 +39,26 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt) ! Regotrack internal variables real(DP) :: rmelt,depthb,dimp,vimp + ! This will be replaced by its own function + real(DP) :: c0,c1,c2,c3,flrad,rh,fld,r + rh = crater%rimheight + fld = -crater%floordepth + flrad = 0.5_DP * crater%floordiam / crater%frad + 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 + !^^^^^^^^^^^^^^^ + ! Executable code ! Get estimate of size of ejb table crater%continuous = RCONT * crater%frad**(EXPCONT) ! Continuous ejecta distance From Moore (1974) eq. 1 crater%ejdis = DISEJB * crater%continuous ! We go out a factor of 3 to get the discontinuous ejecta thickness - domain%ejbres = (log(crater%ejdis) - log(crater%rad)) / EJBTABSIZE - lrad = crater%rad !exp(log(crater%rad) !+ domain%ejbres) - erad = crater%rad + domain%ejbres = (log(crater%ejdis) - log(crater%ejrad)) / EJBTABSIZE + lrad = crater%ejrad !exp(log(crater%rad) !+ domain%ejbres) + erad = crater%ejrad ejtble = EJBTABSIZE firstrun = .true. thick = 0._DP @@ -73,11 +84,14 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt) if (k >= 1) then !call ejecta_thickness(user,crater,eradold,erad,lrad - domain%ejbres,lrad,thick) ejb(k)%lrad = log(lrad) - ! Use McGetchin et al. 1973 for ejecta thickness + + ! This will be replaced + r = lrad / crater%frad if (lrad >= crater%frad) then - thick = 0.14_DP * crater%frad**(0.74_DP) * (lrad / crater%frad)**(-3.0_DP) + thick = crater%rimheight * r**(-3.0_DP) else - thick = 0.14_DP * crater%frad**(0.74_DP) / (crater%frad - crater%rad) * (lrad - crater%rad) + !thick = 0.14_DP * crater%frad**(0.74_DP) / (crater%frad - crater%ejrad) * (lrad - crater%ejrad) + thick = c0 + c1 * r + c2 * r**2 + c3 * r**3 end if ejb(k)%thick = log(thick) ejb(k)%vesq = vejsq diff --git a/src/ejecta/module_ejecta.f90 b/src/ejecta/module_ejecta.f90 index 669ea8a2..944741c4 100644 --- a/src/ejecta/module_ejecta.f90 +++ b/src/ejecta/module_ejecta.f90 @@ -33,7 +33,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot) type(cratertype),intent(inout) :: crater type(domaintype),intent(in) :: domain integer(I4B),intent(in) :: ejtble - type(ejbtype),dimension(ejtble),intent(in) :: ejb + type(ejbtype),dimension(ejtble),intent(inout) :: ejb real(DP),intent(in) :: deltaMtot end subroutine ejecta_emplace end interface diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index abb30406..d21d1748 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -97,6 +97,7 @@ module module_globals real(DP) :: grad ! Strengthless material transient crater radius real(DP) :: frad ! Final crater radius real(DP) :: fcrat ! Final crater diameter + real(DP) :: ejrad ! Radius to begin ejecta/raised rim real(DP) :: vdepth,vrim ! parameter for parabolic crater form real(DP) :: frim ! parameter for parabolic crater form real(DP) :: rimdis ! crater form radius (bowl + upturned rim)