Changeset 1718 in lmdz_wrf

Timestamp:

Dec 14, 2017, 1:32:14 AM (7 years ago)

Author:

lfita

Message:

Adding:

• `var_hur_Uhus': Function to compute relative humidity following 'August - Roche - Magnus' formula using hus

File:

• trunk/tools/diag_tools.py (modified) (4 diffs)

trunk/tools/diag_tools.py

@@ -37 +37 @@
 # var_clt: Function to compute the total cloud fraction following 'newmicro.F90' from LMDZ using 1D vertical column values
 # var_hur: Function to compute relative humidity following 'August - Roche - Magnus' formula
+# var_hur_Uhus: Function to compute relative humidity following 'August - Roche - Magnus' formula using hus
 # var_mslp: Fcuntion to compute mean sea-level pressure
 # var_td: Function to compute dew-point air temperature from temperature and pressure values
@@ -1295 +1296 @@
     hur = q/(ws*1000.)
 
-    print 'q:', q[5,5,5,5], 't:', t[5,5,5,5], 'p:', p[5,5,5,5]
-    print 'es:', es[5,5,5,5], 'ws:', ws[5,5,5,5], 'hur:', hur[5,5,5,5]
+    #print 'q:', q[5,5,5,5], 't:', t[5,5,5,5], 'p:', p[5,5,5,5]
+    #print 'es:', es[5,5,5,5], 'ws:', ws[5,5,5,5], 'hur:', hur[5,5,5,5]
+
+    return hur
+
+def var_hur_Uhus(p, t, hus):
+    """ Function to compute relative humidity following 'August - Roche - Magnus' formula using hus
+      [t]= temperature (assuming [[t],z,y,x] in [K])
+      [p] = pressure field (assuming in [Pa])
+      [hus] = specific humidty [1]
+      ref.: M. G. Lawrence, BAMS, 2005, 225
+    >>> print var_rh(101300., 290., 3.)
+    0.250250256174
+    """
+    fname = 'var_hur_Uhus'
+
+    # Enthalpy of vaporization [Jkg-1]
+    L = 2.453*10.**6
+    # Gas constant for water vapor [JK-1Kg-1]
+    Rv = 461.5 
+
+    # Computing saturated pressure 
+    #es = 10.*0.6112*np.exp(17.67*(t-273.16)/(t-29.65))
+    #es = 6.11*np.exp(L/Rv*(1.-273./t)/273.)
+
+    # August - Roche - Magnus formula
+    es = 6.1094*np.exp(17.625*(t-273.15)/((t-273.15)+243.04))
+
+    # Saturated mixing ratio [g/kg]
+    ws = 0.622*es/(0.01*p-es)
+
+    # Mixing ratio
+    q = hus/(1.+hus)
+
+    hur = q/(ws*1000.)
+
+    #print 'q:', q[5,5,5,5], 't:', t[5,5,5,5], 'p:', p[5,5,5,5]
+    #print 'es:', es[5,5,5,5], 'ws:', ws[5,5,5,5], 'hur:', hur[5,5,5,5]
 
     return hur
@@ -1367 +1404 @@
         if Cdiag == 'hur':
             """ Computing relative humidity
+            """
+            vn = 'hur'
+            CF3Dvars = ['ta', 'plev', 'q']
+            for v3D in CF3Dvars:
+                if not vars3D.has_key(v3D):
+                    print gen.errormsg
+                    print '  ' + fname + ": missing variable '" + v3D +              \
+                      "' attribution to compute '" + vn + "' !!"
+                    print '  Equivalence of 3D variables provided _______'
+                    gen.printing_dictionary(vars3D)
+                    quit(-1)
+                if not self.incvars.has_key(vars3D[v3D]):
+                    print gen.errormsg
+                    print '  ' + fname + ": missing variable '" + vars3D[v3D] +      \
+                      "' in input file to compute '" + vn + "' !!"
+                    print '  available variables:', self.incvars.keys()
+                    print '  looking for variables _______'  
+                    gen.printing_dictionary(vars3D)
+                    quit(-1)
+
+            ta = ncobj.variables[vars3D['ta']][:]
+            p = ncobj.variables[vars3D['plev']][:]
+            q = ncobj.variables[vars3D['q']][:]
+
+            if len(ta.shape) != len(p.shape):
+                p = gen.fill_Narray(p, ta*0., filldim=[0,2,3])
+
+            self.values = var_hur(p, ta, q)
+            self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
+              dictdims['lon']]
+            self.units = '1'
+
+        if Cdiag == 'hur_Uhus':
+            """ Computing relative humidity using hus
             """
             vn = 'hur'
@@ -1394 +1465 @@
                 p = gen.fill_Narray(p, ta*0., filldim=[0,2,3])
 
-            self.values = var_hur(p, ta, hus)
+            self.values = var_hur_Uhus(p, ta, hus)
             self.dims = [dictdims['time'], dictdims['plev'], dictdims['lat'],        \
               dictdims['lon']]