Added debug output of zisodisp

src/VTK/post.f90

@@ -82,7 +82,7 @@ integer iproc,nproc,ierr,nnn,nnne,nnnep,levelmax,ix,iy,iz,nlinks,iordermin
 integer iconmin,iconmax
 integer,dimension(:),allocatable::nstrain,donor,order
 integer,dimension(:),allocatable::levs,li
-integer nz,ni,ijk
+integer nz,ni,ijk,ptcnt,nb,dumpi
 integer,dimension(:),allocatable::invoid,elvoid,rtf,ftr
 integer myicon(100),nstep,ndir,iter,ii,lsf,levmax
 
@@ -100,9 +100,9 @@ double precision,dimension(:,:,:),allocatable::strain,strainn
 double precision,dimension(:),allocatable::xi,yi,zi,ui,vi,wi,si,ei
 double precision,dimension(:),allocatable::s11,s12,s13,s22,s23,s33,str11,str12,str13,str22,str23,str33
 double precision,dimension(:),allocatable::azimuth1,azimuth3,dip1,dip3
-double precision zmax,dz,l1,l2,l3,n11,n12,n13,n21,n22,n23,n31,n32,n33,con
+double precision zmax,dz,l1,l2,l3,n11,n12,n13,n21,n22,n23,n31,n32,n33,con,dxy
 
-double precision dxyz,x,y,z,xcut
+double precision dxyz,x,y,z,xcut,dumpdp
 integer icut 
 
 
@@ -349,6 +349,25 @@ end if
 
 end do
 
+   read (7) (dumpdp,      & 
+                dumpdp,    &
+                dumpdp,    &
+                dumpdp,    &
+                dumpdp,    & 
+                dumpdp,    &
+                dumpdp,    &
+                dumpdp,    &
+                dumpdp,    &
+                dumpdp,    &
+                dumpi,     &
+                i=1,cl%npcl)
+
+   ! read isostasy basal displacement array - dwhipp 11/09
+   dxy=abs(ov%x(1)-ov%x(2))
+   nb=1.d0/dxy
+   allocate(zisodisp(nb+1,nb+1),zisoslx(nb+1,nb+1))
+   if (1==1) read (7) ((zisodisp(i,j),j=1,nb+1),i=1,nb+1)
+
 close(7)
 
 
@@ -1215,6 +1234,66 @@ write(*,*) '--------------------------------------------------------------------
 
 endif
 
+
+!==============================================================================
+!======[produce zisodisp.vtk file]=============================================
+!==============================================================================
+
+if (1==1) then
+  open(unit=30,file='zisodisp.vtk')
+  write(30,'(a)')'# vtk DataFile Version 3.0'
+  write(30,'(a)')'zisodisp'
+  write(30,'(a)')'ASCII'
+  write(30,'(a)')'DATASET UNSTRUCTURED_GRID'
+  write(30,'(a7,i10,a6)')'POINTS ',(nb+1)**2.d0,' float'
+  do j=1,nb+1
+    do i=1,nb+1
+      write(30,'(3f16.11)') ((i-1)*dxy)/nb,((j-1)*dxy)/nb,zisodisp(i,j)
+    enddo
+  enddo
+  write(30,'(a6, 2I10)') 'CELLS ',nb**2.d0,(4+1)*nb**2.d0
+  ptcnt=0
+  do j=1,nb
+    do i=1,nb
+      ptcnt=ptcnt+1
+      write(30,'(9I10)')4,ptcnt,ptcnt+1,ptcnt+1+nb,ptcnt+nb
+    enddo
+    ptcnt=ptcnt+1
+  enddo
+  write(30,'(A11, I10)') 'CELL_TYPES ',nb**2.d0
+  do j=1,nb
+    do i=1,nb
+      write(30,'(I2)')9 ! VTK quad (4 nodes)
+    enddo
+  enddo
+  write(30,'(a11,i10)') 'POINT_DATA ',(nb+1)**2.d0
+  write(30,'(a)')'SCALARS zisodisp float 1'
+  write(30,'(a)')'LOOKUP_TABLE default'
+  do j=1,nb+1
+    do i=1,nb+1
+      write(30,'(e11.4)') zisodisp(i,j)
+    enddo
+  enddo
+  write(30,'(a11,i10)') 'POINT_DATA ',(nb+1)**2.d0
+  write(30,'(a)')'SCALARS zisoslx float 1'
+  write(30,'(a)')'LOOKUP_TABLE default'
+  do j=1,nb+1
+    do i=1,nb+1
+      if (i==1) then
+        write(30,'(e11.4)') (zisodisp(i+1,j)-zisodisp(i,j))/dxy
+      elseif (i==nb+1) then
+        write(30,'(e11.4)') (zisodisp(i,j)-zisodisp(i-1,j))/dxy
+      else
+        write(30,'(e11.4)') ((zisodisp(i,j)-zisodisp(i-1,j))/dxy+&
+                     (zisodisp(i+1,j)-zisodisp(i,j))/dxy)/2.d0
+      endif
+    enddo
+  enddo
+  close(30)
+  write(*,*) '--------------------------------------------------------------------------'
+  write(*,*)'opla I am done producing zisodisp.vtk' 
+end if
+
 write(*,*) '**************************************************************************'
 
 end