Index: codes/icosagcm/devel/Python/dynamico/meshes.py
===================================================================
--- codes/icosagcm/devel/Python/dynamico/meshes.py	(revision 746)
+++ codes/icosagcm/devel/Python/dynamico/meshes.py	(revision 747)
@@ -11,9 +11,8 @@
 import parallel
 from util import list_stencil, Base_class, inverse_list
+from precision import zeros
 
 radian=180/math.pi
-
-def zeros(dims): return np.zeros([n for n in dims if n>1], dtype=unst.np_num)
-
+  
 #------------------- Hybrid mass-based coordinate -------------
 
@@ -170,5 +169,6 @@
                     (indexu(x,y+1),.25),
                     (indexu(x-1,y+1),.25)))
-        Aiv=np.zeros(nx*ny)+dx*dy # Ai=Av=dx*dy
+        Aiv=np.zeros(nx*ny)
+        Aiv[:]=dx*dy # Ai=Av=dx*dy
         unst.init_mesh(llm,nqdyn,2*nx*ny,nx*ny,nx*ny,4,4,4,
                   primal_nb,primal_edge,primal_ne,
@@ -181,12 +181,11 @@
     def field_z(self,n=1): return zeros((n,self.ny,self.nx,self.llm))
     def field_w(self,n=1): return zeros((n,self.ny,self.nx,self.llm+1))
-    def field_u(self,n=1): 
-        if n==1:
-            return np.zeros((self.ny,2*self.nx,self.llm))
-        else:
-            return np.zeros((n,self.ny,2*self.nx,self.llm))
+    def field_u(self,n=1): return zeros((self.ny,2*self.nx,self.llm))
     def field_ps(self,n=1): return zeros((n,self.ny,self.nx))
-    def ucomp(self,u): return u[:,range(0,2*self.nx,2),:]
-    def set_ucomp(self,uv,u): uv[:,range(0,2*self.nx,2),:]=u
+    def ucomp(self,u): 
+        return u[range(0,2*self.nx,2),:] if self.ny==1 else u[:,range(0,2*self.nx,2),:]
+    def set_ucomp(self,uv,u):
+        if self.ny==1 : uv[range(0,2*self.nx,2),:]=u
+        else : uv[:,range(0,2*self.nx,2),:]=u
     def vcomp(self,u): return u[:,range(1,2*self.nx,2),:]
 
@@ -288,5 +287,5 @@
         return self.de*(ulon*np.cos(self.angle_e)+ulat*np.sin(self.angle_e))
     def ucov3D(self, ulon, ulat):
-        ucov = np.zeros((self.edge_num,ulon.shape[1]))
+        ucov = zeros((self.edge_num,ulon.shape[1]))
         for edge in range(self.edge_num):
             angle=self.angle_e[edge]
Index: codes/icosagcm/devel/Python/dynamico/precision.py
===================================================================
--- codes/icosagcm/devel/Python/dynamico/precision.py	(revision 747)
+++ codes/icosagcm/devel/Python/dynamico/precision.py	(revision 747)
@@ -0,0 +1,10 @@
+"Helper function to manage mixed precision"
+import numpy as np
+from unstructured import np_num
+
+zero=np_num(0.)
+
+# cast array to float/double depending on CPP_MIXED_PREC
+def asnum(data): return [ x.astype(np_num) for x in data] if type(data) is list else data.astype(np_num)
+
+def zeros(dims): return np.zeros([n for n in dims if n>1], dtype=np_num)
Index: codes/icosagcm/devel/Python/dynamico/time_step.py
===================================================================
--- codes/icosagcm/devel/Python/dynamico/time_step.py	(revision 746)
+++ codes/icosagcm/devel/Python/dynamico/time_step.py	(revision 747)
@@ -85,10 +85,18 @@
     def cjl(self,b): # converts Butcher tableau b(j,l) to update coefficients c(j,l)
         return np.array([ b[j+1,:]-b[j,:] for j in range(self.nstage)])
-    def __init__(self,bwd_fast_slow, dt,nstage, bsjl,bfjl):
+    def __init__(self,precision,bwd_fast_slow, dt,nstage, bsjl,bfjl):
+        if bwd_fast_slow is None:
+            print 'SIRK with Fortran-side time stepping'
+            precision=np.float64
+        else:
+            print 'SIRK with Python-side time stepping at precision ', precision
+
         self.bwd_fast_slow = bwd_fast_slow
         self.dt, self.nstage = dt, nstage
-        self.csjl = dt*self.cjl(bsjl)
-        self.cfjl = dt*self.cjl(bfjl)
-        self.tauj = np.array([dt*bfjl[j,j] for j in range(self.nstage)])
+        csjl = dt*self.cjl(bsjl)
+        cfjl = dt*self.cjl(bfjl)
+        self.csjl, self.cfjl = [x.astype(precision) for x in csjl, cfjl]
+        self.tauj = np.array([dt*bfjl[j,j] for j in range(self.nstage)], dtype=precision)
+        print 'Types of csjl, cfjl, tauj :', self.csjl.dtype, self.cfjl.dtype, self.tauj.dtype
     def next(self,zj):
         csjl, cfjl, tauj = self.csjl, self.cfjl, self.tauj
@@ -100,4 +108,5 @@
             for l in range(0,j+1): # here zj is initially y(j) and finally z(j+1)
                 zj = [z + csjl[j,l]*s + cfjl[j,l]*f for z,s,f in zip(zj, slowj[l], fastj[l]) ]
+#            print 'SIRK : zj has type ', [ z.dtype for z in zj]
         return zj        
     def advance(self, flow, N):
@@ -113,5 +122,5 @@
     return SIRK(bwd_fast_slow,dt,1,bsjl,bfjl)
         
-def ARK2(bwd_fast_slow, dt, a32 = (3.+2.*np.sqrt(2.))/6.):
+def ARK2(bwd_fast_slow, dt, precision=None, a32 = (3.+2.*np.sqrt(2.))/6.):
     """ ARK2 scheme by Giraldo et al. (2013), noted ARK2(2,3,2) in Weller et al. (2013) """ 
     gamma = 1.-1./np.sqrt(2.)
@@ -120,5 +129,5 @@
     bsjl = np.array([[0,0,0],[2.*gamma,0,0],[1.-a32,a32,0],wj])
     bfjl = np.array([[0,0,0],[gamma,gamma,0],[delta,delta,gamma],wj])
-    return SIRK(bwd_fast_slow,dt,3,bsjl,bfjl)
+    return SIRK(precision,bwd_fast_slow,dt,3,bsjl,bfjl)
         
 def UJ3(bwd_fast_slow, dt):
@@ -135,5 +144,5 @@
 ######### Explicit RK schemes written as ARK with two identical Butcher tableaus ########
 
-def RK4(bwd_fast_slow, dt):
+def RK4(bwd_fast_slow, dt, precision=None):
     """ 4-stage 4th-order RK """ 
     bjl = np.array([[0,0,0,0],
@@ -142,5 +151,5 @@
                     [0.,0.,1.,0],
                     [1/6.,1/3.,1/3.,1/6.]])
-    return SIRK(bwd_fast_slow,dt,4,bjl,bjl)
+    return SIRK(precision, bwd_fast_slow,dt,4,bjl,bjl)
 
 ######### Simplified explicit RK schemes (2nd-order) ########
Index: codes/icosagcm/devel/Python/dynamico/wrap.py
===================================================================
--- codes/icosagcm/devel/Python/dynamico/wrap.py	(revision 746)
+++ codes/icosagcm/devel/Python/dynamico/wrap.py	(revision 747)
@@ -11,5 +11,5 @@
 # defines converter functions depending on type of Python argument
 def noop(x): return x
-py2c = { np.ndarray:loc, np.float64:c_double, float:c_double, int:c_int, 
+py2c = { np.ndarray:loc, np.float64:c_double, np.float32:c_float, int:c_int, 
          c_void_p:noop, str:noop, c_int:noop }
 
Index: codes/icosagcm/devel/Python/src/unstructured.pyx
===================================================================
--- codes/icosagcm/devel/Python/src/unstructured.pyx	(revision 746)
+++ codes/icosagcm/devel/Python/src/unstructured.pyx	(revision 747)
@@ -239,4 +239,5 @@
         # h*s = h => uniform buoyancy s=1 => shallow-water
         dh, du_slow, du_fast, hs, buf = self.fmass(), self.fu(), self.fu(), h.copy(), self.geopot
+        assert type(tau) is np_num, 'tau must be of type unstructured.np_num.'
         ker.dynamico_caldyn_unstructured(tau, self.ms, h, hs, u, self.geopot, buf,
                   self.s, self.ps, self.pk, self.hflux, self.qv,
@@ -296,4 +297,9 @@
               left,right,down,up,trisk_deg,trisk,
               Ai, Av, fv, le_de, Riv2, wee):
+
+    print 'Types of Ai, Av, fv, le_de, Riv2, wee : ', Ai.dtype,Av.dtype,fv.dtype,le_de.dtype,Riv2.dtype,wee.dtype
+    for var,varname in zip((Ai,Av,fv,le_de,Riv2,wee), ('Ai','Av','fv','le_de','Riv2','wee')):
+        assert var.dtype == np.float64, '%s must be double precision'%varname
+
     setvars( ('llm','nqdyn','edge_num','primal_num','dual_num',
               'max_trisk_deg','max_primal_deg','max_dual_deg'),
Index: codes/icosagcm/devel/Python/test/py/DCMIP2008c5.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/DCMIP2008c5.py	(revision 746)
+++ codes/icosagcm/devel/Python/test/py/DCMIP2008c5.py	(revision 747)
@@ -136,8 +136,7 @@
 plt.savefig('fig_DCMIP2008c5/Ai.png'); plt.close()
 
-
-if False: # time stepping in Python
+if False: # time stepping in Python => set precision to np_num
     caldyn = unst.Caldyn_HPE(caldyn_thermo,caldyn_eta, mesh,thermo,params,params.g)
-    scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, a32=0.7)
+    scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, precision=unst.np_num, a32=0.7)
     def next_flow(m,S,u): 
         m,S,u = scheme.advance((m,S,u),nt)
@@ -179,4 +178,7 @@
 z, vol = mesh.field_mass(), mesh.field_mass()
 m,S,u,Phi,W = flow0 
+
+print 'types of m,S,u :', m.dtype, S.dtype, u.dtype
+
 caldyn_step.geopot[:,0]=Phi[:,0]
 plots(0)
Index: codes/icosagcm/devel/Python/test/py/NH_3D_bubble.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/NH_3D_bubble.py	(revision 746)
+++ codes/icosagcm/devel/Python/test/py/NH_3D_bubble.py	(revision 747)
@@ -4,4 +4,5 @@
 from dynamico import DCMIP
 from dynamico import meshes
+from dynamico import precision as prec
 from dynamico.meshes import Cartesian_mesh as Mesh
 import math as math
@@ -51,5 +52,5 @@
     params.rho_bot = 1e6/gas_bot.v
     
-    return thermo, mesh, params, (mass_ik,Sik,ujk,Phi_il,Wil), gas
+    return thermo, mesh, params, prec.asnum([mass_ik,Sik,ujk,Phi_il,Wil]), gas
 
 #Lx, nx, llm, thetac, T, Nslice, courant = 2000., 100, 50, 30., 5., 10, 2.8
@@ -63,4 +64,5 @@
 # compute hybrid coefs from initial distribution of mass
 mass_bl,mass_dak,mass_dbk = meshes.compute_hybrid_coefs(flow0[0])
+print 'Type of mass_bl, mass_dak, mass_dbk : ', [x.dtype for x in mass_bl, mass_dak, mass_dbk]
 unst.ker.dynamico_init_hybrid(mass_bl,mass_dak,mass_dbk)
 
Index: codes/icosagcm/devel/Python/test/py/RSW_2D.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/RSW_2D.py	(revision 746)
+++ codes/icosagcm/devel/Python/test/py/RSW_2D.py	(revision 747)
@@ -18,7 +18,7 @@
 
 x1,x2,yy = mesh.xx[:,:,0]-1., mesh.xx[:,:,0]+1., mesh.yy[:,:,0]
-h0 = 1+0.1*(np.exp(-2.*(x1*x1+yy*yy))+
+h0, u0 = mesh.field_mass(), mesh.field_u()  # flow initally at rest
+h0[:] = 1+0.1*(np.exp(-2.*(x1*x1+yy*yy))+
             np.exp(-2.*(x2*x2+yy*yy)))
-u0 = mesh.field_u()  # flow initally at rest
 flow0=(h0,u0)
 
@@ -32,7 +32,10 @@
 scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt)
 
+print 'types : yy, u0, h0, u0', yy.dtype, u0.dtype, h0.dtype, u0.dtype
+
 flow=flow0
 for i in range(10):
     h,u=flow
+    print 'types : h,u', h.dtype, u.dtype
     caldyn.bwd_fast_slow(flow, 0.)
     plt.figure(); plt.pcolor(mesh.x,mesh.y,caldyn.qv)
@@ -48,4 +51,6 @@
         unst.elapsed=0.
         flow = scheme.advance(flow,N)
+        h,u=flow
+        print 'types : h,u', h.dtype, u.dtype
         # flops
         # mass flux : 1add+1mul per edge          =>  4
Index: codes/icosagcm/devel/Python/test/py/RSW_2D_mesh.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/RSW_2D_mesh.py	(revision 746)
+++ codes/icosagcm/devel/Python/test/py/RSW_2D_mesh.py	(revision 747)
@@ -45,5 +45,5 @@
 dx,dy=Lx/nx,Ly/ny
 
-filename, llm, nqdyn, g, f, radius = 'cart128.nc', 1, 1, 1., 1., 1.
+filename, llm, nqdyn, g, f, radius = 'cart128.nc', 1, 1, 1., 1., None
 unst.setvar('g',g)
 
@@ -64,8 +64,8 @@
 
 x1,x2,yy = xx-1., xx+1., yy
+u0 = mesh.field_u()  # flow initally at rest
 h0 = 1+0.1*(np.exp(-2.*(x1*x1+yy*yy))+
             np.exp(-2.*(x2*x2+yy*yy)))
-u0 = mesh.field_u()  # flow initally at rest
-flow0=(h0,u0)
+flow0=meshes.asnum([h0,u0])
 
 #print('h0 : ', h0[1234]-1.)
@@ -79,5 +79,6 @@
 dt=T/N
 print N,dt,Lx/nx
-scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt)
+#scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, precision=unst.np_num)
+scheme = time_step.RK4(caldyn.bwd_fast_slow, dt, precision=unst.np_num)
 
 flow=flow0
@@ -90,5 +91,5 @@
 for i in range(10):
     h,u=flow
-    flow, fast, slow = caldyn.bwd_fast_slow(flow, 0.)
+    flow, fast, slow = caldyn.bwd_fast_slow(flow, unst.zero_num)
     junk, du_fast = fast
     dh, du_slow = slow
@@ -102,5 +103,6 @@
     plt.figure(); plt.pcolor(x,y,reshape(caldyn.qv))
     plt.colorbar(); plt.title('potential vorticity')
-    plt.savefig('fig_RSW_2D_mesh/%02d.png'%i)
+    plt.savefig('fig_RSW_2D_mesh/%03d.png'%i)
+    plt.close()
 #    plt.figure(); plt.pcolor(mesh.x,mesh.y,h)
 #    plt.colorbar(); plt.title('h')
@@ -112,4 +114,5 @@
         unst.elapsed=0.
         flow = scheme.advance(flow,N)
+#        flow = scheme.advance(flow,5)
         # flops
         # mass flux : 1add+1mul per edge          =>  4
Index: codes/icosagcm/devel/Python/test/py/RSW_MPAS_W02.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/RSW_MPAS_W02.py	(revision 746)
+++ codes/icosagcm/devel/Python/test/py/RSW_MPAS_W02.py	(revision 747)
@@ -1,4 +1,5 @@
 print 'Loading DYNAMICO modules ...'
 from dynamico import unstructured as unst
+from dynamico.precision import asnum
 from dynamico.meshes import MPAS_Format, Unstructured_Mesh
 from dynamico import time_step
@@ -34,5 +35,5 @@
 dt = T/nt
 #scheme = time_step.RKn_simple(1,caldyn.bwd_fast_slow, dt)
-scheme = time_step.RK4(caldyn.bwd_fast_slow, dt)
+scheme = time_step.RK4(caldyn.bwd_fast_slow, dt, precision=unst.np_num)
     
 print dx, dt, dt*c0/dx
@@ -61,5 +62,7 @@
 fu_perp = mesh.ucov2D(0.*ulon,2*Omega*np.sin(mesh.lat_e)*ulon)
 
-flow = (Phi0,mesh.ucov2D(ulon,0.*ulon))
+flow = asnum([Phi0, mesh.ucov2D(ulon,0.*ulon)])
+print 'type of Phi0,u0 : ', flow[0].dtype, flow[1].dtype
+
 for i in range(N):
     h,u=flow
Index: codes/icosagcm/devel/Python/test/py/bubble.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/bubble.py	(revision 746)
+++ codes/icosagcm/devel/Python/test/py/bubble.py	(revision 747)
@@ -1,7 +1,2 @@
-import math as math
-import numpy as np
-import matplotlib.pyplot as plt
-import matplotlib.animation as manimation
-
 from dynamico.meshes import Cartesian_mesh as Mesh
 import dynamico.dyn as dyn
@@ -10,4 +5,10 @@
 import dynamico.advect as advect
 from dynamico import unstructured as unst
+
+import math as math
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.animation as manimation
+
 
 def reldiff(a,b):
Index: codes/icosagcm/devel/Python/test/py/slice_GW_NH.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/slice_GW_NH.py	(revision 746)
+++ codes/icosagcm/devel/Python/test/py/slice_GW_NH.py	(revision 747)
@@ -29,5 +29,5 @@
 
     u0=mesh.field_u() 
-    u0[:,range(0,2*nx,2),:] = u*mesh.dx # Doppler shift by u
+    u0[range(0,2*nx,2),:] = u*mesh.dx # Doppler shift by u
     flow0=(m0ik,S0ik,u0,Phi0_il,0*Phi0_il)
 
@@ -81,5 +81,5 @@
         plt.figure(figsize=(12,3))
         ucomp = mesh.ucomp(u)
-        plt.pcolor(xx,zz,ucomp[0,:,:]/dx)
+        plt.pcolor(xx,zz,ucomp[:,:]/dx)
         plt_axes()
         plt.title("u(t = %f s)" % ((i+1)*dt*nt))
Index: codes/icosagcm/devel/Python/test/py/slice_GW_hydro.py
===================================================================
--- codes/icosagcm/devel/Python/test/py/slice_GW_hydro.py	(revision 746)
+++ codes/icosagcm/devel/Python/test/py/slice_GW_hydro.py	(revision 747)
@@ -23,5 +23,5 @@
 
     u0=mesh.field_u()
-    u0[:,range(0,2*nx,2),:] = 100.*mesh.dx # Doppler shift 100 m/s
+    u0[range(0,2*nx,2),:] = 100.*mesh.dx # Doppler shift 100 m/s
     flow0=(m0ik,S0ik,u0)
 
@@ -60,5 +60,5 @@
         plt.figure(figsize=(12,3))
         ucomp = mesh.ucomp(u)
-        plt.pcolor(xx,zz,ucomp[0,:,:]/dx)
+        plt.pcolor(xx,zz,ucomp[:,:]/dx)
         plt.ylim(0,10.)
         plt.colorbar()
Index: codes/icosagcm/devel/make_icosa
===================================================================
--- codes/icosagcm/devel/make_icosa	(revision 746)
+++ codes/icosagcm/devel/make_icosa	(revision 747)
@@ -227,2 +227,5 @@
   ./build --job $job
 fi
+
+# force recompilation of cython modules
+touch Python/src/*.pyx
Index: codes/icosagcm/devel/make_python
===================================================================
--- codes/icosagcm/devel/make_python	(revision 746)
+++ codes/icosagcm/devel/make_python	(revision 747)
@@ -32,8 +32,10 @@
 	update "$x" "$DYNAMICO_ROOT/src/$2/$x"
     done
+
 }
 
 function cmd_kernels() # this function is invoked by : './make_python kernels'
 {
+    cmd_clean
     cd $KERNELS
     ./codegen hexagonal hex_master unstructured
Index: codes/icosagcm/devel/src/unstructured/data_unstructured.F90
===================================================================
--- codes/icosagcm/devel/src/unstructured/data_unstructured.F90	(revision 746)
+++ codes/icosagcm/devel/src/unstructured/data_unstructured.F90	(revision 747)
@@ -205,12 +205,12 @@
     Riv2(:,:)=Riv2_(:,:)
     wee(:,:) = wee_(:,:)
-    PRINT *, MAXVAL(ABS(Ai))
-    PRINT *, MAXVAL(ABS(Av))
-    PRINT *, MAXVAL(ABS(fv))
-    PRINT *, MAXVAL(ABS(le_de))
-    PRINT *, MAXVAL(ABS(Riv2))
-    PRINT *, MAXVAL(ABS(wee))
-    PRINT *, MINVAL(right), MAXVAL(right)
-    PRINT *, MINVAL(right), MAXVAL(left)
+    PRINT *, 'Max Ai : ',    MAXVAL(ABS(Ai))
+    PRINT *, 'Max Av : ',    MAXVAL(ABS(Av))
+    PRINT *, 'Max fv : ',    MAXVAL(ABS(fv))
+    PRINT *, 'Max le_de : ', MAXVAL(ABS(le_de))
+    PRINT *, 'Max Riv2 : ',  MAXVAL(ABS(Riv2))
+    PRINT *, 'Max wee : ',   MAXVAL(ABS(wee))
+    PRINT *, MINVAL(right),  MAXVAL(right)
+    PRINT *, MINVAL(right),  MAXVAL(left)
     PRINT *,' ... Done.'
     IF(nb_threads==0) nb_threads=OMP_GET_MAX_THREADS()
