Small bug fixes and now report the K value in the tally

src/crater/crater_get_degradation_state.f90

@@ -81,8 +81,8 @@ function crater_get_degradation_state(user,surf,crater,dd) result(Kval)
    rim = rim / nrim 
    bowl = bowl / nbowl 
    outer = outer / nouter
-   dd = (rim - bowl) / crater%fcrat
-   Kval = - 1._SP / Bconst * log(dd / ddinit)
+   dd = max((rim - bowl) / crater%fcrat,SMALLFAC)
+   Kval = - 1._DP / Bconst * log(real(dd,kind=DP) / ddinit)
    return
 
 end function crater_get_degradation_state

src/crater/crater_tally_observed.f90

@@ -18,7 +18,7 @@
 !  Notes       :  
 !
 !**********************************************************************************************************************************
-subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,oposlist,depthdiam)
+subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,oposlist,depthdiam,degradation_state)
    use module_globals
    use module_io
    use module_util
@@ -34,6 +34,7 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o
    real(DP),dimension(:),intent(out),allocatable,optional :: obslist
    real(SP),dimension(:,:),intent(out),allocatable,optional :: oposlist
    real(SP),dimension(:),intent(out),allocatable,optional :: depthdiam
+   real(DP),dimension(:),intent(out),allocatable,optional :: degradation_state
 
    ! Internal variables
    integer(I4B)                               :: i,j,layer,n,m,craternum,obstot
@@ -191,6 +192,7 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o
    ! Bin the observed craters if required 
    if (present(obsdist)) then
       allocate(depthdiam(onum))
+      allocate(degradation_state(onum))
       allocate(obslist(onum))
       allocate(oposlist(3,onum))
       ! Reset all the distribution bins
@@ -214,6 +216,7 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o
             obslist(obstot) = tlist(craternum)
             oposlist(:,obstot) = poslist(:,craternum)
             depthdiam(obstot) = tmp_depthdiam(craternum)
+            degradation_state(obstot)  = Kval(craternum)
             obstot = obstot - 1
          end if
       end do

src/crater/module_crater.f90

@@ -180,7 +180,7 @@ end subroutine crater_tally_true
    end interface
 
    interface
-      subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,oposlist,depthdiam)
+      subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,oposlist,depthdiam,degradation_state)
       use module_globals
       implicit none
       type(usertype),intent(in) :: user
@@ -191,6 +191,7 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o
       real(DP),dimension(:),intent(out),allocatable,optional :: obslist
       real(SP),dimension(:,:),intent(out),allocatable,optional :: oposlist
       real(SP),dimension(:),intent(out),allocatable,optional :: depthdiam
+      real(DP),dimension(:),intent(out),allocatable,optional :: degradation_state
       end subroutine crater_tally_observed
    end interface
 

src/io/io_write_tally.f90

@@ -18,7 +18,7 @@
 !**********************************************************************************************************************************
 
 
-subroutine io_write_tally(tdist,tlist,odist,olist,oposlist,depthdiam)
+subroutine io_write_tally(tdist,tlist,odist,olist,oposlist,depthdiam,degradation_state)
    use module_globals
    use module_io, EXCEPT_THIS_ONE => io_write_tally
    implicit none
@@ -28,6 +28,7 @@ subroutine io_write_tally(tdist,tlist,odist,olist,oposlist,depthdiam)
    real(DP),dimension(:),intent(in) :: olist
    real(SP),dimension(:,:),intent(in) :: oposlist
    real(SP),dimension(:),intent(in) :: depthdiam
+   real(DP),dimension(:),intent(in) :: degradation_state
 
    ! Internals
    integer(I4B)  :: i,distl,distc,ioerr,onum
@@ -42,7 +43,7 @@ subroutine io_write_tally(tdist,tlist,odist,olist,oposlist,depthdiam)
    distl = size(tdist,2)
    allocate(oldtdist(distc,distl))
 1000 format(3F16.4,2F15.0,F15.6)
-2000 format(ES23.15,1X,7(ES14.6,1X,:))
+2000 format(ES23.15,1X,8(ES14.6,1X,:))
    oldtdist = 0._DP
 
    inquire(file=tdistfile, exist=file_exists)
@@ -74,10 +75,10 @@ subroutine io_write_tally(tdist,tlist,odist,olist,oposlist,depthdiam)
    close(LUN)
 
    open(LUN, FILE=OLISTFILE, status='REPLACE')
-   write(LUN,'("#Dcrat(m)                 xpos(m)        ypos(m)        time(y)        depth/diam")'  )
+   write(LUN,'("#Dcrat(m)                 xpos(m)        ypos(m)        time(y)        depth/diam        deg_state(m^2)")' )
    do i=1,onum
 
-      write(LUN,2000) olist(i),oposlist(:,i),depthdiam(i)
+      write(LUN,2000) olist(i),oposlist(:,i),depthdiam(i),degradation_state(i)
    end do
    close(LUN)
 

src/io/module_io.f90

@@ -142,13 +142,14 @@ end subroutine io_write_dist
    end interface
 
    interface
-      subroutine io_write_tally(tdist,tlist,odist,olist,oposlist,depthdiam)
+      subroutine io_write_tally(tdist,tlist,odist,olist,oposlist,depthdiam,degradation_state)
       use module_globals
       implicit none
       real(DP),dimension(:,:),intent(in) :: tdist,tlist,odist
       real(DP),dimension(:),intent(in) :: olist
       real(SP),dimension(:,:),intent(in) :: oposlist
       real(SP),dimension(:),intent(in) :: depthdiam
+      real(DP),dimension(:),intent(in) :: degradation_state
       end subroutine io_write_tally
    end interface
 

src/main/CTEM.f90

@@ -41,6 +41,7 @@ program CTEM
 real(DP),dimension(:,:),allocatable  :: prod,vdist,pdist,crtscl,truedist,obsdist,truelist
 real(DP),dimension(:),allocatable :: obslist
 real(SP),dimension(:),allocatable :: depthdiam
+real(DP),dimension(:),allocatable :: degradation_state
 real(SP),dimension(:,:),allocatable :: oposlist
 integer(I8B),dimension(:),allocatable :: production_list
 ! Miscellaneous variables
@@ -130,10 +131,10 @@ program CTEM
    write(*,*) "Surface-affecting craters generated:   ",ntrue
    write(*,*) "Visible craters generated:             ",vistrue
 end if
-call crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,oposlist,depthdiam)
+call crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,oposlist,depthdiam,degradation_state)
 ntotkilled = ntotkilled + nkilled
 write(*,*) 'Craters killed during tally: ',ntotkilled
-call io_write_tally(truedist,truelist(:,1:ntrue),obsdist,obslist,oposlist,depthdiam)
+call io_write_tally(truedist,truelist(:,1:ntrue),obsdist,obslist,oposlist,depthdiam,degradation_state)
 if (.not.user%tallyonly) then
    write(*,*) "Writing surface files"
    call io_write_surf(user,surf)
@@ -166,7 +167,7 @@ program CTEM
 ! Deallocate all the allocatables
 deallocate(seedarr)
 deallocate(surf,prod,vdist,pdist,crtscl,truedist,truelist,obsdist,obslist,nflux,production_list)
-deallocate(oposlist,depthdiam)
+deallocate(oposlist,depthdiam,degradation_state)
 
 !$ t2 = omp_get_wtime()
 !$ write(*,*) 'Timing information'