diff --git a/src/Makefile.am b/src/Makefile.am index 1bdb5932..5cdc3cb6 100755 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -7,11 +7,11 @@ OPTREPORT = -qopt-report=5 IPRODUCTION = -g -traceback -no-wrap-margin -assume byterecl -O3 -qopt-prefetch=0 -sox $(PAR) $(SIMDVEC) $(HEAPARR) #IDEBUG = -O0 -g -traceback -debug all -nogen-interfaces -assume byterecl -m64 -heap-arrays -FR -no-pie -no-ftz -fpe-all=0 -mp1 -fp-model strict -fpe0 -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margin -init=snan,arrays IDEBUG = -O0 -g -traceback -debug all -nogen-interfaces -assume byterecl -m64 -heap-arrays -FR -no-pie -no-ftz -fpe-all=0 -mp1 -fp-model strict -fpe0 -align all -pad -ip -prec-sqrt -assume protect-parens -CB -no-wrap-margin -init=snan,arrays -AM_FCFLAGS = $(IPRODUCTION) +#AM_FCFLAGS = $(IPRODUCTION) #ifort debug flags #gfortran optimized flags -#AM_FCFLAGS = -O3 -fopenmp -ffree-form -g -fbounds-check -fbacktrace +AM_FCFLAGS = -O3 -fopenmp -ffree-form -g -fbounds-check -fbacktrace #gfortran debug flags #AM_FCFLAGS = -O0 -g -fopenmp -fbounds-check -Wall -Warray-bounds -Warray-temporaries -Wimplicit-interface -ffree-form -fsanitize-address-use-after-scope -fstack-check -fsanitize=bounds-strict -fsanitize=undefined -fsanitize=signed-integer-overflow -fsanitize=object-size -fstack-protector-all @@ -36,8 +36,8 @@ util/util_rootbracketed.f90\ util/util_diffusion_solver.f90\ util/util_area_intersection.f90\ util/util_poisson.f90\ -util/util_pop.f90\ -util/util_push.f90\ +util/util_pop_array.f90\ +util/util_push_array.f90\ util/util_traverse_pop.f90\ util/util_destroy_list.f90\ util/util_init_list.f90\ diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index 7af10cce..9e30e06d 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -91,7 +91,7 @@ module module_globals ! type(regolisttype), pointer :: porolayer => null() ! Pointer to the top of the porosity layer stack !type(regolisttype), pointer :: regolayer => null() ! Pointer to the top of the regolith layer stack !type(regolisttype), pointer :: porolayer => null() ! Pointer to the top of the porosity layer stack - type(regodatatype), dimension(:), allocatable: regolayer + type(regodatatype), dimension(:), allocatable :: regolayer end type surftype ! Derived data type for crater information diff --git a/src/init/init_regolith_stack.f90 b/src/init/init_regolith_stack.f90 index b9fd331e..2a4c35b9 100644 --- a/src/init/init_regolith_stack.f90 +++ b/src/init/init_regolith_stack.f90 @@ -51,7 +51,7 @@ subroutine init_regolith_stack(user,surf) call util_init_list(surf(xp,yp)%regolayer,initstat) if (initstat) then - call util_push(surf(xp,yp)%regolayer,bedrock) + call util_push_array(surf(xp,yp)%regolayer,bedrock) else write(*,*) 'init_regolith_stack: Initialization of regolayer failed.' end if diff --git a/src/io/io_read_porotrack.f90 b/src/io/io_read_porotrack.f90 index 28fb8ffc..357c01e0 100644 --- a/src/io/io_read_porotrack.f90 +++ b/src/io/io_read_porotrack.f90 @@ -95,7 +95,7 @@ subroutine io_read_porotrack(user,surf) do k= 2, stacks_num(i,j) newsurfi%depth = depthk(k) newsurfi%porosity = porok(k) - call util_push(surf(i,j)%porolayer,newsurfi) + call util_push_array(surf(i,j)%porolayer,newsurfi) end do deallocate(depthk, porok) diff --git a/src/io/io_read_regotrack.f90 b/src/io/io_read_regotrack.f90 index 0fa7964a..f49246e2 100644 --- a/src/io/io_read_regotrack.f90 +++ b/src/io/io_read_regotrack.f90 @@ -111,7 +111,7 @@ subroutine io_read_regotrack(user,surf) do q=1,MAXAGEBINS newsurfi%age(q) = agei(MAXAGEBINS*k-(MAXAGEBINS-q)) end do - call util_push(surf(i,j)%regolayer,newsurfi) + call util_push_array(surf(i,j)%regolayer,newsurfi) end do deallocate(regotopi,compi,melti,agei) diff --git a/src/io/io_write_regotrack.f90 b/src/io/io_write_regotrack.f90 index 69b265ce..e21ad5e8 100644 --- a/src/io/io_write_regotrack.f90 +++ b/src/io/io_write_regotrack.f90 @@ -32,7 +32,8 @@ subroutine io_write_regotrack(user,surf) integer(I4B), parameter :: FREGO = 11 integer(I4B), parameter :: FCOMP = 12 integer(I4B), parameter :: FAGE = 13 - type(regolisttype),pointer :: current => null() + !type(regolisttype),pointer :: current => null() + type(regodatatype),dimension(:),allocatable :: current integer(I4B),dimension(user%gridsize,user%gridsize) :: stacks_num real(DP),dimension(:),allocatable :: meltfrac, thickness, comp real(SP),dimension(:,:),allocatable :: age @@ -53,27 +54,34 @@ subroutine io_write_regotrack(user,surf) stacks_num(:,:) = 0 do j=1,user%gridsize do i=1,user%gridsize - current => surf(i,j)%regolayer - do - if (.not. associated(current)) exit ! We've reached the bottom of the linked list - stacks_num(i,j) = stacks_num(i,j) + 1 - current => current%next - end do + !current => surf(i,j)%regolayer + current = surf(i,j)%regolayer(1) + ! do + ! if (.not. associated(current)) exit ! We've reached the bottom of the linked list + ! stacks_num(i,j) = stacks_num(i,j) + 1 + ! current => current%next + ! end do end do end do ! Second pass to get data and save it do j=1,user%gridsize do i=1,user%gridsize - current => surf(i,j)%regolayer + !current => surf(i,j)%regolayer N = stacks_num(i,j) allocate(meltfrac(N),thickness(N),comp(N),age(MAXAGEBINS,N)) do k=1,N - meltfrac(k) = current%regodata%meltfrac - thickness(k) = current%regodata%thickness - comp(k) = current%regodata%comp - age(:,k) = current%regodata%age(:) - current => current%next + current = surf(i,j)%regolayer(k) + ! meltfrac(k) = current%regodata%meltfrac + ! thickness(k) = current%regodata%thickness + ! comp(k) = current%regodata%comp + ! age(:,k) = current%regodata%age(:) + ! current => current%next + meltfrac(k) = current%meltfrac + thickness(k) = current%thickness + comp(k) = current%comp + age(:,k) = current%age(:) + end do write(FMELT) meltfrac(:) write(FREGO) thickness(:) diff --git a/src/porosity/porosity_form_interior.f90 b/src/porosity/porosity_form_interior.f90 index e2c82086..cf56b7a1 100644 --- a/src/porosity/porosity_form_interior.f90 +++ b/src/porosity/porosity_form_interior.f90 @@ -64,7 +64,7 @@ subroutine porosity_form_interior(user, surfi, crater, lradsq) end if else ! Linked list if there is only one porosity layer - call util_push(surfi%porolayer, nlayer) + call util_push_array(surfi%porolayer, nlayer) end if return diff --git a/src/regolith/regolith_mix.f90 b/src/regolith/regolith_mix.f90 index dbd08e75..d535bab7 100644 --- a/src/regolith/regolith_mix.f90 +++ b/src/regolith/regolith_mix.f90 @@ -53,7 +53,7 @@ subroutine regolith_mix(surfi,mixing_depth) newlayer%comp = newlayer%comp / newlayer%thickness newlayer%meltfrac = newlayer%meltfrac / newlayer%thickness - call util_push(surfi%regolayer, newlayer) + call util_push_array(surfi%regolayer, newlayer) call util_destroy_list(poppedlist_top) return diff --git a/src/regolith/regolith_streamtube.f90 b/src/regolith/regolith_streamtube.f90 index 0ae14b88..ae9a1c15 100644 --- a/src/regolith/regolith_streamtube.f90 +++ b/src/regolith/regolith_streamtube.f90 @@ -316,7 +316,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, end if - call util_push(surf(xpi,ypi)%regolayer,newlayer) + call util_push_array(surf(xpi,ypi)%regolayer,newlayer) deallocate(xints,yints) diff --git a/src/regolith/regolith_superdomain.f90 b/src/regolith/regolith_superdomain.f90 index a3faa467..062b0be4 100644 --- a/src/regolith/regolith_superdomain.f90 +++ b/src/regolith/regolith_superdomain.f90 @@ -61,7 +61,7 @@ subroutine regolith_superdomain(user,crater,domain,regolayer,ejdistribution,xpi, vej = cvpg * sqrt(user%gaccel * crater%grad) * (erad / crater%grad)**(-1.0_DP / user%mu_b) !equation 18 in Richardson 2009 lrad = ( vej **2 ) / user%gaccel !assume ejection angle is 45 degree. call regolith_melt_glass(user,crater,age,age_resolution,ebh,rm,erad,lrad,deltar,newlayer,xmints,melt) - call util_push(regolayer,newlayer) + call util_push_array(regolayer,newlayer) return end subroutine regolith_superdomain diff --git a/src/regolith/regolith_transport.f90 b/src/regolith/regolith_transport.f90 index b06314ac..e25f833e 100644 --- a/src/regolith/regolith_transport.f90 +++ b/src/regolith/regolith_transport.f90 @@ -66,7 +66,7 @@ subroutine regolith_transport(user,surfi,crater,domain,ejb,ejtble,lrad,ebh,newla newlayer%meltfrac = melt - call util_push(surfi%regolayer,newlayer) + call util_push_array(surfi%regolayer,newlayer) return end subroutine regolith_transport diff --git a/src/regolith/regolith_traverse_pop.f90 b/src/regolith/regolith_traverse_pop.f90 index afeb9776..41a587a3 100755 --- a/src/regolith/regolith_traverse_pop.f90 +++ b/src/regolith/regolith_traverse_pop.f90 @@ -36,7 +36,7 @@ subroutine util_traverse_pop(elchange,surfi,poppedlist) ! Get initial layer's info and the ! desired info that we want to modify! !======================================= - depth = surfi%regolayer%regodata%thickness + depth = surfi%regolayer%thickness !we don't need '%regodata' anymore dz = 0._DP z = elchange mixedregodata%comp = 0.0_DP @@ -44,25 +44,25 @@ subroutine util_traverse_pop(elchange,surfi,poppedlist) if (z < 0._DP) then do - if (.not. associated(surfi%regolayer)) then - write(*,*) 'Major error in regolith_traverse_pop!' - exit - end if + ! if (.not. associated(surfi%regolayer)) then + ! write(*,*) 'Major error in regolith_traverse_pop!' + ! exit + ! end if if (abs(z) <= depth) then dz = depth - abs(z) - surfi%regolayer%regodata%thickness = dz - mixedregodata%comp = mixedregodata%comp + dz * surfi%regolayer%regodata%comp - mixedregodata%meltfrac = mixedregodata%meltfrac + dz * surfi%regolayer%regodata%meltfrac + surfi%regolayer%thickness = dz + mixedregodata%comp = mixedregodata%comp + dz * surfi%regolayer%comp + mixedregodata%meltfrac = mixedregodata%meltfrac + dz * surfi%regolayer%%meltfrac mixedregodata%thickness = mixedregodata%thickness + dz exit else z = abs(z) - surfi%regolayer%regodata%thickness - call util_pop(surfi,oldregodata) + call util_pop_array(surfi,oldregodata) mixedregodata%comp = mixedregodata%comp + oldregodata%thickness * oldregodata%comp mixedregodata%meltfrac = mixedregodata%meltfrac + oldregodata%thickness * oldregodata%meltfrac mixedregodata%thickness = mixedregodata%thickness + oldregodata%thickness - depth = surfi%regolayer%regodata%thickness + depth = surfi%regolayer%thickness end if end do diff --git a/src/util/module_util.f90 b/src/util/module_util.f90 index eb0236e8..78720ef6 100644 --- a/src/util/module_util.f90 +++ b/src/util/module_util.f90 @@ -33,29 +33,49 @@ module module_util public save +! interface +! subroutine util_push(regolayer,newregodata) +! use module_globals +! implicit none +! type(regolisttype),pointer :: regolayer +! type(regodatatype),intent(in) :: newregodata +! end subroutine util_push +! end interface + interface - subroutine util_push(regolayer,newregodata) + subroutine util_push_array(regolayer,newregodata) use module_globals implicit none - type(regolisttype),pointer :: regolayer - type(regodatatype),intent(in) :: newregodata - end subroutine util_push + type(regodatatype),dimension(:),allocatable :: regolayer + type(regodatatype),dimension(:),allocatable,intent(in) :: newregodata + end subroutine util_push_array end interface interface - subroutine util_pop(regolayer,oldregodata) + subroutine util_pop_array(regolayer,oldregodata) use module_globals implicit none - type(regolisttype),pointer :: regolayer + type(regodatatype),dimension(:),allocatable :: regolayer type(regodatatype),intent(out) :: oldregodata - end subroutine util_pop + 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_traverse_pop(regolayer,traverse_depth,poppedlist) use module_globals implicit none - type(regolisttype),pointer :: regolayer + !type(regolisttype),pointer :: regolayer + type(regodatatype),dimension(:),allocatable :: regolayer real(DP),intent(in) :: traverse_depth type(regolisttype),pointer :: poppedlist end subroutine diff --git a/src/util/util_pop_array.f90 b/src/util/util_pop_array.f90 new file mode 100644 index 00000000..7d148d1e --- /dev/null +++ b/src/util/util_pop_array.f90 @@ -0,0 +1,43 @@ +!********************************************************************************************************************************** +! +! Unit Name : util_pop +! Unit Type : function +! Project : CTEM +! Language : Fortran 2003 +! +! Description : Popped a old regolith block onto an old surface +! +! +! Input +! Arguments : regolayer :: pointer to the top of the regolith stack +! oldlayer :: old layer to pop off of the top of the stack +! +! Output +! Arguments : +! +! +! Notes : +! +!********************************************************************************************************************************** +subroutine util_pop_array(regolayer,oldregodata) + use module_globals + use module_util, EXCEPT_THIS_ONE => util_pop_array + implicit none + + ! Arguments + type(regodatatype),dimension(:),allocatable :: regolayer + type(regodatatype),intent(out) :: oldregodata + + ! Internal variables + type(regodatatype), dimension(:), allocatable :: newlayer + integer(I4B) :: allocstat + integer(I4B) :: nold + + ! Executable code + + nold = size(regolayer) + + allocate(newlayer(nold-1), stat=allocstat) + newlayer(1:nold-1) = regolayer(1:nold-1) ! could also be 2:nold depending on if top or bottom is popped off + call move_alloc(newlayer, regolayer) +end subroutine util_pop_array \ No newline at end of file diff --git a/src/util/util_push_array.f90 b/src/util/util_push_array.f90 new file mode 100644 index 00000000..c1aea892 --- /dev/null +++ b/src/util/util_push_array.f90 @@ -0,0 +1,46 @@ +!********************************************************************************************************************************** +! +! Unit Name : util_push_array +! Unit Type : subroutine +! Project : CTEM +! Language : Fortran 2003 +! +! Description : Pushes a new regolith block onto an old surface +! +! +! Input +! Arguments : regolayer :: allocatable array +! newregodata :: new layer to push onto the top of the stack +! +! Output +! Arguments : +! +! +! Notes : +! +!********************************************************************************************************************************** +subroutine util_push_array(regolayer,newregodata) + use module_globals + use module_util, EXCEPT_THIS_ONE => util_push_array + implicit none + + ! Arguments + type(regodatatype),dimension(:),allocatable :: regolayer + type(regodatatype),dimension(:),allocatable :: newregodata + + ! Internal variables + type(regodatatype), dimension(:), allocatable :: newlayer + integer(I4B) :: allocstat + integer(I4B) :: nold + + ! Executable code + + nold = size(regolayer) + + allocate(newregodata(nold+1), stat=allocstat) + newlayer(1:nold) = regolayer(1:nold) + newlayer(nold+1) = newregodata(:) ! need to see if this adds to the top or bottom of the layer stack + call move_alloc(newlayer, regolayer) +end subroutine util_push_array + + diff --git a/src/util/util_traverse_pop.f90 b/src/util/util_traverse_pop.f90 index cefa3d95..f4725482 100755 --- a/src/util/util_traverse_pop.f90 +++ b/src/util/util_traverse_pop.f90 @@ -60,7 +60,8 @@ subroutine util_traverse_pop(regolayer,traverse_depth,poppedlist) implicit none ! Arguments - type(regolisttype),pointer :: regolayer + !type(regolisttype),pointer :: regolayer + type(regodatatype),dimension(:),allocatable :: regolayer real(DP),intent(in) :: traverse_depth type(regolisttype),pointer :: poppedlist @@ -95,12 +96,12 @@ subroutine util_traverse_pop(regolayer,traverse_depth,poppedlist) recyratio = dz / regolayer%regodata%thickness regolayer%regodata%age(:) = recyratio * regolayer%regodata%age(:) regolayer%regodata%thickness = dz - call util_push(poppedlist,oldregodata) + call util_push_array(poppedlist,oldregodata) exit else z = z - regolayer%regodata%thickness - call util_pop(regolayer,oldregodata) - call util_push(poppedlist,oldregodata) + call util_pop_array(regolayer,oldregodata) + call util_push_array(poppedlist,oldregodata) depth = regolayer%regodata%thickness end if