Laurène,

Your final goal is somewhat unclear but the GPS satellite position computation accuracy generally depends on what kind of ephemeris you use: broadcast or precise. Broadcast ephemerides are transmitted by the satellite itself in real time but their errors are about 1...2 m. Precise ephemerides are available only after post-processing as a SP3 file corresponding to a certain GPS day (see the link below). Precise ephemeris errors, in turn, depend on the post-processing time horizon. The most accurate version (so called "final" ephemeris) does have a centimeter-level error but is available only after a 2 week delay. Of course this post-processing includes interpolation procedures but also takes into account various disturbances acting upon the satellite. Interpolation alone is not enough.    

http://www.ngs.noaa.gov/orbits/orbit_data.shtml

i recommend you a MIT GAMIT/GLOBK software; and using: efemérides sp3 finals, earthquake file, etc. for another hand, do you use CGPS?; and why do you not use Ansi C?

Dear Laurène,

I know that you are specifically dealing with GPS ephemerides but, looking at how some other Fortran programs for interpolating planet ephemerides may give you a guideline to proceed. Please take a look at:

ftp://ssd.jpl.nasa.gov/pub/eph/planets/fortran/userguide.txt

In particular, the programs:

asc2eph.f

testeph.f