diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp
index d8c067450e..e56c2acc93 100644
--- a/src/common/m_derived_types.fpp
+++ b/src/common/m_derived_types.fpp
@@ -57,6 +57,10 @@ module m_derived_types
         real(kind(0d0)) :: ve1
         real(kind(0d0)) :: ve2
         real(kind(0d0)) :: ve3
+        real(kind(0d0)) :: u_in, v_in, w_in, u_out, v_out, w_out, pres_in, pres_out
+        real(kind(0d0)), dimension(num_fluids_max) :: alpha_rho_in,  alpha_in
+        logical :: grcbc_in, grcbc_out, grcbc_vel_out
+
     end type int_bounds_info
 
     !> Derived type adding beginning (beg) and end bounds info as attributes
diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp
index 5b95b6d92f..bc3e78764f 100644
--- a/src/simulation/m_cbc.fpp
+++ b/src/simulation/m_cbc.fpp
@@ -134,6 +134,31 @@ module m_cbc
     integer :: cbc_dir, cbc_loc
 !$acc declare create(dj, bcxb, bcxe, bcyb, bcye, bczb, bcze, cbc_dir, cbc_loc)
 
+    real(kind(0d0)) :: ux_in, ux_out, vx_in, vx_out, wx_in, wx_out, presx_in, presx_out, Delx_in, Delx_out
+     
+    real(kind(0d0)) :: uy_in, uy_out, vy_in, vy_out, wy_in, wy_out, presy_in, presy_out, Dely_in, Dely_out
+     
+    real(kind(0d0)) :: uz_in, uz_out, vz_in, vz_out, wz_in, wz_out, presz_in, presz_out, Delz_in, Delz_out
+    
+!$acc declare create(ux_in, ux_out, vx_in, vx_out, wx_in, wx_out, presx_in, presx_out, Delx_in, Delx_out)
+!$acc declare create(uy_in, uy_out, vy_in, vy_out, wy_in, wy_out, presy_in, presy_out, Dely_in, Dely_out)
+!$acc declare create(uz_in, uz_out, vz_in, vz_out, wz_in, wz_out, presz_in, presz_out, Delz_in, Delz_out)
+
+
+#ifdef CRAY_ACC_WAR
+    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), alpha_rhox_in,  alphax_in)
+    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), alpha_rhoy_in,  alphay_in)
+    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), alpha_rhoz_in,  alphaz_in)
+!$acc declare link(alpha_rhox_in,  alphax_in, alpha_rhoy_in,  alphay_in, alpha_rhoz_in,  alphaz_in) 
+#else
+    real(kind(0d0)), allocatable, dimension(:) :: alpha_rhox_in,  alphax_in
+    real(kind(0d0)), allocatable, dimension(:) :: alpha_rhoy_in,  alphay_in
+    real(kind(0d0)), allocatable, dimension(:) :: alpha_rhoz_in,  alphaz_in
+!$acc declare create(alpha_rhox_in, alphax_in, alpha_rhoy_in, alphay_in, alpha_rhoz_in, alphaz_in) 
+#endif
+
+
+
 #ifndef CRAY_ACC_WAR
 !$acc declare create(q_prim_rsx_vf, q_prim_rsy_vf, q_prim_rsz_vf,  F_rsx_vf, F_src_rsx_vf,flux_rsx_vf, flux_src_rsx_vf, &
 !$acc                 F_rsy_vf, F_src_rsy_vf,flux_rsy_vf, flux_src_rsy_vf, F_rsz_vf, F_src_rsz_vf,flux_rsz_vf, flux_src_rsz_vf,Re, &
@@ -426,6 +451,58 @@ contains
             !$acc update device(bczb, bcze)
         end if
 
+        @:ALLOCATE_GLOBAL(alpha_rhox_in(1:num_fluids), alpha_rhoy_in(1:num_fluids), alpha_rhoz_in(1:num_fluids))
+        @:ALLOCATE_GLOBAL(alphax_in(1:num_fluids), alphay_in(1:num_fluids), alphaz_in(1:num_fluids))
+
+        #:for CBC_DIR, XYZ in [(1, 'x'), (2, 'y'), (3, 'z')]
+
+            if(${CBC_DIR}$ == 1) then 
+                u${XYZ}$_in = bc_${XYZ}$%u_in 
+                v${XYZ}$_in = bc_${XYZ}$%v_in 
+                w${XYZ}$_in = bc_${XYZ}$%w_in 
+                u${XYZ}$_out = bc_${XYZ}$%u_out
+                v${XYZ}$_out = bc_${XYZ}$%v_out 
+                w${XYZ}$_out = bc_${XYZ}$%w_out
+                Del${XYZ}$_in = maxval(dx)
+                Del${XYZ}$_out = maxval(dx)
+            else if(${CBC_DIR}$ == 2) then 
+
+                u${XYZ}$_in = bc_${XYZ}$%v_in 
+                v${XYZ}$_in = bc_${XYZ}$%u_in 
+                w${XYZ}$_in = bc_${XYZ}$%w_in 
+                u${XYZ}$_out = bc_${XYZ}$%v_out
+                v${XYZ}$_out = bc_${XYZ}$%u_out 
+                w${XYZ}$_out = bc_${XYZ}$%w_out 
+                if(n > 0) then 
+                    Del${XYZ}$_in = maxval(dy)
+                    Del${XYZ}$_out = maxval(dy)
+                end if
+            else if(${CBC_DIR}$ == 3) then 
+                u${XYZ}$_in = bc_${XYZ}$%w_in 
+                v${XYZ}$_in = bc_${XYZ}$%u_in 
+                w${XYZ}$_in = bc_${XYZ}$%v_in 
+                u${XYZ}$_out = bc_${XYZ}$%w_out
+                v${XYZ}$_out = bc_${XYZ}$%u_out 
+                w${XYZ}$_out = bc_${XYZ}$%v_out 
+                if(p > 0) then 
+                    Del${XYZ}$_in = maxval(dz)
+                    Del${XYZ}$_out = maxval(dz)
+                end if
+            end if               
+
+            pres${XYZ}$_in = bc_${XYZ}$%pres_in
+            pres${XYZ}$_out = bc_${XYZ}$%pres_out
+            do i = 1, num_fluids
+                alpha_rho${XYZ}$_in(i) = bc_${XYZ}$%alpha_rho_in(i)
+                alpha${XYZ}$_in(i) = bc_${XYZ}$%alpha_in(i)
+            end do
+            !$acc update device(u${XYZ}$_in, v${XYZ}$_in, w${XYZ}$_in, u${XYZ}$_out, v${XYZ}$_out, w${XYZ}$_out)
+            !$acc update device(pres${XYZ}$_in, pres${XYZ}$_out, alpha_rho${XYZ}$_in, alpha${XYZ}$_in)
+            !$acc update device(Del${XYZ}_in, Del${XYZ}_out)
+            
+        #:endfor
+
+
     end subroutine s_initialize_cbc_module
 
     !>  Compute CBC coefficients
@@ -655,6 +732,7 @@ contains
         real(kind(0d0)) :: pi_inf      !< Cell averaged liquid stiffness
         real(kind(0d0)) :: qv          !< Cell averaged fluid reference energy
         real(kind(0d0)) :: c
+        real(kind(0d0)) :: Ma
 
         real(kind(0d0)) :: vel_K_sum, vel_dv_dt_sum
 
@@ -861,14 +939,42 @@ contains
                         lambda(2) = vel(dir_idx(1))
                         lambda(3) = vel(dir_idx(1)) + c
 
+                        Ma = vel(dir_idx(1)) / c                                            
+
                         if ((cbc_loc == -1 .and. bc${XYZ}$b == -5) .or. (cbc_loc == 1 .and. bc${XYZ}$e == -5)) then
                             call s_compute_slip_wall_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
                         else if ((cbc_loc == -1 .and. bc${XYZ}$b == -6) .or. (cbc_loc == 1 .and. bc${XYZ}$e == -6)) then
                             call s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
                         else if ((cbc_loc == -1 .and. bc${XYZ}$b == -7) .or. (cbc_loc == 1 .and. bc${XYZ}$e == -7)) then
                             call s_compute_nonreflecting_subsonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
+                            
+                            if(bc_${XYZ}$%grcbc_in) then 
+                                !$acc loop seq 
+                                do i = 2, momxb
+                                    L(2) = c**3d0 * Ma * (alpha_rho(i-1) - alpha_rho${XYZ}$_in(i-1)) / Del${XYZ}$_in - c * Ma * (pres - pres${XYZ}$_in) / Del${XYZ}$_in
+                                end do
+                                if(n > 0) then 
+                                    L(momxb + 1) = c * Ma * (vel(dir_idx(2)) - v${XYZ}$_in) / Del${XYZ}$_in
+                                    if(p > 0) then 
+                                        L(momxb + 2) = c * Ma * (vel(dir_idx(3)) - w${XYZ}$_in) / Del${XYZ}$_in
+                                    end if 
+                                end if                                
+                                !$acc loop seq 
+                                do i = E_idx, advxe - 1
+                                    L(i) = c * Ma * (adv(i + 1 - E_idx) - alpha${XYZ}$_in(i + 1 - E_idx)) / Del${XYZ}$_in
+                                end do 
+                                L(advxe) = rho*c**2d0*(1d0 + Ma)*(vel(dir_idx(1)) + u${XYZ}$_in * sign(1, cbc_loc)) / Del${XYZ}$_in + c*(1d0 + Ma)*(pres - pres${XYZ}$_in) / Del${XYZ}$_in
+                            end if      
                         else if ((cbc_loc == -1 .and. bc${XYZ}$b == -8) .or. (cbc_loc == 1 .and. bc${XYZ}$e == -8)) then
                             call s_compute_nonreflecting_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
+                            
+                            if(bc_${XYZ}$%grcbc_out) then
+                                L(advxe) =  c*(1d0 - Ma)*(pres - pres${XYZ}$_out) / Del${XYZ}$_out 
+                                if(bc_${XYZ}$%grcbc_vel_out) then
+                                    L(advxe) = L(advxe) +  rho*c**2d0*(1d0 - Ma)*(vel(dir_idx(1)) + u${XYZ}$_out * sign(1, cbc_loc)) / Del${XYZ}$_out
+                                end if
+                            end if
+
                         else if ((cbc_loc == -1 .and. bc${XYZ}$b == -9) .or. (cbc_loc == 1 .and. bc${XYZ}$e == -9)) then
                             call s_compute_force_free_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
                         else if ((cbc_loc == -1 .and. bc${XYZ}$b == -10) .or. (cbc_loc == 1 .and. bc${XYZ}$e == -10)) then
diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp
index 8797a8f45e..b721834fc5 100644
--- a/src/simulation/m_global_parameters.fpp
+++ b/src/simulation/m_global_parameters.fpp
@@ -473,6 +473,7 @@ module m_global_parameters
 
     !$acc declare create(pb_ts, mv_ts)
 #endif
+
     ! ======================================================================
 
 contains
@@ -680,6 +681,18 @@ contains
             integral(i)%ymax = dflt_real
         end do
 
+        bc_x%grcbc_in = .false. 
+        bc_x%grcbc_out = .false.
+        bc_x%grcbc_vel_out = .false.
+
+        bc_y%grcbc_in = .false. 
+        bc_y%grcbc_out = .false.
+        bc_y%grcbc_vel_out = .false.
+
+        bc_z%grcbc_in = .false. 
+        bc_z%grcbc_out = .false.
+        bc_z%grcbc_vel_out = .false.
+
     end subroutine s_assign_default_values_to_user_inputs
 
     !>  The computation of parameters, the allocation of memory,
diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp
index b395fc624a..b4b02a60cb 100644
--- a/src/simulation/m_mpi_proxy.fpp
+++ b/src/simulation/m_mpi_proxy.fpp
@@ -203,6 +203,9 @@ contains
             & 'polydisperse', 'qbmm', 'acoustic_source', 'probe_wrt', 'integral_wrt',   &
             & 'prim_vars_wrt', 'weno_avg', 'file_per_process', 'relax',          &
             & 'adv_n', 'adap_dt', 'ib', 'bodyForces', 'bf_x', 'bf_y', 'bf_z',    &
+            & 'bc_x%grcbc_in', 'bc_x%grcbc_out', 'bc_x%grcbc_vel_out',          &
+            & 'bc_y%grcbc_in', 'bc_y%grcbc_out', 'bc_y%grcbc_vel_out',          &
+            & 'bc_z%grcbc_in', 'bc_z%grcbc_out', 'bc_z%grcbc_vel_out',          &
             & 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt' ]
             call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
         #:endfor
@@ -210,8 +213,14 @@ contains
         #:for VAR in [ 'dt','weno_eps','teno_CT','pref','rhoref','R0ref','Web','Ca', 'sigma', &
             & 'Re_inv', 'poly_sigma', 'palpha_eps', 'ptgalpha_eps', 'pi_fac',    &
             & 'bc_x%vb1','bc_x%vb2','bc_x%vb3','bc_x%ve1','bc_x%ve2','bc_x%ve2', &
+            & 'bc_x%u_in','bc_x%v_in','bc_x%w_in','bc_x%u_out','bc_x%v_out','bc_x%w_out', &
+            & 'bc_x%pres_in','bc_x%pres_out', &
             & 'bc_y%vb1','bc_y%vb2','bc_y%vb3','bc_y%ve1','bc_y%ve2','bc_y%ve3', &
+            & 'bc_y%u_in','bc_y%v_in','bc_y%w_in','bc_y%u_out','bc_y%v_out','bc_y%w_out', &
+            & 'bc_y%pres_in','bc_y%pres_out', & 
             & 'bc_z%vb1','bc_z%vb2','bc_z%vb3','bc_z%ve1','bc_z%ve2','bc_z%ve3', &
+            & 'bc_z%u_in','bc_z%v_in','bc_z%w_in','bc_z%u_out','bc_z%v_out','bc_z%w_out', &
+            & 'bc_z%pres_in','bc_z%pres_out', & 
             & 'x_domain%beg', 'x_domain%end', 'y_domain%beg', 'y_domain%end',    &
             & 'z_domain%beg', 'z_domain%end', 't_stop',  't_save', 'cfl_target']
             call MPI_BCAST(${VAR}$, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
@@ -234,6 +243,12 @@ contains
             call MPI_BCAST(fluid_pp(i)%Re(1), 2, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
         end do
 
+        do i = 1, num_fluids_max
+            #:for VAR in ['bc_x%alpha_rho_in','bc_x%alpha_in','bc_y%alpha_rho_in','bc_y%alpha_in','bc_z%alpha_rho_in','bc_z%alpha_in']
+                call MPI_BCAST(${VAR}$(i), 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr )
+            #:endfor
+        end do 
+
         do i = 1, num_ibs
             #:for VAR in [ 'radius', 'length_x', 'length_y', &
                 & 'x_centroid', 'y_centroid', 'c', 'm', 'p', 't', 'theta', 'slip' ]
diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp
index 4feb49a42a..509e500d07 100644
--- a/src/simulation/m_start_up.fpp
+++ b/src/simulation/m_start_up.fpp
@@ -1474,6 +1474,9 @@ contains
         !$acc update device(bc_y%vb1, bc_y%vb2, bc_y%vb3, bc_y%ve1, bc_y%ve2, bc_y%ve3)
         !$acc update device(bc_z%vb1, bc_z%vb2, bc_z%vb3, bc_z%ve1, bc_z%ve2, bc_z%ve3)
 
+        !$acc update device(bc_x%grcbc_in, bc_x%grcbc_out, bc_x%grcbc_vel_out)
+        !$acc update device(bc_y%grcbc_in, bc_y%grcbc_out, bc_y%grcbc_vel_out)
+        !$acc update device(bc_z%grcbc_in, bc_z%grcbc_out, bc_z%grcbc_vel_out)
 
         !$acc update device(relax, relax_model)
         if (relax) then
diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py
index 95b2a5bcc2..3f883d51a9 100644
--- a/toolchain/mfc/run/case_dicts.py
+++ b/toolchain/mfc/run/case_dicts.py
@@ -237,15 +237,31 @@ def analytic(self):
         SIMULATION[f'patch_ib({ib_id})%{cmp}_centroid'] = ParamType.REAL
         SIMULATION[f'patch_ib({ib_id})%length_{cmp}'] = ParamType.REAL
 
-for cmp in ["x", "y", "z"]:
-    SIMULATION[f'bc_{cmp}%beg'] = ParamType.INT
-    SIMULATION[f'bc_{cmp}%end'] = ParamType.INT
-    SIMULATION[f'bc_{cmp}%vb1'] = ParamType.REAL
-    SIMULATION[f'bc_{cmp}%vb2'] = ParamType.REAL
-    SIMULATION[f'bc_{cmp}%vb3'] = ParamType.REAL
-    SIMULATION[f'bc_{cmp}%ve1'] = ParamType.REAL
-    SIMULATION[f'bc_{cmp}%ve2'] = ParamType.REAL
-    SIMULATION[f'bc_{cmp}%ve3'] = ParamType.REAL
+    for cmp in ["x", "y", "z"]:
+        SIMULATION[f'bc_{cmp}%beg'] = ParamType.INT
+        SIMULATION[f'bc_{cmp}%end'] = ParamType.INT
+        SIMULATION[f'bc_{cmp}%vb1'] = ParamType.REAL
+        SIMULATION[f'bc_{cmp}%vb2'] = ParamType.REAL
+        SIMULATION[f'bc_{cmp}%vb3'] = ParamType.REAL
+        SIMULATION[f'bc_{cmp}%ve1'] = ParamType.REAL
+        SIMULATION[f'bc_{cmp}%ve2'] = ParamType.REAL
+        SIMULATION[f'bc_{cmp}%ve3'] = ParamType.REAL
+        SIMULATION[f'bc_{cmp}%u_in'] = ParamType.REAL
+        SIMULATION[f'bc_{cmp}%v_in'] = ParamType.REAL
+        SIMULATION[f'bc_{cmp}%w_in'] = ParamType.REAL
+        SIMULATION[f'bc_{cmp}%u_out'] = ParamType.REAL
+        SIMULATION[f'bc_{cmp}%v_out'] = ParamType.REAL
+        SIMULATION[f'bc_{cmp}%w_out'] = ParamType.REAL
+        SIMULATION[f'bc_{cmp}%pres_in'] = ParamType.REAL
+        SIMULATION[f'bc_{cmp}%pres_out'] = ParamType.REAL
+        SIMULATION[f'bc_{cmp}%grcbc_in'] = ParamType.LOG
+        SIMULATION[f'bc_{cmp}%grcbc_out'] = ParamType.LOG
+        SIMULATION[f'bc_{cmp}%grcbc_vel_out'] = ParamType.LOG
+
+    for int_id in range(1, 10+1):
+        for cmp in ["x", "y", "z"]:
+            SIMULATION[f"bc_{cmp}%alpha_rho_in({int_id})"] = ParamType.REAL
+            SIMULATION[f"bc_{cmp}%alpha_in({int_id})"] = ParamType.REAL
 
     for var in ["k", "w", "p", "g"]:
         SIMULATION[f'{var}_{cmp}'] = ParamType.REAL
@@ -295,6 +311,7 @@ def analytic(self):
             SIMULATION[f"integral({int_id})%{cmp}max"] = ParamType.REAL
 
 
+
 # Removed: 'fourier_modes%beg', 'fourier_modes%end', 'chem_wrt'
 # Feel free to return them if they are needed once more.
 POST_PROCESS = COMMON.copy()