C+
C NAME:
C	MassFind
C PURPOSE:
C	To find mass of specific heliospheric structures.
C CATEGORY:
C	Data analysis
C CALLING SEQUENCE:
C	call MassFind(XCbe,XCtbeg,XCtend,ALng,nLng,nLat,nMap,nT,
C		      VV,DD,XCshift,Vratio,Dratio,RR,dRR,CLRV,CLRD,CONVT,aNdayV,VLL,VUL,
C		      NSIDE,Vtmp,Dtmp,VZtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt)
C INPUTS:
C	XCbe(2,nT)	real		Boundary values of time maps
C	XCtbeg		real		Beginning time interval
C	XCtend		real		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	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	CONVT		real		Velocity temporal filter
C	aNdayV		real		Velocity numbers of days to combine temporally
C	VLL		real		Lower limit on traceback velocity
C	VUL		real		Upper limit on traceback velocity
C
C	Vtmp(nLng,nLat,nT,2)	real	Internal scratch array
C	Dtmp(nLng,nLat,nT)	real	Internal scratch array
C	VZtmp(nLng,nLat,nT)	real	Internal scratch array
C	XLTtmp(nLng,nT,4)	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) 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) 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) defines the origin of point (I,J,N,K)
C	at the reference surface in terms of a modified Carrington variable.
C	XCshift(I,J,N,K) 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-
	subroutine MkShiftd(XCbe,XCtbeg,XCtend,ALng,nLng,nLat,nMap,nT,
     &			VV,DD,XCshift,Vratio,Dratio,RR,dRR,CLRV,CLRD,CONVT,aNdayV,VLL,VUL,
     &			NSIDE,Vtmp,Dtmp,VZtmp,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),	! 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 
     &			XLT(2),				! Internal array

     &			Vtmp  (nLng,nLat,nT,3),		! Scratch arrays
     &			Dtmp  (nLng,nLat,nT),
     &			VZtmp (nLng,nLat,nT),		! Added 7/30/01 for Timesmooth routine
     &			XLTtmp(nLng,nT,4),
     &			XLtmp (nLng,nT),
     &			XLtmpt(nLng,nT),
     & 			Dtmpt (nLng,nT)

	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

	NR13 = 0
	NR77 = 0
	KPVOLD  = 0

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

	ANt  = XCtend-XCtbeg

C	print *, ' In MKshiftd - XCtbeg XCtend', XCtbeg, XCtend

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

C          call arrR4getminmax(nLngLat,Dtmp(1,1,N),amin,amax)
C          print *, 'DMAP min max',amin,amax
C          call arrR4getminmax(nLngLat,Vtmp(1,1,N,1),amin,amax)
C          print *, 'VMAP 1 min max',amin,amax
C          call arrR4getminmax(nLngLat,Vtmp(1,1,N,2),amin,amax)
C          print *, 'VMAP 2 min max',amin,amax

	call GridSphere2D(ALng,nLng,nLat,1,Vtmp(1,1,N,1),CLRV,4,0.0,0.0)! Fill in bad values
	call ArrR4Copy   (nLngLat,Vtmp(1,1,N,1),Vtmp(1,1,N,3))
	call GridSphere2D(ALng,nLng,nLat,1,Dtmp(1,1,N),CLRD,4,0.0,0.0)

	call ArrR4Zero    (nLngLat,   XCshift(1,1,nRR,N))  ! 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

	    dXC  = Sun__Spiral*(RRhght(KP)-RRhght(KM))	! Shift per (km/s)^(-1) from map KM to map KP
	    Lmin = dXC/VUL*(nLng1/ALng)			! Minimum longitude shift (based on maximum velocity)
	    Lmax = 1 + dXC/VLL*(nLng1/ALng)		! Maximum longitude shift (based on minimum velocity)

	    Lmint = dXC/VUL*((nT-1)/ANt)		! Minimum time shift (based on maximum velocity)
	    Lmaxt = 1 + dXC/VLL*((nT-1)/ANt)		! Maximum time shift (based on minimum velocity)

C	    print *, 'Lmin, Lmax, Lmint, Lmaxt =',Lmin, Lmax, Lmint, Lmaxt

	    do J=1,nLat
		do N=1,nt	
		  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
 		    XLtmp(L,N) = XCvar(XCfull(L)+dXC/Vtmp(L,J,N,KMV)) ! Shifted longitude 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)+dXC/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
		    do L=I+Lmin,min(nLng,I+1+Lmax)	! Search to the right of I on level KM at time N

		    	do M=N+Lmint,min(nT,N+1+Lmaxt)! 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(I,M)
				Ft = abs(Ft)

				if (Ft .lt. 1.) then	! Collect points within one grid spacing		    
				    XCbeg = XCbe(1,M)
				    XCend = XCbe(2,M)
				    XLT(1) = XCvar(XCfull(L)+XCshift(L,J,KM,M))
				    
				    XLT(2) = XCtvar(XCtfull(M)+XCshift(L,J,KM,M))
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), L, J, KM, M
			    	    XCV = FLINT(2,NDLT,Dtmpt,XLT)
			    	    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	  print *, 'Through first loop of nMap', KP

C
C  Put in in 2000 (BVJ, UCSD/CASS)
C
	  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
	      print *, 'Bad middle Vtmp in MkShift, Replace with constant velocity'
	      call ArrR4Constant(nLngLat,VEL1AU,Vtmp(1,1,NNM+1,KPV))
	    end if
	    if(amaxm.eq.Bad) then
	      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) then
	      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

	  PRS  = 2.0
	  FACT = 4.0
	  RSMOOTH = FACT*(RRhght(KP)**PRS)

C	  print *, 'KP, RRhght(KP), FACT, PRS, RSMOOTH =', KP, RRhght(KP), FACT, PRS, RSMOOTH

	  aNgs1S    = 0.0
	  Vtmptgs1S = 0.0
	  aNgs2S    = 0.0
	  Vtmptgs2S = 0.0
	  do N=1,nt

	  aNgs1 = 0.0
	  Vtmptgs1 = 0.0
	  do II=1,nLng
	    do JJ=1,nLat
		if(Vtmp(II,JJ,N,KPV).ne.BadR4()) then
		  aNgs1 = aNgs1 + 1.0
		  Vtmptgs1 = Vtmp(II,JJ,N,KPV) + Vtmptgs1
		end if
	    end do
	  end do
	  aNgs1S    = aNgs1 + aNgs1S
	  Vtmptgs1S = Vtmptgs1 + Vtmptgs1S
	  Vtmptgs1 = Vtmptgs1/aNgs1

C          call arrR4getminmax(nLngLat,Vtmp(1,1,N,KPV),amin,amax)
C          print *, 'Vtmp min max N KP',amin,amax,N,KP

	      call GridSphere2D(ALng,nLng,nLat,1,Vtmp(1,1,N,KPV),90./(nLat-1),4,0.0,90.0)	! Simply fill the holes

C	      call GridSphere2D(ALng,nLng,nLat,1,Vtmp(1,1,N,KPV),CLRV*RSMOOTH,3,0.0,90.0) 	! 5/3/01 Smoothing now used
												! 8/3/01 Now with CLRV*RSMOOTH
	  aNgs2 = 0.0
	  Vtmptgs2 = 0.0
	  do II=1,nLng
	    do JJ=1,nLat
		if(Vtmp(II,JJ,N,KPV).ne.BadR4()) then
		  aNgs2 = aNgs2 + 1.0
		  Vtmptgs2 = Vtmp(II,JJ,N,KPV) + Vtmptgs2
		end if
	    end do
	  end do

	  aNgs2S    = aNgs2 + aNgs2S
	  Vtmptgs2S = Vtmptgs2 + Vtmptgs2S
	  Vtmptgs2  = Vtmptgs2/aNgs2
	  
C	    call arrr4getminmax(nLng*nLat,Vtmp(1,1,N,KPV),amin,amax)
C	    print *, 'MkShift Vtmp BadR4? KPV KP', amin,amax,KPV,KP
C	    call arrr4getminmax(-nLng*nLat,Vtmp(1,1,N,KPV),amin,amax)
C	    print *, 'MkShift Vtmp 0? KPV KP', amin,amax,KPV,KP

C	    call arrr4getminmax(nLng*nLat,Vtmp(1,1,N,KMV),amin,amax)
C	    print *, 'MkShift Vtmp BadR4? KMV KP', amin,amax,KMV,KP
C	    call arrr4getminmax(-nLng*nLat,Vtmp(1,1,N,KMV),amin,amax)
C	    print *, 'MkShift Vtmp 0? KMV KP', amin,amax,KMV,KP

	  end do

	  Vtmptgs1S = Vtmptgs1S/aNgs1S
	  Vtmptgs2S = Vtmptgs2S/aNgs2S

	  RtmptgsS = Vtmptgs1S/Vtmptgs2S

C	  do II=1,nLng
C	    do JJ=1,nLat
C	      do NN=1,nt
C		if(Vtmp(II,JJ,NN,KPV).ne.BadR4()) Vtmp(II,JJ,NN,KPV) = Vtmp(II,JJ,NN,KPV)*RtmptgsS
C	      end do
C	    end do
C	  end do

	  aNs1 = 0.0
	  Vtmpts1 = 0.0
	  do II=1,nLng
	    do JJ=1,nLat
	      do NN=1,nLat
		if(Vtmp(II,JJ,NN,KPV).ne.BadR4()) then
		  aNs1 = aNs1 + 1.0
		  Vtmpts1 = Vtmp(II,JJ,NN,KPV) + Vtmpts1
		end if
	      end do
	    end do
	  end do
	  Vtmpts1 = Vtmpts1/aNs1

C	  call Timesmooth(nLng,nLat,nT,Vtmp(1,1,1,KPV),0,RSMOOTH*CONVT/aNdayV,0.,VZtmp)		! Smoothing put in 5/4/01
												! 8/3/01 now with RSMOOTH*CONVT

	  aNs2 = 0.0
	  Vtmpts2 = 0.0
	  do II=1,nLng
	    do JJ=1,nLat
	      do NN=1,nLat
		if(Vtmp(II,JJ,NN,KPV).ne.BadR4()) then
		  aNs2 = aNs2 + 1.0
		  Vtmpts2 = Vtmp(II,JJ,NN,KPV) + Vtmpts2
		end if
	      end do
	    end do
	  end do
	  Vtmpts2 = Vtmpts2/aNs2
	  Rtmpts = Vtmpts1/Vtmpts2

C	  do II=1,nLng
C	    do JJ=1,nLat
C	      do NN=1,nt
C		if(Vtmp(II,JJ,NN,KPV).ne.BadR4()) Vtmp(II,JJ,NN,KPV) = Vtmp(II,JJ,NN,KPV)*Rtmpts
C	      end do
C	    end do
C	  end do

C	  dXCRR = Sun__Spiral*(RRhght(KP)-RRhght(1))

	  do J=1,nLat
	      do N=1,nt
		  do I=1,nLng
		       XLTtmp(I,N,1) = XCShift(I,J,KM,N)
		       XLTtmp(I,N,2) = Vtmp   (I,J,N, KMV)
		       XLTtmp(I,N,3) = Vratio (I,J,KM,N)
		       XLTtmp(I,N,4) = Dratio (I,J,KM,N)
		       XLtmp (I,N)   = Vtmp   (
I,J,N,3)
		  end do
	      end do
	      do N=1,nT
		  XCbeg = XCbe(1,N)
		  XCend = XCbe(2,N)
	          do I=1,nLng					! Loop over locations at level KP
		      dXC0 = dXC/Vtmp(I,J,N,KPV)		! Shift needed to go back to level KM

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

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

C		      if(J.eq.1) print *, ' Begin'
		      XCshift(I,J,KP,N) = FLINT(2,NDLT,XLTtmp(1,1,1),XLT)-dXC0 ! Accumulated shift

	call arrr4getminmax(nLng*nLat,Vtmp(1,1,N,KPV),amin,amax)
	call arrr4getminmax(nLng*nLat,XLTtmp(1,1,2),amin1,amax1)

C	if(J.eq.5.and.N.eq.20.and.KP.eq.11) print *, ' Begin 2, N XLT(1),XLT(2)',N,NDLT(1),NDLT(2),XLT(1),XLT(2),R,Vtmp(I,J,N,KPV),amin,amax,amin1,amax1

		      dXCshift = XCshift(I,J,KP,N)			! Shift needed to go back to base

		      XLT(1) = min(1.*nLng,XCvar(XCfull(I)+dXCshift)) 	! Longitude on base level on scale [1,nLng]

		      XLT(2) = max(1.,XCtvar(XCtfull(N)+dXCshift))	! Time on base level on scale [1,nT]

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

C		      Vforvratio = -dXCRR/dXCshift

C		      R = Vforvratio/FLINT(2,NDLT,XLtmp,XLT) 	! Velocity ratio between this level from the total shift
								!    and that projected from the base level at this height.

C		      R = R**(1.0/(1.0 + 4.0*RRhght(KP)))

C		      R = 1.0

C	if(J.eq.5.and.N.eq.20.and.KP.eq.11) print *, ' Begin 3, I, J, N, KPV, KP, R', I,J,N,KPV,KP,R

		      Vratio(I,J,KP,N) =  R 				! Accumulated V ratio
 on this level

C		      if(KP.ne. KPVOLD) then
C			print *, 'The number of R greater than 1.3 this ', KP,' level is: ', NR13 
C			NR13 = 0
C		      end if
C		      if(R.ge.1.3)  then
C			NR13 = NR13 + 1
C			R = 1.3
C		      end if
C		      if(KP.ne. KPVOLD) then
C			print *, 'The number of R less than 0.77 this ', KP,' level is: ', NR77 
C			NR77 = 0
C		      end if
C		      If(R.le.0.77) then
C			NR77 = NR77 + 1
C			R = 0.77
C		      end if
C			KPVOLD = KP

C	if(J.eq.5.and.N.eq.20.and.KP.eq.11) print *, ' R, Vratio, Vtmp(I,J,N,KPV),amin,amax,amin1,amax1 =', R, Vratio(I,J,KP,N)
, Vtmp(I,J,N,KPV),amin,amax,amin1,amax1

		      Dratio(I,J,KP,N) =  1.0/R 			! Accumulated D ratio
 this level if mv = constant

C	if(J.eq.5.and.N.eq.20.and.KP.eq.11) print *, ' End 4, Dratio, Vtmpts1/Vtmpts2 =', Dratio(I,J,KP,N)
, Vtmpts1/Vtmpts2, Rtmpts

		  end do
	      end do
C	      print *, ' End nLng, nT latitude J =',J
	  end do

	end do

	return
	end