diff --git a/examples/global-lunar-bombardment/ctem.in b/examples/global-lunar-bombardment/ctem.in index 6e632d39..0a1b45e1 100755 --- a/examples/global-lunar-bombardment/ctem.in +++ b/examples/global-lunar-bombardment/ctem.in @@ -10,14 +10,14 @@ testxoffset 0.0e0 ! x-axis offset of crater center testyoffset 0.0e0 ! y-axis offset of crater center from grid center (m) - Default 0.0 tallyonly F ! Tally the craters without generating any craters testtally F -quasimc F ! MC run constrained by non-MC 'real' craters given in a list +quasimc T ! MC run constrained by non-MC 'real' craters given in a list realcraterlist craterlist.in ! list of 'real' craters for Quasi-MC runs ! IDL driver in uts -interval 3.0 -numintervals 1 ! Total number of intervals +interval 5.86055 +numintervals 100 ! 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)) popupconsole F ! Pop up console window every output interval diff --git a/examples/global-lunar-bombardment/ctem_driver.py b/examples/global-lunar-bombardment/ctem_driver.py index f116cf48..edcf06ef 100644 --- a/examples/global-lunar-bombardment/ctem_driver.py +++ b/examples/global-lunar-bombardment/ctem_driver.py @@ -80,7 +80,7 @@ print("quasi-MC mode is ON") craterlistfile = parameters['workingdir'] + parameters['realcraterlist'] rclist = ctem_io_readers.read_formatted_ascii(craterlistfile, skip_lines = 0) - rclist[:,5] = parameters['interval'] - craterproduction.Tscale(rclist[:,5], 'NPF_Moon') + rclist[:,5] = (parameters['interval'] * parameters['numintervals']) - craterproduction.Tscale(rclist[:,5], 'NPF_Moon') ctem_io_writers.write_realcraters(parameters, rclist) #Create impactor production population diff --git a/src/Makefile.am b/src/Makefile.am index 13b43276..d7ae0072 100755 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -38,6 +38,7 @@ util/util_init_list.f90\ io/io_read_const.f90\ io/io_get_token.f90\ io/io_input.f90\ +io/io_read_craterlist.f90\ io/io_read_prod.f90\ io/io_read_regotrack.f90\ io/io_read_porotrack.f90\ diff --git a/src/Makefile.in b/src/Makefile.in index 6e212687..9a33fe0b 100644 --- a/src/Makefile.in +++ b/src/Makefile.in @@ -115,12 +115,13 @@ am_CTEM_OBJECTS = globals/module_globals.$(OBJEXT) \ util/util_push.$(OBJEXT) util/util_traverse_pop.$(OBJEXT) \ util/util_destroy_list.$(OBJEXT) util/util_init_list.$(OBJEXT) \ io/io_read_const.$(OBJEXT) io/io_get_token.$(OBJEXT) \ - io/io_input.$(OBJEXT) io/io_read_prod.$(OBJEXT) \ - io/io_read_regotrack.$(OBJEXT) io/io_read_porotrack.$(OBJEXT) \ - io/io_read_vdist.$(OBJEXT) io/io_read_surf.$(OBJEXT) \ - io/io_ejecta_table.$(OBJEXT) io/io_write_dist.$(OBJEXT) \ - io/io_write_tally.$(OBJEXT) io/io_write_surf.$(OBJEXT) \ - io/io_write_const.$(OBJEXT) io/io_write_regotrack.$(OBJEXT) \ + io/io_input.$(OBJEXT) io/io_read_craterlist.$(OBJEXT) \ + io/io_read_prod.$(OBJEXT) io/io_read_regotrack.$(OBJEXT) \ + io/io_read_porotrack.$(OBJEXT) io/io_read_vdist.$(OBJEXT) \ + io/io_read_surf.$(OBJEXT) io/io_ejecta_table.$(OBJEXT) \ + io/io_write_dist.$(OBJEXT) io/io_write_tally.$(OBJEXT) \ + io/io_write_surf.$(OBJEXT) io/io_write_const.$(OBJEXT) \ + io/io_write_regotrack.$(OBJEXT) \ io/io_write_porotrack.$(OBJEXT) io/io_crater_profile.$(OBJEXT) \ io/io_updatePbar.$(OBJEXT) io/io_resetPbar.$(OBJEXT) \ io/io_splash.$(OBJEXT) ejecta/ejecta_emplace.$(OBJEXT) \ @@ -339,6 +340,7 @@ util/util_init_list.f90\ io/io_read_const.f90\ io/io_get_token.f90\ io/io_input.f90\ +io/io_read_craterlist.f90\ io/io_read_prod.f90\ io/io_read_regotrack.f90\ io/io_read_porotrack.f90\ @@ -594,6 +596,8 @@ io/io_read_const.$(OBJEXT): io/$(am__dirstamp) \ io/io_get_token.$(OBJEXT): io/$(am__dirstamp) \ io/$(DEPDIR)/$(am__dirstamp) io/io_input.$(OBJEXT): io/$(am__dirstamp) io/$(DEPDIR)/$(am__dirstamp) +io/io_read_craterlist.$(OBJEXT): io/$(am__dirstamp) \ + io/$(DEPDIR)/$(am__dirstamp) io/io_read_prod.$(OBJEXT): io/$(am__dirstamp) \ io/$(DEPDIR)/$(am__dirstamp) io/io_read_regotrack.$(OBJEXT): io/$(am__dirstamp) \ diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index 94b3c4eb..7d519ac7 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -134,6 +134,10 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt icrater = icrater + 1 crater%timestamp = real(curyear + real(icrater,kind=DP) / real(ntotcrat,kind=DP) * user%interval,kind=SP) pbarpos = nint(real(icrater) / real(ntotcrat) * PBARRES) + !if in quasiMC mode: check to see if it's time for a real crater + !if (user%doquasimc) then + + !end if ! generate random crater call crater_generate(user,crater,domain,prod,production_list,vdist,surf) if (user%testflag) write(*,*) 'Dcrat = ',crater%fcrat diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index f2c0c31d..e59a245f 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -137,6 +137,7 @@ module module_globals real(DP) :: ejbres ! Ejecta blanket lookup table resolution integer(I4B) :: pnum ! size of production function array integer(I4B) :: vnum ! size of velocity distribution array + integer(I4B) :: rcnum ! size of real crater list array for quasi-MC real(DP) :: vescsq ! Escape velocity at target integer(I4B) :: vlo ! Index of lowest valid velocity in the velocity distribution file integer(I4B) :: vhi ! Index of highest valid velocity in the velocity distribution file @@ -171,6 +172,7 @@ module module_globals logical :: docollapse ! Set T to use the slope collapse model (turning off speeds up the code for testing) logical :: doangle ! Set to F to only do vertical impacts, otherwise do range of angles (default is T) logical :: doporosity ! Porosity on/off flg. Set to F to turn the model off. Default F. + logical :: doquasimc ! set to T for quasi-MC run. Default F. real(DP) :: basinimp ! Impactor size to switch to multiring basin real(DP) :: maxcrat ! fraction that maximum crater can be relative to grid real(DP) :: deplimit ! complex crater depth limit @@ -212,6 +214,7 @@ module module_globals real(DP) :: testyoffset ! Offset of test crater from center in y direction (m) logical :: tallyonly ! Only run the tally routine (don't generate any new craters) logical :: testtally ! Set to T to count all non-cookie cut craters, regardless of score + real(DP) :: rctime ! time (in interval units) for emplacement of quasi-MC crater. Default 1E30 (aka never emplace real crater) ! IDL driver variables character(STRMAX) :: impfile ! Name of impactor size distribution file (impacts per m^2 per y) @@ -228,6 +231,7 @@ module module_globals real(DP) :: shadedminh ! Minimum height for shaded relief map (m) real(DP) :: shadedmaxh ! Maximum height for shaded relief map (m) character(STRMAX) :: sfdcompare ! Type of run: 0 for normal, 1 for statistical (domain is reset between intervals) + character(STRMAX) :: realcraterlist ! This is only included here so the terminal doesn't return "Unknown parameter" end type usertype ! Derived data type for the ejecta blanket table elements @@ -272,6 +276,7 @@ module module_globals character(*),parameter :: CRTSCLFILE = 'craterscale.dat' character(*),parameter :: DATFILE = 'ctem.dat' character(*),parameter :: MASSFILE = 'impactmass.dat' +character(*),parameter :: RCFILE = 'craterlist.dat' !not sure if this is where this line should go, but putting it here for now... ! Global variables integer(I4B),parameter :: PBCLIM = 1 ! periodic boundary condition limit diff --git a/src/init/init_dist.f90 b/src/init/init_dist.f90 index cbc0af67..9863587c 100644 --- a/src/init/init_dist.f90 +++ b/src/init/init_dist.f90 @@ -69,6 +69,27 @@ subroutine init_dist(user,domain) end if close(LUN) + ! Get size of real crater list array + if (user%doquasimc) then + open(unit=LUN,file=rcfile,status="old",iostat=ierr) + if (ierr /= 0) then + write(*,*) "Unable to open file ", trim(rcfile) + stop + end if + + domain%rcnum = 0 + do + read(LUN,*,iostat=ierr) testreal,testreal,testreal,testreal,testreal,testreal + if (ierr/=0) exit + domain%rcnum = domain%rcnum + 1 + end do + + if (domain%rcnum == 0) then + write(*,*) "No valid entries in ",trim(rcfile) + end if + close(LUN) + end if + return end subroutine init_dist diff --git a/src/io/io_input.f90 b/src/io/io_input.f90 index 258a83ae..83f6c39f 100644 --- a/src/io/io_input.f90 +++ b/src/io/io_input.f90 @@ -90,6 +90,7 @@ subroutine io_input(infile,user) user%psi = 2.0_DP user%fe = 1.0_DP user%ejecta_truncation = 10.0_DP + user%doquasimc = .false. open(unit=LUN,file=infile,status="old",iostat=ierr) if (ierr /= 0) then @@ -251,6 +252,11 @@ subroutine io_input(infile,user) call io_get_token(line, ilength, ifirst, ilast, ierr) token = line(ifirst:ilast) read(token, *) user%docollapse + case ("QUASIMC") + ifirst = ilast + 1 + call io_get_token(line, ilength, ifirst, ilast, ierr) + token = line(ifirst:ilast) + read(token, *) user%doquasimc case ("TESTXOFFSET") ifirst = ilast + 1 call io_get_token(line, ilength, ifirst, ilast, ierr) @@ -426,6 +432,11 @@ subroutine io_input(infile,user) call io_get_token(line, ilength, ifirst, ilast, ierr) token = line(ifirst:ilast) read(token, *) user%sfdcompare + case ("REALCRATERLIST") + ifirst = ilast + 1 + call io_get_token(line, ilength, ifirst, ilast, ierr) + token = line(ifirst:ilast) + read(token, *) user%realcraterlist case ("SHADEDMINH") ifirst = ilast + 1 call io_get_token(line, ilength, ifirst, ilast, ierr) diff --git a/src/io/io_read_craterlist.f90 b/src/io/io_read_craterlist.f90 new file mode 100644 index 00000000..138fc750 --- /dev/null +++ b/src/io/io_read_craterlist.f90 @@ -0,0 +1,52 @@ +!********************************************************************************************************************************** +! +! Unit Name : io_read_craterlist +! Unit Type : subroutine +! Project : CTEM +! Language : Fortran 2003 +! +! Description : Reads in list of 'real' craters for quasi-MC runs +! +! Input +! Arguments : +! +! Output +! Arguments : rclist : List of 'real' craters to read in +! +! Notes : +! +!********************************************************************************************************************************** +subroutine io_read_craterlist(user, domain) + use module_globals + use module_io, EXCEPT_THIS_ONE => io_read_craterlist + implicit none + + ! Arguments + type(usertype),intent(inout) :: user + type(domaintype),intent(in) :: domain + + ! Internals + integer(I4B) :: i,ierr + integer(I4B),parameter :: LUN=7 + + ! Executable code + + ! Read the next crater from the craterlist + + open(unit=LUN,file=rcfile,status='old',iostat=ierr) + if (ierr /= 0) then + write(*,*) "Unable to open file ",trim(rcfile) + stop + end if + + + do i=1,domain%rcnum + 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 +end subroutine io_read_craterlist + + diff --git a/src/io/module_io.f90 b/src/io/module_io.f90 index 0dfcf6d8..d3d36f60 100644 --- a/src/io/module_io.f90 +++ b/src/io/module_io.f90 @@ -121,6 +121,16 @@ subroutine io_read_prod(prod,user,domain) end subroutine io_read_prod end interface + interface + 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 + end interface + interface subroutine io_read_vdist(vdist,user,domain) use module_globals