C+
C NAME:
C	MKVTRACE
C PURPOSE:
C	Calculate velocity traceback array from IPS velocity data
C CALLING SEQUENCE:
	subroutine MKVTRACE(JD,cVel,IPSD,NL,XR,YL,XC,NIND,VTrace,LngV,LatV,VSOL1,VSOL2,NSIDE)
C INPUTS:
C	JD		real*8
C	cVel		character*(*)
C	IPSD		integer
C	NL		integer
C	XR(NL)		real
C	YL(NL)		real
C	XC(NL)		real
C	NIND(NL)	integer
C	LngV		integer
C	LatV		integer
C	VSol1(LngV,LatV,2)real
C	VSol2(LngV,LatV,2)real
C	NSIDE(LngV,LatV,-1:1,-1:1)
C			integer*1
C OUTPUTS:
C	VTRACE(NL)
C CALLS:
C	Julian, N_Carrington, iSearch, ArrR4Copy, iFltArr, ArrR4Mask, GridFill
C	IndexR4, Say, BadR4
C EXTERNAL:
	external	fncWZ
C RESTRICTIONS:
C PROCEDURE:
C MODIFICATION HISTORY:
C	FEB-1996, Paul Hick (UCSD); extracted from main program SANHEL.FOR
C-
	    real*8	JD
	    character	cVel*(*)
	    integer	IPSD
	    integer	NL
	    real	XR(NL)
	    real	YL(NL)
	    real	XC(NL)
	    integer	NIND(NL)
	    real	VTrace(NL)
	    integer	LngV
	    integer	LatV
	    real	VSOL1(LngV,LatV,2)		! Wind velocity (km/s)
	    real	VSOL2(LngV,LatV,2)
	    integer*1	NSIDE(LngV,LatV,-1:1,-1:1)

	real*8		JEpoch
	character	cFile*80
	integer		nPad(4)		/4*0/

	if (NL .eq. 0) return

	Bad = BadR4()

	DXC = 1-LngV
	DYL = (LatV-1)/180.	! 1./Vertical grid spacing

	!-------
	! At the halfway point (JD) the Earth is in Carrington rotation NEarth.

	call Julian(1,iYr,Doy,JD,JEpoch)
	NEarth = N_CARRINGTON(iYr,Doy)

	!-------
	! 3D grids of UCSD IPS velocities are available for a number of Carrington
	! rotation. Find the grid for the Carrington rotation nearest to NEarth.

	J = 1
	I = 0
	write (cVel(8:8+3),'(I4.4)') NEarth+J*I
	do while (iSearch(1,cVel(:8+3)//'*.*',cFile) .ne. 1 .and. I .lt. 7)
	    J = -J
	    if (J .lt. 0) I = I+1
	    write (cVel(8:8+3),'(I4.4)') NEarth+J*I
	end do

	if (cFile .eq. ' ') call Say('MKVTRACE','E',cVel,'no matching UCSD IPS-velocity file')
	call Say('MKVTRACE','I',cFile,'is the UCSD IPS-velocity file')
	call Say('MKVTRACE','I','TRACE','Calculating traceback velocities')

	call IndexR4(1,NL,1,NL,XR,NIND)	! Sort (increasing 
					! heliocentric distance)
	nIPS = -1
	do I4=1,NL				! Loop over all points
	    J4 = NIND(I4)

	    if (XC(J4) .lt. 0) stop 'Negative XC'	! Safety belt

	    nIPSPrev = nIPS
	    nIPS = int(XR(J4)*100/IPSD)*IPSD
	    if (nIPS .ne. nIPSPrev) then		! Get new trace back velocity array

		!-------
		! Read IPS velocity arrays. Fill in missing grid values (i.e. velocity=0 read
		! from file) with GridFill

		if (nIPS .eq. nIPSPrev+IPSD) then
		    call ArrR4Copy(LngV*LatV*2,VSOL2,VSOL1)
		else
		    write (cVel(15:15+2),'(I3.3)') nIPS
		    cFile = cVel
		    I = iFltArr(cFile,140,0.,1.,fncWZ,Bad,LngV*LatV,LngV,LatV,nPad,VSOL1)
		    call GridFill(10,LngV,LatV,VSOL1,NSIDE,VMIN)
		end if

		write (cVel(15:15+2),'(I3.3)') nIPS+IPSD
		cFile = cVel
		I = iFltArr(cFile,140,0.,1.,fncWZ,Bad,LngV*LatV,LngV,LatV,nPad,VSOL2)
		call GridFill(10,LngV,LatV,VSOL2,NSIDE,VMIN)

	    end if

	    DIX = LngV+DXC*(XC(J4)-int(XC(J4)))
	    DIY = 1+DYL*(YL(J4)+90.)

	    IX  = min(LngV-1,int(DIX))	! Nearest grid point with smaller XC and YL
	    IY  = min(LatV-1,int(DIY))

	    DIX = DIX-IX		! Fractional distance to grid point IX,IY
	    DIY = DIY-IY
	    DIZ = (XR(J4)*100-nIPS)/IPSD

	    !-------
	    ! Safety belt

	    if (DIX .lt. 0 .or. DIY .lt. 0 .or. DIZ .lt. 0) stop 'Invalid fractional distance 1'
	    if (DIX .gt. 1 .or. DIY .gt. 1 .or. DIZ .gt. 1) stop 'Invalid fractional distance 2'

	    CIX = 1.-DIX
	    CIY = 1.-DIY
	    CIZ = 1.-DIZ

	    !-------
	    ! Linear interpolation over 8 neighbours (in XC,YL and VSOL).
	    ! The traceback velocity is calculated using the fully filled velocity
	    ! maps (as filled by GridFill). In an earlier version I tried using
	    ! only those data points for which a traceback velocity was available
	    ! from the unfilled maps. As a result so many points were rejected that
	    ! the statistics got significantly worse. Since the heliographic location
	    ! after traceback is fairly insensitive to the actual traceback velocity
	    ! I think the current approach is better.

	    VTrace(J4) =CIX*CIY*CIZ*VSOL1(IX,IY  ,1) + DIX*CIY*CIZ*VSOL1(IX+1,IY  ,1) +
     &			CIX*DIY*CIZ*VSOL1(IX,IY+1,1) + DIX*DIY*CIZ*VSOL1(IX+1,IY+1,1) +
     &			CIX*CIY*DIZ*VSOL2(IX,IY  ,1) + DIX*CIY*DIZ*VSOL2(IX+1,IY  ,1) +
     &			CIX*DIY*DIZ*VSOL2(IX,IY+1,1) + DIX*DIY*DIZ*VSOL2(IX+1,IY+1,1)

	end do

	return
	end