From 4395ac47be7ed9715c6e5600114bf40c29db6f22 Mon Sep 17 00:00:00 2001 From: Austin Blevins Date: Thu, 2 Feb 2023 15:45:43 -0500 Subject: [PATCH] melt tracking WIP sent to cluster for debugging (does not compile) --- src/globals/module_globals.f90 | 4 ++ src/regolith/module_regolith.f90 | 15 +++--- src/regolith/regolith_melt_glass.f90 | 12 ++++- src/regolith/regolith_streamtube.f90 | 25 ++++++++-- src/regolith/regolith_streamtube_head.f90 | 25 +++++++++- src/regolith/regolith_streamtube_lineseg.f90 | 49 +++++++++++++++++-- src/regolith/regolith_subpixel_streamtube.f90 | 24 +++++++-- src/regolith/regolith_traverse_streamtube.f90 | 9 ++-- 8 files changed, 138 insertions(+), 25 deletions(-) diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index 34147a7b..9cb934a1 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -72,6 +72,10 @@ module module_globals real(DP) :: porosity ! Porosity: Maximum 1, Minium 0. real(DP) :: damage ! Damage : Maximum 1, Minium 0. real(DP) :: depth ! Absolute location with respect to the initial surface. + real(DP) :: meltvolume + real(DP) :: totvolume + real(DP) :: ejm !ejected melt + real(DP) :: ejmf !ejected melt fraction real(SP),dimension(:),allocatable :: meltdist !its dimension should be the number of quasimc craters end type regodatatype diff --git a/src/regolith/module_regolith.f90 b/src/regolith/module_regolith.f90 index 9c623b4e..36150d9e 100644 --- a/src/regolith/module_regolith.f90 +++ b/src/regolith/module_regolith.f90 @@ -95,13 +95,13 @@ end subroutine regolith_streamtube interface subroutine regolith_traverse_streamtube(user,surfi,deltar,ri,rip1,eradi,erado,newlayer,vmare,totseb,& - age_collector,xmints,xsfints,rsh,depthb) + age_collector,xmints,xsfints,rsh,depthb,mixedregodata) use module_globals implicit none type(usertype),intent(in) :: user type(surftype),intent(inout) :: surfi real(DP),intent(in) :: deltar,ri,rip1,eradi,erado - type(regodatatype),intent(inout) :: newlayer + type(regodatatype),intent(inout) :: newlayer, mixedregodata real(DP),intent(out) :: vmare,totseb real(SP),dimension(:),intent(inout) :: age_collector real(DP),intent(in) :: xmints @@ -111,13 +111,13 @@ end subroutine regolith_traverse_streamtube interface subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer,vmare,totseb,& - age_collector,xmints,xsfints,vol) + age_collector,xmints,xsfints,vol,mixedregodata) use module_globals implicit none type(usertype),intent(in) :: user type(surftype),intent(in) :: surfi real(DP),intent(in) :: deltar,ri,rip1,eradi - type(regodatatype),intent(inout) :: newlayer + type(regodatatype),intent(inout) :: newlayer, mixedregodata real(DP),intent(out) :: vmare,totseb real(SP),dimension(:),intent(inout) :: age_collector real(DP),intent(in) :: xmints @@ -128,13 +128,13 @@ end subroutine regolith_subpixel_streamtube interface subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad,eradi,deltar,newlayer,vmare,& - totseb,age_collector,xmints,xsfints,depthb) + totseb,age_collector,xmints,xsfints,depthb,mixedregodata) use module_globals implicit none type(usertype),intent(in) :: user type(surftype),intent(in) :: surfi real(DP),intent(in) :: thetast,ri,rip1,zmin,zmax,erad,eradi,deltar - type(regodatatype),intent(inout) :: newlayer + type(regodatatype),intent(inout) :: newlayer,mixedregodata real(DP),intent(inout) :: vmare,totseb real(SP),dimension(:),intent(inout) :: age_collector real(DP),intent(in) :: xmints @@ -143,11 +143,12 @@ end subroutine regolith_streamtube_lineseg end interface interface - subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector) + subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector,mixedregodata) use module_globals implicit none type(usertype),intent(in) :: user type(surftype),intent(in) :: surfi + type(regodatatype),intent(inout) :: mixedregodata real(DP),intent(in) :: deltar real(DP),intent(inout) :: totmare,tots real(SP),dimension(:),intent(inout) :: age_collector diff --git a/src/regolith/regolith_melt_glass.f90 b/src/regolith/regolith_melt_glass.f90 index f25f389d..ad3a60ab 100644 --- a/src/regolith/regolith_melt_glass.f90 +++ b/src/regolith/regolith_melt_glass.f90 @@ -130,6 +130,10 @@ subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad newlayer%meltfrac = 0.0_DP ! default value: no melt (zero) newlayer%thickness = ebh ! default value: stream tube's volume = paraboloid shell's volume newlayer%comp = 0.0_DP + newlayer%meltvolume = 0.0_DP + newlayer%totvolume = 0.0_DP + newlayer%ejm = 0.0_DP + newlayer%ejmf = 0.0_DP rints = sqrt(rm**2 - (crater%imp/2.0)**2) cosvints = min(max(eradi / (crater%imp + eradi), -1.0_DP), 1.0_DP) sinvints = sqrt(1.0 - cosvints**2) @@ -140,6 +144,8 @@ subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad if (eradi <= rints) then volm1 = vst - volv1 melt = volm1 + newlayer%meltvolume = melt + newlayer%totvolume = melt newlayer%meltfrac = 1.0 xmints = rints else if (eradi > rints) then @@ -164,14 +170,18 @@ subroutine regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,erad xmints = eradi * (1.0 - cos(thetaq)) * sin(thetaq) volm1 = regolith_streamtube_volume_func(eradi,0.0_DP,xmints,deltar) melt = volm1 - volv1 + newlayer%meltvolume = melt + newlayer%totvolume = vst-volv1 newlayer%meltfrac = melt/(vst-volv1) + newlayer%ejm = melt + newlayer%ejmf = newlayer%meltfrac end if allocate(newlayer%meltdist((domain%rcnum))) newlayer%meltdist(:) = 0.0_SP if(domain%currentqmc) then - newlayer%meltdist(domain%nqmc) = newlayer%meltfrac + newlayer%meltdist(domain%nqmc) = newlayer%meltfrac !might change this to melt end if n_age = max(ceiling(age / age_resolution), 1) diff --git a/src/regolith/regolith_streamtube.f90 b/src/regolith/regolith_streamtube.f90 index e5e51af0..dd2d916c 100644 --- a/src/regolith/regolith_streamtube.f90 +++ b/src/regolith/regolith_streamtube.f90 @@ -89,7 +89,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, integer(I4B) :: i,j,k,toti,totj,toty,cnt,xstpi,ystpi real(DP) :: vtot,vseg,ri,rip1,xc,yc,thetast real(DP) :: vst,vbody,rbody,vmare,totmare,totseb,tots - type(regodatatype) :: newlayer + type(regodatatype) :: newlayer, mixedregodata ! Constrain the tangital tube's volume with CTEM result real(DP) :: k1,k2,k3,k4,c1,c2 @@ -222,7 +222,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, ! Purpose 2: Once we have the size information of a stream tube, we can ! calculate the distal melt: the precursor of glass spherules within a ! stream tube. The result is contained in a linked list "newlayer". - call regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints, volm) + call regolith_melt_glass(user,crater,domain,age,age_resolution,ebh,rm,eradc,lrad,deltar,newlayer,xmints,volm) ! if (eradc>rm) then ! write(*,*) 'eradc > rm!' ! write(*,*) ebh, exp(ejb(k)%thick) @@ -254,8 +254,11 @@ 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) + age_collector,xmints,xsfints,vol,mixedregodata) newlayer%age(:) = newlayer%age(:) + age_collector(:) @@ -264,6 +267,9 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, newlayer%thickness = ebh newlayer%comp = min(totmare/tots, 1.0_DP) newlayer%age(:) = newlayer%age(:) * min( (ebh * user%pix**2) / tots, 1.0_DP) + newlayer%meltvolume = newlayer%meltvolume + mixedregodata%meltvolume + newlayer%totvolume = newlayer%totvolume + mixedregodata%totvolume + newlayer%meltfrac = newlayer%meltvolume / newlayer%totvolume else rbody = sqrt(xints(2)**2 + yints(2)**2) @@ -275,6 +281,9 @@ 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) totmare = totmare + vmare @@ -292,8 +301,11 @@ 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) + totseb,age_collector,xmints,xsfints,rsh,depthb,mixedregodata) totmare = totmare + vmare tots = tots + totseb end if @@ -302,12 +314,15 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, xstpi = crater%xlpx + nint(eradc*xl/lrad/user%pix) ystpi = crater%ylpx + nint(eradc*yl/lrad/user%pix) call util_periodic(xstpi,ystpi,user%gridsize) - call regolith_streamtube_head(user,surf(xstpi,ystpi),deltar,totmare,tots,age_collector) + call regolith_streamtube_head(user,surf(xstpi,ystpi),deltar,totmare,tots,age_collector,mixedregodata) newlayer%thickness = ebh newlayer%comp = min(totmare/tots, 1.0_DP) newlayer%age(:) = newlayer%age(:) + age_collector(:) newlayer%age(:) = newlayer%age(:) * min( (ebh * user%pix**2) / tots, 1.0_DP) + newlayer%meltvolume = newlayer%meltvolume + mixedregodata%meltvolume + newlayer%totvolume = newlayer%totvolume + mixedregodata%totvolume + newlayer%meltfrac = newlayer%meltvolume / newlayer%totvolume end if diff --git a/src/regolith/regolith_streamtube_head.f90 b/src/regolith/regolith_streamtube_head.f90 index cb5b6d80..97c8570c 100644 --- a/src/regolith/regolith_streamtube_head.f90 +++ b/src/regolith/regolith_streamtube_head.f90 @@ -18,13 +18,14 @@ ! Notes : The stream tube's head is always attached to the surface. ! !********************************************************************************************************************************** -subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector) +subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector,mixedregodata) use module_globals use module_regolith, EXCEPT_THIS_ONE => regolith_streamtube_head implicit none ! arguemnts type(usertype),intent(in) :: user type(surftype),intent(in) :: surfi + type(regodatatype),intent(inout) :: regodatatype real(DP),intent(in) :: deltar real(DP),intent(inout) :: totmare,tots real(SP),dimension(:),intent(inout) :: age_collector @@ -42,6 +43,7 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector ! melt collector real(DP) :: recyratio + real(DP) :: headmeltvol !current => surfi%regolayer !current = surfi%regolayer @@ -63,6 +65,13 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector 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 + mixedregodata%totvolume = mixedregodata%totvolume + vsgly + 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 else ! head is not intersected with layers. do @@ -74,6 +83,13 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector 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)" + end if !current => current%next N = N - 1 z = z + current(N)%thickness @@ -84,6 +100,13 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector 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 diff --git a/src/regolith/regolith_streamtube_lineseg.f90 b/src/regolith/regolith_streamtube_lineseg.f90 index a66d059a..da1fedf8 100644 --- a/src/regolith/regolith_streamtube_lineseg.f90 +++ b/src/regolith/regolith_streamtube_lineseg.f90 @@ -19,7 +19,7 @@ ! !********************************************************************************************************************************** subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad,eradi,deltar,newlayer,vmare,totseb,& - age_collector,xmints,xsfints,depthb) + age_collector,xmints,xsfints,depthb,mixedregodata) use module_globals use module_regolith, EXCEPT_THIS_ONE => regolith_streamtube_lineseg implicit none @@ -27,7 +27,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad type(usertype),intent(in) :: user type(surftype),intent(in) :: surfi real(DP),intent(in) :: thetast,ri,rip1,zmin,zmax,erad,eradi,deltar - type(regodatatype),intent(inout) :: newlayer + type(regodatatype),intent(inout) :: newlayer,mixedregodata real(DP),intent(inout) :: vmare,totseb real(SP),dimension(:),intent(inout) :: age_collector real(DP),intent(in) :: xmints @@ -43,6 +43,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad ! Melt zone real(DP) :: recyratio, xsh, rst real(DP) :: theta1, theta2, r1, r2, vol + real(DP) :: linmelt ! Shock damaged zone real(DP) :: ebh_recyl @@ -81,12 +82,26 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad recyratio = max(ebh_recyl,0.0_DP) / current(N-1)%thickness age_collector(:) = age_collector(:) + current(N-1)%age(:) * recyratio vol = vol + sum(current(N-1)%age(:)) * recyratio + mixedregodata%totvolume = mixedregodata%totvolume + (vsgly-vsh) + 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 + write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)" + end if else if (ri <= xmints .and. rip1 > xmints) then vseg = regolith_streamtube_volume_func(eradi,xmints,rip1,deltar) vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,ri,rip1) recyratio = max((vseg-vsh),0.0_DP) / (user%pix**2) / current(N-1)%thickness age_collector(:) = age_collector(:) + current(N-1)%age(:) * recyratio vol = vol + sum(current(N-1)%age(:)) * recyratio + mixedregodata%totvolume = mixedregodata%totvolume + (vseg-vsh) + 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 + write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)" + end if end if exit else @@ -113,6 +128,13 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad recyratio = max((vseg-vsh),0.0_DP) / (user%pix**2) / current(N)%thickness age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio vol = vol + sum(current(N)%age(:)) * recyratio + mixedregodata%totvolume = mixedregodata%totvolume + (vseg-vsh) + linmelt = current(N)%meltfrac * vseg * recyratio + mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt + mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume + if (mixedregoda%meltfrac > 1.0_DP) then + write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)" + end if end if ! A segment coming from the side of emerging location of a streamtube rip1 @@ -122,6 +144,13 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad recyratio = max((vseg-vsh),0.0_DP) / (user%pix**2) / current(N)%thickness age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio vol = vol + sum(current(N)%age(:)) * recyratio + mixedregodata%totvolume = mixedregodata%totvolume + (vseg-vsh) + linmelt = current(N)%meltfrac * vseg * recyratio + mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt + mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume + if (mixedregoda%meltfrac > 1.0_DP) then + write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)" + end if end if !current => current%next N = N - 1 @@ -143,6 +172,13 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad recyratio = max((vseg-vsh),0.0_DP) / (user%pix**2) / current(N)%thickness age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio vol = vol + sum(current(N)%age(:)) * recyratio + mixedregodata%totvolume = mixedregodata%totvolume + (vseg-vsh) + linmelt = current(N)%meltfrac * vseg * recyratio + mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt + mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume + if (mixedregoda%meltfrac > 1.0_DP) then + write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)" + end if end if exit @@ -151,7 +187,14 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad ! last part of a stream tube vsgly = regolith_streamtube_volume_func(eradi,ri,rip1,deltar) vmare = vmare + (vsgly - totseb) * current(N)%comp - totseb = vsgly + totseb = vsgly + mixedregodata%totvolume = mixedregodata%totvolume + vsgly + linmelt = surfi(N)%meltfrac * vsgly * recyratio + mixedregodata%meltvolume = mixedregodata%meltvolume + linmelt + mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume + if (mixedregoda%meltfrac > 1.0_DP) then + write(*,*) "ERROR! mixedregodata%meltfrac >1! (LINESEG)" + end if exit end if diff --git a/src/regolith/regolith_subpixel_streamtube.f90 b/src/regolith/regolith_subpixel_streamtube.f90 index 3322ce6e..6346d554 100644 --- a/src/regolith/regolith_subpixel_streamtube.f90 +++ b/src/regolith/regolith_subpixel_streamtube.f90 @@ -51,7 +51,7 @@ ! !********************************************************************************************************************************** subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer,vmare,totseb,& - age_collector,xmints,xsfints,vol) + age_collector,xmints,xsfints,vol,mixedregodata) use module_globals use module_regolith, EXCEPT_THIS_ONE => regolith_subpixel_streamtube implicit none @@ -60,7 +60,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer type(usertype),intent(in) :: user type(surftype),intent(in) :: surfi real(DP),intent(in) :: deltar,ri,rip1,eradi - type(regodatatype),intent(inout) :: newlayer + type(regodatatype),intent(inout) :: newlayer, mixedregodata real(DP),intent(out) :: vmare,totseb real(SP),dimension(:),intent(inout) :: age_collector real(DP),intent(in) :: xmints @@ -74,7 +74,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer real(DP),parameter :: b = 1.12368 !type(regolisttype),pointer :: current type(regodatatype),dimension(:),allocatable :: current - real(DP) :: z,zmax,zstart,zend,rlefti,rleftf,rrighti,rrightf,rc,vsgly,vsgly1,vsgly2,x + real(DP) :: z,zmax,zstart,zend,rlefti,rleftf,rrighti,rrightf,rc,vsgly,vsgly1,vsgly2,x,mvl,mvr integer(I4B) :: N ! Stream tube's distance from the edge of a melt zone @@ -88,6 +88,9 @@ 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 @@ -105,7 +108,10 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer 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) + mixedregodata%totvolume = vseg-vsh recyratio = max(vseg-vsh,0.0 )/ (user%pix**2) / (surfi%regolayer(N)%thickness) + mixedregodata%meltvolume = surfi%regolayer(N)%meltfrac * vseg * recyratio + mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume age_collector(:) = age_collector(:) + surfi%regolayer(N)%age(:) * recyratio vol = vol + sum(age_collector(:)) ! write(*,*) '1',eradi, xmints, xsfints, & @@ -181,6 +187,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer recyratio = max((vsgly2 -vsh),0.0)/ (user%pix**2) / current(N)%thickness age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio vol = vol + sum(current(N)%age(:)) * recyratio + mvr = vsgly2 * current(N)%meltfrac * recyratio end if if (rlefti > xmints) then @@ -188,9 +195,16 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer recyratio = max((vsgly1-vsh),0.0) / (user%pix**2) / current(N)%thickness age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio vol = vol + sum(current(N)%age(:)) * recyratio + mvl = vsgly1 * current(N)%meltfrac * recyratio end if !current => current%next + mixedregodata%meltvolume = mixedregodata%meltvolume + mvl + mvr + mixedregodata%totvolume = mixedregodata%totvolume + vsgly + mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume + if (mixedregoda%meltfrac > 1.0_DP) then + write(*,*) "ERROR! mixedregodata%meltfrac >1! (SUBPIXEL)" + end if N = N - 1 z = z + current(N)%thickness zstart = zend @@ -198,7 +212,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer rlefti = rleftf rrightf = rrighti ! final part of a stream tube - else + else !I think this means it's in the melt zone like Ya-Huei said above.. if so, nothing needed here. vsgly = 0.25 * PI * deltar**2 * a**2 * eradi / b * (abs(tan(b) - b)) vmare = vmare + (vsgly - totseb) * current(N)%comp totseb = vsgly @@ -207,7 +221,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer end do end if - call regolith_streamtube_head(user,surfi,deltar,vmare,totseb,age_collector) + call regolith_streamtube_head(user,surfi,deltar,vmare,totseb,age_collector, mixedregodata) return end subroutine regolith_subpixel_streamtube diff --git a/src/regolith/regolith_traverse_streamtube.f90 b/src/regolith/regolith_traverse_streamtube.f90 index b84843dd..6c8636f5 100644 --- a/src/regolith/regolith_traverse_streamtube.f90 +++ b/src/regolith/regolith_traverse_streamtube.f90 @@ -19,7 +19,7 @@ ! !********************************************************************************************************************************** subroutine regolith_traverse_streamtube(user,surfi,deltar,ri,rip1,eradi,erado,newlayer,vmare,totseb,& - age_collector,xmints,xsfints,depthb) + age_collector,xmints,xsfints,depthb,mixedregodata) use module_globals use module_regolith, EXCEPT_THIS_ONE => regolith_traverse_streamtube implicit none @@ -28,7 +28,7 @@ subroutine regolith_traverse_streamtube(user,surfi,deltar,ri,rip1,eradi,erado,ne type(usertype),intent(in) :: user type(surftype),intent(inout) :: surfi real(DP),intent(in) :: deltar,ri,rip1,eradi,erado - type(regodatatype),intent(inout) :: newlayer + type(regodatatype),intent(inout) :: newlayer, mixedregodata real(DP),intent(out) :: vmare,totseb real(SP),dimension(:),intent(inout) :: age_collector real(DP),intent(in) :: xmints @@ -85,12 +85,15 @@ subroutine regolith_traverse_streamtube(user,surfi,deltar,ri,rip1,eradi,erado,ne vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,ri,rip1) recyratio = max(vseg-vsh,0.0_DP) / (user%pix**2) / surfi%regolayer(N)%thickness age_collector(:) = age_collector(:) + surfi%regolayer(N)%age(:) * recyratio + mixedregodata%totvolume = vseg-vsh + mixedregodata%meltvolume = surfi%regolayer(N)%meltfrac * vseg * recyratio + mixedregodata%meltfrac = mixedregodata%meltvolume / mixedregodata%totvolume end if else call regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad,eradi,deltar,& - newlayer,vmare,totseb,age_collector,xmints,xsfints,depthb) + newlayer,vmare,totseb,age_collector,xmints,xsfints,depthb,mixedregodata) end if