Tweaked layering parameters; still unrealistic, but better than before

src/crater/crater_subpixel_diffusion.f90

@@ -49,6 +49,14 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin)
    integer(I4B),dimension(:,:),allocatable :: ejisray
    real(DP) :: xbar,ybar,dD,xp,yp,areafrac,krad,supersize,lrad
    integer(I8B),dimension(user%gridsize,user%gridsize) :: Ngrid
+
+   ! Crater ray parameters
+   real(DP) :: rray = 48_DP ! "L16" in Minton et al. (2019)
+   integer(I4B) :: Nraymax = 14
+   real(DP) :: fpeak = 8000_DP ! narrow ray: rw0 propto 1/4
+   real(DP) :: rayp = 2.0_DP 
+   integer(I4B) :: rayq = 4
+   real(DP) :: rayfmult = (5)**(-4.0_DP / (1.2_DP))
 
 
    ! Create box for soften calculation (will be no bigger than the grid itself)
@@ -176,7 +184,7 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin)
 
             allocate(diffdistribution(imin:imax,jmin:jmax))
             allocate(ejdistribution(imin:imax,jmin:jmax))
-            call ejecta_ray_pattern(user,surf,crater,inc,imin,imax,jmin,jmax,diffdistribution,ejdistribution)
+            call ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution)
             ! Loop over affected matrix area
             !!$OMP PARALLEL DO DEFAULT(SHARED) IF(inc > INCPAR) &
             !!$OMP FIRSTPRIVATE(jmin,jmax,imin,imax) &

src/crater/crater_superdomain.f90

@@ -51,6 +51,14 @@ subroutine crater_superdomain(user,surf,prod,nflux,domain,finterval)
    real(DP) :: rm, depthb, maxgcrat, maxgcratkm
    real(SP) :: gglass, glrad
 
+      ! Crater ray parameters
+   real(DP) :: rray = 48_DP ! "L16" in Minton et al. (2019)
+   integer(I4B) :: Nraymax = 14
+   real(DP) :: fpeak = 8000_DP ! narrow ray: rw0 propto 1/4
+   real(DP) :: rayp = 2.0_DP 
+   integer(I4B) :: rayq = 4
+   real(DP) :: rayfmult = (5)**(-4.0_DP / (1.2_DP))
+
    ! Create box for soften calculation (will be no bigger than the grid itself)
    do j = 0,user%gridsize + 1
       do i = 0,user%gridsize + 1
@@ -148,7 +156,7 @@ subroutine crater_superdomain(user,surf,prod,nflux,domain,finterval)
          allocate(diffdistribution(xi:xf,yi:yf))
          allocate(ejisray(xi:xf,yi:yf))
          ! Now generate ray pattern
-         call ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,diffdistribution,ejdistribution)
+         call ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution)
          ! Now if doregotrack is on, do melt zone calculation
          if (user%doregotrack) call regolith_melt_zone_superdomain(user,crater,domain,rm,depthb)
 

src/ejecta/ejecta_emplace.f90

@@ -128,6 +128,14 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ
    ! Age
    real(SP) :: age_mean
 
+   ! Crater ray parameters
+   real(DP) :: rray = 48_DP ! "L16" in Minton et al. (2019)
+   integer(I4B) :: Nraymax = 12
+   real(DP) :: fpeak = 8000_DP ! narrow ray: rw0 propto 1/4
+   real(DP) :: rayp = 2.0_DP 
+   integer(I4B) :: rayq = 4
+   real(DP) :: rayfmult = (5)**(-4.0_DP / (1.2_DP))
+
 
    ! Executable code
 
@@ -182,7 +190,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ
 
    ! *************************** Layered Ejecta Rays *****************************!
    do i=1,Npatt
-      call ejecta_ray_pattern(user,surf,crater,inc,-inc,inc,-inc,inc,tempdiff,tempej)
+      call ejecta_ray_pattern(user,surf,crater,inc,-inc,inc,-inc,inc,rray,Nraymax+i,fpeak,rayp,rayq,rayfmult,tempdiff,tempej)
       diffdistribution(:,:) = diffdistribution(:,:) + frayreduction**(i-1) * tempdiff(:,:)
       ejdistribution(:,:) = ejdistribution(:,:) + frayreduction**(i-1) * tempej(:,:)
    end do

src/ejecta/ejecta_ray_pattern.f90

@@ -31,7 +31,7 @@
 !***
 
 !**********************************************************************************************************************************
-subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,diffdistribution,ejdistribution)
+subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution)
    use module_globals
    use module_util
    use module_io
@@ -45,6 +45,8 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,diffdistribution,
    type(surftype),dimension(:,:),intent(in) :: surf
    type(cratertype),intent(inout) :: crater
    integer(I4B),intent(in) :: inc,xi,xf,yi,yf
+   real(DP),intent(in) :: rray, fpeak, rayp, rayfmult
+   integer(I4B),intent(in) :: Nraymax, rayq
    real(DP),dimension(xi:xf,yi:yf),intent(out) :: diffdistribution
    real(DP),dimension(xi:xf,yi:yf),intent(out) :: ejdistribution
 
@@ -106,8 +108,8 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,diffdistribution,
             areafrac = util_area_intersection(user%ejecta_truncation * crater%frad,xbar,ybar,user%pix) 
             r = sqrt(xbar**2 + ybar**2) / crater%frad
             theta = mod(atan2(ybar,xbar) + pi + rn * 2 * pi,2 * pi)
-            diffdistribution(i,j) = areafrac * ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,.false.) 
-            ejdistribution(i,j) = areafrac * ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,.true.) 
+            diffdistribution(i,j) = areafrac * ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,rayp,rayq,rayfmult,.false.) 
+            ejdistribution(i,j) = areafrac * ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,rayp,rayq,rayfmult,.true.) 
          end do
       end do
       !!$OMP END PARALLEL DO

src/ejecta/ejecta_ray_pattern_func.f90

@@ -27,13 +27,15 @@
 !***
 
 !**********************************************************************************************************
-function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,ej) result(ans)
+function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,rayp,rayq,rayfmult,ej) result(ans)
     use module_globals
     use module_ejecta, EXCEPT_THIS_ONE => ejecta_ray_pattern_func
     implicit none
     real(DP) :: ans
     real(DP),intent(in) :: r,rmin,rmax,theta
     real(DP),dimension(:),intent(in) :: thetari
+    real(DP),intent(in) :: rray, fpeak, rayp, rayfmult
+    integer(I4B),intent(in) :: Nraymax, rayq
     logical,intent(in) :: ej
     real(DP) :: a,c
     real(DP) :: thetar,rw,rw0,rw1
@@ -42,8 +44,8 @@ function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,ej) result(ans)
     real(DP) :: tmp
 
 
-    minray = rmin * 3 !"L1" in Minton et al. (2019)
-    !minray = 11.0_DP 
+    !minray = rmin * 3 !"L1" in Minton et al. (2019)
+    minray = rray / 2
 
     if (r > rmax) then
        ans = 0._DP

src/ejecta/module_ejecta.f90

@@ -43,25 +43,30 @@ end subroutine ejecta_emplace
    end interface
 
    interface
-      subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,diffdistribution,ejdistribution)
+      subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution)
       use module_globals
       implicit none
       type(usertype),intent(in) :: user
       type(surftype),dimension(:,:),intent(in) :: surf
       type(cratertype),intent(inout) :: crater
       integer(I4B),intent(in) :: inc,xi,xf,yi,yf
+      real(DP),intent(in) :: rray, fpeak, rayp, rayfmult
+      integer(I4B),intent(in) :: Nraymax, rayq
       real(DP),dimension(xi:xf,yi:yf),intent(out) :: diffdistribution
       real(DP),dimension(xi:xf,yi:yf),intent(out) :: ejdistribution
+
       end subroutine ejecta_ray_pattern
    end interface
 
    interface
-      function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,ej) result(ans)
+      function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,rayp,rayq,rayfmult,ej) result(ans)
       use module_globals
       implicit none
       real(DP) :: ans
       real(DP),intent(in) :: r,rmin,rmax,theta
       real(DP),dimension(:),intent(in) :: thetari
+      real(DP),intent(in) :: rray, fpeak, rayp, rayfmult
+      integer(I4B),intent(in) :: Nraymax, rayq
       logical,intent(in) :: ej
       end function ejecta_ray_pattern_func
    end interface

src/globals/module_globals.f90

@@ -383,12 +383,4 @@ module module_globals
 real(DP),parameter :: GFAC = 0.487_DP       ! gravitational acceleration exponent
 real(DP),parameter :: DFAC = 0.556_DP       ! impact distance exponent
 
-! Crater ray parameters
-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
-real(DP),parameter :: rayfmult  = (5)**(-4.0_DP / (1.2_DP))  
-
 end module module_globals