C INTERP    SOURCE    OF166741  25/02/20    21:16:46     12165          
      SUBROUTINE INTERP
C
C=======================================================================
C
C     Opérateur IPOL
C
C     SYNTAXE : voir notice
C
C=======================================================================
C
C    Remarques
C
C       Les listes LISTREE1 et LISTREE2 doivent se correspondre
C
C       L'évolution EVOL1 doit être élémentaire
C
C       ATTENTION : la liste des abscisses donnéee est supposée monotone
C
C=======================================================================
C
C    Auteur   :
C    Création :
C    Modifications : voir fiches d'anomalie
C
C=======================================================================

      IMPLICIT INTEGER(I-N)
      IMPLICIT    REAL*8(A-H,O-Z)
      PARAMETER (NCLE=2)
      CHARACTER*8 TYPE1,TYPE2,TYPEI,TYPO1,TYPO2,TYPOBJ
      CHARACTER*8 CHARIN,CHARRE
      LOGICAL     LOGIN,LOGRE
      CHARACTER*4 MCLE(NCLE)
      DATA MCLE/'TOUS','SPLI'/
      CHARACTER*4 CDER(2)
      DATA CDER/'DGAU','DDRO'/
      LOGICAL LCDER(2)
      REAL*8 XCDER(2)

-INC SMLREEL
      POINTEUR MLREE4.MLREEL,MLDERS.MLREEL,MLDER.MLREEL
-INC CCREEL
-INC SMCHPOI
-INC SMELEME
-INC SMCOORD

-INC PPARAM
-INC CCOPTIO
-INC SMCHAML
-INC SMTABLE
-INC SMEVOLL
      SEGMENT TABR
         REAL*8   TEMA(LOR)
      ENDSEGMENT

      KEVOLL = 0

*     syntaxe 3 ?
      CALL LIROBJ('NUAGE   ',INUA,0,IRETOU)
      IF (IERR.NE.0) RETURN
      IF(IRETOU.EQ.1) THEN
        CALL ACTOBJ('NUAGE   ',INUA,1)
        CALL IPLNU1(INUA)
        RETURN
      ENDIF
*     syntaxe 2 ?
      CALL LIROBJ('TABLE',IPOT,0,IRETOU)
      IF (IERR.NE.0) RETURN
      IF(IRETOU.EQ.1) THEN
        CALL  LIRREE(TEMPS,1,IRETOU)
        IF(IERR.NE.0) RETURN
        GOTO 50
      ENDIF

*     syntaxe 1 (INDIC=1 à 4)
      CALL LIROBJ('CHPOINT ',MTEMP,0,IRETOU)
      IF (IERR.NE.0) RETURN
      IF(IRETOU.EQ.1) THEN
        CALL ACTOBJ('CHPOINT ',MTEMP,1)
        IF (IERR.NE.0) RETURN
        ICHPO1 = MTEMP
C
C       Interpolation d'un point d'abscisse curviligne FLOT1
        CALL LIRREE(FLOT1,0,IRETOU)
        IF (IERR.NE.0) RETURN
        IF (IRETOU.EQ.1) THEN
          CALL IPLCUR(ICHPO1,FLOT1,IPOIN1)
          IF (IERR.NE.0) RETURN
          IF (IPOIN1.EQ.0) CALL ERREUR(251)
          CALL ECROBJ('POINT   ',IPOIN1)
          RETURN
        ELSE
          INDIC=2
        ENDIF
      ELSE
        CALL LIRREE(TEMPS,0,IRETOU)
        IF (IERR.NE.0) RETURN
        IF(IRETOU.EQ.1) THEN
          INDIC=1
        ELSE
          CALL LIROBJ('MCHAML ',IPO1,0,IRETOU)
          IF (IERR.NE.0) RETURN
          IF(IRETOU.EQ.1) THEN
            CALL ACTOBJ('MCHAML ',IPO1,1)
            IF (IERR.NE.0) RETURN
            INDIC=4
          ELSE
            CALL LIROBJ('LISTREEL',MLIST,0,IRETOU)
            IF (IERR.NE.0) RETURN
            IF(IRETOU.EQ.1) THEN
              CALL ACTOBJ('LISTREEL',MLIST,1)
              IF (IERR.NE.0) RETURN
            ELSE
*             Pas d opérande correcte trouvée
              CALL QUETYP(MOTERR(1:8),0,IRETOU)
              IF (IERR.NE.0) RETURN
              IF(IRETOU.NE.0) THEN
*               On ne veut pas d'objet de type %m1:8
                CALL ERREUR(39)
              ELSE
*               Cet opérateur a encore besoin d'un opérande
                CALL ERREUR(533)
              ENDIF
              RETURN
            ENDIF
            INDIC=3
          ENDIF
        ENDIF
      ENDIF

*     Lecture de la fonction à interpoler
      CALL LIROBJ('EVOLUTIO',IEVL,0,IRETOU)
      IF (IERR.NE.0) RETURN
      IF(IRETOU.EQ.1) THEN
        CALL ACTOBJ('EVOLUTIO',IEVL,1)
        IF (IERR.NE.0) RETURN
        MEVOLL=IEVL
*       Vérification que l'évolution est élémentaire
        IF (IEVOLL(/1).NE.1) THEN
*         Opération interdite sur un objet complexe
          CALL ERREUR(25)
          RETURN
        ENDIF
        KEVOLL=IEVOLL(1)
*       de type scalaire
        IF (ITYEVO.NE.'REEL') THEN
*         Opération interdite sur un objet complexe
          CALL ERREUR(25)
          RETURN
        ENDIF
*       peuplée de flottants
        IF ((TYPX.NE.'LISTREEL').OR.(TYPY.NE.'LISTREEL')) THEN
*         Certaines courbes ne sont pas du bon type
          CALL ERREUR(871)
          RETURN
        ENDIF
        KTE=IPROGX
        KFT=IPROGY
      ELSE
        CALL LIROBJ('LISTREEL',KTE,1,IRETOU)
        IF(IERR.NE.0) RETURN
        CALL LIROBJ('LISTREEL',KFT,1,IRETOU)
        IF(IERR.NE.0) RETURN
      ENDIF

*     Longueur des listes décrivant la fonction
      MLREE1=KTE
      MLREE2=KFT
      SEGACT MLREE1,MLREE2
      IF(MLREE1.PROG(/1).NE.MLREE2.PROG(/1)) THEN
*       Les suites n ont pas la même longueur
        CALL ERREUR(212)
        RETURN
      ENDIF
      LON=MLREE1.PROG(/1)

C     erreur 897 2 :
C     "La dimension des LISTREEL doit etre plus grande que 1"
      IF (LON.LT.1) THEN
         CALL ERREUR(897)
         RETURN
      ENDIF

c---- Lecture des options :
C     ITOUS : option TOUS
C     ISPLI : option SPLI
C     IHORS : comportement en dehors de l'intervalle de def. des donnees
C     IHORS = 0 : option ERREur
C     IHORS = 1 : option BORNE (valeur aux bornes)
C     IHORS = 2 : option EXTRapolation lineaire
      ITOUS = 0
      ISPLI = 0
      IHORS = 1
 2    CALL LIRMOT(MCLE,NCLE,ICLE,0)
      IF (IERR.NE.0) RETURN
      IF (ICLE.EQ.1) ITOUS = 1
      IF (ICLE.EQ.2) ISPLI = 1
C     IF (ICLE.EQ.3) IHORS = 1
C     IF (ICLE.EQ.4) IHORS = 2
C     IF (ICLE.EQ.5) IHORS = 0
      IF (ICLE.NE.0) GOTO 2

C     write(6,*) 'ITOUS =',ITOUS
C     write(6,*) 'ISPLI =',ISPLI
C     write(6,*) 'IHORS =',IHORS

C     Petite verification compatibilite options
      IF (ITOUS.EQ.1.AND.ISPLI.EQ.1) THEN
        CALL ERREUR(34)
        RETURN
      ENDIF
C     Option "TOUS" : pas d'extrapolation possible
C     IF (ITOUS.EQ.1.AND.IHORS.EQ.2) THEN
C       CALL ERREUR(34)
C       RETURN
C     ENDIF

C     Verif. donnees option 'TOUS'
      IF(ITOUS.EQ.1) THEN
        IF(INDIC.NE.1) THEN
c         Option %M1:8 incompatible avec les donnees
          MOTERR(1:8)='TOUS'
          CALL ERREUR(803)
c         On desire lire un nombre
          CALL ERREUR(15)
          RETURN
        ENDIF
        ISENS=0
c         GOTO 1
c       On va directement en 10 car SPLINE incompatible avec TOUS
        GOTO 10
      ENDIF

      TDEB=MLREE1.PROG(1)
      TFIN=MLREE1.PROG(LON)
*     Les x sont-ils croissants ou décroissants ?
      IF (TFIN.GE.TDEB) THEN
         ISENS=0
      ELSE
c        si decroissant, on retourne la liste
         ISENS=1
         JG=LON
         SEGINI,MLREE3,MLREE4
         DO ILON=1,LON
            MLREE3.PROG(ILON)=MLREE1.PROG(LON-ILON+1)
            MLREE4.PROG(ILON)=MLREE2.PROG(LON-ILON+1)
         ENDDO
         MLREE1=MLREE3
         MLREE2=MLREE4
         TDEB=MLREE1.PROG(1)
         TFIN=MLREE1.PROG(LON)
      ENDIF
C
C Vérification que la liste est ordonnée
C
       TPRE=TDEB
       DO ILON=2,LON
          TCOU=MLREE1.PROG(ILON)
          IF (TCOU.LT.TPRE) GOTO 6661
          TPRE=TCOU
       ENDDO
c       ENDIF
      GOTO 1
 6661 CONTINUE
C     erreur 249 2 : "La suite de reels doit etre croissante"
cbp : en realite elle doit etre monotone (decroissante possible)
      CALL ERREUR(249)
      RETURN

 1    CONTINUE
*
* Option SPLINE :
*
      IF (ISPLI.EQ.1) THEN
         LCDER(1)=.FALSE.
         LCDER(2)=.FALSE.
*  Lecture des mots clés et valeurs associées
*  On lit les valeurs des dérivées premières à gauche et à droite
*  Si elles ne sont pas données, c'est la condition à la limite
*  naturelle qui s'applique
 77      CONTINUE
         CALL LIRMOT(CDER,2,ICDER,0)
         IF (IERR.NE.0) RETURN
         IF (ICDER.GT.0) THEN
            LCDER(ICDER)=.TRUE.
            CALL LIRREE(XCDER(ICDER),1,IRETOU)
            IF(IERR.NE.0) RETURN
            GOTO 77
         ENDIF
         JG=LON
         SEGINI MLDERS
         SEGINI MLDER
         IF (LCDER(1)) THEN
* Cas où on prescrit la dérivée première à gauche
            MLDERS.PROG(1)=-0.5D0
            DX=MLREE1.PROG(2)-MLREE1.PROG(1)
            DY=MLREE2.PROG(2)-MLREE2.PROG(1)
            MLDER.PROG(1)=(3.D0/DX)*((DY/DX)-XCDER(1))
         ELSE
* Condition de bord naturelle (dérivée seconde nulle)
            MLDERS.PROG(1)=XZERO
            MLDER.PROG(1)=XZERO
         ENDIF
         DO ILON=2,LON-1
            XIM=MLREE1.PROG(ILON-1)
            XI=MLREE1.PROG(ILON)
            XIP=MLREE1.PROG(ILON+1)
            YIM=MLREE2.PROG(ILON-1)
            YI=MLREE2.PROG(ILON)
            YIP=MLREE2.PROG(ILON+1)
            DXIM=XI-XIM
            DXI2=XIP-XIM
            DXIP=XIP-XI
            XRAP=DXIM/DXI2
            XP=XRAP*MLDERS.PROG(ILON-1)+2.D0
            MLDERS.PROG(ILON)=(XRAP-1.D0)/XP
            DYIP=YIP-YI
            DYIM=YI-YIM
            MLDER.PROG(ILON)=(6.D0*(DYIP/DXIP-DYIM/DXIM)/DXI2-XRAP
     $           *MLDER.PROG(ILON-1))/XP
         ENDDO
         IF (LCDER(2)) THEN
            XQN=0.5D0
            DX=MLREE1.PROG(LON)-MLREE1.PROG(LON-1)
            DY=MLREE2.PROG(LON)-MLREE2.PROG(LON-1)
            XUN=(3.D0/DX)*(XCDER(2)-(DY/DX))
         ELSE
* Condition de bord naturelle (dérivée seconde nulle)
            XQN=0.D0
            XUN=0.D0
         ENDIF
         MLDERS.PROG(LON)=(XUN-XQN*MLDER.PROG(LON-1))/
     $        (XQN*MLDERS.PROG(LON-1)+1.D0)
         DO ILON=LON-1,1,-1
            MLDERS.PROG(ILON)=MLDERS.PROG(ILON)*MLDERS.PROG(ILON+1)
     $           +MLDER.PROG(ILON)
         ENDDO
         SEGSUP MLDER
c          write(*,*) 'MLDERS=',(MLDERS.PROG(iou),iou=1,LON)
      ELSE
         MLDERS=0
      ENDIF
*
*     Répartition suivant le type de l'objet fourni
*
C     write(6,*) 'INDIC =',INDIC
      GOTO (10,20,30,40) INDIC

******************   T0  FLOTTANT  *******************************

10    CONTINUE
      IF(ITOUS.EQ.1) THEN
        CALL INTER4(TEMPS,MLREE1,MLREE2,IHORS,IRET)
        IF (IERR.NE.0) RETURN
        IF (IRET.NE.0) THEN
          CALL ACTOBJ('LISTREEL',IRET,1)
          IF (IERR.NE.0) RETURN
          CALL ECROBJ('LISTREEL',IRET)
          IF (IERR.NE.0) RETURN
        ENDIF
      ELSE
        CALL INTER5(TEMPS,MLREE1,MLREE2,FT0,
     &              ISPLI,MLDERS,IHORS,IRET)
        IF (IERR.NE.0) RETURN
        IF (IRET.EQ.1) CALL ECRREE(FT0)
        IF (IERR.NE.0) RETURN
      ENDIF
      GOTO 999

*********************   T0  CHPOINT  *****************************

20    CONTINUE
      MCHPO1=MTEMP
      CALL NBCOMP(MCHPO1,'CHPOINT ',NB_Cmp)
      IF(NB_Cmp .NE. 1)THEN
        CALL ERREUR(180)
        RETURN
      ENDIF
      SEGINI,MCHPOI=MCHPO1
      MFT0=MCHPOI
      NS=IPCHP(/1)
      DO 21 IA=1,NS
        MSOUP1=IPCHP(IA)
        SEGINI,MSOUPO=MSOUP1
        NC=NOHARM(/1)
        IF(KEVOLL .NE. 0)THEN
          MSOUPO.NOCOMP(1)=KEVOLL.NOMEVY
        ELSE
          MSOUPO.NOCOMP(1)='SCAL'
        ENDIF
        IPCHP(IA)=MSOUPO
C        IPT1=IGEOC
C        SEGINI,IPT2=IPT1
C        IGEOC=IPT2
        MPOVA1=IPOVAL
        SEGINI,MPOVAL=MPOVA1
        IPOVAL=MPOVAL
        N=VPOCHA(/1)
        DO IB=1,N
          DO IC=1,NC
            TEMPS=VPOCHA(IB,IC)
            CALL INTER5(TEMPS,MLREE1,MLREE2,FT0,
     &                  ISPLI,MLDERS,IHORS,IRET)
            IF (IERR.NE.0) RETURN
            IF (IRET.EQ.0) THEN
              SEGSUP MCHPOI,MSOUPO,MPOVAL
              GOTO 999
            ENDIF
            VPOCHA(IB,IC)=FT0
          ENDDO
        ENDDO
21    CONTINUE
      CALL ACTOBJ('CHPOINT ',MFT0,1)
      IF (IERR.NE.0) RETURN
      CALL ECROBJ('CHPOINT ',MFT0)
      IF (IERR.NE.0) RETURN
      GOTO 999

*******************   T0 EST UN LISTREEL  ***********************

30    CONTINUE
      MLREE3=MLIST
      SEGACT MLREE3
      LONG=MLREE3.PROG(/1)
      JG=LONG
      SEGINI MLREEL
      MSOL=MLREEL
      DO 31 ILOOP=1,LONG
        TEMPS=MLREE3.PROG(ILOOP)
        CALL INTER5(TEMPS,MLREE1,MLREE2,FTO,
     &              ISPLI,MLDERS,IHORS,IRET)
        IF (IERR.NE.0) RETURN
        IF(IRET.EQ.0) GOTO 999
        PROG(ILOOP)=FTO
31    CONTINUE
      CALL ACTOBJ('LISTREEL',MSOL,1)
      IF (IERR.NE.0) RETURN
      CALL ECROBJ('LISTREEL',MSOL)
      IF (IERR.NE.0) RETURN
      GOTO 999

***********************   T0  MCHAML   ***************************

40    CONTINUE
      IRET=0
      MCHEL1=IPO1

      CALL NBCOMP(MCHEL1,'MCHAML  ',NB_Cmp)
      IF(NB_Cmp .NE. 1)THEN
        CALL ERREUR(320)
        RETURN
      ENDIF
      SEGINI,MCHELM=MCHEL1
      IRET=MCHELM
      NSOUS=IMACHE(/1)
      DO 72 IA=1,NSOUS
        MCHAM1=ICHAML(IA)
        SEGINI,MCHAML=MCHAM1
        IF(KEVOLL .NE. 0)THEN
          MCHAML.NOMCHE(1)=KEVOLL.NOMEVY
        ELSE
          MCHAML.NOMCHE(1)='SCAL'
        ENDIF
        ICHAML(IA)=MCHAML
        DO 75 ICOMP=1,IELVAL(/1)
          MELVA1 = IELVAL(ICOMP)
          SEGINI,MELVAL=MELVA1
          IELVAL(ICOMP) = MELVAL
          SEGACT MELVA1
          IF (TYPCHE(ICOMP).NE.'REAL*8') GOTO 75
          N1PTEL=VELCHE(/1)
          N1EL  =VELCHE(/2)
          N2PTEL=0
          N2EL  =0
          DO IB=1,N1PTEL
            DO ID=1,N1EL
              TEMPS=MELVA1.VELCHE(IB,ID)
              CALL INTER5(TEMPS,MLREE1,MLREE2,FT0,
     &                    ISPLI,MLDERS,IHORS,IREE)
              IF (IERR.NE.0) RETURN
              IF (IREE.EQ.0) THEN
                SEGSUP MCHELM,MCHAML,MELVAL
                GOTO 999
              ENDIF
              VELCHE(IB,ID)=FT0
            ENDDO
          ENDDO
75      CONTINUE
72    CONTINUE
      CALL ACTOBJ('MCHAML  ',IRET,1)
      IF (IERR.NE.0) RETURN
      CALL ECROBJ('MCHAML  ',IRET)
      IF (IERR.NE.0) RETURN
      GOTO 999

************************  OBJET1  TABLE  ******************************

50    CONTINUE
      MTABLE = IPOT
      SEGACT MTABLE
      LO = MLOTAB

*--   Vérification du format de la table
      IF (LO.LE.2) THEN
*       La table n'a pas le format désiré
        CALL ERREUR(647)
        RETURN
      ENDIF

      LOR = LO
      SEGINI TABR

*--   Vérification du sous-type de la table
*     IMOT est son indice dans la table
      IOK = 0
      DO 55 I=1,LO
        TYPE1 = MTABTI(I)
        IF(TYPE1.EQ.'MOT     ') THEN
          CHARIN = 'SOUSTYPE'
          TYPOBJ = ' '
          CALL ACCTAB(MTABLE,TYPE1,IVALIN,XVALIN,CHARIN,LOGIN,
     $               IOBIN,TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
          IF (IERR.NE.0) RETURN
          segact mtable
          IF(CHARIN.EQ.'SOUSTYPE') THEN
            IF(CHARRE.EQ.'RESULTAT') THEN
              IOK  = 1
              IMOT = I
            ENDIF
          ENDIF
        ENDIF
55    CONTINUE

      IF(IOK.EQ.0) THEN
*       Le sous-type de la table est incorrect
        CALL ERREUR(648)
        SEGSUP TABR
        RETURN
      ENDIF

*--   Vérification qu'on a bien des flottants en indice
*     en dehors du sous-type
      J = 0
      DO 56 I=1,LO
        IF (IMOT.EQ.I) GOTO 56
        J=J+1
        TYPEI = MTABTI(I)

        IF (TYPEI.NE.'FLOTTANT') THEN
*         La table n'a pas le format desire
          CALL ERREUR(647)
          SEGSUP TABR
          RETURN
        ENDIF

        TEMA(J) = RMTABI(I)
56    CONTINUE

*--   Vérification de l'ordonnancement des indices
      DO 57 I=1,LOR-2
        TEM1 = TEMA(I)
        TEM2 = TEMA(I+1)

        IF(TEM1.GT.TEM2) THEN
*         la liste des indices n'est pas ordonnee
          CALL ERREUR(211)
          SEGSUP TABR
          RETURN
        ENDIF

57    CONTINUE

      IF(IMOT.EQ.1) THEN
        TEM1   = RMTABI(2)
        TYPO1  = MTABTV(2)
        IVALO1 = MTABIV(2)
        NDEB   = 3
      ENDIF
      IF(IMOT.EQ.2) THEN
        TEM1   = RMTABI(1)
        TYPO1  = MTABTV(1)
        IVALO1 = MTABIV(1)
        NDEB   = 3
      ENDIF
      IF(IMOT.GE.3) THEN
        TEM1   = RMTABI(1)
        TYPO1  = MTABTV(1)
        IVALO1 = MTABIV(1)
        NDEB   = 2
      ENDIF

      DO 58 I=NDEB,LOR
        IF(IMOT.EQ.I) GOTO 58
        TEM2   = RMTABI(I)
        TYPO2  = MTABTV(I)
        IVALO2 = MTABIV(I)
        IF((TEM1.LE.TEMPS).AND.(TEMPS.LE.TEM2)) THEN
          DTEM = (TEMPS-TEM1)/(TEM2-TEM1)
          IF(TYPO1.EQ.'CHPOINT ') THEN
            IF(TYPO2.EQ.'CHPOINT ') THEN
              CALL ACTOBJ('CHPOINT ',IVALO2,1)
              IF (IERR.NE.0) RETURN
              CALL ACTOBJ('CHPOINT ',IVALO1,1)
              IF (IERR.NE.0) RETURN
              CALL ECROBJ('CHPOINT ',IVALO2)
              IF (IERR.NE.0) RETURN
              CALL ECROBJ('CHPOINT ',IVALO1)
              IF (IERR.NE.0) RETURN
            ELSE
              CALL ERREUR(647)
              SEGSUP TABR
              RETURN
            ENDIF
          ENDIF
          IF(TYPO1.EQ.'MCHAML  ') THEN
            IF(TYPO2.EQ.'MCHAML  ') THEN
              CALL ACTOBJ('MCHAML  ',IVALO2,1)
              IF (IERR.NE.0) RETURN
              CALL ACTOBJ('MCHAML  ',IVALO1,1)
              IF (IERR.NE.0) RETURN
              CALL ECROBJ('MCHAML  ',IVALO2)
              IF (IERR.NE.0) RETURN
              CALL ECROBJ('MCHAML  ',IVALO1)
              IF (IERR.NE.0) RETURN
            ELSE
              CALL ERREUR(647)
              SEGSUP TABR
              RETURN
            ENDIF
          ENDIF
          CALL ECRREE(DTEM)
          IF (IERR.NE.0) RETURN
          CALL ECRREE(1.D0-DTEM)
          IF (IERR.NE.0) RETURN
          CALL COLI
          IF (IERR.NE.0) RETURN
          SEGSUP TABR
          GOTO 500
        ENDIF
        TEM1   = TEM2
        TYPO1  = TYPO2
        IVALO1 = IVALO2
58    CONTINUE

*     Le temps est en dehors des limites de la table
      CALL ERREUR(210)
      SEGSUP TABR
      GOTO 500

999   CONTINUE
      IF (ISENS.GT.0) THEN
         SEGSUP MLREE1
         SEGSUP MLREE2
      ENDIF
      IF (ISPLI.EQ.1) THEN
         SEGSUP MLDERS
      ENDIF
*     Sortie
500   CONTINUE
      RETURN
      END