quasiMC routine seems to work until the end when Python errors are thrown

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
+90856.5  15.0e3   90.0  -5.17529e5  1.037286e6  0.4
+18095.0  15.0e3   90.0  -9.38629e5  1.3466084e6 0.2

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))

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)

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
 

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
 
 

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

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)