C KONV2F    SOURCE    MAGN      09/07/30    21:18:32     6458           
      SUBROUTINE KONV2F()
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  KONV2F
C
C DESCRIPTION       :  Subroutine appellée par KONV1
C
C                      Modelling of Two Phase Flow 
C                      (two fluid model, mixtures of water and air)
C
C                      Calcul du residu / jacobien / DELTAT
C     
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C      
C AUTEUR            :  A. BECCANTINI, DRN/DMT/SEMT/LTMF              
C                      Modified for two phase flow by
C                      Jose R. Garcia-Cascales,
C                      Universidad Politecnica de Cartagena,
C                      jr.garcia@upct.es
C************************************************************************
C
C APPELES (Calcul)  :  KONFL8 (2D, two fluid flow)
C                      KONFL9 (3D, two fluid flow)
C
C************************************************************************
C
C*** SYNTAXE
C  FV Discretization of the Two Phase Flow system of equations
C  Unknown: mass flow rate, momentum and total energy per unit volume of
C      each phase (conservative variables)
C 
C
C  RMAT1 RCHPO1 RFLOT1 = 'KONV' 'VF' 'PERFMONO' MOT1 MOT2
C                         TABD LMOT1 
C                         MCHAM1 MCHAM2 MCHAM3 MCHAM4 
C                         MCHAM5 MCHAM6 MCHAM7 MCHAM8
C                         (MOT3 CHPO1 CHPO2 CHPO3 CHPO5) ;
C
C  ENTREES 
C
C
C  MOT1   : objet de type MOT
C           Il vaut 'RESI' si on veut calculer le résidu
C           Il vaut 'FLUX' si on veut calculer le flux
C
C  MOT2   : objet de type MOT
C           Il indique la méthode de décentrement:
C           'AUSMP1' = AUSM +
C           'AUSMP2' = Preconditioned AUSM +
C           'AUSMDV1'= AUSMDV
C           'AUSMDV2'= Preconditioned AUSMDV
C
C LMOT1   : objet de type LISTMOTS
C           Noms de composantes du résultat (RCHPO1)
C           Il contient dans l'ordre suivant: le noms de la densité,
C           de la vitesse, de l'énergie totale par unité de volume
C
C TAB1    : la table domaine
C
C MCHAM1  : MCHAML containing void fraction
C           SPG (support géométrique) l'indice 'FACEL' de 
C           TAB1 (une composante, 'SCAL')
C
C MCHAM2  : MCHAML  vapour velocity and director cosines
C           du repère locale (n,t) dans le repère
C           global (x,y) (dans le cas 2D 6 composantes:
C           * 'UN' = vitesse normale  (SPG = TAB1 . 'FACEL')
C           * 'UT' = vitesse tangentielle  (SPG = 'TAB1 . FACEL')
C           * 'NX' = n.x (SPG = 'FACE')
C           * 'NY' = n.y (SPG = 'FACE')
C           * 'TX' = t.x (SPG = 'FACE')
C           * 'TY' = t.y (SPG = 'FACE')).
C MCHAM3  : MCHAML  liquid velocity and director cosines
C           du repère locale (n,t) dans le repère
C           global (x,y) (dans le cas 2D 6 composantes:
C           * 'UN' = vitesse normale  (SPG = TAB1 . 'FACEL')
C           * 'UT' = vitesse tangentielle  (SPG = 'TAB1 . FACEL')
C           * 'NX' = n.x (SPG = 'FACE')
C           * 'NY' = n.y (SPG = 'FACE')
C           * 'TX' = t.x (SPG = 'FACE')
C           * 'TY' = t.y (SPG = 'FACE')).
C
C MCHAM4  : MCHAML (SPG = TAB1 . 'FACEL') containing pressure
C           une  seule composante, 'SCAL').
C
C MCHAM5  : MCHAML (SPG = TAB1 . 'FACEL') contaning gas/vapour temperature
C           (une seule composante, 'SCAL').
C
C MCHAM6  : MCHAML (SPG = TAB1 . 'FACEL') contaning liquid temperature
C           (une seule composante, 'SCAL').
C
C MCHAM7  : MCHAML (SPG = TAB1 . 'FACEL') contaning gas/vapour density
C           (une seule composante, 'SCAL').
C
C MCHAM8  : MCHAML (SPG = TAB1 . 'FACEL') contaning liquid density
C           (une seule composante, 'SCAL').
C
C RCHPO1  : objet de type CHPOINT (composantes =  LMOT1)
C           Residu  si MOT2 = 'RESI' (SPG = TAB1 . 'CENTRE')
C           Flux si MOT2 = 'FLUX'    (SPG = TAB1 . 'FACE') 
C
C RFLOT1  : objet de type FLOTTANT
C           Il est le temps caracteristique associé à l'onde la plus
C           rapide.
C
C Remarque
C --------
C
C RCHPO1 est égal à:
C * la derivé temporelle des inconnues si l'option 'RESI' est utilisée
C * la projection du flux convectif sur (TAB1 . 'XXNORMAF') si 
C   l'option 'FLUX' est utilisée
C      
C***********************************************************************
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE :
C      
C************************************************************************
C
-INC CCOPTIO
-INC SMLMOTS
-INC SMMATRIK
-INC SMCHPOI
      POINTEUR MLMVIT.MLMOTS,MLMOTY.MLMOTS
C
C**** Variables de COOPTIO
C
      INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
     &        ,IOTER,   IOLEC,  IOIMP,     IOCAR,  IOACQ
     &        ,IOPER,   IOSGB,  IOGRA,     IOSAU,  IORES
     &        ,IECHO,   IIMPI,  IOSPI
     &        ,IDIM,    IFICLE, IPREFI 
     &        ,MCOORD
     &        ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
     &        ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
     &        ,NORINC,NORVAL,NORIND,NORVAD
     &        ,NUCROU, IPSAUV
C
      INTEGER IDOMA, MELEMC, MELEMF, MELEFE, ICHPSU, ICHPDI , ICHPVO
     &      , MELTFA, NBMET, INDMET, IRET, NBJAC
     &      , IAIF, IALPF, IUVF, IULF, IPF, ITVF, ITLF
     &      , IRVF, IRLF, IFRMAF
     &      , ILIINC, NC, ICELL
     &      , INDIC, NBCOMP, IJACO
     &      , NKID, NKMT, NMATRI, NRIGE         
     &      , JGM, JGN, INORM, IIMPL
     &      , ICHFLU, ICHRES

C
      PARAMETER (NBMET=4,NBJAC=5)
      REAL*8 DT
      CHARACTER*8 LMETO(NBMET), TYPE, LJACO(NBJAC)
      CHARACTER*4 MOT, LFLUX(2)
      CHARACTER*(40) MESERR
      LOGICAL LOGNC, LOGAN, LOGRES
C**** Depending on if we want to use AUSM+, AUSM+ with "preconditioning",
C     AUSMDV or AUSMDV with "preconditioning" we have respectively
      DATA LMETO/'AUSMP1','AUSMP2','AUSMDV1','AUSMDV2'/
C**** Option not implemented yet
      DATA LJACO/'VLHJACO1','AUSMJAC1','AUSMLMJ1','VLHJACO2',
     &           'AUSMJAC2'/
C**** If we want flux or residual
      DATA LFLUX/'FLUX','RESI'/
C
C**** Initialisation des variables pour la gestion des erreurs.
C
      LOGNC  = .FALSE.
      LOGAN  = .FALSE.
      MESERR  = '                                        '
C
C******* Flux ou residu???
C
      CALL LIRMOT(LFLUX,2,ICELL,1)
      IF(IERR .NE. 0)GOTO 9999
      IF(ICELL .EQ. 1)THEN
         LOGRES = .FALSE.
      ELSEIF(ICELL .EQ. 2)THEN
         LOGRES = .TRUE.
      ELSE
C
C******** Message d'erreur standard
C         251 2
C         Tentative d'utilisation d'une option non implémentée
C
         CALL ERREUR(251)
      ENDIF
C
C**** Metode utilisée
C
      CALL LIRMOT(LMETO,NBMET,INDMET,1)
      IF(IERR .NE. 0)GOTO 9999
      IF(INDMET .EQ. 0)THEN
C
C******** Message d'erreur standard
C         251 2
C         Tentative d'utilisation d'une option non implémentée
C
         CALL ERREUR(251)
      ENDIF
C     
C*******************************
C**** La table domaine *********
C*******************************
C
      CALL LIROBJ('TABLE',IDOMA,1,IRET)
      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
C     De momento solo 2D
      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********************************
C**** Fin table domaine *********
C********************************
C
C
C**** On va lire les pointeurs des MCHAMLs
C     Reading the MCHAML 'FACEL' containig void fraction
C
      TYPE='MCHAML  '
      CALL  LIROBJ(TYPE,IALPF,1,IRET)
      IF(IERR.NE.0) GOTO 9999
C
C**** Reading the MCHAML 'FACEL' containing vapour velocity
C   
      TYPE='MCHAML  '
      CALL  LIROBJ(TYPE,IUVF,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
C
C**** Reading the MCHAML 'FACEL' containing liquid velocity
C   
      TYPE='MCHAML  '
      CALL  LIROBJ(TYPE,IULF,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
C
C**** Reading the MCHAML 'FACEL' containing pressure
C    
      TYPE='MCHAML  '
      CALL  LIROBJ(TYPE,IPF,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
C
C**** Reading the MCHAML 'FACEL' containing vapour temperature
C
      TYPE='MCHAML  '
      CALL  LIROBJ(TYPE,ITVF,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
C
C**** Reading the MCHAML 'FACEL' containing liquid temperature
C
      TYPE='MCHAML  '
      CALL  LIROBJ(TYPE,ITLF,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
C
C**** Reading the MCHAML 'FACEL' containing vapour density
C
      TYPE='MCHAML  '
      CALL  LIROBJ(TYPE,IRVF,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
C
C**** Reading the MCHAML 'FACEL' containing liquid density
C
      TYPE='MCHAML  '
      CALL  LIROBJ(TYPE,IRLF,1,IRET)
      IF(IERR .NE. 0) GOTO 9999

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(IDIM .EQ. 3)THEN
C*** In 3D we have 10 equations
         NUMEQ = IDIM + 7
      ELSE
C*** In 2D we have 8 equations
         NUMEQ = IDIM + 6
      END IF
      IF(NC .NE. NUMEQ)THEN
         MOTERR(1:40) = 'LISTINCO = ???'
         write(*,*) NC 
         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**** Implicite ?
C
      CALL LIRMOT(LJACO,NBJAC,IIMPL,0)
      IF(IERR .NE. 0)GOTO 9999

C For the moment we only have explicit
C
C     
      IF(IIMPL .EQ. 0)THEN
C
C******** Objet matrik vide
C     
         NRIGE=7
         NMATRI=0
         NKID =9
         NKMT =7
         SEGINI MATRIK
         SEGDES MATRIK
         IJACO = MATRIK
      ELSE
C        Tentative d'utilisation d'une option non implémentée
         CALL ERREUR(251)
         GOTO 9999
      ENDIF
C
C**** Creation des flux aux interfaces
C
      TYPE = 'CHPOINT '
      CALL KRCHP1(TYPE, MELEMF, ICHFLU, MLMOTS)
C
C     
C**** Calcul des flux et du pas du temps.
C
      IF(IDIM .EQ. 2)THEN
         CALL KONF21(INDMET,
     &        IALPF,IUVF,IULF,IPF,
     &        ITVF,ITLF,IRVF,IRLF,
     &        ICHPSU,ICHPDI,
     &        MELEMC,MELEMF,MELEFE,
     &        ICHFLU,
     &        DT,
     &        LOGNC,LOGAN,MESERR)
      ELSEIF(IDIM .EQ. 3)THEN
         CALL KONF22(INDMET, 
     &        IALPF,IUVF,IULF,IPF,
     &        ITVF,ITLF,IRVF,IRLF,
     &        ICHPSU,ICHPDI,
     &        MELEMC,MELEMF,MELEFE,
     &        ICHFLU,DT,
     &        LOGNC,LOGAN,MESERR)
      ELSE
C        Erreur anormale.contactez votre support
C        Dimension has to be 2 or 3
         CALL ERREUR(5)
         GOTO 9999
      ENDIF
C
      IF(LOGAN)THEN
C
C******* Anomalie detectée
C
C
C******* Message d'erreur standard
C        -301 0
C        %m1:40
C
         MOTERR(1:40) =  MESERR(1:40)
         WRITE(IOIMP,*) MOTERR(1:40)
C
C******* Message d'erreur standard
C        5 3
C        Erreur anormale.contactez votre support
C         
         CALL ERREUR(5)
         GOTO 9999
      ENDIF
      IF(LOGNC)THEN
C
C******* Message d'erreur standard
C        -301 0
C        %m1:40
C
         MOTERR(1:40) =  MESERR(1:40)
         WRITE(IOIMP,*) MOTERR(1:40)
C
C******* Message d'erreur standard
C        460 2
C        Pas de convergence dans les itérations internes
C
         CALL ERREUR(460)
         GOTO 9999
      ENDIF
C
C**** Calcul de residu (si LOGRES = .TRUE.)         
C
C      LOGRES = .FALSE.

      IF(LOGRES)THEN
         TYPE = 'CHPOINT '
         CALL KRCHP1(TYPE, MELEMC, ICHRES, MLMOTS)
C
         CALL KONRE8(MELEMC,MELEMF,MELEFE,ICHPVO,
     &        ICHFLU, ICHRES,
     &        LOGAN,MESERR)
         IF(LOGAN)THEN
C
C******* Anomalie detectée
C
C
C******* Message d'erreur standard
C        -301 0
C        %m1:40
C
            MOTERR(1:40) =  MESERR(1:40)
            WRITE(IOIMP,*) MOTERR(1:40)
C
C******* Message d'erreur standard
C        5 3
C        Erreur anormale.contactez votre support
C         
            CALL ERREUR(5)
            GOTO 9999
         ENDIF         
      ELSE
         ICHRES = 0
      ENDIF
C
C**** Ecriture des resultats
C
      CALL ECRREE(DT)
      TYPE = 'CHPOINT '
      IF(ICHRES .NE. 0) CALL ECROBJ(TYPE,ICHRES)
      IF(ICHFLU .NE. 0) CALL ECROBJ(TYPE,ICHFLU)
      TYPE='MATRIK '
C      CALL ECROBJ(TYPE,IJACO)
 9999 CONTINUE
      RETURN
      END