C BSUR      SOURCE    CB215821  16/04/21    21:15:24     8920
      SUBROUTINE BSUR (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)

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION XN(NPP),TN(NPP),EN(NPP),BN(NPP)

C     operateur FUITE
C>>>  superheated vapour
C     version with wall condensation
C
C     QAE:  QA.E  constant air flowrate
C     QEEi: QEi.E water flowrate
C     indice 1,2 : entree et sortie de troncon
C     RECU,XKUL,XKUT1,XKUT2,XKUT3,XKUT4 : coef lois de frot utilisateur
C
C  H correspond au flux d'energie total cede a la paroi
C  HM (convection) est utilise pour le calcul de T2
C
      HW=0.D0
      RINDEX = 1.0
      DQEE = 0.D0

      CALL BPHYS(T0,P0,RA,RV,CA,CV,CL,XLAT0,ROL,XKL,XKT,XREL)

C     interpolation de la temperature de paroi
      CALL BTPAR(XN,TN,X,TP,NPP,ITP)
      TP=TP+T0
C     epaisseur a l'entree
      CALL BTPAR(XN,EN,X,E1,NPP,ITP)
C     etendue a l'entree
      CALL BTPAR(XN,BN,X,B1,NPP,ITP)

      QA1=QAE/E1/B1
      QE1=QEE1/E1/B1
      Q = QA1 + QE1
C     calcul de R et CP
      CALL BCAR(RA,RV,PHI1,CV,CA,Q,RQ,R,QA0,QV,CP)
      IF(KIMP.GE.1) THEN
        write(6,*) 'entree bsur : E1,B1,PHI1= ',E1,B1,PHI1
      ENDIF

      RO=P1/R/T1
      U1=Q/RO

      U=U1
      T=T1
      P=P1
      QE=QE1
      QA=QA1
      QX=QA+QE
      PHI=PHI1

      IF ((X+DX).GT.1.) DX=1-X
      X=X+DX

      DU=0
      T2=0
      P2=0
c     ouverture en sortie de maille
      CALL BTPAR(XN,EN,X,E2,NPP,ITP)
c     etendue en sortie de maille
      CALL BTPAR(XN,BN,X,B2,NPP,ITP)
      QA2=QAE/E2/B2
      E=E1
      B=B1

      NITER=0

  10  CONTINUE

      NITER=NITER+1

c      write(6,*) ' bsur X,NITER,T1,P1,RE',X,NITER,T1,P1,RE
      TT2=T2
      PP2=P2

      XMU=BMUG(T,PHI,T0)
      RE=QX*2*E/XMU
      PR=XMU*CP/BLA(T,T0)
      DPP = XL*QX*QX/4/E
      BK = BKFRO(RE,REL,XKL,XKT,2*E,RUG,RECU,XKUL,XKUT1,XKUT2,XKUT3,
     *XKUT4)
      DP=DX*BK*DPP/RO

      P2=P1-DP
         IF (P2.LT.PSLIM) THEN
            PSQ = -1.
            RETURN
         ENDIF

      S=0.D0
C*
      XLA=BLA(T,T0)

      HM = BHECH(RE,REL,PR,XLA,E)

C if vapour at the inlet

      IF(PHI1.GT.1.D-2) THEN
C
C wall condensation
         D=2.5D-5
         A=2.3D-5
         PP=BPSAT(TP)
         PV=PHI1*P1
         RAPP=(P1-PP)/(P1-PV)

         IF (RAPP.GT.(1.0001)) THEN
         XJ=XW*(HM/CV)*(D/A)**(2./3.)*LOG(RAPP)
         ELSE
         XJ=0.D0
         ENDIF

         DQEE=-2.*XJ*DX*XL*B
         XLAT=XLAT0+(CV-CL)*T
         SW=(XLAT/QX/CP/E)*DQEE
C**
         DELTAT = T-TP
         IF (ABS(DELTAT).GT.1D-3) THEN
         HW=XJ*XLAT/(T-TP)
C        write(6,*) ' T TP HW ',T,TP,HW
         ENDIF

         QEE2=QEE1+DQEE
         QE2=QEE2/E2/B2

         IF ((QEE2/QEE1).LT.(-1.D-5)) THEN
           write(6,*) 'bsur : relative water flowrate ',QE2/QE1
           RINDEX = 0.5
           RETURN
         ELSE
           A=RA*QA2/RV/QE2
           PHI2=1./(1.+A)
         ENDIF
* no vapour at the inlet
      ELSE
         XLAT=XLAT0+(CV-CL)*T
         HW=0.D0
         PHI2=0.D0
         QEE2=0.D0
         QE2=0.D0
      ENDIF
      IF(KIMP.GE.1) THEN
         write(6,*) 'bsur: NITER,PHI2,QA2,QE2= ',NITER,PHI2,QA2,QE2
      ENDIF
      QEEX=(QEE1+QEE2)/2.

      T2 = BT2(T1,TP,S,HM,PR,RE,XLA,XL,DX)
      IF(KIMP.GE.1) WRITE(*,*) 'bsur : T2= ',T2

      U2 = (QA2+QE2)*R*T2/P2
      RO=(P1+P2)/R/(T1+T2)
      T=(T1+T2)/2
      P=(P1+P2)/2
      PHI=(PHI1+PHI2)/2

      E=(E1+E2)/2
      B=(B1+B2)/2
      QA=QAE/E/B
      QE=QEEX/E/B
      QX=QA+QE

      QW2=(QEE2-QEE1)/DX/XL/B

      ERT=ABS((T2-TT2)/T2)
      ERP=ABS((P2-PP2)/P2)

      IF(((ERT.GT.1E-4).OR.(ERP.GT.1E-4))
     $  .AND.(NITER.LE.10))    GOTO 10

         IF(KIMP.GE.2) THEN
         write(6,*) ' HM HW HW/HM ',HM,HW,HW/HM
         ENDIF

C     calcul des puissances echangees

      RAPH=(CV*T+XLAT0)/XLAT
      H=HM+HW*RAPH

C     WRITE(6,*) ' bsur ',H

      HH1=(QAE*CA*T1)+(QEE1*(CV*T1+XLAT0))
      HH2=(QAE*CA*T2)+(QEE2*(CV*T2+XLAT0))
C !   PF = HH2-HH1
      PF = (HH2-HH1) - DQEE*(CV*T1+XLAT0)
      PP = H*(T-TP)*2*DX*XL*B
      IF(KIMP.EQ.-3) THEN
        write(6,*) ' DQEE CV T1 XLAT0 ',DQEE,CV,T1,XLAT0
        write(6,1000)  X,RAPH,PF,PP,HH1,HH2
1000    FORMAT(1X,' bsur X RH DH EC HH1 HH2 ',6E12.5)
      ENDIF
      IF(KIMP.NE.0.AND.X.EQ.1.) THEN
        write(6,2110) X,RE,BK
2110    format(1X,'bsur X RE BKRO ',3E12.5)
      ENDIF
      DPF=PF/DX/XL/B
      DPP=PP/DX/XL/B

      RETURN
      END