Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/mymath.py
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/mymath.py	(revision 341)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/mymath.py	(revision 345)
@@ -5,5 +5,6 @@
 def max (field,axis=None):
         import numpy as np
-        return np.array(field).max(axis=axis)
+        if field is None: return None
+        else: return np.array(field).max(axis=axis)
 
 def mean (field,axis=None):
Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py	(revision 341)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py	(revision 345)
@@ -1,2 +1,3 @@
+## Author: AS
 def errormess(text,printvar=None):
     print text
@@ -5,4 +6,5 @@
     return
 
+## Author: AS
 def getname(var=False,winds=False,anomaly=False):
     if var and winds:     basename = var + '_UV'
@@ -13,4 +15,5 @@
     return basename
 
+## Author: AS
 def localtime(utc,lon):
     ltst = utc + lon / 15.
@@ -19,4 +22,5 @@
     return ltst
 
+## Author: AS
 def whatkindfile (nc):
     if 'controle' in nc.variables:   typefile = 'gcm'
@@ -25,10 +29,13 @@
     elif 'U' in nc.variables:        typefile = 'meso'
     elif 'HGT_M' in nc.variables:    typefile = 'geo'
-    else:                            errormess("whatkindfile: typefile not supported.")
+    #else:                            errormess("whatkindfile: typefile not supported.")
+    else: typefile = 'gcm' # for lslin-ed files from gcm
     return typefile
 
+## Author: AS
 def getfield (nc,var):
     ## this allows to get much faster (than simply referring to nc.variables[var])
     dimension = len(nc.variables[var].dimensions)
+    print "   Opening variable with", dimension, "dimensions ..."
     if dimension == 2:    field = nc.variables[var][:,:]
     elif dimension == 3:  field = nc.variables[var][:,:,:]
@@ -36,8 +43,10 @@
     return field
 
+## Author: AS + TN
 def reducefield (input,d4=None,d3=None,d2=None,d1=None):
     ### we do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x"
     ### it would be actually better to name d4 d3 d2 d1 as t z y x
     import numpy as np
+    from mymath import max
     dimension = np.array(input).ndim
     shape = np.array(input).shape
@@ -45,4 +54,10 @@
     output = input
     error = False
+    ### 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]
+    if d2 is not None and not isinstance(d2, np.ndarray): d2=[d2]
+    if d1 is not None and not isinstance(d1, np.ndarray): d1=[d1]
+    ### now the main part
     if dimension == 2:
         if   d2 >= shape[0]: error = True
@@ -52,19 +67,19 @@
         elif d2 is not None:         output = input[d2,:]
     elif dimension == 3:
-        if   d4 >= shape[0]: error = True
-        elif d2 >= shape[1]: error = True
-        elif d1 >= shape[2]: error = True
+        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:         output = input[d4,d2,:]
-        elif d4 is not None and d1 is not None:         output = input[d4,:,d1]
-        elif d2 is not None and d1 is not None:         output = input[:,d2,d1]
-        elif d1 is not None:                output = input[:,:,d1]
-        elif d2 is not None:                output = input[:,d2,:]
-        elif d4 is not None:                output = input[d4,:,:]
+        elif d4 is not None and d2 is not None:    output=np.mean(input[d4,:,:],axis=0); output=np.mean(output[d2,:],axis=0)
+        elif d4 is not None and d1 is not None:    output=np.mean(input[d4,:,:],0); output=np.mean(output[:,d1],1)
+        elif d2 is not None and d1 is not None:    output=np.mean(input[:,d2,:],axis=1); output=np.mean(output[:,d1],axis=1)
+        elif d1 is not None:                output = np.mean(input[:,:,d1],axis=2)
+        elif d2 is not None:                output = np.mean(input[:,d2,:],axis=1)
+        elif d4 is not None:                output = np.mean(input[d4,:,:],axis=0)
     elif dimension == 4:
-        if   d4 >= shape[0]: error = True
-        elif d3 >= shape[1]: error = True
-        elif d2 >= shape[2]: error = True
-        elif d1 >= shape[3]: error = True
+        if   max(d4) >= shape[0]: error = True
+        elif max(d3) >= shape[1]: error = True
+        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,:]
@@ -72,10 +87,10 @@
         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:                output = input[d4,d3,:,:]
-        elif d4 is not None and d2 is not None:                output = input[d4,:,d2,:]
-        elif d4 is not None and d1 is not None:                output = input[d4,:,:,d1]
-        elif d3 is not None and d2 is not None:                output = input[:,d3,d2,:]
-        elif d3 is not None and d1 is not None:                output = input[:,d3,:,d1]
-        elif d2 is not None and d1 is not None:                output = input[:,:,d2,d1]
+        elif d4 is not None and d3 is not None:  output = np.mean(input[d4,:,:,:],0); output = np.mean(output[d3,:,:],0)
+        elif d4 is not None and d2 is not None:  output = np.mean(input[d4,:,:,:],0); output = np.mean(output[:,d2,:],1)
+        elif d4 is not None and d1 is not None:  output = np.mean(input[d4,:,:,:],0); output = np.mean(output[:,:,d1],2)
+        elif d3 is not None and d2 is not None:  output = np.mean(input[:,d3,:,:],1); output = np.mean(output[:,d2,:],1)
+        elif d3 is not None and d1 is not None:  output = np.mean(input[:,d3,:,:],1); output = np.mean(output[:,:,d1],0)
+        elif d2 is not None and d1 is not None:  output = np.mean(input[:,:,d2,:],2); output = np.mean(output[:,:,d1],2)
         elif d1 is not None:                       output = input[:,:,:,d1]
         elif d2 is not None:                       output = input[:,:,d2,:]
@@ -87,40 +102,60 @@
     return output, error
 
+## Author: AS + TN
 def definesubplot ( numplot, fig ):
     from matplotlib.pyplot import rcParams
     rcParams['font.size'] = 12. ## default (important for multiple calls)
-    if   numplot == 4:
-        sub = 221
-        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
-        rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
+    if numplot <= 0:
+        subv = 99999
+        subh = 99999
+    elif numplot == 1: 
+        subv = 99999
+        subh = 99999
     elif numplot == 2:
-        sub = 121
+        subv = 1
+        subh = 2
         fig.subplots_adjust(wspace = 0.35)
         rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
     elif numplot == 3:
-        sub = 131
+        subv = 3
+        subh = 1
         fig.subplots_adjust(wspace = 0.5)
         rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
-    elif numplot == 6:
-        sub = 231
+    elif   numplot == 4:
+        subv = 2
+        subh = 2
+        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
+        rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
+    elif numplot <= 6:
+        subv = 2
+        subh = 3
         fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
         rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
-    elif numplot == 8:
-        sub = 331 #241
+    elif numplot <= 8:
+        subv = 2
+        subh = 4
         fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
         rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
-    elif numplot == 9:
-        sub = 331
+    elif numplot <= 9:
+        subv = 3
+        subh = 3
         fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
         rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
-    elif numplot == 1:
-        sub = 99999
-    elif numplot <= 0: 
-        sub = 99999
+    elif numplot <= 12:
+        subv = 3
+        subh = 4
+        fig.subplots_adjust(wspace = 0.1, hspace = 0.1)
+        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
+    elif numplot <= 16:
+        subv = 4
+        subh = 4
+        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
+        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
     else:
-        print "supported: 1,2,3,4,6,8,9"
+        print "number of plot supported: 1 to 16"
         exit()
-    return sub
-
+    return subv,subh
+
+## Author: AS
 def getstralt(nc,nvert):
     typefile = whatkindfile(nc)
@@ -140,4 +175,5 @@
     return stralt
 
+## Author: AS
 def getlschar ( namefile ):
     from netCDF4 import Dataset
@@ -169,4 +205,5 @@
     return lschar, zehour, zehourin
 
+## Author: AS
 def getprefix (nc):
     prefix = 'LMD_MMM_'
@@ -175,4 +212,5 @@
     return prefix
 
+## Author: AS
 def getproj (nc):
     typefile = whatkindfile(nc)
@@ -200,4 +238,5 @@
     return proj    
 
+## Author: AS
 def ptitle (name):
     from matplotlib.pyplot import title
@@ -205,4 +244,5 @@
     print name
 
+## Author: AS
 def polarinterv (lon2d,lat2d):
     import numpy as np
@@ -212,8 +252,10 @@
     return [wlon,wlat]
 
+## Author: AS
 def simplinterv (lon2d,lat2d):
     import numpy as np
     return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]]
 
+## Author: AS
 def wrfinterv (lon2d,lat2d):
     nx = len(lon2d[0,:])-1
@@ -231,4 +273,5 @@
     return [wlon,wlat]
 
+## Author: AS
 def makeplotres (filename,res=None,pad_inches_value=0.25,folder='',disp=True,ext='png',erase=False):
     import  matplotlib.pyplot as plt
@@ -246,4 +289,5 @@
     return
 
+## Author: AS
 def dumpbdy (field,n,stag=None):
     nx = len(field[0,:])-1
@@ -253,4 +297,5 @@
     return field[n:ny-n,n:nx-n]
 
+## Author: AS
 def getcoorddef ( nc ):   
     import numpy as np
@@ -267,4 +312,5 @@
     return lon2d,lat2d    
 
+## Author: AS
 def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False):
     import numpy as np
@@ -279,4 +325,5 @@
     return lon2d,lat2d
 
+## Author: AS
 def smooth (field, coeff):
 	## actually blur_image could work with different coeff on x and y
@@ -285,7 +332,7 @@
 	return result
 
+## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
 def gauss_kern(size, sizey=None):
 	import numpy as np
-	## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth	
     	# Returns a normalized 2D gauss kernel array for convolutions
     	size = int(size)
@@ -298,7 +345,7 @@
 	return g / g.sum()
 
+## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
 def blur_image(im, n, ny=None) :
 	from scipy.signal import convolve
-	## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
 	# blurs the image by convolving with a gaussian kernel of typical size n. 
 	# The optional keyword argument ny allows for a different size in the y direction.
@@ -307,4 +354,5 @@
     	return improc
 
+## Author: AS
 def getwinddef (nc):    
     ## getwinds for predefined types
@@ -322,4 +370,5 @@
     return uchar,vchar,metwind
 
+## Author: AS
 def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True):
     ## scale regle la reference du vecteur
@@ -344,4 +393,5 @@
     return 
 
+## Author: AS
 def display (name):
     from os import system
@@ -349,4 +399,5 @@
     return name
 
+## Author: AS
 def findstep (wlon):
     steplon = int((wlon[1]-wlon[0])/4.)  #3
@@ -361,5 +412,6 @@
     return step
 
-def define_proj (char,wlon,wlat,back=None):
+## Author: AS
+def define_proj (char,wlon,wlat,back=None,blat=False):
     from    mpl_toolkits.basemap            import Basemap
     import  numpy                           as np
@@ -368,8 +420,9 @@
     meanlon = 0.5*(wlon[0]+wlon[1])
     meanlat = 0.5*(wlat[0]+wlat[1])
-    if   wlat[0] >= 80.:   blat =  40. 
-    elif wlat[1] <= -80.:  blat = -40.
-    elif wlat[1] >= 0.:    blat = wlat[0] 
-    elif wlat[0] <= 0.:    blat = wlat[1]
+    if not blat:
+        if   wlat[0] >= 80.:   blat =  40. 
+        elif wlat[1] <= -80.:  blat = -40.
+        elif wlat[1] >= 0.:    blat = wlat[0] 
+        elif wlat[0] <= 0.:    blat = wlat[1]
     print "blat ", blat
     h = 50.  ## en km
@@ -407,4 +460,5 @@
     return m
 
+## Author: AS
 #### test temporaire
 def putpoints (map,plot):
@@ -425,4 +479,5 @@
     return
 
+## Author: AS
 def calculate_bounds(field,vmin=None,vmax=None):
     import numpy as np
@@ -448,4 +503,5 @@
     return zevmin, zevmax
 
+## Author: AS
 def bounds(what_I_plot,zevmin,zevmax):
     from mymath import max,min,mean
@@ -461,4 +517,5 @@
     return what_I_plot
 
+## Author: AS
 def nolow(what_I_plot):
     from mymath import max,min
@@ -468,4 +525,5 @@
     return what_I_plot
 
+## Author: AS
 def zoomset (wlon,wlat,zoom):
     dlon = abs(wlon[1]-wlon[0])/2.
@@ -476,4 +534,5 @@
     return wlon,wlat
 
+## Author: AS
 def fmtvar (whichvar="def"):
     fmtvar    =     { \
@@ -489,5 +548,5 @@
              "TAU_ICE":      "%.2f",\
              "VMR_ICE":      "%.1e",\
-             "MTOT":         "%.0f",\
+             "MTOT":         "%.1f",\
              "anomaly":      "%.1f",\
              "W":            "%.1f",\
@@ -497,4 +556,11 @@
              "ALBBARE":      "%.2f",\
              "TAU":          "%.1f",\
+             #### T.N.
+             "TEMP":         "%.0f",\
+             "VMR_H2OICE":   "%.0f",\
+             "VMR_H2OVAP":   "%.0f",\
+             "TAUTES":       "%.2f",\
+             "TAUTESAP":     "%.2f",\
+
                     }
     if whichvar not in fmtvar:
@@ -502,4 +568,5 @@
     return fmtvar[whichvar]
 
+## Author: AS
 ####################################################################################################################
 ### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
@@ -514,5 +581,8 @@
              "HFX":          "RdYlBu",\
              "ICETOT":       "YlGnBu_r",\
-             "MTOT":         "PuBu",\
+             #"MTOT":         "PuBu",\
+             "CCNQ":         "YlOrBr",\
+             "CCNN":         "YlOrBr",\
+             "TEMP":         "Jet",\
              "TAU_ICE":      "Blues",\
              "VMR_ICE":      "Blues",\
@@ -523,4 +593,9 @@
              "ALBBARE":      "spectral",\
              "TAU":          "YlOrBr_r",\
+             #### T.N.
+             "MTOT":         "Jet",\
+             "H2O_ICE_S":    "RdBu",\
+             "VMR_H2OICE":   "PuBu",\
+             "VMR_H2OVAP":   "PuBu",\
                      }
 #W --> spectral ou jet
@@ -531,4 +606,5 @@
     return whichcolorb[whichone]
 
+## Author: AS
 def definecolorvec (whichone="def"):
         whichcolor =    { \
@@ -549,4 +625,5 @@
         return whichcolor[whichone]
 
+## Author: AS
 def marsmap (whichone="vishires"):
         from os import uname
@@ -588,4 +665,5 @@
 #	return whichlink
 
+## Author: AS
 def latinterv (area="Whole"):
     list =    { \
@@ -609,2 +687,135 @@
     return olon,olat
 
+## Author: TN
+def separatenames (name):
+  from numpy import concatenate
+  # look for comas in the input name to separate different names (files, variables,etc ..)
+  if name is None:
+     names = None
+  else:
+    names = []
+    stop = 0
+    currentname = name
+    while stop == 0:
+      indexvir = currentname.find(',')
+      if indexvir == -1:
+        stop = 1
+        name1 = currentname
+      else:
+        name1 = currentname[0:indexvir]
+      names = concatenate((names,[name1]))
+      currentname = currentname[indexvir+1:len(currentname)]
+  return names
+
+## Author: TN [old]
+def getopmatrix (kind,n):
+  import numpy as np
+  # compute matrix of operations between files
+  # the matrix is 'number of files'-square
+  # 1: difference (row minus column), 2: add
+  # not 0 in diag : just plot
+  if n == 1:
+    opm = np.eye(1)
+  elif kind == 'basic':
+    opm = np.eye(n)
+  elif kind == 'difference1': # show differences with 1st file
+    opm = np.zeros((n,n))
+    opm[0,:] = 1
+    opm[0,0] = 0
+  elif kind == 'difference2': # show differences with 1st file AND show 1st file
+    opm = np.zeros((n,n))
+    opm[0,:] = 1
+  else:
+    opm = np.eye(n)
+  return opm
+ 
+## Author: TN [old] 
+def checkcoherence (nfiles,nlat,nlon,ntime):
+  if (nfiles > 1):
+     if (nlat > 1):
+        errormess("what you asked is not possible !")
+  return 1
+
+## Author: TN
+def readslices(saxis):
+  from numpy import empty
+  if saxis == None:
+     zesaxis = None
+  else:
+     zesaxis = empty((len(saxis),2))
+     for i in range(len(saxis)):
+        a = separatenames(saxis[i])
+        if len(a) == 1:
+           zesaxis[i,:] = float(a[0])
+        else:
+           zesaxis[i,0] = float(a[0])
+           zesaxis[i,1] = float(a[1])
+           
+  return zesaxis
+
+## Author: TN
+def  getsindex(saxis,index,axis):
+# input  : all the desired slices and the good index
+# output : all indexes to be taken into account for reducing field
+  import numpy as np
+  if saxis is None:
+      zeindex = None
+  else:
+      aaa = int(np.argmin(abs(saxis[index,0] - axis)))
+      bbb = int(np.argmin(abs(saxis[index,1] - axis)))
+      [imin,imax] = np.sort(np.array([aaa,bbb]))
+      zeindex = np.array(range(imax-imin+1))+imin
+      # because -180 and 180 are the same point in longitude,
+      # we get rid of one for averaging purposes.
+      if axis[imin] == -180 and axis[imax] == 180:
+         zeindex = zeindex[0:len(zeindex)-1]
+         print "whole longitude averaging asked, so last point is not taken into account."
+  return zeindex
+     
+## Author: TN
+def define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,indextime,what_I_plot,dim0,vertmode):
+# Purpose of define_axis is to find x and y axis scales in a smart way
+# x axis priority: 1/time 2/lon 3/lat 4/vertical
+# To be improved !!!...
+   from numpy import array,swapaxes
+   x = None
+   y = None
+   count = 0
+   what_I_plot = array(what_I_plot)
+   shape = what_I_plot.shape
+   if indextime is None:
+      x = time
+      count = count+1
+   if indexlon is None:
+      if count == 0: x = lon
+      else: y = lon
+      count = count+1
+   if indexlat is None:
+      if count == 0: x = lat
+      else: y = lat
+      count = count+1
+   if indexvert is None and dim0 is 4:
+      if vertmode == 0: # vertical axis is as is (GCM grid)
+         if count == 0: x=range(len(vert))
+         else: y=range(len(vert))
+         count = count+1
+      else: # vertical axis is in kms
+         if count == 0: x = vert
+         else: y = vert
+         count = count+1
+   x = array(x)
+   y = array(y)
+   if len(shape) == 1:
+      if shape != len(x):
+         print "WARNING HERE !!!"
+         x = y
+   elif len(shape) == 2:
+      print shape[1], len(y), shape[0], len(x)
+      if shape[1] == len(y) and shape[0] == len(x) and shape[0] != shape[1]:
+         what_I_plot = swapaxes(what_I_plot,0,1)
+         print "swapaxes", what_I_plot.shape, shape
+         #x0 = x
+         #x = y
+         #y = x0
+   #print "define_axis", x, y
+   return what_I_plot,x,y
Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/planetoplot.py
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/planetoplot.py	(revision 345)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/planetoplot.py	(revision 345)
@@ -0,0 +1,528 @@
+#!/usr/bin/env python
+
+### A. Spiga  -- LMD -- 30/06/2011 to 10/07/2011
+### T.Navarro -- LMD -- 10/2011
+### Thanks to A. Colaitis for the parser trick
+
+
+####################################
+####################################
+### The main program to plot vectors
+def planetoplot (namefiles,\
+           vertmode,\
+           proj=None,\
+           back=None,\
+           target=None,
+           stride=3,\
+           numplot=2,\
+           var=None,\
+           colorb="def",\
+           winds=True,\
+           addchar=None,\
+           interv=[0,1],\
+           vmin=None,\
+           vmax=None,\
+           tile=False,\
+           zoom=None,\
+           display=True,\
+           itstep=None,\
+           hole=False,\
+           save="gui",\
+           anomaly=False,\
+           var2=None,\
+           ndiv=10,\
+           first=1,\
+           mult=1.,\
+           zetitle="fill",\
+           slon=None,\
+           slat=None,\
+           svert=None,\
+           stime=None):
+
+    ####################################################################################################################
+    ### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
+
+    #################################
+    ### Load librairies and functions
+    from netCDF4 import Dataset
+    from myplot import getcoord2d,define_proj,makeplotres,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy,\
+                       fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\
+                       zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow,\
+                       getname,localtime,polarinterv,getsindex,define_axis
+    from mymath import deg,max,min,mean
+    from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel
+    from matplotlib.cm import get_cmap
+    import numpy as np
+    from numpy.core.defchararray import find
+        
+    #whichmode = ""
+    nlon = 1 # number of longitudinal slices -- 1 is None
+    nlat = 1
+    nvert = 1
+    ntime = 1
+    nslices = 1
+    if slon is not None:
+        nslices = nslices*len(slon)
+        nlon = len(slon)
+    if slat is not None:
+        nslices = nslices*len(slat)
+        nlat = len(slat)
+    if svert is not None:
+        nslices = nslices*len(svert)
+        nvert = len(svert)
+    if stime is not None:
+        nslices = nslices*len(stime)
+        ntime = len(stime)
+    #else:
+    #    nslices = 2
+    numplot = len(namefiles)*nslices
+    
+    mapmode = 0
+    if slon is None and slat is None:
+       mapmode = 1 # in this case we plot a map, with the given projection
+    elif proj is not None:
+       print "WARNING: you specified a", proj,\
+       "projection but asked for slices", slon,"in longitude and",slat,"in latitude"
+    print "mapmode", mapmode
+       
+    
+    all_var   = [[]]*len(namefiles)
+    all_var2  = [[]]*len(namefiles)
+    all_title = [[]]*len(namefiles)
+    
+    print "len(namefiles), nslices", len(namefiles), nslices
+    print "numplot", numplot
+
+    #########################
+    ### Loop over the files initially separated by comas to be plotted on the same figure
+    k = 0
+    firstfile = True
+    for namefile in namefiles:
+    
+      ######################
+      ### Load NETCDF object
+      nc  = Dataset(namefile)
+      
+      ##################################
+      ### Initial checks and definitions
+      typefile = whatkindfile(nc)                                  ## TYPEFILE
+      if firstfile:
+         typefile0 = typefile
+      elif typefile != typefile0:
+         errormess("Not the same files !", [typefile0, typefile])
+      if var not in nc.variables: var = False                      ## VAR 
+      if winds:                                                    ## WINDS
+         [uchar,vchar,metwind] = getwinddef(nc)             
+         if uchar == 'not found': winds = False
+      if not var and not winds: errormess("please set at least winds or var",printvar=nc.variables)
+      [lon2d,lat2d] = getcoorddef(nc)                              ## COORDINATES, could be moved below
+      if proj == None:   proj = getproj(nc)                        ## PROJECTION
+
+      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!!
+
+      if firstfile:
+         ##########################
+         ### Define plot boundaries
+         ### todo: possible areas in latinterv in argument (ex: "Far_South_Pole")
+         if proj in ["npstere","spstere"]: [wlon,wlat] = polarinterv(lon2d,lat2d)
+         elif proj in ["lcc","laea"]:      [wlon,wlat] = wrfinterv(lon2d,lat2d)
+         else:                             [wlon,wlat] = simplinterv(lon2d,lat2d)
+         if zoom:                          [wlon,wlat] = zoomset(wlon,wlat,zoom) 
+
+         #########################################
+         ### Name for title and graphics save file
+         basename = getname(var=var,winds=winds,anomaly=anomaly)
+         basename = basename + getstralt(nc,nvert)  ## can be moved elsewhere for a more generic routine
+      
+      print "var", var
+      print "var2", var2
+      #print all_var
+      #print nc
+      if var: all_var[k] = getfield(nc,var)
+      if var2: all_var2[k] = getfield(nc,var2)
+      
+      print "k", k
+      print "all_var[k].shape", all_var[k].shape
+      k += 1
+      firstfile = False
+      #### End of for namefile in namefiles
+
+
+    ##################################
+    ### Open a figure and set subplots
+    fig = figure()
+    subv,subh = definesubplot( numplot, fig ) 
+ 
+    #################################
+    ### Time loop for plotting device
+    found_lct = False
+    nplot = 1
+    itime = first
+    error = False
+    if itstep is None and numplot > 0: itstep = int(24./numplot)
+    elif numplot <= 0:                 itstep = 1 
+    
+    ### Map projection
+    if mapmode == 1:
+        m = define_proj(proj,wlon,wlat,back=back)
+        x, y = m(lon2d, lat2d)
+        
+    #for nplot in range(numplot):
+    while error is False:
+       print "nplot", nplot
+       print error
+       
+       ### Which local time ?
+       ltst = localtime ( interv[0]+itime*interv[1], 0.5*(wlon[0]+wlon[1]) )
+
+       ### General plot settings
+       #print itime, int(ltst), numplot, nplot
+       if numplot >= 1: 
+           if nplot > numplot: break
+           if numplot > 1:     
+               if typefile not in ['geo']:  subplot(subv,subh,nplot)
+           found_lct = True
+       ### If only one local time is requested (numplot < 0)
+       elif numplot <= 0: 
+           if int(ltst) + numplot != 0:	
+                 itime += 1 
+                 if found_lct is True: break     ## because it means LT was found at previous iteration
+                 else:                 continue  ## continue to iterate to find the correct LT
+           else: 
+                 found_lct = True
+                 
+                 
+       ## 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)
+
+       #### Contour plot
+       if var2:
+           what_I_plot, error = reducefield(all_var2[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert )
+           what_I_plot = what_I_plot*mult
+           if not error:
+              if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_contour,6)
+              zevmin, zevmax = calculate_bounds(what_I_plot)
+              zelevels = np.linspace(zevmin,zevmax,num=ndiv)
+              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)
+              ### If we plot a 2-D field
+              if len(what_I_plot.shape) is 2:
+                  cs = contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #colors='w' )# , alpha=0.5)
+                  clabel(cs,fmt = '%1.2e')
+              ### 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)  
+           else:
+              errormess("There is an error in reducing field !")
+
+       #### Shaded plot
+       if var:
+           what_I_plot, error = reducefield(all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert )
+           what_I_plot = what_I_plot*mult
+           
+           if not error: 
+               fvar = var
+               ###
+               if anomaly:
+                   what_I_plot = 100. * ((what_I_plot / smooth(what_I_plot,12)) - 1.)
+                   fvar = 'anomaly'
+               #if mult != 1:     
+               #    fvar = str(mult) + "*" + var
+               ###
+               if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot,6)
+               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
+               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)
+                     #contourf(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()         
+                         colorbar(fraction=0.05,pad=0.03,format=fmtvar(fvar.upper()),\
+                                           ticks=np.linspace(zevmin,zevmax,min([ndiv,20])),\
+                                           extend='neither',spacing='proportional')
+                                           # both min max neither
+               ### If we plot a 1-D 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 ??)
+               else:
+                 print "WARNING!!! ",len(what_I_plot.shape),"-D PLOT NOT SUPPORTED !!!"
+                 print "field dimensions: ", what_I_plot.shape
+                 exit()
+               
+           else:
+               errormess("There is an error in reducing field !")
+ 
+    
+       ### Next subplot
+       plottitle = basename+' '+namefiles[index_f]
+       if typefile in ['mesoapi','meso']:
+            if addchar:  plottitle = plottitle + addchar + "_LT"+str(ltst)
+            else:        plottitle = plottitle + "_LT"+str(ltst)
+       if mult != 1:           plottitle = str(mult) + "*" + plottitle
+       if zetitle != "fill":   plottitle = zetitle
+#       if indexlon is not None:
+#         plottitle = plottitle + " lon: " + str(min(lon[indexlon])) +" "+ str(max(lon[indexlon]))
+#       if indexlat is not None:
+#         plottitle = plottitle + " lat: " + str(min(lat[indexlat])) +" "+ str(max(lat[indexlat]))
+#       if indexvert is not None:
+#         plottitle = plottitle + " vert: " + str(min(vert[indexvert])) +" "+ str(max(vert[indexvert]))
+#       if indextime is not None:
+#         plottitle = plottitle + " time: " + str(min(time[indextime])) +" "+ str(max(time[indextime]))
+       ptitle( plottitle )
+       itime += itstep
+       if nplot == numplot:
+          error = True
+       nplot += 1
+
+  
+
+
+
+
+
+      
+    ##########################################################################
+    ### Save the figure in a file in the data folder or an user-defined folder
+    if typefile in ['meso','mesoapi']:   prefix = getprefix(nc)
+    elif typefile in ['gcm']:            prefix = 'LMD_GCM_'
+    else:                                prefix = ''
+    ###
+    zeplot = prefix + basename 
+    if addchar:         zeplot = zeplot + addchar
+    if numplot <= 0:    zeplot = zeplot + "_LT"+str(abs(numplot))
+    ###
+    if not target:      zeplot = namefile[0:find(namefile,'wrfout')] + zeplot
+    else:               zeplot = target + "/" + zeplot  
+    ###
+    if found_lct:     
+        pad_inches_value = 0.35
+        print "save", save
+        if save == 'png': 
+            if display: makeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value) #,erase=True)  ## a miniature
+            makeplotres(zeplot,res=200.,pad_inches_value=pad_inches_value,disp=False)
+        elif save in ['eps','svg','pdf']:
+            makeplotres(zeplot,         pad_inches_value=pad_inches_value,disp=False,ext=save)
+        elif save == 'gui':
+            show()
+        else: 
+            print "save mode not supported. using gui instead."
+            show()
+    else:   print "Local time not found"
+
+    ###############
+    ### Now the end
+    return zeplot
+
+##############################
+### A specific stuff for below
+def adjust_length (tab, zelen):
+    from numpy import ones
+    if tab is None:
+        outtab = ones(zelen) * -999999
+    else:
+        if zelen != len(tab):
+            print "not enough or too much values... setting same values all variables"
+            outtab = ones(zelen) * tab[0]
+        else:
+            outtab = tab
+    return outtab
+
+###########################################################################################
+###########################################################################################
+### What is below relate to running the file as a command line executable (very convenient)
+if __name__ == "__main__":
+    import sys
+    from optparse import OptionParser    ### to be replaced by argparse
+    from api_wrapper import api_onelevel
+    from netCDF4 import Dataset
+    from myplot import getlschar, separatenames, readslices
+    from os import system
+    import numpy as np
+
+    #############################
+    ### Get options and variables
+    parser = OptionParser()
+    parser.add_option('-f', '--file',   action='append',dest='namefile', type="string",  default=None,  help='[NEEDED] name of WRF file (append). Plot files separated by comas in the same figure')
+    parser.add_option('-l', '--level',  action='store',dest='nvert',     type="float",   default=0,     help='level (def=0)(-i 2: p,mbar)(-i 3,4: z,km)')
+    parser.add_option('-p', '--proj',   action='store',dest='proj',      type="string",  default=None,  help='projection')
+    parser.add_option('-b', '--back',   action='store',dest='back',      type="string",  default=None,  help='background image (def: None)')
+    parser.add_option('-t', '--target', action='store',dest='target',    type="string",  default=None,  help='destination folder')
+    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('-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)')
+    parser.add_option('-x', '--no-vect',action='store_false',dest='winds',               default=True,  help='no wind vectors')
+    parser.add_option('-m', '--min',    action='append',dest='vmin',     type="float",   default=None,  help='bounding minimum value (append)')    
+    parser.add_option('-M', '--max',    action='append',dest='vmax',     type="float",   default=None,  help='bounding maximum value (append)') 
+    parser.add_option('-T', '--tiled',  action='store_true',dest='tile',                 default=False, help='draw a tiled plot (no blank zone)')
+    parser.add_option('-z', '--zoom',   action='store',dest='zoom',      type="float",   default=None,  help='zoom factor in %')
+    parser.add_option('-N', '--no-api', action='store_true',dest='nocall',               default=False, help='do not recreate api file')
+    parser.add_option('-d', '--display',action='store_false',dest='display',             default=True,  help='do not pop up created images')
+    parser.add_option('-e', '--itstep', action='store',dest='itstep',    type="int",     default=None,  help='stride time (def=4)')
+    parser.add_option('-H', '--hole',   action='store_true',dest='hole',                 default=False, help='holes above max and below min')
+    parser.add_option('-S', '--save',   action='store',dest='save',      type="string",  default="gui", help='save mode (png,eps,svg,pdf or gui)(def=gui)')
+    parser.add_option('-a', '--anomaly',action='store_true',dest='anomaly',              default=False, help='compute and plot relative anomaly in %')
+    parser.add_option('-w', '--with',   action='store',dest='var2',      type="string",  default=None,  help='variable contoured')
+    parser.add_option('--div',          action='store',dest='ndiv',      type="int",     default=10,    help='number of divisions in colorbar (def: 10)')
+    parser.add_option('-F', '--first',  action='store',dest='first',     type="int",     default=1,     help='first subscript to plot (def: 1)')
+    parser.add_option('--mult',         action='store',dest='mult',      type="float",   default=1.,    help='a multiplicative factor to plotted field')
+    parser.add_option('--title',        action='store',dest='zetitle',   type="string",  default="fill",help='customize the whole title')
+    #parser.add_option('-V', action='store', dest='comb',        type="float",   default=None,  help='a defined combination of variables')
+
+    ############# T.N. changes
+    #parser.add_option('-o','--operation',action='store',dest='operation',type="string",   default=None ,help='matrix of operations between files (for now see code, sorry)')
+    parser.add_option('--lat',          action='append',dest='slat',type="string",   default=None, help='slices along latitude. One value, or two values separated by comas for averaging')
+    parser.add_option('--lon',          action='append',dest='slon', type="string",   default=None, help='slices along longitude. One value, or two values separated by comas for averaging')
+    parser.add_option('--vert',         action='append',dest='svert',type="string",   default=None, help='slices along vertical axis. One value, or two values separated by comas for averaging') # quelles coordonnees ?
+    parser.add_option('--time',         action='append',dest='stime',type="string",   default=None, help='slices along time. One value, or two values separated by comas for averaging') # quelles coordonnees ?
+
+    (opt,args) = parser.parse_args()
+    if opt.namefile is None: 
+        print "I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h"
+        exit()
+    if opt.var is None and opt.anomaly is True:
+        print "Cannot ask to compute anomaly if no variable is set"
+        exit()
+    print "Options:", opt
+    
+    #listvar = '' 
+    #if opt.var is None:
+    #    zerange = [-999999]
+    #else:
+    #    zelen = len(opt.var)
+    #    zerange = range(zelen)
+    #    #if zelen == 1: listvar = opt.var[0] + ','
+    #    #else         : 
+    #    for jjj in zerange: listvar += opt.var[jjj] + ','
+    #    listvar = listvar[0:len(listvar)-1]
+    #    vmintab = adjust_length (opt.vmin, zelen)  
+    #    vmaxtab = adjust_length (opt.vmax, zelen)
+        
+    print "namefile, length", opt.namefile, len(opt.namefile)
+
+    zeslat  = readslices(opt.slat)
+    zeslon  = readslices(opt.slon)
+    zesvert = readslices(opt.svert)
+    zestime = readslices(opt.stime)
+    print "slat,zeslat", opt.slat, zeslat
+    print "slon,zeslon", opt.slon, zeslon
+    print "svert,zesvert", opt.svert, zesvert
+    print "stime,zestime", opt.stime, zestime
+
+      
+    for i in range(len(opt.namefile)):
+      for j in range(len(opt.var)):
+
+        zenamefiles = separatenames(opt.namefile[i])
+        print "zenamefiles", zenamefiles
+        
+        if opt.vmin is not None : zevmin  = opt.vmin[min(i,len(opt.vmin)-1)]
+        else: zevmin = None
+        if opt.vmax is not None : zevmax  = opt.vmax[min(i,len(opt.vmax)-1)]
+        else: zevmax = None
+        print "vmin, zevmin", opt.vmin, zevmin
+        print "vmax, zevmax", opt.vmax, zevmax
+        
+        zevar = separatenames(opt.var[j])
+        zevar = zevar[0]
+        print "var, zevar", opt.var, zevar
+        
+        #checkcoherence(len(zenamefiles),len(opt.slat),len(opt.slon),len(opt.stime))
+        
+        zefile = zenamefiles[0]
+              
+        zelevel = opt.nvert   
+        stralt = None
+        [lschar,zehour,zehourin] = getlschar ( zefile )  ## getlschar from wrfout (or simply return "" if another file)
+        #print "lschar ",lschar
+        #print "zehour ",zehour
+        #print "zehourin ",zehourin
+    
+        #####################################################
+        ### Call Fortran routines for vertical interpolations
+        if opt.interp is not None:
+            if opt.nvert is 0 and opt.interp is 4:  zelevel = 0.010
+            ### winds or no winds
+            if opt.winds            :  zefields = 'uvmet'
+            else                    :  zefields = ''
+            ### var or no var
+            #if opt.var is None      :  pass
+            if zefields == ''       :  zefields = listvar 
+            else                    :  zefields = zefields + "," + listvar 
+            if opt.var2 is not None : zefields = zefields + "," + opt.var2  
+            print zefields
+            zefile = api_onelevel (  path_to_input   = '', \
+                                     input_name      = zefile, \
+                                     fields          = zefields, \
+                                     interp_method   = opt.interp, \
+                                     onelevel        = zelevel, \
+                                     nocall          = opt.nocall )
+            print zefile
+            zelevel = 0 ## so that zelevel could play again the role of nvert
+
+
+        #############
+        ### Main call
+        name = planetoplot (zenamefiles,int(zelevel),\
+                proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=zevar,\
+                numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\
+                addchar=lschar,interv=[zehour,zehourin],vmin=zevmin,vmax=zevmax,\
+                tile=opt.tile,zoom=opt.zoom,display=opt.display,\
+                itstep=opt.itstep,hole=opt.hole,save=opt.save,\
+                anomaly=opt.anomaly,var2=opt.var2,ndiv=opt.ndiv,first=opt.first,\
+                mult=opt.mult,zetitle=opt.zetitle,\
+                slon=zeslon,slat=zeslat,svert=zesvert,stime=zestime)
+        print 'Done: '+name
+        system("rm -f to_be_erased")
+  
+    #########################################################
+    ### Generate a .sh file with the used command saved in it
+    command = ""
+    for arg in sys.argv: command = command + arg + ' '
+    name = 'pycommand'
+    f = open(name+'.sh', 'w')
+    f.write(command)
Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/winds.py
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/winds.py	(revision 341)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/winds.py	(revision 345)
@@ -82,5 +82,5 @@
     ### Open a figure and set subplots
     fig = figure()
-    sub = definesubplot( numplot, fig ) 
+    subv,subh = definesubplot( numplot, fig ) 
  
     #################################
@@ -102,5 +102,5 @@
            if nplot > numplot: break
            if numplot > 1:     
-               if typefile not in ['geo']:  subplot(sub+nplot-1)
+               if typefile not in ['geo']:  subplot(subv,subh,nplot)
            found_lct = True
        ### If only one local time is requested (numplot < 0)
