C+
C NAME:
C	SW_Model_Kinematic
C PURPOSE:
C	Creates a solar wind model from density and velocity map at the source
C	surface. Output is the 3D density and velocity and a 'shift array'
C	to be used in determination of source locations at the source surface.
C	The shift array gives for each grid point the average longitude shift in
C	fractions of a rotation required to get to the source surface.
C CATEGORY:
C	Tomography: solar wind model
C CALLING SEQUENCE:
	subroutine SW_Model_Kinematic(I__VD,XCbegMAT,XCendMAT,VV,G2,XC3D,V3D,D3D,XTtmp)
C INPUTS:
C	I__VD
C	XCbegMAT(nTim)			real		start Carrington variable of matrix (- NCoff)
C	XCendMAT(nTim)			real		end Carrington variable of matrix (- NCoff)
C							(range of matrix typically is three Carrington rotations)
C	VV(nLng,nLat,nTim)		real		Velocity map at at source surface RR
C	G2(nLng,nLat,nTim)		real		g^2 map at at source surface RR
C
C	XTtmp(nLng,nTim,2)		real		Scratch array
C OUTPUTS:
C	XC3D(3,nLng,nLat,nRad,nTim)	real		Shifts in terms of fractions of a Carrington rotation
C	V3D (nLng,nLat,nRad,nTim)	real		3D velocity matrix
C	G23D(nLng,nLat,nRad,nTim)	real		3D g^2 matrix
C CALLS:
C	T3D_get, T3D_iget, T3D_get_grid, Say, ArrR4Zero, ArrR4Copy, CvG2D, CvD2G
C	GridSphere2D, FLINT
C INCLUDE:
	include		'sun.h'
	include		't3d_param.h'
	include		't3d_array.h'
	include		't3d_index.h'

C	include		't3d_grid_fnc.h'
C	include		't3d_loc_fnc.h'
C RESTRICTIONS:
C >	If nTim = 1 (time-independent reconstruction; then dTT should be set to zero (dTT = 0).
C >	I think there is an implicit assumption that XCendMAT(i)-XCbegMAT(i) is the same
C	for all times i=1,nTim.
C >	There is NO CHECK for bad values (BadR4()) in VV and G2. Hence, the caller must make
C	sure that these two input array do not have any holes.
C PROCEDURE:
C >	ClipLng should be set to zero if iand(MODE,TOM__MOD) > 0, and to some non-zero
C	value (90 deg?) if not
C >	XC3D(I,J,N) is the shift needed to trace location (I,J,N) back to the reference
C	surface, i.e. XC = XCvar(I)+XC3D(I,J,N) defines the origin of point (I,J,N)
C	at the reference surface in terms of a modified Carrington variable.
C	XC3D(I,J,N) will be negative above the reference surface (N > 1) and zero
C	at the reference surface (N = 1).
C >	The input array VV and G2 should not have any holos (BadR4() values) in them
C >	G2 (g^2) is converted to normalized density right at the start. 3D density,
C	velocity and shift arrays are then build up stepping from one heliocentric distance
C	level to the next working from the source surface up to the outer boundary.
C	Before returning the density array is converted back to g^2.
C >	Once a level is finished holes at that level are filled with GridSphere2D before
C	continuing to the next level.
C MODIFICATION HISTORY:
C	MAR-1997, B. Jackson (UCSD/CASS) (original name was MAKE_TS)
C	MAY-1998, Paul Hick (UCSD/CASS); serious revision
C	AUG-2001, Paul Hick (UCSD/CASS; pphick@ucsd.edu); modified for time-dependent tomography
C-
	    integer	I__VD

	    real	XCbegMAT(*)
	    real	XCendMAT(*)

	    real	VV      (*)		! Input velocity map at height RR
	    real	G2      (*)		! Input density map at height RR
	    real	XC3D    (*)		! Accumulated shifts at heights (from 1 RR)
	    real	V3D     (*)		! Velocities
	    real	D3D     (*)		! G^2
	    real	XTtmp	(*)

	character	cSay*9		/'SW_Model_Kinematic'/
	real		VLL		/ 100./	! Only used to speed up things
	real		VUL		/1500./

	integer		ND(5)
	real		PP(5)

	logical		bMod
	logical		bNoMod
	logical		bMass
C	logical		bMHD

	real		Tiny		/0.0/

	include		't3d_grid_fnc.h'
	include		't3d_loc_fnc.h'

	locTmp(I,L,N) = I + nLng*( L-1 + nTim*(N-1) )


	call T3D_iget(T3D__MODE,0,MODE)

	if (I__VD .eq. TOM__V) call T3D_iget(T3D__NL_V, 0, NL)
	if (I__VD .eq. TOM__G) call T3D_iget(T3D__NL_G, 0, NL)

	if (I__VD .eq. TOM__V) call T3D_get(T3D__PWN_V, 0, PWN)
	if (I__VD .eq. TOM__G) call T3D_get(T3D__PWN_G, 0, PWN)

	call T3D_get_grid(TT,dTT,RR,dRR, nLng,nLat,nRad,nTim,dTTi,nLng1,nLat1)

	bMass	= iand(MODE,TOM__MASS) .ne. 0
	bMod	= iand(MODE,TOM__MOD ) .ne. 0
	bNoMod	= .not. bMod
	!bMHD	= iand(MODE,TOM__MHD ) .ne. 0

	Bad   = BadR4()

	ND(1) = PRJ__N
	ND(2) = nLng
	ND(3) = nLat
	ND(4) = nRad
	ND(5) = nTim

	PP(1) = float(PRJ__XC)
	Fill  = 90.0/nLat1

	!-------
	! Initialize source surface section in XC3D (set to zero), V3D, D3D
	! (copied from the source surface maps).

	NP = nLng*nLat
	do I=1,nTim
	    iloc = locSHFT(1,1,1,1,I) 
	    call ArrR4Zero(PRJ__N*NP ,XC3D(iloc))	! All shifts at RR are 0

	    jloc = locSS (1,1,I)
	    iloc = locVOL(1,1,1,I)
	    call ArrR4Copy(NP,VV(jloc),V3D(iloc))
	    call CvG2D(PWN,NP,G2(jloc),D3D(iloc))	! g^2 -> normalized density
	end do

	!-------
	! The 3D matrices are filled by an inductive process: given the values at
	! the lower level KM, find the values at the next higher level KP = KM+1.
	! Start out with the known values at the source surface KM = 1.

	do KP=2,nRad
	    KM	 = KP-1					! Radial index of previous map

	    dXC	 = SUN__SPIRAL*(RRhght(KP)-RRhght(KM))	! Shift (fraction of rotation) per (km/s)^(-1) from map KM to map KP

	    !-------
	    ! Imin, Imax, Lmin, Lmax are used only to reduce the range of do-loops.
	    ! The exact values are not important; just make sure to overestimate the range.

	    XCrange = XCendMAT(1)-XCbegMAT(1)

	    tmp  = dXC/XCrange*nLng1			! Shift (longitude grid spacings)
	    Imin = tmp/VUL				! Minimum shift (lng grid spacings) based on max velocity
	    Imax = tmp/VLL+1				! Maximum shift (lng grid spacings) based on min velocity


	    !!!!!!!!!!!  Don't think I trust these !!!

	    tmp	 = dXC*dTTi				! Shift (time grid spacings)
	    Lmin = tmp/VUL				! Minimum shift (time grid spacings) based on max velocity
	    Lmax = tmp/VLL+1				! Maximum shift (time grid spacings) based on min velocity

	    do J=1,nLat

		!-------
		! For each grid point at the lower level KM, calculate where and
		! when mass at that grid point arrives at level KP.

		do LM=1,nTim				! Loop over times on level KM
		    do IM=1,nLng			! Loop over longitudes on level KM
			tmp = dXC/V3D( locVOL(IM,J,KM,LM) ) ! Shift in Carrington variable from KM to KP

			!-------
			! Set up proper XCbeg, XCend, XCrange for time LM

			XCbeg	= XCbegMAT(LM)		! Needed for statement fnc XCvar
			XCend	= XCendMAT(LM)
			XCrange	= XCend-XCbeg

			XCM	= XCvar (IM)+tmp	! Carrington variable of IM traced to level KP
			TTM	= TTtime(LM)+tmp	! Time of IM traced to level KP

			!-------
			! Before converting XCM back to a longitude index, update XCbeg,
			! XCend, XCrange for time TTM.

			XCbeg   = FLINT(1,nTim,XCbegMAT,TTM,Tiny) ! Needed for statement fncs XCmod, XCindx
			XCend	= FLINT(1,nTim,XCendMAT,TTM,Tiny)
			XCrange	= XCend-XCbeg

			if (bMod) XCM = XCmod(XCM)	! Rescale inside [XCbeg,XCend]

			XCM	= XCindx(XCM)		! If ID=TOM__MOD not set these may be < 1
			TTM	= TTindx(TTM)

			XTtmp( locTmp(IM,LM,1) ) = XCM	! Longitude index on level KP on scale [1,nLng]
			XTtmp( locTmp(IM,LM,2) ) = TTM	! Time index on level KP

		    end do
		end do

		!-------
		! Now calculate densities and velocities at the upper level KP using the
		! information from the previous step (stored in XTtmp).

		do LP=1,nTim				! Loop over all times on level KP
		    do IP=1,nLng			! Loop over all longitudes on level KP

			VW = 0.				! Accumulate weighted velocities
			DW = 0.				! Accumulate weights (=normalized density)

			!-------
			! Only points on level KM at larger longitudes (smaller Carrington
			! variable) than IP (on level KP) and earlier times than LP can contribute.
			!
			! If ID=TOM__MOD then IMmax can be larger than nLng (they will be
			! scaled back into the valid index range [1,nLng]).
			! If ID=TOM__MOD not set then IMmax cannot be allowed to be bigger than nLng 

			IMmin = IP+  Imin
			IMmax = IP+1+Imax
			if (bNoMod) IMmax = min(IMmax,nLng)

			LMmin = max(1,LP-  Lmin)
			LMmax = max(1,LP-1-Lmax)

			do LM=LMmax,LMmin		! Search time earlier then LP on level KM
			    do IM0=IMmin,IMmax		! Search to the right of IP on level KM
				IM = IM0
				if (IM .gt. nLng) IM = IM0-nLng1! Make sure IM stays inside [1,nLng]
							! (redundant if ID=TOM__MOD not set)
				F = XTtmp( locTmp(IM,LM,1) )-IP	! Distance between IP and XTtmp(IM,LM,1)
				G = XTtmp( locTmp(IM,LM,2) )-LP	! Time between LP and XTtmp(IM,LM,2)
							! (for nTim=1 always G = 0.0)
				!-------
				! If ID=TOM__MOD then points +/- nLng1 need to be considered only.
				! These are points XCrange rotations away from XTtmp(IM), i.e. if
				! XCrange is an integer value these are location at the same
				! heliographic location but XCrange rotations earlier or later.

				if (bMod) then
				    F = min(abs(F-nLng1),abs(F),abs(F+nLng1))
				else
				    F = abs(F)
				end if

				!-------
				! Collect points within one grid spacing in longitude and time
				! (F < 1 and G < 1). Points within one grid spacing are summed with
				! weights (1-F)*(1-G). This effectively splits up each point into
				! four parcels assigned to four grid points at the upper level
				! KP with a total weight of one.

				if (F .lt. 1.0 .and. G .lt. 1.0) then
				    iloc = locVOL(IM,J,KM,LM)

				    tmp = D3D(iloc)	! Norm. density at IM on level KM

				    if (tmp .eq. Bad) call Say(cSay,'E','Bad','density; should not have happened')

				    F = (1.0-F)*(1.0-G)*tmp ! Weigh with (1-F)*(1-G) and with density at IM on level KM
				    !F = max(0.,F)	! Redundant if input densities are all >=0

				    VW = VW+F*V3D(iloc)	! Mass flux conservation (nvr^2 = constant)

				    !-------
				    ! SW_Model_DevFromR2 models the deviation of the density dropoff from a
				    ! a standard 1/r^2 behavior.

				    tmp = SW_Model_DevFromR2(D3D(iloc),V3D(iloc))
				    if (tmp .eq. 0.0) then
					DW = DW+F
				    else
					DW = DW+F*(RRhght(KM)/RRhght(KP))**tmp
				    end if

				end if

			    end do
			end do

			iloc = locVOL(IP,J,KP,LP)

			if (DW .eq. 0) then
			    V3D(iloc) = Bad
			    D3D(iloc) = Bad
			else
			    V3D(iloc) = VW/DW
			    D3D(iloc) = DW
			end if

		    end do
		end do

	    end do

	    !-------
	    ! Fill holes only using half the binwidth as window.

	    do LP=1,nTim
		XCrange = XCendMAT(LP)-XCbegMAT(LP)

		iloc = locVOL(1,1,KP,LP)
		call GridSphere2D(XCrange,nLng,nLat,1,V3D(iloc),Fill,4,0.0,ClipLng)
		call GridSphere2D(XCrange,nLng,nLat,1,D3D(iloc),Fill,4,0.0,ClipLng)
	    end do

	    !-------
	    ! Then add the new shift to the old shift closer to the reference surface.

	    PP(4) = float(KM)

	    do J=1,nLat					! Loop over locations at level KP
		PP(3) = float(J)

		do LP=1,nTim

		    do IP=1,nLng
							! Shift needed to go back to level KM
			tmp	= dXC/V3D( locVOL(IP,J,KP,LP) )

			XCbeg	= XCbegMAT(LP)		! Needed for statement fnc XCvar
			XCend	= XCendMAT(LP)
			XCrange	= XCend-XCbeg

			XCM	= XCvar (IP)-tmp	! Carrington var on level KM
			TTM	= TTtime(LP)-tmp	! Time on level KM

			XCbeg   = FLINT(1,nTim,XCbegMAT,TTM,Tiny) ! Needed for statement fncs XCmod, XCindx
			XCend	= FLINT(1,nTim,XCendMAT,TTM,Tiny)
			XCrange	= XCend-XCbeg

			if (bMod) XCM = XCmod(XCM)	! Rescale to range [XCbeg,XCend]

			XCM	= XCindx(XCM)		! Longitude index on level KM on scale [1,nLng]
			TTM	= TTindx(TTM)		! Time index on level KM

			!-------
			! TWO KLUDGES:
			!
			! TTM may become an index < 1. In addition, if iand(MODE,TOM__MOD) = 0
			! then XCM may be an index > nLng. If TTM < 0 or XCM > nLng+1 then FLINT
			! will return BadR4() (it extrapolates only over one gridspacing). Since
			! we really do want to fill XC3D with something we don't allow TTM or XCM
			! to go out of range.

			XCM	= min(XCM, nLng+1.0)
			TTM	= max(TTM, 0.0)

			PP(2)	= XCM
			PP(5)	= TTM

			!-------
			! The coordinates for 2nd (latitude) and 3rd (radial) dimension are
			! set to integers, so FLINT does a bi-linear interpolation over 1st
			! (longitude) and 4th (time) dimensions.

			F	= FLINT(5,ND,XC3D,PP,Tiny)! Interpolate longitude shifts at level KM

			!FF = FLINT(1,nLng,XC3D(1,J,KM,LP,PRJ__XC),XCM,Tiny)! Interpolate shifts at level KM
			!if (F .ne. FF) print *, 'SW_Model_Kinematic ',F,FF, F-FF

			if (F .eq. Bad) call Say(cSay,'E','Darn','FLINT returned bad shift')

			iloc = locSHFT(0,IP,J,KP,LP)

			XC3D(iloc+PRJ__XC) = F-tmp	! Accumulated longitude shift (always negative)
			XC3D(iloc+PRJ__XL) = 0.0	! Accumulated latitude shift
			XC3D(iloc+PRJ__TT) = F-tmp	! Accumulated time shift

		    end do
		end do

	    end do

	end do

	if (bMass) call CheckMass(XCbeg,XCend,V3D,D3D)

	call CvD2G(PWN,nLng*nLat*nRad*nTim,D3D,D3D)		! Norm. density -> g^2

	return
	end