C+
C NAME:
C	write3D_infotd3D
C PURPOSE:
C	Writes a 3D file and pertinent information at times
C	XCint(N) + half of the interval. (B Jackson, Mar., 2001)
C
C  XCint(N) values are the beginnings (and ends) of the time intervals used to select data.  When the data selection
C  used is set to a daily cadence with a one day resolution, these XCint(N) values are the beginning and 
C  ending midnight on Earth at the longitude the data is obtained.  This may be different if one site is 
C  used for density and one is used for velocity.  In any case, the inputing XCint(N) values rule in the 
C  selection of the matrix to be output from this routine. The actual start time of the time matrix XCtbeg 
C  is set to the middle of the first set of these values and proceeds at precise nT interval values from this 
C  value until the end value XCtbeg, which is the middle of the last interval. Thus, the time value for the 
C  extraction of the matrix as near as possible should be the middle of the temporal XCint(N) intervals in
C  order to match the Earth location at noon.  The motive here is to select an extract time when most of the 
C  IPS sources are observed.  If temporal intervals are set at a daily cadence, but at higher temporal 
C  resolutions than one day, then the matrix output, should still be set to the middle of the interval. 
C
C CALLING SEQUENCE:
	subroutine write3D_infotd3D(Mo,cPre,Nit,NiterT,Ninter,NCoff,XCintD,XCbe,XCtbeg,XCtend,iYr,RR,dRR,
     &		nLng,nLat,nMap,nT,nTmax,nCar,JDCar,
     &		PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
     &		ClipLng,Scale,VMAP,DMAP,VV,DD,XCshift,Vratio,Dratio,V3D,D3D,LLbeg,LLend,XCTT,XC3DT,V3DT,D3DT)
C INPUTS:
C	Mo		integer		Mode of operation 0 - no additional 1/R beyond 1 AU written
C							  1 - an additional 1/R beyond 1 AU written
C							  2
C	cPre		character	Character prefix for output files
C	Nit		integer 	interation number
C	NiterT		integer		maximum number of iterations
C	Ninter		integer		# of intermediate steps
C	NCoff		integer		Carrington variable offset
C	XCintD(nTmax)	double precision Beginning and ends of the time intervals used to select data
C	XCbe(2,nT)	real		Beginning (1,nT) and end (2,nT) of Carrington rotation times 
C	XCtbeg		real		Start time
C	XCtend		real		End time
C	iYr		integer		year at beginning time
C	RR		real		radial distance of the source surface
C	dRR		real		radial resolution
C	nLng		integer		# longitudes
C	nLat		integer		# latitudes
C	nMap		integer		# radial distances
C	nT		integer		# of time values
C	nTmax		integer		maximum # of time values
C	nCar		integer		# Carrington rotations
C	JDCar		double precision Julian Date at beginning of rotations
C	PWV		real		g-level power to fit velocity
C	PWG		real		g-level power to fit density
C	PWRV		real		radial g-level power to fit g density 
C					(or radial density falloff for Thomson scattering)
C	PWRG		real		radial g-level power to fit V density
C	DEN1AU		real		Density at 1 AU
C	CONRV		real		Velocity spatial smoothing constant
C	CONRD		real		Density spatial smoothing constant
C	CONSTV		real		Velocity spatial hole-filling constant
C	CONSTG		real		Density spatial hole-filling constant
C	CONVT		real		Velocity temporal smoothing constant
C	CONDT		real		Density temporal smoothing constant
C	ClipLng		real		Clip longitude
C	Scale		real		extra scaling factor (set to 1 if no scaling is required)
C	VMAP(nLng,nLat)	real		Velocity map at source surface
C	DMAP(nLng,nLat)	real		Density map at Source Surface
C	VV(nLng,nLat)	real		(scratch array) Velocity map at source surface
C	DD(nLng,nLat)	real		(scratch array) Density map at source surface
C       XCshift(nLng,nLat,nMap,nT,3) real shift amount from array
C       Vratio(nLng,nLat,nMap,nT)  real velocity difference at height
C       Dratio(nLng,nLat,nMap,nT)  real density difference at height
C	V3D(nLng,nLat,nMap)	real	3D velocity matrix (dummy array) written out
C	D3D(nLng,nLat,nMap)	real	3D density matrix (dummy array) written out - normalized to 1 AU
C	iStrAdd		integer		integer number of strings to be added to the file header
C	cfile(*)*(*)	character	strings to be added to the file header
C OUTPUTS:
C	LLbeg		integer		beginning # of volumes to be exported
C	LLend		integer		end # of volumes to be exported
C	XCTT  ((Ninter+1)*nT)	real	times (evenly spaced in UT)
C	XC3DT(nLng,nLat,nMap,(Ninter+1)*nT,3)	real	shifts needed to go back to source surface
C	V3DT(nLng,nLat,nMap,(Ninter+1)*nT) real 3 D velocity over time	
C	D3DT(nLng,nLat,nMap,(Ninter+1)*nT) real Normalized 3 D density over time
C CALLS:
C       T3D_fill_global, T3D_fill_time, T3D_write_nv
C-
	character	cStrAdd(1)*120
	character	cFile*120
	character	cPre*4

	real		XCbe(2,nTmax)

	real		VMAP	(nLng,nLat,nT)		! Velocity source surface map
	real		DMAP	(nLng,nLat,nT)		! Density source surface map
	real		VV	(nLng,nLat)		! Velocity source surface map written 
	real		DD	(nLng,nLat)		! Density source surface map written 
	real		XCshift	(nLng,nLat,nMap,nT,3)
	real		Vratio	(nLng,nLat,nMap,nT)
	real		Dratio	(nLng,nLat,nMap,nT)
	real		V3D	(nLng,nLat,nMap)
	real		D3D	(nLng,nLat,nMap)
	real		V3DT	(nLng,nLat,nMap,(Ninter+1)*nT)
	real		D3DT	(nLng,nLat,nMap,(Ninter+1)*nT)
	real		Scale	(2)
	real		XCsh	(3)

	real		XCTT	((Ninter+1)*nT)		! Contains the carrington-variable times that
							! correspond to the 6-day time-steps.
	real 		XC3DT   (nLng,nLat,nMap,(Ninter+1)*nT,3)!Shift matrix.
							!First 4 dim. refer to position.
							!Last dim stores content, eg,
							!shift-information
							!(xclng,xclat,xctim),	 for that point
	real		PwrR (2)		/0.0, 2.0/

	double precision XCtbeg,XCtend,XCtfull,X,XCtim,XCt,JDCar(nCar)
	double precision XCintD(nTmax)

	external	EARTH8

	include		't3d_array.h'
	include		'mapcoordinates.h'

	Bad = BadR4()

	I = nLng*nLat*nMap*(Ninter+1)*nT
	call ArrR4Bad(I*3,XC3DT)
	call ArrR4Bad(I  ,D3DT )
	call ArrR4Bad(I  ,V3DT )	

	! print *, 'In write3d_infotd'
	
	XCdif  = sngl(XCtend-XCtbeg)/nT
	XCbegt = sngl(XCtbeg)

	if(mo .ge. 0) call T3D_fill_global(t3d,0,Nit,NiterT,NCoff,XCbegt,XCdif,RR,dRR,nLng,nLat,nMap,nT,
     &	      			PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,ClipLng,Scale)
     	call T3D_set(T3D__R_PWR,2,PwrR )
	iset =1
	Ninterp  = Ninter+1
	ANinterp = Ninterp
	RRSDONE2 = RR*RR

	LLbeg = 0
	do N=1,nT

	    XCbeg = XCbe(1,N)
	    XCend = XCbe(2,N)
	    XCdif = sngl(XCintD(N+1)-XCintD(N))
	    AI    = XCdif/ANinterp

	    do L=1,Ninterp

		XCtim = XCintD(N)+dble(AI*(L-1))		! Times the matrix is to be extracted this N

		XCM = XMAP_SC_POS8(EARTH8,iYr,XCtim,nCAR,JDCar)	! Earth location in Carrington coords. at XCtim

		XCbegROI = XCM-0.5 				! Region of interest beginning (-NCoff)
		XCendROI = XCM+0.5				! Region of interest ending (-NCoff)

		TT = XCM
		LL = Ninterp*N-Ninterp+L
								! Store this value in XCTT to pass out of the subroutine
		if (TT .gt. XCbeg+0.75 .and. TT .lt. XCbeg+2.0) XCTT(LL) = TT

		if (TT .gt. XCbeg+1.00 .and. TT .lt. XCbeg+2.0) then
		    if (LLbeg .eq. 0) LLbeg = LL	   
		    LLend = LL
	    	  
!$omp parallel private(I,J,K,RRht,RRSDRR2,XLat,XClng,XCsh,Drat,Vrat,
!$omp&  XC,XLa,XCt,DEN,VEL)
!$omp do schedule(static)
		    do K=1,nMap
			RRht = RRhght(K)			! Height the matrix is to be extracted
			RRSDRR2 = (RRhght(1)/RRht)**2
			do J=1,nLat
			    XLat = XLdeg(J)			! Latitude the matrix is to be extracted
			    do I=1,nLng

				XClng = XCfull(I)		! Longitude the matrix is to be extracted

				call Get4Dval(3,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat,nMap,nT,
     &					XCshift,1,XCtim,XClng,XLat,RRht,XCsh)

				call Get4Dval(1,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat,nMap,nT,
     &					Dratio,1,XCtim,XClng,XLat,RRht,Drat)

				call Get4Dval(1,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat,nMap,nT,
     &					Vratio,1,XCtim,XClng,XLat,RRht,Vrat)

				! print *, 'In write3d_infotd after Get4Dvals'

				XC3DT(I,J,K,LL,1) = XCsh(1)	! storing this value in XC3DT to pass out of the subroutine
				XC3DT(I,J,K,LL,2) = XCsh(2)	! storing this value in XC3DT to pass out of the subroutine
				XC3DT(I,J,K,LL,3) = XCsh(1)	! storing this value in XC3DT to pass out of the subroutine

				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))

				call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,DMAP,1,XC,XLa,XCt,DEN)
				call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,VMAP,1,XC,XLa,XCt,VEL)

				! print *, 'In write3d_infotd after Get3DTvals'

				if (DEN .eq. Bad) then
				    D3D (I,J,K   ) = Bad
				else
				    D3D (I,J,K   ) = Drat*DEN*RRSDONE2
				    D3DT(I,J,K,LL) = Drat*DEN*RRSDRR2
				end if

				if (VEL .eq. Bad) then
				    V3D (I,J,K   ) = Bad
				else
				    V3D (I,J,K   ) = Vrat*VEL
				    V3DT(I,J,K,LL) = Vrat*VEL
				end if

				if (abs(Mo) .eq. 1) then
				    ! If Mo eq 1 remove an additional 1/R density beyond
				    ! 1 AU with the next 3 lines in the file output by
				    ! the write3d subroutine

				    if (K .gt. 11) then
					if (DEN .ne. Bad) then
					    D3D(I,J,K) = D3D(I,J,K)*(10.0)/(K-1.0)
					end if
				    end if
				end if

			    end do	! do I=1,nLng
			end do		! do J=1,nLat
		    end do		! do K=1,nMap
!$omp end do
!$omp do schedule(static)
		    do J=1,nLat
			do I=1,nLng
			    DD(I,J) = D3D(I,J,1)
			    VV(I,J) = V3D(I,J,1)
			end do
		    end do
!$omp end do nowait
!$omp end parallel
		    ! print *, 'mo, TT, XCbeg+1, Xcbeg+2', mo, TT, XCbeg+1, Xcbeg+2

		    if (mo .ge. 0) then
			call T3D_fill_time(t3d,nT,TT,XCbeg,XCend,XCbegROI,XCendROI)
			call T3D_write_nv(cPre,t3d,VV,DD,V3D,D3D,0,cStrAdd,cFile)
		    end if

		end if

	    end do			! do L=1,Ninterp
	end do				! do N=1,nT

	return
	end