C+
C NAME:
C	ExternalRtest.f
C PURPOSE:
C	Set up to read the externally-provided volumetric elements provided by an external 3D-MHD program .
C	These volumetric elements are read from files provided that are interpolated with the cadence 
C	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 the 3D-MHD program's run location
C	by ExternalConvert 
C CATEGORY:
C	I/O
C CALLING SEQUENCE:
C	call ExternalRtest(cWild3DMHD,MHDs,nLng,nLat,nMap,nTmax,LLBegg,ISfile,ISfilef,LLLfstp,LLEnd,NinterD,IINC,IINCD2,IYRBG,XCintDG,TT3D,
C		          RRS,RADMS,RAD1,dRR,RADS,FALLOFFN,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,DDD1,TTT1,VVV3,BBB3,iYes)
C
C INPUTS:
C	cWild3DMHD			character*80	input file name	
C	MHDs				integer		Which 3D MHD program provides the input?   1 - ENLIL, 2 - MSFLUKSS, 3 - HAF3DMHD
C	nLng				integer		Number of Longitude bins	
C	nLat				integer		Number of Latitude bins
C	nMap				integer		Number of Map heights
C	nTmax				integer		Maximum number of times possible
C	LLBegg				integer		Beginning number of the used 6-hour averages
C	ISfile				integer		Beginning file number to write out
C	ISfilef				integer		Ending file number to write out
C	LLLfstp				integer		Beginning increment of daily averages
C	LLEnd				integer		Ending number of 6-hour averages
C	NinterD				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(nTmax)			real*8		Day of Year	
C	TT3D(4*nTmax)			real		Times of the data (in Carrington Variable) for high-resolution output
C	RRS				real		Beginning Height of the IPS modeling
C	RRMS				real		Beginning Height of the Magnetic Field pick-up
C	dRR				real		Outward step height of the modeling
C	FALLOFFN			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	iYes				integer		iYes = 0, check the file
C							iYes = 3, check the file and read the heights input
C							iYes = -3, test the programs with internal arrays
C
C OUTPUTS:
C	RADS(nMap+1)			real		Heights input from the model (The first height is RADMS).
C	DDD1(nLng,nLat,NMap+1,nTmax)	real		Input 1 component density
C	TTT1(nLng,nLat,NMap+1,nTmax)	real		Input 1 component temperature
C	VVV3(nLng,nLat,NMap+1,nTmax,3)	real		Input 3 component velocity (r,t,n)
C	BBB3(nLng,nLat,NMap+1,nTmax,3)	real		Input 3 component magnetic field (r,t,n)
C	iYes				integer		iYes = 0, 3 no file was found
C							iYes = 1, a correct file is found
C							iYes = 2, a new file couldn't be read
C							iYes = 3, a test file was read for header and heights alone
CC							iYes = -3, no files were read
C
C MODIFICATION HISTORY:
C	October-2016, Bernard Jackson (UCSD)
C-
	subroutine ExternalRtest(MHDs,bForecast,cWild3DMHD,nLng,nLat,nMap,nT,nTmax,XCbegg,XCbeg,xInc,NCC,
     &		LLfst,LLEndd,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,
     &		DDD1,TTT1,VVV3,BBB3) 				! Use xcshift3 and xcshiftM matrices

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

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

	real	XCshiftM(nLng,nLat,nT,3),			! these are needed before and after the external read 
     &		BRR2DT  (nLng,nLat,nT),
     &		BRC2DT  (nLng,nLat,nT),
     &		BTC2DT  (nLng,nLat,nT),
     &		BNC2DT  (nLng,nLat,nT),
     &		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,nT),					! Density LOS error source surface map
     &		ANMAPV	(nLng,nLat,nT),					! Velocity LOS error source surface map
     &		DMAP    (nLng,nLat,nT),
     &		VMAP    (nLng,nLat,nT),
     &		VMAPD   (nLng,nLat,nT),
     &		XCSHIFT3(nLng,nLat,nMap,nT,3),
     &		xcbegg  (2,nTmax),
     &		DDfact  (nLng,nLat,nMap,nT),
     &		DVfact  (nLng,nLat,nMap,nT),
     &		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 external read test ExternalRtest.f'
	Nit = -Nit						! Triggers an internal write
	print *, 'NiT =', NiT
C
C	First internally write data from the tomography model for the times not covered by the input reads.
C
	    ALng = 1.*(nLng-1)/(iLng-1)
	    NintLng = 0
	    NintLat = 0
	    NintHt  = 0
	    Ninter = 0
	    nMapm   = 1
	    nMapHRMm = nMap
	    ClipLng  = 90.0
C
	    nLngLat = nLng*nLat
	    RRSCON = (R1AU/RRS)**FALLOFFD
	    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 externalrtest.f. Before INTERNAL writes. By time, sufaces are tested.'
		end if

		  call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax)
		  if(amax.ne.BadR4()) then 
		    do J=1,nLat
		      do I=1,nLng
			if(DMAP(I,J,N).eq.Bad) 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 if any valid
		  end if

		  call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax)
		  if(amax.ne.BadR4()) then 
		    do J=1,nLat
		      do I=1,nLng
			if(VMAP(I,J,N).eq.Bad) 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 if any valid
		  end if

		  call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax)
		  if(amax.ne.BadR4()) then 
		    do J=1,nLat
		      do I=1,nLng
			if(VMAPD(I,J,N).eq.Bad) 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 if any valid
		  end if

		  call arrR4getminmax(nLngLat,BRR2DT(1,1,N),amin,amax)
		  if(amax.ne.BadR4()) then 
		    do J=1,nLat
		      do I=1,nLng
			if(BRR2DT(I,J,N).eq.Bad) 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 if any valid
		  end if

		  call arrR4getminmax(nLngLat,BRC2DT(1,1,N),amin,amax)
		  if(amax.ne.BadR4()) then 
		    do J=1,nLat
		      do I=1,nLng
			if(BRC2DT(I,J,N).eq.Bad) 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 if any valid
		  end if

		  call arrR4getminmax(nLngLat,BTC2DT(1,1,N),amin,amax)
		  if(amax.ne.BadR4()) then 
		    do J=1,nLat
		      do I=1,nLng
			if(BTC2DT(I,J,N).eq.Bad) 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 if any valid
		  end if

		  call arrR4getminmax(nLngLat,BNC2DT(1,1,N),amin,amax)
		  if(amax.ne.BadR4()) then 
		    do J=1,nLat
		      do I=1,nLng
			if(BNC2DT(I,J,N).eq.Bad) 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 if any valid
		  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

	    do ibegend=1,2						! Some day maybe be more efficient.

	      if(ibegend.eq.1) then
		LLfstt = 1
		XCbeggg = 0.0
		print *, 'First XCbegg, XCbeggg =',XCbegg, XCbeggg
		LLEndd = LLfst - 1		! When things work
C		LLEndd = nT
	      end if
	      if(ibegend.eq.2) then
		LLfstt = LLend + 1

		ANinterp = Ninter + 1
		XCdif = sngl(XCintDG(LLfstt)-XCintDG(1))
		AI    = XCdif/ANinterp
		XCbeggg = (XCbegin+xInc)-((0.5*AI)/27.2753)
		print *, 'Second XCbegg, XCbeggg =',XCbegg, XCbeggg

		LLEndd = nT
	      end if

	      RRSCON = (RRS/R1AU)**FALLOFFD
	      do N=LLfstt,LLEndd
	        call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N))
	      end do
C
	      xIncd = xInc
	      if(.not.bForecast) xIncd = 0.0
C
	      nTD   = LLEndd
	      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,XCbeggg,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.
C
	      RRSCON = (R1AU/RRS)**FALLOFFD
	      do N=1,nT
	        call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N))
	      end do

C -Nit makes this a 3DMHD test write.
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

	      call write3D_bbtm_HR_3DMHD(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
C
C Now read in the values from the files
C
	write(*,'(4I3,3F10.4)') 25,5,10,25,BBB3(25,5,10,25,1),BBB3(25,5,10,25,2),BBB3(25,5,10,25,3)

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

	      LLfstt = 1
	      LLendd = nT
	      print *, 'In externalrtest.f Before beginning and end of external read attempt', LLfstt,LLendd

	      if(NiT .lt. 0) then ! Change back to .lt. 0 when the above works
	      call ExternalRead(cWild3DMHD,MHDs,nLng,nLat,nMap,nT,LLfstt,LLEndd,JDCar,nCar,XCEA,NCOFF,iyrBG,XCintDG,			
     &				RRS,dRR,FALLOFFD,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,DDD1,TTT1,VVV3,BBB3,iYes)
	      end if

	      print *, 'In externalrtest.f After External Read'
C  
C	Now remove bad data holes smooth the values between the read data and the iterated tomography model if these exist,
C	and as needed.
C
	      do N=1,nT
		NBDDD1 = 0
		NBTTT1 = 0
		NBVVV1 = 0
		NBVVV2 = 0
		NBVVV3 = 0
		NBBBB1 = 0
		NBBBB2 = 0
	 	NBBBB3 = 0
		do M=1,nMap
		  RRht = RRS*((float(M)*dRR - dRR) + RRS)
		  RRSD = RRht**FALLOFFD
		  RRST = RRht**FALLOFFT
		  if(N.eq.1) write(*,'(A,I4,3F10.5)') 'After External Read. M, RRht, RRSD, RRST =', M, RRht, RRSD, RRST

		  call arrR4getminmax(nLngLat,DDD1(1,1,M,N),amin,amax)
		  if(amax.ne.BadR4()) then 
		    do J=1,nLat
		      do I=1,nLng
			if(DDD1(I,J,M,N).eq.Bad) NBDDD1 = NBDDD1 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,DDD1(1,1,M,N),CONRD,3,0.0,0.0)		! D Should have no holes
		  else
		    if(M.eq.1.and.N.eq.1) print *, 'All densities set to (RRht**FALLOFFD)*4.0/cm^3 (from  1 AU).' 
		    call ArrR4Constant(nLngLat,RRSD*4.0,DDD1(1,1,M,N))				! ,or D Should be 4.0/cm^3 at 1 AU
		  end if

		  call arrR4getminmax(nLngLat,TTT1(1,1,M,N),amin,amax)
		  if(amax.ne.BadR4()) then 
		    do J=1,nLat
		      do I=1,nLng
			if(TTT1(I,J,M,N).eq.Bad) NBTTT1 = NBTTT1 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,TTT1(1,1,M,N),CONRD,3,0.0,0.0)		! T Should have no holes
		  else
		    if(M.eq.1.and.N.eq.1) print *, 'All temperatures set to (RRht**FALLOFFT)*10000.0 (from 1 AU).' 
		    call ArrR4Constant(nLngLat,RRSD*10000.0,DDD1(1,1,M,N))				! ,or D Should be 4.0/cm^3 at 1 AU
		  end if

		  call arrR4getminmax(nLngLat,VVV3(1,1,M,N,1),amin,amax)
		  if(amax.ne.BadR4()) then 
		    do J=1,nLat
		      do I=1,nLng
			if(VVV3(I,J,M,N,1).eq.Bad) NBVVV1 = NBVVV1 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,VVV3(1,1,M,N,1),CONRD,3,0.0,0.0)		! Vr Should have no holes
		  else
		    if(M.eq.1.and.N.eq.1) print *, 'All radial velocities set to 400.0 km/s.' 
		    call ArrR4Constant(nLngLat,400.0,VVV3(1,1,M,N,1))				! ,or Vr Should be 400 km/s
		  end if

		  call arrR4getminmax(nLngLat,VVV3(1,1,M,N,2),amin,amax)
		  if(amax.ne.BadR4()) then 
		    do J=1,nLat
		      do I=1,nLng
			if(VVV3(I,J,M,N,2).eq.Bad) NBVVV2 = NBVVV2 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,VVV3(1,1,M,N,2),CONRD,3,0.0,0.0)		! Vt Should have no holes
		  else
		    if(M.eq.1.and.N.eq.1) print *, 'All tangential velocities set to 0.0 km/s.' 
		    call ArrR4Constant(nLngLat,0.0,VVV3(1,1,M,N,2))				! ,or Vt Should be zero
		  end if
		  call arrR4getminmax(nLngLat,VVV3(1,1,M,N,3),amin,amax)
		  if(amax.ne.BadR4()) then 
		    do J=1,nLat
		      do I=1,nLng
			if(VVV3(I,J,M,N,3).eq.Bad) NBVVV3 = NBVVV3 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,VVV3(1,1,M,N,3),CONRD,3,0.0,0.0)		! Vn Should have no holes
		  else
		    if(M.eq.1.and.N.eq.1) print *, 'All normal velocities set to 0.0 km/s.' 
		    call ArrR4Constant(nLngLat,0.0,VVV3(1,1,M,N,3))				! ,or Vn should be be zero
		  end if

		  call arrR4getminmax(nLngLat,BBB3(1,1,M,N,1),amin,amax)
		  if(amax.ne.BadR4()) then 
		    do J=1,nLat
		      do I=1,nLng
			if(BBB3(I,J,M,N,1).eq.Bad) NBBBB1 = NBBBB1 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,BBB3(1,1,M,N,1),CONRD,3,0.0,0.0)		! Br Should have no holes
		  else
		    if(M.eq.1.and.N.eq.1) print *, 'All radial fields set to 0.0 nT.'
		    call ArrR4Constant(nLngLat,0.0,BBB3(1,1,M,N,1))				! ,or Br should be zero
		  end if
		  call arrR4getminmax(nLngLat,BBB3(1,1,M,N,2),amin,amax)
		  if(amax.ne.BadR4()) then 
		    do J=1,nLat
		      do I=1,nLng
			if(BBB3(I,J,M,N,2).eq.Bad) NBBBB2 = NBBBB2 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,BBB3(1,1,M,N,2),CONRD,3,0.0,0.0)		! Bt Should have no holes
		  else
		    if(M.eq.1.and.N.eq.1) print *, 'All tangential fields set to 0.0 nT.'
		    call ArrR4Constant(nLngLat,0.0,BBB3(1,1,M,N,2))				! ,or Bt should be zero
		  end if
		  call arrR4getminmax(nLngLat,BBB3(1,1,M,N,3),amin,amax)
		  if(amax.ne.BadR4()) then 
		    do J=1,nLat
		      do I=1,nLng
			if(BBB3(I,J,M,N,3).eq.Bad) NBBBB3 = NBBBB3 + 1
		      end do
		    end do
		    call GridSphere2D(ALng,nLng,nLat,1,BBB3(1,1,M,N,3),CONRD,3,0.0,0.0)		! Bn Should have no holes
		  else
		    if(M.eq.1.and.N.eq.1) print *, 'All normal fields set to 0.0 nT.'
		    call ArrR4Constant(nLngLat,0.0,BBB3(1,1,M,N,3))				! ,or Bn should be zero
		  end if

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

	    end do						! some day, may be more efficient
	  print *, 'End of external read test ExternalRtest.f'
	  print *, ' '

	stop

	return
	end