use lapse rates for 1d profile and adapt qvaps func

examples/rainshaft1d/rainshaft1d.py

@@ -73,15 +73,12 @@
 
 ### --- settings for 1-D Thermodynamics --- ###
 PRESS0 = 101315 # [Pa]
-THETA = 298.15  # [K]
+THETA0 = 298.15  # [K]
 qcond = 0.0     # [Kg/Kg]
-WMAX = 0.6      # [m/s]
-VVEL = None     # [m/s]
-Zlength = 1500  # [m]
-Xlength = 1500  # [m]
-qvapmethod = "sratio"
+WVEL = 0.0      # [m/s]
+qvapmethod = "qvap"
 Zbase = 750 # [m]
-sratios = [0.99, 1.0025] # s_ratio [below, above] Zbase
+qvap0 = [0.99, 1.0025] # s_ratio [below, above] Zbase
 moistlayer=False
 
 ### --- settings for initial superdroplets --- ###
@@ -129,10 +126,9 @@
 rgrid.print_domain_info(constsfile, gridfile)
 
 ### ----- write thermodynamics binaries ----- ###
-thermodyngen = thermogen.ConstDryHydrostaticAdiabat(configfile, constsfile, PRESS0, 
-                                                 THETA, qvapmethod, sratios, Zbase,
-                                                 qcond, WMAX, Zlength, Xlength,
-                                                 VVEL, moistlayer)
+thermodyngen = thermogen.ConstHydrostaticLapseRates(configfile, constsfile, PRESS0, 
+                                                    THETA0, qvapmethod, sratios, Zbase,
+                                                    qcond, WVEL, VVEL, UVEL)
 cthermo.write_thermodynamics_binary(thermofile, thermodyngen, configfile,
                                     constsfile, gridfile)
 

pySD/thermobinary_src/thermogen.py

@@ -71,21 +71,26 @@ def sratio2qvap(sratio, press, temp, Mr_ratio):
 
   return qvap
 
-def qparams_to_qvap(method, qparams, Mr_ratio, PRESS, TEMP):
-  ''' determine qvap above and below z (cloud) base '''
+def qparams_to_qvap(method, params, Mr_ratio, PRESS, TEMP):
+  ''' returns qvaps given list of qvaps, supersaturation ratios
+  or relative humidities '''
 
   if method == "qvap":
+    qparams = params
     return qparams
 
   elif method == "sratio":
-    q0 = sratio2qvap(qparams[0], PRESS, TEMP, Mr_ratio)
-    q1 = sratio2qvap(qparams[1], PRESS, TEMP, Mr_ratio)
-    return [q0, q1]
+    qparams = []
+    for sratio in params:
+      qparams.append(sratio2qvap(sratio, PRESS, TEMP, Mr_ratio))
+    return qparams
 
   elif method == "relh":
-    q0 = relh2qvap(PRESS, TEMP, qparams[0], Mr_ratio)  
-    q1 = relh2qvap(PRESS, TEMP, qparams[1], Mr_ratio)     
-    return [q0, q1]
+    qparams = []
+    for relh in params:
+      qparams.append(relh2qvap(PRESS, TEMP, relh, Mr_ratio))
+    return qparams
+
   else:
     raise ValueError("valid method not given to generate qvap")
 
@@ -188,7 +193,7 @@ def __init__(self, configfile, constsfile, PRESS, TEMP, qvapmethod,
     self.TEMP = TEMP                          # temperature [T]
     self.qcond = qcond                        # liquid water content []
 
-    # determine qvap above & below z (cloud) base
+    # determine qvap [below, above] z (cloud) base
     self.Zbase = Zbase
     qvaps = qparams_to_qvap(qvapmethod, qvapparams,
                             inputs["Mr_ratio"], PRESS, TEMP)
@@ -288,7 +293,7 @@ def __init__(self, configfile, constsfile, PRESSz0, THETA,
     self.Xlength = Xlength # wavelength of velocity modulation in x direction [m]
     self.VVEL = VVEL # horizontal (y) velocity
 
-    # determine qvap above & below z (cloud) base
+    # determine qvap [below, above] z (cloud) base
     self.Zbase = Zbase
     self.qvapmethod, self.qvapparams = qvapmethod, qvapparams
     self.qvapz0 = qparams_to_qvap(qvapmethod, qvapparams,
@@ -422,6 +427,21 @@ class ConstHydrostaticLapseRates:
   and in hydrostatic equillibrium and following adiabats
   with constant lapse rates above/below zbase. '''
 
+  def __init__(self, PRESS0, THETA0, qvap0, Zbase,
+               thetalapses, qvlapses, qcond, WVEL, UVEL, VVEL):
+    self.PRESS0 = PRESS0                      # surface pressure [Pa]
+    self.THETA0 = THETA0                      # surface temperature [T]
+    self.qvap0 = qvap0                        # surface water vapour content [Kg/Kg]
+    self.thetalapses = thetalapses            # theta lapse rates [below, above] Zbase
+    self.qvlapses = qvlapses                  # qv lapse rates [below, above] Zbase
+    self.Zbase = Zbase                        # cloud base height [m]
+
+    self.qcond = qcond                        # liquid water content [Kg/Kg]
+    self.WVEL = WVEL                          # vertical (z) velocity [m/s]
+    self.UVEL = UVEL                          # horizontal x velocity [m/s]
+    self.VVEL = VVEL                          # horizontal y velocity [m/s]
+
+
   def generate_thermo(self, gbxbounds, ndims, ntime):
 
     zfulls, xfulls, yfulls = rgrid.fullcoords_forallgridboxes(gbxbounds,