Skip to content

Commit

Permalink
Reverted back to symmetric ejecta blanket
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Feb 1, 2017
1 parent b2e24bf commit 0c6095b
Showing 1 changed file with 34 additions and 58 deletions.
92 changes: 34 additions & 58 deletions src/ejecta/ejecta_emplace.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot)
real(DP) :: lrad,lrad0,lradsq,cdepth,distance,erad,craterslope,landslope,baseline,inneredge,outeredge,lradtrue
integer(I4B),parameter :: MAXLOOP = 10 ! Maximum number of times to loop the ejecta angle correction calculation
integer(I4B) :: xpi,ypi,i,j,k,n,inc,incsq,iradsq
real(DP) :: xp,yp,fradsq,radsq,ebh,ejdissq,continuous,ejbmass
real(DP) :: xp,yp,fradsq,radsq,ebh,ejdissq,continuous,ejbmass,fmasscons
real(DP),dimension(:,:),allocatable :: cumulative_elchange,big_cumulative_elchange
integer(I4B),dimension(:,:,:),allocatable :: indarray,big_indarray
integer(I4B) :: bigi,bigj
Expand Down Expand Up @@ -141,7 +141,9 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot)
if (crater%ejdis <= crater%rad) return

! determine area to effect
continuous = 2.348_DP * crater%frad**(1.006_DP)
inc = max(min(nint(1.5_DP * crater%ejdis / user%pix) + 1,PBCLIM*user%gridsize),1)
if (.not.user%discontinuous) inc = int(continuous / user%pix) + 1
crater%maxinc = max(crater%maxinc,inc)
radsq = crater%rad**2
fradsq = crater%frad**2
Expand All @@ -165,7 +167,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot)
! within a circle of radius irad

! *************************** Continuous Ejecta Formula *****************************!
continuous = 2.3_DP * crater%frad**(1.006_DP)
!^^^^^^^^^^^^^^^^^^^^^^^^

! *************************** Superformula Ray Model ************************************!
! *************************** Part I. Spoke and Skinny Ray Model ************************************!
Expand Down Expand Up @@ -220,7 +222,6 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot)
! inside or outside a ray?
theta = atan2(j * 1._DP,i * 1._DP)
mag = ( (abs(cos(nrays * theta / 4.0_DP)))**n2 + (abs(sin(nrays * theta / 4.0_DP)))**n2 )**(-1.0_DP/n1) !&

magf = rayf * ( (abs(cos(nfrays * theta / 4.0_DP)))**n2f + (abs(sin(nfrays * theta / 4.0_DP)))**n2f )**(-1.0_DP/n1f)
lradp = continuous * mag
lradf = continuous * magf
Expand Down Expand Up @@ -251,13 +252,13 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot)
inneredge = crater%frad + 0.5_DP * domain%ejbres

ejbmass = 0.0_DP
write(*,*) 'NEW EJECTA'
! !$OMP PARALLEL DO DEFAULT(PRIVATE) IF(inc > INCPAR) &
! !$OMP SHARED(user,domain,crater,surf,ejb,ejtble,mvrld,mvrldsc,nrays,n2,nef,ef,n2f,n1f,nfrays,rayf) &
! !$OMP SHARED(inc,incsq,ejdissq,fradsq,radsq,indarray,cumulative_elchange,rn,continuous)
! !$OMP REDUCTION(+:massray,massej,massrayef)
do j = -inc,inc
do i = -inc,inc

! find distance from crater center
iradsq = i*i + j*j
xpi = crater%xlpx + i
Expand All @@ -267,70 +268,41 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot)
xp = xpi * user%pix
yp = ypi * user%pix

!lradsq = (crater%xl - xp)**2 + (crater%yl - yp)**2
!lrad = sqrt(lradsq)


lradsq = (crater%xl - xp)**2 + (crater%yl - yp)**2
lrad = sqrt(lradsq)

! periodic boundary conditions
call util_periodic(xpi,ypi,user%gridsize)

indarray(1,i,j) = xpi
indarray(2,i,j) = ypi
if ((iradsq > incsq).or.(lradsq < fradsq)) cycle

theta = atan2(j * 1._DP,i * 1._DP) + 2.0_DP * PI
mag = ( ( (abs(cos(nrays * theta / 4.0_DP)))**n2 + &
(abs(sin(nrays * theta / 4.0_DP)))**n2 )**(-1.0_DP/n1) )
magf = rayf * ( ( (abs(cos(nfrays * theta / 4.0_DP)))**n2f + &
(abs(sin(nfrays * theta / 4.0_DP)))**n2f )**(-1.0_DP/n1f))
lradp = continuous * mag
lradf = continuous * magf
lradp = max(lradp, lradf)


! Find angle-corrected landing distance
distance = sqrt((crater%xl - xp)**2 + (crater%yl - yp)**2)

baseline = crater%melev + ((i * crater%xslp) + (j * crater%yslp)) * user%pix
craterslope = atan(baseline / distance)


lrad = distance
lrad0 = lrad
do n = 1,MAXLOOP
!k = max(min(1 + int((lrad - inneredge) / (outeredge - inneredge) * (EJBTABSIZE - 1.0_DP)),ejtble),1)
call ejecta_interpolate(crater,domain,lrad,ejb,ejtble,ebh,vsq=vsq,theta=ejtheta,erad=erad)
erad = exp(erad)
baseline = surf(xpi,ypi)%dem - erad * sin(craterslope)
landslope = atan(baseline / (distance - erad * cos(craterslope)))
if (lrad < lradp) then
call ejecta_interpolate(crater,domain,lrad,ejb,ejtble,ebh,vsq=vsq,theta=ejtheta,melt=melt)
ray_pix = ceiling((lrad - crater%rad) / user%pix)
ebh = ebh * ef(ray_pix)

lradtrue = erad * cos(craterslope) + vsq / (user%gaccel * cos(landslope)) * &
(sin(2 * ejtheta + 2 * craterslope - landslope) - sin(landslope))
lrad = lrad * distance / lradtrue
if (lrad > outeredge) exit
if (abs(lrad - lrad0) < domain%small) exit
lrad0 = lrad
end do
lradsq = lrad**2


if ((lradsq <= ejdissq) .and. (lradsq >= fradsq) .and. (lrad > 0.0_DP)) then
theta = atan2(j * 1._DP,i * 1._DP) + 2.0_DP * PI
mag = ( ( (abs(cos(nrays * theta / 4.0_DP)))**n2 + &
(abs(sin(nrays * theta / 4.0_DP)))**n2 )**(-1.0_DP/n1) )
magf = rayf * ( ( (abs(cos(nfrays * theta / 4.0_DP)))**n2f + &
(abs(sin(nfrays * theta / 4.0_DP)))**n2f )**(-1.0_DP/n1f))
lradp = continuous * mag
lradf = continuous * magf
lradp = max(lradp, lradf)


if (lrad < lradp) then
call ejecta_interpolate(crater,domain,lrad,ejb,ejtble,ebh,vsq=vsq,theta=ejtheta,melt=melt)
ray_pix = ceiling((lrad - crater%rad) / user%pix)
ebh = ebh * ef(ray_pix)

if (user%doregotrack .and. ebh>1.0e-8) then
call regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,lrad,ebh,rm)
end if

else
ebh = 0._DP
if (user%doregotrack .and. ebh>1.0e-8) then
call regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,lrad,ebh,rm)
end if

cumulative_elchange(i,j) = cumulative_elchange(i,j) + ebh
ejbmass = ejbmass + ebh
else
ebh = 0._DP
end if

cumulative_elchange(i,j) = cumulative_elchange(i,j) + ebh
ejbmass = ejbmass + ebh

end do
end do
Expand All @@ -339,7 +311,11 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot)
deallocate(ef)

! Do mass conservation by adjusting ejecta thickness
cumulative_elchange = cumulative_elchange * (-deltaMtot)/ ejbmass
write(*,*) 'ejbmass = ',ejbmass
fmasscons = (-deltaMtot)/ ejbmass
write(*,*) 'fmasscons = ',fmasscons

cumulative_elchange = cumulative_elchange * fmasscons

! Create box for soften calculation (will be no bigger than the grid itself)
if (2 * inc + 1 < user%gridsize) then
Expand Down

0 comments on commit 0c6095b

Please sign in to comment.