Index: codes/icosagcm/trunk/src/dynamics/caldyn_gcm.F90
===================================================================
--- codes/icosagcm/trunk/src/dynamics/caldyn_gcm.F90	(revision 899)
+++ codes/icosagcm/trunk/src/dynamics/caldyn_gcm.F90	(revision 953)
@@ -15,4 +15,5 @@
   
   SUBROUTINE init_caldyn
+    USE abort_mod
     CHARACTER(len=255) :: def
     INTEGER            :: ind
@@ -21,5 +22,9 @@
     hydrostatic=.TRUE.
     CALL getin("hydrostatic",hydrostatic)
-  
+
+    IF (.NOT. hydrostatic) THEN
+       CALL abort_acc("hydrostatic /= .TRUE.")
+    END IF
+
     dysl=.NOT.hydrostatic ! dysl defaults to .FALSE. if hydrostatic, .TRUE. if NH                                              
     CALL getin("dysl",dysl)
@@ -34,4 +39,9 @@
     dysl_caldyn_coriolis=dysl
     CALL getin("dysl_caldyn_coriolis",dysl_caldyn_coriolis)
+
+    ! dysl is not supported with openACC
+    IF (dysl_geopot .OR. dysl_caldyn_fast .OR. dysl_slow_hydro .OR. dysl_pvort_only .OR. dysl_caldyn_coriolis) THEN
+       CALL abort_acc("dysl not supported")
+    END IF
 
     def='energy'
@@ -122,8 +132,8 @@
   SUBROUTINE allocate_caldyn
     CALL allocate_field(f_out_u,field_u,type_real,llm) 
-    CALL allocate_field(f_qu,field_u,type_real,llm) 
-    CALL allocate_field(f_qv,field_z,type_real,llm) 
-    CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk')
-    CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu')    
+    CALL allocate_field(f_qu,field_u,type_real,llm, ondevice=.TRUE.) 
+    CALL allocate_field(f_qv,field_z,type_real,llm, ondevice=.TRUE.) 
+    CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk', ondevice=.TRUE.)
+    CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu', ondevice=.TRUE.)    
     CALL allocate_field(f_planetvel,  field_u,type_real, name='planetvel') ! planetary velocity at r=a
     IF(.NOT.hydrostatic) THEN
Index: codes/icosagcm/trunk/src/dynamics/caldyn_hevi.F90
===================================================================
--- codes/icosagcm/trunk/src/dynamics/caldyn_hevi.F90	(revision 953)
+++ codes/icosagcm/trunk/src/dynamics/caldyn_hevi.F90	(revision 953)
@@ -0,0 +1,201 @@
+MODULE caldyn_hevi_mod
+  USE icosa
+  USE transfert_mod
+  USE caldyn_kernels_base_mod
+  USE caldyn_kernels_hevi_mod
+  USE caldyn_gcm_mod
+  IMPLICIT NONE
+  PRIVATE
+  PUBLIC caldyn_hevi
+
+CONTAINS
+  
+  SUBROUTINE caldyn_hevi(write_out,tau, f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
+       f_W, f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, &
+       f_du_slow, f_du_fast, f_dPhi_slow, f_dPhi_fast, f_dW_slow, f_dW_fast) 
+    USE icosa
+    USE observable_mod
+    USE disvert_mod, ONLY : caldyn_eta, eta_mass
+    USE vorticity_mod
+    USE kinetic_mod
+    USE theta2theta_rhodz_mod
+    USE wind_mod
+    USE mpipara
+    USE trace
+    USE omp_para
+    USE output_field_mod
+    USE checksum_mod
+    USE abort_mod
+    IMPLICIT NONE
+    LOGICAL,INTENT(IN)    :: write_out
+    REAL(rstd), INTENT(IN) :: tau
+    TYPE(t_field),POINTER :: f_phis(:)
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_mass(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_u(:)
+    TYPE(t_field),POINTER :: f_q(:)
+    TYPE(t_field),POINTER :: f_W(:)
+    TYPE(t_field),POINTER :: f_geopot(:)
+    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
+    TYPE(t_field) :: f_dps(:)
+    TYPE(t_field) :: f_dmass(:)
+    TYPE(t_field) :: f_dtheta_rhodz(:)
+    TYPE(t_field) :: f_du_slow(:)
+    TYPE(t_field) :: f_du_fast(:)
+    TYPE(t_field) :: f_dW_slow(:)
+    TYPE(t_field) :: f_dW_fast(:)
+    TYPE(t_field) :: f_dPhi_slow(:)
+    TYPE(t_field) :: f_dPhi_fast(:)
+    
+    REAL(rstd),POINTER :: ps(:), dps(:), phis(:)
+    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:,:)
+    REAL(rstd),POINTER :: du(:,:), dW(:,:), dPhi(:,:), hflux(:,:), wflux(:,:)
+    REAL(rstd),POINTER :: u(:,:), w(:,:), qu(:,:), qv(:,:)
+
+! temporary shared variable
+    REAL(rstd),POINTER  :: theta(:,:,:)  
+    REAL(rstd),POINTER  :: pk(:,:)
+    REAL(rstd),POINTER  :: geopot(:,:)
+    REAL(rstd),POINTER  :: convm(:,:) 
+    REAL(rstd),POINTER  :: wwuu(:,:)
+    REAL(rstd),POINTER  :: F_el(:,:), gradPhi2(:,:), w_il(:,:) , W_etadot(:,:), pres(:,:), m_il(:,:)
+    INTEGER :: ind
+    LOGICAL,SAVE :: first=.TRUE.
+!$OMP THREADPRIVATE(first)
+    
+    IF (first) THEN 
+      first=.FALSE.
+      IF(caldyn_eta==eta_mass) THEN
+         CALL init_message(f_ps,req_i1,req_ps)
+      ELSE
+         CALL init_message(f_mass,req_i1,req_mass)
+      END IF
+      CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
+      CALL init_message(f_u,req_e1_vect,req_u)
+      CALL init_message(f_qu,req_e1_scal,req_qu)
+      IF(.NOT.hydrostatic) THEN
+         CALL init_message(f_geopot,req_i1,req_geopot)
+         CALL init_message(f_w,req_i1,req_w)
+      END IF
+    ENDIF
+    
+    CALL trace_start("caldyn")
+    
+    IF(caldyn_eta==eta_mass) THEN
+       CALL send_message(f_ps,req_ps) ! COM00
+       CALL wait_message(req_ps) ! COM00
+    ELSE
+       CALL abort_acc("HEVI_scheme/!eta_mass")
+       CALL send_message(f_mass,req_mass) ! COM00
+       CALL wait_message(req_mass) ! COM00
+    END IF
+    CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01
+    CALL wait_message(req_theta_rhodz) ! COM01 Moved from caldyn_pvort
+
+    IF(.NOT.hydrostatic) THEN
+       CALL abort_acc("HEVI_scheme/!hydrostatic")
+       CALL send_message(f_geopot,req_geopot) ! COM03
+       CALL wait_message(req_geopot) ! COM03
+       CALL send_message(f_w,req_w) ! COM04
+       CALL wait_message(req_w) ! COM04
+    END IF
+    
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       ps=f_ps(ind)
+       theta_rhodz=f_theta_rhodz(ind)
+       mass=f_mass(ind)
+       theta = f_theta(ind)
+       CALL compute_theta(ps,theta_rhodz, mass,theta)
+       pk = f_pk(ind)
+       geopot = f_geopot(ind)
+       du=f_du_fast(ind)
+       IF(hydrostatic) THEN
+          !$acc kernels present(du) async
+          du(:,:)=0.0d0
+          !$acc end kernels
+          CALL compute_geopot(mass,theta, ps,pk,geopot)
+       ELSE
+          CALL abort_acc("HEVI_scheme/!hydrostatic")
+          phis = f_phis(ind)
+          W = f_W(ind)
+          dW = f_dW_fast(ind)
+          dPhi = f_dPhi_fast(ind)
+          ! reuse buffers
+          m_il = f_wil(ind)
+          pres = f_gradPhi2(ind)
+          CALL compute_caldyn_solver(tau,phis, mass,theta,pk,geopot,W, m_il,pres, dPhi,dW,du) ! computes d(Phi,W,du)_fast and updates Phi,W
+       END IF
+       u=f_u(ind)
+       CALL compute_caldyn_fast(tau,u,mass,theta,pk,geopot,du) ! computes du_fast and updates u
+    ENDDO
+    
+    CALL send_message(f_u,req_u) ! COM02
+    CALL wait_message(req_u)   ! COM02
+    
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       u=f_u(ind)
+       mass=f_mass(ind)
+       qu=f_qu(ind)
+       qv=f_qv(ind)
+       CALL compute_pvort_only(u,mass,qu,qv)
+    ENDDO
+    
+    CALL send_message(f_qu,req_qu) ! COM03
+    CALL wait_message(req_qu) ! COM03
+       
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       u=f_u(ind)
+       mass=f_mass(ind)
+       theta = f_theta(ind)
+       qu=f_qu(ind)
+       hflux=f_hflux(ind)
+       convm = f_dmass(ind)
+       dtheta_rhodz=f_dtheta_rhodz(ind)
+       du=f_du_slow(ind)
+
+       IF(hydrostatic) THEN
+          CALL compute_caldyn_slow_hydro(u,mass,hflux,du, .TRUE.)
+       ELSE
+          CALL abort_acc("HEVI_scheme/!hydrostatic")
+          W = f_W(ind)
+          dW = f_dW_slow(ind)
+          geopot = f_geopot(ind)
+          dPhi = f_dPhi_slow(ind)
+          F_el = f_Fel(ind)
+          gradPhi2 = f_gradPhi2(ind)
+          w_il = f_wil(ind)
+          CALL compute_caldyn_slow_NH(u,mass,geopot,W, F_el,gradPhi2,w_il, hflux,du,dPhi,dW)
+       END IF
+       CALL compute_caldyn_Coriolis(hflux,theta,qu,convm,dtheta_rhodz,du)
+       
+       IF(caldyn_eta==eta_mass) THEN
+          wflux=f_wflux(ind)
+          wwuu=f_wwuu(ind)
+          dps=f_dps(ind)
+          CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
+          IF(.NOT.hydrostatic) THEN
+             CALL abort_acc("HEVI_scheme/!hydrostatic")
+             W_etadot=f_Wetadot(ind)
+             CALL compute_caldyn_vert_NH(mass,geopot,W,wflux, W_etadot, du,dPhi,dW)
+          END IF
+       END IF
+    ENDDO
+    
+!$OMP BARRIER
+    !    CALL check_mass_conservation(f_ps,f_dps)
+    CALL trace_end("caldyn_hevi")
+!!$OMP BARRIER
+    
+  END SUBROUTINE caldyn_hevi
+
+END MODULE caldyn_hevi_mod
Index: codes/icosagcm/trunk/src/dynamics/caldyn_hevi.f90
===================================================================
--- codes/icosagcm/trunk/src/dynamics/caldyn_hevi.f90	(revision 899)
+++ 	(revision )
@@ -1,192 +1,0 @@
-MODULE caldyn_hevi_mod
-  USE icosa
-  USE transfert_mod
-  USE caldyn_kernels_base_mod
-  USE caldyn_kernels_hevi_mod
-  USE caldyn_gcm_mod
-  IMPLICIT NONE
-  PRIVATE
-  PUBLIC caldyn_hevi
-
-CONTAINS
-  
-  SUBROUTINE caldyn_hevi(write_out,tau, f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
-       f_W, f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, &
-       f_du_slow, f_du_fast, f_dPhi_slow, f_dPhi_fast, f_dW_slow, f_dW_fast) 
-    USE icosa
-    USE observable_mod
-    USE disvert_mod, ONLY : caldyn_eta, eta_mass
-    USE vorticity_mod
-    USE kinetic_mod
-    USE theta2theta_rhodz_mod
-    USE wind_mod
-    USE mpipara
-    USE trace
-    USE omp_para
-    USE output_field_mod
-    USE checksum_mod
-    IMPLICIT NONE
-    LOGICAL,INTENT(IN)    :: write_out
-    REAL(rstd), INTENT(IN) :: tau
-    TYPE(t_field),POINTER :: f_phis(:)
-    TYPE(t_field),POINTER :: f_ps(:)
-    TYPE(t_field),POINTER :: f_mass(:)
-    TYPE(t_field),POINTER :: f_theta_rhodz(:)
-    TYPE(t_field),POINTER :: f_u(:)
-    TYPE(t_field),POINTER :: f_q(:)
-    TYPE(t_field),POINTER :: f_W(:)
-    TYPE(t_field),POINTER :: f_geopot(:)
-    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
-    TYPE(t_field) :: f_dps(:)
-    TYPE(t_field) :: f_dmass(:)
-    TYPE(t_field) :: f_dtheta_rhodz(:)
-    TYPE(t_field) :: f_du_slow(:)
-    TYPE(t_field) :: f_du_fast(:)
-    TYPE(t_field) :: f_dW_slow(:)
-    TYPE(t_field) :: f_dW_fast(:)
-    TYPE(t_field) :: f_dPhi_slow(:)
-    TYPE(t_field) :: f_dPhi_fast(:)
-    
-    REAL(rstd),POINTER :: ps(:), dps(:), phis(:)
-    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:,:)
-    REAL(rstd),POINTER :: du(:,:), dW(:,:), dPhi(:,:), hflux(:,:), wflux(:,:)
-    REAL(rstd),POINTER :: u(:,:), w(:,:), qu(:,:), qv(:,:)
-
-! temporary shared variable
-    REAL(rstd),POINTER  :: theta(:,:,:)  
-    REAL(rstd),POINTER  :: pk(:,:)
-    REAL(rstd),POINTER  :: geopot(:,:)
-    REAL(rstd),POINTER  :: convm(:,:) 
-    REAL(rstd),POINTER  :: wwuu(:,:)
-    REAL(rstd),POINTER  :: F_el(:,:), gradPhi2(:,:), w_il(:,:) , W_etadot(:,:), pres(:,:), m_il(:,:)
-    INTEGER :: ind
-    LOGICAL,SAVE :: first=.TRUE.
-!$OMP THREADPRIVATE(first)
-    
-    IF (first) THEN 
-      first=.FALSE.
-      IF(caldyn_eta==eta_mass) THEN
-         CALL init_message(f_ps,req_i1,req_ps)
-      ELSE
-         CALL init_message(f_mass,req_i1,req_mass)
-      END IF
-      CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
-      CALL init_message(f_u,req_e1_vect,req_u)
-      CALL init_message(f_qu,req_e1_scal,req_qu)
-      IF(.NOT.hydrostatic) THEN
-         CALL init_message(f_geopot,req_i1,req_geopot)
-         CALL init_message(f_w,req_i1,req_w)
-      END IF
-    ENDIF
-    
-    CALL trace_start("caldyn")
-    
-    IF(caldyn_eta==eta_mass) THEN
-       CALL send_message(f_ps,req_ps) ! COM00
-       CALL wait_message(req_ps) ! COM00
-    ELSE
-       CALL send_message(f_mass,req_mass) ! COM00
-       CALL wait_message(req_mass) ! COM00
-    END IF
-    CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01
-    CALL wait_message(req_theta_rhodz) ! COM01 Moved from caldyn_pvort
-
-    IF(.NOT.hydrostatic) THEN
-       CALL send_message(f_geopot,req_geopot) ! COM03
-       CALL wait_message(req_geopot) ! COM03
-       CALL send_message(f_w,req_w) ! COM04
-       CALL wait_message(req_w) ! COM04
-    END IF
-    
-    DO ind=1,ndomain
-       IF (.NOT. assigned_domain(ind)) CYCLE
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-       ps=f_ps(ind)
-       theta_rhodz=f_theta_rhodz(ind)
-       mass=f_mass(ind)
-       theta = f_theta(ind)
-       CALL compute_theta(ps,theta_rhodz, mass,theta)
-       pk = f_pk(ind)
-       geopot = f_geopot(ind)
-       du=f_du_fast(ind)
-       IF(hydrostatic) THEN
-          du(:,:)=0.
-          CALL compute_geopot(mass,theta, ps,pk,geopot)
-       ELSE
-          phis = f_phis(ind)
-          W = f_W(ind)
-          dW = f_dW_fast(ind)
-          dPhi = f_dPhi_fast(ind)
-          ! reuse buffers
-          m_il = f_wil(ind)
-          pres = f_gradPhi2(ind)
-          CALL compute_caldyn_solver(tau,phis, mass,theta,pk,geopot,W, m_il,pres, dPhi,dW,du) ! computes d(Phi,W,du)_fast and updates Phi,W
-       END IF
-       u=f_u(ind)
-       CALL compute_caldyn_fast(tau,u,mass,theta,pk,geopot,du) ! computes du_fast and updates u
-    ENDDO
-    
-    CALL send_message(f_u,req_u) ! COM02
-    CALL wait_message(req_u)   ! COM02
-    
-    DO ind=1,ndomain
-       IF (.NOT. assigned_domain(ind)) CYCLE
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-       u=f_u(ind)
-       mass=f_mass(ind)
-       qu=f_qu(ind)
-       qv=f_qv(ind)
-       CALL compute_pvort_only(u,mass,qu,qv)
-    ENDDO
-    
-    CALL send_message(f_qu,req_qu) ! COM03
-    CALL wait_message(req_qu) ! COM03
-       
-    DO ind=1,ndomain
-       IF (.NOT. assigned_domain(ind)) CYCLE
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-       u=f_u(ind)
-       mass=f_mass(ind)
-       theta = f_theta(ind)
-       qu=f_qu(ind)
-       hflux=f_hflux(ind)
-       convm = f_dmass(ind)
-       dtheta_rhodz=f_dtheta_rhodz(ind)
-       du=f_du_slow(ind)
-
-       IF(hydrostatic) THEN
-          CALL compute_caldyn_slow_hydro(u,mass,hflux,du, .TRUE.)
-       ELSE
-          W = f_W(ind)
-          dW = f_dW_slow(ind)
-          geopot = f_geopot(ind)
-          dPhi = f_dPhi_slow(ind)
-          F_el = f_Fel(ind)
-          gradPhi2 = f_gradPhi2(ind)
-          w_il = f_wil(ind)
-          CALL compute_caldyn_slow_NH(u,mass,geopot,W, F_el,gradPhi2,w_il, hflux,du,dPhi,dW)
-       END IF
-       CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
-       IF(caldyn_eta==eta_mass) THEN
-          wflux=f_wflux(ind)
-          wwuu=f_wwuu(ind)
-          dps=f_dps(ind)
-          CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
-          IF(.NOT.hydrostatic) THEN
-             W_etadot=f_Wetadot(ind)
-             CALL compute_caldyn_vert_NH(mass,geopot,W,wflux, W_etadot, du,dPhi,dW)
-          END IF
-       END IF
-    ENDDO
-    
-!$OMP BARRIER
-    !    CALL check_mass_conservation(f_ps,f_dps)
-    CALL trace_end("caldyn_hevi")
-!!$OMP BARRIER
-    
-  END SUBROUTINE caldyn_hevi
-
-END MODULE caldyn_hevi_mod
Index: codes/icosagcm/trunk/src/dynamics/caldyn_kernels_base.F90
===================================================================
--- codes/icosagcm/trunk/src/dynamics/caldyn_kernels_base.F90	(revision 899)
+++ codes/icosagcm/trunk/src/dynamics/caldyn_kernels_base.F90	(revision 953)
@@ -5,4 +5,5 @@
   USE omp_para
   USE trace
+  USE abort_mod
   IMPLICIT NONE
   PRIVATE
@@ -25,5 +26,4 @@
 
   !**************************** Geopotential *****************************
-
   SUBROUTINE compute_geopot(rhodz,theta, ps,pk,geopot) 
     REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
@@ -46,4 +46,5 @@
 
     IF(dysl_geopot) THEN
+      CALL abort_acc("HEVI_scheme/!dysl_geopot")
 #include "../kernels/compute_geopot.k90"
     ELSE
@@ -52,26 +53,35 @@
     ! Works also when caldyn_eta=eta_mass          
 
+    !$acc data present(rhodz(:,:), ps(:), pk(:,:), geopot(:,:), theta(:,:,:)) async
+
     IF(boussinesq) THEN ! compute geopotential and pk=Lagrange multiplier
+       CALL abort_acc("HEVI_scheme/boussinesq")
        ! specific volume 1 = dphi/g/rhodz
        !         IF (is_omp_level_master) THEN ! no openMP on vertical due to dependency
        DO l = 1,llm
+          !$acc parallel loop
           !DIR$ SIMD
           DO ij=ij_omp_begin_ext,ij_omp_end_ext         
              geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
           ENDDO
+          !$acc end parallel loop
        ENDDO
        ! use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure) 
        ! uppermost layer
+       !$acc parallel loop
        !DIR$ SIMD
        DO ij=ij_begin_ext,ij_end_ext         
           pk(ij,llm) = ptop + (.5*g)*theta(ij,llm,1)*rhodz(ij,llm)
        END DO
+       !$acc end parallel loop
        ! other layers
        DO l = llm-1, 1, -1
           !          !$OMP DO SCHEDULE(STATIC) 
+          !$acc parallel loop
           !DIR$ SIMD
           DO ij=ij_begin_ext,ij_end_ext         
              pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l,1)*rhodz(ij,l)+theta(ij,l+1,1)*rhodz(ij,l+1))
           END DO
+          !$acc end parallel loop
        END DO
        ! now pk contains the Lagrange multiplier (pressure)
@@ -81,9 +91,15 @@
        SELECT CASE(caldyn_thermo)
           CASE(thermo_theta, thermo_entropy)
+             !$acc parallel loop async
              !DIR$ SIMD
              DO ij=ij_omp_begin_ext,ij_omp_end_ext
                 pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
              END DO
+
              ! other layers
+             ! We use kernels here instead of "loop seq" + "loop gang vector"
+             ! to be sure the code really abides the standard. In practice,
+             ! it seems like the compiler interchanges the loops.
+             !$acc kernels async
              DO l = llm-1, 1, -1
                 !DIR$ SIMD
@@ -92,6 +108,8 @@
                 END DO
              END DO
+             !$acc end kernels
              ! surface pressure (for diagnostics)
              IF(caldyn_eta==eta_lag) THEN
+                !$acc parallel loop async
                 DO ij=ij_omp_begin_ext,ij_omp_end_ext         
                    ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)
@@ -99,10 +117,15 @@
              END IF
           CASE(thermo_moist) ! theta(ij,l,2) = qv = mv/md
+             CALL abort_acc("HEVI_scheme/thermo_moist")
+             !$acc parallel loop
              !DIR$ SIMD
              DO ij=ij_omp_begin_ext,ij_omp_end_ext
                 pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)*(1.+theta(ij,l,2))
              END DO
+             !$acc end parallel loop
+
              ! other layers
              DO l = llm-1, 1, -1
+                !$acc parallel loop
                 !DIR$ SIMD
                 DO ij=ij_omp_begin_ext,ij_omp_end_ext         
@@ -111,51 +134,70 @@
                         rhodz(ij,l+1)*(1.+theta(ij,l+1,2)) )
                 END DO
+                !$acc end parallel loop
              END DO
              ! surface pressure (for diagnostics)
              IF(caldyn_eta==eta_lag) THEN
+                !$acc parallel loop
                 DO ij=ij_omp_begin_ext,ij_omp_end_ext         
                    ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)*(1.+theta(ij,l,2))
                 END DO
+                !$acc end parallel loop
              END IF
           END SELECT
 
-       DO l = 1,llm
-          SELECT CASE(caldyn_thermo)
+        SELECT CASE(caldyn_thermo)
           CASE(thermo_theta)
-             !DIR$ SIMD
-             DO ij=ij_omp_begin_ext,ij_omp_end_ext
-                p_ik = pk(ij,l)
-                exner_ik = cpp * (p_ik/preff) ** kappa
-                pk(ij,l) = exner_ik
-                ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
-                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l,1)*exner_ik/p_ik
-             ENDDO
+            ! We use kernels here instead of "loop seq" + "loop gang vector"
+            ! to be sure the code really abides the standard. In practice,
+            ! it seems like the compiler interchanges the loops.
+            !$acc kernels async
+            DO l = 1,llm
+               !DIR$ SIMD
+               DO ij=ij_omp_begin_ext,ij_omp_end_ext
+                  p_ik = pk(ij,l)
+                  exner_ik = cpp * (p_ik/preff) ** kappa
+                  pk(ij,l) = exner_ik
+                  ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
+                  geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l,1)*exner_ik/p_ik
+               ENDDO
+            ENDDO
+            !$acc end kernels
+
           CASE(thermo_entropy) ! theta is in fact entropy = cpp*log(theta/Treff) = cpp*log(T/Treff) - Rd*log(p/preff)
-             !DIR$ SIMD
-             DO ij=ij_omp_begin_ext,ij_omp_end_ext
-                p_ik = pk(ij,l)
-                temp_ik = Treff*exp((theta(ij,l,1) + Rd*log(p_ik/preff))/cpp)
-                pk(ij,l) = temp_ik
-                ! specific volume v = Rd*T/p = dphi/g/rhodz
-                geopot(ij,l+1) = geopot(ij,l) + (g*Rd)*rhodz(ij,l)*temp_ik/p_ik
+             CALL abort_acc("HEVI_scheme/thermo_entropy")
+             DO l = 1,llm
+               !$acc parallel loop
+               !DIR$ SIMD
+               DO ij=ij_omp_begin_ext,ij_omp_end_ext
+                  p_ik = pk(ij,l)
+                  temp_ik = Treff*exp((theta(ij,l,1) + Rd*log(p_ik/preff))/cpp)
+                  pk(ij,l) = temp_ik
+                  ! specific volume v = Rd*T/p = dphi/g/rhodz
+                  geopot(ij,l+1) = geopot(ij,l) + (g*Rd)*rhodz(ij,l)*temp_ik/p_ik
+               ENDDO
+               !$acc end parallel loop
              ENDDO
           CASE(thermo_moist) ! theta is moist pseudo-entropy per dry air mass
-             DO ij=ij_omp_begin_ext,ij_omp_end_ext
-                p_ik = pk(ij,l)
-                qv = theta(ij,l,2) ! water vaper mixing ratio = mv/md
-                Rmix = Rd+qv*Rv
-                chi = ( theta(ij,l,1) + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff)
-                temp_ik = Treff*exp(chi)
-                pk(ij,l) = temp_ik
-                ! specific volume v = R*T/p = dphi/g/rhodz
-                ! R = (Rd + qv.Rv)/(1+qv)
-                geopot(ij,l+1) = geopot(ij,l) + g*Rmix*rhodz(ij,l)*temp_ik/(p_ik*(1+qv))
+             CALL abort_acc("HEVI_scheme/thermo_moist")
+             DO l = 1,llm
+               !$acc parallel loop
+               DO ij=ij_omp_begin_ext,ij_omp_end_ext
+                  p_ik = pk(ij,l)
+                  qv = theta(ij,l,2) ! water vaper mixing ratio = mv/md
+                  Rmix = Rd+qv*Rv
+                  chi = ( theta(ij,l,1) + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff)
+                  temp_ik = Treff*exp(chi)
+                  pk(ij,l) = temp_ik
+                  ! specific volume v = R*T/p = dphi/g/rhodz
+                  ! R = (Rd + qv.Rv)/(1+qv)
+                  geopot(ij,l+1) = geopot(ij,l) + g*Rmix*rhodz(ij,l)*temp_ik/(p_ik*(1+qv))
+               ENDDO
+               !$acc end parallel loop
              ENDDO
           CASE DEFAULT
              STOP
-          END SELECT
-       ENDDO
+        END SELECT
     END IF
-
+    !$acc end data
     END IF ! dysl
 
@@ -167,5 +209,6 @@
   END SUBROUTINE compute_geopot
 
-  SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
+  SUBROUTINE compute_caldyn_vert(u, theta, rhodz, convm, wflux, wwuu, dps, dtheta_rhodz, du)
+    USE disvert_mod, ONLY : bp
     REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
     REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn)
@@ -193,7 +236,15 @@
     CALL trace_start("compute_caldyn_vert")
 
+    !$acc data async &
+    !$acc present(rhodz(:,:), u(:,:), wwuu(:,:), wflux(:,:), dps(:), convm(:,:), du(:,:), dtheta_rhodz(:,:,:), theta(:,:,:), bp(:))
+
+
     !$OMP BARRIER   
 !!! cumulate mass flux convergence from top to bottom
     !  IF (is_omp_level_master) THEN
+    ! We use kernels here instead of "loop seq" + "loop gang vector"
+    ! to be sure the code really abides the standard. In practice,
+    ! it seems like the compiler interchanges the loops.
+    !$acc kernels async
     DO  l = llm-1, 1, -1
        !    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
@@ -205,4 +256,5 @@
        ENDDO
     ENDDO
+    !$acc end kernels
     !  ENDIF
 
@@ -213,4 +265,5 @@
     ! compute dps
     IF (is_omp_first_level) THEN
+       !$acc parallel loop  async
        !DIR$ SIMD
        DO ij=ij_begin,ij_end         
@@ -221,4 +274,5 @@
 
 !!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC)
+    !$acc parallel loop collapse(2) async
     DO l=ll_beginp1,ll_end
        !    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
@@ -233,5 +287,5 @@
     !--> flush wflux  
     !$OMP BARRIER
-
+    !$acc parallel loop collapse(3) async
     DO iq=1,nqdyn
        DO l=ll_begin,ll_endm1
@@ -242,4 +296,8 @@
           END DO
        END DO
+    END DO
+    
+    !$acc parallel loop collapse(3) async
+    DO iq=1,nqdyn
        DO l=ll_beginp1,ll_end
           !DIR$ SIMD
@@ -252,4 +310,5 @@
 
     ! Compute vertical transport
+    !$acc parallel loop collapse(2) async
     DO l=ll_beginp1,ll_end
        !DIR$ SIMD
@@ -265,4 +324,5 @@
 
     ! Add vertical transport to du
+    !$acc parallel loop collapse(2) async
     DO l=ll_begin,ll_end
        !DIR$ SIMD
@@ -283,4 +343,5 @@
     !   ENDDO
 
+    !$acc end data
     CALL trace_end("compute_caldyn_vert")
 
Index: codes/icosagcm/trunk/src/dynamics/caldyn_kernels_hevi.F90
===================================================================
--- codes/icosagcm/trunk/src/dynamics/caldyn_kernels_hevi.F90	(revision 899)
+++ codes/icosagcm/trunk/src/dynamics/caldyn_kernels_hevi.F90	(revision 953)
@@ -6,4 +6,5 @@
   USE transfert_mod
   USE caldyn_kernels_base_mod
+  USE abort_mod
   IMPLICIT NONE
   PRIVATE
@@ -20,4 +21,5 @@
 
   SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta)
+    USE disvert_mod, ONLY : mass_dbk, mass_dak
     REAL(rstd),INTENT(IN)    :: ps(iim*jjm)
     REAL(rstd),INTENT(IN)    :: theta_rhodz(iim*jjm,llm,nqdyn)
@@ -26,7 +28,12 @@
     INTEGER :: ij,l,iq
     REAL(rstd) :: m
+
     CALL trace_start("compute_theta")
 
+    !$acc data async &
+    !$acc present(ps(:), theta_rhodz(:,:,:), rhodz(:,:), theta(:,:,:), mass_dak(:), mass_dbk(:))
+
     IF(caldyn_eta==eta_mass) THEN ! Compute mass
+       !$acc parallel loop collapse(2) async
        DO l = ll_begin,ll_end
           !DIR$ SIMD
@@ -38,6 +45,7 @@
     END IF
 
-    DO l = ll_begin,ll_end
-       DO iq=1,nqdyn
+    !$acc parallel loop collapse(3) async
+    DO iq=1,nqdyn
+       DO l = ll_begin,ll_end
           !DIR$ SIMD
           DO ij=ij_begin_ext,ij_end_ext
@@ -47,8 +55,10 @@
     END DO
 
+    !$acc end data
     CALL trace_end("compute_theta")
   END SUBROUTINE compute_theta
 
   SUBROUTINE compute_pvort_only(u,rhodz,qu,qv)
+    USE geometry, ONLY : Av, Riv2, fv
     REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
     REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
@@ -61,9 +71,15 @@
     CALL trace_start("compute_pvort_only")  
 !!! Compute shallow-water potential vorticity
+
     IF(dysl_pvort_only) THEN
+      CALL abort_acc("HEVI_scheme/dysl_pvort_only")
 #include "../kernels/pvort_only.k90"
     ELSE
 
+    !$acc data async &
+    !$acc present(rhodz(:,:), u(:,:), qu(:,:), qv(:,:), Av(:), Riv2(:,:), fv(:))
+    
     radius_m2=radius**(-2)
+    !$acc parallel loop collapse(2) async
     DO l = ll_begin,ll_end
        !DIR$ SIMD
@@ -85,5 +101,8 @@
           qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
        ENDDO
-
+    END DO
+
+    !$acc parallel loop collapse(2) async
+    DO l = ll_begin,ll_end
        !DIR$ SIMD
        DO ij=ij_begin,ij_end
@@ -94,6 +113,9 @@
 
     ENDDO
+
+    !$acc end data
    
     END IF ! dysl
+
     CALL trace_end("compute_pvort_only")
 
@@ -409,111 +431,126 @@
 
     IF(dysl_caldyn_fast) THEN
+       CALL abort_acc("HEVI_scheme/dysl_caldyn_fast")
 #include "../kernels/caldyn_fast.k90"
     ELSE
 
-    ! Compute Bernoulli term
-    IF(boussinesq) THEN
-       DO l=ll_begin,ll_end
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-             berni(ij,l) = pk(ij,l)
-             ! from now on pk contains the vertically-averaged geopotential
-             pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
+
+       ! Default case : ported to openacc
+       IF(.not. boussinesq .and. caldyn_thermo /= thermo_moist) THEN
+          !$acc data async &
+          !$acc present(rhodz(:,:), u(:,:),pk(:,:), geopot(:,:), theta(:,:,:), du(:,:))
+
+#define BERNI(ij) .5*(geopot(ij,l)+geopot(ij,l+1))
+          !$acc parallel loop collapse(2) async
+          DO l=ll_begin,ll_end
+             !DIR$ SIMD
+             DO ij=ij_begin,ij_end
+                due_right =  0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))  &
+                     *(pk(ij+t_right,l)-pk(ij,l))        &
+                     +  BERNI(ij+t_right)-BERNI(ij)
+                due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))    &
+                     *(pk(ij+t_lup,l)-pk(ij,l))          &
+                     +  BERNI(ij+t_lup)-BERNI(ij)
+                due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) &
+                     *(pk(ij+t_ldown,l)-pk(ij,l))      &
+                     +   BERNI(ij+t_ldown)-BERNI(ij)
+                du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
+                du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
+                du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
+                u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
+                u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
+                u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
+             END DO
           END DO
-       END DO
-    ELSE ! compressible
-
-       DO l=ll_begin,ll_end
-          SELECT CASE(caldyn_thermo)
-          CASE(thermo_theta) ! vdp = theta.dpi => B = Phi
-             !DIR$ SIMD
-             DO ij=ij_begin,ij_end         
-                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
-             END DO
-          CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s) 
-             !DIR$ SIMD
-             DO ij=ij_begin,ij_end         
-                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
-                     + pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy
-             END DO
-          CASE(thermo_moist) 
-             !DIR$ SIMD
-             DO ij=ij_begin,ij_end         
-                ! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T)
-                ! Bd = Phi + gibbs_d
-                ! Bv = Phi + gibbs_v
-                ! pk=temperature, theta=entropy
-                qv = theta(ij,l,2)
-                temp = pk(ij,l)
-                chi = log(temp/Treff)
-                nu = (chi*(cpp+qv*cppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff)
-                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
-                     + temp*(cpp*(1.-chi)+Rd*nu)
-                berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
-                     + temp*(cppv*(1.-chi)+Rv*nu)
-            END DO
-          END SELECT
-       END DO
-
-    END IF ! Boussinesq/compressible
+
+#undef BERNI
+
+          !$acc end data
+       ELSE ! Not default case : not ported to openacc
+          CALL abort_acc("HEVI_scheme/compute_caldyn_fast")
+          !!$acc parallel loop gang private(berni(:,:), berniv(:,:)) async
+          DO l=ll_begin,ll_end
+             ! Compute Bernoulliterm
+             IF(boussinesq) THEN
+                !!$acc loop vector
+                !DIR$ SIMD
+                DO ij=ij_begin,ij_end
+                   berni(ij,l) = pk(ij,l)
+                   ! from now on pk contains the vertically-averaged geopotential
+                   pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
+                END DO
+             ELSE ! compressible
+                SELECT CASE(caldyn_thermo)
+                CASE(thermo_theta) ! vdp = theta.dpi => B = Phi
+                   CALL dynamico_abort("Case already dealt with")
+                CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s)
+                   !$acc loop vector
+                   !DIR$ SIMD
+                   DO ij=ij_begin,ij_end
+                      berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+                           + pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy
+                   END DO
+                CASE(thermo_moist)
+!!$acc loop vector
+                   !DIR$ SIMD
+                   DO ij=ij_begin,ij_end
+                      ! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T)
+                      ! Bd = Phi + gibbs_d
+                      ! Bv = Phi + gibbs_v
+                      ! pk=temperature, theta=entropy
+                      qv = theta(ij,l,2)
+                      temp = pk(ij,l)
+                      chi = log(temp/Treff)
+                      nu = (chi*(cpp+qv*cppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff)
+                      berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+                           + temp*(cpp*(1.-chi)+Rd*nu)
+                      berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+                           + temp*(cppv*(1.-chi)+Rv*nu)
+                   END DO
+                END SELECT
+             END IF ! Boussinesq/compressible
 
 !!! u:=u+tau*du, du = -grad(B)-theta.grad(pi)
-    DO l=ll_begin,ll_end
-       IF(caldyn_thermo == thermo_moist) THEN
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-             due_right =  berni(ij+t_right,l)-berni(ij,l)       &
-                  + 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))   &
-                       *(pk(ij+t_right,l)-pk(ij,l))             &
-                  + 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2))   &
-                       *(berniv(ij+t_right,l)-berniv(ij,l))
-             
-             due_lup = berni(ij+t_lup,l)-berni(ij,l)            &
-                  + 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))     &
-                       *(pk(ij+t_lup,l)-pk(ij,l))               &
-                  + 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2))     &
-                       *(berniv(ij+t_lup,l)-berniv(ij,l))
-             
-             due_ldown = berni(ij+t_ldown,l)-berni(ij,l)        &
-                  + 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))   &
-                       *(pk(ij+t_ldown,l)-pk(ij,l))             &
-                  + 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2))   &
-                       *(berniv(ij+t_ldown,l)-berniv(ij,l))
-             
-             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
-             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
-             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
-             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
-             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
-             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
-          END DO
-       ELSE
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-             due_right =  0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))  &
-                  *(pk(ij+t_right,l)-pk(ij,l))        &
-                  +  berni(ij+t_right,l)-berni(ij,l)
-             due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))    &
-                  *(pk(ij+t_lup,l)-pk(ij,l))          &
-                  +  berni(ij+t_lup,l)-berni(ij,l)
-             due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) &
-                  *(pk(ij+t_ldown,l)-pk(ij,l))      &
-                  +   berni(ij+t_ldown,l)-berni(ij,l)
-             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
-             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
-             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
-             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
-             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
-             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
+             IF(caldyn_thermo == thermo_moist) THEN
+!!$acc loop vector
+                !DIR$ SIMD
+                DO ij=ij_begin,ij_end
+                   due_right =  berni(ij+t_right,l)-berni(ij,l)       &
+                        + 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))   &
+                        *(pk(ij+t_right,l)-pk(ij,l))             &
+                        + 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2))   &
+                        *(berniv(ij+t_right,l)-berniv(ij,l))
+
+                   due_lup = berni(ij+t_lup,l)-berni(ij,l)            &
+                        + 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))     &
+                        *(pk(ij+t_lup,l)-pk(ij,l))               &
+                        + 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2))     &
+                        *(berniv(ij+t_lup,l)-berniv(ij,l))
+
+                   due_ldown = berni(ij+t_ldown,l)-berni(ij,l)        &
+                        + 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))   &
+                        *(pk(ij+t_ldown,l)-pk(ij,l))             &
+                        + 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2))   &
+                        *(berniv(ij+t_ldown,l)-berniv(ij,l))
+
+                   du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
+                   du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
+                   du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
+                   u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
+                   u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
+                   u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
+                END DO
+             ELSE
+                CALL dynamico_abort("Case already dealt with")
+             END IF
           END DO
        END IF
-    END DO
-
-    END IF ! dysl
+    END IF
     CALL trace_end("compute_caldyn_fast")
 
   END SUBROUTINE compute_caldyn_fast
-
+    
   SUBROUTINE compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
+    USE geometry, ONLY : Ai, wee
     REAL(rstd),INTENT(IN)  :: hflux(3*iim*jjm,llm)  ! hflux in kg/s
     REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn) ! active scalars
@@ -534,10 +571,15 @@
 
     ELSE
-#define FTHETA(ij) Ftheta(ij,1)
-
-    DO l=ll_begin, ll_end
-       ! compute theta flux
-       DO iq=1,nqdyn
-       !DIR$ SIMD
+#define FTHETA(ij) Ftheta(ij,l)
+
+    !$acc data async &
+    !$acc create(Ftheta(3*iim*jjm,llm)) &
+    !$acc present(qu(:,:), hflux(:,:), convm(:,:), du(:,:), dtheta_rhodz(:,:,:), theta(:,:,:), Ai(:), wee(:,:,:))
+
+    ! compute theta flux
+    DO iq=1,nqdyn
+       !$acc parallel loop collapse(2) present(Ftheta, theta) async
+       DO l=ll_begin, ll_end
+          !DIR$ SIMD
           DO ij=ij_begin_ext,ij_end_ext      
              FTHETA(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) &
@@ -548,6 +590,10 @@
                   * hflux(ij+u_ldown,l)
           END DO
+       END DO
+
+       !$acc parallel loop collapse(2) present(Ftheta) async
+       DO l=ll_begin, ll_end
           ! horizontal divergence of fluxes
-       !DIR$ SIMD
+          !DIR$ SIMD
           DO ij=ij_begin,ij_end         
              ! dtheta_rhodz =  -div(flux.theta)
@@ -561,5 +607,8 @@
           END DO
        END DO
-
+    END DO
+
+    !$acc parallel loop collapse(2) present(Ai) async
+    DO l=ll_begin, ll_end
        !DIR$ SIMD
        DO ij=ij_begin,ij_end         
@@ -578,5 +627,5 @@
 
     CASE(energy) ! energy-conserving TRiSK
-
+       !$acc parallel loop collapse(2) present(qu, wee, du) async
        DO l=ll_begin,ll_end
           !DIR$ SIMD
@@ -622,5 +671,6 @@
 
     CASE(enstrophy) ! enstrophy-conserving TRiSK
-
+       
+       !$acc parallel loop collapse(2) present(wee, du) async
        DO l=ll_begin,ll_end
           !DIR$ SIMD
@@ -665,7 +715,11 @@
           END DO
        END DO
+       !$acc end parallel loop
     CASE DEFAULT
        STOP
     END SELECT
+
+    !$acc end data
+
 #undef FTHETA
 
@@ -676,5 +730,6 @@
   END SUBROUTINE compute_caldyn_Coriolis
 
-  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du, zero)
+  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du,zero)
+    USE geometry, ONLY : Ai, le_de
     LOGICAL, INTENT(IN) :: zero
     REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity"
@@ -684,5 +739,4 @@
     
     REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
-    REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function
     REAL(rstd) :: uu_right, uu_lup, uu_ldown, ke, uu
     INTEGER :: ij,l
@@ -698,9 +752,17 @@
      ELSE
 
-#define BERNI(ij) berni1(ij)
+#define BERNI(ij) berni(ij,l)
 
     DO l = ll_begin, ll_end
-       !  Compute mass fluxes
        IF (caldyn_conserv==energy) CALL test_message(req_qu) 
+    END DO
+    
+    !$acc data async &
+    !$acc create(berni(:,ll_begin:ll_end)) &
+    !$acc present(rhodz(:,ll_begin:ll_end), u(:,ll_begin:ll_end), hflux(:,:), du(:,ll_begin:ll_end), Ai(:), le_de(:))
+
+    !  Compute mass fluxes
+    !$acc parallel loop collapse(2) async
+    DO l = ll_begin, ll_end
        !DIR$ SIMD
        DO ij=ij_begin_ext,ij_end_ext
@@ -715,5 +777,9 @@
           hflux(ij+u_ldown,l)=uu_ldown
        ENDDO
-       ! Compute Bernoulli=kinetic energy
+    END DO
+
+    ! Compute Bernoulli=kinetic energy
+    !$acc parallel loop collapse(2) async
+    DO l = ll_begin, ll_end
        !DIR$ SIMD
        DO ij=ij_begin,ij_end         
@@ -726,6 +792,10 @@
                le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
        ENDDO
-       ! Compute du=-grad(Bernoulli)
-       IF(zero) THEN
+    END DO
+
+    ! Compute du=-grad(Bernoulli)
+    IF(zero) THEN
+       !$acc parallel loop collapse(2) async
+       DO l = ll_begin, ll_end
           !DIR$ SIMD
           DO ij=ij_begin,ij_end
@@ -734,5 +804,8 @@
              du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
           END DO
-       ELSE
+       END DO
+    ELSE
+       !$acc parallel loop collapse(2) async
+       DO l = ll_begin, ll_end
           !DIR$ SIMD
           DO ij=ij_begin,ij_end
@@ -744,10 +817,13 @@
                   ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
           END DO
-       END IF
-    END DO
+       END DO
+    END IF
+
+    !$acc end data
 
 #undef BERNI
 
     END IF ! dysl
+
     CALL trace_end("compute_caldyn_slow_hydro")    
   END SUBROUTINE compute_caldyn_slow_hydro
