diff --git a/Source/mass.f90 b/Source/mass.f90
index f92260f7728..3c7975c6e72 100644
--- a/Source/mass.f90
+++ b/Source/mass.f90
@@ -630,7 +630,7 @@ SUBROUTINE CHECK_MASS_DENSITY
 REAL(EB) :: MASS_N(-3:3),CONST,MASS_C,RHO_ZZ_CUT,RHO_CUT,VC(-3:3),SIGN_FACTOR,SUM_MASS_N,VC1(-3:3),&
             RHO_ZZ_MIN,RHO_ZZ_MAX,SUM_RHO_ZZ,RHO_ZZ_TEST
 INTEGER :: IC
-LOGICAL :: CLIP_RHO_ZZ
+LOGICAL :: CLIP_RHO_ZZ(N_TRACKED_SPECIES)
 REAL(EB), POINTER, DIMENSION(:,:,:) :: DELTA_RHO,DELTA_RHO_ZZ,RHOP
 REAL(EB), POINTER, DIMENSION(:,:,:,:) :: RHO_ZZ
 
@@ -716,12 +716,12 @@ SUBROUTINE CHECK_MASS_DENSITY
 ! Correct species mass density
 
 RHO_ZZ_MIN = 0._EB
+CLIP_RHO_ZZ = .FALSE.
 
 SPECIES_LOOP: DO N=1,N_TRACKED_SPECIES
 
    DELTA_RHO_ZZ => WORK5
    DELTA_RHO_ZZ = 0._EB
-   CLIP_RHO_ZZ = .FALSE.
 
    DO K=1,KBAR
       DO J=1,JBAR
@@ -739,7 +739,7 @@ SUBROUTINE CHECK_MASS_DENSITY
 
             RHO_ZZ_MAX = RHOP(I,J,K)
             IF (RHO_ZZ(I,J,K,N)>=RHO_ZZ_MIN .AND. RHO_ZZ(I,J,K,N)<=RHO_ZZ_MAX) CYCLE
-            CLIP_RHO_ZZ = .TRUE.
+            CLIP_RHO_ZZ(N) = .TRUE.
             IF (RHO_ZZ(I,J,K,N)<RHO_ZZ_MIN) THEN
                RHO_ZZ_CUT = RHO_ZZ_MIN
                SIGN_FACTOR = 1._EB
@@ -777,7 +777,7 @@ SUBROUTINE CHECK_MASS_DENSITY
       ENDDO
    ENDDO
 
-   IF (.NOT.CLIP_RHO_ZZ) CYCLE
+   IF (.NOT.CLIP_RHO_ZZ(N)) CYCLE
 
    ! Assign excess/deficit RHO_ZZ neighboring cells
 
@@ -793,7 +793,7 @@ SUBROUTINE CHECK_MASS_DENSITY
 
 ! If nothing has been clipped, return
 
-IF (.NOT.CLIP_RHOMIN .AND. .NOT.CLIP_RHOMAX .AND. .NOT.CLIP_RHO_ZZ) RETURN
+IF (.NOT.CLIP_RHOMIN .AND. .NOT.CLIP_RHOMAX .AND. .NOT. ANY(CLIP_RHO_ZZ)) RETURN
 
 ! Final check of RHO_ZZ to ensure that ZZ(:,:,:,1:N_TRACKED_SPECIES) sums to 1