C+
C NAME:
C	Translate.f
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

	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		CR(1000), 
     &			dBXI(1000), 
     &			dBXT(1000), 
     &			dBYI(1000), 
     &			dBYT(1000), 
     &			dBZI(1000), 
     &			dBZT(1000)

	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,MJDref0,MJDrefB,MJDrefE,
     &			MJDcntr

	real*8		XCintDG(nTmaxG)
	real*8		bad8

	integer		MJDfrst		/ 46100/,
     &			MJDlast		/ 60000/

	character	cdate*8
	character	cfile*80
	character	ccfile*80
	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
    	logical		bENLIL, bMSFLUKSS, bH3DMHD

	external	EARTH8
 	external 	EARTH

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

	bad8 = -1.0d0

	print *, ' '

C	call AskWhat('Use an external input from: ENLIL, MS_FLUKSS, H3D-MHD?$1$1',I)
	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

	   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)

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).
C
	   call AskYN('Forecast mode?$no',bForecast)		! Regular or forecast mode
	   if (bForecast) then
	     call AskI4('Offset between local computer time and UT (hours)?',iOffUT)
	     NLVmin = NLVmin/2
	     NLGmin = NLGmin/2
	   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,'$0$'
	    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

	    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 mode
	   write (cStr,'(A,F7.2,A,F7.2,A)') '$',XClo,'$',XChi,'$0$'
	   XCbeg = max(XClo,NCoff+XCbeg)

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

C
C	************************
C	Begin Quick addition to read and write a new file with year and day of year from Carrington variable	
	   cfile ='./adcc.txt'
	   ccFile = './newadcc.txt'

	   open (13, file=cFile,status='old',recl=1000,access='sequential',form='formatted',iostat=iRead13)
	   open (14, file=ccFile,status='new',recl=1000,access='sequential',form='formatted')

	   if(iRead13.eq.0) write(*,'(A,A40)') ' The file opened was ',cfile
	   if(iRead13.ne.0) then
	     close(13)
	     write(*,'(A,A40)') ' The file not found was ',cfile
	     stop
	   end if
	   do I=1, 150
	     read (13,'(F4.0,6F11.5)',iostat=iRead13) CR(I), dBXI(I), dBXT(I), dBYI(I), dBYT(I), dBZI(I), dBZT(I)	! read values
	print *, 'I, CR(I)', I, CR(I)

	     XCbeg = CR(I)-NCoff

	     call AdjustJDCar(XCbeg, nCar, JDCar, NCoff)

	     call Julian(1,iYr,Doy1,JDCar(int(XCbeg)),JEpoch)
C	     iDoy = Doy
C	     call DATE_DOY(1,iYr,cMon,iMon,iDay,iDoy)
	     call Julian(1,iYr2,Doy2,JDCar(int(XCbeg+1.0)),JEpoch)
	     if(iYr2.gt.iYr) Doy2 = Doy2 + 365.0
	     DOYN = Doy1 + (Doy2 - Doy1)/2.
	     write(*,'(A,2I5,3F11.5)') 'iYr, iYr2, Doy1, Doy2, DOYN', iYr, iYr2, Doy1, Doy2, DOYN

	     write (14,'(F5.0,I6,F9.4,6F9.5)') 
     &		CR(I), iYr, DoYN, dBXI(I), dBXT(I), dBYI(I), dBYT(I), dBZI(I), dBZT(I)
	     write (*,'(F5.0,I6,F9.4,6F9.5)') 
     &		CR(I), iYr, DoYN, dBXI(I), dBXT(I), dBYI(I), dBYT(I), dBZI(I), dBZT(I)

	   end do
	   close(13)
	   close(14)
	stop

C	End Quick addition to read and wrte a new file with year and day of year from Carringto variable	
C	************************
C
	   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 (*,'(/,10X,A,I4,A,I2,3A,I4,A,I3,A,/)') 'Carrington rotation ',NCoff+int(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

	end if

	I = XCtest0
	MJDref0 = JDCar(I)-SUN__MJDtoJD+(XCtest0-I)*(JDCar(I+1)-JDCar(I))
						! Start of MJD for tomographic analysis
	I = XCtest1
	MJDrefB = JDCar(I)-SUN__MJDtoJD+(XCtest1-I)*(JDCar(I+1)-JDCar(I))
						! Start MJD for data search
	I = XCtest2
	MJDrefE = JDCar(I)-SUN__MJDtoJD+(XCtest2-I)*(JDCar(I+1)-JDCar(I))
						! End MJD for data search

	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 = 2007
	  Mo = 4
	  iD = 14
	  iH = 3
	  dRR = 0.1

	  idaynum = NGREST

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

	  call Julian(11,iYr,Doy,MJDrefB,JEpoch)			! Where do the ENLIL 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 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  = 1
	  call ASKI4('Number of ENLIL files input?',     ifiles)

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

	  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
C	    write(*,'(A,I3,F15.8)') 'I, XCintDG(I)', I, XCintDG(I)
	  end do

	  call enlil_read(cWild3DENLIL,INN,iLngE,iLat,nMap-1,nT,LLfst,LLEnd,iY,XCintDG,		! ENLIL has 1 too few distances
     &				RRS,dRR,RADS,DD1E,TT1E,VV3E,BB3E,iNUM)

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

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

C****
C****	    New attampted 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

	    XXXX =  Earth(iYrBG,sngl(XCintDG(I)))
	    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)
	    write(*,'(A,I5,3F12.6)')'iYrBG, DOY, Eclip Aries Lon and Lat',iYrBG,XCintDG(I),AriesLng,AiresLat
	    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)
	    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
	    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 
	    write(*,'(A,I12,F12.5,I6)') 'I, XCEA(I), NCOFF         ',I, XCEA(I), NCOFF
	    write(*,'(A,3F12.5)')       'XCtest0, XCtest2, EarthLng',XCtest0, XCtest2, EarthLng
	    call interpolate_enlil(iLngE,iLat,nMap-1,NNT,nT,nLng,nLat,nMap,EarthLng,XCEA(I),XCtest0,XCtest2, 
     &	 		DD1E,TT1E,VV3E,BB3E,DDD1,TTT1,VVV3,BBB3)
	  end do

	end if

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

	if(bMSFLUKSS) then
	  cWild3DMSFLUKSS =  './MS-FLUKSS_73X10_inputs/msflukssxxxxxxxxxx'
	  cWild3DMSFLUKSSO = './MS-FLUKSS_inputs/nv3oxxxxxxxxx_00001'
	  print *, ' '
	  INN = 34
	  INNO = 24
	  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 ENLIL 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
 	  call AskI4('Beginning MS-FLUKSS test start day?',  iD)
 	  call AskI4('Beginning test start hour?',           iH)

	  ifiles  = 45
	  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	   

	  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,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
	  print *, iYrBG, nMapM, nMap, iMap

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

	  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 = I
	      call SunNewcomb(0,iYrBG,sngl(XCintDG(I)),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(I)),AriesLng,AriesLat)
	      write(*,'(A,I5,3F12.6)')'iYrBG, DOY, Eclip Aries Lon and Lat',iYrBG,XCintDG(I),AriesLng,AiresLat
	      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)
	      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
	      write(*,'(A,I5,3F12.6)')'iYrBG, DOY, HGI Earth Longitude    ',iYrBG,XCintDG(I),HGILng  

	      EarthLng = HGILng

C	      XCEA(I) = XMAP_SC_POS8(EARTH8,iYrBG,XCintDG(I),nCAR,JDCar)	! Earth location in Carrington coords. at DOY XCintDG 
	      write(*,'(A,I12,F12.5,I6)') 'I, XCEA(I), NCOFF         ',I, XCEA(I), NCOFF
	      write(*,'(A,3F12.5)')       'XCtest0, XCtest2, EarthLng',XCtest0, XCtest2, EarthLng

	      call interpolate_enlil(iLngM,iLat,nMap-1,NNT,nT,nLng,nLat,nMap,EarthLng,XCEA(I),XCtest0,XCtest2, 
     &	 		DD1M,TT1M,VV3M,BB3M,DDD1,TTT1,VVV3,BBB3)

	      do I=1,nLng
		do J=1,nLat
		  do N=1,nMap
		    if(N.eq.1.and.J.eq.5.and.L.eq.21) write(*,'(4I4,4F10.4)')  
     &			I,J,N,L, DDD1(I,J,N,L),VVV3(I,J,N,L,1),BBB3(I,J,N,L,1),TTT1(I,J,N,L)
	          end do
		end do
	      end do

	  end do

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

	end if

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


C	 XC = XCEA(1) + float(NCOFF)
C	 write(*,'(A,2I6,F12.5)') 'Before External Write', iYrBG, NCOFF, XC
	
	 print *, ' '
	 if(iNUM.ge.1) then							! At least one file was read from the MHD inputs 
	    	 
	   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