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 pisces mesh and solution files
c
c
      subroutine read_pisces(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 PISCES binary mesh and solution files
c
      logical batchmode, read_error
c
c     locals
c
c
      call read_pisces_binary(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_pisces_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_pisces_binary(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 PISCES binary mesh and solution file
c
      logical batchmode, read_error
c
c     locals
c
      integer i, j, ilen, ioerr, string_len,
     1        nb, nelect, nepb, local_record, ierr, nrecl
      real simtime, bias(25), amp(25)
      real rd1, rd2, rd3, rd4, rd5, rd6, rda1(3), rda2(3), rda3(3)
      integer id1, id2, id3, id4, id5
      logical lsolution_file
      character*(PM_MAX_PATHLEN) expanded_name, datafile
      character*60 ititle
      character*1 ch
c
      double precision dda(7)
      double precision, allocatable :: fv(:), fn(:), fp(:),
     1                                 qfn(:), qfp(:)
c
      write(luttyo,10)
10    format(' ','Reading PISCES mesh file...')
c
c     Open PISCES postprocessing file
      nrecl=0
      datafile = binfilename
      call pm_open_file(binfil, datafile,
     1               'OLD','UNFORMATTED','READ','NONE','NONSHARED',
     2               nrecl, ioerr, ierr)
      if( ierr .ne. 0 ) goto 10000
      local_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)
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
      tri_npts = 0
      tri_ne   = 0
      tri_nreg = 0
c
c     Read header
      read(binfil,end=10010,err=10010,iostat=ioerr) ititle
c
      read(binfil,end=10010,err=10010,iostat=ioerr)
     1    nx, ny,
     1    tri_npts,tri_ne,nb,nelect,tri_nreg,nepb
      if( tri_npts .gt. ARR2DSIZE ) goto 10040
      if( tri_ne   .gt. MAX_TRI ) goto 10080
      if( tri_nreg .gt. MAX_REGIONS ) goto 10050
c     
c     Coordinates, doping
      read(binfil,end=10010,err=10010,iostat=ioerr)
     1             ((cord(j,i),j=1,2), arrdat(i+arrind(1)-1),
     2             rd1, id1, rd2, rd3, rd4, id2, rd5, rd6,i=1,tri_npts)
      in_binary_file(1) = .TRUE.
c
c     Triangle node numbers, adjacency information, region number
      read(binfil,end=10010,err=10010,iostat=ioerr)
     1             ((nop(j,i),j=1,3), (rda1(j),j=1,3),
     2              (rda2(j),j=1,3), (rda3(j),j=1,3),
     3              (tri_nextel(j,i),j=1,3),
     4              tri2reg(i),i=1,tri_ne)
c
c     SKIP FOR NOW: electrode information
      read(binfil,end=10010,err=10010,iostat=ioerr)
      read(binfil,end=10010,err=10010,iostat=ioerr)
c
      read(binfil,end=10010,err=10010,iostat=ioerr)
     1             (imat(i),i=1,tri_nreg)
c     SKIP FOR NOW: timestamp
      read(binfil,end=10010,err=10010,iostat=ioerr)
c
      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
      close(unit=binfil)
c
      if( batchmode ) then
         if( medici_solution_file .eq. 'NOFILE' .or.
     1       medici_solution_file .eq. 'nofile' .or.
     1       medici_solution_file .eq. ' ' ) then
            lsolution_file = .FALSE.
         else
            lsolution_file = .TRUE.            
         endif
      else
         call get_str_with_default(medici_solution_file,
     1   'Pisces solution file? (enter NOFILE to just read mesh file)',
     2   1)
         datafile = medici_solution_file
         call upcase(datafile)
         if( datafile(1:6) .eq. 'NOFILE' .or. datafile .eq. ' ' ) then
            lsolution_file = .FALSE.
         else
            lsolution_file = .TRUE.
         endif
      endif
      if( .NOT. lsolution_file ) then
         read_error = .FALSE.
         datafile_record = local_record
         goto 10200
      endif
c
      write(luttyo,40)
40    format(' ','Reading PISCES data file...')
c
      nrecl=0
      datafile = medici_solution_file
      call pm_open_file(binfil, datafile,
     1               'OLD','UNFORMATTED','READ','NONE','NONSHARED',
     2               nrecl, ioerr, ierr)
      if( ierr .ne. 0 ) goto 10000
c
c     get stem of file name
      call get_filestem(binfilename, binfilestem)
c
c     SKIP FOR NOW: timestamp
      read(binfil,end=10010,err=10010,iostat=ioerr)
c
c     Np, ne, nelect, nmat, rtime, temp, delt, delt0, ncarr0
      read(binfil,end=10010,err=10010,iostat=ioerr)
     1     id1, id2, id3, id4, simtime, temperature, rd1, rd2, id5
c
c     Potential, n, p, qfn, qfp
c     Note: these are double precision, and due to the way they are
c     written by PISCES, we need to allocate a temporary array to read
c     them in.
      allocate(fv(1:tri_npts),fn(1:tri_npts),fp(1:tri_npts),
     1         qfn(1:tri_npts),qfp(1:tri_npts),stat=ierr)
      if( ierr .ne. 0 ) then
c        Failure during array allocation
         write(luttyo,10091)
10091    format(' ','Error: Dynamic allocation of memory ',
     1              '(malloc) failed')
         read_error = .TRUE.
         close(unit=binfil)
         return
      endif
      read(binfil,end=10010,err=10010,iostat=ioerr)
     1    (fv(i), fn(i), fp(i), qfn(i), qfp(i), i=1,tri_npts)
      do i=1,tri_npts
         arrdat(i+arrind(2)-1) = fv(i)
         arrdat(i+arrind(3)-1) = fn(i)
         arrdat(i+arrind(4)-1) = fp(i)
         arrdat(i+arrind(5)-1) = qfn(i)
         arrdat(i+arrind(6)-1) = qfp(i)
      end do
      deallocate(fv,fn,fp,qfn,qfp)
      in_binary_file(2) = .TRUE.
      in_binary_file(3) = .TRUE.
      in_binary_file(4) = .TRUE.
      in_binary_file(5) = .TRUE.
      in_binary_file(6) = .TRUE.
      tri_ndata = 6
c
c     Bias conditions
      read(binfil,end=10010,err=10010,iostat=ioerr)
     1    (bias(i), amp(i), rd1, rd2, rd3, rd4, i=1,nelect)
c     SKIP FOR NOW: dfxpt, dfxpt0
      read(binfil,end=10010,err=10010,iostat=ioerr)
c     SKIP FOR NOW: vsn,vsp,workf,cresis,ccapac (i=1,nelect)
      read(binfil,end=10010,err=10010,iostat=ioerr)
c     SKIP FOR NOW: nintrf,((ssdata(i,j),j=1,7),i=1,nintrf)
      read(binfil,end=10010,err=10010,iostat=ioerr)
c     SKIP for NOW: lconmb,lfldmb,lconlt,lsrh,lauger,lbgn,lboltz,
c                   lincom,limpct,labbas,lcyl, lphgen, abscof,
c                   flux*SNGL(dcscl)*ktq
      read(binfil,end=10010,err=10010,iostat=ioerr)
c     SKIP for NOW: (mattyp(i),epsmat(i),i=1,nmat)
      read(binfil,end=10010,err=10010,iostat=ioerr)
c     SKIP for NOW: semiconductor parameters
      read(binfil,end=10010,err=10010,iostat=ioerr)
c
c     TRY to read values included when 'currents' is specified on solve
c     statement...
      read(binfil,end=109,err=10010,iostat=ioerr) (dda(j),j=1,7)
      do i=7,13
         arrdat(arrind(i)) = dda(i-6)
         in_binary_file(i) = .TRUE.
      end do
      tri_ndata = 13
      do i=2,tri_npts
         read(binfil,end=10010,err=10010,iostat=ioerr)(dda(j),j=1,7)
         do j=7,13
            arrdat(i+arrind(j)-1) = dda(j-6)
         end do
      end do
c
109   write(luttyo,110)
110   format(' ','...done')
      datafile_record = local_record
      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(datafile,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 10200
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')
      read_error = .TRUE.
      close(unit=binfil)
      goto 10200
c
10020 continue
      write(luttyo,10021)
10021 format(' ','Error: Invalid header in PISCES postprocessing file')
      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
10080 continue
      write(luttyo,10081)
10081 format(' ','Error: Number of triangles exceeds internal storage')
      read_error = .TRUE.
      close(unit=binfil)
      goto 10200
c
10200 continue
      return
      end
c
      subroutine load_pisces_specific
      implicit none
c
c     Load descriptive arrays with PISCES 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
c
c     Number of menu items
      max_menu_num = tri_ndata
c
c     (2D) array of net doping
      net_dop_array_index = 1
c     File has only 2D data
      three_d = .FALSE.
c
      call load_desc_arr(1, 'Doping Concentration', 'NET_DOPING',
     1                   'cm^{-3}','cm^{-2}','cm^{-1}','.dconc',
     2                   LOG,NO_OX,NO_JUMP,TWOD,1)
c
      call load_desc_arr(2, 'Potential', 'POTENTIAL',
     1                   'V','V cm','V cm^{2}','.psi',
     2                   LIN,OX,NO_JUMP,TWOD,2)
c
      call load_desc_arr(3, 'Electron Concentration', 'ELECTRONS',
     1                   'cm^{-3}','cm^{-2}','cm^{-1}','.econc',
     2                   LOG,NO_OX,NO_JUMP,TWOD,3)
c
      call load_desc_arr(4, 'Hole Concentration', 'HOLES',
     1                   'cm^{-3}','cm^{-2}','cm^{-1}','.hconc',
     2                   LOG,NO_OX,NO_JUMP,TWOD,4)
c
      call load_desc_arr(5, 'Electron qFermi level', 'QFN',
     1                   'V','V cm','V cm^{2}','.qf_e',
     2                   LIN,NO_OX,NO_JUMP,TWOD,5)
c
      call load_desc_arr(6, 'Hole qFermi level', 'QFP',
     1                   'V','V cm','V cm^{2}','.qf_h',
     2                   LIN,NO_OX,NO_JUMP,TWOD,6)
c
      if( tri_ndata .gt. 6 ) then
         call load_desc_arr(7,'Conduction Current','J_CONDUCTION',
     1                     'A/cm^{2}','A/cm','A','.j_cond',
     2                      LOG,NO_OX,NO_JUMP,TWOD,7)
c
         call load_desc_arr(8,'Electron Current','J_ELEC',
     1                     'A/cm^{2}','A/cm','A','.j_elec',
     2                      LOG,NO_OX,NO_JUMP,TWOD,8)
c
         call load_desc_arr(9,'Hole Current','J_HOLE',
     1                     'A/cm^{2}','A/cm','A','.j_hole',
     2                      LOG,NO_OX,NO_JUMP,TWOD,9)
c
         call load_desc_arr(10,'Displacement Current','J_DISP',
     1                     'A/cm^{2}','A/cm','A','.j_disp',
     2                      LOG,NO_OX,NO_JUMP,TWOD,10)
c
         call load_desc_arr(11,'Total Current','J_TOTAL',
     1                     'A/cm^{2}','A/cm','A','.j_tot',
     2                      LOG,NO_OX,NO_JUMP,TWOD,11)
c
         call load_desc_arr(12,'Electric field','ELEC_FIELD',
     1                     'V/cm','V','V cm','.efld',
     2                      LOG,NO_OX,NO_JUMP,TWOD,12)
c
         call load_desc_arr(13,'Net Recombination','NET_RECOMB',
     1                      'cm^{-3} s^{-1}','cm^{-2} s^{-1})',
     2                      'cm^{-1} s^{-1})','.net_recomb',
     3                      LOG,NO_OX,NO_JUMP,TWOD,13)
      endif
c
      return
      end
