Skip to content

Commit

Permalink
Fixed slope collapse routine to conserve mass
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Feb 16, 2017
1 parent cce37c4 commit f78ab12
Showing 1 changed file with 44 additions and 60 deletions.
104 changes: 44 additions & 60 deletions src/crater/crater_slope_collapse.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,33 +33,33 @@ subroutine crater_slope_collapse(user,surf,crater,domain,deltaMtot)
! Internal variables
real(DP) :: diffmax,difflim
real(DP) :: slpmax,slpsq,oldslpmax
real(DP) :: elchange,elchgmax,critical
real(DP) :: elchgmax,critical
real(DP) :: xslp1,xslp2,yslp1,yslp2
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
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 :: gradarray,slparray
real(DP),dimension(:,:),allocatable :: cumulative_elchange
real(DP),dimension(:,:),allocatable :: cumulative_elchange,kappat
integer(I4B),dimension(:,:,:),allocatable :: indarray
logical :: failflag,oldfail
logical :: failflag
character(len=MESSAGESIZE) :: message ! message for the progress bar
real(DP) :: rfac

! testing
!character(STRMAX) :: filename
!integer(I4B) :: ioerr
character(STRMAX) :: filename
integer(I4B) :: ioerr

! Executable Code

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

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

Expand All @@ -76,38 +76,30 @@ subroutine crater_slope_collapse(user,surf,crater,domain,deltaMtot)
endif

allocate(elevarray(-inc:inc,-inc:inc))
allocate(gradarray(-inc:inc,-inc:inc))
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)

kappat = 0.0_DP
do loopcnt=1,looplim
! reset slope failure flag
oldfail = failflag
failflag = .false.

! loop over affected matrix area
!$OMP PARALLEL DO DEFAULT(PRIVATE) IF(inc > INCPAR) &
!$OMP SHARED(inc,loopcnt,incsq,critical,difflim,failflag,elevarray,gradarray,slparray,cumulative_elchange,indarray) &
!$OMP SHARED(inc,loopcnt,incsq,failflag,elevarray,indarray,kappat,diffmax,critical) &
!$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
if (loopcnt==1) then
mx1 = crater%xlpx+i
my1 = crater%ylpx+j
mx1 = crater%xlpx + i
my1 = crater%ylpx + j
mx2 = mx1 + 1
my2 = my1 + 1
mx3 = mx1 - 1
Expand All @@ -117,12 +109,12 @@ subroutine crater_slope_collapse(user,surf,crater,domain,deltaMtot)
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)
Expand All @@ -132,7 +124,7 @@ subroutine crater_slope_collapse(user,surf,crater,domain,deltaMtot)
mx3 = indarray(5,i,j)
my3 = indarray(6,i,j)
end if
iradsq = i*i+j*j
iradsq = i**2 + j**2
if (iradsq <= incsq) then ! round corners
! define location & distance from crater center in pixel space

Expand All @@ -148,41 +140,35 @@ subroutine crater_slope_collapse(user,surf,crater,domain,deltaMtot)
tslp3sq = (xslp2*xslp2 + yslp1*yslp1)
tslp4sq = (xslp2*xslp2 + yslp2*yslp2)
slpsq = max(tslp1sq,tslp2sq,tslp3sq,tslp4sq)
critical = (crater_critical_slope(user,crater,iradsq)*user%gridsize)**2
if (slpsq > critical) then
kappat(i,j) = diffmax
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
else
kappat(i,j) = 0.0_DP
end if
end if
end do
end do ! end area loopover
!$OMP END PARALLEL DO

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

elevarray = 0._DP

call util_diffusion_solver(user,surf,2 * inc + 1,indarray(1:2,:,:),kappat,elevarray,mult)

! 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)
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 + elevarray(i,j)
surf(xpi,ypi)%ejcov = max(surf(xpi,ypi)%ejcov + elevarray(i,j),0.0_DP)
end do
end do

cumulative_elchange = cumulative_elchange + elevarray

!***************************************************
! TESTING -- makes frames for a slope collapse movie
!***************************************************
Expand All @@ -203,18 +189,16 @@ subroutine crater_slope_collapse(user,surf,crater,domain,deltaMtot)
! end if
!end if
!***************************************************


! if stuck, get out
elchgmax = maxval(elevarray)
if (failflag.and.(abs(elchgmax) < domain%small)) exit

if (slpmax >= oldslpmax) exit
oldslpmax = slpmax


if (.not.failflag) exit
end do ! end downslope motion loops
deltaMtot = deltaMtot + sum(cumulative_elchange)
deallocate(elevarray,gradarray,slparray,cumulative_elchange,indarray)

deltaMtot = deltaMtot + sum(cumulative_elchange)
deallocate(elevarray,cumulative_elchange,indarray)
return
end subroutine crater_slope_collapse

0 comments on commit f78ab12

Please sign in to comment.