C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MAXKPN=3000000)
      DIMENSION SS(9901,2),DD(9901,2)
      DIMENSION S(9901),D(9901),EA(MAXKPN)
	  DIMENSION N(3),P(3)
      READ(5,*) E0,DE,NNEE,ELECN
      DO 1 ISPIN=1,2
      DO 1 I=1,NNEE
      SS(I,ISPIN)=0.0D0
      DD(I,ISPIN)=0.0D0
    1 CONTINUE
      READ(21) NB1,NB2,NSPIN
	  WRITE(6,*) NB1,NB2,NSPIN
      WRITE(13,3001) NNEE,NB1,NB2,NSPIN,E0,DE,ELECN
 3001 FORMAT(4I5,3F10.5)
	  READ(21) P,N,JS
	  WRITE(6,*) N,JS
	  WRITE(6,*) P
	  DO 43 ISPIN=1,NSPIN
	  IB=0
      DO 43 JB=NB1,NB2
	  IB=IB+1
	  NB=JB
      READ(21) N64,NB,ISPIN,NN
	  READ(21) (EA(K),K=1,NN)
      WRITE(6,603) ISPIN,NB,N64,NN,EA(1),EA(NN)
  603 FORMAT(' DENSITY OF STATES FOR BAND ',3I4,I8,2F10.6)
      NXX=N(1)+1
      NYY=N(2)+1
      NZZ=N(3)+1
      CALL SDOS(EA,NNEE,E0,DE,NXX,NYY,NZZ,JS,D,S)
      LASTP=0
      JFIRST=0
      DO 31 I=1,NNEE
         DD(I,ISPIN)=DD(I,ISPIN)+D(I)
         SS(I,ISPIN)=SS(I,ISPIN)+S(I)
         IF(D(I).LT.0.000001) GO TO 31
         IF(JFIRST.EQ.0) JFIRST=I
         LASTP=I
   31 CONTINUE
      WRITE(6,660) ISPIN,NB,JFIRST,LASTP
     & ,E0+(JFIRST-1)*DE,E0+(LASTP-1)*DE
      WRITE(11,660) ISPIN,NB,JFIRST,LASTP
     & ,E0+(JFIRST-1)*DE,E0+(LASTP-1)*DE
  660 FORMAT(4I6,2F8.4)
      DO 32 I=JFIRST,LASTP
	     WRITE(11,600) I,E0+DE*(I-1),D(I),S(I)
   32 CONTINUE
      WRITE(14,660) ISPIN,NB,JFIRST,LASTP
      WRITE(14,3000) (D(I),I=JFIRST,LASTP)
      WRITE(14,3000) (S(I),I=JFIRST,LASTP)
      WRITE(6,*) ' End of SDOS for NB=',NB,' ISPIN=',ISPIN
C	  IF(LASTP.GT.0) STOP
   43 CONTINUE
	  DO 35 ISPIN=1,NSPIN
      DO 33 I=1,NNEE
      IF(DD(I,ISPIN).GT.0.0001)
     &WRITE(12,600) I,E0+DE*(I-1),DD(I,ISPIN),SS(I,ISPIN)
  600 FORMAT(I5,F7.4,F9.4,F9.4)
   33 CONTINUE
   35 CONTINUE
      IPFL=0
      DO 36 I=1,NNEE
      IF(SS(I,1)+SS(I,2).GT.ELECN.AND.IPFL.EQ.0) THEN
        EF=E0+DE*(I-2)+DE*(ELECN-(SS(I-1,1)+SS(I-1,2)))
     &  /((SS(I,1)+SS(I,2))-(SS(I-1,1)+SS(I-1,2)))
        WRITE(12,*) ' electron number=',ELECN,'  EF=',EF
        IPFL=1
      END IF
   36 CONTINUE
	  DO 37 ISPIN=1,NSPIN
      WRITE(13,3000) (DD(I,ISPIN),I=1,NNEE)
 3000 FORMAT(6E12.5)
   37 CONTINUE
      WRITE(13,3000) EF
      STOP
      END
      SUBROUTINE SDOS(EA,NE,E0,DE,IXMAX,IYMAX,IZMAX,JS,D,S)
C     ***  DENSITY OF STATES IS OBTAINED FOR A SINGLE BAND  ***
C     ***  BY LINEAR INTERPOLATION IN TETRAHEDORONS         ***
C     ***  EA  ENERGIES AT MESH POINTS
C     ***  NE  NUMBER OF ENERGY POINTS                      ***
C     ***  E0  STARTING ENERGY                              ***
C     ***  DE  STEP OF ENERGY                               ***
C     ***  IXMAX  NUMBER OF MESH POINTS IN B.Z.             ***
C     ***         IN X, Y, Z-DIRECTION FOR JS=3
C     ***         IN X,Y-DIRECTION FOR JS=2
C     ***         IN X-DIRECTION FOR JS=1
C     ***  IYMAX  NUMBER OF MESH POINTS IN B.Z.
C     ***         IN Y-DIRECTION FOR JS=1
C     ***  IZMAX  NUMBER OF MESH POINTS IN B.Z.             ***
C     ***         IN Z-DIRECTION FOR JS=1 OR 2
C     ***  JS=1  ORTHORHOMBIC
C     ***  JS=2  TETRAGONAL
C     ***  JS=3  CUBIC SYMMETRY
C     ***  D DENSITY OF STATES
C     ***  S NUMBER OF STATES
C     *** BY A. YANASE  JAN.1983                            ***
C     ***    MODIFIED FROM TDOSB
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION D(1),S(1)
      DIMENSION EA(1),EB(4),EC(8),ET(4),IQMAT(6,2)
      SAVE IQMAT
C----------------------------------------------------------
C      DATA IQMAT/2,2,5,3,3,5, 4,6,6,4,7,7/
C The data are changed at 2000.May.18
      DATA IQMAT/5,5,2,3,3,2, 7,6,6,7,4,4/
C----------------------------------------------------------
C      WRITE(6,*) JS,IXMAX,NE,EA(1),EA(20825)
      DO 250 I=1,NE
      D(I)=0.
  250 S(I)=0.
C     **  INTEGRATION OVER B.Z. STARTS    ***
      NTT=0
C      K=0
      JXM=IXMAX-1
      DO 353 IX=1,JXM
      IF(MOD(IX,10).EQ.0.OR.IX.EQ.JXM) write(6,*) ' IN SDOS IX=',IX
      LX=IX-1
      JYM=IX
      IF(JS.EQ.1) JYM=IYMAX-1
      DO 352 IY=1,JYM
      LY=IY-1
      JZM=IY
      IF(JS.NE.3) JZM=IZMAX-1
      DO 351 IZ=1,JZM
      LZ=IZ-1
C     ***  ENERGIES AT CUBE CORNERS  ***
      IF(JS.EQ.3) NI=(LX*(LX+1)*(LX+2))/6+(LY*(LY+1))/2+LZ+1
      IF(JS.EQ.2) NI=IZMAX*((LX*(LX+1))/2+LY)+LZ+1
      IF(JS.EQ.1) NI=IZMAX*(IYMAX*LX+LY)+LZ+1
      NN=NI
      NP=0
      DO 321 KX=1,2
      DO 322 KY=1,2
      DO 323 KZ=1,2
      NP=NP+1
      EC(NP)=EA(NN)
      NN=NN+1
  323 CONTINUE
      IF(JS.EQ.3) NN=NN-2+LY+KY
      IF(JS.EQ.2) NN=NN-2+IZMAX
      IF(JS.EQ.1) NN=NN-2+IZMAX
  322 CONTINUE
      IF(JS.EQ.3) NN=NI+(IX*(IX+1))/2
      IF(JS.EQ.2) NN=NI+IZMAX*(LX+1)
      IF(JS.EQ.1) NN=NI+IZMAX*IYMAX
  321 CONTINUE
C     *** SIX TETRAHEDORONS  ***
      ET(1)=EC(1)
      ET(4)=EC(8)
      IT=1
  325 CONTINUE
      DO 324 IP=1,2
      IQ=IQMAT(IT,IP)
      ET(IP+1)=EC(IQ)
  324 CONTINUE
      NTT=NTT+1
      DO 340 M=1,4
      EB(M)=ET(M)
  340 CONTINUE
C     ***  EB(1).LE.EB(2).LE.EB(3).LE.EB(4)  ***
      DO 342 M=2,4
      J=M
  341 J=J-1
      IF(EB(J).LE.EB(J+1))GOTO 342
      C=EB(J)
      EB(J)=EB(J+1)
      EB(J+1)=C
      IF(J.NE.1)GOTO 341
  342 CONTINUE
C     ***  JEPSON AND ANDERSON FORMULA  ***
      E1=EB(1)
      E2=EB(2)
      E3=EB(3)
      E4=EB(4)
      ISTA=(E1-E0)/DE
      ISTA=ISTA+2
      IF(ISTA.LT.1) ISTA=1
      IEND=(E4-E0)/DE
      IEND=IEND+1
      IF(IEND.GT.NE) IEND=NE
      IF(IEND.LT.0) IEND=0
C      write(6,*) LX,LY,LZ,ISTA,IEND,E1,E2,E3,E4
C      pause
      IF(ISTA.GT.IEND) GO TO 513
      D21=E2-E1
      D31=E3-E1
      D41=E4-E1
      D32=E3-E2
      D42=E4-E2
      D43=E4-E3
      DEL=E4+E3-E2-E1
      EM=(E4*E3-E2*E1)/DEL
      A21=1.0/DEL
      A22=(D31*D41+D32*D42)/(DEL*DEL)
      A20=0.0
      IF(D32.GT.1.0E-4) A20=DEL/(D31*D41*D32*D42)
      A1=0.0
      IF(D21.GT.1.0E-4) A1=1.0/(D21*D31*D41)
      A3=0.0
      IF(D43.GT.1.0E-4) A3=1.0/(D41*D42*D43)
      DO 514 I=ISTA,IEND
      E=E0+DE*FLOAT(I-1)
      IF(E.LE.E2) GO TO 512
      IF(E.GE.E3) GO TO 511
      DM=E-EM
      WD=DM*DM*A20
      WS=WD*DM
      S(I)=S(I)+A22+3.0*DM*A21-WS
      D(I)=D(I)+A21-WD
      GO TO 514
  511 D4=E-E4
      WD=D4*D4*A3
      WS=WD*D4
      S(I)=S(I)+1.0+WS
      D(I)=D(I)+WD
      GO TO 514
  512 D1=E-E1
      WD=D1*D1*A1
      WS=WD*D1
      S(I)=S(I)+WS
      D(I)=D(I)+WD
  514 CONTINUE
  513 IIEND=IEND+1
      IF(IIEND.GT.NE) GO TO 510
      DO 515 I=IIEND,NE
      S(I)=S(I)+1.0
  515 CONTINUE
  510 CONTINUE
C
      IT=IT+1
      IF(IT.EQ.7) GO TO 351
      IF(JS.EQ.2) GO TO 326
      IF(JS.EQ.1) GO TO 325
      IF((IT.EQ.2).AND.(IZ.EQ.IY)) IT=IT+2
      IF((IT.EQ.4).AND.(IY.EQ.IX)) GO TO 351
      IF((IT.EQ.6).AND.(IZ.EQ.IY)) GO TO 351
      GO TO 325
  326 IF((IT.EQ.4).AND.(IY.EQ.IX)) GO TO 351
      GO TO 325
C
  351 CONTINUE
  352 CONTINUE
  353 CONTINUE
      FACTS=1.0/FLOAT(NTT)
      FACT=3.0*FACTS
      DO 530 I=1,NE
      D(I)=FACT*D(I)
      S(I)=S(I)*FACTS
  530 CONTINUE
      RETURN
      END
     
