Index: trunk/UTIL/PYTHON/planetoplot.py
===================================================================
--- trunk/UTIL/PYTHON/planetoplot.py	(revision 760)
+++ trunk/UTIL/PYTHON/planetoplot.py	(revision 763)
@@ -26,5 +26,4 @@
            winds=False,\
            addchar=None,\
-           interv=[0,1],\
            vmin=None,\
            vmax=None,\
@@ -54,4 +53,5 @@
            yaxis=[None,None],\
            ylog=False,\
+           xlog=False,\
            yintegral=False,\
            blat=None,\
@@ -73,5 +73,6 @@
            facwind=1.,\
            trycol=False,\
-           streamflag=False):
+           streamflag=False,\
+           analysis=None):
 
     ####################################################################################################################
@@ -95,4 +96,5 @@
     import numpy as np
     from numpy.core.defchararray import find
+    from scipy.stats import gaussian_kde,describe
     from videosink import VideoSink
     from times import sol2ls
@@ -126,21 +128,25 @@
 ### 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)
+    nlon, nlat, nvert, ntime, mapmode, nslices = determineplot(slon, slat, svert, stime, redope)
     zelen = len(namefiles)*len(var)
+### we correct number of plot fields for possible operation (substract, etc...)
+    if ope is not None:
+        if fileref is not None:       zelen = 3*len(var)*len(namefiles)
+        elif "var" in ope:            zelen = zelen + 1
     numplot = zelen*nslices
     print "********** FILES, SLICES, VARS, TOTAL PLOTS: ", len(namefiles), nslices, len(var), numplot
     print "********** MAPMODE: ", mapmode
-### we correct number of plot fields for possible operation (substract, etc...)
-    if ope is not None:
-        if fileref is not None:       zelen = zelen + 2
-        elif "var" in ope:            zelen = zelen + 1
 ### we define the arrays for plot fields -- which will be placed within multiplot loops
-    all_var  = [[]]*zelen ; all_var2  = [[]]*zelen 
-    all_title = [[]]*zelen ; all_varname = [[]]*zelen ; all_namefile = [[]]*zelen ; all_time = [[]]*zelen ; all_colorb = [[]]*zelen
+    all_var  = [[]]*zelen ; all_var2  = [[]]*zelen
+    all_title = [[]]*zelen ; all_varname = [[]]*zelen ; all_namefile = [[]]*zelen ; all_colorb = [[]]*zelen
+    all_time = [[]]*zelen ; all_vert = [[]]*zelen ; all_lat = [[]]*zelen ; all_lon = [[]]*zelen
     all_windu = [[]]*zelen ; all_windv = [[]]*zelen
-
+    plot_x = [[]]*zelen ; plot_y = [[]]*zelen ; multiplot = [[]]*zelen
+### tool (should be move to mymath)
+    getVar = lambda searchList, ind: [searchList[i] for i in ind]
 #############################
 ### LOOP OVER PLOT FIELDS ###
 #############################
+    
     for nnn in range(len(namefiles)):
      for vvv in range(len(var)): 
@@ -290,5 +296,5 @@
               ftime = np.zeros(len(time))
               for i in range(len(time)): 
-                 ftime[i] = localtime ( interv[0]+time[i]*interv[1], 0.5*(wlon[0]+wlon[1]) )
+                 ftime[i] = localtime ( time[i], slon , nc )
               time=ftime
               time=check_localtime(time)
@@ -312,10 +318,14 @@
       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] = select_getfield(zvarname=var2,znc=nc,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=lon,ylat=lat,yalt=vert,ytime=all_time[k])
+      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] = select_getfield(zvarname=all_varname[k],znc=nc,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=lon,ylat=lat,yalt=vert,ytime=all_time[k])
-    ### we inform the user about the loop then increment the loop. this is the last line of "for namefile in namefiles" 
+      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)
+
+      # last line of for namefile in namefiles
       print "**** GOT SUBDATA:",k," NAMEFILE:",namefile," VAR:",varname, var2 ; k += 1 ; firstfile = False
 
@@ -323,53 +333,92 @@
 
     ##################################
-    ### Operation on files
-    if ope is not None:
-        print "********** OPERATION: ",ope
-        if "var" not in ope:
-             if len(var) > 1: errormess("for this operation... please set only one var !")
-             if ope in ["-","+","-%"]:
-                if fileref is not None:   
-                   ncref = Dataset(fileref)
-                   all_var[k] = select_getfield(zvarname=all_varname[k-1],znc=ncref,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=lon,ylat=lat,yalt=vert,ytime=all_time[k])
-                   all_varname[k] = all_varname[k-1] ; all_time[k] = all_time[k-1] ; all_namefile[k] = all_namefile[k-1] ; all_var2[k] = all_var2[k-1] ; all_colorb[k] = all_colorb[k-1]
-                   if winds: all_windu[k] = getfield(ncref,uchar) ; all_windv[k] = getfield(ncref,vchar)
-                else: errormess("fileref is missing!") 
-                if ope == "-":     all_var[k+1]= all_var[k-1] - all_var[k]
-                elif ope == "+":   all_var[k+1]= all_var[k-1] + all_var[k]
-                elif ope == "-%":
-                    masked = np.ma.masked_where(all_var[k] == 0,all_var[k])
-                    masked.set_fill_value([np.NaN])
-                    all_var[k+1]= 100.*(all_var[k-1] - masked)/masked
-                all_varname[k+1] = all_varname[k] ; all_time[k+1] = all_time[k] ; all_namefile[k+1] = all_namefile[k] ; all_var2[k+1] = all_var2[k] ; numplot = numplot+2 
-                if len(clb) >= zelen:  all_colorb[k+1] = clb[-1]   # last additional user-defined color is for operation plot 
-                else:                  all_colorb[k+1] = "RdBu_r"  # if no additional user-defined color... set a good default one
-                if winds: all_windu[k+1] = all_windu[k-1]-all_windu[k] ; all_windv[k+1] = all_windv[k-1] - all_windv[k]
-             elif ope in ["cat"]:
-                tabtime = all_time[0];tab = all_var[0];k = 1
-                if var2: tab2 = all_var2[0]
-                while k != len(namefiles) and len(all_time[k]) != 0:
-                    if var2: tab2 = np.append(tab2,all_var2[k],axis=0) 
-                    tabtime = np.append(tabtime,all_time[k]) ; tab = np.append(tab,all_var[k],axis=0) ; k += 1
-                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:
-             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   "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")
-             numplot = numplot + 1 ; all_time[k] = all_time[k-1] ; all_namefile[k] = all_namefile[k-1]
-             all_varname[k] = all_varname[k-2] + insert + all_varname[k-1] 
-             if len(clb) >= zelen: all_colorb[k] = clb[-1]   # last additional user-defined color is for operation plot 
-             else:                 all_colorb[k] = all_colorb[k-1]  # if no additional user-defined color... set same as var
-             ### only the operation plot. do not mention colorb so that it is user-defined?
-             if "only" in ope: 
-                 numplot = 1 ; all_var[0] = all_var[k] 
-                 all_time[0] = all_time[k] ; all_namefile[0] = all_namefile[k]
-                 all_varname[0] = all_varname[k-2] + insert + all_varname[k-1]
-
+    ### Operation on files (I) with _var
+    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   "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")
+         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]
+         all_varname[k] = all_varname[k-2] + insert + all_varname[k-1]
+         if len(clb) >= zelen: all_colorb[k] = clb[-1]   # last additional user-defined color is for operation plot 
+         else:                 all_colorb[k] = all_colorb[k-1]  # if no additional user-defined color... set same as var
+         ### only the operation plot. do not mention colorb so that it is user-defined?
+         if "only" in ope:
+             numplot = 1 ; all_var[0] = all_var[k]
+             all_time[0] = all_time[k] ; all_vert[0] = all_vert[k] ; all_lat[0] = all_lat[k] ; all_lon[0] = all_lon[k] ; all_namefile[0] = all_namefile[k]
+             all_varname[0] = all_varname[k-2] + insert + all_varname[k-1]
+
+
+    ##################################
+    ### Operation on files (II) without _var
+    # we re-iterate on the plots to set operation subplots to make it compatible with multifile (using the same ref file)
+    # (k+1)%3==0 is the index of operation plots
+    # (k+2)%3==0 is the index of reference plots
+    # (k+3)%3==0 is the index of first plots
+    opefirstpass=True
+    if ope is not None and "var" not in ope:
+       print "********** OPERATION: ",ope
+       for k in np.arange(zelen):
+               if len(var) > 1: errormess("for this operation... please set only one var !")
+               if ope in ["-","+","-%","-_only","+_only","-%_only","-_histo"]:
+                  if fileref is None: errormess("fileref is missing!")
+                  else:ncref = Dataset(fileref)
+
+                  if opefirstpass: ## first plots
+                     for ll in np.arange(len(namefiles)):
+                        print "SETTING FIRST PLOT"
+                        all_varname[3*ll] = all_varname[ll] ; all_time[3*ll] = all_time[ll] ; all_vert[3*ll] = all_vert[ll] ; all_lat[3*ll] = all_lat[ll] ; all_lon[3*ll] = all_lon[ll] ; all_namefile[3*ll] = all_namefile[ll] ; all_var2[3*ll] = all_var2[ll] ; all_colorb[3*ll] = all_colorb[ll] ; all_var[3*ll] = all_var[ll]
+                        if plot_y[ll] is not None: plot_y[3*ll] = plot_y[ll] ; plot_x[3*ll] = plot_x[ll]
+                        else: plot_y[3*ll] = None ; plot_x[3*ll] = None
+                        opefirstpass=False
+
+                  if (k+2)%3==0: ## reference plots
+                        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)
+
+                  if (k+1)%3==0: ## operation plots
+                     print "SETTING OPERATION 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]
+                     if ope in ["-","-_only","-_histo"]:
+                         all_var[k]= all_var[k-2] - all_var[k-1]
+                         if plot_y[k-1] is not None and plot_y[k-2] is not None: plot_y[k] = plot_y[k-2] - plot_y[k-1]
+                         if plot_y[k-2] is None: plot_y[k] = None; plot_x[k] = None
+                     elif ope in ["+","+_only"]:   
+                         all_var[k]= all_var[k-2] + all_var[k-1]
+                         if plot_y[k-1] is not None and plot_y[k-2] is not None: plot_y[k] = plot_y[k-2] + plot_y[k-1]
+                         if plot_y[k-2] is None: plot_y[k] = None; plot_x[k] = None
+                     elif ope in ["-%","-%_only"]:
+                         masked = np.ma.masked_where(all_var[k-1] == 0,all_var[k-1])
+                         masked.set_fill_value([np.NaN])
+                         all_var[k]= 100.*(all_var[k-2] - masked)/masked
+                         if plot_y[k-1] is not None and plot_y[k-2] is not None: 
+                            masked = np.ma.masked_where(plot_y[k-1] == 0,plot_y[k-1])
+                            masked.set_fill_value([np.NaN])
+                            plot_y[k]= 100.*(plot_y[k-2] - masked)/masked
+                         if plot_y[k-2] is None: plot_y[k] = None; plot_x[k] = None
+                     if len(clb) >= zelen: all_colorb[k] = clb[-1]
+                     else: all_colorb[k] = "RdBu_r" # if no additional user-defined color... set a good default one
+                     if winds: all_windu[k] = all_windu[k-2]-all_windu[k-1] ; all_windv[k] = all_windv[k-2] - all_windv[k-1]
+
+               elif ope in ["cat"]:
+                  tabtime = all_time[0];tab = all_var[0];k = 1
+                  if var2: tab2 = all_var2[0]
+                  while k != len(namefiles) and len(all_time[k]) != 0:
+                      if var2: tab2 = np.append(tab2,all_var2[k],axis=0) 
+                      tabtime = np.append(tabtime,all_time[k]) ; tab = np.append(tab,all_var[k],axis=0) ; k += 1
+                  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")
+       if "only" in ope:
+           numplot = 1 ; all_var[0] = all_var[k]
+           all_time[0] = all_time[k] ; all_vert[0] = all_vert[k] ; all_lat[0] = all_lat[k] ; all_lon[0] = all_lon[k] ; all_namefile[0] = all_namefile[k] ; plot_x[0]=plot_x[k] ; plot_y[0]=plot_y[k]
+           all_varname[0] = all_varname[k]
 
     ##################################
@@ -377,9 +426,10 @@
     fig = figure()
     subv,subh = definesubplot( numplot, fig ) 
-    if ope in ['-','-%']: subv,subh = 2,2
+    if ope in ['-','-%','-_histo'] and len(namefiles) ==1 : subv,subh = 2,2
+    elif ope in ['-','-%'] and len(namefiles)>1 : subv, subh = len(namefiles), 3
  
     #################################
     ### Time loop for plotting device
-    nplot = 1 ; error = False 
+    nplot = 1 ; error = False ; firstpass = True 
     if lstyle is not None: linecycler = cycle(lstyle)
     else: linecycler = cycle(["-","--","-.",":"])
@@ -387,5 +437,5 @@
     while error is False:
      
-       print "********** NPLOT", nplot 
+       print "********** PLOT", nplot, " OF ",numplot
        if nplot > numplot: break
 
@@ -393,17 +443,21 @@
        ## 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 = ((nplot-1)//(nlon*nlat*nvert*ntime))%(3*len(namefiles))  ## OK only 1 var,  see test in the beginning
+           elif "var" in ope:           index_f = ((nplot-1)//(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 firstpass:
+          if len(namefiles) > 1 and len(var) == 1 and which == "unidim": pass
+          else: yeah = len(namefiles)*len(var) ; index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%yeah
+       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) 
-       if ope is not None:
-           if fileref is not None:      index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%(len(namefiles)+2)  ## OK only 1 var,  see test in the beginning
-           elif "var" in ope:           index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%(len(var)+1)        ## OK only 1 file, see test in the beginning
-           elif "cat" in ope:           index_f = 0
-       else:                            yeah = len(namefiles)*len(var) ; index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%yeah
-       time = all_time[index_f]
+       indexvert = 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)
        ltst = None 
-       if typefile in ['meso'] and indextime is not None:  ltst = localtime ( interv[0]+indextime*interv[1], 0.5*(wlon[0]+wlon[1]) )  
+       if typefile in ['meso'] and indextime is not None and len(indextime) < 2: ltst = localtime ( indextime, slon , nc)
        print "********** INDEX LON:",indexlon," LAT:",indexlat," VERT:",indexvert," TIME:",indextime
        ##var = nc.variables["phisinit"][:,:]
@@ -419,4 +473,5 @@
        error = False
        varname = all_varname[index_f]
+       which=''
        if varname:   ### what is shaded.
            what_I_plot, error = reducefield( all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \
@@ -432,5 +487,11 @@
            vecy, error = 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, \
+                                             yint=yintegral, alt=vert, anomaly=anomaly, redope=redope, mesharea=area, unidim=is1d)
+          ploty, error = reducefield( ploty, d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \
+                                             yint=yintegral, alt=vert, anomaly=anomaly, redope=redope, mesharea=area, unidim=is1d)
+          which='xy'
        #####################################################################
        #if truc:
@@ -452,5 +513,5 @@
        ####################################################################
        ### General plot settings
-       changesubplot = (numplot > 1) and (len(what_I_plot.shape) != 1)  ## default for 1D plots: superimposed. to be reworked for better flexibility.
+       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)
        colorb = all_colorb[index_f]
@@ -474,7 +535,7 @@
                        if slat is not None: lonyeah = lon2d[0,:]
                    what_I_plot, x, y = define_axis(lonyeah,latyeah,vert,time,indexlon,indexlat,indexvert,\
-                         itime,what_I_plot, len(all_var[index_f].shape),vertmode)
+                         itime,what_I_plot, len(all_var[index_f].shape),vertmode,redope)
                ###
-               if (fileref is not None) and (index_f == numplot-1):    zevmin, zevmax = calculate_bounds(what_I_plot,vmin=minop,vmax=maxop)
+               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 == numplot-1):    colorb = "RdBu_r"
@@ -493,16 +554,18 @@
                  if mrate is not None: iend=len(time)-1
                  else:                 iend=0
-                 imov = 0 
-                 if len(what_I_plot.shape) == 3:
-                    if var2:               which = "contour" ## have to start with contours rather than shading
-                    else:                  which = "regular"
+                 imov = 0
+                 if analysis in ['density','histo','fft','histodensity']: which="xy"
+                 elif len(what_I_plot.shape) == 3:
+                    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.")
                  elif len(what_I_plot.shape) == 2:
-                    if var2:               which = "contour" ## have to start with contours rather than shading
-                    else:                  which = "regular"
-                    if mrate is not None:  which = "unidim"
-                 elif len(what_I_plot.shape) == 1:
+                    if var2 and which == '':               which = "contour" ## have to start with contours rather than shading
+                    elif which == '':                      which = "regular"
+                    if mrate is not None and which == '':  which = "unidim"
+                 elif len(what_I_plot.shape) == 1 and which == '' :
                     which = "unidim"
                     if what_I_plot.shape[-1] == 1:      print "VALUE VALUE VALUE VALUE ::: ", what_I_plot[0] ; save = 'donothing'
+##                 if which == "unidim" and len(namefiles) > 1: numplot = 1 # this case is similar to several vars from one file 
                  ##### IMOV LOOP #### IMOV LOOP
                  while imov <= iend:
@@ -522,11 +585,12 @@
                     #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
-                            x, y = m(lon2d, lat2d)                                         ## this is dirty, defined above but out of imov loop
+                        m = 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:   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 which == "unidim":
-                        if lbls is not None: lbl=lbls[nplot-1] 
+#                        if lbls is not None: lbl=lbls[nplot-1] 
+                        if lbls is not None: lbl=lbls[index_f]
                         else:
                            lbl = ""
@@ -549,4 +613,50 @@
                         if indextime is None and axtime is not None and xlab is None:    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')
+
+                    elif which == "xy":
+                        if lbls is not None: lbl=lbls[index_f]
+                        else: lbl=None
+                        if analysis is not None:
+                           if analysis == 'histo': 
+                               if zelen == 1:
+                                  mpl.pyplot.hist(ploty.flatten(),bins=ndiv,normed=True, alpha=0.67, facecolor = 'green', label = lbls)
+                                  mpl.pyplot.legend()
+                                  mpl.pyplot.grid(True)
+                                  mpl.pyplot.title(zetitle)
+                               else:
+                                  multiplot[index_f]=ploty.flatten()
+                                  if index_f == zelen-1: 
+                                     if ope is not None: multiplot = getVar(multiplot,3*np.arange((index_f+1)/3)+2) ## we only compute histograms for the operation plots
+                                     mpl.pyplot.hist(multiplot,bins=ndiv,normed=True, alpha=0.75, label = lbls) 
+                                     mpl.pyplot.legend()
+                                     mpl.pyplot.grid(True)
+                                     mpl.pyplot.title(zetitle)
+                           elif analysis in ['density','histodensity']:
+                                  if ope is not None and (index_f+1)%3 !=0: pass
+                                  else:
+                                     plotx = np.linspace(min(ploty.flatten()),max(ploty.flatten()),1000)
+                                     density = 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])
+                                     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
+                                        mpl.pyplot.hist(ploty.flatten(),bins=ndiv,normed=True, alpha=0.30, label = lbl)
+                                     if index_f == zelen-1: mpl.pyplot.legend() ; mpl.pyplot.title(zetitle)
+                           else:
+                              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)
+                           if index_f == zelen-1: mpl.pyplot.legend() ; mpl.pyplot.title(zetitle)
+                        if varname == 'hodograph':
+                            a=0
+                            for ii in np.arange(len(time)): 
+                               if a%6 == 0: mpl.pyplot.text(plotx[ii],ploty[ii],time[ii]) 
+                               a=a+1
+                            mpl.pyplot.grid(True)
 
                     elif which == "regular":
@@ -590,5 +700,5 @@
 
                         if colorb not in ['nobar','onebar']:        
-                            if (fileref is not None) and (index_f == numplot-1):   daformat = "%.3f" 
+                            if (fileref is not None) and ((index_f+1)%3 == 0):   daformat = "%.3f" 
                             elif mult != 1:                                        daformat = "%.1f"
                             else:                                                  daformat = fmtvar(fvar.upper())
@@ -614,5 +724,5 @@
                                                               #200.         ## or csmooth=stride
                         ### THIS IS A QUITE SPECIFIC PIECE (does not work for mesoscale files)
-                        if ope == '-' and nplot == numplot: # this should work as long as ope is '-' guarantees 3 plots for 4 panels without contour
+                        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'
@@ -632,6 +742,7 @@
                             plot(zebins, zegauss, linewidth=1, color='r',label="mean: "+str(zemean)[0:5]+"\nstd: "+str(zestd)[0:5])
                             legend(loc=0,frameon=False)
-                            title("Histogram fig(1) "+ope+" fig(2)")
                             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":
@@ -647,5 +758,5 @@
                         if mapmode == 0:   
                             what_I_plot_frame, x, y = define_axis( lonyeah,latyeah,vert,time,indexlon,indexlat,indexvert,\
-                                                              itime,what_I_plot_frame, len(all_var2[index_f].shape),vertmode )
+                                                              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')
@@ -655,7 +766,7 @@
                             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()))
-                    if which in ["regular","unidim"]:
-
-                        if nplot > 1 and which == "unidim":
+                    if which in ["regular","unidim","xy"]:
+
+                        if nplot > 1 and which in ["unidim","xy"]:
                            pass  ## because we superimpose nplot instances
                         else:
@@ -666,5 +777,9 @@
                            if zymin is not None: mpl.pyplot.ylim(ymin=zymin)
                            if zymax is not None: mpl.pyplot.ylim(ymax=zymax)
-                           if ylog:      mpl.pyplot.semilogy()
+                           if ylog and not xlog:      mpl.pyplot.semilogy()
+                           if xlog and not ylog:      mpl.pyplot.semilogx()
+                           if xlog and ylog: 
+                                mpl.pyplot.xscale('log') 
+                                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)
@@ -715,7 +830,7 @@
        else:
             if fileref is not None:
-                if index_f is numplot-1:     plottitle = basename+' '+"fig(1) "+ope+" fig(2)"
-                elif index_f is numplot-2:   plottitle = basename+' '+fileref
-                else:                        plottitle = basename+' '+all_namefile[index_f]
+               if (index_f+1)%3 == 0:     plottitle = basename+' '+"fig(1) "+ope+" fig(2)"
+               elif (index_f+2)%3 == 0:   plottitle = basename+' '+fileref
+               else:                        plottitle = basename+' '+all_namefile[index_f]
             else:                            plottitle = basename+' '+all_namefile[index_f]
        if mult != 1:                         plottitle = '{:.0e}'.format(mult) + "*" + plottitle
@@ -724,7 +839,7 @@
           else: plottitle = zetitle[0]
           if titleref is "fill":             titleref=zetitle[index_f]
-          if fileref is not None:
-             if index_f is numplot-2:        plottitle = titleref
-             if index_f is numplot-1:        plottitle = "fig(1) "+ope+" fig(2)"
+          if fileref is not None and which != "xy":
+             if (index_f+2)%3 == 0:        plottitle = titleref
+             if (index_f+1)%3 == 0:        plottitle = "fig(1) "+ope+" fig(2)"
 #       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]))
@@ -734,4 +849,7 @@
        if nplot >= numplot: error = True
        nplot += 1
+
+       if len(namefiles) > 1 and len(var) == 1 and which == "unidim": index_f=index_f+1
+       firstpass=False
 
     if colorb == "onebar":
