Re-arranged array allocations in preparation for a more realistic eje…

…cta ray calculation
MintonGroup · Nov 2, 2023 · 649fb61 · 649fb61
1 parent c8356c2
commit 649fb61
Show file tree

Hide file tree

Showing 3 changed files with 22 additions and 11 deletions.
diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90
@@ -168,21 +168,25 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ
        end if
    endif
 
+
+   allocate(ejdistribution(-inc:inc,-inc:inc))
+   allocate(diffdistribution(-inc:inc,-inc:inc))
+   !the other two temporary arrays get allocated here
+
+
+   ! *************************** Continuous Ejecta Formula  *****************************!
+   call ejecta_ray_pattern(user,surf,crater,inc,-inc,inc,-inc,inc,diffdistribution,ejdistribution)
+
+
    allocate(cumulative_elchange(-inc:inc,-inc:inc))
    allocate(cel(-inc:inc,-inc:inc))
    allocate(kdiff(-inc:inc,-inc:inc))
    allocate(indarray(2,-inc:inc,-inc:inc))
-   allocate(ejdistribution(-inc:inc,-inc:inc))
-   allocate(diffdistribution(-inc:inc,-inc:inc))
-
    cumulative_elchange = 0.0_DP
    kdiff = 0.0_DP
    indarray = inc - 1 ! initialize this array to point to a corner (this should have 0 elevation change since we're only doing work
                 ! within a circle of radius irad
 
-   ! *************************** Continuous Ejecta Formula  *****************************!
-   call ejecta_ray_pattern(user,surf,crater,inc,-inc,inc,-inc,inc,diffdistribution,ejdistribution)
-
    ejbmass = 0.0_DP
    nin = 0
    nnot = 0

diff --git a/src/ejecta/ejecta_ray_pattern_func.f90 b/src/ejecta/ejecta_ray_pattern_func.f90
@@ -39,9 +39,11 @@ function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,ej) result(ans)
     real(DP) :: thetar,rw,rw0,rw1
     real(DP) :: f,rtrans,length,rpeak,minray,FF
     integer(I4B) :: n,i
+    real(DP) :: tmp
 
 
-    minray = rmin * 3
+    minray = rmin * 3 !"L1" in Minton et al. (2019)
+    !minray = 11.0_DP 
 
     if (r > rmax) then
        ans = 0._DP
@@ -55,7 +57,12 @@ function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,ej) result(ans)
        rw0 = rmin * pi / Nraymax / 2
        rw1 = 2 * pi / Nraymax
        rw = rw0 * (1._DP - (1.0_DP - rw1 / rw0) * exp(1._DP - (r / rmin)**2)) ! equation 40 Minton et al. 2019
-       n = max(min(floor((Nraymax**rayp - (Nraymax**rayp - 1) * log(r/minray) / log(rray/minray))**(1._DP/rayp)),Nraymax),1) ! Exponential decay of ray number with distance
+       tmp = (Nraymax**rayp - (Nraymax**rayp - 1) * log(r/minray) / log(rray/minray))
+       if (tmp < 0.0_DP) then
+         n = Nraymax ! "Nrays" in Minton et al. (2019)
+       else
+         n = max(min(floor((Nraymax**rayp - (Nraymax**rayp - 1) * log(r/minray) / log(rray/minray))**(1._DP/rayp)),Nraymax),1) ! Exponential decay of ray number with distance
+       end if
        ans = 0._DP
        rtrans = r - 1.0_DP
        c = rw / r

diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90
@@ -221,7 +221,7 @@ module module_globals
 
    ! Ejecta softening variables
    logical           :: dosoftening  ! Set T to use the extra crater softening model
-   real(DP)          :: ejecta_truncation ! Set the number of crater diameters to truncate the ejecta
+   real(DP)          :: ejecta_truncation ! Set the number of crater radii to truncate the ejecta
    logical           :: dorays       ! Set T to use ray model
    logical           :: superdomain  ! Set T to include the superdomain
 
@@ -384,8 +384,8 @@ module module_globals
 real(DP),parameter :: DFAC = 0.556_DP       ! impact distance exponent
 
 ! Crater ray parameters
-real(DP),parameter :: rray = 24_DP ! I think this is "L16" in Minton et al. (2019)
-integer(I4B),parameter :: Nraymax = 16 ! "Nrays" in Minton et al. (2019)
+real(DP),parameter :: rray = 48_DP ! "L16" in Minton et al. (2019)
+integer(I4B),parameter :: Nraymax = 14 ! "Nrays" in Minton et al. (2019)
 real(DP),parameter :: fpeak = 8000_DP ! narrow ray: rw0 propto 1/4
 real(DP),parameter :: rayp = 2.0_DP 
 integer(I4B),parameter :: rayq = 4