Re-wrote tally to explicitly look at degradation state to conform to the equations in Minton et al. (2019)

src/Makefile.am

@@ -6,9 +6,9 @@ bin_PROGRAMS = CTEM
 #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
+AM_FCFLAGS =  -O0 -g -fopenmp -fbounds-check -Wall -Warray-bounds -Warray-temporaries -Wimplicit-interface -ffree-form
 
 CTEM_SOURCES = globals/module_globals.f90\
 util/module_util.f90\
@@ -82,6 +82,8 @@ crater/crater_subpixel_diffusion.f90\
 crater/crater_make_list.f90\
 crater/crater_critical_slope.f90\
 crater/crater_degradation_function.f90\
+crater/crater_visibility.f90\
+crater/crater_get_degradation_state.f90\
 init/init_domain.f90\
 init/init_dist.f90\
 init/init_surf.f90\

src/crater/crater_degradation_function.f90

@@ -36,8 +36,8 @@ function crater_degradation_function(user,r) result(Kd)
    delta = 1.0_DP
 
    A = user%Kd1 * r_break**(alpha_1) / (1_DP + (1.d0 / delta) * (alpha_1 - alpha_2) / alpha_1)**2 
-
    Kd =  A * (r / r_break)**(alpha_1) * (0.5_DP * (1.0_DP + (r / r_break)**(1.0_DP / delta)))**((alpha_2 - alpha_1 ) / delta)
+
    return
 
 end function crater_degradation_function

src/crater/crater_get_degradation_state.f90

@@ -0,0 +1,89 @@
+!  Project     : CTEM
+!  Language    : Fortran 2003
+!
+!  Description : Determines the current degradation state of the crater
+!  
+!
+!  Input
+!    Arguments : 
+!
+!  Output
+!    Arguments :
+!           
+! 
+!  Notes       :  
+!
+!**********************************************************************************************************************************
+
+
+function crater_get_degradation_state(user,surf,crater,dd) result(Kval)
+   use module_globals
+   use module_util
+   use module_crater, EXCEPT_THIS_ONE => crater_get_degradation_state
+   implicit none
+
+   ! Arguments
+   type(usertype),intent(in) :: user
+   type(surftype),dimension(:,:),intent(in) :: surf
+   type(cratertype),intent(in) :: crater
+   real(SP),intent(out) :: dd
+   ! Return value
+   real(DP) :: Kval
+
+   ! Internal variables
+   integer(I4B) :: inc,xpi,ypi,i,j
+   integer(I4B) :: nrim,nbowl,nouter
+   real(DP) :: rim,bowl,outer,rad,baseline,ddinit
+
+   ! Definitions of crater dimensions
+   real(DP),parameter :: OCUTOFF = 5.5e-2_DP
+   real(DP),parameter :: RIMDI = 1.0_DP
+   real(DP),parameter :: RIMDO = 1.2_DP
+   real(DP),parameter :: BOWLD = 0.2_DP
+   real(DP),parameter :: OUTERD = 2.0_DP
+
+   ! Fit for complex crater diffusion model
+   real(DP),parameter :: Bconst = 8.0_DP
+   real(DP),parameter :: dDpikeC = 1.725_DP ! CTEM version of the Pike (1977) depth-to-diam constant
+   real(DP),parameter :: ddsimple = 0.21362307923564000_DP ! CTEM initial depth-to-diameter of simple craterws
+
+   ddinit = min(ddsimple, dDpikeC * (crater%fcrat * 1e-3_DP)**(0.301_DP - 1._DP))
+
+   inc = int(OUTERD * crater%frad / user%pix) + 1 
+
+
+   nrim = 0
+   nbowl = 0
+   nouter = 0
+   bowl = 0.0_DP
+   rim = 0.0_DP
+   outer = 0.0_DP
+
+   do j = -inc, inc
+      do i = -inc, inc
+         rad = sqrt((i**2 + j**2)*1._DP) * user%pix / crater%frad
+         baseline = crater%melev + ((i * crater%xslp) + (j * crater%yslp)) * user%pix 
+         xpi = crater%xlpx + i
+         ypi = crater%ylpx + j
+         call util_periodic(xpi,ypi,user%gridsize)
+         if (rad <= BOWLD) then
+            bowl = bowl + surf(xpi,ypi)%dem - baseline
+            nbowl = nbowl + 1 
+         else if ((rad >= RIMDI).and.(rad < RIMDO)) then
+            rim = rim + surf(xpi,ypi)%dem - baseline
+            nrim = nrim + 1 
+         else if (rad > RIMDO) then
+            outer = outer + surf(xpi,ypi)%dem - baseline
+            nouter = nouter + 1 
+         end if
+      end do
+   end do
+   rim = rim / nrim 
+   bowl = bowl / nbowl 
+   outer = outer / nouter
+   dd = (rim - bowl) / crater%fcrat
+   Kval = - 1._SP / Bconst * log(dd / ddinit)
+   return
+
+end function crater_get_degradation_state
+

src/crater/crater_tally_observed.f90

@@ -45,6 +45,7 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o
    real(SP),dimension(:,:),allocatable :: mposlist
    real(SP),dimension(:),allocatable :: tmp_depthdiam
    logical,dimension(:),allocatable :: countable
+   real(DP),dimension(:),allocatable :: Kval
    integer(I4B),dimension(:,:),allocatable:: mpxlist
    integer(I4B),dimension(:),allocatable :: ind,mlayerlist
    integer(I4B) :: q,mnum,imnum,maxpix,npix,ntot
@@ -55,17 +56,7 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o
    real(DP),dimension(:),allocatable :: dis,elev
    integer(I4B),dimension(:),allocatable :: totpix,istart,iend
    integer(I4B) :: tnum ! True number and observed number
-   integer(I4B) :: inc,xpi,ypi
-   logical :: killable
-   integer(I4B) :: nrim,nbowl,nouter
-   real(DP) :: rim,bowl,outer,rad,baseline,dd
-   ! Counting parameters from Howl study
-   real(DP),parameter :: DDCUTOFF = 5.0e-2_DP
-   real(DP),parameter :: OCUTOFF = 5.5e-2_DP
-   real(DP),parameter :: RIMDI = 1.0_DP
-   real(DP),parameter :: RIMDO = 1.2_DP
-   real(DP),parameter :: BOWLD = 0.2_DP
-   real(DP),parameter :: OUTERD = 2.0_DP
+
 
    ! Executable code
 
@@ -148,17 +139,18 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o
    allocate(tlist(tnum))
    allocate(tmp_depthdiam(tnum))
    allocate(countable(tnum))
+   allocate(Kval(tnum))
    countable = .false.
    nkilled = 0
    onum = 0
+   ! Here we go through the list of craters and take the measures of each one to
+   ! be used by the tally subroutine in the next pass
    !$OMP PARALLEL DEFAULT(PRIVATE) IF(tnum > INCPAR) &
    !$OMP SHARED(user,surf) &
    !$OMP SHARED(domain,maxpix,istart,iend,ind,mlist,mposlist,mpxlist,mlayerlist) &
-   !$OMP SHARED(tnum,totpix,poslist,tlist,tmp_depthdiam,countable) &
+   !$OMP SHARED(tnum,totpix,poslist,tlist,tmp_depthdiam,countable,Kval) &
    !$OMP REDUCTION(+:nkilled) &
    !$OMP REDUCTION(+:onum) 
-   ! Here we go through the list of craters and take the measures of each one to
-   ! be used by the tally subroutine in the next pass
    !$OMP DO 
    do craternum = 1,tnum
       ! This is the first pixel of this crater, so record the appropriate values
@@ -170,71 +162,28 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o
       crater%fradpx = int(crater%frad / user%pix) + 1
       crater%xlpx = int(poslist(1,craternum) / user%pix)
       crater%ylpx = int(poslist(2,craternum) / user%pix)
-      call crater_averages(user,surf,crater)
-      nrim = 0
-      nbowl = 0
-      nouter = 0
-      bowl = 0.0_DP
-      rim = 0.0_DP
-      outer = 0.0_DP
-      inc = int(OUTERD * crater%frad / user%pix) + 1 
-      do j = -inc, inc
-         do i = -inc, inc
-            rad = sqrt((i**2 + j**2)*1._DP) * user%pix / crater%frad
-            baseline = crater%melev + ((i * crater%xslp) + (j * crater%yslp)) * user%pix 
-            xpi = crater%xlpx + i
-            ypi = crater%ylpx + j
-            call util_periodic(xpi,ypi,user%gridsize)
-            if (rad <= BOWLD) then
-               bowl = bowl + surf(xpi,ypi)%dem - baseline
-               nbowl = nbowl + 1 
-            else if ((rad >= RIMDI).and.(rad < RIMDO)) then
-               rim = rim + surf(xpi,ypi)%dem - baseline
-               nrim = nrim + 1 
-            else if (rad > RIMDO) then
-               outer = outer + surf(xpi,ypi)%dem - baseline
-               nouter = nouter + 1 
-            end if
-         end do
-      end do
-      rim = rim / nrim 
-      bowl = bowl / nbowl 
-      outer = outer / nouter
-      tmp_depthdiam(craternum) = (rim - bowl) / crater%fcrat
 
-      !if (crater%fcrat < 20e3_DP) then
-      !   dd = DDCUTOFF
-      !else
-      !   dd = DDCUTOFF - (crater%fcrat - 20e3_DP) * 2e-7
-      !end if
-      dd = min(DDCUTOFF, 24.897 * crater%fcrat ** (-0.632545))
+      call crater_averages(user,surf,crater)
 
-      if (((tmp_depthdiam(craternum) > dd).and.((outer - rim) / crater%fcrat < OCUTOFF)).and.&
-          (nrim /= 0).and.(nbowl /= 0).and.(nouter /= 0)) then
-         countable(craternum) = .true.
-         killable = .false.
-      else  
-         countable(craternum) = .false.
-         killable = .true.
-      end if   
       if ((crater%fcrat < domain%smallest_counted_crater) .or. &
          (totpix(craternum) < int(0.1_DP * 0.25_DP * PI * crater%fcrat / user%pix))) then
          countable(craternum) = .false.
-         killable = .true.
+      else 
+         Kval(craternum) = crater_get_degradation_state(user,surf,crater,tmp_depthdiam(craternum))
+         countable(craternum) = crater_visibility(user,crater,Kval(craternum))
       end if
-      ! TESTING
-      if (user%testtally) then
-         countable = .true.
-         killable = .false.
+
+      if (user%testtally) then ! used for testing the tally system
+         countable(craternum) = .true.
       end if
-      if (killable) then ! Obliterate the crater from the record
+
+      if (.not.countable(craternum)) then ! Obliterate the crater from the record
          do i=istart(craternum),iend(craternum)
             call util_remove_from_layer(surf(mpxlist(1,ind(i)),mpxlist(2,ind(i))),mlayerlist(ind(i)))
          end do
          nkilled = nkilled + 1
       end if
       if (countable(craternum)) onum = onum + 1
-      !deallocate(csurf)
    end do
    !$OMP END DO
    !$OMP END PARALLEL
@@ -293,6 +242,7 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o
    deallocate(tmp_depthdiam)
    deallocate(ind)
    deallocate(countable)
+   deallocate(Kval)
 
    return
 end subroutine crater_tally_observed

src/crater/crater_visibility.f90

@@ -0,0 +1,50 @@
+!  Project     : CTEM
+!  Language    : Fortran 2003
+!
+!  Description : Calculates the visibility function from Minton et al. (2019)
+!  
+!
+!  Input
+!    Arguments : 
+!
+!  Output
+!    Arguments :
+!           
+! 
+!  Notes       :  
+!
+!**********************************************************************************************************************************
+
+
+function crater_visibility(user,crater,Kval) result(iscountable)
+   use module_globals
+   use module_util
+   use module_crater, EXCEPT_THIS_ONE => crater_visibility
+   implicit none
+
+   ! Arguments
+   type(usertype),intent(in) :: user
+   type(cratertype),intent(in) :: crater
+   real(DP),intent(in) :: Kval
+
+   ! Result variable
+   logical :: iscountable 
+
+   ! Internal variables
+   real(DP) :: Kv
+   real(DP),parameter :: Kv1 = 0.17
+   real(DP),parameter :: gam = 2.0
+
+   ! Counting parameter for simple craters from Bryan Howl's study in Minton et al. (2019)
+
+   Kv = Kv1 * (crater%frad)**gam
+   if (Kval >= kv) then
+      iscountable = .false.
+   else
+      iscountable = .true.
+   end if
+
+   return
+
+end function crater_visibility
+

src/crater/module_crater.f90

@@ -271,6 +271,32 @@ function crater_degradation_function(user,r) result(Kd)
       end function crater_degradation_function
    end interface
 
+   interface
+      function crater_get_degradation_state(user,surf,crater,dd) result(Kval)
+      use module_globals
+      implicit none
+      ! Arguments
+      type(usertype),intent(in) :: user
+      type(surftype),dimension(:,:),intent(in) :: surf
+      type(cratertype),intent(in) :: crater
+      real(SP),intent(out) :: dd
+      ! Return value
+      real(DP) :: Kval
+      end function crater_get_degradation_state
+   end interface
+
+   interface
+      function crater_visibility(user,crater,Kval) result(iscountable)
+      use module_globals
+      implicit none
+      ! Arguments
+      type(usertype),intent(in) :: user
+      type(cratertype),intent(in) :: crater
+      real(DP),intent(in) :: Kval
+      ! Result variable
+      logical :: iscountable       
+      end function crater_visibility
+   end interface
 
 
 end module