Took the first step towards melt tracking with the addition of an array 'meltdist'

examples/global-lunar-bombardment/qmctest.in

@@ -0,0 +1,6 @@
+#Dcrat(m)   vel(m/s)   ang(deg)  xoffset(m)  yoffset(m)  t(Ga)
+90000 18300.0 45.0 0 0 0.0009
+50000 18300.0 45.0 50000 50000 0.0007
+40000 18300.0 45.0 -80000 -80000 0.0005
+70000 18300.0 45.0 -80000 50000 0.0003
+10000 18300.0 45.0 50000 -80000 0.0001

examples/global-lunar-bombardment/read_regotrack.py

@@ -3,10 +3,8 @@
 import ctem
 import os
 
-#path = '/Users/owner/Desktop/yh/ctem_20220808'
-#path = '/Users/owner/Desktop/raytest/off'
-#path = '/Users/owner/Desktop/melt_distribution/'
-path = os.getcwd()
+path = '/Users/owner/Documents/git/CTEM/examples/global-lunar-bombardment'
+#path = os.getcwd()
 gridsize = 2000
 pix = 2000
 pixkm = pix / 1000
@@ -60,7 +58,7 @@
 
 fig2 = plt.figure(figsize=(gridsize/dpi, gridsize/dpi), dpi=dpi)
 ax2 = plt.axes([0,0,1,1])
-md = ax2.imshow(meltdepth, interpolation='nearest')#, vmin=0.0, vmax=0.5)
+md = ax2.imshow(meltdepth, interpolation='nearest')#, vmin=0.0, vmax=0.1)
 cbar = fig2.colorbar(md, label = 'Melt Fraction')#, shrink=0.8)
 ax2.invert_yaxis()
 

src/crater/crater_populate.f90

@@ -109,8 +109,8 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
 
    ! read initial quasi-MC position
    if (user%doquasimc) then
-      rccount = 1
-      user%rctime = rclist(6,rccount)
+      domain%rccount = 1
+      user%rctime = rclist(6,domain%rccount)
    end if
 
 
@@ -161,11 +161,11 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
             write(message,*) "Real @ ", crater%timestamp
             call io_updatePbar(message)
             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) 
+            user%testimp = rclist(1, domain%rccount)
+            user%testvel = rclist(2, domain%rccount)
+            user%testang = rclist(3, domain%rccount)
+            user%testxoffset = rclist(4, domain%rccount)
+            user%testyoffset = rclist(5, domain%rccount) 
          end if
       end if
       ! generate random crater
@@ -189,13 +189,13 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
       if (user%doquasimc) then
          if (crater%timestamp > user%rctime) then
             user%testflag = .false.
-            rccount = rccount + 1
-            if (rccount > domain%rcnum) then
+            domain%rccount = domain%rccount + 1
+            if (domain%rccount > domain%rcnum) then
                write(message,*) "Real crater list complete."
                call io_updatePbar(message)
                user%rctime = 1e30
             else
-               user%rctime = rclist(6,rccount)
+               user%rctime = rclist(6,domain%rccount)
             end if
          end if
       end if

src/globals/module_globals.f90

@@ -72,6 +72,7 @@ module module_globals
    real(DP) :: porosity   ! Porosity: Maximum 1, Minium 0.
    real(DP) :: damage     ! Damage  : Maximum 1, Minium 0.
    real(DP) :: depth      ! Absolute location with respect to the initial surface. 
+   real(SP),dimension(:),allocatable :: meltdist !its dimension should be the number of quasimc craters
 end type regodatatype
 
 type regolisttype
@@ -155,6 +156,7 @@ module module_globals
    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
+   integer(I4B) :: rccount   ! Current quasimc crater
    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
@@ -297,9 +299,9 @@ module module_globals
 character(*),parameter :: CRTSCLFILE = 'craterscale.dat'
 character(*),parameter :: SFDFILE    = 'production.dat'
 character(*),parameter :: USERFILE   = 'ctem.in' 
-character(*),parameter :: DATFILE    = 'ctem.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...
+character(*),parameter :: RCFILE     = 'craterlist.dat'
 
 ! Global variables 
 integer(I4B),parameter :: PBCLIM = 1             ! periodic boundary condition limit
@@ -353,7 +355,4 @@ module module_globals
 integer(I4B),parameter :: rayq = 4
 real(DP),parameter :: rayfmult  = (5)**(-4.0_DP / (1.2_DP))  
 
-! quasi-MC test variables
-integer(I4B) :: rccount !start line to be read for quasi-MC (should be 0)
-
 end module module_globals

src/init/init_dist.f90

@@ -88,6 +88,8 @@ subroutine init_dist(user,domain)
          write(*,*) "No valid entries in ",trim(rcfile)
       end if
       close(LUN)
+   else
+      domain%rcnum = 1
    end if
 
 

src/init/init_regolith_stack.f90

@@ -18,7 +18,7 @@
 !  Notes       :  
 !
 !**********************************************************************************************************************************
-subroutine init_regolith_stack(user,surf)
+subroutine init_regolith_stack(user,surf,domain)
    use module_globals
    use module_util
    use module_init, EXCEPT_THIS_ONE => init_regolith_stack
@@ -27,6 +27,7 @@ subroutine init_regolith_stack(user,surf)
    ! Arguments
    type(usertype),intent(in) :: user
    type(surftype),dimension(:,:),intent(inout) :: surf
+   type(domaintype),intent(in)    :: domain
    type(regodatatype) :: bedrock 
    integer(I4B) :: xp,yp,k
 
@@ -49,7 +50,7 @@ subroutine init_regolith_stack(user,surf)
       do xp = 1, user%gridsize
 
          !call util_init_list(surf(xp,yp)%regolayer,initstat)
-         call util_init_array(surf(xp,yp)%regolayer,initstat)
+         call util_init_array(surf(xp,yp)%regolayer,domain,initstat)
 
          if (initstat) then
              call util_push_array(surf(xp,yp)%regolayer,bedrock)

src/init/init_surf.f90

@@ -16,14 +16,15 @@
 !  Notes       :  
 !
 !**********************************************************************************************************************************
-subroutine init_surf(user,surf)
+subroutine init_surf(user,surf,domain)
    use module_globals
    use module_init, EXCEPT_THIS_ONE => init_surf
    implicit none
 
    ! Arguments
    type(usertype),intent(in) :: user
    type(surftype),dimension(:,:),intent(out) :: surf
+   type(domaintype),intent(in)    :: domain
 
    ! Internal variables
    integer(I4B) :: layer
@@ -37,7 +38,7 @@ subroutine init_surf(user,surf)
    end do
    !if (user%docrustal_thinning) surf%mantle = 0._DP
 
-   if (user%doregotrack) call init_regolith_stack(user,surf)
+   if (user%doregotrack) call init_regolith_stack(user,surf,domain)
 
    return
    end subroutine init_surf

src/init/module_init.f90

@@ -16,11 +16,12 @@ module module_init
 save
 
    interface
-      subroutine init_surf(user,surf)
+      subroutine init_surf(user,surf,domain)
       use module_globals
       implicit none
       type(usertype),intent(in) :: user
       type(surftype),dimension(:,:),intent(out) :: surf
+      type(domaintype),intent(in)    :: domain
       end subroutine init_surf
    end interface
 
@@ -48,11 +49,12 @@ end subroutine init_domain
    end interface 
 
    interface
-      subroutine init_regolith_stack(user,surf)
+      subroutine init_regolith_stack(user,surf,domain)
       use module_globals
       implicit none
       type(usertype),intent(in) :: user
       type(surftype),dimension(:,:),intent(inout) :: surf
+      type(domaintype),intent(in)    :: domain
       end subroutine init_regolith_stack
    end interface
 

src/io/io_read_regotrack.f90

@@ -16,7 +16,7 @@
 !  Notes       :  
 !
 !**********************************************************************************************************************************
-subroutine io_read_regotrack(user,surf)
+subroutine io_read_regotrack(user,surf,domain)
    use module_globals
    use module_util
    use module_io, EXCEPT_THIS_ONE => io_read_regotrack
@@ -25,6 +25,7 @@ subroutine io_read_regotrack(user,surf)
    ! Arguments
    type(usertype),intent(in) :: user
    type(surftype),dimension(:,:),intent(out) :: surf
+   type(domaintype),intent(in)    :: domain
 
    ! Internals
    integer(I4B), parameter :: LUN=7
@@ -85,7 +86,7 @@ subroutine io_read_regotrack(user,surf)
       do i=1,user%gridsize
 
          !call util_init_list(surf(i,j)%regolayer,initstat)
-         call util_init_array(surf(i,j)%regolayer,initstat)
+         call util_init_array(surf(i,j)%regolayer,domain,initstat)
 
          allocate(regotopi(stacks_num(i,j)))
          allocate(compi(stacks_num(i,j)))

src/io/io_read_surf.f90

@@ -16,14 +16,15 @@
 !  Notes       :  
 !
 !**********************************************************************************************************************************
-subroutine io_read_surf(user,surf)
+subroutine io_read_surf(user,surf,domain)
    use module_globals
    use module_io, EXCEPT_THIS_ONE => io_read_surf
    implicit none
 
    ! Arguments
    type(usertype),intent(in) :: user
    type(surftype),dimension(:,:),intent(out) :: surf
+   type(domaintype),intent(in)    :: domain
 
    ! Internals
    integer(I4B) :: ioerr,i
@@ -89,7 +90,7 @@ subroutine io_read_surf(user,surf)
    end do
    close(LUN)
 
-   if (user%doregotrack) call io_read_regotrack(user,surf)
+   if (user%doregotrack) call io_read_regotrack(user,surf,domain)
 
    !if (user%docrustal_thinning) then
    !   recsize=sizeof(itmp)*user%gridsize*user%gridsize

src/io/module_io.f90

@@ -49,20 +49,22 @@ end subroutine io_write_const
    end interface
 
    interface   
-      subroutine io_read_surf(user,surf)
+      subroutine io_read_surf(user,surf,domain)
       use module_globals
       implicit none
       type(usertype),intent(in) :: user
       type(surftype),dimension(:,:),intent(out) :: surf
+      type(domaintype),intent(in)    :: domain
       end subroutine io_read_surf
    end interface
 
   interface
-      subroutine io_read_regotrack(user,surf)
+      subroutine io_read_regotrack(user,surf,domain)
       use module_globals
       implicit none
       type(usertype),intent(in) :: user
       type(surftype),dimension(:,:),intent(out) :: surf
+      type(domaintype),intent(in)    :: domain
       end subroutine io_read_regotrack
   end interface  
 

src/main/CTEM.f90

@@ -113,9 +113,9 @@ program CTEM
 
 ! Read in old grid arrays, production function, and velocity distributions
 if (restart .or. user%tallyonly) then
-   call io_read_surf(user,surf)
+   call io_read_surf(user,surf,domain)
 else
-   call init_surf(user,surf)
+   call init_surf(user,surf,domain)
 end if
 
 if (.not.user%tallyonly) then

src/util/module_util.f90

@@ -110,10 +110,11 @@ end subroutine util_destroy_list
 ! end interface
 
 interface
-   subroutine util_init_array(regolayer,initstat)
+   subroutine util_init_array(regolayer,domain,initstat)
    use module_globals
    implicit none
    type(regodatatype),dimension(:),allocatable,intent(inout) :: regolayer
+   type(domaintype),intent(in)    :: domain
    logical, intent(out)     :: initstat
    end subroutine util_init_array
 end interface

src/util/util_init_array.f90

@@ -19,13 +19,14 @@
 !  Notes       : 
 !
 !**********************************************************************************************************************************
-subroutine util_init_array(regolayer,initstat)
+subroutine util_init_array(regolayer,domain,initstat)
     use module_globals
     use module_util, EXCEPT_THIS_ONE => util_init_array
     implicit none
 
     ! Arguments
     type(regodatatype),dimension(:),allocatable,intent(inout) :: regolayer
+    type(domaintype),intent(in)    :: domain
     logical, intent(out)       :: initstat
 
     ! Internal variables
@@ -52,6 +53,8 @@ subroutine util_init_array(regolayer,initstat)
         regolayer(1)%meltfrac = 0.0_DP
         regolayer(1)%porosity = 0.0_DP
         regolayer(1)%age(:)   = 0.0_SP
+        allocate(regolayer(1)%meltdist(domain%rcnum))
+        regolayer(1)%meltdist(:) = 0.0_SP
     end if
     !    else
     !       write(*,*) 'util_init_list: Initialization failed. Exhausted memory.'