c
c     Copyright (c) 1999,2000, Compaq Computer Corporation
c
c     See the file "license.terms" for information on usage and redistribution
c     of this file, and for a DISCLAIMER OF ALL WARRANTIES.
c
      subroutine  e_jcomp(psi,e_conc,h_conc,dop_conc,
     1                    el,et,ez,ej,e_quasi_fermi,h_quasi_fermi,
     2                    e_lat_cur,e_trans_cur,h_lat_cur,h_trans_cur,
     3                    e_mob,h_mob,
     4                    e_charge,h_charge,net_charge,
     5                    e_lat_vel,e_trans_vel,e_vel_mag,
     6                    h_lat_vel,h_trans_vel,h_vel_mag,
     7                    in_binary_file, nintr, ut, temperature, tox,
     8                    vs, vb, xl, model, three_d,
     9                    nonplanar_boundary, chtyp)
c
      implicit none
c
c     patterned after genrat - estimate generation rate.
c        by a.schuetz   1981   tu vienna
c
c
c     3-Jun-85  J. Faricelli
c        stolen for use in postmini to compute electric field,
c        interpolate j from mid-points of mesh to mesh points and
c        calculate E in direction of j at mesh points.
c     10-Jun-85  J. Faricelli
c        array indices changed to use full matrix...
c     31-Dec-85  J. Faricelli
c        compute quasi-Fermi levels
c     03-Feb-87  J. Faricelli
c        compute majority carriers when not explicitly solved
c     04-Feb-87  J. Faricelli
c        compute various charge arrays
c     27_may-87  J. Faricelli
c        convert arrays from min/maj to electron/hole
c     16-Oct-87  J. Faricelli
c        set potential inside gate electrode
c     29-Sep-88  J. Faricelli
c        Si/SiO2 transverse field set to field at midpoint from surface,
c        rather than oxide field refered into silicon. Marden feels this
c        is now the most reasonable thing to do.
c     16-Feb-89  J. Faricelli
c        Modified for 3D binary data
c
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
CC      include 'data.cmn'
      include 'bincodes.cmn'
c
      real psi(ni,nj), e_conc(ni,nj), h_conc(ni,nj), dop_conc(ni,nj),
     1     el(ni,nj), et(ni,nj), ez(ni,nj), ej(ni,nj),
     2     e_quasi_fermi(ni,nj), h_quasi_fermi(ni,nj),
     3     e_lat_cur(ni,nj), e_trans_cur(ni,nj),
     4     h_lat_cur(ni,nj), h_trans_cur(ni,nj),
     5     e_mob(ni,nj), h_mob(ni,nj),
     6     e_charge(ni,nj), h_charge(ni,nj), net_charge(ni,nj),
     7     e_lat_vel(ni,nj), e_trans_vel(ni,nj), e_vel_mag(ni,nj),
     8     h_lat_vel(ni,nj), h_trans_vel(ni,nj), h_vel_mag(ni,nj)
      real nintr, ut, tox, temperature, vs, vb, xl
      integer model
      character*1 chtyp
      logical three_d, nonplanar_boundary, in_binary_file(MENUSIZE)
c
c     locals
      real KBOLTZ, Q
      parameter (KBOLTZ = 1.38066E-23, Q = 1.60218E-19)
c
      real dphmin, dphmax, ccmin, ccmax, h, utch
      real currnt, temp(ARRAYSIZE), logni, cc, vbulk
      integer i, j, imid, nyg, nysem
      logical lmaj
c
c     start of executable code.
c     *************************
c
c     Compute min/max carrier densites
c     This code stolen from MINIMOS4, S. Selberherr et. al.
c
c     Siegi says ccmax is 2.2382/(28.0855*1.6604E-24) = 4.8E22
c     This is the number of silicon atoms per cc.
      ccmax = 4.8E22
c     In Siegi's code, the 1E-20 is variable SMALL. I've used the value
c     he computes for double precision (D) VAX floating point
      ccmin = MAX(1.0E-20, 0.01*MIN(1.0, (nintr/ccmax)*nintr))
      dphmin = ALOG( ccmin )
      dphmax = ALOG( ccmax )
c
c     electrical field at mesh points
c
c     Lateral (X) field
      if( plane_type .ne. YZ_PLANE ) then
         do 10 j=1,nj
c          lateral field at mid points beween mesh...
           do 11 i=1,nxm
11             temp(i) = edxyz(i,1)*( psi(i,j) - psi(i+1,j) )
c
c          lateral field at left and right ends of device...
c          (just use field at left and right mid-points).
           el(1,j) = temp(1)
           el(nxm,j) = temp(nxm)
c          e field at mesh points (simple average as in GENRAT)
           do 12 i=1,nxm-1
c              Introduce special cases for averaging if field is exactly
c              zero on the neighbor point (indicates a conductor).
               if( temp(i) .eq. 0.0 ) then
                  el(i+1,j) = temp(i+1)
               else if( temp(i+1) .eq. 0.0) then
                  el(i+1,j) = temp(i)
               else
                  el(i+1,j) = 0.5*(temp(i) + temp(i+1))
               endif
12             continue
10         continue
         in_binary_file(LAT_E_CODE)   = .TRUE.
      endif
c
c     transverse Y field at mid-points...
      if( plane_type .ne. XZ_PLANE ) then
         do 20 j=1,nytotm
            do 20 i=1,ni
c              store intermediate result one ahead in j...
20             et(i,j+1) = edxyz(j,2)*( psi(i,j) - psi(i,j+1) )
c        transverse field at top and bottom (just use values at inner
c        mid-points). Bottom field is already in et(i,nytot)
c        just do et(i,1)
         do 21 i=1,ni
21          et(i,1) = et(i,2)
c
c        Find gate edge
         nyg = -1
         do 26 j=1,nytot
            if( xyz(j,2) .eq. -tox ) nyg = j
26       continue
c        transverse field at mesh points, use simple average as in GENRAT...
         do 22 j=2,nytotm
c
c           Special case to define field at semiconductor/oxide interface.
c           Use midpoint field as surface field (no average).
c           Marden feels that this is now the better way to compute the
c           surface field (rather than refer the oxide field into the silicon,
c           as he did previously). Also special case transverse field
c           at poly/oxide edge.
c
            if( j .eq. nyoxm+1 ) then
               do 23 i=1,ni
                  et(i,j) = et(i,j+1)
23             continue
            else if( j .eq. nyg ) then
               do 25 i=1,ni
                  et(i,j) = et(i,j+1)
25             continue
            else
               do 24 i=1,ni
                  et(i,j) = 0.5*( et(i,j) + et(i,j+1) )
24             continue
            endif
22       continue
         in_binary_file(TRANS_E_CODE) = .TRUE.
      endif
c
c     Copy transverse component of current at next-to-bottom mesh
c     line of FET to bottom mesh line. Interpolation of current by
c     subroutine INTERY in MINIMOS3 sets bottom value to zero,
c     so we give it a more reasonable value here...
c
      if( plane_type .ne. XZ_PLANE ) then
         do 240 i=1,ni
            e_trans_cur(i,nyoxm+nt) = e_trans_cur(i,nyoxm+nt-1)
            h_trans_cur(i,nyoxm+nt) = h_trans_cur(i,nyoxm+nt-1)
240      continue
      endif
c
c     Z field
      if( plane_type .eq. YZ_PLANE ) then
         do 30 j=1,nytot
c           Z field at mid points beween mesh...
            do 31 i=1,nz-1
31             temp(i) = edxyz(i,3)*( psi(i,j) - psi(i+1,j) )
c
c           Z field at left and right ends of device...
c           (just use field at left and right mid-points).
            ez(1,j)  = temp(1)
            ez(nz,j) = temp(nz-1)
c           E field at mesh points (simple average as in GENRAT)
            do 32 i=1,nz-1
32             ez(i+1,j) = 0.5*(temp(i) + temp(i+1))
30          continue
         in_binary_file(Z_E_CODE)   = .TRUE.
      else if( plane_type .eq. XZ_PLANE ) then
         do 40 j=1,nz-1
            do 40 i=1,nx
c              store intermediate result one ahead in j...
40             ez(i,j+1) = edxyz(j,3)*( psi(i,j) - psi(i,j+1) )
c        transverse field at top and bottom (just use values at inner
c        mid-points). Bottom field is already in ez(i,nz)
c        just do ez(i,1)
         do 41 i=1,nx
41          ez(i,1) = ez(i,2)
c
c        transverse field at mesh points, use simple average as in GENRAT...
         do 42 j=2,nz-1
            do 43 i=1,nx
               ez(i,j) = 0.5*( ez(i,j) + ez(i,j+1) )
43          continue
42       continue
         in_binary_file(Z_E_CODE)   = .TRUE.
      endif
c
c     calculate E in direction of minority carrier j
c     Do for X-Y plane right now...
c
      if( plane_type .eq. XY_PLANE ) then
      do 300 i=2,nxm  ! Si surface of MOSFET
c        count only lateral component, since there is no transverse
c        current under the gate, or along the surface of the FET...
c        (this also avoids having to deal with the definition of the
c        transverse field at the interface where it is discontinuous...)
         if( chtyp .eq.'N' ) then
            currnt=abs(e_lat_cur(i,nyoxm+1))
            if (currnt.eq.0.0) goto 300
            ej(i,nyoxm+1) = el(i,nyoxm+1)*e_lat_cur(i,nyoxm+1)/currnt
         else
            currnt=abs(h_lat_cur(i,nyoxm+1))
            if (currnt.eq.0.0) goto 300
            ej(i,nyoxm+1) = el(i,nyoxm+1)*h_lat_cur(i,nyoxm+1)/currnt
         endif
300   continue

      do 330 j=nyoxm+2,nytot  ! added bottom of MOSFET to computation...
c        at left hand side, el = 0 due to boundary conditions,
c        so only count transverse field in ej
         if( chtyp .eq. 'N' ) then
            currnt=abs(e_trans_cur(1,j))
            if (currnt.eq.0.0) goto 310
            ej(1,j) = et(1,j)*e_trans_cur(1,j)/currnt
         else
            currnt=abs(h_trans_cur(1,j))
            if (currnt.eq.0.0) goto 310
            ej(1,j) = et(1,j)*h_trans_cur(1,j)/currnt
         endif

310      do 320 i=2,nxm
c           interior points (not on boundary)
            if( chtyp .eq. 'N' ) then
               currnt=sqrt(e_lat_cur(i,j)**2+e_trans_cur(i,j)**2)
               if (currnt.eq.0.0) goto 320
               ej(i,j) = (el(i,j)*e_lat_cur(i,j) + 
     1                    et(i,j)*e_trans_cur(i,j))/currnt
            else
               currnt=sqrt(h_lat_cur(i,j)**2+h_trans_cur(i,j)**2)
               if (currnt.eq.0.0) goto 320
               ej(i,j) = (el(i,j)*h_lat_cur(i,j) + 
     1                    et(i,j)*h_trans_cur(i,j))/currnt
            endif
320         continue
c
c        at right hand side, el = 0 due to boundary conditions,
c        so only count transverse field in ej
         if( chtyp .eq. 'N' ) then
            currnt=abs(e_trans_cur(nx,j))
            if (currnt.eq.0.0) goto 330
            ej(nx,j) = et(nx,j)*e_trans_cur(nx,j)/currnt
         else
            currnt=abs(h_trans_cur(nx,j))
            if (currnt.eq.0.0) goto 330
            ej(nx,j) = et(nx,j)*h_trans_cur(nx,j)/currnt
         endif
330      continue
      in_binary_file(EJ_CODE)      = .TRUE.
      endif
c
c     Compute majority carriers when not solved for explicitly
c     by MINIMOS. Assumes constant quasi-Fermi level.
c
c     maj = n_i exp( q(Vbs - psi)/kT )
c
      temperature = ut*Q/KBOLTZ
      logni = ALOG(nintr)
c     Utch flips sign of exponent argument for p-channel devices
      if( chtyp .eq. 'N' ) then
         utch = ut
      else
         utch = - ut
      endif
c     Test for SIMOX (SOI) case where bulk reference is tied to
c     source, not vb
      nysem = nyoxm+ny
      if( (nysem .lt. nytot) .and. (vb .ne. vs) ) then
         h = logni + (- vs)/utch
      else
         h = logni + vb/utch
      endif
c
c     Majority density is included in 3D binary file
c     and is appended to recent (circa 1996) 2D binary files... jvf
      if( chtyp .eq. 'N' ) then
         lmaj = in_binary_file(HOLE_CC_CODE)
      else
         lmaj = in_binary_file(ELEC_CC_CODE)
      endif
      if( model .le. 3 .and. .NOT. three_d .and. .NOT. lmaj ) then
         if( chtyp .eq. 'N' ) then
            do 500 j = nyoxm+1, nysem
               do 510 i = 1,nx
                  h_conc(i,j) = EXP(
     1            MIN( MAX(DPHMIN, (h - psi(i,j)/utch)), DPHMAX))
510            continue
500         continue
c           Zero out concentration in any SOI region
            do 511 j=nysem+1,nytot
               do 511 i=1,nx
                  h_conc(i,j) = 0.0
511         continue
            in_binary_file(HOLE_CC_CODE) = .TRUE.
            in_binary_file(HOLE_QF_CODE) = .TRUE.
         else
            do 520 j = nyoxm+1, nysem
               do 530 i = 1,nx
                  e_conc(i,j) = EXP(
     1            MIN( MAX(DPHMIN, (h - psi(i,j)/utch)), DPHMAX))
530            continue
520         continue
c           Zero out concentration in any SOI region
            do 521 j=nysem+1,nytot
               do 521 i=1,nx
                  e_conc(i,j) = 0.0
521         continue
            in_binary_file(ELEC_CC_CODE) = .TRUE.
            in_binary_file(ELEC_QF_CODE) = .TRUE.
         endif
      endif
c
c     Compute Quasi-Fermi potentials
c
c     n_quasi_fermi = psi - kT/q ln(n/n_i)
c     p_quasi_fermi = psi + kT/q ln(p/n_i)
c
c     where MINIMOS variable ut is kT/q, nintr is n_i
c
c     Note: subtract log of n/p and n_i, to avoid overflow when
c     n_i is very small. Also, guard against small concentrations.
c
c     Use generic array indices (ni,nj) so code is transparent to plane type
c     This also means that we will compute a quasi-Fermi level everywhere
c     the device (even in oxide).
c
c     We also add code to set the quasi-Fermi level for majority carriers
c     to a constant value (Vb) when MODEL < AVL.
c     First, acquire an appropriate "bulk" reference voltage
c     for the SOI (SIMOX) case.
      if( (nysem .lt. nytot) .and. (vb .ne. vs) ) then
         vbulk = vs
      else
         vbulk = vb
      endif
c
      do 400 j = 1,nj
         do 410 i = 1,ni
            if( e_conc(i,j) .lt. ccmin ) then
               cc = ccmin
            else
               cc = e_conc(i,j)
            endif
            if( model .ge. 4 .or. chtyp .eq. 'N' ) then
               e_quasi_fermi(i,j) = psi(i,j) - ut*(ALOG(cc) - logni)
            else
               e_quasi_fermi(i,j) = vbulk
            endif
            if( h_conc(i,j) .lt. ccmin ) then
               cc = ccmin
            else
               cc = h_conc(i,j)
            endif
            if( model .ge. 4 .or. chtyp .eq. 'P' ) then
               h_quasi_fermi(i,j) = psi(i,j) + ut*(ALOG(cc) - logni)
            else
               h_quasi_fermi(i,j) = vbulk
            endif
410      continue
400   continue
      if( .NOT. (three_d .OR. nonplanar_boundary) ) then
c        Set to zero in oxide regions
         do j=1,nyoxm
            do i=1,nx
               e_quasi_fermi(i,j) = 0.0
               h_quasi_fermi(i,j) = 0.0
            enddo
         enddo
      endif
c
c     Compute various charge arrays:
c
c     net charge = (p - n + Nd - Na)
c     electron charge = (Nd - n)
c     hole     charge = (p  - Na)
c     
c     In devices with compensated doping, some charge may get assigned
c     to the wrong catagory (electron/hole) since we only have
c     access to the net doping. Also, in the case
c     of depletion devices, the concept of inversion and bulk charges
c     is murky, since there may be bulk conduction.
c
      do 600 j = 1,nj
         do 610 i = 1,ni
            e_charge(i,j) = MAX(0.0,dop_conc(i,j)) - e_conc(i,j)
            h_charge(i,j) = h_conc(i,j) +  MIN(0.0,dop_conc(i,j))
            net_charge(i,j) = dop_conc(i,j) - e_conc(i,j) + h_conc(i,j)
610       continue
600   continue
      in_binary_file(ELEC_CHARGE_CODE) = .TRUE.
      in_binary_file(HOLE_CHARGE_CODE) = .TRUE.
      in_binary_file(NET_CHARGE_CODE)  = .TRUE.
c
c     Compute carrier velocities (support only XY plane in 3D mode for now)
c
      if( .NOT. (three_d .and. plane_type .ne. XY_PLANE) ) then
c        NOTE WELL: we MUST use the gradient of the quasi-Fermi level
c        when computing the carrier velocity, since the gradient of electric
c        field may lead to "unphysical" velocities in certain regions
c        of the device (e.g. in the source/drain junction depletion
c        regions)
         if( chtyp .eq. 'N' ) then
            call comp_vel(e_lat_vel,e_mob,e_quasi_fermi,xyz(1,1),
     1                    1,nx,nytot,nx)
            call comp_vel(e_trans_vel,e_mob,e_quasi_fermi,xyz(1,2),
     1                    2,nx,nytot,nytot)
            call comp_mag_vel(e_vel_mag,e_trans_vel,e_lat_vel,nx,nytot)
            in_binary_file(LAT_ELEC_VEL_CODE) = .TRUE.
            in_binary_file(TRANS_ELEC_VEL_CODE) = .TRUE.
            in_binary_file(MAG_ELEC_VEL_CODE)  = .TRUE.
         else
            call comp_vel(h_lat_vel,h_mob,h_quasi_fermi,xyz(1,1),
     1                    1,nx,nytot,nx)
            call comp_vel(h_trans_vel,h_mob,h_quasi_fermi,xyz(1,2),
     1                    2,nx,nytot,nytot)
            call comp_mag_vel(h_vel_mag,h_trans_vel,h_lat_vel,nx,nytot)
            in_binary_file(LAT_HOLE_VEL_CODE) = .TRUE.
            in_binary_file(TRANS_HOLE_VEL_CODE) = .TRUE.
            in_binary_file(MAG_HOLE_VEL_CODE)  = .TRUE.
         endif
      endif
      if( model .ge. 4 .and.
     1   .NOT. (three_d .and. plane_type .ne. XY_PLANE) ) then
         if( chtyp .eq. 'N' ) then
            call comp_vel(h_lat_vel,h_mob,h_quasi_fermi,xyz(1,1),
     1                    1,nx,nytot,nx)
            call comp_vel(h_trans_vel,h_mob,h_quasi_fermi,xyz(1,2),
     1                    2,nx,nytot,nytot)
            call comp_mag_vel(h_vel_mag,h_trans_vel,h_lat_vel,nx,nytot)
            in_binary_file(LAT_HOLE_VEL_CODE) = .TRUE.
            in_binary_file(TRANS_HOLE_VEL_CODE) = .TRUE.
            in_binary_file(MAG_HOLE_VEL_CODE)  = .TRUE.
         else
            call comp_vel(e_lat_vel,e_mob,e_quasi_fermi,xyz(1,1),
     1                    1,nx,nytot,nx)
            call comp_vel(e_trans_vel,e_mob,e_quasi_fermi,xyz(1,2),
     1                    2,nx,nytot,nytot)
            call comp_mag_vel(e_vel_mag,e_trans_vel,e_lat_vel,nx,nytot)
            in_binary_file(LAT_ELEC_VEL_CODE) = .TRUE.
            in_binary_file(TRANS_ELEC_VEL_CODE) = .TRUE.
            in_binary_file(MAG_ELEC_VEL_CODE)  = .TRUE.
         endif
      endif
c
c     Determine top of gate electrode
c     Assumes at least two points in the gate
      if( .not. nonplanar_boundary ) then
         do i=1,nx
            if( xyz(i,1) .eq. 0.0 ) nxg1 = i
            if( xyz(i,1) .eq. xl  ) nxg2 = i
         end do
         imid = (nxg2 + nxg1)/2
         do j=1,nytot
            if( xyz(j,2) .eq. -tox ) nyg = j
         end do
         do j=1,nyg
            if( psi(imid,j) .eq. psi(imid,j+1) ) then
               nytop=j
               goto 1000
            endif
         end do
1000     continue
      endif
c
      end
c
      subroutine comp_vel(vel,mob,quasi,mesh,itype,nx,ny,isize)
      implicit none
c
      integer nx, ny, isize, itype
      real vel(nx,ny), quasi(nx,ny), mob(nx,ny),
     1     mesh(isize)
c
      integer i,j
c
c     Compute carrier velocity from mobility and gradient of
c     quasi-Fermi level. Remember: gradient of potential is NEGATIVE
c     of derivative. We use gradient of quasi-Fermi level so we don't
c     get spurious velocites in certain regions of the device
c     (e.g. inside depletion regions in the source/drain junction).
c
      if(itype .eq. 1 ) then
c        Lateral velocity
         do 100 j=1,ny
            do 100 i=2,nx-1
c              Hack for non-planar devices and for dealing with Si/SiO2
c              boundary. If mobility of (i+1) or (i-1) point is exactly zero
c              (in oxide), use one-sided derivative.
               if( mob(i-1,j) .eq. 0.0 ) then
                  vel(i,j) = mob(i,j)*
     1            (quasi(i,j)-quasi(i+1,j))/(mesh(i+1)-mesh(i))
               else if( mob(i+1,j) .eq. 0.0 ) then
                  vel(i,j) = mob(i,j)*
     1            (quasi(i-1,j)-quasi(i,j))/(mesh(i)-mesh(i-1))
               else
                  vel(i,j) = mob(i,j)*
     1            (quasi(i-1,j)-quasi(i+1,j))/(mesh(i+1)-mesh(i-1))
               endif
100      continue
         do 110 j=1,ny
            vel(1,j) = mob(1,j)*
     1         (quasi(1,j)-quasi(2,j))/(mesh(2)-mesh(1))
            vel(nx,j) = mob(nx,j)*
     1         (quasi(nx-1,j)-quasi(nx,j))/(mesh(nx)-mesh(nx-1))
110      continue
      elseif( itype .eq. 2 ) then
c        Transverse velocity
         do 200 j=2,ny-1
            do 200 i=1,nx
c              Hack for non-planar devices and for dealing with Si/SiO2
c              boundary. If mobility of (j-1) point is exactly zero
c              (in oxide), use one-sided derivative.
               if( mob(i,j-1) .eq. 0.0 ) then
                  vel(i,j) = mob(i,j)*
     1            (quasi(i,j)-quasi(i,j+1))/(mesh(j+1)-mesh(j))
               else if( mob(i,j+1) .eq. 0.0 ) then
                  vel(i,j) = mob(i,j)*
     1            (quasi(i,j-1)-quasi(i,j))/(mesh(j)-mesh(j-1))
               else
                  vel(i,j) = mob(i,j)*
     1            (quasi(i,j-1)-quasi(i,j+1))/(mesh(j+1)-mesh(j-1))
               endif
200      continue
         do 210 i=1,nx
            vel(i,1) = mob(i,1)*
     1         (quasi(i,1)-quasi(i,2))/(mesh(2)-mesh(1))
            vel(i,ny) = mob(i,ny)*
     1         (quasi(i,ny-1)-quasi(i,ny))/(mesh(ny)-mesh(ny-1))
210      continue
      endif
c
      return
      end
c
      subroutine comp_mag_vel(vel_mag,vel_lat,vel_trans,nx,ny)
      implicit none
c
      integer nx, ny
      real vel_mag(nx,ny), vel_lat(nx,ny),
     1     vel_trans(nx,ny)
c
      integer i,j
c
      do 100 j=1,ny
         do 100 i=1,nx
            vel_mag(i,j) = sqrt(vel_lat(i,j)**2 + vel_trans(i,j)**2)
100   continue
c
      return
      end
