diff --git a/src_compressible/compressible_functions.cpp b/src_compressible/compressible_functions.cpp
index e2cf9142..3ffe274d 100644
--- a/src_compressible/compressible_functions.cpp
+++ b/src_compressible/compressible_functions.cpp
@@ -3,6 +3,7 @@
 
 AMREX_GPU_MANAGED int compressible::transport_type;
 AMREX_GPU_MANAGED int compressible::membrane_cell;
+AMREX_GPU_MANAGED int compressible::membrane_type;
 AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, MAX_SPECIES> compressible::transmission;
 AMREX_GPU_MANAGED int compressible::do_1D;
 AMREX_GPU_MANAGED int compressible::do_2D;
@@ -39,6 +40,12 @@ void InitializeCompressibleNamespace()
     pp.query("membrane_cell",membrane_cell);
     if (membrane_cell >= 0) amrex::Print() << "Membrane cell is: " << membrane_cell << "\n";
 
+    // get membrane type
+    membrane_type = 0; // membrane type (0: SSA; 1: Langevin)
+    pp.query("membrane_type",membrane_type);
+    if (membrane_type == 0) amrex::Print() << "SSA membrane type" << "\n";
+    if (membrane_type == 1) amrex::Print() << "Langevin membrane type" << "\n";
+
     // get membrane transmission
     for (int i=0; i<MAX_SPECIES; ++i) {
         transmission[i] = 0.0;
diff --git a/src_compressible/compressible_namespace.H b/src_compressible/compressible_namespace.H
index 55d85dc0..b2e04a5a 100644
--- a/src_compressible/compressible_namespace.H
+++ b/src_compressible/compressible_namespace.H
@@ -2,6 +2,7 @@ namespace compressible {
 
     extern AMREX_GPU_MANAGED int transport_type;
     extern AMREX_GPU_MANAGED int membrane_cell;
+    extern AMREX_GPU_MANAGED int membrane_type;
     extern AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, MAX_SPECIES> transmission;
     extern AMREX_GPU_MANAGED int do_1D;
     extern AMREX_GPU_MANAGED int do_2D;
diff --git a/src_compressible_stag/boundaryStag.cpp b/src_compressible_stag/boundaryStag.cpp
index b4c10039..f7bbfa31 100644
--- a/src_compressible_stag/boundaryStag.cpp
+++ b/src_compressible_stag/boundaryStag.cpp
@@ -270,7 +270,7 @@ void setBCStag(MultiFab& prim_in, MultiFab& cons_in,
         Abort("setBC: momentum and velocity need the same number of ghost cells");
     }
 
-    if (membrane_cell >= 0) { // set adiabatic slip BC at the membrane
+    if ((membrane_cell >= 0) and (membrane_type == 1)) { // set adiabatic slip BC at the membrane (Langevin)
         BCMem(prim_in, cons_in, cumom_in, vel_in, geom);
     }
 
diff --git a/src_compressible_stag/compressible_functions_stag.H b/src_compressible_stag/compressible_functions_stag.H
index f68d9e0e..8b123d0c 100644
--- a/src_compressible_stag/compressible_functions_stag.H
+++ b/src_compressible_stag/compressible_functions_stag.H
@@ -150,6 +150,27 @@ void applyEffusion(std::array<MultiFab, AMREX_SPACEDIM>& faceflux,
                    MultiFab& cons_in,
                    std::array< MultiFab, AMREX_SPACEDIM >& cumom);
 
+void 
+ResetMembraneMom(std::array<MultiFab, AMREX_SPACEDIM>& cumom,
+                 const std::array<MultiFab, AMREX_SPACEDIM>& cumom_mem,
+                 const amrex::Geometry& geom);
+
+void 
+ResetMembraneFluxes(const std::array<MultiFab, AMREX_SPACEDIM>& faceflux_mem,
+                    std::array<MultiFab, AMREX_SPACEDIM>& faceflux,
+                    std::array< MultiFab, 2 >& edgeflux_x,
+                    std::array< MultiFab, 2 >& edgeflux_y,
+                    std::array< MultiFab, 2 >& edgeflux_z,
+                    const amrex::Geometry& geom);
+
+void 
+ComputeFluxMomMembrane(const MultiFab& cons0_in, const MultiFab& prim0_in,
+                       const std::array<MultiFab, AMREX_SPACEDIM>& vel0,
+                       std::array<MultiFab, AMREX_SPACEDIM>& cumom_mem,
+                       std::array<MultiFab, AMREX_SPACEDIM>& faceflux_mem,
+                       const amrex::Geometry& geom,
+                       const amrex::Real dt);
+
 void 
 ComputeFluxMomReservoir(const MultiFab& cons0_in, const MultiFab& prim0_in,
                         const std::array<MultiFab, AMREX_SPACEDIM>& vel0,
diff --git a/src_compressible_stag/fluxStag.cpp b/src_compressible_stag/fluxStag.cpp
index c733d2b3..1bf2ef96 100644
--- a/src_compressible_stag/fluxStag.cpp
+++ b/src_compressible_stag/fluxStag.cpp
@@ -806,7 +806,7 @@ void calculateFluxStag(const MultiFab& cons_in, const std::array< MultiFab, AMRE
 
         // Enforce flux boundary conditions
         StochFluxStag(faceflux_in,cenflux_in,edgeflux_x_in,edgeflux_y_in,edgeflux_z_in,geom);
-        if (membrane_cell >= 0) {
+        if ((membrane_cell >= 0) and (membrane_type == 1)) {
             StochFluxMem(faceflux_in,edgeflux_x_in,edgeflux_y_in,edgeflux_z_in);
         }
 
diff --git a/src_compressible_stag/main_driver.cpp b/src_compressible_stag/main_driver.cpp
index af1bcb18..e9d6c24d 100644
--- a/src_compressible_stag/main_driver.cpp
+++ b/src_compressible_stag/main_driver.cpp
@@ -219,6 +219,7 @@ void main_driver(const char* argv)
     if ((all_correl == 1) and (cross_cell > 0) and (cross_cell < n_cells[0]-1)) {
         amrex::Print() << "Correlations will be done at four equi-distant x* because all_correl = 1" << "\n";
     }
+    if ((membrane_cell > 0) and ((membrane_type < 0) or (membrane_type > 1))) Abort("membrane_type needs to be 0 (SSA) or 1 (Langevin)");
 
     // contains yz-averaged running & instantaneous averages of conserved variables (2*nvars) + primitive variables [vx, vy, vz, T, Yk]: 2*4 + 2*nspecies 
     Vector<Real> dataSliceMeans_xcross(2*nvars+8+2*nspecies, 0.0); 
diff --git a/src_compressible_stag/membraneStag.cpp b/src_compressible_stag/membraneStag.cpp
index 4ed2d0a3..87903633 100644
--- a/src_compressible_stag/membraneStag.cpp
+++ b/src_compressible_stag/membraneStag.cpp
@@ -2,6 +2,7 @@
 #include "compressible_functions.H"
 #include "common_functions.H"
 #include "rng_functions.H"
+#include "reservoirStag_K.H"
 #include <math.h>
 
 void doMembraneStag(MultiFab& cons, 
@@ -470,3 +471,239 @@ void applyEffusion(std::array<MultiFab, AMREX_SPACEDIM>& faceflux,
     }
 
 }
+
+///////////////////////////////////////////////////////////////////////////
+// compute momentum and fluxes at the membrane ////////////////////////////
+// for SSA-type membrane; for now membrane is always perpendicular to x ///
+///////////////////////////////////////////////////////////////////////////
+void 
+ComputeFluxMomMembrane(const MultiFab& cons0_in, const MultiFab& prim0_in,
+                       const std::array<MultiFab, AMREX_SPACEDIM>& vel0,
+                       std::array<MultiFab, AMREX_SPACEDIM>& cumom_mem,
+                       std::array<MultiFab, AMREX_SPACEDIM>& faceflux_mem,
+                       const amrex::Geometry& geom,
+                       const amrex::Real dt)
+{
+    BL_PROFILE_VAR("ComputeFluxMomMembrane()",ComputeFluxMomMembrane);
+
+    const GpuArray<Real, AMREX_SPACEDIM> dx = geom.CellSizeArray();
+    Box dom(geom.Domain());
+
+    Real N_A = Runiv/k_B; // Avagadro's number
+    
+    GpuArray<Real,MAX_SPECIES> mass;
+    for (int l=0;l<nspecies;++l) {
+        mass[l] = molmass[l]/(N_A);
+    }
+
+    for (int d=0; d<AMREX_SPACEDIM; ++d) {
+        cumom_mem[d].setVal(0.0);
+        faceflux_mem[d].setVal(0.0);
+    }
+
+    Real area, vol;
+    area = dx[1]*dx[2];
+    vol  = dx[0]*dx[1]*dx[2];
+
+    // face-based flux (mass and energy) and normal momentum /////
+    //////////////////////////////////////////////////////////////
+    for (MFIter mfi(faceflux_mem[0]); mfi.isValid(); ++mfi) {
+        
+        const Box& bx = mfi.validbox();
+        
+        const Array4<const Real> cons0   = cons0_in.array(mfi);
+        const Array4<const Real> prim0   = prim0_in.array(mfi);
+        const Array4<const Real>& xvel0  = (vel0[0]).array(mfi);
+        const Array4<const Real>& yvel0  = (vel0[1]).array(mfi);
+        const Array4<const Real>& zvel0  = (vel0[2]).array(mfi);
+        
+        const Array4<Real>& xmom  = (cumom_mem[0]).array(mfi);
+        const Array4<Real>& xflux = (faceflux_mem[0]).array(mfi);
+
+        // membrane on hi side of bx
+        if ((membrane_cell >= bx.smallEnd(0)) and (membrane_cell <= bx.bigEnd(0))) {
+            amrex::ParallelForRNG(bx, 
+            [=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::RandomEngine const& engine) noexcept
+            {
+                if (i == membrane_cell) {
+                    ////////////////// left cell ///////////////////////
+                    
+                    ////////////////////////////////////////////////////
+                    // number of particles in FHD cell (for correction)
+                    //Real N = 0.0;
+                    //GpuArray<Real,MAX_SPECIES> N_i;
+                    //for (int n=0;n<nspecies;++n) {
+                    //    N_i[n] = vol*(cons0(i-1,j,k,5+n)/mass[n]);
+                    //    N += N_i[n];
+                    //}
+                    ////////////////////////////////////////////////////
+                    
+                    Real T = prim0(i-1,j,k,4);
+                    Real Vx, Vy, Vz;
+                    Vx = 0.5*(xvel0(i-1,j,k) + xvel0(i,j,k)); // normal
+//                    Vx = 0.0;
+                    Vy = 0.5*(yvel0(i-1,j,k) + yvel0(i-1,j+1,k)); // tangential
+                    Vz = 0.5*(zvel0(i-1,j,k) + zvel0(i-1,j,k+1)); // tangential
+                    Real mass_cross; // total mass into membrane
+                    GpuArray<Real,3> mom_cross; // total momentum into membrane
+                    Real en_cross; // total energy into membrane
+                    GpuArray<Real,MAX_SPECIES> spec_mass_cross;
+                    GpuArray<Real,MAX_SPECIES> rhoYk;
+                    for (int n=0;n<nspecies;++n) {
+                        rhoYk[n] = cons0(i-1,j,k,5+n);
+                    }
+                    if (do_1D) {
+                        poisson_process_membrane(mass,rhoYk,T,Vx,nspecies,area,k_B,dt,
+                                                 mass_cross,mom_cross,en_cross,spec_mass_cross,1,Vy,Vz,
+                                                 transmission,engine);
+                    }
+                    else if (do_2D) {
+                        poisson_process_membrane(mass,rhoYk,T,Vx,nspecies,area,k_B,dt,
+                                                 mass_cross,mom_cross,en_cross,spec_mass_cross,2,Vy,Vz,
+                                                 transmission,engine);
+                    }
+                    else {
+                        poisson_process_membrane(mass,rhoYk,T,Vx,nspecies,area,k_B,dt,
+                                                 mass_cross,mom_cross,en_cross,spec_mass_cross,3,Vy,Vz,
+                                                 transmission,engine);
+                    }
+                    xflux(i,j,k,0) += mass_cross/(dt*area); // update mass flux
+                    xflux(i,j,k,4) += en_cross/(dt*area); // update energy flux
+                    for (int n=0;n<nspecies;++n) {
+                        xflux(i,j,k,5+n) += spec_mass_cross[n]/(dt*area); // update species flux
+                    }
+                    
+                    if (do_1D) {
+                        xmom(i,j,k) += mom_cross[0]/vol;
+                    }
+                    else if (do_2D) {
+                        xmom(i,j,k) += mom_cross[0]/vol;
+//                        xflux(i,j,k,1) += mom_cross[1]/dt/area; // tangential flux
+                        xflux(i,j,k,1) += 0.0;
+                    }
+                    else {
+                        xmom(i,j,k) += mom_cross[0]/vol;
+//                        xflux(i,j,k,1) += mom_cross[1]/dt/area; // tangential flux
+//                        xflux(i,j,k,2) += mom_cross[2]/dt/area; // tangential flux
+                        xflux(i,j,k,1) += 0.0;
+                        xflux(i,j,k,2) += 0.0;
+                    }
+                    
+                    ////////////////////////////////////////////////////
+
+                    ////////////////// right cell //////////////////////
+
+                    T = prim0(i,j,k,4);
+                    Vx = -0.5*(xvel0(i,j,k)+xvel0(i+1,j,k)); // normal
+//                    Vx = 0.0;
+                    Vy =  0.5*(yvel0(i,j,k)+yvel0(i,j+1,k)); // tangential
+                    Vz =  0.5*(zvel0(i,j,k)+zvel0(i,j,k+1)); // tangential
+                    for (int n=0;n<nspecies;++n) {
+                        rhoYk[n] = cons0(i,j,k,5+n);
+                    }
+                    if (do_1D) {
+                        poisson_process_membrane(mass,rhoYk,T,Vx,nspecies,area,k_B,dt,
+                                                 mass_cross,mom_cross,en_cross,spec_mass_cross,1,Vy,Vz,
+                                                 transmission,engine);
+                    }
+                    else if (do_2D) {
+                        poisson_process_membrane(mass,rhoYk,T,Vx,nspecies,area,k_B,dt,
+                                                 mass_cross,mom_cross,en_cross,spec_mass_cross,2,Vy,Vz,
+                                                 transmission,engine);
+                    }
+                    else {
+                        poisson_process_membrane(mass,rhoYk,T,Vx,nspecies,area,k_B,dt,
+                                                 mass_cross,mom_cross,en_cross,spec_mass_cross,3,Vy,Vz,
+                                                 transmission,engine);
+                    }
+                    xflux(i,j,k,0) -= mass_cross/(dt*area); // update mass flux
+                    xflux(i,j,k,4) -= en_cross/(dt*area); // update energy flux
+                    for (int n=0;n<nspecies;++n) {
+                        xflux(i,j,k,5+n) -= spec_mass_cross[n]/(dt*area); // update species flux
+                    }
+                    
+                    if (do_1D) {
+                        xmom(i,j,k) -= mom_cross[0]/vol;
+                    }
+                    else if (do_2D) {
+                        xmom(i,j,k) -= mom_cross[0]/vol;
+//                        xflux(i,j,k,1) -= mom_cross[1]/dt/area; // tangential flux
+                        xflux(i,j,k,1) -= 0.0;
+                    }
+                    else {
+                        xmom(i,j,k) -= mom_cross[0]/vol;
+//                        xflux(i,j,k,1) -= mom_cross[1]/dt/area; // tangential flux
+//                        xflux(i,j,k,2) -= mom_cross[2]/dt/area; // tangential flux
+                        xflux(i,j,k,1) -= 0.0;
+                        xflux(i,j,k,2) -= 0.0;
+                    }
+                    ////////////////////////////////////////////////////
+                }
+            });
+        }
+    }
+
+    cumom_mem[0].OverrideSync(geom.periodicity());
+    faceflux_mem[0].OverrideSync(geom.periodicity());
+}
+
+///////////////////////////////////////////////////////////////////////////
+// reset fluxes at the membrane  //////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////
+void 
+ResetMembraneFluxes(const std::array<MultiFab, AMREX_SPACEDIM>& faceflux_mem,
+                    std::array<MultiFab, AMREX_SPACEDIM>& faceflux,
+                    std::array< MultiFab, 2 >& edgeflux_x,
+                    std::array< MultiFab, 2 >& edgeflux_y,
+                    std::array< MultiFab, 2 >& edgeflux_z,
+                    const amrex::Geometry& geom)
+{
+    BL_PROFILE_VAR("ResetMembraneFluxes()",ResetMembraneFluxes);
+
+    for (MFIter mfi(faceflux[0]); mfi.isValid(); ++mfi) {
+        const Box& bx = mfi.validbox();
+        const Array4<Real>& xflux           = (faceflux[0]).array(mfi);
+        const Array4<const Real>& xflux_mem = (faceflux_mem[0]).array(mfi);
+
+        if ((membrane_cell >= bx.smallEnd(0)) and (membrane_cell <= bx.bigEnd(0))) {
+            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                if (i == membrane_cell) {
+                    xflux(i,j,k,0) = xflux_mem(i,j,k,0); // mass flux
+                    xflux(i,j,k,4) = xflux_mem(i,j,k,4); // energy flux
+                    for (int n=0;n<nspecies;++n) {
+                        xflux(i,j,k,5+n) = xflux_mem(i,j,k,5+n); // species flux
+                    }
+                }
+            });
+        }
+    }
+}
+
+///////////////////////////////////////////////////////////////////////////
+// reset momentum at the membrane /////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////
+void 
+ResetMembraneMom(std::array<MultiFab, AMREX_SPACEDIM>& cumom,
+                 const std::array<MultiFab, AMREX_SPACEDIM>& cumom_mem,
+                 const amrex::Geometry& geom)
+{
+    BL_PROFILE_VAR("ResetMembraneMom()",ResetMembraneMom);
+
+    for (MFIter mfi(cumom[0]); mfi.isValid(); ++mfi) {
+        const Box& bx = mfi.validbox();
+        const Array4<Real>& xmom           = (cumom[0]).array(mfi);
+        const Array4<const Real>& xmom_mem = (cumom_mem[0]).array(mfi);
+
+        if ((membrane_cell >= bx.smallEnd(0)) and (membrane_cell <= bx.bigEnd(0))) {
+            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                if (i == membrane_cell) {
+                    xmom(i,j,k) = xmom_mem(i,j,k);
+//                    xmom(i,j,k) = 0.0;
+                }
+            });
+        }
+    }
+}
+
diff --git a/src_compressible_stag/reservoirStag_K.H b/src_compressible_stag/reservoirStag_K.H
index 213cfd92..b1aa962b 100644
--- a/src_compressible_stag/reservoirStag_K.H
+++ b/src_compressible_stag/reservoirStag_K.H
@@ -158,4 +158,124 @@ poisson_process_reservoir(const amrex::GpuArray<amrex::Real,MAX_SPECIES>& mass,
     }
 }
 
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void
+poisson_process_membrane(const amrex::GpuArray<amrex::Real,MAX_SPECIES>& mass,
+                         const amrex::GpuArray<amrex::Real,MAX_SPECIES>& rhoYk,
+                         const amrex::Real& T, // temperature
+                         const amrex::Real& Vnormal, // normal velocity
+                         const int& nspec,
+                         const amrex::Real& area,
+                         const amrex::Real& kB,
+                         const amrex::Real& dt,
+                         amrex::Real& mass_cross,
+                         amrex::GpuArray<amrex::Real,3>& mom_cross,
+                         amrex::Real& en_cross,
+                         amrex::GpuArray<amrex::Real,MAX_SPECIES>& spec_mass_cross,
+                         int dim,         // dimension
+                         amrex::Real VT1, // tangential velocity 1
+                         amrex::Real VT2, // tangential velocity 2
+                         const amrex::GpuArray<amrex::Real, MAX_SPECIES> transmission,
+                         amrex::RandomEngine const& engine)
+
+{
+    constexpr amrex::Real small_number = 1.0e-13;
+    constexpr amrex::Real sqrtPI = 1.77245385091;
+
+    mass_cross = 0.0;
+    en_cross   = 0.0;
+    for (int n=0;n<3;++n) {
+        mom_cross[n] = 0.0;
+    }
+    for (int n=0;n<nspecies;++n) {
+        spec_mass_cross[n] = 0.0;
+    }
+
+    amrex::GpuArray<amrex::Real,MAX_SPECIES> particle_rate;
+    amrex::GpuArray<amrex::Real,MAX_SPECIES> timer;
+    amrex::GpuArray<amrex::Real,MAX_SPECIES> thermal_speed;
+    amrex::GpuArray<amrex::Real,MAX_SPECIES> speed_ratio;
+    amrex::GpuArray<bool,MAX_SPECIES> stop_flag;
+    
+    for (int n=0;n<nspec;++n) {
+        thermal_speed[n] = sqrt(2.0*kB*T/mass[n]);
+        speed_ratio[n] = Vnormal/thermal_speed[n];
+        amrex::Real num_density = rhoYk[n]/mass[n];
+        if (std::abs(speed_ratio[n]) > small_number) {
+            particle_rate[n] = transmission[n]*(1.0/(2.0*sqrtPI))*num_density*area*thermal_speed[n]*
+                (exp(-1.0*speed_ratio[n]*speed_ratio[n]) + speed_ratio[n]*sqrtPI*(1.0 + erf(speed_ratio[n])));
+            if (particle_rate[n] < 0.0) {
+                amrex::Abort("negative particle rate encountered");
+            }
+        }
+        else {
+            particle_rate[n] = transmission[n]*(1.0/(2.0*sqrtPI))*num_density*area*thermal_speed[n];
+            if (particle_rate[n] < 0.0) {
+                amrex::Abort("negative particle rate encountered");
+            }
+        }
+        timer[n] = 0.0;
+        stop_flag[n] = false;
+    }
+
+    // Run Poisson Process
+    bool stop = false;
+    while(!stop) {
+        
+        // add a particle for each species (if stop[n] = false)
+        amrex::Real vel_rand;
+        for (int n=0;n<nspec;++n) {
+
+            if (transmission[n] <= 1.0e-12) {
+                mass_cross += 0.0;
+                spec_mass_cross[n] += 0.0;
+                en_cross += 0.0;
+                mom_cross[0] += 0.0;
+                stop_flag[n] = true;
+            }
+
+            else {
+                if (!stop_flag[n]) {
+                    
+                    // normal component
+                    mass_cross += mass[n]; // add a particle crossing for total mass
+                    spec_mass_cross[n] += mass[n]; // add a particle crossing for species mass
+                    vel_rand = random_cross_vel(speed_ratio[n],thermal_speed[n],false,0.0,engine); // get random velocity
+                    mom_cross[0] += mass[n]*vel_rand; // add a particle crossing for total momentum
+                    en_cross  += 0.5*mass[n]*(vel_rand)*(vel_rand); // add a particle crossing for total energy
+
+                    // tangential component 1
+                    if (dim >= 2) {
+//                        vel_rand = random_cross_vel(speed_ratio[n],thermal_speed[n],true,VT1,engine); // get random velocity
+//                        mom_cross[1] += mass[n]*vel_rand; // add a particle crossing for total momentum
+//                        en_cross  += 0.5*mass[n]*(vel_rand)*(vel_rand); // add a particle crossing for total energy
+                        mom_cross[1] += 0.0;
+                        en_cross  += 0.0;
+                    }
+
+                    // tangential component 2
+                    if (dim == 3) {
+                        //vel_rand = random_cross_vel(speed_ratio[n],thermal_speed[n],true,VT2,engine); // get random velocity
+                        //mom_cross[2] += mass[n]*vel_rand; // add a particle crossing for total momentum
+                        //en_cross  += 0.5*mass[n]*(vel_rand)*(vel_rand); // add a particle crossing for total energy
+                        mom_cross[1] += 0.0;
+                        en_cross  += 0.0;
+                    }
+                }
+            }
+        }
+
+        stop = true;
+        for (int n=0;n<nspec;++n) { // update dt for each species
+            if (!stop_flag[n]) {
+                amrex::Real rn   = amrex::Random(engine);
+                timer[n]         = timer[n] - std::log(rn)/(particle_rate[n]);
+                if (timer[n] >= dt) stop_flag[n] = true;
+            }
+            stop = stop and stop_flag[n];
+        }
+    }
+}
+
+
 #endif
diff --git a/src_compressible_stag/timeStepStag.cpp b/src_compressible_stag/timeStepStag.cpp
index e6663567..07062d47 100644
--- a/src_compressible_stag/timeStepStag.cpp
+++ b/src_compressible_stag/timeStepStag.cpp
@@ -36,12 +36,27 @@ void RK3stepStag(MultiFab& cu,
     // Reservoir stuff
     std::array< MultiFab, AMREX_SPACEDIM > cumom_res; // MFab for storing momentum from reservoir update
     std::array< MultiFab, AMREX_SPACEDIM > faceflux_res; // MFab for storing fluxes (face-based) from reservoir update
-    for (int d=0; d<AMREX_SPACEDIM; ++d) {
-        cumom_res[d].define(convert(cu.boxArray(),nodal_flag_dir[d]), cu.DistributionMap(), 1, 0);
-        faceflux_res[d].define(convert(cu.boxArray(),nodal_flag_dir[d]), cu.DistributionMap(), nvars, 0);
+    if (do_reservoir) {
+        for (int d=0; d<AMREX_SPACEDIM; ++d) {
+            cumom_res[d].define(convert(cu.boxArray(),nodal_flag_dir[d]), cu.DistributionMap(), 1, 0);
+            faceflux_res[d].define(convert(cu.boxArray(),nodal_flag_dir[d]), cu.DistributionMap(), nvars, 0);
+        }
     }
     // Reservoir stuff
     ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Membrane stuff
+    std::array< MultiFab, AMREX_SPACEDIM > cumom_mem; // MFab for storing momentum from membrane update
+    std::array< MultiFab, AMREX_SPACEDIM > faceflux_mem; // MFab for storing fluxes (face-based) from membrane update
+    if ((membrane_cell>=0) and (membrane_type == 0)) {
+        for (int d=0; d<AMREX_SPACEDIM; ++d) {
+            cumom_mem[d].define(convert(cu.boxArray(),nodal_flag_dir[d]), cu.DistributionMap(), 1, 0);
+            faceflux_mem[d].define(convert(cu.boxArray(),nodal_flag_dir[d]), cu.DistributionMap(), nvars, 0);
+        }
+    }
+    // Membrane stuff
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
 
     std::array< MultiFab, AMREX_SPACEDIM > cupmom;
     std::array< MultiFab, AMREX_SPACEDIM > cup2mom;
@@ -335,8 +350,12 @@ void RK3stepStag(MultiFab& cu,
         geom, stoch_weights,dt);
 
     // reservoir timers
-    Real aux1, aux2, aux3, aux4, aux5, aux6; 
+    Real aux1, aux2, aux3, aux4, aux5, aux6;
+    // membrane timers
+    Real auxm1, auxm2, auxm3, auxm4, auxm5, auxm6;
     
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Reservoir stuff
     // add to the total continuum fluxes based on RK3 weight
     if (do_reservoir) {
         aux1 = ParallelDescriptor::second();
@@ -350,6 +369,26 @@ void RK3stepStag(MultiFab& cu,
         aux2 = ParallelDescriptor::second() - aux1;
         ParallelDescriptor::ReduceRealMax(aux2,  ParallelDescriptor::IOProcessorNumber());
     }
+    // Reservoir stuff
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Membrane stuff (SSA)
+    // add to the total continuum fluxes based on RK3 weight
+    if ((membrane_cell>=0) and (membrane_type == 0)) {
+        auxm1 = ParallelDescriptor::second();
+        ComputeFluxMomMembrane(cu,prim,vel,cumom_mem,
+                               faceflux_mem,geom,dt); // compute fluxes and momentum from membrane particle update
+        ResetMembraneFluxes(faceflux_mem, faceflux,
+                            edgeflux_x,
+                            edgeflux_y,
+                            edgeflux_z,
+                            geom); // reset fluxes at the membrane from particle update
+        auxm2 = ParallelDescriptor::second() - auxm1;
+        ParallelDescriptor::ReduceRealMax(auxm2,  ParallelDescriptor::IOProcessorNumber());
+    }
+    // Membrane stuff (SSA)
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
 
     if (nreaction>0) {
         MultiFab::LinComb(ranchem,
@@ -535,9 +574,21 @@ void RK3stepStag(MultiFab& cu,
         BCMomTrans(cupmom[i], vel[i], geom, i);
     }
 
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Reservoir stuff
     if (do_reservoir) {
         ResetReservoirMom(cupmom, cumom_res, geom); // set momentum at the reservoir interface to its value from particle update
     }
+    // Reservoir stuff
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Membrane stuff (SSA)
+    if ((membrane_cell>=0) and (membrane_type == 0)) {
+        ResetMembraneMom(cupmom, cumom_mem, geom); // set momentum at the membrane to its value from particle update
+    }
+    // Membrane stuff (SSA)
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
 
     // Fill boundaries for conserved variables
     for (int d=0; d<AMREX_SPACEDIM; d++) {
@@ -663,6 +714,8 @@ void RK3stepStag(MultiFab& cu,
         stochface, stochedge_x, stochedge_y, stochedge_z, stochcen, 
         geom, stoch_weights,dt);
 
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Reservoir stuff
     // add to the total continuum fluxes based on RK3 weight
     if (do_reservoir) {
         aux3 = ParallelDescriptor::second();
@@ -677,6 +730,27 @@ void RK3stepStag(MultiFab& cu,
         aux4 = ParallelDescriptor::second() - aux3;
         ParallelDescriptor::ReduceRealMax(aux4,  ParallelDescriptor::IOProcessorNumber());
     }
+    // Reservoir stuff
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Membrane stuff (SSA)
+    // add to the total continuum fluxes based on RK3 weight
+    if ((membrane_cell>=0) and (membrane_type == 0)) {
+        auxm3 = ParallelDescriptor::second();
+        ComputeFluxMomMembrane(cup,prim,vel,cumom_mem,
+                               faceflux_mem,
+                               geom,0.25*dt); // compute fluxes and momentum from membrane particle update
+        ResetMembraneFluxes(faceflux_mem, faceflux,
+                            edgeflux_x,
+                            edgeflux_y,
+                            edgeflux_z,
+                            geom); // reset fluxes at the membrane from particle update
+        auxm4 = ParallelDescriptor::second() - auxm3;
+        ParallelDescriptor::ReduceRealMax(auxm4,  ParallelDescriptor::IOProcessorNumber());
+    }
+    // Membrane stuff (SSA)
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
 
     if (nreaction>0) {
         MultiFab::LinComb(ranchem,
@@ -869,9 +943,21 @@ void RK3stepStag(MultiFab& cu,
         BCMomTrans(cup2mom[i], vel[i], geom, i);
     }
 
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Reservoir stuff
     if (do_reservoir) {
         ResetReservoirMom(cup2mom, cumom_res, geom); // set momentum at the reservoir interface to its value from particle update
     }
+    // Reservoir stuff
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Membrane stuff (SSA)
+    if ((membrane_cell>=0) and (membrane_type == 0)) {
+        ResetMembraneMom(cup2mom, cumom_mem, geom); // set momentum at the membrane to its value from particle update
+    }
+    // Membrane stuff (SSA)
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
 
     // Fill  boundaries for conserved variables
     for (int d=0; d<AMREX_SPACEDIM; d++) {
@@ -997,6 +1083,8 @@ void RK3stepStag(MultiFab& cu,
         stochface, stochedge_x, stochedge_y, stochedge_z, stochcen, 
         geom, stoch_weights,dt);
     
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Reservoir stuff
     // add to the total continuum fluxes based on RK3 weight
     if (do_reservoir) {
         aux5 = ParallelDescriptor::second();
@@ -1011,6 +1099,27 @@ void RK3stepStag(MultiFab& cu,
         aux6 = ParallelDescriptor::second() - aux5;
         ParallelDescriptor::ReduceRealMax(aux6,  ParallelDescriptor::IOProcessorNumber());
     }
+    // Reservoir stuff
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Membrane stuff (SSA)
+    // add to the total continuum fluxes based on RK3 weight
+    if ((membrane_cell>=0) and (membrane_type == 0)) {
+        auxm5 = ParallelDescriptor::second();
+        ComputeFluxMomMembrane(cup2,prim,vel,cumom_mem,
+                               faceflux_mem,
+                               geom,(2.0/3.0)*dt); // compute fluxes and momentum from membrane particle update
+        ResetMembraneFluxes(faceflux_mem, faceflux,
+                            edgeflux_x,
+                            edgeflux_y,
+                            edgeflux_z,
+                            geom); // reset fluxes at the membrane from particle update
+        auxm6 = ParallelDescriptor::second() - auxm5;
+        ParallelDescriptor::ReduceRealMax(auxm6,  ParallelDescriptor::IOProcessorNumber());
+    }
+    // Membrane stuff (SSA)
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
 
     if (nreaction>0) {
         MultiFab::LinComb(ranchem,
@@ -1210,9 +1319,21 @@ void RK3stepStag(MultiFab& cu,
         BCMomTrans(cumom[i], vel[i], geom, i);
     }
 
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Reservoir stuff
     if (do_reservoir) {
         ResetReservoirMom(cumom, cumom_res, geom); // set momentum at the reservoir interface to its value from particle update
     }
+    // Reservoir stuff
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Membrane stuff (SSA)
+    if ((membrane_cell>=0) and (membrane_type == 0)) {
+        ResetMembraneMom(cumom, cumom_mem, geom); // set momentum at the membrane to its value from particle update
+    }
+    // Membrane stuff (SSA)
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
     
     // Fill  boundaries for conserved variables
     for (int d=0; d<AMREX_SPACEDIM; d++) {
@@ -1230,17 +1351,34 @@ void RK3stepStag(MultiFab& cu,
     prim.FillBoundary(geom.periodicity());
     cu.FillBoundary(geom.periodicity()); 
 
-    // Membrane setup
-    if (membrane_cell >= 0) {
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Membrane stuff (Langevin)
+    if ((membrane_cell >= 0) and (membrane_type == 1)) {
         doMembraneStag(cu,cumom,prim,vel,faceflux,edgeflux_x,geom,dt);
     }
+    // Membrane stuff (Langevin)
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
 
     // Correctly set momentum and velocity at the walls & temperature, pressure, density & mass/mole fractions in ghost cells
     setBCStag(prim, cu, cumom, vel, geom);
 
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Reservoir stuff
     if (do_reservoir) {
         if (step%100 == 0) {
             amrex::Print() << "Step: " << step << " Reservoir generator time: " << aux2 + aux4 + aux6 << " seconds\n";
         }
     }
+    // Reservoir stuff
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Membrane stuff
+    if ((membrane_cell >= 0) and (membrane_type == 0)) {
+        if (step%100 == 0) {
+            amrex::Print() << "Step: " << step << " Membrane (SSA) time: " << auxm2 + auxm4 + auxm6 << " seconds\n";
+        }
+    }
+    // Membrane
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////
 }