made do loops explicit in an attempt to get debugging to work

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 

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 

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