diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index a04cc516..05651f65 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -71,6 +71,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt character(len=MESSAGESIZE) :: message ! message for the progress bar real(DP) :: ejbmass logical :: makecrater + real(DP),dimension(user%gridsize,user%gridsize) :: kdiff TARGET :: surf ! ejecta blanket array @@ -172,7 +173,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt if (user%doseismic) call seismic_shake(user,surf,crater,domain) ! Generate dynamic diffusion - !if (user%dosoftening) call crater_soften(user,surf,crater,domain) + if (user%dosoftening) call crater_soften_accumulate(user,surf,crater,domain,kdiff) ! find the average height and slope at crater location call crater_averages(user,surf,crater) @@ -253,7 +254,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt call io_updatePbar(message) craters_since_subpixel = icrater - icrater_last_subpixel finterval = craters_since_subpixel / real(ntotcrat,kind=DP) - call crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval) + call crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiff) icrater_last_subpixel = icrater end if ! Intermediate tally step diff --git a/src/crater/crater_soften_accumulate.f90 b/src/crater/crater_soften_accumulate.f90 index cc9fcf12..4b116249 100644 --- a/src/crater/crater_soften_accumulate.f90 +++ b/src/crater/crater_soften_accumulate.f90 @@ -34,9 +34,10 @@ subroutine crater_soften_accumulate(user,surf,crater,domain,kdiff) ! Internal variables real(DP) :: SOFTEN_FACTOR = 1.00e-3_DP ! Constant in topographic diffusion term for crater softening + integer(I4B) :: SOFTEN_SIZE = 100 ! Constant in topographic diffusion term for crater softening integer(I4B) :: inc,incsq,N,xpi,ypi,iradsq,i,j - real(DP) :: kappatmax,lrad,lradsq,xp,yp,fradsq + real(DP) :: kappatmax,lrad,lradsq,xp,yp,fradsq,areafrac,xbar,ybar ! TESTING !open(unit=55,file="SOFTEN_FACTOR.test",status="old") @@ -45,10 +46,8 @@ subroutine crater_soften_accumulate(user,surf,crater,domain,kdiff) !******** - kappatmax = SOFTEN_FACTOR * crater%fcrat**2 - - inc = int(crater%frad/user%pix*(domain%small/crater%rheight)**(-1._DP/RIMDROP)) ! Maximum distance of crater form - inc = max(min(max(crater%rimdispx,inc),PBCLIM*user%gridsize),1) + kappatmax = SOFTEN_FACTOR * crater%frad**2 + inc = SOFTEN_SIZE * crater%fradpx + 2 crater%maxinc = max(crater%maxinc,inc) fradsq = crater%frad**2 incsq = inc**2 @@ -66,19 +65,20 @@ subroutine crater_soften_accumulate(user,surf,crater,domain,kdiff) ! Find distance from crater center to current pixel center in real space xp = xpi * user%pix yp = ypi * user%pix + + xbar = xp - crater%xl + ybar = yp - crater%yl lradsq = (crater%xl - xp)**2 + (crater%yl - yp)**2 ! periodic boundary conditions call util_periodic(xpi,ypi,user%gridsize) + areafrac = util_area_intersection(SOFTEN_SIZE * crater%frad,xbar,ybar,user%pix) - ! interior of the crater should have a constant kappa*t, while the - ! rim should fall away with the power law drop as the rim profile if (lradsq < fradsq) then - kdiff(xpi,ypi) = kdiff(xpi,ypi) + kappatmax + kdiff(xpi,ypi) = 0.0_DP else - lrad = sqrt(lradsq) - kdiff(xpi,ypi) = kdiff(xpi,ypi) + kappatmax * ((crater%frad / lrad)**RIMDROP) + kdiff(xpi,ypi) = kappatmax * areafrac end if end if diff --git a/src/crater/crater_subpixel_diffusion.f90 b/src/crater/crater_subpixel_diffusion.f90 index c8633d22..eca2f03d 100644 --- a/src/crater/crater_subpixel_diffusion.f90 +++ b/src/crater/crater_subpixel_diffusion.f90 @@ -19,7 +19,7 @@ ! Notes : ! !********************************************************************************************************************************** -subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval) +subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiffin) use module_globals use module_util use module_crater, EXCEPT_THIS_ONE => crater_subpixel_diffusion @@ -31,6 +31,7 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval) real(DP),dimension(:,:),intent(in) :: prod,nflux type(domaintype),intent(in) :: domain real(DP),intent(in) :: finterval + real(DP),dimension(:,:),intent(inout) :: kdiffin ! Internal variables real(DP),dimension(0:user%gridsize + 1,0:user%gridsize + 1) :: cumulative_elchange,kdiff @@ -48,6 +49,7 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval) call util_periodic(xpi,ypi,user%gridsize) indarray(1,i,j) = xpi indarray(2,i,j) = ypi + kdiff(i,j) = kdiffin(xpi,ypi) end do end do !goto 100 @@ -68,7 +70,6 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval) end do - kdiff = 0._DP do j = 1,user%gridsize do i = 1,user%gridsize do n = 1,ntot @@ -86,7 +87,7 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval) k = util_poisson(lambda) !rad = 0.5_DP * diam !/ user%pix !Abar = PI * rad**2 / user%pix**2 !2 * sqrt(0.5_DP * PI) * rad**2 / (rad + 1.0_DP / SQRT2)**2 - 0.024_DP * rad**0.682_DP - kdiff(i,j) = kdiff(i,j) + k * kappat + kdiff(i,j) = kdiff(i,j) + k * kappat / user%pix**2 end do end do end do @@ -110,6 +111,7 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval) surf(i,j)%ejcov = max(surf(i,j)%ejcov + cumulative_elchange(i,j),0.0_DP) end do end do + kdiffin = 0.0_DP return end subroutine crater_subpixel_diffusion diff --git a/src/crater/module_crater.f90 b/src/crater/module_crater.f90 index b6af42c8..00b7d167 100644 --- a/src/crater/module_crater.f90 +++ b/src/crater/module_crater.f90 @@ -216,9 +216,19 @@ subroutine crater_soften(user,surf,crater,domain) end subroutine crater_soften end interface - interface - subroutine crater_subpixel_diffusion(user,surf,prod,crtscl,domain,finterval) + subroutine crater_soften_accumulate(user,surf,crater,domain,kdiff) + 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(inout) :: kdiff + end subroutine crater_soften_accumulate + end interface + interface + subroutine crater_subpixel_diffusion(user,surf,prod,crtscl,domain,finterval,kdiffin) use module_globals implicit none type(usertype),intent(in) :: user @@ -226,6 +236,7 @@ subroutine crater_subpixel_diffusion(user,surf,prod,crtscl,domain,finterval) real(DP),dimension(:,:),intent(in) :: prod,crtscl type(domaintype),intent(in) :: domain real(DP),intent(in) :: finterval + real(DP),dimension(:,:),intent(inout) :: kdiffin end subroutine crater_subpixel_diffusion end interface