Télécharger yfpu.eso

Retour à la liste

yfpu
C YFPU      SOURCE    GOUNAND   25/11/12    21:15:48     12399          
      SUBROUTINE YFPU
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C***********************************************************************
C
C     SYNTAXE :
C
C        FPU  NU UET YP <,VPAROI>
C
C       1------2
C       (R1,AL1)  LEF   FLUIDE  NOEUDS 1 2
C
C
C     ANU     VISCOSITE CINEMATIQUE
C     UET     U*
C     YP      DISTANCE A LA PAROI
C     VPAROI  VITESSE DE LA PAROI (PAR DEFAUT 0.)
C
C   CAS TRIDIMENSIONNEL
C                       4 ________ 3
C                       / FLUIDE /
C                    1 /________/2
C
C
C***********************************************************************
 
-INC CCVQUA4
 
-INC PPARAM
-INC CCOPTIO
-INC CCGEOME
-INC SMCOORD
-INC SIZFFB
      POINTEUR IZF1.IZFFM,IZH2.IZHR
-INC SMLREEL
      POINTEUR MLREE4.MLREEL
-INC SMLENTI
      POINTEUR IPADS.MLENTI,IPADT.MLENTI
-INC SMMATRIK
-INC SMELEME
      POINTEUR MELEMS.MELEME,MELEKE.MELEME,MELEMX.MELEME,MLEMST.MELEME
-INC SMCHPOI
      POINTEUR IZD.MPOVAL,IZDS.MPOVAL,IZDK.MPOVAL,IZDE.MPOVAL
      POINTEUR IZG1.MCHPOI, IZGG1.MPOVAL
      POINTEUR IZG2.MCHPOI, IZGG2.MPOVAL
      POINTEUR IZG3.MCHPOI, IZGG3.MPOVAL
      POINTEUR IZGV.MCHPOI, IZGGV.MPOVAL
      POINTEUR IZTU1.MPOVAL
      POINTEUR MUE.MCHPOI,MZUE.MPOVAL
      POINTEUR MYP.MCHPOI,MZYP.MPOVAL
      POINTEUR MZNU.MPOVAL
      POINTEUR MZRO.MPOVAL,MZMU.MPOVAL
      POINTEUR IZVOL.MPOVAL
      POINTEUR MNT.MCHPOI,MZNT.MPOVAL
 
-INC SMLMOTS
      POINTEUR LINCO.MLMOTS
      CHARACTER*8 NOMZ,NOMI,TYPE,TYPC,NOM0
      CHARACTER*(LOCOMP) NOM4(5)
      LOGICAL LOGI
      PARAMETER (NTB=1)
      CHARACTER*8 LTAB(NTB)
      DIMENSION KTAB(NTB),IXV(3),AJ(270),XYZ(90)
      DATA LTAB/'KIZX    '/
C*****************************************************************************
CFPU
      segact mcoord
*pv  ca le fait mieux en initialisant ikel
*pv  si on ne le fait pas et que par hasard il vaut 1
*pv  les resultats ne sont pas ceux attendus
      ikel=0
c      write(6,*)' debut FPU '
      izdk=0
      izde=0
      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)
      IF(NASTOK.EQ.0)THEN
      CALL ZFPU(MTABX,MTAB1)
      RETURN
      ENDIF
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
 
      IKOMP=0
      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 '
      TYPE=' '
      CALL ACMO(KOPTI,'PRECAU',TYPE,IENT)
      IF(TYPE.EQ.'ENTIER')CALL ACME(KOPTI,'PRECAU',IPRE)
      TYPE=' '
      CALL ACMO(KOPTI,'IKOMP',TYPE,IENT)
      IF(TYPE.EQ.'ENTIER')CALL ACME(KOPTI,'IKOMP',IKOMP)
 
      CALL ACME(KOPTI,'KIMPL',KIMPL)
      CALL ACME(KOPTI,'KFORM',KFORM)
 
      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
 
C*****************************************************************************
C
C- Récupération de la table DOMAINE associée au domaine local
C
      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)
 
      CALL LEKTAB(MTABZ,'SOMMET',MELEMS)
      CALL LEKTAB(MTABZ,'MAILLAGE',MELEME)
      CALL LEKTAB(MTABZ,'CENTRE',MELEMC)
      CALL LEKTAB(MTABZ,'XXVOLUM',MCHPOI)
      IF (IERR.NE.0) RETURN
 
      SEGACT MELEME
      CALL LICHT(MCHPOI,IZVOL,TYPC,IGEOM)
      IF (IERR.NE.0) RETURN
C***
 
      TYPE='LISTMOTS'
      CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO)
      SEGACT LINCO
 
      CALL ACME(MTABX,'IARG',IARG)
      IKOMP=0
 
      IF (KFORM.EQ.0.AND.IARG.NE.3)THEN
      write(6,*)' FPU : nb d arguments incorrect :',IARG
      GO TO 90
      ENDIF
      IF (KFORM.EQ.1.AND.IARG.NE.5.AND.IARG.NE.6)THEN
      write(6,*)' FPU : nb d arguments incorrect :',IARG
      GO TO 90
      ENDIF
      IKRS=0
      IKR=1
      N=1
      NC=1
      SEGINI MZRO
      MZRO.VPOCHA(1,1)=1.D0
 
 
C--Cas incompressible
      IF(IARG.EQ.3)THEN
      IF(IKOMP.EQ.0)THEN
C 1er coefficient : nu
      IXV(1)=0
      IXV(2)=1
      IXV(3)=0
      IRET=0
      CALL LEKCOF('Opérateur FPU :',
     & MTABX,KINC,1,IXV,MNU,MZNU,NPT1,NC1,IKM,IRET)
      IF(IRET.EQ.0)RETURN
 
C 2ème coefficient : uet
      IXV(1)=MELEMC
      IF(KFORM.EQ.1)IXV(1)=MELEMS
      IXV(2)=0
      IXV(3)=0
      IRET=1
      CALL LEKCOF('Opérateur FPU :',
     & MTABX,KINC,2,IXV,MUE,MZUE,NPT2,NC2,IK2,IRET)
      IF(IRET.EQ.0)RETURN
      SEGACT MZUE*MOD
 
C 3ème coefficient : yp
      IXV(1)=0
      IXV(2)=1
      IXV(3)=0
      IRET=0
      CALL LEKCOF('Opérateur FPU :',
     & MTABX,KINC,3,IXV,MYP,MZYP,NPT3,NC3,IK3,IRET)
      IF(IRET.EQ.0)RETURN
      ENDIF
 
 
C--Cas compressible
      ELSEIF(IKOMP.EQ.1)THEN
C 1er coefficient : mu
      IXV(1)=MELEMC
      IXV(2)=1
      IXV(3)=0
      IRET=0
      CALL LEKCOF('Opérateur FPU :',
     & MTABX,KINC,1,IXV,MNU,MZNU,NPT1,NC1,IKM,IRET)
      IF(IRET.EQ.0)RETURN
 
C 2ème coefficient : uet
      IXV(1)=MELEMC
      IXV(2)=0
      IXV(3)=0
      IRET=0
      CALL LEKCOF('Opérateur FPU :',
     & MTABX,KINC,2,IXV,MUE,MZUE,NPT2,NC2,IK2,IRET)
      IF(IRET.EQ.0)RETURN
      SEGACT MZUE*MOD
 
      ENDIF
 
      CALL KRIPAD(MELEMS,IPADS)
      NPTS=MELEMS.NUM(/2)
 
C*****************************************************************************
 
C VERIFICATIONS SUR LES INCONNUES
      NBINC=LINCO.MOTS(/2)
      IF(NBINC.NE.1.AND.NBINC.NE.3)THEN
      WRITE(6,*)'Nombre d''inconnues incorrect : ',NBINC,' On attend 3'
      WRITE(6,*)' On attend 1 ou 3'
C           Indice %m1:8 : nombre d'arguments incorrect
            MOTERR(1:8) = 'INCO    '
            CALL ERREUR(804)
      RETURN
      ENDIF
 
C --> 1 ere Inconnue
 
      NOMI=LINCO.MOTS(1)
      DO 15 I=1,IDIM
      WRITE(NOM4(I),FMT='(I1,A3)')I,NOMI(1:3)
 15   CONTINUE
 
      TYPE=' '
      CALL ACMO(KINC,NOMI,TYPE,MCHPOI)
      IF(TYPE.NE.'CHPOINT ')THEN
      WRITE(6,*)' L objet CHPOINT ',NOMI,' n existe pas dans la table'
      RETURN
      ELSE
      CALL LICHT(MCHPOI,IZTU1,TYPC,IGEOM0)
      ENDIF
      CALL KRIPAD(IGEOM0,MLENTI)
 
 
C._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
C--Cas Complet: ro un mu uet yp
      IF((IARG.EQ.5.OR.IARG.EQ.6).AND.KFORM.EQ.1)THEN
 
C 1er coefficient : ro
      IKRS=1
      IXV(1)=IGEOM0
      IXV(2)=1
      IXV(3)=0
      IRET=0
      CALL LEKCOF('Opérateur FPU :',
     & MTABX,KINC,1,IXV,MRO,MZRO,NPT1,NC1,IKR,IRET)
      IF(IRET.EQ.0)RETURN
 
C 2ème coefficient : un
      IXV(1)=-IGEOM0
      IXV(2)=0
      IXV(3)=0
      IRET=1
      CALL LEKCOF('Opérateur FPU :',
     & MTABX,KINC,2,IXV,MUE,MZUE,NPT2,NC2,IKU,IRET)
      IF(IRET.EQ.0)RETURN
      SEGACT MZUE*MOD
 
C 3ème coefficient : mu
      IXV(1)=IGEOM0
      IXV(2)=1
      IXV(3)=0
      IRET=0
      CALL LEKCOF('Opérateur FPU :',
     & MTABX,KINC,3,IXV,MMU,MZMU,NPT3,NC3,IKM,IRET)
      IF(IRET.EQ.0)RETURN
 
C 4ème coefficient : uet
      IXV(1)=MELEMC
      IF(KFORM.EQ.1)IXV(1)=MELEMS
      IXV(2)=0
      IXV(3)=0
      IRET=1
 
      CALL LEKCOF('Opérateur FPU :',
     & MTABX,KINC,4,IXV,MUE,MZUE,NPT2,NC2,IK2,IRET)
      IF(IRET.EQ.0)RETURN
      SEGACT MZUE*MOD
 
      TYPE=' '
      CALL ACMO(MTABX,'ARG5',TYPE,MDTOT)
 
c#######################################################################
c################# IKEL = 1 => Syntaxe II ##############################
      IF(TYPE.EQ.'MMODEL')THEN
      IKEL=1
 
      TYPE=' '
      CALL ACMO(MTABX,'PREFPU',TYPE,MTABP)
      IF(TYPE.NE.'TABLE')THEN
      KPREFPU=0
      write(6,*)'Operateur FPU : On Preconditionne '
      CALL CRTABL(MTABP)
      CALL ECMO(MTABX,'PREFPU','TABLE',MTABP)
      ELSE
      KPREFPU=1
      ENDIF
 
C 6ème coefficient : NUT
      IXV(1)=IGEOM0
      IXV(2)=0
      IXV(3)=0
      IRET=0
      CALL LEKCOF('Opérateur FPU :',
     & MTABX,KINC,6,IXV,MNT,MZNT,NPT3,NC3,IK3,IRET)
      IF(IRET.EQ.0)RETURN
C On calcule le Gradient de U
 
*  IGEOM0 Support géomtrique domaine complet
*  MDTOT modèle pour le maillage complet
*  MTABT table domaine en résultant
*  MTABZ table domaine de la zone frontière
 
       CALL LEKMOD(MDTOT,MTABT,INEFMDT)
       CALL LEKTAB(MTABT,'MAILLAGE',MELEMT)
       CALL LEKTAB(MTABT,'SOMMET'  ,MLEMST)
       CALL LEKTAB(MTABZ,'SOMMET'  ,MELEMS)
 
C._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
c on extrait les élements s'appuyant largement sur la frontière
c -> MELEMX
      IF(KPREFPU.EQ.0)THEN
 
      SEGACT MUE,MZUE
      NSOUPO=MUE.IPCHP(/1)
      IF(NSOUPO.GT.1)THEN
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = 'NSOUPO>1'
            CALL ERREUR(803)
            RETURN
      ENDIF
      NAT=MUE.JATTRI(/1)
      N =MZUE.VPOCHA(/1)
      NC=MZUE.VPOCHA(/2)
      SEGINI MYP,MSOUPO,MZYP
      DO 651 I=1,NAT
      MYP.JATTRI(I)=MUE.JATTRI(I)
 651  CONTINUE
      SEGINI MSOUPO
      MYP.IPCHP(1)=MSOUPO
      MSOUP1=MUE.IPCHP(1)
      SEGACT MSOUP1
      DO 652 I=1,NC
      NOCOMP(I)=MSOUP1.NOCOMP(I)
      NOHARM(I)=MSOUP1.NOHARM(I)
 652  CONTINUE
      IGEOC =MSOUP1.IGEOC
      IPOVAL=MZYP
 
       CALL ECROBJ('MAILLAGE',MELEME)
       CALL ECRCHA('POI1')
       CALL PRCHAN
       CALL ECRCHA('LARGEMENT')
       CALL ECRCHA('APPUYE')
       CALL ECROBJ('MAILLAGE',MELEMT)
       CALL PREXTR
       CALL LIROBJ('MAILLAGE',MELEMX,1,IRET)
       CALL ECROBJ('MAILLAGE',MELEMX)
       CALL NBNO
       CALL LIRENT(NPTA,1,IRET)
 
       call ecmo(mtabx,'MYP'   ,'CHPOINT' ,MYP   )
       call ecmo(mtabp,'MELEMX','MAILLAGE',MELEMX)
       call ecme(mtabp,'NPTA',npta)
      ELSE
       type=' '
       call acmo(mtabx,'MYP',TYPE,MYP)
       SEGACT MYP
       MSOUPO=MYP.IPCHP(1)
       SEGACT MSOUPO
       MZYP=MSOUPO.IPOVAL
       SEGACT MZYP
       type=' '
       call acmo(mtabp,'MELEMX',TYPE,MELEMX)
       call acme(mtabp,'NPTA',npta)
      ENDIF
c -> MELEMX
c on extrait les elements s'appuyant largement sur la frontiere
C._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
 
c     write(6,*)' FPU : NBSOUS=',NBSOUS
c     CALL ECROBJ('MAILLAGE',MELEMT)
c     CALL PRTRAC
c     CALL LIROBJ('MAILLAGE',MELEMT,1,IRET)
c     CALL ECROBJ('MAILLAGE',MELEMX)
c     CALL ECROBJ('MAILLAGE',MELEMS)
c     CALL PRFUSE
c     CALL PRTRAC
c     CALL LIROBJ('MAILLAGE',MELEMX,1,IRET)
 
       CALL KRIPAD(MLEMST,IPADT)
       SEGACT MLEMST
       NPTT=MLEMST.NUM(/2)
       CALL KRIPAD(MELEMX,MLENT2)
       SEGACT MELEMS
       CALL KRIPAD(MELEMS,IPADS)
 
c      nls=MELEMS.NUM(/2)
c      write(6,*)'MELEMS='
c     write(6,1001)(MELEMS.NUM(1,ii),ii=1,nls)
 
c      write(6,*)' I P A D S '
c     npds=ipads.lect(/1)
c     write(6,1001)(ipads.lect(ii),ii=1,npds)
 
       NC=IZTU1.VPOCHA(/2)
       NPTI=IZTU1.VPOCHA(/1)
 
      SEGACT MELEMX
      NBSOUS=MELEMX.LISOUS(/1)
      IF(NBSOUS.EQ.0)NBSOUS=1
      NUTOEL=0
 
      IES=IDIM
 
      SEGACT,MCOORD
      DO 19 L=1,NBSOUS
      IPT1=MELEMX
      IF(NBSOUS.NE.1)IPT1=MELEMX.LISOUS(L)
      SEGACT IPT1
 
      NP  =IPT1.NUM(/1)
      NBEL=IPT1.NUM(/2)
 
      IF(KPREFPU.EQ.0)THEN
      JG=NBEL*NP
      SEGINI MLREE1,MLREE3
      JG=NBEL*NP*IDIM
      SEGINI MLREE2
      CALL ECMO(MTABP,'YP','LISTREEL',MLREE1)
      CALL ECMO(MTABP,'AN','LISTREEL',MLREE2)
      CALL ECMO(MTABP,'MUEF','LISTREEL',MLREE3)
      ELSE
      TYPE=' '
      CALL ACMO(MTABP,'YP',TYPE,MLREE1)
      TYPE=' '
      CALL ACMO(MTABP,'AN',TYPE,MLREE2)
      TYPE=' '
      CALL ACMO(MTABP,'MUEF',TYPE,MLREE3)
      SEGACT MLREE1,MLREE2,MLREE3
      ENDIF
 
      CALL LEKTAB(MTABZ,'XXDIAGSI',MCHPOI)
      CALL LICHT(MCHPOI,IZD,TYPC,IGEOM)
      N=IZD.VPOCHA(/1)
      NC=IZD.VPOCHA(/2)
      SEGINI IZDS
 
      IZGG2=IZDS
      IF(NBINC.EQ.3)THEN
c     write(6,*)' On met les conditions limites sur K et epsilon'
      NOM4(4)=LINCO.MOTS(2)
      NOM4(5)=LINCO.MOTS(3)
      TYPE='SOMMET'
      NC=2
      CALL KRCHPT(TYPE,MELEMS,NC,1,IZG2,NOM4(4))
c     CALL KRCHPT(TYPE,MELEMX,NC,2,IZG2,NOM4(4))
      CALL LICHT(IZG2,IZGG2,TYPC,IGEOM0)
      CALL INITD(IZGG2.VPOCHA,(NC*NPTS),1.D-5)
      ENDIF
 
      CALL FPTAU(IDIM,NP,NBEL,IPT1.NUM,IPADT.LECT,XCOOR,
     & IZTU1.VPOCHA,NPTI,NC,IPADS.LECT,NPTS,
     & MZNT.VPOCHA,MLREE3.PROG,MLREE1.PROG,MLREE2.PROG,MZRO.VPOCHA,IKR,
     & MZMU.VPOCHA,IKM,MZUE.VPOCHA,MZYP.VPOCHA,IPT1.ITYPEL,
     & IZGG2.VPOCHA,NBINC)
 
 19   CONTINUE
      SEGDES,MCOORD
 
c     CALL ECROBJ('MAILLAGE',MELEMX)
c     CALL PRTRAC
c     CALL LIROBJ('MAILLAGE',MELEMX,1,IRET)
      NC=IZTU1.VPOCHA(/2)
      NPTI=IZTU1.VPOCHA(/1)
      TYPE='SOMMET'
      CALL KRCHPT(TYPE,MELEMS,NC,2,IZG1,NOM4)
      CALL LICHT(IZG1,IZGG1,TYPC,IGEOM)
      CALL ECMO(KINC,'FFF','CHPOINT',IZG1)
 
C --> Vitesse seulement on l'impose a 0.
 
      TYPE='SOMMET'
      NC=IDIM
      CALL KRCHPT(TYPE,MELEMS,NC,1,IZGV,NOM4(1))
      CALL LICHT(IZGV,IZGGV,TYPC,IGEOM0)
 
      TYPE=' '
      CALL ACMO(MTAB1,'CLIM',TYPE,KCLIM)
c     write(6,*)' Y F P U :: KCLIM = ',KCLIM
 
      IF(KCLIM.EQ.0)THEN
      CALL ECMO(MTAB1,'CLIM','CHPOINT',IZGV)
 
      ELSE
c On remet à 0  les conditions limites précédentes
      MCHPO4=KCLIM
      SEGACT MCHPO4
      NSP1=MCHPO4.IPCHP(/1)
 
      DO 882 L=1,NSP1
      MSOUPO=MCHPO4.IPCHP(L)
      SEGACT MSOUPO
      NC=NOCOMP(/2)
      MELEKE=IGEOC
      MPOVA4=IPOVAL
      SEGACT MELEKE,MPOVA4*MOD
      NBPKE=MELEKE.NUM(/2)
      DO 883 N=1,NC
      IF(NOCOMP(N).NE.NOM4(1).AND.NOCOMP(N).NE.NOM4(2)
     &                       .AND.NOCOMP(N).NE.NOM4(3))GO TO 883
      DO 884 I=1,NBPKE
      NKE=IPADS.LECT(MELEKE.NUM(1,I))
      IF(NKE.NE.0)MPOVA4.VPOCHA(I,N)=0.D0
 884  CONTINUE
 883  CONTINUE
      SEGDES MSOUPO,MPOVA4,MELEKE
 882  CONTINUE
c      write(6,*)' On remet les CL sur UN CLIM=',MCHPOI
 
      CALL ADCHPO(KCLIM,IZGV,MCHPOI,1D0,1D0)
      CALL ECMO(MTAB1,'CLIM','CHPOINT',MCHPOI)
 
      ENDIF
 
C --> FIN  Vitesse seulement on l'impose a 0.
 
      IZGG2=IZGG1
      IF(NBINC.EQ.3)THEN
C --> On impose K et Epsilon
c On remet à 0 la condition limite précédente pour NOM4(4) et NOM4(5)
      TYPE=' '
      CALL ACMO(MTAB1,'CLIM',TYPE,KCLIM)
 
      MCHPO4=KCLIM
      SEGACT MCHPO4
      NSP1=MCHPO4.IPCHP(/1)
 
      DO 982 L=1,NSP1
      MSOUPO=MCHPO4.IPCHP(L)
      SEGACT MSOUPO
      NC=NOCOMP(/2)
      MELEKE=IGEOC
      MPOVA4=IPOVAL
      SEGACT MELEKE,MPOVA4*MOD
      NBPKE=MELEKE.NUM(/2)
      DO 983 N=1,NC
      IF(NOCOMP(N).NE.NOM4(4).AND.NOCOMP(N).NE.NOM4(5))GO TO 983
      DO 984 I=1,NBPKE
      NKE=IPADS.LECT(MELEKE.NUM(1,I))
      IF(NKE.NE.0)MPOVA4.VPOCHA(I,N)=0.D0
 984  CONTINUE
 983  CONTINUE
      SEGDES MSOUPO,MPOVA4,MELEKE
 982  CONTINUE
 
      CALL ADCHPO(KCLIM,IZG2,MCHPOI,1D0,1D0)
      CALL ECMO(MTAB1,'CLIM','CHPOINT',MCHPOI)
c      write(6,*)' On remet les CL sur KN et EN CLIM=',MCHPOI
 
c     CALL ECROBJ('CHPOINT',MCHPOI)
c     CALL PRLIST
      ENDIF
 
      SEGSUP IZDS
      IF(IKRS.EQ.0)SEGSUP MZRO
      SEGDES IZTU1
      SEGDES IZG1,IZGG1
      IF(NBINC.EQ.3)SEGDES IZG2,IZGG2
      SEGDES IZVOL
      SEGDES LINCO
      SEGSUP MLENTI,IPADS
 
c On crée un MATRIK vide
      NRIGE=7
      NKID =9
      NKMT =7
      NMATRI=0
      SEGINI MATRIK
 
      CALL ECROBJ('MATRIK',MATRIK)
      CALL ECROBJ('CHPOINT',IZG1)
 
c     write(6,*)' Fin FPU Syntaxe II'
      RETURN
c################# IKEL = 1 => Syntaxe II ## FIN #######################
c#######################################################################
      ELSE
 
      IKEL=0
C 5ème coefficient : yp
      IXV(1)=0
      IXV(2)=1
      IXV(3)=0
      IRET=0
      CALL LEKCOF('Opérateur FPU :',
     & MTABX,KINC,5,IXV,MYP,MZYP,NPT3,NC3,IK3,IRET)
      IF(IRET.EQ.0)RETURN
 
      ENDIF
      ENDIF
 
C._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
 
      NC=IZTU1.VPOCHA(/2)
      NPTI=IZTU1.VPOCHA(/1)
      TYPE='SOMMET'
      CALL KRCHPT(TYPE,MELEMS,NC,2,IZG1,NOM4)
      CALL LICHT(IZG1,IZGG1,TYPC,IGEOM)
      CALL ECMO(KINC,'FFF','CHPOINT',IZG1)
 
      IZGG2=IZGG1
      IF(NBINC.GT.1)THEN
C --> 2 eme Inconnue
 
      NOM4(4)=LINCO.MOTS(2)
      NOM4(5)=LINCO.MOTS(3)
      TYPE='SOMMET'
      NC=2
      CALL KRCHPT(TYPE,MELEMS,NC,1,IZG2,NOM4(4))
      CALL LICHT(IZG2,IZGG2,TYPC,IGEOM0)
 
      TYPE=' '
      CALL ACMO(MTAB1,'CLIM',TYPE,KCLIM)
      IF(KCLIM.NE.0)THEN
 
      MCHPO4=KCLIM
      SEGACT MCHPO4
      NSP1=MCHPO4.IPCHP(/1)
 
      DO 782 L=1,NSP1
      MSOUPO=MCHPO4.IPCHP(L)
      SEGACT MSOUPO
      NC=NOCOMP(/2)
      MELEKE=IGEOC
      MPOVA4=IPOVAL
      SEGACT MELEKE,MPOVA4*MOD
      NBPKE=MELEKE.NUM(/2)
      DO 783 N=1,NC
      IF(NOCOMP(N).NE.NOM4(4).AND.NOCOMP(N).NE.NOM4(5))GO TO 783
      DO 784 I=1,NBPKE
      NKE=IPADS.LECT(MELEKE.NUM(1,I))
      IF(NKE.NE.0)MPOVA4.VPOCHA(I,N)=0.D0
 784  CONTINUE
 783  CONTINUE
      SEGDES MSOUPO,MPOVA4,MELEKE
 782  CONTINUE
      ENDIF
      SEGDES MCHPO4
      ENDIF
 
      CALL LEKTAB(MTABZ,'XXDIAGSI',MCHPOI)
      CALL LICHT(MCHPOI,IZD,TYPC,IGEOM)
      N=IZD.VPOCHA(/1)
      NC=IZD.VPOCHA(/2)
      SEGINI IZDS
      IF(KFORM.EQ.1.AND.NBINC.EQ.3)THEN
      SEGINI IZDE,IZDK
      ENDIF
 
 
      SEGACT MELEME
      NBSOUS=LISOUS(/1)
      IF(NBSOUS.EQ.0)NBSOUS=1
      NUTOEL=0
 
      IES=IDIM
 
      DO 1 L=1,NBSOUS
      IPT1=MELEME
      IF(NBSOUS.NE.1)IPT1=LISOUS(L)
      SEGACT IPT1
 
      NOM0=NOMS(IPT1.ITYPEL)//'    '
      CALL KALPBG(NOM0,'FONFORM ',IZFFM)
      SEGACT IZFFM*MOD
      IZHR=KZHR(1)
      SEGACT IZHR*MOD
      IZH2=KZHR(2)
      IZF1=KTP(1)
      NES=GR(/1)
      NPG=GR(/3)
      NP  =IPT1.NUM(/1)
      NBEL=IPT1.NUM(/2)
 
 
 
C SUBROUTINE XCVFPU(NEL,K0,NP,IES,IAXI,IPADL,
C    & LEF,XYZ,               ----> IPT1,COOR
C    & VOLF,                  ----> IZVOL.T,
C    & UN,TK,TE,              ----> IZTU1.T,IZTU2.T,IZTU3.T,
C    & F,                     ----> IZG1,
C    & DK,DE,                 ----> IZD2,IZD3
C    & ANU,IKC,UET,YP,        ----> IZTG1.T,IKM,IZTG2.T,IZTG3.T,
C    & VPAROI,IKV,                  IZTG4.T,IK4,
C    & PORO,NPR,IPOR)         ----> IZPORO,NPOR,IOP7
       if(izde.eq.0) izde=izds
       if(izdk.eq.0) izdk=izds
      SEGACT,MCOORD
      IF(KFORM.EQ.0)THEN
 
      CALL YCVFPU(NBEL,NUTOEL,NP,IES,IAXI,NPTI,LECT,
     & NPTS,IPADS.LECT,
     & IPT1.NUM,XCOOR,
     & IZVOL.VPOCHA,
     & IZTU1.VPOCHA,IZGG2.VPOCHA(1,1),IZGG2.VPOCHA(1,2),
     & IZGG1.VPOCHA,
     & MZNU.VPOCHA,IKM,MZUE.VPOCHA,MZYP.VPOCHA)
      SEGACT,MCOORD
 
      ELSEIF(KFORM.EQ.1.AND.IKEL.EQ.0)THEN
 
      CALL XCVFPU(FN,GR,PG,XYZ,HR,PGSQ,RPG,
     & NES,IDIM,NP,NPG,IAXI,AJ,NBEL,IPT1.NUM,XCOOR,
     & NPTI,LECT,IZTU1.VPOCHA,NPTS,IPADS.LECT,MZUE.VPOCHA,
     & IZGG1.VPOCHA,MZMU.VPOCHA,IKM,MZRO.VPOCHA,IKR,
     & MZYP.VPOCHA,IZDS.VPOCHA,IZDK.VPOCHA,IZDE.VPOCHA,NBINC)
 
      ENDIF
      SEGDES,MCOORD
 
      SEGSUP IZFFM,IZF1,IZHR,IZH2
      NUTOEL=NUTOEL+NBEL
 
 1    CONTINUE
 
      IF(KFORM.EQ.1)THEN
 
      IF(IKEL.EQ.0)THEN
      DO 124 I=1,NPTS
      MZUE.VPOCHA(I,1)=IZDS.VPOCHA(I,1)/IZD.VPOCHA(I,1)
 124  CONTINUE
      ENDIF
 
      IF(NBINC.EQ.3)THEN
      DO 125 I=1,NPTS
      IZGG2.VPOCHA(I,1)=IZDK.VPOCHA(I,1)/IZD.VPOCHA(I,1)
      IZGG2.VPOCHA(I,2)=IZDE.VPOCHA(I,1)/IZD.VPOCHA(I,1)
 125  CONTINUE
      ENDIF
 
      ENDIF
 
      SEGSUP IZDS
      IF(KFORM.EQ.1.AND.NBINC.EQ.3)THEN
      SEGSUP IZDE,IZDK
      ENDIF
      IF(IKRS.EQ.0)SEGSUP MZRO
      SEGDES IZTU1
      SEGDES IZG1,IZGG1
      IF(NBINC.EQ.3)SEGDES IZG2,IZGG2
      SEGDES IZVOL
      SEGDES LINCO
      SEGSUP MLENTI,IPADS
 
      NRIGE=7
      NKID =9
      NKMT =7
      NMATRI=0
      SEGINI MATRIK
 
      IF(NBINC.GT.1)THEN
      IF(KCLIM.EQ.0)THEN
      CALL ECMO(MTAB1,'CLIM','CHPOINT',IZG2)
      ELSE
      CALL ADCHPO(KCLIM,IZG2,MCHPOI,1D0,1D0)
      CALL ECMO(MTAB1,'CLIM','CHPOINT',MCHPOI)
      ENDIF
      ENDIF
 
      CALL ECROBJ('MATRIK',MATRIK)
      CALL ECROBJ('CHPOINT',IZG1)
 
c     write(6,*)' Fin FPU'
      RETURN
 90   CONTINUE
      WRITE(6,*)' Interuption anormale de FPU '
      RETURN
 1001 FORMAT(20(1X,I5))
 1002 FORMAT(10(1X,1PE11.4))
      END
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales