C XDUDW     SOURCE    FANDEUR   13/01/29    21:16:34     7683
      SUBROUTINE XDUDW(FN,FM,GR,PG,XYZ,HR,PGSQ,RPG,AJ,
     & NES,IDIM,NP,MP,NPG,IAXI,NINC,
     & COES,IK1,S,INDGS,IKS,
     & LE,NBEL,K0,XCOOR,AF,AS,CT,PQ,
     & AF1,AF2,AF3,AF4,AF5,AF6,AF7,AF8,AF9,
     & S2,NPT,IPAD,LEP,IPAP,SPQ)

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)

C************************************************************************
C
C   OPERATEUR DUDW
C
C      CALCULE L'OPERATEUR DE PENALISATION  DIV(U)=EPS*P
C
C      OPTIONS : SOURCE    DIV(U)-Q=EPS*P
C
C
C************************************************************************

      DIMENSION XYZ(IDIM,NP),FN(NP,NPG),GR(IDIM,NP,NPG),PG(NPG)
      DIMENSION FM(MP,NPG),HR(NES,NP,NPG)
      DIMENSION PGSQ(NPG),RPG(NPG),AJ(IDIM,IDIM,NPG)
      DIMENSION XCOOR(*)
      DIMENSION COES(1),LE(NP,NBEL),LEP(MP,*)
      DIMENSION AF(NP,NP,NINC,NINC),CT(MP,NP,NINC),AS(NP,NINC),PQ(MP)
      DIMENSION AF1(NBEL,NP,NP),AF2(NBEL,NP,NP),AF3(NBEL,NP,NP)
      DIMENSION AF4(NBEL,NP,NP),AF5(NBEL,NP,NP),AF6(NBEL,NP,NP)
      DIMENSION AF7(NBEL,NP,NP),AF8(NBEL,NP,NP),AF9(NBEL,NP,NP)
      DIMENSION S(*),S2(NPT,NINC),IPAD(*),IPAP(*) ,SPQ(*)
      REAL*8 U,UIM(10),URIM(10)
C*      REAL*8 UJM(10),URJM(10)
C pour pression P0,P1 et P2
-INC CCREEL

C***********************************************************************

      DEUPI=1.D0
      IF(IAXI.NE.0)DEUPI=2.D0*XPI

      NK=K0
      DO 108 KE=1,NBEL
      NK=NK+1
      JC=(1-IK1)*(NK-1)+1
      DO 109 I=1,NP
      J=LE(I,KE)
      DO 109 N=1,IDIM
      XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N)
 109  CONTINUE

      CALL CALJBR
     &(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,AIRE,AJ,SGN)

      DO 31 K=1,NINC
      DO 31 I=1,NP
      DO 39 M=1,MP
      UIM(M)=0.D0
      URIM(M)=0.D0
      DO 33 L=1,NPG
      UIM(M)=UIM(M)+FM(M,L)*HR(K,I,L)*PGSQ(L)*DEUPI*RPG(L)
 33   CONTINUE
      IF(IAXI.NE.0.AND.K.EQ.1)THEN
      DO 34 L=1,NPG
 34   URIM(M)=URIM(M)+FM(M,L)*FN(I,L)*PGSQ(L)*DEUPI
      ENDIF
      CT(M,I,K)=UIM(M)+URIM(M)
C*??    U=U+(UIM(M)+URIM(M))*(UJM(M)+URJM(M))
 39   CONTINUE
 31   CONTINUE

      DO 41 M=1,MP
      PQ(M)=0.D0
      DO 41 L=1,NPG
      PQ(M)=PQ(M)+FM(M,L)*PGSQ(L)*DEUPI*RPG(L)
  41  CONTINUE

      DO 316 I=1,NP
      DO 316 J=1,NP
      DO 316 K=1,NINC
      DO 316 N=1,NINC
      U=0.D0
      DO 315 M=1,MP
C?    U=U+CT(M,I,K)*CT(M,J,N)/PQ(M)
      U=U+CT(M,I,K)*CT(M,J,N)
 315  CONTINUE
      AF(J,I,N,K)=U/COES(JC)
 316  CONTINUE

C
C   CAS DES SOURCES OU PUITS DE MASSE
C
      IF(INDGS.NE.0)THEN
      DO 73 M=1,MP
      J1=IPAP(LEP(M,KE))
      JS=(1-IKS)*(J1-1)+1
      PQ(M)=0.D0
      DO 71 L=1,NPG
      PQ(M)=PQ(M)+FM(M,L)*S(JS)*PGSQ(L)*DEUPI*RPG(L)
  71  CONTINUE
C?    SPQ(J1)=PQ(M)/COES(JC)
      SPQ(J1)=PQ(M)
  73  CONTINUE

      DO 70 K=1,NINC
      DO 70 I=1,NP
      I1=IPAD(LE(I,KE))
      U=0.D0
      DO 72 M=1,MP
      U=U+CT(M,I,K)*PQ(M)
 72   CONTINUE
      S2(I1,K)=S2(I1,K)+U/COES(JC)
  70  CONTINUE
      ENDIF

 107  CONTINUE
C     write(6,*)' KE=',ke,' np=',np,IDIM
      IF(IDIM.EQ.2)THEN
      DO 701 I=1,NP
      DO 701 J=1,NP
      AF1(KE,J,I)=AF(J,I,1,1)
      AF2(KE,J,I)=AF(J,I,2,1)
      AF3(KE,J,I)=AF(J,I,1,2)
      AF4(KE,J,I)=AF(J,I,2,2)
 701  CONTINUE
C     write(6,*)' AF1 '
C     write(6,1002)AF1
C     write(6,*)' AF2 '
C     write(6,1002)AF2
C     write(6,*)' AF3 '
C     write(6,1002)AF3
C     write(6,*)' AF4 '
C     write(6,1002)AF4
      ELSEIF(IDIM.EQ.3)THEN
      DO 702 I=1,NP
      DO 702 J=1,NP
      AF1(KE,J,I)=AF(J,I,1,1)
      AF2(KE,J,I)=AF(J,I,2,1)
      AF3(KE,J,I)=AF(J,I,3,1)
      AF4(KE,J,I)=AF(J,I,1,2)
      AF5(KE,J,I)=AF(J,I,2,2)
      AF6(KE,J,I)=AF(J,I,3,2)
      AF7(KE,J,I)=AF(J,I,1,3)
      AF8(KE,J,I)=AF(J,I,2,3)
      AF9(KE,J,I)=AF(J,I,3,3)
 702  CONTINUE
      ENDIF
 108  CONTINUE

C     write(6,*)' XDUDW FIN '
      RETURN
 1002 FORMAT(8(1X,1PE11.4))
 1001 FORMAT(20(1X,I5))
      END