diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
index 8bd8a4d14006b7538c3644548379b3d866f762a2..61dfb6c6ee6e9df7d44267df68ba805e00631ff3 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
@@ -151,6 +151,44 @@ void Foam::LiquidEvaporation<CloudType>::calculate
     scalarField& dMassPC
 ) const
 {
+    // liquid volume fraction
+    const scalarField X(liquids_.X(Yl));
+
+    // immediately evaporate mass that has reached critical condition
+    if (mag(T - liquids_.Tc(X)) < SMALL)
+    {
+        if (debug)
+        {
+            WarningIn
+            (
+                "void Foam::LiquidEvaporation<CloudType>::calculate"
+                "("
+                    "const scalar, "
+                    "const label, "
+                    "const scalar, "
+                    "const scalar, "
+                    "const scalar, "
+                    "const scalar, "
+                    "const scalar, "
+                    "const scalar, "
+                    "const scalar, "
+                    "const scalar, "
+                    "const scalarField&, "
+                    "scalarField&"
+                ") const"
+            )   << "Parcel reached critical conditions: "
+                << "evaporating all avaliable mass" << endl;
+        }
+
+        forAll(activeLiquids_, i)
+        {
+            const label lid = liqToLiqMap_[i];
+            dMassPC[lid] = GREAT;
+        }
+
+        return;
+    }
+
     // construct carrier phase species volume fractions for cell, cellI
     const scalarField Xc(calcXc(cellI));
 
@@ -242,9 +280,27 @@ Foam::scalar Foam::LiquidEvaporation<CloudType>::dh
 
 
 template<class CloudType>
-Foam::scalar Foam::LiquidEvaporation<CloudType>::TMax(const scalar pIn) const
+Foam::scalar Foam::LiquidEvaporation<CloudType>::Tvap
+(
+    const scalarField& Y
+) const
+{
+    const scalarField X(liquids_.X(Y));
+
+    return liquids_.Tpt(X);
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::LiquidEvaporation<CloudType>::TMax
+(
+    const scalar p,
+    const scalarField& Y
+) const
 {
-    return liquids_.TMax(pIn);
+    const scalarField X(liquids_.X(Y));
+
+    return liquids_.pvInvert(p, X);
 }
 
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H
index fb2e2cdd9c8f38fb7f720e3512cb439ee7aff730..3d5334cd884e759a17eda9f9417d77f523274e29 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H
@@ -131,8 +131,11 @@ public:
             const scalar T
         ) const;
 
+        //- Return vapourisation temperature
+        virtual scalar Tvap(const scalarField& Y) const;
+
         //- Return maximum/limiting temperature
-        virtual scalar TMax(const scalar pIn) const;
+        virtual scalar TMax(const scalar p, const scalarField& Y) const;
 };
 
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C
index 6436d3960538501aa796ac7ce1ad147daeff48af..d78d6b110b4ee13944231af8c0ac42510414e257 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C
@@ -151,18 +151,53 @@ void Foam::LiquidEvaporationBoil<CloudType>::calculate
     scalarField& dMassPC
 ) const
 {
-    // construct carrier phase species volume fractions for cell, cellI
-    const scalarField XcMix(calcXc(cellI));
-
     // liquid volume fraction
     const scalarField X(liquids_.X(Yl));
 
+    // immediately evaporate mass that has reached critical condition
+    if (mag(T - liquids_.Tc(X)) < SMALL)
+    {
+        if (debug)
+        {
+            WarningIn
+            (
+                "void Foam::LiquidEvaporationBoil<CloudType>::calculate"
+                "("
+                    "const scalar, "
+                    "const label, "
+                    "const scalar, "
+                    "const scalar, "
+                    "const scalar, "
+                    "const scalar, "
+                    "const scalar, "
+                    "const scalar, "
+                    "const scalar, "
+                    "const scalar, "
+                    "const scalarField&, "
+                    "scalarField&"
+                ") const"
+            )   << "Parcel reached critical conditions: "
+                << "evaporating all avaliable mass" << endl;
+        }
+
+        forAll(activeLiquids_, i)
+        {
+            const label lid = liqToLiqMap_[i];
+            dMassPC[lid] = GREAT;
+        }
+
+        return;
+    }
+
     // droplet surface pressure assumed to surface vapour pressure
     scalar ps = liquids_.pv(pc, Ts, X);
 
     // vapour density at droplet surface [kg/m3]
     scalar rhos = ps*liquids_.W(X)/(specie::RR*Ts);
 
+    // construct carrier phase species volume fractions for cell, cellI
+    const scalarField XcMix(calcXc(cellI));
+
     // carrier thermo properties
     scalar Hsc = 0.0;
     scalar Hc = 0.0;
@@ -340,13 +375,28 @@ Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::dh
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::Tvap
+(
+    const scalarField& Y
+) const
+{
+    const scalarField X(liquids_.X(Y));
+
+    return liquids_.Tpt(X);
+}
+
+
 template<class CloudType>
 Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::TMax
 (
-    const scalar pIn
+    const scalar p,
+    const scalarField& Y
 ) const
 {
-    return liquids_.TMax(pIn);
+    const scalarField X(liquids_.X(Y));
+
+    return liquids_.pvInvert(p, X);
 }
 
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.H
index 2d11f55f3200869863ffc7632debce83ce2efad0..f3b532dff5827c7bd2a2f6081e774ea60287e1b7 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.H
@@ -141,8 +141,11 @@ public:
             const scalar T
         ) const;
 
+        //- Return vapourisation temperature
+        virtual scalar Tvap(const scalarField& Y) const;
+
         //- Return maximum/limiting temperature
-        virtual scalar TMax(const scalar pIn) const;
+        virtual scalar TMax(const scalar p, const scalarField& Y) const;
 };