C+
C NAME:
C	xc3dtshift.f
C PURPOSE:
C	To shift the XC3DT matrix to incorporate a propagation upward from the lower magnetic surface.
C
C CATEGORY:
C	I/O
C CALLING SEQUENCE:
C	call xc3dtshift(nLng,nLat,nMap,nT,nTmaxG,RADMS,RR,dRR,XCBe,XCtbeg,XCtend,XC3DT)
C
C INPUTS:
C	nLng				integer		Number of Longitude bins	
C	nLat				integer		Number of Latitude bins
C	nMap				integer		Number of Map heights
C	nT				integer		Number of times in written out nson files (used in MAPCOORDINATES.H)		
C	nTmaxG				integer		Number of maximum time allowed in BB3
C	RADMS				real		Beginning Height of the Magnetic Field pick-up
C	RRS				real		Beginning Height of the 3D-MHD modeling
C	dRR				real		Outward step height of the modeling
C	XCbe(2,nTmaxG)			real		Beginning XCbe(1,1) and end XCbe(2,2) of Carrington rotation intervals 
C	XCtbeg				real*8		Start time
C	XCtend				real*8		End time
C
C OUTPUTS:
C	XC3DT(nLng,nLat,nMap,nTmaxG*4,3) real		XC3DT matrix to shift

C MODIFICATION HISTORY:
C	November-2015, Bernard Jackson (UCSD)
C-
	subroutine xc3dtshift(nLng,nLat,nMap,nT,nTmax,RADMS,RR,dRR,XCBe,XCtbeg,XCtend,XC3DT)

	real		XC3DT(nLng,nLat,nMap,nTmax*4,3)
	real		XCbe(2,nTmax)
	real		XCsh(3)
	real		XMAP(nLng,nLat,nTmax*4,3)
	real 		PP(3)								! Location at source surface
	integer 	ND(3)								! Array dimensions at source surface

	real*8		XCtbeg,XCtend,XCtfull,X,XCtim,XCt

	include 'MAPCOORDINATES.H'

	  print *, ' '
	  write(*,'(A,4I5)') 'Inside xc3dtshift, nLng, nLat, nMap, nT',nLng, nLat, nMap, nT

	  BAD = badR4()

	  XCbeg = XCbe(1,1)
	  XCend = XCbe(2,2)

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

	  noprint = 0

	  XBOTFRAC = (RR-RADMS)/dRR			! fraction of way from RADMS to RSS tomography surface relative to dRR

	  do I=1, nLng					! this will place zeros in the two bottom maps
	    do J=1, nLat
	      do N=1, nMap-1
	        do L=1, nT
		 do NX=1,3
		  XC3DT(I,J,nMap+1-N,L,NX) = XC3DT(I,J,nMap-N,L,NX)
		  if(noprint.eq.1) then
		    if(NX.eq.1.and.N.eq.5.and.J.eq.5.and.I.eq.30)
     & 		    write(*,'(A,5I4,2F10.4)'), 'In xc3dtshift.f-0',I,J,N,L,NX,XC3DT(I,J,nMap-N,L,NX),XC3DT(I,J,nMap+1-N,L,NX)
		    if(NX.eq.3.and.N.eq.5.and.J.eq.5.and.I.eq.30)
     & 		    write(*,'(A,5I4,2F10.4)'), 'In xc3dtshift.f-0',I,J,N,L,NX,XC3DT(I,J,nMap-N,L,NX),XC3DT(I,J,nMap+1-N,L,NX)

		    if(NX.eq.1.and.L.eq.20.and.J.eq.5.and.I.eq.30)
     & 		    write(*,'(A,5I4,2F10.4)'), 'In xc3dtshift.f-1',I,J,N,L,NX,XC3DT(I,J,nMap-N,L,NX),XC3DT(I,J,nMap+1-N,L,NX)
		    if(NX.eq.3.and.L.eq.20.and.J.eq.5.and.I.eq.30)
     & 		    write(*,'(A,5I4,2F10.4)'), 'In xc3dtshift.f-2',I,J,N,L,NX,XC3DT(I,J,nMap-N,L,NX),XC3DT(I,J,nMap+1-N,L,NX)
		   end if
		  end do
	        end do
	      end do
	    end do
	  end do
	  if(noprint.eq.1) print *, 'out 1'
	  do I=1, nLng					! Make an approximation that the second map from the bottom is the same as the 
	    do J=1, nLat				! third map from the bottom reduced by fraction XBOTFRAC and replace second row
	      do L=1, nT				! zeros by this
		do NX=1,3
		  XC3DT(I,J,2,L,NX) = XC3DT(I,J,3,L,NX)*XBOTFRAC
		  if(noprint.eq.1) then
		    if(NX.eq.1.and.L.eq.20.and.J.eq.5.and.I.eq.30)
     & 		    write(*,'(A,2I4,A,2I4,2F10.4)'), 'In xc3dtshift.f-3',I,J,' 2-3',L,NX,XC3DT(I,J,2,L,NX),XC3DT(I,J,3,L,NX)
		    if(NX.eq.3.and.L.eq.20.and.J.eq.5.and.I.eq.30)
     & 		    write(*,'(A,2I4,A,2I4,2F10.4)'), 'In xc3dtshift.f-4',I,J,' 2-3',L,NX,XC3DT(I,J,2,L,NX),XC3DT(I,J,3,L,NX)
		  end if
		  XMAP(I,J,L,NX) = XC3DT(I,J,2,L,NX)
	        end do
	      end do
	    end do
	  end do
	  if(noprint.eq.1) then
	    print *, 'out 2'
	    print *, ' '
	  end if
	  do I=1, nLng					! now modify all shift values at the third level and above
	    XClng = XCfull(I)				! Longitude the matrix is to be extracted
	    do J=1, nLat				! by the projected shift values on the second level
	      XLat = XLdeg(J)				! Latitude the matrix is to be extracted
	      do K=3, nMap

	        do L=1, nT
		  XCtim = XCtfull(L)

		  do NX=1,3
		    XCsh(NX) = XC3DT(I,J,K,L,NX)   

		  end do
		  XC = XClng + XCsh(1)				! To get the source surface projected center XCint spatial values
		  XC = max(XCbeg,min(XC,XCend))
		  XLa = XLat + XCsh(2)				! To get the source surface projected heliographic latitude values
		  XCt = XCtim + dble(XCsh(3))			! To get the source surface projected center XCint time values
		  XCt = max(XCtbeg,min(XCt,XCtend))

		  PP(1) = XCvar(XC)
		  PP(2) = XLindx(XLa)
		  PP(3) = XCtvar(XCt)
		  XCshift1 = FLINT(-3,ND,XMAP(1,1,1,1),PP,0.00002)
		  XCshift2 = FLINT(-3,ND,XMAP(1,1,1,2),PP,0.00002)
		  XCshift3 = FLINT(-3,ND,XMAP(1,1,1,3),PP,0.00002)

		  if(XCshift1.ne.BAD) then 
		    XXXX1 = XC3DT(I,J,K,L,1)
		    XC3DT(I,J,K,L,1) = XC3DT(I,J,K,L,1) + XCshift1
		  else
		    XC3DT(I,J,K,L,1) = BAD
		  end if
		  if(XCshift2.ne.BAD) then 
		    XXXX2 = XC3DT(I,J,K,L,2)
		    XC3DT(I,J,K,L,2) = XC3DT(I,J,K,L,2) + XCshift2
		  else
		    XC3DT(I,J,K,L,2) = BAD
		  end if
		  if(XCshift3.ne.BAD) then 
		    XXXX3 = XC3DT(I,J,K,L,3)
		    XC3DT(I,J,K,L,3) = XC3DT(I,J,K,L,3) + XCshift3
		  else
		    XC3DT(I,J,K,L,3) = BAD
		  end if
		  if(noprint.eq.1) then
		   do NX=1,3
		    if(NX.eq.1.and.K.eq.5.and.J.eq.5.and.I.eq.30)
     & 		    write(*,'(A,5I4,3F10.4)'), 'In xc3dtshift.f-5',I,J,K,L,NX,XC3DT(I,J,K,L,NX),XXXX1,XCshift1
		    if(NX.eq.2.and.K.eq.5.and.J.eq.5.and.I.eq.30)
     & 		    write(*,'(A,5I4,3F10.4)'), 'In xc3dtshift.f-6',I,J,K,L,NX,XC3DT(I,J,K,L,NX),XXXX2,XCshift2
		    if(NX.eq.3.and.K.eq.5.and.J.eq.5.and.I.eq.30)
     & 		    write(*,'(A,5I4,3F10.4)'), 'In xc3dtshift.f-7',I,J,K,L,NX,XC3DT(I,J,K,L,NX),XXXX3,XCshift3

		    if(NX.eq.1.and.L.eq.20.and.J.eq.5.and.I.eq.30)
     & 		    write(*,'(A,5I4,3F10.4)'), 'In xc3dtshift.f-8',I,J,K,L,NX,XC3DT(I,J,K,L,NX),XXXX1,XCshift1
		    if(NX.eq.2.and.L.eq.20.and.J.eq.5.and.I.eq.30)
     & 		    write(*,'(A,5I4,3F10.4)'), 'In xc3dtshift.f-9',I,J,K,L,NX,XC3DT(I,J,K,L,NX),XXXX2,XCshift2
		    if(NX.eq.3.and.L.eq.20.and.J.eq.5.and.I.eq.30)
     & 		    write(*,'(A,5I4,3F10.4)'), 'In xc3dtshift.f10',I,J,K,L,NX,XC3DT(I,J,K,L,NX),XXXX3,XCshift3
		    if(NX.eq.3.and.L.eq.20.and.J.eq.5.and.I.eq.30) print *, ' '
	          end do
		 end if
	        end do
	      end do
	    end do
	  end do
	  if(noprint.eq.1) then
	    print *, 'out 3'
	    print *, ' '
	  end if
	  return
	end