      SUBROUTINE STRESS(P,XN,U,TAU)
C
C  ROUTINE FOR THE COMPUTATION OF NORMAL STRESSES
C
      IMPLICIT COMPLEX (P,U,T,E)
      REAL N1,N2,N3
      DIMENSION XN(3),P(3),U(3),TAU(3)
      COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33,
     1               A34,A35,A36,A44,A45,A46,A55,A56,A66,
     1               DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23,
     1               DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44,
     1               DXA45,DXA46,DXA55,DXA56,DXA66,
     1               DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23,
     1               DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44,
     1               DYA45,DYA46,DYA55,DYA56,DYA66,
     1               DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23,
     1               DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44,
     1               DZA45,DZA46,DZA55,DZA56,DZA66,
     1               A2546,A1266,A1355,A1456,A3645,A2344
      COMMON /AUXI/  IANI(20),INTR,INT1,IOUT,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori
      common/dens/rho(20)
C
      P1=P(1)
      P2=P(2)
      P3=P(3)
      U1=U(1)
      U2=U(2)
      U3=U(3)
      N1=XN(1)
      N2=XN(2)
      N3=XN(3)
      IF(IANI(LAY).EQ.0) GO TO 1000
C
C  ANISOTROPIC CASE
C
      RO=1.7+0.2*SQRT(A11)
      IF(IRHO.NE.0) RO=rho(lay)
C
C  DISPLACEMENT TENSOR
C
      E11=U1*P1
      E12=0.5*(U1*P2+U2*P1)
      E13=0.5*(U1*P3+U3*P1)
      E22=U2*P2
      E23=0.5*(U2*P3+U3*P2)
      E33=U3*P3
C
C  STRESS TENSOR
C
      P11=A11*E11+A12*E22+A13*E33+2.*A14*E23+2.*A15*E13+2.*A16*E12
      P22=A12*E11+A22*E22+A23*E33+2.*A24*E23+2.*A25*E13+2.*A26*E12
      P33=A13*E11+A23*E22+A33*E33+2.*A34*E23+2.*A35*E13+2.*A36*E12
      P23=A14*E11+A24*E22+A34*E33+2.*A44*E23+2.*A45*E13+2.*A46*E12
      P13=A15*E11+A25*E22+A35*E33+2.*A45*E23+2.*A55*E13+2.*A56*E12
      P12=A16*E11+A26*E22+A36*E33+2.*A46*E23+2.*A56*E13+2.*A66*E12
C
C  NORMAL STRESS
C
      TAU(1)=P11*N1+P12*N2+P13*N3
      TAU(2)=P12*N1+P22*N2+P23*N3
      TAU(3)=P13*N1+P23*N2+P33*N3
      TAU(1)=TAU(1)*RO
      TAU(2)=TAU(2)*RO
      TAU(3)=TAU(3)*RO
      RETURN
C
C  ISOTROPIC CASE
C
 1000 A12=A11-2.*A44
      RO=1.7+0.2*SQRT(A11)
      IF(IRHO.NE.0) RO=rho(lay)
      IF(ICOEF.EQ.0) GOTO 100
      WRITE(LOUT,'(A,4F12.7)') 'STRESS: A11,A44,A12,RO',A11,A44,A12,RO
      WRITE(LOUT,'(A,/,3(2F12.7,/))') 'STRESS: PSTR',P
      WRITE(LOUT,'(A,/,3(2F12.7,/))') 'STRESS: USTR',U
  100 TETA=U1*P1+U2*P2+U3*P3
      IF(ICOEF.GT.0)
     1WRITE(LOUT,'(A,/,5F12.7)') 'STRESS: UN,TETA',XN,TETA
      P11=A12*TETA+2.*A44*P1*U1
      P12=A44*(P1*U2+P2*U1)
      P13=A44*(P1*U3+P3*U1)
      P22=A12*TETA+2.*A44*P2*U2
      P23=A44*(P2*U3+P3*U2)
      P33=A12*TETA+2.*A44*P3*U3
      IF(ICOEF.GT.0)
     1WRITE(LOUT,'(A,/,6(2F12.7,/))')'STRESS: P11,P12,P13,P22,P23,P33'
     1,P11,P12,P13,P22,P23,P33
      TAU(1)=P11*N1+P12*N2+P13*N3
      TAU(2)=P12*N1+P22*N2+P23*N3
      TAU(3)=P13*N1+P23*N2+P33*N3
      TAU(1)=TAU(1)*RO
      TAU(2)=TAU(2)*RO
      TAU(3)=TAU(3)*RO
      RETURN
      END
