From d382a027ffa5a7e43329204eb80d6b712d392cf9 Mon Sep 17 00:00:00 2001
From: daminton <daminton@purdue.edu>
Date: Mon, 27 Feb 2017 22:27:17 +0000
Subject: [PATCH] Fixed a bug in the superdomain area calculation

---
 src/crater/crater_subpixel_diffusion.f90 | 38 +++++++++++++++++++++---
 1 file changed, 34 insertions(+), 4 deletions(-)

diff --git a/src/crater/crater_subpixel_diffusion.f90 b/src/crater/crater_subpixel_diffusion.f90
index d18fc676..53bfabeb 100644
--- a/src/crater/crater_subpixel_diffusion.f90
+++ b/src/crater/crater_subpixel_diffusion.f90
@@ -38,9 +38,10 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiff
    integer(I4B),dimension(2,0:user%gridsize + 1,0:user%gridsize + 1) :: indarray
    integer(I4B) :: i,j,k,xpi,ypi,n,ntot
    integer(I4B) :: maxhits = 1
-   real(DP) :: lamleft,dburial,lambda,kappat,diam,rad,Abar
+   real(DP) :: lamleft,dburial,lambda,kappat,diam,radius,Area,avgejc
    real(DP),dimension(domain%pnum) :: dN,lambda_regolith,kappat_regolith,lambda_bedrock,kappat_bedrock
 
+
    ! Create box for soften calculation (will be no bigger than the grid itself)
    do j = 0,user%gridsize + 1
       do i = 0,user%gridsize + 1
@@ -55,13 +56,15 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiff
 
    ntot = 1
    ! calculate the subpixel diffusion probability function
+   Area = user%pix**2
+
    do i = 1,domain%pnum
       if ((nflux(1,i) > domain%smallest_crater).and.(nflux(2,i) > domain%smallest_crater)) exit
       ntot = i
       dN(i) = nflux(3,i) * user%interval * finterval
 
-      lambda_bedrock(i) = dN(i) * user%pix**2 
-      lambda_regolith(i) = dN(i) * user%pix**2 
+      lambda_bedrock(i) = dN(i) * Area
+      lambda_regolith(i) = dN(i) * Area
        
       if (user%dosoftening) then
          kappat_bedrock(i) = user%soften_factor * (0.5_DP * nflux(1,i))**(user%soften_slope)
@@ -85,11 +88,38 @@ subroutine crater_subpixel_diffusion(user,surf,prod,nflux,domain,finterval,kdiff
             end if 
             if (diam > domain%smallest_crater) exit
             k = util_poisson(lambda)
-            kdiff(i,j) = kdiff(i,j) + k * kappat / user%pix**2
+            kdiff(i,j) = kdiff(i,j) + k * kappat / Area
          end do
       end do
    end do
+
+   ! Now add the superdomain contribution to crater degradation 
+   if (user%dosoftening) then
+      avgejc = sum(surf%ejcov) / user%gridsize**2
+      do i = 1,domain%pnum
+         if ((user%soften_size * nflux(1,i) <= user%gridsize / SQRT2 ).and.&
+             (user%soften_size * nflux(2,i) <= user%gridsize / SQRT2 )) cycle
+
+         dburial = EXFAC * 0.5_DP * nflux(1,n)
+         if (avgejc > dburial) then
+            radius = nflux(2,i) * 0.5_DP
+         else
+            radius = nflux(1,i) * 0.5_DP
+         end if
+
+         dN(i) = nflux(3,i) * user%interval * finterval
+         Area = min(2 * (user%soften_size * radius)**2 , PI * user%trad**2)
+         Area = max(Area - user%gridsize**2,0.0_DP)
+
+         lambda = dN(i) * Area
+         if (lambda == 0.0_DP) cycle
+         k = util_poisson(lambda)
+          
+         if (k > 0) kdiff = kdiff + k * user%soften_factor / (PI * user%soften_size**2) * radius**(user%soften_slope - 2.0_DP)
+      end do
+   end if
    
+
    call util_diffusion_solver(user,surf,user%gridsize + 2,indarray,kdiff,cumulative_elchange,maxhits)
    do j = 1,user%gridsize
       do i = 1,user%gridsize