C ZCVI      SOURCE    CHAT      05/01/13    04:21:47     5004           
      SUBROUTINE ZCVI
C
C
C     VERSION VECTORISEE
C
C     Les {l{ments sont group{s en paquets de LRV {l{ments, LRV {tant
C     la longueur des registres vectoriels de la machine cible, i.e
C     64 sur Cray, 128 ou 256 sur IBM 3090VF. On prom}ne une fenetre
C     de longueur LRV sur la boucle g{n{rale de longueur NEL.
C
C
     & (HR,RPG,DRR,LE,NEL,K0,NPTD,IES,NP,IAXI,IKOMP,IKAS,
     & COEFF,IK1,RGE,IKG,NELG,TN,IKT,TREF,IKREF,IPADS,
     & UN,IPADU,NPTU,GN,F,IPADI,VF,IPADF,NPTF,
     & VOLU,COTE,NELZ,IDCEN,IPG,
     & DTM1,DT,DTT1,DTT2,NUEL,DIAEL,FN)

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

C***********************************************************************
C
C  CE SP DISCRETISE LES EQUATIONS DE NAVIER STOKES
C    EN 2D SUR LES ELEMENTS QUA4 ET TRI3      PLAN OU AXI
C    EN 3D SUR LES ELEMENTS CUB8 ET PRI6
C  LES OPERATEURS SONT "SOUS-INTEGRES"
C
C  SYNTAXE :
C
C         NS(NU,UN,RGE,DE)  INCO GN :
C
C  COEFFICIENTS :
C  --------------
C
C  UN(NPTD,IES)       CHAMPS DE VITESSE TRANSPORTANT
C  COEFF(SCAL DOMA)  VISCOSITE CINEMATIQE MOLECULAIRE( NU )
C       (SCAL ELEM)
C   RGE(NELG,IES)       TERME SOURCE
C
C
C  INCONNUES :
C  -----------
C
C   UN(NPTD,IES)        CHAMPS DE VITESSE TRANSPORTANT
C   GN(NPTD,IES)        CHAMPS DE VITESSE TRANSPORTE
C   VN(NPTD,IES)        CHAMPS DE VITESSE DU FLUIDE
C
C
C
C
C  TABLEAUX DE TRAVAIL :
C  ---------------------        N     N-1
C                     (D + D1) U - D U              N-1    T  N
C                     -------------------  = F - A U    - C  P
C                               DT
C
C                              N-1
C   F(NPTD,IES)   CONTIENT   A U   - F     (VITESSE)
C
C
C
C***********************************************************************

-INC CCVQUA4
-INC CCREEL
*- 
-INC CCOPTIO
-INC SMCOORD
C
C Longueur des registres vectoriels de la machine cible
C On prend 64 pour ne pas augmenter la taille des tableaux
C n{cessaires @ la vectorisation.
C
      PARAMETER(LRV=64)

      DIMENSION UN(NPTU,IES),GN(NPTD,IES),VF(NPTF,IES)
      DIMENSION TN(*),TREF(*)
      DIMENSION COEFF(*),RGE(NELG,IES)
      DIMENSION COTE(NELZ,IES),VOLU(NELZ),KLIP(100)

      DIMENSION IPADI(*),LE(NP,1),IPADU(*),IPADF(*),IPADS(*)
      DIMENSION HR(NEL,NP,IES),RPG(1),DRR(NP,NEL)

      DIMENSION QGGT(8,8),Q1(8,8),Q2(8,8),Q3(8,8)

      DIMENSION COEF(LRV),AIRE(LRV)
      DIMENSION WX(LRV,9),WY(LRV,9),WZ(LRV,9)
      DIMENSION AL(LRV),AH(LRV),AP(LRV)
C UIX,... vitesse transportante
      DIMENSION UIX(LRV,9),UIY(LRV,9),UIZ(LRV,9)
C GIX,... vitesse massique ou transportée ou inconnue du/dt
      DIMENSION GIX(LRV,9),GIY(LRV,9),GIZ(LRV,9)
C VIX,... vitesse du fluide
      DIMENSION VIX(LRV,9),VIY(LRV,9),VIZ(LRV,9)

      DIMENSION UMI(LRV,3),VMI(LRV,3)
      DIMENSION COEFT(LRV),RGX(LRV),RGY(LRV),RGZ(LRV)
      DIMENSION TO1(LRV),TO2X(LRV),TO2Y(LRV)
      DIMENSION SAF1(LRV,9),SAF2(LRV,9),SAF3(LRV,9)
      DIMENSION CHGLD(LRV),CHGLPX(LRV),CHGLPY(LRV),CHGLPZ(LRV)

      DIMENSION F(NPTD,*),FN(NP,*)

      SAVE IPAS,QGGT,Q1,Q2,Q3
      DATA IPAS/0/

C
C           INITIALISATIONS DIVERSES
C
C     write(6,*)' DEBUT YCVI ',' IKAS=',ikas,' IKOMP=',ikomp,idcen
C     write(6,*)' NPTD=',nptd,IPAS
      IF(IPAS.EQ.0)CALL CALHRH(QGGT,Q1,Q2,Q3,IES)

C                           ********
C                           *  2D  *
C                           ********

      IF(IES.EQ.3)GO TO 10

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                           DIFFERENCES TRIANGLE / QUADRANGLE
      IF(NP.EQ.4)THEN
      QUA4=1.D0
      ELSE
      QUA4=0.D0
      ENDIF
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C Calcul du nombre de paquets de LRV {l{ments
C
      NNN=MOD(NEL,LRV)
      IF(NNN.EQ.0) NPACK=NEL/LRV
      IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV
      KPACKD=1
      KPACKF=NPACK
C
C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS **********
C
C     write(6,*)' DEBUT YCVI 7001'
      DO 7001 KPACK=KPACKD,KPACKF
C
C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS =======
C
C 1. Calcul des limites du paquet courant.
      KDEB=1+(KPACK-1)*LRV
      KFIN=MIN(NEL,KDEB+LRV-1)
C
C Ben voil@, on peut y aller ... i.e. traiter le paquet courant.
C
      DO 7002 K=KDEB,KFIN
      KP=K-KDEB+1
      NK=K+K0
      K1=1+(1-IK1)*(NK-1)
      COEF(KP)=COEFF(K1)
      AIRE(KP)=VOLU(NK)
      AL(KP)=COTE(NK,1)+XPETIT
      AH(KP)=COTE(NK,2)+XPETIT
 7002 CONTINUE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DO 7006 I=1,NP
      DO 7006 K=KDEB,KFIN
      KP=K-KDEB+1
      NU=IPADU(LE(I,K))
      NG=IPADI(LE(I,K))
      NF=IPADF(LE(I,K))
      UIX(KP,I)=UN(NU,1)
      UIY(KP,I)=UN(NU,2)
      GIX(KP,I)=GN(NG,1)
      GIY(KP,I)=GN(NG,2)
      VIX(KP,I)=VF(NF,1)
      VIY(KP,I)=VF(NF,2)
 7006 CONTINUE

      CALL KSUPG1(WX,UMI,CHGLD,CHGLPX,KDEB,KFIN,LRV,
     &GN(1,1),IPADI,UN,COEF,NPTD,NEL,NP,DRR,HR,FN,
     &AIRE,AL,AH,AP,IDCEN,IPADU,LE,QUA4,IKOMP,
     &DTM1,DT,DTT1,DTT2,DIAEL,NUEL)

      CALL KSUPG1(WY,UMI,CHGLD,CHGLPY,KDEB,KFIN,LRV,
     &GN(1,2),IPADI,UN,COEF,NPTD,NEL,NP,DRR,HR,FN,
     &AIRE,AL,AH,AP,IDCEN,IPADU,LE,QUA4,IKOMP,
     &DTM1,DT,DTT1,DTT2,DIAEL,NUEL)

C
C Initialisation des variables d'accumulation SAF1,SAF2,SBF
C

      IF(IAXI.NE.0)THEN
       DO 70050 N=1,IDIM
       DO 70050 K=KDEB,KFIN
       KP=K-KDEB+1
       VMI(KP,N)=XPETIT
70050 CONTINUE
       DO 70051 N=1,IDIM
       DO 70051 I=1,NP
       DO 70051 K=KDEB,KFIN
       KP=K-KDEB+1
       NF=IPADF(LE(I,K))
       VMI(KP,N)=VMI(KP,N)+VF(NF,N)*DRR(I,K)
70051 CONTINUE
       DO 70052 K=KDEB,KFIN
       KP=K-KDEB+1
       VMI(KP,1)=VMI(KP,1)/AIRE(KP)
       VMI(KP,2)=VMI(KP,2)/AIRE(KP)
70052 CONTINUE
      ENDIF

      IF(IKOMP.EQ.0)THEN

       IF(IKAS.EQ.1)THEN
       DO 70061 I=1,NP
       DO 70061 K=KDEB,KFIN
       KP=K-KDEB+1
       SAF1(KP,I)=0.D0
       SAF2(KP,I)=0.D0
70061  CONTINUE
       ELSEIF(IKAS.EQ.2)THEN
       DO 70021 K=KDEB,KFIN
       KP=K-KDEB+1
       NK=K+K0
       NKG=1+(1-IKG)*(NK-1)
       RGX(KP)=RGE(NKG,1)
       RGY(KP)=RGE(NKG,2)
70021  CONTINUE

       IF(IPG.EQ.0)THEN
       DO 70062 I=1,NP
       DO 70062 K=KDEB,KFIN
       KP=K-KDEB+1
       SAF1(KP,I)=(-RGX(KP))*DRR(I,K)
       SAF2(KP,I)=(-RGY(KP))*DRR(I,K)
70062  CONTINUE
       ELSE
       DO 71062 I=1,NP
       DO 71062 K=KDEB,KFIN
       KP=K-KDEB+1
       SAF1(KP,I)=(-RGX(KP))*WX(KP,I)
       SAF2(KP,I)=(-RGY(KP))*WY(KP,I)
71062  CONTINUE
       ENDIF

       ELSEIF(IKAS.EQ.4)THEN
       DO 70022 K=KDEB,KFIN
       KP=K-KDEB+1
       NK=K+K0
       NKG=1+(1-IKG)*(NK-1)
       RGX(KP)=RGE(NKG,1)
       RGY(KP)=RGE(NKG,2)
70022  CONTINUE

       IF(IPG.EQ.0)THEN
       DO 70063 I=1,NP
       DO 70063 K=KDEB,KFIN
       KP=K-KDEB+1
       NF=1+(1-IKT)*(IPADS(LE(I,K))-1)
       NFR=1+(1-IKREF)*(IPADS(LE(I,K))-1)
C?     write(6,*)' NF=',NF,' NFR=',nfr
       SAF1(KP,I)=(-RGX(KP)*(TREF(NFR)-TN(NF)))*DRR(I,K)
       SAF2(KP,I)=(-RGY(KP)*(TREF(NFR)-TN(NF)))*DRR(I,K)
70063  CONTINUE
       ELSE
       DO 71063 I=1,NP
       DO 71063 K=KDEB,KFIN
       KP=K-KDEB+1
       NF=1+(1-IKT)*(IPADS(LE(I,K))-1)
       NFR=1+(1-IKREF)*(IPADS(LE(I,K))-1)
       SAF1(KP,I)=(-RGX(KP)*(TREF(NFR)-TN(NF)))*WX(KP,I)
       SAF2(KP,I)=(-RGY(KP)*(TREF(NFR)-TN(NF)))*WY(KP,I)
71063  CONTINUE
       ENDIF

       ENDIF

      ELSEIF(IKOMP.EQ.1)THEN

       IF(IKAS.EQ.2)THEN
       DO 70064 I=1,NP
       DO 70064 K=KDEB,KFIN
       KP=K-KDEB+1
       SAF1(KP,I)=0.D0
       SAF2(KP,I)=0.D0
70064  CONTINUE

       ELSEIF(IKAS.EQ.3)THEN
       DO 70024 K=KDEB,KFIN
       KP=K-KDEB+1
       NK=K+K0
       NKG=1+(1-IKG)*(NK-1)
       RGX(KP)=RGE(NKG,1)
       RGY(KP)=RGE(NKG,2)
70024  CONTINUE

       IF(IPG.EQ.0)THEN
       DO 70065 I=1,NP
       DO 70065 K=KDEB,KFIN
       KP=K-KDEB+1
       SAF1(KP,I)=-RGX(KP)*DRR(I,K)
       SAF2(KP,I)=-RGY(KP)*DRR(I,K)
70065  CONTINUE
       ELSE
       DO 71065 I=1,NP
       DO 71065 K=KDEB,KFIN
       KP=K-KDEB+1
       SAF1(KP,I)=-RGX(KP)*WX(KP,I)
       SAF2(KP,I)=-RGY(KP)*WY(KP,I)
71065  CONTINUE
       ENDIF

       ENDIF

      ENDIF

C Le coeur du calcul ...

      IF(IKOMP.EQ.0)THEN

      DO 7014 I=1,NP
      DO 7014 J= 1,NP
      DO 7014 K=KDEB,KFIN
      KP=K-KDEB+1

      ZVGGX=AIRE(KP)*CHGLPX(KP)*VGGT(J,I)

      ZVGGY=AIRE(KP)*CHGLPY(KP)*VGGT(J,I)

      ZVGT=AIRE(KP)*(
     &     HR(K,I,1)*HR(K,J,1)*COEF(KP)
     &+    HR(K,I,2)*HR(K,J,2)*COEF(KP)
     &+    CHGLD(KP)*VGGT(J,I) )

      V2=UMI(KP,1)*HR(K,J,1)+UMI(KP,2)*HR(K,J,2)

      SAF1(KP,I)=SAF1(KP,I)+(V2*WX(KP,I)+ZVGGX)*GIX(KP,J)+ZVGT*VIX(KP,J)
      SAF2(KP,I)=SAF2(KP,I)+(V2*WY(KP,I)+ZVGGY)*GIY(KP,J)+ZVGT*VIY(KP,J)

 7014 CONTINUE

      ELSEIF(IKOMP.EQ.1)THEN

      DO 7015 I=1,NP
      DO 7015 J= 1,NP
      DO 7015 K=KDEB,KFIN
      KP=K-KDEB+1
      ZVGGX=AIRE(KP)*CHGLPX(KP)*VGGT(J,I)

      ZVGGY=AIRE(KP)*CHGLPY(KP)*VGGT(J,I)

      ZVGT=AIRE(KP)*(
     &     HR(K,I,1)*HR(K,J,1)*COEF(KP)
     &+    HR(K,I,2)*HR(K,J,2)*COEF(KP)
     &+    CHGLD(KP)*VGGT(J,I) )

      COEFT3=COEF(KP)/3.D0
      ZVGUU=AIRE(KP)*(1.D0/AL(KP)/AL(KP)/12.D0
     &     *VGGT(J,I)*QUA4  )*COEFT3*UIX(KP,J)
      ZVGUV=AIRE(KP)*( HR(K,I,1)*HR(K,J,2))*COEFT3*UIY(KP,J)
      ZVGVU=AIRE(KP)*( HR(K,I,2)*HR(K,J,1))*COEFT3*UIX(KP,J)
      ZVGVV=AIRE(KP)*(1.D0/AH(KP)/AH(KP)/12.D0
     &     *VGGT(J,I)*QUA4  )*COEFT3*UIY(KP,J)

      V2=UIX(KP,J)*HR(K,J,1)+UIY(KP,J)*HR(K,J,2)

      SAF1(KP,I)=SAF1(KP,I)+(V2*WX(KP,I)+ZVGGX)*GIX(KP,J)+ZVGT*UIX(KP,J)
     &+ ZVGUU + ZVGUV
      SAF2(KP,I)=SAF2(KP,I)+(V2*WY(KP,I)+ZVGGY)*GIY(KP,J)+ZVGT*UIY(KP,J)
     &+ ZVGVU + ZVGVV

 7015 CONTINUE

      ENDIF

      IF(IAXI.NE.0) THEN

         DO 7016 I=1,NP
         DO 7016 K=KDEB,KFIN
         KP=K-KDEB+1

         R2=1.D0/RPG(K)/RPG(K)*WX(KP,I)
        SAF1(KP,I)=SAF1(KP,I)+R2*COEF(KP)*VMI(KP,1)

 7016    CONTINUE

         IF(IKOMP.EQ.1)THEN
         DO 7118 I=1,NP
         DO 7118 K=KDEB,KFIN
         KP=K-KDEB+1
         R1=1.D0/RPG(K)*WY(KP,I)
         SAF1(KP,I)=SAF1(KP,I)+R1*UMI(KP,1)*GIX(KP,I)
         SAF2(KP,I)=SAF2(KP,I)+R1*UMI(KP,1)*GIY(KP,I)
 7118    CONTINUE
        ENDIF

      ENDIF
C
C     Fin de l'accumulation dans SAF1,SAF2.
C     On ajoute ces incr{ments @ F.
C
      DO 7017 I=1,NP
      DO 7017 K=KDEB,KFIN
      KP=K-KDEB+1
      NF=IPADI(LE(I,K))
      F(NF,1)=F(NF,1)+SAF1(KP,I)
      F(NF,2)=F(NF,2)+SAF2(KP,I)
 7017 CONTINUE

 1960 FORMAT(/,' ***** SUB XCVTIT : IPAT=',I5,'  K=',I5,' *****')
 1961 FORMAT(2X,I5,' * ',4(1X,I5))
 1962 FORMAT(2X,8(1X,1PE11.4))
 1964 FORMAT(4(1X,1PE11.4))
 7001 CONTINUE

C     WRITE(6,*)' ********** FIN YCVI 2D *****************'
C     CALL ARRET(0)
      IPAS=1
      RETURN

C                           ********
C                           *  3D  *
C                           ********

 10   CONTINUE

C::::::BENET:::SUPPRESION CORRECTION HOURGLASS POUR LES PRISME::29:01:91

      CUB8=0.D0
      IF(NP.EQ.8)CUB8=1.D0

C
C Calcul du nombre de paquets de LRV {l{ments
C
      NNN=MOD(NEL,LRV)
      IF(NNN.EQ.0) NPACK=NEL/LRV
      IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV
      KPACKD=1
      KPACKF=NPACK
C
C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS **********
C
      DO 8001 KPACK=KPACKD,KPACKF
C
C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS =======
C
C 1. Calcul des limites du paquet courant.
      KDEB=1+(KPACK-1)*LRV
      KFIN=MIN(NEL,KDEB+LRV-1)
C
C Ben voil@, on peut y aller ... i.e. traiter le paquet courant.
C
      DO 8002 K=KDEB,KFIN
      KP=K-KDEB+1
      NK=K+K0
      K1=1+(1-IK1)*(NK-1)
      COEF(KP)=COEFF(K1)
      AIRE(KP)=VOLU(NK)
      AL(KP)=COTE(NK,1)+XPETIT
      AH(KP)=COTE(NK,2)+XPETIT
      AP(KP)=COTE(NK,3)+XPETIT
 8002 CONTINUE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Initialisation des UMI avant accumulation

      DO 8006 I=1,NP
      DO 8006 K=KDEB,KFIN
      KP=K-KDEB+1
      NU=IPADU(LE(I,K))
      NF=IPADF(LE(I,K))
      NG=IPADI(LE(I,K))
      UIX(KP,I)=UN(NU,1)
      UIY(KP,I)=UN(NU,2)
      UIZ(KP,I)=UN(NU,3)
      VIX(KP,I)=VF(NF,1)
      VIY(KP,I)=VF(NF,2)
      VIZ(KP,I)=VF(NF,3)
      GIX(KP,I)=GN(NG,1)
      GIY(KP,I)=GN(NG,2)
      GIZ(KP,I)=GN(NG,3)
 8006 CONTINUE

C     write(6,*)'****************************'
C     write(6,*)' DT,DTT1,DTT2=',DT,DTT1,DTT2
      CALL KSUPG1(WX,UMI,CHGLD,CHGLPX,KDEB,KFIN,LRV,
     &GN(1,1),IPADI,UN,COEF,NPTD,NEL,NP,DRR,HR,FN,
     &AIRE,AL,AH,AP,IDCEN,IPADU,LE,CUB8,IKOMP,
     &DTM1,DT,DTT1,DTT2,DIAEL,NUEL)
C     write(6,1002)umi

C     write(6,*)' DT,DTT1,DTT2=',DT,DTT1,DTT2
      CALL KSUPG1(WY,UMI,CHGLD,CHGLPY,KDEB,KFIN,LRV,
     &GN(1,2),IPADI,UN,COEF,NPTD,NEL,NP,DRR,HR,FN,
     &AIRE,AL,AH,AP,IDCEN,IPADU,LE,CUB8,IKOMP,
     &DTM1,DT,DTT1,DTT2,DIAEL,NUEL)

C     write(6,*)' DT,DTT1,DTT2=',DT,DTT1,DTT2
      CALL KSUPG1(WZ,UMI,CHGLD,CHGLPZ,KDEB,KFIN,LRV,
     &GN(1,3),IPADI,UN,COEF,NPTD,NEL,NP,DRR,HR,FN,
     &AIRE,AL,AH,AP,IDCEN,IPADU,LE,CUB8,IKOMP,
     &DTM1,DT,DTT1,DTT2,DIAEL,NUEL)
C     write(6,*)' DT,DTT1,DTT2=',DT,DTT1,DTT2

C
C Initialisation des variables d'accumulation SAF1,SAF2,SBF
C
      IF(IKOMP.EQ.0)THEN

       IF(IKAS.EQ.1)THEN
       DO 80061 I=1,NP
       DO 80061 K=KDEB,KFIN
       KP=K-KDEB+1
       SAF1(KP,I)=0.D0
       SAF2(KP,I)=0.D0
       SAF3(KP,I)=0.D0
80061  CONTINUE
       ELSEIF(IKAS.EQ.2)THEN
       DO 80021 K=KDEB,KFIN
       KP=K-KDEB+1
       NK=K+K0
       NKG=1+(1-IKG)*(NK-1)
       RGX(KP)=RGE(NKG,1)
       RGY(KP)=RGE(NKG,2)
       RGZ(KP)=RGE(NKG,3)
80021  CONTINUE

       IF(IPG.EQ.0)THEN
       DO 80062 I=1,NP
       DO 80062 K=KDEB,KFIN
       KP=K-KDEB+1
       SAF1(KP,I)=(-RGX(KP))*DRR(I,K)
       SAF2(KP,I)=(-RGY(KP))*DRR(I,K)
       SAF3(KP,I)=(-RGZ(KP))*DRR(I,K)
80062  CONTINUE
       ELSE
       DO 81062 I=1,NP
       DO 81062 K=KDEB,KFIN
       KP=K-KDEB+1
       SAF1(KP,I)=(-RGX(KP))*WX(KP,I)
       SAF2(KP,I)=(-RGY(KP))*WY(KP,I)
       SAF3(KP,I)=(-RGZ(KP))*WZ(KP,I)
81062  CONTINUE
       ENDIF

       ELSEIF(IKAS.EQ.4)THEN
       DO 80022 K=KDEB,KFIN
       KP=K-KDEB+1
       NK=K+K0
       NKG=1+(1-IKG)*(NK-1)
       RGX(KP)=RGE(NKG,1)
       RGY(KP)=RGE(NKG,2)
       RGZ(KP)=RGE(NKG,3)
80022  CONTINUE

       IF(IPG.EQ.0)THEN
       DO 80063 I=1,NP
       DO 80063 K=KDEB,KFIN
       KP=K-KDEB+1
       NF=1+(1-IKT)*(IPADS(LE(I,K))-1)
       NFR=1+(1-IKREF)*(IPADS(LE(I,K))-1)
       SAF1(KP,I)=(-RGX(KP)*(TREF(NFR)-TN(NF)))*DRR(I,K)
       SAF2(KP,I)=(-RGY(KP)*(TREF(NFR)-TN(NF)))*DRR(I,K)
       SAF3(KP,I)=(-RGZ(KP)*(TREF(NFR)-TN(NF)))*DRR(I,K)
80063  CONTINUE
       ELSE
       DO 81063 I=1,NP
       DO 81063 K=KDEB,KFIN
       KP=K-KDEB+1
       NF=1+(1-IKT)*(IPADS(LE(I,K))-1)
       NFR=1+(1-IKREF)*(IPADS(LE(I,K))-1)
       SAF1(KP,I)=(-RGX(KP)*(TREF(NFR)-TN(NF)))*WX(KP,I)
       SAF2(KP,I)=(-RGY(KP)*(TREF(NFR)-TN(NF)))*WY(KP,I)
       SAF3(KP,I)=(-RGZ(KP)*(TREF(NFR)-TN(NF)))*WZ(KP,I)
81063  CONTINUE
       ENDIF

       ENDIF

      ELSEIF(IKOMP.EQ.1)THEN

       IF(IKAS.EQ.2)THEN
       DO 80064 I=1,NP
       DO 80064 K=KDEB,KFIN
       KP=K-KDEB+1
       SAF1(KP,I)=0.D0
       SAF2(KP,I)=0.D0
       SAF3(KP,I)=0.D0
80064  CONTINUE

       ELSEIF(IKAS.EQ.3)THEN
       DO 80024 K=KDEB,KFIN
       KP=K-KDEB+1
       NK=K+K0
       NKG=1+(1-IKG)*(NK-1)
       RGX(KP)=RGE(NKG,1)
       RGY(KP)=RGE(NKG,2)
       RGZ(KP)=RGE(NKG,3)
80024  CONTINUE

       IF(IPG.EQ.0)THEN
       DO 80065 I=1,NP
       DO 80065 K=KDEB,KFIN
       KP=K-KDEB+1
       SAF1(KP,I)=-RGX(KP)*DRR(I,K)
       SAF2(KP,I)=-RGY(KP)*DRR(I,K)
       SAF3(KP,I)=-RGZ(KP)*DRR(I,K)
80065  CONTINUE
       ELSE
       DO 81065 I=1,NP
       DO 81065 K=KDEB,KFIN
       KP=K-KDEB+1
       SAF1(KP,I)=-RGX(KP)*WX(KP,I)
       SAF2(KP,I)=-RGY(KP)*WY(KP,I)
       SAF3(KP,I)=-RGZ(KP)*WZ(KP,I)
81065  CONTINUE
       ENDIF

       ENDIF

      ENDIF

C Le coeur du calcul ...

      IF(IKOMP.EQ.0)THEN

      DO 8014 I=1,NP
      DO 8014 J= 1,NP
      DO 8014 K=KDEB,KFIN
      KP=K-KDEB+1

      ZVGGX=AIRE(KP)*CHGLPX(KP)*QGGT(J,I)

      ZVGGY=AIRE(KP)*CHGLPY(KP)*QGGT(J,I)

      ZVGGZ=AIRE(KP)*CHGLPZ(KP)*QGGT(J,I)

      ZVGT=AIRE(KP)*(
     &     HR(K,I,1)*HR(K,J,1)*COEF(KP)
     &+    HR(K,I,2)*HR(K,J,2)*COEF(KP)
     &+    HR(K,I,3)*HR(K,J,3)*COEF(KP) )
     &+    CHGLD(KP)*QGGT(J,I)

      V2=UMI(KP,1)*HR(K,J,1)+UMI(KP,2)*HR(K,J,2)
     &  +UMI(KP,3)*HR(K,J,3)

      SAF1(KP,I)=SAF1(KP,I)+(V2*WX(KP,I)+ZVGGX)*GIX(KP,J)+ZVGT*VIX(KP,J)
      SAF2(KP,I)=SAF2(KP,I)+(V2*WY(KP,I)+ZVGGY)*GIY(KP,J)+ZVGT*VIY(KP,J)
      SAF3(KP,I)=SAF3(KP,I)+(V2*WZ(KP,I)+ZVGGZ)*GIZ(KP,J)+ZVGT*VIZ(KP,J)

 8014 CONTINUE

      ELSEIF(IKOMP.EQ.1)THEN

      DO 8015 I=1,NP
      DO 8015 J= 1,NP
      DO 8015 K=KDEB,KFIN
      KP=K-KDEB+1

      ZVGGX=AIRE(KP)*CHGLPX(KP)*QGGT(J,I)

      ZVGGY=AIRE(KP)*CHGLPY(KP)*QGGT(J,I)

      ZVGGZ=AIRE(KP)*CHGLPZ(KP)*QGGT(J,I)

      ZVGT=AIRE(KP)*(
     &     HR(K,I,1)*HR(K,J,1)*COEF(KP)
     &+    HR(K,I,2)*HR(K,J,2)*COEF(KP)
     &+    HR(K,I,3)*HR(K,J,3)*COEF(KP) )
     &+    CHGLD(KP)*QGGT(J,I)

      COEFT3=COEF(KP)/3.D0
      GEO1=0.D0
      ZVGUU=(AIRE(KP)*(HR(K,I,1)*HR(K,J,1))+GEO1)*COEFT3*UIX(KP,J)
      ZVGUV=AIRE(KP)*( HR(K,I,1)*HR(K,J,2))*COEFT3*UIY(KP,J)
      ZVGUW=AIRE(KP)*( HR(K,I,1)*HR(K,J,3))*COEFT3*UIZ(KP,J)
      ZVGVU=AIRE(KP)*( HR(K,I,2)*HR(K,J,1))*COEFT3*UIX(KP,J)
      ZVGVV=(AIRE(KP)*(HR(K,I,2)*HR(K,J,2))+GEO1)*COEFT3*UIY(KP,J)
      ZVGVW=AIRE(KP)*( HR(K,I,2)*HR(K,J,3))*COEFT3*UIZ(KP,J)
      ZVGWU=AIRE(KP)*( HR(K,I,3)*HR(K,J,1))*COEFT3*UIX(KP,J)
      ZVGWV=AIRE(KP)*( HR(K,I,3)*HR(K,J,2))*COEFT3*UIY(KP,J)
      ZVGWW=(AIRE(KP)*(HR(K,I,3)*HR(K,J,3))+GEO1)*COEFT3*UIZ(KP,J)

      V2=(UIX(KP,J)*HR(K,J,1)+UIY(KP,J)*HR(K,J,2)+UIZ(KP,J)*HR(K,J,3))
     &   *WZ(KP,I)

      SAF1(KP,I)=SAF1(KP,I)+(V2+ZVGGX)*GIX(KP,J)+ZVGT*UIX(KP,J)
     &+ ZVGUU + ZVGUV + ZVGUW
      SAF2(KP,I)=SAF2(KP,I)+(V2+ZVGGY)*GIY(KP,J)+ZVGT*UIY(KP,J)
     &+ ZVGVU + ZVGUV + ZVGVW
      SAF3(KP,I)=SAF3(KP,I)+(V2+ZVGGZ)*GIZ(KP,J)+ZVGT*UIZ(KP,J)
     &+ ZVGWU + ZVGWV + ZVGWW

 8015 CONTINUE

      ENDIF

C
C     Fin de l'accumulation dans SAF1,SAF2.
C     On ajoute ces incr{ments @ F.
C
      DO 8017 I=1,NP
      DO 8017 K=KDEB,KFIN
      KP=K-KDEB+1
      NF=IPADI(LE(I,K))
      F(NF,1)=F(NF,1)+SAF1(KP,I)
      F(NF,2)=F(NF,2)+SAF2(KP,I)
      F(NF,3)=F(NF,3)+SAF3(KP,I)
 8017 CONTINUE

 8001 CONTINUE

C     WRITE(6,*)' ********** FIN YCVI 3D *****************'

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