subroutine subintereuc(glat,glon,syr,smn,sdy,sut,res)
C.................................................................Sep. 2017
C Interpolate foF2 missed data in monthly/annual files From CBK-European foF2 maps
C
C.................................................................Jan. 2015
C Gulyaeva Tamara <gulyaeva@izmiran.ru>
C
C Subroutine to interpolate data from IONEX maps (foF2, hmF2, TEC)
C Input : 
C         glat - geographic latitute
C         glon - geographic longitude
C         syr - YEAR = '99' for 1999,'00' for 2000, etc.
C         smn - month
C         sdy - day of month
C         sut - hour, UT
C.................................................................
C Outpur : res - value of foF2 (0.1 MHz)
C 
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Use fc,hc,tc interpolated from directories dfc, dhc, dtc maps - update addresses for input file
C
	DIMENSION fhdat(51,36) !maps fcF2 glon(1:51)=[-10:2:90],glat(1:36)=[35:1:70]
	CHARACTER*2 syr,smn,sdy,sut
      CHARACTER*45 infile,txt
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Use fc,maps from directory /graf/Emap/
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	read(syr,*) ryr
	iyr=int(ryr)
	read(sdy,*) rdy1
	 ndy1=int(rdy1)
	read(smn,*) rmn1
	imn1=int(rmn1)
C Use fc,hc predictions from directories dfc or dhc
	infile='c:\web\graf\eum\YRMNDDUT.GRD' !path and name of input file !Tamara
	infile(17:18)=SYR  !Tamara
	infile(19:20)=SMN  !Tamara
	infile(21:22)=sdy  !Tamara
	infile(23:24)=sut  !Tamara
C204	format(46(3X,F8.5))
 205	format(A45)
C
	OPEN(10,FILE=INFILE,action='READ')		
C Extra lines:
      do k=1,5
	read(10,*,end=206) txt
	enddo
	goto 207
C Absent map:
  206 res=-1.0
      close(unit=10)
      RETURN
C   
  207	do ii=1,36
	do kk=1,51
	fhdat(kk,ii)=0.
	enddo
	enddo
		do ilat=1,36
	read(10,*,end=1204) (fhdat(k,ilat),k=1,51)	
	enddo
1204	close(unit=10)
C
C      interpolate for glat,glon
C Linear Interpolation   
	aglat=34.
	res=0.
	res11=0.
	res12=0.
	res21=0.
	res22=0.
	do 161 ilat=1,36! select pair of lats-lines (k) & (k+1)
c	aglat=35,36...,70.0
	aglat=34.0+1.0*ilat	   !
	aglat2=aglat+1.0
C
	if ((glat.ge.aglat).and.(glat.le.aglat2)) then  !   <<<<1
C Select pair of longs embracing glong
	 clon=-12.
	along=glon
C	if (along.ge.180.) along=along-360.0
	do j=1,51 ! select pair of longs (j) & (j+1)
	clon=clon+2.0
	clon2=clon+2.0
	if ((along.ge.clon).and.(along.le.clon2)) then  ! +++2
	res11=fhdat(j,ilat)
	res12=fhdat(j+1,ilat)
	res21=fhdat(j,ilat+1)
	res22=fhdat(j+1,ilat+1)
	res1=res11+(res12-res11)/2.*(along-clon)
	res2=res21+(res22-res21)/2.*(along-clon)
	res=res1+(res2-res1)*(glat-aglat)
	endif											 ! +++2
	enddo	! end select pair of longs (j) & (j+1)
	exit
	endif											 !   <<<<1
  161 continue		   ! interpolation end
	close(unit=10)
	RETURN
	END