C BA1D      SOURCE    CB215821  16/04/21    21:15:11     8920
      SUBROUTINE BA1D(SIG0,NSTRS,DEPST,VAR0,XMAT,NCOMAT,XCAR,TRAC,
     &                 NCOURB,SIGF,VARF,DEFP,KERRE)
C
C====&===1=========2=========3=========4=========5=========6=========7==
C Commentaires :   Subroutine permettant de mettre en oeuvre le
C                  modele BA1D (CEA/ENS Cachan) pour representer
C                  le comportement cyclique du beton sous
C                  sollicitations sismiques
C
C Traits       : -
C                -
C                -
C                -
C                -
C
C Auteur       : L. Moutoussamy - CEA, DEN, DANS, DM2S, SEMT, EMSI
C                B. Richard     - CEA, DEN, DANS, DM2S, SEMT, EMSI
C====&===1=========2=========3=========4=========5=========6=========7==
C
C----DECLARATION GENERALES----------------------------------------------
C
      IMPLICIT REAL*8(A-H,O-Z)

      DIMENSION SIG0(*),DEPST(*),VAR0(*),XMAT(*),XCAR(*),TRAC(*),
     &          SIGF(*),VARF(*) ,DEFP(*)
C
C----CARACTERISTIQUES MATERIAUX-----------------------------------------
C
      XYOUN = XMAT(1)
      XNU   = XMAT(2)

      XUEAL = XMAT(5)
      XFPLA = XMAT(6)
      XPFIS = XMAT(7)
      XHCIN = XMAT(8)
      XQFRA = XMAT(9)

      XAPIH = XMAT(10)
      XBPIH = XMAT(11)

C
C----VARIABLES INTERNES-------------------------------------------------
C
      XDP0 = VAR0(1)
      XZP0 = VAR0(2)
      XDM0 = VAR0(3)
      XZM0 = VAR0(4)
      XXPP = VAR0(5)
      XUPP = VAR0(6)
      XALP = VAR0(7)
      XXPM = VAR0(8)
      XUPM = VAR0(9)
      XALM = VAR0(10)
      XXPI = VAR0(11)
      XUPI = VAR0(12)
      XALI = VAR0(13)
      XU1  = VAR0(14)
      XU2  = VAR0(15)
      XU3  = VAR0(16)
C
C----AVANCE EN DEPLACEMENT----------------------------------------------
C
      XU1 = XU1 + DEPST(1)
      XU2 = XU2 + DEPST(2)
      XU3 = XU3 + DEPST(3)
C
C----CALCUL ENERGIE-----------------------------------------------------
C
      XRAI = XMAT(1)*XCAR(2)
      XY0  = 0.5D0*XRAI*(XUEAL)**2.0D0
      XUPL = ((XFPLA - (XRAI*XUEAL))/(XPFIS*XRAI)) + XUEAL
      XDY  = 1.0D0 - ((XUEAL + (XPFIS*(XUPL - XUEAL)))/XUPL)
      XGAM = (2.0D0/XAPIH)

      XENE = 0.5D0*XRAI*(XU3)**2.0D0

C      print*,'xene',xene
C      print*,'XY0',XY0

C
C----ENDOMMAGEMENT------------------------------------------------------
C
      IF (XU3.GE.0.0D0) THEN

         XTEST = XENE - (XY0 + XZP0)

         IF (XTEST.GE.0.0D0) THEN

            XDP0 = (1.0D0 - XPFIS) * (1.0D0 - (XY0/XENE)**XQFRA)
            XZP0 = XZP0 + 0.5D0 * XRAI *
     &             DEPST(3) * (2.0D0 * XU3 - DEPST(3))

         ENDIF

         XD0 = MIN(XDP0,XDY)

      ELSE

         XTEST = XENE - (XY0 + XZM0)

         IF (XTEST.GE.0.0D0) THEN

            XDM0 = (1.0D0 - XPFIS) * (1.0D0 - (XY0/XENE)**XQFRA)
            XZM0 = XZM0 + 0.5D0 * XRAI *
     &             DEPST(3) * (2.0D0 * XU3 - DEPST(3))

         ENDIF

         XD0 = MIN(XDM0,XDY)

      ENDIF

      IF (XU3.GE.0.0D0) THEN
         XD0 = MIN(XDP0,XDY)
      ELSE
         XD0 = MIN(XDM0,XDY)
      ENDIF

C      print*,'XD0',XD0

      XDH = MIN(MAX(XDM0,XDP0),XDY)
      XRAID = (1.0D0 - XD0) * XRAI

      GOTO 999

C
C----GLISSEMENT---------------------------------------------------------
C
      XFPTRI = XDH * XRAI * (XU3 - XUPI)
      XGPTRI = XFPTRI - XXPI
      XFPTRI0= ABS(XGPTRI)

      IF (XFPTRI0.GT.0.0D0) THEN

         IF (XGPTRI.GT.0.0D0) THEN
             XSIG = 1.0D0
         ELSE
             XSIG = -1.0D0
         ENDIF

         XFFPTRI = ABS(XGPTRI)
         XS1 = -1.0D0 * XSIG
         XS2 = (XGAM/2.0D0) * XXPI
         DLAMPI = XFFPTRI/(XDH * XRAI + XBPIH * (XS1 * (XS1 + XS2)))
         XUPI = XUPI + DLAMPI * XS1
         XXPI = XXPI - DLAMPI * XBPIH * (XS1 + XXPI * (XGAM/2.0D0))
         XFPTRI = XFPTRI + (XRAI * XDH * DLAMPI * XS1)
         XALI = XALI + DLAMPI
         XGPTRI = XFPTRI - XXPI

      ENDIF

      XFPI = XDH*XRAI*(XU3-XUPI)

 999  CONTINUE

C
C----PLASTICITE---------------------------------------------------------
C
      IF (XU3.GE.0.0D0) THEN

         FPTRIP = XRAID * (XU3 - XUPP)
         GPTRIP = FPTRIP - XXPP - (XFPLA/2.0D0)
         FFPTRIP= ABS(GPTRIP) - (XFPLA/2.0D0)

         IF (FFPTRIP.GT.0.0D0) THEN

            IF (GPTRIP.GT.0.0D0) THEN
               XSIG = 1.0D0
            ELSE
               XSIG = -1.0D0
            ENDIF

            DLAMPP = XSIG * XRAID *
     &               DEPST(3)/(XRAID + XHCIN)

            XDUPP  = DLAMPP * XSIG
            XUPP   = XUPP + XDUPP
            XXPP   = XXPP + (XHCIN * XDUPP)
            XALP   = XALP + DLAMPP

         ENDIF

         XUP = XUPP
         XUPM= 0.0D0

      ELSE

         FPTRIM = XRAID * (XU3 - XUPM)
         GPTRIM = FPTRIM - XXPM + (XFPLA/2.0D0)
         FFPTRIM= ABS(GPTRIM) - (XFPLA/2.0D0)

         IF (FFPTRIM.GT.0.0D0) THEN

            IF (GPTRIM.GT.0.0D0) THEN
               XSIG = 1.0D0
            ELSE
               XSIG = -1.0D0
            ENDIF

            DLAMPM = XSIG * XRAID *
     &               DEPST(3)/(XRAID + XHCIN)

            XDUPM  = DLAMPM * XSIG
            XUPM   = XUPM + XDUPM
            XXPM   = XXPM + (XHCIN * XDUPM)
            XALM   = XALM + DLAMPM

         ENDIF

         XUP = XUPM
         XUPP= 0.0D0

      ENDIF
C      print*,'xfpi=',xfpi
       XF = XRAID*(XU3 - XUP)
C      XF = XRAID*(XU3 - XUP) + XFPI
C      XF = XRAID*XU3
C
C----VARIABLES INTERNES-------------------------------------------------
C
      VARF(1) = XDP0
      VARF(2) = XZP0
      VARF(3) = XDM0
      VARF(4) = XZM0
      VARF(5) = XXPP
      VARF(6) = XUPP
      VARF(7) = XALP
      VARF(8) = XXPM
      VARF(9) = XUPM
      VARF(10)= XALM
      VARF(11)= XXPI
      VARF(12)= XUPI
      VARF(13)= XALI
      VARF(14)= XU1
      VARF(15)= XU2
      VARF(16)= XU3
C
C----CONTRAINTES--------------------------------------------------------
C
      CISA=0.5D0*XMAT(1)/(1.D0+XMAT(2))
      SIGF(1)=SIG0(1) + XMAT(1)*XCAR(4)*DEPST(1)
      SIGF(2)=SIG0(2) + CISA*(5.0D0/6.0D0)*XCAR(4)*DEPST(2)
C      SIGF(3)=SIG0(3) + XMAT(1)*XCAR(2)*DEPST(3)
      SIGF(3)=XF

      RETURN
      END