Fortran reads list of 'real' craters but doesn't do anything with it yet

MintonGroup · Mar 31, 2021 · 54d4a22 · 54d4a22
1 parent 3f729e9
commit 54d4a22
Show file tree

Hide file tree

Showing 10 changed files with 118 additions and 10 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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