diff --git a/examples/global-lunar-bombardment/craterlist.in b/examples/global-lunar-bombardment/craterlist.in index d7d21187..6c9f00e6 100644 --- a/examples/global-lunar-bombardment/craterlist.in +++ b/examples/global-lunar-bombardment/craterlist.in @@ -1,3 +1,3 @@ # imp vel ang xoff yoff t_Ga -90856.5 15.0e3 90.0 -5.17529e5 1.037286e6 2.00 -18095.0 15.0e3 90.0 -9.38629e5 1.3466084e6 1.00 \ No newline at end of file +90856.5 15.0e3 90.0 -5.17529e5 1.037286e6 0.4 +18095.0 15.0e3 90.0 -9.38629e5 1.3466084e6 0.2 \ No newline at end of file diff --git a/examples/global-lunar-bombardment/ctem.in b/examples/global-lunar-bombardment/ctem.in index c97d8488..986b5672 100755 --- a/examples/global-lunar-bombardment/ctem.in +++ b/examples/global-lunar-bombardment/ctem.in @@ -16,7 +16,7 @@ realcraterlist craterlist.in ! list of 'real' craters for Quas ! IDL driver in uts -interval 4.0 +interval 0.6 numintervals 1 ! Total number of intervals (total time = interval * numintervals) <--when runtype is 'single' restart F ! Restart a previous run impfile NPFextrap.dat ! Impactor SFD rate file (col 1: Dimp (m), col 2: ! impactors > D (m**(-2) y**(-1)) diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index 9858732c..b1863ebe 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,curyear) + fracdone,nflux,ntotcrat,curyear,rclist) use module_globals use module_seismic use module_io @@ -45,6 +45,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt real(DP),dimension(:,:),intent(in) :: nflux integer(I8B),intent(in) :: ntotcrat ! Total number of attempted impacts real(DP),intent(in) :: curyear + real(DP),dimension(:,:), intent(in) :: rclist !array of 'real' craters for quasiMC ! Internal variables real(DP) :: cmin ! Minimum crater diameter (m) @@ -100,10 +101,8 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt ! read initial quasi-MC position if (user%doquasimc) then - call io_read_craterlist(user,domain) rccount = 1 - write(*,*) user%rctime - write(*,*) domain%rcnum + user%rctime = rclist(6,rccount) end if @@ -145,22 +144,31 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt !if in quasiMC mode: check to see if it's time for a real crater if (user%doquasimc) then if (crater%timestamp > user%rctime) then - write(*,*) "real crater @ ", crater%timestamp - ! add the code to emplace the test crater here - call io_read_craterlist(user, domain) + write(*,*) "real crater at time ", rclist(6, rccount) + user%testflag = .true. + user%testimp = rclist(1, rccount) + user%testvel = rclist(2, rccount) + user%testang = rclist(3, rccount) + user%testxoffset = rclist(4, rccount) + user%testyoffset = rclist(5, rccount) + end if + end if + ! generate random crater + call crater_generate(user,crater,domain,prod,production_list,vdist,surf) + if (user%testflag) write(*,*) 'Dcrat = ',crater%fcrat + if (user%testflag) write(*,*) 'Dtrans = ',crater%rad*2 + if (user%doquasimc) then + if (crater%timestamp > user%rctime) then + user%testflag = .false. rccount = rccount + 1 - write(*,*) user%rctime - write(*,*) rccount if (rccount > domain%rcnum) then write(*,*) "real crater list complete." user%rctime = 1e30 + else + user%rctime = rclist(6,rccount) end if end if end if - ! generate random crater - call crater_generate(user,crater,domain,prod,production_list,vdist,surf) - if (user%testflag) write(*,*) 'Dcrat = ',crater%fcrat - if (user%testflag) write(*,*) 'Dtrans = ',crater%rad*2 if (crater%fcrat > domain%biggest_crater) then ! End the run if the crater is too big if (user%killatmaxcrater) then fracdone = real(icrater,kind=DP) / real(ntotcrat,kind=DP) diff --git a/src/crater/module_crater.f90 b/src/crater/module_crater.f90 index ede92236..654e7485 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,curyear) + mass,fracdone,nflux,ntotcrat,curyear, rclist) use module_globals implicit none type(usertype),intent(in) :: user @@ -48,6 +48,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt real(DP),dimension(:,:),intent(in) :: nflux integer(I8B),intent(in) :: ntotcrat real(DP),intent(in) :: curyear + real(DP), dimension(:,:), INTENT(in) :: rclist end subroutine crater_populate end interface diff --git a/src/io/io_read_craterlist.f90 b/src/io/io_read_craterlist.f90 index 86f5880c..ec17ff58 100644 --- a/src/io/io_read_craterlist.f90 +++ b/src/io/io_read_craterlist.f90 @@ -16,18 +16,20 @@ ! Notes : only reads needed line [WIP] ! !********************************************************************************************************************************** -subroutine io_read_craterlist(user, domain) +subroutine io_read_craterlist(rclist, user, domain) use module_globals use module_io, EXCEPT_THIS_ONE => io_read_craterlist implicit none ! Arguments + real(DP),dimension(:,:),intent(out) :: rclist type(usertype),intent(inout) :: user type(domaintype),intent(in) :: domain ! Internals integer(I4B) :: i,ierr integer(I4B),parameter :: LUN=7 + real(DP) :: dumm ! Executable code @@ -39,15 +41,15 @@ subroutine io_read_craterlist(user, domain) stop end if + do i=1,domain%rcnum + read(LUN,*,iostat=ierr) rclist(1:6,i) + if (ierr/=0) then + write(*,*) "Unable to read file ",trim(rcfile) + stop + end if + end do - !do i=0,rccount - !read(LUN,*) - read(LUN,*,iostat=ierr) user%testimp, user%testvel, user%testang, user%testxoffset, user%testyoffset, user%rctime - if (ierr/=0) then - write(*,*) "Unable to read file ",trim(rcfile) - stop - end if - !end do + close(LUN) end subroutine io_read_craterlist diff --git a/src/io/module_io.f90 b/src/io/module_io.f90 index fc1d69e3..edb1548c 100644 --- a/src/io/module_io.f90 +++ b/src/io/module_io.f90 @@ -122,9 +122,10 @@ end subroutine io_read_prod end interface interface - subroutine io_read_craterlist(user,domain) + subroutine io_read_craterlist(rclist, user,domain) use module_globals implicit none + real(DP),dimension(:,:),intent(out) :: rclist type(usertype),intent(in) :: user type(domaintype),intent(in) :: domain end subroutine io_read_craterlist diff --git a/src/main/CTEM.f90 b/src/main/CTEM.f90 index 5d3ca6ea..dc233529 100644 --- a/src/main/CTEM.f90 +++ b/src/main/CTEM.f90 @@ -38,7 +38,7 @@ program CTEM type(cratertype) :: crater ! Distribution arrays -real(DP),dimension(:,:),allocatable :: prod,vdist,pdist,crtscl,truedist,obsdist,truelist +real(DP),dimension(:,:),allocatable :: prod,vdist,pdist,crtscl,truedist,obsdist,truelist,rclist real(DP),dimension(:),allocatable :: obslist real(SP),dimension(:),allocatable :: depthdiam real(SP),dimension(:,:),allocatable :: oposlist @@ -77,6 +77,7 @@ program CTEM allocate(nflux(3,domain%pnum)) allocate(crtscl(2,domain%pnum)) allocate(vdist(3,domain%vnum)) +allocate(rclist(6,domain%rcnum)) allocate(surf(user%gridsize,user%gridsize)) allocate(production_list(domain%pnum)) @@ -86,6 +87,10 @@ program CTEM ! Read in impactor velocity distribution call io_read_vdist(vdist,user,domain) +! Read in real crater list for quasi-MC run +call io_read_craterlist(rclist,user,domain) +write(*,*) rclist + write(*,*) "Initializing simulation domain and determining minimum impactor size" call init_domain(user,crater,domain,prod,pdist,vdist,crtscl,nflux) @@ -115,7 +120,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,curyear) + fracdone,nflux,ntotcrat,curyear,rclist) ! Get the last seed and save it to file call random_seed(get=seedarr)