From 631e291fc9ec9abcd3844a8acbbeca8a41c1234d Mon Sep 17 00:00:00 2001 From: Austin Blevins Date: Fri, 3 Nov 2023 14:44:43 -0400 Subject: [PATCH] first pass at a more realistic ejecta ray model --- src/ejecta/ejecta_emplace.f90 | 38 ++++++++++++++++++++++++++++++++--- 1 file changed, 35 insertions(+), 3 deletions(-) diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index 44372d60..a7a92e98 100644 --- a/src/ejecta/ejecta_emplace.f90 +++ b/src/ejecta/ejecta_emplace.f90 @@ -104,10 +104,13 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ real(DP) :: xp,yp,fradsq,fradpxsq,radsq,ebh,ejdissq,ejbmass,fmasscons,areafrac,xbar,ybar,krad,kdiffmax real(DP),dimension(:,:),allocatable :: big_cumulative_elchange,kdiff,big_kdiff,cel,big_cel integer(I4B),dimension(:,:,:),allocatable :: indarray,big_indarray - real(DP),dimension(:,:),allocatable :: ejdistribution,diffdistribution + real(DP),dimension(:,:),allocatable :: ejdistribution,diffdistribution,maxdiff, maxej + real(DP),dimension(:,:,:),allocatable :: tempdiff,tempej integer(I4B) :: bigi,bigj,maxhits,nin,nnot,dradsq character(len=MESSAGESIZE) :: message ! message for the progress bar real(DP) :: vmelt, totmelt, volm + real(DP), dimension(:), allocatable :: freduction + integer(I4B) :: Npatt = 6 !Number of times to call ray pattern ! Ray mixing model variables @@ -171,12 +174,41 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ allocate(ejdistribution(-inc:inc,-inc:inc)) allocate(diffdistribution(-inc:inc,-inc:inc)) - !the other two temporary arrays get allocated here + allocate(tempdiff(1:Npatt,-inc:inc,-inc:inc)) + allocate(tempej(1:Npatt,-inc:inc,-inc:inc)) + allocate(freduction(Npatt)) + allocate(maxej(-inc:inc,-inc:inc)) + allocate(maxdiff(-inc:inc,-inc:inc)) + ! Set freduction parameters (Manually) <--for now this assumes Npatt=6 + freduction(1) = 0.8_DP + freduction(2) = 0.1_DP + freduction(3) = 0.04_DP + freduction(4) = 0.03_DP + freduction(5) = 0.02_DP + freduction(6) = 0.01_DP ! *************************** Continuous Ejecta Formula *****************************! - call ejecta_ray_pattern(user,surf,crater,inc,-inc,inc,-inc,inc,diffdistribution,ejdistribution) + do i=1,Npatt + call ejecta_ray_pattern(user,surf,crater,inc,-inc,inc,-inc,inc,diffdistribution,ejdistribution) + tempej(i,:,:) = ejdistribution(:,:) + tempdiff(i,:,:) = ejdistribution(:,:) + end do + + maxdiff(:,:) = maxval(tempdiff(:,:,:)) + maxej(:,:) = maxval(tempej(:,:,:)) + + ejdistribution(:,:) = 0.0_DP + diffdistribution(:,:) = 0.0_DP + do i=1,Npatt + ejdistribution(:,:) = ejdistribution(:,:) + (freduction(i) * tempej(i,:,:)) + diffdistribution(:,:) = diffdistribution(:,:) + (freduction(i) * tempdiff(i,:,:)) + end do + + ejdistribution(:,:) = ejdistribution(:,:) / maxej(:,:) + diffdistribution(:,:) = diffdistribution(:,:) / maxdiff(:,:) + deallocate(maxdiff,maxej,tempdiff,tempej,freduction) allocate(cumulative_elchange(-inc:inc,-inc:inc)) allocate(cel(-inc:inc,-inc:inc))