SUBROUTINE interpol(xinp,yinp,ninp,xout,yout,nout)
      IMPLICIT NONE

c=======================================================================
c
c  objet:
c  interpolation lineaire
c
c  ATTENTION!:
c  les tableaux x doivent etre croissants
c
c  Arguments:
c  ----------
c
c  entree:
c  -------
c  ninp        INT  nombre de points des donnes
c  xinp(ninp)  REAL abbcisses des donnees (doivent etre croissantes)
c  yinp(ninp)  REAL valeurs des donnees
c
c  sorties:
c  --------
c  nout        INT  nombre de points interpolles
c  xout(nout)  REAL abbcisses (doivent etre croissantes)
c  yout(nout)  REAL valeurs interpollees
c
c=======================================================================

c-----------------------------------------------------------------------
c    0. Declarations :
c    -----------------

      INTEGER ninp,nout
      INTEGER iinp,iout,iout0,iout1
      REAL xinp(ninp),yinp(ninp)
      REAL xout(nout),yout(nout)

c-----------------------------------------------------------------------
c   traitement des valeurs telles que xout<xinp(1):
c   -----------------------------------------------

      DO 100 iout=1,nout
         IF(xout(iout).lt.xinp(1)) THEN
              yout(iout)=yinp(1)
c   si on veut extrapoler aux bords:
c    1        +(xout(iout)-xinp(1))*(yinp(1)-yinp(2))/(xinp(1)-xinp(2))
          ELSE
              iout0=iout
              GO TO 200
          ENDIF
100   CONTINUE

      RETURN

200   CONTINUE

c-----------------------------------------------------------------------
c   traitement du centre:
c   ---------------------

      iinp=1
      DO 300 iout=iout0,nout
          IF(xout(iout).ge.xinp(ninp)) THEN 
             iout1=iout
             GO TO 400
          ENDIF
320       CONTINUE
          IF(  (xout(iout).ge.xinp(iinp)  ) .AND.
     1         (xout(iout).lt.xinp(iinp+1))        ) THEN
               yout(iout)=yinp(iinp)+
     1         (xout(iout)-xinp(iinp))*
     2         (yinp(iinp+1)-yinp(iinp))/(xinp(iinp+1)-xinp(iinp))
          ELSE
               iinp=iinp+1
               GO TO 320
          ENDIF
300   CONTINUE

      RETURN

400   CONTINUE

c-----------------------------------------------------------------------
c   traitement des points tels que xout>xinp(ninp):
c   -----------------------------------------------

      DO 500 iout=iout1,nout
          yout(iout)=yinp(ninp)
c   si on veut extrapoler aux bords:
c    1    +(xout(iout)-xinp(ninp))*
c    2    (yinp(ninp-1)-yinp(ninp))/(xinp(ninp-1)-xinp(ninp))
500   CONTINUE

      RETURN
      END