        SUBROUTINE flash(NumComp,z,K,MaxError,V,x,y,NumIter,F)

c***********************************************************************
c       This subroutine calculates the compositions of the liquid and
c          vapor phases resulting when a mixture is flashed utilizing
c          specified equilibrium ratios.

c       Author: Scott K. Laudeman
c       Date: June 20, 1987

c       The derivation of the flash calculation equation is described
c          in Amyx, James W.,Bass, Daniel M., Jr., and Whiting,
c          Robert L., Petroleum Reservoir Engineering, McGraw-Hill,
c          1960, pp.342-343.  The resulting equation is solved using
c          Newton's method (see any text on numerical analysis, for
c          example, Gerald, Curtis F. and Wheatly, Patrick O., Applied
c          Numerical Analysis, Addison-Wesley, 1984, pp. 15-16.

c       Input variables.

c          NumComp: Number of components.
c          z: Array containing mole fraction of each component in
c             mixture.
c          K: Array containing equilibrium ratio of each component.
c          MaxError: Error tolerance for total mole fraction of vapor
c             phase.

c       Output variables.

c          V: Total mole fraction of vapor phase.
c          x: Array containing mole fraction of each component in liquid
c             phase.
c          y: Array containing mole fraction of each component in vapor
c             phase.
c          NumIter: Iterations required for calculation.
c          F: Function to be solved for total mole fraction of vapor
c             phase.

c       Other variables.

c          dFdV: Derivative of function to be solved for total mole
c             fraction of vapor phase.
c          i: loop counter
c          MaxIter: Maximum iterations permitted before error message.
c          TermUnit: Number identifying unit for screen input/output.
c          Vold: Value of V from previous iteration.

c***********************************************************************


c***********************************************************************
c       Declarations
c***********************************************************************

c=======================================================================
c       Constants
c=======================================================================

        PARAMETER (MaxIter=20)

c=======================================================================
c       Variables
c=======================================================================


        REAL    K(NumComp),MaxError,x(NumComp),y(NumComp),z(NumComp)

c***********************************************************************
c       Main Program
c***********************************************************************

        NumIter = 0
        V       = 0.5
        Vold    = 0.0

c=======================================================================
c       Loop for calculating vapor phase mole fraction using Newton's
c          method.  Iteration continues until difference between
c          iterations is acceptable or until number of iterations exceeds
c          maximum.
c=======================================================================

10      IF (ABS(V - Vold).GT.MaxError .AND. NumIter.LE.MaxIter) THEN
           NumIter = NumIter + 1
           Vold    = V

           F       = 0.0
           dFdV    = 0.0
           DO 20 i = 1,NumComp
              F    =    F + (K(i) - 1.0)*z(i)/((K(i) - 1.0)*V + 1.0)
              dFdV = dFdV - (K(i) - 1.0)**2*z(i)/((K(i) - 1.0)*V
     1                 + 1.0)**2
20         CONTINUE

           V = Vold - F/dFdV

c----------------------------------------------------------------------
c       Reset vapor phase mole fraction if calculated value is
c          unrealistic.  If calculation fails to converge after maximum
c          iterations then single phase is assumed and the following
c          conditional is used to assign the phase.
c----------------------------------------------------------------------

          IF (V.GT.1.0) V = 1.0
          IF (V.LT.0.0) V = 0.0

          GO TO 10
        ENDIF

c=======================================================================
c       Calculate liquid and vapor compositions.
c=======================================================================

        DO 30 i = 1,NumComp
           x(i) = z(i)/((K(i) - 1.0)*V + 1.0)
           y(i) = K(i)*x(i)
30      CONTINUE

        RETURN

        END
