Index: codes/icosagcm/devel/Python/dynamico/LAM.py
===================================================================
--- codes/icosagcm/devel/Python/dynamico/LAM.py	(revision 791)
+++ codes/icosagcm/devel/Python/dynamico/LAM.py	(revision 791)
@@ -0,0 +1,28 @@
+from numpy import pi, cos, zeros
+
+class Davies:
+    def __init__(self,N1,N2,x_i,y_i,x_e,y_e):
+        self.N1, self.N2 = N1, N2
+        self.beta_i = self.mask(x_i,y_i)
+        self.beta_e = self.mask(x_e,y_e)
+        self.x_e,self.y_e = x_e,y_e
+    def mask0(self,c,c0,delta): # 1D building block for Davies relaxation
+        N1, N2 = self.N1, self.N2
+        m = zeros(c.size)
+        c1,c2 = c0-N1*delta, c0-(N1+N2)*delta
+        for i in range(c.size):
+            ci=c[i]
+            m[i] = (1.+cos((ci-c2)*pi/(N2*delta)))/2.0
+            if ci < c2 : m[i]=1.
+            if ci > c1 : m[i]=0.
+        return m
+    def relax(self, llm, step, flow):
+        beta_i, beta_e = self.beta_i, self.beta_e
+        m,S,u,Phi,W=flow
+        for l in range(llm):
+            step.mass[:,l]        = beta_i*step.mass[:,l]        + (1.-beta_i)*m[:,l]
+            step.theta_rhodz[:,l] = beta_i*step.theta_rhodz[:,l] + (1.-beta_i)*S[:,l]
+            step.u[:,l]           = beta_e*step.u[:,l]           + (1.-beta_e)*u[:,l]
+        for l in range(llm+1):
+            step.geopot[:,l]      = beta_i*step.geopot[:,l]      + (1.-beta_i)*Phi[:,l]
+            step.W[:,l]           = beta_i*step.W[:,l]           + (1.-beta_i)*W[:,l]
Index: codes/icosagcm/devel/Python/dynamico/kernels.py
===================================================================
--- codes/icosagcm/devel/Python/dynamico/kernels.py	(revision 791)
+++ codes/icosagcm/devel/Python/dynamico/kernels.py	(revision 791)
@@ -0,0 +1,49 @@
+from dynamico.precision import zeros
+
+def grad(mesh, s_ik):
+    llm=s_ik.shape[1]
+    edge_num, left, right = mesh.edge_num, mesh.left, mesh.right
+    grad_ek = zeros((edge_num, llm))
+    for e in range(edge_num):
+        left_e = left[e]-1 # Fortran => Python indexing
+        right_e = right[e]-1 # Fortran => Python indexing
+        grad_ek[e,:] = s_ik[right_e,:]-s_ik[left_e,:]
+    return grad_ek
+
+def div(mesh, u_ek):
+    llm=u_ek.shape[1]
+    primal_num, primal_deg, primal_edge = mesh.primal_num, mesh.primal_deg, mesh.primal_edge
+    primal_ne, le_de, Ai = mesh.primal_ne, mesh.le_de, mesh.Ai
+    div_ik = zeros((primal_num, llm))
+    for i in range(primal_num):
+        div_i = 0.
+        for iedge in range(primal_deg[i]):
+            e = primal_edge[i,iedge]-1 # Fortran => Python indexing
+            div_i = div_i + u_ek[e,:]*primal_ne[i,iedge]*le_de[e]
+        div_ik[i,:] = div_i / Ai[i]
+    return div_ik
+
+def KE(mesh, u_ek):
+    llm=u_ek.shape[1]
+    primal_num, primal_deg, primal_edge = mesh.primal_num, mesh.primal_deg, mesh.primal_edge
+    primal_ne, le_de, Ai = mesh.primal_ne, mesh.le_de, mesh.Ai
+    div_ik = zeros((primal_num, llm))
+    for i in range(primal_num):
+        div_i = 0.
+        for iedge in range(primal_deg[i]):
+            e = primal_edge[i,iedge]-1 # Fortran => Python indexing
+            div_i = div_i + le_de[e]*(u_ek[e,:]**2)
+        div_ik[i,:] = div_i * (.25/ Ai[i])
+    return div_ik
+
+def curl(mesh, u_ek):
+    llm=u_ek.shape[1]
+    dual_num, dual_deg, dual_edge, dual_ne, Av = mesh.dual_num, mesh.dual_deg, mesh.dual_edge, mesh.dual_ne, mesh.Av
+    curl_vk = zeros((dual_num, llm))
+    for v in range(dual_num):
+        curl_v = 0.
+        for iedge in range(dual_deg[v]):
+            e = dual_edge[v,iedge]-1 # Fortran => Python indexing
+            curl_v = curl_v + u_ek[e,:]*dual_ne[v,iedge]
+        curl_vk[v,:] = curl_v / Av[v]
+    return curl_vk
Index: codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py	(revision 790)
+++ codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py	(revision 791)
@@ -26,4 +26,6 @@
 from dynamico import precision as prec
 from dynamico.meshes import Cartesian_mesh as Mesh
+from dynamico.kernels import grad, curl, div, KE
+from dynamico.LAM import Davies
 
 import math as math
@@ -140,79 +142,4 @@
     gas = thermo.set_vs(v,s)
     return gas, w, z
-
-def grad(mesh, s_ik):
-    llm=s_ik.shape[1]
-    edge_num, left, right = mesh.edge_num, mesh.left, mesh.right
-    grad_ek = prec.zeros((edge_num, llm))
-    for e in range(edge_num):
-        left_e = left[e]-1 # Fortran => Python indexing
-        right_e = right[e]-1 # Fortran => Python indexing
-        grad_ek[e,:] = s_ik[right_e,:]-s_ik[left_e,:]
-    return grad_ek
-
-def div(mesh, u_ek):
-    llm=u_ek.shape[1]
-    primal_num, primal_deg, primal_edge = mesh.primal_num, mesh.primal_deg, mesh.primal_edge
-    primal_ne, le_de, Ai = mesh.primal_ne, mesh.le_de, mesh.Ai
-    div_ik = prec.zeros((primal_num, llm))
-    for i in range(primal_num):
-        div_i = 0.
-        for iedge in range(primal_deg[i]):
-            e = primal_edge[i,iedge]-1 # Fortran => Python indexing
-            div_i = div_i + u_ek[e,:]*primal_ne[i,iedge]*le_de[e]
-        div_ik[i,:] = div_i / Ai[i]
-    return div_ik
-
-def KE(mesh, u_ek):
-    llm=u_ek.shape[1]
-    primal_num, primal_deg, primal_edge = mesh.primal_num, mesh.primal_deg, mesh.primal_edge
-    primal_ne, le_de, Ai = mesh.primal_ne, mesh.le_de, mesh.Ai
-    div_ik = prec.zeros((primal_num, llm))
-    for i in range(primal_num):
-        div_i = 0.
-        for iedge in range(primal_deg[i]):
-            e = primal_edge[i,iedge]-1 # Fortran => Python indexing
-            div_i = div_i + le_de[e]*(u_ek[e,:]**2)
-        div_ik[i,:] = div_i * (.25/ Ai[i])
-    return div_ik
-
-def curl(mesh, u_ek):
-    llm=u_ek.shape[1]
-    dual_num, dual_deg, dual_edge, dual_ne, Av = mesh.dual_num, mesh.dual_deg, mesh.dual_edge, mesh.dual_ne, mesh.Av
-    curl_vk = prec.zeros((dual_num, llm))
-    for v in range(dual_num):
-        curl_v = 0.
-        for iedge in range(dual_deg[v]):
-            e = dual_edge[v,iedge]-1 # Fortran => Python indexing
-            curl_v = curl_v + u_ek[e,:]*dual_ne[v,iedge]
-        curl_vk[v,:] = curl_v / Av[v]
-    return curl_vk
-
-class Davies:
-    def __init__(self,N1,N2,x_i,y_i,x_e,y_e):
-        self.N1, self.N2 = N1, N2
-        self.beta_i = self.mask(x_i,y_i)
-        self.beta_e = self.mask(x_e,y_e)
-        self.x_e,self.y_e = x_e,y_e
-    def mask0(self,c,c0,delta): # 1D building block for Davies relaxation
-        N1, N2 = self.N1, self.N2
-        m = np.zeros(c.size)
-        c1,c2 = c0-N1*delta, c0-(N1+N2)*delta
-        for i in range(c.size):
-            ci=c[i]
-            m[i] = (1.+cos((ci-c2)*pi/(N2*delta)))/2.0
-            if ci < c2 : m[i]=1.
-            if ci > c1 : m[i]=0.
-        return m
-    def relax(self, llm, step, flow):
-        beta_i, beta_e = self.beta_i, self.beta_e
-        m,S,u,Phi,W=flow
-        for l in range(llm):
-            step.mass[:,l]        = beta_i*step.mass[:,l]        + (1.-beta_i)*m[:,l]
-            step.theta_rhodz[:,l] = beta_i*step.theta_rhodz[:,l] + (1.-beta_i)*S[:,l]
-            step.u[:,l]           = beta_e*step.u[:,l]           + (1.-beta_e)*u[:,l]
-        for l in range(llm+1):
-            step.geopot[:,l]      = beta_i*step.geopot[:,l]      + (1.-beta_i)*Phi[:,l]
-            step.W[:,l]           = beta_i*step.W[:,l]           + (1.-beta_i)*W[:,l]
 
 class myDavies(Davies):
