
MODULE definitions

   type structure
      integer,dimension(:),pointer::octree
      integer noctree,nnode,nleaves
      integer nelemr,nlsf 
      double precision,dimension(:),pointer::x,y,z
      double precision,dimension(:),pointer::u,v,w,wiso,temp,pressure,strain
      double precision,dimension(:),pointer::nodal_pressure,spressure
      double precision,dimension(:),pointer::e2d,eviscosity
      double precision,dimension(:),pointer::crit,dilatr
      integer,dimension(:,:),pointer::icon,iconr
      double precision,dimension(:,:),pointer:: lsf
      integer,dimension(:),pointer::on(:),is_plastic,matnum
      logical, dimension(:), pointer :: whole_leaf_in_fluid
   end type structure

   type sheet
      integer nsurface,nt
      double precision,dimension(:),pointer::r,s,x,y,z,xn,yn,zn,u,v,w
      integer,dimension(:,:),pointer::icon
      integer nhull,levelt,itype,material
   end type sheet

END MODULE definitions


!==============================================================================!
!==============================================================================!
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!==============================================================================!
!==============================================================================!
!                                                                              |
!              POST PROCESSING    Apr. 2008                                    |
!                                                                              |
!==============================================================================!
!==============================================================================!

program data_postprocessing

use definitions

implicit none

type (sheet),dimension(:),allocatable::surface
type(structure),SAVE :: ov
integer :: ie,mer,nn,ner,i,i1,i2,i3,i4,i5,kfx,kfy,kfz,kft
integer :: iunit,ilsf,nne,is,k,npcl,j,np,nl
integer :: output_u_field
integer :: output_v_field
integer :: output_w_field
integer :: output_velo_vect 
integer :: output_preiso_velo_vect 
integer :: output_iso_velo_vect 
integer :: output_disp_field
integer :: output_press_fieldn 
integer :: output_press_fielde 
integer :: output_smooth_press_fielde
integer :: output_countnode_field
integer :: output_temp_field
integer :: output_elem_matnum
integer :: output_e2d_fielde 
integer :: output_e2d_fieldn
integer :: output_eviscosity_fielde 
integer :: output_is_plastic_fielde
integer :: output_dilatr_fielde
integer :: output_crit_field
integer :: output_strain_field
integer :: output_lsf_field
integer :: output_surface_icon
integer :: output_strain_tensor
integer :: output_ps
integer :: output_cubes
integer :: output_rivers
integer :: output_regular
integer :: output_zisodisp
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,ptcnt,nb,dumpi
integer,dimension(:),allocatable::invoid,elvoid,rtf,ftr
integer myicon(100),nstep,ndir,iter,ii,lsf,levmax

logical,dimension(:),allocatable::influid,do_it,subset_leaves,instrain
logical :: nest,nestscl

character clsf*3,c4*4,cc4*4,dir*128
character cs,nestin,nestsclin

double precision :: eps,dil,current_time,activation_time,zmin,xx,yy,zz,maxe2d,dist
double precision, allocatable, dimension (:) :: ov_nodee2d, ov_nodecrit
double precision, allocatable, dimension (:) :: countnode 
double precision,dimension(:,:),allocatable::s,n1,n2,n3,zisodisp
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,dxy
double precision :: sselemx,sselemy,sselemz,xminls,yminls,zminls

double precision dxyz,x,y,z,xcut,dumpdp
integer icut 


!==============================================================================!
!==============================================================================!

call system('clear')

nne=8

print*,'Enter directory (default is OUT)'
read (*,'(a)') dir
ndir=len_trim(dir)
if (ndir.eq.0) then
  ndir=3
  dir(1:ndir)='OUT'
endif

print*,'Enter time step >'
read*,nstep
write (c4,'(i4)') nstep
if (nstep.lt.1000) c4(1:1)='0'
if (nstep.lt.100) c4(1:2)='00'
if (nstep.lt.10) c4(1:3)='000'

do i=1,1
  write (*,*) 'Is this a nested model? (y/n)'
  read (*,'(a)') nestin
  select case (trim(nestin))
  case ('y','Y','yes','Yes')
    nest=.true.
    do j=1,1
      write (*,*) 'Should the values be scaled to fit the large(st) scale model? (y/n)'
      read (*,'(a)') nestsclin
      select case (trim(nestsclin))
      case ('y','Y','yes','Yes')
        nestscl=.true.
      case ('n','N','no','No')
        nestscl=.false.
      case default 
        write (*,*) 'Error: Please input either "y" or "n"'
        j=0
      end select
    enddo
  case ('n','N','no','No')
    nest=.false.
  case default 
    write (*,*) 'Error: Please input either "y" or "n"'
    i=0
  end select
enddo

print*,'Is it a debug file ? (y or n)'
read (*,'(a)') cs
select case (cs)
case ('y')
   print*,'Which debug file: please enter grid iteration number'
   read*,iter
   write (cc4,'(i4)') iter
   if (iter.lt.1000) cc4(1:1)='0'
   if (iter.lt.100) cc4(1:2)='00'
   if (iter.lt.10) cc4(1:3)='000'
case ('n')
   print*,'it is not a debug file'
case default 
   print *,'error: please input y or n'
   stop
end select

!==============================================================================
!==============================================================================

write(*,*) '**************************************************************************'
write(*,*) '--------------------------------------------------------------------------'
write(*,*) '--------------------------------------------------------------------------'
write(*,*) '------------------- VTK file production in progress ----------------------'
write(*,*) '--------------------------------------------------------------------------'
write(*,*) '--------------------------------------------------------------------------'
if (cs=='y') then
   write(*,*) 'output file to be processed: ', '../'//dir(1:ndir)//'/time_'//c4//'_'//cc4//'.bin'
   open(unit=7,file='../'//dir(1:ndir)//'/time_'//c4//'_'//cc4//'.bin',status='old',form='unformatted')
else
   write(*,*) 'output file to be processed: ', '../'//dir(1:ndir)//'/time_'//c4//'.bin'
   open(unit=7,file='../'//dir(1:ndir)//'/time_'//c4//'.bin',status='old',form='unformatted')
end if
write(*,*) '--------------------------------------------------------------------------'
open(unit=77,file='input_of_outputs.txt',status='old')
read(77,*) output_u_field 
read(77,*) output_v_field 
read(77,*) output_w_field
read(77,*) output_velo_vect 
read(77,*) output_preiso_velo_vect 
read(77,*) output_iso_velo_vect 
read(77,*) output_disp_field 
read(77,*) output_press_fieldn 
read(77,*) output_press_fielde 
read(77,*) output_smooth_press_fielde
read(77,*) output_countnode_field
read(77,*) output_temp_field
read(77,*) output_elem_matnum
read(77,*) output_e2d_fielde 
read(77,*) output_e2d_fieldn
read(77,*) output_eviscosity_fielde 
read(77,*) output_is_plastic_fielde
read(77,*) output_dilatr_fielde
read(77,*) output_crit_field 
read(77,*) output_strain_field 
read(77,*) output_lsf_field 
read(77,*) output_surface_icon
read(77,*) output_strain_tensor
read(77,*) output_ps
read(77,*) output_cubes
read(77,*) output_rivers
read(77,*) output_regular
read(77,*) output_zisodisp

write(*,*) 'output u      field       ->',(output_u_field==1)
write(*,*) 'output v      field       ->',(output_v_field==1)
write(*,*) 'output w      field       ->',(output_w_field==1)
write(*,*) 'output velo   vectors     ->',(output_velo_vect==1) 
write(*,*) 'output preiso velo vectors->',(output_preiso_velo_vect==1) 
write(*,*) 'output iso velo vectors->',(output_preiso_velo_vect==1) 
write(*,*) 'output disp   field       ->',(output_disp_field==1)
write(*,*) 'output node press field   ->',(output_press_fieldn==1) 
write(*,*) 'output raw press field    ->',(output_press_fielde==1) 
write(*,*) 'output smooth press field ->',(output_smooth_press_fielde==1)
write(*,*) 'output countnode field    ->',(output_countnode_field==1)
write(*,*) 'output temp   field       ->',(output_temp_field==1)
write(*,*) 'output elem material #    ->',(output_elem_matnum==1)
write(*,*) 'output elemental e2d field->',(output_e2d_fielde==1) 
write(*,*) 'output nodal e2d field    ->',(output_e2d_fieldn==1)
write(*,*) 'output elem eff viscosity ->',(output_eviscosity_fielde==1) 
write(*,*) 'output element is plastic ->',(output_is_plastic_fielde==1)
write(*,*) 'output element dilat rate ->',(output_dilatr_fielde==1) 
write(*,*) 'output crit   field       ->',(output_crit_field==1)
write(*,*) 'output strain field       ->',(output_strain_field==1)
write(*,*) 'output lsf    field       ->',(output_lsf_field==1)
write(*,*) 'output surfaces icon      ->',(output_surface_icon==1)
write(*,*) 'output strain tensor      ->',(output_strain_tensor==1)
write(*,*) 'output principal stresses ->',(output_ps==1)
write(*,*) 'output elemental fields   ->',(output_cubes==1)
write(*,*) 'output rivers             ->',(output_rivers==1)
write(*,*) 'output regular            ->',(output_regular==1)
write(*,*) 'output zisodisp           ->',(output_zisodisp==1)
write(*,*) '--------------------------------------------------------------------------'

!==============================================================================!
!==============================================================================!

read (7) ov%noctree,ov%nnode,ov%nleaves,ner,ov%nlsf,npcl,current_time
nn=ov%nnode
allocate(ov%on(nn))
allocate(ov%x(nn),ov%y(nn),ov%z(nn))
allocate(ov%u(nn),ov%v(nn),ov%w(nn),ov%wiso(nn))
allocate(ov%pressure(ov%nleaves),ov%strain(nn))
allocate(ov%e2d(ov%nleaves),ov%eviscosity(ov%nleaves),ov%is_plastic(ov%nleaves))
allocate(ov%dilatr(ov%nleaves),ov%whole_leaf_in_fluid(ov%nleaves))
allocate(ov%crit(ov%nleaves),ov%matnum(ov%nleaves))
allocate(ov%temp(nn))
allocate(ov%lsf(nn,ov%nlsf))
allocate(influid(nn))
allocate(invoid(nn))
allocate(rtf(nn),ftr(nn),elvoid(nn))
allocate(surface(ov%nlsf))
allocate(ov%nodal_pressure(nn),ov%spressure(ov%nleaves))
allocate(ov%icon(8,ov%nleaves))
allocate(ov%octree(ov%noctree))

!========================
!=====[nodal values]=====
!========================

read(7)(ov%x(i),ov%y(i),ov%z(i),   &
        ov%u(i),ov%v(i),ov%w(i),   &
        ov%wiso(i),                &
        ov%lsf(i,1:ov%nlsf),       &
        ov%temp(i),                &
        ov%nodal_pressure(i),      &
        ov%strain(i),              &
        kfx,kfy,kfz,kft,i=1,nn)

!======================
!=====[icon array]=====
!======================

read(7) (ov%icon(1:8,ie),            &
         ov%pressure(ie),            & 
         ov%spressure(ie),           &
         ov%crit(ie),                &
         ov%e2d(ie),                 &
         ov%eviscosity(ie),          &
         ov%is_plastic(ie),          &
         ov%dilatr(ie),              &
         ov%matnum(ie),              &
         ov%whole_leaf_in_fluid(ie), & 
         ie=1,ov%nleaves)

!==============================
!=====[octree information]=====
!==============================

read(7) (ov%octree(i),i=1,ov%noctree)

!================================
!=====[bad face information]=====
!================================
mer=ner   ! mer=1
if(ner.gt.mer)mer=ner
allocate(ov%iconr(9,mer))
read(7) (ov%iconr(1:9,ie),ie=1,ner)

ov%nelemr=ner

!============================
!=====[void information]=====
!============================

read (7) (invoid(i),elvoid(i),ftr(i),rtf(i),influid(i),i=1,nn)
ov%on=0
do i=1,nn
   if (influid(i)) ov%on(i)=1
enddo

read(7) (i5,i=1,ner)

!==============================================================
!=====[measure work, strain, compute principal directions]=====
!==============================================================

allocate (s (3,ov%nleaves))
allocate (n1(3,ov%nleaves))
allocate (n2(3,ov%nleaves))
allocate (n3(3,ov%nleaves))
allocate (strain(3,3,ov%nleaves))

call find_volumes (ov%icon,ov%octree,ov%noctree,ov%nleaves,ov%lsf,ov%nlsf,ov%nnode,influid, &
                   ov%x,ov%y,ov%z,ov%u,ov%v,ov%w,s,n1,n2,n3,strain)

! rescaling of the values for plotting purposes

s=s/maxval(abs(s))*.01

strain=strain/maxval(abs(strain))*.01

!===============================
!=====[surface information]=====
!===============================

do is=1,ov%nlsf
   read(7) surface(is)%nsurface,activation_time,surface(is)%nt
   allocate (surface(is)%x(surface(is)%nsurface),  &
             surface(is)%y(surface(is)%nsurface),  &
             surface(is)%z(surface(is)%nsurface),  &
             surface(is)%xn(surface(is)%nsurface), &
             surface(is)%yn(surface(is)%nsurface), &
             surface(is)%zn(surface(is)%nsurface), &
             surface(is)%r(surface(is)%nsurface),  &
             surface(is)%s(surface(is)%nsurface),  &
             surface(is)%u(surface(is)%nsurface),  &
             surface(is)%v(surface(is)%nsurface),  &
             surface(is)%w(surface(is)%nsurface),  &
             surface(is)%icon(3,surface(is)%nt))
  
if (cs=='y') then
   read (7) (surface(is)%r(i),surface(is)%s(i), &
               surface(is)%x(i),surface(is)%y(i),surface(is)%z(i), &
               surface(is)%xn(i),surface(is)%yn(i),surface(is)%zn(i), &
               i=1,surface(is)%nsurface)
else
   read (7) (surface(is)%r(i),surface(is)%s(i), &
               surface(is)%x(i),surface(is)%y(i),surface(is)%z(i), &
               surface(is)%xn(i),surface(is)%yn(i),surface(is)%zn(i), &
               surface(is)%u(i),surface(is)%v(i),surface(is)%w(i), &
               i=1,surface(is)%nsurface)
end if

   read (7) (surface(is)%icon(1:3,i),i=1,surface(is)%nt)

end do

!=============================
!=====[cloud information]=====
!=============================

read (7) (dumpdp,      & 
            dumpdp,    &
            dumpdp,    &
            dumpdp,    &
            dumpdp,    & 
            dumpdp,    &
            dumpdp,    &
            dumpdp,    &
            dumpdp,    &
            dumpdp,    &
            dumpi,     &
            i=1,npcl)

!================================
!=====[isostasy information]=====
!================================

read (7) nb
dxy=1.d0/real(nb)
allocate(zisodisp(nb+1,nb+1))
read (7) ((zisodisp(i,j),j=1,nb+1),i=1,nb+1)

!============================
!=====[nest information]=====
!============================

read(7) sselemx,sselemy,sselemz,xminls,yminls,zminls

close(7)

!==============================================================================
!==============================================================================
!=====[C -> N]=================================================================
!==============================================================================
!==============================================================================

allocate(countnode(ov%nnode))
allocate(ov_nodee2d(ov%nnode))
allocate(ov_nodecrit(ov%nnode))

ov_nodee2d       = 0.d0
ov_nodecrit      = 0.d0
countnode        = 0.d0
 
do i=1,ov%nleaves
   do j=1,8
      k=ov%icon(j,i)
      ov_nodee2d(k)       = ov_nodee2d(k)      + ov%e2d(i)    
      ov_nodecrit(k)      = ov_nodecrit(k)     + ov%crit(i)    
      countnode(k)        = countnode(k)       + 1                     
   end do
end do

do i=1,ov%nnode
   if (countnode(i).gt.0.d0) then
      ov_nodee2d(i)      = ov_nodee2d(i)      / countnode(i)
      ov_nodecrit(i)     = ov_nodecrit(i)     / countnode(i)
   end if
end do

! Convert coordinates if using a nest
if (nest) then
  ov%x=ov%x*sselemx+xminls
  ov%y=ov%y*sselemy+yminls
  ov%z=ov%z*sselemz+zminls
endif

!==============================================================================
!==============================================================================

write(*,'(a,i10)') 'octree length              ',ov%noctree
write(*,'(a,i10)') 'nb of nodes                ',ov%nnode
write(*,'(a,i10)') 'nb of leaves               ',ov%nleaves
write(*,'(a,i10)') 'nb of bad faces            ',ner
write(*,'(a,i10)') 'nb of lsf                  ',ov%nlsf
write(*,*) '--------------------------------------------------------------------------'
write(*,'(a,2f30.20)') 'x                 range : ',minval(ov%x(1:nn)),       maxval(ov%x(1:nn))
write(*,'(a,2f30.20)') 'y                 range : ',minval(ov%y(1:nn)),       maxval(ov%y(1:nn))
write(*,'(a,2f30.20)') 'z                 range : ',minval(ov%z(1:nn)),       maxval(ov%z(1:nn))
write(*,'(a,2f30.20)') 'u                 range : ',minval(ov%u(1:nn)),       maxval(ov%u(1:nn))
write(*,'(a,2f30.20)') 'v                 range : ',minval(ov%v(1:nn)),       maxval(ov%v(1:nn))
write(*,'(a,2f30.20)') 'w                 range : ',minval(ov%w(1:nn)),       maxval(ov%w(1:nn))
write(*,'(a,2f30.20)') 'w (iso)           range : ',minval(ov%wiso(1:nn)),    maxval(ov%wiso(1:nn))
write(*,'(a,2f30.20)') 'pressure          range : ',minval(ov%pressure),      maxval(ov%pressure) 
write(*,'(a,2f30.20)') 'smooth pressure   range : ',minval(ov%spressure),     maxval(ov%spressure)
write(*,'(a,2f30.20)') 'pressure (nodes)  range : ',minval(ov%nodal_pressure),maxval(ov%nodal_pressure)
write(*,'(a,2f30.20)') 'e2d               range : ',minval(ov%e2d),           maxval(ov%e2d)
write(*,'(a,2f30.20)') 'e2d (nodes)       range : ',minval(ov_nodee2d),       maxval(ov_nodee2d)
write(*,'(a,2f30.20)') 'dilatation rate   range : ',minval(ov%dilatr),        maxval(ov%dilatr)
write(*,'(a,2f30.20)') 'crit              range : ',minval(ov%crit),          maxval(ov%crit)
write(*,'(a,2f30.20)') 'crit (nodes)      range : ',minval(ov_nodecrit),      maxval(ov_nodecrit)
write(*,'(a,2f30.20)') 'strain            range : ',minval(ov%strain(1:nn)),  maxval(ov%strain(1:nn))
do is=1,ov%nlsf
write(*,*) '--------------------------------------------------------------------------'
write(*,'(a,i3)') 'surface',is
write(*,'(a,i7,a)') 'surface counts ',surface(is)%nsurface, ' points'
write(*,'(a,i7,a)') 'surface counts ',surface(is)%nt, ' triangles'
write(*,'(a,2f30.20)') 'surf%x range : ',minval(surface(is)%x),maxval(surface(is)%x)
write(*,'(a,2f30.20)') 'surf%y range : ',minval(surface(is)%y),maxval(surface(is)%y)
write(*,'(a,2f30.20)') 'surf%z range : ',minval(surface(is)%z),maxval(surface(is)%z)
end do

nnn=0
nnne=0
do i=1,nn
if (invoid(i).eq.0) nnn=nnn+1
enddo
do i=1,ov%nleaves
if (elvoid(i).eq.0) nnne=nnne+1
enddo

!==============================================================================
!=====[producing cubes.vtk]====================================================
!==============================================================================
if (output_cubes==1) then

   !=====[selection of leaves we are interested in]====
   allocate(subset_leaves(ov%nleaves))
   maxe2d=maxval(ov%e2d)
   subset_leaves=.false.
   do i=1,ov%nleaves
      xx=ov%x(ov%icon(1,i))
      yy=ov%y(ov%icon(1,i))
      zz=ov%z(ov%icon(1,i))
      if (.not.(xx.gt..5d0 .and. yy.gt. .5d0) .and.  ov%e2d(i)>0.06 *maxe2d) subset_leaves(i)=.true.
   end do
   !subset_leaves=ov%whole_leaf_in_fluid

   nl=count(subset_leaves)
   !np=8*nl
   np=8*ov%nleaves

   iunit=30
   open(unit=iunit,file='cubes.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 ',np,' float'

   do i=1,ov%nleaves
      !if (subset_leaves(i)) then
         do k=1,8 
           xx=ov%x(ov%icon(k,i))
           yy=ov%y(ov%icon(k,i))
           zz=ov%z(ov%icon(k,i))
           write(iunit,'(3e11.4)') xx,yy,zz
         end do
      !end if
   enddo

   write(iunit,'(A6, 2I10)') 'CELLS ',nl,(1+8)*nl
   do ie=1,ov%nleaves
      if (subset_leaves(ie)) then
         write(iunit,'(9I10)') 8 , ((ie-1)*8+i-1,i=1,8)
      end if
   end do

   write(iunit,'(A11, I10)') 'CELL_TYPES ',nl
   do ie=1,ov%nleaves
      if (subset_leaves(ie)) then
         write(iunit,'(I2)') 11
      end if 
   end do

   write(iunit,'(a11,i10)')'POINT_DATA ',np

   write(iunit,'(a)')'SCALARS e2d float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,ov%nleaves
         do k=1,8 
            write(iunit,'(e11.4)') ov%e2d(i)
         end do
   end do
   
   write(iunit,'(a)')'SCALARS w float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,ov%nleaves
      do k=1,8 
      write(iunit,'(e11.4)')sum(ov%w(ov%icon(:,i)))/8.d0
      end do
   end do
   write(iunit,'(a)')'SCALARS v float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,ov%nleaves
      do k=1,8 
      write(iunit,'(e11.4)')sum(ov%v(ov%icon(:,i)))/8.d0
      end do
   end do
   write(iunit,'(a)')'SCALARS u float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,ov%nleaves
      do k=1,8 
      write(iunit,'(e11.4)')sum(ov%u(ov%icon(:,i)))/8.d0
      end do
   end do

   write(iunit,'(a)')'SCALARS uvw float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,ov%nleaves
      do k=1,8 
      write(iunit,'(e11.4)') sqrt((sum(ov%u(ov%icon(:,i)))/8.d0)**2+(sum(ov%v(ov%icon(:,i)))/8.d0)**2+(sum(ov%w(ov%icon(:,i)))/8.d0)**2)
      end do
   end do

   write(*,*) '--------------------------------------------------------------------------'
   write(*,*) 'opla I am done producing cubes.vtk'
end if




!==============================================================================
!======[producing strain.vtk]==================================================
!==============================================================================

if (output_strain_tensor/=0) then

nnnep=0
levelmax=2
allocate (instrain(ov%nleaves))
instrain=.FALSE.
do i=1,ov%nleaves
   if (elvoid(i).eq.0) then
   xx=sum(ov%x(ov%icon(:,i)))/8.d0
   yy=sum(ov%y(ov%icon(:,i)))/8.d0
   zz=sum(ov%z(ov%icon(:,i)))/8.d0
   xx=xx*(2.d0**levelmax)
   yy=yy*(2.d0**levelmax)
   zz=zz*(2.d0**levelmax)
   ix=int(xx)
   iy=int(yy)
   iz=int(zz)
   if (abs(xx-float(ix)-0.5d0).lt.1.d-6 .and. &
       abs(yy-float(iy)-0.5d0).lt.1.d-6 .and. &
       abs(zz-float(iz)-0.5d0).lt.1.d-6) instrain(i)=.TRUE.
   if (instrain(i)) nnnep=nnnep+1
   endif
enddo 
!print*,'here',nnne,nnnep

iunit=30
open(unit=iunit,file='strain.vtk')
write(iunit,'(a)')'# vtk DataFile Version 3.0'
write(iunit,'(a)')'sigma1'
write(iunit,'(a)')'ASCII'
write(iunit,'(a)')'DATASET UNSTRUCTURED_GRID'
write(iunit,'(a7,i10,a6)')'POINTS ',nnnep,' float'

do i=1,ov%nleaves
   if (elvoid(i).eq.0) then
   xx=sum(ov%x(ov%icon(:,i)))/8.d0
   yy=sum(ov%y(ov%icon(:,i)))/8.d0
   zz=sum(ov%z(ov%icon(:,i)))/8.d0
   if (instrain(i)) write(iunit,'(3f16.11)') xx,yy,zz
   endif
enddo 

write(iunit,'(a11,i10)')'POINT_DATA ',nnnep

write(iunit,'(a)')'TENSORS strain float'
do i=1,ov%nleaves
   if (elvoid(i).eq.0) then
   xx=sum(ov%x(ov%icon(:,i)))/8.d0
   yy=sum(ov%y(ov%icon(:,i)))/8.d0
   zz=sum(ov%z(ov%icon(:,i)))/8.d0
   if (instrain(i)) then
   write(iunit,'(3e11.4)') strain(1,:,i)
   write(iunit,'(3e11.4)') strain(2,:,i)
   write(iunit,'(3e11.4)') strain(3,:,i)
   endif
   endif
enddo
close (iunit)

deallocate (instrain)

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

end if

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

if (output_ps/=0) then

allocate(do_it(ov%nleaves))

np=0
do i=1,ov%nleaves
   zz=sum(ov%z(ov%icon(:,i)))/8.d0  !criterion
   if (zz<.3d0) then
      np=np+1
      do_it(i)=.true.
   end if
end do

write(*,*) '--------------------------------------------------------------------------'
write(*,*) 'subset of points in compressive.vtk:',np

iunit=30
open(unit=iunit,file='compressive.vtk')
write(iunit,'(a)')'# vtk DataFile Version 3.0'
write(iunit,'(a)')'sigma1'
write(iunit,'(a)')'ASCII'
write(iunit,'(a)')'DATASET UNSTRUCTURED_GRID'
write(iunit,'(a7,i10,a6)')'POINTS ',np,' float'

do i=1,ov%nleaves
   if (do_it(i)) then
      xx=sum(ov%x(ov%icon(:,i)))/8.d0
      yy=sum(ov%y(ov%icon(:,i)))/8.d0
      zz=sum(ov%z(ov%icon(:,i)))/8.d0
      write(iunit,'(3f16.11)') xx,yy,zz
   end if
enddo

write(iunit,'(a11,i10)')'POINT_DATA ', np

write(iunit,'(a)')'VECTORS sigma1head+ float'
do i=1,ov%nleaves
   if (do_it(i)) & 
   write(iunit,'(3e15.4)') n1(1:3,i)*s(1,i)
enddo
write(iunit,'(a)')'VECTORS sigma1head- float'
do i=1,ov%nleaves
   if (do_it(i)) & 
   write(iunit,'(3e15.4)') -n1(1:3,i)*s(1,i)
enddo

write(iunit,'(a)')'VECTORS sigma3tail+ float'
do i=1,ov%nleaves
   if (do_it(i)) & 
   write(iunit,'(3e15.4)') n2(1:3,i)*s(2,i)
enddo

write(iunit,'(a)')'VECTORS sigma3tail- float'
do i=1,ov%nleaves
   if (do_it(i)) & 
   write(iunit,'(3e15.4)') -n2(1:3,i)*s(2,i)
enddo

write(iunit,'(a)')'VECTORS sigma2head+ float'
do i=1,ov%nleaves
   if (do_it(i)) & 
   write(iunit,'(3e15.4)') n3(1:3,i)*min(0.d0,s(3,i))
enddo

write(iunit,'(a)')'VECTORS sigma2head- float'
do i=1,ov%nleaves
   if (do_it(i)) & 
   write(iunit,'(3e15.4)') -n3(1:3,i)*min(0.d0,s(3,i))
enddo

write(iunit,'(a)')'VECTORS sigma2tail+ float'
do i=1,ov%nleaves
   if (do_it(i)) & 
   write(iunit,'(3e15.4)') n3(1:3,i)*max(0.d0,s(3,i))
enddo

write(iunit,'(a)')'VECTORS sigma2tail- float'
do i=1,ov%nleaves
   if (do_it(i)) & 
   write(iunit,'(3e15.4)') -n3(1:3,i)*max(0.d0,s(3,i))
enddo

close (iunit)

deallocate(do_it)

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

end if

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


!iunit=30
!open(unit=iunit,file='cut.dat')
!do i=1,nnn              ! loop over nodes
!   x=ov%x(rtf(i))
!   y=ov%y(rtf(i))
!   z=ov%z(rtf(i))
!   icut=3   ! icut between 0 and 64 for uniform level 6
!   ! example , cross section in plane given by x=value
!   xcut=icut*2.d0**(-6)    ! background uniform grid level is 6
!   dxyz=2.d0**(-7)      ! higher level grid in boxes is 7
!   if (abs(x-xcut)<1.d-6 ) then   ! select nodes that are only on the plane
!      if (mod(nint(y/dxyz),2)==0) then  ! les coordonnees des noeuds surla grille de niveau 7 divisees par dxy vont prendre les valeurs 0 a 128
!      if (mod(nint(z/dxyz),2)==0) then
!         write(iunit,*) x,y,z
!      end if
!      end if
!   end if
!enddo
!close(iunit)


iunit=30
open(unit=iunit,file='visu.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 ',nnn,' float'

do i=1,nnn
write(iunit,'(3f16.11)')ov%x(rtf(i)),ov%y(rtf(i)),ov%z(rtf(i))
enddo

write(iunit,'(A6, 2I10)') 'CELLS ',nnne,(1+nne)*nnne
iconmin=nn
iconmax=0
do ie=1,ov%nleaves
   if (elvoid(ie).eq.0) then
   myicon(1:nne)=ftr(ov%icon(1:nne,ie))-1
     do k=1,nne
     iconmin=min(iconmin,myicon(k))
     iconmax=max(iconmax,myicon(k))
     enddo
   write(iunit,'(9I10)')nne,myicon(1:nne)
   endif
enddo
!print*,iconmin,iconmax

write(iunit,'(A11, I10)') 'CELL_TYPES ',nnne
do ie=1,nnne
   write(iunit,'(I2)')11 ! octree  (8 nodes)
end do

write(iunit,'(a11,i10)')'POINT_DATA ',nnn

if(output_u_field==1) then 
   write(iunit,'(a)')'SCALARS u float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      write(iunit,'(e11.4)') ov%u(rtf(i))
   enddo
end if

if(output_v_field==1) then 
   write(iunit,'(a)')'SCALARS v float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      write(iunit,'(e11.4)') ov%v(rtf(i))
   enddo
end if

if(output_w_field==1) then 
   write(iunit,'(a)')'SCALARS w float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      write(iunit,'(e11.4)') ov%w(rtf(i))
   enddo
end if


if(output_disp_field==1) then 
   write(iunit,'(a)')'SCALARS disp float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      write(iunit,'(e11.4)') sqrt(ov%u(rtf(i))**2+ov%v(rtf(i))**2+ov%w(rtf(i))**2)
   enddo
end if

if(output_press_fieldn==1) then 
   write(iunit,'(a)')'SCALARS pressure_n float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      write(iunit,'(e11.4)') ov%nodal_pressure(rtf(i))
   enddo
end if

if(output_countnode_field==1) then 
   write(iunit,'(a)')'SCALARS countnode float 1' 
   write(iunit,'(a)')'LOOKUP_TABLE default' 
   do i=1,nnn 
      write(iunit,'(e11.4)') countnode(rtf(i)) 
   enddo 
end if 
 
if(output_e2d_fieldn==1) then 
   write(iunit,'(a)')'SCALARS e2dn float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      write(iunit,'(e11.4)') ov_nodee2d(rtf(i))
   enddo
end if

if(output_crit_field==1) then 
   write(iunit,'(a)')'SCALARS crit float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      write(iunit,'(e11.4)') ov_nodecrit(rtf(i))
   enddo
end if

if(output_strain_field==1) then 
   write(iunit,'(a)')'SCALARS s float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      write(iunit,'(e11.4)') ov%strain(rtf(i))
   enddo
end if

if(output_temp_field==1) then 
   write(iunit,'(a)')'SCALARS t float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      write(iunit,'(e11.4)') 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//' float 1'
      write(iunit,'(a)')'LOOKUP_TABLE default'
      do i=1,nnn
         write(iunit,'(e11.4)') ov%lsf(rtf(i),ilsf)
      enddo
   enddo
end if

if(output_velo_vect==1) then 
   write(iunit,'(a)')'VECTORS velo float' 
   do i=1,nnn 
      write(iunit,'(3e15.4)')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 float' 
   do i=1,nnn 
      write(iunit,'(3e15.4)')ov%u(rtf(i)),ov%v(rtf(i)),ov%w(rtf(i))-ov%wiso(rtf(i)) 
   enddo 
endif

if(output_iso_velo_vect==1) then 
   write(iunit,'(a)')'VECTORS iso-velo float' 
   do i=1,nnn 
      write(iunit,'(3e15.4)')0.d0*ov%u(rtf(i)),0.d0*ov%v(rtf(i)),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 float'
   do i=1,nnn
      write(iunit,'(3e11.4)') strainn(1,:,i)
      write(iunit,'(3e11.4)') strainn(2,:,i)
      write(iunit,'(3e11.4)') 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 float 1' 
   write(iunit,'(a)')'LOOKUP_TABLE default' 
   do ie=1,ov%nleaves 
      if (elvoid(ie).eq.0) then 
         write(iunit,'(e11.4)') ov%e2d(ie)
      endif 
   enddo 
end if 

if(output_press_fielde==1) then  
   write(iunit,'(a)')'SCALARS pressure_e float 1' 
   write(iunit,'(a)')'LOOKUP_TABLE default' 
   do ie=1,ov%nleaves 
      if (elvoid(ie).eq.0) then 
         write(iunit,'(e11.4)') ov%pressure(ie) 
      endif 
   enddo 
end if 

if(output_smooth_press_fielde==1) then  
   write(iunit,'(a)')'SCALARS smooth_pressure_e float 1' 
   write(iunit,'(a)')'LOOKUP_TABLE default' 
   do ie=1,ov%nleaves 
      if (elvoid(ie).eq.0) then 
         write(iunit,'(e11.4)') ov%spressure(ie)
      endif 
   enddo 
end if

if(output_eviscosity_fielde==1) then  
   write(iunit,'(a)')'SCALARS eff_visc_e float 1' 
   write(iunit,'(a)')'LOOKUP_TABLE default' 
   do ie=1,ov%nleaves 
      if (elvoid(ie).eq.0) then 
         write(iunit,'(e11.4)') 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_dilatr_fielde==1) then
   write(iunit,'(a)')'SCALARS dilatation_rate float 1' 
   write(iunit,'(a)')'LOOKUP_TABLE default' 
   do ie=1,ov%nleaves 
      if (elvoid(ie).eq.0) then 
         write(iunit,'(e11.4)') 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,'(e11.4)') 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,'(e11.4)') -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,'(e11.4)') surface(is)%w(i)
      enddo
endif
      write(30,'(a)') 'VECTORS normal float'
      do i=1,surface(is)%nsurface
         write(30,'(3e15.4)') 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,'(3e15.4)') 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,'(e11.4)') 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,'(e11.4)') 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,'(e11.4)') 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,'(e11.4)') 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,'(e11.4)') 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,'(e11.4)') dip3(i)
!        enddo
      write(30,'(a)')'TENSORS StrainDot float'
        do i=1,(2**levmax+1)*(2**levmax+1)*nz
        write(30,'(3e12.4)') s11(i),s12(i),s13(i)
        write(30,'(3e12.4)') s12(i),s22(i),s23(i)
        write(30,'(3e12.4)') 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,'(e11.4)') 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,'(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
