Ice tests fix

long_test/test_ice_SD.py

@@ -27,7 +27,7 @@ def test_ice_SD():
                      depo = True,
                      ice_nucl = True,
                      time_dep_ice_nucl = time_dep,
-                     aerosol = '{"ammonium_sulfate": {"kappa": 0.61, "sol_frac": 0.1, "mean_r": [0.02e-6], "gstdev": [1.4], "n_tot": [60.0e6]}}',
+                     aerosol = '{"ammonium_sulfate": {"kappa": 0.61, "sol_frac": 6.4e-5, "mean_r": [0.5e-6], "gstdev": [1.4], "n_tot": [60.0e6]}}',
                      outfreq = 4000, 
                      out_bin = '{"liq": {"rght": 1, "moms": [0,3], "drwt": "wet", "nbin": 1, "lnli": "lin", "left": 5e-20}}',
                      outfile = outfile,

long_test/test_plot_ice_SD.py

@@ -17,14 +17,15 @@
 
 def plot_profiles(fnc, output_name):
     plt.clf()
-    plt.rcParams.update({'font.size': 14})
+    plt.rcParams.update({'font.size': 13})
     fig, plots = plt.subplots(1, 2, figsize=(15, 5))
     plots[0].set_xlabel('mixing ratio [g/kg]')
-    plots[1].set_xlabel('T [K]')
+    plots[1].set_xlabel('T [C]')
     for ax in plots:
-        ax.set_ylabel('z [m]')
+        ax.set_ylabel('z [km]')
+        ax.grid()
 
-    z = fnc.variables["z"][:]
+    z = fnc.variables["z"][:] / 1000 #km
     r_v = fnc.variables["r_v"][:] * 1000  #g/kg
     r_liq = np.array([i[0] for i in fnc.variables['liq_m3'][:]]) *4/3 * np.pi * common.rho_w * 1000 #g/kg
     r_ice = fnc.variables["ice_mix_ratio"][:] * 1000 #g/kg
@@ -33,26 +34,27 @@ def plot_profiles(fnc, output_name):
     plots[0].plot(r_v, z)
     plots[0].plot(r_liq, z)
     plots[0].plot(r_ice, z)
-    plots[0].legend(['$r_{tot}$', '$r_v$', '$r_{liq}$', '$r_{ice}$'], loc='best')
-    plots[1].plot(fnc.variables["T"][:], z)
+    plots[0].legend(['total', 'water vapor', 'liquid', 'ice'], loc='best', fontsize=12)
+    plots[1].plot(fnc.variables["T"][:] - 273.15, z)
+    plt.suptitle("Homogeneous ice nucleation" if output_name.startswith("ice_SD_plot_hom") else "Heterogeneous ice nucleation")
 
     if not os.path.exists("plots/outputs/"):
         subprocess.call(["mkdir", "plots/outputs/"])
     plt.savefig(os.path.join("plots/outputs/", output_name))
 
 def test_plot_ice_SD():   
-    for (sol_frac, output_name) in [("0.1", "ice_SD_plot_het.svg"), ("1.0", "ice_SD_plot_hom.svg")]:
+    for (mean_r, sol_frac, output_name) in [("0.5e-6", "6.4e-5", "ice_SD_plot_het.pdf"), ("0.02e-6", "1.0", "ice_SD_plot_hom.pdf")]:
         outfile = "onesim_plot.nc"
         parcel(dt=1.,w=1.,sd_conc=100,
             z_max = 5000.0,
-            T_0 = 263.0,
+            T_0 = 265.0,
             RH_0 = 1.,
             scheme = "lgrngn",
             ice_switch=True,
             depo = True,
             ice_nucl=True,
             time_dep_ice_nucl=True,
-            aerosol = f'{{"ammonium_sulfate": {{"kappa": 0.61, "sol_frac": {sol_frac}, "mean_r": [0.02e-6], "gstdev": [1.4], "n_tot": [60.0e6]}}}}',
+            aerosol = f'{{"ammonium_sulfate": {{"kappa": 0.61, "sol_frac": {sol_frac}, "mean_r": [{mean_r}], "gstdev": [1.4], "n_tot": [60.0e6]}}}}',
             outfreq = 100, 
             out_bin= '{"liq": {"rght": 1, "moms": [0,3], "drwt": "wet", "nbin": 1, "lnli": "lin", "left": 5e-20}}',
             outfile=outfile)

long_test/test_plot_schemes_ice.py

@@ -4,7 +4,7 @@
 - lgrngn with time-dependent freezing (as in Arabas et al., 2025)
 - 1-moment bulk (Grabowski, 1999)
 
-It plots the evolution of ice, liquid and vapor mixing ratios for different schemes. 
+It plots the evolution of ice and liquid mixing ratios for different schemes. 
 """
 
 import sys, os
@@ -21,54 +21,55 @@
 
 def run_scheme(scheme, time_dep, outfile):
     args = dict(dt=1., w=1., sd_conc=100,
-            z_max = 6000.0,
+            z_max = 5000.0,
             T_0 = 265.0,
             RH_0 = 1.,
             scheme = scheme,
             ice_switch=True,
             depo = True,
             ice_nucl=True,
             time_dep_ice_nucl=time_dep,
-            aerosol = f'{{"ammonium_sulfate": {{"kappa": 0.61, "sol_frac": 0.1, "mean_r": [0.02e-6], "gstdev": [1.4], "n_tot": [60.0e6]}}}}',
+            aerosol = f'{{"ammonium_sulfate": {{"kappa": 0.61, "sol_frac": 6.4e-5, "mean_r": [0.5e-6], "gstdev": [1.4], "n_tot": [60.0e6]}}}}',
             outfreq = 100, 
             out_bin= '{"liq": {"rght": 1, "moms": [3], "drwt": "wet", "nbin": 1, "lnli": "lin", "left": 5e-20}}',
             outfile=outfile)
     parcel(**args)
     with netcdf.netcdf_file(outfile, 'r') as f:
-        rv = np.array(f.variables['r_v'][:]).squeeze()
         z = np.array(f.variables['z'][:]).squeeze()
-        T = np.array(f.variables['T'][:]).squeeze()
         if scheme.startswith("blk"):
             r_liq = np.array(f.variables['rc'][:]) + np.array(f.variables['rr'][:]).squeeze()
             r_ice = np.array(f.variables['ria'][:]) + np.array(f.variables['rib'][:]).squeeze()
-        else: 
+        else:
             r_liq = np.array(f.variables['liq_m3'][:]).squeeze() *4/3 * np.pi * common.rho_w
             r_ice = np.array(f.variables['ice_mix_ratio'][:]).squeeze()
-        return z, rv, r_liq, r_ice, T
+        return z, r_liq, r_ice
 
 def test_plot_schemes():
     schemes = [("lgrngn", False), ("lgrngn", True), ("blk_1m", None)]
     plt.rcParams.update({'font.size': 14})
-    fig, ax = plt.subplots(1,4, figsize=(16, 6))
+    fig, ax = plt.subplots(1,2, figsize=(9,6))
     for (scheme, time_dep) in schemes:
-        z, rv, r_liq, r_ice, T = run_scheme(scheme, time_dep, f"test_{scheme}.nc")
+        z, r_liq, r_ice = run_scheme(scheme, time_dep, f"test_{scheme}.nc")
         if scheme == "lgrngn" and time_dep == True:
-            l = "Lagrangian time dependent"
+            l = "Lagrangian time-dependent"
         elif time_dep == False:
             l = "Lagrangian singular"
         else:
             l = "1-moment bulk"
-        ax[0].plot(rv*1e3, z, label=l, linestyle ='--' if scheme=='blk_1m' else '-')
-        ax[1].plot(r_liq*1e3, z, label=l, linestyle ='--' if scheme=='blk_1m' else '-')
-        ax[2].plot(r_ice*1e3, z, label=l, linestyle ='--' if scheme=='blk_1m' else '-')
-        ax[3].plot(T, z, label=l, linestyle ='--' if scheme=='blk_1m' else '-')
+        ax[0].plot(r_liq*1e3, z*1e-3, label=l, linestyle ='--' if scheme=='blk_1m' else '-')
+        ax[1].plot(r_ice*1e3, z*1e-3, label=l, linestyle ='--' if scheme=='blk_1m' else '-')
         os.remove(f"test_{scheme}.nc")
 
-    ax[0].set_ylabel("z [m]")
-    ax[0].set_xlabel("$r_v$ [g/kg]")
-    ax[1].set_xlabel("$r_{liq} [g/kg]$")
-    ax[2].set_xlabel("$r_{ice} [g/kg]$")
-    ax[3].set_xlabel("T [K]")
-    ax[0].legend()
-    plt.tight_layout()
-    plt.savefig("plots/outputs/plot_schemes_ice.svg")
+    ax[0].set_ylabel("z [km]")
+    ax[0].set_xlabel("liquid mix. ratio [g/kg]")
+    ax[1].set_xlabel("ice mix. ratio [g/kg]")
+    handles, labels = ax[0].get_legend_handles_labels()
+    fig.legend(
+        handles, labels,
+        loc="upper center",
+        bbox_to_anchor=(0.5, 1.0),
+        ncol=len(labels),
+        frameon=False,
+    )
+    fig.tight_layout(rect=[0, 0, 1, 0.9]) 
+    plt.savefig("plots/outputs/plot_schemes_ice.pdf", bbox_inches="tight")