From 6783678b9f954db80ef226e1ba01442024bfe09f Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 27 Mar 2020 17:14:50 -0400 Subject: [PATCH] Work in progress on a better terrace model using the scalloped rim model --- src/crater/crater_populate.f90 | 2 +- src/crater/crater_profile.f90 | 55 +++++++++++++ src/crater/crater_realistic_topography.f90 | 91 +++++++++++++--------- src/crater/crater_slope_collapse.f90 | 9 ++- src/crater/module_crater.f90 | 3 +- 5 files changed, 119 insertions(+), 41 deletions(-) create mode 100644 src/crater/crater_profile.f90 diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index fca76d77..79908ab5 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -219,7 +219,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt 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) + if (user%docollapse) call crater_slope_collapse(user,surf,crater,domain,(CRITSLP * user%pix)**2,ejbmass) ! Record crater in an available layer as long as it is above the cutoff call crater_record(user,surf,crater) diff --git a/src/crater/crater_profile.f90 b/src/crater/crater_profile.f90 new file mode 100644 index 00000000..f6dce262 --- /dev/null +++ b/src/crater/crater_profile.f90 @@ -0,0 +1,55 @@ +!********************************************************************************************************************************** +! +! Unit Name : crater_profile +! Unit Type : subroutine +! Project : CTEM +! Language : Fortran 2003 +! +! Description : Function that defines the basic profile of a crater, not including the ejecta blanket +! +! +! Input +! Arguments : +! +! Output +! Arguments : +! +! +! Notes : +! +!********************************************************************************************************************************** +function crater_profile(user,crater,r) result(h) + use module_globals + use module_crater, EXCEPT_THIS_ONE => crater_profile + implicit none + + ! Arguments + type(usertype),intent(in) :: user + type(cratertype),intent(in) :: crater + real(DP),intent(in) :: r + real(DP) :: h + + ! Internal variables + real(DP) :: flrad,c0,c1,c2,c3 + + ! Executable code + 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 = (-crater%floordepth - crater%rimheight) / (flrad + flrad**2 / 3._DP - flrad**3 / 6._DP - 7._DP / 6._DP) + c0 = crater%rimheight - (7._DP / 6._DP) * c1 + c2 = c1 / 3._DP + c3 = -c2 / 2._DP + + if (r < flrad) then + h = -crater%floordepth + if (r > crater%frad) + h = crater%rimheight * ((crater%frad / lrad)**RIMDROP) + else + h = c0 + c1 * r + c2 * r**2 + c3 * r**3 + end if + + + return +end subroutine crater_profile + diff --git a/src/crater/crater_realistic_topography.f90 b/src/crater/crater_realistic_topography.f90 index bb71643e..a28a165e 100644 --- a/src/crater/crater_realistic_topography.f90 +++ b/src/crater/crater_realistic_topography.f90 @@ -105,11 +105,13 @@ end subroutine complex_rim 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_terrace(user,surf,crater,rn,deltaMtot) + call crater_slope_collapse(user,surf,crater,domain,(0.35_DP * user%pix)**2,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) + return end subroutine crater_realistic_topography @@ -132,7 +134,7 @@ subroutine complex_peak(user,surf,crater,rn,deltaMtot) 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 :: xy_noise_fac = 16.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 @@ -272,7 +274,7 @@ subroutine complex_wall_texture(user,surf,crater,rn,deltaMtot) 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-4_DP ! Spatial "size" of noise features at the first octave + real(DP), parameter :: noise_height = 1.0e-5_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 @@ -295,9 +297,10 @@ subroutine complex_wall_texture(user,surf,crater,rn,deltaMtot) 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) + areafrac = areafrac / max(((r - flr) / flr)**2,0.1_DP) ! Add some roughness to the walls noise = 0.0_DP @@ -333,7 +336,7 @@ subroutine complex_rim(user,surf,crater,rn,deltaMtot) real(DP),intent(inout) :: deltaMtot ! Internal variables - real(DP) :: newdem,elchange,rad,xbar,ybar + real(DP) :: newdem,elchange,rad,xbar,ybar,r,areafrac,hprof,tprof integer(I4B) :: xpi,ypi,i,j,inc,maxhits ! Topographic noise parameters @@ -341,7 +344,6 @@ subroutine complex_rim(user,surf,crater,rn,deltaMtot) 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 @@ -352,12 +354,12 @@ subroutine complex_rim(user,surf,crater,rn,deltaMtot) ! First make the interior of the crater rad = 1.25_DP * crater%frad - ! Create diffusion noise for terraces - num_octaves = 4 + ! Create diffusion noise for scalloped rim 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) @@ -373,14 +375,27 @@ subroutine complex_rim(user,surf,crater,rn,deltaMtot) xbar = xpi * user%pix - crater%xl ybar = ypi * user%pix - crater%yl + r = sqrt(xbar**2 + ybar**2) / crater%frad + ! 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 + + + if (r >= 1.0) then + hprof = crater%ejrim * r**(-3) + crater%melev + tprof = 0.5_DP * crater%ejrim + crater%melev !* ((2.0_DP - r)**(-4) - 1.0_DP) + crater%melev + else + hprof = crater%ejrim * (2.0_DP - r)**(-2) + crater%melev + tprof = 0.5_DP * crater%ejrim + crater%melev ! (r**(-4) - 1.0_DP) + crater%melev + end if + + if (crater%melev + crater%ejrim - noise < hprof) then + newdem = tprof - noise + end if elchange = newdem - surf(xpi,ypi)%dem deltaMtot = deltaMtot + elchange @@ -439,15 +454,20 @@ subroutine complex_terrace(user,surf,crater,rn,deltaMtot) ! Create diffusion noise for terraces - num_octaves = 4 + num_octaves = 8 xy_noise_fac = 4.00_DP - diff_height = 6e-5_DP * (crater%fcrat)**2 - noise_height = 0.02_DP + diff_height = 1e-5_DP * (crater%fcrat)**2 + noise_height = 0.01_DP freq = 2.0_DP - pers = 0.80_DP + pers = 0.50_DP terracefac = 8 flr = crater%floordiam / crater%fcrat + + !Borrow from scallops + 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) @@ -491,9 +511,11 @@ subroutine complex_terrace(user,surf,crater,rn,deltaMtot) 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), & + znoise = noise_height * crater%fcrat + dnoise = util_perlin_noise(xynoise * xbar + terrace_num * offset * rn(1), & xynoise * ybar + terrace_num * offset * rn(2))* znoise + + noise = sqrt(sqrt(dnoise**2)) newdem = max(newdem + noise,crater%melev - crater%floordepth) end if elchange = newdem - surf(xpi,ypi)%dem @@ -502,14 +524,13 @@ subroutine complex_terrace(user,surf,crater,rn,deltaMtot) ! 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 + offset = 1500 + do octave = 1, num_octaves + xynoise = 4 * xy_noise_fac * freq ** (octave - 1) / crater%fcrat + znoise = diff_height * pers ** (octave - 1) / max(r,0.3_DP) + noise = noise + util_perlin_noise(xynoise * xbar + offset * rn(1), & + xynoise * ybar + offset * rn(2)) * znoise + end do kdiff(i,j) = noise end do @@ -520,16 +541,16 @@ subroutine complex_terrace(user,surf,crater,rn,deltaMtot) kdiff(:,:) = kdiff(:,:) * areafrac(:,:) !* 0.0_DP maxhits = 1 - call util_diffusion_solver(user,surf,2 * inc + 1,indarray,kdiff,cumulative_elchange,maxhits) - - 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 + !call util_diffusion_solver(user,surf,2 * inc + 1,indarray,kdiff,cumulative_elchange,maxhits) + + !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) diff --git a/src/crater/crater_slope_collapse.f90 b/src/crater/crater_slope_collapse.f90 index 91aba2ae..e9d6860a 100644 --- a/src/crater/crater_slope_collapse.f90 +++ b/src/crater/crater_slope_collapse.f90 @@ -16,7 +16,7 @@ ! Notes : ! !********************************************************************************************************************************** -subroutine crater_slope_collapse(user,surf,crater,domain,deltaMtot) +subroutine crater_slope_collapse(user,surf,crater,domain,critical,deltaMtot) use module_globals use module_util use module_io @@ -28,12 +28,13 @@ subroutine crater_slope_collapse(user,surf,crater,domain,deltaMtot) type(surftype),dimension(:,:),intent(inout) :: surf type(cratertype),intent(inout) :: crater type(domaintype),intent(in) :: domain + real(DP),intent(in) :: critical real(DP),intent(inout) :: deltaMtot ! Internal variables real(DP) :: diffmax,difflim real(DP) :: slpmax,slpsq - real(DP) :: elchgmax,critical + real(DP) :: elchgmax real(DP) :: xslp1,xslp2,yslp1,yslp2 real(DP) :: tslp1sq,tslp2sq,tslp3sq,tslp4sq real(DP) :: xgrad,ygrad,tgrad,gradmax @@ -55,10 +56,10 @@ subroutine crater_slope_collapse(user,surf,crater,domain,deltaMtot) ! Some preliminary setup diffmax = 0.25_DP * user%pix**2 looplim = 100 * crater%fcratpx - critical = (CRITSLP * user%pix)**2 + !critical = (CRITSLP * user%pix)**2 ! determine area to effect - inc = max(min(crater%fcratpx,ceiling(SQRT2*user%gridsize)),1) + 1 + inc = max(min(nint(1.5_DP * crater%frad / user%pix),ceiling(SQRT2*user%gridsize)),1) + 1 crater%maxinc = max(crater%maxinc,inc) incsq = inc**2 diff --git a/src/crater/module_crater.f90 b/src/crater/module_crater.f90 index a573bf97..a3b37775 100644 --- a/src/crater/module_crater.f90 +++ b/src/crater/module_crater.f90 @@ -208,13 +208,14 @@ end subroutine crater_tally_observed end interface interface - subroutine crater_slope_collapse(user,surf,crater,domain,deltaMtot) + subroutine crater_slope_collapse(user,surf,crater,domain,critical,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(in) :: critical real(DP),intent(inout) :: deltaMtot end subroutine crater_slope_collapse end interface