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     Vector Technologies kindly provided the format of their MCP2D
c     simulator output file.
c     Contact Vector Technologies for information on purchase of their
c     TCAD tools.
c
      subroutine read_mcp2d(batchmode, read_error)
      implicit none
c
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
      include 'io.cmn'
      include 'plot.cmn'
c
c     read Vector Technologies MCP2D output file (ASCII)
c
      logical batchmode, read_error
c
c     locals
c
      call read_mcp2d_file(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_mcp2d_specific
c
      if( lfirst_read ) then
         lfirst_read = .FALSE.
         prev_plane_type = plane_type
         prev_binfilename = binfilename
         prev_binfiletype = binfiletype
      endif
c     For FAIM datasets, re-scale plot limits on each read of a new file
      oned_lfirst = .TRUE.
      twod_lfirst = .TRUE.
      threed_lfirst = .TRUE.
c
      datafile_loaded = .TRUE.
c
10000 return
      end
c
      subroutine read_mcp2d_file(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 Vector Technologies MCP2D output file (ASCII)
c
      logical batchmode, read_error
c
c     locals
c
      integer i, j, k, ioffset, npatch, ilen, ioerr, string_len, ierr,
     1        nxm1, nrecl
      logical linit
      character*(PM_MAX_PATHLEN) expanded_name
c
      write(luttyo,10)
10    format(' ','Reading MCP2D file...')
c
      nrecl=0
      call pm_open_file(binfil, binfilename,
     1               'OLD','FORMATTED','READ','FORTRAN','NONSHARED',
     2               nrecl, ioerr, ierr)
      if( ierr .ne. 0 ) goto 10000
c
c     get stem of file name
      call get_filestem(binfilename, binfilestem)
c
      do 51 j=1,MENUSIZE
         in_binary_file(j) = .FALSE.
51       continue
c
c     Split up data storage into chunks
      call allocate_tri_arrays(ARR2DSIZE,1,ierr)
      if( ierr .ne. 0 ) then
c        Failure during array allocation
         write(luttyo,10021)
         read_error = .TRUE.
         close(unit=binfil)
         return
      endif
c     Zero array (not really necessary)
      call zero_real_array(arrdat,ARR2DSIZE,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
      do i=1,MAX_TRI
         do j=1,3
            tri_nextel(j,i) = -1024
         end do
      end do
c
      tri_npts = 0
      tri_ne   = 0
      tri_nreg = 0
c
c     Read number of sections of the file, nx,ny
      read(binfil,*,err=10010,end=10010,iostat=ioerr) npatch, nx, ny
c
      do  k=1,npatch
         ioffset = tri_npts
         do  j=1,ny
            do  i=1,nx
               tri_npts = tri_npts + 1
               if( tri_npts .gt. ARR2DSIZE ) goto 10020
               read(binfil,*,end=10010,err=10010,iostat=ioerr)
     1            cord(1,tri_npts), cord(2,tri_npts), arrdat(tri_npts)
c              Scale coordinates from microns to cm
               cord(1,tri_npts) = cord(1,tri_npts)*MU2CM
               cord(2,tri_npts) = cord(2,tri_npts)*MU2CM
               tri_matpt(tri_npts) = 1
            end do
         end do
c        Generate triangles from trapezoids
         nxm1 = nx - 1
         do j=1,ny-1
               do i=1,nx-1
               tri_ne = tri_ne + 1
               if( tri_ne+1 .gt. MAX_TRI) goto 10030
               nop(1,tri_ne) = ioffset + (j-1)*nx + i
               nop(2,tri_ne) = ioffset + (j-1)*nx + i + 1
               nop(3,tri_ne) = ioffset + j*nx     + i
               tri2reg(tri_ne) = 1
c              Try to precompute most adjacency information
               if( i .ne. 1 ) then
                  tri_nextel(2,tri_ne) = tri_ne - 1
               endif
               tri_nextel(1,tri_ne) = tri_ne + 1
               if( j .ne. 1 ) then
                  tri_nextel(3,tri_ne) = tri_ne - 2*nxm1 + 1
               endif
               tri_ne = tri_ne + 1
               nop(1,tri_ne) = ioffset + (j-1)*nx + i + 1
               nop(2,tri_ne) = ioffset + j*nx     + i + 1
               nop(3,tri_ne) = ioffset + j*nx     + i
               tri2reg(tri_ne) = 1
c              Try to precompute most adjacency information
               if( i .ne. (nx-1) ) then
                  tri_nextel(3,tri_ne) = tri_ne + 1
               endif
               tri_nextel(2,tri_ne) = tri_ne - 1
               if( j .ne. (ny-1) ) then
                  tri_nextel(1,tri_ne) = tri_ne + 2*nxm1 - 1
               endif
            end do
         end do
      end do
c
c     Dummy region/material number
      imat(1) = 1
      tri_nreg = 1
c     One data array read in
      tri_ndata = 1
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     File has only 2D data
      three_d = .FALSE.
c
c     Load generic mesh pointers, depending on plane type
c
      call set_3d_index_limits
      cut_plane_axis = charb
c
c     This variable must be set by all read routines...
      data_is_trigrid = .TRUE.
c
      in_binary_file(1) = .TRUE.
c
c     Compute adjacency information
      linit = .FALSE.
CC    NOTE: because the boundaries of each patch are formed by distinct
CC    nodes, the code below cannot connect the patches together.
CC    So, for now, we commment this call out -- we've already computed
CC    as much adjacency info by inspection...
CC      write(luttyo,*) 'Computing adjacency information...'
CC      call nextel(nop,tri_nextel,tri_ne,MAX_TRI,linit)
c
      write(luttyo,110)
110   format(' ','...done')
      read_error = .FALSE.
c
      close(unit=binfil)
      goto 10100
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 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 10100
c
c     I/O error on read of file
10010 write(luttyo,10011) ioerr
10011 format(/' ','FORTRAN I/O error #',i5,' on read of file')
      if( ioerr .eq. 24 ) then       ! End-of-file during read
         write(luttyo,10012)
10012    format(' ','Unexpected end-of-file during read...')
      else if( ioerr .eq. 67 ) then  ! Not enough data in record
         write(luttyo,10014)
10014    format(' ','Insufficient data in record'/)
      else if( ioerr .eq. -1 ) then  ! End of file at start of read
         write(luttyo,10015)
10015    format(' ','Unexpected end-of-file')
      endif
      read_error = .TRUE.  ! go back and get new name...
      close(unit=binfil)
      goto 10100
c
10020 continue
      write(luttyo,10021)
10021 format(' ','Error: Number of points exceeds internal storage')
      read_error = .TRUE.
      close(unit=binfil)
      goto 10100
c
10030 continue
      write(luttyo,10031)
10031 format(' ','Error: Number of triangles exceeds internal storage')
      read_error = .TRUE.
      close(unit=binfil)
      goto 10100
c
10100 continue
      return
      end
c
      subroutine load_mcp2d_specific
      implicit none
c
c     Load descriptive arrays with FAIM 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
c     Number of menu items
      max_menu_num = 1
c     Array index of net doping (none in this case)
      net_dop_array_index = 0
c
      call load_desc_arr(1,'Resist image','MCP2D',
     1                   ' ',' ',' ','.mcp2d',
     2                   LIN,OX,NO_JUMP,TWOD,1)
      return
      end
c
      subroutine nextel(nop,tri_nextel,tri_ne,max_tri,linit)
c
c     Subroutine to compute triangle adjacency information
c
      implicit none
c
      integer tri_ne, max_tri
c
c     Triangle to node info
      integer nop(3,max_tri)
c
c     Triangle adjacency information
      integer tri_nextel(3,max_tri)
c
c     Logical flag, do not initialize adjacency array (some values may
c     be precomputed)
      logical linit
c
      integer i,j,k,l,m,icnt,imatch, iside,ipt1,ipt2,jpt1,jpt2
c
c     Initialize array to known state
      if( linit ) then
         do i=1,tri_ne
            do j=1,3
               tri_nextel(j,i) = -1024
            end do
         end do
      endif
c
c     Scan down list of triangles, looking for matches
      do i=1,tri_ne
c
c        Count how many adjacencies entries we've already found...
         icnt = 0
         do m=1,3
            if(tri_nextel(m,i) .gt. 0) icnt = icnt + 1
         end do
c
         do k=i+1,tri_ne
c           If we've already found three matches, we're done for this element
            if( icnt .eq. 3 ) goto 1004
            imatch = 0
            do l=1,3
               do m=1,3
                  if( nop(l,i) .eq. nop(m,k) ) then
                     imatch = imatch + 1
                     if( imatch .eq. 1 ) then
                        ipt1 = l
                        jpt1 = m
                     endif
                     if( imatch .eq. 2 ) then
                        ipt2 = l
                        jpt2 = m
                        goto 1001
                     endif
                  endif
               end do
            end do
1001        continue
            if( imatch .eq. 2 ) then
c              Found a triangle with a matching pair of nodes -> adjacent
c              Store triangle number in tri_nextel, update count
               icnt = icnt + 1
c              Look for matches: {1,2} {2,3} {3,1}
               if( ipt1*ipt2 .eq. 2 ) then
                  iside = 3
               else if( ipt1*ipt2 .eq. 6 ) then
                  iside = 1
               else
                  iside = 2
               endif
               tri_nextel(iside,i) = k
c              Look for matches: {1,2} {2,3} {3,1}
               if( jpt1*jpt2 .eq. 2 ) then
                  iside = 3
               else if( jpt1*jpt2 .eq. 6 ) then
                  iside = 1
               else
                  iside = 2
               endif
               tri_nextel(iside,k) = i
            endif
         end do
1004     continue
c
1000     continue
      end do
c
      return
c
      end
