C _____________________________ Density and velocity fields on a mesh
C		                     

       INCLUDE 'PMparameters.h'
       INCLUDE 'PMfields.h'
       REAL     INPUT
C...................................................................
C			Read data and open files
      OPEN(17,FILE='Fields.DAT',STATUS='UNKNOWN')
      OPEN(19,FILE='Fields.small.DAT',STATUS='UNKNOWN')
      CALL RDTAPE

      write (*,*) ' RDTAPE is done'
      WRITE (17,100) HEADER,
     +                  AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,wspecies(1),
     +                  EKIN,
     +                  NROWC,NGRID,NRECL,Om0,Oml0,hubble,
     +                  Ocurv,lspecies(1),Nspecies
      WRITE (19,100) HEADER,
     +                  AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,wspecies(1),
     +                  EKIN,
     +                  NROWC,NGRID,NRECL,Om0,Oml0,hubble,
     +                  Ocurv,lspecies(1),Nspecies
100   FORMAT(1X,'Header=>',A45,/
     +            1X,' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/
     +            1X,' I =',I4,' WEIGHT=',F8.3,' Ekin=',E12.3,/
     +            1X,' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I6,/
     +            1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/
     +            1x,' Omega_curvature=',F7.3,/
     +            1x,' Number of part.1specie=',i8,' Nspecies=',i3)
      Box =INPUT(' Enter box size in Mpc/h   =')
      alpha =INPUT(' Enter angle of rotation   =')
      z0 =INPUT(' Enter z_min   =')
      dz0=INPUT(' Enter dz      =')
      BoxV =Box*100.    ! Box size in km/s
      ScaleV = BoxV/AEXPN/NGRID  ! scale factor for Velocities
      ScaleC = Box/NGRID         ! scale factor for Coordinates
      write (*,*) Box,ScaleV,ScaleC

C
      CALL ReadParticles(ScaleV)
      xsh = 3.2/ScaleC
      ysh = -12.3/ScaleC
      zsh = 0.2/ScaleC
      CALL ShiftParticles(xsh,ysh,zsh)
      CALL RotateYParticles(alpha)
      Write (*,*) Xpb(1),Ypb(1),Zpb(1)
      Write (*,*) Xpb(2000),Ypb(2000),Zpb(2000)
      Write (*,*) Xpb(10000),Ypb(10000),Zpb(10000)
       CALL DENSIT    !   Find density
            write (*,*) ' Density is ok'
            Do i=1,1
              CALL SmoothDENSIT
            Enddo

            
      open(18,file='particlesZ.dat')
      WRITE (18,100) HEADER,
     +                  AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,wspecies(1),
     +                  EKIN,
     +                  NROWC,NGRID,NRECL,Om0,Oml0,hubble,
     +                  Ocurv,lspecies(1),Nspecies
      WRITE (17,*) ' Angle=',alpha,' Z=',z0,z0+dz0
      WRITE (19,*) ' Angle=',alpha,' Z=',z0,z0+dz0
      WRITE (18,*) ' Angle=',alpha,' Z=',z0,z0+dz0
      iz0 =INT(z0*Ndim1/Box+0.5)
      idz0 =INT(dz0/2.*Ndim1/Box+0.5)
      CALL DumpParticles(z0,dz0,ScaleC)
      CALL DumpDENSIT(iz0+idz0,idz0)
      STOP
      END


C---------------------------------------
C                                                       Read Particles 
      SUBROUTINE ReadParticles(ScaleV)
C---------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMfields.h'
      	XN   =FLOAT(NGRID)+1.-1.E-7
	      YN   =FLOAT(NGRID)

      W     = (FLOAT(NGRID)/FLOAT(NROW))**3
      PARTW = W  
      WPAR = wspecies(1)
      N =0
      xmin =1.e+5
      xmax =-1.e+5
      ymin =1.e+5
      ymax =-1.e+5
      zmin =1.e+5
      zmax =-1.e+5
      vxmin =1.e+15
      vxmax =-1.e+15
      vymin =1.e+15
      vymax =-1.e+15
      vzmin =1.e+15
      vzmax =-1.e+15

C				       Loop over particles
         N_particles =lspecies(Nspecies)   ! Total number of particles
         Npages = (N_particles -1)/NPAGE +1
         N_in_last=N_particles -NPAGE*(Npages-1)
      write (*,*) Nspecies,
     + lspecies(Nspecies),wspecies(Nspecies)
      Nsp_current =1
      WPAR     = wspecies(1)
      Lspec_next=lspecies(1)+1
      DO IROW=1,Npages
            In_page =NPAGE
            If(IROW.eq.Npages)In_page =N_in_last
            If(IROW/50*50.eq.IROW)
     +        write (*,*) ' Read page=',IROW,' N=',N,
     +             ' Npages=',Npages,In_page
            iL = NPAGE*(IROW-1)
C				       Loop over particles
         READ  (21,REC=IROW) RECDAT
	   DO   IN=1,In_page
            N   = N +1
            If(N.ge.Lspec_next)Then   ! next specie
               write (*,*) ' Another specie. Current=',Nsp_current
               write (*,*) '     Particle=',N,WPAR,
     &                         ' lspecies=',lspecies(Nsp_current),
     &                         ' next=',lspecies(Nsp_current+1)
               Nsp_current=Nsp_current+1
               Lspec_next =lspecies(Nsp_current)+1
               WPAR     = wspecies(Nsp_current)
            EndIf
	        X=XPAR(IN)
	        Y=YPAR(IN)
	        Z=ZPAR(IN)
           If(N.gt.Nmaxp)Then
              write (*,*) ' Error: not enough space for particles:',N
              stop
           Endif 
           Xp(N) =X                             ! coordinates in 1-NGRID
           Yp(N) =Y
           Zp(N) =Z
           Vxp(N)=ScaleV* VX(IN)    ! velocities in km/s
           Vyp(N)=ScaleV* VY(IN)
           Vzp(N)=ScaleV* VZ(IN)
           Wp(N) =WPAR                   ! particle weight

           xmin =MIN(xmin,X)
           ymin =MIN(ymin,Y)
           zmin =MIN(zmin,Z)
           xmax =MAX(xmax,X)
           ymax =MAX(ymax,Y)
           zmax =MAX(zmax,Z)
           vxmin =MIN(vxmin,Vxp(N))
           vymin =MIN(vymin,Vyp(N))
           vzmin =MIN(vzmin,Vzp(N))
           vxmax =MAX(vxmax,Vxp(N))
           vymax =MAX(vymax,Vyp(N))
           vzmax =MAX(vzmax,Vzp(N))
	        I=INT(X)
	        J=INT(Y)
	        K=INT(Z)
            If(I.le.0)write (*,*) ' X:',X,Y,Z,I,' Irow=',IROW,IN,ifile
            If(J.le.0)write (*,*) ' Y:',X,Y,Z,I,' Irow=',IROW,IN,ifile
            If(K.le.0)write (*,*) ' Z:',X,Y,Z,I,' Irow=',IROW,IN,ifile

	      ENDDO
        ENDDO 

        Np = N
      write (*,*) ' Number of particles read=', N
      write (*,*) ' X:',xmin,xmax,' V=',vxmin,vxmax
      write (*,*) ' Y:',ymin,ymax,' V=',vymin,vymax
      write (*,*) ' Z:',zmin,zmax,' V=',vzmin,vzmax
        RETURN
      END
C---------------------------------------
C                            write particles in a slice
      SUBROUTINE DumpParticles(z0,dz0,ScaleC)
C---------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMfields.h'
C                 Find points close to the center
      Xc =NGRID/2.+1.
      Yc =NGRID/2.+1.
      Zc =NGRID/2.+1.
      Rc =1.5
 10   R2 =Rc**2
C				       Loop over particles
      Nc =0
      Vxc =0.
      Vyc =0.
      Vzc =0.
	   DO   IN=1,Np
	        dX=(Xpb(IN)-Xc)            !coords
	        dY=(Ypb(IN)-Yc)
	        dZ=(Zpb(IN)-Zc)
           dd =dX**2 +dY**2 +dZ**2
             If(dd.lt.R2)Then
                Nc =Nc+1
                Vxc =Vxc +Vxpb(IN)
                Vyc =Vyc +Vypb(IN)
                Vzc =Vzc +Vzpb(IN)
             endif 
      Enddo 
      If(Nc.eq.0)Then
         Rcold =Rc
         Rc       =Rc*1.2
         write (*,*) ' Too small radius=',Rcold,' Increase it to ',Rc
         goto 10
      endif 
      Vxc =Vxc/Nc
      Vyc =Vyc/Nc
      Vzc =Vzc/Nc
      Xc   =ScaleC*(Xc-1.)
      Yc   =ScaleC*(Yc-1.)
      Zc   =ScaleC*(Zc-1.)

      write (*,*) ' There are ',Nc,' particles inside Rc=',Rc*ScaleC,
     &               ' Center(Mpc/h) =',Xc,
     &               Yc,Zc,
     &               ' V(km/s)=', Vxc, Vyc, Vzc
C				       Loop over particles
	   DO   IN=1,Np
	        X=ScaleC*(Xpb(IN)-1.)            !coords
	        Y=ScaleC*(Ypb(IN)-1.)
	        Z=ScaleC*(Zpb(IN)-1.)
             If(Z.gt.z0.and.Z.lt.z0+dz0)Then
	        dX=(X-Xc)            !coords
	        dY=(Y-Yc)
	        dZ=(Z-Zc)
              dd =dX**2 +dY**2 +dZ**2
              d   =sqrt(dd+1.d-10)
                ux =Vxpb(IN)  ! -Vxc
                uy =Vypb(IN)  ! -Vyc
                uz =Vzpb(IN)  ! -Vxc
c                vr  =(dX*ux+dY*uy+dZ*uz)/d
c                Xn =dX*(1.+vr/100./d)+Xc
c                Yn =dY*(1.+vr/100./d)+Yc
c                Zn =dZ*(1.+vr/100./d)+Zc

                write (18,100) X,Y,Z ! ,ux,uy,uz,Wp(IN) !,vr,
c     &                 Xn,Yn,Zn
 100            format(3f8.3,3f8.1,g11.3,4f8.3)
             endif 
      Enddo 
      Return
      End
C---------------------------------------
C                          Rotate around Y-axis by angle alpha
      SUBROUTINE RotateYParticles(alpha)
C---------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMfields.h'
      XN =NGRID+1
      XC =NGRID/2.+1.     ! center of rotation
      ZC =NGRID/2.+1.
      cosA =cos(alpha*3.1415926/180.)
      sinA =sin(alpha*3.1415926/180.)
C				       Loop over particles
	   DO   IN=1,Np
	        X=Xp(IN) -XC            !coords
	        Y=Yp(IN)
	        Z=Zp(IN) -ZC
           Xnew =X*cosA -Z*sinA + XC
           Znew =Z*cosA +X*sinA + ZC
           Vxpb(IN) =Vxp(IN)*cosA-Vzp(IN)*sinA
           Vypb(IN) =Vyp(IN) 
           Vzpb(IN) =Vzp(IN)*cosA+Vxp(IN)*sinA
           
           X =Xnew
              If(X.ge.XN)X=X-NGRID
              If(X.lt.1.)X=X+NGRID
           Z =Znew
              If(Z.ge.XN)Z=Z-NGRID
              If(Z.lt.1.)Z=Z+NGRID
           Xpb(IN) =X
           Ypb(IN) =Y
           Zpb(IN) =Z

      Enddo 
      Return
      End
C---------------------------------------
C                          Shift particles with periodical conditions
      SUBROUTINE ShiftParticles(xsh,ysh,zsh)
C---------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMfields.h'
      XN =NGRID+1
C				       Loop over particles
	   DO   IN=1,Np
	        X=Xp(IN)             !coords
	        Y=Yp(IN)
	        Z=Zp(IN)
           X =X+xsh
              If(X.ge.XN)X=X-NGRID
              If(X.lt.1.)X=X+NGRID
           Y =Y+ysh
              If(Y.ge.XN)Y=Y-NGRID
              If(Y.lt.1.)Y=Y+NGRID
           Z =Z+zsh
              If(Z.ge.XN)Z=Z-NGRID
              If(Z.lt.1.)Z=Z+NGRID
           Xp(IN) =X
           Yp(IN) =Y
           Zp(IN) =Z
c           Vxpb(IN) =Vxp(IN)
c           Vypb(IN) =Vyp(IN) 
c           Vzpb(IN) =Vzp(IN)

      Enddo 
      Return
      End
C---------------------------------------
C                                                       Density Field    
      SUBROUTINE SmoothDENSIT
C---------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMfields.h'
      	XN   =FLOAT(NGRID)+1.-1.E-7
	      YN   =FLOAT(NGRID)
         s    = 1.                 ! gaussian filter
         s2  = 2.*s**2
         w0 = 1.                 ! weights
         w1 = exp(-1./s2)
         w2 = exp(-2./s2)
         w3 = exp(-3./s2)
         sw = w0 +6.*w1 +12.*w2 +8.*w3
         w0 = w0/sw           ! normalize weights
         w1 = w1/sw
         w2 = w2/sw
         w3 = w3/sw
C				       sum contributions
C.... Open_MP
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE (k,km,kp,j,jm,jp,i,im,ip)
       DO k=Nl+1,Nr-1
         km = k -1
         if(km.lt.1)km=km+Ndim1
         kp = k +1
         if(kp.gt.Ndim1)kp=kp-Ndim1
       DO j=1,Ndim1
         jm = j -1
         if(jm.lt.1)jm=jm+Ndim1
         jp = j +1
         if(jp.gt.Ndim1)jp=jp-Ndim1
	      DO i=1,Ndim1
            im = i -1
            if(im.lt.1)im=im+Ndim1
            ip = i +1
            if(ip.gt.Ndim1)ip=ip-Ndim1
	        Dens2(i,j,k) = Dens1(i,j,k)*w0 +
     &            ( Dens1(i,j,km) + Dens1(i,j,kp) +
     &              Dens1(im,j,k) + Dens1(ip,j,k) +                           
     &              Dens1(i,jm,k) + Dens1(i,jp,k) )*w1 +
     &            ( Dens1(i,jm,km) + Dens1(i,jp,km) +
     &              Dens1(im,j,km) + Dens1(ip,j,km) +                           
     &              Dens1(im,jm,k) + Dens1(ip,jm,k) +
     &              Dens1(im,jp,k) + Dens1(ip,jp,k) +
     &              Dens1(i,jm,kp) + Dens1(i,jp,kp) +
     &              Dens1(im,j,kp) + Dens1(ip,j,kp) )*w2 +                           
     &            ( Dens1(im,jm,km) + Dens1(ip,jm,km) +
     &              Dens1(im,jp,km) + Dens1(ip,jp,km) +                           
     &              Dens1(im,jm,kp) + Dens1(ip,jm,kp) +
     &              Dens1(im,jp,kp) + Dens1(ip,jp,kp) )*w3
	      END DO
       END DO
      END DO
c                                restore the density
      dav =0.
      DO k=Nl,Nr
       DO j=1,Ndim1
	      DO i=1,Ndim1
	        Dens1(i,j,k) = Dens2(i,j,k)
              dav =dav + Dens1(i,j,k)
	      END DO
       END DO
      END DO
      write (*,*) ' Done smooting with filter r=',s,
     &     ' Average=',dav/ Ndim1**2/(Nr-Nl+1)     
      Return
      End
C---------------------------------------
C                                                       Density Field    
      SUBROUTINE DENSIT
C---------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMfields.h'
      	XN   =FLOAT(NGRID)+1.-1.E-7
	      YN   =FLOAT(NGRID)
C				       
      DO M3=Nl,Nr
       DO M2=1,Ndim1
	      DO M1=1,Ndim1
	        Dens1(M1,M2,M3) = 0.
	        Vx1(M1,M2,M3) = 0.
	        Vy1(M1,M2,M3) = 0.
	        Vz1(M1,M2,M3) = 0.
	      END DO
       END DO
      END DO


      xmin =1.e+5
      xmax =-1.e+5
      ymin =1.e+5
      ymax =-1.e+5
      zmin =1.e+5
      zmax =-1.e+5
      Scale1 = Float(Ndim1)/float(NGRID)
      AvDens1 =NGRID**3/float(Ndim1)**3
C				       Loop over particles
	   DO   IN=1,Np
	        X=Xpb(IN)             !coords
	        Y=Ypb(IN)
	        Z=Zpb(IN)
              Ux =Vxpb(IN)        ! velocities
              Uy =Vypb(IN)
              Uz =Vzpb(IN)
              W  =Wp(IN)           ! weight
	        X1=(X-1.)*Scale1 +1.     ! scale to mesh 1
	        Y1=(Y-1.)*Scale1 +1.
	        Z1=(Z-1.)*Scale1 +1.
           xmin =MIN(xmin,X1)
           ymin =MIN(ymin,Y1)
           zmin =MIN(zmin,Z1)
           xmax =MAX(xmax,X1)
           ymax =MAX(ymax,Y1)
           zmax =MAX(zmax,Z1)
	        I1=INT(X1)                      ! fields for the first mesh
	        J1=INT(Y1)
	        K1=INT(Z1)
c           write(*,110) IN,X,Y,Z,X1,Y1,Z1,I1,J1,K1
c 110       format(i8,6f9.4,3i4)
            If(I1.le.0)write (*,*) ' Error X1:',X1,Y1,Z1,In,xmin,xmax
            If(J1.le.0)write (*,*) ' Error Y1:',X1,Y1,Z1,In,ymin,ymax
            If(K1.le.0)write (*,*) ' Error  Z1:',X1,Y1,Z1,In,zmin,zmax
	        D1=X1-FLOAT(I1)
	        D2=Y1-FLOAT(J1)
	        D3=Z1-FLOAT(K1)
	        T1=1.-D1
	        T2=1.-D2
	        T3=1.-D3
	        T2W =T2*W
	        D2W =D2*W
	        Ia=I1+1
	           IF(Ia.GT.Ndim1)Ia=1
	        Ja=J1+1
	           IF(Ja.GT.Ndim1)Ja=1
	        Ka=K1+1
	           IF(Ka.GT.Ndim1)Ka=1

              IF(K1.ge.Nl.and.K1.le.Nr)Then
		      Dens1(I1 ,J1 ,K1 ) =Dens1(I1 ,J1 ,K1 ) +T3*T1*T2W
		      Dens1(Ia,J1 ,K1 ) =Dens1(Ia,J1 ,K1 ) +T3*D1*T2W
		      Dens1(I1 ,Ja,K1 ) =Dens1(I1 ,Ja,K1 ) +T3*T1*D2W
		      Dens1(Ia,Ja,K1 ) =Dens1(Ia,Ja,K1 ) +T3*D1*D2W
		      Vx1(I1 ,J1 ,K1 ) =Vx1(I1 ,J1 ,K1 ) +T3*T1*T2W*Ux
		      Vx1(Ia,J1 ,K1 ) =Vx1(Ia,J1 ,K1 ) +T3*D1*T2W*Ux
		      Vx1(I1 ,Ja,K1 ) =Vx1(I1 ,Ja,K1 ) +T3*T1*D2W*Ux
		      Vx1(Ia,Ja,K1 ) =Vx1(Ia,Ja,K1 ) +T3*D1*D2W*Ux
		      Vy1(I1 ,J1 ,K1 ) =Vy1(I1 ,J1 ,K1 ) +T3*T1*T2W*Uy
		      Vy1(Ia,J1 ,K1 ) =Vy1(Ia,J1 ,K1 ) +T3*D1*T2W*Uy
		      Vy1(I1 ,Ja,K1 ) =Vy1(I1 ,Ja,K1 ) +T3*T1*D2W*Uy
		      Vy1(Ia,Ja,K1 ) =Vy1(Ia,Ja,K1 ) +T3*D1*D2W*Uy
		      Vz1(I1 ,J1 ,K1 ) =Vz1(I1 ,J1 ,K1 ) +T3*T1*T2W*Uz
		      Vz1(Ia,J1 ,K1 ) =Vz1(Ia,J1 ,K1 ) +T3*D1*T2W*Uz
		      Vz1(I1 ,Ja,K1 ) =Vz1(I1 ,Ja,K1 ) +T3*T1*D2W*Uz
		      Vz1(Ia,Ja,K1 ) =Vz1(Ia,Ja,K1 ) +T3*D1*D2W*Uz
            EndIf
    
              IF(Ka.ge.Nl.and.Ka.le.Nr)Then
		      Dens1(I1 ,J1 ,Ka) =Dens1(I1 ,J1 ,Ka) +D3*T1*T2W
		      Dens1(Ia,J1 ,Ka) =Dens1(Ia,J1 ,Ka) +D3*D1*T2W
		      Dens1(I1 ,Ja,Ka) =Dens1(I1 ,Ja,Ka) +D3*T1*D2W
		      Dens1(Ia,Ja,Ka) =Dens1(Ia,Ja,Ka) +D3*D1*D2W
    
		      Vx1(I1 ,J1 ,Ka) =Vx1(I1 ,J1 ,Ka) +D3*T1*T2W*Ux
		      Vx1(Ia,J1 ,Ka) =Vx1(Ia,J1 ,Ka) +D3*D1*T2W*Ux
		      Vx1(I1 ,Ja,Ka) =Vx1(I1 ,Ja,Ka) +D3*T1*D2W*Ux
		      Vx1(Ia,Ja,Ka) =Vx1(Ia,Ja,Ka) +D3*D1*D2W*Ux
    
		      Vy1(I1 ,J1 ,Ka) =Vy1(I1 ,J1 ,Ka) +D3*T1*T2W*Uy
		      Vy1(Ia,J1 ,Ka) =Vy1(Ia,J1 ,Ka) +D3*D1*T2W*Uy
		      Vy1(I1 ,Ja,Ka) =Vy1(I1 ,Ja,Ka) +D3*T1*D2W*Uy
		      Vy1(Ia,Ja,Ka) =Vy1(Ia,Ja,Ka) +D3*D1*D2W*Uy

    
		      Vz1(I1 ,J1 ,Ka) =Vz1(I1 ,J1 ,Ka) +D3*T1*T2W*Uz
		      Vz1(Ia,J1 ,Ka) =Vz1(Ia,J1 ,Ka) +D3*D1*T2W*Uz
		      Vz1(I1 ,Ja,Ka) =Vz1(I1 ,Ja,Ka) +D3*T1*D2W*Uz
		      Vz1(Ia,Ja,Ka) =Vz1(Ia,Ja,Ka) +D3*D1*D2W*Uz
            EndIf
            ENDDO 
      write (*,*) '    Min/Max coords for the first mesh'
      write (*,*) ' X:',xmin,xmax
      write (*,*) ' Y:',ymin,ymax
      write (*,*) ' Z:',zmin,zmax
C				                              Normalize the fields
      AvD = 0.
      Dmin =1.e+10
      Dmax =-1.e+10
      RD  =0.
      sx    =0.
      sy    =0.
      sz    =0.
      DO M3=Nl,Nr
       DO M2=1,Ndim1
	      DO M1=1,Ndim1
	        Vx1(M1,M2,M3) = Vx1(M1,M2,M3)/max(1.d-10,Dens1(M1,M2,M3))
	        Vy1(M1,M2,M3) = Vy1(M1,M2,M3)/max(1.d-10,Dens1(M1,M2,M3))
	        Vz1(M1,M2,M3) = Vz1(M1,M2,M3)/max(1.d-10,Dens1(M1,M2,M3))
              If(Dens1(M1,M2,M3).lt.1.d-10)Then
                 Vx1(M1,M2,M3) =0.
                 Vy1(M1,M2,M3) = 0.
      	        Vz1(M1,M2,M3) =0.
              endif 
	        Dens1(M1,M2,M3) = Dens1(M1,M2,M3)/Avdens1
              AvD =AvD + Dens1(M1,M2,M3)
              Dmin =MIN(Dmin,Dens1(M1,M2,M3))
              Dmax =MAX(Dmax,Dens1(M1,M2,M3))
              sx     =sx + Vx1(M1,M2,M3)**2
              sy     =sy + Vy1(M1,M2,M3)**2
              sz     =sz + Vz1(M1,M2,M3)**2
	      END DO
       END DO
      END DO
      AvD =AvD/Ndim1**2/(Nr-Nl+1)
      sx =sqrt(sx/Ndim1**2/(Nr-Nl+1))
      sy =sqrt(sy/Ndim1**2/(Nr-Nl+1))
      sz =sqrt(sz/Ndim1**2/(Nr-Nl+1))
      write (*,*) '<----------- Fields for Mesh=',Ndim1,' ------>' 
      write (*,*) '    Density: Max=',Dmax,' Min=',Dmin,' Average=',AvD
      write (*,*) '    Velocity: rms=',sx,sy,sz
      write (17,*) '<----------- Fields for Mesh=',Ndim1,' ------>' 
      write (17,*) '    Density: Max=',Dmax,' Min=',Dmin,' Average=',AvD
      write (17,*) '    Velocity: rms=',sx,sy,sz

c      enddo
        RETURN
      END

C---------------------------------------
C                                                     Dump  Density Field    
      SUBROUTINE DumpDENSIT(k,kd)
C---------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMfields.h'
      	XN   =FLOAT(NGRID)+1.-1.E-7
	      YN   =FLOAT(NGRID)

         write (*,*) ' dump Density limits:',max(k-kd,Nl),Min(k+kd,Nr)
      do j=1,Ndim1
      do i=1,Ndim1
         Dd =0.
         Wx =0.
         Wy =0.
         Wz =0.
         Do kk=k-kd,k+kd
            If(kk.ge.Nl.and.kk.le.Nr)Then
               Dd =Dd +Dens1(i,j,kk)
               Wx =Wx +Vx1(i,j,kk)*Dens1(i,j,kk)
               Wy =Wy +Vy1(i,j,kk)*Dens1(i,j,kk)
               Wz =Wz +Vz1(i,j,kk)*Dens1(i,j,kk)
            endif
         EndDo
            Dd =MAX(Dd,1.d-10)
        write (17,100) i,j,k,
     &           Dd,Wx/Dd,Wy/Dd,Wz/Dd
        if(i/2*2.eq.i.and.j/2*2.eq.j)
     &            write (19,100) i,j,k,
     &           Dd,Wx/Dd,Wy/Dd,Wz/Dd
 100     format(3i4,4g11.4)
      enddo 
      enddo 
      write (*,*) ' Done density dump for K=',k
        RETURN
      END