!DEC$ FREEFORM

      !INCLUDE 'ABQinterface.f90' !comment this line out when used with MonolithFE2

      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
      RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,&
      TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS,&
      NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,&
      DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)
!
      USE ABQINTERFACE
!
      IMPLICIT NONE
!
      CHARACTER(LEN=80):: MATERL
      REAL(KIND=AbqRK)::  STRESS(NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS),&
                          DDSDDT(NTENS),DRPLDE(NTENS),STRAN(NTENS),DSTRAN(NTENS),&
                          TIME(2),PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3),DROT(3,3),&
                          DFGRD0(3,3),DFGRD1(3,3),DFGRIN(3,3),STRESS_trial(NTENS),&
                          STRESS_trial_dev(NTENS),SSE,SPD,SCD,RPL,DRPLDT,DTEMP,TEMP
!
      REAL(KIND=AbqRK):: NEWTON=15,TOLER=0.000001
      REAL(KIND=AbqRK):: EMOD,ENU,EBULK3,EG2,EG,ELAM
      REAL(KIND=AbqRK):: M,N,A,B,C,L,D,STRAN_CREEP_eq,STRESS_eq,K,&
                         STRESS_trial_eq,SQi,R,DRDSQ,DRDSQT,PNEWDT,CELENT,DTIME
      INTEGER(KIND=ABqIK):: NDI,K1,K2,NSHR,NTENS,NSTATV,NPROPS,NOEL,NPT,KSLAY,KSPT,KSTEP(:),KINC
!
! -----------------------------------------------------------
!     UMAT for isotropic elasticity and creep
!     cannot be be used for plane stress
! -----------------------------------------------------------
!     PROPS(1) - E
!     PROPS(2) - NU
!     PROPS(3) - M
!     PROPS(4) - N
!     PROPS(5) - A
!
!     creep law: DSTRAN_CREEP_eq=B*STRESS_EQ**C*STRAN_CREEP_eq**L*DTIME
!
!     euler-backward
! -----------------------------------------------------------
!
      PNEWDT=1.5_AbqRK

      IF (NDI.NE.3) THEN !output error
         CALL STDB_ABQERR(-3_AbqIK, 'Error: This UMAT may only be used for '//&
         'elements with three direct stress components', 0, 0.0_AbqRK, ' ')
      END IF

!     elastic properties

      EMOD=PROPS(1)
      ENU=PROPS(2)
      IF(ENU.GT.0.4999_AbqRK .AND. ENU .LT. 0.5001_AbqRK) ENU=0.499_AbqRK
      EBULK3=EMOD/(1.0_AbqRK-2.0_AbqRK*ENU)
      EG2=EMOD/(1.0_AbqRK+ENU)
      EG=EG2/2.0_AbqRK
      ELAM=(EBULK3-EG2)/3.0_AbqRK
      K=3.0_AbqRK*ENU/(1.0_AbqRK-2.0_AbqRK*ENU)/(1.0_AbqRK+ENU)*EMOD+2.0_AbqRK*EG

!     creep properties

      M=PROPS(3)
      N=PROPS(4)
      A=PROPS(5)
      B=((M+1)**M/(A)**N)**(1/(M+1))
      C=N/(M+1)
      L=M/(M+1)

!     elastic stiffness

      DO K1=1,NTENS
        DO K2=1,NTENS
           DDSDDE(K2,K1)=0.0_AbqRK
        END DO
      END DO

      DO K1=1,NDI
        DO K2=1,NDI
           DDSDDE(K2,K1)=ELAM
        END DO
        DDSDDE(K1,K1)=EG2+ELAM
      END DO
   
      DO K1=NDI+1,NTENS
        DDSDDE(K1,K1)=EG
      END DO
                                     
!     equivalent creep strain of last increment
         
      STRAN_CREEP_eq=STATEV(1)
!     
!     calculate trial stress
!
      STRESS_trial=STRESS+MATMUL(DDSDDE,DSTRAN)
!
!     calculate deviator of trialstress
!
      STRESS_trial_dev=STRESS_trial
!
      DO K1=1,NDI
          STRESS_trial_dev(K1)=STRESS_trial_dev(K1)-1.0_AbqRK/3.0_AbqRK*SUM(STRESS_trial(1:3))
      END DO

!     compute equivalent trial stress

      STRESS_trial_eq=0.0_AbqRK
!
      DO K1=1,NDI
         STRESS_trial_eq=STRESS_trial_eq+STRESS_trial_dev(K1)**2.0_AbqRK
      END DO
      DO K1=NDI+1,NTENS
         STRESS_trial_eq=STRESS_trial_eq+2.0_AbqRK*STRESS_trial_dev(K1)**2.0_AbqRK
      END DO
      STRESS_trial_eq=SQRT(3.0_AbqRK/2.0_AbqRK*STRESS_trial_eq)

!     initialize equivialent stress of current increment

      STRESS_eq=STRESS_trial_eq

!     newton loop to calculate stress_eq

      IF (STRESS_trial_eq.GT.0.0_AbqRK .AND. STRAN_CREEP_eq.GT.0.0_AbqRK) THEN

!       derivative of residual w.r.t. stress_eq

        DRDSQ=(-L/3.0_AbqRK/EG*(STRAN_CREEP_eq+(STRESS_trial_eq-STRESS_eq)/3/EG)**(L-1.0_AbqRK)+&
               (DTIME*3.0_AbqRK*EG*STRESS_eq**C*B+(STRESS_trial_eq-STRESS_eq)*DTIME*3.0_AbqRK&
               *EG*C*STRESS_eq**(C-1.0_AbqRK)*B)/(DTIME*3*EG*STRESS_eq**C*B)**2.0_AbqRK)
!
        DO K1=1,NEWTON
!
         SQi=STRESS_eq !only needed for termination criterion

!        residual of nonlinear equation  STRESS_eq<->STRESS_trial_eq

         R=((STRAN_CREEP_eq+(STRESS_trial_eq-STRESS_eq)/3.0_AbqRK/EG)**L-&
            (STRESS_trial_eq-STRESS_eq)/(DTIME*3*EG*STRESS_eq**C*B))

!        update stress_eq

         STRESS_eq=STRESS_eq-R/DRDSQ

!        derivative of residual w.r.t. STRESS_eq

         DRDSQ=(-L/3.0_AbqRK/EG*(STRAN_CREEP_eq+(STRESS_trial_eq-STRESS_eq)/3/EG)**(L-1.0_AbqRK)+&
               (DTIME*3.0_AbqRK*EG*STRESS_eq**C*B+(STRESS_trial_eq-STRESS_eq)*DTIME*3.0_AbqRK&
               *EG*C*STRESS_eq**(C-1.0_AbqRK)*B)/(DTIME*3*EG*STRESS_eq**C*B)**2.0_AbqRK)
 

         IF (ABS((SQi-STRESS_eq)/STRESS_trial_eq).LE.TOLER) THEN
         
           PNEWDT=1.5_AbqRK

!          derivative of r w.r.t. STRESS_trial_eq needed for tangent
         
           DRDSQT=L/3.0_AbqRK/EG*(STRAN_CREEP_eq+(STRESS_trial_eq-STRESS_eq)/3.0_AbqRK/EG)**&
                  (L-1.0_AbqRK)-1.0_AbqRK/(3.0_AbqRK*EG*STRESS_eq**C*DTIME*B)

           EXIT
         ELSE IF (K1==NEWTON) THEN
            
            PNEWDT=0.5_AbqRK

         END IF

        END DO

!     in first step, when STRAN_CREEP_eq_m==0.0, the solution STRESS_eq==STRESS_trial_eq .AND.
!     STRAN_CREEP_eq_m+1==0 exists mathematical, which needs to be excluded

      ELSE IF (STRESS_trial_eq.GT.0.0_AbqRK .AND. STRAN_CREEP_eq.EQ.0.0_AbqRK) THEN

         STRESS_eq=0.5*STRESS_trial_eq
     
!        derivative of Residual W.R.T. STRESS_eq

         DRDSQ=1.0_AbqRK-C/(L-1.0_AbqRK)*3.0_AbqRK*EG*STRESS_eq**((-C-1.0_AbqRK)/(L-1.0_AbqRK))&
               /(DTIME*B)**(1/(L-1))

!        derivative of Residual W.R.T. STRESS_trial_eq needed for tangent
           DRDSQT=-1.0

        DO K1=1,NEWTON

         SQi=STRESS_eq

!        residual of nonlinear equation  STRESS_eq<->STRESS_trial_eq

         R=STRESS_eq**(-C/(L-1.0_AbqRK))*3._AbqRK*EG/(DTIME*B)**(1.0_AbqRK/(L-1.0_AbqRK))&
           +STRESS_eq-STRESS_trial_eq

!        update STRESS_eq

         STRESS_eq=STRESS_eq-R/DRDSQ

!        DERIVATIVE OF R W.R.T. STRESS_eq

         DRDSQ=1.0_AbqRK-C/(L-1.0_AbqRK)*3.0_AbqRK*EG*STRESS_eq**((-C-1.0_AbqRK)/&
               (L-1.0_AbqRK))/(DTIME*B)**(1.0_AbqRK/(L-1.0_AbqRK))
!
         IF (ABS((SQi-STRESS_eq)/STRESS_trial_eq).LE.TOLER) THEN
            PNEWDT=1.5_AbqRK
            EXIT
         ELSE IF (K1==NEWTON) THEN
            PNEWDT=0.5_AbqRK
         END IF
 
       END DO

      END IF

!     calculate STRAN_CREEP_eq, STRESS and algorithmic consistent tangent

      IF (STRESS_eq.GT.0.0_AbqRK) THEN

!       STRAN_CREEP_eq

        STRAN_CREEP_eq=((STRESS_trial_eq-STRESS_eq)/(3.0_AbqRK*EG*STRESS_eq**C*DTIME*B))**(1.0_AbqRK/L)

!       STRESS

        STRESS=STRESS_trial_dev*(STRESS_eq/STRESS_trial_eq)
!       
        DO K1=1,NDI
           STRESS(K1)=STRESS(K1)+1.0_AbqRK/3.0_AbqRK*SUM(STRESS_trial(1:3))
        END DO

!       algorithmic consistent tangent DDSDDE

         DO K1=1,NTENS
           DO K2=1,NTENS
              DDSDDE(K2,K1)=0.0_AbqRK
           END DO
          END DO

         DO K1=1,NDI
           DO K2=1,NDI
               IF (K1.EQ.K2) THEN
                   DDSDDE(K1,K2)=K/3.0_AbqRK+4.0_AbqRK/3.0_AbqRK*EG*(STRESS_eq/STRESS_trial_eq)
               ELSE
                   DDSDDE(K1,K2)=K/3.0_AbqRK-2.0_AbqRK/3.0_AbqRK*EG*(STRESS_eq/STRESS_trial_eq)
               END IF
           END DO
         END DO

         DO K1=NDI+1,NTENS
            DDSDDE(K1,K1)=EG*(STRESS_eq/STRESS_trial_eq)
         END DO

         DO K1=1,NTENS
           DO K2=1,NTENS
             DDSDDE(K1,K2)=DDSDDE(K1,K2)+STRESS_trial_dev(K1)*STRESS_trial_dev(K2)*3.0_AbqRK*EG/&
                           STRESS_trial_eq**2.0_AbqRK*(-DRDSQT/DRDSQ-STRESS_eq/STRESS_trial_eq)
           END DO
         END DO

      END IF

!     save equivalent creep strain in STATEV(1)
 
      STATEV(1)=STRAN_CREEP_eq

      RETURN
      END SUBROUTINE UMAT
