C+
C NAME:
C	ExternalWtest.f
C PURPOSE:
C	Set up to write the externally-provided volumetric elements provided by the UCSD 3D-MHD program in order to read 
C	them as a test. 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	Thes 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 ExternalWtest(MHDs,bForecast,nLng,nLat,nMap,nT,nTmax,XCbegg,XCbeg,xInc,NCC,LLfst,LLEnd,JDCar,nCar,
C     &				 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)
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	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,nT)		real		Input 1 component density
C	TTT1(nLng,nLat,NMap,nT)		real		Input 1 component temperature
C	VVV3(nLng,nLat,NMap,nT,3)	real		Input 3 component velocity (r,t,n)
C	BBB3(nLng,nLat,NMap,nT,3)	real		Input 3 component magnetic field (r,t,n)
C
C MODIFICATION HISTORY:
C	September -2016, Bernard Jackson (UCSD)
C-
	subroutine ExternalWtest(MHDs,bForecast,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,
     &		DDD1,TTT1,VVV3,BBB3)

	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),
     &		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 write test ExternalWtest.f'
	Nit = -Nit						! Triggers an internal write
	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
	      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 externalwtest.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
	      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(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

	      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

	      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 externalwtest.f. After INTERNAL write. By time, values are checked.'
	        end if
		nLngLatM = nLngLat*nMap
		do M=1,nMap
		  call arrR4getminmax(nLngLat,DDD1(1,1,M,N),amin,amax)
		  if(amax.ne.BadD) then 
		    do J=1,nLat
		      do I=1,nLng
			if(DDD1(I,J,M,N).eq.BadD) 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
		    NBDLL1 = NBDLL1 + nLngLat
		  end if

		  call arrR4getminmax(nLngLat,VVV3(1,1,M,N,1),amin,amax)
		  if(amax.ne.BadV) then 
		    do J=1,nLat
		      do I=1,nLng
			if(VVV3(I,J,M,N,1).eq.BadV) 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 if any valid
		  else
		    NBVLL1 = NBVLL1 + nLngLat
		  end if

		  call arrR4getminmax(nLngLat,VVV3(1,1,M,N,2),amin,amax)
		  if(amax.ne.BadV) then 
		    do J=1,nLat
		      do I=1,nLng
			if(VVV3(I,J,M,N,2).eq.BadV) 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 if any valid
		  else
		    NBVLL2 = NBVLL2 + nLngLat
		  end if

		  call arrR4getminmax(nLngLat,VVV3(1,1,M,N,3),amin,amax)
		  if(amax.ne.BadV) then 
		    do J=1,nLat
		      do I=1,nLng
			if(VVV3(I,J,M,N,3).eq.BadV) 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 if any valid
		  else
		    NBVLL3 = NBVLL3 + nLngLat
		  end if

		  call arrR4getminmax(nLngLat,BBB3(1,1,M,N,1),amin,amax)
		  if(amax.ne.BadB) then 
		    do J=1,nLat
		      do I=1,nLng
			if(BBB3(I,J,M,N,1).eq.BadB) 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 if any valid
		  else
		    NBBLL1 = NBBLL1 + nLngLat
		  end if

		  call arrR4getminmax(nLngLat,BBB3(1,1,M,N,2),amin,amax)
		  if(amax.ne.BadB) then 
		    do J=1,nLat
		      do I=1,nLng
			if(BBB3(I,J,M,N,2).eq.BadB) NBBBB2 = NBBBB2 + 1
C			if(BBB3(I,J,M,N,2).eq.BadB) print *, 'I =',I,' J =',J,' M =',M,' N =',N,' BadB'
		      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 if any valid
		  else
		    NBBLL2 = NBBLL2 + nLngLat
		  end if

		  call arrR4getminmax(nLngLat,BBB3(1,1,M,N,3),amin,amax)
		  if(amax.ne.BadB) then 
		    do J=1,nLat
		      do I=1,nLng
			if(BBB3(I,J,M,N,3).eq.BadB) 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 if any valid
		  else
		    NBBLL3 = NBBLL3 + nLngLat
		  end if

		  call arrR4getminmax(nLngLat,TTT1(1,1,M,N),amin,amax)
		  if(amax.ne.BadT) then 
		    do J=1,nLat
		      do I=1,nLng
			if(TTT1(I,J,M,N).eq.BadT) 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 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,TTT1(1,1,M,N))				! ,or T Should be 10000.0 at 1 AU
		  end if

		end do										! end distance do loop

		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,VVV3(1,1,1,N,2))				! 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,VVV3(1,1,1,N,3))				! 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,BBB3(1,1,1,N,1))				! 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,BBB3(1,1,1,N,2))				! 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,BBB3(1,1,1,N,3))				! All Bn bad. Set to 0.0
		end if
		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
	print *, 'Before call ExternalWrite', LLfst, LLend
	      if(NiT .lt. 0) then 	! Change back to .lt. 0 when the above works
	        call ExternalWrite(MHDs,nLng,nLat,nMap,nT,LLfst,LLEnd,JDCar,nCar,XCEA,NCOFF,iYrBG,XCintDG,
     &			RRS,dRR,FALLOFFD,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,DDD1,TTT1,VVV3,BBB3)
	      end if 

	    if(NiT.lt.0) Nit = -Nit

	  print *, 'End of external write test ExternalWtest.f'
	  print *, ' '

	return
	end