diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index 1229cd9f..21c96455 100644 --- a/src/crater/crater_populate.f90 +++ b/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 diff --git a/src/crater/crater_tally_observed.f90 b/src/crater/crater_tally_observed.f90 index f6b5a90d..af8211eb 100644 --- a/src/crater/crater_tally_observed.f90 +++ b/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)) diff --git a/src/crater/module_crater.f90 b/src/crater/module_crater.f90 index 00b7d167..d6cddf87 100644 --- a/src/crater/module_crater.f90 +++ b/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 diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index 6eae7272..51ed36c1 100644 --- a/src/globals/module_globals.f90 +++ b/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 diff --git a/src/io/io_read_surf.f90 b/src/io/io_read_surf.f90 index e47da102..3dbade81 100644 --- a/src/io/io_read_surf.f90 +++ b/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 diff --git a/src/io/io_write_surf.f90 b/src/io/io_write_surf.f90 index a8f7813a..9755ed02 100644 --- a/src/io/io_write_surf.f90 +++ b/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 diff --git a/src/io/io_write_tally.f90 b/src/io/io_write_tally.f90 index 64dff416..a8df6774 100644 --- a/src/io/io_write_tally.f90 +++ b/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 diff --git a/src/main/CTEM.f90 b/src/main/CTEM.f90 index 4f64f5c4..5d3ca6ea 100644 --- a/src/main/CTEM.f90 +++ b/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) diff --git a/src/util/util_add_to_layer.f90 b/src/util/util_add_to_layer.f90 index a99dd9dd..04ca287c 100644 --- a/src/util/util_add_to_layer.f90 +++ b/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' diff --git a/src/util/util_remove_from_layer.f90 b/src/util/util_remove_from_layer.f90 index b9e773b7..22bb26d4 100644 --- a/src/util/util_remove_from_layer.f90 +++ b/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 diff --git a/src/util/util_sort_layer.f90 b/src/util/util_sort_layer.f90 index a230c4c5..c3d055b1 100644 --- a/src/util/util_sort_layer.f90 +++ b/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