Successfully changed age interval to Ga from interval time

src/crater/crater_populate.f90

@@ -89,11 +89,11 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
    integer(I4B)                      :: craters_since_subpixel_mix, icrater_last_subpixel_mix
 
    ! doregotrack & age simulation test
-   real(DP)              :: melt, clock, age, thick
+   real(DP)              :: melt, clock, age, thick, maxage
    real(SP),dimension(user%gridsize, user%gridsize)  :: agetop
    real(SP),dimension(60)                            :: agetot
    type(regolisttype),pointer                        :: current => null()
-   real(DP)              :: age_resolution, ageGa, oldGa
+   real(DP)              :: age_resolution, maxageGa, oldGa, agemin 
    integer(I4B)          :: age_counter
 
    nmixingtimes = 0
@@ -144,17 +144,20 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
    ! Reset age
    clock = 0.0_DP
    finterval = 1.0_DP / real(ntotcrat,kind=DP)
-   age     = user%interval * user%numintervals
-   if (age < 0._DP ) then
+   maxage = user%interval * user%numintervals
+   if (maxage < 0._DP ) then
       write(*,*) "MAJOR ERROR: Negative age!"
       stop
-   else if (age < 2330._DP) then
-      ageGa = util_t_from_scale(age,1e-11_DP,4.5_DP)
+   else if (maxage < 2330._DP) then
+      maxageGa = util_t_from_scale(maxage,1e-11_DP,4.5_DP)
    else
-      ageGa = 4.5_DP !util_t_from_scale only supports ages <4.5 Ga
+      maxageGa = 4.5_DP !util_t_from_scale only supports ages <4.5 Ga
    end if
-   age_resolution = ageGa / real(MAXAGEBINS)
+   age_resolution = maxageGa / real(MAXAGEBINS)
    write(*,*) "Age resolution: ", age_resolution, " Ga."
+   do i = 1,MAXAGEBINS
+      domain%age_bin_times(i) = maxageGa-(i*age_resolution)
+   end do
    domain%age_counter = 1
    oldGa = 0._DP
 
@@ -171,24 +174,32 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
       timestamp_old = real(curyear + real(icrater,kind=DP) / real(ntotcrat,kind=DP) * user%interval,kind=DP)
       icrater = icrater + 1
       crater%timestamp = real(curyear + real(icrater,kind=DP) / real(ntotcrat,kind=DP) * user%interval,kind=DP)
+      if (icrater .eq. 1) then
+         agemin = crater%timestamp * 0.9_DP
+      end if
       if (crater%timestamp < 2330._DP) then
          if (oldGa > 0._DP) then 
-            crater%timestampGa = util_t_from_scale(crater%timestamp,oldGa,ageGa)
+            crater%timestampGa = util_t_from_scale(maxage-crater%timestamp,agemin,oldGa)
          else
-            crater%timestampGa = util_t_from_scale(crater%timestamp,1e-11_DP,ageGa)
+            crater%timestampGa = util_t_from_scale(maxage-crater%timestamp,agemin,maxageGa)
          end if
       else
          crater%timestampGa = 4.5_DP
       end if
       pbarpos = nint(real(icrater) / real(ntotcrat) * PBARRES)
-      if (crater%timestampGa > (domain%age_counter*age_resolution)) then
-         domain%age_counter = domain%age_counter + 1
+      if (crater%timestampGa < domain%age_bin_times(domain%age_counter)) then
+         do i = domain%age_counter,MAXAGEBINS
+            if (crater%timestampGa > domain%age_bin_times(i)) then
+               domain%age_counter = i
+               exit
+            end if
+         end do
       end if 
       !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
+            write(message,*) "Real @ ", crater%timestampGa
             call io_updatePbar(message)
             user%testflag = .true.
             user%testimp = rclist(1, domain%rccount)
@@ -292,7 +303,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
                call ejecta_table_define(user,crater,domain,ejb,ejtble)
                !call ejecta_interpolate(crater,domain,crater%frad,ejb(1:ejtble),ejtble,crater%ejrim)
             end if
-            call ejecta_emplace(user,surf,crater,domain,ejb(1:ejtble),ejtble,ejbmass,age,age_resolution,&
+            call ejecta_emplace(user,surf,crater,domain,ejb(1:ejtble),ejtble,ejbmass,&
                ejecta_dem,nmeltsheet,vmeltsheet)
          else
             ejtble = 0
@@ -364,10 +375,9 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
                ! Do superdomain ray deposits
                ! Do sub-pixel craters vertical mixing
                if (user%doregotrack) then
-                     call crater_superdomain(user,surf,age,age_resolution,prod,nflux,domain,finterval)
+                     call crater_superdomain(user,surf,prod,nflux,domain,finterval)
                      call regolith_depth_model(user,domain,finterval,nflux,p)
                      call regolith_subcrater_mix(user,surf,domain,nflux,finterval,p)
-                     age = age - finterval * user%interval
                      nmixingtimes = nmixingtimes + 1
                end if 
 

src/crater/crater_superdomain.f90

@@ -19,7 +19,7 @@
 !  Notes       :  
 !
 !**********************************************************************************************************************************
-subroutine crater_superdomain(user,surf,age,age_resolution,prod,nflux,domain,finterval)
+subroutine crater_superdomain(user,surf,prod,nflux,domain,finterval)
    use module_globals
    use module_util
    use module_ejecta
@@ -30,8 +30,6 @@ subroutine crater_superdomain(user,surf,age,age_resolution,prod,nflux,domain,fin
    ! Arguments
    type(usertype),intent(in)                           :: user
    type(surftype),dimension(:,:),intent(inout)         :: surf
-   real(DP),intent(in)                                 :: age
-   real(DP),intent(in)                                 :: age_resolution
    real(DP),dimension(:,:),intent(in)                  :: prod,nflux 
    type(domaintype),intent(in)                         :: domain
    real(DP),intent(in)                                 :: finterval
@@ -173,7 +171,7 @@ subroutine crater_superdomain(user,surf,age,age_resolution,prod,nflux,domain,fin
                   if ((abs(xpi) > inc) .or. (abs(ypi) > inc)) cycle
                   if (ejisray(xpi,ypi) == 0) cycle  
                   call regolith_superdomain(user,crater,domain,surf(i,j)%regolayer,ejdistribution(xpi,ypi),&
-                       i,j,age,age_resolution,rm,depthb)
+                       i,j,rm,depthb)
                end do
             end do
 

src/crater/module_crater.f90

@@ -328,12 +328,10 @@ end function crater_profile_find_r_inner_wall
    end interface
 
    interface
-      subroutine crater_superdomain(user,surf,age,age_resolution,prod,nflux,domain,finterval)
+      subroutine crater_superdomain(user,surf,prod,nflux,domain,finterval)
       use module_globals
       type(usertype),intent(in)                           :: user
       type(surftype),dimension(:,:),intent(inout)         :: surf
-      real(DP),intent(in)                                 :: age
-      real(DP),intent(in)                                 :: age_resolution
       real(DP),dimension(:,:),intent(in)                  :: prod,nflux
       type(domaintype),intent(in)                         :: domain
       real(DP),intent(in)                                 :: finterval

src/ejecta/ejecta_emplace.f90

@@ -75,7 +75,7 @@
 !                The cutoff of ejecta thickness is still buggy.  
 !
 !**********************************************************************************************************************************
-subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_resolution,cumulative_elchange,&
+subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulative_elchange,&
    nmeltsheet,vmeltsheet)
    use module_globals
    use module_util
@@ -93,8 +93,6 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r
    integer(I4B),intent(in) :: ejtble
    type(ejbtype),dimension(:),intent(inout)    :: ejb
    real(DP),intent(in) :: deltaMtot
-   real(DP),intent(in)  :: age
-   real(DP),intent(in)  :: age_resolution
    real(DP),dimension(:,:),allocatable,intent(out) :: cumulative_elchange
    integer(I4B),intent(in) :: nmeltsheet
    real(DP),intent(out) :: vmeltsheet
@@ -336,7 +334,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r
 
 
                if (user%doregotrack .and. ebh>1.0e-8_DP) then
-                  call regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,lrad,ebh,rm,vsq,age,age_resolution,volm)
+                  call regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,lrad,ebh,rm,vsq,volm)
                   vmelt = vmelt + volm
                   !write(74,*) erad, surf(xpi,ypi)%regolayer%regodata%meltfrac
                end if
@@ -346,7 +344,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r
 
    if (totmelt > vmelt) then
       vmeltsheet = totmelt - vmelt
-   else !give the craters a melt sheet of 1mm
+   else !give the craters a melt sheet of 1m
       vmeltsheet = 1.0_DP * user%pix * user%pix * nmeltsheet
    end if
 

src/ejecta/module_ejecta.f90

@@ -25,7 +25,7 @@ module module_ejecta
 save
 
    interface
-      subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_resolution,cumulative_elchange,&
+      subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulative_elchange,&
          nmeltsheet,vmeltsheet)
       use module_globals
       implicit none
@@ -36,8 +36,6 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r
       integer(I4B),intent(in) :: ejtble
       type(ejbtype),dimension(:),intent(inout)   :: ejb
       real(DP),intent(in) :: deltaMtot
-      real(DP),intent(in)  :: age
-      real(DP),intent(in)  :: age_resolution
       real(DP),dimension(:,:),allocatable,intent(inout) :: cumulative_elchange
       integer(I4B),intent(in) :: nmeltsheet
       real(DP),intent(out) :: vmeltsheet

src/globals/module_globals.f90

@@ -170,6 +170,7 @@ module module_globals
    integer(I4B) :: tallycoverage  ! Estimated areal coverage of craters since the last tally
    integer(I4B) :: subpixelcoverage  ! Estimated areal coverage of craters since the last subpixel step
    real(DP)     :: rmsvel            ! Root mean square impact velocity
+   real(DP),dimension(MAXAGEBINS) :: age_bin_times ! Time in Ga to move to the next age bin
 end type domaintype 
 
 ! Derived data type for user input variables

src/regolith/module_regolith.f90

@@ -77,7 +77,7 @@ end subroutine regolith_transport
 
    interface 
       subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,lrad,ebh,&
-                 rm,vsq,age,age_resolution,volm)
+                 rm,vsq,volm)
       use module_globals 
       implicit none
       type(usertype),intent(in) :: user
@@ -88,7 +88,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,
       type(ejbtype),dimension(:),intent(in)   :: ejb
       real(DP),intent(in)          :: xp,yp,lrad,ebh
       integer(I4B),intent(in)      :: xpi,ypi
-      real(DP),intent(in)          :: rm, vsq, age, age_resolution
+      real(DP),intent(in)          :: rm, vsq
       real(DP),intent(inout)       :: volm
       end subroutine regolith_streamtube
    end interface
@@ -256,13 +256,11 @@ end subroutine regolith_depth_model
    end interface
 
    interface
-      subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints,melt)
+      subroutine regolith_melt_glass(user,crater,domain,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
       real(DP),intent(in)              :: rm
       real(DP),intent(in)              :: eradc
@@ -275,16 +273,14 @@ end subroutine regolith_melt_glass
    end interface
 
    interface 
-     subroutine regolith_superdomain(user,crater,domain,regolayer,ejdistribution,xpi,ypi,age,age_resolution,rm,depthb)
+     subroutine regolith_superdomain(user,crater,domain,regolayer,ejdistribution,xpi,ypi,rm,depthb)
      use module_globals
      type(usertype),intent(in)      :: user
      type(cratertype),intent(inout) :: crater
      type(domaintype),intent(in)    :: domain
      type(regodatatype),dimension(:),allocatable,intent(inout)       :: regolayer
      real(DP),intent(in)            :: ejdistribution 
      integer(I4B),intent(in)        :: xpi, ypi
-     real(DP),intent(in)            :: age
-     real(DP),intent(in)            :: age_resolution
      real(DP),intent(in)            :: rm
      real(DP),intent(in)            :: depthb
      end subroutine regolith_superdomain

src/regolith/regolith_melt_glass.f90

@@ -55,7 +55,7 @@
 !  Notes       :  
 !
 !**********************************************************************************************************************************
-subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints,melt)
+subroutine regolith_melt_glass(user,crater,domain,ebh,rm,eradc,lrad,deltar,newlayer,xmints,melt)
    use module_globals 
    use module_util
    use module_regolith, EXCEPT_THIS_ONE => regolith_melt_glass
@@ -65,8 +65,6 @@ subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad
    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
    real(DP),intent(in)              :: rm
    real(DP),intent(in)              :: eradc
@@ -183,12 +181,5 @@ subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad
       newlayer%age(domain%age_counter) = newlayer%meltvolume
    end if 
 
-   ! n_age = max(ceiling(age / age_resolution), 1)
-   ! if (lrad >= RAD_GP * crater%rad) then
-   !    newlayer%age(n_age) = melt / (user%pix * user%pix) 
-   ! else 
-   !    newlayer%age(n_age) = 0.0_SP   
-   ! end if
-
    return
 end subroutine regolith_melt_glass

src/regolith/regolith_streamtube.f90

@@ -56,7 +56,7 @@
 !  Notes       :  
 !
 !**********************************************************************************************************************************
-subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,lrad,ebh,rm,vsq,age,age_resolution,volm)
+subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,lrad,ebh,rm,vsq,volm)
    use module_globals 
    use module_util
    use module_regolith, EXCEPT_THIS_ONE => regolith_streamtube
@@ -71,7 +71,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,
    type(ejbtype),dimension(:),intent(in)   :: ejb
    real(DP),intent(in)          :: xp,yp,lrad,ebh
    integer(I4B),intent(in)      :: xpi,ypi
-   real(DP),intent(in)          :: rm, vsq, age, age_resolution
+   real(DP),intent(in)          :: rm, vsq
    real(DP),intent(inout)       :: volm
 
    ! Traversing a linked list 
@@ -224,7 +224,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,domain,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints,volm)
+   call regolith_melt_glass(user,crater,domain,ebh,rm,eradc,lrad,deltar,newlayer,xmints,volm)
    ! if (eradc>rm) then
    !    write(*,*) 'eradc > rm!'
    !    write(*,*) ebh, exp(ejb(k)%thick)

src/regolith/regolith_superdomain.f90

@@ -19,7 +19,7 @@
 !  Notes       :  
 !
 !**********************************************************************************************************************************
-subroutine regolith_superdomain(user,crater,domain,regolayer,ejdistribution,xpi,ypi,age,age_resolution,rm,depthb)
+subroutine regolith_superdomain(user,crater,domain,regolayer,ejdistribution,xpi,ypi,rm,depthb)
    use module_globals
    use module_util
    use module_regolith, EXCEPT_THIS_ONE => regolith_superdomain
@@ -32,8 +32,6 @@ subroutine regolith_superdomain(user,crater,domain,regolayer,ejdistribution,xpi,
    type(regodatatype),dimension(:),allocatable,intent(inout)            :: regolayer
    real(DP),intent(in)                   :: ejdistribution
    integer(I4B),intent(in)               :: xpi, ypi
-   real(DP),intent(in)                   :: age
-   real(DP),intent(in)                   :: age_resolution
    real(DP),intent(in)                   :: rm 
    real(DP),intent(in)                   :: depthb
 
@@ -60,7 +58,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,domain,age,age_resolution,ebh,rm,erad,lrad,deltar,newlayer,xmints,melt) 
+   call regolith_melt_glass(user,crater,domain,ebh,rm,erad,lrad,deltar,newlayer,xmints,melt) 
    call util_push_array(regolayer,newlayer)
 
    return

src/util/util_t_from_scale.f90

@@ -61,7 +61,7 @@ function util_t_from_scale(scale,start,finish) result(time)
 
     if (time .lt. 0) then
         write(*,*) "ERROR in util_t_from_scale: Maximum iterations reached!"
-        write(*,*) scale, maxiter
+        write(*,*) scale, start, finish, maxiter
     end if
 end function util_t_from_scale