C _______________________ START 3-D SIMULATIONS
C
C                         Klypin, August 1993
C
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMinitialCR.h'	  
C		   NROW = number of particles in 1 dimension
C		   NGRID= number of cells        in 1 dimension
C		   Npage =  number of particles per page = NROW**2
C


      COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Par(6),ns,qqscaleb,QSCALE
      COMMON / KINENRG/ SKINE,SX,SY,SZ,SX2,SY2,SZ2
      Character  Hd*5,Tail*4
      DATA PI		/3.1415926535/

      Hd  ='PMcrs'
      Tail='.DAT'
     
      CALL InitValues(NBYTE,SCALEL)
      write (*,*) ' InitValues is done. N bytes=',NBYTE
      NspecM   = Levels  ! Nspecies -1
      Nmin_lev = Levels  ! Nspecies
      write (*,*) ' Max resolution Level=',NspecM,
     &               ' Min Level=',Nmin_lev
      write (16,*) ' Max resolution Level=',NspecM,
     &               ' Min Level=',Nmin_lev
      If(Nspecies.ne.0.and.2**(Nspecies-1).ne.Lblock)Then
         write (*,*) ' Wrong number of Levels for mass refinement (',
     &                  Nspecies,') Block size =',Lblock,
     &                  '       2**(Levels-1) must = Block size'
      Endif 
      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')
      write (16,*) ' Seed=',Nseed
       AEXP0 = AEXPN
       AEXPV = AEXPN - ASTEP/2.
C                                    set mask for multiple mass resolution
       If(Nspecies.ne.0)Then
                CALL ZeroMASK(Nmin_lev)
c                If(Nspecies.ge.2)Then
c                   Open(30,file='SelCoords.DAT')
c                   CALL SetMASK(NspecM)
c                EndIf
      EndIf
C			                   get a realization of spectrum in /GRID/
       NRAND =Nseed
c       CALL SPectr_CR(NRAND)
c       If(Levels.gt.1)Then
           CALL SPECTR(NRAND,1)  ! set spectrum for higher levels
           CALL ZERO_SP         ! remove harmonics included at lower levels
c       Endif 
C			                        get the displacement vector by FFT
       IFOUR = INT(ALOG(FLOAT(NROW))/ALOG(2.)+0.5)
       Write (*,*) ' Ifour =',IFOUR,' Nrow=',NROW
       CALL VECTOR(IFOUR)
	    write (*,*) ' SPECTR is done. IFOUR = ',IFOUR,' Alpha=',
     &                   (ALPHA(i),i=1,Levels)
        Fact =  sqrt(Om0+Oml0*AEXPV**3+Ocurv*AEXPV)
        Do Lv =1,Levels
           VCONS(Lv) = -ALPHA(Lv)/(2.*PI/NGRID)*(AEXPV/AEXP0)*
     &                              SQRT(AEXPV)*Fact*2**(Lv-1)
           XCONS(Lv) =	 ALPHA(Lv)/(2.*PI/NGRID)*(AEXPN/AEXP0)
     &                               *2**(Lv-1)
        Enddo 
       QFACT =  FLOAT(NGRID)/FLOAT(NROW)
                      SKINE = 0.
                      SX    = 0.
                      SY    = 0.
                      SZ    = 0.
                      SX2   = 0.
                      SY2   = 0.
                      SZ2   = 0.

C			           find x,v and write to file, page by page
       Do i=1,Nspecies
          lspecies(i) =0
       Enddo 
       If(Nspecies.eq.0)Then     ! old code with one set                    
         DO  KROW = 1, NROW
            CALL SETXV(XCONS(1),VCONS(1),KROW)
c            CALL SETtest(XCONS,VCONS,KROW,NRAND)
	      CALL WRIROW(KROW,1)
         ENDDO
         lspecies(1) =NROW**3
         wspecies(1) =PARTW
         Wtotal = wspecies(1)*lspecies(1)
           EKIN = 0.5*SKINE/AEXPV**2*PARTW
           DP   = NROW**3
           SX   = SX/DP
           SY   = SY/DP
           SZ   = SZ/DP
           SX2  = SQRT(SX2/DP)
           SY2  = SQRT(SY2/DP)
           SZ2  = SQRT(SZ2/DP)
           WRITE (*,'('' Displacement of a particle: MEAN'',3F9.4)')
     .             SX,SY,SZ
           WRITE (*,'('' in cell units               RMS '',3F9.4)')
     .             SX2,SY2,SZ2
           WRITE (16,'('' Displacement of a particle: RMS'',3F9.4)')
     .             SX2,SY2,SZ2
           Veloc=sqrt( SKINE/AEXPV**2/NROW**3)*SCALEL/NGRID*
     .             100.*hubble
           WRITE (*,'('' Ekin='',E12.4,'' Veloc_cold(km/s)='',F8.2,
     .                     '' Weght per cell='',f8.2,'' (must be 1)'')')
     .            EKIN,Veloc,Wtotal/NGRID**3
           WRITE (16,'('' Ekin='',E12.4,'' Veloc_cold(km/s)='',F8.2,
     .                     '' Weght per cell='',f8.2,'' (must be 1)'')')
     .            EKIN,Veloc,Wtotal/NGRID**3
      Else
C                          find x,v for multiple mass resolution code 
         Icurrent    =0                         ! current particle
         Iprevious  =0
         Jspecies   = 0
            Do Lv =Levels,1,-1
               write (*,*) ' Go to BLOCKS. Level =',Lv
                CALL BLOCKS(Lv,Icurrent)   
                If(Icurrent.gt.Iprevious)Then
                   Jspecies             =Jspecies +1
                   lspecies(Jspecies)   =Icurrent
                   wspecies(Jspecies)   =PARTW *8/8**Lv
                   Iprevious            =Icurrent
                Endif 
            Enddo 
            Nspecies = Jspecies
           Vscale =ScaleL*100.*hubble/NGRID
      CALL WriteData(Vscale,AexpV,Wtotal)
      CALL Species_Print(Wtotal_o)      ! rearrange 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
      Endif 
C			        
      CALL WRTAPE   ! write header and control data
C                      write pt.dat file: time-step for particles
      If(Nspecies.eq.0)Then
         Nparticles =lspecies(1)
      Else
         Nparticles =lspecies(Nspecies)
      EndIf
      Do ic1 =1,Nparticles
         XPt(ic1) =astep
      Enddo 
      open ( 60 , file = 'pt.dat' , form = 'unformatted' )
      write(60) (XPt(ic1),ic1=1,Nparticles)
      write (*,*) ' pt=', (XPt(ic1),ic1=1,10)
      STOP
      END
C--------------------------------------------------
C                                                     sqrt(Power spectrum)
C                                                           k = (2pi/L) wk
      REAL  FUNCTION TRUNF(wk)
C--------------------------------------------------
      INCLUDE 'PMparameters.h'
      PARAMETER (NSPEC =NROW/2)
      COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Par(6),ns,qqscaleb,QSCALE
      Real                 ns,k
c      IF (WK.GE.FLOAT(NSPEC)) THEN
c	       TRUNF =0.
c	       RETURN
c      ENDIF
      If(Par(6).ne.0.)Then   ! Holtzman approx
        k = QSCALE*wk
        sk= sqrt(k)
        TRUNF = k**(ns/2.) /
     .         (1.+sk*(Par(2)
     .            +sk*(Par(3)
     .            +sk*(Par(4)
     .            +sk*(Par(5) )))) )**(Par(6))
      Else                   ! BBKS + Sugiyama approx
c        Gammaeff =Om*hsmall/exp(Omb*(1.+sqrt(hsmall/0.5)/Om))
c        Q = wk/hsmall/Gammaeff
        k = QSCALE*wk
        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
      RETURN
      END
C*********************************************************************
C	   Power spectrum for CDM. Wk-  wavenumber(Mpc^-1, h=0.5)
C                         P = wk *T^2(k)
C                         QSCALE = 2pi/Box_size,  QS =hubble**(-2)
c      FUNCTION TRUNF(WK)
c      COMMON / TRUNCOM/ QSCALE,g,gy
c      COMMON / TRUNCDM/ QS,A1,A2,A3,A4,A5,A32,A43,A54
c       Q = QSCALE*QS*WK
c       TRUNF = sqrt( WK*
c     +		 (LOG(1.+A1*Q)/(A1*Q))**2/
c     +		 SQRT(1.+Q*(A2 +Q*(A32   +Q*(A43+Q*A54)))) )
c
c      RETURN
c      END

C*********************************************************************
C			  INITIALIZE CONSTANTS:
C			      Scalel    =length of square in MPC
C			      AEXPN =expansion factor (current)
C			      ASTEP	  =increment of expansion factor: AEXPN =AEXP_0+N*ASTEP
C                     PARTW	=weight of a particle: total 'mass' per cell must be unity
c                                          even for LCDM or OCDM 
C			      AMPLT	=amplitude of perturbations inside the simulation box
C			      NBYTE	=record length in bytes
C
      SUBROUTINE InitValues(NBYTE,SCALEL)
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMinitialCR.h'
      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/

       Write (*,*) 'Would you like to use a model provided in cdm.fit '
       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  =')
       Do i=1,Levels
          AMPLT(i) =INPUT(' Amplitude of density fluctuations=')
       Enddo 
       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(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
       SCALEL=SCALEL/hubble   ! scale it to real megaparsecs      
       NBYTE = NPAGE*6*4
       Nspecies =INPUT(' Number of mass species   = ')
             W     = (FLOAT(NGRID)/FLOAT(NROW))**3
            PARTW = W 
      If(Nspecies.eq.0)Then  ! old constant mass resolution
              write(*,*) ' You use PM code with one particle set '
       Else    ! multiple mass resolution
              write(*,*) 
              write(*,*) ' You use multiple mass resolution'
       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 (*,*) ' Qscale=',QSCALE,SCALEL,ns

      RETURN
      END
C________________________________________Read parameters of the model
C                                 Om       = Omega_0
C                                 Omb     = Omega_baryon
C                                 Omnu  = Omega_neutrino
C                                 hsmall  = hubble constant
C                                 Par      = fitting parameters
C                                 Par(6) = 0  --- bbks-style+Hu&Sugiyama
C                                 Par(6) ne 0 --- holtzman-style
      SUBROUTINE MODEL(hsmall)
C---------------------------------------
      INCLUDE 'PMparameters.h'
      COMMON / TRUNCOM/ Om,Omb,Omc,Omnu,Par(6),ns,qqscaleb,QSCALE
      Real                 ns,INPUT
      Character             Header1*79

      OPEN(2,file='cdm.fit')
      Line1 =INPUT(' Enter Line Number in 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) 
      
      theta = 2.726/2.7  ! = T_cmb/2.7
      Ob0   =Omb/Om      ! Hu&Sugiyama fitting
      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
c      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))
      Return
      End
C------------------------------------------------
C                              Clear MASK for multiple mass resolution                             
C                               Initiat to level Nmin_lev
      SUBROUTINE ZeroMask(Nmin_lev)
C------------------------------------------------
      INCLUDE 'PMparameters.h'	    
      INCLUDE 'PMinitialCR.h'	  

C$omp  parallel do default(shared)
C$omp& private(MK3,MK2,MK1)
         DO  MK3 = 1,Nblocks
         DO  MK2 = 1,Nblocks
         DO  MK1 = 1,Nblocks
	        iMask(MK1,MK2,MK3) =Nmin_lev
         ENDDO
         ENDDO
         ENDDO
      RETURN
      END
C------------------------------------------------
C                       Set MASK for multiple mass resolution  
C                               Read particles with lagrangian coordinates
C                                                to set the mask                           
C                               Block   = 1 = low resolution
C                                       = 2 = higher resolution
C                                       = 3 = even higher resolution
C               After the mask is set, make maximum resolution NspecM
C
      SUBROUTINE SetMask(NspecM)
C------------------------------------------------
      INCLUDE 'PMparameters.h'	  
      INCLUDE 'PMinitialCR.h'	  
      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,ww,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
c             If(iMask(im,jm,km).ne.0)Then
c              write (*,*) ' Occupied block: ',im,jm,km,' Line=',Nlines
c              write (*,*) '                 Block=',iMask(im,jm,km)
c              Stop
c           EndIf
           
           iMask(im,jm,km) =NspecM     ! mark high resolution block
               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 in crest
                                ! 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)=NspecM ! change mask
                      Endif 
                  ENDDO
                ENDDO
                ENDDO
            GoTo 10
C                      Change the mask: iteratevely add layers of lower resolution
 50        write(*,*) ' Total Lines read=',Nlines
       DO kSpec =Nspecies-1,2,-1     !
         DO  MK3 = 1,Nblocks     ! find blocks adjasent to filled blocks
         DO  MK2 = 1,Nblocks     ! mark them with kSpec =Nspecies-1, -2, ..1
         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
                     If(iMask(i,j,k).lt.kSpec)iMask(i,j,k)=kSpec ! change mask
                  ENDDO
               ENDDO
               ENDDO
            EndIf
         ENDDO
         ENDDO
         ENDDO
       Enddo 
       If(NspecM.lt.Nspecies)Then ! change maximum resolution to NspecM

C$omp  parallel do default(shared)
C$omp& private(MK3,MK2,MK1)
         DO  MK3 = 1,Nblocks     
         DO  MK2 = 1,Nblocks     
         DO  MK1 = 1,Nblocks
            iM =      iMask(MK1,MK2,MK3)
            If(iM.gt.NspecM)iMask(MK1,MK2,MK3)=NspecM
         ENDDO
         ENDDO
         ENDDO
       Endif 

       Do i=1,Nspecies            ! 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.Nspecies.and.iM.ge.1)Then
               nLevel(iM) = nLevel(iM) +1 
            Else
               write (*,*) 'wrong mask: ',iM, ' block=',MK1,MK2,MK3
               Stop
            EndIf
         ENDDO
         ENDDO
         ENDDO

         Ntot =0
       Do i=Nspecies,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
C------------------------------------------------
C                             Make a realization of spectrum of perturbations:
C                             ALPHA  = normalization factor for displacements
C                             NRAND = seed for random numbers
C                             Lstart   = start on level Lstart (=1 or 2) 
      SUBROUTINE SPECTR(NRAND,Lstart)
C------------------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMinitialCR.h'
      REAL*8     SUMM, WI3,WJ3,WK3,WD,TS,TRX
c      Dimension SummA(100)
C						    Set spectrum
      If(Lstart.gt.Levels)Then
         write (*,*) ' Wrong range of levels in SPECTR: start/total=',
     &                   Lstart,Levels
         Return
      Endif 
      DO ML =Lstart,Levels
c.... Exemplar
c$DIR LOOP_PARALLEL
c$DIR LOOP_PRIVATE(MK3,MK2,MK1)
c.... Power Challenge
cc$DOACROSS LOCAL (MK3,MK2,MK1)      
C$omp  parallel do default(shared)
C$omp& private(MK3,MK2,MK1)
      DO  MK3 = 1,NROW
      DO  MK2 = 1,NROW
      DO  MK1 = 1,NROW
	     GRX(MK1,MK2,MK3,ML) =0.
	     GRY(MK1,MK2,MK3,ML) =0.
	     GRZ(MK1,MK2,MK3,ML) =0.
      ENDDO
      ENDDO
      ENDDO
      ENDDO

      DO Lv =Lstart,Levels     ! set spectrum for each level
c      DO Lv =2,2     ! set spectrum for each level
         SUMM = 0.

      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
c	        write (*,*) '  K=',k,' K3=',K3,' Nspec=',NSPEC
           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,Lv) = 0.
		         GRY(1,1,1,Lv) = 0.
		         GRZ(1,1,1,Lv) = 0.
	       ELSE
		  WD = (WI3 + WJ3 + WK3)*2**(2*Lv-2)
		  WK = SQRT(WD)
		  TS = TRUNF(WK) * GAUSS(NRAND)
		  TRX = TS / WD
		  GRX(iz, j, k,Lv) = i3 * TRX *isign
		  GRY( i,jz, k,Lv) = j3 * TRX *jsign
		  GRZ( i, j,kz,Lv) = k3 * TRX *ksign
		  SUMM = SUMM + TS**2
	       ENDIF
	    ENDDO
	   ENDDO
      ENDDO
      IF(SUMM.LE.0.)WRITE (*,*) ' Error!!! Summ over spectrum = 0'
         write (*,*) ' ==== SPECTR:  SUMM=',SUMM,' Level=',Lv
      ALPHA(Lv) = AMPLT(Lv) / SQRT(SUMM) *sqrt(8.)*1.087
      ENDDO

      RETURN
      END
C------------------------------------------------
C                             Make spectrum for Constrained Realizations
C                             
C                             
      SUBROUTINE SPectr_CR(NRAND)
C------------------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMinitialCR.h'
      REAL*8     SUMM, WI3,WJ3,WK3,WD,TS,TRX
      Dimension DENS_CR(NROW,NROW,NROW)
C						    Set spectrum
      SUMM = 0.
      WD       =0.
      open(40,file='dens_cr.dat',status='unknown')
      read(40,*) DENS_CR
      Do k=1,NROW
      Do j=1,NROW
      Do i=1,NROW
         SUMM =SUMM +dens_cr(i,j,k)
         WD      =WD +dens_cr(i,j,k)**2
      EndDo
      EndDo
      EndDo
      sigma = sqrt(MAX(1.d-10,WD/NROW**3))
      write (*,*) ' Average density=',SUMM/NROW**3
      write (*,*) ' RMS density=',sigma
       IFOUR = INT(ALOG(FLOAT(NROW))/ALOG(2.)+0.5)
       Write (*,*) ' Ifour =',IFOUR,' Nrow=',NROW

      CALL VECTOR_dens(IFOUR,dens_cr)

C$omp  parallel do default(shared)
C$omp& private(MK3,MK2,MK1)
      DO  MK3 = 1,NROW
      DO  MK2 = 1,NROW
      DO  MK1 = 1,NROW
	     GRX(MK1,MK2,MK3,1) =0.
	     GRY(MK1,MK2,MK3,1) =0.
	     GRZ(MK1,MK2,MK3,1) =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
	        write (*,*) '  K=',k,' K3=',K3,' Nspec=',NSPEC
           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,1) = 0.
		         GRY(1,1,1,1) = 0.
		         GRZ(1,1,1,1) = 0.
	       ELSE
		  WD = WI3 + WJ3 + WK3
		  WK = SQRT(WD)
		  TS = -dens_cr(i,j,k)
		  TRX = TS / WD
		  GRX(iz, j, k,1) = i3 * TRX *isign 
		  GRY( i,jz, k,1) = j3 * TRX *jsign
		  GRZ( i, j,kz,1) = k3 *TRX *ksign
		  SUMM = SUMM + TS**2
	       ENDIF
	    ENDDO
	   ENDDO
      ENDDO

      IF(SUMM.LE.0.)WRITE (*,*) ' Error!!! Summ over spectrum = 0'
         write (*,*) ' ==== SPECTR ===  SUMM=',SUMM
      ALPHA(1) = AMPLT(1) / SQRT(SUMM) *sqrt(8.)
      RETURN
      END
C------------------------------------------------
C                             Set to zero harmonics, which
C                             were included at previous level
C                             of resolution
C                             Lstart   = start on level Lstart (=1 or 2) 
      SUBROUTINE ZERO_SP
C------------------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMinitialCR.h'
      Logical       in_x,in_y,in_z,inside


      If(Levels.lt.2)Then
         write (*,*) ' Wrong range of levels in ZERO_SP: =',
     &                   Levels
         Return
      Endif 
      J_cos = NROW/4 +1   ! limit for cosine harmonics
      J_sin_1 = NROW/2+2  ! limits for sin harmonics
      J_sin_2 = NROW/4*3


C$omp  parallel do default(shared)
C$omp& private(K,J,I,in_z,in_y,in_x,inside,Lv)
      DO  K = 1,NROW
         If(K.le.J_cos.or.(K.ge.J_sin_1.and.K.le.J_sin_2))Then
            in_z =.true.
         else
            in_z =.false.
         Endif 
      DO  J = 1,NROW
         If(J.le.J_cos.or.(J.ge.J_sin_1.and.J.le.J_sin_2))Then
            in_y =.true.
         else
            in_y =.false.
         Endif 
      DO  I = 1,NROW
         If(I.le.J_cos.or.(I.ge.J_sin_1.and.I.le.J_sin_2))Then
            in_x =.true.
         else
            in_x =.false.
         Endif 
         inside = (in_x.and.in_y.and.in_z)
         If(inside)Then
c            write (*,*) I,J,K, in_x,in_y,in_z
            DO Lv =2,Levels
	           GRX(I,J,K,Lv) =0.
	           GRY(I,J,K,Lv) =0.
	           GRZ(I,J,K,Lv) =0.
            ENDDO
         Endif 
      ENDDO
      ENDDO
      ENDDO

      RETURN
      END

C------------------------------------------------
C		                        Define coordinates and velosities for
C                                  particles with resolution = Indx
C                                  Indx =1 = high resolution
C                                  Indx =2 = medium resolution
C                                  Icurrent = number of particles
      SUBROUTINE BLOCKS(Indx,Icurrent)
C------------------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMinitialCR.h'	  
      Dimension nLevel(10)
       Do i=1,Nspecies            ! number of particles at
          nLevel(i)=0                  ! different levels
       Enddo
       If(Indx.eq.Levels)Then
          Do i =1,Levels
             Delta(i)   = Float(NGRID)/NROW *2./2**i
             Delta2(i) = Delta(i)/2.
             iNROW(i) =NROW *2**i/2 
          Enddo 
       Endif 

      QFACT = FLOAT(NGRID)/FLOAT(NROW)
        DO  K = 1,Nblocks
           qq3 = Delta(1)*(K-1.) +1.
        DO  J = 1,Nblocks
           qq2 = Delta(1)*(J-1.) +1.
        DO  I  = 1,Nblocks
           qq1 = Delta(1)*(I-1.) +1.
	      iM =iMask(I,J,K)
C----------------------------------------
c                     
          If(iM.eq.Indx)Then 
          If(iM.eq.1)Then  ! Low resolution: one particle
               DX = GRX(I,J,K,1)
               DY = GRY(I,J,K,1)
               DZ = GRZ(I,J,K,1)
               Q3    = qq3 +Delta2(1)  ! scale them to box size
               Q2    = qq2 +Delta2(1)  
               Q1    = qq1 +Delta2(1) 
               Icurrent = Icurrent +1
               If(Icurrent.gt.Nmaxpart)Then
                    write (*,*)'Too many particles:',Icurrent,' block=',
     &                                            I,J,K, ' Level=',iM
                    STOP
               EndIf
               nLevel(Indx) =nLevel(Indx) +1
	         XPt(Icurrent) = Q1 - XCONS(1)*DX  + 1.E-3  +0.5 
	         YPt(Icurrent) = Q2 - XCONS(1)*DY  + 1.E-3  +0.5
	         ZPt(Icurrent) = Q3 - XCONS(1)*DZ  + 1.E-3  +0.5
	         VXt(Icurrent) =          VCONS(1)*DX
	         VYt(Icurrent) =          VCONS(1)*DY
	         VZt(Icurrent) =          VCONS(1)*DZ
           Else
C------------------------------- High resolution:
C                                          iM = current level ( > 1)
              Nlev =2**(iM-1)   ! number of particles in one dimenstion per Block
              Do kp =1,Nlev
               Q3    = qq3+Delta(iM)*(kp-1.) +Delta2(iM) ! scale them to box size
               kl       = INT((Q3-1.-Delta2(iM))/Delta(iM)) +1 ! inside large box 
                    if(kl.lt.1.or.kl.gt.NROW*2**(iM-1))
     &             write (*,*) ' error kl=',kl,' Lev=',iM,' k=',kp,K
               kll  = kl -INT((kl-1)/NROW)*NROW             ! inside small box of GRZ
                    if(kll.lt.1.or.kll.gt.NROW)
     &             write (*,*) ' error kll=',kll,' Lev=',iM,' kl=',kl

              Do jp =1,Nlev
               Q2    = qq2+Delta(iM)*(jp-1.) +Delta2(iM) 
               jl       = INT((Q2-1.-Delta2(iM))/Delta(iM)) +1 ! inside large box 
                    if(jl.lt.1.or.jl.gt.NROW*2**(iM-1))
     &             write (*,*) ' error jl=',jl,' Lev=',iM,' j=',jp,J
               jll  = jl -INT((jl-1)/NROW)*NROW             ! inside small box of GRY
                    if(jll.lt.1.or.jll.gt.NROW)
     &             write (*,*) ' error jll=',jll,' Lev=',iM,' jl=',jl

              Do ip =1,Nlev
               Q1    = qq1+Delta(iM)*(ip-1.) +Delta2(iM) 
               il       = INT((Q1-1.-Delta2(iM))/Delta(iM)) +1 ! inside large box 
                    if(il.lt.1.or.il.gt.NROW*2**(iM-1))
     &             write (*,*) ' error il=',il,' Lev=',iM,' j=',ip,I
               ill  = il -INT((il-1)/NROW)*NROW             ! inside small box of GRX
                    if(ill.lt.1.or.ill.gt.NROW)
     &             write (*,*) ' error ill=',ill,' Lev=',iM,' il=',il
               DX = XCONS(iM)*GRX(ill,jll,kll,iM)
               DY = XCONS(iM)*GRY(ill,jll,kll,iM)
               DZ = XCONS(iM)*GRZ(ill,jll,kll,iM)
               DvX = VCONS(iM)*GRX(ill,jll,kll,iM)
               DvY = VCONS(iM)*GRY(ill,jll,kll,iM)
               DvZ = VCONS(iM)*GRZ(ill,jll,kll,iM)
                 Do Lv =iM-1,1,-1 ! take contributions from lower levels

                    ii =INT((Q1-1.)/Delta(Lv)+1000.) -999 ! inside large box
                    jm   = INT((Q2-1.)/Delta(Lv)+1000.)-999 ! inside large box
                    km   = INT((Q3-1.)/Delta(Lv)+1000.)-999 ! inside large box

c                    ii =INT((Q1-1.-Delta2(Lv))/Delta(Lv)+1000.) -999 ! inside large box
c                    ii1     = ii +1
c                    dQ1    = Q1 - (Delta(Lv)*(ii-1)+Delta2(Lv)+1.)
c                    if(ii.lt.1)ii =ii+iNROW(Lv)
c                    if(ii.gt.iNROW(Lv))ii =ii - iNROW(Lv)
c                        if(ii.lt.1.or.ii.gt.iNROW(Lv))
c     &                 write (*,*) ' error ii=',ii,' Lev=',Lv,' j=',ip,I
                    iim  = ii  -INT((ii -1)/NROW)*NROW          ! inside small box of GRX
c                    iim1 = ii1-INT((ii1-1)/NROW)*NROW         
                   if(iim.lt.1.or.iim.gt.NROW)
     &             write (*,*) ' error iim=',iim,' Lev=',Lv,' ii=',ii

c                    jm   = INT((Q2-1.-Delta2(Lv))/Delta(Lv)+1000.)-999 ! inside large box
c                    jm1 = jm +1
c                    dQ2= Q2 - (Delta(Lv)*(jm-1)+Delta2(Lv)+1.)
c                    if(jm.lt.1)jm =jm+iNROW(Lv)
c                    if(jm.gt.iNROW(Lv))jm =jm - iNROW(Lv)
c                        if(jm.lt.1.or.jm.gt.iNROW(Lv))
c     &                 write (*,*) ' error jm=',jm,' Lev=',Lv,' j=',jp,J
                    jmm  = jm  -INT((jm -1)/NROW)*NROW          ! inside small box of GRY
c                    jmm1 = jm1-INT((jm1-1)/NROW)*NROW            
                   if(jmm.lt.1.or.jmm.gt.NROW)
     &             write (*,*) ' error jmm=',jmm,' Lev=',Lv,' jm=',jm

c                    km   = INT((Q3-1.-Delta2(Lv))/Delta(Lv)+1000.)-999 ! inside large box
c                    km1 = km +1
c                    dQ3 = Q3 - (Delta(Lv)*(km-1)+Delta2(Lv)+1.)
c                    if(km.lt.1)km =km+iNROW(Lv)
c                    if(km.gt.iNROW(Lv))km =km - iNROW(Lv)
c                        if(km.lt.1.or.km.gt.iNROW(Lv))
c     &                 write (*,*) ' error km=',km,' Lev=',Lv,' j=',jp,J
                    kmm  = km  -INT((km -1)/NROW)*NROW          ! inside small box of GRY
c                    kmm1 = km1-INT((km1-1)/NROW)*NROW            
                   if(kmm.lt.1.or.kmm.gt.NROW)
     &             write (*,*) ' error kmm=',kmm,' Lev=',Lv,' km=',km

c                   write (*,300) I,J,K,ip,jp,kp,Q1,Q2,Q3,
c     &                ill,jll,kll,ii,jm,km,iim,jmm,kmm
c 300               format(3x,' IJK=',3i3,' p=',3i3,' Q=',3f6.1,
c     &                  ' ill=',3i4,' i=',3i4,3i4)

                   uX  = GRX(iim ,jmm ,kmm ,Lv)
                   uY  = GRY(iim ,jmm ,kmm ,Lv)
                   uZ  = GRZ(iim ,jmm ,kmm ,Lv)

c                   dQ1 =dQ1/Delta(Lv)     ! normalize weights
c                   dQ2 =dQ2/Delta(Lv)
c                   dQ3 =dQ3/Delta(Lv)
c                   w111 = (1.-dQ1)  *(1.-dQ2)  *(1.-dQ3)
c                   w211 =  dQ1      *(1.-dQ2)  *(1.-dQ3)
c                   w121 = (1.-dQ1)  * dQ2      *(1.-dQ3)
c                   w221 =  dQ1      * dQ2      *(1.-dQ3)
c                   w112 = (1.-dQ1)  *(1.-dQ2)  * dQ3
c                   w212 =  dQ1      *(1.-dQ2)  * dQ3
c                   w122 = (1.-dQ1)  * dQ2      * dQ3
c                   w222 =  dQ1      * dQ2      * dQ3

c                                                        Interpolate displacements                   
c                   uX  = GRX(iim ,jmm ,kmm ,Lv)*w111 
c     &                   +GRX(iim1,jmm ,kmm ,Lv)*w211
c     &                   +GRX(iim ,jmm1,kmm ,Lv)*w121 
c     &                   +GRX(iim1,jmm1,kmm ,Lv)*w221
c     &                   +GRX(iim ,jmm ,kmm1,Lv)*w112 
c     &                   +GRX(iim1,jmm ,kmm1,Lv)*w212
c     &                   +GRX(iim ,jmm1,kmm1,Lv)*w122 
c     &                   +GRX(iim1,jmm1,kmm1,Lv)*w222
c                   uY  = GRY(iim ,jmm ,kmm ,Lv)*w111 
c     &                   +GRY(iim1,jmm ,kmm ,Lv)*w211
c     &                   +GRY(iim ,jmm1,kmm ,Lv)*w121 
c     &                   +GRY(iim1,jmm1,kmm ,Lv)*w221
c     &                   +GRY(iim ,jmm ,kmm1,Lv)*w112 
c     &                   +GRY(iim1,jmm ,kmm1,Lv)*w212
c     &                   +GRY(iim ,jmm1,kmm1,Lv)*w122 
c     &                   +GRY(iim1,jmm1,kmm1,Lv)*w222
c                   uZ  = GRZ(iim ,jmm ,kmm ,Lv)*w111 
c     &                   +GRZ(iim1,jmm ,kmm ,Lv)*w211
c     &                   +GRZ(iim ,jmm1,kmm ,Lv)*w121 
c     &                   +GRZ(iim1,jmm1,kmm ,Lv)*w221
c     &                   +GRZ(iim ,jmm ,kmm1,Lv)*w112 
c     &                   +GRZ(iim1,jmm ,kmm1,Lv)*w212
c     &                   +GRZ(iim ,jmm1,kmm1,Lv)*w122 
c     &                   +GRZ(iim1,jmm1,kmm1,Lv)*w222

                   DX  = DX  +XCONS(Lv)*uX       ! add contrubutions
                   DY  = DY  +XCONS(Lv)*uY
                   DZ  = DZ  +XCONS(Lv)*uZ
                   DvX = DvX +VCONS(Lv)*uX
                   DvY = DvY +VCONS(Lv)*uY
                   DvZ = DvZ +VCONS(Lv)*uZ
                Enddo                                         ! end Level

                 Icurrent = Icurrent +1
                 If(Icurrent.gt.Nmaxpart)Then
                    write (*,*)'Too many particles:',Icurrent,' block=',
     &                                    I,J,K,' Levels=',iM
                    STOP
                 EndIf
               nLevel(Indx) =nLevel(Indx) +1
	         XPt(Icurrent) = Q1 - DX  - 1.5999E-3  +0.5 
	         YPt(Icurrent) = Q2 - DY  - 1.5999E-3  +0.5
	         ZPt(Icurrent) = Q3 - DZ  - 1.5999E-3  +0.5
	         VXt(Icurrent) =          DvX
	         VYt(Icurrent) =          DvY
	         VZt(Icurrent) =          DvZ
               EndDo
              EndDo
              EndDo
              Endif   ! end iM =1

      EndIf         ! end iM = Indx 
      EndDo
      EndDo
      EndDo

           write (*,*) ' Particles=',Icurrent,' on  levels:',
     &                     (nLevel(i),i=1,Nspecies)
      RETURN
      END
C------------------------------------------------

      SUBROUTINE SETtest(K,NRAND)
C  place particles with small velocities in a sphere 
C  Radius = radius of the sphere in zero-level cell units
C  Xcentr  = position of the sphere
C  
C------------------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMinitialCR.h'	  
      COMMON / KINENRG/ SKINE,SX,SY,SZ,SX2,SY2,SZ2


      QFACT = FLOAT(NGRID)/FLOAT(NROW)
      XMAX  = FLOAT(NGRID) + 1.
      XSHF  = FLOAT(NGRID)
      Xcentr  = NGRID/2.
      Vss       =1.
      Radius = 5.
      write (*,*) ' Xcentr=',Xcentr
      DO  J = 1,NROW
	   DO  I = 1,NROW
            IJ    =      I +(J-1)*NROW
	      DX	      =      Radius*(2.*RANDd(NRAND)-1.)
	      DY	      =      Radius*(2.*RANDd(NRAND)-1.)
	      DZ	      =      Radius*(2.*RANDd(NRAND)-1.)
c	      DX	      =      0.05*I
c	      DY	      =      0.05*J
c	      DZ	      =      0.05*K
C           the factor 0.5 in the following 3 lines places unperturbed
c           particle in the middle of a cell. This reduces the shot noise
C          Small number 1.E-4 is to avoid placing some particles
C          exactly in cell nodes. Does not make much differentce.
	      XPAR(IJ) = Xcentr+ DX    + 1.E-4  +0.5 
	      YPAR(IJ) = Xcentr+ DY    + 1.E-4  +0.5 
	      ZPAR(IJ) = Xcentr+ DZ    + 1.E-4  +0.5 
	      VX(IJ)   =      Vss*GAUSS(NRAND)
	      VY(IJ)   =      Vss*GAUSS(NRAND)
	      VZ(IJ)   =      Vss*GAUSS(NRAND)
C                                Periodical boundary conditions
	       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		                        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(XCONS1,VCONS1,K)
C------------------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMinitialCR.h'	  
      COMMON / KINENRG/ SKINE,SX,SY,SZ,SX2,SY2,SZ2

      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	      =      XCONS1*GRX(I,J,K,1)
	      DY	      =      XCONS1*GRY(I,J,K,1)
	      DZ	      =      XCONS1*GRZ(I,J,K,1)
C           the factor 0.5 in the following 3 lines places unperturbed
c           particle in the middle of a cell. This reduces the shot noise
C          Small number 1.E-4 is to avoid placing some particles
C          exactly in cell nodes. Does not make much differentce.
	      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)   =      VCONS1*GRX(I,J,K,1)
	      VY(IJ)   =      VCONS1*GRY(I,J,K,1)
	      VZ(IJ)   =      VCONS1*GRZ(I,J,K,1)
C                                Periodical boundary conditions
	       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                                              				   FFT of the spectrum
      SUBROUTINE VECTOR_dens(IFOUR,dens_cr)
C------------------------------------------------
      INCLUDE 'PMparameters.h'	  
      INCLUDE 'PMinitialCR.h'	  
      dimension Uf(marr) , Vf(marr)
      dimension Qi(mf67) , jndx(mf67),dens_cr(NROW,NROW,NROW)

      jb1 =3
      jq   =IFOUR
      CALL SETF67(jb1,jq,jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)

C				       x-component
C.... Exemplar
C$DIR LOOP_PARALLEL
C$DIR LOOP_PRIVATE(j,i,Uf,Vf,jp,jsl,ll1,k2,k4)
C.... Power Challenge
Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 )
C$omp  parallel do default(shared)
C$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) =dens_cr(I,J,K)
	    ENDDO
	    CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf)
	    DO I=1,NROW
	       dens_cr(I,J,K) = Vf(I)
	    ENDDO
	   ENDDO
      ENDDO
      write (*,*) ' 1st is done'
C.... Exemplar
C$DIR LOOP_PARALLEL
C$DIR LOOP_PRIVATE(j,i,Uf,Vf,jp,jsl,ll1,k2,k4)
C.... Power Challenge
Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 )
C$omp  parallel do default(shared)
C$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) =dens_cr(I,J,K)
	    ENDDO
	    CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf)
	    DO J=1,NROW
	       dens_cr(I,J,K) = Vf(J)
	    ENDDO
	   ENDDO
      ENDDO
C.... Exemplar
C$DIR LOOP_PARALLEL
C$DIR LOOP_PRIVATE(i,k,Uf,Vf,jp,jsl,ll1,k2,k4)
C.... Power Challenge
Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 )
C$omp  parallel do default(shared)
C$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) =dens_cr(I,J,K)
	    ENDDO
	    CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf)
	    DO K=1,NROW
	       dens_cr(I,J,K) = Vf(K)/8.
	    ENDDO
	   ENDDO
      ENDDO
				       write (*,*) ' x is done'
                                       Return
                                       END
C------------------------------------------------
C                                              				   FFT of the spectrum
      SUBROUTINE VECTOR(IFOUR)
C------------------------------------------------
      INCLUDE 'PMparameters.h'	  
      INCLUDE 'PMinitialCR.h'	  
      dimension Uf(marr) , Vf(marr)
      dimension Qi(mf67) , jndx(mf67)

      jb1 =4
      jq   =IFOUR
      CALL SETF67(jb1,jq,jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)

      DO Lv=1,Levels
C				       x-component
C.... Exemplar
C$DIR LOOP_PARALLEL
C$DIR LOOP_PRIVATE(j,i,Uf,Vf,jp,jsl,ll1,k2,k4)
C.... Power Challenge
Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 )
C$omp  parallel do default(shared)
C$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) =GRX(I,J,K,Lv)
	    ENDDO
	    CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf)
	    DO I=1,NROW
	       GRX(I,J,K,Lv) = Vf(I)
	    ENDDO
	   ENDDO
      ENDDO
      write (*,*) ' 1st is done Level=',Lv
C.... Exemplar
C$DIR LOOP_PARALLEL
C$DIR LOOP_PRIVATE(j,i,Uf,Vf,jp,jsl,ll1,k2,k4)
C.... Power Challenge
Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 )
C$omp  parallel do default(shared)
C$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) =GRX(I,J,K,Lv)
	    ENDDO
	    CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf)
	    DO J=1,NROW
	       GRX(I,J,K,Lv) = Vf(J)
	    ENDDO
	   ENDDO
      ENDDO
C.... Exemplar
C$DIR LOOP_PARALLEL
C$DIR LOOP_PRIVATE(i,k,Uf,Vf,jp,jsl,ll1,k2,k4)
C.... Power Challenge
Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 )
C$omp  parallel do default(shared)
C$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) =GRX(I,J,K,Lv)
	    ENDDO
	    CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf)
	    DO K=1,NROW
	       GRX(I,J,K,Lv) = Vf(K)/8.
	    ENDDO
	   ENDDO
      ENDDO
				       write (*,*) ' x is done ',Lv
C				       Y-component
C.... Exemplar
C$DIR LOOP_PARALLEL
C$DIR LOOP_PRIVATE(j,i,Uf,Vf,jp,jsl,ll1,k2,k4)
C.... Power Challenge
Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 )
C$omp  parallel do default(shared)
C$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) =GRY(I,J,K,Lv)
	    ENDDO
	    CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf)
	    DO I=1,NROW
	       GRY(I,J,K,Lv) = Vf(I)
	    ENDDO
	   ENDDO
      ENDDO
C.... Exemplar
C$DIR LOOP_PARALLEL
C$DIR LOOP_PRIVATE(j,i,Uf,Vf,jp,jsl,ll1,k2,k4)
C.... Power Challenge
Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 )
C$omp  parallel do default(shared)
C$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) =GRY(I,J,K,Lv)
	    ENDDO
	    CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf)
	    DO J=1,NROW
	       GRY(I,J,K,Lv) = Vf(J)
	    ENDDO
	   ENDDO
      ENDDO
C.... Exemplar
C$DIR LOOP_PARALLEL
C$DIR LOOP_PRIVATE(k,i,Uf,Vf,jp,jsl,ll1,k2,k4)
C.... Power Challenge
Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 )
C$omp  parallel do default(shared)
C$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) =GRY(I,J,K,Lv)
	    ENDDO
	    CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf)
	    DO K=1,NROW
	       GRY(I,J,K,Lv) = Vf(K)/8.
	    ENDDO
	   ENDDO
      ENDDO
				       write (*,*) ' y is done'
C				       Z-component
C.... Exemplar
C$DIR LOOP_PARALLEL
C$DIR LOOP_PRIVATE(j,i,Uf,Vf,jp,jsl,ll1,k2,k4)
C.... Power Challenge
Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 )
C$omp  parallel do default(shared)
C$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) =GRZ(I,J,K,Lv)
	    ENDDO
	    CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf)
	    DO I=1,NROW
	       GRZ(I,J,K,Lv) = Vf(I)
	    ENDDO
	   ENDDO
      ENDDO
C.... Exemplar
C$DIR LOOP_PARALLEL
C$DIR LOOP_PRIVATE(j,i,Uf,Vf,jp,jsl,ll1,k2,k4)
C.... Power Challenge
Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 )
C$omp  parallel do default(shared)
C$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) =GRZ(I,J,K,Lv)
	    ENDDO
	    CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf)
	    DO J=1,NROW
	       GRZ(I,J,K,Lv) = Vf(J)
	    ENDDO
	   ENDDO
      ENDDO
C.... Exemplar
C$DIR LOOP_PARALLEL
C$DIR LOOP_PRIVATE(k,i,Uf,Vf,jp,jsl,ll1,k2,k4)
C.... Power Challenge
Cc$DOACROSS LOCAL ( k,j,i,Uf,Vf,jp,jsl,ll1,k2,k4 )
C$omp  parallel do default(shared)
C$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) =GRZ(I,J,K,Lv)
	    ENDDO
	    CALL FOUR67(jb1,jq,jp,jsl,ll1,k2,k7,Qi,jndx,Uf,Vf)
	    DO K=1,NROW
	       GRZ(I,J,K,Lv) = Vf(K)/8.
	    ENDDO
	   ENDDO
      ENDDO
	   write (*,*) ' z is done. Level=',Lv
      ENDDO
      RETURN
      END
C------------------------------------------
c                      
c                                                     Print statistics of spicies
      Subroutine Species_Print
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMinitialCR.h'	  

           write (*,300)Nspecies
           write (16,300)Nspecies
 300       format('---  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
C---------------------------------------------------
C                                  Read  current data from disk/tape,
C                                  Open files
C                                  Nrecl is the number of values in a record
C                                  Npage is the number of particles in a page
      SUBROUTINE RDTAPE
C---------------------------------------------------
      INCLUDE 'PMparameters.h'
      Character  Hd*5,Tail*4,Nm*40
      Hd  ='PMcrs'
      Tail='.DAT'
C                                     Open file on a tape
      OPEN(UNIT=9,FILE='PMcrd.DAT',
     +                FORM='UNFORMATTED', STATUS='UNKNOWN')
C                                 Read control information
C                                 and see whether it has proper
C                                 structure
      READ  (9,err=10,end=10) HEADER,
     +                  AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW,
     +                  TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0,
     +                  NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5
     +                   ,Ocurv,extras
      Ocurv =0.
      WRITE (*,100) HEADER,
     +                  AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW,
     +                  EKIN,EKIN1,EKIN2,
     +                  NROWC,NGRID,NRECL,Om0,Oml0,hubble,
     +                  Ocurv
      WRITE (16,100) HEADER,
     +                  AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW,
     +                  EKIN,EKIN1,EKIN2,
     +                  NROWC,NGRID,NRECL,Om0,Oml0,hubble,
     +                  Ocurv
100   FORMAT(1X,'Header=>',A45,/
     +            1X,' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/
     +            1X,' I =',I4,' WEIGHT=',F8.3,' Ekin=',3E12.3,/
     +            1X,' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I6,/
     +            1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/
     +            1x,' Omega_curvature=',F7.3)
      IF(NROW.NE.NROWC) THEN
         WRITE (*,*)
     +            ' NROW in PARAMETER and in TAPE-FILE are different'
      ENDIF
      IF(NGRID.NE.NGRIDC) THEN
         WRITE (*,*)
     +           ' NGRID in PARAMETER and in TAPE-FILE are different:'
         write (*,*) ' Ngrid=',NGRID,' NgridC=',NGRIDC
      ENDIF
C                                         Open work file on disk
 10   NBYTE = NRECL*4
      NACCES= NBYTE / nbyteword

         Nm =Hd//Char(48)//Tail
         iun=21
         write (*,*) ' File unit=',iun,' Name=',Nm
         OPEN(UNIT=iun,FILE=Nm,ACCESS='DIRECT',
     +	               STATUS='UNKNOWN',RECL=NACCES)

 
      REWIND 9
      RETURN
      END
C---------------------------------------------
C                       Write current data to  disk/tape
C
      SUBROUTINE WRTAPE
C----------------------------------------------
      INCLUDE 'PMparameters.h'
C                                       write header and control data
      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
      RETURN
      END
C--------------------------------------------------
C                             Write current PAGE of particles (x,v) to disk
C                             NRECL - length of ROW block in words
      SUBROUTINE WRIROW(IROW,Ifile)
C--------------------------------------------------
      INCLUDE 'PMparameters.h'
        WRITE (20+Ifile,REC=IROW) RECDAT
      RETURN
      END
cC--------------------------------------------------
cC                             Set Weights of particles for
cC                             fast access
c      SUBROUTINE SetWeight
cC--------------------------------------------------
c      INCLUDE 'PMparameters.h'
c      If(Nspecies.eq.0)Then  ! old  constant weights
c         N_particles =lspecies(1)
c         Do i=1,N_particles
c            iWeight(i) =PARTW
c         EndDo
c      Else
c         N_particles =lspecies(Nspecies)
c         jstart =1
c         Do j=1,Nspecies
c            jend =lspecies(j)
c            Do k=jstart ,jend
c               iWeight(k) =wspecies(j)
c            EndDo
c            jstart =jend
c         EndDo
c      EndIf
c      write (*,*) ' Set Weights for ',N_particles,' particles'
c      RETURN
c      END
C--------------------------------------------
C                            Read current PAGE of particles  from disk
C                            NRECL - length of ROW block in words
      SUBROUTINE GETROW(IROW,Ifile)
C--------------------------------------------
      INCLUDE 'PMparameters.h'
         READ  (20+Ifile,REC=IROW) RECDAT
      RETURN
      END
C------------------------------------------------
C				                                       random number generator
      FUNCTION RANDd(M)
C------------------------------------------------
      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
C--------------------------------------
C				normal random numbers
      FUNCTION GAUSS(M)
C--------------------------------------
      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

C---------------------------------- Read in variables      
      REAL FUNCTION INPUT(text)
C------------------------------------------------
      Character text*(*)
          write (*,'(A,$)')text
          read (*,*) x
          INPUT =x
      Return
      End
C--------------------------------------   Fourier Transform
      SUBROUTINE SETF67(JBC1,JQ1,jbc,jp,jsl,ll1,k2,k3,k4,k7,
     &                  Qi,jndx,Uf,Vf)
c     -----------------------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMinitialCR.h'	  

      dimension Uf(marr) , Vf(marr)
      dimension Qi(mf67) , jndx(mf67)

      PI=DATAN(1.D0)*4.
       JBC=JBC1
       JQ=JQ1
       IF(JBC.LT.3) GO TO 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) GO TO 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) GO TO 102
       GO TO 103
  104     RETURN
       END
C-----------------------------------------
       SUBROUTINE TFOLD(IS,L,ZZZ,k2)
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMinitialCR.h'	  


       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     CONTINUE
       RETURN
       END
C------------------------------------
       SUBROUTINE NEG(I1,I3,I2,jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMinitialCR.h'	  
       dimension Uf(marr) , Vf(marr)
       dimension Qi(mf67) , jndx(mf67)
 
       DO 100 K=I1,I3,I2
       Vf(K)=-Vf(K)
  100     CONTINUE
       RETURN
       END
C-------------------------------------
       SUBROUTINE REVNEG(jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMinitialCR.h'	  
       dimension Uf(marr) , Vf(marr)
       dimension Qi(mf67) , jndx(mf67)

       integer  jbc , jp, jsl , ll1 , k2 , k3 , k4 , k7

       DO 100 I=1,K7
       J=K3+1+I
       K=K4+1-I
       A=Vf(J)
       Vf(J)=-Vf(K)
       Vf(K)=-A
  100     CONTINUE
       RETURN
       END
C-------------------------------------
       SUBROUTINE ZEERO(L,jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMinitialCR.h'	  
       dimension Uf(marr) , Vf(marr)
       dimension Qi(mf67) , jndx(mf67)

       integer  jbc , jp, jsl , ll1 , k2 , k3 , k4 , k7

       DO 100 I=1,K2
       LI=L+I
       Uf(LI-1)=0.0
  100     CONTINUE
       RETURN
       END
C------------------------------------------
       SUBROUTINE TFOLD1(jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMinitialCR.h'	  
       dimension Uf(marr) , Vf(marr)
       dimension Qi(mf67) , jndx(mf67)

       integer  jbc , jp, jsl , ll1 , k2 , k3 , k4 , k7

       II=K2-1
       DO 100 I=1,II
       I1=JSL+I
       I2=LL1-I
       A=Uf(I1)
       Uf(I1)=A+Uf(I2)
       Uf(I2)=A-Uf(I2)
  100     CONTINUE
       RETURN
       END
C-------------------------------------
       SUBROUTINE KFOLD(jbc,jp,jsl,ll1,k2,k3,k4,k7,Qi,jndx,Uf,Vf)
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMinitialCR.h'	  
       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) GO TO 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) GO TO 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) GO TO 210
       I=I+1
  300     IS1=JSL
       IC1=LL1
       JS1=JS1/2
       IF(JS1.EQ.1) GO TO 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) GO TO 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) GO TO 330
       I=I+1
       IF(IS1.EQ.J5) GO TO 300
       GO TO 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) GO TO 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) GO TO 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) GO TO 400
       RETURN
       END

c      ------------------------------------------------------
       SUBROUTINE FOUR67(JBC1,JQ1,jp1,jsl1,ll11,k21,k71,
     &                   Qi,jndx,Uf,Vf)
c      ------------------------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMinitialCR.h'	  
       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
       GO TO (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
       GO TO 103
  102     K3=K3/2
       JQ=JQ-1
  103     N5=K3/4
       K7=K3/2
       N11=3*K7
       K31=K3+1
       GO TO(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

C---------------------------------------------
C                       Write current data to  disk/tape
C                        for PMstartM - multiple masses
      SUBROUTINE WriteData(Vscale,AexpV,Wtotal)
C----------------------------------------------
      INCLUDE 'PMparameters.h'
      COMMON / BUFF /	XPt(Nmaxpart),YPt(Nmaxpart),ZPt(Nmaxpart),
     +			VXt(Nmaxpart),VYt(Nmaxpart),VZt(Nmaxpart)
       COMMON / KINENRG/ SKINE,SX,SY,SZ,SX2,SY2,SZ2
C
      XMAX  = FLOAT(NGRID) + 1.
      XSHF  = FLOAT(NGRID)
      SKINE =0.
      Wtotal =0.
       jstart =1
       write (*,*) ' Nspecies=',Nspecies
       Do j =1,Nspecies
          vvx =0.
          vvy =0.
          vvz =0.
          vx2 =0.
          vy2 =0.
          vz2 =0.
          jend =lspecies(j)
          W     =wspecies(j)
          If(jend.gt.jstart)Then
             write (*,*)'  wrdata: jstart,jend=',jstart,jend
          Do i=jstart,jend
             vvx = vvx + VXt(i)
             vvy = vvy + VYt(i)
             vvz = vvz + VZt(i)
             vx2 = vx2 + VXt(i)**2
             vy2 = vy2 + VYt(i)**2
             vz2 = vz2 + VZt(i)**2
              SKINE    = SKINE + W*(VXt(i)**2 +VYt(i)**2 +VZt(i)**2)
              Wtotal =Wtotal +w 
          EndDo
           npp =max(jend-jstart+1,1)
           v2 = sqrt((vx2 +vy2+vz2)/npp)/AexpV*Vscale
           vvx =vvx/npp/AexpV*Vscale
           vvy =vvy/npp/AexpV*Vscale
           vvz =vvz/npp/AexpV*Vscale
           write (*,350) vvx,vvy,vvz,v2,jend-jstart+1,j
 350       format('       V(km/s)=',4g11.3, ' Npart=',i8,' Species=',i3)
          jstart =  jend+1
          Endif 
       EndDo
      Ibuff =0
      KROW =0
      write (*,*) ' lspecies(Nspecies)=',lspecies(Nspecies)
      Do i=1,lspecies(Nspecies)
C                                Periodical boundary conditions
	       IF(XPt(i).GT.XMAX)	      XPt(i)=XPt(i)-XSHF
	       IF(XPt(i).LE.1.)	            XPt(i)=XPt(i)+XSHF
	       IF(YPt(i).GT.XMAX)	      YPt(i)=YPt(i)-XSHF
	       IF(YPt(i).LE.1.)	            YPt(i)=YPt(i)+XSHF
	       IF(ZPt(i).GT.XMAX)	      ZPt(i)=ZPt(i)-XSHF
	       IF(ZPt(i).LE.1.)	            ZPt(i)=ZPt(i)+XSHF
             Ibuff              = Ibuff +1
             XPAR(Ibuff) = XPt(i)
             YPAR(Ibuff) = YPt(i)
             ZPAR(Ibuff) = ZPt(i)
             VX(Ibuff)      = VXt(i)
             VY(Ibuff)      = VYt(i)
             VZ(Ibuff)      = VZt(i)
             If(  XPAR(Ibuff).lt.1..or.YPAR(Ibuff).lt.1.
     &      .or.ZPAR(Ibuff).lt.1.)write (*,*) ' Error!!!',i, 
     &             XPAR(Ibuff),YPAR(Ibuff),ZPAR(Ibuff)
              If(  XPAR(Ibuff).ge.XMAX.or.YPAR(Ibuff).ge.XMAX
     &      .or.ZPAR(Ibuff).ge.XMAX)write (*,*) ' Error!!!',i, 
     &             XPAR(Ibuff),YPAR(Ibuff),ZPAR(Ibuff)
             If(  XPAR(Ibuff).ge.XMAX)Then
                XPAR(Ibuff) = XPAR(Ibuff) -1.e-4
                write (*,500) XPAR(Ibuff)
             endif 
             If(  YPAR(Ibuff).ge.XMAX)Then
                YPAR(Ibuff) = YPAR(Ibuff) -1.e-4
                write (*,500) YPAR(Ibuff)
             endif 
             If(  ZPAR(Ibuff).ge.XMAX)Then
                ZPAR(Ibuff) = ZPAR(Ibuff) -1.e-4
                write (*,500) ZPAR(Ibuff)
 500            format(' fixed boundary. It is now=',g14.7)
             endif 

            If(Ibuff.ge.NPAGE)Then
                KROW = KROW +1
                if(KROW/10*10.eq.KROW)
     &             write (*,*) ' Write page=',KROW,' i=',i,Ibuff
                CALL WRIROW(KROW,1)
                Ibuff     =0
             EndIf
          EndDo
          If(Ibuff.ne.0)Then    ! write last incomplete page
                KROW = KROW +1
                write (*,*) ' Write page=',KROW,' i_buff=',Ibuff
                CALL WRIROW(KROW,1)
                Ibuff     =0
          EndIf
      RETURN
      END
