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
c   Subroutine to read suprem4 ascii structure files
c
c   Parser routines are patterned after those in parser.for
c
      subroutine read_sup4(batchmode,read_error)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
      include 'io.cmn'
      include 'plot.cmn'
c
c     read SUPREM4 ASCII structure files
c
      logical batchmode, read_error
c
c     locals
c
c
      call read_suprem4_structure(batchmode, read_error)
c
      if ( read_error ) goto 10000
c
c     compute some stuff related to mesh...
c
c     set up boundary limits information
      call coord_limits
c
      call load_suprem4_specific
c
      if( lfirst_read ) then
         lfirst_read = .FALSE.
         prev_plane_type = plane_type
         prev_binfilename = binfilename
         prev_binfiletype = binfiletype
      endif
      if( prev_binfiletype .ne. binfiletype ) then
         oned_lfirst = .TRUE.
         twod_lfirst = .TRUE.
         threed_lfirst = .TRUE.
      endif
      datafile_loaded = .TRUE.
c
10000 return
      end
c
      subroutine read_suprem4_structure(batchmode, read_error)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
      include 'io.cmn'
      include 'plot.cmn'
      include 'bincodes.cmn'
c
c     read SUPREM4 ASCII structure file
c
      logical batchmode, read_error
c
c     locals
c
      integer i, j, ilen, ioerr, string_len,
     1        iline, idum(8), ibnd, ireg, itri, i1, ie, idop, ierr,
     2        nrecl
      real    rdum(NUM_DATA_ARRAYS)
      logical ltsup4_new
      character*(PM_MAX_PATHLEN) expanded_name
      character*1024 buffer
      character*1 ch
c
      write(luttyo,10)
10    format(' ','Reading SUPREM4 structure file...')
c
c     SUPREM4 file can write very long lines (up to Unix BUFSIZ=1024)
      nrecl=1024
      call pm_open_file(binfil, binfilename,
     1                  'OLD','FORMATTED','READ','FORTRAN','NONSHARED',
     2                  nrecl, ioerr, ierr)
      if( ierr .ne. 0 ) goto 10000
      datafile_record = 1
c
c     get stem of binary file name
      call get_filestem(binfilename, binfilestem)
c
      do 51 j=1,NUM_DATA_ARRAYS
         in_binary_file(j) = .FALSE.
51    continue
c
c     Split up data storage into chunks
      call allocate_tri_arrays(ARR2DSIZE,NUM_DATA_ARRAYS,ierr)
      if( ierr .ne. 0 ) then
c        Failure during array allocation
         write(luttyo,10041)
         read_error = .TRUE.
         close(unit=binfil)
         return
      endif
c
c     Clear out arrays. Not really necessary.
      call zero_real_array(arrdat,ARR2DSIZE*NUM_DATA_ARRAYS,1)
c
      do i=1,MAX_REGIONS
         imat(i) = 0
      end do
      do i=1,MAX_TRI
         do j=1,3
           nop(j,i) = 0
         end do
      end do
c
      do i=1,ARR2DSIZE
        tri_matpt(i) = -1
      end do
c
      tri_npts = 0
      tri_ne   = 0
      tri_nreg = 0
      ibnd = 0
      ireg = 0
c
c     Read header
      read(binfil,100,end=10010,err=10010,iostat=ioerr) buffer
100   format(a)
c     Check for valid header
c     Note: only minimal support for TMA structure files.
c     This is mainly a hack so the file can at least be read...
      if( buffer(1:1)  .ne. 'v' .or. (
     1    buffer(3:11) .ne. 'SUPREM-IV' .and.
     2    buffer(3:15) .ne. 'TMA TSUPREM-4') ) goto 10020
c
c     Check for new/old TSUPREM4 structure file (version <= 6.0)
      ltsup4_new = .FALSE.
      if( buffer(3:15) .eq. 'TMA TSUPREM-4' ) then
         if( buffer(18:18) .lt. '6' .or. buffer(18:20) .eq. '6.0' ) then
            ltsup4_new = .FALSE.
         else
            ltsup4_new = .TRUE.
         endif
      endif
c
c     Note well: in SUPREM4 code, points, triangles, etc. are numbered
c     from 0..n-1, but in structure file, they are 1..n
c     We use 1..n numbering which is more natural for FORTRAN.
c     It is arbitrary in any case.
c
      iline = 1
200   read(binfil,100,err=10010,end=300,iostat=ioerr) buffer
         iline = iline + 1
         ch = buffer(1:1)
CC         ilen = string_len(buffer)
         buffer(1:1) = ' '
         if( ch .eq. 'c' ) then
c           Coordinates: c pointnum x y h
            read(buffer,*,err=10030,end=10030)
     1           idum(1), rdum(1), rdum(2)
c           Check for point overflow
            if( idum(1) .gt. ARR2DSIZE ) goto 10040
            cord(1,idum(1)) = MU2CM*rdum(1)
            cord(2,idum(1)) = MU2CM*rdum(2)
            tri_npts=tri_npts+1
         else if( ch .eq. 'e' ) then
c           Edges: e followed by 4 integers
c           Not used, skip
         else if( ch .eq. 'r' ) then
c           Regions: r regionnum mattyp
            if( ltsup4_new ) call fix_tsup4_r_line(buffer)
            read(buffer,*,err=10030,end=10030) idum(1), idum(2)
c           Check for region overflow
            if(idum(1) .le. 0 .or. idum(1) .gt. MAX_REGIONS) goto 10050
c           If this is a new region, save number of boundary points
c           in previous region.
            if( ireg .ne. 0 ) tri_nbndpts(ireg) = ibnd
            ireg = idum(1)
            imat(ireg) = idum(2)
            tri_nreg = tri_nreg + 1
            ibnd = 1
         else if( ch .eq. 'b' ) then
c           Boundaries: b edgenum
c           NOTE WELL: I thought this was the boundary point number.
c           It appears to actually be the edge number.. as such
c           it is not useful yet, so don't save it...
            read(buffer,*,err=10030,end=10030) idum(1)
            if( ibnd .eq. 0 .or. ireg .eq. 0 ) goto 10060
CC            if( ibnd .gt. MAX_BOUNDPTS ) goto 10070
CC            tri_bndpts(ibnd,ireg) = idum(1)
            ibnd = ibnd + 1
         else if( ch .eq. 't' ) then
c           Triangles: t trinum regionnum pointi pointj pointk edge1 edge1 edge3
c           Mark is doing something funny, I think the edge codes are really
c           boundary (< 0) or neighbor triangle? (>0) information?
            read(buffer,*,err=10030,end=10030) (idum(i),i=1,8)
            itri = idum(1)
c           Check for triangle store overflow
            if( itri .gt. MAX_TRI ) goto 10080
            tri2reg(itri) = idum(2)
            nop(1,itri) = idum(3)
            nop(2,itri) = idum(4)
            nop(3,itri) = idum(5)
            tri_nextel(1,itri) = idum(6)
            tri_nextel(2,itri) = idum(7)
            tri_nextel(3,itri) = idum(8)
            tri_ne = tri_ne + 1
         else if( ch .eq. 's' ) then
c           solution: s numimpur {list of impuritity numbers}
c           See suprem4/include/impurity.h for list of impurities and id #'s
            if( ltsup4_new ) then
               call parse_tsup4_s_line(buffer,tri_ndata,id2impur,
     1                                 NUM_DATA_ARRAYS,ierr)
               if( ierr .ne. 0 ) goto 10030
            else
               read(buffer,*,err=10030,end=10030)
     1              tri_ndata, (id2impur(i),i=1,tri_ndata)
            endif
            if( tri_ndata .gt. NUM_DATA_ARRAYS ) goto 10090
         else if( ch .eq. 'n' ) then
c           node: n numpoint material {list of impurity values, 1..numimpur}
c           NOTE WELL: SUPREM4 can output multiple 'n' cards for the
c           same point (since we can have multiple values at a point
c           due to segregation at material boundaries). A workaround
c           is to make these duplicates into new points, and then scan
c           the triangle list later to point to these new points.
            read(buffer,*,err=10030,end=10030) i1,idum(1),
     1          (rdum(i),i=1,tri_ndata)
c           Adjust i1: numbering seems to start at 0,  rather than 1...
            i1 = i1 + 1
            if( i1 .le. 0 .or. i1 .gt. tri_npts ) goto 10100
c           Test for point already defined
            if( tri_matpt(i1) .eq. -1 ) then
               tri_matpt(i1) = idum(1)
               do i=1,tri_ndata
                  arrdat(i1+arrind(i)-1) = rdum(i)
               end do
            else
c              Create new "fake" point and adjust triangles
               if (tri_npts+1 .gt. ARR2DSIZE ) goto 10040
               tri_npts = tri_npts + 1
               cord(1,tri_npts) = cord(1,i1)
               cord(2,tri_npts) = cord(2,i1)
               tri_matpt(tri_npts) = idum(1)
               do i=1,tri_ndata
                  arrdat(tri_npts+arrind(i)-1) = rdum(i)
               end do
               do ie=1,tri_ne
c                 Ick. Another 'new' feature of the TMA SUPREM4 structure file
                  if( ltsup4_new ) then
                     if( tri2reg(ie) .eq. idum(1) ) then
                        do j=1,3
                          if( nop(j,ie) .eq. i1 ) nop(j,ie) = tri_npts
                        end do
                     endif
                  else
                      if( imat(tri2reg(ie)) .eq. idum(1) ) then
                        do j=1,3
                          if( nop(j,ie) .eq. i1 ) nop(j,ie) = tri_npts
                        end do
                     endif
                  endif
               end do
            endif
         else if( ch .eq. 'M' ) then
c           Material information, ignore
         else
c           SUPREM4 ignores other lines, so we do the same
         endif
      goto 200
300   continue
c
c     set up adjacency information, what plots (if any) need this?
c     Adjacency info takes a lot of memory
c     call nxtel(p2t,p2tc)
c
c     Set placeholder values for old grid info (probably not needed)
      nx = 1
      ny = 1
      nt = ny
      nyoxm = 0
      nytot= ny
      nxsp = 1
      nxsoxm = 1
      plane_index = 1
      prev_plane_type = plane_type
      plane_type = XY_PLANE
      xyz(1,3) = -xw/2.0
      nxy=nx*ny
      nxytot=nx*nytot
c
c     Load generic mesh pointers, depending on plane type
c
      call set_3d_index_limits
      cut_plane_axis = charb
c     This variable must be set by all read routines...
      data_is_trigrid = .TRUE.
c
c     Compute net total doping, if not already in SUPREM4 file
      idop = 0
      do i=1,tri_ndata
         if( id2impur(i) .eq. 24 ) idop = i
      end do
c
      if( idop .eq. 0 ) then
         tri_ndata = tri_ndata + 1
         if( tri_ndata .gt. NUM_DATA_ARRAYS ) goto 10090
         id2impur(tri_ndata) = -1 ! Our code for net doping
         do i=1,tri_npts
            arrdat(i+arrind(tri_ndata)-1) = 0.0
         end do
c
         do i=1,tri_npts
c           Scan impurities: 2=As, 3=P, 4=Sb, 5=B
            do j=1,tri_ndata
              if( id2impur(j) .eq. 2 .or. id2impur(j) .eq. 3 .or.
     1             id2impur(j) .eq. 4 ) then
                  arrdat(i+arrind(tri_ndata)-1) =
     1            arrdat(i+arrind(tri_ndata)-1) + arrdat(i+arrind(j)-1)
               else if( id2impur(j) .eq. 5 ) then
                  arrdat(i+arrind(tri_ndata)-1) =
     1            arrdat(i+arrind(tri_ndata)-1) - arrdat(i+arrind(j)-1)
               endif
            end do
         end do
      endif
c
      do i=1,tri_ndata
         in_binary_file(i) = .TRUE.
      end do
c
      write(luttyo,110)
110   format(' ','...done')
      read_error = .FALSE.
      close(unit=binfil)
      goto 10200
c
c     errors of commission/omission
c
c     I/O error on open of file
10000 continue
      call expand_filename(binfilename,expanded_name)
      ilen = string_len(expanded_name)
      write(luttyo,10001) ioerr, expanded_name(1:ilen)
10001 format(/' ','FORTRAN I/O error #',i5,' on open of suprem4 file:'/
     1       ' ',a)
      write(luttyo,10002)
10002 format(' ','Check for these errors:'/
     1       ' ','   - Incorrect file, directory or disk name'/
     1       ' ','   - You must have read permission on the file'/)
      read_error = .TRUE.  ! go back and get new name...
      goto 10200
c
c     I/O error on read of file
10010 write(luttyo,10011) ioerr, iline
10011 format(/' ','FORTRAN I/O error #',i5,' on read of file at line',
     1      i10)
      read_error = .TRUE.
      close(unit=binfil)
      goto 10200
c
10020 continue
      write(luttyo,10021)
10021 format(' ','Error: Invalid header in SUPREM4 structure file')
      read_error = .TRUE.
      close(unit=binfil)
      goto 10200
c
10030 ilen=string_len(buffer)
      buffer(1:1) = ch
      write(luttyo,10031) ioerr, iline, buffer(1:ilen)
10031 format(/' ','FORTRAN I/O error #',i5,' on read of file at line',
     1      i10/' ','Invalid or insufficient data on line:'/
     2      ' ',a)
      read_error = .TRUE.
      close(unit=binfil)
      goto 10200
c
10040 continue
      write(luttyo,10041)
10041 format(' ','Error: Number of points exceeds internal storage')
      read_error = .TRUE.
      close(unit=binfil)
      goto 10200
c
10050 continue
      write(luttyo,10051)
10051 format(' ','Error: Number of regions exceeds internal storage')
      read_error = .TRUE.
      close(unit=binfil)
      goto 10200
c
10060 continue
      write(luttyo,10061)
10061 format(' ','Error: boundary information before region defined')
      read_error = .TRUE.
      close(unit=binfil)
      goto 10200
c
10070 continue
      write(luttyo,10071)
10071 format(' ','Error: Number of boundary pts exceeds internal ',
     1           'storage')
      read_error = .TRUE.
      close(unit=binfil)
      goto 10200
c
10080 continue
      write(luttyo,10081)
10081 format(' ','Error: Number of triangles exceeds internal storage')
      read_error = .TRUE.
      close(unit=binfil)
      goto 10200
c
10090 continue
      write(luttyo,10091)
10091 format(' ','Error: Number of impurities exceeds internal storage')
      read_error = .TRUE.
      close(unit=binfil)
      goto 10200
c
10100 continue
      write(luttyo,10101) iline
10101 format(' ','Error: Invalid point number at line:',i10)
      read_error = .TRUE.
      close(unit=binfil)
      goto 10200
c
10200 continue
      return
      end
c
      subroutine load_suprem4_specific
      implicit none
c
c     Load descriptive arrays with SUPREM4 specific data
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'data.cmn'
c
c     Data varies linear or logarithmic
      logical LOG, LIN
      parameter (LOG = .TRUE., LIN = .FALSE.)
c     Data has value in oxide, or not
      logical OX, NO_OX
      parameter (OX = .TRUE., NO_OX = .FALSE.)
c     Data jumps at oxide interfaces
      logical JUMP, NO_JUMP
      parameter (JUMP = .TRUE., NO_JUMP = .FALSE.)
c     Data is 1D or 2D
      logical ONED, TWOD
      parameter (ONED = .TRUE., TWOD = .FALSE.)
c
      integer i, j
      character*35 eltype
      character*20 nickname
      character*10 ext
c
c     Number of menu items
      max_menu_num = tri_ndata
c
c     (2D) array of net doping
      net_dop_array_index = 0
c     File has only 2D data
      three_d = .FALSE.
c
      j=1
      do i=1,tri_ndata
         if( id2impur(i) .eq. -1 .or. id2impur(i) .eq. 24 ) then
            eltype = 'Net doping'
            ext = '.net'
            nickname = 'NET_DOPING'
c           Array index of net doping
            net_dop_array_index = j
         else if( id2impur(i) .eq. 0 ) then
            eltype = 'Vacancies'
            ext = '.vac'
            nickname = 'VACANCIES'
         elseif( id2impur(i) .eq. 1  ) then
            eltype = 'Interstitials'
            ext = '.int'
            nickname = 'INTERSTITIALS'
         elseif( id2impur(i) .eq. 2  ) then
            eltype = 'Chem. Arsenic Concentration'
            ext = '.ch_as_conc'
            nickname = 'AS_CHEM'
         elseif( id2impur(i) .eq. 3  ) then
            eltype = 'Chem. Phosphorous Concentration'
            ext = '.ch_p_conc'
            nickname = 'P_CHEM'
         elseif( id2impur(i) .eq. 4  ) then
            eltype = 'Chem. Antimony Concentration'
            ext = '.ch_sb_conc'
            nickname = 'SB_CHEM'
         elseif( id2impur(i) .eq. 5  ) then
            eltype = 'Chem. Boron Concentration'
            ext = '.ch_b_conc'
            nickname = 'B_CHEM'
         else if( id2impur(i) .eq. 6 ) then
            eltype = 'Electron concentration'
            ext = '.e_conc'
            nickname = 'ELECTRONS'
         else if( id2impur(i) .eq. 7 ) then
            eltype = 'Hole concentration'
            ext = '.h_conc'
            nickname = 'HOLES'
         else if( id2impur(i) .eq. 8 ) then
            eltype = 'X velocity'
            ext = '.xvel'
            nickname = 'XVELOCITY'
         else if( id2impur(i) .eq. 9 ) then
            eltype = 'Y velocity'
            ext = '.yvel'
            nickname = 'YVELOCITY'
         else if( id2impur(i) .eq. 10 ) then
            eltype = 'O2'
            ext = '.o2'
            nickname = 'O2'
         else if( id2impur(i) .eq. 11 ) then
            eltype = 'H2O'
            ext = '.h2o'
            nickname = 'H2O'
         else if( id2impur(i) .eq. 12 ) then
            eltype = 'Intersitial Traps'
            ext = '.itraps'
            nickname = 'ITRAPS'
         else if( id2impur(i) .eq. 14 ) then
            eltype = 'Potential'
            ext = '.psi'
            nickname = 'POTENTIAL'
         else if( id2impur(i) .eq. 15 ) then
            eltype = 'Sxx'
            ext = '.sxx'
            nickname = 'SXX'
         else if( id2impur(i) .eq. 16 ) then
            eltype = 'Syy'
            ext = '.syy'
            nickname = 'SYY'
         else if( id2impur(i) .eq. 17 ) then
            eltype = 'Sxy'
            ext = '.sxy'
            nickname = 'SXY'
         elseif( id2impur(i) .eq. 19 ) then
            eltype = 'Delta Interface Area'
            ext = '.dela'
            nickname = 'DELA'
         elseif( id2impur(i) .eq. 20 ) then
            eltype = 'Active Arsenic Concentration'
            ext = '.as_conc'
            nickname = 'AS'
         elseif( id2impur(i) .eq. 21 ) then
            eltype = 'Active Phosphorous Concentration'
            ext = '.p_conc'
            nickname = 'P'
         elseif( id2impur(i) .eq. 22 ) then
            eltype = 'Active Antimony Concentration'
            ext = '.sb_conc'
            nickname = 'SB'
         elseif( id2impur(i) .eq. 23 ) then
            eltype = 'Active Boron Concentration'
            ext = '.b_conc'
            nickname = 'B'
         elseif( id2impur(i) .eq. 25 ) then
            eltype = 'Grain size'
            ext = '.lgrain'
            nickname = 'LGRAIN'
         elseif( id2impur(i) .eq. 26 ) then
            eltype = 'Grain orientation'
            ext = '.gorient'
            nickname = 'GORIENT'
         elseif( id2impur(i) .eq. 27 ) then
            eltype = 'T poly'
            ext = '.tpoly'
            nickname = 'TPOLY'
         elseif( id2impur(i) .eq. 28 ) then
            eltype = 'Grain Boundary Boron'
            ext = '.gb_boron'
            nickname = 'GB_BORON'
         elseif( id2impur(i) .eq. 29 ) then
            eltype = 'Grain Boundary Arsenic'
            ext = '.gb_arsenic'
            nickname = 'GB_ARSENIC'
         elseif( id2impur(i) .eq. 30 ) then
            eltype = 'Grain Boundary Phosphorous'
            ext = '.gb_phos'
            nickname = 'GB_PHOS'
         elseif( id2impur(i) .eq. 31 ) then
            eltype = 'Clustered Interstitials'
            ext = '.int_clust'
            nickname = 'INTERSTITIAL_CLUSTER'
         else
            write(eltype,100) id2impur(i)
100         format('SUPREM4 ID number ',i5)
            if( id2impur(i) .lt. 10 ) then
               write(nickname,101) id2impur(i)
101            format('SUPREM_ID_',i1)
            elseif( id2impur(i) .lt. 100) then
               write(nickname,102) id2impur(i)
102            format('SUPREM_ID_',i2)
            else
               write(nickname,103) id2impur(i)
103            format('SUPREM_ID_',i3)
            endif
            ext = '.unknown'
         endif
         if( id2impur(i) .eq. 14 ) then
            call load_desc_arr(j, eltype, nickname,
     1                        'V','V cm','V cm^{2}',ext,
     2                         LIN,OX,NO_JUMP,TWOD,j)
            j = j + 1
         else if( id2impur(i) .eq. 6 .or. id2impur(i) .eq. 7 ) then
            call load_desc_arr(j, eltype, nickname,
     1                        'cm^{-3}','cm^{-2}','cm^{-1}',ext,
     2                         LOG,NO_OX,NO_JUMP,TWOD,j)
            j = j + 1
         else if( id2impur(i) .eq. 8 .or. id2impur(i) .eq. 9 ) then
c           Think the units are micron per second
            call load_desc_arr(j, eltype, nickname,
     1                        '\mu /s','\mu ^{2}/s','\mu ^{3}/s',ext,
     2                         LIN,OX,NO_JUMP,TWOD,j)
            j = j + 1
         else if( id2impur(i) .ge. 15 .and. id2impur(i) .le. 17 ) then
c           Think the units are dynes/cm^2
            call load_desc_arr(j, eltype, nickname,
     1                        'dynes/cm^{2}','dynes/cm','dynes',ext,
     2                         LIN,OX,NO_JUMP,TWOD,j)
            j = j + 1
         else if( id2impur(i) .eq. 19 ) then
c           Not sure the units are square micron
            call load_desc_arr(j, eltype, nickname,
     1              '\mu ^{2}','\mu ^{2} cm','\mu ^{2} cm^{2}',ext,
     2               LIN,OX,NO_JUMP,TWOD,j)
            j = j + 1
         else if( id2impur(i) .eq. 25 ) then
            call load_desc_arr(j, eltype, nickname,
     1                        'cm','cm^{2}','cm^{3}',ext,
     2                         LIN,OX,NO_JUMP,TWOD,j)
            j = j + 1
         else if( id2impur(i) .eq. 26 ) then
            call load_desc_arr(j, eltype, nickname,
     1                        ' ',' ',' ',ext,
     2                         LIN,OX,NO_JUMP,TWOD,j)
            j = j + 1
         else
            call load_desc_arr(j, eltype, nickname,
     1                        'cm^{-3}','cm^{-2}','cm^{-1}',ext,
     2                         LOG,OX,NO_JUMP,TWOD,j)
            j = j + 1
         endif
      end do
c
      return
      end
c
      subroutine coord_limits
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'data.cmn'
      include 'mesh.cmn'
c
      integer i, j
      real dmin, dmax
c
      if( data_is_trigrid ) then
         do 100 j=1,2
            dmin = PM_MAX_FLOAT
            dmax = - PM_MAX_FLOAT
            do 200 i=1,tri_npts
               if( cord(j,i) .gt. dmax ) dmax = cord(j,i)
               if( cord(j,i) .lt. dmin ) dmin = cord(j,i)
200         continue
            xyzlimits(1,j) = dmin
            xyzlimits(2,j) = dmax
100      continue
c        Dummy values for Z coordinate
         xyzlimits(1,3) = -1.0
         xyzlimits(2,3) =  1.0
      else
         xyzlimits(1,1) = xyz(1,1)
         xyzlimits(2,1) = xyz(nx,1)
         xyzlimits(1,2) = xyz(1,2)
         xyzlimits(2,2) = xyz(nytot,2)
         xyzlimits(1,3) = xyz(1,3)
         xyzlimits(2,3) = xyz(nz,3)
      endif
c
      return
      end
c
      subroutine fix_tsup4_r_line(buffer)
      implicit none
      character*(*) buffer
c
      integer ilen, i, istart, iend
      ilen = LEN(buffer)
c
c     Skip initial white space
      i = 1
      do while( buffer(i:i) .eq. ' ' .and. i .lt. ilen )
        i = i + 1
      end do
c
c     Skip region number
      do while( buffer(i:i) .ne. ' ' .and. i .lt. ilen )
        i = i + 1
      end do
c
c     Skip white space
      do while( buffer(i:i) .eq. ' ' .and. i .lt. ilen )
        i = i + 1
      end do
c
c     Replace material name with a numeric value
      istart = i
      do while( buffer(i:i) .ne. ' ' .and. i .lt. ilen )
         i = i + 1
      end do
      iend = i
      if( buffer(istart:iend) .eq. 'silicon' ) then
         buffer(istart:iend) = '1'
      else if( buffer(istart:iend) .eq. 'oxide' ) then
         buffer(istart:iend) = '-1'
      else if( buffer(istart:iend) .eq. 'nitride' ) then
         buffer(istart:iend) = '-2'
      else if( buffer(istart:iend) .eq. 'polysilicon' ) then
         buffer(istart:iend) = '3'
      else if( buffer(istart:iend) .eq. 'aluminum' ) then
         buffer(istart:iend) = '-5'
      else if( buffer(istart:iend) .eq. 'titanium' ) then
         buffer(istart:iend) = '-5'
      else if( buffer(istart:iend) .eq. 'electrode' ) then
         buffer(istart:iend) = '-5'
      else if( buffer(istart:iend) .eq. 'gaas' ) then
         buffer(istart:iend) = '2'
      else if( buffer(istart:iend) .eq. 'semiconductor' ) then
         buffer(istart:iend) = '4'
      else if( buffer(istart:iend) .eq. 'sige' ) then
         buffer(istart:iend) = '5'
      else if( buffer(istart:iend) .eq. 'algaas' ) then
         buffer(istart:iend) = '6'
      else if( buffer(istart:iend) .eq. 'ge' ) then
         buffer(istart:iend) = '7'
      else if( buffer(istart:iend) .eq. 'sic' ) then
         buffer(istart:iend) = '8'
      else if( buffer(istart:iend) .eq. 'TiSi2' ) then
         buffer(istart:iend) = '9'
      else if( buffer(istart:iend) .eq. 'sapphire' ) then
         buffer(istart:iend) = '-3'
      else if( buffer(istart:iend) .eq. 'insulator' ) then
         buffer(istart:iend) = '-4'
      else if( buffer(istart:iend) .eq. 'oxynitride' ) then
         buffer(istart:iend) = '-4'
      else if( buffer(istart:iend) .eq. 'photoresist' ) then
         buffer(istart:iend) = '-4'
      else if( buffer(istart:iend) .eq. 'e_photoresist' ) then
         buffer(istart:iend) = '-4'
      else if( buffer(istart:iend) .eq. 'ambient' ) then
         buffer(istart:iend) = '-4'
      else
c        Punt
         buffer(istart:iend) = '1'
      endif
c
      return
      end
c
      subroutine parse_tsup4_s_line(buffer, ndata, id2impur, impur_size,
     1                              ierr)
      implicit none

      integer ndata, impur_size, id2impur(impur_size)
      character*(*) buffer
      integer ierr
c
      character*50 key
      integer i, ipos, imax
c
      ierr = 0
      ipos = 1
      imax = 32 ! this should be equal to number of known quantities + 1
      call get_keyword(buffer,key,ipos)
      if( ipos .eq. 0 ) goto 1000
      read(key,*,end=1000,err=1000) ndata
c
      do i=1,ndata
c        Note: I can't find documentation on all the possible tsuprem4
c        names used in the new file format. Those I have checked by
c        running tsuprem4 are noted with a *** CHECKED *** comment
         call get_keyword(buffer,key,ipos)
         if( ipos .eq. 0 ) goto 1000
         if( key .eq. 'net' ) then
            id2impur(i) = 24
         else if( key .eq. 'total' ) then
            id2impur(i) = 2
         else if( key .eq. 'donor' ) then
            id2impur(i) = 3
         else if( key .eq. 'acceptor' ) then
            id2impur(i) = 4
         else if( key .eq. 'v' ) then
            id2impur(i) = 14
         else if( key .eq. 'n' ) then
            id2impur(i) = 6
         else if( key .eq. 'p' ) then
            id2impur(i) = 7
         else if( key .eq. 'vacancy' ) then   ! *** CHECKED ***
            id2impur(i) = 0
         else if( key .eq. 'interstitial' ) then   ! *** CHECKED ***
            id2impur(i) = 1
         else if( key .eq. 'arsenic' ) then   ! *** CHECKED ***
            id2impur(i) = 2
         else if( key .eq. 'phosphorus' ) then   ! *** CHECKED ***
            id2impur(i) = 3
         else if( key .eq. 'antimony' ) then   ! *** CHECKED ***
            id2impur(i) = 4
         else if( key .eq. 'boron' ) then   ! *** CHECKED ***
            id2impur(i) = 5
         else if( key .eq. 'vac_s' ) then
            id2impur(i) = 14 ! Vacancy_Steady
         else if( key .eq. 'inter_s' ) then
            id2impur(i) = 15 ! Interstitial_Steady
         else if( key .eq. 'x.vel' ) then   ! *** CHECKED ***
            id2impur(i) = 8
         else if( key .eq. 'y.vel' ) then   ! *** CHECKED ***
            id2impur(i) = 9
         else if( key .eq. 'o2' ) then
            id2impur(i) = 10
         else if( key .eq. 'h2o' ) then   ! *** CHECKED ***
            id2impur(i) = 11
         else if( key .eq. 'trap' ) then
            id2impur(i) = 12 ! Traps
         else if( key .eq. 'sxx' ) then   ! *** CHECKED ***
            id2impur(i) = 15
         else if( key .eq. 'syy' ) then   ! *** CHECKED ***
            id2impur(i) = 16
         else if( key .eq. 'sxy' ) then   ! *** CHECKED ***
            id2impur(i) = 17
         else if( key .eq. 'a_arsenic' ) then   ! *** CHECKED ***
            id2impur(i) = 20 ! Arsenic_Active
         else if( key .eq. 'a_phosphorus' ) then   ! *** CHECKED ***
            id2impur(i) = 21 ! Phosphorus_Active
         else if( key .eq. 'a_antimony' ) then   ! *** CHECKED ***
            id2impur(i) = 22 ! Antimony_Active
         else if( key .eq. 'a_boron' ) then   ! *** CHECKED ***
            id2impur(i) = 23 ! Boron_Active
         else if( key .eq. 'lgrain' ) then
            id2impur(i) = 25 ! grain size
         else if( key .eq. 'gorient' ) then
            id2impur(i) = 26 ! grain orientation
         else if( key .eq. 'tpoly' ) then
            id2impur(i) = 27 ! tpoly ???
         else if( key .eq. 'gb_boron' ) then
            id2impur(i) = 28 ! grain boundary boron
         else if( key .eq. 'gb_arsenic' ) then
            id2impur(i) = 29 ! grain boundary arsenic
         else if( key .eq. 'gb_phosphorus' ) then
            id2impur(i) = 30 ! grain boundary phos
         else if( key .eq. 'cl_interst' ) then
            id2impur(i) = 31 ! Clustered interstitials
         else
c           If you add more quantities, then bump up imax above
            id2impur(i) = imax
            imax = imax + 1
         endif
      end do
c
      return
c
1000  continue
      ierr = 1
      return
      end
