Index: /codes/icosagcm/devel/Python/src/kernels_transport.jin
===================================================================
--- /codes/icosagcm/devel/Python/src/kernels_transport.jin	(revision 791)
+++ /codes/icosagcm/devel/Python/src/kernels_transport.jin	(revision 792)
@@ -312,5 +312,4 @@
 
 KERNEL(scalar_laplacian)
-
   FORALL_CELLS_EXT()
     ON_EDGES
@@ -330,2 +329,21 @@
 
 END_BLOCK
+
+KERNEL(curl_laplacian)
+  FORALL_CELLS_EXT()
+    ON_DUAL
+      etav = 0.d0
+      FORALL_EDGES
+         etav = etav + SIGN*u(EDGE)
+      END_BLOCK
+      qv(DUAL_CELL) = etav/AV
+    END_BLOCK
+  END_BLOCK
+
+  FORALL_CELLS_EXT()
+    ON_EDGES
+       curlcurl(EDGE) = SIGN*(qv(VERTEX1)-qv(VERTEX2))*(1./LE_DE)
+    END_BLOCK
+  END_BLOCK
+
+END_BLOCK
Index: /codes/icosagcm/devel/Python/src/unstructured.pyx
===================================================================
--- /codes/icosagcm/devel/Python/src/unstructured.pyx	(revision 791)
+++ /codes/icosagcm/devel/Python/src/unstructured.pyx	(revision 792)
@@ -78,4 +78,5 @@
                      ['dynamico_caldyn_unstructured', c_num, c_void_p,20],
                      ['dynamico_scalar_laplacian', c_void_p,2],
+                     ['dynamico_curl_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/Python/test/py/Baroclinic_3D_ullrich.py
===================================================================
--- /codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py	(revision 791)
+++ /codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py	(revision 792)
@@ -201,11 +201,15 @@
                 caldyn_step.next()
                 davies.relax(llm, caldyn_step, flow0)
-                m,S = caldyn_step.mass, caldyn_step.theta_rhodz
+                m,S,u = caldyn_step.mass, caldyn_step.theta_rhodz, caldyn_step.u
                 s = S/m
-                laps = mesh.field_mass()
-                bilaps = mesh.field_mass()
+                laps, bilaps = mesh.field_mass(), mesh.field_mass()
+                lapu, bilapu = mesh.field_u(), mesh.field_u()
                 unst.ker.dynamico_scalar_laplacian(s,laps)
                 unst.ker.dynamico_scalar_laplacian(laps,bilaps)
+                unst.ker.dynamico_curl_laplacian(u,lapu)
+                unst.ker.dynamico_curl_laplacian(lapu,bilapu)
                 caldyn_step.theta_rhodz[:] = S - dt*kappa*bilaps*m # Euler step
+                caldyn_step.u[:] = u - dt*kappa*bilapu # Euler step
+
             return (caldyn_step.mass.copy(), caldyn_step.theta_rhodz.copy(), caldyn_step.u.copy(),
                     caldyn_step.geopot.copy(), caldyn_step.W.copy())
@@ -230,5 +234,5 @@
         v = mesh.field_mass() # specific volume (diagnosed)
         
-        for i in range(Nslice):
+        for i in range(Nslice+1):
             context.update_calendar(i+1)
 
Index: /codes/icosagcm/devel/src/kernels_unst/curl_laplacian.k90
===================================================================
--- /codes/icosagcm/devel/src/kernels_unst/curl_laplacian.k90	(revision 792)
+++ /codes/icosagcm/devel/src/kernels_unst/curl_laplacian.k90	(revision 792)
@@ -0,0 +1,64 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- curl_laplacian ----------------------------------
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, dual_num
+      ! this VLOOP iterates over dual cell edges
+      SELECT CASE(dual_deg(ij))
+      CASE(3)
+         edge1 = dual_edge(1,ij)
+         edge2 = dual_edge(2,ij)
+         edge3 = dual_edge(3,ij)
+         sign1 = dual_ne(1,ij)
+         sign2 = dual_ne(2,ij)
+         sign3 = dual_ne(3,ij)
+         !DIR$ SIMD
+         DO l = 1, llm
+            etav = 0.d0
+            etav = etav + sign1*u(l,edge1)
+            etav = etav + sign2*u(l,edge2)
+            etav = etav + sign3*u(l,edge3)
+            qv(l,ij) = etav/Av(ij)
+         END DO
+      CASE(4)
+         edge1 = dual_edge(1,ij)
+         edge2 = dual_edge(2,ij)
+         edge3 = dual_edge(3,ij)
+         edge4 = dual_edge(4,ij)
+         sign1 = dual_ne(1,ij)
+         sign2 = dual_ne(2,ij)
+         sign3 = dual_ne(3,ij)
+         sign4 = dual_ne(4,ij)
+         !DIR$ SIMD
+         DO l = 1, llm
+            etav = 0.d0
+            etav = etav + sign1*u(l,edge1)
+            etav = etav + sign2*u(l,edge2)
+            etav = etav + sign3*u(l,edge3)
+            etav = etav + sign4*u(l,edge4)
+            qv(l,ij) = etav/Av(ij)
+         END DO
+      CASE DEFAULT
+         !DIR$ SIMD
+         DO l = 1, llm
+            etav = 0.d0
+            DO iedge = 1, dual_deg(ij)
+               edge = dual_edge(iedge,ij)
+               etav = etav + dual_ne(iedge,ij)*u(l,edge)
+            END DO
+            qv(l,ij) = etav/Av(ij)
+         END DO
+      END SELECT
+   END DO
+   !$OMP END DO
+   !$OMP DO SCHEDULE(STATIC)
+   DO edge = 1, edge_num
+      ij_up = up(edge)
+      ij_down = down(edge)
+      !DIR$ SIMD
+      DO l = 1, llm
+         curlcurl(l,edge) = 1.*(qv(l,ij_down)-qv(l,ij_up))*(1./le_de(edge))
+      END DO
+   END DO
+   !$OMP END DO
+   !---------------------------- curl_laplacian ----------------------------------
+   !--------------------------------------------------------------------------
Index: /codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90
===================================================================
--- /codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 791)
+++ /codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 792)
@@ -247,3 +247,15 @@
 END SUBROUTINE compute_scalar_laplacian
 
+SUBROUTINE compute_curl_laplacian(u,qv,curlcurl)
+  FIELD_Z   :: qv ! BUF
+  FIELD_U   :: u, curlcurl ! IN,OUT
+  DECLARE_INDICES
+  DECLARE_EDGES
+  DECLARE_VERTICES
+  NUM :: etav
+  START_TRACE(id_scalar_laplacian, 0,0,1)
+#include "../kernels_unst/curl_laplacian.k90"
+  STOP_TRACE
+END SUBROUTINE compute_curl_laplacian
+
 END MODULE caldyn_unstructured_mod
Index: /codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90
===================================================================
--- /codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90	(revision 791)
+++ /codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90	(revision 792)
@@ -29,4 +29,13 @@
      !$OMP END PARALLEL
    END SUBROUTINE scalar_laplacian
+
+   SUBROUTINE curl_laplacian(u,lapu) BINDC(curl_laplacian)
+     FIELD_U   :: u,lapu ! IN, OUT
+     FIELD_Z   :: qv  ! BUF
+     !$OMP PARALLEL NUM_THREADS(nb_threads)
+     CALL update_halo(transfer_edge, u)
+     CALL compute_curl_laplacian(u,qv,lapu)
+     !$OMP END PARALLEL
+   END SUBROUTINE curl_laplacian
 
   !
