C+
C NAME:
C       write3D_infotd3D
C PURPOSE:
C  Writes a 3D file and pertinent information over at times 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:
C       call write3D_infotd3D(Mo,cPre,Nit,NiterT,Ninter,NCoff,XCintD,XCbe,XCtbeg,XCtend,iYr,RR,dRR,
C			nLng,nLat,nMap,nT,nTmax,nCar,JDCar,
C	      		PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
C			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)	real*8		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		real*8		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		Written velocity map at source surface
C	DD(nLng,nLat)	real		Written 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 FUNCTIONS/SUBROUTINES:
C       T3D_fill_global, T3D_fill_time, T3D_write_nv
C-

	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)

	character	cStrAdd(1)*120,
     &			cFile*120,
     &			cPre*4

	real		XCbe(2,nTmax)

	real		VMAP	(nLng,nLat,nT),		! Velocity source surface map
     &			DMAP	(nLng,nLat,nT),		! Density source surface map
     &			VV	(nLng,nLat),		! Velocity source surface map written 
     &			DD	(nLng,nLat),		! Density source surface map written 
     &			XCshift	(nLng,nLat,nMap,nT,3),
     &			Vratio	(nLng,nLat,nMap,nT),
     &			Dratio	(nLng,nLat,nMap,nT),
     &			V3D	(nLng,nLat,nMap),
     &			D3D	(nLng,nLat,nMap),
     &			V3DT	(nLng,nLat,nMap,(Ninter+1)*nT),
     &			D3DT	(nLng,nLat,nMap,(Ninter+1)*nT),
     &			Scale	(2),
     &			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/
     
	real*8		XCtbeg,XCtend,XCtfull,X,XCtim,XCt,JDCar(nCar)
	
	real*8		XCintD(nTmax)
	
	external	EARTH8

	include		't3d_array.h'
	include		'MAPCOORDINATES.H'

	Bad = BadR4()
	
	do L=1,(Ninter+1)*nT
	    XCTT(L) = Bad
	    do K=1,nMap
		do J=1,nLat
		    do I=1, nLng
			XC3DT(I,J,K,L,1) = Bad
			XC3DT(I,J,K,L,2) = Bad
			XC3DT(I,J,K,L,3) = Bad
			D3DT(I,J,K,L)    = Bad
			V3DT(I,J,K,L)    = Bad	
		    end do
		end do
	    end do
	end do

C	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
	   
	   if(TT .gt. (XCbeg+0.75) .and. TT .lt. (XCbeg+2)) then
	     XCTT(LL) = TT			! Store this value in XCTT to pass out of the subroutine
	   end if
	   
	   if(TT .gt. (XCbeg+1) .and. TT .lt. (XCbeg+2)) then
	    if(LLbeg.eq.0) LLbeg = LL
	    LLend = LL	   

	    	  
	    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)
     
C	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)
		  
C	print *, 'In write3d_infotd after Get3DTvals'

		  D3D(I,J,K) = Bad
		  V3D(I,J,K) = Bad
		  if(DEN.ne.Bad) D3D(I,J,K) = Drat*DEN*RRSDONE2
		  if(VEL.ne.Bad) V3D(I,J,K) = Vrat*VEL

		  D3DT(I,J,K,LL) = Bad
		  V3DT(I,J,K,LL) = Bad
		  if(DEN.ne.Bad) D3DT(I,J,K,LL) = Drat*DEN*RRSDRR2
		  if(VEL.ne.Bad) V3DT(I,J,K,LL) = Vrat*VEL

		  if(Mo .eq. 1 .or. Mo .eq. -1) then
C  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.)/(K-1.)
				end if
			end if
		   end if
		  if(K .eq. 1) then
		      DD(I,J) = D3D(I,J,K)
		      VV(I,J) = V3D(I,J,K)
		  end if
		end do
	      end do
	    end do
	    
C	    print *, 'mo, TT, XCbeg+1, Xcbeg+2', mo, TT, XCbeg+1, Xcbeg+2

	    if(mo .ge. 0) call T3D_fill_time(t3d,nT,TT,XCbeg,XCend,XCbegROI,XCendROI)

	    if(mo .ge. 0) call T3D_write_nv(cPre,t3d,VV,DD,V3D,D3D,0,cStrAdd,cFile)

	   end if

	  end do
	end do

      return
      end