C KSUPG1    SOURCE    CHAT      05/01/13    01:09:14     5004           
      SUBROUTINE KSUPG1(WTILDE,UMI,CHGLD,CHGLP,KDEB,KFIN,LRV,
     &     TN,IPADI,UN,COEFT,NPTD,NEL,NP,DRR,HR,FN,
     &     AIRE,AL,AH,AP,IDCEN,IPADU,LE,CHGL,IKOMP,
     &     DTM1,DT,DTT1,DTT2,DIAEL,NUEL)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
-INC CCOPTIO
-INC CCREEL
      DIMENSION IPADU(*),LE(NP,1),IPADI(*)
      DIMENSION UN(NPTD,IDIM),TN(*)
      DIMENSION AIRE(LRV),AL(LRV),AH(LRV),AP(LRV),WTILDE(LRV,NP)
      DIMENSION COEFT(LRV),CHGLD(LRV),CHGLP(LRV)
      DIMENSION HR(NEL,NP,IDIM),DRR(NP,NEL),FN(NP,1)
      PARAMETER (LRV1=64)
      DIMENSION UM(LRV1),UMI(LRV,3),UMJ(LRV1,3)
      DIMENSION GRAD(LRV1,3)
      DIMENSION UP(LRV1),UPI(LRV1,3)
      DIMENSION BM(LRV1),BP(LRV1)
      DIMENSION TO1(LRV1),TO2(LRV1)
      DIMENSION CXT(LRV1),CYT(LRV1),CZT(LRV1)
      DIMENSION CXY(LRV1),CXZ(LRV1),CYZ(LRV1),CFM(LRV1)
      DIMENSION AL2(LRV1),AH2(LRV1),AP2(LRV1),XMB(LRV1)
      DATA IKV/1/,NPG/1/
      CC0=1.D0
      xpeti2 =  sqrt(XPETIT)
      IF(IKOMP.EQ.1)CC0=0.2D0
      IF(LRV.NE.LRV1)STOP
      DO 7005 N=1,IDIM
         DO 70051 K=KDEB,KFIN
            KP=K-KDEB+1
            UMI(KP,N)=xpeti2
            UMJ(KP,N)=xpeti2
            UPI(KP,N)=xpeti2
            GRAD(KP,N)=xpeti2
70051    CONTINUE
 7005 CONTINUE
      IF(IKV.EQ.1)THEN
         DO 70060 N=1,IDIM
            DO 70061 I=1,NP
C*IBMDIR* PREFER VECTOR
               DO 70062 K=KDEB,KFIN
                  KP=K-KDEB+1
                  NF=IPADU(LE(I,K))
                  NG=IPADI(LE(I,K))
                  UMI(KP,N)=UMI(KP,N)+UN(NF,N)*DRR(I,K)
                  UMJ(KP,N)=UMJ(KP,N)+UN(NF,N)*DRR(I,K)
                  GRAD(KP,N)=GRAD(KP,N)+HR(K,I,N)*TN(NG)
70062          CONTINUE
70061       CONTINUE
70060    CONTINUE
         DO 7017 K=KDEB,KFIN
            KP=K-KDEB+1
            UMJ(KP,1)=UMJ(KP,1)/AIRE(KP)
            UMJ(KP,2)=UMJ(KP,2)/AIRE(KP)
            UM(KP)=UMJ(KP,1)*UMJ(KP,1)+UMJ(KP,2)*UMJ(KP,2)
            UM(KP)=SQRT(UM(KP))
 7017    CONTINUE
      ELSE
         DO 70160 N=1,IDIM
            DO 70161 I=1,NP
C*IBMDIR* PREFER VECTOR
               DO 70162 K=KDEB,KFIN
                  KP=K-KDEB+1
                  NF=IPADU(LE(I,K))
                  NG=IPADI(LE(I,K))
                  UMI(KP,N)=UMI(KP,N)+UN(NF,N)*DRR(I,K)
                  GRAD(KP,N)=GRAD(KP,N)+HR(K,I,N)*TN(NG)
                  U=0.D0
                  DO 7026 L=1,NPG
                     U=U+UN(NF,N)*FN(I,L)
 7026             CONTINUE
                  UMJ(KP,N)=UMJ(KP,N)+U
70162          CONTINUE
70161       CONTINUE
70160    CONTINUE
         DO 7027 K=KDEB,KFIN
            KP=K-KDEB+1
            UM(KP)=UMJ(KP,1)*UMJ(KP,1)+UMJ(KP,2)*UMJ(KP,2)
            UM(KP)=SQRT(UM(KP))
 7027    CONTINUE
      ENDIF
C                 ***********
C                 *   2 D   *
C                 ***********
      IF(IDIM.EQ.2)THEN
         QUA4=CHGL
         DO 7004 K=KDEB,KFIN
            KP=K-KDEB+1
            AL2(KP)=1.D0/AL(KP)/AL(KP)
            AH2(KP)=1.D0/AH(KP)/AH(KP)
            XMB(KP)=(AL(KP)+AH(KP))/2.D0
 7004    CONTINUE
         DO 7007 K=KDEB,KFIN
            KP=K-KDEB+1
            UMI(KP,1)=UMI(KP,1)/AIRE(KP)
            UMI(KP,2)=UMI(KP,2)/AIRE(KP)
            AX=GRAD(KP,1)*GRAD(KP,1)+GRAD(KP,2)*GRAD(KP,2)
            AX=SQRT(AX)+xpeti2
            GRAD(KP,1)=GRAD(KP,1)/AX
            GRAD(KP,2)=GRAD(KP,2)/AX
            UP(KP)=(GRAD(KP,1)*UMJ(KP,1)+GRAD(KP,2)*UMJ(KP,2))
            UPI(KP,1)=UP(KP)*GRAD(KP,1)
            UPI(KP,2)=UP(KP)*GRAD(KP,2)
 7007    CONTINUE
         DO 70074 K=KDEB,KFIN
            KP=K-KDEB+1
            BMX=UMJ(KP,1)/AL(KP)
            BMY=UMJ(KP,2)/AH(KP)
            BM(KP)=BMX*BMX+BMY*BMY
            BM(KP)=SQRT(BM(KP))+xpeti2
            BPX=UPI(KP,1)/AL(KP)
            BPY=UPI(KP,2)/AH(KP)
            BP(KP)=BPX*BPX+BPY*BPY
            BP(KP)=SQRT(BP(KP))+xpeti2
70074    CONTINUE
C     
C     DECENTREMENT
C     
         IF(IDCEN.EQ.1)THEN
C CENTREE
            DO 70081 K=KDEB,KFIN
               KP=K-KDEB+1
               TO1(KP)=0.D0
               TO2(KP)=0.D0
               CXT(KP)=0.D0
               CYT(KP)=0.D0
               CXY(KP)=0.D0
70081       CONTINUE
         ELSEIF(IDCEN.EQ.2)THEN
C SUPGDC
            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*UP(KP)/BP(KP)
               ALFA=UP(KP)*HPK/COEFT(KP)/2.D0
               AKSI=ALFA/3.D0
               IF(ALFA.GT.3.D0)AKSI=1.D0
               CCP=AKSI/BP(KP)
               CPT=CCP-CCT
               CC2=0.D0
               IF(CPT.GE.0.D0)CC2=CPT
               CC2=CC2*CC0
               TO1(KP)=CCT
               TO2(KP)=CC2
               CXT(KP)=(UMJ(KP,1)*UMJ(KP,1)*CCT+UPI(KP,1)*UPI(KP,1)*CC2)
               CYT(KP)=(UMJ(KP,2)*UMJ(KP,2)*CCT+UPI(KP,2)*UPI(KP,2)*CC2)
               CXY(KP)=(UMJ(KP,1)*UMJ(KP,2)*CCT+UPI(KP,1)*UPI(KP,2)*CC2)
 7008       CONTINUE
         ELSEIF(IDCEN.EQ.3)THEN
C SUPG
            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)
               TO1(KP)=CCT
               TO2(KP)=0.D0
               CXT(KP)=UMJ(KP,1)*UMJ(KP,1)*CCT
               CYT(KP)=UMJ(KP,2)*UMJ(KP,2)*CCT
               CXY(KP)=UMJ(KP,1)*UMJ(KP,2)*CCT
 7018       CONTINUE
         ELSEIF(IDCEN.EQ.4)THEN
C Tenseur Visqueux
            DT19=DTM1*0.5D0
            DO 7009 K=KDEB,KFIN
               KP=K-KDEB+1
               TO1(KP)=DT19
               TO2(KP)=0.D0
               CXT(KP)=UMJ(KP,1)*UMJ(KP,1)*DT19
               CYT(KP)=UMJ(KP,2)*UMJ(KP,2)*DT19
               CXY(KP)=UMJ(KP,1)*UMJ(KP,2)*DT19
 7009       CONTINUE
         ENDIF
C     
C     CALCUL DU PAS DE TEMPS
C     
C*IBMDIR* PREFER SCALAR
         DO 7010 K=KDEB,KFIN
            KP=K-KDEB+1
            DT0=DT
            DT1=2.D0/
     &           (UMJ(KP,1)*UMJ(KP,1)/(CXT(KP)+COEFT(KP))
     &           +UMJ(KP,2)*UMJ(KP,2)/(CYT(KP)+COEFT(KP)))
            DT2=0.5D0/
     &           ((CXT(KP)+COEFT(KP))*AL2(KP)
     &           +(CYT(KP)+COEFT(KP))*AH2(KP))
            IF(DT1.LT.DT)DT=DT1
            IF(DT2.LT.DT)DT=DT2
            IF(DT.NE.DT0) THEN
               DTT1=DT1
               DTT2=DT2
               DIAEL=XMB(KP)
               NUEL=K
            END IF
 7010    CONTINUE
         DO 7011 I=1,NP
            DO 70111 K=KDEB,KFIN
               KP=K-KDEB+1
               WTILDE(KP,I)=DRR(I,K)
     &              +(TO1(KP)*(UMJ(KP,1)*HR(K,I,1)+UMJ(KP,2)*HR(K,I,2))
     &              + TO2(KP)*(UPI(KP,1)*HR(K,I,1)+UPI(KP,2)*HR(K,I,2))
     &              )*AIRE(KP)
70111       CONTINUE
 7011    CONTINUE
         DO 7012 K=KDEB,KFIN
            KP=K-KDEB+1
            CHGLD(KP)=COEFT(KP)/12.D0*QUA4*(AL2(KP)+AH2(KP))
            CHGLP(KP)=(CXT(KP)*AL2(KP)+CYT(KP)*AH2(KP))/12.D0*QUA4
 7012    CONTINUE
C                 ***********
C                 *   3 D   *
C                 ***********
      ELSEIF(IDIM.EQ.3)THEN
         CUB8=CHGL
         DO 8004 K=KDEB,KFIN
            KP=K-KDEB+1
            AL2(KP)=1.D0/AL(KP)/AL(KP)
            AH2(KP)=1.D0/AH(KP)/AH(KP)
            AP2(KP)=1.D0/AP(KP)/AP(KP)
            XMB(KP)=(AL(KP)+AH(KP)+AP(KP))/3.D0
            CFM(KP)=AL(KP)*AH(KP)/AP(KP)+AL(KP)*AP(KP)/AH(KP)+
     &           AP(KP)*AH(KP)/AL(KP)
 8004    CONTINUE
         DO 8007 K=KDEB,KFIN
            KP=K-KDEB+1
            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=GRAD(KP,1)*GRAD(KP,1)+GRAD(KP,2)*GRAD(KP,2)
     &           +GRAD(KP,3)*GRAD(KP,3)
            AX=SQRT(AX)+xpeti2
            GRAD(KP,1)=GRAD(KP,1)/AX
            GRAD(KP,2)=GRAD(KP,2)/AX
            GRAD(KP,3)=GRAD(KP,3)/AX
            UP(KP)=GRAD(KP,1)*UMI(KP,1)+GRAD(KP,2)*UMI(KP,2)
     &           +GRAD(KP,3)*UMI(KP,3)
            UPI(KP,1)=UP(KP)*GRAD(KP,1)
            UPI(KP,2)=UP(KP)*GRAD(KP,2)
            UPI(KP,3)=UP(KP)*GRAD(KP,3)
 8007    CONTINUE
         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))+xpeti2
            BPX=UPI(KP,1)/AL(KP)
            BPY=UPI(KP,2)/AH(KP)
            BPZ=UPI(KP,3)/AP(KP)
            BP(KP)=BPX*BPX+BPY*BPY+BPZ*BPZ
            BP(KP)=SQRT(BP(KP))+xpeti2
80074    CONTINUE
C     
C     DECENTREMENT
C     
         IF(IDCEN.EQ.1)THEN
C CENTREE
            DO 80081 K=KDEB,KFIN
               KP=K-KDEB+1
               TO1(KP)=0.D0
               TO2(KP)=0.D0
               CXT(KP)=0.D0
               CYT(KP)=0.D0
               CZT(KP)=0.D0
               CXY(KP)=0.D0
               CXZ(KP)=0.D0
               CYZ(KP)=0.D0
80081       CONTINUE
         ELSEIF(IDCEN.EQ.2)THEN
C SUPGDC
            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*UP(KP)/BP(KP)
               ALFA=UP(KP)*HPK/COEFT(KP)/2.D0
               AKSI=ALFA/3.D0
               IF(ALFA.GT.3.D0)AKSI=1.D0
               CCP=AKSI/BP(KP)
               CPT=CCP-CCT
               CC2=0.D0
               IF(CPT.GE.0.D0)CC2=CPT
               CC2=CC2*CC0
               TO1(KP)=CCT
               TO2(KP)=CC2
               CXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT+UPI(KP,1)*UPI(KP,1)*CC2)
               CYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT+UPI(KP,2)*UPI(KP,2)*CC2)
               CZT(KP)=(UMI(KP,3)*UMI(KP,3)*CCT+UPI(KP,3)*UPI(KP,3)*CC2)
               CXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT+UPI(KP,1)*UPI(KP,2)*CC2)
               CXZ(KP)=(UMI(KP,1)*UMI(KP,3)*CCT+UPI(KP,1)*UPI(KP,3)*CC2)
               CYZ(KP)=(UMI(KP,2)*UMI(KP,3)*CCT+UPI(KP,2)*UPI(KP,3)*CC2)
 8008       CONTINUE
         ELSEIF(IDCEN.EQ.3)THEN
C SUPG
            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)
               TO1(KP)=CCT
               TO2(KP)=0.D0
               CXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT
               CYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT
               CZT(KP)=UMI(KP,3)*UMI(KP,3)*CCT
               CXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT
               CXZ(KP)=UMI(KP,1)*UMI(KP,3)*CCT
               CYZ(KP)=UMI(KP,2)*UMI(KP,3)*CCT
 8018       CONTINUE
         ELSEIF(IDCEN.EQ.4)THEN
C Tenseur Visqueux
            DT19=DTM1*0.5D0
            DO 8009 K=KDEB,KFIN
               KP=K-KDEB+1
               TO1(KP)=DT19
               TO2(KP)=0.D0
               CXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19
               CYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19
               CZT(KP)=UMI(KP,3)*UMI(KP,3)*DT19
               CXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19
               CXZ(KP)=UMI(KP,1)*UMI(KP,3)*DT19
               CYZ(KP)=UMI(KP,2)*UMI(KP,3)*DT19
 8009       CONTINUE
         ENDIF
C     
C     CALCUL DU PAS DE TEMPS
C     
C*IBMDIR* PREFER SCALAR
         DO 8010 K=KDEB,KFIN
            KP=K-KDEB+1
            DT0=DT
            DT1=2.D0/
     &           (UMI(KP,1)*UMI(KP,1)/(CXT(KP)+COEFT(KP))
     &           +UMI(KP,2)*UMI(KP,2)/(CYT(KP)+COEFT(KP))
     &           +UMI(KP,3)*UMI(KP,3)/(CZT(KP)+COEFT(KP)))
            DT2=0.5D0/
     &           ((CXT(KP)+COEFT(KP))*AL2(KP)
     &           +(CYT(KP)+COEFT(KP))*AH2(KP)
     &           +(CZT(KP)+COEFT(KP))*AP2(KP))
            IF(DT1.LT.DT)DT=DT1
            IF(DT2.LT.DT)DT=DT2
            IF(DT.NE.DT0) THEN
               DTT1=DT1
               DTT2=DT2
               DIAEL=XMB(KP)
               NUEL=K
            END IF
 8010    CONTINUE
         DO 8011 I=1,NP
            DO 80111 K=KDEB,KFIN
               KP=K-KDEB+1
               WTILDE(KP,I)=DRR(I,K)
     &              +(TO1(KP)*(UMI(KP,1)*HR(K,I,1)+UMI(KP,2)*HR(K,I,2)
     &              +UMI(KP,3)*HR(K,I,3))
     &              + TO2(KP)*(UPI(KP,1)*HR(K,I,1)+UPI(KP,2)*HR(K,I,2)
     &              +UPI(KP,3)*HR(K,I,3))
     &              )*AIRE(KP)
80111       CONTINUE
 8011    CONTINUE
         DO 8012 K=KDEB,KFIN
            KP=K-KDEB+1
            GEO1=CFM(KP)*CUB8
            CHGLD(KP)=COEFT(KP)*GEO1
            CHGLP(KP)=(CXT(KP)+CYT(KP)+CZT(KP))/3.D0*GEO1
 8012    CONTINUE
      ENDIF
      RETURN
 1002 FORMAT(10(1X,1PE11.4))
      END