C   ******************** SEMITIP2 ************************
C
C   E-FIELD COMPUTATIONS FOR TIP IN PROXIMITY TO SEMICONDUCTOR,
C   WITH CIRCULAR SYMMETRY
C   VERSION 1.0 - WRITTEN BY R. M. FEENSTRA, DEC 2002
C           2.0 - OCT/04, USE THIRD ORDER MATCHING AT INTERFACE
C                         USE VARIABLE GRID RADIALLY AND INTO SEMICOND.
C                         INCLUDE SURFACE STATES AND USE TIP MASKING ARRAY
C                         CHECK LAST TWO ITERATIONS FOR STOPPING CRITERION
C                         REVISE CHECKING OF ENERGY VALUES IN TABLE
C           2.1 - DEC/04, REVISE STOPPING CRITERION FOR GOLDEN SECTION SEARCHES
C
C   SEMI COMMON BLOCK:
C      EGAP=BAND GAP (EV)
C      ED=DONOR BINDING ENERGY (EV)
C      EA=ACCEPTOR BINDING ENERGY (EV)
C      ACB=CONDUCTION BAND EFFECTIVE MASS
C      AVB=VALENCE BAND EFFECTIVE MASS
C      CD=DONOR DENSITY (CM^-3)
C      CA=ACCEPTOR DENSITY (CM^-3)
C      EPSIL=DIELECTRIC CONSTANT OF SEMICONDUCTOR
C      TK=TEMPERATURE (EV)
C      IDEG=1 IF DEGENERATE DOPING, 0 OTHERWISE
C
C   INPUT AND OUTPUT PARAMETERS:
C      SEP=SEPARATION BETWEEN SEMICONDUCTOR AND END OF TIP (NM)
C      RAD=RADIUS OF HYPERBOLOID WHICH FORMS TIP (NM)
C      SLOPE=SLOPE OF TIP SHANK
C      ETAT=ETA PARAMETER FOR THE HYPERBOLIC TIP
C      AT=LOCATION OF FOCUS FOR THE HYPERBOLIC TIP
C      ETA1=ETA PARAMETER FOR THE HYPERBOLIC COORDINATE SYSTEM
C      A1=A2=LOCATION OF FOCUS FOR THE HYPERBOLIC COORDINATE SYSTEM
C      VAC=ARRAY OF POTENTIAL VALUES IN VACUUM
C      SEM=ARRAY OF POTENTIAL VALUES IN SEMICONDUCTOR
C      TIP=ARRAY WITH VALUES OF .TRUE. IF INSIDE TIP, .FALSE. OTHERWISE
C      VSINT=ARRAY OF POTENTIAL VALUES AT SEMICONDUCTOR SURFACE
C      R=ARRAY OF R VALUES
C      S=ARRAY OF Z VALUES IN SEMICONDUCTOR
C      DELV=ARRAY OF Z-SPACING VALUES IN VACUUM
C      DELR=R-SPACING (FOR SMALL R) IN VACUUM AND SEMICONDUCTOR
C      DELS=Z-SPACING (FOR SMALL Z) IN SEMICONDUCTOR
C      NRDIM=R-DIMENSION FOR VAC, SEM, AND VSINT ARRAYS
C      NVDIM=Z-DIMENSION FOR VAC ARRAY
C      NSDIM=Z-DIMENSION FOR SEM ARRAY
C      NR=CURRENT VALUE OF R-DIMENSION FOR VAC, SEM, AND VSINT ARRAYS
C      NV=CURRENT VALUE OF Z-DIMENSION FOR VAC ARRAY
C      NS=CURRENT VALUE OF Z-DIMENSION FOR SEM ARRAY
C      BIAS=BIAS OF TIP RELATIVE TO A POINT FAR INSIDE THE SEMICONDUCTOR
C      IWRIT >= 1 IF WANT WRITTEN OUTPUT, 0 OTHERWISE
C      NE=NUMBER OF POINTS IN CHARGE DENSITY TABLE
C      ESTART=STARTING ENERGY IN TABLE
C      DELE=ENERGY STEP SIZE IN TABLE
C      RHOC=CHARGE DENSITY TABLE
C      EF=FERMI-LEVEL POSITION
C      ITMAX=ARRAY OF MAXIMUM ITERATION SPECIFICATIONS
C      EP=ARRAY OF STOPPING CRITERIONS
C      IPMAX=MAXIMUM NUMBER OF SCALING LOOPS FOR GRID DOUBLING
C      PHI0=MAX SURFACE BAND BENDING
C      IERR=ERROR INDICATION (=1 IF GRID SIZE TOO SMALL, =2 IF ENERGY OUTSIDE TABLE)
C
C
C   SEMMIN USED TO FIND OPTIMAL UPDATE EQUATION FOR SOLN IN SEMICONDUCTOR
C
      FUNCTION SEMMIN(PHI,STEMP,DENOM,RHOC,NE,IERR)
      DIMENSION RHOC(NE)
      COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG
      COMMON/STIP/BIAS,DELR0,DELS0,DELETA,A1,ESTART,DELE,EF
      DATA EEP/1.80943E-20/
      ENER=EF-PHI
      IENER=NINT((ENER-ESTART)/DELE)+1
      IF (IENER.LT.1.OR.IENER.GT.NE) GO TO 100
	   RHO=RHOC(IENER)
	   GO TO 110
100	   RHO=RHOTOT(ENER,0.)
	   IERR=2
110	CONTINUE
      TEMP=STEMP-RHO*EEP/EPSIL
      SEMMIN=ABS(PHI-TEMP/DENOM)
      RETURN
      END
C
C   GOLDEN SECTION SEARCH ROUTINE (USED WITH SEMMIN)
C
      SUBROUTINE GSECT2(XMIN,XMAX,EP,STEMP,DENOM,RHOC,NE,IERR)
      DIMENSION RHOC(NE)
      DATA GS/0.3819660/
      IF (XMAX.EQ.XMIN) RETURN
      IF (EP.EQ.0.) RETURN
      IF (XMAX.LT.XMIN) THEN
         TEMP=XMAX
         XMAX=XMIN
         XMIN=TEMP
         END IF
      DELX=XMAX-XMIN
      XA=XMIN+DELX*GS
      FA=SEMMIN(XA,STEMP,DENOM,RHOC,NE,IERR)
      XB=XMAX-DELX*GS
      FB=SEMMIN(XB,STEMP,DENOM,RHOC,NE,IERR)
 100  DELXSAV=DELX
      IF (DELX.LT.EP) RETURN
      IF (FB.LT.FA) GO TO 200
      XMAX=XB
      DELX=XMAX-XMIN
      IF (DELX.EQ.DELXSAV) RETURN
      XB=XA
      FB=FA
      XA=XMIN+DELX*GS
      FA=SEMMIN(XA,STEMP,DENOM,RHOC,NE,IERR)
      GO TO 100
 200  XMIN=XA
      DELX=XMAX-XMIN
      IF (DELX.EQ.DELXSAV) RETURN
      XA=XB
      FA=FB
      XB=XMAX-DELX*GS
      FB=SEMMIN(XB,STEMP,DENOM,RHOC,NE,IERR)
      GO TO 100
      END
C
C   SURFMIN USED TO FIND OPTIMAL UPDATE EQUATION FOR SOLN ON SURFACE
C
      FUNCTION SURFMIN(PHI,RI,STEMP,DENOM,RHOS,NE,IERR)
      DIMENSION RHOS(NE)
      COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG
      COMMON/STIP/BIAS,DELR0,DELS0,DELETA,A1,ESTART,DELE,EF
      DATA EEP/1.80943E-20/
      ENER=EF-PHI
	DELE1=DELE
      IENER=NINT((ENER-ESTART)/DELE)+1
      IF (IENER.LT.1.OR.IENER.GT.NE) GO TO 100
	   RHO=RHOS(IENER)
	   GO TO 110
100	   RHO=SRHO(ENER,DELE1)
	   IERR=2
110	CONTINUE
C   ADD STATEMENT BELOW TO RESTRICT REGION OF SURFACE CHARGE
C	IF (RI.LT.0.) RHO=0.
      TEMP=STEMP-RHO*EEP*1.E7
      SURFMIN=ABS(PHI-TEMP/DENOM)
      RETURN
      END
C
C   GOLDEN SECTION SEARCH ROUTINE (USED WITH SURFMIN)
C
      SUBROUTINE GSECT3(XMIN,XMAX,EP,RI,STEMP,DENOM,RHOS,NE,IERR)
      DIMENSION RHOS(NE)
      DATA GS/0.3819660/
      IF (XMAX.EQ.XMIN) RETURN
      IF (EP.EQ.0.) RETURN
      IF (XMAX.LT.XMIN) THEN
         TEMP=XMAX
         XMAX=XMIN
         XMIN=TEMP
         END IF
      DELX=XMAX-XMIN
      XA=XMIN+DELX*GS
      FA=SURFMIN(XA,RI,STEMP,DENOM,RHOS,NE,IERR)
      XB=XMAX-DELX*GS
      FB=SURFMIN(XB,RI,STEMP,DENOM,RHOS,NE,IERR)
 100  DELXSAV=DELX
      IF (DELX.LT.EP) RETURN
      IF (FB.LT.FA) GO TO 200
      XMAX=XB
      DELX=XMAX-XMIN
      IF (DELX.EQ.DELXSAV) RETURN
      XB=XA
      FB=FA
      XA=XMIN+DELX*GS
      FA=SURFMIN(XA,RI,STEMP,DENOM,RHOS,NE,IERR)
      GO TO 100
 200  XMIN=XA
      DELX=XMAX-XMIN
      IF (DELX.EQ.DELXSAV) RETURN
      XA=XB
      FA=FB
      XB=XMAX-DELX*GS
      FB=SURFMIN(XB,RI,STEMP,DENOM,RHOS,NE,IERR)
      GO TO 100
      END
C
C   MAIN PROGRAM
C
      SUBROUTINE SEMITIP2(SEP,RAD,SLOPE,ETAT,AT,ETA1,A2,VAC,TIP,SEM,
     &VSINT,R,S,DELV,DELR1,DELS1,NRDIM,NVDIM,NSDIM,NR1,NV1,NS1,BIAS1,
     &IWRIT,NE,ESTART1,DELE1,RHOC,RHOS,EF1,ITMAX,EP,IPMAX,PHI0,IERR)
C
      DIMENSION VAC(2,NRDIM,NVDIM),SEM(2,NRDIM,NSDIM),VSINT(2,NRDIM),
     &R(NRDIM),DELR(NRDIM),DELV(NRDIM),S(NSDIM),DELS(NSDIM),ITMAX(10),
     &EP(10),RHOC(NE),RHOS(NE)
	LOGICAL TIP(NRDIM,NVDIM)
      COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG
      COMMON/SIZE/NR,NV,NS
      COMMON/STIP/BIAS,DELR0,DELS0,DELETA,A1,ESTART,DELE,EF
	DATA PI/3.141592654/
      BIAS=BIAS1
      ESTART=ESTART1
      DELE=DELE1
	EF=EF1
      NR=NR1
      NV=NV1
      NS=NS1
	DELR0=DELR1
	DELS0=DELS1
      IERR=0
C
C   CONSTRUCT THE TIP AND VACUUM GRID AND INITIALIZE POTENTIAL
C
      ETAT=1./SQRT(1.+1./SLOPE**2)
      AT=RAD*SLOPE**2/ETAT
      IF (SEP.GT.AT*ETAT) GO TO 50
C   CASE 1 - BASIS TIP HAS SAME ETA & A AS ACTUAL TIP, BUT DIFFERENT SEPARATION
         ETA1=ETAT
         A1=AT
         A2=A1
         GO TO 60
C   CASE 2 - BASIS TIP HAS SAME SEPARATION AS ACTUAL TIP, BUT DIFFERENT ETA & A
50       ETA1=1./SQRT(1.+RAD/SEP)
         A1=SQRT(SEP*(RAD+SEP))
	   A2=A1
60    CONTINUE
      DELETA=ETA1/FLOAT(NV)
      CETA1=ALOG((1.+ETA1)/(1.-ETA1))
      IF (IWRIT.EQ.0) GO TO 100
         WRITE(6,*) 'ACTUAL ETA, A, S =',ETAT,AT,SEP
         WRITE(6,*) 'COORD  ETA, A, S =',ETA1,A1,ETA1*A1
         WRITE(6,*) 'NR,NS,NV =',NR,NS,NV
         WRITE(6,*) 'DELR,DELS,DELV =',DELR0,DELS0,A1*DELETA
         WRITE(16,*) 'ACTUAL ETA, A, S =',ETAT,AT,SEP
         WRITE(16,*) 'COORD  ETA, A, S =',ETA1,A1,ETA1*A1
         WRITE(16,*) 'NR,NS,NV =',NR,NS,NV
         WRITE(16,*) 'DELR,DELS,DELV =',DELR0,DELS0,A1*DELETA
100   CONTINUE
      IF (NR.GT.NRDIM.OR.NV.GT.NVDIM.OR.NS.GT.NSDIM) GO TO 600
      DO 150 I=1,NR
	   R(I)=(2*NR*DELR0/PI)*TAN(PI*(I-0.5)/(2.*NR))
	   DELR(I)=DELR0/(COS(PI*(I-0.5)/(2.*NR)))**2
         DELV(I)=SQRT(A1**2+R(I)**2)*ETA1/FLOAT(NV)
         DO 120 J=1,NV-1
            X2M1=(R(I)/A1)**2
            XSI=DSQRT(1.D0+X2M1)
            ETA=J*ETA1/FLOAT(NV)
            Z=A1*XSI*ETA
            RP=A1*SQRT(X2M1*(1.-ETA**2))
            IF (Z.GT.(ETAT*(SQRT(RP**2/(1.-ETAT**2)+AT**2)-AT)+SEP
     &          -P(RP))) GO TO 110
C   ADD STATEMENT BELOW FOR FLAT TIP APEX
C            IF (Z.GT.(SEP-P(RP)).AND.RP.LT.10.) GO TO 110
            VAC(1,I,J)=0.
            VAC(2,I,J)=0.
	      TIP(I,J)=.FALSE.
            GO TO 120
110         VAC(1,I,J)=BIAS
            VAC(2,I,J)=BIAS
	      TIP(I,J)=.TRUE.
120      CONTINUE
         VAC(1,I,NV)=BIAS
         VAC(2,I,NV)=BIAS
	   TIP(I,NV)=.TRUE.
150   CONTINUE
      DO 160 I=1,NR
         VSINT(1,I)=0.
         VSINT(2,I)=0.
160   CONTINUE
      DO 180 J=1,NS
	   IF (J.EQ.NS) GO TO 165
	   S(J)=(2*NS*DELS0/PI)*TAN(PI*J/(2.*NS))
	   DELS(J)=DELS0/(COS(PI*J/(2.*NS)))**2
165         DO 170 I=1,NR
            SEM(1,I,J)=0.
            SEM(2,I,J)=0.
170      CONTINUE
180   CONTINUE
      IF (.NOT.TIP(1,1)) GO TO 190
         IERR=1
         WRITE(6,*) '*** ERROR - VACUUM GRID SPACING TOO LARGE'
         WRITE(16,*) '*** ERROR - VACUUM GRID SPACING TOO LARGE'
190   CONTINUE
	WRITE(6,*) 'LARGEST RADIUS, DEPTH =',R(NR),S(NS-1)
	WRITE(16,*) 'LARGEST RADIUS, DEPTH =',R(NR),S(NS-1)
      IF (IERR.EQ.1) GO TO 600
C
C   INITIAL GUESS
C
      DO 250 I=1,NR
         DO 210 J=NV,1,-1
            IF (.NOT.TIP(I,J)) GO TO 220
210      CONTINUE
220      JSAV=J
         DO 230 J=JSAV,1,-1
            ETA=J*DELETA*FLOAT(NV)/FLOAT(JSAV+1)
            VAC(1,I,J)=BIAS*ALOG((1.+ETA)/(1.-ETA))/CETA1
            VAC(2,I,J)=VAC(1,I,J)
230      CONTINUE
250   CONTINUE
      PHI0=(9.*VSINT(1,1)-VSINT(1,2))/8.
C
C   ITERATE POISSON'S EQN
C
      DO 590 IP=1,IPMAX
C
         IF (IWRIT.NE.0) WRITE(6,*) 'SOLUTION #',IP
         IF (IWRIT.NE.0) WRITE(16,*) 'SOLUTION #',IP
         ITM=ITMAX(IP)
         EPI=EP(IP)
         CALL ITER2(VAC,TIP,SEM,VSINT,R,DELR,DELV,DELS,NRDIM,NVDIM,
     &              NSDIM,EPI,ITM,PHI0,IWRIT,RHOC,RHOS,NE,IERR)
         IF (IP.EQ.IPMAX) GO TO 600
C
         IF (NR*2.GT.NRDIM.OR.NV*2.GT.NVDIM.OR.NS*2.GT.NSDIM) GO TO 500
C
C   DOUBLE ALL ARRAYS
C
         NR=NR*2
         NS=NS*2
         NV=NV*2
         DELR0=DELR0/2.
         DELS0=DELS0/2.
         DELETA=DELETA/2.
         IF (IWRIT.EQ.0) GO TO 300
            WRITE(6,*) 'NR,NS,NV =',NR,NS,NV
            WRITE(6,*) 'DELR,DELS,DELV =',DELR0,DELS0,A1*DELETA
            WRITE(16,*) 'NR,NS,NV =',NR,NS,NV
            WRITE(16,*) 'DELR,DELS,DELV =',DELR0,DELS0,A1*DELETA
300      CONTINUE
         DO 450 I=1,NR
	      R(I)=(2*NR*DELR0/PI)*TAN(PI*(I-0.5)/(2.*NR))
	      DELR(I)=DELR0/(COS(PI*(I-0.5)/(2.*NR)))**2
            DELV(I)=SQRT(A1**2+R(I)**2)*ETA1/FLOAT(NV)
            DO 430 J=1,NV-1
               X2M1=(R(I)/A1)**2
               XSI=DSQRT(1.D0+X2M1)
               ETA=J*ETA1/FLOAT(NV)
               Z=A1*XSI*ETA
               RP=A1*SQRT(X2M1*(1.-ETA**2))
               IF (Z.GT.(ETAT*(SQRT(RP**2/(1.-ETAT**2)+AT**2)-AT)+SEP
     &            -P(RP))) GO TO 410
C   ADD STATEMENT BELOW FOR FLAT TIP APEX
C               IF (Z.GT.(SEP-P(RP)).AND.RP.LT.10.) GO TO 410
               IF (J.EQ.1) GO TO 390
               IF (MOD(J,2).NE.0) GO TO 380
               TEMP=VAC(2,(I-1)/2+1,J/2)
               GO TO 400
380            TEMP=(VAC(2,(I-1)/2+1,J/2)+VAC(2,(I-1)/2+1,(J+1)/2))/2.
               GO TO 400
390            TEMP=(VAC(2,(I-1)/2+1,1)+VSINT(2,(I-1)/2+1))/2.
400            VAC(1,I,J)=TEMP
	         TIP(I,J)=.FALSE.
               GO TO 430
410            VAC(1,I,J)=BIAS
	         TIP(I,J)=.TRUE.
430         CONTINUE
            VAC(1,I,NV)=BIAS
	      TIP(I,NV)=.TRUE.
450      CONTINUE
         DO 460 I=1,NR
            VSINT(1,I)=VSINT(2,(I-1)/2+1)
460      CONTINUE
	   DO 463 J=1,NS-1
	      S(J)=(2*NS*DELS0/PI)*TAN(PI*J/(2.*NS))
	      DELS(J)=DELS0/(COS(PI*J/(2.*NS)))**2
463	   CONTINUE
         DO 465 I=1,NR
             SEM(1,I,1)=(SEM(2,(I-1)/2+1,1)+VSINT(2,(I-1)/2+1))/2.
465      CONTINUE
         DO 475 J=2,NS,2
            DO 470 I=1,NR
                SEM(1,I,J)=SEM(2,(I-1)/2+1,J/2)
470         CONTINUE
475      CONTINUE
         DO 485 J=3,NS,2
            DO 480 I=1,NR
                SEM(1,I,J)=(SEM(2,(I-1)/2+1,J/2)+
     &                     +SEM(2,(I-1)/2+1,(J+1)/2))/2.
480         CONTINUE
485      CONTINUE
         DO 492 J=1,NV
            DO 490 I=1,NR
               VAC(2,I,J)=VAC(1,I,J)
490         CONTINUE
492      CONTINUE
         DO 494 I=1,NR
            VSINT(2,I)=VSINT(1,I)
494      CONTINUE
         DO 498 J=1,NS
            DO 496 I=1,NR
               SEM(2,I,J)=SEM(1,I,J)
496         CONTINUE
498      CONTINUE
	   WRITE(6,*) 'LARGEST RADIUS, DEPTH =',R(NR),S(NS-1)
	   WRITE(16,*) 'LARGEST RADIUS, DEPTH =',R(NR),S(NS-1)
         GO TO 590
C
500      IF (NS*2.GT.NSDIM) GO TO 600
C
C   DOUBLE ONLY SEM ARRAY ELEMENTS
C
         NS=NS*2
         DELS0=DELS0/2.
         IF (IWRIT.EQ.0) GO TO 505
            WRITE(6,*) 'NS =',NS
            WRITE(6,*) 'DELS =',DELS0
            WRITE(16,*) 'NS =',NS
            WRITE(16,*) 'DELS =',DELS0
505      CONTINUE
	   DO 508 J=1,NS-1
	      S(J)=(2*NS*DELS0/PI)*TAN(PI*J/(2.*NS))
	      DELS(J)=DELS0/COS(PI*J/(2.*NS))
508	   CONTINUE
         DO 510 I=1,NR
             SEM(1,I,1)=(SEM(2,I,1)+VSINT(2,I))/2.
510      CONTINUE
         DO 520 J=2,NS,2
            DO 515 I=1,NR
                SEM(1,I,J)=SEM(2,I,J/2)
515         CONTINUE
520      CONTINUE
         DO 530 J=3,NS,2
            DO 525 I=1,NR
                SEM(1,I,J)=(SEM(2,I,J/2)+
     &                     +SEM(2,I,(J+1)/2))/2.
525         CONTINUE
530      CONTINUE
         DO 550 J=1,NS
            DO 540 I=1,NR
               SEM(2,I,J)=SEM(1,I,J)
540         CONTINUE
550      CONTINUE
	   WRITE(6,*) 'LARGEST DEPTH =',S(NS-1)
	   WRITE(16,*) 'LARGEST DEPTH =',S(NS-1)
590   CONTINUE
C
600   CONTINUE
      IF (IWRIT.NE.0) WRITE(6,*) 'NUMBER OF ITERATIONS = ',ITM
      IF (IWRIT.NE.0) WRITE(6,*) 'BAND BENDING AT MIDPOINT = ',PHI0
      IF (IWRIT.NE.0) WRITE(16,*) 'NUMBER OF ITERATIONS = ',ITM
      IF (IWRIT.NE.0) WRITE(16,*) 'BAND BENDING AT MIDPOINT = ',PHI0
C
      NR1=NR
      NV1=NV
      NS1=NS
      DELR1=DELR0
      DELS1=DELS0
      RETURN
      END
C
      INTEGER FUNCTION IADS(I)
      COMMON/SIZE/NR,NV,NS
      IF (I.EQ.0) GO TO 10
      IF (I.GT.NR) GO TO 20
      IADS=I
      RETURN
10    IADS=1
      RETURN
20    IADS=NR
      RETURN
      END
C
C   3-D SOLUTION FOR POTENTIAL
C
      SUBROUTINE ITER2(VAC,TIP,SEM,VSINT,R,DELR,DELV,DELS,NRDIM,NVDIM,
     &NSDIM,EP,ITMAX,PHI0,IWRIT,RHOC,RHOS,NE,IERR)
C
      DIMENSION VAC(2,NRDIM,NVDIM),SEM(2,NRDIM,NSDIM),VSINT(2,NRDIM),
     &R(NRDIM),DELR(NRDIM),DELV(NRDIM),DELS(NSDIM),RHOC(NE),
     &RHOS(NE)
	LOGICAL TIP(NRDIM,NVDIM)
      COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,DEGEN
      COMMON/SIZE/NR,NV,NS
      COMMON/STIP/BIAS,DELR0,DELS0,DELETA,A1,ESTART,DELE,EF
      DATA EEP/1.80943E-20/
C
	PHISAV=0.
      DO 500 ITER=1,ITMAX
         DO 360 I=1,NR
            X2M1=(R(I)/A1)**2
            XSI=DSQRT(1.D0+X2M1)
	      DELXSI=R(I)*DELR(I)/(XSI*A1**2)
            DO 350 J=1,NV-1
               IF (TIP(I,J)) GO TO 360
               ETA=J*DELETA
               IF (J.EQ.1) GO TO 335
                  TEMP=X2M1*(VAC(1,IADS(I+1),J)+VAC(1,IADS(I-1),J))
     &               /DELXSI**2+
     &              (1.-ETA**2)*(VAC(1,I,J+1)+VAC(1,I,J-1))/DELETA**2+
     &              (VAC(1,IADS(I+1),J)-VAC(1,IADS(I-1),J))*XSI/DELXSI-
     &              (VAC(1,I,J+1)-VAC(1,I,J-1))*ETA/DELETA
               GO TO 340
335               TEMP=X2M1*(VAC(1,IADS(I+1),1)+VAC(1,IADS(I-1),1))
     &            /DELXSI**2+
     &           (1.-ETA**2)*(VAC(1,I,2)+VSINT(1,I))/DELETA**2+
     &           (VAC(1,IADS(I+1),1)-VAC(1,IADS(I-1),1))*XSI/DELXSI-
     &           (VAC(1,I,2)-VSINT(1,I))*ETA/DELETA
340            CONTINUE
               VAC(2,I,J)=TEMP/(2.*X2M1/DELXSI**2+
     &                          2.*(1.-ETA**2)/DELETA**2)
350         CONTINUE
360      CONTINUE
         DO 370 J=1,NV
            DO 365 I=1,NR
               VAC(1,I,J)=VAC(2,I,J)
365         CONTINUE
370      CONTINUE
         DO 400 I=1,NR
	      SURFOLD=VSINT(1,I)
            ENER=EF-VSINT(1,I)
	      IF (TIP(I,3)) GO TO 375
C   3RD ORDER SOLUTION IN SEMICONDUCTOR & VACUUM
            STEMP=
     &      ((3.2*VAC(1,I,1)-1.7*VAC(1,I,2)+0.4*VAC(1,I,3))/DELV(I)+
     &      +EPSIL*(3.2*SEM(1,I,1)-1.7*SEM(1,I,2)+0.4*SEM(1,I,3))/DELS0)
            DENOM=(1.9/DELV(I)+1.9*EPSIL/DELS0)
            GO TO 385
C   3RD ORDER SOLUTION IN SEMICONDUCTOR, 2ND ORDER IN VACUUM
375	      IF (TIP(I,2)) GO TO 380
            STEMP=
     &      ((2.0*VAC(1,I,1)-0.5*VAC(1,I,2))/DELV(I)+
     &      +EPSIL*(3.2*SEM(1,I,1)-1.7*SEM(1,I,2)+0.4*SEM(1,I,3))/DELS0)
            DENOM=(1.5/DELV(I)+1.9*EPSIL/DELS0)
            GO TO 385
C   3RD ORDER SOLUTION IN SEMICONDUCTOR, 1ST ORDER IN VACUUM
380         STEMP=
     &      (VAC(1,I,1)/DELV(I)+
     &      +EPSIL*(3.2*SEM(1,I,1)-1.7*SEM(1,I,2)+0.4*SEM(1,I,3))/DELS0)
            DENOM=(1./DELV(I)+1.9*EPSIL/DELS0)
385	      CONTINUE
            IENER=NINT((ENER-ESTART)/DELE)+1
            IF (IENER.LT.1.OR.IENER.GT.NE) GO TO 388
	         RHO=RHOS(IENER)
	         GO TO 390
388	         RHO=SRHO(ENER,DELE)
	         IERR=2
390	      CONTINUE
C   ADD STATEMENT BELOW TO RESTRICT REGION OF SURFACE CHARGE
C	      IF (R(I).LT.0.) RHO=0.
            TEMP=STEMP-RHO*EEP*1.E7
            SURFNEW=TEMP/DENOM
            DELSURF=ABS(BIAS)/1.E6
            CALL GSECT3(SURFOLD,SURFNEW,DELSURF,R(I),STEMP,DENOM,RHOS,NE
     &                  ,IERR)
            VSINT(2,I)=(SURFOLD+SURFNEW)/2.
400      CONTINUE
         DO 410 I=1,NR
            VSINT(1,I)=VSINT(2,I)
410      CONTINUE
         IF (MOD(ITER,100).EQ.0) THEN
	      PHISAV2=PHISAV
            PHISAV=PHI0
            PHI0=(9.*VSINT(1,1)-VSINT(1,2))/8.
            IF (IWRIT.NE.0) WRITE(6,*) 'ITER,PHI0 =',ITER,PHI0
            IF (IWRIT.NE.0) WRITE(16,*) 'ITER,PHI0 =',ITER,PHI0
         END IF
         DO 450 J=1,NS-1
            DO 440 I=1,NR
               SEMOLD=SEM(1,I,J)
               ENER=EF-SEM(1,I,J)
               IF (J.EQ.1) GO TO 420
               STEMP=(SEM(1,IADS(I+1),J)+SEM(1,IADS(I-1),J))/DELR(I)**2+
     &         (SEM(1,I,J+1)+SEM(1,I,J-1))/DELS(J)**2+
     &         (SEM(1,IADS(I+1),J)-SEM(1,IADS(I-1),J))/(2.*R(I)*DELR(I))
	         GO TO 430
420            STEMP=(SEM(1,IADS(I+1),1)+SEM(1,IADS(I-1),1))/DELR(I)**2+
     &         (SEM(1,I,2)+VSINT(1,I))/DELS(J)**2+
     &         (SEM(1,IADS(I+1),1)-SEM(1,IADS(I-1),1))/(2.*R(I)*DELR(I))
430            CONTINUE
               IENER=NINT((ENER-ESTART)/DELE)+1
               IF (IENER.LT.1.OR.IENER.GT.NE) GO TO 433
	            RHO=RHOC(IENER)
	            GO TO 436
433	            RHO=RHOTOT(ENER,0.)
	            IERR=2
436	         CONTINUE
               TEMP=STEMP-RHO*EEP/EPSIL
	         DENOM=(2./DELR(I)**2+2./DELS(J)**2)
               SEMNEW=TEMP/DENOM
               DELSEM=ABS(BIAS)/1.E6
               CALL GSECT2(SEMOLD,SEMNEW,DELSEM,STEMP,DENOM,RHOC,NE,
     &                     IERR)
               SEM(2,I,J)=(SEMOLD+SEMNEW)/2.
440         CONTINUE
450      CONTINUE
         DO 480 J=1,NS
            DO 470 I=1,NR
               SEM(1,I,J)=SEM(2,I,J)
470         CONTINUE
480      CONTINUE   
         IF ((MOD(ITER,100).EQ.0.AND.ABS(PHI0-PHISAV).LT.EP.AND.
     &        ABS(PHISAV-PHISAV2).LT.2.*EP)) GO TO 510
500   CONTINUE
510   ITMAX=ITER
      PHI0=(9.*VSINT(1,1)-VSINT(1,2))/8.
      RETURN
      END