C+
C NAME:
C	MkShiftdnmam_c   Must use to accomodate variable nT values
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 allows
C	the shifts to be determined from input volumetric data, and goes from
C	the top down. This is unlike the shift maker in the kinematic tomography
C	program that needs to "invent" conservation of mass and mass flux as the 
C	different height surfaces are stepped outward. This program assumes
C	an external program has run to provide the volumes that this program can use,
C	and that these volumes have conserved masss and mass flux, and are otherwise
C	complete with (if possible) accurate 3-component magnetic field data. 
C	Magnetic field data does not need to begin at the tomographic source surface 
C	height, which is constrained to be located at the lower level of the input 
C	volumetric data, but it must begin below that height.
C
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call MkShiftdnmam_c(iPrint,XCbe,XCtbeg,XCtend,ALng,nLng,nLat,nMap,nT,nTmax,
C		      VV3,BB3,DD1,TT1,XCshift3,XCshiftM,Vratio3,Bratio3,Dratio1,Tratio1,RR,RRMS,KPF,
C		      FALLOFFN,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,CLRV,CLRD,VLL,VUL,VLLLL,VLLUL,
C		      Vtmp,Dtmp,XLLTtmp)
C INPUTS:
C	iPrint		integer		if 1, print out lots of information
C	XCbe(2,nT)	real		Boundary values of time maps
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	VV3(nLng,nLat,NMap,nT,3) real	Maps of velocity component values at all heights and times
C	BB3(nLng,nLat,NMap,nT,3) real	Maps of magnetic field component values at all heights and times
C	DD1(nLng,nLat,NMap,nT)   real	Maps of density values at all heights and times
C	TT1(nLng,nLat,NMap,nT)   real	Maps of temperature values at all heights and times
C	RR		real		Distance above sun for reference velocity map
C	RRMS		real		Distance above sun for reference magnetic field map
C	KPF		integer		The height of the first tomographic inversion map. 
C				        If 1, the magnetic field map and the first traceback map have the same values.
C	FALLOFFN	real		Fall off in density
C	FALLOFFT	real		Fall off in temperature
C	FALLOFFBR	real		Fall off in radial field
C	FALLOFFBT	real		Fall off in tangential field
C	FALLOFFBN	real		Fall off in normal field
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	VLLLL		real		Lower limit on traceback longitude, latitude velocity
C	VLLUL		real		Upper limit on traceback longitude, latitude velocity
C
C	Vtmp(nLng,nLat,nT,3)	real	Internal scratch arrayif(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6)  
C	Dtmp(nLng,nLat,nT)	real	Internal scratch array
C	XLLTtmp(nLng,nT,8)	real	Internal scratch array
C OUTPUTS:
C	XCshift3(nLng,nLat,nMap,nT,3) real XCshift contains the final accumulated shifts (longitude at all heights
C					   in terms of fractions of a Carrington rotation), latitude, and time
C	XCshiftM(nLng,nLat,nT,3)      real	shift amount increments from RRMS to RR
C	Vratio3 (nLng,nLat,nMap,nT,3) real	velocity ratios
C	Bratio3 (nLng,nLat,nMap,nT,3) real	mag ratios
C	Dratio1 (nLng,nLat,nMap,nT)   real	density ratios
C	Tratio1 (nLng,nLat,nMap,nT)   real	temperature change
C	
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 MkShiftdnmam_c(iPrint,XCbe,XCtbeg,XCtend,ALng,nLng,nLat,nMap,nT,nTmax,
     &			VV3,BB3,DD1,TT1,XCshift3,XCshiftM,Vratio3,Bratio3,Dratio1,Tratio1,RR,RRMS,dRR,dRRMS,KPF,
     &			FALLOFFN,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,CLRV,CLRD,VLL,VUL,VLLLL,VLLUL,
     &			Vtmp,Dtmp,XLLTtmp)

	real		VV3	(nLng,nLat,NMap,nT,3),	! Input volumetric velocities
     &			BB3	(nLng,nLat,NMap,nT,3),	! Input volumetric magnetic fields
     &			DD1	(nLng,nLat,NMap,nT),		! Input volumetric densities
     &			TT1	(nLng,nLat,NMap,nT),		! Input volumetric temperatures
     &			XCshift3(nLng,nLat,nMap,nT,3),	! Map of accumulated shifts at heights (from 1 RR) (in long, lat, time)
     &			XCshiftM(nLng,nLat,nT,3),		! shift amount increments from RRMS to RR
     &			Vratio3 (nLng,nLat,nMap,nT,3),	! Map of accumulated velocity ratios at heights (from 1 RR)
     &			Bratio3 (nLng,nLat,nMap,nT,3),	! Map of accumulated magnetic field ratios at heights (from 1 RR) 
     &			Dratio1 (nLng,nLat,nMap,nT),		! Map of accumulated density ratios at heights (from 1 RR)
     &			Tratio1 (nLng,nLat,nMap,nT),		! Map of accumulated temperature ratios at heights (from 1 RR)
     &			XCbe   (2,nT),			! Carrington time intervals
     &			XLT (3),				! Internal array
     &			XLTB(3),				! Internal array
     &			BB33(3),				! Internal array
     &			VV33(3),				! Internal array			
     &			Vtmp    (nLng,nLat,nT,3),		! Scratch arrays
     &			Dtmp    (nLng,nLat,nT),
     &			XLLTtmp (nLng,nLat,nT,8)
	
	real*8		XCtbeg,XCtend,XCtfull,X

	integer		NDLT(3)					! Internal array

	include		'MAPCOORDINATES.H'			! Contains executable
C			statements, so must follow last declaration statement
C       iPrint = 0
C       iPrint = 1
	if(iPrint.eq.1) print *, ' '
	if(iPrint.eq.1) print*, 'Into mkshiftdnmam, nRR, iPrint', nRR, iPrint

	if(iPrint.eq.1) write (*,'(A)') 'XCtbeg,XCtend,ALng,nLng,nLat,nMap,nT,nTmax,RR,RRMS,dRR,dRRMS,'
	if(iPrint.eq.1) write (*,'(3F10.5,5I5,4F10.6)') XCtbeg,XCtend,ALng,nLng,nLat,nMap,nT,nTmax,RR,RRMS,dRR,dRRMS


        if(iPrint.eq.1) write (*,'(A)') 'KPF,FALLOFFN,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,CLRV,CLRD,VLL,VUL,VLLLL,VLLUL'
	if(iPrint.eq.1) write (*,'(I5,11F10.2)') KPF,FALLOFFN,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,CLRV,CLRD,VLL,
     &			VUL,VLLLL,VLLUL

	if(iPrint.eq.1) write (*,'(A)') 'XCbe'
	if(iPrint.eq.1) write (*,'(110F10.2)') XCbe

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

	Bad = BadR4()
C	VEL1AU = 425.0

	NR13 = 0
	NR77 = 0
	KPVOLD  = 0

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

	AKMPAUPD  = 1.495979e8/86400.			! Km/(AU*Day)
	AKMPAUPI  = AKMPAUPD				! (Km/(AU*Day))
C	The statement above was changed on 9/2/2011 BVJ
C
C	print *, 'float(nT-1), AKMPAUPI, Sun__Spiral ', float(nT-1), AKMPAUPI, Sun__Spiral
	
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 - NT-1, XCtbeg, XCtend, AKMPAUPI, ERot__Fract', NT-1, XCtbeg, XCtend,  
C     &			AKMPAUPI, ERot__Fract
C-------
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)\
C
C The above was the old way. Now the deconvolution must be done at the first level present 
C in the input model which is the location the inversion needs to take place. This is at
C the input model boundary, which must be defined as the source surface, and the traceback
C for the inversion needs to refer to this location in each higher level of the traceback 
C matrix. This is the first level. This is the place all vratios are
C determined to be 1.0. The magnetic field can come from a lower location, but this in handled 
C by an additional single element shift matrix called XCSHIFTM.
C
C Physically the source surface must start below all of the IPS lines of sight, which is defined 
C to be 11.6 degrees elongation, about what is needed for IPS sources at 327 MHz, is at 0.2 AU.
C
C This scheme presents a problem for magnetic field data that are extrapolated outward from 
C the solar surface since the way to extrapolate these magnetic fields have no way to have 
C their evolution determined below the inversion surface. The approximate solution to this 
C is used here is to have a second shift matrix to get down to the magnetic field source
C surface from the inversion surface by using traceback velocities from the lowest level
C of the modeled velocities. Magnetic fields are traced outward from the magnetic field
C source surface approximately using the lowest modeled velocities, and then from there given 
C a traceback in all of the higher levels to the inversion surface.
C
C 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.
C	
C The ratios on the first levels are only approximate
C
       if(iPrint.eq.1) print *, ' '

       do N=1,nT

C	call GridSphere2D(ALng,nLng,nLat,1,Vtmp(1,1,N,1),CLRV,4,0.0,0.0)! Fill in bad values
C	call GridSphere2D(ALng,nLng,nLat,1,Vtmp(1,1,N,2),CLRV,4,0.0,0.0)
C	call GridSphere2D(ALng,nLng,nLat,1,Vtmp(1,1,N,3),CLRV,4,0.0,0.0)
C	call GridSphere2D(ALng,nLng,nLat,1,Dtmp(1,1,N),CLRD,4,0.0,0.0)
	call ArrR4Zero    (nLngLat,  XCshift3(1,1,1,N,1))  ! Shifts at level 1 are 0
	call ArrR4Zero    (nLngLat,  XCshift3(1,1,1,N,2))  ! Shifts at level 1 are 0
	call ArrR4Zero    (nLngLat,  XCshift3(1,1,1,N,3))  ! Shifts at level 1 are 0
	call ArrR4Constant(nLngLat,1.,Vratio3(1,1,1,N,1))  ! V ratios at 1 are 1.
	call ArrR4Constant(nLngLat,1.,Vratio3(1,1,1,N,2))  ! V ratios at 1 are 1.
	call ArrR4Constant(nLngLat,1.,Vratio3(1,1,1,N,3))  ! V ratios at 1 are 1.
	call ArrR4Constant(nLngLat,1.,Bratio3(1,1,1,N,1))  ! V ratios at 1 are 1.
	call ArrR4Constant(nLngLat,1.,Bratio3(1,1,1,N,2))  ! V ratios at 1 are 1.
	call ArrR4Constant(nLngLat,1.,Bratio3(1,1,1,N,3))  ! V ratios at 1 are 1.
	call ArrR4Constant(nLngLat,1.,Dratio1(1,1,1,N))    ! D ratios at 1 are 1.
	call ArrR4Constant(nLngLat,1.,Tratio1(1,1,1,N))    ! T ratios at 1 are 1.
       end do
C
C 	In this do loop, deal with the shifts from each layer.
C
	if(iPrint.eq.1) write (*,'(A,2F15.5)') , 'ERot__Fract, AKMPAUPI',ERot__Fract, AKMPAUPI

	nlat1 = nLat-1
	nLng1 = nLng-1
	nMap1 = nMap-1
	nT1   = nT1 -1
	KP = 1
C	do K=1,nMap+1
	do K=1,nMap-1

	    KM  = K					! Height index of previous map
	    KP  = K + 1					! Height index of current map

	    if(iPrint.eq.1) print *, ' '
	    if(iPrint.eq.1) write (*,'(A,7I5)') , 'K, KP, KM, nMap, nT, nLat, nLng',K, KP, KM, nMap, nT, nLat, nLng
	    
	    do N=1,nT
	      do J=1,nLat
		do I=1,nLng
		  XLLTtmp(I,J,N,1) = XCShift3(I,J,KM,N,1)	! Longitude
		  XLLTtmp(I,J,N,2) = XCShift3(I,J,KM,N,2)	! Latitude
		  XLLTtmp(I,J,N,3) = XCShift3(I,J,KM,N,3)	! Time
		end do
	      end do
	      if(iPrint.eq.1) write(*,'(A,2I3,3F8.4))'), 
     &		'K,I=20,J=6,XLtp ',K,N,XLLTtmp(20,6,N,1),XLLTtmp(20,6,N,2),XLLTtmp(20,6,N,3) 
	    end do

C	      if(iPrint.eq.1) write(*,'(A,3F8.4))'), 
C     &		'I=20,J=6,N=22,XLtp ',K,N,XLLTtmp(20,6,22,1),XLLTtmp(20,6,22,2),XLLTtmp(20,6,22,3) 

	    dXCL   = ERot__Fract*(RRhght(KP)-RRhght(KM))  	! Long. shift const. /(km/s)^(-1) from map KM to map KP without VT
	    dXCT   = AKMPAUPI*(RRhght(KP)-RRhght(KM))		! Time shift const. per(km/s)^(-1) from map KM to map KP

C	    ALATA    = atan(dXCT/RRhght(KP))	    
	    if(iPrint.eq.1) 
     &		write (*,'(A,2I3,2F8.4,2F8.3)') , 'KP,KM,RRh(KP),RRh(KM),dXCL,dXCT         ',KP,KM,RRhght(KP),RRhght(KM),dXCL,dXCT	

	    if(iPrint.eq.1) write (*,'(A,7I5)') , 'K, KP, KM, nMap, nT, nLat, nLng',K, KP, KM, nMap, nT, nLat, nLng

	    do N=1,nT

	      do J=1,nLat
		ALAT  = 90.0*(-1.0+2.0*(J-1)/nLat1)				! Shift in longitude depends on latitude

		do I=1,nLng

		  if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6) write (*,'(A,2I3,2F8.4,2F8.3)'),
     &		 'KP,KM,RRh(KP),RRh(KM),ALAT (latitude)   ',KP,KM,RRhght(KP),RRhght(KM),ALAT

		  VV331 = VV3(I,J,N,KP,1)
		  VV332 = VV3(I,J,N,KP,2)
		  VV333 = VV3(I,J,N,KP,3)

		  dXCTTM = dXCT/VV331						! Time in days from one level to the next

		  if(VV332.eq.0.0) then						! Shift in Longitude for tangential velocity
		    dXCLOO = 0.0
		  else
		    if(ALAT.ne.-90.0.or.ALAT.ne.90.0) then
C	              ALONGA  = atan((dXCT/VV332)/RRhght(KP))/cos(ALAT)
		      ALONGA  = atan2d(dXCTTM*VV332/(AKMPAUPI*cos(ALAT)),RRhght(KP))	! Shift in longitude in degrees
		      ALONGA  = -ALONGA/360.0						! Shift in longitude in Carrington rotation
		    else
		      ALONGA = 0.0
		    end if

		    dXCLOO = ALONGA
		  end if

		  if(VV333.eq.0.0) then						! Shift in Latitude for normal (n) velocity in degrees
		    dXCLAA = 0.0
		  else
		    dXCLAA = atan2d(dXCTTM*VV333/AKMPAUPI,RRhght(KP))
		  end if

		  if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6) write (*,'(A,3F12.5)'),
     & 			'ERot__Fract, AKMPAUPI, DXCTTM     ',ERot__Fract, AKMPAUPI, dXCTTM
		  if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6) write (*,'(A,4I3,F8.4,F7.3,F9.3,F9.6)'), 
     &		        'KP,N,J,I,dXCT,VV332,VV333,ALONGA  ',KP,N,J,I,dXCT,VV332,VV333,ALONGA
		  if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6) write (*,'(A,4I3,2F7.3,2F9.6)'), 
     &			'KP,N,J,I,RRh(KP),RRh(KM),dXLO,dXLA',KP,N,J,I,RRhght(KP),RRhght(KM),dXCLOO,dXCLAA
		  if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6) print *, ' '

	    	  dXC0L  = dXCL/VV331	- dXCLOO				! XC Shift needed to go back to level KM in longitude
	    	  dXC0LA = dXCLAA						! Deg Shift needed to go back to level KM in latitude
	    	  dXC0T  = dXCT/VV331						! Time Shift needed to go back to level KM in time

		  if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I4,3F8.1))'), 
     &			 'KP,N,J,I,VV331,VV332,VV333',KP,N,J,I,VV331,VV332,VV333

		  if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I4,3F8.4))'), 
     &			 'KP, N, J, I, dL,  dLA, dT ',KP,N,J,I,dXC0L,dXC0LA,dXC0T

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

	    	  II = Itvar(XLT(3))					! Nearest interval time index
	    	  XCbeg = XCbe(1,II)
	    	  XCend = XCbe(2,II)

		  XLT(1)  = min(1.*nLng,XCvar(XCfull(I) - dXC0L)) 	! Carrington variable -> index [1,nLng]
									! Interpolate shifts at level KM

		  XLT(2) = XLindx(XLdeg(J)-dXC0LA)			! Latitude (deg) -> index in [1,nLat]
		  XLT(2) = min(float(nLat),XLT(2))			
		  XLT(2) = max(1.,XLT(2))
									! Interpolate shifts at level KM

		  if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I4,3F8.4))'), 
     &			 'KP,N,J,I, XCfu, XLde, XCtf',KP,N,J,I,XCfull(I),XLdeg(J),XCtfull(N)

		  if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I4,3F8.4))'), 
     &			 'KP,N,J,I, XLT1, XLT2, XLT3',KP,N,J,I,XLT(1),XLT(2),XLT(3)

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

C		  iPrint = 0
C		  if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6) then
		  if(iPrint.eq.1.and.KP.eq.10.and.I.eq.20.and.J.eq.6) then
			  print *, 'After Flint 1'
C			  write(*,'(A,4I4,3F8.4))'), 
C     &			 'KP, N, J, I, XL1, XL2, Xl3',KP,N,J,I,XLT(1),XLT(2),XLT(3)

			  write(*,'(A,4I4,3E9.2))'), 
     &			  'KP, N, J, I, FL1, FL2, FL3',KP,N,J,I,FLINT(3,NDLT,XLLTtmp(1,1,1,1),XLT,0.00002),
     &			  FLINT(3,NDLT,XLLTtmp(1,1,1,2),XLT,0.00002),FLINT(3,NDLT,XLLTtmp(1,1,1,3),XLT,0.00002)
		  end if
C		  iPrint = 0

		  if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6) then
			  print *, 'After Flint 2'
			  write(*,'(A,4I4,3E9.2))'), 
     &			  'KP, N, J, I, FL1, FL2, FL3',KP,N,J,I,FLINT(3,NDLT,XLLTtmp(1,1,1,1),XLT,0.00002),
     &			  FLINT(3,NDLT,XLLTtmp(1,1,1,2),XLT,0.00002),FLINT(3,NDLT,XLLTtmp(1,1,1,3),XLT,0.00002)
			  write(*,'(A,4I4,3F8.4))'), 
     &			  'KP, N, J, I, XS1, XS2, XS3',KP,N,J,I,XCshift3(I,J,KP,N,1),XCshift3(I,J,KP,N,2),XCshift3(I,J,KP,N,3)
			  print *, 'End T Loop'
		  end if

		end do
	      end do

C Maybe here use gridsphere to smooth each xcshift value
C	      call GridSphere2D(ALng,nLng,nLat,1,XCshift3(1,1,KP,N,1),CLRV*2.0/3.0,4,0.0,0.0) ! 4 says smooth and fill a hole if any exists
C	      call GridSphere2D(ALng,nLng,nLat,1,XCshift3(1,1,KP,N,2),CLRV*2.0/3.0,4,0.0,0.0)
C	      call GridSphere2D(ALng,nLng,nLat,1,XCshift3(1,1,KP,N,3),CLRV*2.0/3.0,4,0.0,0.0)

	    end do
	end do
C
C	******The shifts (XCshiftM) needed to get back to the magnetic surface are done seperately below.
C
	call xc3dtshift_rrms(iPrint,0,nLng,nLat,nMap,nT,RR,dRR,RRMS,XCshift3,XCshift3,XCshiftM)
C
C	*******End of XCshiftM creation for fields *********************
C
C
C 	In this do loop, deal with the ratios for each parameter in the input model. Ratios at K=1 are 1.0 except for magnetic field.
C
	if(iPrint.eq.1) print *, ' '
	if(iPrint.eq.1) print *, 'Now begin second loop to provide ratios'
	if(iPrint.eq.1) print *, ' '

	do N=1,nT
	      do J=1,nLat
		do I=1,nLng
		  XLLTtmp(I,J,N,1) = VV3(I,J,1,N,1)					! Base Velocity
		  XLLTtmp(I,J,N,2) = VV3(I,J,1,N,2)
		  XLLTtmp(I,J,N,3) = VV3(I,J,1,N,3)
		  XLLTtmp(I,J,N,4) = BB3(I,J,1,N,1)					! Base Magnetic Field for now
		  XLLTtmp(I,J,N,5) = BB3(I,J,1,N,2)
		  XLLTtmp(I,J,N,6) = BB3(I,J,1,N,3)
		  XLLTtmp(I,J,N,7) = DD1(I,J,1,N)					! Base Density
		  XLLTtmp(I,J,N,8) = TT1(I,J,1,N)					! Base Temperature
		end do
	      end do
	 end do

	 do KP=2,nMap
	    BFR = (RRhght(1)/RRhght(KP))**FALLOFFBR
	    BFT = (RRhght(1)/RRhght(KP))**FALLOFFBT
	    BFN = (RRhght(1)/RRhght(KP))**FALLOFFBN	    
	    DF  = (RRhght(1)/RRhght(KP))**FALLOFFN
	    TF  = (RRhght(1)/RRhght(KP))**FALLOFFT
	    do N=1,nT
	      do J=1,nLat
		do I=1,nLng
	    	  dXC0L  = XCshift3(I,J,KP,N,1)				! Shift needed to go back to level 1 in longitude
	    	  dXC0LA = XCshift3(I,J,KP,N,2)				! Shift needed to go back to level 1 in latitude
	    	  dXC0T  = XCshift3(I,J,KP,N,3)				! Shift needed to go back to level 1 in time

	    	  dXC0LB  = XCshift3(I,J,KP,N,1) + XCSHIFTM(I,J,N,1)	! Shift needed to go back to RRMS in longitude
	    	  dXC0LAB = XCshift3(I,J,KP,N,2) + XCSHIFTM(I,J,N,2)	! Shift needed to go back to RRMS in latitude
	    	  dXC0TB  = XCshift3(I,J,KP,N,3) + XCSHIFTM(I,J,N,3)	! Shift needed to go back to RRMS in time

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

		   XLTB(3) = max(1.,XCtvar(XCtfull(N) - dXC0TB))	! Time index on level KP on scale [1,nT] from RMMS
									! Interpolate shifts at level KP

	    	   II = Itvar(XLT(3))					! Nearest interval time index
	    	   XCbeg = XCbe(1,II)
	    	   XCend = XCbe(2,II)

		   XLT(1)  = min(1.*nLng,XCvar(XCfull(I) - dXC0L)) 	! Carrington variable -> index [1,nLng]
									! Interpolate shifts at level KP

		   II = Itvar(XLTB(3))					! Nearest interval time index from RRMS
	    	   XCbeg = XCbe(1,II)
	    	   XCend = XCbe(2,II)

		   XLTB(1)  = min(1.*nLng,XCvar(XCfull(I) - dXC0LB)) 	! Carrington variable -> index [1,nLng] from RRMS
									! Interpolate shifts at level KP

		   XLT(2) = XLindx(XLdeg(J)-dXC0LA)			! Latitude (deg) -> index in [1,nLat]
		   XLT(2) = min(float(nLat),XLT(2))			
		   XLT(2) = max(1.,XLT(2))

		   XLTB(2) = XLindx(XLdeg(J)-dXC0LAB)			! Latitude (deg) -> index in [1,nLat] from RMMS
		   XLTB(2) = min(float(nLat),XLTB(2))			
		   XLTB(2) = max(1.,XLTB(2))
									! Interpolate shifts at level KP
		iPrint = 0
		   if(iPrint.eq.1.and.KP.eq.10.and.I.eq.20.and.J.eq.6) then
			 write(*,'(A,4I4,3F12.4))'), 
C		   if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6) then
C			 write(*,'(A,4I4,3F8.4))'), 
     &			 'KP, N, J, I, XS1, XS2, XS3   ',KP,N,J,I,dXC0L,dXC0LA,dXC0T
C			 write(*,'(A,4I4,3F8.4))'), 
C     &			 'KP, N, J, I, XS1B, XS2B, XS3B',KP,N,J,I,dXC0LB,dXC0LAB,dXC0TB
			 write(*,'(A,4I4,3F12.5))'), 
C			 write(*,'(A,4I4,3F8.4))'), 
     &			 'KP, N, J, I, XL1, XL2, XL3   ',KP,N,J,I,XLT(1),XLT(2),XLT(3)
C			 write(*,'(A,4I4,3F8.4))'), 
C     &			 'KP, N, J, I, XL1B, XL2B, XL3B',KP,N,J,I,XLTB(1),XLTB(2),XLTB(3)
C			 write(*,'(A,4I4,2F8.2,F8.2))'), 
			 write(*,'(A,4I4,3F12.5))'), 
     &			 'KP, N, J, I, XTV1, XTV2, XTV3',KP,N,J,I,XLLTtmp(I,J,N,1),XLLTtmp(I,J,N,2),XLLTtmp(I,J,N,3)
C			 write(*,'(A,4I4,2F8.2,F8.2))'), 
C     &			 'KP, N, J, I, XTB1, XTB2, XTB3',KP,N,J,I,XLLTtmp(I,J,N,4),XLLTtmp(I,J,N,5),XLLTtmp(I,J,N,6)
			 write(*,'(A,4I4,2F12.5))'), 
     &			 'KP, N, J, I, XTD1, XTT2      ',KP,N,J,I,XLLTtmp(I,J,N,7)*RR*RR,XLLTtmp(I,J,N,8)/1000.
C			 write(*,'(A,4I4,F8.3,F12.1))'), 
C     &			 'KP, N, J, I, XTD1, XTT2      ',KP,N,J,I,XLLTtmp(I,J,N,7),XLLTtmp(I,J,N,8)
		   end if

		iPrint = 0
		   VV33(1) = FLINT(3,NDLT,XLLTtmp(1,1,1,1),XLT,0.00002)	   
		   VV33(2) = FLINT(3,NDLT,XLLTtmp(1,1,1,2),XLT,0.00002)	   
		   VV33(3) = FLINT(3,NDLT,XLLTtmp(1,1,1,3),XLT,0.00002)
		   DD11    = FLINT(3,NDLT,XLLTtmp(1,1,1,7),XLT,0.00002)	   
		   TT11    = FLINT(3,NDLT,XLLTtmp(1,1,1,8),XLT,0.00002)

		   BB33(1) = FLINT(3,NDLT,XLLTtmp(1,1,1,4),XLTB,0.00002)	   
		   BB33(2) = FLINT(3,NDLT,XLLTtmp(1,1,1,5),XLTB,0.00002)	   
		   BB33(3) = FLINT(3,NDLT,XLLTtmp(1,1,1,6),XLTB,0.00002)

C                  call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,XLLTtmp(1,1,1,4),3,	! Interpolated radial field from level 1
C     &                  	 dXC0Lb,dXC0LAB,dXC0TB,BB33)
		  if(BB33(1).gt.-0.1E+30.and.BB3(I,J,KP,N,1).ne.Bad) then
		    Bratio3(I,J,KP,N,1) = BB3(I,J,KP,N,1)/(BB33(1)*BFR)
		  else
		    Bratio3(I,J,KP,N,1) = 1.0
		  end if

C                  call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,XLLTtmp(1,1,1,5),3,	! Interpolated tangential field from level 1
C     &                  	 dXC0LB,dXC0LAB,dXC0TB,BB33)
		  if(BB33(2).gt.-0.1E+30.and.BB3(I,J,KP,N,2).ne.Bad) then
		    Bratio3(I,J,KP,N,2) = BB3(I,J,KP,N,2)/(BB33(2)*BFT)
		  else		iPrint = 1
		    Bratio3(I,J,KP,N,2) = 1.0
		  end if
C
C                  call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,XLLTtmp(1,1,1,6),3,	! Interpolated normal field from level 1
C     &                  	 dXC0LB,dXC0LAB,dXC0TB,BB33)
		  if(BB33(3).gt.-0.1E+30.and.BB3(I,J,KP,N,3).ne.Bad) then
		     Bratio3(I,J,KP,N,3) = BB3(I,J,KP,N,3)/(BB33(3)*BFN)
		  else
		    Bratio3(I,J,KP,N,3) = 1.0
		  end if
C                  call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,XLLTtmp(1,1,1,1),1,	! Interpolated radial velocity from level 1
C     &                  	 dXC0L,dXC0LA,dXC0T,VV33)
		  if(VV33(1).gt.0.0.and.VV3(I,J,KP,N,1).ne.Bad) then
		    Vratio3(I,J,KP,N,1) = VV3(I,J,KP,N,1)/VV33(1)
		  else
		    Vratio3(I,J,KP,N,1) = 1.0
		  end if
C                 call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,XLLTtmp(1,1,1,2),3,	! Interpolated tangential velocity from level 1
C     &                  	 dXC0L,dXC0LA,dXC0T,VV33)
		  if(VV33(2).gt.-0.1E+30.and.VV3(I,J,KP,N,2).ne.Bad) then
		    Vratio3(I,J,KP,N,2) = VV3(I,J,KP,N,2)/VV33(2)
		  else
		    Vratio3(I,J,KP,N,2) = 1.0
		  end if
C                  call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,XLLTtmp(1,1,1,3),3,	! Interpolated normal velocity from level 1
C     &                  	 dXC0L,dXC0LA,dXC0T,VV33)
		  if(VV33(3).gt.-0.1E+30.and.VV3(I,J,KP,N,3).ne.Bad) then
		    Vratio3(I,J,KP,N,3) = VV3(I,J,KP,N,3)/VV33(3)
		  else
		    Vratio3(I,J,KP,N,3) = 1.0
		  end if

		    if(iPrint.eq.1.and.I.eq.22.and.J.eq.5.and.N.eq.22) then
			print *, ' '
		    	write(*,'(A,4I3,3F9.3,F7.3)') 
     &              	'mkshift I,J,KP,N,BB3,BB12,BFR,Brat1',I, J, KP, N, BB3(I,J,KP,N,1),BB33(1),BFR,Bratio3(I,J,KP,N,1)
		    	write(*,'(A,4I3,3F9.3,F7.3)') 
     &              	'mkshift I,J,KP,N,BB3,BB12,BFT,Brat2',I, J, KP, N, BB3(I,J,KP,N,1),BB33(2),BFT,Bratio3(I,J,KP,N,2)
		    	write(*,'(A,4I3,3F9.3,F7.3)') 
     &              	'mkshift I,J,KP,N,BB3,BB13,BFN,Brat3',I, J, KP, N, BB3(I,J,KP,N,1),BB33(3),BFN,Bratio3(I,J,KP,N,3)
		    	write(*,'(A,4I3,2F9.3,9X,F7.3)') 
     &              	'In mkshiftd I,J,KP,N,VV3,VV31,Vrat1',I, J, KP, N, VV3(I,J,KP,N,1),VV33(1),vratio3(I,J,KP,N,1)
		    	write(*,'(A,4I3,2F9.3,9X,F7.3)') 
     &              	'In mkshiftd I,J,KP,N,VV3,VV32,Vrat2',I, J, KP, N, VV3(I,J,KP,N,1),VV33(2),vratio3(I,J,KP,N,2)
		    	write(*,'(A,4I3,2F9.3,9X,F7.3)') 
     &              	'In mkshiftd I,J,KP,N,VV3,VV33,Vrat3',I, J, KP, N, VV3(I,J,KP,N,1),VV33(3),vratio3(I,J,KP,N,3)
		    end if

C                  call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,XLLTtmp(1,1,1,7),1,	! Interpolated density from level 1
C     &                  	 dXC0L,dXC0LA,dXC0T,DD11)

		  if(DD11.gt.0.0.and.DD1(I,J,KP,N).ne.Bad) then
		    Dratio1(I,J,KP,N)   = DD1(I,J,KP,N)/(DD11*DF)
		iPrint = 0
C		    if(iPrint.eq.1.and.I.eq.22.and.J.eq.5.and.N.eq.22) then
C		    	write(*,'(A,4I3,3F9.3,F7.3)') 		    	 
		      if(iPrint.eq.1.and.KP.eq.10.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I3,4E9.2)'), 
     &              	'mkshiftd I,J,KP,N,DD1,DD11,DF,Dra',I, J, KP, N, DD1(I,J,KP,N),DD11,DF,Dratio1(I,J,KP,N)		      
C		    end if
		iPrint = 0
		  else
		    Dratio1(I,J,KP,N) = 1.0
		  end if
C                  call Get3DTvaCl(XCbe,XCtbeg,XCtend,nLng,nLat,nT,XLLTtmp(1,1,1,8),1,	! Interpolated temperature from level 1
C     &                  	 dXC0L,dXC0LA,dXC0T,TT11)
		  if(TT11.ne.0.0.or.TT1(I,J,KP,N).ne.Bad) then
		    Tratio1(I,J,KP,N)   = (TT1(I,J,KP,N)*TF)/TT11
		  else
		    Tratio1(I,J,KP,N) = 1.0
		  end if
C
		      if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I4,3F8.2))'), 
     &			 'KP, N, J, I, V1, V2, V3      ',KP,N,J,I,VV33(1),VV33(2),VV33(3)
		      if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I4,3F8.2))'), 
     &			 'KP, N, J, I, B1, B2, B3      ',KP,N,J,I,BB33(1),BB33(2),BB33(3)
		      if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I4,F8.3,F10.1))'), 
     &			 'KP, N, J, I, D, T            ',KP,N,J,I,DD11,TT11
		      if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I4,2F8.5))'), 
     &			 'KP, N, J, I, VRatio1, VRatio2',KP,N,J,I,Vratio3(I,J,KP,N,1),Vratio3(I,J,KP,N,2)
		iPrint = 0
C		      if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I4,2F8.5))'), 
C     &			 'KP, N, J, I, VRatio3, DRatio ',KP,N,J,I,Vratio3(I,J,KP,N,3),Dratio1(I,J,KP,N)
		      if(iPrint.eq.1.and.KP.eq.10.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I4,2F8.5)'), 
     &			 'KP, N, J, I, VRatio31, DRatio',KP,N,J,I,Vratio3(I,J,KP,N,1),Dratio1(I,J,KP,N)
C		      if(iPrint.eq.1.and.KP.eq.10.and.I.eq.20.and.J.eq.6)  print *, ' '
		iPrint = 0
		      if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I4,2F12.5))'), 
     &			 'KP, N, J, I, Bratio1, BRatio2',KP,N,J,I,Bratio3(I,J,KP,N,1),Bratio3(I,J,KP,N,2)
		      if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I4,2F8.5))'), 
     &			 'KP, N, J, I, Bratio3, TRatio ',KP,N,J,I,Bratio3(I,J,KP,N,3),Tratio1(I,J,KP,N)
		      if(iPrint.eq.1.and.N.eq.22.and.I.eq.20.and.J.eq.6) print *, ' '
		end do
	      end do
	    end do
	end do

	return
	end