C   ************** SURFRHOMULT.F ****************
C
C   CHARGE DENSITIES FOR SURFACE OF AN INHOMOGENEOUS SEMICONDUCTOR
C
C   VERSION 6.0 - FEB/11, DERIVED FROM surfrho-6.0
C
C   CONSTRUCT TABLE OF SURFACE CHARGE DENSITY VALUES
C
      SUBROUTINE SURFRHO(IAR,DELE,ESTART,NE,NEDIM,RHOSTAB)
      PARAMETER(NARDIM=2)
      DIMENSION RHOSTAB(NARDIM,NEDIM)
      COMMON/SURF/ISTK,TK,EN(NARDIM)
      DOUBLE PRECISION SUM
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
            RHOSTAB(IAR,I)=RHOS(IAR,EF1,DELE,EN,ISTK,TK)
200      CONTINUE
      ELSE
         NEN=NINT((EN(IAR)-ESTART)/DELE)+1
         RHOSTAB(IAR,NEN)=0.
         SUM=0.
         DO 300 I=NEN+1,NE
            EF1=(I-1)*DELE+ESTART
            SUM=SUM+SIGSUM(IAR,EF1,0.)
            RHOSTAB(IAR,I)=SUM*DELE
300      CONTINUE
         SUM=0.
         DO 310 I=NEN-1,1,-1
            EF1=(I-1)*DELE+ESTART
            SUM=SUM+SIGSUM(IAR,EF1,0.)
            RHOSTAB(IAR,I)=SUM*DELE
310      CONTINUE
      END IF
      RETURN
      END
C
C   TOTAL DENSITY OF SURFACE CHARGE
C
      FUNCTION RHOS(IAR,ENER,DELE)
      DOUBLE PRECISION SUM
      PARAMETER(NARDIM=2)
      COMMON/SURF/ISTK,TK,EN(NARDIM)
C
      SUM=0.
      E=EN(IAR)
      IF (ENER.LT.EN(IAR)) GO TO 200
100   E=E+DELE
      IF (E.GT.(ENER+10.*TK)) GO TO 900
      IF (ISTK.EQ.0) THEN
         SUM=SUM+SIGSUM(IAR,E,0.)*DELE
      ELSE
         SUM=SUM+SIGSUM(IAR,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(IAR,E,0.)*DELE
      ELSE
         SUM=SUM+SIGSUM(IAR,E,0.)*(1.-fd(e,ener,tk))*DELE
      END IF
      GO TO 200
900   RHOS=SUM
      RETURN
      END
C
      REAL*8 FUNCTION SIGSUM(IAR,EF,Pot)
      REAL*8 SIG
      SIGSUM=SIG(IAR,1,EF,Pot)+SIG(IAR,2,EF,Pot)
      RETURN
      END
C
C   FIND CHARGE-NEUTRALITY LEVEL
C
      SUBROUTINE ENFIND(IAR,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(IAR,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