Changeset 753

Timestamp:

Aug 2, 2012, 11:52:55 AM (12 years ago)

Author:

acolaitis

Message:

Python. Minor updates to enrichment factor computation, made compatible with Yuan Lian et al 2012. (this is old, just commiting uncommited stuff)

Location:

trunk/UTIL/PYTHON

Files:

• myplot.py (modified) (4 diffs)
• planetoplot.py (modified) (1 diff)

trunk/UTIL/PYTHON/myplot.py

@@ -130 +130 @@
 # Compute the enrichment factor of non condensible gases
 # The corresponding variable to call is enfact
-def enrichment_factor(nc):
+# enrichment factor is computed as in Yuan Lian et al. 2012
+# i.e. you need to have VL2 site at LS 135 in your data
+# this only requires co2col so that you can concat.nc at low cost
+def enrichment_factor(nc,lon,lat,time):
     import numpy as np
+    from myplot import reducefield
     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 "co2col" in varinfile: co2col=getfield(nc,'co2col')
+    else: print "error, you need co2col 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
+    dimension = len(nc.variables['co2col'].dimensions)
+    if dimension == 2: 
+      zny,znx = np.array(co2col).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])
+    elif dimension == 3: znt,zny,znx = np.array(co2col).shape
     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)
-    for tt in np.arange(znt):
-      enfact[tt,:,:] =-(meanar[tt] - mmrarcol[tt,:,:])/meanar[tt]
-
+    mmrarcol[:,:,:] = 1. - grav*co2col[:,:,:]/ps[:,:,:]
+# Computation with reference argon mmr at VL2 Ls 135 (as in Yuan Lian et al 2012)
+    lonvl2=np.zeros([1,2])
+    latvl2=np.zeros([1,2])
+    timevl2=np.zeros([1,2])
+    lonvl2[0,0]=-180
+    lonvl2[0,1]=180
+    latvl2[:,:]=48.16
+    timevl2[:,:]=135.
+    indexlon  = getsindex(lonvl2,0,lon)
+    indexlat  = getsindex(latvl2,0,lat)
+    indextime = getsindex(timevl2,0,time)
+    mmrvl2, error = reducefield( mmrarcol, d4=indextime, d1=indexlon, d2=indexlat)
+    print "VL2 Ls 135 mmr arcol:", mmrvl2
+    enfact[:,:,:] = mmrarcol[:,:,:]/mmrvl2
     return enfact
 
@@ -236 +231 @@
     if dimension == 2:
         #### this is needed for 1d-type files (where dim=2 but axes are time-vert and not lat-lon)
-        if unidim == 1: d2=d4 ; d1=d3 ; d4=None ; d3=None
+        if unidim==1: d2=d4 ; d1=d3 ; d4 = None ; d3 = None
         if mesharea is None: mesharea=np.ones(shape)
         if   max(d2) >= shape[0]: error = True
@@ -631 +626 @@
     if typefile in ['meso']:
         if '9999' not in getattr(nc,'START_DATE') :   
-            ## regular mesoscale  
+            ## regular mesoscale 
             [lon2d,lat2d] = getcoord2d(nc)  
         else:                     
@@ -1056 +1051 @@
              "TAU":          "YlOrBr_r",\
              "CO2":          "YlOrBr_r",\
+             "MIXED":        "GnBu",\
              #### T.N.
              "MTOT":         "spectral",\

trunk/UTIL/PYTHON/planetoplot.py

@@ -337 +337 @@
       ### ------------ 4. Enrichment factor
       elif ((typefile in ['gcm']) and (varname in ['enfact'])):
-          all_var[k]=enrichment_factor(nc)
+          all_var[k]=enrichment_factor(nc,lon,lat,time)
       else:
       ### ideally only this line should be here