post.f90

      if (abs(ov%strain(rtf(i))).le.1.d-99) ov%strain(rtf(i))=0.d0
      write(iunit,'(e13.6)') ov%strain(rtf(i))
   enddo
end if

if(output_temp_field==1) then 
   write(iunit,'(a)')'SCALARS t double 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      if (abs(ov%temp(rtf(i))).ge.1.d98) then
        ov%temp(rtf(i))=sign(1.d98,ov%temp(rtf(i)))
      endif
      if (abs(ov%temp(rtf(i))).le.1.d-99) ov%temp(rtf(i))=0.d0
      write(iunit,'(e13.6)') ov%temp(rtf(i))
   enddo
end if

if(output_lsf_field==1) then 
   do ilsf=1,ov%nlsf
      write (clsf,'(i2)') ilsf
      if (ilsf.lt.10) clsf(1:1)='0'
      write(iunit,'(a)')'SCALARS lsf'//clsf//' double 1'
      write(iunit,'(a)')'LOOKUP_TABLE default'
      do i=1,nnn
         if (abs(ov%lsf(rtf(i),ilsf)).ge.1.d98) then
           ov%lsf(rtf(i),ilsf)=sign(1.d98,ov%lsf(rtf(i),ilsf))
         endif
         if (abs(ov%lsf(rtf(i),ilsf)).le.1.d-99) ov%lsf(rtf(i),ilsf)=0.d0
         write(iunit,'(e13.6)') ov%lsf(rtf(i),ilsf)
      enddo
   enddo
end if

if(output_velo_vect==1) then 
   write(iunit,'(a)')'VECTORS velo double' 
   do i=1,nnn 
      if (abs(ov%u(rtf(i))).ge.1.d98) ov%u(rtf(i))=sign(1.d98,ov%u(rtf(i)))
      if (abs(ov%u(rtf(i))).le.1.d-99) ov%u(rtf(i))=0.d0
      if (abs(ov%v(rtf(i))).ge.1.d98) ov%v(rtf(i))=sign(1.d98,ov%v(rtf(i)))
      if (abs(ov%v(rtf(i))).le.1.d-99) ov%v(rtf(i))=0.d0
      if (abs(ov%w(rtf(i))).ge.1.d98) ov%w(rtf(i))=sign(1.d98,ov%w(rtf(i)))
      if (abs(ov%w(rtf(i))).le.1.d-99) ov%w(rtf(i))=0.d0
      write(iunit,'(3e17.6)')ov%u(rtf(i)),ov%v(rtf(i)),ov%w(rtf(i))
   enddo 
endif 

if(output_preiso_velo_vect==1) then 
   write(iunit,'(a)')'VECTORS preiso-velo double' 
   do i=1,nnn 
      if (abs(ov%u(rtf(i))).ge.1.d98) ov%u(rtf(i))=sign(1.d98,ov%u(rtf(i)))
      if (abs(ov%u(rtf(i))).le.1.d-99) ov%u(rtf(i))=0.d0
      if (abs(ov%v(rtf(i))).ge.1.d98) ov%v(rtf(i))=sign(1.d98,ov%v(rtf(i)))
      if (abs(ov%v(rtf(i))).le.1.d-99) ov%v(rtf(i))=0.d0
      wmod=ov%w(rtf(i))-ov%wiso(rtf(i))
      if (abs(wmod).ge.1.d98) wmod=sign(1.d98,wmod)
      if (abs(wmod).le.1.d-99) wmod=0.d0
      write(iunit,'(3e17.6)')ov%u(rtf(i)),ov%v(rtf(i)),wmod
   enddo 
endif

if(output_iso_velo_vect==1) then 
   write(iunit,'(a)')'VECTORS iso-velo double' 
   do i=1,nnn 
      if (abs(ov%wiso(rtf(i))).ge.1.d98) ov%wiso(rtf(i))=sign(1.d98,ov%wiso(rtf(i)))
      if (abs(ov%wiso(rtf(i))).le.1.d-99) ov%wiso(rtf(i))=0.d0
      write(iunit,'(3e17.6)')0.d0,0.d0,ov%wiso(rtf(i))
   enddo 
endif

if (output_strain_tensor/=0) then
   allocate (strainn(3,3,nnn),nstrain(nnn))
   strainn=0.d0
   nstrain=0
   do ie=1,nnne
      if (elvoid(ie).eq.0) then
         do k=1,nne
            strainn(:,:,ftr(ov%icon(k,ie)))=strainn(:,:,ftr(ov%icon(k,ie)))+strain(:,:,ie)
            nstrain(ftr(ov%icon(k,ie)))=nstrain(ftr(ov%icon(k,ie)))+1
         enddo
      endif
   enddo
   do i=1,nnn
      if (nstrain(i).ne.0) strainn(:,:,i)=strainn(:,:,i)/nstrain(i)
   enddo
   write(iunit,'(a)')'TENSORS strain double'
   do i=1,nnn
      !if (strainn(1,:,i).ge.1.d98) strainn(1,:,i)=1.d98
      !if (strainn(1,:,i).le.1.d-99) strainn(1,:,i)=1.d0
      !if (strainn(2,:,i).ge.1.d98) strainn(2,:,i)=1.d98
      !if (strainn(2,:,i).le.1.d-99) strainn(2,:,i)=1.d0
      !if (strainn(3,:,i).ge.1.d98) strainn(3,:,i)=1.d98
      !if (strainn(3,:,i).le.1.d-99) strainn(3,:,i)=1.d0
      write(iunit,'(3e13.6)') strainn(1,:,i)
      write(iunit,'(3e13.6)') strainn(2,:,i)
      write(iunit,'(3e13.6)') strainn(3,:,i)
   enddo
   deallocate (strainn,nstrain)
endif

write(iunit,'(a10,i10)')'CELL_DATA ', nnne 
 
if (output_e2d_fielde==1) then
   write(iunit,'(a)')'SCALARS e2de double 1' 
   write(iunit,'(a)')'LOOKUP_TABLE default' 
   do ie=1,ov%nleaves 
      if (elvoid(ie).eq.0) then 
         if (abs(ov%e2d(ie)).ge.1.d98) ov%e2d(ie)=sign(1.d98,ov%e2d(ie))
         if (abs(ov%e2d(ie)).le.1.d-99) ov%e2d(ie)=0.d0
         write(iunit,'(e13.6)') ov%e2d(ie)
      endif 
   enddo 
end if 

if(output_press_fielde==1) then  
   write(iunit,'(a)')'SCALARS pressure_e double 1' 
   write(iunit,'(a)')'LOOKUP_TABLE default' 
   do ie=1,ov%nleaves 
      if (elvoid(ie).eq.0) then 
         if (abs(ov%pressure(ie)).ge.1.d98) then
           ov%pressure(ie)=sign(1.d98,ov%pressure(ie))
         endif
         if (abs(ov%pressure(ie)).le.1.d-99) ov%pressure(ie)=0.d0
         write(iunit,'(e13.6)') ov%pressure(ie)
      endif 
   enddo 
end if 

if(output_smooth_press_fielde==1) then  
   write(iunit,'(a)')'SCALARS smooth_pressure_e double 1' 
   write(iunit,'(a)')'LOOKUP_TABLE default' 
   do ie=1,ov%nleaves 
      if (elvoid(ie).eq.0) then 
         if (abs(ov%spressure(ie)).ge.1.d98) then
           ov%spressure(ie)=sign(1.d98,ov%spressure(ie))
         endif
         if (abs(ov%spressure(ie)).le.1.d-99) ov%spressure(ie)=0.d0
         write(iunit,'(e13.6)') ov%spressure(ie)
      endif 
   enddo 
end if

if(output_eviscosity_fielde==1) then  
   write(iunit,'(a)')'SCALARS eff_visc_e double 1' 
   write(iunit,'(a)')'LOOKUP_TABLE default' 
   do ie=1,ov%nleaves 
      if (elvoid(ie).eq.0) then 
         if (abs(ov%eviscosity(ie)).ge.1.d98) then
           ov%eviscosity(ie)=sign(1.d98,ov%eviscosity(ie))
         endif
         if (abs(ov%eviscosity(ie)).le.1.d-99) ov%eviscosity(ie)=0.d0
         write(iunit,'(e13.6)') ov%eviscosity(ie)
      endif 
   enddo 
end if 
	 
if(output_is_plastic_fielde==1) then  
   write(iunit,'(a)')'SCALARS is_plastic_e integer 1' 
   write(iunit,'(a)')'LOOKUP_TABLE default' 
   do ie=1,ov%nleaves 
      if (elvoid(ie).eq.0) then 
         write(iunit,'(I10)') ov%is_plastic(ie) 
      endif 
   enddo
end if

if (output_yield_ratio_fielde==1) then
   write(iunit,'(a)')'SCALARS yield_ratio double 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do ie=1,ov%nleaves
      if (elvoid(ie).eq.0) then
         if (abs(ov%yield_ratio(ie)).ge.1.d98) then
           ov%yield_ratio(ie)=sign(1.d98,ov%yield_ratio(ie))
         endif
         if (abs(ov%yield_ratio(ie)).le.1.d-99) ov%yield_ratio(ie)=0.d0
         write(iunit,'(e13.6)') ov%yield_ratio(ie)
      endif
   enddo
end if

if (output_frict_angle_fielde==1) then
   write(iunit,'(a)')'SCALARS frict_angle double 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do ie=1,ov%nleaves
      if (elvoid(ie).eq.0) then
         if (abs(ov%frict_angle(ie)).ge.1.d98) then
           ov%frict_angle(ie)=sign(1.d98,ov%frict_angle(ie))
         endif
         if (abs(ov%frict_angle(ie)).le.1.d-99) ov%frict_angle(ie)=0.d0
         write(iunit,'(e13.6)') ov%frict_angle(ie)
      endif
   enddo
end if

if (output_dilatr_fielde==1) then
   write(iunit,'(a)')'SCALARS dilatation_rate double 1' 
   write(iunit,'(a)')'LOOKUP_TABLE default' 
   do ie=1,ov%nleaves 
      if (elvoid(ie).eq.0) then 
         if (abs(ov%dilatr(ie)).ge.1.d98) then
           ov%dilatr(ie)=sign(1.d98,ov%dilatr(ie))
         endif
         if (abs(ov%dilatr(ie)).le.1.d-99) ov%dilatr(ie)=0.d0
         write(iunit,'(e13.6)') ov%dilatr(ie)
      endif 
   enddo 
end if

if (output_elem_matnum==1) then
   write(iunit,'(a)')'SCALARS material_number integer 1' 
   write(iunit,'(a)')'LOOKUP_TABLE default' 
   do ie=1,ov%nleaves 
      if (elvoid(ie).eq.0) then 
         write(iunit,'(I10)') ov%matnum(ie)
      endif 
   enddo 
end if

close(iunit)
write(*,*) '--------------------------------------------------------------------------'
write(*,*)'opla I am done producing visu.vtk '

deallocate(countnode)

!==============================================================================
!======[produce norm<i>.vtk file]==============================================
!==============================================================================

if (output_surface_icon==1) then 
   do is=1,ov%nlsf
      write(cs,'(i1)') is
      open(unit=30,file='norm'//cs//'.vtk')
      write(30,'(a)')'# vtk DataFile Version 3.0'
      write(30,'(a)')'surface'
      write(30,'(a)')'ASCII'
      write(30,'(a)')'DATASET UNSTRUCTURED_GRID'
      write(30,'(a7,i10,a6)')'POINTS ',surface(is)%nsurface,' float'
      do i=1,surface(is)%nsurface
         write(30,'(3f16.11)')surface(is)%x(i),surface(is)%y(i),surface(is)%z(i)
      enddo
      write(30,'(a6, 2I10)') 'CELLS ',surface(is)%nt,(1+3)*surface(is)%nt
      do i=1,surface(is)%nt
         write(30,'(9I10)')3,surface(is)%icon(1:3,i)-1
      enddo
      write(30,'(A11, I10)') 'CELL_TYPES ',surface(is)%nt
      do i=1,surface(is)%nt
         write(30,'(I2)')5 ! octree  (8 nodes)
      end do
      write(30,'(a11,i10)') 'POINT_DATA ',surface(is)%nsurface
      write(30,'(a)')'SCALARS z float 1'
      write(30,'(a)')'LOOKUP_TABLE default'
      do i=1,surface(is)%nsurface
         write(30,'(e13.6)') surface(is)%z(i)
      enddo
if (is.eq.1) then
      write(30,'(a)')'SCALARS erosion_rate float 1'
      write(30,'(a)')'LOOKUP_TABLE default'
      do i=1,surface(is)%nsurface
         write(30,'(e13.6)') -surface(is)%r(i)
      enddo
      write(30,'(a)')'SCALARS rock_uplift float 1'
      write(30,'(a)')'LOOKUP_TABLE default'
      do i=1,surface(is)%nsurface
         write(30,'(e13.6)') surface(is)%w(i)
      enddo
endif
      write(30,'(a)') 'VECTORS normal float'
      do i=1,surface(is)%nsurface
         write(30,'(3e17.6)') surface(is)%xn(i),surface(is)%yn(i),surface(is)%zn(i)
      end do
      write(30,'(a)') 'VECTORS vel float'
      do i=1,surface(is)%nsurface
         write(30,'(3e17.6)') surface(is)%u(i),surface(is)%v(i),surface(is)%w(i)
      end do
      close(30)
      write(*,*) '--------------------------------------------------------------------------'
      write(*,*)'opla I am done producing norm'//cs//'.vtk' 
   enddo
end if

!==============================================================================
!======[produce rivers.vtk file]==============================================
!==============================================================================

if (output_rivers==1) then
allocate (donor(surface(1)%nsurface),order(surface(1)%nsurface))
call rivers (surface(1)%icon,surface(1)%nt,surface(1)%x,surface(1)%y,surface(1)%z, &
             surface(1)%nsurface,donor,order)
iordermin=3
nlinks=0
do i=1,surface(1)%nsurface
   if (order(i).gt.iordermin.and.donor(i).gt.0) nlinks=nlinks+1
enddo
      open(unit=30,file='rivers.vtk')
      write(30,'(a)')'# vtk DataFile Version 3.0'
      write(30,'(a)')'rivers'
      write(30,'(a)')'ASCII'
      write(30,'(a)')'DATASET UNSTRUCTURED_GRID'
      write(30,'(a7,i10,a6)')'POINTS ',surface(1)%nsurface,' float'
      do i=1,surface(1)%nsurface
         write(30,'(3f16.11)')surface(1)%x(i),surface(1)%y(i),surface(1)%z(i)+.001
      enddo
      write(30,'(a6, 2I10)') 'CELLS ',nlinks,(1+2)*nlinks
      do i=1,surface(1)%nsurface
         if (order(i).gt.iordermin.and.donor(i).gt.0) write(30,'(9I10)')2,i-1,donor(i)-1
      enddo
      write(30,'(A11, I10)') 'CELL_TYPES ',nlinks
      do i=1,nlinks
         write(30,'(I2)')3 ! octree  (2 nodes)
      end do
      close(30)
      write(*,*) '--------------------------------------------------------------------------'
      write(*,*)'opla I am done producing rivers.vtk' 
deallocate (donor,order)
end if

!==============================================================================
!======[produce regular.vtk file]==============================================
!==============================================================================

if (output_regular==1) then

allocate (levs(ov%nleaves))
call octree_find_element_level (ov%octree,ov%noctree,levs,ov%nleaves)
levmax=maxval(levs)

zmax=maxval(surface(1)%z)
dz=1.d0/(2.d0**levmax)
nz=(zmax/dz)+1

allocate (str11(ov%nnode),str12(ov%nnode),str13(ov%nnode), &
          str22(ov%nnode),str23(ov%nnode),str33(ov%nnode),nstrain(ov%nnode))

nstrain=0
str11=0.;str12=0.;str13=0.
str22=0.;str23=0.;str33=0.
  do ie=1,ov%nleaves
    if (elvoid(ie).eq.0) then
    do k=1,nne
    i=ov%icon(k,ie)
    str11(i)=str11(i)+strain(1,1,ie)
    str12(i)=str12(i)+strain(1,2,ie)
    str13(i)=str13(i)+strain(1,3,ie)
    str22(i)=str22(i)+strain(2,2,ie)
    str23(i)=str23(i)+strain(2,3,ie)
    str33(i)=str33(i)+strain(3,3,ie)
    nstrain(i)=nstrain(i)+1
    enddo
    endif
  enddo
do i=1,ov%nnode
  if (nstrain(i).ne.0) then
  str11(i)=str11(i)/nstrain(i)
  str12(i)=str12(i)/nstrain(i)
  str13(i)=str13(i)/nstrain(i)
  str22(i)=str22(i)/nstrain(i)
  str23(i)=str23(i)/nstrain(i)
  str33(i)=str33(i)/nstrain(i)
  endif
enddo

!print*,minval(str11),maxval(str11)
!print*,minval(str12),maxval(str12)
!print*,minval(str13),maxval(str13)
!print*,minval(str22),maxval(str22)
!print*,minval(str23),maxval(str23)
!print*,minval(str33),maxval(str33)

ni=nz*(2**levmax+1)*(2**levmax+1)
allocate (xi(ni),yi(ni),zi(ni),ui(ni),vi(ni),wi(ni),si(ni),ei(ni),li(ni), &
          s11(ni),s12(ni),s13(ni),s22(ni),s23(ni),s33(ni))
do k=1,nz
do j=1,2**levmax+1
do i=1,2**levmax+1
ijk=ijk+1
xi(ijk)=(i-1)*dz
yi(ijk)=(j-1)*dz
zi(ijk)=(k-1)*dz
call octree_interpolate_many (12,ov%octree,ov%noctree,ov%icon,ov%nleaves,ov%nnode, &
                              xi(ijk),yi(ijk),zi(ijk), &
                              ov%u,ui(ijk),ov%v,vi(ijk),ov%w,wi(ijk), &
                              ov%strain,si(ijk),ov_nodee2d,ei(ijk),ov%lsf(1,1),li(ijk), &
                              str11,s11(ijk),str12,s12(ijk),str13,s13(ijk), &
                              str22,s22(ijk),str23,s23(ijk),str33,s33(ijk))
enddo
enddo
enddo

!print*,minval(s11),maxval(s11)
!print*,minval(s12),maxval(s12)
!print*,minval(s13),maxval(s13)
!print*,minval(s22),maxval(s22)
!print*,minval(s23),maxval(s23)
!print*,minval(s33),maxval(s33)

deallocate (str11,str12,str13,str22,str23,str33,strain)

do i=1,2**levmax+1
do j=1,2**levmax+1
do k=1,nz
ijk=ijk+1
!  if (li(ijk).gt.0.) then
!  if (k.eq.1) stop 'error in lsf1'
!  ui(ijk)=ui(ijk-1)
!  vi(ijk)=vi(ijk-1)
!  wi(ijk)=wi(ijk-1)
!  si(ijk)=si(ijk-1)
!  ei(ijk)=ei(ijk-1)
!  endif
enddo
enddo
enddo

allocate (azimuth1(ni),azimuth3(ni),dip1(ni),dip3(ni))

azimuth1=0.d0
azimuth3=0.d0
dip1=0.d0
dip3=0.d0

do i=1,ni
call PS (s11(i),s22(i),s33(i),s12(i),s23(i),s13(i),l1,l2,l3, &
         n11,n12,n13, &
         n21,n22,n23, &
         n31,n32,n33)
if (n11.ne.0.d0 .or. n12.ne.0.d0) azimuth1(i)=atan2(n12,n11)
if (abs(n13).le.1.d0) dip1(i)=asin(n13)
if (n31.ne.0.d0 .or. n32.ne.0.d0) azimuth3(i)=atan2(n32,n31)
if (abs(n33).le.1.d0) dip3(i)=asin(n33)
enddo

con=45.d0/atan(1.d0)
azimuth1=azimuth1*con
azimuth3=azimuth3*con
dip1=dip1*con
dip3=dip3*con

      open(unit=30,file='regular.vtk')
      write(30,'(a)')'# vtk DataFile Version 3.0'
      write(30,'(a)')'regular'
      write(30,'(a)')'ASCII'
      write(30,'(a)')'DATASET STRUCTURED_POINTS'
      write(30,'(a,3i10)') 'DIMENSIONS',2**levmax+1,2**levmax+1,nz
      write(30,'(a,3f16.11)') 'ORIGIN',0.,0.,0.
      write(30,'(a,3f16.11)') 'SPACING',dz,dz,dz
      write(30,'(a11,i10)')'POINT_DATA ',(2**levmax+1)*(2**levmax+1)*nz
      write(30,'(a)')'VECTORS Velocity float'
        do i=1,(2**levmax+1)*(2**levmax+1)*nz
        write(30,'(3e15.4)') ui(i),vi(i),wi(i)
        enddo
      write(30,'(a)')'SCALARS Strain float 1'
      write(30,'(a)')'LOOKUP_TABLE default'
        do i=1,(2**levmax+1)*(2**levmax+1)*nz
        write(30,'(e13.6)') si(i)
        enddo
      write(30,'(a)')'SCALARS StrainRate float 1'
      write(30,'(a)')'LOOKUP_TABLE default'
        do i=1,(2**levmax+1)*(2**levmax+1)*nz
        write(30,'(e13.6)') ei(i)
        enddo
!      write(30,'(a)')'SCALARS azimuth_sigma1 float 1'
!      write(30,'(a)')'LOOKUP_TABLE default'
!        do i=1,(2**levmax+1)*(2**levmax+1)*nz
!        write(30,'(e13.6)') azimuth1(i)
!        enddo
!      write(30,'(a)')'SCALARS dip_sigma1 float 1'
!      write(30,'(a)')'LOOKUP_TABLE default'
!        do i=1,(2**levmax+1)*(2**levmax+1)*nz
!        write(30,'(e13.6)') dip1(i)
!        enddo
!      write(30,'(a)')'SCALARS azimuth_sigma3 float 1'
!      write(30,'(a)')'LOOKUP_TABLE default'
!        do i=1,(2**levmax+1)*(2**levmax+1)*nz
!        write(30,'(e13.6)') azimuth3(i)
!        enddo
!      write(30,'(a)')'SCALARS dip_sigma3 float 1'
!      write(30,'(a)')'LOOKUP_TABLE default'
!        do i=1,(2**levmax+1)*(2**levmax+1)*nz
!        write(30,'(e13.6)') dip3(i)
!        enddo
      write(30,'(a)')'TENSORS StrainDot float'
        do i=1,(2**levmax+1)*(2**levmax+1)*nz
        write(30,'(3e14.6)') s11(i),s12(i),s13(i)
        write(30,'(3e14.6)') s12(i),s22(i),s23(i)
        write(30,'(3e14.6)') s13(i),s23(i),s33(i)
        enddo
      close(30)

deallocate (xi,yi,zi,ui,vi,wi,si,ei,li,s11,s12,s13,s22,s23,s33)

write(*,*) '--------------------------------------------------------------------------'
write(*,*)'opla I am done producing regular.vtk '
write(*,*) '--------------------------------------------------------------------------'

endif


!==============================================================================
!======[produce zisodisp.vtk file]=============================================
!==============================================================================

if (output_zisodisp==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 ',int((nb+1)**2.d0),' float'
  do j=1,nb+1
    do i=1,nb+1
      write(30,'(3f16.11)') (i-1)*dxy,(j-1)*dxy,zisodisp(i,j)
    enddo
  enddo
  write(30,'(a6, 2I10)') 'CELLS ',int(nb**2.d0),int((4+1)*nb**2.d0)
  ptcnt=0
  do j=1,nb
    do i=1,nb
      write(30,'(9I10)')4,ptcnt,ptcnt+1,ptcnt+1+(nb+1),ptcnt+(nb+1)
      ptcnt=ptcnt+1
    enddo
    ptcnt=ptcnt+1
  enddo
  write(30,'(A11, I10)') 'CELL_TYPES ',int(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 ',int((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,'(e13.6)') zisodisp(i,j)
    enddo
  enddo
  !write(30,'(a11,i10)') 'POINT_DATA ',int((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,'(e13.6)') (zisodisp(i+1,j)-zisodisp(i,j))/dxy
      elseif (i==nb+1) then
        write(30,'(e13.6)') (zisodisp(i,j)-zisodisp(i-1,j))/dxy
      else
        write(30,'(e13.6)') ((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