asso_plas76.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    asso_plas76             ../engine/source/materials/mat/mat076/asso_plas76.F
!||--- called by ------------------------------------------------------
!||    sigeps76                ../engine/source/materials/mat/mat076/sigeps76.F
!||--- calls      -----------------------------------------------------
!||    table_mat_vinterp       ../engine/source/materials/tools/table_mat_vinterp.F
!||    vinter                  ../engine/source/tools/curve/vinter.F
!||--- uses       -----------------------------------------------------
!||    table4d_mod             ../common_source/modules/table4d_mod.F
!||    table_mat_vinterp_mod   ../engine/source/materials/tools/table_mat_vinterp.F
!||====================================================================
      SUBROUTINE asso_plas76(
     1     NEL     ,NUPARAM,NUVAR   ,NFUNC   ,IFUNC   ,NVARTMP ,
     2     NPF     ,TF     ,TIME    ,TIMESTEP,UPARAM  ,VARTMP  ,
     3     RHO0    ,PLA    ,DPLA    ,ET      ,NUMTABL ,TABLE   ,
     3     DEPSXX  ,DEPSYY ,DEPSZZ  ,DEPSXY  ,DEPSYZ  ,DEPSZX  ,
     4     SIGOXX  ,SIGOYY ,SIGOZZ  ,SIGOXY  ,SIGOYZ  ,SIGOZX  ,
     5     SIGNXX  ,SIGNYY ,SIGNZZ  ,SIGNXY  ,SIGNYZ  ,SIGNZX  ,
     6     SOUNDSP ,UVAR   ,OFF     ,EPSD    ,YLD     )
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE table4d_mod
      USE table_mat_vinterp_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "tabsiz_c.inc"
C=======================================================================
C   I N P U T   A r g u m e n t s 
C-----------------------------------------------
      INTEGER NEL,NUPARAM,NUVAR,NFUNC,NUMTABL
      INTEGER ,INTENT(IN)    :: NVARTMP
      my_real,INTENT(IN) :: TIME,TIMESTEP
      my_real,DIMENSION(NUPARAM),INTENT(IN) :: UPARAM(NUPARAM)
      my_real,DIMENSION(NEL),INTENT(IN) :: 
     .   RHO0,DEPSXX,DEPSYY,DEPSZZ,
     .   DEPSXY,DEPSYZ,DEPSZX,SIGOXX,SIGOYY,SIGOZZ,
     .   sigoxy,sigoyz,sigozx 
C-----------------------------------------------
C   O U T P U T   A r g u m e n t s
C-----------------------------------------------
      INTEGER ,DIMENSION(NEL,NVARTMP) ,INTENT(INOUT) :: VARTMP
      my_real,DIMENSION(NEL),INTENT(OUT) ::
     .   SIGNXX,SIGNYY,SIGNZZ,SIGNXY,SIGNYZ,SIGNZX,
     .   SOUNDSP,DPLA,ET
C-----------------------------------------------
C   I N P U T   O U T P U T   A r g u m e n t s 
C-----------------------------------------------
      my_real,DIMENSION(NEL),INTENT(INOUT) :: 
     .   off,yld,pla,epsd
      my_real,DIMENSION(NEL,NUVAR),INTENT(INOUT) :: 
     .   uvar
      TYPE(table_4d_),DIMENSION(NUMTABL),INTENT(IN) :: TABLE
C-----------------------------------------------
C   VARIABLES FOR FUNCTION INTERPOLATION 
C-----------------------------------------------
      INTEGER :: NPF(SNPC), IFUNC(NFUNC)
      my_real :: tf(stf)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: I, II, ICONV, NINDX, ICAS, ITER
      INTEGER ,PARAMETER :: FUNC_TRAC  = 1
      INTEGER ,PARAMETER :: FUNC_COMP  = 2
      INTEGER ,PARAMETER :: FUNC_SHEAR = 3
      integer ,DIMENSION(NEL)   :: indx,iad,ilen
      my_real :: g,aa1,aa2,c1,nupc,xfac
      my_real ,DIMENSION(NEL) :: plat,plac,plas,epspt,epspc,epsps,sigt,
     .    sigs,dsigs_dp,sigc,dsigc_dp,sdxx,sdyy,sdzz,sdxy,sdyz,sdzx,dsigt_dp,
     .    nup,p,svm,a0,a1,a2,phi,dplac,dplat,dplas,df
      my_real
     .    normxx,normyy,normzz,normxy,normyz,normzx,norm,aa,cb,
     .    normxx_n,normyy_n,normzz_n,normxy_n,normyz_n,normzx_n,
     .    dsigt_dlam,dsigc_dlam,dsigs_dlam,da1_dsigs,da1_dsigt,
     .    da1_dsigc,da2_dsigs,da2_dsigt,da2_dsigc,da0_dlam,da1_dlam,
     .    da2_dlam,dfdsig2,dphi_dlam,dpxx,dpyy,dpzz,dpxy,dpyz,dpzx,
     .    dpdt_c,dpdt_s,dpdt_t,dpdt,dlam,cc,epdt_min,
     .    epdt_max,epdc_max,epdc_min,epds_min,epds_max,asrate
      my_real ,DIMENSION(NEL,2)   :: xvec
 
      my_real ,DIMENSION(NUMTABL) :: tfac
      my_real ,DIMENSION(NFUNC)   :: yfac
      LOGICAL :: CONV(NEL)
      my_real, PARAMETER :: SFAC = 1.05d0 ! Security factor of ICONV
C
      !=======================================================================
      ! - INITIALIZATION OF COMPUTATION ON TIME STEP
      !=======================================================================
      ! Recovering model parameters
      ! -> Elastic parameters
      g     = uparam(4) ! Shear modulus                
      aa1   = uparam(6) ! First component of elastic matrix
      aa2   = uparam(7) ! Second component of elastic matrix
      c1    = uparam(8) ! Bulk modulus
      ! -> Plastic parameters
      nupc  = uparam(9)  ! Plastic Poisson ratio                   
      ! -> Flags 
      iconv  = nint(uparam(15)) ! Flag to ensure convexity
      asrate = min(one,uparam(16)*timestep)
      icas   = nint(uparam(17)) ! Flag for tabulated case
      !icas      ifunt   | ifunc   | ifuncs
      !  -1         1    |    1    |    1
      !   0         1    |    0    |    0
      !   1         1    |    1    |    0
      !   2         1    |    0    |    1
      xfac  = uparam(18) ! Strain-rate scale factor
      epdt_min = uparam(19)
      epdt_max = uparam(20)
      epdc_min = uparam(21)
      epdc_max = uparam(22)
      epds_min = uparam(23)
      epds_max = uparam(24)
      tfac(1)  = uparam(25)
      tfac(2)  = uparam(26)
      tfac(3)  = uparam(27)
      yfac(1)  = uparam(28)
      yfac(2)  = uparam(29)
c
      ! Recovering internal variables
      DO i=1,nel    
        IF(off(i) < em01) off(i) = zero
        IF(off(i) <  one) off(i) = off(i)*four_over_5    
        plat(i)  = uvar(i,1) ! Uniaxial tension plastic strain    
        plac(i)  = uvar(i,2) ! Uniaxial compression plastic strain    
        plas(i)  = uvar(i,3) ! Shear plastic strain
        epspt(i) = uvar(i,4) ! Uniaxial tension plastic strain-rate
        epspc(i) = uvar(i,5) ! Uniaxial compression plastic strain-rate
        epsps(i) = uvar(i,6) ! Shear plastic strain-rate
        epspt(i) = min(epdt_max, max(epspt(i),epdt_min))
        epspc(i) = min(epdc_max, max(epspc(i),epdc_min))
        epsps(i) = min(epds_max, max(epsps(i),epds_min))
        dpla(i)  = zero
        dplac(i) = zero
        dplat(i) = zero
        dplas(i) = zero
      ENDDO
C     Initialize plastic Poisson ratio
      IF (time == zero) THEN
        nup(1:nel)    = nupc
        uvar(1:nel,7) = nupc
      ELSE
        nup(1:nel) = uvar(1:nel,7)
      END IF   
C
      ! Computation of yield stresses
      xvec(1:nel,1)   = plat(1:nel)
      xvec(1:nel,2)   = epspt(1:nel)*xfac
      CALL table_mat_vinterp(table(func_trac),nel,nel,vartmp(1,1),xvec,sigt,dsigt_dp)
      sigt     = sigt*tfac(1)
      dsigt_dp = dsigt_dp*tfac(1)
      IF (table(func_comp)%NOTABLE > 0) THEN 
        xvec(1:nel,1)   = plac(1:nel)
        xvec(1:nel,2)   = epspc(1:nel)*xfac
        CALL table_mat_vinterp(table(func_comp),nel,nel,vartmp(1,3),xvec,sigc,dsigc_dp)
        sigc     = sigc*tfac(2)
        dsigc_dp = dsigc_dp*tfac(2)
      ENDIF
      IF (table(func_shear)%NOTABLE > 0) THEN 
        xvec(1:nel,1)   = plas(1:nel)
        xvec(1:nel,2)   = epsps(1:nel)*xfac
        CALL table_mat_vinterp(table(func_shear),nel,nel,vartmp(1,5),xvec,sigs,dsigs_dp)
        sigs     = sigs*tfac(3)
        dsigs_dp = dsigs_dp*tfac(3) 
      ENDIF 
      ! Select case for tabulated yield stresses
      IF (icas == 0) THEN 
        sigc(1:nel) = sigt(1:nel)
        sigs(1:nel) = sigt(1:nel)/sqr3
      ELSEIF (icas == 1) THEN 
        DO i = 1,nel
          sigs(i) = sqrt(sigc(i)*sigt(i)/three)
        ENDDO
      ENDIF 
C
      !========================================================================
      ! - COMPUTATION OF TRIAL VALUES
      !========================================================================      
      DO i=1,nel
C        
        ! Computation of the trial stress tensor
        signxx(i) = sigoxx(i) + (aa1*depsxx(i)
     .                        + aa2*(depsyy(i) + depszz(i)))
        signyy(i) = sigoyy(i) + (aa1*depsyy(i)
     .                        + aa2*(depsxx(i) + depszz(i)))
        signzz(i) = sigozz(i) + (aa1*depszz(i)
     .                        + aa2*(depsxx(i) + depsyy(i)))
        signxy(i) = sigoxy(i) + g*depsxy(i)
        signyz(i) = sigoyz(i) + g*depsyz(i)
        signzx(i) = sigozx(i) + g*depszx(i)
C         
        ! Computation of the pressure of the trial stress tensor
        p(i) = -third*(signxx(i) + signyy(i) + signzz(i))
        ! Computation of the Von Mises equivalent stress of the trial stress tensor
        sdxx(i) = signxx(i) + p(i)
        sdyy(i) = signyy(i) + p(i)
        sdzz(i) = signzz(i) + p(i)
        sdxy(i) = signxy(i)
        sdyz(i) = signyz(i)
        sdzx(i) = signzx(i)
        svm(i)  = half*(sdxx(i)**2 + sdyy(i)**2 + sdzz(i)**2)
     .       +           (sdxy(i)**2 + sdzx(i)**2 + sdyz(i)**2)
        svm(i)  = sqrt(three*svm(i))
      ENDDO
C      
      !========================================================================
      ! - COMPUTATION OF YIELD FONCTION
      !======================================================================== 
      ! Compute yield criterion parameters
      IF (iconv == 1) THEN
        ! Ensured convexity 
        DO i = 1,nel
          conv(i) = .false.
          IF (sigs(i) < sfac*sqrt(sigc(i)*sigt(i)/three)) THEN 
            sigs(i) = sfac*sqrt(sigc(i)*sigt(i)/three)
            conv(i) = .true.
          ENDIF
        ENDDO   
      ENDIF
      ! Compute yield criterion parameters A0,A1,A2
      DO i = 1,nel
        aa    = one/sigc(i)/sigt(i)
        a0(i) = three*(sigs(i)**2)
        a1(i) = nine*(sigs(i)**2)*(sigc(i) - sigt(i))*aa
        a2(i) = nine*(sigc(i)*sigt(i) - three*(sigs(i)**2))*aa
      ENDDO   
c
      ! -> Checking plastic behavior for all elements
      nindx = 0
      DO i=1,nel
        phi(i) = (svm(i)**2) - a0(i) - a1(i)*p(i) - a2(i)*p(i)*p(i)
        IF (phi(i) >= zero .AND. off(i) == one) THEN
          nindx = nindx + 1
          indx(nindx) = i
        ENDIF
      ENDDO
C      
      !====================================================================
      ! - PLASTIC RETURN MAPPING WITH CUTTING PLANE METHOD
      !====================================================================
      IF (nindx > 0) THEN
C  
        ! Loop over the iterations   
        DO iter = 1,3
c
          ! Loop over yielding elements
          DO ii = 1,nindx         
            i = indx(ii)
C        
            ! Note     : in this part, the purpose is to compute for each iteration
            ! a plastic multiplier allowing to update internal variables to satisfy
            ! the consistency condition using the cutting plane algorithm
            ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
            ! -> PHI          : current value of yield function (known)
            ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
            !                   into account of internal variables kinetic : 
            !                   plasticity ... 
C        
            ! 1 - Computation of DPHI_DSIG the normal to the yield surface
            !-------------------------------------------------------------
C                  
            ! Computation of normal to yield surface 
            normxx   = three*sdxx(i) + third*(a1(i) + two*a2(i)*p(i)) 
            normyy   = three*sdyy(i) + third*(a1(i) + two*a2(i)*p(i))
            normzz   = three*sdzz(i) + third*(a1(i) + two*a2(i)*p(i))
            normxy   = three*sdxy(i)
            normyz   = three*sdyz(i)
            normzx   = three*sdzx(i)
C
            ! Computation of the normalized normal to the yield surface
            cb       = a1(i) + two*a2(i)*p(i)                                 
            norm     = max(em20,sqrt(six*(svm(i)**2) + third*cb*cb))     
            normxx_n = normxx/norm
            normyy_n = normyy/norm
            normzz_n = normzz/norm
            normxy_n = normxy/norm
            normyz_n = normyz/norm
            normzx_n = normzx/norm
C                     
            ! 2 - Computation of DPHI_DLAMBDA
            !---------------------------------------------------------
C        
            !   a) Derivative with respect stress increments tensor DSIG
            !   --------------------------------------------------------
            dfdsig2 = normxx*(aa1*normxx_n + aa2*(normyy_n + normzz_n)) +
     .                normyy*(aa1*normyy_n + aa2*(normxx_n + normzz_n)) +
     .                normzz*(aa1*normzz_n + aa2*(normxx_n + normyy_n)) + 
     .                two*normxy*two*g*normxy_n + 
     .                two*normyz*two*g*normyz_n + 
     .                two*normzx*two*g*normzx_n
c
            !   b) Derivative of yield criterion parameters
            !   --------------------------------------------
C
            ! Derivative of yield surfaces with respect to yield criterion parameter A0,A1,A2                                                                  
            dsigt_dlam = dsigt_dp(i)*(svm(i)/norm)*(three/(one + nup(i)))
            IF (table(func_comp)%NOTABLE  > 0) dsigc_dlam = dsigc_dp(i)*(svm(i)/norm)*(three/(one + nup(i)))
            IF (table(func_shear)%NOTABLE > 0) dsigs_dlam = dsigs_dp(i)*sqr3*(svm(i)/norm)
            IF (icas == 0) THEN 
              dsigc_dlam = dsigt_dlam 
              dsigs_dlam = (one/sqr3)*dsigt_dlam
            ELSEIF (icas == 1) THEN 
              dsigs_dlam = (one/sqr3)*(one/(two*sqrt(sigt(i)*sigc(i))))*
     .                     (dsigc_dlam*sigt(i) + sigc(i)*dsigt_dlam)
            ENDIF    
            IF (iconv == 1) THEN                                      
              IF (conv(i)) THEN 
                dsigs_dlam = sfac*(one/sqr3)*(one/(two*sqrt(sigt(i)*sigc(i))))*
     .                       (dsigc_dlam*sigt(i) + sigc(i)*dsigt_dlam)             
              ENDIF 
            ENDIF    
C
            ! -> A1 derivatives 
            cc = sigs(i)/sigc(i)/sigt(i)
            !   -> With respect to SIGS
            da1_dsigs = eighteen*(sigc(i) - sigt(i))*cc
            !   -> With respect to SIGC
            da1_dsigc = nine*(sigs(i)/sigc(i))**2
            !   -> With respect to SIGT
            da1_dsigt = -nine*(sigs(i)/sigt(i))**2
C
            ! -> A2 derivatives
            !   -> With respect to SIGS
            da2_dsigs = -fifty4*cc
            !   -> With respect to SIGC                                         
            da2_dsigc = twenty7*cc*sigs(i)/sigc(i)  
            !   -> With respect to SIGT                       
            da2_dsigt = twenty7*cc*sigs(i)/sigt(i)                     
C            
            ! -> A parameters derivatives with respect to lambda        
            !   -> A0 with respect to lambda                                               
            da0_dlam = six*sigs(i)*dsigs_dlam      
            !   -> A1 with respect to lambda                            
            da1_dlam = da1_dsigs*dsigs_dlam + da1_dsigt*dsigt_dlam + da1_dsigc*dsigc_dlam
            !   -> A2 with respect to lambda                     
            da2_dlam = da2_dsigs*dsigs_dlam + da2_dsigt*dsigt_dlam + da2_dsigc*dsigc_dlam                                                
C          
C            
            ! 3 - Computation of plastic multiplier and variables update
            !----------------------------------------------------------
C          
            ! Derivative of PHI with respect to DLAM
            dphi_dlam = - dfdsig2 - da0_dlam - p(i)*da1_dlam - (p(i)**2)*da2_dlam   
            dphi_dlam = sign(max(abs(dphi_dlam),em20),dphi_dlam)                 
            dlam = -phi(i)/dphi_dlam    
C          
            ! Plastic strains tensor update
            dpxx = dlam*normxx_n                          
            dpyy = dlam*normyy_n                          
            dpzz = dlam*normzz_n                          
            dpxy = dlam*normxy_n                          
            dpyz = dlam*normyz_n                          
            dpzx = dlam*normzx_n                 
C          
            ! Elasto-plastic stresses update   
            signxx(i) = signxx(i) - (aa1*dpxx + aa2*(dpyy + dpzz)) 
            signyy(i) = signyy(i) - (aa1*dpyy + aa2*(dpxx + dpzz))
            signzz(i) = signzz(i) - (aa1*dpzz + aa2*(dpxx + dpyy))
            signxy(i) = signxy(i) - two*g*dpxy     
            signyz(i) = signyz(i) - two*g*dpyz  
            signzx(i) = signzx(i) - two*g*dpzx
C          
            ! cumulated plastic strains update       
            plat(i)  = plat(i) + dlam*(svm(i)/norm)*three/(one + nup(i))
            plac(i)  = plat(i)    
            plas(i)  = plas(i) + sqr3*(svm(i)/norm)*dlam   
            pla(i)   = pla(i)  + two*dlam*(svm(i)/norm)
            ! Plastic strain increments update
            dplat(i) = dplat(i) + dlam*(svm(i)/norm)*three/(one + nup(i))
            dplac(i) = dplat(i)
            dplas(i) = dplas(i) + sqr3*(svm(i)/norm)*dlam
            dpla(i)  = dpla(i)  + two*dlam*(svm(i)/norm)
C
            ! Update pressure
            p(i) = -third*(signxx(i) + signyy(i) + signzz(i))
            ! Update Von Mises stress
            sdxx(i) = signxx(i) + p(i)
            sdyy(i) = signyy(i) + p(i)
            sdzz(i) = signzz(i) + p(i)
            sdxy(i) = signxy(i)
            sdyz(i) = signyz(i)
            sdzx(i) = signzx(i)
            svm(i)  = half*(sdxx(i)**2 + sdyy(i)**2 + sdzz(i)**2)
     .         +           (sdxy(i)**2 + sdzx(i)**2 + sdyz(i)**2)
            svm(i)  = sqrt(three*svm(i))
          ENDDO
C
          ! Update yield stresses
          xvec(1:nel,1)   = plat(1:nel)
          xvec(1:nel,2)   = epspt(1:nel)*xfac
          CALL table_mat_vinterp(table(func_trac),nel,nel,vartmp(1,1),xvec,sigt,dsigt_dp)
          sigt     = sigt*tfac(1)
          dsigt_dp = dsigt_dp*tfac(1)
          IF (table(func_comp)%NOTABLE > 0) THEN 
            xvec(1:nel,1)   = plac(1:nel)
            xvec(1:nel,2)   = epspc(1:nel)*xfac
            CALL table_mat_vinterp(table(func_comp),nel,nel,vartmp(1,3),xvec,sigc,dsigc_dp)
            sigc(1:nel)     = sigc*tfac(2)
            dsigc_dp(1:nel) = dsigc_dp*tfac(2)
          ENDIF
          IF (table(func_shear)%NOTABLE > 0) THEN 
            xvec(1:nel,1)   = plas(1:nel)
            xvec(1:nel,2)   = epsps(1:nel)*xfac
            CALL table_mat_vinterp(table(func_shear),nel,nel,vartmp(1,5),xvec,sigs,dsigs_dp)
            sigs(1:nel)     = sigs*tfac(3)
            dsigs_dp(1:nel) = dsigs_dp*tfac(3)
          ENDIF 
          ! Select case for tabulated yield stresses
          IF (icas == 0) THEN 
            DO ii = 1,nindx
              i = indx(ii)
              sigc(i) = sigt(i)
              sigs(i) = sigt(i)/sqr3
            ENDDO
          ELSEIF (icas == 1) THEN 
            DO ii = 1,nindx
              i = indx(ii)
              sigs(i) = sqrt(sigc(i)*sigt(i)/three)
            ENDDO
          ENDIF
C
          ! Update yield function value
          IF (iconv == 1) THEN
            DO ii = 1,nindx         
              i = indx(ii)
              conv(i) = .false.
              IF (sigs(i) < sfac*sqrt(sigc(i)*sigt(i)/three)) THEN 
                sigs(i) = sfac*sqrt(sigc(i)*sigt(i)/three)
                conv(i) = .true.
              ENDIF
            ENDDO
          ENDIF
          DO ii = 1,nindx         
            i = indx(ii) 
            aa = one/sigc(i)/sigt(i)
            a0(i) = three*(sigs(i)**2)
            a1(i) = nine*(sigs(i)**2)*(sigc(i) - sigt(i))*aa
            a2(i) = nine*(sigc(i)*sigt(i) - three*(sigs(i)**2))*aa
            phi(i) = svm(i)**2 - a0(i) - a1(i)*p(i) - a2(i)*p(i)*p(i)  
          ENDDO
        ENDDO 
      ENDIF
c-----------------------------------------------
C     Update plastic Poisson coefficient
c-----------------------------------------------
      IF (ifunc(1) > 0) THEN
        iad(1:nel)  = npf(ifunc(1))   / 2 + 1
        ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad(1:nel) - vartmp(1:nel,8)
!
        CALL vinter(tf,iad,vartmp(1:nel,8),ilen,nel,pla,df,nup)
!        
        uvar(1:nel,7) = yfac(1) * nup(1:nel)
        uvar(1:nel,7) = max(zero, min(nup(1:nel), half))
      END IF
C      
      !====================================================================
      ! - STORING NEW VALUES
      !====================================================================
      DO i=1,nel
        uvar(i,1) = plat(i)          
        uvar(i,2) = plac(i)          
        uvar(i,3) = plas(i)    
        dpdt_t    = dplat(i)/max(timestep,em20)
        uvar(i,4) = asrate*dpdt_t + (one-asrate)*epspt(i)   
        dpdt_c    = dplac(i)/max(timestep,em20)        
        uvar(i,5) = asrate*dpdt_c + (one-asrate)*epspc(i) 
        dpdt_s    = dplas(i)/max(timestep,em20)          
        uvar(i,6) = asrate*dpdt_s + (one-asrate)*epsps(i)
        dpdt      = dpla(i)/max(timestep,em20) 
        epsd(i)   = asrate*dpdt   + (one-asrate)*epsd(i)
        ! Computation of soundspeed
        soundsp(i) = sqrt((c1+four*g/three)/rho0(i)) 
        ! Computation of the yield stress
        yld(i)    = a0(i) + a1(i)*p(i) + a2(i)*p(i)*p(i)
        ! Computation of the hourglass coefficient
        et(i)     = half
      ENDDO         
C         
      END