!DEC$ FREEFORM

SUBROUTINE GETOUTDIR(OUTDIR,LENOUTDIR)
    !returns the current work directory
    
    IMPLICIT NONE
    
    !Name of directory in which UMAT_Driver is stared
    CHARACTER(LEN=256)::  OUTDIR
    !lenght of name of this directory
    INTEGER:: LENOUTDIR
    
    CALL GETCWD(OUTDIR)
    LENOUTDIR=LEN(TRIM(OUTDIR))
    
END SUBROUTINE GETOUTDIR

SUBROUTINE GETJOBNAME(JOBNAME,LENJOBNAME)
    !returns the jobname of the Training Inputfile
    
    USE type_trainingparameters
    
    IMPLICIT NONE
    
    !name of the current UMAT_Driver JOB
    CHARACTER(LEN=256)::  JOBNAME
    !lenght of that name
    INTEGER:: LENJOBNAME
    INTEGER:: error
    
    JOBNAME=training%Jobname
    LENJOBNAME=LEN_TRIM(training%Jobname)
    
END SUBROUTINE GETJOBNAME

SUBROUTINE XIT
    !Abaqus routine version of STOP
    
    STOP

END SUBROUTINE XIT

SUBROUTINE ROTSIG(S,R,SPRIME,LSTR,NDI,NSHR)
    ! Rotate a Stress or Strain Tensor in Voigt Notation
    
    USE ABQinterface
    
    IMPLICIT NONE
    
    INTEGER,INTENT(IN):: NDI,NSHR,LSTR
    INTEGER:: i,j
    REAL(KIND=AbqRK),DIMENSION(NDI+NSHR),INTENT(IN):: S
    REAL(KIND=AbqRK),DIMENSION(NDI+NSHR),INTENT(OUT):: SPRIME
    REAL(KIND=AbqRK),DIMENSION(3,3), INTENT(IN):: R
    REAL(KIND=AbqRK),DIMENSION(NDI+NSHR,NDI+NSHR):: Rotation
    INTEGER,DIMENSION(2,3):: Notation_shear=reshape([1,2,1,3,2,3],[2,3])
    
    FORALL(i=1:NDI,j=1:NDI) Rotation(i,j)=R(j,i)**2
    FORALL(i=1:NSHR,j=1:NDI) Rotation(NDI+i,j)=R(j,Notation_shear(1,i))*R(j,Notation_shear(2,i))
    FORALL(i=1:NDI,j=1:NSHR) Rotation(i,NDI+j)=R(Notation_shear(1,j),i)*R(Notation_shear(2,j),i)
    FORALL(i=1:NSHR,j=1:NSHR) Rotation(NDI+i,NDI+j)=R(Notation_shear(1,j),Notation_shear(1,i))*R(Notation_shear(2,j),Notation_shear(2,i))+&
                                                    R(Notation_shear(2,j),Notation_shear(1,i))*R(Notation_shear(1,j),Notation_shear(2,i))
    
    IF (LSTR==1) THEN !contains stresses
        Rotation(:NDI,NDI+1:)=Rotation(:NDI,NDI+1:)*2.0_AbqRK
        SPRIME=MATMUL(Rotation,S)
    ELSE IF (LSTR==2) THEN !contains strains
        Rotation(NDI+1:,:NDI)=Rotation(NDI+1:,:NDI)*2.0_AbqRK
        SPRIME=MATMUL(Rotation,S)
    ELSE
        STOP 'LSTR in ROTSIG shall be 1 or 2.'
    END IF
    
END SUBROUTINE ROTSIG

SUBROUTINE SPRINC(S,PS,LSTR,NDI,NSHR)
    ! calculate principal values of a stress or strain tensor
    
    USE ABQinterface
    
    IMPLICIT NONE
    
    REAL(KIND=AbqRK),DIMENSION(NDI+NSHR),INTENT(IN):: S
    REAL(KIND=AbqRK),DIMENSION(3,3):: S_TENS
    INTEGER,INTENT(IN):: NDI,NSHR,LSTR
    REAL(KIND=AbqRK),DIMENSION(3),INTENT(OUT):: PS
    INTEGER,DIMENSION(2,3):: Notation_shear=reshape([1,2,1,3,2,3],[2,3])
    INTEGER:: i,k,l,info,lwork
    REAL(KIND=AbqRK),DIMENSION(:),ALLOCATABLE:: work
    REAL(KIND=AbqRK):: shear_factor
    
    !get stress/strain in tensor form
    IF (LSTR==1) THEN !contains stresses 
        shear_factor=1.0_AbqRK
    ELSE IF (LSTR==2) THEN !contains strains
        shear_factor=0.5_AbqRK
    ELSE
        STOP 'LSTR in SPRINC shall be 1 or 2.'
    END IF
    
    !stress in full tensor notation
    S_TENS=0.0_AbqRK
    FORALL(i=1:NDI) S_TENS(i,i)=S(i)
    DO i=1,NSHR
        k=Notation_shear(1,i); l=Notation_shear(2,i)
        S_TENS(k,l)=S(NDI+i)*shear_factor; S_TENS(l,k)=S_TENS(k,l)
    END DO
    
    !compute principal values
    lwork=-1; ALLOCATE(work(1))
    CALL ssyev('N','U',3,S_TENS,3,PS,work,lwork,info)
    lwork=work(1); DEALLOCATE(work); ALLOCATE(work(lwork))
    CALL ssyev('N','U',3,S_TENS,3,PS,work,lwork,info)
    
END SUBROUTINE SPRINC

SUBROUTINE SPRIND(S,PS,AN,LSTR,NDI,NSHR)
    ! calculate principal values and directions of a stress or strain tensor
    
    USE ABQinterface
    
    IMPLICIT NONE
    
    REAL(KIND=AbqRK),DIMENSION(NDI+NSHR),INTENT(IN):: S
    REAL(KIND=AbqRK),DIMENSION(3,3),INTENT(OUT):: AN
    INTEGER,INTENT(IN):: NDI,NSHR,LSTR
    REAL(KIND=AbqRK),DIMENSION(3),INTENT(OUT):: PS
    INTEGER,DIMENSION(2,3):: Notation_shear=reshape([1,2,1,3,2,3],[2,3])
    INTEGER:: i,k,l,info,lwork
    REAL(KIND=AbqRK),DIMENSION(:),ALLOCATABLE:: work
    REAL(KIND=AbqRK):: shear_factor
    
    !get stress/strain in tensor form
    IF (LSTR==1) THEN !contains stresses 
        shear_factor=1.0_AbqRK
    ELSE IF (LSTR==2) THEN !contains strains
        shear_factor=0.5_AbqRK
    ELSE
        STOP 'LSTR in SPRIND shall be 1 or 2.'
    END IF
    
    !stress in full tensor notation
    AN=0.0_AbqRK
    FORALL(i=1:NDI) AN(i,i)=S(i)
    DO i=1,NSHR
        k=Notation_shear(1,i); l=Notation_shear(2,i)
        AN(k,l)=S(NDI+i)*shear_factor; AN(l,k)=AN(k,l)
    END DO
    
    !compute principal values
    lwork=-1; ALLOCATE(work(1))
    CALL ssyev('V','U',3,AN,3,PS,work,lwork,info)
    lwork=work(1); DEALLOCATE(work); ALLOCATE(work(lwork))
    CALL ssyev('V','U',3,AN,3,PS,work,lwork,info)
    AN=TRANSPOSE(AN)
    
END SUBROUTINE SPRIND

SUBROUTINE SINV(STRESS,SINV1,SINV2,NDI,NSHR)
    !compute teh stress invariants
    
    USE ABQinterface
    
    IMPLICIT NONE
    
    REAL(KIND=AbqRK),DIMENSION(NDI+NSHR),INTENT(IN):: STRESS
    INTEGER,INTENT(IN):: NDI,NSHR
    REAL(KIND=AbqRK),INTENT(OUT):: SINV1,SINV2
    INTEGER:: i
    
    !first invariant = 1/3 trace(STRESS)
    SINV1=SUM(STRESS(1:NDI))/3.0_AbqRK
    !second invariant = sqrt(3/2 STRESS_dev:STRESS_dev)
    SINV2=0.0_AbqRK
    DO i=1,NDI
        SINV2=SINV2+(STRESS(i)-SINV1)**2
    END DO
    DO i=NDI+1,NDI+NSHR
        SINV2=SINV2+2.0_AbqRK*(STRESS(i)**2)
    END DO
    SINV2=SQRT(1.5_AbqRK*SINV2)
    
END SUBROUTINE SINV

SUBROUTINE STDB_ABQERR(LOP,STRING,INTV,REALV,CHARV)
    !outputs user warnings; currently only output the information to the file log.msg
    
    USE ABQinterface
    
    IMPLICIT NONE
    
    INTEGER,INTENT(IN):: LOP
    CHARACTER(LEN=*),INTENT(IN):: STRING
    INTEGER(KIND=AbqIK),DIMENSION(*),INTENT(IN):: INTV
    REAL(KIND=AbqRK),DIMENSION(*),INTENT(IN):: REALV
    CHARACTER(LEN=256),DIMENSION(*),INTENT(IN):: CHARV
    CHARACTER(LEN=256),DIMENSION(100):: result_string
    CHARACTER(LEN=10000):: final_result_string
    INTEGER:: i,i_last,n,n_I,n_R,n_S
        
    i_last=1
    n=0;n_I=1;n_R=1;n_S=1
    DO i=1,LEN_TRIM(STRING)
        IF (n==size(result_string)) EXIT
        IF (STRING(i:i)=='%') THEN
            result_string(n+1)=STRING(i_last:i-1)
            i_last=i+2; n=n+1
            IF (STRING(i+1:i+1)=='I') THEN
                WRITE(result_string(n+1),*) INTV(n_I)
                n_I=n_I+1;n=n+1
            ELSE IF (STRING(i+1:i+1)=='R') THEN
                WRITE(result_string(n+1),*) REALV(n_R)
                n_R=n_R+1;n=n+1
            ELSE IF (STRING(i+1:i+1)=='S') THEN
                result_string(n+1)=CHARV(n_S)
                n_S=n_S+1;n=n+1
            ELSE
                n=1
                final_result_string='Error in specifying string in STDB_ABQERR.'
                EXIT
            END IF
        ELSE IF(i==LEN_TRIM(STRING)) THEN
            result_string(n+1)=STRING(i_last:i)
            n=n+1
        END IF
    END DO
    final_result_string=''
    DO i=1,n
        final_result_string=TRIM(final_result_string)//' '//TRIM(ADJUSTL(result_string(i)))
    END DO
    
    IF (LOP==-3) THEN
        WRITE(7,'(a)') 'error found: '//TRIM(final_result_string)
        STOP 'UMAT_Driver exited with errors, look in the .dat or .msg file for possible reasons'
    ELSE
        WRITE(7,'(a)') TRIM(final_result_string)
    END IF
    
END SUBROUTINE STDB_ABQERR

FUNCTION SMAIntArrayCreate(ident,Array_Size,Initval)

    USE SMAArrays
    USE mkl_service
    
    IMPLICIT NONE
    
    INTEGER(KIND=8):: SMAIntArrayCreate
    INTEGER(KIND=8):: alloc_size
    INTEGER(KIND=4),DIMENSION(*)::x
    POINTER(pt_x,x)
    INTEGER(KIND=4):: Array_Size
    INTEGER(KIND=4),INTENT(IN) :: ident
    INTEGER(KIND=4),INTENT(IN) :: INITVAL
    
    alloc_size=KIND(INITVAL)*Array_Size
    pt_x=mkl_malloc(alloc_size,KIND(INITVAL))
    x(1:Array_Size)=INITVAL
    PointerAdressesInt(ident)=pt_x
    SMAIntArrayCreate=pt_x
    
END FUNCTION SMAIntArrayCreate

FUNCTION SMAIntArrayAccess(ident)

    USE SMAArrays
    
    IMPLICIT NONE
    
    INTEGER(KIND=8):: SMAIntArrayAccess
    INTEGER(KIND=4),INTENT(IN) :: ident
    
    SMAIntArrayAccess=PointerAdressesInt(ident)
    
END FUNCTION SMAIntArrayAccess

SUBROUTINE SMAIntArrayDelete(ident)

    USE SMAArrays
    
    IMPLICIT NONE
    
    INTEGER(KIND=4),INTENT(IN) :: ident
    
    CALL mkl_free(PointerAdressesInt(ident))
    
END SUBROUTINE SMAIntArrayDelete

FUNCTION SMAFloatArrayCreateSP(ident,Array_Size,Initval)

    USE SMAArrays
    USE mkl_service
    
    IMPLICIT NONE
    
    INTEGER(KIND=8):: SMAFloatArrayCreateSP
    INTEGER(KIND=8):: alloc_size
    REAL(KIND=4),DIMENSION(*)::x
    POINTER(pt_x,x)
    INTEGER(KIND=4):: Array_Size
    INTEGER(KIND=4),INTENT(IN) :: ident
    REAL(KIND=4),INTENT(IN) :: INITVAL
    
    alloc_size=KIND(INITVAL)*Array_Size
    pt_x=mkl_malloc(alloc_size,KIND(INITVAL))
    x(1:Array_Size)=INITVAL
    PointerAdressesReal(ident)=pt_x
    SMAFloatArrayCreateSP=pt_x
    
END FUNCTION SMAFloatArrayCreateSP

FUNCTION SMAFloatArrayCreateDP(ident,Array_Size,Initval)

    USE SMAArrays
    USE mkl_service
    
    IMPLICIT NONE
    
    INTEGER(KIND=8):: SMAFloatArrayCreateDP
    INTEGER(KIND=8):: alloc_size
    REAL(KIND=8),DIMENSION(*)::x
    POINTER(pt_x,x)
    INTEGER(KIND=4):: Array_Size
    INTEGER(KIND=4),INTENT(IN) :: ident
    REAL(KIND=8),INTENT(IN) :: INITVAL
    
    alloc_size=KIND(INITVAL)*Array_Size
    pt_x=mkl_malloc(alloc_size,KIND(INITVAL))
    x(1:Array_Size)=INITVAL
    PointerAdressesReal(ident)=pt_x
    SMAFloatArrayCreateDP=pt_x
    
END FUNCTION SMAFloatArrayCreateDP

FUNCTION SMAFloatArrayAccess(ident)

    USE SMAArrays
    
    IMPLICIT NONE
    
    INTEGER(KIND=8):: SMAFloatArrayAccess
    INTEGER(KIND=4),INTENT(IN) :: ident
    
    SMAFloatArrayAccess=PointerAdressesReal(ident)
    
END FUNCTION SMAFloatArrayAccess

SUBROUTINE SMAFloatArrayDelete(ident)

    USE SMAArrays
    
    IMPLICIT NONE
    
    INTEGER(KIND=4),INTENT(IN) :: ident
    
    CALL mkl_free(PointerAdressesReal(ident))
    
END SUBROUTINE SMAFloatArrayDelete
