diff --git a/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90 b/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90
index 75a6c3ad22d8ca3a8d2efb1f7da4661c4c59cbed..86d0b0a641f42b3be46723ac3295dcd60263fcff 100644
--- a/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90
+++ b/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90
@@ -744,6 +744,8 @@ end subroutine construct_vertical_coeff
     integer :: ii,ij
     integer :: ize
 
+    integer :: nib,nie,njb,nje,nzb,nze
+
     call boundary_mnh(u)
 
     if (LUseO) then     
@@ -792,37 +794,53 @@ end subroutine construct_vertical_coeff
        zb_k => vert_coeff%b
        zc_k => vert_coeff%c
        zd_k => vert_coeff%d
+
+       nib = Lbound(zv_st,1) ; nie = Ubound(zv_st,1)
+       njb = Lbound(zv_st,2) ; nje = Ubound(zv_st,2)
+       nzb = Lbound(zv_st,3) ; nze = Ubound(zv_st,3)
+
+       call apply_mnh_dim(zv_st,zu_st,zb_k,zc_k,zd_k)
     
-       !$acc kernels present_cr(zb_k,zv_st)
-       iz=1
-       !$mnh_do_concurrent( ii=iib:iie , ij=ijb:ije )
-          zv_st(ii,ij,iz)   = zd_k(iz)* ( (-zb_k(iz)-zc_k(iz))*Tij * zu_st(ii,ij,iz  )  &
-               +zb_k(iz)          *Tij * zu_st(ii,ij,iz+1)  )
-       !$mnh_end_do()
-       
-       !
-!!$       !$acc loop seq
-!!$       do iz=2,ize-1       
-       !$mnh_do_concurrent( ii=iib:iie , ij=ijb:ije , iz=2:ize-1 )
-                zv_st(ii,ij,iz) = zd_k(iz)* ( ((-zb_k(iz)-zc_k(iz))*Tij - 4.0_rl ) * zu_st(ii,ij,iz)    &
-                                                +zb_k(iz)          *Tij            * zu_st(ii,ij,iz+1)  &
-                                                         +zc_k(iz) *Tij            * zu_st(ii,ij,iz-1)  &
-                                           +                                         zu_st(ii+1,ij,iz) &
-                                           +                                         zu_st(ii-1,ij,iz) &
-                                           +                                         zu_st(ii,ij+1,iz) &
-                                           +                                         zu_st(ii,ij-1,iz) &
+    endif
+
+  contains
+
+    subroutine apply_mnh_dim(pzv_st,pzu_st,pzb_k,pzc_k,pzd_k)
+
+      implicit none
+
+      real :: pzv_st(nib:nie,njb:nje,nzb:nze), &
+              pzu_st(nib:nie,njb:nje,nzb:nze), &
+              pzb_k(ize),pzc_k(ize),pzd_k(ize)
+
+      !$acc kernels present_cr(pzb_k,pzv_st,pzu_st)
+      iz=1
+      !$mnh_do_concurrent( ii=iib:iie , ij=ijb:ije )
+         pzv_st(ii,ij,iz) = pzd_k(iz)* ( (-pzb_k(iz)-pzc_k(iz))*Tij * pzu_st(ii,ij,iz  )  &
+                                          +pzb_k(iz)           *Tij * pzu_st(ii,ij,iz+1)  )
+      !$mnh_end_do()
+      
+      !
+      !$mnh_do_concurrent( ii=iib:iie , ij=ijb:ije , iz=2:ize-1 )
+         pzv_st(ii,ij,iz) = pzd_k(iz)* ( ((-pzb_k(iz)-pzc_k(iz))*Tij - 4.0_rl ) * pzu_st(ii,ij,iz)   &
+                                           +pzb_k(iz)           *Tij            * pzu_st(ii,ij,iz+1) &
+                                                     +pzc_k(iz) *Tij            * pzu_st(ii,ij,iz-1) &
+                                     +                                            pzu_st(ii+1,ij,iz) &
+                                     +                                            pzu_st(ii-1,ij,iz) &
+                                     +                                            pzu_st(ii,ij+1,iz) &
+                                     +                                            pzu_st(ii,ij-1,iz) &
                                        )
        !$mnh_end_do()
-!!$          end do
        !
        iz=ize
        !$mnh_do_concurrent( ii=iib:iie , ij=ijb:ije )       
-             zv_st(ii,ij,iz) = zd_k(iz)*  (  (-zb_k(iz)-zc_k(iz))*Tij  * zu_st(ii,ij,iz)    &
-                                                       +zc_k(iz) *Tij  * zu_st(ii,ij,iz-1)  )
+             pzv_st(ii,ij,iz) = pzd_k(iz)*  (  (-pzb_k(iz)-pzc_k(iz))*Tij  * pzu_st(ii,ij,iz)    &
+                                                          +pzc_k(iz) *Tij  * pzu_st(ii,ij,iz-1)  )
        !$mnh_end_do()
-       !$acc end kernels
-    endif
+       !$acc end kernels     
     
+     end subroutine apply_mnh_dim
+
   end subroutine apply_mnh
 
   subroutine apply(u,v)