diff --git a/examples/global-lunar-bombardment/ctem.in b/examples/global-lunar-bombardment/ctem.in index 470aaee4..57c12e29 100644 --- a/examples/global-lunar-bombardment/ctem.in +++ b/examples/global-lunar-bombardment/ctem.in @@ -2,7 +2,7 @@ ! Testing input. These are used to perform non-Monte Carlo tests. -testflag F ! Set to T to create a single crater with user-defined impactor properties +testflag T ! Set to T to create a single crater with user-defined impactor properties testimp 15000.0 ! Diameter of test impactor (m) testvel 18.3e3 ! Velocity of test crater (m/s) testang 45.0 ! Impact angle of test crater (deg) - Default 90.0 @@ -10,19 +10,19 @@ testxoffset 0.0 ! x-axis offset of crater center fr testyoffset 0.0 ! y-axis offset of crater center from grid center (m) - Default 0.0 tallyonly F ! Tally the craters without generating any craters testtally F -quasimc T ! MC run constrained by non-MC 'real' craters given in a list +quasimc F ! MC run constrained by non-MC 'real' craters given in a list realcraterlist craterlist.in ! list of 'real' craters for Quasi-MC runs ! IDL driver in uts -interval 6.280 -numintervals 100 ! Total number of intervals (total time = interval * numintervals) <--when runtype is 'single' +interval 1 +numintervals 1 ! Total number of intervals (total time = interval * numintervals) <--when runtype is 'single' restart F ! Restart a previous run impfile NPFextrap.dat ! Impactor SFD rate file (col 1: Dimp (m), col 2: ! impactors > D (m**(-2) y**(-1)) popupconsole F ! Pop up console window every output interval saveshaded F ! Output shaded relief images -saverego F ! Output regolith map images +saverego T ! Output regolith map images savepres F ! Output simplified console display images (presentation-compatible images) savetruelist T ! Save the true cumulative crater distribution for each interval (large file size) sfdcompare LOLASethCraterCatalogv8gt20-binned.dat ! File name for the SFD comparison file used in the console display @@ -70,7 +70,7 @@ Kd1 0.0001 psi 2.000 fe 4.00 ejecta_truncation 4.0 -doregotrack F -dorays T +doregotrack T +dorays F superdomain F dorealistic F diff --git a/src/Makefile.am b/src/Makefile.am index ada5ce5a..124ad0c9 100755 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -36,6 +36,7 @@ util/util_rootbracketed.f90\ util/util_diffusion_solver.f90\ util/util_area_intersection.f90\ util/util_poisson.f90\ +util/util_pop.f90\ util/util_pop_array.f90\ util/util_push_array.f90\ util/util_traverse_pop_array.f90\ diff --git a/src/crater/crater_form_interior.f90 b/src/crater/crater_form_interior.f90 index 955055e6..aacafde6 100644 --- a/src/crater/crater_form_interior.f90 +++ b/src/crater/crater_form_interior.f90 @@ -36,8 +36,9 @@ subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative ,newele real(DP) :: cform,newdem,elchange,r integer(I4B) :: layer - ! A list for popped data - type(regolisttype),pointer :: poppedlist + ! An array for popped data + !type(regolisttype),pointer :: poppedlist + type(regodatatype),dimension(:),allocatable :: poppedarray ! Empirical crater shape parameters from Fassett et al. (2014) @@ -79,8 +80,9 @@ subroutine crater_form_interior(user,surfi,crater,x_relative, y_relative ,newele surfi%ejcov = max(surfi%ejcov + elchange,0.0_DP) if (user%doregotrack) then - call util_traverse_pop(surfi%regolayer,abs(elchange),poppedlist) - call util_destroy_list(poppedlist) + call util_traverse_pop_array(surfi%regolayer,abs(elchange),poppedarray) + !call util_destroy_list(poppedlist) + deallocate(poppedarray) end if diff --git a/src/init/init_regolith_stack.f90 b/src/init/init_regolith_stack.f90 index 2a4c35b9..14c50519 100644 --- a/src/init/init_regolith_stack.f90 +++ b/src/init/init_regolith_stack.f90 @@ -48,7 +48,8 @@ subroutine init_regolith_stack(user,surf) do yp = 1, user%gridsize do xp = 1, user%gridsize - call util_init_list(surf(xp,yp)%regolayer,initstat) + !call util_init_list(surf(xp,yp)%regolayer,initstat) + call util_init_array(surf(xp,yp)%regolayer,initstat) if (initstat) then call util_push_array(surf(xp,yp)%regolayer,bedrock) diff --git a/src/io/io_read_regotrack.f90 b/src/io/io_read_regotrack.f90 index f49246e2..5327a2aa 100644 --- a/src/io/io_read_regotrack.f90 +++ b/src/io/io_read_regotrack.f90 @@ -84,7 +84,8 @@ subroutine io_read_regotrack(user,surf) do j=1,user%gridsize do i=1,user%gridsize - call util_init_list(surf(i,j)%regolayer,initstat) + !call util_init_list(surf(i,j)%regolayer,initstat) + call util_init_array(surf(i,j)%regolayer,initstat) allocate(regotopi(stacks_num(i,j))) allocate(compi(stacks_num(i,j))) diff --git a/src/main/CTEM.f90 b/src/main/CTEM.f90 index da0ea979..4d4d0519 100644 --- a/src/main/CTEM.f90 +++ b/src/main/CTEM.f90 @@ -160,19 +160,20 @@ program CTEM if (user%doregotrack) then do yp = 1, user%gridsize do xp = 1, user%gridsize - call util_destroy_list(surf(xp,yp)%regolayer) + !call util_destroy_list(surf(xp,yp)%regolayer) + deallocate(surf(xp,yp)%regolayer) end do end do end if ! If doporosity is true, then destroy the linked list for porosity -if (user%doporosity) then - do yp = 1, user%gridsize - do xp = 1, user%gridsize - call util_destroy_list(surf(xp,yp)%porolayer) - end do - end do -end if +! if (user%doporosity) then +! do yp = 1, user%gridsize +! do xp = 1, user%gridsize +! call util_destroy_list(surf(xp,yp)%porolayer) +! end do +! end do +! end if !$ t2 = omp_get_wtime() diff --git a/src/porosity/porosity_form_interior.f90 b/src/porosity/porosity_form_interior.f90 index cf56b7a1..98638849 100644 --- a/src/porosity/porosity_form_interior.f90 +++ b/src/porosity/porosity_form_interior.f90 @@ -52,20 +52,20 @@ subroutine porosity_form_interior(user, surfi, crater, lradsq) depthht = ht - (parabht * lradsq) !The depth at a given pixel. newtrn = crater%melev - depthht !The elevation of the bottom of the transient crater at a given pixel. - nlayer%porosity = 0.2; !The porosity of the first value. This value is an assumed value. - nlayer%depth = newtrn; + ! nlayer%porosity = 0.2; !The porosity of the first value. This value is an assumed value. + ! nlayer%depth = newtrn; - ! Currently consider only the first layer. - if (surfi%porolayer%regodata%porosity == nlayer%porosity) then - ! just update if the newly calculated layer is deeper than the stored one. - if (surfi%porolayer%regodata%depth > nlayer%depth) then - surfi%porolayer%regodata%depth = nlayer%depth - surfi%porolayer%regodata%porosity = nlayer%porosity - end if - else - ! Linked list if there is only one porosity layer - call util_push_array(surfi%porolayer, nlayer) - end if + ! ! ! Currently consider only the first layer. + ! ! if (surfi%porolayer%regodata%porosity == nlayer%porosity) then + ! ! ! just update if the newly calculated layer is deeper than the stored one. + ! ! if (surfi%porolayer%regodata%depth > nlayer%depth) then + ! ! surfi%porolayer%regodata%depth = nlayer%depth + ! ! surfi%porolayer%regodata%porosity = nlayer%porosity + ! ! end if + ! ! else + ! ! ! Linked list if there is only one porosity layer + ! ! call util_push_array(surfi%porolayer, nlayer) + ! ! end if return end subroutine porosity_form_interior diff --git a/src/regolith/regolith_mix.f90 b/src/regolith/regolith_mix.f90 index d535bab7..b2f60cbd 100644 --- a/src/regolith/regolith_mix.f90 +++ b/src/regolith/regolith_mix.f90 @@ -28,25 +28,30 @@ subroutine regolith_mix(surfi,mixing_depth) ! Internal variables type(regodatatype) :: newlayer - type(regolisttype),pointer :: poppedlist,poppedlist_top + !type(regolisttype),pointer :: poppedlist,poppedlist_top + type(regodatatype),dimension(:),allocatable :: poppedarray, poppedarray_top + integer(I4B) :: N !=============================================== ! Add up all layers' info until a desired depth !=============================================== - call util_traverse_pop(surfi%regolayer,mixing_depth,poppedlist_top) + call util_traverse_pop_array(surfi%regolayer,mixing_depth,poppedarray_top) newlayer%thickness = 0.0_DP newlayer%comp = 0.0_DP newlayer%meltfrac = 0.0_DP newlayer%age(:) = 0.0_DP - poppedlist => poppedlist_top - do while(associated(poppedlist%next)) - newlayer%thickness = newlayer%thickness + poppedlist%regodata%thickness - newlayer%comp = newlayer%comp + poppedlist%regodata%thickness * poppedlist%regodata%comp - newlayer%meltfrac = newlayer%meltfrac + poppedlist%regodata%thickness * poppedlist%regodata%meltfrac - newlayer%age(:) = newlayer%age(:) + poppedlist%regodata%age(:) - poppedlist => poppedlist%next + !poppedlist => poppedlist_top + !do while(associated(poppedlist%next)) + N = size(poppedarray) + do + newlayer%thickness = newlayer%thickness + poppedarray(N)%thickness + newlayer%comp = newlayer%comp + poppedarray(N)%thickness * poppedarray(N)%comp + newlayer%meltfrac = newlayer%meltfrac + poppedarray(N)%thickness * poppedarray(N)%meltfrac + newlayer%age(:) = newlayer%age(:) + poppedarray(N)%age(:) + !poppedlist => poppedlist%next + N = N - 1 end do ! Get average values of composition and melt fraction @@ -54,7 +59,8 @@ subroutine regolith_mix(surfi,mixing_depth) newlayer%meltfrac = newlayer%meltfrac / newlayer%thickness call util_push_array(surfi%regolayer, newlayer) - call util_destroy_list(poppedlist_top) + !call util_destroy_list(poppedlist_top) + deallocate(poppedarray_top) return end subroutine regolith_mix diff --git a/src/regolith/regolith_streamtube_head.f90 b/src/regolith/regolith_streamtube_head.f90 index cb7568fc..cb5b6d80 100644 --- a/src/regolith/regolith_streamtube_head.f90 +++ b/src/regolith/regolith_streamtube_head.f90 @@ -38,6 +38,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 ! melt collector real(DP) :: recyratio @@ -45,7 +46,8 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector !current => surfi%regolayer !current = surfi%regolayer allocate(current,source=surfi%regolayer) - z = current%thickness + N = size(current) + z = current(N)%thickness vsgly = vratio * PI * deltar**3 tothead = 0._DP totmarehead = 0._DP @@ -58,9 +60,9 @@ 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%comp - recyratio = vsgly / (user%pix**2) /current%thickness - age_collector(:) = age_collector(:) + current%age(:) * recyratio + totmare = totmare + vsgly * current(N)%comp + recyratio = vsgly / (user%pix**2) /current(N)%thickness + age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio else ! head is not intersected with layers. do @@ -69,18 +71,19 @@ subroutine regolith_streamtube_head(user,surfi,deltar,totmare,tots,age_collector if (zend < zmax) then vhead = regolith_circle_sector_func(deltar,zstart,zend) tothead = tothead + vhead * vratio - totmarehead = totmarehead + vhead * vratio * current%comp - recyratio = vhead * vratio / (user%pix**2) / current%thickness - age_collector(:) = age_collector(:) + current%age(:) * recyratio + totmarehead = totmarehead + vhead * vratio * current(N)%comp + recyratio = vhead * vratio / (user%pix**2) / current(N)%thickness + age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio !current => current%next - z = z + current%thickness + N = N - 1 + z = z + current(N)%thickness zstart = zend zend = z else - totmarehead = totmarehead + (vsgly-tothead) * current%comp + totmarehead = totmarehead + (vsgly-tothead) * current(N)%comp tothead = vsgly - recyratio = (vsgly - tothead) / (user%pix**2) / current%thickness - age_collector(:) = age_collector(:) + current%age(:) * recyratio + recyratio = (vsgly - tothead) / (user%pix**2) / current(N)%thickness + age_collector(:) = age_collector(:) + current(N)%age(:) * recyratio exit end if end do diff --git a/src/regolith/regolith_streamtube_lineseg.f90 b/src/regolith/regolith_streamtube_lineseg.f90 index e859f468..a66d059a 100644 --- a/src/regolith/regolith_streamtube_lineseg.f90 +++ b/src/regolith/regolith_streamtube_lineseg.f90 @@ -34,9 +34,11 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad real(DP),intent(in) :: xsfints, depthb ! internal variables - type(regolisttype),pointer :: current + !type(regolisttype),pointer :: current + type(regodatatype),dimension(:),allocatable :: current real(DP) :: z,zstart,zend,rstart,rend,r real(DP) :: vsgly,x,vseg, vsh + integer(I4B) :: N ! Melt zone real(DP) :: recyratio, xsh, rst @@ -44,8 +46,10 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad ! Shock damaged zone real(DP) :: ebh_recyl - current => surfi%regolayer - z = current%regodata%thickness + !current => surfi%regolayer + allocate(current,source=surfi%regolayer) + N = size(current) + z = current(N)%thickness zstart = 0.0_DP zend = z @@ -61,36 +65,37 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad do - if (.not. associated(current%next)) exit !it should exit until it hit the very bottom. + !if (.not. associated(current%next)) exit !it should exit until it hit the very bottom. if (zend <= zmin) then - if (zmax <= zend + current%next%regodata%thickness) then + if (zmax <= zend + current(N-1)%thickness) then vsgly = newlayer%thickness * user%pix**2 - vmare = vsgly * current%next%regodata%comp + vmare = vsgly * current(N-1)%comp totseb = vsgly if (ri > xmints) then vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,ri,rip1) ebh_recyl = (1.0 - newlayer%meltfrac) * newlayer%thickness - vsh / user%pix**2 - recyratio = max(ebh_recyl,0.0_DP) / current%next%regodata%thickness - age_collector(:) = age_collector(:) + current%next%regodata%age(:) * recyratio - vol = vol + sum(current%next%regodata%age(:)) * recyratio + 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 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%next%regodata%thickness - age_collector(:) = age_collector(:) + current%next%regodata%age(:) * recyratio - vol = vol + sum(current%next%regodata%age(:)) * recyratio + 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 end if exit else - current => current%next - z = z + current%regodata%thickness - zstart = zend - zend = z + !current => current%next + N = N - 1 + z = z + current(N)%thickness + zstart = zend + zend = z end if @@ -98,28 +103,29 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad rend = regolith_quadratic_func(zend,erad,ri,rip1,rstart) vsgly = regolith_streamtube_volume_func(eradi,rstart,rend,deltar) - vmare = vmare + vsgly * current%regodata%comp + vmare = vmare + vsgly * current(N)%comp totseb = totseb + vsgly ! A segment coming from the side of impact site, ri if (thetast>=0_DP .and. rend > xmints) then vseg = regolith_streamtube_volume_func(eradi,max(xmints,rstart),rend,deltar) vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,rstart,rend) - recyratio = max((vseg-vsh),0.0_DP) / (user%pix**2) / current%regodata%thickness - age_collector(:) = age_collector(:) + current%regodata%age(:) * recyratio - vol = vol + sum(current%regodata%age(:)) * recyratio + 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 end if ! A segment coming from the side of emerging location of a streamtube rip1 if (thetast<0._DP .and. rstart > xmints) then vseg = regolith_streamtube_volume_func(eradi,max(xmints,rend),rstart,deltar) vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,rend,rstart) - recyratio = max((vseg-vsh),0.0_DP) / (user%pix**2) / current%regodata%thickness - age_collector(:) = age_collector(:) + current%regodata%age(:) * recyratio - vol = vol + sum(current%regodata%age(:)) * recyratio + 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 end if - current => current%next - z = z + current%regodata%thickness + !current => current%next + N = N - 1 + z = z + current(N)%thickness r = rstart rstart = rend zstart = zend @@ -128,15 +134,15 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad else if (zend >= zmax .and. zstart <= zmin) then vsgly = newlayer%thickness * user%pix**2 - vmare = vsgly * current%regodata%comp + vmare = vsgly * current(N)%comp totseb = vsgly if (rip1 > xmints) then vseg = regolith_streamtube_volume_func(eradi,max(xmints,ri),rip1,deltar) vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,ri,rip1) - recyratio = max((vseg-vsh),0.0_DP) / (user%pix**2) / current%regodata%thickness - age_collector(:) = age_collector(:) + current%regodata%age(:) * recyratio - vol = vol + sum(current%regodata%age(:)) * recyratio + 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 end if exit @@ -144,7 +150,7 @@ subroutine regolith_streamtube_lineseg(user,surfi,thetast,ri,rip1,zmin,zmax,erad else if (zend >= zmax .and. zstart > zmin) then ! last part of a stream tube vsgly = regolith_streamtube_volume_func(eradi,ri,rip1,deltar) - vmare = vmare + (vsgly - totseb) * current%regodata%comp + vmare = vmare + (vsgly - totseb) * current(N)%comp totseb = vsgly exit end if diff --git a/src/regolith/regolith_subcrater_mix.f90 b/src/regolith/regolith_subcrater_mix.f90 index c5f5e66f..5065b96e 100644 --- a/src/regolith/regolith_subcrater_mix.f90 +++ b/src/regolith/regolith_subcrater_mix.f90 @@ -39,6 +39,8 @@ subroutine regolith_subcrater_mix(user,surf,domain,nflux,finterval,p) ! Test age real(SP) :: age_mean + integer(I4B) :: N + ! Find the deepest depth for 100% true saturation klo = 1 ds = p(1,1) @@ -61,7 +63,9 @@ subroutine regolith_subcrater_mix(user,surf,domain,nflux,finterval,p) dd = p(1,domain%smallest_impactor_index) end if - if (surf(i,j)%regolayer%regodata%thickness < dd) then + N = size(surf(i,j)%regolayer) + + if (surf(i,j)%regolayer(N)%thickness < dd) then call regolith_mix(surf(i,j),dd) end if diff --git a/src/regolith/regolith_subpixel_streamtube.f90 b/src/regolith/regolith_subpixel_streamtube.f90 index 57f30744..3322ce6e 100644 --- a/src/regolith/regolith_subpixel_streamtube.f90 +++ b/src/regolith/regolith_subpixel_streamtube.f90 @@ -66,13 +66,16 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer real(DP),intent(in) :: xmints real(DP),intent(in) :: xsfints real(DP),intent(inout) :: vol - + + ! Internal variables ! Traversing a linked list real(DP),parameter :: a = 0.936457 real(DP),parameter :: b = 1.12368 - type(regolisttype),pointer :: current + !type(regolisttype),pointer :: current + type(regodatatype),dimension(:),allocatable :: current real(DP) :: z,zmax,zstart,zend,rlefti,rleftf,rrighti,rrightf,rc,vsgly,vsgly1,vsgly2,x + integer(I4B) :: N ! Stream tube's distance from the edge of a melt zone real(DP) :: zm, recyratio, xmints1, vseg @@ -83,8 +86,10 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer ! The depth that a stream tube dips zmax = rip1/4.0_DP - current => surfi%regolayer - z = surfi%regolayer%regodata%thickness + !current => surfi%regolayer + allocate(current,source=surfi%regolayer) + N = size(current) + z = surfi%regolayer(N)%thickness zstart = 0.0_DP zend = z vmare = 0._DP @@ -95,13 +100,13 @@ 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%regodata%comp + vmare = newlayer%thickness * user%pix**2 * surfi%regolayer(N)%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%regodata%thickness) - age_collector(:) = age_collector(:) + surfi%regolayer%regodata%age(:) * recyratio + recyratio = max(vseg-vsh,0.0 )/ (user%pix**2) / (surfi%regolayer(N)%thickness) + age_collector(:) = age_collector(:) + surfi%regolayer(N)%age(:) * recyratio vol = vol + sum(age_collector(:)) ! write(*,*) '1',eradi, xmints, xsfints, & ! vseg/user%pix**2, (vseg-vsh)/user%pix**2, recyratio @@ -137,7 +142,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer do ! It should hit the bottom layer before it exits, I think. - if (.not. associated(current%next)) exit + !if (.not. associated(current%next)) exit if (zendmax(xmints,sqrt(rm**2 - z**2))) then if (rrighti > xmints) then vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,rrighti,rrightf) - recyratio = max((vsgly2 -vsh),0.0)/ (user%pix**2) / current%regodata%thickness - age_collector(:) = age_collector(:) + current%regodata%age(:) * recyratio - vol = vol + sum(current%regodata%age(:)) * recyratio + 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 end if if (rlefti > xmints) then vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,rlefti,rleftf) - recyratio = max((vsgly1-vsh),0.0) / (user%pix**2) / current%regodata%thickness - age_collector(:) = age_collector(:) + current%regodata%age(:) * recyratio - vol = vol + sum(current%regodata%age(:)) * recyratio + 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 end if - current => current%next - z = z + current%regodata%thickness + !current => current%next + N = N - 1 + z = z + current(N)%thickness zstart = zend zend = z rlefti = rleftf @@ -194,7 +200,7 @@ subroutine regolith_subpixel_streamtube(user,surfi,deltar,ri,rip1,eradi,newlayer ! final part of a stream tube else vsgly = 0.25 * PI * deltar**2 * a**2 * eradi / b * (abs(tan(b) - b)) - vmare = vmare + (vsgly - totseb) * current%regodata%comp + vmare = vmare + (vsgly - totseb) * current(N)%comp totseb = vsgly exit end if diff --git a/src/regolith/regolith_traverse_streamtube.f90 b/src/regolith/regolith_traverse_streamtube.f90 index 4951a8b3..a91684c1 100644 --- a/src/regolith/regolith_traverse_streamtube.f90 +++ b/src/regolith/regolith_traverse_streamtube.f90 @@ -34,9 +34,12 @@ subroutine regolith_traverse_streamtube(user,surfi,deltar,ri,rip1,eradi,erado,ne real(DP),intent(in) :: xmints real(DP),intent(in) :: xsfints, depthb + ! Internal variables + ! Traversing a linked list real(DP) :: zri,zrip1,cosi,coso,rzmax real(DP) :: erad,z,zmin,zmax,thetast,vseg + integer(I4B) :: N real(DP),parameter :: a = 0.936457 real(DP),parameter :: b = 1.12368 @@ -63,19 +66,21 @@ subroutine regolith_traverse_streamtube(user,surfi,deltar,ri,rip1,eradi,erado,ne zmax = max(zmax,erad/4.0) end if - z = surfi%regolayer%regodata%thickness + N = size(surfi%regolayer) + + z = surfi%regolayer(N)%thickness vmare = 0._DP totseb = 0._DP if (z>=zmax) then - vmare = newlayer%thickness * user%pix**2 * surfi%regolayer%regodata%comp + vmare = newlayer%thickness * user%pix**2 * surfi%regolayer(N)%comp totseb = newlayer%thickness * user%pix**2 if (rip1 > xmints .and. ri < xmints) then vseg = regolith_streamtube_volume_func(eradi,max(xmints,ri),rip1,deltar) vsh = regolith_shock_damage(eradi,deltar,xmints,xsfints,ri,rip1) - recyratio = max(vseg-vsh,0.0_DP) / (user%pix**2) / surfi%regolayer%regodata%thickness - age_collector(:) = age_collector(:) + surfi%regolayer%regodata%age(:) * recyratio + recyratio = max(vseg-vsh,0.0_DP) / (user%pix**2) / surfi%regolayer(N)%thickness + age_collector(:) = age_collector(:) + surfi%regolayer(N)%age(:) * recyratio end if else diff --git a/src/util/module_util.f90 b/src/util/module_util.f90 index 75b9731b..a00a84c5 100644 --- a/src/util/module_util.f90 +++ b/src/util/module_util.f90 @@ -61,14 +61,14 @@ end subroutine util_pop_array end interface -! interface -! subroutine util_pop(regolayer,oldregodata) -! use module_globals -! implicit none -! type(regolisttype),pointer :: regolayer -! type(regodatatype),intent(out) :: oldregodata -! end subroutine util_pop -! end interface +interface + subroutine util_pop(regolayer,oldregodata) + use module_globals + implicit none + type(regolisttype),pointer :: regolayer + type(regodatatype),intent(out) :: oldregodata + end subroutine util_pop +end interface ! interface ! subroutine util_traverse_pop(regolayer,traverse_depth,poppedlist)