Télécharger zcvit.eso

Retour à la liste

zcvit
C ZCVIT     SOURCE    GOUNAND   23/07/31    21:15:07     11713          
      SUBROUTINE ZCVIT
     &     (HR,RPG,DRR,LE,NEL,K0,NPTD,IES,NP,IAXI,
     &     IPADL,IKOMP,IKAS,
     &     RO,IKR,
     &     COEFF,IK1,RGE,IKG,NELG,TN,IKT,TREF,IKREF,
     &     AMUT,GN,UN,TK,TE,
     &     F,GK,GE,
     &     GK1,GE1,
     &     DK,DE,
     &     VOLU,COTE,NELZ,IDCEN,IMODEL,
     &     DT,DTT1,DTT2,NUEL,DIAEL,DIAM)
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
C***********************************************************************
C  APPELE PAR NSKE
C  CE SP DISCRETISE LES EQUATIONS DE NAVIER STOKES COUPLEES
C  AUX DEUX EQUATIONS DU MODELE K - EPSILON
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***
C
C                    MODELE K - EPSILON
C
C
C    CONSTANTES :  CNU=0.09D0  C1=1.41D0  C2=1.92D0
C                   SGT=0.9D0  SGK=1.D0   SGE=1.3D0
C
C     P --> TABLEAU P
C
C    LE NU TURBULENT EST CALCULE PAR ELEMENT
C    IL N EST PAS TRANSMIS EN ARGUMENT
C
C    IPRE = 1  ON CONTROLE K ET EPSILON A CHAQUE APS DE TEMPS
C
C
C***********************************************************************
 
      DIMENSION UN(NPTD,IES),GN(NPTD,IES),TK(NPTD),TE(NPTD)
      DIMENSION AMUT(*),RO(*)
      DIMENSION COTE(NELZ,IES),VOLU(NELZ)
 
C NOMBRE DE NOEUDS MAX PAR ELEMENT NPX
      PARAMETER (NPX=9)
 
C LONGUEUR DES REGISTRES VECTORIELS DE LA MACHINE CIBLE
 
      PARAMETER (LRV=64)
C
      DIMENSION IPADL(*),LE(NP,*)
      DIMENSION HR(NEL,NP,IES),RPG(*),DRR(NP,NEL)
      DIMENSION QGGT(8,8),Q1(8,8),Q2(8,8),Q3(8,8)
      SAVE QGGT,Q1,Q2,Q3
 
      REAL*8    F(NPTD,IES)
      REAL*8    GK(NPTD),GK1(NPTD),DK(NPTD)
      REAL*8    GE(NPTD),GE1(NPTD),DE(NPTD)
 
-INC CCVQUA4
-INC CCREEL
*-
 
-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
C* -TABLEAUX ADDITIONNELS POUR LA VECTORISATION   **********************
C
C ---- ILS CONDUISENT A UN SURCOUT DE LRV*(59+25NPX) MOTS. AVEC NPX=9
C      ET LRV=128 ON OBTIENT 37KMOTS (0.3 MEGAOCTETS) ET CECI QUELLE
C      SOIT LA NOMBRE TOTAL D'ELEMENTS DU MAILLAGE
C ------------------------------------------------------------------
      DIMENSION AIRE(LRV),COEF(LRV),AL(LRV),AH(LRV),AP(LRV)
      DIMENSION XMB(LRV),XMH (LRV)
      DIMENSION CFM(LRV)
C     DIMENSION CF1(LRV),CF2(LRV),CF3(LRV)
C
      DIMENSION DR(LRV,NPX)
      DIMENSION UIX(LRV,NPX),UIY(LRV,NPX),UIZ(LRV,NPX)
      DIMENSION GIX(LRV,NPX),GIY(LRV,NPX),GIZ(LRV,NPX)
      DIMENSION AK(LRV,NPX),AE(LRV,NPX)
C
      DIMENSION UM(LRV),UMI(LRV,3),AKM(LRV),AEM(LRV),ARO(LRV)
      DIMENSION P(LRV),PG(LRV),ANUT(LRV),
     ;     DUDX(LRV),DUDY(LRV),DUDZ(LRV),DVDX(LRV),DVDY(LRV),
     ;     DVDZ(LRV),DWDX(LRV),DWDY(LRV),DWDZ(LRV),
     ;     DTDX(LRV),DTDY(LRV),DTDZ(LRV)
C
      DIMENSION C1P(LRV),C1G(LRV),COEFT(LRV),COEFK(LRV),COEFE(LRV)
C
      DIMENSION BE1(LRV,NPX)
C
C
C ---- TABLEAUX POUR L'OPTION RAPIDE -------------------------------
C      DANS LA TRIPLE BOUCLE SUR K,I,J ON PEUT SORTIR UN MAXIMUM
C      D'OPERATIONS DE LA BOUCLES LA PLUS INTERNE. LES TABLEAUX CI-
C      DESSOUS PERMETTENT DE STOCKER LES VALEURS CALCULEES
C ------------------------------------------------------------------
C***
      DIMENSION COEFF(*),RGE(NELG,IES)
      DIMENSION TN(*),TREF(*)
 
      DIMENSION RGX(LRV),RGY(LRV),RGZ(LRV)
      DIMENSION SAF1(LRV,9),SAF2(LRV,9),SAF3(LRV,9)
      DIMENSION SAFK(LRV,9),SAFE(LRV,9)
      DIMENSION AL2(LRV),AH2(LRV),AP2(LRV)
      DIMENSION GRADGX(LRV,3),GRADGY(LRV,3),GRADGZ(LRV,3)
      DIMENSION GRADK (LRV,3),GRADE (LRV,3)
      DIMENSION UPGX(LRV),UPIGX(LRV,3),UPGY(LRV),UPIGY(LRV,3)
      DIMENSION UPGZ(LRV),UPIGZ(LRV,3)
      DIMENSION UPK (LRV),UPIK (LRV,3),UPE (LRV),UPIE (LRV,3)
      DIMENSION BM(LRV),BPGX(LRV),BPGY(LRV),BPGZ(LRV)
      DIMENSION         BPK (LRV),BPE (LRV)
      DIMENSION GXCXT(LRV),GXCYT(LRV),GXCZT(LRV)
      DIMENSION GXCXY(LRV),GXCXZ(LRV),GXCYZ(LRV)
      DIMENSION GYCXT(LRV),GYCYT(LRV),GYCZT(LRV)
      DIMENSION GYCXY(LRV),GYCXZ(LRV),GYCYZ(LRV)
      DIMENSION GZCXT(LRV),GZCYT(LRV),GZCZT(LRV)
      DIMENSION GZCXY(LRV),GZCXZ(LRV),GZCYZ(LRV)
      DIMENSION GKCXT(LRV),GKCYT(LRV),GKCZT(LRV)
      DIMENSION GKCXY(LRV),GKCXZ(LRV),GKCYZ(LRV)
      DIMENSION GECXT(LRV),GECYT(LRV),GECZT(LRV)
      DIMENSION GECXY(LRV),GECXZ(LRV),GECYZ(LRV)
 
      DIMENSION CONST(2,7)
 
      SAVE IPAS
      DATA IPAS/0/
 
C Constantes du modèle RNG
C      DATA CNU,C1,C2,C3/0.0845D0,1.42D0,1.68D0,0.8D0/
C      DATA SGT,SGK,SGE /0.9D0,0.7179D0,0.7179D0/
C
C Modèle standard
C     DATA CNU,C1,C2,C3,SGT,SGE/0.09D0,1.41D0,1.92D0,0.8D0,0.9D0,1.3D0/
C
C IMODEL=1 modele Standard   IMODEL=2 modele RNG
C
C
C*** Variation des constantes en fonction du modèle
      DATA (CONST(1,j),j=1,4) /0.09D0,1.41D0,1.92D0,0.8D0/
      DATA (CONST(1,j),j=5,7) /0.9D0,1.D0,1.3D0/
      DATA (CONST(2,j),j=1,4) /0.0845D0,1.42D0,1.68D0,0.8D0/
      DATA (CONST(2,j),j=5,7) /0.9D0,0.7179D0,0.7179D0/
C***
C
CDIR$ VECTOR
CDIR$ IVDEP
 
C      WRITE(IOIMP,*)' Debut YCVIT'
C      WRITE(IOIMP,*)' NELZ,ies=',nelz,ies
C      WRITE(IOIMP,1002)cote
C      WRITE(IOIMP,*)' Avt precaution ies=',ies
C
C           INITIALISATIONS DIVERSES
C
C     zxzini = sqrt(xpetit)
C un peu plus grand sinon on underflow dans les calculs de normes
       zxzini = xpetit**0.33D0
      DTM1=0.D0
      IPRE=0
 
C     IF(IPAS.EQ.0)CALL CALHRG(QGGT,IES)
      IF(IPAS.EQ.0)CALL CALHRH(QGGT,Q1,Q2,Q3,IES)
      NK=K0
C                           ********
C                           *  2D  *
C                           ********
C
C      IMODEL=3
C Suite de la definition des constantes
C      IF(IMODEL.EQ.3) THEN
C       IMO=2
C      ELSE
C       IMO=1
C      ENDIF
C
C  (IMODEL=1 modèle Standard ; IMODEL=2 modèle RNG)
      CNU=CONST(IMODEL,1)
      C1=CONST(IMODEL,2)
      C2=CONST(IMODEL,3)
      C3=CONST(IMODEL,4)
      SGT=CONST(IMODEL,5)
      SGK=CONST(IMODEL,6)
      SGE=CONST(IMODEL,7)
C
      IF(IES.EQ.3)GO TO 10
C
C           PRECAUTIONS
C
C      WRITE(IOIMP,*)' juste Avt precaution K0=',K0
 
      IF(K0.GT.0)GO TO 52
      IF(IPRE.EQ.0)GO TO 52
      DT3=0.D0
      ECHLM=DIAM
C      WRITE(IOIMP,*)' echlm=',echlm,' nptd=',nptd
 
      DO 51 I=1,NPTD
C.      IF(DK(I).LT.XPETIT)GOTO 51
C.      IF(DE(I).LT.XPETIT)GOTO 51
         UMP=UN(I,1)*UN(I,1)+UN(I,2)*UN(I,2)
         UMP=SQRT(UMP)
C         AKMAX=4.D0*UMP*UMP+zxzini
         TK(I)=ABS(TK(I))+zxzini
         TE(I)=ABS(TE(I))+zxzini
         ECHL=(TK(I)+zxzini)**1.5D0/(TE(I)+zxzini)
C     IF(TK(I).LT.AKMAX.AND.ECHL.LT.ECHLM)GO TO 51
         IF(ECHL.LT.ECHLM)GO TO 51
C         REM=UMP*DIAM/COEFF(1)
         DT3=DT3+1.D0
C         TEMAX=AKMAX**1.5D0/ECHLM
CP    IF(TK(I).GT.AKMAX)TK(I)=((TK(I)-AKMAX)*AKMAX/(TK(I)+AKMAX))+AKMAX
         IF(ECHL.GT.ECHLM)TE(I)=TK(I)**1.5D0/
     &        (((ECHL-ECHLM)*ECHLM/(ECHL+ECHLM))+ECHLM)
C     IF(REM.LT.1.D0)TE(I)=CNU*TK(I)**2.D0*10.D0/COEFF(1)
 51   CONTINUE
 52   CONTINUE
C      WRITE(IOIMP,*)' Apr precaution '
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 
C                           DIFFERENCES TRIANGLE / QUADRANGLE
      QUA4=0.D0
      IF(NP.EQ.4)QUA4=1.D0
 
C
C K  EST LE NUMERO DE L'ELEMENT
C KP PERMET DE SE SITUER A L'INTERIEUR DES TABLEAUX DE LRV ELEMENTS
C --- CALCUL DU NOMBRE DE PAQUETS DE LRV ELEMENTS ---
      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 D'ELEMENTS ***
C
      DO 70001 KPACK=KPACKD,KPACKF
C
C --- POUR CHAQUE PAQUET DE LRV ELEMENTS
         KDEB=1+(KPACK-1)*LRV
         KFIN=MIN(NEL,KDEB+LRV-1)
C      WRITE(IOIMP,*)' KPACK=',KPACK
C
         DO 70002 K=KDEB,KFIN
            KP=K-KDEB+1
            NK=K+K0
            K1=1+(1-IK1)*(NK-1)
            AIRE(KP)=VOLU(NK)
            COEF(KP)=COEFF(K1)
            AL(KP)=COTE(NK,1)
            AH(KP)=COTE(NK,2)
            AL2(KP)=1.D0/AL(KP)/AL(KP)
            AH2(KP)=1.D0/AH(KP)/AH(KP)
            XMB (KP)=(AL(KP)+AH(KP))*0.5D0
70002    CONTINUE
C      WRITE(IOIMP,*)' Apr bcl 70002,,,,,,,,,,,,,,,,,,,,,,,, '
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Initialisation des UMI,UMIT avant accumulation
         DO 7005 K=KDEB,KFIN
            KP=K-KDEB+1
            UMI(KP,1)=zxzini
            UMI(KP,2)=zxzini
            GRADGX(KP,1)=zxzini
            GRADGX(KP,2)=zxzini
            GRADGY(KP,1)=zxzini
            GRADGY(KP,2)=zxzini
            AKM(KP)=zxzini
            AEM(KP)=zxzini
            UPIK(KP,1)=zxzini
            UPIK(KP,2)=zxzini
            UPIE(KP,1)=zxzini
            UPIE(KP,2)=zxzini
            GRADK(KP,1)=zxzini
            GRADK(KP,2)=zxzini
            GRADE(KP,1)=zxzini
            GRADE(KP,2)=zxzini
            DUDX(KP)=0.D0
            DUDY(KP)=0.D0
            DVDX(KP)=0.D0
            DVDY(KP)=0.D0
            DTDX(KP)=0.D0
            DTDY(KP)=0.D0
            C1P(KP)=C1
            C1G(KP)=0.D0
            PG(KP)=0.D0
 7005    CONTINUE
 
         DO 7006 I=1,NP
            DO 6006 K=KDEB,KFIN
               KP=K-KDEB+1
               NF=IPADL(LE(I,K))
               DR(KP,I)=DRR(I,K)
               UMI(KP,1)=UMI(KP,1)+UN(NF,1)*DR(KP,I)
               UMI(KP,2)=UMI(KP,2)+UN(NF,2)*DR(KP,I)
               UIX(KP,I)=UN(NF,1)
               UIY(KP,I)=UN(NF,2)
               GIX(KP,I)=GN(NF,1)
               GIY(KP,I)=GN(NF,2)
               AK(KP,I)=TK(NF)
               AE(KP,I)=TE(NF)
               GRADGX(KP,1)=GRADGX(KP,1)+HR(K,I,1)*GIX(KP,I)
               GRADGX(KP,2)=GRADGX(KP,2)+HR(K,I,2)*GIX(KP,I)
               GRADGY(KP,1)=GRADGY(KP,1)+HR(K,I,1)*GIY(KP,I)
               GRADGY(KP,2)=GRADGY(KP,2)+HR(K,I,2)*GIY(KP,I)
               AKM(KP)=AKM(KP)+TK(NF)*DR(KP,I)
               AEM(KP)=AEM(KP)+TE(NF)*DR(KP,I)
               GRADK(KP,1)=GRADK(KP,1)+HR(K,I,1)*AK(KP,I)
               GRADK(KP,2)=GRADK(KP,2)+HR(K,I,2)*AK(KP,I)
               GRADE(KP,1)=GRADE(KP,1)+HR(K,I,1)*AE(KP,I)
               GRADE(KP,2)=GRADE(KP,2)+HR(K,I,2)*AE(KP,I)
C *** CALCUL DU TERME DE PRODUCTION ****************************
               DUDX(KP) = DUDX(KP)+HR(K,I,1)*UN(NF,1)
               DUDY(KP) = DUDY(KP)+HR(K,I,2)*UN(NF,1)
               DVDX(KP) = DVDX(KP)+HR(K,I,1)*UN(NF,2)
               DVDY(KP) = DVDY(KP)+HR(K,I,2)*UN(NF,2)
 6006       CONTINUE
 7006    CONTINUE
 
C     WRITE(IOIMP,*)'Initialisation SAF1,SAF2,SBF '
C
C Initialisation des variables d'accumulation SAF1,SAF2,SBF
C
         IF(IKOMP.EQ.0)THEN
 
            DO 70060 K=KDEB,KFIN
               KP=K-KDEB+1
               ARO(KP)=1.D0
70060       CONTINUE
 
            IF(IKAS.EQ.2)THEN
               DO 70061 I=1,NP
                  DO 60061 K=KDEB,KFIN
                     KP=K-KDEB+1
                     SAF1(KP,I)=0.D0
                     SAF2(KP,I)=0.D0
                     SAFK(KP,I)=0.D0
                     SAFE(KP,I)=0.D0
60061             CONTINUE
70061          CONTINUE
            ELSEIF(IKAS.EQ.3)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
               DO 70062 I=1,NP
                  DO 60062 K=KDEB,KFIN
                     KP=K-KDEB+1
                     SAF1(KP,I)=(-RGX(KP))*DR(KP,I)
                     SAF2(KP,I)=(-RGY(KP))*DR(KP,I)
                     SAFK(KP,I)=0.D0
                     SAFE(KP,I)=0.D0
                     NKG=1+(1-IKG)*(NK-1)
                     DTDX(KP)=DTDX(KP)+HR(K,I,1)*RGE(NKG,1)
                     DTDY(KP)=DTDY(KP)+HR(K,I,2)*RGE(NKG,2)
60062             CONTINUE
70062          CONTINUE
            ELSEIF(IKAS.EQ.5)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
               DO 70063 I=1,NP
                  DO 60063 K=KDEB,KFIN
                     KP=K-KDEB+1
                     NF=1+(1-IKT)*(IPADL(LE(I,K))-1)
                     NFR=1+(1-IKREF)*(IPADL(LE(I,K))-1)
                     SAF1(KP,I)=(-RGX(KP)*(TREF(NFR)-TN(NF)))*DR(KP,I)
                     SAF2(KP,I)=(-RGY(KP)*(TREF(NFR)-TN(NF)))*DR(KP,I)
                     SAFK(KP,I)=0.D0
                     SAFE(KP,I)=0.D0
                     DTDX(KP)=DTDX(KP)+HR(K,I,1)*TN(NF)
                     DTDY(KP)=DTDY(KP)+HR(K,I,2)*TN(NF)
60063             CONTINUE
70063          CONTINUE
            ENDIF
 
         ELSEIF(IKOMP.EQ.1)THEN
 
            DO 70020 K=KDEB,KFIN
               KP=K-KDEB+1
               NK=K+K0
               NKR=1+(1-IKR)*(NK-1)
               ARO(KP)=RO(NKR)
70020       CONTINUE
 
            IF(IKAS.EQ.4)THEN
               DO 70064 I=1,NP
                  DO 60064 K=KDEB,KFIN
                     KP=K-KDEB+1
                     SAF1(KP,I)=0.D0
                     SAF2(KP,I)=0.D0
                     SAFK(KP,I)=0.D0
                     SAFE(KP,I)=0.D0
60064             CONTINUE
70064          CONTINUE
            ELSEIF(IKAS.EQ.5)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
               DO 70065 I=1,NP
                  DO 60065 K=KDEB,KFIN
                     KP=K-KDEB+1
                     SAF1(KP,I)=-RGX(KP)*DR(KP,I)
                     SAF2(KP,I)=-RGY(KP)*DR(KP,I)
                     SAFK(KP,I)=0.D0
                     SAFE(KP,I)=0.D0
60065             CONTINUE
70065          CONTINUE
            ENDIF
         ENDIF
 
         DO 7007 K=KDEB,KFIN
            KP=K-KDEB+1
            NK=K+K0
            AKM(KP)=ABS(AKM(KP))/AIRE(KP)+zxzini
            AEM(KP)=ABS(AEM(KP))/AIRE(KP)+zxzini
            UMI(KP,1)=UMI(KP,1)/AIRE(KP)
            UMI(KP,2)=UMI(KP,2)/AIRE(KP)
            UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2)
            UM(KP)=SQRT(UM(KP))
            AX=GRADGX(KP,1)*GRADGX(KP,1)+GRADGX(KP,2)*GRADGX(KP,2)
            AX=SQRT(AX)+zxzini
            GRADGX(KP,1)=GRADGX(KP,1)/AX
            GRADGX(KP,2)=GRADGX(KP,2)/AX
            AY=GRADGY(KP,1)*GRADGY(KP,1)+GRADGY(KP,2)*GRADGY(KP,2)
            AY=SQRT(AY)+zxzini
            GRADGY(KP,1)=GRADGY(KP,1)/AY
            GRADGY(KP,2)=GRADGY(KP,2)/AY
            UPGX(KP)=GRADGX(KP,1)*UMI(KP,1)+GRADGX(KP,2)*UMI(KP,2)
            UPIGX(KP,1)=UPGX(KP)*GRADGX(KP,1)
            UPIGX(KP,2)=UPGX(KP)*GRADGX(KP,2)
            UPGY(KP)=GRADGY(KP,1)*UMI(KP,1)+GRADGY(KP,2)*UMI(KP,2)
            UPIGY(KP,1)=UPGY(KP)*GRADGY(KP,1)
            UPIGY(KP,2)=UPGY(KP)*GRADGY(KP,2)
 
            AX=GRADK(KP,1)*GRADK(KP,1)+GRADK(KP,2)*GRADK(KP,2)
            AX=SQRT(AX)+zxzini
            GRADK(KP,1)=GRADK(KP,1)/AX
            GRADK(KP,2)=GRADK(KP,2)/AX
            UPK(KP)=GRADK(KP,1)*UMI(KP,1)+GRADK(KP,2)*UMI(KP,2)
            UPIK(KP,1)=UPK(KP)*GRADK(KP,1)
            UPIK(KP,2)=UPK(KP)*GRADK(KP,2)
 
            AX=GRADE(KP,1)*GRADE(KP,1)+GRADE(KP,2)*GRADE(KP,2)
            AX=SQRT(AX)+zxzini
            GRADE(KP,1)=GRADE(KP,1)/AX
            GRADE(KP,2)=GRADE(KP,2)/AX
            UPE(KP)=GRADE(KP,1)*UMI(KP,1)+GRADE(KP,2)*UMI(KP,2)
            UPIE(KP,1)=UPE(KP)*GRADE(KP,1)
            UPIE(KP,2)=UPE(KP)*GRADE(KP,2)
 
            ANUT(KP)=CNU*AKM(KP)*AKM(KP)/AEM(KP)
            AMUT(NK)=ARO(KP)*ANUT(KP)
            COEFT(KP)=COEF(KP)+AMUT(NK)
            COEFK(KP)=COEF(KP)+ANUT(KP)/SGK
            COEFE(KP)=COEF(KP)+ANUT(KP)/SGE
C --- 2. CALCUL DU TERME DE PRODUCTION PROPREMENT DIT ----------
            PTEMP=( (DUDY(KP)+DVDX(KP))**2
     $            + 2.D0*(DUDX(KP)**2+DVDY(KP)**2) )
            P(KP)=PTEMP
 7007    CONTINUE
C------------ANCIENNE VERSION--------------------
C      IF(IAXI.NE.0) THEN
C        DO 70012 K=KDEB,KFIN
C        KP=K+1-KDEB
C        P(KP)=P(KP)+2.D0*(UMI(KP,1)*UMI(KP,1))/(RPG(K)*RPG(K))
C70012   CONTINUE
C      END IF
C
C      IF(IKAS.EQ.5.AND.IKOMP.EQ.0)THEN
C        IF(IMODEL.EQ.1)THEN
C        DO 70113 K=KDEB,KFIN
C        KP=K+1-KDEB
C        PG(KP)=(RGX(KP)*DTDX(KP)+RGY(KP)*DTDY(KP))/SGT
C        IF(PG(KP).GT.0.D0)C1G(KP)=C1
C70113   CONTINUE
C        ELSEIF(IMODEL.EQ.2)THEN
C        DO 70213 K=KDEB,KFIN
C        KP=K+1-KDEB
C        PG(KP)=(RGX(KP)*DTDX(KP)+RGY(KP)*DTDY(KP))/SGT
C          RF    =-PG(KP)/(P(KP)+PG(KP)+zxzini)
C          IF(PG(KP).GT.0.D0)RF=0.D0
C          C1G(KP)=C1*(1.D0+C3*RF)
C          C1P(KP)=C1*(1.D0+C3*RF)
C70213   CONTINUE
C        ENDIF
C      ELSEIF(IKAS.EQ.2.AND.IKOMP.EQ.0)THEN
C        IF(IMODEL.EQ.1)THEN
C        DO 70114 K=KDEB,KFIN
C        KP=K+1-KDEB
C        PG(KP)=(DTDX(KP)+DTDY(KP))/SGT
C        IF(PG(KP).GT.0.D0)C1G(KP)=C1
C70114   CONTINUE
C        ELSEIF(IMODEL.EQ.2)THEN
C        DO 70214 K=KDEB,KFIN
C        KP=K+1-KDEB
C        PG(KP)=(DTDX(KP)+DTDY(KP))/SGT
C          RF    =-PG(KP)/(P(KP)+PG(KP)+zxzini)
C          IF(PG(KP).GT.0.D0)RF=0.D0
C          C1G(KP)=C1*(1.D0+C3*RF    )
C          C1P(KP)=C1*(1.D0+C3*RF    )
C70214   CONTINUE
C        END IF
C      END IF
C-------------VERSION POUR RNG---------------
         IF(IKAS.EQ.5.AND.IKOMP.EQ.0)THEN
            DO 70213 K=KDEB,KFIN
               KP=K+1-KDEB
               PG(KP)=(RGX(KP)*DTDX(KP)+RGY(KP)*DTDY(KP))/SGT
               RF    =-PG(KP)/(P(KP)+PG(KP)+zxzini)
               IF(PG(KP).GT.0.D0)RF=0.D0
               C1G(KP)=C1*(1.D0+C3*RF)
               C1P(KP)=C1*(1.D0+C3*RF)
70213       CONTINUE
         ELSEIF(IKAS.EQ.2.AND.IKOMP.EQ.0)THEN
            DO 70214 K=KDEB,KFIN
               KP=K+1-KDEB
               PG(KP)=(DTDX(KP)+DTDY(KP))/SGT
               RF    =-PG(KP)/(P(KP)+PG(KP)+zxzini)
               IF(PG(KP).GT.0.D0)RF=0.D0
               C1G(KP)=C1*(1.D0+C3*RF    )
               C1P(KP)=C1*(1.D0+C3*RF    )
70214       CONTINUE
         END IF
 
         IF(IMODEL.EQ.2) THEN
            ETA0=4.377D0
            BETA0=0.012D0
            ETA=(P(KP)**0.5D0)*AKM(KP)/AEM(KP)
            RETA= ETA*(1.D0-(ETA/ETA0))/(1.D0+BETA0*(ETA**3.D0))
            C1P(KP)=C1P(KP)-RETA
         END IF
C --------------------------------------------------------------
C *** FIN DU CALCUL DU TERME DE PRODUCTION *********************
 
         DO 70074 K=KDEB,KFIN
            KP=K-KDEB+1
            BMX=UMI(KP,1)/AL(KP)
            BMY=UMI(KP,2)/AH(KP)
            BM(KP)=BMX*BMX+BMY*BMY
            BM(KP)=SQRT(BM(KP))+zxzini
            BPGXX=UPIGX(KP,1)/AL(KP)
            BPGXY=UPIGX(KP,2)/AH(KP)
            BPGX(KP)=BPGXX*BPGXX+BPGXY*BPGXY
            BPGX(KP)=SQRT(BPGX(KP))+zxzini
            BPGYX=UPIGY(KP,1)/AL(KP)
            BPGYY=UPIGY(KP,2)/AH(KP)
            BPGY(KP)=BPGYX*BPGYX+BPGYY*BPGYY
            BPGY(KP)=SQRT(BPGY(KP))+zxzini
            BPKX=UPIK(KP,1)/AL(KP)
            BPKY=UPIK(KP,2)/AH(KP)
            BPK(KP)=BPKX*BPKX+BPKY*BPKY
            BPK(KP)=SQRT(BPK(KP))+zxzini
            BPEX=UPIE(KP,1)/AL(KP)
            BPEY=UPIE(KP,2)/AH(KP)
            BPE(KP)=BPEX*BPEX+BPEY*BPEY
            BPE(KP)=SQRT(BPE(KP))+zxzini
70074    CONTINUE
C
C     AVANT DECENTREMENT
C
C ALFA,AKSI,CCU utilisées seulement dans la boucle 7008
         IF(IDCEN.EQ.1)THEN
            DO 70081 K=KDEB,KFIN
               KP=K-KDEB+1
               GXCXT(KP)=0.D0
               GXCYT(KP)=0.D0
               GXCXY(KP)=0.D0
               GYCXT(KP)=0.D0
               GYCYT(KP)=0.D0
               GYCXY(KP)=0.D0
               GKCXT(KP)=0.D0
               GKCYT(KP)=0.D0
               GKCXY(KP)=0.D0
               GECXT(KP)=0.D0
               GECYT(KP)=0.D0
               GECXY(KP)=0.D0
70081       CONTINUE
 
         ELSEIF(IDCEN.EQ.2)THEN
            DO 7008 K=KDEB,KFIN
               KP=K-KDEB+1
               HMK=2.D0*UM(KP)/BM(KP)
               XMB(KP)=HMK
               ALFA=UM(KP)*HMK/COEFT(KP)/2.D0
               AKSI=ALFA/3.D0
               IF(ALFA.GT.3.D0)AKSI=1.D0
               CCT=AKSI/BM(KP)
 
               HPK=2.D0*UPGX(KP)/BPGX(KP)
               ALFA=UPGX(KP)*HPK/COEFT(KP)/2.D0
               AKSI=ALFA/3.D0
               IF(ALFA.GT.3.D0)AKSI=1.D0
               CCP=AKSI/BPGX(KP)
               CPT=CCP-CCT
               CC2X=0.D0
               IF(CPT.GE.0.D0)CC2X=CPT
 
               HPK=2.D0*UPGY(KP)/BPGY(KP)
               ALFA=UPGY(KP)*HPK/COEFT(KP)/2.D0
               AKSI=ALFA/3.D0
               IF(ALFA.GT.3.D0)AKSI=1.D0
               CCP=AKSI/BPGY(KP)
               CPT=CCP-CCT
               CC2Y=0.D0
               IF(CPT.GE.0.D0)CC2Y=CPT
 
               HPK=2.D0*UPK(KP)/BPK(KP)
               ALFA=UPK(KP)*HPK/COEFK(KP)/2.D0
               AKSI=ALFA/3.D0
               IF(ALFA.GT.3.D0)AKSI=1.D0
               CCP=AKSI/BPK(KP)
               CPT=CCP-CCT
               CC2K=0.D0
               IF(CPT.GE.0.D0)CC2K=CPT
 
               HPK=2.D0*UPE(KP)/BPE(KP)
               ALFA=UPE(KP)*HPK/COEFE(KP)/2.D0
               AKSI=ALFA/3.D0
               IF(ALFA.GT.3.D0)AKSI=1.D0
               CCP=AKSI/BPE(KP)
               CPT=CCP-CCT
               CC2E=0.D0
               IF(CPT.GE.0.D0)CC2E=CPT
 
               GXCXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT
     $                   +UPIGX(KP,1)*UPIGX(KP,1)*CC2X)
               GXCYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT
     $                   +UPIGX(KP,2)*UPIGX(KP,2)*CC2X)
               GXCXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT
     $                   +UPIGX(KP,1)*UPIGX(KP,2)*CC2X)
 
               GYCXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT
     $                   +UPIGY(KP,1)*UPIGY(KP,1)*CC2Y)
               GYCYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT
     $                   +UPIGY(KP,2)*UPIGY(KP,2)*CC2Y)
               GYCXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT
     $                   +UPIGY(KP,1)*UPIGY(KP,2)*CC2Y)
 
               GKCXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT
     $                   +UPIK(KP,1)*UPIK(KP,1)*CC2K)
               GKCYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT
     $                   +UPIK(KP,2)*UPIK(KP,2)*CC2K)
               GKCXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT
     $                   +UPIK(KP,1)*UPIK(KP,2)*CC2K)
 
               GECXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT
     $                   +UPIE(KP,1)*UPIE(KP,1)*CC2E)
               GECYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT
     $                   +UPIE(KP,2)*UPIE(KP,2)*CC2E)
               GECXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT
     $                   +UPIE(KP,1)*UPIE(KP,2)*CC2E)
 
 7008       CONTINUE
 
         ELSEIF(IDCEN.EQ.3)THEN
            DO 7018 K=KDEB,KFIN
               KP=K-KDEB+1
               HMK=2.D0*UM(KP)/BM(KP)
               XMB(KP)=HMK
               ALFA=UM(KP)*HMK/COEFT(KP)/2.D0
               AKSI=ALFA/3.D0
               IF(ALFA.GT.3.D0)AKSI=1.D0
               CCT=AKSI/BM(KP)
 
               GXCXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT
               GXCYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT
               GXCXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT
 
               GYCXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT
               GYCYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT
               GYCXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT
 
               GKCXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT
               GKCYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT
               GKCXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT
 
               GECXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT
               GECYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT
               GECXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT
 
 7018       CONTINUE
 
         ELSEIF(IDCEN.EQ.4)THEN
            DT19=DTM1*0.5D0
            DO 7009 K=KDEB,KFIN
               KP=K-KDEB+1
               XMB(KP)=(AL(KP)+AH(KP))/2.D0
               GXCXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19
               GXCYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19
               GXCXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19
               GYCXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19
               GYCYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19
               GYCXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19
               GKCXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19
               GKCYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19
               GKCXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19
               GECXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19
               GECYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19
               GECXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19
 7009       CONTINUE
 
         ENDIF
 
C     WRITE(IOIMP,*)' Avt Calcul DT'
C
C     AVANT CALCUL DT
C
         DO 7010 K=KDEB,KFIN
            KP=K-KDEB+1
            DT0=DT
            DT1X=2.D0/
     &           (UMI(KP,1)*UMI(KP,1)/(GXCXT(KP)+COEFT(KP))
     &           +UMI(KP,2)*UMI(KP,2)/(GXCYT(KP)+COEFT(KP)))
 
            DT1Y=2.D0/
     &           (UMI(KP,1)*UMI(KP,1)/(GYCXT(KP)+COEFT(KP))
     &           +UMI(KP,2)*UMI(KP,2)/(GYCYT(KP)+COEFT(KP)))
            DT2X=0.5D0/
     &           ((GXCXT(KP)+COEFT(KP))*AL2(KP)
     &           +(GXCYT(KP)+COEFT(KP))*AH2(KP))
            DT2Y=0.5D0/
     &           ((GYCXT(KP)+COEFT(KP))*AL2(KP)
     &           +(GYCYT(KP)+COEFT(KP))*AH2(KP))
            IF(DT1X.LT.DT)DT=DT1X
            IF(DT1Y.LT.DT)DT=DT1Y
            IF(DT2X.LT.DT)DT=DT2X
            IF(DT2Y.LT.DT)DT=DT2Y
C     IF(DT3.LT.DT)DT=DT3
            IF(DT.NE.DT0) THEN
               DTT1=MIN(DT1X,DT1Y)
               DTT2=MIN(DT2X,DT2Y)
C        DTT3=DT3
               DIAEL=XMB(KP)
               NUEL=K
            END IF
 7010    CONTINUE
 
C     WRITE(IOIMP,*)' le coeur du calcul '
C Le coeur du calcul ...
 
         IF(IKOMP.EQ.0)THEN
 
            DO 7014 I=1,NP
               DO 70141 J=1,NP
                  DO 70142 K=KDEB,KFIN
                     KP=K-KDEB+1
                     ZVGGX=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*GXCXT(KP)
     &               +    HR(K,I,1)*HR(K,J,2)*GXCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,1)*GXCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*GXCYT(KP)
     &               +    (GXCXT(KP)*AL2(KP)+GXCYT(KP)*AH2(KP))
     $                    /12.D0*VGGT(J,I)*QUA4  )
 
                     ZVGGY=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*GYCXT(KP)
     &               +    HR(K,I,1)*HR(K,J,2)*GYCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,1)*GYCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*GYCYT(KP)
     &               +    (GYCXT(KP)*AL2(KP)+GYCYT(KP)*AH2(KP))/12
     $                    .D0*VGGT(J,I)*QUA4  )
 
                     ZVGK =AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*GKCXT(KP)
     &               +    HR(K,I,1)*HR(K,J,2)*GKCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,1)*GKCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*GKCYT(KP)
     &               +    (GKCXT(KP)*AL2(KP)+GKCYT(KP)*AH2(KP))/12
     $                    .D0*VGGT(J,I)*QUA4  )
 
                     ZVGE =AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*GECXT(KP)
     &               +    HR(K,I,1)*HR(K,J,2)*GECXY(KP)
     &               +    HR(K,I,2)*HR(K,J,1)*GECXY(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*GECYT(KP)
     &               +    (GECXT(KP)*AL2(KP)+GECYT(KP)*AH2(KP))/12
     $                    .D0*VGGT(J,I)*QUA4  )
 
                     ZVGD=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*COEFT(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*COEFT(KP)
     &               +    (COEFT(KP)*AL2(KP)+COEFT(KP)*AH2(KP))/12
     $                    .D0*VGGT(J,I)*QUA4  )
 
                     ZVGDK=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*COEFK(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*COEFK(KP)
     &               +    (COEFK(KP)*AL2(KP)+COEFK(KP)*AH2(KP))/12
     $                    .D0*VGGT(J,I)*QUA4  )
 
                     ZVGDE=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*COEFE(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*COEFE(KP)
     &               +    (COEFE(KP)*AL2(KP)+COEFE(KP)*AH2(KP))/12
     $                    .D0*VGGT(J,I)*QUA4  )
 
Cnon  V2=(UIX(KP,I)*HR(K,J,1)+UIY(KP,I)*HR(K,J,2))*DR(KP,I)
 
                     V2=UMI(KP,1)*DR(KP,I)*HR(K,J,1)+UMI(KP,2)*DR(KP,I)
     $                    *HR(K,J,2)
 
                     SAF1(KP,I)=SAF1(KP,I)+(V2+ZVGGX)*GIX(KP,J)+ZVGD
     $                    *UIX(KP,J)
                     SAF2(KP,I)=SAF2(KP,I)+(V2+ZVGGY)*GIY(KP,J)+ZVGD
     $                    *UIY(KP,J)
 
                     SAFK(KP,I)=SAFK(KP,I)+(V2+ZVGK)*AK(KP,J)+ZVGDK
     $                    *AK(KP,J)
                     SAFE(KP,I)=SAFE(KP,I)+(V2+ZVGE)*AE(KP,J)+ZVGDE
     $                    *AE(KP,J)
 
70142             CONTINUE
70141          CONTINUE
 7014       CONTINUE
 
         ELSEIF(IKOMP.EQ.1)THEN
 
            DO 7015 I=1,NP
               DO 70151 J=1,NP
                  DO 70152 K=KDEB,KFIN
                     KP=K-KDEB+1
                     ZVGGX=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*GXCXT(KP)
     &               +    HR(K,I,1)*HR(K,J,2)*GXCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,1)*GXCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*GXCYT(KP)
     &               +    (GXCXT(KP)*AL2(KP)+GXCYT(KP)*AH2(KP))/12
     $                    .D0*VGGT(J,I)*QUA4  )
 
                     ZVGGY=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*GYCXT(KP)
     &               +    HR(K,I,1)*HR(K,J,2)*GYCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,1)*GYCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*GYCYT(KP)
     &               +    (GYCXT(KP)*AL2(KP)+GYCYT(KP)*AH2(KP))/12
     $                    .D0*VGGT(J,I)*QUA4  )
 
                     ZVGK =AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*GKCXT(KP)
     &               +    HR(K,I,1)*HR(K,J,2)*GKCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,1)*GKCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*GKCYT(KP)
     &               +    (GKCXT(KP)*AL2(KP)+GKCYT(KP)*AH2(KP))/12
     $                    .D0*VGGT(J,I)*QUA4  )
 
                     ZVGE =AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*GECXT(KP)
     &               +    HR(K,I,1)*HR(K,J,2)*GECXY(KP)
     &               +    HR(K,I,2)*HR(K,J,1)*GECXY(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*GECYT(KP)
     &               +    (GECXT(KP)*AL2(KP)+GECYT(KP)*AH2(KP))/12
     $                    .D0*VGGT(J,I)*QUA4  )
 
                     ZVGD=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*COEFT(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*COEFT(KP)
     &               +    (COEFT(KP)*AL2(KP)+COEFT(KP)*AH2(KP))/12
     $                    .D0*VGGT(J,I)*QUA4  )
 
                     ZVGDK=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*COEFK(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*COEFK(KP)
     &               +    (COEFK(KP)*AL2(KP)+COEFK(KP)*AH2(KP))/12
     $                    .D0*VGGT(J,I)*QUA4  )
 
                     ZVGDE=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*COEFE(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*COEFE(KP)
     &               +    (COEFE(KP)*AL2(KP)+COEFE(KP)*AH2(KP))/12
     $                    .D0*VGGT(J,I)*QUA4  )
 
                     COEFT3=COEFT(KP)/3.D0
                     ZVGUU=AIRE(KP)*( HR(K,I,1)*HR(K,J,1)
     &                 +    AL2(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)*( HR(K,I,2)*HR(K,J,2)
     &                +    AH2(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))*DR(KP
     $                    ,I)
 
                     SAF1(KP,I)=SAF1(KP,I)+(V2+ZVGGX)*GIX(KP,J)+ZVGD
     $                    *UIX(KP,J)+ ZVGUU + ZVGUV
                     SAF2(KP,I)=SAF2(KP,I)+(V2+ZVGGY)*GIY(KP,J)+ZVGD
     $                    *UIY(KP,J)+ ZVGVU + ZVGVV
 
                     SAFK(KP,I)=SAFK(KP,I)+(V2+ZVGK)*AK(KP,J)+ZVGDK
     $                    *AK(KP,J)
                     SAFE(KP,I)=SAFE(KP,I)+(V2+ZVGE)*AE(KP,J)+ZVGDE
     $                    *AE(KP,J)
 
70152             CONTINUE
70151          CONTINUE
 7015       CONTINUE
 
         ENDIF
 
         IF(IAXI.NE.0) THEN
            DO 7016 I=1,NP
               DO 70161 J=1,NP
                  DO 70162 K=KDEB,KFIN
                     KP=K-KDEB+1
C     IP=LE(I,K)
C     XG=XCOOR((IP-1)*(IDIM+1)    +1)+zxzini
C        R2=1.D0/XG/XG*DR(KP,I)
 
                     R2=1.D0/RPG(K)/RPG(K)*DR(KP,I)
                     AX=(GXCXT(KP)+GXCYT(KP))*0.5D0
                     SAF1(KP,I)=SAF1(KP,I)+R2*(COEFT(KP)*UIX(KP,J)+AX
     $                    *GIX(KP,J))
 
70162             CONTINUE
70161          CONTINUE
 7016       CONTINUE
 
            IF(IKOMP.EQ.1)THEN
 
               DO 7118 I=1,NP
                  DO 71181 K=KDEB,KFIN
                     KP=K-KDEB+1
 
                     R1=1.D0/RPG(K)*DR(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)
 
71181             CONTINUE
 7118          CONTINUE
 
            ENDIF
 
         ENDIF
C
C --------------------------------------------------------------
C
C
         DO 7017 I=1,NP
            DO 70171 K=KDEB,KFIN
               KP=K-KDEB+1
               TSK       =(ANUT(KP)*(P(KP)+PG(KP))-AEM(KP))*DR(KP,I)
               SAFK(KP,I)=SAFK(KP,I)-TSK
 
               TSE       =(CNU*(C1P(KP)*P(KP)+C1G(KP)*PG(KP))*AKM(KP))
     $              *DR(KP,I)
               SAFE(KP,I)=SAFE(KP,I)-TSE
 
               TSE1      =C2*AEM(KP)/AKM(KP)*DR(KP,I)
               BE1(KP,I) =TSE1
70171       CONTINUE
 7017    CONTINUE
 
C
C
C --- ECRITURE DANS LES TABLEAUX F(1,*),F(2,*),G,GK,GE,GE1 (SCALAIRE)
         DO 70019 I=1,NP
            DO 60019 K=KDEB,KFIN
               KP=K+1-KDEB
               NF=IPADL(LE(I,K))
               F(NF,1)=F(NF,1)+SAF1(KP,I)
               F(NF,2)=F(NF,2)+SAF2(KP,I)
               GK(NF) =GK(NF) +SAFK(KP,I)
               GE1(NF)=GE1(NF)+BE1(KP,I)
               GE(NF) =GE (NF)+SAFE(KP,I)
60019       CONTINUE
70019    CONTINUE
70001 CONTINUE
C
C *** FIN DE LA BOUCLE SUR LES PAQUETS D'ELEMENTS ***
C
C     WRITE(IOIMP,*)' FIN YCVIT GK='
C     WRITE(IOIMP,1002)(GK(i),i=1,nptd)
 
      IPAS=1
      RETURN
 
C                           ********
C                           *  3D  *
C                           ********
 
 10   CONTINUE
C
C           PRECAUTIONS
C
      IF(K0.GT.0)GO TO 62
      IF(IPRE.EQ.0)GO TO 62
      DT3=0.D0
      ECHLM=DIAM
      DO 61 I=1,NPTD
         IF(DK(I).LT.XPETIT)GOTO 61
         IF(DE(I).LT.XPETIT)GOTO 61
         UMP=UN(I,1)*UN(I,1)+UN(I,2)*UN(I,2)+UN(I,3)*UN(I,3)
         UMP=SQRT(UMP)
C         AKMAX=4.D0*UMP*UMP+zxzini
         TK(I)=ABS(TK(I))+zxzini
         TE(I)=ABS(TE(I))+zxzini
         ECHL=(TK(I)+zxzini)**1.5D0/(TE(I)+zxzini)
C     IF(TK(I).LT.AKMAX.AND.ECHL.LT.ECHLM)GO TO 61
         IF(ECHL.LT.ECHLM)GO TO 61
C         REM=UMP*DIAM/COEFF(1)
         DT3=DT3+1.D0
C         TEMAX=AKMAX**1.5D0/ECHLM
CP    IF(TK(I).GT.AKMAX)TK(I)=((TK(I)-AKMAX)*AKMAX/(TK(I)+AKMAX))+AKMAX
         IF(ECHL.GT.ECHLM)TE(I)=TK(I)**1.5D0/
     &        (((ECHL-ECHLM)*ECHLM/(ECHL+ECHLM))+ECHLM)
C     IF(REM.LT.1.D0)TE(I)=CNU*TK(I)**2.D0*10.D0/COEFF(1)
 61   CONTINUE
 62   CONTINUE
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 
C                           DIFFERENCES HEXAHEDRES / AUTRES ELEMENTS
      CUB8=0.D0
      IF(NP.EQ.8)CUB8=1.D0
 
C
C K  EST LE NUMERO DE L'ELEMENT
C KP PERMET DE SE SITUER A L'INTERIEUR DES TABLEAUX DE LRV ELEMENTS
C --- CALCUL DU NOMBRE DE PAQUETS DE LRV ELEMENTS ---
      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 D'ELEMENTS ***
C
      DO 80001 KPACK=KPACKD,KPACKF
C
C --- POUR CHAQUE PAQUET DE LRV ELEMENTS
         KDEB=1+(KPACK-1)*LRV
         KFIN=MIN(NEL,KDEB+LRV-1)
C      WRITE(IOIMP,*)' KPACK=',KPACK
C
         DO 80002 K=KDEB,KFIN
            KP=K-KDEB+1
            NK=K+K0
            K1=1+(1-IK1)*(NK-1)
            AIRE(KP)=VOLU(NK)
            COEF(KP)=COEFF(K1)
            AL(KP)=COTE(NK,1)
            AH(KP)=COTE(NK,2)
            AP(KP)=COTE(NK,3)
            AL2(KP)=1.D0/AL(KP)/AL(KP)
            AH2(KP)=1.D0/AH(KP)/AH(KP)
            AP2(KP)=1.D0/AP(KP)/AP(KP)
            CFM(KP)=AL(KP)*AH(KP)/AP(KP)+AL(KP)*AP(KP)/AH(KP)+
     &              AP(KP)*AH(KP)/AL(KP)
C     CF1(KP)=AL(KP)*AH(KP)/AP(KP)
C     CF2(KP)=AL(KP)*AP(KP)/AH(KP)
C     CF3(KP)=AP(KP)*AH(KP)/AL(KP)
            XMH (KP)=(AL(KP)+AH(KP)+AP(KP))/3.D0
            XMB (KP)=XMH (KP)
80002    CONTINUE
C      WRITE(IOIMP,*)' Apr bcl 80002,,,,,,,,,,,,,,,,,,,,,,,, '
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Initialisation des UMI,UMIT avant accumulation
         DO 8005 K=KDEB,KFIN
            KP=K-KDEB+1
            UMI(KP,1)=zxzini
            UMI(KP,2)=zxzini
            UMI(KP,3)=zxzini
            GRADGX(KP,1)=zxzini
            GRADGX(KP,2)=zxzini
            GRADGX(KP,3)=zxzini
            GRADGY(KP,1)=zxzini
            GRADGY(KP,2)=zxzini
            GRADGY(KP,3)=zxzini
            GRADGZ(KP,1)=zxzini
            GRADGZ(KP,2)=zxzini
            GRADGZ(KP,3)=zxzini
            AKM(KP)=zxzini
            AEM(KP)=zxzini
            UPIK(KP,1)=zxzini
            UPIK(KP,2)=zxzini
            UPIK(KP,3)=zxzini
            UPIE(KP,1)=zxzini
            UPIE(KP,2)=zxzini
            UPIE(KP,3)=zxzini
            GRADK(KP,1)=zxzini
            GRADK(KP,2)=zxzini
            GRADK(KP,3)=zxzini
            GRADE(KP,1)=zxzini
            GRADE(KP,2)=zxzini
            GRADE(KP,3)=zxzini
            DUDX(KP)=0.D0
            DUDY(KP)=0.D0
            DUDZ(KP)=0.D0
            DVDX(KP)=0.D0
            DVDY(KP)=0.D0
            DVDZ(KP)=0.D0
            DWDX(KP)=0.D0
            DWDY(KP)=0.D0
            DWDZ(KP)=0.D0
            DTDX(KP)=0.D0
            DTDY(KP)=0.D0
            DTDZ(KP)=0.D0
            C1P(KP)=C1
            C1G(KP)=0.D0
            PG(KP)=0.D0
 8005    CONTINUE
 
         DO 8006 I=1,NP
            DO 5006 K=KDEB,KFIN
               KP=K-KDEB+1
               NF=IPADL(LE(I,K))
               DR(KP,I)=DRR(I,K)
               UMI(KP,1)=UMI(KP,1)+UN(NF,1)*DR(KP,I)
               UMI(KP,2)=UMI(KP,2)+UN(NF,2)*DR(KP,I)
               UMI(KP,3)=UMI(KP,3)+UN(NF,3)*DR(KP,I)
               UIX(KP,I)=UN(NF,1)
               UIY(KP,I)=UN(NF,2)
               UIZ(KP,I)=UN(NF,3)
               GIX(KP,I)=GN(NF,1)
               GIY(KP,I)=GN(NF,2)
               GIZ(KP,I)=GN(NF,3)
               AK(KP,I)=TK(NF)
               AE(KP,I)=TE(NF)
               GRADGX(KP,1)=GRADGX(KP,1)+HR(K,I,1)*GIX(KP,I)
               GRADGX(KP,2)=GRADGX(KP,2)+HR(K,I,2)*GIX(KP,I)
               GRADGX(KP,3)=GRADGX(KP,3)+HR(K,I,3)*GIX(KP,I)
               GRADGY(KP,1)=GRADGY(KP,1)+HR(K,I,1)*GIY(KP,I)
               GRADGY(KP,2)=GRADGY(KP,2)+HR(K,I,2)*GIY(KP,I)
               GRADGY(KP,3)=GRADGY(KP,3)+HR(K,I,3)*GIY(KP,I)
               GRADGZ(KP,1)=GRADGZ(KP,1)+HR(K,I,1)*GIZ(KP,I)
               GRADGZ(KP,2)=GRADGZ(KP,2)+HR(K,I,2)*GIZ(KP,I)
               GRADGZ(KP,3)=GRADGZ(KP,3)+HR(K,I,3)*GIZ(KP,I)
               AKM(KP)=AKM(KP)+TK(NF)*DR(KP,I)
               AEM(KP)=AEM(KP)+TE(NF)*DR(KP,I)
               GRADK(KP,1)=GRADK(KP,1)+HR(K,I,1)*AK(KP,I)
               GRADK(KP,2)=GRADK(KP,2)+HR(K,I,2)*AK(KP,I)
               GRADK(KP,3)=GRADK(KP,3)+HR(K,I,3)*AK(KP,I)
               GRADE(KP,1)=GRADE(KP,1)+HR(K,I,1)*AE(KP,I)
               GRADE(KP,2)=GRADE(KP,2)+HR(K,I,2)*AE(KP,I)
               GRADE(KP,3)=GRADE(KP,3)+HR(K,I,3)*AE(KP,I)
C *** CALCUL DU TERME DE PRODUCTION ****************************
               DUDX(KP) = DUDX(KP)+HR(K,I,1)*UN(NF,1)
               DUDY(KP) = DUDY(KP)+HR(K,I,2)*UN(NF,1)
               DUDZ(KP) = DUDZ(KP)+HR(K,I,3)*UN(NF,1)
               DVDX(KP) = DVDX(KP)+HR(K,I,1)*UN(NF,2)
               DVDY(KP) = DVDY(KP)+HR(K,I,2)*UN(NF,2)
               DVDZ(KP) = DVDZ(KP)+HR(K,I,3)*UN(NF,2)
               DWDX(KP) = DWDX(KP)+HR(K,I,1)*UN(NF,3)
               DWDY(KP) = DWDY(KP)+HR(K,I,2)*UN(NF,3)
               DWDZ(KP) = DWDZ(KP)+HR(K,I,3)*UN(NF,3)
 5006       CONTINUE
 8006    CONTINUE
 
C     WRITE(IOIMP,*)'Initialisation SAF1,SAF2,SBF '
C
C Initialisation des variables d'accumulation SAF1,SAF2,SBF
C
         IF(IKOMP.EQ.0)THEN
 
            DO 80060 K=KDEB,KFIN
               KP=K-KDEB+1
               ARO(KP)=1.D0
80060       CONTINUE
 
            IF(IKAS.EQ.2)THEN
               DO 80061 I=1,NP
                  DO 50061 K=KDEB,KFIN
                     KP=K-KDEB+1
                     SAF1(KP,I)=0.D0
                     SAF2(KP,I)=0.D0
                     SAF3(KP,I)=0.D0
                     SAFK(KP,I)=0.D0
                     SAFE(KP,I)=0.D0
50061             CONTINUE
80061          CONTINUE
            ELSEIF(IKAS.EQ.3)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
               DO 80062 I=1,NP
                  DO 50062 K=KDEB,KFIN
                     KP=K-KDEB+1
                     SAF1(KP,I)=(-RGX(KP))*DR(KP,I)
                     SAF2(KP,I)=(-RGY(KP))*DR(KP,I)
                     SAF3(KP,I)=(-RGZ(KP))*DR(KP,I)
                     SAFK(KP,I)=0.D0
                     SAFE(KP,I)=0.D0
                     NKG=1+(1-IKG)*(NK-1)
                     DTDX(KP)=DTDX(KP)+HR(K,I,1)*RGE(NKG,1)
                     DTDY(KP)=DTDY(KP)+HR(K,I,2)*RGE(NKG,2)
                     DTDZ(KP)=DTDZ(KP)+HR(K,I,3)*RGE(NKG,3)
50062             CONTINUE
80062          CONTINUE
            ELSEIF(IKAS.EQ.5)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
               DO 80063 I=1,NP
                  DO 50063 K=KDEB,KFIN
                     KP=K-KDEB+1
                     NF=1+(1-IKT)*(IPADL(LE(I,K))-1)
                     NFR=1+(1-IKREF)*(IPADL(LE(I,K))-1)
                     SAF1(KP,I)=(-RGX(KP)*(TREF(NFR)-TN(NF)))*DR(KP,I)
                     SAF2(KP,I)=(-RGY(KP)*(TREF(NFR)-TN(NF)))*DR(KP,I)
                     SAF3(KP,I)=(-RGZ(KP)*(TREF(NFR)-TN(NF)))*DR(KP,I)
                     SAFK(KP,I)=0.D0
                     SAFE(KP,I)=0.D0
                     DTDX(KP)=DTDX(KP)+HR(K,I,1)*TN(NF)
                     DTDY(KP)=DTDY(KP)+HR(K,I,2)*TN(NF)
                     DTDZ(KP)=DTDZ(KP)+HR(K,I,3)*TN(NF)
50063             CONTINUE
80063          CONTINUE
            ENDIF
 
         ELSEIF(IKOMP.EQ.1)THEN
 
            DO 80020 K=KDEB,KFIN
               KP=K-KDEB+1
               NK=K+K0
               NKR=1+(1-IKR)*(NK-1)
               ARO(KP)=RO(NKR)
80020       CONTINUE
 
            IF(IKAS.EQ.4)THEN
               DO 80064 I=1,NP
                  DO 50064 K=KDEB,KFIN
                     KP=K-KDEB+1
                     SAF1(KP,I)=0.D0
                     SAF2(KP,I)=0.D0
                     SAF3(KP,I)=0.D0
                     SAFK(KP,I)=0.D0
                     SAFE(KP,I)=0.D0
50064             CONTINUE
80064          CONTINUE
            ELSEIF(IKAS.EQ.5)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
               DO 80065 I=1,NP
                  DO 50065 K=KDEB,KFIN
                     KP=K-KDEB+1
                     SAF1(KP,I)=-RGX(KP)*DR(KP,I)
                     SAF2(KP,I)=-RGY(KP)*DR(KP,I)
                     SAF3(KP,I)=-RGZ(KP)*DR(KP,I)
                     SAFK(KP,I)=0.D0
                     SAFE(KP,I)=0.D0
50065             CONTINUE
80065          CONTINUE
            ENDIF
         ENDIF
 
         DO 8007 K=KDEB,KFIN
            KP=K-KDEB+1
            NK=K+K0
            AKM(KP)=ABS(AKM(KP))/AIRE(KP)+zxzini
            AEM(KP)=ABS(AEM(KP))/AIRE(KP)+zxzini
            UMI(KP,1)=UMI(KP,1)/AIRE(KP)
            UMI(KP,2)=UMI(KP,2)/AIRE(KP)
            UMI(KP,3)=UMI(KP,3)/AIRE(KP)
            UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2)
     &            +UMI(KP,3)*UMI(KP,3)
            UM(KP)=SQRT(UM(KP))
            AX=GRADGX(KP,1)*GRADGX(KP,1)+GRADGX(KP,2)*GRADGX(KP,2)
     &        +GRADGX(KP,3)*GRADGX(KP,3)
            AX=SQRT(AX)+zxzini
            GRADGX(KP,1)=GRADGX(KP,1)/AX
            GRADGX(KP,2)=GRADGX(KP,2)/AX
            GRADGX(KP,3)=GRADGX(KP,3)/AX
            AY=GRADGY(KP,1)*GRADGY(KP,1)+GRADGY(KP,2)*GRADGY(KP,2)
     &        +GRADGY(KP,3)*GRADGY(KP,3)
            AY=SQRT(AY)+zxzini
            GRADGY(KP,1)=GRADGY(KP,1)/AY
            GRADGY(KP,2)=GRADGY(KP,2)/AY
            GRADGY(KP,3)=GRADGY(KP,3)/AY
            AZ=GRADGZ(KP,1)*GRADGZ(KP,1)+GRADGZ(KP,2)*GRADGZ(KP,2)
     &        +GRADGZ(KP,3)*GRADGZ(KP,3)
            AZ=SQRT(AZ)+zxzini
            GRADGZ(KP,1)=GRADGZ(KP,1)/AZ
            GRADGZ(KP,2)=GRADGZ(KP,2)/AZ
            GRADGZ(KP,3)=GRADGZ(KP,3)/AZ
 
            UPGX(KP)=GRADGX(KP,1)*UMI(KP,1)+GRADGX(KP,2)*UMI(KP,2)
     &              +GRADGX(KP,3)*UMI(KP,3)
            UPIGX(KP,1)=UPGX(KP)*GRADGX(KP,1)
            UPIGX(KP,2)=UPGX(KP)*GRADGX(KP,2)
            UPIGX(KP,3)=UPGX(KP)*GRADGX(KP,3)
            UPGY(KP)=GRADGY(KP,1)*UMI(KP,1)+GRADGY(KP,2)*UMI(KP,2)
     &              +GRADGY(KP,3)*UMI(KP,3)
            UPIGY(KP,1)=UPGY(KP)*GRADGY(KP,1)
            UPIGY(KP,2)=UPGY(KP)*GRADGY(KP,2)
            UPIGY(KP,3)=UPGY(KP)*GRADGY(KP,3)
            UPGZ(KP)=GRADGZ(KP,1)*UMI(KP,1)+GRADGZ(KP,2)*UMI(KP,2)
     &              +GRADGZ(KP,3)*UMI(KP,3)
            UPIGZ(KP,1)=UPGZ(KP)*GRADGZ(KP,1)
            UPIGZ(KP,2)=UPGZ(KP)*GRADGZ(KP,2)
            UPIGZ(KP,3)=UPGZ(KP)*GRADGZ(KP,3)
 
            AX=GRADK(KP,1)*GRADK(KP,1)+GRADK(KP,2)*GRADK(KP,2)
     &        +GRADK(KP,3)*GRADK(KP,3)
            AX=SQRT(AX)+zxzini
            GRADK(KP,1)=GRADK(KP,1)/AX
            GRADK(KP,2)=GRADK(KP,2)/AX
            GRADK(KP,3)=GRADK(KP,3)/AX
            UPK(KP)=GRADK(KP,1)*UMI(KP,1)+GRADK(KP,2)*UMI(KP,2)
     &             +GRADK(KP,3)*UMI(KP,3)
            UPIK(KP,1)=UPK(KP)*GRADK(KP,1)
            UPIK(KP,2)=UPK(KP)*GRADK(KP,2)
            UPIK(KP,3)=UPK(KP)*GRADK(KP,3)
 
            AX=GRADE(KP,1)*GRADE(KP,1)+GRADE(KP,2)*GRADE(KP,2)
     &        +GRADE(KP,3)*GRADE(KP,3)
            AX=SQRT(AX)+zxzini
            GRADE(KP,1)=GRADE(KP,1)/AX
            GRADE(KP,2)=GRADE(KP,2)/AX
            GRADE(KP,3)=GRADE(KP,3)/AX
            UPE(KP)=GRADE(KP,1)*UMI(KP,1)+GRADE(KP,2)*UMI(KP,2)
     &             +GRADE(KP,3)*UMI(KP,3)
            UPIE(KP,1)=UPE(KP)*GRADE(KP,1)
            UPIE(KP,2)=UPE(KP)*GRADE(KP,2)
            UPIE(KP,3)=UPE(KP)*GRADE(KP,3)
 
            ANUT(KP)=CNU*AKM(KP)*AKM(KP)/AEM(KP)
            AMUT(NK)=ARO(KP)*ANUT(KP)
            COEFT(KP)=COEF(KP)+AMUT(NK)
            COEFK(KP)=COEF(KP)+ANUT(KP)
            COEFE(KP)=COEF(KP)+ANUT(KP)/SGE
C --- 2.D0 CALCUL DU TERME DE PRODUCTION PROPREMENT DIT ----------
            PTEMP=( (DUDY(KP)+DVDX(KP))**2 + (DUDZ(KP)+DWDX(KP))**2
     &            + (DVDZ(KP)+DWDY(KP))**2
     &            + 2.D0*(DUDX(KP)**2+DVDY(KP)**2+DWDZ(KP)**2) )
            P(KP)=PTEMP
 8007    CONTINUE
 
         IF(IKAS.EQ.5.AND.IKOMP.EQ.0)THEN
            IF(IMODEL.EQ.1)THEN
               DO 80113 K=KDEB,KFIN
                  KP=K+1-KDEB
                  PG(KP)=(RGX(KP)*DTDX(KP)+RGY(KP)*DTDY(KP)
     &                   +RGZ(KP)*DTDZ(KP) )/SGT
                  IF(PG(KP).GT.0.D0)C1G(KP)=C1
80113          CONTINUE
            ELSEIF(IMODEL.EQ.2)THEN
               DO 80213 K=KDEB,KFIN
                  KP=K+1-KDEB
                  PG(KP)=(RGX(KP)*DTDX(KP)+RGY(KP)*DTDY(KP)
     &                   +RGZ(KP)*DTDZ(KP) )/SGT
                  RF    =-PG(KP)/(P(KP)+PG(KP)+zxzini)
                  IF(PG(KP).GT.0.D0)RF=0.D0
                  C1G(KP)=C1*(1.D0+C3*RF)
                  C1P(KP)=C1*(1.D0+C3*RF)
80213          CONTINUE
            ENDIF
         ELSEIF(IKAS.EQ.2.AND.IKOMP.EQ.0)THEN
            IF(IMODEL.EQ.1)THEN
               DO 80114 K=KDEB,KFIN
                  KP=K+1-KDEB
                  PG(KP)=(DTDX(KP)+DTDY(KP)+DTDZ(KP))/SGT
                  IF(PG(KP).GT.0.D0)C1G(KP)=C1
80114          CONTINUE
            ELSEIF(IMODEL.EQ.2)THEN
               DO 80214 K=KDEB,KFIN
                  KP=K+1-KDEB
                  PG(KP)=(DTDX(KP)+DTDY(KP)+DTDZ(KP))/SGT
                  RF    =-PG(KP)/(P(KP)+PG(KP)+zxzini)
                  IF(PG(KP).GT.0.D0)RF=0.D0
                  C1G(KP)=C1*(1.D0+C3*RF    )
                  C1P(KP)=C1*(1.D0+C3*RF    )
80214          CONTINUE
            END IF
         END IF
 
C --------------------------------------------------------------
C *** FIN DU CALCUL DU TERME DE PRODUCTION *********************
 
         DO 80074 K=KDEB,KFIN
            KP=K-KDEB+1
            BMX=UMI(KP,1)/AL(KP)
            BMY=UMI(KP,2)/AH(KP)
            BMZ=UMI(KP,3)/AP(KP)
            BM(KP)=BMX*BMX+BMY*BMY+BMZ*BMZ
            BM(KP)=SQRT(BM(KP))+zxzini
 
            BPGXX=UPIGX(KP,1)/AL(KP)
            BPGXY=UPIGX(KP,2)/AH(KP)
            BPGXZ=UPIGX(KP,3)/AP(KP)
            BPGX(KP)=BPGXX*BPGXX+BPGXY*BPGXY+BPGXZ*BPGXZ
            BPGX(KP)=SQRT(BPGX(KP))+zxzini
            BPGYX=UPIGY(KP,1)/AL(KP)
            BPGYY=UPIGY(KP,2)/AH(KP)
            BPGYZ=UPIGY(KP,3)/AP(KP)
            BPGY(KP)=BPGYX*BPGYX+BPGYY*BPGYY+BPGYZ*BPGYZ
            BPGY(KP)=SQRT(BPGY(KP))+zxzini
            BPGZX=UPIGZ(KP,1)/AL(KP)
            BPGZY=UPIGZ(KP,2)/AH(KP)
            BPGZZ=UPIGZ(KP,3)/AP(KP)
            BPGZ(KP)=BPGZX*BPGZX+BPGZY*BPGZY+BPGZZ*BPGZZ
            BPGZ(KP)=SQRT(BPGZ(KP))+zxzini
 
            BPKX=UPIK(KP,1)/AL(KP)
            BPKY=UPIK(KP,2)/AH(KP)
            BPKZ=UPIK(KP,3)/AP(KP)
            BPK(KP)=BPKX*BPKX+BPKY*BPKY+BPKZ*BPKZ
            BPK(KP)=SQRT(BPK(KP))+zxzini
            BPEX=UPIE(KP,1)/AL(KP)
            BPEY=UPIE(KP,2)/AH(KP)
            BPEZ=UPIE(KP,3)/AP(KP)
            BPE(KP)=BPEX*BPEX+BPEY*BPEY+BPEZ*BPEZ
            BPE(KP)=SQRT(BPE(KP))+zxzini
80074    CONTINUE
C
C     AVANT DECENTREMENT
C
C ALFA,AKSI,CCU utilisées seulement dans la boucle 7008
         IF(IDCEN.EQ.1)THEN
            DO 80081 K=KDEB,KFIN
               KP=K-KDEB+1
               GXCXT(KP)=0.D0
               GXCYT(KP)=0.D0
               GXCZT(KP)=0.D0
               GXCXY(KP)=0.D0
               GXCXZ(KP)=0.D0
               GXCYZ(KP)=0.D0
               GYCXT(KP)=0.D0
               GYCYT(KP)=0.D0
               GYCZT(KP)=0.D0
               GYCXY(KP)=0.D0
               GYCXZ(KP)=0.D0
               GYCYZ(KP)=0.D0
               GZCXT(KP)=0.D0
               GZCYT(KP)=0.D0
               GZCZT(KP)=0.D0
               GZCXY(KP)=0.D0
               GZCXZ(KP)=0.D0
               GZCYZ(KP)=0.D0
               GKCXT(KP)=0.D0
               GKCYT(KP)=0.D0
               GKCZT(KP)=0.D0
               GKCXY(KP)=0.D0
               GKCXZ(KP)=0.D0
               GKCYZ(KP)=0.D0
               GECXT(KP)=0.D0
               GECYT(KP)=0.D0
               GECZT(KP)=0.D0
               GECXY(KP)=0.D0
               GECXZ(KP)=0.D0
               GECYZ(KP)=0.D0
80081       CONTINUE
 
         ELSEIF(IDCEN.EQ.2)THEN
            DO 8008 K=KDEB,KFIN
               KP=K-KDEB+1
               HMK=2.D0*UM(KP)/BM(KP)
               XMB(KP)=HMK
               ALFA=UM(KP)*HMK/COEFT(KP)/2.D0
               AKSI=ALFA/3.D0
               IF(ALFA.GT.3.D0)AKSI=1.D0
               CCT=AKSI/BM(KP)
 
               HPK=2.D0*UPGX(KP)/BPGX(KP)
               ALFA=UPGX(KP)*HPK/COEFT(KP)/2.D0
               AKSI=ALFA/3.D0
               IF(ALFA.GT.3.D0)AKSI=1.D0
               CCP=AKSI/BPGX(KP)
               CPT=CCP-CCT
               CC2X=0.D0
               IF(CPT.GE.0.D0)CC2X=CPT
 
               HPK=2.D0*UPGY(KP)/BPGY(KP)
               ALFA=UPGY(KP)*HPK/COEFT(KP)/2.D0
               AKSI=ALFA/3.D0
               IF(ALFA.GT.3.D0)AKSI=1.D0
               CCP=AKSI/BPGY(KP)
               CPT=CCP-CCT
               CC2Y=0.D0
               IF(CPT.GE.0.D0)CC2Y=CPT
 
               HPK=2.D0*UPGZ(KP)/BPGZ(KP)
               ALFA=UPGZ(KP)*HPK/COEFT(KP)/2.D0
               AKSI=ALFA/3.D0
               IF(ALFA.GT.3.D0)AKSI=1.D0
               CCP=AKSI/BPGZ(KP)
               CPT=CCP-CCT
               CC2Z=0.D0
               IF(CPT.GE.0.D0)CC2Z=CPT
 
               HPK=2.D0*UPK(KP)/BPK(KP)
               ALFA=UPK(KP)*HPK/COEFK(KP)/2.D0
               AKSI=ALFA/3.D0
               IF(ALFA.GT.3.D0)AKSI=1.D0
               CCP=AKSI/BPK(KP)
               CPT=CCP-CCT
               CC2K=0.D0
               IF(CPT.GE.0.D0)CC2K=CPT
 
               HPK=2.D0*UPE(KP)/BPE(KP)
               ALFA=UPE(KP)*HPK/COEFE(KP)/2.D0
               AKSI=ALFA/3.D0
               IF(ALFA.GT.3.D0)AKSI=1.D0
               CCP=AKSI/BPE(KP)
               CPT=CCP-CCT
               CC2E=0.D0
               IF(CPT.GE.0.D0)CC2E=CPT
 
               GXCXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT+UPIGX(KP,1)*UPIGX(KP,1
     $              )*CC2X)
               GXCYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT+UPIGX(KP,2)*UPIGX(KP,2
     $              )*CC2X)
               GXCZT(KP)=(UMI(KP,3)*UMI(KP,3)*CCT+UPIGX(KP,3)*UPIGX(KP,3
     $              )*CC2X)
               GXCXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT+UPIGX(KP,1)*UPIGX(KP,2
     $              )*CC2X)
               GXCXZ(KP)=(UMI(KP,1)*UMI(KP,3)*CCT+UPIGX(KP,1)*UPIGX(KP,3
     $              )*CC2X)
               GXCYZ(KP)=(UMI(KP,2)*UMI(KP,3)*CCT+UPIGX(KP,2)*UPIGX(KP,3
     $              )*CC2X)
 
               GYCXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT+UPIGY(KP,1)*UPIGY(KP,1
     $              )*CC2Y)
               GYCYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT+UPIGY(KP,2)*UPIGY(KP,2
     $              )*CC2Y)
               GYCZT(KP)=(UMI(KP,3)*UMI(KP,3)*CCT+UPIGY(KP,3)*UPIGY(KP,3
     $              )*CC2Y)
               GYCXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT+UPIGY(KP,1)*UPIGY(KP,2
     $              )*CC2Y)
               GYCXZ(KP)=(UMI(KP,1)*UMI(KP,3)*CCT+UPIGY(KP,1)*UPIGY(KP,3
     $              )*CC2Y)
               GYCYZ(KP)=(UMI(KP,2)*UMI(KP,3)*CCT+UPIGY(KP,2)*UPIGY(KP,3
     $              )*CC2Y)
 
               GZCXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT+UPIGZ(KP,1)*UPIGZ(KP,1
     $              )*CC2Z)
               GZCYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT+UPIGZ(KP,2)*UPIGZ(KP,2
     $              )*CC2Z)
               GZCZT(KP)=(UMI(KP,3)*UMI(KP,3)*CCT+UPIGZ(KP,3)*UPIGZ(KP,3
     $              )*CC2Z)
               GZCXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT+UPIGZ(KP,1)*UPIGZ(KP,2
     $              )*CC2Z)
               GZCXZ(KP)=(UMI(KP,1)*UMI(KP,3)*CCT+UPIGZ(KP,1)*UPIGZ(KP,3
     $              )*CC2Z)
               GZCYZ(KP)=(UMI(KP,2)*UMI(KP,3)*CCT+UPIGZ(KP,2)*UPIGZ(KP,3
     $              )*CC2Z)
 
               GKCXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT+UPIK(KP,1)*UPIK(KP,1)
     $              *CC2K)
               GKCYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT+UPIK(KP,2)*UPIK(KP,2)
     $              *CC2K)
               GKCZT(KP)=(UMI(KP,3)*UMI(KP,3)*CCT+UPIK(KP,3)*UPIK(KP,3)
     $              *CC2K)
               GKCXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT+UPIK(KP,1)*UPIK(KP,2)
     $              *CC2K)
               GKCXZ(KP)=(UMI(KP,1)*UMI(KP,3)*CCT+UPIK(KP,1)*UPIK(KP,3)
     $              *CC2K)
               GKCYZ(KP)=(UMI(KP,2)*UMI(KP,3)*CCT+UPIK(KP,2)*UPIK(KP,3)
     $              *CC2K)
 
               GECXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT+UPIE(KP,1)*UPIE(KP,1)
     $              *CC2E)
               GECYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT+UPIE(KP,2)*UPIE(KP,2)
     $              *CC2E)
               GECZT(KP)=(UMI(KP,3)*UMI(KP,3)*CCT+UPIE(KP,3)*UPIE(KP,3)
     $              *CC2E)
               GECXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT+UPIE(KP,1)*UPIE(KP,2)
     $              *CC2E)
               GECXZ(KP)=(UMI(KP,1)*UMI(KP,3)*CCT+UPIE(KP,1)*UPIE(KP,3)
     $              *CC2E)
               GECYZ(KP)=(UMI(KP,2)*UMI(KP,3)*CCT+UPIE(KP,2)*UPIE(KP,3)
     $              *CC2E)
 
 8008       CONTINUE
 
         ELSEIF(IDCEN.EQ.3)THEN
            DO 8018 K=KDEB,KFIN
               KP=K-KDEB+1
               HMK=2.D0*UM(KP)/BM(KP)
               XMB(KP)=HMK
               ALFA=UM(KP)*HMK/COEFT(KP)/2.D0
               AKSI=ALFA/3.D0
               IF(ALFA.GT.3.D0)AKSI=1.D0
               CCT=AKSI/BM(KP)
 
               GXCXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT
               GXCYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT
               GXCZT(KP)=UMI(KP,3)*UMI(KP,3)*CCT
               GXCXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT
               GXCXZ(KP)=UMI(KP,1)*UMI(KP,3)*CCT
               GXCYZ(KP)=UMI(KP,2)*UMI(KP,3)*CCT
 
               GYCXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT
               GYCYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT
               GYCZT(KP)=UMI(KP,3)*UMI(KP,3)*CCT
               GYCXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT
               GYCXZ(KP)=UMI(KP,1)*UMI(KP,3)*CCT
               GYCYZ(KP)=UMI(KP,2)*UMI(KP,3)*CCT
 
               GZCXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT
               GZCYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT
               GZCZT(KP)=UMI(KP,3)*UMI(KP,3)*CCT
               GZCXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT
               GZCXZ(KP)=UMI(KP,1)*UMI(KP,3)*CCT
               GZCYZ(KP)=UMI(KP,2)*UMI(KP,3)*CCT
 
               GKCXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT
               GKCYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT
               GKCZT(KP)=UMI(KP,3)*UMI(KP,3)*CCT
               GKCXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT
               GKCXZ(KP)=UMI(KP,1)*UMI(KP,3)*CCT
               GKCYZ(KP)=UMI(KP,2)*UMI(KP,3)*CCT
 
               GECXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT
               GECYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT
               GECZT(KP)=UMI(KP,3)*UMI(KP,3)*CCT
               GECXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT
               GECXZ(KP)=UMI(KP,1)*UMI(KP,3)*CCT
               GECYZ(KP)=UMI(KP,2)*UMI(KP,3)*CCT
 
 8018       CONTINUE
 
         ELSEIF(IDCEN.EQ.4)THEN
            DT19=DTM1*0.5D0
            DO 8009 K=KDEB,KFIN
               KP=K-KDEB+1
               XMB(KP)=(AL(KP)+AH(KP))/2.D0
 
               GXCXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19
               GXCYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19
               GXCZT(KP)=UMI(KP,3)*UMI(KP,3)*DT19
               GXCXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19
               GXCXZ(KP)=UMI(KP,1)*UMI(KP,3)*DT19
               GXCYZ(KP)=UMI(KP,2)*UMI(KP,3)*DT19
 
               GYCXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19
               GYCYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19
               GYCZT(KP)=UMI(KP,3)*UMI(KP,3)*DT19
               GYCXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19
               GYCXZ(KP)=UMI(KP,1)*UMI(KP,3)*DT19
               GYCYZ(KP)=UMI(KP,2)*UMI(KP,3)*DT19
 
               GZCXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19
               GZCYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19
               GZCZT(KP)=UMI(KP,3)*UMI(KP,3)*DT19
               GZCXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19
               GZCXZ(KP)=UMI(KP,1)*UMI(KP,3)*DT19
               GZCYZ(KP)=UMI(KP,2)*UMI(KP,3)*DT19
 
               GKCXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19
               GKCYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19
               GKCZT(KP)=UMI(KP,3)*UMI(KP,3)*DT19
               GKCXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19
               GKCXZ(KP)=UMI(KP,1)*UMI(KP,3)*DT19
               GKCYZ(KP)=UMI(KP,2)*UMI(KP,3)*DT19
 
               GECXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19
               GECYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19
               GECZT(KP)=UMI(KP,3)*UMI(KP,3)*DT19
               GECXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19
               GECXZ(KP)=UMI(KP,1)*UMI(KP,3)*DT19
               GECYZ(KP)=UMI(KP,2)*UMI(KP,3)*DT19
 
 8009       CONTINUE
 
         ENDIF
 
C     WRITE(IOIMP,*)' Avt Calcul DT'
C
C     AVANT CALCUL DT
C
         DO 8010 K=KDEB,KFIN
            KP=K-KDEB+1
            DT0=DT
            DT1X=2.D0/
     &           (UMI(KP,1)*UMI(KP,1)/(GXCXT(KP)+COEFT(KP))
     &           +UMI(KP,2)*UMI(KP,2)/(GXCYT(KP)+COEFT(KP))
     &           +UMI(KP,3)*UMI(KP,3)/(GXCZT(KP)+COEFT(KP)))
            DT1Y=2.D0/
     &           (UMI(KP,1)*UMI(KP,1)/(GYCXT(KP)+COEFT(KP))
     &           +UMI(KP,2)*UMI(KP,2)/(GYCYT(KP)+COEFT(KP))
     &           +UMI(KP,3)*UMI(KP,3)/(GYCZT(KP)+COEFT(KP)))
            DT1Z=2.D0/
     &           (UMI(KP,1)*UMI(KP,1)/(GZCXT(KP)+COEFT(KP))
     &           +UMI(KP,2)*UMI(KP,2)/(GZCYT(KP)+COEFT(KP))
     &           +UMI(KP,3)*UMI(KP,3)/(GZCZT(KP)+COEFT(KP)))
            DT2X=0.5D0/
     &           ((GXCXT(KP)+COEFT(KP))*AL2(KP)
     &           +(GXCYT(KP)+COEFT(KP))*AH2(KP)
     &           +(GXCZT(KP)+COEFT(KP))*AP2(KP))
            DT2Y=0.5D0/
     &           ((GYCXT(KP)+COEFT(KP))*AL2(KP)
     &           +(GYCYT(KP)+COEFT(KP))*AH2(KP)
     &           +(GYCZT(KP)+COEFT(KP))*AP2(KP))
            DT2Z=0.5D0/
     &           ((GZCXT(KP)+COEFT(KP))*AL2(KP)
     &           +(GZCYT(KP)+COEFT(KP))*AH2(KP)
     &           +(GZCZT(KP)+COEFT(KP))*AP2(KP))
            IF(DT1X.LT.DT)DT=DT1X
            IF(DT1Y.LT.DT)DT=DT1Y
            IF(DT1Z.LT.DT)DT=DT1Z
            IF(DT2X.LT.DT)DT=DT2X
            IF(DT2Y.LT.DT)DT=DT2Y
            IF(DT2Z.LT.DT)DT=DT2Z
C     IF(DT3.LT.DT)DT=DT3
            IF(DT.NE.DT0) THEN
               DTT1=MIN(DT1X,DT1Y,DT1Z)
               DTT2=MIN(DT2X,DT2Y,DT2Z)
C        DTT3=DT3
               DIAEL=XMB(KP)
               NUEL=K
            END IF
 8010    CONTINUE
 
C     WRITE(IOIMP,*)' le coeur du calcul '
C Le coeur du calcul ...
 
         IF(IKOMP.EQ.0)THEN
 
            DO 8014 I=1,NP
               DO 80141 J=1,NP
                  DO 80142 K=KDEB,KFIN
                     KP=K-KDEB+1
C     GEO1=(CF1(KP)*Q1(J,I)+CF2(KP)*Q2(J,I)+CF3(KP)*Q3(J,I))*CUB8
                     GEO1=CFM(KP)*QGGT(J,I)*CUB8
                     ZVGGX=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*GXCXT(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*GXCYT(KP)
     &               +    HR(K,I,3)*HR(K,J,3)*GXCZT(KP)
     &               +    HR(K,I,1)*HR(K,J,2)*GXCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,1)*GXCXY(KP)
     &               +    HR(K,I,1)*HR(K,J,3)*GXCXZ(KP)
     &               +    HR(K,I,3)*HR(K,J,1)*GXCXZ(KP)
     &               +    HR(K,I,2)*HR(K,J,3)*GXCYZ(KP)
     &               +    HR(K,I,3)*HR(K,J,2)*GXCYZ(KP) )
     &               +    (GXCXT(KP)+GXCYT(KP)+GXCZT(KP))/3.D0*GEO1
C    &+    (GXCXT(KP)+GXCYT(KP)+GXCZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8
 
                     ZVGGY=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*GYCXT(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*GYCYT(KP)
     &               +    HR(K,I,3)*HR(K,J,3)*GYCZT(KP)
     &               +    HR(K,I,1)*HR(K,J,2)*GYCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,1)*GYCXY(KP)
     &               +    HR(K,I,1)*HR(K,J,3)*GYCXZ(KP)
     &               +    HR(K,I,3)*HR(K,J,1)*GYCXZ(KP)
     &               +    HR(K,I,2)*HR(K,J,3)*GYCYZ(KP)
     &               +    HR(K,I,3)*HR(K,J,2)*GYCYZ(KP) )
     &               +    (GYCXT(KP)+GYCYT(KP)+GYCZT(KP))/3.D0*GEO1
C    &+    (GYCXT(KP)+GYCYT(KP)+GYCZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8
 
                     ZVGGZ=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*GZCXT(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*GZCYT(KP)
     &               +    HR(K,I,3)*HR(K,J,3)*GZCZT(KP)
     &               +    HR(K,I,1)*HR(K,J,2)*GZCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,1)*GZCXY(KP)
     &               +    HR(K,I,1)*HR(K,J,3)*GZCXZ(KP)
     &               +    HR(K,I,3)*HR(K,J,1)*GZCXZ(KP)
     &               +    HR(K,I,2)*HR(K,J,3)*GZCYZ(KP)
     &               +    HR(K,I,3)*HR(K,J,2)*GZCYZ(KP) )
     &               +    (GZCXT(KP)+GZCYT(KP)+GZCZT(KP))/3.D0*GEO1
C    &+    (GZCXT(KP)+GZCYT(KP)+GZCZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8
 
                     ZVGD=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*COEFT(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*COEFT(KP)
     &               +    HR(K,I,3)*HR(K,J,3)*COEFT(KP) )
     &               +    COEFT(KP)*XMH(KP)*GEO1
C    &+    COEFT(KP)*XMH(KP)*QGGT(J,I)*CUB8
 
                     ZVGK =AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*GKCXT(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*GKCYT(KP)
     &               +    HR(K,I,3)*HR(K,J,3)*GKCZT(KP)
     &               +    HR(K,I,1)*HR(K,J,2)*GKCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,1)*GKCXY(KP)
     &               +    HR(K,I,1)*HR(K,J,3)*GKCXZ(KP)
     &               +    HR(K,I,3)*HR(K,J,1)*GKCXZ(KP)
     &               +    HR(K,I,2)*HR(K,J,3)*GKCYZ(KP)
     &               +    HR(K,I,3)*HR(K,J,2)*GKCYZ(KP) )
     &               +    (GKCXT(KP)+GKCYT(KP)+GKCZT(KP))/3.D0*GEO1
C    &+    (GKCXT(KP)+GKCYT(KP)+GKCZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8
 
                     ZVGE =AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*GECXT(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*GECYT(KP)
     &               +    HR(K,I,3)*HR(K,J,3)*GECZT(KP)
     &               +    HR(K,I,1)*HR(K,J,2)*GECXY(KP)
     &               +    HR(K,I,2)*HR(K,J,1)*GECXY(KP)
     &               +    HR(K,I,1)*HR(K,J,3)*GECXZ(KP)
     &               +    HR(K,I,3)*HR(K,J,1)*GECXZ(KP)
     &               +    HR(K,I,2)*HR(K,J,3)*GECYZ(KP)
     &               +    HR(K,I,3)*HR(K,J,2)*GECYZ(KP) )
     &               +    (GECXT(KP)+GECYT(KP)+GECZT(KP))/3.D0*GEO1
C    &+    (GECXT(KP)+GECYT(KP)+GECZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8
 
                     ZVGDK=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*COEFK(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*COEFK(KP)
     &               +    HR(K,I,3)*HR(K,J,3)*COEFK(KP) )
     &               +    COEFK(KP)*XMH(KP)*GEO1
C    &+    COEFK(KP)*XMH(KP)*QGGT(J,I)*CUB8
 
                     ZVGDE=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*COEFE(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*COEFE(KP)
     &               +    HR(K,I,3)*HR(K,J,3)*COEFE(KP) )
     &               +    COEFE(KP)*XMH(KP)*GEO1
C    &+    COEFE(KP)*XMH(KP)*QGGT(J,I)*CUB8
 
Cnon  V2=(UIX(KP,I)*HR(K,J,1)+UIY(KP,I)*HR(K,J,2)+UIZ(KP,I)*HR(K,J,3))
C    &   *DR(KP,I)
 
                     V2=UMI(KP,1)*DR(KP,I)*HR(K,J,1)+UMI(KP,2)*DR(KP,I)
     $                    *HR(K,J,2)+UMI(KP,3)*DR(KP,I)*HR(K,J,3)
 
                     SAF1(KP,I)=SAF1(KP,I)+(V2+ZVGGX)*GIX(KP,J)+ZVGD
     $                    *UIX(KP,J)
                     SAF2(KP,I)=SAF2(KP,I)+(V2+ZVGGY)*GIY(KP,J)+ZVGD
     $                    *UIY(KP,J)
                     SAF3(KP,I)=SAF3(KP,I)+(V2+ZVGGZ)*GIZ(KP,J)+ZVGD
     $                    *UIZ(KP,J)
 
                     SAFK(KP,I)=SAFK(KP,I)+(V2+ZVGK)*AK(KP,J)+ZVGDK
     $                    *AK(KP,J)
                     SAFE(KP,I)=SAFE(KP,I)+(V2+ZVGE)*AE(KP,J)+ZVGDE
     $                    *AE(KP,J)
 
80142             CONTINUE
80141          CONTINUE
 8014       CONTINUE
 
         ELSEIF(IKOMP.EQ.1)THEN
C->
 
            DO 8015 I=1,NP
               DO 80151 J=1,NP
                  DO 80152 K=KDEB,KFIN
                     KP=K-KDEB+1
C     GEO1=(CF1(KP)*Q1(J,I)+CF2(KP)*Q2(J,I)+CF3(KP)*Q3(J,I))*CUB8
                     GEO1=CFM(KP)*QGGT(J,I)*CUB8
                     ZVGGX=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*GXCXT(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*GXCYT(KP)
     &               +    HR(K,I,3)*HR(K,J,3)*GXCZT(KP)
     &               +    HR(K,I,1)*HR(K,J,2)*GXCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,1)*GXCXY(KP)
     &               +    HR(K,I,1)*HR(K,J,3)*GXCXZ(KP)
     &               +    HR(K,I,3)*HR(K,J,1)*GXCXZ(KP)
     &               +    HR(K,I,2)*HR(K,J,3)*GXCYZ(KP)
     &               +    HR(K,I,3)*HR(K,J,2)*GXCYZ(KP) )
     &               +    (GXCXT(KP)+GXCYT(KP)+GXCZT(KP))/3.D0*GEO1
C    &+    (GXCXT(KP)+GXCYT(KP)+GXCZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8
 
                     ZVGGY=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*GYCXT(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*GYCYT(KP)
     &               +    HR(K,I,3)*HR(K,J,3)*GYCZT(KP)
     &               +    HR(K,I,1)*HR(K,J,2)*GYCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,1)*GYCXY(KP)
     &               +    HR(K,I,1)*HR(K,J,3)*GYCXZ(KP)
     &               +    HR(K,I,3)*HR(K,J,1)*GYCXZ(KP)
     &               +    HR(K,I,2)*HR(K,J,3)*GYCYZ(KP)
     &               +    HR(K,I,3)*HR(K,J,2)*GYCYZ(KP) )
     &               +    (GYCXT(KP)+GYCYT(KP)+GYCZT(KP))/3.D0*GEO1
C    &+    (GYCXT(KP)+GYCYT(KP)+GYCZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8
 
                     ZVGGZ=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*GZCXT(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*GZCYT(KP)
     &               +    HR(K,I,3)*HR(K,J,3)*GZCZT(KP)
     &               +    HR(K,I,1)*HR(K,J,2)*GZCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,1)*GZCXY(KP)
     &               +    HR(K,I,1)*HR(K,J,3)*GZCXZ(KP)
     &               +    HR(K,I,3)*HR(K,J,1)*GZCXZ(KP)
     &               +    HR(K,I,2)*HR(K,J,3)*GZCYZ(KP)
     &               +    HR(K,I,3)*HR(K,J,2)*GZCYZ(KP) )
     &               +    (GZCXT(KP)+GZCYT(KP)+GZCZT(KP))/3.D0*GEO1
C    &+    (GZCXT(KP)+GZCYT(KP)+GZCZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8
 
                     ZVGD=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*COEFT(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*COEFT(KP)
     &               +    HR(K,I,3)*HR(K,J,3)*COEFT(KP) )
     &               +    COEFT(KP)*GEO1
C    &+    COEFT(KP)*XMH(KP)*GEO1
C    &+    COEFT(KP)*XMH(KP)*QGGT(J,I)*CUB8
 
                     ZVGK =AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*GKCXT(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*GKCYT(KP)
     &               +    HR(K,I,3)*HR(K,J,3)*GKCZT(KP)
     &               +    HR(K,I,1)*HR(K,J,2)*GKCXY(KP)
     &               +    HR(K,I,2)*HR(K,J,1)*GKCXY(KP)
     &               +    HR(K,I,1)*HR(K,J,3)*GKCXZ(KP)
     &               +    HR(K,I,3)*HR(K,J,1)*GKCXZ(KP)
     &               +    HR(K,I,2)*HR(K,J,3)*GKCYZ(KP)
     &               +    HR(K,I,3)*HR(K,J,2)*GKCYZ(KP) )
     &               +    (GKCXT(KP)+GKCYT(KP)+GKCZT(KP))/3.D0*GEO1
C    &+    (GKCXT(KP)+GKCYT(KP)+GKCZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8
 
                     ZVGE =AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*GECXT(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*GECYT(KP)
     &               +    HR(K,I,3)*HR(K,J,3)*GECZT(KP)
     &               +    HR(K,I,1)*HR(K,J,2)*GECXY(KP)
     &               +    HR(K,I,2)*HR(K,J,1)*GECXY(KP)
     &               +    HR(K,I,1)*HR(K,J,3)*GECXZ(KP)
     &               +    HR(K,I,3)*HR(K,J,1)*GECXZ(KP)
     &               +    HR(K,I,2)*HR(K,J,3)*GECYZ(KP)
     &               +    HR(K,I,3)*HR(K,J,2)*GECYZ(KP) )
     &               +    (GECXT(KP)+GECYT(KP)+GECZT(KP))/3.D0*GEO1
C    &+    (GECXT(KP)+GECYT(KP)+GECZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8
 
                     ZVGDK=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*COEFK(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*COEFK(KP)
     &               +    HR(K,I,3)*HR(K,J,3)*COEFK(KP) )
     &               +    COEFK(KP)*GEO1
 
                     ZVGDE=AIRE(KP)*(
     &                    HR(K,I,1)*HR(K,J,1)*COEFE(KP)
     &               +    HR(K,I,2)*HR(K,J,2)*COEFE(KP)
     &               +    HR(K,I,3)*HR(K,J,3)*COEFE(KP) )
     &               +    COEFE(KP)*GEO1
C    &+    COEFE(KP)*XMH(KP)*QGGT(J,I)*CUB8
 
                     COEFT3=COEFT(KP)/3.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))*DR(KP,I)
 
C     V2=UMI(KP,1)*DR(KP,I)*HR(K,J,1)+UMI(KP,2)*DR(KP,I)*HR(K,J,2)
 
                     SAF1(KP,I)=SAF1(KP,I)+(V2+ZVGGX)*GIX(KP,J)+ZVGD
     $                    *UIX(KP,J)+ ZVGUU + ZVGUV + ZVGUW
                     SAF2(KP,I)=SAF2(KP,I)+(V2+ZVGGY)*GIY(KP,J)+ZVGD
     $                    *UIY(KP,J)+ ZVGVU + ZVGUV + ZVGVW
                     SAF3(KP,I)=SAF3(KP,I)+(V2+ZVGGZ)*GIZ(KP,J)+ZVGD
     $                    *UIZ(KP,J)+ ZVGWU + ZVGWV + ZVGWW
 
                     SAFK(KP,I)=SAFK(KP,I)+(V2+ZVGK)*AK(KP,J)+ZVGDK
     $                    *AK(KP,J)
                     SAFE(KP,I)=SAFE(KP,I)+(V2+ZVGE)*AE(KP,J)+ZVGDE
     $                    *AE(KP,J)
 
80152             CONTINUE
80151          CONTINUE
 8015       CONTINUE
 
         ENDIF
 
C
C --------------------------------------------------------------
C
C
         DO 8017 I=1,NP
            DO 80171 K=KDEB,KFIN
               KP=K-KDEB+1
               TSK       =(ANUT(KP)*(P(KP)+PG(KP))-AEM(KP))*DR(KP,I)
               SAFK(KP,I)=SAFK(KP,I)-TSK
 
               TSE       =(CNU*(C1P(KP)*P(KP)+C1G(KP)*PG(KP))*AKM(KP))
     $              *DR(KP,I)
               SAFE(KP,I)=SAFE(KP,I)-TSE
 
               TSE1      =C2*AEM(KP)/AKM(KP)*DR(KP,I)
               BE1(KP,I) =TSE1
80171       CONTINUE
 8017    CONTINUE
 
C
C
C --- ECRITURE DANS LES TABLEAUX F(1,*),F(2,*),G,GK,GE,GE1 (SCALAIRE)
         DO 80019 I=1,NP
            DO 50019 K=KDEB,KFIN
               KP=K+1-KDEB
               NF=IPADL(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)
               GK(NF) =GK(NF) +SAFK(KP,I)
               GE1(NF)=GE1(NF)+BE1(KP,I)
               GE(NF) =GE (NF)+SAFE(KP,I)
50019       CONTINUE
80019    CONTINUE
80001 CONTINUE
C
C *** FIN DE LA BOUCLE SUR LES PAQUETS D'ELEMENTS ***
C
C     WRITE(IOIMP,*)' FIN YCVIT GK='
C     WRITE(IOIMP,1002)(GK(i),i=1,nptd)
 
      IPAS=1
      RETURN
 
 1002 FORMAT(10(1X,1PE11.4))
      END