C+
C NAME:
C	GIPSCAST
C PURPOSE:  
C	Calculates aftcast and forecast G values at two day intervals.
C CATEGORY: 
C	Linear interpolation
C CALLING SEQUENCE:
	subroutine GIPSCAST(MJDEarth,nX,Gval,NS)
C INPUTS:
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	NS			integer		# fore- and aftcast steps (2 day
C						steps)
C OUTPUTS: (to file)
C	DOYG(2,-NS:NS)		real		array with current,preceding,  
C		and 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
C		interpolated G values and the G value presently at Earth
C CALLS:
C	Say, iFilePath, bOpenFile, Julian, iFreeLun
C INCLUDE:
	include		'dirspec.h'
	include		'openfile.h'
	include		'sun.h'
C MODIFICATION HISTORY:
C	JUL-1992, Kari Winfield (UCSD)
C	AUG-1992, Susan Rappoport (UCSD)
C-
	parameter	(NSmax=4)

	real*8		MJDEarth
	real*8		JD
	real*8		JEpoch
	real		Gval(-nX/2:nX/2)

	character	APMJD *16		/'apmjd.dat'/
	character	DAYAVE*16		/'3-day_ave.dat'/
	character	SWMJD *16		/'sw_param_mjd.dat'/
	character	KPDST *16		/'kp_dst_ave.dat'/
	character	cStr*80

        integer		AP      (-NSmax:NSmax)
	integer		DEL_AP  (-NSmax:NSmax)
	integer		KP      (-NSmax:NSmax)
	integer		DEL_KP  (-NSmax:NSmax)
	integer		DST     (-NSmax:NSmax)
	integer		DEL_DST (-NSmax:NSmax)

	integer		NO_AP			/  -99/
	integer		NO_DEL_AP		/ -999/
	integer		NO_KP			/  -99/
	integer		NO_DEL_KP		/-9999/
	integer		NO_DST			/ -999/
	integer		NO_DEL_DST		/-9999/

	real		DOYG  (2,-NSmax:NSmax)
	real		DEL_G   (-NSmax:NSmax)
	real		DENS    (-NSmax:NSmax)
	real		DEL_DENS(-NSmax:NSmax)
	real		MAG     (-NSmax:NSmax)
	real		DEL_MAG (-NSmax:NSmax)
	real		VEL     (-NSmax:NSmax)
	real		DEL_VEL (-NSmax:NSmax)

	real		NO_G			/ -99.99/
	real		NO_DEL_G		/-999.99/
	real		NO_DENS			/-1.0   /
	real		NO_DEL_DENS		/-100.0 /
	real		NO_MAG			/ -999.9/
	real		NO_DEL_MAG		/ -999.9/
	real		NO_VEL			/-1.0   /
	real		NO_DEL_VEL		/-9999.0/

	integer		DayStep			/2/
	logical		bOpenFile

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

	iOPN = OPN__REOPEN+OPN__TEXT+OPN__READONLY+OPN__STOP
	call Julian(1,iYr0,Doy0,SUN__MJDtoJD+MJDEarth,JEpoch)

	do I=-NS,NS
	    iYr = iYr0
	    DOYG(1,I) = Doy0+DayStep*I
	    call Julian(0,iYr,DOYG(1,I),JD,JEpoch)

	    DOYG(2,I)	= NO_G
	    DEL_G(I)	= NO_DEL_G
	    AP(I)	= NO_AP
	    DEL_AP(I)	= NO_DEL_AP
	    VEL(I)	= NO_VEL
	    DEL_VEL(I)	= NO_DEL_VEL
	    KP(I)	= NO_KP
	    DEL_KP(I)	= NO_DEL_KP
	    DENS(I)	= NO_DENS
 	    DEL_DENS(I)	= NO_DEL_DENS
	    MAG(I)	= NO_MAG
	    DEL_MAG(I)	= NO_DEL_MAG
	    DST(I)	= NO_DST
	    DEL_DST(I)	= NO_DEL_DST

	end do

	! Interpolate horizontally to obtain g values corresponding to days of the 
	! year (by two day intervals)

	Step = (nX-1)/sngl(SUN__SYNODICP)*DayStep

	do I=-NS,NS
	    StepI = -Step*I
	    K1 = int(StepI)
	    if (StepI .lt. 0) K1 = K1-1
	    K2 = K1+1
	    if (I .eq. 0) K2 = K1

	    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

	! Read file for Julian day and Ap index for that day.

	I = iFilePath(cEnvi//'dat',0,' ',APMJD,cStr)
	if (bOpenFile(iOPN,iU,cStr,iRecl)) then
	    read (iU,'(A)') cStr

	    MJD1 = MJD0-DayStep*NS
	    read (iU,'(2I5)',iostat=I) MJD, IAP
	    do while (I .eq. 0 .and. MJD .lt. MJD1)
		read (iU,'(2I5)',iostat=I) MJD, IAP
	    end do

	    MJD1 = MJD0+DayStep*NS
	    do while (I .eq. 0 .and. MJD .le. MJD1)
		K = MJD-MJD0
		if (mod(K,DayStep) .eq. 0) AP(K/DayStep) = IAP
		read (iU,'(2I5)',iostat=I) MJD, IAP
	    end do
	    iU = iFreeLun(iU)

 	end if

	! Read file for magnetic field, solar wind density, and solar wind velocity
	! for that day.
 
	I = iFilePath(cEnvi//'dat',0,' ',DAYAVE,cStr)
	if (bOpenFile(iOPN,iU,cStr,iRecl)) then
	    read (iU,'(A)') cStr

	    cStr = '(1X,I5,I3,I5,2X,3F6.1,3I4)'

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

	    MJD1 = MJD0+DayStep*NS
	    do while (I .eq. 0 .and. MJD .le. MJD1)
		K = MJD-MJD0
		if (mod(K,DayStep) .eq. 0) then
		    K = K/DayStep
		    DENS(K) = XDENS
		    VEL(K) = XVEL
		    MAG(K) = XMAG
		end if
		read (iU,cStr,iostat=I) MJD,K,K,XMAG,XDENS,XVEL,K,K,K
	    end do
	    iU = iFreeLun(iU)

	end if

	I = iFilePath(cEnvi//'dat',0,' ',SWMJD,cStr)
	if (bOpenFile(iOPN,iU,cStr,iRecl)) then
	    read (iU,'(A)') cStr

	    cStr = '(2I5,F7.1,F6.0)'

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

	    MJD1 = MJD0+DayStep*NS
	    do while (I .eq. 0 .and. MJD .le. MJD1)
		K = MJD-MJD0
		if (mod(K,DayStep) .eq. 0) then
		    K = K/DayStep
		    DENS(K)= XDENS
		    VEL(K) = XVEL
		end if
		read (iU,cStr,iostat=I) MJD,K,XDENS,XVEL
	    end do
	    iU = iFreeLun(iU)

	end if
        
	! Read Kp and Dst indices.

	I = iFilePath(cEnvi//'dat',0,' ',KPDST,cStr)
	if (bOpenFile(iOPN,iU,cStr,iRecl)) then
	    read (iU,'(A)') cStr

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

	    MJD1 = MJD0-DayStep*NS
	    read (iU,cStr,iostat=I) MJD,K,K,XKP,XDST
	    do while (I .eq. 0 .and. MJD .lt. MJD1)
		read (iU,cStr,iostat=I) MJD,K,K,XKP,XDST
	    end do

	    MJD1 = MJD0+DayStep*NS
	    do while (I .eq. 0 .and. MJD .le. MJD1)
		K = MJD-MJD0
		if (mod(K,DayStep) .eq. 0) then
		    K = K/DayStep
		    KP(K) = nint(XKP)
		    DST(K) = nint(XDST)
		end if
		read (iU,cStr,iostat=I) MJD,K,K,XKP,XDST
	    end do
	    iU = iFreeLun(iU)

	end if

	! Calculate differences between values on aftcast and forecast days
	! and value on current day:
	!  AP = Ap index		DENS = density		VEL = velocity
	!  MAG = magnetic field		KP = Kp index		DST = Dst index
	!  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 (AP(I)	  .ne. NO_AP	.and. AP(0)	.ne. NO_AP  ) DEL_AP(I) = AP(I)-AP(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)
	    if (MAG(I)	  .ne. NO_MAG	.and. MAG(0)	.ne. NO_MAG ) DEL_MAG(I) = MAG(I)-MAG(0)
	    if (KP(I)	  .ne. NO_KP	.and. KP(0)	.ne. NO_KP  ) DEL_KP(I) = KP(I)-KP(0)
	    if (DST(I)	  .ne. NO_DST	.and. DST(0)	.ne. NO_DST ) DEL_DST(I) = DST(I)-DST(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,'(3X,3(6X,A4),2(4X,A4),6(6X,A4),4(4X,A4))') 
     &		' DOY','   G','  DG','  AP',' DAP','  Np',' DNp','  Vp',' DVp',
     &		'   B','  DB','  Kp',' DKp',' Dst','DDst'

	    write (iU,'(I3, 3F10.2, 2I8, 6F10.1, 4I8)')
     &		(DayStep*I, (DOYG(J,I),J=1,2), DEL_G(I),
     &		AP(I)  ,DEL_AP(I),
     &		DENS(I),DEL_DENS(I), VEL(I),DEL_VEL(I), MAG(I) ,DEL_MAG(I) ,
     &		KP(I)  ,DEL_KP(I)   , DST(I) ,DEL_DST(I) ,I=-NS,NS)

	    iU = iFreeLun(iU)
	end if

	return
	end