******************************************
	subroutine avgrms( n, x, bar, rms )
* unweighted mean and rms of n data points x(i)
* Input:
*	n	i4 number of data
*	x(n)	r4 data
* output:
*	bar	r4 mean value
*	rms	r4 root-mean-square
* Running update method (B.P.Welford 1962)
* Seems to be more robust to round-off errors.
* 2011 Feb Keith Horne @ St.Andrews
	real*4 x(*)
	real*8 sum, sum2, add
      if( n .le. 1 ) then
	bar = x(1)
	rms = 0.
	return
      end if
	sum1 = x(1)
	sum2 = 0.d0
      do i=2,n
	add = dble( x(i) ) - sum1
	sum1 = sum1 + add / i
	sum2 = sum2 + add * ( dble( x(i) ) - sum1 )
      end do
	bar = sum1
	var = sum2 / max( 1, n - 1 )	
	rms = sqrt( var )
	return
	end

******************************************
	subroutine avgrms1( n, x, bar, rms )
* unweighted mean and rms of n data points x(i)
* Input:
*	n	i4 number of data
*	x(n)	r4 data
* output:
*	bar	r4 mean value
*	rms	r4 root-mean-square
* 1999 Sep Keith Horne @ St.Andrews
* 2011 Feb KDH @ StA - x(i)**2 -> x(i)*x(i)
	real*4 x(*)
	real*8 sum, sum2
	bar = x(1)
	rms = 0.
	if( n.le.1 ) return
	sum = 0.d0
	sum2 = 0.d0
      do i=1,n
	sum = sum + x(i)
	sum2 = sum2 + x(i) * x(i)
      end do
	bar = sum / n
	var = ( sum2 - sum * sum / n ) / max( 1, n-1 )
	rms = sqrt( max( 0., var ) )
	return
	end
********************************************
	subroutine optavg( n, dat, sig, avg, sigavg )
* optimal average
* Input:
*	n	i4 number of data values
*	dat(n)	r4 data values
*	sig(n)	r4 error bars (negative to ignore point)
* Output:
*	avg	r4 optimal average
*	sigavg	r4 error bar of optimal average
* 2000 Sep Keith Horne @ St.Andrews
* 2016 Aug KDH @ St.A - prescaling
	real*4 dat(*), sig(*)
	real*8 sum, sumw, div
	avg = 0.
	sigavg = -1.
	if( n .le. 0 ) return
* prescale divisor
	div = safedivisor( n, dat, sig )
	sum = 0.d0
	sumw = 0.d0
      do i = 1, n
	s = sig( i ) / div
      if( s .gt. 0. ) then
	w = 1. / ( s * s )
	sumw = sumw + w
	sum  = sum  + w * dat( i )
      end if
      end do
	if( sumw .le. 0.d0 ) return
	avg = sum / sumw
	sigavg = div / dsqrt( sumw )
	return
	end
********************************************
	subroutine optavg1( n, dat, sig, avg, sigavg )
* optimal average (old version with no prescaling)
* Input:
*	n	i4 number of data values
*	dat(n)	r4 data values
*	sig(n)	r4 error bars (negative to ignore point)
* Output:
*	avg	r4 optimal average
*	sigavg	r4 error bar of optimal average
* 2000 Sep Keith Horne @ St.Andrews
	real*4 dat(*), sig(*)
	real*8 sum, sumw
	avg = 0.
	sigavg = -1.
	if( n .le. 0 ) return
	sum = 0.d0
	sumw = 0.d0
      do i = 1, n
	s = sig( i )
      if( s .gt. 0. ) then
	w = 1. / ( s * s )
	sumw = sumw + w
	sum = sum + w * dat(i)
      end if
      end do
	if( sumw .le. 0.d0 ) return
	avg = sum / sumw
	sigavg = 1.d0 / dsqrt( sumw )
	return
	end

********************************************
	subroutine optscl( n, dat, sig, pfl, scl, sigscl )
* scale profile to fit data
* Input:
*	n	i4 number of data points
*	dat(n)	r4 data values
*	sig(n)	r4 1-sigma error bars (<0 to omit data point)
*	pfl(n)	r4 profile to be scaled to fit
* Output:
*	scl	r4 scale factor
*	sigscl	r4 1-sigma error bar on scale factor
* 2000 Oct Keith Horne @ St.Andrews
	real*4 dat(*), sig(*), pfl(*)
	real*8 sum, sum1
	scl = 0.
	sigscl = -1.
	if( n.le.0 ) return
	sum = 0.d0
	sum1 = 0.d0
      do i=1,n
      if( sig(i) .gt. 0. ) then
	wgt = pfl(i) / sig(i)**2
	sum = sum + dat(i) * wgt
	sum1 = sum1 + pfl(i) * wgt
      end if
      end do
	if( sum1 .le. 0.d0 ) return
	scl = sum / sum1
	var = 1. / sum1
	sigscl = sqrt( var )
	return
	end
* =================================================
	subroutine avgrmsx( n, dat, sig, good, avg, rms, sigavg, sigrms )
* Maximum likelihood estimator of the mean and intrinsic rms
* given independent samples with error bars
* In:	n	i4 number of data
*	dat(n)	r4 data
*	sig(n)	r4 1-sigma error bars (<0 to skip)
* Out:	good	r4 number of "good" data
*	avg	r4 mean value
*	sigavg	r4 1-sigma error bar
*	rms	r4 excess rms
*	sigrms	r4 1-sigma of excess rms
* 2014 Aug Keith Horne @ St Andrews
* 2014 Nov KDH @ StA - improve estimate of sigrms
* 2016 Mar KDH @ StA - catch crash if no variance
* 2016 Mar KDH @ StA - scale to <1/sig^2>=1, to avoid vmin=0.
* 2018 Mar KDH @ StA - verbose

	real*4 dat( n ), sig( n )
	real*8 sum, sumg, sum1, sum2, sum3
	logical verbose/.false./
* trap no data
	good = 0.
	avg = 0.
	rms = 0.
	sigavg = -1.
	sigrms = -1.
	if( n .le. 0 ) return

* mean of positive error bars
	ng = 0
	sum = 0.d0
      do k = 1, n
	e = sig( k )
      if( e .gt. 0. ) then
	sum = sum + e
	ng = ng + 1
      end if
      end do
	e1 = sum / max( 1, ng ) 

* abort no valid data
      if( ng .le. 0 ) then
	if( verbose ) write(*,*) '** ERROR in AVGRMSX: NO >0 ERROR BARS n', n, ' ng', ng
	return	
      end if

* abort e1 = 0
      if( e1 .le. 0. ) then
	write(*,*) '** ERROR in AVGRMSX. MEAN ERROR BAR e1', e1
	write(*,*) '** sum', sum, ' n', n, ' ng', ng
	return
      end if

* initial estimate
	sum1 = 0.d0
	sum2 = 0.d0
      do k = 1, n
	e = sig( k )
      if( e .gt. 0. ) then
	d = dat( k )
	w = ( e1 / e ) ** 2
	sum1 = sum1 + w * d
	sum2 = sum2 + w	
      end if
      end do

* abort no valid data
      if( sum2 .le. 0. .or. ng .le. 0 ) then
	write(*,*) '** ERROR in AVGRMSX. n', n, ' ng', ng, ' sum2=', sum2
	return	
      end if

* optimal average and its variance
	avg = sum1 / sum2
	sigavg = e1 / dsqrt( sum2 )
* scale factor ( makes <1/sig^2>=1 )
	w = ng / sum2
	e1 = e1 * sqrt( w )
* estimate for rms
	rms = e1

* convergence threshold
	tiny = 2.e-5
* lower limit for extra rms
	rmin = tiny * e1

* KDH : NO LONGER NEEDED ?
      if( rmin .le. 0. ) then
	write(*,*) '** WARNING in AVGRMSX : rmin not positive :('
	write(*,*) ' sum1', sum1, ' sum2', sum2, ' ng', ng
	write(*,*) ' avg', avg, ' sigavg', sigavg
	write(*,*) '  e1', e1, ' rmin', rmin
      end if

* iterate to converge avg and rms
	nloop = 1000
      do loop = 1, nloop
* stow for convergence test
	oldavg = avg
	oldrms = rms
* zero sums
	sum1 = 0.d0
	sum2 = 0.d0
	sumg = 0.d0
	sum  = 0.d0
	sum3 = 0.d0
* scaled avg and rms
	a = avg / e1
	r = rms / e1
	rr = r * r
* data loop
      do k = 1, n
* reject data with negative error bar
	e = sig( k )
      if( e .gt. 0. ) then
* scale data and error bar
	d = dat( k ) / e1
	e = e / e1
* "goodness" of this data point
	g = rr / ( rr + e * e )
	x = g * ( d - a )
	xx = x * x
* increment sums
	sum1 = sum1 + g * d
	sumg = sumg + g
	sum  = sum  + g * g
	sum2 = sum2 + xx
	sum3 = sum3 + g * xx
      end if
* next datum
      end do

* "good" data points ( e.g. with sig < rms )
	good = sumg

* KDH: CATCH SUMG=0
      if( sumg .le. 0.d0 ) then
	write(*,*) '** ERROR in AVGRMSX. NON-POSITIVE SUM(G)=', sumg
	write(*,*) '** loop', loop, ' ndat', n, ' ngood', ng
	write(*,*) '** sumg', sumg, ' sum1', sum1, ' sum', sum
	write(*,*) '**  r', r, ' rr', rr
	write(*,*) '** dat(:)', (dat(i),i=1,n)
	write(*,*) '** sig(:)', (sig(i),i=1,n)
	write(*,*) 'Should we abort here ?'
	stop
      end if

* revise avg and rms
	a = sum1 / sumg
	r = dsqrt( sum2 ) / dsqrt( sumg )
	r = max( rmin, r )
* error bars on scaled avg and rms
	siga = r / sqrt( sumg )
	g = 2.d0 * ( 2.d0 * sum3 / sum2 * sumg - sum )
	sigr = r
	if( g .gt. 0. ) sigr = r / sqrt( g )

* restore scaling
	avg = a * e1
	rms = r * e1
	sigavg = siga * e1
	sigrms = sigr * e1


* KDH: CATCH SUM2=0 (NO VARIANCE)
      if( sum2 .le. 0.d0 ) then
	rms = 0.
	sigrms = -1.
	return
      end if

* converge when test < 1
      if( loop .gt. 1 ) then
* step
	safe = 0.9
	avg = oldavg * ( 1. - safe ) + safe * avg
	rms = oldrms * ( 1. - safe ) + safe * rms
* chi
	chiavg = ( avg - oldavg ) / sigavg
	chirms = ( rms - oldrms ) / sigrms
	test = max( abs( chiavg ), abs( chirms ) ) / tiny

* report on last 5 iterations
      if( loop .gt. nloop - 3 ) then
	write(*,*) 'Loop', loop, ' of', nloop, ' in AVGRMSX'
	write(*,*) 'Ndat', n, ' Ngood', good, ' Neff', g
	write(*,*) ' avg', avg, ' rms', rms
	write(*,*) ' +/-', sigavg, ' +/-', sigrms
	write(*,*) ' chiavg', chiavg, ' chirms', chirms, ' test', test
	if( test .lt. 1. ) write(*,*) '** CONVERGED ** :))'
      end if
* converged
	if( test .lt. 1. ) return
      end if

* quit if rms vanishes
	if( rms .le. 0. ) return

* next loop
      end do

* failed to converge
	write(*,*) '** AVGRMSX FAILED TO CONVERGE :((((('
	return
	end
* =======================================================
	subroutine runoptavg( n, t, d, s, hw, nf, tf, f, sf )
* running optimal average of time series data
* In:	n	i4 number of epochs
*	t(n)	r4 epochs
*	d(n)	r4 data
*	s(n)	r4 1-sigma error bars (<0 to omit)
*	hw	r4 window half-width
*	nf	i4 number of fit epochs
*	tf(nf)	r4 epochs at which to evaluate the fit
* Out:	f(nf)	r4 running optimal average (roa)
*	sf(nf)	r4 1-sigma uncertainty of the roa
* 2003 Jun Keith Horne @ St.Andrews
* 2020 Jun KDH @ StA - call ROAOW - for backwards compatibility
	real*4  t(*), d(*), s(*), tf(*), f(*), sf(*)
	external winb, winc, wing, winl
	call ROAOW( winl, hw, n, t, d, s, nf, tf, f, sf )
	return
	end
* =======================================================
	subroutine ROA( shape, hw, n, t, d, s, nf, tf, f, sf )
* ROA = Running Optimal Average of time series data
* In:	shape	c* window shape (one of B,C,G,L)
*	hw	r4 window half-width
C	n	i4 number of epochs
*	t(n)	r4 epochs
*	d(n)	r4 data
*	s(n)	r4 1-sigma error bars (<0 to omit)
*	nf	i4 number of fit epochs
*	tf(nf)	r4 times at which to evaluate the fit
* Out:	f(nf)	r4 running optimal average (roa)
*	sf(nf)	r4 1-sigma uncertainty of the roa
* 2020 Jun KDH @ St.Andrews - call ROAOW
	real*4  t(*), d(*), s(*), tf(*), f(*), sf(*)
	character*(*) shape
	character*1 w
	external winb, winc, wing, winl
	w = shape
	call upcase( w, w )
* Gaussian
      if(      w .eq. 'G' ) then
	call ROAOW( wing, hw, n, t, d, s, nf, tf, f, sf )
* Boxcar
      else if( w .eq. 'B' ) then
	call ROAOW( winb, hw, n, t, d, s, nf, tf, f, sf )
* Lorentzian
      else if( w .eq. 'L' ) then
	call ROAOW( winl, hw, n, t, d, s, nf, tf, f, sf )
* 1/Cosh
      else if( w .eq. 'C' ) then
	call ROAOW( winc, hw, n, t, d, s, nf, tf, f, sf )
* revert to optimal average
      else
	write(*,*) '** RAO window shape not recognised. ', shape(:len1(shape))
	call optavg( n, d, s, avg, sig )
      do i = 1, nf
	f(i) = avg
	sf(i) = sig
      end do
      end if
	return
	end
* =======================================================
	function wing( x )
* Gaussian window function
* 2020 Jun KDH @ St.Andrews
	wing = exp( - 0.5 * x * x )
	return
	end
* =======================================================
	function winb( x )
* Boxcar window function
* 2020 Jun KDH @ St.Andrews
	winb = 1.
	if( abs(x) .gt. 1. ) winb = 0.
	return
	end
* =======================================================
	function winc( x )
* 1/Cosh window function
* 2020 Jun KDH @ St.Andrews
	winc = 1. / cosh( x )
	return
	end
* =======================================================
	function winl( x )
* Lorentzian window function
* 2020 Jun KDH @ St.Andrews
	winl = 1. / ( 1. + x * x )
	return
	end

* =======================================================
	subroutine ROAP( shape, hw, n, t, s, p )
* ROAP = Running Optimal Average Parameters
* Effective number of parameters fitted to each datum
* for a running optimal average of time series data
* In:	shape	c* window shape (one of B C G L)
*	hw	r4 window half-width ( x = t / hw )
*	n	i4 number of epochs
*	t(n)	r4 epochs
*	s(n)	r4 1-sigma error bars (<0 to omit)
* Out:	p(n)	r4 effective number of parameters
* 2020 Jun Keith Horne @ St Andrews - call ROAPOW
	real*4  t(*), s(*), p(*)
	character*(*) shape
	character*1 w
	external winb, winc, wing, winl
	w = shape
	call upcase( w, w )
* Gaussian
      if(      w .eq. 'G' ) then
	call ROAPOW( wing, hw, n, t, s, p )
* Boxcar
      else if( w .eq. 'B' ) then
	call ROAPOW( winb, hw, n, t, s, p )
* Lorentzian
      else if( w .eq. 'L' ) then
	call ROAPOW( winl, hw, n, t, s, p )
* 1/Cosh
      else if( w .eq. 'C' ) then
	call ROAPOW( winc, hw, n, t, s, p )
* invalid
      else
	write(*,*) '** RAO window shape not recognised. ', shape(:len1(shape))
      do i = 1, n
	p(i) = 0.
      end do
      end if
	return
	end
* =======================================================
	function ROAP1( shape, hw, n, t, s, j )
* ROAP1 = Running Optimal Average Parameters 1 datum
* Effective number of parameters fitted to datum i
* for a running optimal average of time series data
* In:	shape	c* window shape (one of B C G L)
*	hw	r4 window half-width ( x = t / hw )
*	n	i4 number of epochs
*	t(n)	r4 epochs
*	s(n)	r4 1-sigma error bars (<0 to omit)
*	j	i4 index of datum j
* Out:	roap1	r4 effective number of parameters
* 2020 Jun Keith Horne @ St Andrews - call ROAPOW1
	real*4  t(*), s(*)
	character*(*) shape
	character*1 w
	external winb, winc, wing, winl
	w = shape
	call upcase( w, w )
* Gaussian
      if(      w .eq. 'G' ) then
	roap1 = ROAPOW1( wing, hw, n, t, s, j )
* Boxcar
      else if( w .eq. 'B' ) then
	roap1 = ROAPOW1( winb, hw, n, t, s, j )
* Lorentzian
      else if( w .eq. 'L' ) then
	roap1 = ROAPOW1( winl, hw, n, t, s, j )
* 1/Cosh
      else if( w .eq. 'C' ) then
	roap1 =  ROAPOW1( winc, hw, n, t, s, j )
* invalid
      else
	write(*,*) '** RAO window shape not recognised. ', shape(:len1(shape))
	roap1 = 0.
      end if
	return
	end
* =======================================================
	subroutine ROAPOW( win, hw, n, t, s, p )
* ROAPOW = Running Optimal Average Parameters Own Window
* Effective number of parameters fitted to each datum
* for a running optimal average of time series data
* In:	win(x)	r4 window function ( 0 < win(x) < 1 )
*	hw	r4 window half-width ( x = t / hw )
*	n	i4 number of epochs
*	t(n)	r4 epochs
*	s(n)	r4 1-sigma error bars (<0 to omit)
* Out:	p(n)	r4 effective number of parameters
* 2020 Jun Keith Horne @ St Andrews - call ROAPOW1 n times
	real*4  t(*), s(*), p(*)
	external win
      do i = 1, n
	p( i ) = ROAPOW1( win, hw, n, t, s, i )
      end do
	return
	end
* =======================================================
	function ROAPOW1( win, hw, n, t, s, j )
* ROAPOW1 = Running Optimal Average Parameters Own Window 1 Datum
* Effective number of parameters fitted to datum j
* for a running optimal average of time series data
* In:	win(x)	r4 window function ( 0 < win(x) < 1 )
*	hw	r4 window half-width ( x = t / hw )
*	n	i4 number of epochs
*	t(n)	r4 epochs
*	s(n)	r4 1-sigma error bars (<0 to omit)
*	j	i4 index of datum j 
* Out:	roapow1	r4 roa parameters fitted to datum j
* 2020 Jun Keith Horne @ St Andrews 
	real*4  t(*), s(*)
	external win
	real*8 sum
	tj = t( j )
	sum = 0.d0
	wj = 0.
      do i = 1, n
	si = s( i )
      if( si .gt. 0. ) then
	ti = t( i )
	x = ( tj - ti ) / hw
	w = win( x ) / si / si
	sum = sum + w
	if( i .eq. j ) wj = w
      end if
      end do
	roapow1 = wj / sum
	return
	end

* =======================================================
	subroutine ROAOW( win, hw, n, t, d, s, nf, tf, f, sf )
* ROAOW = Running Optimal Average Own Window
* In:	win(x)	r4 window function ( 0 < win(x) < 1 )
*	hw	r4 window half-width ( x = t / hw )
*	n	i4 number of epochs
*	t(n)	r4 epochs
*	d(n)	r4 data
*	s(n)	r4 1-sigma error bars (<0 to omit)
*	nf	i4 number of fit epochs
*	tf(nf)	r4 times at which to evaluate the fit
* Out:	f(nf)	r4 running optimal average
*	sf(nf)	r4 1-sigma uncertainty in running optimal average
* 2003 Jun Keith Horne @ St.Andrews
* 2020 Jun KDH @ StA - adapt from runoptavg
* 2020 Jun KDH @ StA - external function win
* 2020 Jun KDH @ StA - tiny
	real*4  t(*), d(*), s(*), tf(*), f(*), sf(*)
	external win
	real*8 sum, sumw
	tiny = 1.e-12
	tiny = 0.
	tiny = 1.e-20
	tiny = 1.e-15
* nothing to do
	if( n .le. 0 .or. nf .le. 0 ) return
* optimal average
	call optavg( n, d, s, avg, sig )
* validate half-width
      if( hw .le. 0. ) then
	write(*,*) '** ROA INVALID WINDOW HALF-WIDTH', hw
      do i = 1, nf
	 f( i ) = avg
	sf( i ) = sig
      end do
	return
      end if

* fit epoch j
      do j = 1, nf
	tj = tf( j )
	sum = 0.d0
	sumw = 0.d0
* data epoch i
      do i = 1, n
* omit bad data
	si = s( i )
      if( si .gt. 0. ) then
	ti = t( i )
	x = ( ti - tj ) / hw
	w = win( x )
	w = max( w, tiny )
	w = w / si / si
	sumw = sumw + w
	sum  = sum  + w * d( i )
      end if
* next data epoch i
      end do
* stow
      if( sumw .gt. 0.d0 ) then
	 f( j ) = sum / sumw
	sf( j ) = dsqrt( 1.d0 / sumw )
      else
	 f( j ) = avg
	sf( j ) = sig
      end if
* next fit epoch j
      end do
	return
	end
* ==================================
	function safedivisor( n, dat, sig )
* find a "safe" prescale factor
* Input:
*	n	i4 number of data values
*	dat(n)	r4 data values
*	sig(n)	r4 error bars (negative to ignore point)
* Output:
*	safedivisor	r4 "safe" divisor for prescaling
* 2016 Aug Keith Horne @ St.Andrews
	real*4 dat( * )
	real*4 sig( * )
* range of dat and sig where sig > 0
	np = 0
      do i = 1, n
	s = sig( i )
      if( s .gt. 0. ) then
	d = dat( i )
      if( np .eq. 0 ) then
	dmn = d
	dmx = d
	smn = s
	smx = s
      else
	dmn = min( dmn, d )
	dmx = max( dmx, d )
	smn = min( smn, s )
	smx = max( smx, s )
      end if
	np = np + 1	
      end if
      end do
* half-range of good data
	s = ( dmx - dmn ) / 2.
* smallest error bar
	if( np .gt. 0 ) s = smn
* or no safe divisor
	if( s .le. 0. ) s = 1.
	safedivisor = s
	return
	end
*******************************************
	subroutine mnmx( n, x, xmn, xmx )
* find min and max values
* Input:
*	n	i4 number of data values
*	x(n)	r4 data values
* Output:
*	xmn	r4 minimum value
*	xmx	r4 maximum value
	real*4 x(*)
	xmn = x(1)
	xmx = xmn
      do i=1,n
	xmn = min( xmn, x(i) )
	xmx = max( xmx, x(i) )
      end do
	return
	end
*******************************************
	subroutine mnmx_index( n, x, imn, imx )
* find index of min and max values
* Input:
*	n	i4 number of data values
*	x(n)	r4 data values
* Output:
*	imn	r4 index of first min value
*	imx	r4 index of last max value
* 2018 Oct Keith Horne @St.Andrews
	real*4 x(*)
	imn = 1
	imx = 1
      do i=1,n
	if( x(i) .lt. x(imn) ) imn = i
	if( x(i) .ge. x(imx) ) imx = i
      end do
	return
	end
*******************************************
	function minpos( n, x )
* find index of min value
* Input:
*	n	i4 number of data values
*	x(n)	r4 data values
* Output:
*	minpos	r4 index of first min value
* 2018 Oct Keith Horne @St.Andrews
	real*4 x(*)
	m = 1
      do i=1,n
	if( x(i) .lt. x(m) ) m = i
      end do
	minpos = m
	return
	end
*******************************************
	function maxpos( n, x )
* find index of max value
* Input:
*	n	i4 number of data values
*	x(n)	r4 data values
* Output:
*	maxpos	r4 index of last max value
* 2018 Oct Keith Horne @St.Andrews
	real*4 x(*)
	m = 1
      do i=1,n
	if( x(i) .gt. x(m) ) m = i
      end do
	maxpos = m
	return
	end
*******************************************
	subroutine mnmxpos( n, x, xmn, xmx )
* find min and max of positive values
* Input:
*	n	i4 number of data values
*	x(n)	r4 data values
* Output:
*	xmn	r4 minimum value
*	xmx	r4 maximum value
* 2004 Apr Keith Horne @ St.And
* 2012 Sep KDH @ StA fix bug when x(1) < 0
	real*4 x(*), xmn, xmx
	npos = 0
      do i=1,n
      if( x(i) .gt. 0. ) then
	npos = npos + 1
      if( npos .eq. 1 ) then
	xmn = x(i)
	xmx = xmn
      else
	xmn = min( xmn, x(i) )
	xmx = max( xmx, x(i) )
      end if
      end if
      end do
      if( npos .lt. 1 ) then
	xmn = 0.
	xmx = 0.
      end if
	return
	end
*******************************************
	subroutine mnmxi( n, i, imn, imx )
* find min and max values for integer array 
* Input:
*	n	i4 number of data values
*	i(n)	i4 data values
* Output:
*	imn	r4 minimum value
*	imx	r4 maximum value
* 2004 Feb Keith Horne @ St-Andrews
	integer*4 i( * )
	imn = i( 1 )
	imx = imn
      do k=1,n
	imn = min( imn, i( k ) )
	imx = max( imx, i( k ) )
      end do
	return
	end

*******************************************
	subroutine dmnmx( n, x, xmn, xmx )
* find min and max of real*8 values
* Input:
*	n	i4 number of data values
*	x(n)	r8 data values
* Output:
*	xmn	r8 minimum value
*	xmx	r8 maximum value
	real*8 x(*), xmn, xmx
	xmn = x(1)
	xmx = xmn
      do i=1,n
	xmn = min( xmn, x(i) )
	xmx = max( xmx, x(i) )
      end do
	return
	end

************************************************
	function rang( avg, rms, iseed )
* Gaussian random number generator.
* Input:
*	avg	r4	mean value
*	rms	r4	standard deviation
*	iseed	i4	seed integer for ran( iseed ) 
* Output:
*	rang	r4	gaussian random number
*	iseed	i4	seed integer from ran( iseed )
*
* Box-muller transform of 2 independent random variables
* 1978 Peter Young @ Caltech -- original version GAUSS
* 1995 Keith Horne @ St-Andrews -- adapted from GAUSS.
* 2001 Sep KDH @ St-And - tidy up
* 2016 Jul KDH @ Iceland - ran -> rankdh (ran fails w gfortran)
	logical	newpair /.true./
	save y, u, newpair
	rang = avg
	if( rms .le. 0. ) return
      if ( newpair ) then
	r2 = 10.
      do while ( r2 .ge. 1. )
	x = 2. * rankdh( iseed ) - 1.
	y = 2. * rankdh( iseed ) - 1.
	r2 = x * x + y * y
      end do
	u = sqrt( -2. * alog( r2 ) / r2 )
	rang = x * u
      else
	rang = y * u
      end if
	newpair = .not. newpair
	rang = avg + rms * rang
	return
	end
*========================================
        subroutine rangtest( iseed )
* check avg and rms of gaussian random numbers
* 2016 Jul Keith Horne @ Iceland
        real*8 sum, sum2
        logical fail
* sample
        n = 10000
	sum = 0.d0
	sum2 = 0.d0
      do i = 1, n
        add = rang( 0., 1., iseed )
        sum = sum + add
        sum2 = sum2 + add * add
      end do
* mean
        avg = sum / n
	siga = 1. / sqrt( float( n ) )
	chia = avg / siga
* rms
	sum2 = ( sum2 - ( sum / n ) * sum ) / ( n - 1 )
        rms = dsqrt( sum2 )
	sigr = sqrt( 2. / ( n - 1 ) )
	chir = ( rms - 1. ) / sigr
* report
	write(*,*) 'avg', avg, ' =', chia, ' sigma'
	write(*,*) 'rms', rms, ' =', chir, ' sigma'
* test
        fail = max( abs( chia ), abs( chir ) ) .gt. 10.
      if( fail ) then
        write(*,*) '** RANG TEST FAILED  :('
	stop
      end if
        return
        end
*==================================
	real*4 function ranu( a, b, iseed )
* random number uniform on (a,b)
* In:	a	r4 lower limit
* 	b	r4 upper limit
*	iseed	i4 seed integer
* Out:  ranu	r4 random number uniform on (a,b)
*	iseed	i4 seed integer for use on next call
*
* 2020 May Keith Horne @ St.Andrews
	p = rankdh( iseed )
	ranu = a * ( 1. - p ) + p * b
	return
	end
*==================================
	real*4 function rankdh( iseed )
* random number uniform on (0-1)
* hi-level routines call this one, which can
* be used to select different lo-level generators.
* Input:
*	iseed	i4  repeat sequence if iseed < 0, otherwise ignore
*
* 2001 Jul Keith Horne @ St.Andrews
* 2006 Aug KDH @ St.A - use numerical recipes routine ran1
* 2008 Oct KDH @ St.A - try ran3
* 2016 Jan KDH @ Iceland - irand
	logical firstcall/.true./
	save firstcall
	data irand/0/
	is = iseed
      if( irand .ne. 0 ) then
* system-supplied
	ran = rand( is )
      else
* numerical recipes routine 
      if( firstcall ) then
	is = - iabs( is )
	firstcall = .false.
      end if
	ran = ran3( is )
      end if
	iseed = is
	rankdh = ran
	return
	end
*==================================
	FUNCTION RAN3(IDUM)
C			IMPLICIT REAL*4(M)
C			PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG)
	PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1./MBIG)
	DIMENSION MA(55)
	DATA IFF /0/
	IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
		IFF=1
		MJ=MSEED-IABS(IDUM)
		MJ=MOD(MJ,MBIG)
		MA(55)=MJ
		MK=1
		DO I=1,54
			II=MOD(21*I,55)
			MA(II)=MK
			MK=MJ-MK
			IF(MK.LT.MZ)MK=MK+MBIG
			MJ=MA(II)
		ENDDO
		DO K=1,4
			DO I=1,55
				MA(I)=MA(I)-MA(1+MOD(I+30,55))
				IF(MA(I).LT.MZ)MA(I)=MA(I)+MBIG
			ENDDO
		ENDDO
		INEXT=0
		INEXTP=31
		IDUM=1
	ENDIF
	INEXT=INEXT+1
	IF(INEXT.EQ.56)INEXT=1
	INEXTP=INEXTP+1
	IF(INEXTP.EQ.56)INEXTP=1
	MJ=MA(INEXT)-MA(INEXTP)
	IF(MJ.LT.MZ)MJ=MJ+MBIG
	MA(INEXT)=MJ
	RAN3=MJ*FAC
	RETURN
	END
C*LOOKUP -- inverse linear interpolation
C+
	SUBROUTINE LOOKUP( NPIX, DATA, TARGET, ILO, IHI, PART )
*
* Inverse linear interpolator. Finds ILO, IHI, PART such that
*   TARGET = DATA(ILO) * (1.-PART) + DATA(IHI) * PART
*
* This version uses safe but slow method of one-pixel steps
*
* Input:
*	NPIX	= Number of data points
*	DATA	= Data values ( must be monotonic )
*	TARGET	= Target value
* Output:
*	ILO	= Low pixel
*	IHI	= Hi pixel
*	PART	= fraction (0-1) for interpolation
C--
* Nov 1985 KDH @ STScI
* 2001 Aug KDH @ St.And - single point, IBEG
	REAL*4 DATA(*)
	DATA NLAST,IHILAST,ILOLAST/1,1,1/

* Null array
      IF( NPIX.LE.0 ) THEN
	ILO = 0
	IHI = 0
	PART = 0.
	RETURN
      END IF

* Single point
      IF( NPIX .EQ. 1 ) THEN
	ILO = 1
	IHI = 1
	PART = 0.
	RETURN
      END IF

* Use previous location
      IF( NPIX .EQ. NLAST ) THEN
      IF( DATA(IHILAST) .GE. TARGET .AND.
     *	  DATA(ILOLAST) .LE. TARGET ) THEN
	IHI = IHILAST
	ILO = ILOLAST
	GOTO 50
      END IF
	IBEG = ILOLAST
      ELSE
	IBEG = ( 1 + NPIX ) / 2
      END IF

* Determine uphill direction

      IF( DATA(1) .LE. DATA(NPIX) ) THEN
	ILOEND = 1
	IHIEND = NPIX
	IUP = 1
      ELSE
	ILOEND = NPIX
	IHIEND = 1
	IUP = -1
      END IF

* Locate an uphill point

      DO I = IBEG, IHIEND, IUP
	IF( DATA(I) .GE. TARGET ) GOTO 10
      END DO
	I = IHIEND
   10	IHI = I

* Locate nearest downhill point

      DO I = IHI, ILOEND, -IUP
	IF( DATA(I) .LE. TARGET ) GOTO 20
      END DO
	I = ILOEND
   20	ILO = I

* Locate nearest uphill point

      DO I = ILO, IHI, IUP
	IF( DATA(I) .GE. TARGET ) GOTO 30
      END DO
	I = IHI
   30	IHI = I

* Compute fractional part

   50	PART = DATA(IHI) - DATA(ILO)
	IF( PART.NE.0. ) PART = ( TARGET - DATA(ILO) ) / PART

	NLAST = NPIX
	ILOLAST = ILO
	IHILAST = IHI

	RETURN
	END
	function terp1( n, x, i1, i2, p )
* linear interpolation
* 2001 Aug Keith Horne @ St-Andrews
	real*4 x(n)
	terp1 = x(i1) * (1.-p) + x(i2) * p
	return
	end
***********************************************
	SUBROUTINE UPCASE( TEXTIN, TEXTOUT )
* replacement for VMS library routine STR$UPCASE( TEXTOUT, TEXTIN )
	CHARACTER*(*) TEXTIN
	CHARACTER*(*) TEXTOUT
	TEXTOUT = TEXTIN
	JADD = ICHAR('A') - ICHAR('a')
      DO I=1,LEN1(TEXTOUT)
	JTEST = ICHAR(TEXTOUT(I:I) )
      IF( JTEST .GE. ICHAR('a') .AND. JTEST .LE. ICHAR('z') ) THEN
	TEXTOUT(I:I) = CHAR( JTEST + JADD )
      END IF
      END DO
	RETURN
	END
	SUBROUTINE INQ( MESSAGE, DEFAULT, REPLY )
*
* get character string from user
*
* input:
*	MESSAGE	C* message to be sent to user
*	DEFAULT C* default reply
* output:
*	REPLY	C* user reply to the message
*
* 1988 Mar KDH @ STScI
* 2001 Jul KDH @ St.And - avoid char*(*) concatenation
	CHARACTER*(*) MESSAGE
	CHARACTER*(*) DEFAULT
	CHARACTER*(*) REPLY
	CHARACTER*256 LOCAL

      IF( DEFAULT.EQ.' ' ) THEN
	WRITE( *, '(1X,A,A,$)' ) MESSAGE( :LEN1(MESSAGE) ), ' '
      ELSE
	WRITE( *, '(1X,A,A,A,A,$)' ) MESSAGE( :LEN1(MESSAGE) ),
     #		' [' , DEFAULT( :LEN1(DEFAULT) ) , '] '
      END IF
	READ( *, '(A)' ) LOCAL
      IF( LOCAL.EQ.' ' ) THEN
	REPLY = DEFAULT
      ELSE
	REPLY = LOCAL
      END IF
	RETURN
	END
*===========================================================
	SUBROUTINE INQI( MESSAGE, IDEF, IVAL )
*
* get integer from the user
*
* input:
*	MESSAGE	C* message
*	IDEF 	I4 default value
* output:
*	IVAL	I4 returned value
*
* 1988 Mar KDH @ STScI
* 2002 Jul KDH @ St.And - add end=10
*
	CHARACTER*(*) MESSAGE
	CHARACTER*20 NUMBER
   10	WRITE( NUMBER, * ) IDEF
	CALL WORD1( NUMBER, L1, L2 )
	CALL INQ( MESSAGE, NUMBER(L1:L2), NUMBER )
	READ( NUMBER, *, ERR=10, END=10 ) NEW
	IVAL = NEW
	RETURN
	END
*===========================================================
	SUBROUTINE INQR( MESSAGE, DEF, VAL )
*
* get real value from the user
*
* input:
*	MESSAGE	C* message
*	DEF 	R4 default value
* output:
*	VAL	R4 returned value
*
* 1988 Mar KDH @ STScI
* 2002 Jul KDH @ St.And - add end=10
*
	CHARACTER*(*) MESSAGE
	CHARACTER*20 NUMBER
   10	WRITE( NUMBER, * ) DEF
	CALL INQ( MESSAGE, NUMBER, NUMBER )
	READ( NUMBER, *, ERR=10, END=10 ) WANT
	VAL = WANT
	RETURN
	END
*===========================================================
	SUBROUTINE INQD( MESSAGE, DEF, VAL )
*
* get real*8 value from the user
*
* input:
*	MESSAGE	C* message
*	DEF 	R8 default value
* output:
*	VAL	R8 returned value
*
* 1988 Mar KDH @ STScI
* 2002 Jul KDH @ St.And - add end=10
* 2004 Feb KDH @ St.And - adapt from INQR
* 2012 Nov KDH @ St.And - char*20->char*30 number
	REAL*8 DEF, VAL, WANT
	CHARACTER*(*) MESSAGE
	CHARACTER*30 NUMBER
   10	WRITE( NUMBER, * ) DEF
	CALL INQ( MESSAGE, NUMBER, NUMBER )
	READ( NUMBER, *, ERR=10, END=10 ) WANT
	VAL = WANT
	RETURN
	END
*===========================================================
	SUBROUTINE INQ2R( NAME, R1, R2 )
* get 2 reals from the user
* input:
*	NAME	C* name
*	R1 	R4 default value
*	R2 	R4 default value
* output:
*	R1	R4 returned value
*	R2	R4 returned value
* 2002 Jul KDH @ St.And
* 2008 Jul KDH @ SAAO - read v1, then v1,v2
	CHARACTER*(*) NAME
	CHARACTER*256 REPLY
   10	WRITE(*,*) 'Old ', NAME( :LEN1(NAME) ), ' : ', R1, R2
	WRITE(*,'(A,A,A,$)') ' New ', NAME( :LEN1(NAME) ), ' : '
	READ( *, '(A)' ) REPLY
      IF( REPLY .NE. ' ' ) THEN
	READ( REPLY ,*, ERR=10, END=20 ) V1
	R1 = V1
	READ( REPLY ,*, ERR=10, END=20 ) V1, V2
	R2 = V2
      END IF
   20	RETURN
	END
*===========================================================
	SUBROUTINE INQ2D( NAME, R1, R2 )
* get 2 real*8 values from the user
* In:	NAME	C* name
*	R1 	R8 default value
*	R2 	R8 default value
* Out:	R1	R8 returned value
*	R2	R8 returned value
* 2002 Jul KDH @ St.And
* 2008 Jul KDH @ SAAO - read v1, then v1,v2
* 2020 Jul KDH @ ST.And - adapt from INQ2R
	REAL*8 R1, R2, V1, V2
	CHARACTER*(*) NAME
	CHARACTER*256 REPLY
   10	WRITE(*,*) 'Old ', NAME( :LEN1(NAME) ), ' : ', R1, R2
	WRITE(*,'(A,A,A,$)') ' New ', NAME( :LEN1(NAME) ), ' : '
	READ( *, '(A)' ) REPLY
      IF( REPLY .NE. ' ' ) THEN
	READ( REPLY ,*, ERR=10, END=20 ) V1
	R1 = V1
	READ( REPLY ,*, ERR=10, END=20 ) V1, V2
	R2 = V2
      END IF
   20	RETURN
	END
*===========================================================
	SUBROUTINE INQ3R( NAME, R1, R2, R3 )
* get 3 reals from the user
* input:
*	NAME	C* name
*	R1 	R4 default value
*	R2 	R4 default value
*	R3 	R4 default value
* output:
*	R1	R4 returned value
*	R2	R4 returned value
*	R3	R4 returned value
* 2002 Jul KDH @ St.And
* 2008 Jul KDH @ SAAO - read v1, then v1,v2
* 2008 Jul KDH @ StA - adapt from INQ2R
	CHARACTER*(*) NAME
	CHARACTER*256 REPLY
   10	WRITE(*,*) 'Old ', NAME( :LEN1(NAME) ), ' : ', R1, R2, R3
	WRITE(*,'(A,A,A,$)') ' New ', NAME( :LEN1(NAME) ), ' : '
	READ( *, '(A)' ) REPLY
      IF( REPLY .NE. ' ' ) THEN
	READ( REPLY ,*, ERR=10, END=20 ) V1
	R1 = V1
	READ( REPLY ,*, ERR=10, END=20 ) V1, V2
	R2 = V2
	READ( REPLY ,*, ERR=10, END=20 ) V1, V2, V3
	R3 = V3
      END IF
   20	RETURN
	END
*===========================================================
	SUBROUTINE INQ2I( NAME, I1, I2 )
* get 2 integers from the user
* input:
*	NAME	C* name
*	I1 	I4 default value
*	I2 	I4 default value
* output:
*	I1	I4 returned value
*	I2	I4 returned value
* 2002 Jul KDH @ St.And
* 2008 Jul KDH @ SAAO - read j1, then j1,j2
	CHARACTER*(*) NAME
	CHARACTER*256 REPLY
   10	WRITE(*,*) 'Old ', NAME( :LEN1(NAME) ), ' : ', I1, I2
	WRITE(*,'(A,A,A,$)') ' New ', NAME( :LEN1(NAME) ), ' : '
	READ( *, '(A)' ) REPLY
      IF( REPLY .NE. ' ' ) THEN
	read( reply, *, err=10, end=20 ) j1
	i1 = j1
	READ( REPLY ,*, ERR=10, END=20 ) J1, J2
	I2 = J2
      END IF
   20	RETURN
	END

************************************
	integer function len1( string )
* length up to last non-blank character
	character*(*) string
      do i = len( string ), 1, -1
      if( string( i:i ) .ne. ' ' ) then
	len1 = i
	return
      end if
      end do
	len1 = 0
	return
	end
C*WORD1 -- isolate first word of a text string
C+
	SUBROUTINE WORD1( STRING, L1, L2 )
*
* Isolates the first word of a character string 
*
* Input:
*	STRING	= character string
* Output:
*	L1,L2	= Inclusive index limits of first word
C--
* Jan 1986 Keith Horne @ STScI
* Jan 1989 KDH @ STScI - revise to allow BLANK or TAB separators
* Oct 2014 KDH @ StAnd - revise to allow comma separators
	CHARACTER*(*) STRING
	CHARACTER*1 TAB, BLANK, TEST, comma
	comma = ','
	TAB = CHAR( 9 )
	BLANK = ' '
	LMAX = LEN1( STRING )
      DO I1 = 1, LMAX
	TEST = STRING(I1:I1)
      IF(     TEST .NE. BLANK
     &	.AND. TEST .NE. TAB
     &	.and. test .ne. comma ) THEN
	L1 = I1
      DO I2 = L1 + 1, LMAX
	TEST = STRING(I2:I2)
      IF(    TEST .EQ. BLANK
     &	.OR. TEST .EQ. TAB
     &	.or. test .eq. comma ) THEN
	L2 = I2 - 1
	RETURN
      END IF
      END DO
	L2 = LMAX
	RETURN
      END IF
      END DO
	L1 = 0
	L2 = 0
	RETURN
	END
*************************
	subroutine append( str1, str2 )
* append two strings with a space in between
* 2003 Jun Keith Horne @ St.Andrews
	character*(*) str1, str2
	l1 = len1( str1 )
	l2 = len1( str2 )
      if( l1 .le. 0 ) then
	str1 = str2( : l2 )
      else
	str1 = str1( : l1 ) // ' ' // str2( : l2 )
      end if
	return
	end
****************************************
	subroutine append_i( str, i )
* append integer to string
* 2003 Jun Keith Horne @ St.Andrews
	character*(*) str
	character*40 number
	write( number, * ) i
	call word1( number, l1, l2 )
	call append( str, number( l1 : l2 ) )
	return
	end
****************************************
	subroutine append_f( str, x, fmt )
* append formatted real to string
* 2003 Jun Keith Horne @ St.Andrews
* 2008 Aug KDH @ St.A - add '*' format
	character*(*) str, fmt
	character*40 number
      if( fmt(1:1) .eq. '*' ) then
	write( number, * ) x
      else
	write( number, fmt ) x
      end if
	call word1( number, l1, l2 )
	call append( str, number( l1 : l2 ) )
	return
	end
****************************************
	subroutine append_d( str, x, fmt )
* append formatted real*8 to string
* 2004 Feb Keith Horne @ St.Andrews
	real*8 x
	character*(*) str, fmt
	character*40 number
	write( number, fmt ) x
	call word1( number, l1, l2 )
	call append( str, number( l1 : l2 ) )
	return
	end

****************************************
	subroutine append_n( title, val, idig )
* append val to title with idig significant digits
* 2008 Jun Keith Horne @ St.Andrews
* 2008 Jul KDH @ SAAO - idig
	character*(*) title
	character*20 fmt
	i = max( 2, idig )
	a= abs( val ) / 10**(i)
      if( a .eq. 0.0 ) then
	call append_f( title, val, '(f10.1)' )
      else if( a .lt. 9.9e-8) then
	fmt = '(1pe'
	call append_i( fmt, i + 7 )
	call append( fmt, '.' )
	call append_i( fmt, i - 1 )
	call append( fmt, ')' )
	call nospace( fmt )
	call append_f( title, val, fmt )
      else if( a .lt. 9.9e-7) then
	call append_f( title, val, '(f15.6)' )
      else if( a .lt. 9.9e-6) then
	call append_f( title, val, '(f15.5)' )
      else if( a .lt. 9.9e-5 ) then
	call append_f( title, val, '(f15.4)' )
      else if( a .lt. 9.9e-4 ) then
	call append_f( title, val, '(f15.3)' )
      else if( a .lt. 9.9e-3 ) then
	call append_f( title, val, '(f15.2)' )
      else if( a .lt. 0.099 ) then
	call append_f( title, val, '(f15.1)' )
      else
	call append_i( title, nint( val ) )
      end if
	return
	end

********************************************
        subroutine append_datsig( title, dat, sig, kfig )
* Input
*       title   c* string
*       dat     r4 data
*       sig     r4 sigma
*       kfig    i4 number of significant figures on sigma
* Output
*       title   c* string // data(sigma)
* 2004 Feb Keith Horne @ St-Andrews
* 2013 May KDH @ StA trap sig not positive
* 2017 Sep KDH @ StA revise for very small/large numbers
        character*(*) title
        real*8 hjd, big
	dmx = 1000
	dmn = 0.0001

      if( sig .le. 0. ) then
	call append_n( title, dat, max( 2, kfig ) )

      else if( abs(dat) .ge. dmn .and. abs( dat ) .le. dmx ) then
        hjd = dat
        call append_hjdsig( title, hjd, sig, kfig )

      else

* divide by a power of 10
        mpow = nint( log10( abs( dat * 2. ) ) )
        ipow = nint( alog10( sig * 2. ) )
        kpow = ipow - kfig

	kpow = int( alog10( max( abs(dat), sig ) ) ) - 1
        div = 10. ** kpow
	hjd = dat / div

	call append_hjdsig( title, hjd, sig/div, kfig )
* exponent
      if( kpow .ne. 0 ) then
	l = len1( title ) + 1
	title(l:l) = 'E'
	call append_i( title, kpow )
	call nospace( title( l : len1(title) ) )
      end if
      end if
        return
        end

********************************************
        subroutine append_hjdsig( title, hjd, sig, kfig )
* Input
*       title   c* string
*       hjd     r8 data
*       sig     r4 sigma
*       kfig    i4 number of significant figures on sigma
* Output
*       title   c* string // data(sigma)
* 2004 Feb Keith Horne @ St-Andrews
* 2012 Dec KDH @ StA omit E0, skip if sig not positive
* 2013 May KDH @ StA trap sig not positive
        character*(*) title
        real*8 hjd
        character*50 fmt

      if( sig .le. 0. ) then
	call append_n( title, real(hjd), max( 2, kfig ) )
	return
      end if

        l1 = len1( title )
        mpow = nint( log10( abs( hjd * 2. ) ) )
        ipow = nint( alog10( sig * 2. ) )
        kpow = ipow - kfig
        div = 10. ** kpow

      if( kpow .lt. 0 ) then
        fmt = '(f20.'
        call append_i( fmt, -kpow )
        call append( fmt, ')' )
        call nospace( fmt )
        call append_d( title, hjd, fmt )

      if( sig .gt. 0. ) then
        sigdiv = sig / div
        call append( title, '(' )
      if( sig .lt. 1. ) then
        call append_i( title, nint( sigdiv ) )
      else
        call append_f( title, sig, fmt )
      end if
        call append( title, ')' )
      end if

      else if( kpow .gt. 3 ) then
        lpow = mpow
        if( iabs( ipow ) .gt. iabs( mpow ) ) lpow = ipow
        div = 10. ** lpow
        fmt = '(f20.'
        call append_i( fmt, kfig )
        call append( fmt, ')' )
        call nospace( fmt )
        call append_d( title, hjd / div, fmt )

      if( sig .gt. 0. ) then
        call append( title, '(' )
        call append_f( title, sig / div, fmt )
        call append( title, ')' )
      if( lpow .ne. 0 ) then
        call append( title, 'E' )
        call append_i( title, lpow )
      end if
      end if

      else
        idiv = 10 ** kpow
        call append_i( title, nint( hjd / idiv ) * idiv )
      if( sig .gt. 0. ) then
        call append( title, '(' )
        call append_i( title, nint( sig / idiv ) * idiv )
        call append( title, ')' )
      end if

      end if
        call nospace( title(l1+2:) )
        return
        end

****************************************
	subroutine nospace( str )
* remove spaces from string
* 2003 Aug Keith Horne @ St.Andrews
	character*(*) str
	n = len1( str )
	if( n .lt. 1 ) return
	j = 0
      do i=1,n
      if( str( i : i ) .ne. ' ' ) then
	j = j + 1
	str( j : j ) = str( i : i )
      end if
      end do
	str = str( : j )
	return
	end
****************************************
	subroutine notail( str, chr )
* remove specified character from end of string
* 2003 Dec Keith Horne @ St.Andrews
	character*(*) str, chr
	n = len1( str )
	if( n .lt. 1 ) return
      do while( n .gt. 0 .and. str(n:n) .eq. chr(1:1) )
	n = n - 1
      end do
	str = str(:n)
	return
	end
*------------------------------
	subroutine noeqbl( title )
* change '= ' to '='
* 2006 Feb Keith Horne @ St.Andrews
	character*(*) title
    1	l = index( title, '= ' )
	if( l .le. 0 ) return
	title = title(:l) // title(l+2:)
	goto 1
	end
*------------------------------
	subroutine nogreek( title )
* change pgplot greeks to text
* 2007 Aug Keith Horne @ LaSilla
	character*(*) title
	parameter ( nfix = 15 )
	character*3 old( nfix )
	character*5 new( nfix )
	data old/
     &		'\ga',	'\gb',	'\gg',	'\gg',	'\ge'
     &	,	'\gh',	'\gy',	'\gm',	'\gn',	'\gp'
     &	,	'\gr',	'\gs',	'\gx',	'\gD',	'\gW'
     &	/
	data new/
     &		'alp',	'bet',	'gam',	'del',	'eps'
     &	,	'theta', 'eta',	'mu',	'nu',	'pi'
     &	,	'rho',	'sig',	'chi',	'Del',	'Omega'
     &	/

      do k=1,nfix
	i = 1
      do while( i .gt. 0 )
	i = index( title, old(k) )
      if( i .gt. 0 ) then
	n = len1( new(k) )
	title = title(:i-1) // new(k)(:n) // title(i+3:)
      end if
      end do
      end do

	return
	end
*------------------------------
	subroutine noupdn( title )
* change pgplot sub and superscripts to text
* 2007 Aug Keith Horne @ LaSilla
	character*(*) title
	m = 1
      do while( m .gt. 0 )
	i = index( title, '\u' )
	j = index( title, '\d' )
      if( i .gt. 0 .and. j .gt. i ) then
	title = title(:i-1) // '^' // title(i+2:j-1) // title(j+2:)
      else if( j .gt. 0 .and. i .gt. j ) then
	title = title(:j-1) // '_' // title(j+2:i-1) // title(i+2:)
      else
	m = 0
      end if
      end do
	return
	end

*=============================================
        subroutine append_kmb( string, big, n )
* append big with units K=1e3, M=1e6, B=1e9
* In:   string  c* character string
*       big     r4 large number
*       n       i4 number of digits of precision
* Out:  label   c* appended string
* 2020 Apr Keith Horne @ St.Andrews
        character*(*) string
        x = big
        a = abs( x )
      if( a .lt. 1.e4 ) then
        call append_i( string, nint( x ) )
      else if( a .lt. 1.e6 ) then
        call append_n( string, x / 1.e3, n )
        call append( string, 'K' )
      else if( a .lt. 1.e9 ) then
        call append_n( string, x / 1.e6, n )
        call append( string, 'M' )
      else
        call append_n( string, x / 1.e9, n )
        call append( string, 'B' )
      end if
        return
  	end
	subroutine fullpath( s1, s2 )
* evaluate unix path by translating environment variables
* 2001 Jul Keith Horne @ St.Andrews
	character*(*) s1, s2
	parameter ( maxchar = 256 )
	character*( maxchar ) t1, t2
* null string
	call word1( s1, l1, l2 )
      if( l1.lt.1 ) then
	s2 = s1
	return
      end if
* copy to a local string
	t1 = s1( l1:l2 )
      if( l2-l1+1 .gt. maxchar ) then
	write(*,*) '** fullpath truncation:', l1, l2, l2-l1+1, maxchar
	write(*,*) '"', s1( :len1(s1) ), '"'
	write(*,*) '"', t1( :len1(t1) ), '"'
      end if
* isolate the environment variable to be translated
	ldollar = index( t1, '$' )
      do while( ldollar .gt. 0 )
	l1 = ldollar + 1
	lslash = ldollar + index( t1( l1:len1(t1) ) // '/', '/' )
	l2 = lslash - 1
* translate the environment variable
	call getenv( t1( l1:l2 ), t2 )
      if( len1( t2 ) .ge. maxchar ) then
	write(*,*) '** fullpath possible truncation:'
	write(*,*) '"', t1( :len1(t1) ), '"'
	write(*,*) '"', t2( :len1(t2) ), '"'
      end if

	if( ldollar .gt. 1 ) t2 = t1(:ldollar-1) // t2(:len1(t2) )

	t1 = t2(:len1(t2) ) // t1(l2+1:)

	ldollar = index( t1, '$' )
      end do
* copy to output string
	s2 = t1
	l1 = len1( t1 )
	l2 = len1( s2 )
      if( l1 .gt. l2 ) then
	write(*,*) '** fullpath truncation:', l1, l2
	write(*,*) '"', t1(:l1), '"'
	write(*,*) '"', s2(:l2), '"'
      end if
	return
	end
************************
	subroutine spawn( task )
*
* spawn a task to the operating system
*
* 1994 Jan KDH @ Utrecht - dummy version
* 1994 Mar KDH @ Utrecht - call system( task )
* 1994 Mar KDH @ Utrecht - call system( spawn task )
* 2000 Nov KDH @ Austin - call system( task )
* 2001 Jul KDH @ St-And - linux g77 rejects char*(*) concatenations
* 2015 Aug KDH @ St-And - comment out external system
	character*(*) task

* external seems not to be needed and causes problems on some compilers
c	external system

* at Utrecht, this works:
c	write(*,*) ' call system( spawn ' // TASK(:LEN1(TASK)) // ' )'
c	write(*,*) ' '
c	CALL SYSTEM( 'spawn ' // TASK(:LEN1(TASK)) )
* Note: $HOME/bin/spawn is an executable script as follows :
*  #!/bin/csh
*  source $HOME/.login ; $*
* at Austin, this works:
c	write(*,*) ' call system( ' // task(:len1(task)) // ' )'
c	call system( task(:len1(task)) )

* on linux, this works
* Note: linux g77 rejects character*(*) concatenation
	l1 = len1( task )
	write(*,*) 'call system( ', task(:l1), ' ) '
	call system( task( :l1 ) )
	return
	end
********************************************
	SUBROUTINE QUICKSORT( N, X, KEY )
*
* QUICKSORT subroutine (Hoare)
* Produces index such that X(KEY(I)) ascends with I
*
* Input:
*	N	= number of items to be sorted
*	X(N)	= REAL*4 values to be sorted
* Output:
*	KEY(N)	= sort key such that X(KEY(I)) ascends with I
*
* Adapted by T.R.Marsh from a Mike Irwin routine.
* X is left unchanged.  Compatible arguments with
* SHELLSORT and BUBSORT.  About 3 times faster than
* SHELLSORT for more than 5000 points, roughly the
* same below 2000.
*
* July 1985 KDH @ STScI -- typed in and tested.
* 1999 Oct KDH @ St.Andrews - eliminate TYPE
* 2001 Jul KDH @ St.And -- use .eqv. when comparing logicals
* 2018 Feb KDH @ St.And X(*),KEY(*)=>X(N),KEY(N)
	REAL*4 X( N )
	INTEGER*4 KEY( N )
	PARAMETER ( MAXSTACK=300 )
	INTEGER*4 ISTACK(2,MAXSTACK)
	LOGICAL*1 L1, L2, L3, CARRYON
      DO I = 1, N
	KEY(I) = I
      END DO
	CARRYON = .TRUE.
	M = 9
	IF( N.LE.M ) GOTO 500
	IL = 1
	IR = N
	ISP = 0			! Initialize stack pointer
      DO WHILE( CARRYON )
	I = IL
	J = IR + 1
	IM = ( IL + IR ) / 2	! Find median for partition
	L1 = X( KEY(IM) ) .LT. X( KEY(IR) )
	L2 = X( KEY(IR) ) .LT. X( KEY(IL) )
	L3 = X( KEY(IL) ) .LT. X( KEY(IM) )
      IF( L2 .EQV. L3 ) THEN
	IQ = IL
      ELSE IF( L1 .EQV. L3 ) THEN
	IQ = IM
      ELSE
	IQ = IR
      END IF
	XTEST = X( KEY(IQ) )
	ITEST = KEY(IQ)
	KEY(IQ) = KEY(IL)
	KEY(IL) = ITEST
  100	I = I + 1
      DO WHILE( X( KEY(I) ) .LT. XTEST )
	I = I + 1
      END DO
	J = J - 1
      DO WHILE( XTEST .LT. X( KEY(J) ) )
	J = J - 1
      END DO
      IF( J.LE.I ) THEN
	ITEST = KEY(J)
	KEY(J) = KEY(IL)
	KEY(IL) = ITEST
      ELSE
	ITEST = KEY(J)
	KEY(J) = KEY(I)
	KEY(I) = ITEST
	GOTO 100
      END IF

* Update stack

      IF( IR-J.GE.J-IL .AND. J-IL.GT.M ) THEN
	ISP = ISP + 1
	IF( ISP.GT.MAXSTACK ) write(*,*)
     &		CHAR(7),'*** STACK OVERFLOW IN QUICKSORT.'
	ISTACK(1,ISP) = J + 1
	ISTACK(2,ISP) = IR
	IR = J - 1
      ELSE IF( J-IL.GT.IR-J .AND. IR-J.GT.M ) THEN
	ISP = ISP + 1
	IF( ISP.GT.MAXSTACK ) write(*,*)
     &		CHAR(7),'*** STACK OVERFLOW IN QUICKSORT.'
	ISTACK(1,ISP) = IL
	ISTACK(2,ISP) = J - 1
	IL = J+1
      ELSE IF( IR-J.GT.M .AND. M.GE.J-IL ) THEN
	IL = J + 1
      ELSE IF( J-IL.GT.M .AND. M.GE.IR-J ) THEN
	IR = J-1
      ELSE IF( ISP.GT.0 ) THEN
	IL = ISTACK(1,ISP)
	IR = ISTACK(2,ISP)
	ISP = ISP - 1
      ELSE
	CARRYON = .FALSE.
      END IF
      END DO

* Straight insertion sort

  500	IF( M .EQ. 1 ) RETURN
	INT = 1
	IFIN = N - INT
      DO II = 1, IFIN
	I = II
	J = I + INT
      IF( X( KEY(I) ) .GT. X( KEY(J) ) ) THEN
	XTEST = X( KEY(J) )
	ITEST = KEY(J) 
	KEY(J) = KEY(I)
	J = I
	I = I - INT
      DO WHILE( I .GT. 0 .AND. X( KEY(MAX(1,I)) ) .GT. XTEST )
	KEY(J) = KEY(I)
	J = I
	I = I - INT
      END DO
	KEY(J) = ITEST
      END IF
      END DO
	RETURN
	END
C*LOADLIST ... load a list of strings from a disk file
C+
	SUBROUTINE LOADLIST( SOURCE, NLIST, LIST )
*
* load a list of strings from a file
*
* Input:
*	SOURCE	C* "@filename" loads one or more strings from a file
*	NLIST	I4  >0 .. maximum allowed number of strings
*		    <0 .. requests the |NLIST|th string only
* Output:
*	NLIST	I4 total number of strings, or the sequence number
*		closest to the requested sequence number
*		(0 if failed)
*	LIST()	C* the list of strings
C--
* Jul 1988 KDH @ STScI - rename to LOADLIST
* Mar 1989 KDH @ STScI - use SOURCE and 'NONE'.
* Mar 1989 KDH @ STScI - default extension is .lis, then .dat
* May 1989 KDH @ STScI - NLIST<0 option to select specific sequence number
* 			- use ! to start comment
* Jul 1990 KDH @ STScI - use ! or # to start comment
* Aug 1991 KDH @ STScI - fix bug so !,# work when NLIST<0
* Sep 1999 KDH @ StAndrews - eliminate TYPE *
* 2001 Jul KDH @ St.Andrews - no READONLY :(, no char*(*) concatenation
* 2001 Aug KDH @ St.A - openfile
* 2010 Mar KDH @ LCO - avoid ERR=n on statement n
	CHARACTER*(*)	SOURCE
	CHARACTER*(*)	LIST(*)
	PARAMETER ( MAXCHR = 256 )
	CHARACTER*( MAXCHR ) FILE
	CHARACTER*10	REPLY

* check list dimension

	MLIST = NLIST
	NLIST = 0

      IF( MLIST.EQ.0 ) THEN
	write(*,*) '** dimension error ', MLIST
	GOTO 999
      END IF

* an empty list

	IF( LEN1( SOURCE ) .LE. 0 ) GOTO 1000
	CALL WORD1( SOURCE, L1, L2 )
	REPLY = SOURCE( L1:L2 )
	CALL UPCASE( REPLY, REPLY )
	IF( REPLY.EQ. 'NONE' ) GOTO 1000

* a single entry list rather than input from a disk file

	L0 = INDEX( SOURCE, '@' )
      IF( L0.LE.0 ) THEN
	NLIST = 1
	LIST(1) = SOURCE
	GOTO 1000
      END IF

* isolate name of the disk file from which the list is read

	CALL WORD1( SOURCE(L0+1:), L1, L2 )
	IF( L1.LT.0 .OR. L2.LT.L1 ) GOTO 999
	L1 = L0 + L1
	L2 = L0 + L2

* open disk file (default is .lis, then .dat )

	L3 = MIN( L2, MAX( L1, INDEX( SOURCE, ']' ) ) )
      IF( INDEX( SOURCE(L3:L2), '.' ) .GT.0 ) THEN
   25	CALL OPENFILE( 47, SOURCE(L1:L2), 'FORMATTED', 'OLD', IER )
	IF( IER .NE. 0 ) GOTO 90
      ELSE
	FILE = SOURCE( L1:L2 ) // '.lis'
	CALL OPENFILE( 47, FILE, 'FORMATTED', 'OLD', IER )
c	IF( IER .NE. 0 ) GOTO 25
      if( ier .ne. 0 ) then
	CALL OPENFILE( 47, SOURCE(L1:L2), 'FORMATTED', 'OLD', IER )
	IF( IER .NE. 0 ) GOTO 90
      end if
      END IF

* load entire list elements

      IF( MLIST.GT.0 ) THEN
	NLIST = 0
      DO I=1,MLIST
   30	continue
	READ( 47, '(A)', ERR=30, END=100 ) LIST(I)
	L1 = MAX( INDEX( LIST(I), '!' ), INDEX( LIST(I), '#' ) )
	IF( L1 .GT. 0 ) LIST(I) = LIST(I)(:L1-1)
	IF( LIST(I) .EQ. ' ' ) GOTO 30
	NLIST = NLIST + 1
      END DO

* find element closest to the requested sequence number

      ELSE
	MLIST = IABS( MLIST )
	NLIST = 0
   60	continue
	READ( 47, '(A)', ERR=60, END=100 ) LIST(1)
	L1 = MAX( INDEX( LIST(1), '!' ), INDEX( LIST(1), '#' ) )
	IF( L1 .GT. 0 ) LIST(1) = LIST(1)(:L1-1)
	IF( LIST(1) .EQ. ' ' ) GOTO 60
	NLIST = NLIST + 1
	IF( NLIST.LT.MLIST ) GOTO 60
      END IF

* close input file for normal return

  100	CLOSE( UNIT=47 )
	IF( NLIST.LT.1 ) GOTO 999

* normal return

 1000	RETURN

* error return

  999	write(*,*) '** LOADLIST aborted.', NLIST
	RETURN

* handle file open error

   90	write(*,*) '** error opening file ', SOURCE(:LEN1(SOURCE))
	GOTO 999
	END
C* DUMPLC ... Dump HJD and up to 7 arrays to free-format ascii file.
C+
      SUBROUTINE DUMPLCH( FILE, NH, H, NCOL, NROW, HJD
     &	, X1, X2, X3, X4, X5, X6, X7, X8, X9 )
*
* Dump HJD and up to 7 arrays to  free-format ascii disk file.
*
* Input:
*	FILE	C* file name or ' ' for interactive selection
*	NH	I4 Number of header lines
*	H(NH)	C* Header lines
*	NCOL	I4 Number of data arrays
*	NROW	I4 Number of data values in each array
*	HJD	R8 HJDs
*	X2-XNCOL R4(NROW) The data arrays
*
C--
* Oct 1985 by KDH at STScI
* Apr 1988 KDH @ STScI - generalized from DUMPMONGO
* May 1991 KDH @ STScI - install in XCAL
* Jan 1994 KDH @ Utrecht - maxcol=5, DUMP_DAT -> DUMP1D, improve error msgs
* Mar 1994 KDH @ Utrecht - entry points for 1-5 images
* May 1995 KDH @ St.And - fix bug causing no output when file specified
* 1999 Nov KDH @ St.And - maxcol=5->6
* 2000 Nov KDH @ Austin - open with type='unknown'      
* 2001 Aug KDH @ St.And - use OPENFILE
* 2002 Jul KDH @ St.And - add array 7
* 2002 Sep KDH @ St.And - DUMPLC adapted from DUMP1D
* 2002 Sep KDH @ St.And - FMT (linux g77 * format truncates real*8 )	
* 2007 Sep KDH @ Christchurch -- dumplc9 and 10
* 2018 May KDH @ St.And - DUMPLC -> DUMPLCH with header lines
* 2018 Sep KDH @ St.A - parameter( iout=93 )
	REAL*8 HJD(*)
	REAL*4 X1(*), X2(*), X3(*), X4(*), X5(*), X6(*), X7(*)
	REAL*4 X8(*), X9(*)
	CHARACTER*(*) H(*)
	CHARACTER*(*) FILE
	CHARACTER*256 RECORD
	CHARACTER*128 OFILE
	CHARACTER*128 FMT/ '( F18.8, 1P7G14.6 )' /
	CHARACTER*10 REPLY
	PARAMETER ( MAXCOL = 9 )
	parameter ( iout = 93 )
* entry points for dumb linux compiler
	GOTO 100
	ENTRY DUMPLC10H( FILE, NH, H, NCOL, NROW, HJD, X1, X2, X3, X4, X5
     &		, X6, X7, X8, X9 )
	GOTO 100
	ENTRY DUMPLC9H( FILE, NH, H, NCOL, NROW, HJD, X1, X2, X3, X4, X5
     &		, X6, X7, X8 )
	GOTO 100
	ENTRY DUMPLC8H( FILE, NH, H, NCOL, NROW, HJD, X1, X2, X3, X4, X5
     &		, X6, X7 )
	GOTO 100
	ENTRY DUMPLC7H( FILE, NH, H, NCOL, NROW, HJD, X1, X2, X3, X4, X5
     &		, X6 )
	GOTO 100
	ENTRY DUMPLC6H( FILE, NH, H, NCOL, NROW, HJD, X1, X2, X3, X4, X5 )
	GOTO 100
	ENTRY DUMPLC5H( FILE, NH, H, NCOL, NROW, HJD, X1, X2, X3, X4 )
	GOTO 100
	ENTRY DUMPLC4H( FILE, NH, H, NCOL, NROW, HJD, X1, X2, X3 )
	GOTO 100
	ENTRY DUMPLC3H( FILE, NH, H, NCOL, NROW, HJD, X1, X2 )
	GOTO 100
	ENTRY DUMPLC2H( FILE, NH, H, NCOL, NROW, HJD, X1 )
	GOTO 100
	ENTRY DUMPLC1H( FILE, NH, H, NCOL, NROW, HJD )

* check dimensions
  100 IF( NCOL .LT. 1 ) RETURN
      IF( NROW .LT. 1 ) RETURN
	NC = MIN( NCOL, MAXCOL )
      IF( NC .NE. NCOL ) THEN
        WRITE(*,*) '**DUMP1D: skipping columns', NC+1,'-', NCOL
      END IF
	OFILE = FILE
* get filename
   10 IF( OFILE .EQ. ' ' ) THEN
        WRITE(*,'(A,$)') ' Output filename : '
	READ(*,'(A)') OFILE
	REPLY = OFILE
	CALL UPCASE( REPLY, REPLY )
	IF( REPLY.EQ.' ') RETURN
	IF( REPLY.EQ.'NONE' ) RETURN
      END IF
* open file
	CALL OPENFILE( iout, OFILE, 'FORMATTED', 'UNKNOWN', IER )
	IF( IER .NE. 0 ) GOTO 99
	WRITE(*,*) 'File open ', OFILE(:LEN1(OFILE))

* write header lines
      IF( NH .GT. 0 ) THEN
	ITRUNK = 0
      DO I = 1, NH
	IH = I
	RECORD = H( I )

* report truncated header line
	LR = LEN1( RECORD )
	LH = LEN1( H( I ) )
      IF( LR .NE. LH ) THEN
	WRITE(*,*) '** DUMPLCH: HEADER LINE TRUNCATED FROM', LH, ' TO', LR
	WRITE(*,*) '** ', H(I)(:LH)
	ITRUNK = I
      END IF

* write header line
      IF( RECORD(1:1) .NE. '#' .AND. RECORD(1:1) .NE. '!' ) THEN
	RECORD = '# ' // RECORD( : LEN1( RECORD ) )
      END IF
	WRITE( iout, '(A)', ERR=97 ) RECORD( : LEN1( RECORD ) )
      END DO
      END IF

* write data
	NR = 0
      IF( NC .EQ. 1 ) THEN
      DO I=1,NROW
	WRITE( iout, FMT, ERR=98 ) HJD(I)
	NR = I
      END DO

      ELSE IF( NC .EQ. 2 ) THEN
      DO I=1,NROW
	WRITE( iout, FMT, ERR=98 ) HJD(I), X1(I)
	NR = I
      END DO

      ELSE IF( NC .EQ. 3 ) THEN
      DO I=1,NROW
	WRITE( iout, FMT, ERR=98 ) HJD(I), X1(I),X2(I)
	NR = I
      END DO

      ELSE IF( NC.EQ.4 ) THEN
      DO I=1,NROW
	WRITE( iout, FMT, ERR=98 ) HJD(I), X1(I), X2(I), X3(I)
	NR = I
      END DO

      ELSE IF( NC .EQ. 5 ) THEN
      DO I=1,NROW
	WRITE( iout, FMT, ERR=98 ) HJD(I), X1(I), X2(I), X3(I), X4(I)
	NR = I
      END DO

      ELSE IF( NC .EQ. 6 ) THEN
      DO I=1,NROW
	WRITE( iout, FMT, ERR=98 ) HJD(I),
     &	X1(I), X2(I), X3(I), X4(I), X5(I)
	NR = I
      END DO

      ELSE IF( NC .EQ. 7 ) THEN
      DO I=1,NROW
	WRITE( iout, FMT, ERR=98 ) HJD(I),
     &	X1(I), X2(I), X3(I), X4(I), X5(I), X6(I)
	NR = I
      END DO

      ELSE IF( NC .EQ. 8 ) THEN
      DO I=1,NROW
	WRITE( iout, FMT, ERR=98 ) HJD(I),
     &	X1(I), X2(I), X3(I), X4(I), X5(I), X6(I), X7(I)
	NR = I
      END DO

      ELSE IF( NC .EQ. 9 ) THEN
      DO I=1,NROW
	WRITE( iout, FMT, ERR=98 ) HJD(I),
     &	X1(I), X2(I), X3(I), X4(I), X5(I), X6(I), X7(I), X8(I)
	NR = I
      END DO

      ELSE IF( NC .EQ. 10 ) THEN
      DO I=1,NROW
	WRITE( iout, FMT, ERR=98 ) HJD(I)
     &	, X1(I), X2(I), X3(I), X4(I), X5(I)
     &	, X6(I), X7(I), X8(I), X9(I)
	NR = I
      END DO

      END IF

* close file
   50 WRITE(*,*) 'Dump columns:', NC, ' rows:', NR, ' of', NROW
	CLOSE( UNIT=iout )
* normal return
 1000 RETURN


* header write error
   97 WRITE(*,*) '**DUMPLCH: error writing header to file', ofile(:len1(ofile))
	goto 50

* file write error
 98   WRITE(*,*) '**DUMPLCH: error writing to file ', OFILE(:LEN1(OFILE))
	GOTO 50
       
* file open error
 99   WRITE(*,*) '**DUMPLCH: error opening file ', OFILE(:LEN1(OFILE))
	OFILE = ' '
	GOTO 10
	END
	subroutine openfile( iunit, name, form, status, ier )
* Input:
*	iunit	i4 unit number
*	name	c* file name
*	form	c* 'formatted' or 'unformatted'
*	status	c* 'old', 'new', ...
* Output:
*	ier	i4 0 if successful
* 2001 Jul Keith Horne @ St.Andrews -- use FULLPATH
* 2001 Aug KDH @ St-And -- silence
* 2001 Aug KDH @ St-And -- default extension .dat

	character*(*)	name, form, status
	character*256	path
	logical silence
	common/silence/ silence

* get full path
	call fullpath( name, path )
      if( .not. silence ) then
      if( path .ne. name ) then
	write(*,*) 'Path( ', name(:len1(name)), ' ) = ', 
     &		path(:len1(path)) 
      end if
      end if

* open file
	open( unit=iunit, file=path, form=form, status=status, err=2 )
	goto 10

* try again with default extension
    2	path = path( :len1( path ) ) // '.dat'
	open( unit=iunit, file=path, form=form, status=status, err=3 )
	goto 10

* normal return

   10	ier = 0
      if( .not. silence ) then
	write(*,*) 'Open(', iunit, ', ', path(:len1(path)),
     &	', ', form(:len1(form)), ', ', status(:len1(status)), ' )' 
      end if
	return

* error return

    3	ier = 1
	write(*,*) '** OPENFILE FAILED :( '
	write(*,*) '** Path( ', name(:len1(name)), ' ) = ', 
     &		path(:len1(path)) 
	return
	end
C*LOADLC .. load multi-column data from ascii file - 1st column real*8
C+
      SUBROUTINE LOADLC( NAME, NCOL, NROW,
     #	HJD, C2, C3, C4, C5, C6, C7, C8, C9,
     #	C10, C11, C12, C13, C14, C15, C16, C17, C18, C19,
     #	C20, C21, C22, C23, C24, C25, C26, C27, C28, C29,
     #	C30, C31, C32, C33, C34, C35, C36, C37, C38, C39,
     #	C40, C41, C42, C43, C44, C45, C46, C47, C48, C49,
     #	C50, C51, C52, C53, C54, C55, C56, C57, C58, C59,
     #	C60, C61, C62, C63, C64, C65, C66, C67, C68, C69,
     #	C70, C71, C72, C73, C74, C75, C76, C77, C78, C79,
     #	C80, C81, C82, C83, C84, C85, C86, C87, C88, C89,
     #	C90, C91, C92, C93, C94, C95, C96, C97, C98, C99,
     #	C100, C101, C102, C103, C104, C105, C106, C107, C108, C109,
     #	C110, C111, C112, C113, C114, C115, C116, C117, C118, C119,
     #	C120, C121, C122, C123, C124, C125, C126, C127, C128, C129,
     #	C130, C131, C132, C133, C134, C135, C136, C137, C138, C139,
     #	C140 )
*
* Loads up to 140 columns of data from a free-format ASCII data file
* Use as follows:
*       parameter ( maxdat = 10000 )
*       real*4 c1(maxdat), c2(maxdat), c3(maxdat)
*       character*256 file
*       file = 'file.dat'
*       ndat = maxdat
*       call load1d( file, 3, ndat, c1, c2, c3 )
*
*  Input:
*	NAME	= C* Name of disk file, ' ' for interactive selection.
*	NCOL	= I4 Number of columns to be read (if NCOL > 0)
*		     (if NCOL < 0, then only column -NCOL is read)
*	NROW	= I4 Maximum allowed number of rows
*  Output:
*	NAME	= C* Name of disk file (only if selected interactively)
*	NROW	= I4 Number of data points (0 if none loaded)
*	C1(NROW)= R8 first column (REAL*8 TO HOLD HJD)
*	C2(NROW)= R4 second column
*	  ...
*	C140(NROW)= R4 140th column
C--
* May 1987 Keith Horne @ STScI
* Sep 1987 KDH - add NCOL<0 option to read specified column.
* May 1988 KDH - add final read to pick up eof after reading all rows.
* Jul 1988 KDH @ STScI - change name from LOAD_ASCII to LOAD_DAT
* Sep 1991 KDH @ STScI - skip words when COL < 0.
* Jan 1994 KDH @ Utrecht - split LOAD_DAT to LOAD1D,LOAD1D10,...
* Mar 1994 KDH @ Utrecht - entry points LOAD1DX, X=1...20
* Mar 1998 KDH @ St.Andrews - MAXCOL=20->140
* 1999 Oct KDH @ St.And - change OPEN for unix compatibility
* 2001 Jul KDH @ St.And -- no // of character*(*)
* 2001 Aug KDH @ St.And -- use OPENFILE
* 2002 Mar KDH @ St.And -- convert LOAD1D to LOADLC
* 2006 Jul KDH @ StA -- ENTRY LOADPLC for planet format data
* 2009 Jul KDH @ StA -- C(1) -> C(*)
* 2010 Mar KDH @ LCO - continue to avoid err=n branch to same line
	REAL*8 C1(*)
	REAL*4        C2(*), C3(*), C4(*), C5(*)
	REAL*4 C6(*), C7(*), C8(*), C9(*), C10(*)
	REAL*4 C11(*), C12(*), C13(*), C14(*), C15(*)
	REAL*4 C16(*), C17(*), C18(*), C19(*), C20(*)
	REAL*4 C21(*), C22(*), C23(*), C24(*), C25(*)
	REAL*4 C26(*), C27(*), C28(*), C29(*), C30(*)
	REAL*4 C31(*), C32(*), C33(*), C34(*), C35(*)
	REAL*4 C36(*), C37(*), C38(*), C39(*), C40(*)
	REAL*4 C41(*), C42(*), C43(*), C44(*), C45(*)
	REAL*4 C46(*), C47(*), C48(*), C49(*)
	REAL*4 C50(*), C51(*), C52(*), C53(*), C54(*)
	REAL*4 C55(*), C56(*), C57(*), C58(*), C59(*)
	REAL*4 C60(*), C61(*), C62(*), C63(*), C64(*)
	REAL*4 C65(*), C66(*), C67(*), C68(*), C69(*)
	REAL*4 C70(*), C71(*), C72(*), C73(*), C74(*)
	REAL*4 C75(*), C76(*), C77(*), C78(*), C79(*)
	REAL*4 C80(*), C81(*), C82(*), C83(*), C84(*)
	REAL*4 C85(*), C86(*), C87(*), C88(*), C89(*)
	REAL*4 C90(*), C91(*), C92(*), C93(*), C94(*)
	REAL*4 C95(*), C96(*), C97(*), C98(*), C99(*)
	REAL*4 C100(*), C101(*), C102(*), C103(*), C104(*)
	REAL*4 C105(*), C106(*), C107(*), C108(*), C109(*)
	REAL*4 C110(*), C111(*), C112(*), C113(*), C114(*)
	REAL*4 C115(*), C116(*), C117(*), C118(*), C119(*)
	REAL*4 C120(*), C121(*), C122(*), C123(*), C124(*)
	REAL*4 C125(*), C126(*), C127(*), C128(*), C129(*)
	REAL*4 C130(*), C131(*), C132(*), C133(*), C134(*)
	REAL*4 C135(*), C136(*), C137(*), C138(*), C139(*)
	REAL*4 C140(*)
	CHARACTER*(*) NAME
	CHARACTER*256 FILE	! file name
	PARAMETER (INUNIT=49)	! input unit number
	PARAMETER (MAXCOL=140)	! maximum number of columns that can be read
	CHARACTER*(15*MAXCOL) RECORD ! buffer into which data are read
	REAL*8 ROW( MAXCOL )


* entry points since UNIX requires number of arguments to match exactly

	ENTRY LOADLC20( NAME, NCOL, NROW, C1, C2, C3, C4, C5,
     #	C6, C7, C8, C9, C10, C11, C12, C13, C14, C15,
     #	C16, C17, C18, C19, C20 )
        GOTO 100

	ENTRY LOADLC19( NAME, NCOL, NROW, C1, C2, C3, C4, C5,
     #	C6, C7, C8, C9, C10, C11, C12, C13, C14, C15,
     #	C16, C17, C18, C19 )
        GOTO 100

	ENTRY LOADLC18( NAME, NCOL, NROW, C1, C2, C3, C4, C5,
     #	C6, C7, C8, C9, C10, C11, C12, C13, C14, C15,
     #	C16, C17, C18 )
        GOTO 100

	ENTRY LOADLC17( NAME, NCOL, NROW, C1, C2, C3, C4, C5,
     #	C6, C7, C8, C9, C10, C11, C12, C13, C14, C15,
     #	C16, C17 )
        GOTO 100

	ENTRY LOADLC16( NAME, NCOL, NROW, C1, C2, C3, C4, C5,
     #	C6, C7, C8, C9, C10, C11, C12, C13, C14, C15,
     #	C16 )
        GOTO 100

        ENTRY LOADLC15( NAME, NCOL, NROW, C1, C2, C3, C4, C5,
     #	C6, C7, C8, C9, C10, C11, C12, C13, C14, C15 )
        GOTO 100

        ENTRY LOADLC14( NAME, NCOL, NROW, C1, C2, C3, C4, C5,
     #	C6, C7, C8, C9, C10, C11, C12, C13, C14 )
        GOTO 100

        ENTRY LOADLC13( NAME, NCOL, NROW, C1, C2, C3, C4, C5,
     #	C6, C7, C8, C9, C10, C11, C12, C13 )
        GOTO 100

        ENTRY LOADLC12( NAME, NCOL, NROW, C1, C2, C3, C4, C5,
     #	C6, C7, C8, C9, C10, C11, C12 )
        GOTO 100

        ENTRY LOADLC11( NAME, NCOL, NROW, C1, C2, C3, C4, C5,
     #	C6, C7, C8, C9, C10, C11 )
        GOTO 100

        ENTRY LOADLC10( NAME, NCOL, NROW, C1, C2, C3, C4, C5,
     #	C6, C7, C8, C9, C10 )
        GOTO 100

        ENTRY LOADLC9( NAME, NCOL, NROW, C1, C2, C3, C4, C5,
     #	C6, C7, C8, C9 )
        GOTO 100

        ENTRY LOADLC8( NAME, NCOL, NROW, C1, C2, C3, C4, C5,
     #	C6, C7, C8 )
        GOTO 100
        ENTRY LOADLC8T( NAME, NCOL, NHJD, NROW, C1, C2, C3, C4, C5,
     #	C6, C7, C8 )
        GOTO 101

        ENTRY LOADLC7( NAME, NCOL, NROW, C1, C2, C3, C4, C5,
     #	C6, C7 )
        GOTO 100
        ENTRY LOADLC7T( NAME, NCOL, NHJD, NROW, C1, C2, C3, C4, C5,
     #	C6, C7 )
        GOTO 101

        ENTRY LOADLC6( NAME, NCOL, NROW, C1, C2, C3, C4, C5,
     #	C6 )
        GOTO 100
        ENTRY LOADLC6T( NAME, NCOL, NHJD, NROW, C1, C2, C3, C4, C5,
     #	C6 )
        GOTO 101

        ENTRY LOADLC5( NAME, NCOL, NROW, C1, C2, C3, C4, C5 )
        GOTO 100
        ENTRY LOADLC5T( NAME, NCOL, NHJD, NROW, C1, C2, C3, C4, C5 )
        GOTO 101

        ENTRY LOADLC4( NAME, NCOL, NROW, C1, C2, C3, C4 )
        GOTO 100
        ENTRY LOADLC4T( NAME, NCOL, NHJD, NROW, C1, C2, C3, C4 )
        GOTO 101

        ENTRY LOADLC3( NAME, NCOL, NROW, C1, C2, C3 )
        GOTO 100
        ENTRY LOADLC3T( NAME, NCOL, NHJD, NROW, C1, C2, C3 )
        GOTO 101

        ENTRY LOADLC2( NAME, NCOL, NROW, C1, C2 )
        GOTO 100
        ENTRY LOADLC2T( NAME, NCOL, NHJD, NROW, C1, C2 )
        GOTO 101

        ENTRY LOADLC1( NAME, NCOL, NROW, C1 )
        GOTO 100

  100	NT = 1
	GOTO 102
  101	NT = NHJD
  102	NC = NCOL
	write(*,*) 'NCOL', NC, ' HJD COL', NT, ' NROW', nrow

* trap errors
      IF( NROW.LE.0. OR. NC.EQ.0 ) THEN
	WRITE(*,*) '** INVALID DIMENSIONS.', NROW, NC
	GOTO 999
      ELSE IF( NC.GT.MAXCOL ) THEN
	WRITE(*,*) '** WARNING: Too many columns for LOADLC.'
	WRITE(*,*) '** Columns', MAXCOL, ' thru', NC, ' skipped.'
      END IF

	MXT = MIN( NC, 8 )
      IF( NT .GT. MXT ) THEN
	WRITE(*,*) '** ERROR: NT=', NT, ' > MAX =', MXT
	GOTO 999
      END IF

	MAXROW = NROW

* get name of input file

	FILE = NAME
c        WRITE(*,*) 'FILE"', NAME( :LEN1(NAME) ), '"'
      IF( FILE.EQ.' ' ) THEN
	WRITE(*,'(A,I3,A)')
     #	' Enter name of', IABS(NC),'-column free-format ASCII file :'
	READ(*,'(A)') FILE
	CALL WORD1( FILE, L1, L2 )
	IF( L1 .LE. 0 ) GOTO 99
	NAME = FILE( :LEN1( FILE ) )
      END IF

* open input disk file

	CALL OPENFILE( INUNIT, FILE, 'FORMATTED', 'OLD', IER )
	IF( IER .NE. 0 ) GOTO 99

* read desired column by skipping unwanted columns that precede it.

      IF( NC .LE. -2 ) THEN
	NSKIP = IABS(NC)-1
      DO I=1,MAXROW
   51	continue
	READ( INUNIT, '(A)', ERR=51, END=100 ) RECORD
      DO J=1,NSKIP
	CALL WORD1( RECORD, L1, L2 )
	IF( L1.LT.1 ) GOTO 51
	RECORD = RECORD(L2+1:)
      END DO
	READ( RECORD, *, ERR=51, END=51 ) C1(I)
      END DO

* read desired consecutive columns of data from the file

      ELSE IF( IABS(NC) .EQ. 1 ) THEN
      DO I=1,MAXROW
    1	continue
	READ( INUNIT, '(A)', ERR=1, END=200 ) RECORD
	READ( RECORD, *, ERR=1, END=1 ) C1(I)
c	WRITE(*,*) nc, nt, i, C1(I)
      END DO

* 2 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 2 ) THEN
      DO I=1,MAXROW
    2	continue
	READ( INUNIT, '(A)', ERR=2, END=200 ) RECORD
      IF( NT .EQ. 1 ) THEN
	READ( RECORD, *, ERR=2, END=2 ) C1(I), C2(I)
      ELSE IF( NT .EQ. 2 ) THEN
	READ( RECORD, *, ERR=2, END=2 ) C2(I), C1(I)
      END IF
c	WRITE(*,*) nc, nt, i, C1(I), C2(i)
      END DO

* 3 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 3 ) THEN
      DO I=1,MAXROW
    3	continue
	READ( INUNIT, '(A)', ERR=3, END=200 ) RECORD
      IF( NT .EQ. 1 ) THEN
	READ( RECORD, *, ERR=3, END=3 ) C1(I), C2(I), C3(I)
      ELSE IF( NT .EQ. 2 ) THEN
	READ( RECORD, *, ERR=3, END=3 ) C2(I), C1(I), C3(I)
      ELSE IF( NT .EQ. 3 ) THEN
	READ( RECORD, *, ERR=3, END=3 ) C2(I), C3(I), C1(I)
      END IF
c	WRITE(*,*) nc, nt, i, C1(I), C2(i), C3(i)
      END DO

* 4 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 4 ) THEN
      DO I=1,MAXROW
    4	continue
	READ( INUNIT, '(A)', ERR=4, END=200 ) RECORD
      IF( NT .EQ. 1 ) THEN
	READ( RECORD, *, ERR=4, END=4 ) C1(I), C2(I), C3(I), C4(I)
      ELSE IF( NT .EQ. 2 ) THEN
	READ( RECORD, *, ERR=4, END=4 ) C2(I), C1(I), C3(I), C4(I)
      ELSE IF( NT .EQ. 3 ) THEN
	READ( RECORD, *, ERR=4, END=4 ) C2(I), C3(I), C1(I), C4(I)
      ELSE IF( NT .EQ. 4 ) THEN
	READ( RECORD, *, ERR=4, END=4 ) C2(I), C3(I), C4(I), C1(I)
      END IF
c	WRITE(*,*) nc, nt, i, C1(I), C2(i), C3(i), C4(i)
      END DO

* 5 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 5 ) THEN
      DO I=1,MAXROW
    5	continue
	READ( INUNIT, '(A)', ERR=5, END=200 ) RECORD
      IF( NT .EQ. 1 ) THEN
	READ( RECORD, *, ERR=5, END=5 )
     #  C1(I), C2(I), C3(I), C4(I), C5(I)
      ELSE IF( NT .EQ. 2 ) THEN
	READ( RECORD, *, ERR=5, END=5 ) 
     #	C2(I), C1(I), C3(I), C4(I), C5(I)
      ELSE IF( NT .EQ. 3 ) THEN
	READ( RECORD, *, ERR=5, END=5 ) 
     #	C2(I), C3(I), C1(I), C4(I), C5(I)
      ELSE IF( NT .EQ. 4 ) THEN
	READ( RECORD, *, ERR=5, END=5 ) 
     #	C2(I), C3(I), C4(I), C1(I), C5(I)
      ELSE IF( NT .EQ. 5 ) THEN
	READ( RECORD, *, ERR=5, END=5 ) 
     #	C2(I), C3(I), C4(I), C5(I), C1(I)
      END IF

c	WRITE(*,*) nc, nt, i, C1(I), C2(i), C3(i), C4(i), C5(i)
      END DO

* 6 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 6 ) THEN
      DO I=1,MAXROW
    6	continue
	READ( INUNIT, '(A)', ERR=6, END=200 ) RECORD
      IF( NT .EQ. 1 ) THEN
        READ( RECORD, *, ERR=6, END=6 ) 
     #	C1(I), C2(I), C3(I), C4(I), C5(I), C6(I)
      ELSE IF( NT .EQ. 2 ) THEN
        READ( RECORD, *, ERR=6, END=6 ) 
     #	C2(I), C1(I), C3(I), C4(I), C5(I), C6(I)
      ELSE IF( NT .EQ. 3 ) THEN
        READ( RECORD, *, ERR=6, END=6 ) 
     #	C2(I), C3(I), C1(I), C4(I), C5(I), C6(I)
      ELSE IF( NT .EQ. 4 ) THEN
        READ( RECORD, *, ERR=6, END=6 ) 
     #	C2(I), C3(I), C4(I), C1(I), C5(I), C6(I)
      ELSE IF( NT .EQ. 5 ) THEN
        READ( RECORD, *, ERR=6, END=6 ) 
     #	C2(I), C3(I), C4(I), C5(I), C1(I), C6(I)
      ELSE IF( NT .EQ. 6 ) THEN
        READ( RECORD, *, ERR=6, END=6 ) 
     #	C2(I), C3(I), C4(I), C5(I), C6(I), C1(I)
      END IF
c	WRITE(*,*) nc, nt, i, C1(I), C2(i), C3(i), C4(i), C5(i)
c     &	, C6(i)
      END DO

* 7 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 7 ) THEN
      DO I=1,MAXROW
    7	continue
	READ( INUNIT, '(A)', ERR=7, END=200 ) RECORD
      IF( NT  .EQ. 1 ) THEN
        READ( RECORD, *, ERR=7, END=7 ) 
     #	C1(I), C2(I), C3(I), C4(I), C5(I), C6(I), C7(I)
      ELSE IF( NT .EQ. 2 ) THEN
        READ( RECORD, *, ERR=7, END=7 ) 
     #	C2(I), C1(I), C3(I), C4(I), C5(I), C6(I), C7(I)
      ELSE IF( NT .EQ. 3 ) THEN
        READ( RECORD, *, ERR=7, END=7 ) 
     #	C2(I), C3(I), C1(I), C4(I), C5(I), C6(I), C7(I)
      ELSE IF( NT .EQ. 4 ) THEN
        READ( RECORD, *, ERR=7, END=7 ) 
     #	C2(I), C3(I), C4(I), C1(I), C5(I), C6(I), C7(I)
      ELSE IF( NT .EQ. 5 ) THEN
        READ( RECORD, *, ERR=7, END=7 ) 
     #	C2(I), C3(I), C4(I), C5(I), C1(I), C6(I), C7(I)
      ELSE IF( NT .EQ. 6 ) THEN
        READ( RECORD, *, ERR=7, END=7 ) 
     #	C2(I), C3(I), C4(I), C5(I), C6(I), C1(I), C7(I)
      ELSE IF( NT .EQ. 7 ) THEN
        READ( RECORD, *, ERR=7, END=7 ) 
     #	C2(I), C3(I), C4(I), C5(I), C6(I), C7(I), C1(I)
      END IF

c	WRITE(*,*) nc, nt, i, C1(I), C2(i), C3(i), C4(i), C5(i)
c     &	, C6(i), C7(i)
      END DO

* 8 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 8 ) THEN
      DO I=1,MAXROW
    8	continue
	READ( INUNIT, '(A)', ERR=8, END=200 ) RECORD
c	write(*,*) nc, nt, i, ' ', record(:len1(record))
      IF( NT .EQ. 1 ) THEN
        READ( RECORD, *, ERR=8, END=8 ) 
     #	C1(I), C2(I), C3(I), C4(I), C5(I), C6(I), C7(I), C8(I)
      ELSE IF( NT .EQ. 2 ) THEN
        READ( RECORD, *, ERR=8, END=8 ) 
     #	C2(I), C1(I), C3(I), C4(I), C5(I), C6(I), C7(I), C8(I)
      ELSE IF( NT .EQ. 3 ) THEN
        READ( RECORD, *, ERR=8, END=8 ) 
     #	C2(I), C3(I), C1(I), C4(I), C5(I), C6(I), C7(I), C8(I)
      ELSE IF( NT .EQ. 4 ) THEN
        READ( RECORD, *, ERR=8, END=8 ) 
     #	C2(I), C3(I), C4(I), C1(I), C5(I), C6(I), C7(I), C8(I)
      ELSE IF( NT .EQ. 5 ) THEN
        READ( RECORD, *, ERR=8, END=8 ) 
     #	C2(I), C3(I), C4(I), C5(I), C1(I), C6(I), C7(I), C8(I)
      ELSE IF( NT .EQ. 6 ) THEN
        READ( RECORD, *, ERR=8, END=8 ) 
     #	C2(I), C3(I), C4(I), C5(I), C6(I), C1(I), C7(I), C8(I)
      ELSE IF( NT .EQ. 7 ) THEN
        READ( RECORD, *, ERR=8, END=8 ) 
     #	C2(I), C3(I), C4(I), C5(I), C6(I), C7(I), C1(I), C8(I)
      ELSE IF( NT .EQ. 8 ) THEN
        READ( RECORD, *, ERR=8, END=8 ) 
     #	C2(I), C3(I), C4(I), C5(I), C6(I), C7(I), C8(I), C1(I)
      END IF
c	WRITE(*,*) nc, nt, i, C1(I), C2(i), C3(i), C4(i), C5(i)
c     &	, C6(i), C7(i), C8(i)
      END DO

* 9 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 9 ) THEN
      DO I=1,MAXROW
    9   continue
	READ( INUNIT, '(A)', ERR=9, END=200 ) RECORD
	READ( RECORD, *, ERR=9, END=9 ) 
     #	C1(I), C2(I), C3(I), C4(I), C5(I),
     #	C6(I), C7(I), C8(I), C9(I)
      END DO

* 10 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 10 ) THEN
      DO I=1,MAXROW
   10   continue
	READ( INUNIT, '(A)', ERR=10, END=200 ) RECORD
	READ( RECORD, *, ERR=10, END=10 ) 
     #	C1(I), C2(I), C3(I), C4(I), C5(I),
     #	C6(I), C7(I), C8(I), C9(I), C10(I)
      END DO

* 11 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 11 ) THEN
      DO I=1,MAXROW
   11   continue
	READ( INUNIT, '(A)', ERR=11, END=200 ) RECORD
	READ( RECORD, *, ERR=11, END=11 ) 
     #	C1(I), C2(I), C3(I), C4(I), C5(I),
     #	C6(I), C7(I), C8(I), C9(I), C10(I),
     #	C11(I)
      END DO

* 12 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 12 ) THEN
      DO I=1,MAXROW
   12	continue
	READ( INUNIT, '(A)', ERR=12, END=200 ) RECORD
	READ( RECORD, *, ERR=12, END=12 ) 
     #	C1(I), C2(I), C3(I), C4(I), C5(I),
     #	C6(I), C7(I), C8(I), C9(I), C10(I),
     #	C11(I), C12(I)
      END DO

* 13 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 13 ) THEN
      DO I=1,MAXROW
   13   continue
	READ( INUNIT, '(A)', ERR=13, END=200 ) RECORD
	READ( RECORD, *, ERR=13, END=13 ) 
     #	C1(I), C2(I), C3(I), C4(I), C5(I),
     #	C6(I), C7(I), C8(I), C9(I), C10(I),
     #	C11(I), C12(I), C13(I)
      END DO

* 14 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 14 ) THEN
      DO I=1,MAXROW
   14   continue
	READ( INUNIT, '(A)', ERR=14, END=200 ) RECORD
	READ( RECORD, *, ERR=14, END=14 ) 
     #	C1(I), C2(I), C3(I), C4(I), C5(I),
     #	C6(I), C7(I), C8(I), C9(I), C10(I),
     #	C11(I), C12(I), C13(I), C14(I)
      END DO

* 15 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 15 ) THEN
      DO I=1,MAXROW
   15   continue
	READ( INUNIT, '(A)', ERR=15, END=200 ) RECORD
	READ( RECORD, *, ERR=15, END=15 ) 
     #	C1(I), C2(I), C3(I), C4(I), C5(I),
     #	C6(I), C7(I), C8(I), C9(I), C10(I),
     #	C11(I), C12(I), C13(I), C14(I), C15(I)
      END DO

* 16 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 16 ) THEN
      DO I=1,MAXROW
   16   continue
	READ( INUNIT, '(A)', ERR=16, END=200 ) RECORD
	READ( RECORD, *, ERR=16, END=16 ) 
     #	C1(I), C2(I), C3(I), C4(I), C5(I),
     #	C6(I), C7(I), C8(I), C9(I), C10(I),
     #	C11(I), C12(I), C13(I), C14(I), C15(I),
     #	C16(I)
      END DO

* 17 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 17 ) THEN
      DO I=1,MAXROW
   17   continue
	READ( INUNIT, '(A)', ERR=17, END=200 ) RECORD
	READ( RECORD, *, ERR=17, END=17 ) 
     #	C1(I), C2(I), C3(I), C4(I), C5(I),
     #	C6(I), C7(I), C8(I), C9(I), C10(I),
     #	C11(I), C12(I), C13(I), C14(I), C15(I),
     #	C16(I), C17(I)
      END DO
* 18 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 18 ) THEN
      DO I=1,MAXROW
   18   continue
	READ( INUNIT, '(A)', ERR=18, END=200 ) RECORD
	READ( RECORD, *, ERR=18, END=18 ) 
     #	C1(I), C2(I), C3(I), C4(I), C5(I),
     #	C6(I), C7(I), C8(I), C9(I), C10(I),
     #	C11(I), C12(I), C13(I), C14(I), C15(I),
     #	C16(I), C17(I), C18(I)
      END DO

* 19 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 19 ) THEN
      DO I=1,MAXROW
   19   continue
	READ( INUNIT, '(A)', ERR=19, END=200 ) RECORD
	READ( RECORD, *, ERR=19, END=19 ) 
     #	C1(I), C2(I), C3(I), C4(I), C5(I),
     #	C6(I), C7(I), C8(I), C9(I), C10(I),
     #	C11(I), C12(I), C13(I), C14(I), C15(I),
     #	C16(I), C17(I), C18(I), C19(I)
      END DO

* 20 COLUMNS ----------------------------------------------
      ELSE IF( NC .EQ. 20 ) THEN
      DO I=1,MAXROW
   20   continue
	READ( INUNIT, '(A)', ERR=20, END=200 ) RECORD
	READ( RECORD, *, ERR=20, END=20 ) 
     #	C1(I), C2(I), C3(I), C4(I), C5(I),
     #	C6(I), C7(I), C8(I), C9(I), C10(I),
     #	C11(I), C12(I), C13(I), C14(I), C15(I),
     #	C16(I), C17(I), C18(I), C19(I), C20(I)
      END DO

      ELSE IF( NC .LE. MAXCOL ) THEN
      DO I=1,MAXROW
   50   continue
	READ( INUNIT, '(A)', ERR=50, END=200 ) RECORD
	READ( RECORD, *, ERR=50, END=50 ) 
     #		(ROW(J),J=1,NC )
	IF( NC .GE.1 ) C1(I) = ROW(1)
	IF( NC .GE.2 ) C2(I) = ROW(2)
	IF( NC .GE.3 ) C3(I) = ROW(3)
	IF( NC .GE.4 ) C4(I) = ROW(4)
	IF( NC .GE.5 ) C5(I) = ROW(5)
	IF( NC .GE.6 ) C6(I) = ROW(6)
	IF( NC .GE.7 ) C7(I) = ROW(7)
	IF( NC .GE.8 ) C8(I) = ROW(8)
	IF( NC .GE.9 ) C9(I) = ROW(9)
	IF( NC .GE.10 ) C10(I) = ROW(10)
	IF( NC .GE.11 ) C11(I) = ROW(11)
	IF( NC .GE.12 ) C12(I) = ROW(12)
	IF( NC .GE.13 ) C13(I) = ROW(13)
	IF( NC .GE.14 ) C14(I) = ROW(14)
	IF( NC .GE.15 ) C15(I) = ROW(15)
	IF( NC .GE.16 ) C16(I) = ROW(16)
	IF( NC .GE.17 ) C17(I) = ROW(17)
	IF( NC .GE.18 ) C18(I) = ROW(18)
	IF( NC .GE.19 ) C19(I) = ROW(19)
	IF( NC .GE.20 ) C20(I) = ROW(20)
	IF( NC .GE.21 ) C21(I) = ROW(21)
	IF( NC .GE.22 ) C22(I) = ROW(22)
	IF( NC .GE.23 ) C23(I) = ROW(23)
	IF( NC .GE.24 ) C24(I) = ROW(24)
	IF( NC .GE.25 ) C25(I) = ROW(25)
	IF( NC .GE.26 ) C26(I) = ROW(26)
	IF( NC .GE.27 ) C27(I) = ROW(27)
	IF( NC .GE.28 ) C28(I) = ROW(28)
	IF( NC .GE.29 ) C29(I) = ROW(29)
	IF( NC .GE.30 ) C30(I) = ROW(30)
	IF( NC .GE.31 ) C31(I) = ROW(31)
	IF( NC .GE.32 ) C32(I) = ROW(32)
	IF( NC .GE.33 ) C33(I) = ROW(33)
	IF( NC .GE.34 ) C34(I) = ROW(34)
	IF( NC .GE.35 ) C35(I) = ROW(35)
	IF( NC .GE.36 ) C36(I) = ROW(36)
	IF( NC .GE.37 ) C37(I) = ROW(37)
	IF( NC .GE.38 ) C38(I) = ROW(38)
	IF( NC .GE.39 ) C39(I) = ROW(39)
	IF( NC .GE.40 ) C40(I) = ROW(40)
	IF( NC .GE.41 ) C41(I) = ROW(41)
	IF( NC .GE.42 ) C42(I) = ROW(42)
	IF( NC .GE.43 ) C43(I) = ROW(43)
	IF( NC .GE.44 ) C44(I) = ROW(44)
	IF( NC .GE.45 ) C45(I) = ROW(45)
	IF( NC .GE.46 ) C46(I) = ROW(46)
	IF( NC .GE.47 ) C47(I) = ROW(47)
	IF( NC .GE.48 ) C48(I) = ROW(48)
	IF( NC .GE.49 ) C49(I) = ROW(49)
	IF( NC .GE.50 ) C50(I) = ROW(50)
	IF( NC .GE.51 ) C51(I) = ROW(51)
	IF( NC .GE.52 ) C52(I) = ROW(52)
	IF( NC .GE.53 ) C53(I) = ROW(53)
	IF( NC .GE.54 ) C54(I) = ROW(54)
	IF( NC .GE.55 ) C55(I) = ROW(55)
	IF( NC .GE.56 ) C56(I) = ROW(56)
	IF( NC .GE.57 ) C57(I) = ROW(57)
	IF( NC .GE.58 ) C58(I) = ROW(58)
	IF( NC .GE.59 ) C59(I) = ROW(59)
	IF( NC .GE.60 ) C60(I) = ROW(60)
	IF( NC .GE.61 ) C61(I) = ROW(61)
	IF( NC .GE.62 ) C62(I) = ROW(62)
	IF( NC .GE.63 ) C63(I) = ROW(63)
	IF( NC .GE.64 ) C64(I) = ROW(64)
	IF( NC .GE.65 ) C65(I) = ROW(65)
	IF( NC .GE.66 ) C66(I) = ROW(66)
	IF( NC .GE.67 ) C67(I) = ROW(67)
	IF( NC .GE.68 ) C68(I) = ROW(68)
	IF( NC .GE.69 ) C69(I) = ROW(69)
	IF( NC .GE.70 ) C70(I) = ROW(70)
	IF( NC .GE.71 ) C71(I) = ROW(71)
	IF( NC .GE.72 ) C72(I) = ROW(72)
	IF( NC .GE.73 ) C73(I) = ROW(73)
	IF( NC .GE.74 ) C74(I) = ROW(74)
	IF( NC .GE.75 ) C75(I) = ROW(75)
	IF( NC .GE.76 ) C76(I) = ROW(76)
	IF( NC .GE.77 ) C77(I) = ROW(77)
	IF( NC .GE.78 ) C78(I) = ROW(78)
	IF( NC .GE.79 ) C79(I) = ROW(79)
	IF( NC .GE.80 ) C80(I) = ROW(80)
	IF( NC .GE.81 ) C81(I) = ROW(81)
	IF( NC .GE.82 ) C82(I) = ROW(82)
	IF( NC .GE.83 ) C83(I) = ROW(83)
	IF( NC .GE.84 ) C84(I) = ROW(84)
	IF( NC .GE.85 ) C85(I) = ROW(85)
	IF( NC .GE.86 ) C86(I) = ROW(86)
	IF( NC .GE.87 ) C87(I) = ROW(87)
	IF( NC .GE.88 ) C88(I) = ROW(88)
	IF( NC .GE.89 ) C89(I) = ROW(89)
	IF( NC .GE.90 ) C90(I) = ROW(90)
	IF( NC .GE.91 ) C91(I) = ROW(91)
	IF( NC .GE.92 ) C92(I) = ROW(92)
	IF( NC .GE.93 ) C93(I) = ROW(93)
	IF( NC .GE.94 ) C94(I) = ROW(94)
	IF( NC .GE.95 ) C95(I) = ROW(95)
	IF( NC .GE.96 ) C96(I) = ROW(96)
	IF( NC .GE.97 ) C97(I) = ROW(97)
	IF( NC .GE.98 ) C98(I) = ROW(98)
	IF( NC .GE.99 ) C99(I) = ROW(99)
	IF( NC .GE.100 ) C100(I) = ROW(100)
	IF( NC .GE.101 ) C101(I) = ROW(101)
	IF( NC .GE.102 ) C102(I) = ROW(102)
	IF( NC .GE.103 ) C103(I) = ROW(103)
	IF( NC .GE.104 ) C104(I) = ROW(104)
	IF( NC .GE.105 ) C105(I) = ROW(105)
	IF( NC .GE.106 ) C106(I) = ROW(106)
	IF( NC .GE.107 ) C107(I) = ROW(107)
	IF( NC .GE.108 ) C108(I) = ROW(108)
	IF( NC .GE.109 ) C109(I) = ROW(109)
	IF( NC .GE.110 ) C110(I) = ROW(110)
	IF( NC .GE.111 ) C111(I) = ROW(111)
	IF( NC .GE.112 ) C112(I) = ROW(112)
	IF( NC .GE.113 ) C113(I) = ROW(113)
	IF( NC .GE.114 ) C114(I) = ROW(114)
	IF( NC .GE.115 ) C115(I) = ROW(115)
	IF( NC .GE.116 ) C116(I) = ROW(116)
	IF( NC .GE.117 ) C117(I) = ROW(117)
	IF( NC .GE.118 ) C118(I) = ROW(118)
	IF( NC .GE.119 ) C119(I) = ROW(119)
	IF( NC .GE.120 ) C120(I) = ROW(120)
	IF( NC .GE.121 ) C121(I) = ROW(121)
	IF( NC .GE.122 ) C122(I) = ROW(122)
	IF( NC .GE.123 ) C123(I) = ROW(123)
	IF( NC .GE.124 ) C124(I) = ROW(124)
	IF( NC .GE.125 ) C125(I) = ROW(125)
	IF( NC .GE.126 ) C126(I) = ROW(126)
	IF( NC .GE.127 ) C127(I) = ROW(127)
	IF( NC .GE.128 ) C128(I) = ROW(128)
	IF( NC .GE.129 ) C129(I) = ROW(129)
	IF( NC .GE.130 ) C130(I) = ROW(130)
	IF( NC .GE.131 ) C131(I) = ROW(131)
	IF( NC .GE.132 ) C132(I) = ROW(132)
	IF( NC .GE.133 ) C133(I) = ROW(133)
	IF( NC .GE.134 ) C134(I) = ROW(134)
	IF( NC .GE.135 ) C135(I) = ROW(135)
	IF( NC .GE.136 ) C136(I) = ROW(136)
	IF( NC .GE.137 ) C137(I) = ROW(137)
	IF( NC .GE.138 ) C138(I) = ROW(138)
	IF( NC .GE.139 ) C139(I) = ROW(139)
	IF( NC .GE.140 ) C140(I) = ROW(140)
      END DO

      END IF
	I = MAXROW + 1

* final read to pick up the eof

	READ( INUNIT, '(A)', ERR=200, END=200 ) RECORD

* give warning if there was data past the final row

	WRITE(*,*)  '** WARNING: Data beyond row', MAXROW, ' was not read.'
  200	NROW = I - 1

* close input file

	CLOSE( UNIT = INUNIT )

* normal return

 1000	RETURN

* Error return

   99	WRITE(*,*) CHAR(7),'** Error opening file. '
	WRITE(*,*) '** ', FILE(:LEN1(FILE))
	GOTO 999

  999   WRITE(*,*) '** LOADLC aborted.'
	NROW = 0
	RETURN
	END
************************      
      subroutine cpu( tcpu )
*
* cpu time in seconds since previous call
*
* 1994 Mar kdh @ Utrecht - use dtime() to get the answer
* 2010 Mar kdh @ LCO - try cpu_time from gfortran
	real*4 told/-1./
c	real*4 tarray(2)
c	external dtime

	call cpu_time( tnow )
	tcpu = tnow - told
	if( told .lt. 0. ) tcpu = 0.
	told = tnow

c	tcpu = dtime( tarray )
c	tcpu = etime( tarray )
c	tuser = tarray(1)
c	tsystem = tarray(2)
	return
	end
C* PGSTART ... open a plot
C+
      SUBROUTINE PGSTART( NX, NY, NEW )
*
* PGSTART opens and PGSTOP closes a plot.
* Like PGBEGIN and PGEND, but keeping track of previous plot device
* and local plot queues so plots can be spooled at run time.
*
* Input:
*	NX,NY	I4 number of plots per page in X and Y direction
*		selected interactively if input values are 0
* Output:
*	NEW	I4 1 if opening a new plot, otherwise
*	NX,NY	I4 if changed by interaction
*
C--
* 1980 KDH at CIT - Original version
* 1985 Nov KDH @ STScI - added QMS devices
* 1986 Apr KDH @ STScI - added laser printer
* 1987 Jul KDH @ STScI - added Steve Baron's DeAnza driver
* 1987 Nov AWC @ STScI - option not to plot named Q,Z,P,V files.
* 1987 Dec KDH @ STScI - add retrographics option
* 1987 Dec KDH @ STScI - add metafile options
* 1988 Feb KDH @ STScI - call PGIDEN in PGSTOP
* 1988 Nov KDH @ STSCI - PRINT/NOFLAG
* 1989 Apr KDH @ STSCI - CHANGE QUEUE NAMES
* 1989 Jun KDH @ STSCI - put flag page back in
* 1989 Oct KDH @ STSCI - change queue 3C_T800 to 3C_T2400
* 1989 Nov KDH @ STSCI - add VAX WS option, also ?
* 1989 Dec KDH @ STSCI - update STScI printer queues
* 1990 Jan KDH @ STScI - interactive setting of NX, NY if 0 on input
* 1990 Feb KDH @ STScI - pgbegin option, add 3CLPS_POST queue
* 1991 Apr KDH @ STScI - Ant queue
* 1992 Mar Carole Hawell @ STScI - LP5 queue added
* 1992 Mar KDH @ STScI - distinguish VT125 and VT340 drivers
* 1992 Jul KDH @ Amsterdam - add amsterdam postscript queues
* 1992 Jul KDH @ cavad - add cavad canon queues
* 1992 Jul KDH @ cavad - add pericom_mg
* 1992 Oct KDH @ STScI - update STScI postscript queues
* 1993 Feb CAH @ STScI - update P&A postscript queues
* 1993 Apr KDH @ STScI - change SOL331 and ANT423S to 331PS1 and S423PS
* 1994 Jan KDH @ Utrecht - write(*,'(A)') '$...' -> write(*,'(A,$)') ' ...'
* 1994 Jan KDH @ Utrecht - add lpr for use on unix versions
* 1994 Jan KDH @ Utrecht - add XW devicetype
* 1994 Mar KDH @ Utrecht - eliminate de-anzas
* 1994 Mar KDH @ Utrecht - change LIB$SPAWN to SPAWN
* 1994 Mar KDH @ Utrecht - don't upcase filenames
* 1994 May KDH @ STScI - import Utrecht version, set OPSYS->VMS, SITE->STSCI
* 1995 Jan KDH @ St.Andrews - implement St-Andrews queues
* 1999 Sep KDH @ St.And - remove PGIDEN from PGSTOP
* 2000 Nov KDH @ Austin - implement Austin queues
* 2001 Jul KDH @ St.And - new St-Andrews queues
* 2001 Jul KDH @ StA - linux port: must use g77 -fno-automatic
* 2001 Jul KDH @ StA - simplify, strip vms & multi-sites
* 2001 Aug KDH @ StA - FULLPATH plot files, ASPECT
* 2004 Jan KDH @ StA - add err= to read (linux)
* 2005 Feb KDH @ StA - strip @ from plotfile
* 2006 May KDH @ StA - add farside printer
* 2010 Jun KDH @ StA - update StA printers
* 2014 Mar KDH @ StA - kdh -> kdh1
* 2014 Sep KDH @ StA - kdh1 -> kdh

        character*20 site
	character*1 aspect, queue, mode, oldmode
	character*256 reply
	character*64 terminal, devtype
	character*128 prtcom, plotdir, plotfile 
	character*128 psfile, cpsfile, giffile, pltfile, jfile
	character*256 task, path
	logical soft, spool, auto
	integer pgbegin
	common /pgstartmode/ mode
	data icall/0/
	icall = icall + 1

* first call initializations
      if( icall .eq. 1 ) then
	aspect = 'L'
	mode = 'H'
	queue = '?'
	iopen = 0
* page format
	nx0 = 1
	ny0 = 1

* site
	site = 'STANDREWS'
	write(*,*) 'PGSTART at ', SITE( :LEN1(SITE) )

* default directory
	plotdir='$SCR/'
* if system doesn't translate the environment variable $SCR
* correctly try the current directory as the default
c	plotdir = './'
	l1 = len1( plotdir )
        write(*,*) 'Default plot dir: ', plotdir(:l1)
	call fullpath( plotdir, path )
      if( plotdir(:l1) .ne. path(:len1(path)) ) then
	write(*,*) 'Default path: ', path(:len1(path))
      end if

* default plot files
	psfile = plotdir(:l1) // 'pgplot.ps'
	cpsfile = plotdir(:l1) // 'pgplot.cps'
	giffile = plotdir(:l1) // 'pgplot.gif'
	pltfile = plotdir(:l1) // 'pgplot.plt'

* end of first call initializations
      end if

* every call initializations
	oldmode = mode
	new = 0
	goto 1

* menu of device types
   10	write(*,*) 'PGSTART device types:'
	write(*,*) ' X=Xserve W=xWin T=Tektronix G=Gif H=Help'
	write(*,*) ' P=PostScript C=ColourPS N=null B=pgBegin'

* select device type
    1	write(*,'(A,$)')
     &	' PGSTART device type (H=help) [' // mode // '] '
	read(*,'(A)', err=1 ) reply
	call upcase( reply, reply )
	if( reply .eq. 'H' ) goto 10
	if( reply .eq. ' ' ) reply = mode

* soft devices (terminals) 
      if(    reply .eq. 'T'
     &	.or. reply .eq. 'X'
     &	.or. reply .eq. 'W'
     &	.or. reply .eq. 'B'
     &	.or. reply .eq. 'N'
     &	) then
	mode = reply
	soft = .true.
	spool = .false.
	auto = .false.

* hardcopy devices require a disk file
      else if( reply .eq. 'P'
     &	.or.   reply .eq. 'C'
     &	.or.   reply .eq. 'G' ) then
	mode = reply
	soft = .false.
	spool = mode .eq. 'P' .or. mode .eq. 'C'
	auto = .false.

* trap invalid replies
      else if( reply .ne. ' ' ) then
	goto 10
      else if( mode .eq. ' ' .or. mode .eq. 'H' ) then
	goto 10
      end if

* close previous plot if one was open
      if( mode .ne. oldmode .and. iopen .ne. 0 ) then
	call pgend
	iopen = 0
      end if

* terminals
      if( soft ) then
	if( mode .eq. 'N' ) terminal = '/NU'
	if( mode .eq. 'T' ) terminal = '/TE'
	if( mode .eq. 'W' ) terminal = '/XW'
	if( mode .eq. 'X' ) terminal = '/XSERV'
	if( mode .eq. 'B' ) terminal = '?'

* plot file 
      else

	jfile = pltfile
	if( mode .eq. 'P' ) jfile = psfile
	if( mode .eq. 'C' ) jfile = cpsfile
	if( mode .eq. 'G' ) jfile = giffile
      if( mode .ne. oldmode ) then
	plotfile = pltfile
	if( mode .eq. 'P' ) plotfile = psfile
	if( mode .eq. 'C' ) plotfile = cpsfile
	if( mode .eq. 'G' ) plotfile = giffile
      end if

	if( plotfile .ne. ' ' )
     &	write(*,*) '<CR> uses plotfile ', plotfile( :len1(plotfile) )
	write(*,*) 'JUNK autoprints ', jfile( :len1(jfile) )
	if( mode .eq. 'P' ) write(*,'(A,$)') ' PostScript file : '
	if( mode .eq. 'C' ) write(*,'(A,$)') ' ColourPostScript file : '
	if( mode .eq. 'G' ) write(*,'(A,$)') ' GIF file : '
 	read(*,'(A)') reply
      if( reply .eq. ' ' ) then
	reply = plotfile
      else
	if( reply(1:1) .eq. '@' ) reply = reply(2:)
	plotfile = reply
	call upcase( reply, reply )
	auto = .false.
      if( reply .eq. 'JUNK') then
	plotfile = jfile
	auto = .true.
      end if
	if( plotfile .eq. ' ' ) goto 10
      end if

      if( spool ) then

* print queues at each site
      if( site .eq. 'STANDREWS' ) then
   16	write(*,'(A,$)')
     &	' Which Queue (?,U,K,D,T,C) ? [' // queue // '] '
	read(*,'(A)') reply
	if( reply .eq. ' ' ) reply = queue
	call upcase( reply, reply )
	prtcom = ' '
      if( reply .eq. 'K') then
	prtcom = 'lpr -Pkdh1'
c	prtcom = 'lpr -Pkdh'
      else if( reply .eq. 'D') then
	prtcom = 'lpr -Pduplex'
      else if( reply .eq. 'T') then
	prtcom = 'lpr -PT644N'
      else if( reply .eq. 'C') then
	prtcom = 'lpr -PHPCP2025'
      else if( reply .eq. 'U') then
	prtcom = 'lpr -PUniPrint_Ricoh'
      else
        write(*,*) 'St.Andrews Printer queues :'
        write(*,*) ' K  kdh1    (Keith Horne''s office)'
        write(*,*) ' D  duplex  double-sided (Starlink room)'
        write(*,*) ' T  T644N   double-sided (south corridor)'
        write(*,*) ' C HPCP2025 colour (Starlink room)'
        write(*,*) ' U UniPrint_Ricoh'
	goto 16
      end if
	queue = reply

* last site
      else
	write(*,*) '** PGSTART: no queues at ', site( :len1(site) )
	spool = .false.
	auto = .false.
      end if

      end if

* portrait vs landscape
	write(*,'(A,$)')
     &	' P(ortrait) or L(andscape) ? [' // aspect // '] '
	read(*,'(A)') reply
	call upcase( reply, reply )
	if( reply .eq. ' ' ) reply = aspect
	if( reply .eq. 'P' .or. reply .eq. 'L' ) aspect = reply

* device type
      if( aspect .eq. 'L' ) then
	if( mode .eq. 'P' ) devtype  =  'PS'
	if( mode .eq. 'C' ) devtype  =  'CPS'
	if( mode .eq. 'G' ) devtype  =  'GIF'
      else if( aspect .eq. 'P' ) then
	if( mode .eq. 'P' ) devtype  =  'VPS'
	if( mode .eq. 'C' ) devtype  =  'VCPS'
	if( mode .eq. 'G' ) devtype  =  'VGIF'
      end if

* end of hard devices
      end if

* need to open a new plot ?
      if( iopen .le. 0 ) then

* page format
      if( nx .le. 0 .or. ny .le. 0 ) then
   40	write(*,'(A,I3,A,I3,A,$)')
     #	' Plot page format (NX,NY) [', nx0, ' ,', ny0, ' ] '
	read(*,'(A)') reply
      if( reply .eq. ' ' ) then
	nx = nx0
	ny = ny0
      else
	read( reply,*,err=40, end=40 ) nx, ny
	IF( nx .le. 0 .or. ny .le. 0 ) goto 40
      end if
      end if
	nx0 = nx
	ny0 = ny

* open the plot

* terminals
      if( soft ) then
	write(*,*) 'PGBEGIN( 0, ',
     &	'''', terminal( :len1( terminal ) ), ''',',
     &	nx, ',', ny, ' )'
	ierr = pgbegin( 0, terminal, nx, ny )

* hardcopy devices
      else
* get full path
	call fullpath( plotfile, path )
* call pgbegin
	task = '"' // path( :len1( path ) ) // '"'
     &		// '/' // devtype( :len1( devtype ) )
	write(*,*) 'PGBEGIN( 0, ', task( :len1(task) ),
     &	', ', nx, ', ', ny, ' )'
	ierr = pgbegin( 0, task( :len1(task) ), nx, ny )
      end if
* failure
      if( ierr .ne. 1 ) then
	write(*,*) '** FAILED TO OPEN PLOT.'
	goto 1
      end if
* success
	iopen = 1
	new = 1
      end if

c	write(*,*) 'PGSTART spool', spool, ' auto', auto
	return

C* PGSTOP ... close plot, send output to printer queue
C+
	entry pgstop
C--
* write time and username at bottom of plot
c	call pgiden

* close the plot
	call pgend
	iopen = 0

* how to print the plot file
	write(*,*) 'PGSTOP spool', spool, ' auto', auto
      if( spool ) then
	call word1( prtcom, l1, l2 )
	l2 = len1( prtcom )
	write(*,*) 'COMMAND: ', prtcom( l1:l2 )
	call word1( plotfile, l3, l4 )
	write(*,*) 'PLOTFILE: ', plotfile( l3:l4 )
	call fullpath( plotfile, path )
	call word1( path, l5, l6 )
      if( plotfile( l3:l4 ) .ne. path( l5:l6 ) ) then
	write(*,*) 'PATH: ', path( l5:l6 )
      end if
	task = prtcom( l1:l2 ) // ' ' // path( l5:l6 )
	write(*,*) 'TASK: ', task( :len1( task ) )

* decide to spool, or not to spool
      if( .not. auto ) then
	write(*,'(A,$)') ' Queue file for plotting? [Y] '
	read(*,'(A)') reply
	call upcase( reply, reply )
	if( reply .eq. ' ' ) reply = 'Y'
      else
	reply = 'Y'
      end if

* print the plot file
      if( reply .eq. 'Y' ) then
	call spawn( task( : len1( task ) ) )
      end if
      end if

	return
	end
*=====================================
	logical function pgagain()
* plot again ?
* 1997 Aug Keith Horne @ St.Andrews
* 1999 Oct KDH @ St.And - make logical
	character*10 reply
	write(*,'(A,$)') ' Plot again ? [N] '
	read(*,'(A)') reply
	call upcase( reply, reply )
	pgagain = .false.
	if( reply.eq.'Y') pgagain = .true.
	return
	end
*----------------------------------------
	subroutine pgbgw
* set background/foreground colours to white/black
* 2010 Oct Keith Horne @ St-And
	call pgscr( 0, 1., 1., 1. )
	call pgscr( 1, 0., 0., 0. )
	return
	end
