Fixed bug in crater_degradation_function call in subpixel_diffusion

.gitignore

@@ -31,3 +31,5 @@ __pycache__/
 *.DS_Store
 
 alocal.m4
+
+*.m4

src/Makefile.am

@@ -1,12 +1,12 @@
 bin_PROGRAMS = CTEM 
 #ifort optimized flags
-AM_FCFLAGS =  -O3 -qopenmp -parallel -xHost -ipo -assume byterecl -m64 -heap-arrays -FR 
+#AM_FCFLAGS =  -O3 -qopenmp -parallel -xHost -ipo -assume byterecl -m64 -heap-arrays -FR 
 #AM_FCFLAGS =  -O3 -openmp -parallel -assume byterecl -m64 -heap-arrays -FR 
 #ifort debug flags
 #AM_FCFLAGS = -O3 -p -g -openmp -debug all -traceback -CB -assume byterecl -m64 -heap-arrays -FR
 
 #gfortran optimized flags
-#AM_FCFLAGS =  -O3 -fopenmp -ffree-form  -g -fbounds-check -fbacktrace
+AM_FCFLAGS =  -O3 -fopenmp -ffree-form  -g -fbounds-check -fbacktrace
 #gfortran debug flags
 #AM_FCFLAGS =  -O0 -g -fopenmp -fbounds-check -Wall -Warray-bounds -Warray-temporaries -Wimplicit-interface -ffree-form
 

src/crater/crater_degradation_function.f90

@@ -16,20 +16,18 @@
 !**********************************************************************************************************************************
 
 
-function crater_degradation_function(user,crater) result(Kd)
+function crater_degradation_function(user,r) result(Kd)
    use module_globals
    use module_util
    use module_crater, EXCEPT_THIS_ONE => crater_degradation_function
    implicit none
 
   ! Arguments
    type(usertype),intent(in) :: user
-   type(cratertype),intent(in) :: crater
+   real(DP),intent(in) :: r
    real(DP) :: Kd
    real(DP) :: A,r_break,alpha_1,alpha_2,delta,r
 
-   r = crater%frad
-
    ! Testing values that match both mare scale and highlands scale
    alpha_1 = user%psi
    alpha_2 = user%psi2

src/crater/crater_soften.f90

@@ -45,7 +45,7 @@ subroutine crater_soften(user,surf,crater,domain)
 
 
    !kdiffmax = user%Kd1 * crater%frad**user%psi
-   kdiffmax = crater_degradation_function(user,crater)
+   kdiffmax = crater_degradation_function(user,crater%frad)
    inc = int(min(crater%fe * crater%frad / user%pix,SQRT2 * user%gridsize)) + 3 
    maxhits = (1 + (inc / (user%gridsize / 2)))**2
    crater%maxinc = max(crater%maxinc,inc)

src/crater/crater_subpixel_diffusion.f90

@@ -71,7 +71,7 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin)
    if (user%dosoftening) fe = crater%fe
 
    ! Generate both the subpixel and superdomain diffusive degradation
-   do k = 1,domain%pnum - 1
+   superloop: do k = 1,domain%pnum - 1
       dN = nflux(3,k) * user%interval * finterval
 
       diam_bedrock = nflux(1,k)
@@ -93,7 +93,7 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin)
          if (user%dosoftening) then 
          ! User-defined degradation function
             !dKdN = dKdN + user%Kd1 * PI * fe**2 * (radius)**(2.0_DP + user%psi) / domain%parea
-            dKdN = dKdN + PI * fe**2 * radius**2 * crater_degradation_function(user,crater) / domain%parea
+            dKdN = dKdN + PI * fe**2 * radius**2 * crater_degradation_function(user,radius) / domain%parea
          end if
          !Empirically-derived "intrinsic" degradation function from proximal ejecta redistribution
 
@@ -124,11 +124,11 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin)
       else if (user%dosoftening) then
       ! Do the degradation as individual circles
 
-         superlen = fd * diam + domain%side
+         superlen = fe * diam + domain%side
          cutout = 0.0_SP
          crater%continuous = RCONT * radius**(EXPCONT) 
          if (diam > domain%smallest_crater) then
-            if (.not.user%superdomain) exit
+            if (.not.user%superdomain) exit superloop
             ! Superdomain craters 
             cutout = real(domain%side + crater%continuous, kind=SP)
          end if
@@ -164,12 +164,9 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin)
             crater%continuous = RCONT * crater%frad**(EXPCONT) 
             crater%fe = user%fe
             krad = max(fd * crater%frad,crater%fe * crater%frad)
-
-            crater%fe = user%fe
-
 
             !dKdN = user%Kd1 * crater%frad**(user%psi)
-            dKdN = crater_degradation_function(user,crater)
+            dKdN = crater_degradation_function(user,crater%frad)
             inc  = int(krad / user%pix) + 2
             incsq = inc**2
 

src/crater/module_crater.f90

@@ -263,11 +263,10 @@ end function crater_critical_slope
    end interface
 
    interface
-      function crater_degradation_function(user,crater) result(Kd)
+      function crater_degradation_function(user,r) result(Kd)
       use module_globals
       type(usertype),intent(in) :: user
-      type(cratertype),intent(in) :: crater
-      real(DP) :: Kd
+      real(DP),intent(in) :: r
       end function crater_degradation_function
    end interface
 

src/ejecta/ejecta_emplace.f90

@@ -141,7 +141,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot)
 
    if (user%dosoftening) then
       !kdiffmax = user%Kd1 * crater%frad**(user%psi)
-      kdiffmax = crater_degradation_function(user,crater)
+      kdiffmax = crater_degradation_function(user,crater%frad)
    end if
 
    crater%maxinc = max(crater%maxinc,inc)