C+
C NAME:
C	ExternalRWMHD.f
C PURPOSE:
C	Set up to write or read the externally-provided volumetric elements.
C	These volumetric elements are usually read from files provided that are interpolated with  
C	the cadence and in the heliographic coordinate system needed by the UCSD tomography program. 
C	These files needed for the UCSD program are expected to be provided at a subdirectory of the 3D-MHD program's run location.
C CATEGORY:
C	I/O
C CALLING SEQUENCE:
C	CALL ExternalRWMHD(MHDs,bForecast,iYes,cWild3DMHD,nLng,nLat,nMap,nT,nTmax,XCbegg,XCbeg,xInc,NCC,
C     &			   LLfst,LLEnd,JDCar,nCar,RRS,dRR,RRMS,FALLOFFD,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,R1AU,
C     &			   BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,FALLOFFBN,iBFlag,BRR2DT,BRC2DT,BTC2DT,BNC2DT,BR3DHR,BT3DHR,BN3DHR,
C     &			   cPre,-Nit,NiterT,NCoff,XCintDG,XCtbegG,XCtendG,iYrBG,
C     &			   PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
C     &			   bDdenerb,bVdenerb,bVvererb,bDvererb,ERDLOSB,ERVLOSB,ANMAPD,ANMAPV,
C     &			   Scale,VMAPD,VMAP,DMAP,VMHR,DMHR,XCshift3,XCshiftM,DVfact,DDfact,V3DHR,D3DHR,
C     &			   DDD1,TTT1,VVV3,BBB3,DD1,TT1,VV3,BB3,iNum)
C
C INPUTS:
C	MHDs				integer		Which 3D MHD program provides the input?      1 - ENLIL, 2 - MSFLUKSS, 3 - HAF3DMHD
C	bForecast			logical		Yes - this is a forecast analysis, No - this is not a forecast analysis
C	iYes				integer		iYes = 0 	Dont write out any files from the modeling.
C							iYes = 1	Write out test external files from the modeling.
C							iYes = 2	Read in files from an external model.
C	cWild3DMHD			character*80	input file name				
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
C	nTmax				integer		Maximum number of times possible
C	XCbegg  (2,nTmax)		real		Beginning and end of each interval in Carrington rotation value
C	XCbeg				real		Beginning of the Carrington rotation time start
C	LLBeg				integer		Beginning number of 6-hour averages
C	xInc				real			
C	NCC				integer		
C	LLFst				integer		Beginning increment of output data
C	LLEnd				integer		Ending increment of output data
C	JDCar				real*8		
C	nCar				integer
C	ISfile				integer		Beginning file number to write out
C	Ninter				integer		Number of intermediate spaces between tomography intervals
C	IINC				integer		3D-MHD intervals per tomographic time
C	IINCD2				integer		3D-MHD intervals per high resolution tomographic time
C	IYRBG				integer		Beginning year of the tomography sequence
C	XCintDG(4*nTmax)		real*8		Day of Year	
C	TT3D(4*nTmax)			real		Times of data (in Carrington Variable) at the high resolution tomography cadence
C	RRS				real		Beginning Height of the IPS modeling
C	RRMS				real		Beginning Height of the Magnetic Field pick-up
C	RAD1				real		Beginning Height of the 3D-MHD modeling
C	dRR				real		Outward step height of the modeling
C	FALLOFFD			real		Fall off in density
C	FALLOFFT			real		Fall off in temperature
C	FALLOFFBR			real		Fall off in Br
C	FALLOFFBT			real		Fall off in Bt
C	FALLOFFBN			real		Fall off in Bn
C
C OUTPUTS:
C	DDD1(nLng,nLat,NMap,nTmax)	real		Input 1 component density
C	TTT1(nLng,nLat,NMap,nTmax)	real		Input 1 component temperature
C	VVV3(nLng,nLat,NMap,nTmax,3)	real		Input 3 component velocity (r,t,n)
C	BBB3(nLng,nLat,NMap,nTmax,3)	real		Input 3 component magnetic field (r,t,n)
C	DD1(nLng,nLat,NMap,nTmax)	real		Input 1 component density
C	TT1(nLng,nLat,NMap,nTmax)	real		Input 1 component temperature
C	VV3(nLng,nLat,NMap,nTmax,3)	real		Input 3 component velocity (r,t,n)
C	BB3(nLng,nLat,NMap,nTmax,3)	real		Input 3 component magnetic field (r,t,n)
C	iNum				integer		Number of input files read
C
C MODIFICATION HISTORY:
C	September -2016, Bernard Jackson (UCSD)
C-
	subroutine ExternalRWMHD(MHDs,bForecast,iYes,cWild3DMHD,nLng,nLat,nMap,nT,nTmax,XCbegg,XCbeg,xInc,NCC,
     &		LLfst,LLEnd,JDCar,nCar,RRS,dRR,RRMS,FALLOFFD,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,R1AU,
     &		BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,iBFlag,BRR2DT,BRC2DT,BTC2DT,BNC2DT,BR3DHR,BT3DHR,BN3DHR,
     &		cPre,Nit,NiterT,NCoff,XCintDG,XCtbegG,XCtendG,iYrBG,
     &		PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
     &		bDdenerb,bVdenerb,bVvererb,bDvererb,ERDLOSB,ERVLOSB,ANMAPD,ANMAPV,
     &		Scale,VMAPD,VMAP,DMAP,VMHR,DMHR,XCshift3,XCshiftM,DVfact,DDfact,V3DHR,D3DHR,
     &		ConsTMD,ConsTMV,ConstL,aNdayG,aNdayV,
     &		DDD1,TTT1,VVV3,BBB3,DD1,TT1,VV3,BB3,iNum)

	real	XCintG(nTmax),					! these are needed in the external write
     &		XCEA(nT),
     &		DDD1(nLng,nLat,NMap,nTmax),
     &		VVV3(nLng,nLat,NMap,nTmax,3),
     &		BBB3(nLng,nLat,NMap,nTmax,3),
     &		TTT1(nLng,nLat,NMap,nTmax),
     &		DD1(nLng,nLat,NMap,nTmax),
     &		VV3(nLng,nLat,NMap,nTmax,3),
     &		BB3(nLng,nLat,NMap,nTmax,3),
     &		TT1(nLng,nLat,NMap,nTmax)

	real    VVV31(nLng,nLat),			!  Scratch arrays
     &		VVV32(nLng,nLat),
     &		VVV33(nLng,nLat),
     &		BBB31(nLng,nLat),
     &		BBB32(nLng,nLat),
     &		BBB33(nLng,nLat),
     &		DDD11(nLng,nLat),
     &		TTT11(nLng,nLat),
     &		VVVV32(nLng,nLat,nMap),
     &		VVVV33(nLng,nLat,nMap),
     &		BBBB31(nLng,nLat,nMap),
     &		BBBB32(nLng,nLat,nMap),
     &		BBBB33(nLng,nLat,nMap),
     &		DDDD1 (nLng,nLat,nMap),
     &		TTTT1 (nLng,nLat,nMap)

	parameter	(NintLn	= 3,				! High resolution output parameters
     &			 NintLa = 3,
     &			 NintH  = 3)

	real	XCshiftM(nLng,nLat,nTmax,3),
     &		BRR2DT  (nLng,nLat,nTmax),
     &		BRC2DT  (nLng,nLat,nTmax),
     &		BTC2DT  (nLng,nLat,nTmax),
     &		BNC2DT  (nLng,nLat,nTmax),
     &		DMHR	(nLng*(NintLn+1)+1,nLat*(NintLa+1)+1),
     &		VMHR	(nLng*(NintLn+1)+1,nLat*(NintLa+1)+1),
     &		V3DHR	(nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1),
     &		D3DHR	(nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1), 
     &		BR3DHR	(nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1),
     &		BT3DHR	(nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1),
     &		BN3DHR	(nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1),
     &		XCsh    (3),
     &		ANMAPD	(nLng,nLat,nTmax),					! Density LOS error source surface map
     &		ANMAPV	(nLng,nLat,nTmax),					! Velocity LOS error source surface map
     &		DMAP    (nLng,nLat,nTmax),
     &		VMAP    (nLng,nLat,nTmax),
     &		VMAPD   (nLng,nLat,nTmax),
     &		XCSHIFT3(nLng,nLat,nMap,nTmax,3),
     &		xcbegg  (2,nTmax),
     &		ddfact  (nLng,nLat,nMap,nTmax),
     &		dvfact  (nLng,nLat,nMap,nTmax),
     &		Scale	(4)

	real*8		JDCar(nCar)
	real*8		XCintDG(nTmax)
	real*8		XCtbegg,XCtendg

	character	cWild3DMHD*80
	character	cPre*4

	logical		bForecast
	logical		bDdenerb
	logical		bVdenerb
	logical		bVvererb
	logical		bDvererb

	external	EARTH8

	print *, ' '
	print *, 'Into Externalrwmhd.f', nT
	Nit = -Nit						! Triggers an internal write. We will always do this.
	print *, 'NiT =', NiT
C
C	Prepare the external file to write
C
	      BadV = BadR4()
	      BadD = BadR4()		
	      BadB = BadR4()
	      BadT = BadR4()

	      ALng = 1.*(nLng-1)/(iLng-1)
	      NintLng = 0
	      NintLat = 0
	      NintHt  = 0
	      Ninter = 0
	      nMapm   = 1
	      nMapHRMm = nMap
C
	      do I=1,nLng
		do J=1,nLat
		  do K=1,nMap
		    do L=1,nT
		      DDD1(I,J,K,L)   = BadD
		      VVV3(I,J,K,L,1) = BadV
		      VVV3(I,J,K,L,2) = 0.0
		      VVV3(I,J,K,L,3) = 0.0
		      BBB3(I,J,K,L,1) = BadB
		      BBB3(I,J,K,L,2) = BadB
		      BBB3(I,J,K,L,3) = BadB
		      TTT1(I,J,K,L) = 10000.0
		    end do
		  end do
		end do
	      end do

	      nLngLat = nLng*nLat
	      RRSCON = (R1AU/RRS)**FALLOFFD

	      if(iYes.eq.2) go to 9999		! Assumes the MHD files to read are complete amd without holes.

	      do N=1,nT
		  NBDDD1 = 0
		  NBVVV1 = 0
		  NBVVV2 = 0
		  NBBBB1 = 0
		  NBBBB2 = 0
	 	  NBBBB3 = 0
		  NBBBB4 = 0

		  if(N.eq.1) then
		    print *, ' '
		    print *, 'In externalrwmhd.f. Before INTERNAL writes. By time, suface maps are checked.'
		  end if

		  call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax)
		  if(amax.ne.BadD) then 
		    do J=1,nLat
		      do I=1,nLng
			if(DMAP(I,J,N).eq.BadD) NBDDD1 = NBDDD1 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,DMAP(I,J,N),CONRD,3,0.0,0.0)		! DMAP Should have no holes unless all bad
		  else
	            write(*,'(A,I3,A,I5,A)') 'N =',N,'  All', nLngLat,
     &		    ' density points were bad this time'
		  end if

		  call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax)
		  if(amax.ne.BadV) then 
		    do J=1,nLat
		      do I=1,nLng
			if(VMAP(I,J,N).eq.BadV) NBVVV1 = NBVVV1 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,VMAP(1,1,N),CONRV,3,0.0,0.0)		! VMAP Should have no holes unless all bad
		  else
	            write(*,'(A,I3,A,I5,A)') 'N =',N,'  All', nLngLat,
     &		    ' velocity points were bad this time'
		  end if

		  call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax)
		  if(amax.ne.BadV) then 
		    do J=1,nLat
		      do I=1,nLng
			if(VMAPD(I,J,N).eq.BadV) NBVVV2 = NBVVV2 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,VMAPD(1,1,N),CONRD,3,0.0,0.0)		! VMAPD Should have no holes unless all bad
		  else
	            write(*,'(A,I3,A,I5,A)') 'N =',N,'  All', nLngLat,
     &		    ' velocity for density points were bad this time'
		  end if

		  call arrR4getminmax(nLngLat,BRR2DT(1,1,N),amin,amax)
		  if(amax.ne.BadB) then 
		    do J=1,nLat
		      do I=1,nLng
			if(BRR2DT(I,J,N).eq.BadB) NBBBB1 = NBBBB1 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,BRR2DT(1,1,N),CONRD,3,0.0,0.0)		! Br Should have no holes unless all bad
		  else
	            write(*,'(A,I3,A,I5,A)') 'N =',N,'  All', nLngLat,
     &		    ' open radial field points were bad this time'
		  end if

		  call arrR4getminmax(nLngLat,BRC2DT(1,1,N),amin,amax)
		  if(amax.ne.BadB) then 
		    do J=1,nLat
		      do I=1,nLng
			if(BRC2DT(I,J,N).eq.BadB) NBBBB2 = NBBBB2 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,BRC2DT(1,1,N),CONRD,3,0.0,0.0)		! BrC Should have no holes unless all bad
		  else
	            write(*,'(A,I3,A,I5,A)') 'N =',N,'  All', nLngLat,
     &		    ' closed radial field points were bad this time'
		  end if

		  call arrR4getminmax(nLngLat,BTC2DT(1,1,N),amin,amax)
		  if(amax.ne.BadB) then 
		    do J=1,nLat
		      do I=1,nLng
			if(BTC2DT(I,J,N).eq.BadB) NBBBB3 = NBBBB3 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,BTC2DT(1,1,N),CONRD,3,0.0,0.0)		! BtC Should have no holes unless all bad
		  else
	            write(*,'(A,I3,A,I5,A)') 'N =',N,'  All', nLngLat,
     &		    ' closed tangential field points were bad this time'
		  end if

		  call arrR4getminmax(nLngLat,BNC2DT(1,1,N),amin,amax)
		  if(amax.ne.BadB) then 
		    do J=1,nLat
		      do I=1,nLng
			if(BNC2DT(I,J,N).eq.BadB) NBBBB4 = NBBBB4 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,BNC2DT(1,1,N),CONRD,3,0.0,0.0)		! Bn Should have no holes unless all bad
		  else
	            write(*,'(A,I3,A,I5,A)') 'N =',N,'  All', nLngLat,
     &		    ' closed normal field points were bad this time'		    
		  end if

	        if(NBDDD1.ne.0) write(*,'(A,I3,A,I5,A,I5,A)') 'N =',N,'  There were',NBDDD1,' bad points of',
     &		  nLngLat,' this DMAP time'
	        if(NBVVV1.ne.0) write(*,'(A,I3,A,I5,A,I5,A)') 'N =',N,'  There were',NBVVV1,' bad points of',
     &		  nLngLat,' this radial VMAP time'
	        if(NBVVV2.ne.0) write(*,'(A,I3,A,I5,A,I5,A)') 'N =',N,'  There were',NBVVV2,' bad points of',
     &		  nLngLat,' this radial VMAPD time'
	        if(NBBBB1.ne.0) write(*,'(A,I3,A,I5,A,I5,A)') 'N =',N,'  There were',NBBBB1,' bad points of',
     &		  nLngLat,' this open radial field time'
	        if(NBBBB2.ne.0) write(*,'(A,I3,A,I5,A,I5,A)') 'N =',N,'  There were',NBBBB2,' bad points of',
     &		  nLngLat,' this closed radial field time'
	        if(NBBBB3.ne.0) write(*,'(A,I3,A,I5,A,I5,A)') 'N =',N,'  There were',NBBBB3,' bad points of',
     &		  nLngLat,' this closed tangential field time'
	        if(NBBBB4.ne.0) write(*,'(A,I3,A,I5,A,I5,A)') 'N =',N,'  There were',NBBBB4,' bad points of',
     &		  nLngLat,' this closed normal field time'

	    end do										! end time do loop
	    print *, ' '

C
C -Nit makes this a 3DMHD test write.
C
	    write(*,'(3(A,F11.4))') 'In ExtR VMAPD(22,5,22)',VMAPD(22,5,22),' VMAP(22,5,22)',VMAP(22,5,22),' DMAP(22,5,22)',DMAP(22,5,22)

	      nLngLat = nLng*nLat
	      RRSCON = (RRS/R1AU)**FALLOFFD
	      do N=1,nT
	        call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N))
	      end do	

	      xIncd = xInc
	      if(.not.bForecast) xIncd = 0.0
C
	      nTD   = nT
	      if(.not.bForecast) then
		ANinterp = Ninter + 1
		XCdif = sngl(XCintDG(2)-XCintDG(1))
		AI    = XCdif/ANinterp
		nTD   = (xInc*AI)/27.2753 + nT
	      end if
	      print *, ' '
	      write (*,'(A,2F5.2,I5)') 'Beginning and end V, D boundary xInc, xIncd, nTD', xInc, xIncd, nTD

	      call write3D_infotd3DM_HR_3DMHD(0,cPre,Nit,NiterT,NintLng,NintLat,NintHt,Ninter,NCoff,XCintDG,XCbegg,XCtbegG,XCtendG,iYrBG,
     &			RRS,dRR,nLng,nLat,nMapm,nMapHRMm,nMap,nTD,nTmax,nCar,JDCar,bForecast,XCbeg,xIncd,NCC,
     &			PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
     &			bDdenerb,bVdenerb,bVvererb,bDvererb,ERDLOSB,ERVLOSB,ANMAPD,ANMAPV,
     &			90.0,Scale,VMAPD,DMAP,VMHR,DMHR,XCshift3,DVfact,DDfact,V3DHR,D3DHR,
     &			DDD1,VVV3) ! BVJ uses nMapm,nMapHRMm,XCshift3, Writes DDD1, VVV3(i,j,k,n,1). VVV3(i,j,k,n,2) & VVV3(i,j,k,n,3) pass.

	      RRSCON = (R1AU/RRS)**FALLOFFD
	      do N=1,nT
	        call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N))
	      end do

	      ClipLng  = 90.0
	      nMapbm   = 1
	      nMapHRMm = nMap
C
	      xIncm = xInc
	      if(.not.bForecast) xIncm = 0.0
C
	      nTM   = nT
	      if(.not.bForecast) then
		ANinterp = Ninter + 1
		XCdif = sngl(XCintDG(2)-XCintDG(1))
	 	AI    = XCdif/ANinterp
		nTM   = (xInc*AI)/27.2753 + nT
	      end if
	      print *, ' '
	      write (*,'(A,2F5.2,I5)') 'Beginning and end MAG boundary xInc, xIncm, nTM', xInc, xIncm, nTM

C -Nit makes this a 3DMHD test write.

	      call write3D_bbtm_HR_3DMHD(0,bForecast,Nit,NiterT,NintLng,NintLat,NintHt,Ninter,NCoff,XCintDG,XCbeGG,XCtbegG,XCtendG,iYrBG,
     &			RRS,dRR,RRMS,nLng,nLat,nMapbm,nMapHRMm,nMap,nTM,nTmax,nCar,JDCar,XCbeg,xIncm,NCC,
     &	      		PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
     &			ClipLng,VMAP,XCshift3,XCshiftM,DVfact,
     &			BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,FALLOFFBN,iBFlag,BRR2DT,BRC2DT,BTC2DT,BNC2DT,BR3DHR,BT3DHR,BN3DHR,
     &			BBB3) 	! Uses RRMS nMapHRm, nMapbm, XCshift3, XCshiftM, Writes BBB3

	      if(NiT.lt.0) Nit = -Nit   ! If Nit was set negative to provide only an internal output to the tomography program, set positive.

C	      do L=1,nTm
C		XCEA(L) = XMAP_SC_POS8(EARTH8,iYrBG,XCintDG(L),nCAR,JDCar)	! Earth location in Carrington coords. at DOY XCintDG
C	      end do


	      do N=1,nT
		NBDDD1 = 0
		NBTTT1 = 0
		NBVVV1 = 0
		NBVVV2 = 0
		NBVVV3 = 0
		NBBBB1 = 0
		NBBBB2 = 0
	 	NBBBB3 = 0
		NBDLL1 = 0
		NBVLL1 = 0
		NBVLL2 = 0
		NBVLL3 = 0
		NBBLL1 = 0
		NBBLL2 = 0
		NBBLL3 = 0
		NBTLL1 = 0
		if(N.eq.1) then
		  print *, ' '
		  print *, 'In externalrwmhd.f. After INTERNAL write. By time, values are checked.'
	        end if
		nLngLatM = nLngLat*nMap
		do M=1,nMap
		  do J=1,nLat
		    do I=1,nLng
		      VVV31(I,J) = VVV3(I,J,M,N,1)
		      VVV32(I,J) = VVV3(I,J,M,N,2)
		      VVV33(I,J) = VVV3(I,J,M,N,3)
		      BBB31(I,J) = BBB3(I,J,M,N,1)
		      BBB32(I,J) = BBB3(I,J,M,N,2)
		      BBB33(I,J) = BBB3(I,J,M,N,3)
		      DDD11(I,J) = DDD1(I,J,M,N)
		      TTT11(I,J) = TTT1(I,J,M,N)
		    end do
		  end do
		  call arrR4getminmax(nLngLat,DDD11,amin,amax)
		  if(amax.ne.BadD) then 
		    do J=1,nLat
		      do I=1,nLng
			if(DDD11(I,J).eq.BadD) NBDDD1 = NBDDD1 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,DDD11,CONRD,3,0.0,0.0)		! D Should have no holes
		  else
		    NBDLL1 = NBDLL1 + nLngLat
		  end if
		  call arrR4getminmax(nLngLat,VVV31,amin,amax)
		  if(amax.ne.BadV) then 
		    do J=1,nLat
		      do I=1,nLng
			if(VVV31(I,J).eq.BadV) NBVVV1 = NBVVV1 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,VVV31,CONRD,3,0.0,0.0)		! Vr should have no holes if any valid

		  else
		    NBVLL1 = NBVLL1 + nLngLat

		  end if

		  call arrR4getminmax(nLngLat,VVV32,amin,amax)
		  if(amax.ne.BadV) then 
		    do J=1,nLat
		      do I=1,nLng
			if(VVV32(I,J).eq.BadV) NBVVV2 = NBVVV2 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,VVV32,CONRD,3,0.0,0.0)		! Vt should have no holes if any valid
		  else
		    NBVLL2 = NBVLL2 + nLngLat
		  end if

		  call arrR4getminmax(nLngLat,VVV33,amin,amax)
		  if(amax.ne.BadV) then 
		    do J=1,nLat
		      do I=1,nLng
			if(VVV33(I,J).eq.BadV) NBVVV3 = NBVVV3 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,VVV33,CONRD,3,0.0,0.0)		! Vn should have no holes if any valid
		  else
		    NBVLL3 = NBVLL3 + nLngLat
		  end if

		  call arrR4getminmax(nLngLat,BBB31,amin,amax)
		  if(amax.ne.BadB) then 
		    do J=1,nLat
		      do I=1,nLng
			if(BBB31(I,J).eq.BadB) NBBBB1 = NBBBB1 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,BBB31,CONRD,3,0.0,0.0)		! Br should have no holes if any valid
		  else
		    NBBLL1 = NBBLL1 + nLngLat
		  end if

		  call arrR4getminmax(nLngLat,BBB32,amin,amax)
		  if(amax.ne.BadB) then 
		    do J=1,nLat
		      do I=1,nLng
			if(BBB32(I,J).eq.BadB) NBBBB2 = NBBBB2 + 1
C			if(BBB32(I,J).eq.BadB) print *, 'I =',I,' J =',J,' M =',M,' N =',N,' BadB'
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,BBB32,CONRD,3,0.0,0.0)		! Bt should have no holes if any valid
		  else
		    NBBLL2 = NBBLL2 + nLngLat
		  end if

		  call arrR4getminmax(nLngLat,BBB33,amin,amax)
		  if(amax.ne.BadB) then 
		    do J=1,nLat
		      do I=1,nLng
			if(BBB33(I,J).eq.BadB) NBBBB3 = NBBBB3 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,BBB33,CONRD,3,0.0,0.0)		! Bn should have no holes if any valid
		  else
		    NBBLL3 = NBBLL3 + nLngLat
		  end if

		  call arrR4getminmax(nLngLat,TTT11,amin,amax)
		  if(amax.ne.BadT) then 
		    do J=1,nLat
		      do I=1,nLng
			if(TTT11(I,J).eq.BadT) NBTTT1 = NBTTT1 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,TTT11,CONRD,3,0.0,0.0)		! T should have no holes if any valid
		  else
		    NBTLL1 = NBTLL + nLngLat
		    if(M.eq.1.and.N.eq.1) print *, 'All temperatures set to 10000.0 deg (from 1 AU).' 
		    call ArrR4Constant(nLngLat,10000.0,TTT11)				! ,or T Should be 10000.0 at 1 AU
		  end if

		  do J=1,nLat
		    do I=1,nLng
		      VVV3(I,J,M,N,1) = VVV31(I,J)
		      VVV3(I,J,M,N,2) = VVV32(I,J)
		      VVV3(I,J,M,N,3) = VVV33(I,J)
		      BBB3(I,J,M,N,1) = BBB31(I,J)
		      BBB3(I,J,M,N,1) = BBB32(I,J)
		      BBB3(I,J,M,N,1) = BBB33(I,J)
		      DDD1(I,J,M,N)   = DDD11(I,J)
		      TTT1(I,J,M,N)   = TTT11(I,J)	      
		    end do
		  end do

		end do										! end distance do loop

		do M=1,nMap
		  do J=1,nLat
		    do I=1,nLng
		      VVVV32(I,J,M) = VVV3(I,J,M,N,2)
		      VVVV33(I,J,M) = VVV3(I,J,M,N,3)
		      BBBB31(I,J,M) = BBB3(I,J,M,N,1)
		      BBBB32(I,J,M) = BBB3(I,J,M,N,2)
		      BBBB33(I,J,M) = BBB3(I,J,M,N,3)
		    end do
		  end do
		end do

		if((NBDLL1+NBVLL1+NBVLL2+NBVLL3+NBBLL1+NBBLL2+NBBLL3+NBTLL1).ne.0) print *, ' '
		if(NBDLL1.eq.nLngLatM) write(*,'(A,I3,A,I6,A)') 'N =',N,'  All', NBDLL1,' density points were bad this time'
		if(NBVLL1.eq.nLngLatM) write(*,'(A,I3,A,I6,A)') 'N =',N,'  All', NBVLL1,' velocity radial points were bad this time'
		if(NBVLL2.eq.nLngLatM) then
		  write(*,'(A,I3,A,I6,A)') 'N =',N,'  All', NBVLL2,
     &		  ' velocity tangential points were bad this time. Set to 0.0'
		  call ArrR4Constant(nLngLatM,0.0,VVVV32)				! All Vt bad. Set to 0.0
		end if
		if(NBVLL3.eq.nLngLatM) then
	 	  write(*,'(A,I3,A,I6,A)') 'N =',N,'  All', NBVLL3,
     &		  ' velocity normal points were bad this time. Set to 0.0'
		  call ArrR4Constant(nLngLatM,0.0,VVVV33)				! All Vn bad. Set to 0.0
		end if
		if(NBBLL1.eq.nLngLatM) then
		  write(*,'(A,I3,A,I6,A)') 'N =',N,'  All', NBBLL1,
     &		  ' field radial points were bad this time. Set to 0.0'
		  call ArrR4Constant(nLngLatM,0.0,BBBB31)				! All Br bad. Set to 0.0
		end if
		if(NBBLL2.eq.nLngLatM) then
		  write(*,'(A,I3,A,I6,A)') 'N =',N,'  All', NBBLL2,
     &		  ' field tangential points were bad this time. Set to 0.0'
		  call ArrR4Constant(nLngLatM,0.0,BBBB32)				! All Bt bad. Set to 0.0
		end if
		if(NBBLL3.eq.nLngLatM) then
		  write(*,'(A,I3,A,I6,A)') 'N =',N,'  All', NBBLL3,
     &		  ' field normal points were bad this time. Set to 0.0'
		  call ArrR4Constant(nLngLatM,0.0,BBBB33)				! All Bn bad. Set to 0.0
		end if

		do M=1,nMap
		  do J=1,nLat
		    do I=1,nLng
		      VVV3(I,J,M,N,2) = VVVV32(I,J,M) 
		      VVV3(I,J,M,N,3) = VVVV33(I,J,M)
		      BBB3(I,J,M,N,1) = BBBB31(I,J,M)
		      BBB3(I,J,M,N,2) = BBBB32(I,J,M) 
		      BBB3(I,J,M,N,3) = BBBB33(I,J,M)
		    end do
		  end do
		end do

		if(NBTLL1.eq.nLngLatM) write(*,'(A,I3,A,I6,A)') 'N =',N,'  All', NBTLL1,' temperature points were bad this time'

		if((NBDDD1+NBVVV1+NBVVV2+NBVVV3+NBBBB1+NBBBB2+NBBBB3+NBTTT1).ne.0) print *, ' '
	        if(NBDDD1.ne.0) write(*,'(A,I3,A,I6,A,I6,A)') 'N =',N,' There were',NBDDD1,' bad points of',nLngLatM,
     &		  ' this density time'
	        if(NBVVV1.ne.0) write(*,'(A,I3,A,I6,A,I6,A)') 'N =',N,' There were',NBVVV1,' bad points of',nLngLatM,
     &		  ' this velocity radial time'
	        if(NBVVV2.ne.0) write(*,'(A,I3,A,I6,A,I6,A)') 'N =',N,' There were',NBVVV2,' bad points of',nLngLatM,
     &		  ' this velocity tangential time'
	        if(NBVVV3.ne.0) write(*,'(A,I3,A,I6,A,I6,A)') 'N =',N,' There were',NBVVV3,' bad points of',nLngLatM,
     &		  ' this velocity normal time'
	        if(NBBBB1.ne.0) write(*,'(A,I3,A,I6,A,I6,A)') 'N =',N,' There were',NBBBB1,' bad points of',nLngLatM,
     &		  ' this field radial time'
	        if(NBBBB2.ne.0) write(*,'(A,I3,A,I6,A,I6,A)') 'N =',N,' There were',NBBBB2,' bad points of',nLngLatM,
     &		  ' this field tangential time'
	        if(NBBBB3.ne.0) write(*,'(A,I3,A,I6,A,I6,A)') 'N =',N,' There were',NBBBB3,' bad points of',nLngLatM,
     &		  ' this field normal time'
	        if(NBTTT1.ne.0) write(*,'(A,I3,A,I6,A,I6,A)') 'N =',N,' There were',NBTTT1,' bad points of',nLngLatM,
     &		  ' this temperature time'
	      end do										! end time do loop

9999	      continue

	      nTm = nT	  
	      do L=1,nTm
		XCEA(L) = XMAP_SC_POS8(EARTH8,iYrBG,XCintDG(L),nCAR,JDCar)	! Earth location in Carrington coords. at DOY XCintDG
	      end do

	      if(iYes .eq. 1 .or. iYes .eq. 3) then 								! Write out files

	        print *, 'Before call ExternalWrite', LLfst, LLend, nT

	        call ExternalWrite(MHDs,nLng,nLat,nMap,nT,nTmax,LLfst,LLEnd,JDCar,nCar,XCEA,NCOFF,iYrBG,XCintDG,
     &			RRS,dRR,FALLOFFD,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,DDD1,TTT1,VVV3,BBB3)
	      end if 
C

	      if(iYes .eq. 2 .or. iYes .eq. 3) then 								! Read in files
		LLfstt = LLfst
		LLEndd = LLend

	        print *, 'Before call ExternalRead', LLfst, LLend, nT

	        call ExternalRead(cWild3DMHD,MHDs,nLng,nLat,nMap,nT,nTmax,LLfstt,LLEndd,JDCar,nCar,XCEA,NCOFF,iyrBG,XCintDG,	
     &			RRS,dRR,FALLOFFD,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,DDD1,TTT1,VVV3,BBB3,iNum)
C
C		Now we need to fill the spatial holes in the input time series.
C
		ifill = 0
C		ifill = 1
		if(ifill.eq.1) then
		  call Fill_in(0 ,0,NiterT,XCbegg,nLng,nLat,nMap,nT,nTmax,ConsTMD,ConstL,ALng,
     &			CONRD,CONDT,aNdayG,DDD1,VVV3)	! Density
		  call Fill_in(1 ,0,NiterT,XCbegg,nLng,nLat,nMap,nT,nTmax,ConsTMV,ConstL,ALng,
     &			CONRV,CONVT,aNdayV,DDD1,VVV3)	! Velocity
		  call Fill_in(30,0,NiterT,XCbegg,nLng,nLat,nMap,nT,nTmax,ConsTMV,ConstL,ALng,
     &			CONRV,CONVT,aNdayV,DDD1,BBB3)	! Magnetic
		  call Fill_in(10,0,NiterT,XCbegg,nLng,nLat,nMap,nT,nTmax,ConsTMD,ConstL,ALng,
     &			CONRD,CONDT,aNdayG,TTT1,VVV3)	! Temperature
		end if
C
C		Now shift all the external inputs to the internal arrays, and reset DVfact and DDfact to 1.0
C
		do L=1,nT
		  do I=1,nLng
	 	    do J=1,nLat
		      do K=1,nMap

		        DD1   (I,J,K,L)   = DDD1(I,J,K,L)
		        VV3   (I,J,K,L,1) = VVV3(I,J,K,L,1)
		        VV3   (I,J,K,L,2) = VVV3(I,J,K,L,2)
		        VV3   (I,J,K,L,3) = VVV3(I,J,K,L,3)
		        BB3   (I,J,K,L,1) = BBB3(I,J,K,L,1)
		        BB3   (I,J,K,L,2) = BBB3(I,J,K,L,2)
		        BB3   (I,J,K,L,3) = BBB3(I,J,K,L,3)
		        TT1   (I,J,K,L)   = TTT1(I,J,K,L)
		        DVfact(I,J,K,L)   = 1.0	
		        DDfact(I,J,K,L)   = 1.0
C		        if(i.eq.25.and.j.eq.5) write(*,'(4I3,F10.4,3F7.2,3F8.2,F11.1)') i,j,k,L, DD1(i,j,k,L),
C     &			VV3(i,j,k,L,1),VV3(i,j,k,L,2),VV3(i,j,k,L,2),BB3(i,j,k,L,1),BB3(i,j,k,L,2), BB3(i,j,k,L,3), TT1(I,J,K,L)

		      end do
		    end do
		  end do
		end do

	      end if

	      print *, ' '
	      print *, 'End of Externalrwmhd.f'
	      print *, ' '

	 return
	end