C FVEH      SOURCE    AM        11/05/26    21:15:32     6982
      subroutine FVEH (da,force,q,XMKel,Gradgs,distmin_prec,hh,distmin)
C     Routine permettant la fonction vectorielle d'écrouissage

        IMPLICIT INTEGER(I-N)
        IMPLICIT REAL*8(A-H,O-Z)
        DIMENSION da(20), XMKel(5,5), force(5), q(5), Gradgs(5), hh(5)
        DIMENSION Gradf_star(4), force_star(4), tau(4)
*
        a=da(8)
        b=da(9)
        c=da(10)
        d=da(11)
        e=da(12)
        f=da(13)
        a6=da(16)

C       projection du vecteur Gradfs dans l'hyperplan (H*,M*)
        gamma=1.

      den1_star = a*force(1)**c*(gamma-force(1))**d
      den2_star = b*force(1)**e*(gamma-force(1))**f
      force_star(1)=force(2)/den1_star
        force_star(2)=force(3)/den2_star
        force_star(3)=force(4)/den1_star
        force_star(4)=force(5)/den2_star

      Gradf_star(1)=2/q(5)**2*(force_star(1)-q(1))
        Gradf_star(2)=2/q(5)**2*(force_star(2)-q(2))
        Gradf_star(3)=2/q(5)**2*(force_star(3)-q(3))
        Gradf_star(4)=2/q(5)**2*(force_star(4)-q(4))

        qinf=1-q(5)
       DO I=1,4
        tau(I)=q(I)
       ENDDO
        call SCALPR(4,tau,tau,Xnorm_tau)
        Xnorm_tau=sqrt(Xnorm_tau)

        call SCALPR(4,Gradf_star,Gradf_star,Xnorm_Gradf_star)
        Xnorm_Gradf_star=sqrt(Xnorm_Gradf_star)
        if (Xnorm_Gradf_star.ne.0) then

C Test nécessaire pour éviter d'avoir une divergence des variables d'écrouissage
        if (abs(q(1)).gt.abs(Gradf_star(1))/Xnorm_Gradf_star*qinf) then
        q(1)=Gradf_star(1)/Xnorm_Gradf_star*qinf
        endif
        if (abs(q(2)).gt.abs(Gradf_star(2))/Xnorm_Gradf_star*qinf) then
        q(2)=Gradf_star(2)/Xnorm_Gradf_star*qinf
        endif
        if (abs(q(3)).gt.abs(Gradf_star(3))/Xnorm_Gradf_star*qinf) then
        q(3)=Gradf_star(3)/Xnorm_Gradf_star*qinf
        endif
        if (abs(q(4)).gt.abs(Gradf_star(4))/Xnorm_Gradf_star*qinf) then
        q(4)=Gradf_star(4)/Xnorm_Gradf_star*qinf
        endif

        dist=(sqrt((Gradf_star(1)/Xnorm_Gradf_star*qinf-q(1))**2+
     & (Gradf_star(2)/Xnorm_Gradf_star*qinf-q(2))**2+
     & (Gradf_star(3)/Xnorm_Gradf_star*qinf-q(3))**2+
     & (Gradf_star(4)/Xnorm_Gradf_star*qinf-q(4))**2))

        else
        dist=0.9999999
        endif

        if (dist.le.distmin_prec) then
        distmin=dist
        else
        distmin=distmin_prec
        endif


        dist0=2/(1+a6)
        if (((dist+distmin*(1-dist0))/dist0).gt.1) then
        hh(1)=-Gradgs(2)/((1+a6)*den1_star)*1E12*XMKel(2,2)
        hh(2)=-Gradgs(3)/((1+a6)*den2_star)*1E12*XMKel(3,3)
        hh(3)=-Gradgs(4)/((1+a6)*den1_star)*1E12*XMKel(4,4)
        hh(4)=-Gradgs(5)/((1+a6)*den2_star)*1E12*XMKel(5,5)
        hh(5)=-a6*abs((q(1)*Gradgs(2)/((1+a6)*den1_star)*XMKel(2,2)
     &+q(2)*Gradgs(3)/((1+a6)*den2_star)*XMKel(3,3)
     &+q(3)*Gradgs(4)/((1+a6)*den1_star)*XMKel(4,4)
     &+q(4)*Gradgs(5)/((1+a6)*den2_star)*XMKel(5,5))*1E12)/Xnorm_tau
        else
        hh(1)=-Gradgs(2)/((1+a6)*den1_star)*XMKel(2,2)
     & *(dist)/(1-dist/dist0+distmin/dist0*(1-dist0))
        hh(2)=-Gradgs(3)/((1+a6)*den2_star)*XMKel(3,3)
     & *(dist)/(1-dist/dist0+distmin/dist0*(1-dist0))
        hh(3)=-Gradgs(4)/((1+a6)*den1_star)*XMKel(4,4)
     & *(dist)/(1-dist/dist0+distmin/dist0*(1-dist0))
        hh(4)=-Gradgs(5)/((1+a6)*den2_star)*XMKel(5,5)
     & *(dist)/(1-dist/dist0+distmin/dist0*(1-dist0))
        hh(5)=-a6*abs((q(1)*Gradgs(2)/((1+a6)*den1_star)*XMKel(2,2)
     &+q(2)*Gradgs(3)/((1+a6)*den2_star)*XMKel(3,3)
     &+q(3)*Gradgs(4)/((1+a6)*den1_star)*XMKel(4,4)
     &+q(4)*Gradgs(5)/((1+a6)*den2_star)*XMKel(5,5))
     & *(dist)/(1-dist/dist0+distmin/dist0*(1-dist0)))/Xnorm_tau

        endif


        return
        end