Sources Esope | Cast3M

pr12f
C PR12F     SOURCE    CB215821  25/04/23    21:15:33     12247          
      SUBROUTINE PR12f(IM1,  ICH1, ICH2, ICH3,
     &                 ICH4, ICH5, ICH6, ICH7,
     &                 ICH8, ICH9,
     &                 Cp, Cvm,
     &                 OCH1, OCH2, OCH3,
     &                 OCH4, OCH5, OCH6,
     &                 OCH7, OCH8, OCH9,
     &                 LOGNEG, LOGBOR, MESERR,
     &                 VALER, VAL1, VAL2)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  PR12f
C
C DESCRIPTION       :  Primitive variables from conserved variables
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec extensions CISI)
C
C AUTEUR            :  Jose R. Garcia-Cascales,
C                      Universidad Politecnica de Cartagena,
C                      jr.garcia@upct.es
C
C************************************************************************
C
C APPELES (E/S)     :
C
C************************************************************************
C
C ENTREES  :
C
C            IM1     : MELEME contenant les centres des ELTs
C
C            ICH1    : CHPOINT  contenant la masse volumique (mass, gas/vapour)
C
C            ICH2    : CHPOINT  contenant la masse volumique (mass, liquid)
C
C            ICH3    : CHPOINT contenant les dèbits (momentum, gas/vapour)
C                     ( NDIM  composantes);
C
C            ICH4    : CHPOINT contenant les dèbits (momentum, liquid)
C                     ( NDIM  composantes);
C
C            ICH5    : CHPOINT  contenat l'énergie totale per
C                      unité de volume (alpha rho et) (gas/vapour);
C
C            ICH6    : CHPOINT  contenat l'énergie totale per
C                      unité de volume (alpha rho et) (liquid);
C
C            ICH7    : CHPOINT  containing void fraction of the previous
C                      time step, in order to avoid problems when a
C                      phase disappears
C
C            ICH8    : CHPOINT  containing vapour temperature of the previous
C                      time step, in order to avoid problems when a
C                      phase disappears
C
C            ICH9    : CHPOINT  containing liquid temperature of the previous
C                      time step, in order to avoid problems when the liquid
C                      phase disappears
C
C            Cp      : Parameter characterizing  the pressure
C                      correction term
C
C            Cvm     : Parameter characterizing virtual mass  correction,
C                      it is included in the new CATHARE pressure
C                      conrrection term
C
C SORTIES  : OCH1    : CHPOINT  void fraction (alpha);
C
C            OCH2    : CHPOINT  gas/vapour velocity (uvx, uvy, uvz);
C
C            OCH3    : CHPOINT  liquid velocity (ulx, uly, ulz);
C
C            OCH4    : CHPOINT  pressure;
C
C            OCH5    : CHPOINT  vapour temperature (Tv);
C
C            OCH6    : CHPOINT  liquid temperature (Tl);
C
C            OCH7    : CHPOINT  vapour density (rv);
C
C            OCH8    : CHPOINT  liquid density (rl);
C
C            OCH9    : CHPOINT  interfacial pressure correction term (pi)
C
C            LOGNEG  : (LOGICAL): si .TRUE. une pression ou une densité
C                      negative a été detectée -> le programme s'arrete
C                      (sa valeur stockée en MESERR(1) et VALER(1))
C
C            LOGBOR  : (LOGICAL)
C                      gamma a été detecté dehor GAMMIN et GAMMAX
C                      (sa valeur stockée en MESERR(2) et  VALER(2),VAL1,VAL2)
C
C            MESERR,
C            VALER,
C            VAL1,
C            VAL2    : pour message d'erreur
C
C
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE :  Créée le 20.11.05.
C
C************************************************************************
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
 
C**** Les variables
C
      INTEGER NLCE,N1,IGEOMC
C
      INTEGER IM1
     &       ,ICH1, ICH2, ICH3
     &       ,ICH4, ICH5, ICH6
     &       ,ICH7, ICH8, ICH9
     &       ,OCH1,OCH2,OCH3
     &       ,OCH4,OCH5,OCH6
     &       ,OCH7,OCH8,OCH9
C
      REAL*8  VALER(2),VAL1,VAL2
     &       ,alpha, uvx, uvy, uvz, uv2, ulx, uly, ulz, ul2
     &       ,p, Tv, Tl, rv, rl, ev, el
     &       ,oldav, oldTv
     &       ,A, B, delta
     &       ,gamv, cpv, gaml, cpl, pil
     &       ,FFV, FFL
     &       ,rm, rvm, rlm, Cp, cvm, pi
      PARAMETER(gamv = 1.4D0, cpv = 1008.D0,
     &          gaml = 2.8D0, cpl = 4186.D0,  pil = 8.5D8)
C
      CHARACTER*(8) TYPE
      CHARACTER*(40) MESERR(2)
      LOGICAL LOGNEG, LOGBOR
C
      INTEGER  NAT, NSOUPO, N, NC
C
C**** Les includes
C
 
-INC PPARAM
-INC CCOPTIO
-INC SMCHPOI
-INC SMCOORD
      POINTEUR MCHPO5.MCHPOI, MCHPO6.MCHPOI,
     &         MCHPO7.MCHPOI, MCHPO8.MCHPOI, MCHPO9.MCHPOI
C
      POINTEUR MSOUP6.MSOUPO, MSOUP7.MSOUPO, MSOUP8.MSOUPO,
     &         MSOUP9.MSOUPO
C
      POINTEUR MPOW1.MPOVAL, MPOW2.MPOVAL, MPOW3.MPOVAL,
     &         MPOW4.MPOVAL, MPOW5.MPOVAL, MPOW6.MPOVAL,
     &         MPOW7.MPOVAL, MPOW8.MPOVAL, MPOW9.MPOVAL
C
      POINTEUR MPOV1.MPOVAL, MPOV2.MPOVAL, MPOV3.MPOVAL,
     &         MPOV4.MPOVAL, MPOV5.MPOVAL, MPOV6.MPOVAL,
     &         MPOV7.MPOVAL, MPOV8.MPOVAL, MPOV9.MPOVAL
-INC SMELEME
C
C
C**** Initialisation des variables pour la gestion des erreurs pas ici,
C     mais avant, i.e.
C
      LOGNEG  = .FALSE.
      LOGBOR  = .FALSE.
      MESERR(1)  = '                                        '
      MESERR(2)  = '                                        '
      VALER(1) = 0.0D0
      VALER(2) = 0.0D0
      VAL1 = 0.0D0
      VAL2 = 0.0D0
 
C
C**** Activation du MELEME "CENTRE"
C
      IPT1 = IM1
      SEGACT IPT1
      N1 = IPT1.NUM(/2)
      SEGDES IPT1
C
C**** Reading MPOVA of each MCHAMPOIN
C     MPOW1 = alpha*rv
C     MPOW2 = (1 - alpha)*rl
C     MPOW3 = alpha*rv*uvx &  alpha*rv*uvy & alpha*rv*uvz
C     MPOW4 = (1 - alpha)*rl*ulx &  (1 - alpha)*rl*uly & (1 - alpha)*rl*ulz
C     MPOW5 = alpha*rv*Ev
C     MPOW6 = (1 - alpha)*rl*El
C
C     MPOW7: It is the old void fraction, we will use it when
C     alpha approaches to zero
C     MPOW8: It is the old vapour temperature, we will use it when
C     alpha approaches to zero
C     MPOW9: It is the old liquid temperature, we will use it when
C     alpha approaches to one
C
      CALL LICHT(ICH1,MPOW1,TYPE,IGEOMC)
      CALL LICHT(ICH2,MPOW2,TYPE,IGEOMC)
      CALL LICHT(ICH3,MPOW3,TYPE,IGEOMC)
      CALL LICHT(ICH4,MPOW4,TYPE,IGEOMC)
      CALL LICHT(ICH5,MPOW5,TYPE,IGEOMC)
      CALL LICHT(ICH6,MPOW6,TYPE,IGEOMC)
      CALL LICHT(ICH7,MPOW7,TYPE,IGEOMC)
      CALL LICHT(ICH8,MPOW8,TYPE,IGEOMC)
      CALL LICHT(ICH9,MPOW9,TYPE,IGEOMC)
C
C      SEGACT MPOW1
C      SEGACT MPOW2
C      SEGACT MPOW3
C      SEGACT MPOW4
C      SEGACT MPOW5
C      SEGACT MPOW6
C      SEGACT MPOW7
C      SEGACT MPOW8
C      SEGACT MPOW9
CC
C**** LICHT active les MPOVALs en *MO
C
C
C**** Privitive variables
C
C     void fraction: alpha
 
      NAT = 1
      NSOUPO=1
      N = N1
      NC = 1
      SEGINI MCHPO1
      MCHPO1.JATTRI(1)=2
      MCHPO1.IFOPOI=IFOUR
      SEGINI MSOUP1
      MCHPO1.IPCHP(1)=MSOUP1
      SEGINI MPOV1
      MSOUP1.NOCOMP(1)='SCAL'
      MSOUP1.IGEOC = IM1
      MSOUP1.IPOVAL = MPOV1
 
C     vapour velocities: uvx, uvy, uvz
 
 
      NC = IDIM
      SEGINI MCHPO2
      MCHPO2.JATTRI(1)=2
      MCHPO2.IFOPOI=IFOUR
      SEGINI MSOUP2
      MCHPO2.IPCHP(1)=MSOUP2
      SEGINI MPOV2
      MSOUP2.NOCOMP(1)='UVX'
      MSOUP2.NOCOMP(2)='UVY'
      IF(IDIM .EQ. 3)THEN
         MSOUP2.NOCOMP(3)='UVZ'
      ENDIF
      MSOUP2.IGEOC = IM1
      MSOUP2.IPOVAL = MPOV2
 
C     liquid velocities: ulx, uly, ulz
 
 
      NC = IDIM
      SEGINI MCHPO3
      MCHPO3.JATTRI(1)=2
      MCHPO3.IFOPOI=IFOUR
      SEGINI MSOUP3
      MCHPO3.IPCHP(1)=MSOUP3
      SEGINI MPOV3
      MSOUP3.NOCOMP(1)='ULX'
      MSOUP3.NOCOMP(2)='ULY'
      IF(IDIM .EQ. 3)THEN
         MSOUP3.NOCOMP(3)='ULZ'
      ENDIF
      MSOUP3.IGEOC = IM1
      MSOUP3.IPOVAL = MPOV3
 
C     pressure: p
 
      NC = 1
      SEGINI MCHPO4
      MCHPO4.JATTRI(1)=2
      MCHPO4.IFOPOI=IFOUR
      SEGINI MSOUP4
      MCHPO4.IPCHP(1)=MSOUP4
      SEGINI MPOV4
      MSOUP4.NOCOMP(1)='SCAL'
      MSOUP4.IGEOC = IM1
      MSOUP4.IPOVAL =  MPOV4
 
C     vapour temperature: Tv
 
      NC = 1
      SEGINI MCHPO5
      MCHPO5.JATTRI(1)=2
      MCHPO5.IFOPOI=IFOUR
      SEGINI MSOUP5
      MCHPO5.IPCHP(1)=MSOUP5
      SEGINI MPOV5
      MSOUP5.NOCOMP(1)='SCAL'
      MSOUP5.IGEOC = IM1
      MSOUP5.IPOVAL = MPOV5
 
C     liquid temperature: Tl
 
      NC = 1
      SEGINI MCHPO6
      MCHPO6.JATTRI(1)=2
      MCHPO6.IFOPOI=IFOUR
      SEGINI MSOUP6
      MCHPO6.IPCHP(1)=MSOUP6
      SEGINI MPOV6
      MSOUP6.NOCOMP(1)='SCAL'
      MSOUP6.IGEOC = IM1
      MSOUP6.IPOVAL = MPOV6
 
C     vapour density: rv
 
      NC = 1
      SEGINI MCHPO7
      MCHPO7.JATTRI(1)=2
      MCHPO7.IFOPOI=IFOUR
      SEGINI MSOUP7
      MCHPO7.IPCHP(1)=MSOUP7
      SEGINI MPOV7
      MSOUP7.NOCOMP(1)='SCAL'
      MSOUP7.IGEOC = IM1
      MSOUP7.IPOVAL = MPOV7
 
C     liquid density: rl
 
      NC = 1
      SEGINI MCHPO8
      MCHPO8.JATTRI(1)=2
      MCHPO8.IFOPOI=IFOUR
      SEGINI MSOUP8
      MCHPO8.IPCHP(1)=MSOUP8
      SEGINI MPOV8
      MSOUP8.NOCOMP(1)='SCAL'
      MSOUP8.IGEOC = IM1
      MSOUP8.IPOVAL = MPOV8
 
C     pressure correction term, pi
 
      NC = 1
      SEGINI MCHPO9
      MCHPO9.JATTRI(1)=2
      MCHPO9.IFOPOI=IFOUR
      SEGINI MSOUP9
      MCHPO9.IPCHP(1)=MSOUP9
      SEGINI MPOV9
      MSOUP9.NOCOMP(1)='SCAL'
      MSOUP9.IGEOC = IM1
      MSOUP9.IPOVAL = MPOV9
C
C**** "Recapitulando"
C
C     We activate MPOV1 - MPOV9
C
C**** "bucle" for the calculation of the primitives at
C     each centre
 
 
      DO NLCE = 1, N1
 
         oldav = MPOW7.VPOCHA(NLCE,1)
         oldTv = MPOW8.VPOCHA(NLCE,1)
 
         IF(IDIM .EQ. 3)THEN
            if(oldav .LT. 0.01d0)then
            A = (gamv - 1.D0)*cpv*oldTv*MPOW1.VPOCHA(NLCE,1)/gamv
            else
               A = (gamv - 1.D0)*(MPOW5.VPOCHA(NLCE,1) - (MPOW3
     $              .VPOCHA(NLCE,1)**2 + MPOW3.VPOCHA(NLCE,2)**2 +
     &              MPOW3.VPOCHA(NLCE,3)**2)/
     &             (2.d0*(MPOW1.VPOCHA(NLCE,1)+ 1.D-20)))
            endif
            if(oldav .GT. 0.99d0)then
                B = (gaml - 1.D0)*MPOW2.VPOCHA(NLCE,1)*
     &              (cpl*MPOW9.VPOCHA(NLCE,1)/gaml + pil/1000.d0)
            else
                B = (gaml - 1.D0)*(MPOW6.VPOCHA(NLCE,1) - (MPOW4
     $              .VPOCHA(NLCE,1)**2 + MPOW4.VPOCHA(NLCE,2)**2)+
     &              MPOW4.VPOCHA(NLCE,3)**2/(2.D0
     $              *(MPOW2.VPOCHA(NLCE,1)+ 1.D-20)))
            endif
         ELSE
            if(oldav .LT. 0.01d0)then
               A = (gamv - 1.D0)*cpv*oldTv*MPOW1.VPOCHA(NLCE,1)/gamv
            else
               A = (gamv - 1.D0)*(MPOW5.VPOCHA(NLCE,1) - (MPOW3
     $              .VPOCHA(NLCE,1)**2 + MPOW3.VPOCHA(NLCE,2)**2)/(2.d0
     $              *(MPOW1.VPOCHA(NLCE,1)+ 1.D-20)))
            endif
            if(oldav .GT. 0.99d0)then
                B = (gaml - 1.D0)*MPOW2.VPOCHA(NLCE,1)*
     &          (cpl*MPOW9.VPOCHA(NLCE,1)/gaml + pil/1000.d0)
            else
 
               B = (gaml - 1.D0)*(MPOW6.VPOCHA(NLCE,1) - (MPOW4
     $              .VPOCHA(NLCE,1)**2 + MPOW4.VPOCHA(NLCE,2)**2)/(2.D0
     $              *(MPOW2.VPOCHA(NLCE,1)+ 1.D-20)))
            endif
         ENDIF
         delta = (- A - B + gaml*pil)**2.D0 + 4.D0*gaml*pil*A
 
         p = 1.D0/2.D0*(A + B - gaml*pil + sqrt(delta))
 
         alpha = A/p
         call FUNCF(alpha, FFV)
         uvx = FFV*MPOW3.VPOCHA(NLCE,1)/(MPOW1.VPOCHA(NLCE,1) + 1.D-20)
         uvy = FFV*MPOW3.VPOCHA(NLCE,2)/(MPOW1.VPOCHA(NLCE,1) + 1.D-20)
         IF(IDIM .EQ. 3)THEN
            uvz = FFV*MPOW3.VPOCHA(NLCE,3)/(MPOW1.VPOCHA(NLCE,1) + 1.D
     &           -20)
            uv2 = uvx*uvx + uvy*uvy + uvz*uvz
         ELSE
            uv2 = uvx*uvx + uvy*uvy
         ENDIF
 
         call FUNCF((1.d0 - alpha),FFL)
         ulx = FFL*MPOW4.VPOCHA(NLCE,1)/(MPOW2.VPOCHA(NLCE,1) + 1.D-20)
         uly = FFL*MPOW4.VPOCHA(NLCE,2)/(MPOW2.VPOCHA(NLCE,1) + 1.D-20)
         IF(IDIM .EQ. 3)THEN
            ulz = FFL*MPOW4.VPOCHA(NLCE,3)/(MPOW2.VPOCHA(NLCE,1) + 1.D
     &           -20)
            ul2 = ulx*ulx + uly*uly + ulz*ulz
         ELSE
            ul2 = ulx*ulx + uly*uly
         ENDIF
 
         if((alpha .GT. 0.999d0) .OR. (alpha .LT. 0.01d0))then
C
C Water faucet test
C              Tv = 323.15D0
C              Tl = 323.15D0
C
C Shock tube test
C              Tv = 308.15d0
C              Tl = 308.15d0
C
C Phase separation test or slosing of a water column
              Tv = 298.15D0
              Tl = 298.15D0
C For the oscillating manometer should be 323 but to be practical I have left it in 298.15
              ev = cpv*Tv/gamv
              el = cpl*Tl/gaml*(p + gaml*pil)/(p + pil)
         else
              ev = MPOW5.VPOCHA(NLCE,1)/(MPOW1.VPOCHA(NLCE,1) + 1.D-20)
     &             - uv2/2.D0
              el = MPOW6.VPOCHA(NLCE,1)/(MPOW2.VPOCHA(NLCE,1) + 1.D-20)
     &             - ul2/2.D0
              Tv = gamv*ev/cpv
              Tl = gaml*el/(cpl*(1.D0 + pil*(gaml - 1.D0)/(p + pil)))
         endif
         rv = gamv*p/((gamv - 1.D0)*cpv*Tv)
         rl = gaml*(p + pil)/((gaml - 1.D0)*cpl*Tl)
C
C     Bestion´s formula
C
         rm  = alpha*rv + (1.d0 - alpha)*rl
         rvm = rv + Cvm*rm
         rlm = rl + Cvm*rm
 
         pi = p - Cp*(alpha*rvm*(1.D0- alpha)*rlm)/
     &               (alpha*rlm + (1.D0- alpha)*rvm)*
     &               (sqrt(uv2) - sqrt(ul2))**2
 
         MPOV1.VPOCHA(NLCE,1) = alpha
         MPOV2.VPOCHA(NLCE,1) = uvx
         MPOV2.VPOCHA(NLCE,2) = uvy
         MPOV3.VPOCHA(NLCE,1) = ulx
         MPOV3.VPOCHA(NLCE,2) = uly
         IF(IDIM .EQ. 3)THEN
            MPOV2.VPOCHA(NLCE,3) = uvz
            MPOV3.VPOCHA(NLCE,3) = ulz
         ENDIF
         MPOV4.VPOCHA(NLCE,1) = p
         MPOV5.VPOCHA(NLCE,1) = Tv
         MPOV6.VPOCHA(NLCE,1) = Tl
         MPOV7.VPOCHA(NLCE,1) = rv
         MPOV8.VPOCHA(NLCE,1) = rl
         MPOV9.VPOCHA(NLCE,1) = pi
 
      ENDDO
C
      OCH1 = MCHPO1
      OCH2 = MCHPO2
      OCH3 = MCHPO3
      OCH4 = MCHPO4
      OCH5 = MCHPO5
      OCH6 = MCHPO6
      OCH7 = MCHPO7
      OCH8 = MCHPO8
      OCH9 = MCHPO9
C
      SEGDES MPOV1
      SEGDES MPOV2
      SEGDES MPOV3
      SEGDES MPOV4
      SEGDES MPOV5
      SEGDES MPOV6
      SEGDES MPOV7
      SEGDES MPOV8
      SEGDES MPOV9
 
      SEGDES MPOW1
      SEGDES MPOW2
      SEGDES MPOW3
      SEGDES MPOW4
      SEGDES MPOW5
      SEGDES MPOW6
      SEGDES MPOW7
      SEGDES MPOW8
      SEGDES MPOW9
 
      SEGDES MSOUP1
      SEGDES MSOUP2
      SEGDES MSOUP3
      SEGDES MSOUP4
      SEGDES MSOUP5
      SEGDES MSOUP6
      SEGDES MSOUP7
      SEGDES MSOUP8
      SEGDES MSOUP9
 
      SEGDES MCHPO1
      SEGDES MCHPO2
      SEGDES MCHPO3
      SEGDES MCHPO4
      SEGDES MCHPO5
      SEGDES MCHPO6
      SEGDES MCHPO7
      SEGDES MCHPO8
      SEGDES MCHPO9
 
      RETURN
      END