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_2dop(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 MINIMOS 2DOP binary file
c
      logical batchmode, read_error
c
c     locals
c
      logical lint
      integer i
c
      call read_2dop_bin(batchmode,read_error,lint)
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_2dop_specific(lint)
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_2dop_bin(batchmode,read_error,lint)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
      include 'io.cmn'
      include 'plot.cmn'
      include 'units.cmn'
      include 'bincodes.cmn'
c
c     Large arrays used by incore version of code
      include 'incore.cmn'
c
c     Read double precision MINIMOS 2DOP binary file
c
c     For 3D case, ASSUMES that the X and Y mesh for each X-Y plane
c     are identical.
c
c     Also reads associated interface file (.INT), if it exists
c
      logical batchmode, read_error, lint
c
c     locals
c
      integer i, j, k, ilen, ioerr, string_len,
     1        ierr, nrecl, npiece, nx_d, ny_d,
     2        nx_prev, nx_d_prev, nplanes,
     3        plane_cnt, ioff1, ioff2, nx_lhs
      real real_vec(2), ddrain, dsource, offset
      logical lfirst, plane_select_err, alloc_err
      character*10 format_array(2)
      character*(PM_MAX_PATHLEN) expanded_name
c
      logical threed_loaded
      save threed_loaded
      data threed_loaded/.FALSE./
c
c     For interface file read code
      logical lexist, lwidth
      integer nxr, nzr, imdint, nint, ior
      real h, h2
      double precision fxl, fxw, ftox, foxdf, fsgap, fdgap
      character*(PM_MAX_PATHLEN) interface_filename
c
c     Workspace overlapped on double precision work store
      real wk(ARRAYSIZE*(ARRAYSIZE+2))
      equivalence(work_dbl, wk)
c
      write(luttyo,10)
10    format(' ','Reading doping file...',$)
c
      lfirst = .TRUE.
c
c     Append .2DOP if no extension
      call add_extension(binfilename,'.2dop')
c
c     If threed file already loaded, then read from it
      if( threed_loaded .and. binfilename .eq. prev_binfilename .and.
     1    binfiletype .eq. prev_binfiletype ) then
         write(luttyo,60)
60       format(' ','Info: Using 3D doping previously read...')
         call zero_real_array(arrdat,ARR2DSIZE*NUM_DATA_ARRAYS,1)
         nplanes = nz
         goto 2000
      endif
c
c     Attempt to read associated interface file (.INT) file
c     to get actual MOS gate length
      fxl = 0.0
      if( batchmode ) then
         interface_filename = tdop_interface_file
         if( tdop_interface_file .eq. ' ' .or.
     1       tdop_interface_file .eq. 'NOFILE' .or.
     2       tdop_interface_file .eq. 'nofile' ) then
            lexist = .FALSE.
            tdop_interface_file = ' '
         else
            lexist = .TRUE.
         endif
      else
         interface_filename = binfilename
         call replace_extension(interface_filename,'.int')
c        Inquire to see if .int file exists
         inquire(file=interface_filename,exist=lexist)
      endif
c
      if( lexist ) then
         nrecl=0
         call pm_open_file(binfil, interface_filename,
     1               'OLD','UNFORMATTED','READ','NONE','NONSHARED',
     2               nrecl, ioerr, ierr)
         if( ierr .ne. 0 ) goto 10060
         write(luttyo,30)
30       format(' ','Reading interface file to get default gate length')
         read(binfil,err=10070,end=10070) nint, fxl, fxw, ftox,
     1         foxdf, fsgap, fdgap
         close(unit=binfil)
      endif
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
c
c     1. Read NPIECE from binary file
      lwidth = .FALSE.
c        Try to read 3D form first. If there is a read error,
c        then rewind file and try 2D form.
      read(binfil,err=100,end=10010,iostat=ioerr) npiece, nplanes, ior
      three_d = .TRUE.
c     If ior is zero, this a 3DOP file
c     If ior is one, this a length direction 2DOP file
c     If ior is two, this a width direction 2DOP file
      if( ior .eq. 0 ) then
         three_d = .TRUE.
         lwidth = .FALSE.
      else if( ior .eq. 1 ) then
         three_d = .FALSE.
         lwidth = .FALSE.
      else if( ior .eq. 2 ) then
         three_d = .FALSE.
         lwidth = .TRUE.
      endif
c
      if( three_d ) then
c        Free any dynamic memory we've allocated, and
c        set base address of dynamic arrays to NULL value (0),
c        just to have them set to a fixed value. We will pass incore to
c        subroutine read_3d_plane_core, but the array isn't used
c        in 2d mode, so we should be OK.
         call pm_free_incore_arrays
      endif
c
      goto 110
c
100   rewind(binfil)
      read(binfil,err=10010,end=10010,iostat=ioerr) npiece
      three_d = .FALSE.
c
110   continue
c
      plane_cnt = 0
      if( three_d ) then
CC       INCORE
         if( threed_loaded ) then
            threed_loaded = .FALSE.
         endif
CC       END INCORE
      endif
c
c     Read loop for 3D doping
1000  continue
c
c     2. now nx, ny
      read(binfil,err=10010,end=10010,iostat=ioerr) nx,ny
      nz = 1
      nt = ny
      nyoxm = 0
      nytot= ny
      nxsp = 1
      nxsoxm = 1
c
c     Is array storage large enough?
      if( (MAX(nx,nytot,nz) .gt. COMP_ARRAYSIZE) .or.
     1    (MAX(nx*nytot,nx*nz,nytot*nz) .gt. ARR2DSIZE) .or.
     1    (three_d .and. MAX(nx,nytot,nz) .gt. CORESIZE) ) then
         write(luttyo,10050)
10050    format(' ','Error: binary file too large for POSTMINI arrays')
         close(unit=binfil)
         read_error = .TRUE.  ! go back
         return
      endif
c
c     For 3D case, allocate dynamic storage...
      if( three_d ) then
         call pm_allocate_incore_arrays(alloc_err)
         if( alloc_err ) then
            write(luttyo,10052)
10052       format(' ',
     1             'Error: could not allocate array for 3D binary file')
            read_error = .TRUE.  ! go back
         endif
      endif
c
      nxy=nx*ny
      nxytot=nx*nytot
c
c     3. read lateral spacing
      read(binfil,err=10010,end=10010,iostat=ioerr)
     1    (work_dbl(i),i=1,nx)
      do 200 i=1,nx
200      xyz(i,1) = SNGL(work_dbl(i))
c
c     4. read transversal spacing
      read(binfil,err=10010,end=10010,iostat=ioerr)
     1    (work_dbl(i),i=1,nytot)
      do 210 i=1,nytot
210      xyz(i,2) = SNGL(work_dbl(i))
c
c
c     Load generic mesh pointers, depending on plane type
c
      if( .NOT. three_d ) then
         if( lwidth ) then
c           Special case for width direction simulation
            plane_index = 1
            prev_plane_type = plane_type
            plane_type = YZ_PLANE
c           Load "X" coordinate array into "Z"
            do 211 i=1,nx
               xyz(i,3) = xyz(i,1)
211         continue
            xw = -2.0*xyz(1,3)
            nz = nx
            nx = 1
            nxy=ny
            nxytot=nytot
            xyz(1,1) = 0.0
         else
            plane_index = 1
            prev_plane_type = plane_type
            plane_type = XY_PLANE
            xw = 1.0
            xyz(1,3) = -xw/2.0
         endif
         call set_3d_index_limits
         cut_plane_axis = charb
      else
         ni = nx
      endif
      nx_lhs = ni
c
c     Since we create larger arrays that we actually read in
c     (and we don't know exactly how large they will be until they we
c     are finished processing the files), allocate the largest
c     arrays we can. Save the first three arrays for the final doping
c     (donors, acceptors, net).
c     Put the first set of donors/acceptors in the next two slots.
c     Put the second set of donors/acceptors (if any) in next two slots.
c     After creating the final doping in the first three arrays,
c     call allocate_rect_arrays to set arrind() to the final array sizes.
c
      call allocate_rect_arrays(ARR2DSIZE,1,NUM_DATA_ARRAYS,ierr)
      if( ierr .ne. 0 .or. NUM_DATA_ARRAYS .lt. 7 ) then
c        Failure during array allocation
         write(luttyo,10050)
         read_error = .TRUE.
         return
      endif
      call zero_real_array(arrdat,ARR2DSIZE*NUM_DATA_ARRAYS,1)
c
c     5. Load donors
      do j=1,ny
         read(binfil,end=10010,err=10010,iostat=ioerr)
     1       (work_dbl(i+(j-1)*ni),i=1,ni)
      end do
      do k=1,ni*ny
         arrdat(arrind(4)+k-1) = SNGL(work_dbl(k))
      end do
      in_binary_file(1) = .TRUE.
c
c     6. Skip over donor derivative arrays
      if( npiece .gt. 0 ) then
         do j=1,ny
            read(binfil,end=10010,err=10010,iostat=ioerr)
     1          (work_dbl(i+(j-1)*ni),i=1,ni)
         end do
         do j=1,ny
            read(binfil,end=10010,err=10010,iostat=ioerr)
     1          (work_dbl(i+(j-1)*ni),i=1,ni)
         end do
      endif
c
c     7. Load acceptors
      do j=1,ny
         read(binfil,end=10010,err=10010,iostat=ioerr)
     1       (work_dbl(i+(j-1)*ni),i=1,ni)
      end do
      do k=1,ni*ny
         arrdat(arrind(5)+k-1) =  SNGL(work_dbl(k))
      end do
      in_binary_file(2) = .TRUE.
c
      if( ABS(npiece) .eq. 1 ) then
         if( .NOT. lwidth ) then
c        Source/drain doping is symmetric. Mirror source doping to drain side.
c        Is array storage large enough?
         if( 2*nx .gt. COMP_ARRAYSIZE .or.
     1       2*nx*nytot .gt. ARR2DSIZE ) then
            write(luttyo,10051)
10051       format(' ','Warning: combined source/drain profiles are'/
     1             ' ',' too large for POSTMINI arrays. Part of drain '/
     2             ' ',' profile will be truncated.')
         endif
         nx_prev = nx
         nx = MIN(2*nx,COMP_ARRAYSIZE,(ARR2DSIZE/ny))
         nxy=nx*ny
         nxytot=nx*nytot
         if( .NOT. three_d ) call set_3d_index_limits
c
c        Compute offset needed to merge profiles to user specified
c        gate length.
         if( .not. batchmode .and. lfirst) then
            lfirst = .FALSE.
            write(luttyo,70)
70          format(' ','Doping can be scaled to any MOS gate length')
            if( fxl .ne. 0.0 ) then
               real_vec(1) = SNGL(fxl)*pm_scale_cm2plot
            else
               real_vec(1) = 2.0*xyz(nx_prev,1)*pm_scale_cm2plot
            endif
            format_array(1) = 'f9.3'
            call get_reals_with_default(real_vec, 1,
     1          'Enter gate length in microns: ', format_array, 0)
            tdop_gate_length = real_vec(1)*pm_scale_plot2cm
            if( tdop_gate_length .le. 0.0 ) tdop_gate_length = 0.5E-4
         endif
         offset = tdop_gate_length - 2.0*xyz(nx_prev,1)
         if( offset .le. 0.0 ) then
c           Profile extends too far into channel that a simple
c           mirroring of the profile will cause the channel to be
c           longer than the user specified gate length.
c           Reduce points in channel until we get a positive offset.
            i = nx_prev
            do while( offset .le. 0.0 .and. xyz(i,1) .gt. 0.0 )
               i = i - 1
               offset = tdop_gate_length - 2.0*xyz(i,1)
            end do
            nx_prev = i
            nx = MIN(2*nx_prev,COMP_ARRAYSIZE,(ARR2DSIZE/ny))
            nxy=nx*ny
            nxytot=nx*nytot
            if( .NOT. three_d ) call set_3d_index_limits
         endif
         j = nx_prev
         do i=nx_prev+1,nx
            xyz(i,1) = 2*xyz(nx_prev,1) + offset - xyz(j,1)
            j = j - 1
         end do
c        Copy donors to final storage location
         do j=1,ny
            ioff1 = (j-1)*nx
            ioff2 = arrind(4) - 1 + (j-1)*nx_lhs
            do i=1,nx_prev
               arrdat(ioff1+i) = arrdat(ioff2+i)
            end do
c           Mirror donors from source onto drain
            k = nx_prev
            do i=nx_prev+1,nx
               arrdat(ioff1+i) = arrdat(ioff2+k)
               k = k - 1
            end do
         end do
c        Copy acceptors to final storage location
         do j=1,ny
            ioff1 = nx*ny + (j-1)*nx
            ioff2 = arrind(5) - 1 + (j-1)*nx_lhs
            do i=1,nx_prev
               arrdat(ioff1+i) = arrdat(ioff2+i)
            end do
c           Mirror acceptors from source onto drain
            k = nx_prev
            do i=nx_prev+1,nx
               arrdat(ioff1+i) = arrdat(ioff2+k)
               k = k - 1
            end do
         end do
         else
c           Move "width" only doping file to final location
c           Remember: we have switched x and z, since this is a width
c           direction doping file. So use ni rather than nx.
c           Copy donors to final storage location
            do j=1,ny
               ioff1 = (j-1)*ni
               ioff2 = arrind(4) - 1 + (j-1)*nx_lhs
               do i=1,ni
                  arrdat(ioff1+i) = arrdat(ioff2+i)
               end do
            end do
c           Copy acceptors to final storage location
            do j=1,ny
               ioff1 = ni*ny + (j-1)*ni
               ioff2 = arrind(5) - 1 + (j-1)*nx_lhs
               do i=1,ni
                  arrdat(ioff1+i) = arrdat(ioff2+i)
              end do
            end do
         endif
      else
c
c     Read rest of 2DOP file
c
c        8. Skip over acceptor derivative arrays
         if( npiece .gt. 0 ) then
            do j=1,ny
               read(binfil,end=10010,err=10010,iostat=ioerr)
     1             (work_dbl(i),i=1,nx)
            end do
            do j=1,ny
               read(binfil,end=10010,err=10010,iostat=ioerr)
     1             (work_dbl(i),i=1,nx)
            end do
         endif
c
c        9. now nx_d, ny_d for drain side
         read(binfil,err=10010,iostat=ioerr) nx_d,ny_d
c
c        Is array storage large enough?
         if( MAX(nx+nx_d,nytot,nz) .gt. COMP_ARRAYSIZE .or.
     1       MAX(nx*nytot,nx*nz,nytot*nz) .gt. ARR2DSIZE) then
            write(luttyo,10051)
         endif
         nx_prev = nx
         nx_d_prev = nx_d
         nx = MIN(nx+nx_d,COMP_ARRAYSIZE,(ARR2DSIZE/ny))
         nxy=nx*ny
         nxytot=nx*nytot
         if( .NOT. three_d ) call set_3d_index_limits
c
c        10. read lateral spacing
         read(binfil,err=10010,end=10010,iostat=ioerr)
     1       (work_dbl(i),i=1,nx_d)
c        Ask for doping file gate length
         if( .not. batchmode .and. lfirst) then
            lfirst = .FALSE.
            write(luttyo,70)
            if( fxl .ne. 0.0 ) then
               real_vec(1) = SNGL(fxl)*pm_scale_cm2plot
            else
               real_vec(1) = (xyz(nx_prev,1)+work_dbl(nx_d))
     1            *pm_scale_cm2plot
            endif
            format_array(1) = 'f9.3'
            call get_reals_with_default(real_vec, 1,
     1          'Enter gate length in microns: ', format_array, 0)
            tdop_gate_length = real_vec(1)*pm_scale_plot2cm
            if( tdop_gate_length .le. 0.0 ) tdop_gate_length = 0.5E-4
         endif
c        Compute offset needed to merge profiles to user specified
c        gate length.
         dsource = xyz(nx_prev,1)
         ddrain = work_dbl(nx_d)
         offset = tdop_gate_length - (dsource + ddrain)
         if( offset .le. 0.0 ) then
c           Profiles extend too far into channel that we cannot append
c           the two profiles and still have the user specified gate length.
c           Reduce points in channel from each side until this is possible.
            j = 1
            do i=1,nx_prev
               if( xyz(i,1) .ge. 0.5*tdop_gate_length ) goto 215
               j = i
            enddo
215         continue
            nx_prev = j
c
            j = 1
            do i=1,nx_d
               if( work_dbl(i) .ge. 0.5*tdop_gate_length ) goto 216
               j = i
            enddo
216         continue
            nx_d_prev = nx_d
            nx_d = j
            offset = tdop_gate_length -
     1               (xyz(nx_prev,1) + work_dbl(nx_d))
            nx = nx_prev + nx_d
            nxy=nx*ny
            nxytot=nx*nytot
            if( .NOT. three_d ) call set_3d_index_limits
         endif
         j = nx_d
         do 220 i=nx_prev+1,nx
            xyz(i,1) = xyz(nx_prev,1) + offset + 
     1         SNGL(work_dbl(nx_d)) - SNGL(work_dbl(j))
            j = j - 1
220      continue
c
c        11. Read Y spacing. If different than left side, then give up!!
         if( ny .ne. ny_d ) goto 10020
         read(binfil,err=10010,end=10010,iostat=ioerr)
     1       (work_dbl(i),i=1,ny_d)
         do 230 i=1,ny
            if( xyz(i,2) .ne. SNGL(work_dbl(i)) ) goto 10020
230      continue
c
c        12. Load drain donors and append to source donors
c            Note: we load the drain donors "backwards" so the
c            channel sides of both profiles meet.
         do j=1,ny_d
            read(binfil,end=10010,err=10010,iostat=ioerr)
     1          (work_dbl(i+(j-1)*nx_d_prev),i=1,nx_d_prev)
         end do
         do j=1,ny
c           Copy LHS donors to final storage location
            ioff1 = (j-1)*nx
            ioff2 = arrind(4) - 1 + (j-1)*nx_lhs
            do i=1,nx_prev
               arrdat(ioff1+i) = arrdat(ioff2+i)
            end do
c           Copy RHS donors to final storage location
            k = nx_d
            ioff1 = (j-1)*nx
            ioff2 = (j-1)*nx_d_prev
            do i=nx_prev+1,nx
               arrdat(ioff1+i) =  SNGL(work_dbl(ioff2+k))
               k = k - 1
            end do
         end do
c
c        12. Skip over (drain) donor derivative arrays
         if( npiece .gt. 0 ) then
            do j=1,ny_d
               read(binfil,end=10010,err=10010,iostat=ioerr)
     1             (work_dbl(i),i=1,nx_d_prev)
            end do
            do j=1,ny_d
               read(binfil,end=10010,err=10010,iostat=ioerr)
     1             (work_dbl(i),i=1,nx_d_prev)
            end do
         endif
c
c        13. Load drain acceptors with the source acceptors
c            Note: we load the drain acceptors "backwards" so the
c            channel sides of both profiles meet.
         do j=1,ny_d
            read(binfil,end=10010,err=10010,iostat=ioerr)
     1          (work_dbl(i+(j-1)*nx_d_prev),i=1,nx_d_prev)
         end do
         do j=1,ny
c           Copy LHS acceptors to final storage location
            ioff1 = nx*ny + (j-1)*nx
            ioff2 = arrind(5) - 1 + (j-1)*nx_lhs
            do i=1,nx_prev
               arrdat(ioff1+i) = arrdat(ioff2+i)
            end do
c           Copy RHS acceptors to final storage location
            k = nx_d
            ioff1 = nx*ny + (j-1)*nx
            ioff2 = (j-1)*nx_d_prev
            do i=nx_prev+1,nx
               arrdat(ioff1+i) =  SNGL(work_dbl(ioff2+k))
               k = k - 1
            end do
         end do
      endif
c
c     If 3D doping file, write donors and acceptors to scratch file
c     and loop to read next plane.
      if( three_d ) then
         plane_cnt = plane_cnt + 1
         read(binfil,end=10010,err=10010,iostat=ioerr) xyz(plane_cnt,3)
c        Write donors and acceptors to scratch file
c        so we can later select them by an arbitrary cut plane.
CC       INCORE
         do j=1,ny
            ioff1 = (j-1)*nx
            ioff2 = nx*ny + (j-1)*nx
            do i=1,nx
               incore(i,j,plane_cnt,1) = arrdat(ioff1+i)
               incore(i,j,plane_cnt,2) = arrdat(ioff2+i)
            end do
         end do
CC       END INCORE
         if( plane_cnt .lt. nplanes ) goto 1000
      endif
c
      close(unit=binfil)
      if( three_d ) then
         threed_loaded = .TRUE.
      endif
2000  continue
      if( three_d )  then
         nz = nplanes
         prev_plane_type = plane_type
         if( batchmode ) then
            call batch_select_plane(cut_plane_axis,
     1           plane_index, plane_value,
     2           plane_type, plane_actual_value, xyz, nx, ny, nz,
     3           plane_select_err)
            if( plane_select_err ) then
               read_error = .TRUE.
               close(unit=binfil)
               return
            endif
         else
            call select_plane(plane_type, plane_index, plane_value,
     1                        plane_actual_value, xyz, nx, ny, nz)
         endif
c        Load generic mesh pointers, depending on plane type
         call set_3d_index_limits
         cut_plane_axis = charb
c        Allocate arrays to 'final' sizes...
         call allocate_rect_arrays(ni,nj,NUM_DATA_ARRAYS,ierr)
CC       INCORE
         call read_3d_plane_incore(incore(1,1,1,1), arrdat(arrind(1)),
     1                      nx, ny, nz, 0,
     2                      plane_index, plane_type, ni,
     3                      CORESIZE, ierr)
         in_binary_file(1) = .TRUE.
        call read_3d_plane_incore(incore(1,1,1,2), arrdat(arrind(2)),
     1                      nx, ny, nz, 0,
     2                      plane_index, plane_type, ni,
     3                      CORESIZE, ierr)
CC       END INCORE
         in_binary_file(2) = .TRUE.
         if( ierr .ne. 0 ) goto 10010
      else
c        Now, allocate arrays to 'final' sizes...
         call allocate_rect_arrays(ni,nj,NUM_DATA_ARRAYS,ierr)
      endif
c
c     Compute net doping
      call pm_diff_arrays(arrdat(arrind(3)),
     1                    arrdat(arrind(1)),arrdat(arrind(2)),ni,nj,ni)
      in_binary_file(3) = .TRUE.
c     This variable must be set by all read routines...
      data_is_trigrid = .FALSE.
c
      write(luttyo,120)
120   format(' ','...done')
c
c     Read associated interface file (.INT) file
      if( lexist ) then
         nrecl=0
         call pm_open_file(binfil, interface_filename,
     1               'OLD','UNFORMATTED','READ','NONE','NONSHARED',
     2               nrecl, ioerr, ierr)
         if( ierr .ne. 0 ) goto 10060
         write(luttyo,306)
306      format(' ','Reading associated interface file',$)
         read(binfil,err=10070,end=10070) nint, fxl, fxw, ftox,
     1         foxdf, fsgap, fdgap
         do 300 k=1,2
            read(binfil,err=10070,end=10070) nxr, nzr, imdint
            if( MAX(nxr,nzr) .le. ARRAYSIZE ) then
               ioff1 = nxr
               read(binfil,err=10070,end=10070)
     1              (work_dbl(i),i=1,nxr), (work_dbl(ioff1+j),j=1,nzr)
               do 303 i=1,nxr
                  fox_x(i) = work_dbl(i)
303            continue
               do 304 i=1,nzr
                  fox_z(i) = work_dbl(ioff1+i)
304            continue
               read(binfil,err=10070,end=10070)
     1             (work_dbl(i),i=1,nxr*nzr)
               do 307 i=1,nxr*nzr
                   work(i) = SNGL(work_dbl(i))
307            continue
c              Scale to gate length
               if( .NOT. lwidth ) then
                  call modint(work,fox_x,fox_z,nxr,nzr,
     1                        tdop_gate_length, SNGL(fxl))
               endif
c
c              Negate interface information so it plots "upright"
c              on X-Y graphs
               do 305 i=1,nxr*nzr
                  work(i) = - pm_scale_cm2plot*work(i)
305            continue
c              Interpolate onto 2DOP X/Z mesh
c              work_dbl is partitioned into three temporary work spaces
               if( lwidth ) then
                  call int2dm(work, fox_x, nxr, fox_z, nzr,
     1                        xyz(1,3), nz, xyz(1,1), nx,
     2                        wk, wk(ARRAYSIZE*ARRAYSIZE+1),
     3                        wk(ARRAYSIZE*(ARRAYSIZE+1)+1))
               else
                  call int2dm(work, fox_x, nxr, fox_z, nzr,
     1                        xyz(1,1), nx, xyz(1,3), nz,
     2                        wk, wk(ARRAYSIZE*ARRAYSIZE+1),
     3                        wk(ARRAYSIZE*(ARRAYSIZE+1)+1))
               endif
c              Place 1D cuts (for drawing boundaries) in 1d store
               if( plane_type .eq. XY_PLANE ) then
                  do i=1,nx
                     arr2d(i,k) = work(i + (plane_index-1)*nx)
                  end do
               else if( plane_type .eq. YZ_PLANE ) then
                  do i=1,nz
                     arr2d(i,k) = work(plane_index + (i-1)*nx)
                  end do
               else if( plane_type .eq. XZ_PLANE ) then
                  if( k .eq. 2 ) then
c                    Find approximate edge of lower oxide shape
                     do i=1,nx
                        h = xyz(nz,3)
c                       Remember, values are scaled by -pm_scale_cm2plot
                        h2 = - plane_actual_value*pm_scale_cm2plot
                        do j=1,nz-1
                          if( work(i+(j-1)*nx) .gt. h2 .and.
     1                        h2 .ge. work(i+j*nx) )
     2                    h = xyz(j,3) + (xyz(j+1,3)-xyz(j,3))*
     3                        (h2-work(i+(j-1)*nx))/
     4                        (work(i+j*nx)-work(i+(j-1)*nx))
                       enddo
                       arr2d(i,k) = - pm_scale_cm2plot*h
                     end do
                  else
c                    The "upper" shape is assigned a value to make the
c                    fill area work correctly.
                     do i=1,nx
                       arr2d(i,k) = - pm_scale_cm2plot*xyz(nz,3)
                     end do
                  endif
               endif
c              ASSUMES there are three entries preceeding the interface
c              information.
               in_binary_file(3+k) = .TRUE.
            else
               write(luttyo,302)
302            format(' ','Error reading interface file'/
     1                ' ','Data too large for internal arrays')
c              Break out of read loop
               goto 301
            endif
300      continue
301      continue
c
         close(binfil)
         if( in_binary_file(4) .and. in_binary_file(5) ) then
            lint = .TRUE.
            write(luttyo,120)
            tdop_interface_file = interface_filename
         endif
      else
         tdop_interface_file = ' '
         lint = .FALSE.
      endif
c
      read_error = .FALSE.
      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'/
     2       ' ','   - 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 2DOP 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 SUPMIN3_2D'/)
      endif
      read_error = .TRUE.  ! go back and get new name...
      close(unit=binfil)
      goto 10100
10020 write(luttyo,10021)
10021 format(' ','Error: Y mesh for drain profile is different'/
     1       ' ','       than Y mesh for source profile')
      read_error = .TRUE.  ! go back and get new name...
      close(unit=binfil)
      goto 10100
c
10060 call expand_filename(interface_filename,expanded_name)
      ilen = string_len(expanded_name)
      write(luttyo,10061) ioerr, expanded_name(1:ilen)
10061 format(/' ','FORTRAN I/O error #',i5,' on open of interface ',
     1       'file:'/' ',a)
      read_error = .TRUE.
      goto 10100
c
10070 call expand_filename(interface_filename,expanded_name)
      ilen = string_len(expanded_name)
       write(luttyo,10071) ioerr, expanded_name(1:ilen)
10071 format(/' ','FORTRAN I/O error #',i5,' on read of interface ',
     1       'file:'/' ',a/
     2       ' ','Interface information will not be available')
      lint = .FALSE.
      close(binfil)
      goto 10100
c
10100 continue
      return
      end
c
      subroutine load_2dop_specific(lint)
      implicit none
c
c     Load descriptive arrays with MINIMOS 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
      logical lint
c
c     Number of menu items
      if( lint ) then
         max_menu_num = 5
      else
         max_menu_num = 3
      endif
c     Array index of net doping
      net_dop_array_index = 3
c
      call load_desc_arr(1, 'Donor Concentration', 'DONORS',
     1                   'cm^{-3}','cm^{-2}','cm^{-1}','.don',
     2                   LOG,NO_OX,NO_JUMP,TWOD,1)
c
      call load_desc_arr(2, 'Acceptor Concentration', 'ACCEPTORS',
     1                   'cm^{-3}','cm^{-2}','cm^{-1}','.acc',
     2                   LOG,NO_OX,NO_JUMP,TWOD,2)
c
      call load_desc_arr(3, 'Net Concentration', 'NET_DOPING',
     1                   'cm^{-3}','cm^{-2}','cm^{-1}','.net',
     2                   LOG,NO_OX,NO_JUMP,TWOD,3)
      if( lint ) then
         call load_desc_arr(4, 'Upper oxide shape', 'UPPER_OXIDE',
     1                   'microns','microns^{2}','microns^{3}','.foxo',
     2                   LIN,OX,NO_JUMP,ONED,1)
         call load_desc_arr(5, 'Lower oxide shape', 'LOWER_OXIDE',
     1                   'microns','microns^{2}','microns^{3}','.foxu',
     2                   LIN,OX,NO_JUMP,ONED,2)
      endif
      return
      end
c
      subroutine int1dm(z1, zmesh1, npt1, zmesh2, npt2, dest)
      implicit none 
c
c     Linearly interpolates data on two 1d meshes.
c     Use for interpolating oxide information, which varies on a linear
c     scale.
c
c     On input:
c     z1     -- 1d data values to be interpolated
c     zmesh1 -- 1d mesh coordinates of old grid
c     npt1   -- number of points in old grid
c     zmesh2 -- 1d mesh coordinate of new grid
c     npt2   -- number of points in new grid
c
c     On output:
c     dest   -- 1d data values interpolated onto new grid
c
      integer npt1, npt2, max
      real    z1(max(npt1,npt2)), zmesh1(npt1), zmesh2(npt2),
     1          dest(npt2)
c
c     locals
      integer i, j
c
      do 100 j=1,npt2
         if( zmesh2(j) .le. zmesh1(1) ) then
            dest(j) = z1(1)
         else if( zmesh2(j) .ge. zmesh1(npt1) ) then
            dest(j) = z1(npt1)
         else
            do 200 i=1,npt1-1
               if( zmesh2(j) .ge. zmesh1(i) .and.
     1             zmesh2(j) .le. zmesh1(i+1) ) then
                   dest(j)=z1(i)+(z1(i+1)-z1(i))*
     1               (zmesh2(j)-zmesh1(i))/(zmesh1(i+1)-zmesh1(i))
                  goto 150
               endif
200         continue
         endif
150      continue
100   continue
c
      end
c
      subroutine int2dm(z, imesh1, ipt1, jmesh1, jpt1, imesh2, ipt2,
     1                  jmesh2, jpt2, work, temp1, temp2)
      implicit none 
c
c     Linearly interpolates data on two 2d meshes.
c     Use for interpolating oxide information, which varies on a linear
c     scale.
c
c     On input:
c     z      -- 2d data values to be interpolated
c     imesh1 -- 1d mesh coordinates of old grid
c     ipt1   -- number of points in old grid
c     jmesh1 -- 1d mesh coordinates of old grid
c     npt1   -- number of points in old grid
c     imesh2 -- 1d mesh coordinate of new grid
c     jpt2   -- number of points in old grid
c     imesh2 -- 1d mesh coordinate of new grid
c     jpt2   -- number of points in new grid
c
c     On output:
c     z      -- 2d data values interpolated onto new grid
c
c     work   -- 2d working store (must be as large as initial/final z)
c
      integer ipt1, jpt1, ipt2, jpt2, max
      real    z(max(ipt1*jpt1,ipt2*jpt2)), imesh1(ipt1), jmesh1(jpt1),
     1        imesh2(ipt2), jmesh2(jpt2), work(MAX(ipt1*jpt1,ipt2*jpt2))
      real temp1(MAX(jpt1,jpt2)), temp2(MAX(jpt1,jpt2))
c
c     A really disgusting 2D interpolation scheme. First intepolate
c     in 1D strips in one direction, then in the other.
c     Note: we've modified int1dm, so it places its output in the
c     last argument (work)
c
      integer i,j
c
      do 100 j=1,jpt1
          call int1dm(z((j-1)*ipt1+1), imesh1, ipt1, imesh2, ipt2,
     1                work((j-1)*ipt2+1) )
100   continue
c
      do 200 i=1,ipt2
          do 210 j=1,jpt1
             temp1(j) = work((j-1)*ipt2 + i)
210       continue
          call int1dm(temp1, jmesh1, jpt1, jmesh2, jpt2, temp2)
          do 220 j=1,jpt2
             z((j-1)*ipt2 + i) = temp2(j)
220       continue
200   continue
c
      return
      end
c
      subroutine modint(oxida,xox,zox,nxr,nzr,xlgeo,xlint)
      implicit none
c
      include 'sysdep.cmn'
      include 'io.cmn'
c
c     modint - modifies interface information depending on gate length
c
      integer nxr, nzr
      real xlgeo,xlint
      real oxida(nxr*nzr), xox(nxr), zox(nzr)
c
      integer i, j, k, igatel, igatem, igater, nxrp
      real tmp, delta, h1, h2
c
c     find edges of gate. note that x points in interface file do not
c     have to lie exactly on the gate edges.
      igatel=0
      igatem=0
      igater=0
      tmp = xlint/2
      do 1000 i=1,nxr-1
         if( xox(i) .lt. 0   .and. xox(i+1) .ge. 0         ) igatel=i+1
         if( xox(i) .lt. tmp .and. xox(i+1) .ge. tmp       ) igatem=i+1
         if( xox(i) .le. xlint  .and. xox(i+1) .gt. xlint  ) igater=i
1000  continue
      if( igatel .eq. 0 .or. igatem .eq. 0 .or. igater .eq. 0 ) then
         write(luttyo,1020)
1020     format(' ','Error in conversion of interface file'/
     1          ' ','could not find gate edges')
         return
      endif
c
      if( xlgeo .ge. xlint ) then
c       expand gate at middle
        delta = xlgeo - xlint
        do 2000 i=igatem,nxr
           xox(i) = xox(i) + delta
2000    continue
      else
c        subtract off half of gate length difference on each side
         delta = xlint - xlgeo
         h1    = (xlint - delta)/2
         h2    = (xlint + delta)/2
c        determine new number of x points
         k=1
         do 3000 i=1,nxr
            if( xox(i) .le. h1 .or. xox(i) .ge. h2 ) k=k+1
3000     continue
         nxrp = k-1
c        adjust interface
c        note: we can adjust "in place" because the "i" index
c        is always ahead of the "k" index. doing it in place is
c        not very efficient, but we only do this once.
         do 3010 j=1,nzr
            k=1
            do 3020 i=1,nxr
               if( xox(i) .le. h1 ) then
                  oxida(k+nxrp*(j-1)) = oxida(i+nxr*(j-1))
                  k=k+1
               else if ( xox(i) .ge. h2 ) then
                  oxida(k+nxrp*(j-1)) = oxida(i+nxr*(j-1))
                  k=k+1
               endif
3020        continue
3010     continue
c        now adjust mesh
         k=1
         do 3030 i=1,nxr
            if( xox(i) .le. h1 ) then
               k=k+1
            else if ( xox(i) .ge. h2 ) then
               xox(k) = xox(i) - delta
               k=k+1
            endif
3030     continue
         nxr = nxrp
      endif
c
      return
      end
c
      subroutine read_3d_plane_incore(source, dest,
     1                   nx, ny, nz, y_offset,
     2                   plane_index, plane_type, isize, icore,
     3                   ierr)
      implicit none
c
      include 'const.h'
c
      integer nx, ny, nz, y_offset,
     1        plane_type, plane_index,
     2        isize, icore, ierr
      real    source(icore,icore,icore)
      real    dest(isize,1)
c
      integer i, iz, k
c
      ierr = 0
c
      if(plane_type .eq. XY_PLANE) then
         do k=1,ny
            do i=1,nx
               dest(i,k+y_offset) = source(i,k,plane_index)
            end do
         end do
      else if (plane_type .eq. YZ_PLANE) then
         do iz=1,nz
            do k=1,ny
               dest(iz,k+y_offset) = source(plane_index,k,iz)
            end do
         end do
      else if (plane_type .eq. XZ_PLANE) then
         do iz=1,nz
            do k=1,nx
               dest(k,iz) = source(k,plane_index,iz)
            end do
         end do
      else
         ierr = 1
      endif
c
      return
      end
c
      subroutine pm_diff_arrays(dest,a,b,nx,ny,isize)
      implicit none
c
      integer nx,ny,isize
      real dest(isize,1), a(isize,1), b(isize,1)
c
c     locals 
      integer i,j
c
      do j=1,ny
         do i=1,nx
            dest(i,j) = a(i,j) - b(i,j)
         end do
      end do
c
      return
      end

