diff --git a/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90 b/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90
index 86d0b0a641f42b3be46723ac3295dcd60263fcff..8029d6c51375b489d5a6dc771df88ce45347d54d 100644
--- a/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90
+++ b/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90
@@ -663,6 +663,8 @@ end subroutine construct_vertical_coeff
 
     real , dimension(:,:,:) , pointer , contiguous :: zr_st , zb_st
 
+    integer :: nib,nie,njb,nje,nzb,nze
+
     ! r <- A.u
     !call boundary_mnh(u)
     call apply_mnh(u,r)
@@ -694,11 +696,29 @@ end subroutine construct_vertical_coeff
 
        zr_st => r%st
        zb_st => b%st
-       !$acc kernels present(zr_st,zb_st)
-       zr_st(iib:iie,ijb:ije,ikb:ike) = zb_st(iib:iie,ijb:ije,ikb:ike) - zr_st(iib:iie,ijb:ije,ikb:ike)
-       !$acc end kernels
+
+       nib = Lbound(zr_st,1) ; nie = Ubound(zr_st,1)
+       njb = Lbound(zr_st,2) ; nje = Ubound(zr_st,2)
+       nzb = Lbound(zr_st,3) ; nze = Ubound(zr_st,3)
+
+       call calculate_residual_mnh_dim(zr_st,zb_st)
+
     endif
 
+  contains
+
+    subroutine calculate_residual_mnh_dim(pzr_st,pzb_st)
+      implicit none
+
+      real :: pzr_st(nib:nie,njb:nje,nzb:nze), &
+              pzb_st(nib:nie,njb:nje,nzb:nze)
+
+      !$acc kernels present(pzr_st,pzb_st)
+        pzr_st(iib:iie,ijb:ije,ikb:ike) = pzb_st(iib:iie,ijb:ije,ikb:ike) - pzr_st(iib:iie,ijb:ije,ikb:ike)
+      !$acc end kernels
+      
+    end subroutine calculate_residual_mnh_dim
+
   end subroutine calculate_residual_mnh
 
   subroutine calculate_residual(level,m,b,u,r)