!======================================================================!
!                                                                      !
! Software Name : FrontISTR Ver. 3.6                                   !
!                                                                      !
!      Module Name : lib                                               !
!                                                                      !
!                    Written by Xi YUAN (AdavanceSoft)                 !
!                                                                      !
!      Contact address :  IIS,The University of Tokyo, CISS            !
!                                                                      !
!      "Structural Analysis for Large Scale Assembly"                  !
!                                                                      !
!======================================================================!
!======================================================================!
!>  \brief   This module provide functions for elastoplastic calculation
!!
!>  \author     Xi YUAN (AdavanceSoft)
!>  \date       2010/01/12
!>  \version    0.00
!======================================================================!
module m_ElastoPlastic
  use mMaterial
  use m_ElasticLinear

  implicit none
  integer, parameter, private :: kreal = kind(0.0d0)

  contains

   !> This subroutine calculates elastoplastic constitutive relation
   SUBROUTINE calElastoPlasticMatrix( matl, sectType, stress, istat, extval, D, temperature )
     TYPE( tMaterial ), INTENT(IN) :: matl      !< material properties
     INTEGER, INTENT(IN)           :: sectType  !< not used currently
     REAL(KIND=kreal), INTENT(IN)  :: stress(6) !< stress
     REAL(KIND=kreal), INTENT(IN)  :: extval(:) !< plastic strain, back stress
     INTEGER, INTENT(IN)           :: istat     !< plastic state
     REAL(KIND=kreal), INTENT(OUT) :: D(:,:)    !< constitutive relation
     REAL(KIND=kreal), OPTIONAL    :: temperature   !> temprature

     INTEGER :: i,j,ytype
     logical :: kinematic
     REAL(KIND=kreal) :: dum, dj1(6), dj2(6), dj3(6), a(6), De(6,6)
     REAL(KIND=kreal) :: C1,C2,C3, back(6)
     REAL(KIND=kreal) :: J1,J2,J3, fai, sita, harden, khard, da(6), devia(6)

	 ytype = getYieldFunction( matl%mtype )
	 if( ytype==3 ) THEN
       ! Gaku Hashimoto (2015/08/28) <
       ! call uElastoPlasticMatrix( matl, stress, istat, extval, D  )
       ! > Gaku Hashimoto (2015/08/28)
       return
     endif
     IF( sectType /=D3 ) stop "Elastoplastic calculation support only Solid element currently"
     kinematic = isKinematicHarden( matl%mtype )
     khard = 0.d0
     if( kinematic ) then
       back(1:6) = extval(1:6)
       khard = calKinematicHarden( matl, extval(1) )
     endif
     if( present( temperature ) ) then
          call calElasticMatrix( matl, sectTYPE, De, temperature  )
     else
          call calElasticMatrix( matl, sectTYPE, De )
     endif

     J1 = (stress(1)+stress(2)+stress(3))
     devia(1:3) = stress(1:3)-J1/3.d0
     devia(4:6) = stress(4:6)
     if( kinematic ) devia = devia-back
     J2 = 0.5d0* dot_product( devia(1:3), devia(1:3) ) +  &
          dot_product( devia(4:6), devia(4:6) )

     D(:,:) = De(:,:)
     if( istat == 0 ) return   ! elastic state

     !derivative of J2
     dj2(1:3) = devia(1:3)
     dj2(4:6) = 2.d0*devia(4:6)
     dj2 = dj2/( 2.d0*dsqrt(j2) )
     if( present(temperature) ) then
       harden = calHardenCoeff( matl, extval(1), temperature )
     else
       harden = calHardenCoeff( matl, extval(1) )
     endif

     SELECT CASE (yType)
     CASE (0)       ! Mises or. Isotropic
       a(1:6) = dsqrt(3.d0) * dj2
     CASE (1)      ! Mohr-Coulomb
       fai = matl%variables(M_PLCONST3)
       J3 = devia(1)*devia(2)*devia(3)                    &
           +2.d0* devia(4)*devia(5)*devia(6)                     &
           -devia(6)*devia(2)*devia(6)                           &
           -devia(4)*devia(4)*devia(3)                           &
           -devia(1)*devia(5)*devia(5)
       sita = -3.d0*dsqrt(3.d0)*J3/( 2.d0*(J2**1.5d0) )
       if( dabs( dabs(sita)-1.d0 ) <1.d-8 ) THEN
          C1 = 0.d0
          C2 = dsqrt(3.d0)
          C3 = 0.d0
       else
         if( dabs(sita) >1.d0 ) stop "Math Error in Mohr-Coulomb calculation"
         sita = asin( sita )/3.d0
         C2 = cos(sita)*( 1.d0*tan(sita)*tan(3.d0*sita) + sin(fai)* &
             ( tan(3.d0*sita)-tan(sita )/dsqrt(3.d0) ) )
         C1 = sin(fai)/3.d0
         C3 = dsqrt(3.d0)*sin(sita)+cos(sita)*sin(fai)/(2.d0*J2*cos(3.d0*sita))
       endif
       ! deirivative of j1
       dj1(1:3) = 1.d0
       dj1(4:6) = 0.d0
       ! deirivative of j3
       dj3(1) = devia(2)*devia(3)-devia(5)*devia(5)+J2/3.d0
       dj3(2) = devia(1)*devia(3)-devia(6)*devia(6)+J2/3.d0
       dj3(3) = devia(1)*devia(2)-devia(4)*devia(4)+J2/3.d0
       dj3(4) = 2.d0*(devia(5)*devia(6)-devia(3)*devia(4))
       dj3(5) = 2.d0*(devia(4)*devia(6)-devia(1)*devia(5))
       dj3(6) = 2.d0*(devia(4)*devia(5)-devia(2)*devia(6))
       a(:) = C1*dj1 + C2*dj2 + C3*dj3
     CASE (2)      ! Drucker-Prager
       fai = matl%variables(M_PLCONST3)
       ! deirivative of j1
       dj1(1:3) = 1.d0
       dj1(4:6) = 0.d0
       a(:) = fai*dj1(:) + dj2(:)
     END SELECT

     da = matmul( de, a )
     dum = harden + khard+ dot_product( da, a )
     DO i=1,6
     DO j=1,6
       D(i,j) = De(i,j) - da(i)*da(j)/dum
     ENDDO
     ENDDO

   END SUBROUTINE calElastoPlasticMatrix
      
      
      ! Gaku Hashimoto (2015/07/20) <
!####################################################################
      SUBROUTINE calElastoPlasticMatrix_modified(icel, LX, matl, sectType, stress, delta_gamma, apex, istat, extval, strain_e_trial, D, temp)
!####################################################################
      
      INTEGER, INTENT(IN)           :: icel        !> element number
      INTEGER, INTENT(IN)           :: LX          !> Gauss point number
      TYPE( tMaterial ), INTENT(IN) :: matl        !< material properties
      INTEGER, INTENT(IN)           :: sectType    !< not used currently
      REAL(KIND=kreal), INTENT(IN)  :: stress(6)   !< stress
      REAL(KIND=kreal), INTENT(IN)  :: delta_gamma !< incremental plastic multiplier
      LOGICAL, INTENT(IN)           :: apex        !< Apex return
      INTEGER, INTENT(IN)           :: istat       !< plastic state
      REAL(KIND=kreal), INTENT(IN)  :: extval(:)   !< plastic strain, back stress
      ! Gaku Hashimoto (2015/07/20) <
      REAL(KIND=kreal), INTENT(IN)  :: strain_e_trial(6) !> Elastic trial strain
      ! > Gaku Hashimoto (2015/07/20)
      REAL(KIND=kreal), INTENT(OUT) :: D(:,:)      !< constitutive relation
      REAL(KIND=kreal), OPTIONAL    :: temp        !< temprature
      
      INTEGER :: ytype
      INTEGER :: is, js
      INTEGER :: EPFLAG
      
      REAL(KIND=kreal) :: foid(6, 6)
      REAL(KIND=kreal) :: soid(6)
      REAL(KIND=kreal) :: devprj(6, 6)
      REAL(KIND=kreal) :: ina(1), outa(2)
      REAL(KIND=kreal) :: e, nu
      REAL(KIND=kreal) :: kappa, g
      REAL(KIND=kreal) :: p
      REAL(KIND=kreal) :: stess_norm
      REAL(KIND=kreal) :: q_trial
      REAL(KIND=kreal) :: back(6)
      REAL(KIND=kreal) :: j2
      REAL(KIND=kreal) :: harden, khard
      REAL(KIND=kreal) :: stress_dev(6)
      REAL(KIND=kreal) :: coef
      REAL(KIND=kreal) :: coef1, coef2, coef3, coef4
      REAL(KIND=kreal) :: xi, eta
      REAL(KIND=kreal) :: tr_strain_e_trial
      REAL(KIND=kreal) :: strain_e_trial_dev(6)
      REAL(KIND=kreal) :: norm_strain_e_trial_dev, norm_strain_e_trial_dev_inv
      REAL(KIND=kreal) :: unidev(6)
      REAL(KIND=kreal) :: a
      
      LOGICAL :: ierr
      LOGICAL :: kinematic
      
!--------------------------------------------------------------------
      
      IF( sectType .NE. D3 ) STOP "Elastoplastic calculation support only Solid element currently"
      
!--------------------------------------------------------------------
      
      foid = 0.0D0
      
      foid(1, 1) = 1.0D0
      foid(2, 2) = 1.0D0
      foid(3, 3) = 1.0D0
      foid(4, 4) = 0.5D0
      foid(5, 5) = 0.5D0
      foid(6, 6) = 0.5D0
      
      !--------------------------------------------------------------
      
      soid = 0.0D0
      
      soid(1) = 1.0D0
      soid(2) = 1.0D0
      soid(3) = 1.0D0
      
      !--------------------------------------------------------------
      
      DO js = 1, 6
       
       DO is = 1, 6
        
        devprj(is, js) = foid(is, js)-soid(is)*soid(js)/3.0D0
        
       END DO
       
      END DO
      
!--------------------------------------------------------------------
      
      IF( PRESENT( temp ) ) THEN
      
       ina(1) = temp
       
       CALL fetch_TableData(MC_ISOELASTIC, matl%dict, outa, ierr, ina)
       
       IF( ierr ) THEN
        
        e  = matl%variables(M_YOUNGS)
        nu = matl%variables(M_POISSON)
        
       ELSE
        
        e  = outa(1)
        nu = outa(2)
        
       END IF
       
      ELSE
       
       CALL fetch_TableData(MC_ISOELASTIC, matl%dict, outa, ierr)
       
       IF( ierr ) THEN
        
        e  = matl%variables(M_YOUNGS)
        nu = matl%variables(M_POISSON)
        
       ELSE
        
        e  = outa(1)
        nu = outa(2)
        
       END IF
       
      END IF
      
      !--------------------------------------------------------------
      
      kappa = e/( 3.0D0*(1.0D0-2.0D0*nu) )
      g = e/( 2.0D0*(1.0D0+nu) )
      
!--------------------------------------------------------------------
      
      IF( istat .EQ. 0 ) THEN
       
       EPFLAG = 0
       
      ELSE IF( istat .EQ. 1 ) THEN
       
       ! Elasto-plastic tangent flag
       EPFLAG = 1
       
      END IF
      
!--------------------------------------------------------------------
      
      IF( EPFLAG .EQ. 1 ) THEN
       
       !--------------------------------------------------------
       
       IF( PRESENT(temp) ) THEN
        
        harden = calHardenCoeff( matl, extval(1), temp )
        
       ELSE
        
        harden = calHardenCoeff( matl, extval(1) )
        
       END IF
       
       !--------------------------------------------------------
       
       kinematic = isKinematicHarden( matl%mtype )
       
       khard = 0.0D0
       
       IF( kinematic ) THEN
        
        back(1) = extval(1)
        back(2) = extval(2)
        back(3) = extval(3)
        back(4) = extval(4)
        back(5) = extval(5)
        back(6) = extval(6)
        
        khard = calKinematicHarden( matl, extval(1) )
        
       END IF
       
       !--------------------------------------------------------
       
       ytype = getYieldFunction(matl%mtype)
       
       !--------------------------------------------------------
       
       ! von Mises
       IF( ytype .EQ. 0 ) THEN
        
        !--------------------------------------------------
        
        ! Hydrostatic pressure
        p = ( stress(1)+stress(2)+stress(3) )/3.0D0
        
        ! Deviatoric stress components
        stress_dev(1) = stress(1)-p
        stress_dev(2) = stress(2)-p
        stress_dev(3) = stress(3)-p
        stress_dev(4) = stress(4)
        stress_dev(5) = stress(5)
        stress_dev(6) = stress(6)
        
        IF( kinematic ) THEN
         
         stress_dev(1) = stress_dev(1)-back(1)
         stress_dev(2) = stress_dev(2)-back(2)
         stress_dev(3) = stress_dev(3)-back(3)
         stress_dev(4) = stress_dev(4)-back(4)
         stress_dev(5) = stress_dev(5)-back(5)
         stress_dev(6) = stress_dev(6)-back(6)
         
        END IF
        
        !--------------------------------------------------
        
        j2 = 0.5D0*( stress_dev(1)*stress_dev(1)         &
                    +stress_dev(2)*stress_dev(2)         &
                    +stress_dev(3)*stress_dev(3)         &
                    +2.0D0*stress_dev(4)*stress_dev(4)   &
                    +2.0D0*stress_dev(5)*stress_dev(5)   &
                    +2.0D0*stress_dev(6)*stress_dev(6) ) 
        
        stess_norm = DSQRT( 2.0D0*j2 )
        
        q_trial = DSQRT( 3.0D0*j2 )+3.0D0*g*delta_gamma
        
        coef1 = 2.0D0*g*( 1.0D0-3.0D0*g*delta_gamma/q_trial )
        coef2 = 6.0D0*g*g                                             &
                *( delta_gamma/q_trial-1.0D0/(3.0D0*g+harden+khard) ) &
                /( stess_norm*stess_norm )                            
        
        !--------------------------------------------------
        
        DO js = 1, 6
         
         DO is = 1, 6
          
          d(is, js) = kappa*soid(is)*soid(js)             &
                     +coef1*devprj(is, js)                &
                     +coef2*stress_dev(is)*stress_dev(js) 
          
         END DO
         
        END DO
        
        !--------------------------------------------------
        
       ! Drucker-Prager
       ELSE IF( ytype .EQ. 2 ) THEN
        
        !--------------------------------------------------
        
        eta = 3.0D0*matl%variables(M_PLCONST3)
        xi  = matl%variables(M_PLCONST4)
        
        !--------------------------------------------------
        
        ! Apex return
        IF( apex ) THEN
         
         !--------------------------------------------
         
         coef = kappa*( 1.0D0-kappa/( kappa+( xi/eta )*( xi/eta )*harden ) )
         
         DO is = 1, 6
          
          DO js = 1, 6
           
           d(is, js) = coef*soid(is)*soid(js)
           
          END DO
          
         END DO
         
         !--------------------------------------------
         
        ! Smooth cone wall return
        ELSE
         
         !--------------------------------------------
         
         tr_strain_e_trial = strain_e_trial(1)+strain_e_trial(2)+strain_e_trial(3)
         
         strain_e_trial_dev(1) = strain_e_trial(1)-tr_strain_e_trial/3.0D0
         strain_e_trial_dev(2) = strain_e_trial(2)-tr_strain_e_trial/3.0D0
         strain_e_trial_dev(3) = strain_e_trial(3)-tr_strain_e_trial/3.0D0
         strain_e_trial_dev(4) = strain_e_trial(4)
         strain_e_trial_dev(5) = strain_e_trial(5)
         strain_e_trial_dev(6) = strain_e_trial(6)
         
         norm_strain_e_trial_dev                                      &
         = DSQRT( strain_e_trial_dev(1)*strain_e_trial_dev(1)         &
                 +strain_e_trial_dev(2)*strain_e_trial_dev(2)         &
                 +strain_e_trial_dev(3)*strain_e_trial_dev(3)         &
                 +0.5D0*strain_e_trial_dev(4)*strain_e_trial_dev(4)   &
                 +0.5D0*strain_e_trial_dev(5)*strain_e_trial_dev(5)   &
                 +0.5D0*strain_e_trial_dev(6)*strain_e_trial_dev(6) ) 
         
         IF( norm_strain_e_trial_dev .NE. 0.0D0 ) THEN
          
          norm_strain_e_trial_dev_inv = 1.0D0/norm_strain_e_trial_dev
          
         ELSE
          
          norm_strain_e_trial_dev_inv = 0.0D0
          
         END IF
         
         unidev(1) = norm_strain_e_trial_dev_inv*strain_e_trial_dev(1)
         unidev(2) = norm_strain_e_trial_dev_inv*strain_e_trial_dev(2)
         unidev(3) = norm_strain_e_trial_dev_inv*strain_e_trial_dev(3)
         unidev(4) = norm_strain_e_trial_dev_inv*0.5D0*strain_e_trial_dev(4)
         unidev(5) = norm_strain_e_trial_dev_inv*0.5D0*strain_e_trial_dev(5)
         unidev(6) = norm_strain_e_trial_dev_inv*0.5D0*strain_e_trial_dev(6)
         
         a = 1.0D0/( g+kappa*eta*eta+xi*xi*( harden+khard ) )
         
         coef1 = 2.0D0*g                                                           &
                 *( 1.0D0-delta_gamma/( DSQRT( 2.0D0 )*norm_strain_e_trial_dev ) ) 
         coef2 = 2.0D0*g                                                         &
                 *( delta_gamma/( DSQRT( 2.0D0 )*norm_strain_e_trial_dev )-g*a ) 
         coef3 = DSQRT( 2.0D0 )*g*a*kappa
         coef4 = kappa*( 1.0D0-kappa*eta*eta*a )
         
         DO is = 1, 6
          
          DO js = 1, 6
            
           d(is, js) = coef1*devprj(is, js)                                      &
                      +coef2*unidev(is)*unidev(js)                               &
                      -coef3*( eta*unidev(is)*soid(js)+eta*soid(is)*unidev(js) ) &
                      +coef4*soid(is)*soid(js)                                   
           
          END DO
          
         END DO
         
         !--------------------------------------------
         
        END IF
        
        !--------------------------------------------------
        
       ! User subroutine
       ELSE IF( ytype .EQ. 3 ) THEN
        
        !--------------------------------------------------
        
        ! Gaku Hashimoto (2015/08/28) <
        ! CALL uElastoPlasticMatrix(matl, stress, istat, extval, D)
        ! > Gaku Hashimoto (2015/08/28)
        
        !--------------------------------------------------
        
       END IF
       
       !--------------------------------------------------------
       
      ELSE
       
       !--------------------------------------------------------
       
       DO js = 1, 6
        
        DO is = 1, 6
         
         d(is, js) = kappa*soid(is)*soid(js)+2.0D0*g*devprj(is, js)
         
        END DO
          
       END DO
       
       !--------------------------------------------------------
       
      END IF
      
!--------------------------------------------------------------------
      
      RETURN
      
!####################################################################
      END SUBROUTINE calElastoPlasticMatrix_modified
!####################################################################
      ! > Gaku Hashimoto (2015/07/20)
      
      
   !> This subrouitne calculate equivalent stress
   REAL(KIND=kreal) FUNCTION cal_equivalent_stress(matl, stress, extval)
     TYPE( tMaterial ), INTENT(IN) :: matl        !< material property
     REAL(kind=kreal), INTENT(IN)  :: stress(6)   !< stress
     REAL(KIND=kreal), INTENT(IN)  :: extval(:)   !< plastic strain, back stress

     INTEGER :: ytype
     logical :: kinematic
     REAL(kind=kreal) :: eqvs, sita, fai, J1,J2,J3, devia(6)
     REAL(KIND=kreal) :: back(6)
     kinematic = isKinematicHarden( matl%mtype )
     if( kinematic ) back(1:6) = extval(1:6)

     ytype = getYieldFunction( matl%mtype )
     J1 = (stress(1)+stress(2)+stress(3))
     devia(1:3) = stress(1:3)-J1/3.d0
     devia(4:6) = stress(4:6)
     if( kinematic ) devia = devia-back
     J2 = 0.5d0* dot_product( devia(1:3), devia(1:3) ) +  &
          dot_product( devia(4:6), devia(4:6) )

     SELECT CASE (yType)
     CASE (0)       ! Mises or. Isotropic
        eqvs = dsqrt( 3.d0*J2 )
     CASE (1)      ! Mohr-Coulomb
       fai = matl%variables(M_PLCONST1)
       J3 = devia(1)*devia(2)*devia(3)                    &
           +2.d0* devia(4)*devia(5)*devia(6)                     &
           -devia(6)*devia(2)*devia(6)                           &
           -devia(4)*devia(4)*devia(3)                           &
           -devia(1)*devia(5)*devia(5)
       sita = -3.d0*dsqrt(3.d0)*J3/( 2.d0*(J2**1.5d0) )
       if( dabs( dabs(sita)-1.d0 ) <1.d-8 ) sita=SIGN(1.d0, sita)
       if( dabs(sita) >1.d0 ) stop "Math Error in Mohr-Coulomb calculation"
       sita = asin( sita )/3.d0
       eqvs = (cos(sita)-sin(sita)*sin(fai)/dsqrt(3.d0))*dsqrt(J2)  &
             +J1*sin(fai)/3.d0
     CASE (2)      ! Drucker-Prager
       eqvs = dsqrt(J2)
     CASE DEFAULT
       eqvs = -1.d0
     END SELECT

     cal_equivalent_stress = eqvs
   END FUNCTION

   !> This subrouitne calculate equivalent stress
   REAL(KIND=kreal) FUNCTION cal_mises_strain( strain )
     REAL(kind=kreal), INTENT(IN)  :: strain(6)        !< strain
	 cal_mises_strain = 2.d0*dot_product( strain(1:3), strain(1:3) )
	 cal_mises_strain = cal_mises_strain+ dot_product( strain(4:6), strain(4:6) )
	 cal_mises_strain = dsqrt( cal_mises_strain/3.d0 )
   END FUNCTION

   !> This function calcualtes hardening coefficient
   REAL(KIND=kreal) FUNCTION calHardenCoeff( matl, pstrain, temp )
     TYPE( tMaterial ), INTENT(IN)          :: matl    !< material property
     REAL(KIND=kreal), INTENT(IN)           :: pstrain !< plastic strain
     REAL(KIND=kreal), INTENT(IN), OPTIONAL :: temp !< temprature

     INTEGER :: i, nc, htype
     logical :: ierr
     REAL(KIND=kreal) :: s0, s1,s2, ef, ina(2)

     calHardenCoeff = -1.d0
     htype = getHardenType( matl%mtype )
     SELECT CASE (htype)
     CASE (0)  ! Linear hardening
       calHardenCoeff = matl%variables(M_PLCONST2)
     CASE (1)  ! Multilinear approximation
       if( present(temp) ) then
          ina(1) = temp;  ina(2)=pstrain
          call fetch_TableGrad( MC_YIELD, ina, matl%dict, calHardenCoeff, ierr )
	!	  print *, ina, calHardenCoeff; pause
       else
         ina(1)=pstrain
         call fetch_TableGrad( MC_YIELD, ina(1:1), matl%dict, calHardenCoeff, ierr )
       endif
	 !  print *,1, calHardenCoeff; pause
     CASE (2)  ! Swift
       s0= matl%variables(M_PLCONST1)
       s1= matl%variables(M_PLCONST2)
       s2= matl%variables(M_PLCONST3)
       calHardenCoeff = s1*s2*( s0+pstrain )**(s2-1)
     CASE (3)  ! Ramberg-Osgood
       s0= matl%variables(M_PLCONST1)
       s1= matl%variables(M_PLCONST2)
       s2= matl%variables(M_PLCONST3)
       if( present(temp) ) then
         ef = calCurrYield( matl, pstrain, temp )
       else
         ef = calCurrYield( matl, pstrain )
       endif
       calHardenCoeff = s1*(ef/s1)**(1.d0-s2) /(s0*s2)
     CASE(4)   ! Prager
       calHardenCoeff = 0.d0
     CASE(5)   ! Prager+linear
       calHardenCoeff = matl%variables(M_PLCONST2)
     END SELECT
   END FUNCTION

   !> This function calcualtes kinematic hardening coefficient
   REAL(KIND=kreal) FUNCTION calKinematicHarden( matl, pstrain )
     TYPE( tMaterial ), INTENT(IN) :: matl    !< material property
     REAL(KIND=kreal), INTENT(IN)  :: pstrain !< plastic strain

     INTEGER :: htype
     htype = getHardenType( matl%mtype )
     SELECT CASE (htype)
     CASE(4, 5)   ! Prager
       calKinematicHarden = matl%variables(M_PLCONST3)
     CASE DEFAULT
       calKinematicHarden = 0.d0
     END SELECT
   END FUNCTION

   !> This function calcualtes state of kinematic hardening
   REAL(KIND=kreal) FUNCTION calCurrKinematic( matl, pstrain )
     TYPE( tMaterial ), INTENT(IN) :: matl    !< material property
     REAL(KIND=kreal), INTENT(IN)  :: pstrain !< plastic strain

     INTEGER :: htype
     htype = getHardenType( matl%mtype )
     SELECT CASE (htype)
     CASE(4, 5)   ! Prager
       calCurrKinematic = matl%variables(M_PLCONST3)*pstrain
     CASE DEFAULT
       calCurrKinematic = 0.d0
     END SELECT
   END FUNCTION

   !> This function calcualtes current yield stress
   REAL(KIND=kreal) FUNCTION calCurrYield( matl, pstrain, temp )
     TYPE( tMaterial ), INTENT(IN) :: matl    !< material property
     REAL(KIND=kreal), INTENT(IN)  :: pstrain !< plastic strain
     REAL(kind=kreal), INTENT(IN), optional :: temp  !< temperature

     INTEGER :: i, nc, htype
     REAL(KIND=kreal) :: s0, s1,s2, ina(2), outa(1)
	 logical :: ierr
     calCurrYield = -1.d0
     htype = getHardenType( matl%mtype )

     SELECT CASE (htype)
     CASE (0, 5)  ! Linear hardening, Linear+Parger hardening
       calCurrYield = matl%variables(M_PLCONST1)+matl%variables(M_PLCONST2)*pstrain
     CASE (1)  ! Multilinear approximation
       if( present(temp) ) then
         ina(1) = temp;  ina(2)=pstrain
         CALL fetch_TableData(MC_YIELD, matl%dict, outa, ierr, ina)
       else
         ina(1) = pstrain
         CALL fetch_TableData(MC_YIELD, matl%dict, outa, ierr, ina(1:1))
       endif
       if( ierr ) stop "Fail to get yield stress!"
       calCurrYield = outa(1)
     CASE (2)  ! Swift
       s0= matl%variables(M_PLCONST1)
       s1= matl%variables(M_PLCONST2)
       s2= matl%variables(M_PLCONST3)
       calCurrYield = s1*( s0+pstrain )**s2
     CASE (3)  ! Ramberg-Osgood
       s0= matl%variables(M_PLCONST1)
       s1= matl%variables(M_PLCONST2)
       s2= matl%variables(M_PLCONST3)
       if( pstrain<=s0 ) then
         calCurrYield = s1
       else
         calCurrYield = s1*( pstrain/s0 )**(1.d0/s2)
       endif
     CASE (4)  ! Parger hardening
       calCurrYield = matl%variables(M_PLCONST1)
     END SELECT
   END FUNCTION

   !> This function calcualtes yield state
   REAL(KIND=kreal) FUNCTION calYieldFunc( matl, stress, extval, temp )
     TYPE( tMaterial ), INTENT(IN) :: matl        !< material property
     REAL(KIND=kreal), INTENT(IN)  :: stress(6)   !< stress
     REAL(KIND=kreal), INTENT(IN)  :: extval(:)   !< plastic strain, back stress
     REAL(kind=kreal), INTENT(IN), optional :: temp  !< temperature

     INTEGER :: ytype
     logical :: kinematic
     REAL(kind=kreal) :: eqvs, sita, eta, fai, J1,J2,J3, f, devia(6)
     REAL(KIND=kreal) :: pstrain, back(6)

     kinematic = isKinematicHarden( matl%mtype )
     if( kinematic ) back(1:6) = extval(1:6)

     pstrain = extval(1)
     ytype = getYieldFunction( matl%mtype )
     J1 = (stress(1)+stress(2)+stress(3))
     devia(1:3) = stress(1:3)-J1/3.d0
     devia(4:6) = stress(4:6)
     if( kinematic ) devia = devia-back

     J2 = 0.5d0* dot_product( devia(1:3), devia(1:3) ) +  &
          dot_product( devia(4:6), devia(4:6) )
     if( present(temp) ) then
       eqvs = calCurrYield( matl, pstrain, temp )
     else
       eqvs = calCurrYield( matl, pstrain )
     endif

     SELECT CASE (yType)
     CASE (0)       ! Mises or. Isotropic
        f = dsqrt( 3.d0*J2 ) - eqvs
     CASE (1)      ! Mohr-Coulomb
       fai = matl%variables(M_PLCONST3)
       J3 = devia(1)*devia(2)*devia(3)                    &
           +2.d0* devia(4)*devia(5)*devia(6)                     &
           -devia(2)*devia(6)*devia(6)                           &
           -devia(3)*devia(4)*devia(4)                           &
           -devia(1)*devia(5)*devia(5)
       sita = -3.d0*dsqrt(3.d0)*J3/( 2.d0*(J2**1.5d0) )
       if( dabs( dabs(sita)-1.d0 ) <1.d-8 ) sita=SIGN(1.d0, sita)
       if( dabs(sita) >1.d0 ) stop "Math Error in Mohr-Coulomb calculation"
       sita = asin( sita )/3.d0
       f = (cos(sita)-sin(sita)*sin(fai)/dsqrt(3.d0))*dsqrt(J2)  &
             +J1*sin(fai)/3.d0 - eqvs*cos(fai)
     CASE (2)      ! Drucker-Prager
       eta = matl%variables(M_PLCONST3)
       f = dsqrt(J2) + eta*J1 - eqvs*matl%variables(M_PLCONST4)
     END SELECT

     calYieldFunc = f
   END FUNCTION

   !> This subroutine does backward-Euler return calculation
   subroutine BackwardEuler( matl, stress, plstrain, istat, fstat, temp )
      use m_utilities, only : eigen3
      type( tMaterial ), intent(in)    :: matl        !< material properties
      real(kind=kreal), intent(inout)  :: stress(6)   !< trial->real stress
      real(kind=kreal), intent(in)     :: plstrain    !< plastic strain till current substep
      integer, intent(inout)           :: istat       !< plastic state
      real(kind=kreal), intent(inout)  :: fstat(:)    !< plastic strain, back stress
      REAL(kind=kreal), INTENT(IN), optional :: temp  !< temperature

      real(kind=kreal), parameter :: tol =1.d-3
      integer, parameter          :: MAXITER = 5
      real(kind=kreal) :: dlambda, f, mat(3,3)
      integer :: i,j,ytype, maxp(1), minp(1), mm
      real(kind=kreal) :: youngs, poisson, pstrain, dum, ina(1), ee(2)
      real(kind=kreal) :: J1,J2,J3, H, KH, KK, dd, yd, G, K, devia(6)
      real(kind=kreal) :: prnstre(3), prnprj(3,3), tstre(3,3)
      real(kind=kreal) :: sita, fai, dep, trialprn(3)
      real(kind=kreal) :: a,b,siga,sigb,lamab(2),fab(2)
      real(kind=kreal) :: resi(2,2), invd(2,2), fstat_bak(6)
      logical          :: right, kinematic, ierr
      REAL(KIND=kreal) :: ftrial, betan, back(6)

      ytype = getYieldFunction( matl%mtype )
      if( ytype==3 ) THEN
       ! Gaku Hashimoto (2015/08/28) <
       ! call uBackwardEuler( matl, stress, istat, fstat )
       ! > Gaku Hashimoto (2015/08/28)
       return
      endif

      pstrain = plstrain
      fstat_bak = plstrain
      if( present(temp) ) then
        f = calYieldFunc( matl, stress, fstat_bak, temp )
      else
        f = calYieldFunc( matl, stress, fstat_bak )
      endif
      if( dabs(f)<tol ) then  ! yielded
        istat = 1
        return
      elseif( f<0.d0 ) then   ! not yielded or unloading
          istat =0
          return
      endif
      istat = 1           ! yielded
      KH = 0.d0; KK=0.d0; betan=0.d0; back(:)=0.d0

      kinematic = isKinematicHarden( matl%mtype )
      if( kinematic ) then
        back(1:6) = fstat(1:6)
        betan = calCurrKinematic( matl, pstrain )
      endif

      J1 = (stress(1)+stress(2)+stress(3))/3.d0
      devia(1:3) = stress(1:3)-J1
      devia(4:6) = stress(4:6)
      if( kinematic ) devia = devia-back
      yd = cal_equivalent_stress(matl, stress, fstat)

      if( present(temp) ) then
        ina(1) = temp
        CALL fetch_TableData(MC_ISOELASTIC, matl%dict, ee, ierr, ina)
      else
        call fetch_TableData(MC_ISOELASTIC, matl%dict, ee, ierr )
      endif
      if( ierr ) then
         stop " fail to fetch young's modulus in elastoplastic calculation"
      else
        youngs = ee(1)
        poisson = ee(2)
      endif
      if( youngs==0.d0 ) stop "YOUNG's ratio==0"
      G = youngs/ ( 2.d0*(1.d0+poisson) )
      K = youngs/ ( 3.d0*(1.d0-2.d0*poisson) )
      dlambda = 0.d0

      if( yType==0 ) then    ! Mises or. Isotropic
        do i=1,MAXITER
          if( present(temp) ) then
            H= calHardenCoeff( matl, pstrain+dlambda, temp )
          else
            H= calHardenCoeff( matl, pstrain+dlambda )
          endif
          if( kinematic ) then
            KH = calKinematicHarden( matl, pstrain+dlambda )
          endif
          dd= 3.d0*G+H+KH
          dlambda = dlambda+f/dd
          if( dlambda<0.d0 ) then
            dlambda = 0.d0
            istat=0; exit
          endif
          if( present(temp) ) then
            dum = calCurrYield( matl, pstrain+dlambda, temp )
          else
            dum = calCurrYield( matl, pstrain+dlambda )
          endif
          if( kinematic ) then
            KK = calCurrKinematic( matl, pstrain+dlambda )
          endif
          f = yd-3.d0*G*dlambda-dum -(KK-betan)
          if( dabs(f)<tol*tol ) exit
        enddo
        pstrain = pstrain+dlambda
        if( kinematic ) then
          KK = calCurrKinematic( matl, pstrain )
          fstat(1:6) = back(:)+(KK-betan)*devia(:)/yd
        endif
        devia(:) = (1.d0-3.d0*dlambda*G/yd)*devia(:)
        stress(1:3) = devia(1:3)+J1
        stress(4:6) = devia(4:6)
        stress(:)= stress(:)+back(:)
      elseif(yType==1) then    ! Mohr-Coulomb
        fai = matl%variables(M_PLCONST3)

     !   do j=1,MAXITER
        J2 = 0.5d0* dot_product( devia(1:3), devia(1:3) ) +  &
          dot_product( devia(4:6), devia(4:6) )
        J3 = devia(1)*devia(2)*devia(3)                    &
           +2.d0* devia(4)*devia(5)*devia(6)                     &
           -devia(6)*devia(2)*devia(6)                           &
           -devia(4)*devia(4)*devia(3)                           &
           -devia(1)*devia(5)*devia(5)
        sita = -3.d0*dsqrt(3.d0)*J3/( 2.d0*(J2**1.5d0) )
        if( dabs( dabs(sita)-1.d0 ) <1.d-8 ) sita=SIGN(1.d0, sita)
        if( dabs(sita) >1.d0 ) stop "Math Error in Mohr-Coulomb calculation"
        sita = asin( sita )/3.d0
        do mm=1,6
          if( dabs(stress(mm))<1.d-10 ) stress(mm)=0.d0
        enddo
        call eigen3( stress, prnstre, prnprj )
        trialprn = prnstre
        maxp = MAXLOC( prnstre )
        minp = MINLOC( prnstre )
        mm = 1
        if( maxp(1)==1 .or. minp(1)==1 ) mm =2
        if( maxp(1)==2 .or. minp(1)==2 ) mm =3
        do i=1,MAXITER
          if( present(temp) ) then
            H= calHardenCoeff( matl, pstrain, temp )
          else
            H= calHardenCoeff( matl, pstrain )
          endif
          dd= 4.d0*G*( 1.d0+sin(fai)*sin(sita)/3.d0 )+4.d0*K         &
             *sin(fai)*sin(sita)+4.d0*H*cos(fai)*cos(fai)
          dlambda = dlambda+f/dd
          if( 2.d0*dlambda*cos(fai)<0.d0 ) then
            if( cos(fai)==0.d0 ) STOP "Math error in return mapping"
            dlambda = 0.d0
            istat=0; exit
          endif
          dum = pstrain + 2.d0*dlambda*cos(fai)
          if( present(temp) ) then
            yd = calCurrYield( matl, dum, temp )
          else
            yd = calCurrYield( matl, dum )
          endif
          f = prnstre(maxp(1))-prnstre(minp(1))+                     &
            (prnstre(maxp(1))+prnstre(minp(1)))*sin(fai)-            &
            (4.d0*G*(1.d0+sin(fai)*sin(sita)/3.d0)+4.d0*K*sin(fai)   &
            *sin(sita))*dlambda-2.d0*yd*cos(fai)
          if( dabs(f)<tol ) exit
        enddo
        pstrain = pstrain + 2.d0*dlambda*cos(fai)
        prnstre(maxp(1)) = prnstre(maxp(1))-(2.d0*G*(1.d0+sin(fai)/3.d0)  &
                 + 2.d0*K*sin(fai) )*dlambda
        prnstre(minp(1)) = prnstre(minp(1))+(2.d0*G*(1.d0-sin(fai)/3.d0)  &
                 - 2.d0*K*sin(fai) )*dlambda
        prnstre(mm) = prnstre(mm)+(4.d0*G/3.d0-2.d0*K)*sin(fai)*dlambda

        tstre(:,:) = 0.d0
        tstre(1,1)= prnstre(1); tstre(2,2)=prnstre(2); tstre(3,3)=prnstre(3)
        mat= matmul( prnprj, tstre )
		mat= matmul( mat, transpose(prnprj) )
        stress(1) = mat(1,1)
        stress(2) = mat(2,2)
        stress(3) = mat(3,3)
        stress(4) = mat(1,2)
        stress(5) = mat(2,3)
        stress(6) = mat(3,1)
      elseif(yType==2) then    ! Drucker-Prager
        fai = matl%variables(M_PLCONST3)
        dum = matl%variables(M_PLCONST4)
        do i=1,MAXITER
          if( present(temp) ) then
            H= calHardenCoeff( matl, pstrain, temp )
          else
            H= calHardenCoeff( matl, pstrain )
          endif
          dd= G+K*fai*fai+H*dum*dum
          dlambda = dlambda+f/dd
          if( dum*dlambda<0.d0 ) then
            if( dum==0.d0 ) STOP "Math error in return mapping"
            dlambda = 0.d0
            istat=0; exit
          endif
          if( present(temp) ) then
            f = calCurrYield( matl, pstrain+dum*dlambda, temp  )
          else
            f = calCurrYield( matl, pstrain+dum*dlambda  )
          endif
          f = yd-G*dlambda+fai*(J1-K*fai*dlambda)- dum*f
          if( dabs(f)<tol*tol ) exit
        enddo
        pstrain = pstrain+dum*dlambda
        devia(:) = (1.d0-G*dlambda/yd)*devia(:)
        J1 = J1-K*fai*dlambda
        stress(1:3) = devia(1:3)+J1
        stress(4:6) = devia(4:6)
      end if

      fstat(1) = pstrain
   end subroutine BackwardEuler
      
      
      ! Gaku Hashimoto (2015/07/20) <
!####################################################################
      SUBROUTINE BackwardEuler_modified(icel, lx, matl, stress, strain, plstrain, delta_gamma, apex, istat, fstat, temp)
!####################################################################
      
      USE m_utilities, ONLY : eigen3
      
!--------------------------------------------------------------------
      
      INTEGER, INTENT(IN)              :: icel
      INTEGER, INTENT(IN)              :: lx
      TYPE( tMaterial ), INTENT(IN)    :: matl        !< material properties
      REAL(kind=kreal), INTENT(INOUT)  :: stress(6)   !< trial->real stress
      REAL(kind=kreal), INTENT(INOUT)  :: strain(6)   !< elastic strain
      REAL(kind=kreal), INTENT(IN)     :: plstrain    !< plastic strain till current substep
      REAL(kind=kreal), INTENT(INOUT)  :: delta_gamma !< incremental plastic multiplier
      LOGICAL, INTENT(INOUT)           :: apex        !< Apex return
      INTEGER, INTENT(INOUT)           :: istat       !< plastic state
      REAL(kind=kreal), INTENT(INOUT)  :: fstat(:)    !< plastic strain, back stress
      REAL(kind=kreal), INTENT(IN), OPTIONAL :: temp  !< temperature
      
!--------------------------------------------------------------------
      
      INTEGER :: MAXITER
      INTEGER :: iter
      INTEGER :: iter1, iter2
      INTEGER :: ytype
      INTEGER :: maxp(1), minp(1)
      INTEGER :: mm
      
      REAL(kind=kreal) :: tol
      REAL(kind=kreal) :: ina(1), outa(2)
      REAL(kind=kreal) :: e, nu
      REAL(KIND=kreal) :: kappa, g
      REAL(kind=kreal) :: pstrain
      REAL(KIND=kreal) :: pstrain_trial
      REAL(KIND=kreal) :: fstat_trial(6)
      REAL(KIND=kreal) :: tr_strain_e
      REAL(KIND=kreal) :: tr_strain_e_trial
      REAL(KIND=kreal) :: strain_e_trial_dev(6)
      REAL(kind=kreal) :: p_trial
      REAL(kind=kreal) :: p
      REAL(kind=kreal) :: stress_trial_dev(6)
      REAL(kind=kreal) :: stress_trial(6)
      REAL(kind=kreal) :: tr_delta_strain_p
      REAL(kind=kreal) :: dtr_delta_strain_p
      REAL(kind=kreal) :: q_trial
      REAL(kind=kreal) :: harden, khard
      REAL(kind=kreal) :: KK
      REAL(kind=kreal) :: coef
      REAL(kind=kreal) :: f
      REAL(kind=kreal) :: sigma_y, c
      REAL(kind=kreal) :: j1, j2, j3
      REAL(KIND=kreal) :: prnstre(3), prnprj(3,3), tstre(3,3)
      REAL(KIND=kreal) :: sita, fai, dep, trialprn(3)
      REAL(KIND=kreal) :: xi, eta
      REAL(KIND=kreal) :: a, b, siga, sigb, lamab(2), fab(2)
      REAL(KIND=kreal) :: resi(2,2), invd(2,2)
      REAL(KIND=kreal) :: ftrial, betan, back(6)
      REAL(KIND=kreal) :: mat(3, 3)
      REAL(KIND=kreal) :: ddelta_gamma
      
      LOGICAL :: right, kinematic, ierr
      
!--------------------------------------------------------------------
      
      MAXITER = 50
      tol = 1.0D-6
      
!--------------------------------------------------------------------
      
      IF( PRESENT( temp ) ) THEN
       
       ina(1) = temp
       
       CALL fetch_TableData(MC_ISOELASTIC, matl%dict, outa, ierr, ina)
       
       IF( ierr ) THEN
        
        e  = matl%variables(M_YOUNGS)
        nu = matl%variables(M_POISSON)
        
       ELSE
        
        e  = outa(1)
        nu = outa(2)
        
       END IF
       
      ELSE
       
       CALL fetch_TableData(MC_ISOELASTIC, matl%dict, outa, ierr )
       
       IF( ierr ) THEN
        
        e  = matl%variables(M_YOUNGS)
        nu = matl%variables(M_POISSON)
        
       ELSE
        
        e  = outa(1)
        nu = outa(2)
        
       END IF
       
      END IF
      
!--------------------------------------------------------------------
      
      kappa = e/( 3.0D0*(1.0D0-2.0D0*nu) )
      g = e/( 2.0D0*(1.0D0+nu) )
      
!--------------------------------------------------------------------
      
      pstrain_trial = plstrain
      
      fstat_trial(1) = pstrain_trial
      fstat_trial(2) = pstrain_trial
      fstat_trial(3) = pstrain_trial
      fstat_trial(4) = pstrain_trial
      fstat_trial(5) = pstrain_trial
      fstat_trial(6) = pstrain_trial
      
!--------------------------------------------------------------------
      
      kinematic = isKinematicHarden( matl%mtype )
      
      IF( kinematic ) THEN
       
       back(1) = fstat(1)
       back(2) = fstat(2)
       back(3) = fstat(3)
       back(4) = fstat(4)
       back(5) = fstat(5)
       back(6) = fstat(6)
       
       betan = calCurrKinematic(matl, pstrain_trial)
       
      END IF

      !--------------------------------------------------------------
      
      ! Elastic trial volumetric strain
      tr_strain_e_trial = strain(1)+strain(2)+strain(3)
      
      ! Elastic trial deviatoric strain
      strain_e_trial_dev(1) = strain(1)-tr_strain_e_trial/3.0D0
      strain_e_trial_dev(2) = strain(2)-tr_strain_e_trial/3.0D0
      strain_e_trial_dev(3) = strain(3)-tr_strain_e_trial/3.0D0
      strain_e_trial_dev(4) = strain(4)
      strain_e_trial_dev(5) = strain(5)
      strain_e_trial_dev(6) = strain(6)
      
!--------------------------------------------------------------------
      
      ! Elastic trial pressure
      p_trial = kappa*tr_strain_e_trial
      
      ! Elastic trial deviatoric stress
      stress_trial_dev(1) = 2.0D0*g*strain_e_trial_dev(1)
      stress_trial_dev(2) = 2.0D0*g*strain_e_trial_dev(2)
      stress_trial_dev(3) = 2.0D0*g*strain_e_trial_dev(3)
      stress_trial_dev(4) =       g*strain_e_trial_dev(4)
      stress_trial_dev(5) =       g*strain_e_trial_dev(5)
      stress_trial_dev(6) =       g*strain_e_trial_dev(6)
      
      IF( kinematic ) THEN
       
       stress_trial_dev(1) = stress_trial_dev(1)-back(1)
       stress_trial_dev(2) = stress_trial_dev(2)-back(2)
       stress_trial_dev(3) = stress_trial_dev(3)-back(3)
       stress_trial_dev(4) = stress_trial_dev(4)-back(4)
       stress_trial_dev(5) = stress_trial_dev(5)-back(5)
       stress_trial_dev(6) = stress_trial_dev(6)-back(6)
       
      END IF
      
      ! Elastic trial stress
      stress_trial(1) = p_trial+stress_trial_dev(1)
      stress_trial(2) = p_trial+stress_trial_dev(2)
      stress_trial(3) = p_trial+stress_trial_dev(3)
      stress_trial(4) =         stress_trial_dev(4)
      stress_trial(5) =         stress_trial_dev(5)
      stress_trial(6) =         stress_trial_dev(6)
      
!--------------------------------------------------------------------
      
      ! Uniaxial yield stress
      q_trial = cal_equivalent_stress(matl, stress_trial, fstat)
      
!--------------------------------------------------------------------
      
      ! Yield function
      IF( PRESENT(temp) ) THEN
       
       f = calYieldFunc(matl, stress_trial, fstat_trial, temp)
       
      ELSE
       
       f = calYieldFunc(matl, stress_trial, fstat_trial)
       
      END IF
      
!--------------------------------------------------------------------
      
      IF( DABS( f ) .LT. tol ) THEN
       
       ! Yielded
       istat = 1
       
       RETURN
       
      ELSE IF( f .LE. -tol ) THEN
       
       ! Not yielded or unloading
       istat =0
       
       RETURN
       
      END IF
      
      ! Yielded
      istat = 1
      
!--------------------------------------------------------------------
      
      ytype = getYieldFunction(matl%mtype)
      
!--------------------------------------------------------------------
      
      ! von Mises
      IF( yType .EQ. 0 ) THEN
       
       !--------------------------------------------------------
       
       delta_gamma = 0.0D0
       
       pstrain = pstrain_trial
       
       !--------------------------------------------------------
       
       DO iter = 1, MAXITER
        
        !--------------------------------------------------
        
        IF( PRESENT(temp) ) THEN
         
         harden = calHardenCoeff(matl, pstrain, temp)
         
        ELSE
         
         harden = calHardenCoeff(matl, pstrain)
         
        END IF
        
        !--------------------------------------------------
        
        khard = 0.0D0
        
        IF( kinematic ) THEN
         
         khard = calKinematicHarden(matl, pstrain)
         
        END IF
        
        !--------------------------------------------------
        
        ddelta_gamma = -f/( -3.0D0*g-( harden+khard ) )
        delta_gamma = delta_gamma+ddelta_gamma
        
        pstrain = pstrain_trial+delta_gamma
        
        !--------------------------------------------------
        
        IF( delta_gamma/sigma_y .LT. 0.0D0 ) THEN
         
         ! Not yielded or unloading
         istat = 0
         
         delta_gamma = 0.0D0
         
         pstrain = pstrain_trial
         
         EXIT
         
        END IF
        
        !--------------------------------------------------
        
        IF( PRESENT(temp) ) THEN
         
         sigma_y = calCurrYield(matl, pstrain, temp)
         
        ELSE
         
         sigma_y = calCurrYield(matl, pstrain)
         
        END IF
        
        !--------------------------------------------------
        
        IF( kinematic ) THEN
         
         KK = calCurrKinematic(matl, pstrain)
         
         ! Yield function
         f = q_trial-3.d0*g*delta_gamma-sigma_y-(KK-betan)
         
        ELSE
         
         ! Yield function
         f = q_trial-3.0d0*g*delta_gamma-sigma_y
         
        END IF
        
        !--------------------------------------------------
        
        IF( DABS( f ) < tol ) THEN
         
         IF( kinematic ) THEN
          
          KK = calCurrKinematic( matl, pstrain )
          
          fstat(1) = back(1)+( KK-betan )*stress_trial_dev(1)/q_trial
          fstat(2) = back(2)+( KK-betan )*stress_trial_dev(2)/q_trial
          fstat(3) = back(3)+( KK-betan )*stress_trial_dev(3)/q_trial
          fstat(4) = back(4)+( KK-betan )*stress_trial_dev(4)/q_trial
          fstat(5) = back(5)+( KK-betan )*stress_trial_dev(5)/q_trial
          fstat(6) = back(6)+( KK-betan )*stress_trial_dev(6)/q_trial
          
         END IF
         
         coef = 1.0D0-3.0D0*g*delta_gamma/q_trial
         
         stress(1) = p_trial+coef*stress_trial_dev(1)
         stress(2) = p_trial+coef*stress_trial_dev(2)
         stress(3) = p_trial+coef*stress_trial_dev(3)
         stress(4) =         coef*stress_trial_dev(4)
         stress(5) =         coef*stress_trial_dev(5)
         stress(6) =         coef*stress_trial_dev(6)
         
         IF( kinematic ) THEN
          
          stress(1)= stress(1)+back(1)
          stress(2)= stress(2)+back(2)
          stress(3)= stress(3)+back(3)
          stress(4)= stress(4)+back(4)
          stress(5)= stress(5)+back(5)
          stress(6)= stress(6)+back(6)
          
         END IF
         
         strain(1) = tr_strain_e_trial/3.0D0+coef*strain_e_trial_dev(1)
         strain(2) = tr_strain_e_trial/3.0D0+coef*strain_e_trial_dev(2)
         strain(3) = tr_strain_e_trial/3.0D0+coef*strain_e_trial_dev(3)
         strain(4) =                         coef*strain_e_trial_dev(4)
         strain(5) =                         coef*strain_e_trial_dev(5)
         strain(6) =                         coef*strain_e_trial_dev(6)
         
         RETURN
         
        END IF
        
        !--------------------------------------------------
        
       END DO
       
       !--------------------------------------------------------
       
      ! Drucker-Prager
      ELSE IF( yType .EQ. 2 ) THEN
       
       !--------------------------------------------------------
       
       eta = 3.0D0*matl%variables(M_PLCONST3)
       xi  = matl%variables(M_PLCONST4)
       
       !--------------------------------------------------------
       
       apex = .FALSE.
       
       delta_gamma = 0.0D0
       
       pstrain = pstrain_trial
       
       !--------------------------------------------------------
       
       ! Smooth cone wall return
       DO iter1 = 1, MAXITER
        
        !--------------------------------------------------
        
        IF( PRESENT( temp ) ) THEN
         
         harden = calHardenCoeff( matl, pstrain, temp )
         
        ELSE
         
         harden = calHardenCoeff( matl, pstrain )
         
        END IF
        
        !--------------------------------------------------
        
        ddelta_gamma = -f/( -g-kappa*eta*eta-xi*xi*harden )
        
        delta_gamma = delta_gamma+ddelta_gamma
        
        pstrain = pstrain_trial+xi*delta_gamma
        
        p = p_trial-delta_gamma*kappa*eta
        
        ! IF( xi*delta_gamma .LT. 0.0D0 ) THEN
         
        !  ! Not yielded or unloading
        !  istat = 0
         
        !  delta_gamma = 0.0D0
         
        !  pstrain = pstrain_trial
         
        !  EXIT
         
        ! END IF
        
        IF( PRESENT( temp ) ) THEN
         
         sigma_y = calCurrYield(matl, pstrain, temp)
         
        ELSE
         
         sigma_y = calCurrYield(matl, pstrain)
         
        END IF
        
        f = q_trial-g*delta_gamma+eta*p-xi*sigma_y
        
        !--------------------------------------------------
        
        IF( DABS( f ) .LT. tol ) THEN
         
         !--------------------------------------------
         
         IF( q_trial-g*delta_gamma .GE. 0.0D0 ) THEN
          
          !--------------------------------------
          
          IF( q_trial .EQ. 0.0D0 ) THEN
           
           coef = 0.0D0
           
          ELSE
           
           coef = 1.0D0-delta_gamma*g/q_trial
           
          END IF
          
          stress(1) = p+coef*stress_trial_dev(1)
          stress(2) = p+coef*stress_trial_dev(2)
          stress(3) = p+coef*stress_trial_dev(3)
          stress(4) =   coef*stress_trial_dev(4)
          stress(5) =   coef*stress_trial_dev(5)
          stress(6) =   coef*stress_trial_dev(6)
          
          tr_strain_e = p/kappa
          
          strain(1) = tr_strain_e/3.0D0+coef*strain_e_trial_dev(1)
          strain(2) = tr_strain_e/3.0D0+coef*strain_e_trial_dev(2)
          strain(3) = tr_strain_e/3.0D0+coef*strain_e_trial_dev(3)
          strain(4) =                   coef*strain_e_trial_dev(4)
          strain(5) =                   coef*strain_e_trial_dev(5)
          strain(6) =                   coef*strain_e_trial_dev(6)
          
          fstat(1) = pstrain
          
          RETURN
          
          !--------------------------------------
          
         ELSE
          
          !--------------------------------------
          
          apex = .TRUE.
          
          !--------------------------------------
          
          tr_delta_strain_p = 0.0D0
          
          pstrain = pstrain_trial
          
          p = p_trial
          
          !--------------------------------------
          
          IF( PRESENT( temp ) ) THEN
           
           c = calCurrYield(matl, pstrain, temp)
           
          ELSE
           
           c = calCurrYield(matl, pstrain)
           
          END IF
          
          !--------------------------------------
          
          ! Apex return
          DO iter2 = 1, MAXITER
           
           IF( PRESENT( temp ) ) THEN
            
            harden = calHardenCoeff( matl, pstrain, temp )
            
           ELSE
            
            harden = calHardenCoeff( matl, pstrain )
            
           END IF
           
           dtr_delta_strain_p = -( ( xi/eta )*c-p )/( ( xi/eta )*( xi/eta )*harden+kappa )
           
           tr_delta_strain_p = tr_delta_strain_p+dtr_delta_strain_p
           
           pstrain = pstrain_trial+( xi/eta )*tr_delta_strain_p
           
           IF( PRESENT( temp ) ) THEN
            
            c = calCurrYield(matl, pstrain, temp)
            
           ELSE
            
            c = calCurrYield(matl, pstrain)
            
           END IF
           
           p = p_trial-kappa*tr_delta_strain_p
           
           IF( DABS( ( ( xi/eta )*c-p )/c ) .LE. tol ) THEN
            
            delta_gamma = tr_delta_strain_p/eta
            
            stress(1) = p
            stress(2) = p
            stress(3) = p
            stress(4) = 0.0D0
            stress(5) = 0.0D0
            stress(6) = 0.0D0
            
            tr_strain_e = p/kappa
            
            strain(1) = tr_strain_e/3.0D0
            strain(2) = tr_strain_e/3.0D0
            strain(3) = tr_strain_e/3.0D0
            strain(4) = 0.0D0
            strain(5) = 0.0D0
            strain(6) = 0.0D0
            
            fstat(1) = pstrain
            
            RETURN
            
           END IF
           
          END DO
          
          !--------------------------------------
          
         END IF
         
         !--------------------------------------------
         
        END IF
        
        !--------------------------------------------------
        
       END DO
       
      ! User subroutine
      ELSE IF( ytype .EQ. 3 ) THEN
       
       !--------------------------------------------------------
       
       ! Gaku Hashimoto (2015/08/28) <
       ! CALL uBackwardEuler(matl, stress, istat, fstat)
       ! > Gaku Hashimoto (2015/08/28)
       
       !--------------------------------------------------------
       
      END IF
      
!--------------------------------------------------------------------
      
      RETURN
      
!####################################################################
      END SUBROUTINE BackwardEuler_modified
!####################################################################
      ! > Gaku Hashimoto (2015/07/20)
      
      
   !> Clear elatoplastic state
   subroutine updateEPState( gauss )
      use mMechGauss
      type(tGaussStatus), intent(inout) :: gauss  ! status of curr gauss point
      gauss%plstrain= gauss%fstatus(1)
   end subroutine

end module m_ElastoPlastic
