Index: trunk/UTIL/PYTHON/mcd/examples/dustdevil.py
===================================================================
--- trunk/UTIL/PYTHON/mcd/examples/dustdevil.py	(revision 971)
+++ trunk/UTIL/PYTHON/mcd/examples/dustdevil.py	(revision 1007)
@@ -7,4 +7,7 @@
 
 from math import isnan
+
+from ppplot import writeascii
+
 
 dafile = "/home/aymeric/Work/submitted/coauthorship/2012reiss/reiss_data.txt"
@@ -48,8 +51,15 @@
 #exit()
 
-query = mcd()
-mcdwind = [] ; mcdwindle = [] ; mcdwindhe = [] ; mcdangle = [] ; mcdanglele = [] ; mcdanglehe = [] 
-
-for i in yorgl[0]:
+## constant_heights
+#for heights in [20.,100.,500.,1000.,2000.,5000.,7000.,10000.,20000.]:
+for heights in [10.]:
+
+## multiple_obsheight
+#for heights in [0.1,0.2,0.5,0.75,1.0,1.5,2.0,10.]:
+
+ query = mcd()
+ mcdwind = [] ; mcdwindle = [] ; mcdwindhe = [] ; mcdangle = [] ; mcdanglele = [] ; mcdanglehe = [] 
+
+ for i in yorgl[0]:
 
    if not isnan(ddwind[i-1]):
@@ -58,6 +68,13 @@
     query.xdate = yorgl[3][i-1]
     query.loct = yorgl[4][i-1]
-    query.xz = ddheight[i-1]
-    #query.xz = 1000.  #ddheight is really the one that works best
+
+    query.xz = ddheight[i-1] #ddheight is really the one that works best
+
+    ## constant_heights
+    query.xz = heights
+
+    ## multiple_obsheight
+    #query.xz = ddheight[i-1]*heights
+
     query.update() ; uu = query.zonwind ; vv = query.merwind
     if uu == -999. or vv == -999.: 
@@ -70,6 +87,9 @@
        if angle < 0.: angle = angle + 360.
 
-       query.xzs = max(ddheight[i-1] - 2.*ddheighterror[i-1],10.)
-       query.xze = ddheight[i-1] + 2.*ddheighterror[i-1]
+       #query.xzs = max(ddheight[i-1] - 2.*ddheighterror[i-1],10.)
+       #query.xze = ddheight[i-1] + 2.*ddheighterror[i-1]
+       query.xzs = max(query.xz - 0.1*query.xz,10.)
+       query.xze = query.xz + 0.1*query.xz
+
        query.profile()
        tab = np.sqrt(query.zonwindtab*query.zonwindtab + query.merwindtab*query.merwindtab)
@@ -87,25 +107,40 @@
     speed = -999. ; speedmin = -1000. ; speedmax = -998. ; angle = -999. ; anglemin = -1000. ; anglemax = -998.
 
-   mcdwind = np.append(mcdwind,speed)
-   mcdwindle = np.append(mcdwindle,speed-speedmin)
-   mcdwindhe = np.append(mcdwindhe,speedmax-speed)
-   mcdangle = np.append(mcdangle,angle)
-   mcdanglele = np.append(mcdanglele,angle-anglemin)
-   mcdanglehe = np.append(mcdanglehe,anglemax-angle)
-
-mpl.figure(figsize=(12,10))
-mpl.plot([0,30], [stddev,stddev+30], 'g--')
-mpl.plot([0,30], [-stddev,-stddev+30], 'g--')
-mpl.plot([0,30], [0,30], 'g-')
-
-mpl.errorbar(mcdwind, ddwind, yerr=[ddwinderror,ddwinderror], xerr=[mcdwindle,mcdwindhe], fmt='bo')
-mpl.xlim(xmin=0.,xmax=17.)
-mpl.ylim(ymin=0.,ymax=30.)
-
-mpl.xlabel('MCD horizontal wind speed (m/s)')
-mpl.ylabel('Dust devil observed drift speed (m/s)')
-
-
-for i in yorgl[0]:
+   mcdwind.append(speed)
+   mcdwindle.append(speed-speedmin)
+   mcdwindhe.append(speedmax-speed)
+   mcdangle.append(angle)
+   mcdanglele.append(angle-anglemin)
+   mcdanglehe.append(anglemax-angle)
+
+ mcdwind = np.array(mcdwind)
+ mcdwindle = np.array(mcdwindle)
+ mcdwindhe = np.array(mcdwindhe)
+
+ mcdangle = np.array(mcdangle)
+ mcdanglele = np.array(mcdanglele)
+ mcdanglehe = np.array(mcdanglehe)
+
+
+ mpl.figure(figsize=(12,10))
+ mpl.plot([0,30], [stddev,stddev+30], 'g--')
+ mpl.plot([0,30], [-stddev,-stddev+30], 'g--')
+ mpl.plot([0,30], [0,30], 'g-')
+
+ writeascii(field=mcdwind,absc=ddwind,name=str(heights)+'.txt')
+
+
+ mpl.errorbar(mcdwind, ddwind, yerr=[ddwinderror,ddwinderror], xerr=[mcdwindle,mcdwindhe], fmt='bo')
+# mpl.xlim(xmin=0.,xmax=17.)
+# mpl.ylim(ymin=0.,ymax=30.)
+ mpl.xlim(xmin=0.,xmax=30.)
+ mpl.ylim(ymin=0.,ymax=30.)
+
+
+ mpl.xlabel('MCD horizontal wind speed (m/s)')
+ mpl.ylabel('Dust devil observed drift speed (m/s)')
+
+
+ for i in yorgl[0]:
     ## on multiplie stddev par 1.2 de facon a enlever les points tres au bord
     ##  qui sont in si on considere l error bar
@@ -118,32 +153,32 @@
 
 
-##mpl.show()
-mpl.savefig("comp.eps")
-mpl.savefig("comp.png")
-#mpl.savefig("comp4.eps")
-#mpl.savefig("comp4.png")
-
-ddwinddirerror = mcdanglele*0. + 15.
-
-mpl.figure(figsize=(12,10))
-mpl.errorbar(mcdangle,ddwinddir,yerr=[ddwinddirerror,ddwinddirerror],xerr=[mcdanglele,mcdanglehe],fmt='bo')
-mpl.plot(mcdangle,ddwinddir,'bo')
-mpl.xlim(xmin=0.,xmax=360.)
-mpl.ylim(ymin=0.,ymax=360.)
-
-mpl.xlabel('MCD horizontal wind orientation (deg)')
-mpl.ylabel('Dust devil observed drift orientation (deg)')
-
-
-stddev = 50.  ## voir e.g. 8 9 10
-mpl.plot([0,360], [stddev,stddev+360], 'g--')
-mpl.plot([0,360], [-stddev,-stddev+360], 'g--')
-mpl.plot([0,360], [0,360], 'g-')
-
-
-for i in yorgl[0]:
+ #mpl.show()
+ mpl.savefig("speed"+str(heights)+".eps")
+ mpl.savefig("speed"+str(heights)+".png")
+ ##mpl.savefig("comp4.eps")
+ ##mpl.savefig("comp4.png")
+
+ ddwinddirerror = mcdanglele*0. + 15.
+
+ mpl.figure(figsize=(12,10))
+ mpl.errorbar(mcdangle,ddwinddir,yerr=[ddwinddirerror,ddwinddirerror],xerr=[mcdanglele,mcdanglehe],fmt='bo')
+ mpl.plot(mcdangle,ddwinddir,'bo')
+ mpl.xlim(xmin=0.,xmax=360.)
+ mpl.ylim(ymin=0.,ymax=360.)
+
+ mpl.xlabel('MCD horizontal wind orientation (deg)')
+ mpl.ylabel('Dust devil observed drift orientation (deg)')
+
+
+ stddev2 = 50.  ## voir e.g. 8 9 10
+ mpl.plot([0,360], [stddev2,stddev2+360], 'g--')
+ mpl.plot([0,360], [-stddev2,-stddev2+360], 'g--')
+ mpl.plot([0,360], [0,360], 'g-')
+
+
+ for i in yorgl[0]:
     ## on multiplie stddev par 1.2 de facon a enlever les points tres au bord
     ##  qui sont in si on considere l error bar
-    if abs(ddwinddir[i-1] - mcdangle[i-1]) > stddev*1.2:
+    if abs(ddwinddir[i-1] - mcdangle[i-1]) > stddev2*1.2:
         mpl.annotate(str(int(i)), xy=(mcdangle[i-1], ddwinddir[i-1]), xytext=(10,2),
             textcoords='offset points', ha='center', va='bottom')#,
@@ -152,131 +187,130 @@
             #color='red'))
 
-##xytext=(18,30)
-
-
-mpl.savefig("comp2.eps")
-mpl.savefig("comp2.png")
-#mpl.savefig("comp3.eps")
-#mpl.savefig("comp3.png")
-
-
-exit()
-
-
-dd = mcd()
-
-
-
-
-
-#dd.zkey = 4
-#dd.xzs = 600.
-#dd.xze = 0.1
-#dd.profile()
-#dd.plot1d(["t","p"])
-#mpl.show()
+ ##xytext=(18,30)
+
+ #mpl.show()
+ mpl.savefig("angle"+str(heights)+".eps")
+ mpl.savefig("angle"+str(heights)+".png")
+ #mpl.savefig("comp3.eps")
+ #mpl.savefig("comp3.png")
+
 #exit()
-
-## CASE 1
-dd.lat   = -14.6584
-dd.lon   = 175.54
-dd.xdate = 265.278792
-dd.loct  = 14.8487
-dd.xz    = 130.4
-
-dd.update()
-dd.printmcd()
-
-dd.lats  = -5.  
-dd.late  = -25.
-dd.lons  = 160. 
-dd.lone  = 190.
-dd.map2d(["wind"],incwind=True,fixedlt=True)
-## il y a un piege pour les cartes locales il faut fixedlt=True
-mpl.savefig("case1.png",dpi=110,bbox_inches='tight',pad_inches=0.4)
-
-## CASE 2
-dd.lat   = -68.6125
-dd.lon   = 11.4303
-dd.xdate = 286.507293
-dd.loct  = 15.129030
-dd.xz    = 872.8
-#dd.xz    = 10000.0
-
-dd.update()
-dd.printmcd()
-
-dd.lats  = -75.
-dd.late  = -60.
-dd.lons  = 0.
-dd.lone  = 20.
-dd.map2d(["wind"],incwind=True,fixedlt=True)
-mpl.savefig("case2.png",dpi=110,bbox_inches='tight',pad_inches=0.4)
-
-##3	25,1306	314,4842	60,8946	14,9971	30,89575722	232	250,217961736111
-
-## CASE 3
-dd.lat   = 25.1306
-dd.lon   = 314.4842
-dd.xdate = 60.8946
-dd.loct  = 14.9971
-dd.xz    = 250.218
-
-dd.update()
-dd.printmcd()
-
-##4	35,564	199,486	60,5883	14,9949	12,2972688101353	125	895,247012611329
-
-## CASE 4
-dd.lat   = 35.564
-dd.lon   = 199.486
-dd.xdate = 60.5883
-dd.loct  = 14.9949
-dd.xz    = 895.247
-
-dd.update()
-dd.printmcd()
-
-##5	68,346	234,396	142,828	15,1777	13,4079899322222	83	1128,06581216829
-##6	68,316	234,462	142,828	15,1819	16,1939071396022	85	582,624074786015
-##7	68,302	234,144	142,828	15,1607	17,3354121022885	112	262,299040101764
-
-## CASE 5
-dd.lat   = 68.32
-dd.lon   = 234.3
-dd.xdate = 142.828
-dd.loct  = 15.2
-
-dd.xz    = 1128.066
-dd.update()
-dd.printmcd()
-
-dd.xz    = 582.624
-dd.update()
-dd.printmcd()
-
-dd.xz    = 262.299
-dd.update()
-dd.printmcd()
-
-
-##19	27,055	129,614	13,5818	14,444	42,6488205391137	302	1395,96076100437
-
-## CASE 19
-dd.lat   = 27.055
-dd.lon   = 129.614
-dd.xdate = 13.5818
-dd.loct  = 14.444
-
-dd.xz    = 1395.961
-dd.update()
-dd.printmcd()
-
-dd.lats  = 20.
-dd.late  = 30.
-dd.lons  = 125.
-dd.lone  = 135.
-dd.map2d(["wind"],incwind=True,fixedlt=True)
-mpl.savefig("case19.png",dpi=110,bbox_inches='tight',pad_inches=0.4)
-
-
+#
+#
+#dd = mcd()
+#
+#
+#
+#
+#
+##dd.zkey = 4
+##dd.xzs = 600.
+##dd.xze = 0.1
+##dd.profile()
+##dd.plot1d(["t","p"])
+##mpl.show()
+##exit()
+#
+### CASE 1
+#dd.lat   = -14.6584
+#dd.lon   = 175.54
+#dd.xdate = 265.278792
+#dd.loct  = 14.8487
+#dd.xz    = 130.4
+#
+#dd.update()
+#dd.printmcd()
+#
+#dd.lats  = -5.  
+#dd.late  = -25.
+#dd.lons  = 160. 
+#dd.lone  = 190.
+#dd.map2d(["wind"],incwind=True,fixedlt=True)
+### il y a un piege pour les cartes locales il faut fixedlt=True
+#mpl.savefig("case1.png",dpi=110,bbox_inches='tight',pad_inches=0.4)
+#
+### CASE 2
+#dd.lat   = -68.6125
+#dd.lon   = 11.4303
+#dd.xdate = 286.507293
+#dd.loct  = 15.129030
+#dd.xz    = 872.8
+##dd.xz    = 10000.0
+#
+#dd.update()
+#dd.printmcd()
+#
+#dd.lats  = -75.
+#dd.late  = -60.
+#dd.lons  = 0.
+#dd.lone  = 20.
+#dd.map2d(["wind"],incwind=True,fixedlt=True)
+#mpl.savefig("case2.png",dpi=110,bbox_inches='tight',pad_inches=0.4)
+#
+###3	25,1306	314,4842	60,8946	14,9971	30,89575722	232	250,217961736111
+#
+### CASE 3
+#dd.lat   = 25.1306
+#dd.lon   = 314.4842
+#dd.xdate = 60.8946
+#dd.loct  = 14.9971
+#dd.xz    = 250.218
+#
+#dd.update()
+#dd.printmcd()
+#
+###4	35,564	199,486	60,5883	14,9949	12,2972688101353	125	895,247012611329
+#
+### CASE 4
+#dd.lat   = 35.564
+#dd.lon   = 199.486
+#dd.xdate = 60.5883
+#dd.loct  = 14.9949
+#dd.xz    = 895.247
+#
+#dd.update()
+#dd.printmcd()
+#
+###5	68,346	234,396	142,828	15,1777	13,4079899322222	83	1128,06581216829
+###6	68,316	234,462	142,828	15,1819	16,1939071396022	85	582,624074786015
+###7	68,302	234,144	142,828	15,1607	17,3354121022885	112	262,299040101764
+#
+### CASE 5
+#dd.lat   = 68.32
+#dd.lon   = 234.3
+#dd.xdate = 142.828
+#dd.loct  = 15.2
+#
+#dd.xz    = 1128.066
+#dd.update()
+#dd.printmcd()
+#
+#dd.xz    = 582.624
+#dd.update()
+#dd.printmcd()
+#
+#dd.xz    = 262.299
+#dd.update()
+#dd.printmcd()
+#
+#
+###19	27,055	129,614	13,5818	14,444	42,6488205391137	302	1395,96076100437
+#
+### CASE 19
+#dd.lat   = 27.055
+#dd.lon   = 129.614
+#dd.xdate = 13.5818
+#dd.loct  = 14.444
+#
+#dd.xz    = 1395.961
+#dd.update()
+#dd.printmcd()
+#
+#dd.lats  = 20.
+#dd.late  = 30.
+#dd.lons  = 125.
+#dd.lone  = 135.
+#dd.map2d(["wind"],incwind=True,fixedlt=True)
+#mpl.savefig("case19.png",dpi=110,bbox_inches='tight',pad_inches=0.4)
+#
+#
Index: trunk/UTIL/PYTHON/mcd/examples/gale.py
===================================================================
--- trunk/UTIL/PYTHON/mcd/examples/gale.py	(revision 971)
+++ trunk/UTIL/PYTHON/mcd/examples/gale.py	(revision 1007)
@@ -35,4 +35,8 @@
 mpl.savefig("temp.png",dpi=85,bbox_inches='tight',pad_inches=0.25)
 
+exit()
+
+
+
 
 #gale.xdate = 270.
Index: trunk/UTIL/PYTHON/mcd/inimeso/inimeso.py
===================================================================
--- trunk/UTIL/PYTHON/mcd/inimeso/inimeso.py	(revision 971)
+++ trunk/UTIL/PYTHON/mcd/inimeso/inimeso.py	(revision 1007)
@@ -37,2 +37,7 @@
 sounding.close() ; additional.close()
 query.plot1d(["p","t","u","v"]) ; mpl.show()
+
+### Add information about dust opacity
+dod = open("dustopacity.def", "w")
+dod.write( "%4.2f"  %(query.extvartab[iz,36]) )
+dod.close()
