Index: trunk/UTIL/PYTHON/mymath.py
===================================================================
--- trunk/UTIL/PYTHON/mymath.py	(revision 642)
+++ trunk/UTIL/PYTHON/mymath.py	(revision 647)
@@ -58,4 +58,13 @@
            else:
               return np.array(field).sum(axis=axis)
+              
+def getmask (field):
+        import numpy as np
+        if field is None: return None
+        if type(field).__name__=='MaskedArray':
+              return np.ma.getmask(field)
+        else:
+              return np.isnan(field)
+        
 
 def deg ():
@@ -66,5 +75,8 @@
     myfile = open(filename, 'w')
     for line in mydata:
-        myfile.write(str(line) + '\n')
+        zeline = str(line)
+        zeline = zeline.replace('[','')
+        zeline = zeline.replace(']','')
+        myfile.write(zeline + '\n')
     myfile.close()
     return
Index: trunk/UTIL/PYTHON/myplot.py
===================================================================
--- trunk/UTIL/PYTHON/myplot.py	(revision 642)
+++ trunk/UTIL/PYTHON/myplot.py	(revision 647)
@@ -55,4 +55,5 @@
 ## Author: AS, AC, JL
 def whatkindfile (nc):
+    typefile = 'gcm' # default
     if 'controle' in nc.variables:             typefile = 'gcm'
     elif 'phisinit' in nc.variables:           typefile = 'gcm'
@@ -62,5 +63,4 @@
     elif hasattr(nc,'institution'):
       if "European Centre" in getattr(nc,'institution'):  typefile = 'ecmwf'
-    else:                                      typefile = 'gcm' # for lslin-ed files from gcm
     return typefile
 
@@ -165,5 +165,5 @@
     ### ... note, anomaly is only computed over d1 and d2 for the moment
     import numpy as np
-    from mymath import max,mean,min,sum
+    from mymath import max,mean,min,sum,getmask
     csmooth = 12 ## a fair amount of grid points (too high results in high computation time)
     if redope is not None:
@@ -194,8 +194,12 @@
         elif max(d1) >= shape[1]: error = True
         elif d1 is not None and d2 is not None:
-	  totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0)
-	  output = output*mesharea; output = sum(output[d2,:],axis=0); output = sum(output[d1],axis=0)/totalarea
+          totalarea = np.ma.masked_where(getmask(output),mesharea)
+          totalarea = mean(totalarea[d2,:],axis=0);totalarea = mean(totalarea[d1])
+          output = output*mesharea; output = mean(output[d2,:],axis=0); output = mean(output[d1])/totalarea
         elif d1 is not None:         output = mean(input[:,d1],axis=1)
-        elif d2 is not None:         totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[d2,:],axis=0)/totalarea
+        elif d2 is not None:
+          totalarea = np.ma.masked_where(getmask(output),mesharea)
+          totalarea = mean(totalarea[d2,:],axis=0)
+          output = output*mesharea; output = mean(output[d2,:],axis=0)/totalarea
     elif dimension == 3:
         if mesharea is None: mesharea=np.ones(shape[[1,2]])
@@ -203,20 +207,29 @@
         elif max(d2) >= shape[1]: error = True
         elif max(d1) >= shape[2]: error = True
-        elif d4 is not None and d2 is not None and d1 is not None:  
-	  totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0)
+        elif d4 is not None and d2 is not None and d1 is not None:
           output = mean(input[d4,:,:],axis=0)
-	  output = output*mesharea; output = sum(output[d2,:],axis=0); output = sum(output[d1],axis=0)/totalarea
+          totalarea = np.ma.masked_where(getmask(output),mesharea)
+          totalarea = mean(totalarea[d2,:],axis=0);totalarea = mean(totalarea[d1])
+          output = output*mesharea; output = mean(output[d2,:],axis=0); output = mean(output[d1])/totalarea
         elif d4 is not None and d2 is not None:
 	  output = mean(input[d4,:,:],axis=0)
-	  totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[d2,:],axis=0)/totalarea
+          totalarea = np.ma.masked_where(getmask(output),mesharea)
+          totalarea = mean(totalarea[d2,:],axis=0)
+          output = output*mesharea; output = mean(output[d2,:],axis=0)/totalarea
         elif d4 is not None and d1 is not None:    output = mean(input[d4,:,:],axis=0); output=mean(output[:,d1],axis=1)
         elif d2 is not None and d1 is not None:
-	  totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0)
-	  output = output*mesharea; output = sum(output[:,d2,:],axis=1); output = sum(output[:,d1],axis=1)/totalarea
+          totalarea = np.tile(mesharea,(output.shape[0],1,1))
+          totalarea = np.ma.masked_where(getmask(output),totalarea)
+          totalarea = mean(totalarea[:,d2,:],axis=1);totalarea = mean(totalarea[:,d1],axis=1)
+          output = output*mesharea; output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)/totalarea
         elif d1 is not None:   output = mean(input[:,:,d1],axis=2)
-        elif d2 is not None:   totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[:,d2,:],axis=1)/totalarea
+        elif d2 is not None:   
+          totalarea = np.tile(mesharea,(output.shape[0],1,1))
+          totalarea = np.ma.masked_where(getmask(output),totalarea)
+          totalarea = mean(totalarea[:,d2,:],axis=1)
+          output = output*mesharea; output = mean(output[:,d2,:],axis=1)/totalarea
         elif d4 is not None:   output = mean(input[d4,:,:],axis=0)
     elif dimension == 4:
-        if mesharea is None: mesharea=np.ones(shape[[2,3]])
+        if mesharea is None: mesharea=np.ones(shape[[2,3]]) # mesharea=np.random.random_sample(shape[[2,3]])*5. + 2. # pour tester
         if   max(d4) >= shape[0]: error = True
         elif max(d3) >= shape[1]: error = True
@@ -227,14 +240,14 @@
              output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
              if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
-	     totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0)
-	     output = output*mesharea
-             output = sum(output[d2,:],axis=0)
-             output = sum(output[d1],axis=0)/totalarea
+             totalarea = np.ma.masked_where(np.isnan(output),mesharea)
+             totalarea = mean(totalarea[d2,:],axis=0); totalarea = mean(totalarea[d1])
+             output = output*mesharea; output = mean(output[d2,:],axis=0); output = mean(output[d1])/totalarea
         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,indice=d3)
              if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
-	     totalarea=sum(mesharea[d2,:],axis=0)
-	     output = output*mesharea; output = sum(output[d2,:],axis=0)/totalarea
+             totalarea = np.ma.masked_where(np.isnan(output),mesharea)
+             totalarea = mean(totalarea[d2,:],axis=0)
+             output = output*mesharea; output = mean(output[d2,:],axis=0)/totalarea
         elif d4 is not None and d3 is not None and d1 is not None: 
              output = mean(input[d4,:,:,:],axis=0)
@@ -246,15 +259,19 @@
              if anomaly:
                  for k in range(output.shape[0]):  output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.)
-	     totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0)
-	     output = output*mesharea; output = sum(output[:,d2,:],axis=1); output = sum(output[:,d1],axis=1)/totalarea
+             totalarea = np.tile(mesharea,(output.shape[0],1,1))
+             totalarea = np.ma.masked_where(getmask(output),mesharea)
+             totalarea = mean(totalarea[:,d2,:],axis=1); totalarea = mean(totalarea[:,d1],axis=1)
+             output = output*mesharea; output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)/totalarea
              #noperturb = smooth1d(output,window_len=7)
              #lenlen = len(output) ; output = output[1:lenlen-7] ; yeye = noperturb[4:lenlen-4]
              #plot(output) ; plot(yeye) ; show() ; plot(output-yeye) ; show()
-        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,indice=d3) 
+        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,indice=d3)
              if anomaly:
                  for k in range(output.shape[0]):  output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.)
-	     totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0)
-	     output = output*mesharea; output = sum(output[:,d2,:],axis=1); output = sum(output[:,d1],axis=1)/totalarea
+             totalarea = np.tile(mesharea,(output.shape[0],1,1))
+             totalarea = np.ma.masked_where(getmask(output),totalarea)
+             totalarea = mean(totalarea[:,d2,:],axis=1); totalarea = mean(totalarea[:,d1],axis=1)
+             output = output*mesharea; output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)/totalarea
         elif d4 is not None and d3 is not None:  
              output = mean(input[d4,:,:,:],axis=0) 
@@ -262,21 +279,32 @@
              if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.) 
         elif d4 is not None and d2 is not None:  
-             output = mean(input[d4,:,:,:],axis=0) 
-	     totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[:,d2,:],axis=1)/totalarea
+             output = mean(input[d4,:,:,:],axis=0)
+             totalarea = np.tile(mesharea,(output.shape[0],1,1))
+             totalarea = np.ma.masked_where(getmask(output),mesharea)
+             totalarea = mean(totalarea[:,d2,:],axis=1)
+             output = output*mesharea; output = mean(output[:,d2,:],axis=1)/totalarea
         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:  
+        elif d3 is not None and d2 is not None:
              output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
-	     totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[:,d2,:],axis=1)/totalarea
+             totalarea = np.tile(mesharea,(output.shape[0],1,1))
+             totalarea = np.ma.masked_where(getmask(output),mesharea)
+             totalarea = mean(totalarea[:,d2,:],axis=1)
+             output = output*mesharea; output = mean(output[:,d2,:],axis=1)/totalarea
         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=2)
-        elif d2 is not None and d1 is not None:  
- 	     totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0)
-	     output = output*mesharea; output = sum(output[:,:,d2,:],axis=2); output = sum(output[:,:,d1],axis=2)/totalarea
+        elif d2 is not None and d1 is not None:
+             totalarea = np.tile(mesharea,(output.shape[0],output.shape[1],1,1))
+             totalarea = np.ma.masked_where(getmask(output),totalarea)
+             totalarea = mean(totalarea[:,:,d2,:],axis=2); totalarea = mean(totalarea[:,:,d1],axis=1)
+             output = output*mesharea; output = mean(output[:,:,d2,:],axis=2); output = mean(output[:,:,d1],axis=1)/totalarea
         elif d1 is not None:        output = mean(input[:,:,:,d1],axis=3)
-        elif d2 is not None: 
-             totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[:,:,d2,:],axis=2)/totalarea
+        elif d2 is not None:
+             totalarea = np.tile(mesharea,(output.shape[0],output.shape[1],1,output.shape[3]))
+             totalarea = np.ma.masked_where(getmask(output),totalarea)
+             totalarea = mean(totalarea[:,:,d2,:],axis=2)
+             output = output*mesharea; output = mean(output[:,:,d2,:],axis=2)/totalarea
         elif d3 is not None:        output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
         elif d4 is not None:        output = mean(input[d4,:,:,:],axis=0)
@@ -739,7 +767,7 @@
     #zecolor = 'r'
     #zelinewidth = 1
-    #step = 5.625
-    #steplon = 5.625
-    #zelatmax = 89.9
+    #step = 180./48.
+    #steplon = 360./64.
+    #zelatmax = 90. - step/3
     if char not in ["moll"]:
         m.drawmeridians(np.r_[-180.:180.:steplon], labels=[0,0,0,1], color=zecolor, linewidth=zelinewidth, fontsize=fontsizemer, latmax=zelatmax)
@@ -932,5 +960,5 @@
              "CO2":          "YlOrBr_r",\
              #### T.N.
-             "MTOT":         "Jet",\
+             "MTOT":         "spectral",\
              "H2O_ICE_S":    "RdBu",\
              "VMR_H2OICE":   "PuBu",\
@@ -1056,32 +1084,4 @@
   return names
 
-## Author: TN [old]
-def getopmatrix (kind,n):
-  import numpy as np
-  # compute matrix of operations between files
-  # the matrix is 'number of files'-square
-  # 1: difference (row minus column), 2: add
-  # not 0 in diag : just plot
-  if n == 1:
-    opm = np.eye(1)
-  elif kind == 'basic':
-    opm = np.eye(n)
-  elif kind == 'difference1': # show differences with 1st file
-    opm = np.zeros((n,n))
-    opm[0,:] = 1
-    opm[0,0] = 0
-  elif kind == 'difference2': # show differences with 1st file AND show 1st file
-    opm = np.zeros((n,n))
-    opm[0,:] = 1
-  else:
-    opm = np.eye(n)
-  return opm
- 
-## Author: TN [old] 
-def checkcoherence (nfiles,nlat,nlon,ntime):
-  if (nfiles > 1):
-     if (nlat > 1):
-        errormess("what you asked is not possible !")
-  return 1
 
 ## Author: TN
Index: trunk/UTIL/PYTHON/myscript.py
===================================================================
--- trunk/UTIL/PYTHON/myscript.py	(revision 642)
+++ trunk/UTIL/PYTHON/myscript.py	(revision 647)
@@ -74,4 +74,5 @@
     ### SPECIAL
     parser.add_option('--tsat',         action='store_true',dest='tsat',               default=False,help='convert temperature field T in Tsat-T using pressure')
+    parser.add_option('--stream',       action='store_true',dest='stream',             default=False,help='plot streamlines from streamfunction.e for vert/lat slices.')
 
     return parser
Index: trunk/UTIL/PYTHON/planetoplot.py
===================================================================
--- trunk/UTIL/PYTHON/planetoplot.py	(revision 642)
+++ trunk/UTIL/PYTHON/planetoplot.py	(revision 647)
@@ -66,5 +66,6 @@
            lbls=None,\
            lstyle=None,\
-           cross=None):
+           cross=None,\
+           streamflag=False):
 
     ####################################################################################################################
@@ -91,4 +92,5 @@
     #from singlet import singlet
     from itertools import cycle
+    import os
 
 #########################
@@ -414,4 +416,5 @@
            if mult != 2718.:  what_I_plot = what_I_plot*mult 
            else:              what_I_plot = np.log10(what_I_plot) ; print "log plot"
+                      
        if var2:      ### what is contoured.
            what_I_plot_contour, error = reducefield( all_var2[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert, \
@@ -523,5 +526,5 @@
                            if indexvert is not None: lbl = lbl + " iz" + str(indexvert[0])
                            if indextime is not None: lbl = lbl + " it" + str(indextime[0])
-                           if lbl == "": lbl = namefiles[index_f]
+                           if lbl == "": lbl = all_namefile[index_f]
 
                         if mrate is not None: x = y  ## because swapaxes...
@@ -539,4 +542,25 @@
 
                     elif which == "regular":
+                    
+                        # plot stream lines if there is a stream file and a vert/lat slice. Might not work with movies ??
+                        if streamflag and sslat is None and svert is None:
+                             streamfile = all_namefile[index_f].replace('.nc','_stream.nc')
+                             teststream = os.path.exists(streamfile)
+                             if teststream:
+                                print 'INFO: Using stream file',streamfile, 'for stream lines'
+                                ncstream = Dataset(streamfile)
+                                psi = getfield(ncstream,'psi')
+                                psi = psi[0,:,:,0] # time and longitude are dummy dimensions
+                                if psi.shape[1] != len(x) or psi.shape[0] != len(y):
+                                    errormess('stream file does not fit! Dimensions: '+str(psi.shape)+' '+str(x.shape)+' '+str(y.shape))
+                                zelevels = np.arange(-1.e10,1.e10,1.e9)
+                                zemin = np.min(abs(zelevels))
+                                zemax = np.max(abs(zelevels))
+                                zewidth  =  (abs(zelevels)-zemin)*(5.- 0.5)/(zemax - zemin) + 0.5 # linewidth ranges from 5 to 0.5
+                                cs = contour( x,y,psi, zelevels, colors='k', linewidths = zewidth)
+                                clabel(cs, inline=True, fontsize = 4.*rcParams['font.size']/5., fmt="%1.1e")
+                             else:
+                                print 'WARNING: STREAM FILE',streamfile, 'DOES NOT EXIST !'
+                             
                         if hole:         what_I_plot_frame = hole_bounds(what_I_plot_frame,zevmin,zevmax)
                         else:            what_I_plot_frame = bounds(what_I_plot_frame,zevmin,zevmax)
@@ -554,4 +578,5 @@
                             if mapmode == 1:       m.pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax, alpha=trans)
                             elif mapmode == 0:     pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax, alpha=trans)
+                        
 
                         if colorb not in ['nobar','onebar']:        
@@ -591,6 +616,5 @@
                             zebins = np.linspace(zevmin,zevmax,num=30)
                             hist(toplot,bins=zebins,histtype='step',linewidth=2,color='k',normed=True)
-                            zestd = np.std(toplot)
-                            zemean = np.mean(toplot)
+                            zestd = np.std(toplot);zemean = np.mean(toplot)
                             zebins = np.linspace(zevmin,zevmax,num=300)
                             zegauss = (1./(zestd * np.sqrt(2 * np.pi)) * np.exp( - (zebins - zemean)**2 / (2 * zestd**2) ) )
@@ -675,6 +699,6 @@
                 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+' '+namefiles[0]#index_f]
-            else:                            plottitle = basename+' '+namefiles[0]#index_f]
+                else:                        plottitle = basename+' '+all_namefile[index_f]
+            else:                            plottitle = basename+' '+all_namefile[index_f]
        if mult != 1:                         plottitle = '{:.0e}'.format(mult) + "*" + plottitle
        if zetitle[0] != "fill":                 
Index: trunk/UTIL/PYTHON/pp.py
===================================================================
--- trunk/UTIL/PYTHON/pp.py	(revision 642)
+++ trunk/UTIL/PYTHON/pp.py	(revision 647)
@@ -172,5 +172,5 @@
                 mrate=opt.rate,mquality=opt.quality,trans=opt.trans,zarea=opt.area,axtime=opt.axtime,\
                 redope=opt.redope,seevar=opt.seevar,xlab=opt.xlab,ylab=opt.ylab,lbls=separatenames(opt.labels),\
-                lstyle=separatenames(opt.linestyle),cross=readslices(opt.mark))
+                lstyle=separatenames(opt.linestyle),cross=readslices(opt.mark),streamflag=opt.stream)
         print 'DONE: '+name
         system("rm -f to_be_erased")
