Changeset 758

Timestamp:

10/10/18 18:12:01 (6 years ago)

Author:

dubos

Message:

devel/Python : added ncwrite to Cartesian_mesh + bugfixes + cleanup

Location:

codes/icosagcm/devel/Python

Files:

• dynamico/meshes.py (modified) (7 diffs)
• test/py/write_Cartesian_mesh.py (modified) (2 diffs)

codes/icosagcm/devel/Python/dynamico/meshes.py

@@ -81 +81 @@
     deg[ij]=k
 
-class Cartesian_mesh:
+def concat(x,y):
+    z = np.asarray([x,y])
+    return z.transpose()
+
+class Cartesian_mesh(Base_class):
     def __init__(self,nx,ny,llm,nqdyn,Lx,Ly,f):
         dx,dy = Lx/nx, Ly/ny
@@ -103 +107 @@
         indexu=lambda x,y : 2*index(x,y)-1
         indexv=lambda x,y : 2*index(x,y)
-        indices = lambda shape : np.zeros(shape,np.int32)
-
-        primal_nb = indices(nx*ny)
-        primal_edge = indices((nx*ny,4))
-        primal_ne = indices((nx*ny,4))
+        def zeros(shape,n=1): return [np.zeros(shape) for i in range(n)] if n>1 else np.zeros(shape)
+        def indices(shape,n=1): return [np.zeros(shape,np.int32) for i in range(n)] if n>1 else np.zeros(shape,np.int32)
+
+        Aiv, lon_i,lon_v,lat_i,lat_v   = zeros( nx*ny,   5)
+        lon_e,lat_e,le,de,angle_e      = zeros( 2*nx*ny, 5)
+        Riv2                           = zeros((nx*ny,4)  )
+        wee                            = zeros((2*nx*ny,4))
+
+        primal_deg,dual_deg                 = indices( nx*ny,   2)
+        primal_edge,primal_ne,primal_vertex = indices((nx*ny,4),3)
+        dual_edge,dual_ne,dual_vertex       = indices((nx*ny,4),3)
+        trisk_deg,left,right,up,down        = indices( 2*nx*ny, 5)
+        trisk                               = indices((2*nx*ny,4)) # 4 TRiSK coefs per edge
     
-        dual_nb = indices(nx*ny)
-        dual_edge = indices((nx*ny,4))
-        dual_ne = indices((nx*ny,4))
-        dual_vertex = indices((nx*ny,4))
-        Riv2 = np.zeros((nx*ny,4))
-    
-        left = indices(2*nx*ny)
-        right = indices(2*nx*ny)
-        up = indices(2*nx*ny)
-        down = indices(2*nx*ny)
-        le_de = np.zeros(2*nx*ny)
-    
-        trisk_deg = indices(2*nx*ny)
-        trisk = indices((2*nx*ny,4)) # 4 TRiSK coefs per edge
-        wee = np.zeros((2*nx*ny,4))
-    
+        def putval(ij, arrays, vals): 
+            for array,v in zip(arrays,vals): array[ij]=v
         for x in range(nx):
             for y in range(ny):
                 # NB : Fortran indices start at 1
                 # primal cells
-                put(index(x,y)-1,primal_nb,(primal_edge,primal_ne),
-                    ((indexu(x,y),1),
-                     (indexv(x,y),1),
-                     (indexu(x-1,y),-1),
-                     (indexv(x,y-1),-1) ))
+                ij=index(x,y)-1
+                lon_i[ij], lat_i[ij] = x, y
+                put(index(x,y)-1, primal_deg,(primal_edge,primal_vertex,primal_ne),
+                    ((indexu(x,y),index(x,y),1), 
+                    (indexv(x,y),index(x-1,y),1),
+                    (indexu(x-1,y),index(x-1,y-1),-1),
+                    (indexv(x,y-1),index(x,y-1),-1) ))
                 # dual cells
-                put(index(x,y)-1,dual_nb,(dual_edge,dual_vertex,dual_ne,Riv2),
+                ij=index(x,y)-1
+                lon_v[ij], lat_v[ij] = x+.5, y+.5
+                put(ij,dual_deg,(dual_edge,dual_vertex,dual_ne,Riv2),
                    ((indexv(x+1,y),index(x,y),1,.25),
                     (indexu(x,y+1),index(x+1,y),-1,.25),
@@ -147 +149 @@
                 # u-points
                 ij =indexu(x,y)-1 
-                left[ij]=index(x,y)
-                right[ij]=index(x+1,y)
-                down[ij]=index(x,y-1)
-                up[ij]=index(x,y)
-                le_de[ij]=dy/dx
+                putval(ij, (left, right, down, up, le, de, angle_e, lon_e, lat_e),
+                       (index(x,y), index(x+1,y), index(x,y-1), index(x,y), dy, dx, 0., x+.5, y ) )
                 put(ij,trisk_deg,(trisk,wee),(
-                    (indexv(x,y),-.25),
-                    (indexv(x+1,y),-.25),
-                    (indexv(x,y-1),-.25),
-                    (indexv(x+1,y-1),-.25)))
+                    (indexv(x,y),.25),
+                    (indexv(x+1,y),.25),
+                    (indexv(x,y-1),.25),
+                    (indexv(x+1,y-1),.25)))
                 # v-points
                 ij = indexv(x,y)-1
-                left[ij]=index(x,y)
-                right[ij]=index(x,y+1)
-                down[ij]=index(x,y)
-                up[ij]=index(x-1,y)
-                le_de[ij]=dx/dy
+                putval(ij, (left, right, down, up, le, de, angle_e, lon_e, lat_e),
+                       (index(x,y), index(x,y+1), index(x,y), index(x-1,y), dx, dy, 5.*math.pi, x, y+.5 ) )
                 put(ij,trisk_deg,(trisk,wee),(
-                    (indexu(x,y),.25),
-                    (indexu(x-1,y),.25),
-                    (indexu(x,y+1),.25),
-                    (indexu(x-1,y+1),.25)))
-        Aiv=np.zeros(nx*ny)
+                    (indexu(x,y),-.25),
+                    (indexu(x-1,y),-.25),
+                    (indexu(x,y+1),-.25),
+                    (indexu(x-1,y+1),-.25)))
+
+        primal_i, primal_j = [ x.astype(np.int32) for x in lon_i, lat_i]
+        lon_i, lon_v, lon_e = [x*dx-Lx/2 for x in lon_i, lon_v, lon_e]
+        lat_i, lat_v, lat_e = [y*dy-Ly/2 for y in lat_i, lat_v, lat_e]
+
         Aiv[:]=dx*dy # Ai=Av=dx*dy
+        le_de=le/de
         unst.init_mesh(llm,nqdyn,2*nx*ny,nx*ny,nx*ny,4,4,4,
-                  primal_nb,primal_edge,primal_ne,
-                  dual_nb,dual_edge,dual_ne,dual_vertex,
+                  primal_deg,primal_edge,primal_ne,
+                  dual_deg,dual_edge,dual_ne,dual_vertex,
                   left,right,down,up,trisk_deg,trisk,
                   Aiv,Aiv,f+0.*Aiv,le_de,Riv2,-wee)
+        self.set_members(locals(), 'primal_deg', 'primal_edge', 'primal_ne', 'primal_vertex',
+                         'dual_deg', 'dual_edge', 'dual_ne', 'dual_vertex',
+                         'left','right','down','up',
+                         'trisk_deg','trisk','Riv2','wee',
+                         'Aiv','lon_i','lat_i','lon_v','lat_v',
+                         'le_de','le','de','lon_e','lat_e','angle_e')
+    def ncwrite(self, name):
+        """The method writes Cartesian mesh on the disc.
+            Args: Mesh parameters"""
+
+        #----------------OPENING NETCDF OUTPUT FILE------------
+
+        try:
+            f = cdf.Dataset(name, 'w', format='NETCDF4')
+        except:
+            print("Cartesian_mesh.ncwrite : Error occurred while opening new netCDF mesh file")
+            raise
+
+        #----------------DEFINING DIMENSIONS--------------------
+
+        for dimname, dimsize in [("primal_cell", self.primal_deg.size),
+                                 ("dual_cell", self.dual_deg.size),
+                                 ("edge", self.left.size),
+                                 ("primal_edge_or_vertex", self.primal_edge.shape[1]),
+                                 ("dual_edge_or_vertex", self.dual_edge.shape[1]),
+                                 ("TWO", 2),
+                                 ("trisk_edge", 4)]:
+            f.createDimension(dimname,dimsize)
+
+        #----------------DEFINING VARIABLES---------------------
+        
+        f.description = "Cartesian_mesh"
+        f.nx, f.ny = self.nx, self.ny
+
+        def create_vars(dimname, info):
+            for vname, vtype, vdata in info:
+                var = f.createVariable(vname,vtype,dimname)
+                var[:] = vdata
+                
+        create_vars("primal_cell", 
+                    [("primal_deg","i4",self.primal_deg), 
+                     ("Ai","f8",self.Aiv),
+                     ("lon_i","f8",self.lon_i),
+                     ("lat_i","f8",self.lat_i)] )
+        create_vars("dual_cell", 
+                    [("dual_deg","i4",self.dual_deg), 
+                     ("Av","f8",self.Aiv),
+                     ("lon_v","f8",self.lon_v),
+                     ("lat_v","f8",self.lat_v),
+                     ] )
+
+        create_vars( ("primal_cell","primal_edge_or_vertex"), 
+                     [("primal_edge", "i4", self.primal_edge),
+                      ("primal_ne", "i4", self.primal_ne),
+                      ("primal_vertex", "i4", self.primal_vertex)] )
+
+        create_vars( ("dual_cell","dual_edge_or_vertex"), 
+                     [("dual_edge", "i4", self.dual_edge),
+                      ("dual_vertex","i4",self.dual_vertex),
+                      ("dual_ne","i4",self.dual_ne),
+                      ("Riv2","f8",self.Riv2)] )
+
+        create_vars("edge",
+                    [("trisk_deg","i4",self.trisk_deg),
+                     ("le","f8",self.le),
+                     ("de","f8",self.de),
+                     ("lon_e","f8", self.lon_e),
+                     ("lat_e","f8", self.lat_e),
+                      ("angle_e","f8", self.angle_e)] )
+
+        create_vars( ("edge","TWO"),
+                     [("edge_left_right","i4", concat(self.left,self.right)),
+                      ("edge_down_up","i4", concat(self.down,self.up))] )
+                      
+
+        create_vars( ("edge","trisk_edge"),
+                    [("trisk","i4",self.trisk),
+                     ("wee","f8",self.wee)] )
+
+        f.close()
 
     def field_theta(self,n=1): return zeros((n,self.nqdyn,self.ny,self.nx,self.llm))
@@ -371 +452 @@
         
 class Unstructured_PMesh(Abstract_Mesh): # Mesh data distributed across MPI processes
-    def __init__(self,comm, gridfile):
+    def __init__(self,comm, gridfile, meshtype='unstructured'):
         self.comm, self.gridfile, get_pvars = comm, gridfile, gridfile.get_pvars
         dim_primal, dim_edge, dim_dual = gridfile.get_pdims(comm,
-            'primal_num','edge_num','dual_num')
+            'primal_cell','edge','dual_cell')
         primal_deg, primal_vertex, primal_edge, lon_i, lat_i, Ai = get_pvars(
             dim_primal, 'primal_deg', 'primal_vertex', 'primal_edge', 'lon_i', 'lat_i', 'Ai' )
@@ -397 +478 @@
         primal_owner = partition_from_stencil(edge_owner, primal_deg, primal_edge)
         primal_owner = parallel.LocPArray1D(dim_primal, primal_owner)
-        self.set_members(locals(), 'dim_primal', 'dim_dual', 'dim_edge',
+        self.set_members(locals(), 'meshtype', 'dim_primal', 'dim_dual', 'dim_edge',
                          'edge_owner', 'primal_owner', 'primal_deg', 'primal_vertex', 
                          'lon_v', 'lat_v', 'Av', 'lon_i', 'lat_i', 'Ai', 
                          'le', 'de', 'lon_e', 'lat_e', 'angle_e')
+#        if meshtype == 'curvilinear' : # read extra information from mesh file
+#            self.primal_i, self.primal_j = get_pvars(dim_primal, 'primal_i','primal_j')
+            
     def get(self, getter, *names): return [ getter(self.__dict__[x]) for x in names ]
 
@@ -445 +529 @@
             get_all_edges, 'le', 'de', 'lon_e', 'lat_e', 'angle_e')
         primal_num, dual_num, edge_num = map(len, (cells_C1, vertices_V1, edges_E2))
+
+        if pmesh.meshtype == 'curvilinear' :
+            self.primal_i, self.primal_j = pmesh.get(get_own_cells, 'primal_i', 'primal_j')
 
         # construct local stencils from global stencils
@@ -655 +742 @@
     # example : edge_list = progressive_list(E0,E1,E2) with E0,E1,E2 increasing lists of cell edges
     return list(progressive_iter([], cell_lists))
+from dynamico.meshes import zeros
+import numpy as np
+import netCDF4 as cdf
+import argparse
+
+
+if False:
+    def __init__(self,nx,ny,llm,nqdyn,Lx,Ly):
+        dx,dy = Lx/nx, Ly/ny
+        self.dx, self.dy = dx,dy
+        self.nx, self.ny, self.llm, self.nqdyn = nx,ny,llm,nqdyn
+        self.field_z = self.field_mass
+        # 1D coordinate arrays
+        x=(np.arange(nx)-nx/2.)*Lx/nx
+        y=(np.arange(ny)-ny/2.)*Ly/ny
+        lev=np.arange(llm)
+        levp1=np.arange(llm+1)
+        self.x, self.y, self.lev, self.levp1 = x,y,lev,levp1
+        # 3D coordinate arrays
+        self.xx,self.yy,self.ll = np.meshgrid(x,y,lev, indexing='ij')
+        self.xxp1,self.yyp1,self.llp1 = np.meshgrid(x,y,levp1, indexing='ij')
+        # beware conventions for indexing
+        # Fortran order : llm,nx*ny,nqdyn  / indices start at 1
+        # Python order : nqdyn,ny,nx,llm   / indices start at 0
+        # indices below follow Fortran while x,y follow Python/C
+        index=lambda x,y : ((x+(nx*(y+2*ny)))%(nx*ny))+1
+        indexu=lambda x,y : 2*index(x,y)-1
+        indexv=lambda x,y : 2*index(x,y)
+        indices = lambda shape : np.zeros(shape,np.int32)
+
+        for x in range(nx):
+            for y in range(ny):
+                # NB : Fortran indices start at 1
+                # primal cells
+                put(index(x,y)-1, primal_deg,(primal_edge,primal_vertex,primal_ne),
+                    ((indexu(x,y),index(x,y),1), 
+                    (indexv(x,y),index(x-1,y),1),
+                    (indexu(x-1,y),index(x-1,y-1),-1),
+                    (indexv(x,y-1),index(x,y-1),-1) ))
+                # dual cells
+                put(index(x,y)-1,dual_deg,(dual_edge,dual_vertex,dual_ne,Riv2),
+                   ((indexv(x+1,y),index(x,y),1,.25),
+                    (indexu(x,y+1),index(x+1,y),-1,.25),
+                    (indexv(x,y),index(x+1,y+1),-1,.25),
+                    (indexu(x,y),index(x,y+1),1,.25) ))
+                # edges :
+                # left and right are adjacent primal cells
+                # flux is positive when going from left to right
+                # up and down are adjacent dual cells (vertices)
+                # circulation is positive when going from down to up
+
+        Aiv=np.zeros(nx*ny)+dx*dy  
+        Ai=Av=np.zeros(nx*ny)+dx*dy
+
+        self.llm=llm
+        self.nqdyn=nqdyn
+        self.nx=nx
+        self.ny=ny
+        self.primal_deg=primal_deg
+        self.primal_edge=primal_edge
+        self.primal_ne=primal_ne
+        self.dual_deg=dual_deg
+        self.dual_edge=dual_edge
+        self.dual_ne=dual_ne
+        self.dual_vertex=dual_vertex
+        self.primal_vertex=primal_vertex
+        self.left=left
+        self.right=right
+        self.down=down
+        self.up=up
+        self.trisk_deg=trisk_deg
+        self.trisk=trisk
+        self.Aiv=Aiv
+        self.Ai=Ai
+        self.Av=Av
+        self.le=le
+        self.de=de
+        self.angle_e=angle_e
+        self.Riv2=Riv2
+        self.wee=wee
+        self.lon_i = lon_i
+        self.lon_v = lon_v
+        self.lon_e = lon_e
+        #self.lon_ev = indices(2*nx*ny)
+        self.lat_i = lat_i
+        self.lat_v = lat_v
+        self.lat_e = lat_e
+        #self.lat_ev = indices(2*nx*ny)
+      
+
+    def field_theta(self,n=1): return zeros((n,self.nqdyn,self.ny,self.nx,self.llm))
+    def field_mass(self,n=1): return zeros((n,self.ny,self.nx,self.llm))
+    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_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 vcomp(self,u): return u[:,range(1,2*self.nx,2),:]
+
+

codes/icosagcm/devel/Python/test/py/write_Cartesian_mesh.py

@@ -1 +1 @@
 from dynamico.meshes import zeros
+from dynamico import meshes
 import numpy as np
 import netCDF4 as cdf
 import argparse
 
-#----------------------- Cartesian mesh -----------------------
-
-def concat(x,y):
-    z = np.asarray([x,y])
-    return z.transpose()
-
-# arrays is a list of arrays
-# vals is a list of tuples
-# each tuple is stored in each array
-def put(ij, deg, arrays, vals):
-    k=0
-    for vv in vals: # vv is a tuple of values to be stored in arrays
-        for array,v in zip(arrays,vv):
-            array[ij,k]=v
-        k=k+1
-    deg[ij]=k
-
-class Cartesian_mesh:
-    def __init__(self,nx,ny,llm,nqdyn,Lx,Ly):
-        dx,dy = Lx/nx, Ly/ny
-        self.dx, self.dy = dx,dy
-        self.nx, self.ny, self.llm, self.nqdyn = nx,ny,llm,nqdyn
-        self.field_z = self.field_mass
-        # 1D coordinate arrays
-        x=(np.arange(nx)-nx/2.)*Lx/nx
-        y=(np.arange(ny)-ny/2.)*Ly/ny
-        lev=np.arange(llm)
-        levp1=np.arange(llm+1)
-        self.x, self.y, self.lev, self.levp1 = x,y,lev,levp1
-        # 3D coordinate arrays
-        self.xx,self.yy,self.ll = np.meshgrid(x,y,lev, indexing='ij')
-        self.xxp1,self.yyp1,self.llp1 = np.meshgrid(x,y,levp1, indexing='ij')
-        # beware conventions for indexing
-        # Fortran order : llm,nx*ny,nqdyn  / indices start at 1
-        # Python order : nqdyn,ny,nx,llm   / indices start at 0
-        # indices below follow Fortran while x,y follow Python/C
-        index=lambda x,y : ((x+(nx*(y+2*ny)))%(nx*ny))+1
-        indexu=lambda x,y : 2*index(x,y)-1
-        indexv=lambda x,y : 2*index(x,y)
-        indices = lambda shape : np.zeros(shape,np.int32)
-
-        primal_deg = indices(nx*ny)
-        primal_edge = indices((nx*ny,4))
-        primal_ne = indices((nx*ny,4))
-
-        dual_deg = indices(nx*ny)
-        dual_edge = indices((nx*ny,4))
-        dual_ne = indices((nx*ny,4))
-        primal_vertex = indices((nx*ny,4))
-        dual_vertex = indices((nx*ny,4))
-        Riv2 = np.zeros((nx*ny,4))
-
-        left = indices(2*nx*ny)
-        right = indices(2*nx*ny)
-        up = indices(2*nx*ny)
-        down = indices(2*nx*ny)
-        le_de = np.zeros(2*nx*ny)
-        le = np.zeros(2*nx*ny)
-        de = np.zeros(2*nx*ny)
-        angle_e = np.zeros(2*nx*ny)
-
-        trisk_deg = indices(2*nx*ny)
-        trisk = indices((2*nx*ny,4)) # 4 TRiSK coefs per edge
-        wee = np.zeros((2*nx*ny,4))
-
-        lon_i = indices(nx*ny)
-        lon_v = indices(nx*ny)
-        lon_e = indices(2*nx*ny)
-        lat_i = indices(nx*ny)
-        lat_v = indices(nx*ny)
-        lat_e = indices(2*nx*ny)
-
-        for x in range(nx):
-            for y in range(ny):
-                # NB : Fortran indices start at 1
-                # primal cells
-                put(index(x,y)-1, primal_deg,(primal_edge,primal_vertex,primal_ne),
-                    ((indexu(x,y),index(x,y),1), 
-                    (indexv(x,y),index(x-1,y),1),
-                    (indexu(x-1,y),index(x-1,y-1),-1),
-                    (indexv(x,y-1),index(x,y-1),-1) ))
-                # dual cells
-                put(index(x,y)-1,dual_deg,(dual_edge,dual_vertex,dual_ne,Riv2),
-                   ((indexv(x+1,y),index(x,y),1,.25),
-                    (indexu(x,y+1),index(x+1,y),-1,.25),
-                    (indexv(x,y),index(x+1,y+1),-1,.25),
-                    (indexu(x,y),index(x,y+1),1,.25) ))
-                # edges :
-                # left and right are adjacent primal cells
-                # flux is positive when going from left to right
-                # up and down are adjacent dual cells (vertices)
-                # circulation is positive when going from down to up
-                # u-points
-                ij =indexu(x,y)-1
-                left[ij]=index(x,y)
-                left[ij]=index(x,y)
-                right[ij]=index(x+1,y)
-                down[ij]=index(x,y-1)
-                up[ij]=index(x,y)
-                #le_de[ij]=dy/dx
-                le[ij]=dy
-                de[ij]=dx
-                le_de[ij]=le[ij]/de[ij]
-                angle_e[ij]=0.
-                put(ij,trisk_deg,(trisk,wee),(
-                    (indexv(x,y),.25),
-                    (indexv(x+1,y),.25),
-                    (indexv(x,y-1),.25),
-                    (indexv(x+1,y-1),.25)))
-                # v-points
-                ij = indexv(x,y)-1
-                left[ij]=index(x,y)
-                right[ij]=index(x,y+1)
-                down[ij]=index(x,y)
-                up[ij]=index(x-1,y)
-                #le_de[ij]=dx/dy
-                le[ij]=dx
-                de[ij]=dy
-                le_de[ij]=le[ij]/de[ij]
-                angle_e[ij]=.5*np.pi
-                put(ij,trisk_deg,(trisk,wee),(
-                    (indexu(x,y),-.25),
-                    (indexu(x-1,y),-.25),
-                    (indexu(x,y+1),-.25),
-                    (indexu(x-1,y+1),-.25)))
-                ij = index(x,y)-1
-                lon_i[ij], lat_i[ij] = x, y
-                ij = index(x,y)-1
-                lon_v[ij], lat_v[ij] = x+.5, y+.5
-                ij = indexu(x,y)-1
-                lon_e[ij], lat_e[ij] = x+.5, y
-                ij = indexv(x,y)-1
-                lon_e[ij], lat_e[ij] = x, y+.5
-
-        Aiv=np.zeros(nx*ny)+dx*dy  
-        Ai=Av=np.zeros(nx*ny)+dx*dy
-
-        self.llm=llm
-        self.nqdyn=nqdyn
-        self.nx=nx
-        self.ny=ny
-        self.primal_deg=primal_deg
-        self.primal_edge=primal_edge
-        self.primal_ne=primal_ne
-        self.dual_deg=dual_deg
-        self.dual_edge=dual_edge
-        self.dual_ne=dual_ne
-        self.dual_vertex=dual_vertex
-        self.primal_vertex=primal_vertex
-        self.left=left
-        self.right=right
-        self.down=down
-        self.up=up
-        self.trisk_deg=trisk_deg
-        self.trisk=trisk
-        self.Aiv=Aiv
-        self.Ai=Ai
-        self.Av=Av
-        self.le=le
-        self.de=de
-        self.angle_e=angle_e
-        self.Riv2=Riv2
-        self.wee=wee
-        self.lon_i = lon_i
-        self.lon_v = lon_v
-        self.lon_e = lon_e
-        #self.lon_ev = indices(2*nx*ny)
-        self.lat_i = lat_i
-        self.lat_v = lat_v
-        self.lat_e = lat_e
-        #self.lat_ev = indices(2*nx*ny)
-      
-
-    def field_theta(self,n=1): return zeros((n,self.nqdyn,self.ny,self.nx,self.llm))
-    def field_mass(self,n=1): return zeros((n,self.ny,self.nx,self.llm))
-    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_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 vcomp(self,u): return u[:,range(1,2*self.nx,2),:]
-
-
-    #----------------------------WRITING MESH----------------------
-    def ncwrite(self, name):
-        """The method writes Cartesian mesh on the disc.
-            Args: Mesh parameters"""
-
-        #----------------OPENING NETCDF OUTPUT FILE------------
-
-        try:
-            f = cdf.Dataset(name, 'w', format='NETCDF4')
-        except:
-            print("CartesianMesh.ncwrite : Error occurred while opening new netCDF mesh file")
-            raise
-
-        #----------------DEFINING DIMENSIONS--------------------
-
-        for dimname, dimsize in [("primal_cell", self.primal_deg.size),
-                                 ("dual_cell", self.dual_deg.size),
-                                 ("edge", self.left.size),
-                                 ("primal_edge_or_vertex", self.primal_edge.shape[1]),
-                                 ("dual_edge_or_vertex", self.dual_edge.shape[1]),
-                                 ("TWO", 2),
-                                 ("trisk_edge", 4)]:
-            f.createDimension(dimname,dimsize)
-
-        #----------------DEFINING VARIABLES---------------------
-        
-        f.description = "Cartesian_mesh"
-        f.nx, f.ny = self.nx, self.ny
-
-        def create_vars(dimname, info):
-            for vname, vtype, vdata in info:
-#                print('create_vars', dimname, vname, vdata.shape)
-                var = f.createVariable(vname,vtype,dimname)
-                var[:] = vdata
-                
-        create_vars("primal_cell", 
-                    [("primal_deg","i4",self.primal_deg), 
-                     ("Ai","f8",self.Aiv),
-                     ("lon_i","f8",self.lon_i),
-                     ("lat_i","f8",self.lat_i)] )
-        create_vars("dual_cell", 
-                    [("dual_deg","i4",self.dual_deg), 
-                     ("Av","f8",self.Aiv),
-                     ("lon_v","f8",self.lon_v),
-                     ("lat_v","f8",self.lat_v),
-                     ] )
-
-        create_vars( ("primal_cell","primal_edge_or_vertex"), 
-                     [("primal_edge", "i4", self.primal_edge),
-                      ("primal_ne", "i4", self.primal_ne),
-                      ("primal_vertex", "i4", self.primal_vertex)] )
-
-        create_vars( ("dual_cell","dual_edge_or_vertex"), 
-                     [("dual_edge", "i4", self.dual_edge),
-                      ("dual_vertex","i4",self.dual_vertex),
-                      ("dual_ne","i4",self.dual_ne),
-                      ("Riv2","f8",self.Riv2)] )
-
-        create_vars("edge",
-                    [("trisk_deg","i4",self.trisk_deg),
-                     ("le","f8",self.le),
-                     ("de","f8",self.de),
-                     ("lon_e","f8", self.lon_e),
-                     ("lat_e","f8", self.lat_e),
-                      ("angle_e","f8", self.angle_e)] )
-
-        create_vars( ("edge","TWO"),
-                     [("edge_left_right","i4", concat(self.left,self.right)),
-                      ("edge_down_up","i4", concat(self.down,self.up))] )
-                      
-
-        create_vars( ("edge","trisk_edge"),
-                    [("trisk","i4",self.trisk),
-                     ("wee","f8",self.wee)] )
-
-        f.close()
-
-
-#--------------------------------- Main program -------------------------------------
-
 parser = argparse.ArgumentParser()
 
 parser.add_argument("-nx", type=int,
-                    default=128, choices=None,
+                    default=64, choices=None,
                     help="number of x points")
 parser.add_argument("-ny", type=int,
-                    default=128, choices=None,
+                    default=64, choices=None,
                     help="number of y points")
 parser.add_argument("-Lx", type=float,
@@ -293 +26 @@
 
 dx,dy=Lx/nx,Ly/ny
-mesh = Cartesian_mesh(nx,ny,llm,nqdyn,Lx,Ly)
+mesh = meshes.Cartesian_mesh(nx,ny,llm,nqdyn,Lx,Ly,0.)
 print('Successfully initialized Cartesian mesh')
-mesh.ncwrite('cart_'+str(nx)+'.nc')
+mesh.ncwrite('cart_%03d_%03d.nc'%(nx,ny))
 print('Successfully written Cartesian mesh to NetCDF File')