Changed softening function and added the kappa*t accumulator to preve…

…nt having to run the diffusion solver on every single crater

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 

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

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

src/crater/module_crater.f90

@@ -216,16 +216,27 @@ 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
       type(surftype),dimension(:,:),intent(inout) :: surf
       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