Non-functional version of the nest stitching code from end of Grenoble visit in June 2011

src/NEST/Makefile

@@ -11,26 +11,36 @@ OBJECTS = \
 ../module_definitions.o \
 ../module_random.o \
 ../module_threads.o \
+../check_delaunay.o \
 ../compute_normals.o \
+../compute_positive_volume.o \
 ../create_surfaces.o \
 ../define_cloud.o \
 ../define_ov.o \
 ../define_surface.o \
+../DoRuRe.o \
+../erosion.o \
 ../heap.o \
 ../initialize_temperature.o \
+../move_surface.o \
 ../octreev_shrink_xyz.o \
 ../read_controlling_parameters.o \
 ../read_input_file.o \
+read_restart.o \
+../remove_point.o \
 ../scanfile.o \
 stitch_nest.o \
-../toolbox.o
+../toolbox.o \
+../write_global_output.f90
 
 #OCTREE = ../OCTREE/libOctree$(BITS).a
 
 LIBS = \
- -L../OCTREE -lOctree$(BITS) \
+ -L../CASCADE -lcascade$(BITS) \
  -L../NN -lnn_f$(BITS) \
  -L../NN -lnn_c$(BITS) \
+ -L../OCTREE -lOctree$(BITS) \
+ -L../RESAMPLE -lresample$(BITS) \
  $(MPI_LIB)
 
 # make the program

src/NEST/read_restart.f90

@@ -21,14 +21,10 @@ type (ziso) :: zi
 type (nest_info) :: nest
 
 ! Internal variables
-double precision :: current_time
-integer :: err;i,is,k
-
-osolve%noctree=params%noctreemax
-allocate (osolve%octree(osolve%noctree),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%octree in read_restart$')
+integer :: err,i,is,j,k,zisosize
 
 open (9,file=trim(params%restartfile),status='old',form='unformatted')
-read (9) osolve%octree(3), &
+read (9) osolve%noctree,   &
          osolve%nnode,     &
          osolve%nleaves,   &
          osolve%nface,     &
@@ -40,10 +36,10 @@ read (9) osolve%octree(3), &
 allocate (osolve%x(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%x in read_restart$')
 allocate (osolve%y(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%y in read_restart$')
 allocate (osolve%z(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%z in read_restart$')
-allocate (osolve%unode(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%unode in read_restart$')
-allocate (osolve%vnode(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%vnode in read_restart$')
-allocate (osolve%wnode(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%wnode in read_restart$')
-allocate (osolve%wnodeiso(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%wnodeiso in read_restart$')
+allocate (osolve%u(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%u in read_restart$')
+allocate (osolve%v(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%v in read_restart$')
+allocate (osolve%w(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%w in read_restart$')
+allocate (osolve%wiso(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%wiso in read_restart$')
 allocate (osolve%lsf(osolve%nnode,osolve%nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%lsf in read_restart$')
 allocate (osolve%temp(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%temp in read_restart$')
 allocate (osolve%strain(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%strain in read_restart$')
@@ -58,16 +54,29 @@ allocate (osolve%eviscosity(osolve%nleaves),stat=err) ; if (err.ne.0) call stop_
 allocate (osolve%is_plastic(osolve%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%is_plastic in read_restart$')
 allocate (osolve%dilatr(osolve%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%dilatr in read_restart$')
 allocate (osolve%matnum(osolve%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%matnum in read_restart$')
+allocate (osolve%octree(osolve%noctree),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%octree in read_restart$')
 allocate (osolve%iface(9,osolve%nface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc osolve%iface in read_restart$')
-allocate (vo%node(osolve%nnode),stat=threadinfo%err) ; if (err.ne.0) call stop_run ('Error alloc vo%node in read_restart$')
-allocate (vo%leaf(osolve%nnode),stat=threadinfo%err) ; if (err.ne.0) call stop_run ('Error alloc vo%leaf in read_restart$')
-allocate (vo%ftr(osolve%nnode),stat=threadinfo%err) ; if (err.ne.0) call stop_run ('Error alloc vo%ftr in read_restart$')
-allocate (vo%rtf(osolve%nnode),stat=threadinfo%err) ; if (err.ne.0) call stop_run ('Error alloc vo%rtf in read_restart$')
-allocate (vo%influid(osolve%nnode),stat=threadinfo%err) ; if (err.ne.0) call stop_run ('Error alloc vo%influid in read_restart$')
-allocate (vo%face(osolve%nface),stat=threadinfo%err) ; if (err.ne.0) call stop_run ('Error alloc vo%face in read_restart$')
-
-allocate (ov%temporary_nodal_pressure(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%temp_nodal_pressure in define_ov$')
-allocate (ov%whole_leaf_in_fluid(ov%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%whole_leaf_in_fluid in define_ov$')
+allocate (vo%node(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vo%node in read_restart$')
+allocate (vo%leaf(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vo%leaf in read_restart$')
+allocate (vo%ftr(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vo%ftr in read_restart$')
+allocate (vo%rtf(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vo%rtf in read_restart$')
+allocate (vo%influid(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vo%influid in read_restart$')
+allocate (vo%face(osolve%nface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vo%face in read_restart$')
+allocate (ov%temporary_nodal_pressure(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%temp_nodal_pressure in define_ov$')
+allocate (ov%whole_leaf_in_fluid(osolve%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%whole_leaf_in_fluid in define_ov$')
+allocate (cl%x(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%x in read_restart$')
+allocate (cl%y(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%y in read_restart$')
+allocate (cl%z(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%z in read_restart$')
+allocate (cl%x0(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%x0 in read_restart$')
+allocate (cl%y0(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%x0 in read_restart$')
+allocate (cl%z0(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%x0 in read_restart$')
+allocate (cl%strain(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%strain in read_restart$')
+allocate (cl%lsf0(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%lsf0 in read_restart$')
+allocate (cl%temp(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%temp in read_restart$')
+allocate (cl%press(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%press in read_restart$')
+allocate (cl%tag(cl%np), stat=err) ; if (err.ne.0) call stop_run ('Error alloc cl%tag in read_restart$')
+
+write (*,*) 'Read osolve nodes'
 
 ! Read info for octree solve nodes (x,y,z,u,v,w,lsf,temp)
 read (9) (osolve%x(i),                         &
@@ -77,7 +86,7 @@ read (9) (osolve%x(i),                         &
           osolve%v(i),                         & 
           osolve%w(i),                         &
           osolve%wiso(i),                      &
-          (osolve%lsf(i,j),j=1,osolve%nlsf),   &
+          osolve%lsf(i,1:osolve%nlsf),         &
           osolve%temp(i),                      &
           ov%temporary_nodal_pressure(i),      &
           osolve%strain(i),                    &
@@ -87,8 +96,10 @@ read (9) (osolve%x(i),                         &
           osolve%kfixt(i),                     &
           i=1,osolve%nnode)
 
+write (*,*) 'Read icon' 
+
 ! Read icon array
-read (9) ((osolve%icon(k,i),k=1,8),   &
+read (9) (osolve%icon(1:8,i),         &
            osolve%pressure(i),        &
            osolve%spressure(i),       &
            osolve%crit(i),            &
@@ -100,12 +111,15 @@ read (9) ((osolve%icon(k,i),k=1,8),   &
            ov%whole_leaf_in_fluid(i), &
            i=1,osolve%nleaves)
 
+write (*,*) 'Read octree' 
 ! Read octree information
 read (9) (osolve%octree(i),i=1,osolve%octree(3))
 
+write (*,*) 'Read bad faces' 
 ! Read bad face information
-read (9) ((osolve%iface(k,i),k=1,9),i=1,osolve%nface)
+read (9) (osolve%iface(1:9,i),i=1,osolve%nface)
 
+write (*,*) 'Read void' 
 ! Read void information
 read (9) (vo%node(i),    &
           vo%leaf(i),    &
@@ -114,11 +128,14 @@ read (9) (vo%node(i),    &
           vo%influid(i), &
           i=1,osolve%nnode)
 
+write (*,*) 'Read bad faces 2' 
 ! Read bad faces information
 read (9) (vo%face(i),i=1,osolve%nface)
 
+write (*,*) 'Read surfaces' 
 ! Read surface information (r,s,x,y,z,xn,yn,zn)
 do is=1,osolve%nlsf
+write (*,*) 'Read surface ',is 
   read (9) surface(is)%nsurface,        &
            surface(is)%activation_time, &
            surface(is)%nt
@@ -148,13 +165,14 @@ do is=1,osolve%nlsf
             surface(is)%v(i),  &
             surface(is)%w(i),  &
             i=1,surface(is)%nsurface)
-  read (9) ((surface(is)%icon(k,i),k=1,3),i=1,surface(is)%nt)
+  read (9) (surface(is)%icon(1:3,i),i=1,surface(is)%nt)
 
   ! correct surface for vertical exageration
   surface(is)%z=surface(is)%z/params%vex
   surface(is)%zn=surface(is)%zn/params%vex
 enddo
 
+write (*,*) 'Read cloud' 
 ! Read cloud information
 read (9) (cl%x(i),             &
           cl%y(i),             &
@@ -173,17 +191,18 @@ read (9) (cl%x(i),             &
 cl%z=cl%z/params%vex
 cl%z0=cl%z0/params%vex
 
+write (*,*) 'Read ziso' 
 ! Read isostasy basal displacement array
-read (9) 2**params%levelmax_oct
-allocate (zi%zisodisp(2**params%levelmax_oct+1,2**params%levelmax_oct+1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc zi%zisodisp in read_restart$')
+read (9) zisosize
+allocate (zi%zisodisp(zisosize+1,zisosize+1), stat=err) ; if (err.ne.0) call stop_run ('Error alloc zi%zisodisp in read_restart$')
 if (params%isobc) then
-  read (9) ((zi%zisodisp(i,j)+surface(osolve%nlsf)%sp01, &
-             j=1,2**params%levelmax_oct+1),i=1,2**params%levelmax_oct+1)
+  read (9) ((zi%zisodisp(i,j),j=1,zisosize+1),i=1,zisosize+1)
+  zi%zisodisp=zi%zisodisp-surface(osolve%nlsf)%sp01
 else
-  read (9) ((0.d0+surface(osolve%nlsf)%sp01,j=1,2**params%levelmax_oct+1), &
-             i=1,2**params%levelmax_oct+1)
+  zi%zisodisp=0.d0
 endif
 
+write (*,*) 'Read nest' 
 ! Read nested model info
 if (params%nest) then
   read (9) nest%sselemx,nest%sselemy,nest%sselemz,nest%xminls,nest%yminls,&
@@ -198,3 +217,5 @@ else
 endif
 
 close (9)
+
+end subroutine read_restart

src/NEST/stitch_nest.f90

@@ -17,27 +17,33 @@ include 'mpif.h'
 ! Variable type declaration
 character (len=4) :: cistep,cnnest
 
-double precision :: current_time,inc,step,total
+double precision :: current_time,inc,step,total,xss,yss,zss
+double precision :: xmin,xmax,ymin,ymax,zmin,zmax
 
-integer :: err,i,ierr,iproc,istep,material0,nnest,nproc
+integer :: err,i,ierr,iproc,istep,iter,j,material0,nnest,nproc
 
 type (box),dimension(:),allocatable :: boxes
 type (cloud),dimension(:),allocatable :: ss_cl
-type (cloud) :: ls_cl
+type (cloud) :: ls0_cl,ls_cl
 type (cross_section),dimension(:),allocatable :: sections
 type (face),dimension(6) :: cube_faces
-type (material),dimension(:),allocatable::mat
-type (nest_info) nest
+type (material),dimension(:),allocatable :: mat
+type (nest_info),dimension(:),allocatable :: nest
+type (octreelsf) :: olsf
+type (octreesolve) :: ls_osolve
 type (octreev),dimension(:),allocatable :: ss_ov
-type (octreev) ls_ov
+type (octreev) ls0_ov,ls_ov
 type (parameters) :: params
-type (sheet),dimension(:),allocatable :: ls_surf,ls_surf0,ss_surf,ss_surf0
+type (sheet),dimension(:,:),allocatable :: ss_surf,ss_surf0
+type (sheet),dimension(:),allocatable :: ls0_surf,ls0_surf0,ls_surf,ls_surf0
 type (thread) :: threadinfo
-type (ziso) zi
+type (void) :: ls_vo
+type (ziso) :: zi
 
 !-------------------------------------------------------------------------------
 
-write (*,*) 'Program started'
+!write (*,*) 'Program started'
+call show_time (total,step,inc,0,'Program started$')
 
 call mpi_init(ierr)
 call mpi_comm_size (mpi_comm_world,nproc,ierr)
@@ -45,9 +51,6 @@ call mpi_comm_rank (mpi_comm_world,iproc,ierr)
 
 ! Init some variables
 current_time=0.d0
-inc=0.d0
-step=0.d0
-total=0.d0
 
 ! Read/echo current time step and # of nests from command line
 call getarg (1,cistep)
@@ -57,37 +60,178 @@ write (*,*) 'istep = ',trim(cistep)
 write (*,*) 'nnest = ',trim(cnnest)
 
 ! Convert character values to integers
-write (istep,'(i)') cistep
-write (nnest,'(i)') cnnest
+read (cistep,'(i4)') istep
+read (cnnest,'(i4)') nnest
 
 ! Read the coarse model input file
 params%infile=trim('input_ls.txt')
 call read_controlling_parameters (params,threadinfo,'main')
-allocate (ls_surf  (params%ns),stat=err); if (err.ne.0) call stop_run ('Error alloc ls_surf in main$')
-allocate (ls_surf0 (params%ns),stat=err); if (err.ne.0) call stop_run ('Error alloc ls_surf0 in main$')
-call read_input_file (params,threadinfo,material0,mat,ls_surf,boxes,sections,  &
-                      cube_faces,nest)
+allocate (ls0_surf  (params%ns),stat=err); if (err/=0) call stop_run ('Error alloc ls0_surf in main$')
+allocate (ls0_surf0 (params%ns),stat=err); if (err/=0) call stop_run ('Error alloc ls0_surf0 in main$')
+allocate (mat (0:params%nmat),stat=err); if (err/=0) call stop_run ('Error alloc mat in main$')
+allocate (nest (nnest),stat=err); if (err/=0) call stop_run ('Error alloc nest in main$')
+if (params%nboxes.gt.0) then
+  allocate (boxes(params%nboxes),stat=err) ; if (err/=0) call stop_run ('Error alloc boxes arrays$')
+else
+  allocate (boxes(1),stat=err) ! necessary to avoid nil size argument
+endif
+if (params%nsections.gt.0) then
+  allocate (sections(params%nsections),stat=err) ; if (err/=0) call stop_run ('Error alloc cross_section arrays$')
+else
+  allocate (sections(1),stat=err) ! necessary to avoid nil size argument
+end if
+
+call read_input_file (params,threadinfo,material0,mat,ls0_surf,boxes,sections,  &
+                      cube_faces,nest(1))
 
+params%irestart=0
 ! Modify params%restartfile to read previous restart file
 if (istep > 1) then
   call int_to_char (cistep,4,istep-1)
-  params%restartfile='OUT/time_'//cistep//'.bin'
+  ! This may not always be the OUT directory...need to keep it general
+  params%irestart=istep
+  params%restartfile='LSOUT/time_'//cistep//'.bin'
 endif
 
 ! Read the previous coarse model restart file if istep > 1, define surface,
 ! cloud and ov
-call define_surface (params,ls_surf,threadinfo,total,step,inc,current_time)
-call define_cloud (ls_cl,params,zi)
-call define_ov (ls_ov,params,threadinfo)
+call define_surface (params,ls0_surf,threadinfo,total,step,inc,current_time)
+call define_cloud (ls0_cl,params,zi)
+call define_ov (ls0_ov,params,threadinfo)
 
 ! Read the coarse model restart file from the current time step
+call int_to_char (cistep,4,istep)
+! As above, this may now always be in the OUT directory
+params%restartfile='LSOUT/time_'//cistep//'.bin'
+
+! Read the current coarse model restart file, define surface, cloud and ov
+allocate (ls_surf  (params%ns),stat=err); if (err/=0) call stop_run ('Error alloc ls0_surf in main$')
+allocate (ls_surf0 (params%ns),stat=err); if (err/=0) call stop_run ('Error alloc ls0_surf0 in main$')
+call read_restart (params,istep,iter,current_time,ls_osolve,ls_ov,ls_vo,ls_surf,ls_cl,zi,nest(1))
+
+deallocate(mat,boxes,sections)
+
+allocate (ss_surf(params%ns,nnest),stat=err); if (err/=0) call stop_run ('Error alloc ss_surf in main$')
+allocate (ss_surf0(params%ns,nnest),stat=err); if (err/=0) call stop_run ('Error alloc ss_surf in main$')
+allocate (ss_cl(nnest),stat=err); if (err/=0) call stop_run ('Error alloc ss_cl in main$')
+allocate (ss_ov(nnest),stat=err); if (err/=0) call stop_run ('Error alloc ss_ov in main$')
 
 ! Loop through fine model(s) and read their current restart files
-do i = 1,nnest
-  ! Change params%infile
+do i=1,nnest
+  ! Read the coarse model input file
+  params%infile=trim('input_ss.txt')
+  call read_controlling_parameters (params,threadinfo,'main')
+  allocate (mat (0:params%nmat),stat=err); if (err/=0) call stop_run ('Error alloc mat in main$')
+  if (params%nboxes.gt.0) then
+    allocate (boxes(params%nboxes),stat=err) ; if (err/=0) call stop_run ('Error alloc boxes arrays$')
+  else
+    allocate (boxes(1),stat=err) ! necessary to avoid nil size argument
+  endif
+  if (params%nsections.gt.0) then
+    allocate (sections(params%nsections),stat=err) ; if (err/=0) call stop_run ('Error alloc cross_section arrays$')
+  else
+    allocate (sections(1),stat=err) ! necessary to avoid nil size argument
+  endif
+
+  call read_input_file (params,threadinfo,material0,mat,ss_surf(:,i),boxes,sections,  &
+                        cube_faces,nest(i))
+
+  ! Modify params%restartfile to read previous restart file
+  params%restartfile='SSOUT/time_'//cistep//'.bin'
+
+  call define_surface (params,ss_surf(:,i),threadinfo,total,step,inc,current_time)
+  call define_cloud (ss_cl(i),params,zi)
+  call define_ov (ss_ov(i),params,threadinfo)
+  !call read_restart (params,istep,iter,current_time,ls_osolve,ls_ov,ls_vo,ls_surf,ls_cl,zi,nest(1))
+enddo
+
+do i=1,nnest
+  xmin=nest(i)%xminls
+  xmax=nest(i)%xminls+nest(i)%sselemx
+  ymin=nest(i)%yminls
+  ymax=nest(i)%yminls+nest(i)%sselemy
+  zmin=nest(i)%zminls
+  zmax=nest(i)%zminls+nest(i)%sselemz
+  do j=1,ls_osolve%nnode
+    xss=(ls_osolve%x(i)-nest(i)%xminls)/nest(i)%sselemx
+    yss=(ls_osolve%y(i)-nest(i)%yminls)/nest(i)%sselemy
+    zss=(ls_osolve%z(i)-nest(i)%zminls)/nest(i)%sselemz
+    if(ls_osolve%x(i)<=xmax .and. ls_osolve%x(i)>=xmin .and. ls_osolve%y(i)<=ymax .and. ls_osolve%y(i)>=ymin .and. ls_osolve%z(i)<=zmax .and. ls_osolve%z(i)>=zmin) then
+      call octree_interpolate_three (3,ss_ov(i)%octree,ss_ov(i)%noctree,ss_ov(i)%icon,       &
+                                   ss_ov(i)%nleaves,ss_ov(i)%nnode,xss,yss,zss,         &
+                                   ss_ov(i)%unode,ls_osolve%u(i),ss_ov(i)%vnode,ls_osolve%v(i),ss_ov(i)%wnode, &
+                                   ls_osolve%w(i))
+    endif
+  enddo
+enddo
+
+! Create ov arrays for move_surface
+ls_ov%noctree=ls_osolve%noctree
+ls_ov%nnode=ls_osolve%nnode
+ls_ov%nleaves=ls_osolve%nleaves
+
+allocate (ls_ov%octree(ls_ov%noctree),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ls_ov%octree in stitch_nest$')
+allocate (ls_ov%x(ls_ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ls_ov%x in stitch_nest$')
+allocate (ls_ov%y(ls_ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ls_ov%y in stitch_nest$')
+allocate (ls_ov%z(ls_ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ls_ov%z in stitch_nest$')
+allocate (ls_ov%unode(ls_ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ls_ov%unode in stitch_nest$')
+allocate (ls_ov%vnode(ls_ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ls_ov%vnode in stitch_nest$')
+allocate (ls_ov%wnode(ls_ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ls_ov%wnode in stitch_nest$')
+allocate (ls_ov%wnodeiso(ls_ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ls_ov%wnodeiso in stitch_nest$')
+allocate (ls_ov%icon(8,ls_ov%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ls_ov%icon in stitch_nest$')
+
+ls_ov%octree=ls_osolve%octree
+ls_ov%x=ls_osolve%x
+ls_ov%y=ls_osolve%y
+ls_ov%z=ls_osolve%z
+ls_ov%unode=ls_osolve%u
+ls_ov%vnode=ls_osolve%v
+ls_ov%wnode=ls_osolve%w
+ls_ov%wnodeiso=ls_osolve%wiso
+ls_ov%icon=ls_osolve%icon
+
+write (*,*) 'Before calls to move surface'
+
+do i=1,params%ns
+  if (current_time+tiny(current_time).ge. ls_surf(i)%activation_time) then
+    write (*,*) 'Before call to move surface ',i
+    ! Need to make sure to use params%dt from the ls model here!!!
+    call move_surface (ls0_surf(i),ls0_surf(i),0,0,ls_ov,params%dt,params,istep,i)
+    write (*,*) 'After call to move surface ',i
+  else
+    ls_surf(i)%x=ls_surf(1)%x
+    ls_surf(i)%y=ls_surf(1)%y
+    ls_surf(i)%z=ls_surf(1)%z
+    ls_surf(i)%xn=ls_surf(1)%xn
+    ls_surf(i)%yn=ls_surf(1)%yn
+    ls_surf(i)%zn=ls_surf(1)%zn
+    ls_surf(i)%r=ls_surf(1)%r
+    ls_surf(i)%s=ls_surf(1)%s
+  endif
+
+  ! We apply erosion
+  if (material0.eq.0.and.params%erosion) then
+    if (current_time+tiny(current_time).ge. ls_surf(i)%activation_time) then    
+      call erosion (ls_surf(i),olsf,i,params)
+    else
+      if (iproc.eq.0) write(8,*) 'Surface ',i,' is attached to Surface0'
+      ls_surf(i)%z=ls_surf(1)%z   
+    endif
+  end if
+  if (iproc.eq.0) write (8,*) 'Min-Max z surf ',i,':',minval(ls_surf(i)%z),maxval(ls_surf(i)%z)
+
+  allocate(ls_surf(i)%u(ls_surf(i)%nsurface))
+  allocate(ls_surf(i)%v(ls_surf(i)%nsurface))
+  allocate(ls_surf(i)%w(ls_surf(i)%nsurface))
+
 enddo
 
-deallocate(ls_surf,ls_surf0)
+! Need to do isostasy???
+
+! Write new output file
+call write_global_output (params,istep,iter,current_time,ls_osolve,ls_ov,ls_vo,ls_surf,ls_cl,zi,nest(1),'final')
+
+deallocate(ls_surf,ls_surf0,ss_surf,ss_surf0)
 
 call mpi_finalize (ierr)