Index: codes/icosagcm/devel/src/diagnostics/check_conserve.f90
===================================================================
--- codes/icosagcm/devel/src/diagnostics/check_conserve.f90	(revision 1040)
+++ codes/icosagcm/devel/src/diagnostics/check_conserve.f90	(revision 1052)
@@ -52,5 +52,5 @@
 
   SUBROUTINE check_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it) 
-  USE vorticity_mod
+  USE compute_vorticity_mod
   USE caldyn_gcm_mod
   USE exner_mod 
@@ -153,8 +153,6 @@
 
   SUBROUTINE check_conserve_detailed(it,tag, f_ps,f_dps,f_ue,f_theta_rhodz,f_phis) 
-    USE vorticity_mod
     USE caldyn_gcm_mod
     USE exner_mod 
-!    USE mpipara, ONLY : is_mpi_root, comm_icosa
     USE omp_para, ONLY : is_master
     INTEGER, INTENT(IN) :: tag
@@ -268,5 +266,4 @@
   SUBROUTINE check_energy(f_ue,f_theta_rhodz,f_phis, etot, &
                       stot, AAM_mass_tot, AAM_vel_tot, AAM_velp_tot, AAM_velm_tot, rmsvtot)  
-  USE vorticity_mod
     TYPE(t_field), POINTER :: f_ue(:)
     TYPE(t_field), POINTER :: f_theta_rhodz(:)
Index: codes/icosagcm/devel/src/diagnostics/compute_diagnostics.f90
===================================================================
--- codes/icosagcm/devel/src/diagnostics/compute_diagnostics.f90	(revision 1040)
+++ codes/icosagcm/devel/src/diagnostics/compute_diagnostics.f90	(revision 1052)
@@ -27,4 +27,16 @@
     END SUBROUTINE comp_temperature
 
+    SUBROUTINE comp_vorticity(ue, vort)
+      IMPORT
+      REAL(rstd),INTENT(IN)  :: ue(:,:)
+      REAL(rstd),INTENT(OUT) :: vort(:,:)
+    END SUBROUTINE comp_vorticity
+
+    SUBROUTINE comp_divergence(ue, div)
+      IMPORT
+      REAL(rstd),INTENT(IN)  :: ue(:,:)
+      REAL(rstd),INTENT(OUT) :: div(:,:)
+    END SUBROUTINE comp_divergence
+
     SUBROUTINE comp_hydro_press(rhodz, theta_rhodz, ps, p)
       IMPORT
@@ -48,9 +60,11 @@
   PROCEDURE(comp_pression),    POINTER, SAVE :: compute_pression => NULL(), compute_pression_mid => NULL()
   PROCEDURE(comp_temperature), POINTER, SAVE :: compute_temperature => NULL()
+  PROCEDURE(comp_vorticity),   POINTER, SAVE :: compute_vorticity => NULL()
+  PROCEDURE(comp_divergence),  POINTER, SAVE :: compute_divergence => NULL()
   PROCEDURE(comp_hydro_press), POINTER, SAVE :: compute_hydrostatic_pressure => NULL()
   PROCEDURE(comp_vert_interp), POINTER, SAVE :: compute_vertical_interp => NULL()
 
-  PUBLIC :: compute_rhodz, compute_pression, compute_pression_mid, compute_temperature, compute_hydrostatic_pressure, &
-       compute_vertical_interp
+  PUBLIC :: compute_rhodz, compute_pression, compute_pression_mid, compute_temperature, &
+       compute_vorticity, compute_divergence, compute_hydrostatic_pressure, compute_vertical_interp
 
 END MODULE compute_diagnostics_mod
Index: codes/icosagcm/devel/src/diagnostics/compute_divergence.F90
===================================================================
--- codes/icosagcm/devel/src/diagnostics/compute_divergence.F90	(revision 1052)
+++ codes/icosagcm/devel/src/diagnostics/compute_divergence.F90	(revision 1052)
@@ -0,0 +1,94 @@
+MODULE compute_divergence_mod
+  USE grid_param
+  USE prec, ONLY : rstd
+  IMPLICIT NONE
+  PRIVATE
+
+#include "../unstructured/unstructured.h90"
+
+  PUBLIC :: divergence, compute_divergence_hex, compute_divergence_unst
+
+CONTAINS
+
+#ifdef BEGIN_DYSL
+
+KERNEL(compute_divergence)
+
+  FORALL_CELLS_EXT()
+    ON_PRIMAL
+      div_ij = 0.d0
+      FORALL_EDGES
+         div_ij = div_ij + SIGN*ue(EDGE)*LE
+      END_BLOCK
+      div(CELL) = div_ij / AI
+    END_BLOCK
+  END_BLOCK
+
+END_BLOCK
+
+#endif END_DYSL
+
+  SUBROUTINE divergence(f_ue,f_div)
+    USE icosa
+    USE compute_diagnostics_mod, ONLY : compute_divergence
+    TYPE(t_field), POINTER :: f_ue(:)
+    TYPE(t_field), POINTER :: f_div(:)
+  
+    REAL(rstd), POINTER :: ue(:,:)
+    REAL(rstd), POINTER :: div(:,:)
+    INTEGER :: ind
+
+    CALL transfert_request(f_ue,req_e1_vect)
+   
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind)) CYCLE
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      ue=f_ue(ind)
+      div=f_div(ind)
+      CALL compute_divergence(ue, div)
+    ENDDO
+    
+  END SUBROUTINE divergence
+
+!-------------- Wrappers for F2008 conformity -----------------
+
+  SUBROUTINE compute_divergence_unst(ue,div)
+    REAL(rstd),INTENT(IN)  :: ue(:,:)
+    REAL(rstd),INTENT(OUT) :: div(:,:)
+    CALL compute_divergence_unst_(ue,div)
+  END SUBROUTINE compute_divergence_unst
+
+  SUBROUTINE compute_divergence_hex(ue,div)
+    REAL(rstd),INTENT(IN)  :: ue(:,:)
+    REAL(rstd),INTENT(OUT) :: div(:,:)
+    CALL compute_divergence_hex_(ue,div)
+  END SUBROUTINE compute_divergence_hex
+
+!--------------------------------------------------------------
+
+  SUBROUTINE compute_divergence_unst_(ue,div)
+    USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
+    USE geometry, ONLY : le, Ai
+    USE data_unstructured_mod, ONLY : primal_deg, primal_edge, primal_ne
+    FIELD_U    :: ue
+    FIELD_MASS :: div
+    DECLARE_INDICES
+    DECLARE_EDGES
+    DECLARE_VERTICES
+    NUM :: div_ij
+#include "../kernels_unst/compute_divergence.k90"
+  END SUBROUTINE compute_divergence_unst_
+
+  SUBROUTINE compute_divergence_hex_(ue,div)
+    USE icosa
+    USE omp_para
+    REAL(rstd),INTENT(IN)  :: ue(3*iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: div(iim*jjm,llm)
+    REAL(rstd) :: div_ij
+    INTEGER :: ij,l
+
+#include "../kernels_hex/compute_divergence.k90"
+  END SUBROUTINE compute_divergence_hex_
+
+END MODULE compute_divergence_mod
Index: codes/icosagcm/devel/src/diagnostics/compute_vorticity.F90
===================================================================
--- codes/icosagcm/devel/src/diagnostics/compute_vorticity.F90	(revision 1052)
+++ codes/icosagcm/devel/src/diagnostics/compute_vorticity.F90	(revision 1052)
@@ -0,0 +1,123 @@
+MODULE compute_vorticity_mod
+  USE grid_param
+  USE prec, ONLY : rstd
+  IMPLICIT NONE
+  PRIVATE
+
+#include "../unstructured/unstructured.h90"
+
+  PUBLIC :: vorticity, compute_vorticity_manual, compute_vorticity_hex, compute_vorticity_unst
+
+CONTAINS
+
+#ifdef BEGIN_DYSL
+
+KERNEL(compute_vorticity)
+
+  FORALL_CELLS_EXT()
+    ON_DUAL
+      etav = 0.d0
+      FORALL_EDGES
+         etav = etav + SIGN*ue(EDGE)*DE
+      END_BLOCK
+      vort(DUAL_CELL) = etav / AV
+    END_BLOCK
+
+  END_BLOCK
+
+END_BLOCK
+
+#endif END_DYSL
+
+  SUBROUTINE vorticity(f_ue,f_vort)
+    USE icosa
+    USE compute_diagnostics_mod, ONLY : compute_vorticity
+    TYPE(t_field), POINTER :: f_ue(:)
+    TYPE(t_field), POINTER :: f_vort(:)
+  
+    REAL(rstd), POINTER :: ue(:,:)
+    REAL(rstd), POINTER :: vort(:,:)
+    INTEGER :: ind
+
+    CALL transfert_request(f_ue,req_e1_vect)
+   
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind)) CYCLE
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      ue=f_ue(ind)
+      vort=f_vort(ind)
+      CALL compute_vorticity(ue, vort)
+    ENDDO
+  
+  END SUBROUTINE vorticity
+
+!-------------- Wrappers for F2008 conformity -----------------
+
+  SUBROUTINE compute_vorticity_unst(ue,vort)
+    REAL(rstd),INTENT(IN)  :: ue(:,:)
+    REAL(rstd),INTENT(OUT) :: vort(:,:)
+    CALL compute_vorticity_unst_(ue,vort)
+  END SUBROUTINE compute_vorticity_unst
+
+  SUBROUTINE compute_vorticity_hex(ue,vort)
+    REAL(rstd),INTENT(IN)  :: ue(:,:)
+    REAL(rstd),INTENT(OUT) :: vort(:,:)
+    CALL compute_vorticity_hex_(ue,vort)
+  END SUBROUTINE compute_vorticity_hex
+
+!--------------------------------------------------------------
+
+  SUBROUTINE compute_vorticity_unst_(ue,vort)
+    USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
+    USE geometry, ONLY : de, Av
+    USE data_unstructured_mod, ONLY : dual_deg, dual_edge, dual_ne
+    FIELD_U    :: ue
+    FIELD_Z    :: vort
+    DECLARE_INDICES
+    DECLARE_EDGES
+    DECLARE_VERTICES
+    NUM :: etav
+#include "../kernels_unst/compute_vorticity.k90"
+  END SUBROUTINE compute_vorticity_unst_
+
+  SUBROUTINE compute_vorticity_hex_(ue,vort)
+    USE icosa
+    USE omp_para
+    REAL(rstd),INTENT(IN)  :: ue(3*iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: vort(2*iim*jjm,llm)
+    REAL(rstd) :: etav
+    INTEGER :: ij,l
+
+#include "../kernels_hex/compute_vorticity.k90"
+  END SUBROUTINE compute_vorticity_hex_
+
+  SUBROUTINE compute_vorticity_manual(ue,vort)
+  USE icosa
+  USE disvert_mod
+  USE omp_para
+  IMPLICIT NONE
+    REAL(rstd),INTENT(IN)  :: ue(3*iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: vort(2*iim*jjm,llm)
+    INTEGER :: i,j,ij,l
+
+    DO l = ll_begin,ll_end
+      DO j=jj_begin-1,jj_end+1
+        DO i=ii_begin-1,ii_end+1
+          ij=(j-1)*iim+i
+        
+           vort(ij+z_up,l) = 1./Av(ij+z_up)*(  ne(ij,rup)        * ue(ij+u_rup,l)        * de(ij+u_rup)         &
+                                + ne(ij+t_rup,left) * ue(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)               &
+                                - ne(ij,lup)        * ue(ij+u_lup,l)        * de(ij+u_lup) )                               
+
+           vort(ij+z_down,l) = 1./Av(ij+z_down)*(  ne(ij,ldown)         * ue(ij+u_ldown,l)         * de(ij+u_ldown)          &
+                                   + ne(ij+t_ldown,right) * ue(ij+t_ldown+u_right,l) * de(ij+t_ldown+u_right)                &
+                                   - ne(ij,rdown)         * ue(ij+u_rdown,l)          * de(ij+u_rdown) )
+
+        ENDDO
+      ENDDO
+    ENDDO
+    
+  END SUBROUTINE compute_vorticity_manual
+
+END MODULE compute_vorticity_mod
Index: codes/icosagcm/devel/src/diagnostics/observable.f90
===================================================================
--- codes/icosagcm/devel/src/diagnostics/observable.f90	(revision 1040)
+++ codes/icosagcm/devel/src/diagnostics/observable.f90	(revision 1052)
@@ -46,4 +46,6 @@
     USE compute_temperature_mod
     USE compute_velocity_mod
+    USE compute_vorticity_mod
+    USE compute_divergence_mod
     USE vertical_interp_mod
     USE theta2theta_rhodz_mod
@@ -109,4 +111,5 @@
 
        CALL output_field("temp_init",f_buf_i)
+       CALL output_field("vort_init",f_buf_v)
     ELSE
        CALL output_field("p",f_pmid)
@@ -123,4 +126,14 @@
        CALL vertical_interp(f_pmid,f_buf_i,f_buf_s,preff)
        CALL output_field("SST",f_buf_s)       
+    END IF
+
+    CALL vorticity(f_u, f_buf_v)
+    CALL divergence(f_u, f_buf_i)
+    IF(init) THEN
+       CALL output_field("vort_init",f_buf_v)
+       CALL output_field("div_init",f_buf_i)
+    ELSE
+       CALL output_field("vort",f_buf_v)
+       CALL output_field("div",f_buf_i)
     END IF
     
