c************
c  Description:	XMIDAS record oriented utility routines (3000/5000/6000)
c
c  Author:
c
c  $Id: type5000.for,v 1.5 2014/11/05 21:14:15 jtm Exp $
c  $Revision: 1.5 $
c  $Modtime:   11 Aug 1997 06:35:34  $
c              DR 536469 jlr 10Nov2005 - In M$ADDT6SUBR properly account for
c                                        multiple atoms when updating record
c                                        length
c
c************

	function M$SUB_FIND (hcb,name)
	implicit none
	include 'headers.inc'
	include 'types.inc'
	record /HEADER/ hcb
	character*(*) name
	integer*4 n

	M$SUB_FIND = 0

	do n = 1,hcb.subrecords
	  if (M$SUB_NAME(hcb, n) .eq. name) then
	    M$SUB_FIND = n
	    return
	  endif
	enddo

	if (M$TYPE(name).eq.1) then
	  read(name(1:M$LENGTH(name)),*) n
	  if (n.gt.0 .and. n.le.hcb.subrecords) M$SUB_FIND = n
	endif
	return
	end


c	character*24 function M$SUB_NAME (hcb,n)
	function M$SUB_NAME (hcb,n)
	implicit none
	include 'headers.inc'
	include 'types.inc'
	record /HEADER/ hcb

	record /M_PTR/ name_ptr			! address pointer
	integer*4 n
	character*24 bname

	if (hcb.type/1000 .eq. 3) then
	  bname = hcb.subr(n).name
	elseif (hcb.type/1000 .eq. 5) then
	  bname = hcb.comp(n).name
	elseif (hcb.type/1000 .eq. 6) then
	  name_ptr.p = hcb.t6subr_loc.p + (n-1)*L_t6subr_len + T6O_name
	  call M$MOVE (%val(name_ptr.p), %ref(bname), 24)
	else
	  call M$ERROR ('Illegal type class in M$SUB_NAME')
	endif
	M$SUB_NAME = bname
	return
	end

c	character*2 function M$SUB_FORMAT (hcb,n)
	function M$SUB_FORMAT (hcb,n)
	implicit none
	include 'headers.inc'
	include 'types.inc'
	integer*4 n
	record /HEADER/ hcb

	record /M_PTR/ format_ptr			! address pointer
	character*2 bformat

	if (hcb.type/1000 .eq. 3) then
	  bformat = hcb.subr(n).format
	elseif (hcb.type/1000 .eq. 5) then
	  bformat = hcb.comp(n).format
	elseif (hcb.type/1000 .eq. 6) then
	  format_ptr.p = hcb.t6subr_loc.p + (n-1)*L_t6subr_len + T6O_format
	  call M$MOVE (%val(format_ptr.p), %ref(bformat), 2)
	else
	  call M$ERROR ('Illegal type class in M$SUB_FORMAT')
	endif
	M$SUB_FORMAT = bformat
	return
	end



	function M$SUB_OFFSET (hcb,n)
	implicit none
	include 'headers.inc'
	include 'types.inc'
	record /HEADER/ hcb
	record /T6SUBR_STRUCT/ t6subr
	integer*4 boff,i,n

	if (hcb.type.ge.3000 .and. hcb.type.le.3999) then
	  boff = hcb.subr(n).offset
	elseif (hcb.type.ge.5000 .and. hcb.type.le.5999) then
	  boff = 0
	  do i = 1,n-1
	    boff = boff + M$BPA(hcb.comp(i).format)
	  enddo
#if !_FDS
	elseif (hcb.type/1000 .eq. 6) then
	  call M$GET_T6SUBRECORD (hcb, n, t6subr)
	  boff = t6subr.offset
#endif
	else
	  call M$ERROR ('Illegal type class in M$SUB_OFFSET')
	endif
	M$SUB_OFFSET = boff
	return
	end


	subroutine M$ADDCOMP (hcb,name,format,type,units)
	implicit none
c
	include 'types.inc'
	include 'headers.inc'
	record /HEADER/ hcb
	character*(*) name,format
	integer*4 type,units,n
	integer*4 bpa
c
	n = hcb.components+1
	hcb.comp(n).name = name
	hcb.comp(n).format = M$TEST_FORMAT(format,' ')
	hcb.comp(n).type = type
	hcb.comp(n).units = units
	hcb.components = n
	if (n.eq.1) then
	  hcb.record_length = 0
	  hcb.format = hcb.subr(1).format
	else
	  if (hcb.format .ne. hcb.subr(n).format) hcb.format = 'NH'
	endif
	bpa = M$BPA( hcb.subr(n).format )
	if (bpa .gt. 0) then
	  hcb.record_length = hcb.record_length + bpa
	else
	  hcb.record_length = hcb.record_length + (-bpa + 7)/8
	endif
c
	return
	end


	subroutine M$ADDSUBR (hcb,name,format)
	implicit none
	include 'types.inc'
	include 'headers.inc'
	record /HEADER/ hcb
	character*(*) name,format
	integer*4 n
	integer*4 bpa
c
	entry M$ADDSREC (hcb,name,format)
c
	n = hcb.subrecords+1
	hcb.subr(n).name = name
	hcb.subr(n).format = M$TEST_FORMAT(format,' ')
	hcb.subr(n).offset = hcb.record_length
	hcb.subrecords = n
	if (n.eq.1) then
	  hcb.record_length = 0
	  hcb.format = hcb.subr(1).format
	else
	  if (hcb.format .ne. hcb.subr(n).format) hcb.format = 'NH'
	endif
	bpa = M$BPA( hcb.subr(n).format )
	if (bpa .gt. 0) then
	  hcb.record_length = hcb.record_length + bpa
	else
	  hcb.record_length = hcb.record_length + (-bpa + 7)/8
	endif
c
	return
	end


#if !_FDS
	subroutine M$ADDT6SUBR
     &	(hcb,name,format,minval,maxval,nel,units,uprefix)
	implicit none
c
	include 'types.inc'
	include 'headers.inc'
	record /HEADER/ hcb
c parameters
	character*(*) name,format
	real*8    minval, maxval
	integer*4 nel, units, uprefix, junk
c local
	record /T6SUBR_STRUCT/ t6subr
	integer*4 n, old_alloc
	record /M_PTR/ old_loc		! Address pointer to previous
					! T6000 subrecord memory buffer
	integer*4 bpa
c
	entry M$ADDT6SREC
     &	(hcb,name,format,minval,maxval,nel,units,uprefix)

c if this is a brand new file, make sure that the SUBREC_DEF keyword
c exists, and that memory has been allocated to hold the subrecords info
	if (hcb.t6subr_alloc .eq. 0) then
	   junk = M$REPL_KEYWORD(hcb, 'SUBREC_DEF', ' ')
	   call M$OPEN_T6SUBRECORDS (hcb)
	endif

	n = hcb.subrecords+1

c check to see if we need to add more memory to hold this subrecord info
c (this logic is identical to a block of code in HEADERMOD;
c is there a better place for it?)
	if (hcb.t6subr_alloc .lt. n*L_t6subr_len) then
	  old_loc.p   = hcb.t6subr_loc.p
	  old_alloc = hcb.t6subr_alloc

	  hcb.t6subr_alloc = (n*L_t6subr_len + 511)/512 * 512
	  call M$STMALLOC (hcb.t6subr_alloc, hcb.t6subr_loc.p)

	  if (old_alloc .gt. 0) then
	    call M$MOVE (%val(old_loc.p), %val(hcb.t6subr_loc.p), old_alloc)
	    call M$STMFREE (old_alloc, old_loc)
	  endif
	endif

	t6subr.name     = name
	t6subr.format   = M$TEST_FORMAT(format,' ')
	t6subr.minval   = minval
	t6subr.maxval   = maxval
	t6subr.num_elts = nel
	t6subr.units    = units
	t6subr.uprefix  = uprefix
	t6subr.offset   = hcb.record_length
	call M$PUT_T6SUBRECORD (hcb, n, t6subr)

	hcb.subrecords = n
	if (n.eq.1) then
	  hcb.record_length = 0
	  hcb.format = t6subr.format
	else
	  if (hcb.format .ne. t6subr.format) hcb.format = 'NH'
	endif
	bpa = M$BPA( t6subr.format )
	if (bpa .ge. 0) then
	  hcb.record_length = hcb.record_length + bpa * nel
	else
	  hcb.record_length = hcb.record_length + (((-bpa * nel) + 7) / 8)
	endif
	call M$UPDATE_T6SUBRECORDS (hcb)

	return
	end
#endif

c************
c  Description:	Type3/5000 IO routines
c************

	subroutine M$GRAB5N (hcb,buf,start,nget,ngot,list)
	implicit none
	include 'headers.inc'
	record /HEADER/ hcb,hcbt
	real*8		start
	integer*4	nget, ngot
	character*(*) list
	byte		buf(*)
	hcbt.main = hcb.main
	call M$TRANSFORM (hcbt,list,'ECR','CARTESIAN','METRIC')
	call M$GRABT (hcb,buf,start,nget,ngot,hcbt)
	return
	end

	subroutine M$GRABT (hcb,buf,start,nget,ngot,hcbt)
	implicit none
	include 'headers.inc'
	include 'types.inc'
	record /HEADER/ hcb,hcbt
	real*8		start
	integer*4	nget, ngot, nrec
	byte		buf(*)
	record /M_PTR/  lt
	call M$MALLOC (hcb.record_len*nget,lt.p)
	call M$GRAB (hcb,%val(lt.p),start,nget,ngot)
	nrec = ngot
	if (ngot.eq.-1) nrec = nget
	call M$TRANSFER (hcb,%val(lt.p),hcbt,buf,nrec,50)
	call M$FREE (lt)
	return
	end

	subroutine M$GRABXT (hcb,buf,ngot,hcbt)
	implicit none
	include 'headers.inc'
	include 'types.inc'
	record /HEADER/ hcb,hcbt
	integer*4	ngot,nrec
	byte		buf(*)
	record /M_PTR/  lt
	call M$MALLOC (hcb.record_len*hcb.xfer_len,lt.p)
	call M$GRABX (hcb,%val(lt.p),ngot)
	nrec = ngot
	if (ngot.eq.-1) nrec = hcb.xfer_len
	call M$TRANSFER (hcb,%val(lt.p),hcbt,buf,nrec,50)
	call M$FREE (lt)
	return
	end


	subroutine M$TRANSFORM (H,list,frame,type,units)
	implicit none
	include 'headers.inc'
	include 'types.inc'
	character*(*) frame,type,units,list
	character*80 clist
	character*5 name
	integer*4 it,iu,ll,n,jt,ju,i
	record /HEADER/ H,HL

	HL.main = H.main

	if (frame.eq.' ') then
	  ! leave alone
	else
	  H.frame_of_ref = frame
	endif

	if (type.eq.' ') then
	  it = -1	! leave alone
	elseif (type(1:3).eq.'CAR') then
	  it = 2
	elseif (type(1:3).eq.'SPH') then
	  it = 3
	elseif (type(1:3).eq.'CYL') then
	  it = 4
	elseif (type(1:3).eq.'ELL') then
	  it = 5
	elseif (type(1:3).eq.'GEO') then
	  it = 6
	else
	  call M$ERROR('Unknown coordinate system in M$TRANSFORM: '//type)
	endif

	if (units.eq.' ') then
	  iu = -1	! leave alone
	elseif (units(1:2).eq.'ME') then
	  iu = 1
	elseif (units(1:2).eq.'EN') then
	  iu = 2
	elseif (units(1:2).eq.'NA') then
	  iu = 3
	else
	  call M$ERROR('Unknown units system in M$TRANSFORM: '//units)
	endif

	ll = M$LENGTH(list)+1
	if (ll.gt.1) then		! create new comp structures
	  clist = list(1:ll-1)//','
	  n = 0
	  H.components = 0
	  H.record_len = 0
	  do while (M$PARSE(clist,n+1,name).gt.0)
	    n = n+1
	    if (name(1:1).eq.'+') then
	      i = M$SUB_FIND (HL,name(2:))
	    else
	      i = M$SUB_FIND (HL,name)
	      if (i.eq.0) call M$ERROR ('Component not found: '//name)
	    endif
	    if (i.gt.0) then
	      H.comp(n) = HL.comp(i)
	      H.record_len = H.record_len + M$BPA(H.comp(n).format)
 	      H.components = H.components+1
	    elseif (name.eq.'+TIME') then
	      call M$ADDCOMP (H,'TIME','SD',1,1)
	    elseif (name.eq.'+POS') then
	      call M$ADDCOMP (H,'POS ','VD',2,5)
	    elseif (name.eq.'+VEL') then
	      call M$ADDCOMP (H,'VEL ','VD',2,6)
	    elseif (name.eq.'+ACC') then
	      call M$ADDCOMP (H,'ACC ','VD',2,7)
	    else
	      call M$ERROR ('Dont know how create component: '//name)
	    endif
	  enddo
	endif

	do n = 1,H.components
	  jt = H.comp(n).type
	  ju = H.comp(n).units
	  if (it.gt.0 .and. jt.ge.2 .and. jt.le.6) jt = it
	  if (iu.gt.0) call M$CONVERT_UNITS (ju,iu,ju)
	  H.comp(n).type = jt
	  H.comp(n).units = ju
	enddo

	return
	end



	subroutine M$TRANSFER (HA,A,HB,B,nrec,mode)
	implicit none
c
	include 'midas.inc'
	include 'headers.inc'
	include 'types.inc'
	record /HEADER/ HA,HB
	record /T6SUBR_STRUCT/ t6subr
	include 'type5000.inc'
	record /TRANSMATRIX/ M
	byte A(*),B(*)
	integer*4 nrec,mode,rla,rlb,offa,offb,bpa,bpb
	integer*4 type,mpva,nc,i,n,ii,ua,ub,spa
	character forma*2, formb*2, blanks*8
	real*8 scale
	real*8 tpcquad(3) / -90.d0, 0.d0, 0.d0 /
        integer*4 NULLATOM(L_MAXSPA_) / L_MAXSPA_*0 /
	record /M_PTR/ la,lb,lt,lgha,lat

	la.p = %loc(A)
	lb.p = %loc(B)
	if (nrec.eq.0) return

	if (mode.eq.1) then		! move subrecord into
	  call M$MOVE (HA,n,4)
	  offb = M$SUB_OFFSET (HB,n)
	  formb = M$SUB_FORMAT(HB,n)
	  bpb = M$BPA(formb)
	  rlb = HB.record_length
	  call M$VMOVNB (%val(la.p),bpb,%val(lb.p+offb),rlb,nrec,bpb)
	  return
	elseif (mode.eq.-1) then	! extract subrecord from
	  call M$MOVE (HB,n,4)
	  offa = M$SUB_OFFSET (HA,n)
	  forma = M$SUB_FORMAT(HA,n)
	  bpa = M$BPA(forma)
	  rla = HA.record_length
	  call M$VMOVNB (%val(la.p+offa),rla,%val(lb.p),bpa,nrec,bpa)
	  return
	elseif (mode.eq.9 .or. mode.eq.-9) then	! convert machine rep
	  if (HA.data_rep .eq. Mc.machine) return
	  lt.p = 0
	  rla = HA.record_length
	  do n = 1,HA.subrecords
	    forma = M$SUB_FORMAT(HA,n)
	    if (M$SEARCH('ABOP',forma(2:2)).eq.0) then
	      offa = M$SUB_OFFSET (HA,n)
	      bpa = M$BPA(forma)
	      spa = M$NPM(forma(1:1))
	      call M$REMALLOC(nrec*bpa,lt.p)
	      call M$VMOVNB (%val(la.p+offa),rla,%val(lt.p),bpa,nrec,bpa)
	      if (mode.gt.0) then
		call M$CONVERT (HA.data_rep,%val(lt.p),forma(2:2),nrec*spa)
	      else
		call M$CONVERT_TO (HA.data_rep,%val(lt.p),forma(2:2),nrec*spa)
	      endif
	      call M$VMOVNB (%val(lt.p),bpa,%val(la.p+offa),rla,nrec,bpa)
	    endif
	  enddo
	  call M$FREE (lt)
	  return
	elseif ((mode.eq.5 .or. mode.eq.50) .and. (HA.type/1000.ne.5)) then
	  call M$ERROR('Cannot perform transformation on non-5000 file')
	endif

	lt.p = 0
	lgha.p = 0
	lat.p = 0
	rla = HA.record_length
	rlb = HB.record_length
	if (la.p.eq.lb.p) then  ! in place - copy buffer
	  call M$MALLOC (rla*nrec,lat.p)
	  call M$MOVE (%val(la.p),%val(lat.p),rla*nrec)
	  la.p = lat.p
	endif

	n = HB.subrecords
	do while (n.gt.0)
	  offb = M$SUB_OFFSET (HB,n)
	  formb = M$SUB_FORMAT(HB,n)
	  bpb = M$BPA(formb)

#if !_FDS
	  if (hb.class .eq. 6) then
	    call M$GET_T6SUBRECORD (HB, n, t6subr)
	    i = M$SUB_FIND (HA, t6subr.name)
	  else
#endif
	    i = M$SUB_FIND (HA,HB.subr(n).name)
#if !_FDS
	  endif
#endif

	  if (i.eq.0) then
	    if (mode.lt.10) then
	      ! no zero fill mode
	    elseif (HB.comp(n).name.eq.'TIME') then
	      call M$REMALLOC (nrec*max(8,bpb),lt.p)
	      scale = HB.tstart + HB.offset*HB.tdelta   ! tstart
	      call DRAMP (%val(lt.p),scale,HB.tdelta,nrec)
	      call M$REFORMAT (%val(lt.p),'SD',%val(lt.p),formb(2:2),nrec)
	      call M$VMOVNB (%val(lt.p),bpb,%val(lb.p+offb),rlb,nrec,bpb)
	    elseif (formb(2:2).eq.'A') then
	      call M$REMALLOC (nrec*bpb,lt.p)
	      call M$VMOVNB (%ref(' '),0,%val(lt.p),1,bpb,1)
	      call M$VMOVNB (%val(lt.p),0,%val(lb.p+offb),rlb,nrec,bpb)
	    else
	      call M$REMALLOC (nrec*bpb,lt.p)
c avl DR# 6361  call M$REFORMAT (0,'SL',%val(lt.p),formb(2:2),M$SPA(formb))
	      call M$REFORMAT (NULLATOM,'SL',%val(lt.p),formb(2:2),M$SPA(formb))
	      call M$VMOVNB (%val(lt.p),0,%val(lb.p+offb),rlb,nrec,bpb)
	    endif
	  else
	    offa = M$SUB_OFFSET (HA,i)
	    forma = M$SUB_FORMAT(HA,i)
	    bpa = M$BPA(forma)
	    if ((mode.eq.5 .or. mode.eq.50) .and.
     &	         (HA.comp(i).type .ne. 0) ) then
c type 5000 transformation section
c special case for non-vectors
	      ua = HA.comp(i).units
	      ub = HB.comp(n).units
	      if (HA.comp(i).type .eq. P_SCA) then
	        call M$CONVERSION_FACTOR (ua,ub,scale)
	        if (scale.eq.1.d0) goto 51	! just copy
	        call M$REMALLOC(nrec*max(bpa,bpb,8),lt.p)
	        call M$VMOVNB (%val(la.p+offa),rla,%val(lt.p),bpa,nrec,bpa)
	        call M$REFORMAT(%val(lt.p),forma,%val(lt.p),'D',nrec)
	        call M$DVSCALE (%val(lt.p),1,scale,nrec)
	        call M$REFORMAT(%val(lt.p),'SD',%val(lt.p),formb(2:2),nrec)
	        call M$VMOVNB (%val(lt.p),bpb,%val(lb.p+offb),rlb,nrec,bpb)
	        goto 52
	      endif

c determine special vector mode for (POS,VEL,ACC)
	      mpva = M$SEARCH ('POS,VEL,ACC',HB.comp(n).name(1:3))
	      if (mpva.ne.0) mpva = (mpva+3)/4  ! modes 0,1,2,3
	      nc = max(1,mpva)  ! number of components in group
	      if (nc.gt.1) then
	        n = n - nc + 1    ! base calcs on POS component
	        i = i - nc + 1	  ! assume input vectors are there also
	        if (mpva.eq.3 .and. HB.comp(n+1).name(1:3).ne.'VEL')
     +	        call M$ERROR('ACC vector must be preceded by VEL vector')
	        if (mpva.ge.2 .and. HB.comp(n).name(1:3).ne.'POS')
     +	        call M$ERROR('VEL vector must be preceded by POS vector')
	        offa = offa - bpa*(nc-1)
	        offb = offb - bpa*(nc-1)
	      endif
	      if (forma(1:1).ne.'V') call M$ERROR
     +	      ('Non-vector assigned to a vector coordinate system type')

c reformat to double precision data format
	      call M$REMALLOC(nrec*nc*max(bpa,bpb,24),lt.p)
	      call M$VMOVNB (%val(la.p+offa),rla,%val(lt.p),bpa*nc,nrec,bpa*nc)
	      call M$REFORMAT(%val(lt.p),forma,%val(lt.p),'D',nrec*nc)

c complete TOPOCENT quadword fields
	      if (HA.frame_of_ref.eq.'TOPOCENT')
     +	        call M$MOVE (tpcquad,HA.quadwords(5),24)
	      if (HB.frame_of_ref.eq.'TOPOCENT')
     +	        call M$MOVE (tpcquad,HB.quadwords(5),24)

c transform to CARTESIAN,METRIC
	      call M$CONVERSION_FACTOR (ua,0,scale)
	      if     (HA.comp(i).type .eq. P_CAR) then
	        call M$DVSCALE (%val(lt.p),1,scale,nc*nrec*3)
	      elseif (HA.comp(i).type .eq. P_SPH) then
	        call M$DVSCALE (%val(lt.p),3,scale,nc*nrec)
	        call M$SPH_TO_CAR (%val(lt.p),mpva,nrec)
	      elseif (HA.comp(i).type .eq. P_GEO) then
	        call M$DVSCALE (%val(lt.p),3,scale,nc*nrec)
	        call M$GEO_TO_CAR (%val(lt.p),mpva,nrec)
	      else
	        call M$ERROR('Unsupported coordinate system conversion')
	      endif

c flip reference frames
	      if     (HA.frame_of_ref .eq. 'ECR') then
	        if     (HB.frame_of_ref .eq. 'ECR') then
	        elseif (HB.frame_of_ref .eq. 'ECI') then
	          call GET_GHA (HA,A,lgha,nrec,HB.quadwords(10))
	          call M$ECR_TO_ECI (%val(lt.p),mpva,nrec,%val(lgha.p))
	        elseif (HB.frame_of_ref .eq. 'TOP'  .or.
     +	      	  HB.frame_of_ref .eq. 'TOPOCENT') then
	          call M$ECR_TO_TOP (%val(lt.p),mpva,nrec,HB.quadwords(2))
	        else
	          call M$ECR_TO_CRF (%val(lt.p),mpva,nrec,HB.quadwords)
	        endif
	      elseif (HA.frame_of_ref .eq. 'ECI') then
	        if     (HB.frame_of_ref .eq. 'ECR') then
	          call GET_GHA (HA,A,lgha,nrec,HA.quadwords(10))
	          call M$ECI_TO_ECR (%val(lt.p),mpva,nrec,%val(lgha.p))
	        elseif (HB.frame_of_ref .eq. 'ECI') then
	        elseif (HB.frame_of_ref .eq. 'TOP'  .or.
     +	      	  HB.frame_of_ref .eq. 'TOPOCENT') then
	          call GET_GHA (HA,A,lgha,nrec,HA.quadwords(10))
	          call M$ECI_TO_ECR (%val(lt.p),mpva,nrec,%val(lgha.p))
	          call M$ECR_TO_TOP (%val(lt.p),mpva,nrec,HB.quadwords(2))
	        else
	          call GET_GHA (HA,A,lgha,nrec,HA.quadwords(10))
	          call M$ECI_TO_ECR (%val(lt.p),mpva,nrec,%val(lgha.p))
	          call M$ECR_TO_CRF (%val(lt.p),mpva,nrec,HB.quadwords)
	        endif
	      elseif (HA.frame_of_ref .eq. 'TOP'  .or.
     +	      	HA.frame_of_ref .eq. 'TOPOCENT') then
	        if     (HB.frame_of_ref .eq. 'ECR') then
	          call M$TOP_TO_ECR (%val(lt.p),mpva,nrec,HA.quadwords(2))
	        elseif (HB.frame_of_ref .eq. 'ECI') then
	          call M$TOP_TO_ECR (%val(lt.p),mpva,nrec,HA.quadwords(2))
	          call GET_GHA (HA,A,lgha,nrec,HB.quadwords(10))
	          call M$ECR_TO_ECI (%val(lt.p),mpva,nrec,%val(lgha.p))
	        elseif (HB.frame_of_ref .eq. 'TOP'  .or.
     +	      	  HB.frame_of_ref .eq. 'TOPOCENT') then
	          ! only if TOPO parameters are not the same
	          if (M$CHK(HA.quadwords,HB.quadwords,7*8).ne.1)
     &	          call M$TOP_TO_TOP (%val(lt.p),mpva,nrec,
     &	      	HA.quadwords(2),HB.quadwords(2))
	        else
	          call M$TOP_TO_ECR (%val(lt.p),mpva,nrec,HA.quadwords(2))
	          call M$ECR_TO_CRF (%val(lt.p),mpva,nrec,HB.quadwords)
	        endif
	      else
	        if     (HB.frame_of_ref .eq. 'ECR') then
	          call M$CRF_TO_ECR (%val(lt.p),mpva,nrec,HA.quadwords)
	        elseif (HB.frame_of_ref .eq. 'ECI') then
	          call M$CRF_TO_ECR (%val(lt.p),mpva,nrec,HA.quadwords)
	          call GET_GHA (HA,A,lgha,nrec,HB.quadwords(10))
	          call M$ECR_TO_ECI (%val(lt.p),mpva,nrec,%val(lgha.p))
	        elseif (HB.frame_of_ref .eq. 'TOP'  .or.
     +	      	  HB.frame_of_ref .eq. 'TOPOCENT') then
	          call M$CRF_TO_ECR (%val(lt.p),mpva,nrec,HA.quadwords)
	          call M$ECR_TO_TOP (%val(lt.p),mpva,nrec,HB.quadwords(2))
	        else
	          ! only if custom frames are not the same
	          if (M$CHK(HA.quadwords,HB.quadwords,8).ne.1)
     &	          call M$CRF_TO_CRF (%val(lt.p),mpva,nrec,
     &	      	HA.quadwords,HB.quadwords)
	        endif
	      endif

c transform from CARTESIAN,METRIC
	      call M$CONVERSION_FACTOR (0,ub,scale)
	      if     (HB.comp(n).type .eq. P_CAR) then
	        call M$DVSCALE (%val(lt.p),1,scale,nc*nrec*3)
	      elseif (HB.comp(n).type .eq. P_SPH) then
	        call M$CAR_TO_SPH (%val(lt.p),mpva,nrec)
	        call M$DVSCALE (%val(lt.p),3,scale,nc*nrec)
	      elseif (HB.comp(n).type .eq. P_GEO) then
	        call M$CAR_TO_GEO (%val(lt.p),mpva,nrec)
	        call M$DVSCALE (%val(lt.p),3,scale,nc*nrec)
	      else
	        call M$ERROR('Unsupported coordinate system conversion')
	      endif

c reformat to desired output format
	      call M$REFORMAT(%val(lt.p),'VD',%val(lt.p),formb(2:2),nrec*nc)
	      call M$VMOVNB (%val(lt.p),bpb*nc,%val(lb.p+offb),rlb,nrec,bpb*nc)
	      goto 52
	    endif

 51	    if (forma.eq.formb) then	! copy subrecord
	      call M$VMOVNB (%val(la.p+offa),rla,%val(lb.p+offb),rlb,nrec,bpa)
	    elseif (forma(2:2).eq.'A' .and. formb(2:2).eq.'A') then
	      call M$VMOVNB (%val(la.p+offa),rla,%val(lb.p+offb),rlb,
     &				nrec,min(bpa,bpb))
	      blanks = ' '	! fill strings
	      do i = bpa, bpb-8, 8
		call M$VMOVNB (%ref(blanks),0,%val(lb.p+offb+i),rlb,nrec,8)
	      enddo
	    else			! reformat subrecord
	      call M$REMALLOC(nrec*max(bpa,bpb),lt.p)
	      call M$VMOVNB (%val(la.p+offa),rla,%val(lt.p),bpa,nrec,bpa)
	      call M$REFORMAT(%val(lt.p),forma,%val(lt.p),formb(2:2),nrec)
	      call M$VMOVNB (%val(lt.p),bpb,%val(lb.p+offb),rlb,nrec,bpb)
 	    endif
 52	    continue
	  endif
	  n = n-1
	enddo

	call M$FREE (lt)
	call M$FREE (lat)
	call M$FREE (lgha)
	return


	end



c************
c  Description:	Type5000 coordinate system transformations
c************

	subroutine M$SPH_TO_CAR (V,mode,nrec)
	! radius,azimuth,elevation -> x,y,z
	implicit none
	include 'constants.inc'
	integer*4 mode,nrec,i,isk,n
	real*8 V(*),ca,cb,sa,sb,r,a,b,rr,aa,bb,ri,rrr,aaa,bbb,rri
	i = 0
	isk = mode*3
	do n = 1,nrec
	  r = v(i+1)
	  a = v(i+2)*ddeg2rad
	  b = v(i+3)*ddeg2rad
	  ca = dcos(a)
	  sa = dsin(a)
	  cb = dcos(b)
	  sb = dsin(b)
	  v(i+1) = r*cb*ca
	  v(i+2) = r*cb*sa
	  v(i+3) = r*sb
	  if (mode.le.1) goto 10
	  rr = v(i+4)
	  aa = v(i+5)*ddeg2rad
	  bb = v(i+6)*ddeg2rad
	  ri = rr/r
	  v(i+4) = v(i+1)*ri - v(i+3)*bb*ca - v(i+2)*aa
	  v(i+5) = v(i+2)*ri - v(i+3)*bb*sa + v(i+1)*aa
	  v(i+6) = v(i+3)*ri + r*bb*cb
	  if (mode.le.2) goto 10
	  rrr = v(i+7)
	  aaa = v(i+8)*ddeg2rad
	  bbb = v(i+9)*ddeg2rad
	  rri = rrr/r
	  v(i+7) = (v(i+4)-v(i+1)*ri)*ri + v(i+1)*rri - (v(i+6)*bb
     &	    +v(i+3)*bbb)*ca + v(i+3)*bb*aa*sa - v(i+2)*aaa - v(i+5)*aa
	  v(i+8) = (v(i+5)-v(i+2)*ri)*ri + v(i+2)*rri - (v(i+6)*bb
     &	    +v(i+3)*bbb)*sa - v(i+3)*bb*aa*ca - v(i+1)*aaa + v(i+4)*aa
	  v(i+9) = (v(i+6)-v(i+3)*ri)*ri + v(i+2)*rri + (rr*bb+r*bbb)*ca
     &	    + v(i+3)*bb*bb
  10	  i = i+isk
	enddo
	return
	end

	subroutine M$CAR_TO_SPH (V,mode,nrec)
	! x,y,z -> radius,azimuth,elevation
	implicit none
	include 'constants.inc'
	integer*4 mode,nrec,i,isk,n
	real*8 V(*),x,y,z,xx,yy,zz,xxx,yyy,zzz,ri,si,si2
	i = 0
	isk = mode*3
	do n = 1,nrec
	  x = v(i+1)
	  y = v(i+2)
	  z = v(i+3)
	  v(i+1) = dsqrt(x*x+y*y+z*z)
	  ri = 1.d0/v(i+1)
	  v(i+2) = datan2d(y,x)
	  v(i+3) = dasind(z*ri)
	  if (mode.le.1) goto 10
	  xx = v(i+4)
	  yy = v(i+5)
	  zz = v(i+6)
	  si2 = 1.d0/(x*x+y*y)
	  si = dsqrt(si2)
	  v(i+4) = (x*xx+y*yy+z*zz)*ri
	  v(i+5) = (x*yy-xx*y)*si2*drad2deg
	  v(i+6) = (zz-z*v(i+4)*ri)*si*drad2deg
	  if (mode.le.2) goto 10
	  xxx = v(i+7)
	  yyy = v(i+8)
	  zzz = v(i+9)
	  v(i+7) = (xx*xx+yy*yy+zz*zz + x*xxx+y*yyy+z*zzz - v(i+4)**2)*ri
	  v(i+8) = ((x*yyy-xxx*y)*drad2deg-2*v(i+5)*(x*xx+y*yy))*si2
	  v(i+9) = (zzz-zz*v(i+4)*ri-z*(v(i+7)*ri-(v(i+4)*ri)**2))*si*drad2deg
     &		 - v(i+6)*(x*xx+y*yy)*si2
  10	  i = i+isk
	enddo
	return
	end

	subroutine M$GEO_TO_CAR (V,mode,nrec)
	! altitude,latitude,longitude -> x,y,z
	implicit none
	include 'constants.inc'
	integer*4 mode,nrec,i,isk,n
	real*8 V(*),ca,cb,sa,sb,r,a,b,ee,t
	ee = earth_eccentricity**2
	i = 0
	isk = mode*3
	do n = 1,nrec
	  r = v(i+1)
	  a = v(i+2)*ddeg2rad
	  b = v(i+3)*ddeg2rad
	  ca = dcos(a)
	  sa = dsin(a)
	  cb = dcos(b)
	  sb = dsin(b)
	  t = earth_equatorial_radius/dsqrt(1.d0-ee*sa**2)
	  v(i+1) = (r+t)*ca*cb
	  v(i+2) = (r+t)*ca*sb
	  v(i+3) = (r+(1-ee)*t)*sa
	  if (mode.gt.1) call M$ERROR
     &	  ('VEL,ACC Geodetic transforms not supported yet')
  10	  i = i+isk
	enddo
	return
	end

	subroutine M$CAR_TO_GEO (V,mode,nrec)
	! x,y,z, -> altitude,latitude,longitude
	implicit none
	include 'constants.inc'
	integer*4 mode,nrec,i,isk,n
	real*8 V(*),x,y,z,a,b,c,eer,ee,ee1,ee1i,r,t,f1,f2,sa,ca,h,da,dh
        real*8 oldda

	if (mode.gt.1) call M$ERROR
     &      ('VEL,ACC Geodetic transforms not supported yet')

c++DR 15149 JLR Initializing da which is saved to oldda through each iteration.
        da = 0.0
        ee = earth_eccentricity**2
	ee1 = (1-ee)
	ee1i = 1./ee1
	eer = earth_equatorial_radius
	i = 0
	isk = mode*3
	do n = 1,nrec
c++DR 15149 JLR Initializing da which is saved to oldda through each iteration.
          da = 0.0

	  x = v(i+1)
	  y = v(i+2)
	  z = v(i+3)
	  if (x.eq.0 .and. y.eq.0) then
	    v(i+1) = abs(z) - eer*dsqrt(1-ee)
	    v(i+2) = sign (90.d0, z)
	    v(i+3) = 0
	    goto 10
	  endif
	  r = dsqrt(x*x+y*y+z*z)
	  b = dsqrt(x*x+y*y)
	  c = dasin(z/r)
	  a = datan( dtan(c)*ee1i )
	  sa = dsin(a)
	  h = r - eer/dsqrt(1+ee*sa*sa*ee1i)
	  goto 2
 1	  sa = dsin(a)			! top of iteration loop
 2	  ca = dcos(a)
	  t  = 1./dsqrt(1-ee*sa*sa)
	  f1 = (h + eer*t)*ca - b
	  f2 = (h + eer*ee1*t)*sa - z
          oldda = da
	  da = (-f1*sa+f2*ca)/(h+eer*ee1*t*t*t)
	  dh = f1*ca+f2*sa
          a = a-da
          h = h-dh

          if ((oldda.eq.da) .and. (oldda.ne.0)) then
	     write(*,201) n, nrec, x, y, z, oldda, da
	     call M$WARNING
     &	  ('Invalid data - possibly a corrupted type 5000 file')
	     call M$WARNING('Skipping record')
	     goto 10                  ! ignore this record; skip to next
	  endif

	  if (da.gt.0.0001) goto 1	! iterate on lat tolerance
	  if (dh.gt.0.0001) goto 1	! iterate on height tolerance
	  v(i+1) = h
	  v(i+2) = a*drad2deg
	  v(i+3) = datan2d(y,x)
  10	  i = i+isk
	enddo
	return
 201	format('M$CAR_TO_GEO, input record ', I4, ' of ', I4, ':'/
     &         ' x = ', 1PG20.12, ';'/
     &         ' y = ', 1PG20.12, ';'/
     &         ' z = ', 1PG20.12, ';'/
     &         'da = ', 1PG20.12, ', oldda = ', 1PG20.12)
	end


c************
c  Description:	Type5000 frame of reference transformations
c************

	subroutine M$CRF_TO_CRF (V,mode,nrec,a,b)
	implicit none
	integer*4 mode,nrec
	include 'type5000.inc'
	record /TRANSMATRIX/ M,MA,MB
	character*80 filename
	character*8 name
	byte		a(*)
	byte		b(*)
	byte		V(*)
	call M$MOVE (b,%ref(name),8)
	if (name.ne.'NULL') call M$ERROR
     &	('Custom reference frames not supported yet')
	entry M$CRF_TO_ECR (V,mode,nrec,a)
	entry M$ECR_TO_CRF (V,mode,nrec,a)
	call M$MOVE (a,%ref(name),8)
	if (name.ne.'NULL') call M$ERROR
     &	('Custom reference frames not supported yet')
	return
	end

	subroutine M$TOP_TO_TOP (V,mode,nrec,a,b)
	implicit none
	include 'type5000.inc'
	integer*4 mode,nrec
	record /TRANSMATRIX/ M,MA,MB
	byte		a(*)
	byte		b(*)
	byte		V(*)
	call M$MATRIX_TOP2ECR (MA,a)
	call M$MATRIX_TOP2ECR (MB,b)
	call M$MATRIX_TRN (MB)
	call M$MATRIX_MUL (MB,MA,M)
	call M$MATRIX_MULPVA (M,V,mode,nrec)
	return
	entry M$TOP_TO_ECR (V,mode,nrec,a)
	call M$MATRIX_TOP2ECR (M,a)
	call M$MATRIX_MULPVA (M,V,mode,nrec)
	return
	entry M$ECR_TO_TOP (V,mode,nrec,a)
	call M$MATRIX_TOP2ECR (M,a)
	call M$MATRIX_TRN (M)
	call M$MATRIX_MULPVA (M,V,mode,nrec)
	return
	end

	subroutine M$MATRIX_TOP2ECR (M,a)
	implicit none
	include 'type5000.inc'
	record /TRANSMATRIX/ M
	include 'constants.inc'
	real*8 t1,t2,a(*),alt,lat,lon,slat,clat,slon,clon
	if (a(2).gt.90.0 .or. a(2).lt.-90.0) call M$WARNING
     &	('Latitude out of allowable range -90 < lat < 90 in TOP2ECR')
	alt = a(1)
	lat = a(2)*ddeg2rad
	lon = a(3)*ddeg2rad
	slat = dsin(lat)
	clat = dcos(lat)
	slon = dsin(lon)
	clon = dcos(lon)
	t1 = earth_eccentricity**2
	t2 = earth_equatorial_radius/dsqrt(1.d0-t1*slat**2)
	call M$MATRIX_IDN (M)
	call M$MATRIX_ROT (M,a(6),P_ROLL)	! local ROLL
	call M$MATRIX_ROT (M,a(5),P_ELEV)	! local ELEVation
	call M$MATRIX_ROT (M,a(4),P_AZIM)	! local AZIMuth
	call M$MATRIX_ROT (M,90-a(2),P_ROLL)	! topo roll rot
	call M$MATRIX_ROT (M,a(3)+90,P_AZIM)	! topo azimuth rot
	M.b1 = (alt+t2)*clat*clon		! topo origin
	M.b2 = (alt+t2)*clat*slon
	M.b3 = (alt+(1-t1)*t2)*slat
	return
	end

	subroutine M$ECI_TO_ECR (V,mode,nrec,gha)
	implicit none
	include 'constants.inc'
	integer*4 mode,nrec,i,isk,dir,n
	real*8 v(*),gha(*),eav,eavt,eav2,x,y,xx,yy,xxx,yyy,ca,sa
	dir = 1
	goto 1
	entry M$ECR_TO_ECI (V,mode,nrec,gha)
	dir = -1
  1	eav = dir*earth_angular_velocity
	eavt = eav+eav
	eav2 = eav*eav
	i = 0
	isk = 3*mode
	do n = 1,nrec
	  sa = dsin(gha(n))
	  ca = dcos(gha(n))
	  if (dir.lt.0) sa = -sa
	  x = v(i+1)
	  y = v(i+2)
	  v(i+1) =  ca*x + sa*y
	  v(i+2) = -sa*x + ca*y
	  if (mode.le.1) goto 10
	  xx = v(i+4)
	  yy = v(i+5)
	  v(i+4) =  ca*xx + sa*yy + eav*(-sa*x+ca*y)
	  v(i+5) = -sa*xx + ca*yy + eav*(-ca*x-sa*y)
	  if (mode.le.2) goto 10
	  xxx = v(i+7)
	  yyy = v(i+8)
	  v(i+7) =  ca*xxx + sa*yyy + eavt*(-sa*xx+ca*yy) + eav2*(-ca*x-sa*y)
	  v(i+8) = -sa*xxx + ca*yyy + eavt*(-ca*xx-sa*yy) + eav2*(+sa*x-ca*y)
  10	  i = i + isk
	enddo
	return
	end

	subroutine GET_GHA (H,A,lgha,nrec,quads)
	implicit none
	include 'constants.inc'
	include 'headers.inc'
	include 'types.inc'
	record /HEADER/ H
	integer*4 nrec,bpa,i
	real*8 A(*),ghas,ghad,quads(3)
	record /M_PTR/ la, lgha
	call M$REMALLOC (nrec*8,lgha.p)
	i = M$SUB_FIND(H,'TIME')
	if (i.eq.0) then
	  ghad = earth_angular_velocity*H.xdelta
	  ghas = quads(3) + earth_angular_velocity*H.xstart
	  ghas = ghas + (H.offset-nrec)*ghad
	  call DRAMP (%val(lgha.p),ghas,ghad,nrec)
	else
	  bpa = M$BPA(H.comp(i).format)
	  la.p = %loc(A) + M$SUB_OFFSET(H,i)
	  call M$VMOVNB (%val(la.p),H.record_len,%val(lgha.p),bpa,nrec,bpa)
	  call M$REFORMAT(%val(lgha.p),H.comp(i).format,%val(lgha.p),'D',nrec)
	  call M$DVSCALE (%val(lgha.p),1,earth_angular_velocity,nrec)
	  ghas = quads(3) + earth_angular_velocity*H.xstart
	  call M$DVADD   (%val(lgha.p),1,ghas,nrec)
	endif
	return
	end

	subroutine DRAMP (A,S,D,nrec)
	implicit none
	integer*4 n,nrec
	real*8 A(*),S,D
	A(1) = S
	do n = 2,nrec
	  A(n) = A(n-1) + D
	enddo
	return
	end

c************
c  Description:	Type5000 matrix utilities
c************

	subroutine M$MATRIX_ROT (M,angle,type)
	implicit none
	include 'type5000.inc'
	include 'constants.inc'
	record /TRANSMATRIX/ M,MA,MB
	real*8 angle,a,sa,ca
	integer*4 type
c
	if (angle.eq.0.d0) return
	a = angle*ddeg2rad
	sa = dsin(a)
	ca = dcos(a)
	MB = M
	call M$MATRIX_IDN (MA)
	if (type.eq.P_AZIM .or. type.eq.P_YAW) then
	  MA.a11 = ca
	  MA.a12 = -sa
	  MA.a21 = sa
	  MA.a22 = ca
	elseif (type.eq.P_ELEV .or. type.eq.P_PITCH) then
	  if (type.eq.P_PITCH) sa = -sa
	  MA.a11 = ca
	  MA.a13 = -sa
	  MA.a31 = sa
	  MA.a33 = ca
	elseif (type.eq.P_ROLL) then
	  MA.a22 = ca
	  MA.a23 = -sa
	  MA.a32 = sa
	  MA.a33 = ca
	else
	  call M$ERROR('Illegal type in M$MATRIX_ROT')
	endif
	call M$MATRIX_MUL (MA,MB,M)
	return
	end


	subroutine M$MATRIX_CPY (a,b)
	implicit none
	byte a(*)
	byte b(*)
	call M$MOVE (a,b,15*8)
	return
	end

	subroutine M$MATRIX_IDN (a)
	implicit none
	byte a(*)
	call M$VMOVN(0.d0,0,a,1,15,8)	! zero all
	call M$VMOVN(1.d0,0,a,4,3,8)	! set diagonal to 1
	return
	end

	subroutine M$MATRIX_MUL (a,b,c)
	implicit none
	integer*4 i,j
	real*8 a(3,5),b(3,5),c(3,5)
	do i = 1,3
	  do j = 1,3
	    c(i,j) = a(i,1)*b(1,j) + a(i,2)*b(2,j) + a(i,3)*b(3,j)
	  enddo
	  c(i,4) = a(i,4) + a(i,1)*b(1,4) + a(i,2)*b(2,4) + a(i,3)*b(3,4)
	enddo
	return
	end

	subroutine M$MATRIX_TRN (a)
	implicit none
	real*8 a(3,5),b(3),t
	! calculate transpose inplace
	t = a(1,2)
	a(1,2) = a(2,1)
	a(2,1) = t
	t = a(1,3)
	a(1,3) = a(3,1)
	a(3,1) = t
	t = a(3,2)
	a(3,2) = a(2,3)
	a(2,3) = t
	! invert scalar array
	call M$MOVE (a(1,4),b,24)
	a(1,4) = -( a(1,1)*b(1) + a(1,2)*b(2) + a(1,3)*b(3) )
	a(2,4) = -( a(2,1)*b(1) + a(2,2)*b(2) + a(2,3)*b(3) )
	a(3,4) = -( a(3,1)*b(1) + a(3,2)*b(2) + a(3,3)*b(3) )
	return
	end

	subroutine M$MATRIX_INV (a,b)
	implicit none
	real*8 a(3,5),b(3,5),det
	integer i,j
	! calculate sub-determinants
	b(1,1) = +( a(2,2)*a(3,3) - a(2,3)*a(3,2) )
	b(1,2) = -( a(1,2)*a(3,3) - a(1,3)*a(3,2) )
	b(1,3) = +( a(1,2)*a(2,3) - a(1,3)*a(2,2) )
	b(2,1) = -( a(2,1)*a(3,3) - a(2,3)*a(3,1) )
	b(2,2) = +( a(1,1)*a(3,3) - a(1,3)*a(3,1) )
	b(2,3) = -( a(1,1)*a(2,3) - a(1,3)*a(2,1) )
	b(3,1) = +( a(2,1)*a(3,2) - a(2,2)*a(3,1) )
	b(3,2) = -( a(1,1)*a(3,2) - a(1,2)*a(3,1) )
	b(3,3) = +( a(1,1)*a(2,2) - a(1,2)*a(2,1) )
	! calculate determinant
	det = a(1,1)*b(1,1) + a(1,2)*b(2,1) + a(1,3)*b(3,1)
	! scale matrix
	if (det.ne.0.d0) det = 1.d0/det
D	write(*,*) 'det = ',det
	do i = 1,3
	  do j = 1,3
	    b(i,j) = det*b(i,j)
	  enddo
	enddo
	! invert scalar array
	b(1,4) = -( b(1,1)*a(1,4) + b(1,2)*a(2,4) + b(1,3)*a(3,4) )
	b(2,4) = -( b(2,1)*a(1,4) + b(2,2)*a(2,4) + b(2,3)*a(3,4) )
	b(3,4) = -( b(3,1)*a(1,4) + b(3,2)*a(2,4) + b(3,3)*a(3,4) )
	return
	end

	subroutine M$MATRIX_DUMP (a)
	implicit none
	integer*4 i
	real*8 a(3,5)
	write(*,*) 'A1x = ',(a(1,i),i=1,3)
	write(*,*) 'A2x = ',(a(2,i),i=1,3)
	write(*,*) 'A3x = ',(a(3,i),i=1,3)
	write(*,*) 'Bx  = ',(a(i,4),i=1,3)
	write(*,*) 'Sx  = ',(a(i,5),i=1,3)
	return
	end

	subroutine M$MATRIX_MULPVA (A,V,mode,nrec)
	implicit none
	real*8 V(*),A(*)
	integer*4 nrec,mode
	if (mode.eq.0) then		! std vector
	  call M$MATRIX_MULA(A,V,3,nrec)
	elseif (mode.eq.1) then		! POS
	  call M$MATRIX_MULAB(A,V,3,nrec)
	elseif (mode.eq.2) then		! POS,VEL
	  call M$MATRIX_MULAB(A,V,6,nrec)
	  call M$MATRIX_MULA(A,V(4),6,nrec)
	elseif (mode.eq.3) then		! POS,VEL,ACC
	  call M$MATRIX_MULAB(A,V,9,nrec)
	  call M$MATRIX_MULA(A,V(4),9,nrec)
	  call M$MATRIX_MULA(A,V(7),9,nrec)
	else
	  call M$ERROR ('Illegal mode in MATRIX_MULPVA')
	endif
	return
	end

	subroutine M$MATRIX_MULA (A,V,isk,nrec)
	implicit none
	integer*4 i, nrec, isk, n
	real*8 V(*),A(*),vx,vy,vz
	i = 0
	do n = 1,nrec
	  vx = V(i+1)
	  vy = V(i+2)
	  vz = V(i+3)
	  V(i+1) = vx*A(1) + vy*A(4) + vz*A(7)
	  V(i+2) = vx*A(2) + vy*A(5) + vz*A(8)
	  V(i+3) = vx*A(3) + vy*A(6) + vz*A(9)
	  i = i+isk
	enddo
	return
	end

	subroutine M$MATRIX_MULAB (A,V,isk,nrec)
	implicit none
	integer*4 i, nrec, isk, n
	real*8 V(*),A(*),vx,vy,vz
	i = 0
	do n = 1,nrec
	  vx = V(i+1)
	  vy = V(i+2)
	  vz = V(i+3)
	  V(i+1) = vx*A(1) + vy*A(4) + vz*A(7) + A(10)
	  V(i+2) = vx*A(2) + vy*A(5) + vz*A(8) + A(11)
	  V(i+3) = vx*A(3) + vy*A(6) + vz*A(9) + A(12)
	  i = i+isk
	enddo
	return
	end

c************
c  Description:	Type5000 transform utilities
c************

	subroutine M$DVSCALE (a,isk,scale,nrec)
	implicit none
	real*8 a(*),scale
	integer*4 nrec,n,isk
	if (scale.eq.1.d0) return
	do n = 1,nrec*isk,isk
	  a(n) = a(n)*scale
	enddo
	return
	end

	subroutine M$DVADD (a,isk,add,nrec)
	implicit none
	real*8 a(*),add
	integer*4 isk,nrec,n
	if (add.eq.0.d0) return
	do n = 1,nrec*isk,isk
	  a(n) = a(n) + add
	enddo
	return
	end


	subroutine M$CONVERT_UNITS (units1,type,units2)
	implicit none
	include 'types.inc'
	integer*4 units1,units2,type
	character*4 s
	byte una(5,3) / 5,6,7,8,63,
     &			40,42,45,49,62,
     &			41,43,46,49,62 /
	units2 = units1
	if (type.eq.0) return
	if (type.lt.0 .or. type.gt.3) call M$ERROR
     &	('Units type parameter out of range in M$CONVERT_UNITS')
	s = M$UNITS_NAME (units1)
	if     (s.eq.'Dist') then
	  units2 = una(1,type)
	elseif (s.eq.'Spee') then
	  units2 = una(2,type)
	elseif (s.eq.'Acce') then
	  units2 = una(3,type)
	elseif (s.eq.'Jerk') then
	  units2 = una(4,type)
	elseif (s.eq.'Alti') then
	  units2 = una(5,type)
	endif
	return
	end



	subroutine M$CONVERSION_FACTOR (units1,units2,factor)
	implicit none
	include 'constants.inc'
	integer*4 units1,units2,i1,i2
	real*8 factor
	save
	real*8 fm(8)
	byte bfi(0:64) / 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,   ! 0-20
     &			1,1,2,1,3,1,1,1,1,1,1,1,2,1,1,1,1,1,1,4,5,4,  ! 21-42
     &			5,6,4,5,6,7,7,3,8,2,1,2,1,1,1,1,1,1,1,4,1,1 / ! 43-64
	factor = 1.d0
	if (units1.eq.units2) return
	i1 = bfi(units1)
	i2 = bfi(units2)
	if (i1.eq.i2) return
	if (fm(1).ne.1.d0) then
	  fm(1) = 1.d0
	  fm(2) = drad2deg
	  fm(3) = 360.d0
	  fm(4) = feet2meters
	  fm(5) = nm2meters
	  fm(6) = nm2meters/60.d0
	  fm(7) = 9.8d0
	  fm(8) = 3600.d0
	endif
	factor = fm(i1)/fm(i2)
	return
	end


	subroutine M$EPOCH_HOUR_ANGLE (time,year,second,rag)
	implicit none
	real*8 time,year,second,rag,ds50,ragw
	integer*4 mjd1950_jan1,mjd_jan1,nd,nm,ny,mjd
	mjd = (time+1d-7)/86400.d0
	if (time.lt.0.d0) mjd = mjd-1
	call M$MJD2YMD (mjd,ny,nm,nd)
	call M$YMD2MJD (ny,1,1,mjd_jan1)
	call M$YMD2MJD (1950,1,1,mjd1950_jan1)
	year = ny
	second = time - mjd_jan1*86400.d0
	ds50 = time/86400.d0 - mjd1950_jan1
	rag = RAGW(ds50)
	return
	end



      FUNCTION RAGW(DS50)
C     CAL. R. ASC. OF  GREENWICH (FROM NORAD ALGORITHM DEC. 1986)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      CR = 1.72027915249490159D-2
      THGR70 = 1.73213060363866D0
      TPI = 6.28318530717958648D0
      FK4R = 5.06494548928754D-15
C
      TS70 = DS50   -   7304.D0
      DS70 = IDINT( TS70 )
      RAGW = THGR70 + CR*DS70 + (CR+TPI)*(TS70-DS70) + TS70*TS70*FK4R
      RAGW = DMOD( RAGW, TPI)
      IF (RAGW.LT.-.5*TPI) RAGW = RAGW + TPI
      IF (RAGW.GT.+.5*TPI) RAGW = RAGW - TPI
      RETURN
      END

#if !_FDS
c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	subroutine M$OPEN_T6SUBRECORDS (hcb)

	implicit none

	include 'midas.inc'
	include 'headers.inc'
	include 'types.inc'
c Arguments
	record /HEADER/ hcb

c Locals
	record /T6SUBR_STR_STRUCT/ t6subr_str
	record /T6SUBR_STRUCT/     t6subr

	record /M_PTR/ t6recorddef_loc		! address pointers
	integer*4 startpos, n, junk, ls
	real*8    dval

	character*1 type
        character*2 fmt

c if not Type 6000, don't allocate any memory for subrecord info
	if (hcb.type/1000 .ne. 6) then
	  hcb.t6subr_alloc = 0
	  return
	endif

	hcb.t6subr_alloc = (hcb.subrecords*L_t6subr_len + 511)/512 * 512
	if (hcb.t6subr_alloc .eq. 0) then
	  hcb.t6subr_alloc = 512
	endif
	call M$STMALLOC (hcb.t6subr_alloc, hcb.t6subr_loc.p)

	type = 'A'
	ls = M$FIND_KEYDATA
     &	     (hcb, 'SUBREC_DEF', t6recorddef_loc.p, junk, 0 , type)

	if (ls .gt. 0) then
	  startpos = 0
	  do n = 1,hcb.subrecords
	    call M$MOVE
     &	    (%val(t6recorddef_loc.p+startpos),t6subr_str,L_t6subr_str_len)

	    t6subr.name		= t6subr_str.name
	    t6subr.format	= t6subr_str.format
	    ls			= M$A2D(t6subr_str.offset, dval)
	    t6subr.offset	= dval
	    ls			= M$A2D(t6subr_str.num_elts, dval)
	    t6subr.num_elts	= dval
	    ls			= M$A2D(t6subr_str.units, dval)
	    t6subr.units	= dval
	    ls			= M$A2D(t6subr_str.uprefix, dval)
	    t6subr.uprefix	= dval
	    ls			= M$A2D(t6subr_str.minval, t6subr.minval)
	    ls			= M$A2D(t6subr_str.maxval, t6subr.maxval)

	    call M$PUT_T6SUBRECORD (hcb, n, t6subr)

	    startpos = startpos + L_t6subr_str_len
	  enddo

	else
	  call M$WARNING ('Type 6000 subrecord definition keyword not found')
	endif

	return
	end

c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	subroutine M$UPDATE_T6SUBRECORDS (hcb)

	implicit none

	include 'midas.inc'
	include 'headers.inc'
	include 'types.inc'
c Arguments
	record /HEADER/ hcb

c Locals
	record /T6SUBR_STR_STRUCT/ t6subr_str
	record /T6SUBR_STRUCT/     t6subr

	integer*4 t6recorddef_alloc
	integer*4 startpos, n, junk, ls
	real*8    dval
	record /M_PTR/ t6recorddef_loc

	character*1 type

	if (hcb.type/1000 .ne. 6) then
	  return
	endif

	if (hcb.subrecords .gt. 0) then

	  t6recorddef_alloc = (hcb.subrecords*L_t6subr_str_len)
	  call M$STMALLOC (t6recorddef_alloc, t6recorddef_loc.p)

	  startpos = 0
	  do n = 1,hcb.subrecords
	    call M$GET_T6SUBRECORD (hcb, n, t6subr)

	    t6subr_str.name	= t6subr.name
	    t6subr_str.format	= t6subr.format
	    write (t6subr_str.offset,'(I8.8)')		t6subr.offset
	    write (t6subr_str.units,'(I4.4)')		t6subr.units
	    write (t6subr_str.uprefix,'(''+'',I2.2)')	t6subr.uprefix
	    write (t6subr_str.num_elts,'(I4.4)')	t6subr.num_elts
	    write (t6subr_str.minval,'(e24.17)')	t6subr.minval
	    write (t6subr_str.maxval,'(e24.17)')	t6subr.maxval
	    if (t6subr_str.minval(1:1) .NE. '-') t6subr_str.minval(1:1) = '+'
	    if (t6subr_str.maxval(1:1) .NE. '-') t6subr_str.maxval(1:1) = '+'

	    call M$MOVE
     &	    (t6subr_str,%val(t6recorddef_loc.p+startpos),L_t6subr_str_len)

	    startpos = startpos + L_t6subr_str_len
	  enddo

	  type = 'A'
	  ls = M$PUT_KEYDATA (hcb, 'SUBREC_DEF', %val(t6recorddef_loc.p),
     &         hcb.subrecords*L_t6subr_str_len, MK_replace, type)
	  call M$STMFREE (t6recorddef_alloc, t6recorddef_loc)

	endif

	return
	end

c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	subroutine M$PUT_T6SUBRECORD (hcb, n, t6subr)

	implicit none

	include 'midas.inc'
	include 'headers.inc'
c Arguments
	record /HEADER/ hcb
	record /T6SUBR_STRUCT/     t6subr
	integer*4 n

	call M$MOVE (t6subr.definition,
     &	             %val(hcb.t6subr_loc.p+L_t6subr_len*(n-1)),
     &	             L_t6subr_len)

c also, save corresponding info in Type 3000 subrecord fields
	if (n .le. L_subrecords) then
	  hcb.subr(n).name   = t6subr.name(1:4)
	  hcb.subr(n).format = t6subr.format
	  hcb.subr(n).offset = t6subr.offset
	endif

	return
	end

c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	subroutine M$GET_T6SUBRECORD (hcb, n, t6subr)

	implicit none

	include 'midas.inc'
	include 'headers.inc'
c Arguments
	record /HEADER/ hcb
	record /T6SUBR_STRUCT/     t6subr
	integer*4 n

	call M$MOVE (%val(hcb.t6subr_loc.p+L_t6subr_len*(n-1)),
     &	             t6subr.definition,
     &	             L_t6subr_len)

	return
	end

c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
#endif
