Index: trunk/UTIL/PYTHON/mymath.py
===================================================================
--- trunk/UTIL/PYTHON/mymath.py	(revision 518)
+++ trunk/UTIL/PYTHON/mymath.py	(revision 525)
@@ -38,4 +38,24 @@
            else: 
               return np.array(field).mean(axis=axis)
+
+def sum (field,axis=None):
+        import numpy as np
+        if field is None: return None
+        else:
+           if type(field).__name__=='MaskedArray':
+              field.set_fill_value(np.NaN)
+              zout=np.ma.array(field).sum(axis=axis)
+              if axis is not None:
+                 zout.set_fill_value(np.NaN)
+                 return zout.filled()
+              else:return zout
+           elif (np.isnan(np.sum(field)) and (type(field).__name__ not in 'MaskedArray')):
+              zout=np.ma.masked_invalid(field).sum(axis=axis)
+              if axis is not None:
+                 zout.set_fill_value([np.NaN])
+                 return zout.filled()
+              else:return zout
+           else:
+              return np.array(field).sum(axis=axis)
 
 def deg ():
Index: trunk/UTIL/PYTHON/myplot.py
===================================================================
--- trunk/UTIL/PYTHON/myplot.py	(revision 518)
+++ trunk/UTIL/PYTHON/myplot.py	(revision 525)
@@ -36,8 +36,9 @@
     return ltst
 
-## Author: AS, AC
+## Author: AS, AC, JL
 def whatkindfile (nc):
     if 'controle' in nc.variables:             typefile = 'gcm'
     elif 'phisinit' in nc.variables:           typefile = 'gcm'
+    elif 'time_counter' in nc.variables:       typefile = 'earthgcm'
     elif hasattr(nc,'START_DATE'):
       if '9999' in getattr(nc,'START_DATE') :  typefile = 'mesoideal'
@@ -77,10 +78,10 @@
 
 ## Author: AS + TN + AC
-def reducefield (input,d4=None,d3=None,d2=None,d1=None,yint=False,alt=None,anomaly=False,redope=None):
+def reducefield (input,d4=None,d3=None,d2=None,d1=None,yint=False,alt=None,anomaly=False,redope=None,mesharea=None):
     ### we do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x"
     ### it would be actually better to name d4 d3 d2 d1 as t z y x
     ### ... note, anomaly is only computed over d1 and d2 for the moment
     import numpy as np
-    from mymath import max,mean,min
+    from mymath import max,mean,min,sum
     csmooth = 12 ## a fair amount of grid points (too high results in high computation time)
     if redope is not None:
@@ -95,5 +96,5 @@
        #elif redope == "maxx":     input = max(input,axis=3) ; d4 = None
     dimension = np.array(input).ndim
-    shape = np.array(input).shape
+    shape = np.array(np.array(input).shape)
     #print 'd1,d2,d3,d4: ',d1,d2,d3,d4
     if anomaly: print 'ANOMALY ANOMALY'
@@ -107,22 +108,33 @@
     ### now the main part
     if dimension == 2:
-        if   d2 >= shape[0]: error = True
-        elif d1 >= shape[1]: error = True
-        elif d1 is not None and d2 is not None:  output = mean(input[d2,:],axis=0); output = mean(output[d1],axis=0)
+        if mesharea is None: mesharea=np.ones(shape)
+        if   max(d2) >= shape[0]: error = True
+        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
         elif d1 is not None:         output = mean(input[:,d1],axis=1)
-        elif d2 is not None:         output = mean(input[d2,:],axis=0)
+        elif d2 is not None:         totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[d2,:],axis=0)/totalarea
     elif dimension == 3:
+        if mesharea is None: mesharea=np.ones(shape[[1,2]])
         if   max(d4) >= shape[0]: error = True
         elif max(d2) >= shape[1]: error = True
         elif max(d1) >= shape[2]: error = True
         elif d4 is not None and d2 is not None and d1 is not None:  
-            output = mean(input[d4,:,:],axis=0); output = mean(output[d2,:],axis=0); output = mean(output[d1],axis=0)
-        elif d4 is not None and d2 is not None:    output = mean(input[d4,:,:],axis=0); output=mean(output[d2,:],axis=0)
+	  totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0)
+          output = mean(input[d4,:,:],axis=0)
+	  output = output*mesharea; output = sum(output[d2,:],axis=0); output = sum(output[d1],axis=0)/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
         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:    output = mean(input[:,d2,:],axis=1); output=mean(output[:,d1],axis=1)
-        elif d1 is not None:                       output = mean(input[:,:,d1],axis=2)
-        elif d2 is not None:                       output = mean(input[:,d2,:],axis=1)
-        elif d4 is not None:                       output = mean(input[d4,:,:],axis=0)
+        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
+        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 d4 is not None:   output = mean(input[d4,:,:],axis=0)
     elif dimension == 4:
+        if mesharea is None: mesharea=np.ones(shape[[2,3]])
         if   max(d4) >= shape[0]: error = True
         elif max(d3) >= shape[1]: error = True
@@ -133,11 +145,14 @@
              output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
              if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
-             output = mean(output[d2,:],axis=0)
-             output = mean(output[d1],axis=0)
+	     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
         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.) 
-             output = mean(output[d2,:],axis=0)
+             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
         elif d4 is not None and d3 is not None and d1 is not None: 
              output = mean(input[d4,:,:,:],axis=0)
@@ -149,6 +164,6 @@
              if anomaly:
                  for k in range(output.shape[0]):  output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.)
-             output = mean(output[:,d2,:],axis=1)
-             output = mean(output[:,d1],axis=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
              #noperturb = smooth1d(output,window_len=7)
              #lenlen = len(output) ; output = output[1:lenlen-7] ; yeye = noperturb[4:lenlen-4]
@@ -158,6 +173,6 @@
              if anomaly:
                  for k in range(output.shape[0]):  output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.)
-             output = mean(output[:,d2,:],axis=1)
-             output = mean(output[:,d1],axis=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
         elif d4 is not None and d3 is not None:  
              output = mean(input[d4,:,:,:],axis=0) 
@@ -166,5 +181,5 @@
         elif d4 is not None and d2 is not None:  
              output = mean(input[d4,:,:,:],axis=0) 
-             output = mean(output[:,d2,:],axis=1)
+	     totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[:,d2,:],axis=1)/totalarea
         elif d4 is not None and d1 is not None:  
              output = mean(input[d4,:,:,:],axis=0) 
@@ -172,13 +187,14 @@
         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)
+	     totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(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:  
-             output = mean(input[:,:,d2,:],axis=2)
-             output = mean(output[:,:,d1],axis=2)
+ 	     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 d1 is not None:        output = mean(input[:,:,:,d1],axis=3)
-        elif d2 is not None:        output = mean(input[:,:,d2,:],axis=2)
+        elif d2 is not None: 
+             totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(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)
@@ -427,4 +443,6 @@
     elif typefile in ['gcm']: 
         [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True)
+    elif typefile in ['earthgcm']: 
+        [lon2d,lat2d] = getcoord2d(nc,nlat="lat",nlon="lon",is1d=True)
     elif typefile in ['geo']:
         [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M')
Index: trunk/UTIL/PYTHON/planetoplot.py
===================================================================
--- trunk/UTIL/PYTHON/planetoplot.py	(revision 518)
+++ trunk/UTIL/PYTHON/planetoplot.py	(revision 525)
@@ -9,5 +9,5 @@
 ### A. Colaitis  -- LMD --    12/2011 -- Added movie capability [mencoder must be installed]
 ### A. Spiga     -- LMD --    12/2011 -- Added HTML animated page capability + general tests of consistency [winds, etc...] + consistent generic movie loop
-
+### J. Leconte   -- LMD --    02/2012 -- Added area weighted averaging. Compatibility with terrestrial gcm.
 def planetoplot (namefiles,\
            level=0,\
@@ -126,8 +126,9 @@
       ### ... TYPEFILE
       typefile = whatkindfile(nc)
-      if typefile in ['mesoideal']:   mapmode=0;winds=False
-      elif typefile in ['gcm'] and len(nc.variables["longitude"][:]) is 1 and len(nc.variables["latitude"][:]) is 1:       mapmode=0;winds=False
+      if typefile in ['mesoideal']:   mapmode=0
+      elif typefile in ['gcm'] and len(nc.variables["longitude"][:]) is 1 and len(nc.variables["latitude"][:]) is 1:       mapmode=0
+      elif typefile in ['earthgcm'] and len(nc.variables["lon"][:]) is 1 and len(nc.variables["lat"][:]) is 1:       mapmode=0
       if mapmode == 0:       winds=False
-      elif mapmode == 1:
+      elif mapmode == 1 and redope is None:
           if svert is None:  svert = readslices(str(level)) ; nvert=1
           if stime is None and mrate is None:
@@ -166,9 +167,15 @@
 ############ LOAD 4D DIMENSIONS : x, y, z, t #############
 ##########################################################
-      if typefile == "gcm":
-          lat = nc.variables["latitude"][:] ; lon = nc.variables["longitude"][:] ; vert = nc.variables["altitude"][:]
-          if "Time" in nc.variables:      time = nc.variables["Time"][:]
-          elif "time" in nc.variables:    time = nc.variables["time"][:]
-          else:                           errormess("no time axis found.")
+      if typefile in ["gcm","earthgcm"]:
+          ### SPACE
+          if typefile == "gcm":         lat = nc.variables["latitude"][:] ; lon = nc.variables["longitude"][:] ; vert = nc.variables["altitude"][:]
+          elif typefile == "earthgcm":  lat = nc.variables["lat"][:] ; lon = nc.variables["lon"][:] ; vert = nc.variables["Alt"][:]
+          if "aire" in nc.variables:      area = nc.variables["aire"][:,:]  #JL to weight means with the area
+          else:                           area = None
+          ### TIME
+          if "Time" in nc.variables:      	time = nc.variables["Time"][:]
+          elif "time" in nc.variables:    	time = nc.variables["time"][:]
+          elif "time_counter" in nc.variables:	time = nc.variables["time_counter"][:]/86400. #### time counter cinverstion from s-> days
+          else:                           	errormess("no time axis found.")
           if axtime in ["ls","sol"]:   errormess("not supported. should not be too difficult though.")
           # for 1D plots (no need for longitude computation):
@@ -178,4 +185,5 @@
               print "LOCAL TIMES.... ", time
       elif typefile in ['meso','mesoapi','geo','mesoideal']:
+          area = None ## not active for the moment
           ###### STUFF TO GET THE CORRECT LAT/LON FROM MESOSCALE FILES WITH 2D LAT/LON ARRAYS
           ###### principle: calculate correct indices then repopulate slon and slat
@@ -341,5 +349,5 @@
        if varname:   ### what is shaded.
            what_I_plot, error = reducefield( all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \
-                                             yint=yintegral, alt=vert, anomaly=anomaly, redope=redope )
+                                             yint=yintegral, alt=vert, anomaly=anomaly, redope=redope, mesharea=area )
            if mult != 2718.:  what_I_plot = what_I_plot*mult 
            else:              what_I_plot = np.log10(what_I_plot) ; print "log plot"
