C********************************************************************
C               (STRAIN DRIVEN ALGORITHM)                           *
C      								    *
C  EFFECTIVE ELASTIC BEHAVIOR OF 2DRFPCs CONSIDERING                *
C  INTERFACIAL DAMAGE: STATISTICAL (WEIBUALL) APPROACH              *
C                      						    *
C               HAENG-KI LEE,	LAST EDITED: 7/30/01	            *
C                                                                   *
C********************************************************************
      subroutine umshl53(cm,ldcmw,capa,dt1,time,
     $     deps,ldeps,
     $     sig,ldsig,
     $     hisv,ldhv,
     $     tsh,
     $     ndsh)
      implicit double precision (a-h,o-z)                                   
c
c     user material model for belytschko-tsay shell
c
c     Haeng-Ki Lee
c     Srdan Simunovic
c     Oak Ridge National Laboratory
c     P.O. Box 2008, MS-6359
c     Tel: (865) 574-4837
c     Fax: (865) 574-0680
c     E-mail: leeh@ornl.gov
c
c     Variables
c
c     Micromechanical Damage Model
c
c     deps(1) - local x  strain increment
c     deps(2) - local y  strain increment
c     deps(3) - local z  strain increment
c     deps(4) - local xy strain increment
c     deps(5) - local yz strain increment
c     deps(6) - local zx strain increment
c
c     sig(1) - local x  stress
c     sig(2) - local y  stress
c     sig(3) - local z  stress
c     sig(4) - local xy stress
c     sig(5) - local yz stress
c     sig(6) - local zx stress
c
c     hisv(i) - ith history variable
c
c     dt1 - current time step size
c     capa - reduction factor for transverse shear
c     time - current problem time
c     fail - element failure parameter
c     tsh - time increment for the element
c     ndsh(4) - node array for the element
c
      dimension cm(ldcmw),deps(ldeps),sig(ldsig),hisv(ldhv),ndsh(4)
      dimension dfl(20)
c
c     material parameters
c
c     ee - Young's modulus of composites
c     enu - Poisson's ratio of composites
c     et1 - Young's modulus of the matrix
c     et2 - Young's modulus of fibers
c     anu1 - Poisson's ratio of the matrix
c     anu2 - Poisson's ratio of fibers
c     et3 - Young's modulus of damaged fibers
c     anu3 - Poisson's ratio of damaged fibers
c     sigy - Yiels stress of composites
c     h - Linear isotropic hardening parameter
c     q - Exponential isotropic hardening parameter
c
      ee = cm(1)
      enu = cm(2)
      et1 = cm(3)
      et2 = cm(4)
      anu1 = cm(5)
      anu2 = cm(6)
      et3 = .7*et2
      anu3 = anu2
      sigy = cm(7)
      h = cm(8)
      q = cm(11)
c      
c     s0 - Weibull parameter related to the shape of the function
c     am - Weibull parameter related to the wavelength of the function 
c     ph - Initial volume fraction of fibers
c     a - aspect ratio of fibers
c
      s0 = cm(12)
      am = cm(13)
      ph = cm(14)
      a = cm(15)
c
c --- Compute total strain increments
c
      v1 = deps(1)
      v2 = deps(2)                                                           
      v3 = deps(3)
      v12 = deps(4)                                                           
c
c     v12= 2.0*deps(4) already engineering strain
c
      ph1 = hisv(1)
      ph2 = hisv(2)
      oph2 = hisv(8)
c
c --- Compute new strains                                           
c
      eps1 = hisv(4) + v1
      eps2 = hisv(5) + v2
      eps3 = hisv(6) + v3
      eps12 = hisv(7) + v12
c
      hisv(4) = eps1
      hisv(5) = eps2
      hisv(6) = eps3
      hisv(7) = eps12
c
c --- NEW CONSTITUTIVE MODEL
c
      aka1 = et1/(3. *(1. - 2.*anu1))
      aka2 = et2/(3. *(1. - 2.*anu2))      
      aka3 = et3/(3. *(1. - 2.*anu3))
      amu1 = et1/(2. *(1. + anu1))
      amu2 = et2/(2. *(1. + anu2))
      amu3 = et3/(2. *(1. + anu3))
      alam1 = et1*anu1/((1. + anu1)*(1. -2.*anu1))  
      alam2 = et2*anu2/((1. + anu2)*(1. -2.*anu2))       
      alam3 = et3*anu3/((1. + anu3)*(1. -2.*anu3))       
c
c --- Call input
c
      z11 = 1./3.*(aka1/(aka2 - aka1) - amu1/(amu2 - amu1))
      z12 = 1./3.*(aka1/(aka3 - aka1) - amu1/(amu3 - amu1))
      z21 = amu1/(2.*(amu2 - amu1))
      z22 = amu1/(2.*(amu3 - amu1))
c
      if (a. gt. 1.) then
        g0 = a/(a**2 - 1.)**1.5*(log(a + sqrt(a**2 - 1.)) - 
     x       a*(a**2 - 1.)**.5)
      else 
        g0 = a/(1. - a**2)**1.5*(a*(1. -a**2)**.5 - acos(a))
      endif
c
c --- Eshelby's tensor for spheroidal inclusions 
c        (generalized fourth-rank tensor)
c
      s111 = (4.*anu1 + 2./(a**2 - 1.))*g0 + 4.*anu1 + 
     x       4./(3.*(a**2 - 1.))
      s112 = (4.*anu1 - (2.*a**2 + 1.)/(a**2 - 1.))*g0 + 
     x       4.*anu1 - 2.*a**2/(a**2 - 1.)
      s113 = s112
      s121 = (-2.*anu1 - (2.*a**2 + 1.)/(a**2 - 1.))*g0 - 
     x       2.*a**2/(a**2 - 1.)
      s131 = s121
      s122 = (-2.*anu1 + (4.*a**2 - 1.)/(4.*(a**2 - 1.)))*g0
     x       +a**2/(2.*(a**2 - 1.))
      s123 = s122
      s132 = s122
      s133 = s122
      s211 = (-4.*anu1 + (4.*a**2 - 2.)/(a**2 - 1.))*g0
     x       -4.*anu1 + (12.*a**2 - 8.)/(3.*(a**2 - 1.))
      s212 = (-anu1 - (a**2 + 2.)/(a**2 - 1.))*g0 - 2.*anu1 - 
     x       2./(a**2 - 1.)
      s213 = s212
      s221 = s212
      s231 = s212
      s222 = (2.*anu1 - (4.*a**2 - 7.)/(4.*(a**2 - 1.)))*g0
     x       +a**2/(2.*(a**2 - 1.))
      s223 = s222
      s232 = s222
      s233 = s222
c
c
c --- Eshelby's tensor for spherical inclusions
c       (anisotropic fourth-rank tensor)
c
      si1 = (16. + 45.*g0 + 54.*a**2 + 60.*g0*a**2)/
     x      (16.*(anu1 - 1.)*(1. - a**2))
      si2 = .75*(2./3. + 5./4.*g0 - 2./3.*anu1 - g0*anu1 + 
     x      1./6.*a**2 + 2./3.*anu1*a**2 + g0*anu1*a**2)/
     x      (1. - anu1 - a**2 + anu1**a**2)
      si3 = (3.*g0 + 10.*a**2 + 12.*g0*a**2)/(16.*(anu1 - 1.)*
     x      (a**2 - 1.))
      si4 = (3.*g0 + 16.*anu1 + 24.*g0*anu1 + 10.*a**2 + 12.*g0*a**2
     x      -16.*anu1*a**2 - 24.*g0*anu1*a**2)/(16.*(anu1 - 1.)
     x      *(a**2 - 1.))
      si5 = (g0 - 8.*g0*anu1 - 2.*a**2 - 4.*g0*a**2 + 
     x      8.*g0*anu1*a**2)/(16.*(anu1 - 1.)*(a**2 - 1.))
      si6 = .0625*(-7.*g0 + 8.*g0*anu1 - 2.*a**2 + 4.*g0*a**2-
     x      8.*g0*anu1*a**2)/(1. - anu1 - a**2 + anu1*a**2)
c
c --- Elastic behavior
c
      da11 = si1 + 4.*si2 + si3 + 2.*(.5*amu1/(amu2 - amu1) + si6)
      da12 = si1 + 4.*si2 + si3 + 2.*(.5*amu1/(amu3 - amu1) + si6)
      da21 = si1 + 4.*si2 + 3.*si3
      da22 = ds21
      da31 = si4 + 1./3.*(aka1/(aka2 - aka1) - amu1/(amu2 - amu1)) + 
     x       si5        
      da32 = si4 + 1./3.*(aka1/(aka3 - aka1) - amu1/(amu3 - amu1)) + 
     x       si5
      da41 = si4 + 3.*(1./3.*(aka1/(aka2 - aka1) - amu1/
     x       (amu2 - amu1)) + si5)+ 2.*(.5*amu1/(amu2 - amu1) + si6)
      da42 = si4 + 3.*(1./3.*(aka1/(aka3 - aka1) - amu1/
     x       (amu3 - amu1)) + si5) + 2.*(.5*amu1/(amu3 - amu1) + si6)
c
      dain11 = da41/(da11*da41 - da21*da31)
      dain12 = da42/(da12*da42 - da22*da32)
      dain21 = -da21/(da11*da41 - da21*da31)
      dain22 = -da22/(da12*da42 - da22*da32)
      dain31 = -da31/(da11*da41 - da21*da31)
      dain32 = -da32/(da12*da42 - da22*da32)        
      dain41 = da11/(da11*da41 - da21*da31)
      dain42 = da12/(da12*da42 - da22*da32)
c    
      pp61 = 1./(4.*(.5*amu1/(amu2 - amu1) + si6))        
      pp62 = 1./(4.*(.5*amu1/(amu3 - amu1) + si6))
      pp21 = si2/(-4.*(.5*amu1/(amu2 - amu1) + si6)*(si2 + 
     x       .5*amu1/(amu2 - amu1) + si6))
      pp22 = si2/(-4.*(.5*amu1/(amu3 - amu1) + si6)*(si2 +
     x       .5*amu1/(amu3 - amu1) + si6))
      pp11 = dain11*(2.*pp61*si1 - 4.*pp21*(si1 + 2.*
     x       si2 + si3)) + dain21*(-2.*pp61*si4 - 4.*
     x       pp21 * (si4 + 1./3.*(aka1/(aka2 - aka1)-
     x       amu1/(amu2 - amu1)) + si5))        
      pp12 = dain12*(2.*pp62*si1 - 4.*pp22*(si1 + 2.*
     x       si2 + si3)) + dain22*(-2.*pp62*si4 - 4.*
     x       pp22*(si4 + 1./3.*(aka1/(aka3 - aka1) -
     x       amu1/(amu3 - amu1)) + si5))
      pp41 = dain31*(2.*pp61*si1 - 4.*pp21*(si1 + 2.*
     x       si2 + si3)) + dain41*(-2.*pp61*si4 - 4.*
     x       pp21*(si4 + 1./3.*(aka1/(aka2 - aka1)-
     x       amu1/(amu2 - amu1)) + si5))
      pp42 = dain32*(2.*pp62*si1 - 4.*pp22*(si1 + 2.*
     x       si2 + si3)) + dain42*(-2.*pp62*si4 - 4.*
     x       pp22*(si4 + 1./3.*(aka1/(aka3 - aka1)-
     x       amu1/(amu3 - amu1)) + si5))     
      pp31 = dain11*(-2.*pp61*si3) + dain21*(-2.*pp61*
     x       (1./3.*(aka1/(aka2 - aka1) - amu1/(amu2 - amu1)) + si5))
      pp32 = dain12*(-2.*pp62*si3) + dain22*(-2.*pp62*
     x       (1./3.*(aka1/(aka3 - aka1) - amu1/(amu3 - amu1)) + si5))
      pp51 = dain31*(-2.*pp61*si3) + dain41*(-2.*pp61*
     x       (1./3.*(aka1/(aka2 - aka1) - amu1/(amu2 - amu1)) + si5))
      pp52 = dain32*(-2.*pp62*si3) + dain42*(-2.*pp62*
     x       (1./3.*(aka1/(aka3 - aka1) - amu1/(amu3 - amu1)) + si5))
c
c --- Plastic behavior
c
      if (a .lt. 1.) then
        bb = 5./(2.*(1. - a**4)**2)*(2. + a**4 - 3.*a**4*
     x       acos(a**2)/(a**2*sqrt(1. - a**4)))
        cc = 15.*a**4/(4.*(1.-a**4)**2)*(-3. + (1. + 2.*a**4)*
     x       acos(a**2)/(a**2*sqrt(1. - a**4)))
        dd = 1./8.*(15. - 3.*bb - 4.*cc)
        delt1 = 3. *(1. - a**4*acos(a**2)/(a**2*sqrt(1. - a**4)))/
     x          (1. - a**4)
        delt2 = .5*(3.- delt1)
        delt3 = delt2
        delta11 = bb
        delta12 = cc
        delta13 = cc
        delta21 = cc
        delta22 = dd
        delta23 = dd
        delta31 = cc
        delta32 = dd
        delta33 = dd
      else
        bb = 5./(2.*(1. - a**4)**2)*(2. + a**4 - 3.*a**4
     x       *log(a**2 + sqrt(a**4 - 1.))/(a**2*sqrt(a**4 - 1.)))
        cc = 15.*a**4/(4.*(1. - a**4)**2)*(-3. + (1. + 2.*a**4)
     x       *log(a**2 + sqrt(a**4 - 1.))/(a**2*sqrt(a**4 - 1.)))
        dd = 1./8.*(15. - 3.*bb - 4.*cc)
        delt1 = 3.*(1. - a**4*log(a**2 + sqrt(a**4 - 1.))
     x          /(a**2*sqrt(a**4 - 1.)))/(1. - a**4)
        delt2 = .5*(3. - delt1)
        delt3 = delt2
        delta11 = bb
        delta12 = cc
        delta13 = cc
        delta21 = cc
        delta22 = dd
        delta23 = dd
        delta31 = cc
        delta32 = dd
        delta33 = dd            
      endif
c
      b111 = 2.*(z21 + s211)
      b112 = 2.*(z21 + s212)
      b113 = 2.*(z21 + s213)
      b121 = 2.*(z21 + s221)
      b122 = 2.*(z21 + s222)
      b123 = 2.*(z21 + s223)
      b131 = 2.*(z21 + s231)
      b132 = 2.*(z21 + s232)
      b133 = 2.*(z21 + s233)
      b211 = 2.*(z22 + s211)
      b212 = 2.*(z22 + s212)
      b213 = 2.*(z22 + s213)
      b221 = 2.*(z22 + s221)
      b222 = 2.*(z22 + s222)
      b223 = 2.*(z22 + s223)
      b231 = 2.*(z22 + s231)
      b232 = 2.*(z22 + s232)
      b233 = 2.*(z22 + s233)     
c 
      gamm111 = ((2.*z11 + 2.*s122 + b122)*(z11 + s111)
     x          -2.*(z11 + s121)*(z11 + s112))/((2.*z11 + 2.
     x          *s122 + b122)*(z11 + s111 + b111) - 2.*
     x          (z11 + s112)*(z11 + s121))
      gamm121 = ((2.*z11 + 2.*s122 + b122)*(z11 + s121)
     x          -2.*(z11 + s121)*(z11 + s122))/((2.*z11 + 2.
     x          *s122 + b122)*(z11 + s111 + b111) - 2.*
     x          (z11 + s112)*(z11 + s121))
      gamm131 = ((2.*z11 + 2.*s122 + b122)*(z11 + s131)
     x          -2.*(z11 + s121)*(z11 + s132))/((2.*z11 + 2.
     x          *s122 + b122)*(z11 + s111 + b111) - 2.*
     x          (z11 + s112)*(z11 + s121))
      gamm211 = ((2.*z12 + 2.*s122 + b222)*(z12 + s111)
     x          -2.*(z12 + s121)*(z12 + s112))/((2.*z12 + 2.
     x          *s122 + b222)*(z12 + s111 + b211) - 2.*
     x          (z12 + s112)*(z12 + s121))
      gamm221 = ((2.*z12 + 2.*s122 + b222)*(z12 + s121)
     x          -2.*(z12 + s121)*(z12 + s122))/((2.*z12 + 2.
     x          *s122 + b222)*(z12 + s111 + b211) - 2.*
     x          (z12 + s112)*(z12 + s121))     
      gamm231 = ((2.*z12 + 2.*s122 + b222)*(z12 + s131)
     x          -2.*(z12 + s121)*(z12 + s132))/((2.*z12 + 2.
     x          *s122 + b222)*(z12 + s111 + b211) - 2.*
     x          (z12 + s112)*(z12 + s121))
      gamm112 = ((z11 + s111 + b111)*(z11 + s112) - (z11 + s112)*
     x          (z11 + s111))/((2.*z11 + 2.*s122 + b122)*
     x          (z11 + s111 + b111) - 2.*(z11 + s112)*(z11 + s121))
      gamm122 = ((z11 + s111 + b111)*(z11 + s122) - (z11 + s112)*
     x          (z11 + s121))/((2.*z11 + 2.*s122 + b122)*
     x          (z11 + s111 + b111) - 2.*(z11 + s112)*(z11 + s121))
      gamm132 = ((z11 + s111 + b111)*(z11 + s132) - (z11 + s112)*
     x          (z11 + s131))/((2.*z11 + 2.*s122 + b122)*
     x          (z11 + s111 + b111) - 2.*(z11 + s112)*(z11 + s121))          
      gamm212 = ((z12 + s111 + b211)*(z12 + s112) - (z12 + s112)*
     x          (z12 + s111))/((2.*z12 + 2.*s122 + b222)*
     x          (z12 + s111 + b211) - 2.*(z12 + s112)*(z12 + s121))  
      gamm222 = ((z12 + s111 + b211)*(z12 + s122) - (z12 + s112)*
     x          (z12 + s121))/((2.*z12 + 2.*s122 + b222)*
     x          (z12 + s111 + b211) - 2.*(z12 + s112)*(z12 + s121))   
      gamm232 = ((z12 + s111 + b211)*(z12 + s132) - (z12 + s112)*
     x          (z12 + s131))/((2.*z12 + 2.*s122 + b222)*
     x          (z12 + s111 + b211) - 2.*(z12 + s112)*(z12 + s121))                
      gamm113 = gamm112
      gamm123 = gamm122
      gamm133 = gamm132
      gamm213 = gamm212
      gamm223 = gamm222
      gamm233 = gamm232
c
c --- Internal stresses of fibers
c
      aj11 = si1*(pp11 + 4.*pp21 + pp31 + 2.*pp61) +
     x       4.*si2*(pp11 + 2.*pp21 + pp31) + si4*
     x       (pp11 + 4.*pp21 + 3.*pp31) + 2.*si6*pp11
      aj12 = si1*(pp12 + 4.*pp22 + pp32 + 2.*pp62) +
     x       4.*si2*(pp12 + 2.*pp22 + pp32) + si4*
     x       (pp12 + 4.*pp22 + 3.*pp32) + 2.*si6*pp12
      aj21 = 2.*si2*(pp21 + pp61) + 2.*si6*pp21
      aj22 = 2.*si2*(pp22 + pp62) + 2.*si6*pp22
      aj31 = si3*(pp11 + 4.*pp21 + pp31 + 2.*pp61) + 
     x       si5*(pp11 + 4.*pp21 + 3.*pp31) + 2.*si6*pp31
      aj32 = si3*(pp12 + 4.*pp22 + pp32 + 2.*pp62) + 
     x       si5*(pp12 + 4.*pp22 + 3.*pp32) + 2.*si6*pp32
      aj41 = si1*(pp41 + pp51) + 4.*si2*(pp41 + pp51) + 
     x       si4*(pp41 + 3.*pp51 + 2.*pp61) + 2.*si6*pp41
      aj42 = si1*(pp42 + pp52) + 4.*si2*(pp42 + pp52) + 
     x       si4*(pp42 + 3.*pp52 + 2.*pp62) + 2.*si6*pp42
      aj51 = si3*(pp41 + pp51) + si5*(pp41 + 3.*pp51 + 2.*pp61) + 
     x       2.*si6*pp51
      aj52 = si3*(pp42 + pp52) + si5*(pp42 + 3.*pp52 + 2.*pp62) + 
     x       2.*si6*pp52
      aj61 = 2.*si6*pp61
      aj62 = 2.*si6*pp62
c       
      ajj11 = -aj11
      ajj12 = -aj12
      ajj21 = -aj21
      ajj22 = -aj22
      ajj31 = -aj31
      ajj32 = -aj32
      ajj41 = -aj41 
      ajj42 = -aj42
      ajj51 = -aj51
      ajj52 = -aj52
      ajj61 = .5 - aj61    
      ajj62 = .5 - aj62   
c       
c --- Call avgt
c
      t111 = -1./3.+2.*ph1/(4725.*(1.-anu1**2)*b111*b111)*(1575.*
     x     (1.-2.*anu1)**2*gamm111*gamm111+21.*(25.*anu1-23.)*(1.-2.*
     x     anu1)*(gamm111*delt1+gamm111*delt1)+21.*(25.*anu1-2.)*
     x     (1.-2.*anu1)*(gamm111+gamm111)+3.*(35.*anu1**2-70.*anu1+
     x     36.)*delta11+7.*(50.*anu1**2-59.*anu1+8.)*(delt1+
     x     delt1)-2.*(175.*anu1**2-343.*anu1+103.))+2.*ph2/(4725.*
     x     (1.-anu1**2)*b211*b211)*(1575.*(1.-2.*anu1)**2*gamm211*
     x     gamm211+21.*(25.*anu1-23.)*(1.-2.*anu1)*(gamm211*delt1+
     x     gamm211*delt1)+21.*(25.*anu1-2.)*(1.-2.*anu1)*(gamm211+
     x     gamm211)+3.*(35.*anu1**2-70.*anu1+36.)*delta11+7.*(50.*
     x     anu1**2-59.*anu1+8.)*(delt1+delt1)-2.*(175.*anu1**2-
     x     343.*anu1+103.))
      t211 = .5+ph1/(1575.*(1.-anu1)**2*b111*b111)*((72.-140.*anu1+
     x     70.*anu1**2)*delta11-(75.-266.*anu1+175.*anu1**2)*
     x     (delt1+delt1)/2.+164.-476.*anu1+350.*anu1**2)+ph2/
     x     (1575.*(1.-anu1)**2*b211*b211)*((72.-140.*anu1+70.*anu1**2)*
     x     delta11-(75.-266.*anu1+175.*anu1**2)*(delt1+delt1)/2.+
     x     164.-476.*anu1+350.*anu1**2)
c
      t112 = -1./3.+2.*ph1/(4725.*(1.-anu1**2)*b111*b122)*(1575.*
     x     (1.-2.*anu1)**2*gamm111*gamm122+21.*(25.*anu1-23.)*(1.-2.*
     x     anu1)*(gamm111*delt2+gamm122*delt1)+21.*(25.*anu1-2.)*
     x     (1.-2.*anu1)*(gamm111+gamm122)+3.*(35.*anu1**2-70.*anu1+
     x     36.)*delta12+7.*(50.*anu1**2-59.*anu1+8.)*(delt1+
     x     delt2)-2.*(175.*anu1**2-343.*anu1+103.))+2.*ph2/(4725.*
     x     (1.-anu1**2)*b211*b222)*(1575.*(1.-2.*anu1)**2*gamm211*
     x     gamm222+21.*(25.*anu1-23.)*(1.-2.*anu1)*(gamm211*delt2+
     x     gamm222*delt1)+21.*(25.*anu1-2.)*(1.-2.*anu1)*(gamm211+
     x     gamm222)+3.*(35.*anu1**2-70.*anu1+36.)*delta12+7.*(50.*
     x     anu1**2-59.*anu1+8.)*(delt1+delt2)-2.*(175.*anu1**2-
     x     343.*anu1+103.))
      t212 = .5+ph1/(1575.*(1.-anu1)**2*b112*b112)*((72.-140.*anu1+
     x     70.*anu1**2)*delta12-(75.-266.*anu1+175.*anu1**2)*
     x     (delt1+delt2)/2.+164.-476.*anu1+350.*anu1**2)+ph2/
     x     (1575.*(1.-anu1)**2*b212*b212)*((72.-140.*anu1+70.*anu1**2)*
     x     delta12-(75.-266.*anu1+175.*anu1**2)*(delt1+delt2)/2.+
     x     164.-476.*anu1+350.*anu1**2)
c
      t113 = -1./3.+2.*ph1/(4725.*(1.-anu1**2)*b111*b133)*(1575.*
     x     (1.-2.*anu1)**2*gamm111*gamm133+21.*(25.*anu1-23.)*(1.-2.*
     x     anu1)*(gamm111*delt3+gamm133*delt1)+21.*(25.*anu1-2.)*
     x     (1.-2.*anu1)*(gamm111+gamm133)+3.*(35.*anu1**2-70.*anu1+
     x     36.)*delta13+7.*(50.*anu1**2-59.*anu1+8.)*(delt1+
     x     delt3)-2.*(175.*anu1**2-343.*anu1+103.))+2.*ph2/(4725.*
     x     (1.-anu1**2)*b211*b233)*(1575.*(1.-2.*anu1)**2*gamm211*
     x     gamm233+21.*(25.*anu1-23.)*(1.-2.*anu1)*(gamm211*delt3+
     x     gamm233*delt1)+21.*(25.*anu1-2.)*(1.-2.*anu1)*(gamm211+
     x     gamm233)+3.*(35.*anu1**2-70.*anu1+36.)*delta13+7.*(50.*
     x     anu1**2-59.*anu1+8.)*(delt1+delt3)-2.*(175.*anu1**2-
     x     343.*anu1+103.))
      t213 = .5+ph1/(1575.*(1.-anu1)**2*b113*b113)*((72.-140.*anu1+
     x     70.*anu1**2)*delta13-(75.-266.*anu1+175.*anu1**2)*
     x     (delt1+delt3)/2.+164.-476.*anu1+350.*anu1**2)+ph2/
     x     (1575.*(1.-anu1)**2*b213*b213)*((72.-140.*anu1+70.*anu1**2)*
     x     delta13-(75.-266.*anu1+175.*anu1**2)*(delt1+delt3)/2.+
     x     164.-476.*anu1+350.*anu1**2)
c
      t121 = -1./3.+2.*ph1/(4725.*(1.-anu1**2)*b122*b111)*(1575.*
     x     (1.-2.*anu1)**2*gamm122*gamm111+21.*(25.*anu1-23.)*(1.-2.*
     x     anu1)*(gamm122*delt1+gamm111*delt2)+21.*(25.*anu1-2.)*
     x     (1.-2.*anu1)*(gamm122+gamm111)+3.*(35.*anu1**2-70.*anu1+
     x     36.)*delta21+7.*(50.*anu1**2-59.*anu1+8.)*(delt2+
     x     delt1)-2.*(175.*anu1**2-343.*anu1+103.))+2.*ph2/(4725.*
     x     (1.-anu1**2)*b222*b211)*(1575.*(1.-2.*anu1)**2*gamm222*
     x     gamm211+21.*(25.*anu1-23.)*(1.-2.*anu1)*(gamm222*delt1+
     x     gamm211*delt2)+21.*(25.*anu1-2.)*(1.-2.*anu1)*(gamm222+
     x     gamm211)+3.*(35.*anu1**2-70.*anu1+36.)*delta21+7.*(50.*
     x     anu1**2-59.*anu1+8.)*(delt2+delt1)-2.*(175.*anu1**2-
     x     343.*anu1+103.))
      t221 = .5+ph1/(1575.*(1.-anu1)**2*b121*b121)*((72.-140.*anu1+
     x     70.*anu1**2)*delta21-(75.-266.*anu1+175.*anu1**2)*
     x     (delt2+delt1)/2.+164.-476.*anu1+350.*anu1**2)+ph2/
     x     (1575.*(1.-anu1)**2*b221*b221)*((72.-140.*anu1+70.*anu1**2)*
     x     delta21-(75.-266.*anu1+175.*anu1**2)*(delt2+delt1)/2.+
     x     164.-476.*anu1+350.*anu1**2)
c
      t122 = -1./3.+2.*ph1/(4725.*(1.-anu1**2)*b122*b122)*(1575.*
     x     (1.-2.*anu1)**2*gamm122*gamm122+21.*(25.*anu1-23.)*(1.-2.*
     x     anu1)*(gamm122*delt2+gamm122*delt2)+21.*(25.*anu1-2.)*
     x     (1.-2.*anu1)*(gamm122+gamm122)+3.*(35.*anu1**2-70.*anu1+
     x     36.)*delta22+7.*(50.*anu1**2-59.*anu1+8.)*(delt2+
     x     delt2)-2.*(175.*anu1**2-343.*anu1+103.))+2.*ph2/(4725.*
     x     (1.-anu1**2)*b222*b222)*(1575.*(1.-2.*anu1)**2*gamm222*
     x     gamm222+21.*(25.*anu1-23.)*(1.-2.*anu1)*(gamm222*delt2+
     x     gamm222*delt2)+21.*(25.*anu1-2.)*(1.-2.*anu1)*(gamm222+
     x     gamm222)+3.*(35.*anu1**2-70.*anu1+36.)*delta22+7.*(50.*
     x     anu1**2-59.*anu1+8.)*(delt2+delt2)-2.*(175.*anu1**2-
     x     343.*anu1+103.))
      t222 = .5+ph1/(1575.*(1.-anu1)**2*b122*b122)*((72.-140.*anu1+
     x     70.*anu1**2)*delta22-(75.-266.*anu1+175.*anu1**2)*
     x     (delt2+delt2)/2.+164.-476.*anu1+350.*anu1**2)+ph2/
     x     (1575.*(1.-anu1)**2*b222*b222)*((72.-140.*anu1+70.*anu1**2)*
     x     delta22-(75.-266.*anu1+175.*anu1**2)*(delt2+delt2)/2.+
     x     164.-476.*anu1+350.*anu1**2)
c
      t123 = -1./3.+2.*ph1/(4725.*(1.-anu1**2)*b122*b133)*(1575.*
     x     (1.-2.*anu1)**2*gamm122*gamm133+21.*(25.*anu1-23.)*(1.-2.*
     x     anu1)*(gamm122*delt3+gamm133*delt2)+21.*(25.*anu1-2.)*
     x     (1.-2.*anu1)*(gamm122+gamm133)+3.*(35.*anu1**2-70.*anu1+
     x     36.)*delta23+7.*(50.*anu1**2-59.*anu1+8.)*(delt2+
     x     delt3)-2.*(175.*anu1**2-343.*anu1+103.))+2.*ph2/(4725.*
     x     (1.-anu1**2)*b222*b233)*(1575.*(1.-2.*anu1)**2*gamm222*
     x     gamm233+21.*(25.*anu1-23.)*(1.-2.*anu1)*(gamm222*delt3+
     x     gamm233*delt2)+21.*(25.*anu1-2.)*(1.-2.*anu1)*(gamm222+
     x     gamm233)+3.*(35.*anu1**2-70.*anu1+36.)*delta23+7.*(50.*
     x     anu1**2-59.*anu1+8.)*(delt2+delt3)-2.*(175.*anu1**2-
     x     343.*anu1+103.))
      t223 = .5+ph1/(1575.*(1.-anu1)**2*b123*b123)*((72.-140.*anu1+
     x     70.*anu1**2)*delta23-(75.-266.*anu1+175.*anu1**2)*
     x     (delt2+delt3)/2.+164.-476.*anu1+350.*anu1**2)+ph2/
     x     (1575.*(1.-anu1)**2*b223*b223)*((72.-140.*anu1+70.*anu1**2)*
     x     delta23-(75.-266.*anu1+175.*anu1**2)*(delt2+delt3)/2.+
     x     164.-476.*anu1+350.*anu1**2)
c 
      t131 = -1./3.+2.*ph1/(4725.*(1.-anu1**2)*b133*b111)*(1575.*
     x     (1.-2.*anu1)**2*gamm133*gamm111+21.*(25.*anu1-23.)*(1.-2.*
     x     anu1)*(gamm133*delt1+gamm111*delt3)+21.*(25.*anu1-2.)*
     x     (1.-2.*anu1)*(gamm133+gamm111)+3.*(35.*anu1**2-70.*anu1+
     x     36.)*delta31+7.*(50.*anu1**2-59.*anu1+8.)*(delt1+
     x     delt1)-2.*(175.*anu1**2-343.*anu1+103.))+2.*ph2/(4725.*
     x     (1.-anu1**2)*b233*b211)*(1575.*(1.-2.*anu1)**2*gamm233*
     x     gamm211+21.*(25.*anu1-23.)*(1.-2.*anu1)*(gamm233*delt1+
     x     gamm211*delt1)+21.*(25.*anu1-2.)*(1.-2.*anu1)*(gamm233+
     x     gamm211)+3.*(35.*anu1**2-70.*anu1+36.)*delta31+7.*(50.*
     x     anu1**2-59.*anu1+8.)*(delt3+delt1)-2.*(175.*anu1**2-
     x     343.*anu1+103.))
      t231 = .5+ph1/(1575.*(1.-anu1)**2*b131*b131)*((72.-140.*anu1+
     x     70.*anu1**2)*delta31-(75.-266.*anu1+175.*anu1**2)*
     x     (delt3+delt1)/2.+164.-476.*anu1+350.*anu1**2)+ph2/
     x     (1575.*(1.-anu1)**2*b231*b231)*((72.-140.*anu1+70.*anu1**2)*
     x     delta31-(75.-266.*anu1+175.*anu1**2)*(delt3+delt1)/2.+
     x     164.-476.*anu1+350.*anu1**2)
c
      t132 = -1./3.+2.*ph1/(4725.*(1.-anu1**2)*b133*b122)*(1575.*
     x     (1.-2.*anu1)**2*gamm133*gamm122+21.*(25.*anu1-23.)*(1.-2.*
     x     anu1)*(gamm133*delt2+gamm122*delt3)+21.*(25.*anu1-2.)*
     x     (1.-2.*anu1)*(gamm133+gamm122)+3.*(35.*anu1**2-70.*anu1+
     x     36.)*delta32+7.*(50.*anu1**2-59.*anu1+8.)*(delt3+
     x     delt2)-2.*(175.*anu1**2-343.*anu1+103.))+2.*ph2/(4725.*
     x     (1.-anu1**2)*b233*b222)*(1575.*(1.-2.*anu1)**2*gamm233*
     x     gamm222+21.*(25.*anu1-23.)*(1.-2.*anu1)*(gamm233*delt2+
     x     gamm222*delt3)+21.*(25.*anu1-2.)*(1.-2.*anu1)*(gamm233+
     x     gamm222)+3.*(35.*anu1**2-70.*anu1+36.)*delta32+7.*(50.*
     x     anu1**2-59.*anu1+8.)*(delt3+delt2)-2.*(175.*anu1**2-
     x     343.*anu1+103.))
      t232 = .5+ph1/(1575.*(1.-anu1)**2*b132*b132)*((72.-140.*anu1+
     x     70.*anu1**2)*delta32-(75.-266.*anu1+175.*anu1**2)*
     x     (delt3+delt2)/2.+164.-476.*anu1+350.*anu1**2)+ph2/
     x     (1575.*(1.-anu1)**2*b232*b232)*((72.-140.*anu1+70.*anu1**2)*
     x     delta32-(75.-266.*anu1+175.*anu1**2)*(delt3+delt2)/2.+
     x     164.-476.*anu1+350.*anu1**2)
c
      t133 = -1./3.+2.*ph1/(4725.*(1.-anu1**2)*b133*b133)*(1575.*
     x     (1.-2.*anu1)**2*gamm133*gamm133+21.*(25.*anu1-23.)*(1.-2.*
     x     anu1)*(gamm133*delt3+gamm133*delt3)+21.*(25.*anu1-2.)*
     x     (1.-2.*anu1)*(gamm133+gamm133)+3.*(35.*anu1**2-70.*anu1+
     x     36.)*delta33+7.*(50.*anu1**2-59.*anu1+8.)*(delt3+
     x     delt3)-2.*(175.*anu1**2-343.*anu1+103.))+2.*ph2/(4725.*
     x     (1.-anu1**2)*b233*b233)*(1575.*(1.-2.*anu1)**2*gamm233*
     x     gamm233+21.*(25.*anu1-23.)*(1.-2.*anu1)*(gamm233*delt3+
     x     gamm233*delt3)+21.*(25.*anu1-2.)*(1.-2.*anu1)*(gamm233+
     x     gamm233)+3.*(35.*anu1**2-70.*anu1+36.)*delta33+7.*(50.*
     x     anu1**2-59.*anu1+8.)*(delt3+delt3)-2.*(175.*anu1**2-
     x     343.*anu1+103.))
      t233 = .5+ph1/(1575.*(1.-anu1)**2*b133*b133)*((72.-140.*anu1+
     x     70.*anu1**2)*delta33-(75.-266.*anu1+175.*anu1**2)*
     x     (delt3+delt3)/2.+164.-476.*anu1+350.*anu1**2)+ph2/
     x     (1575.*(1.-anu1)**2*b233*b233)*((72.-140.*anu1+70.*anu1**2)*
     x     delta33-(75.-266.*anu1+175.*anu1**2)*(delt3+delt3)/2.+
     x     164.-476.*anu1+350.*anu1**2)   
c
      bt1 = t111 + 2.*t211 - t112 - t121 + t122 + 2.*t222 - 4.*t231
      bt2 = -.5*(t122 + 2.*t222 - t123) + t231
      bt3 = t121 - t123
      bt4 = t112 - t123
      bt5 = t123
      bt6 = .5*(t122 + 2.*t222 - t123)
c          
      avgt1 = 3./8.*bt1
      avgt2 = -1./8.*(bt1 + 4.*bt2)
      avgt3 = -1./8.*(bt1 + 4.*bt3)
      avgt4 = -1./8.*(bt1 + 4.*bt4)
      avgt5 = 1./8.*(bt1 + 4.*(bt3 + bt4 + 2.*bt5))
      avgt6 = 1./8.*bt1 + bt2 + bt6
c
      gam111 = ph1*(gamm111*s111/b111+gamm121*s112/b122+gamm131*s113/
     x             b133-s111/b111-gamm111*(1.-2.*s211)/b111)+
     x         ph2*(gamm211*s111/b211+gamm221*s112/b222+gamm231*s113/
     x             b233-s111/b211-gamm211*(1.-2.*s211)/b211)
      gam211 = .5+ph1*(1.-2.*s211)/(2.*b111)+ph2*(1.-2.*s211)/(2.*b211)
c
      gam112 = ph1*(gamm112*s111/b111+gamm122*s112/b122+gamm132*s113/
     x             b133-s112/b122-gamm112*(1.-2.*s211)/b111)+
     x         ph2*(gamm212*s111/b211+gamm222*s112/b222+gamm232*s113/
     x             b233-s112/b222-gamm212*(1.-2.*s211)/b211)
      gam212 = .5+ph1*(1.-2.*s212)/(2.*b112)+ph2*(1.-2.*s212)/(2.*b212)
c
      gam113 = ph1*(gamm113*s111/b111+gamm123*s112/b122+gamm133*s113/
     x             b133-s113/b133-gamm113*(1.-2.*s211)/b111)+
     x         ph2*(gamm213*s111/b211+gamm223*s112/b222+gamm233*s113/
     x             b233-s113/b233-gamm213*(1.-2.*s211)/b211)
      gam213 = .5+ph1*(1.-2.*s213)/(2.*b113)+ph2*(1.-2.*s213)/(2.*b213)
c
      gam121 = ph1*(gamm111*s121/b111+gamm121*s122/b122+gamm131*s123/
     x             b133-s121/b111-gamm121*(1.-2.*s222)/b122)+
     x         ph2*(gamm211*s121/b211+gamm221*s122/b222+gamm231*s123/
     x             b233-s121/b211-gamm221*(1.-2.*s222)/b222)
      gam221 = .5+ph1*(1.-2.*s221)/(2.*b121)+ph2*(1.-2.*s221)/(2.*b221)
c
      gam122 = ph1*(gamm112*s121/b111+gamm122*s122/b122+gamm132*s123/
     x             b133-s122/b122-gamm122*(1.-2.*s222)/b122)+
     x         ph2*(gamm212*s121/b211+gamm222*s122/b222+gamm232*s123/
     x             b233-s122/b222-gamm222*(1.-2.*s222)/b222)
      gam222 = .5+ph1*(1.-2.*s222)/(2.*b122)+ph2*(1.-2.*s222)/(2.*b222)
c      
      gam123 = ph1*(gamm113*s121/b111+gamm123*s122/b122+gamm133*s123/
     x             b133-s123/b133-gamm123*(1.-2.*s222)/b122)+
     x         ph2*(gamm213*s121/b211+gamm223*s122/b222+gamm233*s123/
     x             b233-s123/b233-gamm223*(1.-2.*s222)/b222)
      gam223 = .5+ph1*(1.-2.*s223)/(2.*b123)+ph2*(1.-2.*s223)/(2.*b223)
c      
      gam131 = ph1*(gamm111*s131/b111+gamm121*s132/b122+gamm131*s133/
     x             b133-s131/b111-gamm131*(1.-2.*s233)/b133)+
     x         ph2*(gamm211*s131/b211+gamm221*s132/b222+gamm231*s133/
     x             b233-s131/b211-gamm231*(1.-2.*s233)/b233)
      gam231 = .5+ph1*(1.-2.*s231)/(2.*b131)+ph2*(1.-2.*s231)/(2.*b231)
c      
      gam132 = ph1*(gamm112*s131/b111+gamm122*s132/b122+gamm132*s133/
     x             b133-s132/b122-gamm132*(1.-2.*s233)/b133)+
     x         ph2*(gamm212*s131/b211+gamm222*s132/b222+gamm232*s133/
     x             b233-s132/b222-gamm232*(1.-2.*s233)/b233)
      gam232 = .5+ph1*(1.-2.*s232)/(2.*b132)+ph2*(1.-2.*s232)/(2.*b232)
c      
      gam133 = ph1*(gamm113*s131/b111+gamm123*s132/b122+gamm133*s133/
     x             b133-s133/b133-gamm133*(1.-2.*s233)/b133)+
     x         ph2*(gamm213*s131/b211+gamm223*s132/b222+gamm233*s133/
     x             b233-s133/b233-gamm233*(1.-2.*s233)/b233)
      gam233 = .5+ph1*(1.-2.*s233)/(2.*b133)+ph2*(1.-2.*s233)/(2.*b233)
c
      gam311 = ((gam122+gam222)*gam111-gam121*gam112)/((gam111+
     x         2.*gam211)*(gam122+gam222)-gam112*gam121)        
      gam312 = ((gam111+2.*gam211)*gam112-gam112*gam111)/(2.*(gam111+
     x         2.*gam211)*(gam122+gam222)-2.*gam112*gam121)
      gam313 = gam312
c          
      gam321 = ((gam122+gam222)*gam121-gam121*gam122)/((gam111+
     x         2.*gam211)*(gam122+gam222)-gam112*gam121)        
      gam322 = ((gam111+2.*gam211)*gam122-gam112*gam121)/(2.*(gam111+
     x         2.*gam211)*(gam122+gam222)-2.*gam112*gam121)
      gam323 = gam322
c 
      gam331 = ((gam122+gam222)*gam131-gam121*gam132)/((gam111+
     x         2.*gam211)*(gam122+gam222)-gam112*gam121)        
      gam332 = ((gam111+2.*gam211)*gam132-gam112*gam131)/(2.*(gam111+
     x         2.*gam211)*(gam122+gam222)-2.*gam112*gam121)
      gam333 = gam332
c  
      tp111 = -gam311/(2.*gam211)
      tp211 = 1./(4.*gam211)
      tp112 = -gam312/(2.*gam212)
      tp212 = 1./(4.*gam212)                     
      tp113 = -gam313/(2.*gam213)
      tp213 = 1./(4.*gam231)     
      tp121 = -gam321/(2.*gam221)
      tp221 = 1./(4.*gam221)
      tp122 = -gam322/(2.*gam222)
      tp222 = 1./(4.*gam222)                     
      tp123 = -gam323/(2.*gam223)
      tp223 = 1./(4.*gam223)     
      tp131 = -gam331/(2.*gam231)
      tp231 = 1./(4.*gam231)
      tp132 = -gam332/(2.*gam232)
      tp232 = 1./(4.*gam232)                     
      tp133 = -gam333/(2.*gam233)
      tp233 = 1./(4.*gam233)     
c
      bp1 = tp111 + 2.*tp211 - tp112 - tp121 + tp122 + 2.*tp222 - 
     x      4.*tp231
      bp2 = -.5*(tp122 + 2.*tp222 - tp123) + tp231
      bp3 = tp121 - tp123
      bp4 = tp112 - tp123
      bp5 = tp123
      bp6 = .5*(tp122 + 2.*tp222 - tp123)
c          
      avgtp1 = 3./8.*bp1
      avgtp2 = -1./8.*(bp1 + 4.*bp2)
      avgtp3 = -1./8.*(bp1 + 4.*bp3)
      avgtp4 = -1./8.*(bp1 + 4.*bp4)
      avgtp5 = 1./8.*(bp1 + 4.*(bp3 + bp4 + 2.*bp5))
      avgtp6 = 1./8.*bp1 + bp2 + bp6
c          
      ab1 = 4.*avgtp1*avgtp6*avgt6+8.*(avgtp1+2.*
     x       avgtp2+avgtp3)*(avgtp2*avgt2+avgtp6*
     x       avgt2+avgtp2*avgt6)+(avgtp1+4.*avgtp2+
     x       3.*avgtp3)*(avgtp1*avgt4+4.*avgtp2*
     x       avgt4+avgtp4*avgt4+2.*avgtp6*avgt4+
     x       avgtp1*avgt5+4.*avgtp2*avgt5+3.*avgtp4*
     x       avgt5+2.*avgtp4*avgt6)+(avgtp1+4.*avgtp2+
     x       avgtp3+2.*avgtp6)*(2.*avgtp6*avgt1+4.*
     x       avgtp2*(avgt1+2.*avgt2+avgt3)+avgtp4*
     x       (avgt1+4.*avgt2+3.*avgt3)+avgtp1*(avgt1+
     x       4.*avgt2+avgt3+2.*avgt6))
      ab2 = 4.*(avgtp2*avgtp2*avgt2+2.*avgtp2*avgtp6*
     x       avgt2+avgtp6*avgtp6*avgt2+avgtp2*
     x       avgtp2*avgt6+2.*avgtp2*avgtp6*avgt6)
      ab3 = 4.*avgtp3*avgtp6*avgt6+(avgtp1+4.*avgtp2+
     x       3.*avgtp3)*(avgtp3*avgt4+avgtp5*avgt4+
     x       avgtp3*avgt5+3.*avgtp5*avgt5+2.*avgtp6*
     x       avgt5+2.*avgtp5*avgt6)+(avgtp1+4.*avgtp2+
     x       avgtp3+2.*avgtp6)*(2.*avgtp6*avgt3+
     x       avgtp5*(avgt1+4.*avgt2+3.*avgt3)+avgtp3*
     x       (avgt1+4.*avgt2+avgt3+2.*avgt6))
      ab4 = 4.*avgtp4*avgtp6*avgt6+8.*(avgtp4+
     x       avgtp5)*(avgtp2*avgt2+avgtp6*avgt2+
     x       avgtp2*avgt6)+(avgtp4+3.*avgtp5+2.*
     x       avgtp6)*(avgtp1*avgt4+4.*avgtp2*avgt4+
     x       avgtp4*avgt4+2.*avgtp6*avgt4+avgtp1*
     x       avgt5+4.*avgtp2*avgt5+3.*avgtp4*avgt5+
     x       2.*avgtp4*avgt6)+(avgtp4+avgtp5)*(2.*
     x       avgtp6*avgt1+4.*avgtp2*(avgt1+2.*avgt2+
     x       avgt3)+avgtp4*(avgt1+4.*avgt2+3.*avgt3)+
     x       avgtp1*(avgt1+4.*avgt2+avgt3+2.*avgt6))
      ab5 = 4.*avgtp5*avgtp6*avgt6+(avgtp4+3.*avgtp5+
     x       2.*avgtp6)*(avgtp3*avgt4+avgtp5*avgt4+
     x       avgtp3*avgt5+3.*avgtp5*avgt5+2.*avgtp6*
     x       avgt5+2.*avgtp5*avgt6)+(avgtp4+avgtp5)*
     x       (2.*avgtp6*avgt3+avgtp5*(avgt1+4.*
     x       avgt2+3.*avgt3)+avgtp3*(avgt1+4.*avgt2+
     x       avgt3+2.*avgt6))
      ab6 = 4.*avgtp6*avgtp6*avgt6
c
c --- Call elast
c
      up1 = -ph1*(si1*(pp11+4.*pp21+2.*pp61)+4.*si2*(pp11+2.*pp21+
     x           pp31)+si4*(pp11+4.*pp21+3.*pp31)+2.*si6*pp11)
     x      -ph2*(si1*(pp12+4.*pp22+2.*pp62)+4.*si2*(pp12+2.*pp22+
     x           pp32)+si4*(pp12+4.*pp22+3.*pp32)+2.*si6*pp12)
      up2 = -ph1*(2.*si2*(pp21+pp61)+2.*si6*pp21)
     x      -ph2*(2.*si2*(pp22+pp62)+2.*si6*pp22)
      up3 = -ph1*(si3*(pp11+4.*pp21+pp31+2.*pp61)+si5*(pp11+4.*pp21+
     x           3.*pp31)+2.*si6*pp31)
     x      -ph2*(si3*(pp12+4.*pp22+pp32+2.*pp62)+si5*(pp12+4.*pp22+
     x           3.*pp32)+2.*si6*pp32)
      up4 = -ph1*(si1*(pp41+pp51)+4.*si2*(pp41+pp51)+si4*(pp41+3.*
     x           pp51+2.*pp61)+2.*si6*pp41)
     x      -ph2*(si1*(pp42+pp52)+4.*si2*(pp42+pp52)+si4*(pp42+3.*
     x           pp52+2.*pp62)+2.*si6*pp42)
      up5 = -ph1*(si3*(pp41+pp51)+si5*(pp41+3.*pp51+2.*pp61)+
     x           2.*si6*pp51)
     x      -ph2*(si3*(pp42+pp52)+si5*(pp42+3.*pp52+2.*pp62)+
     x           2.*si6*pp52)
      up6 = .5-ph1*2.*si6*pp61-ph2*2.*si6*pp62
c
      db1 = up1 + 4.*up2 + up3 + 2.*up6
      db2 = up1 + 4.*up2 + 3.*up3
      db3 = up4 + up5
      db4 = up4 + 3.*up5 + 2.*up6
c
      dbin1 = db4/(db1*db4 - db2*db3)
      dbin2 = -db2/(db1*db4 - db2*db3)
      dbin3 = -db3/(db1*db4 - db2*db3)
      dbin4 = db1/(db1*db4 - db2*db3)
c    
      upp6 = 1./(4.*up6)
      upp2 = -up2/(4.*up6*(up2 + up6))
      upp1 = dbin1*(-2.*upp6*up1 - 4.*upp2*(up1 + 2.*up2 + up3)) +
     x       dbin2*(-2.*upp6*up4 - 4.*upp2*(up4 + up5))
      upp4 = dbin3*(-2.*upp6*up1 - 4.*upp2*(up1 + 2.*up2 + up3)) + 
     x       dbin4*(-2.*upp6*up4 - 4.*upp2*(up4 + up5))
      upp3 = dbin1*(-2.*upp6*up3) + dbin2*(-2.*upp6*up5)
      upp5 = dbin3*(-2.*upp6*up3) + dbin4*(-2.*upp6*up5)
c
      all1 = ph1*(pp11*(upp1 + 4.*upp2 + 2.*upp6) + 4.*pp21*(upp1 + 
     x       2.*upp2 + upp3) + pp41*(upp1 + 4.*upp2 + 3.*upp3) + 
     x       2.*pp61*upp1) + ph2*(pp12*(upp1 + 4.*upp2 + 2.*upp6) + 
     x       4.*pp22*(upp1 + 2.*upp2 + upp3) + pp42*(upp1 + 4.*upp2 + 
     x       3.*upp3) + 2.*pp62*upp1)
      all2 = ph1*(2.*pp21*(upp2 + upp6) + 2.*pp61*upp2) +
     x       ph2*(2.*pp22*(upp2 + upp6) + 2.*pp62*upp2)
      all3 = ph1*(pp31*(upp1 + 4.*upp2 + upp3 + 2.*upp6) + 
     x       pp51*(upp1 + 4.*upp2 + 3.*upp3) + 2.*pp61*upp3) +
     x       ph2*(pp32*(upp1 + 4.*upp2 + upp3 + 2.*upp6) + 
     x       pp52*(upp1 + 4.*upp2 + 3.*upp3) + 2.*pp62*upp3)
      all4 = ph1*(pp11*(upp4 + upp5) + 4.*pp21*(upp4 + upp5) + 
     x       pp41*(upp4 + 3.*upp5 + 2.*upp6) + 2.*pp61*upp4) +
     x       ph2*(pp12*(upp4 + upp5) + 4.*pp22*(upp4 + upp5) + 
     x       pp42*(upp4 + 3.*upp5 + 2.*upp6) + 2.*pp62*upp4)
      all5 = ph1*(pp31*(upp4 + upp5) + pp51*(upp4 + 3.*upp5 + 
     x       2.*upp6) + 2.*pp61*upp5) + ph2*(pp32*(upp4 + upp5) + 
     x       pp52*(upp4 + 3.*upp5 + 2.*upp6) + 2.*pp62*upp5)
      all6 = .5 + ph1*2.*pp61*upp6 + ph2*2.*pp62*upp6
c
      aji1 = 2.*amu1*all1
      aji2 = 2.*amu1*all2
      aji3 = alam1*(all1 + 4.*all2 + 3.*all3) + 2.*amu1*all3
      aji4 = 2.*amu1*all4
      aji5 = alam1*(all4 + 3.*all5 + 2.*all6) + 2.*amu1*all5
      aji6 = 2.*amu1*all6
c
      abar1 = 3./8.*aji1
      abar2 = -1./8.*(aji1 + 4.*aji2)
      abar3 = -1./8.*(aji1 + 4.*aji3)
      abar4 = -1./8.*(aji1 + 4.*aji4)
      abar5 = 1./8.*(aji1 + 4.*(aji3 + aji4 + 2.*aji5))
      abar6 = 1./8.*aji1 + aji2 + aji6 
c
c --- For 3D random distribution
c      td1=1./15.*(aji1+5.*(aji3+aji4+3.*aji5))
c      td2=1./15.*(aji1+10.*aji2+15.*aji6)     
c 
c      deps(3)=-enu/(1.-enu)*(deps(1)+deps(2))
c
c --- For 2D random distribution
c
      eeps = -1./3.*(eps1 + eps2 + eps3)
      zz1 = -3.*aji5
      zz2 = 2.*aji6
      zz3 = .5*zz2
c     
      sig(1) = zz1*eeps + zz2*eps1
      sig(2) = zz1*eeps + zz2*eps2
      sig(3) = 0.
      sig(4) = zz3*eps12
      sig(5) = sig(5) + capa*zz3*deps(5)
      sig(6) = sig(6) + capa*zz3*deps(6)
c      sig(1)=(aji5+2.*aji6)*eps1+aji5*eps2
c      sig(2)=aji5*eps1+(aji5+2.*aji6)*eps2
c      sig(1)=(aji5+2.*aji6)*eps1+aji5*eps2
c    x       +(aji3+aji5)*eps3
c      sig(2)=aji5*eps1+(aji5+2.*aji6)*eps2
c    x       +(aji3+aji5)*eps3
c      sig(3)=0.
c      sig(4)=aji6*eps12
c      sig(5)=sig(5)+ee/(2.*(1.+enu))*capa*deps(5)
c      sig(6)=sig(6)+ee/(2.*(1.+enu))*capa*deps(5)    
c      write(*,*)zz1,zz2
c
c --- Call internal
c  
      hh1 = -ph1*aj11 - ph2*aj12
      hh2 = -ph1*aj21 - ph2*aj22
      hh3 = -ph1*aj31 - ph2*aj32
      hh4 = -ph1*aj41 - ph2*aj42
      hh5 = -ph1*aj51 - ph2*aj52
      hh6 = .5 - ph1*aj61 - ph2*aj62
c
      dhh1 = hh1 + 4.*hh2 + hh3 + 2.*hh6
      dhh2 = hh1 + 4.*hh2 + 3.*hh3
      dhh3 = hh4 + hh5
      dhh4 = hh4 + 3.*hh5 + 2.*hh6
c 
      dhhin1 = dhh4/(dhh1*dhh4 - dhh2*dhh3)
      dhhin2 = -dhh2/(dhh1*dhh4 - dhh2*dhh3)
      dhhin3 = -dhh3/(dhh1*dhh4 - dhh2*dhh3)
      dhhin4 = dhh1/(dhh1*dhh4 - dhh2*dhh3)
c    
      aij6 = 1./(4.*hh6)
      aij2 = -hh2/(4.*hh6*(hh2 + hh6))
      aij1 = dhhin1*(2.*aij6*hh1 - 4.*aij2*(hh1 + 2.*hh2 + hh3)) +
     x       dhhin2*(-2.*aij6*hh4 - 4.*aij2*(hh4 + hh5))
      aij4 = dhhin3*(2.*aij6*hh1 - 4.*aij2*(hh1 + 2.*hh2 + hh3)) + 
     x       dhhin4*(-2.*aij6*hh4 - 4.*aij2*(hh4 + hh5))
      aij3 = dhhin1*(-2.*aij6*hh3) + dhhin2*(-2.*aij6*hh5)
      aij5 = dhhin3*(-2.*aij6*hh3) + dhhin4*(-2.*aij6*hh5)
c
      bk1 = ajj11*(aij1 + 4.*aij2 + aij3 + 2.*aij6) + 4.*ajj21*
     x      (aij1 + 2.*aij2 + aij3) + ajj41*(aij1 + 4.*aij2 + 
     x      3.*aij3) + 2.*ajj61*aij1
      bk2 = 2.*ajj21*(aij2 + aij6) + 2.*ajj61*aij2
      bk3 = ajj31*(aij1 + 4.*aij2 + aij3 + 2.*aij6) + ajj51*
     x      (aij1 + 4.*aij2 + 3.*aij3) + 2.*ajj61*aij3
      bk4 = ajj11*(aij4 + aij5) + 4.*ajj21*(aij4 + aij5) + 
     x      ajj41*(aij4 + 3.*aij5 + 2.*aij6) + 2.*ajj61*aij4
      bk5 = ajj31*(aij4 + aij5) + ajj51*(aij4 + 3.*aij5 + 
     x      2.*aij6) + 2.*ajj61*aij5
      bk6 = 2.*ajj61*aij6
c
      rr1 = 2.*amu2*bk1
      rr2 = 2.*amu2*bk2
      rr3 = alam2*(bk1 + 4.*bk2 + 3.*bk3) + 2.*amu2*bk3
      rr4 = 2.*amu2*bk4
      rr5 = alam2*(bk4 + 3.*bk5 + 2.*bk6) + 2.*amu2*bk5
      rr6 = 2.*amu2*bk6
c
      upi1 = 3./8.*rr1
      upi2 = -1./8.*(rr1 + 4.*rr2)
      upi3 = -1./8.*(rr1 + 4.*rr3)
      upi4 = -1./8.*(rr1 + 4.*rr4)
      upi5 = 1./8.*(rr1 + 4.*(rr3 + rr4 + 2.*rr5))
      upi6 = 1./8.*rr1 + rr2 + rr6
c
      psig1 = (upi5 + 2.*upi6)*eps1 + upi5*eps2
      psig2 = upi5*eps1 + (upi5 + 2.*upi6)*eps2
      psig12 = upi6*eps12
c
      epsig = sqrt(((psig1 - psig2)**2 + psig2**2 + psig1**2 +
     x        6.*psig12**2)/2.)
      proba = 1. - exp(-(epsig/s0)**am)
c
c --- Damage variable
c
      pph2 = ph*proba
      if (pph2 .ge. oph2) then
        ph2 = pph2
      else
        ph2 = oph2
      endif
c     
      if (ph2 .le. oph2) then 
        ph2 = oph2
      endif
      oph2 = ph2
      ph1 = ph-ph2
c
      hisv(1) = ph1
      hisv(2) = ph2
      hisv(8) = oph2  
      return
      end 
