To merge

README.md

@@ -12,7 +12,7 @@ For more information on the [libcloudph++](http://libcloudphxx.igf.fuw.edu.pl/)
   [project repository](https://github.com/igfuw/libcloudphxx) and 
   the [project documentation](http://www.geosci-model-dev.net/8/1677/2015/).
 
-The parcel model is written in Python 2.7. 
+The parcel model was written in Python 2.7 and now works with Python 3.# . 
 
 # installation
 

long_test/test_timestep.py

@@ -1,5 +1,6 @@
 import sys
 sys.path.insert(0, "../")
+#sys.path.insert(0, "path to /parcel/plots/comparison")
 sys.path.insert(0, "./")
 sys.path.insert(0, "plots/comparison/")
 
@@ -15,33 +16,35 @@
 import pdb
 
 """
-This set of tests checks how the timestep of the simulation affects 
-the maximum supersaturation (RH) and the concentration of the activated 
-particles (N).                                      
-The expected result is to see nearly constant RH and N for small timesteps 
-(.001 - .03 for this setup) and then decrease of RH and increase of N 
-for bigger timesteps. 
+This set of tests checks how the timestep of the simulation affects
+the maximum supersaturation (RH) and the concentration of the activated
+particles (N).
+The expected result is to see nearly constant RH and N for small timesteps
+(.001 - .03 for this setup) and then decrease of RH and increase of N
+for bigger timesteps.
 """
 
 Dt_list = [1e-3, 2e-3, 4e-3, 8e-3, 1e-2, 2e-2, 4e-2, 1e-1, 2e-1, 1.]
 
-# runs all simulations 
+# runs all simulations
 # returns data with values of RH_max and N at the end of simulations
 # keep the netcdf alive for all tests
 @pytest.fixture(scope="module")
 def data(request):
-    # initial values                                        
+    # initial values
     RH_init = .99999
     T_init  = 280.
     p_init  = 100000.
     r_init  = common.eps * RH_init * common.p_vs(T_init) / (p_init - RH_init * common.p_vs(T_init))
-    # lists to store RH_max and N at the end of the simulation from each test run 
+    # lists to store RH_max and N at the end of the simulation from each test run
     RH_list = []
     N_list  = []
 
     for dt in Dt_list:
+
         print("\nt time step", dt)
         outfile_nc = "timesteptest_dt=" + str(dt) + ".nc" 
+
         parcel(dt=dt, outfreq = int(100/dt),   outfile = outfile_nc,\
                 w = 1., T_0 = T_init, p_0 = p_init, r_0 = r_init, z_max = 200, \
                 sd_conc = 1000, \
@@ -52,9 +55,9 @@ def data(request):
         f_out  = netcdf.netcdf_file(outfile_nc, "r")
         RH_max = f_out.RH_max
         N_end  = f_out.variables["radii_m0"][-1,0] # concentration of drops > 1e-6 m
-    
-        RH_list.append((RH_max - 1)*100)  # [%]                                      
-        N_list.append(N_end / 1e6)        # [1/mg]           
+
+        RH_list.append((RH_max - 1)*100)  # [%]
+        N_list.append(N_end / 1e6)        # [1/mg]
 
     data = {"RH" : RH_list, "N" : N_list, "dt" : Dt_list}
 
@@ -68,35 +71,35 @@ def removing_files():
 
 def test_timestep_eps(data, eps=0.01, dt_lim=0.01):
     """
-    checking if the results obtained from simulations with different timesteps   
-    do not differ from the referential one (the one with the smallest timestep) 
-    more than eps times 
-    (Unitill we think of a better convergence test, the check is done for the 
+    checking if the results obtained from simulations with different timesteps
+    do not differ from the referential one (the one with the smallest timestep)
+    more than eps times
+    (Unitill we think of a better convergence test, the check is done for the
     smallest timesteps. This is done in order to avoid too big epsilon.
     """
     for var, val in data.items():
         # check for RH and N
         if var in ["RH", "N"]:
             # for simulations with small timesteps
             for idx in range(len(data["dt"])):
-                if data["dt"][idx] < dt_lim: 
+                if data["dt"][idx] < dt_lim:
                     # assert that the results are close to the one with the smallest timestep
-                    assert np.isclose(val[idx], val[0], atol=0, rtol=eps), str(val[idx]) + str(val[0]) 
+                    assert np.isclose(val[idx], val[0], atol=0, rtol=eps), str(val[idx]) + str(val[0])
 
 
 @pytest.mark.parametrize("dt", Dt_list)
 def test_timestep_diff(data, dt, eps=2.e-4):
     """
-    checking if the results are close to the referential ones 
+    checking if the results are close to the referential ones
     (stored in refdata folder)
     """
 
     filename = "timesteptest_dt=" + str(dt) + ".nc"
     f_test = netcdf.netcdf_file(filename, "r")
-    f_ref  = netcdf.netcdf_file(os.path.join("long_test/refdata", filename), "r")
+    f_ref  = netcdf.netcdf_file(os.path.join("/home/piotr/Piotr/IGF/parcel/long_test/refdata", filename), "r")
     for var in ["t", "z", "th_d", "T", "p", "r_v", "rhod"]:
          assert np.isclose(f_test.variables[var][:], f_ref.variables[var][:], atol=0, rtol=eps).all(), "differs e.g. " + str(var) + "; max(ref diff) = " + str(np.where(f_ref.variables[var][:] != 0., abs((f_test.variables[var][:]-f_ref.variables[var][:])/f_ref.variables[var][:]), 0.).max())
-        
+
 
 def test_timestep_plot(data):
     timestep_plot(data, output_folder="plots/outputs")

parcel.py

@@ -5,6 +5,7 @@
 #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, "path_to/lib/python3/dist-packages")
 # TEMP TODO TEMP TODO !!!
 
 from argparse import ArgumentParser, RawTextHelpFormatter
@@ -155,17 +156,19 @@ def _output_bins(fout, t, micro, opts, spectra):
   for dim, dct in spectra.items():
     for bin in range(dct["nbin"]):
       if dct["drwt"] == 'wet':
-        micro.diag_wet_rng(
-          fout.variables[dim+"_r_wet"][bin],
-          fout.variables[dim+"_r_wet"][bin] + fout.variables[dim+"_dr_wet"][bin]
-        )
+         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':
-        micro.diag_dry_rng(
-          fout.variables[dim+"_r_dry"][bin],
-          fout.variables[dim+"_r_dry"][bin] + fout.variables[dim+"_dr_dry"][bin]
-        )
+          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")
 
+
       for vm in dct["moms"]:
         if type(vm) == int:
           # calculating moments
@@ -175,6 +178,7 @@ def _output_bins(fout, t, micro, opts, spectra):
             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())
+
         else:
           # calculate chemistry
           micro.diag_chem(_Chem_a_id[vm])
@@ -214,9 +218,14 @@ def _output_init(micro, opts, spectra):
       	fout.createVariable(name+'_'+vm, 'd', ('t',name))
       	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'
+
+          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('number_of_rc_m0', 'd', ('t'))
+          fout.createVariable('number_of_rc_m1', 'd', ('t'))
+          fout.variables['number_of_rc_m0'].unit = 'm^'+str(vm)+' (kg of dry air)^-1'
+          fout.variables['number_of_rc_m1'].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"
@@ -242,6 +251,12 @@ def _save_attrs(fout, dictnr):
     setattr(fout, var, val)
 
 def _output(fout, opts, micro, state, rec, spectra):
+  micro.diag_rw_ge_rc()
+  micro.diag_wet_mom(0)
+  fout.variables['number_of_rc_m0'][int(rec)] = np.frombuffer(micro.outbuf())
+  micro.diag_rw_ge_rc()
+  micro.diag_wet_mom(1)
+  fout.variables['number_of_rc_m1'][int(rec)] = np.frombuffer(micro.outbuf())
   _output_bins(fout, rec, micro, opts, spectra)
   _output_save(fout, state, rec)
 
@@ -336,7 +351,8 @@ def parcel(dt=.1, z_max=200., w=1., T_0=300., p_0=101300.,
   args, _, _, _ = inspect.getargvalues(inspect.currentframe())
   opts = dict()
   for k in args:
-    opts[k] = locals()[k]
+      opts[k] = locals()[k]
+
 
   # parsing json specification of output spectra
   spectra = json.loads(opts["out_bin"])
@@ -350,7 +366,7 @@ def parcel(dt=.1, z_max=200., w=1., T_0=300., p_0=101300.,
     r_0 = .022
   # water coontent specified with RH
   if ((opts["r_0"] < 0) and (opts["RH_0"] >= 0)):
-    r_0 = common.eps * opts["RH_0"] * common.p_vs(T_0) / (p_0 - opts["RH_0"] * common.p_vs(T_0))
+      r_0 = common.eps * opts["RH_0"] * common.p_vs(T_0) / (p_0 - opts["RH_0"] * common.p_vs(T_0))
 
   # sanity checks for arguments
   _arguments_checking(opts, spectra, aerosol)
@@ -490,8 +506,8 @@ def _arguments_checking(opts, spectra, aerosol):
       if gstdev <= 0:
         raise Exception("standard deviation should be > 0 for aerosol[" + name + "]")
     # necessary?
-      if gstdev == 1.:
-        raise Exception("standard deviation should be != 1 to avoid monodisperse distribution for aerosol[" + name + "]")
+     # if gstdev == 1.:
+       # raise Exception("standard deviation should be != 1 to avoid monodisperse distribution for aerosol[" + name + "]")
 
   for name, dct in spectra.items():
     # TODO: check if name is valid netCDF identifier

plots/one_simulat/init_spectrum_plot.py

@@ -27,9 +27,12 @@ def plot_init_spectrum(data, outfolder):
     import Gnuplot
 
     # size distribution parameters from Kreidenweis 2003
-    n_tot   = 566e6
-    mean_r  = 0.04e-6
-    gstdev  = 2
+    n_tot   = 125e6#, 15e6]#566e6
+    mean_r  = 0.011e-6#, 0.14e-6]#0.04e-6
+    gstdev  = 1.2#, 1.75]#2
+    n_tot2   = 65e6#566e6
+    mean_r2  = 0.06e-6#0.04e-6
+    gstdev2  = 1.7#2
 
     # from ncdf file attributes read out_bin parameters as a dictionary ...
     out_bin = eval(getattr(data, "out_bin"))
@@ -49,9 +52,10 @@ def plot_init_spectrum(data, outfolder):
     # variables for plotting theoretical solution
     radii = np.logspace(-3, 1, 100) * 1e-6
     theor = np.empty(radii.shape)
+    theor2 = np.empty(radii.shape)
     for it in range(radii.shape[0]):
         theor[it] = fn.log10_size_of_lnr(n_tot, mean_r, math.log(radii[it], 10), gstdev)
-
+        theor2[it] = fn.log10_size_of_lnr(n_tot2, mean_r2, math.log(radii[it], 10), gstdev2)
     g = Gnuplot.Gnuplot()
     g('set term svg dynamic enhanced')
     g('reset')
@@ -64,9 +68,10 @@ def plot_init_spectrum(data, outfolder):
     g('set yrange [0:800]')
 
     theory_r = Gnuplot.PlotItems.Data(radii * 2 * 1e6,  theor * 1e-6, with_="lines", title="theory")
+    theory_r2 = Gnuplot.PlotItems.Data(radii * 2 * 1e6,  theor2 * 1e-6, with_="lines", title="theory2")
     plot     = Gnuplot.PlotItems.Data(rd    * 2 * 1e6,  model * 1e-6, with_="steps", title="model" )
 
-    g.plot(theory_r, plot)
+    g.plot(theory_r, theory_r2, plot)
 
 def main():
     # copy options from chem_conditions ...
@@ -83,7 +88,9 @@ def main():
     opts_dict['sd_conc']  = 1024 * 44
     opts_dict['outfreq']  = 1
 
-    opts_dict['out_bin'] = '{"drad": {"rght": 1e-6, "left": 1e-10, "drwt": "dry", "lnli": "log", "nbin": 26, "moms": [0,3]}}'
+    opts_dict['out_bin'] = '{"drad": {"rght": 2.5e-06, "left": 5e-09, "drwt": "dry", "lnli": "log", "nbin": 26, "moms": [0,1,3]}}'
+
+                            #{"rght": 2.5e-05, "moms": [0,1,3], "drwt": "wet", "nbin": 26, "lnli": "log", "left": 5e-07}}
 
     # run parcel
     parcel(**opts_dict)

plots/one_simulat/profiles_plot.py

@@ -22,10 +22,10 @@ def plot_profiles(fnc, output_folder="../outputs"):
 
     plots[0].set_xlabel('p [hPa]')
 
-    plots[1].ticklabel_format(useOffset=False) 
+    plots[1].ticklabel_format(useOffset=False)
     plots[1].set_xlabel('th_d [K]')
     plots[2].set_xlabel('T [K]')
-    plots[3].set_xlabel('kappa(rho_d :)) [kg/m3]')  
+    plots[3].set_xlabel('kappa(rho_d :)) [kg/m3]')
     plots[4].set_xlabel('rv [g/kg]')
     plots[5].set_xlabel('RH')
 
@@ -39,18 +39,18 @@ def plot_profiles(fnc, output_folder="../outputs"):
     plots[3].plot(fnc.variables["rhod"][:]       , z)
     plots[4].plot(fnc.variables["r_v"][:] * 1000 , z)
     plots[5].plot(
-        fnc.variables["RH"][:]                     , z, 
+        fnc.variables["RH"][:]                     , z,
         [fnc.variables["RH"][:].max()] * z.shape[0], z
         )
-   
+
     if not os.path.exists(output_folder):
         subprocess.call(["mkdir", output_folder])
     plt.savefig(os.path.join(output_folder, "plot_profiles_onesim.svg"))
 
 def main(dt=1):
-    # running parcel model for different ways to solve for pressure  ...            
+    # running parcel model for different ways to solve for pressure  ...
     outfile = "onesim_plot.nc"
-    parcel(dt=dt, outfreq = 10, outfile=outfile)
+    parcel(dt=dt, outfreq = 10, w=5, outfile=outfile)
     fnc = netcdf.netcdf_file(outfile)
     plot_profiles(fnc)
     fnc.close()