C+
C NAME:
C	external_fill_vol
C PURPOSE:
C	To fill the tomography with externally provided volumes of data over times of the volumetric data set.
C	This program also reads the input volumetric data file that has been prepared and converted for use in the current 
C	UCSD tomography. These input files are in heliographic coordinates, and at evenly spaced heights (usually 0.1 AU) upward
C	from the source surface. Input are 3 velocity components, 3 magnetic field components, density, and temperature.  
C	The volumetric input data can be between limits that may be extended to those required by the kinematic tomography in time; 
C 	the kinematic tomography can provide an extension it time on either end of the input volumetric data.
C	This program also allows written volumetric output of the type expected from the input volumetric data.
C
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	external_fill_vol(    )
C INPUTS:
C	nTTVa		   integer	0  - cycle through all values of nT
C					>0 - use nTTVa as the one and only value of nT to correct
C	XCbe		   real		Beginning (and ending) rotation values
C	nLng		   integer	# Longitude points
C	nLat		   integer	# Latitude points
C	nT		   integer	# DMap times
C	CONST		   real		Time filter constant (in terms of the time interval)
C	DMap(nLng,nLat,nT) real		Input map
C OUTPUTS:
C	DMap(nLng,nLat,nT) real		Output map
C SCRATCH ARRAY:
C	DMp(nLng)	   real
C FUNCTIONS/SUBROUTINES:
C PROCEDURE:
C	Bad values (indicated by BadR4()) are not processed
C MODIFICATION HISTORY:
C	FEB, 2016 B. Jackson (UCSD)
C-
	subroutine external_fill_vol(nTTVa,XCbe,nLng,nLat,nT,CONST,DMap,DMp)

	real		DMap(nLng,nLat,nT)	! Map
	real		XCbe(2,nT)		! Map longitude index

	real		DMp(nLng)
	real		dRR	/.1/
	
	real*8		XCtbeg,XCtend,XCtfull,X		! Probably not needed since XCTfull is never called
C
	include		'MAPCOORDINATES.H'
C
C  If no common time values are present, then average the interpolated
C  non-Bad time values.
C
	Bad = BadR4()
	    write(*,'(A,2I4)') 'Before filling test arrays; time LLbeg, LLend         ',LLbeg,LLend
	    print *, 'Fill array D & V NOT between LLbeg and LLend from the kinematic tomography'
	    print *, ' '

C 	Fill the portion not filled from an external source by normal means from DMAP and VMAP (needed anyway)

	    do N=1,nTG
	      do K=1,nMap
	        do J=1,nLat
		  do I=1,nLng
		    XCshift(I,J,K,N,2) = 0.0		    
		    if(K.eq.1) then
		      BMAPR(I,J,N) = 0.0
		      BMAPT(I,J,N) = 0.0
		      BMAPN(I,J,N) = 0.0
		    end if
		  end do
		end do
	      end do
	    end do

	    VT       = 100.0*Sun__Omega*Sun__AU
	    BMAPN    = 0.0
	    XCtimmf  = 0.0

C	    ND(1)    = nLng
C	    ND(2)    = nLat
C	    ND(3)    = nMap
C	    ND(4)    = NTG

	    XCbegfst = XCbegG(1,1)
	    XCendfst = XCbegG(2,1) 
	    nLat1    = nLat-1
	    nLng1    = nLng-1
	    nMap1    = nMap-1
	    nT1      = NTG -1
	    nfilesEXT = ((ISFilef - ISFile + IINC)*IINCD2)/IINC
	    write(*,'(A,I5)' ) 'External files numbering system, some may be skipped. nfilesEXT =', nfilesEXT
	    print *, ' '

	    do N=1,NTG
	      XCtimm = XCintG(N) 						! Times the matrix is to be extracted this N
	      if(XCtimm.ge.(TT3D(LLbeg+ifn)+0.0001).and.XCTimmf.eq.0.0) then
		XCtimmf = sngl(XCtimm)									! fill in one ifn beyond beginning
		write(*,'(A,F16.12)') '             Beginning of External fill, XCtimmf =',XCtimmf
	      end if

	      write(*,'(A,I5,F16.12)') 'Value of Carrington time              N, XCtimm', N, XCtimm
	      if(XCtimm.lt.(TT3D(LLbeg+ifn)-0.0001).or.XCtimm.gt.(TT3D(LLbeg+ifn+nfilesEXT)+0.0001)) then ! fill in 1 ifn beyond beginning

	      if(XCtimm.gt.(TT3D(LLbeg+ifn+nfilesEXT)+0.0001).and.XCtimm.lt.(TT3D(LLbeg+ifn+nfilesEXT+1)-0.0001)) print *, ' '
	      write(*,'(A,I5,F16.12)') 'Fill this Carringtion time internally N, XCtimm', N, XCtimm

		do K=1, nMap
		  if(K.lt.KPF) then			! This sets the first surface to the value of the magnetic surface
		    if(K.eq.1) RHIGHT = RRMS - dRRMS
		    RHIGHT = RHIGHT + dRRMS		! K < KPF are the Magnetic Map heights
		    FBFAC = RHIGHT
		  else
		    if(K.eq.KPF) RHIGHT = RRS - dRR
	            RHIGHT = RHIGHT + dRR		! K >= KPF are the deconvolution heights 
		    FBFAC  = RHIGHT
		  end if
	    	  FN  = FBFAC**FALLOFFN							! For the internal source
		  FT  = FBFAC**FALLOFFT
	    	  FNDMAP  = (RRS/RRMS)**FALLOFFN
	          FTDMAP  = (RRS/RRMS)**FALLOFFT
		  FNLMAP  = (RRS/RHIGHT)**FALLOFFN
		  FTLMAP  = (RRS/RHIGHT)**FALLOFFT
		  FBR = FBFAC**FALLOFFBR
		  FBT = FBFAC**FALLOFFBT
		  FBN = FBFAC**FALLOFFBN

		  do J=1,nLat
		    XLatt  = 90.0*(-1.0+2.0*(J-1)/(nLat1))				! Latitude the matrix is to be extracted
		    do I=1,nLng
		      nLng1 = nLng -1
		      XCfraci = 1.0*(nLng-I)/nLng1					! Index [1,nLng] -> fraction of rotation
		      XClng = XCbegfst + XCfraci*(XCendfst-XCbegfst)			! Longitude the matrix is to be extracted

C		      PP(1) = 1+(XCendfst-XCLng)/(XCendfst-XCbegfst)*nLng1		! Carrington variable -> index [1,nLng]
C		      PP(2) = 1+((XLatt+90.0)/180.0)*nLat1				! Latitude (deg) -> index in [1,nLat]
C		      PP(3) = 1+(RHIGHT-RRS)/dRR
C		      PP(4) = 1+sngl(XCtimm-XCtbegG)/sngl(XCtendG-XCtbegG)*nT1		! CV -> time interval
		      if(K.ge.KPF) then
			do NX=1,3
C			  XCsh(NX) = FLINT(-4,ND,XCshift(1,1,1,1,NX),PP,0.00002)		! Minus sign forces check for bad values
			  XCsh(NX) = XCshift(I,J,K-KPF+1,N,NX)   
			end do
C			Drat       = FLINT(-4,ND,DDfact (1,1,1,1   ),PP,0.00002)		! Minus sign forces check for bad values
			Drat	 = DDfact(I,J,K-KPF+1,N) 
C			Vrat       = FLINT(-4,ND,DVfact (1,1,1,1   ),PP,0.00002)		! Minus sign forces check for bad values
			Vrat	 = DVfact(I,J,K-KPF+1,N) 
		      end if
		      XC = XClng + XCsh(1)				! To get the source surface projected center XCint spatial values
		      XC = max(XCbegfst,min(XC,XCendfst))

		      XLa = XLatt + XCsh(2)				! To get the source surface projected heliographic latitude values

		      XCtt = XCintDG(N) + dble(XCsh(3))			! To get the source surface projected center XCint time values
		      XCtt = max(XCtbegG,min(XCtt,XCtendG))

		      if(K.eq.3.and.I.eq.20.and.J.eq.6)  write(*,'(A,I4)') 'Test values at a given distance. K =',K
		      if(K.eq.3.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I4,3F10.5)'),
     &			 'N, K, J, I,  XCsh1,2,3  ',N,K,J,I,XCsh(1),XCsh(2),XCsh(3)

C		      if(K.eq.3.and.I.eq.20.and.J.eq.6)  write(*,'(A,F11.5,1F12.5)'),
C     &			 'Other parameter 22',Vrat,Drat

		      if(K.eq.3.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I4,3F10.5)'),
     &			 'N, K, J, I,XC, Xla, XCtt',N,K,J,I,XC,Xla,XCtt

		      call Get3DTval(XCbeGG,XCtbegG,XCtendG,nLng,nLat,nTG,DMAP, 1,XC,XLa,XCtt,DEN)
		      call Get3DTval(XCbeGG,XCtbegV,XCtendV,nLng,nLat,nTV,VMAPD,1,XC,XLa,XCtt,VEL)
C		      call Get3DTval(XCbeGG,XCtbegV,XCtendV,nLng,nLat,nT,BMAPR,1,XC,XLa,XCtt,BMAPRR)	! Not needed because from write
C		      call Get3DTval(XCbeGG,XCtbegV,XCtendV,nLng,nLat,nT,BMAPT,1,XC,XLa,XCtt,BMAPTT)
C		      call Get3DTval(XCbeGG,XCtbegV,XCtendV,nLng,nLat,nT,BMAPN,1,XC,XLa,XCtt,BMAPNN)
C
		      BMAPRR = 0.0					! BMAPRR will be needed from Write3D_bbtm4D some day.
		      BMAPTT = 0.0
		      BMAPNN = 0.0

		      if(K.lt.KPF) then					! Deal with the maps below the source surface approximately
			VV3(I,J,K,N,1) = VMAP(I,J,N)
			VV3(I,J,K,N,2) = 0.0
			VV3(I,J,K,N,3) = 0.0
			DD1(I,J,K,N)   = DMAP(I,J,N)*FNDMAP		! Set to higher density, but never used 
			TT1(I,J,K,N)   = 10000.0*FTDMAP
			if(RRS.eq.RRMS) then				! This will not happen
			  BB3(I,J,K,N,1) = BMAPR(I,J,N)
			  BB3(I,J,K,N,2) = BMAPT(I,J,N)
			  BB3(I,J,K,N,3) = BMAPN(I,J,N)
			else
C			  BB3(I,J,K,N,1) = BR2DT(I,J,N*ifn)
			  BB3(I,J,K,N,1) = BADR4()		! This will happen. No maps at magnetic surface at K=1, no T or N yet
			  BB3(I,J,K,N,2) = BADR4()			! No BR2DT(I,J,N*ifn) available beyond region near center
			  BB3(I,J,K,N,3) = BADR4()	
C			  BB3(I,J,K,N,1) = BMAPR(I,J,N)*FBR
C			  BB3(I,J,K,N,2) = BMAPT(I,J,N)*FBT
C			  BB3(I,J,K,N,3) = BMAPN(I,J,N)*FBN			  
			end if
		      else						! Maps at all other solar distances
			if(VEL.eq.BadR4().or.Vrat.eq.BadR4()) then
			  VV3(I,J,K,N,1) = BadR4()
			else
			  VV3(I,J,K,N,1) = VEL*Vrat
			end if
			VV3(I,J,K,N,2) = 0.0
			VV3(I,J,K,N,3) = 0.0
			if(Den.eq.BadR4().or.Drat.eq.BadR4()) then
			  DD1(I,J,K,N)   = BadR4()
			else
			  DD1(I,J,K,N)   = (DEN*Drat)*FNLMAP
			end if
			if(I.eq.20.and.J.eq.6.and.N.eq.10.and.K.le.4) write(*,'(A,I4)') 'Test values over distance at time N =',N
			if(I.eq.20.and.J.eq.6.and.N.eq.10) write(*,'(A,4I4,F10.7,2F10.2,F10.5))'), 
     &			 'N, K, J, I, Drat DEN D11',N,K,J,I,Drat,DEN,DD1(I,J,K,N),FNLMAP

			if(K.eq.3.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I4,F10.5,F10.4))'), 
     &			 'N, K, J, I, FN  DD1     ',N,K,J,I,FN,DD1(I,J,K,N)
			if(K.eq.3.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I4,F10.7,2F10.4))'), 
     &			 'N, K, J, I, Drat DEN DD1',N,K,J,I,Drat,DEN,DD1(I,J,K,N)
			if(K.eq.3.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I4,F10.7,2F10.2))'), 
     &			 'N, K, J, I, Vrat VEL VV3',N,K,J,I,Vrat,VEL,VV3(I,J,K,N,1)
			if(K.eq.3.and.I.eq.20.and.J.eq.6)  write(*,'(A,4I4,F10.5,F10.4))'), 
     &			 'N, K, J, I,  FN  DD1    ',N,K,J,I,FN,DD1(I,J,K,N)
			if(K.eq.3.and.I.eq.20.and.J.eq.6) print *, ' '

			TT1(I,J,K,N)   = 10000.0*FTLMAP
			  BB3(I,J,K,N,1) = BADR4()
			  BB3(I,J,K,N,2) = BADR4()
			  BB3(I,J,K,N,3) = BADR4()
C			  BB3(I,J,K,N,1) = BR3D(I,J,K,N*ifn,1)*FBR			! From write above. This not available before.
C			  BB3(I,J,K,N,2) = 0.0*FBT					! N needs to be fixed. No values of B beyond 
C			  BB3(I,J,K,N,3) = BT3D(I,J,K,N*ifn,3)*FBN			! immediate time region, sorry. Use gridsphere.
		      end if

			  if(K.eq.1.and.I.eq.20.and.J.eq.6.and.N.eq.10) print *, ' '
			  if(K.eq.1.and.I.eq.20.and.J.eq.6.and.N.eq.10) write(*,'(A)') 'Check values at distance K = 1'
			  if(K.eq.1.and.I.eq.20.and.J.eq.6.and.N.eq.10) write(*,'(A,4I4,F10.7,2F10.2,F10.5))'), 
     &			   'N,K,J,I,  Drat DEN DD1',N,K,J,I,Drat,DEN,DD1(I,J,K,N),RHIGHT
			  if(K.eq.1.and.I.eq.20.and.J.eq.6.and.N.eq.10) write(*,'(A,4I4,2E10.3,2E10.3))'), 
     &			   'N,K,J,I,BB31 BB32 BB33',N,K,J,I,BB3(I,J,K,N,1),BB3(I,J,K,N,2),BB3(I,J,K,N,3)
			  if(K.eq.nMap.and.I.eq.20.and.J.eq.6.and.N.eq.10) print *, ' '
			  if(K.eq.1.and.I.eq.20.and.J.eq.6.and.N.eq.10) print *, ' '

		    end do
		  end do
		end do
	      end if
	    end do 

	    print *, 'Provide internal source parameters for the new traceback matrix'
	    print *, ' '

	    do N=1,NTG
	      XCtimm = XCintG(N) 						! Times the matrix is to be extracted this N
C	      if(XCtimm.lt.(TT3D(LLbeg+ifn)+0.0001) .or. XCtimm.gt.(TT3D(LLend)-0.0001)) then		    ! fill in one beyond beginning
	      if(XCtimm.lt.(TT3D(LLbeg+ifn)+0.0001) .or. XCtimm.gt.(TT3D(LLbeg+ifn+nfilesEXT)-0.0001)) then ! fill in 1 beyond beginning
		do J=1,nLat
		  do I=1,nLng
		    DMAPV(I,J,N) = DMAP(I,J,N)
		    DMAP (I,J,N) = DMAP(I,J,N)
		    VMAP (I,J,N) = VMAP(I,J,N)
		    if(I.eq.20.and.J.eq.6) 
     &			write(*,'(A,3I4,2F12.1)'), 'N , J, I, DMAPV, DMAP',N,J,I,DMAPV(I,J,N),DMAP(I,J,N)
		    do KK =1,nMap
		      if(KK.ge.KPF) then
		        XCshift3(I,J,KK,N,1) = XCshift(I,J,KK-KPF+1,N,1)	! fill in all shifts at and above KPF with old tomography
		        XCshift3(I,J,KK,N,2) = XCshift(I,J,KK-KPF+1,N,2)
		        XCshift3(I,J,KK,N,3) = XCshift(I,J,KK-KPF+1,N,3)
		      end if
		      if(KK.lt.KPF) then
		        XCshift3(I,J,KK,N,1) = XCshift(I,J,KK+1,N,1)		! fill in the shifts below KPF with the tomography
		        XCshift3(I,J,KK,N,2) = XCshift(I,J,KK+1,N,2)		! (the shifts from above are on the same level as the
		        XCshift3(I,J,KK,N,3) = XCshift(I,J,KK+1,N,3)		!  			parameter values below, below KPF)
			VV3(I,J,KK,N,1)    = VV3(I,J,KK+KPF-1,N,1)		! fill in the velocities below KPF with the old tomography
			VV3(I,J,KK,N,2)    = VV3(I,J,KK+KPF-1,N,2)	
			VV3(I,J,KK,N,3)    = VV3(I,J,KK+KPF-1,N,3)	
			DD1(I,J,KK,N)      = DD1(I,J,KK+KPF-1,N)*FNDMAP		! fill in the densities below KPF with the old tomography
			TT1(I,J,KK,N)      = TT1(I,J,KK+KPF-1,N)*FTDMAP		! fill in the temperatures below KPF in an ad hoc way
			if(RRMS.eq.RRS) then
			  BB3(I,J,K,N,1)   = BB3(I,J,KK+KPF-1,N,1)
			  BB3(I,J,K,N,2)   = BB3(I,J,KK+KPF-1,N,2)
			  BB3(I,J,K,N,3)   = BB3(I,J,KK+KPF-1,N,3)
			else
			  BB3(I,J,K,N,1)   = BADR4()
			  BB3(I,J,K,N,2)   = BADR4()
			  BB3(I,J,K,N,3)   = BADR4()
			end if
		      end if
		    end do
		  end do
		end do
	      end if
	    end do 

	    print *, ' '
	    do N=1,NTG
	      XCtimm = XCintG(N) 						! Times the matrix is to be extracted this N
	      iDoy = int(XCintDG(N))
	      call DATE_DOY(1,IYRBG,cMon,mo,id,iDoy)
	      ih = nint((XCintDG(N) - float(iDoy))*24.0)
	      if(XCtimm.gt.(TT3D(LLbeg+ifn)-0.0001) .and. XCtimm.lt.(TT3D(LLbeg+ifn+nfilesEXT)+0.0001))
     &		write (*,'(A,5I6)'), 'Arrays filled in externally, N, Times', N, IYRBG, mo, id, ih
	      if(XCtimm.lt.(TT3D(LLbeg+ifn)+0.0001) .or. XCtimm.gt.(TT3D(LLbeg+ifn+nfilesEXT)-0.0001)) then  ! fill in 1 beyond beginning
		write(*,'(A,I4,F18.14)') 'Arrays filled in from tomography - N, XCtimm', N,  XCtimm
		do J=1,nLat
		  do I=1,nLng
		    do KK =1,nMap
		      if(KK.ge.1.and.I.eq.20.and.J.eq.6.and.N.eq.10)  then
			if(KK.eq.1.and.I.eq.20.and.J.eq.6.and.N.eq.10) print *, ' '
			if(KK.eq.1.and.I.eq.20.and.J.eq.6.and.N.eq.10) print *, 'Check of parameters over distance at time N = 10'
			write(*,'(A,4I4,3F10.2))'), 
     &			  'KK, N, J, I, VV3, DD1, TT1',KK,N,J,I,VV3(I,J,KK,N,1),DD1(I,J,KK,N),TT1(I,J,KK,N)
			if(KK.eq.nMap+1.and.I.eq.20.and.J.eq.6.and.N.eq.10) print *, ' '
		      end if
		    end do
		  end do
		end do
	      end if
	    end do

C	Fill as if from an external source (could be)

	    print *, ' '
	    write(*,'(A,2I5)') 'Fill in between LLbeg+ifn and LLend from an "external" source', LLbeg+ifn, LLend
	    print *, ' '

	    LLfstp = 10000
	    LL = -10000
	    do N=1,LLend,ifn					! go through time arrays from 1 to LLend at steps of ifn
	      LLfbeg = N/ifn
	      XCtimm = XCintG(LLfbeg)				! Carrington variable times the matrix is to be extracted this N

	      if(XCtimm.gt.(TT3D(LLbeg+ifn)-0.0001) .and. XCtimm.lt.(TT3D(LLbeg+ifn+nfilesEXT)+0.0001)) then
		if(LL.lt.0) then
		  LLfstp = N - LLfbeg*ifn 
		  LLbegg = N - ifn - 2					! Strange this it what is needed
		  LBG = LLbegg - ifn
		  LLLfstp = ((N + LLfstp)/ifn)
		  LL = LLLfstp - 1
	        end if
	        LBG = LBG + ifn
	        LL  = LL + 1

C		print *, 'Before the external read 0, LL =', 
C     &			LL, N, LLbeg, LLend, LLbegg, LLfstp, LLfbeg, XCtimm, XCtimmf, TT3D(N), ISfile 
	
C	************************************

C		print *, 'Before the external read 1, LL =', LL 

		if(iYes.eq.0) then
		  call ExternalRead(cWild3DMHD,MHDs,nLng,nLat,nMap,nTmaxG,LLBegg,ISfile,ISfilef,LLLfstp,LLEnd,NinterD,IINC,IINCD2,
     &			      IYRBG,XCintDG,TT3D,
     &			      RRS,RADMS,RAD1,dRR,RADS,FALLOFFN,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,DDD1,TTT1,VVV3,BBB3,iYes) 
		end if
		
		if(iYes.eq.1) then

C		  print *, 'Before external read 2, LL =', LL
		  call ExternalRead(cWild3DMHD,MHDs,nLng,nLat,nMap,nTmaxG,LLBegg,ISfile,ISfilef,LLLfstp,LLEnd,NinterD,IINC,IINCD2,
     &			      IYRBG,XCintDG,TT3D,
     &			      RRS,RADMS,RAD1,dRR,RADS,FALLOFFN,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,DDD1,TTT1,VVV3,BBB3,iYes) 
		    print *, 'Was this a successful read? iYes =', iYes
		  if(iYes.eq.2) then
		    print *, 'This was a successful read. iYes =', iYes
		  end if
		end if

C	************************************

		write(*,'(A,4I5)') 'Filling external arrays N, LL, LBG ',N,LL,LBG

	        do K=1, nMAP+1				! go through heights
		  if(K.lt.KPF) then			! This sets the first surface to the value of the magnetic surface
		    KK = K
		    if(RRS.eq.RRMS) then
		      FBFAC = (RRS + dRRMS*(float(K-1)))/RRMS
		      RHIGHT = RRS + dRRMS*(float(K-KPF+1))
		    else
		      if(K.eq.1) RHIGHT = RRMS - dRRMS		! This begins K .lt. KPF below the deconvolution surface
		      RHIGHT = RRMS + dRRMS
		      FBFAC  = RHIGHT
		    end if
		  else
		    if(K.eq.KPF) RHIGHT = RRS - dRR		! Will this work now? BVJ
		    KK = K - KPF + 1
	            RHIGHT = RHIGHT + dRR
		    FBFAC  = RHIGHT
		  end if 
		  FN  = FBFAC**FALLOFFN
		  FT  = FBFAC**FALLOFFT
		  FBR = FBFAC**FALLOFFBR
		  FBT = FBFAC**FALLOFFBT
		  FBN = FBFAC**FALLOFFBN

C		  print *, 'N, K, KK',N, K, KK 

		  do J=1,nLat
		    do I=1,nLng
		      if(K.eq.1) then					! Imported values of DMAPV, etc. possibly above magnetic surface
			if(iYes.eq.3.or.iYes.eq.-3) then					! If iYes=3, -3 these are tests
			  DMAPV(I,J,LL) = D3DT(I,J,1,LBG)
			  DMAP (I,J,LL) = D3DT(I,J,1,LBG)
		     	  VMAP (I,J,LL) = V3DT(I,J,1,LBG)
		        end if
		      end if
		        if(iYes.eq.3.or.iYes.eq.-3) then				! If iYes=3 these are test writes
		          if(K.eq.2.and.I.eq.20.and.J.eq.6) write(*,'(A,4I4,3F10.2)'), 
     &			    'LL, LBG, J, I, DMAPV,DMAP,VMAP',LL,LBG,J,I,DMAPV(I,J,LL),DMAP(I,J,LL),VMAP(I,J,LL)
		          if(I.eq.20.and.J.eq.6) write(*,'(A,4I4,3F10.2)'), 
     &			    'I,J,K,LBG,BB3(1),BB3(2),BB3(3)',I, J, K, LBG, BB3(I,J,K,LBG,1),BB3(I,J,K,LBG,2),BB3(I,J,K,LBG,3) 

			  VV3(I,J,K,LL,1) = V3DT(I,J,KK,LL)				! Test Imported values of V over height (KK)
			  VV3(I,J,K,LL,2) = 0.0
			  VV3(I,J,K,LL,3) = 0.0
			  if(D3DT(I,J,KK,LL).ne.BADR4()) then
			    DD1(I,J,K,LL)   = D3DT(I,J,KK,LBG)*FN			! Test Imported values of D over height (KK)
			  else
			    DD1(I,J,K,LL)   = BADR4()
			  end if
			  TT1(I,J,K,LL)   = 10000.0/FT					! Test Imported values of T over height (KK)
			  BB3(I,J,K,LL,1) = BB3(I,J,K,LBG,1)				! Imported values of B over height (KK)
		 	  BB3(I,J,K,LL,2) = BB3(I,J,K,LBG,2)
C			  BB3(I,J,K,LL,3) = BB3(I,J,K,LBG,3)
			  BB3(I,J,K,LL,3) = 0.0
			else								! Reads from external input
			  if(K.eq.1) then
			    if(VVV3(I,J,2,LL,1).lt.9999.9) then
		              VV3(I,J,K,LL,1) = VVV3(I,J,2,LL,1)			! Imported values of V over height (KK)
			    else
			      VV3(I,J,K,LL,1) = BADR4()
			    end if
		            VV3(I,J,K,LL,2) = VVV3(I,J,2,LL,2)
		            VV3(I,J,K,LL,3) = VVV3(I,J,2,LL,3)
			    DD1(I,J,K,LL)   = DDD1(I,J,2,LL)				! Imported values of D over height (KK)
			    TT1(I,J,K,LL)   = TTT1(I,J,2,LL)
			    BB3(I,J,K,LL,1) = BBB3(I,J,1,LL,1)				! Imported values of B over height (KK)
		 	    BB3(I,J,K,LL,2) = BBB3(I,J,1,LL,2)
C			    BB3(I,J,K,LL,3) = BBB3(I,J,1,LL,3)
			    BB3(I,J,K,LL,3) = BADR4()
			      VMAP (I,J,LL) = VVV3(I,J,2,LL,1)			      
			      DMAP (I,J,LL) = DDD1(I,J,2,LL)
		     	      DMAPV(I,J,LL) = DDD1(I,J,2,LL)
			  else
		            VV3(I,J,K,LL,1) = VVV3(I,J,K,LL,1)				! Imported values of V over height (KK)
		            VV3(I,J,K,LL,2) = VVV3(I,J,K,LL,2)
		            VV3(I,J,K,LL,3) = VVV3(I,J,K,LL,3)
			    DD1(I,J,K,LL)   = DDD1(I,J,K,LL)				! Imported values of D over height (KK)
			    TT1(I,J,K,LL)   = TTT1(I,J,K,LL)				! Imported values of T over height (KK)
			    BB3(I,J,K,LL,1) = BBB3(I,J,K,LL,1)				! Imported values of B over height (KK)
		 	    BB3(I,J,K,LL,2) = BBB3(I,J,K,LL,2)
C			    BB3(I,J,K,LL,3) = BBB3(I,J,K,LL,3)
			    BB3(I,J,K,LL,3) = BADR4()
			  end if
			end if

		        if(K.eq.2.and.I.eq.20.and.J.eq.6) write(*,'(A,4I4,2F10.1,F10.5)'), 
     &			  'LL, LBG, J, I, DMAPV,DMAP,RHIGHT',LL,LBG,J,I,DMAPV(I,J,LL),DMAP(I,J,LL),RHIGHT

		    end do
		  end do  

		  if(iYes.eq.3.or.iYes.eq.-3)	then						   ! If iYes=3, -3 these are test writes
	            call GridSphere2D(ALng,nLng,nLat,1,BB3(1,1,KK,LL,1),CONRD,4,0.0,0.0)	   ! Smooth Br polar holes, elsewhere
	            call GridSphere2D(ALng,nLng,nLat,1,BB3(1,1,KK,LL,2),CONRD,4,0.0,0.0)	   ! Smooth Bt polar holes, elsewhere
	            call GridSphere2D(ALng,nLng,nLat,1,BB3(1,1,KK,LL,3),CONRD,4,0.0,0.0)	   ! Smooth Bn polar holes, elsewhere
		    do J=1,nLat
		      do I=1,nLng
			BB3(I,J,K,LL,1) = BB3(I,J,KK,LL,1)				! Test Imported values of B over height (KK)
		 	BB3(I,J,K,LL,2) = BB3(I,J,KK,LL,2)				! Test First height assumed the same as next step
			BB3(I,J,K,LL,3) = BB3(I,J,KK,LL,3)

		        if(K.eq.2.and.I.eq.20.and.J.eq.6) write(*,'(A,3I4,A,F12.1)'), 
     &			 'LL, J, I,        DD1  ',LL,J,I,'            ',DD1(I,J,K,LL)
		        if(K.eq.2.and.I.eq.20.and.J.eq.6) write(*,'(A,3I4,A,F12.1)'), 
     &			 'LL, J, I,        DD1  ',LL,J,I,'            ',TT1(I,J,K,LL)
		        if(K.eq.2.and.I.eq.20.and.J.eq.6) write(*,'(A,3I4,A,F12.1)'), 
     &			 'LL, J, I,        VV3  ',LL,J,I,'            ',VV3(I,J,K,LL,1)
		        if(K.eq.2.and.I.eq.20.and.J.eq.6) write(*,'(A,3I4,A,F12.1)'), 
     &			 'LL, J, I,        VV3  ',LL,J,I,'            ',BB3(I,J,K,LL,1)

		        if(K.ge.1.and.I.eq.20.and.J.eq.6.and.LL.eq.22)  then
			  if(K.eq.1.and.I.eq.20.and.J.eq.6.and.LL.eq.22) print *, ' '
			  write(*,'(A,5I4,4F7.2))'), 
     &			   'K,KK,N,J,I, DD1, TT1, VV3, BB3',K, KK, N,J,I,DD1(I,J,K,LL),TT1(I,J,K,LL)/10000.0,VV3(I,J,K,LL,1),
     &			  BB3(I,J,K,LL,1)
			  if(K.eq.nMap+1.and.I.eq.20.and.J.eq.6.and.LL.eq.22) print *, ' '
		        end if
		      end do
		    end do
		  end if		
	        end do
	      end if
	    end do

	    print *, 'Before External Write iYes =', iYes, LLLfstp  

	    if(iyes.eq.3) then
	      print *, NTG,NTV,NTP,LLbegg,LLLfstp,LLend
	      call ExternalWrite(MHDs,nLng,nLat,nMap,nTmaxG,LLbegg,ISfile,ISfilef,LLLfstp,LLend,NinterD,IINC,IINCD2,
     &				IYRBG,XCintDG,TT3D,
     &			     	RRS,RADMS,RAD1,dRR,FALLOFFN,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,DD1,TT1,VV3,BB3)	       
	    end if

	    print *, ' '

	    if(bmagex) then

C	      print *, 'DATA EXTRACTED FROM FINAL MATRIX INCLUDING Magnetic fields'

	      print *, ' Before Final Extractdvdm',XCbegG(1,1),XCtbegG,XCstrt

	      print *, ' '
	      print *, 'Check the array disposition at the end of the tomography - XCshift 1, 2, 3'
	      print *, ' '
		  do N=1,nTG
		    do K=1, nMap
		      do J=1,nLat
			do I=1,nLng

		      if(I.eq.20.and.J.eq.6.and.N.eq.22) write(*,'(A,4I4,3F10.4)'), 
     &			 'N, K, J, I,      XCshift 1-3  ',N, K, J, I,XCshift(I,J,K,N,1),XCshift(I,J,K,N,2),XCshift(I,J,K,N,3)
C			  XCshift(I,J,K,N,1) = XCshift3(I,J,K+1,N,1)
C			  XCshift(I,J,K,N,2) = 0.0
C			  XCshift(I,J,K,N,3) = XCshift3(I,J,K+1,N,3)
C			  DVfact (I,J,K,N) = Vratio3 (I,J,K+1,N,1)
C			  DDfact (I,J,K,N) = Dratio1 (I,J,K+1,N)
			end do
		      end do
		    end do
		  end do
	      print *, ' '
	      print *, 'Check the array disposition at the end of the tomography - DVfact, DDfact'
	      print *, ' '
		  do N=1,nTG
		    do K=1, nMap
		      do J=1,nLat
			do I=1,nLng
		      if(I.eq.20.and.J.eq.6.and.N.eq.22) write(*,'(A,4I4,2F10.4)'), 
     &			 'N, K, J, I, DVfact,  DDFact   ',N, K, J, I,DVfact (I,J,K,N),DDfact (I,J,K,N)
			end do
		      end do
		    end do
		  end do

	      print *, ' '
	      print *, 'Check the array disposition at the end of the tomography - VV3, DD1'
	      print *, ' '
		  do N=1,nTG
		    do K=1, nMap
		      do J=1,nLat
			do I=1,nLng
		      if(I.eq.20.and.J.eq.6.and.N.eq.22) write(*,'(A,4I4,2F10.4)'), 
     &			 'N, K, J, I,VV3(I,J,K,N,1), DD1',N, K, J, I,VV3 (I,J,K,N,1),DD1(I,J,K,N)
			end do
		      end do
		    end do
		  end do
	      print *, ' '
	      print *, 'Check the array disposition at the end of the tomography - DMAPV, DMAP'
	      print *, ' '
		  do N=1,nTG
			I=20
			J=6
		      write(*,'(A,3I4,2F10.2,F10.5)'), 
     &			 'N, J, I, DMAPV(I,J,N) DMAP    ',N, J, I,DMAPV(I,J,N),DMAP (I,J,N),RADS(2)
		  end do
	      print *, ' '

	  write(*,'(A,I3,A,I3,A,I3)') 'The values of Nintert =',NiterT,' NiterT2 =',NiterT2,' and Nit2 =',Nit2
	STOP
	if(Nit2.eq.2) then
	  NiterT2 = NiterT2 + 1
	  if(Nitert2.eq.19) go to 9997
	  Nit2 = 1	! Provide an old kinematic model traceback with the new base maps 		  
	  NiterT = 1
	  Nit    = 1
	  NitV   = 1
	  NitG   = 1
C        	  call MkShiftdn0n(XCbeGG,XCintG,XCtbegG,XCtendG,ALng,nLng,nLat,nMap,nTG,nTmaxG,Speed,
C     &		  nCar,JDCar,NCoff,VMAPD,DMAP,XCshift,DVfact,DDfact,RRS,dRR,FALLOFF,CONRV,CONRD,VLL,VUL,
C    &		  NSIDE,Vtmp,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt)

	  go to 9988		! Provide a revised set of base maps and a 4D density matrix for a new traceback iteration
	end if

	if(Nit2.eq.0.or.Nit2.eq.1) then
	  Nit2 = 2	! After resetting NiterT, go around again with the new analysis input and begin iteration with new traceback.
	  NiterT = 1
	  Nit    = 0
	  NitV   = 0
	  NitG   = 0

	if(MHDs.eq.1) then
	  write(*,'(A,I3,A)') 'New number of iterations NiterT =',NiterT,' with array input' 
	  write(*,'(A)') 'The program is ready to continue with the input external data from ENLIL'
	  write(*,'(A)') 'almost'
	  stop
	end if

		dRRMS	  = RAD1 - RADMS
		RRMS      = RADMS
		KPF       = 2
		FALLOFFN  = 2.0
		FALLOFFT  = 2.0
		FALLOFFBR = 2.0
		FALLOFFBT = 1.0
		FALLOFFBN = 1.34
		VLLLL	  = 50.0
		VLLUL	  = -50.0

		  call MkShiftdnma(XCbeGG,XCtbegG,XCtendG,ALng,nLng,nLat,nMap,nTG,
     &		  VV3,BB3,DD1,TT1,XCshift3,Vratio3,Bratio3,Dratio1,Tratio1,RRS,RRMS,dRR,dRRMS,KPF,
     &		  FALLOFFN,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,CONRV,CONRD,VLL,VUL,VLLLL,VLLUL,
     &		  Vtmp,Dtmp,XLLTtmp)

C		  iprint = 0
		  iprint = 1
		  if(iprint.eq.1) then
	          KPFM1 = KPF-1
	          print *, ' '
	          write(*,'(A,I4)'), ' Check the array disposition after MkShiftdnma KPFM1 =', KPF-1
	          print *, ' '
		  do N=1,nTG
		    do K=1, nMap
		      do J=1,nLat
			do I=1,nLng

		        if(I.eq.20.and.J.eq.6.and.N.eq.22) write(*,'(A,4I4,3F10.4)'), 
     &			 'N, K, J, I,      XCshift 1-3  ',N, K, J, I,XCshift(I,J,K,N,1),XCshift(I,J,K,N,2),XCshift(I,J,K,N,3)
		         if(I.eq.20.and.J.eq.6.and.N.eq.22) write(*,'(A,4I4,3F10.4)'), 
     &			 'N, K, J, I,      XCshift3 1-3  ',N, K, J, I,XCshift3(I,J,K+KPFM1,N,1),
     &				XCshift3(I,J,K+KPFM1,N,2),XCshift3(I,J,K+KPFM1,N,3)
C			  XCshift(I,J,K,N,1) = XCshift3(I,J,K+KPFM1,N,1)
C			  XCshift(I,J,K,N,2) = 0.0
C			  XCshift(I,J,K,N,3) = XCshift3(I,J,K+KPFM1,N,3)
C			  DVfact (I,J,K,N) = Vratio3 (I,J,K+KPFM1,N,1)
C			  DDfact (I,J,K,N) = Dratio1 (I,J,K+KPFM1,N)
			end do
		      end do
		    end do
		  end do
		K = 1
		      if(I.eq.20.and.J.eq.6.and.N.eq.22.and.K.eq.1) write(*,'(A,4I4,3F10.4)'), 
     &			 'N, K, J, I,      XCshift3 1-3  ',N, K, J, I,XCshift3(I,J,K,N,1),XCshift3(I,J,K,N,2),XCshift3(I,J,K,N,3)

		print *, ' '
		  do N=1,nTG
		    do K=1, nMap
		      do J=1,nLat
			do I=1,nLng
		      if(I.eq.20.and.J.eq.6.and.N.eq.22) write(*,'(A,4I4,2F10.4)'), 
     &			 'N, K, J, I, DVfact,  DDFact   ',N, K, J, I,DVfact (I,J,K,N),DDfact (I,J,K,N)
		      if(I.eq.20.and.J.eq.6.and.N.eq.22) write(*,'(A,4I4,2F10.4)'), 
     &			 'N, K, J, I, Vratio3, Dratio1   ',N, K, J, I,Vratio3(I,J,K+KPFM1,N,1),Dratio1(I,J,K+KPFM1,N)
			end do
		      end do
		    end do
		  end do		  

		  end if				! End of print loop

C		The following shifts over the current matrix and ratios for use in the UCSD tomography

		  print*, 'Transferring to UCSD shift system after MkShiftdnma to return with new shifts'

	          KPFM1 = KPF-1
		  do N=1,nTG
		    do K=1, nMap
		      do J=1,nLat
			do I=1,nLng
			  XCshift(I,J,K,N,1) = XCshift3(I,J,K+KPFM1,N,1)
			  XCshift(I,J,K,N,2) = 0.0
			  XCshift(I,J,K,N,3) = XCshift3(I,J,K+KPFM1,N,3)
			  DVfact (I,J,K,N) = Vratio3 (I,J,K+KPFM1,N,1)
			  DDfact (I,J,K,N) = Dratio1 (I,J,K+KPFM1,N)
			end do
		      end do
		    end do
		  end do
	end if

C	print *, 'out of fillmaptn'
	return
	end