program compiles but gives a lot of errors

examples/global-lunar-bombardment/ctem.in

@@ -2,27 +2,27 @@
 
 
 ! 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
 testxoffset  0.0                             ! x-axis offset of crater center from grid center (m) - Default 0.0
 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

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\

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
 
 

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)

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

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

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

src/regolith/regolith_mix.f90

@@ -28,33 +28,39 @@ 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
    newlayer%comp = newlayer%comp / newlayer%thickness 
    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

src/regolith/regolith_streamtube_head.f90

@@ -38,14 +38,16 @@ 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
 
    !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