C+
C NAME:
C	MkShiftdn0
C PURPOSE:
C	To make XC shift maps at even heights to use in determining projections.
C	The maps give the average longitude shift in fractions of a rotation
C	required to get to the input reference velocity map at each latitude
C	and longitude point on the map. This version of the subroutine underwent 
C	considerable modification on about 4/1/02 when the way of determining ratios
C	were not accumulated with each height, but were accumulated relative to the base
C	at the source surface.  This version, which is where the majority of the 
C	time is spent in the tomography program, is more than twice as slow as the 
C	old version of the mkshiftd subroutine.
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call MkShiftdn0(XCbe,XCint,XCtbeg,XCtend,ALng,nLng,nLat,nMap,nT,nTmax
C		      nCar,JDCar,NCoff,VV,DD,XCshift,Vratio,Dratio,RR,dRR,CLRV,CLRD,VLL,VUL,
C		      NSIDE,Vtmp,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt)
C INPUTS:
C	XCbe(2,nT)	real		Boundary values of time maps
C	XCint(nTmax)	real		Boundaries of time intervals 
C	XCtbeg		real*8		Beginning time interval
C	XCtend		real*8		Ending time interval
C	Alng		real		Rotations per nLng1
C	nLng		integer		# longitude points
C	nLat		integer		# latitude points
C	nMap		integer		# heights (height begins at Sun__RAu AU)
C	nT		integer		# time intervals
C	nTmax		integer		Maximum # of time intervals
C	nCar		integer		# Carrington rotations
C	JDCar		real*8		Julian Date at beginning of rotations
C	NCoff		integer		# number of rotation offset values from the beginning rotation number
C	VV(nLng,nLat,nT)real		Map of velocity values at a given height RR
C	DD(nLng,nLat,nT)real		Map of density values at a given height RR
C	RR		real		Distance above sun for reference velocity map
C	dRR		real		Distance between individual maps (in AU)
C	CLRV		real		Velocity map radius filter distance
C	CLRD		real		Density map radius filter distance
C	VLL		real		Lower limit on traceback velocity
C	VUL		real		Upper limit on traceback velocity
C
C	Vtmp(nLng,nLat,nT,3)	real	Internal scratch array
C	Dtmp(nLng,nLat,nT)	real	Internal scratch array
C	XLTtmp(nLng,nT,5)	real	Internal scratch array
C	XLtmp(nLng,nT)		real	Internal scratch array
C	XLtmpt(nLng,nT)		real	Internal scratch array
C	Dtmpt(nT)		real	Internal scratch array
C OUTPUTS:
C	XCshift(nLng,nLat,nMap,nT,3) real XCshift contains the final shifts at all heights
C					in terms of fractions of a Carrington rotation
C	Vratio (nLng,nLat,nMap,nT) real	velocity ratios
C	Dratio (nLng,nLat,nMap,nT) real	density ratios
C FUNCTIONS/SUBROUTINES:
C PROCEDURE:
C >	XCshift(I,J,N,K,3) is the shift needed to trace location (I,J,N,K) back 
C	to the reference surface at a given time, i.e. 
C	XC = XCfrac(I)+XCshift(I,J,N,K,3) defines the origin of point (I,J,N,K,3)
C	at the reference surface in terms of a modified Carrington variable.
C	XCshift(I,J,N,K,3) will be negative above the reference surface (N>nRR), 
C	zero at the reference surface (N=nRR) and will point to an earlier
C	or the same time. Bad grid points in VV and DD are filled with 
C	GridSphere2D
C MODIFICATION HISTORY:
C	MAR-1997, B. Jackson (UCSD) (original name was MAKE_TS)
C	MAY-1998, Paul Hick (UCSD; pphick@ucsd.edu); revision
C	APR-1999, Paul Hick, B. Jackson (UCSD); revision
C	APR-2001, B. Jackson (UCSD); revision
C-
	subroutine MkShiftdn0(XCbe,XCint,XCtbeg,XCtend,ALng,nLng,nLat,nMap,nT,nTmax,
     &			nCar,JDCar,NCoff,VV,DD,XCshift,Vratio,Dratio,RR,dRR,FALLOFF,CLRV,CLRD,VLL,VUL,
     &			NSIDE,Vtmp,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt)

	real		VV(nLng,nLat,nT),		! Input velocity map at height RR
     &			DD(nLng,nLat,nT),		! Input density map at height RR
     &			XCshift(nLng,nLat,nMap,nT,3),	! Map of accumulated shifts at heights (from 1 RR)
     &			Vratio (nLng,nLat,nMap,nT),	! Map of accumulated velocity ratios at heights (from 1 RR)
     &			Dratio (nLng,nLat,nMap,nT),	! Map of accumulated density ratios at heights (from 1 RR)
     &			XCbe   (2,nT),			! Carrington time intervals
     &			XCint  (nTmax),			! Boundaries of time intervals
     &			XLT(2),				! Internal array

     &			Vtmp  (nLng,nLat,nT,3),		! Scratch arrays
     &			Dtmp  (nLng,nLat,nT),
     &			XLTtmp(nLng,nT,5),
     &			XLtmp (nLng,nT),
     &			XLtmpt(nLng,nT),
     & 			Dtmpt (nLng,nT)

	real*8		JDCar(nCar)
	
	real*8		XCtbeg,XCtend,XCtfull,X

	integer		NDLT(2)				! Internal array

	byte		NSIDE(*)			! (9*nLng*nLat)

	include		'MAPCOORDINATES.H'		! Contains executable
C			statements, so must follow last declaration statement

	if (nRR .lt. 1 .or. nRR .gt. nMap) call Say('MkShiftd','F','Illegal','reference distance')

	Bad = BadR4()
	VEL1AU = 425.0

C	print*, 'Into this mkshiftdn'

	NR13 = 0
	NR77 = 0
	KPVOLD  = 0

	NDLT(1) = nLng
	NDLT(2) = nT

	ANt  	  = sngl(XCtend-XCtbeg)			! No. Days over whole interval
	AKMPAUPD  = 1.495979e8/86400.			! Km*Day/AU
	AKMPAUPI  = AKMPAUPD*(float(nT-1)/ANt)		! Km*time_interval/AU
	
C	print *, 'ANt, float(nT-1), AKMPAUPI, Sun__Spiral ', ANt, float(nT-1), AKMPAUPI, Sun__Spiral, XCint(1)
	
C	ERot__Fract = Sun__Spiral*Sun__SiderealP/Sun__SynodicP	! 68.21347545*25.38/27.2753 
C								  (Rotates a little slower than inertial from Earth)
	ERot__Fract = AKMPAUPI/Sun__SiderealP		! What longitude has this material come from on the Sun?

C	print *, ' In MKshiftd - XCtbeg XCtend', XCtbeg, XCtend, ERot__Fract
C-------
C Store input arrays in Vtmp and Dtmp
	call ArrR4Copy2    ( nLngLatnt,VV,DD,Vtmp,Dtmp)	! Copy input arrays into scratch arrays
	call ArrR4SetMinMax(-nLngLatnt,Vtmp,VLL,VUL)	! Constrain Vtmp to [VLL,VUL] (exclude bad values)
C-------
C in the main program calls to this subroutine are preceeded by INTERP calls. Hence there are
C no bad values at this point. The GridSphere2D calls should be redundant.

       do N=1,nT

	call GridSphere2D(ALng,nLng,nLat,1,Vtmp(1,1,N,1),CLRV,4,0.0,0.0)! Fill in bad values
	call GridSphere2D(ALng,nLng,nLat,1,Dtmp(1,1,N),CLRD,4,0.0,0.0)
	call ArrR4Zero    (nLngLat,   XCshift(1,1,nRR,N,1))  ! Shifts at RR are 0
	call ArrR4Zero    (nLngLat,   XCshift(1,1,nRR,N,2))  ! Shifts at RR are 0
	call ArrR4Zero    (nLngLat,   XCshift(1,1,nRR,N,3))  ! Shifts at RR are 0
	call ArrR4Constant(nLngLat,1.,Vratio (1,1,nRR,N))    ! V ratios at RR are 1.
	call ArrR4Constant(nLngLat,1.,Dratio (1,1,nRR,N))    ! D ratios at RR are 1.
       end do

C-------
C The known velocities for the previous level (KM) are stored in Vtmp(*,*,*,KMV)
C The velocities to be computed for the current level (KP) are stored in Vtmp(*,*,*,KPV)

	KPV = 1
	do KP=nRR+1,nMap
	    KM = KP-1					! Height index of previous map

	    KMV = KPV
	    KPV = 3-KMV
	    
	    dXCL = ERot__Fract*(RRhght(KP)-RRhght(KM))	! Long. shift const. /(km/s)^(-1) from map KM to map KP
	    	    
	    dXCT = AKMPAUPI*(RRhght(KP)-RRhght(KM))	! Time shift const. per(km/s)^(-1) from map KM to map KP 
	    
	    Lmin = dXCL/VUL*(nLng1/ALng)		! Minimum longitude shift (based on maximum velocity)
	    Lmax = 1 + dXCL/VLL*(nLng1/ALng)		! Maximum longitude shift (based on minimum velocity)

	    Lmint = dXCT/VUL*((nT-1)/ANt)		! Minimum time shift (based on maximum velocity)
	    Lmaxt = 1 + dXCT/VLL*((nT-1)/ANt)		! Maximum time shift (based on minimum velocity)
	    
C	    print *, 'Lmin, Lmax, Lmint, Lmaxt =',Lmin, Lmax, Lmint, Lmaxt
C	    print *, 'dXCL, VUL, VLL, nLng1, ALng',dXCL, VUL, VLL, nLng1, ALng
C	    print *, 'dXCT, VUL, VLL, nT, ANt',   dXCT, VUL, VLL, nT, ANt

	    do J=1,nLat
		do N=1,nT

C		  dXCL = ERot__Fract*(RRhght(KP)-RRhght(KM))/(Sun__SynodicP/(XCint(N+1)-XCint(N)))
		  
C		  if(N.eq.2) print *, ' If ', dXCL, N, (XCint(N+1)-XCint(N))*(nT-1)/ANt, 1./Sun__SynodicP

		  XCbeg = XCbe(1,N)
		  XCend = XCbe(2,N)
		  do L=1,nLng
C
C  		    Longitude at longitude L at time N on level KM shifts to
C		    longitude XLtmp(L,N) at level KP
C
		    X          = (XCfull(L)+dXCL/Vtmp(L,J,N,KMV))		    
 		    XLtmp(L,N) = XCvar(X) 		! Shifted longitude (+ C.R.) on scale [1,nLng]
C
C		    Time at longitude L at time N on level KM shifts to time
C		    XLtmpt(L,N) at level KP
C
 		    XLtmpt(L,N) = XCtvar(XCtfull(N)+dXCT/Vtmp(L,J,N,KMV)) ! Shifted time (+) on scale [1,nT]
		    Dtmpt(L,N) = Dtmp(L,J,N)
		  end do
		end do

		do I=1,nLng				! Loop over all longitudes on level KP
		  do N=1,nT				! Loop over all times on level KP
		    VW  = 0.				! Accumulate longitude weighted velocities
		    W   = 0.				! Accumulate longitude weights
C
C  Calculate the weighted velocities and their weights at each longitude
C
		    LLmax = I+1+Lmax
		    if(LLmax.gt.nLng) LLmax = nLng
C		    do L=I+Lmin,min(nLng,I+1+Lmax)	! Search to the right of I on level KM at time N
		    do L=I+Lmin,LLmax			! Search to the right of I on level KM at time N

C		    	do M=N+Lmint,min(nT,N+1+Lmaxt)	! Search to lower times of N on level KM at longitude I
							! There was an error found here by BVJ (8/4/04)
			LTmin = N-1-Lmaxt
			if(LTmin.lt.1) LTmin = 1
			LTmax = N-Lmint
			if(LTmax.lt.1) LTmax = 1
		    	do M=LTmin,LTmax		! Search to lower times of N on level KM at longitude I

			    F = XLtmp(L,M)-I		! Distance between I and XLtmp(L,N)
			    F = abs(F)

			    if (F .lt. 1.) then		! Collect points within one grid spacing

				Ft = XLtmpt(L,M)-N	! Distance between N and XLtmpt(L,M)
				Ft = abs(Ft)

C	print *, 'L, M, F, XLtmp(L,M), I, Ft, XLtmpt(L,M), N', L, M, F, XLtmp(L,M), I, Ft, XLtmpt(L,M), N

				if (Ft .lt. 1.) then	! Collect points within one grid spacing
				    XCbeg = XCbe(1,M)
				    XCend = XCbe(2,M)
				    X      = (XCfull(L)+XCshift(L,J,KM,M,1))
				    XLT(1) = XCvar(X)
				    XLT(2) = XCtvar(XCtfull(M)+XCshift(L,J,KM,M,3))

C	print *, ' '
C	print *, ' Before FLINT - XLtmpt XLT(1) XLT(2) XCfull XCtFull XCshift
C     &	L J KM M', XLtmpt(L,M), XLT(1), XLT(2), XCfull(L), XCtfull(M), XCshift(L,J,KM,M,1),XCshift(L,J,KM,M,3),L, J, KM, M

			    	    XCV = FLINT(2,NDLT,Dtmpt,XLT,0.00002)

			    	    if (XCV .ne. Bad) then
					F = (1.-Ft)*(1.-F)*Dratio(L,J,KM,M)*XCV ! Weight with 1-F and with density at M on level KM
					F = max(0.,F)	! Redundant if input densities are all >=0

				    	VW = VW + F*Vtmp(L,J,M,KMV)
					W  = W  + F
				    else		! Density beyond time (or longitude) domain
					F = (1.-Ft)*(1.-F)*1.0
					VW = VW + F*Vtmp(L,J,M,KMV)
					W  = W  + F
				    end if
C	print *, ' After FLINTt - XCV F I J KM M VW W', XCV, F, I, J, KM, M, VW, W
				end if
			    end if
			end do
 		    end do
		    if (W .eq. 0) then
			Vtmp(I,J,N,KPV) = Bad
		    else
			Vtmp(I,J,N,KPV) = VW/W
		    end if
C	print *, 'V, I J N KPV', Vtmp(I,J,N,KPV), I, J, N, KPV
		  end do
		end do
	    end do
C-------
C Fill holes and smooth V maps using half the binwidth as window.
C Then add the new shift to the old shift closer to the reference surface.
C Amount of new shift needed to get to the last location is:
C
C	  print *, 'Through first loop of nMap', KP
C
C  Put in in 2000 (BVJ, UCSD/CASS)
C
	  NCM = 0
	  NCPP = 0
	  NCMM = 0
	  NTT  = nT
	  NTT2 = NTT/2
	  do N=NTT2,NTT-1
	    NNM = NTT - N
	    NNP = N   + 1
            call ArrR4getminmax(nLngLat,Vtmp(1,1,NNM,KPV),aminm,amaxm)
            call ArrR4getminmax(nLngLat,Vtmp(1,1,NNP,KPV),aminp,amaxp)
	    if(amaxp.eq.Bad.and.N.eq.NTT2) then
	      NCM = NCM + 1
C	      print *, 'Bad middle Vtmp in MkShift, Replace with constant velocity'
	      call ArrR4Constant(nLngLat,VEL1AU,Vtmp(1,1,NNP,KPV))
	    end if
	    if(amaxm.eq.Bad) then
	      NCMM = NCMM + 1
C	      print *, 'All bad values in Vtmp in MkShift NNM =',NNM,' at Ht. =',KP
	      call ArrR4Copy(nLngLat,Vtmp(1,1,NNM+1,KPV),Vtmp(1,1,NNM,KPV))
	    end if
	    if(amaxp.eq.Bad.and.N.ne.NTT2) then
	      NCPP = NCPP + 1	
C	      print *, 'All bad values in Vtmp in MkShift NNP =',NNP,' at Ht. =',KP
	      call ArrR4Copy(nLngLat,Vtmp(1,1,NNP-1,KPV),Vtmp(1,1,NNP,KPV))
	    end if
	  end do

	  if(NCM.gt.0)  print *, 'Bad middle Vtmp map in MkShift replaced with constant velocity'
	  if(NCMM.gt.0) print *, NCMM,' of',NTT2,' bad Vtmp - maps in MkShift at Ht. =',KP
	  if(NCPP.gt.0) print *, NCPP,' of',NTT2-1,' bad Vtmp + maps in MkShift at Ht. =',KP

	  do N=1,nt
	      call GridSphere2D(ALng,nLng,nLat,1,Vtmp(1,1,N,KPV),CLRV*2.0/3.0,4,0.0,0.0)	! Simply fill the rest of the holes
	  end do

	  NCM = 0
	  do N=1, nT
            call ArrR4getminmax(nLngLat,Vtmp(1,1,N,KPV),amin,amax)
	    if(amax.eq.Bad) NCM = 2
	  end do

	  if(NCM.eq.2) then
	      print *, 'Gridsphere in mkshift unable to fill all Vtmp holes. Attempt to fill with timesmooth'
	      call Timesmooth(nLng,nLat,nT,Vtmp(1,1,1,KPV),2,3.0,0.0,Dtmp)
	  end if

	  do J=1,nLat
	      do N=1,nt
		  do I=1,nLng
		       XLTtmp(I,N,1) = XCShift(I,J,KM,N,1)		! Longitude
		       XLTtmp(I,N,3) = XCShift(I,J,KM,N,3)		! Time
		       XLTtmp(I,N,2) = Vtmp   (I,J,N,KMV)		! Lower level velocity
		       XLTtmp(I,N,4) = Vratio (I,J,KM,N)		! Lower level Vratio
		       XLTtmp(I,N,5) = Dratio (I,J,KM,N)		! Lower level Dratio
		  end do
	      end do
	      do N=1,nT
C	          dXCL = AKMPAUPI*(XCint(N+1)-XCint(N))*(RRhght(KP)-RRhght(KP-1))*((nT-1)/ANt) ! Approximate this time

C		  dXCL = ERot__Fract*(RRhght(KP)-RRhght(KM))/(Sun__SynodicP*(XCint(N+1)-XCint(N))
C		  if(N.eq.2) print *, ' If ', dXCL, N, (XCint(N+1)-XCint(N))*(nT-1)/ANt, 1./Sun__SynodicP

		  XCbeg = XCbe(1,N)
		  XCend = XCbe(2,N)
	          do I=1,nLng					! Loop over locations at level KP
		      dXC0L = dXCL/Vtmp(I,J,N,KPV)		! Shift needed to go back to level KM in longitude
		      dXC0T = dXCT/Vtmp(I,J,N,KPV)		! Shift needed to go back to level KM in time

		      X      = (XCfull(I)-dXC0L)
		      XLT(1) = min(1.*nLng,XCvar(X)) 		! Longitude on level KM on scale [1,nLng]
								! Interpolate shifts at level KM

		      XLT(2) = max(1.,XCtvar(XCtfull(N)-dXC0T))	! Time on level KM on scale [1,nT]
								! Interpolate shifts at level KM

		      XCshift(I,J,KP,N,1) = FLINT(2,NDLT,XLTtmp(1,1,1),XLT,0.00002)-dXC0L ! Accumulated longitude shift
		      XCshift(I,J,KP,N,2) = 0.0						  ! Accumulated latitude shift
		      XCshift(I,J,KP,N,3) = FLINT(2,NDLT,XLTtmp(1,1,3),XLT,0.00002)-dXC0T ! Accumulated time shift

		      R = Vtmp(I,J,N,KPV)/FLINT(2,NDLT,XLTtmp(1,1,2),XLT,0.00002) ! Velocity ratio between this level and that
									!  projected from the base level at this height

		      Vratio(I,J,KP,N) =  FLINT(2,NDLT,XLTtmp(1,1,4),XLT,0.00002)*R ! Accumulated V ratio	     old way 9/17/03

		      Dratio(I,J,KP,N) =  FLINT(2,NDLT,XLTtmp(1,1,5),XLT,0.00002)/R ! Accumulated D ratio	     old way 9/17/03
		      		      
C		      if(I.ge.25.and.I.le.30.and.J.eq.5.and.N.eq.25) then		      
C		        print *, I,J,KP,N, dXCshiftL, dXCshiftT, Vratio(I,J,KP,N), Dratio(I,J,KP,N), Vtmp(I,J,N,KPV),
C     &				 FLINT(2,NDLT,XLTtmp(1,1,2),XLT,0.00002) 		      
C		      end if

		  end do
	      end do
	  end do

	end do

C	    call sim_MOD(XCbe,XCtbeg,XCtend,ALng,nLng,nLat,nMap,nT,
C     &			nCar,JDCar,NCoff,VV,DD,XCshift,Vratio,Dratio,RR,dRR,FALLOFF,CLRV,CLRD)

	return
	end