Added timestamp feature. Simulation time of formation now saved.

src/crater/crater_populate.f90

@@ -18,7 +18,7 @@
 !
 !**********************************************************************************************************************************
 subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,ntrue,vistrue,ntotkilled,truelist,mass, &
-                           fracdone,nflux,ntotcrat)
+                           fracdone,nflux,ntotcrat,curyear)
    use module_globals
    use module_seismic
    use module_io
@@ -44,6 +44,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
    real(DP),intent(out)                            :: fracdone
    real(DP),dimension(:,:),intent(in)              :: nflux 
    integer(I8B),intent(in)                         :: ntotcrat  ! Total number of attempted impacts
+   real(DP),intent(in)                             :: curyear
 
    ! Internal variables
    real(DP)                :: cmin     ! Minimum crater diameter (m)
@@ -160,7 +161,8 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
       if (crater%fcrat < domain%smallest_crater) makecrater = .false.
 
       if (makecrater) then
-
+         ! Stamp the current time onto the crater
+         crater%timestamp = real(curyear + real(icrater,kind=DP) / real(ntotcrat,kind=DP) * user%interval,kind=SP)
          ! Find the visible crater parameters
          call crater_find_visible(user,crater,domain)
 
@@ -181,6 +183,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
          truelist(4,ntrue) = crater%yl
          truelist(5,ntrue) = crater%impvel
          truelist(6,ntrue) = crater%sinimpang
+         truelist(7,ntrue) = crater%timestamp
          mass = mass + crater%impmass
 
          crater%maxinc = 0

src/crater/crater_tally_observed.f90

@@ -88,7 +88,7 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o
 
 
    allocate(mlist(max(mnum,1)))
-   allocate(mposlist(2,max(mnum,1)))
+   allocate(mposlist(3,max(mnum,1)))
    allocate(mpxlist(2,max(mnum,1)))
    allocate(mlayerlist(max(mnum,1)))
    allocate(ind(max(mnum,1)))
@@ -111,6 +111,7 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o
                ! Save the crater center coordinates
                mposlist(1,imnum) = surf(m,n)%xl(layer) 
                mposlist(2,imnum) = surf(m,n)%yl(layer) 
+               mposlist(3,imnum) = surf(m,n)%timestamp(layer) 
                ! Save the pixel index location
                mpxlist(1,imnum) = m 
                mpxlist(2,imnum) = n 
@@ -143,7 +144,7 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o
    iend(tnum) = mnum
 
    allocate(totpix(tnum))
-   allocate(poslist(2,tnum))
+   allocate(poslist(3,tnum))
    allocate(tlist(tnum))
    allocate(tmp_depthdiam(tnum))
    allocate(countable(tnum))
@@ -242,7 +243,7 @@ subroutine crater_tally_observed(user,surf,domain,nkilled,onum,obsdist,obslist,o
    if (present(obsdist)) then
       allocate(depthdiam(onum))
       allocate(obslist(onum))
-      allocate(oposlist(2,onum))
+      allocate(oposlist(3,onum))
       ! Reset all the distribution bins
       do i=1,domain%distl
          obsdist(1,i) = 1e3_DP*SQRT2**(domain%plo+(i-1))

src/crater/module_crater.f90

@@ -30,7 +30,7 @@ end subroutine crater_mass_conservation
 
    interface
       subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,ntrue,vistrue,ntotkilled,truelist,&
-                                 mass,fracdone,nflux,ntotcrat)
+                                 mass,fracdone,nflux,ntotcrat,curyear)
       use module_globals
       implicit none
       type(usertype),intent(in) :: user
@@ -47,6 +47,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
       real(DP),intent(out)                         :: fracdone
       real(DP),dimension(:,:),intent(in)           :: nflux 
       integer(I8B),intent(in)                      :: ntotcrat
+      real(DP),intent(in)                          :: curyear
       end subroutine crater_populate
    end interface
 

src/globals/module_globals.f90

@@ -79,6 +79,7 @@ module module_globals
 type surftype
    real(DP),dimension(MAXLAYER) :: diam   ! Crater diameter
    real(SP),dimension(MAXLAYER) :: xl,yl ! Crater center 
+   real(SP),dimension(MAXLAYER) :: timestamp ! Crater formation time
    real(DP) :: ejcov                ! Ejecta coverage
    real(DP) :: dem                  ! Digital elevation model
    real(DP) :: mantle               ! Height of mantle (should be smaller than dem)
@@ -91,6 +92,7 @@ module module_globals
    real(DP) :: imp,imprad,impvel,sinimpang,impmass ! Impactor properties
    ! Real domain properties
    real(SP) :: xl,yl               ! Crater center in simulated surface size units 
+   real(SP) :: timestamp           ! The formation time of the crater
    real(DP) :: rad                 ! Transient radius
    real(DP) :: grad                ! Strengthless material transient crater radius
    real(DP) :: frad                ! Final crater radius
@@ -249,6 +251,7 @@ module module_globals
 
 ! Grid array file names
 character(*),parameter :: DIAMFILE   = 'surface_diam.dat'
+character(*),parameter :: TIMEFILE   = 'surface_time.dat'
 character(*),parameter :: EJCOVFILE  = 'surface_ejc.dat'
 character(*),parameter :: DEMFILE    = 'surface_dem.dat'
 character(*),parameter :: REGOFILE   = 'surface_regotop.dat'
@@ -274,7 +277,7 @@ module module_globals
 integer(I4B),parameter :: PBCLIM = 1             ! periodic boundary condition limit
 integer(I4B),parameter :: SMALLESTCOUNTABLE = 5 ! Minimum pixel diameter for a crater to be considered countable
 real(DP),parameter     :: SMALLESTEJECTA = 1.5  ! Minimum number of pixels from center of crater for an ejecta to have any surface effects
-integer(I4B),parameter :: TRUECOLS = 6 ! Number of columns in the true crater count array
+integer(I4B),parameter :: TRUECOLS = 7 ! Number of columns in the true crater count array
 integer(I4B)           :: NTHREADS = 1 ! Number of OpenMP threads (reset by OpenMP if a parallel environment is detected)
 integer(I4B),parameter :: INCPAR = 1   ! Minimum size of inc variables before parallelization kicks in
 

src/io/io_read_surf.f90

@@ -48,6 +48,7 @@ subroutine io_read_surf(user,surf)
    read(LUN,rec=1) surf%dem
    close(LUN)
 
+
    open(LUN,file=EJCOVFILE,status='old',form='unformatted',recl=recsize,access='direct',iostat=ioerr)
    if (ioerr/=0) then
       write(*,*) 'Error! Cannot read file ',trim(adjustl(EJCOVFILE))
@@ -66,6 +67,17 @@ subroutine io_read_surf(user,surf)
    end do
    close(LUN)
 
+   open(LUN,file=TIMEFILE,status='old',form='unformatted',recl=recsize,access='direct',iostat=ioerr)
+   if (ioerr/=0) then
+      write(*,*) 'Error! Cannot read file ',trim(adjustl(DIAMFILE))
+      stop
+   end if
+   do i=1,user%numlayers 
+      read(LUN,rec=i) surf%timestamp(i)
+   end do
+   close(LUN)
+
+
    recsize = sizeof(stmp) * user%gridsize * user%gridsize
    open(LUN,file=POSFILE,status='old',form='unformatted',recl=recsize,access='direct',iostat=ioerr)
    if (ioerr/=0) then

src/io/io_write_surf.f90

@@ -53,6 +53,12 @@ subroutine io_write_surf(user,surf)
    end do
    close(LUN)
 
+   open(LUN,file=TIMEFILE,status='replace',form='unformatted',recl=recsize,access='direct')
+   do i=1,user%numlayers 
+      write(LUN,rec=i) surf%timestamp(i)
+   end do
+   close(LUN)
+
    recsize = sizeof(stmp) * user%gridsize * user%gridsize
    open(LUN,file=POSFILE,status='replace',form='unformatted',recl=recsize,access='direct')
    do i=1,user%numlayers 

src/io/io_write_tally.f90

@@ -42,7 +42,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,6(ES14.6,1X,:))
+2000 format(ES23.15,1X,7(ES14.6,1X,:))
    oldtdist = 0._DP
 
    inquire(file=tdistfile, exist=file_exists)
@@ -74,7 +74,7 @@ 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)        depth/diam")'  )
+   write(LUN,'("#Dcrat(m)                 xpos(m)        ypos(m)        time(y)        depth/diam")'  )
    do i=1,onum
 
       write(LUN,2000) olist(i),oposlist(:,i),depthdiam(i)
@@ -83,7 +83,7 @@ subroutine io_write_tally(tdist,tlist,odist,olist,oposlist,depthdiam)
 
 
    open(LUN, FILE=TLISTFILE, status='REPLACE')
-   write(LUN,'("#Dcrat(m)                 Dimp(m)        xpos(m)        ypos(m)        vimp(m/s)      sinang")')
+   write(LUN,'("#Dcrat(m)                 Dimp(m)        xpos(m)        ypos(m)        vimp(m/s)      sinang         time(y)")')
    do i=1,size(tlist,2)
       write(LUN,2000) tlist(:,i)
    end do

src/main/CTEM.f90

@@ -48,7 +48,7 @@ program CTEM
 logical                 :: restart  ! F = new run (start with a fresh surface)
 integer(I8B)            :: totalimpacts ! Total number of impacts ever produced 
 integer(I4B)            :: ncount   ! Current count in ctem_driver IDL run
-integer(I4B)            :: n, xp, yp, i        ! Size of random number generator seed array
+integer(I4B)            :: n, xp, yp        ! Size of random number generator seed array
 integer(I4B),dimension(:),allocatable :: seedarr ! Random number generator seed array
 real(DP)                :: curyear
 real(DP)                :: mass
@@ -60,7 +60,6 @@ program CTEM
 integer(I4B)            :: ntotkilled 
 integer(I8B)            :: ntotcrat
 integer(I4B)            :: onum
-real(DP)                :: lambda
 !$ real(DP)             :: t1,t2
 real(DP),dimension(:,:),allocatable :: nflux
 
@@ -116,7 +115,7 @@ program CTEM
       call crater_make_list(domain,prod,ntotcrat,production_list)
    end if
    call crater_populate(user,surf,crater,domain,prod,production_list,vdist,ntrue,vistrue,ntotkilled,truelist,mass,&
-                        fracdone,nflux,ntotcrat)
+                        fracdone,nflux,ntotcrat,curyear)
 
    ! Get the last seed and save it to file
    call random_seed(get=seedarr)

src/util/util_add_to_layer.f90

@@ -47,6 +47,7 @@ subroutine util_add_to_layer(user,surfi,crater)
    surfi%diam(layer) = crater%fcrat
    surfi%xl(layer) = crater%xl
    surfi%yl(layer) = crater%yl
+   surfi%timestamp(layer) = crater%timestamp
 else
    write(*,*)
    write(*,*) 'WARNING! No free layer to add crater pixel. Consider increasing NUMLAYERS'

src/util/util_remove_from_layer.f90

@@ -32,5 +32,6 @@ subroutine util_remove_from_layer(surfi,layer)
 surfi%diam(layer) = 0.0_DP
 surfi%xl(layer) = 0.0_SP
 surfi%yl(layer) = 0.0_SP
+surfi%timestamp(layer) = 0.0_SP
 
 end subroutine util_remove_from_layer

src/util/util_sort_layer.f90

@@ -31,7 +31,7 @@ subroutine util_sort_layer(user,surf,crater)
    ! Internal variables
    integer(I4B),dimension(user%numlayers) :: ind
    real(DP),dimension(user%numlayers) :: tempdiam
-   real(SP),dimension(user%numlayers) :: tempxpos,tempypos
+   real(SP),dimension(user%numlayers) :: tempxpos,tempypos,temptime
    integer(I4B) :: i,j,k,inc,incsq,mx,my,iradsq
 
    inc = min(crater%rimdispx,(user%gridsize - 1)/2)
@@ -54,6 +54,7 @@ subroutine util_sort_layer(user,surf,crater)
             tempdiam = surf(mx,my)%diam(1:user%numlayers)
             tempxpos = surf(mx,my)%xl(1:user%numlayers)
             tempypos = surf(mx,my)%yl(1:user%numlayers)
+            temptime = surf(mx,my)%timestamp(1:user%numlayers)
 
             ! Sort the layers by crater diameter
             call util_mrgrnk(surf(mx,my)%diam(1:user%numlayers),ind)
@@ -62,6 +63,7 @@ subroutine util_sort_layer(user,surf,crater)
                surf(mx,my)%diam(k)=tempdiam(ind(k))
                surf(mx,my)%xl(k)=tempxpos(ind(k))
                surf(mx,my)%yl(k)=tempypos(ind(k))
+               surf(mx,my)%timestamp(k)=temptime(ind(k))
             end do
          end if
       end do