!DEC$ FREEFORM

!=====================================================================
! - general empty Interface for UMAT to be used inside UelXCU
! - contains general plane stress algorithm
!=====================================================================

INCLUDE 'UMATUtilityModules.f'  !include file with USER Modules

MODULE UXMATModule
    USE ABQINTERFACE
    IMPLICIT NONE
    PRIVATE
    PUBLIC::UXMAT, CheckMaterialParameters, GetNSTATV
    CONTAINS
                                                                                                                                                                   
        INCLUDE 'UMAT.f'    !include the actual UMAT
        INCLUDE 'CheckMaterialParameters.f' !CheckMaterialParameters
        INCLUDE 'GetNSTATV.f' !GetNSTATV

        SUBROUTINE UXMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,DRPLDE,&
                         DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,&
                         CMNAME,NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,&
                         PNEWDT,CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
            ! Declare interface
            REAL(KIND=AbqRK),INTENT(IN):: TIME(:),PREDEF(:),STRAN(:),DSTRAN(:),&
                                          DPRED(:),PROPS(:),COORDS(:),DROT(:,:),&
                                          DFGRD0(:,:),DFGRD1(:,:),DTIME,TEMP,&
                                          DTEMP,CELENT
            INTEGER(KIND=AbqIK), INTENT(IN):: NTENS,NSTATV,NOEL,NPT,LAYER,KSPT,&
                                            KINC,NPROPS,NDI,NSHR,JSTEP(:)
            CHARACTER*80 CMNAME
            REAL(KIND=AbqRK)::STRESS(:),STATEV(:),DDSDDE(:,:),DDSDDT(:),DRPLDE(:),&
                            SSE,SPD,SCD,RPL,DRPLDT,PNEWDT,STRESS_PS(4),DSTRAN_PS(4),&
                            DDSDDE_PS(4,4),STRAN_PS(4),PNEWDT_PS,DDSDDT_PS(4),&
                            DRPLDE_PS(4),STATEV_PS(NSTATV-1),DFGRD0_PS(3,3),DFGRD1_PS(3,3)
            INTEGER(KIND=AbqIK),PARAMETER:: NDI_PS=3,NTENS_PS=4,NSHR_PS=1
            INTEGER(KIND=AbqIK):: I
            REAL(KIND=AbqRK),PARAMETER:: TOLER=0.00001
            INTEGER:: NEWTON=15
            
            IF (NTENS==3) THEN !plane stress case
                !STRAN of last increment
                STRAN_PS(1:2)=STRAN(1:2)
                STRAN_PS(3)=STATEV(1)
                STRAN_PS(4)=STRAN(3)
                !deformation gradient of last increment
                DFGRD0_PS=DFGRD0
                IF (JSTEP(3)==1) DFGRD0_PS(3,3)=EXP(STRAN_PS(3))
                !deformation gradient of current increment
                DFGRD1_PS=DFGRD1
                IF (JSTEP(3)==1) DFGRD1_PS(3,3)=EXP(STRAN_PS(3)) !initial guess
                !DSTRAN
                DSTRAN_PS(1:2)=DSTRAN(1:2)
                DSTRAN_PS(4)=DSTRAN(3)
                DSTRAN_PS(3)=0.0_AbqRK !initial guess

                DO I=1,NEWTON !iterate for DSTRAN_PS(3) to unti STRESS_PS(3)==0.0
                   !STATEV of last increment
                   STATEV_PS=STATEV(2:NSTATV)
                   !STRESS of last increment
                   STRESS_PS(1:2)=STRESS(1:2)
                   STRESS_PS(4)=STRESS(3)
                   STRESS_PS(3)=0.0_AbqRK
                   
                   PNEWDT_PS=1.0_AbqRK
                   ! Call UMAT in FORTRAN77 style  
                   CALL UMAT(STRESS_PS,STATEV_PS,DDSDDE_PS,SSE,SPD,SCD, RPL,DDSDDT_PS,&
                            DRPLDE_PS,DRPLDT, STRAN_PS,DSTRAN_PS,TIME,DTIME,TEMP,DTEMP,&
                            PREDEF,DPRED,CMNAME,NDI_PS,NSHR_PS,NTENS_PS,NSTATV-1,&
                            PROPS,NPROPS,COORDS,DROT,PNEWDT_PS, CELENT,DFGRD0_PS,&
                            DFGRD1_PS,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
                             
                    IF (ABS(STRESS_PS(3)/MAXVAL(ABS(STRESS_PS)))<=TOLER .OR.&
                        MAXVAL(ABS(STRESS_PS))==0.0_AbqRK) THEN
                    
                        IF (I<=NEWTON/2 .AND. PNEWDT_PS>=1.0_AbqRK) THEN
                            PNEWDT=1.5_AbqRK
                        ELSE
                            PNEWDT=PNEWDT_PS
                        END IF
                        
                        EXIT
                        
                    ELSE IF (I==NEWTON) THEN
                    
                        PNEWDT=0.5_AbqRK
                    
                    END IF
                    
                    !update strain in x3 direction
                    DSTRAN_PS(3)=DSTRAN_PS(3)-STRESS_PS(3)/DDSDDE_PS(3,3)
                    !update deformation gradient
                    IF (JSTEP(3)==1) DFGRD1_PS(3,3)=EXP(STRAN_PS(3)+DSTRAN_PS(3))
                    
                  END DO

                  !update STATEV
                  STATEV(1)=STATEV(1)+DSTRAN_PS(3) !updated stran in x3 direction
                  STATEV(2:NSTATV)=STATEV_PS
                  
                  !output reduced STRESS and reduced DDSDDE
                  
                  STRESS(1:2)=STRESS_PS(1:2)
                  STRESS(3)=STRESS_PS(4)
                  
                  DDSDDE(1,1)=DDSDDE_PS(1,1)-1/DDSDDE_PS(3,3)*DDSDDE_PS(1,3)**2
                  DDSDDE(1,2)=DDSDDE_PS(1,2)-1/DDSDDE_PS(3,3)*DDSDDE_PS(1,3)*DDSDDE_PS(2,3)
                  DDSDDE(2,1)=DDSDDE(1,2)
                  DDSDDE(1,3)=DDSDDE_PS(1,4)-1/DDSDDE_PS(3,3)*DDSDDE_PS(1,3)*DDSDDE_PS(4,3)
                  DDSDDE(3,1)=DDSDDE(1,3)
                  DDSDDE(2,2)=DDSDDE_PS(2,2)-1/DDSDDE_PS(3,3)*DDSDDE_PS(2,3)**2
                  DDSDDE(2,3)=DDSDDE_PS(2,4)-1/DDSDDE_PS(3,3)*DDSDDE_PS(2,3)*DDSDDE_PS(4,3)
                  DDSDDE(3,2)=DDSDDE(2,3)
                  DDSDDE(3,3)=DDSDDE_PS(4,4)-1/DDSDDE_PS(3,3)*DDSDDE_PS(4,3)**2
                     
            ELSE !3D, plane strain, axissymmetric
                ! Call UMAT in FORTRAN77 style
                CALL UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, RPL,DDSDDT,DRPLDE,DRPLDT,&
                        STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,NDI,&
                        NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,&
                        DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
            END IF
            
      END SUBROUTINE UXMAT
      
END MODULE UXMATModule
