Index: /codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90
===================================================================
--- /codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90	(revision 916)
+++ /codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90	(revision 917)
@@ -5,5 +5,5 @@
   USE caldyn_kernels_hevi_mod
   USE caldyn_kernels_base_mod
-  USE compute_theta_mod, ONLY : compute_theta
+  USE compute_theta_mod, ONLY : compute_theta => compute_theta_manual
   USE compute_caldyn_kv_mod, ONLY : compute_caldyn_kv
   USE compute_caldyn_Coriolis_mod, ONLY : compute_caldyn_Coriolis
Index: /codes/icosagcm/devel/src/dynamics/compute_theta.F90
===================================================================
--- /codes/icosagcm/devel/src/dynamics/compute_theta.F90	(revision 916)
+++ /codes/icosagcm/devel/src/dynamics/compute_theta.F90	(revision 917)
@@ -1,4 +1,5 @@
 MODULE compute_theta_mod
-  USE grid_param, ONLY : llm
+  USE grid_param, ONLY : llm, nqdyn
+  USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass, ptop
   IMPLICIT NONE
   PRIVATE
@@ -6,7 +7,19 @@
 #include "../unstructured/unstructured.h90"
 
-  PUBLIC :: compute_theta
+  PUBLIC :: compute_theta_unst, compute_theta_hex, compute_theta_manual
 
 CONTAINS
+
+  ! Python/Fortran differences to be resolved at some point :
+
+  ! the Fortran driver prognoses ps/mass_col or rhodz depending on caldyn_eta
+  ! the Python driver prognoses rhodz even if caldyn_eta==eta_mass
+  ! so that the DOFs are the same whatever caldyn_eta
+
+  ! in the Fortran driver dak, dbk are 1D and based on pressure
+  ! => m = mass_dak(l)+(mass_col(ij)*g+ptop)*mass_dbk(l)
+  !    rhodz(CELL) = m/g 
+  ! in the Python driver dak, dbk are 2D and based on mass
+  ! => rhodz(CELL) = MASS_DAK(CELL) + mass_col(HIDX(CELL))*MASS_DBK(CELL)
 
 #ifdef BEGIN_DYSL
@@ -14,7 +27,5 @@
 KERNEL(theta)
   IF(caldyn_eta==eta_mass) THEN ! Compute mass
-    ! FIXME : here mass_col is computed from rhodz
-    ! so that the DOFs are the same whatever caldyn_eta
-    ! in DYNAMICO mass_col is prognosed rather than rhodz
+    ! compute mass_col from rhodz
     SEQUENCE_C1
       PROLOGUE(0)
@@ -27,8 +38,4 @@
     FORALL_CELLS_EXT()
       ON_PRIMAL
-        ! FIXME : formula below (used in DYNAMICO) is for dak, dbk based on
-        ! pressure rather than mass
-        !        m = mass_dak(CELL)+(mass_col(HIDX(CELL))*g+ptop)*mass_dbk(CELL)
-        !        rhodz(CELL) = m/g 
         rhodz(CELL) = MASS_DAK(CELL) + mass_col(HIDX(CELL))*MASS_DBK(CELL)
       END_BLOCK
@@ -44,18 +51,34 @@
 END_BLOCK
 
+KERNEL(compute_theta)
+  IF(caldyn_eta==eta_mass) THEN ! Compute mass
+     FORALL_CELLS_EXT()
+       ON_PRIMAL
+         m = MASS_DAK(CELL)+(mass_col(HIDX(CELL))*g+ptop)*MASS_DBK(CELL)
+         rhodz(CELL) = m/g
+        END_BLOCK
+     END_BLOCK
+  END IF  
+  DO iq=1,nqdyn
+     FORALL_CELLS_EXT()
+       ON_PRIMAL
+         theta(CELL,iq) = theta_rhodz(CELL,iq)/rhodz(CELL)
+       END_BLOCK
+     END_BLOCK
+  END DO
+END_BLOCK
+
 #endif END_DYSL
 
-SUBROUTINE compute_theta_unst(mass_col,rhodz,theta_rhodz, theta)
-  USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
-  USE data_unstructured_mod, ONLY : id_theta,primal_num,dual_num,edge_num, &
-    enter_trace, exit_trace
-  USE grid_param, ONLY : nqdyn
-  USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass, ptop
-  FIELD_PS :: mass_col
-  FIELD_MASS :: rhodz
-  FIELD_THETA :: theta, theta_rhodz
-  DECLARE_INDICES
-  NUM :: m
-  START_TRACE(id_theta, 3,0,0) ! primal, dual, edge  
+  SUBROUTINE compute_theta_unst(mass_col,theta_rhodz, rhodz,theta)
+    USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
+    USE data_unstructured_mod, ONLY : id_theta,primal_num,dual_num,edge_num, &
+         enter_trace, exit_trace
+    FIELD_PS :: mass_col
+    FIELD_MASS :: rhodz
+    FIELD_THETA :: theta, theta_rhodz
+    DECLARE_INDICES
+    NUM :: m
+    START_TRACE(id_theta, 3,0,0) ! primal, dual, edge  
 #define MASS_DAK(l,ij) mass_dak(l)
 #define MASS_DBK(l,ij) mass_dbk(l)
@@ -63,12 +86,30 @@
 #undef MASS_DAK
 #undef MASS_DBK
-  STOP_TRACE
-END SUBROUTINE
+    STOP_TRACE
+  END SUBROUTINE compute_theta_unst
 
-  SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta)
+  SUBROUTINE compute_theta_hex(mass_col,theta_rhodz, rhodz,theta)
     USE icosa
     USE trace, ONLY : trace_start, trace_end
     USE omp_para, ONLY : ll_begin, ll_end
-    USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass, ptop
+    REAL(rstd),INTENT(IN)    :: mass_col(iim*jjm)
+    REAL(rstd),INTENT(IN)    :: theta_rhodz(iim*jjm,llm,nqdyn)
+    REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT)   :: theta(iim*jjm,llm,nqdyn)
+    INTEGER :: ij,l,iq
+    REAL(rstd) :: m
+    CALL trace_start("compute_theta")
+#define MASS_DAK(ij,l) mass_dak(l)
+#define MASS_DBK(ij,l) mass_dbk(l)
+#include "../kernels_hex/compute_theta.k90"
+#undef MASS_DAK
+#undef MASS_DBK
+    CALL trace_end("compute_theta")
+  END SUBROUTINE compute_theta_hex
+  
+  SUBROUTINE compute_theta_manual(ps,theta_rhodz, rhodz,theta)
+    USE icosa
+    USE trace, ONLY : trace_start, trace_end
+    USE omp_para, ONLY : ll_begin, ll_end
     REAL(rstd),INTENT(IN)    :: ps(iim*jjm)
     REAL(rstd),INTENT(IN)    :: theta_rhodz(iim*jjm,llm,nqdyn)
@@ -99,5 +140,5 @@
 
     CALL trace_end("compute_theta")
-  END SUBROUTINE compute_theta
+  END SUBROUTINE compute_theta_manual
   
 END MODULE compute_theta_mod
Index: /codes/icosagcm/devel/src/kernels_hex/compute_theta.k90
===================================================================
--- /codes/icosagcm/devel/src/kernels_hex/compute_theta.k90	(revision 917)
+++ /codes/icosagcm/devel/src/kernels_hex/compute_theta.k90	(revision 917)
@@ -0,0 +1,21 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- compute_theta ----------------------------------
+   IF(caldyn_eta==eta_mass) THEN ! Compute mass
+      DO l = ll_begin, ll_end
+         !DIR$ SIMD
+         DO ij=ij_begin_ext, ij_end_ext
+            m = MASS_DAK(ij,l)+(mass_col(ij)*g+ptop)*MASS_DBK(ij,l)
+            rhodz(ij,l) = m/g
+         END DO
+      END DO
+   END IF
+   DO iq=1,nqdyn
+      DO l = ll_begin, ll_end
+         !DIR$ SIMD
+         DO ij=ij_begin_ext, ij_end_ext
+            theta(ij,l,iq) = theta_rhodz(ij,l,iq)/rhodz(ij,l)
+         END DO
+      END DO
+   END DO
+   !---------------------------- compute_theta ----------------------------------
+   !--------------------------------------------------------------------------
Index: /codes/icosagcm/devel/src/kernels_unst/theta.k90
===================================================================
--- /codes/icosagcm/devel/src/kernels_unst/theta.k90	(revision 916)
+++ /codes/icosagcm/devel/src/kernels_unst/theta.k90	(revision 917)
@@ -2,7 +2,5 @@
    !---------------------------- theta ----------------------------------
    IF(caldyn_eta==eta_mass) THEN ! Compute mass
-      ! FIXME : here mass_col is computed from rhodz
-      ! so that the DOFs are the same whatever caldyn_eta
-      ! in DYNAMICO mass_col is prognosed rather than rhodz
+      ! compute mass_col from rhodz
       !$OMP DO SCHEDULE(STATIC)
       DO ij=1,primal_num
@@ -17,8 +15,4 @@
          !DIR$ SIMD
          DO l = 1, llm
-            ! FIXME : formula below (used in DYNAMICO) is for dak, dbk based on
-            ! pressure rather than mass
-            ! m = mass_dak(l,ij)+(mass_col(ij)*g+ptop)*mass_dbk(l,ij)
-            ! rhodz(l,ij) = m/g
             rhodz(l,ij) = MASS_DAK(l,ij) + mass_col(ij)*MASS_DBK(l,ij)
          END DO
