Sources Esope | Cast3M

clmi
C CLMI      SOURCE    CB215821  23/01/25    21:15:07     11573          
      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 PPARAM
-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
       IFOPOI = IFOUR
       MTYPOI = '        '
       MOCHDE = 'cree par CLIM'
       JATTRI(1)=2
       NC=1
       SEGINI MSOUPO
       IPCHP(1)=MSOUPO
       NOCOMP(1)=NOMI
       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.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.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.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.IGEOC=MELEMS
       SEGINI MPOVAB
       MSOUPB.IPOVAL=MPOVAB
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
C     MTYP='MAILLAGE'
C     CALL LIROBJ(MTYP,MELEME,1,IRETOU)
C     CALL ACTOBJ(MTYP,MELEME,1)
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
 
      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
 
C     CB215821 : DT n'est pas utilise mais doit etre initialise sinon NAN
      DT=0.D0
 
      SEGACT,MCOORD
      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)
      SEGDES,MCOORD
 
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 LINCO
C     SEGDES MTEL
 
      CALL ECROBJ('MATRIK  ',MATRIK)
      CALL ACTOBJ('CHPOINT ',MCHPOI,1)
      CALL ECROBJ('CHPOINT ',MCHPOI)
      RETURN
 1003 FORMAT(6(1X,1PE11.4))
 1002 FORMAT(10(1X,1PE11.4))
 1001 FORMAT(20(1X,I5))
      END