!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              VRM    Feb. 2007                                                |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine vrm (params,J2d,J3d,compressibility,viscosity,pressure,             &
               friction_angle,plasticity_parameters,plasticity_type,           &
               plasticity_ss_type,is_plastic,flag_vrm_pb,straintot,yield)

use constants
use definitions
use invariants

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
!This routine computes the rescaled viscosity in the element in the case it is
!above the failure criterion

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

implicit none

type (parameters) params
double precision :: J2d,J3d,viscosity,pressure,compressibility,straintot
double precision :: friction_angle
character (len=8) plasticity_type
character (len=12) plasticity_ss_type
double precision plasticity_parameters(9)
logical ::  is_plastic,flag_vrm_pb

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

double precision :: yield,fail,theta,c,e2d,zeta,alpha,k
double precision :: sin_theta,cos_theta,sin_phi,cos_phi
double precision :: strain_soft_in,strain_soft_out,strain_soft_phi,phi,fact
double precision :: strain_soft_yield,eps,strain_soft_ref
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------

flag_vrm_pb = .false.
eps = 1.d-10
c = 0.d0
phi = 0.d0

if (params%init_e2d) then
  e2d=params%e2d0
else
  e2d=sqrt(J2d)
endif

select case (trim(plasticity_type))

  case ('vM')
    yield = plasticity_parameters(1)
    ! Set initial cohesion (yield) value for cohesion (yield) weakening
    strain_soft_in_c = plasticity_parameters(6)
    ! Call strain softening routine if initial cohesion (yield) is positive
    if (strain_soft_in_c.gt.0.d0) then
      call strain_soften (yield,phi,plasticity_ss_type_coh,                    &
                          plasticity_ss_type_phi,plasticity_parameters,        &
                          straintot,e2d,flag_vrm_pb)
    endif
    fail = 2.d0*viscosity*e2d-yield
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d
      is_plastic=.true.
    else
      is_plastic=.false.
    end if

  case ('DPI')
    ! Strain softening/hardening
    phi = plasticity_parameters(1)*pi/180.d0
    strain_soft_in_phi = plasticity_parameters(3)
    c = plasticity_parameters(2)
    strain_soft_in_c = plasticity_parameters(6)
    if (strain_soft_in_phi.gt.0.d0 .or. strain_soft_in_c.gt.0.d0) then
      call strain_soften (c,phi,plasticity_ss_type_coh,plasticity_ss_type_phi, &
                          plasticity_parameters,straintot,e2d,flag_vrm_pb)
    endif
    friction_angle = phi*180.d0/pi
    alpha = (2.d0*sin(phi))  /(sqrt(3.d0)*(3.d0-sin(phi)))
    k     = (6.d0*c*cos(phi))/(sqrt(3.d0)*(3.d0-sin(phi)))
    if (params%plastic_stress_min.lt.eps) then
      yield = max(c/compressibility,k-alpha*pressure)
    else
      yield = max(params%plastic_stress_min,k-alpha*pressure)
    endif
    fail = 2.d0*viscosity*e2d-yield
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d
      is_plastic=.true.
    else
      is_plastic=.false.
    end if

  case ('DPII')
    ! Strain softening/hardening
    phi = plasticity_parameters(1)*pi/180.d0
    strain_soft_in_phi = plasticity_parameters(3)
    c = plasticity_parameters(2)
    strain_soft_in_c = plasticity_parameters(6)
    if (strain_soft_in_phi.gt.0.d0 .or. strain_soft_in_c.gt.0.d0) then
      call strain_soften (c,phi,plasticity_ss_type_coh,plasticity_ss_type_phi, &
                          plasticity_parameters,straintot,e2d,flag_vrm_pb)
    endif
    friction_angle = phi*180.d0/pi
    alpha = (2.d0*sin(phi))  /(sqrt(3.d0)*(3.d0+sin(phi)))
    k     = (6.d0*c*cos(phi))/(sqrt(3.d0)*(3.d0+sin(phi)))
    if (params%plastic_stress_min.lt.eps) then
      yield = max(c/compressibility,k-alpha*pressure)
    else
      yield = max(params%plastic_stress_min,k-alpha*pressure)
    endif
    fail  = 2.d0*viscosity*e2d-yield 
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d
      is_plastic=.true.
    else
      is_plastic=.false.
    end if

  case ('DPIII')
    ! Strain softening/hardening
    phi = plasticity_parameters(1)*pi/180.d0
    strain_soft_in_phi = plasticity_parameters(3)
    c = plasticity_parameters(2)
    strain_soft_in_c = plasticity_parameters(6)
    if (strain_soft_in_phi.gt.0.d0 .or. strain_soft_in_c.gt.0.d0) then
      call strain_soften (c,phi,plasticity_ss_type_coh,plasticity_ss_type_phi, &
                          plasticity_parameters,straintot,e2d,flag_vrm_pb)
    endif
    friction_angle = phi*180.d0/pi
    alpha = (2.d0*sin(phi))  /(sqrt(3.d0)*3.d0)
    k     = (6.d0*c*cos(phi))/(sqrt(3.d0)*3.d0)
    if (params%plastic_stress_min.lt.eps) then
	  yield = max(c/compressibility,k-alpha*pressure)
	else
	  yield = max(params%plastic_stress_min,k-alpha*pressure)
	endif
    fail  = 2.d0*viscosity*e2d-yield 
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d
      is_plastic=.true.
    else
      is_plastic=.false.
    end if

  case ('DPIV')
    ! Strain softening/hardening
    phi = plasticity_parameters(1)*pi/180.d0
    strain_soft_in_phi = plasticity_parameters(3)
    c = plasticity_parameters(2)
    strain_soft_in_c = plasticity_parameters(6)
    if (strain_soft_in_phi.gt.0.d0 .or. strain_soft_in_c.gt.0.d0) then
      call strain_soften (c,phi,plasticity_ss_type_coh,plasticity_ss_type_phi, &
                          plasticity_parameters,straintot,e2d,flag_vrm_pb)
    endif
    friction_angle = phi*180.d0/pi
    alpha = tan(phi)/(sqrt(9.d0+12.d0*(tan(phi))**2.d0))
    k     = (3.d0*c)/(sqrt(9.d0+12.d0*(tan(phi))**2.d0))
    if (params%plastic_stress_min.lt.eps) then
      yield = max(c/compressibility,k-alpha*pressure)
    else
      yield = max(params%plastic_stress_min,k-alpha*pressure)
    endif
    fail  = 2.d0*viscosity*e2d-yield 
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d
      is_plastic=.true.
    else
      is_plastic=.false.
    end if

  case ('DPV') 
    alpha = plasticity_parameters(1)
    k     = plasticity_parameters(2)
    if (params%plastic_stress_min.lt.eps) then
      yield = max(k/compressibility,k-alpha*pressure)
    else
      yield = max(params%plastic_stress_min,k-alpha*pressure)
    endif
    fail  = 2.d0*viscosity*e2d-yield 
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d
      is_plastic=.true.
    else
      is_plastic=.false.
    end if
    
  case ('DPVI')
    ! Strain softening/hardening
    phi = plasticity_parameters(1)*pi/180.d0
    strain_soft_in_phi = plasticity_parameters(3)
    c = plasticity_parameters(2)
    strain_soft_in_c = plasticity_parameters(6)
    if (strain_soft_in_phi.gt.0.d0 .or. strain_soft_in_c.gt.0.d0) then
      call strain_soften (c,phi,plasticity_ss_type_coh,plasticity_ss_type_phi, &
                          plasticity_parameters,straintot,e2d,flag_vrm_pb)
    endif
    friction_angle = phi*180.d0/pi
    alpha = tan(phi)
    k     = c
    if (params%plastic_stress_min.lt.eps) then
      yield = max(c/compressibility,k-alpha*pressure)
    else
      yield = max(params%plastic_stress_min,k-alpha*pressure)
    endif
    fail  = 2.d0*viscosity*e2d-yield
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d
      is_plastic=.true.
    else
      is_plastic=.false.
    end if

  case ('DPVII')
    ! Strain softening/hardening
    phi = plasticity_parameters(1)*pi/180.d0
    strain_soft_in_phi = plasticity_parameters(3)
    c = plasticity_parameters(2)
    strain_soft_in_c = plasticity_parameters(6)
    if (strain_soft_in_phi.gt.0.d0 .or. strain_soft_in_c.gt.0.d0) then
      call strain_soften (c,phi,plasticity_ss_type_coh,plasticity_ss_type_phi, &
                          plasticity_parameters,straintot,e2d,flag_vrm_pb)
    endif
    friction_angle = phi*180.d0/pi
    alpha = sin(phi)
    k     = c*cos(phi)
    if (params%plastic_stress_min.lt.eps) then
      yield = max(c/compressibility,k-alpha*pressure)
    else
      yield = max(params%plastic_stress_min,k-alpha*pressure)
    endif
    fail  = 2.d0*viscosity*e2d-yield 
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d
      is_plastic=.true.
    else
      is_plastic=.false.
    end if

  case ('MC')
    if (params%init_e2d) then
      theta = 0.d0
    else
      theta = lode_angle(J2d,J3d)
    endif
    sin_theta = sin(theta)
    cos_theta = cos(theta)
    ! Strain softening/hardening
    phi = plasticity_parameters(1)*pi/180.d0
    strain_soft_in_phi = plasticity_parameters(3)
    c = plasticity_parameters(2)
    strain_soft_in_c = plasticity_parameters(6)
    if (strain_soft_in_phi.gt.0.d0 .or. strain_soft_in_c.gt.0.d0) then
      call strain_soften (c,phi,plasticity_ss_type_coh,plasticity_ss_type_phi, &
                          plasticity_parameters,straintot,e2d,flag_vrm_pb)
    endif
    friction_angle = phi*180.d0/pi
    sin_phi = sin(phi)
    cos_phi = cos(phi)
    zeta    = cos_theta-sin_phi*sin_theta/sqrt3
    if (params%plastic_stress_min.lt.eps) then
      yield = max(c/compressibility,(-pressure*sin_phi+c*cos_phi)/zeta)
    else
      yield = max(params%plastic_stress_min,(-pressure*sin_phi+c*cos_phi)/zeta)
    endif
    fail  = 2.d0*viscosity*e2d-yield
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d
      !if (viscosity.le.params%viscositymin) then
      !  write(*,*) 'VISCOSITY BELOW MINIMUM'
      !  write(*,*) 'viscosity: ',viscosity
      !  write(*,*) 'yield: ',yield
      !  write(*,*) 'fail: ',fail
      !  write(*,*) '-pressure: ',-pressure
      !  write(*,*) 'phi: ',phi
      !  write(*,*) 'sin_phi: ',sin_phi
      !  write(*,*) 'c: ',c
      !  write(*,*) 'cos_phi: ',cos_phi
      !  write(*,*) 'zeta: ',zeta
      !  write(*,*) 'e2d: ',e2d
      !endif
      is_plastic=.true.
    else
      is_plastic=.false.
    end if

  case ('Tresca')
    if (params%init_e2d) then
      theta = 0.d0
    else
      theta = lode_angle(J2d,J3d)
    endif
    yield = plasticity_parameters(1)
    ! Set initial cohesion (yield) value for cohesion (yield) weakening
    strain_soft_in_c = plasticity_parameters(6)
    ! Call strain softening routine if initial cohesion (yield) is positive
    if (strain_soft_in_c.gt.0.d0) then
      call strain_soften (yield,phi,plasticity_ss_type_coh,                    &
                          plasticity_ss_type_phi,plasticity_parameters,        &
                          straintot,e2d,flag_vrm_pb)
    endif
    cos_theta = cos(theta)
    fail  = 2.d0*viscosity*e2d*cos_theta-yield 
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d/cos_theta
      is_plastic=.true.
    else
      is_plastic=.false.
    end if

  case default
    print *,plasticity_type
    call stop_run('error in vrm$')
    flag_vrm_pb = .true.

end select

return
end subroutine vrm

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

subroutine strain_soften (c,phi,plasticity_ss_type_coh,plasticity_ss_type_phi, &
                          plasticity_parameters,straintot,e2d,flag_vrm_pb)

use constants

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! This routine calculates strain softened values of either the angle of internal
! friction or cohesion/yield strength as a function of total strain or strain
! rate

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

implicit none
character (len=12) plasticity_ss_type_coh,plasticity_ss_type_phi
double precision plasticity_parameters(9)
double precision :: c,phi,straintot,e2d
logical :: flag_vrm_pb

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

double precision :: strain_soft_phi_in,strain_soft_phi_out,strain_soft_phi
double precision :: strain_soft_c_in,strain_soft_c_out,strain_soft_phi

! Set strain softening onset/termination and min values
strain_soft_c_in    = plasticity_parameters(6)
strain_soft_c_out   = plasticity_parameters(7)
strain_soft_c       = plasticity_parameters(8)
strain_soft_phi_in  = plasticity_parameters(3)
strain_soft_phi_out = plasticity_parameters(4)
strain_soft_phi     = plasticity_parameters(5)*pi/180.d0

! Apply cohesion softening if onset cohesion value is positive
if (strain_soft_c_in.gt.0.d0) then
  ! Set type of strain softening
  ss_type_c: select case (trim(plasticity_ss_type_coh))
    case ('tot_strain')
      strain_soft_ref=straintot
    case ('strain_rate')
      strain_soft_ref=e2d
    case default
      write (*,*) 'plasticity_ss_type_coh: ',plasticity_ss_type_coh
      call stop_run('error in vrm$')
      flag_vrm_pb = .true.
  end select ss_type_c
  ! Soften cohesion if beyond onset value
  if (strain_soft_ref.gt.strain_soft_in_c) then
    fact=(strain_soft_ref-strain_soft_in_c)/(strain_soft_out_c-strain_soft_in_c)
    fact=min(fact,1.d0)
    c=c+(strain_soft_c-c)*fact
  endif
endif

! Apply internal angle softening if onset phi value is positive
if (strain_soft_phi_in.gt.0.d0) then
  ! Set type of strain softening
  ss_type_phi: select case (trim(plasticity_ss_type_phi))
    case ('tot_strain')
      strain_soft_ref=straintot
    case ('strain_rate')
      strain_soft_ref=e2d
    case default
      write (*,*) 'plasticity_ss_type_phi: ',plasticity_ss_type_phi
      call stop_run('error in vrm$')
      flag_vrm_pb = .true.
  end select ss_type_phi
  ! Soften phi if beyond onset value
  if (strain_soft_ref.gt.strain_soft_in_phi) then
    fact=(strain_soft_ref-strain_soft_in_phi)/(strain_soft_out_phi-strain_soft_in_phi)
    fact=min(fact,1.d0)
    phi=phi+(strain_soft_phi-phi)*fact
  endif
endif

end subroutine strain_soften

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