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 readbinfile
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
      include 'io.cmn'
      include 'plot.cmn'
      include 'user.cmn'
      include 'callable.cmn'
      include 'analytical.cmn'
c
      logical read_error
      character*(PM_MAX_PATHLEN) savefile
      character*10 file_type
      character*1 answer
      integer type, rc
c
      if( binfiletype .eq. DT_MINIMOS ) then
         file_type = 'MINIMOS'
      else if( binfiletype .eq. DT_PROMIS ) then
         file_type = 'PROMIS'
      else if( binfiletype .eq. DT_SUPREM3S ) then
         file_type = 'SUPREM3S'
      else if( binfiletype .eq. DT_SUPREM3P ) then
         file_type = 'SUPREM3P'
      else if( binfiletype .eq. DT_VLSICAP ) then
         file_type = 'VLSICAP'
      else if( binfiletype .eq. DT_2DOP ) then
         file_type = '2DOP'
      else if( binfiletype .eq. DT_USEOUT ) then
         file_type = 'USEOUT'
      else if( binfiletype .eq. DT_BAMBI ) then
         file_type = 'BAMBI'
      else if( binfiletype .eq. DT_CALLABLE ) then
         file_type = 'RELOAD'
      else if( binfiletype .eq. DT_SUPREM4 ) then
         file_type = 'SUPREM4'
      else if( binfiletype .eq. DT_PISCES ) then
         file_type = 'PISCES'
      else if( binfiletype .eq. DT_ASCII ) then
         file_type = 'ASCII'
      else if( binfiletype .eq. DT_TIF ) then
         file_type = 'TIF'
      else if( binfiletype .eq. DT_FAIM ) then
         file_type = 'FAIM'
      else if( binfiletype .eq. DT_MCP2D ) then
         file_type = 'MCP2D'
      else if( binfiletype .eq. DT_MEDICI ) then
         file_type = 'MEDICI'
      else if( binfiletype .eq. DT_GRAPES ) then
         file_type = 'GRAPES'
      else if( binfiletype .eq. DT_SPICE3 ) then
         file_type = 'SPICE3'
      else if( binfiletype .eq. DT_ANALYTICAL ) then
         file_type = 'ANALYTICAL'
      else
         file_type = 'PMI'
      endif
c
100   continue
      write(luttyo,110)
110   FORMAT(' ','Data file types:'/
     1      ' ','ASCII (also SPICE output file), PMI, ANALYTICAL,'/
     1      ' ','SPICE3,'/
     1      ' ','MINIMOS, PISCES, TIF, 2DOP, USEOUT,'/
     1      ' ','SUPREM3S, SUPREM3P, SUPREM4, PROMIS,'/
     1      ' ','VLSICAP, FAIM, MCP2D, BAMBI')
      if( callable_loaded ) then
         write(luttyo,111)
111      format(' ',
     1   ' RELOAD -- reload data from application (e.g. VIDE)')
      endif
      call get_str_with_default(file_type,
     1   'Enter data file type: ',0)
      call upcase(file_type)
      if( file_type(1:1) .eq. 'M' .and.
     1    file_type(1:2) .ne. 'MC' .and.
     1    file_type(1:2) .ne. 'ME' ) then
         file_type = 'MINIMOS'
         prev_binfiletype = binfiletype
         binfiletype = DT_MINIMOS
      else if( file_type(1:1) .eq. 'P' ) then
         if( file_type(2:2) .eq. 'M' ) then
            file_type = 'PMI'
         else if( file_type(2:2) .eq. 'R' ) then
            file_type = 'PROMIS'
            prev_binfiletype = binfiletype
            binfiletype = DT_PROMIS
         else if( file_type(2:2) .eq. 'I' ) then
            file_type = 'PISCES'
            prev_binfiletype = binfiletype
            binfiletype = DT_PISCES
         else
            file_type = 'PROMIS'
            prev_binfiletype = binfiletype
            binfiletype = DT_PROMIS
         endif
      else if( file_type(1:1) .eq. 'V' ) then
         file_type = 'VLSICAP'
         prev_binfiletype = binfiletype
         binfiletype = DT_VLSICAP
      else if( file_type(1:1) .eq. 'B' ) then
         file_type = 'BAMBI'
         prev_binfiletype = binfiletype
         binfiletype = DT_BAMBI
      else if( file_type(1:1) .eq. 'U' ) then
         file_type = 'USEOUT'
         prev_binfiletype = binfiletype
         binfiletype = DT_USEOUT
      else if( file_type(1:1) .eq. 'S' ) then
         if( file_type(2:2) .eq. 'P' ) then
            file_type = 'SPICE3'
            prev_binfiletype = binfiletype
            binfiletype = DT_SPICE3
         else if( file_type(7:7) .eq. '4') then
            file_type = 'SUPREM4'
            prev_binfiletype = binfiletype
            binfiletype = DT_SUPREM4
         else if( file_type(8:8) .eq. 'P') then
            file_type = 'SUPREM3S'
            prev_binfiletype = binfiletype
            binfiletype = DT_SUPREM3P
         else if( file_type(8:8) .eq. 'S') then
            file_type = 'SUPREM3P'
            prev_binfiletype = binfiletype
            binfiletype = DT_SUPREM3S
         endif
      else if( file_type(1:1) .eq. '2' ) then
         file_type = '2DOP'
         prev_binfiletype = binfiletype
         binfiletype = DT_2DOP
      else if( file_type(1:1) .eq. 'R' ) then
         file_type = 'RELOAD'
         prev_binfiletype = binfiletype
         binfiletype = DT_CALLABLE
      else if( file_type(1:2) .eq. 'AS' ) then
         file_type = 'ASCII'
         prev_binfiletype = binfiletype
         binfiletype = DT_ASCII
      else if( file_type(1:2) .eq. 'AN' ) then
         file_type = 'ANALYTICAL'
         prev_binfiletype = binfiletype
         binfiletype = DT_ANALYTICAL
      else if( file_type(1:1) .eq. 'G' ) then
         file_type = 'GRAPES'
         prev_binfiletype = binfiletype
         binfiletype = DT_GRAPES
      else if( file_type(1:1) .eq. 'T' ) then
         file_type = 'TIF'
         prev_binfiletype = binfiletype
         binfiletype = DT_TIF
      else if( file_type(1:1) .eq. 'F' ) then
         file_type = 'FAIM'
         prev_binfiletype = binfiletype
         binfiletype = DT_FAIM
      else if( file_type(1:2) .eq. 'MC' ) then
         file_type = 'MCP2D'
         prev_binfiletype = binfiletype
         binfiletype = DT_MCP2D
      else if( file_type(1:2) .eq. 'ME' ) then
         file_type = 'MEDICI'
         prev_binfiletype = binfiletype
         binfiletype = DT_MEDICI
      else
         write(luttyo,101) file_type
101      format(' ','Unknown data file type: ',a)
         file_type = 'MINIMOS'
         goto 100
      endif
c
      if( file_type .eq. 'PMI' ) then
         call get_str_with_default(savefilename,
     1     'Please enter PMI file name: ' , 1)
         savefile = savefilename
         call do_restore_plot(savefile,read_error)
         if( read_error ) then
            call get_char(answer,'Hit return to continue... ',0)
         endif
         return
      endif
c
      if( binfiletype .eq. DT_CALLABLE ) then
         call load_callable_data(.FALSE.,read_error)
         if( read_error ) then
            call get_char(answer,
     1      'Error loading data from memory, hit return to continue ',0)
         else
            lsavefile_read = .FALSE.
         endif
         return
      endif
c
      prev_binfilename = binfilename
      if( binfiletype .eq. DT_MEDICI ) then
        call get_str_with_default(binfilename,
     1        'Please enter the name of mesh file: ' , 1)
      else if( binfiletype .eq. DT_ANALYTICAL ) then
c       No filename needed
      else
        call get_str_with_default(binfilename,
     1        'Please enter the name of data file: ' , 1)
      endif
      write(luttyo,220)
220   format(' ')
c
      if( binfiletype .eq. DT_MINIMOS ) then
         call read_mini(.FALSE., read_error)
      else if( binfiletype .eq. DT_PROMIS ) then
         call read_promis(.FALSE., read_error)
      else if( binfiletype .eq. DT_VLSICAP ) then
         call read_vlsicap(.FALSE., read_error)
      else if( binfiletype .eq. DT_BAMBI ) then
         call read_bambi(read_error)
      else if( binfiletype .eq. DT_USEOUT ) then
         call read_useout(.FALSE., read_error)
      else if( binfiletype .eq. DT_SUPREM3S ) then
         type = 2
         call read_suprem3(type, read_error)
      else if( binfiletype .eq. DT_SUPREM3P ) then
         type = 1
         call read_suprem3(type, read_error)
      else if( binfiletype .eq. DT_SEDAN ) then
         call read_sedan(read_error)
      else if( binfiletype .eq. DT_2DOP ) then
         call read_2dop(.FALSE., read_error)
      else if( binfiletype .eq. DT_SUPREM4 ) then
         call read_sup4(.FALSE., read_error)
      else if( binfiletype .eq. DT_PISCES ) then
         call read_pisces(.FALSE., read_error)
      else if( binfiletype .eq. DT_TIF ) then
         call read_tif(.FALSE., read_error)
      else if( binfiletype .eq. DT_FAIM ) then
         call read_faim(.FALSE., read_error)
      else if( binfiletype .eq. DT_MCP2D ) then
         call read_mcp2d(.FALSE., read_error)
      else if( binfiletype .eq. DT_MEDICI ) then
         call read_medici(.FALSE., read_error)
      else if( binfiletype .eq. DT_ASCII ) then
         call read_ascii_file
         read_error = .NOT. have_user_1d_data
      else if( binfiletype .eq. DT_GRAPES ) then
         call read_grapesfile(.FALSE., read_error)
      else if( binfiletype .eq. DT_SPICE3 ) then
         call read_spice3file(read_error)
      else if( binfiletype .eq. DT_ANALYTICAL ) then
         call create_analytical
         read_error = .NOT. have_analytical_data
      endif
c
      if( read_error ) then
         call get_char(answer,
     1        'Do you wish to try again (y/n) [Y]: ',0)
         call upcase(answer)
         if( answer .ne. 'N' ) goto 100
      else
         lsavefile_read = .FALSE.
      endif
c
10000 continue
      return
      end
c
      subroutine read_sedan(read_error)
      implicit none
c
      include 'sysdep.cmn'
      include 'io.cmn'
c
      logical read_error
c
      write(luttyo,100)
100   format(' ','Sorry, not implemented yet...')
      read_error = .TRUE.
      return
      end
c
      subroutine zero_real_array(array,n,m)
      implicit none
c
c     clears real array dimensioned n by m.
      real array(1)
      integer n, m
c
      integer i
c
      do 10 i=1,n*m
10       array(i) = 0.0
c
      return
      end
c
      subroutine zero_dbl_array(array,n,m)
      implicit none
c
c     clears double precision array dimensioned n by m.
      double precision array(1)
      integer n, m
c
      integer i
c
      do 10 i=1,n*m
10       array(i) = 0.0
c
      return
      end
c
      subroutine copy_dbl_to_single(dbl_array,single_array,n,m)
      implicit none
c
c     copies double precision array dimensioned n by m into single
c     precision array.
      double precision dbl_array(1)
      real single_array(1)
      integer n, m
c
      integer i
c
      do 10 i=1,n*m
10       single_array(i) = SNGL( dbl_array(i) )
c
      return
      end
c
      subroutine select_line_from_2d(work, dest,
     1           x, z, nx, nz, plane_type, plane_index, isize)
      implicit none
c
      include 'const.h'
c
      integer nx, nz, plane_type, plane_index, isize
      double precision  work(isize*isize)
      real    dest(isize), x(isize), z(isize)
c
      integer i
c
      if( plane_type .eq. XY_PLANE ) then
         do i=1,nx
            dest(i) = SNGL(work(i + (plane_index-1)*nx))
         end do
      else if( plane_type .eq. YZ_PLANE ) then
         do i=1,nz
            dest(i) = SNGL(work(plane_index + (i-1)*nx))
         end do
      else
c        In the X-Z plane, the oxide or contact lines are contours
c        in the 2D plane. So do nothing here.
      endif
c
      return
      end
c
      subroutine read_3d_plane(binfil, work, dest,
     1                   nx, ny, nz, y_offset,
     2                   plane_index, plane_type, ni, nj,
     3                   isize, ierr, ioerr)
      implicit none
c
      include 'const.h'
c
      integer binfil, nx, ny, nz, y_offset,
     1        plane_type, plane_index,
     2        ni, nj, isize, ierr, ioerr
      double precision  work(isize*isize)
      real    dest(ni,nj)
c
      integer i, iz, k, zero_offset
c
      ierr = 0
      zero_offset = y_offset*nx
c
      if(plane_type .eq. XY_PLANE) then
         do iz=1,plane_index
            read(binfil,err=1000,end=2000,iostat=ioerr)
     1          (work(i),i=1,nx*ny)
         end do
         do k=1,ny
            do i=1,nx
               dest(i,k+y_offset) = SNGL(work(i + nx*(k-1)))
            end do
         end do
c        Skip over rest of planes
         do iz=plane_index+1,nz
            read(binfil,err=1000,end=2000,iostat=ioerr)
     1          (work(i),i=1,nx*ny)
         end do
      else if (plane_type .eq. YZ_PLANE) then
         do iz=1,nz
            read(binfil,err=1000,end=2000,iostat=ioerr)
     1          (work(i),i=1,nx*ny)
            do k=1,ny
               dest(iz,k+y_offset) = 
     1            SNGL(work(plane_index + nx*(k-1)))
            end do
         end do
      else if (plane_type .eq. XZ_PLANE) then
         do iz=1,nz
            read(binfil,err=1000,end=2000,iostat=ioerr)
     1          (work(i + zero_offset),i=1,nx*ny)
            do k=1,nx
               dest(k,iz) = 
     1            SNGL(work(nx*(plane_index-1) + k))
            end do
         end do
      else
         ierr = 1
      endif
c
      return
c
1000  ierr = 1
      return
c
2000  ierr = -1
      end
c
c
      subroutine read_3d_plane_sng(binfil, work, dest,
     1                   nx, ny, nz, y_offset,
     2                   plane_index, plane_type, isize, ierr, ioerr)
      implicit none
c
      include 'const.h'
c
      integer binfil, nx, ny, nz, y_offset,
     1        plane_type, plane_index,
     2        isize, ierr, ioerr
      real    work(isize*isize)
      real    dest(isize,isize)
c
      integer i, iz, k, zero_offset
c
      ierr = 0
      zero_offset = y_offset*nx
c
      if(plane_type .eq. XY_PLANE) then
         do iz=1,plane_index
            read(binfil,err=1000,end=2000,iostat=ioerr)
     1          (work(i),i=1,nx*ny)
         end do
         do k=1,ny
            do i=1,nx
               dest(i,k+y_offset) = work(i + nx*(k-1))
            end do
         end do
c        Skip over rest of planes
         do iz=plane_index+1,nz
            read(binfil,err=1000,end=2000,iostat=ioerr)
     1          (work(i),i=1,nx*ny)
         end do
      else if (plane_type .eq. YZ_PLANE) then
         do iz=1,nz
            read(binfil,err=1000,end=2000,iostat=ioerr)
     1          (work(i),i=1,nx*ny)
            do k=1,ny
               dest(iz,k+y_offset) = work(plane_index + nx*(k-1))
            end do
         end do
      else if (plane_type .eq. XZ_PLANE) then
         do iz=1,nz
            read(binfil,err=1000,end=2000,iostat=ioerr)
     1          (work(i + zero_offset),i=1,nx*ny)
            do k=1,nx
               dest(k,iz) = work(nx*(plane_index-1) + k)
            end do
         end do
      else
         ierr = 1
      endif
c
      return
c
1000  ierr = 1
      return
c
2000  ierr = -1
      end
c
      subroutine select_plane(plane_type, plane_index, plane_value,
     1                        plane_actual_value, xyz, nx, ny, nz)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'io.cmn'
c
      integer plane_type, plane_index, nx, ny, nz
      real    plane_value, plane_actual_value, xyz(COMP_ARRAYSIZE,3)
c
      integer bindex, nb, int_vec(2), getindex
      real real_vec(2), h
      character*1 ch, ans
c
      write(luttyo,100)
100   format(' ','Binary file contains 3D data'/
     1       ' ','You may view data along the X-Y (length) plane'/
     2       ' ','                            Y-Z (width)  plane'/
     3       ' ','                            X-Z (wafer)  plane')
200   call get_ints(int_vec,1,
     1   'Enter view direction (1=length; 2=width; 3=wafer) ',0)
      plane_type = int_vec(1)
      if( plane_type .lt. 1 .or. plane_type .gt. 3 ) then
         write(luttyo,210)
210      format(' ','Illegal plane type... re-enter')
         goto 200
      endif
c
      if( plane_type .eq. 1 ) then
         bindex = 3
         nb = nz
         ch = 'Z'
      else if( plane_type .eq. 2 ) then
         bindex = 1
         nb = nx
         ch = 'X'
      else
         bindex = 2
         nb = ny
         ch = 'Y'
      endif
c
      call get_char(ans,
     1   'Select plane by coordinate or plane index? (c/p) [C] ',0)
      call upcase(ans)
      if( ans .ne. 'P' ) then
         write(luttyo,300) xyz(1,bindex)*CM2MU, xyz(nb,bindex)*CM2MU, ch
300      format(' ','Device extends from ',f9.4,' to ', f9.4,
     1          ' microns in ',a1,' direction')
         call get_reals(real_vec,1,
     1      'Select cut plane coordinate to plot (in microns): ',0)
         h = real_vec(1)*MU2CM
         if( h .le. xyz(1,bindex) ) then
            plane_index = 1
         else if( h .ge. xyz(nb,bindex) ) then
            plane_index = nb
         else
            plane_index = getindex(h,xyz(1,bindex),nb)
            if( (xyz(plane_index+1,bindex) - h) .lt.
     1          (h - xyz(plane_index,bindex)))
     1         plane_index = plane_index + 1
         endif
         plane_value = real_vec(1)
         plane_actual_value = xyz(plane_index,bindex)
         write(luttyo,301) ch, plane_actual_value*CM2MU
301      format(' ','Closest plane was at ',a1,' = ',f9.4)
      else
         write(luttyo,400) nb
400      format(' ','Plane indices are 1 to ',i5)
         call get_reals(real_vec,1,
     1        'Enter plane index: ',0)
         plane_index = NINT(real_vec(1))
         if( plane_index .lt. 1 ) plane_index = 1
         if( plane_index .gt. nb ) plane_index = nb
         plane_actual_value = xyz(plane_index,bindex)
         plane_value = plane_actual_value
         write(luttyo,401) plane_index, ch, plane_actual_value*CM2MU
401      format(' ','Plane ',i5,' is located at ',a1,' = ',f9.4)
      endif
c
      return
      end
c
      subroutine batch_select_plane(cut_plane_axis,
     1           plane_index, plane_value,
     2           plane_type, plane_actual_value, xyz, nx, ny, nz,
     3           err)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'units.cmn'
      include 'io.cmn'
c
      integer plane_type, plane_index, nx, ny, nz
      real    plane_value, plane_actual_value, xyz(COMP_ARRAYSIZE,3)
      logical err
      character*1 cut_plane_axis
c
      integer bindex, nb, getindex
      real h
      character*1 ch
c
      err = .FALSE.
c
      if( cut_plane_axis .eq. 'Z' ) then
         plane_type = XY_PLANE
      elseif( cut_plane_axis .eq. 'X' ) then
         plane_type = YZ_PLANE
      elseif( cut_plane_axis .eq. 'Y' ) then
         plane_type = XZ_PLANE
      else
         write(luttyo,100) cut_plane_axis
100      format(' ','Error: Illegal cut plane axis: ',a)
         err = .TRUE.
         return
      endif
c
      if( plane_type .eq. XY_PLANE ) then
         bindex = 3
         nb = nz
         ch = 'Z'
      else if( plane_type .eq. YZ_PLANE ) then
         bindex = 1
         nb = nx
         ch = 'X'
      else
         bindex = 2
         nb = ny
         ch = 'Y'
      endif
c
      if( plane_index .eq. 0 ) then
         h = plane_value
         if( h .le. xyz(1,bindex) ) then
            plane_index = 1
         else if( h .ge. xyz(nb,bindex) ) then
            plane_index = nb
         else
            plane_index = getindex(h,xyz(1,bindex),nb)
            if( (xyz(plane_index+1,bindex) - h) .lt.
     1          (h - xyz(plane_index,bindex)))
     1         plane_index = plane_index + 1
         endif
         plane_actual_value = xyz(plane_index,bindex)
         write(luttyo,200) ch, plane_actual_value*pm_scale_cm2plot
200      format(' ','Closest plane was at ',a1,' = ',f9.4)
      else
         if( plane_index .lt. 1 ) plane_index = 1
         if( plane_index .gt. nb ) plane_index = nb
         plane_actual_value = xyz(plane_index,bindex)
         plane_value = plane_actual_value
         write(luttyo,401) plane_index, ch, plane_actual_value*CM2MU
401      format(' ','Plane ',i5,' is located at ',a1,' = ',f9.4)
      endif
c
      return
      end
c
      subroutine set_3d_index_limits
      implicit none
c
      include 'const.h'
      include 'mesh.cmn'
c
c     ni, nj are number of data points in the horizontal and vertical
c        cut plane
c     idisti,idistj are second indices of xyz (edxyz, dxyzm) arrays
c        which point to mesh spacings for horizontal and vertical
c        axes of the cut plane
c     charh, charv, charb are characters which identify absissa, ordinate
c        and "depth into the board" of the cut plane
c
      if( plane_type .eq. XY_PLANE ) then
         ni = nx
         nj = nytot
         idisti = 1
         idistj = 2
         charh = 'X'
         charv = 'Y'
         charb = 'Z'
      else if( plane_type .eq. YZ_PLANE ) then
         ni = nz
         nj = nytot
         idisti = 3
         idistj = 2
         charh = 'Z'
         charv = 'Y'
         charb = 'X'
      else if( plane_type .eq. XZ_PLANE ) then
         ni = nx
         nj = nz
         idisti = 1
         idistj = 3
         charh = 'X'
         charv = 'Z'
         charb = 'Y'
      else if( plane_type .eq. XT_PLANE ) then
         ni = nx
         nj = ny
         idisti = 1
         idistj = 2
         charh = 'T'
         charv = 'X'
c        Logic says charb should be "Y", but compare_plot needs this
c        charb to be 'T' to determine axis labels
         charb = 'T'
      else if( plane_type .eq. YT_PLANE ) then
         ni = nx
         nj = ny
         idisti = 1
         idistj = 2
         charh = 'T'
         charv = 'Y'
         charb = 'T'
      endif
c
      return
      end

c
      subroutine allocate_rect_arrays(nx,ny,ndata,ierr)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'data.cmn'
c
      integer nx, ny, ndata, ierr
c
c     locals
      integer i
c
c     Define the arrind() array entries so that it splits up the arrdat
c     array into appropriate chunks...
c
      if( nx .le. 0 .or. ny .le. 0) then
         ierr = 1
         return
      endif
c
      if( nx*ny*ndata .gt. ARR2DSIZE*NUM_DATA_ARRAYS ) then
c        Not enough static storage to satisfy request
         ierr = 2
         return
      endif
c
c     Compute offsets into arrdat which partition it into ndata pieces
c     of nx*ny each...
      do i=1,ndata
         arrind(i) = 1 + (i-1)*nx*ny
      end do
      ierr = 0
c
      return
      end
c
      subroutine allocate_tri_arrays(ntri,ndata,ierr)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'data.cmn'
c
      integer ntri, ndata, ierr
c
c     locals
      integer i
c
c     Define the arrind() array entries so that it splits up the arrdat
c     array into appropriate chunks...
c
      if( ntri .le. 0 ) then
         ierr = 1
         return
      endif
c
      if( ntri*ndata .gt. ARR2DSIZE*NUM_DATA_ARRAYS ) then
c        Not enough static storage to satisfy request
         ierr = 2
         return
      endif
c
c     Compute offsets into arrdat which partition it into ndata pieces
c     of ntri each...
      do i=1,ndata
         arrind(i) = 1 + (i-1)*ntri
      end do
      ierr = 0
c
      return
      end
