first pass at allocatable arrays (does not compile)

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\

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

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

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)

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)

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(:)

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

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

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)
 

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

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

src/regolith/regolith_traverse_pop.f90

@@ -36,33 +36,33 @@ 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
 
    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

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 

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

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
+
+

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