Changeset 701 for trunk/UTIL

Timestamp:

Jun 11, 2012, 4:17:57 PM (13 years ago)

Author:

acolaitis

Message:

PYTHON. Added possibility to call for variable <<-v enfact>> - enrichment factor - given that you have co2, ap, bp and ps fields in your file. This will give you the enrichment factor of non condensible gases, ranging from 0 (no enrichment) to X (enrichment by a factor X compared to the mean non condensible gases mass mixing ratio (which is also computed for the relevant time step)).

Location:

trunk/UTIL/PYTHON

Files:

• myplot.py (modified) (2 diffs)
• planetoplot.py (modified) (4 diffs)

trunk/UTIL/PYTHON/myplot.py

@@ -125 +125 @@
        zvint=zv
     return np.sqrt(zuint**2 + zvint**2)
+
+## Author: AC
+# Compute the enrichment factor of non condensible gases
+# The corresponding variable to call is enfact
+def enrichment_factor(nc):
+    import numpy as np
+    varinfile = nc.variables.keys()
+    if "co2" in varinfile: co2=getfield(nc,'co2')
+    else: print "error, you need co2 var in your file"
+    if "ap" in varinfile: ap=getfield(nc,'ap')
+    else: print "error, you need ap var in your file"
+    if "bp" in varinfile: bp=getfield(nc,'bp')
+    else: print "error, you need bp var in your file"
+    if "ps" in varinfile: ps=getfield(nc,'ps')
+    else: print "error, you need ps var in your file"
+    dimension = len(nc.variables['co2'].dimensions)
+    if dimension == 3: 
+      znz,zny,znx = np.array(co2).shape
+      znt=1
+    elif dimension == 4: znt,znz,zny,znx = np.array(co2).shape
+    co2col = np.zeros([znt,zny,znx])
+    arcol = np.zeros([znt,zny,znx])
+    mmrarcol = np.zeros([znt,zny,znx])
+    meanar = np.zeros([znt])
+    pplev = np.zeros([znt,znz+1,zny,znx])
+    enfact = np.zeros([znt,zny,znx])
+    grav=3.72
+    for zz in np.arange(znz):
+      pplev[:,zz,:,:] = ap[zz]+bp[zz]*ps[:,:,:]
+    pplev[:,znz,:,:]=0.
+
+    for zz in np.arange(znz):
+      co2col[:,:,:] = co2col[:,:,:] + co2[:,zz,:,:]*(pplev[:,zz,:,:]-pplev[:,zz+1,:,:])/grav
+      arcol[:,:,:] = arcol[:,:,:] + (1.-co2[:,zz,:,:])*(pplev[:,zz,:,:]-pplev[:,zz+1,:,:])/grav
+
+    mmrarcol = arcol/(arcol + co2col)
+    for xx in np.arange(znx):
+      for yy in np.arange(zny):
+          meanar[:] = meanar[:] + mmrarcol[:,yy,xx]
+    meanar = meanar/(znx*zny)
+
+    enfact =-(meanar - mmrarcol)/meanar
+    return enfact
 
 ## Author: AC
@@ -947 +990 @@
              "TAU":          "%.1f",\
              "CO2":          "%.2f",\
+             "ENFACT":       "%.1f",\
              #### T.N.
              "TEMP":         "%.0f",\

trunk/UTIL/PYTHON/planetoplot.py

@@ -79 +79 @@
                        zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow,\
                        getname,localtime,check_localtime,polarinterv,getsindex,define_axis,determineplot,readslices,bidimfind,getlschar,hole_bounds,\
-                       windamplitude,slopeamplitude,deltat0t1
+                       windamplitude,slopeamplitude,deltat0t1,enrichment_factor
     from mymath import deg,max,min,mean,get_tsat,writeascii,fig2data,fig2img
     import matplotlib as mpl
@@ -155 +155 @@
     ### we get the names of variables to be read. in case only one available, we choose this one.
       varname=var[vvv] 
-      if ((varname not in nc.variables) and not ((typefile in ['meso']) and (varname in ['UV','uv','uvmet','slopexy','SLOPEXY','deltat','DELTAT']))): 
+      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']))): 
           if len(varinfile) == 1:   varname = varinfile[0] 
           else:                     varname = False
@@ -299 +299 @@
       #####################
       ##### SPECIFIC CASES
+      ##### 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:
@@ -310 +311 @@
       elif ((varname in ['slopexy','SLOPEXY']) and (typefile in ['meso']) and (varname not in nc.variables)):
           all_var[k]=slopeamplitude(nc)
+      #### 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
+      elif ((typefile in ['gcm']) and (varname in ['enfact'])):
+          all_var[k]=enrichment_factor(nc)
       else:
       ##### GENERAL STUFF HERE