Index: trunk/UTIL/PYTHON/mcd/exampleascii.py
===================================================================
--- trunk/UTIL/PYTHON/mcd/exampleascii.py	(revision 937)
+++ 	(revision )
@@ -1,14 +1,0 @@
-#! /usr/bin/python
-
-from mcd import mcd
-
-
-yeah = mcd()
-
-yeah.diurnal()
-yeah.getascii(["u","t"],filename="diurnal.txt")
-
-
-yeah.zonal()
-yeah.getascii(["u","t"],filename="zonal.txt")
-
Index: trunk/UTIL/PYTHON/mcd/examples/dustdevil.py
===================================================================
--- trunk/UTIL/PYTHON/mcd/examples/dustdevil.py	(revision 943)
+++ trunk/UTIL/PYTHON/mcd/examples/dustdevil.py	(revision 943)
@@ -0,0 +1,282 @@
+#! /usr/bin/env python
+
+from mcd import mcd
+from myplot import maplatlon
+import numpy as np
+import matplotlib.pyplot as mpl
+
+from math import isnan
+
+dafile = "/home/aymeric/Work/submitted/coauthorship/2012reiss/reiss_data.txt"
+
+yorgl = np.loadtxt(dafile)
+yorgl = np.transpose(yorgl)
+ddwind = yorgl[5]
+ddwinddir = yorgl[6]
+ddwinderror = yorgl[9]
+
+ddheight = yorgl[7]
+ddheighterror = yorgl[8]
+
+#ddheight = yorgl[7]*0. + 10.
+#ddheighterror = yorgl[8]*0. + 0.5
+
+stddev = np.std(ddwind[np.where(ddwind != -999.)])
+print stddev
+
+## missing points are NaN
+ind = np.where(ddwind == -999.) ; ddwind[ind] = np.NaN
+print "1. wind not available", yorgl[0][ind]
+ind = np.where(ddheight == -999.) ; ddwind[ind] = np.NaN
+print "1. height not available", yorgl[0][ind]
+## no data point with rel error more than x%
+#ind = np.where(100. * ddwinderror / ddwinddir > 5.)
+#there was an error here!
+tol = 50.
+ind = np.where(100. * ddwinderror / ddwind > tol) ; ddwind[ind] = np.NaN 
+print "2. error on wind too large"
+print yorgl[0][ind]
+print ddwinderror[ind]
+ind = np.where(100. * ddheighterror / ddheight > tol) ; ddwind[ind] = np.NaN 
+print "3. error on height too large"
+print yorgl[0][ind]
+print ddheighterror[ind]
+## one record is doubtful according to Reiss
+ind = np.where( yorgl[0] == 3 ) ; ddwind[ind] = np.NaN #; ddwinddir[ind] = np.NaN ; ddwinderror[ind] = np.NaN
+print "4. doubtful", yorgl[0][ind]
+
+#exit()
+
+query = mcd()
+mcdwind = [] ; mcdwindle = [] ; mcdwindhe = [] ; mcdangle = [] ; mcdanglele = [] ; mcdanglehe = [] 
+
+for i in yorgl[0]:
+
+   if not isnan(ddwind[i-1]):
+    query.lat = yorgl[1][i-1]
+    query.lon = yorgl[2][i-1]
+    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.update() ; uu = query.zonwind ; vv = query.merwind
+    if uu == -999. or vv == -999.: 
+       speed = np.NaN ; speedmin = np.NaN ; speedmax = np.NaN
+       angle = np.NaN ; anglemin = np.NaN ; anglemax = np.NaN
+    else:
+       speed = np.sqrt(uu*uu+vv*vv) ; speedmin = speed ; speedmax = speed
+
+       angle = 90. - np.arctan2(vv,uu) * 180. / np.pi
+       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.profile()
+       tab = np.sqrt(query.zonwindtab*query.zonwindtab + query.merwindtab*query.merwindtab)
+       tab2 = 90. - np.arctan2(query.merwindtab,query.zonwindtab) * 180. / np.pi
+       tab2[np.where(tab2 < 0.)] = tab2[np.where(tab2 < 0.)] + 360.
+       speedmin = min(tab)
+       speedmax = max(tab)
+       anglemin = min(tab2)-5.
+       anglemax = max(tab2)+5.
+       if anglemin < 0.: anglemin = angle - (anglemax - angle)
+       print int(i), speed, speedmin, speedmax 
+       print int(i), angle, anglemin, anglemax
+   else:
+    #print "removed ",i,ddwinderror[i-1],ddheighterror[i-1]
+    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]:
+    ## 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(ddwind[i-1] - mcdwind[i-1]) > stddev*1.2:
+        mpl.annotate(str(int(i)), xy=(mcdwind[i-1], ddwind[i-1]), xytext=(10,2), 
+            textcoords='offset points', ha='center', va='bottom' )#,
+            #bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.3),
+            #arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=-0.3', 
+            #color='red'))
+
+
+##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]:
+    ## 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:
+        mpl.annotate(str(int(i)), xy=(mcdangle[i-1], ddwinddir[i-1]), xytext=(10,2),
+            textcoords='offset points', ha='center', va='bottom')#,
+            #bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.3),
+            #arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=-0.3',
+            #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()
+#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/exampleascii.py
===================================================================
--- trunk/UTIL/PYTHON/mcd/examples/exampleascii.py	(revision 943)
+++ trunk/UTIL/PYTHON/mcd/examples/exampleascii.py	(revision 943)
@@ -0,0 +1,14 @@
+#! /usr/bin/python
+
+from mcd import mcd
+
+
+yeah = mcd()
+
+yeah.diurnal()
+yeah.getascii(["u","t"],filename="diurnal.txt")
+
+
+yeah.zonal()
+yeah.getascii(["u","t"],filename="zonal.txt")
+
Index: trunk/UTIL/PYTHON/mcd/examples/gale.py
===================================================================
--- trunk/UTIL/PYTHON/mcd/examples/gale.py	(revision 943)
+++ trunk/UTIL/PYTHON/mcd/examples/gale.py	(revision 943)
@@ -0,0 +1,47 @@
+#! /usr/bin/env python
+
+#import sys
+#
+#class NullDevice():
+#    def write(self, s):
+#        pass
+#
+#original_stdout = sys.stdout  # keep a reference to STDOUT
+#
+#sys.stdout = NullDevice()  # redirect the real STDOUT
+#sys.stderr = NullDevice()  # redirect the real STDOUT
+
+from mcd import mcd
+
+gale = mcd()
+
+
+gale.lat = -4.6
+gale.lon = 137.4
+gale.loct = 15.
+gale.xz = 1.
+
+gale.xdate = 150.6
+
+gale.update()
+gale.printmcd()
+gale.printallextvar()
+#gale.printextvar("rho")
+
+gale.seasonal()
+gale.plot1d(["tsurf","u","v"])
+
+import matplotlib.pyplot as mpl
+mpl.savefig("temp.png",dpi=85,bbox_inches='tight',pad_inches=0.25)
+
+
+#gale.xdate = 270.
+#gale.diurnal()
+#gale.plot1d(["u","v"])
+
+gale.latlon()
+gale.map2d("tsurf")
+
+import matplotlib.pyplot as mpl
+#mpl.show()
+mpl.savefig("temp.png",dpi=110,bbox_inches='tight',pad_inches=0.4)
Index: trunk/UTIL/PYTHON/mcd/examples/mcdtest.ipython
===================================================================
--- trunk/UTIL/PYTHON/mcd/examples/mcdtest.ipython	(revision 943)
+++ trunk/UTIL/PYTHON/mcd/examples/mcdtest.ipython	(revision 943)
@@ -0,0 +1,7 @@
+example of use
+--------------
+
+from mcd import mcd
+myrequest = mcd()
+myrequest.latlon()
+myrequest.map2d(22)
Index: trunk/UTIL/PYTHON/mcd/examples/test_mcd.py
===================================================================
--- trunk/UTIL/PYTHON/mcd/examples/test_mcd.py	(revision 943)
+++ trunk/UTIL/PYTHON/mcd/examples/test_mcd.py	(revision 943)
@@ -0,0 +1,99 @@
+#! /usr/bin/env python
+
+# This python program illustrates how the Mars Climate Database
+# fortran routines can be called interactively from python
+
+from fmcd import call_mcd,julian
+import numpy as np
+
+# 1. Inputs:
+dset = '' # default to 'MCD_DATA'
+perturkey = 1 # default to no perturbation
+seedin = 0 # perturbation seed (unused if perturkey=1)
+gwlength = 0 # Gravity Wave length for perturbations (unused if perturkey=1)
+
+# 1.1 Dates
+choice_date=raw_input('Use Earth date (e) or Mars date (m)?')
+if (choice_date == "e") :
+  datekey=0 # Earth date
+  loct=0 # local time must then also be set to zero
+  day,month,year,hour,minute,second = \
+  raw_input("Enter date: day/month/year/hour/minute/second: ").split('/')
+  day=int(day)
+  month=int(month)
+  year=int(year)
+  hour=int(hour)
+  minute=int(minute)
+  second=int(second)
+  # now call Julian routine to convert to julian date
+  (ier,xdate)=julian(month,day,year,hour,minute,second)
+  print " Julian date %16.8f" % xdate
+else :
+  datekey=1 # Mars date
+  xdate=float(raw_input("Enter solar longitude Ls (deg.):"))
+  loct=float(raw_input("Local time (0 < time < 24)?"))
+
+# 1.2 Vertical coordinate
+zkey=int(raw_input("Select verical coordinate type (1: distance to center of planet, 2: height above areoid, 3: height above surface, 4: Pressure) "))
+if (zkey == 1) :
+  xz=float(raw_input("Enter distance to planet center (m) "))
+if (zkey == 2) :
+  xz=float(raw_input("Enter altitude above areoid (m) "))
+if (zkey == 3) :
+  xz=float(raw_input("Enter altitude above surface (m) "))
+if (zkey == 4) :
+  xz=float(raw_input("Enter pressure value (Pa) "))
+
+# high resolution mode
+hrkey=int(raw_input("High resolution? (1: yes, 0: no) "))
+
+# 1.3 Position
+lat = float(raw_input('Latitude (deg)?'))
+lon = float(raw_input('Longitude (deg)?'))
+
+# 1.4 Dust and solar scenario
+print "Dust scenario?"
+print "1= Climatology       typical Mars year dust scenario"
+print "                     average solar EUV conditions"
+print "2= Climatology       typical Mars year dust scenario"
+print "                     minimum solar EUV conditions"
+print "3= Climatology       typical Mars year dust scenario"
+print "                     maximum solar EUV conditions"
+print "4= dust storm        constant dust opacity = 5 (dark dust)"
+print "                     minimum solar EUV conditions"
+print "5= dust storm        constant dust opacity = 5 (dark dust)"
+print "                     average solar EUV conditions"
+print "6= dust storm        constant dust opacity = 5 (dark dust)"
+print "                     maximum solar EUV conditions"
+print "7= warm scenario     dustier than Climatology scenario"
+print "                     maximum solar EUV conditions"
+print "8= cold scenario     clearer than Climatology scenario"
+print "                     minimum solar EUV conditions"
+dust=int(raw_input(''))
+
+# 1.5 perturbations
+perturkey=int(raw_input("Perturbation? (1:none, 2: large scale, 3: small scale, 4: small+large, 5: n sigmas) "))
+if (perturkey > 1) :
+  seedin=int(raw_input("seedin? (only matters if adding perturbations) "))
+if ((perturkey == 3) or (perturkey == 4)) :
+  gwlength=float(raw_input("Gravity wave length? (for small scale perturbations) "))
+
+# 1.6 extra outputs
+# here we only implement an all-or-nothing case
+extvarkey=int(raw_input("Output the extra variables? (yes==1; no==0) "))
+if (extvarkey == 0) :
+  extvarkeys = np.zeros(100)
+else :
+  extvarkeys = np.ones(100)
+
+# 2. Call MCD
+(pres, dens, temp, zonwind, merwind, \
+ meanvar, extvar, seedout, ierr) \
+ = \
+call_mcd(zkey,xz,lon,lat,hrkey, \
+ datekey,xdate,loct,dset,dust, \
+ perturkey,seedin,gwlength,extvarkeys )
+
+# 3. Write outputs
+print "temperature is %.0f K, pressure is %.0f Pa, density is %5.3e kg/m3, zonal wind is %.1f m/s, meridional wind is %.1f m/s" % (temp,pres,dens,zonwind,merwind)
+
Index: trunk/UTIL/PYTHON/mcd/examples/test_mcd_alt.py
===================================================================
--- trunk/UTIL/PYTHON/mcd/examples/test_mcd_alt.py	(revision 943)
+++ trunk/UTIL/PYTHON/mcd/examples/test_mcd_alt.py	(revision 943)
@@ -0,0 +1,32 @@
+#! /usr/bin/env python
+
+from fmcd import call_mcd
+import numpy as np
+
+lon = float(raw_input('Longitude?'))
+lat = float(raw_input('Latitude?'))
+xdate = float(raw_input('Ls?'))
+loct = float(raw_input('Local time?'))
+xz = float(raw_input('Altitude?'))
+
+dset = '/home/aymeric/Science/MCD_v4.3/data/'
+
+zkey      = 3  # specify that xz is the altitude above surface (m)
+datekey   = 1  # 1 = "Mars date": xdate is the value of Ls
+dust      = 2  #our best guess MY24 scenario, with solar average conditions
+hrkey     = 1  #set high resolution mode on (hrkey=0 to set high resolution off)
+perturkey = 0  #integer perturkey ! perturbation type (0: none)
+seedin    = 1  #random number generator seed (unused if perturkey=0)
+gwlength  = 0. #gravity Wave wavelength (unused if perturkey=0)
+extvarkey = 1
+#extvarkeys = np.ones(100)
+
+(pres, dens, temp, zonwind, merwind, \
+ meanvar, extvar, seedout, ierr) \
+ = \
+call_mcd(zkey,xz,lon,lat,hrkey, \
+ datekey,xdate,loct,dset,dust, \
+ perturkey,seedin,gwlength,extvarkey )
+
+print "temperature is %.0f K, pressure is %.0f Pa, density is %5.3e kg/m3, zonal wind is %.1f m/s, meridional wind is %.1f m/s" % (temp,pres,dens,zonwind,merwind)
+
Index: trunk/UTIL/PYTHON/mcd/fmcd5_gfortran.sh
===================================================================
--- trunk/UTIL/PYTHON/mcd/fmcd5_gfortran.sh	(revision 943)
+++ trunk/UTIL/PYTHON/mcd/fmcd5_gfortran.sh	(revision 943)
@@ -0,0 +1,57 @@
+#! /bin/bash
+
+###
+### fmcd_gfortran.sh
+###
+### ---> This should be used to make MCD fortran stuff directly accessible in python 
+### ---> This script is for gfortran, but it is easy to adapt to your other compilers
+### ---> A file fmcd.so should be created
+### ---> See mcd.py for use in python. Very easy!
+###
+### AS. 17/04/2012. 
+###
+
+### LINKS
+NETCDF=/home/marshttp/NETCDF/netcdf64-4.0.1_gfortran_fPIC/
+wheremcd=/home/marshttp/MCD_v5.0//mcd/
+
+### LOG FILE
+touch fmcd5.log
+\rm fmcd5.log
+
+### COPY/PREPARE SOURCES
+### perform changes that makes f2py not to fail
+sed s/"\!\!'"/"'"/g $wheremcd/call_mcd.F         | sed s/"\!'"/"'"/g | sed s/"\!"/"\n\!"/g > tmp.call_mcd.F
+sed s/"\!\!'"/"'"/g $wheremcd/julian.F           | sed s/"\!'"/"'"/g | sed s/"\!"/"\n\!"/g > tmp.julian.F
+sed s/"\!\!'"/"'"/g $wheremcd/heights.F          | sed s/"\!'"/"'"/g | sed s/"\!"/"\n\!"/g > tmp.heights.F
+sed s/"\!\!'"/"'"/g $wheremcd/constants_mcd.inc  | sed s/"\!'"/"'"/g | sed s/"\!"/"\n\!"/g > constants_mcd.inc
+
+### BUILD THROUGH f2py WHAT IS NECESSARY TO CREATE THE PYTHON FUNCTIONS
+touch fmcd5.pyf
+\rm fmcd5.pyf
+f2py -h fmcd5.pyf -m fmcd5 tmp.call_mcd.F tmp.julian.F tmp.heights.F > fmcd5.log 2>&1
+
+#### IMPORTANT: we teach f2py about variables in the call_mcd subroutines which are intended to be out
+sed s/"real :: pres"/"real, intent(out) :: pres"/g fmcd5.pyf | \
+sed s/"real :: dens"/"real, intent(out) :: dens"/g | \
+sed s/"real :: temp"/"real, intent(out) :: temp"/g | \
+sed s/"real :: zonwind"/"real, intent(out) :: zonwind"/g | \
+sed s/"real :: merwind"/"real, intent(out) :: merwind"/g | \
+sed s/"real dimension(5) :: meanvar"/"real dimension(5),intent(out) :: meanvar"/g | \
+sed s/"real dimension(100) :: extvar"/"real dimension(100),intent(out) :: extvar"/g | \
+sed s/"real :: seedout"/"real, intent(out) :: seedout"/g | \
+sed s/"integer :: ier"/"integer, intent(out) :: ier"/g > fmcd5.pyf.modif
+mv fmcd5.pyf.modif fmcd5.pyf
+
+### BUILD
+f2py -c fmcd5.pyf tmp.call_mcd.F tmp.julian.F tmp.heights.F --fcompiler=gnu95 \
+  -L$NETCDF/lib -lnetcdf \
+  -lm -I$NETCDF/include \
+  --f90flags="-fPIC" \
+  --f77flags="-fPIC" \
+  --verbose \
+  > fmcd5.log 2>&1
+#  --include-paths $NETCDF/include:$NETCDF/lib \ ---> makes it fail
+
+### CLEAN THE PLACE
+\rm tmp.call_mcd.F tmp.julian.F tmp.heights.F constants_mcd.inc
Index: trunk/UTIL/PYTHON/mcd/gale.py
===================================================================
--- trunk/UTIL/PYTHON/mcd/gale.py	(revision 937)
+++ 	(revision )
@@ -1,47 +1,0 @@
-#! /usr/bin/env python
-
-#import sys
-#
-#class NullDevice():
-#    def write(self, s):
-#        pass
-#
-#original_stdout = sys.stdout  # keep a reference to STDOUT
-#
-#sys.stdout = NullDevice()  # redirect the real STDOUT
-#sys.stderr = NullDevice()  # redirect the real STDOUT
-
-from mcd import mcd
-
-gale = mcd()
-
-
-gale.lat = -4.6
-gale.lon = 137.4
-gale.loct = 15.
-gale.xz = 1.
-
-gale.xdate = 150.6
-
-gale.update()
-gale.printmcd()
-gale.printallextvar()
-#gale.printextvar("rho")
-
-gale.seasonal()
-gale.plot1d(["tsurf","u","v"])
-
-import matplotlib.pyplot as mpl
-mpl.savefig("temp.png",dpi=85,bbox_inches='tight',pad_inches=0.25)
-
-
-#gale.xdate = 270.
-#gale.diurnal()
-#gale.plot1d(["u","v"])
-
-gale.latlon()
-gale.map2d("tsurf")
-
-import matplotlib.pyplot as mpl
-#mpl.show()
-mpl.savefig("temp.png",dpi=110,bbox_inches='tight',pad_inches=0.4)
Index: trunk/UTIL/PYTHON/mcd/inimeso.py
===================================================================
--- trunk/UTIL/PYTHON/mcd/inimeso.py	(revision 937)
+++ 	(revision )
@@ -1,38 +1,0 @@
-#! /usr/bin/env python
-
-### AS 10/05/2012. A python script to prepare initial state for idealized mesoscale runs.
-###                use : ensure mcd class is working. fill in input_coord. execute inimeso.
-
-from string import split ; import numpy as np ; import matplotlib.pyplot as mpl
-from mcd import mcd
-
-### MCD INSTANCE and SETTINGS (actually, default. but one never knows)
-query = mcd() ; query.zkey = 3 ; query.dust = 2 ; query.hrkey = 1
-
-### GET COORDINATES
-lines = open("input_coord", "r").readlines()
-query.lon   = float(split(lines[0])[0]) ; query.lat   = float(split(lines[1])[0]) 
-query.xdate = float(split(lines[2])[0]) ; query.loct  = float(split(lines[3])[0])
-query.printcoord()
-
-### OPEN FILES TO BE WRITTEN
-sounding = open("input_sounding", "w") ; additional = open("input_therm", "w") ; more = open("input_more", "w")
-
-### GET and WRITE SURFACE VALUES
-query.xz = 0. ; query.update() ; query.printmeanvar()
-sounding.write( "%10.2f%12.2f%12.2f\n" % (query.pres/100.,query.temp*(610./query.pres)**(1.0/3.9),0.) )
-more.write( "%10.2f%10.2f" % (query.extvar[1],query.extvar[14]) ) ; more.close()
-
-### GET and WRITE VERTICAL PROFILE
-query.profile( tabperso = np.append([0,1,5,10,20,50,100],np.linspace(200.,float(split(lines[4])[0])*1000.,float(split(lines[5])[0]))) )
-for iz in range(len(query.prestab)):
-    sounding.write(   "%10.2f%12.2f%12.2f%12.2f%12.2f\n" % ( \
-                      query.extvartab[iz,2],query.temptab[iz]*(610./query.prestab[iz])**(1.0/3.9),\
-                      0.,query.zonwindtab[iz],query.merwindtab[iz]) )
-    additional.write( "%12.2f%12.2f%18.6e%18.6e%12.2f\n" % ( \
-                      query.extvartab[iz,49],query.extvartab[iz,8],\
-                      query.prestab[iz],query.denstab[iz],query.temptab[iz]) )
-
-### FINISH
-sounding.close() ; additional.close()
-query.plot1d(["p","t","u","v"]) ; mpl.show()
Index: trunk/UTIL/PYTHON/mcd/inimeso/inimeso.py
===================================================================
--- trunk/UTIL/PYTHON/mcd/inimeso/inimeso.py	(revision 943)
+++ trunk/UTIL/PYTHON/mcd/inimeso/inimeso.py	(revision 943)
@@ -0,0 +1,38 @@
+#! /usr/bin/env python
+
+### AS 10/05/2012. A python script to prepare initial state for idealized mesoscale runs.
+###                use : ensure mcd class is working. fill in input_coord. execute inimeso.
+
+from string import split ; import numpy as np ; import matplotlib.pyplot as mpl
+from mcd import mcd
+
+### MCD INSTANCE and SETTINGS (actually, default. but one never knows)
+query = mcd() ; query.zkey = 3 ; query.dust = 2 ; query.hrkey = 1
+
+### GET COORDINATES
+lines = open("input_coord", "r").readlines()
+query.lon   = float(split(lines[0])[0]) ; query.lat   = float(split(lines[1])[0]) 
+query.xdate = float(split(lines[2])[0]) ; query.loct  = float(split(lines[3])[0])
+query.printcoord()
+
+### OPEN FILES TO BE WRITTEN
+sounding = open("input_sounding", "w") ; additional = open("input_therm", "w") ; more = open("input_more", "w")
+
+### GET and WRITE SURFACE VALUES
+query.xz = 0. ; query.update() ; query.printmeanvar()
+sounding.write( "%10.2f%12.2f%12.2f\n" % (query.pres/100.,query.temp*(610./query.pres)**(1.0/3.9),0.) )
+more.write( "%10.2f%10.2f" % (query.extvar[1],query.extvar[14]) ) ; more.close()
+
+### GET and WRITE VERTICAL PROFILE
+query.profile( tabperso = np.append([0,1,5,10,20,50,100],np.linspace(200.,float(split(lines[4])[0])*1000.,float(split(lines[5])[0]))) )
+for iz in range(len(query.prestab)):
+    sounding.write(   "%10.2f%12.2f%12.2f%12.2f%12.2f\n" % ( \
+                      query.extvartab[iz,2],query.temptab[iz]*(610./query.prestab[iz])**(1.0/3.9),\
+                      0.,query.zonwindtab[iz],query.merwindtab[iz]) )
+    additional.write( "%12.2f%12.2f%18.6e%18.6e%12.2f\n" % ( \
+                      query.extvartab[iz,49],query.extvartab[iz,8],\
+                      query.prestab[iz],query.denstab[iz],query.temptab[iz]) )
+
+### FINISH
+sounding.close() ; additional.close()
+query.plot1d(["p","t","u","v"]) ; mpl.show()
Index: trunk/UTIL/PYTHON/mcd/inimeso/inimeso5.py
===================================================================
--- trunk/UTIL/PYTHON/mcd/inimeso/inimeso5.py	(revision 943)
+++ trunk/UTIL/PYTHON/mcd/inimeso/inimeso5.py	(revision 943)
@@ -0,0 +1,40 @@
+#! /usr/bin/env python
+
+### AS 10/05/2012. A python script to prepare initial state for idealized mesoscale runs.
+###                use : ensure mcd class is working. fill in input_coord. execute inimeso.
+
+from string import split ; import numpy as np ; import matplotlib.pyplot as mpl
+from mcd import mcd
+
+### MCD INSTANCE and SETTINGS (actually, default. but one never knows)
+query = mcd() ; query.zkey = 3 ; query.dust = 2 ; query.hrkey = 1
+
+query.toversion5()
+
+### GET COORDINATES
+lines = open("input_coord", "r").readlines()
+query.lon   = float(split(lines[0])[0]) ; query.lat   = float(split(lines[1])[0]) 
+query.xdate = float(split(lines[2])[0]) ; query.loct  = float(split(lines[3])[0])
+query.printcoord()
+
+### OPEN FILES TO BE WRITTEN
+sounding = open("input_sounding", "w") ; additional = open("input_therm", "w") ; more = open("input_more", "w")
+
+### GET and WRITE SURFACE VALUES
+query.xz = 0. ; query.update() ; query.printmeanvar()
+sounding.write( "%10.2f%12.2f%12.2f\n" % (query.pres/100.,query.temp*(610./query.pres)**(1.0/3.9),0.) )
+more.write( "%10.2f%10.2f" % (query.extvar[1],query.extvar[14]) ) ; more.close()
+
+### GET and WRITE VERTICAL PROFILE
+query.profile( tabperso = np.append([0,5,10,20,50,100],np.linspace(200.,float(split(lines[4])[0])*1000.,float(split(lines[5])[0]))) )
+for iz in range(len(query.prestab)):
+    sounding.write(   "%10.2f%12.2f%12.2f%12.2f%12.2f\n" % ( \
+                      query.extvartab[iz,2],query.temptab[iz]*(610./query.prestab[iz])**(1.0/3.9),\
+                      0.,query.zonwindtab[iz],query.merwindtab[iz]) )
+    additional.write( "%12.2f%12.2f%18.6e%18.6e%12.2f\n" % ( \
+                      query.extvartab[iz,53],query.extvartab[iz,8],\
+                      query.prestab[iz],query.denstab[iz],query.temptab[iz]) )
+
+### FINISH
+sounding.close() ; additional.close()
+query.plot1d(["p","t","u","v"]) ; mpl.show()
Index: trunk/UTIL/PYTHON/mcd/inimeso/input_coord
===================================================================
--- trunk/UTIL/PYTHON/mcd/inimeso/input_coord	(revision 943)
+++ trunk/UTIL/PYTHON/mcd/inimeso/input_coord	(revision 943)
@@ -0,0 +1,6 @@
+353.87  lon
+-1.88   lat
+244.    ls
+06.     lct
+20      ztop
+100.    nz
Index: trunk/UTIL/PYTHON/mcd/inimeso5.py
===================================================================
--- trunk/UTIL/PYTHON/mcd/inimeso5.py	(revision 937)
+++ 	(revision )
@@ -1,40 +1,0 @@
-#! /usr/bin/env python
-
-### AS 10/05/2012. A python script to prepare initial state for idealized mesoscale runs.
-###                use : ensure mcd class is working. fill in input_coord. execute inimeso.
-
-from string import split ; import numpy as np ; import matplotlib.pyplot as mpl
-from mcd import mcd
-
-### MCD INSTANCE and SETTINGS (actually, default. but one never knows)
-query = mcd() ; query.zkey = 3 ; query.dust = 2 ; query.hrkey = 1
-
-query.toversion5()
-
-### GET COORDINATES
-lines = open("input_coord", "r").readlines()
-query.lon   = float(split(lines[0])[0]) ; query.lat   = float(split(lines[1])[0]) 
-query.xdate = float(split(lines[2])[0]) ; query.loct  = float(split(lines[3])[0])
-query.printcoord()
-
-### OPEN FILES TO BE WRITTEN
-sounding = open("input_sounding", "w") ; additional = open("input_therm", "w") ; more = open("input_more", "w")
-
-### GET and WRITE SURFACE VALUES
-query.xz = 0. ; query.update() ; query.printmeanvar()
-sounding.write( "%10.2f%12.2f%12.2f\n" % (query.pres/100.,query.temp*(610./query.pres)**(1.0/3.9),0.) )
-more.write( "%10.2f%10.2f" % (query.extvar[1],query.extvar[14]) ) ; more.close()
-
-### GET and WRITE VERTICAL PROFILE
-query.profile( tabperso = np.append([0,5,10,20,50,100],np.linspace(200.,float(split(lines[4])[0])*1000.,float(split(lines[5])[0]))) )
-for iz in range(len(query.prestab)):
-    sounding.write(   "%10.2f%12.2f%12.2f%12.2f%12.2f\n" % ( \
-                      query.extvartab[iz,2],query.temptab[iz]*(610./query.prestab[iz])**(1.0/3.9),\
-                      0.,query.zonwindtab[iz],query.merwindtab[iz]) )
-    additional.write( "%12.2f%12.2f%18.6e%18.6e%12.2f\n" % ( \
-                      query.extvartab[iz,53],query.extvartab[iz,8],\
-                      query.prestab[iz],query.denstab[iz],query.temptab[iz]) )
-
-### FINISH
-sounding.close() ; additional.close()
-query.plot1d(["p","t","u","v"]) ; mpl.show()
Index: trunk/UTIL/PYTHON/mcd/input_coord
===================================================================
--- trunk/UTIL/PYTHON/mcd/input_coord	(revision 937)
+++ 	(revision )
@@ -1,6 +1,0 @@
-353.87  lon
--1.88   lat
-244.    ls
-06.     lct
-20      ztop
-100.    nz
Index: trunk/UTIL/PYTHON/mcd/mcd.py
===================================================================
--- trunk/UTIL/PYTHON/mcd/mcd.py	(revision 937)
+++ trunk/UTIL/PYTHON/mcd/mcd.py	(revision 943)
@@ -28,6 +28,6 @@
         self.name      = "MCD v4.3"
         self.ack       = "Mars Climate Database (c) LMD/OU/IAA/ESA/CNES"
-        #self.dset      = '/home/aymeric/Science/MCD_v4.3/data/'
-        self.dset      = '/home/marshttp/MCD_v4.3/data/'
+        self.dset      = '/home/aymeric/Science/MCD_v4.3/data/'
+        #self.dset      = '/home/marshttp/MCD_v4.3/data/'
         ## 1. spatio-temporal coordinates
         self.lat       = 0.
Index: trunk/UTIL/PYTHON/mcd/mcdtest.ipython
===================================================================
--- trunk/UTIL/PYTHON/mcd/mcdtest.ipython	(revision 937)
+++ 	(revision )
@@ -1,7 +1,0 @@
-example of use
---------------
-
-from mcd import mcd
-myrequest = mcd()
-myrequest.latlon()
-myrequest.map2d(22)
Index: trunk/UTIL/PYTHON/mcd/test_mcd.py
===================================================================
--- trunk/UTIL/PYTHON/mcd/test_mcd.py	(revision 937)
+++ 	(revision )
@@ -1,99 +1,0 @@
-#! /usr/bin/env python
-
-# This python program illustrates how the Mars Climate Database
-# fortran routines can be called interactively from python
-
-from fmcd import call_mcd,julian
-import numpy as np
-
-# 1. Inputs:
-dset = '' # default to 'MCD_DATA'
-perturkey = 1 # default to no perturbation
-seedin = 0 # perturbation seed (unused if perturkey=1)
-gwlength = 0 # Gravity Wave length for perturbations (unused if perturkey=1)
-
-# 1.1 Dates
-choice_date=raw_input('Use Earth date (e) or Mars date (m)?')
-if (choice_date == "e") :
-  datekey=0 # Earth date
-  loct=0 # local time must then also be set to zero
-  day,month,year,hour,minute,second = \
-  raw_input("Enter date: day/month/year/hour/minute/second: ").split('/')
-  day=int(day)
-  month=int(month)
-  year=int(year)
-  hour=int(hour)
-  minute=int(minute)
-  second=int(second)
-  # now call Julian routine to convert to julian date
-  (ier,xdate)=julian(month,day,year,hour,minute,second)
-  print " Julian date %16.8f" % xdate
-else :
-  datekey=1 # Mars date
-  xdate=float(raw_input("Enter solar longitude Ls (deg.):"))
-  loct=float(raw_input("Local time (0 < time < 24)?"))
-
-# 1.2 Vertical coordinate
-zkey=int(raw_input("Select verical coordinate type (1: distance to center of planet, 2: height above areoid, 3: height above surface, 4: Pressure) "))
-if (zkey == 1) :
-  xz=float(raw_input("Enter distance to planet center (m) "))
-if (zkey == 2) :
-  xz=float(raw_input("Enter altitude above areoid (m) "))
-if (zkey == 3) :
-  xz=float(raw_input("Enter altitude above surface (m) "))
-if (zkey == 4) :
-  xz=float(raw_input("Enter pressure value (Pa) "))
-
-# high resolution mode
-hrkey=int(raw_input("High resolution? (1: yes, 0: no) "))
-
-# 1.3 Position
-lat = float(raw_input('Latitude (deg)?'))
-lon = float(raw_input('Longitude (deg)?'))
-
-# 1.4 Dust and solar scenario
-print "Dust scenario?"
-print "1= Climatology       typical Mars year dust scenario"
-print "                     average solar EUV conditions"
-print "2= Climatology       typical Mars year dust scenario"
-print "                     minimum solar EUV conditions"
-print "3= Climatology       typical Mars year dust scenario"
-print "                     maximum solar EUV conditions"
-print "4= dust storm        constant dust opacity = 5 (dark dust)"
-print "                     minimum solar EUV conditions"
-print "5= dust storm        constant dust opacity = 5 (dark dust)"
-print "                     average solar EUV conditions"
-print "6= dust storm        constant dust opacity = 5 (dark dust)"
-print "                     maximum solar EUV conditions"
-print "7= warm scenario     dustier than Climatology scenario"
-print "                     maximum solar EUV conditions"
-print "8= cold scenario     clearer than Climatology scenario"
-print "                     minimum solar EUV conditions"
-dust=int(raw_input(''))
-
-# 1.5 perturbations
-perturkey=int(raw_input("Perturbation? (1:none, 2: large scale, 3: small scale, 4: small+large, 5: n sigmas) "))
-if (perturkey > 1) :
-  seedin=int(raw_input("seedin? (only matters if adding perturbations) "))
-if ((perturkey == 3) or (perturkey == 4)) :
-  gwlength=float(raw_input("Gravity wave length? (for small scale perturbations) "))
-
-# 1.6 extra outputs
-# here we only implement an all-or-nothing case
-extvarkey=int(raw_input("Output the extra variables? (yes==1; no==0) "))
-if (extvarkey == 0) :
-  extvarkeys = np.zeros(100)
-else :
-  extvarkeys = np.ones(100)
-
-# 2. Call MCD
-(pres, dens, temp, zonwind, merwind, \
- meanvar, extvar, seedout, ierr) \
- = \
-call_mcd(zkey,xz,lon,lat,hrkey, \
- datekey,xdate,loct,dset,dust, \
- perturkey,seedin,gwlength,extvarkeys )
-
-# 3. Write outputs
-print "temperature is %.0f K, pressure is %.0f Pa, density is %5.3e kg/m3, zonal wind is %.1f m/s, meridional wind is %.1f m/s" % (temp,pres,dens,zonwind,merwind)
-
