Add ports of craterproduction.py functions to convert to true age internally

MintonGroup · Jul 19, 2023 · aa92e35 · aa92e35
1 parent 87dd833
commit aa92e35
Show file tree

Hide file tree

Showing 7 changed files with 183 additions and 6 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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)<tol) then
+        ans = 0.0_DP
+    else
+        do while(i .lt. maxiter)
+            c = (a+b)/2
+            temp = util_tscale(c)
+
+            if (abs(temp-scale)<tol) then
+                time = c
+                exit
+            else if ((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)<tol) then
+                write(*,*) "ERROR in util_t_from_scale: Convergence failed!"
+                write(*,*) scale, start, finish
+                exit
+            end if
+        end do
+    end if
+
+    if (ans .lt. 0) then
+        write(*,*) "ERROR in util_t_from_scale: Maximum iterations reached!"
+        write(*,*) scale, maxiter
+    end if
+end function util_t_from_scale
+
+
diff --git a/src/util/util_tscale.f90 b/src/util/util_tscale.f90
@@ -0,0 +1,35 @@
+!**********************************************************************************************************************************
+!
+!  Unit Name   : util_tscale
+!  Unit Type   : subroutine
+!  Project     : CTEM
+!  Language    : Fortran 2003
+!
+!  Description : converts age from true time in Ga to "interval time" 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_scale" from craterproduction.py
+!
+!**********************************************************************************************************************************
+function util_tscale(t) result(tscale)
+    use module_globals
+    use module_util, EXCEPT_THIS_ONE => 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
+
+
+