diff --git a/src/regolith/regolith_streamtube_head.f90 b/src/regolith/regolith_streamtube_head.f90 index 38735c85..9cef521b 100644 --- a/src/regolith/regolith_streamtube_head.f90 +++ b/src/regolith/regolith_streamtube_head.f90 @@ -18,7 +18,7 @@ ! Notes : The stream tube's head is always attached to the surface. ! !********************************************************************************************************************************** -subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots) +subroutine regolith_streamtube_head(user,surfi,deltar,newlayer,eradi,rm) !subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,turnover) use module_globals use module_regolith, EXCEPT_THIS_ONE => regolith_streamtube_head @@ -27,7 +27,8 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots) type(usertype),intent(in) :: user type(surftype),intent(in) :: surfi real(DP),intent(in) :: deltar - real(DP),intent(inout) :: totmare,tots + type(regodatatype),intent(inout) :: newlayer + real(DP),intent(in) :: eradi,rm ! internal variables type(regolisttype),pointer :: current @@ -36,52 +37,49 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots) ! from the analytical function that we use here to approximate the stream tube's ! head intersected with underlying layers, about 30% of volume difference. real(DP) :: z,zstart,zend,zmin,zmax - real(DP) :: tothead,totmarehead,marehead,vhead,vsgly - - ! * Mixing - real(DP) :: zmix + real(DP) :: headtot,headcomp,headmelt,vhead,vsgly current => surfi%regolayer z = current%regodata%thickness - zmix = z vsgly = vratio * PI * deltar**3 - tothead = 0._DP - totmarehead = 0._DP - vhead = 0._DP - marehead = 0.0_DP zstart = 0._DP zend = z zmin = zstart zmax = 2.0 * deltar if (zend >= zmax) then ! Stream tube's head is inside the 1st layer. - tots = tots + vsgly - totmare = totmare + vsgly * current%regodata%comp - !write(*,*) '0',zstart,zend,zmin,zmax,current%comp,totmare/2500.0,tots/2500.0 + + newlayer%thickness = newlayer%thickness + vsgly + newlayer%comp = newlayer%comp + vsgly * current%regodata%comp + else ! head is not intersected with layers. + headtot = 0._DP + headcomp = 0._DP + headmelt = 0._DP + vhead = 0._DP + 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%regodata%comp - !write(*,*) '1',zstart,zend,current%comp,vhead*vratio/2500.0,totmarehead/2500.0,tothead/2500.0 + headtot = headtot + vhead * vratio + headcomp = headcomp + vhead * vratio * current%regodata%comp current => current%next z = z + current%regodata%thickness zstart = zend zend = z else - totmarehead = totmarehead + (vsgly-tothead) * current%regodata%comp - tothead = vsgly + headcomp = headcomp + (vsgly-headtot) * current%regodata%comp + headtot = vsgly !write(*,*) '2',zstart,zend,current%comp,(vsgly-tothead)/2500.0,totmarehead/2500.0,tothead/2500.0 exit end if end do - tots = tots + tothead - totmare = totmare + totmarehead + newlayer%thickness = newlayer%thickness + headtot + newlayer%comp = newlayer%comp + headcomp end if