C+
C NAME:
C	Translate_c.f This is a test translate program that checks the ability to provide the next day output for the interative ENLIL system
C PURPOSE:
C	Read the externally-provided volumetric elements provided by an external 3D-MHD program, and translate these into 
C	the volumetric elements used by the tomography program 
C CATEGORY:
C	I/O
C CALLING SEQUENCE:
C	run translate
C
C INPUTS:
C	3-D data file from 3D-MHD program in a specified coordinate system
C
C OUTPUTS:
C	3-D data file that can be read by the UCSD tomography program
C
C MODIFICATION HISTORY:
C	November-2015, Bernard Jackson (UCSD)
C-
	program translate_c

	parameter      (NLmaxG	=  30000,	! Max # G data points  (13000) (400,000)
     &			NGREST	=  1,		! G temporal resolution factor better than 1 day
     &			NGRESS	=  1,		! Spatial G resolution factor better than 20 degrees
     &			NF	=  2,		! LOS resolution factor (Should be 1 if Spatial resolution is >= 2)
     &			nTmaxG	=  55*NGREST)	! Max # G times	(50) (100) (150)

	parameter	(nLng	= 54*NGRESS+1,	! # Longitudes, resolution is 360/(nLng-1) 	(55) (109) (217) (325)
     &			 iLngE	= 18+1,  	! # Longitudes to use in input ENLIL map     	(19)  (37)  (73) (109)
     &			 iLngM	= 72+1,  	! # Longitudes to use in input MS-FLUKSS map    (73)  (37)  (73) (109)
     &			 nLat	=  9*NGRESS+1,	! # Latitudes , resolution is 180/(nLat-1)      (10)  (19)  (37)  (55)
     &			 iLat	=  9+1,  	! # Latitudes , resolution is 180/(nLat-1)      (10)  (19)  (37)  (55)
     &			 nMap	=  31,		! # Radial heights. Range covered is [0,(nMap-1)*dRR) (31)
     &			 nLngLat = nLng*nLat,
     &			 nLngLatnTG =nLng*nLat*nTMaxG)

	real		XCEA(nTmaxG),
     &			DD1E(iLngE,iLat,nMap,nTmaxG),
     &			TT1E(iLngE,iLat,nMap,nTmaxG),
     &			VV3E(iLngE,iLat,nMap,nTmaxG,3),
     &			BB3E(iLngE,iLat,nMap,nTmaxG,3),
     &			DD1M(iLngM,iLat,nMap,nTmaxG),
     &			TT1M(iLngM,iLat,nMap,nTmaxG),
     &			VV3M(iLngM,iLat,nMap,nTmaxG,3),
     &			BB3M(iLngM,iLat,nMap,nTmaxG,3),
     &			DDD1(nLng,nLat,nMap,nTmaxG),
     &			TTT1(nLng,nLat,nMap,nTmaxG),
     &			VVV3(nLng,nLat,nMap,nTmaxG,3),
     &			BBB3(nLng,nLat,nMap,nTmaxG,3),
     &			RADS(nMap)

	real		FALLOFFD	/2.0/,
     &			FALLOFFV	/0.0/,
     &			FALLOFFBR	/2.0/,
     &			FALLOFFBT	/1.0/,
     &			FALLOFFBN	/1.34/

	real		XCMargin	/.5/

	parameter	(xInc	= 0.5)
	parameter	(nCar	= 2000)

	real*8		JDCar(nCar),
     &			JD, JEpoch, JEpoch1, JEpoch2,JDdate,
     &			dLngSun,dLatSun,dDisSun,
     &			DDoy,XC,
     &			MJDref,MJDrefB,MJDfile,
     &			MJDcntr

	real*8		XCintDG(nTmaxG)
	real*8		bad8
	real		NTBAD(nTmaxG)

	integer		MJDfrst		/ 46100/,
     &			MJDlast		/ 60000/

	character	cdate*8
	character	cfile*80
	character	ccfile*80
	character	cfileMJD*19
	character	cfileNUM*19
	character	ctimestp*21
	character	ctime*18
	character	cWild*80
	character	cheader*132
	character	cWild3DENLIL*80
	character	cWild3DMSFLUKSS*80
	character	cWild3DMSFLUKSSO*80
	character	cFmtF4*9		/'(55F10.4)'/
	character	cMon*3
	character	cStr*80

	logical		bForecast, bMJDtimefile, bMJDfile, btimefile
    	logical		bENLIL, bMSFLUKSS, bH3DMHD

	external	EARTH8
 	external 	EARTH

	include		'dirspec.h'
	include		'sun.h'

	bad8 = -1.0d0

	print *, ' '

	call AskWhat('Use an external input from: ENLIL, MS-FLUKSS, H3D-MHD?$1$1',I)
C	call AskWhat('Use an external input from: ENLIL, MS-FLUKSS, H3D-MHD?$1$2',I)

	    bENLIL    = I .eq. 1
	    bMSFLUKSS = I .eq. 2
	    bH3DMHD   = I .eq. 3

	 if(bENLIL)     MHDS = 1
	 if(bMSFLUKSS)  MHDS = 2
	 if(bH3DMHD)    MHDS = 3

C	Set up program
	   iOffUT = 7			! If PDT
	   iOffUT = 8			! If PST

	   NmidHRV = -9
	   call AskI4('Number of hours of tomography reconst. from UT midnight?',NmidHRV)

	   Doy1985 = 1.
	   iY1985 = 1985
	   call MAP_TZERO(EARTH,iY1985,Doy1985,.01,nCar,JDCar)
	   NCoff = N_CARRINGTON(iY1985,Doy1985)

	   call Julian(11,iYr1985,Doy1985,dble(MJDfrst),JEpoch)
	   XClo = NCoff+XMAP_SC_POS(EARTH,iYr1985,Doy1985,nCar,JDCar)
	   XClo = XClo-.5
C	   call Local2UT(iOffUT,iYr1985,Doy1985)
	   call Local2UT(NmidHRV,iYr1985,Doy1985)					! Changed 5/2/11 BVJ
	   XChi = max(NCoff+XMAP_SC_POS(EARTH,iYr1985,Doy1985,nCar,JDCar)-1.,XClo)
	   XChi = XChi + 1.0								! Changed 4/23/18 BVJ
C------
C iOffUT is an offset time (in hours) between the local computer time and UT. This is used only if the program
C is operated in forecast mode. Note that the local time is assumed to be the same as the
C computer time. The default offset is 7 hours (the difference between PDT and GST). If PST, the offset is 8 hours.
C
	   call AskYN('Forecast mode?$yes',bForecast)		! Regular or forecast mode
	   if (bForecast) then
	     call AskYN('Write MJD and timestamp values and stop?$no',bMJDtimefile)
	     call AskI4('Offset between local computer time and UT (hours)?',iOffUT)
	     if(bMJDtimefile) call AskYN('When write out/stop write actual MJD value?$no',bMJDfile)
	     NLVmin = NLVmin/2
	     NLGmin = NLGmin/2
	   else
	     call AskYN('Write out timestamp values and stop?$no',btimefile)
	   end if

C***************
	  if (bForecast) then			! Forecast mode

	    call Local2UT(iOffUT,iYr,Doy)	! calls the local computer time clock. With the offset calculates the Doy = universal time
	    call Julian(10,iYr,Doy,JD,JEpoch)	! Current time in MJD

	    write (*,'(/,10X,A,I6,A,I6,A,F11.3)') 'First MJD =',MJDfrst,
     &			'   Last MJD =',MJDlast,'   Current MJD =',JD

	    MJDref = JD					! MJD at center of plot (Earth position)

	    write (cStr,'(A,I5,A,F7.1,A)') '$',MJDfrst,'$',JD+1.5d0,'$0$' ! Allow the forecast to operate 1.5 day early (BVJ 2109/02/17) 
	    print *, ' '
	    call AskR8('Modified Julian Day (0 = STOP)'//cStr,MJDref)
	    if (MJDref .le. 0) call Say('IPSD2020','I','StopPlay','StopPlay')
C
C-------   Make sure iYr,Doy match MJDref
C	   Get heliographic coordinates for center of plot (Earth position)
C	   MJDref = MJDref + 7
C
	    call Julian(11,iYr,Doy,MJDref,JEpoch)
	    Idoy = int(Doy)
	    HR   = (Doy-float(Idoy))*24.0
	    IH   = int(HR)
	    IMN   = nint((HR-float(IH))*60.0)
	    call DATE_DOY(1,iYr,cMon,iM,iDay,iDoy)
	    write(*,'(A,I2.2,A,I2.2,A,I2.2,A,A,I5,A,F8.3,A)') 
     &		'      (This MJD is ',IH,':',IMN,' UT ',iDay,' ',cMon,iYr,' -- DOY',Doy,')' 

	    XCEarth = XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar)
	    call AdjustJDCar(XCEarth, nCar, JDCar, NCoff)

	    XCbeg   = XCEarth-.5

C            if(bMJDtimefile) then
	      IHfile   = 3
	      Idoyfile = Idoy
	      iYrfile  = iYr
	      iMfile   = iM
	      IDayfile = iDay
	      MJDfile  = MJDref
	      IMJD     = MJDref

C	      if(IH.lt.3) then
C		MJDfile = MJDref - 1.0d+0
C	        call Julian(11,iYrfile,Doyfile,MJDfile,JEpoch)
C		IDoyfile = int(Doyfile)
C	        call DATE_DOY(1,iYrfile,cMon,iMfile,IDayfile,IDoyfile)
C	      end if


	      if(.not.bMJDfile) then
	        if((MJDref-IMJD).gt.0.625) then		! add one day to the output file MJD if greater than a half day beyond 3 UT
		  MJDfile = MJDref + 1.0d+0
	          call Julian(11,iYrfile,Doyfile,MJDfile,JEpoch)
		  IDoyfile = int(Doyfile)
	          call DATE_DOY(1,iYrfile,cMon,iMfile,IDayfile,IDoyfile)
	        end if
	      end if

	    print*, ' '
	      write(*,'(A,F15.3)') 'MJDRef = ', MJDref
	      write(*,'(A,F15.3)') 'MJDfile = ', MJDfile

            if(bMJDtimefile) then	! Only do this in Forecast if ENLIL time, file start time, and CR are output in script, then stop.
C	      IDayfile = iDay
	      print*, '  '

	      write(*,'(A,F15.3)') 'MJDfile = ', MJDfile	! Modified Julian date of the run (MJDref or MJDfile).

	      write(*,'(A,I4,2I2.2,A,I2.2)') 'timestamp = ', iYrfile,iMfile,Idayfile,'T',IHfile

	      XCvalue = NCOFF + XCbeg				! XC value of the Carrington rotations output (CRNumber)
	      write(*,'(A,F15.3)') 'XCvalue = ', XCvalue

	      cfileMJD = 'MJD_value_         '
	      write(cfileMJD(11:19),'(F9.3)') MJDfile
	      cfileNUM = 'CRNumber__         '
	      write(cfileNUM(11:18),'(F8.3)') XCvalue
	      ctimestp = 'timestamp_            '
	      write(ctimestp(11:21),'(I4,2I2.2,A,I2.2)') iYrfile,iMfile,Idayfile,'T',IHfile

	      print*, cfileNUM,'  ', cfileMJD,'  ', ctimestp
C	 stop
	      close(13)
	      open (13, file=cFileNum,status='new',recl=550,access='sequential',form='formatted',iostat=iWrtNUM)
		write (13,'(A)',iostat=iwrtNUM) cFileNum								! write cfileNUM
	      close(13)
	      open (13, file=cFileMJD,status='new',recl=550,access='sequential',form='formatted',iostat=iWrtMJD)
		write (13,'(A)',iostat=iwrtNUM) cFileMJD								! write cfileMJD
	      close(13)
	      open (13, file=ctimestp,status='new',recl=550,access='sequential',form='formatted',iostat=iWrttim)
		write (13,'(A)',iostat=iwrtNUM) ctimestp								! write ctimestp
	      close(13)
	      stop

	    end if

	    XCtest0 = XCbeg-1.0			! Beginning of three rotations

	    XCtest1 = XCbeg-XCmargin		! Beginning of needed input files

	    XCtest2 = XCbeg+NC+1.0		! Ending of three rotations

C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

	  else						! Normal (non Forecast) mode where data are read from ENLIL inputs
	    write (cStr,'(A,F7.2,A,F7.2,A)') '$',XClo,'$',XChi,'$0$'
	    XCbeg = max(XClo,NCoff+XCbeg)  

	    if(bENLIL)    XCbeg = 2114.
	    if(bMSFLUKSS) XCbeg = 2114.

	    call AskR4('Carrington rotation (0 = STOP)?'//cStr,XCbeg)
	    if (XCbeg .le. 0) call Say('IPSD2020','I','StopPlay','StopPlay')

	    NintF = 1
C	    call AskI4('# Intermediate rotations interpolated',NintF)	

	    NC = 1
C	    call AskI4('# Rotations averaged',NC)
	    NC = min(NC,1+int(XChi-XCbeg) )	! Stay within available data range
	    NCC = NC

	    XCbeg = XCbeg-NCoff

	    call AdjustJDCar(XCbeg, nCar, JDCar, NCoff)

C-------	Info only: start time of Carrington rotation

	    call Julian(1,iYr,Doy,JDCar(int(XCbeg)),JEpoch)
	    iDoy = Doy
	    call DATE_DOY(1,iYr,cMon,iMon,iDay,iDoy)

	    write (*,'(/,8X,A,F7.2,A,I2,3A,I4,A,I3,A,/)') 'Carrington rotation ',NCoff+XCbeg,
     &		' starts on date ',iDay,'-',cMon,'-',iYr,' ( Doy ',iDoy,' )'

C-------	Get modified Carrington variable for center of plot
C		Get heliographic coordinates for center of plot (Earth position)

	    XCEarth = XCbeg+.5			! Center of plot
	    I = XCEarth				! Rotation containing plot center
	    MJDcntr = JDCar(I)+(XCEarth-I)*(JDCar(I+1)-JDCar(I))
						! Julian day for plot center
	    call Julian(1,iYr,Doy,MJDcntr,JEpoch) ! iYr,Doy for plot center
	    MJDcntr = MJDcntr-SUN__MJDtoJD	! Modified Julian day for plot center

	    XCtest0 = XCbeg-1.0			! Beginning of three rotations

	    XCtest1 = XCbeg-XCmargin		! Beginning of needed input files

	    XCtest2 = XCbeg+NC+1.0		! Ending of three rotations

C *******

	    if(btimefile) then			! Only do this in Normal mode if ENLIL time, and start time are output in script, then stop.

C	      print *, "Start of input files ",XCtest1

	      I = XCbeg
	      MJDrefB = JDCar(I)-SUN__MJDtoJD+(XCbeg-I)*(JDCar(I+1)-JDCar(I))

	      call Julian(11,iYr,Doy,MJDrefB,JEpoch)			! Where do the data files begin in year and day of year?
	      iDoy = int(Doy)
	      iH = nint((Doy - float(iDoY))*24.)
	      if(iH.lt.3.0) then 
		iDoy = iDoy - 1
		if(iDoy.lt.1) then
		  MJDrefB = MJDrefb - 1.0d0
	          call Julian(11,iYr,Doy,MJDrefB,JEpoch)		  
		  iDoy = int(Doy)
		end if
	      end if

C	      print *, 'XCbeg is ',iYr, Doy, iDoy 

	      call DATE_DOY(1,iYr,cMon,iM,iDay,iDoy)
C	      print *, iYr, iM, iDay     

	      ctimestp = 'time_year-iM-iDT03_00'
	      write(ctimestp(1:15),'(A,I4,A,I2.2,A,I2.2)') 'time_',iYr,'-',iM,'-',iday

C	      I = XCbeg - XCmargin
C	      MJDrefB = JDCar(I)-SUN__MJDtoJD+(XCbeg-XCmargin-I)*(JDCar(I+1)-JDCar(I))

	      MJDrefB = MJDrefB - 14.0d0
	      call Julian(11,iYr,Doy,MJDrefB,JEpoch)			! Where do the data files begin in year and day of year?
	      iDoy = int(Doy)
	      iH = nint((Doy - float(iDoY))*24.)
C	      if(iH.lt.3.0) then 
C		iDoy = iDoy - 1
		if(iDoy.lt.1) then
		  MJDrefB = MJDrefB - 1.0d0
	          call Julian(11,iYr,Doy,MJDrefB,JEpoch)		  
		  iDoy = int(Doy)
		end if
C	      end if
C	      print *, 'time above minus 14.0 days is ',iYr, Doy, iDoy

	      call DATE_DOY(1,iYr,cMon,iM,iDay,iDoy)

	      iH = 3
C	      print *, 'Start time is ', iYr, iM, IDay, iH  

	      ctime = 'strt_year-iM-iD-iH'
	      write(ctime(1:18),'(A,I4,A,I2.2,A,I2.2,A,I2.2)') 'strt_',iYr,'-',iM,'-',iday,'-',iH

	      print *, ctimestp, '   ',ctime

	      close(13)
	      open (13, file=ctimestp,status='new',recl=550,access='sequential',form='formatted',iostat=iWrttim)
		write (13,'(A)',iostat=iwrtNUM) ctimestp								! write ctimestp
	      close(13)
	      close(13)
	      open (13, file=ctime,status='new',recl=550,access='sequential',form='formatted',iostat=iWrttim)
		write (13,'(A)',iostat=iwrtNUM) ctime									! write ctime
	      close(13)
	      stop

	    end if
	  end if

	if(bforecast) then
	  print *, ' '
	  write(*,'(A,5F10.5)') 'Forecast Carrington variables', XCbeg, XCtest0, XCtest1, XCtest2, XCEarth
	else
	  write(*,'(A,4F10.5)') 'Non-forecast Carrington variables', XCbeg, XCtest0, XCtest1, XCtest2
	end if

	if(bENLIL) then
	  cWild3DENLIL = './ENLIL_inputs/enlil_ips_xxxxxxxxxx.dat'
	  print *, ' '
	  INN = 26
C	  Write(*,'(A,A50)'), 'ENLIL source ', cWild3DENLIL
c	  print *, ' '
	  iY = 2011
	  Mo = 8
	  iD = 13
	  iH = 3
	  dRR = 0.1

	  idaynum = NGREST

	  call ASKI4('Number of ENLIL values per day?',  idaynum)
	  print *, ' '

	  if(bForecast) then
	    Ifst = XCtest1 + 1 
            Ilst = Ifst -1
	    if(XCtest1.lt.XCMargin) then
              Ifst = XCtest1
	      Ilst = Ifst - 1
	    end if

C  Start from input start date
	    print *, 'MJDrefB = ', MJDfile, - (JDCar(Ifst)-JDCar(Ilst))

	    MJDrefB = MJDfile - (JDCar(Ifst)-JDCar(Ilst))

	    print *, 'MJDrefB = ', MJDrefB
						! Start MJD for data search
	  else
 
	    I = XCtest1
	    MJDrefB = JDCar(I)-SUN__MJDtoJD+(XCtest1-I)*(JDCar(I+1)-JDCar(I))
						! Start MJD for data search
	  end if

	  call Julian(11,iYr,Doy,MJDrefB,JEpoch)			! Where do the ENLIL data files begin in year and day of year?

	  print *, iYr, Doy

	  idoy = int(Doy)
	  call DATE_DOY(1,iYr,cMon,iM,iDay,iDoy)

	  iY = iYr
	  call AskI4('Beginning ENLIL test start year?', iY)
	  Mo = iM
 	  call AskI4('Beginning ENLIL test start month?',Mo)
	  iD = iDay
 	  call AskI4('Beginning ENLIL test start day?',  iD)
 	  call AskI4('Beginning test start hour?',       iH)

	  ifiles  = 48
	  if(bForecast) ifiles = 34
	  call ASKI4('Number of ENLIL files input?',     ifiles)



C Determine beginning start date from input year, month and start date. Suggestions above are given

	  call DATE_DOY(0,iY,cMon,Mo,iD,iDoy)		! determine the actual day of year (idoy) from above (BVJ Jan 03, 2019)
	  Doy = float(iDoy) + float(iH)/24.0

	  call Julian(10,iY,Doy,MJDrefB,JEpoch)		! Where do the ENLIL data files begin in MJD? (BVJ Jan 03, 2019)

	 print *, 'Here 0 - iY,cMon,Mo,iD',iY,cMon,Mo,iD

	  write(cWild3DENLIL(INN:INN+9),'(I4,3I2.2)') iY,mo,iD,iH
	  nT = ifiles
	  LLfst = 1
	  LLEnd = nT
	  write(*,'(A,A50)') 'ENLIL begining file ',cWild3DENLIL	   

	  cMon = '   '				! Needed because cMon is set above and if not blank provides this.
	  call DATE_DOY(0,iY,cMon,Mo,iD,iDoy)

	  do I=1,nT
	    DDOY       = dble(iDOY) + dble(iH + ((24.0d0/idaynum)*(I-1)))/24.0d0
	    XCintDG(I) = DDOY
	    write(*,'(A,I3,F15.8)') 'In Translate I, XCintDG(I)', I, XCintDG(I)
	  end do
	   
	  call enlil_read(cWild3DENLIL,INN,iLngE,iLat,nMap-1,nT,nTmaxG,LLfst,LLEnd,iY,XCintDG,MJDRefB, ! ENLIL has 1 too few distances
     &				RRS,dRR,RADS,DD1E,TT1E,VV3E,BB3E,iNUM,NTBAD)

C	  print *, ' '
	  do L=1,0
C	  do L=1,nT
	    do I=1,iLngE
	      do J=1,iLat
		do N=1,nMap-1
		  if(N.eq.10.and.J.eq.5.and.I.eq.12) write(*,'(A,4I3,F6.2,F7.1,5F7.3,F7.2)')  
     &		     'After ERead',I,J,N,L, DD1E(I,J,N,L),VV3E(I,J,N,L,1),VV3E(I,J,N,L,2),VV3E(I,J,N,L,3),
     &						BB3E(I,J,N,L,1),BB3E(I,J,N,L,2),BB3E(I,J,N,L,3),TT1E(I,J,N,L)
		end do
	      end do
	    end do
	  end do
C	stop

C	    call interpolate_enlil_m(iPr,iLngE,iLat,nMap-1,NNT,nTmaxG,nLng,nLat,nMap,iYrBG,XCintDG(I),nCar,JDCar,EarthLng,XCtest0,XCtest2, 
C     &	 		DD1E,TT1E,VV3E,BB3E,DDD1,TTT1,VVV3,BBB3)

	  RRS = RADS(1)
	  iYrBG = iY
	  EarthLngL = 0
	  do I=1,nT

	    NNT = I
	    call SunNewcomb(0,iYrBG,sngl(XCintDG(I)),dLngSun,dLatSun,dDisSun)

C****
C****	    New attempted fix
	    iskip = 1
	    if(iskip.eq.1) go to 999
           
	    call Julian(0,iYrBG,sngl(XCintDG(I)),JDdate,JEpoch1)
	    call Julian(0,2000.0,0.0,JDdate,JEpoch2)
	    sunlng19 = sngl(dLngSun)
	    sunlat19 = sngl(dLatSun)
	    call Precession_rot(JEpoch1,JEpoch2,sunlng19,sunlat19)
	    EarthLng = mod(sngl(sunLng19)+180.,360.)
	    EarthLat = -sunLat19
	    go to 998

999	    continue
C	    End New attempted fix
C****
C
C****
C    As in ipstd 3D
C	    call SunNewcomb(0,iYr,Doy,dLngSun,dLatSun,dDisSun)
C	    EarthLng = mod(sngl(dLngSun)+180.,360.)
C	    EarthLat = -dLatSun
C	    call ECLIPTIC_HELIOGRAPHIC(0,iYr,Doy,EarthLng,EarthLat)
C****
	    EarthLng = mod(sngl(dLngSun)+180.,360.)
	    EarthLat = sngl(-dLatSun)

998	    continue

	    iPr = 1
C	    iPr = 0

	    XXXX =  Earth(iYrBG,sngl(XCintDG(I)))
	    if(ipr.eq.1) write(*,'(A,I5,3F12.6)')'Earth Heliographic sub-Earth point ',iYrBG,XCintDG(I),XXXX

	    AriesLng = 0.0
	    AriesLat = 0.0
	    call ECLIPTIC_HELIOGRAPHIC(0,iYrBG,sngl(XCintDG(I)),AriesLng,AriesLat)
	    if(iPr.eq.1) write(*,'(A,I5,3F12.6)')'iYrBG, DOY, Eclip Aries Lon and Lat',iYrBG,XCintDG(I),AriesLng,AiresLat
	    if(iPr.eq.1) write(*,'(A,I5,3F12.6)')'iYrBG, DOY, Eclip Earth Lon and Lat',iYrBG,XCintDG(I),EarthLng,EarthLat
  
	    call ECLIPTIC_HELIOGRAPHIC(0,iYrBG,sngl(XCintDG(I)),EarthLng,EarthLat)
	    if(iPr.eq.1) write(*,'(A,I5,3F12.6)')'iYrBG, DOY, Helio Earth Lon and Lat',iYrBG,XCintDG(I),EarthLng,EarthLat

	    HGILng = EarthLng - AriesLng - 75.7570			! Why 75.7570 constant, I don't know
	    if(HGILng.lt.0.0) HGILng = HGILng + 360.0
	    if(iPr.eq.1) write(*,'(A,I5,3F12.6)')'iYrBG, DOY, HGI Earth Longitude    ',iYrBG,XCintDG(I),HGILng  
C****
	    EarthLng = HGILng
C
	    XCEA(I) = XMAP_SC_POS8(EARTH8,iYrBG,XCintDG(I),nCAR,JDCar)		! Earth location in Carrington coords. at DOY XCintDG
	    if(iPr.eq.1) write(*,'(A,I12,F12.5,I6)') 'I, XCEA(I), NCOFF         ',I, XCEA(I), NCOFF
	    if(iPr.eq.1) write(*,'(A,3F12.5)')       'XCtest0, XCtest2, EarthLng',XCtest0, XCtest2, EarthLng

C	    if(EarthLngL.ne.0.0) then
C	      COvsINTI = (EarthLng-EarthLngL)			! Amount of IHG ENLIL rotation per UCSD interval in deg.
C	      COvsINTE = (XCEA(I) - XCEAL)*1080./float(nLng-1)	! Amount of EHG UCSD rotation per UCSD interval in deg.
C	      COvsINT  = COvsINTE/COvsINTI			! Degrees amount of ENG relative to IHG in a UCSD interval at this time.
C	    else
C	      COvsINTI = (360.0*24.25/365.25)/float(idaynum)	! Average amount of IHG ENLIL rotation per UCSD interval in deg.
C	      COvsINTE = (360.0*27.232/365.25)/float(idaynum)	! Average amount of EHG UCSD rotation per UCSD interval in deg.
C	      COvsINT  = COvsINTE/COvsINTI			! Degrees amount of ENG relative to IHG in a UCSD interval at this time.
C	    end if

	    if(iPr.eq.1) print *, ' '
	    if(iPr.eq.1) write(*,'(I4,F12.4,5F12.5)') I, XCintDG(I), EarthLng, EarthLngL, COvsINT, COvsINTI, COvsINTE

	    iPr = 0
C	    iPr = 1
C	    iPr = 2

C	  NTVAL = 48
C
C	    call interpolate_enlil(iLngE,iLat,nMap-1,NNT,nTmaxG,nLng,nLat,nMap,EarthLng,COvsINT,XCEA(I),XCtest0,XCtest2, 
C     &	 		DD1E,TT1E,VV3E,BB3E,DDD1,TTT1,VVV3,BBB3)
C
C	    call interpolate_enlil_n(iPr,iLngE,iLat,nMap-1,NNT,nTmaxG,nLng,nLat,nMap,iYrBG,XCintDG(I),nCar,JDCar,EarthLng,XCtest0,XCtest2, 
C     &	 		DD1E,TT1E,VV3E,BB3E,DDD1,TTT1,VVV3,BBB3)
C
C	    call interpolate_enlil_m(iPr,iLngE,iLat,nMap-1,NNT,nTmaxG,nLng,nLat,nMap,iYrBG,XCintDG(I),nCar,JDCar,EarthLng,XCtest0,XCtest2, 
C     &	 		DD1E,TT1E,VV3E,BB3E,DDD1,TTT1,VVV3,BBB3)
C

C	    if(I.eq.1.or.I.eq.30.or.I.eq.47) then
C	      iPr = 1
C	      print *, ' '
C	      write(*,'(A,I3)') '***** Print ENLIL output in interpolate_enlil_nm **** at I =', I
C	      print *, ' '
C	    end if
C

	print *, 'Just before interpolate_enlil'
C 	if(I.eq.2) stop

	    call interpolate_enlil_nm(iPr,iLngE,iLat,nMap-1,NNT,nTmaxG,nLng,nLat,nMap,iYrBG,XCintDG(I),nCar,JDCar,EarthLng,
     &	 		XCtest0,XCtest2,DD1E,TT1E,VV3E,BB3E,DDD1,TTT1,VVV3,BBB3)
C
C	    print *, 'after interpolate_enlil ENLIL L, nt', L, nt

C	  if(NNT.eq.NTVAL) stop
	    iPr = 0
	    EarthLngL = EarthLng
	    XCEAL     = XCEA(I)
	  end do

C	stop

	end if

C******************************************

	if(bMSFLUKSS) then
	  cWild3DMSFLUKSS =  './MS-FLUKSS_inputs/msflukssxxxxxxxxxx'
	  cWild3DMSFLUKSSO = './MS-FLUKSS_nv3o_inputs/nv3oxxxxxxxxx_00001'
	  print *, ' '
	  INN = 28
	  INNO = 29
	  write(*,'(A,A50)'), 'MS-FLUKSS source ', cWild3DMSFLUKSS
	  print *, ' '
	  iY = 2011
	  Mo = 8
	  iD = 13
	  iH = 3
	  RRS = 0.1
	  dRR = 0.1

	  idaynum = NGREST

	  call ASKI4('Number of MS-FLUKSS values per day?',  idaynum)
	  print *, ' '

	  call Julian(11,iYr,Doy,MJDrefB,JEpoch)			! Where do the MS-FLUKSS data files begin in year and day of year?
	  idoy = int(Doy)
	  call DATE_DOY(1,iYr,cMon,iM,iDay,iDoy)

	  iY = iYr
	  call AskI4('Beginning MS-FLUKSS test start year?', iY)
	  Mo = iM
 	  call AskI4('Beginning MS-FLUKSS test start month?',Mo)
	  iD = iDay
	  if(iY.eq.2011.and.Mo.eq.8) iD = iD + 1			! This sets the MS-FLUKSS analysis right for CR 2114.0
 	  call AskI4('Beginning MS-FLUKSS test start day?',  iD)
 	  call AskI4('Beginning test start hour?',           iH)

	  ifiles  = 48
	  call ASKI4('Number of MS-FLUKSS files input?', ifiles)

	  write(cWild3DMSFLUKSS(INN:INN+9),'(I4,3I2.2)') iY,mo,iD,iH

	  nT = ifiles
	  LLfst = 1
	  LLEnd = nT
C	  write(*,'(A,A50)') 'MS-FLUKSS begining file ',cWild3DMSFLUKSS	   

	  cMon = '   '				! Needed because cMon is set above and if not blank provides this.
	  call DATE_DOY(0,iY,cMon,Mo,iD,iDoy)

	  iYrBG = iY
	  iDoyy = iDoy
	  do I=1,nT

	    DDOY       = dble(iDoyy) + dble(iH + ((24.0d0/idaynum)*(I-1)))/24.0d0
	    XCintDG(I) = DDOY

	    XCEA(I) = XMAP_SC_POS8(EARTH8,iYrBG,XCintDG(I),nCAR,JDCar)		! Earth location in Carrington coords. at DOY XCintDG 

	    write(*,'(A,I3,2F15.8)') 'I, DO XCintDG(I) XCEA(I)', I, XCintDG(I), XCEA(I)

	    if(I.eq.2) then
	      print *, ' '
	      write(cwild3DMSFLUKSS(INN:INN+8),'(F9.4)') XCEA(I) + NCOFF	   
	      if(I.eq.2) write(*,'(I4,A,A50)') I, ' First MS-FLUKSS file ',cWild3DMSFLUKSS	  
	      iDoy = int(XCintDG(I))
	      call DATE_DOY(1,iY,cMon,Mo,iD,iDoy)
	      write(*,'(A,5I6)') '     iYrBG, mo, iD, iH, iDoy:',iYrBG, mo, iD, iH, iDoy
	    end if
	  end do

C		!MS-FLUXSS inputs at UCSD Time locations
	  call msflukss_readt(cWild3DMSFLUKSS,INN,iLngM,iLat,nMap-1,nT,nTmaxG,LLfst,LLEnd,iY,XCintDG,  	! MS-FLUKSS has 1 too few distances
     &				RRS,dRR,RADS,DD1M,TT1M,VV3M,BB3M,iNUM)

	  iYrBG = iY
	  iMap  = nMap-1
	  nMapm = nMap-iMap
C	  print *, iYrBG, nMapM, nMap, iMap

C*******************
	print *, 'before interpolate nt', nt
	  do L=1,nT
	      do I=1,0
C	      do I=1,iLngM
		do J=1,iLat
		  do N=1,nMap-1
		    if(N.eq.10.and.J.eq.5) 
     &		    write(*,'(4I4,4F10.4)') L,I,J,N, DD1M(I,J,N,L),VV3M(I,J,N,L,1),BB3M(I,J,N,L,1),TT1M(I,J,N,L)
	          end do
		end do
	      end do		

	      NNT = L
	      call SunNewcomb(0,iYrBG,sngl(XCintDG(L)),dLngSun,dLatSun,dDisSun)

	      EarthLng = mod(sngl(dLngSun)+180.,360.)
	      EarthLat = sngl(-dLatSun)

	      AriesLng = 0.0
	      AriesLat = 0.0

	      call ECLIPTIC_HELIOGRAPHIC(0,iYrBG,sngl(XCintDG(L)),AriesLng,AriesLat)
C	      write(*,'(A,I5,3F12.6)')'iYrBG, DOY, Eclip Aries Lon and Lat',iYrBG,XCintDG(L),AriesLng,AiresLat
C	      write(*,'(A,I5,3F12.6)')'iYrBG, DOY, Eclip Earth Lon and Lat',iYrBG,XCintDG(L),EarthLng,EarthLat
  
	      call ECLIPTIC_HELIOGRAPHIC(0,iYrBG,sngl(XCintDG(L)),EarthLng,EarthLat)
C	      write(*,'(A,I5,3F12.6)')'iYrBG, DOY, Helio Earth Lon and Lat',iYrBG,XCintDG(L),EarthLng,EarthLat

	      HGILng = EarthLng - AriesLng - 75.7570			! Why 75.7570 constant, I don't know
	      if(HGILng.lt.0.0) HGILng = HGILng + 360.0
C	      write(*,'(A,I5,3F12.6)')'iYrBG, DOY, HGI Earth Longitude    ',iYrBG,XCintDG(L),HGILng  

	      EarthLng = HGILng

C	      XCEA(I) = XMAP_SC_POS8(EARTH8,iYrBG,XCintDG(L),nCAR,JDCar)	! Earth location in Carrington coords. at DOY XCintDG 
C	      write(*,'(A,I12,F12.5,I6)') 'L, XCEA(L), NCOFF         ',L, XCEA(L), NCOFF
C	      write(*,'(A,3F12.5)')       'XCtest0, XCtest2, EarthLng',XCtest0, XCtest2, EarthLng
C
	      iPr = 0
C	      iPr = 1
C
C	      call interpolate_enlil(iLngM,iLat,nMap-1,NNT,nTmaxG,nLng,nLat,nMap,EarthLng,XCEA(L),XCtest0,XCtest2, 
C     &	 		DD1M,TT1M,VV3M,BB3M,DDD1,TTT1,VVV3,BBB3)
C

C	      if(L.eq.1.or.L.eq.30.or.L.eq.47) then
C	        iPr = 1
C	        print *, ' '
C	        write(*,'(A,I3)') '***** Print MS-FLUKSS output in interpolate_enlil_nm **** at L =', L
C	        print *, ' '
C	      end if

	      call interpolate_enlil_nm(iPr,iLngM,iLat,nMap-1,NNT,nTmaxG,nLng,nLat,nMap,iYrBG,XCintDG(L),nCar,JDCar,EarthLng,
     &	 		XCtest0,XCtest2,DD1M,TT1M,VV3M,BB3M,DDD1,TTT1,VVV3,BBB3)
C
	      if(iPr.eq.1) print *, 'after interpolate_enlil_nm MS-FLUKSS L, nt', L, nt

	  end do

C*******************

	end if

C******************************************

C	  print *, ' '
	  do L=1,0
C	  do L=1,nT
	    do I=1,nLng
	      do J=1,nLat
		do N=1,nMap
		  if(N.eq.10.and.J.eq.5.and.I.eq.22) write(*,'(A,4I3,F6.2,F7.1,5F7.3,F7.2)')  
     &		     'Before ExWrite',I,J,N,L, DDD1(I,J,N,L),VVV3(I,J,N,L,1),VVV3(I,J,N,L,2),VVV3(I,J,N,L,3),
     &						BBB3(I,J,N,L,1),BBB3(I,J,N,L,2),BBB3(I,J,N,L,3),TTT1(I,J,N,L)
		end do
	      end do
	    end do
	  end do

C	 XC = XCEA(1) + float(NCOFF)
C	 write(*,'(A,2I6,F12.5)') 'Before External Write', iYrBG, NCOFF, XC
C	 stop
	 print *, ' '
	 if(iNUM.ge.1) then							! At least one file was read from the MHD inputs 
	   do I=1, nT
	     if(NTBAD(I).eq.99) XCintDG(I) = -XCintDG(I)			! If NTBAD from enlil_read don't write file and stop
	   end do	    	 
	   call ExternalWrite(MHDs,nLng,nLat,nMap,nT,nTmaxG,LLfst,LLEnd,JDCar,nCar,XCEA,NCOFF,iYrBG,XCintDG,
     &			RRS,dRR,FALLOFFD,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,DDD1,TTT1,VVV3,BBB3)
	 end if

	stop
	end