	Subroutine SIGINI(SIGMA,COORDS,NTENS,NCRDS,NOEL,NPT,LAYER,
     1 KSPT,LREBAR,NAMES)
!-----------------------------------------------------------------------------
! Subroutine used to define the initial stress state. In this version
! The initial stress is linearly dependent on the vertical coordinate z, which
! is assumed to be zero at the surface and negative below. It is intended
! for use on with a homogenous soil with horizontal surface. The vertical stress
! is found from the weight of the soil, and the horizontal stresses are found
! from the K0 factor, Specific mass, gravity acceleration and K0 must be defined
! in the subroutine, i.e. they are not given as input parameters.
!
! INPUT (most explanations taken from ABAQUS help)
!
!	COORDS(ncrds)	An array containing the initial coordinates of this point.
!	NTENS			Number of stresses to be defined, which depends on the element type.
!	NCRDS			Number of coordinates.
!	NOEL			Element number.
!	NPT				Integration point number in the element.
!	LAYER			Layer number (for composite shells and layered solids).
!	KSPT			Section point number within the current layer.
!	LREBAR			Rebar flag. If LREBAR=1, the current integration point is associated with element rebar. Otherwise, LREBAR=0.
!	NAMES(1)		Name of the rebar to which the current integration point belongs, which is the name given in the rebar or rebar layer definition (“Defining reinforcement,” Section 2.2.3 of the ABAQUS Analysis User's Manual, or “Defining rebar as an element property,” Section 2.2.4 of the ABAQUS Analysis User's Manual). If no name was given in the rebar or rebar layer definition, this variable will be blank. This variable is relevant only when LREBAR=1.
!	NAMES(2)		Element type name (see Section I.1, “ABAQUS/Standard Element Index,” of the ABAQUS Analysis User's Manual).
!
! OUTPUT
!
!	Sigma(ntens)	Stress vector. z is taken to be the vertical coordinate,
!					  while x and y are horizontal.
!					  Some examples of the ordering of the stress components are
!					  3D: x, y, z, xy, xz, yz
!
!-----------------------------------------------------------------------------
! Johan Clausen
! Department of Civil Engineering
! Aalborg University
! January 2008
!-----------------------------------------------------------------------------
      INCLUDE 'ABA_PARAM.INC'
	
	integer(4) ntens, ncrds
	DIMENSION SIGMA(NTENS),COORDS(NCRDS)
      CHARACTER NAMES(2)*80
	real(8), parameter :: g	= 9.82_8   ! [m/s^2] Gravitational accelaration
      real(8), parameter :: rho = 2036.7_8 ! [kg/m^3] Mass density
      real(8), parameter :: K0  = 1.0_8	   ! [-] Earth pressure coefficient at rest
	
!------------------------------------------------------------
!		User coding
!------------------------------------------------------------
	Sigma = 0.0_8
	
	if (ntens == 6) then ! General 3D with z (or 3) as the vertical coordinate
		Sigma(3) = g*rho*Coords(3)
		Sigma(1) = K0*Sigma(3)
		Sigma(2) = Sigma(1)
	elseif (ntens == 4) then ! Plane situations with y (or 2) as the vertical coordinate
		Sigma(2) = g*rho*Coords(2)
		Sigma(1) = K0*Sigma(2)
		Sigma(3) = Sigma(1)
	end if
	
!-----------------------------------------------------------------------------      
	end subroutine SigIni !-------------------------------------------------
!-----------------------------------------------------------------------------

      subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
     1 RPL,DDSDDT,DRPLDE,DRPLDT,
     2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
     3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
     4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)

!-------------------------------------------------------------------------
! Plastic strain update for a linearly elastic - perfectly plastic material
! which obeys the Hoek-Brown yield criterion (and plastic potential).
!
! This version is intended for use with ABAQUS as a user programmable 
! feature. Because of this, the dummy arguments have the the names of the 
! example on user plasticity provided by ABAQUS.
!
! Input arguments:
!       
!   STRAN(NTENS)         -  total strain at beginning of increment
!   DSTRAN(NTENS)        -  strain increments
!   TIME(1)              -  step time at beginning of increment
!   TIME(2)              -  total time at beginning of increment
!   DTIME                -  time increment
!   TEMP                 -  temperature at beginning of increment
!   DTEMP                -  increment of temperature
!   PREDEF               -  interp. values of predef. field variables
!   DPRED                -  increment in predef. field variables
!   CMNAME               -  name given on *material option.
!   NDI                  -  no. of direct stress/strain components
!   NSHR                 -  no. of shear stress/strain components
!   NTENS                -  total no. of stress/strain components
!   NSTATV               -  no. of state variables
!   PROPS(NPROPS)        -  array with material properties
!   NPROPS               -  no. of material properties
!   COORDS(3)            -  coordinates at point
!   DROT(3,3)            -  rotation increment vector
!   CELENT               -  characteristic element length
!   DFGRD0(3,3)          -  deformation increment at start of step
!   DFGRD1(3,3)          -  deformation increment at end of step
!   NOEL                 -  element number
!   NPT                  -  integration point number
!   LAYER                -  layer number (not used)
!   KSPT                 -  section point number (not used)
!   KSTEP                -  step number
!   KINC                 -  increment number
!       
! Arguments that can be updated:
!
!   PNEWDT               -  ratio of suggested new time increment
!       
! Output arguments:
!
!   DDSDDE(NTENS,NTENS)  -  the elastoplastic constitutive tensor
!   STRESS(NTENS)        -  array with stresses
!   STATEV(NSTATEV)      -  statevariables**
!   SSE                  -  specific elastic strain energy
!   SPE                  -  specific plastic dissipation
!   SCD                  -  specific "creep" dissipation
!   RPL                  -  volumetric heat generation*
!   DDSDDT(NTENS)        -  stress variation due to temperature*
!   DRPLDE(NTENS)        -  variation of RPL due to strain*
!   DRPLDT               -  variation of RPL due to temperature*
!
! *:  Only required in fully coupled temperature-displacement analysis
! **: Region of the updated stress point (see [1]) is saved in STATEV(1):
!     
!	--------------------------------------------
!	 region  |  Return to    | from region no.
!     --------------------------------------------
!        0       no yielding
!        1       yield surface	region I
!        2       Edge 1			region II
!        3       Edge 2			region III
!        4       Apex point		region IV
!     --------------------------------------------
! REFERENCE:
! [1] Johan Clausen and Lars Damkilde: "An exact implementation of the Hoek-Brown
!	criterion for elasto-plastic finite element calculations", 
!	International Journal of Rock Mechanics & Mining Sciences, 2008,
!	volume 45, pages 831-847.
!-----------------------------------------------------------------------
! Johan Clausen
! Department of Civil Engineering
! Aalborg University
! September 2009
!-----------------------------------------------------------------------

      INCLUDE 'ABA_PARAM.INC'
      !implicit none

!-----------------------------------------------------------------------
!     Variables defined by ABAQUS input/output
!-----------------------------------------------------------------------
	
	CHARACTER*80 CMNAME
      integer(4) NTENS,NDI,NSHR,NSTATV,NPROPS,NOEL,NPT,LAYER,KSPT,
     &	KSTEP,KINC
      real(8) SPD,SCD,RPL,DRPLDT,DTIME,DTEMP,PNEWDT,CELENT
      real(8) STRESS(NTENS),STATEV(NSTATV),
     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
     3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)

!-----------------------------------------------------------------------
!     Variables defined by this subroutine
!-----------------------------------------------------------------------

      integer(4) region
      real(8) SigB(NTENS), SigC(NTENS),
     1 D(NTENS,NTENS), Dinv(NTENS,NTENS), Depc(NTENS,NTENS), K0,
     2 StrainElas(NTENS), DelStrainPlas(NTENS), 
     3 comp,m,s,a,mg,sg,ag, E, nu,
	4 DSTRES(NTENS), SSE, SPE, DEE, TDE, DPE_eq, temp, tens, tid

!---- Explicit interfaces ----------------------------------
	interface
		pure subroutine HoekBrownStressReturn(Sigma,nsigma,ElasPar,PlasPar,D,
     &												Dinv,Sigma_up,Depc,region)
			integer(4), intent(in) :: nsigma
			real(8), intent(in)    :: Sigma(nsigma), PlasPar(7), ElasPar(2)
			real(8), intent(in)    :: D(nsigma,nsigma), Dinv(nsigma,nsigma)
			
			real(8), intent(out)   :: Sigma_up(nsigma), Depc(nsigma,nsigma)
			integer(4), intent(out):: region
		end subroutine HoekBrownStressReturn
		!---
		pure subroutine DlinElas(E,nu,nsigma,D,Dinv)
			integer(4), intent(in) :: nsigma
			real(8), intent(in) :: E, nu

			real(8), intent(out) :: D(nsigma,nsigma)
			real(8), intent(out), optional :: Dinv(nsigma,nsigma)
		end subroutine DlinElas
	end interface
!-----------------------------------------------------------

!-----------------------------------------------------------------------
!	Some info about output to ABAQUS data and message files
!-----------------------------------------------------------------------
!
!	write(6,*) 'FILE ID# 6' ! .DAT-file
!	write(7,*) 'FILE ID# 7' ! .MSG-file
!
!-----------------------------------------------------------------------
	
!-----------------------------------------------------------------------
!     Material Properties
!-----------------------------------------------------------------------
	!--- Elasticity parameters --
      E      =  PROPS(1)  !  E    - Young's modulus
      nu     =  PROPS(2)  !  nu   - Poisson's ratio
      !----------------------------
      
      !--- Yield criterion paramters ----
      !--- f = sig1 - sig3 - comp*(s - m*sig1/comp)^a = 0
      comp   =  PROPS(3)  !  Uniaxial compressive strength of the intact rock
	m		 =  PROPS(4)  !  "friction" parameter of the rock mass
	s		 =  PROPS(5)  !  Hoek-Brown parameter
	a		 =  PROPS(6)  !  Curvature parameter in the Hoek-Brown criterion
	!----------------------------------
	
	!--- Plastic potential parameters ---
	!--- g = sig1 - sig3 - comp*(sg - mg*sig1/comp)^ag
	if (NPROPS > 6) then ! Possibly non-associated flow rule
		mg	 =  PROPS(7)  !  "dilation" parameter of the rock mass
		sg	 =  PROPS(8)  !  Hoek-Brown plastic potential parameter
		ag	 =  PROPS(9)  !  Curvature parameter in the Hoek-Brown plastic potential
	else ! Associated flow rule
		mg = m
		sg = s
		ag = a
	end if
	!------------------------------------
	
!-----------------------------------------------------------------------
!     Predictor stress
!-----------------------------------------------------------------------

      call DlinElas(E,nu,NTENS,D,Dinv)

	DSTRES = matmul(D,DSTRAN) ! stress increment
	SigB   = STRESS + DSTRES  ! elastic predictor stress

!-----------------------------------------------------------------------
!     Stress update and creation of consistent constitutive matrix
!-----------------------------------------------------------------------
	
	call HoekBrownStressReturn(SigB,NTENS,(/E,nu/),
     &					(/comp, m, s, a, mg, sg, ag/),
     &					D,Dinv,		SigC,Depc,region)
	
!-----------------------------------------------------------------------
!     Postconditioner for ABAQUS
!-----------------------------------------------------------------------

      DDSDDE    = Depc          ! consistent constitutive matrix
      STRESS    = SigC          ! updated stress at the Gauss point

!-----------------------------------------------------------------------
!	Saving of stress point region
!-----------------------------------------------------------------------

      STATEV(1) = real(region,8)  ! type of stress return is saved as a state variable

!-----------------------------------------------------------------------
      end subroutine UMAT
!-----------------------------------------------------------------------


	pure subroutine HoekBrownStressReturn(Sigma,nsigma,ElasPar,PlasPar,D,
     &								Dinv,Sigma_up,   Depc,region)
!--------------------------------------------------------------------------
! Wrapping subroutine for return mapping with a linearly elastic perfectly
! plastic Hoek-Brown material. Performed in principal coordinates according
! to Reference [1]
!
! The Hoek-Brown criterion is defined by
!		f = sigP(1) - sigP(3) - comp*(s - m*sigP(1)/comp)^a = 0
!
! and the corresponding plastic potential
!       g = sigP(1) - sigP(3) - comp*(sg - mg*sigP(1)/comp)^ag
! 
! INPUT
!  - Name -		-type and description --
!	Sigma		real array(nsigma). Stress vector. Size is dependent on the type of stress.
!				  The ordering of the components must be as follows:
!				    Plane situations (plane stress, plane strain and axisymmetry):
!					Sigma = [sig_x sig_y sig_z tau_xy] (i.e. sig_z is out-of-plane)
!					General 3D:  Sigma = [sig_x sig_y sig_z tau_xy tau_xz tau_yz]
!	nsigma		integer,scalar. Number of stress components, i.e. the size of Sigma
!	ElasPar		real, array(2).	[E nu] Elastic parameters
!	PlasPar		real, array(7). Vector of plastic material parameters [comp, m, s, a, mg, sg, ag]
!	D			real, array(nsigma,nsigma) Elastic constitutive matrix.
!	Dinv		real, array(nsigma,nsigma) Inverted elastic constitutive matrix (compliance matrix).
!
! OUTPUT
!	Sigma_up	real, array(nsigma) Stress vector containg the updated stresses.
!	Depc		real, array(nsigma,nsigma) Consistent constitutive elasto-plastic matrix
!				  If the material is not yielding Depc = D
!	region		integer, scalar. Number of the region of the returned stress
!					region = 0: the elastic region, i.e. no yielding
!					region = 1: Return to the Hoek-Brown surface
!					region = 2: Return to the triaxial compressive line, sigp_1 = sigp_2
!					region = 3: Return to the triaxial tensile line, sigp_2 = sigp_3
!					region = 4: Return to the apex, sigp_1 = sigp_2 = sigp_3.
!
! REFERENCES:
! [1] Johan Clausen and Lars Damkilde: "An exact implementation of the Hoek-Brown
!	  criterion for elasto-plastic finite element calculations", 
!	  International Journal of Rock Mechanics & Mining Sciences, 2008,
!	  volume 45, pages 831-847.
!------------------------------------------------------------------------
! Johan Clausen
! Department of Civil Engineering
! Aalborg University
! April 2008
!------------------------------------------------------------------------
	implicit none

	!-----------------------
	integer(4), intent(in) :: nsigma
	real(8), intent(in)    :: Sigma(nsigma), PlasPar(7), ElasPar(2)
	real(8), intent(in)    :: D(nsigma,nsigma), Dinv(nsigma,nsigma)

	real(8), intent(out)	 :: Sigma_up(nsigma), Depc(nsigma,nsigma)
	integer(4), intent(out):: region
	!-----------------------

	integer(4) hoop,ouplP, s1, s2, nshear, Lt(3,3)
	real(8) psi(3,3), f, SigP(3), SigP_up(3), nu
	real(8), dimension(nsigma,nsigma) :: DepP, DepcP, A, Atrans
	real(8), dimension(nsigma,nsigma) :: MatProd, T, Dc, Dcinv, Tinv
	real(8) Temp(3,3), Norm(3), Edir(3), NormS(3), EdirS(3), R(3), Rg(3)
	real(8) Hessian(3,3), Gdir(3), Fdir(3), SiPla(3), Tnorm(3,3)
	real(8) Tinvnorm(3,3)
	real(8) Tshear(nsigma-3,nsigma-3),Tshearinv(nsigma-3,nsigma-3)
!------------------------------------------------------------------------

	!--- Settings regarding calculation of the constitutive matrix, DepP -------------
	integer(4), parameter :: DepFormLine = 1 ! Form of the consistent constitutive matrix on a line
					! DepFormLine = 0: Double singular constitutive matrix
					! DepFormLine = 1: Modified "almost double singular" matrix on a line ("Tilnæ_b=axr" in the MatLab code)
					! DepFormLine = 2: A single singular matrix in the Koiter direction, i.e. the direction of the plastic strain.
	real(8), parameter :: beta = 100.0_8 ! The value used in the modification of DepP when DepFormLine = 1
	
	integer(4), parameter :: DepFormApex = 2 ! Form of the consistent constitutive matrix on the apex
					! DepFormApex = 0: DepP is the zero matrix
					! DepFormApex = 1: A modified single singular DepP in the Koiter direction ("DepKo/fak" in the MatLab code)
					! DepFormApex = 2: A modified double singular DepP in the Koiter and the hydrostatic direction ("DepKoiteP" in the MatLab code)
	real(8), parameter :: alpha = 1000.0_8 ! ! Factor in the modified apex-Depc-formulation. Used when DepFormApex = 1 or DepFormApex = 2
	!--------------------------------------------------------------------------------

!---- Explicit interfaces ----------------------------------
	interface
		pure subroutine PrinStressAna(Sigma,nsigma,SigP)
			integer(4), intent(in) :: nsigma
			real(8), intent(in) :: Sigma(nsigma)

			real(8), intent(out) :: SigP(3)
		end subroutine PrinStressAna
		!---
		pure subroutine PrinRetHoekBrown(SigP,nu,PlasPar,Dfull,Dinvfull,
     &			nsigma,	  SigP_up,region,Gdir,Fdir)
			integer(4), intent(in) :: nsigma
			real(8), intent(in)	   :: SigP(3), PlasPar(7), nu
			real(8), dimension(nsigma,nsigma), intent(in) :: Dfull, Dinvfull

			real(8), intent(out)	:: SigP_up(3), Gdir(3), Fdir(3)
			integer(4), intent(out) :: region
		end subroutine PrinRetHoekBrown
		!---
		pure subroutine TnormHoekBrownPrin(SigP_up,SiPla,PlasPar,Dn,Dinvn,
     &							region,Gdir,Fdir,  Tnorm,Tnorminv)
			real(8), intent(in) :: SigP_up(3), PlasPar(7), Dn(3,3), Dinvn(3,3)
			real(8), intent(in) :: Sipla(3), Gdir(3), Fdir(3)
			integer(4), intent(in) :: region

			real(8), intent(out) :: Tnorm(3,3), Tnorminv(3,3)
		end subroutine TnormHoekBrownPrin
		!---
		pure subroutine TshearPrinPerfect(SigP,SigP_up,nshear,s1,s2,
     &												Tshear,Tshearinv)
			real(8), intent(in) :: SigP(3), SigP_up(3)
			integer(4), intent(in) :: nshear, s1, s2

			real(8), dimension(nshear,nshear), intent(out) :: Tshear, Tshearinv
		end subroutine TshearPrinPerfect
		!---
		pure subroutine FormDepPerfect(D,Norm,Edir,nsigma,Dep)
			integer(4), intent(in) :: nsigma
			real(8), intent(in) :: D(nsigma,nsigma), Norm(nsigma)
			real(8), intent(in) :: Edir(nsigma)
			real(8), intent(out) :: Dep(nsigma,nsigma)
		end subroutine FormDepPerfect
		!---
		pure subroutine FormModDepLine(mod_type,beta,SiPla,Fline,Gline,Dninv,
     1											Dc,Dcinv,nsigma,Depc)
			integer(4), intent(in) :: mod_type, nsigma
			real(8), intent(in) :: Dninv(3,3),Dc(nsigma,nsigma),Gline(3)
			real(8), intent(in) :: beta,Sipla(3),Fline(3),Dcinv(nsigma,nsigma)
			real(8), intent(out) :: Depc(nsigma,nsigma)
		end subroutine FormModDepLine
		!---
		pure subroutine FormModDepApex(mod_type,alpha,SiPla,Dninv,Dc,Dcinv,
     1												nsigma,ishydro,  Depc)
			integer, intent(in) :: mod_type, nsigma, ishydro
			real(8), intent(in) :: alpha,Sipla(3),Dninv(3,3),Dc(nsigma,nsigma),
     1										Dcinv(nsigma,nsigma)
			real(8), intent(out) :: Depc(nsigma,nsigma)
		end subroutine FormModDepApex
		!---
		pure subroutine PrinDirect(Sigma,nsigma,SigP,psi)
			integer(4), intent(in) :: nsigma
			real(8), intent(in) :: Sigma(nsigma), SigP(3)
			real(8), intent(out)  :: psi(3,3)
		end subroutine PrinDirect
		!---
		pure subroutine TransMatrix(psi,nsigma,ouplP,A)
			integer(4), intent(in) :: nsigma, ouplP
			real(8), intent(in)    :: psi(3,3)
			real(8), intent(out)   :: A(nsigma,nsigma)
		end subroutine TransMatrix
	end interface
!-----------------------------------------------------------
	
	call PrinStressAna(Sigma,nsigma,SigP) ! Principal stresses SigP

	!--- Value of Yield function f ---  
		! PlasPar = [comp, m, s, a, mg, sg, ag]
	if (SigP(1) >= PlasPar(3)*PlasPar(1)/PlasPar(2)) then ! i.e. if SigP(1) >= apex
		f = 1.0_8
	else								
		! i.e.:	f = sigP1 - sigP3 - comp*(s-m*sigP1/comp)^a
		f = SigP(1) - SigP(3) - PlasPar(1)*
     &			(PlasPar(3) - PlasPar(2)*SigP(1)/PlasPar(1))**PlasPar(4) ! Hoek-Brown yield criterion
	end if
	
	!---------------------------------

	if (f > 1.0e-7_8) then
!	--- Position of out-of plane or hoop stress in principal stress vector SigP ----
		if (nsigma == 4) then
			if (Sigma(3) == SigP(1)) then
				ouplP = 1 ; s1 = 2 ; s2 = 3
			elseif (Sigma(3) == SigP(2)) then
				ouplP = 2 ; s1 = 1 ; s2 = 3
			elseif (Sigma(3) == SigP(3)) then
				ouplP = 3 ; s1 = 1 ; s2 = 2
			end if
		end if
!	--------------------------------------------------------------------------------


!	--- Stress return in principal coordinates --
		nu = ElasPar(2) ! Poissons ratio
		call PrinRetHoekBrown(SigP,nu,PlasPar,D,Dinv,nsigma,
     &									SigP_up,region,Gdir,Fdir)
     		SiPla = SigP - SigP_up ! Plastic corrector stress in principal stress space
!	---------------------------------------------

!	--- Consistent constitutive matrix -------
		call TnormHoekBrownPrin(SigP_up,SiPla,PlasPar,D(1:3,1:3),Dinv(1:3,1:3)
     &											,region,Gdir,Fdir,Tnorm,Tinvnorm) ! Forms the normal components of the T-matrix
		nshear = nsigma - 3 ! Number of shear components
		call TshearPrinPerfect(SigP,SigP_up,nshear,s1,s2,Tshear,Tshearinv) ! Forms the shear part modification matrix T
	
		T = 0.0_8 ; Tinv = 0.0_8
		T(1:3,1:3) = Tnorm ; T(4:nsigma,4:nsigma) = Tshear
		T	inv(1:3,1:3) = Tinvnorm ; Tinv(4:nsigma,4:nsigma) = Tshearinv

		Dc = matmul(T,D) ! Modified elastic stiffness
	
		DepcP = 0.0_8 ! Initiation of the consistent constitutive matrix in principal stress space
		if (region == 1) then ! Return to the yield surface -------------------------------
			DepcP(4:nsigma,4:nsigma) = Dc(4:nsigma,4:nsigma)  ! Shear components of consistent constitutive matrix in principal stress space
			call FormDepPerfect(Dc(1:3,1:3),Fdir,Gdir,3,Temp) ! Normal components of consistent constitutive matrix in principal stress space
			DepcP(1:3,1:3) = Temp   ! Consistent constitutive matrix in principal stress space

		elseif (region == 2 .or. region == 3) then ! return to a line ---------------------
			Dcinv = matmul(Dinv,Tinv) ! Inverse of Dc
			call FormModDepLine(DepFormLine,beta,SiPla,Fdir,Gdir,Dinv(1:3,1:3),
     &													    Dc,Dcinv,nsigma,DepcP) ! Consistent constitutive matrix in principal stress space
		else ! apex return ----------------------------------------------------------------
			Dcinv = Matmul(Dinv,Tinv) ! Inverse modified elastic stiffness
			call FormModDepApex(DepFormApex,alpha,SiPla,Dinv(1:3,1:3),Dc,Dcinv,
     &																nsigma,1,DepcP)
		end if ! region = ... -------------------------------------------------------------
		
!	--- Tranformation matrix A ------------------
		call PrinDirect(Sigma,nsigma,SigP,psi)
		call TransMatrix(psi,nsigma,ouplP,A)
!	---------------------------------------------

!	--- Coordinate transformation ---------------
		Atrans = transpose(A)
		Sigma_up = 0.0_8
		Sigma_up = matmul(Atrans(:,1:3),SigP_up)
		MatProd = matmul(Atrans,DepcP)
		Depc = matmul(MatProd,A)
!	---------------------------------------------

	else ! no yielding
		Sigma_up = Sigma
		Depc = D
		region = 0
	end if ! f > 0
			
!--------------------------------------------------------------------------------
	end subroutine HoekBrownStressReturn !-------------------------------------
!--------------------------------------------------------------------------------

	pure subroutine PrinRetHoekBrown(SigP,nu,PlasPar,Dfull,Dinvfull,nsigma,
     &												SigP_up,region,Gdir,Fdir)
!-----------------------------------------------------------------------
! Stress return for Hoek-Brown plasticity in principal stress
! space. Includes non-associated plasticity, but no hardening is allowed.
! The Hoek-Brown criterion is written as
!				f = sigP_1 - sigP_3 - comp*(s-m*sigP_1/comp)^a = 0
! The plastic potential is
!				g = sigP_1 - sigP_3 - comp*(sg-m*sigP_1/comp)^ag
!
! INPUT
!  - Name -		-type and description --
!	SigP		(real,array(3))	Principal predictor stresses in descending order
!				  SigP = [sigP_1 sigP_2 sigP_3]
!	nu			(real,scalar) Poisson's ratio
!	PlasPar		(real,array(7))	Vector of plastic material parameters
!				  PlasPar = [comp, m, s, a, mg, sg, ag]
!				  comp: yp, Uniaxial compressive strength of the inctact rock.
!				  m: yp, Parameter which has a large saying in the amount of material friction
!				  s: yp, Has a large saying in determination of tension strength
!				  a: yp, Curvature parameter. For rock usually 0.5 <= a <= 0.65
!				  mg:pp, The potential equivalent of the yield parameter m
!				  sg:pp, The potential equivalent of the yield parameter s
!				  ag:pp, Curvature parameter of plastic potential
!	Dfull		(real,array(nsigma,nsigma)) Elastic constitutive matrix
!	Dinvfull	(real,array(nsigma,nsigma)) The inverse of the elastic constitutive matrix (the compliance matrix)
!	nsigma		(integer,scalar) Number of stress components
!
! OUTPUT
!	SigP_up	(real,array(3))	Updated principal stresses, i.e. they obey the
!			  Hoek-Brown yield criterion
!	region	(integer,scalar) Type of return, i.e. the region of the
!			  predictor stress
!	Gdir	(real,array(3)) A characteristic plastic potential direction
!			  vector in principal stress space, depending on the return:
!				region = 1 (surface): Plastic potential normal
!				region = 2 (line 1, tri. comp.): Potential line 1 direction
!				region = 3 (line 2, tri. tens.): Potential line 2 direction
!				region = 4 (apex return): Hydrostatic line direction
!	Fdir	(real,array(3)) A characteristic yield surface direction
!			  vector in principal stress space, depending on the return:
!				region = 1 (surface): Yield surface normal
!				region = 2 (line 1, tri. comp.): Yield line 1 direction
!				region = 3 (line 2, tri. tens.): Yield line 2 direction
!				region = 4 (apex return): Hydrostatic line direction
!
! REFERENCES:
! [1] Johan Clausen and Lars Damkilde: "An exact implementation of the Hoek-Brown
!	  criterion for elasto-plastic finite element calculations", 
!	  International Journal Rock Mechanics & Mining Sciences, 2008,
!	  volume 45, pages 831-847.
!----------------------------------------------------------------------
! Johan Clausen
! Department of Civil Engineering
! Aalborg University
! April 2008
!----------------------------------------------------------------------
	implicit none

	!-------------------
	integer(4), intent(in) :: nsigma
	real(8), intent(in)	   :: SigP(3), PlasPar(7), nu
	real(8), dimension(nsigma,nsigma), intent(in) :: Dfull, Dinvfull

	real(8), intent(out)	:: SigP_up(3), Gdir(3), Fdir(3)
	integer(4), intent(out) :: region
	!-------------------

	real(8) comp, m, a, s, mg, ag, sg, apex, apexg
	real(8) s1_0, s1_1, parf, parg, k, kg, dkg, alfas, alfar
	real(8) dalfar, dalfas, h0, dh, tol, t0, kga
	real(8) s3_0, s3_1
	real(8), dimension(3) :: Norm, Edir, R1, R2, Rg, DRg
	real(8), dimension(3) :: N_IV_II, N_IV_III, Ra_II, Ra_III
	real(8), dimension(3) :: SiPla, PlasStrain
	real(8), dimension(3) :: Rp_up, ApexVec, Sig0, Rpa
	real(8), dimension(3,3) :: D, Dinv
	integer(4) i
	real(8), parameter :: tolfak = 1.0e-8_8 ! tolerance factor
	integer(4), parameter :: imax = 50		! Maximum number of iterations
	integer(4), parameter :: bisec = 12		! Maximum number of Newton-Raphson iterations before bisection is invoked

!----- Explicit interfaces -----------------------------------
	interface
		pure function Sigma1HoekBrown(sigma3,comp,m,s,a)
			real(8), intent(in) :: sigma3, comp, m, s, a

			real(8) Sigma1HoekBrown
		end function Sigma1HoekBrown
		!---
		pure function Cross(A,B)
			real(8), intent(in)	::	A(3), B(3)
			real(8)				::  Cross(3)
		end function Cross
	end interface
!--------------------------------------------------------------
!----- Stress return ----------------------------------------------------------------

    !--- Preliminary parameters ----------------
	!PlasPar = [comp, m, s, a, mg, sg, ag]
	comp = PlasPar(1) ! Compressive strength parameter
	m	 = PlasPar(2) ! Parameter which has a large saying in the amount of material friction
	s	 = PlasPar(3) ! Has a large saying in determination of tension strength
	a	 = PlasPar(4) ! Curvature yield parameter. For rock usually 0.5 <= a <= 0.65
	apex = s*comp/m	  ! Stress coordinate of the criterion's apex

	mg	 = PlasPar(5) ! The potential equivalent of the yield parameter m
	sg	 = PlasPar(6) ! The potential equivalent of the yield parameter s
	ag	 = PlasPar(7) ! Curvature parameter of plastic potential
	apexg= sg*comp/mg ! Stress coordinate of the apex of the plastic potential

	D	 = Dfull(1:3,1:3)	 ! Relates to normal stresses
	Dinv = Dinvfull(1:3,1:3) ! Relates to normal stresses


	ApexVec = (/apex, apex, apex/) ! Apex coordinate vector


	if (sg-mg*s/m > 0) then
		kga = 1.0_8 + ag*mg*(sg-mg*s/m)**(ag-1.0_8) ! Value of kg at the apex
    		Rpa = matmul(D,(/kga, 0.0_8, -1.0_8/))		! Plastic corrector stress at apex by return from region I
    		Ra_II = matmul(D,(/1.0_8, 1.0_8, -2/kga/))	! Plastic corrector stress at apex by return along the plane sigma1 = sigma2
    		Ra_III = matmul(D,(/2*kga, -1.0_8, -1.0_8/))! Plastic corrector stress at apex by return along the plane sigma2 = sigma3
    		N_IV_II = Cross(Rpa,Ra_II)			! Normal of the region II and IV boundary plane
    		N_IV_III = Cross(Ra_III,Rpa)		! Normal of the region III and IV boundary plane
	else
		N_IV_II = (/-nu, -nu, 1.0_8/)   ! Normal of the region II and IV boundary plane
    		N_IV_III = (/0.0_8, 0.0_8, 0.0_8/)		! Normal of the region III and IV boundary plane
	end if
    !-------------------------------------------

    !--- Region IV - apex return ------------------
	if (dot_product(N_IV_II,(SigP-ApexVec))  >= 0.0_8  .and.
     &	  dot_product(N_IV_III,(SigP-ApexVec)) >= 0.0_8) then
     		region = 4
		SigP_up = ApexVec		! Updated stress
		SiPla = SigP - SigP_up	! Plastic corrector stress
    !----------------------------------------------

    !--- Region I - surface return -----------------------------------
	else
		region = 1

		if (SigP(3) < apex) then ! Initial guess
			s1_1 = Sigma1HoekBrown(SigP(3),comp,m,s,a)
    		else
			s1_1 = 0.95*apex
    		end if
            
    		if (SigP(3) /= 0.0_8) then
			tol = tolfak*abs(SigP(3)) ! Tolerance
    		else
    			tol = tolfak*comp ! Tolerance
    		end if
		if (abs(s1_1-SigP(1)) > tol) then
			s1_0 = s1_1 + 2*tol
		else
			s1_0 = s1_1
		end if	
    	
    		i = 0
    		do while (abs(s1_1 - s1_0) > tol) ! Newton-Raphson iteration
			s1_0 = s1_1
			
			parf = s - m*s1_0/comp			! Value of f-parenthesis
			parg = sg - mg*s1_0/comp		! Value of g-parenthesis
			
			s3_0 = s1_0 - comp*parf**a		! Corresponding sigmap3-value
			k =  1.0_8 + a*m*parf**(a-1.0_8)		! First component of yield surface normal
			kg = 1.0_8 + ag*mg*parg**(ag-1.0_8)		! First component of plastic potential gradient
			dkg = (1.0_8-ag)/comp*ag*mg**2*parg**(ag-2.0_8) ! Derivative of gradients first component
			
			dalfar = (nu*dkg*((1-nu)*kg-nu) - (kg*nu-1+nu)*(1-nu)*dkg)/
     &														((1-nu)*kg-nu)**2 ! Derivative of alfar
     			alfar = (nu*kg-1+nu)/(kg*(1-nu)-nu)		  ! Slope of the gradient
			
			dalfas = (k*(s1_0 - SigP(1)) + SigP(3) - s3_0)/(s1_1 - SigP(1))**2	! Derivative of alfas
			alfas = (s3_0 - SigP(3))/(s1_1 - SigP(1)) ! Slope based on stress points
			
			h0 = alfas - alfar   ! function value at s1_0
			dh = dalfas - dalfar ! function derivative at s1_0
			
			s1_1 = dble(s1_0 - h0/dh)  ! New value for sigmap1 (first principal stress, s1_1)
			   
			if (s1_1 > apex) then
				s1_1 = 0.9_8*apex + 0.1_8*s1_0
			end if
			
			i = i + 1
			if (i > imax) then
				s1_0 = apex
				s1_1 = apex
				Norm = (/1.0_8, 0.0_8, 0.0_8/)
				if (apex == apexg) then
					Edir = (/1.0_8, 0.0_8, 0.0_8/)
					parg = 0.0_8		
				else
					parg = sg - mg*s/m
					Edir = 0.0_8
					Edir(1) = 1.0_8 + ag*mg*parg**(ag-1.0_8)
				end if
			end if
    		end do ! while abs(s1_1 - s1_0) > tol % Newton-Raphson iteration
    		
    		SigP_up(1) = dble(s1_1) ! Updated largest principal stress
    		SigP_up(3) = dble(s1_1 - comp*(s - m*s1_1/comp)**a) ! Updated smallest principal stress
            
    		if (i <= imax) then
			k = 1.0_8 + a*m*(s - m*SigP_up(1)/comp)**(a-1.0_8) ! first component of the yield normal
			parg = sg - mg*SigP_up(1)/comp  ! Value of the parenthesis
			kg = 1.0_8 + ag*mg*parg**(ag-1.0_8) ! First component of the gradient
        
			Norm = dble((/k,  0.0_8, -1.0_8/)) ! Yield surface normal at updated stress
			Edir = dble((/kg, 0.0_8, -1.0_8/)) ! Gradient at updated stress
        
			Rp_up = matmul(D,Edir)				! Direction of plastic corrector
			t0 = (SigP_up(1) - SigP(1))/Rp_up(1)! Parameter
			SigP_up(2) = t0*Rp_up(2) + SigP(2)	! Updated intermediate principal stress
			Gdir = Edir ; Fdir = Norm
    		else
			SigP_up = ApexVec
			Gdir = (/1.0_8, 1.0_8, 1.0_8/) ; Fdir = (/1.0_8, 1.0_8, 1.0_8/)
    		end if ! (i <= imax)
	end if ! Region IV or I stress return
    !-----------------------------------------------------------------

    !---- Region II ------------------------------ 
	if (SigP_up(2) > SigP_up(1) .and.
     &		min(SigP_up(2),apex) 
     &		- comp*(s-m*min(SigP_up(2),apex)/comp)**a >= SigP_up(3)) then ! Region II, return to curve 1
    		region = 2
    		
    		s1_1 = 0.98*SigP_up(1) ! Initial estimate
       
    		s1_0 = s1_1 + 2*tol
    		i = 0
    		do while (abs(s1_1 - s1_0) > tol) ! Newton-Raphson iteration
			s1_0 = s1_1

			parf = s - m*s1_0/comp			 ! Value of the f-parenthesis
			parg = sg - mg*s1_0/comp		 ! Value of the g-parenthesis
        
			k =  1.0_8 + a*m*parf**(a-1.0_8)         ! Component of yield surface normal
			kg = 1.0_8 + ag*mg*parg**(ag-1.0_8)      ! Component of plastic potential gradient
		
			dkg = (1.0_8-ag)*ag*mg**2/comp*parg**(ag-2.0_8) ! Derivative of potential gradient

			Sig0 = (/s1_0, s1_0, s1_0-comp*parf**a/) ! Stress point of the curve
			Rg = (/1.0_8, 1.0_8, kg/) ! Direction vector of the potential curve
		
			SiPla = SigP - Sig0 ! Temporary value of the plastic corrector
			PlasStrain = matmul(Dinv,SiPla) ! Current value of plastic strain increment

			h0 = dot_product(Rg,PlasStrain) ! function value at s1_0

			DRg = (/0.0_8, 0.0_8, dkg/) ! Derivative of the line 1 direction

			R1 = -(/1.0_8, 1.0_8, k/) ! Negative direction vector of curve 1 at s1_0

			dh = dot_product(DRg,PlasStrain) + dot_product(Rg,matmul(Dinv,R1)) ! Derivative of h in s1_0

			s1_1 = dble(s1_0 - h0/dh) ! New estimate of s1_1

			if (s1_1 > apex) then
				s1_1 = 0.95_8*apex + 0.05_8*dble(s1_0)
			end if

			i = i + 1
			if (i > imax) then
				s1_0 = apex
				s1_1 = apex
				region = 4
			end if

    		end do ! while (abs(s1_1 - s1_0) > tol) ! Newton-Raphson iteration
		parf = s - m*s1_1/comp
		parg = sg - mg*s1_1/comp
    		SigP_up = dble((/s1_1, s1_1, s1_1-comp*parf**a /)) ! Updated principal stress

		k =  1 + a*m*parf**(a-1)     ! Updated component of yield surface normal
		kg = 1 + ag*mg*parg**(ag-1)	 ! Updated component of plastic potential gradient
		Gdir = (/1.0_8, 1.0_8, kg/) ; Fdir = (/1.0_8, 1.0_8, k/) ! Updated line directions
    !--------------------------------------------------------

    !---- Region III --------------------------------- 
	elseif (SigP_up(2) < SigP_up(3)) then ! Region III, return to curve 2
    		region = 3
		s1_1 = 0.99*SigP_up(1) ! Initial estimate

    		s1_0 = s1_1 + 2*tol
    		i = 0
    		do while (abs(s1_1 - s1_0) > tol) ! Newton-Raphson iteration
			s1_0 = s1_1

			parf = s - m*s1_0/comp			 ! Value of the f-parenthesis
			parg = sg - mg*s1_0/comp		 ! Value of the g-parenthesis
        
			k =  1.0_8 + a*m*parf**(a-1.0_8)    ! Component of yield surface normal
			kg = 1.0_8 + ag*mg*parg**(ag-1.0_8) ! Component of plastic potential gradient
		
			dkg = (1.0_8-ag)*ag*mg**2/comp*parg**(ag-2.0_8) ! Derivative of potential gradient

			Sig0 = (/s1_0, s1_0-comp*parf**a, s1_0-comp*parf**a/) ! Stress point of the curve
			Rg = (/1.0_8, kg, kg/) ! Direction vector of the potential curve
		
			SiPla = SigP - Sig0 ! Temporary value of the plastic corrector
			PlasStrain = matmul(Dinv,SiPla) ! Current value of plastic strain increment

			h0 = dot_product(Rg,PlasStrain) ! function value at s1_0

			DRg = (/0.0_8, dkg, dkg/) ! Derivative of the line 2 direction

			R2 = -(/1.0_8, k, k/) ! Negative direction vector of curve 2 at s1_0

			dh = dot_product(DRg,PlasStrain) + dot_product(Rg,matmul(Dinv,R2)) ! Derivative of h in s1_0

			s1_1 = dble(s1_0 - h0/dh) ! New estimate of s1_1

			if (s1_1 > apex) then
				s1_1 = 0.95_8*apex + 0.05_8*dble(s1_0)
			end if

			i = i + 1
			if (i > imax) then
				s1_0 = apex
				s1_1 = apex
				region = 4
			end if
    		end do ! while (abs(s1_1 - s1_0) > tol) ! Newton-Raphson iteration
		parf = s - m*s1_1/comp
		parg = sg - mg*s1_1/comp
    		SigP_up = dble((/s1_1, s1_1-comp*parf**a, s1_1-comp*parf**a /)) ! Updated principal stress

		k =  1.0_8 + a*m*parf**(a-1.0_8)     ! Updated component of yield surface normal
		kg = 1.0_8 + ag*mg*parg**(ag-1.0_8)	 ! Updated component of plastic potential gradient
		Gdir = (/1.0_8, kg, kg/) ; Fdir = (/1.0_8, k, k/) ! Updated line directions
	end if ! (SigP_up(2) > SigP_up(1)) Region II or III
	!--------------------------------------------------------

!--------------------------------------------------------------------------------
	end subroutine PrinRetHoekBrown !------------------------------------------
!--------------------------------------------------------------------------------

	pure function Sigma1HoekBrown(sigma3,comp,m,s,a)
!--------------------------------------------------------------------------
! Determines the largest principal stress in a Hoek-Brown material when
! the smallest principal stress sigma3 is given (tension is positive)
! The Hoek-Brown criterion with is given as
! f = sigma1 - sigma3 - comp(s - m*sigma1/comp)^a = 0
! The solution is based on Newton-Raphson iteration
!
! INPUT
!  - Name -	-type,size - - Description -
!	sigma3	(real,sc)	 Smallest principal stress
!	comp	(real,sc)	 Compressive strength of intact rock
!	m		(real,sc)	 Hoek-Brown parameter
!	s		(real,sc)	 Hoek-Brown parameter
!	a		(real,sc)	 Curvature parameter. For rock usually 0.5 <= a <= 0.65
!
! OUTPUT
!	Sigma1HoekBrown	(real,sc)	Corresponding largest principal stress
!
!------------------------------------------------------------------------
! Johan Clausen
! Department of Civil Engineering
! Aalborg University
! April 2008
!------------------------------------------------------------------------
	implicit none
	real(8), intent(in) :: sigma3, comp, m, s, a

	real(8) Sigma1HoekBrown
	!----------------------
	real(8) apex, tol, x0, x1, par, f0, df
	integer(4) i

!------------------------------------------------------------------------
	apex = s*comp/m ! apex stress
	if (sigma3 < apex) then
		if (sigma3 /= 0.0_8) then
				tol = abs(1.0e-10*sigma3) !Tolerance
    		else 
			tol = 1.0e-10
    		end if

		x0 = 0.0_8
    		x1 = sigma3 + (1+1.0e-7_8)*tol ! First estimate
    
    		i = 0
    		do while (abs(x1 - x0) > tol)
			x0 = x1
			par = s - m*x0/comp            ! Value in parenthesis
			f0 = x0 - sigma3 - comp*par**a ! Value of the function at x0
			df = 1.0_8 + a*m*par**(a-1.0_8)! Derivative at x0
			x1 = x0 - f0/df				   ! New approximative value of x1 at f0 = 0
			i = i + 1
    
			if (x1 > apex) then
				x1 = (apex+x0)/2
			end if
			if (i > 50) then
				x1 = apex
				x0 = apex
			end if
    		end do
    		Sigma1HoekBrown = x1
	else
		Sigma1HoekBrown = apex
	end if 
!----------------------------------------------------------------------------------
	end function Sigma1HoekBrown !-----------------------------------------------
!----------------------------------------------------------------------------------

	pure subroutine TnormHoekBrownPrin(SigP_up,SiPla,PlasPar,Dn,Dinvn,
     &							region,Gdir,Fdir,  Tnorm,Tnorminv)
!--------------------------------------------------------------------------
! Forms the upper left 3 x 3 components of the modification matrix T for use
! with the Hoek-Brown criterion i.e. the components of
!
!	[Tnorm] = ([I] + lambda*D*d^2g/dsigma^2)^-1
!
! where g is the plastic potential and sigma are the normal stresses. The
! modification matrix, T, is used in the construction of the consistent
! constitutive matrix, Depc.
! All stresses are expressed in principal stress space. In addition, the
! derivative of the plastic potential, i.e. the strain direction vector
! Edir and the yield surface normal Norm are also calculated for the case
! of a yield surface return (i.e. region 1).
! The Hoek-Brown criterion is written in the form
!			f = sigP_1 - sigP_3 - comp*(s-m*sigP_1/comp)^a = 0
! The plastic potential is g = sigP_1 - sigP_3 - mg/m*comp*(sg-m*sigP_1/comp)^ag
!
!
! INPUT
!  - Name -	-type and description --
!	SigP_up	(real,array(3))	Updated principal stresses, i.e. they obey the
!			  Hoek-Brown yield criterion
!	SiPla	(real,array(3)) Plastic corrector stress
!	PlasPar	(real,array(7))	Vector of plastic material parameters
!			  PlasPar = [comp, m, s, a, mg, sg, ag]
!			  comp:  Uniaxial compressive strength of the inctact rock.
!			  m:  Parameter which has a large saying in the amount of material friction
!			  s:  Has a large saying in determination of tension strength
!			  a:  Curvature parameter. For rock usually 0.5 <= a <= 0.65
!			  mg: The potential equivalent of the yield parameter m
!			  sg: The potential equivalent of the yield parameter s
!			  ag: Curvature parameter of plastic potential
!	Dn		(real,array(3,3)) Normal part of the elastic constitutive matrix
!	Dinvn	(real,array(3,3)) Normal part of the inverse elastic constitutive matrix
!	region	(integer,scalar) Indicates on which region the updated stress
!			  point is located:
!			  region = 1 : The yield surface
!			  region = 2 : Line 1; Triaxial compression
!			  region = 3 : Line 2; Triaxial extension
!			  region = 4 : The apex
!	Gdir	(real,array(3)) A characteristic plastic potential direction
!			  vector in principal stress space, depending on the return:
!				region = 1 (surface): Plastic potential normal
!				region = 2 (line 1, tri. comp.): Potential line 1 direction
!				region = 3 (line 2, tri. tens.): Potential line 2 direction
!				region = 4 (apex return): Hydrostatic line direction
!	Fdir	(real,array(3)) A characteristic yield surface direction
!			  vector in principal stress space, depending on the return:
!				region = 1 (surface): Yield surface normal
!				region = 2 (line 1, tri. comp.): Yield line 1 direction
!				region = 3 (line 2, tri. tens.): Yield line 2 direction
!				region = 4 (apex return): Hydrostatic line direction
!
! OUTPUT
!	Tnorm	(real,array(3,3)) Normal elements of the modification matrix T
!	Tnorminv(real,array(3,3)) The inverse of Tnorm
!
! REFERENCES:
! [1] Johan Clausen and Lars Damkilde: "An exact implementation of the Hoek-Brown
!	  criterion for elasto-plastic finite element calculations", 
!	  International Journal Rock Mechanics & Mining Sciences, 2008,
!	  volume 45, pages 831-847.
! [2] Johan Clausen: Efficient Non-Linear Finite Element Implementation of
!	  Elasto-Plasticity for Geotechnical Problems. Ph.D. Thesis. 2006.
!	  Esbjerg Technical Institute, Aalborg University. Can be downloaded
!	  from http://vbn.aau.dk/fbspretrieve/14058639/JCthesis.pdf
!------------------------------------------------------------------------
! Johan Clausen
! Department of Civil Engineering
! Aalborg University
! April 2008
!------------------------------------------------------------------------
	implicit none

	real(8), intent(in) :: SigP_up(3), PlasPar(7), Dn(3,3), Dinvn(3,3)
	real(8), intent(in) :: Sipla(3), Gdir(3), Fdir(3)
	integer(4), intent(in) :: region

	real(8), intent(out) :: Tnorm(3,3), Tnorminv(3,3)
	!----------------------
	real(8) comp, m, s, a, mg, sg, ag, parf, parg, dellam, kg, dellamsek
	real(8), dimension(3,3) :: Hessian, HessianSek, Eye
	real(8), dimension(3) :: PlasStrain

!----- Explicit interfaces --------------------------------
	interface
		pure subroutine MatInv_2x2_3x3(A,n,Ainv)
			integer(4), intent(in) :: n
			real(8), dimension(n,n), intent(in) :: A

			real(8), dimension(n,n), intent(out) :: Ainv
		end subroutine MatInv_2x2_3x3
	end interface
!----------------------------------------------------------

!------------------------------------------------------------------------
	
!----- Initiation of variables -----------------------------
	comp = PlasPar(1) ! Uniaxial compressive strength
	m	 = PlasPar(2) ! Parameter which has a large saying in the amount of material friction
	s	 = PlasPar(3) ! Has a large saying in determination of tension strength
	a	 = PlasPar(4) ! Curvature yield parameter. For rock usually 0.5 <= a <= 0.65

	mg = PlasPar(5) ! The potential equivalent of the yield parameter m
	sg = PlasPar(6) ! The potential equivalent of the yield parameter s
	ag = PlasPar(7) ! Curvature parameter of plastic potential
!-----------------------------------------------------------
	Eye = reshape((/1.0_8, 0.0_8, 0.0_8,
     &				  0.0_8, 1.0_8, 0.0_8,
     &				  0.0_8, 0.0_8, 1.0_8/), (/3,3/)) ! Unit Matrix
      Tnorminv = Eye
      if (region < 4) then
		parg = sg - mg*SigP_up(1)/comp
		Hessian = 0.0_8
		Hessian(1,1) = (1-ag)*ag*mg**2/comp*parg**(ag-2)	
		PlasStrain = matmul(Dinvn,SiPla) ! Plastic strain increment
	
	
		if (region == 1) then ! return to the yield surface
			kg = Gdir(1) ! Plastic potential normal = [kg 0 -1]
			dellam = PlasStrain(1)/kg ! Plastic multiplier

			Tnorminv = Eye + dellam*matmul(Dn,Hessian)
	
		elseif (region == 2) then ! return to line 1, i.e. triaxial compression
			kg = Gdir(3) ! Gdir = [1 1 kg]' on line 1
			dellam = PlasStrain(1)/kg ! Plastic multiplier
			dellamsek = PlasStrain(2)/kg ! Plastic multiplier on secondary yield surface

			HessianSek = 0.0_8
			HessianSek(2,2) = Hessian(1,1)

			Tnorminv = Eye + dellam*matmul(Dn,Hessian) + dellamsek
     &											*matmul(Dn,HessianSek)

		elseif (region == 3) then ! return to line 2, i.e. triaxial extension
			dellam = -PlasStrain(3) ! Plastic multiplier
			dellamsek = -PlasStrain(2) ! Plastic multiplier on secondary yield surface
	
			HessianSek = 0.0_8
			HessianSek(1,1) = Hessian(1,1)

			Tnorminv = Eye + dellam*matmul(Dn,Hessian) + 
     &							dellamsek*matmul(Dn,HessianSek)
		end if ! region == ...
	end if ! (region < 4)
	call MatInv_2x2_3x3(Tnorminv,3,Tnorm) ! calculates the inverse of Tnorminv, which is Tnorm
	
!--------------------------------------------------------------------------------
	end subroutine TnormHoekBrownPrin !----------------------------------------
!--------------------------------------------------------------------------------

	pure subroutine FormModDepLine(mod_type,beta,SiPla,Fline,Gline,Dninv,
     1											Dc,Dcinv,nsigma,  Depc)
!--------------------------------------------------------------------------
! Calculates the modified elasto-plastic constitutive matrix at on a line
! for a perfectly plastic material, It is modified in the sence that is not
! the double singular matrix that theory prescribes. Only one type of
! modification is possible. The calculation are carried out in principal
! stress space.
!
! INPUT
! - Name -		-type and description --
!	mod_type	(integer,scalar) Detemines the modified type
!				  mod_type = 0: A Doubly singular matrix, as given in the
!							    References [1] and [2].
!				  mod_type = 1: A single singular matrix in the plastic
!							    strain direction, but almost singular in
!							    any direction perpendicular to the potential
!							    line. See Reference [3] for a detailed explanation.
!				  mod_type = 2: A single singular matrix in the Koiter direction,
!								i.e. the direction of the plastic strain.
!	beta		(real,scalar) The factor by which to divide the stiffness in
!				  the directions perpendicular to the line.
!	SiPla		(real,array(3)) The plastic corrector principal stresses.
!	Fline		(real,array(3)) Direction of the yield surface line in principal
!				  stress space.
!	Gline		(real,array(3)) Direction of the plastic potential line in principal
!				  stress space.
!	Dninv		(real,array(3,3)) Normal components of the elastic
!				  compliance matrix.
!	Dc			(real,array(nsigma,nsigma)) Modified elastic stiffness.
!				  If the yield criterion is linear Dc = D.
!	Dcinv		(real,array(nsigma,nsigma)) Inverse of the modified
!				  elastic stiffness. If the yield criterion is linear,
!				  Dcinv = Dinv.
!	nsigma		(integer, scalar) Number of stress components
!
! OUTPUT
!	Depc	(real,array(nsigma,nsigma)) Modified Elasto-plastic constitutive
!			  matrix on a line. If the criterion is linear Depc is the
!			  infinitesimal version, otherwise it is the consistent. In
!			  both cases it is expressed in principal stres space.
!
! REFERENCES:
! [1] Johan Clausen, Lars Damkilde and Lars Andersen: "Efficient return
!	  algorithms for associated plasticity with multiple yield planes", 
!	  International Journal for Numerical Methods in Engineering, 2006,
!	  volume 66, pages 1036-1059.
! [2] Johan Clausen, Lars Damkilde and Lars Andersen: An efficient return
!	  algorithm for non-associated plasticity with linear yield criteria
!	  in principal stress space. Computers & Structures, 2007, volume 85,
!	  pages 1795-1807.
! [3] Johan Clausen: Efficient Non-Linear Finite Element Implementation of
!	  Elasto-Plasticity for Geotechnical Problems. Ph.D. Thesis. 2006.
!	  Esbjerg Technical Institute, Aalborg University. Can be downloaded
!	  from http://vbn.aau.dk/fbspretrieve/14058639/JCthesis.pdf
!
!------------------------------------------------------------------------
! Johan Clausen
! Department of Civil Engineering
! Aalborg University
! January 2009
!------------------------------------------------------------------------
	implicit none
	integer, parameter :: dp = kind(1.0D0) ! Double precision kind

	!----- INOUT/OUTPUT ------
	integer, intent(in) :: mod_type, nsigma
	real(dp), intent(in) :: Dninv(3,3),Dc(nsigma,nsigma),
     1		Dcinv(nsigma,nsigma),beta,Sipla(3),Fline(3),Gline(3)

	real(dp), intent(out) :: Depc(nsigma,nsigma)
	!-------------------------

	real(dp) KoitDir(3),KoitPerDir(3), Dper(3,3)
	
	interface
		pure subroutine FormDepLinePerfect(D,Dinv,Ra,Rb,nsigma,DepLine)
			integer(4), intent(in) :: nsigma
			real(8), intent(in), dimension(nsigma,nsigma) :: D, Dinv
			real(8), intent(in), dimension(3) :: Ra, Rb
			real(8), intent(out) :: DepLine(nsigma,nsigma)
		end subroutine FormDepLinePerfect
		
		pure function Cross(A,B)
			real(8), intent(in)	::	A(3), B(3)
			real(8)				::  Cross(3)
		end function Cross
		
		pure subroutine FormDepPerfect(D,Norm,Edir,nsigma,Dep)
			integer(4), intent(in) :: nsigma
			real(8), intent(in) :: D(nsigma,nsigma), Norm(nsigma)
			real(8), intent(in) :: Edir(nsigma)
			real(8), intent(out) :: Dep(nsigma,nsigma)
		end subroutine FormDepPerfect
	end interface
!------------------------------------------------------------------
	Depc = 0.0_dp

	if (mod_type == 0) then
		call FormDepLinePerfect(Dc,Dcinv,Fline,Gline,nsigma,Depc)
	
	elseif (mod_type == 1) then
		KoitDir = matmul(Dninv,SiPla) ! The plastic strain direction in principal stress space
		KoitPerDir = Cross(KoitDir,Gline) ! Direction perpendicular to the strain direction
											 ! and the plastic potential line
	
		call FormDepLinePerfect(Dc,Dcinv,Fline,Gline,nsigma,Depc)
		call FormDepLinePerfect(Dc(1:3,1:3),Dcinv(1:3,1:3),KoitPerDir,
     1												KoitPerDir,3,Dper)
		Depc(1:3,1:3) = Depc(1:3,1:3) + Dper/beta
	elseif (mod_type == 2) then
		KoitDir = matmul(Dninv,SiPla) ! The plastic strain direction in principal stress space
		call FormDepPerfect(Dc(1:3,1:3),KoitDir,KoitDir,3,Depc(1:3,1:3))
		Depc(4:nsigma,4:nsigma) = Dc(4:nsigma,4:nsigma)
	end if
!--------------------------------------------------------------------------------
	end subroutine FormModDepLine !--------------------------------------------
!--------------------------------------------------------------------------------

	pure subroutine FormModDepApex(mod_type,alpha,SiPla,Dninv,Dc,Dcinv,
     1												nsigma,ishydro,  Depc)
!--------------------------------------------------------------------------
! Calculates the modified elasto-plastic constitutive matrix at the apex
! for a perfectly plastic material, It is modified in the sence that is not
! the zero matrix. Two different modifications are possible, and they are
! determined by the variable mod_type. The calculations are carried out in
! principal stress space.
!
! INPUT
! - Name -		-type and description --
!	mod_type	(integer,scalar) Detemines the modified type
!				  mod_type= 0: Depc is the zero matrix. See Refernce [1]
!				  mod_type= 1: A single singular matrix in the plastic
!							   strain direction. See Reference [2]
!				  mod_type= 2: A double singular matrix in the plastic
!							   strain direction and in the direction of
!							   the hydrostatic stress line.
!	alpha		(real,scalar) The facto by which to divide the stiffness
!	SiPla		(real,array(3)) The plastic corrector principal stresses.
!	Dninv		(real,arrey(3,3)) Normal components of the elastic
!				  compliance matrix.
!	Dc			(real,array(nsigma,nsigma)) Modified elastic stiffness.
!				  If the yield criterion is linear Dc = D.
!	Dcinv		(real,array(nsigma,nsigma)) Inverse of the modified
!				  elastic stiffness. If the yield criterion is linear,
!				  Dcinv = Dinv.
!	nsigma		(integer, scalar) Number of stress components
!	ishydro		(integer, scalar). Determines whether or not the point is located
!				  on the hydrostatic axis.
!				  ishydro = 0: The point is not on the hydrostatic axis.
!				  ishydro = 1: The point is on the hydrostatic axis.
!
! OUTPUT
!	Depc	(real,array(nsigma,nsigma) Modified Elasto-plastic constitutive
!			  matrix on an apex. If the criterion is linear Depc is the
!			  infinitesimal version, otherwise it is the consistent. In
!			  both cases it is expressed in principal stres space.
!
! REFERENCES:
! [1] Johan Clausen, Lars Damkilde and Lars Andersen: "Efficient return
!	  algorithms for associated plasticity with multiple yield planes", 
!	  International Journal for Numerical Methods in Engineering, 2006,
!	  volume 66, pages 1036-1059.
! [2] Johan Clausen: Efficient Non-Linear Finite Element Implementation of
!	  Elasto-Plasticity for Geotechnical Problems. Ph.D. Thesis. 2006.
!	  Esbjerg Technical Institute, Aalborg University. Can be downloaded
!	  from http://vbn.aau.dk/fbspretrieve/14058639/JCthesis.pdf
!------------------------------------------------------------------------
! Johan Clausen
! Department of Civil Engineering
! Aalborg University
! January 2009
!------------------------------------------------------------------------
	implicit none
	integer, parameter :: dp = kind(1.0D0) ! Double precision kind

	!----- INPUT/OUTPUT ------
	integer, intent(in) :: mod_type, nsigma, ishydro
	real(dp), intent(in) :: alpha,Sipla(3),Dninv(3,3),Dc(nsigma,nsigma),
     1						  Dcinv(nsigma,nsigma)

	real(dp), intent(out) :: Depc(nsigma,nsigma)
	!-------------------------

	real(dp) KoitDir(3), Dnorm(3,3), Pdir(3), PerDir(3)
	interface
		pure subroutine FormDepLinePerfect(D,Dinv,Ra,Rb,nsigma,DepLine)
		  integer(4), intent(in) :: nsigma
		  real(8), intent(in), dimension(nsigma,nsigma) :: D, Dinv
		  real(8), intent(in), dimension(3) :: Ra, Rb
		  real(8), intent(out) :: DepLine(nsigma,nsigma)
		end subroutine FormDepLinePerfect
		!---
		pure function Cross(A,B)
		  real(8), intent(in)	::	A(3), B(3)
		  real(8)				::  Cross(3)
		end function Cross
		!---
		pure subroutine FormDepPerfect(D,Norm,Edir,nsigma,Dep)
			integer(4), intent(in) :: nsigma
			real(8), intent(in) :: D(nsigma,nsigma), Norm(nsigma)
			real(8), intent(in) :: Edir(nsigma)
			real(8), intent(out) :: Dep(nsigma,nsigma)
		end subroutine FormDepPerfect
	end interface
!--------------------------------------------------------------
	Depc = 0.0_dp

	KoitDir = matmul(Dninv,SiPla)

	if (mod_type == 0) then ! Zero matrix at the apex
		if (ishydro == 0) then ! the stress point is not located on the hydrostatic line
			Depc(4:nsigma,4:nsigma) = Dc(4:nsigma,4:nsigma) ! Shear part
		end if	
	elseif (mod_type == 1) then ! Single singular matrix (DepKo/fak)
		call FormDepPerfect(Dc(1:3,1:3),KoitDir,KoitDir,3,Dnorm)
		Depc(1:3,1:3) = Dnorm/alpha ! Normal part
		Depc(4:nsigma,4:nsigma) = Dc(4:nsigma,4:nsigma)/alpha ! Shear part
	elseif (mod_type == 2) then ! Double singular matrix (KoiteP)
		Pdir = (/1.0_8, 1.0_8, 1.0_8/) ! Hydrostatic direction
		PerDir = Cross(KoitDir,Pdir)
	
		call FormDepLinePerfect(Dc,Dcinv,Perdir,PerDir,nsigma,Depc)
		Depc = Depc/alpha
	end if
!--------------------------------------------------------------------------------
	end subroutine FormModDepApex !--------------------------------------------
!--------------------------------------------------------------------------------

	pure subroutine TshearPrinPerfect(SigP,SigP_up,nshear,s1,s2,
     &												Tshear,Tshearinv)
!--------------------------------------------------------------------------
! Forms the shear part of the modification matrix T in principal stress space
! for at perfectly plastic material. T is used when calculation the modified
! elastic stiffness matrix Dc which is, in turn, used to form the consistent
! constitutive matrix. The shear part of T is unaffected whether the yield
! criterion is linear or not.
! 
!
! INPUT
!  - Name -	-type,size -	-- Description --
!	SigP	 real,ar(3)		Principal predictor stresses in descending order
!							  SigP = [sigP_1 sigP_2 sigP_3]
!	SigP_up	 real,ar(3)		Updated principal stresses in descending order
!	nshear	 int,sc			Number of shear stress components. 1 in plane problems
!							  and 3 in full 3D
!	s1		 int,sc			Position of largest principal in-plane stress in SigP.
!							  Only used in plane situations
!	s2		 int,sc			Position of smallest principal in-plane stress in SigP.
!							  Only used in plane situations
!
! OUTPUT
!	Tshear		real,ar(nshear,nshear) Shear part of the modification matrix T
!	Tshearinv	real,ar(nshear,nshear) Inverse of Tshear
!
!------------------------------------------------------------------------
! Johan Clausen
! Department of Civil Engineering
! Aalborg University
! April 2008
!------------------------------------------------------------------------
	implicit none
	real(8), intent(in) :: SigP(3), SigP_up(3)
	integer(4), intent(in) :: nshear, s1, s2

	real(8), dimension(nshear,nshear), intent(out) :: Tshear, Tshearinv
!------------------------------------------------------------------------

	Tshear = 0.0_8
	Tshearinv = 0.0_8

	if (nshear == 1) then
		if (SigP_up(s1) - SigP_up(s2) > 0.0_8) then
			Tshear = (SigP_up(s1)-SigP_up(s2)) / (SigP(s1)-SigP(s2))
			Tshearinv = 1.0_8/Tshear
		end if
	elseif (nshear == 3) then
		if (SigP_up(1) - SigP_up(2) > 0.0_8) then
			Tshear(1,1) = (SigP_up(1)-SigP_up(2)) / (SigP(1)-SigP(2))
			Tshearinv(1,1) = 1.0_8/Tshear(1,1)
		end if
		if (SigP_up(1) - SigP_up(3) > 0.0_8) then
			Tshear(2,2) = (SigP_up(1)-SigP_up(3)) / (SigP(1)-SigP(3))
			Tshearinv(2,2) = 1.0_8/Tshear(2,2)
		end if
		if (SigP_up(2) - SigP_up(3) > 0.0_8) then
			Tshear(3,3) = (SigP_up(2)-SigP_up(3)) / (SigP(2)-SigP(3))
			Tshearinv(3,3) = 1.0_8/Tshear(3,3)
		end if
	end if

!--------------------------------------------------------------------------------
	end subroutine TshearPrinPerfect !-----------------------------------------
!--------------------------------------------------------------------------------


	pure subroutine PrinStressAna(Sigma,nsigma,SigP)
!--------------------------------------------------------------------------
! Calculates the principal stress of a stress vector Sigma. Analytical 
! expressions are used in the calculations.
!
! INPUT
!  - Name - -type,size -  -- description --
!     Sigma (real,ar(nsigma)Stress vector. Size is dependent on the type of stress
!                   nsigma = 4 in plane situations and nsigma = 6 in 3D
!                   The ordering of the components must be as follows:
!                   Plane situations (plane stress, plane strain and axisymmetry):
!                   Sigma = [sig_x sig_y sig_z tau_xy] (i.e. sig_z is out-of-plane)
!                   General 3D:  Sigma = [sig_x sig_y sig_z tau_xy tau_xz tau_yz]
!     nsigma  (integer,sc)  Number of stress components, i.e. the size of Sigma
!
! OUTPUT
!     SigP  (real,ar(3))  Principal stresses in descending order
!
!---------------------------------------------------------------------------
! Johan Clausen
! Esbjerg Institute of Technology
! Aalborg University
! November 2005
!--------------------------------------------------------------------------
	implicit none
      integer(4), intent(in) :: nsigma
      real(8), intent(in) :: Sigma(nsigma)

      real(8), intent(out) :: SigP(3)

!----- Internal variables -------------------------------------------------
      real(8) sig_av, sig_hj, I1, J2, J3, lode, sigm
      real(8) sqJ2, sin3lode
      real(8) S(6)

!----- Explicit interfaces ------------------------------------------
	interface
		pure subroutine Invariants(Sigma,nsigma,S,I1,J2,J3,lode,sin3lode)
			integer(4), intent(in) :: nsigma
			real(8), intent(in) :: Sigma(nsigma)
	
			real(8), intent(out) :: S(nsigma), I1, J2
			real(8), intent(out) :: J3, lode, sin3lode
		end subroutine Invariants
	end interface
!--------------------------------------------------------------------

!----- Initialization of SigP ---------------------------------------------
      SigP = 0.0_8

!----- Plane situations including axisymmetry -----------------------------

      if (nsigma == 4) then
        
        sig_av = 0.5_8*(Sigma(1) + Sigma(2)) ;
        sig_hj = sqrt((0.5_8*(Sigma(1) - Sigma(2)))**2 + Sigma(4)**2)
    
        SigP(1) = sig_av + sig_hj
        SigP(2) = sig_av - sig_hj
        SigP(3) = Sigma(3) ! The out-of-plane stress

      !-- Sorting: SigP(1) >= SigP(2) >= SigP(3) --
        if (Sigma(3) > SigP(1)) then
          SigP(3) = SigP(2)  ;  SigP(2) = SigP(1)  ; SigP(1) = Sigma(3)
        elseif (Sigma(3) > SigP(2) ) then
          SigP(3) = SigP(2)  ;  SigP(2) = Sigma(3) 
        end if
      !--------------------------------------------
  
!------------------------------------------------------------------------

!----- General stress state ---------------------------------------------
      elseif (nsigma == 6) then
      
!   ----- Invariants ----------------------------------------------------
		call Invariants(Sigma,nsigma,S,I1,J2,J3,lode,sin3lode)
		sigm = I1/3.0_8 ! Hydrostatic stress

		if (sqrt(J2) < 1.0D-12*abs(I1)) then ! A hydrostatic stress state
			SigP(1) = Sigma(1)
			SigP(2) = Sigma(2)
			SigP(3) = Sigma(3)
		else
			
!   ----- Principal stresses -----------------------------------------
			sqJ2 = 1.154700538379252_8*sqrt(J2)
      !				2/sqrt(3)
      !                                    2*pi/3
			SigP(1) = sqJ2 * dsin(lode + 2.094395102393195_8) + sigm
			SigP(2) = sqJ2 * dsin(lode)                       + sigm
			SigP(3) = sqJ2 * dsin(lode - 2.094395102393195_8) + sigm
		endif
      end if ! nsigma == ...

!----------------------------------------------------------------------------------
      end subroutine PrinStressAna !-----------------------------------------------
!----------------------------------------------------------------------------------


      pure subroutine FormDepPerfect(D,Norm,Edir,nsigma,Dep)
!--------------------------------------------------------------------------
! Calculates elasto-plastic constitutive matrix of a material with 
! non-associated and perfect plasticity. If Norm = Edir the flow rule
! is associated
!
! INPUT
!     D   (real) elastic constitutive matrix. Size = nsigma x nsigma
!     Norm  (real) Gradient (normal) of the yield surface. Size = nsigma
!     Edir  (real) Gradient (normal) of the plastic potential. Size = nsigma
!     nsigma  (integer) Number of stress components
!
! OUTPUT
!     Dep   (real) Elasto-plastic infinitesimal constitutive matrix
!
!--------------------------------------------------------------------------
! Johan Clausen
! Esbjerg Department of Engineering
! Aalborg University
! July 2005
!--------------------------------------------------------------------------
	implicit none
      integer(4), intent(in) :: nsigma
	real(8), intent(in) :: D(nsigma,nsigma), Norm(nsigma)
	real(8), intent(in) :: Edir(nsigma)
      real(8), intent(out) :: Dep(nsigma,nsigma)
!--------------------------------------------------------------------------
      real(8), dimension(nsigma) :: Num1, Num2
      real(8) den, Num(nsigma,nsigma)
      integer(4) i,j
!--------------------------------------------------------------------------
      Num1 = matmul(D,Edir)
      Num2 = matmul(Norm,D)
      forall(i = 1:nsigma, j = 1:nsigma) Num(i,j) = Num1(i)*Num2(j)
      den  = dot_product(Norm,Num1)
      Dep = D - Num/den
!--------------------------------------------------------------------------
      end subroutine FormDepPerfect !--------------------------------------
!--------------------------------------------------------------------------

	pure subroutine FormDepLinePerfect(D,Dinv,Ra,Rb,nsigma,DepLine)
!--------------------------------------------------------------------------
! Calculates double singular elasto-plastic constitutive matrix of a
! material with non-associated and perfect plasticity. If Ra = Rb the flow
! rule is associated. The formula is only valid in principal stress space
!
! INPUT
!     D   (real) elastic constitutive matrix. Size = nsigma x nsigma
!     Dinv  (real) Inverted elastic constitutive matrix. Size = nsigma x nsigma
!     Ra    (real) Direction of the line on the yield surface in principal
!           stress space. If the line is curved Ra is the tangent
!     Rb    (real) Direction of the line on the plastic potential in principal
!           stress space. If the line is curved Rb is the tangent
!
! OUTPUT
!     Depline (real) Double singular elasto-plastic infinitesimal
!           constitutive matrix on a line. Size = nsigma x nsigma
!
!------------------------------------------------------------------------
! Johan Clausen
! Esbjerg Department of Engineering
! Aalborg University
! July 2005
!------------------------------------------------------------------------
	implicit none
	!----- INPUT/OUTPUT -----
      integer(4), intent(in) :: nsigma
      real(8), intent(in), dimension(nsigma,nsigma) :: D, Dinv
      real(8), intent(in), dimension(3) :: Ra, Rb

      real(8), intent(out) :: DepLine(nsigma,nsigma)
      !-------------------------

      real(8) den, Den1(3), Num(3,3), Dprin(3,3)
      integer(4) i,j
!-----------------------------------------------------------------

      forall(i = 1:3, j = 1:3) Num(i,j) = Ra(i)*Rb(j)

      Den1 = matmul(Dinv(1:3,1:3),Rb)
      den = dot_product(Ra,Den1) ! denomimator, Ra'*Dinv*Rb

      Dprin = Num/den

      DepLine = 0.0_8
      DepLine(1:3,1:3) = Dprin
      if (nsigma > 3) then
		DepLine(4,4) = D(4,4) ! In plane situations nsigma == 4
		if (nsigma == 6) then ! General three dimensional stress state 
			DepLine(5,5) = D(5,5)
			DepLine(6,6) = D(6,6)
		end if
	end if

!-------------------------------------------------------------------------------
      end subroutine FormDepLinePerfect !---------------------------------------
!-------------------------------------------------------------------------------

				
	pure subroutine PrinDirect(Sigma,nsigma,SigP,psi)
!--------------------------------------------------------------------------
! Calculates the principal directions of a principal stress state
! Analytical expressions are used in the calculations.
!
! INPUT
!  - Name -	-type,size -	-- description --
!	Sigma	(real,ar(nsigma)Stress vector. Size is dependent on the type of stress
!							  nsigma = 4 in plane situations and nsigma = 6 in 3D
!							  The ordering of the components must be as follows:
!							  Plane situations (plane stress, plane strain and axisymmetry):
!							  Sigma = [sig_x sig_y sig_z tau_xy] (i.e. sig_z is out-of-plane)
!							  General 3D:  Sigma = [sig_x sig_y sig_z tau_xy tau_xz tau_yz]
!	nsigma	(integer,sc)	Number of stress components, i.e. the size of Sigma
!	SigP	(real,ar(3))	Principal stresses in descending order
!
! OUTPUT
!	psi		(real,ar(3,3))	Principal angle if nsigma = 4 or
!							  vector of normalized eigenvectors (principal directions)
!							  if nsigma = 6. When nsigma = 4 the angle is stored
!							  in the upper left corner, i.e. psi(1,1).
!
!---------------------------------------------------------------------------
! Johan Clausen
! Esbjerg Institute of Technology
! Aalborg University
! November 2005
!--------------------------------------------------------------------------
	integer(4), intent(in) :: nsigma
	real(8), intent(in) :: Sigma(nsigma), SigP(3)

	real(8), intent(out)  :: psi(3,3)

!----- Internal variables -------------------------------------------------
	real(8), parameter :: pi = 3.14159265358979323_8
	real(8) tol, nl(3), normal(3,3), leng, tol1, tol2
	integer(4) i, hj1, hj2, nr

	interface
		pure function Cross(A,B)
			real(8), intent(in)	::	A(3), B(3)
			real(8)				::  Cross(3)
		end function Cross
	end interface
!----- Plane situations including axisymmetry -----------------------------
	if (nsigma == 4) then
      	if (Sigma(1) > Sigma(2)  .and.  Sigma(4) >= 0.0_8) then
			psi = 0.5_8*atan(2*Sigma(4)/(Sigma(1)-Sigma(2)))
		elseif (Sigma(1) < Sigma(2)  .and.  Sigma(4) >= 0.0_8) then
			psi = 0.5_8*(pi - atan(2*Sigma(4)/(Sigma(2)-Sigma(1))))
		elseif (Sigma(1) < Sigma(2)  .and.  Sigma(4) < 0.0_8) then
			psi = 0.5_8*(atan(-2*Sigma(4)/(Sigma(2)-Sigma(1))) + pi)
		elseif (Sigma(1) > Sigma(2)  .and.  Sigma(4) < 0.0_8) then
			psi = 0.5_8*(2*pi - atan(-2*Sigma(4)/(Sigma(1)-Sigma(2))))
		elseif (Sigma(1) == Sigma(2) .and. Sigma(4) > 0.0_8) then
			psi = 0.25_8*pi
		elseif (Sigma(1) == Sigma(2) .and. Sigma(4) < 0.0_8) then
			psi = 0.75_8*pi
		elseif (Sigma(1) == Sigma(2) .and. Sigma(4) == 0.0_8) then
			psi = 0.0_8 ;
		end if

!------------------------------------------------------------------------


!----- General 3D ----------------------------------------------------
	elseif (nsigma == 6) then
		tol1 = abs(1.0e-10_8*(SigP(1) - SigP(3))) ! tolerance value 1
		tol2 = abs(1.0e-12_8*(SigP(1) + SigP(2) + SigP(3))) ! tolerance
		tol = maxval((/tol1,tol2,1.0e-13_8/))
		
		nl = 0.0_8
		normal = 0.0_8

!	-- Sigma is already expressed in principal coordinates ---
		if (abs(Sigma(4)) < tol  .and.
     &		abs(Sigma(5)) < tol  .and.
     &		abs(Sigma(6)) < tol)		then
     			if (Sigma(1) >= Sigma(2) .and. Sigma(1) >= Sigma(3)) then
				normal(1,1) = 1.0_8
				if (Sigma(2) >= Sigma(3)) then
					normal(2,2) = 1.0_8 ; normal(3,3) = 1.0_8
				else
					normal(3,2) = 1.0_8 ; normal(2,3) = 1.0_8
				end if
			elseif (Sigma(2) >= Sigma(1) .and. Sigma(2) >= Sigma(3))
     &															 then
     				normal(2,1) = 1.0_8
				if (Sigma(1) >= Sigma(3)) then
					normal(1,2) = 1.0_8 ; normal(3,3) = 1.0_8
				else
					normal(3,2) = 1.0_8 ; normal(1,3) = 1.0_8
				end if
			elseif (Sigma(3) >= Sigma(1) .and. Sigma(3) >= Sigma(2))
     &															 then
     				normal(3,1) = 1.0_8
				if (Sigma(1) >= Sigma(2)) then
					normal(1,2) = 1.0_8 ; normal(2,3) = 1.0_8
				else
					normal(2,2) = 1.0_8 ; normal(1,3) = 1.0_8
				end if
			end if
!	-- Three distinct eigenvalues -----
		elseif (abs(SigP(1)-SigP(2)) > tol  .and.
     &			abs(SigP(2)-SigP(3)) > tol  .and.
     &			abs(SigP(1)-SigP(3)) > tol)			then
     			do i = 1,2
				nl(1) = (Sigma(2)-SigP(i))
     &						* (Sigma(3)-SigP(i)) - Sigma(6)**2
				nl(2) = -Sigma(4) * (Sigma(3)-SigP(i))
     &						+ Sigma(6) * Sigma(5)
				nl(3) = Sigma(4) * Sigma(6)
     &						- (Sigma(2)-SigP(i)) * Sigma(5)
				leng = sqrt(nl(1)**2 + nl(2)**2 + nl(3)**2)
				if (leng < tol) then
					normal(i,i) = 1 ;
				else
					normal(:,i) = nl/leng;
				end if
			end do
			normal(:,3) = Cross(Normal(:,1),Normal(:,2))
			
!	-- Three equal eigenvalues ----
		elseif (abs(SigP(1)-SigP(2)) < tol  .and.
     &			abs(SigP(2)-SigP(3)) < tol  .and.
     &			abs(SigP(1)-SigP(3)) < tol)			then
     			
			normal(1,1) = 1.0_8
			normal(2,2) = 1.0_8
			normal(3,3) = 1.0_8

!	-- two equal eigenvalues -----
		else
			if (abs(SigP(1) - SigP(2)) <= tol) then
				hj1 = 1 ; hj2 = 2 ; nr = 3
			elseif (abs(SigP(1) - SigP(3)) <= tol) then
				hj1 = 1 ; nr = 2 ; hj2 = 3
			else 
				nr = 1 ; hj1 = 2 ; hj2 = 3
			end if
            
!       --- First principal direction -------------------------            
			nl(1) = (Sigma(2)-SigP(nr)) * (Sigma(3)-SigP(nr))
     &				- Sigma(6)**2
			nl(2) = -Sigma(4) * (Sigma(3)-SigP(nr))
     &				+ Sigma(6) * Sigma(5)
			nl(3) = Sigma(4) * Sigma(6)
     &				- (Sigma(2)-SigP(nr)) * Sigma(5)
			leng = sqrt(nl(1)**2 + nl(2)**2 + nl(3)**2)
			normal(:,nr) = nl/leng ;

!       --- Second principal direction ------------------------            
			if (normal(3,nr) > 0.0_8  .or. normal(3,nr) < 0.0_8) then
				nl(1) = 1.0_8 ; nl(2) = 1.0_8 ;
				nl(3) = -(normal(1,nr)+normal(2,nr))/normal(3,nr) ;
			elseif (normal(2,nr) > 0.0_8  .or.
     &								normal(2,nr) < 0.0_8) then
				nl(1) = 1.0_8 ; nl(3) = 1.0_8 ;
				nl(2) = -(normal(1,nr)+normal(3,nr))/normal(2,nr) ;
			elseif (normal(1,nr) > 0.0_8  .or.
     &								normal(1,nr) < 0.0_8) then
				nl(2) = 1.0_8 ; nl(3) = 1.0_8
				nl(1) = -(normal(2,nr)+normal(3,nr))/normal(1,nr) ;             
			end if
			leng = sqrt(nl(1)**2 + nl(2)**2 + nl(3)**2)
			normal(:,hj1) = nl/leng ;
!       --- Third principal direction -------------------------
			normal(:,hj2) = Cross(normal(:,nr),Normal(:,hj1))
		end if
		psi = normal
	end if ! nsigma == ....

!----------------------------------------------------------------------------------
	end subroutine PrinDirect !--------------------------------------------------
!----------------------------------------------------------------------------------



      pure subroutine TransMatrix(psi,nsigma,ouplP,A)
!--------------------------------------------------------------------
! Calculates transformation matrix, A, for stress transformation
! based on either an angle (plane and axisymmetry) or direction
! cosines (general 3D). The ordering of the shear stresses
! must be specified, as it influences how A is built
!
! INPUT
!  - Name- -type,size -   -- Description --
!     psi   (real,ar(3,3))  Matrix of direction cosines in general
!                   three-dimensional stress state or the rotation
!                   angle otherwise. In 3D we have 
!                   sig_tranformed = psi'*sig*psi, where sig
!                   is a stress tensor
!     nsigma  (integer,sc)  Number of stress components, i.e. size of A
!     ouplP (integer,sc)  Only used in plan situations in which it is similar to
!                   hoop but for the principal stress vector SigP, i.e.
!                   it is the position of the out-of-plane principal stress.
!
! OUTPUT
!     A Transformation matrix, with elements ordered according to hoop and ouplP
!         A stress is transformed according to ("'" signifies matrix transpose)
!         Sigma_transformed = inv(A')*Sigma and conversely
!         Sigma = A'*Sigma_transformed
!         A strain vector, Epsilon, is tranformed according to
!         Epsilon_transformed = A*Epsilon and conversely
!         Epsilon = inv(A)*Epsilon_transformed
!         See e.g. [Cook, Malkus and Plesha, 1989] for further details
!--------------------------------------------------------------------
! Johan Clausen
! Esbjerg Department of Engineering
! Aalborg University
! July 2005
!--------------------------------------------------------------------

      !-----------------
      integer(4), intent(in) :: nsigma, ouplP
      real(8), intent(in)    :: psi(3,3)

      real(8), intent(out)   :: A(nsigma,nsigma)
      !-----------------

      real(8) sin_psi, cos_psi, sin_psi2, cos_psi2, sin_2psi
      real(8), dimension(3,3) :: A1, A2, A3, A4, Asmall
      real(8), dimension(3,4) :: Aint
      integer(4), dimension(3):: Hj1, Hj2
      integer(4) Ext(3,4), ExtP(4,3)
      integer(4) i, j

!---- Plane situations --------------------------------------------
! Assumes that the stress components are ordered according to
! Sigma = [sig_x sig_y sig_z tau_xy]
!------------------------------------------------------------------ 
      if (nsigma == 4) then
        sin_psi  = dsin(psi(1,1)) ; cos_psi = dcos(psi(1,1))
        sin_psi2 = sin_psi**2     ; cos_psi2 = cos_psi**2
        sin_2psi = dsin(2.0_8*psi(1,1))
      
      
        Asmall(1,1) = cos_psi2  ; Asmall(1,2) = sin_psi2
        Asmall(2,1) = sin_psi2  ; Asmall(2,2) = cos_psi2
        Asmall(3,1) = -sin_2psi ; Asmall(3,2) = sin_2psi
        
        Asmall(1,3) =  0.5_8*sin_2psi
        Asmall(2,3) = -0.5_8*sin_2psi
        Asmall(3,3) = cos_psi2 - sin_psi2

        ExtP = 0.0_8
        if (ouplP == 2) then
          ExtP(1,1) = 1  ;  ExtP(3,2) = 1  ;  ExtP(4,3) = 1
        elseif (ouplP == 1) then
          ExtP(2,1) = 1  ;  ExtP(3,2) = 1  ;  ExtP(4,3) = 1
        elseif (ouplP == 3) then
          ExtP(1,1) = 1  ;  ExtP(2,2) = 1  ;  ExtP(4,3) = 1
        end if

				  ! IMPORTANT! This order of
        Ext = 0.0_8 ! Ext assumes that the out-of-plane stress i in position 3!!
        Ext(1,1) = 1 ; Ext(2,2) = 1 ; Ext(3,4) = 1

        Aint = matmul(Asmall,Ext)
        A = matmul(ExtP,Aint) ! Transformation matrix
        A(ouplP,3) = 1.0_8

!---- General three dimensional stress state ------------------------
! Assumes that the stress components are ordered according to
! Sigma = [sig_x sig_y sig_z tau_xy tau_xz tau_yz]
!------------------------------------------------------------------ 
      elseif (nsigma == 6) then
      
        Hj1 = (/1, 3, 2/)
        Hj2 = (/2, 1, 3/)
      
        A1 = psi**2
      
        do i = 1,3
          do j = 1,3
            A2(i,j) = psi(i,Hj1(j))*psi(i,Hj2(j))
            A3(i,j) = psi(Hj1(i),j)*psi(Hj2(i),j)
            A4(i,j) = psi(Hj1(i),Hj1(j))*psi(Hj2(i),Hj2(j))
     &              + psi(Hj2(i),Hj1(j))*psi(Hj1(i),Hj2(j))
          end do
        end do
        A2 = 2*A2
        
        A(1:3,1:3) = A1 ; A(1:3,4:6) = A2
        A(4:6,1:3) = A3 ; A(4:6,4:6) = A4
		A = transpose(A)
      end if ! nsigma == ....

!----------------------------------------------------------------------------------
      end subroutine TransMatrix !-------------------------------------------------
!----------------------------------------------------------------------------------

      pure subroutine DlinElas(E,nu,nsigma,D,Dinv)
!--------------------------------------------------------------------------
! Calculates linear elastic stiffness matrix D. If nsigma = 3, i.e.
! only three stress components plane stress is assumed
!
! INPUT
!     E   (real) Youngs module, module of elasticity
!     nu    (real) Poisson's ratio
!     nsigma  (int)  Number of stress components. 
!             nsigma = 4 in plane problems in general
!             if nsigma = 3 plane stress is assumed
!             nsigma = 6 in general stress states
! OUTPUT
!     D   (real) Linear elastic stiffness matrix D
!           Size is dependent on stress type
!     Dinv  (real) Inverse of D, optional
!------------------------------------------------------------------------
! Johan Clausen
! Esbjerg Department of Engineering
! Aalborg University
! July 2005
!------------------------------------------------------------------------
	implicit none
	!--- input/output -----
      integer(4), intent(in) :: nsigma
      real(8), intent(in) :: E, nu

      real(8), intent(out) :: D(nsigma,nsigma)
      real(8), intent(out), optional :: Dinv(nsigma,nsigma)
	!---------------------

      real(8) c, g
!-------------------------------------------------------------
      D = 0.0_8

!----- Plane stress with three stress components -----------------------
      if (nsigma == 3) then
      
        D(1,1) = 1.0_8 ; D(1,2) = nu
        D(2,1) = nu    ; D(2,2) = 1.0_8
                         D(3,3) = (1.0_8-nu)/2
        D = E/(1.0_8-nu*nu) * D

        if (present(Dinv)) then
          Dinv = 0.0_8
          Dinv(1,1) = 1.0_8 ; Dinv(1,2) = -nu
          Dinv(2,1) = -nu   ; Dinv(2,2) = 1.0_8
                              Dinv(3,3) = 2.0_8*(1+nu)
          Dinv = Dinv/E
        end if

!----- Plane stress, strain or axisymmetry (four stress components) ---
      elseif(nsigma == 4) then
		D(1,1) = 1.0_8 - nu ; D(1,2) = nu         ; D(1,3) = nu
		D(2,1) = nu         ; D(2,2) = 1.0_8 - nu ; D(2,3) = nu
		D(3,1) = nu         ; D(3,2) = nu         ; D(3,3) = 1.0_8- nu
		D(4,4) = (1.0_8-2.0_8*nu)/2
        D = E/((1+nu)*(1.0_8-2.0_8*nu)) * D

        if (present(Dinv)) then
          Dinv = 0.0_8
          Dinv(1,1) = 1.0_8 ; Dinv(1,2) = -nu     ; Dinv(1,3) = -nu
          Dinv(2,1) = -nu   ; Dinv(2,2) = 1.0_8   ; Dinv(2,3) = -nu
          Dinv(3,1) = -nu   ; Dinv(3,2) = -nu     ; Dinv(3,3) = 1.0_8
          Dinv(4,4) = 2.0_8*(1.0_8+nu)
          Dinv = Dinv/E
        end if
!----- General stress state ---------------------------------------------
      elseif(nsigma == 6) then
        c = E/((1.0_8+nu)*(1.0_8-2.0_8*nu))
        g = E/(2.0_8*(1.0_8+nu))

        D(1,1) = (1.0_8-nu)*c ; D(1,2) = nu*c   ; D(1,3) = D(1,2)
        D(2,1) = D(1,2)   ; D(2,2) = D(1,1) ; D(2,3) = D(1,2)
        D(3,1) = D(1,2)   ; D(3,2) = D(1,2) ; D(3,3) = D(1,1)
        D(4,4) = g
        D(5,5) = g
        D(6,6) = g

        if (present(Dinv)) then
          Dinv = 0.0_8
          Dinv(1,1) = 1.0_8 ; Dinv(1,2) = -nu   ; Dinv(1,3) = -nu
          Dinv(2,1) = -nu   ; Dinv(2,2) = 1.0_8 ; Dinv(2,3) = -nu
          Dinv(3,1) = -nu   ; Dinv(3,2) = -nu   ; Dinv(3,3) = 1.0_8
          Dinv(4,4) = 2.0_8*(1.0_8+nu)
          Dinv(5,5) = 2.0_8*(1.0_8+nu)
          Dinv(6,6) = 2.0_8*(1.0_8+nu)
          Dinv = Dinv/E
        end if
      end if ! nsigma = ...
!--------------------------------------------------------------------------------
      end subroutine  DlinElas !-------------------------------------------------
!--------------------------------------------------------------------------------


	pure subroutine Invariants(Sigma,nsigma,S,I1,J2,J3,lode,sin3lode)
!--------------------------------------------------------------------------
! Calculates the deviator stress vector of a stress vector Sigma, along
! with several stress invariants.
!
! INPUT
!  - Name -	-- description --
!	Sigma	(real,ar(nsigma): Stress vector. Size is dependent on the type of stress
!			  as given by nsigma.
!			  nsigme = 3: It is assumed that Sigma is principal stresses,
!						  Sigma = [SigP_1 SigP_2 SigP_3]
!			  nsigma = 4: Plane situations, i.e. plane stress, plane strain
!						  and axisymmetry. The ordering of the stresses must
!						  follow Sigma = [sig_x sig_y sig_z tau_xy] , i.e.
!						  sig_z is the out-of-plane stress.
!			  nsigma = 6  General 3D. The stresses must be ordered according
!						  to Sigma = [sig_x sig_y sig_z tau_xy tau_xz tau_yz]
!	nsigma	(integer,sc): Number of stress components, i.e. the size of Sigma.
!
! OUTPUT
!	S		(real, ar(nsigma)): Deviator stress vector. Same size as Sigma
!	I1		(real,scalar): First stress invariant.
!	J2		(real,scalar): Second deviator stress invariant.
!	J3		(real,scalar): Third deviator stress invariant.
!	lode	(real,scalar): Lode angle in radians. In the 
!			  range -pi/6 <= lode <= pi/6,    (pi/6 = 30 deg)
!	sin3lode (real,scalar): sin(3*lode)
!
!---------------------------------------------------------------------------
! Johan Clausen
! Deparment of Civil Engineering
! Aalborg University
! May 2008
!--------------------------------------------------------------------------
	implicit none
	integer(4), intent(in) :: nsigma
	real(8), intent(in) :: Sigma(nsigma)

	real(8), intent(out) :: S(nsigma), I1, J2
	real(8), intent(out) :: J3, lode, sin3lode

!--------------------------------------------------------------------------

	I1 = sum(Sigma(1:3))
	if (nsigma == 3) then ! Sigma is the principal stresses
		S = Sigma - I1*0.333333333333333333_8 ! Deviator stress vector
		J2 = 0.5_8*sum(S**2)
		J3 = S(1)*S(2)*S(3)
    
	elseif (nsigma == 4) then ! Plane situation, including axisymmetry
		S(1:3) = Sigma(1:3) - I1*0.333333333333333333_8 ! Deviator stress vector
		S(4) = Sigma(4)
		J2 = 0.5_8 * sum(S(1:3)**2) + S(4)**2
    
		J3 = (S(1)**3 + S(2)**3 + S(3)**3 +3.0_8*S(4)**2 *(S(1)+S(2)))
     &								*0.3333333333333333_8

	elseif (nsigma == 6) then ! General three-dimensional stress state
		S(1:3) = Sigma(1:3) - I1*0.33333333333333333_8 ! Deviator stress vector
		S(4:6) = Sigma(4:6)

		J2 = 0.5_8 * sum(S(1:3)**2) + sum(S(4:6)**2)

		J3 = ( sum(S(1:3)**3) + 6.0_8*S(4)*S(5)*S(6) +
     &		 3.0_8*(S(1)*(S(4)**2 + S(5)**2) +
     &			    S(2)*(S(4)**2 + S(6)**2) +
     &				S(3)*(S(5)**2 + S(6)**2)) )*0.3333333333333333_8
    
	end if

	if (J2 > 0.0_8) then
		sin3lode = -2.598076211353316_8 * J3 / J2**(1.5_8)
	else
		sin3lode = 0.0_8
	endif

	if (sin3lode <= -1.0_8) then
		lode = -0.5235987755982988_8 ! -pi/6
	elseif (sin3lode >= 1.0_8) then
		lode = +0.5235987755982988_8 ! +pi/6
	else
		lode = asin(sin3lode)/3.0_8
	endif

!----------------------------------------------------------------------------------
	end subroutine Invariants !--------------------------------------------------
!----------------------------------------------------------------------------------

	pure function Cross(A,B)
!--------------------------------------------------------------------------
! Calculates the cross product (vector product) between the two
! three dimensional vectors A and B
!
! INPUT
!	A, B	(real) Vectors with three components
!
! OUTPUT
!	Cross	(real) A vector perpendicular to both A and B
!
!------------------------------------------------------------------------
! Johan Clausen
! Esbjerg Department of Engineering
! Aalborg University
! July 2005
!------------------------------------------------------------------------

	real(8), intent(in) :: A(3), B(3)
!-------------------------

	real(8) Cross(3)

!---------------------------------------------

	Cross(1) = A(2)*B(3) - A(3)*B(2)
	Cross(2) = A(3)*B(1) - A(1)*B(3)
	Cross(3) = A(1)*B(2) - A(2)*B(1)

!--------------------------------------------------------------------------------
	end function Cross !-------------------------------------------------------
!--------------------------------------------------------------------------------

	pure subroutine MatInv_2x2_3x3(A,n,Ainv)
!-----------------------------------------------------
! Calculates the inverse of a square 1x1, 2x2 or 3x2 matrix
!
! INPUT
!  - Name- -type,size -		-- Description --
!	A		real,(n,n)		Square matrix
!	n		int,sc			size of A. n = 1, n = 2 or n = 3
!
! OUTPUT
!	Ainv	real,(n,n)		Inverse of A
!
!-----------------------------------------------------
! Johan Clausen
! Esbjerg Institute of Technology
! Aalborg University
! February 2006
!-----------------------------------------------------
	integer(4), intent(in) :: n
	real(8), dimension(n,n), intent(in) :: A

	real(8), dimension(n,n), intent(out) :: Ainv
	!-----------
	real(8) D(n,n), det

!---------------------------------
	if (n == 1) then
		Ainv = 1/A
	elseif (n == 2) then
		det = A(1,1)*A(2,2) - A(1,2)*A(2,1) ! Determinant
    
    		D(1,1) =  A(2,2) ;  D(1,2) = -A(1,2)
    		D(2,1) = -A(2,1) ;  D(2,2) =  A(1,1)

		Ainv = D/det

	elseif (n == 3) then
		D(1,1) = A(2,2)*A(3,3) - A(3,2)*A(2,3)
    		D(1,2) = A(1,3)*A(3,2) - A(3,3)*A(1,2)
    		D(1,3) = A(1,2)*A(2,3) - A(2,2)*A(1,3)
    		D(2,1) = A(2,3)*A(3,1) - A(3,3)*A(2,1)
    		D(2,2) = A(1,1)*A(3,3) - A(3,1)*A(1,3)
    		D(2,3) = A(1,3)*A(2,1) - A(2,3)*A(1,1)
    		D(3,1) = A(2,1)*A(3,2) - A(3,1)*A(2,2)
    		D(3,2) = A(1,2)*A(3,1) - A(3,2)*A(1,1)
    		D(3,3) = A(1,1)*A(2,2) - A(2,1)*A(1,2)

		det = A(1,1)*D(1,1) + A(1,2)*D(2,1) + A(1,3)*D(3,1) ! A's determinant

		Ainv = D/det
	end if
!----------------------------------------------------------------------
	end subroutine MatInv_2x2_3x3 !----------------------------------
!----------------------------------------------------------------------
