Sources Esope | Cast3M

kfmm2
C KFMM2     SOURCE    CB215821  20/11/25    13:31:21     10792          
      SUBROUTINE KFMM2(IRN,IGN,IRETN,IGAMN,
     &     ICHPSU,ICHPVO,INORM,
     &     MELEMC,MELEMF,MELEFE,
     &     MELLIM,ICHLIM,
     &     ICHRES)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  KFMM2
C
C DESCRIPTION       :  Voir KFMM
C
C                      Cas deux dimensions, gaz "calorically perfect"
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C
C AUTEUR            :  A. BECCANTINI, DRN/DMT/SEMT/TTMF
C
C************************************************************************
C ENTREES
C
C
C
C  1) Pointeurs des CHPOINTs
C
C     IRN     : CHPOINT 'CENTRE' contenant la masse volumique
C
C     IGN     : CHPOINT 'CENTRE' contenant la q.d.m.
C
C     IRETN   : CHPOINT 'CENTRE' contenant l'energie totale
C
C     IGAMN   : CHPOINT 'CENTRE' contenant le gamma
C
C     ICHLIM  : CHPOINT conditions aux bords
C
C  3) Pointeurs de CHPOINTs de la table DOMAINE
C
C     ICHPSU  : CHPOINT "FACE" contenant la surface des faces
C
C     ICHPVO  : CHPOINT "CENTRE" contenant le volume
C               de chaque element
C
C
C  4) Pointeurs de MELEME de la table DOMAINE
C
C     MELEMC  : MELEME 'CENTRE' du SPG des CENTRES
C
C     MELEMF  : MELEME 'FACE' du SPG des FACES
C
C     MELEFE  : MELEME 'FACEL' du connectivité FACES -> ELEM
C
C     MELLIM  : MELEME SPG conditions aux bords
C
C SORTIES
C
C     ICHRES  : resultat
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE : Avril 2002  : creation
C              Janvier 2003: implementation de condition aux limites
C
C************************************************************************
C
C
C N.B.: On suppose qu'on a déjà controllé RO, P > 0
C                                         GAMMA \in (1,3)
C                                         Y \in (0,1)
C       Si non il faut le faire!!!
C
C************************************************************************
C
C
C**** Variables de COOPTIO
C
C      INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
C     &        ,IOTER,   IOLEC,  IOIMP,     IOCAR,  IOACQ
C     &        ,IOPER,   IOSGB,  IOGRA,     IOSAU,  IORES
C     &        ,IECHO,   IIMPI,  IOSPI
C     &        ,IDIM
C     &        ,MCOORD
C     &        ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
C     &        ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
C     &        ,NORINC,NORVAL,NORIND,NORVAD
C     &        ,NUCROU, IPSAUV, IFICLE, IPREFI
C
      IMPLICIT INTEGER(I-N)
      INTEGER IRN,IGN,IRETN,IGAMN,ICHPSU,ICHPVO,INORM
     &     ,MELEMC,MELEMF,MELEFE,ICHRES,ISPG1,ISPG2,NFAC
     &     ,IGEOMF,IGEOMC, NGCEG, NGCED, NGCF, NLCF, NLCF1, NLCEG, NLCED
     &     ,MELLIM,ICHLIM,NLLIM
      REAL*8 ROG, RUXG, RUYG, UNG, RETG, GAMG, REC, PG, VOLG, UTG
     &     , ROD, RUXD, RUYD, UND, RETD, GAMD, PD, VOLD, UTD
     &     , CNX, CNY, CTX, CTY
     &     , SURF, SIGMAF
     &     , UXG, UYG, RHTG
     &     , UXD, UYD, RHTD, COEF
     &     ,FR, FRUX, FRUY, FRET
 
      CHARACTER*8 TYPE
C
C**** LES INCLUDES
C
 
-INC PPARAM
-INC CCOPTIO
-INC SMCHPOI
      POINTEUR MPOVSU.MPOVAL
     &        , MPORES.MPOVAL, MPOVOL.MPOVAL
     &       , MPRN.MPOVAL, MPGN.MPOVAL, MPRETN.MPOVAL, MPGAMN.MPOVAL
     &       , MPNORM.MPOVAL, MPLIM.MPOVAL
-INC SMELEME
-INC SMLMOTS
-INC SMLENTI
      POINTEUR MLELIM.MLENTI
C
C**** Initialisation des MLENTI des conditions aux limites
C
C
      CALL KRIPAD(MELLIM,MLELIM)
C     SEGINI MLELIM
C
C**** Initialisation des MELEMEs
C
C     'CENTRE', 'FACEL'
C
      IPT2 = MELEFE
      SEGACT IPT2
      NFAC = IPT2.NUM(/2)
C
C**** KRIPAD pour la correspondance global/local de centre
C
      CALL KRIPAD(MELEMC,MLENT1)
C
C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements
C
C     Si i est le numero global d'un noeud de ICEN,
C     MLENT1.LECT(i) contient sa position, i.e.
C
C     I              = numero global du noeud centre
C     MLENT1.LECT(i) = numero local du noeud centre
C
C     MLENT1 déjà activé, i.e.
C
C     SEGINI MLENT1
C
C
C**** KRIPAD pour la correspondance global/local de 'FACE'
C
      CALL KRIPAD(MELEMF,MLENT2)
C     SEGINI MLENT2
C
C**** CHPOINTs de la table DOMAINE
C
      CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
      CALL LICHT(ICHPVO,MPOVOL,TYPE,IGEOMC)
      CALL LICHT(INORM,MPNORM,TYPE,IGEOMC)
C
C**** LICHT active les MPOVALs en *MOD
C
C     i.e.
C
C     SEGACT MPOVSU*MOD
C     SEGACT MPOVOL*MOD
C     SEGACT MPNORM*MOD
C
      CALL LICHT(ICHRES,MPORES,TYPE,IGEOMC)
C     SEGACT MPORES*MOD
C
C
      CALL LICHT(IRN,MPRN,TYPE,IGEOMC)
      CALL LICHT(IGN,MPGN,TYPE,IGEOMC)
      CALL LICHT(IRETN,MPRETN,TYPE,IGEOMC)
      CALL LICHT(IGAMN,MPGAMN,TYPE,IGEOMC)
      CALL LICHT(ICHLIM,MPLIM,TYPE,MELLIM)
C
C     SEGACT MPRN*MOD
C     SEGACT MPGN*MOD
C     SEGACT MPRETN*MOD
C     SEGACT MPGAMN*MOD
C     SEGACT MPLIM*MOD
C
      IF(IGEOMF .NE. MELEMF)THEN
         WRITE(IOIMP,*) 'Anomalie dedans kfmm1.eso'
         WRITE(IOIMP,*) 'Probleme de SPG'
C        21 2
C        Données incompatibles
         CALL ERREUR(21)
         GOTO 9999
      ENDIF
      IF(IGEOMC .NE. MELEMC)THEN
         WRITE(IOIMP,*) 'Anomalie dedans kfmm1.eso'
         WRITE(IOIMP,*) 'Probleme de SPG'
C        21 2
C        Données incompatibles
         CALL ERREUR(21)
         GOTO 9999
      ENDIF
C
C
C**** BOUCLE SUR FACEL pour le calcul du FLUX
C
      DO NLCF = 1, NFAC
C
C******* NLCF  = numero local du centre de facel
C        NGCF  = numero global du centre de facel
C        NLCF1 = numero local du centre de face
C        NGCEG = numero global du centre ELT "gauche"
C        NLCEG = numero local du centre ELT "gauche"
C        NGCED = numero global du centre ELT "droite"
C        NLCED = numero local du centre ELT "droite"
C
         NGCEG = IPT2.NUM(1,NLCF)
         NGCED = IPT2.NUM(3,NLCF)
         NGCF  = IPT2.NUM(2,NLCF)
         NLCF1 = MLENT2.LECT(NGCF)
         NLCEG = MLENT1.LECT(NGCEG)
         NLCED = MLENT1.LECT(NGCED)
C
C******* NLCF != NLCF1 -> l'auteur (MOI) n'a rien compris.
C
         IF(NLCF .NE. NLCF1)THEN
            WRITE(IOIMP,*) 'Anomalie dedans kfmm1.eso'
            WRITE(IOIMP,*) 'Probleme de SPG'
C           21 2
C           Données incompatibles
            CALL ERREUR(21)
            GOTO 9999
         ENDIF
C
         CNX = MPNORM.VPOCHA(NLCF,1)
         CNY = MPNORM.VPOCHA(NLCF,2)
         CTX = -1.0D0 * CNY
         CTY = CNX
C
C
C******* Recuperation des Etats "gauche" et "droite"
C
         ROG = MPRN.VPOCHA(NLCEG,1)
         RUXG = MPGN.VPOCHA(NLCEG,1)
         RUYG = MPGN.VPOCHA(NLCEG,2)
         UXG = RUXG / ROG
         UYG = RUYG / ROG
         UNG = (UXG * CNX) + (UYG * CNY)
         UTG = (UXG * CTX) + (UYG * CTY)
         RETG = MPRETN.VPOCHA(NLCEG,1)
         GAMG = MPGAMN.VPOCHA(NLCEG,1)
         REC = 0.5D0 * ((RUXG * RUXG) + (RUYG*RUYG))
         REC = REC / ROG
         PG = (GAMG - 1.0D0) * (RETG - REC)
         VOLG = MPOVOL.VPOCHA(NLCEG,1)
C
         IF(NLCEG .NE. NLCED)THEN
            ROD = MPRN.VPOCHA(NLCED,1)
            RUXD = MPGN.VPOCHA(NLCED,1)
            RUYD = MPGN.VPOCHA(NLCED,2)
            UXD = RUXD / ROD
            UYD = RUYD / ROD
            UND = (UXD * CNX) + (UYD * CNY)
            RETD = MPRETN.VPOCHA(NLCED,1)
            GAMD = MPGAMN.VPOCHA(NLCED,1)
            REC = 0.5D0 * ((RUXD * RUXD) + (RUYD*RUYD))
            REC = REC / ROD
            PD = (GAMD - 1.0D0) * (RETD - REC)
            VOLD = MPOVOL.VPOCHA(NLCED,1)
         ELSE
C           Murs au condition aux limite
            NLLIM=MLELIM.LECT(NGCF)
            IF(NLLIM .EQ.0)THEN
C              Mur
               ROD = ROG
               UND = -1.0D0 * UNG
               UTD = UTG
               UXD = UND * CNX + UTD * CTX
               UYD = UND * CNY + UTD * CTY
               RUXD = UXD * ROD
               RUYD = UYD * ROD
               GAMD = GAMG
               RETD = RETG
               PD = PG
               VOLD = VOLG
            ELSE
C              Boundary
               ROD = MPLIM.VPOCHA(NLLIM,1)
               UXD = MPLIM.VPOCHA(NLLIM,2)
               UYD = MPLIM.VPOCHA(NLLIM,3)
               PD = MPLIM.VPOCHA(NLLIM,4)
               GAMD = GAMG
               UND = (UXD * CNX) + (UYD * CNY)
               RETD = ((1.0D0/(GAMD - 1.0D0))*PD)+
     &         (0.5D0*ROD*((UXD*UXD)+(UYD*UYD)))
               VOLD = VOLG
            ENDIF
         ENDIF
C
         RHTG=RETG+PG
         RHTD=RETD+PD
C
         SURF   = MPOVSU.VPOCHA(NLCF,1)
         SIGMAF = (0.5D0 * (GAMG + GAMD)) * (PG + PD) / (ROG + ROD)
         SIGMAF = SIGMAF ** 0.5D0
         SIGMAF = (0.5D0 * (ABS(UNG) + ABS(UND))) + SIGMAF
C
         FR=(ROG*UNG)+(ROD*UND)
         FRUX=(ROG*UNG*UXG)+(ROD*UND*UXD)+((PG+PD)*CNX)
         FRUY=(ROG*UNG*UYG)+(ROD*UND*UYD)+((PG+PD)*CNY)
         FRET=(UNG*RHTG)+(UND*RHTD)
C
         COEF=SURF/(2.0D0*VOLG)
         MPORES.VPOCHA(NLCEG,1) = MPORES.VPOCHA(NLCEG,1) +
     &        (COEF * (FR -  (SIGMAF * ROD)))
         MPORES.VPOCHA(NLCEG,2) = MPORES.VPOCHA(NLCEG,2) +
     &        (COEF * (FRUX - (SIGMAF * ROD * UXD)))
         MPORES.VPOCHA(NLCEG,3) = MPORES.VPOCHA(NLCEG,3) +
     &        (COEF *  (FRUY - (SIGMAF * ROD * UYD)))
         MPORES.VPOCHA(NLCEG,4) = MPORES.VPOCHA(NLCEG,4) +
     &        (COEF *  (FRET - (SIGMAF * RETD)))
C
         IF(NLCEG .NE. NLCED)THEN
C
C********** Domain interieur
C
            COEF=SURF/(2.0D0*VOLD)
            MPORES.VPOCHA(NLCED,1) = MPORES.VPOCHA(NLCED,1) -
     &           (COEF * (FR +  (SIGMAF * ROG)))
            MPORES.VPOCHA(NLCED,2) = MPORES.VPOCHA(NLCED,2) -
     &           (COEF * (FRUX + (SIGMAF * ROG * UXG)))
            MPORES.VPOCHA(NLCED,3) = MPORES.VPOCHA(NLCED,3) -
     &           (COEF *  (FRUY + (SIGMAF * ROG * UYG)))
            MPORES.VPOCHA(NLCED,4) = MPORES.VPOCHA(NLCED,4) -
     &           (COEF *  (FRET + (SIGMAF * RETG)))
         ENDIF
      ENDDO
C
      SEGDES MPOVSU
      SEGDES MPOVOL
      SEGDES MPNORM
C
      SEGDES MPORES
C
      SEGDES MPRN
      SEGDES MPRETN
      SEGDES MPGN
      SEGDES MPGAMN
C
      SEGSUP MLENT1
      SEGSUP MLENT2
      SEGDES IPT2
C
      SEGSUP MLELIM
      IF(MPLIM .GT. 0) SEGDES MPLIM
C
C
 9999 CONTINUE
C
      RETURN
      END
C