Skip to content

Commit

Permalink
Fixed bugs arising from not converting the logspace values in the eje…
Browse files Browse the repository at this point in the history
…cta table to real space.
  • Loading branch information
daminton committed Feb 11, 2024
1 parent 56054d9 commit b0e64a2
Show file tree
Hide file tree
Showing 3 changed files with 8 additions and 3 deletions.
4 changes: 2 additions & 2 deletions src/ejecta/ejecta_emplace.f90
Original file line number Diff line number Diff line change
Expand Up @@ -282,10 +282,10 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ

! Find out where in the table this new velocity corresponds to
ind = 1
call util_search(ejb%vesq,ind,ejtble,vsq,klo)
call util_search(ejb%vesq,ind,ejtble,log(vsq),klo)
klo = min(max(klo,1),ejtble-1)
! Interpolate on the table to find the flat plane equivalent landing distance for this velocity
frac = (vsq - ejb(klo)%vesq) / (ejb(klo+1)%vesq - ejb(klo)%vesq)
frac = (vsq - exp(ejb(klo)%vesq)) / (exp(ejb(klo+1)%vesq) - exp(ejb(klo)%vesq))
distance = exp(ejb(klo)%lrad) + frac * (exp(ejb(klo+1)%lrad) - exp(ejb(klo)%lrad))
end do

Expand Down
2 changes: 1 addition & 1 deletion src/regolith/regolith_streamtube.f90
Original file line number Diff line number Diff line change
Expand Up @@ -154,7 +154,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,
end if

if (eradc<=0.0_DP) then
write(*,*) k,ebh,crater%ejdis,lrad,exp(ejb(k-1)%lrad),exp(ejb(k)%lrad),eradc,ejb(k-1)%erad,ejb(k)%erad
write(*,*) k,ebh,crater%ejdis,lrad,exp(ejb(k-1)%lrad),exp(ejb(k)%lrad),eradc,exp(ejb(k-1)%erad),exp(ejb(k)%erad)
stop
end if

Expand Down
5 changes: 5 additions & 0 deletions src/regolith/regolith_transport.f90
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,11 @@ subroutine regolith_transport(user,surfi,crater,domain,ejb,ejtble,lrad,ebh,newla
frac = (loglrad - logtablerad) / logdelta
melt = ejb(k)%meltfrac - ((ejb(k)%meltfrac - ejb(k+1)%meltfrac) * frac)
end if
if (melt < LOGVSMALL) then
melt = 0.0_DP
else
melt = exp(melt)
end if

newlayer%meltvolume = melt * newlayer%totvolume

Expand Down

0 comments on commit b0e64a2

Please sign in to comment.