From 145881aeaa0d8c2a441f01c3279dac0ef4cf6107 Mon Sep 17 00:00:00 2001 From: Austin Blevins Date: Fri, 10 Nov 2023 16:33:41 -0500 Subject: [PATCH] Attempt to add length relations based on fits to Jake's data; currently bugged --- src/crater/crater_subpixel_diffusion.f90 | 9 ++++++--- src/crater/crater_superdomain.f90 | 8 ++++++-- src/ejecta/ejecta_emplace.f90 | 10 ++++++++-- src/ejecta/ejecta_ray_pattern.f90 | 9 ++++++--- src/ejecta/ejecta_ray_pattern_func.f90 | 6 ++++-- src/ejecta/module_ejecta.f90 | 6 ++++-- 6 files changed, 34 insertions(+), 14 deletions(-) diff --git a/src/crater/crater_subpixel_diffusion.f90 b/src/crater/crater_subpixel_diffusion.f90 index 05c9e05f..cdba3408 100644 --- a/src/crater/crater_subpixel_diffusion.f90 +++ b/src/crater/crater_subpixel_diffusion.f90 @@ -51,13 +51,16 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin) 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) :: 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)) + real(DP) :: l1 + rray = 11.95*crater%frad**1.32 + l1 = 5.32*crater%frad**1.27 ! Create box for soften calculation (will be no bigger than the grid itself) do j = 0,user%gridsize + 1 @@ -184,7 +187,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,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution) + call ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution,l1) ! Loop over affected matrix area !!$OMP PARALLEL DO DEFAULT(SHARED) IF(inc > INCPAR) & !!$OMP FIRSTPRIVATE(jmin,jmax,imin,imax) & diff --git a/src/crater/crater_superdomain.f90 b/src/crater/crater_superdomain.f90 index 8cc47ff9..65c94e09 100644 --- a/src/crater/crater_superdomain.f90 +++ b/src/crater/crater_superdomain.f90 @@ -52,12 +52,16 @@ subroutine crater_superdomain(user,surf,prod,nflux,domain,finterval) real(SP) :: gglass, glrad ! Crater ray parameters - real(DP) :: rray = 48_DP ! "L16" in Minton et al. (2019) + 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)) + real(DP) :: l1 + + rray = 11.95*crater%frad**1.32 + l1 = 5.32*crater%frad**1.27 ! Create box for soften calculation (will be no bigger than the grid itself) do j = 0,user%gridsize + 1 @@ -156,7 +160,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,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution) + call ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution,l1) ! Now if doregotrack is on, do melt zone calculation if (user%doregotrack) call regolith_melt_zone_superdomain(user,crater,domain,rm,depthb) diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index cc5f73d0..bbbea701 100644 --- a/src/ejecta/ejecta_emplace.f90 +++ b/src/ejecta/ejecta_emplace.f90 @@ -129,12 +129,18 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ real(SP) :: age_mean ! Crater ray parameters - real(DP) :: rray = 48_DP ! "L16" in Minton et al. (2019) + 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)) + real(DP) :: l1 + + + rray = 11.95*crater%frad**1.32 + l1 = 5.32*crater%frad**1.27 + write(*,*) "L16 = ", rray, "; L1 = ", l1 ! Executable code @@ -190,7 +196,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,rray,Nraymax+i,fpeak,rayp,rayq,rayfmult,tempdiff,tempej) + call ejecta_ray_pattern(user,surf,crater,inc,-inc,inc,-inc,inc,rray,Nraymax+i,fpeak,rayp,rayq,rayfmult,tempdiff,tempej,l1) diffdistribution(:,:) = diffdistribution(:,:) + frayreduction**(i-1) * tempdiff(:,:) ejdistribution(:,:) = ejdistribution(:,:) + frayreduction**(i-1) * tempej(:,:) end do diff --git a/src/ejecta/ejecta_ray_pattern.f90 b/src/ejecta/ejecta_ray_pattern.f90 index f9da323f..7a25bd50 100644 --- a/src/ejecta/ejecta_ray_pattern.f90 +++ b/src/ejecta/ejecta_ray_pattern.f90 @@ -31,7 +31,7 @@ !*** !********************************************************************************************************************************** -subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution) +subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution,l1) use module_globals use module_util use module_io @@ -49,6 +49,7 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpea 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 + real(DP),intent(in) :: l1 ! Internal variables integer(I4B) :: nrays,i,j,k,n,nef,incsq,iradsq,xpi,ypi,ejpxsq @@ -87,6 +88,8 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpea call random_number(rn) ! randomize the ray orientation rmax = user%ejecta_truncation rmin = crater%continuous / crater%frad + !rmax = rray / crater%frad + !rmin = l1 / crater%frad crater%fe = 10.0_DP ! Estimate the equivalent degradation radius !ejdistribution = 0.0_DP !diffdistribution = 0.0_DP @@ -108,8 +111,8 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpea 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,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.) + diffdistribution(i,j) = areafrac * ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,rayp,rayq,rayfmult,l1,.false.) + ejdistribution(i,j) = areafrac * ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,rayp,rayq,rayfmult,l1,.true.) end do end do !!$OMP END PARALLEL DO diff --git a/src/ejecta/ejecta_ray_pattern_func.f90 b/src/ejecta/ejecta_ray_pattern_func.f90 index 61e173fa..d5d81354 100644 --- a/src/ejecta/ejecta_ray_pattern_func.f90 +++ b/src/ejecta/ejecta_ray_pattern_func.f90 @@ -27,7 +27,7 @@ !*** !********************************************************************************************************** -function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,rayp,rayq,rayfmult,ej) result(ans) +function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,rayp,rayq,rayfmult,l1,ej) result(ans) use module_globals use module_ejecta, EXCEPT_THIS_ONE => ejecta_ray_pattern_func implicit none @@ -36,6 +36,7 @@ function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,ra real(DP),dimension(:),intent(in) :: thetari real(DP),intent(in) :: rray, fpeak, rayp, rayfmult integer(I4B),intent(in) :: Nraymax, rayq + real(DP),intent(in) :: l1 logical,intent(in) :: ej real(DP) :: a,c real(DP) :: thetar,rw,rw0,rw1 @@ -45,7 +46,8 @@ function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,ra !minray = rmin * 3 !"L1" in Minton et al. (2019) - minray = rray / 2 + !minray = rray / 2 + minray = l1 if (r > rmax) then ans = 0._DP diff --git a/src/ejecta/module_ejecta.f90 b/src/ejecta/module_ejecta.f90 index 2c21a963..e2d8747a 100644 --- a/src/ejecta/module_ejecta.f90 +++ b/src/ejecta/module_ejecta.f90 @@ -43,7 +43,7 @@ end subroutine ejecta_emplace end interface interface - subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution) + subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution,l1) use module_globals implicit none type(usertype),intent(in) :: user @@ -54,12 +54,13 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpea 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 + real(DP),intent(in) :: l1 end subroutine ejecta_ray_pattern end interface interface - function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,rayp,rayq,rayfmult,ej) result(ans) + function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,rayp,rayq,rayfmult,l1,ej) result(ans) use module_globals implicit none real(DP) :: ans @@ -67,6 +68,7 @@ function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,ra real(DP),dimension(:),intent(in) :: thetari real(DP),intent(in) :: rray, fpeak, rayp, rayfmult integer(I4B),intent(in) :: Nraymax, rayq + real(DP),intent(in) :: l1 logical,intent(in) :: ej end function ejecta_ray_pattern_func end interface