diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index 6bef8315..e946e00b 100644 --- a/src/crater/crater_populate.f90 +++ b/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 diff --git a/src/crater/crater_superdomain.f90 b/src/crater/crater_superdomain.f90 index a153bdea..1d1f6092 100644 --- a/src/crater/crater_superdomain.f90 +++ b/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 diff --git a/src/crater/module_crater.f90 b/src/crater/module_crater.f90 index 307d0b78..c0dd5dfe 100644 --- a/src/crater/module_crater.f90 +++ b/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 diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index 72193d42..20706f1c 100644 --- a/src/ejecta/ejecta_emplace.f90 +++ b/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 diff --git a/src/ejecta/module_ejecta.f90 b/src/ejecta/module_ejecta.f90 index c36daf22..06f8d9e9 100644 --- a/src/ejecta/module_ejecta.f90 +++ b/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 diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index 2fa31110..d4ae2f8b 100644 --- a/src/globals/module_globals.f90 +++ b/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 diff --git a/src/regolith/module_regolith.f90 b/src/regolith/module_regolith.f90 index f3245e9e..2ab24cd9 100644 --- a/src/regolith/module_regolith.f90 +++ b/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,7 +273,7 @@ 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 @@ -283,8 +281,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 end subroutine regolith_superdomain diff --git a/src/regolith/regolith_melt_glass.f90 b/src/regolith/regolith_melt_glass.f90 index 10b69fd6..b7d63d6e 100644 --- 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,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 diff --git a/src/regolith/regolith_streamtube.f90 b/src/regolith/regolith_streamtube.f90 index 4088d894..28c28d08 100644 --- a/src/regolith/regolith_streamtube.f90 +++ b/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) diff --git a/src/regolith/regolith_superdomain.f90 b/src/regolith/regolith_superdomain.f90 index 3c3d76ce..ce73a3a6 100644 --- a/src/regolith/regolith_superdomain.f90 +++ b/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 diff --git a/src/util/util_t_from_scale.f90 b/src/util/util_t_from_scale.f90 index 216c7154..e45fcfd0 100644 --- a/src/util/util_t_from_scale.f90 +++ b/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