diff --git a/src/crater/crater_get_degradation_state.f90 b/src/crater/crater_get_degradation_state.f90 index 37356104..e94f293e 100644 --- a/src/crater/crater_get_degradation_state.f90 +++ b/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 diff --git a/src/crater/crater_tally_observed.f90 b/src/crater/crater_tally_observed.f90 index e163bc73..face1cba 100644 --- a/src/crater/crater_tally_observed.f90 +++ b/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 diff --git a/src/crater/module_crater.f90 b/src/crater/module_crater.f90 index 719b96b8..ee421e05 100644 --- a/src/crater/module_crater.f90 +++ b/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 diff --git a/src/io/io_write_tally.f90 b/src/io/io_write_tally.f90 index a8df6774..57654235 100644 --- a/src/io/io_write_tally.f90 +++ b/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) diff --git a/src/io/module_io.f90 b/src/io/module_io.f90 index 0dfcf6d8..5bcafa61 100644 --- a/src/io/module_io.f90 +++ b/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 diff --git a/src/main/CTEM.f90 b/src/main/CTEM.f90 index 5d3ca6ea..5858beb7 100644 --- a/src/main/CTEM.f90 +++ b/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'