diff --git a/src/vrm.f90 b/src/vrm.f90
index 7414aefd4f5b64628bb5d07f191b46c0a3f1e940..98011dc27f86136f19e543e8885b9dcb5969a2a8 100644
--- a/src/vrm.f90
+++ b/src/vrm.f90
@@ -18,7 +18,7 @@
 
 subroutine vrm (params,J2d,J3d,compressibility,viscosity,pressure,             &
                plasticity_parameters,plasticity_type,is_plastic,flag_vrm_pb,   &
-               straintot)
+               straintot,call_cnt)
 
 use constants
 use definitions
@@ -42,6 +42,7 @@ double precision :: J2d,J3d,viscosity,pressure,compressibility,straintot
 character (len=8) plasticity_type
 double precision plasticity_parameters(9)
 logical ::  is_plastic,flag_vrm_pb
+integer call_cnt
 
 !------------------------------------------------------------------------------|
 !(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
@@ -108,8 +109,19 @@ select case (trim(plasticity_type))
     c     = plasticity_parameters(2)
     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)))
-    yield = max(c/compressibility,k+alpha*pressure)
-    fail  = 2.d0*viscosity*e2d-yield 
+    yield = max(c/compressibility,k-alpha*pressure)
+    fail  = 2.d0*viscosity*e2d-yield
+    !if (call_cnt.le.10) then
+    !  write(*,*) '******************************************'
+    !  write(*,*) 'phi: ',phi
+    !  write(*,*) 'c: ',c
+    !  write(*,*) 'alpha: ',alpha
+    !  write(*,*) 'k: ',k
+    !  write(*,*) 'pressure: ',pressure
+    !  write(*,*) 'yield: ',yield
+    !  write(*,*) 'fail: ',fail
+    !  write(*,*) '******************************************'
+    !endif
     if (fail.gt.0.d0) then
       viscosity=0.5d0*yield/e2d
       is_plastic=.true.
@@ -133,7 +145,7 @@ select case (trim(plasticity_type))
     c     = plasticity_parameters(2)
     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)))
-    yield = max(c/compressibility,k+alpha*pressure)
+    yield = max(c/compressibility,k-alpha*pressure)
     fail  = 2.d0*viscosity*e2d-yield 
     if (fail.gt.0.d0) then
       viscosity=0.5d0*yield/e2d
@@ -158,7 +170,7 @@ select case (trim(plasticity_type))
     c     = plasticity_parameters(2)
     alpha = (2.d0*sin(phi))  /(sqrt(3.d0)*3.d0)
     k     = (6.d0*c*cos(phi))/(sqrt(3.d0)*3.d0)
-    yield = max(c/compressibility,k+alpha*pressure)
+    yield = max(c/compressibility,k-alpha*pressure)
     fail  = 2.d0*viscosity*e2d-yield 
     if (fail.gt.0.d0) then
       viscosity=0.5d0*yield/e2d
@@ -183,7 +195,7 @@ select case (trim(plasticity_type))
     c     = plasticity_parameters(2)
     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))
-    yield = max(c/compressibility,k+alpha*pressure)
+    yield = max(c/compressibility,k-alpha*pressure)
     fail  = 2.d0*viscosity*e2d-yield 
     if (fail.gt.0.d0) then
       viscosity=0.5d0*yield/e2d
@@ -195,7 +207,7 @@ select case (trim(plasticity_type))
   case ('DPV') 
     alpha = plasticity_parameters(1)
     k     = plasticity_parameters(2)
-    yield = max(k/compressibility,k+alpha*pressure)
+    yield = max(k/compressibility,k-alpha*pressure)
     fail  = 2.d0*viscosity*e2d-yield 
     if (fail.gt.0.d0) then
       viscosity=0.5d0*yield/e2d
@@ -218,7 +230,7 @@ select case (trim(plasticity_type))
       endif
     endif
     c     = plasticity_parameters(2)
-    yield = max(c/compressibility,c+tan(phi)*pressure)
+    yield = max(c/compressibility,c-tan(phi)*pressure)
     fail  = 2.d0*viscosity*e2d-yield
     if (fail.gt.0.d0) then
       viscosity=0.5d0*yield/e2d
@@ -243,7 +255,7 @@ select case (trim(plasticity_type))
     c = plasticity_parameters(2)
     sin_phi=sin(phi)
     cos_phi=cos(phi)
-    yield = max(plasticity_parameters(2)/compressibility,pressure*sin_phi+c*cos_phi)
+    yield = max(plasticity_parameters(2)/compressibility,-pressure*sin_phi+c*cos_phi)
     fail  = 2.d0*viscosity*e2d-yield 
     if (fail.gt.0.d0) then
       viscosity=0.5d0*yield/e2d