C   ******************** SEMICURRTIP_V4 ************************
C
C   CALLING PROGRAM FOR E-FIELD AND TUNNEL CURRENT COMPUTATIONS FOR 
C   TIP-SEMICONDUCTOR SYSTEM; TAKES INPUT FROM FILE FORT.9
C
C   VERSION 4.0 - AUG/09
C           4.1 - DEC/10 - ELIMINATE 2-BAND PARAMETERS, AND MODIFY
C                          OUTPUT OF CURRENTS AND CONDUCTANCES
C           4.2 - JAN/11 - NEW GRID FOR Z COORDINATES IN SEMICONDUCTOR
C
      PARAMETER(NRDIM=2048,NVDIM=2048,NSDIM=2048,NVDIM1=NVDIM+1,
     &NVDIM2=2048,NSDIM1=NSDIM+1,NSDIM2=20000,NEDIM=50000)
C      
      DIMENSION VAC(2,NRDIM,NVDIM),SEM(2,NRDIM,NSDIM),VSINT(2,NRDIM),
     &R(NRDIM),S(NSDIM),DELV(NRDIM),ITMAX(10),EP(10),RHOB(NEDIM),
     &RHOS(NEDIM),BBIAS(1000),NLOC(4),CR(4,NRDIM),CX(4,NRDIM),
     &CS(4,NSDIM)
      LOGICAL TIP(NRDIM,NVDIM)
      CHARACTER*1 ANS
      COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG,IINV,ISTK
      COMMON/PROTRU/RAD2
      COMMON/SURF/EN0,EN(2),DENS(2),FWHM(2),ECENT(2)
      COMMON/FD/EFDUM
      DATA EPSIL0/8.854185E-12/E/1.60210E-19/
      PI=4.*ATAN(1.)
C
C   SEMICONDUCTOR AND TIP PARAMETERS
C
      READ(9,*) NPARM
      DO 900 IPARM=1,NPARM
      READ(9,*) SLOPE
      READ(9,*) SEPIN
      READ(9,*) RAD
      READ(9,*) RAD2
      READ(9,*) DELPHI
      READ(9,*) CD
      READ(9,*) CA
      THETA=360.*ATAN(1./SLOPE)/PI
      WRITE(6,*) ' '
      WRITE(16,*) ' '
      WRITE(6,*) 'SEP, CONTACT POTENTIAL =',SEPIN,DELPHI
      WRITE(16,*) 'SEP, CONTACT POTENTIAL =',SEPIN,DELPHI
      WRITE(6,*) 'DOPING =',CD,CA
      WRITE(16,*) 'DOPING =',CD,CA
      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,*) AVBH
      read(9,*) AVBL
      read(9,*) AVBSO
      read(9,*) ESO
C   PARAMETERS BELOW ELIMINATED FROM VERSION 4.1      
C      read(9,*) E2HH
C      read(9,*) E2SO
      AVB=exp(2.*alog(sqrt(AVBH**3)+sqrt(AVBL**3))/3.)
      READ(9,*) EPSIL
      READ(9,*) TEM
      TK=TEM*8.617E-5
      READ(9,*) BMOD
      READ(9,*) IDEG
      READ(9,*) IINV
      IF ((CA.GT.CD.AND.(IINV.EQ.1.OR.IINV.EQ.3)).OR.
     &    (CD.GT.CA.AND.(IINV.EQ.2.OR.IINV.EQ.3))) THEN
         WRITE(6,*) '****** WARNING - LIKELY INCOMPATIBLE DOPING AND ',
     &   'INVERSION (IINV) PARAMETER'
         WRITE(6,*) 'CONTINUE (y/n) ?'
         READ(5,50) ANS
50       FORMAT(1A1)         
         IF (ANS.NE.'y'.AND.ANS.NE.'Y') STOP
      END IF
      READ(9,*) ISTK
      READ(9,*) DENS(1)
      READ(9,*) EN(1)
      READ(9,*) FWHM(1)
      READ(9,*) ECENT(1)
      READ(9,*) DENS(2)
      READ(9,*) EN(2)
      READ(9,*) FWHM(2)
      READ(9,*) ECENT(2)
      WRITE(6,*) 'FIRST DISTRIBUTION OF SURFACE STATES:'
      WRITE(16,*) 'FIRST DISTRIBUTION OF SURFACE STATES:'
      WRITE(6,*) 'SURFACE STATE DENSITY, EN =',DENS(1),EN(1)
      WRITE(16,*) 'SURFACE STATE DENSITY, EN =',DENS(1),EN(1)
      WRITE(6,*) 'FWHM, ECENT =',FWHM(1),ECENT(1)
      WRITE(16,*) 'FWHM, ECENT =',FWHM(1),ECENT(1)
      WRITE(6,*) 'SECOND DISTRIBUTION OF SURFACE STATES:'
      WRITE(16,*) 'SECOND DISTRIBUTION OF SURFACE STATES:'
      WRITE(6,*) 'SURFACE STATE DENSITY, EN =',DENS(2),EN(2)
      WRITE(16,*) 'SURFACE STATE DENSITY, EN =',DENS(2),EN(2)
      WRITE(6,*) 'FWHM, ECENT =',FWHM(2),ECENT(2)
      WRITE(16,*) 'FWHM, ECENT =',FWHM(2),ECENT(2)
      READ(9,*) CHI
      READ(9,*) EFTIP
C
C   GRID SIZES, STEP SIZES, AND ITERATION LIMITS
C
      READ(9,*) NRIN
      READ(9,*) NVIN
      READ(9,*) NSIN
      READ(9,*) SIZE
      READ(9,*) IPMAX
      READ(9,*) (ITMAX(IP),IP=1,IPMAX)
      READ(9,*) (EP(IP),IP=1,IPMAX)
      READ(9,*) NE
      READ(9,*) NWK
      READ(9,*) NEE
      READ(9,*) EXPANI
      READ(9,*) FRACZ
C
C    OUTPUT PARAMETER
C
      READ(9,*) IWRIT
C
C   VOLTAGES, AND SPECTRUM PARAMETERS
C
      READ(9,*) NBIAS
      READ(9,*)(BBIAS(IBIAS),IBIAS=1,NBIAS)
      READ(9,*) ANEG
      READ(9,*) APOS
      READ(9,*) VSTART
      if (vstart.lt.0.) then
         dels1=abs(vstart)*aneg
      else
         dels1=abs(vstart)*apos
      end if
c   sep0 is separation at V=0
      sep0=sepin-dels1
C
C   PARAMETERS FOR CONTOUR PLOT
C
      READ(9,*) NUMC
      READ(9,*) DELP
C
C   FIND FERMI-LEVEL POSITION
C
      CALL EFFIND(EF)
      WRITE(6,*) 'FERMI-LEVEL =',EF
      WRITE(16,*) 'FERMI-LEVEL =',EF
      RHOCC=RHOC(EF,0.)
      RHOVV=RHOV(EF,0.)
      WRITE(6,*) 'CARRIER DENSITY IN CB, VB =',RHOCC,RHOVV
      WRITE(16,*)'CARRIER DENSITY IN CB, VB =',RHOCC,RHOVV
C
C   FIND OVERALL CHARGE NEUTRALITY LEVEL, EN0
C
      IF (DENS(1).EQ.0.AND.DENS(2).EQ.0.) THEN
         EN0=0.
      ELSE
         IF (DENS(1).EQ.0.) THEN
            EN0=EN(2)
         ELSE
            IF (DENS(2).EQ.0.) THEN
               EN0=EN(1)
            ELSE
               CALL ENFIND(EN(1),EN(2),EN0)
            END IF
         END IF
      END IF
      WRITE(6,*) 'CHARGE-NEUTRALITY LEVEL =',EN0
      PHIS=EF-EN0
C
C   ****************** LOOP OVER BIAS VOLTAGES ****************************
C
      DO 600 IBIAS=1,NBIAS
      write(6,*) ' '
      write(16,*) ' '
      BIAS0=BBIAS(IBIAS)
      if (bias0.le.0) then
         sep=sep0+aneg*abs(bias0)
      else
         sep=sep0+apos*abs(bias0)
      end if
      write(6,*) 'SEPARATION =',sep
      write(16,*) 'SEPARATION =',sep
C
      IMODMAX=1
      IF (BMOD.EQ.0.) IMODMAX=-1
      do 550 imod=-1,IMODMAX,2
      BIASSAV=BIAS
      bias=bias0+imod*bmod*sqrt(2.)
      PHITIP=BIAS+DELPHI
      write(6,*) ' '
      write(16,*) ' '
      WRITE(6,*) 'BIAS, TIP POTENTIAL =',BIAS,PHITIP
      WRITE(16,*) 'BIAS, TIP POTENTIAL =',BIAS,PHITIP
C
C   1-D SOLUTION FOR BAND BENDING
C
      IF ((CD-CA).EQ.0.) GO TO 105
      W=1.E9*SQRT(2.*EPSIL*EPSIL0*AMAX1(1.,ABS(PHITIP))/
     &   (ABS(CD-CA)*1.E6*E))
      WRITE(6,*) '1-D ESTIMATE OF DEPLETION WIDTH (NM) =',W
      WRITE(16,*) '1-D ESTIMATE OF DEPLETION WIDTH (NM) =',W
      GO TO 110
105   W=1.E10
C
C   CONSTRUCT TABLES OF CHARGE DENSITY VALUES
C
110   ESTART=AMIN1(EF,EF-PHITIP,EF-PHIS)
      EEND=AMAX1(EF,EF-PHITIP,EF-PHIS)
      ETMP=EEND-ESTART
      ESTART=ESTART-2.*ETMP
      EEND=EEND+2.*ETMP
      DELE=(EEND-ESTART)/FLOAT(NE-1)
C   PLACE ONE OF THE TABLE VALUES FOR ENERGY AT EF +/- (DELE/2)
      NETMP=NINT((EF-ESTART)/DELE)
      ESTART=EF-(NETMP-0.5)*DELE
      EEND=ESTART+(NE-1)*DELE
      WRITE(6,*) 'ESTART,EEND,NE =',ESTART,EEND,NE
      WRITE(16,*) 'ESTART,EEND,NE =',ESTART,EEND,NE
      WRITE(6,*) 'COMPUTING TABLE OF BULK CHARGE DENSITIES'
      WRITE(16,*) 'COMPUTING TABLE OF BULK CHARGE DENSITIES'
      CALL SEMIRHO(DELE,ESTART,NE,NEDIM,RHOB)
      WRITE(6,*) 'COMPUTING TABLE OF SURFACE CHARGE DENSITIES'
      WRITE(16,*) 'COMPUTING TABLE OF SURFACE CHARGE DENSITIES'
      IF (DENS(1).EQ.0.AND.DENS(2).EQ.0.) THEN
         DO 115 IE=1,NE
            RHOS(IE)=0.
115      CONTINUE
      ELSE
         CALL SURFRHO(DELE,ESTART,NE,NEDIM,RHOS)
      END IF
C
C   SEMICONDUCTOR GRID SIZE AND SPACING
C
      NR=NRIN
      DELR=RAD*SIZE
      IF (RAD2.NE.0.) DELR=AMIN1(RAD2,DELR)
      DELR=AMIN1(DELR,W/NR)
      NV=NVIN
      NS=NSIN
      DELS=RAD*SIZE
      IF (RAD2.NE.0.) DELS=AMIN1(RAD2,DELS)
      DELS=AMIN1(DELS,W/NS)
C
C   SOLVE THE POTENTIAL PROBLEM
C
      IWRIT1=IWRIT
      IF (IWRIT.GT.5) IWRIT1=MOD(IWRIT,5)
      CALL SEMITIP2(SEP,RAD,SLOPE,ETAT,A,Z0,C,VAC,TIP,SEM,VSINT,
     &R,S,DELV,DELR,DELS,NRDIM,NVDIM,NSDIM,NR,NV,NS,PHITIP,IWRIT1,
     &NE,ESTART,DELE,RHOB,RHOS,EF,ITMAX,EP,IPMAX,PHI0,IERR,CR,CX,CS)
      WRITE(6,120) NR,NS,NV,IERR
      WRITE(16,120) NR,NS,NV,IERR
120   FORMAT(' RETURN FROM SEMTIP2, NR,NS,NV,IERR =',4I7)
C
C   CURRENT COMPUTATION
C
      WRITE(10,430) RAD,SEP,BIAS,DELPHI,PHI0
430   FORMAT(' ',5G12.4)
      WRITE(6,*) ' '
      WRITE(16,*)' '
      WRITE(6,*) 'COMPUTATION OF CURRENT:'
      WRITE(16,*)'COMPUTATION OF CURRENT:'
      NSP=NINT(FRACZ*NS)
      SDEPTH=(2*NS*DELS/PI)*TAN(PI*NSP/(2.*NS))
      WRITE(6,*) '# GRID POINTS INTO SEMICONDUCTOR USED FOR INTEGRATION'
     &,' =',NSP
      WRITE(16,*)'# GRID POINTS INTO SEMICONDUCTOR USED FOR INTEGRATION'
     &,' =',NSP
      WRITE(6,*) 'DEPTH INTO SEMICONDUCTOR USED FOR INTEGRATION =',
     &SDEPTH
      WRITE(16,*)'DEPTH INTO SEMICONDUCTOR USED FOR INTEGRATION =',
     &SDEPTH
      IWRIT1=MOD(IWRIT,5)
      CALL CURR2(0,VAC,TIP,SEM,VSINT,NV,NR,NS,NVDIM,NRDIM,NSDIM,
     &S,SEP,DELV,BIAS,EF,CHI,EFTIP,DELPHI,AVBH,AVBL,AVBSO,
     &ESO,E2HH,E2SO,nee,nwk,PHI0,NVDIM1,NVDIM2,NSDIM1,NSDIM2,EXPANI,
     &FRACZ,NLOC,CURRVE,CURRVL,CURRCE,CURRCL,CURRE,CURRL,IWRIT1)
      CURR=CURRE+CURRL
      CURRV=CURRVE+CURRVL
      CURRC=CURRCE+CURRCL
C
C   OUTPUT CURRENT AND CONDUCTANCE
C
      if (imod.eq.-1) then
         CSAV=CURR
         CSAVE=CURRE
         CSAVL=CURRL
         CSAVV=CURRV
         CSAVVE=CURRVE
         CSAVVL=CURRVL
         CSAVC=CURRC
         CSAVCE=CURRCE
         CSAVCL=CURRCL
      else
         COND=(CURR-CSAV)
         CONDE=(CURRE-CSAVE)
         CONDL=(CURRL-CSAVL)
         CONDV=(CURRV-CSAVV)
         CONDVE=(CURRVE-CSAVVE)
         CONDVL=(CURRVL-CSAVVL)
         CONDC=(CURRC-CSAVC)
         CONDCE=(CURRCE-CSAVCE)
         CONDCL=(CURRCL-CSAVCL)
         COND=COND*exp(2.*10.*(sep-sep0-dels1))/(2.*bmod*sqrt(2.))
         CONDE=CONDE*exp(2.*10.*(sep-sep0-dels1))/(2.*bmod*sqrt(2.))
         CONDL=CONDL*exp(2.*10.*(sep-sep0-dels1))/(2.*bmod*sqrt(2.))
         CONDV=CONDV*exp(2.*10.*(sep-sep0-dels1))/(2.*bmod*sqrt(2.))
         CONDVE=CONDVE*exp(2.*10.*(sep-sep0-dels1))/(2.*bmod*sqrt(2.))
         CONDVL=CONDVL*exp(2.*10.*(sep-sep0-dels1))/(2.*bmod*sqrt(2.))
         CONDC=CONDC*exp(2.*10.*(sep-sep0-dels1))/(2.*bmod*sqrt(2.))
         CONDCE=CONDCE*exp(2.*10.*(sep-sep0-dels1))/(2.*bmod*sqrt(2.))
         CONDCL=CONDCL*exp(2.*10.*(sep-sep0-dels1))/(2.*bmod*sqrt(2.))
         write(15,*) bias0,COND,CONDE,CONDL
         write(93,*) bias0,CONDV,CONDVE,CONDVL
         write(94,*) bias0,CONDC,CONDCE,CONDCL
      end if
      CURR=CURR*exp(2.*10.*(sep-sep0-dels1))
      CURRE=CURRE*exp(2.*10.*(sep-sep0-dels1))
      CURRL=CURRL*exp(2.*10.*(sep-sep0-dels1))
      CURRV=CURRV*exp(2.*10.*(sep-sep0-dels1))
      CURRVE=CURRVE*exp(2.*10.*(sep-sep0-dels1))
      CURRVL=CURRVL*exp(2.*10.*(sep-sep0-dels1))
      CURRC=CURRC*exp(2.*10.*(sep-sep0-dels1))
      CURRCE=CURRCE*exp(2.*10.*(sep-sep0-dels1))
      CURRCL=CURRCL*exp(2.*10.*(sep-sep0-dels1))
      write(14,*) bias,CURR,CURRE,CURRL
      write(91,*) bias,CURRV,CURRVE,CURRVL
      write(92,*) bias,CURRC,CURRCE,CURRCL
C
C   PLOT CROSS-SECTIONAL PROFILE
C
      IF (IWRIT.GE.1) THEN
         DO 505 J=NV,1,-1
           WRITE(11,*) -J*DELV(1),VAC(1,1,J),VAC(1,NR,J)
505      CONTINUE
         WRITE(11,*) 0.,VSINT(1,1),VSINT(1,NR)
         DO 510 J=1,NS
            WRITE(11,*) S(J),SEM(1,1,J),SEM(1,NR,J)
510      CONTINUE
C
C   PLOT SURFACE POTENTIAL AND SURFACE CHARGE DENSITY
C
         WRITE(12,*) -R(1),VSINT(1,1)
         DO 520 I=1,NR
            WRITE(12,*) R(I),VSINT(1,I)
520      CONTINUE
         CLOSE(11)
         CLOSE(12)
      END IF
C
550   continue
600   CONTINUE
C
C   PLOT CONTOURS
C
      IF (IWRIT.GE.2) THEN
         CALL CONTR(ETAT,VAC,TIP,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=36+4*nr*(nv+ns+3)+4*ns
         open(unit=13,file='fort.13',access='direct',recl=nrecl)                                            
         write(13,rec=1) nr,nv,ns,sep,rad,rad2,slope,bias,epsil,
     &   ((vac(1,i,j),i=1,nr),j=1,nv),
     &   ((sem(1,i,j),i=1,nr),j=1,ns),(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,*)
C
      stop
      end
C
      real function p(r)
      common/protru/rad2
      p=0.
      if (r.lt.rad2) p=sqrt(rad2**2-r**2)
      return
      end
c
      real*8 function sig(ID,EF,PHI)
      COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG
      common/SURF/EN0,EN(2),DENS(2),FWHM(2),ECENT(2)
      PI=4.*ATAN(1.)
c
c     sig=0.
c     if ((ef-phi).lt.0.) return
c     if ((ef-phi).gt.egap) return
c
      if (fwhm(ID).eq.0.) go to 200
      width=fwhm(ID)/(2.*sqrt(2.*alog(2.)))
      sig=-dexp(-1.d0*(EF-PHI-(EN(ID)+ecent(ID)))**2/(2.*width**2))+
     &     dexp(-1.d0*(EF-PHI-(EN(ID)-ecent(ID)))**2/(2.*width**2))
      sig=sig*dens(ID)/(SQRT(2.*pi)*width)
      return
c     
c     sig=exp(-(EF-PHI-EN)**2/(2.*width**2))
c     sig=sig*dens/(SQRT(2.*pi)*width)
c
c     sig=exp(-(EF-PHI)/width)-
c     &    exp((EF-PHI-egap)/width)
c     sig=sig*dens/WIDTH
c
200   sig=dens(ID)
      if ((ef-phi).gt.en(ID)) sig=-sig
      return
      end