C CLMI      SOURCE    MAGN      10/05/18    21:15:13     6675           
      SUBROUTINE CLMI
      PARAMETER (NTB=1)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 ctype,ferm
C***********************************************************************
C
C
C***********************************************************************

-INC CCOPTIO
-INC CCGEOME
-INC SMELEME
-INC SMCOORD
-INC SIZFFB
-INC SMMATRIK
      POINTEUR IMAT1.IMATRI,IMAT2.IMATRI,IMAT3.IMATRI
-INC SMCHPOI
      POINTEUR IZTU1.MPOVAL,MCHPO9.MCHPOI,MSOUP9.MSOUPO,MPOVA9.MPOVAL
      POINTEUR MCHPO8.MCHPOI,MSOUP8.MSOUPO,MPOVA8.MPOVAL
      POINTEUR MCHPOA.MCHPOI,MSOUPA.MSOUPO,MPOVAA.MPOVAL
      POINTEUR MCHPOB.MCHPOI,MSOUPB.MSOUPO,MPOVAB.MPOVAL
      SEGMENT TRAV
      REAL*8 UN(NPTI,IDIM),TN(NPTI,1)
      REAL*8 WT(NP,NPG),WS(NP,NPG)
      REAL*8 UIL(IDIM,NPG),DUIL(IDIM,NPG),ANUK(1),AIRE(1)
      REAL*8 HNM(NPTI),BNM(NPTI),Coef1(NPTI),Coef2(NPTI),DNM(NPTI)
      REAL*8 H32NM(NPTI),CFN(NPTI),CEN(NPTI),DeN(NPTI),D1N(NPTI)
      REAL*8 HHN(NPTI)
      REAL*8 UCL(NPTI),Delta(NPTI),SM1(NPTI),SM2(NPTI)
      REAL*8 FP(NPTI),PHIP(NPTI),VT(NPTI)
      ENDSEGMENT
      
-INC SMLENTI
-INC SMLMOTS
      POINTEUR LINCO.MLMOTS
      real*8 nu
      CHARACTER*8 MTYP,NOME,NOMI,TYPE,NOMZ,TYPC
      CHARACTER*8 LTAB(NTB)
      DIMENSION KTAB(NTB),IXV(9)
      DATA LTAB/'KIZX    '/

C*****************************************************************************
C
C     write(6,*)' Debut CLMI '

      CALL LITABS(LTAB,KTAB,NTB,1,IRET)
      IF (IERR.NE.0) RETURN
      MTABX=KTAB(1)
C
C- Récupération de la table EQEX (pointeur MTAB1)
C
      CALL LEKTAB(MTABX,'EQEX',MTAB1)
      IF(MTAB1.EQ.0)THEN
C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
         MOTERR( 1: 8) = '  EQEX  '
         MOTERR( 9:16) = '  EQEX  '
         MOTERR(17:24) = '  KIZX  '
         CALL ERREUR(786)
         RETURN
      ENDIF
      CALL ACME(MTAB1,'NAVISTOK',NASTOK)
C
C- Récupération de la table INCO (pointeur KINC)
C
      CALL LEKTAB(MTAB1,'INCO',KINC)
      IF(KINC.EQ.0)THEN
C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
         MOTERR( 1: 8) = '  INCO  '
         MOTERR( 9:16) = '  INCO  '
         MOTERR(17:24) = '  EQEX  '
         CALL ERREUR(786)
         RETURN
      ENDIF

C*****************************************************************************
C OPTIONS
C KIMPL = 0 -> EXPL     1 -> IMPL     2 -> SEMI
C KFORM = 0 -> SI       1 -> EF       2 -> VF    3 -> EFMC
C IDCEN = 0-> rien  1-> CENTREE  2-> SUPGDC  3-> SUPG  4-> TVISQUEU  5-> CNG

      IAXI=0
      IF(IFOMOD.EQ.0)IAXI=2
C
C- Récupération de la table des options KOPT (pointeur KOPTI)
C
      CALL LEKTAB(MTABX,'KOPT',KOPTI)
      IF (KOPTI.EQ.0) THEN
C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
         MOTERR( 1: 8) = '  KOPT  '
         MOTERR( 9:16) = '  KOPT  '
         MOTERR(17:24) = '  KIZX  '
         CALL ERREUR(786)
         RETURN
      ENDIF

C     write(6,*)' Avant les options '
      CALL ACME(KOPTI,'IDCEN',IDCEN)
      KDIM=1
      IF(IDCEN.EQ.2)KDIM=IDIM
      CALL ACME(KOPTI,'IKOMP',IKOMP)
      CALL ACME(KOPTI,'KIMPL',KIMPL)
      CALL ACME(KOPTI,'KPOIN',KPRE)
      CALL ACME(KOPTI,'KFORM',KFORM)
      CALL ACMF(KOPTI,'CMD',CMD)

      IF(KFORM.NE.0.AND.KFORM.NE.1)THEN
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = 'EF/EFM1 '
            CALL ERREUR(803)
            RETURN
        ENDIF
      CALL ACMF(KOPTI,'AIMPL',AIMPL)
      AG=AIMPL
      AD=AIMPL-1.D0
      IF (IERR.NE.0) RETURN

C     write(6,*)' Apres les options '
C*****************************************************************************
C
C- Récupération de la table DOMAINE associée au domaine local
C
c      write(6,*)' Recuperation du modele '
      CALL ACMM(MTABX,'NOMZONE',NOMZ)
      TYPE=' '
      CALL ACMO(MTABX,'DOMZ',TYPE,MMDZ)
      IF(TYPE.NE.'MMODEL')THEN
C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
         MOTERR( 1: 8) = '  DOMZ  '
         MOTERR( 9:16) = '  DOMZ  '
         MOTERR(17:24) = '  KIZX  '
         CALL ERREUR(786)
         RETURN
      ENDIF

C E/  MMODEL : Pointeur de la table contenant l'information cherchée
C  /S IPOINT : Pointeur sur la table DOMAINE
C  /S INEFMD : Type formulation INEFMD=1 LINE,=2 MACRO,=3 QUADRATIQUE
C                               INEFMD=4 LINB

      CALL LEKMOD(MMDZ,MTABZ,INEFMD)
C verifier que la formulation est bonne  INEFMD = 1

      CALL LEKTAB(MTABZ,'MAILLAGE',MELEME)
      CALL LEKTAB(MTABZ,'SOMMET',MELEMS)
C     CALL LEKTAB(MTABZ,'CENTRE',MELEMC)
      CALL KRIPAD(MELEMS,MLENTI)
      IF (IERR.NE.0) RETURN


C*************************************************************************
C VERIFICATIONS SUR LES INCONNUES
C
C- Récupération du nombre d'inconnues et du nom de l'inconnue NOMI
C
c      write(6,*)' Recuperation inconnues '
      TYPE='LISTMOTS'
      CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO)
      IF (IERR.NE.0) RETURN
      SEGACT LINCO
      NBINC=LINCO.MOTS(/2)
      IF(NBINC.NE.1)THEN
C        Indice %m1:8 : contient plus de %i1 %m9:16
         MOTERR( 1:8) = 'LISTINCO'
         INTERR(1) = 1
         MOTERR(9:16) = ' MOTS   '
         CALL ERREUR(799)
         RETURN
      ENDIF

      NOMI=LINCO.MOTS(1)
c      write(6,*) 'Nb inco=',NBINC
c      write(6,*) 'NOMI=',NOMI
C
C- Récupération de l'inconnue
c      write(6,*) 'Récupération de l''inconnue'
      TYPE=' '
      CALL ACMO(KINC,NOMI,TYPE,MCHPO1)
      IF(TYPE.NE.'CHPOINT ')THEN
C        Indice %m1:8 : ne contient pas un objet de type %m9:16
         MOTERR( 1: 8) = 'INC '//NOMI
         MOTERR( 9:16) = 'CHPOINT '
         CALL ERREUR(800)
         RETURN
      ELSE
      CALL LICHT(MCHPO1,IZTU1,TYPC,MELEMI)
         NINKO = IZTU1.VPOCHA(/2)
         NPTI  = IZTU1.VPOCHA(/1)
c         write (6,*) 'NINKO =',NINKO 
c      write(6,*)' MCHPOI,MELEMI=',MCHPOI,MELEMI
C On fait pointer ces deux tableaux sur le champ U inconu (tjs présent) pour
C eviter de les enlever lors de l'appel FORTRAN si les options sont absentes
      ENDIF

C*****************************************************************************
C Le domaine de definition est donne par le SPG  de la premiere inconnue
C Les inconnues suivantes devront posseder ce meme pointeur
C On verifie que les points de la zone sont tous inclus dans ce SPG

      CALL KRIPAD(MELEMI,IPADI)

      CALL VERPAD(IPADI,MELEME,IRET)
      IF(IRET.NE.0)THEN
C     Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
         MOTERR(1: 8) = 'INC '//NOMI
         MOTERR(9:16) = 'CHPOINT '
         CALL ERREUR(788)
      RETURN
      ENDIF

C*****************************************************************************

      NU=1.d-6

C*************************************************************************
C Lecture des coefficients
C Type du coefficient :
C IK1=0 CHPOINT IK1=1 scalaire  IK1=2 vecteur

c     write(6,*)' lecture coef '
      CALL ACME(MTABX,'IARG',IARG)
      IF(IARG.NE.6)THEN
      WRITE(6,*)' Operateur DETO : '
      WRITE(6,*)' Nombre d''arguments incorrect : ',IARG
      WRITE(6,*)' On attend 6  '
C           Indice %m1:8 : nombre d'arguments incorrect
            MOTERR(1:8) = 'IARG    '
            CALL ERREUR(804)
      RETURN
      ENDIF


C--Cas incompréssible
c-------1er coefficient/type de fermeture     
      IXV(1)=0
      IXV(2)=1
      IXV(3)=0      
      CALL LEKCOF('Opérateur DETO :',
     & MTABX,KINC,1,IXV,IZTG1,MPOVA1,NPT1,NC1,IK1,IRET)
      IF(IRET.EQ.0)RETURN
      ferm = MPOVA1.VPOCHA(1,1)
c------ 2eme coef/n°de l'équation traitée
      IXV(1)=0
      IXV(2)=1
      IXV(3)=0     
      CALL LEKCOF('Opérateur DETO :',
     & MTABX,KINC,2,IXV,IZTG2,MPOVA2,NPT2,NC2,IK2,IRET)
      IF(IRET.EQ.0)RETURN
      ctype = MPOVA2.VPOCHA(1,1)
c------ Ue, vitesse extérieure
      IXV(1)=MELEMS
      IXV(2)=0
      IXV(3)=0     
      CALL LEKCOF('Opérateur DETO :',
     & MTABX,KINC,3,IXV,IZTG3,MPOVA3,NPT3,NC3,IK3,IRET)
      IF(IRET.EQ.0)RETURN
c------ dUe/dX, gradient de vitesse extérieure
      IXV(1)=MELEMS
      IXV(2)=0
      IXV(3)=0     
      CALL LEKCOF('Opérateur DETO :',
     & MTABX,KINC,4,IXV,IZTG4,MPOVA4,NPT4,NC4,IK4,IRET)
      IF(IRET.EQ.0)RETURN
c------ lecture de l'inconnue D2 au temps n-1
      IXV(1)=MELEMS
      IXV(2)=0
      IXV(3)=0     
      CALL LEKCOF('Opérateur DETO :',
     & MTABX,KINC,5,IXV,IZTG5,MPOVA5,NPT5,NC5,IK5,IRET)
      IF(IRET.EQ.0)RETURN 
c------ lecture de l'inconnue D3 au temps n-1(non obligatoire)
      IXV(1)=MELEMS
      IXV(2)=0
      IXV(3)=0     
      CALL LEKCOF('Opérateur DETO :',
     & MTABX,KINC,6,IXV,IZTG6,MPOVA6,NPT6,NC6,IK6,IRET)
      IF(IRET.EQ.0)RETURN

      
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C On cree  le champoint second membre
       NAT=2
       NSOUPO=1
       SEGINI MCHPOI
       JATTRI(1)=2
       NC=1
       SEGINI MSOUPO
       IPCHP(1)=MSOUPO
       NOCOMP(1)=NOMI(1:4)
       NOCONS(1)='QDM '
       IGEOC=MELEMS
       N=NPTI
       SEGINI MPOVAL
       IPOVAL=MPOVAL
       
c       write (6,*) 'longueur1 chpoint = ',VPOCHA(/1)
c       write (6,*) 'longueur2 chpoint = ',VPOCHA(/2)
c       write (6,*) 'NPTI = ',NPTI

C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C On crée le champoint qui reçoit les valeurs du facteur de forme     
       SEGINI MCHPO9
       MCHPO9.JATTRI(1)=2
       NC=1
       SEGINI MSOUP9
       MCHPO9.IPCHP(1)=MSOUP9
       MSOUP9.NOCOMP(1)='SCAL'
       MSOUP9.NOCONS(1)='H   '
       MSOUP9.IGEOC=MELEMS
       SEGINI MPOVA9                                                        
       MSOUP9.IPOVAL=MPOVA9
C On crée le champoint qui reçoit les valeurs du coefficient de frottement     
       SEGINI MCHPO8
       MCHPO8.JATTRI(1)=2
       NC=1
       SEGINI MSOUP8
       MCHPO8.IPCHP(1)=MSOUP8
       MSOUP8.NOCOMP(1)='SCAL'
       MSOUP8.NOCONS(1)='Cf  ' 
       MSOUP8.IGEOC=MELEMS
       SEGINI MPOVA8                                                        
       MSOUP8.IPOVAL=MPOVA8
C On crée le champoint qui reçoit les valeurs del'épaisseur de déplacement
       SEGINI MCHPOA
       MCHPOA.JATTRI(1)=2
       NC=1
       SEGINI MSOUPA
       MCHPOA.IPCHP(1)=MSOUPA
       MSOUPA.NOCOMP(1)='SCAL'
       MSOUPA.NOCONS(1)='D1  ' 
       MSOUPA.IGEOC=MELEMS
       SEGINI MPOVAA                                                        
       MSOUPA.IPOVAL=MPOVAA
C On crée le champoint qui reçoit les valeur de d(UE*delta1)/dx
       SEGINI MCHPOB
       MCHPOB.JATTRI(1)=2
       NC=1
       SEGINI MSOUPB
       MCHPOB.IPCHP(1)=MSOUPB
       MSOUPB.NOCOMP(1)='SCAL'
       MSOUPB.NOCONS(1)='Vtra' 
       MSOUPB.IGEOC=MELEMS
       SEGINI MPOVAB                                                        
       MSOUPB.IPOVAL=MPOVAB
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

C     MTYP='MAILLAGE'
C     CALL LIROBJ(MTYP,MELEME,1,IRETOU)
C     IF (IRETOU.NE.1)THEN
C     MOTERR( 1: 8) ='MAILLAGE'
C     CALL ERREUR(471)
C     RETURN
C     ENDIF

c      write(6,*)' On peut commencer a travailler '
c      write(6,*)' La dimension de l espace IDIM=',idim

      SEGACT MELEME
      NBSOUS=LISOUS(/1)
c      write(6,*)' Nb de sous maillage NBSOUS=',NBSOUS
      IF (NBSOUS.NE.0)RETURN
      IF (NBSOUS.EQ.0)NBSOUS=1

c      write(6,*)' ITYPEL=',ITYPEL,NOMS(ITYPEL)
      NOME=NOMS(ITYPEL)//'    '
      CALL KALPBG(NOME,'FONFORM ',IZFFM)

      SEGACT IZFFM*MOD
      IZHR=KZHR(1)
      SEGACT IZHR*MOD

      NES=GR(/1)
      NPG=GR(/3)
      IAXI=0

      NBELEM=NUM(/2)
      NBNN  =NUM(/1)
      write(6,*)' NBELEM=',nbelem,' NBNN=',nbnn
      NBEL = NBELEM
      NP=NBNN
      MP=NBNN

c      SEGINI MTEL
      SEGINI TRAV
      CALL INITD(UN,(NPTI*IDIM),1.D0)
      CALL INITD(TN,NPTI,1.D0)

C creation de l'objet MATRIK

      NRIGE=7
      NKID =9
      NKMT =7
      NMATRI=2

      SEGINI MATRIK
      CALL ECROBJ('MATRIK',MATRIK)
      CALL ECROBJ('CHPOINT',MCHPOI)

      NBME=1

c matrice masse élémentaire
      SEGINI IMAT1
      IRIGEL(1,1)=MELEME
      IRIGEL(2,1)=MELEME
      IRIGEL(4,1)=IMAT1
      IRIGEL(7,1)=0

      IMAT1.LISPRI(1)=NOMI
      IMAT1.LISDUA(1)=NOMI
      IMAT1.KSPGP=MELEMS
      IMAT1.KSPGD=MELEMS
      SEGINI IPM1
      IMAT1.LIZAFM(1,1)=IPM1
c matrice de convection
      SEGINI IMAT2
      IRIGEL(1,2)=MELEME
      IRIGEL(2,2)=MELEME
      IRIGEL(4,2)=IMAT2
      IRIGEL(7,2)=2

      IMAT2.LISPRI(1)=NOMI
      IMAT2.LISDUA(1)=NOMI
      IMAT2.KSPGP=MELEMS
      IMAT2.KSPGD=MELEMS
      SEGINI IPM2
      IMAT2.LIZAFM(1,1)=IPM2
      nd1=ipm2.am(/1)
      nd2=ipm2.am(/2)
      nd3=ipm2.am(/3)

      
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C>>>>>>>>>>Relations de fermeture >>>>>>>>>>>>>>>>>>>>>>>>>>
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      IF (ferm.EQ.1.) THEN
         write(6,*) 'CAS LAMINAIRE: Méthode simplifiée'
         IF(ctype.EQ.1.) THEN       
C----Cas laminaire, méthode simplifiée
            write(6,*) 'EQUATION DE QDM'
            DO  I=1,VPOCHA(/1)
               HNM(I) = 2.5911
               BNM(I) = 0.22052
               MPOVA9.VPOCHA(I,1) = HNM(I)
            ENDDO
         ENDIF
      ENDIF
   
      IF(ferm.EQ.2.) THEN         
C----Cas laminaire, méthode de Karman Pohlhausen
         write(6,*) 'CAS LAMINAIRE: Méthode de Von Karman-Pohlhausen'
         IF(ctype.EQ.1.) THEN
            write(6,*) 'EQUATION DE QDM'
            CALL Karpohl (MPOVA5.VPOCHA,MPOVA4.VPOCHA,MPOVA3.VPOCHA,
     &                       VPOCHA(/1),HNM,BNM)
            DO I=1,VPOCHA(/1)
               Coef1(I)= (HNM(I)+2)*MPOVA4.VPOCHA(I,1)
     &                     *((MPOVA3.VPOCHA(I,1))**(-1))
               Coef2(I)=(BNM(I)*1e-6)*((MPOVA3.VPOCHA(I,1)
     &                     *MPOVA5.VPOCHA(I,1))**(-1))
               MPOVA9.VPOCHA(I,1) = HNM(I)

            ENDDO
         ENDIF
      ENDIF
   
      IF(ferm.EQ.3.) THEN
C----Cas laminaire, méthode à deux équations
         write(6,*) 'CAS LAMINAIRE: Méthode à 2 équations'
         IF(ctype.EQ.1.) THEN
C--------Equation de Qdm 
            write(6,*) 'EQUATION DE QDM'
            CALL LAM2(MPOVA5.VPOCHA,MPOVA6.VPOCHA,H32NM,HNM,BNM,DNM,D1N,
     &                VPOCHA(/1))
            DO I=1,VPOCHA(/1)
               Coef1(I)= (HNM(I)+2)*MPOVA4.VPOCHA(I,1)
     &                     *((MPOVA3.VPOCHA(I,1))**(-1))
               Coef2(I)=((BNM(I)*1e-6)*((MPOVA3.VPOCHA(I,1)
     &                     *MPOVA5.VPOCHA(I,1))**(-1)))
               MPOVA9.VPOCHA(I,1)=HNM(I)
               MPOVA8.VPOCHA(I,1)=2*Coef2(I)
               MPOVAA.VPOCHA(I,1)=D1N(I)
            ENDDO
         ENDIF         
         IF(ctype.EQ.2.) THEN
c--------Equation d'énegie cinétique
            write(6,*) 'EQUATION D''ENERGIE CINETIQUE'
            CALL LAM2(MPOVA6.VPOCHA,MPOVA5.VPOCHA,H32NM,HNM,BNM,DNM,D1N,
     &                VPOCHA(/1))
            DO I=1,VPOCHA(/1)
               Coef1(I)=3*MPOVA4.VPOCHA(I,1)*
     &              ((MPOVA3.VPOCHA(I,1))**(-1))
               Coef2(I)=(2*DNM(I)*1e-6*((MPOVA3.VPOCHA(I,1)
     &              *MPOVA5.VPOCHA(I,1))**(-1))*H32NM(I))
               MPOVA9.VPOCHA(I,1)=HNM(I)
               MPOVAA.VPOCHA(I,1)=D1N(I)
            ENDDO
         ENDIF
      ENDIF

      IF(ferm.EQ.4.) THEN
c----Cas Turbulent, méthode à deux équations (Cousteix)
         write(6,*) 'CAS TURBULENT: Méthode à 2 équations(Cousteix)'
         IF(ctype.EQ.1.) THEN
c--------Equation de qdm
            write(6,*) 'EQUATION DE QDM'
            CALL TURB5(MPOVA5.VPOCHA,MPOVA6.VPOCHA,MPOVA3.VPOCHA,
     &           HHN,HNM,CFN,CEN,D1N,DeN,VPOCHA(/1))
            DO I=1,VPOCHA(/1)
               Coef1(I)= (HNM(I)+2)*MPOVA4.VPOCHA(I,1)
     &                     *((MPOVA3.VPOCHA(I,1))**(-1))
               Coef2(I)=CFN(I)/2
               MPOVA9.VPOCHA(I,1)=HNM(I)
               MPOVA8.VPOCHA(I,1)=CFN(I)
               MPOVAA.VPOCHA(I,1)=D1N(I)               
            ENDDO
         ENDIF
         IF(ctype.EQ.3.) THEN
c--------Equation d'entrainement
            write(6,*) 'EQUATION D''ENTRAINEMENT'
            CALL TURB5(MPOVA6.VPOCHA,MPOVA5.VPOCHA,MPOVA3.VPOCHA,
     &           HHN,HNM,CFN,CEN,D1N,DeN,VPOCHA(/1))
            DO I=1,VPOCHA(/1)
               Coef1(I)=(MPOVA4.VPOCHA(I,1))
     &                  *((MPOVA3.VPOCHA(I,1))**(-1))
               Coef2(I)=CEN(I)
               MPOVA9.VPOCHA(I,1)=HNM(I)
               MPOVA8.VPOCHA(I,1)=CFN(I)
               MPOVAA.VPOCHA(I,1)=D1N(I)               
            ENDDO   
         ENDIF
      ENDIF
      IF(ferm.EQ.5.) THEN
c----Cas Turbulent, méthode à deux équations relations de Head
         write(6,*) 'CAS TURBULENT: Méthode à 2 équations(Head)'
         IF(ctype.EQ.1.) THEN
c--------Equation de qdm
            write(6,*) 'EQUATION DE QDM'
c            write(6,*)MPOVA3.VPOCHA(/1),MPOVAL.VPOCHA(/1)
            CALL TURB6(MPOVA5.VPOCHA,MPOVA6.VPOCHA,MPOVA3.VPOCHA,
     &           HHN,HNM,CFN,CEN,D1N,VPOCHA(/1))
            DO I=1,VPOCHA(/1)
               Coef1(I)= (HNM(I)+2)*MPOVA4.VPOCHA(I,1)
     &                     *((MPOVA3.VPOCHA(I,1))**(-1))
               Coef2(I)=CFN(I)/2
               MPOVA9.VPOCHA(I,1)=HNM(I)
               MPOVA8.VPOCHA(I,1)=CFN(I)
               MPOVAA.VPOCHA(I,1)=D1N(I)
            ENDDO
         ENDIF
         IF(ctype.EQ.3.) THEN
c--------Equation d'entrainement
            write(6,*) 'EQUATION D''ENTRAINEMENT'
            CALL TURB6(MPOVA6.VPOCHA,MPOVA5.VPOCHA,MPOVA3.VPOCHA,
     &           HHN,HNM,CFN,CEN,D1N,VPOCHA(/1))
            DO I=1,VPOCHA(/1)
               Coef1(I)=(MPOVA4.VPOCHA(I,1))
     &                  *((MPOVA3.VPOCHA(I,1))**(-1))
               Coef2(I)=CEN(I)
               MPOVA9.VPOCHA(I,1)=HNM(I)
               MPOVA8.VPOCHA(I,1)=CFN(I)
               MPOVAA.VPOCHA(I,1)=D1N(I)
            ENDDO   
         ENDIF
      ENDIF
c      IF(ferm.EQ.6) THEN
C---- Convection naturelle laminaire
c         write(6,*) 'CONVECTION NATURELLE'
c         IF(ctype.EQ.1) THEN
C--------Equation de quantité de mouvement
c            write(6,*) 'EQUATION DE QDM'
c            CALL CNAT(MPOVA5.VPOCHA,MPOVA6.VPOCHA,MPOVA3.VPOCHA,
c     &              MPOVA4.VPOCHA,UCL,Delta,SM1,SM2,FP,PHIP,VPOCHA(/1))        
c            DO I=1,VPOCHA(/1)
c               Coef1(I) = 0
c               Coef2(I) = SM1(I)
c               MPOVA9.VPOCHA(I,1)=FP(I)
c               MPOVA8.VPOCHA(I,1)=PHIP(I)
c               MPOVAA.VPOCHA(I,1)=Delta(I)
c            ENDDO   
c         ENDIF
c         IF(ctype.EQ.4.) THEN
C--------Equation de la chaleur
c            write(6,*) 'EQUATION DE LA CHALEUR'
c            CALL CNAT(MPOVA5.VPOCHA,MPOVA6.VPOCHA,MPOVA3.VPOCHA,
c     &              MPOVA4.VPOCHA,UCL,Delta,SM1,SM2,FP,PHIP,VPOCHA(/1)) 
c            DO I=1,VPOCHA(/1)
c               Coef1(I) = 0
c               Coef2(I) = SM2(I)
c               MPOVA9.VPOCHA(I,1)=FP(I)
c               MPOVA8.VPOCHA(I,1)=PHIP(I)
c               MPOVAA.VPOCHA(I,1)=Delta(I)
c            ENDDO
c        ENDIF
c      ENDIF

C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c      write(6,*)' nd1,nd2,nd3=',nd1,nd2,nd3
c      write(6,*)' NBELEM=',NBELEM


      DO 101 K=1,NBELEM
c       write(6,*)' Element n ',k
c      write(6,1001)(num(ii,k),ii=1,nbnn)

C?    DO 102 I=1,NBNN
C?    IP1=num(I,K)
C?    XYZ(1,I)=XCOOR((IP1-1)*(IDIM+1)    +1)
C?    XYZ(2,I)=XCOOR((IP1-1)*(IDIM+1)    +2)
c      write(6,*)'Coordonnees element ',K,' noeud ',ip1
c      write(6,1002)(XYZ(M,I),M=1,IDIM)
C? 102  CONTINUE

C?    CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,
C?   &   NES,IDIM,NBNN,NPG,IAXI,AIRE,AJ)

      KDEB=1
      KFIN=1
      LRV =1
      IKOMP=0
      ANUK(1)=1.D0

      CALL SUPGEF(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ,
     & NES,NP,NPG,IAXI,XCOOR,
     & NUM(1,K),KDEB,KFIN,LRV,IDCEN,CMD,IKOMP,
     & TN,LECT,UN,LECT,NPTI,ANUK,
     & WT,WS,HR,PGSQ,RPG,AJ,AIRE,UIL,DUIL,
     & DTM1,DT,DTT1,DTT2,DIAEL,NUEL)

C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C>>>>>>>>>>Calcul des matrices élémentaires>>>>>>>>>>>>>>>>>
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c calcul des matrices masses élémentaires
c      write(6,*)' NBNN=',nbnn
      DO 103 I=1,NBNN
         DO 104 J=1,NBNN
            U=0.D0
            DO 105 L=1,NPG
              V=0.D0
              DO 1051 L1=1,NBNN
                 J1=lect(num(L1,K))
                 V=V+FN(L1,L)*Coef1(J1)
 1051         CONTINUE
              U=U+(WT(J,L)*FN(I,L)*PGSQ(L))*V
 105        CONTINUE
            IPM1.AM(K,I,J)=U
c            write(6,*) 'M(',I,',',J,',',K,')=',IPM1.AM(K,I,J)
 104     CONTINUE
 103  CONTINUE

      ax=1.d0
      ay=1.d0
c calcul des matrices de convection élémentaires      
      DO 106 I=1,NBNN
         DO 107 J=1,NBNN
c            J1 = lect(num(J,K))
            U=0.D0
            DO 108 L=1,NPG
               V=0.D0
               DO 1081 L1=1,NBNN
                  J1 = lect(num(L1,K))
                  V = V+FN(L1,L)*1.
 1081          CONTINUE
               U=U+(WT(J,L)*((ax*HR(1,I,L)) )*PGSQ(L))*V
c              write(6,*) 'hr(',J,',',L,')=',HR(1,J,L)
c              write(6,*) 'FN(',I,',',L,',',K,')=',FN(I,L)
c              write(6,*) 'WT(',I,',',L,',',K,')=',WT(I,L)
c              U=U+(WT(I,L)*((ax*HR(1,J,L)) + (ay*HR(1,J,L)) )*PGSQ(L))
c               U=U+(WT(J,L)*((ax*HR(1,I,L)) )*PGSQ(L))
c     &*MPOVA2.VPOCHA(J1,1)
 108        CONTINUE
C           write(6,*)'K,i,j=', K,i,j
            IPM2.AM(K,I,J) = U
 107     CONTINUE
 106  CONTINUE

C calcul du terme source
      DO 109 I=1,NBNN
         I1= lect(num(I,K))
         U=0.D0
         DO 110 L=1,NPG
            V=0.D0
            v2=0.D0
            DO 1101 J=1,NBNN
               J1 =lect(num(J,K))
               V=V+FN(J,L)*Coef2(J1)
 1101       CONTINUE
            U=U+V*WT(I,L)*PGSQ(L)          
 110    CONTINUE
        VPOCHA(I1,1)=VPOCHA(I1,1)+ U
 109  CONTINUE

C calcul de la vitesse de transpiration
      DO 111 I=1,NBNN
         I1 = lect(num(I,K))
         U=0.D0
         DO 112 J=1,NBNN
            J1 = lect(num(J,K))
            U = U+HR(1,J,1)*(D1N(J1)*MPOVA3.VPOCHA(J1,1))
 112     CONTINUE
         VT(I1) = U
c         write(6,*) 'VT(',I1,')=',VT(I1)
         MPOVAB.VPOCHA(I1,1)= VT(I1)
 111  CONTINUE 

101   CONTINUE


C Sortie des valeurs des coefficients H, Cf, Ce et VT

      CALL ECMO(KINC,'H','CHPOINT',MCHPO9)      
      CALL ECMO(KINC,'CF','CHPOINT',MCHPO8)
      CALL ECMO(KINC,'D1','CHPOINT',MCHPOA)
      CALL ECMO(KINC,'Vtra','CHPOINT',MCHPOB)
     

      SEGDES IPM2,IMAT2,IPM1,IMAT1
      SEGDES MATRIK
      SEGDES TRAV
      SEGDES IZFFM,IZHR
      SEGDES MELEME
      SEGDES MPOVAL,MPOVA8,MPOVA9,MPOVAA,MPOVAB
      SEGDES MSOUPO,MSOUP8,MSOUP9,MSOUPA,MSOUPB
      SEGDES MCHPOI,MCHPO9,MCHPO8,MCHPOA,MCHPOB
      SEGDES LINCO

C     SEGDES MTEL

      RETURN
 1003 FORMAT(6(1X,1PE11.4))
 1002 FORMAT(10(1X,1PE11.4))
 1001 FORMAT(20(1X,I5))
      END