C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
      SUBROUTINE FUEINT
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'prmfit.f'
      COMMON/SPG2/IL,NG,IG(48),JV(2,3,48)
      COMMON/SCL/NATOM,NKAT,IPR,ISO,E0,DE,NNEE,NELEC,ENEND
      COMMON/ENR/KXM(3,MAXKPT),ICM(MAXKPT),IUDM(MAXKPT)
     &   ,NEIGM(MAXKPT),EIGM(MAXEIG,MAXKPT),NKPINT,NBAND,NB1,NB2
      COMMON/LCM/COMP(MAXEIG,MAXKPT,MAXCMP),NCOMP,
     $           NKCOMP(MAXCMP),NCCAR(MAXCMP)
      COMMON/NXYZ/NX,NY,NZ,N64,NBLOCK,JS
      COMMON/LATVC/LATOM(3,MAXLP),ALM(MAXLP),MP(MAXLP),NLAT,NSHEL,ICD
	  COMMON/FITB/B(MAXKPT,MAXEIG,2),FF(MAXKPT,MAXKPT),JPMAT(MAXKPT)
     &   ,NFSHEL
	  DIMENSION CRK(0:600),ENRR(MAXMESH,10,2),EK(25,2)
	  DIMENSION N(3),P(3)
	  WRITE(6,*) MAXEIG,NBAND,NB1,NB2,NX,NY,NZ,N64,NBLOCK
	  KRMESH=2*N64
	  IXM=NX
	  IYM=NY
	  IZM=NZ
	  P(1)=1.0D0
	  P(2)=1.0D0
	  P(3)=1.0D0
	  IF(IL.EQ.4.OR.IL.EQ.0) THEN
	     IZM=NZ/2
		 P(3)=0.5D0
      END IF
	  IF(IL.EQ.1) THEN
	     IXM=NX/2
		 IYM=NY/2
		 IZM=NZ/2
		 P(1)=0.5D0
		 P(2)=0.5D0
		 P(3)=0.5D0
		 KRMESH=N64
      END IF
	  N(1)=IXM
	  N(2)=IYM
	  N(3)=IZM
	  HPAI=DATAN(1.0D0)*2.0D0
	  MPX=(N64)/NX
	  MPY=(N64)/NY
	  MPZ=(N64)/NZ
	  write(6,*) NX,NY,NZ,MPX,MPY,MPZ,N64,KRMESH
	  IF(MOD(KRMESH,4).NE.0.OR.KRMESH.GT.1200) THEN
	     WRITE(6,*) ' THE NUMBER OF MESH INCORRECT',N64,KRMESH,ICD
		 STOP ' STOP IN FIEINT'
      END IF
	  MESH=KRMESH/4
	  CRK(0)=1.0D0
	  CRK(MESH*2)=-1.0D0
	  CRK(KRMESH)=1.0D0
	  CRK(MESH)=0.0D0
	  CRK(MESH*3)=0.0D0
	  write(6,*) MESH,KRMESH
      DO 11 I=1,MESH-1
         CRK(I)=DCOS((HPAI*I)/MESH)
		 CRK(MESH*2-I)=-CRK(I)
		 CRK(MESH*2+I)=-CRK(I)
		 CRK(KRMESH-I)=CRK(I)
   11 CONTINUE
      WRITE(6,601) (CRK(I),I=0,KRMESH)
  601 FORMAT(10F8.4)
	  WRITE(6,*) NBAND,NB1,NB2,NX,NY,NZ,N64,NBLOCK
	  WRITE(6,*) MPX,MPY,MPZ
      NSPIN=1
	  IF(ISO.EQ.2) NSPIN=2
	  write(6,*) ISO,NSPIN,JS,NB1,NB2
	  NN=0
      DO 21 IX=0,IXM
	     IF(IPR.GE.3) WRITE(6,*) ' in FUEINT IX=',IX,N64
		 IYMM=IYM
		 IF(JS.EQ.3.OR.JS.EQ.2) IYMM=IX
	  DO 22 IY=0,IYMM
C	     IF(IL.EQ.0) write(6,*) IX,IY,IX*MPX-IY*MPY,IX*MPX+IY*MPY
	     IZMM=IZM
	     IF(JS.EQ.3) IZMM=IY
	  DO 23 IZ=0,IZMM
	     NN=NN+1
		 IF(NN.GT.MAXMESH) THEN
		    WRITE(6,*) ' NN=',NN,' MAXMESH=',MAXMESH
			STOP ' AT MAXMESH'
		 END IF
C		 IF(IZ.EQ.IZMM.OR.(IX.EQ.IXM.AND.IY.EQ.IYMM)) THEN
C		 WRITE(6,*) IX*MPX,IY*MPY,IZ*MPZ,N64
C		 END IF
		 DO 41 ISPIN=1,NSPIN
		 IB=0
		 DO 41 JB=NB1,NB2
		 IB=IB+1
	     EK(IB,ISPIN)=0.0D0
   41    CONTINUE		 
		 IP=1
	     DO 31 JP=1,NFSHEL
		 J=JPMAT(JP)
		 IF(J.NE.1) IP=MP(J-1)+1
		 WW=0.0D0
		 DO 32 I=IP,MP(J)
		 IF(IL.GE.1) THEN
		    IRK=LATOM(1,I)*IX*MPX+LATOM(2,I)*IY*MPY
     &     +LATOM(3,I)*IZ*MPZ
         ELSE IF(IL.EQ.0) THEN
		    IRK=LATOM(1,I)*(IX*MPX-IY*MPY)
     &         +LATOM(2,I)*(IX*MPX+IY*MPY)
     &         +LATOM(3,I)*IZ*MPZ*2
         END IF
		 IF(IRK.LT.0) IRK=-IRK
   33    IF(IRK.LE.KRMESH) GO TO 34 
		    IRK=IRK-KRMESH
			GO TO 33
   34		WW=WW+CRK(IRK)   
   32    CONTINUE
         WW=WW/(MP(J)-IP+1)
C		 WRITE(6,*) WW
		 DO 42 ISPIN=1,NSPIN
		 IB=0
		 DO 42 JB=NB1,NB2
		 IB=IB+1
C		 write(6,*) JP,J,IB,ISPIN
		 EK(IB,ISPIN)=EK(IB,ISPIN)+WW*B(JP,IB,ISPIN)
   42    CONTINUE
   31    CONTINUE
         DO 46 ISPIN=1,NSPIN
         IB=0
         DO 46 JB=NB1,NB2
		 IB=IB+1
		 ENRR(NN,IB,ISPIN)=EK(IB,ISPIN)
   46    CONTINUE
C   		 WRITE(6,667) IX,IY,IZ,NN,IB,NSPIN
C     &			,(ENRR(NN,IIB,NSPIN),IIB=1,IB)
C  667     FORMAT(3I4,I7,2I4,10F7.4)
   23 CONTINUE
   22 CONTINUE	 
   21 CONTINUE
      WRITE(6,*) ' INTERPORATION N64=',N64,' NN=',NN
	  WRITE(6,*) N
	  WRITE(6,*) P
	  WRITE(6,*) ' JS=',JS
	  WRITE(21) NB1,NB2,NSPIN
	  WRITE(21) P,N,JS
	  DO 43 ISPIN=1,NSPIN
	  IB=0
      DO 43 JB=NB1,NB2
	  IB=IB+1
	  NB=JB
      WRITE(21) N64,NB,ISPIN,NN
	  WRITE(21) (ENRR(K,IB,ISPIN),K=1,NN)
C	  WRITE(6,*) (ENRR(K,IB,ISPIN),K=1,5)
   43 CONTINUE
      RETURN
      END
	  