! _______________________ START 3-D SIMULATIONS                         
!                                                                       
!                         Klypin, August 1993                           
!                                                                       
!-------------------------------------------------------------
MODULE Random
Contains
!--------------------------------------                                 
!				normal random numbers                                              
      FUNCTION GAUSS (M) 
!--------------------------------------                                 
      X = 0. 
      DO I = 1, 5 
      X = X + RANDd (M) 
      EndDo 
      X2 = 1.5491933 * (X - 2.5) 
      GAUSS = X2 * (1. - 0.01 * (3. - X2 * X2) ) 
      RETURN 
      END FUNCTION GAUSS                            

!------------------------------------------------                       
!				                                       random number generator     
      FUNCTION RANDd (M) 
!------------------------------------------------                       
      DATA LC, AM, KI, K1, K2, K3, K4, L1, L2, L3, L4	 / 453815927,     &
      2147483648., 2147483647, 536870912, 131072, 256, 	16777216, 4,    &
      16384, 8388608, 128 /
      ML = M / K1 * K1 
      M1 = (M - ML) * L1 
      ML = M / K2 * K2 
      M2 = (M - ML) * L2 
      ML = M / K3 * K3 
      M3 = (M - ML) * L3 
      ML = M / K4 * K4 
      M4 = (M - ML) * L4 
      M5 = KI - M 
      IF (M1.GE.M5) M1 = M1 - KI - 1 
      ML = M + M1 
      M5 = KI - ML 
      IF (M2.GE.M5) M2 = M2 - KI - 1 
      ML = ML + M2 
      M5 = KI - ML 
      IF (M3.GE.M5) M3 = M3 - KI - 1 
      ML = ML + M3 
      M5 = KI - ML 
      IF (M4.GE.M5) M4 = M4 - KI - 1 
      ML = ML + M4 
      M5 = KI - ML 
      IF (LC.GE.M5) ML = ML - KI - 1 
      M = ML + LC 
      RANDd = M / AM 
      RETURN 
      END FUNCTION RANDd                            
end module Random
!
!
!------------------------------------------------------------------
MODULE   setInitialConditions
PUBLIC
Character (LEN=12), DIMENSION(2) :: NamesDir=(/ &
       'PMstartX.DAT',&                                                                 
       'PMstartY.DAT'/)                                                                 

INCLUDE 'PMinitialp.h'
      INTEGER, PARAMETER :: NPAGE = NROW**2,  & ! # particles in a record
                            NMAX  = NGRID/2,     & 
                            NRECL = NPAGE*6,   & ! # particles in a record
                            NARR  = MAX(2*NROW+1,NGRID+1), & ! FFT
                            NF67  = MAX(NROW,NGRID/2), &  
                            Nblocks=NROW/Lblock,       &
                            NSPEC = NROW/2,            &  ! No waves shorter than Ny
                            marr  = NROW+1,            &  ! Arrays for FFT
                            mf67  = NROW/2

      REAL, DIMENSION(10) ::   timing =0.
      COMMON / CONTROL/ AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & 
                       TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0, &
                       NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5 &
                       ,Ocurv,extras(100)
      COMMON / HEADDR/  HEADER
      CHARACTER*45      HEADER
      COMMON /FOURAR/Zf(NARR),Yf(NARR)
      COMMON /F67COM/ &
                      IBC,      IP,       ISL,     L1,     N2, &
                      N3,       N4,        N7, &
                      SI(NF67),    INDEX(NF67)

      COMMON / ROW /	XPAR(NPAGE),YPAR(NPAGE),ZPAR(NPAGE),&
     			VX(NPAGE),VY(NPAGE),VZ(NPAGE)
      DIMENSION     RECDAT(NRECL),wspecies(10),lspecies(10)
      EQUIVALENCE    (RECDAT(1),XPAR(1)) &
                                    ,(wspecies(1),extras(1)), &
                                    (lspecies(1),extras(11))
      COMMON / MASK/	iMask(Nblocks,Nblocks,Nblocks) ! Refinement Mask
      REAL, ALLOCATABLE  :: GR(:,:,:)
      COMMON / BUFFx /	XPt(Nmaxpart)
      COMMON / BUFFv/   VXt(Nmaxpart)
      COMMON /SeedsPage/ NseedP(NROW)
      COMMON / LevelsM/ Max_lev_abs,Nmin_lev,Nmax_lev

Contains

!--------------------------------------------------                     
!                    : Store and retrieve seeds for parallelization 
!                      of random number generator
!--------------------------------------------------                     
      SUBROUTINE SetRandom(iCont)
use     Random

INCLUDE 'luxuryp.h' 
Character (LEN=70) :: FileName
Logical            :: FileExists =.false., RNGluxury =.false. 
Logical            :: iCont
REAL               ::    t0=0,t1=0.        ! counters for timing

      CALL CPU_TIME(t0)                    ! initialize time counter
      write(FileName,'(a,i1.1,".",i4.4,".",i8.8,".dat")')        &
                     'Seeds.',INT(extras(80)),NROW,Nseed
INQUIRE(file=FileName,EXIST=FileExists)

If(abs(extras(80)-1.) < 1.e-5)RNGluxury =.true.

If(FileExists==.true.)Then            ! read the data
      write(*,'(10x,2a)')'Open Existing FileName=',FileName
      Open(62,file=FileName,form='unformatted')
   read(62)NROWdummy,Nseeddummy
   if(NROWdummy /= NROW .or. Nseeddummy /= Nseed)Then
      write(*,*) ' File with seeds does not match parameters'
      write(*,*) ' Read Nrow =',NROWdummy,' Internal=',NROW
      write(*,*) ' Read Nseed=',Nseeddummy,' Internal=',Nseed
      STOP
   endIf
   If(RNGluxury)Then ! 
      read(62)i24A,j24A,in24A,kountA,carryA,gsetA,iFlagA
      read(62)SeedsPage
      Close(62)
   Else
      read(62)i24A
   EndIf
   Close(62)
Else
                       ! open file with seeds 
      write(*,'(10x,2a)')'Create FileName=',FileName
      Open(62,file=FileName,form='unformatted')
                                  ! generate seeds 
   If(RNGluxury)Then !  Detect random generator: luxury or nrandd 
      iFlag = 0
      gSet  = 0.
      Do k=1,NROW                               ! luxury
          i24A(k)    = i24
          j24A(k)    = j24
          in24A(k)  = in24
          kountA(k) = kount
          carryA(k) = carry
          gsetA(k)  = gSet
          iFlagA(k) = iFlag
             Do m=1,24
                 SeedsPage(m,k) = seeds(m)
             EndDo    
             Do j=1,NROW
                Do i=1,NROW
                  If(i+j+k.ne.3)Then
                    x = GAUSS3(gSet,iFlag)
                  EndIf
                EndDo
             EndDo
      EndDo
      write(62)NROW,Nseed
      write(62)i24A,j24A,in24A,kountA,carryA,gsetA,iFlagA
      write(62)SeedsPage
      Close(62)
   Else 
   write(*,*) ' Nseed=',Nseed
       NRAND = Nseed                  
      Do k=1,NROW                               ! gauss and randd
             i24A(k)  = NRAND
             Do j=1,NROW
                Do i=1,NROW
                   If(i+j+k.ne.3)Then
                      x = GAUSS(NRAND)
                   EndIf
                EndDo
             EndDo
      EndDo
      write(62)NROW,Nseed
      write(62)i24A
      Close(62)

   EndIf                                       ! end luxury/randd switch
      write(*,'(/2(20("-"),a/))')'File with seeds was created. STOP',  &
         'Do not change anything. Just restart the code to produce Initial Conditions'
      iCont = .false.
EndIf                                          ! file exists switch
      CALL CPU_TIME(t1)
      timing(2) = t1-t0             ! time for setting random numbers 

        Return
        End SUBROUTINE SetRandom
!--------------------------------------------------                     
!                                                     sqrt(Power spectrum)
!                                                           k = (2pi/L) 
!--------------------------------------------------                     
      FUNCTION TRUNF (WK) 
!-------------------------------------------------
      PARAMETER (NtabM = 100000) 
      COMMON / PTABLE / xkt (0:NtabM), Pkt (0:NtabM), Ntab, StepK,      &
                                           alog0, iFlag                       
      COMMON / TRUNCOM / Om, Omb, Omc, Omnu, Par (6), ns, qqscaleb,     &
                                           QSCALE
      Real ns, k 
                                                                        
      IF (WK.GE.FLOAT (NSPEC) ) THEN 
         TRUNF = 0. 
         RETURN 
      ENDIF 
      k = QSCALE * wk 
                                                                        
                               ! use approximations                     
      If (iFlag.gt.0) Then 
                             ! Holtzman approx                          
      If (Par (6) .ne.0.) Then 
      sk = sqrt (k) 
      TRUNF = k** (ns / 2.) / (1. + sk * (Par (2) + sk * (Par (3)       &
      + sk * (Par (4) + sk * (Par (5) ) ) ) ) ) ** (Par (6) )           
                             ! BBKS + Sugiyama approx                   
      Else 
!        Gammaeff =Om*hsmall/exp(Omb*(1.+sqrt(hsmall/0.5)/Om))          
!        Q = wk/hsmall/Gammaeff                                         
      Q = k * qqscaleb 
      TRUNF = k** (ns / 2.) * LOG (1. + Par (1) * Q) / (Par (1) * Q)    &
      / sqrt (sqrt (1. + Q * (Par (2) + Q * (Par (3) **2 + Q * (Par (4) &
      **3 + Q * (Par (5) **4) ) ) ) ) )                                 
      EndIf 
                         ! use table for P(k)                           
      Else 
         TRUNF = sqrt (Ppk (k) ) 
      EndIf 
      RETURN 
      END FUNCTION TRUNF                            
!---------------------------------------                                
      FUNCTION Ppk (x) 
!                     interpolate table with p(k)                       
!                       x is in real 1/Mpc                              
!---------------------------------------                                
      PARAMETER (NtabM = 100000) 
      COMMON / PTABLE / xkt (0:NtabM), Pkt (0:NtabM), Ntab, StepK,      &
      alog0, iFlag                                                      
                                                                        
                              ! slope is ns =-3                         
      If (x.ge.xkt (Ntab) ) Then 
         Ppk = Pkt (Ntab) / (x / xkt (Ntab) ) **3 
         Return 
      EndIf 
                                         
      If (x.lt.xkt (1) ) Then    ! slope is ns=1                
         Ppk = Pkt (1) * (x / xkt (1) ) 
         Return 
      EndIf 
      ind = INT ( (log10 (x) - alog0) / StepK) + 1 
      dk = xkt (ind+1) - xkt (ind) 
      Ppk = (Pkt (ind) * (xkt (ind+1) - x) + Pkt (ind+1) * (x - xkt (   &
      ind) ) ) / dk                                                     
!        write (*,'(5x," k=",g12.4," table=",2g12.4,"  P=",3g12.4)')    
!     &    x,xkt(ind),xkt(ind+1),                                       
!     &    Ppk,Pkt(ind),Pkt(ind+1)                                      
      Return 
      END FUNCTION Ppk                              
!---------------------------------------                                
      SUBROUTINE ReadTable (hsmall) 
!                     interpolate table with p(k)                       
!---------------------------------------                                
      PARAMETER (NtabM = 100000) 
      COMMON / PTABLE / xkt (0:NtabM), Pkt (0:NtabM), Ntab, StepK,      &
      alog0, iFlag                                                      
      CHARACTER(80) Line 
                                                                        
                           ! W.Hu table                                 
      If (iFlag.eq. - 1) Then 
         Open (50, file = 'WHUpk.dat') 
         twop = 2. * 3.14145926**2 
         Ntab = 0 
   10    READ (50, *, end = 30, err = 30) xx, pp 
            Ntab = Ntab + 1 
            xkt (Ntab) = xx * hsmall     ! scale to real Mpc      
            Pkt (Ntab) = pp / xx**3 * twop ! k^3Pk -> Pk                
         GOTO 10 
   30    Write (*,*) 
         Write (*,*) ' Read ', Ntab, ' lines from file WHUpk.dat' 
         Write (*,*) ' hsmall =', hsmall 
         Goto 50
      EndIf 
                           ! Eis.Hu table                               
      If (iFlag.eq. - 2) Then 
         Open (50, file = 'pk_EHu_Om0.3_Omb0.043_s8=0.9.dat') 
      Else 
         Open (50, file = 'pk_EHu_WMAP06.dat') 
      EndIf 
                                 ! read header                          
      Do i = 1, 9 
         READ (50, '(a)', end = 100, err = 100) Line 
      EndDo 
      Ntab = 0 
   12 READ (50, *, end = 32, err = 32) xx, pp 
      Ntab = Ntab + 1 
      xkt (Ntab) = xx * hsmall   ! scale to real Mpc      
      Pkt (Ntab) = pp                !  Pk 
      GOTO 12 
   32 Write (*,*) 
      Write (*,*) ' Read ', Ntab, ' lines from file pk_EHu_Om0.3...'
      Write (*,*) ' hsmall =', hsmall 
                                                                        
                                                                        
   50 If (Ntab.le.1) stop'wrong table for p(k)' 
      StepK = log10 (xkt (2) / xkt (1) ) 
      alog0 = log10 (xkt (1) ) 
                                                                        
                            ! test that spacing is the same             
      Do k = 2, Ntab - 1 
        ss = log10 (xkt (k + 1) / xkt (k) ) 
        If (abs (ss / StepK - 1.) .gt.2.e-2) Then 
        Write (*,*) ' error in K spacing. k=', k, xkt (k + 1),xkt (k)
        STOP 
      EndIf 
      EndDo 
      CLOSE (50) 
      Return 
  100 Write (*,*) ' Error reading table of P(k)' 
      STOP 
                                                                        
      END SUBROUTINE ReadTable                      
                                                                        
                                                                        
!*********************************************************************  
!	   Power spectrum for CDM. Wk-  wavenumber(Mpc^-1, h=0.5)             
!                         P = wk *T^2(k)                                
!                         QSCALE = 2pi/Box_size,  QS =hubble**(-2)      
!      FUNCTION TRUNF(WK)                                               
!      COMMON / TRUNCOM/ QSCALE,g,gy                                    
!      COMMON / TRUNCDM/ QS,A1,A2,A3,A4,A5,A32,A43,A54                  
!       Q = QSCALE*QS*WK                                                
!       TRUNF = sqrt( WK*                                               
!     +		 (LOG(1.+A1*Q)/(A1*Q))**2/                                     
!     +		 SQRT(1.+Q*(A2 +Q*(A32   +Q*(A43+Q*A54)))) )                   
!                                                                       
!      RETURN                                                           
!      END                                                              
                                                                        
!*********************************************************************  
!			  INITIALIZE CONSTANTS:                                             
!			      Scalel    =length of square in MPC                            
!			      AEXPN =expansion factor (current)                             
!			      ASTEP	  =increment of expansion factor: AEXPN =AEXP_0+N*ASTEP 
!                     PARTW	=weight of a particle: total 'mass' per cell
!                                          even for LCDM or OCDM        
!			      AMPLT	=amplitude of perturbations inside the simulation box   
!			      NBYTE	=record length in bytes                                 
!                                                                       
      SUBROUTINE InitValues (NBYTE, SCALEL) 
!--------------------------------------------
      PARAMETER (NtabM = 100000) 
      COMMON / PTABLE / xkt (0:NtabM), Pkt (0:NtabM), Ntab, StepK,      &
                                           alog0, iFlag                                                      
      COMMON / TRUNCOM / Om, Omb, Omc, Omnu, Par (6), ns, qqscaleb,     &
                                                QSCALE                                                            
      COMMON / FERMI / Vscale 
      Real           ns 
      Character Answer*1 
      Real          INPUT 
      DATA PI  / 3.1415926535 / 
                                                                        
      Do i = 1, 100 
         extras (i) = 0. 
      enddo 
      Write (*,  * ) 'Would you like to use a model provided in cdm.fi&
     &t '                                                               
      Write (*,*) 'or your own model?' 
      Write (*, '(A,$)') ' Enter Yes for cdm.fit;  No for your model =&
     &'                                                                 
      Read (*, '(A)') Answer 
      HEADER = 'N=128x256 L=20h-1CDMsig=0.239 z=30-----------' 
      Write (*,  * ) '------ Enter Header for the run up to 45 characters'
      Write (*,  * ) '       Example:' 
      Write (*, '(A)') HEADER 
      READ (*, '(A)') HEADER 
      Write (*, '(A)') HEADER 
      AEXPN = INPUT (' Initial expansion parameter (0-1)=') 
      ASTEP = INPUT (' Step of the expansion parameter  =') 
      AMPLT = INPUT (' Amplitude of density fluctuations=') 
      ns = INPUT (' Slope of the Power spectrum     n=') 
      SCALEL = INPUT (' Box size (Mpc/h)                 =') 
      Ocurv = INPUT (' Omega_curvature   at present     =') 
                                                                        
      Nseed = INPUT (' Random seed (Integer  1-2^31)    =') 
      If (Nseed.le.0) Then 
         lux             = 2                              ! level for luxury 
         extras (80) = 1. 
         extras (81) = lux 
         If (Nseed.eq.0) Nseed = 121071 
         Nseed        = ABS (Nseed) 
         CALL rluxgo (lux, Nseed, 0, 0)      ! initialize luxury 
      Else                                                    ! use randd and gauss
          extras (80) = 0. 
      EndIf 
!                       extras(80) is used as switch for                
!                                        random numbers generators      
!                         = 0 - old Randd+Gauss                         
!                         = 1 - luxury+Gauss3                           
!                        extras(81) = lux - level of luxury             
                                                                        
                                              ! store the box size      
      extras (100) = SCALEL 
      If (Answer.eq.'Y'.or.Answer.eq.'y') Then 
         CALL MODEL (hsmall) 
         hubble = hsmall 
         Om0 = Om 
         Oml0 = 1. - Om0 - Ocurv 
      Else 
         Write (*,*) ' You use your cosmological model.' 
         Write (*,*) ' Be sure you provide routine TRUNF' 
         hubble = INPUT (' The Hubble constant (h=H/100)    =') 
         Om0 = INPUT (' Omega_0 matter    at present     =') 
         Oml0 = INPUT (' Omega_lambda      at present     =') 
      EndIf 
                              ! scale it to real megaparsecs            
      SCALEL = SCALEL / hubble 
      NBYTE = NPAGE * 6 * 4 
                                                                        
      Write (*, '(A,$)') ' Enter Min and Max number of mass species =' 
      READ ( *, * ) Nmin_lev, Nmax_lev 
      W = (FLOAT (NGRID) / FLOAT (NROW) ) **3 
      PARTW = W 
                                       ! multiple mass resolution       
      If (Nmin_lev.lt.Nmax_lev) Then 
         Write ( *, * ) 
         Write (*,*) ' You use multiple mass resolution' 
      ENDIF 
                                    ! change the order of levels        
       If(Nmin_lev.gt.Nmax_lev)Then                                     
          i = Nmax_lev                                                  
          Nmax_lev = Nmin_lev                                           
          Nmin_lev = i                                                  
       endif                                                            
                                                                        
       ISTEP = 0                                                        
       TINTG = 0.                                                       
       AU0   = 0.                                                       
       AEU0  = 0.                                                       
       EKIN  = 0.                                                       
       EKIN1 = 0.                                                       
       EKIN2 = 0.                                                       
       NROWC = NROW                                                     
       NGRIDC= NGRID                                                    
           QSCALE = 2.*PI/SCALEL                                        
           QS     = hubble**(-2)                                        
        write (*,'(/3x,3(a,g12.3))') ' Qscale=',QSCALE,                &
                                     ' Box=',SCALEL,' slope=',ns                         
                                                                        
      RETURN                                                            
      END  SUBROUTINE InitValues                                         
!________________________________________Read parameters of the model   
!                                 Om       = Omega_0                    
!                                 Omb     = Omega_baryon                
!                                 Omnu  = Omega_neutrino                
!                                 hsmall  = hubble constant             
!                                 Par      = fitting parameters         
!                                 Par(6) = 0  --- bbks-style+Hu&Sugiyama
!                                 Par(6) ne 0 --- holtzman-style        
      SUBROUTINE MODEL (hsmall) 
!---------------------------------------                                
      COMMON / TRUNCOM / Om, Omb, Omc, Omnu, Par (6), ns, qqscaleb,     &
      QSCALE                                                            
      PARAMETER (NtabM = 100000) 
      COMMON / PTABLE / xkt (0:NtabM), Pkt (0:NtabM), Ntab, StepK,      &
      alog0, iFlag                                                      
      Real ns, INPUT 
      Character Header1*79 
                                                                        
      Line1 = INPUT (' Enter Line Number in cdm.fit     =') 
      If (Line1.gt.0) Then                 ! use analytical approximation
         iFlag = Line1 
         OPEN (2, file = 'cdm.fit') 
         Read (2, '(A)') Header1 
         If (Line1.gt.2) Then 
            Do i = 1, Line1 - 2 
               Read (2, * ) a 
            EndDo 
         EndIf 
         Read (2, * ) Om, Omb, Omc, Omnu, hsmall, Par 
         CLOSE (2) 
                                                                           
                            ! = T_cmb/2.7                                  
         theta = 2.726 / 2.7 
                            ! Hu&Sugiyama fitting                          
         Ob0 = Omb / Om 
         Omh2 = Om * hsmall**2 
         a1 = (46.9 * Omh2) **0.670 * (1. + (32.1 * Omh2) ** ( - 0.532) ) 
         a2 = (12.0 * Omh2) **0.424 * (1. + (45. * Omh2) ** ( - 0.582) ) 
         alpha = 1. / a1**Ob0 / a2** (Ob0**3) 
         qqscaleb = theta**2 / (1. - Ob0) **0.60 / sqrt (alpha) / Omh2 
                                                                           
         Write ( *, 20) Om, Omb, Omc, Omnu, hsmall, Par 
!         Write (16,20)Om,Omb,Omc,Omnu,hsmall,Par                          
   20    Format (' Model: Om_0=', F5.3, ' Om_baryon=', F5.3, ' Om_cold=',  &
                F5.3, ' O_nu=', F5.3, ' hsmall=', F4.2, / 8x, 'Parameters=',    &
               (6G9.3) )                                                         
                                  
      Else                      ! read spectrum from table               
         iFlag = Line1 
         If (iFlag.lt. - 3) Stop' Error: model Line # less than -3' 
                                 ! W.Hu table                              
         If (iFlag.eq. - 1) Then 
            hsmall = 0.7 
            Omb = 0.024 / hsmall**2 
            Om = 0.3 
            Ocurv = 0. 
            Omc = Om 
         EndIf 
                                 ! Eis.Hu approximation                    
         If (iFlag.eq. - 2) Then 
            hsmall = 0.7 
            Omb = 0.04 
            Om = 0.3 
            Ocurv = 0. 
            Omc = Om 
         EndIf 
                                 ! Eis.Hu approximation: WMAP              
         If (iFlag.eq. - 3) Then 
            hsmall = 0.73 
            Omb = 0.04 
            Om = 0.24 
            Ocurv = 0. 
            Omc = Om 
         EndIf 
                                                                              
         CALL ReadTable (hsmall) 
                                                                           
         Write ( *, 22) Om, Omb, hsmall 
         Write (1, 22) Om, Omb, hsmall 
   22    Format (' Model: Om_0=', F5.3, ' Om_baryon=', F5.3, ' hsmall=',   &
         F4.2, ' <== Using table for P(k)')                                
      EndIf 
                                                                        
      Return 
      END SUBROUTINE MODEL                          
!------------------------------------------------                       
!                              Clear MASK for multiple mass resolution  
!                               Initiat to level Nmin_lev               
      SUBROUTINE ZeroMask 
!------------------------------------------------                       
!      INCLUDE 'PMinitialp.h'	 
                                                                        
      DO MK3 = 1, Nblocks 
      DO MK2 = 1, Nblocks 
      DO MK1 = 1, Nblocks 
         iMask (MK1, MK2, MK3) = Nmin_lev 
      ENDDO 
      ENDDO 
      ENDDO 
      RETURN 
      END SUBROUTINE ZeroMask                       
!------------------------------------------------                       
!                       Set MASK for multiple mass resolution           
!                               Read particles with lagrangian coordinat
!                                                to set the mask        
!                               Block   = 1 = low resolution            
!                                       = 2 = higher resolution         
!                                       = 3 = even higher resolution    
!               After the mask is set, make maximum resolution Nmax_lev 
!                                                                       
      SUBROUTINE SetMask 
!------------------------------------------------                       
      Dimension nLevel (10) 
      Character *70 Line 
                                                                        
      Do i = 1, 5 
         Read (30, '(A)') Line 
         Write (*, '(A)') Line 
      EnDDo 
      Nlines = 0 
   10 Read (30, *, end = 50, err = 50) i, xw, yw, zw, vwx, vwy, vwz, im,&
      jm, km                                                            
      Nlines = Nlines + 1 
!           write (*,*) ' Block=',im,jm,km                              
      If (im.lt.0.or.im.gt.Nblocks) Then 
         Write (*,*) ' wrong x_block: ', im, jm, km, ' Line=', Nlines 
         Stop 
      EndIf 
      If (jm.lt.0.or.jm.gt.Nblocks) Then 
         Write (*,*) ' wrong y_block: ', im, jm, km, ' Line=', Nlines 
         Stop 
      EndIf 
      If (km.lt.0.or.km.gt.Nblocks) Then 
         Write (*,*) ' wrong z_block: ', im, jm, km, ' Line=', Nlines 
         Stop 
      EndIf 
!             If(iMask(im,jm,km).ne.0)Then                              
!              write (*,*) ' Occupied block: ',im,jm,km,' Line=',Nlines 
!              write (*,*) '                 Block=',iMask(im,jm,km)    
!              Stop                                                     
!           EndIf                                                       
                                                                        
                                        ! mark high resolution block    
      iMask (im, jm, km) = Nmax_lev 
      Do ik = - 1, 1 
         k = km + ik 
         If (k.le.0) k = k + Nblocks 
         If (k.gt.Nblocks) k = k - Nblocks 
         Do ij = - 1, 1 
            j = jm + ij 
            If (j.le.0) j = j + Nblocks 
            If (j.gt.Nblocks) j = j - Nblocks 
            Do ii = - 1, 1 
               i = im + ii 
               If (i.le.0) i = i + Nblocks 
               If (i.gt.Nblocks) i = i - Nblocks 
                                                            
               ind = abs (ik) + abs (ij) + abs (ii) ! only seven blocks 
                                                                   ! are marked out                        
!                               If(iMask(i,j,k).le.0.or.iMask(i,j,k).gt.Nspecies)
!              +                     Then                                        
!                                  Stop                                          
!                               Endif                                            
               If (ind.le.1) Then 
                   iMask (i, j, k) = Nmax_lev ! change mask           
               Endif 
            ENDDO 
         ENDDO 
      ENDDO 
      GoTo 10 

   50 Write (*,*) ' Total Lines read=', Nlines 
!                      Change the mask: iteratevely add layers of lower resolution. 
!                       Find blocks adjasent to filled blocks and
!                       mark them with kSpec =Nspecies-1, -2,
!                      
                                     !                                  
   DO kSpec = Nmax_lev - 1, 2, - 1 
      DO MK3 = 1, Nblocks 
      DO MK2 = 1, Nblocks 
      DO MK1 = 1, Nblocks 
         iM = iMask (MK1, MK2, MK3) 
         If (iM.eq.kSpec + 1) Then           !  change to kSpec if not high res
            Do ik = - 1, 1 
               k = MK3 + ik 
               If (k.le.0) k = k + Nblocks 
               If (k.gt.Nblocks) k = k - Nblocks 
                  Do ij = - 1, 1 
                     j = MK2 + ij 
                     If (j.le.0) j = j + Nblocks 
                     If (j.gt.Nblocks) j = j - Nblocks 
                     Do ii = - 1, 1 
                        i = MK1 + ii 
                        If (i.le.0) i = i + Nblocks 
                        If (i.gt.Nblocks) i = i - Nblocks 
                                                                                   ! chang
                        If (iMask (i, j, k) .lt.kSpec) iMask (i, j, k) = kSpec 
                     ENDDO 
                  ENDDO 
            ENDDO 
         EndIf                   ! end iM=kSpec+1 test
      ENDDO                  ! end MK1
      ENDDO                  ! end MK2
      ENDDO                  ! end MK3
  Enddo                        ! end kSpec 
                                                                        
      Do i = 1, Max_lev_abs  ! count marked blocks              
         nLevel (i) = 0 
      Enddo 
      DO MK3 = 1, Nblocks 
      DO MK2 = 1, Nblocks 
      DO MK1 = 1, Nblocks 
         iM = iMask (MK1, MK2, MK3) 
         If (iM.le.Nmax_lev.and.iM.ge.Nmin_lev) Then 
            nLevel (iM) = nLevel (iM) + 1 
         Else 
            Write (*,*) 'wrong mask: ', iM, ' block=', MK1, MK2, MK3 
            Write (*,  * ) '      it should be within limits:', Nmin_lev, Nmax_lev
            Stop 
         EndIf 
      ENDDO 
      ENDDO 
      ENDDO 
                                                                        
      Ntot = 0 
      Do i = Max_lev_abs, 1, - 1 
         icurrent = nLevel (i) * 8** (i - 1) 
         Write (*,*) ' Level=', i, ' Number of Blocks=', nLevel (i) ,  &
         ' N_Particles=', icurrent                                         
         Ntot = Ntot + icurrent 
      Enddo 
      Write (*,*) ' Total expected particles =', Ntot 
      RETURN 
      END SUBROUTINE SetMask                        
!------------------------------------------------                       
!                             Make a realization of spectrum of perturbations
!                             ALPHA  = normalization factor for displacement
!                             NRAND = seed for random numbers           
      SUBROUTINE SPECTR (ALPHA, NRAND, iDirection) 
!------------------------------------------------                       
use Random

INCLUDE 'luxuryp.h' 
      REAL(8)   :: SUMM=0., Wi3, Wj3, Wk3, WD, TS, TRX 
      REAL      ::    t0=0,t1=0.        ! counters for timing
      INTEGER, DIMENSION(NROW) :: mapz, map3

      CALL CPU_TIME(t0)                    ! initialize time counter
      iRandom = INT (extras (80) + 1.e-5) ! set type of Random Numbers       
                                          ! 0 - old GAUSS 1-GAUSS3(ranlux)
      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 
         mapz(j) = jz
         map3(j) = j3*jsign  
      EndDo
!.... Open_MP                                                           
!$OMP PARALLEL DO DEFAULT(SHARED)  &                                     
!$OMP PRIVATE(Mk3,Mk2,Mk1)                                            
      DO Mk3 = 1, NROW 
      DO Mk2 = 1, NROW 
      DO Mk1 = 1, NROW 
          GR (Mk1, Mk2, Mk3) = 0. 
      ENDDO 
      ENDDO 
      ENDDO 
!                                      Set Spectrum                     
!.... Open_MP                                                           
!$OMP PARALLEL DO DEFAULT(SHARED)  &                                     
!$OMP PRIVATE(k,ksign,kz,k3,Wk3,j,jsign,jz,j3,Wj3,i,isign,iz,i3,Wi3)  &
!$OMP PRIVATE(WD,Wk,TS,TRX,m,NRAND,gSet,iFlag)           &
!$OMP REDUCTION(+:SUMM)
      DO k = 1, NROW 
         If (iRandom.eq.0) Then 
            NRAND = i24A(k)
         Else 
           i24      = i24A(k) 
           j24      = j24A(k) 
           in24    = in24A(k) 
           kount   = kountA(k)
           mkount  = 0 
           carry  = carryA(k)
           iflag  = iFlagA(k)
           gset   = gSetA(k)
              Do m=1,24
                  seeds(m) =SeedsPage(m,k) 
              EndDo    
         EndIf         ! end iRandom

         Wk3 = map3(k)**2 

         DO j = 1, NROW 
            Wj3 = map3(j)**2 
            DO i = 1, NROW 
               Wi3 = map3(i)**2 
                                                                                 
               IF (i + j + k.EQ.3) THEN 
                  GR (1, 1, 1) = 0. 
               ELSE 
                  WD = Wi3 + Wj3 + Wk3 
                  Wk = SQRT (WD) 
                  If (iRandom.eq.0) Then 
                     TS = TRUNF (Wk) * GAUSS (NRAND) 
                  Else 
                     TS = TRUNF (Wk) * GAUSS3 (gSet,iFlag) 
                  EndIf 
                  TRX = TS / WD 
                                                                   ! start switch
                  If (iDirection.eq.1) Then 
                     GR (mapz(i), j, k) = TRX * map3(i) 
                  ElseIF (iDirection.eq.2) Then 
                     GR (i, mapz(j), k) = TRX * map3(j) 
                  Else 
                     GR (i, j, mapz(k)) = TRX * map3(k)
                  EndIf       ! end iDirection
                  SUMM = SUMM + TS**2 
               ENDIF        ! end kx=ky=kz=0 switch
            ENDDO         ! i
         ENDDO            ! j
      ENDDO               ! k
                                                                        
      IF (SUMM.LE.0.) Write (*,  * ) ' Error!!! Summ over spectrum = 0'
      ALPHA = AMPLT / SQRT (SUMM) * sqrt (8.) 
      Write (*,'(10x,a20,3g12.4)') 'SPECTR ==>  SUMM=', SUMM, ALPHA 
      CALL CPU_TIME(t1)
      timing(3) = timing(3) +t1-t0        ! time for setting SPECTR
                                                                        
      RETURN 
      END SUBROUTINE SPECTR                         
!------------------------------------------------                       
!		                        Define coordinates and velosities for        
!                                  particles with resolution = Indx     
!                                  Indx =1 = high resolution            
!                                  Indx =2 = medium resolution          
!                                  Icurrent = number of particles       
      SUBROUTINE BLOCKS (XCONS, VCONS, Indx, Icurrent, iDirection) 
!------------------------------------------------                       
      Dimension nLevel (10) 
      Real(8)       :: sDispl =0. 
      REAL          ::    t0=0,t1=0.        ! counters for timing
                                 
      CALL CPU_TIME(t0)                    ! initialize time counter   
      Do i = 1, Max_lev_abs   ! number of particles at           
         nLevel (i) = 0               ! different levels               
      Enddo 
      Boff = Lblock / 2 + 0.5 
      sDispl = 0. 
                                                                        
      QFACT = FLOAT (NGRID) / FLOAT (NROW) 
      DO Mk3 = 1, Nblocks 
         M3 = (Mk3 - 1) * Lblock       ! part of offset
      DO Mk2 = 1, Nblocks 
         M2 = (Mk2 - 1) * Lblock      ! part of offset
      DO Mk1 = 1, Nblocks 
         M1 = (Mk1 - 1) * Lblock     ! part of offset
         iM = iMask (Mk1, Mk2, Mk3) 
!                                             --------------------------
!                     M1 + Boff  = lagrangian coordinate of 'mean' particle
!                     Lblock**3 = total number of averaged cells        
         If (iM.eq.Indx) Then 

         If (iM.eq.1) Then        !-------------- Low resolution: one particle               
         kp = M3 + (Lblock + 1) / 2 
         jp = M2 + (Lblock + 1) / 2 
         ip = M1 + (Lblock + 1) / 2 
         DX = GR (ip, jp, kp) 
                                                                           
         Icurrent = Icurrent + 1 
         If (Icurrent.gt.Nmaxpart) Then 
             Write (*,*) 'Too many particles:', Icurrent, ' block=',      &
                                   Mk1,Mk2, Mk3                                                          
             STOP 
         EndIf 
         nLevel (Indx) = nLevel (Indx) + 1 
         If (iDirection.eq.1) Then 
            Q = QFACT * (M1 + Boff - 1.) + 1. 
         ElseIF (iDirection.eq.2) Then 
            Q = QFACT * (M2 + Boff - 1.) + 1. 
         Else 
            Q = QFACT * (M3 + Boff - 1.) + 1. 
         EndIf	 
         XPt (Icurrent)  = Q - XCONS * DX + 1.E-4 + 0.5 
         VXt (Icurrent) = VCONS * DX 
          sDispl            = sDispl + DX**2                                                                 
      Else 
                                      ! --------------- Highest resolution:             
      If (iM.eq.Max_lev_abs) Then 
      Do kp = M3 + 1, M3 + Lblock 
      Do jp = M2 + 1, M2 + Lblock 
      Do ip = M1 + 1, M1 + Lblock 
         DX = GR (ip, jp, kp)              !  find displacements 
         Icurrent = Icurrent + 1 
         If (Icurrent.gt.Nmaxpart) Then 
            Write (*,*) 'Too many particles:', Icurrent, ' block=',      &
                                 Mk1,Mk2, Mk3                                                          
            STOP 
         EndIf 
         nLevel (Indx) = nLevel (Indx) + 1 
         If (iDirection.eq.1) Then 
            Q = QFACT * (ip - 1.) + 1. 
         ElseIF (iDirection.eq.2) Then 
            Q = QFACT * (jp - 1.) + 1. 
         Else 
            Q = QFACT * (kp - 1.) + 1. 
         EndIf 
         XPt (Icurrent) = Q - XCONS * DX + 1.E-4 + 0.5 
         VXt (Icurrent) = VCONS * DX 
         sDispl = sDispl + DX**2
      EndDo        ! end ip
      EndDo        ! end jp
      EndDo        ! end kp
                                                                        
      Else 
            !-------------------            Medium resolution:                                         
            !                   KBlock = 2**(iM-1)        = number of particles in 1
            !                   Nbcells= Lblock/KBlock = number of cells per particl
            !                   (MK-1)*Lblock +(i-1)*Nbcells = offset of cells in 1D
            !                   (MK-1)*Lblock +(i-1)*Nbcells +Nbcells/2 +1/2 = lagra
            !                                                                 coordi
            !                    M1   = (Mk1-1)*Lbloc                               
      KBlock    = 2** (iM - 1)          ! number of particles in 1D per block
      Nbcells   = Lblock / KBlock   ! number of cells per particle in  1D
      Xlag        = Nbcells / 2 + 0.5   ! offset for lagrangian coordinates
      Nbcells3 = Nbcells**3            ! total number of cells per particle
                                                                        
      Do i = 1, KBlock                                ! loop over particles
         i1 = M1 + (i - 1) * Nbcells               ! offset for cells in X 
        Q1 = QFACT * (i1 + Xlag - 1.) + 1. ! scale it to box size 
           Do j = 1, KBlock 
               j1 = M2 + (j - 1) * Nbcells         ! offset for cells in Y
              Q2 = QFACT * (j1 + Xlag - 1.) + 1. 
              Do k = 1, KBlock 
                 k1 = M3 + (k - 1) * Nbcells     ! offset for cells in Z
                 Q3 = QFACT * (k1 + Xlag - 1.) + 1. 
                                                                                   
                 kp = k1 + (Nbcells + 1) / 2 
                 jp = j1 + (Nbcells + 1) / 2 
                 ip = i1 + (Nbcells + 1) / 2 
                 DX = GR (ip, jp, kp) 
                                                                                   
                 Icurrent = Icurrent + 1 
                 If (Icurrent.gt.Nmaxpart) Then 
                    Write (*,*) 'Too many particles:', Icurrent, ' block=', Mk1,  &
                    Mk2, Mk3                                                          
                    STOP 
                 EndIf 
                 nLevel (Indx) = nLevel (Indx) + 1 
                 If (iDirection.eq.1) Then 
                    Q = Q1 
                 ElseIF (iDirection.eq.2) Then 
                    Q = Q2 
                 Else 
                    Q = Q3 
                 EndIf	 
                 XPt (Icurrent) = Q - XCONS * DX + 1.E-4 + 0.5 
                 VXt (Icurrent) = VCONS * DX 
                 sDispl = sDispl + DX**2
            EndDo            ! end k
         EndDo               ! end j
      EndDo                  ! end i
      EndIf                    ! end iM=Max_lev_abs
      EndIf                    ! end iM = 1
      EndIf                    ! end iM = Indx  
      EndDo                 ! end Mk1 =1,Nblocks
      EndDo                 ! end Mk2 =1,Nblocks
      EndDo                 ! end Mk3 =1,Nblocks
                                                                        
      If (nLevel (indx) .gt.0) Then 
         Nspecies = Nspecies + 1 
         lspecies (Nspecies) = Icurrent 
         wspecies (Nspecies) = PARTW * Lblock**3 / 8** (Indx - 1) 
         Write (*, '(10x,a20,g12.4)') ' RMS 3d diplacement=',              &
                        XCONS*sqrt (sDispl/nLevel (indx) )   
         write (*,'(10x,a20,i4,a,i10,a,6i10)') ' Species=',Nspecies,   &
                         ' Current Particles=',Icurrent,' on  all levels:',         &
                          (nLevel(i),i=1,Max_lev_abs   )     
      EndIf 

      CALL CPU_TIME(t1)
      timing(4) = timing(4) +t1-t0        ! time for setting SPECTR
                                                                        
      RETURN                                                            
      END  SUBROUTINE BLOCKS                                         
!------------------------------------------------                       
!                                              				   FFT of the spectru
      SUBROUTINE VECTOR (IFOUR) 
!------------------------------------------------                       
      DIMENSION Uf (marr), Vf (marr) 
      DIMENSION Qi (mf67), jndx (mf67) 
      REAL               ::    t0=0,t1=0.        ! counters for timing
                          
      CALL CPU_TIME(t0)                                        
      jb1 = 4 
      jq = IFOUR 
      CALL SETF67 (jb1, jq, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx,&
      Uf, Vf)                                                           
                                                                        
!			                                                                    
!.... Open_MP                                                   -1-     
!$OMP PARALLEL DO DEFAULT(SHARED)   &                                    
!$OMP PRIVATE ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 )                          
      DO K = 1, NROW 
      DO J = 1, NROW 
      DO I = 1, NROW 
         Uf (I) = GR (I, J, K) 
      ENDDO 
      CALL FOUR67 (jb1, jq, jp, jsl, ll1, k2, k7, Qi, jndx, Uf, Vf) 
      DO I = 1, NROW 
         GR (I, J, K) = Vf (I) 
      ENDDO 
      ENDDO 
      ENDDO 
      Write (*,'(10x,a)') 'Xfft is done' 
!.... Open_MP                                                  -2-      
!$OMP PARALLEL DO DEFAULT(SHARED)   &                                    
!$OMP PRIVATE ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 )                          
      DO K = 1, NROW 
      DO I = 1, NROW 
      DO J = 1, NROW 
         Uf (J) = GR (I, J, K) 
      ENDDO 
      CALL FOUR67 (jb1, jq, jp, jsl, ll1, k2, k7, Qi, jndx, Uf, Vf) 
      DO J = 1, NROW 
         GR (I, J, K) = Vf (J) 
      ENDDO 
      ENDDO 
      ENDDO 
      Write (*,'(10x,a)') 'Yfft is done' 

      write (*,'(10x,a)') 'Swap k<->i'
!.... Open_MP
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE ( k,j,i,aa)
      DO J=1,NROW
      DO K=1,NROW-1
            DO I=K+1,NROW
               aa = GR(I,J,K)
               GR(I,J,K) =GR(K,J,I)
               GR(K,J,I) =aa
            ENDDO
         ENDDO
      ENDDO

!.... Open_MP                                                  -3-      
!$OMP PARALLEL DO DEFAULT(SHARED)   &                                    
!$OMP PRIVATE ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 )                          
      DO K = 1, NROW 
      DO J = 1, NROW 
      DO I = 1, NROW 
         Uf (I) = GR (I, J, K) 
      ENDDO 
      CALL FOUR67 (jb1, jq, jp, jsl, ll1, k2, k7, Qi, jndx, Uf, Vf) 
      DO I = 1, NROW 
         GR (I, J, K) = Vf (I) / 8. 
      ENDDO 
      ENDDO 
      ENDDO 

      write (*,'(10x,a)') 'Swap i<->k'
!.... Open_MP
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE ( k,j,i,aa)
      DO J=1,NROW
      DO K=1,NROW-1
            DO I=K+1,NROW
               aa = GR(I,J,K)
               GR(I,J,K) =GR(K,J,I)
               GR(K,J,I) =aa
            ENDDO
         ENDDO
      ENDDO

!!.... Open_MP                                                  -3-      
!!$OMP PARALLEL DO DEFAULT(SHARED)   &                                    
!!$OMP PRIVATE ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 )                          
!      DO J = 1, NROW 
!      DO I = 1, NROW 
!      DO K = 1, NROW 
!         Uf (K) = GR (I, J, K) 
!      ENDDO 
!      CALL FOUR67 (jb1, jq, jp, jsl, ll1, k2, k7, Qi, jndx, Uf, Vf) 
!      DO K = 1, NROW 
!         GR (I, J, K) = Vf (K) / 8. 
!      ENDDO 
!      ENDDO 
!      ENDDO 

      Write (*,'(10x,a)') ' FFT is done' 
      CALL CPU_TIME(t1)
      timing(5) = timing(5) +t1-t0        ! time for VECTOR

      RETURN 
      END SUBROUTINE VECTOR                         
!---------------------------------------------                          
!               1)    Update SKINE and Wtotal
!               2)    Write current data to  disk/tape                
!                   do not write files for 3rd direction
      SUBROUTINE WriteData (Vscale, AexpV, Wtotal,SKINE,iDirection) 
!----------------------------------------------                         
DOUBLE PRECISION    Wtotal,SKINE,vvx,vx2
REAL               ::    t0=0,t1=0.        ! counters for timing


!        write(30,*) ' XXX'       ! test data 
!        write(30,'(10g12.5)') (XPt(i),i=1,lspecies(Nspecies),5000)
!        write(30,*) ' VVV'
!        write(30,'(10g12.5)') (VXt(i),i=1,lspecies(Nspecies),5000)

      CALL CPU_TIME(t0)                    ! initialize time counter
      jstart = 1 
Do j = 1, Nspecies 
      vvx  = 0. 
      vx2  = 0. 
      jend = lspecies (j) 
      W    = wspecies (j) 
      If (jend.gt.jstart) Then 
            Do i = jstart, jend 
               vvx = vvx + VXt (i) 
               vx2 = vx2 + VXt (i) **2 
               SKINE = SKINE+W * VXt (i) **2
               Wtotal = Wtotal + w 
            EndDo                             ! end current species
         npp  = max (jend-jstart + 1, 1) 
         v2    = sqrt ( vx2 / npp) / AexpV * Vscale 
         vvx  = vvx / npp / AexpV * Vscale 
         Write ( *,'(22x,"V(km/s)=",2g11.3," Npart=",i10," Species=",i3)')   & 
                                        vvx, v2, jend-jstart + 1, j 
         jstart = jend+1 
      Endif 
EndDo                                       ! end all species
If(iDirection == 3)Then
      CALL CPU_TIME(t1)
      timing(6) = timing(6) +t1-t0        ! time for VECTOR
   Return                   ! do not write the last page
EndIf

      Ibuff = 0                              ! write data page by page
      KROW = 0 
Do i = 1, lspecies (Nspecies) 
        Ibuff = Ibuff + 1          ! buffer counter
        XPAR (Ibuff) = XPt (i) 
        VX (Ibuff)     = VXt (i) 
        If(Ibuff.ge.NPAGE)Then                                      
            KROW = KROW +1                                          
            write (*,*) ' Write page=',KROW,' i=',i,Ibuff           
            CALL WRIROW                                             
            Ibuff     =0                 ! initialize buffer counter                                           
         EndIf                                                      
EndDo                            ! end lspecies loop                             
          If(Ibuff.ne.0)Then          ! write leftover                                  
                KROW = KROW +1                                          
                write (*,*) ' Write page=',KROW,' i_buff=',Ibuff        
                CALL WRIROW                                             
                Ibuff     =0                                            
          EndIf                                                         
      RETURN                                                            
      CALL CPU_TIME(t1)
      timing(6) = timing(6) +t1-t0        ! time for VECTOR

      END  SUBROUTINE WriteData                                         
!--------------------------------------------------                     
!                    Write current  page of particles (x,v) to disk 
!                                                                       
      SUBROUTINE WRIROW 
         Write (1) XPAR 
         Write (1) VX
      RETURN 
      END SUBROUTINE WRIROW  
!--------------------------------------------------
!                             Write current PAGE of particles (x,v) to disk
!                             NRECL - length of ROW block in words
      SUBROUTINE WRIROWd(IROW)
        WRITE (21,REC=IROW) RECDAT
      RETURN
      END SUBROUTINE WRIROWd
!-------------------------------------------------------                
!                                                                       
!                                     Print statistics of spicies       
      Subroutine Species_Change (Wtotal) 

DOUBLE PRECISION Wtotal	 
                                                                        
      Write ( *, 300) Nspecies 
      Write (16, 300) Nspecies 
  300 FORMAT     ('--- New number of species=',i3,/                     &
               3x,'Specie',2x,                                         &
               'Nparticles Ntotal  Weight/part Weight_Total ')         
  310 FORMAT     (T5,i3,T13,i8,T22,i8,T30,2g10.3) 
      Wtotal = 0. 
      Do i = 1, Nspecies 
          If (i.eq.1) Then 
             Wtotal = wspecies (i) * lspecies (i) 
             Write ( *, 310) i, lspecies (i), lspecies (i), wspecies (i),      &
                                       Wtotal                                                            
            Write (16, 310) i, lspecies (i), lspecies (i), wspecies (i),      &
                                      Wtotal                                                            
         Else 
            ww = wspecies (i) * (lspecies (i) - lspecies (i - 1) ) 
           Write ( *, 310) i, lspecies (i) - lspecies (i - 1), lspecies (i), &
                                      wspecies (i), ww                                                  
            Write (16, 310) i, lspecies (i) - lspecies (i - 1), lspecies (i), &
                                     wspecies (i), ww                                                  
            Wtotal = Wtotal + ww 
      Endif 
      Enddo 
      Write (*,*) ' Total weight =', Wtotal 
                                                                        
      Return 
      END  Subroutine Species_Change                                         
                                                                        
!---------------------------------------------------                    
!                                  Read  current data from disk/tape,   
!                                  Open files                           
!                                  Nrecl is the number of values in a re
!                                  Npage is the number of particles in a
      SUBROUTINE FilesOpen 

Character (LEN=50) :: FileName
Logical            :: FileExists 

      write(*,'(10x,a/,2(30x,a/))')'Create temporary Files====>',NamesDir



!                                     Open file on a tape               
      OPEN (UNIT = 9, FILE = 'PMcrd.DAT', form = 'unformatted') 
                            ! this clears old header file               
      Write (9) 
      CLOSE (9) 
      OPEN (9, FILE = 'PMcrd.DAT', FORM = 'UNFORMATTED', STATUS =       &
         'UNKNOWN')                                                        
      OPEN (UNIT = 16, FILE = 'RESNEW.DAT') 
      OPEN (60, file = 'pt.dat', form = 'unformatted') 
                              

      RETURN 
      END SUBROUTINE FilesOpen                      
!---------------------------------------------                          
!                       Write current data to  disk/tape                
!                                                                       
      SUBROUTINE FilesWrite 
!----------------------------------------------                         
!                                       write header and control data   
REAL               ::    t0=0,t1=0.        ! counters for timing

      CALL CPU_TIME(t0)                    ! initialize time counter

      Write (9) HEADER, AEXPN, AEXP0, AMPLT, ASTEP, ISTEP, PARTW, TINTG,&
      EKIN, EKIN1, EKIN2, AU0, AEU0, NROWC, NGRIDC, Nspecies, Nseed,    &
      Om0, Oml0, hubble, Wp5, Ocurv, extras                             
      REWIND 9 

      NBYTE = NRECL*4
      NACCES= NBYTE / nbyteword

      OPEN(UNIT=21,FILE='PMcrs0.DAT',ACCESS='DIRECT', &
     	               FORM='unformatted',STATUS='UNKNOWN',RECL=NACCES)
      OPEN(31,FILE=NamesDir(1),form='unformatted')  ! x coord and veloc
      OPEN(32,FILE=NamesDir(2),form='unformatted')  ! y
      XMAX = FLOAT (NGRID) + 1. 
      XSHF = FLOAT (NGRID) 
      Ibuff =0
      KROW =0


   Npages  = (lspecies (Nspecies)-1)/NPAGE+1
   If(Npages.eq.0)stop ' Zero number of pages requested. Stop'
Do i=1,Npages
   If(i==Npages)Then
      NinPage = lspecies (Nspecies) -(i-1)*NPAGE  ! # particles in the current page
   Else
      NinPage = NPAGE
   EndIf
   jfirst  = (i-1)*NPAGE +1
   jlast   = jfirst + NinPage-1
   write(*,'(10x,a,i5,a,4i9)')'Write page=',i,' Particles=',NinPage,jfirst,jlast
   read(31)XPAR
   read(31)VX
   read(32)YPAR
   read(32)VY
   Do j = jfirst,jlast
      ZPAR(j-jfirst+1) =XPt(j)
      VZ (j-jfirst+1) =Vxt(j)
   EndDo
   Do j=1,NinPage
!                       Impose  Periodical boundary conditions         
        IF (XPAR(j) .GT.XMAX) XPAR(j) = XPAR(j) - XSHF 
        IF (XPAR(j) .LE.1.)   XPAR(j) = XPAR(j) + XSHF 
        If (XPAR (j) .lt.1.)write(*,*) ' Boundary Error X!!', i,j,XPAR(j)
        If (XPAR (j) .ge.XMAX)Then
            write (*,*)' Fixing Boundary:', j, XPAR(j)  
            XPAR (j) = XPAR (j) - 1.e-4 
        ENDIF
        IF (YPAR(j) .GT.XMAX) YPAR(j) = YPAR(j) - XSHF 
        IF (YPAR(j) .LE.1.)   YPAR(j) = YPAR(j) + XSHF 
        If (YPAR (j) .lt.1.)write(*,*) ' Boundary Error Y!!', i,j,YPAR(j)
        If (YPAR (j) .ge.XMAX)Then
            write (*,*)' Fixing Boundary:', j, YPAR(j)  
            YPAR (j) = YPAR (j) - 1.e-4 
        ENDIF
        IF (ZPAR(j) .GT.XMAX) ZPAR(j) = ZPAR(j) - XSHF 
        IF (ZPAR(j) .LE.1.)   ZPAR(j) = ZPAR(j) + XSHF 
        If (ZPAR (j) .lt.1.)write(*,*) ' Boundary Error Z!!', i,j,ZPAR(j)
        If (ZPAR (j) .ge.XMAX)Then
            write (*,*)' Fixing Boundary:', j, ZPAR(j)  
            ZPAR (j) = ZPAR (j) - 1.e-4 
        ENDIF
   EndDo                                    
   WRITE (21,REC=i) RECDAT
EndDo                            ! end lspecies loop                             
rewind(31)                       ! erase temporary files
write(31)NROW
rewind(32)
write(32)NROW
      CALL CPU_TIME(t1)
      timing(7) = timing(7) +t1-t0        ! time for VECTOR

      RETURN 
      END SUBROUTINE FilesWrite                     
!--------------------------------------                                 
!		normal random numbers                                                
      FUNCTION GAUSS3 (gSet,iFlag) 
!                        Uses ranlux for homogenous rand numbers        
!                        Uses Box-Muller inversion for gaussian         
!                 Excellent quality and speed                           
!  N=       100000000  GAUSS3+ranlux ------                             
!  sigma= 0.99997     mean= 0.18692E-03                                 
!  sigma Frac(>sig) Frac(<-sig)   True       n(>)    n(<)               
!  1.00 0.1587     0.1586     0.1587     15869095 15857724              
!  1.50 0.6682E-01 0.6680E-01 0.6681E-01  6681652  6679768              
!  2.00 0.2277E-01 0.2273E-01 0.2275E-01  2276651  2273197              
!  2.50 0.6222E-02 0.6206E-02 0.6210E-02   622190   620610              
!  3.00 0.1356E-02 0.1347E-02 0.1350E-02   135628   134674              
!  4.00 0.3242E-04 0.3145E-04 0.3170E-04     3242     3145              
!                                                                       
!--------------------------------------                                 
      DIMENSION RanNum (2)
      If (iFlag.eq.0) Then 
    1    CALL ranlux (RanNum, 2 ) 
         x1 = 2. * RanNum (1) - 1. 
         x2 = 2. * RanNum (2) - 1. 
         R = x1**2 + x2**2 
         If (R.ge.1.) GoTo 1 
         F              = sqrt ( - 2. * LOG (R) / R) 
         gSet        = x1 * F 
         GAUSS3 = x2 * F 
         iFlag        = 1 
      Else 
         GAUSS3 = gSet 
         iFlag        = 0 
      EndIf 
      RETURN 
      END FUNCTION GAUSS3                           
                                                                        
                                                                        
!--------------------------------------------------------------------   
                                                                        
!   LUXURY LEVELS.                                                      
!  level 0  (p=24): equivalent to the original RCARRY of Marsaglia      
!           and Zaman, very long period, but fails many tests.          
!  level 1  (p=48): considerable improvement in quality over level 0,   
!           now passes the gap test, but still fails spectral test.     
!  level 2  (p=97): passes all known tests, but theoretically still     
!           defective.                                                  
!  level 3  (p=223): DEFAULT VALUE.  Any theoretically possible         
!           correlations have very small chance of being observed.      
!  level 4  (p=389): highest possible luxury, all 24 bits chaotic.      
!---------------------------------------------                          
      SUBROUTINE ranlux (rvec, lenv) 
!                    returns a vector RVEC of LEN                       
!                   32-bit random floating point numbers between        
!                   zero (not included) and one (also not incl.).       
!           Next call to ranlux gives next LEN of random numbers        
!---------------------------------------------                          
      INTEGER lenv 
      REAL rvec (lenv)
      INTEGER iseeds (24) 

INCLUDE 'luxuryp.h' 

      DATA ndskip / 0, 24, 73, 199, 365 / 
      DATA i24, j24, luxlev / 24, 10, lxdflt / 
      DATA notyet / .true. / 
      DATA in24, kount, mkount, carry / 0, 0, 0, 0. / 
                                                                        
                                                                        
!  NOTYET is .TRUE. if no initialization has been performed yet.        
!              Default Initialization by Multiplicative Congruential    
                                                                        
      IF (notyet) THEN 
         stop ' Access to runlux before initialization STOP'
      ENDIF 
                                                                        
!          The Generator proper: "Subtract-with-borrow",                
!          as proposed by Marsaglia and Zaman,                          
!          Florida State University, March, 1989                        
                                                                        
      DO ivec = 1, lenv 
      uni = seeds (j24) - seeds (i24) - carry 
      IF (uni.LT.0.) THEN 
         uni = uni + 1.0 
         carry = twom24 
      ELSE 
         carry = 0. 
      ENDIF 
      seeds (i24) = uni 
      i24 = next (i24) 
      j24 = next (j24) 
      rvec (ivec) = uni 
      !  small numbers (with less than 12 "significant" bits) are "padde
      IF (uni.LT.twom12) THEN 
         rvec (ivec) = rvec (ivec) + twom24 * seeds (j24) 
      !        and zero is forbidden in case someone takes a logarithm  
         IF (rvec (ivec) .EQ.0.) rvec (ivec) = twom24 * twom24 
      ENDIF 
      !        Skipping to luxury.  As proposed by Martin Luscher.      
      in24 = in24 + 1 
      IF (in24.EQ.24) THEN 
         in24 = 0 
         kount = kount + nskip 
         DO isk = 1, nskip 
         uni = seeds (j24) - seeds (i24) - carry 
         IF (uni.LT.0.) THEN 
            uni = uni + 1.0 
            carry = twom24 
         ELSE 
            carry = 0. 
         ENDIF 
         seeds (i24) = uni 
         i24 = next (i24) 
         j24 = next (j24) 
         ENDDO 
      ENDIF 
      ENDDO 
      kount = kount + lenv 
      IF (kount.GE.igiga) THEN 
         mkount = mkount + 1 
         kount = kount - igiga 
      ENDIF 
      RETURN 
                                                                        
          ! SUBROUTINE ranlux                                           
      END SUBROUTINE ranlux                         
                                                                        
!---------------------------------------------                          
!                    Subroutine to initialize from one or three integers
      SUBROUTINE rluxgo (lux, ins, k1, k2) 
!                initializes the generator from                         
!               one 32-bit integer INT and sets Luxury Level LUX        
!               which is integer between zero and MAXLEV, or if         
!               LUX .GT. 24, it sets p=LUX directly.  K1 and K2         
!               should be set to zero unless restarting at a break      
!               point given by output of RLUXAT                         
!---------------------------------------------                          
INCLUDE'luxuryp.h' 
      INTEGER iseeds (24)
                                                                        
      IF (lux.LT.0) THEN 
         luxlev = lxdflt 
      ELSEIF (lux.LE.maxlev) THEN 
         luxlev = lux 
      ELSEIF (lux.LT.24.OR.lux.GT.2000) THEN 
         luxlev = maxlev 
         Write (6, '(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ', lux 
      ELSE 
         luxlev = lux 
         DO ilx = 0, maxlev 
         IF (lux.EQ.ndskip (ilx) + 24) luxlev = ilx 
         ENDDO 
      ENDIF 
      IF (luxlev.LE.maxlev) THEN 
         nskip = ndskip (luxlev) 
      Write (6, '(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :', &
                       luxlev, '     P=', nskip + 24                                      
      ELSE 
         nskip = luxlev - 24 
         Write (6, '(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',       &
         luxlev                                                         
      ENDIF 
      in24 = 0 
      IF (ins.LT.0) Write (6,'(A)')' Illegal initialization by RLUXGO,'&
                               'negative input seed'                                             
      IF (ins.GT.0) THEN 
         jseed = ins 
      Write (6, '(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS', &
                                        jseed, k1, k2                                                     
      ELSE 
         jseed = jsdflt 
      Write (6, '(A)') ' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED' 
      ENDIF 
      inseed = jseed 
      notyet = .false. 
      twom24 = 1. 
      DO i = 1, 24 
      twom24 = twom24 * 0.5 
      k = jseed / 53668 
      jseed = 40014 * (jseed-k * 53668) - k * 12211 
      IF (jseed.LT.0) jseed = jseed+icons 
      iseeds (i) = MOD (jseed, itwo24) 
      ENDDO 
      twom12 = twom24 * 4096. 
      DO i = 1, 24 
      seeds (i) = REAL (iseeds (i) ) * twom24 
      next (i) = i - 1 
      ENDDO 
      next (1) = 24 
      i24 = 24 
      j24 = 10 
      carry = 0. 
      IF (seeds (24) .EQ.0.) carry = twom24 
      !        If restarting at a break point, skip K1 + IGIGA*K2       
      !        Note that this is the number of numbers delivered to     
      !        the user PLUS the number skipped (if luxury .GT. 0).     
      kount = k1 
      mkount = k2 
      IF (k1 + k2.NE.0) THEN 
         DO iouter = 1, k2 + 1 
         inner = igiga 
         IF (iouter.EQ.k2 + 1) inner = k1 
         DO isk = 1, inner 
         uni = seeds (j24) - seeds (i24) - carry 
         IF (uni.LT.0.) THEN 
            uni = uni + 1.0 
            carry = twom24 
         ELSE 
            carry = 0. 
         ENDIF 
         seeds (i24) = uni 
         i24 = next (i24) 
         j24 = next (j24) 
         ENDDO 
         ENDDO 
      !         Get the right value of IN24 by direct calculation       
         in24 = MOD (kount, nskip + 24) 
         IF (mkount.GT.0) THEN 
            izip = MOD (igiga, nskip + 24) 
            izip2 = mkount * izip + in24 
            in24 = MOD (izip2, nskip + 24) 
         ENDIF 
      !       Now IN24 had better be between zero and 23 inclusive      
         IF (in24.GT.23) THEN 
      Write (6, '(A/A,3I11,A,I5)') '  Error in RESTARTING with RLUXGO:',&
           '  The values', ins, k1, k2, ' cannot occur at luxury level', luxlev 
            in24 = 0 
         ENDIF 
      ENDIF 
      RETURN 
                                                                        
          ! SUBROUTINE rluxgo                                           
      END SUBROUTINE rluxgo                         
                                                                        
!---------------------------------- Read in variables                   
      REAL FUNCTION INPUT (text) 
!------------------------------------------------                       
      Character text * (*) 
         Write (*,'(A,$)') text 
         READ (*,*) x 
         INPUT = x 
      Return 
      END FUNCTION INPUT                            
!--------------------------------------   Fourier Transform             
      SUBROUTINE SETF67 (JBC1, JQ1, jbc, jp, jsl, ll1, k2, k3, k4, k7,  &
      Qi, jndx, Uf, Vf)                                                 
!     -----------------------------------------------------             
                                                                        
      DIMENSION Uf (marr), Vf (marr) 
      DIMENSION Qi (mf67), jndx (mf67) 
                                                                        
      PI = DATAN (1.D0) * 4. 
      JBC = JBC1 
      JQ = JQ1 
      IF (JBC.LT.3) GOTO 101 
      JQ = JQ - 1 
  101 K3 = 2**JQ 
      K7 = K3 / 2 
      N5 = K3 / 4 
      I = 1 
      JNDX (I) = N5 
      QI (I) = 0.5 * SQRT (2.) 
      K = I 
      I = I + 1 
  102 IL = I 
      IF (I.EQ.K7) GOTO 104 
  103 K1 = JNDX (K) / 2 
      JNDX (I) = K1 
      QI (I) = SIN (PI * FLOAT (K1) / FLOAT (K3) ) 
      K1 = K7 - K1 
      I = I + 1 
      JNDX (I) = K1 
      QI (I) = SIN (PI * FLOAT (K1) / FLOAT (K3) ) 
      K = K + 1 
      I = I + 1 
      IF (K.EQ.IL) GOTO 102 
      GOTO 103 
  104 RETURN 
      END SUBROUTINE SETF67                         
!-----------------------------------------                              
      SUBROUTINE TFOLD (IS, L, ZZZ, k2) 
                                                                        
      DIMENSION ZZZ (MARR) 
      IH2 = K2 / 2 - 1 
      DO 100 I = IS, IH2 
         I1 = I + L 
         I2 = K2 - I + L 
         A = ZZZ (I1) 
         ZZZ (I1) = A - ZZZ (I2) 
         ZZZ (I2) = A + ZZZ (I2) 
  100 END DO 
      RETURN 
      END SUBROUTINE TFOLD                          
!------------------------------------                                   
      SUBROUTINE NEG (I1, I3, I2, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi,&
      jndx, Uf, Vf)                                                     

      DIMENSION Uf (marr), Vf (marr) 
      DIMENSION Qi (mf67), jndx (mf67) 
                                                                        
      DO  K = I1, I3, I2 
         Vf (K) = - Vf (K) 
      END DO 
      RETURN 
      END SUBROUTINE NEG                            
!-------------------------------------                                  
      SUBROUTINE REVNEG (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx,   &
      Uf, Vf)                                                           

      DIMENSION Uf (marr), Vf (marr) 
      DIMENSION Qi (mf67), jndx (mf67) 
                                                                        
      INTEGER jbc, jp, jsl, ll1, k2, k3, k4, k7 
                                                                        
      DO  I = 1, K7 
         J = K3 + 1 + I 
         K = K4 + 1 - I 
         A = Vf (J) 
         Vf (J) = - Vf (K) 
         Vf (K) = - A 
      END DO 
      RETURN 
      END SUBROUTINE REVNEG                         
!-------------------------------------                                  
      SUBROUTINE ZEERO (L, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, &
      Uf, Vf)                                                           

      DIMENSION Uf (marr), Vf (marr) 
      DIMENSION Qi (mf67), jndx (mf67) 
                                                                        
      INTEGER jbc, jp, jsl, ll1, k2, k3, k4, k7 
                                                                        
      DO I = 1, K2 
         LI = L + I 
         Uf (LI - 1) = 0.0 
      END DO 
      RETURN 
      END SUBROUTINE ZEERO                          
!------------------------------------------                             
      SUBROUTINE TFOLD1 (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx,   &
      Uf, Vf)                                                           

      DIMENSION Uf (marr), Vf (marr) 
      DIMENSION Qi (mf67), jndx (mf67) 
                                                                        
      INTEGER jbc, jp, jsl, ll1, k2, k3, k4, k7 
                                                                        
      II = K2 - 1 
      DO  I = 1, II 
         I1 = JSL + I 
         I2 = LL1 - I 
         A = Uf (I1) 
         Uf (I1) = A + Uf (I2) 
         Uf (I2) = A - Uf (I2) 
      END DO 
      RETURN 
      END SUBROUTINE TFOLD1                         
!-------------------------------------                                  
      SUBROUTINE KFOLD (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf,&
      Vf)                                                               

      DIMENSION Uf (marr), Vf (marr) 
      DIMENSION Qi (mf67), jndx (mf67) 
                                                                        
      INTEGER jbc, jp, jsl, ll1, k2, k3, k4, k7 
                                                                        
      JS1 = K2 
      I = 1 
      J5 = JSL + K2 
      IS1 = JSL 
      IC1 = LL1 
      JS1 = JS1 / 2 
      IF (JS1.NE.1) GOTO 200 
      K1 = JNDX (I) 
      SN = QI (I) 
      IS0 = IS1 
      IS1 = IS1 + JS1 
      IC0 = IC1 
      IC1 = IC1 + JS1 
      ODD1 = SN * (Uf (IC1) - Uf (IS1) ) 
      Vf (K1 + 1) = Uf (IC0) + ODD1 
      K3MK1 = K3 - K1 
      Vf (K3MK1 + 1) = Uf (IC0) - ODD1 
      IF (JBC.LT.3) GOTO 110 
      ODD2 = SN * (Uf (IC1) + Uf (IS1) ) 
      K3PK1 = K3 + K1 
      Vf (K3PK1 + 1) = Uf (IS0) + ODD2 
      K4MK1 = K4 - K1 
      Vf (K4MK1 + 1) = - Uf (IS0) + ODD2 
  110 RETURN 
  200 SN = QI (I) 
      IS1 = IS1 + JS1 
      IC1 = IC1 + JS1 
      J3 = IS1 + JS1 
  210 IS0 = IS1 - JS1 
      IC0 = IC1 - JS1 
      ODD1 = SN * (Uf (IC1) - Uf (IS1) ) 
      ODD2 = SN * (Uf (IC1) + Uf (IS1) ) 
      Uf (IC1) = Uf (IC0) - ODD1 
      Uf (IS1) = - Uf (IS0) + ODD2 
      Uf (IC0) = Uf (IC0) + ODD1 
      Uf (IS0) = Uf (IS0) + ODD2 
      IS1 = IS1 + 1 
      IC1 = IC1 + 1 
      IF (IS1.NE.J3) GOTO 210 
      I = I + 1 
  300 IS1 = JSL 
      IC1 = LL1 
      JS1 = JS1 / 2 
      IF (JS1.EQ.1) GOTO 400 
  310 SN = QI (I) 
      I = I + 1 
      CS = QI (I) 
      IS1 = IS1 + JS1 
      IC1 = IC1 + JS1 
      J3 = IS1 + JS1 
  320 IS0 = IS1 - JS1 
      IC0 = IC1 - JS1 
      ODD1 = CS * Uf (IC1) - SN * Uf (IS1) 
      ODD2 = SN * Uf (IC1) + CS * Uf (IS1) 
      Uf (IC1) = Uf (IC0) - ODD1 
      Uf (IC0) = Uf (IC0) + ODD1 
      Uf (IS1) = - Uf (IS0) + ODD2 
      Uf (IS0) = Uf (IS0) + ODD2 
      IS1 = IS1 + 1 
      IC1 = IC1 + 1 
      IF (IS1.NE.J3) GOTO 320 
      IS1 = IS1 + JS1 
      IC1 = IC1 + JS1 
      J3 = IS1 + JS1 
  330 IS0 = IS1 - JS1 
      IC0 = IC1 - JS1 
      ODD1 = SN * Uf (IC1) - CS * Uf (IS1) 
      ODD2 = CS * Uf (IC1) + SN * Uf (IS1) 
      Uf (IC1) = Uf (IC0) - ODD1 
      Uf (IC0) = Uf (IC0) + ODD1 
      Uf (IS1) = - Uf (IS0) + ODD2 
      Uf (IS0) = Uf (IS0) + ODD2 
      IS1 = IS1 + 1 
      IC1 = IC1 + 1 
      IF (IS1.NE.J3) GOTO 330 
      I = I + 1 
      IF (IS1.EQ.J5) GOTO 300 
      GOTO 310 
  400 K1 = JNDX (I) 
      SN = QI (I) 
      I = I + 1 
      CS = QI (I) 
      IS0 = IS1 
      IS1 = IS1 + JS1 
      IC0 = IC1 
      IC1 = IC1 + JS1 
      ODD1 = CS * Uf (IC1) - SN * Uf (IS1) 
      Vf (K1 + 1) = Uf (IC0) + ODD1 
      K3MK1 = K3 - K1 
      Vf (K3MK1 + 1) = Uf (IC0) - ODD1 
      IF (JBC.LT.3) GOTO 410 
      ODD2 = SN * Uf (IC1) + CS * Uf (IS1) 
      K3PK1 = K3 + K1 
      Vf (K3PK1 + 1) = Uf (IS0) + ODD2 
      K4MK1 = K4 - K1 
      Vf (K4MK1 + 1) = - Uf (IS0) + ODD2 
  410 IS1 = IS1 + 1 
      IC1 = IC1 + 1 
      K1 = JNDX (I) 
      IS0 = IS1 
      IS1 = IS1 + JS1 
      IC0 = IC1 
      IC1 = IC1 + JS1 
      ODD1 = SN * Uf (IC1) - CS * Uf (IS1) 
      Vf (K1 + 1) = Uf (IC0) + ODD1 
      K3MK1 = K3 - K1 
      Vf (K3MK1 + 1) = Uf (IC0) - ODD1 
      IF (JBC.LT.3) GOTO 420 
      ODD2 = CS * Uf (IC1) + SN * Uf (IS1) 
      K3PK1 = K3 + K1 
      Vf (K3PK1 + 1) = Uf (IS0) + ODD2 
      K4MK1 = K4 - K1 
      Vf (K4MK1 + 1) = - Uf (IS0) + ODD2 
  420 IS1 = IS1 + 1 
      IC1 = IC1 + 1 
      I = I + 1 
      IF (IS1.NE.J5) GOTO 400 
      RETURN 
      END SUBROUTINE KFOLD                          
                                                                        
!      ------------------------------------------------------           
      SUBROUTINE FOUR67 (JBC1, JQ1, jp1, jsl1, ll11, k21, k71, Qi, jndx,&
      Uf, Vf)                                                           
!      ------------------------------------------------------           

      DIMENSION Uf (marr), Vf (marr) 
      DIMENSION Qi (mf67), jndx (mf67) 
                                                                        
      INTEGER jbc, jp, jsl, ll1, k2, k3, k4, k7 
                                                                        
      JBC = JBC1 
      JQ = JQ1 
      jp = jp1 
      jsl = jsl1 
      ll1 = ll11 
      k2 = k21 
      k7 = k71 
      A5 = 0.5 * SQRT (2.0) 
      K4 = 2**JQ 
      K3 = K4 
      GOTO (103, 103, 101, 102), JBC 
  101 Uf (1) = Uf (1) / 2.0 
      Uf (K3 + 1) = Uf (1) 
      K2 = K3 
                                                                        
      CALL TFOLD (0, 1, Uf, k2) 
                                                                        
      K3 = K3 / 2 
      JQ = JQ - 1 
      GOTO 103 
  102 K3 = K3 / 2 
      JQ = JQ - 1 
  103 N5 = K3 / 4 
      K7 = K3 / 2 
      N11 = 3 * K7 
      K31 = K3 + 1 
      GOTO (300, 400, 500, 600), JBC 
  300 Uf (K31) = 0.0 
      Uf (1) = 0.0 
      K2 = K3 
      DO 301 I = 2, JQ 
                                                                        
         CALL TFOLD (1, 1, Uf, k2) 
                                                                        
  301 K2 = K2 / 2 
      Vf (K7 + 1) = Uf (2) 
      JF = N5 
      JSL = 1 
      DO 302 JP = 2, JQ 
         LL1 = K2 + 1 
                                                                        
         CALL ZEERO (1, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf,&
         Vf)                                                            
         CALL KFOLD (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf,   &
         Vf)                                                            
                                                                        
         I1 = 3 * JF + 1 
         I2 = 4 * JF 
         I3 = I1 + (K2 / 2 - 1) * I2 
                                                                        
         CALL NEG (I1, I3, I2, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi,   &
         jndx, Uf, Vf)                                                  
                                                                        
         K2 = K2 + K2 
  302 JF = JF / 2 
      RETURN 
  400 Uf (1) = Uf (1) / 2.0 
      Uf (K31) = Uf (K31) / 2.0 
      K2 = K3 
      DO 401 I = 2, JQ 
                                                                        
         CALL TFOLD (0, K31 - K2, Uf, k2) 
                                                                        
  401 K2 = K2 / 2 
      LL1 = K31 - K2 
      A = Uf (LL1) + Uf (LL1 + 2) 
      Vf (1) = - A - Uf (LL1 + 1) 
      Vf (K31) = - A + Uf (LL1 + 1) 
      Vf (K7 + 1) = Uf (LL1) - Uf (LL1 + 2) 
      DO 402 JP = 2, JQ 
         JSL = K31 - K2 
         LL1 = JSL - K2 
         idum = jsl 
                                                                        
         CALL ZEERO (idum, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, &
         Uf, Vf)                                                        
         CALL KFOLD (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf,   &
         Vf)                                                            
                                                                        
  402 K2 = K2 + K2 
                                                                        
      CALL NEG (1, K31, 2, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, &
      Uf, Vf)                                                           
                                                                        
      RETURN 
  500 K2 = K3 
      L2 = K4 
      DO 501 JP = 2, JQ 
                                                                        
         CALL TFOLD (1, 1, Uf, k2) 
         CALL TFOLD (0, L2 - K2 + 1, Uf, k2) 
                                                                        
  501 K2 = K2 / 2 
      LL1 = L2 - K2 + 1 
      A = Uf (LL1) + Uf (LL1 + 2) 
      Vf (K7 + 1) = 2.0 * ( - Uf (LL1) + Uf (LL1 + 2) ) 
      Vf (1) = 2.0 * (A + Uf (LL1 + 1) ) 
      Vf (K31) = 2.0 * (A - Uf (LL1 + 1) ) 
      Vf (N11 + 1) = 2.0 * Uf (2) 
      DO 502 JP = 2, JQ 
         Uf (K2 + 1) = 2.0 * Uf (K2 + 1) 
         JSL = K2 + 1 
                                                                        
         CALL TFOLD1 (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf,  &
         Vf)                                                            
                                                                        
         LL1 = LL1 - K2 
         Uf (LL1) = - 2.0 * Uf (LL1) 
                                                                        
         CALL KFOLD (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf,   &
         Vf)                                                            
                                                                        
  502 K2 = K2 + K2 
      Vf (1) = Vf (1) * A5 
      Vf (K31) = Vf (K31) * A5 
      RETURN 
  600 Uf (1) = Uf (1) * A5 
      Uf (K31) = Uf (K31) * A5 
      K2 = K3 
      DO 601 JP = 2, JQ 
                                                                        
         CALL TFOLD (0, K31 - K2, Uf, k2) 
         CALL TFOLD (1, K31, Uf, k2) 
                                                                        
  601 K2 = K2 / 2 
      LL1 = K31 - K2 
      A = Uf (LL1) + Uf (LL1 + 2) 
      Vf (1) = 2.0 * (A + Uf (LL1 + 1) ) 
      Vf (K31) = 2.0 * (A - Uf (LL1 + 1) ) 
      Vf (K7 + 1) = 2.0 * ( - Uf (LL1) + Uf (LL1 + 2) ) 
      Vf (N11 + 1) = 2.0 * Uf (K3 + 2) 
      DO 602 JP = 2, JQ 
         JSL = K31 + K2 
         Uf (JSL) = 2.0 * Uf (JSL) 
                                                                        
         CALL TFOLD1 (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf,  &
         Vf)                                                            
                                                                        
         LL1 = LL1 - K2 
         Uf (LL1) = - 2.0 * Uf (LL1) 
                                                                        
         CALL KFOLD (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf,   &
         Vf)                                                            
                                                                        
  602 K2 = K2 + K2 
                                                                        
      CALL REVNEG (jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx, Uf, Vf) 
                                                                        
      idum = k3 
                                                                        
      CALL NEG (2, idum, 2, jbc, jp, jsl, ll1, k2, k3, k4, k7, Qi, jndx,&
      Uf, Vf)                                                           
                                                                        
      K2 = K4 
                                                                        
      CALL TFOLD (1, 1, Vf, k2) 
                                                                        
      Vf (K4 + 1) = 0.0 
                                                                        
      RETURN 
      END SUBROUTINE FOUR67                         
end module setInitialConditions
!
!
!
!-------------------------------------------------------------------------
PROGRAM  PMstartMp

use      setInitialConditions

!		   NROW = number of particles in 1 dimension                         
!		   NGRID= number of cells        in 1 dimension                      
!		   Npage =  number of particles per page = NROW**2               
!                                                                       
COMMON / TRUNCOM / Om, Omb, Omc, Omnu, Par (6), ns, qqscaleb,QSCALE
DOUBLE PRECISION ::    Wtotal=0., SKINE=0.
REAL, PARAMETER  ::    PI=3.1415926535
REAL             ::    t0=0,t1=0.        ! counters for timing
Logical          ::    iCont =.true.     ! continue (true) after setting seeds

      CALL CPU_TIME(t0)
      Open(30,file='TestDataParallel.dat')
      Wtotal = (float(NROW)**3*4+2.*Nmaxpart*4. &
                + 31.*NROW*4 +float(Nblocks)**3*4 &
                + 6. *NPAGE*4)/1024.**3
      write (*,'(//10x,a,g13.4,a/)') ' Memory required=',Wtotal,'Gb'
      CALL InitValues (NBYTE, SCALEL) 
      Write (*,*) ' InitValues is done. N bytes=', NBYTE 
      Max_lev_abs = INT (ALOG (REAL (Lblock) ) / ALOG (2.) + 1.1) 
      Write (*,'(10x,a,2i4)') 'Absolute Limits on mass resolution levels: ', 1,&
                            Max_lev_abs
      Write (*,'(10x,2(a,i4))') 'Actual Max resolution Level=', Nmax_lev, & 
                              ' Min Level=', Nmin_lev
      Write (16,'(10x,2(a,i4))') 'Actual Max resolution Level=', Nmax_lev, &
                              ' Min Level=', Nmin_lev
      If (Nmax_lev.gt.Max_lev_abs) Then 
               Write (*,*) ' Wrong number of Levels for mass refinement (',  &
                                  Nmax_lev, ') Block size =', Lblock, ' Use Absolute value=', &
                Max_lev_abs                                                       
                Nmax_lev = Max_lev_abs 
      Endif 
                                                                        
      CALL FilesOpen                   ! this opens files on disk        
         Write (16, * ) ' Seed=', Nseed 
         AEXP0 = AEXPN 
         AEXPV = AEXPN - ASTEP / 2. 

      CALL ZeroMASK                ! set mask for multiple mass resolution
      If (Nmax_lev.ne.Nmin_lev) Then 
         Open (30, file = 'SelCoords.DAT', status = 'unknown') 
         CALL SetMASK 
      EndIf 
                                     ! set scaling constants for coordinates and velocities
         IFOUR = INT (ALOG (FLOAT (NROW) ) / ALOG (2.) + 0.5) 
            Write (*,'(3x,a12,i4,a,i4)') ' Ifour =', IFOUR, ' Nrow=', NROW 
         Fact   = sqrt (Om0 + Oml0 * AEXPV**3 + Ocurv * AEXPV) 
         QFACT  = FLOAT (NGRID) / FLOAT (NROW) 
         Vscale = ScaleL * 100. * hubble / NGRID 

			      !  get a realization of spectrum in /GRID/          
         NRAND = Nseed 
      CALL CPU_TIME(t1)
      timing(1) = t1-t0       ! time for intialization

      CALL SetRandom(iCont)            ! set seeds for ramdom numbers

         If(.NOT.iCont)Then
            write(*,'(2(10x,a,f10.2/))') &
                     'cpu_time for init =',timing(1), &
                     'cpu_time for seeds=',timing(2)
            Stop
         EndIf
         SKINE   = 0.                              ! counter for  kinetic energy
         Wtotal  = 0.                             ! counter for  total weight

      ALLOCATE(GR(NROW,NROW,NROW))     ! allocate memory for spectrum

      Do iDirection = 1, 3                  ! ----------------- loop over x,y,z directions 
          Write(*,'(//3x,a12,i3,15("-"))' )'- Direction=',iDirection
          If(iDirection < 3)  &          
                 OPEN(1,FILE=NamesDir(iDirection), FORM='UNFORMATTED',        &
                        STATUS ='UNKNOWN')                                  

          CALL SPECTR (ALPHA, NRAND, iDirection) 
			       !   get the displacement vector by FFT 
          VCONS = - ALPHA/(2.*PI/NGRID)*(AEXPV/AEXP0)*SQRT(AEXPV)*Fact
          XCONS = 	ALPHA / (2. * PI / NGRID) * (AEXPN / AEXP0) 
            Write(*,'(3x,a12,g12.4,a,g12.4)' ) 'Scaling:(x)=', XCONS, ' (v)=', VCONS 

          CALL VECTOR (IFOUR) 
             Write (*,'(3x,a12,i4,a,g12.4)' ) 'SPECTR done:', IFOUR, ' Alpha=', ALPHA 
            write(*,'(2(a,f8.3))')' Spectr t=',timing(3),' FFT=',timing(5)
			       !    find x,v and write to file, page by page
          Do i = 1, Max_lev_abs 
                lspecies (i) = 0 
          Enddo 
          Nspecies = 0 
          Icurrent = 0                            ! current particle   
          Do kSpec = Max_lev_abs, 1, - 1 
             CALL BLOCKS (XCONS, VCONS, kSpec, Icurrent, iDirection) 
          Enddo 
          CALL WriteData (Vscale, AexpV, Wtotal,SKINE,iDirection)
          CLOSE(1) 
      EndDo 
      
      DEALLOCATE (GR)                                   
      CALL Species_Change (Wtotal)     ! print info on species
      EKIN = 0.5 * SKINE / AEXPV**2 
         Write (*, '('' Ekin='',E12.4,'' Weght per cell='',g12.5,        &
                 '' (must be 1)'')') EKIN, Wtotal / NGRID**3            
         Write (16, '('' Ekin='',E12.4,'' Weght per cell='',g12.5,         &
     &         '' (must be 1)'')') EKIN, Wtotal / NGRID**3              

                       
      CALL FilesWrite           ! write header and control data
      Write (60) astep        ! write pt.dat file: time-step for particles             

      CALL CPU_TIME(t1)
            write(*,'(10x,a/,(10x,a15,f10.2))') &
                     'Wall cpu_time (seconds) ', &
                     'Init        =',timing(1), &
                     'Seeds       =',timing(2), &
                     'Spectr      =',timing(3), &
                     'Coords      =',timing(4), &
                     'FFT         =',timing(5), &
                     'Write temp  =',timing(6), &
                     'Write final =',timing(7), &
                     'Total       =',t1-t0

      STOP 
      END                                           
