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 read_suprem3(type, 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 SUPREM3 save files (structure and ASCII save file formats)
c
      integer type
      logical read_error
c
c     locals
c
      integer i
c
      if( type .eq. 1 ) then
         call read_suprem3_ascii(read_error)
      else
         call read_suprem3_str(read_error)
      endif
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
      nxm = nx - 1
      nym = ny - 1
      nytotm = nytot - 1
c
c     inverse mesh spacing...
c
      do 210 i=1,nxm
210      edxyz(i,1) = 1.0/(xyz(i+1,1) - xyz(i,1))
      do 211 i=1,nytotm
211      edxyz(i,2) = 1.0/(xyz(i+1,2) - xyz(i,2))
      do 212 i=1,nz-1
212      edxyz(i,3) = 1.0/(xyz(i+1,3) - xyz(i,3))
c
c     average mesh spacing...
c
      dxyzm(1,1)  = (xyz(2,1)  - xyz(1,1))/2.0
      dxyzm(nx,1) = (xyz(nx,1) - xyz(nx-1,1))/2.0
      do 220 i=2,nxm
220      dxyzm(i,1) = (xyz(i+1,1) - xyz(i-1,1))/2.0
c
      dxyzm(1,2)     = (xyz(2,2)  - xyz(1,2))/2.0
      dxyzm(nytot,2) = (xyz(nytot,2) - xyz(nytot-1,2))/2.0
      do 221 i=2,nytotm
221      dxyzm(i,2) = (xyz(i+1,2) - xyz(i-1,2))/2.0
c
      dxyzm(1,3)     = (xyz(2,3)  - xyz(1,3))/2.0
      dxyzm(nz,3) = (xyz(nz,3) - xyz(nz-1,3))/2.0
      do 222 i=2,nz-1
222      dxyzm(i,3) = (xyz(i+1,3) - xyz(i-1,3))/2.0
c
      call load_suprem3_specific(type)
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_suprem3_ascii(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 SUPREM3 ASCII save file (e.g. PRINT MINIMOS)
c
      logical read_error
c
c     locals
c
      integer i, j, ilen, ioerr, ierr, nrecl, string_len
      character*(PM_MAX_PATHLEN) expanded_name
      character*132 buffer
c
      write(luttyo,10)
10    format(' ','Reading SUPREM "PRINT MINIMOS" 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
      datafile_record = 1
c
c     get stem of binary file name
      call get_filestem(binfilename, binfilestem)
c
c     clear out arrays
      do 45 j=1,3
         do 45 i=1,ARRAYSIZE
            xyz(i,j) = 0.0
            edxyz(i,j) = 0.0
            dxyzm(i,j) = 0.0
45    continue
      do 51 j=1,MENUSIZE
         in_binary_file(j) = .FALSE.
51       continue
      call zero_real_array(arr2d,COMP_ARRAYSIZE,MENUSIZE)
c
c     Skip header
      read(binfil,100,end=10010,err=10010,iostat=ioerr) buffer
      read(binfil,100,end=10010,err=10010,iostat=ioerr) buffer
c     Check for valid header
      if( buffer(2:12) .ne. '!     depth' ) goto 10020
100   format(a)
c
c     
      i = 0
200   read(binfil,100,err=10010,end=10010,iostat=ioerr) buffer
         i = i + 1
         if( buffer(2:15) .eq. '! Total Active' ) goto 300
         if( i .gt. COMP_ARRAYSIZE ) then
            write(luttyo,201) COMP_ARRAYSIZE
201         format(' ','Save file size exceeds POSTMINI array storage'/
     1             ' ','Data will be truncated after ',i5,' points')
            goto 300
         endif
         read(buffer,*) xyz(i,1), arr2d(i,5), arr2d(i,1), arr2d(i,2),
     1                  arr2d(i,3), arr2d(i,4)
         xyz(i,1) = xyz(i,1)*MU2CM
      goto 200
300   continue
      nx = i - 1
c     Fake second dimension
      ny = 1
      xyz(1,2) = 0.0
c
      nz = 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
c     Since SUPREM3 has only 1D data in the Y direction, flip axes
c     so labeling comes out right.
      charv = 'X'
      charh = 'Y'
      cut_plane_axis = charb
c     This variable must be set by all read routines...
      data_is_trigrid = .FALSE.
c
      do i=1,5
         in_binary_file(i) = .TRUE.
      end do
c
      write(luttyo,110)
110   format(' ','...done')
      read_error = .FALSE.
      close(unit=binfil)
      goto 10100
c
c     errors of commission/omission
c
c     I/O error on open of binary 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 binary 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
c     I/O error on read of binary 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...'/
     1          ' ','File may not be a SUPREM3 ASCII save file'/)
         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'/
     1     ' ','File may not have been written completely by SUPREM3'/)
      endif
      read_error = .TRUE.  ! go back and get new name...
      close(unit=binfil)
      goto 10100
c
10020 continue
      write(luttyo,10021)
10021 format(' ','Error: Invalid header in Suprem-III doping file')
      read_error = .TRUE.  ! go back and get new name...
      close(unit=binfil)
      goto 10100
c
10100 continue
      return
      end
c
      subroutine load_suprem3_specific(type)
      implicit none
c
c     Load descriptive arrays with SUPREM3 specific data
c     If type = 1, then input was from "PRINT MINIMOS" file (5 entries)
c     Otherwise, input was from structure file (14 entries)
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 type
c
c     Number of menu items
      if( type .eq. 1 ) then
         max_menu_num = 5
      else
         max_menu_num = 14
      endif
c
c     (2D) array of net doping, not applicable to SUPREM3.
      net_dop_array_index = 0
c     File has only 1D data
      three_d = .FALSE.
c
      call load_desc_arr(1, 'Active Boron Concentration','B',
     1                   'cm^{-3}','cm^{-2}','cm^{-1}','.b_conc',
     2                   LOG,OX,NO_JUMP,ONED,1)
c
      call load_desc_arr(2, 'Active Phosphorous Concentration','P',
     1                   'cm^{-3}','cm^{-2}','cm^{-1}','.p_conc',
     2                   LOG,OX,NO_JUMP,ONED,2)
c
      call load_desc_arr(3, 'Active Arsenic Concentration','AS',
     1                   'cm^{-3}','cm^{-2}','cm^{-1}','.as_conc',
     2                   LOG,OX,NO_JUMP,ONED,3)
c
      call load_desc_arr(4, 'Active Antimony Concentration','SB',
     1                   'cm^{-3}','cm^{-2}','cm^{-1}','.sb_conc',
     2                   LOG,OX,NO_JUMP,ONED,4)
c
      if( type .eq. 1 ) then
         call load_desc_arr(5, 'Net Concentration','NET_DOPING',
     1                      'cm^{-3}','cm^{-2}','cm^{-1}','.net_conc',
     2                      LOG,OX,NO_JUMP,ONED,5)
      else
         call load_desc_arr(5, 'Chem. Boron Concentration','B_CHEM',
     1                      'cm^{-3}','cm^{-2}','cm^{-1}','.ch_b_conc',
     2                      LOG,OX,NO_JUMP,ONED,5)
         call load_desc_arr(6, 'Chem. Phosphorous Concentration',
     1                      'P_CHEM',
     2                      'cm^{-3}','cm^{-2}','cm^{-1}','.ch_p_conc',
     3                      LOG,OX,NO_JUMP,ONED,6)
         call load_desc_arr(7, 'Chem. Arsenic Concentration', 'AS_CHEM',
     1                      'cm^{-3}','cm^{-2}','cm^{-1}','.ch_as_conc',
     2                      LOG,OX,NO_JUMP,ONED,7)
         call load_desc_arr(8, 'Chem. Antimony Concentration','SB_CHEM',
     1                      'cm^{-3}','cm^{-2}','cm^{-1}','.ch_sb_conc',
     2                      LOG,OX,NO_JUMP,ONED,8)
         call load_desc_arr(9, 'Inactive Boron Conc.','B_INACTIVE',
     1                      'cm^{-3}','cm^{-2}','cm^{-1}','.in_b_conc',
     2                      LOG,OX,NO_JUMP,ONED,9)
         call load_desc_arr(10, 'Inactive Phosphorous Conc.',
     1                      'P_INACTIVE',
     1                      'cm^{-3}','cm^{-2}','cm^{-1}','.in_p_conc',
     2                      LOG,OX,NO_JUMP,ONED,10)
         call load_desc_arr(11, 'Inactive Arsenic Conc.', 'AS_INACTIVE',
     1                      'cm^{-3}','cm^{-2}','cm^{-1}','.in_as_conc',
     2                      LOG,OX,NO_JUMP,ONED,11)
         call load_desc_arr(12, 'Inactive Antimony Conc.','SB_INACTIVE',
     1                      'cm^{-3}','cm^{-2}','cm^{-1}','.in_sb_conc',
     2                      LOG,OX,NO_JUMP,ONED,12)
         call load_desc_arr(13, 'Net Active Concentration','NET_DOPING',
     1                      'cm^{-3}','cm^{-2}','cm^{-1}','.net_conc',
     2                      LOG,OX,NO_JUMP,ONED,13)
         call load_desc_arr(14, 'Net Chem. Concentration','NET_CHEM',
     1                      'cm^{-3}','cm^{-2}','cm^{-1}','.net_chem',
     2                      LOG,OX,NO_JUMP,ONED,14)
      endif
      return
      end
c
      subroutine read_suprem3_str(read_error)
      implicit none
c
c     Read SUPREM3 binary structure file
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
      include 'io.cmn'
      include 'plot.cmn'
c
      logical read_error
c
c     locals
c
      integer MAXIMP, MAXGRID,  MAXGRIDE, MAXLAYERS, MAXMAT, MAXMATP1
      parameter (MAXIMP = 4)
      parameter (MAXGRID = 500)
      parameter (MAXGRIDE = 600)
      parameter (MAXLAYERS = 10)
      parameter (MAXMAT = 10)
      parameter (MAXMATP1 = MAXMAT + 1)
c
      integer revcod, nimps, iused(MAXIMP), mattyp(MAXLAYERS),
     1        top(MAXLAYERS),bottom(MAXLAYERS),layrsp(MAXLAYERS), nlayrs
      logical materl(MAXMATP1), impur(MAXIMP)
      real layrth(MAXLAYERS), h(MAXGRIDE), dum1, dum2, dum3, dum4, dum5
      character*2 fdatyp
c
      integer i, j, k, ilen, ioerr, string_len, npts, isi_layer, ierr,
     1        nrecl
      real thick, xtotal, imptyp(4)
      character*(PM_MAX_PATHLEN) expanded_name
      data imptyp/-1.0, 1.0, 1.0, 1.0/
c
      write(luttyo,10)
10    format(' ','Reading SUPREM3 structure file...',$)
c
      nrecl=0
      call pm_open_file(binfil, binfilename,
     1               'OLD','UNFORMATTED','READ','NONE','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
c     clear out arrays
      do 45 j=1,3
         do 45 i=1,ARRAYSIZE
            xyz(i,j) = 0.0
            edxyz(i,j) = 0.0
            dxyzm(i,j) = 0.0
45    continue
      do 51 j=1,MENUSIZE
         in_binary_file(j) = .FALSE.
51       continue
      call zero_real_array(arr2d,COMP_ARRAYSIZE,MENUSIZE)
c
c     Header
      read(binfil,end=10010,err=10010,iostat=ioerr) fdatyp
      if( fdatyp .NE. 'A=' .AND. fdatyp .NE. 'S=' ) goto 10020
c
      read(binfil,end=10010,err=10010,iostat=ioerr) revcod
      if( revcod .NE. 8628 ) goto 10030
c
      if( MAXGRID .gt. COMP_ARRAYSIZE ) then
         write(luttyo,10050)
10050    format(' ','SUPREM save file too large for POSTMINI arrays')
         read_error = .TRUE.  ! go back and get new name...
         close(unit=binfil)
         goto 10100
      endif
c
c     The number of records in the structure definition.
      read(binfil, err=10010, end=10010, iostat=ioerr)  j
c
c     The number of layers and the materials in use array.
      read(binfil, err=10010, end=10010, iostat=ioerr) nlayrs,
     1   (materl(i),i=1,MAXMATP1)
c    
c     For each layer, the material type, orientation, top pointer, bottom
c     pointer, number of spaces, thickness, nominal grid spacing, location
c     of nominal grid spacing, minimum grid spacing, grain size.
      read(binfil, err=10010, end=10010, iostat=ioerr)
     1   (mattyp(i),dum1,top(i),bottom(i),layrsp(i),
     2    layrth(i),dum2,dum3,dum4,dum5,i=1,MAXLAYERS)
c    
c     For each impurity and each layer, the integrated dopant, interior
c     concentration of polysilicon grains.
      read(binfil, err=10010, end=10010, iostat=ioerr)
     1    ((dum1,dum2,j=1,MAXIMP),i=1,MAXLAYERS)
c    
c     More stuff that we can skip (pqd,impas,tempC,otempK,tempK,kT)
      read(binfil, err=10010, end=10010, iostat=ioerr)
c    
c     The impurities used.
      read(binfil, err=10010, end=10010, iostat=ioerr)
     1   (impur(i),i=1,MAXIMP)
c    
c     The grid spacing (NOT position!).
      read(binfil, err=10010, end=10010, iostat=ioerr)
     1   (h(i),i=top(nlayrs),MAXGRID)
c    
c     For each impurity, the concentration distribution.
      nimps = 0
      do i = 1, MAXIMP 
         if (impur(i)) then
            nimps = nimps + 1 
            iused(nimps) = i 
         endif
      end do
c     Load chemical and active concentrations
      npts = MAXGRID - top(nlayrs) + 1
      do j = 1, nimps
        read(binfil, err=10010, end=10010, iostat=ioerr)
     1  (arr2d(i,iused(j)+4),arr2d(i,iused(j)),i=1,npts)
c        Compute inactive concentration
         do i = 1,npts
           arr2d(i,iused(j)+8) = arr2d(i,iused(j)+4) - arr2d(i,iused(j))
         end do
      end do
c
      nx = MAXGRID - top(nlayrs) + 1
c
c     Compute distance array from mesh spacings
c     Start at layer 1 (deepest layer) and look for a silicon layer.
c     Offset rest of thickness of layers so silicon surface is at X=0.0
c     If no silicon layer is found, then do not offset coordinates.
      isi_layer = nlayrs
      do 60 i=1,nlayrs
         if( mattyp(i) .eq. 1 ) then
c           Found silicon layer (mattyp == 1)
            isi_layer = i
            goto 61
         endif
60    continue
61    continue
c
      thick = 0.0
      do i=isi_layer+1,nlayrs
         thick = thick - layrth(i)
      end do
      xtotal = 0.0
      j = 0
      sup3_nlayers = nlayrs
      sup3_interface_coord(1) = thick*MU2CM
      k = 1
      do i=top(nlayrs),bottom(1)
         j = j + 1
         if( h(i) .eq. 0.0 ) then
c           Offset mesh line to fake discontinuity at interface
            xyz(j,1) = (xtotal + thick - 1.0E-4)*MU2CM
c           Record position of interface
            k = k + 1
            sup3_interface_coord(k) = (xtotal + thick)*MU2CM
         else
            xyz(j,1) = (xtotal + thick)*MU2CM
         endif
         xtotal = xtotal + h(i)
      end do
c     Compute net active concentration
      do i=1,npts
         do j=1,MAXIMP
            arr2d(i,13) = arr2d(i,13) + imptyp(j)*arr2d(i,j)
         end do
      end do
c     Compute net total concentration
      do i=1,npts
         do j=1,MAXIMP
            arr2d(i,14) = arr2d(i,14) + imptyp(j)*arr2d(i,j+4)
         end do
      end do
c
c     Store material types in reverse order (from surface down)
c     Material types: Si=1, SiO2=2, Poly=3
      do i=1,nlayrs
         sup3_mattyp(i) = mattyp(nlayrs-i+1)
      end do
c
      do i=1,4
         if( impur(i) ) then
            in_binary_file(i)   = .TRUE.
            in_binary_file(i+4) = .TRUE.
            in_binary_file(i+8) = .TRUE.
         endif
      end do
      in_binary_file(13) = .TRUE.
      in_binary_file(14) = .TRUE.
c
c     Fake second dimension...
      ny = 1
      nz = 1
      nt = ny
      nyoxm = 0
      nytot= ny
      nxsp = 1
      nxsoxm = 1
      plane_index = 1
      plane_type = XY_PLANE
      xyz(1,3) = -xw/2.0
c
      nxy=nx*ny
      nxytot=nx*nytot
c
c     Load generic mesh pointers, depending on plane type
c
      call set_3d_index_limits
c     Since SUPREM3 has only 1D data in the Y direction, flip axes
c     so labeling comes out right.
      charv = 'X'
      charh = 'Y'
      cut_plane_axis = charb
c
      write(luttyo,110)
110   format(' ','...done')
      read_error = .FALSE.
      close(unit=binfil)
      goto 10100
c
c     errors of commission/omission
c
c     I/O error on open of binary 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 binary 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
c     I/O error on read of binary file
10010 write(luttyo,10011) ioerr
10011 format(/' ','FORTRAN I/O error #',i5,' on read of binary file')
      if( ioerr .eq. 24 ) then       ! End-of-file during read
         write(luttyo,10012)
10012    format(' ','Unexpected end-of-file during read...'/
     1          ' ','File may not be a MINIMOS binary file'/)
      else if( ioerr .eq. 35 ) then  ! Not binary file format
         write(luttyo,10013)
10013    format(' ','File is probably not a FORTRAN binary file...'/)
      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'/
     1     ' ','File may not have been written completely by SUPREM3'/)
      endif
      read_error = .TRUE.  ! go back and get new name...
      close(unit=binfil)
      goto 10100
c
10020 continue
      write(luttyo,10021)
10021 format(' ','Error: file does not contain structure information')
      read_error = .TRUE.
      close(binfil)
      goto 10100
c
10030 continue
      write(luttyo,10031) revcod
10031 format(' ','Error: old structure file: revision code ',i5)
      read_error = .TRUE.
      close(binfil)
      goto 10100
c
10100 continue
      return
      end
