C CCGMUS    SOURCE    GOUNAND   06/12/08    21:15:18     5597           
      SUBROUTINE CCGMUS(LCOF,NOMLOI,
     $     FC,
     $     IMPR,IRET)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C***********************************************************************
C NOM         : CCGMUS
C DESCRIPTION : Calcul de la loi de comportement aux points de Gauss : 
C               Taille d'un élément suivant une direction donnée
C     
C LANGAGE     : ESOPE
C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
C               mél : gounand@semt2.smts.cea.fr
C***********************************************************************
C APPELES          : 
C APPELE PAR       : 
C***********************************************************************
C ENTREES            :
C ENTREES/SORTIES    :
C SORTIES            : -
C TRAVAIL            :
C***********************************************************************
C VERSION    : v1, 30/08/06, version initiale
C HISTORIQUE : v1, 30/08/06, création
C HISTORIQUE :
C HISTORIQUE :
C***********************************************************************
C Prière de PRENDRE LE TEMPS de compléter les commentaires
C en cas de modification de ce sous-programme afin de faciliter
C la maintenance !
C***********************************************************************
-INC CCOPTIO
-INC CCREEL
CBEGININCLUDE SMCHAEL
      SEGMENT MCHAEL
      POINTEUR IMACHE(N1).MELEME
      POINTEUR ICHEVA(N1).MCHEVA
      ENDSEGMENT
      SEGMENT MCHEVA
      REAL*8 VELCHE(NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM)
      ENDSEGMENT
      SEGMENT LCHEVA
      POINTEUR LISCHE(NBCHE).MCHEVA
      ENDSEGMENT
CENDINCLUDE SMCHAEL
      INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM,N1
      POINTEUR FC.MCHEVA
      POINTEUR LCOF.LCHEVA
      POINTEUR JMAJAC.MCHEVA
      POINTEUR JMIJAC.MCHEVA
      POINTEUR JDTJAC.MCHEVA
      POINTEUR JMAREG.MCHEVA
      POINTEUR JDIAMA.MCHEVA
      POINTEUR JPC.MCHEVA
      POINTEUR JVIT.MCHEVA
      POINTEUR JRHO.MCHEVA
      POINTEUR JMU.MCHEVA
      POINTEUR JPEC.MCHEVA
      CHARACTER*8 NOMLOI
      INTEGER ICOF
*
-INC TMXMAT      
      POINTEUR A.MXMAT
      POINTEUR JMA.MXMAT
      POINTEUR KJMA.MXMAT
      POINTEUR JM1.MXMAT
      POINTEUR K.MXMAT
*
      SEGMENT MVIT
      POINTEUR MVCOMP(IDIM).MCHEVA
      ENDSEGMENT
      LOGICAL LCRIT
*
      INTEGER IMPR,IRET
*     
* Executable statements
*     
      IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ccgmus'
      XPETI=SQRT(XPETIT)
      CALL CH2INT(NOMLOI(7:7),IMETH,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      CALL CH2INT(NOMLOI(8:8),IDIMI,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      NLFC=FC.VELCHE(/6)
      NPFC=FC.VELCHE(/5)
      ICOF=0 
*
* Récupération de RHO et MU
*
      ICOF=ICOF+1
      JRHO=LCOF.LISCHE(ICOF)
      NLRO=JRHO.VELCHE(/6)
      NPRO=JRHO.VELCHE(/5)
      ICOF=ICOF+1
      JMU=LCOF.LISCHE(ICOF)
      NLMU=JMU.VELCHE(/6)
      NPMU=JMU.VELCHE(/5)
*
*     Récupération des coefficients de la metrique
*
      ICOF1=ICOF+1
      SEGINI MVIT
      DO IIDIM=1,IDIM
         ICOF=ICOF+1
         JVIT=LCOF.LISCHE(ICOF)
         IF (ICOF.EQ.ICOF1) THEN 
            NLJV=JVIT.VELCHE(/6)
            NPJV=JVIT.VELCHE(/5)
         ELSE
            NLJV2=JVIT.VELCHE(/6)
            NPJV2=JVIT.VELCHE(/5)
            IF (NLJV2.NE.NLJV.OR.NPJV2.NE.NPJV) THEN 
               WRITE(IOIMP,*) 'Erreur grave dims JVIT'
               GOTO 9999
            ENDIF
         ENDIF
         MVIT.MVCOMP(IIDIM)=JVIT
      ENDDO
*
* Récupération du Péclet critique
*
      ICOF=ICOF+1
      JPEC=LCOF.LISCHE(ICOF)
      NLPEC=JPEC.VELCHE(/6)
      NPPEC=JPEC.VELCHE(/5)
*
      ICOF=ICOF+1
      JMAJAC=LCOF.LISCHE(ICOF)
C      NLJA=JMAJAC.VELCHE(/6)
C      NPJA=JMAJAC.VELCHE(/5)
      IREF=JMAJAC.VELCHE(/4)
      IREL=JMAJAC.VELCHE(/3)
*
      IF (IREL.NE.IDIM.OR.IREF.NE.IDIM) THEN
         WRITE(IOIMP,*) 'Erreur dims JMAJAC'
         GOTO 9999
      ENDIF
*
      ICOF=ICOF+1
      JMIJAC=LCOF.LISCHE(ICOF)
      IF (JMIJAC.EQ.0) THEN
         WRITE(IOIMP,*) 'Erreur JMIJAC=0'
         GOTO 9999
      ENDIF
      NLJI=JMIJAC.VELCHE(/6)
      NPJI=JMIJAC.VELCHE(/5)
      IREL2=JMIJAC.VELCHE(/4)
      IREF2=JMIJAC.VELCHE(/3)
*
      IF (IREL2.NE.IREL.OR.IREF2.NE.IREF) THEN
         WRITE(IOIMP,*) 'Erreur dims JMIJAC'
         GOTO 9999
      ENDIF
*
      ICOF=ICOF+1
      ICOF=ICOF+1
      JMAREG=LCOF.LISCHE(ICOF)
      NLJR=JMAREG.VELCHE(/6)
      NPJR=JMAREG.VELCHE(/5)
      I1  =JMAREG.VELCHE(/4)
      I2  =JMAREG.VELCHE(/3)
      IF ((NLJR.NE.1).OR.(NPJR.NE.1).OR.(I1.NE.IREF).OR.(I2.NE.IREF))
     $     THEN
         WRITE(IOIMP,*) 'Erreur dims JMAREG'
         GOTO 9999
      ENDIF
      ICOF=ICOF+1
      JDIAMA=LCOF.LISCHE(ICOF)
      NLJD=JDIAMA.VELCHE(/6)
      NPJD=JDIAMA.VELCHE(/5)
      I1  =JDIAMA.VELCHE(/4)
      I2  =JDIAMA.VELCHE(/3)
      IF ((NLJD.NE.1).OR.(NPJD.NE.1).OR.(I1.NE.1).OR.(I2.NE.1))
     $     THEN
         WRITE(IOIMP,*) 'Erreur dims JDIAMA'
         GOTO 9999
      ENDIF
      XDIAMA=JDIAMA.VELCHE(1,1,1,1,1,1)
*
* Objets temporaires 
*
      LDIM1=IREL
      LDIM2=1
      SEGINI,A
      LDIM1=IREF
      LDIM2=1
      SEGINI,JMA
      LDIM1=IREF
      LDIM2=1
      SEGINI,KJMA
      LDIM1=IREF
      LDIM2=IREL
      SEGINI,JM1
      LDIM1=IREF
      LDIM2=IREF
      SEGINI,K
*
* Calcul de la métrique des éléments réguliers
*
      CALL MAMA(JMAREG.VELCHE,IREF,IREF,
     $     'COPIE   ',K.XMAT,IREF,IREF,
     $     IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
*      SEGPRT,H
*      
      DO ILFC=1,NLFC
         IF (NLRO.EQ.1) THEN
            ILRO=1
         ELSE
            ILRO=ILFC
         ENDIF
         IF (NLMU.EQ.1) THEN
            ILMU=1
         ELSE
            ILMU=ILFC
         ENDIF
         IF (NLJV.EQ.1) THEN
            ILJV=1
         ELSE
            ILJV=ILFC
         ENDIF
         IF (NLPEC.EQ.1) THEN
            ILPEC=1
         ELSE
            ILPEC=ILFC
         ENDIF
C         IF (NLJA.EQ.1) THEN
C            ILJA=1
C         ELSE
C            ILJA=ILFC
C         ENDIF
         IF (NLJI.EQ.1) THEN
            ILJI=1
         ELSE
            ILJI=ILFC
         ENDIF
*
         DO IPFC=1,NPFC
            IF (NPRO.EQ.1) THEN
               IPRO=1
            ELSE
               IPRO=IPFC
            ENDIF
            IF (NPMU.EQ.1) THEN
               IPMU=1
            ELSE
               IPMU=IPFC
            ENDIF
            IF (NPJV.EQ.1) THEN
               IPJV=1
            ELSE
               IPJV=IPFC
            ENDIF
            IF (NPPEC.EQ.1) THEN
               IPPEC=1
            ELSE
               IPPEC=IPFC
            ENDIF
C            IF (NPJA.EQ.1) THEN
C               IPJA=1
C            ELSE
C               IPJA=IPFC
C            ENDIF
            IF (NPJI.EQ.1) THEN
               IPJI=1
            ELSE
               IPJI=IPFC
            ENDIF
            XRHO=JRHO.VELCHE(1,1,1,1,IPRO,ILRO)
            XMU=JMU.VELCHE(1,1,1,1,IPMU,ILMU)
            XPEC=JPEC.VELCHE(1,1,1,1,IPPEC,ILPEC)
*
* Copie de la vitesse
*
            DO IIDIM=1,IDIM
               JVIT=MVIT.MVCOMP(IIDIM)
               A.XMAT(IIDIM,1)=JVIT.VELCHE(1,1,1,1,IPJV,ILJV)
            ENDDO
*            SEGPRT,A
*
* Copie de l'inverse du jacobien ref->reel
*
            CALL MAMA(JMIJAC.VELCHE(1,1,1,1,IPJI,ILJI),IREF,IREL,
     $           'COPIE   ',
     $           JM1.XMAT,IREF,IREL,
     $           IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
C            SEGPRT,JM1
*
* Vecteur A dans le repère de référence
*
            CALL MAMAMA(JM1.XMAT,IREF,IREL,A.XMAT,IREL,1,
     $           'FOIS    ',JMA.XMAT,IREF,1,
     $           IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
*
* Vecteur A dans le repère de l'élément régulier
*
            CALL MAMAMA(K.XMAT,IREF,IREF,JMA.XMAT,IREF,1,
     $           'FOIS    ',KJMA.XMAT,IREF,1,
     $           IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
*
* Normes de A
*
            XAREG2=XZERO
            DO IIREF=1,IREF
               XAREG2=XAREG2+(KJMA.XMAT(IIREF,1)**2)
            ENDDO
            XA2=XZERO
            DO IIDIM=1,IDIM
               XA2=XA2+(A.XMAT(IIDIM,1)**2)
            ENDDO
*
            IF (XA2.LT.XPETI) THEN
               CONTRI=XZERO
            ELSE
               XV=SQRT(XA2)
* Pondération par un diamètre externe !
*               XH=SQRT(XA2/XAREG2)
               XH=SQRT(XA2/XAREG2)*XDIAMA
*
               IF (IMETH.EQ.1) THEN
* Diffusion artificielle UPWIND
                  CONTRI=((XRHO*XV*XH)/XPEC)
                  CONTRI=SQRT(CONTRI)
                  IF (IDIMI.GT.0) THEN
                     CONTRI=CONTRI*
     $                    (A.XMAT(IDIMI,1)/XV)
                  ENDIF
               ELSEIF (IMETH.EQ.2) THEN
* Diffusion artificielle comme Alberto (equiv critical approx)
* Doit-on décentrer ?
                  LCRIT=((XRHO*XV*XH).GT.(XPEC*XMU))
                  IF (LCRIT) THEN
                     CONTRI=((XRHO*XV*XH)/XPEC)-XMU
                     CONTRI=SQRT(CONTRI)
                     IF (IDIMI.GT.0) THEN
                        CONTRI=CONTRI*
     $                       (A.XMAT(IDIMI,1)/XV)
                     ENDIF
                  ELSE 
                     CONTRI=XZERO
                  ENDIF
               ELSEIF (IMETH.EQ.3) THEN
* Diffusion artificielle (equiv doubly asymptotic)
                  LCRIT=((XRHO*XV*XH).GT.(XPEC*XMU*3.D0))
                  IF (LCRIT) THEN
                     CONTRI=((XRHO*XV*XH)/XPEC)
                  ELSE
                     CONTRI=(((XRHO*XV*XH)/XPEC)**2)/(3.D0*XMU)
                  ENDIF
                  CONTRI=SQRT(CONTRI)
                  IF (IDIMI.GT.0) THEN
                     CONTRI=CONTRI*
     $                    (A.XMAT(IDIMI,1)/XV)
                  ENDIF
               ELSE
                  WRITE(IOIMP,*) 'IMETH=',IMETH,' unknown'
                  GOTO 9999
               ENDIF
            ENDIF
            FC.VELCHE(1,1,1,1,IPFC,ILFC)=
     $           FC.VELCHE(1,1,1,1,IPFC,ILFC)+
     $           CONTRI
         ENDDO
      ENDDO
      SEGSUP,K
      SEGSUP,JM1
      SEGSUP,KJMA
      SEGSUP,JMA
      SEGSUP,A
      SEGSUP,MVIT
*     
* Normal termination
*     
      IRET=0
      RETURN
*     
* Format handling
*     
*     
* Error handling
*     
 9999 CONTINUE
      IRET=1
      WRITE(IOIMP,*) 'An error was detected in subroutine ccgmus'
      RETURN
*     
* End of subroutine CCGMUS
*     
      END