	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 - add end=10
*
	REAL*8 DEF, VAL, WANT
	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 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 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
	CHARACTER*(*) STRING
	CHARACTER*1 TAB, BLANK, TEST
	TAB = CHAR(9)
	BLANK = ' '
	LMAX = LEN1(STRING)
      DO I1=1,LMAX
	TEST = STRING(I1:I1)
      IF( TEST.NE.BLANK .AND. TEST.NE.TAB ) THEN
	L1 = I1
      DO I2=L1+1,LMAX
	TEST = STRING(I2:I2)
      IF( TEST.EQ.BLANK .OR. TEST.EQ.TAB ) THEN
	L2 = I2-1
	RETURN
      END IF
      END DO
	L2 = LMAX
	RETURN
      END IF
      END DO
	L1 = 0
	L2 = 0
	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 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
        character*(*) title
        real*8 hjd
        hjd = dat
        call append_hjdsig( title, hjd, sig, kfig )
        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
        character*(*) title
        real*8 hjd
        character*50 fmt
        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 )
        sigdiv = sig / div
        call append_d( title, hjd, fmt )
        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, ')' )
      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 )
        call append( title, '(' )
        call append_f( title, sig/ div, fmt )
        call append( title, ')' )
        call append( title, 'E' )
        call append_i( title, lpow )
      else
        idiv = 10 ** kpow
        call append_i( title, nint( hjd / idiv ) * idiv )
        call append( title, '(' )
        call append_i( title, nint( sig / idiv ) * idiv )
        call append( title, ')' )
      end if
        call nospace( title(l1+1:) )
        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

	real*4 function ran( 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
	logical firstcall/.true./
	save firstcall
* system-supplied
c	ran = rand( iseed )
* numerical recipes routine 
      if( firstcall ) then
	is = - iabs( iseed )
c	ran = ran1( is )
	ran = ran3( is )
	firstcall = .false.
      else
	is = iseed
c	ran = ran1( is )
	ran = ran3( is )
      end if
	iseed = is
	return
	end
*==================================
	FUNCTION RAN1(IDUM)
* random number uniform on (0-1).
* Numerical Recipes routine.
*
* NOTE: RAN1 SEEMS TO HAVE PROBLEMS FOR LARGE NUMBERS OF CALLS.
* THERE IS A SYSTEMATIC BIAS REVEALED BY TESTS OF CENTRAL LIMIT THEOREM
* RAN3 DOES A BETTER JOB.
* 2008 KEITH HORNE @ ST-ANDREWS
* 2010 Jun KDH @ St.A - pause -> call pause
	DIMENSION R(97)
	PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=1./M1)
	PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=1./M2)
	PARAMETER (M3=243000,IA3=4561,IC3=51349)
	DATA IFF /0/
	IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
		IFF=1
		IX1=MOD(IC1-IDUM,M1)
		IX1=MOD(IA1*IX1+IC1,M1)
		IX2=MOD(IX1,M2)
		IX1=MOD(IA1*IX1+IC1,M1)
		IX3=MOD(IX1,M3)
		DO J=1,97
			IX1=MOD(IA1*IX1+IC1,M1)
			IX2=MOD(IA2*IX2+IC2,M2)
			R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
		ENDDO
		IDUM=1
	ENDIF
	IX1=MOD(IA1*IX1+IC1,M1)
	IX2=MOD(IA2*IX2+IC2,M2)
	IX3=MOD(IA3*IX3+IC3,M3)
	J=1+(97*IX3)/M3
	IF(J.GT.97.OR.J.LT.1) call PAUSE( 'RAN1 FAILED. :(' )
	RAN1=R(J)
	R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
	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
************************      
      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
************************
	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
	character*(*) task
	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 pause( message )
* substitute for pause statement (withdrawn from gfortran)
* 2010 Jun Keith Horne @ St.Andrews
	character*(*) message
	character*50 reply
	l = len1( message )
	if( l .gt. 0 ) write(*,*) '** pause ', message(:l)
   10	call inq( '** pause C(ontinue) or D(ie) ?', 'C', reply )
	call upcase( reply, reply )
	if( reply .eq. 'D' ) stop 'bye for now.'
	if( reply .ne. 'C' ) goto 10
	return
	end
*---------------------------------------------
      subroutine  get_dmy( id, im, iy )

* 2006 Feb Keith Horne @ St.Andrews -- from ian bonnell
* 2010 Jun KDH @ StA - idate -> date_and_time 
* 2010 Jun KDH @ StA - itime -> time() and gmtime

      INTEGER*4 iarray(3), istatb(13)
      INTEGER*4 ival1, ival2
      CHARACTER*11 job, inname
      CHARACTER*7 file
	character*20 date, localtime
	integer*4 time, stime, tarray(9)
c
c--Get date
c
* idate fails with ifort compiler
c      CALL idate(iarray)
c      id = iarray(1)
c      im = iarray(2)
c      iy = iarray(3)

	call date_and_time( date, localtime )
	read( date, '(i4,i2,i2)') iy, im, id
	read( localtime, '(i2,i2,i2,1x,i3)') ihr, imin, isec, ims

      RETURN
c
c--Get time
c
      ENTRY get_hmst( ih, im, is, fhour )

* itime returns local time, not UT
c      CALL itime(iarray)
c      ih = iarray(1)
c      im = iarray(2)
c      is = iarray(3)

	stime = time()
	call gmtime( stime, tarray )
	ih = tarray(3)
	im = tarray(2)
	is = tarray(1)

	fhour = ih + ( im + is / 60.d0 ) / 60.d0
      RETURN
c
c--Get time used since begining
c
c      ENTRY getused(tused)
c      tused = etime(tarray)
c      RETURN
c
c--Check for file status
c
      ENTRY statfile( file, ifsize )
c      ifile = lstat( file, istatb )
      ifsize = istatb(8)
      RETURN
c
c--Get argument on command line
c
      ENTRY getcom(job, inname)
      ival1 = 1
      ival2 = 2
      CALL getarg( ival1, job )
      CALL getarg( ival2, inname )
      WRITE (*,*) 'Job: ', job, inname
      RETURN

      END
      function jdate( iy, im, id )
C  Julian Date at 12 hrs UT
* pre-1983 code from Peter Young.
* 1999 Oct Keith Horne @ St.Andrews - eliminate type
* 2002 Mar KDH @ St.And - mdy -> ymd
      if(iy.lt.1801.or.iy.gt.2099) then
	write(*,*) '** JDATE: Year out of range 1801 - 2099', iy
	stop
      end if
      if(im.lt.1.or.im.gt.12) then
	write(*,*) '** JDATE: Month out of range', im
	stop
      end if
      if(id.lt.1.or.id.gt.31) then
	write(*,*) '** JDATE: Day out of range', id
	stop
      end if
	jj = 367*iy - (7*(iy+(im+9)/12))/4 + (275*im)/9 + id + 1721014
	if( 100*iy + im .lt. 190002.5 ) jj = jj + 1
	jdate = jj
	return
	end

      SUBROUTINE CALENDAR( DAYJULIAN, IYR, IMO, IDY, HOURS )
*
*  Find calendar date and time from Julian date.
*
*  Input:
*	DAYJULIAN	= (REAL*8) JULIAN DATE
*  Output:
*	IYR,IMO,IDY	= (INTEGER*4) CALENDAR DATE
*	HOURS		= (REAL*8) TIME (HOURS)
*
* 1995 Aug KDH @ St-Andrews -- declare HOURS real*8
* 2000 Jan KDH @ St-Andrews -- make millennium leap years
* 2002 Mar KDH @ St-And -- mdy -> ymd
	REAL*8 DAYJULIAN, DAY, HOURS
	INTEGER*4 NDAY(12)
	DATA NDAY/0,31,59,90,120,151,181,212,243,273,304,334/
*  Find the year
	IYR = DAYJULIAN / 365.2475D0 - 4712.D0
   10	JAN1 = JDATE( IYR, 1, 1 )
	DAY = DAYJULIAN - JAN1 + 0.5D0
*	write(*,*) DAYJULIAN, JAN1, IYR, DAY
      IF( DAY.LT.0.D0 ) THEN
	IYR=IYR-1
	GOTO 10
      END IF
	LEAP=0
	IF( (IYR/4)*4 .EQ. IYR ) LEAP=1
	IF( (IYR/100)*100 .EQ. IYR ) LEAP=0
	IF( (IYR/1000)*1000 .EQ. IYR ) LEAP=1
      IF( DAY.GT. 365+LEAP ) THEN
	IYR=IYR+1
	GOTO 10
      END IF
*  Find the month
	IMO=1
	IF( DAY .LT. NDAY(2) ) GOTO 30
      DO IMO = 2, 11
	IF( DAY .LT. NDAY(IMO+1) + LEAP ) GOTO 30
      END DO
	IMO=12
*  Find the day
   30	IF(IMO.LT.3) LEAP=0
	NDY = NDAY(IMO) + LEAP
	IDY = DAY - NDY + 1
*  Find the time
	HOURS = (DAY - NDY - IDY + 1) * 24.D0
	RETURN
	END
      SUBROUTINE TSPLIT( TIME, IH, IM, IS, MS )
*
*  Decomposes time into hours, minutes, seconds and milliseconds.
*
* Input:
*	TIME	= (R8) time in hours
* Output:
*	IH,IM,IS,MS	= (I4) hours, minutes, seconds and milliseconds
*
* Jan 1987 KDH @ STScI - correct 0.5 ms roundoff
* Oct 1987 KDH @ STScI - use mod 24 hours.
* Jul 1990 KDH @ STScI - take out mod(24), routine is used for ra and dec.
*
	REAL*8 TIME, DUMMY
c	DUMMY = MOD( TIME, 24.D0 )
	DUMMY = TIME
     #	 + 0.5D0/3600.D3	! rounds time to nearest ms
	IH = DUMMY
	DUMMY = (DUMMY-IH)*60.D0
	IM = DUMMY
	DUMMY = (DUMMY-IM)*60.D0
	IS = DUMMY
	MS = (DUMMY-IS)*1000.D0
	RETURN
	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
	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 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 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
	real*4 x(*), xmn, xmx
	xmn = x(1)
	xmx = xmn
	npos = 0
      do i=1,n
      if( x(i) .gt. 0. ) then
	npos = npos + 1
	xmn = min( xmn, x(i) )
	xmx = max( xmx, x(i) )
      end if
      end do
	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

******************************************
	subroutine medmad( n, x, xmed, xmad )
* median and mean absolute deviation of n data points
* Input:
*	n	i4 number of data
*	x(n)	r4 data
* output:
*	xmed	r4 median value
*	xmad	r4 mean absolute deviation
* 2009 May Keith Horne @ St.Andrews
	real*4 x(*)
	real*8 sum
	xmed = xmedian( n, x )
	sum = 0.d0
      do i=1,n
	sum = sum + abs( x(i) - xmed )
      end do
	xmad = sum / max( 1, n )
	return
	end
	function rad2deg( rad )
* convert radians to degrees
	logical firstcall/.true./
      if( firstcall ) then
	f = 45. / atan2( 1., 1. )
	firstcall = .false.
      end if
	rad2deg = rad * f
	return
	end
* 
	function deg2rad( deg )
* convert degrees to radians
	logical firstcall/.true./
      if( firstcall ) then
	f = atan2( 1., 1. ) / 45.
	firstcall = .false.
      end if
	deg2rad = deg * f
	return
	end

******************************************
	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 avgrmspdf( n, p, x, avg, rms )
* mean and rms averaged over a pdf
* 2003 Jul Keith Horne @ St.Andrews
* 2008 Aug KDH @ StA - rms was wrong, fixed
	real*4 p(*), x(*)
	real*8 sum, sump
	avg = x(1)
	if( n .le. 0 ) return
	sump = 0.d0
	sum = 0.d0
      do i=1,n
	sump = sump + p(i)
	sum  =  sum + p(i) * x(i)
      end do
	if( sump .le. 0. ) return
	avg = sum / sump
	sum = 0.d0
      do i = 1, n
	dx = x(i) - avg
	sum = sum + p(i) * dx * dx
      end do
	var = sum / sump
	rms = sqrt( 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
	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
      if( sig(i) .gt. 0. ) then
	wgt = 1. / sig(i)**2
	sum = sum + wgt * dat(i)
	sumw = sumw + wgt
      end if
      end do
	if( sumw .le. 0.d0 ) return
	avg = sum / sumw
	var = 1. / sumw
	sigavg = sqrt( var )
	return
	end
********************************************
	subroutine optavgrms( n, dat, sig, bar, rms )
* optimal average
	real*4 dat(*), sig(*)
	real*8 sum
	bar = 0.
	rms = -1.
	if( n.le.0 ) return
	call optavg( n, dat, sig, bar, sigbar )
	nsum = 0
	sum = 0.d0
      do i=1,n
      if( sig(i) .gt. 0. ) then
	nsum = nsum + 1
	sum = sum + ( dat(i) - bar )**2
      end if
      end do
	var = sum / max( 1, nsum - 1 )
	rms = sqrt( var )
	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 fitline( n, t, dat, sig, t0, a, b, siga, sigb )
*  fit y(t)=a+b(t-t0) to data points dat(i)+/-sig(i) at times t(i)
*  input:
*	n	i4 number of data points
*	t(n)	r4 independent variable
*	dat(n)	r4 data values
*	sig(n)	r4 error bars
*  output:
*	t0,a,b	r4 best fit y(t) = a + b * ( t - t0 )
*	siga,sigb r4 fit variance is siga^2 + ((t-t0)*sigb)^2
*
* 2002 May Keith Horne @ St.Andrews
	real*4 t(*), dat(*), sig(*)
	real*8 w, wx, wy, wxx, wxy, det
* sanity check
	t0 = 0.
	a = 0.
	b = 0.
	siga = -1.
	sigb = -1.
	if( n .le. 0 ) return
* centroids and scales
	call optavg( n, t, sig, t0, tscale )
	call optavg( n, dat, sig, y0, yscale )
	call mnmx( n, t, tmn, tmx )
	call mnmx( n, dat, ymn, ymx )
	tscale = abs( tmx - tmn )
	yscale = abs( ymx - ymn )
	if( tscale. le. 0. ) tscale = 1.
	if( yscale .le. 0. ) yscale = 1.
* matrix elements
	w = 0.d0
	wx = 0.d0
	wy = 0.d0
	wxx = 0.d0
	wxy = 0.d0
      do i=1,n
      if( sig(i) .gt. 0. ) then
	wgt = ( yscale / sig(i) )**2
	x = ( t(i) - t0 ) / tscale
	y = ( dat(i) - y0 ) / yscale
	w  = w  + wgt
	wx  = wx  + wgt * x
	wy  = wy  + wgt * y
	wxx = wxx + wgt * x * x
	wxy = wxy + wgt * y * x
      end if
      end do
	if( w .le. 0.d0 ) return
* invert matrix
	det = w * wxx - wx * wx
      if( det .gt. 0.d0 ) then
	a = ( wy * wxx - wx * wxy ) / det
	b = (  w * wxy - wx * wy ) / det
	siga = 1.d0 / dsqrt( w )
c	sigb = 1.d0 / dsqrt( det )
	sigb = 1.d0 / dsqrt( wxx )
      else
	a = wy / w
	siga = 1.d0 / dsqrt( w )
	b = 0.
	sigb = 0.
      end if
* de-scale results
	a = a * yscale + y0
	siga = siga * yscale
	bscale = yscale / tscale
	b = b * bscale
	sigb = sigb * bscale
	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
	REAL*4 X(*)
	INTEGER*4 KEY(*)
	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(I) ) .GT. XTEST )
	KEY(J) = KEY(I)
	J = I
	I = I - INT
      END DO
	KEY(J) = ITEST
      END IF
      END DO
	RETURN
	END
********************************************
	SUBROUTINE QUICKSORT8( N, X, KEY )
*
* REAL*8 VERSION OF QUICKSORT
* QUICKSORT subroutine (Hoare)
* Produces index such that X(KEY(I)) ascends with I
*
* Input:
*	N	= number of items to be sorted
*	X(N)	= REAL*8 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
* 2009 Jul KDH @ St.And -- real*8 version
	REAL*8 X(*), XTEST
	INTEGER*4 KEY(*)
	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(I) ) .GT. XTEST )
	KEY(J) = KEY(I)
	J = I
	I = I - INT
      END DO
	KEY(J) = ITEST
      END IF
      END DO
	RETURN
	END
***************************************************
	FUNCTION XMEDIAN( NDATA, DATA )
* find median of an array of real numbers
*  Input:
*	NDATA	= I4 number of data values
*	DATA	= R4 array of data values
*  Output:
*	MEDVAL	= R4 median value of data
*
* Mar 1988 Keith Horne @ STScI - add special cases N=1,2, trap invalid inputs
* Jan 1995 KDH @ St.And - fix bug when NPTS = 2, eliminate type.
* Mar 2000 KDH @ St.And - MAXDATA 10000->20000
* Mar 2000 KDH @ St.And - medval(x,n) -> xmedian(n,x)
* Sep 2000 KDH @ St.And - shellsort -> quicksort
* Sep 2003 KDH @ St.And - MAXDATA 20,000->100,000

	REAL*4 DATA(*)
	PARAMETER ( MAXDATA = 100000 )
	INTEGER*4 IRANK( MAXDATA )

* trap invalid dimensions

      IF( NDATA.LT.1 ) THEN
	WRITE(*,*) '** INVALID DIMENSION', NDATA 
	GOTO 999
      ELSE IF( NDATA.GT.MAXDATA ) THEN
	WRITE(*,*) '** WARNING: buffer overflow in XMEDIAN', NDATA, MAXDATA
      END IF
	NPTS = MIN( NDATA, MAXDATA )

      IF( NPTS.EQ.1 ) THEN
	XMEDIAN = DATA(1)
	GOTO 1000
      ELSE IF( NPTS.EQ.2 ) THEN
	XMEDIAN = 0.5 * ( DATA(1) + DATA(2) )
	GOTO 1000
      END IF

* compute rank of each data value using a sort

	CALL QUICKSORT( NPTS, DATA, IRANK )
	MEDRANK = ( NPTS + 1 ) / 2 
	XMEDIAN = DATA( IRANK( MEDRANK ) )

* average two values if the number of points is even 

      IF( MEDRANK + MEDRANK .EQ.NPTS ) THEN
	XMEDIAN = 0.5 * ( XMEDIAN + DATA( IRANK( MEDRANK+1 ) ) )
      END IF

* normal return

 1000	RETURN

* error return

  999	XMEDIAN = 0.
	WRITE(*,*) '** XMEDIAN aborted.'
	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 SUBROUTINE MATINV.F 
C
C SOURCE								       
C   BEVINGTON, PAGES 302-303.
C
C Modified 6/26/89 by G. Kriss for use in specfit.
C	Changed array sizes to MAXPAR=500.
C Modified 5/31/95 by G. Kriss for use in specfit
C	Changed array sizes to MAXFREE=200
C
C PURPOSE
C   INVERT A SYMMETRIC MATRIX AND CALCULATE ITS DETERMINANT
C
C USAGE 
C   CALL MATINV (ARRAY, NORDER, DET)
C
C DESCRIPTION OF PARAMETERS
C   ARRAY  - INPUT MATRIX WHICH IS REPLACED BY ITS INVERSE
C   NORDER - DEGREE OF MATRIX (ORDER OF DETERMINANT)
C   DET    - DETERMINANT OF INPUT MATRIX
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
C   NONE
C
C COMMENT
C   DIMENSION STATEMENT VALID FOR NORDER UP TO 200
C
	SUBROUTINE MATINV (ARRAY,NORDER,DET)
	DOUBLE PRECISION ARRAY,AMAX,SAVE
	parameter ( maxfree = 500 )
	DIMENSION ARRAY(norder,norder)
	integer*4 IK(maxfree), JK(maxfree)
	DET = 0.
	IF( NORDER .LE. 0 ) RETURN
	IF( NORDER .GT. MAXFREE ) RETURN
C
10	DET=1.
11	DO 100 K=1,NORDER
C
C FIND LARGEST ELEMENT ARRAY(I,J) IN REST OF MATRIX
C
	AMAX=0. 
21	DO 30 I=K,NORDER
	DO 30 J=K,NORDER
23	IF (DABS(AMAX)-DABS(ARRAY(I,J))) 24,24,30
24	AMAX=ARRAY(I,J) 
	IK(K)=I 
	JK(K)=J 
30	CONTINUE
C
C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K,K)
C
31	IF (AMAX) 41,32,41
32	DET=0.
	GOTO 140
41	I=IK(K) 
	IF (I-K) 21,51,43
43	DO 50 J=1,NORDER
	SAVE=ARRAY(K,J) 
	ARRAY(K,J)=ARRAY(I,J)
50	ARRAY(I,J)=-SAVE
51	J=JK(K) 
	IF (J-K) 21,61,53
53	DO 60 I=1,NORDER
	SAVE=ARRAY(I,K) 
	ARRAY(I,K)=ARRAY(I,J)
60	ARRAY(I,J)=-SAVE
C
C ACCUMULATE ELEMENTS OF INVERSE MATRIX 
C
61	DO 70 I=1,NORDER
	IF (I-K) 63,70,63
63	ARRAY(I,K)=-ARRAY(I,K)/AMAX
70	CONTINUE
71	DO 80 I=1,NORDER
	DO 80 J=1,NORDER
	IF (I-K) 74,80,74
74	IF (J-K) 75,80,75
75	ARRAY(I,J)=ARRAY(I,J)+ARRAY(I,K)*ARRAY(K,J)
80	CONTINUE
81	DO 90 J=1,NORDER
	IF (J-K) 83,90,83
83	ARRAY(K,J)=ARRAY(K,J)/AMAX
90	CONTINUE
	ARRAY(K,K)=1./AMAX
100	DET=DET*AMAX
C
C RESTORE ORDERING OF MATRIX
C
101	DO 130 L=1,NORDER
	K=NORDER-L+1
	J=IK(K) 
	IF (J-K) 111,111,105
105	DO 110 I=1,NORDER
	SAVE=ARRAY(I,K) 
	ARRAY(I,K)=-ARRAY(I,J)
110	ARRAY(I,J)=SAVE 
111	I=JK(K) 
	IF (I-K) 130,130,113
113	DO 120 J=1,NORDER
	SAVE=ARRAY(K,J) 
	ARRAY(K,J)=-ARRAY(I,J)
120	ARRAY(I,J)=SAVE 
130	CONTINUE
140	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
C* DUMP1D ... Dump NCOL arrays to a free-format ascii disk file.
C+
      SUBROUTINE DUMP1D( FILE, NCOL, NROW
     &	, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10
     &	, X11, X12, X13, X14, X15, X16, X17, X18, X19, X20 )
* Input:
*	FILE	C* file name or ' ' for interactive selection
*	NCOL	I4 Number of data arrays
*	NROW	I4 Number of data values in each array
*	X1-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
* 2007 Oct KDH @ Christchurch - arrays 8-20
* 2007 Oct KDH @ Christchurch - PARAMETER( IOUT = 93 )
* 2007 Oct KDH @ Christchurch -- WRITE(RECORD,*), then WRITE(IOUT,'(A)')
*	       since WRITE(IOUT,*) sometimes writes more than 1 line
	REAL*4 X1(*), X2(*), X3(*), X4(*), X5(*)
	REAL*4 X6(*), X7(*), X8(*), X9(*), X10(*)
	REAL*4 X11(*), X12(*), X13(*), X14(*), X15(*)
	REAL*4 X16(*), X17(*), X18(*), X19(*), X20(*)
	CHARACTER*(*) FILE
	CHARACTER*128 OFILE
	CHARACTER*512 RECORD
	CHARACTER*10 REPLY
	PARAMETER ( MAXCOL = 10 )
	PARAMETER ( IOUT = 93 )
	GOTO 100
* multiple entries because unix/linux can't handle more arrays in
* the subroutine argument list than in the calling program
	ENTRY DUMP1D20( FILE, NCOL, NROW
     &		, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10
     &		, X11, X12, X13, X14, X15, X16, X17, X18, X19, X20 )
	GOTO 100
	ENTRY DUMP1D19( FILE, NCOL, NROW
     &		, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10
     &		, X11, X12, X13, X14, X15, X16, X17, X18, X19 )
	GOTO 100
	ENTRY DUMP1D18( FILE, NCOL, NROW
     &		, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10
     &		, X11, X12, X13, X14, X15, X16, X17, X18 )
	GOTO 100
	ENTRY DUMP1D17( FILE, NCOL, NROW
     &		, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10
     &		, X11, X12, X13, X14, X15, X16, X17 )
	GOTO 100
	ENTRY DUMP1D16( FILE, NCOL, NROW
     &		, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10
     &		, X11, X12, X13, X14, X15, X16 )
	GOTO 100
	ENTRY DUMP1D15( FILE, NCOL, NROW
     &		, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10
     &		, X11, X12, X13, X14, X15 )
	GOTO 100
	ENTRY DUMP1D14( FILE, NCOL, NROW
     &		, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10
     &		, X11, X12, X13, X14 )
	GOTO 100
	ENTRY DUMP1D13( FILE, NCOL, NROW
     &		, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10
     &		, X11, X12, X13 )
	GOTO 100
	ENTRY DUMP1D12( FILE, NCOL, NROW
     &		, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10, X11, X12 )
	GOTO 100
	ENTRY DUMP1D11( FILE, NCOL, NROW
     &		, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10, X11 )
	GOTO 100
	ENTRY DUMP1D10( FILE, NCOL, NROW
     &		, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10 )
	GOTO 100
	ENTRY DUMP1D9( FILE, NCOL, NROW
     &		, X1, X2, X3, X4, X5, X6, X7, X8, X9 )
	GOTO 100
	ENTRY DUMP1D8( FILE, NCOL, NROW
     &		, X1, X2, X3, X4, X5, X6, X7, X8 )
	GOTO 100
	ENTRY DUMP1D7( FILE, NCOL, NROW, X1, X2, X3, X4, X5, X6, X7 )
	GOTO 100
	ENTRY DUMP1D6( FILE, NCOL, NROW, X1, X2, X3, X4, X5, X6 )
	GOTO 100
	ENTRY DUMP1D5( FILE, NCOL, NROW, X1, X2, X3, X4, X5 )
	GOTO 100
	ENTRY DUMP1D4( FILE, NCOL, NROW, X1, X2, X3, X4 )
	GOTO 100
	ENTRY DUMP1D3( FILE, NCOL, NROW, X1, X2, X3 )
	GOTO 100
	ENTRY DUMP1D2( FILE, NCOL, NROW, X1, X2 )
	GOTO 100
	ENTRY DUMP1D1( FILE, NCOL, NROW, X1 )

* 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 data
      NR = 0
      DO I=1,NROW
      IF( NC .EQ. 1 ) THEN
	WRITE( RECORD, *, ERR=98 ) X1(I)
      ELSE IF( NC .EQ. 2 ) THEN
	WRITE( RECORD, *, ERR=98 ) X1(I), X2(I)
      ELSE IF( NC .EQ. 3 ) THEN
	WRITE( RECORD, *, ERR=98 ) X1(I), X2(I), X3(I)
      ELSE IF( NC .EQ. 4 ) THEN
	WRITE( RECORD, *, ERR=98 ) X1(I), X2(I), X3(I), X4(I)
      ELSE IF( NC .EQ. 5 ) THEN
	WRITE( RECORD, *, ERR=98 ) X1(I), X2(I), X3(I), X4(I), X5(I)
      ELSE IF( NC .EQ. 6 ) THEN
	WRITE( RECORD, *, ERR=98 ) X1(I), X2(I), X3(I), X4(I), X5(I)
     &	, X6(I)
      ELSE IF( NC .EQ. 7 ) THEN
	WRITE( RECORD, *, ERR=98 ) X1(I), X2(I), X3(I), X4(I), X5(I)
     &	, X6(I), X7(I)
      ELSE IF( NC .EQ. 8 ) THEN
	WRITE( RECORD, *, ERR=98 ) X1(I), X2(I), X3(I), X4(I), X5(I)
     &	, X6(I), X7(I), X8(I)
      ELSE IF( NC .EQ. 9 ) THEN
	WRITE( RECORD, *, ERR=98 ) X1(I), X2(I), X3(I), X4(I), X5(I)
     &	, X6(I), X7(I), X8(I), X9(I)
      ELSE IF( NC .EQ. 10 ) THEN
	WRITE( RECORD, *, ERR=98 ) X1(I), X2(I), X3(I), X4(I), X5(I)
     &	, X6(I), X7(I), X8(I), X9(I), X10(I)
      END IF
	WRITE( IOUT, '(A)', ERR=98 ) RECORD( : LEN1( RECORD ) )
        NR = I
      END DO

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

* file write error
 98   WRITE(*,*) '**DUMP1D: error writing to file ', OFILE(:LEN1(OFILE))
      GOTO 50
       
* file open error
 99   WRITE(*,*) '**DUMP1D: error opening file ', OFILE(:LEN1(OFILE))
      OFILE = ' '
      GOTO 10
      END
      SUBROUTINE DUMP2D( FILE, NIM, NX, NY, X, Y,
     #  D1, D2, D3, D4, D5 )
*    
* Load one or more 2D images to a disk file
*
* Input:
*	FILE	C* file into which to dump data
*	NIM	I4 number of images
*	NX,NY	I4 actual number of pixels
*	X(NX)	R4 X values
*	Y(NY)	R4 Y values
*	D1(NX,NY)	R4 data values in first image
*	D2,D3...	R4 data values in subsequent images
*
* Apr 1988 Keith Horne @ STScI
* Aug 1990 KDH @ STScI - add nim and multi-image capability
* Apr 1991 KDH @ STScI - fix bug on file close found by J.Krolik
* Jan 1994 KDH @ Utrecht - MAXIM=5, DUMP_2D -> DUMP2D
* Mar 1994 KDH @ Utrecht - entry points for multiple image dumps
* 2000 Nov KDH @ Austin - open type=unknown
* 2001 Aug KDH @ Austin - use OPENFILE
	CHARACTER*(*) FILE
	CHARACTER*256 OUTFILE
	PARAMETER (MAXIM=5)
	REAL*4 X(*), Y(*)
	REAL*4 D1(*), D2(*), D3(*), D4(*), D5(*)

* entry points since UNIX needs exact match to number of arguments
	GOTO 100
	ENTRY DUMP2D5( FILE, NIM, NX,NY, X,Y, D1, D2, D3, D4, D5 )
	GOTO 100
	ENTRY DUMP2D4( FILE, NIM, NX,NY, X,Y, D1, D2, D3, D4 )
	GOTO 100
	ENTRY DUMP2D3( FILE, NIM, NX,NY, X,Y, D1, D2, D3 )
	GOTO 100
	ENTRY DUMP2D2( FILE, NIM, NX,NY, X,Y, D1, D2 )
	GOTO 100
	ENTRY DUMP2D1( FILE, NIM, NX,NY, X,Y, D1 )

* check dimensions

  100   IF( NIM.LE.0 .OR. NX.LE.0 .OR. NY.LE.0 ) THEN
	WRITE(*,*), '** Invalid dimensions', NIM, NX, NY
	GOTO 999
      ELSE IF( NIM.GT.MAXIM ) THEN
	WRITE(*,*) '** Too many images', NIM, MAXIM
      END IF

	OUTFILE = FILE
      IF( FILE.EQ.' ' ) THEN
	WRITE(*,'(A,$)') ' Output file for 2D data ? '
	READ(*,'(A)') OUTFILE
      END IF

* open the file

c	IUNIT = 57
c	OPEN( UNIT=IUNIT,
c    #		NAME=FILE,
c    #		FORM='UNFORMATTED',
c    #		TYPE='UNKNOWN',
c    #		ERR=996 )

	CALL OPENFILE( IUNIT, FILE, 'UNFORMATTED', 'UNKNOWN', IER )
	IF( IER .NE. 0 ) GOTO 996

* write the data

	WRITE( IUNIT, ERR=997 ) NX, NY, NIM
	WRITE( IUNIT, ERR=997 ) (X(I),I=1,NX)
	WRITE( IUNIT, ERR=997 ) (Y(I),I=1,NY)
	IF( NIM.GE.1 )
     #	WRITE( IUNIT, ERR=997 ) (D1(I),I=1,NX*NY)
	IF( NIM.GE.2 )
     #	WRITE( IUNIT, ERR=997 ) (D2(I),I=1,NX*NY)
	IF( NIM.GE.3 )
     #	WRITE( IUNIT, ERR=997 ) (D3(I),I=1,NX*NY)
	IF( NIM.GE.4 )
     #	WRITE( IUNIT, ERR=997 ) (D4(I),I=1,NX*NY)
	IF( NIM.GE.5 )
     #	WRITE( IUNIT, ERR=997 ) (D5(I),I=1,NX*NY)

        CLOSE( IUNIT )
 1000	RETURN

* error messages

  996	WRITE(*,*) '** Error on file open.'
	GOTO 999
  997	WRITE(*,*) '** Error on file write.'
	CLOSE( UNIT=IUNIT )
	GOTO 999
  999	WRITE(*,*) '** DUMP2D aborted.'
	RETURN
	END
***********************************
	subroutine peqc( n, x, c )
* x(i) = c for i=1...n
	real*4 x(n)
      do i=1,n
	x(i) = c
      end do
	return
	end
***********************************
	subroutine peqp( n, x, y )
* x(i) = y(i) for i=1...n
	real*4 x(n), y(n)
      do i=1,n
	x(i) = y(i)
      end do
	return
	end
	subroutine peqpac( n, x, y, c )
* x(i) = y(i) + c for i=1...n
	real*4 x(n), y(n)
	if( c .eq. 0. ) return
      do i=1,n
	x(i) = y(i) + c
      end do
	return
	end
***********************************
	subroutine peqpmc( n, x, y, c )
* x(i) = y(i) for i=1...n
	real*4 x(n), y(n)
	if( c .eq. 1.0 ) return
      do i=1,n
	x(i) = y(i) * c
      end do
	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 @ St.And - linux port: must use g77 -fno-automatic
* 2001 Jul KDH @ St.And - simplify, strip vms & multi-sites
* 2001 Aug KDH @ St.And - FULLPATH plot files, ASPECT
* 2004 Jan KDH @ St.And - add err= to read (linux)
* 2005 Feb KDH @ St.And - strip @ from plotfile
* 2006 May KDH @ St.And - add farside printer
* 2010 Jun KDH @ St.And - update StA printers

        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 (?,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 -Pkdh'
*      else if( reply .eq. 'S') then
*	prtcom = 'lpr -Pstar-lw'
      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 -PC544'
      else
        write(*,*) 'St.Andrews Printer queues :'
        write(*,*) ' K  kdh     (Keith Horne''s office)'
        write(*,*) ' D  duplex  double-sided (Starlink room)'
        write(*,*) ' T  T644N   double-sided (south corridor)'
        write(*,*) ' C  C544    colour (Starlink room)'
	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

	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
      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
C*AMOEBA ... downhill simplex fitting program
C+
	SUBROUTINE AMOEBA(P,Y,MP,NP,NDIM,FTOL,FUNK,ITER)
*
* adjust parameters to minimize value of a function
*
* input:
*	P	R4(MP,NP)	simplex vertices
*	Y	R4(NDIM+1)	function value at each vertex
*	MP	I4		max number of vertices
*	NP	I4		max number of parameters
*	NDIM	I4		actual number of parameters
*	FTOL	R4		relative tolerance for convergence
*	ITER	I4		maximum iterations
* output:
*	P	R4(MP,NP)	final vertices of simplex
*	Y	R4(NDIM+1)	FUNK values at each vertex
*	ITER	I4		number of iterations
C--
*
* May 1989 Keith Horne @ STScI - adapted from Numerical Recipes
* May 1989 KDH @ STScI - modified convergence test ala Kip Kuntz.
* Jun 1989 KDH @ STScI - convergence messages to terminal.
* 2002 Nov KDH @ St.And - report iter
* 2006 May KDH @ St.And - protect RTOL divide by zero
* 2007 Apr KDH @ St.And - quit if RTOL = 0.
* 2007 Apr KDH @ St.And - param span convergence test bugs fixed.
* 2008 Jul KDH @ SAAO - common/muteamoeba/
* 2009 Jul KDH @ St.A - entry amoebaplot convergence plot
* 2009 Jul KDH @ St.A - converge on (Ptol < Ftol or Rtol=0) + extra iterations
	PARAMETER ( NMAX = 1000 )
	DIMENSION P(MP,NP), Y(MP)
	DIMENSION PR( NMAX ), PRR( NMAX ), PBAR( NMAX )
	LOGICAL CONTRACT
	common/muteamoeba/ mute
* plot
	parameter ( maxplot = 10000 )
	real*4 plot1( maxplot )
	real*4 plot2( maxplot )
	real*4 plot3( maxplot )
	real*4 plot4( maxplot )
	real*4 plot5( maxplot )
	real*4 plot6( maxplot )
	character*(*) atitle
	character*75 title, xlabel, ylabel
	logical pgagain
	data kount / 0 /

* amoeba plot ?
	makeplot = 1

* sanity check inputs
	if( np .lt. 1 ) return
	if( ndim .lt. 1 ) return
	if( mp .lt. ndim+1 ) return
	if( iter .lt. 1 ) return

* parameters
	ALPHA = 1.0
	BETA = 0.5
	GAMMA = 2.0
	MPTS = NDIM + 1
	ITMAX = ITER
	ITER = 0

* extra iterations
	NEXTRA = MIN( ITMAX, MPTS )
	IEXTRA = NEXTRA

* needed for diagnostic plot
	nplot = 0
	needlogs = 1
	ctol = ftol
	npar = ndim

* find lowest, hightest, and next-highest
1	ILO=1
	IF(Y(1).GT.Y(2))THEN
		IHI=1
		INHI=2
	ELSE
		IHI=2
		INHI=1
	ENDIF
	DO I=1,MPTS
		IF(Y(I).LT.Y(ILO)) ILO=I
		IF(Y(I).GT.Y(IHI))THEN
			INHI=IHI
			IHI=I
		ELSE IF(Y(I).GT.Y(INHI))THEN
			IF(I.NE.IHI) INHI=I
		ENDIF
	ENDDO

* RTOL is fractional range of Chi^2 ( or -2 ln( P ) )
* note this assumes that Y is positive.
* true for Chi^2 but not true in general, e.g. -2lnL.
* not sure how to do better.
	YHI = Y(IHI)
	YLO = Y(ILO)
	DIF = ABS( YHI - YLO )
	SUM = ABS( YHI ) + ABS( YLO )
      IF( SUM .NE. 0. ) THEN
	RTOL = DIF / SUM
      ELSE
	WRITE(*,*) '** AMOEBA DIVIDE BY ZERO AVERTED. RTOL=', DIF, ' / ',  SUM
	RTOL = 2. * ABS( FTOL ) 
      END IF

* PTOL is max fractional parameter span
	PTOL = 0.
	IPWORST = 0
	TINY = MIN( 1.E-10, FTOL / 100. )
      DO I=1,NDIM
	JHI = 1
	JLO = 1
      DO J=1,MPTS
	IF( P(J,I) .GT. P(JHI,I) ) JHI = J
	IF( P(J,I) .LT. P(JLO,I) ) JLO = J
      END DO
	PHI = P(JHI,I)
	PLO = P(JLO,I)
	SUM = ABS( PHI ) + ABS( PLO )
	DIF = ABS( PHI - PLO )
	TEST = DIF / MAX( SUM, TINY )
      IF( TEST .GT. PTOL ) THEN
	PTOL = MAX( PTOL, TEST )
	IPWORST = I
      END IF
      END DO

* stow for plot
      if( makeplot .ne. 0 ) then
	nplot = nplot + 1
	if( nplot .gt. maxplot ) nplot = nplot - maxplot
	plot1( nplot ) = iter
	plot2( nplot ) = yhi
	plot3( nplot ) = ylo
	plot4( nplot ) = yhi - ylo
	plot5( nplot ) = rtol
	plot6( nplot ) = ptol
      if( nplot .eq. 1 ) then
	bofhi = yhi
	boflo = ylo
      else
	bofhi = max( bofhi, yhi )
	boflo = min( boflo, ylo )
      end if
      end if

* converge if small chi^2 range and last operation is contraction

      IF( RTOL .LT. FTOL ) then
      IF( CONTRACT .OR. RTOL .LE. 0. ) THEN
	if( mute .eq. 0 ) WRITE(*,*) 'Amoeba CONVERGED.'
     &		, ' PTOL', PTOL, ' RTOL', RTOL, ' <', FTOL, ' :)'
c	GOTO 1001
      ELSE
C	if( mute .eq. 0  ) WRITE(*,*) 'Amoeba its', ITER,
C     &		, ' PTOL', PTOL, ' RTOL', RTOL, ' <', FTOL, ' NOT contracting. :('
      END IF
      END IF

* converge if small range in all parameters

c      IF( PTOL .LT. FTOL .AND. RTOL .LT. FTOL ) THEN
      IF( PTOL .LT. FTOL .OR. RTOL .EQ. 0. ) THEN
	if( mute .eq. 0 ) WRITE(*,*) 'Amoeba CONVERGED.'
     &		, ' RTOL', RTOL, ' PTOL', PTOL, ' <', FTOL, ' :)'
	GOTO 1001
      END IF

* converge if maximum iterations exceeded

      IF( ITER .GE. ITMAX ) THEN
      IF( ITMAX .GT. 1 ) THEN
	if( mute .eq. 0 ) WRITE(*,*) 'Amoeba Iter', ITER, ' max', ITMAX
     &		, ' PTOL', PTOL, ' RTOL', RTOL, ' >', FTOL, ' :('
	GOTO 1001
      END IF
	GOTO 1000
      END IF

* iterate
   10	CONTRACT = .FALSE.
	ITER=ITER+1
	DO J=1,NDIM
		PBAR(J)=0.
	END DO
	DO I=1,MPTS
		IF(I.NE.IHI)THEN
			DO J=1,NDIM
				PBAR(J)=PBAR(J)+P(I,J)
			ENDDO
		ENDIF
	END DO
	DO J=1,NDIM
		PBAR(J)=PBAR(J)/NDIM
		PR(J)=(1.+ALPHA)*PBAR(J)-ALPHA*P(IHI,J)
	END DO
	YPR=FUNK(PR)
	IF(YPR.LE.Y(ILO))THEN
		DO J=1,NDIM
			PRR(J)=GAMMA*PR(J)+(1.-GAMMA)*PBAR(J)
		END DO
		YPRR=FUNK(PRR)
		IF(YPRR.LT.Y(ILO))THEN
			DO J=1,NDIM
				P(IHI,J)=PRR(J)
			END DO
			Y(IHI)=YPRR
		ELSE
			DO J=1,NDIM
				P(IHI,J)=PR(J)
			END DO
			Y(IHI)=YPR
		END IF
	ELSE IF(YPR.GE.Y(INHI))THEN
		IF(YPR.LT.Y(IHI))THEN
			DO J=1,NDIM
				P(IHI,J)=PR(J)
			END DO
			Y(IHI)=YPR
		END IF
		DO J=1,NDIM
			PRR(J)=BETA*P(IHI,J)+(1.-BETA)*PBAR(J)
		END DO
		YPRR=FUNK(PRR)
		IF(YPRR.LT.Y(IHI))THEN
			DO J=1,NDIM
				P(IHI,J)=PRR(J)
			END DO
			Y(IHI)=YPRR
		ELSE
			contract = .true.
			DO I=1,MPTS
				IF(I.NE.ILO)THEN
					DO J=1,NDIM
						PR(J)=0.5*(P(I,J)+P(ILO,J))
						P(I,J)=PR(J)
					END DO
					Y(I)=FUNK(PR)
				END IF
			END DO
		END IF
	ELSE
		DO J=1,NDIM
			P(IHI,J)=PR(J)
		END DO
		Y(IHI)=YPR
	END IF

	GOTO 1


* report vertices

 1001 if( mute .eq. 0 ) then
	WRITE(*,*) 'Amoeba its', ITER, ' RTOL', rtol, ' PTOL', ptol
c      do i=1,mpts
c	write(*,*) i, y(i), (p(i,j),j=1,ndim)
c      end do
      end if

* extra iterations
	iextra = iextra - 1
	if( iextra .ge. 0 ) goto 10

* normal return

 1000	RETURN

* diagnostic plot

	entry amoebaplot( atitle )

	write(*,*) '** amoebaplot **********', nplot
	if( nplot .le. 0 ) return

* logs of plot diagnostics

      if( needlogs .ne. 0 ) then
	kount = kount + 1
	floor = max( 1.e-12, ctol / 100. )
      do i=1,nplot
	plot2(i) = alog10( max( floor, abs( plot2(i) ) ) )
	plot3(i) = alog10( max( floor, abs( plot3(i) ) ) )
	plot4(i) = alog10( max( floor, abs( plot4(i) ) ) )
	plot5(i) = alog10( max( floor, abs( plot5(i) ) ) )
	plot6(i) = alog10( max( floor, abs( plot6(i) ) ) )
      end do
	needlogs = 0
      end if

* format
  100	nx = 1
	ny = 1
* open
	call pgstart( nx, ny, new )
	csize = 1.5
	call pgsch( csize )
* window
	call mnmx( nplot, plot1, xmn, xmx ) 
	call mnmx( nplot, plot2, ymn, ymx ) 
	call mnmx( nplot, plot3, yymn, yymx ) 
	call mnmx( nplot, plot4, yyymn, yyymx ) 
	ymn = min( ymn, yymn, yyymn, alog10( ctol ) )
	split = 2.
	ymn = nint( ymn/split - 0.5 )*split
	ymx = nint( ymx/split + 0.5 )*split
	call pgsci( 5 )
	logy = 0
	if( ymx - ymn .le. 10. ) logy = 20
	call pgenv( xmn, xmx, ymn, ymx, 0, logy )
* mart itmax
	call pgsls( 4 )
	x = itmax
	call pgmove( x, ymn )
	call pgdraw( x, ymx )
* mark Ftol
	yy = alog10( ctol )
	call pgmove( xmn, yy )
	call pgdraw( xmx, yy )
	call pgsls( 1 )
	xpc = ( xmx - xmn ) / 100.
	x = xmn + 2. * xpc
	yy = alog10( ctol ) + 0.01 * (ymx - ymn )
	call pgtext( x, yy, 'Ftol' )
* labels
	title = atitle
      if( title .eq. ' ' ) then
	title = 'Amoeba'
	call append_i( title, kount )
      end if
	xlabel = 'Iterations'
	i = nint( plot1( nplot ) )
	call append_i( xlabel, i )
	ylabel = 'Badness of Fit'
	if( logy .eq. 0 ) ylabel = 'log[ ' // ylabel(:len1(ylabel)) // ' ]' 
	call pgsci( 7 )
	call pglabel( xlabel, ylabel, title )
* number of params
	title = 'N\dP\u'
	call append_i( title, npar )
	call pgmtxt( 'T', 0.5, 0.0, 0.0, title )
* badness of fit
	title = 'BoF'
	call append_f( title, bofhi, '(f15.2)' )
	call pgmtxt( 'T', 0.5, 1.0, 1.0, title )
	title = ' '
	call append( title, '=>' )
	call append_f( title, boflo, '(f15.2)' )
	call pgmtxt( 'T', -1.7, 0.93, 1.0, title )
* BoF (max and min)
	call pgsci( 1 )
	call pgline( nplot, plot1, plot2 )
	call pgline( nplot, plot1, plot3 )
	i = max( 1, nplot - nextra - 2 )
	x = xmx + xpc
	yy = max( plot3(i), ymn )
	call pgtext( x , yy, 'BoF' )

* Bof range
	call pgsci( 7 )
	call pgline( nplot, plot1, plot4 )
	yy = max( plot4(i), ymn )
	call pgtext( x , yy, '\gDBoF' )
* Ptol
	call pgsci( 3 )
	call pgsls( 2 )
	call pgline( nplot, plot1, plot5 )
	yy = max( plot5(i), ymn )
	call pgtext( x , yy, 'Rtol' )
* Rtol
	call pgsci( 4 )
	call pgsls( 3 )
	call pgline( nplot, plot1, plot6 )
	yy = max( plot6(i), ymn )
	title = 'Ptol'
	if( ipworst .gt. 0 ) call append_i( title, ipworst )
	call pgtext( x , yy, title )
	call pgsls( 1 )
* close
	call pgstop
* again
	if( pgagain() ) goto 100
	goto 1000

	END
*-------------------------------------------------------
	SUBROUTINE dZROOTS( A, M, ROOTS, POLISH )
* 2008 Keith Horne @ St.Andrews -- double precision
	implicit double precision ( a-h, o-z )
	PARAMETER ( EPS = 1.d-15, MAXM = 101 )
	double COMPLEX A( M + 1 ), ROOTS( M ), AD( MAXM ), X, B, C
	LOGICAL POLISH
      DO J = 1, M + 1
	AD( J ) = A( J )
      END DO
      DO J = M, 1, -1
	X = dCMPLX( 0., 0. )
	CALL dLAGUER( AD, J, X, EPS, .FALSE. )
	big = 2.d0 * eps * eps * abs( dreal( x ) )
	IF( ABS( dIMAG( X ) ) .LE. big ) X = dCMPLX( dREAL( X ), 0.d0 )
	ROOTS( J ) = X
	B = AD( J + 1 )
      DO JJ = J, 1, -1
	C = AD( JJ )
	AD( JJ ) = B
	B = X * B + C
      END DO
      END DO
      IF( POLISH ) THEN
      DO J = 1, M
	CALL dLAGUER( A, M, ROOTS( J ), EPS, .TRUE. )
      END DO
      END IF
c don't re-order the roots
c      DO J = 2, M
c	X = ROOTS( J )
c      DO I = J-1, 1, -1
c	IF( dREAL( ROOTS( I ) ) .LE. dREAL( X ) ) GO TO 10
c	ROOTS( I + 1 ) = ROOTS( I )
c      END DO
c	I = 0
c   10	ROOTS( I + 1 ) = X
c      END DO
	RETURN
	END
*---------------------------------------------------------
	SUBROUTINE dLAGUER( A, M, X, EPS, POLISH )
* 2008 Keith Horne @ St.Andrews -- double precision
	implicit double precision (a-h,o-z)
	double COMPLEX A(M+1),X,DX,X1,B,D,F,G,H,SQ,GP,GM,G2,ZERO
	LOGICAL POLISH
	PARAMETER( ZERO = ( 0.d0, 0.d0 ) )
	parameter( EPSS = 6.d-15, MAXIT = 100 )
	DXOLD = CdABS( X )
      DO ITER = 1, MAXIT
	B = A( M + 1 )
	ERR = CdABS( B )
	D = ZERO
	F = ZERO
	ABX = CdABS( X )
      DO J = M, 1, -1
	F = X * F + D
	D = X * D + B
	B = X * B + A( J )
	ERR = CdABS( B ) + ABX * ERR
      END DO
	ERR = EPSS * ERR
      IF( CdABS( B ) .LE. ERR ) THEN
	DX = ZERO
	RETURN
      ELSE
	G = D / B
	G2 = G * G
	H = G2 - ( F + F ) / B
	SQ = CdSQRT( ( M - 1 ) * ( M * H - G2 ) )
	GP = G + SQ
	GM = G - SQ
* trick from Dave Bennett's version
	abp = cdabs( gp )
	abm = cdabs( gm )
      IF( abp .LT. abm ) GP = GM
      if( max( apb, abm ) .gt. 0.) then
	DX = M / GP
      else
	dx = exp( dcmplx( log( 1.d0 + abx ), dble( iter ) ) ) 
      end if
      END IF
	X1 = X - DX
	itns = iter
	IF( X .EQ. X1 ) RETURN
	X = X1
	CDX = CdABS( DX )
      IF( ITER .GT. 6. AND. CDX .GE. DXOLD ) THEN
c	write(*,*) 'LAGUER diverging.  :(((.'
	RETURN
      END IF
	DXOLD = CDX
      IF( .NOT. POLISH ) THEN
	IF( CdABS( DX ) .LE. EPS * CdABS( X ) ) RETURN
      END IF
      END DO

	write(*,*) 'its:', itns, ' x:', x, ' dx:', dx
	write(*,*) 'm:', m, ' b:', b
	write(*,*) 'a:', (a(i),i=1,m+1)
	STOP 'LAGUER too many iterations'
	RETURN
	END
