      SUBROUTINE contrd(iunit, chimin, stepsubx, stepsuby, dset,
     &                  dset_new)
 
c	fwj haberl   15-MAR-1989 10:31:44
c  Read STEPPAR package from XSPEC history file for plotting
c  of chi-square contours
 
c	implicit none
      REAL xsq(256, 256)
      LOGICAL*4 qndone
      INTEGER*4 iunit, len, ind, nx, ny, i, j, ij, isub, nsub, ierrsf,
     &          infoar(4)
      INTEGER*4 dset, dset_new
      CHARACTER*12 pkgtyp
      REAL*4 chimin
      REAL*4 delxin, delyin, xsqmax, xsqmin, xmin, xmax, ymin, ymax
      INCLUDE '../inc/historysf.inc'
      INCLUDE 'xplot.inc'
      RECORD /his_steppar/step
      RECORD /his_steppar_sub1/stepsubx, stepsuby
      RECORD /his_steppar_sub2/stepsub
 
      ind = 0
      pkgtyp = 'steppar'
      len = len_his_steppar
      CALL nxpksf(iunit, pkgtyp, ind, infoar, nsub, step, len, .TRUE.,
     &            ierrsf)
      IF (ierrsf.NE.0) THEN
         IF (ierrsf.EQ.7) THEN
            WRITE (*, *)
     &                  ' Unable to find a ''steppar'' package in file.'
         ELSE
            WRITE (*, *) ' SF read error in looking for package.'
         ENDIF
         GOTO 999
      ENDIF
      IF (step.nspar.le.1) THEN
         WRITE (*, *) ' Need two dimensions to be stepped'
         GOTO 999
      ELSEIF (step.nspar.gt.2) THEN
         WRITE (*, *)
     &               ' There are more than two dimensions being stepped'
     &               , ' Using the first two.'
      ENDIF
      chimin = step.chimin
      WRITE (*, '(2a,g15.5)') indent, 'Minimum chi square from fit: ',
     &                        chimin
      len = len_his_steppar_sub1
      CALL rsubsf(iunit, stepsubx, len, ierrsf)
      len = len_his_steppar_sub1
      CALL rsubsf(iunit, stepsuby, len, ierrsf)
      DO i = 3, step.nspar
         CALL sksbsf(iunit, ierrsf)
      ENDDO
      CALL stepdl(stepsubx, nx, delxin)
      CALL stepdl(stepsuby, ny, delyin)
c
c  Set bad-data flag in array.
c
      DO j = 1, ny
         DO i = 1, nx
            xsq(i, j) = -1.
         ENDDO
      ENDDO
      nsub = nsub - step.nspar
      IF (nsub.LE.0) nsub = 1000000
c 	** fill in the grid
      qndone = .TRUE.
      isub = 0
      DO WHILE (qndone)
         isub = isub + 1
         len = len_his_steppar_sub2_ini + 4*step.nspar
         CALL rsubsf(iunit, stepsub, len, ierrsf)
cd
c	    write (*,*) stepsub.chisqr, stepsub.par(1), stepsub.par(2), ierrsf
cd
         IF (ierrsf.EQ.0) THEN
            CALL getpos(stepsub.par(1), stepsubx, nx, delxin, i)
            CALL getpos(stepsub.par(2), stepsuby, ny, delyin, j)
            IF (xsq(i,j).EQ.-1.) THEN
               xsq(i, j) = stepsub.chisqr
            ELSE
               xsq(i, j) = min(xsq(i,j), stepsub.chisqr)
            ENDIF
            qndone = isub.LT.nsub
         ELSE
            qndone = .FALSE.
         ENDIF
      ENDDO
c
c  Find minimum and maximum.
c
      xsqmin = 1.e30
      xsqmax = -1.e30
      DO i = 1, nx
         DO j = 1, ny
            IF (xsq(i,j).GT.xsqmax) xsqmax = xsq(i, j)
            IF (xsq(i,j).LT.xsqmin) xsqmin = xsq(i, j)
         ENDDO
      ENDDO
      IF (xsqmin.LT.chimin) THEN
         WRITE (*, '(2a,g15.5)') indent,
     &                     'Lower chi square in steppar package found: '
     &                     , xsqmin
      ENDIF
c
c  Replace any points which still have -1 with maximum value.
c
      ij = 0
      DO i = 1, dset
         ij = ij + nbinz(i, 1)*nbinz(i, 2)
      ENDDO
      DO j = 1, ny
         DO i = 1, nx
            ij = ij + 1
            IF (xsq(i,j).EQ.-1) xsq(i, j) = xsqmax
            datax(ij) = xsq(i, j)
         ENDDO
      ENDDO
cd
c	do i=1, ny
c	    write (*,*) (xsq(j,i),j=1,nx)
c	end do
cd
      dset_new = 1
      nbinz(dset+1, 1) = nx
      nbinz(dset+1, 2) = ny
      RETURN
c	**destination for errors
 999  CONTINUE
      dset_new = 0
      RETURN
      END
