C XLAPL     SOURCE    CHAT      05/01/13    04:14:37     5004
      SUBROUTINE XLAPL(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,
     & COEF,COG,COES,KITT,KJTT,IK1,
     & LE,NBEL,K0,XCOOR,AIMPL,IKOMP,
     & AF1,AF2,AF3,
     & AS1,AS2,AS3,
     & NINC,IHV,IARG,S2)
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C************************************************************************
C
C     CALCUL MATRICE ELEMENTAIRE DU LAPLACIEN
C
C
C************************************************************************

C     DIMENSION FN(NP,NPG),HR(KES,NP,NPG),PGSQ(NPG),RPG(NPG)
      DIMENSION LE(NP,NBEL)
      DIMENSION XYZ(IDIM,NP),FN(NP,NPG),GR(IDIM,NP,NPG),PG(NPG)
      DIMENSION HR(IDIM,NP,NPG),PGSQ(NPG),RPG(NPG)
      DIMENSION AF1(NBEL,NP,NP),AF2(NBEL,NP,NP),AF3(NBEL,NP,NP)
      DIMENSION AS1(NBEL,NP,NP),AS2(NBEL,NP,NP),AS3(NBEL,NP,NP)
      DIMENSION CL(9),COEF(1),COES(IDIM,IDIM,*)
      DIMENSION AF(4,4,2,2),XCOOR(*)
C
      DIMENSION COE(3,3),S2(IDIM,1),COG(2,1)

C     write(6,*)' XLAPL  KITT,KJTT,IK1 =',KITT,KJTT,IK1
      NK=K0
      DO 108 KE=1,NBEL
      NK=NK+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 CALJBC(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,
     *IDIM,NP,NPG,IAXI,AIRE)
C
      IAX=0
C     WRITE(6,*)
C     WRITE(6,*)' KITT ',KITT,' KJTT ',KJTT
C     if(iax.eq.0)go to 108
      IF(IHV.EQ.1)IAX=IAXI
      IF(KITT.EQ.4)GO TO 70
      IF(KITT.EQ.3)GO TO 80
      IF(KITT.EQ.1)GO TO 80
      IF(KJTT.EQ.4)GO TO 5
      IC=1+(1-IK1)*(NK-1)
C     WRITE(6,*)' IC ',IC ,COEF(IC)
      DO 1 I=1,NP
      DO 1 J=1,NP
      U=0.D0
      DO 2 L=1,NPG
      V=0.D0
      DO 3 N=1,IDIM
      V=V+HR(N,I,L)*HR(N,J,L)
 3    CONTINUE
      U=U+V*PGSQ(L)
 2    CONTINUE
      U=U*COEF(IC)
      AF1(KE,J,I)=AF1(KE,J,I)+U*AIMPL
      AS1(KE,J,I)=AS1(KE,J,I)+U*(AIMPL-1.D0)
 1    CONTINUE

      IF(NINC.GE.2)THEN
      DO 1762 I=1,NP
      DO 1762 J=1,NP
      AF2(KE,I,J)=AF1(KE,I,J)
      AS2(KE,I,J)=AS1(KE,I,J)
 1762 CONTINUE
      ENDIF
      IF(NINC.GE.3)THEN
      DO 1763 I=1,NP
      DO 1763 J=1,NP
      AF3(KE,I,J)=AF1(KE,I,J)
      AS3(KE,I,J)=AS1(KE,I,J)
 1763 CONTINUE
      ENDIF

      IF(IHV.EQ.1)THEN
      IF(IAXI.EQ.1)THEN
      DO 41 I=1,NP
      DO 41 J=1,NP
      U=0.D0
      DO 42 L=1,NPG
      U=U+FN(I,L)*FN(J,L)/RPG(L)/RPG(L)*PGSQ(L)
 42   CONTINUE
      U=U*COEF(IC)
      AF2(KE,J,I)=AF2(KE,J,I)+U*AIMPL
      AS2(KE,J,I)=AS2(KE,J,I)+U*(AIMPL-1.D0)
 41   CONTINUE
      ELSEIF(IAXI.EQ.2)THEN
      DO 43 I=1,NP
      DO 43 J=1,NP
      U=0.D0
      DO 44 L=1,NPG
      U=U+FN(I,L)*FN(J,L)/RPG(L)/RPG(L)*PGSQ(L)
 44   CONTINUE
      U=U*COEF(IC)
      AF1(KE,J,I)=AF1(KE,J,I)+U*AIMPL
      AS1(KE,J,I)=AS1(KE,J,I)+U*(AIMPL-1.D0)
 43   CONTINUE
      ENDIF
      ENDIF
      GO TO 108

 5    CONTINUE
      DO 51 L=1,NPG
      C=0.D0
      DO 52 I=1,NP
      IU=LE(I,KE)
      C=C+COEF(IU)*FN(I,L)
 52   CONTINUE
      CL(L)=C
 51   CONTINUE
C
C     write(6,*)' BCL 15 '
      DO 15 I=1,NP
      DO 15 J=I,NP
      U=0.D0
      DO 12 L=1,NPG
      V=0.D0
      DO 13 N=1,IDIM
      V=V+HR(N,I,L)*HR(N,J,L)
 13   CONTINUE
      U=U+V*PGSQ(L)*CL(L)
 12   CONTINUE
      DO 14 K=1,NINC
      AF(J,I,K,K)=AF(J,I,K,K)+U
      IF(I.NE.J)AF(I,J,K,K)=AF(I,J,K,K)+U
 14   CONTINUE
      IF(IAX.EQ.0)GO TO 15
      U=0.D0
      DO 410 L=1,NPG
      U=U+FN(I,L)*FN(J,L)/RPG(L)/RPG(L)*PGSQ(L)*CL(L)
 410  CONTINUE
      KU=3-IAX
      AF(J,I,KU,KU)=AF(J,I,KU,KU)+U
      IF(I.NE.J)AF(I,J,KU,KU)=AF(I,J,KU,KU)+U
 15   CONTINUE
      GO TO 108
C
C
 70   CONTINUE
      IF(NINC.NE.1)CALL ARRET(0)
      IF(KJTT.EQ.4)CALL ARRET(0)
      IC=1+(1-IK1)*(NK-1)
      DO 71 I=1,NP
      DO 71 J=1,NP
      U=0.D0
      DO 72 L=1,NPG
      UL=0.D0
      DO 73 M=1,IDIM
      UN=0.D0
      DO 74 N=1,IDIM
      UN=UN+COES(N,M,IC)*HR(N,J,L)
 74   CONTINUE
      UL=UL+UN*HR(M,I,L)
 73   CONTINUE
      U=U+UL*PGSQ(L)
 72   CONTINUE
      AF(J,I,1,1)=AF(J,I,1,1)+U
 71   CONTINUE
      GO TO 108
C
C
   80 CONTINUE
      IF(KJTT.EQ.4)CALL ARRET(0)
      IF(IARG.EQ.0)CALL ARRET(0)
      IC=1+(1-IK1)*(NK-1)
      ALT=COG(1,IC)-COG(2,IC)
      AT=COG(2,IC)
C     WRITE(6,*)' ALT ',ALT,' AT ',AT
      IF(IDIM.EQ.2) THEN
      CO=SQRT(S2(1,IC)*S2(1,IC)+S2(2,IC)*S2(2,IC))
      ENDIF
      IF(IDIM.EQ.3) THEN
      CO=SQRT(S2(1,IC)*S2(1,IC)+S2(2,IC)*S2(2,IC)+S2(3,IC)*S2(3,IC))
      ENDIF
C     WRITE(6,*)' CO ',CO
C     WRITE(6,*)' S2 ',S2(1,IC),S2(2,IC)
      IF(IDIM.EQ.2) THEN
      COE(1,1)=ALT*S2(1,IC)*S2(1,IC)/CO+AT*CO
      COE(2,2)=ALT*S2(2,IC)*S2(2,IC)/CO+AT*CO
      COE(1,2)=ALT*S2(1,IC)*S2(2,IC)/CO
      COE(2,1)=COE(1,2)
      ENDIF
      IF(IDIM.EQ.3) THEN
      COE(1,1)=ALT*S2(1,IC)*S2(1,IC)/CO+AT*CO
      COE(2,2)=ALT*S2(2,IC)*S2(2,IC)/CO+AT*CO
      COE(3,3)=ALT*S2(3,IC)*S2(3,IC)/CO+AT*CO
      COE(1,2)=ALT*S2(1,IC)*S2(2,IC)/CO
      COE(1,3)=ALT*S2(1,IC)*S2(3,IC)/CO
      COE(2,3)=ALT*S2(2,IC)*S2(3,IC)/CO
      COE(2,1)=COE(1,2)
      COE(3,1)=COE(1,3)
      COE(3,2)=COE(2,3)
      ENDIF
C     WRITE(6,*)' COE ',COE(1,1),COE(1,2),COE(1,3)
C     WRITE(6,*)'     ',COE(2,1),COE(2,2),COE(2,3)
C     WRITE(6,*)'     ',COE(3,1),COE(3,2),COE(3,3)
C
      DO 81 I=1,NP
      DO 81 J=1,NP
      U=0.D0
      DO 82 L=1,NPG
      UL=0.D0
      DO 83 M=1,IDIM
      UN=0.D0
      DO 84 N=1,IDIM
      UN=UN+COE(N,M)*HR(N,J,L)
 84   CONTINUE
      UL=UL+UN*HR(M,I,L)
 83   CONTINUE
      U=U+UL*PGSQ(L)
 82   CONTINUE
      AF(J,I,1,1)=AF(J,I,1,1)+U
 81   CONTINUE

 108  CONTINUE

C     write(6,*)' FIN XLAPL'

      RETURN
 1002 FORMAT(10(1X,1PE11.4))
 1001 FORMAT(20(1X,I5))
      END