Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/gcm.py
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/gcm.py	(revision 349)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/gcm.py	(revision 350)
@@ -26,5 +26,5 @@
     parser.add_option('-s', '--stride', action='store',dest='stride',    type="int",     default=3,     help='stride vectors (def=3)')
     parser.add_option('-v', '--var',    action='append',dest='var',      type="string",  default=None,  help='variable color-shaded (append)')
-    parser.add_option('-n', '--num',    action='store',dest='numplot',   type="int",     default=2,     help='plot number (def=2)(<0: plot LT -*numplot*)')
+    parser.add_option('-n', '--num',    action='store',dest='numplot',   type="int",     default=None,  help='plot number (def=2)(<0: plot LT -*numplot*)')
     parser.add_option('-i', '--interp', action='store',dest='interp',    type="int",     default=None,  help='interpolation (2: p, 3: z-amr, 4:z-als)')
     parser.add_option('-c', '--color',  action='store',dest='colorb',    type="string",  default="def", help='change colormap (nobar: no colorbar)')
@@ -138,5 +138,5 @@
         #############
         ### Main call
-        name = planetoplot (zenamefiles,int(zelevel),\
+        name = planetoplot (zenamefiles,nvert=int(zelevel),vertmode=opt.interp,\
                 proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=zevar,\
                 numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/meso.py
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/meso.py	(revision 349)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/meso.py	(revision 350)
@@ -26,5 +26,5 @@
     parser.add_option('-s', '--stride', action='store',dest='stride',    type="int",     default=3,     help='stride vectors (def=3)')
     parser.add_option('-v', '--var',    action='append',dest='var',      type="string",  default=None,  help='variable color-shaded (append)')
-    parser.add_option('-n', '--num',    action='store',dest='numplot',   type="int",     default=2,     help='plot number (def=2)(<0: plot LT -*numplot*)')
+    parser.add_option('-n', '--num',    action='store',dest='numplot',   type="int",     default=None,  help='plot number (def=2)(<0: plot LT -*numplot*)')
     parser.add_option('-i', '--interp', action='store',dest='interp',    type="int",     default=None,  help='interpolation (2: p, 3: z-amr, 4:z-als)')
     parser.add_option('-c', '--color',  action='store',dest='colorb',    type="string",  default="def", help='change colormap (nobar: no colorbar)')
@@ -133,5 +133,5 @@
             #############
             ### Main call
-            name = planetoplot (zefile,int(zelevel),\
+            name = planetoplot (zefile,nvert=int(zelevel),vertmode=opt.interp,\
                 proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=argvar,\
                 numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py	(revision 349)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py	(revision 350)
@@ -68,5 +68,5 @@
     output = input
     error = False
-    ### this is needed to cope the case where d4,d3,d2,d1 are single integers and not arrays
+    #### this is needed to cope the case where d4,d3,d2,d1 are single integers and not arrays
     if d4 is not None and not isinstance(d4, np.ndarray): d4=[d4]
     if d3 is not None and not isinstance(d3, np.ndarray): d3=[d3]
@@ -77,12 +77,13 @@
         if   d2 >= shape[0]: error = True
         elif d1 >= shape[1]: error = True
-        elif d1 is not None and d2 is not None:  output = input[d2,d1]
-        elif d1 is not None:         output = input[:,d1]
-        elif d2 is not None:         output = input[d2,:]
+        elif d1 is not None and d2 is not None:  output = mean(input[d2,:],axis=0); output = mean(output[d1],axis=0)
+        elif d1 is not None:         output = mean(input[:,d1],axis=1)
+        elif d2 is not None:         output = mean(input[d2,:],axis=0)
     elif dimension == 3:
         if   max(d4) >= shape[0]: error = True
         elif max(d2) >= shape[1]: error = True
         elif max(d1) >= shape[2]: error = True
-        elif d4 is not None and d2 is not None and d1 is not None:  output = input[d4,d2,d1]
+        elif d4 is not None and d2 is not None and d1 is not None:  
+            output = mean(input[d4,:,:],axis=0); output = mean(output[d2,:],axis=0); output = mean(output[d1],axis=0)
         elif d4 is not None and d2 is not None:    output = mean(input[d4,:,:],axis=0); output=mean(output[d2,:],axis=0)
         elif d4 is not None and d1 is not None:    output = mean(input[d4,:,:],axis=0); output=mean(output[:,d1],axis=1)
@@ -96,9 +97,14 @@
         elif max(d2) >= shape[2]: error = True
         elif max(d1) >= shape[3]: error = True
-        elif d4 is not None and d3 is not None and d2 is not None and d1 is not None:  output = input[d4,d3,d2,d1]
-        elif d4 is not None and d3 is not None and d2 is not None:         output = input[d4,d3,d2,:]
-        elif d4 is not None and d3 is not None and d1 is not None:         output = input[d4,d3,:,d1]
-        elif d4 is not None and d2 is not None and d1 is not None:         output = input[d4,:,d2,d1]
-        elif d3 is not None and d2 is not None and d1 is not None:         output = input[:,d3,d2,d1]
+        elif d4 is not None and d3 is not None and d2 is not None and d1 is not None:  
+            output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0); output = mean(output[d2,:],axis=0); output = mean(output[d1],axis=0)
+        elif d4 is not None and d3 is not None and d2 is not None: 
+            output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0); output = mean(output[d2,:],axis=0)
+        elif d4 is not None and d3 is not None and d1 is not None: 
+            output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0); output = mean(output[:,d1],axis=1)
+        elif d4 is not None and d2 is not None and d1 is not None: 
+            output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)
+        elif d3 is not None and d2 is not None and d1 is not None: 
+            output = mean(input[:,d3,:,:],axis=1); output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1) 
         elif d4 is not None and d3 is not None:  output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0)
         elif d4 is not None and d2 is not None:  output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,d2,:],axis=1)
@@ -107,8 +113,8 @@
         elif d3 is not None and d1 is not None:  output = mean(input[:,d3,:,:],axis=1); output = mean(output[:,:,d1],axis=0)
         elif d2 is not None and d1 is not None:  output = mean(input[:,:,d2,:],axis=2); output = mean(output[:,:,d1],axis=2)
-        elif d1 is not None:                       output = input[:,:,:,d1]
-        elif d2 is not None:                       output = input[:,:,d2,:]
-        elif d3 is not None:                       output = input[:,d3,:,:]
-        elif d4 is not None:                       output = input[d4,:,:,:]
+        elif d1 is not None:                     output = mean(input[:,:,:,d1],axis=3)
+        elif d2 is not None:                     output = mean(input[:,:,d2,:],axis=2)
+        elif d3 is not None:                     output = mean(input[:,d3,:,:],axis=1)
+        elif d4 is not None:                     output = mean(input[d4,:,:,:],axis=0)
     dimension = np.array(output).ndim
     shape = np.array(output).shape
@@ -801,15 +807,19 @@
    shape = what_I_plot.shape
    if indextime is None:
+      print "axis is time"
       x = time
       count = count+1
    if indexlon is None:
+      print "axis is lon"
       if count == 0: x = lon
       else: y = lon
       count = count+1
    if indexlat is None:
+      print "axis is lat"
       if count == 0: x = lat
       else: y = lat
       count = count+1
    if indexvert is None and dim0 is 4:
+      print "axis is vert"
       if vertmode == 0: # vertical axis is as is (GCM grid)
          if count == 0: x=range(len(vert))
@@ -822,6 +832,9 @@
    x = array(x)
    y = array(y)
+   print "what_I_plot.shape", what_I_plot.shape
+   print "x.shape, y.shape", x.shape, y.shape
    if len(shape) == 1:
-      if shape != len(x):
+      print shape[0]
+      if shape[0] != len(x):
          print "WARNING HERE !!!"
          x = y
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/planetoplot.py
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/planetoplot.py	(revision 349)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/planetoplot.py	(revision 350)
@@ -7,10 +7,11 @@
 
 def planetoplot (namefiles,\
-           vertmode,\
+           nvert=0,\
+           vertmode=0,\
            proj=None,\
            back=None,\
            target=None,
            stride=3,\
-           numplot=2,\
+           numplot=None,\
            var=None,\
            colorb="def",\
@@ -59,5 +60,5 @@
     if not isinstance(namefiles, np.ndarray): namefiles = [namefiles]
     zelen = len(namefiles)
-    numplot = zelen*nslices
+    if numplot == None:  numplot = zelen*nslices
     print "len(namefiles), nslices, numplot: ", zelen, nslices, numplot
     all_var   = [[]]*zelen
@@ -97,15 +98,27 @@
 ##########################################################
       if typefile == "gcm":
-          lat = nc.variables["latitude"][:]
+          lat = nc.variables["latitude"][:] 
           lon = nc.variables["longitude"][:]
           time = nc.variables["Time"][:]
           vert = nc.variables["altitude"][:]
-          #if firstfile:
-          #   lat0 = lat
-          #elif len(lat0) != len(lat):
-          #   errormess("Not the same latitude lengths !", [len(lat0), len(lat)])
-          #elif sum((lat == lat0) == False) != 0:
-          #   errormess("Not the same latitudes !", [lat,lat0])
-          ## Faire d'autre checks sur les compatibilites entre fichiers!!
+      elif typefile in ['meso','mesoapi']:
+          if mapmode == 0:
+              if var in ['PHTOT','W']:  vertdim='BOTTOM-TOP_PATCH_END_STAG'
+              else:                     vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
+              if var in ['V']:  latdim='SOUTH-NORTH_PATCH_END_STAG'
+              else:             latdim='SOUTH-NORTH_PATCH_END_UNSTAG'
+              if var in ['U']:  londim='WEST-EAST_PATCH_END_STAG'
+              else:             londim='WEST-EAST_PATCH_END_UNSTAG'
+              lon = np.arange(0,getattr(nc,londim),1)
+              lat = np.arange(0,getattr(nc,latdim),1)
+              vert = np.arange(0,getattr(nc,vertdim),1)
+              time = np.arange(0,len(nc.variables["Times"]),1)
+       #if firstfile:
+       #   lat0 = lat
+       #elif len(lat0) != len(lat):
+       #   errormess("Not the same latitude lengths !", [len(lat0), len(lat)])
+       #elif sum((lat == lat0) == False) != 0:
+       #   errormess("Not the same latitudes !", [lat,lat0])
+       ## Faire d'autre checks sur les compatibilites entre fichiers!!
 ##########################################################
 
@@ -122,5 +135,5 @@
          ### Name for title and graphics save file
          basename = getname(var=var,winds=winds,anomaly=anomaly)
-         basename = basename + getstralt(nc,vertmode)  ## can be moved elsewhere for a more generic routine
+         basename = basename + getstralt(nc,nvert)  ## can be moved elsewhere for a more generic routine
       
       print "var, var2: ", var, var2
@@ -179,26 +192,25 @@
 
 ####################################################################
-       if typefile in ['mesoapi','meso']: 
-           indexlon = None
-           indexlat = None
-           indexvert = vertmode  
-           indextime = itime
-           nlon = 1
-           nlat = 1
-           nvert = 1
-           ntime = 1  
-           index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%len(namefiles)
-           print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime)
-       elif typefile in ['gcm']:
+       if typefile in ['meso','mesoapi'] and mapmode == 1: 
+               indextime = itime
+               indexlon = None
+               indexlat = None
+               indexvert = nvert  ## ou svert ??? 
+               nlon = 1
+               nlat = 1
+               nvert = 1
+               ntime = 1 
+       else:
            ## get all indexes to be taken into account for this subplot and then reduce field
            ## We plot 1) all lon slices 2) all lat slices 3) all vert slices 4) all time slices and then go to the next slice
-           #print "(nplot-1)%nlon", (nplot-1)%nlon
            indexlon  = getsindex(slon,(nplot-1)%nlon,lon)
-           print "indexlon", indexlon
            indexlat  = getsindex(slat,((nplot-1)//nlon)%nlat,lat)
            indexvert = getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert)  
            indextime = getsindex(stime,((nplot-1)//(nlon*nlat*nvert))%ntime,time)  
-           index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%len(namefiles)
-           #print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime)
+       index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%len(namefiles)
+       print nlon, nlat, nvert, ntime  
+       print slon, slat, svert, stime       
+       print index_f
+       print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime)
 ####################################################################
 
@@ -208,5 +220,6 @@
            #what_I_plot = what_I_plot*mult
            if not error:
-              if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot,6)
+              if mapmode == 1:
+                  if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot,6)
               zevmin, zevmax = calculate_bounds(what_I_plot)
               zelevels = np.linspace(zevmin,zevmax,num=ndiv+1) #20)
@@ -214,5 +227,5 @@
               if mapmode == 0:
                   x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
-                        indextime,what_I_plot.shape, len(all_var[index_f].shape),vertmode)
+                        indextime,what_I_plot.shape, len(all_var2[index_f].shape),vertmode)
               ### If we plot a 2-D field
               if len(what_I_plot.shape) is 2:
@@ -221,5 +234,5 @@
               ### If we plot a 1-D field
               elif len(what_I_plot.shape) is 1:
-                  plot(x,y,what_I_plot, zelevels, colors='g', linewidths = 0.33 ) #colors='w' )# , alpha=0.5)  
+                  plot(what_I_plot,x) 
            else:
               errormess("There is an error in reducing field !")
@@ -239,43 +252,32 @@
                #    fvar = str(mult) + "*" + var
                ###
-               if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot,6)
+               if mapmode == 1:
+                   if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot,6)
+               elif mapmode == 0:
+                   what_I_plot, x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
+                         indextime,what_I_plot, len(all_var[index_f].shape),vertmode)
                zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
                if colorb in ["def","nobar"]:   palette = get_cmap(name=defcolorb(fvar.upper()))
                else:                           palette = get_cmap(name=colorb)
-               ### If we plot a 2-D field
+               ##### 2D field
                if len(what_I_plot.shape) is 2:
-                 if mapmode == 0:
-                     what_I_plot, x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
-                            indextime,what_I_plot, len(all_var[index_f].shape),vertmode)
                  if not tile:
                      if not hole: what_I_plot = bounds(what_I_plot,zevmin,zevmax)
                      #zelevels = np.linspace(zevmin*(1. + 1.e-7),zevmax*(1. - 1.e-7)) #,num=20)
-                     print "what_I_plot.shape", what_I_plot.shape
-                     print "x.shape, y.shape", x.shape, y.shape
-                     zelevels = np.linspace(zevmin,zevmax,num=ndiv+1)
+                     zelevels = np.linspace(zevmin,zevmax)#,num=ndiv+1)
                      #contourf(what_I_plot, zelevels, cmap = palette )
-                     contourf(x,y,what_I_plot, zelevels, cmap = palette )
+                     contourf( x, y, what_I_plot, zelevels, cmap = palette )
                  else:
                      if hole:  what_I_plot = nolow(what_I_plot)
-                     pcolor(x,y,what_I_plot, cmap = palette, \
-                           vmin=zevmin, vmax=zevmax )
-                 if colorb != 'nobar' and var != 'HGT' :
-                         #subplot(111)
-                         #colorbar()         
+                     pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax )
+                 if colorb != 'nobar' and var != 'HGT' :        
                          colorbar(fraction=0.05,pad=0.03,format=fmtvar(fvar.upper()),\
                                            ticks=np.linspace(zevmin,zevmax,min([ndiv+1,20])),\
                                            extend='neither',spacing='proportional')
                                            # both min max neither
-               ### If we plot a 1-D field
+               ##### 1D field
                elif len(what_I_plot.shape) is 1:
-                 what_I_plot,x,y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
-                       indextime,what_I_plot,len(all_var[index_f].shape),vertmode)
-                # x = np.array(x)
-                # y = np.array(y)
-                 print "what_I_plot.shape", what_I_plot.shape
-                 print "x.shape", x.shape
                  plot(x,what_I_plot)
-               ### If we plot something that is not a 2-D or 1-D field
-               ### (maybe plot 3-D field one day or movie ??)
+               #### Other cases: (maybe plot 3-D field one day or movie ??)
                else:
                  print "WARNING!!! ",len(what_I_plot.shape),"-D PLOT NOT SUPPORTED !!!"
