Index: trunk/MESOSCALE/PLOT/PYTHON/domain.py
===================================================================
--- trunk/MESOSCALE/PLOT/PYTHON/domain.py	(revision 179)
+++ 	(revision )
@@ -1,25 +1,0 @@
-#!/usr/bin/env python
-
-### A. Spiga LMD 29/05/2011
-
-def usage():
-    print 'USAGE:  plot.py file '
-    print 'file  : name of input netcdf file.'
-    print 'Example:  domain.py /tmp7/aslmd/HELLAS/wrfout_d02_2024-10-03_06:00:00'
-
-def domain (namefile): 
-    from netCDF4 import Dataset
-    from myplot import getcoord2d,define_proj,makeplotpng,simplinterv
-    from matplotlib.pyplot import contourf
-    nc  = Dataset(namefile)
-    [lon2d,lat2d] = getcoord2d(nc)
-    [wlon,wlat] = simplinterv(lon2d,lat2d)
-    m = define_proj("ortho",wlon,wlat,back="vishires")
-    x, y = m(lon2d, lat2d)
-    contourf(x, y, nc.variables['HGT'][0,:,:] / 1000.,50)
-    makeplotpng(namefile+".domain",pad_inches_value=0.15)
-
-if __name__ == "__main__":
-    import sys
-    domain(str(sys.argv[1]))
-
Index: trunk/MESOSCALE/PLOT/PYTHON/map_ecmwf.py
===================================================================
--- trunk/MESOSCALE/PLOT/PYTHON/map_ecmwf.py	(revision 179)
+++ trunk/MESOSCALE/PLOT/PYTHON/map_ecmwf.py	(revision 180)
@@ -11,9 +11,9 @@
 ##########################################################################
 proj = "cyl"  #"moll" "ortho" "lcc"
-proj = "ortho"
+#proj = "ortho"
 area = "Europe"
-area = "Central_America"
-area = "Southern_Hemisphere"
-area = "Northern_Hemisphere"
+#area = "Central_America"
+#area = "Southern_Hemisphere"
+#area = "Northern_Hemisphere"
 ##########################################################################
 
Index: trunk/MESOSCALE/PLOT/PYTHON/mylib/ecmwf.py
===================================================================
--- trunk/MESOSCALE/PLOT/PYTHON/mylib/ecmwf.py	(revision 180)
+++ trunk/MESOSCALE/PLOT/PYTHON/mylib/ecmwf.py	(revision 180)
@@ -0,0 +1,101 @@
+import urllib
+import urllib2
+import time
+import datetime
+
+class ECMWFDataServer:
+    def __init__(self,portal,token,email):
+        self.version = '0.3'
+        self.portal  = portal
+        self.token   = token
+        self.email   = email
+
+    def _call(self,action,args):
+
+        params = {'_token'   : self.token,
+                  '_email'   : self.email,
+                  '_action'  : action,
+                  '_version' : self.version}
+        params.update(args)
+            
+        data = urllib.urlencode(params)
+        req = urllib2.Request(self.portal, data)
+        response = urllib2.urlopen(req)
+
+        json = response.read();
+
+        undef = None;
+        json = eval(json)
+
+        if json != None:
+            if 'error' in json:
+                raise RuntimeError(json['error'])
+            if 'message' in json:
+                self.put(json['message'])
+
+        return json
+
+    def put(self,*args):
+        print datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S'),
+        for a in args:
+            print a,
+        print
+
+        
+    def retrieve(self,args):
+        self.put("ECMWF data server batch tool version",self.version);
+        user = self._call("user_info",{});
+        self.put("Welcome to",user['name'], "from" , user['organisation']);
+
+        r = self._call('retrieve',args)
+        rid = r['request']
+
+        last  = ''
+        sleep = 0
+        while r['status'] != 'complete' and r['status'] != 'aborted':
+            text = r['status'] + '.'
+            if 'info' in r and r['info'] != None:
+                text = text + ' ' + r['info'] 
+
+            if text != last:
+                self.put("Request",text)
+                last = text
+
+            time.sleep(sleep)
+            r = self._call('status',{'request':rid})
+            if sleep < 60:
+                sleep = sleep + 1
+
+        if r['status'] != last:
+            self.put("Request",r['status'])
+
+        if 'reason' in r:
+            for m in r['reason']:
+                self.put(m)
+
+        if 'result' in r:
+            size = long(r['size'])
+            self.put("Downloading",self._bytename(size))
+            done = self._transfer(r['result'],args['target'])
+            self.put("Done")
+            if done != size:
+                raise RuntimeError("Size mismatch: " + str(done) + " and " + str(size))
+
+        self._call('delete',{'request':rid})
+        
+        if r['status'] == 'aborted':
+            raise RuntimeError("Request aborted")
+        
+    def _transfer(self,url,path):
+        result =  urllib.urlretrieve(url,path)
+        return long(result[1]['content-length'])
+        
+    def _bytename(self,size):   
+        next = {'':'K','K':'M','M':'G','G':'T','T':'P'}
+        l    = ''
+        size = size*1.0
+        while 1024 < size:
+            l = next[l]
+            size = size / 1024
+        return "%g %sbyte%s" % (size,l,'s')
+
Index: trunk/MESOSCALE/PLOT/PYTHON/mylib/myecmwf.py
===================================================================
--- trunk/MESOSCALE/PLOT/PYTHON/mylib/myecmwf.py	(revision 180)
+++ trunk/MESOSCALE/PLOT/PYTHON/mylib/myecmwf.py	(revision 180)
@@ -0,0 +1,51 @@
+def ecmwf_title_field (var):
+	### http://www.ecmwf.int/services/archive/d/table/grib_table_2_versions
+	### http://www.ecmwf.int/services/archive/d/parameters/mars=1/order=grib_parameter/table=128/
+	### http://www.atmos.albany.edu/facstaff/rmctc/ecmwfGrib/ecmwfgrib128.tbl
+	if   var == "151": name="Mean Sea Level Pressure (Pa)"          # ["151", "fc", "3", "sfc"]
+	elif var == "146": name="Sensible Heat Flux (W m^-2)"           # ["146", "fc", "3", "sfc"]
+	elif var == "130": name="Atmospheric Temperature (K)"           # ["130", "an", "0", "pl" ]
+	elif var == "167": name="2m Atmospheric Temperature (K)"        # ["167", "fc", "3", "sfc"]
+	return name
+
+def split_char_add (var,add):
+	import  numpy as np
+	varchar = ""
+	for i in range( np.array(var).size ):
+		varchar = varchar + str(var[i]) + add + "/"
+	return varchar
+
+def split_char (var):
+	varchar = split_char_add (var,"")
+	return varchar
+
+def get_ecmwf (var, dataset, wlat, wlon, lev, date, tim):
+	from ecmwf 	import ECMWFDataServer
+	from os		import system
+	from netCDF4	import Dataset
+	gbfile = 'output.grib'
+	ncfile = 'output.nc'
+	######
+	timchar = split_char (tim)
+	levchar = split_char (lev)
+	varchar = split_char_add (var,".128")
+	##########################################################################
+	## First registrer at http://data-portal.ecmwf.int/data/d/license/interim/
+	## Then get your token at http://data-portal.ecmwf.int/data/d/token/interim_daily/
+	server = ECMWFDataServer('http://data-portal.ecmwf.int/data/d/dataserver/','6948746482e9e3e29d64211e06329a2e','spiga@lmd.jussieu.fr')
+	server.retrieve({
+		'dataset'  : "interim_daily",
+		'date'     : date[2] + date[1] + date[0] + "/to/" + date[5] + date[4] + date[3],
+		'time'     : timchar,
+		'step'     : dataset[2],
+		'levtype'  : dataset[1],
+		'type'     : dataset[0],
+		'param'    : varchar,
+		'levelist' : levchar,   
+		'area'     : str(wlat[0])+"/"+str(wlon[0])+"/"+str(wlat[1])+"/"+str(wlon[1]),
+		'target'   : gbfile
+	})
+	system("cdo -f nc copy "+gbfile+" "+ncfile+" ; rm -f "+gbfile)
+	nc = Dataset(ncfile)
+	#system("rm -f "+ncfile)
+	return nc
Index: trunk/MESOSCALE/PLOT/PYTHON/mylib/mymath.py
===================================================================
--- trunk/MESOSCALE/PLOT/PYTHON/mylib/mymath.py	(revision 180)
+++ trunk/MESOSCALE/PLOT/PYTHON/mylib/mymath.py	(revision 180)
@@ -0,0 +1,15 @@
+def min (field,axis=None): 
+        import numpy as np
+        return np.array(field).min(axis=axis)
+
+def max (field,axis=None):
+        import numpy as np
+        return np.array(field).max(axis=axis)
+
+def mean (field,axis=None):
+        import numpy as np
+        return np.array(field).mean(axis=axis)
+
+def deg ():
+        return u'\u00b0'
+
Index: trunk/MESOSCALE/PLOT/PYTHON/mylib/myplot.py
===================================================================
--- trunk/MESOSCALE/PLOT/PYTHON/mylib/myplot.py	(revision 180)
+++ trunk/MESOSCALE/PLOT/PYTHON/mylib/myplot.py	(revision 180)
@@ -0,0 +1,175 @@
+def latinterv (area):
+	if   area == "Europe": 
+		wlat = [20.,80.]
+		wlon = [-50.,50.]
+	elif area == "Central_America":
+		wlat = [-10.,40.]
+		wlon = [230.,300.]
+	elif area == "Africa":
+		wlat = [-20.,50.]
+		wlon = [-50.,50.]
+	elif area == "Whole":
+		wlat = [-90.,90.]
+		wlon = [-180.,180.]
+	elif area == "Southern_Hemisphere":
+		wlat = [-90.,60.]
+		wlon = [-180.,180.]
+        elif area == "Northern_Hemisphere":
+                wlat = [-60.,90.]
+                wlon = [-180.,180.]
+        elif area == "Tharsis":
+		wlat = [-30.,60.]
+		wlon = [-170.,-10.]
+	elif area == "Whole_No_High":
+		wlat = [-60.,60.]
+		wlon = [-180.,180.]
+	elif area == "Chryse":
+                wlat = [-60.,60.]
+                wlon = [-60.,60.]
+        elif area == "North_Pole":
+                wlat = [60.,90.]
+                wlon = [-180.,180.]
+        elif area == "Close_North_Pole":
+                wlat = [75.,90.]
+                wlon = [-180.,180.]
+	return wlon,wlat
+
+def ptitle (name):
+    from matplotlib.pyplot import title
+    title(name)
+    print name
+
+def simplinterv (lon2d,lat2d):
+    import numpy as np
+    return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]]
+
+def makeplotpngres (filename,res,pad_inches_value=0.25,folder='',disp=True):
+    import  matplotlib.pyplot as plt
+    res = int(res)
+    if folder != '':	name = folder+'/'+filename+str(res)+".png"
+    else:           	name = filename+str(res)+".png"
+    plt.savefig(name,dpi=res,bbox_inches='tight',pad_inches=pad_inches_value)
+    if disp:		display(name)		
+    return
+
+def makeplotpng (filename,pad_inches_value=0.25,minres=100.,folder=''):
+    makeplotpngres(filename,minres,     pad_inches_value=pad_inches_value,folder=folder)
+    makeplotpngres(filename,minres+200.,pad_inches_value=pad_inches_value,folder=folder,disp=False)
+    return
+
+def getcoord2d (nc,nlat='XLAT',nlon='XLONG'):
+	import numpy as np
+	lat = nc.variables[nlat][0,:,:]
+	lon = nc.variables[nlon][0,:,:]
+	if np.array(lat).ndim != 2:	[lon2d,lat2d] = np.meshgrid(lon,lat)
+	else:				[lon2d,lat2d] = [lon,lat]
+	return lon2d,lat2d
+
+def smooth (field, coeff):
+	## actually blur_image could work with different coeff on x and y
+	if coeff > 1:	result = blur_image(field,int(coeff))
+	else:		result = field
+	return result
+
+def gauss_kern(size, sizey=None):
+	import numpy as np
+	## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth	
+    	# Returns a normalized 2D gauss kernel array for convolutions
+    	size = int(size)
+    	if not sizey:
+	        sizey = size
+	else:
+	        sizey = int(sizey)
+	x, y = np.mgrid[-size:size+1, -sizey:sizey+1]
+	g = np.exp(-(x**2/float(size)+y**2/float(sizey)))
+	return g / g.sum()
+
+def blur_image(im, n, ny=None) :
+	from scipy.signal import convolve
+	## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
+	# blurs the image by convolving with a gaussian kernel of typical size n. 
+	# The optional keyword argument ny allows for a different size in the y direction.
+    	g = gauss_kern(n, sizey=ny)
+    	improc = convolve(im, g, mode='same')
+    	return improc
+
+def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1):
+	## scale regle la reference du vecteur
+	## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs.
+	import  matplotlib.pyplot               as plt
+	import  numpy                           as np
+	posx = np.max(x)*0.90 
+	posy = np.mean(y)
+	u = smooth(u,csmooth)
+	v = smooth(v,csmooth)
+	q = plt.quiver(	x[::stride,::stride],\
+			y[::stride,::stride],\
+			u[::stride,::stride],\
+			v[::stride,::stride],\
+			angles='xy',color=color,\
+			scale=factor,width=0.003 )
+	if color=='white':	kcolor='black'
+	elif color=='yellow':	kcolor=color
+	else:			kcolor=color
+	p = plt.quiverkey(q,posx,posy,scale,\
+			str(int(scale)),coordinates='data',color=kcolor)
+	return p
+
+def display (name):
+	from 	os 				import system
+	system("eog "+name+" > /dev/null 2> /dev/null &")
+	return name
+
+def findstep (wlon):
+    steplon = int((wlon[1]-wlon[0])/3.)
+    step = 60.
+    if steplon < 60.: step = 30.
+    if steplon < 30.: step = 15.
+    if steplon < 15.: step = 10.
+    if steplon < 10.: step = 5.
+    if steplon < 5.:  step = 1.
+    return step
+
+def define_proj (char,wlon,wlat,back="."):
+    from    mpl_toolkits.basemap            import Basemap
+    import  numpy                           as np
+    import  matplotlib                      as mpl
+    meanlon = 0.5*(wlon[0]+wlon[1])
+    meanlat = 0.5*(wlat[0]+wlat[1])
+    h = 2000.
+    if   char == "cyl":     m = Basemap(projection='cyl',llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])#,suppress_ticks=False)
+    elif char == "moll":    m = Basemap(projection='moll',lon_0=meanlon)
+    elif char == "ortho":   m = Basemap(projection='ortho',lon_0=meanlon,lat_0=meanlat)
+    elif char == "lcc":     m = Basemap(projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\
+                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
+    elif char == "npstere": m = Basemap(projection='npstere', boundinglat=wlat[0], lon_0=0.)
+    elif char == "spstere": m = Basemap(projection='spstere', boundinglat=wlat[0], lon_0=0.)
+    elif char == "nsper":   m = Basemap(projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=h*1000.)
+    fontsizemer = int(mpl.rcParams['font.size']*2./3.)
+    if char in ["cyl","lcc"]:   step = findstep(wlon)
+    else:                       step = 10.
+    m.drawmeridians(np.r_[-180.:180.:step*2.], 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)
+    if back == ".":      m.warpimage(marsmap(),scale=0.75)
+    elif back == None:   pass 
+    else:                m.warpimage(marsmap(back),scale=0.75)
+    return m
+
+def marsmap (whichone="vishires"):
+	whichlink = 	{ \
+		"vis":		"http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\
+		"vishires":	"http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/MarsMap_2500x1250.jpg",\
+		"mola":		"http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\
+		"molabw":	"http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/MarsElevation_2500x1250.jpg",\
+			}
+	if whichone not in whichlink: 
+		print "marsmap: choice not defined... you'll get the default one... "
+		whichone = "vishires"  
+        return whichlink[whichone]
+
+def earthmap (whichone):
+	if   whichone == "contrast":	whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg"
+	elif whichone == "bw":		whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg"
+	elif whichone == "nice":	whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg"
+	return whichlink
+
Index: trunk/MESOSCALE/PLOT/PYTHON/scripts/domain.py
===================================================================
--- trunk/MESOSCALE/PLOT/PYTHON/scripts/domain.py	(revision 180)
+++ trunk/MESOSCALE/PLOT/PYTHON/scripts/domain.py	(revision 180)
@@ -0,0 +1,25 @@
+#!/usr/bin/env python
+
+### A. Spiga LMD 29/05/2011
+
+def usage():
+    print 'USAGE:  plot.py file '
+    print 'file  : name of input netcdf file.'
+    print 'Example:  domain.py /tmp7/aslmd/HELLAS/wrfout_d02_2024-10-03_06:00:00'
+
+def domain (namefile,proj="ortho",back="molabw"): #vishires
+    from netCDF4 import Dataset
+    from myplot import getcoord2d,define_proj,makeplotpng,simplinterv
+    from matplotlib.pyplot import contourf
+    nc  = Dataset(namefile)
+    [lon2d,lat2d] = getcoord2d(nc)
+    [wlon,wlat] = simplinterv(lon2d,lat2d)
+    m = define_proj(proj,wlon,wlat,back=back)
+    x, y = m(lon2d, lat2d)
+    contourf(x, y, nc.variables['HGT'][0,:,:] / 1000., 50)
+    makeplotpng(namefile+".domain",pad_inches_value=0.05)
+
+if __name__ == "__main__":
+    import sys
+    domain(str(sys.argv[1]))
+
Index: trunk/MESOSCALE/PLOT/PYTHON/scripts/winds.py
===================================================================
--- trunk/MESOSCALE/PLOT/PYTHON/scripts/winds.py	(revision 180)
+++ trunk/MESOSCALE/PLOT/PYTHON/scripts/winds.py	(revision 180)
@@ -0,0 +1,42 @@
+#!/usr/bin/env python
+
+### A. Spiga LMD 29/05/2011
+
+def usage():
+    print 'USAGE:  winds.py ntime nvert file'
+    print 'file  : name of input netcdf file.'
+    print 'nvert : vertical level.'
+    print 'Example: winds.py 0 /d5/aslmd/LMD_MM_MARS_SIMUS/OM/OM6_TI85/wrfout_d01_2024-06-43_06:00:00_zabg'
+
+def winds (namefile,nvert,proj="cyl",back=None): 
+    from netCDF4 import Dataset
+    from myplot import getcoord2d,define_proj,makeplotpng,simplinterv,vectorfield,ptitle
+    from matplotlib.pyplot import contourf, subplot, figure
+    nc  = Dataset(namefile)
+    [lon2d,lat2d] = getcoord2d(nc)
+    [wlon,wlat] = simplinterv(lon2d,lat2d)
+    [u,v] = getwinds(nc)
+    nt = len(u[:,0,0,0])
+    fig = figure()
+    fig.subplots_adjust(wspace = 0.0, hspace = 0.3)
+    sub = 221
+    for i in range(0,nt-1,int(nt/4.)):
+       subplot(sub)
+       ptitle("Winds time"+str(i)+" level"+str(nvert))
+       m = define_proj(proj,wlon,wlat,back=back)
+       x, y = m(lon2d, lat2d)
+       contourf(x, y, nc.variables['HGT'][0,:,:] / 1000., 30)
+       vectorfield(u[i,nvert,:,:], v[i,nvert,:,:],\
+                      x, y, stride=5, csmooth=3,\
+                      scale=20., factor=300., color='k')
+       sub += 1
+    makeplotpng(namefile+".winds"+str(nvert),pad_inches_value=0.35)   
+
+def getwinds (nc,charu='Um',charv='Vm'):
+    u = nc.variables[charu]
+    v = nc.variables[charv]
+    return u,v
+
+if __name__ == "__main__":
+    import sys
+    winds(str(sys.argv[2]),int(sys.argv[1]))
