C
C ADJUSTING PARAMETERS TO FIT GIVEN MATERIAL TEST DATA
C 
C Using the free parameters to fit given data may sometimes lead to less
C than perfect fits of test data. In that case, it is necessary to
C adjust some of the fixed, or hard-to-change, parameters c_i. This is
C quite challenging. Improving the fit of one type of test by changing
C one parameter, c_j, may spoil the fit of another type of test. It may
C also necessitate changing several parameters c_j, which is a rather
C difficult nonlinear optimization problem. Attempts by several
C researchers to automatically adjust these parameters using
C optimization software has been unsuccessful in the past because of the
C complexity of the response. To succeed requires considerable
C experience. In that case, you can obtain assistance from
C Prof. F.C. Caner ( <ferhun.caner@upc.edu>, Tel. +34659816715, Skype
C ferhun.caner). He has ample experience in the fitting of test data
C with M7 using the finite element method.
C 
C The program is based on the following papers, which can be freely
C downloaded as 527.pdf, 528.pdf, 519.pdf and 547.pdf, respectively,
C from ZP Bazant's website:
C http://www.civil.northwestern.edu/people/bazant/PDFs/Papers/
C 
C [1] Caner, F.C., and Ba\v zant, Z.P. (2013). ``Microplane model M7 for
C     plain concrete: I. formulation." {\em ASCE J. of Engrg. Mechanics}
C     139 (12), Dec., 1714--1723.
C 
C [2] Caner, F.C., and Ba\v zant, Z.P. (2013). ``Microplane model M7 for
C     plain concrete: II. calibration and verification." {\em ASCE J. of
C     Engrg. Mechanics} 139 (12), Dec., 1724--1735.
C 
C [3] Caner, F.C., Ba\v zant, Z.P., and Wendner, R. (2013). ``Microplane
C     model M7f for fiber reinforced concrete." {\em Engrg. Fracture
C     Mechanics} 105, 41--57.
C 
C [4] Kirane, K., and Ba\v zant, Z.P. (2014). ``Microplane damage model
C     for fatigue of quasibrittle materials: Sub-critical crack growth,
C     lifetime and residual strength." {\em International Journal of
C     Fatigue} 70, 93--105.
C      
      subroutine vumat( 
     $     nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
     $     stepTime, totalTime, dt, cmname, coordMp, charLength,
     $     props, density, strainInc, relSpinInc,
     $     tempOld, stretchOld, defgradOld, fieldOld,
     $     stressOld, stateOld, enerInternOld, enerInelasOld,
     $     tempNew, stretchNew, defgradNew, fieldNew,
     $     stressNew, stateNew, enerInternNew, enerInelasNew )
      include 'vaba_param.inc'
      dimension props(nprops), density(nblock),
     $coordMp(nblock,*),
     $charLength(nblock), strainInc(nblock,ndir+nshr),
     $relSpinInc(nblock,nshr), tempOld(nblock),
     $stretchOld(nblock,ndir+nshr), defgradOld(nblock,ndir+nshr+nshr),
     $fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
     $stateOld(nblock,nstatev), enerInternOld(nblock),
     $enerInelasOld(nblock), tempNew(nblock),
     $stretchNew(nblock,ndir+nshr), defgradNew(nblock,ndir+nshr+nshr), 
     $fieldNew(nblock,nfieldv),
     $stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
     $enerInternNew(nblock), enerInelasNew(nblock)
      character*80 cmname

      IF(CMNAME(1:11).EQ."M7FCONCRETE") THEN
         CALL M7fMaterial( 
     $ nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
     $ stepTime, totalTime, dt, cmname, coordMp, charLength,
     $ props, density, strainInc, relSpinInc,
     $ tempOld, stretchOld, defgradOld, fieldOld,
     $ stressOld, stateOld, enerInternOld, enerInelasOld,
     $ tempNew, stretchNew, defgradNew, fieldNew,
     $ stressNew, stateNew, enerInternNew, enerInelasNew )
      END IF
      return
      END subroutine vumat

C +--------------------------------------------------------------------+
C |                 SUBROUTINE C_NORM_ELASTIC                          |
C +--------------------------------------------------------------------+
      subroutine c_norm_elastic(ef,def,sn0,eps_N0_pos,eps_N0_neg,sv0,
     $snf_no_split,E_N,zeta) 
      include 'vaba_param.inc'
      real, parameter :: PI=3.1415926535897932384626433832795d0
      integer, parameter ::  np=37, nvhm=5, nvhf=2, nvhi=np*nvhm+nvhf
      common /KCommon_Vars/young, poisson,k_1, C0,dtime,th_del,unit_conv  
      real :: young, poisson, k_1, C0, dtime, th_del, unit_conv, zeta
      common /KM7fBounds_M7fIO/k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9,
     $     sv0_p, c_1, c_2, c_3, c_4, c_5,
     $     c_6, c_7, c_8, c_9, c_10, c_11,c_12,c_13,c_14,c_15,c_16,
     $     c_17,c_18,c_19,c_20,c_21,C_R0,C_R2
      real :: k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9, sv0_p
      real :: c_1, c_2, c_3, c_4, c_5, c_6, c_7, c_8,
     $     c_9, c_10
      real :: c_11,c_12,c_13,c_14,c_15,c_16,c_17,c_18,
     $     c_19,c_20,c_21,C_R0,C_R2
      common /KM7fBounds_M7fIO_1/ c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1,
     $     c_6_2, c_6_M4, c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4,
     $     c_9_0, c_9_1, c_9_2, c_9_M4
      real :: c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1, c_6_2, c_6_M4
      real :: c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4
      real :: c_9_0, c_9_1, c_9_2, c_9_M4
      common /KM7fBounds_M7fIO_2/sig_fib_0, p_1, p_2, p_5, p_6,xp_1,xp_2,
     $     xp_3, V_f, cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3, cf_4,
     $     cf_5, cf_4_4, ap1, bp1, cp2, dp2, ep2, d_1, d_2, d_3, d_4,
     $     d_5, d_6
      real :: sig_fib_0, p_1, p_2, p_5, p_6
      real :: xp_1, xp_2, xp_3, V_f
      real :: cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3
      real :: cf_4, cf_5, cf_4_4
      real :: ap1, bp1, cp2, dp2, ep2
      real :: d_1, d_2, d_3, d_4, d_5, d_6
      save :: /KCommon_Vars/, /KM7fBounds_M7fIO/, /KM7fBounds_M7fIO_1/,
     $     /KM7fBounds_M7fIO_2/
      real, dimension(1:np) :: ef, def, sn0
      real, dimension(1:np) :: eps_N0_pos, eps_N0_neg
      real :: sv0
      real, dimension(1:np) :: snf_no_split 
      real, dimension(1:np) :: E_N

      real :: E_N_0
      real :: weight_factor, t0, t1, t2, t3, t4, tsum, fzeta
      integer :: isize, i 
      E_N_0= young / (1.0d0 - 2.0d0*poisson)

      isize=size(ef)

      do i=1,isize
         if ((sn0(i) > 0.d0)) then 
            weight_factor=4.0d1
            t0=1.d0
            t1=(weight_factor*zeta)**2.d0
            t2=(weight_factor*zeta)**4.d0
            t3=(weight_factor*zeta)**6.d0
            t4=(weight_factor*zeta)**8.d0
            tsum=t0+t1+t2
            fzeta=1.d0/tsum
            E_N(i) = E_N_0*fzeta*exp(-c_19*eps_N0_pos(i))
            if (sn0(i) > E_N_0*ef(i) .and. sn0(i)*def(i)<0.0d0) then 
               E_N(i)=E_N_0                
            end if                         
         else 
            E_N(i) = E_N_0*(exp(-c_20*abs(eps_N0_neg(i))/
     $           (1.d0+c_18*max(-sv0,0.d0)/E_N_0))+
     $           c_21*max(-sv0,0.d0)/E_N_0)
         end if
      end do

      snf_no_split = sn0*C0 + E_N*def

      return
      end subroutine c_norm_elastic

C +--------------------------------------------------------------------+
C |                         SUBROUTINE C_NORM                          |
C +--------------------------------------------------------------------+
      subroutine c_norm(ef,sv0,R_N,snb) 
      include 'vaba_param.inc'
      real, parameter :: PI=3.1415926535897932384626433832795d0
      integer, parameter ::  np=37, nvhm=5, nvhf=2, nvhi=np*nvhm+nvhf
      common /KCommon_Vars/young, poisson,k_1, C0,dtime,th_del,unit_conv  
      real :: young, poisson, k_1, C0, dtime, th_del, unit_conv
      common /KM7fBounds_M7fIO/k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9,
     $     sv0_p, c_1, c_2, c_3, c_4, c_5,
     $     c_6, c_7, c_8, c_9, c_10, c_11,c_12,c_13,c_14,c_15,c_16,
     $     c_17,c_18,c_19,c_20,c_21,C_R0,C_R2
      real :: k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9, sv0_p
      real :: c_1, c_2, c_3, c_4, c_5, c_6, c_7, c_8,
     $     c_9, c_10
      real :: c_11,c_12,c_13,c_14,c_15,c_16,c_17,c_18,
     $     c_19,c_20,c_21,C_R0,C_R2
      common /KM7fBounds_M7fIO_1/ c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1,
     $     c_6_2, c_6_M4, c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4,
     $     c_9_0, c_9_1, c_9_2, c_9_M4
      real :: c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1, c_6_2, c_6_M4
      real :: c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4
      real :: c_9_0, c_9_1, c_9_2, c_9_M4
      common /KM7fBounds_M7fIO_2/sig_fib_0, p_1, p_2, p_5,p_6,xp_1,xp_2,
     $     xp_3, V_f, cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3, cf_4,
     $     cf_5, cf_4_4, ap1, bp1, cp2, dp2, ep2, d_1, d_2, d_3, d_4,
     $     d_5, d_6
      real :: sig_fib_0, p_1, p_2, p_5, p_6
      real :: xp_1, xp_2, xp_3, V_f
      real :: cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3
      real :: cf_4, cf_5, cf_4_4
      real :: ap1, bp1, cp2, dp2, ep2
      real :: d_1, d_2, d_3, d_4, d_5, d_6
      save :: /KCommon_Vars/, /KM7fBounds_M7fIO/, /KM7fBounds_M7fIO_1/,
     $     /KM7fBounds_M7fIO_2/
      real, dimension(1:np) :: ef 
      real, dimension(1:np) :: snb 
      real :: sv0, R_N 
      real, dimension(1:np)  :: aux1
      real :: C_N0, RATE_FACTOR, fstar, beta_N
      real :: E_V, eb_N0, eb_NK, eb_N
      integer :: allocstat, isize

      isize = size(ef)

      E_V    = young/(1.0d0-2.0d0*poisson)
      eb_N0  = c_3*k_1 
      eb_NK  = c_4 

      c_1 = d_1*tanh(d_2*V_f - d_3)+d_4*exp(-max(-sv0-d_6,0.d0)/E_V*d_5) 

      if (sv0.lt.0.0d0) then 
         eb_N = eb_N0 - eb_NK / E_V * sv0  
      else 
         eb_N = eb_N0 
      end if       

      fstar  = k_1*young*c_1 
      beta_N = c_2*c_1*k_1 

      C_N0=C_R2 

      aux1 = max(ef-beta_N,0.0d0)/eb_N
      snb  = fstar*exp(-aux1)

      RATE_FACTOR = C_N0*R_N
      snb = snb * (1.d0 + RATE_FACTOR)

      return 
      end subroutine c_norm

! +--------------------------------------------------------------------+
! |                       SUBROUTINE C_NORM_FIB                        |
! +--------------------------------------------------------------------+
      subroutine c_norm_fib(ef,R_N,snb_fib) 
      include 'vaba_param.inc'
      real, parameter :: PI=3.1415926535897932384626433832795d0
      integer, parameter ::  np=37, nvhm=5, nvhf=2, nvhi=np*nvhm+nvhf
      common /KCommon_Vars/young, poisson,k_1, C0,dtime,th_del,unit_conv  
      real :: young, poisson, k_1, C0, dtime, th_del, unit_conv
      common /KM7fBounds_M7fIO/k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9,
     $     sv0_p, c_1, c_2, c_3, c_4, c_5,
     $     c_6, c_7, c_8, c_9, c_10, c_11,c_12,c_13,c_14,c_15,c_16,
     $     c_17,c_18,c_19,c_20,c_21,C_R0,C_R2
      real :: k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9, sv0_p
      real :: c_1, c_2, c_3, c_4, c_5, c_6, c_7, c_8,
     $     c_9, c_10
      real :: c_11,c_12,c_13,c_14,c_15,c_16,c_17,c_18,
     $     c_19,c_20,c_21,C_R0,C_R2
      common /KM7fBounds_M7fIO_1/ c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1,
     $     c_6_2, c_6_M4, c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4,
     $     c_9_0, c_9_1, c_9_2, c_9_M4
      real :: c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1, c_6_2, c_6_M4
      real :: c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4
      real :: c_9_0, c_9_1, c_9_2, c_9_M4
      common /KM7fBounds_M7fIO_2/sig_fib_0, p_1, p_2, p_5,p_6,xp_1,xp_2,
     $     xp_3, V_f, cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3, cf_4,
     $     cf_5, cf_4_4, ap1, bp1, cp2, dp2, ep2, d_1, d_2, d_3, d_4,
     $     d_5, d_6
      real :: sig_fib_0, p_1, p_2, p_5, p_6
      real :: xp_1, xp_2, xp_3, V_f
      real :: cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3
      real :: cf_4, cf_5, cf_4_4
      real :: ap1, bp1, cp2, dp2, ep2
      real :: d_1, d_2, d_3, d_4, d_5, d_6
      save :: /KCommon_Vars/, /KM7fBounds_M7fIO/, /KM7fBounds_M7fIO_1/,
     $     /KM7fBounds_M7fIO_2/
      real, dimension(1:np) :: ef 
      real, dimension(1:np) :: snb_fib 

      real :: R_N 

      real :: C_N0, RATE_FACTOR 

      real :: e_N_max, e_N_0, e_N_hat, eps_shift, pow 
      integer :: isize, i 
      
      isize = size(ef)

      eps_shift = xp_2
      pow = xp_3

      e_N_max = pow/p_2 + eps_shift

      e_N_0 = xp_1  

      e_N_hat = sig_fib_0 

      C_N0=C_R2 

      do i=1, isize
         if ( ef(i)/k_1 < e_N_max ) then
            snb_fib(i) = young*p_1*k_1*
     $           max(ef(i)/k_1 - eps_shift,0.0d0)**pow*
     $           exp(-p_2*max(ef(i)/k_1 - eps_shift,0.0d0)) 
         else if ((ef(i)/k_1 >= e_N_max) .and. (ef(i)/k_1 < e_N_0)) then
            snb_fib(i) = young*p_1*k_1*(pow/p_2)**pow*exp(-1.0d0*pow)
         else if (ef(i)/k_1 >= e_N_0) then 
            snb_fib(i) = young*p_1*k_1*
     $           (ef(i)/k_1 - e_N_0  + e_N_max - eps_shift)**pow*
     $           exp(-p_2*(ef(i)/k_1 - e_N_0 + e_N_max - eps_shift))  

         end if
      end do

      RATE_FACTOR = C_N0*R_N
      snb_fib = snb_fib * (1.d0 + RATE_FACTOR)
      return
      end subroutine c_norm_fib

! +--------------------------------------------------------------------+
! |                      SUBROUTINE C_DEV                              |
! +--------------------------------------------------------------------+
      subroutine c_dev(ded, ed0, dev, ev0, sv0, sd0, sdf, C_d, R_D,
     $     sdneg, sdpos) 
      include 'vaba_param.inc'
      real, parameter :: PI=3.1415926535897932384626433832795d0
      integer, parameter ::  np=37, nvhm=5, nvhf=2, nvhi=np*nvhm+nvhf
      common /KCommon_Vars/young, poisson,k_1, C0,dtime,th_del,unit_conv  
      real :: young, poisson, k_1, C0, dtime, th_del, unit_conv
      common /KM7fBounds_M7fIO/k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9,
     $     sv0_p, c_1, c_2, c_3, c_4, c_5,
     $     c_6, c_7, c_8, c_9, c_10, c_11,c_12,c_13,c_14,c_15,c_16,
     $     c_17,c_18,c_19,c_20,c_21,C_R0,C_R2
      real :: k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9, sv0_p
      real :: c_1, c_2, c_3, c_4, c_5, c_6, c_7, c_8,
     $     c_9, c_10
      real :: c_11,c_12,c_13,c_14,c_15,c_16,c_17,c_18,
     $     c_19,c_20,c_21,C_R0,C_R2
      common /KM7fBounds_M7fIO_1/ c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1,
     $     c_6_2, c_6_M4, c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4,
     $     c_9_0, c_9_1, c_9_2, c_9_M4
      real :: c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1, c_6_2, c_6_M4
      real :: c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4
      real :: c_9_0, c_9_1, c_9_2, c_9_M4
      common /KM7fBounds_M7fIO_2/sig_fib_0, p_1, p_2,p_5,p_6,xp_1,xp_2,
     $     xp_3, V_f, cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3, cf_4,
     $     cf_5, cf_4_4, ap1, bp1, cp2, dp2, ep2, d_1, d_2, d_3, d_4,
     $     d_5, d_6
      real :: sig_fib_0, p_1, p_2, p_5, p_6
      real :: xp_1, xp_2, xp_3, V_f
      real :: cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3
      real :: cf_4, cf_5, cf_4_4
      real :: ap1, bp1, cp2, dp2, ep2
      real :: d_1, d_2, d_3, d_4, d_5, d_6
      save :: /KCommon_Vars/, /KM7fBounds_M7fIO/, /KM7fBounds_M7fIO_1/,
     $     /KM7fBounds_M7fIO_2/
      real, dimension(1:np) :: ded, ed0, sd0
      real, dimension(1:np) :: sdf, C_d

      real, dimension(1:np) :: sdneg, sdpos

      real :: dev, ev0, sv0, R_D
      real :: usi_C, usi_T, par5, par6, E_s, Cd0, Ev
      real, dimension(1:np) :: edf, sde, sdlower, sdupper
      real, dimension(1:np) :: sdlower_fib, sdupper_fib
      real :: RATE_FACTOR_C, RATE_FACTOR_T, C_DC0, C_DT0, evf, E_v

      real :: f_c0, E_0, f_cp, c_40, beta_15, beta_16, beta_17, beta_18,
     $     beta_19 
      logical :: l0,l1,l2,l3,l4
      integer :: j, isize, allocstat
      
      isize = size(ded)

      evf = ev0 + dev

      Ev=young/(1.0d0 - 2.0d0*poisson)

      f_c0=15.08d0
      E_0=20000.d0
      c_40=1.0d+0

      f_cp=90.3
      beta_15=c_5_1*exp(-c_40*(f_cp/young-f_c0/E_0))
      beta_16=c_8_1*exp(-c_40*(f_cp/young-f_c0/E_0))
      beta_17=c_7_1*exp(-c_40*(f_cp/young-f_c0/E_0))
      beta_18=c_6_0*exp(-c_40*(f_cp/young-f_c0/E_0))
      beta_19=c_9_0*exp(-c_40*(f_cp/young-f_c0/E_0))

      c_5 = (beta_15*tanh(c_5_0*max(-evf,0.0d0)/k_1))*
     $     (1.0d0 + cf_1*tanh(cf_1_1*V_f))+cf_1*tanh(cf_1_1*V_f)+c_5_M4

      c_8 = (beta_16*tanh(c_8_0*max(-evf,0.0d0)/k_1))*
     $     (1.0d0 + cf_3*tanh(cf_3_3*V_f))+cf_3*tanh(cf_3_3*V_f)+c_8_M4

      c_7 = beta_17*tanh(c_7_0*max(-evf,0.0d0)/k_1)+
     $     cf_2*tanh(cf_2_2*V_f) + c_7_M4 

      if (beta_18*max(-evf/k_1-c_6_1,0.0d0)>=log(c_6_2)) then
         c_6 = c_6_M4*c_6_2
      else
         c_6 = c_6_M4*exp(beta_18*max(-evf/k_1-c_6_1,0.0d0))
      end if
      if (beta_19*max(-evf/k_1-c_9_1,0.0d0)>=log(c_9_2)) then
         c_9 = c_9_M4*c_9_2
      else
         c_9 = c_9_M4*exp(beta_19*max(-evf/k_1-c_9_1,0.0d0))
      end if

      par5 = k_1*young*c_8
      par6 = young*c_5*k_1
      Cd0=young/(1.d0+poisson)*(1.d0-4.d0*poisson)/(1.d0-2.d0*poisson)
      E_v  = young/(1.d0-2*poisson)

      C_DC0  = C_R2

      C_DT0  = C_R2

      C_d=Cd0

      edf = ed0+ded 

      call dev_comp(edf, sdlower) 

      call dev_tens(edf, sdupper)
      call dev_tens_fib(edf, sdupper_fib)

      RATE_FACTOR_C = C_DC0*R_D
      RATE_FACTOR_T = C_DT0*R_D

      sdlower = sdlower * (1.d0 + RATE_FACTOR_C)      

      sdupper = sdupper * (1.d0 + RATE_FACTOR_T)

      sdneg=sdlower
      sdpos=sdupper

      return
      end subroutine c_dev

! +--------------------------------------------------------------------+
! |                        SUBROUTINE DEV_COMP                         |
! +--------------------------------------------------------------------+
      subroutine dev_comp(edf, sdlower) 
      include 'vaba_param.inc'
      real, parameter :: PI=3.1415926535897932384626433832795d0
      integer, parameter ::  np=37, nvhm=5, nvhf=2, nvhi=np*nvhm+nvhf
      common /KCommon_Vars/young, poisson,k_1, C0,dtime,th_del,unit_conv  
      real :: young, poisson, k_1, C0, dtime, th_del, unit_conv
      common /KM7fBounds_M7fIO/k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9,
     $     sv0_p, c_1, c_2, c_3, c_4, c_5,
     $     c_6, c_7, c_8, c_9, c_10, c_11,c_12,c_13,c_14,c_15,c_16,
     $     c_17,c_18,c_19,c_20,c_21,C_R0,C_R2
      real :: k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9, sv0_p
      real :: c_1, c_2, c_3, c_4, c_5, c_6, c_7, c_8,
     $     c_9, c_10
      real :: c_11,c_12,c_13,c_14,c_15,c_16,c_17,c_18,
     $     c_19,c_20,c_21,C_R0,C_R2
      common /KM7fBounds_M7fIO_1/ c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1,
     $     c_6_2, c_6_M4, c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4,
     $     c_9_0, c_9_1, c_9_2, c_9_M4
      real :: c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1, c_6_2, c_6_M4
      real :: c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4
      real :: c_9_0, c_9_1, c_9_2, c_9_M4
      common /KM7fBounds_M7fIO_2/sig_fib_0, p_1, p_2, p_5,p_6,xp_1,xp_2,
     $     xp_3, V_f, cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3, cf_4,
     $     cf_5, cf_4_4, ap1, bp1, cp2, dp2, ep2, d_1, d_2, d_3, d_4,
     $     d_5, d_6
      real :: sig_fib_0, p_1, p_2, p_5, p_6
      real :: xp_1, xp_2, xp_3, V_f
      real :: cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3
      real :: cf_4, cf_5, cf_4_4
      real :: ap1, bp1, cp2, dp2, ep2
      real :: d_1, d_2, d_3, d_4, d_5, d_6
      save :: /KCommon_Vars/, /KM7fBounds_M7fIO/, /KM7fBounds_M7fIO_1/,
     $     /KM7fBounds_M7fIO_2/
      real, dimension(1:np) :: edf

      real, dimension(1:np) :: sdlower
      real :: par1, par2, par5

      par1 = c_9*c_8*k_1
      par2 = c_7*k_1 
      par5 = k_1*young*c_8

      sdlower=-par5/(1.d0+(max(-edf-par1,0.0d0)/par2)**2.0d0) 

      return 
      end subroutine dev_comp 
C
C
C
C +--------------------------------------------------------------------+
C |                         SUBROUTINE DEV_TENS                        |
C +--------------------------------------------------------------------+

      subroutine dev_tens(edf, sdupper) 
      include 'vaba_param.inc'
      real, parameter :: PI=3.1415926535897932384626433832795d0
      integer, parameter ::  np=37, nvhm=5, nvhf=2, nvhi=np*nvhm+nvhf
      common /KCommon_Vars/young, poisson,k_1, C0,dtime,th_del,unit_conv  
      real :: young, poisson, k_1, C0, dtime, th_del, unit_conv
      common /KM7fBounds_M7fIO/k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9,
     $     sv0_p, c_1, c_2, c_3, c_4, c_5,
     $     c_6, c_7, c_8, c_9, c_10, c_11,c_12,c_13,c_14,c_15,c_16,
     $     c_17,c_18,c_19,c_20,c_21,C_R0,C_R2
      real :: k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9, sv0_p
      real :: c_1, c_2, c_3, c_4, c_5, c_6, c_7, c_8,
     $     c_9, c_10
      real :: c_11,c_12,c_13,c_14,c_15,c_16,c_17,c_18,
     $     c_19,c_20,c_21,C_R0,C_R2
      common /KM7fBounds_M7fIO_1/ c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1,
     $     c_6_2, c_6_M4, c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4,
     $     c_9_0, c_9_1, c_9_2, c_9_M4
      real :: c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1, c_6_2, c_6_M4
      real :: c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4
      real :: c_9_0, c_9_1, c_9_2, c_9_M4
      common /KM7fBounds_M7fIO_2/sig_fib_0, p_1, p_2, p_5,p_6,xp_1,xp_2,
     $     xp_3, V_f, cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3, cf_4,
     $     cf_5, cf_4_4, ap1, bp1, cp2, dp2, ep2, d_1, d_2, d_3, d_4,
     $     d_5, d_6
      real :: sig_fib_0, p_1, p_2, p_5, p_6
      real :: xp_1, xp_2, xp_3, V_f
      real :: cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3
      real :: cf_4, cf_5, cf_4_4
      real :: ap1, bp1, cp2, dp2, ep2
      real :: d_1, d_2, d_3, d_4, d_5, d_6
      save :: /KCommon_Vars/, /KM7fBounds_M7fIO/, /KM7fBounds_M7fIO_1/,
     $     /KM7fBounds_M7fIO_2/
      real, dimension(1:np) :: edf

      real, dimension(1:np) :: sdupper
      real :: par2, par3, par4, par6

      par2 = c_7*k_1 
      par3 = c_20
      par4 = c_6*c_5*k_1
      par6 = young*c_5*k_1

      sdupper=par6/(1.0d0+(max(edf-par4,0.0d0)/par2/par3)**2.0d0)

      return 
      end subroutine dev_tens 

! +--------------------------------------------------------------------+
! |                         SUBROUTINE DEV_TENS_FIB                    |
! +--------------------------------------------------------------------+
      subroutine dev_tens_fib(edf, sdupper_fib) 
      include 'vaba_param.inc'
      real, parameter :: PI=3.1415926535897932384626433832795d0
      integer, parameter ::  np=37, nvhm=5, nvhf=2, nvhi=np*nvhm+nvhf
      common /KCommon_Vars/young, poisson,k_1, C0,dtime,th_del,unit_conv  
      real :: young, poisson, k_1, C0, dtime, th_del, unit_conv
      common /KM7fBounds_M7fIO/k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9,
     $     sv0_p, c_1, c_2, c_3, c_4, c_5,
     $     c_6, c_7, c_8, c_9, c_10, c_11,c_12,c_13,c_14,c_15,c_16,
     $     c_17,c_18,c_19,c_20,c_21,C_R0,C_R2
      real :: k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9, sv0_p
      real :: c_1, c_2, c_3, c_4, c_5, c_6, c_7, c_8,
     $     c_9, c_10
      real :: c_11,c_12,c_13,c_14,c_15,c_16,c_17,c_18,
     $     c_19,c_20,c_21,C_R0,C_R2
      common /KM7fBounds_M7fIO_1/ c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1,
     $     c_6_2, c_6_M4, c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4,
     $     c_9_0, c_9_1, c_9_2, c_9_M4
      real :: c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1, c_6_2, c_6_M4
      real :: c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4
      real :: c_9_0, c_9_1, c_9_2, c_9_M4
      common /KM7fBounds_M7fIO_2/sig_fib_0, p_1, p_2, p_5,p_6,xp_1,xp_2,
     $     xp_3, V_f, cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3, cf_4,
     $     cf_5, cf_4_4, ap1, bp1, cp2, dp2, ep2, d_1, d_2, d_3, d_4,
     $     d_5, d_6
      real :: sig_fib_0, p_1, p_2, p_5, p_6
      real :: xp_1, xp_2, xp_3, V_f
      real :: cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3
      real :: cf_4, cf_5, cf_4_4
      real :: ap1, bp1, cp2, dp2, ep2
      real :: d_1, d_2, d_3, d_4, d_5, d_6
      save :: /KCommon_Vars/, /KM7fBounds_M7fIO/, /KM7fBounds_M7fIO_1/,
     $     /KM7fBounds_M7fIO_2/
      real, dimension(1:np) :: edf
      real, dimension(1:np) :: sdupper_fib

      sdupper_fib=young*p_5*k_1*max(edf/k_1,0.0d0)*
     $     exp(-p_6*max(edf/k_1,0.0d0))

      return 
      end subroutine dev_tens_fib 

      subroutine c_shear_tens(sig_o, s_0, fsp_0)
      include 'vaba_param.inc'
      real, parameter :: PI=3.1415926535897932384626433832795d0
      integer, parameter ::  np=37, nvhm=5, nvhf=2, nvhi=np*nvhm+nvhf
      common /KCommon_Vars/young, poisson,k_1, C0,dtime,th_del,unit_conv  
      real :: young, poisson, k_1, C0, dtime, th_del, unit_conv
      common /KM7fBounds_M7fIO/k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9,
     $     sv0_p, c_1, c_2, c_3, c_4, c_5,
     $     c_6, c_7, c_8, c_9, c_10, c_11,c_12,c_13,c_14,c_15,c_16,
     $     c_17,c_18,c_19,c_20,c_21,C_R0,C_R2
      real :: k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9, sv0_p
      real :: c_1, c_2, c_3, c_4, c_5, c_6, c_7, c_8,
     $     c_9, c_10
      real :: c_11,c_12,c_13,c_14,c_15,c_16,c_17,c_18,
     $     c_19,c_20,c_21,C_R0,C_R2
      common /KM7fBounds_M7fIO_1/ c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1,
     $     c_6_2, c_6_M4, c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4,
     $     c_9_0, c_9_1, c_9_2, c_9_M4
      real :: c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1, c_6_2, c_6_M4
      real :: c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4
      real :: c_9_0, c_9_1, c_9_2, c_9_M4
      common /KM7fBounds_M7fIO_2/sig_fib_0, p_1, p_2, p_5,p_6,xp_1,xp_2,
     $     xp_3, V_f, cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3, cf_4,
     $     cf_5, cf_4_4, ap1, bp1, cp2, dp2, ep2, d_1, d_2, d_3, d_4,
     $     d_5, d_6
      real :: sig_fib_0, p_1, p_2, p_5, p_6
      real :: xp_1, xp_2, xp_3, V_f
      real :: cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3
      real :: cf_4, cf_5, cf_4_4
      real :: ap1, bp1, cp2, dp2, ep2
      real :: d_1, d_2, d_3, d_4, d_5, d_6
      save :: /KCommon_Vars/, /KM7fBounds_M7fIO/, /KM7fBounds_M7fIO_1/,
     $     /KM7fBounds_M7fIO_2/
      real :: sig_o, s_0
      real :: fsp_0
      real :: c_6_p

      c_6_p = c_10
      fsp_0 = c_6_p * max(sig_o, 0.0d0) / (1.d0 + c_6_p / s_0 *
     $     max(sig_o, 0.0d0)) 
      return
      end subroutine c_shear_tens

! +--------------------------------------------------------------------+
! |                          SUBROUTINE C_SHEAR2                       |
! +--------------------------------------------------------------------+
      subroutine c_shear2(eps_L, eps_M, C_d, snf, R_S, del, dem, sl0,
     $     sm0, slf, smf, eps_V) 
      include 'vaba_param.inc'
      real, parameter :: PI=3.1415926535897932384626433832795d0
      integer, parameter ::  np=37, nvhm=5, nvhf=2, nvhi=np*nvhm+nvhf
      common /KCommon_Vars/young, poisson,k_1, C0,dtime,th_del,unit_conv  
      real :: young, poisson, k_1, C0, dtime, th_del, unit_conv
      common /KM7fBounds_M7fIO/k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9,
     $     sv0_p, c_1, c_2, c_3, c_4, c_5,
     $     c_6, c_7, c_8, c_9, c_10, c_11,c_12,c_13,c_14,c_15,c_16,
     $     c_17,c_18,c_19,c_20,c_21,C_R0,C_R2
      real :: k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9, sv0_p
      real :: c_1, c_2, c_3, c_4, c_5, c_6, c_7, c_8,
     $     c_9, c_10
      real :: c_11,c_12,c_13,c_14,c_15,c_16,c_17,c_18,
     $     c_19,c_20,c_21,C_R0,C_R2
      common /KM7fBounds_M7fIO_1/ c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1,
     $     c_6_2, c_6_M4, c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4,
     $     c_9_0, c_9_1, c_9_2, c_9_M4
      real :: c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1, c_6_2, c_6_M4
      real :: c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4
      real :: c_9_0, c_9_1, c_9_2, c_9_M4
      common /KM7fBounds_M7fIO_2/sig_fib_0, p_1, p_2, p_5,p_6,xp_1,xp_2,
     $     xp_3, V_f, cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3, cf_4,
     $     cf_5, cf_4_4, ap1, bp1, cp2, dp2, ep2, d_1, d_2, d_3, d_4,
     $     d_5, d_6
      real :: sig_fib_0, p_1, p_2, p_5, p_6
      real :: xp_1, xp_2, xp_3, V_f
      real :: cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3
      real :: cf_4, cf_5, cf_4_4
      real :: ap1, bp1, cp2, dp2, ep2
      real :: d_1, d_2, d_3, d_4, d_5, d_6
      save :: /KCommon_Vars/, /KM7fBounds_M7fIO/, /KM7fBounds_M7fIO_1/,
     $     /KM7fBounds_M7fIO_2/
      real, dimension(1:np)  :: eps_L, eps_M, C_d, snf, del,
     $     dem, sl0, sm0 
      real  :: R_S, eps_V
      real, dimension(1:np) :: slf, smf  
      real :: E_T, c_12_p, c_13_p, c_6_p, s_0
      real, dimension(1:np) :: C_t, fsp
      real, dimension(1:np) :: stf, fsp_0, epsf, depsf
      real :: RATE_FACTOR_1, C_S0, sig_o
      integer :: isize, allocstat, j
      
      isize = size(snf)

      E_T=young/(1.d0+poisson)*(1.d0-4.d0*poisson)/(1.d0-2.d0*poisson) 
      c_6_p  = c_10
      c_12_p = c_11 
      c_13_p = c_12*3.371d-4

      s_0    = k_1*k_2*E_T

      C_S0 = C_R2

      RATE_FACTOR_1 = C_S0*R_S  

      sig_o = max(E_T * k_1 *
     $     (c_12_p - c_13_p*max(eps_V, 0.0d0)/k_1),0.0d0)

      fsp = c_6_p * max(-snf + sig_o, 0.0d0) / (1.d0 + c_6_p / s_0 *
     $     max(-snf + sig_o, 0.0d0)) * (1.d0 + RATE_FACTOR_1)

      epsf = sqrt(eps_L*eps_L + eps_M*eps_M)
      depsf = sqrt(del*del + dem*dem)
      do j=1, isize
         call c_shear_tens(sig_o, s_0, fsp_0(j))
         fsp_0(j)= fsp_0(j) * (1.d0 + RATE_FACTOR_1)
      end do

      do j=1,isize
        if(sl0(j)*del(j) < 0.0d0 .or. sm0(j)*dem(j) < 0.0d0) then
          C_t(j) = C_d(j)
        else                 
          C_t(j) = E_T
        end if  
      end do

      slf = sl0*C0 + C_t*del
      smf = sm0*C0 + C_t*dem
      stf = sqrt((slf*slf + smf*smf))
      do j = 1, isize
            if (stf(j) /= 0.0d0) then 
               slf(j) = slf(j) / stf(j)
               smf(j) = smf(j) / stf(j)
            end if
      end do

      do j=1, isize
         if (snf(j) < 0.0 ) then 
            if (stf(j) > fsp(j)) then
               stf(j) = fsp(j)
            end if
         else
            if (stf(j) > fsp_0(j)) then
               stf(j) = fsp_0(j)
            end if
         end if
         slf(j) = stf(j) * slf(j)
         smf(j) = stf(j) * smf(j)
      end do

      return 
      end subroutine c_shear2 

C +--------------------------------------------------------------------+
C |                        SUBROUTINE C_VOL                            |
C +--------------------------------------------------------------------+
      subroutine c_vol(dev, ev0, sv0, deps_N, eps_N, R_N, svneg) 
      include 'vaba_param.inc'
      external fb
      real, parameter :: PI=3.1415926535897932384626433832795d0
      integer, parameter ::  np=37, nvhm=5, nvhf=2, nvhi=np*nvhm+nvhf
      common /KCommon_Vars/young, poisson,k_1, C0,dtime,th_del,unit_conv  
      real :: young, poisson, k_1, C0, dtime, th_del, unit_conv
      common /KM7fBounds_M7fIO/k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9,
     $     sv0_p, c_1, c_2, c_3, c_4, c_5,
     $     c_6, c_7, c_8, c_9, c_10, c_11,c_12,c_13,c_14,c_15,c_16,
     $     c_17,c_18,c_19,c_20,c_21,C_R0,C_R2
      real :: k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9, sv0_p
      real :: c_1, c_2, c_3, c_4, c_5, c_6, c_7, c_8,
     $     c_9, c_10
      real :: c_11,c_12,c_13,c_14,c_15,c_16,c_17,c_18,
     $     c_19,c_20,c_21,C_R0,C_R2
      common /KM7fBounds_M7fIO_1/ c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1,
     $     c_6_2, c_6_M4, c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4,
     $     c_9_0, c_9_1, c_9_2, c_9_M4
      real :: c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1, c_6_2, c_6_M4
      real :: c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4
      real :: c_9_0, c_9_1, c_9_2, c_9_M4
      common /KM7fBounds_M7fIO_2/sig_fib_0, p_1, p_2, p_5,p_6,xp_1,xp_2,
     $     xp_3, V_f, cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3, cf_4,
     $     cf_5, cf_4_4, ap1, bp1, cp2, dp2, ep2, d_1, d_2, d_3, d_4,
     $     d_5, d_6
      real :: sig_fib_0, p_1, p_2, p_5, p_6
      real :: xp_1, xp_2, xp_3, V_f
      real :: cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3
      real :: cf_4, cf_5, cf_4_4
      real :: ap1, bp1, cp2, dp2, ep2
      real :: d_1, d_2, d_3, d_4, d_5, d_6
      save :: /KCommon_Vars/, /KM7fBounds_M7fIO/, /KM7fBounds_M7fIO_1/,
     $     /KM7fBounds_M7fIO_2/
      real :: dev, ev0, sv0
      real, dimension(1:np) :: eps_N, deps_N
      real :: svneg 
      real :: R_N 
      real :: Cv0, evf, svb_neg 
      real :: xk0, e0_V

      real :: prStrainDiff, xk4 
      
      xk0  = k_3*k_1*young 
      e0_V = k_4*k_1 
      Cv0  = young / (1.0d0 - 2.0d0 * poisson) 

      prStrainDiff = maxval(eps_N+deps_N,1) - minval(eps_N+deps_N,1)

      xk4 = (k_6*(prStrainDiff/k_1)**k_7)
     $     /(1.0d0 + min(max(-sv0,0.d0),sv0_p)/Cv0) + k_4

      e0_V = xk4*k_1

      evf = ev0+dev 

      svb_neg = fb(0,evf,xk0,e0_V) 
      svneg=svb_neg*(1.0d0 + C_R2*R_N)

      end subroutine c_vol
 
C +--------------------------------------------------------------------+
C |                            FUNCTION FB                             |    
C +--------------------------------------------------------------------+
      function fb(i,ef,xk0,e0_V) 
      include 'vaba_param.inc'
      real, parameter :: PI=3.1415926535897932384626433832795d0
      integer, parameter ::  np=37, nvhm=5, nvhf=2, nvhi=np*nvhm+nvhf
      common /KCommon_Vars/young, poisson,k_1, C0,dtime,th_del,unit_conv  
      real :: young, poisson, k_1, C0, dtime, th_del, unit_conv
      common /KM7fBounds_M7fIO/k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9,
     $     sv0_p, c_1, c_2, c_3, c_4, c_5,
     $     c_6, c_7, c_8, c_9, c_10, c_11,c_12,c_13,c_14,c_15,c_16,
     $     c_17,c_18,c_19,c_20,c_21,C_R0,C_R2
      real :: k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9, sv0_p
      real :: c_1, c_2, c_3, c_4, c_5, c_6, c_7, c_8,
     $     c_9, c_10
      real :: c_11,c_12,c_13,c_14,c_15,c_16,c_17,c_18,
     $     c_19,c_20,c_21,C_R0,C_R2
      common /KM7fBounds_M7fIO_1/ c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1,
     $     c_6_2, c_6_M4, c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4,
     $     c_9_0, c_9_1, c_9_2, c_9_M4
      real :: c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1, c_6_2, c_6_M4
      real :: c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4
      real :: c_9_0, c_9_1, c_9_2, c_9_M4
      common /KM7fBounds_M7fIO_2/sig_fib_0, p_1, p_2, p_5,p_6,xp_1,xp_2,
     $     xp_3, V_f, cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3, cf_4,
     $     cf_5, cf_4_4, ap1, bp1, cp2, dp2, ep2, d_1, d_2, d_3, d_4,
     $     d_5, d_6
      real :: sig_fib_0, p_1, p_2, p_5, p_6
      real :: xp_1, xp_2, xp_3, V_f
      real :: cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3
      real :: cf_4, cf_5, cf_4_4
      real :: ap1, bp1, cp2, dp2, ep2
      real :: d_1, d_2, d_3, d_4, d_5, d_6
      save :: /KCommon_Vars/, /KM7fBounds_M7fIO/, /KM7fBounds_M7fIO_1/,
     $     /KM7fBounds_M7fIO_2/

      integer ::  i
      real :: ef,xk0,e0_V
      real :: fb, aux, aux2
      
      aux=-xk0*exp(-ef/e0_V) 
      if (i.eq.0) fb = aux 

      if (i.eq.1) fb = -aux/e0_V 
      return
      end function fb

C **********************************************************************
C *** SUBROUTINE SETSYSTEM *********************************************
C **********************************************************************
 
      subroutine setsystem() 
      include 'vaba_param.inc'
      real, parameter :: PI=3.1415926535897932384626433832795d0
      integer, parameter ::  np=37, nvhm=5, nvhf=2, nvhi=np*nvhm+nvhf
      common /KCommon_Vars/young, poisson,k_1, C0,dtime,th_del,unit_conv  
      real :: young, poisson, k_1, C0, dtime, th_del, unit_conv, zeta
      real :: zeta0
      common /KM7f_Microplane_1/ qn, ql, qm, w 
      common /KM7f_Microplane_2/ deps_N, deps_L, deps_M, eps_N,  eps_L,
     $     eps_M,  sig_N
      common /KM7f_Microplane_3/ sd0, sn0, ded, ed0, edf, enf, snf, snb,
     $     sdf
      common /KM7f_Microplane_4/ sl0, sm0, slf, smf, del, dem, vh_ini,
     $     vh_fin, C_d, E_N
      common /KM7f_Microplane_5/ sdneg, sdpos, eps_N0_neg, eps_N0_pos,
     $     snb_fib, Cv, zeta, zeta0
      real, dimension(1:6,1:np) :: qn, ql, qm
      real, dimension(1:np)     :: w
      real, dimension(1:np)     :: deps_N,
     $     deps_L, deps_M, eps_N,  eps_L,  eps_M,  sig_N
      real, dimension(1:np)     :: sd0, sn0,
     $     ded, ed0, edf
      real, dimension(1:np)     :: enf, snf,
     $     snb, sdf, sdneg, sdpos, snb_fib
      real, dimension(1:np)     :: sl0, sm0, slf,
     $     smf, del, dem, eps_N0_neg, eps_N0_pos
      real, dimension(1:nvhi+2)    :: vh_ini, vh_fin 
      real, dimension(1:np)     :: C_d, E_N
      real :: Cv
      save    :: /KCommon_Vars/
      save    :: /KM7f_Microplane_1/,/KM7f_Microplane_2/,
     $     /KM7f_Microplane_3/,/KM7f_Microplane_4/,/KM7f_Microplane_5/
      integer :: jp, ij(1:2,1:6), i, j, k, allocstat
      integer :: rs_size
      real, dimension(1:4,1:np) :: te
      real, dimension(1:3) :: xn, xm, xl, rand_vec 
      real :: lengthn, lengthm, lengthl
C

         vh_ini = 0.0d0  

         vh_fin = 0.0d0  

         qn = 0.0d0

         ql = 0.0d0

         qm = 0.0d0

         w = 0.0d0

         deps_N = 0.0d0

         deps_M = 0.0d0

         deps_L = 0.0d0

         eps_N = 0.0d0

         eps_M = 0.0d0

         eps_L = 0.0d0

         sig_N = 0.0d0

         sd0 = 0.0d0
         sn0 = 0.0d0

         ded = 0.0d0

         ed0 = 0.0d0

         edf = 0.0d0

         enf = 0.0d0

         snf = 0.0d0

         snb = 0.0d0

         snb_fib = 0.0d0

         sdf = 0.0d0

         sdpos = 0.0d0 
         sdneg=0.0d0

         eps_N0_pos = 0.0d0
         zeta = 0.d0
         zeta0 = 0.d0
         eps_N0_neg=0.0d0

         sl0 = 0.0d0
         sm0=0.0d0
         slf=0.0d0
         smf=0.0d0

         del = 0.0d0

         dem = 0.0d0

         C_d=0.0d0
         E_N=0.0d0

      te=0.0d0

      ij=reshape((/1,1,2,2,3,3,1,2,2,3,3,1/),(/2,6/))

      if(np.eq.21) then

      te = reshape(
     $(/0.000000000000D+00,0.000000000000D+00,1.000000000000D+00,
     $        2.652142440930D-02,
     $0.000000000000D+00,1.000000000000D+00,0.000000000000D+00,
     $        2.652142440930D-02,
     $1.000000000000D+00,0.000000000000D+00,0.000000000000D+00,
     $        2.652142440930D-02,
     $0.000000000000D+00,7.071067811870D-01,7.071067811870D-01,
     $        1.993014763120D-02,
     $0.000000000000D+00,-7.071067811870D-01,7.071067811870D-01,
     $        1.993014763120D-02,
     $7.071067811870D-01,0.000000000000D+00,7.071067811870D-01,
     $        1.993014763120D-02,
     $-7.071067811870D-01,0.000000000000D+00,7.071067811870D-01,
     $        1.993014763120D-02,
     $7.071067811870D-01,7.071067811870D-01,0.000000000000D+00,
     $        1.993014763120D-02,
     $-7.071067811870D-01,7.071067811870D-01,0.000000000000D+00,
     $        1.993014763120D-02,
     $8.360955967490D-01,3.879073040670D-01,3.879073040670D-01,
     $        2.507123674870D-02,
     $-8.360955967490D-01,3.879073040670D-01,3.879073040670D-01,
     $        2.507123674870D-02,
     $8.360955967490D-01,-3.879073040670D-01,3.879073040670D-01,
     $        2.507123674870D-02,
     $-8.360955967490D-01,-3.879073040670D-01,3.879073040670D-01,
     $        2.507123674870D-02,
     $3.879073040670D-01,8.360955967490D-01,3.879073040670D-01,
     $        2.507123674870D-02,
     $-3.879073040670D-01,8.360955967490D-01,3.879073040670D-01,
     $        2.507123674870D-02,
     $3.879073040670D-01,-8.360955967490D-01,3.879073040670D-01,
     $        2.507123674870D-02,
     $-3.879073040670D-01,-8.360955967490D-01,3.879073040670D-01,
     $        2.507123674870D-02,
     $3.879073040670D-01,3.879073040670D-01,8.360955967490D-01,
     $        2.507123674870D-02,
     $-3.879073040670D-01,3.879073040670D-01,8.360955967490D-01,
     $        2.507123674870D-02,
     $3.879073040670D-01,-3.879073040670D-01,8.360955967490D-01,
     $        2.507123674870D-02,
     $-3.879073040670D-01,-3.879073040670D-01,8.360955967490D-01,
     $        2.507123674870D-02/)
     $,(/4,np/))
      else if(np.eq.28) then

      te = reshape(
     $(/.577350259D0, .577350259D0, .577350259D0, .0160714276D0,
     $.577350259D0, .577350259D0,-.577350259D0, .0160714276D0,
     $.577350259D0,-.577350259D0, .577350259D0, .0160714276D0,
     $.577350259D0,-.577350259D0,-.577350259D0, .0160714276D0,
     $.935113132D0, .250562787D0, .250562787D0, .0204744730D0,
     $.935113132D0, .250562787D0,-.250562787D0, .0204744730D0,
     $.935113132D0,-.250562787D0, .250562787D0, .0204744730D0,
     $.935113132D0,-.250562787D0,-.250562787D0, .0204744730D0,
     $.250562787D0, .935113132D0, .250562787D0, .0204744730D0,
     $.250562787D0, .935113132D0,-.250562787D0, .0204744730D0,
     $.250562787D0,-.935113132D0, .250562787D0, .0204744730D0,
     $.250562787D0,-.935113132D0,-.250562787D0, .0204744730D0,
     $.250562787D0, .250562787D0, .935113132D0, .0204744730D0,
     $.250562787D0, .250562787D0,-.935113132D0, .0204744730D0,
     $.250562787D0,-.250562787D0, .935113132D0, .0204744730D0,
     $.250562787D0,-.250562787D0,-.935113132D0, .0204744730D0,
     $.186156720D0, .694746614D0, .694746614D0, .0158350505D0,
     $.186156720D0, .694746614D0,-.694746614D0, .0158350505D0,
     $.186156720D0,-.694746614D0, .694746614D0, .0158350505D0,
     $.186156720D0,-.694746614D0,-.694746614D0, .0158350505D0,
     $.694746614D0, .186156720D0, .694746614D0, .0158350505D0,
     $.694746614D0, .186156720D0,-.694746614D0, .0158350505D0,
     $.694746614D0,-.186156720D0, .694746614D0, .0158350505D0,
     $.694746614D0,-.186156720D0,-.694746614D0, .0158350505D0,
     $.694746614D0, .694746614D0, .186156720D0, .0158350505D0,
     $.694746614D0, .694746614D0,-.186156720D0, .0158350505D0,
     $.694746614D0,-.694746614D0, .186156720D0, .0158350505D0,
     $.694746614D0,-.694746614D0,-.186156720D0, .0158350505D0/)
     $,(/4,np/))     
      else if(np.eq.37) then

      te = reshape(
     $(/0.000000000000D+00,0.000000000000D+00,1.000000000000D+00,
     $        1.072388573030D-02,
     $0.000000000000D+00,1.000000000000D+00,0.000000000000D+00,
     $        1.072388573030D-02,
     $1.000000000000D+00,0.000000000000D+00,0.000000000000D+00,
     $        1.072388573030D-02,
     $0.000000000000D+00,7.071067811870D-01,7.071067811870D-01,
     $        2.114160951980D-02,
     $0.000000000000D+00,-7.071067811870D-01,7.071067811870D-01,
     $        2.114160951980D-02,
     $7.071067811870D-01,0.000000000000D+00,7.071067811870D-01,
     $        2.114160951980D-02,
     $-7.071067811870D-01,0.000000000000D+00,7.071067811870D-01,
     $        2.114160951980D-02,
     $7.071067811870D-01,7.071067811870D-01,0.000000000000D+00,
     $        2.114160951980D-02,
     $-7.071067811870D-01,7.071067811870D-01,0.000000000000D+00,
     $        2.114160951980D-02,
     $0.000000000000D+00,3.089512677750D-01,9.510778696510D-01,
     $        5.355055908370D-03,
     $0.000000000000D+00,-3.089512677750D-01,9.510778696510D-01,
     $        5.355055908370D-03,
     $0.000000000000D+00,9.510778696510D-01,3.089512677750D-01,
     $        5.355055908370D-03,
     $0.000000000000D+00,-9.510778696510D-01,3.089512677750D-01,
     $        5.355055908370D-03,
     $3.089512677750D-01,0.000000000000D+00,9.510778696510D-01,
     $        5.355055908370D-03,
     $-3.089512677750D-01,0.000000000000D+00,9.510778696510D-01,
     $        5.355055908370D-03,
     $9.510778696510D-01,0.000000000000D+00,3.089512677750D-01,
     $        5.355055908370D-03,
     $-9.510778696510D-01,0.000000000000D+00,3.089512677750D-01,
     $        5.355055908370D-03,
     $3.089512677750D-01,9.510778696510D-01,0.000000000000D+00,
     $        5.355055908370D-03,
     $-3.089512677750D-01,9.510778696510D-01,0.000000000000D+00,
     $        5.355055908370D-03,
     $9.510778696510D-01,3.089512677750D-01,0.000000000000D+00,
     $        5.355055908370D-03,
     $-9.510778696510D-01,3.089512677750D-01,0.000000000000D+00,
     $        5.355055908370D-03,
     $8.805355183100D-01,3.351545919390D-01,3.351545919390D-01,
     $        1.677709091560D-02,
     $-8.805355183100D-01,3.351545919390D-01,3.351545919390D-01,
     $        1.677709091560D-02,
     $8.805355183100D-01,-3.351545919390D-01,3.351545919390D-01,
     $        1.677709091560D-02,
     $-8.805355183100D-01,-3.351545919390D-01,3.351545919390D-01,
     $        1.677709091560D-02,
     $3.351545919390D-01,8.805355183100D-01,3.351545919390D-01,
     $        1.677709091560D-02,
     $-3.351545919390D-01,8.805355183100D-01,3.351545919390D-01,
     $        1.677709091560D-02,
     $3.351545919390D-01,-8.805355183100D-01,3.351545919390D-01,
     $        1.677709091560D-02,
     $-3.351545919390D-01,-8.805355183100D-01,3.351545919390D-01,
     $        1.677709091560D-02,
     $3.351545919390D-01,3.351545919390D-01,8.805355183100D-01,
     $        1.677709091560D-02,
     $-3.351545919390D-01,3.351545919390D-01,8.805355183100D-01,
     $        1.677709091560D-02,
     $3.351545919390D-01,-3.351545919390D-01,8.805355183100D-01,
     $        1.677709091560D-02,
     $-3.351545919390D-01,-3.351545919390D-01,8.805355183100D-01,
     $        1.677709091560D-02,
     $5.773502691900D-01,5.773502691900D-01,5.773502691900D-01,
     $        1.884823095080D-02,
     $-5.773502691900D-01,5.773502691900D-01,5.773502691900D-01,
     $        1.884823095080D-02,
     $5.773502691900D-01,-5.773502691900D-01,5.773502691900D-01,
     $        1.884823095080D-02,
     $-5.773502691900D-01,-5.773502691900D-01,5.773502691900D-01,
     $        1.884823095080D-02/)
     $,(/4,np/))

      end if   

C     ---------------------------------------
C ... Assemble tensors from direction cosines
C     ---------------------------------------
      qn = 0.0d0
      ql = 0.0d0
      qm = 0.0d0
       w = 0.0d0
      call random_seed(size=rs_size)
      do jp=1,np 
        w(jp) = te(4,jp)*6.0D0 
        xn(1) = te(3,jp) 
        xn(2) = te(2,jp) 
        xn(3) = te(1,jp) 

        lengthm = 0.0d0
        do while (lengthm .lt. epsilon(lengthm))
           call random_number(rand_vec)
           xm = rand_vec - dot_product(xn,rand_vec)*xn
           lengthm = sqrt(dot_product(xm,xm))
        end do
        xm = xm/lengthm

        xl(1) = xn(2)*xm(3)-xn(3)*xm(2) 
        xl(2) = xn(3)*xm(1)-xn(1)*xm(3) 
        xl(3) = xn(1)*xm(2)-xn(2)*xm(1) 
        lengthl = sqrt(dot_product(xl,xl))
        xl=xl/lengthl
        lengthn = sqrt(dot_product(xn,xn))
        lengthm = sqrt(dot_product(xm,xm))
        lengthl = sqrt(dot_product(xl,xl))
        do k=1,6 
          i=ij(1,k) 
          j=ij(2,k) 
          qn(k,jp) = xn(i)*xn(j) 
          qm(k,jp) = 0.5D0*(xn(i)*xm(j)+xn(j)*xm(i)) 
          ql(k,jp) = 0.5D0*(xn(i)*xl(j)+xn(j)*xl(i)) 
        end do 
      end do 
      return
      end subroutine setsystem 

C *******************************************************************
C *** SUBROUTINE M7fMATERIAL *****************************************
C *******************************************************************

      subroutine M7fMaterial(
     $     nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
     $     stepTime, totalTime, dt, cmname, coordMp, charLength,
     $     props, density, strainInc, relSpinInc,
     $     tempOld, stretchOld, defgradOld, fieldOld,
     $     stressOld, stateOld, enerInternOld, enerInelasOld,
     $     tempNew, stretchNew, defgradNew, fieldNew,
     $     stressNew, stateNew, enerInternNew, enerInelasNew )

      include 'vaba_param.inc'

      external RateFunc
C
      dimension props(nprops), density(nblock),
     $coordMp(nblock,*),
     $charLength(nblock), strainInc(nblock,ndir+nshr),
     $relSpinInc(nblock,nshr), tempOld(nblock),
     $stretchOld(nblock,ndir+nshr), defgradOld(nblock,ndir+nshr+nshr),
     $fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
     $stateOld(nblock,nstatev), enerInternOld(nblock),
     $enerInelasOld(nblock), tempNew(nblock),
     $stretchNew(nblock,ndir+nshr), defgradNew(nblock,ndir+nshr+nshr), 
     $fieldNew(nblock,nfieldv),
     $stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
     $enerInternNew(nblock), enerInelasNew(nblock)

      character*80 cmname

C
      real :: equivStress, fractureWorkInc, smean, stressPower, two

      real, parameter :: PI=3.1415926535897932384626433832795d0
      integer, parameter ::  np=37, nvhm=5, nvhf=2, nvhi=np*nvhm+nvhf
      common /KCommon_Vars/young, poisson,k_1, C0,dtime,th_del,unit_conv  
      real :: young, poisson, k_1, C0, dtime, th_del, unit_conv, zeta
      real :: sum_bulk, bulk, zeta0
      common /KIFED_M7f/ StrainRateTensorInv2, DevStrainRateTensorInv2,
     $     old_dtime, lruncreep, lfirstpass 
      real :: StrainRateTensorInv2,
     $     DevStrainRateTensorInv2, old_dtime
      common /KM7f_M7fIO/ on_rate_dep_frac, on_creep
      common /KM7f_M7fIO_1/ C_R1, q_1
      common /KM7f_M7fIO_2/ t_0, t_ch, f_c, a_to_c, g_c, w_to_c, E28
      integer      :: on_rate_dep_frac, on_creep
      real :: C_R1, q_1  
      real :: t_0, t_ch, f_c, a_to_c, g_c, w_to_c, E28
      common /KM7f_Microplane_1/ qn, ql, qm, w 
      common /KM7f_Microplane_2/ deps_N, deps_L, deps_M, eps_N,  eps_L,
     $     eps_M,  sig_N
      common /KM7f_Microplane_3/ sd0, sn0, ded, ed0, edf, enf, snf, snb,
     $     sdf
      common /KM7f_Microplane_4/ sl0, sm0, slf, smf, del, dem, vh_ini,
     $     vh_fin, C_d, E_N
      common /KM7f_Microplane_5/ sdneg, sdpos, eps_N0_neg, eps_N0_pos,
     $     snb_fib, Cv, zeta, zeta0
      real, dimension(1:6,1:np) :: qn, ql, qm
      real, dimension(1:np)     :: w
      real, dimension(1:np)     :: deps_N,
     $     deps_L, deps_M, eps_N,  eps_L,  eps_M,  sig_N
      real, dimension(1:np)     :: sd0, sn0,
     $     ded, ed0, edf
      real, dimension(1:np)     :: enf, snf,
     $     snb, sdf, sdneg, sdpos, snb_fib
      real, dimension(1:np)     :: sl0, sm0, slf,
     $     smf, del, dem, eps_N0_neg, eps_N0_pos

      real, dimension(1:nvhi+2)    :: vh_ini, vh_fin 
      real, dimension(1:np)     :: C_d, E_N
      real :: Cv
      save    :: /KCommon_Vars/, /KIFED_M7f/
      save    :: /KM7f_M7fIO/, /KM7f_M7fIO_1/, /KM7f_M7fIO_2/
      save    :: /KM7f_Microplane_1/, /KM7f_Microplane_2/,
     $     /KM7f_Microplane_3/, /KM7f_Microplane_4/, /KM7f_Microplane_5/

      real, dimension(1:6) :: eps, deps, sig_old 
      real, dimension(1:6) :: sig
      real, dimension(1:6) :: deps_s
      real :: R_N, dfb
      real :: sv0, ev0, dev, evf, sum_snf, devv
      real :: svfB, svf, phi0, phi, sig_I_old, eps_I
      integer ::  jp, i 
      real :: svneg, sum1, strain_rate, ek(1:6,1:6)

      NTENS = 6

      if (totalTime<=dt) then
         call inputparams()
         call setsystem() 
      end if

      do i = 1, nblock

         vh_ini = stateOld(i,:) 

         if (totalTime<=dt) then
            vh_ini(2)=1.0d0
         end if

         DO IC = 1, NTENS
            sig(IC)=stressOld(i,IC)
            sig_old(IC)=stressOld(i,IC)
            sig(IC)=sig(IC)*unit_conv
            sig_old(IC)=sig_old(IC)*unit_conv
            IF (IC > 3) THEN    
               deps(IC)=2.0*strainInc(i,IC)
            ELSE
               deps(IC)=strainInc(i,IC)
            END IF
         END DO

         eps(1)=(defgradOld(i,1)*defgradOld(i,1) + 
     $           defgradOld(i,7)*defgradOld(i,7) +
     $           defgradOld(i,6)*defgradOld(i,6) - 1.0)/2.0
         eps(2)=(defgradOld(i,4)*defgradOld(i,4) + 
     $           defgradOld(i,2)*defgradOld(i,2) +
     $           defgradOld(i,8)*defgradOld(i,8) - 1.0)/2.0
         eps(3)=(defgradOld(i,9)*defgradOld(i,9) + 
     $           defgradOld(i,5)*defgradOld(i,5) +
     $           defgradOld(i,3)*defgradOld(i,3) - 1.0)/2.0
         eps(4)=(defgradOld(i,1)*defgradOld(i,4) +
     $           defgradOld(i,2)*defgradOld(i,7) +
     $           defgradOld(i,6)*defgradOld(i,8))/2.0
         eps(5)=(defgradOld(i,9)*defgradOld(i,4) + 
     $           defgradOld(i,2)*defgradOld(i,5) +
     $           defgradOld(i,3)*defgradOld(i,8))/2.0
         eps(6)=(defgradOld(i,1)*defgradOld(i,9) + 
     $           defgradOld(i,7)*defgradOld(i,5) +
     $           defgradOld(i,3)*defgradOld(i,6))/2.0
         DO IC = 1, NTENS
            IF (IC > 3) THEN    ! SHEAR STRAINS ARE SHEAR ANGLES IN MICROPLANE MODEL
               eps(IC)=2.0*eps(IC)
            END IF
         END DO

         deps_s(1)=deps(1)*deps(1)+deps(4)*deps(4)+deps(6)*deps(6)
         deps_s(2)=deps(4)*deps(4)+deps(2)*deps(2)+deps(5)*deps(5)
         deps_s(3)=deps(6)*deps(6)+deps(5)*deps(5)+deps(3)*deps(3)
         deps_s(4)=deps(1)*deps(4)+deps(2)*deps(4)+deps(5)*deps(6)
         deps_s(5)=deps(4)*deps(6)+deps(2)*deps(5)+deps(5)*deps(3)
         deps_s(6)=deps(1)*deps(6)+deps(4)*deps(5)+deps(6)*deps(3)
         sum1=0.d0
         do jp = 1, np
            sum1=sum1+(qn(1,jp)*deps_s(1)+qn(4,jp)*deps_s(4)+
     $           qn(6,jp)*deps_s(6)+qn(4,jp)*deps_s(4)+
     $           qn(2,jp)*deps_s(2)+qn(5,jp)*deps_s(5)+
     $           qn(6,jp)*deps_s(6)+qn(5,jp)*deps_s(5)+
     $           qn(3,jp)*deps_s(3))*w(jp)
         end do
         strain_rate=sqrt(sum1/2.d0)/dt
 
C     -------------------------------------------------------------------
C     Compute material response :
C     Bounding curve formulation for volumetric,
C ... deviatoric and shear variables.
C     Shear strains are shear angles.
C     -------------------------------------------------------------------

C     ---------------
C ... Initializations
C     ---------------

      jp = size(w) 
      
      sig = 0.0d0
       
C     -------------------------------------------------------------------
C ... Calculate the rate coefficients using the rate function
C     -------------------------------------------------------------------
      
      if(on_rate_dep_frac.eq.0) then 
        R_N=0.0d0
      else  

        R_N = RateFunc(strain_rate/C_R1) 

      end if  

C     --------------------------------------------
C ... Evaluate local strains (eps_N, eps_L, eps_M)
C     --------------------------------------------
  
       eps_N =  matmul( eps,qn) 
       eps_L =  matmul( eps,ql)
       eps_M =  matmul( eps,qm)
      deps_N =  matmul(deps,qn)
      deps_L =  matmul(deps,ql)
      deps_M =  matmul(deps,qm)

      if (maxval(eps_N)>th_del) then
         vh_fin(188)=0.
      else
         vh_fin(188)=1.
      end if

C     -----------------------------------------
C ... Volumetric law (same for all microplanes)
C     -----------------------------------------
 
      sv0 = vh_ini(1)
      phi0 = vh_ini(2)
      
      ev0 =  (eps(1)+ eps(2)+ eps(3))/3.0d0
      dev = (deps(1)+deps(2)+deps(3))/3.0d0
      evf = ev0+dev 
      zeta0 = vh_ini(189)
      zeta = zeta0

      call c_vol(dev, ev0, sv0, deps_N, eps_N, R_N, svneg)

C     ---------------------
C ... Normal deviatoric law
C     ---------------------

      sn0 = vh_ini(3:nvhi:nvhm)

      ded = deps_N - dev
      ed0 = eps_N - ev0     
      edf = ed0 + ded

      call c_dev(ded, ed0, dev, ev0, sv0, sd0, sdf, C_d, R_N,
     $     sdneg, sdpos)

      enf = evf + edf 

      eps_N0_pos=vh_ini(6:nvhi:nvhm)
      eps_N0_neg=vh_ini(7:nvhi:nvhm)
      call c_norm_elastic(eps_N, deps_N, sn0, eps_N0_pos,
     $     eps_N0_neg, sv0, snf, E_N, zeta)

      call c_norm(enf, sv0, R_N, snb)
      call c_norm_fib(enf, R_N, snb_fib)

      phi0 = vh_ini(2)

      sig_N = max(min(snf,snb+snb_fib),svneg+sdneg) 

      do ic = 1, jp
         if (snf(ic) > snb(ic)+snb_fib(ic)) then 
               eps_N0_pos(ic) = enf(ic)
         else if (snf(ic) < svneg+sdneg(ic)) then 
            eps_N0_neg(ic) = enf(ic)
         end if
      end do
      
      sum_snf = dot_product(sig_N,w)              
      svfB = sum_snf/3.0d0 

      svf = svfB

      sum_bulk = dot_product(E_N,w)   
      bulk = sum_bulk/3.0d0

      if(sv0 > 0.d0) then
         if(svf > 0.d0) then
            devv = abs(dev - ((svf-sv0)/bulk))
            zeta = zeta0 + abs(devv)

         end if
      end if
C     -----------------
C ... History variables
C     -----------------

      sl0 = vh_ini(4:nvhi:nvhm) 
      sm0 = vh_ini(5:nvhi:nvhm) 
 
      snf = sig_N

      enf = eps_N + deps_N
      del = deps_L
      dem = deps_M

C     ---------
C ... Shear law
C     ---------

      call c_shear2(eps_L, eps_M, C_d, snf, R_N, del, dem, sl0, sm0,
     $     slf, smf, evf)

C     -------------------------------------------------
C     Final Stress Vector
C     -------------------------------------------------
      
      sig = matmul(qn,snf*w) + matmul(qm,smf*w) + matmul(ql,slf*w)

      sig = sig/unit_conv
C
C     -------------------------------------------
C ... Update microplane normal and shear stresses
C     -------------------------------------------
      
      vh_fin(1) = svf

      vh_fin(3:nvhi:nvhm) = snf 
      vh_fin(4:nvhi:nvhm) = slf 
      vh_fin(5:nvhi:nvhm) = smf 
      vh_fin(6:nvhi:nvhm) = eps_N0_pos
      vh_fin(7:nvhi:nvhm) = eps_N0_neg 
      vh_fin(189)=zeta

      stateNew(i,:)=vh_fin
      stressNew(i,:)=sig

      stressPower = (( stressOld(i,1)+stressNew(i,1) )*
     $        strainInc(i,1) + ( stressOld(i,2)+stressNew(i,2) )*
     $        strainInc(i,2) + ( stressOld(i,3)+stressNew(i,3) )*
     $        strainInc(i,3) + 2.0*( stressOld(i,4)+stressNew(i,4) )*
     $        strainInc(i,4) )/2.0

      enerInternNew(i) = enerInternOld(i) + stressPower / density(i)

      smean = ( stressNew(i,1) + stressNew(i,2) +
     $        stressNew(i,3) )/3.0
      equivStress = sqrt( 3.0/2.0 * ( (stressNew(i,1)-smean)**2 +
     $        (stressNew(i,2)-smean)**2 + (stressNew(i,3)-smean)**2 +
     $        two * stressNew(i,4)**2 ) )

      fractureWorkInc = stressNew(i,1)*deps(1) +
     $        stressNew(i,2)*deps(2) + stressNew(i,3)*deps(3) +
     $        stressNew(i,4)*deps(4) + stressNew(i,5)*deps(5) +
     $        stressNew(i,6)*deps(6) 
      enerInelasNew(i) = enerInelasOld(i)+fractureWorkInc/density(i)

      end do

      return
      end subroutine M7fMaterial

C***********************************************************************
C**** This function characterizes the theoretical rate effect      *****
C***********************************************************************
      function RateFunc(arg)
      include 'vaba_param.inc'
      real :: RateFunc
      real :: arg

      RateFunc = log(arg + sqrt(arg * arg + 1.0d0))
      return
      end function RateFunc

C **********************************************************************
C *** SUBROUTINE CREEP *************************************************
C **********************************************************************
      subroutine M7f_Creep(dt)
C+---------------------------------------------------------------------+
C|  Replaces the Elastic Modulus E by the new one Ep which accounts for|
C|  creep effects and stores the coefficient C0                        |
C+---------------------------------------------------------------------+
      include 'vaba_param.inc'
      real, parameter :: PI=3.1415926535897932384626433832795d0
      integer, parameter ::  np=37, nvhm=5, nvhf=2, nvhi=np*nvhm+nvhf
      common /KCommon_Vars/young, poisson,k_1, C0,dtime,th_del,unit_conv  
      real :: young, poisson, k_1, C0, dtime, th_del, unit_conv
      common /KM7f_M7fIO/ on_rate_dep_frac, on_creep
      common /KM7f_M7fIO_1/ C_R1, q_1
      common /KM7f_M7fIO_2/ t_0, t_ch, f_c, a_to_c, g_c, w_to_c, E28
      integer      :: on_rate_dep_frac, on_creep
      real :: C_R1, q_1  
      real :: t_0, t_ch, f_c, a_to_c, g_c, w_to_c, E28
      real :: dt
      real, parameter  :: DaysPerSec = 1.1574074074074074d-5 
      real :: t_age, t_ch_p, density, q1_opt
      real :: q1, q2, q3, q4, n, m, t, J, r, Z, Qf, Q, J_dot
      real :: Ep, eta, tau1, C1

      t_age=t_0*DaysPerSec
      t_ch_p =t_ch*DaysPerSec/2.0d0 
      density=g_c
      q1_opt=q_1

      if(f_c.eq.0.0d0) then
        f_c=E28/4734.0d0*E28/4734.0d0 
      else
        E28=4734.0d0*sqrt(f_c)
      end if

      q2  = 185.4d0*sqrt(density)*f_c**(-0.9d0)
      q4  = 20.3d0*a_to_c**(-0.7d0)
      if(w_to_c.eq.0.0d0) w_to_c=1.0d0/(f_c/22.8d0+0.535d0)

      q1 = 0.60d+6/E28
      q3 = 0.29d0*w_to_c**4*q2
      n=0.1d0
      m=0.5d0
      t=t_age+t_ch_p

      r=1.7d0*t_age**0.12d0+8.0d0
      Z=t_age**(-m)*log(1.0d0+(t-t_age)**n)
      Qf=1.0d0/(0.086d0*t_age**(2.d0/9.d0)+1.21d0*t_age**(4.0d0/9.0d0))
      Q=Qf*Z/(Z**r+Qf**r)**(1.0d0/r)
      J=q1+q2*Q+q3*log(1.0d0+(t-t_age)**n)+q4*log(t/t_age)
      J=q1_opt*(J/q1)

      J_dot=n*(q2*t**(-m)+q3)/(t-t_age+(t-t_age)**(1.0d0-n))+q4/t
      J_dot=q1_opt*(J_dot/q1)
      Ep=1.0d+6/(J-t_ch_p*J_dot)

      eta=1.0d+6/J_dot
      tau1=eta/Ep

      if(dt/tau1.ge.1.0d+2) then
        C0=0.0d0
      else
        C0=exp(-dt/tau1*DaysPerSec)
      end if
      if(dt/tau1.le.1.0d-3) then
        C1=1.0d0 - dt/tau1/2.0d0*DaysPerSec
      else
        C1=(1.0d0-C0)/(dt/tau1*DaysPerSec)
      end if
      Ep=C1*Ep
      young = Ep

      return
      end subroutine M7f_Creep

C **********************************************************************
C *** SUBROUTINE INPUTPARAMS *******************************************
C **********************************************************************
      subroutine inputparams() 
      include 'vaba_param.inc'
      real, parameter :: PI=3.1415926535897932384626433832795d0
      integer, parameter :: np=37, nvhm=5, nvhf=2, nvhi=np*nvhm+nvhf
      common /KCommon_Vars/young, poisson,k_1, C0,dtime,th_del,unit_conv  
      real :: young, poisson, k_1, C0, dtime, th_del, unit_conv
      common /KM7f_M7fIO/ on_rate_dep_frac, on_creep
      common /KM7f_M7fIO_1/ C_R1, q_1
      common /KM7f_M7fIO_2/ t_0, t_ch, f_c, a_to_c, g_c, w_to_c, E28
      integer      :: on_rate_dep_frac, on_creep
      real :: C_R1, q_1  
      real :: t_0, t_ch, f_c, a_to_c, g_c, w_to_c, E28
      common /KM7fBounds_M7fIO/k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9,
     $     sv0_p, c_1, c_2, c_3, c_4, c_5,
     $     c_6, c_7, c_8, c_9, c_10, c_11,c_12,c_13,c_14,c_15,c_16,
     $     c_17,c_18,c_19,c_20,c_21,C_R0,C_R2
      real :: k_2, k_3, k_4, k_5, k_6, k_7, k_8, k_9, sv0_p
      real :: c_1, c_2, c_3, c_4, c_5, c_6, c_7, c_8,
     $     c_9, c_10
      real :: c_11,c_12,c_13,c_14,c_15,c_16,c_17,c_18,
     $     c_19,c_20,c_21,C_R0,C_R2
      common /KM7fBounds_M7fIO_1/ c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1,
     $     c_6_2, c_6_M4, c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4,
     $     c_9_0, c_9_1, c_9_2, c_9_M4
      real :: c_5_0, c_5_1, c_5_M4, c_6_0, c_6_1, c_6_2, c_6_M4
      real :: c_7_0, c_7_1, c_7_M4, c_8_0, c_8_1, c_8_M4
      real :: c_9_0, c_9_1, c_9_2, c_9_M4
      common /KM7fBounds_M7fIO_2/sig_fib_0, p_1, p_2, p_5,p_6,xp_1,xp_2,
     $     xp_3, V_f, cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3, cf_4,
     $     cf_5, cf_4_4, ap1, bp1, cp2, dp2, ep2, d_1, d_2, d_3, d_4,
     $     d_5, d_6
      real :: sig_fib_0, p_1, p_2, p_5, p_6
      real :: xp_1, xp_2, xp_3, V_f
      real :: cf_1, cf_2, cf_3, cf_1_1, cf_2_2, cf_3_3
      real :: cf_4, cf_5, cf_4_4
      real :: ap1, bp1, cp2, dp2, ep2
      real :: d_1, d_2, d_3, d_4, d_5, d_6
      save :: /KCommon_Vars/,/KM7f_M7fIO/,/KM7f_M7fIO_1/,/KM7f_M7fIO_2/,
     $     /KM7fBounds_M7fIO/,/KM7fBounds_M7fIO_1/,/KM7fBounds_M7fIO_2/

      th_del = 0.005d10 

      unit_conv=1.                

      young  = 38300.0d0 
      poisson= 0.18000d0
      k_1    = 150.00d-6 
      k_2    = 110.000d0
      k_3    = 30.d0 
      k_4    = 100.d0 
      c_5_0=1.3d-2
      c_5_1=4.0d0 
      c_7_0  = 1.2d-2 
      c_7_1=35.0d2 
      c_8_0  = 1.2d-2 
      c_8_1=20.0d0 
      c_6_0=4.0d2 
      c_6_1  = 4.0d1 
      c_6_2  = 13.0d0
      c_9_0=4.0d2 
      c_9_1  = 4.0d1 
      c_9_2  = 13.0d0
      c_5_M4=3.0d0    
      c_6_M4 = 1.30d0 
      c_7_M4=10.d-1   
      c_8_M4 = 8.0d0  
      c_9_M4 = 0.d-1 

      d_1=0.095 
      d_2=35.   
      d_3=1.7   
      d_4=1.7d0 
      d_5 = 1.d3
      d_6 = 25.d0
     
      sv0_p=250.d0
      k_8=3.d0
      k_9=5.d-1
      k_5=1.0d0
      k_6=1.0d-4 
      k_7=1.8d0 
     
      c_2 = 1.76d-1 
      c_3 = 4.d0 
      c_4 = 50.d0 
      c_10 = 3.3d-1 
      c_11 = 5.d-1 
      c_12 = 7.00d3 
      c_16 = 10.0d0 
      c_17 = 1.00d-2       
      c_18 = 4.d3 
      c_19 = 4.5d3 
      c_20 = 3.d2 
      c_21 = 6.d1  
     
      sig_fib_0 = -4.5d0 
      p_1    = 0.0d0     
      p_2    = 300.d-3   
      p_5    = 0.0d0     
      p_6    = 0.0d0     
      xp_1   = 3.34d0    
      xp_2   = 0.        
      xp_3   = 1.0d0     
      V_f    = 0.0d0     
      cf_1=104.0d-2
      cf_2=54.0d+1
      cf_3=104.0d-2
      cf_1_1=3.5d+1
      cf_2_2=1.5d+1
      cf_3_3=3.5d+1
      cf_4_4=4.4d2 
      cf_4  = 3.2d-2 
      cf_5  = 5.0d-2 
     
      on_rate_dep_frac = 0
      C_R1 = 4.0d-6
      C_R2 = 22.0d-3
     
      on_creep = 0
      t_0  = 241.920d+4
      t_ch = 180.000d+0
      f_c  =  48.000d+0
      a_to_c = 7.000d+0
      g_c  = 304.340d+0
      w_to_c = 0.000d+0
     
      q_1=21.00d0        
      C0=1.0d0           

      E28 = young
      if(t_ch.eq.0.0d0) then
         write(*,*) 'Error: Duration of loading is not allowed to'//
     $        ' be zero.'
         write(*,*) 'If you want to turn off creep effects, use the '//
     $        'provided switch.'
         stop
      end if
      return
      end subroutine inputparams
