Changeset 1343 for trunk/UTIL

Timestamp:

Aug 26, 2014, 11:53:30 PM (10 years ago)

Author:

aslmd

Message:

python. mcd. bug fixes to improve ticks intervals, and max/min, and add warning if no latitude (tbd for others)

Location:

trunk/UTIL/PYTHON/mcd

Files:

• frozen_myplot.py (modified) (1 diff)
• inimeso/inimeso5.py (modified) (3 diffs)
• mcd.py (modified) (5 diffs)
• proto/cgi-bin/mcdcgi.py (modified) (1 diff)

trunk/UTIL/PYTHON/mcd/frozen_myplot.py

@@ -987 +987 @@
                       zevmax = mean(fieldcalc) + dev  ### for continuity            
     ###
-    if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0.
+    if zevmin < min(fieldcalc): zevmin = min(fieldcalc)
+    if zevmax > max(fieldcalc): zevmax = max(fieldcalc)
+    #if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0.
     #print "BOUNDS field ", min(fieldcalc), max(fieldcalc), " //// adopted", zevmin, zevmax
     return zevmin, zevmax

trunk/UTIL/PYTHON/mcd/inimeso/inimeso5.py

@@ -7 +7 @@
 from mcd import mcd
 
+rho_dust = 2500.  # Mars dust density (kg.m-3)
+grav = 3.72
+ksi = 3. / 4. / rho_dust / grav
+nueff = 0.5
+
 ### MCD INSTANCE and SETTINGS (actually, default. but one never knows)
 query = mcd() ; query.zkey = 3 ; query.dust = 2 ; query.hrkey = 1
-
 query.toversion5()
 
@@ -20 +24 @@
 ### OPEN FILES TO BE WRITTEN
 sounding = open("input_sounding", "w") ; additional = open("input_therm", "w") ; more = open("input_more", "w")
+dust = open("input_dust", "w")
 
 ### GET and WRITE SURFACE VALUES
 query.xz = 0.1 ; 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.) )
+query.extvar[42] = 1.5*1e-3
+query.extvar[44] = 0.0
+sounding.write( "%10.2f%12.2f%12.2f\n" % (query.pres/100.,query.temp*(610./query.pres)**(1.0/3.9),(query.extvar[42]+query.extvar[44])*1e3) )
 more.write( "%10.2f%10.2f" % (query.extvar[1],query.extvar[14]) ) ; more.close()
 
@@ -29 +36 @@
 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)):
+
+    query.extvartab[iz,42] = 1.5*1e-3
+    query.extvartab[iz,44] = 0.0
+
     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]) )
+                      (query.extvartab[iz,42]+query.extvartab[iz,44])*1e3,\
+                      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]) )
 
+    ### DUST PROFILES
+    q = query.extvartab[iz,38] # extvar(38)= Dust mass mixing ratio (kg/kg)
+    reff = query.extvartab[iz,39] # extvar(39)= Dust effective radius (m)
+    print q,reff
+    N = (grav*ksi*((1+nueff)**3)/np.pi)*q/(reff**3)
+    dust.write( "%18.6e%18.6e\n" % (q,N) )
+
 ### FINISH
-sounding.close() ; additional.close()
-query.plot1d(["p","t","u","v"]) ; mpl.show()
+sounding.close() ; additional.close() ; dust.close()
+query.plot1d(["p","t","u","v","h2ovap","h2oice"]) ; mpl.show()
+

trunk/UTIL/PYTHON/mcd/mcd.py

@@ -512 +512 @@
 
     def vertlabel(self):
+#      if self.zkey == 1:   self.xlabel = "radius from centre of planet (m)"
+#      elif self.zkey == 2: self.xlabel = "height above areoid (m) (MOLA zero datum)"
+#      elif self.zkey == 3: self.xlabel = "height above surface (m)"
+#      elif self.zkey == 4: self.xlabel = "pressure level (Pa)"
+#      elif self.zkey == 5: self.xlabel = "altitude above mean Mars Radius(=3396000m) (m)"
       if self.zkey == 1:   self.xlabel = "radius from centre of planet (m)"
-      elif self.zkey == 2: self.xlabel = "height above areoid (m) (MOLA zero datum)"
+      elif self.zkey == 2: self.xlabel = "altitude above MOLA$_0$ (m)"
       elif self.zkey == 3: self.xlabel = "height above surface (m)"
-      elif self.zkey == 4: self.xlabel = "pressure level (Pa)"
-      elif self.zkey == 5: self.xlabel = "altitude above mean Mars Radius(=3396000m) (m)"
+      elif self.zkey == 4: self.xlabel = "pressure (Pa)"
+      elif self.zkey == 5: self.xlabel = "altitude above mean Mars radius (m)"
 
     def vertunits(self):
@@ -665 +670 @@
         if not self.vertplot and self.islog: ax.set_yscale('log')
         if self.vertplot and self.islog: ax.set_xscale('log')
+
+        #ax.ticklabel_format(useOffset=False,axis='x')
+        #ax.ticklabel_format(useOffset=False,axis='y')
 
         if self.lats is not None:      ax.set_xticks(np.arange(-90,91,15)) ; ax.set_xbound(lower=self.lats, upper=self.late)
@@ -815 +823 @@
       for i in range(len(tabtodo)): mpl.subplot(subv,subh,i+1) ; self.makemap2d(tabtodo[i],incwind=incwind,proj=proj)
 
+    def makeinterv(self):
+      self.latinterv = 30.
+      self.loninterv = 45.
+      if self.lats is not None:
+        if (abs(self.late-self.lats) < 90.): self.latinterv = 10.
+        if (abs(self.late-self.lats) < 10.): self.latinterv = 1.
+      if self.lons is not None:
+        if (abs(self.lone-self.lons) < 135.): self.loninterv = 15.
+        if (abs(self.lone-self.lons) < 15.): self.loninterv = 1.
+
     def htmlmap2d(self,tabtodo,incwind=False,figname="temp.png",back="zMOL"):
     ### complete 2D figure with possible multiplots
@@ -890 +908 @@
           [x2d,y2d] = np.meshgrid(x,y)
           yeah.quiver(x2d,y2d,np.transpose(windx),np.transpose(windy))
+
         ax = fig.gca() ; ax.set_ylabel("Latitude") ; ax.set_xlabel("Longitude")
-        ax.set_xticks(np.arange(-360,361,45)) ; ax.set_xbound(lower=self.lons, upper=self.lone)
-        ax.set_yticks(np.arange(-90,91,30)) ; ax.set_ybound(lower=self.lats, upper=self.late)
+
+        # make intervals 
+        self.makeinterv()
+        ax.set_xticks(np.arange(-360,361,self.loninterv)) ; ax.set_xbound(lower=self.lons, upper=self.lone)
+        ax.set_yticks(np.arange(-90,91,self.latinterv)) ; ax.set_ybound(lower=self.lats, upper=self.late)
       self.gettitle()
       fig.text(0.5, 0.95, self.title, ha='center')
@@ -952 +974 @@
         ax = fig.gca() ; ax.set_ylabel(self.ylabel) ; ax.set_xlabel(self.xlabel)
 
-        if self.lons is not None:   ax.set_xticks(np.arange(-360,361,45)) ; ax.set_xbound(lower=self.lons, upper=self.lone)
-        elif self.lats is not None: ax.set_xticks(np.arange(-90,91,30)) ; ax.set_xbound(lower=self.lats, upper=self.late)
+        self.makeinterv()
+        if self.lons is not None:   ax.set_xticks(np.arange(-360,361,self.loninterv)) ; ax.set_xbound(lower=self.lons, upper=self.lone)
+        elif self.lats is not None: ax.set_xticks(np.arange(-90,91,self.latinterv)) ; ax.set_xbound(lower=self.lats, upper=self.late)
 
         if self.locts is not None: 

trunk/UTIL/PYTHON/mcd/proto/cgi-bin/mcdcgi.py

@@ -98 +98 @@
 
 # Get data from user-defined fields and define free dimensions
-islatfree,  query.lat,  query.lats,  query.late  = gethtmlcoord( form.getvalue("latitude"),   -90.,  90. )
+getlat = form.getvalue("latitude")
+if getlat is None: 
+  errormess = errormess+"<li>A value or interval for latitude is missing. Please correct."
+  islatfree = False ; query.lat = 0. ; query.lats = 0. ; query.late = 0.
+else:
+  islatfree,  query.lat,  query.lats,  query.late  = gethtmlcoord( getlat,   -90.,  90. )
 islonfree,  query.lon,  query.lons,  query.lone  = gethtmlcoord( form.getvalue("longitude"), -180., 180. )
 isloctfree, query.loct, query.locts, query.locte = gethtmlcoord( form.getvalue("localtime"),    0.,  24. )