Checked out from main branch

src/crater/crater_soften.f90

@@ -31,74 +31,115 @@ subroutine crater_soften(user,surf,crater,domain)
    type(domaintype),intent(in) :: domain
 
    ! Internal variables
-   real(DP),dimension(:,:),allocatable :: cumulative_elchange,kappat
+   real(DP),dimension(:,:),allocatable :: cumulative_elchange,kdiff
    integer(I4B),dimension(:,:,:),allocatable :: indarray
+   real(DP),dimension(0:user%gridsize+1,0:user%gridsize+1) :: big_cel,big_kdiff
+   integer(I4B),dimension(2,0:user%gridsize+1,0:user%gridsize+1) :: big_indarray
    integer(I4B) :: maxhits 
 
-   integer(I4B) :: inc,incsq,N,xpi,ypi,iradsq,i,j,ss,ioerr
-   real(DP) :: xp,yp,fradsq,xbar,ybar,areafrac,kappatmax,sf
+   integer(I4B) :: inc,incsq,N,xpi,ypi,iradsq,i,j,ss,ioerr,bigi,bigj
+   real(DP) :: xp,yp,fradsq,xbar,ybar,areafrac,kdiffmax,sf
 
-   !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
+   kdiffmax = user%Kd1 * crater%frad**user%psi
+   inc = int(min(user%fe * crater%frad / user%pix,SQRT2 * user%gridsize)) + 3 
    maxhits = (1 + (inc / (user%gridsize / 2)))**2
    crater%maxinc = max(crater%maxinc,inc)
-   incsq = inc**2
-   kappatmax = (sf / PI) * crater%frad**2 
+   incsq = (inc - 2)**2
 
-   allocate(kappat(-inc:inc,-inc:inc))
+   allocate(kdiff(-inc:inc,-inc:inc))
    allocate(indarray(2,-inc:inc,-inc:inc))
    allocate(cumulative_elchange(-inc:inc,-inc:inc))
-   cumulative_elchange = 0._DP
+   cumulative_elchange = 0.0_DP
    indarray = inc
-   kappat = 0.0_DP
+   kdiff = 0.0_DP
 
    ! Loop over affected matrix area
    do j=-inc,inc  ! Do the loop in pixel space
       do i=-inc,inc
+         ! find elevation and grid point
+         xpi = crater%xlpx + i
+         ypi = crater%ylpx + j
+
+         ! periodic boundary conditions
+         call util_periodic(xpi,ypi,user%gridsize)
+
+         indarray(1,i,j) = xpi
+         indarray(2,i,j) = ypi
+
          ! find distance from crater center
          iradsq = i*i + j*j
-         if (iradsq <= incsq) then
-            ! find elevation and grid point
-            xpi = crater%xlpx + i
-            ypi = crater%ylpx + j
 
+         ! We set the diffusion constant to be proportional to the fraction of pixel area covered by the interior of the crater
+         !if (iradsq <= incsq) then
             ! 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
 
-            ! periodic boundary conditions
-            call util_periodic(xpi,ypi,user%gridsize)
-
-            indarray(1,i,j) = xpi
-            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*0.98_DP,xbar,ybar,user%pix)
+            areafrac = util_area_intersection(user%fe * crater%frad,xbar,ybar,user%pix)
 
-            kappat(i,j) = kappatmax * areafrac  ! This is the extra per-crater diffusion required to match equilibrium
-         end if
+            kdiff(i,j) = kdiffmax * areafrac  ! Apply user-defined diffusive degradation to degradation region
+         !end if
       end do
    end do !end area loopover 
+   !write(*,*) 'crater_soften: ',crater%frad,maxval(kdiff)
 
+   if (2 * inc + 1 < user%gridsize) then
+      write(*,*) 'crater inc = ',inc
+      write(*,*) 'crater kdiff '
+      write(*,*) maxval(kdiff)
+      call util_diffusion_solver(user,surf,2 * inc + 1,indarray,kdiff,cumulative_elchange,maxhits)
+      write(16,*) crater%frad
+      do j = -inc,inc
+         do i = -inc,inc
+            xpi = indarray(1,i,j)
+            ypi = indarray(2,i,j)
+            write(16,*) i,j,xpi,ypi,kdiff(i,j),surf(xpi,ypi)%dem,cumulative_elchange(i,j)
+            !surf(xpi,ypi)%dem = surf(xpi,ypi)%dem + cumulative_elchange(i,j)
+            !surf(xpi,ypi)%ejcov = max(surf(xpi,ypi)%ejcov + cumulative_elchange(i,j),0.0_DP)
+         end do
+      end do
+   else
+      maxhits = 1
+      big_kdiff = 0.0_DP
+      big_cel = 0.0_DP
+      do bigj = 0,user%gridsize + 1
+         do bigi = 0,user%gridsize + 1
+            xpi = bigi
+            ypi = bigj
+            call util_periodic(xpi,ypi,user%gridsize)
+            big_indarray(1,bigi,bigj) = xpi
+            big_indarray(2,bigi,bigj) = ypi
+         end do
+      end do
+
+      do j = -inc,inc
+         do i = -inc,inc
+            xpi = indarray(1,i,j) 
+            ypi = indarray(2,i,j)
+            big_kdiff(xpi,ypi) = big_kdiff(xpi,ypi) + kdiff(i,j)
+         end do
+      end do
+
+      big_kdiff(0,0) = big_kdiff(user%gridsize,user%gridsize)
+      big_kdiff(user%gridsize + 1,user%gridsize + 1) = big_kdiff(1,1)
+      big_kdiff(0,:) = big_kdiff(user%gridsize,:)
+      big_kdiff(:,0) = big_kdiff(:,user%gridsize)
+      big_kdiff(user%gridsize + 1,:) = big_kdiff(1,:)
+      big_kdiff(:,user%gridsize + 1) = big_kdiff(:,1)
 
-   call util_diffusion_solver(user,surf,2 * inc + 1,indarray,kappat,cumulative_elchange,maxhits)
-   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 + cumulative_elchange(i,j)
-         surf(xpi,ypi)%ejcov = max(surf(xpi,ypi)%ejcov + cumulative_elchange(i,j),0.0_DP)
+      call util_diffusion_solver(user,surf,user%gridsize + 2,big_indarray,big_kdiff,big_cel,maxhits)
+      do j = 1,user%gridsize
+         do i = 1,user%gridsize
+            surf(i,j)%dem = surf(i,j)%dem + big_cel(i,j)
+            surf(i,j)%ejcov = max(surf(i,j)%ejcov + big_cel(i,j),0.0_DP)
+         end do
       end do
-   end do
-   deallocate(kappat,indarray,cumulative_elchange)
+   end if
+
+   deallocate(kdiff,indarray,cumulative_elchange)
 return
 end subroutine crater_soften