C   ******************** SEMITIP_V5 ************************
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 5.0 - AUG/09
C           5.1 - DEC/10 - ELIMINATE 2-BAND PARAMETERS, AND MODIFY
C                          OUTPUT OF CURRENTS AND CONDUCTANCES
C           5.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),CDESEM(NSDIM),CDLSEM(NSDIM),
     &CDEVAC(NVDIM1),CDLVAC(NVDIM1),CDSEM(NRDIM,NSDIM),CDSURF(NRDIM),
     &CDSEM1(NRDIM,NSDIM),CDSURF1(NRDIM),RHOVB(NEDIM),RHOCB(NEDIM),
     &plot1(NRDIM,NSDIM),plot2(NRDIM,NSDIM),plot3(NRDIM,NSDIM),
     &img(512,512)
      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,*) FRACZ1
      READ(9,*) FRACZ2
      READ(9,*) FRACR
      IF (FRACZ2.GE.FRACZ1) THEN
         WRITE(6,*) 'ERROR - FRACZ2 MUST BE LESS THAN FRACZ1'
         WRITE(16,*)'ERROR - FRACZ2 MUST BE LESS THAN FRACZ1'
         WRITE(6,*) 'PRESS THE ENTER KEY TO CONTINUE'
         WRITE(16,*) 'PRESS THE ENTER KEY TO CONTINUE'
         STOP
      END IF
      READ(9,*) NCD
      READ(9,*) EPSC
      READ(9,*) FILT
      NFILT=MAX0(NINT(FILT+0.5),1)
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,RHOCB,RHOVB)
      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 112 IE=1,NE
            RHOS(IE)=0.
112      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 INITIAL POTENTIAL PROBLEM
C
      IQUAN=0
      IINIT=1
      DO 117 I=1,NR*2**(IPMAX-1)
         CDSURF(I)=0.
         DO 115 J=1,NS*2**(IPMAX-1)
            CDSEM(I,J)=0.
115       CONTINUE
117   CONTINUE
      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,RHOVB,RHOCB,RHOS,EF,ITMAX,EP,IPMAX,PHI0,IERR,
     &CDSEM,CDSURF,FRACZ2,FRACR,IINIT,IQUAN)
      WRITE(6,120) NR,NS,NV,IERR
      WRITE(16,120) NR,NS,NV,IERR
120   FORMAT(' RETURN FROM SEMTIP2, NR,NS,NV,IERR =',4I7)
      IF (NCD.EQ.0) GO TO 420
C
C   INITIALIZE THE FILTERED CHARGE DENSITY ARRAY; ZERO VALUES ARE
C   EMPLOYED FOR INITCD=0, OR SEMI-CLASSICAL VALUES FOR INITCD=1
C
      INITCD=1
      NSP=NINT(FRACZ1*NS)
      NRP=NINT(FRACR*NR)
      IF (FILT.NE.0.) THEN
         DO 160 I=1,NRP
            IF (INITCD.EQ.0) THEN
               CDSURF1(I)=0.
            ELSE
               ENER=EF-VSINT(1,I)
               IENER=NINT((ENER-ESTART)/DELE)+1
               IF (IENER.GE.1.AND.IENER.LE.NE) THEN
                  RHO=-RHOCB(IENER)+RHOVB(IENER)
               ELSE
                  RHO=-RHOC(ENER,0.)+RHOV(ENER,0.)
               END IF
               CDSURF1(I)=RHO*1.E-21*(DELS/2.)/(1-EXP(-1./FILT))
            END IF
            DO 150 J=1,NSP
               IF (INITCD.EQ.0) THEN
                  CDSEM1(I,J)=0.
               ELSE
                  ENER=EF-SEM(1,I,J)
                  IENER=NINT((ENER-ESTART)/DELE)+1
                  IF (IENER.GE.1.AND.IENER.LE.NE) THEN
                     RHO=-RHOCB(IENER)+RHOVB(IENER)
                  ELSE
                     RHO=-RHOC(ENER,0.)+RHOV(ENER,0.)
                  END IF
                  CDSEM1(I,J)=RHO*1.E-21/(1-EXP(-1./FILT))
               END IF
150         CONTINUE
160      CONTINUE
      END IF
C
      IQUAN=1
      IINIT=0
      PHISAV1=0.
      PHISAV2=0.
      NSP=NINT(FRACZ1*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
      NSP=NINT(FRACZ2*NS)
      SDEPTH=(2*NS*DELS/PI)*TAN(PI*NSP/(2.*NS))
      WRITE(6,*) '# GRID POINTS INTO SEMICONDUCTOR USED FOR QUANTUM',
     &' CHARGE =',NSP
      WRITE(16,*)'# GRID POINTS INTO SEMICONDUCTOR USED FOR QUANTUM',
     &' CHARGE =',NSP
      WRITE(6,*) 'DEPTH INTO SEMICONDUCTOR USED FOR QUANTUM CHARGE =',
     &SDEPTH
      WRITE(16,*)'DEPTH INTO SEMICONDUCTOR USED FOR QUANTUM CHARGE =',
     &SDEPTH
      NRP=NINT(FRACR*NR)
      RDIST=(2*NR*DELR/PI)*TAN(PI*(NRP-0.5)/(2.*NR))
      WRITE(6,*) '# RADIAL GRID POINTS USED FOR QUANTUM CHARGE =',
     &NRP
      WRITE(16,*)'# RADIAL GRID POINTS USED FOR QUANTUM CHARGE =',
     &NRP
      WRITE(6,*) 'RADIAL DISTANCE USED FOR QUANTUM CHARGE =',
     &RDIST
      WRITE(16,*)'RADIAL DISTANCE USED FOR QUANTUM CHARGE =',
     &RDIST
      WRITE(6,170)
      WRITE(16,170)
170   FORMAT(/'*** SELF-CONSISTENCY LOOP *** '/4x,
     &'BIAS       ITER      PHI0       NLOC    CD SURF        CD BULK')
C
C   GET CHARGE DENSITY (BY EVALUATING STATES AT SINGLE K-PARALLEL POINT)
C   AND RE-SOLVE FOR THE POTENTIAL, ETC.
C
      NWK1=1
      DO 400 ICD=1,NCD
         DO 245 I=0,NRP
            IWRIT1=0
            IF (I.EQ.0.AND.IWRIT.GT.5) IWRIT1=MOD(IWRIT,5)
            CALL CURR2(I,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,NWK1,PHI0,NVDIM1,NVDIM2,NSDIM1,NSDIM2,
     &      EXPANI,FRACZ1,NLOC,CDESEM,CDESURF,CDLSEM,CDLSURF,CDEVAC,
     &      CDLVAC,CURRVE,CURRVL,CURRCE,CURRCL,CURRE,CURRL,IWRIT1)
            IF (I.EQ.0) THEN
               NSUM=NLOC(1)+NLOC(2)+NLOC(3)+NLOC(4)
               CDS0=CDESURF+CDLSURF
               CDB0=CDESEM(1)+CDLSEM(1)
            ELSE
               CDSURF(I)=CDESURF+CDLSURF
               DO 230 J=1,NSP
                  CDSEM(I,J)=CDESEM(J)+CDLSEM(J)
                  plot1(I,J)=abs(CDLSEM(J))
                  plot2(I,J)=abs(CDESEM(J))
                  plot3(I,J)=abs(CDSEM(I,J))
230            CONTINUE            
            END IF
            IF (IWRIT.GE.7) THEN
               IF (I.EQ.1) THEN
                  DO 235 J=NV+1,1,-1
                     WRITE(17,*) -SEP*(J-1)/(NV+1),CDEVAC(J)+CDLVAC(J),
     &                     CDEVAC(J),CDLVAC(J)
235               CONTINUE
                  DO 240 J=1,NSP
                     WRITE(17,*) S(J),CDESEM(J)+CDLSEM(J),CDESEM(J),
     &                     CDLSEM(J)
240               CONTINUE
               END IF
               IF (I.EQ.0) THEN
                  WRITE(18,*) 0.,CDEVAC(1)+CDLVAC(1),CDEVAC(1),CDLVAC(1)
               ELSE
                  WRITE(18,*) R(I),CDEVAC(1)+CDLVAC(1),CDEVAC(1),
     &                  CDLVAC(1)
               END IF
            END IF
245      CONTINUE
         IF (IWRIT.GE.7) THEN 
            CLOSE(17)
            CLOSE(18)
         END IF
C
C   FILTER THE CHARGE DENSITY
C
         IF (FILT.NE.0) THEN
            DO 330 I=1,NRP
               CDSURF1(I)=CDSURF(I)+EXP(-1./FILT)*CDSURF1(I)
               DO 320 J=1,NSP
                  CDSEM1(I,J)=CDSEM(I,J)+EXP(-1./FILT)*CDSEM1(I,J)
320            CONTINUE
330         CONTINUE
C     
            DO 360 I=1,NRP
               CDSURF(I)=CDSURF1(I)*(1-EXP(-1./FILT))
               DO 350 J=1,NSP
                  CDSEM(I,J)=CDSEM1(I,J)*(1-EXP(-1./FILT))
                  plot3(I,J)=abs(CDSEM(I,J))
350            CONTINUE
360         CONTINUE
         END IF
C
C   OUTPUT GRAY-SCALE IMAGES OF CHARGE DENSITIES
C
         IF (IWRIT.GE.7) THEN 
            ilog=1
            iexpan=3
            shade=0.
            igraycut=1
            call PlotGray(plot1,nrdim,nsdim,nrp,nsp,ilog,iexpan,shade,
     &      ibin,ixbin,iybin,px,py,pz,ipow,ycorrec,nfilt,ixfilt,iyfilt,
     &      irfilt,disc,nscal,zsmin,zsmax,igraycut,grang,smin,smax,img)
            call dispir(1,img)
            call PlotGray(plot2,nrdim,nsdim,nrp,nsp,ilog,iexpan,shade,
     &      ibin,ixbin,iybin,px,py,pz,ipow,ycorrec,nfilt,ixfilt,iyfilt,
     &      irfilt,disc,nscal,zsmin,zsmax,igraycut,grang,smin,smax,img)
            call dispir(2,img)
            call PlotGray(plot3,nrdim,nsdim,nrp,nsp,ilog,iexpan,shade,
     &      ibin,ixbin,iybin,px,py,pz,ipow,ycorrec,nfilt,ixfilt,iyfilt,
     &      irfilt,disc,nscal,zsmin,zsmax,igraycut,grang,smin,smax,img)
            call dispir(3,img)
         END IF
C
C   RE-SOLVE FOR POTENTIAL & STATES
C
         IPMAX2=1
         ITMAX2=1
         IWRIT1=0
         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,RHOVB,RHOCB,RHOS,EF,ITMAX2,EP,IPMAX2,PHI0,
     &   IERR,CDSEM,CDSURF,FRACZ2,FRACR,IINIT,IQUAN)
         WRITE(6,365) BIAS,ICD,PHI0,NSUM,CDS0*1.E14,CDB0*1.E21
         WRITE(16,365)BIAS,ICD,PHI0,NSUM,CDS0*1.E14,CDB0*1.E21
365      FORMAT(G13.5,I5,3X,G13.6,I4,2X,2G14.5)
         IF (ICD.EQ.NCD) GO TO 400
         IF (ICD.GE.NFILT.AND.(ABS(PHI0-PHISAV1).LT.EPSC).AND.
     &       (ABS(PHISAV1-PHISAV2).LT.2.*EPSC)) GO TO 420
         PHISAV2=PHISAV1
         PHISAV1=PHI0
C
400   CONTINUE
C
C   FINAL CURRENT COMPUTATION
C
420   WRITE(10,430) RAD,SEP,BIAS,DELPHI,PHI0
430   FORMAT(' ',5G12.4)
      WRITE(6,*) ' '
      WRITE(16,*)' '
      WRITE(6,*) 'FINAL COMPUTATION OF CURRENT:'
      WRITE(16,*)'FINAL COMPUTATION OF CURRENT:'
      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,FRACZ1,NLOC,CDESEM,CDESURF,CDLSEM,CDLSURF,CDEVAC,CDLVAC,
     &CURRVE,CURRVL,CURRCE,CURRCL,CURRE,CURRL,IWRIT1)
      IF (IWRIT.GE.2)
     &   WRITE(19,*) BIAS,CDEVAC(1)+CDLVAC(1),CDEVAC(1),CDLVAC(1)
      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
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
C   OUTPUT CHARGE DENSITIES
C
      IF (IWRIT.GE.2) THEN
         DO 535 J=NV+1,1,-1
            WRITE(17,*) -SEP*(J-1)/(NV+1),CDEVAC(J)+CDLVAC(J),
     &            CDEVAC(J),CDLVAC(J)
535      CONTINUE
         DO 540 J=1,NSP
            WRITE(17,*) S(J),CDESEM(J)+CDLSEM(J),CDESEM(J),
     &            CDLSEM(J)
540      CONTINUE
         WRITE(18,*) 0.,CDEVAC(1)+CDLVAC(1),CDEVAC(1),CDLVAC(1)
         CLOSE(17)
         CLOSE(18)
         ilog=1
         iexpan=3
         shade=0.
         igraycut=1
         callPlotGray(plot1,nrdim,nsdim,nrp,nsp,ilog,iexpan,shade,
     &   ibin,ixbin,iybin,px,py,pz,ipow,ycorrec,nfilt,ixfilt,iyfilt,
     &   irfilt,disc,nscal,zsmin,zsmax,igraycut,grang,smin,smax,img)
         call dispir(1,img)
         callPlotGray(plot2,nrdim,nsdim,nrp,nsp,ilog,iexpan,shade,
     &   ibin,ixbin,iybin,px,py,pz,ipow,ycorrec,nfilt,ixfilt,iyfilt,
     &   irfilt,disc,nscal,zsmin,zsmax,igraycut,grang,smin,smax,img)
         call dispir(2,img)
         callPlotGray(plot3,nrdim,nsdim,nrp,nsp,ilog,iexpan,shade,
     &   ibin,ixbin,iybin,px,py,pz,ipow,ycorrec,nfilt,ixfilt,iyfilt,
     &   irfilt,disc,nscal,zsmin,zsmax,igraycut,grang,smin,smax,img)
         call dispir(3,img)
      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,IINV,ISTK
      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