diff --git a/incompact3d.f90 b/incompact3d.f90
index f220fd7..9df8755 100644
--- a/incompact3d.f90
+++ b/incompact3d.f90
@@ -266,7 +266,7 @@ PROGRAM incompact3d
       poissiter = 0
       do while(converged.eqv..FALSE.)
         if (ilmn.ne.0) then
-          if (nrhoscheme.eq.0) then
+          if (ivarcoeff.ne.0) then
             if (poissiter.ne.0) then
               !! Compute correction term
               call divergence_corr(rho1, px1, py1, pz1, ta1, tb1, tc1, td1, te1, tf1, di1, tg1, &
@@ -309,7 +309,7 @@ PROGRAM incompact3d
         
 !!! CM call test_min_max('pp3  ','In main deco   ',pp3,size(pp3))
 
-        if ((nrhoscheme.ne.0).or.(ilmn.eq.0)) then
+        if ((ilmn.eq.0).or.(ivarcoeff.eq.0)) then
           converged = .TRUE.
         endif
 
diff --git a/incompact3d.prm b/incompact3d.prm
index e40249b..99f4ab4 100644
--- a/incompact3d.prm
+++ b/incompact3d.prm
@@ -30,6 +30,7 @@
 0       #istret # y mesh refinement (0:no, 1:center, 2:both sides, 3:bottom)
 0.3     #beta   # Refinement parameter (beta)
 1       #ilmn   # (0:incompressible, 1: Low Mach Number)
+1       #ivarcoeff# (0:constant-coefficient Poisson equation, 1:variable-coefficient Poisson equation)
 1       #npoisscheme# How to compute rho0 in var-coeff poisson equation? (0: rho0 = MIN(rho), Anything else: harmonic average)
 1       #iskew  # (0:urotu, 1:skew, for the convective terms)
 1       #iprops # (0: constant properties, 1: variable properties)
diff --git a/module_param.f90 b/module_param.f90
index 964c915..5224b55 100644
--- a/module_param.f90
+++ b/module_param.f90
@@ -138,7 +138,7 @@ module param
   integer :: ifft, ivirt,istret,iforc_entree,iturb
   integer :: itype, iskew, iin, nscheme, ifirst, ilast, iles
   integer :: isave,ilit,idebmod, imodulo, idemarre, icommence, irecord
-  integer :: ilmn, nrhoscheme, npoissscheme
+  integer :: ilmn, nrhoscheme, npoissscheme, ivarcoeff
   integer :: iscalar
   integer :: iprops
   integer :: nxboite, istat,iread,iadvance_time 
diff --git a/navier.f90 b/navier.f90
index 53344c2..16dc041 100644
--- a/navier.f90
+++ b/navier.f90
@@ -274,7 +274,7 @@ SUBROUTINE inttdensity(rho1, rhos1, rhoss1, rhos01, tg1, drhodt1, ux1, uy1, uz1,
     phi1(:,:,:) = phi1(:,:,:) / rho1(:,:,:)
   endif
   
-  if (nrhoscheme.eq.0) then
+  if (ivarcoeff.ne.0) then
     !---------------------------------------------------------------------------------
     ! XXX Using variable-coefficient Poisson equation, convert momentum back to
     !     velocity.
@@ -317,7 +317,7 @@ subroutine corgp (ux,gx,uy,uz,px,py,pz,rho)
   nxyz=xsize(1)*xsize(2)*xsize(3)
 
   if (ilmn.ne.0) then
-    if (nrhoscheme.ne.0) then
+    if (ivarcoeff.eq.0) then
       !! We are solving constant-coefficient Poisson equation,
       !! first convert momentum->velocity
       if (iskew.ne.2) then
@@ -947,7 +947,7 @@ subroutine divergence (ux1,uy1,uz1,ep1,ta1,tb1,tc1,di1,td1,te1,tf1,&
 
   !WORK X-PENCILS
   call decx6(td1,ta1,di1,sx,cfx6,csx6,cwx6,xsize(1),nxmsize,xsize(2),xsize(3),0)
-  if (nrhoscheme.eq.0) then
+  if (ivarcoeff.ne.0) then
     !! Solving variable-coefficient Poisson equation
     !  Get divu to x pencils and interpolate to pressure points
     call transpose_z_to_y(divu3, divu2)
@@ -1053,14 +1053,6 @@ SUBROUTINE extrapol_rhotrans(rho1, rhos1, rhoss1, rhos01, drhodt1)
   INTEGER :: subitr
 
   INTEGER :: ijk, nxyz
-  INTEGER :: nrhoscheme_tmp
-
-  IF (nrhoscheme.LE.0) THEN
-    nrhoscheme_tmp = nrhoscheme
-    nrhoscheme = 3
-  ELSE
-    nrhoscheme_tmp = nrhoscheme
-  ENDIF
 
   nxyz = xsize(1) * xsize(2) * xsize(3)
 
@@ -1115,8 +1107,6 @@ SUBROUTINE extrapol_rhotrans(rho1, rhos1, rhoss1, rhos01, drhodt1)
     ENDIF
   ENDIF
 
-  nrhoscheme = nrhoscheme_tmp
-
 ENDSUBROUTINE extrapol_rhotrans
 
 !********************************************************************
diff --git a/parameters.f90 b/parameters.f90
index ac628e4..ed5931f 100644
--- a/parameters.f90
+++ b/parameters.f90
@@ -90,6 +90,7 @@ subroutine parameter()
 read (10,*) istret
 read (10,*) beta
 read (10,*) ilmn
+read (10,*) ivarcoeff
 read (10,*) npoissscheme
 read (10,*) iskew
 read (10,*) iprops
@@ -198,7 +199,8 @@ subroutine parameter()
 if (ilmn.ne.0) then
   print *, "Low Mach Number: Enabled"
 
-  if (nrhoscheme.le.0) then
+  if (ivarcoeff.ne.0) then
+    print *, "Var-coeff Poisson: ENABLED"
     print *, "Poisson tolerance: ", tol
 
     if (npoissscheme.eq.0) then