C BCALQ2    SOURCE    CB215821  16/04/21    21:15:15     8920
      SUBROUTINE BCALQ2(PE,PVE,TE,PSO,XL,DX0,RUG,Q,XW,NPP,XN,TN,EN,BN,
     $             KIMP,NT,NX,XX,XP,XT,XY,XU,XHF,XQ,XQW,
     $             QAE,XRE,XDH,PSQ,RECU,XKUL,XKUT1,XKUT2,XKUT3,XKUT4)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)

C     operateur FUITE
C     cf. sub. BECALC2
C     calcul de la solution pour un debit Q d'injection impose
C     QAE, XQ : debit air, debit eau,
C     XQW: flux d'eau a la paroi
C     RECU,XKUL,XKUT1,XKUT2,XKUT3,XKUT4 : coef lois de frot utilisateur
C
      DIMENSION XX(NT),XP(NT),XT(NT),XY(NT),XU(NT),XN(NPP),TN(NPP)
      DIMENSION EN(NPP),BN(NPP)
      DIMENSION XHF(NT),XQ(NT),XQW(NT),XRE(NT),XDH(NT)

      IF(KIMP.GE.1) THEN
        WRITE(6,*)
        WRITE(6,*) 'entree bcalq2'
      ENDIF
C     variation enthalpie et puissance echangée avec la paroi
      DHTOT = 0.D0
      PPTOT = 0.D0

      DTMX=10
      PSLIM=PSO
      NITMAX = 100

      PHI = PVE/PE

C     IMPLEMENTATION DES PROPRIETES PHYSIQUES
      CALL BPHYS(T0,P0,RA,RV,CA,XCV,XCL,XLAT0,XROL,XKL,XKT,REL)

C     Calcul des debits en entree Qa et Qv

C     Calcul de la masse volumique de vapeur entree

       RHOVE=ROVAP0(PVE,TE)

C     Calcul du facteur Z Compressibilite en entree

       ZE=ZVAP0(RHOVE,TE)


C     Calcul des debits air et eau, le debit total Q etant connu

       IF(PHI.LT.1D-3) THEN
         RQ = 1.E8
       ELSE
         RQ = (1-PHI)*(RV*ZE)/PHI/RA
       ENDIF
C
       QA = Q*RQ/(1+RQ)
       QV = Q/(1+RQ)

C !!
C     QAE invariant

      QAE=QA*EN(1)*BN(1)

C      'constante de gaz' equivalente a l'entree

       R = PHI/RV/ZE+(1-PHI)/RA
       R=1/R

      ROE = PE/R/TE
      UE=Q/ROE

      IF(KIMP.NE.0) THEN
      write(6,*)
      write(6,*) ' bcalq2: conditions en entree'
      write(6,2000) TE,RHOVE,ZE,R,QA,QV
 2000 FORMAT(1X,'   TE ROV ZE R QA QV', 6E12.5)
      ENDIF

C>>>  positionnement a l'entree de la fissure

      PSQ=1
      IX=1

      P2=PE
      PV2=PVE
      T2=TE
      U2=UE

      QEE2=QV*EN(1)*BN(1)
      PHI2=PHI

      IF(KIMP.NE.0) THEN
      write(6,*) '>BCALQ2 Q PHI R ZE ',Q,PHI,R,ZE
      ENDIF

      X=0.D0

      TIT=1
      ITP=1

C>>>  le gaz a l'entree est-il surchauffe ou sature

      PS2=PSATT0(T2)
      R2=ROVAP0(PS2,T2)
      ZV =ZVAP0(R2,T2)
      IF (PVE.GT.PS2) THEN
         TIT=0
         Y2=(PS2/(P2-PS2))*(1-PHI)/(PHI*ZV)
      ELSE
         Y2=1
      ENDIF

      QTOT=QEE2+QAE

      CALL BSTOK(IX,X,P2,T2-T0,Y2,U2,0.D0,QEE2,0.D0,0.D0,0.D0,
     &               XX,XP,XT,XY,XU,XHF,XQ,XQW,XRE,XDH,NT)

C>>>  boucle "tant que" sur la longueur de la fissure


  10  CONTINUE
      IF ((X.LT.0.999).AND.(PSQ.NE.-1.)) THEN
      IF(KIMP.GE.1) THEN
        write(6,*)
        write(6,*) 'bcalq2 X',X
      ENDIF

C>>>  positionnement au point local 1

         P1=P2
         PS1=PSATT0(T2)
         T1=T2
         Y1=Y2

         QEE1=QEE2
         PHI1=PHI2

         DX=DX0

C>>>  test sur le titre

         IF (TIT.EQ.1) THEN

             Y2=1
             NIT = 0
  20         CONTINUE
             NIT = NIT + 1
             IF (NIT.GT.99) THEN
                write(6,*) 'bcalq2: NIT20=100 X',X
             ENDIF

           CALL BSUR2 (X,DX,XL,RUG,XW,XN,TN,EN,BN,KIMP,PSLIM,REL,RINDEX,
     &     P1,T1,QAE,QEE1,PHI1,P2,T2,U2,QEE2,PHI2,QW2,RE,H,PSQ,
     &     NPP,ITP,PF,PP,DPF,DPP,RECU,XKUL,XKUT1,XKUT2,XKUT3,XKUT4)

             IF (PSQ.EQ.-1.) THEN
                RETURN
             ENDIF

             IF ((RINDEX.LT.0.999).AND.
     &       (NIT.LE.NITMAX)) THEN
                X=X-DX
                DX=DX*0.5
                GO TO 20
             ENDIF


C>>>  y-a-t'il condensation ?

             PS2=PSATT0(T2)
C            WRITE(6,*) 'PS2',PS2
             PSI=PS2/P2
C            write(6,*) ' PSI ',PSI
             AL=PHI2/PSI

             IF(KIMP.GE.2) THEN
             WRITE(6,998) X,AL,P2/P0,T2-T0,QEE1/EN(1),QEE2/EN(1),U2
 998         FORMAT(1X,'sur -XALP2T2 ',F9.6,F9.4,E12.5,4F9.2)
             WRITE(6,*)
             ENDIF

             IF ((AL.GT.(1.02)).OR.(ABS(T2-T1).GT.DTMX).AND.
     &       (NIT.LE.NITMAX)) THEN
                 X=X-DX
                 DX=DX/2
                 GOTO 20
             ENDIF

             IF (AL.GT.1.) THEN
               TIT=0
C
               PS2=PSATT0(T2)
               RHO2=ROVAP0(PS2,T2)
               ZVAP2=ZVAP0(RHO2,T2)
               Y2 = (RA/RV)*(QAE/QEE2)*(PSI/(1-PSI))
               Y2 = Y2/ZVAP2

              IF(KIMP.NE.0) THEN
                write(6,*) 'bcalq2: QAE QEE2 ',QAE,QEE2
                write(6,*) 'bcalq2: transition vers condensation X= ',X
                write(6,2100) PS2,AL,RHO2,ZVAP2,Y2
 2100           FORMAT(1X,'bcalq2: PS2 AL RO2 Z2 Y2 ', 6E12.5)
              ENDIF

C
             ENDIF


         ELSE

             NIT = 0
  30         CONTINUE
             NIT = NIT + 1
             IF (NIT.GT.99) THEN
                write(6,*) 'bcalq2: NIT30=100 X',X
             ENDIF

             CALL BCOND2 (X,DX,XL,RUG,XW,XN,TN,EN,BN,KIMP,PSLIM,REL,
     &       P1,PS1,T1,Y1,QAE,QEE1,PHI1,
     &       P2,T2,Y2,U2,QEE2,PHI2,QW2,RE,H,PSQ,RINDEX,
     &       NPP,ITP,PF,PP,DPF,DPP,RECU,XKUL,XKUT1,XKUT2,XKUT3,XKUT4)
C
C            WRITE(6,*)'****** RINDEX,X= ',RINDEX,X

             IF(KIMP.GE.2) THEN
             WRITE(6,999) X,Y2,P2/P0,T2-T0,QEE1/EN(1)/BN(1),
     &       QEE2/EN(1)/BN(1),U2
 999         FORMAT(1X,'bcalq2: XY2P2T2 ',F9.6,F9.4,E12.5,4F9.2)
             ENDIF

             IF (PSQ.EQ.-1.) THEN
                RETURN
             ENDIF

             IF ((RINDEX.LT.0.999).AND.
     &       (NIT.LE.NITMAX)) THEN
                X=X-DX
                DX=DX*0.5
                GO TO 30
             ENDIF

C>>>        y-a-t'il evaporation totale ?

             IF ((Y2.GT.(1.01)).OR.(ABS(T2-T1).GT.DTMX).AND.
     &       (NIT.LE.NITMAX)) THEN
                X=X-DX
                DX=DX/2
                GOTO 30
             ENDIF

* vapeur surchauffee
             IF (Y2.GT.0.999) THEN
                TIT=1
* attention si E variable : OK
               IF((QEE2/QAE).LT.1.D-7) THEN
C                write(6,*) ' QEE2/QAE ',QEE2/QAE
                 PHI2=0.D0
               ELSE


                 ZTRAN = 1.D0
                 NZ=0
 31              CONTINUE
                 NZ = NZ + 1

                 AA = (QEE2/QAE)* (RV*ZTRAN/RA)
                 PHI2 = AA/(1.+AA)
                 PV2 = PHI2 * P2

                 RHO2 = ROVAP0(PV2,T2)
                 Z2  = ZVAP0(RHO2,T2)

                 DELTAZ = ABS(Z2 - ZTRAN)

                 IF(KIMP.GE.2) THEN
                 WRITE(6,*) ' NZ ZVAP ', NZ,DELTAZ
                 ENDIF
                 ZTRAN=Z2

                 IF ((DELTAZ.GT.1D-5).AND.(NZ.LE.20)) GOTO 31


              IF(KIMP.NE.0) THEN
                 write(6,*) 'bcalq2: transition vers evaporation  '
                 write(6,2200) PHI2,RHO2,Z2
 2200            FORMAT(1X,'bcalq2 PHI2 RHO2 Z2 ', 3E12.5)
              ENDIF

               ENDIF
             ENDIF

         ENDIF

         QTOT=QEE2+QAE

         CALL BSTOK(IX,X,P2,T2-T0,Y2,U2,H,QEE2,QW2,RE,DPF,
     &                XX,XP,XT,XY,XU,XHF,XQ,XQW,XRE,XDH,NT)

      DHTOT = DHTOT + PF
      PPTOT = PPTOT + PP

      GOTO 10
      ENDIF

      NX=IX-1
      PSQ=P2

      XX(NX)=1.D0
C     write(6,*) ' bcalq  NX ',NX
C     CALL UTPRIM(XX,NX)
C**
      IF(KIMP.GT.0) THEN
      WRITE(*,*) 'bcalq: X,PSQ= ',X,PSQ
      ENDIF
      IF(KIMP.GT.1) THEN
      write(6,1000) Q,PHI,DHTOT,PPTOT
 1000 FORMAT(1X,'bcalq2 Q phi DH H*T',4E12.5)
      ENDIF
C**

      RETURN
      END