C
C Subroutine file 'wanpfa.for' with subroutines used in 'crtpfv.for' and
C 'wan.for'
C
C Version: 7.00
C Date: 2013, June 12
C
C Coded by: Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.cz/staff/bulant.htm
C
C=======================================================================
C
      SUBROUTINE WACHRI(P1,P2,P3,B11,B12,B22,B13,B23,B33,
     *                  B14,B24,B34,B44,B15,B25,B35,B45,B55,
     *                  B16,B26,B36,B46,B56,B66,
     *                  G1,G2,G3,EE,DER)
C
C----------------------------------------------------------------------
C Subroutine to compute the Christoffel matrix, its eigenvalues
C and eigenvectors.
      REAL P1,P2,P3
      REAL B11,B12,B22,B13,B23,B33,B14,B24,B34,B44
      REAL B15,B25,B35,B45,B55,B16,B26,B36,B46,B56,B66
      REAL G1,G2,G3,EE(9),DER(9)
C
C Input:
C     P1,P2,P3... Slowness vector.
C     Bii ... Values of real parts of 21 reduced
C             (divided by the density) elastic parameters.
C Output:
C     G1,G2,G3 ... Eigenvalues of the Christoffel matrix.
C     EE  ... Eigenvectors of the Christoffel matrix.
C     DER ... Derivatives dx/dt=dH/dp=Aijkl Ej Ek pl stored columnwise
C             for qP, qS1 and qS2 waves.
C
C-----------------------------------------------------------------------
C
      REAL A11,A12,A13,A14,A21,A22,A23,A24,A31,A32,A33,A34,A44
      REAL A15,A25,A35,A45,A55,A16,A26,A36,A46,A56,A66
      REAL A1111,A2111,A3111,A1211,A2211,A3211,A1311,A2311,A3311
      REAL A1121,A2121,A3121,A1221,A2221,A3221,A1321,A2321,A3321
      REAL A1131,A2131,A3131,A1231,A2231,A3231,A1331,A2331,A3331
      REAL A1112,A2112,A3112,A1212,A2212,A3212,A1312,A2312,A3312
      REAL A1122,A2122,A3122,A1222,A2222,A3222,A1322,A2322,A3322
      REAL A1132,A2132,A3132,A1232,A2232,A3232,A1332,A2332,A3332
      REAL A1113,A2113,A3113,A1213,A2213,A3213,A1313,A2313,A3313
      REAL A1123,A2123,A3123,A1223,A2223,A3223,A1323,A2323,A3323
      REAL A1133,A2133,A3133,A1233,A2233,A3233,A1333,A2333,A3333
      EQUIVALENCE (A11,A1111)
      EQUIVALENCE (A22,A2222)
      EQUIVALENCE (A33,A3333)
      EQUIVALENCE (A16,A1112,A1121,A1211,A2111)
      EQUIVALENCE (A26,A2221,A2212,A2122,A1222)
      EQUIVALENCE (A15,A1113,A1131,A1311,A3111)
      EQUIVALENCE (A35,A3331,A3313,A3133,A1333)
      EQUIVALENCE (A24,A2223,A2232,A2322,A3222)
      EQUIVALENCE (A34,A3332,A3323,A3233,A2333)
      EQUIVALENCE (A23,A2233,A3322)
      EQUIVALENCE (A13,A1133,A3311)
      EQUIVALENCE (A12,A1122,A2211)
      EQUIVALENCE (A44,            A2323,A3232,A2332,A3223)
      EQUIVALENCE (A55,            A1313,A3131,A1331,A3113)
      EQUIVALENCE (A66,            A1212,A2121,A1221,A2112)
      EQUIVALENCE (A14,A1123,A1132,A2311,A3211)
      EQUIVALENCE (A25,A2213,A2231,A1322,A3122)
      EQUIVALENCE (A36,A3312,A3321,A1233,A2133)
      EQUIVALENCE (A56,A1321,A3112,A2113,A1231,A1213,A2131,A1312,A3121)
      EQUIVALENCE (A46,A2312,A3221,A1223,A2132,A2123,A1232,A2321,A3212)
      EQUIVALENCE (A45,A3213,A2331,A1332,A3123,A3132,A1323,A3231,A2313)
      REAL GAMMA(6),E11,E21,E31,E12,E22,E32,E13,E23,E33
C     EQUIVALENCE (GAMMA(1),G1),(GAMMA(3),G2),(GAMMA(6),G3)
C     EQUIVALENCE (EE(1),E11),(EE(4),E12),(EE(7),E13)
C     EQUIVALENCE (EE(2),E21),(EE(5),E22),(EE(8),E23)
C     EQUIVALENCE (EE(3),E31),(EE(6),E32),(EE(9),E33)
      REAL A111,A112,A121,A122,A113,A123,A131,A132,A133
      REAL A211,A212,A221,A222,A213,A223,A231,A232,A233
      REAL A311,A312,A322,A313,A321,A323,A331,A332,A333
      REAL AUX
      INTEGER KQICHM
      DATA    KQICHM/-1/
C
C     GAMMA,G1,G2,G3...Christoffel matrix, later its eigenvalues.
C        (E11,E12,E13)
C     EE=(E21,E22,E23)... Eigenvectors of the christoffel matrix.
C        (E31,E32,E33)
C     A111,A211,A311,A112,A212,A312,A122,A222,A322,A113,A213,A313,A123,
C     A223,A323,A133,A233,A333... A(I,J,K,L)*P(L) summed over L.
C     A11,A21,A31,A12,A22,A32,A13,A23,A33 ... Aijk*Ek
C
C.......................................................................
C
      A11=B11
      A22=B22
      A33=B33
      A16=B16
      A26=B26
      A15=B15
      A35=B35
      A24=B24
      A34=B34
      A23=B23
      A13=B13
      A12=B12
      A44=B44
      A55=B55
      A66=B66
      A14=B14
      A25=B25
      A36=B36
      A56=B56
      A46=B46
      A45=B45
C     Christoffel matrix:
      A111=A1111*P1+A1112*P2+A1113*P3
      A112=A1121*P1+A1122*P2+A1123*P3
      A121=A1211*P1+A1212*P2+A1213*P3
      A122=A1221*P1+A1222*P2+A1223*P3
      A113=A1131*P1+A1132*P2+A1133*P3
      A123=A1231*P1+A1232*P2+A1233*P3
      A131=A1311*P1+A1312*P2+A1313*P3
      A132=A1321*P1+A1322*P2+A1323*P3
      A133=A1331*P1+A1332*P2+A1333*P3
      A211=A2111*P1+A2112*P2+A2113*P3
      A212=A2121*P1+A2122*P2+A2123*P3
      A221=A2211*P1+A2212*P2+A2213*P3
      A222=A2221*P1+A2222*P2+A2223*P3
      A213=A2131*P1+A2132*P2+A2133*P3
      A223=A2231*P1+A2232*P2+A2233*P3
      A231=A2311*P1+A2312*P2+A2313*P3
      A232=A2321*P1+A2322*P2+A2323*P3
      A233=A2331*P1+A2332*P2+A2333*P3
      A311=A3111*P1+A3112*P2+A3113*P3
      A312=A3121*P1+A3122*P2+A3123*P3
      A322=A3221*P1+A3222*P2+A3223*P3
      A313=A3131*P1+A3132*P2+A3133*P3
      A321=A3211*P1+A3212*P2+A3213*P3
      A323=A3231*P1+A3232*P2+A3233*P3
      A331=A3311*P1+A3312*P2+A3313*P3
      A332=A3321*P1+A3322*P2+A3323*P3
      A333=A3331*P1+A3332*P2+A3333*P3
      GAMMA(1)=P1*A111+P2*A211+P3*A311
      GAMMA(2)=P1*A112+P2*A212+P3*A312
      GAMMA(3)=P1*A122+P2*A222+P3*A322
      GAMMA(4)=P1*A113+P2*A213+P3*A313
      GAMMA(5)=P1*A123+P2*A223+P3*A323
      GAMMA(6)=P1*A133+P2*A233+P3*A333
C
C     Quasi-isotropic modification of the Christoffel matrix:
      IF(KQICHM.LT.0) THEN
        CALL RSEP3I('QICHM',KQICHM,0)
      END IF
      IF(KQICHM.GT.0) THEN
        AUX=SQRT(P1*P1+P2*P2+P3*P3)
        E11=P1/AUX
        E21=P2/AUX
        E31=P3/AUX
        E12=GAMMA(1)*E11+GAMMA(2)*E21+GAMMA(4)*E31
        E22=GAMMA(2)*E11+GAMMA(3)*E21+GAMMA(5)*E31
        E32=GAMMA(4)*E11+GAMMA(5)*E21+GAMMA(6)*E31
        AUX=2.*(E11*E12+E21*E22+E31*E32)
        GAMMA(1)=GAMMA(1)-E11*E12-E12*E11+AUX*E11*E11
        GAMMA(2)=GAMMA(2)-E11*E22-E12*E21+AUX*E11*E21
        GAMMA(3)=GAMMA(3)-E21*E22-E22*E21+AUX*E21*E21
        GAMMA(4)=GAMMA(4)-E11*E32-E12*E31+AUX*E11*E31
        GAMMA(5)=GAMMA(5)-E21*E32-E22*E31+AUX*E21*E31
        GAMMA(6)=GAMMA(6)-E31*E32-E32*E31+AUX*E31*E31
      END IF
C
C     Selecting eigenvalue and eigenvector of the Christoffel matrix:
      CALL EIGEN(GAMMA,EE,3,0)
      G1=GAMMA(1)
      G2=GAMMA(3)
      G3=GAMMA(6)
      E11=EE(1)
      E21=EE(2)
      E31=EE(3)
      E12=EE(4)
      E22=EE(5)
      E32=EE(6)
      E13=EE(7)
      E23=EE(8)
      E33=EE(9)
      IF (G3.LT.0.) THEN
C       WAN-11
        CALL ERROR('WAN-11: Negative eigenvalue of Christoffel matrix')
C       This error should not appear.
      END IF
      AUX=E11*E11+E21*E21+E31*E31
      IF (ABS(AUX-1.).GT.0.00001) THEN
C       WAN-12
        CALL ERROR('WAN-12: Eigenvector is not normalized')
C       This error should not appear.
      ENDIF
      AUX=E12*E12+E22*E22+E32*E32
      IF (ABS(AUX-1.).GT.0.00001) THEN
C       WAN-13
        CALL ERROR('WAN-13: Eigenvector is not normalized')
C       This error should not appear.
      ENDIF
      AUX=E13*E13+E23*E23+E33*E33
      IF (ABS(AUX-1.).GT.0.00001) THEN
C       WAN-14
        CALL ERROR('WAN-14: Eigenvector is not normalized')
C       This error should not appear.
      ENDIF
C
C     Computation of derivatives dx/dt:
      A11=  A111*E11+A112*E21+A113*E31
      A21=  A211*E11+A212*E21+A213*E31
      A31=  A311*E11+A312*E21+A313*E31
      A12=  A121*E11+A122*E21+A123*E31
      A22=  A221*E11+A222*E21+A223*E31
      A32=  A321*E11+A322*E21+A323*E31
      A13=  A131*E11+A132*E21+A133*E31
      A23=  A231*E11+A232*E21+A233*E31
      A33=  A331*E11+A332*E21+A333*E31
      DER(1)=A11*E11+ A12*E21+ A13*E31
      DER(2)=A21*E11+ A22*E21+ A23*E31
      DER(3)=A31*E11+ A32*E21+ A33*E31
      A11=  A111*E12+A112*E22+A113*E32
      A21=  A211*E12+A212*E22+A213*E32
      A31=  A311*E12+A312*E22+A313*E32
      A12=  A121*E12+A122*E22+A123*E32
      A22=  A221*E12+A222*E22+A223*E32
      A32=  A321*E12+A322*E22+A323*E32
      A13=  A131*E12+A132*E22+A133*E32
      A23=  A231*E12+A232*E22+A233*E32
      A33=  A331*E12+A332*E22+A333*E32
      DER(4)=A11*E12+ A12*E22+ A13*E32
      DER(5)=A21*E12+ A22*E22+ A23*E32
      DER(6)=A31*E12+ A32*E22+ A33*E32
      A11=  A111*E13+A112*E23+A113*E33
      A21=  A211*E13+A212*E23+A213*E33
      A31=  A311*E13+A312*E23+A313*E33
      A12=  A121*E13+A122*E23+A123*E33
      A22=  A221*E13+A222*E23+A223*E33
      A32=  A321*E13+A322*E23+A323*E33
      A13=  A131*E13+A132*E23+A133*E33
      A23=  A231*E13+A232*E23+A233*E33
      A33=  A331*E13+A332*E23+A333*E33
      DER(7)=A11*E13+ A12*E23+ A13*E33
      DER(8)=A21*E13+ A22*E23+ A23*E33
      DER(9)=A31*E13+ A32*E23+ A33*E33
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE WAMAT(A,B,C,D)
C
C----------------------------------------------------------------------
C Subroutine to compute the product of two 2x2 complex matrices.
C The second matrix (C+iD) is destroyed in the computation.
      REAL A(2,2),B(2,2),C(2,2),D(2,2)
C Input:
C     A,B,C,D ... Real and imaginary parts of the two matrices.
C Output:
C     C,D ... Real and imaginary parts of resulting matrix.
C
C.......................................................................
C     Auxiliary storage locations:
      REAL CR11,CR21,CR12,CR22,CI11,CI21,CI12,CI22
C.......................................................................
      CR11=A(1,1)*C(1,1)-B(1,1)*D(1,1)+A(1,2)*C(2,1)-B(1,2)*D(2,1)
      CR21=A(2,1)*C(1,1)-B(2,1)*D(1,1)+A(2,2)*C(2,1)-B(2,2)*D(2,1)
      CR12=A(1,1)*C(1,2)-B(1,1)*D(1,2)+A(1,2)*C(2,2)-B(1,2)*D(2,2)
      CR22=A(2,1)*C(1,2)-B(2,1)*D(1,2)+A(2,2)*C(2,2)-B(2,2)*D(2,2)
C
      CI11=A(1,1)*D(1,1)+B(1,1)*C(1,1)+A(1,2)*D(2,1)+B(1,2)*C(2,1)
      CI21=A(2,1)*D(1,1)+B(2,1)*C(1,1)+A(2,2)*D(2,1)+B(2,2)*C(2,1)
      CI12=A(1,1)*D(1,2)+B(1,1)*C(1,2)+A(1,2)*D(2,2)+B(1,2)*C(2,2)
      CI22=A(2,1)*D(1,2)+B(2,1)*C(1,2)+A(2,2)*D(2,2)+B(2,2)*C(2,2)
C
      C(1,1)=CR11
      C(2,1)=CR21
      C(1,2)=CR12
      C(2,2)=CR22
      D(1,1)=CI11
      D(2,1)=CI21
      D(1,2)=CI12
      D(2,2)=CI22
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE WAVECR(G,H,U,V,RU,RV)
C
C----------------------------------------------------------------------
C Subroutine to reorganize the vectors U and V in such way,
C that the pair U,V is equal to pair G,H rotated with angle smaller
C than 45 degree (clockwise or anticlockwise).
C The real numbers RU (associated with U) and RV (associated with V),
C are reorganized in the same way.
C
      REAL G(3),H(3),U(3),V(3),RU,RV
C Input:
C     G,H ... A pair of orthonormal vectors.
C     U,V ... A pair of orthonormal vectors.
C     RU,RV ... Real numbers associated to U and V.
C Output:
C     U,V ... Selection from input vectors U,V,-U,-V in such way, that
C             the output pair U,V is equal to pair G,H rotated
C             with angle smaller than 45 degree.
C     RU,RV ... Real numbers associated to U and V.
C
C.......................................................................
C
C     Auxiliary storage locations:
      REAL E1,E2,E3,AUX
C
C-----------------------------------------------------------------------
        AUX=G(1)*U(1)+G(2)*U(2)+G(3)*U(3)
        IF(ABS(AUX).LT.0.707107) THEN
          E1=U(1)
          E2=U(2)
          E3=U(3)
          U(1)=V(1)
          U(2)=V(2)
          U(3)=V(3)
          V(1)=E1
          V(2)=E2
          V(3)=E3
          AUX =RU
          RU  =RV
          RV  =AUX
        END IF
        AUX=G(1)*U(1)+G(2)*U(2)+G(3)*U(3)
        IF(AUX.LT.0.0) THEN
          U(1)=-U(1)
          U(2)=-U(2)
          U(3)=-U(3)
        END IF
        AUX=H(1)*V(1)+H(2)*V(2)+H(3)*V(3)
        IF(AUX.LT.0.0) THEN
          V(1)=-V(1)
          V(2)=-V(2)
          V(3)=-V(3)
        END IF
      RETURN
      END
C
C
C=======================================================================
C
      REAL FUNCTION WAANGL(G,H,U,V)
C
C----------------------------------------------------------------------
C Subroutine to compute the rotation angle between the pair of
C orthonormal vectors G and H and the pair of orthonormal vectors U and
C V according to equation 18 of Bulant & Klimes, 2002, Pageoph.
C
      REAL G(3),H(3),U(3),V(3)
C Input:
C     G,H ... A pair of orthonormal vectors.
C     U,V ... A pair of orthonormal vectors.
C Output:
C     WAANGL  ... The rotation angle.
C
C.......................................................................
C Subroutines and external functions required:
      EXTERNAL WASUM
      REAL WASUM
C     WASUM ... This file.
C-----------------------------------------------------------------------
      WAANGL=ATAN2((WASUM(U,H)-WASUM(V,G)),(WASUM(U,G)+WASUM(V,H)))
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE WAPROJ(P1,P2,P3,E1,E2,E3,G1,G2,G3,H1,H2,H3,
     *                  GE1,GE2,HE1,HE2)
C
C----------------------------------------------------------------------
C Subroutine to project vectors G and H to the plane defined by vector E
C and vector PxE.
C
      REAL P1,P2,P3,E1,E2,E3,G1,G2,G3,H1,H2,H3,GE1,GE2,HE1,HE2
C Input:
C     P1,P2,P3,E1,E2,E3 ... Vectors defining the plane.
C     G1,G2,G3,H1,H2,H3 ... Vectors to be projected.
C Output:
C     GE1,GE2,HE1,HE2 ... Projected vectors.
C
C.......................................................................
C
C     Auxiliary storage locations:
      REAL F1,F2,F3,AUX
C
C-----------------------------------------------------------------------
C     Second vector defining the plane:
      F1=P2*E3-E2*P3
      F2=P3*E1-E3*P1
      F3=P1*E2-E1*P2
      AUX=SQRT(F1*F1+F2*F2+F3*F3)
      F1=F1/AUX
      F2=F2/AUX
      F3=F3/AUX
C     Projecting vectors G and H to the plane defined by E and F:
      GE1=E1*G1+E2*G2+E3*G3
      GE2=F1*G1+F2*G2+F3*G3
      HE1=E1*H1+E2*H2+E3*H3
      HE2=F1*H1+F2*H2+F3*H3
C     Normalizing vectors GE and HE:
      AUX=SQRT(GE1*GE1+GE2*GE2)
      GE1=GE1/AUX
      GE2=GE2/AUX
      AUX=SQRT(HE1*HE1+HE2*HE2)
      HE1=HE1/AUX
      HE2=HE2/AUX
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE WAIJKL(A,B,I)
C
C----------------------------------------------------------------------
      REAL A(10,21),B(3,3,3,3)
      INTEGER I
      B(1,1,1,1)=A(I,1)
      B(1,1,2,2)=A(I,2)
      B(2,2,1,1)=A(I,2)
      B(2,2,2,2)=A(I,3)
      B(1,1,3,3)=A(I,4)
      B(3,3,1,1)=A(I,4)
      B(2,2,3,3)=A(I,5)
      B(3,3,2,2)=A(I,5)
      B(3,3,3,3)=A(I,6)
      B(1,1,2,3)=A(I,7)
      B(1,1,3,2)=A(I,7)
      B(2,3,1,1)=A(I,7)
      B(3,2,1,1)=A(I,7)
      B(2,2,2,3)=A(I,8)
      B(2,2,3,2)=A(I,8)
      B(2,3,2,2)=A(I,8)
      B(3,2,2,2)=A(I,8)
      B(3,3,3,2)=A(I,9)
      B(3,3,2,3)=A(I,9)
      B(3,2,3,3)=A(I,9)
      B(2,3,3,3)=A(I,9)
      B(2,3,2,3)=A(I,10)
      B(3,2,3,2)=A(I,10)
      B(2,3,3,2)=A(I,10)
      B(3,2,2,3)=A(I,10)
      B(1,1,1,3)=A(I,11)
      B(1,1,3,1)=A(I,11)
      B(1,3,1,1)=A(I,11)
      B(3,1,1,1)=A(I,11)
      B(2,2,1,3)=A(I,12)
      B(2,2,3,1)=A(I,12)
      B(1,3,2,2)=A(I,12)
      B(3,1,2,2)=A(I,12)
      B(3,3,3,1)=A(I,13)
      B(3,3,1,3)=A(I,13)
      B(3,1,3,3)=A(I,13)
      B(1,3,3,3)=A(I,13)
      B(3,2,1,3)=A(I,14)
      B(2,3,3,1)=A(I,14)
      B(1,3,3,2)=A(I,14)
      B(3,1,2,3)=A(I,14)
      B(3,1,3,2)=A(I,14)
      B(1,3,2,3)=A(I,14)
      B(3,2,3,1)=A(I,14)
      B(2,3,1,3)=A(I,14)
      B(1,3,1,3)=A(I,15)
      B(3,1,3,1)=A(I,15)
      B(1,3,3,1)=A(I,15)
      B(3,1,1,3)=A(I,15)
      B(1,1,1,2)=A(I,16)
      B(1,1,2,1)=A(I,16)
      B(1,2,1,1)=A(I,16)
      B(2,1,1,1)=A(I,16)
      B(2,2,2,1)=A(I,17)
      B(2,2,1,2)=A(I,17)
      B(2,1,2,2)=A(I,17)
      B(1,2,2,2)=A(I,17)
      B(3,3,1,2)=A(I,18)
      B(3,3,2,1)=A(I,18)
      B(1,2,3,3)=A(I,18)
      B(2,1,3,3)=A(I,18)
      B(2,3,1,2)=A(I,19)
      B(3,2,2,1)=A(I,19)
      B(1,2,2,3)=A(I,19)
      B(2,1,3,2)=A(I,19)
      B(2,1,2,3)=A(I,19)
      B(1,2,3,2)=A(I,19)
      B(2,3,2,1)=A(I,19)
      B(3,2,1,2)=A(I,19)
      B(1,3,2,1)=A(I,20)
      B(3,1,1,2)=A(I,20)
      B(2,1,1,3)=A(I,20)
      B(1,2,3,1)=A(I,20)
      B(1,2,1,3)=A(I,20)
      B(2,1,3,1)=A(I,20)
      B(1,3,1,2)=A(I,20)
      B(3,1,2,1)=A(I,20)
      B(1,2,1,2)=A(I,21)
      B(2,1,2,1)=A(I,21)
      B(1,2,2,1)=A(I,21)
      B(2,1,1,2)=A(I,21)
      RETURN
      END
C
C=======================================================================
C
      REAL FUNCTION WASUM(A,B)
C
C----------------------------------------------------------------------
      REAL A(3),B(3)
      WASUM=A(1)*B(1)+A(2)*B(2)+A(3)*B(3)
      RETURN
      END
C
C=======================================================================
C
      REAL FUNCTION WASUM3(AA,B,C)
C
C----------------------------------------------------------------------
      REAL AA(3,3),B(3),C(3)
      INTEGER I,J
      WASUM3=0.
      DO 20, I=1,3
        DO 10, J=1,3
          WASUM3=WASUM3+AA(I,J)*B(I)*C(J)
  10    CONTINUE
  20  CONTINUE
      RETURN
      END
C
C=======================================================================
C
      REAL FUNCTION WASUM5(A,B,C,D,E)
C
C----------------------------------------------------------------------
      REAL A(3,3,3,3),B(3),C(3),D(3),E(3)
      INTEGER I,J,K,L
      WASUM5=0.
      DO 13, I=1,3
        DO 12, J=1,3
          DO 11, K=1,3
            DO 10, L=1,3
              WASUM5=WASUM5+A(I,J,K,L)*B(I)*C(J)*D(K)*E(L)
  10        CONTINUE
  11      CONTINUE
  12    CONTINUE
  13  CONTINUE
      RETURN
      END
C
C=======================================================================
C
      REAL FUNCTION WASUM4(A,B,C,D,J)
C
C----------------------------------------------------------------------
      REAL A(3,3,3,3),B(3),C(3),D(3)
      INTEGER I,J,K,L
      WASUM4=0.
      DO 13, I=1,3
        DO 11, K=1,3
          DO 10, L=1,3
            WASUM4=WASUM4+A(I,J,K,L)*B(I)*C(K)*D(L)
  10      CONTINUE
  11    CONTINUE
  13  CONTINUE
      RETURN
      END
C
C=======================================================================
C