c************
c
c   Description :  Utility routines for Midas libraries
c
c   Author :       JGS  4/4/90
c
c   $Revision: 1.5 $
c   6-Jun-98 RDArendt: Added Mc.visClass, Mc.visDepth, and Mc.visID
c
c
c   Modified:  [July 31 2002 ajb]
c     DR# 10552 - CALC result truncated in display
c
c     2020-06026 Updates from Duy contributed as part of adding GFC support NXM-3531
c
c
c************

	function M$SEARCH (string,sub)
c
	character*(*) string,sub
c
c use fortran intrinsic to find placement of substring
	M$SEARCH = INDEX(string,sub)
c
	return
	end



	function M$SEARCHB (string,sub)

	implicit none

	include 'types.inc'
c Arguments
	character*(*)	string, sub
c Locals
	integer*4	i, L1, L2
c Code
c
c find placement of substring searching from the back
	entry INDEXB (string,sub)
	L1 = LEN(string)
	L2 = LEN(sub)-1
	do i = L1-L2, 1, -1
	  if (string(i:i+L2).eq.sub) goto 10
	enddo
c
	INDEXB    = 0
	M$SEARCHB = 0
	return
c
 10	INDEXB    = i
	M$SEARCHB = i
	return
c
	end



	function M$LENGTHQ (string)
c
	character*(*) string
c
c find length of string up to the trailing blanks - fast
c start at front looking for at least 16 blanks in succession
	i = 1
	ls = LEN(string)
  1	if (i.eq.ls) goto 9
	i = min(ls,i+8)
	if (string(i:i).gt.' ') goto 1
	i1 = i
	i2 = min(ls,i+8)
	do i = i2,i1,-1
	  if (string(i:i).gt.' ') goto 1
	enddo
	do i = i1-1,1,-1
	  if (string(i:i).gt.' ') goto 9
	enddo
	i = 0
  9	M$LENGTHQ = i
	return
	end



	function M$PARSE (string,ind,ostring)
c
	character*(*) string,ostring
c
c skip by ind-1 delimitted strings
	i = 1
	n = 1
	il = 0
	ls = LEN(string)
	do while (n.lt.ind)
	  i = i + M$LENGTHD(string(i:))
#if _IRIX
	  if (i.gt.ls) goto 10
	  if (string(i:i).eq.CHAR(0)) goto 10
#else
	  if (i.gt.ls .or. string(i:i).eq.CHAR(0)) goto 10
#endif
	  if (string(i:i).eq.',') n = n+1
	  i = i+1
	enddo
c
c pull out INDth delimitted string
	il = M$LENGTHD(string(i:))
	do while (string(i+il:i+il).eq.' ')
	  if (i+il.ge.ls) then
	    if (string(i:ls).eq.' ') il = 0
	    goto 10
	  endif
	  il = il+1 + M$LENGTHD(string(i+il+1:))
	enddo
c
c fill in the output string and return value
  10    if (il.eq.0) then
	  ostring = ' '
	elseif (string(i:i).eq.'"') then
	  ! tag quoted strings by negating return length (trim off quotes)
	  il = M$TRIM_QUOTES (string(i:i+il-1),ostring)
	  if (il .lt. LEN(ostring)) ostring(il+1:) = ' '
	  il = -il
	else
	  ostring = string(i:i+il-1)
	endif
	M$PARSE = il
	return
	end

	function M$TRIM_QUOTES (sq,sr)
c
	character*(*) sq,sr
c
	iq = 0
	ir = 0
	lq = LEN(sq)
	lr = LEN(sr)
	if (sq(1:1) .ne. '"') then
	  call M$WARNING('Assuming leading quote')
	else
	  iq = iq + 1
	endif
	do while (iq.lt.lq)
	  iq = iq+1
	  if (sq(iq:iq).eq.'"') then
	    if (iq.lt.lq .and. sq(iq+1:iq+1).eq.'"') then
	      ir = ir+1
	      if (ir.le.lr) sr(ir:ir) = '"'
	      iq = iq+1
	    endif
	  else
	    ir = ir+1
	    if (ir.le.lr) sr(ir:ir) = sq(iq:iq)
	  endif
	enddo
	M$TRIM_QUOTES = ir
	return
	end



	subroutine M$LOWERCASE(string)
	include 'types.inc'
	character*(*) string
	integer*4 i, j, iascii

	i = 1
	j = M$LENGTHQ(string)
	do i = 1, j
	  iascii = ICHAR(string(i:i))
	  if ((iascii .ge. 65) .and. (iascii .le. 90))
     &		string(i:i) = CHAR(iascii + 32)
	enddo

	return
	end


	subroutine M$UPPERCASE(string)
	include 'types.inc'
	character*(*) string
	integer*4 i, j, iascii

	i = 1
	j = M$LENGTHQ(string)
	do i = 1, j
	  iascii = ICHAR(string(i:i))
	  if ((iascii .ge. 97) .and. (iascii .le. 122))
     &		string(i:i) = CHAR(iascii - 32)
	enddo

	return
	end




	subroutine M$WORDWRAP(str, start, margin, consume, wrap, new)

	implicit none
c Arguments
	character*(*)	str
	integer*4	start, margin
	logical*4	consume
	integer*4	wrap, new		! (OUT parameters)
c Locals
	integer*4	wrapmargin, i, lastOK, col, lstr

	wrapmargin = start + margin - 1
	lstr = LEN(str)
	col  = 1

c Find the point of wrapping
	lastOK = -1
	wrap = wrapmargin + 1
	i    = start
#if _IRIX
	i = MIN(i,lstr)
#endif

	if (str(i:i).eq.CHAR(9) .and. MOD(col, 8).ne.0)
     &	    col = col + (8 - MOD(col, 8))
	do while (i.le.lstr .and. col.lt.margin)
	  if      (str(i:i) .eq. ',') then
	    wrap   = i
	    lastOK = i
	  else if (str(i:i) .gt. ' ') then
	    lastOK = i

c Special case:  if we're at the last character in the string, and it's not
c 		 whitespace, allow the wrap point to occur here.

	    if (i .eq. lstr) wrap = lastOK
	  else
	    if (.not. consume) then
	      wrap = i
	    else if (lastOK .gt. 0) then
	      wrap = lastOK
	    endif
	  endif
	  i   = i   + 1
	  col = col + 1
	  if (str(i:i).eq.CHAR(9) .and. MOD(col, 8).ne.0)
     &	    col = col + (8 - MOD(col, 8))
	enddo

	if ( col.lt.margin .or.		! no wrapping required
     &	    (col.eq.margin .and. wrapmargin.ge.lstr) ) then
	  wrap = lstr
	  new  = -1
	  return
	endif

	new = wrap + 1

c If wrap is > wrapmargin, we have an overlength line.
	if (wrap .gt. wrapmargin) then
	  wrap = wrapmargin - 1
	  new  = wrapmargin
	endif

	if (new .gt. lstr) new = -1
	if (new .lt. 0 .or. .not. consume) return

c The user has asked that intermediate white space be consumed,
c and furthermore, 'new' is presently .le. lstr.
c Now find where the next string should begin.
	do while (new .le. lstr)
	  if (str(new:new) .gt. ' ') return
	  new = new + 1
	enddo

c We could not find a new place to start, because all the rest of the string
c is white space.  Since we're here, we know the user has asked to consume
c it... therefore tell them there is nothing left to wrap.
	new = -1

	return
	end




	subroutine M$WRAPNEXTLINES(lun, str)

	implicit none
c Arguments
	integer*4	lun
	character*(*)	str
c Locals
	integer*4	orig_margin
	parameter(orig_margin=72)
	integer*4	hang_margin
	parameter(hang_margin=orig_margin-8)

	character*8	fmt

	integer*4	j, new, wrap, margin
c Code
c
c Trick with M$WORDWRAP:  saying 'consume' will cause 'new' to be advanced
c to the next nonwhite.  But we want to leave that white space BEFORE the
c ampersand, so we write to new-1 if wrapped.
	j      = 1
	fmt    = '(a)'
	margin = orig_margin
	do while (j .gt. 0)
	  call M$WORDWRAP(str, j, margin,
     &				.TRUE., wrap, new)
	  if (new .le. 0) then
	    write (lun, fmt) str(j:wrap)
	  else
c Careful here!  Make sure that the ampersand M$WRAPNEXTLINES writes will
c not be immediately preceded by one in the wrapped string.  It will be
c ignored by M$NEXTLINE on the way back in.
c This won't be too hard since we know that 'new' is sitting on a non-white
c character.  Either we're on white space or in the middle of an unbreakable
c sequence.
c It is possible to have a degenerate case of nothing but ampersands.
c We will try to prevent it by walking all the way back to j +1 if we have to.
	    do while (new-1.gt.j+1 .and. str(new-1:new-1) .eq. '&')
	      new = new - 1
	    enddo
	    write (lun, fmt) str(j:new-1)//'&'
	  endif
	  if (j .eq. 1) then
	    margin = hang_margin
	    fmt    = '(1h'//CHAR(9)//',a)'
	  endif
	  j = new
	enddo

	return
	end




	function M$CHK (n1,n2,n)
c
	byte n1(1),n2(1)
	integer*4 n
c
	M$CHK = 0
	do i = 1,n
	  if (n1(i).ne.n2(i)) return
	enddo
	M$CHK = 1
c
	return
	end


	function M$D2A (dval,string)
	include 'midas.inc'
	character*(*) string
	real*8 dval
	M$D2A = M$D2A_FORM (dval,string,Mc.nformat)
	return
	end


	function M$D2A_FORM (dval,string,form)
c
	implicit none
c
	include 'types.inc'
c Locals
	character string*(*),form*(*),tmp*21,tmp2*21,tmp3*(L_max_ascii)
	real*8 dval
	integer*4 lf, ls, i
	integer*4 iexp, lexp, kexp
c Functions Called
	integer*4 M$CFORMAT
c Code
	ls = LEN(string)
	lf = M$LENGTH(form)

c If format string starts with paren, assume it's a FORTRAN
c I/O format specifier and immediately write it into <string>.
	if (form(1:1) .eq. '(') then
	  tmp3   = ' '
	  string = ' '
	  i = M$LENGTHQ(form)
	  if (form(i:i).eq.'L') then
	    write(tmp3, form(1:i-1)) INT (dval)
	  else
	    write(tmp3, form) dval
	  endif
	  i = min(M$LENGTHQ(tmp3), ls)
	  string = tmp3(1:i)
	  M$D2A_FORM = i
	  return

c If format string starts with square brace [, assume it's a C
c I/O format specifier and immediately write it into <string>.
	elseif (form(1:1) .eq. '[') then
	  tmp3   = ' '
	  string = ' '
	  i = M$SEARCHB(form, ']')
	  if (i.gt.0) then
#if _IRIX
	    i = min(i,lf)
#endif
	    if (i.lt.lf .and. form(i+1:i+1).eq.'L') then
	      i = M$CFORMAT (INT(dval), 'L', tmp3, form(2:i-1))
	    else
	      i = M$CFORMAT (dval, 'D', tmp3, form(2:i-1))
	    endif
	    i = min(M$LENGTHQ(tmp3), ls)
	    if (i.gt.0) string = tmp3(1:i)
	  else
	    i = -1
	  endif
	  M$D2A_FORM = i
	  return

	elseif (form(1:lf) .eq. 'DMS') then
	  i = min(M$DEG2DMS(dval, tmp3), ls)
	  string = tmp3(1:i)
	  M$D2A_FORM = i
	  return

	elseif (form(1:lf) .eq. 'HMS') then
	  i = min(M$SEC2TOD(dval, tmp3), ls)
	  string = tmp3(1:i)
	  M$D2A_FORM = i
	  return

	endif		! special cases

c create numeric string in exponential notation
	write(tmp,'(1pe21.14)') dval
c
c find the exponent (without using READ for speed)
	iexp = ICHAR(tmp(21:21)) - ICHAR('0')
	if (tmp(20:20).ne.'0') iexp = iexp +
     &		10 * ( ICHAR(tmp(20:20)) - ICHAR('0') )
	if (tmp(19:19).eq.'-') iexp = -iexp
c
c determine the output exponent based on form
	if (form.eq.'GEN') then
	  lexp = 0
c##################
c DR# 10552
c OLD	  if (iexp.gt.18 .or. iexp.lt.-3) lexp = iexp
c######ajb
c NEW
	  if (iexp.gt.14 .or. iexp.lt.-3) lexp = iexp
c##################
	elseif (form.eq.'VIS') then
	  lexp = 0
	  if (iexp.gt.7 .or. iexp.lt.-3) lexp = iexp
	  lexp = lexp/3*3
	  if (lexp.lt.0) lexp = lexp-3
	elseif (form.eq.'ENG') then
	  lexp = iexp/3*3
	  if (lexp.ne.iexp .and. iexp.lt.0) lexp = lexp-3
	elseif (form.eq.'SCI') then
	  lexp = iexp
	elseif (form.eq.'MAN') then
	  lexp = 0
	else
	  call M$ERROR ('Unknown format in M$D2A_FORM: '//form)
	endif
c
c adjust exponent
	kexp = lexp - iexp
	if (kexp.lt.0) then
	  tmp(3:3-kexp) = tmp(4:3-kexp)//'.'
	elseif (kexp.gt.0) then
	  tmp2 = tmp			! Sun cannot copy strings backwards
	  tmp(4+kexp:17) = tmp2(4:17)
	  tmp(3+kexp:3+kexp) = tmp2(2:2)
	  tmp(2:2+kexp) = '0.000000000000'
	endif
	if (lexp.eq.iexp) then
	  ! exponent already formed
	elseif (lexp.gt.0) then
	  write(tmp(19:21),'(1h+,i2.2)') lexp
	elseif (lexp.lt.0) then
	  write(tmp(19:21),'(1h-,i2.2)') -lexp
	endif
c
c remove redundant zeros
	i = 17
	do while (tmp(i-1:i).eq.'00')
	  i = i-1
	enddo
c
c remove trailing period
	if (tmp(i-1:i).eq.'.0') i = i-2
c
c remove leading space if positive number
	if (tmp(1:1).eq.' ') then
	  i = i-1
	  tmp(1:i) = tmp(2:i+1)
	endif
c
c put into final string
	if (lexp.eq.0) then
	  i = min(i,ls)
	  string = tmp(1:i)
	else
	  i = min(i+4,ls)
	  string = tmp(1:i-4)//tmp(18:21)
	endif
	M$D2A_FORM = i
	return
	end



	function M$F2A (fval,string)
c
	character*(*) string
	real*4 fval
	real*8 dval
c
	call M$F2D (fval,dval)
	M$F2A = M$D2A(dval,string)
	return
c
	end


	function M$L2A (ival,string)
c
	character*(*) string*(*),tmp*20
	integer*4 ival
c
	write(tmp,'(i20)') ival
	i = 20
	do while (tmp(i:i).ne.' ')
	  i = i-1
	enddo
	string = tmp(i+1:20)
	M$L2A = 20-i
	return
	end


	subroutine M$F2D (fval,dval)
c use only 1st seven digits of precision
	real*8 dval
	real*4 fval
	character*13 string
c
	write(string,'(e13.7)') fval
	read(string,*) dval
c
	return
	end


	subroutine M$D2F (dval,fval)
c nothing fancy here yet
	real*8 dval
	real*4 fval
	fval = dval
c
	return
	end



	function M$DMS2DEG (dms_str, deg)

	character*(*)	dms_str		!deg/min/sec string
	real*8		deg		!equivalent degrees
	real*8		min, sec	!intermediate min, sec
	integer*4	lt, lt1		!tic locators

	deg = 0
	lt = M$SEARCH (dms_str, '''')
	read (dms_str(1:lt-1), *, err=99) deg
	lt1 = M$SEARCH (dms_str(lt+1:), '''') + lt
	if (lt1.gt.lt) then
	  read (dms_str(lt+1:lt1-1), *, err=99) min	!ddd'mm'ss.sss format
	  read (dms_str(lt1+1:), *, err=99) sec
	  deg = SIGN(ABS(deg) + min/60 + sec/3600, deg)
	else
	  read (dms_str(lt1+1:), *, err=99) min		!ddd'mm.mmm format
	  deg = SIGN(ABS(deg) + min/60, deg)
	endif
	M$DMS2DEG = M$LENGTH(dms_str)
	return
   99	M$DMS2DEG = -1
	return
	end


	function M$DEG2DM (deg, dm_str)

	real*8		deg		!degrees
	character*(*)	dm_str		!deg/min string
	real*8		ddeg, dmin	!intermediate deg,min
	integer*4	ideg, i		!intermediate deg; temp index

	dm_str = ' 000''00.000'
c Get degrees, minutes, checking for overflow at 60
	ddeg = ABS(deg)
	ideg = INT(ddeg)
	dmin = DNINT(((ddeg - ideg)*60)*1000) * 0.001d0	!round to 3 decimals
	if (dmin .ge. 60) then
	  dmin = dmin - 60
	  ideg = ideg + 1
	endif
c Write to string
	write (dm_str(2:4), '(i3)', err=99) ideg
	write (dm_str(6:), '(f6.3)',err=99) dmin
c Zero-fill if needed
	if (dm_str(2:2) .eq. ' ') dm_str(2:2) = '0'
	if (dm_str(3:3) .eq. ' ') dm_str(3:3) = '0'
	if (dm_str(6:6) .eq. ' ') dm_str(6:6) = '0'
c Suppress trailing zeros in fraction
	i = 11
	do while (dm_str(i:i) .eq. '0')
	  i = i-1
	enddo
	if (i .eq. 8) i = i-1	!skip decimal if nothing follows
	M$DEG2DM = i
c Add sign, if needed
	if (deg .lt. 0) dm_str(1:1) = '-'
	return
   99	M$DEG2DM = -1
	return
	end


	function M$DEG2DMS (deg, dms_str)

	real*8		deg		!degrees
	character*(*)	dms_str		!deg/min/sec string
	real*8		ddeg,dmin,dsec	!intermediate deg,min,sec
	integer*4	ideg, imin, i	!intermediate deg, min; temp index

	dms_str = ' 000''00''00.000'
c Get degrees, minutes, seconds, checking for overflow at 60
	ddeg = ABS(deg)
	ideg = INT(ddeg)
	dmin = (ddeg - ideg)*60
	imin = INT(dmin)
	dsec = DNINT(((dmin - imin)*60)*1000) * 0.001d0	!round to 3 decimals
	if (dsec .ge. 60) then
	  dsec = dsec - 60
	  imin = imin + 1
	  if (imin .ge. 60) then
	    imin = imin - 60
	    ideg = ideg + 1
	  endif
	endif
c Write to string
	write (dms_str(2:4), '(i3)', err=99) ideg
	write (dms_str(6:7), '(i2)', err=99) imin
	write (dms_str(9:), '(f6.3)',err=99) dsec
c Zero-fill if needed
	if (dms_str(2:2) .eq. ' ') dms_str(2:2) = '0'
	if (dms_str(3:3) .eq. ' ') dms_str(3:3) = '0'
	if (dms_str(6:6) .eq. ' ') dms_str(6:6) = '0'
	if (dms_str(9:9) .eq. ' ') dms_str(9:9) = '0'
c Suppress trailing zeros in fraction
	i = 14
	do while (dms_str(i:i) .eq. '0')
	  i = i-1
	enddo
	if (i .eq. 11) i = i-1	!skip decimal if nothing follows
	M$DEG2DMS = i
c Add sign, if needed
	if (deg .lt. 0) dms_str(1:1) = '-'
	return
   99	M$DEG2DMS = -1
	return
	end



	function M$SPA (format)
c
	character*2 format
c
c determine scalars per atom from from mode
	M$SPA = M$NPT(format(1:1))

	if (M$SPA.eq.0) then
	  if (format.ne.'NH'.and.format.ne.'KW') call M$WARNING
     &	    ('Unknown format mode in M$SPA, format = '//format)
	  M$SPA = 1
	endif
c
	return
	end


	function M$BPS (format)
c
	character*2 format
c
c determine bytes per scalar from format type
	if (format.eq.'NH' .or. format.eq.'KW') then
	  M$BPS = 1
	elseif (format(1:1).eq.'U') then
	  M$BPS = ICHAR(format(2:2))		! in bytes
	  if (M$BPS.ge.128) M$BPS = M$BPS - 256	! in -bits
	else
	  M$BPS = M$NPT(format(2:2))
	  if (M$BPS.eq.0) then
	    call M$WARNING ('Illegal format type in M$BPS, format = '//format)
	    M$BPS = 1
	  endif
	endif
c
	return
	end


	function M$BPA (format)
c
	character*2 format
	M$BPA = M$SPA(format) * M$BPS(format)
c
	return
	end


	function M$BPE (hcb)
c
	include 'headers.inc'
	record /HEADER/ hcb
c
c bytes per element depends on type class
	if      (hcb.class.eq.1) then
	  M$BPE = M$BPA(hcb.format)
	else if (hcb.class.eq.2) then
	  M$BPE = hcb.subsize * M$BPA(hcb.format)
	else if (hcb.class.eq.3) then
	  M$BPE = hcb.record_length
	else if (hcb.class.eq.4) then
	  M$BPE = 8
	else if (hcb.class.eq.5) then
	  M$BPE = hcb.record_length
	else
	  call M$WARNING ('Unknown type class in M$BPE - treating as 1000')
	  M$BPE = M$BPA(hcb.format)
	endif
	if (M$BPE.eq.0) then
	  call M$WARNING('Encountered 0 bytes per element in M$BPE - set to 1')
	  M$BPE = 1
	endif
c
	return
	end



	function M$TYPE (char)
c
c typing a string depends on first character only
c
c Return	Means character is one of:
c   2		any character in the alphabet and the underscore
c   1		the digits 0 1 2 3 4 5 6 7 8 9
c   0		plus, minus, and period (i.e. numeric punctuation)
c  -2		special delimiters (comma, double quotes)
c  -1		all other characters (non-special delimiters)
c
c Valid results names start with a type 2 and thereafter contain types 1 or 2
c
	character*1 char
	byte	ascii_map(0:127) / 34*-1,-2,8*-1,0,-2,2*0,-1,10*1,7*-1,
     &				   26*2,4*-1,2,-1,26*2,5*-1 /
	save ascii_map
c
	M$TYPE = ascii_map( ICHAR(char) )
	return
c
	end


	function M$NPT (type)
c
c number of bytes per type is used for format modes as well
c
c  modes: S,C,V,Q,M,X,T,U =  1,2,3,4,9,10,16,1
c  modes: 0-9             =  0-9
c  modes: ICHAR(128-255)  =  1-128
c
c  types: P,B,O,I,L,F,D,A = -1,1,1,2,4,4,8,8
c
	character*1 type
	integer*4   itype
	byte	ascii_map(0:127) / 48*0,0,1,2,3,4,5,6,7,8,9,7*0,8,1,2,8,
     &		0,4,2*0,2,2*0,4,9,0,1,-1,4,0,1,16,1,3,0,10,39*0 /
	save ascii_map
c
	itype = ICHAR(type)
	if (itype.le.127) then
	  M$NPT = ascii_map( itype )
	else
	  M$NPT = itype - 127
	endif
	return
c
	end


	subroutine M$INDEX (hcb,abscissa,start,end)
c
	include 'headers.inc'
	record /HEADER/ hcb
c
	real*8 dstart,dend,xstart,xdelta,start,end
	integer*4 abscissa
c
	dstart = start
	dend   = end
	if (hcb.class .eq. 2) then
	  xstart = hcb.ystart
	  xdelta = hcb.ydelta
	else
	  xstart = hcb.xstart
	  xdelta = hcb.xdelta
	endif
c
	if (abscissa.ne.0) then
	  if (abscissa.eq.3) then
	    if (hcb.timecode.eq.0.d0) then
	      call M$WARNING (hcb.file_name//'has no time code field')
	    else
	      dstart = dstart - hcb.timecode
	      dend   = dend   - hcb.timecode
	    endif
	  endif
	  start = DNINT( (dstart-xstart)/xdelta ) + 1.d0
	  start = dmax1 (1.d0,start)
	  if (start .gt. hcb.size) call M$ERROR
     +	     ('Abscissa past end of file in M$INDEX')
	  end = DNINT( (dend-xstart)/xdelta ) + 1.d0
	  if (dend .eq. dstart) end = start
	  if (dend .lt. dstart) end = hcb.size
	  end = dmin1 (end,hcb.size)
c
	else
	  start = DNINT(dstart)
	  start = dmax1 (1.d0,start)
	  if (start .gt. hcb.size) call M$ERROR
     +	     ('Index past end of file in M$INDEX')
	  end = DNINT(dend)
	  if (dend .eq. dstart) end = start
	  if (dend .lt. dstart) end = hcb.size
	  end = dmin1 (end,hcb.size)
	endif
c
	return
	end


	subroutine M$ALIGN (hcb1,hcb2,abscissa,size,start)
c
	include 'headers.inc'
	record /HEADER/ hcb1,hcb2
c
	integer*4 abscissa
	real*8 size,start
	real*8 start1,delta1,off1,start2,delta2,off2
	real*8 tc1,tc1x,tc2,tc2x,tcoff
c
c do special 2000 changeover of x/y fields
	if (hcb1.class .eq. 2) then
	  start1 = hcb1.ystart
	  delta1 = hcb1.ydelta
	else
	  start1 = hcb1.xstart
	  delta1 = hcb1.xdelta
	endif
	if (hcb2.class .eq. 2) then
	  start2 = hcb2.ystart
	  delta2 = hcb2.ydelta
	else
	  start2 = hcb2.xstart
	  delta2 = hcb2.xdelta
	endif
c
c handle abscissa registration case
	if (abscissa.ne.0) then
	  if (delta1.eq.0.d0 .or. delta2.eq.0.d0) call M$ERROR
     &	  ('Cannot abscissa M$ALIGN with delta = 0')
	  if ( abs(delta2-delta1)/delta1 .gt. 1.d-6) call M$ERROR
     &	  ('Cannot abscissa M$ALIGN with delta mismatch > 1e-6')
	  if (abscissa.eq.3) then
	    if (hcb1.timecode.eq.0.d0) then
	      call M$WARNING (hcb1.file_name//'has no time code field')
	    elseif (hcb2.timecode.eq.0.d0) then
	      call M$WARNING (hcb2.file_name//'has no time code field')
	    else
	      call M$GET_EPOCH (hcb1,tc1,tc1x)
	      call M$GET_EPOCH (hcb2,tc2,tc2x)
	      tcoff = tc2x-tc1x
	      tcoff = (tc2-tc1)+tcoff
	      start2 = start2 + tcoff
	    endif
	  endif
	  if (start1.ge.start2) then
	    off1 = 0.d0
	    off2 = dnint ((start1-start2)/delta1)
	    start = start1
	  else
	    off2 = 0.d0
	    off1 = dnint ((start2-start1)/delta1)
	    start = start2
	  endif
	  start1 = start1 + off1*delta1
	  start2 = start2 + off2*delta2
	  if ( abs(start2-start1)/delta1 .gt. 1.d-3) call M$WARNING
     &	  ('Start abscissa mismatch > 1e-3 * delta in M$ALIGN')
	else
	  off1 = 0.d0
	  off2 = 0.d0
	  start = start1
	endif
c
	if (off1.ge.hcb1.size) call M$ERROR
     &    ('Index past end of file #1 in M$ALIGN')
	if (off2.ge.hcb2.size) call M$ERROR
     &    ('Index past end of file #2 in M$ALIGN')
c
c set internals header offset fields
	hcb1.offset = off1
	hcb2.offset = off2
c
c return mutual size (mutual start set above)
	size = min (hcb1.size-off1, hcb2.size-off2)
c
	return
	end


	subroutine M$TALIGN (hcb1,hcb2,astart,aoff,size)
c
	include 'headers.inc'
	record /HEADER/ hcb1,hcb2
c
	real*8 astart,aoff,start1,delta1,start2,delta2
	real*8 off1,off2,size,start
	character*1 address
c
	if (hcb1.class .ne. 1) call M$ERROR
     &	('Only type 1000 files supported in M$TALIGN')
	if (hcb2.class .ne. 1) call M$ERROR
     &	('Only type 1000 files supported in M$TALIGN')
c
	start1 = hcb1.xstart
	delta1 = hcb1.xdelta
	start2 = hcb2.xstart
	delta2 = hcb2.xdelta
c
	if (delta1.eq.0.d0) call M$ERROR
     &	('Cannot abscissa M$TALIGN with delta = 0')
	if ( abs(delta2-delta1)/delta1 .gt. 1.d-6) call M$ERROR
     &	('Cannot abscissa M$TALIGN with delta mismatch > 1e-6')
c
	start = astart
c
	noff = nint(aoff/delta1)
	off1 = dnint( (start-start1)/delta1 )
	if (off1.lt.0.d0) then
	  off1 = 0.d0
	  start = start1
	endif
	off2 = dnint( (start-start2)/delta2 ) + DBLE(noff)
	if (off2.lt.0.d0) then
	  off2 = 0.d0
	  start  = start2 - DBLE(noff)*delta1
	  off1 = dnint( (start-start1)/delta1 )
	endif
c
	start1 = start1 + off1*delta1
	start2 = start2 + off2*delta2
	if ( abs(start2-start1)/delta1 .gt. 1.d-3) call M$WARNING
     &	('Start abscissa mismatch > 1e-3 * delta in M$TALIGN')
c
	if (off1.ge.hcb1.size) call M$ERROR
     &    ('Index past end of file #1 in M$TALIGN')
	if (off2.ge.hcb2.size) call M$ERROR
     &    ('Index past end of file #2 in M$TALIGN')
c
c set internals header offset fields
	hcb1.offset = off1
	hcb2.offset = off2
c
c return mutual size
	size = min (hcb1.size-off1, hcb2.size-off2)
c
c return actual start
	astart = start1
c
c return actual offset if not desired zero (probably passed as 0.d0)
	if (aoff.ne.0.d0) aoff = noff*delta1
c
	return
	end


	function M$UNITS_NAME (units)
c
	include 'types.inc'
	integer*4 units
	character*14 resname
	character*35 aval
	real*8	     dval
c
	i = units
c
c branch on units code
	if (i.eq.0) then
	  M$UNITS_NAME = 'None (U)'
	else if (i.eq.1) then
	  M$UNITS_NAME = 'Time (sec)'
	else if (i.eq.2) then
	  M$UNITS_NAME = 'Delay (sec)'
	else if (i.eq.3) then
	  M$UNITS_NAME = 'Frequency (Hz)'
	else if (i.eq.4) then
	  M$UNITS_NAME = 'Time code format'
	else if (i.eq.5) then
	  M$UNITS_NAME = 'Distance (m)'
	else if (i.eq.6) then
	  M$UNITS_NAME = 'Speed (m/sec)'
	else if (i.eq.7) then
	  M$UNITS_NAME = 'Acceleration (m/sec^2)'
	else if (i.eq.8) then
	  M$UNITS_NAME = 'Jerk (m/sec^3)'
	else if (i.eq.9) then
	  M$UNITS_NAME = 'Doppler (Hz)'
	else if (i.eq.10) then
	  M$UNITS_NAME = 'Doppler rate (Hz/sec)'
	else if (i.eq.11) then
	  M$UNITS_NAME = 'Energy (J)'
	else if (i.eq.12) then
	  M$UNITS_NAME = 'Power (W)'
	else if (i.eq.13) then
	  M$UNITS_NAME = 'Mass (g)'
	else if (i.eq.14) then
	  M$UNITS_NAME = 'Volume (l)'
	else if (i.eq.15) then
	  M$UNITS_NAME = 'Angular power density (W/ster)'
	else if (i.eq.16) then
	  M$UNITS_NAME = 'Integrated power density (W/rad)'
	else if (i.eq.17) then
	  M$UNITS_NAME = 'Spatial power density (W/m^2)'
	else if (i.eq.18) then
	  M$UNITS_NAME = 'Integrated power density (W/m)'
	else if (i.eq.19) then
	  M$UNITS_NAME = 'Spectral power density (W/MHz)'
	else if (i.eq.20) then
	  M$UNITS_NAME = 'Amplitude (U)'
	else if (i.eq.21) then
	  M$UNITS_NAME = 'Real (U)'
	else if (i.eq.22) then
	  M$UNITS_NAME = 'Imaginary (U)'
	else if (i.eq.23) then
	  M$UNITS_NAME = 'Phase (rad)'
	else if (i.eq.24) then
	  M$UNITS_NAME = 'Phase (deg)'
	else if (i.eq.25) then
	  M$UNITS_NAME = 'Phase (cycles)'
	else if (i.eq.26) then
	  M$UNITS_NAME = '10*Log (U)'
	else if (i.eq.27) then
	  M$UNITS_NAME = '20*Log (U)'
	else if (i.eq.28) then
	  M$UNITS_NAME = 'Magnitude (U)'
	else if (i.eq.29) then
	  M$UNITS_NAME = 'Frequency difference (Hz)'
	else if (i.le.30) then
	  M$UNITS_NAME = 'Unknown (U)'
	else if (i.eq.31) then
	  M$UNITS_NAME = 'General dimensionless'
	else if (i.eq.32) then
	  M$UNITS_NAME = 'Counts'
	else if (i.eq.33) then
	  M$UNITS_NAME = 'Angle (rad)'
	else if (i.eq.34) then
	  M$UNITS_NAME = 'Angle (deg)'
	else if (i.eq.35) then
	  M$UNITS_NAME = 'Relative power (dB)'
	else if (i.eq.36) then
	  M$UNITS_NAME = 'Relative power (dBm)'
	else if (i.eq.37) then
	  M$UNITS_NAME = 'Relative power (dBW)'
	else if (i.eq.38) then
	  M$UNITS_NAME = 'Solid angle (ster)'
	else if (i.eq.40) then
	  M$UNITS_NAME = 'Distance (ft)'
	else if (i.eq.41) then
	  M$UNITS_NAME = 'Distance (nmi)'
	else if (i.eq.42) then
	  M$UNITS_NAME = 'Speed (ft/sec)'
	else if (i.eq.43) then
	  M$UNITS_NAME = 'Speed (nmi/sec)'
	else if (i.eq.44) then
	  M$UNITS_NAME = 'Speed (knots=nmi/hr)'
	else if (i.eq.45) then
	  M$UNITS_NAME = 'Acceleration (ft/sec^2)'
	else if (i.eq.46) then
	  M$UNITS_NAME = 'Acceleration (nmi/sec^2)'
	else if (i.eq.47) then
	  M$UNITS_NAME = 'Acceleration (knots/sec)'
	else if (i.eq.48) then
	  M$UNITS_NAME = 'Acceleration (G)'
	else if (i.eq.49) then
	  M$UNITS_NAME = 'Jerk (G/sec)'
	else if (i.eq.50) then
	  M$UNITS_NAME = 'Rotation (rps)'
	else if (i.eq.51) then
	  M$UNITS_NAME = 'Rotation (rpm)'
	else if (i.eq.52) then
	  M$UNITS_NAME = 'Angular velocity (rad/sec)'
	else if (i.eq.53) then
	  M$UNITS_NAME = 'Angular velocity (deg/sec)'
	else if (i.eq.54) then
	  M$UNITS_NAME = 'Angular acceleration (rad/sec^2)'
	else if (i.eq.55) then
	  M$UNITS_NAME = 'Angular acceleration (deg/sec^2)'
	else if (i.eq.60) then
	  M$UNITS_NAME = 'Latitude (deg)'
	else if (i.eq.61) then
	  M$UNITS_NAME = 'Longitude (deg)'
	else if (i.eq.62) then
	  M$UNITS_NAME = 'Altitude (ft)'
	else if (i.eq.63) then
	  M$UNITS_NAME = 'Altitude (m)'

	else if (i.ge.100 .and. i.le.127) then
	  write (resname,'(11hUNITS_NAME_,i3)') i
	  if (M$RFIND(resname,aval,dval).eq.'A') then
	    M$UNITS_NAME = aval
	  else
	    M$UNITS_NAME = 'User Undefined (U)'
	  endif
	else
	  M$UNITS_NAME = 'Unknown (U)'
	endif
c
	return
	end


	function M$FORMAT_NAME (format)
c
	include 'types.inc'
	character*2 format
	character*1 mode,rep
	character*10 cmode
	character*12 crep
	byte ubpa
	equivalence (rep,ubpa)
c
	mode = format(1:1)
	rep  = format(2:2)
c
c check for standard/special format types
	if      (format.eq.'SF') then
	  M$FORMAT_NAME = 'Scalar Float - (R*04)'
	else if (format.eq.'SI') then
	  M$FORMAT_NAME = 'Scalar Integer - (I*02)'
	else if (format.eq.'SL') then
	  M$FORMAT_NAME = 'Scalar Long Int - (I*04)'
	else if (format.eq.'SD') then
	  M$FORMAT_NAME = 'Scalar Double - (R*08)'
	else if (format.eq.'CF') then
	  M$FORMAT_NAME = 'Complex Float - (C*08)'
	else if (format.eq.'NH') then
	  M$FORMAT_NAME = 'NonHomogeneous'
	else if (format.eq.'KW') then
	  M$FORMAT_NAME = 'KeyWord Struct'
	else
c
c handle generic case mode
	  if      (mode.eq.'S') then
	    cmode = 'Scalar'
	  else if (mode.eq.'C') then
	    cmode = 'Complex'
	  else if (mode.eq.'V') then
	    cmode = 'Vector'
	  else if (mode.eq.'Q') then
	    cmode = 'Quad'
	  else if (mode.eq.'M') then
	    cmode = 'Matrix'
	  else if (mode.eq.'T') then
	    cmode = 'TransM'
	  else if (mode.eq.'U') then
	    cmode = 'UserDef'
	  else if (mode.eq.'X') then
	    cmode = '10 *'
	  else if (mode.ge.'1' .and. mode.le.'9') then
	    cmode = mode//' *'
	  else
	    cmode = 'Undefined'
	  endif
c
c handle generic case format representation
	  if (mode.eq.'U') then
	    if (ubpa.ge.0) write (crep,'(i3,3h by)') ubpa
	    if (ubpa.lt.0) write (crep,'(i3,4h bit)') -ubpa
	  else if (rep.eq.'F') then
	    crep = 'Float'
	  else if (rep.eq.'L') then
	    crep = 'Long Int'
	  else if (rep.eq.'I') then
	    crep = 'Integer'
	  else if (rep.eq.'B') then
	    crep = 'Byte'
	  else if (rep.eq.'O') then
	    crep = 'Offset Byte'
	  else if (rep.eq.'P') then
	    crep = 'Bit Pack'
	  else if (rep.eq.'D') then
	    crep = 'Double'
	  else if (rep.eq.'A') then
	    crep = 'Ascii'
	  else
	    crep = 'Undefined'
	  endif
c
c combine mode and rep
	  M$FORMAT_NAME = cmode(1:M$LENGTH(cmode)) // ' ' //
     +			crep(1:M$LENGTH(crep))
c
	endif
c
	return
	end



	function M$COMPTYPE_NAME (type)
c
	include 'types.inc'
	integer*4 type
c
	i = type
c
c branch on component type code
	if (i.eq.1) then
	  M$COMPTYPE_NAME = 'Scalar'
	else if (i.eq.2) then
	  M$COMPTYPE_NAME = 'Cartesian'
	else if (i.eq.3) then
	  M$COMPTYPE_NAME = 'Spherical'
	else if (i.eq.4) then
	  M$COMPTYPE_NAME = 'Cylindric'
	else if (i.eq.5) then
	  M$COMPTYPE_NAME = 'Ellipsoid'
	else if (i.eq.6) then
	  M$COMPTYPE_NAME = 'Geodetic'
	else if (i.eq.10) then
	  M$COMPTYPE_NAME = 'Matrix'
	else
	  M$COMPTYPE_NAME = 'Undefined'
	endif
c
	return
	end



	function M$PROMOTE_FORMAT (form1, form2, form3)
c
	implicit none
c
	include 'types.inc'
	character    form1*2, form2*2, form3*(*)
	character*1  mode, type, t1, t2, t3

	character*1	M$PROMOTE_TYPE3
c
c promote format mode
	t1 = form1(1:1)
	t2 = form2(1:1)
	t3 = form3(1:1)
	if (t1.eq.'V' .or. t2.eq.'V' .or. t3.eq.'V') then
	  mode = 'V'
	elseif (t1.eq.'C' .or. t2.eq.'C' .or. t3.eq.'C') then
	  mode = 'C'
	else
	  mode = 'S'
	endif
c
c promote format type
	t1 = form1(2:2)
	t2 = form2(2:2)
	if (form3(1:1) .eq. '(') then
	  type = M$PROMOTE_TYPE3(t1, t2, form3)
	else
	  t3 = form3(2:2)
	  type = M$PROMOTE_TYPE3(t1, t2, t3)
	endif
c
	M$PROMOTE_FORMAT = mode//type
c
	return
	END	! M$PROMOTE_FORMAT



	function M$PROMOTE_TYPE (t1, t2)
c
	implicit none
c
	include 'types.inc'
	character*1  t1, t2

	character*1	M$PROMOTE_TYPE3
c
c promote format type

	M$PROMOTE_TYPE = M$PROMOTE_TYPE3(t1, t2, 'B')

	return
	END	! M$PROMOTE_TYPE



	character*1 function M$PROMOTE_TYPE3 (t1, t2, t3)
	implicit none
	include 'types.inc'
	character  t1, t2, t3*(*)

	character  pfostr*(L_max_ascii), type
	integer*4  stat

	integer*4  M$PROMOTE_FORMAT_OVERRIDE

	if (t1.eq.'A' .OR. t2.eq.'A' .OR. t3.eq.'A')
     +	  call M$ERROR('Non-numeric format types cannot be promoted')

	if (M$GET_PSWITCH('DP')) then
	  M$PROMOTE_TYPE3 = 'D'
	  return
	endif

	if (M$GET_USWITCH('PFO', pfostr) .gt. 0) then
	  stat = M$PROMOTE_FORMAT_OVERRIDE(t1, t2, pfostr, type)
	  if (stat .lt. 0) then
	    if (stat .eq. -1) then
	      call M$ERROR('Illegal value in /PFO switch: '//pfostr)
	    else
	      call M$ERROR('Non-numeric or invalid type in /PFO switch: '
     +			   //type)
	    endif
	  endif
	  if (type .ne. ' ') then
	    M$PROMOTE_TYPE3 = type
	    return
	  else
	    type = 'B'	! 'B' is harmless: anything higher will override
	  endif
	endif

	if (t3(1:1) .eq. '(') then
	  stat = M$PROMOTE_FORMAT_OVERRIDE(t1, t2, t3, type)
	  if (stat .lt. 0) then
	    if (stat .eq. -1) then
	      call M$ERROR('Illegal value in 3rd argument to '
     +			   //'M$PROMOTE_FORMAT: '//t3)
	    else
	      call M$ERROR('Non-numeric or invalid type in 3rd arg to '
     +			   //'M$PROMOTE_FORMAT: '//type)
	    endif
	  endif
	  if (type .ne. ' ') then
	    M$PROMOTE_TYPE3 = type
	    return
	  else
	    type = 'B'	! 'B' is harmless: anything higher will override
	  endif
	endif

	! Determine type the old fashioned way: 'DFLIB'
	if (t1.eq.'D' .or. t2.eq.'D' .or. t3.eq.'D') then
	  type = 'D'
	elseif (t1.eq.'F' .or. t2.eq.'F' .or. t3.eq.'F') then
	  type = 'F'
	elseif (t1.eq.'L' .or. t2.eq.'L' .or. t3.eq.'L') then
	  type = 'L'
	elseif (t1.eq.'I' .or. t2.eq.'I' .or. t3.eq.'I') then
	  type = 'I'
	else
	  type = 'B'
	endif

	M$PROMOTE_TYPE3 = type
	return
	END	! M$PROMOTE_TYPE3



C M$PROMOTE_FORMAT_OVERRIDE
C
C Parses and interprets the value of the /PFO switch or the 3rd
C argument to M$PROMOTE_FORMAT
C
	integer*4 function M$PROMOTE_FORMAT_OVERRIDE (t1, t2, pfostr, type)
	character	t1, t2, pfostr*(*), type

	include 'types.inc'
	integer*4   returnstat, pfolen, len, ii, jj, nn
	character   str*16

	returnstat = 0
	type = ' '

	pfolen = M$LENGTH(pfostr)
	ii = 1
	jj = pfolen
	if (pfostr(1:1) .eq. '(') then
	  ii = 2
	  jj = pfolen - 1
	endif

	nn = 1
	do while (.TRUE.)
	  len = ABS(M$PARSE(pfostr(ii:jj), nn, str))
	  if (len .eq. 0) then
	    goto 50	! all done - no more entries
	  elseif (len .eq. 1) then
	    type = str(1:1)
	    if (M$SEARCH('DFLIB', type) .eq. 0)
     +	      returnstat = -2	! Bad type
	    goto 50
	  elseif (len .eq. 3) then
	    if (((t1.eq.str(1:1) .OR. str(1:1) .eq. '#') .AND.
     +		 (t2.eq.str(2:2) .OR. str(2:2) .eq. '#'))
     +				 .OR.
     +	        ((t2.eq.str(1:1) .OR. str(1:1) .eq. '#') .AND.
     +		 (t1.eq.str(2:2) .OR. str(2:2) .eq. '#')))
     +	    then
	      type = str(3:3)
	      if (M$SEARCH('DFLIB', type) .eq. 0)
     +		returnstat = -2	! Bad type
	      goto 50
	    endif
	  else
	    ! Illegal length
	    returnstat = -1
	    goto 50
	  endif
	  nn = nn + 1
	enddo

  50	M$PROMOTE_FORMAT_OVERRIDE = returnstat
	return
	END	! M$PROMOTE_FORMAT_OVERRIDE



c M$NEXTLINE
c
c Description:
c   Reads the next relevant line from the given file unit into the given
c string, up to the given maximum number of characters.
c   Blank lines and comment lines are ignored.
c   In-line comments are removed.
c   Lines ending with a continuation character ('&') are concatenated up
c to the maximum length, except when there are two final cont. chars ('&&'),
c in which case they are treated as a single ampersand.
c
c Arguments:
c   f_in (IN)  --   logical unit of the file to read from
c   string (OUT) -- string to hold the next relevant line
c   verbose (IN) -- causes each input line to print if TRUE; stays
c		    quiet during reading if FALSE
c
c Returns:
c   1 if everything went OK.
c  -1 if the end of file was found.
c  -2 if a fatal I/O error was encountered.
c
c Revisions:
c   DTJ 3/10/94 - the logical 'contchar' needs to be reinitialized to .FALSE.
c		  for each time through the loop

	function M$NEXTLINE(f_in, string, verbose)

	implicit none

	include 'types.inc'
c Arguments
	integer*4	f_in, srcln
	character*(*)	string
	logical*4	verbose
c Locals
	logical*4	countsrc
	integer*4	maxchars, strlen
	character*128	line
	logical*4	inquote, contchar
	integer*4	i, j, iscan, lline
c
c Code
c
	countsrc = .FALSE.
	goto 5		! don't run over the entry statement

	entry	M$NEXTLINE_SOURCE(f_in, string, verbose, srcln)

	countsrc = .TRUE.

c Initializations
5	maxchars = LEN(string)
	string   = ' '
	strlen   = 1
	inquote  = .FALSE.
	if (countsrc) srcln    = 0

20	read (f_in,'(a)',err=90,end=80) line	! top of the "loop"
	if (countsrc) srcln = srcln + 1

c Remove trailing blanks; check for blank lines
	j = M$LENGTH(line)
	if (j .eq. 0) then
	  string(strlen:strlen) = ' '
	  goto 10	! just leave with what <string> has so far
	endif

c Show the user what line was just read
	if (verbose) write(*,*) line(1:j)

c Scan line for blanks, in-line comments, etc.
c Leading blanks
	i = 1
	do while (line(i:i) .le. ' ')
	  i = i + 1
	enddo

c Comment line
	if (line(i:i) .eq. '!') then
	  string(strlen:strlen) = ' '
	  goto 10 	! just leave with <string> so far
	endif

c In-line comments
	iscan = i
	do while (iscan .le. j)

	  if      (line(iscan:iscan) .eq. '"') then
	    inquote = .not. inquote
	  else if (line(iscan:iscan) .eq. '!') then
	    if (.not. inquote) j = iscan - 1
	  endif

	  iscan = iscan + 1

	enddo

c Trailing blanks
	do while (line(j:j) .le. ' ')
	  j = j - 1
	enddo

	lline = j - i + 1
	string(strlen:) = line(i:j)
	strlen = strlen + lline - 1

c Check for the double ampersand
c DTJ 3/10/94 - initialize contchar for each time thru loop
	contchar = .FALSE.
	if (string(strlen-1:strlen) .eq. '&&') then
	  string(strlen:strlen) = ' '
	  strlen = strlen - 1
	else
	  contchar = string(strlen:strlen) .eq. '&'
	endif

c Getting ready to return - do we have anything?
10	if (contchar .or.
     &	    (string .eq. ' ').or.(strlen .le. 0)) goto 20  ! back for more

c check for string overflow
	if (strlen .gt. maxchars) then
	  write(line,'(A, i4.0, A)')
     &		'Line too long; exceeded maximum length of ', maxchars,
     &		' characters'
	  call M$ERROR(line)
	endif

	M$NEXTLINE        = 1          ! normal exit
	M$NEXTLINE_SOURCE = 1
	return

 80     M$NEXTLINE        = -1         ! end-of-file
	M$NEXTLINE_SOURCE = -1
	return

 90     M$NEXTLINE        = -2         ! fatal I/O error
	M$NEXTLINE_SOURCE = -1
	return

	end



	subroutine M$EXPAND_IOTYPES(types, ntypes)

	implicit none
	include 'types.inc'
	character*(*)	types
	character*80	newtypes
	integer*4	ntypes, i, j, k, rep

	j = 1
	k = 0
	do while (j .le. ntypes)

	  i = j
	  do while (M$TYPE(types(j:j)) .eq. 1 .and. j .lt. ntypes)
	    j = j + 1
	  enddo
	  if (i .ne. j) then
	    read(types(i:j-1), *) rep
	  else
	    rep = 1
	  endif

	  do i = 1, rep
	    k = k + 1
	    newtypes(k:k) = types(j:j)
	  enddo

	  j = j + 1

	enddo

	ntypes = k
	types = newtypes

	return
	end



	subroutine M$PUT_EPOCH (hcb,time1,time2,toQuads)
c
	include 'headers.inc'
	record /HEADER/ hcb
	real*8 timeu,time1,time2,timep
	character*32 psec
	logical*4 toQuads
c
	timeu = DNINT( (time2)*1d6 ) * 1d-6
	timep = DNINT( (time2-timeu)*1d12 ) * 1d-12
	hcb.timecode = time1 + timeu
	if (timep.ne.0) then
	  lp = M$D2A(timep,psec)
	  mode = MK_main + MK_replace
	  i = M$PUT_KEYDATA (hcb,'TC_PREC',%val(%loc(psec)),lp,mode,'A')
	  if (i.le.0) call M$WARNING
     &	  ('No room in header for extra timecode precision')
	else
	  i = M$DEL_KEYDATA (hcb,'TC_PREC',MK_main)
	endif
c
c also save equivalent timecode in Type 5000 Epoch quadwords
	if (toQuads) call M$EPOCH_HOUR_ANGLE (time1+time2,
     &		hcb.quad(10), hcb.quad(11), hcb.quad(12))
	return
	end


	subroutine M$GET_EPOCH (hcb,time1,time2)
c
	include 'headers.inc'
	record /HEADER/ hcb
	real*8 timep,time1,time2
	character*32 psec
c
	time1 = DNINT(hcb.timecode)
	time2 = DNINT((hcb.timecode-time1)*1d6)*1d-6
	i = M$GET_KEYDATA (hcb,'TC_PREC',%val(%loc(psec)),32,MK_main,'A')
	if (i.gt.0) then
	  if (M$A2D(psec(1:i),timep).gt.0) time2 = time2+timep
	endif
c
	return
	end




c ***********************************************************************
c			Time Utilities
c ***********************************************************************

	function M$TIME_FORMAT (str, type, strout, sec, frac, warn)
	implicit none
	include 'types.inc'
c Arguments:
	character*(*)	str, strout
	character*1	type
	real*8		sec, frac
	logical*4	warn
c Locals:
	integer*4	i,j, ls, year, month, day
	real*8		time
	character*40	tod_str, expect
	character*2	anum(31) /
     &	'01', '02', '03', '04', '05', '06', '07', '08', '09', '10',
     &	'11', '12', '13', '14', '15', '16', '17', '18', '19', '20',
     &	'21', '22', '23', '24', '25', '26', '27', '28', '29', '30',
     &	'31'/
	character*3	amonth(12) /
     &		'JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN',
     &		'JUL', 'AUG', 'SEP', 'OCT', 'NOV', 'DEC'/

	M$TIME_FORMAT = -1
	sec = 0
	frac = 0

c create TOD string from input
	if (type .eq. ' ' .or. type .eq. '_') then	! standard format
	  expect = '[yy]yy:mm:dd::hh:mm:ss'
	  tod_str = str

	elseif (type .eq. 'V') then		! VAX format
	  expect = 'dd-MMM-yy[yy]:[hh:mm:ss]'
	  j = INDEX (str, ':')  !handle 2-char or 4-char year
	  if (j .gt. 0) then
	    tod_str = str(8:j-1)//':00:'//str(1:2)//'::'//str(j+1:)
	  else
	    j = M$LENGTH (str)
	    tod_str = str(8:j-1)//':00:'//str(1:2)//'::0'
	  endif
	  call M$UPPERCASE (str)
	  do i = 1,12
	    if (str(4:6) .eq. amonth(i)) tod_str(j-6:j-5) = anum(i)
	  enddo

	elseif (type .eq. 'A') then		! ACQ format
	  expect = '[yy.ddd]:[hh:mm:ss]'
	  j = INDEX (str, '.')		!yy.ddd optional and may be yyyy.ddd
	  if (j .gt. 0) then
	    if (M$A2D (str(1:j+3), time) .le. 0) goto 99
	    year = INT (time)
	    day = NINT ((time-year)*1000)
	    call M$YDY2YMD (year, day, month, day)
	    tod_str = str(1:j-1)//':'//anum(month)//':'//anum(day)//'::'
	    ls = j+7
	  else
	    tod_str = '0::'
	    ls = 3
	  endif
	  j = INDEX (str, ':')		!hh:mm:ss optional
	  if (j .gt. 0) then
	    tod_str(ls+1:) = str(j+1:)
	  else
	    tod_str(ls+1:) = '0'
	  endif

	elseif (type .eq. 'T') then		! Timecode reader format
	  expect = 'ddd:hh:mm:ss'	!strictly positional
#if _VMS
	  call LIB$DATE_TIME (tod_str)         ! gives dd-mmm-yyyy hh:mm:ss.hh
	  tod_str(8:9) = tod_str(10:11)        ! gives dd-mmm-yy
	  tod_str(10:) = ' '
#else
	  call DATE (tod_str)		!use this year for year, dd-mmm-yy
#endif
	  tod_str(1:2) = tod_str(8:9)
	  if (M$A2D (tod_str(1:2)//str(1:3), time) .le. 0) goto 99
	  year = time/1000
	  day = time - year*1000
	  call M$YDY2YMD (year, day, month, day)
	  tod_str(3:) = ':'//anum(month)//':'//anum(day)//'::'//str(5:)

	elseif (type .eq. 'N') then		! NORAD format
	  expect = 'yyddd.ssssss'	!strictly positional
	  if (M$A2D (str(1:5), time) .le. 0) goto 99
	  year = time/1000
	  day = time - year*1000
	  call M$YDY2YMD (year, day, month, day)
	  tod_str = str(1:2)//':'//anum(month)//':'//anum(day)//'::'
	  if (str(6:6) .ne. '.') goto 99
	  if (M$A2D (str(6:), time) .le. 0) goto 99
	  if (M$D2A (time*86400, tod_str(11:)) .le. 0) goto 99

	elseif (type .eq. 'E') then 		! Epoch quadword format
	  expect = '[yy]yy:sec_of_year'
	  j = INDEX (str, ':')
	  if (j .le. 0) goto 99
	  if (M$A2D (str(1:j-1), time) .le. 0) goto 99
	  year = time
	  if (M$A2D (str(j+1:), time) .le. 0) goto 99
	  day = time / 86400 + 1
	  time = time - (day-1)*86400
	  call M$YDY2YMD (year, day, month, day)
	  tod_str = str(1:j-1)//':'//anum(month)//':'//anum(day)//'::'
	  if (M$D2A (time, tod_str(j+8:)) .le. 0) goto 99

	else
	  if (warn) call M$WARNING (type//
     &	    ' format not supported in M$TIME_FORMAT')
	  return
	endif

c get times from TOD string
	ls = M$LENGTH (tod_str)
	ls = M$TOD2TIMES (tod_str(1:ls), sec, frac)
	if (ls .le. 0) goto 99
	ls = M$TIMES2TOD (sec, frac, tod_str)
	strout = tod_str
	M$TIME_FORMAT = ls
	return

   99	if (.not. warn) return
	ls = M$LENGTH (expect)
	call M$WARNING (type//' time expect '//expect(1:ls)//
     &	  ', got '//str)
	return
	end


	subroutine M$SECNDS (ref_time, elapsed)
	implicit none
	include 'types.inc'
	real*4 ref_time, elapsed
	elapsed = SECNDS(ref_time)
	return
	end


	function M$TOD2SEC (tod_str, sec)
c
c Time-of-day format is [[yy]yy:mm:dd::][h]h:mm[:ss[.ssssss]],
c                    or [[yy]yy:mm:dd::]sssssss...
c Year assumed between 1950 and 2049 if first two chars suppressed
c
	implicit none
	include 'types.inc'
	character*(*)	tod_str		! time-of-day string
	real*8		sec, frac	! seconds, fractional seconds
	integer*4	mode		! 0 = single output, 1 = two outputs
	integer*4	ih, im		! intermediate hrs, min
	integer*4	ny, nm, nd, mjd	! year month day
	integer*4	lc, len, i	! colon locator, length

#if !_PFC
	character*80 fmthack
#endif

	mode = 0
	goto 1
	entry M$TOD2TIMES (tod_str, sec, frac)
	mode = 1

c Handle date
    1	len = M$LENGTH(tod_str)
	lc  = INDEX(tod_str,'::')
	if (lc.gt.0) then
	  if (INDEX(tod_str(1:lc-1),':') .eq.0) then
	    read(tod_str(1:lc-1), *, err=99,end=99) mjd
	  else
	    read(tod_str(lc-2:lc-1), '(i2)', err=99,end=99) nd		! day
	    read(tod_str(lc-5:lc-4), '(i2)', err=99,end=99) nm		! month
#if _PFC
	    read(tod_str(1   :lc-7), 100, err=99,end=99) ny		! year
100	    format (i<lc-7>)
#else
	    write(fmthack, "('(I',I12,')')") (lc-7)
	    read(tod_str(1   :lc-7), fmthack, err=99) ny		! year
#endif
	    if (ny.lt.100.and.ny.ge.50) ny = ny+1900
	    if (ny.lt.100.and.ny.lt.50) ny = ny+2000
	    call M$YMD2MJD (ny,nm,nd,mjd)
	  endif
	  i = lc+1
	else
	  mjd = 0
	  i = 0
	endif

c Handle time
c Must use explicit format, else expressions like 'T1' read as valid numbers!
	lc  = INDEX(tod_str(i+1:),':') + i
	if (lc.eq.i) then	! only seconds follow date
	  ih = 0
	  im = 0
	  sec = 0
	  if (len.gt.i) then
#if _PFC
             read (tod_str(i+1:len), 200, err=99,end=99) sec	!seconds
200          format (f<len-i>.0)
#else
	     write(fmthack,"('(F',I12,'.0)')") (len-i)
	     read (tod_str(i+1:len), fmthack, err=99) sec	!seconds
#endif
	  endif
	else			! hh:mm:ss follow date
#if _PFC
          read (tod_str(i+1:lc-1),  300, err=99,end=99) ih	!hours
300       format (i<lc-1-i>)
#else
	  write(fmthack, "('(I',I12,')')") (lc-1-i)
	  read (tod_str(i+1:lc-1),  fmthack, err=99) ih 	!hours
#endif
	  read (tod_str(lc+1:lc+2), '(i2)', err=99,end=99) im	!minutes
	  if (im.lt.0 .or. im.gt.59) goto 99
	  if (len .gt. lc+2) then
	    if (tod_str(lc+3:lc+3) .ne. ':') goto 99
#if _PFC
            read (tod_str(lc+4:len), 400, err=99,end=99) sec	!seconds
400         format (f<len-lc-3>.0)
#else
	    write(fmthack, "('(F',I12,'.0)')") (len-lc-3)
	    read (tod_str(lc+4:len), fmthack, err=99) sec	!seconds
#endif

	    if (sec.lt.0 .or. sec.ge.60) goto 99
	    if (mode .eq. 0) then
	      sec = DNINT(sec*1d6)*1d-6		!round to usec
	    else
	      sec = DNINT(sec*1d12)*1d-12	!round to psec
	    endif
	  else
	    sec = 0
	  endif
	endif

c Set outputs
	if (mode .eq. 0) then
	  sec = mjd*86400d0 + ih*3600 + im*60 + sec
 	  M$TOD2SEC = len
	else
	  frac = DMOD (sec, 1d0)
	  sec = mjd*86400d0 + ih*3600 + im*60 + DINT (sec)
	  M$TOD2TIMES = len
	endif
	return

   99	M$TOD2SEC = -1
	M$TOD2TIMES = -1
	return
	end


	function M$SEC2TOD (inputsec, out_str)
c
c Time-of-day format is [[yy]yy:mm:dd::][h]h:mm[:ss[.ssssss]],
c                    or [[yy]yy:mm:dd::]sssssss...
c Year assumed between 1950 and 2049 if first two chars suppressed
c
	implicit none
	include  'types.inc'
	real*8		inputsec, inputfrac ! what the user gave
	real*8		sec, frac	  ! seconds, fractional seconds
	character*(*)	out_str		  ! time-of-day string output
	character*33	tod_str		  ! time-of-day string
	real*8		s, dis, psec	  ! intermediate sec, psec
	integer*4	i, ih, im, is, ip ! intermediate hrs, mins, sec, psec
	integer*4	ny, nm, nd, mjd	  ! year, month, day, J2000 date

	sec  = dabs(inputsec)

	s = DNINT (sec*1d6) * 1d-6	! round to nearest usec
	ip = 0
	goto 1

	entry M$TIMES2TOD (inputsec, inputfrac, out_str)

	sec  = dabs(inputsec)
	frac = inputfrac

	psec = DNINT (frac*1d12)	! round to nearest psec
	if (psec .lt. 0) sec = sec - 1
	if (psec .lt. 0) psec = psec + 1d12
	ip = DMOD (psec, 1d6)
	s = sec + (psec-ip)*1d-12

c Handle date
    1	if (s.gt.86400 .or. inputsec.lt.0.d0) then
	  mjd = DINT (s/86400)
	  if (inputsec.lt.0.d0) then	! convert neg times to -day:+hh:mm:ss
	    dis = mjd*86400d0
	    if (s .ne. dis) then
	      dis = dis + 86400d0	! perform FLOOR operation on negative
	      mjd = mjd + 1		! times that aren't on day boundaries
	    endif
	    s = dis - s		! convert: e.g, -1 sec == -1::23:59:59
	    mjd = - mjd		! i.e. 1 sec before midnight yesterday
	  else
	    s = s - mjd*86400d0
	  endif
	  if (iabs(mjd).lt.365) then
	    i = M$L2A(mjd,tod_str(1:4)) + 2
	    tod_str(i-1:i) = '::'
	  else
	    call M$MJD2YMD (mjd, ny, nm, nd)
	    i = 12
	    write (tod_str(1:i),'(i4,1h:,i2.2,1h:,i2.2,2h::)',err=99) ny,nm,nd
	  endif
	else
	  i = 0
	endif

c Handle time
	im = s / 60
	s  = s - im*60
	ih = im / 60
	im = im - ih*60
	is = s
	write (tod_str(i+1:i+8),'(i2.2,1h:,i2.2,1h:,i2.2)',err=99) ih,im,is
	i = i+8
	if (ip .eq. 0) then		! fractional seconds to usec
	  is = NINT ((s-is)*1d6)
	  if (is.ne.0) then
	    write (tod_str(i+1:i+7),'(1h.,i6.6)',err=99) is
	    i = i+7
	  endif
	else				! fractional seconds to psec
	  is = psec*1d-6
	  write (tod_str(i+1:i+13),'(1h.,i6.6,i6.6)',err=99) is, ip
	  i = i+13
	endif
	M$SEC2TOD = i
	M$TIMES2TOD = i
	out_str = tod_str(1:i)
	return

   99	M$SEC2TOD = -1
	M$TIMES2TOD = -1
	out_str = ' '
	return
	end



	function M$TIMES2STR (sec, frac, type, strout)

	implicit none

	include 'types.inc'

c Arguments
	real*8		sec, frac
	character*1	type
	character*(*)	strout
c Locals
	integer*4	year, month, day
	character*(L_max_ascii)	tod_str, str
	integer*4	i, j, ls, lstr, mjd
	integer*4	doy
	integer*4	yoc		! year-of-century
	integer*4	sid		! seconds-in-day
	integer*4	siy		! seconds in the year
	real*8		time, fracday
	character*2	anum(31) /
     &	'01', '02', '03', '04', '05', '06', '07', '08', '09', '10',
     &	'11', '12', '13', '14', '15', '16', '17', '18', '19', '20',
     &	'21', '22', '23', '24', '25', '26', '27', '28', '29', '30',
     &	'31'/
	character*3	amonth(12) /
     &		'JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN',
     &		'JUL', 'AUG', 'SEP', 'OCT', 'NOV', 'DEC'/

	M$TIMES2STR = -1

	tod_str = ' '
	if (type .eq. ' ' .or. type .eq. '_'. or. type .eq. 'X') then
!
! [yy]yy:mm:dd::hh:mm:ss
!
	  ls = M$TIMES2TOD(sec, frac, tod_str)

	else

c Need to do this crazy SEARCHB thing because there might be trailing
c fractional seconds.  Otherwise we could count on the time-of-day being
c 8 characters, always.
	  lstr = M$TIMES2TOD(sec, frac, str)
	  i = M$SEARCHB(str(1:lstr), '::')
	  if (i .eq. 0) then
	    i = 1
	  else
	    i = i + 2
	  endif
	  str = str(i:lstr)//' '
	  lstr = lstr - i + 1

	  if (type .eq. 'V') then		! VAX format
!
! dd-MMM-yy[yy]:[hh:mm:ss]
!
	    call M$MJD2YMD(INT(sec/86400), year, month, day)
	    write(tod_str,
     &		'(i2.2,1h-,a3,1h-,i4.4,1h:,a)')
     &		day, amonth(month), year, str(1:lstr)

	  else		! 'T', 'A', 'N', or 'E' : all need the day of year

	    mjd = DINT (sec / 86400)
	    call M$MJD2YDY(mjd, year, doy)
	    yoc = MOD  (year, 100)
	    sid = DMOD (sec, 86400.0d0)

	    if (type .eq. 'T') then		! Timecode reader format
!
! ddd:hh:mm:ss
!
	      write (tod_str, '(i3.3,1h:,a)') doy, str(1:lstr)

	    elseif (type .eq. 'A') then		! ACQ format
!
! [yy.ddd]:[hh:mm:ss]
!
	      write (tod_str, '(i2.2,1h.,i3.3,1h:,a)')
     &		yoc, doy, str(1:lstr)


	    elseif (type .eq. 'N') then		! NORAD format
!
! yyddd.<fractional day>
!
c We cannot use M$D2A_FORM here because NORAD format may need the full 16
c digits of fractional precision in order to retain microsecond precision
c on the time.  We can get away with a lot of assumptions because we know
c what we're formatting will be a double precision number between 0 and 1.
	      write(str, '(f17.16)') (sid + frac)/86400.0d0
	      i = 17
	      do while (str(i:i) .eq. '0')
		i = i - 1
	      enddo
	      if (str(i:i) .eq. '.') i = i + 1
	      write (tod_str, '(i2.2,i3.3,a)')
     &				yoc, doy, str(1:i)

	    elseif (type .eq. 'E') then 		! Epoch quadword format
!
! [yy]yy:sec_of_year
!
!	      lstr = M$D2A((doy - 1) * 86400 + sid, str)
	      lstr = M$D2A_FORM(DMOD(sec, 86400.0d0) +
     &				     (doy - 1) * 86400.0d0 + frac,
     &		     str, 'MAN')
	      write (tod_str, '(i4.4,1h:,a)') year, str(1:lstr)

	    else
	      call M$WARNING (type//
     &	    ' format not supported in M$TIMES2STR')
	      return
	    endif	    ! another format
	  endif		  ! VAX format
	endif		! std. format

	ls = M$LENGTH(tod_str)
	strout = tod_str(1:ls)
	M$TIMES2STR = ls
	return

	end



	subroutine M$MJDCONVERSIONS_ ()
	! J1950 = Jan 1 1950 Midnight
	implicit none
	integer*4 ny, nm, nd, mjd, ndy, ndoy, ly, ndm(24)
	data ndm /0,31,59,90,120,151,181,212,243,273,304,334,
     &	          0,31,60,91,121,152,182,213,244,274,305,335/

	entry M$YMD2MJD (ny,nm,nd,mjd)
c Converts Year Month Day to modified Julian day J1950
	ly = 0
	if (MOD(ny,4).eq.0) ly = 12
	ndy = nd + ndm(nm+ly)
	mjd = (ny-1950)*365 + (ny+3)/4 - 488 + (ndy-1)
	return

	entry M$MJD2YMD (mjd,ny,nm,nd)
c Converts modified Julian day J1950 to Year Month Day
	ny = (mjd+712237.5d0)/365.25d0 + .002
	ndy = mjd - ( (ny-1950)*365 + (ny+3)/4 - 488 ) + 1
	ly = 0
	if (MOD(ny,4).eq.0) ly = 12
	do nm = 12,1,-1
	  if (ndy.gt.ndm(nm+ly)) then
	    nd = ndy - ndm(nm+ly)
	    return
	  endif
	enddo
	return

c <<< DTJ 3/23/93 -- 3rd argument was not used, causing incorrect results
	entry M$MJD2YDY (mjd,ny,ndoy)
c Converts modified Julian day J1950 to Year and DayofYear
	ny   = (mjd+712237.5d0)/365.25d0 + .002
	ndoy = mjd - ( (ny-1950)*365 + (ny+3)/4 - 488 ) + 1
	return

	entry M$YDY2YMD (ny,ndoy,nm,nd)
c Converts Year and DayofYear to Year Month Day
	ly = 0
	if (MOD(ny,4).eq.0) ly = 12
	do nm = 12,1,-1
	  if (ndoy.gt.ndm(nm+ly)) then
	    nd = ndoy - ndm(nm+ly)
	    return
	  endif
	enddo
	call M$WARNING ('Day of year <= 0 in M$YDY2YMD')
	nm = 1
	nd = ndoy
	return

	end

	function M$SEC2TSPEC (dval, string, mode)
	implicit none
	include 'types.inc'
	real*8 dval, dvalmod
	character*(*) string
	integer*4 mode
	real*8 sec_per_day
	data sec_per_day/86400.D0/

	if (dval.ge.0.d0 .and. dval.le.sec_per_day) then
	  M$SEC2TSPEC = M$SEC2TOD(dval,string)
	else
	  dvalmod = DMOD(dval,sec_per_day)
	  if (mode.eq.0 .and. dvalmod.le.0.d0) then
	    M$SEC2TSPEC = M$SEC2TOD(dvalmod+sec_per_day,string)
	  elseif (dvalmod .lt. 0.d0) then
	    string(1:1) = '-'
	    M$SEC2TSPEC = M$SEC2TOD(-dvalmod,string(2:)) + 1
	  else
	    M$SEC2TSPEC = M$SEC2TOD(dvalmod,string)
	  endif
	endif

	return
	end



