C HBMFNL    SOURCE    OF166741  26/05/11    21:15:11     12538          

*=======================================================================
*      Calcul de FNL(X) par AFT (Alternating Frequency Time)
*=======================================================================

      SUBROUTINE HBMFNL(NT,NHBM,NDDL,NFFT,KTQ,KTKAM,KTPHI,KTLIAA,
     &                  KTLIAB,KTPAS,KTFEX,KOCLFA,KOCLB1,mFNL)

*----- Declarations ----------------------------------------------------
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)

-INC SMLREEL
      POINTEUR mFNL.MLREEL

-INC TMDYNC

      SEGMENT mwork
        REAL*8 V1(NT)
        REAL*8 XT(NDDL,NFFT),FT(NDDL,NFFT),VT(NDDL,NFFT)
        REAL*8 XAUX(NDDL,4), VAUX(NDDL,4)
      ENDSEGMENT

      INTEGER NDDL,NT,NFFT
      REAL*8 PDT
      LOGICAL RIGIDE, GETJAC

      REAL*8 ZERO,ONE
      PARAMETER (ZERO=0.D0, ONE=1.D0)

*   Variables generalisees: coefficients de Fourier et frequence
      MTQ=KTQ
*   Matrices XK,KASM,XM,GAM,IGAM,DL
      MTKAM=KTKAM
*   Deformees modales
      MTPHI = KTPHI
      NSB   = XPHILB(/1)
      NPLSB = XPHILB(/2)
      NA2   = XPHILB(/3)
      IDIMB = XPHILB(/4)
*   Liaisons sur base A
      MTLIAA = KTLIAA
      NLIAA = IPALA(/1)
*   Liaisons sur base B
      MTLIAB = KTLIAB
      NLIAB  = IPALB(/1)
      NIP    = XABSCI(/2)
      NPLB   = JPLIB(/1)
*   Forces externes
      MTFEX=KTFEX
*   Truc local base A
      LOCLFA = KOCLFA
*   Truc local base B
      LOCLB1 = KOCLB1
*   Inconnues sur un pas de temps (AFT)
      MTPAS=KTPAS
      RIGIDE =.FALSE.
c   Initialisation (pas d'erreur)
      IERRD=0
      GETJAC=.FALSE.

      SEGINI,MWORK

*------ Time-domain displacements and velocities: XT, VT ---------------

      CALL IFT(Q1,NT,NDDL,NHBM,NFFT,GAM,XT)
      CALL HBMDVEC(NT,NHBM,NDDL,Q1,OMEG,V1)
      CALL IFT(V1,NT,NDDL,NHBM,NFFT,GAM,VT)

*------ Time-domain nonlinear forces, tangent matrix: FT -------

c   >>> Loop over time steps >>>
      PDT = ONE/NFFT
      DO I = 1,NFFT

c          initialisation des vecteurs du pas de temps
        DO L=1,NDDL
          XAUX(L,1) = XT(L,I)
*     VAUX pour la vitesse des liaisons A
          VAUX(L,1) = VT(L,I)
*     XAUX(:,2) permettait de retrouver le bonne vitesse pour les liaisons A
*     quand D2VLFA utilisait XVIT=
          XAUX(L,2) = XAUX(L,1)-PDT*VAUX(L,1)
          FTOTA(L,1)=ZERO
        ENDDO
        DO JJ=1,NDDL
          DO IJ=1,NDDL
            KTOTXA(IJ,JJ) = ZERO
            KTOTVA(IJ,JJ) = ZERO
          ENDDO
        ENDDO

c   Fnl(x,v) en base A (ddls modaux)
        IF (NLIAA.NE.0) THEN
          CALL D2VLFA(XAUX,VAUX,FTOTA,NDDL,IPALA,IPLIA,XPALA,XVALA,
     &                NLIAA,PDT,ZERO,I,1,FINERT,IVINIT,LOCLFA.FTESTA,
     &                KTOTXA,KTOTVA,GETJAC)
        ENDIF

c   Fnl(x,-) en base B (ddls physique)
        IF (NLIAB.NE.0) THEN
*         pour les liaisons *_frottement, on impose le glissement a vrai
          do il=1,NLIAB
            IPALB(il,2)=1
          enddo
          CALL D2VLFB(XAUX,VAUX,FTOTA,NDDL,IPALB,IPLIB,XPALB,XVALB,
     &                NLIAB,XPHILB,JPLIB,NPLB,IDIMB,FTOTB,FTOTBA,XPTB,
     &                PDT,ZERO,I,IBASB,IPLSB,INMSB,IORSB,NSB,NPLSB,NA2,
     &                1,FEXPSM,NPC1,IERRD,LOCLB1.FTESTB,
     &                XABSCI,XORDON,NIP,FEXB,RIGIDE,IAROTA,XCHPFB,
     &                KTOTXA,KTOTVA,KTOTXB,KTOTVB,GETJAC)
        ENDIF

c   Stockage de la force au pas de temps I
        DO J=1,NDDL
          FT(J,I) = FTOTA(J,1)
        ENDDO

cbp,2020-09 c          Decalage pour le pas de temps suivant (utile pour la vitesse tangentielle)
cbp,2020-09 *          remarque/TODO : avec HBM il faudrait utiliser Q2 directement pour la vitesse
cbp,2020-09              DO 22 ID = 1,IDIMB
cbp,2020-09              DO 20 IP = 1,NPLB
cbp,2020-09                 XPTB(IP,2,ID) = XPTB(IP,1,ID)
cbp,2020-09  20          CONTINUE
cbp,2020-09  22          CONTINUE

      ENDDO
c   <<< end of Loop over time steps <<<


c------ Fourier coefficients of nonlinear forces: FNL ------------------
      CALL DFT(FT,NDDL,NHBM,NFFT,IGAM,mFNL.prog(1))

*        DO IM = 1,NDDL
*          DO KI = 1,NFFT
*            XX51(KI) = FT(IM,KI)
*          ENDDO
*          CALL rfft1f(NFFT,1,XX51,NFFT,wwsave,lensav,work,lenwrk,ier)
*          IA = 0
*          DO KI = IM,NT-(NDDL-IM),NDDL
*           IA = IA+1
*           mFNL.PROG(KI) = XX51(IA)
*          ENDDO
*        ENDDO

c------ Menage ------------------
      SEGSUP,MWORK

      RETURN
      END