diff --git a/src/Makefile.am b/src/Makefile.am index 2ea4799b..2e8c1902 100755 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -43,6 +43,9 @@ util/util_traverse_pop_array.f90\ util/util_destroy_list.f90\ util/util_init_array.f90\ util/util_perlin_noise.f90\ +util/util_tscale.f90\ +util/util_t_from_scale.f90\ +util/util_npf_timefunc.f90\ io/io_read_const.f90\ io/io_get_token.f90\ io/io_input.f90\ diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index 118cdd5d..409dfe58 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -93,7 +93,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt real(SP),dimension(user%gridsize, user%gridsize) :: agetop real(SP),dimension(60) :: agetot type(regolisttype),pointer :: current => null() - real(DP) :: age_resolution + real(DP) :: age_resolution, ageGa, oldGa integer(I4B) :: age_counter nmixingtimes = 0 @@ -145,8 +145,18 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt clock = 0.0_DP finterval = 1.0_DP / real(ntotcrat,kind=DP) age = user%interval * user%numintervals - age_resolution = age / real(MAXAGEBINS) + if (age < 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 + ageGa = 4.5_DP !util_t_from_scale only supports ages <4.5 Ga + end if + age_resolution = ageGa / real(MAXAGEBINS) + write(*,*) "Age resolution: ", age_resolution, " Ga." domain%age_counter = 1 + oldGa = 0._DP ! Reset coverage map domain%tallycoverage = 0 @@ -158,11 +168,20 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt 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) + 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=SP) + crater%timestamp = real(curyear + real(icrater,kind=DP) / real(ntotcrat,kind=DP) * user%interval,kind=DP) + if (crater%timestamp < 2330._DP) then + if (oldGa > 0._DP) then + crater%timestampGa = util_t_from_scale(crater%timestamp,oldGa,ageGa) + else + crater%timestampGa = util_t_from_scale(crater%timestamp,1e-11_DP,ageGa) + end if + else + crater%timestampGa = 4.5_DP + end if pbarpos = nint(real(icrater) / real(ntotcrat) * PBARRES) - if (crater%timestamp > (domain%age_counter*age_resolution)) then + if (crater%timestampGa > (domain%age_counter*age_resolution)) then domain%age_counter = domain%age_counter + 1 end if !if in quasiMC mode: check to see if it's time for a real crater diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index e105c279..2fa31110 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -102,7 +102,8 @@ module module_globals real(DP) :: imp,imprad,impvel,sinimpang,impmass ! Impactor properties ! Real domain properties real(SP) :: xl,yl ! Crater center in simulated surface size units - real(SP) :: timestamp ! The formation time of the crater + real(DP) :: timestamp ! The formation time of the crater + real(DP) :: timestampGa ! Crater formation time in Ga real(DP) :: rad ! Transient radius real(DP) :: grad ! Strengthless material transient crater radius real(DP) :: frad ! Final crater radius diff --git a/src/util/module_util.f90 b/src/util/module_util.f90 index 26ecb045..a15dfaf6 100644 --- a/src/util/module_util.f90 +++ b/src/util/module_util.f90 @@ -278,5 +278,29 @@ function util_perlin_noise(x,y,z) result (noise) end function util_perlin_noise end interface +interface + function util_npf_timefunc(T) result(N1) + use module_globals + real(DP), intent(in) :: T + real(DP) :: N1 + end function util_npf_timefunc +end interface + +interface + function util_tscale(t) result(tscale) + use module_globals + real(DP), intent(in) :: t + real(DP) :: tscale + end function util_tscale +end interface + +interface + function util_t_from_scale(scale,start,finish) result(time) + use module_globals + real(DP), intent(in) :: scale, start, finish + real(DP) :: time + end function util_t_from_scale +end interface + end module diff --git a/src/util/util_npf_timefunc.f90 b/src/util/util_npf_timefunc.f90 new file mode 100644 index 00000000..cd4c287d --- /dev/null +++ b/src/util/util_npf_timefunc.f90 @@ -0,0 +1,27 @@ +!********************************************************************************************************************************** +! +! Unit Name : util_npf_timefunc +! Unit Type : subroutine +! Project : CTEM +! Language : Fortran 2003 +! +! Description : The time function of the lunar Neukum production function (Neukum et al., 2001) +! +! Input +! Arguments : +! +! Output +! Arguments : +! +! +! Notes : +! +!********************************************************************************************************************************** +function util_npf_timefunc(T) result(N1) + use module_globals + use module_util, EXCEPT_THIS_ONE => util_npf_timefunc + real(DP), intent(in) :: T + real(DP) :: N1 + + N1 = 5.44e-14 * (exp(6.93*T)-1) + 8.17e-4*T +end function util_npf_timefunc diff --git a/src/util/util_t_from_scale.f90 b/src/util/util_t_from_scale.f90 new file mode 100644 index 00000000..13d17953 --- /dev/null +++ b/src/util/util_t_from_scale.f90 @@ -0,0 +1,68 @@ +!********************************************************************************************************************************** +! +! Unit Name : util_t_from_scale +! Unit Type : subroutine +! Project : CTEM +! Language : Fortran 2003 +! +! Description : converts age from "interval time" to true time in Ga based on the lunar Neukum production function (Neukum et al., 2001) +! +! Input +! Arguments : +! +! Output +! Arguments : +! +! +! Notes : This is a Fortran port of the function "T_from_scale" from craterproduction.py. It uses the bisect method so is simple but inefficient. +! +!********************************************************************************************************************************** +function util_t_from_scale(scale,start,finish) result(time) + use module_globals + use module_util, EXCEPT_THIS_ONE => util_t_from_scale + real(DP), intent(in) :: scale, start, finish + real(DP) :: time + + real(DP) :: tol = 1e-11_DP + integer(I4B) :: maxiter = 1000 + integer(I4B) :: i + real(DP) :: a, b, c + + + a = start + b = finish + time = -1. + i = 0 + temp = 0._DP + + if (abs(temp-scale)0) then + b = c + i = i + 1 + else if ((temp-scale)<0) then + a = c + i = i + 1 + else if (abs(a-b) util_tscale + real(DP), intent(in) :: t + real(DP) :: tscale + + real(DP) :: N1 + real(DP) :: CSFD = 0.0008173348179780709 !this is the result of the shape function at T=1 Ga for the lunar NPF + + N1 = util_npf_timefunc(t) + + tscale = N1 / CSFD +end function util_tscale + + +