diff --git a/src/BIO/BIO_mem.f90 b/src/BIO/BIO_mem.f90
index fb5e3bd2..b3544716 100644
--- a/src/BIO/BIO_mem.f90
+++ b/src/BIO/BIO_mem.f90
@@ -18,6 +18,10 @@ MODULE BIO_mem
       double precision, allocatable :: ogstm_ph(:,:,:) ! GUESS for FOLLOWS algorithm
       double precision, allocatable :: NPPF2(:,:,:)
       double precision, allocatable :: ogstm_co2(:,:), co2_IO(:,:,:)
+      double precision, allocatable :: sediPPY(:,:)
+      double precision, allocatable :: local_D3DIAGNOS(:,:)
+      double precision, allocatable :: local_D2DIAGNOS(:,:)
+      double precision, allocatable :: er(:,:)
       double precision:: ice
 
 
@@ -38,7 +42,7 @@ subroutine myalloc_BIO()
        allocate(co2_IO(jpj,jpi,2))    
         
        co2_IO    = huge(co2_IO(1,1,1))
-       allocate(ogstm_sedipi(jpk,jpj,jpi,4)) 
+       allocate(ogstm_sedipi(jpk,jpj,jpi,4))
        ogstm_sedipi    = huge(ogstm_sedipi(1,1,1,1))
        allocate(ogstm_ph(jpk,jpj,jpi))       
        ogstm_ph        = huge(ogstm_ph(1,1,1))
@@ -48,6 +52,12 @@ subroutine myalloc_BIO()
                  ! and used in hard_tissue_pump.F also in land points
        ice=0
 
+       allocate(sediPPY(jpi * jpj * jpk, 4))
+       allocate(local_D3DIAGNOS(jpi * jpj * jpk, jptra_dia))
+       allocate(local_D2DIAGNOS(jpi * jpj, jptra_dia_2d))
+       allocate(er(jpi * jpj * jpk, 11))
+       !$acc enter data create(ogstm_co2,ogstm_sedipi,ogstm_ph,sediPPY,local_D3DIAGNOS,local_D2DIAGNOS,er)
+
 #ifdef Mem_Monitor
        mem_all=get_mem(err) - aux_mem
 #endif
@@ -63,6 +73,11 @@ subroutine clean_memory_bio()
             deallocate(ogstm_sedipi)
             deallocate(ogstm_ph)
             deallocate(NPPF2)
+            deallocate(sediPPY)
+            deallocate(local_D3DIAGNOS)
+            deallocate(local_D2DIAGNOS)
+            deallocate(er)
+            !$acc exit data delete(ogstm_co2,ogstm_sedipi,ogstm_ph,sediPPY,local_D3DIAGNOS,local_D2DIAGNOS,er)
 
       end subroutine clean_memory_bio
 
diff --git a/src/BIO/OPT_mem.f90 b/src/BIO/OPT_mem.f90
index 60af75dc..a08e2228 100644
--- a/src/BIO/OPT_mem.f90
+++ b/src/BIO/OPT_mem.f90
@@ -18,8 +18,9 @@ MODULE OPT_mem
 
 
       INTEGER, allocatable :: itabe(:),imaske(:,:) 
-      double precision, allocatable :: zpar(:,:),xEPS_ogstm(:,:)
-      double precision, allocatable :: zpar0m(:),zpar100(:) 
+      ! double precision, allocatable :: zpar(:,:)
+      double precision, allocatable :: xEPS_ogstm(:,:)
+      ! double precision, allocatable :: zpar0m(:),zpar100(:)
       double precision, allocatable :: kef(:,:)
       double precision, allocatable :: kextIO(:,:,:)
       real, allocatable :: zkef_f (:,:)
@@ -40,17 +41,21 @@ subroutine myalloc_OPT()
        allocate(imaske(jpk,jpi))   
        imaske  = huge(imaske(1,1))
 !!!$omp parallel default (none) shared(jpk,jpi)
-       allocate(zpar(jpk,jpi))     
-       zpar    = huge(zpar(1,1))
-       allocate(xEPS_ogstm(jpk,jpi))     
+       ! allocate(zpar(jpk,jpi))
+       ! zpar    = huge(zpar(1,1))
+       allocate(xEPS_ogstm(jpk,jpi))
+       !$acc enter data create(xEPS_ogstm)
+       !$acc kernels default(present)
        xEPS_ogstm    = huge(xEPS_ogstm(1,1))
-       allocate(zpar0m(jpi))        
-       zpar0m  = huge(zpar0m(1))
-       allocate(zpar100(jpi))       
-       zpar100 = huge(zpar100(1))
+       !$acc end kernels
+       ! allocate(zpar0m(jpi))
+       ! zpar0m  = huge(zpar0m(1))
+       ! allocate(zpar100(jpi))
+       ! zpar100 = huge(zpar100(1))
 !!!$omp end parallel
 
-       allocate(kef(jpj,jpi))       
+       allocate(kef(jpj,jpi))
+       !$acc enter data create(kef)
        kef     = huge(kef(1,1))
        allocate(kextIO(jpj,jpi,2))  
        kextIO  = huge(kextIO(1,1,1))
@@ -71,10 +76,12 @@ subroutine clean_memory_opt
 
           deallocate(itabe)
           deallocate(imaske)
-          deallocate(zpar)
+          ! deallocate(zpar)
+          !$acc exit data delete(xEPS_ogstm)
           deallocate(xEPS_ogstm)
-          deallocate(zpar0m)
-          deallocate(zpar100)
+          ! deallocate(zpar0m)
+          ! deallocate(zpar100)
+          !$acc exit data delete(kef)
           deallocate(kef)
           deallocate(kextIO)
 
diff --git a/src/BIO/SED_mem.f90 b/src/BIO/SED_mem.f90
index 0062dce6..4dae9d15 100644
--- a/src/BIO/SED_mem.f90
+++ b/src/BIO/SED_mem.f90
@@ -46,7 +46,7 @@ subroutine myalloc_SED()
 #endif
        dimen_jvsed=0
 
-       allocate(sed_idx(nsed))  
+       allocate(sed_idx(nsed))
        sed_idx = huge(sed_idx(1))
 
        sed_idx(1)  = ppR6c
@@ -79,10 +79,13 @@ subroutine myalloc_SED()
        jarr_sed     = huge(jarr_sed(1,1))
        allocate(jarr_sed_flx(jpk,jpi*jpj)) 
        jarr_sed_flx = huge(jarr_sed_flx(1,1)) 
-       allocate( ztra(nsed,ntids))         
+#ifndef _OPENACC
+       allocate( ztra(nsed,ntids))
        ztra         = huge(ztra(1,1))
        allocate(zwork(jpk,nsed, ntids))   
        zwork        = huge(zwork(1,1,1))
+#endif
+       !$acc enter data create(sed_idx,jarr_sed,jarr_sed_flx)
 
 
 #ifdef Mem_Monitor
@@ -90,11 +93,26 @@ subroutine myalloc_SED()
 #endif
 
       END subroutine myalloc_SED
+
+#ifdef _OPENACC
+      subroutine myalloc_SED_gpu()
+
+        allocate(ztra(nsed,dimen_jvsed))
+        allocate(zwork(jpk,nsed,dimen_jvsed))
+        !$acc enter data create(ztra,zwork)
+        !$acc kernels default(present)
+        ztra = huge(ztra(1,1))
+        zwork = huge(zwork(1,1,1))
+        !$acc end kernels
+
+      end subroutine myalloc_SED_gpu
+#endif
       
       
       
       subroutine clean_memory_sed
 
+          !$acc exit data delete(ztra,zwork,sed_idx,jarr_sed,jarr_sed_flx)
           deallocate(sed_idx)
           deallocate(jarr_sed)
           deallocate(jarr_sed_flx)
diff --git a/src/BIO/trcbio.f90 b/src/BIO/trcbio.f90
index c6e52bcc..93c0939b 100644
--- a/src/BIO/trcbio.f90
+++ b/src/BIO/trcbio.f90
@@ -48,12 +48,7 @@ SUBROUTINE trcbio
 
   IMPLICIT NONE
 
-  double precision, dimension(jpi * jpj * jpk, 4) :: sediPPY
-  double precision, dimension(jpi * jpj * jpk, jptra_dia) :: local_D3DIAGNOS
-  double precision, dimension(jpi * jpj, jptra_dia_2d) :: local_D2DIAGNOS
-  double precision, dimension(jpi * jpj * jpk, 11) :: er
-
-  integer :: jk, jj, ji, jn, jlinear2d, jlinear3d, bottom
+  integer :: jk, jj, ji, jn, jlinear2d, jlinear3d, bottom, queue
   double precision :: correct_fact, gdept_local, gdeptmax_local
   
   integer :: year, month, day
@@ -62,6 +57,10 @@ SUBROUTINE trcbio
   BIOparttime = MPI_WTIME()
   
   ! Initialization
+
+  queue=1
+
+  !$acc kernels default(present) async(queue)
   D3STATE = 1.0
   er = 1.0
   er(:,10) = 8.1
@@ -70,120 +69,141 @@ SUBROUTINE trcbio
 
   ! ogstm_sediPI appear to be unused
   ogstm_sediPI = 0.
+  !$acc end kernels
+
+  ! NOTE: this kernel *needs* to be executed synchronously as we need the
+  ! reduced `bottom` scalar value on host before launching the next kernel.
+  bottom=0
+  !$acc parallel loop gang vector reduction(max:bottom) collapse(2) default(present)
+  do ji = 1, jpi
+     do jj = 1, jpj
+        bottom = max(bottom,mbathy(jj,ji))
+     end do
+  end do
+  !$acc end parallel loop
 
   ! Set D3STATE (pass state to BFM)
+  !$acc parallel loop gang vector default(present) collapse(4) async(queue)
   do jn = 1, jptra
     do ji = 1, jpi
       do jj = 1, jpj
-        if (bfmmask(1, jj, ji) == 0) then
-          cycle
-        else
-          bottom = mbathy(jj, ji)
-          do jk = 1, bottom
-            jlinear3d = jk + (jj - 1) * jpk + (ji - 1) * jpk * jpj
-            D3STATE(jlinear3d, jn) = trn(jk, jj, ji, jn)
+         do jk = 1, bottom
+            if (.not. bfmmask(1, jj, ji) == 0) then
+               if (jk <= mbathy(jj, ji)) then
+                  jlinear3d = jk + (jj - 1) * jpk + (ji - 1) * jpk * jpj
+                  D3STATE(jlinear3d, jn) = trn(jk, jj, ji, jn)
+               endif
+            endif
           end do
-        end if
       end do
     end do
   end do
+  !$acc end parallel loop
 
   call read_date_string(COMMON_DATEstring, year, month, day, sec)
 
-  ! Set er 
+  ! Set er
+  !$acc parallel loop gang vector default(present) collapse(3) async(queue)
   do ji = 1, jpi
-    do jj = 1, jpj
-      if (bfmmask(1, jj, ji) == 0) then
-        cycle
-      else
-        bottom = mbathy(jj, ji)
+     do jj = 1, jpj
         do jk = 1, bottom
-          jlinear3d = jk + (jj - 1) * jpk + (ji - 1) * jpk * jpj
-          if (jk .eq. 1) then
-            er(jlinear3d, 4) = ice ! from 0 to 1 adimensional
-            er(jlinear3d, 5) = ogstm_co2(jj, ji) ! CO2 Mixing Ratios (ppm)  390
-            er(jlinear3d, 7) = DAY_LENGTH(jj, ji) ! fotoperiod expressed in hours
-            er(jlinear3d, 9) = vatm(jj, ji) ! wind speed (m/s)
-          end if
-          er(jlinear3d, 1) = tn(jk, jj, ji) ! Temperature (Celsius)
-          er(jlinear3d, 2) = sn(jk, jj, ji) ! Salinity PSU
-          er(jlinear3d, 3) = rho(jk, jj, ji) ! Density Kg/m3
-          er(jlinear3d, 6) = instant_par_from_sec(sec, xpar(jk, jj, ji)) ! PAR umoles/m2/s | Watt to umoles photons W2E=1./0.217
-          er(jlinear3d, 8) = e3t(jk, jj, ji) ! depth in meters of the given cell
-          er(jlinear3d, 10) = ogstm_PH(jk, jj, ji) ! 8.1
+           if (.not. bfmmask(1, jj, ji) == 0) then
+              if (jk <= mbathy(jj, ji)) then
+                 jlinear3d = jk + (jj - 1) * jpk + (ji - 1) * jpk * jpj
+                 if (jk .eq. 1) then
+                    er(jlinear3d, 4) = ice ! from 0 to 1 adimensional
+                    er(jlinear3d, 5) = ogstm_co2(jj, ji) ! CO2 Mixing Ratios (ppm)  390
+                    er(jlinear3d, 7) = DAY_LENGTH(jj, ji) ! fotoperiod expressed in hours
+                    er(jlinear3d, 9) = vatm(jj, ji) ! wind speed (m/s)
+                 end if
+                 er(jlinear3d, 1) = tn(jk, jj, ji) ! Temperature (Celsius)
+                 er(jlinear3d, 2) = sn(jk, jj, ji) ! Salinity PSU
+                 er(jlinear3d, 3) = rho(jk, jj, ji) ! Density Kg/m3
+                 er(jlinear3d, 6) = instant_par_from_sec(sec, xpar(jk, jj, ji)) ! PAR umoles/m2/s | Watt to umoles photons W2E=1./0.217
+                 er(jlinear3d, 8) = e3t(jk, jj, ji) ! depth in meters of the given cell
+                 er(jlinear3d, 10) = ogstm_PH(jk, jj, ji) ! 8.1
 #ifdef gdept1d
-          gdept_local = gdept(jk)
-          gdeptmax_local = gdept(jpk)
+                 gdept_local = gdept(jk)
+                 gdeptmax_local = gdept(jpk)
 #else
-          gdept_local = gdept(jk, jj, ji)
-          gdeptmax_local = gdept(jpk, jj, ji)
+                 gdept_local = gdept(jk, jj, ji)
+                 gdeptmax_local = gdept(jpk, jj, ji)
 #endif
-          if (gdept_local .lt. 1000.0D0) then
-            correct_fact = 1.0D0
-          else if (gdept_local .lt. 2000.0D0) then
-            correct_fact = 0.25D0
-          else
-            correct_fact = 0.0D0
-          end if
-          er(jlinear3d, 11) = correct_fact * (gdeptmax_local - gdept_local) / gdept_local
+                 if (gdept_local .lt. 1000.0D0) then
+                    correct_fact = 1.0D0
+                 else if (gdept_local .lt. 2000.0D0) then
+                    correct_fact = 0.25D0
+                 else
+                    correct_fact = 0.0D0
+                 end if
+                 er(jlinear3d, 11) = correct_fact * (gdeptmax_local - gdept_local) / gdept_local
+              endif
+           end if
         end do
-      end if
-    end do
+     end do
   end do
-   
+  !$acc end parallel loop
+
+  !$acc wait(queue)
+
   call BFM1D_Input_EcologyDynamics(mbathy, er) ! here mbathy was bottom
   call BFM1D_reset()
   call EcologyDynamics()
   call BFM1D_Output_EcologyDynamics(sediPPY, local_D3DIAGNOS, local_D2DIAGNOS)
 
-  ! The following copies could be avoided
+  !$acc parallel loop gang vector collapse(4) default(present) async(queue)
   do jn = 1, max(4, jptra, jptra_dia)
     do ji = 1, jpi
       do jj = 1, jpj
-        if (bfmmask(1, jj, ji) == 0) then
-          cycle
-        else
-          bottom = mbathy(jj, ji)
-          do jk = 1, bottom
-            jlinear3d = jk + (jj - 1) * jpk + (ji - 1) * jpk * jpj
-            if (jn .le. jptra) then
-             tra(jk, jj, ji, jn) = tra(jk, jj, ji, jn) + D3SOURCE(jlinear3d, jn) ! trend
-            end if
-            if (jn .le. jptra_dia) then 
-              tra_DIA(jk, jj, ji, jn) = local_D3DIAGNOS(jlinear3d, jn)
-            end if
-            if (jn .le. 4) then
-              ogstm_sediPI(jk, jj, ji, jn) = sediPPY(jlinear3d, jn)   ! BFM output of sedimentation speed (m/d)
+        do jk = 1, bottom
+          if (.not. bfmmask(1, jj, ji) == 0) then
+            if (jk <= mbathy(jj, ji)) then
+               jlinear3d = jk + (jj - 1) * jpk + (ji - 1) * jpk * jpj
+               if (jn .le. jptra) then
+                  tra(jk, jj, ji, jn) = tra(jk, jj, ji, jn) + D3SOURCE(jlinear3d, jn) ! trend
+               endif
+               if (jn .le. jptra_dia) then
+                  tra_DIA(jk, jj, ji, jn) = local_D3DIAGNOS(jlinear3d, jn)
+               endif
+               if (jn .le. 4) then
+                  ogstm_sediPI(jk, jj, ji, jn) = sediPPY(jlinear3d, jn)   ! BFM output of sedimentation speed (m/d)
+               endif
             end if
-          end do
-        end if
+          end if
+        end do
       end do
     end do
   end do
-  
+  !$acc end parallel loop
+
+  !$acc parallel loop gang vector collapse(3) default(present) async(queue)
   do ji = 1, jpi
     do jj = 1, jpj
-      if (bfmmask(1, jj, ji) == 0) then
-        cycle
-      else
-        bottom = mbathy(jj, ji)
-        jlinear2d = jj + (ji - 1) * jpj
-        tra_DIA_2d(:, jj, ji) = local_D2DIAGNOS(jlinear2d, :)
-        do jk = 1, bottom
-          jlinear3d = jk + (jj - 1) * jpk + (ji - 1) * jpk * jpj
-          ogstm_PH(jk, jj, ji) = local_D3DIAGNOS(jlinear3d, pppH) ! Follows solver guess, put 8.0 if pppH is not defined
-        end do
-      end if
+      do jk = 1, bottom
+         if (.not. bfmmask(1, jj, ji) == 0) then
+            if (jk <= mbathy(jj, ji)) then
+               jlinear2d = jj + (ji - 1) * jpj
+               tra_DIA_2d(:, jj, ji) = local_D2DIAGNOS(jlinear2d, :)
+               jlinear3d = jk + (jj - 1) * jpk + (ji - 1) * jpk * jpj
+               ogstm_PH(jk, jj, ji) = local_D3DIAGNOS(jlinear3d, pppH) ! Follows solver guess, put 8.0 if pppH is not defined
+            end if
+         end if
+      end do
     end do
   end do
+  !$acc end parallel loop
 
   !---------------------------------------------------------------------
   ! BEGIN BC_REFACTORING SECTION
   !---------------------------------------------------------------------
+  ! XXX: when should we care about this ?
   call boundaries%fix_diagnostic_vars()
   !----------------------------------------------------------------------
   ! END BC_REFACTORING SECTION
   !---------------------------------------------------------------------
+
+  !$acc wait(queue)
+
   BIOparttime = MPI_WTIME() - BIOparttime
   BIOtottime = BIOtottime + BIOparttime
 END SUBROUTINE trcbio
diff --git a/src/BIO/trclec.f90 b/src/BIO/trclec.f90
index d6af01f4..b99be8a3 100644
--- a/src/BIO/trclec.f90
+++ b/src/BIO/trclec.f90
@@ -141,6 +141,7 @@ SUBROUTINE trclec
           if (ctr_hf(ji).eq.1) jptra_high = jptra_high + 1
       enddo
       allocate(highfreq_table(jptra_HIGH))
+      !$acc enter data create(highfreq_table)
        highfreq_table = huge(highfreq_table(1))
 
       jptra_high = 0
@@ -165,6 +166,7 @@ SUBROUTINE trclec
 
       if (lwp) write(*,*) 'High freq diagnostics number :', jptra_dia_HIGH
       allocate(highfreq_table_dia(jptra_dia_HIGH))
+      !$acc enter data create(highfreq_table_dia)
 
       jptra_dia_high = 0
 
@@ -188,6 +190,7 @@ SUBROUTINE trclec
 
       if (lwp) write(*,*) 'High freq diagnostics number 2d:', jptra_dia2d_HIGH
       allocate(highfreq_table_dia2d(jptra_dia2d_HIGH))
+      !$acc enter data create(highfreq_table_dia2d)
 
       jptra_dia2d_high = 0
 
diff --git a/src/BIO/trcopt.f90 b/src/BIO/trcopt.f90
index b0f21d82..e76308f6 100644
--- a/src/BIO/trcopt.f90
+++ b/src/BIO/trcopt.f90
@@ -18,7 +18,7 @@ SUBROUTINE trcopt
       double precision :: conversion
 #if defined key_trc_nnpzddom || defined key_trc_npzd || key_trc_bfm
 
-      INTEGER :: jk,jj,ji
+      INTEGER :: jk,jj,ji,queue
       INTEGER :: mytid, ntids! omp variables
 
       trcoptparttime = MPI_WTIME() ! cronometer-start
@@ -33,48 +33,39 @@ SUBROUTINE trcopt
 !!!$omp&                       shared(jj,jpk,jpj,jpi,xpar,conversion,kef,e3t,qsr)
 
 ! 1. determination of surface irradiance
-        DO ji = 1,jpi
-      DO jj = 1,jpj
-
-          zpar0m(ji)          = qsr(jj,ji)*conversion
-          zpar100(ji)         = zpar0m(ji)*0.01
-          xpar(1,jj,ji) = zpar0m(ji)
-          zpar(1,ji)          = zpar0m(1)
-          xEPS_ogstm(1,ji)          = kef(jj,ji)
-
-      !   ENDDO
-      ! ENDDO
-!! 2. determination of xpar
-!! ------------------------
-            !       DO ji = 1,jpi
-            ! DO jj = 1,jpj
-
-        DO jk = 2,jpk
-            xEPS_ogstm(jk,ji) = max(kef(jj,ji),1.D-15) ! avoid denormalized number
+      queue=1
+      ! XXX: this can be improved but no way to test as it doesn't seem to
+      ! impact numerical results
+      !$acc kernels default(present) async(queue)
+      DO ji = 1,jpi
+         DO jj = 1,jpj
+
+            ! zpar0m(ji)          = qsr(jj,ji)*conversion
+            ! zpar100(ji)         = zpar0m(ji)*0.01
+            xpar(1,jj,ji) = qsr(jj,ji)*conversion
+            ! zpar(1,ji)          = zpar0m(1)
+            xEPS_ogstm(1,ji)          = kef(jj,ji)
+
+            !! 2. determination of xpar
+            !! ------------------------
+
+            DO jk = 2,jpk
+               xEPS_ogstm(jk,ji) = max(kef(jj,ji),1.D-15) ! avoid denormalized number
+            END DO
+
+            DO jk = 2,jpk
+               xpar(jk,jj,ji) = max( xpar(jk-1,jj,ji) *exp(-1. * xEPS_ogstm(jk-1,ji)* e3t(jk-1,jj,ji)) ,1.D-15) ! avoid denormalized number
+
+            END DO
+
+            DO jk = 1,jpk
+               xpar(jk,jj,ji) = max( xpar(jk,jj,ji) * exp(- xEPS_ogstm(jk,ji)* 0.5D+00* e3t(jk,jj,ji) ) ,1.D-15)
+            END DO
+
          END DO
-      !   END DO
-      ! ENDDO
-
-            !       DO ji = 1,jpi
-            ! DO jj = 1,jpj
-      DO jk = 2,jpk
-            !print * ,"CHECK",xpar(jk-1,jj,ji),xEPS_ogstm(jk-1,ji),
-            xpar(jk,jj,ji) = max( xpar(jk-1,jj,ji) *exp(-1. * xEPS_ogstm(jk-1,ji)* e3t(jk-1,jj,ji)) ,1.D-15) ! avoid denormalized number
-
-          END DO
-      !   END DO
-      ! ENDDO
-
-            !       DO ji = 1,jpi
-            ! DO jj = 1,jpj
-        DO jk = 1,jpk
-            xpar(jk,jj,ji) = max( xpar(jk,jj,ji) * exp(- xEPS_ogstm(jk,ji)* 0.5D+00* e3t(jk,jj,ji) ) ,1.D-15)
-          END DO
-
-        END DO
       ENDDO
-
-    
+      !$acc end kernels
+      !$acc wait(queue)
 
        trcoptparttime = MPI_WTIME() - trcoptparttime ! cronometer-stop
        trcopttottime = trcopttottime + trcoptparttime
diff --git a/src/BIO/trcsed.f90 b/src/BIO/trcsed.f90
index b9c5c484..2295e64b 100644
--- a/src/BIO/trcsed.f90
+++ b/src/BIO/trcsed.f90
@@ -61,8 +61,8 @@ SUBROUTINE trcsed
 #ifdef key_trc_bfm
 
       LOGICAL :: l1,l2,l3
-      INTEGER :: ji,jj,jk,jv,jf,js
-      INTEGER :: bottom
+      INTEGER :: ji,jj,jk,jv,jf,js,ntx
+      INTEGER :: bottom,queue
       double precision :: ze3tr,d2s
 ! omp variables
     
@@ -71,7 +71,7 @@ SUBROUTINE trcsed
 !! ===================
 
 
-
+      queue=1
 
       d2s=1./3600./24.  ! speed from (m/day) to  (m/s)
 
@@ -100,107 +100,200 @@ SUBROUTINE trcsed
                END DO
             END DO
          END DO
-
+#ifdef _OPENACC
+         call myalloc_SED_gpu()
+#endif
+         !$acc update device(jarr_sed,jarr_sed_flx)
       ENDIF ! End initialization phase (once at the beginning)
 
 
 ! vertical slab
 ! =============
 
+      !      if( mytid + jv <=  dimen_jvsed) then
+      ! 1. sedimentation of detritus  : upstream scheme
+      ! -----------------------------------------------
+      ! 1.1 initialisation needed for bottom and surface value
+
+      !$acc parallel loop gang vector collapse(3) default(present) async(queue)
+      DO jv=1,dimen_jvsed
+         DO js = 1, nsed
+            DO  jk = 1,jpk
+#ifdef _OPENACC
+               ntx=jv
+#else
+               ntx=1
+#endif
+               ji = jarr_sed(2,jv)
+               jj = jarr_sed(1,jv)
+               zwork(jk,js,ntx) = 0.
+
+            END DO
+         END DO
+      END DO
+
+      ! 1.2 tracer flux at w-point: we use -vsed (downward flux)
+      ! with simplification : no e1*e2
+
+      !                Particulate
+      !$acc parallel loop gang vector collapse(3) default(present) async(queue)
+      DO jv=1,dimen_jvsed
+         DO js =1,4
+            DO  jk = 2,jpkm1
+#ifdef _OPENACC
+               ntx=jv
+#else
+               ntx=1
+#endif
+               ji = jarr_sed(2,jv)
+               jj = jarr_sed(1,jv)
+               zwork(jk,js,jv) = -vsed * trn(jk-1,jj,ji, sed_idx(js))
+            END DO
+         END DO
+      END DO
+
+      !                Diatoms
+      !$acc parallel loop gang vector collapse(3) default(present) async(queue)
+      DO jv=1,dimen_jvsed
+         DO js =5,9
+            DO  jk = 2,jpkm1
+#ifdef _OPENACC
+               ntx=jv
+#else
+               ntx=1
+#endif
+               ji = jarr_sed(2,jv)
+               jj = jarr_sed(1,jv)
+               zwork(jk,js,ntx) = -ogstm_sedipi(jk-1,jj,ji,1) * trn(jk-1,jj,ji, sed_idx(js))
+            END DO
+         END DO
+      END DO
+
+      !                Flagellates
+      !$acc parallel loop gang vector collapse(3) default(present) async(queue)
+      DO jv=1,dimen_jvsed
+         DO js =10,13
+            DO  jk = 2,jpkm1
+#ifdef _OPENACC
+               ntx=jv
+#else
+               ntx=1
+#endif
+               ji = jarr_sed(2,jv)
+               jj = jarr_sed(1,jv)
+               zwork(jk,js,ntx) = -ogstm_sedipi(jk-1,jj,ji,2) * trn(jk-1,jj,ji, sed_idx(js))
+            END DO
+         END DO
+      END DO
+
+      !                Picophytoplankton
+      !$acc parallel loop gang vector collapse(3) default(present) async(queue)
+      DO jv=1,dimen_jvsed
+         DO js =14,17
+            DO  jk = 2,jpkm1
+#ifdef _OPENACC
+               ntx=jv
+#else
+               ntx=1
+#endif
+               ji = jarr_sed(2,jv)
+               jj = jarr_sed(1,jv)
+               zwork(jk,js,ntx) = -ogstm_sedipi(jk-1,jj,ji,3) * trn(jk-1,jj,ji, sed_idx(js))
+            END DO
+         END DO
+      END DO
+
+      !                Dinoflagellates
+      !$acc parallel loop gang vector collapse(3) default(present) async(queue)
+      DO jv=1,dimen_jvsed
+         DO js =18,21
+            DO  jk = 2,jpkm1
+#ifdef _OPENACC
+               ntx=jv
+#else
+               ntx=1
+#endif
+               ji = jarr_sed(2,jv)
+               jj = jarr_sed(1,jv)
+               zwork(jk,js,ntx) = -ogstm_sedipi(jk-1,jj,ji,4) * trn(jk-1,jj,ji, sed_idx(js))
+            END DO
+         END DO
+      END DO
+
+      !                Calcite
+      !$acc parallel loop gang vector collapse(3) default(present) async(queue)
+      DO jv=1,dimen_jvsed
+         DO js =22,22
+            DO  jk = 2,jpkm1
+#ifdef _OPENACC
+               ntx=jv
+#else
+               ntx=1
+#endif
+               ji = jarr_sed(2,jv)
+               jj = jarr_sed(1,jv)
+               zwork(jk,js,ntx) = - vsedO5c * trn(jk-1,jj,ji, sed_idx(js))
+            END DO
+         END DO
+      END DO
+
+      !$acc parallel loop gang vector collapse(2) default(present) async(queue)
+      DO jv=1,dimen_jvsed
+         DO js = 1,nsed
+#ifdef _OPENACC
+            ntx=jv
+#else
+            ntx=1
+#endif
+            ji = jarr_sed(2,jv)
+            jj = jarr_sed(1,jv)
+            bottom = mbathy(jj,ji) + 1
+            zwork(bottom,js,ntx) = bottom_flux * zwork(bottom,js,ntx) ! bottom_flux = 0 -> no flux in the sea floor
+         END DO
+      END DO
+
+      ! 1.3 tracer flux divergence at t-point added to the general trend
 
+      !$acc parallel loop gang vector collapse(2) default(present) async(queue)
+      DO jv=1,dimen_jvsed
+         DO  jk = 1,jpkm1
+#ifdef _OPENACC
+            ntx=jv
+#else
+            ntx=1
+#endif
+            ji = jarr_sed(2,jv)
+            jj = jarr_sed(1,jv)
+
+            jf=  jarr_sed_flx(jk,jV)
+
+            ze3tr = 1./e3t(jk,jj,ji)
+
+            DO js =1,nsed
+               ztra(js,ntx) = -ze3tr * (zwork(jk,js,ntx) - zwork(jk+1,js,ntx))
+               IF ((Fsize .GT. 0) .AND. (jf .GT. 0)) THEN
+                  diaflx(4,jf,sed_idx(js)) = diaflx(4, jf,sed_idx(js)) + zwork(jk,js,ntx)*rdt
+               ENDIF
 
-      MAIN_LOOP: DO jv=1,dimen_jvsed
-
-!      if( mytid + jv <=  dimen_jvsed) then
-! 1. sedimentation of detritus  : upstream scheme
-! -----------------------------------------------
-! 1.1 initialisation needed for bottom and surface value
-
-           ji = jarr_sed(2,jv)
-           jj = jarr_sed(1,jv)
-
-
-                 DO js = 1, nsed
-              DO  jk = 1,jpk
-
-                    zwork(jk,js,1) = 0.
-
-                 END DO
-              END DO
-
-! 1.2 tracer flux at w-point: we use -vsed (downward flux)
-! with simplification : no e1*e2
-
-             
-
-!                Particulate
-              DO js =1,4
-                 DO  jk = 2,jpkm1
-                    zwork(jk,js,1) = -vsed * trn(jk-1,jj,ji, sed_idx(js))
-                 END DO
-              END DO
-!                Diatoms
-              DO js =5,9
-                 DO  jk = 2,jpkm1
-                    zwork(jk,js,1) = -ogstm_sedipi(jk-1,jj,ji,1) * trn(jk-1,jj,ji, sed_idx(js))
-                 END DO
-              END DO
-!                Flagellates
-              DO js =10,13
-                 DO  jk = 2,jpkm1
-                    zwork(jk,js,1) = -ogstm_sedipi(jk-1,jj,ji,2) * trn(jk-1,jj,ji, sed_idx(js))
-                 END DO
-              END DO
-!                Picophytoplankton
-              DO js =14,17
-                 DO  jk = 2,jpkm1
-                    zwork(jk,js,1) = -ogstm_sedipi(jk-1,jj,ji,3) * trn(jk-1,jj,ji, sed_idx(js))
-                 END DO
-              END DO
-!                Dinoflagellates
-              DO js =18,21
-                 DO  jk = 2,jpkm1
-                    zwork(jk,js,1) = -ogstm_sedipi(jk-1,jj,ji,4) * trn(jk-1,jj,ji, sed_idx(js))
-                 END DO
-              END DO
-
-!                Calcite
-              DO js =22,22
-                 DO  jk = 2,jpkm1
-                    zwork(jk,js,1) = - vsedO5c * trn(jk-1,jj,ji, sed_idx(js))
-                 END DO
-              END DO
-               bottom = mbathy(jj,ji) + 1
-               zwork(bottom,:,1) = bottom_flux * zwork(bottom,:,1) ! bottom_flux = 0 -> no flux in the sea floor
-
-! 1.3 tracer flux divergence at t-point added to the general trend
-
-              DO  jk = 1,jpkm1
-                  jf=  jarr_sed_flx(jk,jV)
-
-                 ze3tr = 1./e3t(jk,jj,ji)
-
-                 DO js =1,nsed
-                    ztra(js,1) = -ze3tr * (zwork(jk,js,1) - zwork(jk+1,js,1))
-                    IF ((Fsize .GT. 0) .AND. (jf .GT. 0)) THEN
-                         diaflx(4,jf,sed_idx(js)) = diaflx(4, jf,sed_idx(js)) + zwork(jk,js,1)*rdt
-                    ENDIF
-                 END DO
-
-                 DO js =1,nsed
 !!!  d2s convert speed from (m/day) to  (m/s)
-                    tra(jk,jj,ji,sed_idx(js)) = tra(jk,jj,ji,sed_idx(js)) + ztra(js,1)*d2s
-                 END DO
+               tra(jk,jj,ji,sed_idx(js)) = tra(jk,jj,ji,sed_idx(js)) + ztra(js,ntx)*d2s
+            END DO
+         END DO
+      END DO
+      !$acc end parallel loop
+      !$acc wait(queue)
 
 #ifdef key_trc_diabio
-                  trbio(jk,jj,ji,8) = ztra
-#endif
-
-              END DO
+      DO jv=1,dimen_jvsed
+         ji = jarr_sed(2,jv)
+         jj = jarr_sed(1,jv)
+         DO  jk = 1,jpkm1
+            trbio(jk,jj,ji,8) = ztra
+         END DO
+      END DO
 #endif
 
-
-      END DO MAIN_LOOP
+#endif ! key_trc_bfm
 
 !!!$omp    end parallel do
 
diff --git a/src/BIO/trcsms.f90 b/src/BIO/trcsms.f90
index 7c463d54..1b4faff9 100644
--- a/src/BIO/trcsms.f90
+++ b/src/BIO/trcsms.f90
@@ -22,6 +22,7 @@ SUBROUTINE trcsms
 
        USE myalloc
        USE mpi
+
        IMPLICIT NONE
 
 
@@ -29,11 +30,10 @@ SUBROUTINE trcsms
 
        trcsmsparttime = MPI_WTIME() ! cronometer-start
 
-
 !! this first routines are parallelized on vertical slab
 
        CALL trcopt ! tracers: optical model
-       
+
        CALL trcbio ! tracers: biological model
 
 !! trcsed no updated for time step advancing
diff --git a/src/General/TimeManager.f90 b/src/General/TimeManager.f90
index e095e8de..ba36df25 100644
--- a/src/General/TimeManager.f90
+++ b/src/General/TimeManager.f90
@@ -799,6 +799,7 @@ END FUNCTION INSTANT_PAR
 
 
        double precision FUNCTION INSTANT_PAR_from_sec(sec, MEAN_PAR)
+       !$acc routine seq
        IMPLICIT NONE
        double precision, INTENT(IN) :: MEAN_PAR ! daily integral 
        ! LOCAL
diff --git a/src/General/memory.f90 b/src/General/memory.f90
index 2af71f75..15a6a3cf 100644
--- a/src/General/memory.f90
+++ b/src/General/memory.f90
@@ -516,7 +516,7 @@ subroutine alloc_tot()
       gphiv    = huge(gphiv(1,1))
       allocate(gphif(jpj,jpi))          
       gphif    = huge(gphif(1,1))
-      allocate(e1t(jpj,jpi))            
+      allocate(e1t(jpj,jpi))
       e1t      = huge(e1t(1,1))
       allocate(e1u(jpj,jpi))            
       e1u      = huge(e1u(1,1))
@@ -542,6 +542,7 @@ subroutine alloc_tot()
      allocate(gdept(jpk,jpj,jpi))
        gdept = huge(gdept(1,1,1))
 #endif
+      !$acc enter data create(gdept)
 
       allocate(gdepw(jpk)) 
         gdepw = huge(gdepw(1))
@@ -567,7 +568,8 @@ subroutine alloc_tot()
         e3w = huge(e3w(1,1,1))
 
       allocate(mbathy(jpj,jpi)) 
-       mbathy = huge(mbathy(1,1))
+      mbathy = huge(mbathy(1,1))
+      !$acc enter data create(mbathy)
 
       allocate(tmask(jpk,jpj,jpi)) 
       tmask = huge(tmask(1,1,1))
@@ -586,6 +588,7 @@ subroutine alloc_tot()
 
       allocate(bfmmask(jpk, jpj, jpi))
       bfmmask = huge(bfmmask(1, 1, 1))
+      !$acc enter data create(bfmmask)
 
        allocate(un(jpk,jpj,jpi))    
       un     = huge(un(1,1,1))
@@ -603,6 +606,7 @@ subroutine alloc_tot()
       rhopn  = huge(rhopn(1,1,1))
        allocate(rho(jpk,jpj,jpi))   
       rho    = huge(rho(1,1,1))
+      !$acc enter data create(rho)
       
       allocate(ahtu(jpk)) 
       ahtu = huge(ahtu(1))
@@ -781,7 +785,8 @@ subroutine alloc_tot()
 #    endif
 
 #ifdef key_trc_bfm
-      allocate(xpar(jpk,jpj,jpi))   
+      allocate(xpar(jpk,jpj,jpi))
+      !$acc enter data create(xpar)
        xpar = huge(xpar(1,1,1))
 
 #endif
@@ -789,11 +794,19 @@ subroutine alloc_tot()
 !!    photoperiod
         allocate(DAY_LENGTH(jpj,jpi))   
        DAY_LENGTH = huge(DAY_LENGTH(1,1))
+       !$acc enter data create(DAY_LENGTH)
        forcing_phys_initialized = .false.
+
+       !$acc enter data create(e1t,e2t,e3t,e3w,e3t_back,tra,trb,tmask,umask,&
+       !$acc& vmask,avt,e1u,e2u,e3u,e1v,e2v,e3v,un,vn,wn,trn,ahtt,traIO,traIO_HIGH,&
+       !$acc& snIO,tnIO,wnIO,avtIO,e3tIO,unIO,vnIO,vatmIO,empIO,qsrIO,sn,tn,&
+       !$acc& tra_DIA,tra_DIA_IO,tra_DIA_2d_IO,tra_DIA_2d,tra_DIA_IO_HIGH,&
+       !$acc& vatm,emp,qsr,tra_DIA_2d_IO_HIGH)
+
 #ifdef Mem_Monitor
       mem_all=get_mem(err) - aux_mem
 #endif
-  
+
         END subroutine alloc_tot
 
 
@@ -802,6 +815,8 @@ subroutine clean_memory()
 
             ! myalloc (memory.f90)
 
+            !$acc exit data delete(e1t,e2t,e3t,e3w,e3t_back,tra,trb,tmask,umask,vmask,avt)
+
 #ifdef key_mpp
 
             !$acc exit data delete(te_send, tw_send, tn_send, ts_send) finalize
@@ -860,6 +875,7 @@ subroutine clean_memory()
             deallocate(ff)
             
             deallocate(gdept)
+            !$acc exit data delete(gdept)
             deallocate(gdepw)
             deallocate(e3t_0)
             deallocate(e3u_0)
@@ -873,6 +889,7 @@ subroutine clean_memory()
             deallocate(e3w)
             
             deallocate(mbathy)
+            !$acc exit data delete(mbathy)
             
             deallocate(tmask)
             deallocate(h_column)
@@ -884,6 +901,7 @@ subroutine clean_memory()
             deallocate(vmask)
             
             deallocate(bfmmask)
+            !$acc exit data delete(bfmmask)
             
             deallocate(un)
             deallocate(vn)
@@ -893,6 +911,7 @@ subroutine clean_memory()
             deallocate(rdn)
             deallocate(rhopn)
             deallocate(rho)
+            !$acc exit data delete(rho)
             
             deallocate(ahtu)
             deallocate(ahtv)
@@ -953,7 +972,7 @@ subroutine clean_memory()
             deallocate(tra_DIA_2d_IO_HIGH)
             deallocate(tra_PHYS_2d_IO)
             deallocate(tra_PHYS_2d_IO_HIGH)
-            
+
 
             if(lwp) then
                 deallocate(tottrn)
@@ -985,16 +1004,23 @@ subroutine clean_memory()
 
 #ifdef key_trc_bfm
             deallocate(xpar)
+            !$acc exit data delete(xpar)
 #endif
 
             deallocate(DAY_LENGTH)
+            !$acc exit data delete(DAY_LENGTH)
 
             ! trclec
 
+            !$acc exit data delete(highfreq_table,highfreq_table_dia,highfreq_table_dia2d)
             deallocate(highfreq_table)
             deallocate(highfreq_table_dia)
             deallocate(highfreq_table_dia2d)
 
+            !$acc exit data delete(trn, e1u, e2u, e3u, e1v, e2v, e3v, un, vn, wn,&
+            !$acc& ahtt, traio,traIO_HIGH,snIO,tnIO,wnIO,avtIO,e3tIO,unIO,vnIO,vatmIO,empIO,qsrIO,sn,tn,&
+            !$acc& tra_DIA,tra_DIA_IO,tra_DIA_2d_IO,tra_DIA_2d,tra_DIA_IO_HIGH,vatm,emp,qsr,tra_DIA_2d_IO_HIGH)
+
         end subroutine clean_memory
 
     INTEGER FUNCTION find_index_var(string)
diff --git a/src/General/step.f90 b/src/General/step.f90
index 80c8ef31..66f96be6 100644
--- a/src/General/step.f90
+++ b/src/General/step.f90
@@ -265,6 +265,13 @@ SUBROUTINE trcstp
 !         with surface boundary condition
 !         with IMPLICIT vertical diffusion
 
+      ! XXX: to be removed
+      use DIA_mem, only: diaflx,flx_ridxt
+      use myalloc, only: tra,trb,e1t,e3t_back,e2t,e3t,e3w,umask,vmask,tmask,avt,ahtt
+      use BIO_mem, only: ogstm_sediPI,ogstm_PH,ogstm_co2
+      USE OPT_mem, only: kef
+      USE SED_mem
+
        IMPLICIT NONE
       integer jn,jk,ji,jj
       trcstpparttime = MPI_WTIME() ! cronometer-start
@@ -289,14 +296,26 @@ SUBROUTINE trcstp
 ! tracers: horizontal diffusion IF namelist flags are activated
 ! -----------------------------
 
+      !$acc update device(umask,vmask,tmask,trb,ahtt,tra,diaflx,flx_ridxt) if(lhdf)
       IF (lhdf) CALL trchdf
+      !$acc update host(diaflx,tra) if(lhdf)
 
 ! tracers: sink and source (must be  parallelized on vertical slab)
       IF (lsbc) CALL trcsbc ! surface cell processes, default lsbc = False
 
+       !$acc update device(kef,qsr,mbathy,bfmmask,trn,DAY_LENGTH,vatm,tn,sn,rho,e3t,gdept,ogstm_PH,ogstm_co2) if(lbfm)
+#if  defined key_trc_sed
+       !$acc update device(sed_idx,diaflx,e3t,tra,ogstm_sedipi,mbathy) if(lbfm)
+#endif
       IF (lbfm) CALL trcsms
+       !$acc update host(tra,tra_DIA,tra_DIA_2d,ogstm_sediPI,ogstm_PH) if(lbfm)
+#if  defined key_trc_sed
+       !$acc update host(diaflx,zwork) if(lbfm)
+#endif
 
+      !$acc update device(e1t,diaflx,e3t_back,e2t,trb,tmask,e3t,tra,avt,e3w) if (lzdf)
       IF (lzdf) CALL trczdf ! tracers: vertical diffusion
+      !$acc update host(diaflx,tra) if (lzdf)
 
       IF (lsnu) CALL snutel
 
diff --git a/src/IO/DIA_mem.f90 b/src/IO/DIA_mem.f90
index e385f32e..4b4a896a 100644
--- a/src/IO/DIA_mem.f90
+++ b/src/IO/DIA_mem.f90
@@ -50,6 +50,7 @@ SUBROUTINE alloc_DIA_local_flx()
        INDflxDUMP = huge(INDflxDUMP(1))
        allocate(diaflx        (7, Fsize, jptra   ))
        diaflx = 0
+       !$acc enter data create(flx_ridxt,diaflx)
       END SUBROUTINE alloc_DIA_local_flx
 
 
@@ -95,6 +96,7 @@ subroutine clean_memory_dia()
 
           if (allocated(diaflx)) then
               deallocate(diaflx)
+              !$acc exit data delete(diaflx)
           endif
 
           if (allocated(INDflxDUMPZERO)) then
diff --git a/src/MPI/ogstm_mpi.f90 b/src/MPI/ogstm_mpi.f90
index b6c0bcd4..1cf8a3e2 100644
--- a/src/MPI/ogstm_mpi.f90
+++ b/src/MPI/ogstm_mpi.f90
@@ -118,168 +118,7 @@ SUBROUTINE mynode
 !!                    nono   : number for local neighboring processors
 !!
 !!----------------------------------------------------------------------
- SUBROUTINE mpplnk_my(ptab)
-
-      double precision ptab(jpk,jpj,jpi)
-
-
-#ifdef key_mpp_mpi
-
-      INTEGER jk,jj,ji
-      INTEGER reqs1, reqs2, reqr1, reqr2
-      INTEGER reqs3, reqs4, reqr3, reqr4
-      INTEGER jw, packsize
-
-!!     trcadvparttime = MPI_WTIME()
-
-!!
-!!2. East and west directions exchange
-!!------------------------------------
-
-
-
-!!
-!!2.2 Migrations
-!!
-!!
-!               3     4
-!               |     ^
-!               |     |
-!               v     |
-!           ________________
-!          |                |
-!     1<-- |                | 1 <--
-!     2--> |                | 2 -->
-!          |________________|
-!               3     4
-!               |     ^
-!               |     |
-!               v     |
-
-    packsize=jpk*jpj
-
-      IF(nbondi.eq.-1) THEN ! We are at the west side of the domain
-
-          CALL mppsend(2,ptab(:,:,jpi-1),packsize,noea,0,reqs1)
-          CALL mpprecv(1,ptab(:,:,  jpi),packsize,reqr1)
-
-      ELSE IF(nbondi.eq.0) THEN
-          CALL mppsend(1, ptab(:,:    ,2),packsize,nowe,0,reqs1)
-          CALL mppsend(2, ptab(:,:,jpi-1),packsize,noea,0,reqs2)
-
-          CALL mpprecv(1,ptab(:,:,jpi),packsize,reqr1)
-          CALL mpprecv(2,ptab(:,:,  1),packsize,reqr2)
-
-      ELSE IF(nbondi.eq.1) THEN ! We are at the east side of the domain
-
-          CALL mppsend(1,ptab(:,:,2), packsize, nowe,0, reqs1)
-          CALL mpprecv(2,ptab(:,:,1), packsize, reqr1)
-
-
-      ENDIF
-
-
-!!
-!!
-!!3. North and south directions
-!!-----------------------------
-!!
-!!3.1 Read Dirichlet lateral conditions
-!!
-
-
-      IF(nbondj.eq.0.or.nbondj.eq.-1) THEN
-         DO jw=1,NORTH_count_send
-              ji = NORTHpoints_send(1,jw)
-              jk = NORTHpoints_send(2,jw)
-              tn_send(jw) = ptab(jk,jpj-1,ji)
-         ENDDO
-     ENDIF
-     IF(nbondj.eq.0.or.nbondj.eq.1) THEN
-         DO jw=1,SOUTH_count_send
-             ji = SOUTHpoints_send(1,jw)
-             jk = SOUTHpoints_send(2,jw)
-             ts_send(jw) = ptab(jk,2,ji)
-         ENDDO
-
-
-      ENDIF!    PACK_LOOP4
-
-
-!!
-!!2.2 Migrations
-!!
-!!
-
-      IF(nbondj.eq.-1) THEN ! We are at the south side of the domain
-          CALL mppsend(4,tn_send,NORTH_count_send,nono,0,reqs4)
-          CALL mpprecv(3,tn_recv,NORTH_count_recv,reqr3)
-          CALL mppwait(reqs4)
-          CALL mppwait(reqr3)
-      ELSE IF(nbondj.eq.0) THEN
-          CALL mppsend(4, tn_send,NORTH_count_send,nono,0,reqs4)
-          CALL mppsend(3, ts_send,SOUTH_count_send,noso,0,reqs3)
-          CALL mpprecv(3,tn_recv,NORTH_count_recv,reqr3)
-          CALL mpprecv(4,ts_recv,SOUTH_count_recv,reqr4)
-
-          CALL mppwait(reqs4)
-          CALL mppwait(reqs3)
-          CALL mppwait(reqr3)
-          CALL mppwait(reqr4)
-      ELSE IF(nbondj.eq.1) THEN ! We are at the north side of the domain
-          CALL mppsend(3,ts_send, SOUTH_count_send, noso,0, reqs3)
-          CALL mpprecv(4,ts_recv, SOUTH_count_recv, reqr4)
-          CALL mppwait(reqs3)
-          CALL mppwait(reqr4)
-      ENDIF
-
-
-
-!!
-!!2.3 Write Dirichlet lateral conditions
-!!
-
-      IF(nbondj.eq.0.or.nbondj.eq.1) THEN ! All but south boundary, we received from south
-
-         DO jw=1,SOUTH_count_recv
-              ji = SOUTHpoints_recv(1,jw)
-              jk = SOUTHpoints_recv(2,jw)
-              ptab(jk,1,ji)= ts_recv(jw)
-         ENDDO
-
-      ENDIF
-
-      IF(nbondj.eq.-1.or.nbondj.eq.0) THEN ! All but north boundary, we received from north
-
-        DO jw=1,NORTH_count_recv
-              ji = NORTHpoints_recv(1,jw)
-              jk = NORTHpoints_recv(2,jw)
-             ptab(jk,jpj,ji)= tn_recv(jw)
-        ENDDO
-
-      ENDIF ! PACK_LOOP5
-
-
-!!!  East - West waits
-
-      IF(nbondi.eq.-1) THEN ! We are at the west side of the domain
-          CALL mppwait(reqs1)
-          CALL mppwait(reqr1)
-      ELSE IF(nbondi.eq.0) THEN
-          CALL mppwait(reqs1)
-          CALL mppwait(reqs2)
-          CALL mppwait(reqr1)
-          CALL mppwait(reqr2)
-      ELSE IF(nbondi.eq.1) THEN ! We are at the east side of the domain
-          CALL mppwait(reqs1)
-          CALL mppwait(reqr1)
-      ENDIF
-
-#endif
-
-END SUBROUTINE
-
-SUBROUTINE mpplnk_my_openacc(ptab,gpu)
+SUBROUTINE mpplnk_my(ptab,gpu)
 
       double precision ptab(jpk,jpj,jpi)
 #ifdef key_mpp_mpi
@@ -413,6 +252,7 @@ SUBROUTINE mpplnk_my_openacc(ptab,gpu)
       !$acc kernels default(present) if(use_gpu)
       IF(nbondj.eq.0.or.nbondj.eq.1) THEN ! All but south boundary, we received from south
 
+         !$acc loop independent
          DO jw=1,SOUTH_count_recv
               ji = SOUTHpoints_recv(1,jw)
               jk = SOUTHpoints_recv(2,jw)
@@ -423,6 +263,7 @@ SUBROUTINE mpplnk_my_openacc(ptab,gpu)
 
       IF(nbondj.eq.-1.or.nbondj.eq.0) THEN ! All but north boundary, we received from north
 
+        !$acc loop independent
         DO jw=1,NORTH_count_recv
               ji = NORTHpoints_recv(1,jw)
               jk = NORTHpoints_recv(2,jw)
diff --git a/src/PHYS/ADV_mem.f90 b/src/PHYS/ADV_mem.f90
index 9d201de7..d1bba762 100644
--- a/src/PHYS/ADV_mem.f90
+++ b/src/PHYS/ADV_mem.f90
@@ -114,6 +114,9 @@ subroutine myalloc_ADV()
       mem_all=get_mem(err) - aux_mem
 #endif
 
+      !$acc enter data create(advmask,zaa,zbb,zcc,inv_eu,inv_ev,inv_et,&
+      !$acc& big_fact_zaa,big_fact_zbb,big_fact_zcc,zbtr_arr)
+
       END subroutine myalloc_ADV
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
 
@@ -133,6 +136,9 @@ subroutine clean_memory_adv()
           deallocate(big_fact_zcc)
           deallocate(zbtr_arr)
 
+          !$acc exit data delete(advmask,zaa,zbb,zcc,inv_eu,inv_ev,inv_et,&
+          !$acc& big_fact_zaa,big_fact_zbb,big_fact_zcc,zbtr_arr)
+
       end subroutine clean_memory_adv
 
 
diff --git a/src/PHYS/ZDF_mem.f90 b/src/PHYS/ZDF_mem.f90
index cb3fe00b..265e7004 100644
--- a/src/PHYS/ZDF_mem.f90
+++ b/src/PHYS/ZDF_mem.f90
@@ -46,6 +46,9 @@ subroutine myalloc_ZDF()
         jarr_zdf     = huge(jarr_zdf(1,1))
        allocate(jarr_zdf_flx(jpi*jpj,jpk)) 
         jarr_zdf_flx = huge(jarr_zdf_flx(1,1))
+        !$acc enter data create(jarr_zdf,jarr_zdf_flx)
+        !$acc update device(jarr_zdf,jarr_zdf_flx)
+#ifndef _OPENACC
        allocate(zwd(jpk, ntids))           
         zwd          = huge(zwd(1,1))
        allocate(zws(jpk, ntids))           
@@ -60,6 +63,7 @@ subroutine myalloc_ZDF()
         zwz          = huge(zwz(1,1))
        allocate(zwt(jpk, ntids))           
         zwt          = huge(zwt(1,1))
+#endif
 
 #ifdef Mem_Monitor
       mem_all=get_mem(err) - aux_mem
@@ -67,7 +71,28 @@ subroutine myalloc_ZDF()
 
 
       END subroutine myalloc_ZDF
-      
+
+#ifdef _OPENACC
+      subroutine myalloc_ZDF_gpu()
+        allocate(zwd(jpk, dimen_jvzdf))
+        zwd          = huge(zwd(1,1))
+        allocate(zws(jpk, dimen_jvzdf))
+        zws          = huge(zws(1,1))
+        allocate(zwi(jpk, dimen_jvzdf))
+        zwi          = huge(zwi(1,1))
+        allocate(zwx(jpk, dimen_jvzdf))
+        zwx          = huge(zwx(1,1))
+        allocate(zwy(jpk, dimen_jvzdf))
+        zwy          = huge(zwy(1,1))
+        allocate(zwz(jpk, dimen_jvzdf))
+        zwz          = huge(zwz(1,1))
+        allocate(zwt(jpk, dimen_jvzdf))
+        zwt          = huge(zwt(1,1))
+
+        !$acc enter data create(zwd,zwi,zwx,zws,zwz,zwy,zwt)
+        !$acc update device(zwd,zwi,zwx,zws,zwz,zwy,zwt)
+      END subroutine myalloc_ZDF_gpu
+#endif
       
       
       subroutine clean_memory_zdf()
@@ -82,6 +107,8 @@ subroutine clean_memory_zdf()
           deallocate(zwz)
           deallocate(zwt)
 
+          !$acc exit data delete(jarr_zdf,jarr_zdf_flx,zwd,zwi,zwx,zws,zwz,zwy,zwt)
+
       end subroutine clean_memory_zdf
 
 
diff --git a/src/PHYS/trcadv.f90 b/src/PHYS/trcadv.f90
index 129613be..9fa07701 100644
--- a/src/PHYS/trcadv.f90
+++ b/src/PHYS/trcadv.f90
@@ -1,3 +1,8 @@
+#ifdef _OPENACC
+! BUG?: the fsx routine causes additional H2D copies
+#define fsx(pfx1, pfx2, pfu) ((((pfu) + abs(pfu)) * (pfx1) + ((pfu) - abs(pfu)) * (pfx2)) * 0.5)
+#endif
+
 SUBROUTINE trcadv
 !!!---------------------------------------------------------------------
 !!!
@@ -68,7 +73,7 @@ SUBROUTINE trcadv
   !!                               monthly weather review, pp 479-486
   !!
   LOGICAL :: MPI_CHECK,l1,l2,l3
-  INTEGER :: jk,jj,ji,jt,jn,jf,ju
+  INTEGER :: jk,jj,ji,jt,jn,jf,ju,queue
   double precision :: zbtr,zdt
   double precision :: junk, junki, junkj, junkk
   double precision :: timer
@@ -78,10 +83,16 @@ SUBROUTINE trcadv
   double precision, allocatable,dimension(:,:,:) :: zx,zy,zz,zbuf
   double precision, allocatable,dimension(:,:,:) :: zkx,zky,zkz
   logical :: use_gpu
+
+  queue=1
   
   trcadvparttime = MPI_WTIME()
 
+#ifdef _OPENACC
   use_gpu=.true.
+#else
+  use_gpu=.false.
+#endif
 
   !-------------------------------------------------------------------
 
@@ -159,8 +170,22 @@ SUBROUTINE trcadv
 
 
      adv_initialized=.true.
+
   endif
 
+  !!OpenMP compatibility broken. Possibility to use ifndef OpenMP + rename the file in trcadv.F90 to keep it
+  allocate(zy(jpk,jpj,jpi))
+  allocate(zx(jpk,jpj,jpi))
+  allocate(zz(jpk,jpj,jpi))
+  allocate(ztj(jpk,jpj,jpi))
+  allocate(zti(jpk,jpj,jpi))
+  allocate(zkx(jpk,jpj,jpi))
+  allocate(zky(jpk,jpj,jpi))
+  allocate(zkz(jpk,jpj,jpi))
+  allocate(zbuf(jpk,jpj,jpi))
+
+  !$acc enter data create(zy,zx,zz,ztj,zti,zkx,zky,zkz,zbuf)
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! end initialization phase
 
   jk=1
@@ -168,30 +193,24 @@ SUBROUTINE trcadv
   zdt = rdt*ndttrc
   !$OMP TASK private(ji,jj) firstprivate(jpim1,jpjm1) shared(zbtr_arr,e1t,e2t,e3t) default(none)
 
+  !$acc update device( zaa(1:jpk,1:jpj,1:jpi), zbb(1:jpk,1:jpj,1:jpi), zcc(1:jpk,1:jpj,1:jpi) )
+  !$acc update device( inv_eu(1:jpk,1:jpj,1:jpi), inv_ev(1:jpk,1:jpj,1:jpi), inv_et(1:jpk,1:jpj,1:jpi) )
+  !$acc update device( big_fact_zaa (1:jpk,1:jpj,1:jpi), big_fact_zbb(1:jpk,1:jpj,1:jpi), big_fact_zcc(1:jpk,1:jpj,1:jpi) )
+  !$acc update device( zbtr_arr(1:jpk,1:jpj,1:jpi) )
 
-  !$acc enter data create( zaa(1:jpk,1:jpj,1:jpi), zbb(1:jpk,1:jpj,1:jpi), zcc(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc enter data create( inv_eu(1:jpk,1:jpj,1:jpi), inv_ev(1:jpk,1:jpj,1:jpi), inv_et(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc enter data create( big_fact_zaa (1:jpk,1:jpj,1:jpi), big_fact_zbb(1:jpk,1:jpj,1:jpi), big_fact_zcc(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc enter data create( zbtr_arr(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-
-  !$acc enter data create( e1t(1:jpj,1:jpi), e2t(1:jpj,1:jpi), e3t(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc enter data create( e1u(1:jpj,1:jpi), e2u(1:jpj,1:jpi), e3u(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc enter data create( e1v(1:jpj,1:jpi), e2v(1:jpj,1:jpi), e3v(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc enter data create( e3w(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc enter data create( un(1:jpk,1:jpj,1:jpi), vn(1:jpk,1:jpj,1:jpi), wn(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
+  !$acc update device( e1t(1:jpj,1:jpi), e2t(1:jpj,1:jpi), e3t(1:jpk,1:jpj,1:jpi) )
+  !$acc update device( e1u(1:jpj,1:jpi), e2u(1:jpj,1:jpi), e3u(1:jpk,1:jpj,1:jpi) )
+  !$acc update device( e1v(1:jpj,1:jpi), e2v(1:jpj,1:jpi), e3v(1:jpk,1:jpj,1:jpi) )
+  !$acc update device( e3w(1:jpk,1:jpj,1:jpi) )
+  !$acc update device( un(1:jpk,1:jpj,1:jpi), vn(1:jpk,1:jpj,1:jpi), wn(1:jpk,1:jpj,1:jpi) )
 
-  !$acc update device( zaa(1:jpk,1:jpj,1:jpi), zbb(1:jpk,1:jpj,1:jpi), zcc(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc update device( inv_eu(1:jpk,1:jpj,1:jpi), inv_ev(1:jpk,1:jpj,1:jpi), inv_et(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc update device( big_fact_zaa (1:jpk,1:jpj,1:jpi), big_fact_zbb(1:jpk,1:jpj,1:jpi), big_fact_zcc(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc update device( zbtr_arr(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
+  !$acc update device(tra(1:jpk,1:jpj,1:jpi,1:jptra))
+  !$acc update device(trn(1:jpk,1:jpj,1:jpi,1:jptra))
+  !$acc update device(advmask(1:jpk,1:jpj,1:jpi))
+  !$acc update device(flx_ridxt(1:Fsize,1:4))
+  !$acc update device( diaflx(1:7, 1:Fsize, 1:jptra))
 
-  !$acc update device( e1t(1:jpj,1:jpi), e2t(1:jpj,1:jpi), e3t(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc update device( e1u(1:jpj,1:jpi), e2u(1:jpj,1:jpi), e3u(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc update device( e1v(1:jpj,1:jpi), e2v(1:jpj,1:jpi), e3v(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc update device( e3w(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc update device( un(1:jpk,1:jpj,1:jpi), vn(1:jpk,1:jpj,1:jpi), wn(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-
-  !$acc kernels default(present) if(use_gpu)
+  !$acc kernels default(present) async(queue)
   DO ji = 1,jpi
      DO jj = 1,jpj
         !dir$ vector aligned
@@ -206,7 +225,7 @@ SUBROUTINE trcadv
   !$OMP TASK private(ji,jj) firstprivate(jpim1,jpjm1,jpi,jpj,jpk) default(none) &
   !$OMP shared(zdt,zaa,inv_eu,e1u,e2u,e3u,un,big_fact_zaa)
 
-  !$acc kernels default(present) if(use_gpu)
+  !$acc kernels default(present) async(queue)
   DO ji = 1,jpi
      DO jj = 1,jpj
         !dir$ vector aligned
@@ -218,7 +237,7 @@ SUBROUTINE trcadv
   !$acc end kernels
 
 
-  !$acc kernels default(present) if(use_gpu)
+  !$acc kernels default(present) async(queue)
   DO ji = 1,jpi
      DO jj = 1,jpj
         !dir$ vector aligned
@@ -229,7 +248,7 @@ SUBROUTINE trcadv
   END DO
   !$acc end kernels
 
-  !$acc kernels default(present) if(use_gpu)
+  !$acc kernels default(present) async(queue)
   DO ji = 1,jpi
      DO jj = 1,jpj
         !dir$ vector aligned
@@ -246,7 +265,7 @@ SUBROUTINE trcadv
   !$OMP TASK private(ji,jj) firstprivate(jpim1,jpjm1,jpi,jpj,jpk)  default(none) &
   !$OMP shared(inv_ev,e1v,e2v,e3v,vn,zdt,zbb,big_fact_zbb)
 
-  !$acc kernels default(present) if(use_gpu)
+  !$acc kernels default(present) async(queue)
   DO ji = 1,jpi
      DO jj = 1,jpj
         !dir$ vector aligned
@@ -257,7 +276,7 @@ SUBROUTINE trcadv
   END DO
   !$acc end kernels
 
-  !$acc kernels default(present) if(use_gpu)
+  !$acc kernels default(present) async(queue)
   DO ji = 1,jpi
      DO jj = 1,jpj
         !dir$ vector aligned
@@ -268,7 +287,7 @@ SUBROUTINE trcadv
   END DO
   !$acc end kernels
 
-  !$acc kernels default(present) if(use_gpu)
+  !$acc kernels default(present) async(queue)
   DO ji = 1,jpi
      DO jj = 1,jpj
         !dir$ vector aligned
@@ -283,7 +302,7 @@ SUBROUTINE trcadv
   !$OMP TASK private(ji,jj) firstprivate(jpim1,jpjm1,jpi,jpj,jpk) default(none) &
   !$OMP shared(inv_et,e1t,e2t,e3w,wn,zcc,zdt,big_fact_zcc)
 
-  !$acc kernels default(present) if(use_gpu)
+  !$acc kernels default(present) async(queue)
   DO ji = 1,jpi
      DO jj = 1,jpj
         !dir$ vector aligned
@@ -294,7 +313,7 @@ SUBROUTINE trcadv
   END DO
   !$acc end kernels
 
-  !$acc kernels default(present) if(use_gpu)
+  !$acc kernels default(present) async(queue)
   DO ji = 1,jpi
      DO jj = 1,jpj
         !dir$ vector aligned
@@ -305,7 +324,7 @@ SUBROUTINE trcadv
   END DO
   !$acc end kernels
 
-  !$acc kernels default(present) if(use_gpu)
+  !$acc kernels default(present) async(queue)
   DO ji = 1,jpi
      DO jj = 1,jpj
         !dir$ vector aligned
@@ -319,57 +338,10 @@ SUBROUTINE trcadv
 
   !$OMP TASKWAIT
 
-  !$acc update host( zaa(1:jpk,1:jpj,1:jpi), zbb(1:jpk,1:jpj,1:jpi), zcc(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc update host( inv_eu(1:jpk,1:jpj,1:jpi), inv_ev(1:jpk,1:jpj,1:jpi), inv_et(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc update host( big_fact_zaa (1:jpk,1:jpj,1:jpi), big_fact_zbb(1:jpk,1:jpj,1:jpi), big_fact_zcc(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc update host( zbtr_arr(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-
-
   !!     tracer loop parallelized (macrotasking)
   !!     =======================================
 
-
-
-  !!OpenMP compatibility broken. Possibility to use ifndef OpenMP + rename the file in trcadv.F90 to keep it
-  allocate(zy(jpk,jpj,jpi))
-  allocate(zx(jpk,jpj,jpi))
-  allocate(zz(jpk,jpj,jpi))
-  allocate(ztj(jpk,jpj,jpi))
-  allocate(zti(jpk,jpj,jpi))
-  allocate(zkx(jpk,jpj,jpi))
-  allocate(zky(jpk,jpj,jpi))
-  allocate(zkz(jpk,jpj,jpi))
-  allocate(zbuf(jpk,jpj,jpi))
-
-  zy(:,:,:) = 0
-  zz(:,:,:) = 0
-  zx(:,:,:) = 0
-  ztj(:,:,:)= 0
-  zti(:,:,:)= 0
-  zbuf(:,:,:) = 0.
-  zkx(:,:,:)=0.
-  zky(:,:,:)=0.
-  zkz(:,:,:)=0.
-
-  !!trn could be allocate earlier
-  !$acc enter data create(trn(1:jpk,1:jpj,1:jpi,1:jptra)) if(use_gpu)
-  !$acc enter data create(tra(1:jpk,1:jpj,1:jpi,1:jptra)) if(use_gpu)
-  !$acc enter data create(advmask(1:jpk,1:jpj,1:jpi)) if(use_gpu)
-  !$acc enter data create(flx_ridxt(1:Fsize,1:4)) if(use_gpu)
-  !$acc enter data create( diaflx(1:7, 1:Fsize, 1:jptra)) if(use_gpu)
-
-  !$acc enter data create( zy(1:jpk,1:jpj,1:jpi), zx(1:jpk,1:jpj,1:jpi), zz(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc enter data create( ztj(1:jpk,1:jpj,1:jpi), zti(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc enter data create( zkx(1:jpk,1:jpj,1:jpi), zky(1:jpk,1:jpj,1:jpi), zkz(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc enter data create( zbuf(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-
-  !$acc update device(tra(1:jpk,1:jpj,1:jpi,1:jptra)) if(use_gpu)
-  !$acc update device(trn(1:jpk,1:jpj,1:jpi,1:jptra)) if(use_gpu)
-  !$acc update device(advmask(1:jpk,1:jpj,1:jpi)) if(use_gpu)
-  !$acc update device(flx_ridxt(1:Fsize,1:4)) if(use_gpu)
-  !$acc update device( diaflx(1:7, 1:Fsize, 1:jptra)) if(use_gpu)
-
-  !$acc kernels default(present) if(use_gpu)
+  !$acc kernels default(present) async(queue)
   DO ji = 1, jpi
      DO jj = 1, jpj
         DO jk = 1, jpk
@@ -386,6 +358,7 @@ SUBROUTINE trcadv
      ENDDO
   ENDDO
   !$acc end kernels
+  !$acc wait(queue)
 
   !$omp taskloop default(none) private(jf,junk,junki,junkj,junkk,zbtr) &
   !$omp private(zkx,zky,zkz,zti,ztj,zx,zy,zz,zbuf) shared(diaflx,jarrt,tra,zdt) &
@@ -404,8 +377,8 @@ SUBROUTINE trcadv
      !!           and mass fluxes calculated above
      !!       calcul of tracer flux in the i and j direction
 
-     !$acc kernels async default(present) if(use_gpu)
-     !$acc loop independent
+
+     !$acc kernels default(present) async(queue)
      DO ji = 2,jpim1
         !dir$ vector aligned
         DO jj = 2,jpjm1
@@ -414,7 +387,7 @@ SUBROUTINE trcadv
      END DO
      !$acc end kernels
 
-     !$acc kernels async default(present) if(use_gpu)
+     !$acc kernels default(present) async(queue)
      DO ji = 2,jpim1
         !dir$ vector aligned
         !$acc loop independent
@@ -424,7 +397,7 @@ SUBROUTINE trcadv
      END DO
      !$acc end kernels
 
-     !$acc kernels async default(present) if(use_gpu)
+     !$acc kernels default(present) async(queue)
      DO ji = 1,jpi
         !dir$ vector aligned
         !$acc loop independent
@@ -434,7 +407,7 @@ SUBROUTINE trcadv
      ENDDO
      !$acc end kernels
 
-     !$acc kernels async default(present) if(use_gpu)
+     !$acc kernels default(present) async(queue)
      DO ji = 1,jpi
         !dir$ vector aligned
         !$acc loop independent
@@ -444,7 +417,7 @@ SUBROUTINE trcadv
      END DO
      !$acc end kernels
      ! loop unfusion
-     !$acc kernels async default(present) if(use_gpu)
+     !$acc kernels default(present) async(queue)
      DO jj = 2,jpjm1
         !dir$ vector aligned
         !$acc loop independent
@@ -454,7 +427,7 @@ SUBROUTINE trcadv
      END DO
      !$acc end kernels
 
-     !$acc kernels async default(present) if(use_gpu)
+     !$acc kernels default(present) async(queue)
      DO jj = 2,jpjm1
         !dir$ vector aligned
         !$acc loop independent
@@ -464,7 +437,7 @@ SUBROUTINE trcadv
      END DO
      !$acc end kernels
 
-     !$acc kernels async default(present) if(use_gpu)
+     !$acc kernels default(present) async(queue)
      !$acc loop independent
      DO  ji = 2,jpim1
         DO jj = 2,jpjm1
@@ -476,7 +449,7 @@ SUBROUTINE trcadv
      END DO
      !$acc end kernels
 
-     !$acc kernels async default(present) if(use_gpu)
+     !$acc kernels default(present) async(queue)
      DO  ji = 2,jpim1
         !$acc loop independent
         DO jj = 2,jpjm1
@@ -488,7 +461,7 @@ SUBROUTINE trcadv
      END DO
      !$acc end kernels
 
-     !$acc parallel loop async collapse(3) gang vector default(present) if(use_gpu)
+     !$acc parallel loop collapse(3) gang vector default(present) async(queue)
      DO  ji = 2,jpim1
         DO jj = 2,jpjm1
            !dir$ vector aligned
@@ -498,21 +471,15 @@ SUBROUTINE trcadv
         END DO
      END DO
      !$acc end parallel loop
-     !$acc wait
+     !$acc wait(queue)
 
      ! ... Lateral boundary conditions on zk[xy]
 #ifdef key_mpp
 
      !  ... Mpp : export boundary values to neighboring processors
 
-#ifndef _OPENACC
-     CALL mpplnk_my(zkx)
-     CALL mpplnk_my(zky)
-#else
-     CALL mpplnk_my_openacc(zkx,gpu=use_gpu)
-     CALL mpplnk_my_openacc(zky,gpu=use_gpu)
-
-#endif
+     CALL mpplnk_my(zkx,gpu=use_gpu)
+     CALL mpplnk_my(zky,gpu=use_gpu)
 
 #else
 
@@ -526,7 +493,7 @@ SUBROUTINE trcadv
      !! 2. calcul of after field using an upstream advection scheme
      !! -----------------------------------------------------------
 
-     !$acc kernels default(present) if(use_gpu)
+     !$acc kernels default(present) async(queue)
      DO ji =2,jpim1
         DO jj =2,jpjm1
            DO jk =1,jpkm1
@@ -539,7 +506,7 @@ SUBROUTINE trcadv
      ENDDO
      !$acc end kernels
 
-     !$acc kernels default(present) if(use_gpu)
+     !$acc kernels default(present) async(queue)
      DO jf=1,Fsize
         jk = flx_ridxt(jf,2)
         jj = flx_ridxt(jf,3)
@@ -565,7 +532,7 @@ SUBROUTINE trcadv
         if(jt .EQ. 1) then
 
            if(ncor .EQ. 1) then
-              !$acc kernels default(present) if(use_gpu)
+              !$acc kernels default(present)
               DO ji = 2,jpim1
                  DO jj = 2,jpjm1
                     !dir$ vector aligned
@@ -577,7 +544,7 @@ SUBROUTINE trcadv
               !$acc end kernels
 
            else
-              !$acc kernels default(present) if(use_gpu)
+              !$acc kernels default(present) async(queue)
               DO ji = 2,jpim1
                  DO jj = 2,jpjm1
                     !dir$ vector aligned
@@ -589,7 +556,7 @@ SUBROUTINE trcadv
               END DO
               !$acc end kernels
 
-              !$acc kernels default(present) if(use_gpu)
+              !$acc kernels default(present) async(queue)
               DO ji = 2,jpim1
                  DO jj = 2,jpjm1
                     !dir$ vector aligned
@@ -604,7 +571,7 @@ SUBROUTINE trcadv
            endif
 
         else
-           !$acc kernels default(present) if(use_gpu)
+           !$acc kernels default(present) async(queue)
            DO ji = 2,jpim1
               DO jj = 2,jpjm1
                  !dir$ vector aligned
@@ -615,7 +582,7 @@ SUBROUTINE trcadv
            END DO
            !$acc end kernels
 
-           !$acc kernels default(present) if(use_gpu)
+           !$acc kernels default(present) async(queue)
            DO ji = 2,jpim1
               DO jj = 2,jpjm1
                  !dir$ vector aligned
@@ -627,17 +594,12 @@ SUBROUTINE trcadv
            !$acc end kernels
         endif
 
-
-
+        !$acc wait(queue)
 
         !! ... Lateral boundary conditions on zti
 #ifdef key_mpp
         ! ... Mpp : export boundary values to neighboring processors
-#ifndef _OPENACC
-        CALL mpplnk_my(zti)
-#else
-        CALL mpplnk_my_openacc(zti,gpu=use_gpu)
-#endif
+        CALL mpplnk_my(zti,gpu=use_gpu)
 #else
         ! ... T-point, 3D array, full local array zti is initialised
         CALL lbc( zti(:,:,:), 1, 1, 1, 1, jpk, 1, gpu=use_gpu )
@@ -648,7 +610,7 @@ SUBROUTINE trcadv
 
         !jk = 1
         !          DO jk = 1,jpkm1
-        !$acc kernels default(present) if(use_gpu)
+        !$acc kernels default(present) async(queue)
         DO ji = 2,jpim1
            DO jj = 2,jpjm1
               junk  = zti(1,jj,ji )
@@ -661,7 +623,7 @@ SUBROUTINE trcadv
         !$acc end kernels
 
         !DO ju=1, dimen_jarr2
-        !$acc kernels default(present) if(use_gpu)
+        !$acc kernels default(present) async(queue)
         DO ji = 2,jpim1
            DO jj = 2,jpjm1
               !dir$ vector aligned
@@ -687,20 +649,15 @@ SUBROUTINE trcadv
         !$acc end kernels
         !                 endif
 
+        !$acc wait(queue)
 
         ! ... Lateral boundary conditions on z[xyz]
 #ifdef key_mpp
 
         ! ... Mpp : export boundary values to neighboring processors
-#ifndef _OPENACC
-        CALL mpplnk_my(zx)
-        CALL mpplnk_my(zy)
-        CALL mpplnk_my(zz)
-#else
-        CALL mpplnk_my_openacc(zx,gpu=use_gpu)
-        CALL mpplnk_my_openacc(zy,gpu=use_gpu)
-        CALL mpplnk_my_openacc(zz,gpu=use_gpu)
-#endif
+        CALL mpplnk_my(zx,gpu=use_gpu)
+        CALL mpplnk_my(zy,gpu=use_gpu)
+        CALL mpplnk_my(zz,gpu=use_gpu)
 #else
 
         !  ... T-point, 3D array, full local array z[xyz] are initialised
@@ -713,7 +670,7 @@ SUBROUTINE trcadv
         !!            2.5 calcul of the final field:
         !!                advection by antidiffusive mass fluxes and an upstream scheme
 
-        !$acc kernels default(present) if(use_gpu)
+        !$acc kernels default(present) async(queue)
         !$acc loop independent
         DO ji = 2,jpim1
            !dir$ vector aligned
@@ -723,7 +680,7 @@ SUBROUTINE trcadv
         END DO
         !$acc end kernels
 
-        !$acc kernels default(present) if(use_gpu)
+        !$acc kernels default(present) async(queue)
         DO ji = 2,jpim1
            !dir$ vector aligned
            !$acc loop independent
@@ -733,7 +690,7 @@ SUBROUTINE trcadv
         END DO
         !$acc end kernels
 
-        !$acc kernels default(present) if(use_gpu)
+        !$acc kernels default(present) async(queue)
         DO ji = 1,jpi
            !dir$ vector aligned
            !$acc loop independent
@@ -743,7 +700,7 @@ SUBROUTINE trcadv
         ENDDO
         !$acc end kernels
 
-        !$acc kernels default(present) if(use_gpu)
+        !$acc kernels default(present) async(queue)
         DO ji = 1,jpi
            !dir$ vector aligned
            !$acc loop independent
@@ -753,7 +710,7 @@ SUBROUTINE trcadv
         ENDDO
         !$acc end kernels
 
-        !$acc kernels default(present) if(use_gpu)
+        !$acc kernels default(present) async(queue)
         DO jj = 2,jpjm1
            !dir$ vector aligned
            !$acc loop independent
@@ -763,7 +720,7 @@ SUBROUTINE trcadv
         ENDDO
         !$acc end kernels
 
-        !$acc kernels default(present) if(use_gpu)
+        !$acc kernels default(present) async(queue)
         DO jj = 2,jpjm1
            !dir$ vector aligned
            !$acc loop independent
@@ -773,7 +730,7 @@ SUBROUTINE trcadv
         END DO
         !$acc end kernels
 
-        !$acc kernels default(present) if(use_gpu)
+        !$acc kernels default(present) async(queue)
         !$acc loop independent
         DO  ji = 2,jpim1
            DO jj = 2,jpjm1
@@ -787,7 +744,7 @@ SUBROUTINE trcadv
         END DO
         !$acc end kernels
 
-        !$acc kernels default(present) if(use_gpu)
+        !$acc kernels default(present) async(queue)
         DO  ji = 2,jpim1
            !$acc loop independent
            DO jj = 2,jpjm1
@@ -801,7 +758,7 @@ SUBROUTINE trcadv
         END DO
         !$acc end kernels
 
-        !$acc kernels default(present) if(use_gpu)
+        !$acc kernels default(present) async(queue)
         DO  ji = 2,jpim1
            DO jj = 2,jpjm1
               !dir$ vector aligned
@@ -815,17 +772,13 @@ SUBROUTINE trcadv
         END DO
         !$acc end kernels
 
+        !$acc wait(queue)
 
         !... Lateral boundary conditions on zk[xy]
 #ifdef key_mpp
         !  ... Mpp : export boundary values to neighboring processors
-#ifndef _OPENACC
-        CALL mpplnk_my(zkx)
-        CALL mpplnk_my(zky)
-#else
-        CALL mpplnk_my_openacc(zkx,gpu=use_gpu)
-        CALL mpplnk_my_openacc(zky,gpu=use_gpu)
-#endif
+        CALL mpplnk_my(zkx,gpu=use_gpu)
+        CALL mpplnk_my(zky,gpu=use_gpu)
 #else
         ! ... T-point, 3D array, full local array zk[xy] are initialised
         CALL lbc( zkx(:,:,:), 1, 1, 1, 1, jpk, 1, gpu=use_gpu )
@@ -837,7 +790,7 @@ SUBROUTINE trcadv
 
 
         if(ncor .EQ. 1) then
-           !$acc kernels default(present) if(use_gpu)
+           !$acc kernels default(present) async(queue)
            DO ji =2,jpim1
               DO jj =2,jpjm1
                  DO jk =1,jpkm1
@@ -848,7 +801,7 @@ SUBROUTINE trcadv
            ENDDO
            !$acc end kernels
 
-           !$acc kernels default(present) if(use_gpu)
+           !$acc kernels default(present) async(queue)
            DO jf=1,Fsize
               jk = flx_ridxt(jf,2)
               jj = flx_ridxt(jf,3)
@@ -863,7 +816,7 @@ SUBROUTINE trcadv
 
         else
 
-           !$acc kernels default(present) if(use_gpu)
+           !$acc kernels default(present) async(queue)
            DO ji =2,jpim1
               DO jj =2,jpjm1
                  DO jk =1,jpkm1
@@ -874,7 +827,7 @@ SUBROUTINE trcadv
            ENDDO
            !$acc end kernels
 
-           !$acc kernels default(present) if(use_gpu)
+           !$acc kernels default(present) async(queue)
            DO jf=1,Fsize
               jk = flx_ridxt(jf,2)
               jj = flx_ridxt(jf,3)
@@ -886,7 +839,6 @@ SUBROUTINE trcadv
            ENDDO
            !$acc end kernels
 
-
         endif
 
      ENDDO ANTIDIFF_CORR
@@ -898,7 +850,7 @@ SUBROUTINE trcadv
 
      if(ncor .EQ. 1) then
 
-        !$acc kernels default(present) if(use_gpu)
+        !$acc kernels default(present) async(queue)
         do ji=1,jpi
            do jj=1,jpj
               do jk=1,jpk
@@ -910,7 +862,7 @@ SUBROUTINE trcadv
 
      else
 
-        !$acc kernels default(present) if(use_gpu)
+        !$acc kernels default(present) async(queue)
         do ji=1,jpi
            do jj=1,jpj
               do jk=1,jpk
@@ -936,21 +888,24 @@ SUBROUTINE trcadv
 
 
   END DO TRACER_LOOP
+  !$acc wait(queue)
+
   !$OMP end taskloop
 
-  !$acc update host( diaflx(1:7, 1:Fsize, 1:jptra) ) if(use_gpu)
-  !$acc update host( tra(1:jpk,1:jpj,1:jpi,1:jptra) ) if(use_gpu)
+  !$acc update host( zaa(1:jpk,1:jpj,1:jpi), zbb(1:jpk,1:jpj,1:jpi), zcc(1:jpk,1:jpj,1:jpi) )
+  !$acc update host( inv_eu(1:jpk,1:jpj,1:jpi), inv_ev(1:jpk,1:jpj,1:jpi), inv_et(1:jpk,1:jpj,1:jpi) )
+  !$acc update host( big_fact_zaa (1:jpk,1:jpj,1:jpi), big_fact_zbb(1:jpk,1:jpj,1:jpi), big_fact_zcc(1:jpk,1:jpj,1:jpi) )
+  !$acc update host( zbtr_arr(1:jpk,1:jpj,1:jpi) )
 
-  !$acc update host( zy(1:jpk,1:jpj,1:jpi), zx(1:jpk,1:jpj,1:jpi), zz(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc update host( ztj(1:jpk,1:jpj,1:jpi), zti(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc update host( zkx(1:jpk,1:jpj,1:jpi), zky(1:jpk,1:jpj,1:jpi), zkz(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
-  !$acc update host( zbuf(1:jpk,1:jpj,1:jpi) ) if(use_gpu)
+  !$acc update host( diaflx(1:7, 1:Fsize, 1:jptra) )
+  !$acc update host( tra(1:jpk,1:jpj,1:jpi,1:jptra) )
 
-  !$acc exit data delete( tra) finalize if(use_gpu)
-  !$acc exit data delete( trn, advmask ) finalize if(use_gpu)
-  !$acc exit data delete( flx_ridxt, diaflx ) finalize if(use_gpu)
-  !$acc exit data delete( zy, zx, zz, ztj, zti, zkx, zky, zkz, zbuf ) finalize if(use_gpu)
+  !$acc update host( zy(1:jpk,1:jpj,1:jpi), zx(1:jpk,1:jpj,1:jpi), zz(1:jpk,1:jpj,1:jpi) )
+  !$acc update host( ztj(1:jpk,1:jpj,1:jpi), zti(1:jpk,1:jpj,1:jpi) )
+  !$acc update host( zkx(1:jpk,1:jpj,1:jpi), zky(1:jpk,1:jpj,1:jpi), zkz(1:jpk,1:jpj,1:jpi) )
+  !$acc update host( zbuf(1:jpk,1:jpj,1:jpi) )
 
+  !$acc exit data delete( zy, zx, zz, ztj, zti, zkx, zky, zkz, zbuf ) finalize
   !!OpenMP compatibility broken. Possibility to use ifndef OpenMP + rename the file in trcadv.F90 to keep it
   deallocate(zy )
   deallocate(zx )
@@ -962,24 +917,22 @@ SUBROUTINE trcadv
   deallocate(zkz )
   deallocate(zbuf )
 
-  !$acc exit data delete( zaa, zbb, zcc, inv_eu, inv_ev, inv_et, big_fact_zaa , big_fact_zbb, big_fact_zcc, zbtr_arr ) finalize if(use_gpu)
-  !$acc exit data delete( e1t, e2t, e3t, e1u, e2u, e3u, e1v, e2v, e3v, e3w, un, vn, wn ) finalize if(use_gpu)
-
   trcadvparttime = MPI_WTIME() - trcadvparttime
   trcadvtottime = trcadvtottime + trcadvparttime
 !!!!
 
 contains
 
+#ifndef _OPENACC
   double precision function fsx(pfx1, pfx2, pfu)
     !dir$ attributes vector :: fsx
-    !$acc routine seq
     IMPLICIT NONE
     double precision, INTENT(IN) :: pfx1, pfx2, pfu
     double precision :: abspfu
     abspfu = abs(pfu)
     fsx = ((pfu + abspfu) * pfx1 + (pfu - abspfu) * pfx2) * 0.5
   end function fsx
+#endif
 
 
 END SUBROUTINE trcadv
diff --git a/src/PHYS/trcave.f90 b/src/PHYS/trcave.f90
index 8d8313a5..617a0212 100644
--- a/src/PHYS/trcave.f90
+++ b/src/PHYS/trcave.f90
@@ -10,14 +10,19 @@ SUBROUTINE trcave
       integer :: jn_high, jn_on_all
       double precision ::  Miss_val =1.e20
       double precision :: elapsed_time, inv_incremented_time
+      integer :: queue
 
       ave_partTime = MPI_WTIME()
 
+      queue=1
 
 !     FIRST, LOW FREQUENCY
       elapsed_time         =     elapsed_time_2
       inv_incremented_time = 1./(elapsed_time_2 + rdt)
 
+      !$acc update device(traIO,trn,umask,vmask,tmask,traIO_HIGH,highfreq_table,snIO,tnIO,wnIO,avtIO,e3tIO,unIO,vnIO,sn,tn,wn,avt,e3t,un,vn,tra_DIA_IO,tra_DIA,tra_DIA_2d_IO,tra_DIA_2d,vatmIO,empIO,qsrIO,vatm,emp,qsr,highfreq_table_dia,tra_DIA_IO_HIGH,tra_DIA_2d_IO_HIGH,highfreq_table_dia2d)
+
+      !$acc parallel loop gang vector collapse(4) default(present) async(queue)
       DO jn=1 ,jptra
 
           DO ji=1, jpi
@@ -33,6 +38,7 @@ SUBROUTINE trcave
           END DO
 
       END DO
+      !$acc end parallel loop
 
        
 
@@ -40,15 +46,21 @@ SUBROUTINE trcave
 
       elapsed_time         =     elapsed_time_1
       inv_incremented_time = 1./(elapsed_time_1 + rdt)! ****************** HIGH FREQUENCY
+      !$acc parallel loop gang vector collapse(4) default(present) async(queue)
       DO jn_high=1 ,jptra_high
 
 
 
        
+#ifndef _OPENACC
           jn_on_all = highfreq_table(jn_high)
+#endif
           DO ji=1, jpi
           DO jj=1, jpj
           DO jk=1, jpk
+#ifdef _OPENACC
+               jn_on_all = highfreq_table(jn_high)
+#endif
                IF(tmask(jk,jj,ji) .NE. 0.) THEN
                   traIO_HIGH(jk,jj,ji,jn_high  )= &
      &           (traIO_HIGH(jk,jj,ji,jn_high )*elapsed_time+trn(jk,jj,ji,jn_on_all)*rdt)*inv_incremented_time
@@ -61,6 +73,7 @@ SUBROUTINE trcave
    
 
       END DO
+      !$acc end parallel loop
 
 
 !     *****************  PHYS *****************************************************
@@ -73,6 +86,7 @@ SUBROUTINE trcave
       endif
 
 
+          !$acc parallel loop gang vector collapse(3) default(present) async(queue)
           DO ji=1, jpi
        DO jj=1, jpj
       DO jk=1, jpk
@@ -107,7 +121,9 @@ SUBROUTINE trcave
           END DO
        END DO
       END DO
+      !$acc end parallel loop
 
+      !$acc parallel loop gang vector collapse(2) default(present) async(queue)
       DO jj=1, jpj
         DO ji=1, jpi
            IF (tmask(1,jj,ji) .NE. 0.) THEN
@@ -121,6 +137,7 @@ SUBROUTINE trcave
            ENDIF
         END DO
       END DO
+      !$acc end parallel loop
 
 
 !     *****************  END PHYS *************************************************
@@ -135,6 +152,7 @@ SUBROUTINE trcave
       inv_incremented_time = 1./(elapsed_time_2 + rdt)
 
 
+         !$acc parallel loop gang vector collapse(4) default(present) async(queue)
          DO jn = 1,jptra_dia
          DO ji=1, jpi
          DO jj=1, jpj
@@ -148,10 +166,12 @@ SUBROUTINE trcave
          END DO
          END DO
          ENDDO
+         !$acc end parallel loop
 
 !     *********************  DIAGNOSTICS 2D **********
 
 
+          !$acc parallel loop gang vector collapse(2) default(present) async(queue)
           DO ji=1, jpi
           DO jj=1, jpj
               IF(tmask(1,jj,ji) .NE. 0.) THEN ! Warning ! Tested only for surface
@@ -162,6 +182,7 @@ SUBROUTINE trcave
                   ENDIF
           END DO
           END DO
+          !$acc end parallel loop
 
 
 
@@ -172,11 +193,17 @@ SUBROUTINE trcave
 
 
       if (jptra_dia_high.gt.0) THEN
+         !$acc parallel loop gang vector collapse(4) default(present) async(queue)
          DO jn_high=1, jptra_dia_high
+#ifndef _OPENACC
          jn_on_all = highfreq_table_dia(jn_high )
+#endif
          DO ji=1, jpi
          DO jj=1, jpj
          DO jk=1, jpk
+#ifdef _OPENACC
+         jn_on_all = highfreq_table_dia(jn_high )
+#endif
          IF(tmask(jk,jj,ji) .NE. 0.) THEN
                tra_DIA_IO_HIGH(jk,jj,ji,jn_high )= &
      &         (tra_DIA_IO_HIGH(jk,jj,ji,jn_high )*elapsed_time+tra_DIA(jk,jj,ji,jn_on_all)*rdt)*inv_incremented_time
@@ -187,6 +214,7 @@ SUBROUTINE trcave
          END DO
          END DO
          END DO
+         !$acc end parallel loop
       endif
 
 
@@ -194,24 +222,29 @@ SUBROUTINE trcave
 
 
        if (jptra_dia2d_high.gt.0) THEN
+             !$acc parallel loop gang vector collapse(3) default(present) async(queue)
              DO ji=1, jpi
              DO jj=1, jpj
-                IF(tmask(1,jj,ji) .NE. 0.) THEN
                 DO jn_high=1, jptra_dia2d_high
-                   jn_on_all = highfreq_table_dia2d(jn_high)
-                   tra_DIA_2d_IO_HIGH(jn_high,jj,ji)= &
-     &            (tra_DIA_2d_IO_HIGH(jn_high,jj,ji)*elapsed_time+tra_DIA_2d(jn_on_all,jj,ji)*rdt)*inv_incremented_time
+                  IF(tmask(1,jj,ji) .NE. 0.) THEN
+                     jn_on_all = highfreq_table_dia2d(jn_high)
+                       tra_DIA_2d_IO_HIGH(jn_high,jj,ji)= &
+         &            (tra_DIA_2d_IO_HIGH(jn_high,jj,ji)*elapsed_time+tra_DIA_2d(jn_on_all,jj,ji)*rdt)*inv_incremented_time
+                  ELSE
+                     tra_DIA_2d_IO_HIGH(jn_high,jj,ji)=Miss_val
+                  ENDIF
                 END DO
-                ELSE
-                   tra_DIA_2d_IO_HIGH(:,jj,ji)=Miss_val
-                ENDIF
              END DO
              END DO
+             !$acc end parallel loop
         endif
 
 
       endif ! lfbm
 
+      !$acc wait(queue)
+      !$acc update host(traIO,traIO_HIGH,snIO,tnIO,wnIO,avtIO,e3tIO,unIO,vnIO,vatmIO,empIO,qsrIO,tra_DIA_IO,tra_DIA_2d_IO,tra_DIA_2d,tra_DIA_IO_HIGH,tra_DIA_2d_IO_HIGH)
+
       ave_partTime = MPI_WTIME() - ave_partTime
       ave_TotTime = ave_TotTime  + ave_partTime
 
diff --git a/src/PHYS/trchdf.f90 b/src/PHYS/trchdf.f90
index a70bfc9d..f43b8198 100644
--- a/src/PHYS/trchdf.f90
+++ b/src/PHYS/trchdf.f90
@@ -109,10 +109,19 @@ SUBROUTINE trchdf
       INTEGER :: locsum,jklef,jjlef,jilef,jkrig,jjrig,jirig
       !INTEGER, allocatable :: jarr_hdf(:,:,:),jarr_hdf_flx(:)
       double precision, allocatable,dimension(:,:,:) :: zlt, ztu, ztv
+      integer :: queue
+      logical :: use_gpu
 !!----------------------------------------------------------------------
 !! statement functions
 !! ===================
 
+
+#ifdef _OPENACC
+      use_gpu=.true.
+#else
+      use_gpu=.false.
+#endif
+
     !  #include "BFM_var_list.h"
        trcbilaphdfparttime = MPI_WTIME()
 
@@ -132,8 +141,9 @@ SUBROUTINE trchdf
              allocate(zeev   (jpk,jpj,jpi      )) 
              zeev         = huge(zeev(1,1,1))
              allocate(zbtr   (jpk,jpj,jpi      )) 
-             zbtr         = huge(zbtr(1,1,1)) 
+             zbtr         = huge(zbtr(1,1,1))
 
+             !$acc enter data create(hdfmask,zeeu,zeev,zbtr)
 
                 DO ji = 1,jpi
              DO jj = 1,jpj
@@ -168,10 +178,20 @@ SUBROUTINE trchdf
              END DO
           END DO
 
+        queue=1
+        !$acc update device(hdfmask)
         hdf_initialized=.true.
+
         ENDIF
+
+        allocate(zlt    (jpk,jpj,jpi))
+        allocate(ztu    (jpk,jpj,jpi))
+        allocate(ztv    (jpk,jpj,jpi))
+        !$acc enter data create(zlt,ztu,ztv)
+
 !  Metric arrays calculated out of the initialisation phase(for z- or s-coordinates)
 ! !! ----------------------------------
+                !$acc parallel loop gang vector collapse(3) default(present) async(queue)
                 DO ji = 1, jpi
               DO jj = 1, jpj
            DO jk=1,jpk
@@ -182,13 +202,11 @@ SUBROUTINE trchdf
                 END DO
               END DO
            END DO
+           !$acc end parallel loop
 
 
 
      
-             allocate(zlt    (jpk,jpj,jpi)) 
-             allocate(ztu    (jpk,jpj,jpi)) 
-             allocate(ztv    (jpk,jpj,jpi)) 
 !! tracer slab
 !! =============
 
@@ -199,9 +217,11 @@ SUBROUTINE trchdf
       TRACER_LOOP: DO  jn = 1, jptra
 
 
+             !$acc kernels default(present) async(queue)
              zlt = 0.
              ztu = 0.
              ztv = 0.
+             !$acc end kernels
             
 !! 1. Laplacian
 !! ------------
@@ -215,6 +235,7 @@ SUBROUTINE trchdf
             !  jk = jarr_hdf(1,jv,1a)
             
             ! $OMP TASK default(shared) private(ji,jj,jk)
+                  !$acc parallel loop gang vector collapse(3) default(present) async(queue)
                   DO ji = 1,jpi-1
               DO jj = 1,jpj-1
            DO jk = 1,jpk
@@ -228,8 +249,10 @@ SUBROUTINE trchdf
             END DO
           END DO
           ! $OMP END TASK
+          !$acc end parallel loop
 
                 ! $OMP TASK default(shared) private(ji,jj,jk)
+                !$acc parallel loop gang vector collapse(3) default(present) async(queue)
                 DO ji = 1,jpi-1
               DO jj = 1,jpj-1
            DO jk = 1,jpk
@@ -242,6 +265,7 @@ SUBROUTINE trchdf
              END DO
             END DO
           END DO
+          !$acc end parallel loop
           ! $OMP END TASK
 
           ! $OMP TASKWAIT
@@ -254,6 +278,7 @@ SUBROUTINE trchdf
       !        ji = jarr_hdf(3,jv,2)
       !        jj = jarr_hdf(2,jv,2)
       !        jk = jarr_hdf(1,jv,2)
+                  !$acc parallel loop gang vector collapse(3) default(present) async(queue)
                   DO ji = 2,jpi-1
               DO jj = 2,jpj-1
            DO jk = 1,jpk
@@ -266,6 +291,7 @@ SUBROUTINE trchdf
              END DO
             END DO
           END DO
+          !$acc end parallel loop
 
 
 
@@ -274,10 +300,10 @@ SUBROUTINE trchdf
 !! ... Lateral boundary conditions on the laplacian (zlt,zls)
 
 
+        !$acc wait(queue)
 #ifdef key_mpp
 !  ... Mpp : export boundary values to neighboring processors
-
-        CALL mpplnk_my(zlt)
+        CALL mpplnk_my(zlt,gpu=use_gpu)
 #else
                CALL lbc( zlt(:,:,:), 1, 1, 1, 1, jpk, 1 )
 #endif
@@ -291,6 +317,7 @@ SUBROUTINE trchdf
 !!!&omp&                               dimen_jvhdf3,zta,zbtr,tra,jarr_hdf_flx,diaflx,Fsize)
 
             ! $OMP TASK default(shared) private(ji,jj,jk)
+                  !$acc parallel loop gang vector collapse(3) default(present) async(queue)
                   DO ji = 1,jpi-1
               DO jj = 1,jpj-1
            DO jk = 1,jpk
@@ -301,9 +328,11 @@ SUBROUTINE trchdf
               END DO
             END DO
           END DO
+          !$acc end parallel loop
           ! $OMP END TASK
 
           ! $OMP TASK default(shared) private(ji,jj,jk)
+               !$acc parallel loop gang vector collapse(3) default(present) async(queue)
                DO ji = 1,jpi-1
               DO jj = 1,jpj-1
            DO jk = 1,jpk
@@ -315,6 +344,7 @@ SUBROUTINE trchdf
               END DO
             END DO
           END DO
+          !$acc end parallel loop
           ! $OMP END TASK
 
 !! ... fourth derivative (divergence) and add to the general tracer trend
@@ -326,7 +356,8 @@ SUBROUTINE trchdf
       !        jk = jarr_hdf(1,jv,2)
       !        jf = jarr_hdf_flx(jv)
           
-      ! $OMP TASKWAIT 
+      ! $OMP TASKWAIT
+                  !$acc parallel loop gang vector collapse(3) default(present) async(queue)
                   DO ji = 2,jpim1
               DO jj = 2,jpjm1
            DO jk = 1,jpk
@@ -339,9 +370,10 @@ SUBROUTINE trchdf
               !tra(jk,jj,ji,jn ) = tra(jk,jj,ji,jn ) + zta
                END DO
             END DO
-          END DO
-
+         END DO
+         !$acc end parallel loop
 
+         !$acc parallel loop gang vector default(present) async(queue)
          DO jf=1,Fsize
              jk = flx_ridxt(jf,2)
              jj = flx_ridxt(jf,3)
@@ -351,6 +383,7 @@ SUBROUTINE trchdf
              diaflx(6,jf, jn) = diaflx(6,jf, jn) - ztv(jk,jj,ji)*rdt
 
           ENDDO
+          !$acc end parallel loop
 
  
       
@@ -359,16 +392,18 @@ SUBROUTINE trchdf
 
         END DO TRACER_LOOP
         ! $OMP END TASKLOOP
+        !$acc wait(queue)
 
         ! deallocate(hdfmask)
         ! deallocate(zeeu)
         ! deallocate(zeev)
         ! deallocate(zbtr)
 
-             deallocate(zlt) 
-             deallocate(ztu) 
-             deallocate(ztv) 
-            
+        !$acc exit data delete(zlt,ztu,ztv)
+        deallocate(zlt)
+        deallocate(ztu)
+        deallocate(ztv)
+
 
        trcbilaphdfparttime = MPI_WTIME() - trcbilaphdfparttime
        trcbilaphdftottime = trcbilaphdftottime + trcbilaphdfparttime
diff --git a/src/PHYS/trczdf.f90 b/src/PHYS/trczdf.f90
index 0a02d892..615883f7 100644
--- a/src/PHYS/trczdf.f90
+++ b/src/PHYS/trczdf.f90
@@ -69,7 +69,7 @@ SUBROUTINE trczdf
 
 
       LOGICAL :: l1,l2,l3
-      INTEGER :: jk,jj,ji, jn, jv, jf
+      INTEGER :: jk,jj,ji, jn, jv, jf, ntx
 ! omp variables
       
 
@@ -113,8 +113,12 @@ SUBROUTINE trczdf
                END DO
             END DO
          END DO
-
+         !$acc update device(jarr_zdf,jarr_zdf_flx)
+#ifdef _OPENACC
+         call myalloc_ZDF_gpu()
+#endif
       ENDIF
+      !$acc enter data create(delta_tra,int_tra)
 
 
 !! passive tracer slab
@@ -129,12 +133,21 @@ SUBROUTINE trczdf
         ztavg = 0.e0
 !! vertical slab
 
+        ! NOTE: kernel is too big, should be split by adding a new jv dimension
+        ! on zwi zws zwd zwy zwt zwz zwx
+        !$acc parallel loop gang vector default(present) async vector_length(32)
         DO jv = 1, dimen_jvzdf
 
            ji  = jarr_zdf(2,jv)
            jj  = jarr_zdf(1,jv)
            Aij = e1t(jj,ji) * e2t(jj,ji)
 
+#ifdef _OPENACC
+           ntx=jv
+#else
+           ntx=1
+#endif
+
 !! I. Vertical trends associated with lateral mixing
 !! -------------------------------------------------
 !!    (excluding the vertical flux proportional to dk[t] )
@@ -155,14 +168,14 @@ SUBROUTINE trczdf
 !!   ... Euler time stepping when starting from rest
         DO jk = 1, jpkm1
           z2dtt = zdt * rdt
-          zwi(jk, 1) = - z2dtt * avt(jk,jj,ji  )/( e3t(jk,jj,ji) * e3w(jk,jj,ji  ) )
-          zws(jk, 1) = - z2dtt * avt(jk+1,jj,ji)/( e3t(jk,jj,ji) * e3w(jk+1,jj,ji) )
-          zwd(jk, 1) = 1. - zwi(jk, 1) - zws(jk, 1)
+          zwi(jk, ntx) = - z2dtt * avt(jk,jj,ji  )/( e3t(jk,jj,ji) * e3w(jk,jj,ji  ) )
+          zws(jk, ntx) = - z2dtt * avt(jk+1,jj,ji)/( e3t(jk,jj,ji) * e3w(jk+1,jj,ji) )
+          zwd(jk, ntx) = 1. - zwi(jk, ntx) - zws(jk, ntx)
         END DO
 
 !! Surface boundary conditions
-          zwi(1,1) = 0.e0
-          zwd(1,1) = 1. - zws(1,1)
+          zwi(1,ntx) = 0.e0
+          zwd(1,ntx) = 1. - zws(1,ntx)
 
 !! II.1. Vertical diffusion on tr
 !! ------------------------------
@@ -170,7 +183,7 @@ SUBROUTINE trczdf
 !!   ... Euler time stepping when starting from rest
         DO jk = 1, jpkm1
           z2dtt = zdt * rdt
-          zwy(jk,1) = trb(jk,jj,ji,jn)*e3t_back(jk,jj,ji)/e3t(jk,jj,ji) + z2dtt * tra(jk,jj,ji,jn)
+          zwy(jk,ntx) = trb(jk,jj,ji,jn)*e3t_back(jk,jj,ji)/e3t(jk,jj,ji) + z2dtt * tra(jk,jj,ji,jn)
         END DO
 
 !! Matrix inversion from the first level
@@ -208,22 +221,22 @@ SUBROUTINE trczdf
         ikstp1=ikst+1
         ikenm2=jpk-2
 
-        zwt(ikst,1)=zwd(ikst,1)
+        zwt(ikst,ntx)=zwd(ikst,ntx)
 
         DO jk=ikstp1,jpkm1
-            zwt(jk,1)=zwd(jk,1)-zwi(jk,1)*zws(jk-1,1)/zwt(jk-1,1)
+            zwt(jk,ntx)=zwd(jk,ntx)-zwi(jk,ntx)*zws(jk-1,ntx)/zwt(jk-1,ntx)
         END DO
 
-          zwz(ikst,1)=zwy(ikst,1)
+          zwz(ikst,ntx)=zwy(ikst,ntx)
 
         DO jk=ikstp1,jpkm1
-            zwz(jk,1)=zwy(jk,1)-zwi(jk, 1)/zwt(jk-1, 1)*zwz(jk-1, 1)
+            zwz(jk,ntx)=zwy(jk,ntx)-zwi(jk, ntx)/zwt(jk-1, ntx)*zwz(jk-1, ntx)
         END DO
 
-        zwx(jpkm1, 1)=zwz(jpkm1, 1)/zwt(jpkm1, 1)
+        zwx(jpkm1, ntx)=zwz(jpkm1, ntx)/zwt(jpkm1, ntx)
 
         DO jk=ikenm2,ikst,-1
-            zwx(jk, 1)=( zwz(jk, 1)-zws(jk, 1)*zwx(jk+1, 1) )/zwt(jk, 1)
+            zwx(jk, ntx)=( zwz(jk, ntx)-zws(jk, ntx)*zwx(jk+1, ntx) )/zwt(jk, ntx)
         END DO
 
 ! calculate flux due to vertical diffusion (on top face of the grid cell jk)
@@ -232,7 +245,7 @@ SUBROUTINE trczdf
       DO jk=1,jpkm1
 
          z2dtt = zdt * rdt
-         delta_tra(jk) = ( zwx(jk,1) - zwy(jk,1) ) / z2dtt * Aij * e3t(jk,jj,ji)! or trn(jk,jj,ji,jn+mytid)
+         delta_tra(jk) = ( zwx(jk,ntx) - zwy(jk,ntx) ) / z2dtt * Aij * e3t(jk,jj,ji)! or trn(jk,jj,ji,jn+mytid)
 
          IF (jk .EQ. 1) THEN
              int_tra(1)  = 0
@@ -258,16 +271,18 @@ SUBROUTINE trczdf
 !! (c a u t i o n:  tracer not its trend, Leap-frog scheme done
 !!                  it will not be done in trcnxt)
          DO jk = 1, jpkm1
-            tra(jk,jj,ji,jn) = zwx(jk,1) * tmask(jk,jj,ji)
+            tra(jk,jj,ji,jn) = zwx(jk,ntx) * tmask(jk,jj,ji)
          END DO
 
         END DO ! jv
+        !$acc end parallel loop
 
 !      end if
 
        END DO TRACER_LOOP
 !!!$omp    end parallel do
 
+       !$acc wait
 
 
        trczdfparttime = MPI_WTIME() - trczdfparttime   !cronometer-stop