C+
C NAME:
C	GIPSIMP
C PURPOSE:  
C	Calculates aftcast and forecast G values at two day intervals.
C CATEGORY: 
C	Linear interpolation
C CALLING SEQUENCE:
	subroutine GIPSIMP(cVFile,MJDEarth,nX,Gval,NS)
C INPUTS:
C	cVFile		character*(*)	name of IMP data file
C	MJDEarth	real*8		modified julian day at the
C						position of the Earth
C	nX		integer		# latititudes
C	Gval(-nX/2:nX/2)real		G values at latitude of Earth
C					0 = no G-value is available
C	NS		integer		# fore- and aftcast steps (2 day steps)
C OUTPUTS: (to file)
C	DOYG(2,-NS:NS)	real		array with current, preceding, and
C		subsequent days of the year (at two day intervals) and the
C		corresponding G values for those days
C	DEL_G(-NS:NS)	real		differences between interpolated G values
C					and the G value presently at Earth
C CALLS:
C	bOpenFile, iFreeLun, Say, Julian
C INCLUDE:
	include		'openfile.h'
	include		'sun.h'
C MODIFICATION HISTORY:
C	July 1992, Kari Winfield (UCSD)
C	August 1992, Susan Rappoport (UCSD)
C-
	    character	cVFile*(*)
	    real*8	MJDEarth
	    integer	nX
	    real	Gval(-nX/2:nX/2)
	    integer	NS

	parameter	(NSmax=4)

	real*8		JD
	real*8		JEpoch

	character	cStr*80

	!-------
	! The number following the array declarations will be used to indicate that
	! no data values are available.

	real		DOYG    (2,-NSmax:NSmax),	NO_G		/ -99.99/,
     &			DEL_G   (-NSmax:NSmax),		NO_DEL_G	/-999.99/,
     &			DENS    (-NSmax:NSmax),		NO_DENS		/-1./,
     &			DEL_DENS(-NSmax:NSmax),		NO_DEL_DENS	/-100./,
     &			VEL     (-NSmax:NSmax),		NO_VEL		/-1./,
     &			DEL_VEL (-NSmax:NSmax),		NO_DEL_VEL	/-9999./

	integer		DayStep			/2/
	logical		bOpenFile


	if (NS .gt. NSmax) call Say('GIPSIMP','E','NSTOOSMALL',
     &		'increase value of parameter NS')

	!-------
	! Initialize DOYG(1,*) with times at fore- and aft- cast positions.
	! Initialize data arrays with values indicating 'NO DATA AVAILABLE'

	call Julian(1,iYr0,Doy0,SUN__MJDtoJD+MJDEarth,JEpoch)
					! iYr0,Doy0 refers to center of map
	do I=-NS,NS
	    iYr = iYr0
	    DOYG(1,I) = Doy0+DayStep*I	! Time at fore- and aft-cast positions
	    call Julian(0,iYr,DOYG(1,I),JD,JEpoch)

	    DOYG(2,I)	= NO_G
	    DEL_G(I)	= NO_DEL_G
	    VEL(I)	= NO_VEL
	    DEL_VEL(I)	= NO_DEL_VEL
	    DENS(I)	= NO_DENS
 	    DEL_DENS(I)	= NO_DEL_DENS

	end do

	!-------
	!  Interpolate horizontally to obtain g values corresponding to days of the 
	!  year (by intervals of DayStep days)
	!
	! (nX-1)/SUN__SYNODICP =Speed of Earth relative to solar surface in grid
	!			spacings per day
	! Step = # grid spacings traveled by Earth in DayStep days

	Step = (nX-1)/sngl(SUN__SYNODICP)*DayStep
	do I=-NS,NS		! Neg. I is aft-cast; Pos. I is forecast
	    StepI = -Step*I	! Displacement relative to map center (grid spacings)
				! For aft-casting, StepI is positive
	    K1 = int(StepI)	! K1,K2 are the enclosing grid points
	    if (StepI .lt. 0) K1 = K1-1
	    K2 = K1+1
	    if (I .eq. 0) K2 = K1
				! For I=0 use the grid value of map center
				! Linear interpolation between enclosing grid points
	    if (Gval(K1) .ne. 0 .and. Gval(K2) .ne. 0) then
		WeightK1 = 1.-(StepI-K1)
		DOYG(2,I) = WeightK1*Gval(K1)+(1.-WeightK1)*Gval(K2)
	    end if

	end do

	MJD0 = int(MJDEarth)	! Truncation of real*8 MJD at Earth
				! MJD0 refers to a day starting at midnight
	!-------
	! Read file for solar wind density and velocity.
 
	cStr = cVFile
	if (.not. bOpenFile(OPN__REOPEN+OPN__TEXT+OPN__READONLY+OPN__STOP,iU,cStr,iRecl)) I = I

	read (iU,'(A)') cStr

	cStr = '(1X,I5,2F7.1,I5,I4,I7)'

	!-------
	! Read records until the first requested modified Julian day is located.
	! The use of the DO WHILE loop assumes that the input file is
	! chronologically ordered.

	!-------
	! Read records until the earliest MJD time needed is located

	MJD1 = MJD0-DayStep*NS		! Earliest MJD time requested
	read (iU,cStr,iostat=I) MJD,XDENS,XVEL,K,K,K
	do while (I .eq. 0 .and. MJD .lt. MJD1)
	    read (iU,cStr,iostat=I) MJD,XDENS,XVEL,K,K,K
	end do

	!-------
	! Process records until the last MJD time needed is located. Extract only
	! data values at DayStep days intervals.
	! Note that invalid data in XDENS and XVEL should be indicated by the values
	! NO_DENS and NO_VEL. If this is not the case the calculation of DEL_DENS and
	! DEL_VEL doesn't work properly.

	MJD1 = MJD0+DayStep*NS		! Latest MJD time requested
	do while (I .eq. 0 .and. MJD .le. MJD1)
	    K = MJD-MJD0
	    if (mod(K,DayStep) .eq. 0) then	! Only DayStep offsetts are used
		K = K/DayStep
		DENS(K) = XDENS
		VEL(K) = XVEL
	    end if
	    read (iU,cStr,iostat=I) MJD,XDENS,XVEL,K,K,K
	end do

	iU = iFreeLun(iU)

	!-------
	! Calculate differences between values on aftcast and forecast days
	! and value on current day:
	!	DENS = density		VEL = velocity		DEL_G = IPS G value

	do I=-NS,NS
	    if (DOYG(2,I) .ne. NO_G    .and. DOYG(2,0) .ne. NO_G   ) DEL_G(I) = DOYG(2,I)-DOYG(2,0)
	    if (DENS(I)   .ne. NO_DENS .and. DENS(0)   .ne. NO_DENS) DEL_DENS(I)= DENS(I)-DENS(0)
	    if (VEL(I)	  .ne. NO_VEL  .and. VEL(0)    .ne. NO_VEL ) DEL_VEL(I) = VEL(I)-VEL(0)
	end do

	!-------
	! Write results to file

	write (cStr,'(I5,A)') MJD0,'.dat'
	iRecl = 0
	if (bOpenFile(OPN__TEXT+OPN__NEW+OPN__TRYINPUT+OPN__ONEPASS,iU,cStr,iRecl)) then
	    write (iU,'(8X,A,7X,A,10X,A,4(6X,A))')
     &		'DOY','G','DG','Np','DNp','Vp','DVp'

	    write (iU,'(1X,I2,F9.1,F10.2,F11.2,F7.1,F10.1,F8.1,F9.1)')
     &		(DayStep*I, (DOYG(J,I),J=1,2), DEL_G(I),
     &		DENS(I),DEL_DENS(I), VEL(I),DEL_VEL(I) ,I=-NS,NS)

	    iU = iFreeLun(iU)
	end if

	return
	end