Jon GCCN comparison

parcel.py

@@ -5,6 +5,9 @@
 #sys.path.insert(0, "../libcloudphxx/build/bindings/python/")
 #sys.path.insert(0, "../../../libcloudphxx/build/bindings/python/")
 #sys.path.insert(0, "/usr/local/lib/site-python/")
+#sys.path.insert(0, "/mnt/local/pdziekan/usr/local/lib")
+#sys.path.insert(0, "/mnt/local/pdziekan/lib/python/site-packages/")
+sys.path.insert(0, "/usr/local/lib/python2.7/dist-packages/")
 # TEMP TODO TEMP TODO !!!
 
 from argparse import ArgumentParser, RawTextHelpFormatter
@@ -67,13 +70,17 @@ def __call__(self, lnr):
       res += lognormal(lnr)
     return res
 
-def _micro_init(aerosol, opts, state, info):
+def _micro_init(aerosol, aerosol_sizes, opts, state, info):
 
   # lagrangian scheme options
   opts_init = lgrngn.opts_init_t()
   for opt in ["dt", "sd_conc", "chem_rho", "sstp_cond"]:
     setattr(opts_init, opt, opts[opt])
-  opts_init.n_sd_max = opts_init.sd_conc
+  #opts_init.n_sd_max = int(1.1*opts_init.sd_conc)
+  opts_init.n_sd_max = opts_init.sd_conc + 38  # some more space for the dry_sizes GCCNs
+
+  opts_init.kernel = lgrngn.kernel_t.hall_davis_no_waals
+  opts_init.terminal_velocity = lgrngn.vt_t.khvorostyanov_nonspherical;
 
   # read in the initial aerosol size distribution
   dry_distros = {}
@@ -84,14 +91,29 @@ def _micro_init(aerosol, opts, state, info):
     dry_distros[dct["kappa"]] = sum_of_lognormals(lognormals)
   opts_init.dry_distros = dry_distros
 
+  # read in the initial aerosol dry sizes
+  if (aerosol_sizes != None):
+    print aerosol_sizes
+    dry_sizes = {}
+    for name, dct in aerosol_sizes.iteritems(): # loop over kappas
+      print dct["kappa"]
+      print dct["r"]
+      print dct["n_tot"]
+      print dct["multi"]
+      sizes = {}
+      for i in range(len(dct["r"])):
+        sizes[dct["r"][i]] = [dct["n_tot"][i], dct["multi"][i]]
+      dry_sizes[dct["kappa"]] = sizes
+    opts_init.dry_sizes = dry_sizes
+
   # better resolution for the SD tail
   if opts["large_tail"]:
       opts_init.sd_conc_large_tail = 1
       opts_init.n_sd_max = int(1e6)  # some more space for the tail SDs
 
   # switch off sedimentation and collisions
   opts_init.sedi_switch = False
-  opts_init.coal_switch = False
+  opts_init.coal_switch = True
 
   # switching on chemistry if either dissolving, dissociation or reactions are chosen
   opts_init.chem_switch = False
@@ -100,22 +122,22 @@ def _micro_init(aerosol, opts, state, info):
     opts_init.sstp_chem = opts["sstp_chem"]
 
   # initialisation
-  micro = lgrngn.factory(lgrngn.backend_t.serial, opts_init)
+  micro = lgrngn.factory(lgrngn.backend_t.OpenMP, opts_init)
   ambient_chem = {}
   if micro.opts_init.chem_switch:
     ambient_chem = dict((v, state[k]) for k,v in _Chem_g_id.iteritems())
   micro.init(state["th_d"], state["r_v"], state["rhod"], ambient_chem=ambient_chem)
 
   # sanity check
-  _stats(state, info)
+  _stats(state, info, micro)
   if (state["RH"] > 1): raise Exception("Please supply initial T,p,r_v below supersaturation")
 
   return micro
 
 def _micro_step(micro, state, info, opts, it, fout):
   libopts = lgrngn.opts_t()
   libopts.cond = True
-  libopts.coal = False
+  libopts.coal = True
   libopts.adve = False
   libopts.sedi = False
 
@@ -136,7 +158,7 @@ def _micro_step(micro, state, info, opts, it, fout):
   micro.step_async(libopts)
 
   # update state after microphysics (needed for below update for chemistry)
-  _stats(state, info)
+  _stats(state, info, micro)
 
   # update in state for aqueous chem (TODO do we still want to have aq chem in state?)
   if micro.opts_init.chem_switch:
@@ -146,25 +168,27 @@ def _micro_step(micro, state, info, opts, it, fout):
       micro.diag_chem(id_int)
       state[id_str.replace('_g', '_a')] = np.frombuffer(micro.outbuf())[0]
 
-def _stats(state, info):
+def _stats(state, info, micro):
   state["T"] = np.array([common.T(state["th_d"][0], state["rhod"][0])])
   state["RH"] = state["p"] * state["r_v"] / (state["r_v"] + common.eps) / common.p_vs(state["T"][0])
+  micro.diag_RH()
+  state["RH_libcloud"] =  np.frombuffer(micro.outbuf())[0]
   info["RH_max"] = max(info["RH_max"], state["RH"])
 
 def _output_bins(fout, t, micro, opts, spectra):
   for dim, dct in spectra.iteritems():
     for bin in range(dct["nbin"]):
-      if dct["drwt"] == 'wet':
+      if dct["slct"] == 'wet':
 	micro.diag_wet_rng(
 	  fout.variables[dim+"_r_wet"][bin],
 	  fout.variables[dim+"_r_wet"][bin] + fout.variables[dim+"_dr_wet"][bin]
 	)
-      elif dct["drwt"] == 'dry':
+      elif dct["slct"] == 'dry':
 	micro.diag_dry_rng(
 	  fout.variables[dim+"_r_dry"][bin],
 	  fout.variables[dim+"_r_dry"][bin] + fout.variables[dim+"_dr_dry"][bin]
 	)
-      else: raise Exception("drwt should be wet or dry")
+      else: raise Exception("slct should be wet or dry")
 
       for vm in dct["moms"]:
         if type(vm) == int:
@@ -174,7 +198,7 @@ def _output_bins(fout, t, micro, opts, spectra):
           elif dct["drwt"] == 'dry':
             micro.diag_dry_mom(vm)
           else: raise Exception("drwt should be wet or dry")
-          fout.variables[dim+'_m'+str(vm)][int(t), int(bin)] = np.frombuffer(micro.outbuf())
+          fout.variables[dim+'_'+dct["drwt"]+'_mom'+str(vm)][int(t), int(bin)] = np.frombuffer(micro.outbuf())
         else:
           # calculate chemistry
           micro.diag_chem(_Chem_a_id[vm])
@@ -187,12 +211,12 @@ def _output_init(micro, opts, spectra):
   for name, dct in spectra.iteritems():
     fout.createDimension(name, dct["nbin"])
 
-    tmp = name + '_r_' + dct["drwt"]
+    tmp = name + '_r_' + dct["slct"]
     fout.createVariable(tmp, 'd', (name,))
     fout.variables[tmp].unit = "m"
-    fout.variables[tmp].description = "particle wet radius (left bin edge)"
+    fout.variables[tmp].description = "particle "+dct["slct"]+" radius (left bin edge)"
 
-    tmp = name + '_dr_' + dct["drwt"]
+    tmp = name + '_dr_' + dct["slct"]
     fout.createVariable(tmp, 'd', (name,))
     fout.variables[tmp].unit = "m"
     fout.variables[tmp].description = "bin width"
@@ -201,12 +225,12 @@ def _output_init(micro, opts, spectra):
       from math import exp, log
       dlnr = (log(dct["rght"]) - log(dct["left"])) / dct["nbin"]
       allbins = np.exp(log(dct["left"]) + np.arange(dct["nbin"]+1) * dlnr)
-      fout.variables[name+'_r_'+dct["drwt"]][:] = allbins[0:-1]
-      fout.variables[name+'_dr_'+dct["drwt"]][:] = allbins[1:] - allbins[0:-1]
+      fout.variables[name+'_r_'+dct["slct"]][:] = allbins[0:-1]
+      fout.variables[name+'_dr_'+dct["slct"]][:] = allbins[1:] - allbins[0:-1]
     elif dct["lnli"] == 'lin':
       dr = (dct["rght"] - dct["left"]) / dct["nbin"]
-      fout.variables[name+'_r_'+dct["drwt"]][:] = dct["left"] + np.arange(dct["nbin"]) * dr
-      fout.variables[name+'_dr_'+dct["drwt"]][:] = dr
+      fout.variables[name+'_r_'+dct["slct"]][:] = dct["left"] + np.arange(dct["nbin"]) * dr
+      fout.variables[name+'_dr_'+dct["slct"]][:] = dr
     else: raise Exception("lnli should be log or lin")
 
     for vm in dct["moms"]:
@@ -215,11 +239,11 @@ def _output_init(micro, opts, spectra):
       	fout.variables[name+'_'+vm].unit = 'kg of chem species dissolved in cloud droplets (kg of dry air)^-1'
       else:
         assert(type(vm)==int)
-	fout.createVariable(name+'_m'+str(vm), 'd', ('t',name))
-	fout.variables[name+'_m'+str(vm)].unit = 'm^'+str(vm)+' (kg of dry air)^-1'
+	fout.createVariable(name+'_'+dct["drwt"]+'_mom'+str(vm), 'd', ('t',name))
+	fout.variables[name+'_'+dct["drwt"]+'_mom'+str(vm)].unit = 'm^'+str(vm)+' (kg of dry air)^-1'
 
   units = {"z"  : "m",     "t"   : "s",     "r_v"  : "kg/kg", "th_d" : "K", "rhod" : "kg/m3",
-           "p"  : "Pa",    "T"   : "K",     "RH"   : "1"
+           "p"  : "Pa",    "T"   : "K",     "RH"   : "1"    , "RH_libcloud" : "1", "curr_w" : "m/s"
   }
 
   if micro.opts_init.chem_switch:
@@ -253,25 +277,29 @@ def _p_hydro_const_th_rv(z_lev, p_0, th_std, r_v, z_0=0.):
   # hydrostatic pressure assuming constatnt theta and r_v
   return common.p_hydro(z_lev, th_std, r_v, z_0, p_0)
 
-def parcel(dt=.1, z_max=200., w=1., T_0=300., p_0=101300.,
+def parcel(dt=.1, nt=0, z_max=200., z_min=-1., w=1., T_0=300., p_0=101300.,
   r_0=-1., RH_0=-1., #if none specified, the default will be r_0=.022,
   outfile="test.nc",
   pprof="pprof_piecewise_const_rhod",
   outfreq=100, sd_conc=64,
   aerosol = '{"ammonium_sulfate": {"kappa": 0.61, "mean_r": [0.02e-6], "gstdev": [1.4], "n_tot": [60.0e6]}}',
-  out_bin = '{"radii": {"rght": 0.0001, "moms": [0], "drwt": "wet", "nbin": 1, "lnli": "log", "left": 1e-09}}',
+  aerosol_sizes = 'null', 
+  out_bin = '{"radii": {"rght": 0.0001, "moms": [0], "drwt": "wet", "slct": "wet", "nbin": 1, "lnli": "log", "left": 1e-09}}',
   SO2_g = 0., O3_g = 0., H2O2_g = 0., CO2_g = 0., HNO3_g = 0., NH3_g = 0.,
   chem_dsl = False, chem_dsc = False, chem_rct = False,
   chem_rho = 1.8e3,
   sstp_cond = 1,
   sstp_chem = 1,
   wait = 0,
+  top_stop = 0.,
   large_tail = False
 ):
   """
   Args:
     dt      (Optional[float]):    timestep [s]
+    nt      (Optional[int]):      number of timesteps, only for updraft
     z_max   (Optional[float]):    maximum vertical displacement [m]
+    z_min   (Optional[float]):    if non-negative, parcel will descend to z_min after reaching z_max [m]
     w       (Optional[float]):    updraft velocity [m/s]
     T_0     (Optional[float]):    initial temperature [K]
     p_0     (Optional[float]):    initial pressure [Pa]
@@ -284,6 +312,7 @@ def parcel(dt=.1, z_max=200., w=1., T_0=300., p_0=101300.,
                                   valid options are: pprof_const_th_rv, pprof_const_rhod, pprof_piecewise_const_rhod
     wait (Optional[float]):       number of timesteps to run parcel model with vertical velocity=0 at the end of simulation
                                   (added for testing)
+    top_stop (Optional[float]):   time the parcel waits at the top of the cloud before starting the descent [s]
     sd_conc (Optional[int]):      number of moving bins (super-droplets)
 
     aerosol (Optional[json str]): dict of dicts defining aerosol distribution, e.g.:
@@ -297,12 +326,21 @@ def parcel(dt=.1, z_max=200., w=1., T_0=300., p_0=101300.,
                                         n_tot  - lognormal distribution total concentration under standard
                                                  conditions (T=20C, p=1013.25 hPa, rv=0) [m^-3]            (list if multimodal distribution)
 
+    aerosol_sizes (Optional[json str]): dict of dicts defining aerosol by size-concentration pair, e.g.:
+
+                                  {gccn"            : {"kappa": 1.28, "r": [2e-6, 4e-6],      "n_tot": [1e2, 5e1]}}
+
+                                  where kappa  - hygroscopicity parameter (see doi:10.5194/acp-7-1961-2007)
+                                        r -      dry radius [m]                    (list if multimodal distribution)
+                                        n_tot  - total concentration under standard
+                                                 conditions (T=20C, p=1013.25 hPa, rv=0) [m^-3]            (list if multimodal distribution)
+
     large_tail (Optional[bool]) : use more SD to better represent the large tail of the initial aerosol distribution
 
     out_bin (Optional[json str]): dict of dicts defining spectrum diagnostics, e.g.:
 
-                                  {"radii": {"rght": 0.0001,  "moms": [0],          "drwt": "wet", "nbin": 26, "lnli": "log", "left": 1e-09},
-                                   "cloud": {"rght": 2.5e-05, "moms": [0, 1, 2, 3], "drwt": "wet", "nbin": 49, "lnli": "lin", "left": 5e-07}}
+                                  {"radii": {"rght": 0.0001,  "moms": [0],          "drwt": "wet", "slct": "wet", "nbin": 26, "lnli": "log", "left": 1e-09},
+                                   "cloud": {"rght": 2.5e-05, "moms": [0, 1, 2, 3], "drwt": "wet", "slct": "wet", "nbin": 49, "lnli": "lin", "left": 5e-07}}
                                   will generate five output spectra:
                                   - 0-th spectrum moment for 26 bins spaced logarithmically between 0 and 1e-4 m for dry radius
                                   - 0,1,2 & 3-rd moments for 49 bins spaced linearly between .5e-6 and 25e-6 for wet radius
@@ -341,6 +379,7 @@ def parcel(dt=.1, z_max=200., w=1., T_0=300., p_0=101300.,
 
   # parsing json specification of init aerosol spectra
   aerosol = json.loads(opts["aerosol"])
+  aerosol_sizes = json.loads(opts["aerosol_sizes"])
 
   # default water content
   if ((opts["r_0"] < 0) and (opts["RH_0"] < 0)):
@@ -354,13 +393,22 @@ def parcel(dt=.1, z_max=200., w=1., T_0=300., p_0=101300.,
   _arguments_checking(opts, spectra, aerosol)
 
   th_0 = T_0 * (common.p_1000 / p_0)**(common.R_d / common.c_pd)
-  nt = int(z_max / (w * dt))
+  if(z_max > 0):
+    nt_ud = int(round(z_max / (w * dt)))
+    nt_dd = 0
+    if(z_min >= 0):
+      nt_dd = int(round((z_max - z_min) / (w * dt)))
+    nt_topstop = int(round(top_stop / dt))
+    nt = nt_ud + nt_topstop + nt_dd
+
+
   state = {
     "t" : 0, "z" : 0,
     "r_v" : np.array([r_0]), "p" : p_0,
     "th_d" : np.array([common.th_std2dry(th_0, r_0)]),
     "rhod" : np.array([common.rhod(p_0, th_0, r_0)]),
-    "T" : None, "RH" : None
+    "T" : None, "RH" : None, "RH_libcloud" : None,
+    "curr_w" : w
   }
 
   if opts["chem_dsl"] or opts["chem_dsc"] or opts["chem_rct"]:
@@ -370,7 +418,7 @@ def parcel(dt=.1, z_max=200., w=1., T_0=300., p_0=101300.,
   info = { "RH_max" : 0, "libcloud_Git_revision" : libcloud_version,
            "parcel_Git_revision" : parcel_version }
 
-  micro = _micro_init(aerosol, opts, state, info)
+  micro = _micro_init(aerosol, aerosol_sizes, opts, state, info)
 
   with _output_init(micro, opts, spectra) as fout:
     # adding chem state vars
@@ -391,7 +439,16 @@ def parcel(dt=.1, z_max=200., w=1., T_0=300., p_0=101300.,
       # the reasons to use analytic solution:
       # - independent of dt
       # - same as in 2D kinematic model
-      state["z"] += w * dt
+
+      state["curr_w"] = w # for runs with nt specified and z_max = 0
+      if(it < nt_ud):
+        state["curr_w"] = w
+      elif(it < nt_ud + nt_topstop):
+        state["curr_w"] = 0.
+      else:
+        state["curr_w"] = -w
+
+      state["z"] += state["curr_w"] * dt
       state["t"] = it * dt
 
       # pressure
@@ -406,7 +463,7 @@ def parcel(dt=.1, z_max=200., w=1., T_0=300., p_0=101300.,
       elif pprof == "pprof_piecewise_const_rhod":
         # as in Grabowski and Wang 2009 but calculating pressure
         # for rho piecewise constant per each time step
-        state["p"] = _p_hydro_const_rho(w*dt, state["p"], state["rhod"][0])
+        state["p"] = _p_hydro_const_rho(state["curr_w"]*dt, state["p"], state["rhod"][0])
 
       else: raise Exception("pprof should be pprof_const_th_rv, pprof_const_rhod, or pprof_piecewise_const_rhod")
 
@@ -455,6 +512,18 @@ def _arguments_checking(opts, spectra, aerosol):
     raise Exception("temperature should be larger than 0C - microphysics works only for warm clouds")
   elif ((opts["r_0"] >= 0) and (opts["RH_0"] >= 0)):
     raise Exception("both r_0 and RH_0 specified, please use only one")
+  elif ((opts["z_max"] > 0) and (opts["nt"] > 0)):
+    raise Exception("both z_max and nt specified, please use only one")
+  elif ((opts["z_min"] > 0) and (opts["nt"] > 0)):
+    raise Exception("both z_min and nt specified, please use only one")
+  elif ((opts["top_stop"] > 0) and (opts["nt"] > 0)):
+    raise Exception("both top_stop and nt specified, please use only one")
+  elif ((opts["top_stop"] < 0)):
+    raise Exception("top_stop < 0")
+  elif ((opts["z_max"] <= 0) and (opts["nt"] <= 0)):
+    raise Exception("either z_max or nt should be larger than 0")
+  elif ((opts["z_max"] > 0) and (opts["z_min"] >= opts["z_max"])):
+    raise Exception("z_min >= z_max !")
   if opts["w"] < 0:
     raise Exception("vertical velocity should be larger than 0")
 
@@ -494,7 +563,7 @@ def _arguments_checking(opts, spectra, aerosol):
   for name, dct in spectra.iteritems():
     # TODO: check if name is valid netCDF identifier
     # (http://www.unidata.ucar.edu/software/thredds/current/netcdf-java/CDM/Identifiers.html)
-    keys = ["left", "rght", "nbin", "drwt", "lnli", "moms"]
+    keys = ["left", "rght", "nbin", "drwt", "slct", "lnli", "moms"]
     for key in keys:
       if key not in dct:
         raise Exception(">>" + key + "<< is missing in out_bin[" + name + "]")
@@ -509,6 +578,8 @@ def _arguments_checking(opts, spectra, aerosol):
         raise Exception(">>left<< is greater than >>rght<< in out_bin["+ name +"]")
     if dct["drwt"] not in ["dry", "wet"]:
         raise Exception(">>drwt<< key in out_bin["+ name +"] must be either >>dry<< or >>wet<<")
+    if dct["slct"] not in ["dry", "wet"]:
+        raise Exception(">>slct<< key in out_bin["+ name +"] must be either >>dry<< or >>wet<<")
     if dct["lnli"] not in ["lin", "log"]:
         raise Exception(">>lnli<< key in out_bin["+ name +"] must be either >>lin<< or >>log<<")
     if type(dct["nbin"]) != int: