C   ******************** SEMITIP_V1 ************************
C
C   CALLING PROGRAM FOR E-FIELD COMPUTATIONS FOR TIP-SEMICONDUCTOR
C   TAKES INPUT FROM FILE FORT.9
C
      DIMENSION VAC(2,1024,4096),SEM(2,1024,2048),VSINT(2,2048),
     &DELV(2048),ITMAX(10),EP(10),RHOC(50000),R(2048),S(2048)
      COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG
      COMMON/PROTRU/RAD2
      DATA PI/3.141592654/EPSIL0/8.854185E-12/E/1.60210E-19/
C
      NRDIM=1024
      NVDIM=4096
      NSDIM=2048
      NE=50000
C
C   SEMICONDUCTOR AND TIP PARAMETERS
C
      READ(9,*) NPARM
      DO 900 IPARM=1,NPARM
      READ(9,*) SLOPE
      READ(9,*) SEP
      READ(9,*) RAD
      READ(9,*) RAD2
      READ(9,*) BIAS
      READ(9,*) CD
      READ(9,*) CA
      THETA=360.*ATAN(1./SLOPE)/PI
      WRITE(6,*) ' '
      WRITE(16,*) ' '
      WRITE(6,*) 'SEP, BIAS, DOPING =',SEP,BIAS,CD
      WRITE(16,*) 'SEP, BIAS, DOPING =',SEP,BIAS,CD
      WRITE(6,*) 'RAD, SLOPE, ANGLE =',RAD,SLOPE,THETA 
      WRITE(16,*) 'RAD, SLOPE, ANGLE =',RAD,SLOPE,THETA 
      READ(9,*) EGAP
      READ(9,*) ED
      READ(9,*) EA
      READ(9,*) ACB
      READ(9,*) AVB
      READ(9,*) EPSIL
      READ(9,*) TEM
      TK=TEM*8.617E-5
      READ(9,*) IDEG
C
C   FIND FERMI-LEVEL POSITION
C
      CALL EFFIND(EF)
      WRITE(6,*) 'DOPING, FERMI-LEVEL =',CD,CA,EF
      WRITE(16,*) 'DOPING, FERMI-LEVEL =',CD,CA,EF
C
C   1-D SOLUTION FOR BAND BENDING
C
      W=1.E9*SQRT(2.*EPSIL*EPSIL0*ABS(BIAS)/(ABS(CD-CA)*1.E6*E))
      WRITE(6,*) 'ABSOLUTE MAX DEPLETION WIDTH (nm) =',W
      WRITE(16,*) 'ABSOLUTE MAX DEPLETION WIDTH (nm) =',W
      ATMP=1.
      BTMP=2.*EPSIL*SEP
      CTMP=-1.E18*2.*ABS(BIAS)*EPSIL*EPSIL0/(ABS(CD-CA)*1.E6*E)
      W=(-BTMP+SQRT(BTMP**2-4.*ATMP*CTMP))/(2.*ATMP)
      WRITE(6,*) '1-D SOLN FOR DEPLETION WIDTH (nm) =',W
      WRITE(16,*) '1-D SOLN FOR DEPLETION WIDTH (nm) =',W
      PHI00=(W*1.E-9)**2*E*ABS(CD-CA)*1.E6/(2.*EPSIL*EPSIL0)
     &     *(BIAS/ABS(BIAS))
      WRITE(6,*) '1-D SOLN FOR BAND BENDING (V) =',PHI00
      WRITE(16,*) '1-D SOLN FOR BAND BENDING (V) =',PHI00
C
C   CONSTRUCT TABLE OF CHARGE DENSITY VALUES
C
      ESTART=EF
      IF ((EF-PHI00).LT.ESTART) ESTART=EF-PHI00
      ESTART=ESTART-0.1*ABS(PHI00)
      DELE=(1.2*ABS(PHI00))/FLOAT(NE)
      WRITE(6,*) 'ESTART,DELE,NE =',ESTART,DELE,NE
      WRITE(16,*) 'ESTART,DELE,NE =',ESTART,DELE,NE
      CALL SEMIRHO(DELE,ESTART,NE,RHOC)
C
C   SEMICONDUCTOR GRID SIZE AND SPACING
C
      READ(9,*) NR
      READ(9,*) NV
      NV=MAX0(NV,NINT(2.*RAD*SLOPE**2/SEP))
      IF (NV.GT.NVDIM) GO TO 900
      READ(9,*) NS
      READ(9,*) ISIZE
      DELR=ISIZE*W/FLOAT(NR)
      DELS=ISIZE*W/FLOAT(NS)
C
C   ITERATION LIMITS
C
      READ(9,*) IPMAX
      READ(9,*) (ITMAX(IP),IP=1,IPMAX)
      READ(9,*) (EP(IP),IP=1,IPMAX)
      READ(9,*) IWRIT
C
C   SOLVE THE PROBLEM !
C
      CALL SEMITIP2(SEP,RAD,SLOPE,ETAT,AT,ETA1,A1,VAC,SEM,VSINT,
     &DELV,DELR,DELS,NRDIM,NVDIM,NSDIM,NR,NV,NS,BIAS,IWRIT,
     &NE,ESTART,DELE,RHOC,EF,ITMAX,EP,IPMAX,PHI0,IERR)
      WRITE(6,500) NR,NS,NV,IERR
      WRITE(16,500) NR,NS,NV,IERR
500   FORMAT(' RETURN FROM SEMTIP2, NR,NS,NV,IERR =',4I7)
C
C   OUTPUT
C
      WRITE(10,510) RAD,SEP,BIAS,THETA,PHI0,PHI00
510   FORMAT(' ',6G12.4)
      IF (IWRIT.GE.1) THEN
      DO 610 I=1,NR
         R(I)=(I-0.5)*DELR
610   CONTINUE
      DO 620 J=1,NS
         S(J)=J*DELS
620   CONTINUE
C
C   PLOT CROSS-SECTIONAL PROFILE
C
      DO 705 J=NV,1,-1
        WRITE(11,*) -J*DELV(1),VAC(1,1,J)
705   CONTINUE
      WRITE(11,*) 0.,VSINT(1,1)
      DO 710 J=1,NS-1
          WRITE(11,*) S(J),SEM(1,1,J)
710   CONTINUE

C
C   PLOT SURFACE POTENTIAL
C
      WRITE(12,*) -R(1),VSINT(1,1)
      DO 720 I=1,NR-1
          WRITE(12,*) R(I),VSINT(1,I)
720   CONTINUE
      END IF
C
C   PLOT CONTOURS
C
      IF (IWRIT.GE.2) THEN
      READ(9,*) NUMC
      READ(9,*) DELP
      CALL CONTR(ETA1,VAC,BIAS,SEM,VSINT,R,S,DELV,NRDIM,NVDIM,
     &NSDIM,NR,NV,NS,NUMC,DELP)
      END IF
C
C   OUTPUT ENTIRE POTENTIAL
C
      IF (IWRIT.GE.3) THEN
      nrecl=5+2*nr*(nv+ns+3)+2*ns
      open(unit=13,file='fort.13',access='direct',recl=nrecl)                                            
      write(13,rec=1) nr,nv,ns,eta1,((vac(1,i,j),i=1,nr),j=1,nv),
     &((sem(1,i,j),i=1,nr),j=1,nv),(vsint(1,i),i=1,nr),
     &(r(i),i=1,nr),(s(j),j=1,ns),(delv(i),i=1,nr)
      END IF
C
900   CONTINUE
      WRITE(6,*) 'PRESS THE ENTER KEY TO EXIT'
      WRITE(16,*) 'PRESS THE ENTER KEY TO EXIT'
      READ(5,*)
      STOP
      END
C
      REAL FUNCTION P(R)
      COMMON/PROTRU/RAD2
      P=0.
      IF (R.LT.RAD2) P=SQRT(RAD2**2-R**2)
      RETURN
      END