Index: codes/icosagcm/devel/src/base/init_grid_param.f90
===================================================================
--- codes/icosagcm/devel/src/base/init_grid_param.f90	(revision 905)
+++ codes/icosagcm/devel/src/base/init_grid_param.f90	(revision 906)
@@ -1,4 +1,15 @@
 MODULE init_grid_param_mod
+
+  USE compute_diagnostics_mod
+  USE compute_rhodz_mod
+
+  USE compute_caldyn_mod
+  USE compute_pvort_only_mod
+
   IMPLICIT NONE
+  PRIVATE
+  SAVE
+  
+  PUBLIC :: init_grid_param
 
 CONTAINS
@@ -20,4 +31,5 @@
        jjm_glo=iim_glo
        IF(is_mpi_root) PRINT *, 'GETIN nbp = ',iim_glo
+       CALL select_compute_hex
     CASE('unstructured')
        grid_type = grid_unst
@@ -26,4 +38,5 @@
        CALL open_local_mesh_file
        IF (is_mpi_root) PRINT *,'DYNAMICO mesh is unstructured/LAM.'
+       CALL select_compute_unst
     CASE DEFAULT
        PRINT *, 'Invalid value of grid_type :',TRIM(grid_type_var)
@@ -41,5 +54,16 @@
     
   END SUBROUTINE  init_grid_param
+
+
+  SUBROUTINE select_compute_hex
+    compute_rhodz      => compute_rhodz_hex
+    compute_pvort_only => compute_pvort_only_hex
+  END SUBROUTINE select_compute_hex
+
+  SUBROUTINE select_compute_unst
+    compute_rhodz      => compute_rhodz_unst
+    compute_pvort_only => compute_pvort_only_unst
+  END SUBROUTINE select_compute_unst
+
   
 END MODULE init_grid_param_mod
-    
Index: codes/icosagcm/devel/src/diagnostics/compute_diagnostics.f90
===================================================================
--- codes/icosagcm/devel/src/diagnostics/compute_diagnostics.f90	(revision 906)
+++ codes/icosagcm/devel/src/diagnostics/compute_diagnostics.f90	(revision 906)
@@ -0,0 +1,18 @@
+MODULE compute_diagnostics_mod
+  IMPLICIT NONE
+  SAVE
+
+  INTERFACE
+
+    SUBROUTINE comp_rhodz(flag,ps,rhodz)
+      USE prec, ONLY : rstd
+      LOGICAL, INTENT(IN) :: flag
+      REAL(rstd),INTENT(IN)  :: ps(1)
+      REAL(rstd),INTENT(OUT) :: rhodz(1,1)
+    END SUBROUTINE comp_rhodz
+
+  END INTERFACE
+
+  PROCEDURE(comp_rhodz), POINTER :: compute_rhodz
+
+END MODULE compute_diagnostics_mod
Index: codes/icosagcm/devel/src/diagnostics/compute_rhodz.F90
===================================================================
--- codes/icosagcm/devel/src/diagnostics/compute_rhodz.F90	(revision 906)
+++ codes/icosagcm/devel/src/diagnostics/compute_rhodz.F90	(revision 906)
@@ -0,0 +1,120 @@
+MODULE compute_rhodz_mod
+  USE earth_const, ONLY : g
+  USE disvert_mod, ONLY : ap, bp
+  USE grid_param, ONLY : llm
+  USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
+  IMPLICIT NONE
+  PRIVATE
+  SAVE
+
+#include "../unstructured/unstructured.h90"
+
+  PUBLIC :: compute_rhodz_hex, compute_rhodz_unst
+
+CONTAINS
+
+  SUBROUTINE check_interface
+    USE compute_diagnostics_mod
+    compute_rhodz => compute_rhodz_hex
+    compute_rhodz => compute_rhodz_unst
+  END SUBROUTINE check_interface
+
+#if BEGIN_DYSL
+
+KERNEL(compute_rhodz)
+IF(comp) THEN
+   FORALL_CELLS()
+     ON_PRIMAL
+       m = ( AP(CELL)-AP(UP(CELL)) + (BP(CELL)-BP(UP(CELL)))*ps(HIDX(CELL)) )/g
+       rhodz(CELL)=m
+     END_BLOCK
+   END_BLOCK
+ELSE
+   err=0.
+   FORALL_CELLS_EXT()
+     ON_PRIMAL
+       m = ( AP(CELL)-AP(UP(CELL)) + (BP(CELL)-BP(UP(CELL)))*ps(HIDX(CELL)) )/g
+       err = MAX(err, ABS(m-rhodz(CELL)))
+     END_BLOCK
+   END_BLOCK
+   IF(err>1e-10) THEN
+      PRINT *, 'Discrepancy between ps and rhodz detected', err
+      STOP
+   END IF
+END IF
+END_BLOCK
+
+#endif END_DYSL
+
+  SUBROUTINE compute_rhodz_unst(comp, ps, rhodz)
+    USE data_unstructured_mod, ONLY : primal_num
+    LOGICAL, INTENT(IN)  :: comp
+    FIELD_PS, INTENT(IN) :: ps
+    FIELD_MASS, INTENT(INOUT) :: rhodz
+    DECLARE_INDICES
+    NUM :: m, err
+#define AP(l,ij) ap(l)
+#define BP(l,ij) bp(l)
+#include "../kernels_unst/compute_rhodz.k90"
+#undef AP
+#undef BP
+  END SUBROUTINE compute_rhodz_unst
+  
+  SUBROUTINE compute_rhodz_hex(comp, ps, rhodz)
+    USE icosa
+    USE omp_para
+    LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check
+    REAL(rstd), INTENT(IN) :: ps(iim*jjm)
+    REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm)
+    REAL(rstd) :: m, err
+    INTEGER :: ij,l
+#define AP(ij,l) ap(l)
+#define BP(ij,l) bp(l)
+#include "../kernels_hex/compute_rhodz.k90"
+#undef AP
+#undef BP
+  END SUBROUTINE compute_rhodz_hex
+
+  SUBROUTINE compute_rhodz_handmade(comp, ps, rhodz)
+    USE icosa
+    USE omp_para
+    LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check
+    REAL(rstd), INTENT(IN) :: ps(iim*jjm)
+    REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm)
+    REAL(rstd) :: m, err
+    INTEGER :: l,i,j,ij,dd
+    err=0.
+
+    IF(comp) THEN
+       dd=1
+    ELSE
+       dd=0
+    END IF
+
+    DO l = ll_begin, ll_end
+       DO j=jj_begin-dd,jj_end+dd
+          DO i=ii_begin-dd,ii_end+dd
+             ij=(j-1)*iim+i
+             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
+             IF(comp) THEN
+                rhodz(ij,l) = m
+             ELSE
+                err = MAX(err,abs(m-rhodz(ij,l)))
+             END IF
+          ENDDO
+       ENDDO
+    ENDDO
+
+    IF(.NOT. comp) THEN
+       IF(err>1e-10) THEN
+          PRINT *, 'Discrepancy between ps and rhodz detected', err
+          STOP
+       ELSE
+!          PRINT *, 'No discrepancy between ps and rhodz detected'
+       END IF
+    END IF
+
+  END SUBROUTINE compute_rhodz_handmade
+ 
+
+END MODULE compute_rhodz_mod
Index: codes/icosagcm/devel/src/diagnostics/observable.f90
===================================================================
--- codes/icosagcm/devel/src/diagnostics/observable.f90	(revision 905)
+++ codes/icosagcm/devel/src/diagnostics/observable.f90	(revision 906)
@@ -58,6 +58,4 @@
 !    IF (is_master) PRINT *,'CALL write_output_fields_basic'
 
-    IF(grid_type==grid_unst) RETURN ! FIXME
-
     CALL transfert_request(f_ps,req_i1)
     
@@ -102,6 +100,18 @@
 
     IF(init) THEN
+       CALL output_field("p_init",f_pmid)
+       CALL output_field("ps_init",f_ps)
+       CALL output_field("mass_init",f_mass)
+       CALL output_field("geopot_init",f_geopot)
+       CALL output_field("q_init",f_q)
+
        CALL output_field("temp_init",f_buf_i)
     ELSE
+       CALL output_field("p",f_pmid)
+       CALL output_field("ps",f_ps)
+       CALL output_field("mass",f_mass)
+       CALL output_field("geopot",f_geopot)
+       CALL output_field("q",f_q)
+
        CALL output_field("temp",f_buf_i)
        CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,85000.)
@@ -113,4 +123,6 @@
     END IF
     
+    IF(grid_type == grid_unst) RETURN
+
     CALL progonostic_vel_to_horiz(f_geopot, f_ps, f_mass, f_u, f_W, f_buf_uh, f_buf_i)
     CALL transfert_request(f_buf_uh,req_e1_vect) 
@@ -120,9 +132,4 @@
        CALL output_field("ulon_init",f_buf_ulon)
        CALL output_field("ulat_init",f_buf_ulat)
-       CALL output_field("p_init",f_pmid)
-       CALL output_field("ps_init",f_ps)
-       CALL output_field("mass_init",f_mass)
-       CALL output_field("geopot_init",f_geopot)
-       CALL output_field("q_init",f_q)
        IF(is_master) PRINT *, 'Done writing initial condition ...'
     ELSE
@@ -130,9 +137,4 @@
        CALL output_field("ulon",f_buf_ulon)
        CALL output_field("ulat",f_buf_ulat)
-       CALL output_field("p",f_pmid)
-       CALL output_field("ps",f_ps)
-       CALL output_field("mass",f_mass)
-       CALL output_field("geopot",f_geopot)
-       CALL output_field("q",f_q)
 
        !       CALL output_field("exner",f_pk)
@@ -211,5 +213,6 @@
 
   SUBROUTINE progonostic_vel_to_horiz(f_geopot, f_ps, f_rhodz, f_u, f_W, f_uh, f_uz)
-    USE disvert_mod
+    USE disvert_mod, ONLY : caldyn_eta, eta_mass
+    USE compute_diagnostics_mod, ONLY : compute_rhodz
     TYPE(t_field), POINTER :: f_geopot(:), f_ps(:), f_rhodz(:), &
          f_u(:), f_W(:), f_uz(:), &  ! IN
@@ -294,5 +297,5 @@
     USE pression_mod
     IMPLICIT NONE
-    TYPE(t_field), POINTER :: f_pmid(:)           ! IN
+    TYPE(t_field), POINTER :: f_pmid(:)         ! IN
     TYPE(t_field), POINTER :: f_q(:)            ! IN
     TYPE(t_field), POINTER :: f_temp(:)         ! INOUT
@@ -329,5 +332,5 @@
           p_ik = pmid(ij,l)
           theta_ik = temp(ij,l)
-          qv = q(ij,l,1) ! water vaper mixing ratio = mv/md
+          qv = q(ij,l,1) ! water vapor mixing ratio = mv/md
           SELECT CASE(caldyn_thermo)
           CASE(thermo_theta)
Index: codes/icosagcm/devel/src/dynamics/caldyn_gcm.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_gcm.F90	(revision 905)
+++ codes/icosagcm/devel/src/dynamics/caldyn_gcm.F90	(revision 906)
@@ -16,6 +16,4 @@
   
   SUBROUTINE init_caldyn
-    USE compute_caldyn_mod, ONLY : compute_pvort_only
-    USE compute_pvort_only_mod, ONLY : compute_pvort_only_hex, compute_pvort_only_unst
 
     CHARACTER(len=255) :: def
@@ -160,7 +158,5 @@
           CALL compute_planetvel(planetvel)
        END DO
-       compute_pvort_only => compute_pvort_only_hex
     CASE(grid_unst)
-       compute_pvort_only => compute_pvort_only_unst
     END SELECT
 
Index: codes/icosagcm/devel/src/initial/etat0.f90
===================================================================
--- codes/icosagcm/devel/src/initial/etat0.f90	(revision 905)
+++ codes/icosagcm/devel/src/initial/etat0.f90	(revision 906)
@@ -63,5 +63,6 @@
     USE omp_para, ONLY : is_omp_level_master
     USE mpipara, ONLY : is_mpi_root
-    USE disvert_mod
+    USE disvert_mod, ONLY : caldyn_eta, eta_lag
+    USE compute_diagnostics_mod, ONLY : compute_rhodz
     ! Generic interface
     USE etat0_isothermal_mod, ONLY : getin_etat0_isothermal=>getin_etat0
Index: codes/icosagcm/devel/src/kernels_hex/compute_rhodz.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_hex/compute_rhodz.k90	(revision 906)
+++ codes/icosagcm/devel/src/kernels_hex/compute_rhodz.k90	(revision 906)
@@ -0,0 +1,26 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- compute_rhodz ----------------------------------
+   IF(comp) THEN
+      DO l = ll_begin, ll_end
+         !DIR$ SIMD
+         DO ij=ij_begin, ij_end
+            m = ( AP(ij,l)-AP(ij,l+1) + (BP(ij,l)-BP(ij,l+1))*ps(ij) )/g
+            rhodz(ij,l)=m
+         END DO
+      END DO
+   ELSE
+      err=0.
+      DO l = ll_begin, ll_end
+         !DIR$ SIMD
+         DO ij=ij_begin_ext, ij_end_ext
+            m = ( AP(ij,l)-AP(ij,l+1) + (BP(ij,l)-BP(ij,l+1))*ps(ij) )/g
+            err = MAX(err, ABS(m-rhodz(ij,l)))
+         END DO
+      END DO
+      IF(err>1e-10) THEN
+         PRINT *, 'Discrepancy between ps and rhodz detected', err
+         STOP
+      END IF
+   END IF
+   !---------------------------- compute_rhodz ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/compute_rhodz.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/compute_rhodz.k90	(revision 906)
+++ codes/icosagcm/devel/src/kernels_unst/compute_rhodz.k90	(revision 906)
@@ -0,0 +1,30 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- compute_rhodz ----------------------------------
+   IF(comp) THEN
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij = 1, primal_num
+         !DIR$ SIMD
+         DO l = 1, llm
+            m = ( AP(l,ij)-AP(l+1,ij) + (BP(l,ij)-BP(l+1,ij))*ps(ij) )/g
+            rhodz(l,ij)=m
+         END DO
+      END DO
+      !$OMP END DO
+   ELSE
+      err=0.
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij = 1, primal_num
+         !DIR$ SIMD
+         DO l = 1, llm
+            m = ( AP(l,ij)-AP(l+1,ij) + (BP(l,ij)-BP(l+1,ij))*ps(ij) )/g
+            err = MAX(err, ABS(m-rhodz(l,ij)))
+         END DO
+      END DO
+      !$OMP END DO
+      IF(err>1e-10) THEN
+         PRINT *, 'Discrepancy between ps and rhodz detected', err
+         STOP
+      END IF
+   END IF
+   !---------------------------- compute_rhodz ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/time/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/devel/src/time/timeloop_gcm.f90	(revision 905)
+++ codes/icosagcm/devel/src/time/timeloop_gcm.f90	(revision 906)
@@ -5,4 +5,5 @@
   USE trace
   USE omp_para
+  USE compute_diagnostics_mod, ONLY : compute_rhodz
   USE euler_scheme_mod
   USE explicit_scheme_mod
@@ -247,4 +248,6 @@
     !$OMP END MASTER   
 
+    IF(grid_type == grid_unst) RETURN
+
     DO it=itau0+1,itau0+itaumax
        IF (is_master) CALL print_iteration(it, itau0, itaumax, start_clock, rate_clock)
Index: codes/icosagcm/devel/src/vertical/disvert.f90
===================================================================
--- codes/icosagcm/devel/src/vertical/disvert.f90	(revision 905)
+++ codes/icosagcm/devel/src/vertical/disvert.f90	(revision 906)
@@ -148,46 +148,4 @@
   END SUBROUTINE init_disvert  
   
-  SUBROUTINE compute_rhodz(comp, ps, rhodz)
-    USE icosa
-    USE omp_para
-    LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check
-    REAL(rstd), INTENT(IN) :: ps(iim*jjm)
-    REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm)
-    REAL(rstd) :: m, err
-    INTEGER :: l,i,j,ij,dd
-    err=0.
-
-    IF(comp) THEN
-       dd=1
-    ELSE
-       dd=0
-    END IF
-
-    DO l = ll_begin, ll_end
-       DO j=jj_begin-dd,jj_end+dd
-          DO i=ii_begin-dd,ii_end+dd
-             ij=(j-1)*iim+i
-             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
-             IF(comp) THEN
-                rhodz(ij,l) = m
-             ELSE
-                err = MAX(err,abs(m-rhodz(ij,l)))
-             END IF
-          ENDDO
-       ENDDO
-    ENDDO
-
-    IF(.NOT. comp) THEN
-       IF(err>1e-10) THEN
-          PRINT *, 'Discrepancy between ps and rhodz detected', err
-          STOP
-       ELSE
-!          PRINT *, 'No discrepancy between ps and rhodz detected'
-       END IF
-    END IF
-
-  END SUBROUTINE compute_rhodz
-
- 
   SUBROUTINE write_apbp
   USE icosa
