C DYPOL1    SOURCE    CB215821  25/04/22    21:15:03     12245          
      SUBROUTINE DYPOL1(Q1,Q2,NA1,IPLIA,XPALA,XVALA,NLIAA,IND,XPDT,
     &                  NPAS,NUMLIA,NMOD,FTOTA,IVINIT,IANNUL)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
*--------------------------------------------------------------------*
*                                                                    *
*     Operateur DYNE : algorithme de Fu - de Vogelaere               *
*     ________________________________________________               *
*                                                                    *
*     Calcul de la force elementaire pour la liaison POLYNOMIALE     *
*                                                                    *
*     Parametres:                                                    *
*                                                                    *
* e   Q1        Vecteur des deplacements generalises                 *
* e   Q2        Vecteur des vitesses generalisees                    *
* e   NA1       Nombre total d'inconnues en base A                   *
* e   IPLIA     Tableau des points de types support et origine       *
* es  XPALA     Tableau des parametres de la liaison                 *
* es  XVALA     Tableau des variables internes de la liaison         *
* e   NLIAA     Nombre de liaisons sur la base A                     *
* e   IND       Indice du demi-pas                                   *
* e   XPDT       pas de temps                             *
* e   NPAS      Numero du pas de temps                               *
* e   NUMLIA    Numero de la liaison consideree                      *
* e   NMOD      Nombre de modes "origine" de la liaison              *
* es  FTOTA     Forces exterieures totalisees sur la base A          *
* e   IVINIT    =1 si vitesses initiales, =0 sinon                   *
*                                                                    *
*     Auteur, date de creation:                                      *
*                                                                    *
*     Denis ROBERT, le 30 avril 1992.                                *
*                                                                    *
*--------------------------------------------------------------------*
*

-INC PPARAM
-INC CCOPTIO
*
      INTEGER IPLIA(NLIAA,*)
      REAL*8 XPALA(NLIAA,*),XVALA(NLIAA,4,*),Q1(NA1,*),FTOTA(NA1,*)
      REAL*8 Q2(NA1,*)
      PARAMETER (XZERO = 0.D0)
      PARAMETER (XUN = 1.D0 )
*
c       XPDT = XDT(NPAS)
      XPDTS2 = 0.5 * XPDT
*
*     Boucle sur les modes "origine"
*     Organisation de XPALA et indices :
*     COEF [ ED / RD / XJD / EV / RV / XJV / ... NMAX val de depl ... ]
* i=    1   IDEB                            IDEB+6
*
      IF (IIMPI.EQ.666) THEN
         WRITE(IOIMP,*)'       '
         WRITE(IOIMP,*)'NPAS   = ',NPAS
         WRITE(IOIMP,*)'IVINIT = ',IVINIT
         WRITE(IOIMP,*)'IND    = ',IND
      ENDIF
      XCOEF = XPALA(NUMLIA,1)
      IF (IANNUL.EQ.1) XCOEF = 0.D0
      FX = 1.D0
      IDEB = 2
      DO 10 I = 1,NMOD
         INA2  = IPLIA(NUMLIA,I+1)
         RD    = XPALA(NUMLIA,IDEB+1)
*        NRD : decalage pour trouver le deplacement lie au retard RD
         NRD   = INT(RD/XPDTS2) + 1
         RV    = XPALA(NUMLIA,IDEB+4)
*        NRV : decalage pour trouver la vitesse liee au retard RV
         NRV   = INT(RV/XPDTS2) + 1
*        NRV2 : un deplacement de plus pour calculer la vitesse
         NRV2  = NRV + 1
         NMAX  = MAX(NRD,NRV2)
         ED    = XPALA(NUMLIA,IDEB)
         EV    = XPALA(NUMLIA,IDEB+3)
         XJD   = XPALA(NUMLIA,IDEB+2)
         XJV   = XPALA(NUMLIA,IDEB+5)
*
*        Les deplacements sont reactualises en First In Last Out
*
         XDNEW = Q1(INA2,IND)
         J1 = IDEB + 6
         J2 = IDEB + NMAX + 5
         DO 15 J = J1,J2
            XAUX = XPALA(NUMLIA,J)
            XPALA(NUMLIA,J) = XDNEW
            XDNEW = XAUX
 15      CONTINUE
*
*        Boucle d'impression de XPALA
*
         IF (IIMPI.EQ.666) THEN
            WRITE(IOIMP,*)'DYPOL1 : IMOD = ',I
            DO 200 II = IDEB,J2
               WRITE(IOIMP,*)'XPALA(',NUMLIA,',',II,')=' 
     &                               ,XPALA(NUMLIA,II)
 200        CONTINUE
            WRITE(IOIMP,*)'NRD = ',NRD
            WRITE(IOIMP,*)'NRV = ',NRV
         ENDIF
*
*          faut-il calculer ( cad t > tau)
*
         jmax =  ideb + 5+Nmax
         Xtest = XPALA(NUMLIA,jmax)
         itest = int (xtest)
       IF ( itest .eq. 123456) then
*           on n'a pas ce qu'il faut pour calculer : on met zero
            fx = 0d0
       ELSE

         XDEPL = XPALA(NUMLIA,IDEB+5+NRD)
*
*            Calcul de la vitesse
*
         IF ((IND.EQ.3).AND.(IVINIT.EQ.1).and. (nrd.eq.1)) THEN
            XVIT = Q2(INA2,IND)
         ELSE
            XD1   = XPALA(NUMLIA,IDEB+5+NRV)
            XD2   = XPALA(NUMLIA,IDEB+5+NRV2)
               XVIT = (XD1 - XD2) / XPDTS2
         ENDIF
         IF (IIMPI.EQ.666) THEN
            WRITE(IOIMP,*)'---> XVIT = ',XVIT
         ENDIF
*
         IF ((XJD.NE.XZERO).AND.(ED.EQ.XZERO)) THEN
            ISEUIL = 1
         ELSE
            ISEUIL = 0
         ENDIF
         IF (XJD.GE.XZERO) THEN
*
*           partie positive de (XDEPL - XJD)
*
            IF (XDEPL.GE.XJD) THEN
               XD = ABS(XDEPL-XJD)
            ELSE
               XD = XZERO
            ENDIF
         ELSE
            IF (XDEPL.LT.XJD) THEN
               XD = ABS(XJD-XDEPL)
            ELSE
               XD = XZERO
            ENDIF
         ENDIF
         IF (IIMPI.EQ.666) THEN
            WRITE(IOIMP,*)'XD = ',XD
         ENDIF
         IF (XJV.GE.XZERO) THEN
            IF (XVIT.GE.XJV) THEN
               XV = ABS(XVIT-XJV)
            ELSE
               XV = XZERO
            ENDIF
         ELSE
            IF (XVIT.LT.XJV) THEN
               XV = ABS(XJV-XVIT)
            ELSE
               XV = XZERO
            ENDIF
         ENDIF
         IF (IIMPI.EQ.666) THEN
            WRITE(IOIMP,*)'XV = ',XV
         ENDIF
         IF (XD.GT.XZERO) THEN
            FXD = EXP(ED*LOG(XD))
         ELSE
            IF (ED.EQ.XZERO) THEN
               FXD = XUN
            ELSE
               FXD = XZERO
            ENDIF
         ENDIF
         IF (XV.GT.XZERO) THEN
            IF (ISEUIL.EQ.0) THEN
               FXV = EXP(EV*LOG(XV))
            ELSE
               IF (XD.GT.XZERO) THEN
               FXV = EXP(EV*LOG(XV))
               ELSE
                  FXV = XZERO
               ENDIF
            ENDIF
         ELSE
            IF (EV.EQ.XZERO) THEN
               FXV = XUN
            ELSE
               FXV = XZERO
            ENDIF
         ENDIF
         FX = FX * FXV * FXD
         IF (IIMPI.EQ.666) THEN
            WRITE(IOIMP,*)'FXD = ',FXD
            WRITE(IOIMP,*)'FXV = ',FXV
            WRITE(IOIMP,*)'FX  = ',FX
         ENDIF
         IDEB = IDEB + 6 + NMAX
      endif
 10   CONTINUE
      INA1 = IPLIA(NUMLIA,1)
      FXTOT = XCOEF * FX
      XVALA(NUMLIA,IND,1) = FXTOT
      FTOTA(INA1,IND) = FTOTA(INA1,IND) + FXTOT
*
      END