Restore pT equilibrium on phase depletion

When the pTg relaxation depletes one of the reacting phases, the Gibbs
equality no longer constrains the state, so the pT solver runs again
on the final masses. The seven Phase Change model 6 goldens and the
phasechange_bubble example change.

Fixes #1820.

Author: Billy Cao

src/common/m_phase_change.fpp

@@ -136,6 +136,12 @@ contains
                         q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = m2
 
                         call s_infinite_ptg_relaxation_k(j, k, l, pS, rhoe, q_cons_vf, TS)
+
+                        ! A depleted phase leaves Gibbs equality unconstrained; solve pT equilibrium again for the final masses
+                        if (q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, &
+                            & l) == 0.0_wp .or. q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) == 0.0_wp) then
+                            call s_infinite_pt_relaxation_k(j, k, l, 2, pS, p_infpT, q_cons_vf, rhoe, TS)
+                        end if
                     end if
 
                     ! Calculations AFTER equilibrium