diff --git a/src/ZSOLVER/tensorproductmultigrid_Source/conjugategradient.f90 b/src/ZSOLVER/tensorproductmultigrid_Source/conjugategradient.f90
index 30277d8625ad6735444ad7335b761ab6ed369b7a..34eacb6b86a7cf47e41a6b14a08ea2685e8c7206 100644
--- a/src/ZSOLVER/tensorproductmultigrid_Source/conjugategradient.f90
+++ b/src/ZSOLVER/tensorproductmultigrid_Source/conjugategradient.f90
@@ -206,10 +206,13 @@ contains
     allocate(Ap%st(1-halo_size:nlocal+halo_size,   &
                    1-halo_size:nlocal+halo_size,   &
                    0:nz+1) )
+    !$acc enter data create (r%st,z%st,p%st,Ap%st)
+    !$acc kernels
     r%st = 0.0_rl
     z%st = 0.0_rl
     p%st = 0.0_rl
     Ap%st = 0.0_rl
+    !$acc end kernels
     endif
 
 
@@ -307,6 +310,7 @@ contains
     end if
 
     if (LUseT) then
+    !$acc exit data delete(r%st,z%st,p%st,Ap%st)
     deallocate(r%st)
     deallocate(z%st)
     deallocate(p%st)
@@ -419,10 +423,13 @@ contains
     allocate(Ap%st(1-halo_size:nlocal+halo_size,  &
                    1-halo_size:nlocal+halo_size,  &
                    0:nz+1) )
+    !$acc enter data create (r%st,z%st,p%st,Ap%st)
+    !$acc kernels
     r%st = 0.0_rl
     z%st = 0.0_rl
     p%st = 0.0_rl
     Ap%st = 0.0_rl
+    !$acc end kernels
     end if
 
     ! Initialise
@@ -518,6 +525,7 @@ contains
     deallocate(Ap%s)
     end if
     if (LUseT) then
+    !$acc exit data delete(r%st,z%st,p%st,Ap%st)
     deallocate(r%st)
     deallocate(z%st)
     deallocate(p%st)
diff --git a/src/ZSOLVER/tensorproductmultigrid_Source/datatypes.f90 b/src/ZSOLVER/tensorproductmultigrid_Source/datatypes.f90
index d0ddee6daae5fa031349db044389cb510c92b144..e4c2a370395eca6f76468299e4b9d01213ab4862 100644
--- a/src/ZSOLVER/tensorproductmultigrid_Source/datatypes.f90
+++ b/src/ZSOLVER/tensorproductmultigrid_Source/datatypes.f90
@@ -195,6 +195,7 @@ private
     allocate(zphi_st(1-halo_size:nlocal+halo_size, &
                     1-halo_size:nlocal+halo_size, &
                     0:grid_param%nz+1))
+    !$acc enter data create (zphi_st)
     phi%st => zphi_st
     !$acc kernels
     zphi_st(:,:,:) = 0.0_rl
@@ -211,7 +212,10 @@ private
     type(scalar3d), intent(inout) :: phi
 
     if (LUseO) deallocate(phi%s)
-    if (LUseT) deallocate(phi%st)
+    if (LUseT) then
+       !$acc exit data delete(phi%st)
+       deallocate(phi%st)
+    end if
 
   end subroutine destroy_scalar3d
 
diff --git a/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90 b/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90
index 11a793ce2fcf4bf89123a2a837f4c3d055868604..1bd6d51d85bcd52b6a93103136e06caf8881a182 100644
--- a/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90
+++ b/src/ZSOLVER/tensorproductmultigrid_Source/discretisation.f90
@@ -204,6 +204,7 @@ contains
     allocate(vert_coeff%c(grid_param%nz))
     allocate(vert_coeff%d(grid_param%nz))
     call construct_vertical_coeff_mnh(PA_K,PB_K,PC_K,PD_K)
+    !$acc enter data copyin(vert_coeff%a,vert_coeff%b,vert_coeff%c,vert_coeff%d)
   end subroutine discretisation_initialise_mnh
 
   subroutine discretisation_initialise(grid_param_in, &
@@ -817,6 +818,7 @@ end subroutine construct_vertical_coeff
        allocate(r%st(1-halo_size:nlocal+halo_size, &
                      1-halo_size:nlocal+halo_size, &
                      0:nz+1))
+       !$acc enter data create (r%st)
     end if
 
     call calculate_residual(level,m,b,u,r)
@@ -843,7 +845,10 @@ end subroutine construct_vertical_coeff
                             + (log_res_final - log_res_initial)
     nsmooth_total(level) = nsmooth_total(level) + nsmooth
     if (LUseO) deallocate(r%s)
-    if (LUseT) deallocate(r%st)
+    if (LUseT) then
+       !$acc exit data delete(r%st)
+       deallocate(r%st)       
+    end if
 #endif
   end subroutine smooth_mnh
 !==================================================================
@@ -888,6 +893,7 @@ end subroutine construct_vertical_coeff
     allocate(r%st(1-halo_size:nlocal+halo_size,   &
                   1-halo_size:nlocal+halo_size,   &
                   0:nz+1) )
+    !$acc enter data create (r%st)
     endif
     log_res_initial = log(l2norm(r))
 #endif
@@ -912,7 +918,10 @@ end subroutine construct_vertical_coeff
                             + (log_res_final - log_res_initial)
     nsmooth_total(level) = nsmooth_total(level) + nsmooth
     if (LUseO) deallocate(r%s)
-    if (LUseT) deallocate(r%st)
+    if (LUseT) then
+       !$acc exit data delete(r%st)
+       deallocate(r%st)
+    end if
 #endif
   end subroutine smooth
 !==================================================================
@@ -1929,6 +1938,7 @@ end subroutine line_Jacobi_mnh
          !call construct_alpha_T_cst_mnh(b%grid_param,alpha_T,Tij)
          Tij = ( b%grid_param%L/b%grid_param%n ) ** 2
          alpha_div_Tij = 4.0_rl / Tij
+         !$acc enter data copyin(Tij,alpha_div_Tij)
          !print*,"level=",level," Tij=",Tij," alpha_div_Tij=",alpha_div_Tij
 
          allocate(Ttridiag_mnh(level)%tmp_k(size(zb_k)))
diff --git a/src/ZSOLVER/tensorproductmultigrid_Source/multigrid.f90 b/src/ZSOLVER/tensorproductmultigrid_Source/multigrid.f90
index d1eefcc42910aa2eabacfcfac27d46a5e4a5a240..e2a0b94a9f31087bb2ffe8b1173fd72963bab935 100644
--- a/src/ZSOLVER/tensorproductmultigrid_Source/multigrid.f90
+++ b/src/ZSOLVER/tensorproductmultigrid_Source/multigrid.f90
@@ -275,18 +275,20 @@ contains
          allocate(zxu_mg_st(1-halo_size:nlocal+halo_size, &
               1-halo_size:nlocal+halo_size, &
               0:nz+1))
+         !$acc enter data create (zxu_mg_st)
          xu_mg(level,m)%st => zxu_mg_st
          
          allocate(zxb_mg_st(1-halo_size:nlocal+halo_size, &
               1-halo_size:nlocal+halo_size, &
               0:nz+1))
+         !$acc enter data create (zxb_mg_st)
          xb_mg(level,m)%st => zxb_mg_st
          
          allocate(zxr_mg_st(1-halo_size:nlocal+halo_size, &
               1-halo_size:nlocal+halo_size, &
               0:nz+1))
-          xr_mg(level,m)%st => zxr_mg_st
-
+         !$acc enter data create (zxr_mg_st)
+         xr_mg(level,m)%st => zxr_mg_st
       
       !$acc kernels
       zxu_mg_st(:,:,:) = 0.0_rl
@@ -429,9 +431,14 @@ contains
       endif
 
       if (LUseT) then
-      deallocate(xu_mg(level,m)%st)
-      deallocate(xb_mg(level,m)%st)
-      deallocate(xr_mg(level,m)%st)
+         !$acc exit data delete(xu_mg(level,m)%st)  
+         deallocate(xu_mg(level,m)%st)
+         !
+         !$acc exit data delete(xb_mg(level,m)%st)
+         deallocate(xb_mg(level,m)%st)
+         !
+         !$acc exit data delete(xr_mg(level,m)%st)
+         deallocate(xr_mg(level,m)%st)
       endif
 
       ! If we are below L_split, split data