C CCGMUS    SOURCE    GOUNAND   26/01/09    21:15:07     12441          
      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 : 09/06/2015 : SG ajout de tests de positivite pour eviter
C              des pbs dans les SQRTs
C              Ajout du cas ou l'espace de reference est de dimension
C              strictement inferieure a l'espace reel
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 PPARAM
-INC CCOPTIO
-INC CCREEL
-INC TNLIN
*-INC 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 AP.MXMAT
      POINTEUR JMA.MXMAT
      POINTEUR KJMA.MXMAT
      POINTEUR JM1.MXMAT,J.MXMAT
      POINTEUR K.MXMAT
*
      SEGMENT MVIT
      POINTEUR MVCOMP(IDIM).MCHEVA
      ENDSEGMENT
      LOGICAL LCRIT,LFIRST
      LOGICAL LRRHO,LRMU,LRPEC
*
      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.WELCHE(/6)
      NPFC=FC.WELCHE(/5)
      ICOF=0
*
* Récupération de RHO et MU
*
      ICOF=ICOF+1
      JRHO=LCOF.LISCHE(ICOF)
      NLRO=JRHO.WELCHE(/6)
      NPRO=JRHO.WELCHE(/5)
      ICOF=ICOF+1
      JMU=LCOF.LISCHE(ICOF)
      NLMU=JMU.WELCHE(/6)
      NPMU=JMU.WELCHE(/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.WELCHE(/6)
            NPJV=JVIT.WELCHE(/5)
         ELSE
            NLJV2=JVIT.WELCHE(/6)
            NPJV2=JVIT.WELCHE(/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.WELCHE(/6)
      NPPEC=JPEC.WELCHE(/5)
*
* sg : 09/06/2015
* Verification que rho, mu et pec sont positifs.
* Meme si les champs sont positifs, il peut y avoir un souci
* avec les quadratiques aux points de Gauss
* 220 0
*Attention: une des composantes du champ est négative ou nulle
* 220 0
*sa puissance est mise à zéro
      LRRHO=.FALSE.
      DO ILRO=1,NLRO
         DO IPRO=1,NPRO
            XRHO=JRHO.WELCHE(1,1,1,1,IPRO,ILRO)
            IF (XRHO.LT.XZERO) THEN
               IF (.NOT.LRRHO) THEN
                  LRRHO=.TRUE.
                  WRITE(IOIMP,*) 'ccgmus.eso : RHO'
                  CALL ERREUR(220)
               ENDIF
            ENDIF
         ENDDO
      ENDDO
      LRMU=.FALSE.
      DO ILMU=1,NLMU
         DO IPMU=1,NPMU
            XMU=JMU.WELCHE(1,1,1,1,IPMU,ILMU)
            IF (XMU.LT.XZERO) THEN
               IF (.NOT.LRMU) THEN
                  LRMU=.TRUE.
                  WRITE(IOIMP,*) 'ccgmus.eso : MU'
                  CALL ERREUR(220)
               ENDIF
            ENDIF
         ENDDO
      ENDDO
      LRPEC=.FALSE.
      DO ILPEC=1,NLPEC
         DO IPPEC=1,NPPEC
            XPEC=JPEC.WELCHE(1,1,1,1,IPPEC,ILPEC)
            IF (XPEC.LT.XZERO) THEN
               IF (.NOT.LRPEC) THEN
                  LRPEC=.TRUE.
                  WRITE(IOIMP,*) 'ccgmus.eso : PEC'
                  CALL ERREUR(220)
               ENDIF
            ENDIF
         ENDDO
      ENDDO
*dbg                  IF (CONTRI.LT.0.D0) THEN
*dbg                     WRITE(IOIMP,*) 'XRHO=',XRHO,' XMU=',XMU
*dbg     $                    ,' XV=',XV,' XH=',XH,' XPEC=',XPEC
*dbg                     WRITE(IOIMP,*) 'NPRO=',NPRO,'ILRO=',ILRO
*dbg                     do ii=1,npro
*dbg                        XRHO=JRHO.WELCHE(1,1,1,1,Ii,ILRO)
*dbg                        WRITE(IOIMP,*) 'ii XRHO=',ii,XRHO
*dbg                     enddo
*dbg                     GOTO 9999
*dbg                  ENDIF
*
      ICOF=ICOF+1
      JMAJAC=LCOF.LISCHE(ICOF)
C      NLJA=JMAJAC.WELCHE(/6)
C      NPJA=JMAJAC.WELCHE(/5)
      IREF=JMAJAC.WELCHE(/4)
      IREL=JMAJAC.WELCHE(/3)
*dbg      WRITE(IOIMP,*) 'IREF=',IREF,' IREL=',IREL
*dbg      IF (IREL.NE.IDIM.OR.IREF.NE.IDIM) THEN
      IF (IREL.NE.IDIM) THEN
         WRITE(IOIMP,*) 'Erreur dims JMAJAC'
         GOTO 9999
      ENDIF
*
      ICOF=ICOF+1
      JMIJAC=LCOF.LISCHE(ICOF)
      IF (JMIJAC.LE.0) THEN
         write(ioimp,*) 'Jacobien non inversible'
         GOTO 9999
      ENDIF
      NLJI=JMIJAC.WELCHE(/6)
      NPJI=JMIJAC.WELCHE(/5)
      IREL2=JMIJAC.WELCHE(/4)
      IREF2=JMIJAC.WELCHE(/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.WELCHE(/6)
      NPJR=JMAREG.WELCHE(/5)
      I1  =JMAREG.WELCHE(/4)
      I2  =JMAREG.WELCHE(/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.WELCHE(/6)
      NPJD=JDIAMA.WELCHE(/5)
      I1  =JDIAMA.WELCHE(/4)
      I2  =JDIAMA.WELCHE(/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.WELCHE(1,1,1,1,1,1)
*dbg
*dbg      WRITE(IOIMP,*) 'XDIAMA=',XDIAMA
*
* Objets temporaires
*
      LDIM1=IREL
      LDIM2=1
      SEGINI,A
      SEGINI,AP
      LDIM1=IREF
      LDIM2=1
      SEGINI,JMA
      LDIM1=IREF
      LDIM2=1
      SEGINI,KJMA
      LDIM1=IREL
      LDIM2=IREF
      SEGINI,J
      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.WELCHE,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.WELCHE(1,1,1,1,IPRO,ILRO)
            XMU=JMU.WELCHE(1,1,1,1,IPMU,ILMU)
            XPEC=JPEC.WELCHE(1,1,1,1,IPPEC,ILPEC)
            IF (LRRHO) XRHO=MAX(XRHO,XZERO)
            IF (LRMU)  XMU=MAX(XMU,XZERO)
            IF (LRPEC) XPEC=MAX(XPEC,XZERO)

*dbg            LFIRST=(ILFC.EQ.1).AND.(IPFC.EQ.1)
            LFIRST=.FALSE.
            IF (LFIRST) WRITE(IOIMP,*) 'XRHO=',XRHO,' XMU=',XMU,' XPEC='
     $           ,XPEC
*
* Copie de la vitesse
*
            DO IIDIM=1,IDIM
               JVIT=MVIT.MVCOMP(IIDIM)
               A.XMAT(IIDIM,1)=JVIT.WELCHE(1,1,1,1,IPJV,ILJV)
            ENDDO
            IF (LFIRST) WRITE(IOIMP,*) 'A ',(A.XMAT(ii,1),ii=1,IDIM)
*            SEGPRT,A
*
* Copie de l'inverse du jacobien ref->reel
*
            CALL MAMA(JMIJAC.WELCHE(1,1,1,1,IPJI,ILJI),IREF,IREL,
     $           'COPIE   ',
     $           JM1.XMAT,IREF,IREL,
     $           IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
C            SEGPRT,JM1
*
* Copie du jacobien ref->reel
*
            CALL MAMA(JMAJAC.WELCHE(1,1,1,1,IPJI,ILJI),IREL,IREF,
     $           'COPIE   ',
     $           J.XMAT,IREL,IREF,
     $           IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
C            SEGPRT,J
*
* 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 réel (projection sur le domaine considéré)
*
            CALL MAMAMA(J.XMAT,IREL,IREF,JMA.XMAT,IREF,1,
     $           'FOIS    ',AP.XMAT,IREL,1,
     $           IMPR,IRET)
            IF (IRET.NE.0) GOTO 9999
            IF (LFIRST) WRITE(IOIMP,*) 'AP ',(AP.XMAT(ii,1),ii=1,IDIM)
*
* 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
            XAP2=XZERO
            DO IIDIM=1,IDIM
               XAP2=XAP2+(AP.XMAT(IIDIM,1)**2)
            ENDDO
            IF (LFIRST) WRITE(IOIMP,*) 'XAREG2=',XAREG2,' XA2=',XA2,
     $           ' XAP2=',XAP2
*
*old            IF (XA2.LT.XPETI) THEN
            IF (XAP2.LT.XPETI) THEN
               CONTRI=XZERO
            ELSE
*old               XV=SQRT(XA2)
*old* Pondération par un diamètre externe !
*old*               XH=SQRT(XA2/XAREG2)
*old               XH=SQRT(XA2/XAREG2)*XDIAMA
* Pour que le décentrement marche sur une surface
               XV=SQRT(XAP2)
               XH=SQRT(XAP2/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*
*old     $                    (A.XMAT(IDIMI,1)/XV)
     $                    (AP.XMAT(IDIMI,1)/XV)
                  ENDIF
                  IF (LFIRST) WRITE(IOIMP,*) 'CONTRI=',CONTRI
               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*
*old     $                       (A.XMAT(IDIMI,1)/XV)
     $                       (AP.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*
*old     $                    (A.XMAT(IDIMI,1)/XV)
     $                    (AP.XMAT(IDIMI,1)/XV)
                  ENDIF
               ELSE
                  WRITE(IOIMP,*) 'IMETH=',IMETH,' unknown'
                  GOTO 9999
               ENDIF
            ENDIF
            FC.WELCHE(1,1,1,1,IPFC,ILFC)=
     $           FC.WELCHE(1,1,1,1,IPFC,ILFC)+
     $           CONTRI
         ENDDO
      ENDDO
      SEGSUP,K
      SEGSUP,J
      SEGSUP,JM1
      SEGSUP,KJMA
      SEGSUP,JMA
      SEGSUP,AP
      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