From 52d2b5dc9d23eb0a2e77d932885959bc05950d00 Mon Sep 17 00:00:00 2001
From: holger <holger-grosshans@gmx.de>
Date: Wed, 13 May 2020 10:44:39 +0200
Subject: [PATCH] partial merge japan

	modified:   input/input_example.dat
	modified:   src/bc.f90
	modified:   src/electrostatics.f90
	modified:   src/makefile
	modified:   src/particles.f90
	modified:   src/post.f90
	new file:   src/writedat_fluid_xy.f90
	new file:   src/writedat_fluid_xz.f90
	new file:   src/writevtk_fluid_xy.f90
	modified:   src/writevtk_fluid_xyz.f90
	modified:   src/writevtk_fluid_xz.f90
	modified:   src/writevtk_fluid_yz.f90
---
 input/input_example.dat    |   2 +-
 src/bc.f90                 |   9 +--
 src/electrostatics.f90     |  33 ++++-----
 src/makefile               |  17 +++--
 src/particles.f90          |   2 +-
 src/post.f90               |  10 ++-
 src/writedat_fluid_xy.f90  |  93 ++++++++++++++++++++++++
 src/writedat_fluid_xz.f90  |  93 ++++++++++++++++++++++++
 src/writevtk_fluid_xy.f90  | 144 +++++++++++++++++++++++++++++++++++++
 src/writevtk_fluid_xyz.f90 |   2 +-
 src/writevtk_fluid_xz.f90  |   2 +-
 src/writevtk_fluid_yz.f90  |   5 +-
 12 files changed, 376 insertions(+), 36 deletions(-)
 create mode 100644 src/writedat_fluid_xy.f90
 create mode 100644 src/writedat_fluid_xz.f90
 create mode 100644 src/writevtk_fluid_xy.f90

diff --git a/input/input_example.dat b/input/input_example.dat
index f8e75a6..0aa6a27 100644
--- a/input/input_example.dat
+++ b/input/input_example.dat
@@ -14,7 +14,7 @@ interval of writing result files, restart files (-,-) =
 10,10000
 CFL number (-) =
 0.4
-initial bulk flow (m/s), pressure gradient (N/m**3) =
+in/initial flow (m/s), pressure gradient (N/m**3) =
 3.65,-2.9
 fluid density (kg/m**3), fluid kinematic viscosity (m**2/s) =
 1.2,1.46E-5
diff --git a/src/bc.f90 b/src/bc.f90
index 091ea2e..a27ae75 100644
--- a/src/bc.f90
+++ b/src/bc.f90
@@ -50,7 +50,7 @@
             if (bcy.eq.'w'.and.bcz.eq.'p') dist=delta-abs(yc(j))
             if (bcy.eq.'p'.and.bcz.eq.'w') dist=delta-abs(zc(l))
             do m=1,gc
-              u(imin-m,j,l)=1._pr*ubulk0*dist/delta + (random(mm)-.5_pr)*alpha*ubulk
+              u(imin-m,j,l)=ubulk0*dist/delta + (random(mm)-.5_pr)*alpha*ubulk
               v(imin-m,:,:)=0._pr
               w(imin-m,:,:)=0._pr
             enddo
@@ -58,9 +58,10 @@
         endif
         if (myid.eq.nrprocs-1) then
            do 8 m=1,gc
-           u(imax+m,:,:)  =u(imax,:,:)
-           v(imax+m,:,:)  =v(imax,:,:)
-           w(imax+m,:,:)  =w(imax,:,:)
+           u(imax+m-1,:,:)=u(imax-1,:,:)
+           v(imax+m,:,:)=  v(imax,:,:)
+           w(imax+m,:,:)=  w(imax,:,:)
+           p(imax+m,:,:)=  p(imax,:,:)
 8          enddo
         endif
       endif
diff --git a/src/electrostatics.f90 b/src/electrostatics.f90
index 47a3fb3..7b67d3c 100644
--- a/src/electrostatics.f90
+++ b/src/electrostatics.f90
@@ -292,7 +292,10 @@
 !        ut2=up(n)**2
 !        un2=vp(n)**2+wp(n)**2
 
-      if (direction.eq.2) then
+      if (direction.eq.1) then
+        ut2=vp(n)**2+wp(n)**2
+        un2=up(n)**2
+      elseif (direction.eq.2) then
         ut2=up(n)**2+wp(n)**2
         un2=vp(n)**2
       elseif (direction.eq.3) then
@@ -301,9 +304,10 @@
       endif
 
 ! John (1979)
-      alpha1=(0.625_pr*pi*rhop*(1._pr+restRatio)*(un2+ut2)* &
+      alpha1=(0.625_pr*pi*rhop*(1._pr+restRatio)*un2* &
             ((1._pr-nyw**2)/Ew+(1._pr-nyp**2)/Ep))**(0.4_pr)
       AoAtot=alpha1/4._pr
+
       dqp=Qaccfactor*AoAtot*(qpmax-q_el(n))
       q_el(n)=min(dqp+q_el(n),qpmax)
       
@@ -323,26 +327,23 @@
 !> @brief compute particle-particle charging
       subroutine chargeParticleParticle(n1,n2)
       use var
-      real(kind=pr) uax2,urad2,alpha1,dqp
+      real(kind=pr) urel2,alpha1,dqp
       integer :: n1,n2
       character(70) :: filename
 
-      uax2=(up(n1)-up(n2))**2
-      urad2=(vp(n1)-vp(n2))**2+(wp(n1)-wp(n2))**2
-
-      if(urad2.eq.0.) then
-        dqp=0.5*(q_el(n2)-q_el(n1))
-        goto 10
-      endif
+      urel2=(up(n1)-up(n2))**2+(vp(n1)-vp(n2))**2+(wp(n1)-wp(n2))**2
 
+      if (urel2.eq.0._pr) then
+        dqp=0.5_pr*(q_el(n2)-q_el(n1))
+      else
 ! Soo (1971)
-      alpha1=radp(n1)*radp(n2)*(0.625*pi*rhop*(1+restratio)*(uax2+urad2) &
-             *(radp(n1)+radp(n2))**(0.5)/(radp(n1)**3+radp(n2)**3) &
-             *(1-nyp**(2.0))/Ep)**(0.4)
-
-      dqp=alpha1/8._pr/radp(n1)*(q_el(n2)-q_el(n1))
+        alpha1=radp(n1)*radp(n2)*(0.625_pr*pi*rhop*(1._pr+restRatio)*urel2 &
+               *(radp(n1)+radp(n2))**(0.5_pr)/(radp(n1)**3+radp(n2)**3) &
+               *(1._pr-nyp**2)/Ep)**(0.4_pr)
+        dqp=alpha1/8._pr/radp(n1)*(q_el(n2)-q_el(n1))
+      endif
 
-10    q_el(n1)=q_el(n1)+dqp
+      q_el(n1)=q_el(n1)+dqp
       q_el(n2)=q_el(n2)-dqp
 
 
diff --git a/src/makefile b/src/makefile
index 148b351..f0288b7 100644
--- a/src/makefile
+++ b/src/makefile
@@ -3,21 +3,22 @@
 OBJ =	var.o main.o pre.o bc.o restart.o fluid.o \
       momentum.o post.o schemes.o \
       pressure.o writevtk_fluid_xz.o \
-      writevtk_fluid_yz.o writedat_fluid.o \
       parallel.o writevtk_grid.o particles.o \
       mom1.o mom5.o particlesTransport.o \
       mass.o electrostatics.o \
       writevtk_particles.o writedat_particles.o \
-	writevtk_fluid_xyz.o
+ 	    writevtk_fluid_xy.o writevtk_fluid_xyz.o \
+      writevtk_fluid_yz.o writedat_fluid_xz.o \
+      writedat_fluid_xy.o
 INC =
 
 
 # gcc:
 CMP = mpifort
 #O3: optimization
-FLAGS   = -O3 -mcmodel=medium -pg
+#FLAGS   = -O3 -mcmodel=medium
 #debugging:
-#FLAGS   = -g3 -O0 -fbounds-check -mcmodel=medium -fimplicit-none -fcheck=all -fbacktrace -floop-nest-optimize -ffpe-trap=invalid,zero,overflow -Wconversion -fno-tree-vectorize
+FLAGS   = -g3 -O0 -fbounds-check -mcmodel=medium -fimplicit-none -fcheck=all -fbacktrace -floop-nest-optimize -ffpe-trap=invalid,zero,overflow -Wconversion -fno-tree-vectorize
 FFLAGS  = -c $(FLAGS)
 
 
@@ -43,14 +44,18 @@ pressure.o: pressure.f90 $(INC)
 	$(CMP) $(FFLAGS) pressure.f90
 mass.o: mass.f90 $(INC)
 	$(CMP) $(FFLAGS) mass.f90
+writevtk_fluid_xy.o: writevtk_fluid_xy.f90 $(INC)
+	$(CMP) $(FFLAGS) writevtk_fluid_xy.f90
 writevtk_fluid_xz.o: writevtk_fluid_xz.f90 $(INC)
 	$(CMP) $(FFLAGS) writevtk_fluid_xz.f90
 writevtk_fluid_yz.o: writevtk_fluid_yz.f90 $(INC)
 	$(CMP) $(FFLAGS) writevtk_fluid_yz.f90
 writevtk_fluid_xyz.o: writevtk_fluid_xyz.f90 $(INC)
 	$(CMP) $(FFLAGS) writevtk_fluid_xyz.f90
-writedat_fluid.o: writedat_fluid.f90 $(INC)
-	$(CMP) $(FFLAGS) writedat_fluid.f90
+writedat_fluid_xy.o: writedat_fluid_xy.f90 $(INC)
+	$(CMP) $(FFLAGS) writedat_fluid_xy.f90
+writedat_fluid_xz.o: writedat_fluid_xz.f90 $(INC)
+	$(CMP) $(FFLAGS) writedat_fluid_xz.f90
 writevtk_particles.o: writevtk_particles.f90 $(INC)
 	$(CMP) $(FFLAGS) writevtk_particles.f90
 writedat_particles.o: writedat_particles.f90 $(INC)
diff --git a/src/particles.f90 b/src/particles.f90
index 7b1586d..82ffade 100644
--- a/src/particles.f90
+++ b/src/particles.f90
@@ -13,7 +13,7 @@
         if (elForceScheme.ne.2) call forcesGauss
         if (elForceScheme.ne.1) call forcesCoulomb
         call particlesVelocity
-        call particlesCollide
+!        call particlesCollide
         call particlesTransport    !particles from neighbour gc added
         call chargeDensity
         call momentumCoupling !particles from neighbour gc removed
diff --git a/src/post.f90 b/src/post.f90
index 88746c3..515f085 100644
--- a/src/post.f90
+++ b/src/post.f90
@@ -6,11 +6,15 @@
       if (mod(nt,int_results).eq.0.or.(nt).eq.1) then
 
       call monitor
+      call writevtk_fluid_xy
       call writevtk_fluid_xz
       call writevtk_fluid_yz
-      call writedat_ufluid_xz
-      call writedat_vfluid_xz
-      call writedat_wfluid_xz
+      call writedat_ufluid_xy
+      call writedat_vfluid_xy
+      call writedat_wfluid_xy
+!      call writedat_ufluid_xz
+!      call writedat_vfluid_xz
+!      call writedat_wfluid_xz
 !      call writedat_coulomb
 !      call writevtk_fluid_xyz
       call writevtk_particles
diff --git a/src/writedat_fluid_xy.f90 b/src/writedat_fluid_xy.f90
new file mode 100644
index 0000000..a3c0e89
--- /dev/null
+++ b/src/writedat_fluid_xy.f90
@@ -0,0 +1,93 @@
+!####################################################################
+!> @author Holger Grosshans
+      subroutine writedat_ufluid_xy
+      use var
+      character(70) :: filename,rowfmt
+      integer :: i,j,l
+
+100   format(es11.4e2)
+      write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))'
+
+      l=int(ll/2)
+
+      write(filename,'(a,i3.3,a,i6.6,a)') &
+        'results/fluid_u_xy_p',myid,'_',nt,'.dat'
+      open(10,file=filename)
+
+      write(10,'(a,es14.6e2,a,x,a)') '# u: t=',t,' s --',version
+
+      write(10,'(a11)',advance='no') '#       y/x'
+      write(10,fmt=rowfmt) (xf(i),i=imin,imax+1)
+      do j=jmin-1,jmax+1
+        write(10,fmt=100,advance='no') yc(j) 
+        write(10,fmt=rowfmt) (u(i,j,l),i=imin,imax+1)
+      enddo
+ 
+      close(10)
+
+      return
+
+      end
+
+!####################################################################
+!> @author Holger Grosshans
+      subroutine writedat_vfluid_xy
+      use var
+      character(70) :: filename,rowfmt
+      integer :: i,j,l
+
+100   format(es11.4e2)
+      write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))'
+
+      l=int(ll/2)
+
+      write(filename,'(a,i3.3,a,i6.6,a)') &
+        'results/fluid_v_xy_p',myid,'_',nt,'.dat'
+      open(10,file=filename)
+
+      write(10,'(a,es14.6e2,a,x,a)') '# v: t=',t,' s --',version
+
+      write(10,'(a11)',advance='no') '#       y/x'
+      write(10,fmt=rowfmt) (xc(i),i=imin,imax+1)
+      do j=jmin-1,jmax+1
+        write(10,fmt=100,advance='no') yf(j) 
+        write(10,fmt=rowfmt) (v(i,j,l),i=imin,imax+1)
+      enddo
+ 
+      close(10)
+
+      return
+
+      end
+
+!####################################################################
+!> @author Holger Grosshans
+      subroutine writedat_wfluid_xy
+      use var
+      character(70) :: filename,rowfmt
+      integer :: i,j,l
+
+100   format(es11.4e2)
+      write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))'
+
+      l=int(ll/2)
+
+      write(filename,'(a,i3.3,a,i6.6,a)') &
+        'results/fluid_w_xy_p',myid,'_',nt,'.dat'
+      open(10,file=filename)
+
+      write(10,'(a,es14.6e2,a,x,a)') '# w: t=',t,' s --',version
+
+      write(10,'(a11)',advance='no') '#       y/x'
+      write(10,fmt=rowfmt) (xc(i),i=imin,imax+1)
+      do j=jmin-1,jmax+1
+        write(10,fmt=100,advance='no') yc(j) 
+        write(10,fmt=rowfmt) (w(i,j,l),i=imin,imax+1)
+      enddo
+ 
+      close(10)
+
+      return
+
+      end
+
diff --git a/src/writedat_fluid_xz.f90 b/src/writedat_fluid_xz.f90
new file mode 100644
index 0000000..7798953
--- /dev/null
+++ b/src/writedat_fluid_xz.f90
@@ -0,0 +1,93 @@
+!####################################################################
+!> @author Holger Grosshans
+      subroutine writedat_ufluid_xz
+      use var
+      character(70) :: filename,rowfmt
+      integer :: i,j,l
+
+100   format(es11.4e2)
+      write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))'
+
+      j=int(jj/2)
+
+      write(filename,'(a,i3.3,a,i6.6,a)') &
+        'results/fluid_u_xz_p',myid,'_',nt,'.dat'
+      open(10,file=filename)
+
+      write(10,'(a,es14.6e2,a,x,a)') '# u: t=',t,' s --',version
+
+      write(10,'(a11)',advance='no') '#       z/x'
+      write(10,fmt=rowfmt) (xf(i),i=imin,imax+1)
+      do l=lmin-1,lmax+1
+        write(10,fmt=100,advance='no') zc(l) 
+        write(10,fmt=rowfmt) (u(i,j,l),i=imin,imax+1)
+      enddo
+ 
+      close(10)
+
+      return
+
+      end
+
+!####################################################################
+!> @author Holger Grosshans
+      subroutine writedat_vfluid_xz
+      use var
+      character(70) :: filename,rowfmt
+      integer :: i,j,l
+
+100   format(es11.4e2)
+      write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))'
+
+      j=int(jj/2)
+
+      write(filename,'(a,i3.3,a,i6.6,a)') &
+        'results/fluid_v_xz_p',myid,'_',nt,'.dat'
+      open(10,file=filename)
+
+      write(10,'(a,es14.6e2,a,x,a)') '# v: t=',t,' s --',version
+
+      write(10,'(a11)',advance='no') '#       z/x'
+      write(10,fmt=rowfmt) (xc(i),i=imin,imax+1)
+      do l=lmin-1,lmax+1
+        write(10,fmt=100,advance='no') zc(l) 
+        write(10,fmt=rowfmt) (v(i,j,l),i=imin,imax+1)
+      enddo
+ 
+      close(10)
+
+      return
+
+      end
+
+!####################################################################
+!> @author Holger Grosshans
+      subroutine writedat_wfluid_xz
+      use var
+      character(70) :: filename,rowfmt
+      integer :: i,j,l
+
+100   format(es11.4e2)
+      write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))'
+
+      j=int(jj/2)
+
+      write(filename,'(a,i3.3,a,i6.6,a)') &
+        'results/fluid_w_xz_p',myid,'_',nt,'.dat'
+      open(10,file=filename)
+
+      write(10,'(a,es14.6e2,a,x,a)') '# w: t=',t,' s --',version
+
+      write(10,'(a11)',advance='no') '#       z/x'
+      write(10,fmt=rowfmt) (xc(i),i=imin,imax+1)
+      do l=lmin-1,lmax+1
+        write(10,fmt=100,advance='no') zf(l) 
+        write(10,fmt=rowfmt) (w(i,j,l),i=imin,imax+1)
+      enddo
+ 
+      close(10)
+
+      return
+
+      end
+
diff --git a/src/writevtk_fluid_xy.f90 b/src/writevtk_fluid_xy.f90
new file mode 100644
index 0000000..3ef111e
--- /dev/null
+++ b/src/writevtk_fluid_xy.f90
@@ -0,0 +1,144 @@
+!> @author Holger Grosshans
+      subroutine writevtk_fluid_xy
+      use var
+      integer :: i,j,l,m
+      character(70) :: filename,filename2,rowfmt,vecfmt
+      logical :: file_ex
+
+
+100   format(es11.4e2)
+      write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))'
+      write(vecfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,3(es11.4e2)))'
+
+      l=int(ll/2.)+1
+
+      write(filename,'(a,i3.3,a,i6.6,a)') 'results/fluid_xy_p',myid,'_',nt,'.vtk'
+      open(10,file=filename)
+
+! write visit container file
+      filename2='results/fluid_xy.visit'
+      if(myid.eq.0) then
+        inquire(file=filename2,exist=file_ex)
+        if (file_ex.and.nt.gt.1) then
+          open(11,file=filename2,access='append')
+        else
+          open(11,file=filename2)
+          write(11,'(a8,x,i3)') '!NBLOCKS',nrprocs
+        endif
+        do m=1,nrprocs
+          write(11,'(a,i3.3,a,i6.6,a)') 'fluid_xy_p',(m-1),'_',nt,'.vtk'
+        enddo
+        close(11)
+      endif
+
+
+      write(10,'(a)') '# vtk DataFile Version 3.0'
+      write(10,'(a,es14.6e2,a,x,a)') 'fluid data: t=',t,' s --',version
+      write(10,'(a)') 'ASCII'
+      write(10,'(a)') 'DATASET RECTILINEAR_GRID'
+      write(10,'(a10,i4,i4,a4)') &
+            'DIMENSIONS',(imax-imin+2),(jmax-jmin+3),'1'
+
+      write(10,*) 'X_COORDINATES',(imax-imin+2),'float'
+      write(10,fmt=rowfmt) (xc(i),i=imin,imax+1)
+
+      write(10,*) 'Y_COORDINATES',(jmax-jmin+3),'float'
+      write(10,fmt=rowfmt) (yc(j),j=jmin-1,jmax+1)
+
+      write(10,*) 'Z_COORDINATES 1 float'
+      write(10,100) zc(l)
+
+      write(10,*) 'POINT_DATA ',(imax-imin+2)*(jmax-jmin+3)
+
+      write(10,'(a)') 'SCALARS u float 1'
+      write(10,'(a)') 'LOOKUP_TABLE default '
+      do j=jmin-1,jmax+1
+        write(10,fmt=rowfmt) (0.5_pr*(u(i,j,l)+u(max(i-1,1),j,l)),i=imin,imax+1)
+      enddo
+
+      write(10,'(a)') 'SCALARS v float 1'
+      write(10,'(a)') 'LOOKUP_TABLE default '
+      do j=jmin-1,jmax+1
+        write(10,fmt=rowfmt) (0.5_pr*(v(i,j,l)+v(i,max(j-1,1),l)),i=imin,imax+1)
+      enddo
+
+      write(10,'(a)') 'SCALARS w float 1'
+      write(10,'(a)') 'LOOKUP_TABLE default '
+      do j=jmin-1,jmax+1
+        write(10,fmt=rowfmt) (0.5_pr*(w(i,j,l)+w(i,j,max(l-1,1))),i=imin,imax+1)
+      enddo
+
+!      write(10,'(a)') 'SCALARS uraw float 1'
+!      write(10,'(a)') 'LOOKUP_TABLE default '
+!      do j=jmin-1,jmax+1
+!        write(10,fmt=rowfmt) (u(i,j,l),i=imin,imax+1)
+!      enddo
+!
+!      write(10,'(a)') 'SCALARS vraw float 1'
+!      write(10,'(a)') 'LOOKUP_TABLE default '
+!      do j=jmin-1,jmax+1
+!        write(10,fmt=rowfmt) (v(i,j,l),i=imin,imax+1)
+!      enddo
+!
+!      write(10,'(a)') 'SCALARS wraw float 1'
+!      write(10,'(a)') 'LOOKUP_TABLE default '
+!      do j=jmin-1,jmax+1
+!        write(10,fmt=rowfmt) (w(i,j,l),i=imin,imax+1)
+!      enddo
+
+      write(10,'(a)') 'SCALARS p float 1'
+      write(10,'(a)') 'LOOKUP_TABLE default '
+      do j=jmin-1,jmax+1
+        write(10,fmt=rowfmt) (p(i,j,l),i=imin,imax+1)
+      enddo
+
+      write(10,'(a)') 'SCALARS rho_el float 1'
+      write(10,'(a)') 'LOOKUP_TABLE default '
+      do j=jmin-1,jmax+1
+        write(10,fmt=rowfmt) (rho_el(i,j,l),i=imin,imax+1)
+      enddo
+
+      write(10,'(a)') 'SCALARS phi_el float 1'
+      write(10,'(a)') 'LOOKUP_TABLE default '
+      do j=jmin-1,jmax+1
+        write(10,fmt=rowfmt) (phi_el(i,j,l),i=imin,imax+1)
+      enddo
+
+      write(10,'(a)') 'SCALARS Ex_el float 1'
+      write(10,'(a)') 'LOOKUP_TABLE default '
+      do j=jmin-1,jmax+1
+        write(10,fmt=rowfmt) (Ex_el(i,j,l),i=imin,imax+1)
+      enddo
+
+      write(10,'(a)') 'SCALARS Ey_el float 1'
+      write(10,'(a)') 'LOOKUP_TABLE default '
+      do j=jmin-1,jmax+1
+        write(10,fmt=rowfmt) (Ey_el(i,j,l),i=imin,imax+1)
+      enddo
+
+      write(10,'(a)') 'SCALARS Ez_el float 1'
+      write(10,'(a)') 'LOOKUP_TABLE default '
+      do j=jmin-1,jmax+1
+        write(10,fmt=rowfmt) (Ez_el(i,j,l),i=imin,imax+1)
+      enddo
+
+!      write(10,'(a)') 'SCALARS Emag_el float 1'
+!      write(10,'(a)') 'LOOKUP_TABLE default '
+!      do j=jmin-1,jmax+1
+!        write(10,fmt=rowfmt) &
+!        (sqrt(Ex_el(i,j,l)**2+Ey_el(i,j,l)**2+Ez_el(i,j,l)**2)**(0.5),i=imin,imax+1)
+!      enddo
+!
+!
+!      write(10,'(a)') 'VECTORS uvw float'
+!      do j=jmin-1,jmax+1
+!        write(10,fmt='(3(es12.2e2))') (0.5_pr*(u(i,j,l)+u(max(i-1,1),j,l)), &
+!                              0.5_pr*(v(i,j,l)+v(i,max(j-1,1),l)), &
+!                              0.5_pr*(w(i,j,l)+w(i,j,max(l-1,1))),i=imin,imax+1)
+!      enddo
+
+
+      close(10)
+
+      return
+      end
diff --git a/src/writevtk_fluid_xyz.f90 b/src/writevtk_fluid_xyz.f90
index 00ee61a..22a9429 100644
--- a/src/writevtk_fluid_xyz.f90
+++ b/src/writevtk_fluid_xyz.f90
@@ -20,7 +20,7 @@
       filename2='results/fluid_xyz.visit'
       if(myid.eq.0) then
         inquire(file=filename2,exist=file_ex)
-        if (file_ex.and.nt.ne.0) then
+        if (file_ex.and.nt.gt.1) then
           open(11,file=filename2,access='append')
         else
           open(11,file=filename2)
diff --git a/src/writevtk_fluid_xz.f90 b/src/writevtk_fluid_xz.f90
index 9c81d44..832508e 100644
--- a/src/writevtk_fluid_xz.f90
+++ b/src/writevtk_fluid_xz.f90
@@ -21,7 +21,7 @@
       filename2='results/fluid_xz.visit'
       if(myid.eq.0) then
         inquire(file=filename2,exist=file_ex)
-        if (file_ex.and.nt.ne.0) then
+        if (file_ex.and.nt.gt.1) then
           open(11,file=filename2,access='append')
         else
           open(11,file=filename2)
diff --git a/src/writevtk_fluid_yz.f90 b/src/writevtk_fluid_yz.f90
index 50052ae..59f40d5 100644
--- a/src/writevtk_fluid_yz.f90
+++ b/src/writevtk_fluid_yz.f90
@@ -2,7 +2,7 @@
       subroutine writevtk_fluid_yz
       use var
       real(kind=pr), dimension(jj) :: temp4
-      real(kind=pr) :: temp,temp1,temp2,temp3,pp
+      real(kind=pr) :: temp,temp1,temp2,temp3
       integer :: i,j,l,m
       character(70) :: filename,filename2,rowfmt
       logical :: file_ex
@@ -12,7 +12,6 @@
       write(rowfmt,'(a,i4,a)') '(',(jmax-jmin+3),'(1x,es11.4e2))'
 
       i=int(ii/2.)+1
-      pp=0._pr
 
       write(filename,'(a,i3.3,a,i6.6,a)') 'results/fluid_yz_p',myid,'_',nt,'.vtk'
       open(10,file=filename)
@@ -21,7 +20,7 @@
       filename2='results/fluid_yz.visit'
       if(myid.eq.0) then
         inquire(file=filename2,exist=file_ex)
-        if (file_ex.and.nt.ne.1) then
+        if (file_ex.and.nt.gt.1) then
           open(11,file=filename2,access='append')
         else
           open(11,file=filename2)
-- 
GitLab