C ______________________________	   POWER SPECTRUM
C		                      25 april 1990
C		                     R. Kates (Potsdam), A. Klypin (NMSU)
C	     NROW = number of particles per row
C	     NGRID= number of cells on a side (normal grid)
C
C           Warning: comment out lines with declarations of
C                         common blocks /FOURAR/ and /F67COM/
C                         in files PMparameters.h        
C
       INCLUDE 'PMparameters.h'
       INCLUDE 'PMmesh.h'
       REAL     INPUT
C...................................................................
C			Read data and open files
      OPEN(17,FILE='Spectrum.DAT',STATUS='UNKNOWN')
      CALL RDTAPE

      write (*,*) ' RDTAPE is done'
      WRITE (17,100) HEADER,
     +                  AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,wspecies(1),
     +                  EKIN,
     +                  NROWC,NGRID,NRECL,Om0,Oml0,hubble,
     +                  Ocurv,lspecies(1),Nspecies
100   FORMAT(1X,'Header=>',A45,/
     +            1X,' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/
     +            1X,' I =',I4,' WEIGHT=',F8.3,' Ekin=',E12.3,/
     +            1X,' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I9,/
     +            1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/
     +            1x,' Omega_curvature=',F7.3,/
     +            1x,' Number of part.1specie=',i10,' Nspecies=',i3)
      Box = extras(100)
      If(Box.le.1.e-3)Then
          Box =INPUT(' Enter box size in Mpc/h   =')
      EndIf 
C
         CALL DENSIT    !   Find density
            write (*,*) ' Density is ok'
c      Call DensCR
         CALL DENTES(DELR) ! Statistics
         CALL POWER(Box)   ! Power spectrum

      STOP
      END

C-----------------------------------------------------
C		                                                  Compute mean density and rms
      SUBROUTINE DENTES(DELR)
C-----------------------------------------------------
       INCLUDE 'PMparameters.h'
       INCLUDE 'PMmesh.h'
      Dimension x(44),amass(44),imass(44)
      real*8 SUM,SUM2

      Data x/  0.,  0.2, 0.4, 0.6, 0.7,  0.8,  0.9,  1.,  1.1,   1.2,
     .         1.3, 1.5, 1.75, 2., 2.25, 2.5,  3.0,  3.5,  4.,    5.,
     .         7.,  10., 12., 15., 20.,  25.,  30.,  40.,  50.,  70.,
     .       100., 120.,150.,200.,300., 500., 700.,1000.,1200.,1500.,
     .       2.e3, 3.e3,5.e3,1.e4/
      Do i=1,44
         amass(i) =0.
         imass(i) =0
      EndDo
      
      SUM = 0.
      SUM2= 0.
      Nn  = 0
      Am  = 0.
C.... Open_MP
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE ( k,j,i)
C$OMP+REDUCTION(+:SUM,SUM2)
      DO K=1,NGRID
      DO J=1,NGRID
      DO I=1,NGRID
       SUM = SUM + FI(I,J,K)
       SUM2= SUM2+ FI(I,J,K)**2
c       Do l=44,1,-1
c          Dens =FI(I,J,K) +1.
c          If(Dens .GT. x(l))Then
c             Nn       =Nn +1
c             Am       =Am +Dens
c             imass(l) =imass(l) +1
c             amass(l) =amass(l) +Dens
c             GoTo 10
c          EndIf
c       EndDo
c 10    Continue
      ENDDO
      ENDDO
      ENDDO
      Total =(FLOAT(NGRID))**3
      DENM = SUM/Total
      DELR = DSQRT(SUM2/Total-DENM**2)
      DensAv=PARTW*float(NROW)**3/(float(NGRID))**3
      DensOne= float(NGRID)**3/float(NROW)**3
         WRITE (*,150)  DELR,DENM,DensAv,DensOne
         WRITE (17,150) DELR,DENM,DensAv,DensOne
 150     format(20('-'),' Density Distribution',20('-'),/
     &   '              Density is in units of average density',
     &                       ' in the Box',/20x,
     &                     ' RMS Delta Rho/Rho   =',G11.4,/20x,
     &                     ' Mean Delta Rho/Rho  =',G11.4,/20x,
     &                     ' Average density     =',G11.4,/20x,
     &                     ' One particle density=',G11.4)
c     &,/' bin    Rho     Rho   Sum(rho_i)/dRho/TotMass        Sum',/
c     &'      <center>< left > Normalized density     Number of cells ',/
c     &'      < bin  ><border>   Distribution         With this density')  
c      Do k=1,44-1
c         drho =x(k+1)-x(k)
c         rhom =x(k)+drho/2.
c         Sn   =drho*Total
c            If(Sn .lt.1.e-7)Sn=1.
c         write (*,20)  k,rhom,x(k),amass(k)/Sn,imass(k)
c         write (17,20) k,rhom,x(k),amass(k)/Sn,imass(k)
c 20      Format(i3,2f8.1,2x,G12.3,12X,i12)
c      EndDo
c         k =44
c         write (*,20)  k,x(k),x(k),amass(k),imass(k)
c         write (17,20) k,x(k),x(k),amass(k),imass(k)
      RETURN
      END

C----------------------------------------------------
C                                                     power spectrum   P(k)
C                                                     for given density field FI
      SUBROUTINE POWER(Box)
C----------------------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMmesh.h'
      PARAMETER ( Pi=3.141592653, NPOWER =NGRID)
      PARAMETER (NPROCS =64)
      DIMENSION DPOWER(NPOWER),iVolume(NPOWER),dk(NPOWER)
      DIMENSION DPW(NPOWER,NPROCS),iVol(NPOWER,NPROCS),
     &          dkk(NPOWER,NPROCS)

      dimension Zf(narr) , Yf(narr)
      dimension si(nf67) , indx(nf67)
      REAL*8    SIGMA,DPOWER,DPW

      NSPEC = NGRID/2
      If(Box.le.0..or.Box.gt.1.e+4)Then
         write (*,*) ' Wrong Box size=',Box,' STOP'
         Return
      EndIf
      SIGMA =0.
      Scale = 2.*Pi/Box
      DO M =1,NPOWER
         iVolume(M) =0
         DPOWER(M)  =0.
         dk(M)      =0.
         Do i=1,NPROCS
           iVol(M,i) =0
           DPW(M,i)     =0.
           dkk(M,i)     =0.
         EndDo 
      ENDDO
      ib1 = 3  
      iq  = int(alog(float(NGRID))/alog(2.)+0.5)
      FT  = (2.*NGRID)**3*float(NGRID)**3   
C                                               FFT of density field
      write (*,*) '     Do fft of rho along x-axis'
      call setf67(ib1,iq,ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf)

C.... Open_MP
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE ( k,j,i,zf,yf,ip,isl,l1,n2,n4 )
      DO K=1,NGRID
         DO J=1,NGRID
            DO I=1,NGRID
               Zf(I) =FI(I,J,K)
            ENDDO
	    call four67(ib1,iq,ip,isl,l1,n2,n7,
     &                 si,indx,Zf,Yf)
            DO I=1,NGRID
               FI(I,J,K) = Yf(I)
            ENDDO
c      If(J.lt.12.and.K.eq.1)write (*,'(2i4,16g12.3)')J,K,(Zf(I),I=1,16)
         ENDDO
      ENDDO

C.... Open_MP
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE ( k,j,i,zf,yf,ip,isl,l1,n2,n4 )
      DO K=1,NGRID
         DO I=1,NGRID
            DO J=1,NGRID
               Zf(J) =FI(I,J,K)
            ENDDO
            call four67(ib1,iq,ip,isl,l1,n2,n7,si,indx,
     &                        Zf,Yf)
            DO J=1,NGRID
               FI(I,J,K) = Yf(J)
            ENDDO
         ENDDO
      ENDDO

      write (*,*) '  Swap k<->i'
C.... Open_MP
C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP+PRIVATE ( k,j,i,aa)
      DO J=1,NGRID
      DO K=1,NGRID-1
            DO I=K+1,NGRID
               aa = FI(I,J,K)
               FI(I,J,K) =FI(K,J,I)
               FI(K,J,I) =aa
            ENDDO
         ENDDO
      ENDDO


      write (*,*) '     ....z'
C.... Open_MP
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE ( k,j,i,zf,yf,ip,isl,l1,n2,n4 )
      DO K=1,NGRID
         DO J=1,NGRID
            DO I=1,NGRID
               Zf(I) =FI(I,J,K)
            ENDDO
            call four67(ib1,iq,ip,isl,l1,n2,n7,si,indx,
     &                        Zf,Yf)
            DO I=1,NGRID
                  FI(I,J,K) =  Yf(I)**2/FT
            ENDDO
         ENDDO
      ENDDO

      FI(1,1,1) =0.
      kStep     =NGRID/NPROCS
      If(kstep*NPROCS.ne.NGRID)STOP '-- Wrong number of NPROCS'
      write (*,*) '    Power .... kStep=',kStep

C.... Open_MP
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE (iP,kStart,kEnd, k,j,i,MJ,MK,AMP,RK,MP)
C$OMP+REDUCTION(+:SIGMA)
      Do iP=1,NPROCS
         kStart =(iP-1)*kStep+1
         kEnd   =iP*kStep
      Do K=kStart,kEnd
         IF(K.LE.(NSPEC+1)) THEN
               MK =(K-1)**2 
         ELSE  
               MK =(K-1-NSPEC)**2
         ENDIF                                                                                              
      DO J=1,NGRID
         IF(J.LE.(NSPEC+1)) THEN
            MJ =(J-1)**2
         ELSE
            MJ =(J-1-NSPEC)**2
         ENDIF
         DO I=1,NGRID
               AMP = FI(I,J,K)
                  SIGMA = SIGMA +AMP
C                                          Find Abs(K)
                  IF(I.LE.(NSPEC+1)) THEN
                     RK =SQRT( FLOAT((I-1)**2     +MK+MJ) )
                  ELSE
                     RK =SQRT( FLOAT((I-1-NSPEC)**2+MK+MJ) )
                  ENDIF
C                                          Power Spectrum
                  MP =INT(RK +0.5)
                  IF(MP.LE.NPOWER.and.MP.gt.0)THEN
                     DPW(MP,iP)  =DPW(MP,iP) +AMP
                     iVol(MP,iP) =iVol(MP,iP)+1
                     dkk(MP,iP)  =dkk(MP,iP) +RK
                  ENDIF
            ENDDO
         ENDDO
      ENDDO
      ENDDO
      Do MP=1,NPOWER
      Do iP=1,NPROCS
          DPOWER(MP)  =DPOWER(MP) +DPW(MP,iP)
          iVolume(MP) =iVolume(MP)+iVol(MP,iP)
          dk(MP)      =dk(MP)     +dkk(MP,iP)
      EndDo 
      EndDo 

      WRITE (*,100) SQRT(SIGMA),Scale*NROW/2.
      WRITE (17,100) SQRT(SIGMA),Scale*NROW/2.
100   FORMAT(5x,'RMS DRho/Rho=',F7.3,' k_Nyquist for Particles=',f6.2,
     &        /4X,'n',2X,'k(h/Mpc)',1X,
     &       'Harmonics ','Power Spectrum(h^3)')
110   FORMAT(1x,I4,F8.3,I10,3X,5G12.4)
      DO I=1,NPOWER
         IF(iVolume(I).GT.0)THEN
            wk  = dk(I)/iVolume(I)*Scale
            Pk  = Box**3*DPOWER(I)/iVolume(I)
             
            WRITE (*,110)  I,wk,iVolume(I),Pk
            WRITE (17,110) I,wk,iVolume(I),Pk
         ENDIF
      ENDDO
      RETURN
      END
C*********************************************************************
C                                      Corr.function for given R
      FUNCTION FCORR(R)
       INCLUDE 'PMparameters.h'
       INCLUDE 'PMmesh.h'
      PARAMETER (  MMAX=NMAX-2)
      PARAMETER (PI=3.141592653)
 
      IF(R.LE.1.E-5)THEN
            WRITE (*,*) ' Fcorr(R): R can not be zero'
            STOP
      ENDIF
      QF    =2.*PI*R/FLOAT(NGRID)
      QSHF  =(FLOAT(NGRID)/2.)**2
      M     =NMAX +1
 
      SIGMA =0.
      FK    =SQRT(3.)*NMAX
      SIGMA =SIGMA + SIN(QF*FK) * FI(M,M,M)
      FK    =SQRT(2.)*NMAX
      SIGMA =SIGMA + SIN(QF*FK) *(FI(M,M,1)+FI(M,1,M)+FI(1,M,M))
      FK    =         NMAX
      SIGMA =SIGMA + SIN(QF*FK) *(FI(M,1,1)+FI(1,M,1)+FI(1,1,M))
      DO K1 =1,MMAX
         M1    =K1 +1
         L1    =K1 +1+NMAX
         K12   =K1**2
         SIGMA =SIGMA + SIN(QF*K1)*( FI(M1,1,1)+FI(L1,1,1)+
     .                               FI(1,M1,1)+FI(1,L1,1)+
     .                               FI(1,1,M1)+FI(1,1,L1))
         FK    =SQRT(K12+2.*QSHF)
         SIGMA =SIGMA + SIN(QF*FK)*( FI(M1,M,M)+FI(L1,M,M)+
     .                               FI(M,M1,M)+FI(M,L1,M)+
     .                               FI(M,M,M1)+FI(M,M,L1))
         FK    =SQRT(K12+   QSHF)
         SIGMA =SIGMA + SIN(QF*FK)*( FI(M1,M,1)+FI(L1,M,1)+
     .                               FI(M,M1,1)+FI(M,L1,1)+
     .                               FI(M1,1,M)+FI(L1,1,M)+
     .                               FI(M,1,M1)+FI(M,1,L1)+
     .                               FI(1,M1,M)+FI(1,L1,M)+
     .                               FI(1,M,M1)+FI(1,M,L1))
         DO K2 =1,MMAX
            M2    =K2 +1
            L2    =K2 +1+NMAX
            K22   =K2**2
            FK    =SQRT(Float(K12+ K22))
            SIGMA =SIGMA + SIN(QF*FK)*( FI(M1,M2,1)+FI(M1,L2,1)+
     .                                  FI(L1,M2,1)+FI(L1,L2,1)+
     .                                  FI(1,M1,M2)+FI(1,M1,L2)+
     .                                  FI(1,L1,M2)+FI(1,L1,L2)+
     .                                  FI(M1,1,M2)+FI(M1,1,L2)+
     .                                  FI(L1,1,M2)+FI(L1,1,L2))
            FK    =SQRT(K12+ K22 +QSHF)
            SIGMA =SIGMA + SIN(QF*FK)*( FI(M1,M2,M)+FI(M1,L2,M)+
     .                                  FI(L1,M2,M)+FI(L1,L2,M)+
     .                                  FI(M,M1,M2)+FI(M,M1,L2)+
     .                                  FI(M,L1,M2)+FI(M,L1,L2)+
     .                                  FI(M1,M,M2)+FI(M1,M,L2)+
     .                                  FI(L1,M,M2)+FI(L1,M,L2))
            DO K3 =1,MMAX
               M3    =K3 +1
               L3    =K3 +1+NMAX
               K32   =K3**2
               FK    =SQRT(Float(K12+ K22 +K32))
               SIGMA =SIGMA + SIN(QF*FK)*( FI(M1,M2,M3)+FI(M1,L2,M3)+
     .                                     FI(L1,M2,M3)+FI(L1,L2,M3)+
     .                                     FI(M1,M2,L3)+FI(M1,L2,L3)+
     .                                     FI(L1,M2,L3)+FI(L1,L2,L3))
            ENDDO
         ENDDO
      ENDDO
 
      FCORR  = SIGMA/QF
 
      RETURN
      END
C---------------------------------------
C                                                       Density Field    
      SUBROUTINE DENSIT
C---------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMmesh.h'
      	XN   =FLOAT(NGRID)+1.-1.E-7
	      YN   =FLOAT(NGRID)
C				       Subtract mean density
      DO M3=1,NGRID
       DO M2=1,NGRID
	      DO M1=1,NGRID
	        FI(M1,M2,M3) = -1.
	      END DO
       END DO
      END DO

      W     = (FLOAT(NGRID)/FLOAT(NROW))**3
      PARTW = W  
      WPAR = wspecies(1)
      N =0
      xmin =1.e+5
      xmax =-1.e+5
      ymin =1.e+5
      ymax =-1.e+5
      zmin =1.e+5
      zmax =-1.e+5

C				       Loop over particles
       If(Nspecies.eq.0)Then ! old case: 1 complete set
         Npages =NROW          ! Number of pages
         N_in_last=NPAGE       ! Number of particles in last page
      Else                             ! multiple masses
         N_particles =lspecies(Nspecies)   ! Total number of particles
         Npages = (N_particles -1)/NPAGE +1
         N_in_last=N_particles -NPAGE*(Npages-1)
      EndIf     
      write (*,*) Nspecies,
     + lspecies(Nspecies),wspecies(Nspecies)
      Nsp_current =1
      ff       = float(NGRIDC)
      ind = INT(log(float(NGRID)/ff)/log(2.)+10.)-10
      WPAR     = wspecies(1) *8.**ind
      Lspec_next=lspecies(1)+1

      DO IROW=1,Npages
            In_page =NPAGE
            If(IROW.eq.Npages)In_page =N_in_last
            If(IROW/50*50.eq.IROW)
     +        write (*,*) ' Read page=',IROW,' N=',N,
     +             ' Npages=',Npages,In_page
            iL = NPAGE*(IROW-1)
C				       Loop over particles
         READ  (21,REC=IROW) RECDAT
	   DO   IN=1,In_page
            N   = N +1
            If(N.ge.Lspec_next)Then   ! next specie
               write (*,*) ' Another specie. Current=',Nsp_current
               write (*,*) '     Particle=',N,WPAR,
     &                         ' lspecies=',lspecies(Nsp_current),
     &                         ' next=',lspecies(Nsp_current+1)
               Nsp_current=Nsp_current+1
               Lspec_next =lspecies(Nsp_current)+1
               WPAR     = wspecies(Nsp_current)
            EndIf
            X=(XPAR(IN)-1.)*NGRID/ff +1.
            Y=(YPAR(IN)-1.)*NGRID/ff +1.
            Z=(ZPAR(IN)-1.)*NGRID/ff +1.
           IF(X.ge.NGRID+1.)X=X-NGRID
           IF(Y.ge.NGRID+1.)Y=Y-NGRID
           IF(Z.ge.NGRID+1.)Z=Z-NGRID
           xmin =MIN(xmin,X)
           ymin =MIN(ymin,Y)
           zmin =MIN(zmin,Z)
           xmax =MAX(xmax,X)
           ymax =MAX(ymax,Y)
           zmax =MAX(zmax,Z)
	        I=INT(X)
	        J=INT(Y)
	        K=INT(Z)
            If(I.le.0)write (*,*) ' X:',X,Y,Z,I,' Irow=',IROW,IN,ifile
            If(J.le.0)write (*,*) ' Y:',X,Y,Z,I,' Irow=',IROW,IN,ifile
            If(K.le.0)write (*,*) ' Z:',X,Y,Z,I,' Irow=',IROW,IN,ifile
	        D1=X-FLOAT(I)
	        D2=Y-FLOAT(J)
	        D3=Z-FLOAT(K)
	        T1=1.-D1
	        T2=1.-D2
	        T3=1.-D3
	        T2W =T2*WPAR
	        D2W =D2*WPAR
	        I1=I+1
	           IF(I1.GT.NGRID)I1=1
	        J1=J+1
	           IF(J1.GT.NGRID)J1=1
	        K1=K+1
	           IF(K1.GT.NGRID)K1=1
		      FI(I ,J ,K ) =FI(I ,J ,K ) +T3*T1*T2W
		      FI(I1,J ,K ) =FI(I1,J ,K ) +T3*D1*T2W
		      FI(I ,J1,K ) =FI(I ,J1,K ) +T3*T1*D2W
		      FI(I1,J1,K ) =FI(I1,J1,K ) +T3*D1*D2W
    
		      FI(I ,J ,K1) =FI(I ,J ,K1) +D3*T1*T2W
		      FI(I1,J ,K1) =FI(I1,J ,K1) +D3*D1*T2W
		      FI(I ,J1,K1) =FI(I ,J1,K1) +D3*T1*D2W
		      FI(I1,J1,K1) =FI(I1,J1,K1) +D3*D1*D2W
	      ENDDO
        ENDDO 

      write (*,*) ' Number of particles read=', N
      write (*,*) ' X:',xmin,xmax
      write (*,*) ' Y:',ymin,ymax
      write (*,*) ' Z:',zmin,zmax
        RETURN
      END

C------------------------------------
c      Subroutine densCR
c      INCLUDE 'PMparameters.h'
c      INCLUDE 'PMmesh.h'
c      	XN   =FLOAT(NGRID)+1.-1.E-7
c	      YN   =FLOAT(NGRID)

c      SUMM = 0.
c      WD       =0.
c      open(40,file='dens_cr.dat',status='unknown')
c      read(40,*) FI
c      Do k=1,NROW
c      Do j=1,NROW
c      Do i=1,NROW
c         SUMM =SUMM +FI(i,j,k)
c         WD      =WD +FI(i,j,k)**2
c      EndDo
c      EndDo
c      EndDo
c      sigma = sqrt(MAX(1.d-10,WD/NROW**3))
c      write (*,*) ' Average density=',SUMM/NROW**3
c      write (*,*) ' RMS density=',sigma
c       IFOUR = INT(ALOG(FLOAT(NROW))/ALOG(2.)+0.5)
c       Write (*,*) ' Ifour =',IFOUR,' Nrow=',NROW

c       Return
c       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'
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=20,end=20) HEADER,
     +                  AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW,
     +                  TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0,
     +                  NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5
     +                   ,Ocurv,extras
      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=',I9,/
     +            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
      write (*,*)  ' start openning PMcrs0'
 10   NBYTE = NRECL*4
      NACCES= NBYTE / nbyteword

      OPEN(UNIT=21,FILE='PMcrs0.DAT',ACCESS='DIRECT',
     +	               STATUS='UNKNOWN',RECL=NACCES)


      REWIND 9
       write (*,*)  ' done openning PMcrs0'
      RETURN
 20   write (*,*) ' Error reading the header file: Abort'
      stop
      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(IBC1,IQ1,ibc,ip,isl,l1,n2,n3,n4,n7,
     &                  si,indx,Zf,Yf)
c     -----------------------------------------------------
       INCLUDE 'PMparameters.h'
       include 'PMmesh.h'
      dimension Zf(narr) , Yf(narr)
      dimension si(nf67) , indx(nf67)

      integer  ibc , ip, isl , l1 , n2 , n3 , n4 , n7
c      INCLUDE 'Control.param'
      PI=DATAN(1.D0)*4.
       IBC=IBC1
       IQ=IQ1
       IF(IBC.LT.3) GO TO 101
       IQ=IQ-1
  101      N3=2**IQ
       N7=N3/2
       N5=N3/4
       I=1
       INDX(I)=N5
       SI(I)=0.5*SQRT(2.)
       K=I
       I=I+1
  102     IL=I
       IF(I.EQ.N7) GO TO 104
  103     K1=INDX(K)/2
       INDX(I)=K1
       SI(I)=SIN(PI*FLOAT(K1)/FLOAT(N3))
       K1=N7-K1
       I=I+1
       INDX(I)=K1
       SI(I)=SIN(PI*FLOAT(K1)/FLOAT(N3))
       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,n2)
       INCLUDE 'PMparameters.h'
       include 'PMmesh.h'

       integer  n2

       DIMENSION ZZZ(NARR)
       IH2=N2/2-1
       DO 100 I=IS,IH2
       I1=I+L
       I2=N2-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,ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf)
       INCLUDE 'PMparameters.h'
       include 'PMmesh.h'
       dimension Zf(narr) , Yf(narr)
       dimension si(nf67) , indx(nf67)
 
       integer  ibc , ip, isl , l1 , n2 , n3 , n4 , n7

       DO 100 K=I1,I3,I2
       Yf(K)=-Yf(K)
  100     CONTINUE
       RETURN
       END
C-------------------------------------
       SUBROUTINE REVNEG(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf)
       INCLUDE 'PMparameters.h'
       include 'PMmesh.h'
       dimension Zf(narr) , Yf(narr)
       dimension si(nf67) , indx(nf67)

       integer  ibc , ip, isl , l1 , n2 , n3 , n4 , n7

       DO 100 I=1,N7
       J=N3+1+I
       K=N4+1-I
       A=Yf(J)
       Yf(J)=-Yf(K)
       Yf(K)=-A
  100     CONTINUE
       RETURN
       END
C-------------------------------------
       SUBROUTINE ZEERO(L,ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf)
       INCLUDE 'PMparameters.h'
       include 'PMmesh.h'
       dimension Zf(narr) , Yf(narr)
       dimension si(nf67) , indx(nf67)

       integer  ibc , ip, isl , l1 , n2 , n3 , n4 , n7

       DO 100 I=1,N2
       LI=L+I
       Zf(LI-1)=0.0
  100     CONTINUE
       RETURN
       END
C------------------------------------------
       SUBROUTINE TFOLD1(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf)
       INCLUDE 'PMparameters.h'
       include 'PMmesh.h'
       dimension Zf(narr) , Yf(narr)
       dimension si(nf67) , indx(nf67)

       integer  ibc , ip, isl , l1 , n2 , n3 , n4 , n7

       II=N2-1
       DO 100 I=1,II
       I1=ISL+I
       I2=L1-I
       A=Zf(I1)
       Zf(I1)=A+Zf(I2)
       Zf(I2)=A-Zf(I2)
  100     CONTINUE
       RETURN
       END
C-------------------------------------
       SUBROUTINE KFOLD(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf)
       INCLUDE 'PMparameters.h'
       include 'PMmesh.h'
       dimension Zf(narr) , Yf(narr)
       dimension si(nf67) , indx(nf67)

       integer  ibc , ip, isl , l1 , n2 , n3 , n4 , n7

       JS1=N2
       I=1
       J5=ISL+N2
       IS1=ISL
       IC1=L1
       JS1=JS1/2
       IF(JS1.NE.1) GO TO 200
       K1=INDX(I)
       SN=SI(I)
       IS0=IS1
       IS1=IS1+JS1
       IC0=IC1
       IC1=IC1+JS1
       ODD1=SN*(Zf(IC1)-Zf(IS1))
       Yf(K1+1)=Zf(IC0)+ODD1
       N3MK1=N3-K1
       Yf(N3MK1+1)=Zf(IC0)-ODD1
       IF(IBC.LT.3) GO TO 110
       ODD2=SN*(Zf(IC1)+Zf(IS1))
       N3PK1=N3+K1
       Yf(N3PK1+1)=Zf(IS0)+ODD2
       N4MK1=N4-K1
       Yf(N4MK1+1)=-Zf(IS0)+ODD2
  110     RETURN
  200     SN=SI(I)
       IS1=IS1+JS1
       IC1=IC1+JS1
       J3=IS1+JS1
  210     IS0=IS1-JS1
       IC0=IC1-JS1
       ODD1=SN*(Zf(IC1)-Zf(IS1))
       ODD2=SN*(Zf(IC1)+Zf(IS1))
       Zf(IC1)=Zf(IC0)-ODD1
       Zf(IS1)=-Zf(IS0)+ODD2
       Zf(IC0)=Zf(IC0)+ODD1
       Zf(IS0)=Zf(IS0)+ODD2
       IS1=IS1+1
       IC1=IC1+1
       IF(IS1.NE.J3) GO TO 210
       I=I+1
  300     IS1=ISL
       IC1=L1
       JS1=JS1/2
       IF(JS1.EQ.1) GO TO 400
  310     SN=SI(I)
       I=I+1
       CS=SI(I)
       IS1=IS1+JS1
       IC1=IC1+JS1
       J3=IS1+JS1
  320     IS0=IS1-JS1
       IC0=IC1-JS1
       ODD1=CS*Zf(IC1)-SN*Zf(IS1)
       ODD2=SN*Zf(IC1)+CS*Zf(IS1)
       Zf(IC1)=Zf(IC0)-ODD1
       Zf(IC0)=Zf(IC0)+ODD1
       Zf(IS1)=-Zf(IS0)+ODD2
       Zf(IS0)=Zf(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*Zf(IC1)-CS*Zf(IS1)
       ODD2=CS*Zf(IC1)+SN*Zf(IS1)
       Zf(IC1)=Zf(IC0)-ODD1
       Zf(IC0)=Zf(IC0)+ODD1
       Zf(IS1)=-Zf(IS0)+ODD2
       Zf(IS0)=Zf(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=INDX(I)
       SN=SI(I)
       I=I+1
       CS=SI(I)
       IS0=IS1
       IS1=IS1+JS1
       IC0=IC1
       IC1=IC1+JS1
       ODD1=CS*Zf(IC1)-SN*Zf(IS1)
       Yf(K1+1)=Zf(IC0)+ODD1
       N3MK1=N3-K1
       Yf(N3MK1+1)=Zf(IC0)-ODD1
       IF(IBC.LT.3) GO TO 410
       ODD2=SN*Zf(IC1)+CS*Zf(IS1)
       N3PK1=N3+K1
       Yf(N3PK1+1)=Zf(IS0)+ODD2
       N4MK1=N4-K1
       Yf(N4MK1+1)=-Zf(IS0)+ODD2
  410     IS1=IS1+1
       IC1=IC1+1
       K1=INDX(I)
       IS0=IS1
       IS1=IS1+JS1
       IC0=IC1
       IC1=IC1+JS1
       ODD1=SN*Zf(IC1)-CS*Zf(IS1)
       Yf(K1+1)=Zf(IC0)+ODD1
       N3MK1=N3-K1
       Yf(N3MK1+1)=Zf(IC0)-ODD1
       IF(IBC.LT.3) GO TO 420
       ODD2=CS*Zf(IC1)+SN*Zf(IS1)
       N3PK1=N3+K1
       Yf(N3PK1+1)=Zf(IS0)+ODD2
       N4MK1=N4-K1
       Yf(N4MK1+1)=-Zf(IS0)+ODD2
  420     IS1=IS1+1
       IC1=IC1+1
       I=I+1
       IF(IS1.NE.J5) GO TO 400
       RETURN
       END

c      ------------------------------------------------------
       SUBROUTINE FOUR67(IBC1,IQ1,ip1,isl1,l11,n21,n71,
     &                   si,indx,Zf,Yf)
c      ------------------------------------------------------
       INCLUDE 'PMparameters.h'
       include 'PMmesh.h'
       dimension Zf(narr) , Yf(narr)
       dimension si(nf67) , indx(nf67)

       integer  ibc , ip, isl , l1 , n2 , n3 , n4 , n7

c       INCLUDE 'Control.param'
       IBC=IBC1
       IQ=IQ1
       ip = ip1
       isl = isl1
       l1 = l11
       n2 = n21
       n7 = n71
       A5=0.5*SQRT(2.0)
       N4=2**IQ
       N3=N4
       GO TO (103,103,101,102),IBC
  101     Zf(1)=Zf(1)/2.0
       Zf(N3+1)=Zf(1)
       N2=N3

       CALL TFOLD(0,1,Zf,n2)

       N3=N3/2
       IQ=IQ-1
       GO TO 103
  102     N3=N3/2
       IQ=IQ-1
  103     N5=N3/4
       N7=N3/2
       N11=3*N7
       N31=N3+1
       GO TO(300,400,500,600),IBC
  300     Zf(N31)=0.0
       Zf(1)=0.0
       N2=N3
       DO 301 I=2,IQ

       CALL TFOLD(1,1,Zf,n2)

  301     N2=N2/2
       Yf(N7+1)=Zf(2)
       JF=N5
       ISL=1
       DO 302 IP=2,IQ
       L1=N2+1

       CALL ZEERO(1,ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf)
       CALL KFOLD(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf)

       I1=3*JF+1
       I2=4*JF
       I3=I1+(N2/2-1)*I2

       CALL NEG(I1,I3,I2,ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf)

       N2=N2+N2
  302     JF=JF/2
       RETURN
  400     Zf(1)=Zf(1)/2.0
       Zf(N31)=Zf(N31)/2.0
       N2=N3
       DO 401 I=2,IQ

       CALL TFOLD(0,N31-N2,Zf,n2)

  401     N2=N2/2
       L1=N31-N2
       A=Zf(L1)+Zf(L1+2)
       Yf(1)=-A-Zf(L1+1)
       Yf(N31)=-A+Zf(L1+1)
       Yf(N7+1)=Zf(L1)-Zf(L1+2)
       DO 402 IP=2,IQ
       ISL=N31-N2
       L1=ISL-N2
       idum = isl

       CALL ZEERO(idum,ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf)
       CALL KFOLD(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf)

  402     N2=N2+N2

       CALL NEG(1,N31,2,ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf)

       RETURN
  500     N2=N3
       L2=N4
       DO 501 IP=2,IQ

       CALL TFOLD(1,1,Zf,n2)
       CALL TFOLD(0,L2-N2+1,Zf,n2)

  501     N2=N2/2
       L1=L2-N2+1
       A=Zf(L1)+Zf(L1+2)
       Yf(N7+1)=2.0*(-Zf(L1)+Zf(L1+2))
       Yf(1)=2.0*(A+Zf(L1+1))
       Yf(N31)=2.0*(A-Zf(L1+1))
       Yf(N11+1)=2.0*Zf(2)
       DO 502 IP=2,IQ
       Zf(N2+1)=2.0*Zf(N2+1)
       ISL=N2+1

       CALL TFOLD1(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf)

       L1=L1-N2
       Zf(L1)=-2.0*Zf(L1)

       CALL KFOLD(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf)

  502     N2=N2+N2
       Yf(1)=Yf(1)*A5
       Yf(N31)=Yf(N31)*A5
       RETURN
  600     Zf(1)=Zf(1)*A5
       Zf(N31)=Zf(N31)*A5
       N2=N3
       DO 601 IP=2,IQ

       CALL TFOLD(0,N31-N2,Zf,n2)
       CALL TFOLD(1,N31,Zf,n2)

  601     N2=N2/2
       L1=N31-N2
       A=Zf(L1)+Zf(L1+2)
       Yf(1)=2.0*(A+Zf(L1+1))
       Yf(N31)=2.0*(A-Zf(L1+1))
       Yf(N7+1)=2.0*(-Zf(L1)+Zf(L1+2))
       Yf(N11+1)=2.0*Zf(N3+2)
       DO 602 IP=2,IQ
       ISL=N31+N2
       Zf(ISL)=2.0*Zf(ISL)

       CALL TFOLD1(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf)

       L1=L1-N2
       Zf(L1)=-2.0*Zf(L1)

       CALL KFOLD(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf)

  602     N2=N2+N2

       CALL REVNEG(ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf)

       idum = n3

       CALL NEG(2,idum,2,ibc,ip,isl,l1,n2,n3,n4,n7,si,indx,Zf,Yf)

       N2=N4

       CALL TFOLD(1,1,Yf,n2)

       Yf(N4+1)=0.0

       RETURN
       END

