diff --git a/.dockerignore b/.dockerignore index e725e276f..42155c2ca 100644 --- a/.dockerignore +++ b/.dockerignore @@ -25,3 +25,4 @@ scenarios/4_chamber/out scenarios/5_coag_brownian/out !scenarios/6_camp scenarios/6_camp/out +!scenarios/7_freezing diff --git a/.github/workflows/tchem.yml b/.github/workflows/tchem.yml.temporarily_off similarity index 100% rename from .github/workflows/tchem.yml rename to .github/workflows/tchem.yml.temporarily_off diff --git a/.gitignore b/.gitignore index e0ec53304..f8b4ede53 100644 --- a/.gitignore +++ b/.gitignore @@ -17,5 +17,7 @@ scenarios/3_condense/spec/ scenarios/3_condense/temp/ scenarios/3_condense/out/ scenarios/4_chamber/out/ +scenarios/7_freezing/out/ +scenarios/7_freezing *.o *.mod diff --git a/CMakeLists.txt b/CMakeLists.txt index 788097ae3..30e13b04f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -198,12 +198,13 @@ if(ENABLE_SUNDIALS) endif() add_test(test_emission ${CMAKE_BINARY_DIR}/test_run/run_test_directory.sh emission) add_test(test_fractal ${CMAKE_BINARY_DIR}/test_run/run_test_directory.sh fractal) +add_test(test_freezing ${CMAKE_BINARY_DIR}/test_run/run_test_directory.sh freezing) add_test(test_loss ${CMAKE_BINARY_DIR}/test_run/run_test_directory.sh loss) -add_test(test_nucleate ${CMAKE_BINARY_DIR}/test_run/run_test_directory.sh nucleate) add_test(test_mixing_state ${CMAKE_BINARY_DIR}/test_run/run_test_directory.sh mixing_state) if(ENABLE_MOSAIC) add_test(test_mosaic ${CMAKE_BINARY_DIR}/test_run/run_test_directory.sh mosaic) endif() +add_test(test_nucleate ${CMAKE_BINARY_DIR}/test_run/run_test_directory.sh nucleate) if(ENABLE_MPI) add_test(test_parallel ${CMAKE_BINARY_DIR}/test_run/run_test_directory.sh parallel) endif() @@ -213,6 +214,7 @@ if (ENABLE_TCHEM) add_test(test_tchem ${CMAKE_BINARY_DIR}/test_run/run_test_directory.sh tchem) endif() add_test(test_weighting ${CMAKE_BINARY_DIR}/test_run/run_test_directory.sh weighting) +#>>>>>>> 2341ef410d6f49f3169b8461b5fa8c89dbd3c7a2 ###################################################################### # partmc library @@ -234,6 +236,7 @@ add_library(partmclib src/aero_state.F90 src/integer_varray.F90 src/netcdf.F90 src/aero_info.F90 src/aero_info_array.F90 src/nucleate.F90 src/condense.F90 src/fractal.F90 src/chamber.F90 src/camp_interface.F90 src/photolysis.F90 src/sys.F90 + src/ice_nucleation.F90 src/tchem_interface.F90 src/aero_component.F90 ${SUNDIALS_SRC} ${GSL_SRC} ${TCHEM_SRC} @@ -286,6 +289,22 @@ add_executable(test_bidisperse_extract target_link_libraries(test_bidisperse_extract ${NETCDF_LIBS}) +###################################################################### +# test_freezing_extract + +add_executable(test_freezing_extract + test/freezing/test_freezing_extract.F90 src/getopt.F90) + +target_link_libraries(test_freezing_extract ${NETCDF_LIBS} partmclib) + +###################################################################### +# test_freezing_theoretical + +add_executable(test_freezing_theoretical + test/freezing/test_freezing_theoretical.F90) + +target_link_libraries(test_freezing_theoretical partmclib) + ###################################################################### # test_nucleate_ode @@ -405,6 +424,13 @@ add_executable(urban_plume_process scenarios/1_urban_plume/urban_plume_process.F90) target_link_libraries(urban_plume_process partmclib) +###################################################################### +# scenarios/7_freezing/freezing_process + +add_executable(freezing_process + scenarios/7_freezing/freezing_process.F90 src/getopt.F90) +target_link_libraries(freezing_process partmclib) + ###################################################################### # scenarios/4_chamber/chamber_process diff --git a/doc/freezing/freezing.tex b/doc/freezing/freezing.tex new file mode 100644 index 000000000..9e5b95812 --- /dev/null +++ b/doc/freezing/freezing.tex @@ -0,0 +1,90 @@ +\documentclass{article} +\usepackage{amsmath} +\usepackage[margin=2cm]{geometry} +\usepackage{longtable} +\usepackage{comment} + +\begin{document} + +\section{Symbol list} +\label{sec:symbol-list} + +\newcommand{\rr}{\raggedright} +\newcommand{\tn}{\tabularnewline\hline} +\renewcommand{\arraystretch}{1.5} +\begin{longtable}{|l|p{5.5cm}|l|l|p{4.5cm}|} +\hline \textbf{Symbol} & \textbf{Meaning} & \textbf{Type} & \textbf{Units} & \textbf{Reference} \tn +\hline \endhead +\(i\) & species index. & index & discrete & \tn +\( A_j \) & The total surface area of the INP. & function & $m^2$ & \tn +\( A_j^\text{(i)} \) & Surface area of \( i \)th species. & function & $m^2$ & \tn +\( a_{\text{w}} \) & Water activity of the droplet in which the INP is immersed. & function & 1 & eq.~\ref{eq:4}\tn +\( a_\text{w}^\text{ice} \) & Water activity in equilibrium with ice. & function & 1 & eq.~\ref{eq:5} \tn +\(a_\text{INAS}\) & INAS parameter for singular scheme. & constant & $K^{-1}$ & Niemand et al., 2012 \tn +\(b_\text{INAS}\) & INAS parameter for singular scheme. & constant & 1 & Niemand et al., 2012 \tn +\( \Delta t \) & Time interval. & function & $s$ & \tn +\( e_\text{s} \) & Saturated vapor pressure with respect to water. & function & Pa & eq.~\ref{eq:2} \tn +\( e_\text{s}^\text{ice} \) & Saturated vapor pressure with respect to ice. & function & Pa & eq.~\ref{eq:3} \tn +\( H \) & Relative humidity with respect to water. & dynamics & 1 & \tn +\( J_{\text{het}} \) & A constant heterogeneous freezing rate for all particles (for the constant nucleation rate scheme only) & constant & $m^{-2}s^{-1}$ & \tn +\( J_{\text{het}}^{(i)} \) & Heterogeneous freezing rate coefficients of the \( i \)th species. & function & $m^{-2}s^{-1}$ & eq.~\ref{eq:6} \tn +\( m_i, c_i \) & ABIFM parameters of \( i^\text{th} \) species. & constant & 1 & Knopf et al., 2013 \tn +\(n_0\) & An unit constant (=1 $m^{-2}$). & constant & $m^{-2}$ & \tn +\( N_\text{s} \) & The total number of species covering the surface of the INP. & constant & 1 & \tn +\( P_{\text{frz,}j} (\Delta t)\) & The freezing probability of $j^\text{th}$ particle within a time interval of $\Delta t$. & function & 1 & eq.~\ref{eq:1},eq.~\ref{eq:7} \tn +\(p_j\) & A random number uniformly distributed in [0,1] regarding the sampling of $j^\text{th}$ particle. & function & 1 & eq.~\ref{eq:8} \tn +\(T\) & Current environmental temperature & prescribed & $K$ & \tn +\(T_0\) & Water freezing temperature & constant & K & \verb+env_state%temp+ \tn +\(T_{\text{m,}j}\) & The freezing temperature of the $j^\text{th}$ particle. & function & K & eq.~\ref{eq:9} \tn + +\end{longtable} + +\begin{comment} +\( J_{\text{het}}^{(i)} \) & Heterogeneous freezing rate coefficients of the \( i \)th species. \\ + \( A_i \) & Surface area of \( i \)th species. \\ + \( A \) & The total surface area of the INP. \\ + \( \Delta t \) & Time interval. \\ + \( a_{\text{w}} \) & Water activity of the droplet in which the INP is immersed. \\ + \( a_\text{w}^\text{ice} \) & Water activity in equilibrium with ice. \\ + \( RH \) & Relative humidity with respect to water. \\ + \( e_\text{s} \) & Saturated vapor pressure with respect to water. \\ + \( e_\text{s}^\text{ice} \) & Saturated vapor pressure with respect to ice. \\ + \( m_i, c_i \) & ABIFM parameters of \( i^\text{th} \) species. \\ + \( N_\text{s} \) & The total number of species covering the surface of the INP. \\ + \hline +\end{comment} + +\section{Equations} +\label{sec:equations} + +\subsection{Time-dependent scheme with a constant nucleation rate} +\begin{equation} + P_{\text{frz,}j}(\Delta t) = 1 - \exp\left(J_{\text{het}} \cdot A_j \cdot \Delta t \right) \label{eq:1} +\end{equation} + +\subsection{Time-dependent ABIFM scheme} +\begin{equation} + \begin{split} + e_\text{s}(T) &= \exp \left\{ 54.842763 - \frac{6763.22}{T} - 4.21 * \ln{T} + 0.000367T \right.\\ + &+ \left. \tanh \left[{0.0415 \cdot (T - 218.8) \cdot (53.878 - \frac{1331.22}{T} - 9.44523 \ln{T} + 0.014025T)}\right] \right\} \label{eq:2} + \end{split} +\end{equation} +\begin{equation} + e_\text{s}^\text{ice} (T) = \exp \left( 9.550426 - \frac{5723.265}{T} + 3.53068 \cdot \ln{T} - 0.00728332T\right). \label{eq:3} +\end{equation} +\begin{align} + a_{\text{w}} &= 1 \label{eq:4} \\ + a_\text{w}^\text{ice} &= \frac{e_\text{s}^\text{ice}(T)}{e_\text{s}(T)} \label{eq:5} \\ + \log_{10}J_{\text{het}}^{(i)} &= m_{\text{i}} \cdot (a_{\text{w}} - a_\text{w}^\text{ice}) + c_i \label{eq:6} \\ + P_{\text{frz,}j}(\Delta t) &= 1 - \exp\left[ -\sum_{i=1}^{N_{\text{s}}} (A_j^\text{(i)} \cdot J_{\text{het}}^{(i)}) \cdot \Delta t \right] \label{eq:7} +\end{align} + +\subsection{Singular scheme} +\begin{equation} + p_j = random() \in [0,1] \label{eq:8} +\end{equation} + +\begin{equation} + T_{\text{m,}j} = T_0 + \frac{1}{a_\text{INAS}} \left[\ln{\left( \frac{\ln{(1-p_j)} + \exp{\left(-A_j n_0 \exp{\left(-a_\text{INAS} T_0 + b_\text{INAS}\right)}\right)}}{-A_j n_0} \right)} - b_\text{INAS}\right] \label{eq:9} +\end{equation} +\end{document} diff --git a/scenarios/1_urban_plume/aero_data.dat b/scenarios/1_urban_plume/aero_data.dat index 43422f091..26fcd7e32 100644 --- a/scenarios/1_urban_plume/aero_data.dat +++ b/scenarios/1_urban_plume/aero_data.dat @@ -1,21 +1,21 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -SO4 1800 0 96d-3 0.65 -NO3 1800 0 62d-3 0.65 -Cl 2200 0 35.5d-3 1.28 -NH4 1800 0 18d-3 0.65 -MSA 1800 0 95d-3 0.53 -ARO1 1400 0 150d-3 0.1 -ARO2 1400 0 150d-3 0.1 -ALK1 1400 0 140d-3 0.1 -OLE1 1400 0 140d-3 0.1 -API1 1400 0 184d-3 0.1 -API2 1400 0 184d-3 0.1 -LIM1 1400 0 200d-3 0.1 -LIM2 1400 0 200d-3 0.1 -CO3 2600 0 60d-3 0.53 -Na 2200 0 23d-3 1.28 -Ca 2600 0 40d-3 0.53 -OIN 2600 0 1d-3 0.1 -OC 1000 0 1d-3 0.001 -BC 1800 0 1d-3 0 -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) +SO4 1800 0 96d-3 0.65 0 0 +NO3 1800 0 62d-3 0.65 0 0 +Cl 2200 0 35.5d-3 1.28 0 0 +NH4 1800 0 18d-3 0.65 0 0 +MSA 1800 0 95d-3 0.53 0 0 +ARO1 1400 0 150d-3 0.1 0 0 +ARO2 1400 0 150d-3 0.1 0 0 +ALK1 1400 0 140d-3 0.1 0 0 +OLE1 1400 0 140d-3 0.1 0 0 +API1 1400 0 184d-3 0.1 0 0 +API2 1400 0 184d-3 0.1 0 0 +LIM1 1400 0 200d-3 0.1 0 0 +LIM2 1400 0 200d-3 0.1 0 0 +CO3 2600 0 60d-3 0.53 0 0 +Na 2200 0 23d-3 1.28 0 0 +Ca 2600 0 40d-3 0.53 0 0 +OIN 2600 0 1d-3 0.1 0 0 +OC 1000 0 1d-3 0.001 0 0 +BC 1800 0 1d-3 0 0 0 +H2O 1000 0 18d-3 0 0 0 diff --git a/scenarios/1_urban_plume/urban_plume.spec b/scenarios/1_urban_plume/urban_plume.spec index 4f9290bbd..a44a95cab 100644 --- a/scenarios/1_urban_plume/urban_plume.spec +++ b/scenarios/1_urban_plume/urban_plume.spec @@ -42,6 +42,7 @@ do_condensation no # whether to do condensation (yes/no) do_mosaic yes # whether to do MOSAIC (yes/no) do_optical yes # whether to compute optical props (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to use time) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/scenarios/2_urban_plume2/aero_data.dat b/scenarios/2_urban_plume2/aero_data.dat index cd98619d5..26fcd7e32 100644 --- a/scenarios/2_urban_plume2/aero_data.dat +++ b/scenarios/2_urban_plume2/aero_data.dat @@ -1,21 +1,21 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -SO4 1800 0 96d-3 0.65 -NO3 1800 0 62d-3 0.65 -Cl 2200 0 35.5d-3 1.28 -NH4 1800 0 18d-3 0.65 -MSA 1800 0 95d-3 0.53 -ARO1 1400 0 150d-3 0.1 -ARO2 1400 0 150d-3 0.1 -ALK1 1400 0 140d-3 0.1 -OLE1 1400 0 140d-3 0.1 -API1 1400 0 184d-3 0.1 -API2 1400 0 184d-3 0.1 -LIM1 1400 0 200d-3 0.1 -LIM2 1400 0 200d-3 0.1 -CO3 2600 0 60d-3 0.53 -Na 2200 0 23d-3 1.28 -Ca 2600 0 40d-3 0.53 -OIN 2600 0 1d-3 0.1 -OC 1000 0 1d-3 0.001 -BC 1800 0 1d-3 0 -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) +SO4 1800 0 96d-3 0.65 0 0 +NO3 1800 0 62d-3 0.65 0 0 +Cl 2200 0 35.5d-3 1.28 0 0 +NH4 1800 0 18d-3 0.65 0 0 +MSA 1800 0 95d-3 0.53 0 0 +ARO1 1400 0 150d-3 0.1 0 0 +ARO2 1400 0 150d-3 0.1 0 0 +ALK1 1400 0 140d-3 0.1 0 0 +OLE1 1400 0 140d-3 0.1 0 0 +API1 1400 0 184d-3 0.1 0 0 +API2 1400 0 184d-3 0.1 0 0 +LIM1 1400 0 200d-3 0.1 0 0 +LIM2 1400 0 200d-3 0.1 0 0 +CO3 2600 0 60d-3 0.53 0 0 +Na 2200 0 23d-3 1.28 0 0 +Ca 2600 0 40d-3 0.53 0 0 +OIN 2600 0 1d-3 0.1 0 0 +OC 1000 0 1d-3 0.001 0 0 +BC 1800 0 1d-3 0 0 0 +H2O 1000 0 18d-3 0 0 0 diff --git a/scenarios/2_urban_plume2/urban_plume2.spec b/scenarios/2_urban_plume2/urban_plume2.spec index 6a22fa81e..5e1ecd27e 100644 --- a/scenarios/2_urban_plume2/urban_plume2.spec +++ b/scenarios/2_urban_plume2/urban_plume2.spec @@ -42,6 +42,7 @@ do_condensation no # whether to do condensation (yes/no) do_mosaic yes # whether to do MOSAIC (yes/no) do_optical yes # whether to compute optical props (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to use time) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/scenarios/3_condense/aero_data.dat b/scenarios/3_condense/aero_data.dat index 66960e1c5..d0cd57cfa 100644 --- a/scenarios/3_condense/aero_data.dat +++ b/scenarios/3_condense/aero_data.dat @@ -1,21 +1,21 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -SO4 1800 1 96d-3 0 -NO3 1800 1 62d-3 0 -Cl 2200 1 35.5d-3 0 -NH4 1800 1 18d-3 0 -MSA 1800 0 95d-3 0.53 -ARO1 1400 0 150d-3 0.1 -ARO2 1400 0 150d-3 0.1 -ALK1 1400 0 140d-3 0.1 -OLE1 1400 0 140d-3 0.1 -API1 1400 0 184d-3 0.1 -API2 1400 0 184d-3 0.1 -LIM1 1400 0 200d-3 0.1 -LIM2 1400 0 200d-3 0.1 -CO3 2600 1 60d-3 0 -Na 2200 1 23d-3 0 -Ca 2600 1 40d-3 0 -OIN 2600 0 1d-3 0.1 -OC 1000 0 1d-3 0.1 -BC 1800 0 1d-3 0 -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) +SO4 1800 1 96d-3 0 0 0 +NO3 1800 1 62d-3 0 0 0 +Cl 2200 1 35.5d-3 0 0 0 +NH4 1800 1 18d-3 0 0 0 +MSA 1800 0 95d-3 0.53 0 0 +ARO1 1400 0 150d-3 0.1 0 0 +ARO2 1400 0 150d-3 0.1 0 0 +ALK1 1400 0 140d-3 0.1 0 0 +OLE1 1400 0 140d-3 0.1 0 0 +API1 1400 0 184d-3 0.1 0 0 +API2 1400 0 184d-3 0.1 0 0 +LIM1 1400 0 200d-3 0.1 0 0 +LIM2 1400 0 200d-3 0.1 0 0 +CO3 2600 1 60d-3 0 0 0 +Na 2200 1 23d-3 0 0 0 +Ca 2600 1 40d-3 0 0 0 +OIN 2600 0 1d-3 0.1 0 0 +OC 1000 0 1d-3 0.1 0 0 +BC 1800 0 1d-3 0 0 0 +H2O 1000 0 18d-3 0 0 0 diff --git a/scenarios/3_condense/cond_template.spec b/scenarios/3_condense/cond_template.spec index 699fb6725..64235625e 100644 --- a/scenarios/3_condense/cond_template.spec +++ b/scenarios/3_condense/cond_template.spec @@ -34,6 +34,7 @@ do_condensation yes # whether to do condensation (yes/no) do_init_equilibrate yes # whether to initially equilibrate water (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to use time) allow_doubling no # whether to allow doubling (yes/no) diff --git a/scenarios/4_chamber/aero_data.dat b/scenarios/4_chamber/aero_data.dat index 43422f091..26fcd7e32 100644 --- a/scenarios/4_chamber/aero_data.dat +++ b/scenarios/4_chamber/aero_data.dat @@ -1,21 +1,21 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -SO4 1800 0 96d-3 0.65 -NO3 1800 0 62d-3 0.65 -Cl 2200 0 35.5d-3 1.28 -NH4 1800 0 18d-3 0.65 -MSA 1800 0 95d-3 0.53 -ARO1 1400 0 150d-3 0.1 -ARO2 1400 0 150d-3 0.1 -ALK1 1400 0 140d-3 0.1 -OLE1 1400 0 140d-3 0.1 -API1 1400 0 184d-3 0.1 -API2 1400 0 184d-3 0.1 -LIM1 1400 0 200d-3 0.1 -LIM2 1400 0 200d-3 0.1 -CO3 2600 0 60d-3 0.53 -Na 2200 0 23d-3 1.28 -Ca 2600 0 40d-3 0.53 -OIN 2600 0 1d-3 0.1 -OC 1000 0 1d-3 0.001 -BC 1800 0 1d-3 0 -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) +SO4 1800 0 96d-3 0.65 0 0 +NO3 1800 0 62d-3 0.65 0 0 +Cl 2200 0 35.5d-3 1.28 0 0 +NH4 1800 0 18d-3 0.65 0 0 +MSA 1800 0 95d-3 0.53 0 0 +ARO1 1400 0 150d-3 0.1 0 0 +ARO2 1400 0 150d-3 0.1 0 0 +ALK1 1400 0 140d-3 0.1 0 0 +OLE1 1400 0 140d-3 0.1 0 0 +API1 1400 0 184d-3 0.1 0 0 +API2 1400 0 184d-3 0.1 0 0 +LIM1 1400 0 200d-3 0.1 0 0 +LIM2 1400 0 200d-3 0.1 0 0 +CO3 2600 0 60d-3 0.53 0 0 +Na 2200 0 23d-3 1.28 0 0 +Ca 2600 0 40d-3 0.53 0 0 +OIN 2600 0 1d-3 0.1 0 0 +OC 1000 0 1d-3 0.001 0 0 +BC 1800 0 1d-3 0 0 0 +H2O 1000 0 18d-3 0 0 0 diff --git a/scenarios/4_chamber/chamber.spec b/scenarios/4_chamber/chamber.spec index be35e8ed8..2300c21f8 100644 --- a/scenarios/4_chamber/chamber.spec +++ b/scenarios/4_chamber/chamber.spec @@ -49,6 +49,7 @@ coag_kernel brown # coagulation kernel do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 7 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/scenarios/5_coag_brownian/1_run.sh b/scenarios/5_coag_brownian/1_run.sh old mode 100644 new mode 100755 diff --git a/scenarios/5_coag_brownian/aero_data.dat b/scenarios/5_coag_brownian/aero_data.dat index 43422f091..26fcd7e32 100644 --- a/scenarios/5_coag_brownian/aero_data.dat +++ b/scenarios/5_coag_brownian/aero_data.dat @@ -1,21 +1,21 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -SO4 1800 0 96d-3 0.65 -NO3 1800 0 62d-3 0.65 -Cl 2200 0 35.5d-3 1.28 -NH4 1800 0 18d-3 0.65 -MSA 1800 0 95d-3 0.53 -ARO1 1400 0 150d-3 0.1 -ARO2 1400 0 150d-3 0.1 -ALK1 1400 0 140d-3 0.1 -OLE1 1400 0 140d-3 0.1 -API1 1400 0 184d-3 0.1 -API2 1400 0 184d-3 0.1 -LIM1 1400 0 200d-3 0.1 -LIM2 1400 0 200d-3 0.1 -CO3 2600 0 60d-3 0.53 -Na 2200 0 23d-3 1.28 -Ca 2600 0 40d-3 0.53 -OIN 2600 0 1d-3 0.1 -OC 1000 0 1d-3 0.001 -BC 1800 0 1d-3 0 -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) +SO4 1800 0 96d-3 0.65 0 0 +NO3 1800 0 62d-3 0.65 0 0 +Cl 2200 0 35.5d-3 1.28 0 0 +NH4 1800 0 18d-3 0.65 0 0 +MSA 1800 0 95d-3 0.53 0 0 +ARO1 1400 0 150d-3 0.1 0 0 +ARO2 1400 0 150d-3 0.1 0 0 +ALK1 1400 0 140d-3 0.1 0 0 +OLE1 1400 0 140d-3 0.1 0 0 +API1 1400 0 184d-3 0.1 0 0 +API2 1400 0 184d-3 0.1 0 0 +LIM1 1400 0 200d-3 0.1 0 0 +LIM2 1400 0 200d-3 0.1 0 0 +CO3 2600 0 60d-3 0.53 0 0 +Na 2200 0 23d-3 1.28 0 0 +Ca 2600 0 40d-3 0.53 0 0 +OIN 2600 0 1d-3 0.1 0 0 +OC 1000 0 1d-3 0.001 0 0 +BC 1800 0 1d-3 0 0 0 +H2O 1000 0 18d-3 0 0 0 diff --git a/scenarios/5_coag_brownian/example.spec b/scenarios/5_coag_brownian/example.spec index fff0dae5c..90e63fb4b 100644 --- a/scenarios/5_coag_brownian/example.spec +++ b/scenarios/5_coag_brownian/example.spec @@ -41,6 +41,7 @@ coag_kernel brown # coagulation kernel do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 66 # random initialization (0 to use time) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/scenarios/6_camp/camp.spec b/scenarios/6_camp/camp.spec index e51383447..c1db38bab 100644 --- a/scenarios/6_camp/camp.spec +++ b/scenarios/6_camp/camp.spec @@ -39,6 +39,7 @@ do_coagulation no # whether to do coagulation (yes/no) do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to use time) allow_doubling no # whether to allow doubling (yes/no) diff --git a/scenarios/6_camp/monarch_mod37/custom_species.json b/scenarios/6_camp/monarch_mod37/custom_species.json index df8673d47..0ef9e68ca 100644 --- a/scenarios/6_camp/monarch_mod37/custom_species.json +++ b/scenarios/6_camp/monarch_mod37/custom_species.json @@ -28,7 +28,9 @@ "molecular weight [kg mol-1]" : 0.1, "density [kg m-3]" : 2650.0, "num_ions" : 0, - "kappa" : 0.1 + "kappa" : 0.1, + "abifm_m": 0, + "abifm_c": 0 }, { "monarch name" : "lumped low-density dust species", @@ -39,7 +41,9 @@ "molecular weight [kg mol-1]" : 0.1, "density [kg m-3]" : 2650.0, "num_ions" : 0, - "kappa" : 0.1 + "kappa" : 0.1, + "abifm_m": 0, + "abifm_c": 0 }, { "monarch name" : "lumped sea salt species", @@ -50,7 +54,9 @@ "molecular weight [kg mol-1]" : 0.1, "density [kg m-3]" : 2160.0, "num_ions" : 0, - "kappa" : 0.53 + "kappa" : 0.53, + "abifm_m": 0, + "abifm_c": 0 }, { "name" : "BC_phob", @@ -61,7 +67,9 @@ "molecular weight [kg mol-1]" : 0.1, "description" : "hydrophobic black carbon", "num_ions" : 0, - "kappa" : 0 + "kappa" : 0, + "abifm_m": 0, + "abifm_c": 0 }, { "name" : "BC_phil", @@ -72,7 +80,9 @@ "molecular weight [kg mol-1]" : 0.1, "description" : "hydrophilic black carbon", "num_ions" : 0, - "kappa" : 0.001 + "kappa" : 0.001, + "abifm_m": 0, + "abifm_c": 0 }, { "name" : "other_PM", @@ -83,7 +93,9 @@ "density [kg m-3]" : 1800.0, "description" : "unspecified particulate matter FIXME", "num_ions" : 0, - "kappa" : 0 + "kappa" : 0, + "abifm_m": 0, + "abifm_c": 0 }, { "name" : "other_other_PM", @@ -94,7 +106,9 @@ "density [kg m-3]" : 1800.0, "description" : "second unspecified particulate matter species FIXME", "num_ions" : 0, - "kappa" : 0 + "kappa" : 0, + "abifm_m": 0, + "abifm_c": 0 } ] } diff --git a/scenarios/6_camp/monarch_mod37/tsigaridis_2_product_SOA_scheme/species.json b/scenarios/6_camp/monarch_mod37/tsigaridis_2_product_SOA_scheme/species.json index e46bc5511..72bfcc6b0 100644 --- a/scenarios/6_camp/monarch_mod37/tsigaridis_2_product_SOA_scheme/species.json +++ b/scenarios/6_camp/monarch_mod37/tsigaridis_2_product_SOA_scheme/species.json @@ -58,6 +58,8 @@ "description" : "lumped hydrophobic particulate matter", "num_ions" : 0, "kappa" : 0.0, + "abifm_m": 0, + "abifm_c": 0, "note" : "Using C30H54 for molecular weight. TODO find best surrogate" }, { @@ -70,6 +72,8 @@ "description" : "First isoprene SOA species in 2-product scheme", "num_ions" : 0, "kappa" : 0.1, + "abifm_m": 0, + "abifm_c": 0, "note" : [ "Using ketopropanoic acid for molecular weight. TODO find best surrogate", "TODO update SIMPOL parameters based on MW of new surrogate" @@ -85,6 +89,8 @@ "description" : "Second isoprene SOA species in 2-product scheme", "num_ions" : 0, "kappa" : 0.1, + "abifm_m": 0, + "abifm_c": 0, "note" : [ "Using oxalic acid for molecular weight. TODO find best surrogate", "TODO update SIMPOL parameters based on MW of new surrogate" @@ -100,6 +106,8 @@ "description" : "First monoterpene SOA species in 2-product scheme", "num_ions" : 0, "kappa" : 0.1, + "abifm_m": 0, + "abifm_c": 0, "note" : [ "Using 2-hydroxy-3-isopropyl-6-methyl-cyclohexanone for molecular weight", "TODO find best surrogate", @@ -116,6 +124,8 @@ "description" : "Second monoterpene SOA species in 2-product scheme", "num_ions" : 0, "kappa" : 0.1, + "abifm_m": 0, + "abifm_c": 0, "note" : [ "Using 2-methyl-5-carboxy-2,4-hexodienoic acid. TODO find best surrogate", "TODO update SIMPOL parameters based on MW of new surrogate" diff --git a/src/aero_data.F90 b/src/aero_data.F90 index 5ff634c97..d0fe85384 100644 --- a/src/aero_data.F90 +++ b/src/aero_data.F90 @@ -70,6 +70,10 @@ module pmc_aero_data real(kind=dp), allocatable :: molec_weight(:) !> Length \c aero_data_n_spec(aero_data), kappas (1). real(kind=dp), allocatable :: kappa(:) + !> Length \c aero_data_n_spec(aero_data), abifm_m (1). + real(kind=dp), allocatable :: abifm_m(:) + !> Length \c aero_data_n_spec(aero_data), abifm_c (1). + real(kind=dp), allocatable :: abifm_c(:) !> Length \c aero_data_n_source(aero_data), source names. character(len=AERO_SOURCE_NAME_LEN), allocatable :: source_name(:) !> Length \c aero_data_n_weight_classes, weight class names. @@ -455,6 +459,8 @@ subroutine spec_file_read_aero_data(file, aero_data) !! - molecular weight (real, unit kg/mol) !! - kappa hygroscopicity parameter (real, dimensionless) - if !! zero, then inferred from the ions value + !! - abifm_m: The m parameter for ABIFM algorithm + !! - abifm_c: The c parameter for ABIFM algorithm !! !! This specifies both which species are to be recognized as !! aerosol consituents, as well as their physical properties. For @@ -481,9 +487,9 @@ subroutine spec_file_read_aero_data(file, aero_data) ! check the data size n_species = size(species_data, 1) - if (.not. ((size(species_data, 2) == 4) .or. (n_species == 0))) then + if (.not. ((size(species_data, 2) == 6) .or. (n_species == 0))) then call die_msg(428926381, 'each line in ' // trim(file%name) & - // ' should contain exactly 5 values') + // ' should contain exactly 7 values') end if ! allocate and copy over the data @@ -494,12 +500,16 @@ subroutine spec_file_read_aero_data(file, aero_data) call ensure_integer_array_size(aero_data%num_ions, n_species) call ensure_real_array_size(aero_data%molec_weight, n_species) call ensure_real_array_size(aero_data%kappa, n_species) + call ensure_real_array_size(aero_data%abifm_m, n_species) + call ensure_real_array_size(aero_data%abifm_c, n_species) do i = 1,n_species aero_data%name(i) = species_name(i)(1:AERO_NAME_LEN) aero_data%density(i) = species_data(i,1) aero_data%num_ions(i) = nint(species_data(i,2)) aero_data%molec_weight(i) = species_data(i,3) aero_data%kappa(i) = species_data(i,4) + aero_data%abifm_m(i) = species_data(i,5) + aero_data%abifm_c(i) = species_data(i,6) call assert_msg(232362742, & (aero_data%num_ions(i) == 0) .or. (aero_data%kappa(i) == 0d0), & "ions and kappa both non-zero for species " & @@ -572,6 +582,8 @@ integer function pmc_mpi_pack_size_aero_data(val) + pmc_mpi_pack_size_integer_array(val%num_ions) & + pmc_mpi_pack_size_real_array(val%molec_weight) & + pmc_mpi_pack_size_real_array(val%kappa) & + + pmc_mpi_pack_size_real_array(val%abifm_m) & + + pmc_mpi_pack_size_real_array(val%abifm_c) & + pmc_mpi_pack_size_string_array(val%source_name) & + pmc_mpi_pack_size_string_array(val%weight_class_name) & + pmc_mpi_pack_size_fractal(val%fractal) @@ -602,6 +614,8 @@ subroutine pmc_mpi_pack_aero_data(buffer, position, val) call pmc_mpi_pack_integer_array(buffer, position, val%num_ions) call pmc_mpi_pack_real_array(buffer, position, val%molec_weight) call pmc_mpi_pack_real_array(buffer, position, val%kappa) + call pmc_mpi_pack_real_array(buffer, position, val%abifm_m) + call pmc_mpi_pack_real_array(buffer, position, val%abifm_c) call pmc_mpi_pack_string_array(buffer, position, val%source_name) call pmc_mpi_pack_string_array(buffer, position, val%weight_class_name) call pmc_mpi_pack_fractal(buffer, position, val%fractal) @@ -635,6 +649,8 @@ subroutine pmc_mpi_unpack_aero_data(buffer, position, val) call pmc_mpi_unpack_integer_array(buffer, position, val%num_ions) call pmc_mpi_unpack_real_array(buffer, position, val%molec_weight) call pmc_mpi_unpack_real_array(buffer, position, val%kappa) + call pmc_mpi_unpack_real_array(buffer, position, val%abifm_m) + call pmc_mpi_unpack_real_array(buffer, position, val%abifm_c) call pmc_mpi_unpack_string_array(buffer, position, val%source_name) call pmc_mpi_unpack_string_array(buffer, position, val%weight_class_name) call pmc_mpi_unpack_fractal(buffer, position, val%fractal) @@ -886,6 +902,8 @@ subroutine aero_data_output_netcdf(aero_data, ncid) !! weights of aerosol species !! - \b aero_kappa (unit kg/mol, dim \c aero_species): hygroscopicity !! parameters of aerosol species + !! - \b aero_abifm_m (dim \c aero_species): m parameter of ABIFM + !! - \b aero_abifm_c (dim \c aero_species): c parameter of ABIFM !! - \b fractal parameters, see \ref output_format_fractal !! !! See also: @@ -918,6 +936,12 @@ subroutine aero_data_output_netcdf(aero_data, ncid) call pmc_nc_write_real_1d(ncid, aero_data%kappa, & "aero_kappa", (/ dimid_aero_species /), unit="1", & long_name="hygroscopicity parameters (kappas) of aerosol species") + call pmc_nc_write_real_1d(ncid, aero_data%abifm_m, & + "aero_abifm_m", (/ dimid_aero_species /), unit="1", & + long_name="m parameter of ABIFM") + call pmc_nc_write_real_1d(ncid, aero_data%abifm_c, & + "aero_abifm_c", (/ dimid_aero_species /), unit="1", & + long_name="c parameter of ABIFM") call pmc_nc_write_integer(ncid, aero_data%i_water, & "aero_i_water", long_name="Index of aerosol water or " & // "0 if water does not exist.") @@ -967,6 +991,8 @@ subroutine aero_data_input_netcdf(aero_data, ncid) call pmc_nc_read_integer_1d(ncid, aero_data%num_ions, "aero_num_ions") call pmc_nc_read_real_1d(ncid, aero_data%molec_weight, "aero_molec_weight") call pmc_nc_read_real_1d(ncid, aero_data%kappa, "aero_kappa") + call pmc_nc_read_real_1d(ncid, aero_data%abifm_m, "aero_abifm_m") + call pmc_nc_read_real_1d(ncid, aero_data%abifm_c, "aero_abifm_c") call pmc_nc_check(nf90_inq_varid(ncid, "aero_species", & varid_aero_species)) @@ -1090,6 +1116,8 @@ subroutine aero_data_initialize(aero_data, camp_core) allocate(aero_data%num_ions(num_spec)) allocate(aero_data%molec_weight(num_spec)) allocate(aero_data%kappa(num_spec)) + allocate(aero_data%abifm_m(num_spec)) + allocate(aero_data%abifm_c(num_spec)) allocate(aero_data%camp_particle_spec_id(num_spec)) ! Assume no aerosol water @@ -1127,6 +1155,18 @@ subroutine aero_data_initialize(aero_data, camp_core) call die_msg(944207343, "Missing kappa for aerosol species " & // spec_names(i_spec)%string) end if + prop_name = "abifm_m" + if (.not. property_set%get_real(prop_name, & + aero_data%abifm_m(i_spec))) then + call die_msg(944207345, "Missing abifm_m for aerosol species " & + // spec_names(i_spec)%string) + end if + prop_name = "abifm_c" + if (.not. property_set%get_real(prop_name, & + aero_data%abifm_c(i_spec))) then + call die_msg(944207346, "Missing abifm_c for aerosol species " & + // spec_names(i_spec)%string) + end if prop_name = "PartMC name" if (property_set%get_string(prop_name, str_val)) then if (str_val == "H2O") then diff --git a/src/aero_particle.F90 b/src/aero_particle.F90 index acd60af9f..fb6712af2 100644 --- a/src/aero_particle.F90 +++ b/src/aero_particle.F90 @@ -56,6 +56,14 @@ module pmc_aero_particle real(kind=dp) :: least_create_time !> Last time a constituent was created (s). real(kind=dp) :: greatest_create_time + !> Ice-phase flag. + logical :: frozen + !> Immersion freezing temperature (K). + real(kind=dp) :: imf_temperature + !> Ice density (kg m^{-3}). + real(kind=dp) :: den_ice + !> Ice shape. + real(kind=dp) :: ice_shape_phi !> True number of primary particles contributing to particle. integer :: n_primary_parts end type aero_particle_t @@ -96,6 +104,10 @@ subroutine aero_particle_shift(aero_particle_from, aero_particle_to) aero_particle_to%least_create_time = aero_particle_from%least_create_time aero_particle_to%greatest_create_time = & aero_particle_from%greatest_create_time + aero_particle_to%frozen = aero_particle_from%frozen + aero_particle_to%imf_temperature = aero_particle_from%imf_temperature + aero_particle_to%den_ice = aero_particle_from%den_ice + aero_particle_to%ice_shape_phi = aero_particle_from%ice_shape_phi aero_particle_to%n_primary_parts = aero_particle_from%n_primary_parts end subroutine aero_particle_shift @@ -131,6 +143,10 @@ subroutine aero_particle_zero(aero_particle, aero_data) allocate(aero_particle%component(0)) aero_particle%least_create_time = 0d0 aero_particle%greatest_create_time = 0d0 + aero_particle%frozen = .false. + aero_particle%imf_temperature = 0d0 + aero_particle%den_ice = const%nan + aero_particle%ice_shape_phi = const%nan aero_particle%n_primary_parts = 0 end subroutine aero_particle_zero @@ -899,15 +915,19 @@ end function aero_particle_crit_diameter !> Coagulate two particles together to make a new one. The new !> particle will not have its ID set. subroutine aero_particle_coagulate(aero_particle_1, & - aero_particle_2, aero_particle_new) + aero_particle_2, aero_particle_new, aero_data) !> First particle. type(aero_particle_t), intent(in) :: aero_particle_1 !> Second particle. type(aero_particle_t), intent(in) :: aero_particle_2 + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data !> Result particle. type(aero_particle_t), intent(inout) :: aero_particle_new + + real(kind=dp) :: ice_vol_1, ice_vol_2 integer :: n_comp_1, n_comp_2, n_comp_1_new, n_comp_2_new, i type(aero_component_t), allocatable :: new_aero_component(:) integer :: n_swbands @@ -977,6 +997,35 @@ subroutine aero_particle_coagulate(aero_particle_1, & aero_particle_new%greatest_create_time = & max(aero_particle_1%greatest_create_time, & aero_particle_2%greatest_create_time) + aero_particle_new%frozen = aero_particle_1%frozen .OR. & + aero_particle_2%frozen + + ! Assumption: If two particles coagulate, the overall freezing temperature + ! is determined by the particle with higher freezing temperature. + aero_particle_new%imf_temperature = max(aero_particle_1%imf_temperature, & + aero_particle_2%imf_temperature) + + if (aero_particle_new%frozen) then + ice_vol_1 = aero_particle_1%vol(aero_data%i_water) + ice_vol_2 = aero_particle_2%vol(aero_data%i_water) + + if (aero_particle_1%frozen .and. aero_particle_2%frozen) then + ice_vol_1 = aero_particle_1%vol(aero_data%i_water) * & + const%water_density / aero_particle_1%den_ice + ice_vol_2 = aero_particle_2%vol(aero_data%i_water) * & + const%water_density / aero_particle_2%den_ice + aero_particle_new%den_ice = (aero_particle_1%den_ice * ice_vol_1 + & + aero_particle_2%den_ice * ice_vol_2) / (ice_vol_1 + ice_vol_2) + else if(aero_particle_1%frozen) then + aero_particle_new%den_ice = aero_particle_1%den_ice + else + aero_particle_new%den_ice = aero_particle_2%den_ice + end if + aero_particle_new%ice_shape_phi = 1d0 + else + aero_particle_new%den_ice = const%nan + end if + aero_particle_new%n_primary_parts = aero_particle_1%n_primary_parts & + aero_particle_2%n_primary_parts @@ -1024,6 +1073,10 @@ integer function pmc_mpi_pack_size_aero_particle(val) + pmc_mpi_pack_size_integer(aero_particle_n_components(val)) & + pmc_mpi_pack_size_real(val%least_create_time) & + pmc_mpi_pack_size_real(val%greatest_create_time) & + + pmc_mpi_pack_size_logical(val%frozen) & + + pmc_mpi_pack_size_real(val%imf_temperature) & + + pmc_mpi_pack_size_real(val%den_ice) & + + pmc_mpi_pack_size_real(val%ice_shape_phi) & + pmc_mpi_pack_size_integer(val%n_primary_parts) do i = 1,aero_particle_n_components(val) @@ -1067,6 +1120,10 @@ subroutine pmc_mpi_pack_aero_particle(buffer, position, val) end do call pmc_mpi_pack_real(buffer, position, val%least_create_time) call pmc_mpi_pack_real(buffer, position, val%greatest_create_time) + call pmc_mpi_pack_logical(buffer, position, val%frozen) + call pmc_mpi_pack_real(buffer, position, val%imf_temperature) + call pmc_mpi_pack_real(buffer, position, val%den_ice) + call pmc_mpi_pack_real(buffer, position, val%ice_shape_phi) call pmc_mpi_pack_integer(buffer, position, val%n_primary_parts) call assert(810223998, position - prev_position & <= pmc_mpi_pack_size_aero_particle(val)) @@ -1110,6 +1167,10 @@ subroutine pmc_mpi_unpack_aero_particle(buffer, position, val) end do call pmc_mpi_unpack_real(buffer, position, val%least_create_time) call pmc_mpi_unpack_real(buffer, position, val%greatest_create_time) + call pmc_mpi_unpack_logical(buffer, position, val%frozen) + call pmc_mpi_unpack_real(buffer, position, val%imf_temperature) + call pmc_mpi_unpack_real(buffer, position, val%den_ice) + call pmc_mpi_unpack_real(buffer, position, val%ice_shape_phi) call pmc_mpi_unpack_integer(buffer, position, val%n_primary_parts) call assert(287447241, position - prev_position & <= pmc_mpi_pack_size_aero_particle(val)) diff --git a/src/aero_state.F90 b/src/aero_state.F90 index 90720111d..e86f741f1 100644 --- a/src/aero_state.F90 +++ b/src/aero_state.F90 @@ -1345,6 +1345,33 @@ function aero_state_num_concs_by_source(aero_state, aero_data) end function aero_state_num_concs_by_source +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Returns the frozen fraction (unitless) for whole aerosol population.. + !> frozen fraction = number concentration of frozen particles / + !> total number concentration. + real(kind=dp) function aero_state_frozen_fraction(aero_state, aero_data) + !> Aerosol state. + type(aero_state_t), intent(in) :: aero_state + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + + real(kind=dp) :: particle_num_concs(aero_state_n_part(aero_state)) + logical :: particle_frozen(aero_state_n_part(aero_state)) + + integer :: i_part + + particle_num_concs = aero_state_num_concs(aero_state, aero_data) + + do i_part = 1,aero_state_n_part(aero_state) + particle_frozen(i_part) = aero_state%apa%particle(i_part)%frozen + end do + + aero_state_frozen_fraction = sum(pack(particle_num_concs, particle_frozen)) & + / sum(particle_num_concs) + + end function aero_state_frozen_fraction + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Returns the number concentration of a given weight group and class. @@ -2643,6 +2670,10 @@ subroutine aero_state_output_netcdf(aero_state, ncid, aero_data, & integer :: aero_water_hyst_leg(aero_state_n_part(aero_state)) real(kind=dp) :: aero_num_conc(aero_state_n_part(aero_state)) integer(kind=8) :: aero_id(aero_state_n_part(aero_state)) + integer :: aero_frozen(aero_state_n_part(aero_state)) + real(kind=dp) :: aero_imf_temperature(aero_state_n_part(aero_state)) + real(kind=dp) :: aero_ice_density(aero_state_n_part(aero_state)) + real(kind=dp) :: aero_ice_shape_phi(aero_state_n_part(aero_state)) real(kind=dp) :: aero_least_create_time(aero_state_n_part(aero_state)) real(kind=dp) :: aero_greatest_create_time(aero_state_n_part(aero_state)) integer(kind=8) :: aero_removed_id( & @@ -2802,6 +2833,16 @@ subroutine aero_state_output_netcdf(aero_state, ncid, aero_data, & = aero_state_particle_num_conc(aero_state, & aero_state%apa%particle(i_part), aero_data) aero_id(i_part) = aero_state%apa%particle(i_part)%id + if (aero_state%apa%particle(i_part)%frozen) then + aero_frozen(i_part) = 1 + else + aero_frozen(i_part) = 0 + end if + aero_imf_temperature(i_part) & + = aero_state%apa%particle(i_part)%imf_temperature + aero_ice_density(i_part) = aero_state%apa%particle(i_part)%den_ice + aero_ice_shape_phi(i_part) & + = aero_state%apa%particle(i_part)%ice_shape_phi aero_least_create_time(i_part) & = aero_state%apa%particle(i_part)%least_create_time aero_greatest_create_time(i_part) & @@ -2863,6 +2904,18 @@ subroutine aero_state_output_netcdf(aero_state, ncid, aero_data, & call pmc_nc_write_integer64_1d(ncid, aero_id, & "aero_id", (/ dimid_aero_particle /), & long_name="unique ID number of each aerosol particle") + call pmc_nc_write_integer_1d(ncid, aero_frozen, & + "aero_frozen", (/ dimid_aero_particle /), & + long_name="flag indicating ice phase of each aerosol particle") + call pmc_nc_write_real_1d(ncid, aero_imf_temperature, & + "aero_imf_temperature", (/ dimid_aero_particle /), & + long_name="immersion freezing temperature (Singular)") + call pmc_nc_write_real_1d(ncid, aero_ice_density, & + "aero_ice_density", (/ dimid_aero_particle /), & + long_name="Ice density if the particle nucleates to ice, -9999 indicates the particle is not an ice.") + call pmc_nc_write_real_1d(ncid, aero_ice_shape_phi, & + "aero_ice_shape_phi", (/ dimid_aero_particle /), & + long_name="Ice shape parameter phi") call pmc_nc_write_real_1d(ncid, aero_least_create_time, & "aero_least_create_time", (/ dimid_aero_particle /), unit="s", & long_name="least creation time of each aerosol particle", & @@ -3057,6 +3110,10 @@ subroutine aero_state_input_netcdf(aero_state, ncid, aero_data) real(kind=dp), allocatable :: aero_core_vol(:) integer, allocatable :: aero_water_hyst_leg(:) real(kind=dp), allocatable :: aero_num_conc(:) + integer, allocatable :: aero_frozen(:) + real(kind=dp), allocatable :: aero_imf_temperature(:) + real(kind=dp), allocatable :: aero_ice_density(:) + real(kind=dp), allocatable :: aero_ice_shape_phi(:) integer(kind=8), allocatable :: aero_id(:) real(kind=dp), allocatable :: aero_least_create_time(:) real(kind=dp), allocatable :: aero_greatest_create_time(:) @@ -3123,6 +3180,14 @@ subroutine aero_state_input_netcdf(aero_state, ncid, aero_data) "aero_num_conc") call pmc_nc_read_integer64_1d(ncid, aero_id, & "aero_id") + call pmc_nc_read_integer_1d(ncid, aero_frozen, & + "aero_frozen") + call pmc_nc_read_real_1d(ncid, aero_imf_temperature, & + "aero_imf_temperature") + call pmc_nc_read_real_1d(ncid, aero_ice_density, & + "aero_ice_density", must_be_present=.true.) + call pmc_nc_read_real_1d(ncid, aero_ice_shape_phi, & + "aero_ice_shape_phi", must_be_present=.true.) call pmc_nc_read_real_1d(ncid, aero_least_create_time, & "aero_least_create_time") call pmc_nc_read_real_1d(ncid, aero_greatest_create_time, & @@ -3174,6 +3239,14 @@ subroutine aero_state_input_netcdf(aero_state, ncid, aero_data) end if aero_particle%water_hyst_leg = aero_water_hyst_leg(i_part) aero_particle%id = aero_id(i_part) + if (aero_frozen(i_part) == 1) then + aero_particle%frozen = .true. + else + aero_particle%frozen = .false. + end if + aero_particle%imf_temperature = aero_imf_temperature(i_part) + aero_particle%den_ice = aero_ice_density(i_part) + aero_particle%ice_shape_phi = aero_ice_shape_phi(i_part) aero_particle%least_create_time = aero_least_create_time(i_part) aero_particle%greatest_create_time = aero_greatest_create_time(i_part) diff --git a/src/coagulation.F90 b/src/coagulation.F90 index 8925bc9b7..891efb6fd 100644 --- a/src/coagulation.F90 +++ b/src/coagulation.F90 @@ -438,7 +438,7 @@ subroutine sample_source_particle(aero_state, aero_data, env_state, & if (pmc_random() < prob_coag) then n_avg = n_avg + 1 call aero_particle_coagulate(source_particle, & - aero_state%apa%particle(i_part), source_particle) + aero_state%apa%particle(i_part), source_particle, aero_data) vol_sq = vol_sq + aero_state%apa%particle(i_part)%vol**2 if (i_samp <= n_samp_remove) then num_conc_i = aero_weight_array_num_conc(aero_state%awa, & @@ -547,7 +547,7 @@ subroutine coag_target_with_source(aero_state, aero_data, bt, ct, & old_num_conc_target = aero_weight_array_num_conc(aero_state%awa, & aero_state%apa%particle(target_part), aero_data) call aero_particle_coagulate(aero_state%apa%particle(target_part), & - source_particle, aero_state%apa%particle(target_part)) + source_particle, aero_state%apa%particle(target_part), aero_data) aero_state%apa%particle(target_part)%id = target_id ! assign to a randomly chosen group new_group = aero_weight_array_rand_group(aero_state%awa, cc, & @@ -893,7 +893,7 @@ subroutine coagulate_weighting(pt1, pt2, ptc, c1, c2, cc, aero_data, & ! create a new particle and set its ID if (create_new) then - call aero_particle_coagulate(pt1, pt2, ptc) + call aero_particle_coagulate(pt1, pt2, ptc, aero_data) call aero_particle_set_weight(ptc, new_group, cc) if (remove_1 .and. (.not. id_1_lost)) then ptc%id = pt1%id diff --git a/src/constants.F90 b/src/constants.F90 index 477725848..041187ca0 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -63,6 +63,12 @@ module pmc_constants real(kind=dp) :: air_std_press = 101325d0 !> Dynamic viscosity of air (kg m^{-1} s^{-1}). real(kind=dp) :: air_dyn_visc = 1.78d-5 + !> Reference ice density (kg m^{-3}). + real(kind=dp) :: reference_ice_density = 920d0 + !> Immersion freezing water mass ratio threshold. + real(kind=dp) :: imf_water_threshold = 1d-2 + !> NaN value. + real(kind=dp) :: nan = -9999d0 end type const_t !> Fixed variable for accessing the constant's values. diff --git a/src/env_state.F90 b/src/env_state.F90 index 5daf17491..db02d2b55 100644 --- a/src/env_state.F90 +++ b/src/env_state.F90 @@ -643,6 +643,58 @@ subroutine env_state_input_netcdf(env_state, ncid) end subroutine env_state_input_netcdf + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> compute saturated vapor pressure (units : Pa) with respective to water + !> Formula (10) from [Murphy & Koop, 2004] (https://doi.org/10.1256/qj.04.94) + real(kind=dp) function env_state_saturated_vapor_pressure_wrt_water(T) + + !> Temperature (K) + real(kind=dp), intent(in) :: T + + real(kind=dp) :: tmp + + + call warn_assert_msg(571128376, (T > 123) .and. (T < 332), & + "The environment temperature is less then 123K or larger than "& + "332K, the subroutine env_state_saturated_vapor_pressure_wrt_water"& + " isn't applicable") + + tmp = 54.842763 & + - 6763.22 / T & + - 4.210 * log(T) & + + 0.000367 * T & + + tanh( 0.0415 * (T - 218.8)) & + * (53.878 - 1331.22 / T - 9.44523 * log(T) + 0.014025 * T) + env_state_saturated_vapor_pressure_wrt_water = exp(tmp) + + end function env_state_saturated_vapor_pressure_wrt_water + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Compute saturated vapor pressure (units : Pa) with respective to ice + !> Formula (7) from [Murphy & Koop, 2004] (https://doi.org/10.1256/qj.04.94) + real(kind=dp) function env_state_saturated_vapor_pressure_wrt_ice(T) + + !> Temperature (K) + real(kind=dp), intent(in) :: T + + real(kind=dp) :: tmp + + call warn_assert_msg(482130832, T > 110, & + "The environment temperature is less then 110K, "& + "the subroutine env_state_saturated_vapor_pressure_wrt_ice"& + " isn't applicable") + + tmp = 9.550426 & + - 5723.265 / T & + + 3.53068 * log(T) & + - 0.00728332 * T + env_state_saturated_vapor_pressure_wrt_ice = exp(tmp) + + end function env_state_saturated_vapor_pressure_wrt_ice + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! end module pmc_env_state diff --git a/src/extract_aero_particles.F90 b/src/extract_aero_particles.F90 index cda2f2115..936dc3689 100644 --- a/src/extract_aero_particles.F90 +++ b/src/extract_aero_particles.F90 @@ -83,6 +83,7 @@ program extract_aero_particles trim(aero_data%name(i_spec)), ' mass (kg) - density = ', & aero_data%density(i_spec), ' (kg/m^3)' end do + write(*,'(a,i2,a)') ' column ', aero_data_n_spec(aero_data) + 4 + 1,': frozen flag' call open_file_write(out_filename, out_unit) do i_part = 1,aero_state_n_part(aero_state) @@ -97,6 +98,11 @@ program extract_aero_particles aero_particle_species_mass(aero_state%apa%particle(i_part), & i_spec, aero_data) end do + if (aero_state%apa%particle(i_part)%frozen) then + write(out_unit, '(a,i1)', advance='no') ' ', 1 + else + write(out_unit, '(a,i1)', advance='no') ' ', 0 + end if write(out_unit, '(a)') '' end do call close_file(out_unit) diff --git a/src/ice_nucleation.F90 b/src/ice_nucleation.F90 new file mode 100644 index 000000000..0cfcd913a --- /dev/null +++ b/src/ice_nucleation.F90 @@ -0,0 +1,515 @@ +! Copyright (C) 2024 University of Illinois at Urbana-Champaign +! Licensed under the GNU General Public License version 2 or (at your +! option) any later version. See the file COPYING for details. + +!> \file +! The pmc_ice_nucleation module. + +module pmc_ice_nucleation + use pmc_aero_state + use pmc_env_state + use pmc_aero_data + use pmc_util + use pmc_aero_particle + use pmc_constants + use pmc_rand + + implicit none + + !> Type code for an undefined or invalid immersion freezing scheme. + integer, parameter :: IMMERSION_FREEZING_SCHEME_INVALID = 0 + !> Type code for constant ice nucleation rate (J_het) immersion freezing scheme. + integer, parameter :: IMMERSION_FREEZING_SCHEME_CONST = 1 + !> Type code for the singular (INAS) immersion freezing scheme. + integer, parameter :: IMMERSION_FREEZING_SCHEME_SINGULAR = 2 + !> Type code for the ABIFM immersion freezing scheme. + integer, parameter :: IMMERSION_FREEZING_SCHEME_ABIFM = 3 + +contains + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Main subroutine for immersion freezing simulation. + subroutine ice_nucleation_immersion_freezing(aero_state, aero_data, & + env_state, del_t, immersion_freezing_scheme_type, & + freezing_rate, do_freezing_naive) + + !> Aerosol state. + type(aero_state_t), intent(inout) :: aero_state + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + !> Environment state. + type(env_state_t), intent(inout) :: env_state + !> Total time to integrate. + real(kind=dp), intent(in) :: del_t + !> Immersion freezing scheme type. + integer, intent(in) :: immersion_freezing_scheme_type + !> Freezing rate (only used for the constant rate scheme). + real(kind=dp), intent(in) :: freezing_rate + !> Whether to use the naive algorithm for time-dependent scheme. + !> (If false, use the binned tau-leaping algorithm.) + logical, intent(in) :: do_freezing_naive + + !> Call the immersion freezing subroutine according to the immersion + !> freezing scheme. + if (env_state%temp <= const%water_freeze_temp) then + if ((immersion_freezing_scheme_type == IMMERSION_FREEZING_SCHEME_ABIFM) & + .OR. (immersion_freezing_scheme_type == IMMERSION_FREEZING_SCHEME_CONST)) then + if (do_freezing_naive) then + call ice_nucleation_immersion_freezing_time_dependent_naive( & + aero_state, aero_data, env_state, del_t, & + immersion_freezing_scheme_type, freezing_rate) + else + call ice_nucleation_immersion_freezing_time_dependent( & + aero_state, aero_data, env_state, del_t, & + immersion_freezing_scheme_type, freezing_rate) + end if + else if (immersion_freezing_scheme_type == & + IMMERSION_FREEZING_SCHEME_SINGULAR) then + call ice_nucleation_immersion_freezing_singular(aero_state, & + aero_data, env_state) + else + call assert_msg(121370299, .false., & + 'Error type of immersion freezing scheme') + endif + endif + + end subroutine ice_nucleation_immersion_freezing + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Initialization for the sigular scheme, sampling the freezing temperature + !> for each particles. + subroutine ice_nucleation_singular_initialize(aero_state, aero_data, & + INAS_a, INAS_b) + implicit none + !> Aerosol state. + type(aero_state_t), intent(inout) :: aero_state + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + real(kind=dp), intent(in) :: INAS_a, INAS_b + + integer :: i_part + real(kind=dp) :: p, S, T0, temp + real(kind=dp) :: aerosol_diameter + + T0 = const%water_freeze_temp + do i_part = 1, aero_state_n_part(aero_state) + aerosol_diameter = aero_particle_dry_diameter( & + aero_state%apa%particle(i_part), aero_data) + S = const%pi * aerosol_diameter **2 + p = pmc_random() + temp = (log(1 - p) + exp(-S * exp(-INAS_a * T0 + INAS_b))) / (-S) + aero_state%apa%particle(i_part)%imf_temperature = T0 + (log(temp) & + - INAS_b) / INAS_a + end do + end subroutine ice_nucleation_singular_initialize + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Simulation for singular scheme, deciding whether to freeze for each + !> particle. Run in each time step. + subroutine ice_nucleation_immersion_freezing_singular(aero_state, aero_data,& + env_state) + + !> Aerosol state. + type(aero_state_t), intent(inout) :: aero_state + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + !> Environment state. + type(env_state_t), intent(inout) :: env_state + + real(kind=dp), allocatable :: H2O_masses(:), total_masses(:), & + H2O_frac(:) + integer :: i_part + + ! FIXME: Do this to avoid compiler warning/error, fix it in the future. + allocate(total_masses(aero_state_n_part(aero_state))) + allocate(H2O_masses(aero_state_n_part(aero_state))) + allocate(H2O_frac(aero_state_n_part(aero_state))) + + total_masses = aero_state_masses(aero_state, aero_data) + H2O_masses = aero_state_masses(aero_state, aero_data, include=["H2O"]) + H2O_frac = H2O_masses / total_masses + do i_part = 1, aero_state_n_part(aero_state) + if (aero_state%apa%particle(i_part)%frozen) then + cycle + end if + if (H2O_frac(i_part) < const%imf_water_threshold) then + cycle + end if + if (env_state%temp <= & + aero_state%apa%particle(i_part)%imf_temperature) then + aero_state%apa%particle(i_part)%frozen = .TRUE. + aero_state%apa%particle(i_part)%den_ice = & + const%reference_ice_density + aero_state%apa%particle(i_part)%ice_shape_phi = 1d0 + end if + end do + + deallocate(total_masses) + deallocate(H2O_masses) + deallocate(H2O_frac) + + end subroutine ice_nucleation_immersion_freezing_singular + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Simulation for time-dependent scheme (e.g., ABIFM, constant rate), + !> deciding whether to freeze for each particle. Run in each time step. + !> This subroutine applys the binned-tau leaping algorithm for speeding up. + subroutine ice_nucleation_immersion_freezing_time_dependent(aero_state, & + aero_data, env_state, del_t, immersion_freezing_scheme_type, & + freezing_rate) + + !> Aerosol state. + type(aero_state_t), intent(inout) :: aero_state + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + !> Environment state. + type(env_state_t), intent(inout) :: env_state + !> Total time to integrate. + real(kind=dp), intent(in) :: del_t + !> Freezing rate (only used for the constant rate scheme). + real(kind=dp), intent(in) :: freezing_rate + integer, intent(in) :: immersion_freezing_scheme_type + + integer :: i_part, i_bin, i_class, n_bins, n_class + real(kind=dp) :: a_w_ice, pis, pvs + real(kind=dp) :: p_freeze + + real(kind=dp) :: p_freeze_max, radius_max, diameter_max + + integer :: k_th, n_parts_in_bin + real(kind=dp) :: rand + real(kind=dp), allocatable :: H2O_masses(:), total_masses(:), & + H2O_frac(:) + integer :: i_spec_max + real(kind=dp) :: j_het_max + integer :: rand_geo + + + allocate(total_masses(aero_state_n_part(aero_state))) + allocate(H2O_masses(aero_state_n_part(aero_state))) + allocate(H2O_frac(aero_state_n_part(aero_state))) + + call aero_state_sort(aero_state, aero_data) + + total_masses = aero_state_masses(aero_state, aero_data) + H2O_masses = aero_state_masses(aero_state, aero_data, include=["H2O"]) + H2O_frac = H2O_masses / total_masses + pvs = env_state_saturated_vapor_pressure_wrt_water(env_state%temp) + pis = env_state_saturated_vapor_pressure_wrt_ice(env_state%temp) + a_w_ice = pis / pvs + if (immersion_freezing_scheme_type == IMMERSION_FREEZING_SCHEME_ABIFM) then + call ABIFM_max_spec(aero_data, a_w_ice, i_spec_max, j_het_max) + endif + + n_bins = aero_sorted_n_bin(aero_state%aero_sorted) + n_class = aero_sorted_n_class(aero_state%aero_sorted) + + loop_bins: do i_bin = 1, n_bins + loop_classes: do i_class = 1, n_class + n_parts_in_bin = integer_varray_n_entry(& + aero_state%aero_sorted%size_class%inverse(i_bin, i_class)) + radius_max = aero_state%aero_sorted%bin_grid%edges(i_bin + 1) + diameter_max = radius_max * 2 + if (immersion_freezing_scheme_type == & + IMMERSION_FREEZING_SCHEME_ABIFM) then + p_freeze_max = ABIFM_Pfrz_max(diameter_max, aero_data, j_het_max, & + del_t) + else if (immersion_freezing_scheme_type == & + IMMERSION_FREEZING_SCHEME_CONST) then + p_freeze_max = 1 - exp(freezing_rate * del_t) + endif + + k_th = n_parts_in_bin + 1 + loop_choosed_particles: do while(.TRUE.) + rand_geo = rand_geometric(p_freeze_max) + k_th = k_th - rand_geo + if (k_th <= 0) then + EXIT loop_choosed_particles + endif + i_part = aero_state%aero_sorted%size_class & + %inverse(i_bin, i_class)%entry(k_th) + if (aero_state%apa%particle(i_part)%frozen) then + cycle + end if + if (H2O_frac(i_part) < const%imf_water_threshold) then + cycle + end if + if (immersion_freezing_scheme_type == & + IMMERSION_FREEZING_SCHEME_ABIFM) then + p_freeze = ABIFM_Pfrz_particle( & + aero_state%apa%particle(i_part), aero_data, a_w_ice, del_t) + call warn_assert_msg(301184565, p_freeze <= p_freeze_max,& + "p_freeze > p_freeze_max.") + rand = pmc_random() + if (rand < p_freeze / p_freeze_max) then + aero_state%apa%particle(i_part)%frozen = .TRUE. + aero_state%apa%particle(i_part)%den_ice = & + const%reference_ice_density + aero_state%apa%particle(i_part)%ice_shape_phi = 1d0 + + endif + else + aero_state%apa%particle(i_part)%frozen = .TRUE. + aero_state%apa%particle(i_part)%den_ice = & + const%reference_ice_density + aero_state%apa%particle(i_part)%ice_shape_phi = 1d0 + endif + + enddo loop_choosed_particles + enddo loop_classes + enddo loop_bins + + deallocate(total_masses) + deallocate(H2O_masses) + deallocate(H2O_frac) + + end subroutine ice_nucleation_immersion_freezing_time_dependent + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Simulation for time-dependent scheme (e.g., ABIFM, constant rate), + !> deciding whether to freeze for each particle. Run in each time step. + !> This subroutine applies the naive algorithm that checks each particle. + subroutine ice_nucleation_immersion_freezing_time_dependent_naive( & + aero_state, aero_data, env_state, del_t, & + immersion_freezing_scheme_type, freezing_rate) + + !> Aerosol state. + type(aero_state_t), intent(inout) :: aero_state + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + !> Environment state. + type(env_state_t), intent(inout) :: env_state + !> Total time to integrate. + real(kind=dp), intent(in) :: del_t + !> Type of the immersion freezing scheme + integer, intent(in) :: immersion_freezing_scheme_type + !> Freezing rate (only used for the constant rate scheme). + real(kind=dp), intent(in) :: freezing_rate + + integer :: i_part + real(kind=dp) :: a_w_ice, pis, pvs + real(kind=dp) :: p_freeze = 0 + real(kind=dp), allocatable :: H2O_masses(:), total_masses(:), & + H2O_frac(:) + real(kind=dp) :: rand + + + ! FIXME: Do this to avoid compiler warning/error, fix it in the future. + allocate(total_masses(aero_state_n_part(aero_state))) + allocate(H2O_masses(aero_state_n_part(aero_state))) + allocate(H2O_frac(aero_state_n_part(aero_state))) + + total_masses = aero_state_masses(aero_state, aero_data) + H2O_masses = aero_state_masses(aero_state, aero_data, include=["H2O"]) + H2O_frac = H2O_masses / total_masses + + pvs = env_state_saturated_vapor_pressure_wrt_water(env_state%temp) + pis = env_state_saturated_vapor_pressure_wrt_ice(env_state%temp) + a_w_ice = pis / pvs + + do i_part = 1, aero_state_n_part(aero_state) + if (aero_state%apa%particle(i_part)%frozen) cycle + if (H2O_frac(i_part) < const%imf_water_threshold) cycle + rand = pmc_random() + + if (immersion_freezing_scheme_type == & + IMMERSION_FREEZING_SCHEME_ABIFM) then + p_freeze = ABIFM_Pfrz_particle(aero_state%apa%particle(i_part), & + aero_data, a_w_ice, del_t) + else if (immersion_freezing_scheme_type == & + IMMERSION_FREEZING_SCHEME_CONST) then + p_freeze = 1 - exp(freezing_rate * del_t) + end if + + if (rand < p_freeze) then + aero_state%apa%particle(i_part)%frozen = .TRUE. + aero_state%apa%particle(i_part)%den_ice = & + const%reference_ice_density + aero_state%apa%particle(i_part)%ice_shape_phi = 1d0 + end if + end do + + deallocate(total_masses) + deallocate(H2O_masses) + deallocate(H2O_frac) + + end subroutine ice_nucleation_immersion_freezing_time_dependent_naive + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Simulates melting: if the environmental temperature is above the freezing + !> temperature of water, all particles are set to be unfrozen. + subroutine ice_nucleation_melting(aero_state, aero_data, env_state) + + !> Aerosol state. + type(aero_state_t), intent(inout) :: aero_state + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + !> Environment state. + type(env_state_t), intent(inout) :: env_state + + integer :: i_part + + if (env_state%temp > const%water_freeze_temp) then + do i_part = 1, aero_state_n_part(aero_state) + aero_state%apa%particle(i_part)%frozen = .false. + aero_state%apa%particle(i_part)%den_ice = const%nan + aero_state%apa%particle(i_part)%ice_shape_phi = const%nan + end do + end if + + end subroutine ice_nucleation_melting + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Calculating the freezing probability for the particle (i_part) using ABIFM + !> method (Knopf et al.,2013) + real(kind=dp) function ABIFM_Pfrz_particle(aero_particle, aero_data, & + a_w_ice, del_t) + + !> Aerosol particle. + type(aero_particle_t), intent(in) :: aero_particle + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + !> The water activity w.r.t. ice. + real(kind=dp), intent(in) :: a_w_ice + !> Time interval. + real(kind=dp), intent(in) :: del_t + + real(kind=dp) :: aerosol_diameter + real(kind=dp) :: immersed_surface_area + real(kind=dp) :: total_vol + real(kind=dp) :: surface_ratio + real(kind=dp) :: abifm_m, abifm_c + real(kind=dp) :: j_het, j_het_x_area + integer :: i_spec + + aerosol_diameter = aero_particle_dry_diameter(aero_particle, aero_data) + immersed_surface_area = const%pi * aerosol_diameter **2 + + total_vol = 0d0 + total_vol = aero_particle_dry_volume(aero_particle, aero_data) + + j_het_x_area = 0d0 + do i_spec = 1,aero_data_n_spec(aero_data) + if (i_spec == aero_data%i_water) cycle + abifm_m = aero_data%abifm_m(i_spec) + abifm_c = aero_data%abifm_c(i_spec) + j_het = 10 ** (abifm_m * (1 - a_w_ice) + abifm_c) * 10000 + surface_ratio = aero_particle%vol(i_spec) / total_vol + j_het_x_area = j_het_x_area + j_het * immersed_surface_area * & + surface_ratio + end do + + ABIFM_Pfrz_particle = 1 - exp(-j_het_x_area * del_t) + + end function ABIFM_Pfrz_particle + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Finding the maximum heterogeneous ice nucleation rate coefficient. + subroutine ABIFM_max_spec(aero_data, a_w_ice, i_spec_max, j_het_max) + + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + !> The water activity w.r.t. ice. + real(kind=dp), intent(in) :: a_w_ice + !> The index of the maximum J_het species. + integer, intent(out) :: i_spec_max + !> The maximum value of J_het among all species. + real(kind=dp), intent(out) :: j_het_max + + real(kind=dp) :: abifm_m, abifm_c + real(kind=dp) :: j_het + integer :: i_spec + + j_het_max = const%nan + do i_spec = 1,aero_data_n_spec(aero_data) + if (i_spec == aero_data%i_water) cycle + abifm_m = aero_data%abifm_m(i_spec) + abifm_c = aero_data%abifm_c(i_spec) + j_het = 10 ** (abifm_m * (1 - a_w_ice) + abifm_c) * 10000 + if (j_het > j_het_max) then + j_het_max = j_het + i_spec_max = i_spec + end if + end do + + end subroutine ABIFM_max_spec + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Calculating the maximum freezing probability for particles in + !> one bin using ABIFM method (Knopf et al.,2013). Only used by + !> the binned-tau leaping algorithm. + real(kind=dp) function ABIFM_Pfrz_max(diameter_max, aero_data, j_het_max, & + del_t) + + !> Aerosol data. + type(aero_data_t), intent(in) :: aero_data + !> Maximum diameter. + real(kind=dp), intent(in) :: diameter_max + !> Time interval. + real(kind=dp), intent(in) :: del_t + !> Maximum J_het among all species. + real(kind=dp), intent(in) :: j_het_max + + real(kind=dp) :: immersed_surface_area + + immersed_surface_area = const%pi * diameter_max **2 + ABIFM_Pfrz_max = 1 - exp(-j_het_max * immersed_surface_area * del_t) + + end function ABIFM_Pfrz_max + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Read the specification for immersion freezing scheme. + subroutine spec_file_read_immersion_freezing_scheme_type(file, & + immersion_freezing_scheme_type) + + !> Spec file. + type(spec_file_t), intent(inout) :: file + !> Kernel type. + integer, intent(out) :: immersion_freezing_scheme_type + + character(len=SPEC_LINE_MAX_VAR_LEN) :: imf_scheme + + !> \page input_format_immersion_freezing_scheme Input File Format: Immersion freezing scheme + !! + !! The immersion freezing scheme is specified by the parameter: + !! - \b immersion_freezing_scheme (string): the type of + !! immersion freezing scheme must be one of: + !! \c const for freezing with a constant freezing rate; + !! \c singular for the INAS scheme based on singular perspective; + !! \c ABIFM for the ABIFM scheme based on CNT perspetive. + !! + !! If \c immersion_freezing_scheme is \c singular, the parameters INAS_a + !! and INAS_b need to be provided. + !! If \c immersion_freezing_scheme is \c const, the parameter freezing_rate + !! need to be provided. + !! If \c immersion_freezing_scheme is \c ABIFM or const, a logical variable + !! \c do_freezing_naive need to be provided. + !! + !! See also: + !! - \ref spec_file_format --- the input file text format + + call spec_file_read_string(file, 'immersion_freezing_scheme', imf_scheme) + if (trim(imf_scheme) == 'const') then + immersion_freezing_scheme_type = IMMERSION_FREEZING_SCHEME_CONST + elseif (trim(imf_scheme) == 'singular') then + immersion_freezing_scheme_type = IMMERSION_FREEZING_SCHEME_SINGULAR + elseif (trim(imf_scheme) == 'ABIFM') then + immersion_freezing_scheme_type = IMMERSION_FREEZING_SCHEME_ABIFM + else + call spec_file_die_msg(920761229, file, & + "Unknown immersion freezing scheme: " // trim(imf_scheme)) + end if + + end subroutine spec_file_read_immersion_freezing_scheme_type + +end module pmc_ice_nucleation diff --git a/src/partmc.F90 b/src/partmc.F90 index b873970f2..f7e34947e 100644 --- a/src/partmc.F90 +++ b/src/partmc.F90 @@ -417,6 +417,8 @@ subroutine partmc_part(file) !! nucleation. If \c do_nucleation is \c yes, then the following !! parameters must also be provided: !! - \subpage input_format_nucleate + !! - \b do_freezing (logical): whether to perform particle + !! freezing. !! - \b rand_init (integer): if greater than zero then use as !! the seed for the random number generator, or if zero then !! generate a random seed for the random number generator --- diff --git a/src/rand.F90 b/src/rand.F90 index 6e390da22..7a596bea1 100644 --- a/src/rand.F90 +++ b/src/rand.F90 @@ -697,4 +697,75 @@ end subroutine uuid4_str !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Generate a Geometric-distributed random number with the given + !> probability. + !! + !! See https://en.wikipedia.org/wiki/Geometric_distribution and + !! https://www.ucl.ac.uk/~ucakarc/work/software/randgen.f for + !! more details + + integer function rand_geometric(p) + implicit none + !> Sample probability. + real(kind=dp), intent(in) :: p + +#ifdef PMC_USE_GSL + real(kind=c_double) :: p_c + integer(kind=c_int), target :: harvest + type(c_ptr) :: harvest_ptr + logical :: acceptable +#else + real(kind=dp) :: U, TINY +#endif + +#ifdef PMC_USE_GSL + interface + integer(kind=c_int) function pmc_rand_geometric_gsl(p, harvest) bind(c) + use iso_c_binding + real(kind=c_double), value :: p + type(c_ptr), value :: harvest + end function pmc_rand_geometric_gsl + end interface +#endif + call assert(386975305, p >= 0d0) + call assert(667756333, p <= 1d0) +#ifdef PMC_USE_GSL + p_c = real(p, kind=c_double) + harvest_ptr = c_loc(harvest) + acceptable = .false. + do while(.not. acceptable) + call rand_check_gsl(218328792, & + pmc_rand_geometric_gsl(p_c, harvest_ptr)) + rand_geometric = int(harvest) + if (rand_geometric .ge. 1) then + acceptable = .true. + end if + end do +#else + TINY = 1.0D-12 + rand_geometric = 0 + + call assert_msg(927129543, (P.ge.0.0D0).and.(P.le.1.0D0),& + "Range error") + + if (P.gt.0.9D0) then + rand_geometric = rand_geometric + 1 + U = pmc_random() + do while( U.gt.P ) + rand_geometric = rand_geometric + 1 + U = pmc_random() + enddo + else + U = pmc_random() + + if (P.gt.TINY) then + rand_geometric = 1 + INT( DLOG(U)/DLOG(1.0D0-P) ) + else + rand_geometric = 1 + INT(-DLOG(U)/P) + endif + endif +#endif + end function rand_geometric +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + end module pmc_rand diff --git a/src/rand_gsl.c b/src/rand_gsl.c index db200abbf..179b927b6 100644 --- a/src/rand_gsl.c +++ b/src/rand_gsl.c @@ -150,3 +150,18 @@ int pmc_rand_binomial_gsl(int n, double p, int *harvest) *harvest = gsl_ran_binomial(pmc_rand_gsl_rng, p, u); return PMC_RAND_GSL_SUCCESS; } + +/** \brief Generate a Geometric-distributed random integer. + * + * \param p The sample probability for the distribution. + * \param harvest A pointer to the generated random number. + * \return PMC_RAND_GSL_SUCCESS on success, otherwise an error code. + */ +int pmc_rand_geometric_gsl(double p, int *harvest) +{ + if (!pmc_rand_gsl_rng) { + return PMC_RAND_GSL_NOT_INIT; + } + *harvest = gsl_ran_geometric(pmc_rand_gsl_rng, p); + return PMC_RAND_GSL_SUCCESS; +} diff --git a/src/run_part.F90 b/src/run_part.F90 index 2c380e713..76a2ddcfe 100644 --- a/src/run_part.F90 +++ b/src/run_part.F90 @@ -22,6 +22,7 @@ module pmc_run_part use pmc_coagulation_dist use pmc_coag_kernel use pmc_nucleate + use pmc_ice_nucleation use pmc_mpi use pmc_camp_interface use pmc_photolysis @@ -69,6 +70,18 @@ module pmc_run_part logical :: do_coagulation !> Whether to do nucleation. logical :: do_nucleation + !> Whether to do freezing. + logical :: do_immersion_freezing + !> The immersion freezing scheme options. + integer :: immersion_freezing_scheme_type + !> The INAS parameters for "singular" scheme. + real(kind=dp) :: INAS_a + real(kind=dp) :: INAS_b + !> The freezing rate parameter for "const" scheme. + real(kind=dp) :: freezing_rate + !> Whether to use the naive algorithm for time-dependent scheme. + !> (If false, use the binned tau-leaping algorithm.) + logical :: do_freezing_naive !> Allow doubling if needed. logical :: allow_doubling !> Allow halving if needed. @@ -232,10 +245,16 @@ subroutine run_part(scenario, env_state, aero_data, aero_state, gas_data, & call print_part_progress(run_part_opt%i_repeat, time, & global_n_part, 0, 0, 0, 0, 0, t_wall_elapsed, t_wall_remain) end if + ! initialize the immersion freezing temperature for Singular scheme + if (run_part_opt%do_immersion_freezing .and. & + (run_part_opt%immersion_freezing_scheme_type .eq. & + IMMERSION_FREEZING_SCHEME_SINGULAR)) then + call ice_nucleation_singular_initialize(aero_state, aero_data, & + run_part_opt%INAS_a, run_part_opt%INAS_b) + end if i_cur = 1 i_next = n_time - call run_part_timeblock(scenario, env_state, aero_data, aero_state, & gas_data, gas_state, run_part_opt, & #ifdef PMC_USE_CAMP @@ -316,6 +335,12 @@ integer function pmc_mpi_pack_size_run_part_opt(val) + pmc_mpi_pack_size_integer(val%nucleate_weight_class) & + pmc_mpi_pack_size_logical(val%do_coagulation) & + pmc_mpi_pack_size_logical(val%do_nucleation) & + + pmc_mpi_pack_size_logical(val%do_immersion_freezing) & + + pmc_mpi_pack_size_integer(val%immersion_freezing_scheme_type) & + + pmc_mpi_pack_size_real(val%INAS_a) & + + pmc_mpi_pack_size_real(val%INAS_b) & + + pmc_mpi_pack_size_real(val%freezing_rate) & + + pmc_mpi_pack_size_logical(val%do_freezing_naive) & + pmc_mpi_pack_size_logical(val%allow_doubling) & + pmc_mpi_pack_size_logical(val%allow_halving) & + pmc_mpi_pack_size_logical(val%do_condensation) & @@ -367,6 +392,14 @@ subroutine pmc_mpi_pack_run_part_opt(buffer, position, val) call pmc_mpi_pack_integer(buffer, position, val%nucleate_weight_class) call pmc_mpi_pack_logical(buffer, position, val%do_coagulation) call pmc_mpi_pack_logical(buffer, position, val%do_nucleation) + call pmc_mpi_pack_logical(buffer, position, val%do_immersion_freezing) + call pmc_mpi_pack_integer(buffer, position, & + val%immersion_freezing_scheme_type) + + call pmc_mpi_pack_real(buffer, position, val%INAS_a) + call pmc_mpi_pack_real(buffer, position, val%INAS_b) + call pmc_mpi_pack_real(buffer, position, val%freezing_rate) + call pmc_mpi_pack_logical(buffer, position, val%do_freezing_naive) call pmc_mpi_pack_logical(buffer, position, val%allow_doubling) call pmc_mpi_pack_logical(buffer, position, val%allow_halving) call pmc_mpi_pack_logical(buffer, position, val%do_condensation) @@ -421,6 +454,12 @@ subroutine pmc_mpi_unpack_run_part_opt(buffer, position, val) call pmc_mpi_unpack_integer(buffer, position, val%nucleate_weight_class) call pmc_mpi_unpack_logical(buffer, position, val%do_coagulation) call pmc_mpi_unpack_logical(buffer, position, val%do_nucleation) + call pmc_mpi_unpack_logical(buffer, position, val%do_immersion_freezing) + call pmc_mpi_unpack_integer(buffer, position, val%immersion_freezing_scheme_type) + call pmc_mpi_unpack_real(buffer, position, val%INAS_a) + call pmc_mpi_unpack_real(buffer, position, val%INAS_b) + call pmc_mpi_unpack_real(buffer, position, val%freezing_rate) + call pmc_mpi_unpack_logical(buffer, position, val%do_freezing_naive) call pmc_mpi_unpack_logical(buffer, position, val%allow_doubling) call pmc_mpi_unpack_logical(buffer, position, val%allow_halving) call pmc_mpi_unpack_logical(buffer, position, val%do_condensation) @@ -673,6 +712,36 @@ subroutine spec_file_read_run_part(file, run_part_opt, aero_data, & else run_part_opt%nucleate_type = NUCLEATE_TYPE_INVALID end if + call spec_file_read_logical(file, 'do_immersion_freezing', & + run_part_opt%do_immersion_freezing) + + if (run_part_opt%do_immersion_freezing) then + + call spec_file_read_immersion_freezing_scheme_type(file, & + run_part_opt%immersion_freezing_scheme_type) + + if (run_part_opt%immersion_freezing_scheme_type .eq. & + IMMERSION_FREEZING_SCHEME_CONST) then + call spec_file_read_real(file, 'freezing_rate', & + run_part_opt%freezing_rate) + endif + + if ((run_part_opt%immersion_freezing_scheme_type .eq. & + IMMERSION_FREEZING_SCHEME_ABIFM) .or. & + (run_part_opt%immersion_freezing_scheme_type .eq. & + IMMERSION_FREEZING_SCHEME_CONST)) then + call spec_file_read_logical(file, 'do_freezing_naive', & + run_part_opt%do_freezing_naive) + endif + + + if (run_part_opt%immersion_freezing_scheme_type .eq.& + IMMERSION_FREEZING_SCHEME_SINGULAR) then + call spec_file_read_real(file, 'INAS_a', run_part_opt%INAS_a) + call spec_file_read_real(file, 'INAS_b', run_part_opt%INAS_b) + endif + + endif call spec_file_read_integer(file, 'rand_init', rand_init) call spec_file_read_logical(file, 'allow_doubling', & @@ -843,7 +912,12 @@ subroutine run_part_timestep(scenario, env_state, aero_data, aero_state, & - n_part_before progress_n_nuc = progress_n_nuc + n_nuc end if - + if (run_part_opt%do_immersion_freezing) then + call ice_nucleation_immersion_freezing(aero_state, aero_data, env_state, & + run_part_opt%del_t, run_part_opt%immersion_freezing_scheme_type, & + run_part_opt%freezing_rate, run_part_opt%do_freezing_naive) + call ice_nucleation_melting(aero_state, aero_data, env_state) + end if if (run_part_opt%do_coagulation) then if (run_part_opt%parallel_coag_type & == PARALLEL_COAG_TYPE_LOCAL) then diff --git a/test/additive/aero_data.dat b/test/additive/aero_data.dat index 41cdc37fc..b38ac540e 100644 --- a/test/additive/aero_data.dat +++ b/test/additive/aero_data.dat @@ -1,2 +1,2 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) +H2O 1000 0 18d-3 0 0 0 diff --git a/test/additive/run_part.spec b/test/additive/run_part.spec index b2da7e0ab..0b28eb747 100644 --- a/test/additive/run_part.spec +++ b/test/additive/run_part.spec @@ -43,6 +43,7 @@ additive_kernel_coeff 1000 # additive kernel constant do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/average/aero_data.dat b/test/average/aero_data.dat index adc790ba7..65d773f46 100644 --- a/test/average/aero_data.dat +++ b/test/average/aero_data.dat @@ -1,21 +1,21 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa -H2O 1000 0 18d-3 0 -SO4 1800 0 96d-3 0.53 -NO3 1800 0 62d-3 0.53 -Cl 2200 0 35.5d-3 1.28 -NH4 1800 0 18d-3 0.53 -CO3 2600 0 60d-3 0.53 -MSA 1800 0 95d-3 0.53 -Na 2200 0 23d-3 1.28 -Ca 2600 0 40d-3 0.53 -OC 1000 0 1d-3 0.1 -BC 1800 0 1d-3 0 -OIN 2600 0 1d-3 0.1 -ARO1 1000 0 150d-3 0.1 -ARO2 1000 0 150d-3 0.1 -ALK1 1000 0 140d-3 0.1 -OLE1 1000 0 140d-3 0.1 -API1 1000 0 184d-3 0.1 -API2 1000 0 184d-3 0.1 -LIM1 1000 0 200d-3 0.1 -LIM2 1000 0 200d-3 0.1 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa abifm_m (1) abifm_c (1) +H2O 1000 0 18d-3 0 0 0 +SO4 1800 0 96d-3 0.53 0 0 +NO3 1800 0 62d-3 0.53 0 0 +Cl 2200 0 35.5d-3 1.28 0 0 +NH4 1800 0 18d-3 0.53 0 0 +CO3 2600 0 60d-3 0.53 0 0 +MSA 1800 0 95d-3 0.53 0 0 +Na 2200 0 23d-3 1.28 0 0 +Ca 2600 0 40d-3 0.53 0 0 +OC 1000 0 1d-3 0.1 0 0 +BC 1800 0 1d-3 0 0 0 +OIN 2600 0 1d-3 0.1 0 0 +ARO1 1000 0 150d-3 0.1 0 0 +ARO2 1000 0 150d-3 0.1 0 0 +ALK1 1000 0 140d-3 0.1 0 0 +OLE1 1000 0 140d-3 0.1 0 0 +API1 1000 0 184d-3 0.1 0 0 +API2 1000 0 184d-3 0.1 0 0 +LIM1 1000 0 200d-3 0.1 0 0 +LIM2 1000 0 200d-3 0.1 0 0 diff --git a/test/average/run_part.spec b/test/average/run_part.spec index 22a7fb06e..1dad42d99 100644 --- a/test/average/run_part.spec +++ b/test/average/run_part.spec @@ -40,6 +40,7 @@ do_coagulation no # whether to do coagulation (yes/no) do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling no # whether to allow doubling (yes/no) diff --git a/test/bidisperse/aero_data.dat b/test/bidisperse/aero_data.dat index 41cdc37fc..cc430aaee 100644 --- a/test/bidisperse/aero_data.dat +++ b/test/bidisperse/aero_data.dat @@ -1,2 +1,2 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) +H2O 1000 0 18d-3 0 0 0 diff --git a/test/bidisperse/run_part.spec b/test/bidisperse/run_part.spec index 260689cd6..95dbb750d 100644 --- a/test/bidisperse/run_part.spec +++ b/test/bidisperse/run_part.spec @@ -41,6 +41,7 @@ coag_kernel sedi # coagulation kernel do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/brownian/aero_data.dat b/test/brownian/aero_data.dat index 41cdc37fc..cc430aaee 100644 --- a/test/brownian/aero_data.dat +++ b/test/brownian/aero_data.dat @@ -1,2 +1,2 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) +H2O 1000 0 18d-3 0 0 0 diff --git a/test/brownian/run_part.spec b/test/brownian/run_part.spec index a2ba909fe..164bc7c39 100644 --- a/test/brownian/run_part.spec +++ b/test/brownian/run_part.spec @@ -41,6 +41,7 @@ coag_kernel brown # coagulation kernel do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/camp/camp.spec b/test/camp/camp.spec index f563bbbf9..30e6444f2 100644 --- a/test/camp/camp.spec +++ b/test/camp/camp.spec @@ -39,6 +39,8 @@ do_coagulation no # whether to do coagulation (yes/no) do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) + rand_init 0 # random initialization (0 to use time) allow_doubling no # whether to allow doubling (yes/no) diff --git a/test/camp/monarch_mod37/tsigaridis_2_product_SOA_scheme/species.json b/test/camp/monarch_mod37/tsigaridis_2_product_SOA_scheme/species.json index e46bc5511..e0c09ac1c 100644 --- a/test/camp/monarch_mod37/tsigaridis_2_product_SOA_scheme/species.json +++ b/test/camp/monarch_mod37/tsigaridis_2_product_SOA_scheme/species.json @@ -9,6 +9,8 @@ "type" : "CHEM_SPEC", "absolute integration tolerance" : 1.0E-03, "molecular weight [kg mol-1]" : 0.08806, + "abifm_m": 0, + "abifm_c": 0, "description" : "gas-phase product from isoprene oxidation by OH", "notes" : [ "using diffusion coefficient from dry deposition", @@ -20,6 +22,8 @@ "type" : "CHEM_SPEC", "absolute integration tolerance" : 1.0E-03, "molecular weight [kg mol-1]" : 0.09003, + "abifm_m": 0, + "abifm_c": 0, "description" : "gas-phase product from isoprene oxidation by O3", "notes" : [ "using diffusion coefficient from dry deposition", @@ -31,6 +35,8 @@ "type" : "CHEM_SPEC", "absolute integration tolerance" : 1.0E-03, "molecular weight [kg mol-1]" : 0.17025, + "abifm_m": 0, + "abifm_c": 0, "description" : "gas-phase product from monoterpene oxidation by OH", "notes" : [ "using diffusion coefficient from dry deposition", @@ -42,6 +48,8 @@ "type" : "CHEM_SPEC", "absolute integration tolerance" : 1.0E-03, "molecular weight [kg mol-1]" : 0.202162, + "abifm_m": 0, + "abifm_c": 0, "description" : "gas-phase product from monoterpene oxidation by O3", "notes" : [ "using diffusion coefficient from dry deposition", @@ -55,6 +63,8 @@ "absolute integration tolerance" : 1.0E-05, "density [kg m-3]" : 1800.0, "molecular weight [kg mol-1]" : 0.41475, + "abifm_m": 0, + "abifm_c": 0, "description" : "lumped hydrophobic particulate matter", "num_ions" : 0, "kappa" : 0.0, @@ -67,6 +77,8 @@ "absolute integration tolerance" : 1.0E-05, "density [kg m-3]" : 1800.0, "molecular weight [kg mol-1]" : 0.08806, + "abifm_m": 0, + "abifm_c": 0, "description" : "First isoprene SOA species in 2-product scheme", "num_ions" : 0, "kappa" : 0.1, @@ -82,6 +94,8 @@ "absolute integration tolerance" : 1.0E-05, "density [kg m-3]" : 1800.0, "molecular weight [kg mol-1]" : 0.09003, + "abifm_m": 0, + "abifm_c": 0, "description" : "Second isoprene SOA species in 2-product scheme", "num_ions" : 0, "kappa" : 0.1, @@ -97,6 +111,8 @@ "absolute integration tolerance" : 1.0E-05, "density [kg m-3]" : 1800.0, "molecular weight [kg mol-1]" : 0.17025, + "abifm_m": 0, + "abifm_c": 0, "description" : "First monoterpene SOA species in 2-product scheme", "num_ions" : 0, "kappa" : 0.1, @@ -113,6 +129,8 @@ "absolute integration tolerance" : 1.0E-05, "density [kg m-3]" : 1800.0, "molecular weight [kg mol-1]" : 0.202162, + "abifm_m": 0, + "abifm_c": 0, "description" : "Second monoterpene SOA species in 2-product scheme", "num_ions" : 0, "kappa" : 0.1, diff --git a/test/condense/aero_data.dat b/test/condense/aero_data.dat index af22101b2..9d466783c 100644 --- a/test/condense/aero_data.dat +++ b/test/condense/aero_data.dat @@ -1,21 +1,21 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -SO4 1800 1 96d-3 0 -NO3 1800 1 62d-3 0 -Cl 2200 1 35.5d-3 0 -NH4 1800 1 18d-3 0 -MSA 1800 0 95d-3 0.53 -ARO1 1400 0 150d-3 0.1 -ARO2 1400 0 150d-3 0.1 -ALK1 1400 0 140d-3 0.1 -OLE1 1400 0 140d-3 0.1 -API1 1400 0 184d-3 0.1 -API2 1400 0 184d-3 0.1 -LIM1 1400 0 200d-3 0.1 -LIM2 1400 0 200d-3 0.1 -CO3 2600 1 60d-3 0 -Na 2200 1 23d-3 0 -Ca 2600 1 40d-3 0 -OIN 2600 0 1d-3 0.1 -OC 1400 0 1d-3 0.1 -BC 1800 0 1d-3 0 -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) +SO4 1800 1 96d-3 0 0 0 +NO3 1800 1 62d-3 0 0 0 +Cl 2200 1 35.5d-3 0 0 0 +NH4 1800 1 18d-3 0 0 0 +MSA 1800 0 95d-3 0.53 0 0 +ARO1 1400 0 150d-3 0.1 0 0 +ARO2 1400 0 150d-3 0.1 0 0 +ALK1 1400 0 140d-3 0.1 0 0 +OLE1 1400 0 140d-3 0.1 0 0 +API1 1400 0 184d-3 0.1 0 0 +API2 1400 0 184d-3 0.1 0 0 +LIM1 1400 0 200d-3 0.1 0 0 +LIM2 1400 0 200d-3 0.1 0 0 +CO3 2600 1 60d-3 0 0 0 +Na 2200 1 23d-3 0 0 0 +Ca 2600 1 40d-3 0 0 0 +OIN 2600 0 1d-3 0.1 0 0 +OC 1400 0 1d-3 0.1 0 0 +BC 1800 0 1d-3 0 0 0 +H2O 1000 0 18d-3 0 0 0 diff --git a/test/condense/run_part.spec b/test/condense/run_part.spec index 23cf93f4b..d0e0bd20c 100644 --- a/test/condense/run_part.spec +++ b/test/condense/run_part.spec @@ -41,6 +41,7 @@ do_condensation yes # whether to do condensation (yes/no) do_init_equilibrate yes # whether to initially equilibrate water (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to use time) allow_doubling no # whether to allow doubling (yes/no) diff --git a/test/emission/aero_data.dat b/test/emission/aero_data.dat index 83d7a7091..d52f270c3 100644 --- a/test/emission/aero_data.dat +++ b/test/emission/aero_data.dat @@ -1,4 +1,4 @@ # dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -SI 10000 0 18d-3 0 -SB 100 0 18d-3 0 -SE 1500 0 18d-3 0 +SI 10000 0 18d-3 0 0 0 +SB 100 0 18d-3 0 0 0 +SE 1500 0 18d-3 0 0 0 diff --git a/test/emission/run_part.spec b/test/emission/run_part.spec index 016496f3a..5947e9ed2 100644 --- a/test/emission/run_part.spec +++ b/test/emission/run_part.spec @@ -40,6 +40,7 @@ do_coagulation no # whether to do coagulation (yes/no) do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/fractal/aero_data.dat b/test/fractal/aero_data.dat index 21e82154f..9d5459403 100644 --- a/test/fractal/aero_data.dat +++ b/test/fractal/aero_data.dat @@ -1,2 +1,2 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -BC 4200 0 1d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) +BC 4200 0 1d-3 0 0 0 diff --git a/test/fractal/run_part_brown_cont_df_1_8_restart.spec b/test/fractal/run_part_brown_cont_df_1_8_restart.spec index 173d16530..176ef3cb2 100644 --- a/test/fractal/run_part_brown_cont_df_1_8_restart.spec +++ b/test/fractal/run_part_brown_cont_df_1_8_restart.spec @@ -34,6 +34,7 @@ coag_kernel brown_cont # coagulation kernel do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/fractal/run_part_brown_cont_df_1_8_upto1000s.spec b/test/fractal/run_part_brown_cont_df_1_8_upto1000s.spec index 6df294249..ba06c265c 100644 --- a/test/fractal/run_part_brown_cont_df_1_8_upto1000s.spec +++ b/test/fractal/run_part_brown_cont_df_1_8_upto1000s.spec @@ -44,6 +44,7 @@ coag_kernel brown_cont # coagulation kernel do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/fractal/run_part_brown_free_df_2_4_restart.spec b/test/fractal/run_part_brown_free_df_2_4_restart.spec index fb5356fbd..375381264 100644 --- a/test/fractal/run_part_brown_free_df_2_4_restart.spec +++ b/test/fractal/run_part_brown_free_df_2_4_restart.spec @@ -34,6 +34,7 @@ coag_kernel brown_free # coagulation kernel do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/fractal/run_part_brown_free_df_2_4_upto1000s.spec b/test/fractal/run_part_brown_free_df_2_4_upto1000s.spec index 1051a29e3..a209b132e 100644 --- a/test/fractal/run_part_brown_free_df_2_4_upto1000s.spec +++ b/test/fractal/run_part_brown_free_df_2_4_upto1000s.spec @@ -44,6 +44,7 @@ coag_kernel brown_free # coagulation kernel do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/freezing/README b/test/freezing/README new file mode 100644 index 000000000..5e3cc4ebc --- /dev/null +++ b/test/freezing/README @@ -0,0 +1,13 @@ + +Immersion freezing Test-case +================================== +The initial condition is an internally mixed INPs in a constant +temperature process. All particles initially immersed in water, +with wet diameter of ~2.15 um in average, and dry diameter +of 1 um in average. The log-std of wet and dry diameter is 0.5. +All particles contain an insoluble INP, with 50% of surface +covered by illite, and 50% of surface by Fe2O3. The temperature +is -20 degrees celsius. + +The theoretical frozen fraction results are calculated using +the equation (3.40) in Tang 2024 (https://www.ideals.illinois.edu/items/131711). diff --git a/test/freezing/aero_back.dat b/test/freezing/aero_back.dat new file mode 100644 index 000000000..3879b8f6b --- /dev/null +++ b/test/freezing/aero_back.dat @@ -0,0 +1,6 @@ +# time (s) +# rate (s^{-1}) +# aerosol distribution filename +time 0 +rate 0 +dist aero_back_dist.dat diff --git a/test/freezing/aero_back_dist.dat b/test/freezing/aero_back_dist.dat new file mode 100644 index 000000000..c56ed75ed --- /dev/null +++ b/test/freezing/aero_back_dist.dat @@ -0,0 +1 @@ +# empty aerosol state with no particles diff --git a/test/freezing/aero_data.dat b/test/freezing/aero_data.dat new file mode 100644 index 000000000..3735cb147 --- /dev/null +++ b/test/freezing/aero_data.dat @@ -0,0 +1,30 @@ +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) +H2O 1000 0 18d-3 0 0 0 +#OIN 2600 0 1d-3 0.1 28.13797 -2.92414 + +ILT 2750 0 389.34d-3 0.003 54.48075 -10.66873 +# Illite: +# dens, molec wght: http://www.webmineral.com/data/Illite.shtml#.ZA0JpezMJ_Q +# kappa: 0.1039/b901585j + +#Al2O3 3970 0 101.96d-3 0.01 14.96639 1.60671 +# Aluminium oxide: +# dens, molec wght: https://en.wikipedia.org/wiki/Aluminium_oxide +# kappa: [unknown] + +Fe2O3 5240 0 159.69d-3 0.01 17.62106 1.42411 +# Iron(III) oxide: +# den, molec wght: https://en.wikipedia.org/wiki/Iron(III)_oxide +# kappa: [unknown] + +#KLN 2650 0 258.16d-3 0.023 54.58834 -10.54758 +# Kaolinite +# den: https://cameochemicals.noaa.gov/chemical/25036 +# molec wght: http://webmineral.com/data/Kaolinite.shtml#.ZA0TTezMJ_Q +# kappa: 10.1039/b901585j +# + +#NVF 5240 0 1d-3 0.01 0.0 0.0 +# NeVer Freeze (other parameters are same as OIN) +#VOLC 2600 0 1d-3 0.01 28.13797 -2.92414 +# Volcano (Steinke et al., 2011) diff --git a/test/freezing/aero_emit.dat b/test/freezing/aero_emit.dat new file mode 100644 index 000000000..34fb5fb66 --- /dev/null +++ b/test/freezing/aero_emit.dat @@ -0,0 +1,6 @@ +# time (s) +# rate (s^{-1}) +# aerosol distribution filename +time 0 +rate 0 +dist aero_emit_dist.dat diff --git a/test/freezing/aero_emit_dist.dat b/test/freezing/aero_emit_dist.dat new file mode 100644 index 000000000..c56ed75ed --- /dev/null +++ b/test/freezing/aero_emit_dist.dat @@ -0,0 +1 @@ +# empty aerosol state with no particles diff --git a/test/freezing/aero_init_comp.dat b/test/freezing/aero_init_comp.dat new file mode 100644 index 000000000..014decb58 --- /dev/null +++ b/test/freezing/aero_init_comp.dat @@ -0,0 +1,5 @@ +# mass fractions +H2O 0.69257407 +Fe2O3 0.20161601 +ILT 0.10580993 + diff --git a/test/freezing/aero_init_dist.dat b/test/freezing/aero_init_dist.dat new file mode 100644 index 000000000..8b4c874dd --- /dev/null +++ b/test/freezing/aero_init_dist.dat @@ -0,0 +1,7 @@ +# +mode_name init1 # name of mode source +mass_frac aero_init_comp.dat # species mass fractions +diam_type geometric # type of diameter specified +mode_type mono # type of distribution +num_conc 1e8 # particle number concentration (#/m^3) +diam 2.154434690031885e-06 # geometric mean diameter (m) diff --git a/test/freezing/gas_back.dat b/test/freezing/gas_back.dat new file mode 100644 index 000000000..38737461f --- /dev/null +++ b/test/freezing/gas_back.dat @@ -0,0 +1,7 @@ +# time (s) +# rate (s^{-1}) +# mixing ratios (ppb) +time 0 +rate 1 +# these lines intentionally left empty +# to serve as an empty gas state diff --git a/test/freezing/gas_data.dat b/test/freezing/gas_data.dat new file mode 100644 index 000000000..5bf86180b --- /dev/null +++ b/test/freezing/gas_data.dat @@ -0,0 +1,7 @@ +# list of gas species +H2SO4 +HNO3 +HCl +NH3 +NO +NO2 diff --git a/test/freezing/gas_emit.dat b/test/freezing/gas_emit.dat new file mode 100644 index 000000000..80db9bcf2 --- /dev/null +++ b/test/freezing/gas_emit.dat @@ -0,0 +1,7 @@ +# time (s) +# rate (s^{-1}) +# emissions (mol m^{-2} s^{-1}) +time 0 +rate 1 +# these lines intentionally left empty +# to serve as an empty gas state diff --git a/test/freezing/gas_init.dat b/test/freezing/gas_init.dat new file mode 100644 index 000000000..5a02fb933 --- /dev/null +++ b/test/freezing/gas_init.dat @@ -0,0 +1,3 @@ +# species init mixing ratio (ppb) +# these lines intentionally left empty +# to serve as an empty gas state diff --git a/test/freezing/height.dat b/test/freezing/height.dat new file mode 100644 index 000000000..d4927c431 --- /dev/null +++ b/test/freezing/height.dat @@ -0,0 +1,4 @@ +# time (s) +# height (m) +time 0 +height 1000 diff --git a/test/freezing/pressure.dat b/test/freezing/pressure.dat new file mode 100644 index 000000000..20df8e87a --- /dev/null +++ b/test/freezing/pressure.dat @@ -0,0 +1,4 @@ +# time (s) +# pressure (Pa) +time 0 +pressure 50000.000 diff --git a/test/freezing/run_part.ABIFM.naive.spec b/test/freezing/run_part.ABIFM.naive.spec new file mode 100644 index 000000000..c2a5876f0 --- /dev/null +++ b/test/freezing/run_part.ABIFM.naive.spec @@ -0,0 +1,53 @@ +run_type particle # particle-resolved run +output_prefix out/ABIFM_naive/freezing_part # prefix of output files +n_repeat 5 # number of Monte Carlo repeats +n_part 1000 # total number of particles +restart no # whether to restart from saved state (yes/no) +do_select_weighting yes # whether to select weighting explicitly (yes/no) +weight_type flat + + +t_max 600 # total simulation time (s) +del_t 1 # timestep (s) +t_output 60 # output interval (0 disables) (s) +t_progress 120 # progress printing interval (0 disables) (s) + +do_camp_chem no # whether to run the campible chemistry module +do_tchem no # whether to use TChem for chemistry + +gas_data gas_data.dat # file containing gas data +gas_init gas_init.dat # initial gas mixing ratios + +aerosol_data aero_data.dat # file containing aerosol data +do_fractal no # whether to do fractal treatment +aerosol_init aero_init_dist.dat # aerosol initial condition file + +temp_profile temp.dat # temperature profile file +pressure_profile pressure.dat # pressure profile file +height_profile height.dat # height profile file +gas_emissions gas_emit.dat # gas emissions file +gas_background gas_back.dat # background gas mixing ratios file +aero_emissions aero_emit.dat # aerosol emissions file +aero_background aero_back.dat # aerosol background file +loss_function none # particle loss function + +rel_humidity 1.000 # initial relative humidity (1) +latitude 40 # latitude (degrees, -90 to 90) +longitude 0 # longitude (degrees, -180 to 180) +altitude 0 # altitude (m) +start_time 0 # start time (s since 00:00 UTC) +start_day 1 # start day of year (UTC) + +do_coagulation no # whether to do coagulation (yes/no) +do_condensation no # whether to do condensation (yes/no) +do_mosaic no # whether to do MOSAIC (yes/no) +do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing yes # whether to do freezing (yes/no) +immersion_freezing_scheme ABIFM +do_freezing_naive yes # whether to use naive algorithm (yes/no) + +rand_init 0 # random initialization (0 to auto-generate) +allow_doubling yes # whether to allow doubling (yes/no) +allow_halving yes # whether to allow halving (yes/no) +record_removals no # whether to record particle removals (yes/no) +do_parallel no # whether to run in parallel (yes/no) diff --git a/test/freezing/run_part.ABIFM.spec b/test/freezing/run_part.ABIFM.spec new file mode 100644 index 000000000..6b025a404 --- /dev/null +++ b/test/freezing/run_part.ABIFM.spec @@ -0,0 +1,53 @@ +run_type particle # particle-resolved run +output_prefix out/ABIFM/freezing_part # prefix of output files +n_repeat 5 # number of Monte Carlo repeats +n_part 1000 # total number of particles +restart no # whether to restart from saved state (yes/no) +do_select_weighting yes # whether to select weighting explicitly (yes/no) +weight_type flat + + +t_max 600 # total simulation time (s) +del_t 1 # timestep (s) +t_output 60 # output interval (0 disables) (s) +t_progress 120 # progress printing interval (0 disables) (s) + +do_camp_chem no # whether to run the campible chemistry module +do_tchem no # whether to use TChem for chemistry + +gas_data gas_data.dat # file containing gas data +gas_init gas_init.dat # initial gas mixing ratios + +aerosol_data aero_data.dat # file containing aerosol data +do_fractal no # whether to do fractal treatment +aerosol_init aero_init_dist.dat # aerosol initial condition file + +temp_profile temp.dat # temperature profile file +pressure_profile pressure.dat # pressure profile file +height_profile height.dat # height profile file +gas_emissions gas_emit.dat # gas emissions file +gas_background gas_back.dat # background gas mixing ratios file +aero_emissions aero_emit.dat # aerosol emissions file +aero_background aero_back.dat # aerosol background file +loss_function none # particle loss function + +rel_humidity 1.000 # initial relative humidity (1) +latitude 40 # latitude (degrees, -90 to 90) +longitude 0 # longitude (degrees, -180 to 180) +altitude 0 # altitude (m) +start_time 0 # start time (s since 00:00 UTC) +start_day 1 # start day of year (UTC) + +do_coagulation no # whether to do coagulation (yes/no) +do_condensation no # whether to do condensation (yes/no) +do_mosaic no # whether to do MOSAIC (yes/no) +do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing yes # whether to do freezing (yes/no) +immersion_freezing_scheme ABIFM +do_freezing_naive no # whether to use naive algorithm (yes/no) + +rand_init 0 # random initialization (0 to auto-generate) +allow_doubling yes # whether to allow doubling (yes/no) +allow_halving yes # whether to allow halving (yes/no) +record_removals no # whether to record particle removals (yes/no) +do_parallel no # whether to run in parallel (yes/no) diff --git a/test/freezing/run_part.const.naive.spec b/test/freezing/run_part.const.naive.spec new file mode 100644 index 000000000..76f0558b7 --- /dev/null +++ b/test/freezing/run_part.const.naive.spec @@ -0,0 +1,54 @@ +run_type particle # particle-resolved run +output_prefix out/const_naive/freezing_part # prefix of output files +n_repeat 5 # number of Monte Carlo repeats +n_part 1000 # total number of particles +restart no # whether to restart from saved state (yes/no) +do_select_weighting yes # whether to select weighting explicitly (yes/no) +weight_type flat + + +t_max 600 # total simulation time (s) +del_t 1 # timestep (s) +t_output 60 # output interval (0 disables) (s) +t_progress 120 # progress printing interval (0 disables) (s) + +do_camp_chem no # whether to run the campible chemistry module +do_tchem no # whether to use TChem for chemistry + +gas_data gas_data.dat # file containing gas data +gas_init gas_init.dat # initial gas mixing ratios + +aerosol_data aero_data.dat # file containing aerosol data +do_fractal no # whether to do fractal treatment +aerosol_init aero_init_dist.dat # aerosol initial condition file + +temp_profile temp.dat # temperature profile file +pressure_profile pressure.dat # pressure profile file +height_profile height.dat # height profile file +gas_emissions gas_emit.dat # gas emissions file +gas_background gas_back.dat # background gas mixing ratios file +aero_emissions aero_emit.dat # aerosol emissions file +aero_background aero_back.dat # aerosol background file +loss_function none # particle loss function + +rel_humidity 1.000 # initial relative humidity (1) +latitude 40 # latitude (degrees, -90 to 90) +longitude 0 # longitude (degrees, -180 to 180) +altitude 0 # altitude (m) +start_time 0 # start time (s since 00:00 UTC) +start_day 1 # start day of year (UTC) + +do_coagulation no # whether to do coagulation (yes/no) +do_condensation no # whether to do condensation (yes/no) +do_mosaic no # whether to do MOSAIC (yes/no) +do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing yes # whether to do freezing (yes/no) +immersion_freezing_scheme const +freezing_rate -.01123456789 +do_freezing_naive yes # whether to use naive algorithm (yes/no) + +rand_init 0 # random initialization (0 to auto-generate) +allow_doubling yes # whether to allow doubling (yes/no) +allow_halving yes # whether to allow halving (yes/no) +record_removals no # whether to record particle removals (yes/no) +do_parallel no # whether to run in parallel (yes/no) diff --git a/test/freezing/run_part.const.spec b/test/freezing/run_part.const.spec new file mode 100644 index 000000000..d45d0b89f --- /dev/null +++ b/test/freezing/run_part.const.spec @@ -0,0 +1,54 @@ +run_type particle # particle-resolved run +output_prefix out/const/freezing_part # prefix of output files +n_repeat 5 # number of Monte Carlo repeats +n_part 1000 # total number of particles +restart no # whether to restart from saved state (yes/no) +do_select_weighting yes # whether to select weighting explicitly (yes/no) +weight_type flat + + +t_max 600 # total simulation time (s) +del_t 1 # timestep (s) +t_output 60 # output interval (0 disables) (s) +t_progress 120 # progress printing interval (0 disables) (s) + +do_camp_chem no # whether to run the campible chemistry module +do_tchem no # whether to use TChem for chemistry + +gas_data gas_data.dat # file containing gas data +gas_init gas_init.dat # initial gas mixing ratios + +aerosol_data aero_data.dat # file containing aerosol data +do_fractal no # whether to do fractal treatment +aerosol_init aero_init_dist.dat # aerosol initial condition file + +temp_profile temp.dat # temperature profile file +pressure_profile pressure.dat # pressure profile file +height_profile height.dat # height profile file +gas_emissions gas_emit.dat # gas emissions file +gas_background gas_back.dat # background gas mixing ratios file +aero_emissions aero_emit.dat # aerosol emissions file +aero_background aero_back.dat # aerosol background file +loss_function none # particle loss function + +rel_humidity 1.000 # initial relative humidity (1) +latitude 40 # latitude (degrees, -90 to 90) +longitude 0 # longitude (degrees, -180 to 180) +altitude 0 # altitude (m) +start_time 0 # start time (s since 00:00 UTC) +start_day 1 # start day of year (UTC) + +do_coagulation no # whether to do coagulation (yes/no) +do_condensation no # whether to do condensation (yes/no) +do_mosaic no # whether to do MOSAIC (yes/no) +do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing yes # whether to do freezing (yes/no) +immersion_freezing_scheme const +freezing_rate -.01123456789 +do_freezing_naive no # whether to use naive algorithm (yes/no) + +rand_init 0 # random initialization (0 to auto-generate) +allow_doubling yes # whether to allow doubling (yes/no) +allow_halving yes # whether to allow halving (yes/no) +record_removals no # whether to record particle removals (yes/no) +do_parallel no # whether to run in parallel (yes/no) diff --git a/test/freezing/run_part.singular.spec b/test/freezing/run_part.singular.spec new file mode 100644 index 000000000..eb3e31e06 --- /dev/null +++ b/test/freezing/run_part.singular.spec @@ -0,0 +1,54 @@ +run_type particle # particle-resolved run +output_prefix out/singular/freezing_part # prefix of output files +n_repeat 5 # number of Monte Carlo repeats +n_part 1000 # total number of particles +restart no # whether to restart from saved state (yes/no) +do_select_weighting yes # whether to select weighting explicitly (yes/no) +weight_type flat + + +t_max 600 # total simulation time (s) +del_t 1 # timestep (s) +t_output 60 # output interval (0 disables) (s) +t_progress 120 # progress printing interval (0 disables) (s) + +do_camp_chem no # whether to run the campible chemistry module +do_tchem no # whether to use TChem for chemistry + +gas_data gas_data.dat # file containing gas data +gas_init gas_init.dat # initial gas mixing ratios + +aerosol_data aero_data.dat # file containing aerosol data +do_fractal no # whether to do fractal treatment +aerosol_init aero_init_dist.dat # aerosol initial condition file + +temp_profile temp.dat # temperature profile file +pressure_profile pressure.dat # pressure profile file +height_profile height.dat # height profile file +gas_emissions gas_emit.dat # gas emissions file +gas_background gas_back.dat # background gas mixing ratios file +aero_emissions aero_emit.dat # aerosol emissions file +aero_background aero_back.dat # aerosol background file +loss_function none # particle loss function + +rel_humidity 1.000 # initial relative humidity (1) +latitude 40 # latitude (degrees, -90 to 90) +longitude 0 # longitude (degrees, -180 to 180) +altitude 0 # altitude (m) +start_time 0 # start time (s since 00:00 UTC) +start_day 1 # start day of year (UTC) + +do_coagulation no # whether to do coagulation (yes/no) +do_condensation no # whether to do condensation (yes/no) +do_mosaic no # whether to do MOSAIC (yes/no) +do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing yes # whether to do freezing (yes/no) +immersion_freezing_scheme singular +INAS_a -0.517 +INAS_b 8.934 + +rand_init 0 # random initialization (0 to auto-generate) +allow_doubling yes # whether to allow doubling (yes/no) +allow_halving yes # whether to allow halving (yes/no) +record_removals no # whether to record particle removals (yes/no) +do_parallel no # whether to run in parallel (yes/no) diff --git a/test/freezing/temp.dat b/test/freezing/temp.dat new file mode 100644 index 000000000..45e691790 --- /dev/null +++ b/test/freezing/temp.dat @@ -0,0 +1,4 @@ +# time (s) +# temp (K) +time 0.000 +temp 243.150 diff --git a/test/freezing/test_freezing_1.sh b/test/freezing/test_freezing_1.sh new file mode 100755 index 000000000..3ee5d9ef5 --- /dev/null +++ b/test/freezing/test_freezing_1.sh @@ -0,0 +1,23 @@ +#!/bin/bash + +# exit on error +set -e +# turn on command echoing +set -v +# make sure that the current directory is the one where this script is +cd ${0%/*} + +# make the output directory if it doesn't exist +mkdir -p out/ABIFM + +# run PartMC freezing test +../../partmc run_part.ABIFM.spec + +# extract frozen fraction results to a txt file +../../test_freezing_extract out/ABIFM/freezing_part + +# calculate the freezing theoretical results +../../test_freezing_theoretical ABIFM + +# compare the PartMC versus theoretical results +../../numeric_diff --rel-tol 0.02 out/ABIFM/freezing_part_frozen_fraction_ensemble_mean.txt out/freezing_theoretical_ABIFM_data.txt diff --git a/test/freezing/test_freezing_2.sh b/test/freezing/test_freezing_2.sh new file mode 100755 index 000000000..9a63e2e19 --- /dev/null +++ b/test/freezing/test_freezing_2.sh @@ -0,0 +1,23 @@ +#!/bin/bash + +# exit on error +set -e +# turn on command echoing +set -v +# make sure that the current directory is the one where this script is +cd ${0%/*} + +# make the output directory if it doesn't exist +mkdir -p out/singular + +# run PartMC freezing test +../../partmc run_part.singular.spec + +# extract frozen fraction results to a txt file +../../test_freezing_extract out/singular/freezing_part + +# calculate the freezing theoretical results +../../test_freezing_theoretical singular + +# compare the PartMC versus theoretical results +../../numeric_diff --rel-tol 0.04 out/singular/freezing_part_frozen_fraction_ensemble_mean.txt out/freezing_theoretical_singular_data.txt diff --git a/test/freezing/test_freezing_3.sh b/test/freezing/test_freezing_3.sh new file mode 100755 index 000000000..0dac7d6f8 --- /dev/null +++ b/test/freezing/test_freezing_3.sh @@ -0,0 +1,23 @@ +#!/bin/bash + +# exit on error +set -e +# turn on command echoing +set -v +# make sure that the current directory is the one where this script is +cd ${0%/*} + +# make the output directory if it doesn't exist +mkdir -p out/const + +# run PartMC freezing test +../../partmc run_part.const.spec + +# extract frozen fraction results to a txt file +../../test_freezing_extract out/const/freezing_part + +# calculate the freezing theoretical results +../../test_freezing_theoretical const + +# compare the PartMC versus theoretical results +../../numeric_diff --rel-tol 0.02 out/const/freezing_part_frozen_fraction_ensemble_mean.txt out/freezing_theoretical_const_data.txt diff --git a/test/freezing/test_freezing_4.sh b/test/freezing/test_freezing_4.sh new file mode 100755 index 000000000..6431fc809 --- /dev/null +++ b/test/freezing/test_freezing_4.sh @@ -0,0 +1,23 @@ +#!/bin/bash + +# exit on error +set -e +# turn on command echoing +set -v +# make sure that the current directory is the one where this script is +cd ${0%/*} + +# make the output directory if it doesn't exist +mkdir -p out/ABIFM_naive + +# run PartMC freezing test +../../partmc run_part.ABIFM.naive.spec + +# extract frozen fraction results to a txt file +../../test_freezing_extract out/ABIFM_naive/freezing_part + +# calculate the freezing theoretical results +#../../test_freezing_theoretical ABIFM + +# compare the PartMC versus theoretical results +../../numeric_diff --rel-tol 0.02 out/ABIFM_naive/freezing_part_frozen_fraction_ensemble_mean.txt out/freezing_theoretical_ABIFM_data.txt diff --git a/test/freezing/test_freezing_5.sh b/test/freezing/test_freezing_5.sh new file mode 100755 index 000000000..f2b60b2b8 --- /dev/null +++ b/test/freezing/test_freezing_5.sh @@ -0,0 +1,23 @@ +#!/bin/bash + +# exit on error +set -e +# turn on command echoing +set -v +# make sure that the current directory is the one where this script is +cd ${0%/*} + +# make the output directory if it doesn't exist +mkdir -p out/const_naive + +# run PartMC freezing test +../../partmc run_part.const.naive.spec + +# extract frozen fraction results to a txt file +../../test_freezing_extract out/const_naive/freezing_part + +# calculate the freezing theoretical results +#../../test_freezing_theoretical const + +# compare the PartMC versus theoretical results +../../numeric_diff --rel-tol 0.02 out/const_naive/freezing_part_frozen_fraction_ensemble_mean.txt out/freezing_theoretical_const_data.txt diff --git a/test/freezing/test_freezing_extract.F90 b/test/freezing/test_freezing_extract.F90 new file mode 100644 index 000000000..eb7d584d9 --- /dev/null +++ b/test/freezing/test_freezing_extract.F90 @@ -0,0 +1,144 @@ +!! Copyright (C) 2009-2012 Matthew West +! Licensed under the GNU General Public License version 2 or (at your +! option) any later version. See the file COPYING for details. + +!> \file +!> The extract_freezing program. + +!> Read NetCDF output files and write out the time evolution of +!> frozen fraction average among ensembles in text format. +program extract_freezing + + use pmc_aero_state + use pmc_aero_particle + use pmc_output + use pmc_mpi + use getopt_m + + integer :: n_ensemble + character(len=PMC_MAX_FILENAME_LEN) :: in_prefix, in_prefix_ensemble, out_filename + character(len=PMC_MAX_FILENAME_LEN), allocatable :: filename_list(:) + character(len=1000) :: tmp_str + type(aero_data_t) :: aero_data + type(aero_state_t) :: aero_state + type(env_state_t) :: env_state + integer :: index, i_repeat, i_spec, out_unit + integer :: i_file, n_file, i_ensemble + real(kind=dp) :: time, del_t + character(len=PMC_UUID_LEN) :: uuid, run_uuid + real(kind=dp), allocatable :: frozen_fractions(:),& + frozen_fractions_mean(:) + type(option_s) :: opts(2) + + call pmc_mpi_init() + + opts(1) = option_s("help", .false., 'h') + opts(2) = option_s("output", .true., 'o') + + out_filename = "" + + do + select case(getopt("ho:", opts)) + case(char(0)) + exit + case('h') + call print_help() + stop + case('o') + out_filename = optarg + case( '?' ) + call print_help() + call die_msg(514364550, 'unknown option: ' // trim(optopt)) + case default + call print_help() + call die_msg(603100341, 'unhandled option: ' // trim(optopt)) + end select + end do + + if (optind /= command_argument_count()) then + call print_help() + call die_msg(967032896, 'expected exactly one non-option prefix argument') + end if + + call get_command_argument(optind, in_prefix) + + if (out_filename == "") then + out_filename = trim(in_prefix) // "_frozen_fraction_ensemble_mean.txt" + end if + + + call input_n_files(in_prefix, n_ensemble, n_file) + + allocate(filename_list(0)) + + write(in_prefix_ensemble, "(a,a,i4.4)") trim(in_prefix), "_", 1 + call input_filename_list(in_prefix_ensemble, filename_list) + + allocate(frozen_fractions(n_file)) + allocate(frozen_fractions_mean(n_file)) + frozen_fractions_mean = 0d0 + + do i_ensemble = 1,N_ensemble + write(in_prefix_ensemble, "(a,a,i4.4)") trim(in_prefix), '_', i_ensemble + call input_filename_list(in_prefix_ensemble, filename_list) + n_file = size(filename_list) + call assert_msg(323514871, n_file > 0, & + "no NetCDF files found with prefix: " // trim(in_prefix_ensemble)) + + call input_state(filename_list(1), index, time, del_t, i_repeat, uuid, & + aero_data=aero_data, aero_state=aero_state, env_state=env_state) + run_uuid = uuid + + + + do i_file = 1,n_file + call input_state(filename_list(i_file), index, time, del_t, i_repeat, & + uuid, aero_data=aero_data, aero_state=aero_state, env_state=env_state) + + call assert_msg(397906326, uuid == run_uuid, & + "UUID mismatch between " // trim(filename_list(1)) // " and " & + // trim(filename_list(i_file))) + + frozen_fractions(i_file) = aero_state_frozen_fraction(aero_state, aero_data) + end do + frozen_fractions_mean = frozen_fractions_mean + frozen_fractions + + end do + frozen_fractions_mean = frozen_fractions_mean / N_ensemble + + write(*,'(a,a)') "Output file: ", trim(out_filename) + write(*,'(a)') " Each row of output is one time." + write(*,'(a)') " The columns of output are:" + write(*,'(a)') " column 1: frozen fraction ensemble mean (unitless)" + call open_file_write(out_filename, out_unit) + do i_file = 1,n_file + write(out_unit, '(e30.15e3)', advance='no') frozen_fractions_mean(i_file) + write(out_unit, '(a)') '' + end do + call close_file(out_unit) + + deallocate(frozen_fractions) + deallocate(frozen_fractions_mean) + deallocate(filename_list) + + call pmc_mpi_finalize() + +contains + + subroutine print_help() + + write(*,'(a)') 'Usage: extract_freezing [options] ' + write(*,'(a)') '' + write(*,'(a)') 'options are:' + write(*,'(a)') ' -h, --help Print this help message.' + write(*,'(a)') ' -o, --out Output filename.' + write(*,'(a)') '' + write(*,'(a)') 'Examples:' + write(*,'(a)') ' extract_freezing freezing_part' + write(*,'(a)') '' + + end subroutine print_help + +end program extract_freezing + + diff --git a/test/freezing/test_freezing_theoretical.F90 b/test/freezing/test_freezing_theoretical.F90 new file mode 100644 index 000000000..d585866de --- /dev/null +++ b/test/freezing/test_freezing_theoretical.F90 @@ -0,0 +1,141 @@ +! --- +! This program is used to calculate the theoretical frozen fraction of +! mono-disperse INPs composed of two species in internal mixing at a +! constant temperature during the immersion freezing proces. +! All INPs are assumed to be immersed in supercooled water. +! Calculation of J_het is based on the ABIFM method. + +! References: +! Tang, W. (2024). Particle-resolved simulations of immersion freezing with multi-species ice nucleating particles +! (Master’s thesis, University of Illinois at Urbana-Champaign). https://hdl.handle.net/2142/124611 + +! Knopf, D. A., & Alpert, P. A. (2013). A water activity based model of +! heterogeneous ice nucleation kinetics for freezing of water and aqueous +! solution droplets. Faraday discussions, 165, 513-534. Doi:10.1039/C3FD00035D. + + +program theoretical_freezing + use pmc_env_state + implicit none + + ! Temperature (unit: Kelvin) + real(kind=dp), parameter :: temperature = 243.15 + ! Melting point (unit: Kelvin) + real(kind=dp), parameter :: T0 = 273.15 + ! Time duration (unit: s) + real(kind=dp), parameter :: total_time = 600 + ! Output time interval (unit: s) + real(kind=dp), parameter :: out_dt = 60 + ! ABIFM parameters (m and c) for species 1. + real(kind=dp), parameter :: abifm_m_1 = 17.62106, abifm_c_1 = 1.42411 + ! ABIFM parameters (m and c) for species 2. + real(kind=dp), parameter :: abifm_m_2 = 54.48075, abifm_c_2 = -10.66873 + ! The ratio of total surface for species 1 and 2 (for ABIFM only). + real(kind=dp), parameter :: s_ratio_1 = 0.5, s_ratio_2 = 0.5 + ! INAS parameters (a and b) for mineral dust. (for singular only) + real(kind=dp), parameter :: inas_a = -0.517, inas_b = 8.934 + ! Constant freezing rate. (for const only) + real(kind=dp), parameter :: freezing_rate = -0.01123456789 + ! Dry diameter of INPs (unit: m) + real(kind=dp), parameter :: Dp_dry = 1d-6 + ! The name of output file + !character(len=*), parameter :: out_filename = "out/freezing_theoretical_data.txt" + character(len=100) :: out_filename, imf_scheme + !character(len=*), parameter :: out_dir = "out" + !character(len=*), parameter :: filename = "freezing_theoretical_data.txt" + integer, parameter :: out_unit = 65 + integer, parameter :: IMMERSION_FREEZING_SCHEME_CONST = 0 + integer, parameter :: IMMERSION_FREEZING_SCHEME_SINGULAR = 1 + integer, parameter :: IMMERSION_FREEZING_SCHEME_ABIFM = 2 + real(kind=dp) :: time, Jhet_1, Jhet_2, Jhet_mean, Phi_mean, frozen_fraction + real(kind=dp) :: ns + integer :: ios, immersion_freezing_scheme_type + + call getarg(1, imf_scheme) + if (trim(imf_scheme) == 'const') then + immersion_freezing_scheme_type = IMMERSION_FREEZING_SCHEME_CONST + elseif (trim(imf_scheme) == 'singular') then + immersion_freezing_scheme_type = IMMERSION_FREEZING_SCHEME_SINGULAR + elseif (trim(imf_scheme) == 'ABIFM') then + immersion_freezing_scheme_type = IMMERSION_FREEZING_SCHEME_ABIFM + else + print*, "Unknown immersion freezing scheme: " // trim(imf_scheme) + immersion_freezing_scheme_type = -1 + endif + out_filename = "out/freezing_theoretical_" // trim(imf_scheme) // & + "_data.txt" + ! open output + open(unit=out_unit, file=out_filename, iostat=ios) + if (ios /= 0) then + write(0,'(a,a,a,i4)') 'ERROR: unable to open file ', & + trim(out_filename), ' for writing: ', ios + stop 1 + end if + + if (immersion_freezing_scheme_type == & + IMMERSION_FREEZING_SCHEME_ABIFM) then + call compute_Jhet_ABIFM(temperature, abifm_m_1, abifm_c_1, Jhet_1) + call compute_Jhet_ABIFM(temperature, abifm_m_2, abifm_c_2, Jhet_2) + Jhet_mean = Jhet_1 * s_ratio_1 + Jhet_2 * s_ratio_2 + + time = 0 + do while(time .le. total_time) + Phi_mean = Jhet_mean * time + frozen_fraction = 1 - exp(- const%pi * Dp_dry**2 * Phi_mean) + write(out_unit,'(e20.10)') frozen_fraction + time = time + out_dt + end do + elseif (immersion_freezing_scheme_type == & + IMMERSION_FREEZING_SCHEME_SINGULAR) then + time = 0 + do while(time .le. total_time) + if (time .eq. 0) then + frozen_fraction = 0.0 + else + call compute_INAS(temperature, inas_a, inas_b, ns) + frozen_fraction = 1 - exp(- const%pi * Dp_dry**2 * ns) + end if + write(out_unit,'(e20.10)') frozen_fraction + time = time + out_dt + end do + elseif (immersion_freezing_scheme_type == & + IMMERSION_FREEZING_SCHEME_CONST) then + time = 0 + do while(time .le. total_time) + frozen_fraction = 1 - exp(freezing_rate * time) + write(out_unit,'(e20.10)') frozen_fraction + time = time + out_dt + end do + endif + + close(out_unit) + + contains + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine compute_Jhet_ABIFM(T, abifm_m, abifm_c, Jhet) + ! Calculate J_het value using ABIFM method (Knopf et al., 2013) + use pmc_env_state + use pmc_constants + implicit none + real(kind=dp), intent(in) :: T, abifm_m, abifm_c + real(kind=dp), intent(out) :: Jhet + real(kind=dp) :: es, ei, a_w_ice + es = env_state_saturated_vapor_pressure_wrt_water(T) + ei = env_state_saturated_vapor_pressure_wrt_ice(T) + a_w_ice = ei / es + Jhet = 10 ** (abifm_m * (1 - a_w_ice) + abifm_c) * 10000 + + end subroutine compute_Jhet_ABIFM + + subroutine compute_INAS(T, inas_a, inas_b, ns) + implicit none + ! Calculate the INAS density (Niemand et al., 2012) + real(kind=dp), intent(in) :: T, inas_a, inas_b + real(kind=dp), intent(out) :: ns + ns = exp(inas_a * (T - T0) + inas_b) + end subroutine + +end program theoretical_freezing + + diff --git a/test/loss/aero_data.dat b/test/loss/aero_data.dat index 41cdc37fc..cc430aaee 100644 --- a/test/loss/aero_data.dat +++ b/test/loss/aero_data.dat @@ -1,2 +1,2 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) +H2O 1000 0 18d-3 0 0 0 diff --git a/test/loss/run_chamber_part.spec b/test/loss/run_chamber_part.spec index 3675a6872..59f6c0a5f 100644 --- a/test/loss/run_chamber_part.spec +++ b/test/loss/run_chamber_part.spec @@ -45,6 +45,7 @@ do_coagulation no # whether to do coagulation (yes/no) do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/loss/run_constant_part.spec b/test/loss/run_constant_part.spec index cc469ed6c..06790780d 100644 --- a/test/loss/run_constant_part.spec +++ b/test/loss/run_constant_part.spec @@ -40,6 +40,7 @@ do_coagulation no # whether to do coagulation (yes/no) do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/loss/run_drydep_part.spec b/test/loss/run_drydep_part.spec index 54c24a261..d35c5eab2 100644 --- a/test/loss/run_drydep_part.spec +++ b/test/loss/run_drydep_part.spec @@ -40,6 +40,7 @@ do_coagulation no # whether to do coagulation (yes/no) do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/loss/run_volume_part.spec b/test/loss/run_volume_part.spec index 864127e09..7d8a515f1 100644 --- a/test/loss/run_volume_part.spec +++ b/test/loss/run_volume_part.spec @@ -40,6 +40,7 @@ do_coagulation no # whether to do coagulation (yes/no) do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/mixing_state/aero_data.dat b/test/mixing_state/aero_data.dat index 740f9933b..9ee382e85 100644 --- a/test/mixing_state/aero_data.dat +++ b/test/mixing_state/aero_data.dat @@ -1,5 +1,5 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -A 1000 0 100d-3 0.1 -B 1500 0 100d-3 0.1 -C 2000 0 100d-3 0.1 -D 2500 0 100d-3 0.1 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) +A 1000 0 100d-3 0.1 0 0 +B 1500 0 100d-3 0.1 0 0 +C 2000 0 100d-3 0.1 0 0 +D 2500 0 100d-3 0.1 0 0 diff --git a/test/mixing_state/run_part.spec b/test/mixing_state/run_part.spec index f4832e215..af3df0392 100644 --- a/test/mixing_state/run_part.spec +++ b/test/mixing_state/run_part.spec @@ -41,6 +41,7 @@ coag_kernel brown # coagulation kernel do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to use time) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/mosaic/aero_data.dat b/test/mosaic/aero_data.dat index d6d5f3538..65d773f46 100644 --- a/test/mosaic/aero_data.dat +++ b/test/mosaic/aero_data.dat @@ -1,21 +1,21 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa -H2O 1000 0 18d-3 0 -SO4 1800 0 96d-3 0.53 -NO3 1800 0 62d-3 0.53 -Cl 2200 0 35.5d-3 1.28 -NH4 1800 0 18d-3 0.53 -CO3 2600 0 60d-3 0.53 -MSA 1800 0 95d-3 0.53 -Na 2200 0 23d-3 1.28 -Ca 2600 0 40d-3 0.53 -OC 1000 0 1d-3 0.1 -BC 1800 0 1d-3 0 -OIN 2600 0 1d-3 0.1 -ARO1 1000 0 150d-3 0.1 -ARO2 1000 0 150d-3 0.1 -ALK1 1000 0 140d-3 0.1 -OLE1 1000 0 140d-3 0.1 -API1 1000 0 184d-3 0.1 -API2 1000 0 184d-3 0.1 -LIM1 1000 0 200d-3 0.1 -LIM2 1000 0 200d-3 0.1 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa abifm_m (1) abifm_c (1) +H2O 1000 0 18d-3 0 0 0 +SO4 1800 0 96d-3 0.53 0 0 +NO3 1800 0 62d-3 0.53 0 0 +Cl 2200 0 35.5d-3 1.28 0 0 +NH4 1800 0 18d-3 0.53 0 0 +CO3 2600 0 60d-3 0.53 0 0 +MSA 1800 0 95d-3 0.53 0 0 +Na 2200 0 23d-3 1.28 0 0 +Ca 2600 0 40d-3 0.53 0 0 +OC 1000 0 1d-3 0.1 0 0 +BC 1800 0 1d-3 0 0 0 +OIN 2600 0 1d-3 0.1 0 0 +ARO1 1000 0 150d-3 0.1 0 0 +ARO2 1000 0 150d-3 0.1 0 0 +ALK1 1000 0 140d-3 0.1 0 0 +OLE1 1000 0 140d-3 0.1 0 0 +API1 1000 0 184d-3 0.1 0 0 +API2 1000 0 184d-3 0.1 0 0 +LIM1 1000 0 200d-3 0.1 0 0 +LIM2 1000 0 200d-3 0.1 0 0 diff --git a/test/mosaic/run_part.spec b/test/mosaic/run_part.spec index 6d5426ec8..f065d826f 100644 --- a/test/mosaic/run_part.spec +++ b/test/mosaic/run_part.spec @@ -41,6 +41,7 @@ do_condensation no # whether to do condensation (yes/no) do_mosaic yes # whether to do MOSAIC (yes/no) do_optical yes # whether to compute optical props (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/mosaic/run_part_restarted.spec b/test/mosaic/run_part_restarted.spec index c6f2c5986..c9a9c149c 100644 --- a/test/mosaic/run_part_restarted.spec +++ b/test/mosaic/run_part_restarted.spec @@ -34,6 +34,7 @@ do_condensation no # whether to do condensation (yes/no) do_mosaic yes # whether to do MOSAIC (yes/no) do_optical yes # whether to compute optical props (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/nucleate/aero_data.dat b/test/nucleate/aero_data.dat index 1f0f968fa..3b62e55ae 100644 --- a/test/nucleate/aero_data.dat +++ b/test/nucleate/aero_data.dat @@ -1,3 +1,3 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -SO4 1800 0 96d-3 0.65 -NO3 1800 0 62d-3 0.65 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) +SO4 1800 0 96d-3 0.65 0 0 +NO3 1800 0 62d-3 0.65 0 0 diff --git a/test/nucleate/run_part.spec b/test/nucleate/run_part.spec index b806c2b81..28686527d 100644 --- a/test/nucleate/run_part.spec +++ b/test/nucleate/run_part.spec @@ -41,6 +41,7 @@ do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation yes # whether to do nucleation (yes/no) nucleate sulf_acid # nucleation parameterization +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling no # whether to allow doubling (yes/no) diff --git a/test/parallel/aero_data.dat b/test/parallel/aero_data.dat index 41cdc37fc..cc430aaee 100644 --- a/test/parallel/aero_data.dat +++ b/test/parallel/aero_data.dat @@ -1,2 +1,2 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) +H2O 1000 0 18d-3 0 0 0 diff --git a/test/parallel/run_part_parallel_dist.spec b/test/parallel/run_part_parallel_dist.spec index 210e729ad..182cb5c3e 100644 --- a/test/parallel/run_part_parallel_dist.spec +++ b/test/parallel/run_part_parallel_dist.spec @@ -42,6 +42,7 @@ coag_kernel brown # coagulation kernel do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/parallel/run_part_parallel_dist_single.spec b/test/parallel/run_part_parallel_dist_single.spec index b1a803233..ecd46580b 100644 --- a/test/parallel/run_part_parallel_dist_single.spec +++ b/test/parallel/run_part_parallel_dist_single.spec @@ -42,6 +42,7 @@ coag_kernel brown # coagulation kernel do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/parallel/run_part_parallel_mix.spec b/test/parallel/run_part_parallel_mix.spec index 819299b18..4d8ded8cd 100644 --- a/test/parallel/run_part_parallel_mix.spec +++ b/test/parallel/run_part_parallel_mix.spec @@ -41,6 +41,7 @@ coag_kernel brown # coagulation kernel do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/parallel/run_part_serial.spec b/test/parallel/run_part_serial.spec index c9fb42ea0..0097a6c33 100644 --- a/test/parallel/run_part_serial.spec +++ b/test/parallel/run_part_serial.spec @@ -41,6 +41,7 @@ coag_kernel brown # coagulation kernel do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/sedi/aero_data.dat b/test/sedi/aero_data.dat index 41cdc37fc..5a703b6c9 100644 --- a/test/sedi/aero_data.dat +++ b/test/sedi/aero_data.dat @@ -1,2 +1,2 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -H2O 1000 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) +H2O 1000 0 18d-3 0 0 0 diff --git a/test/sedi/run_part.spec b/test/sedi/run_part.spec index 53e565105..6ecfb120e 100644 --- a/test/sedi/run_part.spec +++ b/test/sedi/run_part.spec @@ -41,6 +41,7 @@ coag_kernel sedi # coagulation kernel do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/tchem/aero_data.dat b/test/tchem/aero_data.dat index 83d7a7091..892067a0d 100644 --- a/test/tchem/aero_data.dat +++ b/test/tchem/aero_data.dat @@ -1,4 +1,4 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -SI 10000 0 18d-3 0 -SB 100 0 18d-3 0 -SE 1500 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) +SI 10000 0 18d-3 0 0 0 +SB 100 0 18d-3 0 0 0 +SE 1500 0 18d-3 0 0 0 diff --git a/test/tchem/run_part_cb05cl_ae5.spec b/test/tchem/run_part_cb05cl_ae5.spec index 296463b2e..0526191c0 100644 --- a/test/tchem/run_part_cb05cl_ae5.spec +++ b/test/tchem/run_part_cb05cl_ae5.spec @@ -42,6 +42,7 @@ do_coagulation no # whether to do coagulation (yes/no) do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/tchem/run_part_chapman.spec b/test/tchem/run_part_chapman.spec index e28a9dae5..91a449619 100644 --- a/test/tchem/run_part_chapman.spec +++ b/test/tchem/run_part_chapman.spec @@ -42,6 +42,7 @@ do_coagulation no # whether to do coagulation (yes/no) do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/weighting/aero_data.dat b/test/weighting/aero_data.dat index 83d7a7091..892067a0d 100644 --- a/test/weighting/aero_data.dat +++ b/test/weighting/aero_data.dat @@ -1,4 +1,4 @@ -# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) -SI 10000 0 18d-3 0 -SB 100 0 18d-3 0 -SE 1500 0 18d-3 0 +# dens (kg/m^3) ions in soln (1) molec wght (kg/mole) kappa (1) abifm_m (1) abifm_c (1) +SI 10000 0 18d-3 0 0 0 +SB 100 0 18d-3 0 0 0 +SE 1500 0 18d-3 0 0 0 diff --git a/test/weighting/run_part_flat.spec b/test/weighting/run_part_flat.spec index 4ae85f318..3e0d24570 100644 --- a/test/weighting/run_part_flat.spec +++ b/test/weighting/run_part_flat.spec @@ -41,6 +41,7 @@ do_coagulation no # whether to do coagulation (yes/no) do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/weighting/run_part_flat_source.spec b/test/weighting/run_part_flat_source.spec index 5ff999a45..1f9839264 100644 --- a/test/weighting/run_part_flat_source.spec +++ b/test/weighting/run_part_flat_source.spec @@ -41,6 +41,7 @@ do_coagulation no # whether to do coagulation (yes/no) do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/weighting/run_part_flat_specified.spec b/test/weighting/run_part_flat_specified.spec index 3bb804688..1e0901c64 100644 --- a/test/weighting/run_part_flat_specified.spec +++ b/test/weighting/run_part_flat_specified.spec @@ -41,6 +41,7 @@ do_coagulation no # whether to do coagulation (yes/no) do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/weighting/run_part_nummass.spec b/test/weighting/run_part_nummass.spec index eb6a5d208..c402c93d5 100644 --- a/test/weighting/run_part_nummass.spec +++ b/test/weighting/run_part_nummass.spec @@ -41,6 +41,7 @@ do_coagulation no # whether to do coagulation (yes/no) do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/weighting/run_part_nummass_source.spec b/test/weighting/run_part_nummass_source.spec index fd107ebab..3d70d4cdd 100644 --- a/test/weighting/run_part_nummass_source.spec +++ b/test/weighting/run_part_nummass_source.spec @@ -41,6 +41,7 @@ do_coagulation no # whether to do coagulation (yes/no) do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no) diff --git a/test/weighting/run_part_nummass_specified.spec b/test/weighting/run_part_nummass_specified.spec index 0dc3e5cca..ab8129f67 100644 --- a/test/weighting/run_part_nummass_specified.spec +++ b/test/weighting/run_part_nummass_specified.spec @@ -41,6 +41,7 @@ do_coagulation no # whether to do coagulation (yes/no) do_condensation no # whether to do condensation (yes/no) do_mosaic no # whether to do MOSAIC (yes/no) do_nucleation no # whether to do nucleation (yes/no) +do_immersion_freezing no # whether to do freezing (yes/no) rand_init 0 # random initialization (0 to auto-generate) allow_doubling yes # whether to allow doubling (yes/no)