diff --git a/src/regolith/regolith_streamtube_head.f90 b/src/regolith/regolith_streamtube_head.f90 index ac95a08c..4c6d57fc 100644 --- a/src/regolith/regolith_streamtube_head.f90 +++ b/src/regolith/regolith_streamtube_head.f90 @@ -40,7 +40,7 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector ! head intersected with underlying layers, about 30% of volume difference. real(DP) :: z,zstart,zend,zmin,zmax real(DP) :: tothead,totmarehead,marehead,vhead,vsgly - integer(I4B) :: N + integer(I4B) :: N,M ! melt collector real(DP) :: recyratio @@ -49,8 +49,8 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector !current => surfi%regolayer !current = surfi%regolayer allocate(current,source=surfi%regolayer) - N = size(current) - z = current(N)%thickness + M = size(current) + z = current(M)%thickness vsgly = vratio * PI * deltar**3 tothead = 0._DP totmarehead = 0._DP @@ -63,16 +63,16 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector if (zend >= zmax) then ! Stream tube's head is inside the 1st layer. tots = tots + vsgly - totmare = totmare + vsgly * current(N)%comp - recyratio = vsgly / (user%pix**2) /current(N)%thickness - age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio - headmeltvol = current(N)%meltfrac * recyratio * vsgly + totmare = totmare + vsgly * current(M)%comp + recyratio = vsgly / (user%pix**2) /current(M)%thickness + age_collector(:) = age_collector(:) + current(M)%age(:) * recyratio + headmeltvol = current(M)%meltfrac * recyratio * vsgly meltinejecta = meltinejecta + headmeltvol - distvol(:) = distvol + (current(N)%meltdist(:)*recyratio*vsgly) + distvol(:) = distvol + (current(M)%meltdist(:)*recyratio*vsgly) totvol = totvol + vsgly else ! head is not intersected with layers. - do + do N=M,2,-1 ! if (.not. associated(current%next)) exit if (zend < zmax) then @@ -86,8 +86,8 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector distvol(:) = distvol + (current(N)%meltdist(:)*recyratio*tothead) totvol = totvol + tothead !current => current%next - N = N - 1 - z = z + current(N)%thickness + !N = N - 1 + z = z + current(N-1)%thickness zstart = zend zend = z else diff --git a/src/regolith/regolith_streamtube_lineseg.f90 b/src/regolith/regolith_streamtube_lineseg.f90 index 77646e4f..24137c2f 100644 --- a/src/regolith/regolith_streamtube_lineseg.f90 +++ b/src/regolith/regolith_streamtube_lineseg.f90 @@ -40,7 +40,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad type(regodatatype),dimension(:),allocatable :: current real(DP) :: z,zstart,zend,rstart,rend,r real(DP) :: vsgly,x,vseg, vsh - integer(I4B) :: N + integer(I4B) :: N,M ! Melt zone real(DP) :: recyratio, xsh, rst @@ -51,8 +51,8 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad !current => surfi%regolayer allocate(current,source=surfi%regolayer) - N = size(current) - z = current(N)%thickness + M = size(current) + z = current(M)%thickness zstart = 0.0_DP zend = z @@ -66,7 +66,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad vol = 0.0_DP - do + do N=N,2,-1 !if (.not. associated(current%next)) exit !it should exit until it hit the very bottom. @@ -103,8 +103,8 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad else !current => current%next - N = N - 1 - z = z + current(N)%thickness + !N = N - 1 + z = z + current(N-1)%thickness zstart = zend zend = z @@ -143,8 +143,8 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad totvol = totvol + vseg end if !current => current%next - N = N - 1 - z = z + current(N)%thickness + !N = N - 1 + z = z + current(N-1)%thickness r = rstart rstart = rend zstart = zend diff --git a/src/regolith/regolith_subpixel_streamtube.f90 b/src/regolith/regolith_subpixel_streamtube.f90 index 46758758..fc42f514 100644 --- a/src/regolith/regolith_subpixel_streamtube.f90 +++ b/src/regolith/regolith_subpixel_streamtube.f90 @@ -77,7 +77,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer !type(regolisttype),pointer :: current type(regodatatype),dimension(:),allocatable :: current real(DP) :: z,zmax,zstart,zend,rlefti,rleftf,rrighti,rrightf,rc,vsgly,vsgly1,vsgly2,x,mvl,mvr - integer(I4B) :: N + integer(I4B) :: N,M ! Stream tube's distance from the edge of a melt zone real(DP) :: zm, recyratio, xmints1, vseg, recyratio2 @@ -90,8 +90,8 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer !current => surfi%regolayer allocate(current,source=surfi%regolayer) - N = size(current) - z = surfi%regolayer(N)%thickness + M = size(current) + z = surfi%regolayer(M)%thickness zstart = 0.0_DP zend = z vmare = 0._DP @@ -107,16 +107,16 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer ! 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(N)%comp + vmare = newlayer%thickness * user%pix**2 * surfi%regolayer(M)%comp totseb = newlayer%thickness * user%pix**2 if (eradi>xmints) then vseg = regolith_streamtube_volume_func(eradi,xmints,eradi,deltar) vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,0.0_DP,eradi) - recyratio = max(vseg-vsh,0.0 )/ (user%pix**2) / (surfi%regolayer(N)%thickness) - meltinejecta = surfi%regolayer(N)%meltfrac * vseg * recyratio - distvol = distvol + (surfi%regolayer(N)%meltdist(:)*vseg*recyratio) + recyratio = max(vseg-vsh,0.0 )/ (user%pix**2) / (surfi%regolayer(M)%thickness) + meltinejecta = surfi%regolayer(M)%meltfrac * vseg * recyratio + distvol = distvol + (surfi%regolayer(M)%meltdist(:)*vseg*recyratio) totvol = vseg - age_collector(:) = age_collector(:) + surfi%regolayer(N)%age(:) * recyratio + age_collector(:) = age_collector(:) + surfi%regolayer(M)%age(:) * recyratio vol = vol + sum(age_collector(:)) ! write(*,*) '1',eradi, xmints, xsfints, & ! vseg/user%pix**2, (vseg-vsh)/user%pix**2, recyratio @@ -149,7 +149,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer rrightf = eradi rc = rip1 * sqrt(3.0) / 4.0 - do + do N=M,2,-1 ! It should hit the bottom layer before it exits, I think. !if (.not. associated(current%next)) exit @@ -207,8 +207,8 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer meltinejecta = meltinejecta + mvl + mvr distvol = distvol + (current(N)%meltdist(:)*((vsgly1*recyratio)+(vsgly2*recyratio2))) totvol = totvol + vsgly1 + vsgly2 - N = N - 1 - z = z + current(N)%thickness + !N = N - 1 + z = z + current(N-1)%thickness zstart = zend zend = z rlefti = rleftf