attempt to add meltdist array; doesn't work yet

MintonGroup · Dec 12, 2022 · 2aad6a7 · 2aad6a7
1 parent c903522
commit 2aad6a7
Show file tree

Hide file tree

Showing 13 changed files with 115 additions and 19 deletions.
diff --git a/examples/global-lunar-bombardment/ctem.in b/examples/global-lunar-bombardment/ctem.in
@@ -11,12 +11,12 @@ testyoffset  0.0                            ! y-axis offset of crater center fro
 tallyonly    F                                 ! Tally the craters without generating any craters
 testtally    F   
 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
+realcraterlist    qmctest.in                ! list of 'real' craters for Quasi-MC runs
 
 
 
 ! IDL driver in uts
-interval     628
+interval     0.1
 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/examples/global-lunar-bombardment/temp.in b/examples/global-lunar-bombardment/temp.in
@@ -0,0 +1,76 @@
+! CTEM Input file
+
+
+! Testing input. These are used to perform non-Monte Carlo tests.
+testflag T!     F                                 ! Set to T to create a single crater with user-defined impactor properties
+testimp 10 !      15000.0       ! Diameter of test impactor (m)
+testvel      18.3e3                            ! Velocity of test crater (m/s)
+testang      45.0                              ! Impact angle of test crater (deg) - Default 90.0
+testxoffset  0.0                             ! x-axis offset of crater center from grid center (m) - Default 0.0
+testyoffset  0.0                            ! 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!      T                                 ! MC run constrained by non-MC 'real' craters given in a list
+realcraterlist    qmctest.in                ! list of 'real' craters for Quasi-MC runs
+
+
+
+! IDL driver in uts
+interval 1 !     0.1
+numintervals 1 ! 1                                  ! Total number of interval 1 !s (total time = interval 1 ! * numintervals 1 !) <--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 1 ! 
+saveshaded   F                                  ! Output shaded relief images  
+saverego     T                                  ! Output regolith map images 
+savepres     F                                  ! Output simplified console display images (presentation-compatible images)
+savetruelist T                                  ! Save the true cumulative crater distribution for each interval 1 ! (large file size)
+sfdcompare   LOLASethCraterCatalogv8gt20-binned.dat         ! File name for the SFD comparison file used in the console display
+shadedminh   -85.0                              ! Minimum height for shaded relief map (m) (Default - automatic)
+shadedmaxh   85.0                               ! Maximum height for shaded relief map (m) (Default - automatic)
+runtype      single                             ! Run type: options are normal / statistical
+                                                !    single:      craters accumulate in successive interval 1 !s
+                                                !    statistical: surface is reset between interval 1 !s
+
+! CTEM required inputs
+seed         76535                    ! Random number generator seed
+gridsize     2000                          ! Size of grid in pixels
+numlayers    10                                 ! Number of perched layers
+pix          3.08e3                           ! Pixel size (m)
+mat          rock                               ! Material (rock or ice)
+! Bedrock scaling parameters
+mu_b         0.55e0                             ! Experimentally derived parameter for bedrock crater scaling law 
+kv_b         0.20e0                             ! Experimentally derived parameter for bedrock crater scaling law 
+trho_b       2250.0e0                           ! Target density (bedrock) (kg/m**3)
+ybar_b       0.0e6                               ! Bedrock strength (Pa)
+! Regolith scaling parameters
+mu_r         0.55e0                             ! Experimentally derived parameter for regolith crater scaling law 
+kv_r         0.20e0                             ! Experimentally derived parameter for regolith crater scaling law
+trho_r       2250.0e0                           ! Target density (regolith) (kg/m**3)
+ybar_r       0.00e6                              ! Regolith strength (Pa)
+! Body parameters
+gaccel       1.62e0                             ! Gravitational acceleration at target (m/s**2)
+trad         1737.35e3                          ! Target radius (m)
+prho         2500.0e0                           ! Projectile density (kg/m**3)
+sfdfile      production.dat                     ! Impactor SFD file (col 1: Dimp (m), col 2: ! impactors > D
+velfile      lunar-MBA-impactor-velocities.dat  ! Impactor velocity dist file
+
+! Seismic shaking input (required if seismic shaking is set to T, otherwise ignored)
+doseismic    F                                  ! Perform seismic shaking calculations with each impact - Default F
+
+! Optional inputF These have internally set default values that work reasonable well. Comment them out with
+deplimit        9e99                            ! Depth limit for craters (m) - Default is to ignore.
+maxcrat         3.15e-2                          ! Fraction of gridsize that maximum crater can be - Default 1.0
+killatmaxcrater F                               ! Stop the run if a crater larger than the maximum is produced - Default F      
+basinimp        35.0e3                          ! Size of impactor to switch to lunar basin scaling law - Default is to ignore
+docollapse      T                               ! Do slope collapse - Default T
+dosoftening     T                               ! Do ejecta softening - Default T
+doangle         T                               ! Vary the impact angle. Set to F to have only vertical impacts - Default T
+Kd1             0.0001
+psi             2.000
+fe              4.00
+ejecta_truncation 4.0
+doregotrack T
+dorays  F 
+superdomain F
+dorealistic F
diff --git a/src/Makefile.am b/src/Makefile.am
@@ -11,7 +11,7 @@ AM_FCFLAGS =  $(IDEBUG)
 #ifort debug flags
 
 #gfortran optimized flags
-#AM_FCFLAGS =  -O3 -fopenmp -ffree-form  -g -fbounds-check -fbacktrace
+AM_FCFLAGS =  -O3 -fopenmp -ffree-form  -g -fbounds-check -fbacktrace
 #gfortran debug flags
 #AM_FCFLAGS =  -O0 -g -fopenmp -fbounds-check -Wall -Warray-bounds -Warray-temporaries -Wimplicit-interface -ffree-form -fsanitize-address-use-after-scope -fstack-check  -fsanitize=bounds-strict -fsanitize=undefined -fsanitize=signed-integer-overflow -fsanitize=object-size -fstack-protector-all
 

diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90
@@ -151,13 +151,15 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
    oldpbarpos = 0
    do while (icrater < ntotcrat)
       makecrater = .true.
+      domain%currentqmc = .false.
       timestamp_old = real(curyear + real(icrater,kind=DP) / real(ntotcrat,kind=DP) * user%interval,kind=SP)
       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
          if ((user%rctime > timestamp_old) .and. (user%rctime < crater%timestamp)) then
+            domain%currentqmc = .true.
             write(message,*) "Real @ ", crater%timestamp
             call io_updatePbar(message)
             user%testflag = .true.
@@ -189,6 +191,7 @@ 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.
+            domain%nqmc = domain%rccount
             domain%rccount = domain%rccount + 1
             if (domain%rccount > domain%rcnum) then
                write(message,*) "Real crater list complete."

diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90
@@ -155,8 +155,10 @@ 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
-   integer(I4B) :: rccount   ! Current quasimc crater
+   integer(I4B) :: rcnum     ! size of real crater list array for quasi-MC
+   integer(I4B) :: rccount   ! quasimc crater for iterating through the list of craters
+   integer(I4B) :: nqmc      ! current Quasi-MC crater for regolith distribution
+   logical      :: currentqmc ! is the current crater a Quasi-MC 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
@@ -302,6 +304,7 @@ module module_globals
 character(*),parameter :: DATFILE    = 'ctem.dat' 
 character(*),parameter :: MASSFILE   = 'impactmass.dat'
 character(*),parameter :: RCFILE     = 'craterlist.dat'
+character(*),parameter :: MDFILE     = 'surface_meltdist.dat'
 
 ! Global variables 
 integer(I4B),parameter :: PBCLIM = 1             ! periodic boundary condition limit

diff --git a/src/io/io_write_regotrack.f90 b/src/io/io_write_regotrack.f90
@@ -16,14 +16,15 @@
 !  Notes       :  
 !
 !**********************************************************************************************************************************
-subroutine io_write_regotrack(user,surf)
+subroutine io_write_regotrack(user,surf,domain)
    use module_globals
    use module_io, EXCEPT_THIS_ONE => io_write_regotrack
    implicit none
 
    ! Arguments
    type(usertype),intent(in) :: user
    type(surftype),dimension(:,:),intent(in) :: surf
+   type(domaintype),intent(in) :: domain
 
    ! Regotrack Internals
    integer(I4B) :: i,j,k
@@ -32,11 +33,12 @@ subroutine io_write_regotrack(user,surf)
    integer(I4B), parameter :: FREGO = 11
    integer(I4B), parameter :: FCOMP = 12
    integer(I4B), parameter :: FAGE = 13
+   integer(I4B), parameter :: FMD = 14
    !type(regolisttype),pointer :: current => null()
    type(regodatatype),dimension(:),allocatable :: current
    integer(I4B),dimension(user%gridsize,user%gridsize) :: stacks_num
    real(DP),dimension(:),allocatable :: meltfrac, thickness, comp
-   real(SP),dimension(:,:),allocatable :: age
+   real(SP),dimension(:,:),allocatable :: age, meltdist
    integer(kind=8) :: recsize
    real(DP) :: dtmp
    real(SP) :: stmp
@@ -49,6 +51,7 @@ subroutine io_write_regotrack(user,surf)
    open(FREGO,file=REGOFILE,status='replace',form='unformatted')
    open(FCOMP,file=COMPFILE,status='replace',form='unformatted')
    open(FAGE,file=AGEFILE,status='replace',form='unformatted')
+   open(FMD,file=MDFILE,status='replace',form='unformatted')
 
    ! First pass to get stack numbers
    stacks_num(:,:) = 0
@@ -71,7 +74,7 @@ subroutine io_write_regotrack(user,surf)
       do i=1,user%gridsize
          !current => surf(i,j)%regolayer
          N = stacks_num(i,j)
-         allocate(meltfrac(N),thickness(N),comp(N),age(MAXAGEBINS,N))
+         allocate(meltfrac(N),thickness(N),comp(N),age(MAXAGEBINS,N),meltdist(domain%rcnum,N))
          allocate(current,source=surf(i,j)%regolayer)
          do k=1,N
             ! meltfrac(k) = current%regodata%meltfrac
@@ -83,20 +86,25 @@ subroutine io_write_regotrack(user,surf)
             thickness(k) = current(k)%thickness
             comp(k) = current(k)%comp
             age(:,k) = current(k)%age(:)
+            write(*,*) size(current(k)%age)
+            write(*,*) size(current(k)%meltdist)
+            meltdist(:,k) = current(k)%meltdist(:)
 
          end do
          deallocate(current)
          write(FMELT) meltfrac(:)
          write(FREGO) thickness(:)
          write(FCOMP) comp(:)
          write(FAGE) age(:,:)
-         deallocate(meltfrac,thickness,comp,age)
+         write(FMD) meltdist(:,:)
+         deallocate(meltfrac,thickness,comp,age,meltdist)
       end do 
    end do
    close(FMELT)
    close(FREGO)
    close(FCOMP)
    close(FAGE)
+   close(FMD)
 
    recsize = sizeof(itmp) * user%gridsize * user%gridsize
    open(LUN,file=STACKNUMFILE,status='replace',form='unformatted',recl=recsize,access='direct')

diff --git a/src/io/io_write_surf.f90 b/src/io/io_write_surf.f90
@@ -16,14 +16,15 @@
 !  Notes       :  
 !
 !**********************************************************************************************************************************
-subroutine io_write_surf(user,surf)
+subroutine io_write_surf(user,surf,domain)
    use module_globals
    use module_io, EXCEPT_THIS_ONE => io_write_surf
    implicit none
 
    ! Arguments
    type(usertype),intent(in) :: user
    type(surftype),dimension(:,:),intent(in) :: surf
+   type(domaintype),intent(in) :: domain
 
    ! Internals
    integer(I4B) :: i
@@ -67,7 +68,7 @@ subroutine io_write_surf(user,surf)
    end do
    close(LUN)
 
-   if (user%doregotrack) call io_write_regotrack(user,surf) 
+   if (user%doregotrack) call io_write_regotrack(user,surf,domain) 
 
 !   if (user%docrustal_thinning) then
 !      recsize = sizeof(itmp) * user%gridsize * user%gridsize

diff --git a/src/io/module_io.f90 b/src/io/module_io.f90
@@ -69,20 +69,22 @@ end subroutine io_read_regotrack
   end interface  
 
   interface   
-      subroutine io_write_surf(user,surf)
+      subroutine io_write_surf(user,surf,domain)
       use module_globals
       implicit none
       type(usertype),intent(in) :: user
       type(surftype),dimension(:,:),intent(in) :: surf
+      type(domaintype),intent(in) :: domain
       end subroutine io_write_surf
    end interface
 
   interface 
-      subroutine io_write_regotrack(user,surf)
+      subroutine io_write_regotrack(user,surf,domain)
       use module_globals
       implicit none
       type(usertype),intent(in) :: user
-      type(surftype),dimension(:,:),intent(in) :: surf      
+      type(surftype),dimension(:,:),intent(in) :: surf   
+      type(domaintype),intent(in) :: domain   
       end subroutine io_write_regotrack
   end interface 
 

diff --git a/src/main/CTEM.f90 b/src/main/CTEM.f90
@@ -147,7 +147,7 @@ program CTEM
 call io_write_tally(truedist,truelist(:,1:ntrue),obsdist,obslist,oposlist,depthdiam,degradation_state)
 if (.not.user%tallyonly) then
    write(*,*) "Writing surface files"
-   call io_write_surf(user,surf)
+   call io_write_surf(user,surf,domain)
 end if
 
 if (user%testflag) then ! Draw a profile across the crater

diff --git a/src/regolith/module_regolith.f90 b/src/regolith/module_regolith.f90
@@ -246,10 +246,11 @@ end subroutine regolith_depth_model
    end interface
 
    interface
-      subroutine regolith_melt_glass(user,crater,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints,melt)
+      subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints,melt)
       use module_globals
       type(usertype),intent(in)        :: user
       type(cratertype),intent(in)      :: crater
+      type(domaintype),intent(in)      :: domain
       real(DP),intent(in)              :: age
       real(DP),intent(in)              :: age_resolution
       real(DP),intent(in)              :: ebh

diff --git a/src/regolith/regolith_melt_glass.f90 b/src/regolith/regolith_melt_glass.f90
@@ -55,7 +55,7 @@
 !  Notes       :  
 !
 !**********************************************************************************************************************************
-subroutine regolith_melt_glass(user,crater,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints,melt)
+subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints,melt)
    use module_globals 
    use module_util
    use module_regolith, EXCEPT_THIS_ONE => regolith_melt_glass
@@ -64,6 +64,7 @@ subroutine regolith_melt_glass(user,crater,age,age_resolution,ebh,rm,eradc,lrad,
    ! Arguments
    type(usertype),intent(in)        :: user
    type(cratertype),intent(in)      :: crater
+   type(domaintype),intent(in)      :: domain
    real(DP),intent(in)              :: age
    real(DP),intent(in)              :: age_resolution
    real(DP),intent(in)              :: ebh
@@ -164,6 +165,7 @@ subroutine regolith_melt_glass(user,crater,age,age_resolution,ebh,rm,eradc,lrad,
            volm1        = regolith_streamtube_volume_func(eradi,0.0_DP,xmints,deltar)
            melt         = volm1 - volv1
            newlayer%meltfrac = melt/(vst-volv1)
+           newlayer%meltdist(domain%nqmc) = newlayer%meltfrac
    end if
 
    n_age = max(ceiling(age / age_resolution), 1)

diff --git a/src/regolith/regolith_streamtube.f90 b/src/regolith/regolith_streamtube.f90
@@ -226,7 +226,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,
    ! Purpose 2: Once we have the size information of a stream tube, we can
    ! calculate the distal melt: the precursor of glass spherules within a
    ! stream tube. The result is contained in a linked list "newlayer".
-   call regolith_melt_glass(user,crater,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints, volm)
+   call regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints, volm)
    ! if (eradc>rm) then
    !    write(*,*) 'eradc > rm!'
    !    write(*,*) ebh, exp(ejb(k)%thick)

diff --git a/src/regolith/regolith_superdomain.f90 b/src/regolith/regolith_superdomain.f90
@@ -60,7 +60,7 @@ subroutine regolith_superdomain(user,crater,domain,regolayer,ejdistribution,xpi,
    erad        = cvpg**(user%mu_b) * (lrad**(-0.5_DP * user%mu_b)) * (crater%rad)**(user%mu_b * 0.5_DP + 1.0_DP)
    vej         = cvpg * sqrt(user%gaccel * crater%grad) * (erad / crater%grad)**(-1.0_DP / user%mu_b) !equation 18 in Richardson 2009
    lrad        = ( vej **2 ) / user%gaccel !assume ejection angle is 45 degree.
-   call regolith_melt_glass(user,crater,age,age_resolution,ebh,rm,erad,lrad,deltar,newlayer,xmints,melt) 
+   call regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad,lrad,deltar,newlayer,xmints,melt) 
    call util_push_array(regolayer,newlayer)
 
    return