Index: codes/icosagcm/devel/Python/src/kernels_transport.jin
===================================================================
--- codes/icosagcm/devel/Python/src/kernels_transport.jin	(revision 783)
+++ codes/icosagcm/devel/Python/src/kernels_transport.jin	(revision 784)
@@ -310,2 +310,22 @@
     END_BLOCK
 END_BLOCK
+
+KERNEL(scalar_laplacian)
+
+  FORALL_CELLS_EXT()
+    ON_EDGES
+       grad(EDGE) = SIGN*(b(CELL2)-b(CELL1))
+    END_BLOCK
+  END_BLOCK
+
+  FORALL_CELLS_EXT()
+    ON_PRIMAL
+      div_ij=0.
+      FORALL_EDGES
+        div_ij = div_ij + SIGN*LE_DE*grad(EDGE)
+      END_BLOCK
+      divu(CELL) = div_ij / AI
+    END_BLOCK
+  END_BLOCK
+
+END_BLOCK
Index: codes/icosagcm/devel/Python/src/unstructured.pyx
===================================================================
--- codes/icosagcm/devel/Python/src/unstructured.pyx	(revision 783)
+++ codes/icosagcm/devel/Python/src/unstructured.pyx	(revision 784)
@@ -77,4 +77,5 @@
                      ['dynamico_init_hybrid', c_void_p,3],
                      ['dynamico_caldyn_unstructured', c_num, c_void_p,20],
+                     ['dynamico_scalar_laplacian', c_void_p,2],
                      ['dynamico_partition_graph', c_int,2, c_void_p,3, c_int, c_void_p],
                      ['dynamico_init_transfer', c_int, c_int,2,c_void_p,3, c_int,2,c_void_p,3],
Index: codes/icosagcm/devel/src/kernels_unst/scalar_laplacian.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/scalar_laplacian.k90	(revision 784)
+++ codes/icosagcm/devel/src/kernels_unst/scalar_laplacian.k90	(revision 784)
@@ -0,0 +1,109 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- scalar_laplacian ----------------------------------
+   !$OMP DO SCHEDULE(STATIC)
+   DO edge = 1, edge_num
+      ij_left = left(edge)
+      ij_right = right(edge)
+      !DIR$ SIMD
+      DO l = 1, llm
+         grad(l,edge) = 1.*(b(l,ij_right)-b(l,ij_left))
+      END DO
+   END DO
+   !$OMP END DO
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, primal_num
+      ! this VLOOP iterates over primal cell edges
+      SELECT CASE(primal_deg(ij))
+      CASE(4)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         le_de1 = le_de(edge1)
+         le_de2 = le_de(edge2)
+         le_de3 = le_de(edge3)
+         le_de4 = le_de(edge4)
+         sign1 = primal_ne(1,ij)
+         sign2 = primal_ne(2,ij)
+         sign3 = primal_ne(3,ij)
+         sign4 = primal_ne(4,ij)
+         !DIR$ SIMD
+         DO l = 1, llm
+            div_ij=0.
+            div_ij = div_ij + sign1*le_de1*grad(l,edge1)
+            div_ij = div_ij + sign2*le_de2*grad(l,edge2)
+            div_ij = div_ij + sign3*le_de3*grad(l,edge3)
+            div_ij = div_ij + sign4*le_de4*grad(l,edge4)
+            divu(l,ij) = div_ij / Ai(ij)
+         END DO
+      CASE(5)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         edge5 = primal_edge(5,ij)
+         le_de1 = le_de(edge1)
+         le_de2 = le_de(edge2)
+         le_de3 = le_de(edge3)
+         le_de4 = le_de(edge4)
+         le_de5 = le_de(edge5)
+         sign1 = primal_ne(1,ij)
+         sign2 = primal_ne(2,ij)
+         sign3 = primal_ne(3,ij)
+         sign4 = primal_ne(4,ij)
+         sign5 = primal_ne(5,ij)
+         !DIR$ SIMD
+         DO l = 1, llm
+            div_ij=0.
+            div_ij = div_ij + sign1*le_de1*grad(l,edge1)
+            div_ij = div_ij + sign2*le_de2*grad(l,edge2)
+            div_ij = div_ij + sign3*le_de3*grad(l,edge3)
+            div_ij = div_ij + sign4*le_de4*grad(l,edge4)
+            div_ij = div_ij + sign5*le_de5*grad(l,edge5)
+            divu(l,ij) = div_ij / Ai(ij)
+         END DO
+      CASE(6)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         edge5 = primal_edge(5,ij)
+         edge6 = primal_edge(6,ij)
+         le_de1 = le_de(edge1)
+         le_de2 = le_de(edge2)
+         le_de3 = le_de(edge3)
+         le_de4 = le_de(edge4)
+         le_de5 = le_de(edge5)
+         le_de6 = le_de(edge6)
+         sign1 = primal_ne(1,ij)
+         sign2 = primal_ne(2,ij)
+         sign3 = primal_ne(3,ij)
+         sign4 = primal_ne(4,ij)
+         sign5 = primal_ne(5,ij)
+         sign6 = primal_ne(6,ij)
+         !DIR$ SIMD
+         DO l = 1, llm
+            div_ij=0.
+            div_ij = div_ij + sign1*le_de1*grad(l,edge1)
+            div_ij = div_ij + sign2*le_de2*grad(l,edge2)
+            div_ij = div_ij + sign3*le_de3*grad(l,edge3)
+            div_ij = div_ij + sign4*le_de4*grad(l,edge4)
+            div_ij = div_ij + sign5*le_de5*grad(l,edge5)
+            div_ij = div_ij + sign6*le_de6*grad(l,edge6)
+            divu(l,ij) = div_ij / Ai(ij)
+         END DO
+      CASE DEFAULT
+         !DIR$ SIMD
+         DO l = 1, llm
+            div_ij=0.
+            DO iedge = 1, primal_deg(ij)
+               edge = primal_edge(iedge,ij)
+               div_ij = div_ij + primal_ne(iedge,ij)*le_de(edge)*grad(l,edge)
+            END DO
+            divu(l,ij) = div_ij / Ai(ij)
+         END DO
+      END SELECT
+   END DO
+   !$OMP END DO
+   !---------------------------- scalar_laplacian ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90
===================================================================
--- codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 783)
+++ codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 784)
@@ -236,3 +236,14 @@
 END SUBROUTINE
 
+SUBROUTINE compute_scalar_laplacian(b,grad,divu)
+  FIELD_MASS   :: b, divu ! IN, OUT
+  FIELD_U      :: grad ! BUF
+  DECLARE_INDICES
+  DECLARE_EDGES
+  NUM :: div_ij
+  START_TRACE(id_scalar_laplacian, 0,0,1)
+#include "../kernels_unst/scalar_laplacian.k90"
+  STOP_TRACE
+END SUBROUTINE compute_scalar_laplacian
+
 END MODULE caldyn_unstructured_mod
Index: codes/icosagcm/devel/src/unstructured/data_unstructured.F90
===================================================================
--- codes/icosagcm/devel/src/unstructured/data_unstructured.F90	(revision 783)
+++ codes/icosagcm/devel/src/unstructured/data_unstructured.F90	(revision 784)
@@ -46,5 +46,6 @@
   INTEGER, PARAMETER :: id_dev1=1, id_dev2=2, &
        id_pvort_only=3, id_slow_hydro=4, id_fast=5, id_coriolis=6, id_theta=7, id_geopot=8, id_vert=9, &
-       id_solver=10, id_slow_NH=11, id_NH_geopot=12, id_vert_NH=13, id_update=14, id_halo=15, nb_routines=15 
+       id_solver=10, id_slow_NH=11, id_NH_geopot=12, id_vert_NH=13, id_update=14, id_halo=15, &
+       id_scalar_laplacian=16, nb_routines=16 
   TIME, PRIVATE :: start_time, time_spent(nb_routines) ! time spent in each kernel
   INTEGER, PRIVATE :: current_id, nb_calls(nb_routines)
@@ -53,5 +54,5 @@
        (/'dev1      ', 'dev2      ', &
        'pvort_only', 'slow_hydro', 'fast      ', 'coriolis  ', 'theta     ', 'geopot    ', 'vert      ', &
-       'solver    ', 'slow_NH   ', 'NH_geopot ', 'vert_NH   ',  'update    ', 'halo_xchg ' /)
+       'solver    ', 'slow_NH   ', 'NH_geopot ', 'vert_NH   ',  'update    ', 'halo_xchg ', 'scalar_lap' /)
 
   INTEGER, PARAMETER ::transfer_primal=1, transfer_edge=2, transfer_dual=3, transfer_max=3
Index: codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90
===================================================================
--- codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90	(revision 783)
+++ codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90	(revision 784)
@@ -21,4 +21,14 @@
 #define HASNAN(field) (ANY(.NOT.ABS(field)<1e20))
 
+   SUBROUTINE scalar_laplacian(s,laps) BINDC(scalar_laplacian)
+     FIELD_MASS   :: s,laps ! IN, OUT
+     FIELD_U      :: grads  ! BUF
+     !$OMP PARALLEL NUM_THREADS(nb_threads)
+     CALL update_halo(transfer_primal, s)
+     CALL compute_scalar_laplacian(s,grads,laps)
+     !$OMP END PARALLEL
+   END SUBROUTINE scalar_laplacian
+
+  !
   SUBROUTINE caldyn_unstructured(tau, mass_col,rhodz,theta_rhodz,u,geopot,w, & ! IN : flow state
                                  theta,ps,pk,hflux,qv, &    ! OUT : diags (except surface geopot : IN)
