diff --git a/tensorproductmultigrid_Source/communication.f90 b/tensorproductmultigrid_Source/communication.f90
index 3f7c44e2f9b2ab950ee59c5343aa05fabddd01fe..0739f83a05e43c53ece21292a7ce1e95dab31c63 100644
--- a/tensorproductmultigrid_Source/communication.f90
+++ b/tensorproductmultigrid_Source/communication.f90
@@ -48,6 +48,7 @@ public::comm_preinitialise
 public::comm_initialise
 public::comm_finalise
 public::scalarprod
+public::boundary_mnh
 public::haloswap_mnh
 public::haloswap
 public::ihaloswap_mnh
@@ -425,6 +426,39 @@ contains
     s = global_sum
   end subroutine scalarprod
 !==================================================================
+! Boundary Neumann
+!==================================================================
+  subroutine boundary_mnh(a)              ! data field
+
+    implicit none
+ 
+    type(scalar3d), intent(inout) :: a
+    !local var
+   
+
+    integer :: n, ix_min,ix_max,iy_min,iy_max
+
+
+    ! Update Real Boundary for Newman case   u(0) = u(1) , etc ...
+    n      = a%grid_param%n
+    ix_min = a%ix_min
+    ix_max = a%ix_max
+    iy_min = a%iy_min
+    iy_max = a%iy_max
+    if ( ix_min == 1 ) then
+       a%s(:,:,0) = a%s(:,:,1)
+    endif
+    if ( ix_max == n ) then
+       a%s(:,:,a%icompx_max+1) = a%s(:,:,a%icompx_max)
+    endif
+    if ( iy_min == 1 ) then
+       a%s(:,0,:) = a%s(:,1,:)
+    endif
+    if ( iy_max == n ) then
+       a%s(:,a%icompy_max+1,:) = a%s(:,a%icompy_max,:)
+    endif
+  end subroutine boundary_mnh
+!==================================================================
 !  Initiate asynchronous halo exchange
 !
 !  For all processes with horizontal indices that are multiples
@@ -666,6 +700,7 @@ contains
     integer, intent(in) :: level
     integer, intent(in) :: m
     type(scalar3d), intent(inout) :: a
+    !local var
     integer :: a_n  ! horizontal grid size
     integer :: nz   ! vertical grid size
     integer, dimension(2) :: p_horiz
@@ -769,6 +804,7 @@ contains
         call finish_timer(t_haloswap(level,m))
       end if
     end if
+
   end subroutine haloswap_mnh
 !==================================================================
   subroutine haloswap(level,m, &  ! multigrid- and processor- level
diff --git a/tensorproductmultigrid_Source/discretisation.f90 b/tensorproductmultigrid_Source/discretisation.f90
index 4962aeb9eed0143dae41ced013951a369b21fd86..92a850e13399ba3af17cfb7b7e303e6f273e55d2 100644
--- a/tensorproductmultigrid_Source/discretisation.f90
+++ b/tensorproductmultigrid_Source/discretisation.f90
@@ -302,7 +302,7 @@ contains
     h = grid_param%L/grid_param%n
     ! Cartesian coefficients
     Tij = h**2
-    xcoef = 0.0
+    xcoef = 0.5_rl ! 0.0
     l2nd = .false. ! .true. ! .false.
     alpha_T(1) = 1.0
     alpha_T(2) = 1.0
@@ -531,6 +531,7 @@ end subroutine construct_vertical_coeff
     integer :: ix,iy,iz
 
     ! r <- A.u
+    call boundary_mnh(u)
     call apply_mnh(u,r)
     ! r <- b - r = b - A.u
     do ix=u%icompx_min,u%icompx_max
@@ -1147,6 +1148,8 @@ end subroutine construct_vertical_coeff
     integer :: iymin(5), iymax(5)
     integer :: iblock, ierr
 
+    call boundary_mnh(u)
+
     ! Set optimal smoothing parameter on each level
     rho = 2.0_rl/(2.0_rl+4.0_rl*model_param%omega2*u%grid_param%n**2/(1.0_rl+4.0_rl*model_param%omega2*u%grid_param%n**2))
 
diff --git a/tensorproductmultigrid_Source/multigrid.f90 b/tensorproductmultigrid_Source/multigrid.f90
index 06f61526ac9bc55a1013795270c3ed6cd6ff31fc..d026a69a7c120b528c81bb8a10bc722951ab62f0 100644
--- a/tensorproductmultigrid_Source/multigrid.f90
+++ b/tensorproductmultigrid_Source/multigrid.f90
@@ -660,7 +660,7 @@ contains
       !local var
       integer :: ix,iy,iz
       real (kind=rl) :: xcoef
-      xcoef = 0.0
+      xcoef = 0.5_rl ! 0.0
       do ix=ixmin(iblock),ixmax(iblock)
         do iy=iymin(iblock),iymax(iblock)
 #ifdef PIECEWISELINEAR