Index: /trunk/UTIL/PYTHON/planetoplot.py
===================================================================
--- /trunk/UTIL/PYTHON/planetoplot.py	(revision 724)
+++ /trunk/UTIL/PYTHON/planetoplot.py	(revision 725)
@@ -11,4 +11,8 @@
 ### J. Leconte   -- LMD --    02/2012 -- Added area weighted averaging. Compatibility with terrestrial gcm.
 ### A. Spiga     -- LMD --    03/2012 -- Cleaning and improved comments 
+### T. Navarro   -- LMD --    04/2012 -- Added capabilities (e.g. histograms for difference maps)
+### A. Colaitis  -- LMD -- 05~06/2012 -- Added capabilities for analysis of mesoscale files (wind speed, etc...)
+### A. Spiga     -- LMD -- 04~07/2012 -- Added larger support of files + planets. Corrected a few bugs. Cleaning and improved comments
+
 def planetoplot (namefiles,\
            level=0,\
@@ -83,5 +87,7 @@
     from mymath import deg,max,min,mean,get_tsat,writeascii,fig2data,fig2img
     import matplotlib as mpl
-    from matplotlib.pyplot import contour,contourf,hist, text,subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel, title, close, legend, xlabel, axis, ylabel, subplots_adjust, axes, clabel
+    from matplotlib.pyplot import contour,contourf,hist, text,subplot, figure, rcParams, savefig, colorbar, \
+                                  pcolor, show, plot, clabel, title, close, legend, xlabel, axis, ylabel, \
+                                  subplots_adjust, axes, clabel
     from matplotlib.cm import get_cmap
     #from mpl_toolkits.basemap import cm
@@ -156,6 +162,12 @@
              print "WELL... nothing about time axis. I took default: first time reference stored in file."
     ### we get the names of variables to be read. in case only one available, we choose this one.
+    ### (we take out of the test specific names e.g. UV is not in the file but used to ask a wind speed computation)
       varname=var[vvv] 
-      if ((varname not in nc.variables) and not ((typefile in ['meso']) and (varname in ['UV','uv','uvmet','slopexy','SLOPEXY','deltat','DELTAT'])) and not ((typefile in ['gcm']) and (varname in ['enfact']))): 
+      logical_novarname = varname not in nc.variables
+      specificname_meso = ['UV','uv','uvmet','slopexy','SLOPEXY','deltat','DELTAT']
+      logical_nospecificname_meso = not ((typefile in ['meso']) and (varname in specificname_meso))
+      specificname_gcm = ['enfact']
+      logical_nospecificname_gcm = not ((typefile in ['gcm']) and (varname in specificname_gcm))
+      if ( logical_novarname and logical_nospecificname_meso and logical_nospecificname_gcm ): 
           if len(varinfile) == 1:   varname = varinfile[0] 
           else:                     varname = False
@@ -178,33 +190,35 @@
          elif zarea is not None:           [wlon,wlat] = latinterv(area=zarea)
 
-##########################################################
-############ LOAD 4D DIMENSIONS : x, y, z, t #############
-##########################################################
+#############################################################
+############ WE LOAD 4D DIMENSIONS : x, y, z, t #############
+#############################################################
+
+    ### TYPE 1 : GCM files or simple files
       if typefile in ["gcm","earthgcm","ecmwf"]:
+      ### this is needed for continuity 
           if slon is not None: sslon = slon  
           if slat is not None: sslat = slat  
-          ### LAT : done in getcoorddef
-          lat = lat2d[0,:] 
-          ### LON : done in getcoorddef
-          lon = lon2d[:,0]
-          ### ALT
+      ### we define lat/lon vectors. we get what was done in getcoorddef.
+          lat = lat2d[0,:] ; lon = lon2d[:,0]
+      ### we define areas. this is needed for calculate means and weight with area. this is not compulsory (see reduce_field).
+          if "aire" in nc.variables:      area = nc.variables["aire"][:,:]
+          ### --> add a line here if your reference is not present
+          else:                           area = None
+      ### we define altitude vector. either it is referenced or it is guessed based on last variable's dimensions.
           if "altitude" in nc.variables:   vert = nc.variables["altitude"][:]
           elif "Alt" in nc.variables:      vert = nc.variables["Alt"][:]
           elif "lev" in nc.variables:      vert = nc.variables["lev"][:]
           elif "presnivs" in nc.variables: vert = nc.variables["presnivs"][:]
+          ### --> add a line here if your reference is not present
           else: 
-              print "No altitude found. Try to build a simple axis." 
-              dadim = getdimfromvar(nc) 
+              print "No altitude found. Try to build a simple axis." ; dadim = getdimfromvar(nc) 
               if   len(dadim) == 4:  print "-- 4D field. Assume z is dim 2." ; vert = np.arange(dadim[-3])
               elif len(dadim) == 3:  print "-- 3D field. Assume z is dim 1." ; vert = [0.]
-              else:                  errormess("-- 2D or 1D field. Unclear. Stop.")
-          #else:                            vert = [0.]
-          ### AIRE (to weight means with the area)
-          if "aire" in nc.variables:      area = nc.variables["aire"][:,:]  
-          else:                           area = None
-          ### TIME
+              else:                  vert = [0.]
+      ### we define time vector. either it is referenced or it is guessed based on last variable's dimensions.
           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
+          elif "time_counter" in nc.variables:	time = nc.variables["time_counter"][:]/86400. #### convert from s to days
+          ### --> add a line here if your reference is not present
           else: 
               print "No time found. Try to build a simple axis. Assume t is dim 1." 
@@ -212,6 +226,6 @@
               if   len(dadim) == 4:  time = np.arange(dadim[-4])
               elif len(dadim) == 3:  time = np.arange(dadim[-3])
-              else:                  errormess("-- I am sorry this is 2D or 1D field. I failed.")
-          #else:                           	time = [0.] #errormess("no time axis found.")
+              else:                  time = [0.] #errormess("no time axis found.")
+          ### (SPECIFIC. convert to Ls for Martian GCM files.)
           if axtime in ["ls"]:
               print "converting to Ls ..."
@@ -219,18 +233,20 @@
                 time[iii] = sol2ls(time[iii])
                 if iii > 0:
-                  while abs(time[iii]-time[iii-1]) > 300:
-                    time[iii] = time[iii]+360
-          # for 1D plots (no need for longitude computation):
+                  while abs(time[iii]-time[iii-1]) > 300: time[iii] = time[iii]+360
+          ### (a case where the user would like to set local times. e.g. : 1D plot with no longitude reference)
           if axtime in ["lt"]:
               if initime == -1: initime=input("Please type initial local time:")
-              time = (initime+time*24)%24
-              print "LOCAL TIMES.... ", time
+              time = (initime+time*24)%24 ; print "LOCAL TIMES.... ", time
+
+    ### TYPE 2 : MESOSCALE FILES
       elif typefile in ['meso','geo']:
-          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
+          ### area not active with mesoscale files
+          area = None 
+          ### HACK TO GET THE CORRECT LAT/LON FROM MESOSCALE FILES WITH 2D LAT/LON ARRAYS
+          ### principle: calculate correct indices then repopulate slon and slat
           if slon is not None or slat is not None:
-              if firstfile and save == 'png' and typefile == 'meso' and "HGT" in varinfile and display:   iwantawhereplot = nc     #show a topo map with a cross on the chosen point
-              else:                                                                           iwantawhereplot = None   #do not show anything, just select indices
+              show_topo_map_user = ((save == 'png') and (typefile == 'meso') and ("HGT" in varinfile) and display)
+              if firstfile and show_topo_map_user:  iwantawhereplot = nc     #show a topo map with a cross on the chosen point
+              else:                                 iwantawhereplot = None   #do not show anything, just select indices
               numlon = 1 ; numlat = 1 
               if slon is not None:   numlon = slon.shape[1]    
@@ -249,11 +265,12 @@
                  if slat is not None: sslat[0][jjj] = indices[0,jjj,0] #...this is idy
               lonp,latp = ( lon2d[indices[0,0,0],indices[0,0,1]] , lat2d[indices[0,0,0],indices[0,0,1]] )
-          ######
-          if typefile in ['meso'] and mapmode == 1: lon2d = dumpbdy(lon2d,6) ; lat2d = dumpbdy(lat2d,6)  ### important to do that now and not before
-          ######
+          ### we get rid of boundary relaxation zone for plots. important to do that now and not before.
+          if typefile in ['meso'] and mapmode == 1: lon2d = dumpbdy(lon2d,6) ; lat2d = dumpbdy(lat2d,6)  
+          ### we read the keyword for vertical dimension. we take care for vertical staggering.
           if varname in ['PHTOT','W']:    vertdim='BOTTOM-TOP_PATCH_END_STAG'
           else:                           vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
           if (var2 is not None and var2 not in ['PHTOT','W']): dumped_vert_stag=True ; vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
           else:                                                dumped_vert_stag=False
+          ### we read the keyword for horizontal dimensions. we take care for horizontal staggering.
           if varname in ['V']:  latdim='SOUTH-NORTH_PATCH_END_STAG'
           else:                 latdim='SOUTH-NORTH_PATCH_END_UNSTAG'
@@ -261,5 +278,5 @@
           else:                 londim='WEST-EAST_PATCH_END_UNSTAG'
           lon = np.arange(0,getattr(nc,londim),1) ; lat = np.arange(0,getattr(nc,latdim),1)
-          ###
+          ### we define the time axis and take care of various specificities (lt, ls, sol) or issues (concatenation)
           if axtime in ["ls","sol"]:
               lstab, soltab, lttab = getlschar ( namefile, getaxis = True )
@@ -279,6 +296,5 @@
               time=check_localtime(time)
               print "LOCAL TIMES.... ", time
-
-          ###
+          ### we define the vertical axis (or lack thereof) and cover possibilities for it to be altitude, pressure, geopotential. quite SPECIFIC.
           if typefile in ['geo']:   vert = [0.] ; stime = readslices(str(0))
           else:
@@ -290,15 +306,9 @@
               elif vertmode == 1 or vertmode == 2:  vert = nc.variables["vert"][:]        ## pressure in Pa
               else:                                 vert = nc.variables["vert"][:]/1000.  ## altitude in km
-       #if firstfile:
-       #   lat0 = lat
-       #elif len(lat0) != len(lat):
-       #   errormess("Not the same latitude lengths !", [len(lat0), len(lat)])
-       #elif sum((lat == lat0) == False) != 0:
-       #   errormess("Not the same latitudes !", [lat,lat0])
-       ## Faire d'autre checks sur les compatibilites entre fichiers!!
-##########################################################
-##########################################################
-##########################################################
-
+####################################################################
+############ END of WE LOAD 4D DIMENSIONS : x, y, z, t #############
+####################################################################
+
+    ### we fill the arrays of varname, namefile, time at the current step considered (NB: why use both k and nnn ?)
       all_varname[k] = varname
       all_namefile[k] = namefile
@@ -306,12 +316,10 @@
       if var2: 
          all_var2[k] = getfield(nc,var2)
+         ### v--- too SPECIFIC (see below)
          if ((var2 in ['slopexy','SLOPEXY']) and (typefile in ['meso']) and (var2 not in nc.variables)): all_var2[k] = slopeamplitude(nc)
       if winds: all_windu[k] = getfield(nc,uchar) ; all_windv[k] = getfield(nc,vchar)
-      #####################
-      ##### GETFIELDS #####
-      #####################
-      ##### SPECIFIC CASES
-      ##### note : we could probably call those via a "toolbox" in myplot.
-      ##### 1. saturation temperature
+    ### we fill the arrays of fields to be plotted at the current step considered
+      ### ----------- SPECIFIC CASES. NOT HAPPY WITH THIS. note : we could probably call those via a "toolbox" in myplot.
+      ### ----------- 1. saturation temperature
       if varname in ["temp","t","T_nadir_nit","T_nadir_day","temp_day","temp_night"] and tsat:
           tt=getfield(nc,varname) ; print "computing Tsat-T, I ASSUME Z-AXIS IS PRESSURE"
@@ -319,20 +327,22 @@
           else:                                 tinput=tt
           all_var[k]=get_tsat(vert,tinput,zlon=lon,zlat=lat,zalt=vert,ztime=time)
-      ##### 2. wind amplitude
+      ### ----------- 2. wind amplitude
       elif ((varname in ['UV','uv','uvmet']) and (typefile in ['meso']) and (varname not in nc.variables)):
           all_var[k]=windamplitude(nc)
       elif ((varname in ['slopexy','SLOPEXY']) and (typefile in ['meso']) and (varname not in nc.variables)):
           all_var[k]=slopeamplitude(nc)
-      #### 3. Near surface instability
+      ### ------------ 3. Near surface instability
       elif ((varname in ['DELTAT','deltat']) and (typefile in ['meso']) and (varname not in nc.variables)):
           all_var[k]=deltat0t1(nc)
-      #### 4. Enrichment factor
+      ### ------------ 4. Enrichment factor
       elif ((typefile in ['gcm']) and (varname in ['enfact'])):
           all_var[k]=enrichment_factor(nc)
       else:
-      ##### GENERAL STUFF HERE
+      ### ideally only this line should be here
           all_var[k] = getfield(nc,varname)
+    ### we inform the user about the loop then increment the loop. this is the last line of "for namefile in namefiles" 
       print "**** GOT SUBDATA:",k," NAMEFILE:",namefile," VAR:",varname, var2 ; k += 1 ; firstfile = False
-      #### End of for namefile in namefiles
+
+    ### note to self : stopped comment rewriting here.
 
     ##################################
