diff --git a/src/regolith/regolith_subpixel_streamtube.f90 b/src/regolith/regolith_subpixel_streamtube.f90 index 9b873c0b..5b2458ef 100644 --- a/src/regolith/regolith_subpixel_streamtube.f90 +++ b/src/regolith/regolith_subpixel_streamtube.f90 @@ -50,7 +50,7 @@ ! Notes : ! !********************************************************************************************************************************** -subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer,vmare,totseb) +subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,vseg,newlayer,rm) use module_globals use module_regolith, EXCEPT_THIS_ONE => regolith_subpixel_streamtube implicit none @@ -58,9 +58,9 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer ! Arguments type(usertype),intent(in) :: user type(surftype),intent(inout) :: surfi - real(DP),intent(in) :: deltar,ri,rip1,eradi + real(DP),intent(in) :: deltar,ri,rip1,eradi,vseg type(regodatatype),intent(inout) :: newlayer - real(DP),intent(out) :: vmare,totseb + real(DP),intent(in) :: rm ! Traversing a linked list real(DP),parameter :: a = 0.936457 @@ -68,24 +68,20 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer type(regolisttype),pointer :: current real(DP) :: z,zmax,zstart,zend,rlefti,rleftf,rrighti,rrightf,rc,vsgly,x - ! * Mixing - real(DP) :: zmix - ! The depth that a stream tube dips zmax = rip1/4.0_DP current => surfi%regolayer z = surfi%regolayer%regodata%thickness - zmix = z zstart = 0.0_DP zend = z - vmare = 0._DP - totseb = 0._DP ! Two cases: subpixel is inside the first layer, and its volume is simply the landing ejecta blanket. if (zend>=zmax) then - vmare = newlayer%thickness * user%pix**2 * surfi%regolayer%regodata%comp - totseb = newlayer%thickness * user%pix**2 + + newlayer%thickness = vseg + newlayer%comp = vseg * surfi%regolayer%regodata%comp + else ! The subpixel stream tube may dip deeper than the first layer. And we use the line segments approximation, but ! the layer now intersects with both sides of a stream tube, so two intersection points between the layer and the @@ -112,6 +108,8 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer rrighti = eradi rrightf = eradi rc = rip1 * sqrt(3.0) / 4.0 + newlayer%thickness = 0._DP + newlayer%comp = 0._DP do @@ -126,8 +124,8 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer - tan(b/eradi * rlefti) - tan(b/eradi * rrighti)) - abs(b/eradi * rleftf + b/eradi * rrightf & - b/eradi * rlefti - b/eradi * rrighti)) - vmare = vmare + vsgly * current%regodata%comp - totseb = totseb + vsgly + newlayer%comp = newlayer%comp + vsgly * current%regodata%comp + newlayer%thickness = newlayer%thickness + vsgly current => current%next z = z + current%regodata%thickness @@ -138,14 +136,16 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer ! final part of a stream tube else vsgly = 0.25 * PI * deltar**2 * a**2 * eradi / b * (abs(tan(b) - b)) - vmare = vmare + (vsgly - totseb) * current%regodata%comp - totseb = vsgly + newlayer%comp = newlayer%comp + (vsgly - newlayer%thickness) * current%regodata%comp + newlayer%thickness = vsgly exit end if + end do + end if - call regolith_streamtube_head(user,surfi,deltar,vmare,totseb) + call regolith_streamtube_head(user,surfi,deltar,newlayer,eradi,rm) return end subroutine regolith_subpixel_streamtube