Fixed typos

src/calculate_density_profile.f90

@@ -178,7 +178,7 @@ do i=size_str,1,-1
        do k=1,osolve%nlsf
          ! unsquished z into octree_interpolate
          call octree_interpolate (osolve%octree,osolve%noctree,osolve%icon,         &
-                                  osolve%nleaves,osolve%lsf(:,i),osolve%nnode,      &
+                                  osolve%nleaves,osolve%lsf(:,k),osolve%nnode,      &
                                   x,y,z/params%vex,cur_lsf)
          if (cur_lsf < eps) then
            matnum = k

src/isostasy.f90

-subroutine isostasy (params,weightel,ov,surface,mat,flag0,bcdef,istep)
+subroutine isostasy (params,weightel,ov,surface,mat,flag0,bcdef,istep,iter)
 
 ! first go at imposing isostasy
 ! this routine assumes that the original surface are flat
@@ -22,6 +22,7 @@ type(material) mat(0:params%nmat)
 integer flag0
 type(bc_definition) bcdef
 integer istep
+integer iter
 
 double precision xc,yc,zc,eps,dtcur
 integer i,j,ii,jj,kk,ie,dlev,is,iproc,nproc,ierr
@@ -88,8 +89,8 @@ do ie=1,ov%nleaves
 !  w(leafb)=w(leafb)+weightel(leaf)/(4.d0**dlev)*(4.d0**levelb)
 enddo
 
-! Store minimum column weight as reference weight in first isostasy call
-if (params%first_iter) then
+! Store minimum column weight as reference weight in first isostasy call of run
+if (istep==params%irestart+1 .and. iter==1) then
   weightrefloc=minloc(abs(weight))
   params%weightref = weight(weightrefloc(1),weightrefloc(2))
 endif
@@ -111,7 +112,6 @@ weightref=params%weightref
 !write (*,*) 'weightref mod: ',weightref
 
 if (.not.params%flexure) then
-
   ! apply local isostasy
   ! the deflection is stored in array weight
   do j=1,nb
@@ -120,9 +120,7 @@ if (.not.params%flexure) then
       weight(i,j)=(weight(i,j)-weightref)/mat(surface(params%nmat)%material)%density
     enddo
   enddo
-
 else
-
   ! flexural isostasy
   nflex=8*nb
   allocate(work(nflex,nflex),flex(nflex,nflex))
@@ -237,7 +235,7 @@ enddo
 !  write (*,*) 'displacement max: ',maxval(diso)
 !endif
 
-if (bcdef%bctype == 'iso_only' .and. params%first_iter) then
+if (bcdef%bctype == 'iso_only' .and. istep==params%irestart+1 .and. iter==1) then
   if (maxval(diso) < (params%dt * isovelomax)) then
     params%nstep = 1
     !write (*,*) 'Maxval(diso) less than dt*isovelomax'