Sources Esope | Cast3M

kon19
C KON19     SOURCE    CB215821  20/11/25    13:31:57     10792          
      SUBROUTINE KON19
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  KON19
C
C DESCRIPTION       :  Subroutine appellée par KON2
C
C                      Modelisation 2D/3D des equations d'Euler
C                       MULTIESPECES!!!!!!!!!!!!!!!
C
C                      Calcul du jacobien
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C
C AUTEUR            :  S. KUDRIAKOV, DRN/DMT/SEMT/LTMF
C
C*********************************************************************
C
C APPELES (Calcul)  :  KONMSP (calcul du jacobien, gaz "calorically
C                              perfect",multiespeces, 2D, AUSM+)
C                   :  KONMS3 (calcul du jacobien, gaz "calorically
C                              perfect",multiespeces, 3D, AUSM+)
C*********************************************************************
C
C*** SYNTAXE
C
C  Discrétisation en VF "cell-centered" des équations d'Euler pour
C  un gaz parfait mono-constituent polytropique
C  Inconnues: densité, quantité de mouvement, énergie totale par
C  unité de volumes (variables conservatives)
C
C  RMAT1 = 'KONV' 'VF' 'PERFMULT' 'JACOCONV' MOT3 TABD TABP LMOT1
C          (MAILIM) CHPO1 CHPO2 CHPO3 CHPO4 ;
C-----------------------
C  ENTREES
C-----------------------
C  MOT3   : objet de type MOT
C           'AUSMPLUS' : jacobien du residu pour la methode AUSM+
C
C  TABD   : la table domaine
C
C  TABP    : TABLE qui contient :
C              * les noms des especes qui apparessent
C                explicitement dans les equations d'Euler en
C                TABP . 'ESPEULE' (list de mots);
C              * le nom de l'espece qui n'y est pas (mots);
C              * les CP et les CV du gas en
C                TABP . 'CP' (table)
C                TABP . 'CV' (table)
C
C  LMOT1   : objet de type LISTMOTS
C           Noms de composantes des variable primales et duales de RMAT1.
C           Il contient dans l'ordre suivant: le noms de la densité,
C           du momentum, de l'énergie totale et les fraction massiques
C           par unité de volume
C
C  MAILIM : MAIILAGE de POI1 ou on ne veut pas calculer le FLUX
C           convective
C
C  CHPO1  : CHPOINT contenant la masse volumique
C           (SPG = TAB1 . 'CENTRE', une seule composante,
C           'SCAL').
C
C  CHPO2  : CHPOINT contenant la vitesse
C           (SPG = TAB1 . 'CENTRE', deux/trois composantes
C           'UX', 'UY', 'UZ')
C
C  CHPO3  : CHPOINT contenant la pression du gaz
C           (SPG = TAB1 . 'CENTRE', une seule composante,
C            'SCAL').
C
C  CHPO4  : CHPOINT contenant les fraction massiques
C           des differentes especes
C           (SPG = TAB1 . 'CENTRE', nombre des composantes = NESP-1).
C
C---------------------
C SORTIES
C---------------------
C  RMAT1  : objet de type MATRIK
C          (SPG =  TAB1 . 'CENTRE')
C          (inconnues primales = inconnues duales = LMOT1)
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE :
C
C************************************************************************
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
 
-INC PPARAM
-INC CCOPTIO
-INC SMLMOTS
-INC SMCHPOI
-INC SMELEME
 
      POINTEUR MLMVIT.MLMOTS
C
      INTEGER NBJAC, IRET, INDIC, NBCOMP, NESP, JGN, JGM
     &     ,IDOMA, MELEMC, MELEMF, MELEFE, MELTFA, ICHPSU, ICHPDI
     &     ,ICHPVO, INORM, ICP, MMODEL
     &     ,IJACO, ILIINC, NC, IFLIM, MELLIM, ICACCA
     &     ,IIMPL, IRN, IVN, IPN, IYN, IUINF, IUPRI
     &     ,I1,  JG, ICOND, IRETOU, ICV, INEFMD, IPGAS
C
      PARAMETER (NBJAC=2)
      CHARACTER*8 TYPE, LJACO(NBJAC), MTYPR
      CHARACTER*4 MOT
      CHARACTER*(4) MOT1(3), CELLCH
      CHARACTER*(40) MESERR
      REAL*8 CP, CV
C----------------------------------------------------
C************ Variables en ACCTAB  ******************
C----------------------------------------------------
      INTEGER IVALI, IRETI,IVALR, IRETR
      REAL*8 XVALI,XVALR
      LOGICAL LOGII, LOGIR
      CHARACTER*(8) CHARR,MTYPI
      DATA LJACO/'AUSMPLUS','AUSMPLM '/
cC----------------------------------------------------
cC************ Variables en ACCTAB  ******************
cC----------------------------------------------------
c      INTEGER IVALI, IRETI,IVALR, IRETR
c      REAL*8 XVALI,XVALR
c      LOGICAL LOGII, LOGIR
c      CHARACTER*(8) CHARR,MTYPI
C----------------------------------------------------
-INC SMLREEL
      POINTEUR MLRECP.MLREEL, MLRECV.MLREEL
C-----------------------------------------
C Initialisation des variables en ACCTAB
C-----------------------------------------
      IVALI = 0
      IVALR = 0
      XVALI = 0.0D0
      XVALR = 0.0D0
      LOGII = .FALSE.
      LOGIR = .FALSE.
      IRETI = 0
      IRETR = 0
      CHARR = '        '
C-----------------------------
C**** Type of Jacobian
C-----------------------------
      CALL LIRMOT(LJACO,NBJAC,IIMPL,1)
      IF(IERR .NE. 0)GOTO 9999
C*******************************
C**** La table domaine *********
C*******************************
      ICOND = 1
      CALL QUETYP(TYPE,ICOND,IRET)
C
      IF(IRET.EQ.0.AND.TYPE.NE.'MMODEL')THEN
        WRITE(IOIMP,*)' On attend un objet MMODEL'
        GOTO 9999
      ENDIF
      CALL LIROBJ('MMODEL',MMODEL,ICOND,IRET)
      IF(IERR.NE.0)GOTO 9999
      CALL LEKMOD(MMODEL,IDOMA,INEFMD)
C     INEFMD inutilisé
      IF(IERR.NE.0)GOTO 9999
C--------------------------------
C**** Centre, FACE, FACEL, ELTFA
C--------------------------------
      CALL LEKTAB(IDOMA,'CENTRE',MELEMC)
      IF(IERR .NE. 0) GOTO 9999
C
      CALL LEKTAB(IDOMA,'FACE',MELEMF)
      IF(IERR .NE. 0) GOTO 9999
C
      CALL LEKTAB(IDOMA,'FACEL',MELEFE)
      IF(IERR .NE. 0) GOTO 9999
C
      CALL LEKTAB(IDOMA,'ELTFA',MELTFA)
      IF(IERR .NE. 0) GOTO 9999
C---------------------------------------------------------
C**** Lecture du CHPOINT contenant les surfaces des faces.
C---------------------------------------------------------
      CALL LEKTAB(IDOMA,'XXSURFAC',ICHPSU)
      IF(IERR .NE. 0) GOTO 9999
      INDIC = 1
      NBCOMP = 1
      MOT = 'SCAL'
      CALL QUEPOI(ICHPSU, MELEMF, INDIC, NBCOMP, MOT)
      IF(IERR .NE. 0) GOTO 9999
C---------------------------------------------------------
C**** Lecture du CHPOINT contenant les diametres minimums.
C---------------------------------------------------------
      CALL LEKTAB(IDOMA,'XXDIEMIN',ICHPDI)
      IF(IERR .NE. 0) GOTO 9999
      INDIC = 1
      NBCOMP = 1
      MOT = 'SCAL'
      CALL QUEPOI(ICHPDI, MELEMC, INDIC, NBCOMP, MOT)
      IF(IERR .NE. 0) GOTO 9999
C----------------------------------------------
C**** Lecture du CHPOINT contenant les volumes
C----------------------------------------------
      CALL LEKTAB(IDOMA,'XXVOLUM',ICHPVO)
      IF(IERR .NE. 0) GOTO 9999
      INDIC = 1
      NBCOMP = 1
      MOT = 'SCAL'
      CALL QUEPOI(ICHPVO, MELEMC, INDIC, NBCOMP, MOT)
      IF(IERR .NE. 0) GOTO 9999
C----------------------------
C**** Les normales aux faces
C----------------------------
      IF(IDIM .EQ. 2)THEN
C        Que les normales
         CALL LEKTAB(IDOMA,'XXNORMAF',INORM)
         IF(IERR .NE. 0) GOTO 9999
         JGN = 4
         JGM = 2
         SEGINI MLMVIT
         MLMVIT.MOTS(1) = 'UX  '
         MLMVIT.MOTS(2) = 'UY  '
         CALL QUEPO1(INORM, MELEMF, MLMVIT)
         SEGSUP MLMVIT
         IF(IERR .NE. 0) GOTO 9999
      ELSE
C        Les normales et les tangentes
         TYPE = ' '
         CALL ACMO(IDOMA,'MATROT',TYPE,INORM)
         IF (TYPE .NE. 'CHPOINT ') THEN
            CALL MATRAN(IDOMA,INORM)
            IF(IERR .NE. 0) GOTO 9999
         ENDIF
         JGN = 4
         JGM = 9
         SEGINI MLMVIT
         MLMVIT.MOTS(1) = 'UX  '
         MLMVIT.MOTS(2) = 'UY  '
         MLMVIT.MOTS(3) = 'UZ  '
         MLMVIT.MOTS(4) = 'RX  '
         MLMVIT.MOTS(5) = 'RY  '
         MLMVIT.MOTS(6) = 'RZ  '
         MLMVIT.MOTS(7) = 'MX  '
         MLMVIT.MOTS(8) = 'MY  '
         MLMVIT.MOTS(9) = 'MZ  '
         CALL QUEPO1(INORM, MELEMF, MLMVIT)
         SEGSUP MLMVIT
      ENDIF
C********************************
C**** Fin table domaine *********
C********************************
      ICOND = 1
      CALL QUETYP(MTYPR,ICOND,IRETOU)
      IF(IERR .NE. 0)GOTO 9999
      IF(MTYPR .NE. 'TABLE   ')THEN
C---------------------------------------
C Message d'erreur standard
C 37 2
C On ne trouve pas d'objet de type %m1:8
C---------------------------------------
         MOTERR(1:8) = 'TABLE   '
         CALL ERREUR(37)
         GOTO 9999
      ELSE
         ICOND = 1
         CALL LIROBJ(MTYPR,IPGAS,ICOND,IRETOU)
         IF(IERR .NE. 0)GOTO 9999
      ENDIF
C-------------------------------------------
C********  Reading the CPs *****************
C-------------------------------------------
      MTYPR = '        '
      CALL ACMO(IPGAS,'CP',MTYPR,ICP)
      IF(MTYPR .NE. 'TABLE   ')THEN
C--------------------------------
C Message d'erreur standard
C -301 0 %m1:40
C--------------------------------
         MOTERR(1:40) = 'TAB1 . CP = ???                     '
         WRITE(IOIMP,*) MOTERR(1:40)
C---------------------------------
C Message d'erreur standard
C 21 2
C Données incompatibles
C--------------------------------
         CALL ERREUR(21)
         GOTO 9999
      ENDIF
C-------------------------------------------
C********  Reading the CVs *****************
C-------------------------------------------
      MTYPR = '        '
      CALL ACMO(IPGAS,'CV',MTYPR,ICV)
      IF(MTYPR .NE. 'TABLE   ')THEN
C--------------------------------
C Message d'erreur standard
C -301 0 %m1:40
C--------------------------------
         MOTERR(1:40) = 'TAB1 . CV = ???                     '
         WRITE(IOIMP,*) MOTERR(1:40)
C---------------------------------
C Message d'erreur standard
C 21 2
C Données incompatibles
C--------------------------------
         CALL ERREUR(21)
         GOTO 9999
      ENDIF
C---------------------------------------------------------
C Reading Les especes qui sont dans les Equations d'Euler
C---------------------------------------------------------
      MTYPR = '        '
      CALL ACMO(IPGAS,'ESPEULE',MTYPR,MLMOT1)
      IF(MTYPR .NE. 'LISTMOTS')THEN
C---------------------------
C Message d'erreur standard
C -301 0 %m1:40
C---------------------------
         MOTERR(1:40) = 'TAB1 . ESPEULE = ???                '
         WRITE(IOIMP,*) MOTERR(1:40)
C---------------------------
C Message d'erreur standard
C 21 2
C Données incompatibles
C--------------------------
         CALL ERREUR(21)
         GOTO 9999
      ENDIF
C-----------------------------------------------------------------
C Reading Nom de l'espece qui n'est pas dans les equations d'Euler
C-----------------------------------------------------------------
      MTYPI = 'MOT     '
      MTYPR = '        '
      CALL ACCTAB(IPGAS,MTYPI,IVALI,XVALI,'ESPNEULE', LOGII,IRETI,
     &                  MTYPR,IVALR,XVALR ,CELLCH,LOGIR,IRETR)
      IF((IERR .NE. 0) .OR. (MTYPR .NE. 'MOT     ')) THEN
C-------------------------------
C Message d'erreur standard
C -301 0 %m1:40
C-------------------------------
         MOTERR = 'TAB1 . ESPNEULE = ???        '
         WRITE(IOIMP,*) MOTERR(1:40)
C------------------------------
C Message d'erreur standard
C 21 2
C Données incompatibles
C------------------------------
         CALL ERREUR(21)
         GOTO 9999
      ENDIF
C--------------------------------------------------------------
C  Creating the lists of real numbers of CP and CV
C--------------------------------------------------------------
      SEGACT MLMOT1
      NESP = MLMOT1.MOTS(/2)
C----------------------
C** List de CP et CV **
C----------------------
      JG = NESP+1
      SEGINI MLRECP
      SEGINI MLRECV
      DO I1 = 1, NESP
C-------------------------------
C N.B. MOT1 est un CHARACTER*(4)
C-------------------------------
         MOT1(1) = MLMOT1.MOTS(I1)
         MTYPI = 'MOT     '
         MTYPR = '        '
         CALL ACCTAB(ICP,MTYPI,IVALI,XVALI,MOT1(1), LOGII,IRETI,
     &               MTYPR,IVALR,CP   ,CHARR,LOGIR,IRETR)
         IF((IERR .NE. 0) .OR. (MTYPR .NE. 'FLOTTANT')) THEN
C------------------------------
C Message d'erreur standard
C -301 0 %m1:40
C------------------------------
            MOTERR = 'TAB1 . CP , TAB1 . ESPEULE = ???        '
            WRITE(IOIMP,*) MOTERR(1:40)
C------------------------------
C Message d'erreur standard
C 21 2
C Données incompatibles
C------------------------------
            CALL ERREUR(21)
            GOTO 9999
         ENDIF
         MLRECP.PROG(I1) = CP
C-------------------------------------------
         CALL ACCTAB(ICV,MTYPI,IVALI,XVALI,MOT1(1), LOGII,IRETI,
     &               MTYPR,IVALR,CV   ,CHARR,LOGIR,IRETR)
         IF((IERR .NE. 0) .OR. (MTYPR .NE. 'FLOTTANT')) THEN
C---------------------------
C Message d'erreur standard
C -301 0 %m1:40
C---------------------------
            MOTERR = 'TAB1  . CV , TAB1 . ESPEULE = ???       '
            WRITE(IOIMP,*) MOTERR(1:40)
C---------------------------
C Message d'erreur standard
C 21 2
C Données incompatibles
C---------------------------
            CALL ERREUR(21)
            GOTO 9999
         ENDIF
         MLRECV.PROG(I1) = CV
      ENDDO
      MTYPI = 'MOT     '
      MTYPR = '        '
      CALL ACCTAB(ICP,MTYPI,IVALI,XVALI,CELLCH, LOGII,IRETI,
     &            MTYPR,IVALR,CP   ,CHARR,LOGIR,IRETR)
      IF((IERR .NE. 0) .OR. (MTYPR .NE. 'FLOTTANT')) THEN
C---------------------------
C Message d'erreur standard
C -301 0 %m1:40
C---------------------------
         MOTERR = 'TAB1  . CP , TAB1 . ESPNEULE = ???      '
         WRITE(IOIMP,*) MOTERR(1:40)
C---------------------------
C Message d'erreur standard
C 21 2
C Données incompatibles
C---------------------------
         CALL ERREUR(21)
         GOTO 9999
      ENDIF
      MLRECP.PROG(JG) = CP
C
      MTYPI = 'MOT     '
      MTYPR = '        '
      CALL ACCTAB(ICV,MTYPI,IVALI,XVALI,CELLCH, LOGII,IRETI,
     &            MTYPR,IVALR,CV   ,CHARR,LOGIR,IRETR)
      IF((IERR .NE. 0) .OR. (MTYPR .NE. 'FLOTTANT')) THEN
C----------------------------
C Message d'erreur standard
C -301 0 %m1:40
C----------------------------
         MOTERR = 'TAB1 . CV , TAB1 . ESPNEULE = ???       '
         WRITE(IOIMP,*) MOTERR(1:40)
C----------------------------
C Message d'erreur standard
C 21 2
C Données incompatibles
C----------------------------
         CALL ERREUR(21)
         GOTO 9999
      ENDIF
      MLRECV.PROG(JG) = CV
C********************************************
C**** Fin de la table de gas ****************
C********************************************
C--------------------------
C**** La list des inconnues
C--------------------------
      TYPE='LISTMOTS'
      CALL LIROBJ(TYPE,ILIINC,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
      MLMOTS = ILIINC
      SEGACT MLMOTS
      NC = MLMOTS.MOTS(/2)
      SEGDES MLMOTS
      IF(NC .NE. (IDIM+2+NESP))THEN
         MOTERR(1:40) = 'LISTINCO = ???'
         WRITE(IOIMP,*) MOTERR
C--------------------------
C Message d'erreur standard
C  21 2
C Données incompatibles
C--------------------------
         CALL ERREUR(21)
         GOTO 9999
      ENDIF
C-------------------------------------
C******** Boundary condition *********
C-------------------------------------
      IRET=0
      TYPE='MAILLAGE'
      CALL LIROBJ(TYPE,IFLIM,0,IRET)
      IF(IERR.NE.0)GOTO 9999
      IF(IRET .EQ. 0)THEN
         MELLIM = 0
      ELSE
         MELEME=IFLIM
         SEGACT MELEME
         ICACCA=MELEME.NUM(/2)
         IF(ICACCA .EQ. 0)THEN
            MELLIM = 0
         ELSE
            MELLIM = IFLIM
         ENDIF
         SEGDES MELEME
      ENDIF
C--------------------------
C**** La densité au centre
C--------------------------
      TYPE = 'CHPOINT '
      CALL LIROBJ(TYPE,IRN,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
      INDIC = 1
      NBCOMP = 1
      MOT = 'SCAL'
      CALL QUEPOI(IRN, MELEMC, INDIC, NBCOMP, MOT)
      IF(IERR .NE. 0) GOTO 9999
C---------------------------
C*** La vitesse au centre
C---------------------------
      TYPE = 'CHPOINT '
      CALL LIROBJ(TYPE,IVN,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
      JGN = 4
      JGM = IDIM
      SEGINI MLMVIT
      MLMVIT.MOTS(1) = 'UX  '
      MLMVIT.MOTS(2) = 'UY  '
      IF(IDIM .EQ. 3) MLMVIT.MOTS(3) = 'UZ  '
      CALL QUEPO1(IVN, MELEMC, MLMVIT)
      SEGSUP MLMVIT
      IF(IERR .NE. 0) GOTO 9999
C----------------------------
C*** La pression au centre
C----------------------------
      TYPE = 'CHPOINT '
      CALL LIROBJ(TYPE,IPN,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
      INDIC = 1
      NBCOMP = 1
      MOT = 'SCAL'
      CALL QUEPOI(IPN, MELEMC, INDIC, NBCOMP, MOT)
      IF(IERR .NE. 0) GOTO 9999
C-----------------------------
C*** Fraction massiques
C*** des especes "splittees"
C-----------------------------
      TYPE = 'CHPOINT '
      CALL LIROBJ(TYPE,IYN,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
      CALL QUEPO1(IYN, MELEMC, MLMOT1)
      IF(IERR .NE. 0) THEN
 
       MOTERR = 'CHPO FR.MAS. = ?????'
       WRITE(IOIMP,*) MOTERR(1:40)
       CALL ERREUR(21)
       GOTO 9999
      ENDIF
C-----------------------------------
C**** Bas Mach
C-----------------------------------
      IF(IIMPL .EQ. 2)THEN
         TYPE = 'CHPOINT '
C------------------
C******* Cut off 1
C------------------
         CALL LIROBJ(TYPE,IUINF,1,IRET)
         IF(IERR .NE. 0) GOTO 9999
         INDIC = 1
         NBCOMP = 1
         MOT = 'SCAL'
         CALL QUEPOI(IUINF, MELEMC, INDIC, NBCOMP, MOT)
         IF(IERR .NE. 0) GOTO 9999
C-----------------
C******* Cut off 2
C-----------------
         TYPE = 'CHPOINT '
         CALL LIROBJ(TYPE,IUPRI,1,IRET)
         IF(IERR .NE. 0) GOTO 9999
         INDIC = 1
         NBCOMP = 1
         MOT = 'SCAL'
         CALL QUEPOI(IUPRI, MELEMC, INDIC, NBCOMP, MOT)
         IF(IERR .NE. 0) GOTO 9999
C
      ELSE
         IUINF=0
         IUPRI=0
      ENDIF
C***************************************
C*****   Fin reading chpoints  *********
C***************************************
      NESP=NESP+1
      IF(IIMPL .EQ. 1)THEN
c------------------------------
c   AUSM+
c-----------------------------
       IF(IDIM .EQ. 2) THEN
        CALL KONMSP(NESP,ILIINC,IRN,IVN,IPN,IYN,MLRECP,
     &    MLRECV,INORM,ICHPVO,ICHPSU,MELEMC,MELEFE,MELLIM,IJACO)
       ELSE
        CALL KONMS3(NESP,ILIINC,IRN,IVN,IPN,IYN,MLRECP,
     &    MLRECV,INORM,ICHPVO,ICHPSU,MELEMC,MELEFE,MELLIM,IJACO)
       ENDIF
C------------------------------------------------------------
      ELSEIF(IIMPL .EQ. 2)THEN
C-----------------------------
C********** AUSM+ low Mach
C-----------------------------
         IF(IDIM .EQ. 2)THEN
          CALL KBMMS2(NESP,ILIINC,IRN,IVN,IPN,IYN,MLRECP,MLRECV,
     $     INORM,ICHPVO,ICHPSU, IUINF, IUPRI, MELEMC, MELEFE,
     $           MELLIM,IJACO)
         ENDIF
      ENDIF
      TYPE='MATRIK '
      CALL ECROBJ(TYPE,IJACO)
 9999 CONTINUE
      RETURN
      END