C _______________________ START 3-D SIMULATIONS
C
C                         Klypin, Jenuary 1997
C
      INCLUDE 'PMparameters.h'
C		   NROW = number of particles per row
C		   NGRID= number of cells on a side (normal grid)
C		   Npage =  number of points per page of the grid
C
      COMMON / GRID/	GRX(NROW,NROW,NROW),GRY(NROW,NROW,NROW),
     +			GRZ(NROW,NROW,NROW)
      COMMON / GRIV/	GVX(NROW,NROW,NROW),GVY(NROW,NROW,NROW),
     +			GVZ(NROW,NROW,NROW)
      COMMON / ROW2/	XPAR2(NPAGE),YPAR2(NPAGE),ZPAR2(NPAGE),
     +			VX2(NPAGE),VY2(NPAGE),VZ2(NPAGE)
      COMMON / KINENRG/ SKINE,SX,SY,SZ,SX2,SY2,SZ2
      DATA PI		/3.1415926535/
      CALL InitValues(NBYTE,SCALEL)  ! Set parameters of the run
         write (*,*) '  InitValues is done. N bytes=',NBYTE
         write (*,*) '  Nspecies =',Nspecies
         write (16,*) ' Seed=',Nseed
      OPEN(UNIT=9,FILE='PMcrd.DAT',form='unformatted')
      write (9)                ! this   clears old header file
      CLOSE (9)
      CALL RDTAPE                  ! this opens files on disk         
      OPEN(UNIT=16,FILE='RESNEW.DAT') ! open protocol
C
        AEXP0 = AEXPN            ! expansion parameter for coords
        AEXPV = AEXPN - ASTEP/2. !                     for veloc
C----                        Cold Particles --------------------------      
        NRAND =Nseed             ! store the initial seed
      CALL SPECTR(ALPHA,NRAND) ! get a realization of spectrum
        IFOUR = INT(ALOG(FLOAT(NROW))/ALOG(2.)+0.5)
        write (*,*) ' Ifour =',IFOUR,' Nrow=',NROW
      CALL VECTOR(IFOUR)       ! get displacement vector by FFT
   	write (*,*) ' SPECTR is done. IFOUR = ',IFOUR,' Alpha=',ALPHA
        VCONS = -ALPHA/(2.*PI/NGRID)*(AEXPV/AEXP0)*SQRT(AEXPV)
        XCONS =	ALPHA/(2.*PI/NGRID)*(AEXPN/AEXP0)
        QFACT = FLOAT(NGRID)/FLOAT(NROW)
        write (*,*) ' Aexp=',AEXPN,AEXP0,AEXPV,ASTEP
                      SKINE = 0.   ! set counters to zero
                      SX    = 0.
                      SY    = 0.
                      SZ    = 0.
                      SX2   = 0.
                      SY2   = 0.
                      SZ2   = 0.
      DO  KROW = 1, NROW  ! find x,v and write to file, page by page
	     CALL SETXV(XCONS,VCONS,KROW)
	     CALL WRIROW(KROW,1)
      ENDDO
      EKIN = 0.5*SKINE/AEXPV**2*PARTW !kinet.energy at A(i-1/2)
      DP   = NROW**3
      SX   = SX/DP                    ! mean displacement
      SY   = SY/DP
      SZ   = SZ/DP
      SX2  = SQRT(SX2/DP)             ! rms displacement
      SY2  = SQRT(SY2/DP)
      SZ2  = SQRT(SZ2/DP)
      WRITE (*,'('' Displacement of a particle: MEAN'',3F9.4)')
     .	     SX,SY,SZ
      WRITE (*,'(''                             RMS '',3F9.4)')
     .	     SX2,SY2,SZ2
      WRITE (16,'('' Displacement of a particle: MEAN'',3F9.4)')
     .	     SX2,SY2,SZ2
      Veloc=sqrt( SKINE/AEXPV**2/NROW**3)*SCALEL/NGRID*50.
      WRITE (*,'('' Ekin='',E12.4,'' Veloc_cold(km/s)='',F8.2)')
     .       EKIN,Veloc
      WRITE (16,'('' Ekin='',E12.4,'' Veloc_cold(km/s)='',F8.2)')
     .       EKIN,Veloc

C----                        Neutrino ------------------------------      
        NRAND =Nseed       !  Restore initial seed
      CALL SPECNU(NRAND) !  Renormalize the spectrum for neutrino
      CALL VECTOR(IFOUR) !   get displacement vector by FFT
	     write (*,*) ' Neutrino SPECTRUM is done.'
C			     find x,v and write to file, page by page
         Svel =SKINE
         SKINE=0.
         DO  KROW = 1, NROW
           Do js=1,Nspecies,2
              CALL SETXVN(XCONS,VCONS,KROW,NRAND)
              CALL WRIROW(KROW,js+1)
              If(Nspecies.ne.js)Then ! write another pair of
                 Do I=1,NPAGE        ! particles
                    XPAR(I)=XPAR2(I)
                    YPAR(I)=YPAR2(I)
                    ZPAR(I)=ZPAR2(I)
                    VX(I)  =VX2(I)
                    VY(I)  =VY2(I)
                    VZ(I)  =VZ2(I)
                 EndDo
                 CALL WRIROW(KROW,js+2)
              EndIf
           EndDo
         ENDDO
      W    = (FLOAT(NGRID)/FLOAT(NROW))**3         
      EKIN = 0.5*(Svel*PARTW+SKINE*(W-PARTW)/Nspecies)/AEXPV**2
      Veloc= sqrt( SKINE/AEXPV**2/NROW**3/Float(Nspecies))
     .           *SCALEL/NGRID*50. 
      WRITE (16,'('' Ekin='',E12.4,'' Veloc_hot(km/s)='',F8.2)')
     .       EKIN,Veloc
      WRITE (*,'('' Ekin='',E12.4,'' Veloc_hot(km/s)='',F8.2)')
     .       EKIN,Veloc
      CALL WRTAPE    !   write header and control data
      STOP
      END
C************************** MDM:    0.1-0.3-0.6    7eV
c      FUNCTION TRUNF(wk)
c      INCLUDE 'PMparameters.h'
c      PARAMETER (NSPEC =NROW/2)
c      COMMON / TRUNCOM/ QSCALE, ns
c      DATA A1,A2,A3,A4,A5 
c     .         /3.93, - 60.2, 409., -17.1, 12.5/
c       IF (WK.GE.FLOAT(NSPEC)) THEN
c	  TRUNF =0.
c	  RETURN
c       ENDIF
c       Q = QSCALE *WK
c       sk=sqrt(Q)
c       TRUNF = sqrt(LOG(1.+A5*g*Q) / (A5*g)) /
c     .         (1.+sk*(A1*(1.-g/3.2)
c     .            +sk*(A2*(1.-g/3.5)
c     .            +sk*(A3*(1.-g/5.)
c     .            +sk*(A4*(1.-g**2/15.) )))) )
c      RETURN
c      END
C************************** MDM: Om_nu=0.2, z=30. 
      FUNCTION TRUNF(wk)      
      INCLUDE 'PMparameters.h'
      PARAMETER (NSPEC =NROW/2)
      COMMON / TRUNCOM/ QSCALE, ns
      DIMENSION          PP(6)
      REAL                ns
      DATA   PP/1.,-1.464,19.81,154.5,143.6,0.961/
       IF (WK.GE.FLOAT(NSPEC)) THEN
	      TRUNF =0. ! truncate above nyquist frequency
	      RETURN
       ENDIF
       Q = QSCALE*WK
       sk= sqrt(Q)
       TRUNF= sqrt(wk**ns/
     .        (1.+sk*(PP(2) +sk*(PP(3) +sk*(PP(4)
     .           +sk*PP(5) ))))**(2.*PP(6)) )
      RETURN
      END
C**************************:  Approx for Growth rate for MDM cold
      FUNCTION GrPcold(wk)
c      PARAMETER (C1=0.805,C2=0.195,C3=0.3,C4=0.70803)
      PARAMETER (C1=0.8736,C2=1.-C1,C3=0.0891,C4=0.7347)
      COMMON / TRUNCOM/ QSCALE, ns
      REAL                ns
         Q = QSCALE*WK
         GrPcold =C1+C2/(1. +exp(LOG(Q/C3)/C4))
      Return
      End
C**************************:  Approx for Growth rate for MDM hot
      FUNCTION GrPhot(wk)
c      PARAMETER (C1=1.80,C2=-0.8,C3=1.4,C4=1.065)
      PARAMETER (C1=1.75,C2=-0.75,C3=0.2165,C4=0.88)
      COMMON / TRUNCOM/ QSCALE, ns
      REAL                ns
         Q = QSCALE*WK
         GrPhot =C1+C2/(1. +exp(LOG(Q/C3)/C4))
      Return
      End
C************************** Neutrino:   Cold * correction factor
      FUNCTION TRUNU(WK)
      INCLUDE 'PMparameters.h'
      COMMON / TRUNCOM/ QSCALE, ns
      DIMENSION PP(6)
      REAL                ns
      DATA PP/1.442,-0.651,4.03,1.553,-35.50,141.7/
         Q = QSCALE*WK
         sk= (Q)**0.33333
         TRUNU =TRUNF(WK)  *
     .                    sqrt( exp(-PP(1)*Q)
     .                   /(1.+sk*(PP(2) +sk*(PP(3)
     .                       +sk*(PP(4) +sk*(PP(5) 
     .                       +sk*PP(6)    ))))) )
      RETURN
      END
C*************************************** Velocity distribution
      SUBROUTINE VELCT(Wx,Wy,Wz,NRAND)
      COMMON / FERMI  / Vscale
C                                       get random unit vector
 1    x =(RANDd(NRAND)-0.5)*2.
      y =(RANDd(NRAND)-0.5)*2.
      z =(RANDd(NRAND)-0.5)*2.
      r2 =x**2 +y**2 +z**2
            If(r2 .gt. 1.)goto1
C                                      Fermi distribution            
      p  =RANDd(NRAND)
      v  =Vscale *Vfermi(p)/sqrt(r2)
      Wx =x *v
      Wy =y *v
      Wz =z *v

      Return
      End
C*********************************************************************
C			  INITIALIZE CONSTANTS:
C			 Scalel -length of square in MPC
C			 AEXPN	-expansion factor (current)
C			 ASTEP	-increment of expansion factor
C			 PARTW	-weight of a particle
C			 AMPLT	-amplitude of perturbations
C			 NBYTE	-record length in bytes
C
      SUBROUTINE InitValues(NBYTE,SCALEL)
      INCLUDE 'PMparameters.h'
      COMMON / TRUNCOM/ QSCALE, ns
      COMMON / FERMI  / Vscale
      REAL                INPUT, ns
      DATA PI		 /3.1415926535/
C                          '123456789012345678901234567890123456789012345'   
       HEADER='N=256x512 L=50/hsig=0.174totCHDM2nu z0=30Ns=3'
       Om0   = 1.         ! \Omega_0 at z=0
       Oml0  = 0.         ! \Omega_lambda at z=0
       Ocurv = 0.
       hubble= 0.5          
       write (*,*) '------ Enter Header for the run up to 45 characters'
       write (*,*) '       Example:'
       write (*,'(A)') HEADER
       read  (*,'(A)') HEADER
       write (*,'(A)') HEADER
       AEXPN =1./(1.+30.)   !  z_init = 30.
       write (*,*) ' Initial expansion parameter     =',AEXPN
       ASTEP =INPUT(' Step of the expansion parameter  =')
       AMPLT =INPUT(' Amplitude of density fluctuations=')
       ns    =INPUT(' Slope of the power spectrum      =')
       SCALEL=INPUT(' Box size (Mpc/h)                 =')
       Nseed =INPUT(' Random seed (Integer  1-2^31)    =') 
       Nspecies =INPUT(' Number of neutrino/cold_part(1-6)=') 
       SCALEL=SCALEL/hubble   ! scale it to real megaparsecs  
       write (*,*) ' Aexpn=',AEXPN,' dA=',ASTEP
C                          Velocity scale for neutrino thermal vel.
C                          Vrms = 3.59*Vscale
C                          No dependance on Aexpn (P, not V)       
c       Vscale= 7.2/(SCALEL*50./NGRID)
c       Vscale= 8.82/(SCALEL*50./NGRID)  ! m_nu=5.72eV, Omega_nu=0.25
       Vscale= 22.0/(SCALEL*50./NGRID) ! m_nu=2.29eV, Omega_nu=0.1x2
       NBYTE = NPAGE*6*4
       W     = (FLOAT(NGRID)/FLOAT(NROW))**3
       PARTW = W *0.80       ! this is for 20% neutrino
c                                 W1    = (W -Partw)/Nspecies
       ISTEP = 0
       TINTG = 0.
       AU0   = 0.
       AEU0  = 0.
       EKIN  = 0.
       EKIN1 = 0.
       EKIN2 = 0.
       NROWC = NROW
       NGRIDC= NGRID
	    QSCALE = 2.*PI/SCALEL
      RETURN
      END
C*********************************************************************
C		    Set spectrum of perturbations:
      SUBROUTINE SPECTR(ALPHA,NRAND)
      INCLUDE 'PMparameters.h'
      PARAMETER (NSPEC = NROW/2 )
      COMMON / GRID/	GRX(NROW,NROW,NROW),GRY(NROW,NROW,NROW),
     +			GRZ(NROW,NROW,NROW)
      COMMON / GRIV/	GVX(NROW,NROW,NROW),GVY(NROW,NROW,NROW),
     +			GVZ(NROW,NROW,NROW)
C						    Set spectrum
      SUMM = 0.
      DO  MK3 = 1,NROW
      DO  MK2 = 1,NROW
      DO  MK1 = 1,NROW
	       GRX(MK1,MK2,MK3) =0.
	       GRY(MK1,MK2,MK3) =0.
	       GRZ(MK1,MK2,MK3) =0.
      ENDDO
      ENDDO
      ENDDO

      DO  k = 1,NROW
         ksign = -1
         kz    = k +NSPEC
         k3    = k - 1
         If(k3.GT.NSPEC)k3=k3-NSPEC
           IF(kz.gt.NROW)THEN
              kz    =kz -NROW
              ksign =1
           ENDIF
           WK3= K3**2
           DO  j = 1,NROW
              jsign = -1
              jz    = j +NSPEC
              j3    = j - 1
              If(j3.GT.NSPEC)j3=j3-NSPEC
                IF(jz.gt.NROW)THEN
                   jz    =jz -NROW
                   jsign =1
                ENDIF
                WJ3= j3**2
                DO  i = 1,NROW
                   isign = -1
                   iz    = i +NSPEC
                   i3    = i - 1
                   If(i3.GT.NSPEC)i3=i3-NSPEC
                     IF(iz.gt.NROW)THEN
                        iz    =iz -NROW
                        isign =1
                     ENDIF
                     WI3= i3**2

	       IF(i3+j3+k3.EQ.0)THEN
		  GRX(1,1,1) = 0.
		  GRY(1,1,1) = 0.
		  GRZ(1,1,1) = 0.
	       ELSE
		  WD = WI3 + WJ3 + WK3
		  WK = SQRT(WD)

		  TS = TRUNF(WK) * GAUSS(NRAND)
		  TRX = TS / WD
                  TV = GrPcold(WK)
		  GRX(iz, j, k) = i3 * TRX *isign
		  GRY( i,jz, k) = j3 * TRX *jsign
		  GRZ( i, j,kz) = k3 * TRX *ksign
		  GVX(iz, j, k) = TV * GRX(iz, j, k)
		  GVY( i,jz, k) = TV * GRY( i,jz, k)
		  GVZ( i, j,kz) = TV * GRZ( i, j,kz)
		  SUMM = SUMM + TS**2
	       ENDIF
	    ENDDO
	   ENDDO
      ENDDO
      IF(SUMM.LE.0.)WRITE (*,*) ' Error!!! Summ over spectrum = 0'
         write (*,*) ' ==== SPECTR ===  SUMM=',SUMM
      ALPHA = AMPLT / SQRT(SUMM) *sqrt(8.)
      RETURN
      END
C*********************************************************************
C		    Set spectrum of perturbations:
      SUBROUTINE SPECNU(NRAND)
      INCLUDE 'PMparameters.h'
      PARAMETER (NSPEC = NROW/2 )
      COMMON / GRID/	GRX(NROW,NROW,NROW),GRY(NROW,NROW,NROW),
     +			GRZ(NROW,NROW,NROW)
      COMMON / GRIV/	GVX(NROW,NROW,NROW),GVY(NROW,NROW,NROW),
     +			GVZ(NROW,NROW,NROW)
C						    Set spectrum
      SUMM = 0.
      DO  MK3 = 1,NROW
      DO  MK2 = 1,NROW
      DO  MK1 = 1,NROW
	      GRX(MK1,MK2,MK3) =0.
	      GRY(MK1,MK2,MK3) =0.
	      GRZ(MK1,MK2,MK3) =0.
      ENDDO
      ENDDO
      ENDDO

      DO  k = 1,NROW
         ksign = -1
         kz    = k +NSPEC
         k3    = k - 1
         If(k3.GT.NSPEC)k3=k3-NSPEC
           IF(kz.gt.NROW)THEN
              kz    =kz -NROW
              ksign =1
           ENDIF
           WK3= K3**2
           DO  j = 1,NROW
              jsign = -1
              jz    = j +NSPEC
              j3    = j - 1
              If(j3.GT.NSPEC)j3=j3-NSPEC
                IF(jz.gt.NROW)THEN
                   jz    =jz -NROW
                   jsign =1
                ENDIF
                WJ3= j3**2
                DO  i = 1,NROW
                   isign = -1
                   iz    = i +NSPEC
                   i3    = i - 1
                   If(i3.GT.NSPEC)i3=i3-NSPEC
                     IF(iz.gt.NROW)THEN
                        iz    =iz -NROW
                        isign =1
                     ENDIF
                     WI3= i3**2

	       IF(i3+j3+k3.EQ.0)THEN
		  GRX(1,1,1) = 0.
		  GRY(1,1,1) = 0.
		  GRZ(1,1,1) = 0.
	       ELSE
		  WD = WI3 + WJ3 + WK3
		  WK = SQRT(WD)

		  TS = TRUNU(WK) * GAUSS(NRAND)
		  TRX = TS / WD
                  TV = GrPhot(WK)
		  GRX(iz, j, k) = i3 * TRX *isign
		  GRY( i,jz, k) = j3 * TRX *jsign
		  GRZ( i, j,kz) = k3 * TRX *ksign
		  GVX(iz, j, k) = TV * GRX(iz, j, k)
		  GVY( i,jz, k) = TV * GRY( i,jz, k)
		  GVZ( i, j,kz) = TV * GRZ( i, j,kz)
		  SUMM = SUMM + TS**2
	       ENDIF
	    ENDDO
	   ENDDO
      ENDDO
      IF(SUMM.LE.0.)WRITE (*,*) ' Error!!! Summ over spectrum = 0'
         write (*,*) ' ==== SPECTR ===  SUMM=',SUMM
      RETURN
      END
C*********************************************************************
C		   Define coordinates and velosities for
C		   all particles of current row (given q2)
C			  x = q - amplt*b(t)/b(t0)*S(q1,q2)
C			  momentum = a**2*deriv x
C
C		   Twidled variables:
C		   x= q - {1/(2 pi Ngrid)}{Alpha*a(t)/a(0)}*S
C		   P=	- {1/(2 pi Ngrid)}{Alpha*a(t)**1.5/a(0)}*S
C
C		   q(i) = Ngrid/Nrow*(i-1) + 1
C		   x = q - xcons*S
C		   P =	   vcons*S
      SUBROUTINE SETXV(XCONS,VCONS,K)
      INCLUDE 'PMparameters.h'
      COMMON / KINENRG/ SKINE,SX,SY,SZ,SX2,SY2,SZ2
      COMMON / GRID/	GRX(NROW,NROW,NROW),GRY(NROW,NROW,NROW),
     +			GRZ(NROW,NROW,NROW)
      COMMON / GRIV/	GVX(NROW,NROW,NROW),GVY(NROW,NROW,NROW),
     +			GVZ(NROW,NROW,NROW)

      QFACT = FLOAT(NGRID)/FLOAT(NROW)
      XMAX  = FLOAT(NGRID) + 1.
      XSHF  = FLOAT(NGRID)
      Q3    = QFACT*(K-1.) +1.
      DO  J = 1,NROW
	     Q2 = QFACT*(J-1.) +1.
	     DO  I = 1,NROW
	       Q1       = QFACT*(I-1.) +1.
          IJ       =      I +(J-1)*NROW
	       DX	    =      XCONS*GRX(I,J,K)
	       DY	    =      XCONS*GRY(I,J,K)
	       DZ	    =      XCONS*GRZ(I,J,K)
	       XPAR(IJ) = Q1 - DX       + 1.E-4 +0.5
	       YPAR(IJ) = Q2 - DY       + 1.E-4 +0.5
	       ZPAR(IJ) = Q3 - DZ       + 1.E-4 +0.5
	       VX(IJ)   =      VCONS*GVX(I,J,K)
	       VY(IJ)   =      VCONS*GVY(I,J,K)
	       VZ(IJ)   =      VCONS*GVZ(I,J,K)
	       IF(XPAR(IJ).GT.XMAX)	XPAR(IJ)=XPAR(IJ)-XSHF
	       IF(XPAR(IJ).LE.1.)	XPAR(IJ)=XPAR(IJ)+XSHF
	       IF(YPAR(IJ).GT.XMAX)	YPAR(IJ)=YPAR(IJ)-XSHF
	       IF(YPAR(IJ).LE.1.)	YPAR(IJ)=YPAR(IJ)+XSHF
	       IF(ZPAR(IJ).GT.XMAX)	ZPAR(IJ)=ZPAR(IJ)-XSHF
	       IF(ZPAR(IJ).LE.1.)	ZPAR(IJ)=ZPAR(IJ)+XSHF
	       SX	    = SX    + DX
	       SY	    = SY    + DY
	       SZ	    = SZ    + DZ
	       SX2      = SX2   + DX**2
	       SY2      = SY2   + DY**2
	       SZ2      = SZ2   + DZ**2
	       SKINE    = SKINE + VX(IJ)**2 +VY(IJ)**2 +VZ(IJ)**2
	     ENDDO
      ENDDO
      RETURN
      END
C*********************************************************************
C		   Define coordinates and velosities for
C		   neutrino in the current row 
C
      SUBROUTINE SETXVN(XCONS,VCONS,K,NRAND)
      INCLUDE 'PMparameters.h'
      COMMON / KINENRG/ SKINE,SX,SY,SZ,SX2,SY2,SZ2
      COMMON / GRID/	GRX(NROW,NROW,NROW),GRY(NROW,NROW,NROW),
     +			GRZ(NROW,NROW,NROW)
      COMMON / GRIV/	GVX(NROW,NROW,NROW),GVY(NROW,NROW,NROW),
     +			GVZ(NROW,NROW,NROW)
      COMMON / ROW2/	XPAR2(NPAGE),YPAR2(NPAGE),ZPAR2(NPAGE),
     +			VX2(NPAGE),VY2(NPAGE),VZ2(NPAGE)

      QFACT = FLOAT(NGRID)/FLOAT(NROW)
      XMAX  = FLOAT(NGRID) + 1.
      XSHF  = FLOAT(NGRID)
      Q3    = QFACT*(K-1.) +1.
            
      DO  J = 1,NROW
	     Q2 = QFACT*(J-1.) +1.
	     DO  I = 1,NROW
	       Q1 = QFACT*(I-1.) +1.
          IJ       =      I +(J-1)*NROW
	       DX	    =      XCONS*GRX(I,J,K)
	       DY	    =      XCONS*GRY(I,J,K)
	       DZ	    =      XCONS*GRZ(I,J,K)
	       XPAR(IJ) = Q1 - DX       + 1.E-4 +0.5
	       YPAR(IJ) = Q2 - DY       + 1.E-4 +0.5
	       ZPAR(IJ) = Q3 - DZ       + 1.E-4 +0.5
          CALL VELCT(Wx,Wy,Wz,NRAND)
	       VX(IJ)   =      VCONS*GVX(I,J,K) +Wx
	       VY(IJ)   =      VCONS*GVY(I,J,K) +Wy
	       VZ(IJ)   =      VCONS*GVZ(I,J,K) +Wz
	       VX2(IJ)  =      VCONS*GVX(I,J,K) -Wx
	       VY2(IJ)  =      VCONS*GVY(I,J,K) -Wy
	       VZ2(IJ)  =      VCONS*GVZ(I,J,K) -Wz
	          IF(XPAR(IJ).GT.XMAX)	XPAR(IJ)=XPAR(IJ)-XSHF
	          IF(XPAR(IJ).LE.1.)	XPAR(IJ)=XPAR(IJ)+XSHF
	          IF(YPAR(IJ).GT.XMAX)	YPAR(IJ)=YPAR(IJ)-XSHF
	          IF(YPAR(IJ).LE.1.)	YPAR(IJ)=YPAR(IJ)+XSHF
	          IF(ZPAR(IJ).GT.XMAX)	ZPAR(IJ)=ZPAR(IJ)-XSHF
	          IF(ZPAR(IJ).LE.1.)	ZPAR(IJ)=ZPAR(IJ)+XSHF
          XPAR2(IJ)=      XPAR(IJ)
          YPAR2(IJ)=      YPAR(IJ)
          ZPAR2(IJ)=      ZPAR(IJ)
	       SKINE    = SKINE + VX(IJ)**2 +VY(IJ)**2 +VZ(IJ)**2
     .                     + VX2(IJ)**2+VY2(IJ)**2+VZ2(IJ)**2
   	  ENDDO
      ENDDO
      RETURN
      END
C*********************************************************************
C				   FFT of the spectrum
      SUBROUTINE VECTOR(IFOUR)
      INCLUDE 'PMparameters.h'
      COMMON / GRID/	GRX(NROW,NROW,NROW),GRY(NROW,NROW,NROW),
     +			GRZ(NROW,NROW,NROW)
      COMMON / GRIV/	GVX(NROW,NROW,NROW),GVY(NROW,NROW,NROW),
     +			GVZ(NROW,NROW,NROW)
      CALL SETF67(4,IFOUR)
C				       x-component
      DO K=1,NROW
	   DO J=1,NROW
	    DO I=1,NROW
	       Zf(I) =GRX(I,J,K)
	    ENDDO
	    CALL FOUR67(4,IFOUR)
	    DO I=1,NROW
	       GRX(I,J,K) = Yf(I)
	    ENDDO
	    DO I=1,NROW
	       Zf(I) =GVX(I,J,K)
	    ENDDO
	    CALL FOUR67(4,IFOUR)
	    DO I=1,NROW
	       GVX(I,J,K) = Yf(I)
	    ENDDO
   	ENDDO
      ENDDO
      DO K=1,NROW
	   DO I=1,NROW
	    DO J=1,NROW
	       Zf(J) =GRX(I,J,K)
	    ENDDO
	    CALL FOUR67(4,IFOUR)
	    DO J=1,NROW
	       GRX(I,J,K) = Yf(J)
	    ENDDO
	    DO J=1,NROW
	       Zf(J) =GVX(I,J,K)
	    ENDDO
	    CALL FOUR67(4,IFOUR)
	    DO J=1,NROW
	       GVX(I,J,K) = Yf(J)
	    ENDDO
	   ENDDO
      ENDDO
      DO J=1,NROW
	   DO I=1,NROW
	    DO K=1,NROW
	       Zf(K) =GRX(I,J,K)
	    ENDDO
	    CALL FOUR67(4,IFOUR)
	    DO K=1,NROW
	       GRX(I,J,K) = Yf(K)/8.
	    ENDDO
	    DO K=1,NROW
	       Zf(K) =GVX(I,J,K)
	    ENDDO
	    CALL FOUR67(4,IFOUR)
	    DO K=1,NROW
	       GVX(I,J,K) = Yf(K)/8.
	    ENDDO
	   ENDDO
      ENDDO
				       write (*,*) ' x is done'
C				       Y-component
      DO K=1,NROW
	   DO J=1,NROW
	    DO I=1,NROW
	       Zf(I) =GRY(I,J,K)
	    ENDDO
	    CALL FOUR67(4,IFOUR)
	    DO I=1,NROW
	       GRY(I,J,K) = Yf(I)
	    ENDDO
	    DO I=1,NROW
	       Zf(I) =GVY(I,J,K)
	    ENDDO
	    CALL FOUR67(4,IFOUR)
	    DO I=1,NROW
	       GVY(I,J,K) = Yf(I)
	    ENDDO
	   ENDDO
      ENDDO
      DO K=1,NROW
	   DO I=1,NROW
	    DO J=1,NROW
	       Zf(J) =GRY(I,J,K)
	    ENDDO
	    CALL FOUR67(4,IFOUR)
	    DO J=1,NROW
	       GRY(I,J,K) = Yf(J)
	    ENDDO
	    DO J=1,NROW
	       Zf(J) =GVY(I,J,K)
	    ENDDO
	    CALL FOUR67(4,IFOUR)
	    DO J=1,NROW
	       GVY(I,J,K) = Yf(J)
	    ENDDO
	   ENDDO
      ENDDO
      DO J=1,NROW
	   DO I=1,NROW
	    DO K=1,NROW
	       Zf(K) =GRY(I,J,K)
	    ENDDO
	    CALL FOUR67(4,IFOUR)
	    DO K=1,NROW
	       GRY(I,J,K) = Yf(K)/8.
	    ENDDO
	    DO K=1,NROW
	       Zf(K) =GVY(I,J,K)
	    ENDDO
	    CALL FOUR67(4,IFOUR)
	    DO K=1,NROW
	       GVY(I,J,K) = Yf(K)/8.
	    ENDDO
	   ENDDO
      ENDDO
				       write (*,*) ' y is done'
C				       Z-component
      DO K=1,NROW
	   DO J=1,NROW
	    DO I=1,NROW
	       Zf(I) =GRZ(I,J,K)
	    ENDDO
	    CALL FOUR67(4,IFOUR)
	    DO I=1,NROW
	       GRZ(I,J,K) = Yf(I)
	    ENDDO
	    DO I=1,NROW
	       Zf(I) =GVZ(I,J,K)
	    ENDDO
	    CALL FOUR67(4,IFOUR)
	    DO I=1,NROW
	       GVZ(I,J,K) = Yf(I)
	    ENDDO
	   ENDDO
      ENDDO
      DO K=1,NROW
	   DO I=1,NROW
	    DO J=1,NROW
	       Zf(J) =GRZ(I,J,K)
	    ENDDO
	    CALL FOUR67(4,IFOUR)
	    DO J=1,NROW
	       GRZ(I,J,K) = Yf(J)
	    ENDDO
	    DO J=1,NROW
	       Zf(J) =GVZ(I,J,K)
	    ENDDO
	    CALL FOUR67(4,IFOUR)
	    DO J=1,NROW
	       GVZ(I,J,K) = Yf(J)
	    ENDDO
	   ENDDO
      ENDDO
      DO J=1,NROW
	   DO I=1,NROW
	    DO K=1,NROW
	       Zf(K) =GRZ(I,J,K)
	    ENDDO
	    CALL FOUR67(4,IFOUR)
	    DO K=1,NROW
	       GRZ(I,J,K) = Yf(K)/8.
	    ENDDO
	    DO K=1,NROW
	       Zf(K) =GVZ(I,J,K)
	    ENDDO
	    CALL FOUR67(4,IFOUR)
	    DO K=1,NROW
	       GVZ(I,J,K) = Yf(K)/8.
	    ENDDO
	   ENDDO
      ENDDO
				       write (*,*) ' z is done'
      RETURN
      END

C_________________________________________ Fermi distribution
C                                y - random number (0-1)
C                            dn = v^2dv/(exp(v)+1)/1.8036
C                                 \int_0^\infty dn =1      
      FUNCTION Vfermi(y)
         z  = y**0.333
         d  = MAX(1.-y,1.e-6)
            Vfermi = -1.15*log(d) +3.7e-4+z*(2.219+0.396*z) 
         IF( y.gt. 0.35)Then
            IF(y .gt. 0.55)Then
               Vfermi =Vfermi+2.04-1.937*sqrt(d)-1.24*y
            Else
               Vfermi =Vfermi+(2.04-1.937*sqrt(d)-1.24*y)*
     .                        (y-0.35)/0.2               
            EndIf
         Endif
      RETURN
      END
