C SMTP      SOURCE    CHAT      11/03/16    21:31:37     6902           
      SUBROUTINE SMTP
C-----------------------------------------------------------------------
C 1) Calcul du second membre du systeme en trace de charge dans le cas
C de la résolution des équations de Darcy par EFMH avec le modèle DARCY.
C 2) Calcul du second membre du systeme en trace de charge et/ou d'une
C contribution matricielle provenant de la convection dans le cas de la
C résolution d'une équation de diffusion-convection par EFMH avec le
C modèle DARCY.
C-----------------------------------------------------------------------
C
C---------------------------
C Phrase d'appel (GIBIANE) :
C---------------------------
C
C  (RIG2) (CHP3) = 'SMTP' MMODEL TABLE1 RIG1 (TABLE2) (CHP1) (CHP2) ;
C
C------------------------
C Operandes et resultat :
C------------------------
C
C  MMODEL  : Objet modele décrivant la formulation.
C  TABLE1  : TABLE DOMAINE contenant les maillages et les connectivités.
C  RIG1    : Matrices hybrides elementaires de DARCY crees par MHYB.
C  TABLE2  : TABLE DARCY_TRANSITOIRE contenant les infos transitoires.
C  CHP1    : CHPO centre de composante SOUR, source moyenne par élément.
C  CHP2    : CHPO face de composante FLUX, flux de vitesse convective.
C  CHP3    : CHPO face de composante FLUX, second membre du pb en TH.
C  RIG2    : RIGIDITE, contribution du terme convectif au LHS en TH.
C
C-----------------------------------------------------------------------
C
C Langage : ESOPE + FORTRAN77
C
C Auteurs : 09/93 F.DABBENE   - Cas permanent
C           09/94 X.NOUVELLON - Extension au transitoire
C           12/95 F.DABBENE   - Extension à la convection
C
C-----------------------------------------------------------------------
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
*
-INC CCOPTIO                         
-INC SMCHAML
-INC SMCHPOI
-INC SMELEME
-INC SMMODEL
-INC SMRIGID
-INC SMTABLE
*                                                                 
      SEGMENT IPMAHY
         INTEGER MAHYBR(NSOUS)
      ENDSEGMENT
*       
      LOGICAL LOGRE,LOGIN
      CHARACTER*4 NOMTOT(1)
      CHARACTER*6 CHAR6
      CHARACTER*8 TAPIND,TYPOBJ,CHARIN,CHARRE,LETYPE
*                                                 
* Initialisations
*
      ISOU1  = 0
      IHN1   = 0
      ITPN1  = 0
      IVALIN = 0
      XVALIN = 0.D0
      LOGIN  = .TRUE.
      IOBIN  = 0
      TAPIND = 'MOT     '
*                                                 
* Lecture du MMODEL
*
      CALL LIROBJ('MMODEL',IPMODE,1,IRET)
      IF (IERR.NE.0) RETURN                     
      MMODEL = IPMODE
*                                                
* Lecture de la TABLE domaine
*                                                        
      IPTABL = 0
      CALL LEKMOD(IPMODE,IPTABL,IRET)
      IF (IERR.NE.0) RETURN
      CHARIN = 'MAILLAGE'
      TYPOBJ = 'MAILLAGE'
      CALL LEKTAB(IPTABL,CHARIN,IOBRE)
      IF (IERR.NE.0) RETURN
      IPGEOM = IOBRE
      CALL LEKTAB(IPTABL,'ELTFA',IOBRE)
      IF (IERR.NE.0) RETURN
      IELTFA = IOBRE
      CALL LEKTAB(IPTABL,'CENTRE',IOBRE)
      IF (IERR.NE.0) RETURN
      ICENTR = IOBRE
      CALL LEKTAB(IPTABL,'FACE',IOBRE)
      IF (IERR.NE.0) RETURN
      IPFACE = IOBRE
*                                                                      
* Lecture de RIGIDITE
*                                                                     
      CALL LIROBJ('RIGIDITE',IPRIGI,1,IRET)
      IF (IERR.NE.0) RETURN                                           
      MRIGID = IPRIGI
*
* Lecture éventuelle de la TABLE Darcy_transitoire
*
      IPTAB1 = 0
      CALL LIRTAB('DARCY_TRANSITOIRE',IPTAB1,0,IRET)
      IF (IERR.NE.0) RETURN
      IF (IRET.EQ.1) THEN
         THEMIN = -1.D-12
         THEMAX = 1.D0 - THEMIN
         TYPOBJ = 'FLOTTANT'
         CALL ACCTAB(IPTAB1,TAPIND,IVALIN,XVALIN,'THETA',LOGIN,IOBIN,
     S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
         IF (IERR.NE.0) RETURN 
         THETA  = XVALRE
         IF ((THETA.LT.THEMIN).OR.(THETA.GT.THEMAX)) THEN
            REAERR(1)   = REAL(THETA)
C            REAERR(2)   = REAL(0.D0)
C            REAERR(3)   = REAL(1.D0)
            MOTERR(1:8) = '  THETA '
            CALL ERREUR(42)
            RETURN
         ENDIF
         IF (THETA.LT.0.D0) THETA=0.D0
         IF (THETA.GT.1.D0) THETA=1.D0
         TYPOBJ = '      '
         CALL ACCTAB(IPTAB1,TAPIND,IVALIN,XVALIN,'THETA_CONVECTION',
     S                                              LOGIN,IOBIN,
     S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
         IF (IERR.NE.0) RETURN 
         IF (TYPOBJ.EQ.'FLOTTANT') THEN
           THETAC  = XVALRE
           IF ((THETAC.LT.THEMIN).OR.(THETAC.GT.THEMAX)) THEN
              REAERR(1)   = REAL(THETAC)
              REAERR(2)   = REAL(0.D0)
              REAERR(3)   = REAL(1.D0)
              MOTERR(1:8) = 'TETACONV'
              CALL ERREUR(42)
              RETURN
           ENDIF
           IF (THETAC.LT.0.D0) THETA=0.D0
           IF (THETAC.GT.1.D0) THETA=1.D0
         ELSE
           THETAC  = THETA
         ENDIF
         CALL ACCTAB(IPTAB1,TAPIND,IVALIN,XVALIN,'PAS',LOGIN,IOBIN,
     S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
         IF (IERR.NE.0) RETURN 
         DELTAT = XVALRE
         IF (DELTAT.LE.0.D0) THEN
            REAERR(1)   = REAL(DELTAT)
            REAERR(2)   = REAL(0.D0)
            MOTERR(1:8) = ' DELTAT '
            CALL ERREUR(41)
            RETURN
         ENDIF
         TYPOBJ = 'MCHAML  '
         CALL ACCTAB(IPTAB1,TAPIND,IVALIN,XVALIN,'SURF',LOGIN,IOBIN,
     S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
         IF (IERR.NE.0) RETURN 
         IPCK = IOBRE
         TYPOBJ = 'CHPOINT '
         CALL ACCTAB(IPTAB1,TAPIND,IVALIN,XVALIN,'TRACE',LOGIN,IOBIN,
     S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
         IF (IERR.NE.0) RETURN 
         ITP = IOBRE
         CALL ACCTAB(IPTAB1,TAPIND,IVALIN,XVALIN,'CHARGE',LOGIN,IOBIN,
     S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
         IF (IERR.NE.0) RETURN 
         IH = IOBRE
      ELSE
         THETA  = 0.D0
         DELTAT = 0.D0
         IPCK   = 0
         ITP    = 0
         IH     = 0
      ENDIF                                                           
*
* Lecture des champoints facultatifs SOUR et FLUX
*     prise en compte des termes sources
*     prise en compte de la convection
*
      ISOUR = 0
      CALL LIROBJ('CHPOINT',ISOUR,0,IRET)
      IF (IERR.NE.0) RETURN
*
      IFLUX = 0
      CALL LIROBJ('CHPOINT',IFLUX,0,IRET)
      IF (IERR.NE.0) RETURN
*
* Tri des deux champoints facultatifs
*
      IFAC = ISOUR + IFLUX
      IF (IFAC.EQ.0.AND.IPTAB1.EQ.0) THEN
         MOTERR(1:8) = 'CHAMPOIN'
         CALL ERREUR(37)
         RETURN
      ENDIF
      IF (IFAC.NE.0) THEN
         NBCOMP = 0
         INDIC  = 0
         CALL QUEPOI(ISOUR,ICENTR,INDIC,NBCOMP,NOMTOT)
         IF (NOMTOT(1).NE.'SOUR'.AND.NOMTOT(1).NE.'FLUX') THEN
            MOTERR(1:4) = NOMTOT(1)
            CALL ERREUR(197)
            RETURN
         ELSE IF (NOMTOT(1).EQ.'FLUX') THEN
            IFAC  = IFLUX
            IFLUX = ISOUR
            ISOUR = IFAC
         ENDIF
      ENDIF
*
* Dans le cas de la convection, récup. MCHAML orientation normale
*
      IF (IFLUX.NE.0) THEN
         CALL LEKTAB(IPTABL,'XXNORMAE',IORIE)
         IF (IERR.NE.0) RETURN
      ENDIF
*                                                                     
* Test du CHAMPOINT de composante SOUR de spg CENTRE
*                                                                   
      IF (ISOUR.NE.0) THEN
         NBCOMP    = 1
         INDIC     = 1
         NOMTOT(1) = 'SOUR'
         CALL QUEPOI(ISOUR,ICENTR,INDIC,NBCOMP,NOMTOT)
         IF (IERR.NE.0) RETURN
      ENDIF
*                                                                     
* Test du CHAMPOINT de composante FLUX de spg FACE
*                                                                   
      IF (IFLUX.NE.0) THEN
         NBCOMP    = 1
         INDIC     = 1
         NOMTOT(1) = 'FLUX'
         CALL QUEPOI(IFLUX,IPFACE,INDIC,NBCOMP,NOMTOT)
         IF (IERR.NE.0) RETURN
      ENDIF
*                                                                     
* Test du CHAMPOINT de composante TH de spg FACES
*                                                                   
      IF (ITP.NE.0) THEN
         NBCOMP    = 1
         INDIC     = 1
         NOMTOT(1) = 'TH'
         CALL QUEPOI(ITP,IPFACE,INDIC,NBCOMP,NOMTOT)
         IF (IERR.NE.0) RETURN
      ENDIF
*                                                                     
* Test du CHAMPOINT de composante H de spg CENTRE
*                                                                   
      IF (IH.NE.0) THEN
         NBCOMP = 1
         INDIC  = 1
         NOMTOT(1) = 'H'
         CALL QUEPOI(IH,ICENTR,INDIC,NBCOMP,NOMTOT)
         IF (IERR.NE.0) RETURN
      ENDIF
*                                                 
* Test de la formulation
*
      SEGACT MMODEL
      NSOUS  = KMODEL(/1)
      SEGINI IPMAHY
      IDARCY = 0
      IPT1=IPGEOM
      IPT2=IPGEOM
      DO 10 ISOUS=1,NSOUS
         IF(NSOUS.GT.1)THEN
              SEGACT IPT2
              IPT1=IPT2.LISOUS(ISOUS)
              SEGDES IPT2
         ENDIF
         IMODEL = KMODEL(ISOUS)
         SEGACT IMODEL
         LETYPE = FORMOD(1)
         IF (LETYPE.EQ.'DARCY') THEN
            IDARCY = IDARCY + 1
            MAHYBR(ISOUS) = IPT1
         ENDIF
         SEGDES IMODEL
 10   CONTINUE
      SEGDES MMODEL
      IF (IDARCY.EQ.0) THEN
         MOTERR = LETYPE
         CALL ERREUR(193)
         GOTO 100
      ENDIF                                         
*
* Recuperation des pointeurs ELTFA pour les zones ou DARCY est defini
*
      MELEME = IELTFA
      SEGACT MELEME
      LZONES = LISOUS(/1)
      IF (LZONES.EQ.0) LZONES = 1
      IPT1   = IPGEOM
      SEGACT IPT1
      DO 30 ISOUS=1,NSOUS
         IMAMEL = MAHYBR(ISOUS)
         IF (IMAMEL.NE.0) THEN
            DO 20 ISZ=1,LZONES
               IF (LZONES.EQ.1) THEN
                  IPT2 = IPT1
                  IPT3 = MELEME
               ELSE
                  IPT2 = IPT1.LISOUS(ISZ)
                  IPT3 = LISOUS(ISZ)
               ENDIF
               IF (IPT2.EQ.IMAMEL) THEN
                  MAHYBR(ISOUS) = IPT3
                  GOTO 30
               ENDIF
 20         CONTINUE
            IF (IMAMEL.EQ.MAHYBR(ISOUS)) THEN
               INTERR(1) = ISOUS
               CALL ERREUR(698)
               SEGDES IPT1
               SEGDES MELEME
               GOTO 100
            ENDIF
         ENDIF
 30   CONTINUE
      SEGDES IPT1
      SEGDES MELEME
*                                                                      
* Test du sous-type de la matrice de rigiditée récupérée
*                                                                     
      SEGACT MRIGID
      LETYPE = MTYMAT
      IF (LETYPE.NE.'DARCY') THEN
         MOTERR(1:8)  = 'RIGIDITE'
         MOTERR(9:16) = 'DARCY   '
         CALL ERREUR(79)
         SEGDES MRIGID
         GOTO 100
      ENDIF
*                                                                      
* Controle des pointeurs de MELEME de la rigidité
*                                                                     
      DO 40 ISOUS=1,NSOUS
         IMAMEL = MAHYBR(ISOUS)
         IF (IMAMEL.NE.0) THEN
            IPTEST = IRIGEL(1,ISOUS)
            IF (MAHYBR(ISOUS).NE.IRIGEL(1,ISOUS)) THEN
               MOTERR(1:8)  = 'DARCY   '
               MOTERR(9:16) = 'ELTFA   '
               INTERR(1)    = ISOUS
               CALL ERREUR(698)
               SEGDES MRIGID
               GOTO 100
            ENDIF
         ENDIF
 40   CONTINUE
      SEGDES MRIGID
*
* Vérification du support du MCHAML % au MMODEL
* Controle du nombre de composantes reelles
*
      IF (IPCK.NE.0) THEN
         ISUP  = 2
         ICOND = 1
         CALL QUESUP(IPMODE,IPCK,ISUP,ICOND,IRET,IRET2)
         IF (IRET.GT.1) GOTO 100
         MCHELM = IPCK
         SEGACT MCHELM
         DO 50 ISOUS=1,NSOUS
            IF (MAHYBR(ISOUS).NE.0) THEN
               MCHAML = ICHAML(ISOUS)
               SEGACT MCHAML
               N2     = TYPCHE(/2)
               IF (N2.GT.1) THEN
                  CALL ERREUR(320)
                  SEGDES MCHELM
                  SEGDES MCHAML
                  GOTO 100
               ENDIF
               CHAR6  = TYPCHE(1)(1:6)
               IF (CHAR6.NE.'REAL*8') THEN
                  MOTERR(1:8) = NOMCHE(1)
                  CALL ERREUR(679)
                  SEGDES MCHELM
                  SEGDES MCHAML
                  GOTO 100
               ENDIF
               SEGDES MCHAML
            ENDIF
 50      CONTINUE
         SEGDES MCHELM
      ENDIF
*                                                                    
* Construction de MCHPOI
*                                                                    
      SEGDES IPMAHY
      CALL SMTP1(IPMAHY,IPFACE,IPRIGI,THETA,THETAC,DELTAT,IPCK,
     S           ITP,IH,ISOUR,IFLUX,IORIE,MCHPOI,IRI1)
*
      IF (MCHPOI.NE.0) CALL ECROBJ('CHPOINT',MCHPOI)
      IF (IRI1.NE.0)   CALL ECROBJ('RIGIDITE',IRI1)
*
* Ménage
*
 100  CONTINUE
      SEGSUP IPMAHY
      RETURN                                                          
      END