Added storage of zisodisp array

src/isostasy.f90

-subroutine isostasy (params,weightel,ov,surface,mat)
+subroutine isostasy (params,weightel,ov,surface,mat,zisodisp)
 
 ! first go at imposing isostasy
 ! this routine assumes that the original surface are flat
@@ -15,6 +15,7 @@ type(octreev) ov
 double precision weightel(ov%nleaves)
 type(sheet) surface(params%ns)
 type(material) mat(0:params%nmat)
+double precision zisodisp(2**params%levelmax_oct+1,2**params%levelmax_oct+1)
 
 double precision xc,yc,zc,eps
 integer i,j,ii,jj,kk,ie,dlev,is,iproc,nproc,ierr,ic
@@ -107,7 +108,7 @@ xk=-mat(surface(params%ns)%material)%density*params%density_difference*9.81d0
 d=1.d11/12.d0/(1.d0-0.25d0**2)*params%elastic_plate_thickness**3
 dh=hx/(nflex-1)
 
-if (iproc.eq.0) print*,'hx,dh,xk,d',hx,dh,xk,d
+!if (iproc.eq.0) print*,'hx,dh,xk,d',hx,dh,xk,d
 
 ii=nflex/2-nb/2-1
 jj=ii
@@ -118,7 +119,7 @@ flex=0.d0
     enddo
   enddo
 
-if (iproc.eq.0) print*,'poids',minval(flex),maxval(flex)
+!if (iproc.eq.0) print*,'poids',minval(flex),maxval(flex)
 
   do j=1,nflex
   call dsinft (flex(1,j),nflex)
@@ -165,7 +166,7 @@ pihx=pi/hx
   call dsinft (flex(1,j),nflex)
   enddo
 
-if (iproc.eq.0) print*,'def in m',minval(flex),maxval(flex)
+!if (iproc.eq.0) print*,'def in m',minval(flex),maxval(flex)
 
 ii=nflex/2-nb/2-1
 jj=ii
@@ -177,7 +178,7 @@ jj=ii
 
 deallocate (flex,work)
 
-if (iproc.eq.0) print*,'dimensionless def ',minval(weight),maxval(weight)
+!if (iproc.eq.0) print*,'dimensionless def ',minval(weight),maxval(weight)
 
 endif
 
@@ -206,59 +207,62 @@ if ((ii-1)*(ii-nb-1).gt.0 .or. (jj-1)*(jj-nb-1).gt.0) stop 'error in isostasy'
 diso(i)=work(ii,jj)
 enddo
 
-goto 3333
+!goto 3333
 
 ! old version : to be neglected new is much faster
 
-diso=0.d0
-np=0
-  do ie=1,ov%nleaves
-  xc=0.d0;yc=0.d0;zc=0.d0
-    do i=1,8
-    xc=xc+ov%x(ov%icon(i,ie))
-    yc=yc+ov%y(ov%icon(i,ie))
-    zc=zc+ov%z(ov%icon(i,ie))
-    enddo
-  xc=xc/8.d0
-  yc=yc/8.d0
-  zc=zc/8.d0
-  call octree_find_leaf (ov%octree,ov%noctree,xc,yc,zc,leaf,level,loc,x0,y0,z0,dxyz)
-!  call octree_find_leaf (ov%octree,ov%noctree,xc,yc,eps,leafb,levelb,locb,x0b,y0b,z0b,dxyzb)
-!  disoel=(w(leafb)-weightref)/(mat(surface(params%ns)%material)%density)
-  w=0.d0
-  dlev=params%levelmax_oct-level
-  call find_integer_coordinates (x0+dd,y0+dd,z0+dd,ii,jj,kk,params%levelmax_oct)
-  ii=min(ii,2**(params%levelmax_oct)-2**dlev)
-  jj=min(jj,2**(params%levelmax_oct)-2**dlev)
-    do i=1,2**dlev
-      do j=1,2**dlev 
-      w=w+weight(ii+i,jj+j)
-      enddo
-    enddo
-  disoel=w/(2**dlev*2.d0**dlev)
-!  disoel=(w-weightref)/(mat(surface(params%ns)%material)%density)
-    do i=1,8
-    ic=ov%icon(i,ie)
-    np(ic)=np(ic)+1
-    diso(ic)=diso(ic)+disoel
-    enddo
-  enddo
-
-  do i=1,ov%nnode
-    if (np(i).ne.0) then
-    diso(i)=diso(i)/np(i)
-    else
-    diso(i)=0
-    print*,'error np',i,ov%x(i),ov%y(i),ov%z(i)
-    endif
-  enddo
-
-3333 continue
+!diso=0.d0
+!np=0
+!  do ie=1,ov%nleaves
+!  xc=0.d0;yc=0.d0;zc=0.d0
+!    do i=1,8
+!    xc=xc+ov%x(ov%icon(i,ie))
+!    yc=yc+ov%y(ov%icon(i,ie))
+!    zc=zc+ov%z(ov%icon(i,ie))
+!    enddo
+!  xc=xc/8.d0
+!  yc=yc/8.d0
+!  zc=zc/8.d0
+!  call octree_find_leaf (ov%octree,ov%noctree,xc,yc,zc,leaf,level,loc,x0,y0,z0,dxyz)
+!!  call octree_find_leaf (ov%octree,ov%noctree,xc,yc,eps,leafb,levelb,locb,x0b,y0b,z0b,dxyzb)
+!!  disoel=(w(leafb)-weightref)/(mat(surface(params%ns)%material)%density)
+!  w=0.d0
+!  dlev=params%levelmax_oct-level
+!  call find_integer_coordinates (x0+dd,y0+dd,z0+dd,ii,jj,kk,params%levelmax_oct)
+!  ii=min(ii,2**(params%levelmax_oct)-2**dlev)
+!  jj=min(jj,2**(params%levelmax_oct)-2**dlev)
+!    do i=1,2**dlev
+!      do j=1,2**dlev 
+!      w=w+weight(ii+i,jj+j)
+!      enddo
+!    enddo
+!  disoel=w/(2**dlev*2.d0**dlev)
+!!  disoel=(w-weightref)/(mat(surface(params%ns)%material)%density)
+!    do i=1,8
+!    ic=ov%icon(i,ie)
+!    np(ic)=np(ic)+1
+!    diso(ic)=diso(ic)+disoel
+!    enddo
+!  enddo
+!
+!  do i=1,ov%nnode
+!    if (np(i).ne.0) then
+!    diso(i)=diso(i)/np(i)
+!    else
+!    diso(i)=0
+!    print*,'error np',i,ov%x(i),ov%y(i),ov%z(i)
+!    endif
+!  enddo
+!
+!3333 continue
 
   if (params%debug>=1.and.iproc.eq.0) then
   write(*,'(a,f10.7,1x,f10.7)') shift//'Isostatic velocity range = ',minval(diso)/params%dt,maxval(diso)/params%dt
   endif
 
+! Store the total displacement of the basal node plane 
+if (params%isobc) zisodisp=zisodisp-work
+
 ! transforms the displacements into velocities and add this contribution to the velocity
 ! solution of the Stokes equation before interpolation onto the surfaces and their advection
 
@@ -267,47 +271,47 @@ np=0
   ov%wnode(i)=ov%wnode(i)-diso(i)/params%dt
   enddo
 
-goto 1111
-
-! display results into a vtk file for debugging
-
-iunit=30
-open(unit=iunit,file='isostasy.vtk')
-write(iunit,'(a)')'# vtk DataFile Version 3.0'
-write(iunit,'(a)')'velocities'
-write(iunit,'(a)')'ASCII'
-write(iunit,'(a)')'DATASET UNSTRUCTURED_GRID'
-write(iunit,'(a7,i10,a6)')'POINTS ',ov%nnode,' float'
-
-do i=1,ov%nnode
-write(iunit,'(3f16.11)')ov%x(i),ov%y(i),ov%z(i)
-enddo
-
-write(iunit,'(A6, 2I10)') 'CELLS ',ov%nleaves,(1+8)*ov%nleaves
-do ie=1,ov%nleaves
-write(iunit,'(9I10)')8,ov%icon(1:8,ie)-1
-enddo
-
-write(iunit,'(A11, I10)') 'CELL_TYPES ',ov%nleaves
-do ie=1,ov%nleaves
-   write(iunit,'(I2)')11 ! octree  (8 nodes)
-end do
-
-write(iunit,'(a11,i10)')'POINT_DATA ',ov%nnode
-
-write(iunit,'(a)')'VECTORS iso float'
-do i=1,ov%nnode
-   write(iunit,'(3e11.4)')0.,0.,-diso(i)/params%dt
-enddo
-
-write(iunit,'(a)')'VECTORS velo float'
-do i=1,ov%nnode
-   write(iunit,'(3e11.4)')0.,0.,ov%wnode(i)
-enddo
-
-close (iunit)
-
-1111 continue
+!goto 1111
+!
+!! display results into a vtk file for debugging
+!
+!iunit=30
+!open(unit=iunit,file='isostasy.vtk')
+!write(iunit,'(a)')'# vtk DataFile Version 3.0'
+!write(iunit,'(a)')'velocities'
+!write(iunit,'(a)')'ASCII'
+!write(iunit,'(a)')'DATASET UNSTRUCTURED_GRID'
+!write(iunit,'(a7,i10,a6)')'POINTS ',ov%nnode,' float'
+!
+!do i=1,ov%nnode
+!write(iunit,'(3f16.11)')ov%x(i),ov%y(i),ov%z(i)
+!enddo
+!
+!write(iunit,'(A6, 2I10)') 'CELLS ',ov%nleaves,(1+8)*ov%nleaves
+!do ie=1,ov%nleaves
+!write(iunit,'(9I10)')8,ov%icon(1:8,ie)-1
+!enddo
+!
+!write(iunit,'(A11, I10)') 'CELL_TYPES ',ov%nleaves
+!do ie=1,ov%nleaves
+!   write(iunit,'(I2)')11 ! octree  (8 nodes)
+!end do
+!
+!write(iunit,'(a11,i10)')'POINT_DATA ',ov%nnode
+!
+!write(iunit,'(a)')'VECTORS iso float'
+!do i=1,ov%nnode
+!   write(iunit,'(3e11.4)')0.,0.,-diso(i)/params%dt
+!enddo
+!
+!write(iunit,'(a)')'VECTORS velo float'
+!do i=1,ov%nnode
+!   write(iunit,'(3e11.4)')0.,0.,ov%wnode(i)
+!enddo
+!
+!close (iunit)
+!
+!1111 continue
 
 deallocate (np)
 deallocate (diso)