--- a/+scheme/Utux2D.m	Sat Oct 14 22:36:31 2017 -0700
+++ b/+scheme/Utux2D.m	Sat Oct 14 22:38:42 2017 -0700
@@ -25,13 +25,19 @@
         % Other:   'centered'
         coupling_type 
 
+        % String, type of interpolation operators
+        % Default: 'AWW' (Almquist Wang Werpers)
+        % Other:   'MC' (Mattsson Carpenter)
+        interpolation_type
+
         
     end
 
 
     methods 
-         function obj = Utux2D(g ,order, opSet, a, coupling_type)
+         function obj = Utux2D(g ,order, opSet, a, coupling_type, interpolation_type)
             
+            default_arg('interpolation_type','AWW'); 
             default_arg('coupling_type','upwind'); 
             default_arg('a',1/sqrt(2)*[1, 1]); 
             default_arg('opSet',@sbp.D2Standard);
@@ -84,6 +90,7 @@
             obj.order = order;
             obj.a = a;
             obj.coupling_type = coupling_type;
+            obj.interpolation_type = interpolation_type;
             obj.D = -(a(1)*obj.Dx + a(2)*obj.Dy);
 
         end
@@ -117,12 +124,16 @@
              switch neighbour_boundary
                  case {'e','E','east','East'}
                      e_neighbour = neighbour_scheme.e_e;
+                     m_neighbour = neighbour_scheme.m(2);
                  case {'w','W','west','West'}
                      e_neighbour = neighbour_scheme.e_w;
+                     m_neighbour = neighbour_scheme.m(2);
                  case {'n','N','north','North'}
                      e_neighbour = neighbour_scheme.e_n;
+                     m_neighbour = neighbour_scheme.m(1);
                  case {'s','S','south','South'}
                      e_neighbour = neighbour_scheme.e_s;
+                     m_neighbour = neighbour_scheme.m(1);
              end
              
              switch obj.coupling_type
@@ -140,22 +151,74 @@
              otherwise
                 error(['Interface coupling type ' coupling_type ' is not available.'])
              end
+
+             % Check grid ratio for interpolation
+             switch boundary
+                 case {'w','W','west','West','e','E','east','East'}
+                     m = obj.m(2);       
+                 case {'s','S','south','South','n','N','north','North'}
+                     m = obj.m(1);
+             end
+             grid_ratio = m/m_neighbour;
+             if grid_ratio ~= 1
+
+                [ms, index] = sort([m, m_neighbour]);
+                orders = [obj.order, neighbour_scheme.order];
+                orders = orders(index);
+
+                switch obj.interpolation_type
+                case 'MC'
+                    interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2));
+                    if grid_ratio < 1
+                        I_neighbour2local_us = interpOpSet.IF2C;
+                        I_neighbour2local_ds = interpOpSet.IF2C;
+                    elseif grid_ratio > 1
+                        I_neighbour2local_us = interpOpSet.IC2F;
+                        I_neighbour2local_ds = interpOpSet.IC2F;
+                    end
+                case 'AWW'
+                    %String 'C2F' indicates that ICF2 is more accurate.
+                    interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C');
+                    interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); 
+                    if grid_ratio < 1 
+                        % Local is coarser than neighbour
+                        I_neighbour2local_us = interpOpSetC2F.IF2C;
+                        I_neighbour2local_ds = interpOpSetF2C.IF2C;
+                    elseif grid_ratio > 1
+                        % Local is finer than neighbour 
+                        I_neighbour2local_us = interpOpSetF2C.IC2F;
+                        I_neighbour2local_ds = interpOpSetC2F.IC2F;
+                    end
+                otherwise
+                    error(['Interpolation type ' obj.interpolation_type ...
+                         ' is not available.' ]);
+                end
+
+             else 
+                % No interpolation required
+                I_neighbour2local_us = speye(m,m);
+                I_neighbour2local_ds = speye(m,m);
+            end    
              
              switch boundary
                  case {'w','W','west','West'}
                      tau = sigma_ds*obj.a(1)*obj.e_w*obj.H_y;
-                     closure = obj.Hi*tau*obj.e_w';       
+                     closure = obj.Hi*tau*obj.e_w';
+                     penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';     
                  case {'e','E','east','East'}
                      tau = sigma_us*obj.a(1)*obj.e_e*obj.H_y;
                      closure = obj.Hi*tau*obj.e_e';
+                     penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
                  case {'s','S','south','South'}
                      tau = sigma_ds*obj.a(2)*obj.e_s*obj.H_x;
                      closure = obj.Hi*tau*obj.e_s'; 
+                     penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
                  case {'n','N','north','North'}
                      tau = sigma_us*obj.a(2)*obj.e_n*obj.H_x;
                      closure = obj.Hi*tau*obj.e_n';
+                     penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
              end
-             penalty = -obj.Hi*tau*e_neighbour';
+             
                  
          end