Changeset 2151 in lmdz_wrf for trunk

Timestamp:

Sep 20, 2018, 4:47:43 PM (7 years ago)

Author:

lfita

Message:

Trying to fix moist-adiabts in `SkewT_Plog'

Location:

trunk/tools

Files:

• drawing.py (modified) (1 diff)
• drawing_tools.py (modified) (8 diffs)

trunk/tools/drawing.py

@@ -8999 +8999 @@
     return
 
-#filen='/home/lluis/estudios/ChemGBsAs/tests/199807/obs/snd/UWyoming_snd_87576.nc'
-#values='time|0:auto:auto:Sounding!at!Ezeiza!airport!on!3rd!July!1998:png:yes'
-#varn='ta,tda,pres'
-#draw_SkewT(filen, values, varn)
+filen='/home/lluis/estudios/ChemGBsAs/tests/199807/obs/snd/UWyoming_snd_87576.nc'
+values='time|0:auto:auto:Sounding!at!Ezeiza!airport!on!3rd!July!1998:png:yes'
+varn='ta,tda,pres'
+draw_SkewT(filen, values, varn)
+
 
 def draw_multi_SkewT(ncfiles, values):

trunk/tools/drawing_tools.py

@@ -11697 +11697 @@
     return
 
-def SaturationPressure(temp):
+def es(temp):
     """ Provides de saturation pressure [hPa] at a given temperature following The 
       August-Roche-Magnus equation
@@ -11704 +11704 @@
     23.334406231
     """
-    fname = 'SaturationPressure'
+    fname = 'es'
 
     ARM1 = 6.1094 # hPa
@@ -11711 +11711 @@
     
     tempC = temp - 273.15  
-    satpres = ARM1*exp(ARM2*tempC/(tempC+ARM3))
-
-    return satpres
+    esv = ARM1*exp(ARM2*tempC/(tempC+ARM3))
+
+    return esv
 
 #print 'satpres:', SaturationPressure(293.15)
 
-def SatMixingRatio(ta,pres):
+def qs(ta,pres):
     """ Provides de saturation mixing ratio [kgkg-1]
       ta: temperature [K]
@@ -11724 +11724 @@
     0.014531424608
     """
-    fname = 'SatMixingRatio'
-
-    es = SaturationPressure(ta)
-    rs = (0.6222*es)/(0.01*pres-es)
-
-    return rs
-
-#print 'saturation mixing ratio:', SatMixingRatio(293.,101300.)
+    fname = 'qs'
+
+    esv = es(ta)
+    qsv = (0.6222*es)/(0.01*pres-es)
+
+    return qsv
+
+#for t in range(260,310):
+#    print t, 'qs:', SatMixingRatio(t*1.,101300.)
+#quit()
+
+def e(ta, pres):
+    """ Function to compute water vapour pressure [hPa]
+      ta: air temperature [K]
+      pres: pressure [Pa]
+    """
+    fname = 'e'
+    tk = ta
+    data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
+    ev = 0.622*data1/(0.01*pres-(1.-0.622)*data1)
+
+    return ev
+
+def hur(ta, pres, hus):
+    """ Function to compute relative humidtiy [1]
+      ta: air temperature [K]
+      pres: pressure [Pa]
+      hus: water vapor mixing ratio [kgkg-1]
+    """
+    fname = 'hur'
+    
+    ev = e(ta,pres)
+    hurv = hus/ev
+
+    return hurv
+
+def td(ta,pres,hus):
+    """ Function to compute dew-point temperature [C]
+      ta: air temperature [K]
+      pres: pressure [Pa]
+      hus: water vapor mixing ratio [kgkg-1]
+    """
+    fname = 'td'
+    
+    hurv = hur(ta,pres,hus)
+    esv = es(ta)
+
+    pa = hurv*esv
+    tdv = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
+
+    return tdv
+
+def tl(ta, pres, hus):
+    """ Temperature at the Lifted condensation Level [C]
+      ta: air temperature [K]
+      pres: pressure [Pa]
+      hus: water vapor mixing ratio [kgkg-1]
+    """
+    fname = 'tl'
+
+    tdv = td(ta, pres, hus)
+
+    tC= ta- 273.15
+
+    frac = 1./(tdv-56.) + np.log(tC/tdv)/800.
+
+    tlv = 1./frac + 56.
+
+    return tlv
+
+def equivpot(ta, pres, hus):
+    """ Equivalent potenatial temperature as in WikiPedia
+    https://en.wikipedia.org/wiki/Equivalent_potential_temperature
+      ta: temperature of the air [K]
+      pres: pressure of the air [Pa]
+      hus: mixing ratio of the air [kgkg-1]
+    """
+    fname = 'equivpot'
+
+    ev = e(ta, pres)
+    tlv = tl(ta,pres,hus) + 273.15
+
+    # Dry potential temperature at LCL
+    thetal_1 = ta*(pres-ev)**0.2854
+    thetal_2 = (ta/tlv)**(0.28*hus)
+    thetal = thetal_1*thetal_2
+
+    ethetav_exp = (3036./tlv-1.78)*hus*(1.+0.448*hus)
+    ethetav = thetal*np.exp(ethetav_exp)
+
+    return ethetav
+
+print equivpot(300., 98000., 0.02)
+quit()
 
 def ethetapotct(tsfc, dz=50, pmin=5000.):
@@ -11758 +11844 @@
     dp = (psfc-pmin)/(dz-1)
     
-    ws = SatMixingRatio(tsfc+273.15,psfc)
+    wsv = ws(tsfc+273.15,psfc)
     thetasfc = (tsfc+273.15)*(psfc/p0)**(R/Cp)
     for iz in range(dz):
         ethetapotv[iz,0] = psfc-dp*iz
-        thetapotv = (tsfc+273.15)*((psfc-dp*iz)/p0)**(R/Cp)
-        tfromthetasfc = thetasfc*(p0/(psfc-dp*iz))**(-R/Cp)
-        wstfromthetasfc = SatMixingRatio(tfromthetasfc,psfc-dp*iz)
-        #ethetapotv[iz,1] = thetapotv*np.exp(L*1000.*ws/(Cp*(tsfc+273.15))) -273.15
-        ethetapotsfcv = thetapotv*np.exp(L*1000.*wstfromthetasfc/(Cp*(tsfc+273.15))) -273.15
-        ethetapotv[iz,1] = ethetapotsfcv
-        print iz,psfc-dp*iz,thetapotv-273.15,tfromthetasfc-273.15,':',ethetapotv[iz,1],wstfromthetasfc,ethetapotsfcv
+        thetav = (tsfc+273.15)*((psfc-dp*iz)/p0)**(R/Cp)
+        temp = thetav*(p0/(psfc-dp*iz))**(-R/Cp)
+        wsv = ws(temp,psfc-dp*iz)
+        ethetapotval = thetav*np.exp(L*1000.*wsv/(Cp*(tsfc+273.15)))-273.15
+        ethetapotv[iz,1] = ethetapotval
+        print iz,psfc-dp*iz,thetav-273.15,temp-273.15,':',ethetapotv[iz,1],ws
 
     print thetasfc,ws
-    #quit()
+    quit()
     return ethetapotv
 
@@ -11809 +11894 @@
     """
     import skewt
+
+    lbck = 1.
 
     fname = 'plot_SkewT'
@@ -11832 +11919 @@
     ax = fig.add_subplot(111, projection='skewx')
 
-    plt.grid(True)
+    plt.grid(True, linewidth=lbck)
 
     # Quantity of dryadiabats
-    xmax = taxtrm[1]
+    xmax = 150.
     xmin = taxtrm[0]
     if len(np.arange(xmin,xmax,10.)) < 10: xfrqtk = 10.
@@ -11844 +11931 @@
     for itt in range(ddt):
         dryt = thetapotct(xmin+itt*xfrqtk, pmin=10000.)
-        ax.semilogy(dryt[:,1], dryt[:,0]/100., color='#FFCCCC', linewidth=0.75)
+        ax.semilogy(dryt[:,1], dryt[:,0]/100., '-.', color='#FFCCCC', linewidth=lbck)
 
     # Including moist- adiabats
-    #for itt in range(ddt):
-    #    moist = ethetapotct(xmin+itt*xfrqtk, pmin=10000.)
-    #    if xmin+itt*xfrqtk == 10.: print moist[:,1], moist[:,0]/100.
-    #    ax.semilogy(moist[:,1], moist[:,0]/100., color='#CCFFCC', linewidth=0.75)   
+    for itt in range(ddt):
+        moist = ethetapotct(xmin+itt*xfrqtk, pmin=10000.)
+        if xmin+itt*xfrqtk == 10.: print moist[:,1], moist[:,0]/100.
+        ax.semilogy(moist[:,1], moist[:,0]/100., '--', color='#CCFFCC', linewidth=lbck)
 
     # Plot the data using normal plotting functions, in this case using