C   ************** SURFRHO.F ****************
C
C   CHARGE DENSITIES FOR SURFACE OF A SEMICONDUCTOR
C   VERSION 1.0 - WRITTEN BY R. M. FEENSTRA, JAN 2004
C   VERSION 2.0 - AUG 2009, INCORPORATES TWO CHARGE DENSITY ARRAYS,
C                 FINDS CHARGE-NEUTRALITY LEVEL FOR COMBINED DISTRIBUTION
C                 USES DOUBLE PRECISION FOR INDIVIDUAL CHARGE DENSITIES
C
C   CONSTRUCT TABLE OF SURFACE CHARGE DENSITY VALUES
C
      SUBROUTINE SURFRHO(DELE,ESTART,NE,NEDIM,RHOS)
      DIMENSION RHOS(NEDIM)
      DOUBLE PRECISION SUM
      COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG,IINV,ISTK
      COMMON/SURF/EN
C
      IF (NE.GT.NEDIM) THEN
         WRITE(6,*) '*** ERROR - NE > NEDIM; PROGRAM HALTED'
         WRITE(6,*) 'TYPE ENTER TO CONTINUE'
         READ(5,*)
         STOP
      END IF
      IF (ISTK.EQ.1) THEN
         DO 200 I=1,NE
            EF1=(I-1)*DELE+ESTART
            RHOS(I)=SRHO(EF1,DELE)
200      CONTINUE
      ELSE
         NEN=NINT((EN-ESTART)/DELE)+1
         RHOS(NEN)=0.
         SUM=0.
         DO 300 I=NEN+1,NE
            EF1=(I-1)*DELE+ESTART
            SUM=SUM+SIGSUM(EF1,0.)
            RHOS(I)=SUM*DELE
300      CONTINUE
         SUM=0.
         DO 310 I=NEN-1,1,-1
            EF1=(I-1)*DELE+ESTART
            SUM=SUM+SIGSUM(EF1,0.)
            RHOS(I)=SUM*DELE
310      CONTINUE
      END IF
      RETURN
      END
C     
      FUNCTION SRHO(ENER,DELE)
      DOUBLE PRECISION SUM
      COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG,IINV,ISTK
      COMMON/SURF/EN
C
      SUM=0.
      E=EN
      IF (ENER.LT.EN) GO TO 200
100   E=E+DELE
      IF (E.GT.(ENER+10.*TK)) GO TO 900
      IF (ISTK.EQ.0) THEN
         SUM=SUM+SIGSUM(E,0.)*DELE
      ELSE
         SUM=SUM+SIGSUM(E,0.)*fd(e,ener,tk)*DELE
      END IF
      GO TO 100
200   E=E-DELE
      IF (E.LT.(ENER-10.*TK)) GO TO 900
      IF (ISTK.EQ.0) THEN
         SUM=SUM+SIGSUM(E,0.)*DELE
      ELSE
         SUM=SUM+SIGSUM(E,0.)*(1.-fd(e,ener,tk))*DELE
      END IF
      GO TO 200
900   SRHO=SUM
      RETURN
      END
C
      REAL*8 FUNCTION SIGSUM(EF,PHI)
      REAL*8 SIG
      SIGSUM=SIG(1,EF,PHI)+SIG(2,EF,PHI)
      RETURN
      END
C
C   FIND CHARGE-NEUTRALITY LEVEL
C
      SUBROUTINE ENFIND(EN1,EN2,EN0)
      REAL*8 SIGTMP,SIGSUM
C
      EN0=EN1
      ESTART=EN1
      NE=100000
      DELE=ABS(EN1-EN2)/FLOAT(NE)
      IF (DELE.EQ.0.) RETURN
      IF (EN2.LT.EN1) DELE=-DELE
      DO 100 IE=0,NE
         ENER=ESTART+IE*DELE
         SIGTMP=SIGSUM(ENER,0.)
         IF (DELE.GT.0.) THEN
            IF (SIGTMP.LE.0.D0) GO TO 200
         ELSE
            IF (SIGTMP.GE.0.D0) GO TO 200
         END IF
100   CONTINUE
200   EN0=ENER
      RETURN
      END