diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index 20706f1c..44372d60 100644 --- 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 index 8f83a704..fa4b712f 100644 --- 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 index f0014af0..3a325840 100644 --- 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