From b0c96e29d34f91e5b1735ae02bb62442614f7371 Mon Sep 17 00:00:00 2001
From: daminton <daminton@purdue.edu>
Date: Tue, 7 Mar 2017 02:29:15 +0000
Subject: [PATCH] Updated test code for crater_soften

---
 src/crater/crater_soften.f90 | 20 +++++++++++++-------
 1 file changed, 13 insertions(+), 7 deletions(-)

diff --git a/src/crater/crater_soften.f90 b/src/crater/crater_soften.f90
index 36a84a1c..d6943b76 100644
--- a/src/crater/crater_soften.f90
+++ b/src/crater/crater_soften.f90
@@ -33,16 +33,22 @@ subroutine crater_soften(user,surf,crater,domain)
    ! Internal variables
    real(DP),dimension(:,:),allocatable :: cumulative_elchange,kappat
    integer(I4B),dimension(:,:,:),allocatable :: indarray
-   integer(I4B),parameter :: MAXHITS = 1
+   integer(I4B) :: maxhits 
 
-   integer(I4B) :: inc,incsq,N,xpi,ypi,iradsq,i,j
-   real(DP) :: xp,yp,fradsq,xbar,ybar,areafrac,kappatmax,ss,sf
+   integer(I4B) :: inc,incsq,N,xpi,ypi,iradsq,i,j,ss,ioerr
+   real(DP) :: xp,yp,fradsq,xbar,ybar,areafrac,kappatmax,sf
 
-   kappatmax = user%soften_factor * crater%frad**user%soften_slope
-   
-   inc = max(min(int(crater%fradpx + 2),PBCLIM * user%gridsize),2) 
+   !kappatmax = user%soften_factor * crater%frad**user%soften_slope
+   open(unit=22,file='soften.dat',status='old',iostat=ioerr)
+   if (ioerr /= 0) return
+   read(22,*) sf
+   close(22)
+   !inc = max(min(int(crater%fradpx + 2),PBCLIM * user%gridsize),2) 
+   inc = crater%fradpx + 2
+   maxhits = (1 + (inc / (user%gridsize / 2)))**2
    crater%maxinc = max(crater%maxinc,inc)
    incsq = inc**2
+   kappatmax = (sf / PI) * crater%frad**2 
 
    allocate(kappat(-inc:inc,-inc:inc))
    allocate(indarray(2,-inc:inc,-inc:inc))
@@ -75,7 +81,7 @@ subroutine crater_soften(user,surf,crater,domain)
             indarray(2,i,j) = ypi
 
             ! We set the diffusion constant to be proportional to the fraction of pixel area covered by the interior of the crater
-            areafrac = util_area_intersection(crater%frad,xbar,ybar,user%pix)
+            areafrac = util_area_intersection(crater%frad*0.98_DP,xbar,ybar,user%pix)
              
             kappat(i,j) = kappatmax * areafrac  ! This is the extra per-crater diffusion required to match equilibrium
          end if