diff --git a/src/regolith/regolith_streamtube.f90 b/src/regolith/regolith_streamtube.f90 index dd2d916c..69ffe66c 100644 --- a/src/regolith/regolith_streamtube.f90 +++ b/src/regolith/regolith_streamtube.f90 @@ -117,6 +117,10 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, ! Executalbe code + mixedregodata%totvolume = 0 + mixedregodata%meltvolume = 0 + mixedregodata%meltfrac = 0 + ! ****** Interpolate radial distance, erad, for a given pixel ******* ! outeredge = crater%frad + domain%ejbres * (EJBTABSIZE - 0.5_DP) ! inneredge = crater%frad + 0.5_DP * domain%ejbres @@ -254,9 +258,6 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, vseg = regolith_streamtube_volume_func(eradi,0.0_DP,eradi,deltar) newlayer%thickness = vseg/(user%pix**2) call util_periodic(xstpi,ystpi,user%gridsize) - mixedregodata%meltfrac = surf(xstpi,ystpi)%regolayer%meltfrac - mixedregodata%meltvolume = surf(xstpi,ystpi)%regolayer%meltvolume - mixedregodata%totvolume = surf(xstpi,ystpi)%regolayer%totvolume call regolith_subpixel_streamtube(user,surf(xstpi,ystpi),deltar,ri,rip1,eradi,newlayer,vmare,totseb,& age_collector,xmints,xsfints,vol,mixedregodata) @@ -281,11 +282,8 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, vseg = regolith_streamtube_volume_func(eradi,rbody,eradi,deltar) newlayer%thickness = vseg/(user%pix**2) call util_periodic(xstpi,ystpi,user%gridsize) - mixedregodata%meltfrac = surf(xstpi,ystpi)%regolayer%meltfrac - mixedregodata%meltvolume = surf(xstpi,ystpi)%regolayer%meltvolume - mixedregodata%totvolume = surf(xstpi,ystpi)%regolayer%totvolume call regolith_traverse_streamtube(user,surf(xstpi,ystpi),deltar,rbody,eradi,eradi,erado,newlayer,vmare,& - totseb,age_collector,xmints,xsfints,rsh,depthb) + totseb,age_collector,xmints,xsfints,rsh,depthb,mixedregodata) totmare = totmare + vmare tots = tots + totseb end if @@ -301,9 +299,6 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, vseg = regolith_streamtube_volume_func(eradi,ri,rip1,deltar) newlayer%thickness = vseg/(user%pix**2) call util_periodic(xstpi,ystpi,user%gridsize) - mixedregodata%meltfrac = surf(xstpi,ystpi)%regolayer%meltfrac - mixedregodata%meltvolume = surf(xstpi,ystpi)%regolayer%meltvolume - mixedregodata%totvolume = surf(xstpi,ystpi)%regolayer%totvolume call regolith_traverse_streamtube(user,surf(xstpi,ystpi),deltar,ri,rip1,eradi,erado,newlayer,vmare,& totseb,age_collector,xmints,xsfints,rsh,depthb,mixedregodata) totmare = totmare + vmare diff --git a/src/regolith/regolith_streamtube_head.f90 b/src/regolith/regolith_streamtube_head.f90 index 97c8570c..a13cb3bd 100644 --- a/src/regolith/regolith_streamtube_head.f90 +++ b/src/regolith/regolith_streamtube_head.f90 @@ -25,7 +25,7 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector ! arguemnts type(usertype),intent(in) :: user type(surftype),intent(in) :: surfi - type(regodatatype),intent(inout) :: regodatatype + type(regodatatype),intent(inout) :: mixedregodata real(DP),intent(in) :: deltar real(DP),intent(inout) :: totmare,tots real(SP),dimension(:),intent(inout) :: age_collector @@ -69,50 +69,51 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector mixedregodata%totvolume = mixedregodata%totvolume + vsgly mixedregodata%meltvolume = mixedregodata%meltvolume + headmeltvol mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume - if (mixedregoda%meltfrac > 1.0_DP) then + if (mixedregodata%meltfrac > 1.0_DP) then write(*,*) "ERROR! mixedregodata%meltfrac >1! (HEAD)" end if else ! head is not intersected with layers. - do - ! if (.not. associated(current%next)) exit - - if (zend < zmax) then - vhead = regolith_circle_sector_func(deltar,zstart,zend) - tothead = tothead + vhead * vratio - totmarehead = totmarehead + vhead * vratio * current(N)%comp - recyratio = vhead * vratio / (user%pix**2) / current(N)%thickness - age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio - mixedregodata%totvolume = mixedregodata%totvolume + tothead - headmeltvol = current(N)%meltfrac * recyratio * tothead - mixedregodata%meltvolume = mixedregodata%meltvolume + headmeltvol - mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume - if (mixedregoda%meltfrac > 1.0_DP) then - write(*,*) "ERROR! mixedregodata%meltfrac >1! (HEAD)" + do + ! if (.not. associated(current%next)) exit + + if (zend < zmax) then + vhead = regolith_circle_sector_func(deltar,zstart,zend) + tothead = tothead + vhead * vratio + totmarehead = totmarehead + vhead * vratio * current(N)%comp + recyratio = vhead * vratio / (user%pix**2) / current(N)%thickness + age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio + mixedregodata%totvolume = mixedregodata%totvolume + tothead + headmeltvol = current(N)%meltfrac * recyratio * tothead + mixedregodata%meltvolume = mixedregodata%meltvolume + headmeltvol + mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume + if (mixedregodata%meltfrac > 1.0_DP) then + write(*,*) "ERROR! mixedregodata%meltfrac >1! (HEAD)" + end if + !current => current%next + N = N - 1 + z = z + current(N)%thickness + zstart = zend + zend = z + else + totmarehead = totmarehead + (vsgly-tothead) * current(N)%comp + tothead = vsgly + recyratio = (vsgly - tothead) / (user%pix**2) / current(N)%thickness + age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio + mixedregodata%totvolume = mixedregodata%totvolume + tothead + headmeltvol = current(N)%meltfrac * recyratio * tothead + mixedregodata%meltvolume = mixedregodata%meltvolume + headmeltvol + mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume + if (mixedregodata%meltfrac > 1.0_DP) then + write(*,*) "ERROR! mixedregodata%meltfrac >1! (HEAD)" + end if + exit end if - !current => current%next - N = N - 1 - z = z + current(N)%thickness - zstart = zend - zend = z - else - totmarehead = totmarehead + (vsgly-tothead) * current(N)%comp - tothead = vsgly - recyratio = (vsgly - tothead) / (user%pix**2) / current(N)%thickness - age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio - mixedregodata%totvolume = mixedregodata%totvolume + tothead - headmeltvol = current(N)%meltfrac * recyratio * tothead - mixedregodata%meltvolume = mixedregodata%meltvolume + headmeltvol - mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume - if (mixedregoda%meltfrac > 1.0_DP) then - write(*,*) "ERROR! mixedregodata%meltfrac >1! (HEAD)" - end if - exit - end if - end do + end do + - tots = tots + tothead - totmare = totmare + totmarehead + tots = tots + tothead + totmare = totmare + totmarehead end if diff --git a/src/regolith/regolith_streamtube_lineseg.f90 b/src/regolith/regolith_streamtube_lineseg.f90 index da1fedf8..fe99112c 100644 --- a/src/regolith/regolith_streamtube_lineseg.f90 +++ b/src/regolith/regolith_streamtube_lineseg.f90 @@ -86,7 +86,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad linmelt = current(N-1)%meltfrac * vsgly * recyratio mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume - if (mixedregoda%meltfrac > 1.0_DP) then + if (mixedregodata%meltfrac > 1.0_DP) then write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)" end if else if (ri <= xmints .and. rip1 > xmints) then @@ -99,7 +99,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad linmelt = current(N-1)%meltfrac * vseg * recyratio mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume - if (mixedregoda%meltfrac > 1.0_DP) then + if (mixedregodata%meltfrac > 1.0_DP) then write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)" end if end if @@ -132,7 +132,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad linmelt = current(N)%meltfrac * vseg * recyratio mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume - if (mixedregoda%meltfrac > 1.0_DP) then + if (mixedregodata%meltfrac > 1.0_DP) then write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)" end if end if @@ -148,7 +148,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad linmelt = current(N)%meltfrac * vseg * recyratio mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume - if (mixedregoda%meltfrac > 1.0_DP) then + if (mixedregodata%meltfrac > 1.0_DP) then write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)" end if end if @@ -176,7 +176,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad linmelt = current(N)%meltfrac * vseg * recyratio mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume - if (mixedregoda%meltfrac > 1.0_DP) then + if (mixedregodata%meltfrac > 1.0_DP) then write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)" end if end if @@ -189,10 +189,10 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad vmare = vmare + (vsgly - totseb) * current(N)%comp totseb = vsgly mixedregodata%totvolume = mixedregodata%totvolume + vsgly - linmelt = surfi(N)%meltfrac * vsgly * recyratio + linmelt = current(N)%meltfrac * vsgly * recyratio mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume - if (mixedregoda%meltfrac > 1.0_DP) then + if (mixedregodata%meltfrac > 1.0_DP) then write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)" end if exit diff --git a/src/regolith/regolith_subpixel_streamtube.f90 b/src/regolith/regolith_subpixel_streamtube.f90 index 6346d554..aba3df9f 100644 --- a/src/regolith/regolith_subpixel_streamtube.f90 +++ b/src/regolith/regolith_subpixel_streamtube.f90 @@ -88,9 +88,6 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer !current => surfi%regolayer allocate(current,source=surfi%regolayer) - mixedregodata%totvolume = 0 - mixedregodata%meltvolume = 0 - mixedregodata%meltfrac = 0 N = size(current) z = surfi%regolayer(N)%thickness zstart = 0.0_DP @@ -202,7 +199,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer mixedregodata%meltvolume = mixedregodata%meltvolume + mvl + mvr mixedregodata%totvolume = mixedregodata%totvolume + vsgly mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume - if (mixedregoda%meltfrac > 1.0_DP) then + if (mixedregodata%meltfrac > 1.0_DP) then write(*,*) "ERROR! mixedregodata%meltfrac >1! (SUBPIXEL)" end if N = N - 1