From f78ab1273d379d94edf4defb5dd3bdf5afeefd8e Mon Sep 17 00:00:00 2001 From: daminton Date: Thu, 16 Feb 2017 19:14:46 +0000 Subject: [PATCH] Fixed slope collapse routine to conserve mass --- src/crater/crater_slope_collapse.f90 | 104 ++++++++++++--------------- 1 file changed, 44 insertions(+), 60 deletions(-) diff --git a/src/crater/crater_slope_collapse.f90 b/src/crater/crater_slope_collapse.f90 index 7989bb6b..7fda521b 100644 --- a/src/crater/crater_slope_collapse.f90 +++ b/src/crater/crater_slope_collapse.f90 @@ -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 @@ -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 @@ -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) @@ -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 @@ -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 !*************************************************** @@ -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