Skip to content

Commit

Permalink
Reverted back to old style of crater slope collapse
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Feb 3, 2017
1 parent 24f33cb commit 2bce361
Showing 1 changed file with 56 additions and 78 deletions.
134 changes: 56 additions & 78 deletions src/crater/crater_slope_collapse.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,15 +37,14 @@ subroutine crater_slope_collapse(user,surf,crater,domain)
real(DP) :: tslp1sq,tslp2sq,tslp3sq,tslp4sq
real(DP) :: xgrad,ygrad,tgrad,gradmax
integer(I4B) :: iradsq
integer(I4B) :: i,j,inc,incsq,mx1,mx2,mx3,my1,my2,my3,xpi,ypi
integer(I4B) :: i,j,inc,incsq,mx1,mx2,mx3,my1,my2,my3
integer(I4B) :: loopcnt,looplim,mult
real(DP),dimension(:,:),allocatable :: elevarray
real(DP),dimension(:,:),allocatable :: gradarray,slparray
real(DP),dimension(:,:),allocatable :: cumulative_elchange,kappat
real(DP),dimension(:,:),allocatable :: cumulative_elchange
integer(I4B),dimension(:,:,:),allocatable :: indarray
logical :: failflag
logical :: failflag,oldfail
character(len=MESSAGESIZE) :: message ! message for the progress bar
real(DP) :: rfac

! testing
!character(STRMAX) :: filename
Expand All @@ -54,11 +53,12 @@ subroutine crater_slope_collapse(user,surf,crater,domain)
! Executable Code

! Some preliminary setup
diffmax = 0.25_DP * user%pix**2
looplim = 1000 * crater%fcratpx
diffmax = 0.25_DP
looplim = 10*crater%fcratpx
critical = (CRITSLP*user%pix)**2

! determine area to effect
inc = max(min(crater%fcratpx,ceiling(SQRT2*user%gridsize)),1) + 1
inc = max(min(crater%fcratpx,ceiling(SQRT2*user%gridsize)),1)
crater%maxinc = max(crater%maxinc,inc)
incsq = inc*inc

Expand All @@ -79,33 +79,31 @@ subroutine crater_slope_collapse(user,surf,crater,domain)
allocate(slparray(-inc:inc,-inc:inc))
allocate(cumulative_elchange(-inc:inc,-inc:inc))
allocate(indarray(6,-inc:inc,-inc:inc))
allocate(kappat(-inc:inc,-inc:inc))

!mult = inc/(user%gridsize/2) + 1 ! Reduce max diffusion if the loopover area exceeds the periodic boundary conditions
mult = (1 + (inc/(user%gridsize/2)))**2
!difflim = diffmax/mult
difflim = diffmax/mult

cumulative_elchange = 0._DP

! begin downslope motion loop
failflag=.true.
oldslpmax = huge(oldslpmax)
slparray = 0.0_DP
kappat = 0.0_DP

do loopcnt=1,looplim
! reset slope failure flag
!if (mod(loopcnt,100)==0) write(*,*) loopcnt,looplim
oldfail = failflag
failflag = .false.

! loop over affected matrix area
!$OMP PARALLEL DO DEFAULT(PRIVATE) IF(inc > INCPAR) &
!$OMP SHARED(inc,loopcnt,incsq,failflag,elevarray,gradarray,slparray,cumulative_elchange,indarray,kappat) &
!$OMP SHARED(inc,loopcnt,incsq,critical,difflim,failflag,elevarray,gradarray,slparray,cumulative_elchange,indarray) &
!$OMP SHARED(crater,user,surf)
do j=-inc,inc
do i=-inc,inc
elevarray(i,j) = 0._DP
gradarray(i,j) = 0._DP
slparray(i,j) = 0._DP
elevarray(i,j)=0._DP
gradarray(i,j)=0._DP
slparray(i,j)=0._DP
if (loopcnt==1) then
mx1 = crater%xlpx+i
my1 = crater%ylpx+j
Expand All @@ -118,20 +116,20 @@ subroutine crater_slope_collapse(user,surf,crater,domain)
call util_periodic(mx1,my1,user%gridsize)
call util_periodic(mx2,my2,user%gridsize)
call util_periodic(mx3,my3,user%gridsize)
indarray(1,i,j) = mx1
indarray(2,i,j) = my1
indarray(3,i,j) = mx2
indarray(4,i,j) = my2
indarray(5,i,j) = mx3
indarray(6,i,j) = my3
indarray(1,i,j)=mx1
indarray(2,i,j)=my1
indarray(3,i,j)=mx2
indarray(4,i,j)=my2
indarray(5,i,j)=mx3
indarray(6,i,j)=my3

else
mx1 = indarray(1,i,j)
my1 = indarray(2,i,j)
mx2 = indarray(3,i,j)
my2 = indarray(4,i,j)
mx3 = indarray(5,i,j)
my3 = indarray(6,i,j)
mx1=indarray(1,i,j)
my1=indarray(2,i,j)
mx2=indarray(3,i,j)
my2=indarray(4,i,j)
mx3=indarray(5,i,j)
my3=indarray(6,i,j)
end if
iradsq = i*i+j*j
if (iradsq <= incsq) then ! round corners
Expand All @@ -148,54 +146,41 @@ subroutine crater_slope_collapse(user,surf,crater,domain)
tslp2sq = (xslp1*xslp1 + yslp2*yslp2)
tslp3sq = (xslp2*xslp2 + yslp1*yslp1)
tslp4sq = (xslp2*xslp2 + yslp2*yslp2)
slpsq = max(tslp1sq,tslp2sq,tslp3sq,tslp4sq)
slparray(i,j) = slpsq
critical = crater_critical_slope(user,crater,iradsq) * user%pix**2
kappat(i,j) = max(slpsq - critical,0.0_DP)
if (slpsq > critical) failflag = .true.
! ! find elevation change
! xgrad = (xslp1 - xslp2)
! ygrad = (yslp1 - yslp2)
! tgrad = xgrad + ygrad
! elchange = difflim * tgrad
!
! ! change digital elevation map
! elevarray(i,j) = elchange
! gradarray(i,j) = tgrad
! cumulative_elchange(i,j) = cumulative_elchange(i,j) + elchange
! end if
slpsq=max(tslp1sq,tslp2sq,tslp3sq,tslp4sq)

if (slpsq > critical) then
failflag = .true.
! find elevation change
xgrad = (xslp1 - xslp2)
ygrad = (yslp1 - yslp2)
tgrad = xgrad + ygrad
elchange = difflim * tgrad

! change digital elevation map
elevarray(i,j)=elchange
gradarray(i,j)=tgrad
slparray(i,j)=slpsq
cumulative_elchange(i,j)=cumulative_elchange(i,J) + elchange
end if
end if
end do
end do ! end area loopover
!$OMP END PARALLEL DO

if (.not.failflag) exit
kappat = kappat / maxval(kappat) * diffmax * 2_DP

call util_diffusion_solver(user,surf,2 * inc + 1,indarray(1:2,:,:),kappat,cumulative_elchange,mult)
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

! See if we have reached maximums
!elchgmax = maxval(elevarray)
slpmax = sqrt(maxval(slparray)) / user%pix
gradmax = maxval(gradarray)
elchgmax=maxval(elevarray)
slpmax=sqrt(maxval(slparray))/user%pix
gradmax=maxval(gradarray)

! Add the elevation changes back to the dem
!do j=-inc,inc
! do i=-inc,inc
! mx1=indarray(1,i,j)
! my1=indarray(2,i,j)
! surf(mx1,my1)%dem = surf(mx1,my1)%dem + elevarray(i,j)
! surf(mx1,my1)%ejcov = max(surf(mx1,my1)%ejcov + elevarray(i,j),0.0_DP)
! end do
!end do
do j=-inc,inc
do i=-inc,inc
mx1=indarray(1,i,j)
my1=indarray(2,i,j)
surf(mx1,my1)%dem = surf(mx1,my1)%dem + elevarray(i,j)
surf(mx1,my1)%ejcov = max(surf(mx1,my1)%ejcov + elevarray(i,j),0.0_DP)
end do
end do

!***************************************************
! TESTING -- makes frames for a slope collapse movie
Expand All @@ -219,21 +204,14 @@ subroutine crater_slope_collapse(user,surf,crater,domain)
!***************************************************

! if stuck, get out
!if (failflag.and.(abs(elchgmax) < 1.0e-8*domain%small)) exit
if (failflag.and.(abs(elchgmax) < domain%small)) exit

!if (slpmax >= oldslpmax) exit
!oldslpmax = slpmax
if (slpmax >= oldslpmax) exit
oldslpmax = slpmax


if (.not.failflag) exit
end do ! end downslope motion loops
!j = 0
!do i = -inc,inc
! iradsq = i**2 + j**2
! critical = crater_critical_slope(user,crater,iradsq)
! write(25,*) i*user%pix/crater%frad,slparray(i,j),critical*user%pix**2
!end do


deallocate(elevarray,gradarray,slparray,cumulative_elchange,indarray)

Expand Down

0 comments on commit 2bce361

Please sign in to comment.