diff --git a/src/crater/crater_realistic_topography.f90 b/src/crater/crater_realistic_topography.f90 index 938c2359..3919ea75 100644 --- a/src/crater/crater_realistic_topography.f90 +++ b/src/crater/crater_realistic_topography.f90 @@ -43,85 +43,88 @@ subroutine crater_realistic_topography(user,surf,crater,domain,ejecta_dem) ! Internal variables real(DP) :: deltaMtot - real(DP), dimension(2) :: rn integer(I4B) :: inc interface - subroutine complex_floor(user,surf,crater,rn,deltaMtot) + subroutine complex_floor(user,surf,crater,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_floor - subroutine complex_peak(user,surf,crater,rn,deltaMtot) + subroutine complex_peak(user,surf,crater,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_peak - subroutine complex_wall_texture(user,surf,crater,domain,rn,deltaMtot) + subroutine complex_wall_texture(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),dimension(:),intent(in) :: rn real(DP),intent(inout) :: deltaMtot end subroutine complex_wall_texture - subroutine complex_terrace(user,surf,crater,rn,deltaMtot) + subroutine complex_terrace(user,surf,crater,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 ejecta_texture(user,surf,crater,rn,deltaMtot,inc,ejecta_dem) + subroutine ejecta_texture(user,surf,crater,deltaMtot,inc,ejecta_dem) 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(in) :: deltaMtot integer(I4B),intent(in) :: inc real(DP),dimension(-inc:inc,-inc:inc),intent(inout) :: ejecta_dem end subroutine ejecta_texture + subroutine crater_realistic_slope_texture(user,critical,inc,critarray) + use module_globals + implicit none + type(usertype),intent(in) :: user + real(DP),intent(in) :: critical + integer(I4B),intent(in) :: inc + real(DP),dimension(-inc:inc,-inc:inc),intent(out) :: critarray + end subroutine crater_realistic_slope_texture + + end interface ! Executable code - call random_number(rn) - - call complex_terrace(user,surf,crater,rn,deltaMtot) - call complex_wall_texture(user,surf,crater,domain,rn,deltaMtot) - call complex_floor(user,surf,crater,rn,deltaMtot) - call complex_peak(user,surf,crater,rn,deltaMtot) - + 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) ! Retrieve the size of the ejecta dem and correct for indexing inc = (size(ejecta_dem,1) - 1) / 2 - - call ejecta_texture(user,surf,crater,rn,deltaMtot,inc,ejecta_dem) + call ejecta_texture(user,surf,crater,deltaMtot,inc,ejecta_dem) + + + ! 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 -subroutine complex_peak(user,surf,crater,rn,deltaMtot) +subroutine complex_peak(user,surf,crater,deltaMtot) ! Makes the central peak or peak ring, depending on the crater size use module_globals use module_util @@ -130,12 +133,12 @@ subroutine complex_peak(user,surf,crater,rn,deltaMtot) 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 + real(DP), dimension(2) :: rn ! Complex crater peak integer(I4B), parameter :: num_octaves = 10 ! Number of Perlin noise octaves @@ -151,6 +154,10 @@ subroutine complex_peak(user,surf,crater,rn,deltaMtot) integer(I4B) :: octave real(DP) :: noise + + ! Executable code + call random_number(rn) + !Preliminary empirical fits based on a handful of complex craters !More work needs to be done to make these models more robust FWHM = 0.1_DP @@ -198,7 +205,7 @@ subroutine complex_peak(user,surf,crater,rn,deltaMtot) end subroutine complex_peak -subroutine complex_floor(user,surf,crater,rn,deltaMtot) +subroutine complex_floor(user,surf,crater,deltaMtot) ! Makes the flat floor with the appearance of blocks and smooth melts use module_globals use module_util @@ -207,12 +214,12 @@ subroutine complex_floor(user,surf,crater,rn,deltaMtot) 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 + real(DP), dimension(2) :: rn ! Complex crater floor parameters integer(I4B), parameter :: num_octaves = 8 ! Number of Perlin noise octaves @@ -222,12 +229,14 @@ subroutine complex_floor(user,surf,crater,rn,deltaMtot) 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 + !Executable code + call random_number(rn) + inc = max(min(nint(0.5_DP * crater%floordiam / user%pix),PBCLIM*user%gridsize),1) + 1 crater%maxinc = max(crater%maxinc,inc) @@ -265,7 +274,7 @@ subroutine complex_floor(user,surf,crater,rn,deltaMtot) end subroutine complex_floor -subroutine complex_wall_texture(user,surf,crater,domain,rn,deltaMtot) +subroutine complex_wall_texture(user,surf,crater,domain,deltaMtot) ! Gives the terraced slopes some rough texture to increase realism use module_globals use module_util @@ -277,12 +286,12 @@ subroutine complex_wall_texture(user,surf,crater,domain,rn,deltaMtot) type(surftype),dimension(:,:),intent(inout) :: surf type(cratertype),intent(inout) :: crater type(domaintype),intent(in) :: domain - real(DP),dimension(:),intent(in) :: rn real(DP),intent(inout) :: deltaMtot ! Internal variables real(DP) :: newdem,elchange,xbar,ybar,areafrac,flr,r integer(I4B) :: i,j,inc,xpi,ypi + real(DP), dimension(2) :: rn ! Topographic noise parameters real(DP) :: xynoise, znoise @@ -298,6 +307,8 @@ subroutine complex_wall_texture(user,surf,crater,domain,rn,deltaMtot) real(DP), parameter :: freq = 2.0_DP ! Spatial size scale factor multiplier at each octave level real(DP), parameter :: pers = 1.20_DP ! The relative size scaling at each octave level + !Executable code + call random_number(rn) inc = max(min(nint(2.1_DP * crater%frad / user%pix),PBCLIM*user%gridsize),1) + 1 crater%maxinc = max(crater%maxinc,inc) @@ -345,7 +356,7 @@ subroutine complex_wall_texture(user,surf,crater,domain,rn,deltaMtot) end subroutine complex_wall_texture -subroutine complex_terrace(user,surf,crater,rn,deltaMtot) +subroutine complex_terrace(user,surf,crater,deltaMtot) !Makes terraced walls by applying noisy topographic diffusion in discrete radial zones along the wall use module_globals use module_util @@ -356,12 +367,12 @@ subroutine complex_terrace(user,surf,crater,rn,deltaMtot) 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,rad,xbar,ybar,r,flr,hprof,tprof integer(I4B) :: xpi,ypi,i,j,inc,maxhits + real(DP), dimension(2) :: rn ! Topographic noise parameters real(DP) :: xynoise, znoise @@ -382,6 +393,8 @@ subroutine complex_terrace(user,surf,crater,rn,deltaMtot) real(DP) :: rinner ! The radius of the terrace outer edge real(DP) :: upshift,dfloor + !Executable code + call random_number(rn) nscallops = 16 rad = 2.0_DP * crater%frad @@ -496,7 +509,8 @@ subroutine complex_terrace(user,surf,crater,rn,deltaMtot) end subroutine complex_terrace -subroutine ejecta_texture(user,surf,crater,rn,deltaMtot,inc,ejecta_dem) +subroutine ejecta_texture(user,surf,crater,deltaMtot,inc,ejecta_dem) + ! Adds realistic texture to the ejecta use module_globals use module_util use module_crater, EXCEPT_THIS_ONE => crater_realistic_topography @@ -506,7 +520,6 @@ subroutine ejecta_texture(user,surf,crater,rn,deltaMtot,inc,ejecta_dem) 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(in) :: deltaMtot integer(I4B),intent(in) :: inc real(DP),dimension(-inc:inc,-inc:inc),intent(inout) :: ejecta_dem @@ -516,6 +529,7 @@ subroutine ejecta_texture(user,surf,crater,rn,deltaMtot,inc,ejecta_dem) integer(I4B) :: i,j,xpi,ypi integer(I4B),dimension(2,-inc:inc,-inc:inc) :: indarray real(DP),dimension(-inc:inc,-inc:inc) :: noise_dem + real(DP), dimension(2) :: rn ! Topographic noise parameters real(DP) :: xynoise, znoise, hprof,xysplat,zsplat,dsplat,splatnoise @@ -543,6 +557,9 @@ subroutine ejecta_texture(user,surf,crater,rn,deltaMtot,inc,ejecta_dem) logical :: insplat real(DP),parameter :: splatmag = 0.15_DP ! The magnitude of the splat features relative to the ejecta thickness + !Executable code + call random_number(rn) + nsplats = 32 ! Get the ejecta mass @@ -646,3 +663,56 @@ subroutine ejecta_texture(user,surf,crater,rn,deltaMtot,inc,ejecta_dem) return end subroutine ejecta_texture + +subroutine crater_realistic_slope_texture(user,critical,inc,critarray) + ! Adds noise to the critical slope to give texture to regions that undergo slope collapse + use module_globals + use module_util + use module_crater, EXCEPT_THIS_ONE => crater_realistic_topography + implicit none + + ! Arguments + type(usertype),intent(in) :: user + real(DP),intent(in) :: critical + integer(I4B),intent(in) :: inc + real(DP),dimension(-inc:inc,-inc:inc),intent(out) :: critarray + + ! Internal variables + real(DP) :: xynoise, znoise,xbar,ybar,r,noise + integer(I4B) :: octave,i,j + real(DP), dimension(2) :: rn + + ! Topographic noise parameters + 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 = 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 :: 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 + + ! Executable code + call random_number(rn) + + do j = -inc,inc + do i = -inc,inc + + xbar = real(i,kind=DP) + ybar = real(j,kind=DP) + + r = sqrt(xbar**2 + ybar**2) + + noise = 0.0_DP + do octave = 1, num_octaves + xynoise = xy_noise_fac * freq ** (octave - 1) + znoise = noise_height * pers ** (octave - 1) + noise = noise + util_perlin_noise(xynoise * xbar + offset * rn(1), & + xynoise * ybar + offset * rn(2)) * znoise + end do + !write(*,*) i,j,noise + critarray(i,j) = max(critical * (1.0_DP + noise),0.0_DP) + end do + end do + + return +end subroutine crater_realistic_slope_texture + diff --git a/src/crater/crater_slope_collapse.f90 b/src/crater/crater_slope_collapse.f90 index baf029c8..71944fd2 100644 --- a/src/crater/crater_slope_collapse.f90 +++ b/src/crater/crater_slope_collapse.f90 @@ -41,19 +41,17 @@ subroutine crater_slope_collapse(user,surf,crater,domain,critical,deltaMtot) integer(I4B) :: iradsq integer(I4B) :: i,j,inc,incsq,mx1,mx2,mx3,my1,my2,my3,xpi,ypi integer(I4B) :: loopcnt,looplim,mult - real(DP),dimension(:,:),allocatable :: elevarray + real(DP),dimension(:,:),allocatable :: elevarray,critarray real(DP),dimension(:,:),allocatable :: cumulative_elchange,kappat integer(I4B),dimension(:,:,:),allocatable :: indarray logical :: failflag character(len=MESSAGESIZE) :: message ! message for the progress bar - ! Executable Code ! Some preliminary setup diffmax = 0.25_DP * user%pix**2 looplim = 100 * crater%fcratpx - !critical = (CRITSLP * user%pix)**2 ! determine area to effect inc = max(min(nint(1.5_DP * crater%frad / user%pix),ceiling(SQRT2*user%gridsize)),1) + 1 @@ -73,6 +71,7 @@ subroutine crater_slope_collapse(user,surf,crater,domain,critical,deltaMtot) endif allocate(elevarray(-inc:inc,-inc:inc)) + allocate(critarray(-inc:inc,-inc:inc)) allocate(cumulative_elchange(-inc:inc,-inc:inc)) allocate(indarray(6,-inc:inc,-inc:inc)) allocate(kappat(-inc:inc,-inc:inc)) @@ -80,14 +79,24 @@ subroutine crater_slope_collapse(user,surf,crater,domain,critical,deltaMtot) cumulative_elchange = 0._DP + + + if (user%dorealistic) then + call crater_realistic_slope_texture(user,critical,inc,critarray) + write(*,*) 'critical, max min, ',critical,minval(critarray),maxval(critarray) + else + critarray(:,:) = critical + end if + ! begin downslope motion loop do loopcnt=1,looplim failflag = .false. - kappat = 0.0_DP + kappat(:,:) = 0.0_DP + ! loop over affected matrix area !$OMP PARALLEL DO DEFAULT(PRIVATE) IF(inc > INCPAR) & - !$OMP SHARED(inc,loopcnt,incsq,failflag,elevarray,indarray,kappat,diffmax,critical) & + !$OMP SHARED(inc,loopcnt,incsq,failflag,elevarray,indarray,kappat,diffmax,critarray) & !$OMP SHARED(crater,user,surf) do j=-inc,inc do i=-inc,inc @@ -134,7 +143,7 @@ subroutine crater_slope_collapse(user,surf,crater,domain,critical,deltaMtot) tslp3sq = (xslp2*xslp2 + yslp1*yslp1) tslp4sq = (xslp2*xslp2 + yslp2*yslp2) slpsq = max(tslp1sq,tslp2sq,tslp3sq,tslp4sq) - if (slpsq > critical) then + if (slpsq > critarray(i,j)) then kappat(i,j) = diffmax failflag = .true. else @@ -147,7 +156,7 @@ subroutine crater_slope_collapse(user,surf,crater,domain,critical,deltaMtot) if (.not.failflag) exit - elevarray = 0._DP + elevarray(:,:) = 0._DP call util_diffusion_solver(user,surf,2 * inc + 1,indarray(1:2,:,:),kappat,elevarray,mult) @@ -193,6 +202,6 @@ subroutine crater_slope_collapse(user,surf,crater,domain,critical,deltaMtot) end do ! end downslope motion loops deltaMtot = deltaMtot + sum(cumulative_elchange) - deallocate(elevarray,cumulative_elchange,indarray,kappat) + deallocate(elevarray,critarray,cumulative_elchange,indarray,kappat) return end subroutine crater_slope_collapse