SUBROUTINE SUPDOS(ENER,IDIM1,IDIM2,N1,N2,NB,DK1,DK2,X,Y,Z,NX)
C    ************************************************************
C    *            DENSITY OF STATES OF A SUPERLATTICE           *
C    * OR ANY OTHER SYSTEMS WITH CYLINDRICAL BRILOUIN ZONE      *
C    * ENER : BAND STRUCTURE, FUNCTION OF K1 (= K PERPENDICULAR)*
C    * AND K2 (= MODULUS OF K PARALLEL TO THE INTERFACES)       *
C    * K1 AND K2 ARE DISCRETIZED ON A REGULAR MESH, SUCH THAT   *
C    * ENER(I1,I2,IB) = ENERGY OF BAND IB FOR                   *
C    * K1 = (I1-1)*DK1 , I1 = 1,...,N1 , (N1-1)*DK1 = PI/PERIOD *
C    * K2 = (I2-1)*DK2 , I2 = 1,...,N2                          *
C    * IDIM1, IDIM2 ARE THE TWO FIRST DIMENSIONS OF ENER        *
C    * NB IS THE NUMBER OF BANDS                                *
C    * X(1),...,X(NX) : TABULATED ENERGIES FOR DOS CALCULATION  *
C    * Y(1),...,Y(NX) : COMPUTED DOS (OUTPUT)                   *
C    * Z(1),...,Z(NX) : COMPUTED INTEGRATED DOS (OUTPUT)        *
C    * NORMALIZATION                                            *
C    * Z : NUMBER OF STATES PER CUBE OF THE RECIPROCAL UNITS OF *
C    *     DK1 AND DK2                                          *
C    * Y : SAME AS Z, PER UNIT OF ENERGY USED FOR ENER AND X    *
C    ************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 KPAR
      DIMENSION ENER(IDIM1,IDIM2,NB),X(NX),Y(NX),Z(NX)
      DIMENSION KPAR(3),ZONE(3)
      EQUIVALENCE (ZONE(1),E1),(ZONE(2),E2),(ZONE(3),E3)
      IF(NX.LE.0 .OR. NB.LE.0) RETURN
      M1=N1-1
      M2=N2-1
      IF(M1.LE.0 .OR. M2.LE.0) RETURN
C *** INITIALIZATION
      DO 17 IX=1,NX
      Y(IX)=0.0D0
   17 Z(IX)=0.0D0
C *** LOOP OVER THE BANDS
      DO 20 IB=1,NB
      DO 16 I1=1,M1
      DO 15 I2=1,M2
C *** CONSIDER TWO KINDS OF TRIANGULAR DECOMPOSITION OF THE CELLS
      IPASS=1
    1 GOTO(2,3,4,5),IPASS
      GOTO 15
    2 ZONE(1)=ENER(I1,I2,IB)
      ZONE(2)=ENER(I1,I2+1,IB)
      ZONE(3)=ENER(I1+1,I2+1,IB)
      KPAR(1)=(I2-1)*DK2
      KPAR(2)=I2*DK2
      KPAR(3)=KPAR(2)
      GOTO 6
    3 ZONE(1)=ENER(I1,I2,IB)
      ZONE(2)=ENER(I1+1,I2,IB)
      ZONE(3)=ENER(I1+1,I2+1,IB)
      KPAR(1)=(I2-1)*DK2
      KPAR(2)=KPAR(1)
      KPAR(3)=I2*DK2
      GOTO 6
    4 ZONE(1)=ENER(I1,I2,IB)
      ZONE(2)=ENER(I1+1,I2,IB)
      ZONE(3)=ENER(I1,I2+1,IB)
      KPAR(1)=(I2-1)*DK2
      KPAR(2)=KPAR(1)
      KPAR(3)=I2*DK2
      GOTO 6
    5 ZONE(1)=ENER(I1+1,I2+1,IB)
      ZONE(2)=ENER(I1+1,I2,IB)
      ZONE(3)=ENER(I1,I2+1,IB)
      KPAR(1)=I2*DK2
      KPAR(2)=(I2-1)*DK2
      KPAR(3)=KPAR(1)
C *** SORT THE CORNER ENERGIES
    6 DO 8 K=1,2
      DO 7 L=K+1,3
      IF(ZONE(L).GE.ZONE(K)) GOTO 7
      SAVE=ZONE(L)
      ZONE(L)=ZONE(K)
      ZONE(K)=SAVE
      SAVE=KPAR(L)
      KPAR(L)=KPAR(K)
      KPAR(K)=SAVE
    7 CONTINUE
    8 CONTINUE
      E21=E2-E1
      E31=E3-E1
      E32=E3-E2
C *** LOOP OVER THE ENERGIES
      DO 14 IX=1,NX
      E=X(IX)
      EE1=E-E1
      EE3=E3-E
      IF(EE1.EQ.0.0D0 .AND. EE3.EQ.0.0D0) GOTO 10
      IF(EE1.LE.0.0D0) GOTO 14
      IF(EE3.LE.0.0D0) GOTO 11
      IF(E.GT.E2) GOTO 9
      EE2=E2-E
      SAVE=EE2/E21+EE3/E31
      C1=EE1/(E21*E31)
      D1=C1*EE1*(1.0D0+SAVE)
      C1=C1*SAVE
      C2=(EE1/E21)**2/E31
      D2=C2*EE1
      C3=(EE1/E31)**2/E21
      D3=C3*EE1
      GOTO 12
    9 EE2=E-E2
      C1=(EE3/E31)**2/E32
      D1=1.0D0-EE3*C1
      C2=(EE3/E32)**2/E31
      D2=1.0D0-EE3*C2
      SAVE=EE1/E31+EE2/E32
      C3=EE3/(E31*E32)
      D3=1.0D0-C3*EE3*(1.0D0+SAVE)
      C3=C3*SAVE
      GOTO 12
   10 C1=1.0D+15
      C2=C1
      C3=C1
      D1=0.5D0
      D2=D1
      D3=D1
      GOTO 12
   11 D1=1.0D0
      D2=D1
      D3=D1
      GOTO 13
   12 Y(IX)=Y(IX)+C1*KPAR(1)+C2*KPAR(2)+C3*KPAR(3)
   13 Z(IX)=Z(IX)+D1*KPAR(1)+D2*KPAR(2)+D3*KPAR(3)
   14 CONTINUE
      IPASS=IPASS+1
      GOTO 1
   15 CONTINUE
   16 CONTINUE
   20 CONTINUE
C *** NORMALIZATION (SPIN IS NOT INCLUDED)
      PI=3.141592653589793238D0
      SURF=0.5D0*DK1*DK2
      FACC=0.25D0*SURF/PI**2
      FACD=FACC/3.0D0
      DO 25 IX=1,NX
      Y(IX)=Y(IX)*FACC
   25 Z(IX)=Z(IX)*FACD
      RETURN
      END