Index: unk/UTIL/PYTHON/README.GCM
===================================================================
--- /trunk/UTIL/PYTHON/README.GCM	(revision 391)
+++ 	(revision )
@@ -1,80 +1,0 @@
-
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-PYTHON COMMAND LINE PLOTS EXAMPLES
-
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-
-***********************************************************************************
-Simple 2D plot: Zonal mean.
-**********************************************************************************
-
-Goal:
-
-Plot the zonal mean temperature from a netcdf fiel representing one month.
-
-Command:
-
-gcm.py -f POLAR_NIGHT_RUN/diagfi16.nc --var temp --lon 180,-180 --time 0,65
-
-Note:
-
-The --time, --lat, --lon and --vert command takes in input values corresponding to the unit stored in the netcdf file, and not indices ! For example, if the "Time" unit is in sol (which is common for a gcm output), --time 2 means sol 2.0 and not index 2 along the time direction. Consequently, one can ask --time 2.5 for temperatures at 12:00 on sol 2.
-
-Means are easy to perform by specifying a range. Here: --time 0,65 and --lat -180,180.
-
-***********************************************************************************
-Vertical interpolation of the field.
-***********************************************************************************
-
-Goal: 
-
-Calls to zrecast and api are built-in the python functions. One can call them using -i with the appropriate argument (see meso.py -h or gcm.py -h). Here is an example that re-interpolates data using zrecast before plotting it in a 2D contour.
-
-Command:
-
-gcm.py -f POLAR_NIGHT_RUN/diagfi16.nc --var temp --lon -180,180 --time 0,65 -i 4
-
-Note:
-
-All interpolation modes available in zrecast and api (pressure, AGL, distance from planet center, etc...) are theoretically possible, but may not be coded yet in the routine. See gcm.py -h or meso.py -h.
-
-For this example, the default behavior of zrecast for -i 4 is to interpolate in (m) from the local surface, between 0 and 150 km. The command will generate a reinterpolated netcdf file "POLAR_NIGHT_RUN/diagfi16_S.nc" with only the requested field, which is not deleted afterward.
-
-
-***********************************************************************************
-2D plot of the difference between two files.
-***********************************************************************************
-
-Goal:
-
-Comparing two .nc files with similar dimension axis can be done in a single command, by specifying which files to compare and the comparison operator (i.e. is it a difference, an addition, etc...). When comparing data along a vertical axis, it can be wise to also ask for an interpolation of the fields to make sure the comparison is correct.
-
-Command:
-
-gcm.py -f POLAR_NIGHT_RUN/stats16.nc --var temp --lon -180,180 --time 1 -i 4 --fref POLAR_NIGHT_REF/stats16.nc --operation - --mope -2 --Mope 2 --title "Polar temperatures with new parametrizations" --titleref "Reference run"
-
-Note:
-
-The command will output 3 plots: the field from file 1, the field from file 2, and the comparison between the two. One can specify specific names for the title of these plots by using --title and --titleref for the titles of file 1 and file 2, and can specify different plotting range for the normal field (-m -M) and the compared field (--mope --Mope).
-
-One can combine this command with projections and means, so that for example, to compare co2 depletion at the south pole:
-
-gcm.py -f POLAR_NIGHT_RUN/start16.nc --var co2 --vert 0,150 --proj spstere --time 1 -i 4 --fref POLAR_NIGHT_REF/start16.nc --operation - --mope -0.5 --Mope 0.5 --title "Polar co2 with new parametrizations" --titleref "Reference run"
-
-***********************************************************************************
-2D plot of data with missing values, along a pressure axis (decreasing with height).
-***********************************************************************************
-
-Goal:
-
-By default, python will force the y-axis of a 2D plot to be ordered by increasing values. Here is how to force it otherwise.
-
-Command:
-
-gcm.py -f TES.MappedClimatology.nadir.MY25.nc --var T_nadir_day --lat -76. --time 90 --ymin 500 --ymax 1 -m 128 -M 148 -H
-
-Note:
-
-The axis reversal is done by specifying ymin and ymax in the right order. One can also simply use --inverty and not specify (ymin,ymax). Missing values (out of range values) are replaced by holes by the option "-H".
Index: /trunk/UTIL/PYTHON/README.PP
===================================================================
--- /trunk/UTIL/PYTHON/README.PP	(revision 392)
+++ /trunk/UTIL/PYTHON/README.PP	(revision 392)
@@ -0,0 +1,137 @@
+**************************************
+**************************************
+**************************************
+    PLANETOPLOT TUTORIAL EXAMPLES
+**************************************
+         Authors : AC + AS
+**************************************
+  DON'T FORGET YOUR BEST FRIEND IS
+     pp.py -h [or] pp.py --help
+**************************************
+**************************************
+**************************************
+
+
+*************************************
+MAPMODE 1 
+MAPPING MODE
+SIMPLE EXAMPLES on a SAMPLE GCM FILE
+*************************************
+
+Goal: The simplest, most minimal example. Mapping near-surface wind vectors over topography.
+pp.py -f diagfired.nc
+
+Goal: I would like finer contours.
+pp.py -f diagfired.nc --div 30
+
+Goal: I would like more vectors [i.e. lower the stride].
+pp.py -f diagfired.nc --div 30 -s 1
+
+Goal: I want to get rid of those damn vectors.
+pp.py -f diagfired.nc -x
+
+Goal: I want to map a given field (surface temperature).
+pp.py -f diagfired.nc -x -v tsurf
+
+Goal: I want to map two fields next to one another (topography and tauice).
+pp.py -f diagfired.nc -x -v phisinit,tauice
+
+Goal: I want to map two fields, tauice shaded, topography contoured, same plot.
+pp.py -f diagfired.nc -x -v tauice -w phisinit
+
+Goal: I want to map a field but projected on the sphere.
+pp.py -f diagfired.nc -x -v tauice -p ortho
+
+Goal: I want to redefine the minimum and maximum values shown.
+pp.py -f diagfired.nc -x -v tauice -m 0.2 -M 0.9
+
+Goal: I want to insert holes wherever values are lower than 0.2 and higher than 0.9
+pp.py -f diagfired.nc -x -v tauice -m 0.2 -M 0.9 -H
+
+Goal: I want to fill holes with an background image of Mars [you have to be connected to Internet]
+pp.py -f diagfired.nc -x -v tauice -m 0.2 -M 0.9 -H -b vishires
+
+Goal: I want the same map, but projected on the sphere
+pp.py -f diagfired.nc -x -v tauice -m 0.2 -M 0.9 -H -b vishires -p ortho
+
+Goal: I want the same map, but projected with north polar stereographic view
+pp.py -f diagfired.nc -x -v tauice -m 0.2 -M 0.9 -H -b vishires -p npstere
+
+Goal: I want to save this in PNG format
+pp.py -f diagfired.nc -x -v tauice -m 0.2 -M 0.9 -H -b vishires -p ortho -S png
+
+Goal: I want to plot results from two simulation files next to one another
+pp.py -f diagfired.nc,diagfired.nc -x -v tsurf 
+
+***********************************************************************************
+Simple 2D plot: Zonal mean.
+**********************************************************************************
+
+Goal:
+
+Plot the zonal mean temperature from a netcdf fiel representing one month.
+
+Command:
+
+gcm.py -f POLAR_NIGHT_RUN/diagfi16.nc --var temp --lon 180,-180 --time 0,65
+
+Note:
+
+The --time, --lat, --lon and --vert command takes in input values corresponding to the unit stored in the netcdf file, and not indices ! For example, if the "Time" unit is in sol (which is common for a gcm output), --time 2 means sol 2.0 and not index 2 along the time direction. Consequently, one can ask --time 2.5 for temperatures at 12:00 on sol 2.
+
+Means are easy to perform by specifying a range. Here: --time 0,65 and --lat -180,180.
+
+***********************************************************************************
+Vertical interpolation of the field.
+***********************************************************************************
+
+Goal: 
+
+Calls to zrecast and api are built-in the python functions. One can call them using -i with the appropriate argument (see meso.py -h or gcm.py -h). Here is an example that re-interpolates data using zrecast before plotting it in a 2D contour.
+
+Command:
+
+gcm.py -f POLAR_NIGHT_RUN/diagfi16.nc --var temp --lon -180,180 --time 0,65 -i 4
+
+Note:
+
+All interpolation modes available in zrecast and api (pressure, AGL, distance from planet center, etc...) are theoretically possible, but may not be coded yet in the routine. See gcm.py -h or meso.py -h.
+
+For this example, the default behavior of zrecast for -i 4 is to interpolate in (m) from the local surface, between 0 and 150 km. The command will generate a reinterpolated netcdf file "POLAR_NIGHT_RUN/diagfi16_S.nc" with only the requested field, which is not deleted afterward.
+
+
+***********************************************************************************
+2D plot of the difference between two files.
+***********************************************************************************
+
+Goal:
+
+Comparing two .nc files with similar dimension axis can be done in a single command, by specifying which files to compare and the comparison operator (i.e. is it a difference, an addition, etc...). When comparing data along a vertical axis, it can be wise to also ask for an interpolation of the fields to make sure the comparison is correct.
+
+Command:
+
+gcm.py -f POLAR_NIGHT_RUN/stats16.nc --var temp --lon -180,180 --time 1 -i 4 --fref POLAR_NIGHT_REF/stats16.nc --operation - --mope -2 --Mope 2 --title "Polar temperatures with new parametrizations" --titleref "Reference run"
+
+Note:
+
+The command will output 3 plots: the field from file 1, the field from file 2, and the comparison between the two. One can specify specific names for the title of these plots by using --title and --titleref for the titles of file 1 and file 2, and can specify different plotting range for the normal field (-m -M) and the compared field (--mope --Mope).
+
+One can combine this command with projections and means, so that for example, to compare co2 depletion at the south pole:
+
+gcm.py -f POLAR_NIGHT_RUN/start16.nc --var co2 --vert 0,150 --proj spstere --time 1 -i 4 --fref POLAR_NIGHT_REF/start16.nc --operation - --mope -0.5 --Mope 0.5 --title "Polar co2 with new parametrizations" --titleref "Reference run"
+
+***********************************************************************************
+2D plot of data with missing values, along a pressure axis (decreasing with height).
+***********************************************************************************
+
+Goal:
+
+By default, python will force the y-axis of a 2D plot to be ordered by increasing values. Here is how to force it otherwise.
+
+Command:
+
+gcm.py -f TES.MappedClimatology.nadir.MY25.nc --var T_nadir_day --lat -76. --time 90 --ymin 500 --ymax 1 -m 128 -M 148 -H
+
+Note:
+
+The axis reversal is done by specifying ymin and ymax in the right order. One can also simply use --inverty and not specify (ymin,ymax). Missing values (out of range values) are replaced by holes by the option "-H".
Index: /trunk/UTIL/PYTHON/api_wrapper.py
===================================================================
--- /trunk/UTIL/PYTHON/api_wrapper.py	(revision 391)
+++ /trunk/UTIL/PYTHON/api_wrapper.py	(revision 392)
@@ -37,5 +37,5 @@
         onelevel = -99999.
 
-    #print input_name, output_name
+    print input_name, output_name
 
     if nocall:     pass
Index: unk/UTIL/PYTHON/gcm.py
===================================================================
--- /trunk/UTIL/PYTHON/gcm.py	(revision 391)
+++ 	(revision )
@@ -1,188 +1,0 @@
-#!/usr/bin/env python
-
-### T. Navarro + A. Spiga + A. Colaitis
-
-###########################################################################################
-###########################################################################################
-### 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 zrecast_wrapper import call_zrecast
-    from netCDF4 import Dataset
-    from myplot import getlschar, separatenames, readslices, adjust_length
-    from os import system
-    from planetoplot import planetoplot
-    from myscript import getparseroptions
-    import numpy as np
-
-    #############################
-    ### Get options and variables
-    parser = OptionParser()
-    getparseroptions(parser)
-    (opt,args) = parser.parse_args()
-
-    if opt.file 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)
-
-
-    ################################
-    ### General check
- 
-    if opt.fref is not None:
-       if opt.operat is None:
-          print "you must specify an operation when using a reference file"
-          exit()
-    if opt.operat in ["+","-"]:
-       if opt.fref is None:
-          print "you must specifiy a reference file when using inter-file operations"
-          exit()
-
-    interpref=False
-    if opt.fref is not None:
-       if opt.operat is not None:
-          if opt.itp is not None:
-             interpref=True
-
-    ################################
-
-
-    print "file, length", opt.file, len(opt.file)
-
-    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.file)):
-      for j in range(len(opt.var)):
-
-        zefiles = separatenames(opt.file[i])
-        print "zefiles", zefiles
-        
-        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])
-        zevars = zevar[j]
-        zevar = zevar[0]
-        print "var, zevar", opt.var, zevar
-        
-        #checkcoherence(len(zefiles),len(opt.slat),len(opt.slon),len(opt.stime))
-        
-        zefile = zefiles[0]
-              
-        #zelevel = opt.lvl   
-        stralt = None
-        #[lschar,zehour,zehourin] = getlschar ( zefile )  ## getlschar from wrfout (or simply return "" if another file)
-        [lschar,zehour,zehourin] = ["",0,0]  ## dummy
-    
-        inputnvert = separatenames(opt.lvl)
-        if np.array(inputnvert).size == 1:
-            zelevel = float(inputnvert[0])
-            ze_interp_levels = [-9999.]
-        else:
-            zelevel = -99.
-            ze_interp_levels = np.linspace(float(inputnvert[0]),float(inputnvert[1]),float(inputnvert[2]))
-        print 'level: ', zelevel
-        print 'interp_levels: ',ze_interp_levels
-
-        #####################################################
-        ### Call Fortran routines for vertical interpolations  -->  zrecast for GCM ?       
-#        if opt.itp is not None:
-#            if zelevel == 0. and opt.itp == 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.itp, \
-#                                     interp_level    = ze_interp_levels, \
-#                                     onelevel        = zelevel, \
-#                                     nocall          = opt.nocall )
-#            print zefile
-#            zelevel = 0 ## so that zelevel could play again the role of nvert
-
-# A.C.  #####################################################
-        ### Call Fortran routines for vertical interpolations  -->  zrecast
-
-        if opt.itp is not None:
-           interpolated_files=""
-           interpolated_files=call_zrecast(interp_mode=opt.itp,\
-                    input_name=zefiles,\
-                    fields=zevars)
-
-           zefiles=interpolated_files
-        if interpref:
-           interpolated_ref=""
-           interpolated_ref=call_zrecast(interp_mode=opt.itp,\
-                    input_name=[opt.fref],\
-                    fields=zevars)
-
-           reffile=interpolated_ref[0]
-        else:
-           reffile=opt.fref
-# Divers ####################################################
-    
-        zexaxis=[opt.xmin,opt.xmax]
-        zeyaxis=[opt.ymin,opt.ymax]
-
-        #############
-        ### Main call
-        name = planetoplot (zefiles,level=int(zelevel),vertmode=opt.itp,\
-                proj=opt.proj,back=opt.back,target=opt.tgt,stride=opt.ste,var=zevar,\
-                numplot=opt.num,colorb=opt.clb,winds=opt.winds,\
-                addchar=lschar,interv=[zehour,zehourin],vmin=zevmin,vmax=zevmax,\
-                tile=opt.tile,zoom=opt.zoom,display=opt.display,\
-                itstep=opt.it,hole=opt.hole,save=opt.save,\
-                anomaly=opt.anomaly,var2=opt.var2,ndiv=opt.ndiv,first=opt.frt,\
-                mult=opt.mult,zetitle=opt.zetitle,\
-                slon=zeslon,slat=zeslat,svert=zesvert,stime=zestime,\
-                outputname=opt.out,resolution=opt.res,\
-                ope=opt.operat,fileref=reffile,minop=opt.vminope,maxop=opt.vmaxope,titleref=opt.titref,\
-                invert_y=opt.inverty,xaxis=zexaxis,yaxis=zeyaxis,ylog=opt.logy,yintegral=opt.column,\
-                blat=opt.blat,tsat=opt.tsat)
-        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: unk/UTIL/PYTHON/meso.py
===================================================================
--- /trunk/UTIL/PYTHON/meso.py	(revision 391)
+++ 	(revision )
@@ -1,184 +1,0 @@
-#!/usr/bin/env python
-
-### A. Spiga
-
-###########################################################################################
-###########################################################################################
-### 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, adjust_length
-    from os import system
-    from planetoplot import planetoplot
-    from myscript import getparseroptions
-    import numpy as np
-
-    #############################
-    ### Get options and variables
-    parser = OptionParser()
-    getparseroptions(parser)
-    (opt,args) = parser.parse_args()
-
-    if opt.file 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)
-
-
-    ################################
-    ### General check
- 
-    if opt.fref is not None:
-       if opt.operat is None:
-          print "you must specify an operation when using a reference file"
-          exit()
-    if opt.operat in ["+","-"]:
-       if opt.fref is None:
-          print "you must specifiy a reference file when using inter-file operations"
-          exit()
-
-    interpref=False
-    if opt.fref is not None:
-       if opt.operat is not None:
-          if opt.itp is not None:
-             interpref=True
-
-    ################################
-
-
-    print "file, length", opt.file, len(opt.file)
-
-    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.file)):
-
-        zefiles = separatenames(opt.file[i])
-        print "zefiles", zefiles
-
-        #zefile = opt.file[i]
-        #print zefile    
-        zefile = zefiles[0]
-              
-        #zelevel = opt.lvl   
-        stralt = None
-        [lschar,zehour,zehourin] = getlschar ( zefile )  ## getlschar from wrfout (or simply return "" if another file)
-    
-        inputnvert = separatenames(opt.lvl)
-        if np.array(inputnvert).size == 1:
-            zelevel = float(inputnvert[0])
-            ze_interp_levels = [-9999.]
-        else:
-            zelevel = -99.
-            ze_interp_levels = np.linspace(float(inputnvert[0]),float(inputnvert[1]),float(inputnvert[2]))
-        print 'level: ', zelevel
-        print 'interp_levels: ',ze_interp_levels
-
-        #####################################################
-        ### Call Fortran routines for vertical interpolations        
-        if opt.itp is not None:
-            if zelevel == 0. and opt.itp == 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.itp, \
-                                     interp_level    = ze_interp_levels, \
-                                     onelevel        = zelevel, \
-                                     nocall          = opt.nocall )
-            print zefile
-            zelevel = 0 ## so that zelevel could play again the role of nvert
-
-        if opt.var is None: zerange = [-999999]
-        else:               zerange = range(zelen) 
-
-
-
-
-
-
-
-
-        if interpref:
-           interpolated_ref=""
-           interpolated_ref=call_zrecast(interp_mode=opt.itp,\
-                    input_name=[opt.fref],\
-                    fields=zevars)
-
-           reffile=interpolated_ref[0]
-        else:
-           reffile=opt.fref
-# Divers ####################################################
-    
-        zexaxis=[opt.xmin,opt.xmax]
-        zeyaxis=[opt.ymin,opt.ymax]
-
-        for jjj in zerange:
-            if jjj == -999999: 
-                argvar  = None
-                zevmin = None
-                zevmax = None
-            else:
-                argvar = opt.var[jjj]
-                if vmintab[jjj] != -999999:  zevmin = vmintab[jjj]
-                else:                        zevmin = None
-                if vmaxtab[jjj] != -999999:  zevmax = vmaxtab[jjj] 
-                else:                        zevmax = None
-
-            #############
-            ### Main call
-            name = planetoplot (zefiles,level=int(zelevel),vertmode=opt.itp,\
-                proj=opt.proj,back=opt.back,target=opt.tgt,stride=opt.ste,var=argvar,\
-                numplot=opt.num,colorb=opt.clb,winds=opt.winds,\
-                addchar=lschar,interv=[zehour,zehourin],vmin=zevmin,vmax=zevmax,\
-                tile=opt.tile,zoom=opt.zoom,display=opt.display,\
-                itstep=opt.it,hole=opt.hole,save=opt.save,\
-                anomaly=opt.anomaly,var2=opt.var2,ndiv=opt.ndiv,first=opt.frt,\
-                mult=opt.mult,zetitle=opt.zetitle,\
-                slon=zeslon,slat=zeslat,svert=zesvert,stime=zestime,\
-                outputname=opt.out,resolution=opt.res,\
-                ope=opt.operat,fileref=reffile,minop=opt.vminope,maxop=opt.vmaxope,titleref=opt.titref,\
-                invert_y=opt.inverty,xaxis=zexaxis,yaxis=zeyaxis,ylog=opt.logy,yintegral=opt.column,\
-                blat=opt.blat,tsat=opt.tstat)
-        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 + ' '
-    f = open(name+'.sh', 'w')
-    f.write(command)
Index: /trunk/UTIL/PYTHON/myplot.py
===================================================================
--- /trunk/UTIL/PYTHON/myplot.py	(revision 391)
+++ /trunk/UTIL/PYTHON/myplot.py	(revision 392)
@@ -43,5 +43,5 @@
     elif 'HGT_M' in nc.variables:    typefile = 'geo'
     #else:                            errormess("whatkindfile: typefile not supported.")
-    else: typefile = 'gcm' # for lslin-ed files from gcm
+    else:                            typefile = 'gcm' # for lslin-ed files from gcm
     return typefile
 
@@ -50,5 +50,5 @@
     ## 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 ..."
+    #print "   Opening variable with", dimension, "dimensions ..."
     if dimension == 2:    field = nc.variables[var][:,:]
     elif dimension == 3:  field = nc.variables[var][:,:,:]
@@ -65,5 +65,5 @@
     shape = np.array(input).shape
     #print 'd1,d2,d3,d4: ',d1,d2,d3,d4
-    print 'dim,shape: ',dimension,shape
+    print 'IN  REDUCEFIELD dim,shape: ',dimension,shape
     output = input
     error = False
@@ -98,39 +98,61 @@
         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 = mean(input[d4,:,:,:],axis=0); output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt[d3]); output = mean(output[d2,:],axis=0); output = mean(output[d1],axis=0)
+             output = mean(input[d4,:,:,:],axis=0)
+             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
+             output = mean(output[d2,:],axis=0)
+             output = mean(output[d1],axis=0)
         elif d4 is not None and d3 is not None and d2 is not None: 
-            output = mean(input[d4,:,:,:],axis=0); output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt[d3]); output = mean(output[d2,:],axis=0)
+             output = mean(input[d4,:,:,:],axis=0)
+             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
+             output = mean(output[d2,:],axis=0)
         elif d4 is not None and d3 is not None and d1 is not None: 
-            output = mean(input[d4,:,:,:],axis=0); output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt[d3]); output = mean(output[:,d1],axis=1)
+             output = mean(input[d4,:,:,:],axis=0)
+             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
+             output = mean(output[:,d1],axis=1)
         elif d4 is not None and d2 is not None and d1 is not None: 
-            output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)
+             output = mean(input[d4,:,:,:],axis=0)
+             output = mean(output[:,d2,:],axis=1)
+             output = mean(output[:,d1],axis=1)
         elif d3 is not None and d2 is not None and d1 is not None: 
-            output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt[d3]); output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1) 
-        elif d4 is not None and d3 is not None:  output = mean(input[d4,:,:,:],axis=0); output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt[d3])
-        elif d4 is not None and d2 is not None:  output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,d2,:],axis=1)
-        elif d4 is not None and d1 is not None:  output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,:,d1],axis=2)
-        elif d3 is not None and d2 is not None:  output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt[d3]); output = mean(output[:,d2,:],axis=1)
-        elif d3 is not None and d1 is not None:  output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt[d3]); output = mean(output[:,:,d1],axis=0)
-        elif d2 is not None and d1 is not None:  output = mean(input[:,:,d2,:],axis=2); output = mean(output[:,:,d1],axis=2)
-        elif d1 is not None:                     output = mean(input[:,:,:,d1],axis=3)
-        elif d2 is not None:                     output = mean(input[:,:,d2,:],axis=2)
-        elif d3 is not None:                     output = reduce_zaxis(input[:,d3,:,:],ax=0,yint=yint,vert=alt[d3])
-        elif d4 is not None:                     output = mean(input[d4,:,:,:],axis=0)
+             output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3) 
+             output = mean(output[:,d2,:],axis=1)
+             output = mean(output[:,d1],axis=1) 
+        elif d4 is not None and d3 is not None:  
+             output = mean(input[d4,:,:,:],axis=0) 
+             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
+        elif d4 is not None and d2 is not None:  
+             output = mean(input[d4,:,:,:],axis=0) 
+             output = mean(output[:,d2,:],axis=1)
+        elif d4 is not None and d1 is not None:  
+             output = mean(input[d4,:,:,:],axis=0) 
+             output = mean(output[:,:,d1],axis=2)
+        elif d3 is not None and d2 is not None:  
+             output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
+             output = mean(output[:,d2,:],axis=1)
+        elif d3 is not None and d1 is not None:  
+             output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
+             output = mean(output[:,:,d1],axis=0)
+        elif d2 is not None and d1 is not None:  
+             output = mean(input[:,:,d2,:],axis=2)
+             output = mean(output[:,:,d1],axis=2)
+        elif d1 is not None:        output = mean(input[:,:,:,d1],axis=3)
+        elif d2 is not None:        output = mean(input[:,:,d2,:],axis=2)
+        elif d3 is not None:        output = reduce_zaxis(input[:,d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
+        elif d4 is not None:        output = mean(input[d4,:,:,:],axis=0)
     dimension = np.array(output).ndim
     shape = np.array(output).shape
-    print 'dim,shape: ',dimension,shape
+    print 'OUT REDUCEFIELD dim,shape: ',dimension,shape
     return output, error
 
-## Author: AC
-
-def reduce_zaxis (input,ax=None,yint=False,vert=None):
+## Author: AC + AS
+def reduce_zaxis (input,ax=None,yint=False,vert=None,indice=None):
     from mymath import max,mean
     from scipy import integrate
-    if yint:
+    if yint and vert is not None and indice is not None:
        if type(input).__name__=='MaskedArray':
           input.set_fill_value([np.NaN])
-          output = integrate.trapz(input.filled(),x=vert,axis=ax)
+          output = integrate.trapz(input.filled(),x=vert[indice],axis=ax)
        else:
-          output = integrate.trapz(input.filled(),x=vert,axis=ax)
+          output = integrate.trapz(input.filled(),x=vert[indice],axis=ax)
     else:
        output = mean(input,axis=ax)
@@ -257,13 +279,13 @@
             if cen_lat > 10.:    
                 proj="npstere"
-                print "NP stereographic polar domain" 
+                #print "NP stereographic polar domain" 
             else:            
                 proj="spstere"
-                print "SP stereographic polar domain"
+                #print "SP stereographic polar domain"
         elif map_proj == 1: 
-            print "lambert projection domain" 
+            #print "lambert projection domain" 
             proj="lcc"
         elif map_proj == 3: 
-            print "mercator projection"
+            #print "mercator projection"
             proj="merc"
         else:
@@ -460,5 +482,5 @@
         elif wlat[1] >= 0.:    blat = wlat[0] 
         elif wlat[0] <= 0.:    blat = wlat[1]
-    print "blat ", blat
+    #print "blat ", blat
     h = 50.  ## en km
     radius = 3397200.
@@ -484,5 +506,4 @@
     #    step = np.min([5.,step])
     #    steplon = step
-    print step
     m.drawmeridians(np.r_[-180.:180.:steplon], labels=[0,0,0,1], color='grey', fontsize=fontsizemer)
     m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color='grey', fontsize=fontsizemer)
@@ -534,6 +555,6 @@
     ###
     if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0.
-    print "field ", min(fieldcalc), max(fieldcalc)
-    print "bounds ", zevmin, zevmax
+    print "BOUNDS field ", min(fieldcalc), max(fieldcalc)
+    print "BOUNDS adopted ", zevmin, zevmax
     return zevmin, zevmax
 
@@ -545,9 +566,8 @@
     if zevmin < 0: what_I_plot[ what_I_plot < zevmin*(1. - 1.e-7) ] = zevmin*(1. - 1.e-7)
     else:          what_I_plot[ what_I_plot < zevmin*(1. + 1.e-7) ] = zevmin*(1. + 1.e-7)
-    print "new min ", min(what_I_plot)
+    print "NEW MIN ", min(what_I_plot)
     what_I_plot[ what_I_plot > 9e+35  ] = -9e+35
     what_I_plot[ what_I_plot > zevmax ] = zevmax
-    print "new max ", max(what_I_plot)
-    
+    print "NEW MAX ", max(what_I_plot)
     return what_I_plot
 
@@ -556,5 +576,5 @@
     from mymath import max,min
     lim = 0.15*0.5*(abs(max(what_I_plot))+abs(min(what_I_plot)))
-    print "on vire en dessous de ", lim
+    print "NO PLOT BELOW VALUE ", lim
     what_I_plot [ abs(what_I_plot) < lim ] = 1.e40 
     return what_I_plot
@@ -566,5 +586,5 @@
     [wlon,wlat] = [ [wlon[0]+zoom*dlon/100.,wlon[1]-zoom*dlon/100.],\
                     [wlat[0]+zoom*dlat/100.,wlat[1]-zoom*dlat/100.] ]
-    print "zoom %",zoom,wlon,wlat
+    print "ZOOM %",zoom,wlon,wlat
     return wlon,wlat
 
@@ -642,5 +662,5 @@
 #W --> spectral ou jet
 #spectral BrBG RdBu_r
-    print "predefined colorbars"
+    #print "predefined colorbars"
     if whichone not in whichcolorb:
         whichone = "def"
@@ -813,5 +833,5 @@
       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."
+         print "INFO: whole longitude averaging asked, so last point is not taken into account."
   return zeindex
      
@@ -828,19 +848,19 @@
    shape = what_I_plot.shape
    if indextime is None:
-      print "axis is time"
+      print "AXIS is time"
       x = time
       count = count+1
    if indexlon is None:
-      print "axis is lon"
+      print "AXIS is lon"
       if count == 0: x = lon
       else: y = lon
       count = count+1
    if indexlat is None:
-      print "axis is lat"
+      print "AXIS is lat"
       if count == 0: x = lat
       else: y = lat
       count = count+1
    if indexvert is None and dim0 is 4:
-      print "axis is vert"
+      print "AXIS is vert"
       if vertmode == 0: # vertical axis is as is (GCM grid)
          if count == 0: x=range(len(vert))
@@ -853,6 +873,6 @@
    x = array(x)
    y = array(y)
-   print "what_I_plot.shape", what_I_plot.shape
-   print "x.shape, y.shape", x.shape, y.shape
+   print "CHECK: what_I_plot.shape", what_I_plot.shape
+   print "CHECK: x.shape, y.shape", x.shape, y.shape
    if len(shape) == 1:
       print shape[0]
@@ -864,5 +884,5 @@
       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
+         print "INFO: swapaxes", what_I_plot.shape, shape
          #x0 = x
          #x = y
@@ -896,8 +916,4 @@
     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
 
     return nlon, nlat, nvert, ntime, mapmode, nslices
Index: /trunk/UTIL/PYTHON/myscript.py
===================================================================
--- /trunk/UTIL/PYTHON/myscript.py	(revision 391)
+++ /trunk/UTIL/PYTHON/myscript.py	(revision 392)
@@ -34,5 +34,5 @@
     parser.add_option('-s', '--stride', action='store',dest='ste',       type="int",     default=3,     help='stride vectors [3]')
     parser.add_option('-x', '--no-vect',action='store_false',dest='winds',               default=True,  help='no wind vectors')
-    parser.add_option('-n', '--num',    action='store',dest='num',       type="int",     default=None,  help='plot number (<0: plot LT -*numplot*) [2]')
+    parser.add_option('-n', '--num',    action='store',dest='num',       type="int",     default=None,  help='plot number (<0: plot LT -*numplot*) [1]')
     parser.add_option('-z', '--zoom',   action='store',dest='zoom',      type="float",   default=None,  help='zoom factor in %')
     parser.add_option('-e', '--itstep', action='store',dest='it',        type="int",     default=None,  help='stride time [4]')
@@ -44,5 +44,5 @@
     parser.add_option('--lon',          action='append',dest='slon',   type="string",  default=None, help='slices along lon. 2 comma-separated values: averaging')
     parser.add_option('--vert',         action='append',dest='svert',  type="string",  default=None, help='slices along vert. 2 comma-separated values: averaging') 
-    parser.add_option('--column',         action='store_true',dest='column',               default=False,help='changes the function of --vert z1,z2 from MEAN to INTEGRATE along z')
+    parser.add_option('--column',         action='store_true',dest='column',           default=False,help='changes the function of --vert z1,z2 from MEAN to INTEGRATE along z')
     parser.add_option('--time',         action='append',dest='stime',  type="string",  default=None, help='slices along time. 2 comma-separated values: averaging') 
     parser.add_option('--xmax',         action='store',dest='xmax',    type="float",   default=None, help='max value for x-axis in contour-plots [max(xaxis)]')
Index: /trunk/UTIL/PYTHON/planetoplot.py
===================================================================
--- /trunk/UTIL/PYTHON/planetoplot.py	(revision 391)
+++ /trunk/UTIL/PYTHON/planetoplot.py	(revision 392)
@@ -5,5 +5,6 @@
 ### A. Spiga     -- LMD -- 06~09/2011 -- General building and mapping capabilities
 ### T. Navarro   -- LMD -- 10~11/2011 -- Improved use for GCM and added sections + 1Dplot capabilities 
-### A. Colaitis  -- LMD --            -- Mostly minor improvements and inter-plot operation capabilities + zrecast interpolation for gcm
+### A. Colaitis  -- LMD --    11/2011 -- Mostly minor improvements and inter-plot operation capabilities + zrecast interpolation for gcm
+### A. Spiga     -- LMD --    11/2011 -- Extended multivar subplot capabilities + cosmetic changes + general cleaning and tests
 
 def planetoplot (namefiles,\
@@ -63,8 +64,8 @@
                        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,determineplot
+                       getname,localtime,polarinterv,getsindex,define_axis,determineplot,readslices
     from mymath import deg,max,min,mean,get_tsat
     import matplotlib as mpl
-    from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel
+    from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel, title
     from matplotlib.cm import get_cmap
     import numpy as np
@@ -72,26 +73,40 @@
 
     ################################
+    ### Preliminary stuff
+    ################################
+    print "********************************************"
+    print "********** WELCOME TO PLANETOPLOT **********"
+    print "********************************************"
+    if not isinstance(namefiles, np.ndarray): namefiles = [namefiles]
+    if not isinstance(var, np.ndarray):       var = [var]
+    if ope is not None and len(var) > 1:      errormess("you are using an operation... please set only one var !")
+
+    ################################
     ### Which plot needs to be done?
+    ################################
     nlon, nlat, nvert, ntime, mapmode, nslices = determineplot(slon, slat, svert, stime)
-    if mapmode == 0: winds=False
-    if not isinstance(namefiles, np.ndarray): namefiles = [namefiles]
-    zelen = len(namefiles)
-    if numplot == None:  numplot = zelen*nslices
-    print "len(namefiles), nslices, numplot: ", zelen, nslices, numplot
-    if fileref is not None:
-        all_var  = [[]]*(zelen+2)
-    else:
-        all_var   = [[]]*zelen
-        all_var2  = [[]]*zelen
-        all_title = [[]]*zelen
-    
-    #########################
-    ### Loop over the files initially separated by comas to be plotted on the same figure
+    if mapmode == 0:       winds=False
+    elif mapmode == 1:     
+        if svert is None:  svert = readslices(str(level)) ; nvert=1
+    zelen = len(namefiles)*len(var)
+    if numplot is None:  numplot = zelen*nslices
+    print "********** FILES, SLICES, VARS, TOTAL PLOTS: ", len(namefiles), nslices, len(var), numplot
+    print "********** MAPMODE: ", mapmode
+    if fileref is not None:       all_var  = [[]]*(zelen+2) ;  all_varname = [[]]*(zelen+2)
+    else:                         all_var  = [[]]*zelen     ;  all_var2  = [[]]*zelen ;  all_title = [[]]*zelen ;  all_varname = [[]]*zelen
+ 
+    #################################################################################################
+    ### Loop over the files + vars initially separated by commas to be plotted on the same figure ###
+    #################################################################################################
     k = 0
     firstfile = True
-    for namefile in namefiles:
-      print namefile 
+    for nnn in range(len(namefiles)):
+     for vvv in range(len(var)): 
+
+      print "********** LOOP..... THIS IS SUBPLOT NUMBER.....",k
+
       ######################
       ### Load NETCDF object
+      namefile = namefiles[nnn] ; print "********** THE NAMEFILE IS....", namefile 
       nc  = Dataset(namefile)
 
@@ -100,20 +115,19 @@
       ### ... TYPEFILE
       typefile = whatkindfile(nc)                                  
-      if firstfile:
-         typefile0 = typefile
-      elif typefile != typefile0:
-         errormess("Not the same kind of files !", [typefile0, typefile])
+      if firstfile:                 typefile0 = typefile
+      elif typefile != typefile0:   errormess("Not the same kind of files !", [typefile0, typefile])
       ### ... VAR
-      varname=var
-      if var not in nc.variables: var = False
+      varname=var[vvv]
+      print "********** THE VAR IS....",varname, var2
+      if varname not in nc.variables: varname = False
       ### ... WINDS
       if 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)
+      if not varname and not winds: errormess("please set at least winds or var",printvar=nc.variables)
       ### ... COORDINATES, could be moved below
       [lon2d,lat2d] = getcoorddef(nc)                       
       ### ... PROJECTION
-      if proj == None:   proj = getproj(nc)                  
+      if proj == None:   proj = getproj(nc)                 
 
 ##########################################################
@@ -131,16 +145,18 @@
       elif typefile in ['meso','mesoapi']:
           if mapmode == 0:
-              if var in ['PHTOT','W']:    vertdim='BOTTOM-TOP_PATCH_END_STAG'
-              else:                       vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
-              if var in ['V']:  latdim='SOUTH-NORTH_PATCH_END_STAG'
-              else:             latdim='SOUTH-NORTH_PATCH_END_UNSTAG'
-              if var in ['U']:  londim='WEST-EAST_PATCH_END_STAG'
-              else:             londim='WEST-EAST_PATCH_END_UNSTAG'
+              if varname in ['PHTOT','W']:    vertdim='BOTTOM-TOP_PATCH_END_STAG'
+              else:                           vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
+              if varname in ['V']:  latdim='SOUTH-NORTH_PATCH_END_STAG'
+              else:                 latdim='SOUTH-NORTH_PATCH_END_UNSTAG'
+              if varname in ['U']:  londim='WEST-EAST_PATCH_END_STAG'
+              else:                 londim='WEST-EAST_PATCH_END_UNSTAG'
               lon = np.arange(0,getattr(nc,londim),1)
               lat = np.arange(0,getattr(nc,latdim),1)
-              if vertmode is None: vertmode=0
-              if vertmode == 0:   vert = np.arange(0,getattr(nc,vertdim),1)
-              else:               vert = nc.variables["vert"][:]
+              if vertmode is None:  vertmode=0
+              if vertmode == 0:     vert = np.arange(0,getattr(nc,vertdim),1)
+              else:                 vert = nc.variables["vert"][:]
               time = np.arange(0,len(nc.variables["Times"]),1)
+          else: 
+              lon=None ; lat=None ; vert=None ; time=None
        #if firstfile:
        #   lat0 = lat
@@ -160,26 +176,21 @@
          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)
+         ### ... NAME FOR TITLE and FILE (have to do something for multiple vars...)
+         basename = getname(var=varname,winds=winds,anomaly=anomaly)
          basename = basename + getstralt(nc,level)  ## can be moved elsewhere for a more generic routine
 
-      print "var, var2: ", var, var2
-      if varname in ["temp","t","T_nadir_nit","T_nadir_day"] and var and tsat:
-          print "computing Tsat-T, I ASSUME Z-AXIS IS PRESSURE"
-          tt=getfield(nc,var)
-          if type(tt).__name__=='MaskedArray':
-             tt.set_fill_value([np.NaN])
-             tinput=tt.filled()
-          else:
-             tinput=tt
+      ##### SPECIFIC
+      if varname in ["temp","t","T_nadir_nit","T_nadir_day"] and tsat:
+          tt=getfield(nc,varname) ; print "computing Tsat-T, I ASSUME Z-AXIS IS PRESSURE"
+          if type(tt).__name__=='MaskedArray':  tt.set_fill_value([np.NaN]) ; tinput=tt.filled()
+          else:                                 tinput=tt
           all_var[k]=get_tsat(vert,tinput,zlon=lon,zlat=lat,zalt=vert,ztime=time)
       else:
-         if var: all_var[k] = getfield(nc,var)
+      ##### GENERAL STUFF HERE
+          all_var[k] = getfield(nc,varname)
+          all_varname[k] = varname
       if var2: all_var2[k] = getfield(nc,var2)
       
-      print "k", k
-      print "all_var[k].shape", all_var[k].shape
+      print "********** all_var[k].shape", all_var[k].shape
       k += 1
       firstfile = False
@@ -188,13 +199,13 @@
     ##################################
     ### Operation on files
+    ### ... for the moment only OK when 1 var only
     if ope is not None:
-       if var not in nc.variables: var = False
-       if var: 
              print ope
              if ope in ["-","+"]:
-                if fileref is not None:   all_var[k] = getfield(Dataset(fileref),var)
+                if fileref is not None:   all_var[k] = getfield(Dataset(fileref),all_varname[k-1]) ; all_varname[k] = all_varname[k-1]
                 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]
+                all_varname[k+1] = all_varname[k]
                 numplot = numplot+2
              elif ope in ["cat"]:
@@ -219,8 +230,7 @@
     elif numplot <= 0:                 itstep = 1 
     
-    #for nplot in range(numplot):
+    print "********************************************"
     while error is False:
-       print "nplot", nplot
-       print error
+       print "********** nplot", nplot, "itime",itime,"error",error
        
        ### Which local time ?
@@ -249,31 +259,28 @@
 
 ####################################################################
-       if typefile in ['meso','mesoapi'] and mapmode == 1: 
-               indextime = itime
-               indexlon = None
-               indexlat = None
-               indexvert = level  ## ou svert ??? 
-               nlon = 1
-               nlat = 1
-               nvert = 1
-               ntime = 1 
+       ## 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
+       indexlon  = getsindex(slon,(nplot-1)%nlon,lon)
+       indexlat  = getsindex(slat,((nplot-1)//nlon)%nlat,lat)
+       indexvert = getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert) 
+       if mapmode == 1 and stime is None: 
+           indextime = itime
+           if len(var) != 1 or len(namefiles) != 1: indextime = first
        else:
-           ## get all indexes to be taken into account for this subplot and then reduce field
-           ## We plot 1) all lon slices 2) all lat slices 3) all vert slices 4) all time slices and then go to the next slice
-           indexlon  = getsindex(slon,(nplot-1)%nlon,lon)
-           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)  
 
        if fileref is not None:
-           index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%(len(namefiles)+2)
+           index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%(len(namefiles)+2)  ## OK only 1 var, see test in the beginning
        else:
-           index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%len(namefiles)
-
-       print nlon, nlat, nvert, ntime  
-       print slon, slat, svert, stime       
-       print index_f
-       print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime)
+           yeah = len(namefiles)*len(var)
+           index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%yeah
+       #print nlon, nlat, nvert, ntime  
+       #print slon, slat, svert, stime       
+       #print index_f
+       #print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime)
+       #if varname != all_varname[index_f]:  indextime = first
 ####################################################################
+
+       ticks = ndiv + 1
 
        #### Contour plot
@@ -285,5 +292,5 @@
                   if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot,6)
               zevmin, zevmax = calculate_bounds(what_I_plot)
-              zelevels = np.linspace(zevmin,zevmax,num=ndiv+1) #20)
+              zelevels = np.linspace(zevmin,zevmax,ticks) #20)
               if var2 == 'HGT':  zelevels = np.arange(-10000.,30000.,2000.)
               if mapmode == 0:
@@ -293,5 +300,5 @@
               if len(what_I_plot.shape) is 2:
                   cs = contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #0.33 colors='w' )# , alpha=0.5)
-                  if typefile in ['gcm']: clabel(cs,fmt = '%1.2e')
+                  #if typefile in ['gcm']: clabel(cs,fmt = '%1.2e')
               ### If we plot a 1-D field
               elif len(what_I_plot.shape) is 1:
@@ -301,9 +308,10 @@
 
        #### Shaded plot
-       if var:
+       varname = all_varname[index_f]
+       if varname:
            what_I_plot, error = reducefield(all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert , yint=yintegral, alt=vert)
            what_I_plot = what_I_plot*mult
            if not error: 
-               fvar = var
+               fvar = varname
                ###
                if anomaly:
@@ -329,5 +337,5 @@
                      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)
-                     zelevels = np.linspace(zevmin,zevmax,num=2.*ndiv+1)
+                     zelevels = np.linspace(zevmin,zevmax,num=ticks)
                      #contourf(what_I_plot, zelevels, cmap = palette )
 
@@ -336,32 +344,25 @@
                      elif mapmode == 0:
                        contourf( x, y, what_I_plot, zelevels, cmap = palette)
+
                      zxmin, zxmax = xaxis
                      zymin, zymax = yaxis
-                     if zxmin is not None:
-                        mpl.pyplot.xlim(xmin=zxmin)
-                     if zxmax is not None:
-                        mpl.pyplot.xlim(xmax=zxmax)
-                     if zymin is not None:
-                        mpl.pyplot.ylim(ymin=zymin)
-                     if zymax is not None:
-                        mpl.pyplot.ylim(ymax=zymax)
-
-                     if invert_y:
-                        lima,limb = mpl.pyplot.ylim()
-                        mpl.pyplot.ylim(limb,lima)
-                     if ylog:
-                        mpl.pyplot.semilogy()
+                     if zxmin is not None: mpl.pyplot.xlim(xmin=zxmin)
+                     if zxmax is not None: mpl.pyplot.xlim(xmax=zxmax)
+                     if zymin is not None: mpl.pyplot.ylim(ymin=zymin)
+                     if zymax is not None: mpl.pyplot.ylim(ymax=zymax)
+                     if invert_y:     lima,limb = mpl.pyplot.ylim() ; mpl.pyplot.ylim(limb,lima)
+                     if ylog:         mpl.pyplot.semilogy()
 
                  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' :        
+                 if colorb != 'nobar' and varname != 'HGT' :        
                      if (fileref is not None) and (index_f is numplot-1):
                         colorbar(fraction=0.05,pad=0.03,format="%.2f",\
-                                           ticks=np.linspace(zevmin,zevmax,min([ndiv+1,20])),\
+                                           ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,20])),\
                                            extend='neither',spacing='proportional')
                      else:
                         colorbar(fraction=0.05,pad=0.03,format=fmtvar(fvar.upper()),\
-                                           ticks=np.linspace(zevmin,zevmax,min([ndiv+1,10])),\
+                                           ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,20])),\
                                            extend='neither',spacing='proportional')
 
@@ -391,5 +392,5 @@
                    key = False
                if metwind:  [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d)
-               if var == False:       colorvec = definecolorvec(back)
+               if varname == False:       colorvec = definecolorvec(back)
                else:                  colorvec = definecolorvec(colorb)
                vectorfield(vecx, vecy,\
@@ -400,4 +401,5 @@
     
        ### Next subplot
+       basename = getname(var=varname,winds=winds,anomaly=anomaly)
        if typefile in ['mesoapi','meso']:
             plottitle = basename
@@ -411,7 +413,7 @@
                    plottitle = basename+' '+fileref
                 else:
-                   plottitle = basename+' '+namefiles[index_f]
+                   plottitle = basename+' '+namefiles[0]#index_f]
             else:
-                plottitle = basename+' '+namefiles[index_f]
+                plottitle = basename+' '+namefiles[0]#index_f]
        if mult != 1:           plottitle = str(mult) + "*" + plottitle
        if zetitle != "fill":   
@@ -432,5 +434,5 @@
 #       if indextime is not None:
 #         plottitle = plottitle + " time: " + str(min(time[indextime])) +" "+ str(max(time[indextime]))
-       ptitle( plottitle )
+       title( plottitle )
        itime += itstep
        if nplot >= numplot:
@@ -464,5 +466,5 @@
     if found_lct:     
         pad_inches_value = 0.35
-        print "save", save
+        print "********** SAVE ", save
         if save == 'png': 
             if display: makeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value) #,erase=True)  ## a miniature
@@ -473,7 +475,7 @@
             show()
         else: 
-            print "save mode not supported. using gui instead."
+            print "INFO: save mode not supported. using gui instead."
             show()
-    else:   print "Local time not found"
+    else:   print "!!! Local time not found"
 
     ###############
Index: /trunk/UTIL/PYTHON/pp.py
===================================================================
--- /trunk/UTIL/PYTHON/pp.py	(revision 392)
+++ /trunk/UTIL/PYTHON/pp.py	(revision 392)
@@ -0,0 +1,148 @@
+#!/usr/bin/env python
+
+### A. Spiga + T. Navarro + A. Colaitis
+
+###########################################################################################
+###########################################################################################
+### 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 zrecast_wrapper import call_zrecast
+    from netCDF4 import Dataset
+    from myplot import getlschar, separatenames, readslices, adjust_length, whatkindfile, errormess
+    from os import system
+    from planetoplot import planetoplot
+    from myscript import getparseroptions
+    import numpy as np
+
+    #############################
+    ### Get options and variables
+    parser = OptionParser() ; getparseroptions(parser) ; (opt,args) = parser.parse_args()
+    if opt.file is None: 				errormess("I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h")
+    if opt.var is None and opt.anomaly is True: 	errormess("Cannot ask to compute anomaly if no variable is set")
+    if opt.fref is not None and opt.operat is None: 	errormess("you must specify an operation when using a reference file")
+    if opt.operat in ["+","-"] and opt.fref is None: 	errormess("you must specifiy a reference file when using inter-file operations")
+    if opt.fref is not None and opt.operat is not None and opt.itp is not None: interpref=True
+    else:   interpref=False
+
+    #############################
+    ### Get infos about slices
+    zeslat  = readslices(opt.slat) ; zeslon  = readslices(opt.slon) ; zesvert = readslices(opt.svert) ; zestime = readslices(opt.stime)
+    reffile = opt.fref
+    zexaxis = [opt.xmin,opt.xmax] ; zeyaxis=[opt.ymin,opt.ymax]
+
+    #############################
+    ### 1. LOOP ON FILE LISTS TO BE PUT IN DIFFERENT FIGURES
+    for i in range(len(opt.file)):
+
+      zefiles = separatenames(opt.file[i])
+
+      typefile = whatkindfile(Dataset(zefiles[0]))
+      stralt = None
+      if typefile in ["meso","mesoapi"]:          
+          [lschar,zehour,zehourin] = getlschar ( zefiles[0] )
+          if opt.var is None:  opt.var = ["HGT"] 
+      else:                                       
+          [lschar,zehour,zehourin] = ["",0,0]
+          if opt.var is None:  opt.var = ["phisinit"]
+
+      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
+
+      #############################
+      ### 2. LOOP ON VAR LISTS TO BE PUT IN DIFFERENT FIGURES
+      for j in range(len(opt.var)):
+
+        zevars = separatenames(opt.var[j])
+
+        inputnvert = separatenames(opt.lvl)
+        if np.array(inputnvert).size == 1:
+            zelevel = float(inputnvert[0])
+            ze_interp_levels = [-9999.]
+        else:
+            zelevel = -99.
+            ze_interp_levels = np.linspace(float(inputnvert[0]),float(inputnvert[1]),float(inputnvert[2]))
+
+        #########################################################
+        ### Call Fortran routines for vertical interpolations ###      
+        ######################################################### 
+        if opt.itp is not None:
+          #####
+          ##### MESOSCALE : written by AS
+          #####
+          if typefile == "meso":
+            if zelevel == 0. and opt.itp == 4:  zelevel = 0.010
+            ### winds or no winds
+            if opt.winds            :  zefields = 'uvmet'
+            else                    :  zefields = ''
+            ### var or no var
+            if zefields == ''       :  zefields = opt.var[j] 
+            else                    :  zefields = zefields + "," + opt.var[j]
+            if opt.var2 is not None :  zefields = zefields + "," + opt.var2  
+            ### call fortran routines
+            for fff in range(len(zefiles)):
+                newname = api_onelevel (  path_to_input   = '', \
+                                               input_name      = zefiles[fff], \
+                                               fields          = zefields, \
+                                               interp_method   = opt.itp, \
+                                               interp_level    = ze_interp_levels, \
+                                               onelevel        = zelevel, \
+                                               nocall          = opt.nocall )
+                if fff == 0: zetab = newname 
+                else:        zetab = np.append(zetab,newname)
+            zefiles = zetab #; print zefiles
+            zelevel = 0 ## so that zelevel could play again the role of nvert
+          #####
+          ##### GCM : written by AC
+          #####
+          elif typefile == "gcm":
+            interpolated_files=""
+            interpolated_files=call_zrecast(interp_mode=opt.itp,\
+                    input_name=zefiles,\
+                    fields=zevars)
+
+            zefiles=interpolated_files
+            if interpref:
+               interpolated_ref=""
+               interpolated_ref=call_zrecast(interp_mode=opt.itp,\
+                    input_name=[opt.fref],\
+                    fields=zevars)
+
+               reffile=interpolated_ref[0]
+          else:
+            print "not supported"
+            exit()
+
+        #############
+        ### Main call
+        name = planetoplot (zefiles,level=int(zelevel),vertmode=opt.itp,\
+                proj=opt.proj,back=opt.back,target=opt.tgt,stride=opt.ste,var=zevars,\
+                numplot=opt.num,colorb=opt.clb,winds=opt.winds,\
+                addchar=lschar,interv=[zehour,zehourin],vmin=zevmin,vmax=zevmax,\
+                tile=opt.tile,zoom=opt.zoom,display=opt.display,\
+                itstep=opt.it,hole=opt.hole,save=opt.save,\
+                anomaly=opt.anomaly,var2=opt.var2,ndiv=opt.ndiv,first=opt.frt,\
+                mult=opt.mult,zetitle=opt.zetitle,\
+                slon=zeslon,slat=zeslat,svert=zesvert,stime=zestime,\
+                outputname=opt.out,resolution=opt.res,\
+                ope=opt.operat,fileref=reffile,minop=opt.vminope,maxop=opt.vmaxope,titleref=opt.titref,\
+                invert_y=opt.inverty,xaxis=zexaxis,yaxis=zeyaxis,ylog=opt.logy,yintegral=opt.column,\
+                blat=opt.blat,tsat=opt.tsat)
+        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 + ' '
+    if typefile not in ["meso","mesoapi"]: name = 'pycommand'
+    f = open(name+'.sh', 'w')
+    f.write(command)
+
+    print "********** OPTIONS: ", opt
+    print "********************************************************** END"
