Index: trunk/UTIL/PYTHON/planetoplot.py
===================================================================
--- trunk/UTIL/PYTHON/planetoplot.py	(revision 874)
+++ trunk/UTIL/PYTHON/planetoplot.py	(revision 876)
@@ -78,5 +78,6 @@
            streamflag=False,\
            nocolorb=False,\
-           analysis=None):
+           analysis=None,\
+           monster=False):
 
     ####################################################################################################################
@@ -85,27 +86,24 @@
     #################################
     ### 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,check_localtime,polarinterv,getsindex,define_axis,determineplot,readslices,bidimfind,getlschar,hole_bounds,\
-                       getdimfromvar,select_getfield,find_devil
-    from mymath import deg,max,min,mean,writeascii,fig2data,fig2img
+    import netCDF4
+
     import matplotlib as mpl
-    from matplotlib.pyplot import contour,contourf,hist, text,subplot, figure, rcParams, savefig, colorbar, \
-                                  pcolor, show, plot, clabel, title, close, legend, xlabel, axis, ylabel, \
-                                  subplots_adjust, axes, clabel
-    from matplotlib.cm import get_cmap
+    import matplotlib.pyplot
+    import matplotlib.cm
     from mpl_toolkits.basemap import cm
+    if mrate is not None: from videosink import VideoSink
+
     import numpy as np
     from numpy.core.defchararray import find
-    from scipy.stats import gaussian_kde,describe
-    from videosink import VideoSink
+    import scipy
+    if analysis in ['laplace']: import scipy.ndimage.laplace as laplace
+
+    import itertools
+    import os
+
+    import myplot
+    from mymath import deg,max,min,mean,writeascii,fig2data,fig2img
     from times import sol2ls
-    import subprocess
     #from singlet import singlet
-    from itertools import cycle
-    import os
-    from scipy import ndimage
 
 
@@ -117,4 +115,7 @@
     print "********** WELCOME TO PLANETOPLOT **********"
     print "********************************************"
+    if monster: 
+        print "********* SPEED MODE. EXPERIMENTAL. ********"
+        print "********************************************"
 ### we ensure namefiles and var are arrays
     if not isinstance(namefiles, np.ndarray): namefiles = [namefiles]
@@ -132,9 +133,10 @@
     if trycol: clb = ["Greys","Blues","YlOrRd","jet","spectral","hot","RdBu","RdYlBu","Paired"] ; zetitle = clb ; var = [var[0]]*9
 ### we avoid specific cases not yet implemented
-    if mrate is not None and len(var) > 1: errormess("multivar not allowed in movies. should be fixed soon!")
+    if mrate is not None and len(var) > 1: myplot.errormess("multivar not allowed in movies. should be fixed soon!")
 ### we prepare number of plot fields [zelen] and number of plot instances [numplot] according to user choices
 ### --> we support multifile and multivar plots : files + vars separated by commas are on the same figure
-    nlon, nlat, nvert, ntime, mapmode, nslices = determineplot(slon, slat, svert, stime, redope)
+    nlon, nlat, nvert, ntime, mapmode, nslices = myplot.determineplot(slon, slat, svert, stime, redope)
     zelen = len(namefiles)*len(var)
+    if (nslices > 1 and monster): myplot.errormess("multislice + monster not supported yet. to be done soon")
 ### we have a special mode obtained by -p noproj in which lat/lon plots are just flat 2D plots
     if proj == "noproj": mapmode = 0
@@ -157,5 +159,5 @@
 ### LOOP OVER PLOT FIELDS ###
 #############################
-    
+   
     for nnn in range(len(namefiles)):
      for vvv in range(len(var)): 
@@ -163,43 +165,43 @@
     ### we load each NETCDF objects in namefiles
       namefile = namefiles[nnn] 
-      nc  = Dataset(namefile)
+      nc  = netCDF4.Dataset(namefile)
     ### we explore the variables in the file
       varinfile = nc.variables.keys()
       if seevar: print varinfile ; exit() 
     ### we define the type of file we have (gcm, meso, etc...)
-      typefile = whatkindfile(nc)
+      typefile = myplot.whatkindfile(nc)
       if firstfile:                 typefile0 = typefile
-      elif typefile != typefile0:   errormess("Not the same kind of files !", [typefile0, typefile])
+      elif typefile != typefile0:   myplot.errormess("Not the same kind of files !", [typefile0, typefile])
     ### we care for input file being 1D simulations
       is1d=999 
-      if "longitude" in nc.dimensions and "latitude" in nc.dimensions: is1d = len(nc.variables["longitude"][:])*len(nc.variables["latitude"][:])
-      elif "lon" in nc.dimensions and "lat" in nc.dimensions: is1d = len(nc.variables["lon"][:])*len(nc.variables["lat"][:])
+      if "longitude" in nc.dimensions and "latitude" in nc.dimensions: is1d = len(nc.dimensions["longitude"])*len(nc.dimensions["latitude"])
+      elif "lon" in nc.dimensions and "lat" in nc.dimensions: is1d = len(nc.dimensions["lon"])*len(nc.dimensions["lat"])
       if typefile in ['gcm','earthgcm'] and is1d == 1:       mapmode=0 ; winds=False
     ### we create default vert and time prescriptions if not here in case mapping mode is on (lat/lon maps)
       if redope is None and mapmode == 1:
-          if svert is None:  svert = readslices(str(level)) ; nvert=1
+          if svert is None:  svert = myplot.readslices(str(level)) ; nvert=1
           if stime is None and mrate is None:
-             stime = readslices(str(0)) ; ntime=1 ## this is a default choice
+             stime = myplot.readslices(str(0)) ; ntime=1 ## this is a default choice
              print "WELL... nothing about time axis. I took default: first time reference stored in file."
     ### we get the names of variables to be read. in case only one available, we choose this one.
     ### (we take out of the test specific names e.g. UV is not in the file but used to ask a wind speed computation)
-      varname = select_getfield(zvarname=var[vvv],znc=nc,ztypefile=typefile,mode='check')
+      varname = myplot.select_getfield(zvarname=var[vvv],znc=nc,ztypefile=typefile,mode='check')
     ### we get the names of wind variables to be read (if any)
       if winds:                                                    
-         [uchar,vchar,metwind] = getwinddef(nc)             
+         [uchar,vchar,metwind] = myplot.getwinddef(nc)             
          if uchar == 'not found': winds = False
     ### we tell the user that either no var or no wind is not acceptable
-      if not varname and not winds: errormess("please set at least winds or var",printvar=nc.variables)
+      if not varname and not winds: myplot.errormess("please set at least winds or var",printvar=nc.variables)
     ### we get the coordinates lat/lon to be used
-      [lon2d,lat2d] = getcoorddef(nc)
+      [lon2d,lat2d] = myplot.getcoorddef(nc)
     ### we get an adapted map projection if none is provided by the user
-      if proj == None:   proj = getproj(nc)   
+      if proj == None:   proj = myplot.getproj(nc)   
     ### we define plot boundaries given projection or user choices
       if firstfile:
-         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)
-         elif zarea is not None:           [wlon,wlat] = latinterv(area=zarea)
+         if proj in ["npstere","spstere"]: [wlon,wlat] = myplot.polarinterv(lon2d,lat2d)
+         elif proj in ["lcc","laea"]:      [wlon,wlat] = myplot.wrfinterv(lon2d,lat2d)
+         else:                             [wlon,wlat] = myplot.simplinterv(lon2d,lat2d)
+         if zoom:                          [wlon,wlat] = myplot.zoomset(wlon,wlat,zoom)
+         elif zarea is not None:           [wlon,wlat] = myplot.latinterv(area=zarea)
 
 #############################################################
@@ -228,19 +230,33 @@
           ### --> add a line here if your reference is not present
           else: 
-              dadim = getdimfromvar(nc) ; print "No altitude found. Try to build a simple axis.",dadim
+              dadim = myplot.getdimfromvar(nc) ; print "No altitude found. Try to build a simple axis.",dadim
               if   len(dadim) == 4:  print "-- 4D field. Assume z is dim 2." ; vert = np.arange(dadim[-3])
               elif len(dadim) == 3:  print "-- 3D field. Assume z is dim 1." ; vert = [0.]
               else:                  vert = [0.]
       ### we define time vector. either it is referenced or it is guessed based on last variable's dimensions.
-          if "Time" in nc.variables:      	time = nc.variables["Time"][:]
-          elif "time_counter" in nc.variables:  time = nc.variables["time_counter"][:]/86400. #### convert from s to days
-          elif "time" in nc.variables:    	time = nc.variables["time"][:]
+          if "Time" in nc.variables:            letime = "Time"
+          elif "time_counter" in nc.variables:  letime = "time_counter"
+          elif "time" in nc.variables:          letime = "time"
           ### --> add a line here if your reference is not present
+          else:                                 letime = None
+          if letime is not None:          
+              if monster:
+                  ### very long simulation... we just retrieve 3 values for time
+                  timelength = len(nc.dimensions[letime])
+                  dafirst = nc.variables[letime][0]
+                  dalast = nc.variables[letime][timelength-1]
+                  step = nc.variables[letime][1] - dafirst
+                  time = np.arange(dafirst,dalast,step)
+              else:
+                  ### if the time record is not too long, what follows is pretty quick
+                  time = nc.variables[letime][:]
           else: 
               print "No time found. Try to build a simple axis. Assume t is dim 1." 
-              dadim = getdimfromvar(nc)
+              dadim = myplot.getdimfromvar(nc)
               if   len(dadim) == 4:  time = np.arange(dadim[-4])
               elif len(dadim) == 3:  time = np.arange(dadim[-3])
-              else:                  time = [0.] #errormess("no time axis found.")
+              else:                  time = [0.] #myplot.errormess("no time axis found.")
+          ### (SPECIFIC?)
+          if "time_counter" in nc.variables:  time = time / 86400. #### convert from s to days
           ### (SPECIFIC. convert to Ls for Martian GCM files.)
           if axtime in ["ls"]:
@@ -256,5 +272,5 @@
           ### (simply ask for subscript)
           if axtime in ["ind"]:
-              dadim = getdimfromvar(nc)
+              dadim = myplot.getdimfromvar(nc)
               if   len(dadim) == 4:  time = np.arange(dadim[-4])
               elif len(dadim) == 3:  time = np.arange(dadim[-3])
@@ -280,5 +296,5 @@
                  if slon is not None:  vlon = slon[0][iii]  ### note: slon[:][0] does not work
                  if slat is not None:  vlat = slat[0][jjj]  ### note: slon[:][0] does not work
-                 indices[iii,jjj,:] = bidimfind(lon2d,lat2d,vlon,vlat,file=iwantawhereplot)  
+                 indices[iii,jjj,:] = myplot.bidimfind(lon2d,lat2d,vlon,vlat,file=iwantawhereplot)  
                  lonp,latp = ( lon2d[indices[iii,jjj,0],indices[iii,jjj,1]] , lat2d[indices[iii,jjj,0],indices[iii,jjj,1]] )
               ### possible bug here if several --lat
@@ -290,5 +306,5 @@
           ### we get rid of boundary relaxation zone for plots. important to do that now and not before.
           if (typefile in ['meso'] and mapmode == 1):
-             if '9999' not in getattr(nc,'START_DATE'): lon2d = dumpbdy(lon2d,6) ; lat2d = dumpbdy(lat2d,6)  
+             if '9999' not in getattr(nc,'START_DATE'): lon2d = myplot.dumpbdy(lon2d,6) ; lat2d = myplot.dumpbdy(lat2d,6)  
           ### we read the keyword for vertical dimension. we take care for vertical staggering.
           if varname in ['PHTOT','W']:    vertdim='BOTTOM-TOP_PATCH_END_STAG'
@@ -304,5 +320,5 @@
           ### we define the time axis and take care of various specificities (lt, ls, sol) or issues (concatenation)
           if axtime in ["ls","sol"]:
-              lstab, soltab, lttab = getlschar ( namefile, getaxis = True )
+              lstab, soltab, lttab = myplot.getlschar ( namefile, getaxis = True )
               if axtime == "ls":      time = lstab 
               elif axtime == "sol":   time = soltab
@@ -315,9 +331,9 @@
           if axtime in ["lt"]:
               ftime = np.zeros(len(time))
-              for i in range(len(time)): ftime[i] = localtime ( time[i], slon , namefile )
-              time=ftime ; time=check_localtime(time)
+              for i in range(len(time)): ftime[i] = myplot.localtime ( time[i], slon , namefile )
+              time=ftime ; time = myplot.check_localtime(time)
               print "LOCAL TIMES.... ", time
           ### we define the vertical axis (or lack thereof) and cover possibilities for it to be altitude, pressure, geopotential. quite SPECIFIC.
-          if typefile in ['geo']:   vert = [0.] ; stime = readslices(str(0))
+          if typefile in ['geo']:   vert = [0.] ; stime = myplot.readslices(str(0))
           else:
               if vertmode is None:  vertmode=0
@@ -332,16 +348,69 @@
 ####################################################################
 
+
     ### we fill the arrays of varname, namefile, time, colorbar at the current step considered (NB: why use both k and nnn ?)
       all_varname[k] = varname
       all_namefile[k] = namefile
-      all_time[k] = time
-      all_vert[k] = vert 
-      all_lat[k] = lat
-      all_lon[k] =  lon
       all_colorb[k] = clb[vvv]
-      if var2:  all_var2[k], plot_x[k], plot_y[k] = select_getfield(zvarname=var2,znc=nc,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=all_lon[k],ylat=all_lat[k],yalt=all_vert[k],ytime=all_time[k],analysis=analysis)
-      if winds: all_windu[k] = getfield(nc,uchar) ; all_windv[k] = getfield(nc,vchar)
-    ### we fill the arrays of fields to be plotted at the current step considered
-      all_var[k], plot_x[k], plot_y[k] = select_getfield(zvarname=all_varname[k],znc=nc,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=all_lon[k],ylat=all_lat[k],yalt=all_vert[k],ytime=all_time[k],analysis=analysis)
+
+##############################################################################
+##### LARGE FILE. WE'D BETTER NOT FILL THE MEMORY WITH THE STUPID WHOLE ARRAY!
+      if monster:
+        ####################################################################
+        ## 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
+        if ope is not None:
+            if fileref is not None:      index_f = (k//(nlon*nlat*nvert*ntime))%(3*len(namefiles))  ## OK only 1 var,  see test in the beginning
+            elif "var" in ope:           index_f = (k//(nlon*nlat*nvert*ntime))%(len(var)+1)        ## OK only 1 file, see test in the beginning
+            elif "cat" in ope:           index_f = 0
+        elif not (k == 0):
+           #if len(namefiles) > 1 and len(var) == 1 and which == "unidim": pass ## TROUVER UN MOYEN POUR unidim
+           if len(namefiles) > 1 and len(var) == 1: pass
+           else: yeah = len(namefiles)*len(var) ; index_f = (k//(nlon*nlat*nvert*ntime))%yeah
+        else: yeah = len(namefiles)*len(var) ; index_f = (k//(nlon*nlat*nvert*ntime))%yeah
+
+        ilon  = myplot.getsindex(sslon,k%nlon,lon)
+        ilat  = myplot.getsindex(sslat,(k//nlon)%nlat,lat)
+        ivert = myplot.getsindex(svert,(k//(nlon*nlat))%nvert,vert)
+ 
+        if mrate is not None:                 itime = None
+        else:                                 itime = myplot.getsindex(stime,(k//(nlon*nlat*nvert))%ntime,time)
+        ltst = None
+        if typefile in ['meso'] and itime is not None and len(itime) < 2: ltst = myplot.localtime ( itime, slon , all_namefile[index_f] )
+ 
+        if ilat  is not None:  print "********** LAT : INDEX",ilat[0],  ilat[-1],  "VALUE",lat[ilat[0]],   lat[ilat[-1]]
+        else:                  ilat = range(len(lat))
+        if ilon  is not None:  print "********** LON : INDEX",ilon[0],  ilon[-1],  "VALUE",lon[ilon[0]],   lon[ilon[-1]]
+        else:                  ilon = range(len(lon))
+        if ivert is not None:  print "********** VERT: INDEX",ivert[0], ivert[-1], "VALUE",vert[ivert[0]], vert[ivert[-1]]
+        else:                  ivert = range(len(vert))
+        if itime is not None:  print "********** TIME: INDEX",itime[0], itime[-1], "VALUE",time[itime[0]], time[itime[-1]]
+        else:                  itime = range(len(time))
+
+        all_time[k] = time[itime]
+        all_vert[k] = vert[ivert]
+        all_lat[k]  = lat[ilat]
+        all_lon[k]  = lon[ilon]
+
+        all_var[k] = myplot.getfieldred(nc,all_varname[k],ilon,ilat,ivert,itime)
+        if var2:   all_var2[k] = myplot.getfieldred(nc,var2,ilon,ilat,ivert,itime)
+        if winds:  
+            all_windu[k] = myplot.getfield(nc,uchar,ilon,ilat,ivert,itime) 
+            all_windv[k] = myplot.getfield(nc,vchar,ilon,ilat,ivert,itime)
+        plot_x[k] = None ; plot_y[k] = None
+
+      else:
+        ## regular stuff: not large array.
+        all_time[k] = time
+        all_vert[k] = vert
+        all_lat[k]  = lat
+        all_lon[k]  = lon
+        if var2:  all_var2[k], plot_x[k], plot_y[k] = myplot.select_getfield(zvarname=var2,znc=nc,ztypefile=typefile,\
+                         mode='getvar',ztsat=tsat,ylon=all_lon[k],ylat=all_lat[k],yalt=all_vert[k],ytime=all_time[k],analysis=analysis)
+        if winds: all_windu[k] = myplot.getfield(nc,uchar) ; all_windv[k] = myplot.getfield(nc,vchar)
+        ### we fill the arrays of fields to be plotted at the current step considered
+        all_var[k], plot_x[k], plot_y[k] = myplot.select_getfield(zvarname=all_varname[k],znc=nc,ztypefile=typefile,\
+                         mode='getvar',ztsat=tsat,ylon=all_lon[k],ylat=all_lat[k],yalt=all_vert[k],ytime=all_time[k],analysis=analysis)
+####################################################################
 
       # last line of for namefile in namefiles
@@ -354,11 +423,11 @@
     if ope is not None and "var" in ope:
          print "********** OPERATION: ",ope
-         if len(namefiles) > 1: errormess("for this operation... please set only one file !")
-         if len(var) > 2:       errormess("not sure this works for more than 2 vars... please check.")
+         if len(namefiles) > 1: myplot.errormess("for this operation... please set only one file !")
+         if len(var) > 2:       myplot.errormess("not sure this works for more than 2 vars... please check.")
          if   "div_var" in ope: all_var[k] = all_var[k-2] / all_var[k-1] ; insert = '_div_'
          elif "mul_var" in ope: all_var[k] = all_var[k-2] * all_var[k-1] ; insert = '_mul_'
          elif "add_var" in ope: all_var[k] = all_var[k-2] + all_var[k-1] ; insert = '_add_'
          elif "sub_var" in ope: all_var[k] = all_var[k-2] - all_var[k-1] ; insert = '_sub_'
-         else:                    errormess(ope+" : non-implemented operation. Check pp.py --help")
+         else:                    myplot.errormess(ope+" : non-implemented operation. Check pp.py --help")
          plot_x[k] = None ; plot_y[k] = None
          all_time[k] = all_time[k-1] ; all_vert[k] = all_vert[k-1] ; all_lat[k] = all_lat[k-1] ; all_lon[k] = all_lon[k-1] ; all_namefile[k] = all_namefile[k-1]
@@ -383,7 +452,7 @@
        print "********** OPERATION: ",ope
        for k in np.arange(zelen):
-               if len(var) > 1: errormess("for this operation... please set only one var !")
+               if len(var) > 1: myplot.errormess("for this operation... please set only one var !")
                if ope in ["-","+","-%","-_only","+_only","-%_only","-_histo"]:
-                  if fileref is None: errormess("fileref is missing!")
+                  if fileref is None: myplot.errormess("fileref is missing!")
                   else:ncref = Dataset(fileref)
 
@@ -399,6 +468,6 @@
                         print "SETTING REFERENCE PLOT"
                         all_varname[k] = all_varname[k-1] ; all_time[k] = all_time[k-1] ; all_vert[k] = all_vert[k-1] ; all_lat[k] = all_lat[k-1] ; all_lon[k] = all_lon[k-1] ; all_namefile[k] = all_namefile[k-1] ; all_var2[k] = all_var2[k-1] ; all_colorb[k] = all_colorb[k-1]
-                        all_var[k], plot_x[k], plot_y[k] = select_getfield(zvarname=all_varname[k-1],znc=ncref,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=all_lon[k],ylat=all_lat[k],yalt=all_vert[k],ytime=all_time[k],analysis=analysis)
-                        if winds: all_windu[k] = getfield(ncref,uchar) ; all_windv[k] = getfield(ncref,vchar)
+                        all_var[k], plot_x[k], plot_y[k] = myplot.select_getfield(zvarname=all_varname[k-1],znc=ncref,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=all_lon[k],ylat=all_lat[k],yalt=all_vert[k],ytime=all_time[k],analysis=analysis)
+                        if winds: all_windu[k] = myplot.getfield(ncref,uchar) ; all_windv[k] = myplot.getfield(ncref,vchar)
 
                   if (k+1)%3==0: ## operation plots
@@ -434,5 +503,5 @@
                   all_time[0] = np.array(tabtime) ; all_var[0] = np.array(tab) ; numplot = 1
                   if var2: all_var2[0] = np.array(tab2)
-               else: errormess(ope+" : non-implemented operation. Check pp.py --help")
+               else: myplot.errormess(ope+" : non-implemented operation. Check pp.py --help")
        if "only" in ope:
            numplot = 1 ; all_var[0] = all_var[k]
@@ -442,6 +511,6 @@
     ##################################
     ### Open a figure and set subplots
-    fig = figure()
-    subv,subh = definesubplot( numplot, fig, ipreferline = (mapmode == 1) ) 
+    fig = mpl.pyplot.figure()
+    subv,subh = myplot.definesubplot( numplot, fig, ipreferline = (mapmode == 1) ) 
     if ope in ['-','-%','-_histo'] and len(namefiles) ==1 : subv,subh = 2,2
     elif ope in ['-','-%'] and len(namefiles)>1 : subv, subh = len(namefiles), 3
@@ -450,6 +519,6 @@
     ### Time loop for plotting device
     nplot = 1 ; error = False ; firstpass = True 
-    if lstyle is not None: linecycler = cycle(lstyle)
-    else: linecycler = cycle(["-","--","-.",":"])
+    if lstyle is not None: linecycler = itertools.cycle(lstyle)
+    else: linecycler = itertools.cycle(["-","--","-.",":"])
     print "********************************************"
     while error is False:
@@ -470,26 +539,22 @@
        else: yeah = len(namefiles)*len(var) ; index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%yeah
        time = all_time[index_f] ; vert = all_vert[index_f] ; lat = all_lat[index_f] ; lon = all_lon[index_f]
-       indexlon  = getsindex(sslon,(nplot-1)%nlon,lon)
-       indexlat  = getsindex(sslat,((nplot-1)//nlon)%nlat,lat)
-       indexvert = getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert)
+       indexlon  = myplot.getsindex(sslon,(nplot-1)%nlon,lon)
+       indexlat  = myplot.getsindex(sslat,((nplot-1)//nlon)%nlat,lat)
+       indexvert = myplot.getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert)
        plotx=plot_x[index_f] ; ploty=plot_y[index_f]
        if mrate is not None:                 indextime = None 
-       else:                                 indextime = getsindex(stime,((nplot-1)//(nlon*nlat*nvert))%ntime,time)
+       else:                                 indextime = myplot.getsindex(stime,((nplot-1)//(nlon*nlat*nvert))%ntime,time)
        ltst = None 
-       if typefile in ['meso'] and indextime is not None and len(indextime) < 2: ltst = localtime ( indextime, slon , all_namefile[index_f])
-
-
-       #if mapmode == 0:
-       # if xlab is None:
-       #  if indexlon is None:  xlab = "Longitude"
-       #  elif indexlat is None: xlab = "Latitude"
-       # if ylab is None:
-       #  if indexvert is None: ylab = "Altitude"
-
-
-       print "********** INDEX LON:",indexlon," LAT:",indexlat," VERT:",indexvert," TIME:",indextime
+       if typefile in ['meso'] and indextime is not None and len(indextime) < 2: ltst = myplot.localtime ( indextime, slon , all_namefile[index_f] )
+
+       if not monster:
+         if indexlat  is not None: print "********** LAT : INDEX",indexlat[0],  indexlat[-1],  "VALUE",lat[indexlat[0]],   lat[indexlat[-1]]
+         if indexlon  is not None: print "********** LON : INDEX",indexlon[0],  indexlon[-1],  "VALUE",lon[indexlon[0]],   lon[indexlon[-1]]
+         if indexvert is not None: print "********** VERT: INDEX",indexvert[0], indexvert[-1], "VALUE",vert[indexvert[0]], vert[indexvert[-1]]
+         if indextime is not None: print "********** TIME: INDEX",indextime[0], indextime[-1], "VALUE",time[indextime[0]], time[indextime[-1]]
+
        ##var = nc.variables["phisinit"][:,:]
-       ##contourf(np.transpose(var),30,cmap = get_cmap(name="Greys_r") ) ; axis('off') ; plot(indexlat,indexlon,'mx',mew=4.0,ms=20.0)
-       ##show()
+       ##contourf(np.transpose(var),30,cmap = get_cmap(name="Greys_r") ) ; mpl.pyplot.axis('off') ; plot(indexlat,indexlon,'mx',mew=4.0,ms=20.0)
+       ##mpl.pyplot.show()
        ##exit()
        #truc = True 
@@ -503,5 +568,5 @@
        which=''
        if varname:   ### what is shaded.
-           what_I_plot, error = reducefield( all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \
+           what_I_plot, error = myplot.reducefield( all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \
                                              yint=yintegral, alt=vert, anomaly=anomaly, redope=redope, mesharea=area, unidim=is1d)
            if add != 0.:      what_I_plot = what_I_plot + add
@@ -510,15 +575,15 @@
 
        if var2:      ### what is contoured.
-           what_I_plot_contour, error = reducefield( all_var2[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert, \
+           what_I_plot_contour, error = myplot.reducefield( all_var2[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert, \
                                                      yint=yintegral, alt=vert, redope=redope )
        if winds:     ### what is plot as vectors.
-           vecx, error = reducefield( all_windu[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert)
-           vecy, error = reducefield( all_windv[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert)
+           vecx, error = myplot.reducefield( all_windu[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert)
+           vecy, error = myplot.reducefield( all_windv[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert)
            if varname in [uchar,vchar]: what_I_plot = np.sqrt( np.square(vecx) + np.square(vecy) ) ; varname = "wind"
  
        if plotx is not None:
-          plotx, error = reducefield( plotx, d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \
+          plotx, error = myplot.reducefield( plotx, d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \
                                              yint=yintegral, alt=vert, anomaly=anomaly, redope=redope, mesharea=area, unidim=is1d)
-          ploty, error = reducefield( ploty, d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \
+          ploty, error = myplot.reducefield( ploty, d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \
                                              yint=yintegral, alt=vert, anomaly=anomaly, redope=redope, mesharea=area, unidim=is1d)
           which='xy'
@@ -526,5 +591,5 @@
        #if truc:
        #   nx = what_I_plot.shape[2] ; ny = what_I_plot.shape[1] ; nz = what_I_plot.shape[0] 
-       #   for k in range(nz): print k,' over ',nz ; what_I_plot[k,:,:] = what_I_plot[k,:,:] / smooth(what_I_plot[k,:,:],12)
+       #   for k in range(nz): print k,' over ',nz ; what_I_plot[k,:,:] = what_I_plot[k,:,:] / myplot.smooth(what_I_plot[k,:,:],12)
        #   for iii in range(nx):
        #    for jjj in range(ny):
@@ -534,6 +599,6 @@
        #     if np.abs(what_I_plot[0,jjj,iii]) > 1.5: 
        #         print iii,jjj,what_I_plot[0,jjj,iii],int(abs(1.-mx)*100.),int(abs(1.-mn)*100.)
-       #         plot(rel)
-       #   show()
+       #         mpl.pyplot.plot(rel)
+       #   mpl.pyplot.show()
        #   anomaly = True ### pour avoir les bons reglages plots
        #   what_I_plot = what_I_plot[0,:,:]  
@@ -543,9 +608,9 @@
        ### General plot settings
        changesubplot = (numplot > 1) and (len(what_I_plot.shape) != 1) and (which != "xy")  ## default for 1D plots: superimposed. to be reworked for better flexibility.
-       if changesubplot: subplot(subv,subh,nplot) #; subplots_adjust(wspace=0,hspace=0)
+       if changesubplot: mpl.pyplot.subplot(subv,subh,nplot) #; mpl.pyplot.subplots_adjust(wspace=0,hspace=0)
        colorb = all_colorb[index_f]
        ####################################################################
        if error:
-               errormess("There is an error in reducing field !")
+               myplot.errormess("There is an error in reducing field !")
        else:
                ticks = ndiv + 1 
@@ -565,16 +630,16 @@
                        if slon is not None or proj == "noproj": latyeah = lat2d[:,milieux]
                        if slat is not None or proj == "noproj": lonyeah = lon2d[milieuy,:]
-                   what_I_plot, x, y = define_axis(lonyeah,latyeah,vert,time,indexlon,indexlat,indexvert,\
+                   what_I_plot, x, y = myplot.define_axis(lonyeah,latyeah,vert,time,indexlon,indexlat,indexvert,\
                          itime,what_I_plot, len(all_var[index_f].shape),vertmode,redope)
                ###
-               if analysis in ['laplace']: what_I_plot = ndimage.laplace(what_I_plot)
+               if analysis in ['laplace']: what_I_plot = laplace(what_I_plot)
                ###
-               if (fileref is not None) and ((index_f+1)%3 == 0):    zevmin, zevmax = calculate_bounds(what_I_plot,vmin=minop,vmax=maxop)
-               else:                                                   zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
+               if (fileref is not None) and ((index_f+1)%3 == 0):    zevmin, zevmax = myplot.calculate_bounds(what_I_plot,vmin=minop,vmax=maxop)
+               else:                                                   zevmin, zevmax = myplot.calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
                #if (fileref is not None) and (index_f == numplot-1):    colorb = "RdBu_r"
-               if colorb in ["def","nobar","onebar"]:                  palette = get_cmap(name=defcolorb(fvar.upper()))
+               if colorb in ["def","nobar","onebar"]:                  palette = mpl.cm.get_cmap(name=myplot.defcolorb(fvar.upper()))
                elif colorb == "relief":                                palette = cm.GMT_relief
                elif colorb == "haxby":                                 palette = cm.GMT_haxby
-               else:                                                   palette = get_cmap(name=colorb)
+               else:                                                   palette = mpl.cm.get_cmap(name=colorb)
                #palette = cm.GMT_split
                #palette = cm.GMT_globe
@@ -585,5 +650,5 @@
                elif len(what_I_plot.shape) >= 4:
                  print "WARNING!!! ",len(what_I_plot.shape),"-D PLOT NOT SUPPORTED !!! dimensions: ",what_I_plot.shape
-                 errormess("Are you sure you did not forget to prescribe a dimension ?")
+                 myplot.errormess("Are you sure you did not forget to prescribe a dimension ?")
                ##### 2. HANDLE simple 1D/2D field and movies of 1D/2D fields
                else:
@@ -595,5 +660,5 @@
                     if var2 and which == '':               which = "contour" ## have to start with contours rather than shading
                     elif which == '':                      which = "regular"
-                    if mrate is None:      errormess("3D field. Use --rate RATE for movie or specify --time TIME. Exit.")
+                    if mrate is None:      myplot.errormess("3D field. Use --rate RATE for movie or specify --time TIME. Exit.")
                  elif len(what_I_plot.shape) == 2:
                     if var2 and which == '':               which = "contour" ## have to start with contours rather than shading
@@ -621,9 +686,9 @@
                     #if mrate is not None:     
                     if mapmode == 1: 
-                        m = define_proj(proj,wlon,wlat,back=back,blat=blat,blon=blon)  ## this is dirty, defined above but out of imov loop
+                        m = myplot.define_proj(proj,wlon,wlat,back=back,blat=blat,blon=blon)  ## this is dirty, defined above but out of imov loop
                         x, y = m(lon2d, lat2d)                                         ## this is dirty, defined above but out of imov loop
                     if (typefile in ['meso'] and mapmode == 1):
-                       if '9999' not in getattr(nc,'START_DATE'): what_I_plot_frame = dumpbdy(what_I_plot_frame,6,condition=True)
-#                   if typefile in ['mesoideal']:    what_I_plot_frame = dumpbdy(what_I_plot_frame,0,stag='W',condition=dumped_vert_stag)
+                       if '9999' not in getattr(nc,'START_DATE'): what_I_plot_frame = myplot.dumpbdy(what_I_plot_frame,6,condition=True)
+#                   if typefile in ['mesoideal']:    what_I_plot_frame = myplot.dumpbdy(what_I_plot_frame,0,stag='W',condition=dumped_vert_stag)
 
                     if which == "unidim":
@@ -645,9 +710,9 @@
                         else:         zemarker = None 
                         this_is_a_regular_plot = (indexvert is not None) or (indextime is None) or (indexlat is None) or (indexlon is None)
-                        if this_is_a_regular_plot:   plot(x,what_I_plot_frame,zeline,label=lbl,marker=zemarker)  ## vertical profile
-                        else:                        plot(what_I_plot_frame,x,zeline,label=lbl,marker=zemarker)  ## regular plot
+                        if this_is_a_regular_plot:   mpl.pyplot.plot(x,what_I_plot_frame,zeline,label=lbl,marker=zemarker)  ## vertical profile
+                        else:                        mpl.pyplot.plot(what_I_plot_frame,x,zeline,label=lbl,marker=zemarker)  ## regular plot
                         mpl.pyplot.grid(True)
-                        if nplot > 1: legend(loc='best')
-                        if indextime is None and axtime is not None and xlab is None:    xlabel(axtime.upper()) ## define the right label
+                        if nplot > 1: mpl.pyplot.legend(loc='best')
+                        if indextime is None and axtime is not None and xlab is None:    mpl.pyplot.xlabel(axtime.upper()) ## define the right label
                         if save == 'txt':  writeascii(np.transpose([x,np.transpose(what_I_plot)]),'profile'+str(nplot*1000+imov)+'.txt')
                         if axtime == "lt" and indextime is None:
@@ -682,10 +747,10 @@
                                   else:
                                      plotx = np.linspace(min(ploty.flatten()),max(ploty.flatten()),1000)
-                                     density = gaussian_kde(ploty.flatten())
+                                     density = scipy.stats.gaussian_kde(ploty.flatten())
    #                                  density.covariance_factor = lambda : .25  # adjust the covariance factor to change the bandwidth if needed
    #                                  density._compute_covariance()
                                         # display the mean and variance of the kde:
                                      sample = density.resample(size=20000)
-                                     n, (smin, smax), sm, sv, ss, sk = describe(sample[0])
+                                     n, (smin, smax), sm, sv, ss, sk = scipy.stats.describe(sample[0])
                                      mpl.pyplot.plot(plotx,density(plotx), label = lbl+'\nmean: '+str(sm)[0:5]+'   std: '+str(np.sqrt(sv))[0:5]+'\nskewness: '+str(ss)[0:5]+'   kurtosis: '+str(sk)[0:5])
                                      if analysis == 'histodensity':  # plot both histo and density (to assess the rightness of the kernel density estimate for exemple) and display the estimated variance
@@ -693,9 +758,9 @@
                                      if index_f == zelen-1: mpl.pyplot.legend() ; mpl.pyplot.title(zetitle)
                            else:
-                              plot(plotx,ploty,label = lbl)
+                              mpl.pyplot.plot(plotx,ploty,label = lbl)
                               if index_f == zelen-1: mpl.pyplot.legend() ; mpl.pyplot.title(zetitle)
                               mpl.pyplot.grid(True)
                         else:
-                           plot(plotx,ploty,label = lbl)
+                           mpl.pyplot.plot(plotx,ploty,label = lbl)
                            if index_f == zelen-1: mpl.pyplot.legend() ; mpl.pyplot.title(zetitle)
                         if varname == 'hodograph':
@@ -715,29 +780,29 @@
                                 print 'INFO: Using stream file',streamfile, 'for stream lines'
                                 ncstream = Dataset(streamfile)
-                                psi = getfield(ncstream,'psi')
+                                psi = myplot.getfield(ncstream,'psi')
                                 psi = psi[0,:,:,0] # time and longitude are dummy dimensions
                                 if psi.shape[1] != len(x) or psi.shape[0] != len(y):
-                                    errormess('stream file does not fit! Dimensions: '+str(psi.shape)+' '+str(x.shape)+' '+str(y.shape))
+                                    myplot.errormess('stream file does not fit! Dimensions: '+str(psi.shape)+' '+str(x.shape)+' '+str(y.shape))
                                 zelevels = np.arange(-1.e10,1.e10,1.e9)
                                 zemin = np.min(abs(zelevels))
                                 zemax = np.max(abs(zelevels))
                                 zewidth  =  (abs(zelevels)-zemin)*(5.- 0.5)/(zemax - zemin) + 0.5 # linewidth ranges from 5 to 0.5
-                                cs = contour( x,y,psi, zelevels, colors='k', linewidths = zewidth)
-                                clabel(cs, inline=True, fontsize = 4.*rcParams['font.size']/5., fmt="%1.1e")
+                                cs = mpl.pyplot.contour( x,y,psi, zelevels, colors='k', linewidths = zewidth)
+                                mpl.pyplot.clabel(cs, inline=True, fontsize = 4.*mpl.pyplot.rcParams['font.size']/5., fmt="%1.1e")
                              else:
                                 print 'WARNING: STREAM FILE',streamfile, 'DOES NOT EXIST !'
                              
-                        if hole:         what_I_plot_frame = hole_bounds(what_I_plot_frame,zevmin,zevmax)
-                        else:            what_I_plot_frame = bounds(what_I_plot_frame,zevmin,zevmax)
-                        if flagnolow:    what_I_plot_frame = nolow(what_I_plot_frame)
+                        if hole:         what_I_plot_frame = myplot.hole_bounds(what_I_plot_frame,zevmin,zevmax)
+                        else:            what_I_plot_frame = myplot.bounds(what_I_plot_frame,zevmin,zevmax)
+                        if flagnolow:    what_I_plot_frame = myplot.nolow(what_I_plot_frame)
                         if not tile:
                             #zelevels = np.linspace(zevmin*(1. + 1.e-7),zevmax*(1. - 1.e-7)) #,num=20)
                             zelevels = np.linspace(zevmin,zevmax,num=ticks)
-                            #what_I_plot_frame = smooth(what_I_plot_frame,100)
+                            #what_I_plot_frame = myplot.smooth(what_I_plot_frame,100)
                             if mapmode == 1:       m.contourf( x, y, what_I_plot_frame, zelevels, cmap = palette, alpha=trans)
-                            elif mapmode == 0:     contourf( x, y, what_I_plot_frame, zelevels, cmap = palette, alpha=trans)
+                            elif mapmode == 0:     mpl.pyplot.contourf( x, y, what_I_plot_frame, zelevels, cmap = palette, alpha=trans)
                         else:
                             if mapmode == 1:       m.pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax, alpha=trans)
-                            elif mapmode == 0:     pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax, alpha=trans)
+                            elif mapmode == 0:     mpl.pyplot.pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax, alpha=trans)
 
                         if (cross is not None or markdevil) and mapmode == 1:
@@ -748,5 +813,5 @@
                                   mpl.pyplot.plot([idx],[idy],'+k',mew=0.5,ms=18)
                             elif markdevil:
-                                idx,idy=find_devil(nc,indextime)
+                                idx,idy=myplot.find_devil(nc,indextime)
                                 idx,idy=x[idx,idy],y[idx,idy]
                                 mpl.pyplot.plot([idx],[idy],'+k',mew=0.5,ms=18)
@@ -756,10 +821,10 @@
                             if (fileref is not None) and ((index_f+1)%3 == 0):   daformat = "%.3f" 
                             elif mult != 1:                                        daformat = "%.1f"
-                            else:                                                  daformat = fmtvar(fvar.upper())
+                            else:                                                  daformat = myplot.fmtvar(fvar.upper())
                             if proj in ['moll']:  zeorientation="horizontal" ; zepad = 0.07 ; zefrac = 0.1 #zepad=0.05
                             elif proj in ['cyl']: zeorientation="vertical" ; zepad = 0.03 ; zefrac = 0.023
                             else:                 zeorientation="vertical" ; zepad = 0.03 ; zefrac = 0.05
                             if mapmode == 0:      zefrac = 0.1
-                            zecb = colorbar( fraction=zefrac,pad=zepad,format=daformat,orientation=zeorientation,\
+                            zecb = mpl.pyplot.colorbar( fraction=zefrac,pad=zepad,format=daformat,orientation=zeorientation,\
                                       ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,21])),extend='neither',spacing='proportional' ) 
                             #if zeorientation == "horizontal": 
@@ -769,5 +834,5 @@
                         if winds:
                             if typefile in ['meso']:
-                                if '9999' not in getattr(nc,'START_DATE') : [vecx_frame,vecy_frame] = [dumpbdy(vecx_frame,6,stag=uchar,condition=True), dumpbdy(vecy_frame,6,stag=vchar,condition=True)]
+                                if '9999' not in getattr(nc,'START_DATE') : [vecx_frame,vecy_frame] = [myplot.dumpbdy(vecx_frame,6,stag=uchar,condition=True), myplot.dumpbdy(vecy_frame,6,stag=vchar,condition=True)]
                                 key = True
                                 if fvar in ['UV','uv','uvmet']: key = False
@@ -775,7 +840,7 @@
                                 key = False
                             if metwind and mapmode == 1:   [vecx_frame,vecy_frame] = m.rotate_vector(vecx_frame, vecy_frame, lon2d, lat2d)
-                            if trans != 0.0:   colorvec = definecolorvec(colorb) 
-                            else:              colorvec = definecolorvec(back) 
-                            vectorfield(vecx_frame, vecy_frame, x, y, stride=stride, csmooth=2,\
+                            if trans != 0.0:   colorvec = myplot.definecolorvec(colorb) 
+                            else:              colorvec = myplot.definecolorvec(back) 
+                            myplot.vectorfield(vecx_frame, vecy_frame, x, y, stride=stride, csmooth=2,\
                                              #scale=15., factor=300., color=colorvec, key=key)
                                              scale=20., factor=250./facwind, color=colorvec, key=key)
@@ -783,6 +848,6 @@
                         ### THIS IS A QUITE SPECIFIC PIECE (does not work for mesoscale files)
                         if ope == '-_histo' and nplot == numplot: # this should work as long as ope is '-' guarantees 3 plots for 4 panels without contour
-                            subplot(subv,subh,nplot+1)
-                            rcParams["legend.fontsize"] = 'xx-large'
+                            mpl.pyplot.subplot(subv,subh,nplot+1)
+                            mpl.pyplot.rcParams["legend.fontsize"] = 'xx-large'
                             if indexlat is None:
                                 latmin = -50.; latmax = 50. # latitude range for histogram of difference
@@ -794,17 +859,17 @@
                             toplot = np.ravel(what_I_plot_frame[np.isnan(what_I_plot_frame)==False])
                             zebins = np.linspace(zevmin,zevmax,num=30)
-                            hist(toplot,bins=zebins,histtype='step',linewidth=2,color='k',normed=True)
+                            mpl.pyplot.hist(toplot,bins=zebins,histtype='step',linewidth=2,color='k',normed=True)
                             zestd = np.std(toplot);zemean = np.mean(toplot)
                             zebins = np.linspace(zevmin,zevmax,num=300)
                             zegauss = (1./(zestd * np.sqrt(2 * np.pi)) * np.exp( - (zebins - zemean)**2 / (2 * zestd**2) ) )
-                            plot(zebins, zegauss, linewidth=1, color='r',label="mean: "+str(zemean)[0:5]+"\nstd: "+str(zestd)[0:5])
-                            legend(loc=0,frameon=False)
-                            subplot(subv,subh,nplot) # go back to last plot for title of contour difference
+                            mpl.pyplot.plot(zebins, zegauss, linewidth=1, color='r',label="mean: "+str(zemean)[0:5]+"\nstd: "+str(zestd)[0:5])
+                            mpl.pyplot.legend(loc=0,frameon=False)
+                            mpl.pyplot.subplot(subv,subh,nplot) # go back to last plot for title of contour difference
                         if ope is not None and "only" not in ope: title("fig(1) "+ope+" fig(2)")
                         elif ope is not None and "only" in ope: title("fig(1) "+ope[0]+" fig(2)")
                             
                     elif which == "contour":
-                        rcParams['contour.negative_linestyle'] = 'solid' # no dashed line for negative values
-                        zevminc, zevmaxc = calculate_bounds(what_I_plot_frame, vmin=min(what_I_plot_frame), vmax=max(what_I_plot_frame))
+                        mpl.pyplot.rcParams['contour.negative_linestyle'] = 'solid' # no dashed line for negative values
+                        zevminc, zevmaxc = myplot.calculate_bounds(what_I_plot_frame, vmin=min(what_I_plot_frame), vmax=max(what_I_plot_frame))
                         zelevels = np.linspace(zevminc,zevmaxc,ticks/2) #20)
                         ### another dirty specific stuff in the wall
@@ -814,15 +879,16 @@
                         elif var2 == 'wstar':    zelevels = np.arange(0,10,1.0)
                         elif var2 == 'zmax_th':  zelevels = np.arange(0,10,2.0) ; what_I_plot_frame = what_I_plot_frame / 1000.
+                        elif var2 in ['u']:      zelevels = np.arange(-400.,400.,10.)
                         ###
                         if mapmode == 0:   
-                            what_I_plot_frame, x, y = define_axis( lonyeah,latyeah,vert,time,indexlon,indexlat,indexvert,\
+                            what_I_plot_frame, x, y = myplot.define_axis( lonyeah,latyeah,vert,time,indexlon,indexlat,indexvert,\
                                                               itime,what_I_plot_frame, len(all_var2[index_f].shape),vertmode,redope)
                             ## is this needed? only if len(all_var2[index_f].shape) != len(all_var[index_f].shape)
-                            cs = contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 0.33)#, alpha=0.5, linestyles='solid')
-                            ##cs = contour( x,y,what_I_plot_frame, zelevels, colors='w', linewidths = 0.5)#, alpha=0.5, linestyles='solid')
-                            #clabel(cs, inline=1, fontsize = 4.*rcParams['font.size']/5., fmt=fmtvar(var2.upper()))
+                            cs = mpl.pyplot.contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 0.33)#, alpha=0.5, linestyles='solid')
+                            ##cs = mpl.pyplot.contour( x,y,what_I_plot_frame, zelevels, colors='w', linewidths = 0.5)#, alpha=0.5, linestyles='solid')
+                            #mpl.pyplot.clabel(cs, inline=1, fontsize = 4.*rcParams['font.size']/5., fmt=fmtvar(var2.upper()))
                         elif mapmode == 1:  
                             cs = m.contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 0.33)#, alpha=0.5, linestyles='solid')
-                            #clabel(cs, inline=0, fontsize = rcParams['font.size'], fmt="%.0f") #fmtvar(var2.upper()))
+                            #mpl.pyplot.clabel(cs, inline=0, fontsize = mpl.pyplot.rcParams['font.size'], fmt="%.0f") #myplot.fmtvar(var2.upper()))
                     if which in ["regular","unidim","xy"]:
 
@@ -842,6 +908,6 @@
                                 mpl.pyplot.yscale('log')
                            if invert_y:  ax = mpl.pyplot.gca() ; ax.set_ylim(ax.get_ylim()[::-1])
-                           if xlab is not None: xlabel(xlab)
-                           if ylab is not None: ylabel(ylab)
+                           if xlab is not None: mpl.pyplot.xlabel(xlab)
+                           if ylab is not None: mpl.pyplot.ylabel(ylab)
 
                         if mrate is not None:
@@ -855,10 +921,10 @@
                                 moviename='movie' ;W,H = figframe.canvas.get_width_height()
                                 video = VideoSink((H,W), moviename, rate=mrate, byteorder="rgba")
-                             video.run(mframe) ; close()
+                             video.run(mframe) ; mpl.pyplot.close()
                              if imov == iend: video.close()                            
                            ### THIS IS A WEBPAGE MOVIE 
                            else:
                              nameframe = "image"+str(1000+imov)
-                             makeplotres(nameframe,res=100.,disp=False) ; close()
+                             myplot.makeplotres(nameframe,res=100.,disp=False) ; mpl.pyplot.close()
                              if imov == 0: myfile = open("zepics", 'w')
                              myfile.write("modImages["+str(imov)+"] = '"+nameframe+"_100.png';"+ '\n')
@@ -875,7 +941,7 @@
        zevarname = varname
        if redope is not None: zevarname = zevarname + "_" + redope 
-       basename = getname(var=zevarname,var2=var2,winds=winds,anomaly=anomaly)
+       basename = myplot.getname(var=zevarname,var2=var2,winds=winds,anomaly=anomaly)
        if len(what_I_plot.shape) > 3:
-           basename = basename + getstralt(nc,level) 
+           basename = basename + myplot.getstralt(nc,level) 
        if mrate is not None: basename = "movie_" + basename 
        if typefile in ['meso']:
@@ -884,5 +950,5 @@
             plottitle = basename
             ### dans le nouveau systeme time=ls,sol,lt cette ligne pourrait ne servir a rien (ou deplacer au dessus)
-            if addchar and indextime is not None:   [addchar,gogol,gogol2] = getlschar ( all_namefile[index_f] )  ;  plottitle = plottitle + addchar
+            if addchar and indextime is not None:   [addchar,gogol,gogol2] = myplot.getlschar ( all_namefile[index_f] )  ;  plottitle = plottitle + addchar
             ### en fait redope is None doit etre remplace par : n'est ni maxt ni mint
             if redope is None and ltst is not None and ( (mapmode == 0) or (proj in ["lcc","laea","merc","nsper"]) ):  plottitle = plottitle + "_LT" + str(ltst)
@@ -905,5 +971,5 @@
 #       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]))
-       if colorb != "onebar": title( plottitle )
+       if colorb != "onebar": mpl.pyplot.title( plottitle )
        if nplot >= numplot: error = True
        nplot += 1
@@ -913,6 +979,6 @@
 
     if colorb == "onebar":
-        cax = axes([0.1, 0.2, 0.8, 0.03]) # a ameliorer
-        zecb = colorbar(cax=cax, orientation="horizontal", format=fmtvar(fvar.upper()),\
+        cax = mpl.pyplot.axes([0.1, 0.2, 0.8, 0.03]) # a ameliorer
+        zecb = mpl.pyplot.colorbar(cax=cax, orientation="horizontal", format=myplot.fmtvar(fvar.upper()),\
                  ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,21])),extend='neither',spacing='proportional')
         if zetitle[0] != "fill": zecb.ax.set_xlabel(zetitle[index_f]) ; zetitle[index_f]=""
@@ -922,5 +988,5 @@
     ### Save the figure in a file in the data folder or an user-defined folder
     if outputname is None:
-       if typefile in ['meso']:   prefix = getprefix(nc)
+       if typefile in ['meso']:   prefix = myplot.getprefix(nc)
        elif typefile in ['gcm']:            prefix = 'LMD_GCM_'
        else:                                prefix = ''
@@ -942,13 +1008,13 @@
         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=resolution,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()
+            if display: myplot.makeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value) #,erase=True)  ## a miniature
+            myplot.makeplotres(zeplot,res=resolution,pad_inches_value=pad_inches_value,disp=False)
+        elif save in ['eps','svg','pdf']:     myplot.makeplotres(zeplot,pad_inches_value=pad_inches_value,disp=False,ext=save)
+        elif save == 'gui':                   mpl.pyplot.show()
         elif save == 'donothing':             pass
         elif save == 'txt':                   print "Saved results in txt file." 
         else: 
             print "INFO: save mode not supported. using gui instead."
-            show()
+            mpl.pyplot.show()
 
     ###################################
