From d10c23b39b28b77d8586c5dcd0839a99138bade7 Mon Sep 17 00:00:00 2001 From: Agnieszka Makulska Date: Thu, 17 Oct 2024 12:52:15 +0200 Subject: [PATCH 01/17] formulas for ice A nucleation --- include/libcloudph++/blk_1m/formulae.hpp | 49 ++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/include/libcloudph++/blk_1m/formulae.hpp b/include/libcloudph++/blk_1m/formulae.hpp index 4375b705..53d57645 100644 --- a/include/libcloudph++/blk_1m/formulae.hpp +++ b/include/libcloudph++/blk_1m/formulae.hpp @@ -92,6 +92,55 @@ namespace libcloudphxx * sqrt(rhod_0 / rhod) ); } + + //ice A nucleation + template + quantity::type, real_t> ice_A_nucleation( + const quantity &ra, + const quantity &rc, + const quantity &rv, + const quantity &rvs, + const quantity &rvsi, + const quantity &T, + const quantity &rhod_0, + const quantity &dt + ) { + const quantity taunuc = dt; // time scale for nucleation + //homogeneous part + if (T < real_t(233.16) * si::kelvins) { // only if T < -40 C + real_t beta = (T > real_t(213.16)*si::kelvins) ? real_t(0.1)+real_t(0.9)*(T-real_t(213.16)*si::kelvins)/(real_t(20)*si::kelvins) : real_t(0.1); + quantity rv_adj = beta * rvs + (1-beta) * rvsi; + //nucleation rates: + quantity::type, real_t> hom1 = (real_t(1)-std::exp(-dt/taunuc)) * std::max(real_t(0), rv-rv_adj) / dt; + quantity::type, real_t> hom2 = (real_t(1)-std::exp(-dt/taunuc)) * rc / dt; + } + else { + quantity::type, real_t> hom1 = real_t(0); + quantity::type, real_t> hom2 = real_t(0); + } + //heterogeneous part + //calculating the mass of small ice A particle: + real_t IWCS = std::min(real_t(1e-3), rhod_0*ra/si::mass_density, real_t(2.52e-4)*std::pow(real_t(1e3)*rhod_0*ra/si::mass_density,real_t(0.837))); + real_t rho_i = real_t(916.8); //ice density + real_t alpha = std::max(real_t(2.2e4), real_t(-4.99e3)-real_t(4.94e4)*std::log10(real_t(1e3)*IWCS)); + quantity m_as = real_t(2*std::numbers::pi)*rho_i/std::pow(alpha,3) *si::mass; + //calculating the mass of large ice A particle: + real_t IWCL = rhod_0*ra - IWCS; + real_t a_mu = real_t(5.2) + real_t(1.3e-3)*(T/si::kelvin-273.16); + real_t b_mu = real_t(0.026) - real_t(1.2e-3)*(T/si::kelvin-273.16); + real_t a_sigma = real_t(0.47) + real_t(2.1e-3)*(T/si::kelvin-273.16); + real_t b_sigma = real_t(0.018) - real_t(2.1e-4)*(T/si::kelvin-273.16); + real_t mu = a_mu + b_mu*std::log10(real_t(1e3)*IWCL); + real_t sigma = a_sigma + b_sigma*std::log10(real_t(1e3)*IWCL); + quantity m_al = real_t(1.67e17*std::numbers::pi)*rho_i*std::exp(3*mu +9/2*std::pow(sigma,2))*si::mass; + //concentration of ice nuclei: + quantity::type, real_t> N_in = std::min(real_t(1e5), real_t(1e-2)*std::exp(real_t(0.6)*(real_t(273.16)-T/si::kelvins))); + //nucleation rate: + quantity::type, real_t> het = (real_t(1)-std::exp(-dt/taunuc)) * std::min(rc, std::max(real_t(0), N_in*std::max(real_t(1e-12)*si::mass, ma)/rhod_0 - ra)) / dt; + inia=het+hom1+hom2; + return inia; + } + }; }; }; From 18a2e4b5d138b9ba563ca95ea6f6aecfb33f4445 Mon Sep 17 00:00:00 2001 From: Agnieszka Makulska Date: Mon, 21 Oct 2024 10:39:17 +0200 Subject: [PATCH 02/17] nucleation rates homa1, homa1, heta --- include/libcloudph++/blk_1m/formulae.hpp | 47 +++++++++++++++++++----- 1 file changed, 38 insertions(+), 9 deletions(-) diff --git a/include/libcloudph++/blk_1m/formulae.hpp b/include/libcloudph++/blk_1m/formulae.hpp index 53d57645..d5f7f61c 100644 --- a/include/libcloudph++/blk_1m/formulae.hpp +++ b/include/libcloudph++/blk_1m/formulae.hpp @@ -93,16 +93,13 @@ namespace libcloudphxx ); } - //ice A nucleation + //ice A homogeneous nucleation 1 template - quantity::type, real_t> ice_A_nucleation( - const quantity &ra, - const quantity &rc, + quantity::type, real_t> hom_A_nucleation_1( const quantity &rv, const quantity &rvs, const quantity &rvsi, const quantity &T, - const quantity &rhod_0, const quantity &dt ) { const quantity taunuc = dt; // time scale for nucleation @@ -110,14 +107,45 @@ namespace libcloudphxx if (T < real_t(233.16) * si::kelvins) { // only if T < -40 C real_t beta = (T > real_t(213.16)*si::kelvins) ? real_t(0.1)+real_t(0.9)*(T-real_t(213.16)*si::kelvins)/(real_t(20)*si::kelvins) : real_t(0.1); quantity rv_adj = beta * rvs + (1-beta) * rvsi; - //nucleation rates: + //nucleation rate: quantity::type, real_t> hom1 = (real_t(1)-std::exp(-dt/taunuc)) * std::max(real_t(0), rv-rv_adj) / dt; - quantity::type, real_t> hom2 = (real_t(1)-std::exp(-dt/taunuc)) * rc / dt; } else { quantity::type, real_t> hom1 = real_t(0); + } + return hom1; + } + + + //ice A homogeneous nucleation 2 + template + quantity::type, real_t> hom_A_nucleation_2( + const quantity &rc, + const quantity &T, + const quantity &dt + ) { + const quantity taunuc = dt; // time scale for nucleation + //homogeneous part + if (T < real_t(233.16) * si::kelvins) { // only if T < -40 C + //nucleation rate: + quantity::type, real_t> hom2 = (real_t(1)-std::exp(-dt/taunuc)) * rc / dt; + } + else { quantity::type, real_t> hom2 = real_t(0); } + return hom2; + } + + //ice A heterogeneous nucleation + template + quantity::type, real_t> het_A_nucleation( + const quantity &ra, + const quantity &rc, + const quantity &T, + const quantity &rhod_0, + const quantity &dt + ) { + const quantity taunuc = dt; // time scale for nucleation //heterogeneous part //calculating the mass of small ice A particle: real_t IWCS = std::min(real_t(1e-3), rhod_0*ra/si::mass_density, real_t(2.52e-4)*std::pow(real_t(1e3)*rhod_0*ra/si::mass_density,real_t(0.837))); @@ -133,12 +161,13 @@ namespace libcloudphxx real_t mu = a_mu + b_mu*std::log10(real_t(1e3)*IWCL); real_t sigma = a_sigma + b_sigma*std::log10(real_t(1e3)*IWCL); quantity m_al = real_t(1.67e17*std::numbers::pi)*rho_i*std::exp(3*mu +9/2*std::pow(sigma,2))*si::mass; + real_t delta = IWCS/(ra*rhod_0/si::mass_density); + quantity ma = delta*m_as + (real_t(1)-delta)*m_al; //concentration of ice nuclei: quantity::type, real_t> N_in = std::min(real_t(1e5), real_t(1e-2)*std::exp(real_t(0.6)*(real_t(273.16)-T/si::kelvins))); //nucleation rate: quantity::type, real_t> het = (real_t(1)-std::exp(-dt/taunuc)) * std::min(rc, std::max(real_t(0), N_in*std::max(real_t(1e-12)*si::mass, ma)/rhod_0 - ra)) / dt; - inia=het+hom1+hom2; - return inia; + return het; } }; From 6321d9ffb900626c2e400abdff7dd20b3e220445 Mon Sep 17 00:00:00 2001 From: Agnieszka Makulska Date: Wed, 23 Oct 2024 16:18:47 +0200 Subject: [PATCH 03/17] ice A nucleation --- bindings/python/blk_1m.hpp | 4 +- cmake/FindThrust.cmake | 66 ----------------- cmake/LICENSE_VTKm.txt | 78 -------------------- include/libcloudph++/blk_1m/formulae.hpp | 16 ++-- include/libcloudph++/blk_1m/rhs_cellwise.hpp | 66 ++++++++++++++--- 5 files changed, 69 insertions(+), 161 deletions(-) delete mode 100644 cmake/FindThrust.cmake delete mode 100644 cmake/LICENSE_VTKm.txt diff --git a/bindings/python/blk_1m.hpp b/bindings/python/blk_1m.hpp index b20aa3f9..6329071f 100644 --- a/bindings/python/blk_1m.hpp +++ b/bindings/python/blk_1m.hpp @@ -106,8 +106,10 @@ namespace libcloudphxx const b1m::opts_t &opts, bp_array &dot_rc, bp_array &dot_rr, + bp_array &dot_ra, const bp_array &rc, - const bp_array &rr + const bp_array &rr, + const bp_array &ra, ) { arr_t diff --git a/cmake/FindThrust.cmake b/cmake/FindThrust.cmake deleted file mode 100644 index 38e8fc47..00000000 --- a/cmake/FindThrust.cmake +++ /dev/null @@ -1,66 +0,0 @@ -##============================================================================= -## -## Copyright (c) Kitware, Inc. -## All rights reserved. -## See LICENSE_VTKm.txt for details. -## -## This software is distributed WITHOUT ANY WARRANTY; without even -## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -## PURPOSE. See the above copyright notice for more information. -## -## Copyright 2012 Sandia Corporation. -## Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, -## the U.S. Government retains certain rights in this software. -## -##============================================================================= - -# -# FindThrust -# -# This module finds the Thrust header files and extrats their version. It -# sets the following variables. -# -# THRUST_INCLUDE_DIR - Include directory for thrust header files. (All header -# files will actually be in the thrust subdirectory.) -# THRUST_VERSION - Version of thrust in the form "major.minor.patch". -# - -find_path( THRUST_INCLUDE_DIR - HINTS - ${THRUST_INCLUDE_DIR} - ${CUDA_INCLUDE_DIRS} - /usr/include/cuda - /usr/local/include - /usr/local/cuda/include - NAMES thrust/version.h - DOC "Thrust headers" - ) -if( THRUST_INCLUDE_DIR ) - list( REMOVE_DUPLICATES THRUST_INCLUDE_DIR ) - include_directories(BEFORE ${THRUST_INCLUDE_DIR} ) -endif() - -# Find thrust version -file( STRINGS ${THRUST_INCLUDE_DIR}/thrust/version.h - version - REGEX "#define THRUST_VERSION[ \t]+([0-9x]+)" - ) -string( REGEX REPLACE - "#define THRUST_VERSION[ \t]+" - "" - version - "${version}" - ) - -string( REGEX MATCH "^[0-9]" major ${version} ) -string( REGEX REPLACE "^${major}00" "" version "${version}" ) -string( REGEX MATCH "^[0-9]" minor ${version} ) -string( REGEX REPLACE "^${minor}0" "" version "${version}" ) -set( THRUST_VERSION "${major}.${minor}.${version}") - -# Check for required components -include( FindPackageHandleStandardArgs ) -find_package_handle_standard_args( Thrust - REQUIRED_VARS THRUST_INCLUDE_DIR - VERSION_VAR THRUST_VERSION - ) diff --git a/cmake/LICENSE_VTKm.txt b/cmake/LICENSE_VTKm.txt deleted file mode 100644 index 9750502b..00000000 --- a/cmake/LICENSE_VTKm.txt +++ /dev/null @@ -1,78 +0,0 @@ -VTKm License Version 1.0 -======================================================================== - -Copyright (c) 2014, -Sandia Corporation, Los Alamos National Security, LLC., -UT-Battelle, LLC., Kitware Inc., University of California Davis -All rights reserved. - -Sandia National Laboratories, New Mexico -PO Box 5800 -Albuquerque, NM 87185 -USA - -UT-Battelle -1 Bethel Valley Rd -Oak Ridge, TN 37830 - -Los Alamos National Security, LLC -105 Central Park Square -Los Alamos, NM 87544 - -Kitware Inc. -28 Corporate Drive -Clifton Park, NY 12065 -USA - -University of California, Davis -One Shields Avenue -Davis, CA 95616 -USA - -Under the terms of Contract DE-AC04-94AL85000, there is a -non-exclusive license for use of this work by or on behalf of the -U.S. Government. - -Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National -Laboratory (LANL), the U.S. Government retains certain rights in -this software. - -Redistribution and use in source and binary forms, with or without -modification, are permitted provided that the following conditions are -met: - - * Redistributions of source code must retain the above copyright - notice, this list of conditions and the following disclaimer. - - * Redistributions in binary form must reproduce the above copyright - notice, this list of conditions and the following disclaimer in the - documentation and/or other materials provided with the - distribution. - - * Neither the name of Kitware nor the names of any contributors may - be used to endorse or promote products derived from this software - without specific prior written permission. - -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR -A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR -CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -======================================================================== - -The following files and directories come from third parties. Check the -contents of these for details on the specifics of their respective -licenses. -- - - - - - - - - - - - - - - - - - - - - - - - do not remove this line -CMake/CheckCXX11Features.cmake -CMake/FindBoostHeaders.cmake -CMake/FindTBB.cmake -CMake/FindGLEW.cmake -vtkm/cont/tbb/internal/parallel_sort.h -vtkm/testing/OptionParser.h diff --git a/include/libcloudph++/blk_1m/formulae.hpp b/include/libcloudph++/blk_1m/formulae.hpp index d5f7f61c..c1b6d54c 100644 --- a/include/libcloudph++/blk_1m/formulae.hpp +++ b/include/libcloudph++/blk_1m/formulae.hpp @@ -99,9 +99,11 @@ namespace libcloudphxx const quantity &rv, const quantity &rvs, const quantity &rvsi, - const quantity &T, + const quantity &th, + const quantity &p, const quantity &dt ) { + const quantity T = std::pow((p/real_t(100000)/si::pascal), 0.286); const quantity taunuc = dt; // time scale for nucleation //homogeneous part if (T < real_t(233.16) * si::kelvins) { // only if T < -40 C @@ -116,14 +118,15 @@ namespace libcloudphxx return hom1; } - //ice A homogeneous nucleation 2 template quantity::type, real_t> hom_A_nucleation_2( const quantity &rc, - const quantity &T, + const quantity &th, + const quantity &p, const quantity &dt ) { + const quantity T = std::pow((p/real_t(100000)/si::pascal), 0.286); const quantity taunuc = dt; // time scale for nucleation //homogeneous part if (T < real_t(233.16) * si::kelvins) { // only if T < -40 C @@ -141,10 +144,12 @@ namespace libcloudphxx quantity::type, real_t> het_A_nucleation( const quantity &ra, const quantity &rc, - const quantity &T, + const quantity &th, + const quantity &p, const quantity &rhod_0, const quantity &dt ) { + const quantity T = std::pow((p/real_t(100000)/si::pascal), 0.286); const quantity taunuc = dt; // time scale for nucleation //heterogeneous part //calculating the mass of small ice A particle: @@ -160,7 +165,8 @@ namespace libcloudphxx real_t b_sigma = real_t(0.018) - real_t(2.1e-4)*(T/si::kelvin-273.16); real_t mu = a_mu + b_mu*std::log10(real_t(1e3)*IWCL); real_t sigma = a_sigma + b_sigma*std::log10(real_t(1e3)*IWCL); - quantity m_al = real_t(1.67e17*std::numbers::pi)*rho_i*std::exp(3*mu +9/2*std::pow(sigma,2))*si::mass; + quantity m_al = real_t(1.67e17*pi())*rho_i*std::exp(3*mu +9/2*std::pow(sigma,2))*si::mass; + //mass of average ice A particle real_t delta = IWCS/(ra*rhod_0/si::mass_density); quantity ma = delta*m_as + (real_t(1)-delta)*m_al; //concentration of ice nuclei: diff --git a/include/libcloudph++/blk_1m/rhs_cellwise.hpp b/include/libcloudph++/blk_1m/rhs_cellwise.hpp index c87dc5ce..e3c54a45 100644 --- a/include/libcloudph++/blk_1m/rhs_cellwise.hpp +++ b/include/libcloudph++/blk_1m/rhs_cellwise.hpp @@ -19,25 +19,37 @@ namespace libcloudphxx const opts_t &opts, cont_t &dot_rc_cont, cont_t &dot_rr_cont, + cont_t &dot_ra_cont, const cont_t &rc_cont, - const cont_t &rr_cont + const cont_t &rr_cont, + const cont_t &ra_cont, + const cont_t &th_cont, + const cont_t &p_cont, + const cont_t &rhod_cont, + const real_t &dt ) // { - for (auto tup : zip(dot_rc_cont, dot_rr_cont, rc_cont, rr_cont)) + for (auto tup : zip(dot_rc_cont, dot_rr_cont, dot_ra_cont, rc_cont, rr_cont, ra_cont, th_cont, p_cont, rhod_cont)) { real_t - tmp = 0, + rc_to_rr = 0, + rc_to_ra = 0, &dot_rc = std::get<0>(tup), - &dot_rr = std::get<1>(tup); + &dot_rr = std::get<1>(tup), + &dot_ra = std::get<2>(tup); const real_t - &rc = std::get<2>(tup), - &rr = std::get<3>(tup); + &rc = std::get<3>(tup), + &rr = std::get<4>(tup), + &ra = std::get<5>(tup), + &th = std::get<6>(tup), + &p = std::get<7>(tup), + &rhod = std::get<8>(tup); // autoconversion if (opts.conv) { - tmp += ( + rc_to_rr += ( formulae::autoconversion_rate( rc * si::dimensionless(), opts.r_c0 * si::dimensionless(), @@ -49,7 +61,7 @@ namespace libcloudphxx // collection if (opts.accr) { - tmp += ( + rc_to_rr += ( formulae::collection_rate( rc * si::dimensionless(), rr * si::dimensionless() @@ -57,8 +69,38 @@ namespace libcloudphxx ); } - dot_rr += tmp; - dot_rc -= tmp; + // ice A heterogeneous nucleation + //if (opts.hetA) + //{ + rc_to_ra += ( + formulae::het_A_nucleation( + ra * si::dimensionless(), + rc * si::dimensionless(), + th * si::kelvin, + p * si::pascal, + rhod * si::mass_density, + dt * si::time + ) * si::seconds // to make it dimensionless + ); + //} + + // ice A homogeneous nucleation 2 + //if (opts.homA2) + //{ + rc_to_ra += ( + formulae::hom_A_nucleation_2( + rc * si::dimensionless(), + th * si::kelvin, + p * si::pascal, + dt * si::time + ) * si::seconds // to make it dimensionless + ); + //} + + dot_rr += rc_to_rr; + dot_rc -= rc_to_rr; + dot_rc -= rc_to_ra; + dot_ra += rc_to_ra; } } @@ -69,16 +111,18 @@ namespace libcloudphxx cont_t &dot_rv_cont, cont_t &dot_rc_cont, cont_t &dot_rr_cont, + cont_t &dot_ra_cont, const cont_t &rhod_cont, const cont_t &p_cont, const cont_t &th_cont, const cont_t &rv_cont, const cont_t &rc_cont, const cont_t &rr_cont, + const cont_t &ra_cont, const real_t &dt // time step in seconds ) { - rhs_cellwise(opts, dot_rc_cont, dot_rr_cont, rc_cont, rr_cont); + rhs_cellwise(opts, dot_rc_cont, dot_rr_cont, dot_ra_cont, rc_cont, rr_cont, ra_cont, th_cont, p_cont, rhod_cont, dt); //calculate T from theta // rain evaporation treated as a force in Newthon-Raphson saturation adjustment for (auto tup : zip(dot_th_cont, dot_rv_cont, dot_rr_cont, rhod_cont, p_cont, th_cont, rv_cont, rr_cont)) From 54512a85b106c017572776cf0ad0429809c7482d Mon Sep 17 00:00:00 2001 From: Agnieszka Makulska Date: Mon, 28 Oct 2024 14:14:22 +0100 Subject: [PATCH 04/17] rhs_cellwise_ice function and testing ice nucleation --- CMakeLists.txt | 130 +++++++------- bindings/python/blk_1m.hpp | 58 ++++++- bindings/python/lib.cpp | 7 +- include/libcloudph++/blk_1m/adj_cellwise.hpp | 8 +- include/libcloudph++/blk_1m/formulae.hpp | 68 +++----- include/libcloudph++/blk_1m/options.hpp | 6 +- include/libcloudph++/blk_1m/rhs_cellwise.hpp | 169 +++++++++++++------ include/libcloudph++/common/const_cp.hpp | 40 ++++- include/libcloudph++/common/moist_air.hpp | 1 + tests/python/unit/api_blk_1m.py | 21 ++- 10 files changed, 324 insertions(+), 184 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 00d9fa58..7f8f1433 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -260,20 +260,20 @@ endif() ############################################################################################ # Thrust, location of Thrust can be hinted by setting THRUST_INCLUDE_DIR -find_package(Thrust) - -if (NOT THRUST_INCLUDE_DIR) - message(FATAL_ERROR "Thrust library not found. - -* To install Thrust, please try: -* Debian/Ubuntu: sudo apt-get install libthrust-dev -* Fedora: TODO -* Homebrew: TODO - ") -endif() +find_package(Thrust REQUIRED) + +#if (NOT THRUST_INCLUDE_DIR) +# message(FATAL_ERROR "Thrust library not found. +# +#* To install Thrust, please try: +#* Debian/Ubuntu: sudo apt-get install libthrust-dev +#* Fedora: TODO +#* Homebrew: TODO +# ") +#endif() # include thrust found here instead of the one shipped with cuda -target_include_directories(cloudphxx_lgrngn PRIVATE ${THRUST_INCLUDE_DIR}) +#target_include_directories(cloudphxx_lgrngn PRIVATE ${THRUST_INCLUDE_DIR}) ############################################################################################ # Boost libraries @@ -296,59 +296,59 @@ endif() ############################################################################################ # BOOST ODEINT VERSION TEST -message(STATUS "Testing if Boost ODEINT version >= 1.58") -set(pfx "boost odeint check") -execute_process(COMMAND "mktemp" "-d" "/tmp/tmp.XXX" RESULT_VARIABLE status OUTPUT_VARIABLE tmpdir) -if (NOT status EQUAL 0) - message(FATAL_ERROR "${pfx}: mkdtemp failed") -endif() -file(WRITE "${tmpdir}/test.cpp" " - #define THRUST_DEVICE_SYSTEM THRUST_DEVICE_SYSTEM_CPP - - #include - - #include - #include - - struct rhs - { - void operator()( - const thrust::cpp::vector &psi, - thrust::cpp::vector &dot_psi, - const float /* t */ - ) - { - assert(psi.size() == dot_psi.size()); - } - }; - - int main() - { - boost::numeric::odeint::euler< - thrust::cpp::vector, // state_type - float, // value_type - thrust::cpp::vector, // deriv_type - float, // time_type - boost::numeric::odeint::thrust_algebra, - boost::numeric::odeint::thrust_operations - > chem_stepper; - - thrust::cpp::vector v(2); - chem_stepper.do_step(rhs(), v, 0, 1); - } -") -execute_process(COMMAND "${CMAKE_CXX_COMPILER}" "-std=c++14" "test.cpp" "-I${Boost_INCLUDE_DIR}" "-I${THRUST_INCLUDE_DIR}" WORKING_DIRECTORY ${tmpdir} RESULT_VARIABLE status ERROR_VARIABLE msg) -if (NOT status EQUAL 0) - message(FATAL_ERROR "${pfx}: c++ compiler failed\n ${msg}") -endif() -execute_process(COMMAND "./a.out" WORKING_DIRECTORY ${tmpdir} RESULT_VARIABLE status OUTPUT_VARIABLE msg) -if (NOT status EQUAL 0) - message(FATAL_ERROR "${pfx}: test program failed, install Boost odeint version >= 1.58") -endif() -unset(pfx) -unset(tmpdir) -unset(msg) -unset(status) +#message(STATUS "Testing if Boost ODEINT version >= 1.58") +#set(pfx "boost odeint check") +#execute_process(COMMAND "mktemp" "-d" "/tmp/tmp.XXX" RESULT_VARIABLE status OUTPUT_VARIABLE tmpdir) +#if (NOT status EQUAL 0) +# message(FATAL_ERROR "${pfx}: mkdtemp failed") +#endif() +#file(WRITE "${tmpdir}/test.cpp" " +# #define THRUST_DEVICE_SYSTEM THRUST_DEVICE_SYSTEM_CPP +# +# #include +# +# #include +# #include +# +# struct rhs +# { +# void operator()( +# const thrust::cpp::vector &psi, +# thrust::cpp::vector &dot_psi, +# const float /* t */ +# ) +# { +# assert(psi.size() == dot_psi.size()); +# } +# }; +# +# int main() +# { +# boost::numeric::odeint::euler< +# thrust::cpp::vector, // state_type +# float, // value_type +# thrust::cpp::vector, // deriv_type +# float, // time_type +# boost::numeric::odeint::thrust_algebra, +# boost::numeric::odeint::thrust_operations +# > chem_stepper; +# +# thrust::cpp::vector v(2); +# chem_stepper.do_step(rhs(), v, 0, 1); +# } +#") +#execute_process(COMMAND "${CMAKE_CXX_COMPILER}" "-std=c++14" "test.cpp" "-I${Boost_INCLUDE_DIR}" "-I${THRUST_INCLUDE_DIR}" WORKING_DIRECTORY ${tmpdir} RESULT_VARIABLE status ERROR_VARIABLE msg) +#if (NOT status EQUAL 0) +# message(FATAL_ERROR "${pfx}: c++ compiler failed\n ${msg}") +#endif() +#execute_process(COMMAND "./a.out" WORKING_DIRECTORY ${tmpdir} RESULT_VARIABLE status OUTPUT_VARIABLE msg) +#if (NOT status EQUAL 0) +# message(FATAL_ERROR "${pfx}: test program failed, install Boost odeint version >= 1.58") +#endif() +#unset(pfx) +#unset(tmpdir) +#unset(msg) +#unset(status) # generate a header file with git revision id if (EXISTS "${CMAKE_SOURCE_DIR}/.git") diff --git a/bindings/python/blk_1m.hpp b/bindings/python/blk_1m.hpp index 6329071f..1d918900 100644 --- a/bindings/python/blk_1m.hpp +++ b/bindings/python/blk_1m.hpp @@ -43,7 +43,7 @@ namespace libcloudphxx np2bz_th, np2bz_rv, np2bz_rc, - np2bz_rr, + np2bz_rr, dt ); } @@ -100,27 +100,64 @@ namespace libcloudphxx dt ); } + + template + void rhs_cellwise( + const b1m::opts_t &opts, + bp_array &dot_rc, + bp_array &dot_rr, + const bp_array &rc, + const bp_array &rr + ) + { + arr_t + np2bz_dot_rc(np2bz(dot_rc)), + np2bz_dot_rr(np2bz(dot_rr)); + b1m::rhs_cellwise( + opts, + np2bz_dot_rc, + np2bz_dot_rr, + np2bz(rc), + np2bz(rr) + ); + } template - void rhs_cellwise( + void rhs_cellwise_ice( const b1m::opts_t &opts, bp_array &dot_rc, bp_array &dot_rr, - bp_array &dot_ra, + bp_array &dot_rv, + bp_array &dot_ria, const bp_array &rc, const bp_array &rr, - const bp_array &ra, + const bp_array &rv, + const bp_array &ria, + const bp_array &theta, + const bp_array &p, + const bp_array &rhod, + const typename arr_t::T_numtype &dt ) { arr_t np2bz_dot_rc(np2bz(dot_rc)), - np2bz_dot_rr(np2bz(dot_rr)); - b1m::rhs_cellwise( + np2bz_dot_rr(np2bz(dot_rr)), + np2bz_dot_rv(np2bz(dot_rv)), + np2bz_dot_ria(np2bz(dot_ria)); + b1m::rhs_cellwise_ice( opts, np2bz_dot_rc, np2bz_dot_rr, + np2bz_dot_rv, + np2bz_dot_ria, np2bz(rc), - np2bz(rr) + np2bz(rr), + np2bz(rv), + np2bz(ria), + np2bz(theta), + np2bz(p), + np2bz(rhod), + dt ); } @@ -131,16 +168,19 @@ namespace libcloudphxx bp_array &dot_rv, bp_array &dot_rc, bp_array &dot_rr, + bp_array &dot_ria, const bp_array &rhod, const bp_array &p, const bp_array &th, const bp_array &rv, const bp_array &rc, const bp_array &rr, + const bp_array &ria, const typename arr_t::T_numtype &dt ) { arr_t + np2bz_dot_ria(np2bz(dot_ria)), np2bz_dot_rc(np2bz(dot_rc)), np2bz_dot_rr(np2bz(dot_rr)), np2bz_dot_rv(np2bz(dot_rv)), @@ -151,13 +191,15 @@ namespace libcloudphxx np2bz_dot_rv, np2bz_dot_rc, np2bz_dot_rr, + np2bz_dot_ria, np2bz(rhod), np2bz(p), np2bz(th), np2bz(rv), np2bz(rc), np2bz(rr), - dt + np2bz(ria), + dt ); } diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index aaed65dd..d7adc5f4 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -155,6 +155,10 @@ BOOST_PYTHON_MODULE(libcloudphxx) .def_readwrite("conv", &b1m::opts_t::conv) .def_readwrite("accr", &b1m::opts_t::accr) .def_readwrite("sedi", &b1m::opts_t::sedi) + .def_readwrite("ice", &b1m::opts_t::ice) + .def_readwrite("homA1", &b1m::opts_t::homA1) + .def_readwrite("homA2", &b1m::opts_t::homA2) + .def_readwrite("hetA", &b1m::opts_t::hetA) .def_readwrite("r_c0", &b1m::opts_t::r_c0) .def_readwrite("r_eps", &b1m::opts_t::r_eps) .def_readwrite("nwtrph_iters", &b1m::opts_t::nwtrph_iters) @@ -162,7 +166,8 @@ BOOST_PYTHON_MODULE(libcloudphxx) bp::def("adj_cellwise", blk_1m::adj_cellwise); bp::def("adj_cellwise_constp", blk_1m::adj_cellwise_constp); bp::def("adj_cellwise_nwtrph", blk_1m::adj_cellwise_nwtrph); - bp::def("rhs_cellwise", blk_1m::rhs_cellwise); + bp::def("rhs_cellwise", blk_1m::rhs_cellwise); + bp::def("rhs_cellwise_ice", blk_1m::rhs_cellwise_ice); bp::def("rhs_cellwise_nwtrph", blk_1m::rhs_cellwise_nwtrph); bp::def("rhs_columnwise", blk_1m::rhs_columnwise); // TODO: handle the returned flux } diff --git a/include/libcloudph++/blk_1m/adj_cellwise.hpp b/include/libcloudph++/blk_1m/adj_cellwise.hpp index 8ccd68d5..5b8f22e6 100644 --- a/include/libcloudph++/blk_1m/adj_cellwise.hpp +++ b/include/libcloudph++/blk_1m/adj_cellwise.hpp @@ -106,7 +106,7 @@ namespace libcloudphxx &rc = std::get<3>(tup); // double-checking.... - assert(th >= 273.15); // TODO: that's theta, not T! + //assert(th >= 273.15); // TODO: that's theta, not T! assert(rc >= 0); assert(rv >= 0); @@ -150,7 +150,7 @@ namespace libcloudphxx th += L0 / (moist_air::c_pd() * exner_p) * drc / si::kelvins; // triple-checking.... - assert(th >= 273.15); // that is theta, not T ! TODO + //assert(th >= 273.15); // that is theta, not T ! TODO assert(rc >= 0); assert(rv >= 0); } @@ -199,7 +199,7 @@ namespace libcloudphxx &rr = std::get<5>(tup); // double-checking.... - assert(th >= 273.15); // TODO: that's theta, not T! + //assert(th >= 273.15); // TODO: that's theta, not T! assert(rc >= 0); assert(rv >= 0); assert(rr >= 0); @@ -281,7 +281,7 @@ namespace libcloudphxx // hopefully true for RK4 assert(F.r == rv); // triple-checking.... - assert(th >= 273.15); // that is theta, not T ! TODO + //assert(th >= 273.15); // that is theta, not T ! TODO assert(rc >= 0); assert(rv >= 0); assert(rr >= 0); diff --git a/include/libcloudph++/blk_1m/formulae.hpp b/include/libcloudph++/blk_1m/formulae.hpp index c1b6d54c..13c6da9f 100644 --- a/include/libcloudph++/blk_1m/formulae.hpp +++ b/include/libcloudph++/blk_1m/formulae.hpp @@ -99,22 +99,16 @@ namespace libcloudphxx const quantity &rv, const quantity &rvs, const quantity &rvsi, - const quantity &th, - const quantity &p, + const quantity &th, + const quantity &p, const quantity &dt ) { - const quantity T = std::pow((p/real_t(100000)/si::pascal), 0.286); - const quantity taunuc = dt; // time scale for nucleation - //homogeneous part - if (T < real_t(233.16) * si::kelvins) { // only if T < -40 C - real_t beta = (T > real_t(213.16)*si::kelvins) ? real_t(0.1)+real_t(0.9)*(T-real_t(213.16)*si::kelvins)/(real_t(20)*si::kelvins) : real_t(0.1); + const quantity T = th * std::pow((p/real_t(100000)/si::pascal), 0.286); + const quantity taunuc = dt; // timescale for nucleation + real_t beta = (T/si::kelvins > real_t(213.16)*si::dimensionless()) ? real_t(0.1)+real_t(0.9)*(T-real_t(213.16)*si::kelvins)/(real_t(20)*si::kelvins) : real_t(0.1); quantity rv_adj = beta * rvs + (1-beta) * rvsi; - //nucleation rate: - quantity::type, real_t> hom1 = (real_t(1)-std::exp(-dt/taunuc)) * std::max(real_t(0), rv-rv_adj) / dt; - } - else { - quantity::type, real_t> hom1 = real_t(0); - } + //homogeneous nucleation rate (only if T < -40 C): + quantity::type, real_t> hom1 = (T/si::kelvins < real_t(233.16)*si::dimensionless()) ? (real_t(1)-std::exp(-dt/taunuc)) * std::max(real_t(0)*si::dimensionless(), rv-rv_adj) / dt : real_t(0)/si::seconds; return hom1; } @@ -122,57 +116,51 @@ namespace libcloudphxx template quantity::type, real_t> hom_A_nucleation_2( const quantity &rc, - const quantity &th, - const quantity &p, + const quantity &th, + const quantity &p, const quantity &dt ) { - const quantity T = std::pow((p/real_t(100000)/si::pascal), 0.286); + const quantity T = th * std::pow((p/real_t(100000)/si::pascal), 0.286); const quantity taunuc = dt; // time scale for nucleation - //homogeneous part - if (T < real_t(233.16) * si::kelvins) { // only if T < -40 C - //nucleation rate: - quantity::type, real_t> hom2 = (real_t(1)-std::exp(-dt/taunuc)) * rc / dt; - } - else { - quantity::type, real_t> hom2 = real_t(0); - } + //homogeneous nucleation rate (only if T < -40 C) + quantity::type, real_t> hom2 = (T/si::kelvins < real_t(233.16)*si::dimensionless()) ? (real_t(1)-std::exp(-dt/taunuc)) * rc / dt : real_t(0)/si::seconds; return hom2; } //ice A heterogeneous nucleation template quantity::type, real_t> het_A_nucleation( - const quantity &ra, + const quantity &ria, const quantity &rc, - const quantity &th, - const quantity &p, + const quantity &th, + const quantity &p, const quantity &rhod_0, const quantity &dt ) { - const quantity T = std::pow((p/real_t(100000)/si::pascal), 0.286); - const quantity taunuc = dt; // time scale for nucleation + const quantity T = th * std::pow((p/real_t(100000)/si::pascal), 0.286); + const quantity taunuc = dt; // timescale for nucleation //heterogeneous part //calculating the mass of small ice A particle: - real_t IWCS = std::min(real_t(1e-3), rhod_0*ra/si::mass_density, real_t(2.52e-4)*std::pow(real_t(1e3)*rhod_0*ra/si::mass_density,real_t(0.837))); - real_t rho_i = real_t(916.8); //ice density - real_t alpha = std::max(real_t(2.2e4), real_t(-4.99e3)-real_t(4.94e4)*std::log10(real_t(1e3)*IWCS)); - quantity m_as = real_t(2*std::numbers::pi)*rho_i/std::pow(alpha,3) *si::mass; + quantity IWCS = std::min(std::min(real_t(1e-3)*si::dimensionless(), rhod_0*ria/(si::kilogram/si::cubic_metres)), real_t(2.52e-4)*std::pow(real_t(1e3)*rhod_0*ria/(si::kilogram/si::cubic_metres),real_t(0.837)) *si::dimensionless()) * si::kilograms / si::cubic_meters; + quantity rho_i = real_t(916.8) *si::kilograms / si::cubic_meters; //ice density + quantity::type, real_t> alpha = std::max(real_t(2.2e4), real_t(-4.99e3)-real_t(4.94e4)*std::log10(real_t(1e3)*IWCS/si::kilograms*si::cubic_meters)) / si::meters; + quantity m_as = real_t(2)*pi()*rho_i/(std::pow(alpha*si::meters,3)/si::cubic_meters); //calculating the mass of large ice A particle: - real_t IWCL = rhod_0*ra - IWCS; + quantity IWCL = rhod_0*ria - IWCS; real_t a_mu = real_t(5.2) + real_t(1.3e-3)*(T/si::kelvin-273.16); real_t b_mu = real_t(0.026) - real_t(1.2e-3)*(T/si::kelvin-273.16); real_t a_sigma = real_t(0.47) + real_t(2.1e-3)*(T/si::kelvin-273.16); real_t b_sigma = real_t(0.018) - real_t(2.1e-4)*(T/si::kelvin-273.16); - real_t mu = a_mu + b_mu*std::log10(real_t(1e3)*IWCL); - real_t sigma = a_sigma + b_sigma*std::log10(real_t(1e3)*IWCL); - quantity m_al = real_t(1.67e17*pi())*rho_i*std::exp(3*mu +9/2*std::pow(sigma,2))*si::mass; + real_t mu = a_mu + b_mu*std::log10(real_t(1e3)*IWCL/si::kilograms*si::cubic_meters); + real_t sigma = a_sigma + b_sigma*std::log10(real_t(1e3)*IWCL/si::kilograms*si::cubic_meters); + quantity m_al = real_t(1.67e17)*pi()*rho_i*std::exp(3*mu +9/2*std::pow(sigma,2))*si::cubic_meters; //mass of average ice A particle - real_t delta = IWCS/(ra*rhod_0/si::mass_density); + real_t delta = IWCS/(ria*rhod_0); quantity ma = delta*m_as + (real_t(1)-delta)*m_al; //concentration of ice nuclei: - quantity::type, real_t> N_in = std::min(real_t(1e5), real_t(1e-2)*std::exp(real_t(0.6)*(real_t(273.16)-T/si::kelvins))); + quantity::type, real_t> N_in = std::min(real_t(1e5), real_t(1e-2)*std::exp(real_t(0.6)*(real_t(273.16)-T/si::kelvins))) / si::cubic_meters; //nucleation rate: - quantity::type, real_t> het = (real_t(1)-std::exp(-dt/taunuc)) * std::min(rc, std::max(real_t(0), N_in*std::max(real_t(1e-12)*si::mass, ma)/rhod_0 - ra)) / dt; + quantity::type, real_t> het = (real_t(1)-std::exp(-dt/taunuc)) * std::min(rc, std::max(real_t(0)*si::dimensionless(), N_in*std::max(real_t(1e-12)*si::dimensionless(), ma/si::kilograms)/rhod_0*si::kilograms - ria)) / dt; return het; } diff --git a/include/libcloudph++/blk_1m/options.hpp b/include/libcloudph++/blk_1m/options.hpp index 8e703680..38526171 100644 --- a/include/libcloudph++/blk_1m/options.hpp +++ b/include/libcloudph++/blk_1m/options.hpp @@ -20,7 +20,11 @@ namespace libcloudphxx revp = true, // evaporation of rain conv = true, // autoconversion accr = true, // accretion - sedi = true; // sedimentation + sedi = true, // sedimentation + ice = true, // enable ice processes + homA1 = true, // homogeneous nucleation 1 of ice A + homA2 = true, // homogeneous nucleation 2 of ice A + hetA = true; // heterogeneous nucleation of ice A real_t r_c0 = 5e-4, // autoconv. threshold k_acnv = 0.001, // Kessler autoconversion (eq. 5a in Grabowski & Smolarkiewicz 1996) diff --git a/include/libcloudph++/blk_1m/rhs_cellwise.hpp b/include/libcloudph++/blk_1m/rhs_cellwise.hpp index e3c54a45..1daf00e2 100644 --- a/include/libcloudph++/blk_1m/rhs_cellwise.hpp +++ b/include/libcloudph++/blk_1m/rhs_cellwise.hpp @@ -1,6 +1,7 @@ /** @file * @copyright University of Warsaw - * @brief Autoconversion and collection righ-hand side terms using Kessler formulae + * @brief Autoconversion and collection righ-hand side terms using Kessler formulae. + * Ice nucleation, growth by deposition and riming, and melting from Grabowski (1999). * @section LICENSE * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/) */ @@ -13,94 +14,149 @@ namespace libcloudphxx { namespace blk_1m { -//+ // template void rhs_cellwise( const opts_t &opts, cont_t &dot_rc_cont, cont_t &dot_rr_cont, - cont_t &dot_ra_cont, const cont_t &rc_cont, - const cont_t &rr_cont, - const cont_t &ra_cont, - const cont_t &th_cont, - const cont_t &p_cont, - const cont_t &rhod_cont, - const real_t &dt + const cont_t &rr_cont ) // { - for (auto tup : zip(dot_rc_cont, dot_rr_cont, dot_ra_cont, rc_cont, rr_cont, ra_cont, th_cont, p_cont, rhod_cont)) + for (auto tup : zip(dot_rc_cont, dot_rr_cont, rc_cont, rr_cont)) { real_t rc_to_rr = 0, - rc_to_ra = 0, &dot_rc = std::get<0>(tup), - &dot_rr = std::get<1>(tup), - &dot_ra = std::get<2>(tup); + &dot_rr = std::get<1>(tup); const real_t - &rc = std::get<3>(tup), - &rr = std::get<4>(tup), - &ra = std::get<5>(tup), - &th = std::get<6>(tup), - &p = std::get<7>(tup), - &rhod = std::get<8>(tup); + &rc = std::get<2>(tup), + &rr = std::get<3>(tup); // autoconversion if (opts.conv) { - rc_to_rr += ( - formulae::autoconversion_rate( - rc * si::dimensionless(), - opts.r_c0 * si::dimensionless(), - opts.k_acnv / si::seconds - ) * si::seconds // to make it dimensionless - ); + rc_to_rr += ( + formulae::autoconversion_rate( + rc * si::dimensionless(), + opts.r_c0 * si::dimensionless(), + opts.k_acnv / si::seconds + ) * si::seconds // to make it dimensionless + ); } // collection if (opts.accr) { - rc_to_rr += ( - formulae::collection_rate( - rc * si::dimensionless(), - rr * si::dimensionless() - ) * si::seconds // to make it dimensionless - ); + rc_to_rr += ( + formulae::collection_rate( + rc * si::dimensionless(), + rr * si::dimensionless() + ) * si::seconds // to make it dimensionless + ); } + dot_rr += rc_to_rr; + dot_rc -= rc_to_rr; + } + } + +//+ template + void rhs_cellwise_ice( + const opts_t &opts, + cont_t &dot_rc_cont, + cont_t &dot_rr_cont, + cont_t &dot_rv_cont, + cont_t &dot_ria_cont, + const cont_t &rc_cont, + const cont_t &rr_cont, + const cont_t &rv_cont, + const cont_t &ria_cont, + const cont_t &th_cont, + const cont_t &p_cont, + const cont_t &rhod_cont, + const real_t &dt + ) +// + { + for (auto tup : zip(dot_rc_cont, dot_rr_cont, dot_rv_cont, dot_ria_cont, rc_cont, rr_cont, rv_cont, ria_cont, th_cont, p_cont, rhod_cont)) + { + using namespace common; + real_t + rc_to_ria = 0, + rv_to_ria = 0, + &dot_rc = std::get<0>(tup), + &dot_rr = std::get<1>(tup), + &dot_rv = std::get<2>(tup), + &dot_ria = std::get<3>(tup); + const real_t + &rc = std::get<4>(tup), + &rr = std::get<5>(tup), + &rv = std::get<6>(tup), + &ria = std::get<7>(tup); + + const quantity + th = std::get<8>(tup) * si::kelvins; + + const quantity + p = std::get<9>(tup) * si::pascals; + + const quantity + rhod = std::get<10>(tup) * si::kilograms / si::cubic_metres; + + quantity T = th * theta_std::exner(p); + quantity rvs = const_cp::r_vs(T, p); + quantity rvsi = const_cp::r_vsi(T, p); + // ice A heterogeneous nucleation - //if (opts.hetA) - //{ - rc_to_ra += ( + if (opts.hetA) + { + rc_to_ria += ( formulae::het_A_nucleation( - ra * si::dimensionless(), + ria * si::dimensionless(), rc * si::dimensionless(), - th * si::kelvin, - p * si::pascal, - rhod * si::mass_density, - dt * si::time + th, + p, + rhod, + dt * si::seconds ) * si::seconds // to make it dimensionless ); - //} + } + + // ice A homogeneous nucleation 1 + if (opts.homA1) + { + rv_to_ria += ( + formulae::hom_A_nucleation_1( + rv * si::dimensionless(), + rvs * si::dimensionless(), + rvsi * si::dimensionless(), + th, + p, + dt * si::seconds + ) * si::seconds // to make it dimensionless + ); + } // ice A homogeneous nucleation 2 - //if (opts.homA2) - //{ - rc_to_ra += ( + if (opts.homA2) + { + rc_to_ria += ( formulae::hom_A_nucleation_2( rc * si::dimensionless(), - th * si::kelvin, - p * si::pascal, - dt * si::time + th, + p, + dt * si::seconds ) * si::seconds // to make it dimensionless ); - //} + } - dot_rr += rc_to_rr; - dot_rc -= rc_to_rr; - dot_rc -= rc_to_ra; - dot_ra += rc_to_ra; + dot_rc -= rc_to_ria; + dot_rv -= rv_to_ria; + dot_ria += rc_to_ria + rv_to_ria; } } @@ -111,18 +167,21 @@ namespace libcloudphxx cont_t &dot_rv_cont, cont_t &dot_rc_cont, cont_t &dot_rr_cont, - cont_t &dot_ra_cont, + cont_t &dot_ria_cont, const cont_t &rhod_cont, const cont_t &p_cont, const cont_t &th_cont, const cont_t &rv_cont, const cont_t &rc_cont, const cont_t &rr_cont, - const cont_t &ra_cont, + const cont_t &ria_cont, const real_t &dt // time step in seconds ) { - rhs_cellwise(opts, dot_rc_cont, dot_rr_cont, dot_ra_cont, rc_cont, rr_cont, ra_cont, th_cont, p_cont, rhod_cont, dt); //calculate T from theta + rhs_cellwise(opts, dot_rc_cont, dot_rr_cont, rc_cont, rr_cont); + if (opts.ice) { + rhs_cellwise_ice(opts, dot_rc_cont, dot_rr_cont, dot_rv_cont, dot_ria_cont, rc_cont, rr_cont, rv_cont, ria_cont, th_cont, p_cont, rhod_cont, dt); + } // rain evaporation treated as a force in Newthon-Raphson saturation adjustment for (auto tup : zip(dot_th_cont, dot_rv_cont, dot_rr_cont, rhod_cont, p_cont, th_cont, rv_cont, rr_cont)) diff --git a/include/libcloudph++/common/const_cp.hpp b/include/libcloudph++/common/const_cp.hpp index 1029ba3d..4fb387b0 100644 --- a/include/libcloudph++/common/const_cp.hpp +++ b/include/libcloudph++/common/const_cp.hpp @@ -11,6 +11,7 @@ namespace libcloudphxx namespace const_cp { using moist_air::c_pw; + using moist_air::c_pi; using moist_air::c_pv; using moist_air::R_v; using moist_air::eps; @@ -20,10 +21,11 @@ namespace libcloudphxx // water triple point parameters libcloudphxx_const(si::pressure, p_tri, 611.73, si::pascals) // pressure libcloudphxx_const(si::temperature, T_tri, 273.16, si::kelvins) // temperature - libcloudphxx_const(energy_over_mass, l_tri, 2.5e6, si::joules / si::kilograms) // latent heat of evaporation + libcloudphxx_const(energy_over_mass, l_tri_evap, 2.5e6, si::joules / si::kilograms) // latent heat of evaporation + libcloudphxx_const(energy_over_mass, l_tri_sublim, 2.83e6, si::joules / si::kilograms) // latent heat of sublimation - // saturation vapour pressure for water assuming constant c_p_v and c_p_w - // with constants taken at triple point + // saturation vapour pressure with respect to liquid water + // assuming constant c_p_v and c_p_w with constants taken at triple point // (solution to the Clausius-Clapeyron equation assuming rho_vapour << rho_liquid) // template @@ -34,12 +36,26 @@ namespace libcloudphxx // { return p_tri() * exp( - (l_tri() + (c_pw() - c_pv()) * T_tri()) / R_v() * (real_t(1) / T_tri() - real_t(1) / T) + (l_tri_evap() + (c_pw() - c_pv()) * T_tri()) / R_v() * (real_t(1) / T_tri() - real_t(1) / T) - (c_pw() - c_pv()) / R_v() * std::log(T / T_tri()) ); } - // saturation vapour mixing ratio for water as a function of pressure and temperature + // saturation vapour pressure with respect to ice assuming constant c_p_v and c_p_i + // with constants taken at triple point + template + BOOST_GPU_ENABLED + quantity p_vsi( + const quantity &T + ) + { + return p_tri() * exp( + (l_tri_sublim() + (c_pi() - c_pv()) * T_tri()) / R_v() * (real_t(1) / T_tri() - real_t(1) / T) + - (c_pi() - c_pv()) / R_v() * std::log(T / T_tri()) + ); + } + + // saturation vapour mixing ratio with respect to liquid water as a function of pressure and temperature template BOOST_GPU_ENABLED quantity r_vs( @@ -49,13 +65,23 @@ namespace libcloudphxx return eps() / (p / p_vs(T) - 1); } - // latent heat for constant c_p + // saturation vapour mixing ratio with respect to ice as a function of pressure and temperature + template + BOOST_GPU_ENABLED + quantity r_vsi( + const quantity &T, + const quantity &p + ) { + return eps() / (p / p_vsi(T) - 1); + } + + // latent heat of evaporation for constant c_p template BOOST_GPU_ENABLED quantity::type , real_t> l_v( const quantity &T ) { - return l_tri() + (c_pv() - c_pw()) * (T - T_tri()); + return l_tri_evap() + (c_pv() - c_pw()) * (T - T_tri()); } }; }; diff --git a/include/libcloudph++/common/moist_air.hpp b/include/libcloudph++/common/moist_air.hpp index 536f3361..e78dcea1 100644 --- a/include/libcloudph++/common/moist_air.hpp +++ b/include/libcloudph++/common/moist_air.hpp @@ -26,6 +26,7 @@ namespace libcloudphxx libcloudphxx_const(energy_over_temperature_mass, c_pd, 1005, si::joules / si::kilograms / si::kelvins) // dry air libcloudphxx_const(energy_over_temperature_mass, c_pv, 1850, si::joules / si::kilograms / si::kelvins) // water vapour libcloudphxx_const(energy_over_temperature_mass, c_pw, 4218, si::joules / si::kilograms / si::kelvins) // liquid water + libcloudphxx_const(energy_over_temperature_mass, c_pi, 2107, si::joules / si::kilograms / si::kelvins) // ice // molar masses libcloudphxx_const(mass_over_amount, M_d, 0.02897, si::kilograms / si::moles) // dry air (Curry & Webster / Seinfeld & Pandis) diff --git a/tests/python/unit/api_blk_1m.py b/tests/python/unit/api_blk_1m.py index a084397d..6e87a18d 100644 --- a/tests/python/unit/api_blk_1m.py +++ b/tests/python/unit/api_blk_1m.py @@ -11,16 +11,21 @@ print("revp =", opts.revp) print("conv =", opts.conv) print("accr =", opts.accr) -print("sedi =", opts.sedi) +print("sedi =", opts.sedi) +print("ice =", opts.ice) +print("homA1 =", opts.homA1) +print("homA2 =", opts.homA2) +print("hetA =", opts.hetA) print("r_c0 =", opts.r_c0) print("r_eps =", opts.r_eps) rhod = arr_t([1. ]) p = arr_t([1.e5]) th = arr_t([300.]) -rv = arr_t([0. ]) rc = arr_t([0.01]) +rv = arr_t([0. ]) rr = arr_t([0. ]) +ria = arr_t([0. ]) dt = 1 dz = 1 @@ -58,6 +63,7 @@ dot_rc = arr_t([0.]) dot_rr = arr_t([0.]) +dot_ria = arr_t([0.]) blk_1m.rhs_cellwise(opts, dot_rc, dot_rr, rc, rr) assert dot_rc != 0 # some water should have coalesced assert dot_rr != 0 @@ -69,7 +75,7 @@ dot_rc = arr_t([0.]) dot_rr = arr_t([0.]) -blk_1m.rhs_cellwise_nwtrph(opts, dot_th, dot_rv, dot_rc, dot_rr, rhod, p, th, rv, rc, rr, dt) +blk_1m.rhs_cellwise_nwtrph(opts, dot_th, dot_rv, dot_rc, dot_rr, dot_ria, rhod, p, th, rv, rc, rr, ria, dt) assert dot_rc != 0 # some water should have coalesced assert dot_rr != 0 assert dot_th != 0 # some rain should have evaporated @@ -80,3 +86,12 @@ flux = blk_1m.rhs_columnwise(opts, dot_rr, rhod, rr, dz) assert flux == 0 assert dot_rr == dot_rr_old # no rain water -> no precip + +th = arr_t([230.]) #testing ice physics +dot_rc = arr_t([0.]) +dot_rr = arr_t([0.]) +dot_rv = arr_t([0.]) +dot_ria = arr_t([0.]) +rv = arr_t([0.01]) +blk_1m.rhs_cellwise_ice(opts, dot_rc, dot_rr, dot_rv, dot_ria, rc, rr, rv, ria, th, p, rhod, dt) +assert dot_ria != 0 From 9ae4f9636d580696e1e81fc100894418916ed500 Mon Sep 17 00:00:00 2001 From: Agnieszka Makulska Date: Tue, 29 Oct 2024 12:00:43 +0100 Subject: [PATCH 05/17] rhs_cellwise_nwtrhp_ice function --- bindings/python/blk_1m.hpp | 79 ++++--- bindings/python/lib.cpp | 5 +- include/libcloudph++/blk_1m/options.hpp | 1 - include/libcloudph++/blk_1m/rhs_cellwise.hpp | 222 ++++++++++--------- include/libcloudph++/common/const_cp.hpp | 21 +- tests/python/unit/api_blk_1m.py | 9 +- 6 files changed, 184 insertions(+), 153 deletions(-) diff --git a/bindings/python/blk_1m.hpp b/bindings/python/blk_1m.hpp index 1d918900..db1c227f 100644 --- a/bindings/python/blk_1m.hpp +++ b/bindings/python/blk_1m.hpp @@ -121,48 +121,47 @@ namespace libcloudphxx np2bz(rr) ); } - - template - void rhs_cellwise_ice( - const b1m::opts_t &opts, - bp_array &dot_rc, - bp_array &dot_rr, - bp_array &dot_rv, - bp_array &dot_ria, - const bp_array &rc, - const bp_array &rr, - const bp_array &rv, - const bp_array &ria, - const bp_array &theta, - const bp_array &p, - const bp_array &rhod, - const typename arr_t::T_numtype &dt - ) + + template + void rhs_cellwise_nwtrph( + const b1m::opts_t &opts, + bp_array &dot_th, + bp_array &dot_rv, + bp_array &dot_rc, + bp_array &dot_rr, + const bp_array &rhod, + const bp_array &p, + const bp_array &th, + const bp_array &rv, + const bp_array &rc, + const bp_array &rr, + const typename arr_t::T_numtype &dt + ) { - arr_t - np2bz_dot_rc(np2bz(dot_rc)), - np2bz_dot_rr(np2bz(dot_rr)), - np2bz_dot_rv(np2bz(dot_rv)), - np2bz_dot_ria(np2bz(dot_ria)); - b1m::rhs_cellwise_ice( - opts, - np2bz_dot_rc, - np2bz_dot_rr, - np2bz_dot_rv, - np2bz_dot_ria, - np2bz(rc), - np2bz(rr), - np2bz(rv), - np2bz(ria), - np2bz(theta), - np2bz(p), - np2bz(rhod), - dt - ); - } + arr_t + np2bz_dot_rc(np2bz(dot_rc)), + np2bz_dot_rr(np2bz(dot_rr)), + np2bz_dot_rv(np2bz(dot_rv)), + np2bz_dot_th(np2bz(dot_th)); + b1m::rhs_cellwise_nwtrph( + opts, + np2bz_dot_th, + np2bz_dot_rv, + np2bz_dot_rc, + np2bz_dot_rr, + np2bz(rhod), + np2bz(p), + np2bz(th), + np2bz(rv), + np2bz(rc), + np2bz(rr), + dt + ); + } + template - void rhs_cellwise_nwtrph( + void rhs_cellwise_nwtrph_ice( const b1m::opts_t &opts, bp_array &dot_th, bp_array &dot_rv, @@ -185,7 +184,7 @@ namespace libcloudphxx np2bz_dot_rr(np2bz(dot_rr)), np2bz_dot_rv(np2bz(dot_rv)), np2bz_dot_th(np2bz(dot_th)); - b1m::rhs_cellwise_nwtrph( + b1m::rhs_cellwise_nwtrph_ice( opts, np2bz_dot_th, np2bz_dot_rv, diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index d7adc5f4..dfe3b84c 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -155,7 +155,6 @@ BOOST_PYTHON_MODULE(libcloudphxx) .def_readwrite("conv", &b1m::opts_t::conv) .def_readwrite("accr", &b1m::opts_t::accr) .def_readwrite("sedi", &b1m::opts_t::sedi) - .def_readwrite("ice", &b1m::opts_t::ice) .def_readwrite("homA1", &b1m::opts_t::homA1) .def_readwrite("homA2", &b1m::opts_t::homA2) .def_readwrite("hetA", &b1m::opts_t::hetA) @@ -167,8 +166,8 @@ BOOST_PYTHON_MODULE(libcloudphxx) bp::def("adj_cellwise_constp", blk_1m::adj_cellwise_constp); bp::def("adj_cellwise_nwtrph", blk_1m::adj_cellwise_nwtrph); bp::def("rhs_cellwise", blk_1m::rhs_cellwise); - bp::def("rhs_cellwise_ice", blk_1m::rhs_cellwise_ice); - bp::def("rhs_cellwise_nwtrph", blk_1m::rhs_cellwise_nwtrph); + bp::def("rhs_cellwise_nwtrph", blk_1m::rhs_cellwise_nwtrph); + bp::def("rhs_cellwise_nwtrph_ice", blk_1m::rhs_cellwise_nwtrph_ice); bp::def("rhs_columnwise", blk_1m::rhs_columnwise); // TODO: handle the returned flux } diff --git a/include/libcloudph++/blk_1m/options.hpp b/include/libcloudph++/blk_1m/options.hpp index 38526171..792b7e44 100644 --- a/include/libcloudph++/blk_1m/options.hpp +++ b/include/libcloudph++/blk_1m/options.hpp @@ -21,7 +21,6 @@ namespace libcloudphxx conv = true, // autoconversion accr = true, // accretion sedi = true, // sedimentation - ice = true, // enable ice processes homA1 = true, // homogeneous nucleation 1 of ice A homA2 = true, // homogeneous nucleation 2 of ice A hetA = true; // heterogeneous nucleation of ice A diff --git a/include/libcloudph++/blk_1m/rhs_cellwise.hpp b/include/libcloudph++/blk_1m/rhs_cellwise.hpp index 1daf00e2..af5b2fa6 100644 --- a/include/libcloudph++/blk_1m/rhs_cellwise.hpp +++ b/include/libcloudph++/blk_1m/rhs_cellwise.hpp @@ -63,105 +63,72 @@ namespace libcloudphxx } } -// template - void rhs_cellwise_ice( + void rhs_cellwise_nwtrph( const opts_t &opts, + cont_t &dot_th_cont, + cont_t &dot_rv_cont, cont_t &dot_rc_cont, cont_t &dot_rr_cont, - cont_t &dot_rv_cont, - cont_t &dot_ria_cont, + const cont_t &rhod_cont, + const cont_t &p_cont, + const cont_t &th_cont, + const cont_t &rv_cont, const cont_t &rc_cont, const cont_t &rr_cont, - const cont_t &rv_cont, - const cont_t &ria_cont, - const cont_t &th_cont, - const cont_t &p_cont, - const cont_t &rhod_cont, const real_t &dt ) -// { - for (auto tup : zip(dot_rc_cont, dot_rr_cont, dot_rv_cont, dot_ria_cont, rc_cont, rr_cont, rv_cont, ria_cont, th_cont, p_cont, rhod_cont)) + rhs_cellwise(opts, dot_rc_cont, dot_rr_cont, rc_cont, rr_cont); + + // rain evaporation treated as a force in Newthon-Raphson saturation adjustment + for (auto tup : zip(dot_th_cont, dot_rv_cont, dot_rr_cont, rhod_cont, p_cont, th_cont, rv_cont, rr_cont)) { using namespace common; + real_t - rc_to_ria = 0, - rv_to_ria = 0, - &dot_rc = std::get<0>(tup), - &dot_rr = std::get<1>(tup), - &dot_rv = std::get<2>(tup), - &dot_ria = std::get<3>(tup); - const real_t - &rc = std::get<4>(tup), - &rr = std::get<5>(tup), - &rv = std::get<6>(tup), - &ria = std::get<7>(tup); + rr_to_rv = 0, + &dot_th = std::get<0>(tup), + &dot_rv = std::get<1>(tup), + &dot_rr = std::get<2>(tup); - const quantity - th = std::get<8>(tup) * si::kelvins; + const quantity + rhod = std::get<3>(tup) * si::kilograms / si::cubic_metres; const quantity - p = std::get<9>(tup) * si::pascals; - - const quantity - rhod = std::get<10>(tup) * si::kilograms / si::cubic_metres; + p = std::get<4>(tup) * si::pascals; - quantity T = th * theta_std::exner(p); - quantity rvs = const_cp::r_vs(T, p); - quantity rvsi = const_cp::r_vsi(T, p); + const quantity + th = std::get<5>(tup) * si::kelvins; - // ice A heterogeneous nucleation - if (opts.hetA) - { - rc_to_ria += ( - formulae::het_A_nucleation( - ria * si::dimensionless(), - rc * si::dimensionless(), - th, - p, - rhod, - dt * si::seconds - ) * si::seconds // to make it dimensionless - ); - } + const real_t + &rv = std::get<6>(tup), + &rr = std::get<7>(tup); - // ice A homogeneous nucleation 1 - if (opts.homA1) - { - rv_to_ria += ( - formulae::hom_A_nucleation_1( - rv * si::dimensionless(), - rvs * si::dimensionless(), - rvsi * si::dimensionless(), - th, - p, - dt * si::seconds - ) * si::seconds // to make it dimensionless - ); - } + quantity T = th * theta_std::exner(p); + real_t r_vs = const_cp::r_vs(T, p); - // ice A homogeneous nucleation 2 - if (opts.homA2) - { - rc_to_ria += ( - formulae::hom_A_nucleation_2( - rc * si::dimensionless(), - th, - p, - dt * si::seconds - ) * si::seconds // to make it dimensionless + rr_to_rv += ( + formulae::evaporation_rate( + rv * si::dimensionless(), + r_vs * si::dimensionless(), + rr * si::dimensionless(), + rhod, + p + ) * si::seconds * dt ); - } - dot_rc -= rc_to_ria; - dot_rv -= rv_to_ria; - dot_ria += rc_to_ria + rv_to_ria; + // limiting + rr_to_rv = std::min(rr, rr_to_rv) / dt; + + dot_rv += rr_to_rv; + dot_rr -= rr_to_rv; + dot_th -= const_cp::l_v(T) / (moist_air::c_pd() * theta_std::exner(p)) * dot_rv / si::kelvins; } } template - void rhs_cellwise_nwtrph( + void rhs_cellwise_nwtrph_ice( const opts_t &opts, cont_t &dot_th_cont, cont_t &dot_rv_cont, @@ -175,58 +142,109 @@ namespace libcloudphxx const cont_t &rc_cont, const cont_t &rr_cont, const cont_t &ria_cont, - const real_t &dt // time step in seconds + const real_t &dt ) { rhs_cellwise(opts, dot_rc_cont, dot_rr_cont, rc_cont, rr_cont); - if (opts.ice) { - rhs_cellwise_ice(opts, dot_rc_cont, dot_rr_cont, dot_rv_cont, dot_ria_cont, rc_cont, rr_cont, rv_cont, ria_cont, th_cont, p_cont, rhod_cont, dt); - } // rain evaporation treated as a force in Newthon-Raphson saturation adjustment - for (auto tup : zip(dot_th_cont, dot_rv_cont, dot_rr_cont, rhod_cont, p_cont, th_cont, rv_cont, rr_cont)) + for (auto tup : zip(dot_th_cont, dot_rv_cont, dot_rc_cont, dot_rr_cont, dot_ria_cont, rhod_cont, p_cont, th_cont, rv_cont, rc_cont, rr_cont, ria_cont)) { using namespace common; real_t - tmp = 0, - &dot_th = std::get<0>(tup), - &dot_rv = std::get<1>(tup), - &dot_rr = std::get<2>(tup); + rr_to_rv = 0, + rv_to_ria = 0, + rc_to_ria = 0, + &dot_th = std::get<0>(tup), + &dot_rv = std::get<1>(tup), + &dot_rc = std::get<2>(tup), + &dot_rr = std::get<3>(tup), + &dot_ria = std::get<4>(tup); const quantity - rhod = std::get<3>(tup) * si::kilograms / si::cubic_metres; + rhod = std::get<5>(tup) * si::kilograms / si::cubic_metres; const quantity - p = std::get<4>(tup) * si::pascals; + p = std::get<6>(tup) * si::pascals; const quantity - th = std::get<5>(tup) * si::kelvins; + th = std::get<7>(tup) * si::kelvins; const real_t - &rv = std::get<6>(tup), - &rr = std::get<7>(tup); + &rv = std::get<8>(tup), + &rc = std::get<9>(tup), + &rr = std::get<10>(tup), + &ria = std::get<11>(tup); quantity T = th * theta_std::exner(p); - real_t r_vs = const_cp::r_vs(T, p); - - tmp = ( - formulae::evaporation_rate( - rv * si::dimensionless(), - r_vs * si::dimensionless(), - rr * si::dimensionless(), - rhod, - p - ) * si::seconds * dt + real_t rvs = const_cp::r_vs(T, p); + real_t rvsi = const_cp::r_vsi(T, p); + + // rain evaporation + rr_to_rv += ( + formulae::evaporation_rate( + rv * si::dimensionless(), + rvs * si::dimensionless(), + rr * si::dimensionless(), + rhod, + p + ) * si::seconds * dt ); // limiting - tmp = std::min(rr, tmp) / dt; + rr_to_rv = std::min(rr, rr_to_rv) / dt; - dot_rv += tmp; - dot_rr -= tmp; + // ice A heterogeneous nucleation + if (opts.hetA) + { + rc_to_ria += ( + formulae::het_A_nucleation( + ria * si::dimensionless(), + rc * si::dimensionless(), + th, + p, + rhod, + dt * si::seconds + ) * si::seconds // to make it dimensionless + ); + } + + // ice A homogeneous nucleation 1 + if (opts.homA1) + { + rv_to_ria += ( + formulae::hom_A_nucleation_1( + rv * si::dimensionless(), + rvs * si::dimensionless(), + rvsi * si::dimensionless(), + th, + p, + dt * si::seconds + ) * si::seconds // to make it dimensionless + ); + } + + // ice A homogeneous nucleation 2 + if (opts.homA2) + { + rc_to_ria += ( + formulae::hom_A_nucleation_2( + rc * si::dimensionless(), + th, + p, + dt * si::seconds + ) * si::seconds // to make it dimensionless + ); + } - dot_th -= const_cp::l_v(T) / (moist_air::c_pd() * theta_std::exner(p)) * tmp / si::kelvins; + dot_rr -= rr_to_rv; + dot_rc -= rc_to_ria; + dot_rv += rr_to_rv - rv_to_ria; + dot_ria += rc_to_ria + rv_to_ria; + dot_th -= const_cp::l_v(T) / (moist_air::c_pd() * theta_std::exner(p)) * rr_to_rv / si::kelvins; //heat of vaporisation + dot_th += const_cp::l_sublim(T) / (moist_air::c_pd() * theta_std::exner(p)) * rv_to_ria / si::kelvins; //heat of sublimation + dot_th += const_cp::l_freez(T) / (moist_air::c_pd() * theta_std::exner(p)) * rc_to_ria / si::kelvins; //heat of freezing } } } diff --git a/include/libcloudph++/common/const_cp.hpp b/include/libcloudph++/common/const_cp.hpp index 4fb387b0..71151795 100644 --- a/include/libcloudph++/common/const_cp.hpp +++ b/include/libcloudph++/common/const_cp.hpp @@ -23,6 +23,7 @@ namespace libcloudphxx libcloudphxx_const(si::temperature, T_tri, 273.16, si::kelvins) // temperature libcloudphxx_const(energy_over_mass, l_tri_evap, 2.5e6, si::joules / si::kilograms) // latent heat of evaporation libcloudphxx_const(energy_over_mass, l_tri_sublim, 2.83e6, si::joules / si::kilograms) // latent heat of sublimation + libcloudphxx_const(energy_over_mass, l_tri_freez, 3.34e5, si::joules / si::kilograms) // latent heat of freezing/melting // saturation vapour pressure with respect to liquid water // assuming constant c_p_v and c_p_w with constants taken at triple point @@ -51,7 +52,6 @@ namespace libcloudphxx { return p_tri() * exp( (l_tri_sublim() + (c_pi() - c_pv()) * T_tri()) / R_v() * (real_t(1) / T_tri() - real_t(1) / T) - - (c_pi() - c_pv()) / R_v() * std::log(T / T_tri()) ); } @@ -83,6 +83,25 @@ namespace libcloudphxx ) { return l_tri_evap() + (c_pv() - c_pw()) * (T - T_tri()); } + + // latent heat of sublimation for constant c_p + template + BOOST_GPU_ENABLED + quantity::type , real_t> l_sublim( + const quantity &T + ) { + return l_tri_sublim() + (c_pv() - c_pi()) * (T - T_tri()); + } + + // latent heat of freezing/melting for constant c_p + template + BOOST_GPU_ENABLED + quantity::type , real_t> l_freez( + const quantity &T + ) { + return l_tri_freez() + (c_pw() - c_pi()) * (T - T_tri()); + } + }; }; }; diff --git a/tests/python/unit/api_blk_1m.py b/tests/python/unit/api_blk_1m.py index 6e87a18d..1c36b53a 100644 --- a/tests/python/unit/api_blk_1m.py +++ b/tests/python/unit/api_blk_1m.py @@ -12,7 +12,6 @@ print("conv =", opts.conv) print("accr =", opts.accr) print("sedi =", opts.sedi) -print("ice =", opts.ice) print("homA1 =", opts.homA1) print("homA2 =", opts.homA2) print("hetA =", opts.hetA) @@ -63,7 +62,6 @@ dot_rc = arr_t([0.]) dot_rr = arr_t([0.]) -dot_ria = arr_t([0.]) blk_1m.rhs_cellwise(opts, dot_rc, dot_rr, rc, rr) assert dot_rc != 0 # some water should have coalesced assert dot_rr != 0 @@ -75,7 +73,7 @@ dot_rc = arr_t([0.]) dot_rr = arr_t([0.]) -blk_1m.rhs_cellwise_nwtrph(opts, dot_th, dot_rv, dot_rc, dot_rr, dot_ria, rhod, p, th, rv, rc, rr, ria, dt) +blk_1m.rhs_cellwise_nwtrph(opts, dot_th, dot_rv, dot_rc, dot_rr, rhod, p, th, rv, rc, rr, dt) assert dot_rc != 0 # some water should have coalesced assert dot_rr != 0 assert dot_th != 0 # some rain should have evaporated @@ -92,6 +90,5 @@ dot_rr = arr_t([0.]) dot_rv = arr_t([0.]) dot_ria = arr_t([0.]) -rv = arr_t([0.01]) -blk_1m.rhs_cellwise_ice(opts, dot_rc, dot_rr, dot_rv, dot_ria, rc, rr, rv, ria, th, p, rhod, dt) -assert dot_ria != 0 +blk_1m.rhs_cellwise_nwtrph_ice(opts, dot_th, dot_rv, dot_rc, dot_rr, dot_ria, rhod, p, th, rv, rc, rr, ria, dt) +assert dot_ria != 0 \ No newline at end of file From 7e1595cfd4298b8a2655d8cb443f551c47644195 Mon Sep 17 00:00:00 2001 From: Agnieszka Makulska Date: Thu, 31 Oct 2024 12:19:44 +0100 Subject: [PATCH 06/17] fixes and adding hetB nucleation in formulae --- bindings/python/blk_1m.hpp | 372 +++++++++---------- bindings/python/lib.cpp | 1 + cmake/FindThrust.cmake | 66 ++++ cmake/LICENSE_VTKm.txt | 78 ++++ include/libcloudph++/blk_1m/formulae.hpp | 182 ++++++--- include/libcloudph++/blk_1m/options.hpp | 7 +- include/libcloudph++/blk_1m/rhs_cellwise.hpp | 38 +- include/libcloudph++/common/const_cp.hpp | 23 +- include/libcloudph++/common/moist_air.hpp | 4 +- tests/python/unit/api_blk_1m.py | 3 + 10 files changed, 497 insertions(+), 277 deletions(-) create mode 100644 cmake/FindThrust.cmake create mode 100644 cmake/LICENSE_VTKm.txt diff --git a/bindings/python/blk_1m.hpp b/bindings/python/blk_1m.hpp index db1c227f..00ce7de1 100644 --- a/bindings/python/blk_1m.hpp +++ b/bindings/python/blk_1m.hpp @@ -13,213 +13,209 @@ #include #include -namespace libcloudphxx -{ - namespace python - { +namespace libcloudphxx { + namespace python { namespace b1m = libcloudphxx::blk_1m; - namespace blk_1m - { - template + namespace blk_1m { + template void adj_cellwise( - const b1m::opts_t& opts, - const bp_array &rhod, - bp_array &th, - bp_array &rv, - bp_array &rc, - bp_array &rr, - const typename arr_t::T_numtype &dt - ) - { - arr_t - np2bz_th(np2bz(th)), - np2bz_rv(np2bz(rv)), - np2bz_rc(np2bz(rc)), - np2bz_rr(np2bz(rr)); - b1m::adj_cellwise( - opts, - np2bz(rhod), // since it is const, it may be a temporary object - np2bz_th, - np2bz_rv, - np2bz_rc, - np2bz_rr, - dt - ); + const b1m::opts_t &opts, + const bp_array &rhod, + bp_array &th, + bp_array &rv, + bp_array &rc, + bp_array &rr, + const typename arr_t::T_numtype &dt + ) { + arr_t + np2bz_th(np2bz(th)), + np2bz_rv(np2bz(rv)), + np2bz_rc(np2bz(rc)), + np2bz_rr(np2bz(rr)); + b1m::adj_cellwise( + opts, + np2bz(rhod), // since it is const, it may be a temporary object + np2bz_th, + np2bz_rv, + np2bz_rc, + np2bz_rr, + dt + ); } - template + template void adj_cellwise_constp( - const b1m::opts_t& opts, - const bp_array &rhod, - const bp_array &p, - bp_array &th, - bp_array &rv, - bp_array &rc, - bp_array &rr, - const typename arr_t::T_numtype &dt - ) - { - arr_t - np2bz_th(np2bz(th)), - np2bz_rv(np2bz(rv)), - np2bz_rc(np2bz(rc)), - np2bz_rr(np2bz(rr)); - b1m::adj_cellwise_constp( - opts, - np2bz(rhod), // since it is const, it may be a temporary object - np2bz(p), - np2bz_th, - np2bz_rv, - np2bz_rc, - np2bz_rr, - dt - ); + const b1m::opts_t &opts, + const bp_array &rhod, + const bp_array &p, + bp_array &th, + bp_array &rv, + bp_array &rc, + bp_array &rr, + const typename arr_t::T_numtype &dt + ) { + arr_t + np2bz_th(np2bz(th)), + np2bz_rv(np2bz(rv)), + np2bz_rc(np2bz(rc)), + np2bz_rr(np2bz(rr)); + b1m::adj_cellwise_constp( + opts, + np2bz(rhod), // since it is const, it may be a temporary object + np2bz(p), + np2bz_th, + np2bz_rv, + np2bz_rc, + np2bz_rr, + dt + ); } - - template + + template void adj_cellwise_nwtrph( - const b1m::opts_t& opts, - const bp_array &p, - bp_array &th, - bp_array &rv, - bp_array &rc, - const typename arr_t::T_numtype &dt - ) - { - arr_t - np2bz_th(np2bz(th)), - np2bz_rv(np2bz(rv)), - np2bz_rc(np2bz(rc)); - b1m::adj_cellwise_nwtrph( - opts, - np2bz(p), - np2bz_th, - np2bz_rv, - np2bz_rc, - dt - ); + const b1m::opts_t &opts, + const bp_array &p, + bp_array &th, + bp_array &rv, + bp_array &rc, + const typename arr_t::T_numtype &dt + ) { + arr_t + np2bz_th(np2bz(th)), + np2bz_rv(np2bz(rv)), + np2bz_rc(np2bz(rc)); + b1m::adj_cellwise_nwtrph( + opts, + np2bz(p), + np2bz_th, + np2bz_rv, + np2bz_rc, + dt + ); } - template - void rhs_cellwise( - const b1m::opts_t &opts, - bp_array &dot_rc, - bp_array &dot_rr, - const bp_array &rc, - const bp_array &rr - ) - { - arr_t - np2bz_dot_rc(np2bz(dot_rc)), - np2bz_dot_rr(np2bz(dot_rr)); - b1m::rhs_cellwise( - opts, - np2bz_dot_rc, - np2bz_dot_rr, - np2bz(rc), - np2bz(rr) - ); + template + void rhs_cellwise( + const b1m::opts_t &opts, + bp_array &dot_rc, + bp_array &dot_rr, + const bp_array &rc, + const bp_array &rr + ) { + arr_t + np2bz_dot_rc(np2bz(dot_rc)), + np2bz_dot_rr(np2bz(dot_rr)); + b1m::rhs_cellwise( + opts, + np2bz_dot_rc, + np2bz_dot_rr, + np2bz(rc), + np2bz(rr) + ); } - template - void rhs_cellwise_nwtrph( - const b1m::opts_t &opts, - bp_array &dot_th, - bp_array &dot_rv, - bp_array &dot_rc, - bp_array &dot_rr, - const bp_array &rhod, - const bp_array &p, - const bp_array &th, - const bp_array &rv, - const bp_array &rc, - const bp_array &rr, - const typename arr_t::T_numtype &dt - ) - { - arr_t - np2bz_dot_rc(np2bz(dot_rc)), - np2bz_dot_rr(np2bz(dot_rr)), - np2bz_dot_rv(np2bz(dot_rv)), - np2bz_dot_th(np2bz(dot_th)); - b1m::rhs_cellwise_nwtrph( - opts, - np2bz_dot_th, - np2bz_dot_rv, - np2bz_dot_rc, - np2bz_dot_rr, - np2bz(rhod), - np2bz(p), - np2bz(th), - np2bz(rv), - np2bz(rc), - np2bz(rr), - dt - ); + template + void rhs_cellwise_nwtrph( + const b1m::opts_t &opts, + bp_array &dot_th, + bp_array &dot_rv, + bp_array &dot_rc, + bp_array &dot_rr, + const bp_array &rhod, + const bp_array &p, + const bp_array &th, + const bp_array &rv, + const bp_array &rc, + const bp_array &rr, + const typename arr_t::T_numtype &dt + ) { + arr_t + np2bz_dot_rc(np2bz(dot_rc)), + np2bz_dot_rr(np2bz(dot_rr)), + np2bz_dot_rv(np2bz(dot_rv)), + np2bz_dot_th(np2bz(dot_th)); + b1m::rhs_cellwise_nwtrph( + opts, + np2bz_dot_th, + np2bz_dot_rv, + np2bz_dot_rc, + np2bz_dot_rr, + np2bz(rhod), + np2bz(p), + np2bz(th), + np2bz(rv), + np2bz(rc), + np2bz(rr), + dt + ); } - template + template void rhs_cellwise_nwtrph_ice( - const b1m::opts_t &opts, - bp_array &dot_th, - bp_array &dot_rv, - bp_array &dot_rc, - bp_array &dot_rr, - bp_array &dot_ria, - const bp_array &rhod, - const bp_array &p, - const bp_array &th, - const bp_array &rv, - const bp_array &rc, - const bp_array &rr, - const bp_array &ria, - const typename arr_t::T_numtype &dt - ) - { - arr_t - np2bz_dot_ria(np2bz(dot_ria)), - np2bz_dot_rc(np2bz(dot_rc)), - np2bz_dot_rr(np2bz(dot_rr)), - np2bz_dot_rv(np2bz(dot_rv)), - np2bz_dot_th(np2bz(dot_th)); - b1m::rhs_cellwise_nwtrph_ice( - opts, - np2bz_dot_th, - np2bz_dot_rv, - np2bz_dot_rc, - np2bz_dot_rr, - np2bz_dot_ria, - np2bz(rhod), - np2bz(p), - np2bz(th), - np2bz(rv), - np2bz(rc), - np2bz(rr), - np2bz(ria), - dt - ); - } + const b1m::opts_t &opts, + bp_array &dot_th, + bp_array &dot_rv, + bp_array &dot_rc, + bp_array &dot_rr, + bp_array &dot_ria, + //bp_array &dot_rib, + const bp_array &rhod, + const bp_array &p, + const bp_array &th, + const bp_array &rv, + const bp_array &rc, + const bp_array &rr, + const bp_array &ria, + //const bp_array &rib, + const typename arr_t::T_numtype &dt + ) { + arr_t + np2bz_dot_ria(np2bz(dot_ria)), + //np2bz_dot_rib(np2bz(dot_rib)), + np2bz_dot_rc(np2bz(dot_rc)), + np2bz_dot_rr(np2bz(dot_rr)), + np2bz_dot_rv(np2bz(dot_rv)), + np2bz_dot_th(np2bz(dot_th)); + b1m::rhs_cellwise_nwtrph_ice( + opts, + np2bz_dot_th, + np2bz_dot_rv, + np2bz_dot_rc, + np2bz_dot_rr, + np2bz_dot_ria, + //np2bz_dot_rib, + np2bz(rhod), + np2bz(p), + np2bz(th), + np2bz(rv), + np2bz(rc), + np2bz(rr), + np2bz(ria), + //np2bz(rib), + dt + ); + } - template + template typename arr_t::T_numtype rhs_columnwise( - const b1m::opts_t &opts, - bp_array &dot_rr, - const bp_array &rhod, - const bp_array &rr, - const typename arr_t::T_numtype &dz + const b1m::opts_t &opts, + bp_array &dot_rr, + const bp_array &rhod, + const bp_array &rr, + const typename arr_t::T_numtype &dz ) { - arr_t - np2bz_dot_rr(np2bz(dot_rr)); - return b1m::rhs_columnwise( - opts, - np2bz_dot_rr, - np2bz(rhod), - np2bz(rr), - dz - ); - } + arr_t + np2bz_dot_rr(np2bz(dot_rr)); + return b1m::rhs_columnwise( + opts, + np2bz_dot_rr, + np2bz(rhod), + np2bz(rr), + dz + ); + } }; }; }; diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index dfe3b84c..76bde5ac 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -158,6 +158,7 @@ BOOST_PYTHON_MODULE(libcloudphxx) .def_readwrite("homA1", &b1m::opts_t::homA1) .def_readwrite("homA2", &b1m::opts_t::homA2) .def_readwrite("hetA", &b1m::opts_t::hetA) + .def_readwrite("hetB", &b1m::opts_t::hetB) .def_readwrite("r_c0", &b1m::opts_t::r_c0) .def_readwrite("r_eps", &b1m::opts_t::r_eps) .def_readwrite("nwtrph_iters", &b1m::opts_t::nwtrph_iters) diff --git a/cmake/FindThrust.cmake b/cmake/FindThrust.cmake new file mode 100644 index 00000000..38e8fc47 --- /dev/null +++ b/cmake/FindThrust.cmake @@ -0,0 +1,66 @@ +##============================================================================= +## +## Copyright (c) Kitware, Inc. +## All rights reserved. +## See LICENSE_VTKm.txt for details. +## +## This software is distributed WITHOUT ANY WARRANTY; without even +## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +## PURPOSE. See the above copyright notice for more information. +## +## Copyright 2012 Sandia Corporation. +## Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +## the U.S. Government retains certain rights in this software. +## +##============================================================================= + +# +# FindThrust +# +# This module finds the Thrust header files and extrats their version. It +# sets the following variables. +# +# THRUST_INCLUDE_DIR - Include directory for thrust header files. (All header +# files will actually be in the thrust subdirectory.) +# THRUST_VERSION - Version of thrust in the form "major.minor.patch". +# + +find_path( THRUST_INCLUDE_DIR + HINTS + ${THRUST_INCLUDE_DIR} + ${CUDA_INCLUDE_DIRS} + /usr/include/cuda + /usr/local/include + /usr/local/cuda/include + NAMES thrust/version.h + DOC "Thrust headers" + ) +if( THRUST_INCLUDE_DIR ) + list( REMOVE_DUPLICATES THRUST_INCLUDE_DIR ) + include_directories(BEFORE ${THRUST_INCLUDE_DIR} ) +endif() + +# Find thrust version +file( STRINGS ${THRUST_INCLUDE_DIR}/thrust/version.h + version + REGEX "#define THRUST_VERSION[ \t]+([0-9x]+)" + ) +string( REGEX REPLACE + "#define THRUST_VERSION[ \t]+" + "" + version + "${version}" + ) + +string( REGEX MATCH "^[0-9]" major ${version} ) +string( REGEX REPLACE "^${major}00" "" version "${version}" ) +string( REGEX MATCH "^[0-9]" minor ${version} ) +string( REGEX REPLACE "^${minor}0" "" version "${version}" ) +set( THRUST_VERSION "${major}.${minor}.${version}") + +# Check for required components +include( FindPackageHandleStandardArgs ) +find_package_handle_standard_args( Thrust + REQUIRED_VARS THRUST_INCLUDE_DIR + VERSION_VAR THRUST_VERSION + ) diff --git a/cmake/LICENSE_VTKm.txt b/cmake/LICENSE_VTKm.txt new file mode 100644 index 00000000..9750502b --- /dev/null +++ b/cmake/LICENSE_VTKm.txt @@ -0,0 +1,78 @@ +VTKm License Version 1.0 +======================================================================== + +Copyright (c) 2014, +Sandia Corporation, Los Alamos National Security, LLC., +UT-Battelle, LLC., Kitware Inc., University of California Davis +All rights reserved. + +Sandia National Laboratories, New Mexico +PO Box 5800 +Albuquerque, NM 87185 +USA + +UT-Battelle +1 Bethel Valley Rd +Oak Ridge, TN 37830 + +Los Alamos National Security, LLC +105 Central Park Square +Los Alamos, NM 87544 + +Kitware Inc. +28 Corporate Drive +Clifton Park, NY 12065 +USA + +University of California, Davis +One Shields Avenue +Davis, CA 95616 +USA + +Under the terms of Contract DE-AC04-94AL85000, there is a +non-exclusive license for use of this work by or on behalf of the +U.S. Government. + +Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +Laboratory (LANL), the U.S. Government retains certain rights in +this software. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the + distribution. + + * Neither the name of Kitware nor the names of any contributors may + be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +======================================================================== + +The following files and directories come from third parties. Check the +contents of these for details on the specifics of their respective +licenses. +- - - - - - - - - - - - - - - - - - - - - - - - do not remove this line +CMake/CheckCXX11Features.cmake +CMake/FindBoostHeaders.cmake +CMake/FindTBB.cmake +CMake/FindGLEW.cmake +vtkm/cont/tbb/internal/parallel_sort.h +vtkm/testing/OptionParser.h diff --git a/include/libcloudph++/blk_1m/formulae.hpp b/include/libcloudph++/blk_1m/formulae.hpp index 13c6da9f..bd1d6040 100644 --- a/include/libcloudph++/blk_1m/formulae.hpp +++ b/include/libcloudph++/blk_1m/formulae.hpp @@ -43,7 +43,7 @@ namespace libcloudphxx // Kessler evaporation rate // eq. 5c in Grabowski & Smolarkiewicz 1996 (multiplied by rho!) template - quantity::type, real_t> evaporation_rate( + quantity::type, real_t> evaporation_rate( quantity rv, quantity rvs, quantity rr, @@ -51,27 +51,27 @@ namespace libcloudphxx quantity p ) { - return + return (1 - rv / rvs) / rhod * ( - real_t(1.6) - + real_t(124.9) * std::pow( + real_t(1.6) + + real_t(124.9) * std::pow( real_t(1e-3) * rhod * rr * si::cubic_metres / si::kilograms, real_t(.2046) - ) + ) ) // ventilation factor TODO- move to ventil.hpp * std::pow( - real_t(1e-3) * rhod * rr * si::cubic_metres / si::kilograms, + real_t(1e-3) * rhod * rr * si::cubic_metres / si::kilograms, real_t(.525) - ) - / (real_t(5.4e2) - + real_t(2.55e5) - * (real_t(1) / (p / si::pascals) / rvs)) + ) + / (real_t(5.4e2) + + real_t(2.55e5) + * (real_t(1) / (p / si::pascals) / rvs)) / si::seconds * si::kilograms / si::cubic_metres; } // Kessler/Beard terminal velocity - // eq. 5d in Grabowski & Smolarkiewicz 1996 + // eq. 5d in Grabowski & Smolarkiewicz 1996 libcloudphxx_const(si::velocity, vterm_A, 36.34, si::metre_per_second) using inverse_density = divide_typeof_helper::type; @@ -83,87 +83,175 @@ namespace libcloudphxx const quantity &rhod, const quantity &rhod_0 ) { - return - vterm_A() + return + vterm_A() * real_t(std::pow( - (rhod * rr * vterm_B()), + (rhod * rr * vterm_B()), real_t(.1346) - ) + ) * sqrt(rhod_0 / rhod) ); } + // mean mass of ice A particle + // eq. A.15a in Grabowski 1999 + template + quantity mass_a( + const quantity &ria, + const quantity &T, + const quantity &rhod_0 + ) { + using namespace common; + //calculating the mass of small ice A particle: + quantity IWCS = std::min(std::min(real_t(1e-3)*si::dimensionless(), rhod_0*ria/(si::kilogram/si::cubic_metres)), real_t(2.52e-4)*std::pow(real_t(1e3)*rhod_0*ria/(si::kilogram/si::cubic_metres),real_t(0.837)) *si::dimensionless()) * si::kilograms / si::cubic_meters; + quantity::type, real_t> alpha = std::max(real_t(2.2e4), real_t(-4.99e3)-real_t(4.94e4)*std::log10(real_t(1e3)*IWCS/si::kilograms*si::cubic_meters)) / si::meters; + quantity m_as = real_t(2)*pi()*moist_air::rho_i()/(std::pow(alpha*si::meters,3)/si::cubic_meters); + //calculating the mass of large ice A particle: + quantity IWCL = rhod_0*ria - IWCS; + real_t a_mu = real_t(5.2) + real_t(1.3e-3)*(T/si::kelvin-273.16); + real_t b_mu = real_t(0.026) - real_t(1.2e-3)*(T/si::kelvin-273.16); + real_t a_sigma = real_t(0.47) + real_t(2.1e-3)*(T/si::kelvin-273.16); + real_t b_sigma = real_t(0.018) - real_t(2.1e-4)*(T/si::kelvin-273.16); + real_t mu = a_mu + b_mu*std::log10(real_t(1e3)*IWCL/si::kilograms*si::cubic_meters); + real_t sigma = a_sigma + b_sigma*std::log10(real_t(1e3)*IWCL/si::kilograms*si::cubic_meters); + quantity m_al = real_t(1.67e17)*pi()*moist_air::rho_i()*std::exp(3*mu +9/2*std::pow(sigma,2))*si::cubic_meters; + //mass of the average ice A particle: + real_t delta = IWCS/(ria*rhod_0); + quantity m_a = delta*m_as + (real_t(1)-delta)*m_al; + return m_a; + } + + // mean velocity of ice A particle + // eq. A.15b in Grabowski 1999 + template + quantity::type, real_t> velocity_a( + const quantity &ria, + const quantity &rhod_0 + ) { + using namespace common; + //velocity of small ice A particle: + quantity::type, real_t> v_as = real_t(0.1)*si::meters/si::seconds; + //velocity of large ice A particle: + quantity IWCS = std::min(std::min(real_t(1e-3)*si::dimensionless(), rhod_0*ria/(si::kilogram/si::cubic_metres)), real_t(2.52e-4)*std::pow(real_t(1e3)*rhod_0*ria/(si::kilogram/si::cubic_metres),real_t(0.837)) *si::dimensionless()) * si::kilograms / si::cubic_meters; + quantity IWCL = rhod_0*ria - IWCS; + quantity::type, real_t> v_al = real_t(0.9)+real_t(0.1)*std::log10(real_t(1e3)*IWCL); + //average velocity of ice A particle: + real_t delta = IWCS/(ria*rhod_0); + quantity::type, real_t> v_a = (delta*v_as + (real_t(1)-delta)*v_al)*std::pow(real_t(0.3)/rhod_0,real_t(0.5)); + return v_a; + } + //ice A homogeneous nucleation 1 + // eq. A.21a in Grabowski 1999 template quantity::type, real_t> hom_A_nucleation_1( const quantity &rv, const quantity &rvs, const quantity &rvsi, - const quantity &th, - const quantity &p, + const quantity &T, const quantity &dt ) { - const quantity T = th * std::pow((p/real_t(100000)/si::pascal), 0.286); const quantity taunuc = dt; // timescale for nucleation - real_t beta = (T/si::kelvins > real_t(213.16)*si::dimensionless()) ? real_t(0.1)+real_t(0.9)*(T-real_t(213.16)*si::kelvins)/(real_t(20)*si::kelvins) : real_t(0.1); + real_t beta = (T > real_t(213.16)*si::kelvins) ? real_t(0.1)+real_t(0.9)*(T-real_t(213.16)*si::kelvins)/(real_t(20)*si::kelvins) : real_t(0.1); quantity rv_adj = beta * rvs + (1-beta) * rvsi; - //homogeneous nucleation rate (only if T < -40 C): - quantity::type, real_t> hom1 = (T/si::kelvins < real_t(233.16)*si::dimensionless()) ? (real_t(1)-std::exp(-dt/taunuc)) * std::max(real_t(0)*si::dimensionless(), rv-rv_adj) / dt : real_t(0)/si::seconds; + //homogeneous nucleation rate (only if T < -40 C) + quantity::type, real_t> hom1 = (T < real_t(233.16)*si::kelvins) ? (real_t(1)-std::exp(-dt/taunuc)) * std::max(real_t(0)*si::dimensionless(), rv-rv_adj) / dt : real_t(0)/si::seconds; return hom1; } //ice A homogeneous nucleation 2 + // eq. A.21b in Grabowski 1999 template quantity::type, real_t> hom_A_nucleation_2( const quantity &rc, - const quantity &th, - const quantity &p, + const quantity &T, const quantity &dt ) { - const quantity T = th * std::pow((p/real_t(100000)/si::pascal), 0.286); const quantity taunuc = dt; // time scale for nucleation //homogeneous nucleation rate (only if T < -40 C) - quantity::type, real_t> hom2 = (T/si::kelvins < real_t(233.16)*si::dimensionless()) ? (real_t(1)-std::exp(-dt/taunuc)) * rc / dt : real_t(0)/si::seconds; + quantity::type, real_t> hom2 = (T < real_t(233.16)*si::kelvins) ? (real_t(1)-std::exp(-dt/taunuc)) * rc / dt : real_t(0)/si::seconds; return hom2; } //ice A heterogeneous nucleation + // eq. A.19 in Grabowski 1999 template quantity::type, real_t> het_A_nucleation( const quantity &ria, const quantity &rc, - const quantity &th, - const quantity &p, + const quantity &T, const quantity &rhod_0, const quantity &dt ) { - const quantity T = th * std::pow((p/real_t(100000)/si::pascal), 0.286); + using namespace common; const quantity taunuc = dt; // timescale for nucleation - //heterogeneous part - //calculating the mass of small ice A particle: - quantity IWCS = std::min(std::min(real_t(1e-3)*si::dimensionless(), rhod_0*ria/(si::kilogram/si::cubic_metres)), real_t(2.52e-4)*std::pow(real_t(1e3)*rhod_0*ria/(si::kilogram/si::cubic_metres),real_t(0.837)) *si::dimensionless()) * si::kilograms / si::cubic_meters; - quantity rho_i = real_t(916.8) *si::kilograms / si::cubic_meters; //ice density - quantity::type, real_t> alpha = std::max(real_t(2.2e4), real_t(-4.99e3)-real_t(4.94e4)*std::log10(real_t(1e3)*IWCS/si::kilograms*si::cubic_meters)) / si::meters; - quantity m_as = real_t(2)*pi()*rho_i/(std::pow(alpha*si::meters,3)/si::cubic_meters); - //calculating the mass of large ice A particle: - quantity IWCL = rhod_0*ria - IWCS; - real_t a_mu = real_t(5.2) + real_t(1.3e-3)*(T/si::kelvin-273.16); - real_t b_mu = real_t(0.026) - real_t(1.2e-3)*(T/si::kelvin-273.16); - real_t a_sigma = real_t(0.47) + real_t(2.1e-3)*(T/si::kelvin-273.16); - real_t b_sigma = real_t(0.018) - real_t(2.1e-4)*(T/si::kelvin-273.16); - real_t mu = a_mu + b_mu*std::log10(real_t(1e3)*IWCL/si::kilograms*si::cubic_meters); - real_t sigma = a_sigma + b_sigma*std::log10(real_t(1e3)*IWCL/si::kilograms*si::cubic_meters); - quantity m_al = real_t(1.67e17)*pi()*rho_i*std::exp(3*mu +9/2*std::pow(sigma,2))*si::cubic_meters; - //mass of average ice A particle - real_t delta = IWCS/(ria*rhod_0); - quantity ma = delta*m_as + (real_t(1)-delta)*m_al; + //mean ice A particle mass + quantity m_a = mass_a(ria, T, rhod_0); //concentration of ice nuclei: quantity::type, real_t> N_in = std::min(real_t(1e5), real_t(1e-2)*std::exp(real_t(0.6)*(real_t(273.16)-T/si::kelvins))) / si::cubic_meters; //nucleation rate: - quantity::type, real_t> het = (real_t(1)-std::exp(-dt/taunuc)) * std::min(rc, std::max(real_t(0)*si::dimensionless(), N_in*std::max(real_t(1e-12)*si::dimensionless(), ma/si::kilograms)/rhod_0*si::kilograms - ria)) / dt; + quantity::type, real_t> het = (real_t(1)-std::exp(-dt/taunuc)) * std::min(rc, std::max(real_t(0)*si::dimensionless(), N_in*std::max(real_t(1e-12)*si::dimensionless(), m_a/si::kilograms)/rhod_0*si::kilograms - ria)) / dt; return het; } + // ice B heterogeneous nucleation 1 + // eq. A.23 in Grabowski 1999 + template + quantity::type, real_t> het_B_nucleation_1( + const quantity &rr, + const quantity &ria, + const quantity &T, + const quantity &rhod_0 + ) { + using namespace common; + //parameters of Marshall & Palmer distribution for rain: + real_t N_0r = real_t(1e7); + real_t gamma_r = std::pow(pi()*moist_air::rho_w()*N_0r / (rhod_0*rr) , real_t(0.25)); + //mean raindrop mass + quantity m_r = pi()*moist_air::rho_w() / (real_t(6)*std::pow(gamma_r,3)) * si::cubic_meters; + //mean raindrop velocity + real_t v_r = real_t(251)*std::pow(gamma_r*rhod_0/si::kilograms*si::cubic_meters,real_t(-0.5)); + //mean raindrop radius + quantity R_r = real_t(0.5)/gamma_r * si::meters; + //mean ice A particle mass + quantity m_a = mass_a(ria, T, rhod_0); + //mean ice A velocity + real_t v_a = velocity_a(ria, rhod_0) *si::seconds/si::meters; + //rate of collisions between raindrops and ice A + real_t N_ra = N_0r/gamma_r*std::abs(v_r-v_a)*pi()*R_r*R_r*rhod_0*ria/m_a * si::meters; + //nucleation rate: + quantity::type, real_t> het1 = N_ra * m_r /si::kilograms /si::seconds; + return het1; + } + + //ice B heterogeneous nucleation 2 + // eq. A.23 in Grabowski 1999 + template + quantity::type, real_t> het_B_nucleation_2( + const quantity &rr, + const quantity &ria, + const quantity &T, + const quantity &rhod_0 + ) { + using namespace common; + //parameters of Marshall & Palmer distribution for rain: + real_t N_0r = real_t(1e7); + real_t gamma_r = std::pow(pi()*moist_air::rho_w()*N_0r / (rhod_0*rr) , real_t(0.25)); + //mean raindrop velocity + real_t v_r = real_t(251)*std::pow(gamma_r*rhod_0/si::kilograms*si::cubic_meters,real_t(-0.5)); + //mean raindrop radius + quantity R_r = real_t(0.5)/gamma_r * si::meters; + //mean ice A particle mass + quantity m_a = mass_a(ria, T, rhod_0); + //mean ice A velocity + real_t v_a = velocity_a(ria, rhod_0) *si::seconds/si::meters ; + //rate of collisions between raindrops and ice A + real_t N_ra = N_0r/gamma_r*abs(v_r-v_a)*pi()*R_r*R_r*rhod_0*ria/m_a * si::meters; + //nucleation rate: + quantity::type, real_t> het2 = N_ra * m_a /si::kilograms /si::seconds; + return het2; + } + }; }; }; diff --git a/include/libcloudph++/blk_1m/options.hpp b/include/libcloudph++/blk_1m/options.hpp index 792b7e44..cbcedafb 100644 --- a/include/libcloudph++/blk_1m/options.hpp +++ b/include/libcloudph++/blk_1m/options.hpp @@ -21,9 +21,10 @@ namespace libcloudphxx conv = true, // autoconversion accr = true, // accretion sedi = true, // sedimentation - homA1 = true, // homogeneous nucleation 1 of ice A - homA2 = true, // homogeneous nucleation 2 of ice A - hetA = true; // heterogeneous nucleation of ice A + homA1 = true, // homogeneous nucleation 1 of ice A + homA2 = true, // homogeneous nucleation 2 of ice A + hetA = true, // heterogeneous nucleation of ice A + hetB = true; // heterogeneous nucleation of ice B real_t r_c0 = 5e-4, // autoconv. threshold k_acnv = 0.001, // Kessler autoconversion (eq. 5a in Grabowski & Smolarkiewicz 1996) diff --git a/include/libcloudph++/blk_1m/rhs_cellwise.hpp b/include/libcloudph++/blk_1m/rhs_cellwise.hpp index af5b2fa6..24145305 100644 --- a/include/libcloudph++/blk_1m/rhs_cellwise.hpp +++ b/include/libcloudph++/blk_1m/rhs_cellwise.hpp @@ -135,6 +135,7 @@ namespace libcloudphxx cont_t &dot_rc_cont, cont_t &dot_rr_cont, cont_t &dot_ria_cont, + //cont_t &dot_rib_cont, const cont_t &rhod_cont, const cont_t &p_cont, const cont_t &th_cont, @@ -142,10 +143,12 @@ namespace libcloudphxx const cont_t &rc_cont, const cont_t &rr_cont, const cont_t &ria_cont, + //const cont_t &rib_cont, const real_t &dt ) { - rhs_cellwise(opts, dot_rc_cont, dot_rr_cont, rc_cont, rr_cont); + // autoconversion, collection and rain evaporation: + rhs_cellwise_nwtrph(opts, dot_th_cont, dot_rv_cont, dot_rc_cont, dot_rr_cont, rhod_cont, p_cont, th_cont, rv_cont, rc_cont, rr_cont, dt); // rain evaporation treated as a force in Newthon-Raphson saturation adjustment for (auto tup : zip(dot_th_cont, dot_rv_cont, dot_rc_cont, dot_rr_cont, dot_ria_cont, rhod_cont, p_cont, th_cont, rv_cont, rc_cont, rr_cont, ria_cont)) @@ -153,7 +156,6 @@ namespace libcloudphxx using namespace common; real_t - rr_to_rv = 0, rv_to_ria = 0, rc_to_ria = 0, &dot_th = std::get<0>(tup), @@ -181,20 +183,6 @@ namespace libcloudphxx real_t rvs = const_cp::r_vs(T, p); real_t rvsi = const_cp::r_vsi(T, p); - // rain evaporation - rr_to_rv += ( - formulae::evaporation_rate( - rv * si::dimensionless(), - rvs * si::dimensionless(), - rr * si::dimensionless(), - rhod, - p - ) * si::seconds * dt - ); - - // limiting - rr_to_rv = std::min(rr, rr_to_rv) / dt; - // ice A heterogeneous nucleation if (opts.hetA) { @@ -202,8 +190,7 @@ namespace libcloudphxx formulae::het_A_nucleation( ria * si::dimensionless(), rc * si::dimensionless(), - th, - p, + T, rhod, dt * si::seconds ) * si::seconds // to make it dimensionless @@ -218,8 +205,7 @@ namespace libcloudphxx rv * si::dimensionless(), rvs * si::dimensionless(), rvsi * si::dimensionless(), - th, - p, + T, dt * si::seconds ) * si::seconds // to make it dimensionless ); @@ -231,21 +217,19 @@ namespace libcloudphxx rc_to_ria += ( formulae::hom_A_nucleation_2( rc * si::dimensionless(), - th, - p, + T, dt * si::seconds ) * si::seconds // to make it dimensionless ); } - dot_rr -= rr_to_rv; dot_rc -= rc_to_ria; - dot_rv += rr_to_rv - rv_to_ria; + dot_rv -= rv_to_ria; dot_ria += rc_to_ria + rv_to_ria; - dot_th -= const_cp::l_v(T) / (moist_air::c_pd() * theta_std::exner(p)) * rr_to_rv / si::kelvins; //heat of vaporisation - dot_th += const_cp::l_sublim(T) / (moist_air::c_pd() * theta_std::exner(p)) * rv_to_ria / si::kelvins; //heat of sublimation - dot_th += const_cp::l_freez(T) / (moist_air::c_pd() * theta_std::exner(p)) * rc_to_ria / si::kelvins; //heat of freezing + dot_th += const_cp::l_s(T) / (moist_air::c_pd() * theta_std::exner(p)) * rv_to_ria / si::kelvins; //heat of sublimation + dot_th += const_cp::l_f(T) / (moist_air::c_pd() * theta_std::exner(p)) * rc_to_ria / si::kelvins; //heat of freezing } } + } } diff --git a/include/libcloudph++/common/const_cp.hpp b/include/libcloudph++/common/const_cp.hpp index 71151795..d2433437 100644 --- a/include/libcloudph++/common/const_cp.hpp +++ b/include/libcloudph++/common/const_cp.hpp @@ -21,9 +21,9 @@ namespace libcloudphxx // water triple point parameters libcloudphxx_const(si::pressure, p_tri, 611.73, si::pascals) // pressure libcloudphxx_const(si::temperature, T_tri, 273.16, si::kelvins) // temperature - libcloudphxx_const(energy_over_mass, l_tri_evap, 2.5e6, si::joules / si::kilograms) // latent heat of evaporation - libcloudphxx_const(energy_over_mass, l_tri_sublim, 2.83e6, si::joules / si::kilograms) // latent heat of sublimation - libcloudphxx_const(energy_over_mass, l_tri_freez, 3.34e5, si::joules / si::kilograms) // latent heat of freezing/melting + libcloudphxx_const(energy_over_mass, lv_tri, 2.5e6, si::joules / si::kilograms) // latent heat of evaporation + libcloudphxx_const(energy_over_mass, ls_tri, 2.834e6, si::joules / si::kilograms) // latent heat of sublimation + libcloudphxx_const(energy_over_mass, lf_tri, 3.34e5, si::joules / si::kilograms) // latent heat of freezing // saturation vapour pressure with respect to liquid water // assuming constant c_p_v and c_p_w with constants taken at triple point @@ -37,7 +37,7 @@ namespace libcloudphxx // { return p_tri() * exp( - (l_tri_evap() + (c_pw() - c_pv()) * T_tri()) / R_v() * (real_t(1) / T_tri() - real_t(1) / T) + (lv_tri() + (c_pw() - c_pv()) * T_tri()) / R_v() * (real_t(1) / T_tri() - real_t(1) / T) - (c_pw() - c_pv()) / R_v() * std::log(T / T_tri()) ); } @@ -51,7 +51,8 @@ namespace libcloudphxx ) { return p_tri() * exp( - (l_tri_sublim() + (c_pi() - c_pv()) * T_tri()) / R_v() * (real_t(1) / T_tri() - real_t(1) / T) + (ls_tri() + (c_pi() - c_pv()) * T_tri()) / R_v() * (real_t(1) / T_tri() - real_t(1) / T) + - (c_pi() - c_pv()) / R_v() * std::log(T / T_tri()) ); } @@ -81,25 +82,25 @@ namespace libcloudphxx quantity::type , real_t> l_v( const quantity &T ) { - return l_tri_evap() + (c_pv() - c_pw()) * (T - T_tri()); + return lv_tri() + (c_pv() - c_pw()) * (T - T_tri()); } // latent heat of sublimation for constant c_p template BOOST_GPU_ENABLED - quantity::type , real_t> l_sublim( + quantity::type , real_t> l_s( const quantity &T ) { - return l_tri_sublim() + (c_pv() - c_pi()) * (T - T_tri()); + return ls_tri() + (c_pv() - c_pi()) * (T - T_tri()); } - // latent heat of freezing/melting for constant c_p + // latent heat of freezing for constant c_p template BOOST_GPU_ENABLED - quantity::type , real_t> l_freez( + quantity::type , real_t> l_f( const quantity &T ) { - return l_tri_freez() + (c_pw() - c_pi()) * (T - T_tri()); + return lf_tri() + (c_pw() - c_pi()) * (T - T_tri()); } }; diff --git a/include/libcloudph++/common/moist_air.hpp b/include/libcloudph++/common/moist_air.hpp index e78dcea1..159b04cf 100644 --- a/include/libcloudph++/common/moist_air.hpp +++ b/include/libcloudph++/common/moist_air.hpp @@ -26,7 +26,7 @@ namespace libcloudphxx libcloudphxx_const(energy_over_temperature_mass, c_pd, 1005, si::joules / si::kilograms / si::kelvins) // dry air libcloudphxx_const(energy_over_temperature_mass, c_pv, 1850, si::joules / si::kilograms / si::kelvins) // water vapour libcloudphxx_const(energy_over_temperature_mass, c_pw, 4218, si::joules / si::kilograms / si::kelvins) // liquid water - libcloudphxx_const(energy_over_temperature_mass, c_pi, 2107, si::joules / si::kilograms / si::kelvins) // ice + libcloudphxx_const(energy_over_temperature_mass, c_pi, 2114, si::joules / si::kilograms / si::kelvins) // ice // ice // molar masses libcloudphxx_const(mass_over_amount, M_d, 0.02897, si::kilograms / si::moles) // dry air (Curry & Webster / Seinfeld & Pandis) @@ -48,6 +48,8 @@ namespace libcloudphxx // water density libcloudphxx_const(si::mass_density, rho_w, 1e3, si::kilograms / si::cubic_metres) + // ice density + libcloudphxx_const(si::mass_density, rho_i, 910, si::kilograms / si::cubic_metres) // mixing rule for extensive quantitites (i.e. using mass mixing ratio) template diff --git a/tests/python/unit/api_blk_1m.py b/tests/python/unit/api_blk_1m.py index 1c36b53a..c6939fa9 100644 --- a/tests/python/unit/api_blk_1m.py +++ b/tests/python/unit/api_blk_1m.py @@ -15,6 +15,7 @@ print("homA1 =", opts.homA1) print("homA2 =", opts.homA2) print("hetA =", opts.hetA) +print("hetB =", opts.hetB) print("r_c0 =", opts.r_c0) print("r_eps =", opts.r_eps) @@ -25,6 +26,7 @@ rv = arr_t([0. ]) rr = arr_t([0. ]) ria = arr_t([0. ]) +rib = arr_t([0. ]) dt = 1 dz = 1 @@ -90,5 +92,6 @@ dot_rr = arr_t([0.]) dot_rv = arr_t([0.]) dot_ria = arr_t([0.]) +dot_rib = arr_t([0.]) blk_1m.rhs_cellwise_nwtrph_ice(opts, dot_th, dot_rv, dot_rc, dot_rr, dot_ria, rhod, p, th, rv, rc, rr, ria, dt) assert dot_ria != 0 \ No newline at end of file From 06319fe202a72b9150e20c8aafc5006b4eebfafc Mon Sep 17 00:00:00 2001 From: Agnieszka Makulska Date: Thu, 31 Oct 2024 15:11:41 +0100 Subject: [PATCH 07/17] adding hetB to rhs_cellwise_nwtrph_ice --- bindings/python/blk_1m.hpp | 54 +++---- bindings/python/util.hpp | 1 + include/libcloudph++/blk_1m/formulae.hpp | 148 ++++++++++++------- include/libcloudph++/blk_1m/rhs_cellwise.hpp | 54 +++++-- tests/python/unit/api_blk_1m.py | 5 +- 5 files changed, 164 insertions(+), 98 deletions(-) diff --git a/bindings/python/blk_1m.hpp b/bindings/python/blk_1m.hpp index 00ce7de1..b9f9928c 100644 --- a/bindings/python/blk_1m.hpp +++ b/bindings/python/blk_1m.hpp @@ -29,10 +29,10 @@ namespace libcloudphxx { const typename arr_t::T_numtype &dt ) { arr_t - np2bz_th(np2bz(th)), - np2bz_rv(np2bz(rv)), - np2bz_rc(np2bz(rc)), - np2bz_rr(np2bz(rr)); + np2bz_th(np2bz(th)), + np2bz_rv(np2bz(rv)), + np2bz_rc(np2bz(rc)), + np2bz_rr(np2bz(rr)); b1m::adj_cellwise( opts, np2bz(rhod), // since it is const, it may be a temporary object @@ -56,10 +56,10 @@ namespace libcloudphxx { const typename arr_t::T_numtype &dt ) { arr_t - np2bz_th(np2bz(th)), - np2bz_rv(np2bz(rv)), - np2bz_rc(np2bz(rc)), - np2bz_rr(np2bz(rr)); + np2bz_th(np2bz(th)), + np2bz_rv(np2bz(rv)), + np2bz_rc(np2bz(rc)), + np2bz_rr(np2bz(rr)); b1m::adj_cellwise_constp( opts, np2bz(rhod), // since it is const, it may be a temporary object @@ -82,9 +82,9 @@ namespace libcloudphxx { const typename arr_t::T_numtype &dt ) { arr_t - np2bz_th(np2bz(th)), - np2bz_rv(np2bz(rv)), - np2bz_rc(np2bz(rc)); + np2bz_th(np2bz(th)), + np2bz_rv(np2bz(rv)), + np2bz_rc(np2bz(rc)); b1m::adj_cellwise_nwtrph( opts, np2bz(p), @@ -104,8 +104,8 @@ namespace libcloudphxx { const bp_array &rr ) { arr_t - np2bz_dot_rc(np2bz(dot_rc)), - np2bz_dot_rr(np2bz(dot_rr)); + np2bz_dot_rc(np2bz(dot_rc)), + np2bz_dot_rr(np2bz(dot_rr)); b1m::rhs_cellwise( opts, np2bz_dot_rc, @@ -131,10 +131,10 @@ namespace libcloudphxx { const typename arr_t::T_numtype &dt ) { arr_t - np2bz_dot_rc(np2bz(dot_rc)), - np2bz_dot_rr(np2bz(dot_rr)), - np2bz_dot_rv(np2bz(dot_rv)), - np2bz_dot_th(np2bz(dot_th)); + np2bz_dot_rc(np2bz(dot_rc)), + np2bz_dot_rr(np2bz(dot_rr)), + np2bz_dot_rv(np2bz(dot_rv)), + np2bz_dot_th(np2bz(dot_th)); b1m::rhs_cellwise_nwtrph( opts, np2bz_dot_th, @@ -160,7 +160,7 @@ namespace libcloudphxx { bp_array &dot_rc, bp_array &dot_rr, bp_array &dot_ria, - //bp_array &dot_rib, + bp_array &dot_rib, const bp_array &rhod, const bp_array &p, const bp_array &th, @@ -168,16 +168,16 @@ namespace libcloudphxx { const bp_array &rc, const bp_array &rr, const bp_array &ria, - //const bp_array &rib, + const bp_array &rib, const typename arr_t::T_numtype &dt ) { arr_t - np2bz_dot_ria(np2bz(dot_ria)), - //np2bz_dot_rib(np2bz(dot_rib)), - np2bz_dot_rc(np2bz(dot_rc)), - np2bz_dot_rr(np2bz(dot_rr)), - np2bz_dot_rv(np2bz(dot_rv)), - np2bz_dot_th(np2bz(dot_th)); + np2bz_dot_ria(np2bz(dot_ria)), + np2bz_dot_rib(np2bz(dot_rib)), + np2bz_dot_rc(np2bz(dot_rc)), + np2bz_dot_rr(np2bz(dot_rr)), + np2bz_dot_rv(np2bz(dot_rv)), + np2bz_dot_th(np2bz(dot_th)); b1m::rhs_cellwise_nwtrph_ice( opts, np2bz_dot_th, @@ -185,7 +185,7 @@ namespace libcloudphxx { np2bz_dot_rc, np2bz_dot_rr, np2bz_dot_ria, - //np2bz_dot_rib, + np2bz_dot_rib, np2bz(rhod), np2bz(p), np2bz(th), @@ -193,7 +193,7 @@ namespace libcloudphxx { np2bz(rc), np2bz(rr), np2bz(ria), - //np2bz(rib), + np2bz(rib), dt ); } diff --git a/bindings/python/util.hpp b/bindings/python/util.hpp index f6039857..78a617bd 100644 --- a/bindings/python/util.hpp +++ b/bindings/python/util.hpp @@ -12,6 +12,7 @@ #include // otherwise Clang fails in debug mode #include +#define BOOST_PYTHON_MAX_ARITY 20 // max number of arguments in a function #include #ifdef BPNUMPY #include diff --git a/include/libcloudph++/blk_1m/formulae.hpp b/include/libcloudph++/blk_1m/formulae.hpp index bd1d6040..ca9d9994 100644 --- a/include/libcloudph++/blk_1m/formulae.hpp +++ b/include/libcloudph++/blk_1m/formulae.hpp @@ -95,55 +95,76 @@ namespace libcloudphxx // mean mass of ice A particle // eq. A.15a in Grabowski 1999 - template + template quantity mass_a( - const quantity &ria, - const quantity &T, - const quantity &rhod_0 + const quantity &ria, + const quantity &T, + const quantity &rhod_0 ) { using namespace common; //calculating the mass of small ice A particle: - quantity IWCS = std::min(std::min(real_t(1e-3)*si::dimensionless(), rhod_0*ria/(si::kilogram/si::cubic_metres)), real_t(2.52e-4)*std::pow(real_t(1e3)*rhod_0*ria/(si::kilogram/si::cubic_metres),real_t(0.837)) *si::dimensionless()) * si::kilograms / si::cubic_meters; - quantity::type, real_t> alpha = std::max(real_t(2.2e4), real_t(-4.99e3)-real_t(4.94e4)*std::log10(real_t(1e3)*IWCS/si::kilograms*si::cubic_meters)) / si::meters; - quantity m_as = real_t(2)*pi()*moist_air::rho_i()/(std::pow(alpha*si::meters,3)/si::cubic_meters); + quantity IWCS = std::min( + std::min(real_t(1e-3) * si::dimensionless(), + rhod_0 * ria / (si::kilogram / si::cubic_metres)), + real_t(2.52e-4) * std::pow( + real_t(1e3) * rhod_0 * ria / (si::kilogram / si::cubic_metres), + real_t(0.837)) * si::dimensionless()) * si::kilograms / + si::cubic_meters; + quantity::type, real_t> alpha = std::max(real_t(2.2e4), + real_t(-4.99e3) - real_t(4.94e4) * std::log10(real_t(1e3) * IWCS / si::kilograms * si::cubic_meters)) / + si::meters; + quantity m_as = real_t(2) * pi() * moist_air::rho_i() / ( + std::pow(alpha * si::meters, 3) / si::cubic_meters); //calculating the mass of large ice A particle: - quantity IWCL = rhod_0*ria - IWCS; - real_t a_mu = real_t(5.2) + real_t(1.3e-3)*(T/si::kelvin-273.16); - real_t b_mu = real_t(0.026) - real_t(1.2e-3)*(T/si::kelvin-273.16); - real_t a_sigma = real_t(0.47) + real_t(2.1e-3)*(T/si::kelvin-273.16); - real_t b_sigma = real_t(0.018) - real_t(2.1e-4)*(T/si::kelvin-273.16); - real_t mu = a_mu + b_mu*std::log10(real_t(1e3)*IWCL/si::kilograms*si::cubic_meters); - real_t sigma = a_sigma + b_sigma*std::log10(real_t(1e3)*IWCL/si::kilograms*si::cubic_meters); - quantity m_al = real_t(1.67e17)*pi()*moist_air::rho_i()*std::exp(3*mu +9/2*std::pow(sigma,2))*si::cubic_meters; + quantity IWCL = rhod_0 * ria - IWCS; + real_t a_mu = real_t(5.2) + real_t(1.3e-3) * (T / si::kelvin - 273.16); + real_t b_mu = real_t(0.026) - real_t(1.2e-3) * (T / si::kelvin - 273.16); + real_t a_sigma = real_t(0.47) + real_t(2.1e-3) * (T / si::kelvin - 273.16); + real_t b_sigma = real_t(0.018) - real_t(2.1e-4) * (T / si::kelvin - 273.16); + real_t mu = a_mu + b_mu * std::log10(real_t(1e3) * IWCL / si::kilograms * si::cubic_meters); + real_t sigma = a_sigma + b_sigma * std::log10(real_t(1e3) * IWCL / si::kilograms * si::cubic_meters); + quantity m_al = real_t(1.67e17) * pi() * moist_air::rho_i() * std::exp( + 3 * mu + 9 / 2 * std::pow(sigma, 2)) * si::cubic_meters; //mass of the average ice A particle: - real_t delta = IWCS/(ria*rhod_0); - quantity m_a = delta*m_as + (real_t(1)-delta)*m_al; + real_t delta = IWCS / (ria * rhod_0); + quantity m_a = delta * m_as + (real_t(1) - delta) * m_al; return m_a; } // mean velocity of ice A particle // eq. A.15b in Grabowski 1999 - template + template quantity::type, real_t> velocity_a( - const quantity &ria, - const quantity &rhod_0 + const quantity &ria, + const quantity &rhod_0 ) { using namespace common; //velocity of small ice A particle: - quantity::type, real_t> v_as = real_t(0.1)*si::meters/si::seconds; + quantity::type, real_t> v_as = + real_t(0.1) * si::meters / si::seconds; //velocity of large ice A particle: - quantity IWCS = std::min(std::min(real_t(1e-3)*si::dimensionless(), rhod_0*ria/(si::kilogram/si::cubic_metres)), real_t(2.52e-4)*std::pow(real_t(1e3)*rhod_0*ria/(si::kilogram/si::cubic_metres),real_t(0.837)) *si::dimensionless()) * si::kilograms / si::cubic_meters; - quantity IWCL = rhod_0*ria - IWCS; - quantity::type, real_t> v_al = real_t(0.9)+real_t(0.1)*std::log10(real_t(1e3)*IWCL); + quantity IWCS = std::min( + std::min(real_t(1e-3) * si::dimensionless(), + rhod_0 * ria / (si::kilogram / si::cubic_metres)), + real_t(2.52e-4) * std::pow( + real_t(1e3) * rhod_0 * ria / (si::kilogram / si::cubic_metres), + real_t(0.837)) * si::dimensionless()) * si::kilograms / + si::cubic_meters; + quantity IWCL = rhod_0 * ria - IWCS; + quantity::type, real_t> v_al = + (real_t(0.9) + real_t(0.1) * std::log10(real_t(1e3) * IWCL / si::kilograms * si::cubic_meters)) * si::meters + / si::seconds; //average velocity of ice A particle: - real_t delta = IWCS/(ria*rhod_0); - quantity::type, real_t> v_a = (delta*v_as + (real_t(1)-delta)*v_al)*std::pow(real_t(0.3)/rhod_0,real_t(0.5)); + real_t delta = IWCS / (ria * rhod_0); + quantity::type, real_t> v_a = + (delta * v_as + (real_t(1) - delta) * v_al) * std::pow( + real_t(0.3) / rhod_0 * si::kilograms / si::cubic_meters, real_t(0.5)); return v_a; } - //ice A homogeneous nucleation 1 + //ice A homogeneous nucleation 1 (rv -> ria) // eq. A.21a in Grabowski 1999 - template + template quantity::type, real_t> hom_A_nucleation_1( const quantity &rv, const quantity &rvs, @@ -152,16 +173,21 @@ namespace libcloudphxx const quantity &dt ) { const quantity taunuc = dt; // timescale for nucleation - real_t beta = (T > real_t(213.16)*si::kelvins) ? real_t(0.1)+real_t(0.9)*(T-real_t(213.16)*si::kelvins)/(real_t(20)*si::kelvins) : real_t(0.1); - quantity rv_adj = beta * rvs + (1-beta) * rvsi; + real_t beta = (T > real_t(213.16) * si::kelvins) + ? real_t(0.1) + real_t(0.9) * (T - real_t(213.16) * si::kelvins) / (real_t(20) * si::kelvins) + : real_t(0.1); + quantity rv_adj = beta * rvs + (1 - beta) * rvsi; //homogeneous nucleation rate (only if T < -40 C) - quantity::type, real_t> hom1 = (T < real_t(233.16)*si::kelvins) ? (real_t(1)-std::exp(-dt/taunuc)) * std::max(real_t(0)*si::dimensionless(), rv-rv_adj) / dt : real_t(0)/si::seconds; + quantity::type, real_t> hom1 = + (T < real_t(233.16) * si::kelvins) + ? (real_t(1) - std::exp(-dt / taunuc)) * std::max(real_t(0) * si::dimensionless(), rv - rv_adj) / dt + : real_t(0) / si::seconds; return hom1; } - //ice A homogeneous nucleation 2 + //ice A homogeneous nucleation 2 (rc -> ria) // eq. A.21b in Grabowski 1999 - template + template quantity::type, real_t> hom_A_nucleation_2( const quantity &rc, const quantity &T, @@ -169,13 +195,16 @@ namespace libcloudphxx ) { const quantity taunuc = dt; // time scale for nucleation //homogeneous nucleation rate (only if T < -40 C) - quantity::type, real_t> hom2 = (T < real_t(233.16)*si::kelvins) ? (real_t(1)-std::exp(-dt/taunuc)) * rc / dt : real_t(0)/si::seconds; + quantity::type, real_t> hom2 = + (T < real_t(233.16) * si::kelvins) + ? (real_t(1) - std::exp(-dt / taunuc)) * rc / dt + : real_t(0) / si::seconds; return hom2; } - //ice A heterogeneous nucleation + //ice A heterogeneous nucleation (rc-> ria) // eq. A.19 in Grabowski 1999 - template + template quantity::type, real_t> het_A_nucleation( const quantity &ria, const quantity &rc, @@ -188,15 +217,21 @@ namespace libcloudphxx //mean ice A particle mass quantity m_a = mass_a(ria, T, rhod_0); //concentration of ice nuclei: - quantity::type, real_t> N_in = std::min(real_t(1e5), real_t(1e-2)*std::exp(real_t(0.6)*(real_t(273.16)-T/si::kelvins))) / si::cubic_meters; + quantity::type, real_t> N_in = std::min(real_t(1e5), + real_t(1e-2) * std::exp(real_t(0.6) * (real_t(273.16) - T / si::kelvins))) / si::cubic_meters; //nucleation rate: - quantity::type, real_t> het = (real_t(1)-std::exp(-dt/taunuc)) * std::min(rc, std::max(real_t(0)*si::dimensionless(), N_in*std::max(real_t(1e-12)*si::dimensionless(), m_a/si::kilograms)/rhod_0*si::kilograms - ria)) / dt; + quantity::type, real_t> het = + (real_t(1) - std::exp(-dt / taunuc)) * std::min(rc, std::max(real_t(0) * si::dimensionless(), + N_in * std::max( + real_t(1e-12) * si::dimensionless(), + m_a / si::kilograms) / rhod_0 * si::kilograms + - ria)) / dt; return het; } - // ice B heterogeneous nucleation 1 + // ice B heterogeneous nucleation 1 (rr -> rib) // eq. A.23 in Grabowski 1999 - template + template quantity::type, real_t> het_B_nucleation_1( const quantity &rr, const quantity &ria, @@ -206,27 +241,29 @@ namespace libcloudphxx using namespace common; //parameters of Marshall & Palmer distribution for rain: real_t N_0r = real_t(1e7); - real_t gamma_r = std::pow(pi()*moist_air::rho_w()*N_0r / (rhod_0*rr) , real_t(0.25)); + real_t gamma_r = std::pow(pi() * moist_air::rho_w() * N_0r / (rhod_0 * rr), real_t(0.25)); //mean raindrop mass - quantity m_r = pi()*moist_air::rho_w() / (real_t(6)*std::pow(gamma_r,3)) * si::cubic_meters; + quantity m_r = pi() * moist_air::rho_w() / (real_t(6) * std::pow(gamma_r, 3)) + * si::cubic_meters; //mean raindrop velocity - real_t v_r = real_t(251)*std::pow(gamma_r*rhod_0/si::kilograms*si::cubic_meters,real_t(-0.5)); + real_t v_r = real_t(251) * std::pow(gamma_r * rhod_0 / si::kilograms * si::cubic_meters, real_t(-0.5)); //mean raindrop radius - quantity R_r = real_t(0.5)/gamma_r * si::meters; + quantity R_r = real_t(0.5) / gamma_r * si::meters; //mean ice A particle mass quantity m_a = mass_a(ria, T, rhod_0); //mean ice A velocity - real_t v_a = velocity_a(ria, rhod_0) *si::seconds/si::meters; + real_t v_a = velocity_a(ria, rhod_0) * si::seconds / si::meters; //rate of collisions between raindrops and ice A - real_t N_ra = N_0r/gamma_r*std::abs(v_r-v_a)*pi()*R_r*R_r*rhod_0*ria/m_a * si::meters; + real_t N_ra = N_0r / gamma_r * std::abs(v_r - v_a) * pi() * R_r * R_r * rhod_0 * ria / m_a * si::meters; //nucleation rate: - quantity::type, real_t> het1 = N_ra * m_r /si::kilograms /si::seconds; + quantity::type, real_t> het1 = + N_ra * m_r / si::kilograms / si::seconds; return het1; } - //ice B heterogeneous nucleation 2 + //ice B heterogeneous nucleation 2 (ria -> rib) // eq. A.23 in Grabowski 1999 - template + template quantity::type, real_t> het_B_nucleation_2( const quantity &rr, const quantity &ria, @@ -236,19 +273,20 @@ namespace libcloudphxx using namespace common; //parameters of Marshall & Palmer distribution for rain: real_t N_0r = real_t(1e7); - real_t gamma_r = std::pow(pi()*moist_air::rho_w()*N_0r / (rhod_0*rr) , real_t(0.25)); + real_t gamma_r = std::pow(pi() * moist_air::rho_w() * N_0r / (rhod_0 * rr), real_t(0.25)); //mean raindrop velocity - real_t v_r = real_t(251)*std::pow(gamma_r*rhod_0/si::kilograms*si::cubic_meters,real_t(-0.5)); + real_t v_r = real_t(251) * std::pow(gamma_r * rhod_0 / si::kilograms * si::cubic_meters, real_t(-0.5)); //mean raindrop radius - quantity R_r = real_t(0.5)/gamma_r * si::meters; + quantity R_r = real_t(0.5) / gamma_r * si::meters; //mean ice A particle mass quantity m_a = mass_a(ria, T, rhod_0); //mean ice A velocity - real_t v_a = velocity_a(ria, rhod_0) *si::seconds/si::meters ; + real_t v_a = velocity_a(ria, rhod_0) * si::seconds / si::meters; //rate of collisions between raindrops and ice A - real_t N_ra = N_0r/gamma_r*abs(v_r-v_a)*pi()*R_r*R_r*rhod_0*ria/m_a * si::meters; + real_t N_ra = N_0r / gamma_r * abs(v_r - v_a) * pi() * R_r * R_r * rhod_0 * ria / m_a * si::meters; //nucleation rate: - quantity::type, real_t> het2 = N_ra * m_a /si::kilograms /si::seconds; + quantity::type, real_t> het2 = + N_ra * m_a / si::kilograms / si::seconds; return het2; } diff --git a/include/libcloudph++/blk_1m/rhs_cellwise.hpp b/include/libcloudph++/blk_1m/rhs_cellwise.hpp index 24145305..ca6f88cb 100644 --- a/include/libcloudph++/blk_1m/rhs_cellwise.hpp +++ b/include/libcloudph++/blk_1m/rhs_cellwise.hpp @@ -135,7 +135,7 @@ namespace libcloudphxx cont_t &dot_rc_cont, cont_t &dot_rr_cont, cont_t &dot_ria_cont, - //cont_t &dot_rib_cont, + cont_t &dot_rib_cont, const cont_t &rhod_cont, const cont_t &p_cont, const cont_t &th_cont, @@ -143,41 +143,44 @@ namespace libcloudphxx const cont_t &rc_cont, const cont_t &rr_cont, const cont_t &ria_cont, - //const cont_t &rib_cont, + const cont_t &rib_cont, const real_t &dt ) { // autoconversion, collection and rain evaporation: rhs_cellwise_nwtrph(opts, dot_th_cont, dot_rv_cont, dot_rc_cont, dot_rr_cont, rhod_cont, p_cont, th_cont, rv_cont, rc_cont, rr_cont, dt); - // rain evaporation treated as a force in Newthon-Raphson saturation adjustment - for (auto tup : zip(dot_th_cont, dot_rv_cont, dot_rc_cont, dot_rr_cont, dot_ria_cont, rhod_cont, p_cont, th_cont, rv_cont, rc_cont, rr_cont, ria_cont)) + for (auto tup : zip(dot_th_cont, dot_rv_cont, dot_rc_cont, dot_rr_cont, dot_ria_cont, dot_rib_cont, rhod_cont, p_cont, th_cont, rv_cont, rc_cont, rr_cont, ria_cont, rib_cont)) { using namespace common; real_t rv_to_ria = 0, rc_to_ria = 0, + rr_to_rib = 0, + ria_to_rib = 0, &dot_th = std::get<0>(tup), &dot_rv = std::get<1>(tup), &dot_rc = std::get<2>(tup), &dot_rr = std::get<3>(tup), - &dot_ria = std::get<4>(tup); + &dot_ria = std::get<4>(tup), + &dot_rib = std::get<5>(tup); const quantity - rhod = std::get<5>(tup) * si::kilograms / si::cubic_metres; + rhod = std::get<6>(tup) * si::kilograms / si::cubic_metres; const quantity - p = std::get<6>(tup) * si::pascals; + p = std::get<7>(tup) * si::pascals; const quantity - th = std::get<7>(tup) * si::kelvins; + th = std::get<8>(tup) * si::kelvins; const real_t - &rv = std::get<8>(tup), - &rc = std::get<9>(tup), - &rr = std::get<10>(tup), - &ria = std::get<11>(tup); + &rv = std::get<9>(tup), + &rc = std::get<10>(tup), + &rr = std::get<11>(tup), + &ria = std::get<12>(tup), + &rib = std::get<13>(tup); quantity T = th * theta_std::exner(p); real_t rvs = const_cp::r_vs(T, p); @@ -223,11 +226,34 @@ namespace libcloudphxx ); } + // ice B heterogeneous nucleation + if (opts.hetB) + { + rr_to_rib += ( + formulae::het_B_nucleation_1( + rr * si::dimensionless(), + ria * si::dimensionless(), + T, + rhod + ) * si::seconds // to make it dimensionless + ); + ria_to_rib += ( + formulae::het_B_nucleation_2( + rr * si::dimensionless(), + ria * si::dimensionless(), + T, + rhod + ) * si::seconds // to make it dimensionless + ); + } + dot_rc -= rc_to_ria; dot_rv -= rv_to_ria; - dot_ria += rc_to_ria + rv_to_ria; + dot_rr -= rr_to_rib; + dot_ria += rc_to_ria + rv_to_ria - ria_to_rib; + dot_rib += rr_to_rib + ria_to_rib; dot_th += const_cp::l_s(T) / (moist_air::c_pd() * theta_std::exner(p)) * rv_to_ria / si::kelvins; //heat of sublimation - dot_th += const_cp::l_f(T) / (moist_air::c_pd() * theta_std::exner(p)) * rc_to_ria / si::kelvins; //heat of freezing + dot_th += const_cp::l_f(T) / (moist_air::c_pd() * theta_std::exner(p)) * (rc_to_ria+rr_to_rib) / si::kelvins; //heat of freezing } } diff --git a/tests/python/unit/api_blk_1m.py b/tests/python/unit/api_blk_1m.py index c6939fa9..917fd3f1 100644 --- a/tests/python/unit/api_blk_1m.py +++ b/tests/python/unit/api_blk_1m.py @@ -93,5 +93,6 @@ dot_rv = arr_t([0.]) dot_ria = arr_t([0.]) dot_rib = arr_t([0.]) -blk_1m.rhs_cellwise_nwtrph_ice(opts, dot_th, dot_rv, dot_rc, dot_rr, dot_ria, rhod, p, th, rv, rc, rr, ria, dt) -assert dot_ria != 0 \ No newline at end of file +blk_1m.rhs_cellwise_nwtrph_ice(opts, dot_th, dot_rv, dot_rc, dot_rr, dot_ria, dot_rib, rhod, p, th, rv, rc, rr, ria, rib, dt) +assert dot_ria != 0 +assert dot_rib != 0 \ No newline at end of file From f89dc38c1c65c7702c3c14a6a35f8f027998d5da Mon Sep 17 00:00:00 2001 From: Agnieszka Makulska Date: Mon, 4 Nov 2024 13:04:38 +0100 Subject: [PATCH 08/17] renaming latent heat of evaporation back to l_tri --- include/libcloudph++/common/const_cp.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/libcloudph++/common/const_cp.hpp b/include/libcloudph++/common/const_cp.hpp index d2433437..ab9a2abc 100644 --- a/include/libcloudph++/common/const_cp.hpp +++ b/include/libcloudph++/common/const_cp.hpp @@ -21,7 +21,7 @@ namespace libcloudphxx // water triple point parameters libcloudphxx_const(si::pressure, p_tri, 611.73, si::pascals) // pressure libcloudphxx_const(si::temperature, T_tri, 273.16, si::kelvins) // temperature - libcloudphxx_const(energy_over_mass, lv_tri, 2.5e6, si::joules / si::kilograms) // latent heat of evaporation + libcloudphxx_const(energy_over_mass, l_tri, 2.5e6, si::joules / si::kilograms) // latent heat of evaporation libcloudphxx_const(energy_over_mass, ls_tri, 2.834e6, si::joules / si::kilograms) // latent heat of sublimation libcloudphxx_const(energy_over_mass, lf_tri, 3.34e5, si::joules / si::kilograms) // latent heat of freezing @@ -37,7 +37,7 @@ namespace libcloudphxx // { return p_tri() * exp( - (lv_tri() + (c_pw() - c_pv()) * T_tri()) / R_v() * (real_t(1) / T_tri() - real_t(1) / T) + (l_tri() + (c_pw() - c_pv()) * T_tri()) / R_v() * (real_t(1) / T_tri() - real_t(1) / T) - (c_pw() - c_pv()) / R_v() * std::log(T / T_tri()) ); } @@ -82,7 +82,7 @@ namespace libcloudphxx quantity::type , real_t> l_v( const quantity &T ) { - return lv_tri() + (c_pv() - c_pw()) * (T - T_tri()); + return l_tri() + (c_pv() - c_pw()) * (T - T_tri()); } // latent heat of sublimation for constant c_p From d41918a704e7aa9cf6ef1f2194aa10553c89447d Mon Sep 17 00:00:00 2001 From: Agnieszka Makulska Date: Mon, 4 Nov 2024 19:39:03 +0100 Subject: [PATCH 09/17] melting of ice A --- bindings/python/blk_1m.hpp | 48 +++---- bindings/python/lib.cpp | 2 + include/libcloudph++/blk_1m/formulae.hpp | 124 +++++++++++-------- include/libcloudph++/blk_1m/options.hpp | 4 +- include/libcloudph++/blk_1m/rhs_cellwise.hpp | 19 ++- tests/python/unit/api_blk_1m.py | 2 + 6 files changed, 120 insertions(+), 79 deletions(-) diff --git a/bindings/python/blk_1m.hpp b/bindings/python/blk_1m.hpp index b9f9928c..64b907aa 100644 --- a/bindings/python/blk_1m.hpp +++ b/bindings/python/blk_1m.hpp @@ -29,10 +29,10 @@ namespace libcloudphxx { const typename arr_t::T_numtype &dt ) { arr_t - np2bz_th(np2bz(th)), - np2bz_rv(np2bz(rv)), - np2bz_rc(np2bz(rc)), - np2bz_rr(np2bz(rr)); + np2bz_th(np2bz(th)), + np2bz_rv(np2bz(rv)), + np2bz_rc(np2bz(rc)), + np2bz_rr(np2bz(rr)); b1m::adj_cellwise( opts, np2bz(rhod), // since it is const, it may be a temporary object @@ -56,10 +56,10 @@ namespace libcloudphxx { const typename arr_t::T_numtype &dt ) { arr_t - np2bz_th(np2bz(th)), - np2bz_rv(np2bz(rv)), - np2bz_rc(np2bz(rc)), - np2bz_rr(np2bz(rr)); + np2bz_th(np2bz(th)), + np2bz_rv(np2bz(rv)), + np2bz_rc(np2bz(rc)), + np2bz_rr(np2bz(rr)); b1m::adj_cellwise_constp( opts, np2bz(rhod), // since it is const, it may be a temporary object @@ -82,9 +82,9 @@ namespace libcloudphxx { const typename arr_t::T_numtype &dt ) { arr_t - np2bz_th(np2bz(th)), - np2bz_rv(np2bz(rv)), - np2bz_rc(np2bz(rc)); + np2bz_th(np2bz(th)), + np2bz_rv(np2bz(rv)), + np2bz_rc(np2bz(rc)); b1m::adj_cellwise_nwtrph( opts, np2bz(p), @@ -104,8 +104,8 @@ namespace libcloudphxx { const bp_array &rr ) { arr_t - np2bz_dot_rc(np2bz(dot_rc)), - np2bz_dot_rr(np2bz(dot_rr)); + np2bz_dot_rc(np2bz(dot_rc)), + np2bz_dot_rr(np2bz(dot_rr)); b1m::rhs_cellwise( opts, np2bz_dot_rc, @@ -131,10 +131,10 @@ namespace libcloudphxx { const typename arr_t::T_numtype &dt ) { arr_t - np2bz_dot_rc(np2bz(dot_rc)), - np2bz_dot_rr(np2bz(dot_rr)), - np2bz_dot_rv(np2bz(dot_rv)), - np2bz_dot_th(np2bz(dot_th)); + np2bz_dot_rc(np2bz(dot_rc)), + np2bz_dot_rr(np2bz(dot_rr)), + np2bz_dot_rv(np2bz(dot_rv)), + np2bz_dot_th(np2bz(dot_th)); b1m::rhs_cellwise_nwtrph( opts, np2bz_dot_th, @@ -172,12 +172,12 @@ namespace libcloudphxx { const typename arr_t::T_numtype &dt ) { arr_t - np2bz_dot_ria(np2bz(dot_ria)), - np2bz_dot_rib(np2bz(dot_rib)), - np2bz_dot_rc(np2bz(dot_rc)), - np2bz_dot_rr(np2bz(dot_rr)), - np2bz_dot_rv(np2bz(dot_rv)), - np2bz_dot_th(np2bz(dot_th)); + np2bz_dot_ria(np2bz(dot_ria)), + np2bz_dot_rib(np2bz(dot_rib)), + np2bz_dot_rc(np2bz(dot_rc)), + np2bz_dot_rr(np2bz(dot_rr)), + np2bz_dot_rv(np2bz(dot_rv)), + np2bz_dot_th(np2bz(dot_th)); b1m::rhs_cellwise_nwtrph_ice( opts, np2bz_dot_th, @@ -207,7 +207,7 @@ namespace libcloudphxx { const typename arr_t::T_numtype &dz ) { arr_t - np2bz_dot_rr(np2bz(dot_rr)); + np2bz_dot_rr(np2bz(dot_rr)); return b1m::rhs_columnwise( opts, np2bz_dot_rr, diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index 76bde5ac..b363a184 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -159,6 +159,8 @@ BOOST_PYTHON_MODULE(libcloudphxx) .def_readwrite("homA2", &b1m::opts_t::homA2) .def_readwrite("hetA", &b1m::opts_t::hetA) .def_readwrite("hetB", &b1m::opts_t::hetB) + .def_readwrite("melA", &b1m::opts_t::melA) + .def_readwrite("melB", &b1m::opts_t::melB) .def_readwrite("r_c0", &b1m::opts_t::r_c0) .def_readwrite("r_eps", &b1m::opts_t::r_eps) .def_readwrite("nwtrph_iters", &b1m::opts_t::nwtrph_iters) diff --git a/include/libcloudph++/blk_1m/formulae.hpp b/include/libcloudph++/blk_1m/formulae.hpp index ca9d9994..b941ce54 100644 --- a/include/libcloudph++/blk_1m/formulae.hpp +++ b/include/libcloudph++/blk_1m/formulae.hpp @@ -9,6 +9,7 @@ #pragma once #include "../common/moist_air.hpp" +#include "../common/vterm.hpp" namespace libcloudphxx { @@ -25,7 +26,7 @@ namespace libcloudphxx const quantity rc_thresh, const quantity k_autoconv ) { - return k_autoconv * std::max( real_t(0) * si::dimensionless(), rc - rc_thresh); + return k_autoconv * std::max(real_t(0) * si::dimensionless(), rc - rc_thresh); } //Kessler collection @@ -49,8 +50,7 @@ namespace libcloudphxx quantity rr, quantity rhod, quantity p - ) - { + ) { return (1 - rv / rvs) / rhod * ( @@ -65,8 +65,8 @@ namespace libcloudphxx real_t(.525) ) / (real_t(5.4e2) - + real_t(2.55e5) - * (real_t(1) / (p / si::pascals) / rvs)) + + real_t(2.55e5) + * (real_t(1) / (p / si::pascals) / rvs)) / si::seconds * si::kilograms / si::cubic_metres; } @@ -74,7 +74,7 @@ namespace libcloudphxx // eq. 5d in Grabowski & Smolarkiewicz 1996 libcloudphxx_const(si::velocity, vterm_A, 36.34, si::metre_per_second) - using inverse_density = divide_typeof_helper::type; + using inverse_density = divide_typeof_helper::type; libcloudphxx_const(inverse_density, vterm_B, 1e-3, si::cubic_metres / si::kilograms) template @@ -86,16 +86,16 @@ namespace libcloudphxx return vterm_A() * real_t(std::pow( - (rhod * rr * vterm_B()), - real_t(.1346) - ) - * sqrt(rhod_0 / rhod) - ); + (rhod * rr * vterm_B()), + real_t(.1346) + ) + * sqrt(rhod_0 / rhod) + ); } // mean mass of ice A particle // eq. A.15a in Grabowski 1999 - template + template quantity mass_a( const quantity &ria, const quantity &T, @@ -104,17 +104,17 @@ namespace libcloudphxx using namespace common; //calculating the mass of small ice A particle: quantity IWCS = std::min( - std::min(real_t(1e-3) * si::dimensionless(), - rhod_0 * ria / (si::kilogram / si::cubic_metres)), - real_t(2.52e-4) * std::pow( - real_t(1e3) * rhod_0 * ria / (si::kilogram / si::cubic_metres), - real_t(0.837)) * si::dimensionless()) * si::kilograms / - si::cubic_meters; + std::min(real_t(1e-3) * si::dimensionless(), + rhod_0 * ria / (si::kilogram / si::cubic_metres)), + real_t(2.52e-4) * std::pow( + real_t(1e3) * rhod_0 * ria / (si::kilogram / si::cubic_metres), + real_t(0.837)) * si::dimensionless()) * si::kilograms / + si::cubic_meters; quantity::type, real_t> alpha = std::max(real_t(2.2e4), real_t(-4.99e3) - real_t(4.94e4) * std::log10(real_t(1e3) * IWCS / si::kilograms * si::cubic_meters)) / si::meters; quantity m_as = real_t(2) * pi() * moist_air::rho_i() / ( - std::pow(alpha * si::meters, 3) / si::cubic_meters); + std::pow(alpha * si::meters, 3) / si::cubic_meters); //calculating the mass of large ice A particle: quantity IWCL = rhod_0 * ria - IWCS; real_t a_mu = real_t(5.2) + real_t(1.3e-3) * (T / si::kelvin - 273.16); @@ -124,7 +124,7 @@ namespace libcloudphxx real_t mu = a_mu + b_mu * std::log10(real_t(1e3) * IWCL / si::kilograms * si::cubic_meters); real_t sigma = a_sigma + b_sigma * std::log10(real_t(1e3) * IWCL / si::kilograms * si::cubic_meters); quantity m_al = real_t(1.67e17) * pi() * moist_air::rho_i() * std::exp( - 3 * mu + 9 / 2 * std::pow(sigma, 2)) * si::cubic_meters; + 3 * mu + 9 / 2 * std::pow(sigma, 2)) * si::cubic_meters; //mass of the average ice A particle: real_t delta = IWCS / (ria * rhod_0); quantity m_a = delta * m_as + (real_t(1) - delta) * m_al; @@ -133,7 +133,7 @@ namespace libcloudphxx // mean velocity of ice A particle // eq. A.15b in Grabowski 1999 - template + template quantity::type, real_t> velocity_a( const quantity &ria, const quantity &rhod_0 @@ -141,30 +141,30 @@ namespace libcloudphxx using namespace common; //velocity of small ice A particle: quantity::type, real_t> v_as = - real_t(0.1) * si::meters / si::seconds; + real_t(0.1) * si::meters / si::seconds; //velocity of large ice A particle: quantity IWCS = std::min( - std::min(real_t(1e-3) * si::dimensionless(), - rhod_0 * ria / (si::kilogram / si::cubic_metres)), - real_t(2.52e-4) * std::pow( - real_t(1e3) * rhod_0 * ria / (si::kilogram / si::cubic_metres), - real_t(0.837)) * si::dimensionless()) * si::kilograms / - si::cubic_meters; + std::min(real_t(1e-3) * si::dimensionless(), + rhod_0 * ria / (si::kilogram / si::cubic_metres)), + real_t(2.52e-4) * std::pow( + real_t(1e3) * rhod_0 * ria / (si::kilogram / si::cubic_metres), + real_t(0.837)) * si::dimensionless()) * si::kilograms / + si::cubic_meters; quantity IWCL = rhod_0 * ria - IWCS; quantity::type, real_t> v_al = - (real_t(0.9) + real_t(0.1) * std::log10(real_t(1e3) * IWCL / si::kilograms * si::cubic_meters)) * si::meters - / si::seconds; + (real_t(0.9) + real_t(0.1) * std::log10(real_t(1e3) * IWCL / si::kilograms * si::cubic_meters)) * si::meters + / si::seconds; //average velocity of ice A particle: real_t delta = IWCS / (ria * rhod_0); quantity::type, real_t> v_a = - (delta * v_as + (real_t(1) - delta) * v_al) * std::pow( - real_t(0.3) / rhod_0 * si::kilograms / si::cubic_meters, real_t(0.5)); + (delta * v_as + (real_t(1) - delta) * v_al) * std::pow( + real_t(0.3) / rhod_0 * si::kilograms / si::cubic_meters, real_t(0.5)); return v_a; } //ice A homogeneous nucleation 1 (rv -> ria) // eq. A.21a in Grabowski 1999 - template + template quantity::type, real_t> hom_A_nucleation_1( const quantity &rv, const quantity &rvs, @@ -179,15 +179,15 @@ namespace libcloudphxx quantity rv_adj = beta * rvs + (1 - beta) * rvsi; //homogeneous nucleation rate (only if T < -40 C) quantity::type, real_t> hom1 = - (T < real_t(233.16) * si::kelvins) - ? (real_t(1) - std::exp(-dt / taunuc)) * std::max(real_t(0) * si::dimensionless(), rv - rv_adj) / dt - : real_t(0) / si::seconds; + (T < real_t(233.16) * si::kelvins) + ? (real_t(1) - std::exp(-dt / taunuc)) * std::max(real_t(0) * si::dimensionless(), rv - rv_adj) / dt + : real_t(0) / si::seconds; return hom1; } //ice A homogeneous nucleation 2 (rc -> ria) // eq. A.21b in Grabowski 1999 - template + template quantity::type, real_t> hom_A_nucleation_2( const quantity &rc, const quantity &T, @@ -196,15 +196,15 @@ namespace libcloudphxx const quantity taunuc = dt; // time scale for nucleation //homogeneous nucleation rate (only if T < -40 C) quantity::type, real_t> hom2 = - (T < real_t(233.16) * si::kelvins) - ? (real_t(1) - std::exp(-dt / taunuc)) * rc / dt - : real_t(0) / si::seconds; + (T < real_t(233.16) * si::kelvins) + ? (real_t(1) - std::exp(-dt / taunuc)) * rc / dt + : real_t(0) / si::seconds; return hom2; } //ice A heterogeneous nucleation (rc-> ria) // eq. A.19 in Grabowski 1999 - template + template quantity::type, real_t> het_A_nucleation( const quantity &ria, const quantity &rc, @@ -221,17 +221,17 @@ namespace libcloudphxx real_t(1e-2) * std::exp(real_t(0.6) * (real_t(273.16) - T / si::kelvins))) / si::cubic_meters; //nucleation rate: quantity::type, real_t> het = - (real_t(1) - std::exp(-dt / taunuc)) * std::min(rc, std::max(real_t(0) * si::dimensionless(), - N_in * std::max( - real_t(1e-12) * si::dimensionless(), - m_a / si::kilograms) / rhod_0 * si::kilograms - - ria)) / dt; + (real_t(1) - std::exp(-dt / taunuc)) * std::min(rc, std::max(real_t(0) * si::dimensionless(), + N_in * std::max( + real_t(1e-12) * si::dimensionless(), + m_a / si::kilograms) / rhod_0 * si::kilograms + - ria)) / dt; return het; } // ice B heterogeneous nucleation 1 (rr -> rib) // eq. A.23 in Grabowski 1999 - template + template quantity::type, real_t> het_B_nucleation_1( const quantity &rr, const quantity &ria, @@ -244,7 +244,7 @@ namespace libcloudphxx real_t gamma_r = std::pow(pi() * moist_air::rho_w() * N_0r / (rhod_0 * rr), real_t(0.25)); //mean raindrop mass quantity m_r = pi() * moist_air::rho_w() / (real_t(6) * std::pow(gamma_r, 3)) - * si::cubic_meters; + * si::cubic_meters; //mean raindrop velocity real_t v_r = real_t(251) * std::pow(gamma_r * rhod_0 / si::kilograms * si::cubic_meters, real_t(-0.5)); //mean raindrop radius @@ -257,13 +257,13 @@ namespace libcloudphxx real_t N_ra = N_0r / gamma_r * std::abs(v_r - v_a) * pi() * R_r * R_r * rhod_0 * ria / m_a * si::meters; //nucleation rate: quantity::type, real_t> het1 = - N_ra * m_r / si::kilograms / si::seconds; + N_ra * m_r / si::kilograms / si::seconds; return het1; } //ice B heterogeneous nucleation 2 (ria -> rib) // eq. A.23 in Grabowski 1999 - template + template quantity::type, real_t> het_B_nucleation_2( const quantity &rr, const quantity &ria, @@ -286,10 +286,32 @@ namespace libcloudphxx real_t N_ra = N_0r / gamma_r * abs(v_r - v_a) * pi() * R_r * R_r * rhod_0 * ria / m_a * si::meters; //nucleation rate: quantity::type, real_t> het2 = - N_ra * m_a / si::kilograms / si::seconds; + N_ra * m_a / si::kilograms / si::seconds; return het2; } + //melting of ice A (ria -> rr) + // eq. A.26 in Grabowski 1999 + template + quantity::type, real_t> melting_A( + const quantity &ria, + const quantity &T, + const quantity &rhod_0 + ) { + using namespace common; + //mass of average ice A particle: + quantity m_a = mass_a(ria, T, rhod_0); + //diameter of average ice A particle: + quantity D_a = std::pow(m_a / real_t(0.025) / si::kilograms, real_t(0.5)) * si::meters; + //fall velocity of avg. ice A particle: + quantity::type, real_t> v_a = real_t(4) * std::pow(D_a / si::meters, real_t(0.25)) *si::meters/si::seconds; + //Reynolds number: + real_t Re = D_a * v_a * rhod_0 / vterm::visc(T); + real_t F_a = real_t(0.78) + real_t(0.27) * std::pow(Re, real_t(0.5)); + quantity::type, real_t> mela = real_t(4.5e-7) * rhod_0 * ria / + m_a * D_a * (273.16 - T / si::kelvins) * F_a * si::square_meters / si::seconds; + return mela; + } }; }; }; diff --git a/include/libcloudph++/blk_1m/options.hpp b/include/libcloudph++/blk_1m/options.hpp index cbcedafb..d14eee6d 100644 --- a/include/libcloudph++/blk_1m/options.hpp +++ b/include/libcloudph++/blk_1m/options.hpp @@ -24,7 +24,9 @@ namespace libcloudphxx homA1 = true, // homogeneous nucleation 1 of ice A homA2 = true, // homogeneous nucleation 2 of ice A hetA = true, // heterogeneous nucleation of ice A - hetB = true; // heterogeneous nucleation of ice B + hetB = true, // heterogeneous nucleation of ice B + melA = true, // melting of ice A + melB = true; // melting of ice B real_t r_c0 = 5e-4, // autoconv. threshold k_acnv = 0.001, // Kessler autoconversion (eq. 5a in Grabowski & Smolarkiewicz 1996) diff --git a/include/libcloudph++/blk_1m/rhs_cellwise.hpp b/include/libcloudph++/blk_1m/rhs_cellwise.hpp index ca6f88cb..26e53624 100644 --- a/include/libcloudph++/blk_1m/rhs_cellwise.hpp +++ b/include/libcloudph++/blk_1m/rhs_cellwise.hpp @@ -159,6 +159,7 @@ namespace libcloudphxx rc_to_ria = 0, rr_to_rib = 0, ria_to_rib = 0, + ria_to_rr = 0, &dot_th = std::get<0>(tup), &dot_rv = std::get<1>(tup), &dot_rc = std::get<2>(tup), @@ -247,13 +248,25 @@ namespace libcloudphxx ); } + // melting of ice A + if (opts.melA) + { + ria_to_rr += ( + formulae::melting_A( + ria * si::dimensionless(), + T, + rhod + ) * si::seconds // to make it dimensionless + ); + } + dot_rc -= rc_to_ria; dot_rv -= rv_to_ria; - dot_rr -= rr_to_rib; - dot_ria += rc_to_ria + rv_to_ria - ria_to_rib; + dot_rr += ria_to_rr - rr_to_rib; + dot_ria += rc_to_ria + rv_to_ria - ria_to_rib - ria_to_rr; dot_rib += rr_to_rib + ria_to_rib; dot_th += const_cp::l_s(T) / (moist_air::c_pd() * theta_std::exner(p)) * rv_to_ria / si::kelvins; //heat of sublimation - dot_th += const_cp::l_f(T) / (moist_air::c_pd() * theta_std::exner(p)) * (rc_to_ria+rr_to_rib) / si::kelvins; //heat of freezing + dot_th += const_cp::l_f(T) / (moist_air::c_pd() * theta_std::exner(p)) * (rc_to_ria+rr_to_rib-ria_to_rr) / si::kelvins; //heat of freezing } } diff --git a/tests/python/unit/api_blk_1m.py b/tests/python/unit/api_blk_1m.py index 917fd3f1..a57e73af 100644 --- a/tests/python/unit/api_blk_1m.py +++ b/tests/python/unit/api_blk_1m.py @@ -16,6 +16,8 @@ print("homA2 =", opts.homA2) print("hetA =", opts.hetA) print("hetB =", opts.hetB) +print("melA =", opts.melA) +print("melB =", opts.melB) print("r_c0 =", opts.r_c0) print("r_eps =", opts.r_eps) From 30f1ccae73e02abf303a87ca5b1f172d2542363a Mon Sep 17 00:00:00 2001 From: Agnieszka Makulska Date: Wed, 6 Nov 2024 13:10:45 +0100 Subject: [PATCH 10/17] melting of ice B and limiting the changes of variables --- bindings/python/lib.cpp | 4 ++ include/libcloudph++/blk_1m/formulae.hpp | 52 +++++++++++++++++++- include/libcloudph++/blk_1m/options.hpp | 14 ++++-- include/libcloudph++/blk_1m/rhs_cellwise.hpp | 18 +++++++ tests/python/unit/api_blk_1m.py | 6 ++- 5 files changed, 86 insertions(+), 8 deletions(-) diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index b363a184..6f3253f7 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -159,6 +159,10 @@ BOOST_PYTHON_MODULE(libcloudphxx) .def_readwrite("homA2", &b1m::opts_t::homA2) .def_readwrite("hetA", &b1m::opts_t::hetA) .def_readwrite("hetB", &b1m::opts_t::hetB) + .def_readwrite("depA", &b1m::opts_t::depA) + .def_readwrite("depB", &b1m::opts_t::depB) + .def_readwrite("rimA", &b1m::opts_t::rimA) + .def_readwrite("rimB", &b1m::opts_t::rimB) .def_readwrite("melA", &b1m::opts_t::melA) .def_readwrite("melB", &b1m::opts_t::melB) .def_readwrite("r_c0", &b1m::opts_t::r_c0) diff --git a/include/libcloudph++/blk_1m/formulae.hpp b/include/libcloudph++/blk_1m/formulae.hpp index b941ce54..1b07c948 100644 --- a/include/libcloudph++/blk_1m/formulae.hpp +++ b/include/libcloudph++/blk_1m/formulae.hpp @@ -239,6 +239,10 @@ namespace libcloudphxx const quantity &rhod_0 ) { using namespace common; + if (ria == 0) // ice B is formed only if ice A is present + { + return real_t(0) / si::seconds; + } //parameters of Marshall & Palmer distribution for rain: real_t N_0r = real_t(1e7); real_t gamma_r = std::pow(pi() * moist_air::rho_w() * N_0r / (rhod_0 * rr), real_t(0.25)); @@ -271,6 +275,10 @@ namespace libcloudphxx const quantity &rhod_0 ) { using namespace common; + if (ria == 0) // ice B is formed only if ice A is present + { + return real_t(0) / si::seconds; + } //parameters of Marshall & Palmer distribution for rain: real_t N_0r = real_t(1e7); real_t gamma_r = std::pow(pi() * moist_air::rho_w() * N_0r / (rhod_0 * rr), real_t(0.25)); @@ -299,19 +307,61 @@ namespace libcloudphxx const quantity &rhod_0 ) { using namespace common; + if (ria == 0) // calculating melting rate only if ice A is present + { + return real_t(0) / si::seconds; + } //mass of average ice A particle: quantity m_a = mass_a(ria, T, rhod_0); //diameter of average ice A particle: quantity D_a = std::pow(m_a / real_t(0.025) / si::kilograms, real_t(0.5)) * si::meters; //fall velocity of avg. ice A particle: - quantity::type, real_t> v_a = real_t(4) * std::pow(D_a / si::meters, real_t(0.25)) *si::meters/si::seconds; + quantity::type, real_t> v_a = real_t(4) * std::pow( + D_a / si::meters, real_t(0.25)) * si::meters / si::seconds; //Reynolds number: real_t Re = D_a * v_a * rhod_0 / vterm::visc(T); + //ventilation factor: real_t F_a = real_t(0.78) + real_t(0.27) * std::pow(Re, real_t(0.5)); quantity::type, real_t> mela = real_t(4.5e-7) * rhod_0 * ria / m_a * D_a * (273.16 - T / si::kelvins) * F_a * si::square_meters / si::seconds; return mela; } + + //melting of ice B (rib -> rr) + // eq. A.26 in Grabowski 1999 + template + quantity::type, real_t> melting_B( + const quantity &rib, + const quantity &T, + const quantity &rhod_0 + ) { + using namespace common; + if (rib == 0) // calculating melting rate only if ice B is present + { + return real_t(0) / si::seconds; + } + quantity rho_b = real_t(400) * si::kilograms / si::cubic_meters; + //graupel density (from Grabowski 1999) + // ice B distribution parameters + real_t N_0b = real_t(4e6); + real_t gamma_b = std::pow(pi() * rho_b * N_0b / (rhod_0 * rib), real_t(0.25)); + //mass of average ice B particle: + quantity m_b = pi() * rho_b / (real_t(6) * std::pow(gamma_b, real_t(3))) * + si::cubic_meters; + //diameter of average ice B particle: + quantity D_b = real_t(1) / gamma_b * si::meters; + //fall velocity of avg. ice B particle: + quantity::type, real_t> v_b = real_t(31.2) * + std::pow(gamma_b, real_t(-0.37)) * std::pow(rhod_0 / si::kilograms * si::cubic_meters, real_t(-0.5)) * + si::meters / si::seconds; + //Reynolds number: + real_t Re = D_b * v_b * rhod_0 / vterm::visc(T); + //ventilation factor: + real_t F_b = real_t(0.78) + real_t(0.27) * std::pow(Re, real_t(0.5)); + quantity::type, real_t> melb = real_t(4.5e-7) * rhod_0 * rib / + m_b * D_b * (273.16 - T / si::kelvins) * F_b * si::square_meters / si::seconds; + return melb; + } }; }; }; diff --git a/include/libcloudph++/blk_1m/options.hpp b/include/libcloudph++/blk_1m/options.hpp index d14eee6d..553e7510 100644 --- a/include/libcloudph++/blk_1m/options.hpp +++ b/include/libcloudph++/blk_1m/options.hpp @@ -21,12 +21,16 @@ namespace libcloudphxx conv = true, // autoconversion accr = true, // accretion sedi = true, // sedimentation - homA1 = true, // homogeneous nucleation 1 of ice A - homA2 = true, // homogeneous nucleation 2 of ice A + homA1 = true, // homogeneous nucleation of ice A from water vapor + homA2 = true, // homogeneous nucleation of ice A from cloud droplets hetA = true, // heterogeneous nucleation of ice A - hetB = true, // heterogeneous nucleation of ice B - melA = true, // melting of ice A - melB = true; // melting of ice B + hetB = true, // heterogeneous nucleation of ice B + depA = true, // depositional growth of ice A + depB = true, // depositional growth of ice B + rimA = true, // growth of ice A by riming + rimB = true, // growth of ice B by riming + melA = true, // melting of ice A + melB = true; // melting of ice B real_t r_c0 = 5e-4, // autoconv. threshold k_acnv = 0.001, // Kessler autoconversion (eq. 5a in Grabowski & Smolarkiewicz 1996) diff --git a/include/libcloudph++/blk_1m/rhs_cellwise.hpp b/include/libcloudph++/blk_1m/rhs_cellwise.hpp index 26e53624..b9033783 100644 --- a/include/libcloudph++/blk_1m/rhs_cellwise.hpp +++ b/include/libcloudph++/blk_1m/rhs_cellwise.hpp @@ -260,6 +260,24 @@ namespace libcloudphxx ); } + // melting of ice B + if (opts.melB) + { + rr_to_rib -= ( + formulae::melting_B( + rib * si::dimensionless(), + T, + rhod + ) * si::seconds // to make it dimensionless + ); + } + //limiting + rv_to_ria = std::min(rv, rv_to_ria) / dt; + rc_to_ria = std::min(rc, rc_to_ria) / dt; + rr_to_rib = std::min(rr, rr_to_rib) / dt; + ria_to_rib = std::min(ria, ria_to_rib) / dt; + ria_to_rr = std::min(ria, ria_to_rr) / dt; + dot_rc -= rc_to_ria; dot_rv -= rv_to_ria; dot_rr += ria_to_rr - rr_to_rib; diff --git a/tests/python/unit/api_blk_1m.py b/tests/python/unit/api_blk_1m.py index a57e73af..c34b71d3 100644 --- a/tests/python/unit/api_blk_1m.py +++ b/tests/python/unit/api_blk_1m.py @@ -16,6 +16,10 @@ print("homA2 =", opts.homA2) print("hetA =", opts.hetA) print("hetB =", opts.hetB) +print("depA =", opts.depA) +print("depB =", opts.depB) +print("rimA =", opts.rimA) +print("rimB =", opts.rimB) print("melA =", opts.melA) print("melB =", opts.melB) print("r_c0 =", opts.r_c0) @@ -27,8 +31,6 @@ rc = arr_t([0.01]) rv = arr_t([0. ]) rr = arr_t([0. ]) -ria = arr_t([0. ]) -rib = arr_t([0. ]) dt = 1 dz = 1 From e9270f45e49575d2e05f9842cd03f96e256b3727 Mon Sep 17 00:00:00 2001 From: Agnieszka Makulska Date: Tue, 12 Nov 2024 13:27:37 +0100 Subject: [PATCH 11/17] growth by deposition and riming --- include/libcloudph++/blk_1m/formulae.hpp | 332 +++++++++++++++++-- include/libcloudph++/blk_1m/rhs_cellwise.hpp | 90 ++++- include/libcloudph++/common/moist_air.hpp | 2 + tests/python/unit/api_blk_1m.py | 2 + 4 files changed, 389 insertions(+), 37 deletions(-) diff --git a/include/libcloudph++/blk_1m/formulae.hpp b/include/libcloudph++/blk_1m/formulae.hpp index 1b07c948..0983de26 100644 --- a/include/libcloudph++/blk_1m/formulae.hpp +++ b/include/libcloudph++/blk_1m/formulae.hpp @@ -1,7 +1,7 @@ /** @file * @copyright University of Warsaw * @brief single-moment bulk parameterisation formulae (Kessler) - * from @copybrief bib::Grabowski_and_Smolarkiewicz_1996 + * from @copybrief bib::Grabowski_and_Smolarkiewicz_1996, Grabowski 1999 * @section LICENSE * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/) */ @@ -162,6 +162,93 @@ namespace libcloudphxx return v_a; } + // slope of Marshall Palmer distribution for rain + // eq. A.1 in Grabowski 1999 + template + quantity::type, real_t> lambda_rain( + const quantity &rr, + const quantity &rhod_0 + ) { + using namespace common; + real_t N_0r = real_t(1e7); + return std::pow(pi() * moist_air::rho_w() * N_0r / (rhod_0 * rr), real_t(0.25)) / si::meters; + } + + // slope of Marshall Palmer distribution for ice B + // eq. A.4 in Grabowski 1999 + template + quantity::type, real_t> lambda_ice_b( + const quantity &rib, + const quantity &rhod_0 + ) { + using namespace common; + real_t N_0b = real_t(4e6); + return std::pow(pi() * moist_air::rho_ib() * N_0b / (rhod_0 * rib), real_t(0.25)) / si::meters; + } + + + // mean mass of ice B particle + // eq. A.5 in Grabowski 1999 + template + quantity mass_b( + const quantity &rib, + const quantity &rhod_0 + ) { + using namespace common; + return pi() * moist_air::rho_ib() / (real_t(6) * std::pow( + lambda_ice_b(rib, rhod_0) * si::meters, real_t(3))) * si::cubic_meters; + } + + //coefficients for ice growth by deposition/riming + // table 2 from Koenig 1972 + template + real_t coeff_alpha(const quantity &T) { + real_t alpha_list[32] = { + real_t(0.), real_t(0.7939e-7), real_t(0.7841e-6), real_t(0.3369e-5), real_t(0.4336e-5), real_t(0.5285e-5), + real_t(0.3728e-5), real_t(0.1852e-5), real_t(0.2991e-6), real_t(0.4248e-6), real_t(0.7434e-6), + real_t(0.1812e-5), + real_t(0.4394e-5), real_t(0.9145e-5), real_t(0.1725e-4), real_t(0.3348e-4), real_t(0.1725e-4), + real_t(0.9175e-5), + real_t(0.4412e-5), real_t(0.2252e-5), real_t(0.9115e-6), real_t(0.4876e-6), real_t(0.3473e-6), + real_t(0.4758e-6), + real_t(0.6306e-6), real_t(0.8573e-6), real_t(0.7868e-6), real_t(0.7192e-6), real_t(0.6515e-6), + real_t(0.5956e-6), + real_t(0.533e-6), real_t(0.4834e-6) + }; + //temperature in Celsius: + real_t Tc = T / si::kelvins - real_t(273.16); + //limiting the temerature to (-31, 0): + real_t ttcoe = std::min(real_t(-0.), std::max(real_t(-31), Tc)); + int index1 = static_cast(std::trunc(-ttcoe)); //index between 0 and 31 + int index2 = index1 + 1; // index between 1 and 32 + real_t del = -ttcoe - real_t(index1); + return (real_t(1.) - del) * alpha_list[index1] + del * alpha_list[index2]; //interpolation + } + + //coefficients for ice growth by deposition/riming + // table 2 from Koenig 1972 + template + real_t coeff_beta(const quantity &T) { + real_t beta_list[32] = { + real_t(0.), real_t(0.4006), real_t(0.4831), real_t(0.5320), real_t(0.5307), real_t(0.5319), real_t(0.5249), + real_t(0.4888), + real_t(0.3894), real_t(0.4047), real_t(0.4318), real_t(0.4771), real_t(0.5183), real_t(0.5463), + real_t(0.5651), real_t(0.5813), + real_t(0.5655), real_t(0.5478), real_t(0.5203), real_t(0.4906), real_t(0.4447), real_t(0.4126), + real_t(0.3960), real_t(0.4149), + real_t(0.4320), real_t(0.4506), real_t(0.4483), real_t(0.4460), real_t(0.4433), real_t(0.4413), + real_t(0.4382), real_t(0.4361) + }; + //temperature in Celsius: + real_t Tc = T / si::kelvins - real_t(273.16); + //limiting the temerature to (-31, 0): + real_t ttcoe = std::min(real_t(-0.), std::max(real_t(-31), Tc)); + int index1 = static_cast(std::trunc(-ttcoe)); //index between 0 and 31 + int index2 = index1 + 1; // index between 1 and 32 + real_t del = -ttcoe - real_t(index1); + return (real_t(1.) - del) * beta_list[index1] + del * beta_list[index2]; //interpolation + } + //ice A homogeneous nucleation 1 (rv -> ria) // eq. A.21a in Grabowski 1999 template @@ -223,8 +310,8 @@ namespace libcloudphxx quantity::type, real_t> het = (real_t(1) - std::exp(-dt / taunuc)) * std::min(rc, std::max(real_t(0) * si::dimensionless(), N_in * std::max( - real_t(1e-12) * si::dimensionless(), - m_a / si::kilograms) / rhod_0 * si::kilograms + real_t(1e-12) * si::kilograms, + m_a) / rhod_0 - ria)) / dt; return het; } @@ -243,25 +330,25 @@ namespace libcloudphxx { return real_t(0) / si::seconds; } - //parameters of Marshall & Palmer distribution for rain: - real_t N_0r = real_t(1e7); - real_t gamma_r = std::pow(pi() * moist_air::rho_w() * N_0r / (rhod_0 * rr), real_t(0.25)); + quantity::type, real_t> lambda_r = lambda_rain(rr, rhod_0); //mean raindrop mass - quantity m_r = pi() * moist_air::rho_w() / (real_t(6) * std::pow(gamma_r, 3)) - * si::cubic_meters; + quantity m_r = pi() * moist_air::rho_w() / (real_t(6) * std::pow( + lambda_r * si::meters, 3)) * si::cubic_meters; //mean raindrop velocity - real_t v_r = real_t(251) * std::pow(gamma_r * rhod_0 / si::kilograms * si::cubic_meters, real_t(-0.5)); + real_t v_r = real_t(251) * std::pow(lambda_r * si::meters * rhod_0 / si::kilograms * si::cubic_meters, + real_t(-0.5)); //mean raindrop radius - quantity R_r = real_t(0.5) / gamma_r * si::meters; + quantity R_r = real_t(0.5) / lambda_r; //mean ice A particle mass quantity m_a = mass_a(ria, T, rhod_0); //mean ice A velocity real_t v_a = velocity_a(ria, rhod_0) * si::seconds / si::meters; //rate of collisions between raindrops and ice A - real_t N_ra = N_0r / gamma_r * std::abs(v_r - v_a) * pi() * R_r * R_r * rhod_0 * ria / m_a * si::meters; + quantity::type, real_t> N_ra = real_t(1e7) / si::cubic_meters / + si::meters / lambda_r * std::abs(v_r - v_a) * si::meters / si::seconds * pi() * R_r * R_r * ria / m_a; //nucleation rate: quantity::type, real_t> het1 = - N_ra * m_r / si::kilograms / si::seconds; + N_ra * m_r; return het1; } @@ -279,22 +366,22 @@ namespace libcloudphxx { return real_t(0) / si::seconds; } - //parameters of Marshall & Palmer distribution for rain: - real_t N_0r = real_t(1e7); - real_t gamma_r = std::pow(pi() * moist_air::rho_w() * N_0r / (rhod_0 * rr), real_t(0.25)); + quantity::type, real_t> lambda_r = lambda_rain(rr, rhod_0); //mean raindrop velocity - real_t v_r = real_t(251) * std::pow(gamma_r * rhod_0 / si::kilograms * si::cubic_meters, real_t(-0.5)); + real_t v_r = real_t(251) * std::pow(lambda_r * si::meters * rhod_0 / si::kilograms * si::cubic_meters, + real_t(-0.5)); //mean raindrop radius - quantity R_r = real_t(0.5) / gamma_r * si::meters; + quantity R_r = real_t(0.5) / lambda_r; //mean ice A particle mass quantity m_a = mass_a(ria, T, rhod_0); //mean ice A velocity real_t v_a = velocity_a(ria, rhod_0) * si::seconds / si::meters; //rate of collisions between raindrops and ice A - real_t N_ra = N_0r / gamma_r * abs(v_r - v_a) * pi() * R_r * R_r * rhod_0 * ria / m_a * si::meters; + quantity::type, real_t> N_ra = real_t(1e7) / si::cubic_meters / + si::meters / lambda_r * abs(v_r - v_a) * si::meters / si::seconds * pi() * R_r * R_r * ria / m_a; //nucleation rate: quantity::type, real_t> het2 = - N_ra * m_a / si::kilograms / si::seconds; + N_ra * m_a; return het2; } @@ -322,8 +409,8 @@ namespace libcloudphxx real_t Re = D_a * v_a * rhod_0 / vterm::visc(T); //ventilation factor: real_t F_a = real_t(0.78) + real_t(0.27) * std::pow(Re, real_t(0.5)); - quantity::type, real_t> mela = real_t(4.5e-7) * rhod_0 * ria / - m_a * D_a * (273.16 - T / si::kelvins) * F_a * si::square_meters / si::seconds; + quantity::type, real_t> mela = real_t(4.5e-7) * ria / + m_a * D_a * (273.16 - T / si::kelvins) * F_a * si::kilograms / si::meters / si::seconds; return mela; } @@ -340,28 +427,209 @@ namespace libcloudphxx { return real_t(0) / si::seconds; } - quantity rho_b = real_t(400) * si::kilograms / si::cubic_meters; - //graupel density (from Grabowski 1999) - // ice B distribution parameters - real_t N_0b = real_t(4e6); - real_t gamma_b = std::pow(pi() * rho_b * N_0b / (rhod_0 * rib), real_t(0.25)); + quantity::type, real_t> lambda_b = + lambda_ice_b(rib, rhod_0); //mass of average ice B particle: - quantity m_b = pi() * rho_b / (real_t(6) * std::pow(gamma_b, real_t(3))) * - si::cubic_meters; + quantity m_b = mass_b(rib, rhod_0); //diameter of average ice B particle: - quantity D_b = real_t(1) / gamma_b * si::meters; + quantity D_b = real_t(1) / lambda_b; //fall velocity of avg. ice B particle: quantity::type, real_t> v_b = real_t(31.2) * - std::pow(gamma_b, real_t(-0.37)) * std::pow(rhod_0 / si::kilograms * si::cubic_meters, real_t(-0.5)) * + std::pow(lambda_b * si::meters, real_t(-0.37)) * std::pow(rhod_0 / si::kilograms * si::cubic_meters, + real_t(-0.5)) * si::meters / si::seconds; //Reynolds number: real_t Re = D_b * v_b * rhod_0 / vterm::visc(T); //ventilation factor: real_t F_b = real_t(0.78) + real_t(0.27) * std::pow(Re, real_t(0.5)); - quantity::type, real_t> melb = real_t(4.5e-7) * rhod_0 * rib / - m_b * D_b * (273.16 - T / si::kelvins) * F_b * si::square_meters / si::seconds; + quantity::type, real_t> melb = real_t(4.5e-7) * rib / + m_b * D_b * (273.16 - T / si::kelvins) * F_b * si::kilograms / si::meters / si::seconds; return melb; } + + + // growth of ice A by deposition (rv->ria) + // eq. A.24a in Grabowski 1999, eq. 27 in Koenig 1976 + template + quantity::type, real_t> deposition_A( + const quantity &ria, + const quantity &rv, + const quantity &rvs, + const quantity &rvsi, + const quantity &T, + const quantity &rhod_0 + ) { + if (ria == 0 || T > 273.16 * si::kelvins) + { + return real_t(0) / si::seconds; + } + //mass of average ice A particle: + quantity m_a = mass_a(ria, T, rhod_0); + real_t alpha = coeff_alpha(T); + real_t beta = coeff_beta(T); + // growth rate for a single particle: + quantity::type, real_t> dm_dt_AE = (rv - rvsi) / (rvs - rvsi) * alpha * + std::pow(m_a / si::kilograms, beta) * si::kilograms / si::seconds; + quantity::type, real_t> depa = ria / m_a * dm_dt_AE; + return depa; + } + + + //growth of ice A by riming (rc->ria) + // eq. A.24b in Grabowski 1999, eq. 27-34 in Koenig 1976 + template + quantity::type, real_t> riming_A( + const quantity &ria, + const quantity &rc, + const quantity &rv, + const quantity &rvs, + const quantity &rvsi, + const quantity &T, + const quantity &rhod_0 + ) { + if (ria == 0 || T > 273.16 * si::kelvins) + { + return real_t(0) / si::seconds; + } + //mass of average ice A particle: + quantity m_a = mass_a(ria, T, rhod_0); + real_t alpha = coeff_alpha(T); + real_t beta = coeff_beta(T); + // regime AE + quantity::type, real_t> dm_dt_AE = (rv - rvsi) / (rvs - rvsi) * alpha * + std::pow(m_a / si::kilograms, beta) * si::kilograms / si::seconds; + // regime BC + real_t tan_theta = real_t(1.) + real_t(0.1) * std::log( + rhod_0 * rc * real_t(1e3) / si::kilograms * si::cubic_meters); + real_t gamma = alpha * std::pow(real_t(5e-8), beta); + quantity::type, real_t> dm_dt_BC = gamma * std::pow( + m_a / real_t(5e-11) / si::kilograms, tan_theta) * si::kilograms / si::seconds; + //regime CD + real_t dzeta = gamma * std::pow(real_t(2e3), tan_theta); + real_t xi = std::log(rc * rhod_0 / si::kilograms * si::cubic_meters * real_t(1e9) / dzeta) / std::log( + real_t(1e4)); + quantity::type, real_t> dm_dt_CD = dzeta * std::pow( + m_a * real_t(1e7) / si::kilograms, xi) * si::kilograms / si::seconds; + //total growth + quantity::type, real_t> rima = real_t(0) / si::seconds; + if (m_a > real_t(5e-11) * si::kilograms && m_a < real_t(1e-7) * si::kilograms) + { + rima += std::max(real_t(0) * si::kilograms / si::seconds, dm_dt_BC - dm_dt_AE) * ria / m_a; + } + if (m_a >= real_t(1e-7) * si::kilograms) + { + rima += std::max(real_t(0) * si::kilograms / si::seconds, dm_dt_CD - dm_dt_AE) * ria / m_a; + } + return rima; + } + + // growth of ice B by deposition (rv->rib) + // eq. A.24c in Grabowski 1999, eq. 27 in Koenig 1976 + template + quantity::type, real_t> deposition_B( + const quantity &rib, + const quantity &rv, + const quantity &rvs, + const quantity &rvsi, + const quantity &T, + const quantity &rhod_0 + ) { + if (rib == 0 || T > 273.16 * si::kelvins) // calculating growth rate only if ice A is present + { + return real_t(0) / si::seconds; + } + //mass of average ice B particle: + quantity m_b = mass_b(rib, rhod_0); + real_t alpha = coeff_alpha(T); + real_t beta = coeff_beta(T); + // growth rate for a single particle: + quantity::type, real_t> dm_dt_AE = (rv - rvsi) / (rvs - rvsi) * alpha * + std::pow(m_b / si::kilograms, beta) * si::kilograms / si::seconds; + quantity::type, real_t> depb = rib / m_b * dm_dt_AE; + return depb; + } + + //growth of ice B by riming (rc, rr ->rib) + // eq. A.24d in Grabowski 1999, eq. 27-34 in Koenig 1976 + template + quantity::type, real_t> riming_B( + const quantity &rib, + const quantity &rc, + const quantity &rv, + const quantity &rvs, + const quantity &rvsi, + const quantity &T, + const quantity &rhod_0 + ) { + if (rib == 0 || T > 273.16 * si::kelvins) + { + return real_t(0) / si::seconds; + } + //mass of average ice B particle: + quantity m_b = mass_a(rib, T, rhod_0); + real_t alpha = coeff_alpha(T); + real_t beta = coeff_beta(T); + // regime AE + quantity::type, real_t> dm_dt_AE = (rv - rvsi) / (rvs - rvsi) * alpha * + std::pow(m_b / si::kilograms, beta) * si::kilograms / si::seconds; + // regime BC + real_t tan_theta = real_t(1.) + real_t(0.1) * std::log( + rhod_0 * rc * real_t(1e3) / si::kilograms * si::cubic_meters); + real_t gamma = alpha * std::pow(real_t(5e-8), beta); + quantity::type, real_t> dm_dt_BC = gamma * std::pow( + m_b / real_t(5e-11) / si::kilograms, tan_theta) * si::kilograms / si::seconds; + //regime CD + real_t dzeta = gamma * std::pow(real_t(2e3), tan_theta); + real_t xi = std::log(rc * rhod_0 / si::kilograms * si::cubic_meters * real_t(1e9) / dzeta) / std::log( + real_t(1e4)); + quantity::type, real_t> dm_dt_CD = dzeta * std::pow( + m_b * real_t(1e7) / si::kilograms, xi) * si::kilograms / si::seconds; + //total growth + quantity::type, real_t> rimb = real_t(0) / si::seconds; + if (m_b > real_t(5e-11) * si::kilograms && m_b < real_t(1e-7) * si::kilograms) + { + rimb += std::max(real_t(0) * si::kilograms / si::seconds, dm_dt_BC - dm_dt_AE) * rib / m_b; + } + if (m_b >= real_t(1e-7) * si::kilograms) + { + rimb += std::max(real_t(0) * si::kilograms / si::seconds, dm_dt_CD - dm_dt_AE) * rib / m_b; + } + return rimb; + } + + //growth of ice B by riming (only rc->rib) + // eq. A.24d in Grabowski 1999, eq. 27-34 in Koenig 1976 + template + quantity::type, real_t> riming_B_1( + const quantity &rib, + const quantity &rc, + const quantity &rr, + const quantity &rv, + const quantity &rvs, + const quantity &rvsi, + const quantity &T, + const quantity &rhod_0 + ) { + real_t coeff_rc = rc / (rc + rr + real_t(1e-10)); + return coeff_rc * riming_B(rib, rc, rv, rvs, rvsi, T, rhod_0); + } + + //growth of ice B by riming (only rr->rib) + // eq. A.24d in Grabowski 1999, eq. 27-34 in Koenig 1976 + template + quantity::type, real_t> riming_B_2( + const quantity &rib, + const quantity &rc, + const quantity &rr, + const quantity &rv, + const quantity &rvs, + const quantity &rvsi, + const quantity &T, + const quantity &rhod_0 + ) { + real_t coeff_rc = rc / (rc + rr + real_t(1e-10)); + return (real_t(1) - coeff_rc) * riming_B(rib, rc, rv, rvs, rvsi, T, rhod_0); + } }; }; }; diff --git a/include/libcloudph++/blk_1m/rhs_cellwise.hpp b/include/libcloudph++/blk_1m/rhs_cellwise.hpp index b9033783..e0628597 100644 --- a/include/libcloudph++/blk_1m/rhs_cellwise.hpp +++ b/include/libcloudph++/blk_1m/rhs_cellwise.hpp @@ -156,7 +156,9 @@ namespace libcloudphxx real_t rv_to_ria = 0, + rv_to_rib = 0, rc_to_ria = 0, + rc_to_rib = 0, rr_to_rib = 0, ria_to_rib = 0, ria_to_rr = 0, @@ -271,20 +273,98 @@ namespace libcloudphxx ) * si::seconds // to make it dimensionless ); } + + // depositional growth of ice A + if (opts.depA) + { + rv_to_ria += ( + formulae::deposition_A( + ria * si::dimensionless(), + rv * si::dimensionless(), + rvs * si::dimensionless(), + rvsi * si::dimensionless(), + T, + rhod + ) * si::seconds // to make it dimensionless + ); + } + + // growth of ice A by riming + if (opts.rimA) + { + rc_to_ria += ( + formulae::riming_A( + ria * si::dimensionless(), + rc * si::dimensionless(), + rv * si::dimensionless(), + rvs * si::dimensionless(), + rvsi * si::dimensionless(), + T, + rhod + ) * si::seconds // to make it dimensionless + ); + } + + // depositional growth of ice B + if (opts.depB) + { + rv_to_rib += ( + formulae::deposition_B( + rib * si::dimensionless(), + rv * si::dimensionless(), + rvs * si::dimensionless(), + rvsi * si::dimensionless(), + T, + rhod + ) * si::seconds // to make it dimensionless + ); + } + + // growth of ice B by riming + if (opts.rimB) + { + rc_to_rib += ( + formulae::riming_B_1( + rib * si::dimensionless(), + rc * si::dimensionless(), + rr * si::dimensionless(), + rv * si::dimensionless(), + rvs * si::dimensionless(), + rvsi * si::dimensionless(), + T, + rhod + ) * si::seconds // to make it dimensionless + ); + rr_to_rib += ( + formulae::riming_B_2( + rib * si::dimensionless(), + rc * si::dimensionless(), + rr * si::dimensionless(), + rv * si::dimensionless(), + rvs * si::dimensionless(), + rvsi * si::dimensionless(), + T, + rhod + ) * si::seconds // to make it dimensionless + ); + } + //limiting rv_to_ria = std::min(rv, rv_to_ria) / dt; + rv_to_rib = std::min(rv, rv_to_rib) / dt; rc_to_ria = std::min(rc, rc_to_ria) / dt; + rc_to_rib = std::min(rc, rc_to_rib) / dt; rr_to_rib = std::min(rr, rr_to_rib) / dt; ria_to_rib = std::min(ria, ria_to_rib) / dt; ria_to_rr = std::min(ria, ria_to_rr) / dt; - dot_rc -= rc_to_ria; - dot_rv -= rv_to_ria; + dot_rc += - rc_to_ria - rc_to_rib; + dot_rv += - rv_to_ria - rv_to_rib; dot_rr += ria_to_rr - rr_to_rib; dot_ria += rc_to_ria + rv_to_ria - ria_to_rib - ria_to_rr; - dot_rib += rr_to_rib + ria_to_rib; - dot_th += const_cp::l_s(T) / (moist_air::c_pd() * theta_std::exner(p)) * rv_to_ria / si::kelvins; //heat of sublimation - dot_th += const_cp::l_f(T) / (moist_air::c_pd() * theta_std::exner(p)) * (rc_to_ria+rr_to_rib-ria_to_rr) / si::kelvins; //heat of freezing + dot_rib += rr_to_rib + ria_to_rib + rv_to_rib + rc_to_rib; + dot_th += const_cp::l_s(T) / (moist_air::c_pd() * theta_std::exner(p)) * (rv_to_ria + rv_to_rib) / si::kelvins; //heat of sublimation + dot_th += const_cp::l_f(T) / (moist_air::c_pd() * theta_std::exner(p)) * (rc_to_ria + rc_to_rib + rr_to_rib - ria_to_rr) / si::kelvins; //heat of freezing } } diff --git a/include/libcloudph++/common/moist_air.hpp b/include/libcloudph++/common/moist_air.hpp index 159b04cf..74c5a318 100644 --- a/include/libcloudph++/common/moist_air.hpp +++ b/include/libcloudph++/common/moist_air.hpp @@ -50,6 +50,8 @@ namespace libcloudphxx libcloudphxx_const(si::mass_density, rho_w, 1e3, si::kilograms / si::cubic_metres) // ice density libcloudphxx_const(si::mass_density, rho_i, 910, si::kilograms / si::cubic_metres) + // graupel density used for ice B (from Grabowski 1999) + libcloudphxx_const(si::mass_density, rho_ib, 400, si::kilograms / si::cubic_metres) // mixing rule for extensive quantitites (i.e. using mass mixing ratio) template diff --git a/tests/python/unit/api_blk_1m.py b/tests/python/unit/api_blk_1m.py index c34b71d3..7c6af568 100644 --- a/tests/python/unit/api_blk_1m.py +++ b/tests/python/unit/api_blk_1m.py @@ -92,6 +92,8 @@ assert dot_rr == dot_rr_old # no rain water -> no precip th = arr_t([230.]) #testing ice physics +ria = arr_t([0.1]) +rib = arr_t([0.1]) dot_rc = arr_t([0.]) dot_rr = arr_t([0.]) dot_rv = arr_t([0.]) From ffebf8aa9aa5be5b7371292c998c35a263a1a2d2 Mon Sep 17 00:00:00 2001 From: Agnieszka Makulska Date: Tue, 12 Nov 2024 15:41:11 +0100 Subject: [PATCH 12/17] sedimentation of ice --- bindings/python/blk_1m.hpp | 16 ++- include/libcloudph++/blk_1m/formulae.hpp | 129 +++++++++--------- .../libcloudph++/blk_1m/rhs_columnwise.hpp | 113 ++++++++++----- include/libcloudph++/common/moist_air.hpp | 2 - tests/python/unit/api_blk_1m.py | 10 +- 5 files changed, 162 insertions(+), 108 deletions(-) diff --git a/bindings/python/blk_1m.hpp b/bindings/python/blk_1m.hpp index 64b907aa..ba912128 100644 --- a/bindings/python/blk_1m.hpp +++ b/bindings/python/blk_1m.hpp @@ -201,19 +201,21 @@ namespace libcloudphxx { template typename arr_t::T_numtype rhs_columnwise( const b1m::opts_t &opts, - bp_array &dot_rr, + bp_array &dot_r, const bp_array &rhod, - const bp_array &rr, - const typename arr_t::T_numtype &dz + const bp_array &r, + const typename arr_t::T_numtype &dz, + const std::string& precip_type = "rain" ) { arr_t - np2bz_dot_rr(np2bz(dot_rr)); + np2bz_dot_r(np2bz(dot_r)); return b1m::rhs_columnwise( opts, - np2bz_dot_rr, + np2bz_dot_r, np2bz(rhod), - np2bz(rr), - dz + np2bz(r), + dz, + precip_type ); } }; diff --git a/include/libcloudph++/blk_1m/formulae.hpp b/include/libcloudph++/blk_1m/formulae.hpp index 0983de26..cd07c282 100644 --- a/include/libcloudph++/blk_1m/formulae.hpp +++ b/include/libcloudph++/blk_1m/formulae.hpp @@ -41,7 +41,7 @@ namespace libcloudphxx return k_2() * rc * std::pow(rr, real_t(.875)); } - // Kessler evaporation rate + // Kessler rain evaporation rate // eq. 5c in Grabowski & Smolarkiewicz 1996 (multiplied by rho!) template quantity::type, real_t> evaporation_rate( @@ -70,7 +70,7 @@ namespace libcloudphxx / si::seconds * si::kilograms / si::cubic_metres; } - // Kessler/Beard terminal velocity + // Kessler/Beard rain terminal velocity // eq. 5d in Grabowski & Smolarkiewicz 1996 libcloudphxx_const(si::velocity, vterm_A, 36.34, si::metre_per_second) @@ -93,6 +93,18 @@ namespace libcloudphxx ); } + // slope of Marshall Palmer distribution for rain + // eq. A.1 in Grabowski 1999 + template + quantity::type, real_t> lambda_rain( + const quantity &rr, + const quantity &rhod_0 + ) { + using namespace common; + real_t N_0r = real_t(1e7); + return std::pow(pi() * moist_air::rho_w() * N_0r / (rhod_0 * rr), real_t(0.25)) / si::meters; + } + // mean mass of ice A particle // eq. A.15a in Grabowski 1999 template @@ -127,20 +139,19 @@ namespace libcloudphxx 3 * mu + 9 / 2 * std::pow(sigma, 2)) * si::cubic_meters; //mass of the average ice A particle: real_t delta = IWCS / (ria * rhod_0); - quantity m_a = delta * m_as + (real_t(1) - delta) * m_al; - return m_a; + return delta * m_as + (real_t(1) - delta) * m_al; } - // mean velocity of ice A particle + // mean terminal velocity of ice A particle // eq. A.15b in Grabowski 1999 template - quantity::type, real_t> velocity_a( + quantity velocity_iceA( const quantity &ria, const quantity &rhod_0 ) { using namespace common; //velocity of small ice A particle: - quantity::type, real_t> v_as = + quantity v_as = real_t(0.1) * si::meters / si::seconds; //velocity of large ice A particle: quantity IWCS = std::min( @@ -151,28 +162,18 @@ namespace libcloudphxx real_t(0.837)) * si::dimensionless()) * si::kilograms / si::cubic_meters; quantity IWCL = rhod_0 * ria - IWCS; - quantity::type, real_t> v_al = + quantity v_al = (real_t(0.9) + real_t(0.1) * std::log10(real_t(1e3) * IWCL / si::kilograms * si::cubic_meters)) * si::meters / si::seconds; //average velocity of ice A particle: real_t delta = IWCS / (ria * rhod_0); - quantity::type, real_t> v_a = - (delta * v_as + (real_t(1) - delta) * v_al) * std::pow( - real_t(0.3) / rhod_0 * si::kilograms / si::cubic_meters, real_t(0.5)); - return v_a; + return (delta * v_as + (real_t(1) - delta) * v_al) * std::pow( + real_t(0.3) / rhod_0 * si::kilograms / si::cubic_meters, real_t(0.5)); } - // slope of Marshall Palmer distribution for rain - // eq. A.1 in Grabowski 1999 - template - quantity::type, real_t> lambda_rain( - const quantity &rr, - const quantity &rhod_0 - ) { - using namespace common; - real_t N_0r = real_t(1e7); - return std::pow(pi() * moist_air::rho_w() * N_0r / (rhod_0 * rr), real_t(0.25)) / si::meters; - } + + // graupel density used for ice B (from Grabowski 1999) + libcloudphxx_const(si::mass_density, rho_ib, 400, si::kilograms / si::cubic_metres) // slope of Marshall Palmer distribution for ice B // eq. A.4 in Grabowski 1999 @@ -183,7 +184,7 @@ namespace libcloudphxx ) { using namespace common; real_t N_0b = real_t(4e6); - return std::pow(pi() * moist_air::rho_ib() * N_0b / (rhod_0 * rib), real_t(0.25)) / si::meters; + return std::pow(pi() * rho_ib() * N_0b / (rhod_0 * rib), real_t(0.25)) / si::meters; } @@ -195,8 +196,22 @@ namespace libcloudphxx const quantity &rhod_0 ) { using namespace common; - return pi() * moist_air::rho_ib() / (real_t(6) * std::pow( - lambda_ice_b(rib, rhod_0) * si::meters, real_t(3))) * si::cubic_meters; + return pi() * rho_ib() / (real_t(6) * std::pow(lambda_ice_b(rib, rhod_0) * si::meters, + real_t(3))) * si::cubic_meters; + } + + // mean terminal velocity of ice B particle + // eq. A.6 in Grabowski 1999 + template + quantity velocity_iceB( + const quantity &rib, + const quantity &rhod_0 + ) { + return real_t(31.2) * + std::pow(lambda_ice_b(rib, rhod_0) * si::meters, real_t(-0.37)) * std::pow( + rhod_0 / si::kilograms * si::cubic_meters, + real_t(-0.5)) * + si::meters / si::seconds; } //coefficients for ice growth by deposition/riming @@ -218,7 +233,7 @@ namespace libcloudphxx //temperature in Celsius: real_t Tc = T / si::kelvins - real_t(273.16); //limiting the temerature to (-31, 0): - real_t ttcoe = std::min(real_t(-0.), std::max(real_t(-31), Tc)); + real_t ttcoe = std::min(real_t(0.), std::max(real_t(-31), Tc)); int index1 = static_cast(std::trunc(-ttcoe)); //index between 0 and 31 int index2 = index1 + 1; // index between 1 and 32 real_t del = -ttcoe - real_t(index1); @@ -265,11 +280,9 @@ namespace libcloudphxx : real_t(0.1); quantity rv_adj = beta * rvs + (1 - beta) * rvsi; //homogeneous nucleation rate (only if T < -40 C) - quantity::type, real_t> hom1 = - (T < real_t(233.16) * si::kelvins) - ? (real_t(1) - std::exp(-dt / taunuc)) * std::max(real_t(0) * si::dimensionless(), rv - rv_adj) / dt - : real_t(0) / si::seconds; - return hom1; + return (T < real_t(233.16) * si::kelvins) + ? (real_t(1) - std::exp(-dt / taunuc)) * std::max(real_t(0) * si::dimensionless(), rv - rv_adj) / dt + : real_t(0) / si::seconds; } //ice A homogeneous nucleation 2 (rc -> ria) @@ -282,11 +295,9 @@ namespace libcloudphxx ) { const quantity taunuc = dt; // time scale for nucleation //homogeneous nucleation rate (only if T < -40 C) - quantity::type, real_t> hom2 = - (T < real_t(233.16) * si::kelvins) - ? (real_t(1) - std::exp(-dt / taunuc)) * rc / dt - : real_t(0) / si::seconds; - return hom2; + return (T < real_t(233.16) * si::kelvins) + ? (real_t(1) - std::exp(-dt / taunuc)) * rc / dt + : real_t(0) / si::seconds; } //ice A heterogeneous nucleation (rc-> ria) @@ -307,13 +318,11 @@ namespace libcloudphxx quantity::type, real_t> N_in = std::min(real_t(1e5), real_t(1e-2) * std::exp(real_t(0.6) * (real_t(273.16) - T / si::kelvins))) / si::cubic_meters; //nucleation rate: - quantity::type, real_t> het = - (real_t(1) - std::exp(-dt / taunuc)) * std::min(rc, std::max(real_t(0) * si::dimensionless(), - N_in * std::max( - real_t(1e-12) * si::kilograms, - m_a) / rhod_0 - - ria)) / dt; - return het; + return (real_t(1) - std::exp(-dt / taunuc)) * std::min(rc, std::max(real_t(0) * si::dimensionless(), + N_in * std::max( + real_t(1e-12) * si::kilograms, + m_a) / rhod_0 + - ria)) / dt; } // ice B heterogeneous nucleation 1 (rr -> rib) @@ -342,14 +351,12 @@ namespace libcloudphxx //mean ice A particle mass quantity m_a = mass_a(ria, T, rhod_0); //mean ice A velocity - real_t v_a = velocity_a(ria, rhod_0) * si::seconds / si::meters; + real_t v_a = velocity_iceA(ria, rhod_0) * si::seconds / si::meters; //rate of collisions between raindrops and ice A quantity::type, real_t> N_ra = real_t(1e7) / si::cubic_meters / si::meters / lambda_r * std::abs(v_r - v_a) * si::meters / si::seconds * pi() * R_r * R_r * ria / m_a; //nucleation rate: - quantity::type, real_t> het1 = - N_ra * m_r; - return het1; + return N_ra * m_r; } //ice B heterogeneous nucleation 2 (ria -> rib) @@ -375,14 +382,12 @@ namespace libcloudphxx //mean ice A particle mass quantity m_a = mass_a(ria, T, rhod_0); //mean ice A velocity - real_t v_a = velocity_a(ria, rhod_0) * si::seconds / si::meters; + real_t v_a = velocity_iceA(ria, rhod_0) * si::seconds / si::meters; //rate of collisions between raindrops and ice A quantity::type, real_t> N_ra = real_t(1e7) / si::cubic_meters / si::meters / lambda_r * abs(v_r - v_a) * si::meters / si::seconds * pi() * R_r * R_r * ria / m_a; //nucleation rate: - quantity::type, real_t> het2 = - N_ra * m_a; - return het2; + return N_ra * m_a; } //melting of ice A (ria -> rr) @@ -403,15 +408,13 @@ namespace libcloudphxx //diameter of average ice A particle: quantity D_a = std::pow(m_a / real_t(0.025) / si::kilograms, real_t(0.5)) * si::meters; //fall velocity of avg. ice A particle: - quantity::type, real_t> v_a = real_t(4) * std::pow( - D_a / si::meters, real_t(0.25)) * si::meters / si::seconds; + quantity v_a = velocity_iceA(ria, rhod_0); //Reynolds number: real_t Re = D_a * v_a * rhod_0 / vterm::visc(T); //ventilation factor: real_t F_a = real_t(0.78) + real_t(0.27) * std::pow(Re, real_t(0.5)); - quantity::type, real_t> mela = real_t(4.5e-7) * ria / + return real_t(4.5e-7) * ria / m_a * D_a * (273.16 - T / si::kelvins) * F_a * si::kilograms / si::meters / si::seconds; - return mela; } //melting of ice B (rib -> rr) @@ -434,17 +437,13 @@ namespace libcloudphxx //diameter of average ice B particle: quantity D_b = real_t(1) / lambda_b; //fall velocity of avg. ice B particle: - quantity::type, real_t> v_b = real_t(31.2) * - std::pow(lambda_b * si::meters, real_t(-0.37)) * std::pow(rhod_0 / si::kilograms * si::cubic_meters, - real_t(-0.5)) * - si::meters / si::seconds; + quantity v_b = velocity_iceB(rib, rhod_0); //Reynolds number: real_t Re = D_b * v_b * rhod_0 / vterm::visc(T); //ventilation factor: real_t F_b = real_t(0.78) + real_t(0.27) * std::pow(Re, real_t(0.5)); - quantity::type, real_t> melb = real_t(4.5e-7) * rib / + return real_t(4.5e-7) * rib / m_b * D_b * (273.16 - T / si::kelvins) * F_b * si::kilograms / si::meters / si::seconds; - return melb; } @@ -470,8 +469,7 @@ namespace libcloudphxx // growth rate for a single particle: quantity::type, real_t> dm_dt_AE = (rv - rvsi) / (rvs - rvsi) * alpha * std::pow(m_a / si::kilograms, beta) * si::kilograms / si::seconds; - quantity::type, real_t> depa = ria / m_a * dm_dt_AE; - return depa; + return ria / m_a * dm_dt_AE; } @@ -545,8 +543,7 @@ namespace libcloudphxx // growth rate for a single particle: quantity::type, real_t> dm_dt_AE = (rv - rvsi) / (rvs - rvsi) * alpha * std::pow(m_b / si::kilograms, beta) * si::kilograms / si::seconds; - quantity::type, real_t> depb = rib / m_b * dm_dt_AE; - return depb; + return rib / m_b * dm_dt_AE; } //growth of ice B by riming (rc, rr ->rib) diff --git a/include/libcloudph++/blk_1m/rhs_columnwise.hpp b/include/libcloudph++/blk_1m/rhs_columnwise.hpp index 8e410e14..bb0fd956 100644 --- a/include/libcloudph++/blk_1m/rhs_columnwise.hpp +++ b/include/libcloudph++/blk_1m/rhs_columnwise.hpp @@ -20,69 +20,120 @@ namespace libcloudphxx template real_t rhs_columnwise( const opts_t &opts, - cont_t &dot_rr_cont, + cont_t &dot_r_cont, const cont_t &rhod_cont, - const cont_t &rr_cont, - const real_t &dz + const cont_t &r_cont, + const real_t &dz, + const std::string& precip_type = "rain" ) // { using flux_t = quantity::type, real_t>; - auto dot_rr_unit = si::hertz; + auto dot_r_unit = si::hertz; if (!opts.sedi) return 0; // flux_t flux_in = 0 * si::kilograms / si::cubic_metres / si::seconds; - real_t *dot_rr = NULL; + real_t *dot_r = NULL; const real_t zero = 0; // this should give zero flux from above the domain top - const real_t *rhod = &*(--(rhod_cont.end())), *rr = &zero; + const real_t *rhod = &*(--(rhod_cont.end())), *r = &zero; - auto iter = zip(dot_rr_cont, rhod_cont, rr_cont); + auto iter = zip(dot_r_cont, rhod_cont, r_cont); for (auto tup_ptr = iter.end(); tup_ptr != iter.begin();) { --tup_ptr; const real_t *rhod_below = &std::get<1>(*tup_ptr), - *rr_below = &std::get<2>(*tup_ptr); + *r_below = &std::get<2>(*tup_ptr); - if (dot_rr != NULL) // i.e. all but first (top) grid cell + if (dot_r != NULL) // i.e. all but first (top) grid cell { - // terminal momenta at grid-cell edge (to assure precip mass conservation) - flux_t flux_out = -real_t(.5) * ( // averaging + axis orientation - (*rhod_below * si::kilograms / si::cubic_metres) * formulae::v_term( - *rr_below * si::kilograms / si::kilograms, - *rhod_below * si::kilograms / si::cubic_metres, - *rhod_cont.begin() * si::kilograms / si::cubic_metres - ) + - (*rhod * si::kilograms / si::cubic_metres) * formulae::v_term( - *rr * si::kilograms / si::kilograms, - *rhod * si::kilograms / si::cubic_metres, - *rhod_cont.begin() * si::kilograms / si::cubic_metres - ) - ) * (*rr * si::kilograms / si::kilograms) / (dz * si::metres); + flux_t flux_out = real_t(0) * si::kilograms / si::cubic_metres / si::seconds; + if (precip_type == "rain") + { + // terminal momenta at grid-cell edge (to assure precip mass conservation) + flux_out += -real_t(.5) * ( // averaging + axis orientation + (*rhod_below * si::kilograms / si::cubic_metres) * formulae::v_term( + *r_below * si::kilograms / si::kilograms, + *rhod_below * si::kilograms / si::cubic_metres, + *rhod_cont.begin() * si::kilograms / si::cubic_metres + ) + + (*rhod * si::kilograms / si::cubic_metres) * formulae::v_term( + *r * si::kilograms / si::kilograms, + *rhod * si::kilograms / si::cubic_metres, + *rhod_cont.begin() * si::kilograms / si::cubic_metres + ) + ) * (*r * si::kilograms / si::kilograms) / (dz * si::metres); + } + else if (precip_type == "iceA") + { + // terminal momenta at grid-cell edge (to assure precip mass conservation) + flux_out += -real_t(.5) * ( // averaging + axis orientation + (*rhod_below * si::kilograms / si::cubic_metres) * formulae::velocity_iceA( + *r_below * si::kilograms / si::kilograms, + *rhod_below * si::kilograms / si::cubic_metres + ) + + (*rhod * si::kilograms / si::cubic_metres) * formulae::velocity_iceA( + *r * si::kilograms / si::kilograms, + *rhod * si::kilograms / si::cubic_metres + ) + ) * (*r * si::kilograms / si::kilograms) / (dz * si::metres); + } + else if (precip_type == "iceB") + { + // terminal momenta at grid-cell edge (to assure precip mass conservation) + flux_out += -real_t(.5) * ( // averaging + axis orientation + (*rhod_below * si::kilograms / si::cubic_metres) * formulae::velocity_iceB( + *r_below * si::kilograms / si::kilograms, + *rhod_below * si::kilograms / si::cubic_metres + ) + + (*rhod * si::kilograms / si::cubic_metres) * formulae::velocity_iceB( + *r * si::kilograms / si::kilograms, + *rhod * si::kilograms / si::cubic_metres + ) + ) * (*r * si::kilograms / si::kilograms) / (dz * si::metres); + } - *dot_rr -= (flux_in - flux_out) / (*rhod * si::kilograms / si::cubic_metres) / dot_rr_unit; + *dot_r -= (flux_in - flux_out) / (*rhod * si::kilograms / si::cubic_metres) / dot_r_unit; flux_in = flux_out; // inflow = outflow from above } - dot_rr = &std::get<0>(*tup_ptr); + dot_r = &std::get<0>(*tup_ptr); rhod = rhod_below; - rr = rr_below; + r = r_below; } // the bottom grid cell (with mid-cell vterm approximation) - flux_t flux_out = - (*rhod * si::kilograms / si::cubic_metres) * formulae::v_term( - *rr * si::kilograms / si::kilograms, - *rhod * si::kilograms / si::cubic_metres, - *rhod_cont.begin() * si::kilograms / si::cubic_metres - ) * (*rr * si::kilograms / si::kilograms) / (dz * si::metres); - *dot_rr -= (flux_in - flux_out) / (*rhod * si::kilograms / si::cubic_metres) / dot_rr_unit; + flux_t flux_out = real_t(0) * si::kilograms / si::cubic_metres / si::seconds; + if (precip_type == "rain") + { + flux_out += - (*rhod * si::kilograms / si::cubic_metres) * formulae::v_term( + *r * si::kilograms / si::kilograms, + *rhod * si::kilograms / si::cubic_metres, + *rhod_cont.begin() * si::kilograms / si::cubic_metres + ) * (*r * si::kilograms / si::kilograms) / (dz * si::metres); + } + else if (precip_type == "iceA") + { + flux_out += - (*rhod * si::kilograms / si::cubic_metres) * formulae::velocity_iceA( + *r * si::kilograms / si::kilograms, + *rhod * si::kilograms / si::cubic_metres + ) * (*r * si::kilograms / si::kilograms) / (dz * si::metres); + } + else if (precip_type == "iceB") + { + flux_out += - (*rhod * si::kilograms / si::cubic_metres) * formulae::velocity_iceB( + *r * si::kilograms / si::kilograms, + *rhod * si::kilograms / si::cubic_metres + ) * (*r * si::kilograms / si::kilograms) / (dz * si::metres); + } + *dot_r -= (flux_in - flux_out) / (*rhod * si::kilograms / si::cubic_metres) / dot_r_unit; // outflow from the domain return real_t(flux_out / (si::kilograms / si::cubic_metres / si::seconds)); } diff --git a/include/libcloudph++/common/moist_air.hpp b/include/libcloudph++/common/moist_air.hpp index 74c5a318..159b04cf 100644 --- a/include/libcloudph++/common/moist_air.hpp +++ b/include/libcloudph++/common/moist_air.hpp @@ -50,8 +50,6 @@ namespace libcloudphxx libcloudphxx_const(si::mass_density, rho_w, 1e3, si::kilograms / si::cubic_metres) // ice density libcloudphxx_const(si::mass_density, rho_i, 910, si::kilograms / si::cubic_metres) - // graupel density used for ice B (from Grabowski 1999) - libcloudphxx_const(si::mass_density, rho_ib, 400, si::kilograms / si::cubic_metres) // mixing rule for extensive quantitites (i.e. using mass mixing ratio) template diff --git a/tests/python/unit/api_blk_1m.py b/tests/python/unit/api_blk_1m.py index 7c6af568..d36a8a49 100644 --- a/tests/python/unit/api_blk_1m.py +++ b/tests/python/unit/api_blk_1m.py @@ -87,7 +87,7 @@ rr = arr_t([0.]) dot_rr_old = dot_rr.copy() -flux = blk_1m.rhs_columnwise(opts, dot_rr, rhod, rr, dz) +flux = blk_1m.rhs_columnwise(opts, dot_rr, rhod, rr, dz, "rain") assert flux == 0 assert dot_rr == dot_rr_old # no rain water -> no precip @@ -101,4 +101,10 @@ dot_rib = arr_t([0.]) blk_1m.rhs_cellwise_nwtrph_ice(opts, dot_th, dot_rv, dot_rc, dot_rr, dot_ria, dot_rib, rhod, p, th, rv, rc, rr, ria, rib, dt) assert dot_ria != 0 -assert dot_rib != 0 \ No newline at end of file +assert dot_rib != 0 + +#testing sedimentation of ice +flux_iceA = blk_1m.rhs_columnwise(opts, dot_ria, rhod, ria, dz, "iceA") +flux_iceB = blk_1m.rhs_columnwise(opts, dot_rib, rhod, rib, dz, "iceB") +assert flux_iceA != 0 +assert flux_iceB != 0 \ No newline at end of file From c225df4faabb900ef72f5495d35fe0227608e0e6 Mon Sep 17 00:00:00 2001 From: Agnieszka Makulska Date: Thu, 14 Nov 2024 11:21:53 +0100 Subject: [PATCH 13/17] fixing ice sedimentation --- bindings/python/blk_1m.hpp | 35 +++- bindings/python/lib.cpp | 1 + include/libcloudph++/blk_1m/formulae.hpp | 62 +++---- .../libcloudph++/blk_1m/rhs_columnwise.hpp | 152 ++++++++++++------ tests/python/unit/api_blk_1m.py | 6 +- 5 files changed, 163 insertions(+), 93 deletions(-) diff --git a/bindings/python/blk_1m.hpp b/bindings/python/blk_1m.hpp index ba912128..3f6c2dab 100644 --- a/bindings/python/blk_1m.hpp +++ b/bindings/python/blk_1m.hpp @@ -201,21 +201,40 @@ namespace libcloudphxx { template typename arr_t::T_numtype rhs_columnwise( const b1m::opts_t &opts, - bp_array &dot_r, + bp_array &dot_rr, const bp_array &rhod, - const bp_array &r, - const typename arr_t::T_numtype &dz, - const std::string& precip_type = "rain" + const bp_array &rr, + const typename arr_t::T_numtype &dz ) { arr_t - np2bz_dot_r(np2bz(dot_r)); + np2bz_dot_rr(np2bz(dot_rr)); return b1m::rhs_columnwise( opts, - np2bz_dot_r, + np2bz_dot_rr, + np2bz(rhod), + np2bz(rr), + dz + ); + } + + template + typename arr_t::T_numtype rhs_columnwise_ice( + const b1m::opts_t &opts, + bp_array &dot_ri, + const bp_array &rhod, + const bp_array &ri, + const typename arr_t::T_numtype &dz, + const std::string ice_type + ) { + arr_t + np2bz_dot_ri(np2bz(dot_ri)); + return b1m::rhs_columnwise_ice( + opts, + np2bz_dot_ri, np2bz(rhod), - np2bz(r), + np2bz(ri), dz, - precip_type + ice_type ); } }; diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index 6f3253f7..4bc0e9c0 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -176,6 +176,7 @@ BOOST_PYTHON_MODULE(libcloudphxx) bp::def("rhs_cellwise_nwtrph", blk_1m::rhs_cellwise_nwtrph); bp::def("rhs_cellwise_nwtrph_ice", blk_1m::rhs_cellwise_nwtrph_ice); bp::def("rhs_columnwise", blk_1m::rhs_columnwise); // TODO: handle the returned flux + bp::def("rhs_columnwise_ice", blk_1m::rhs_columnwise_ice); } // blk_2m stuff diff --git a/include/libcloudph++/blk_1m/formulae.hpp b/include/libcloudph++/blk_1m/formulae.hpp index cd07c282..ce12917e 100644 --- a/include/libcloudph++/blk_1m/formulae.hpp +++ b/include/libcloudph++/blk_1m/formulae.hpp @@ -106,7 +106,7 @@ namespace libcloudphxx } // mean mass of ice A particle - // eq. A.15a in Grabowski 1999 + // eq. A.7 - A.15a in Grabowski 1999 template quantity mass_a( const quantity &ria, @@ -257,7 +257,7 @@ namespace libcloudphxx //temperature in Celsius: real_t Tc = T / si::kelvins - real_t(273.16); //limiting the temerature to (-31, 0): - real_t ttcoe = std::min(real_t(-0.), std::max(real_t(-31), Tc)); + real_t ttcoe = std::min(real_t(0.), std::max(real_t(-31), Tc)); int index1 = static_cast(std::trunc(-ttcoe)); //index between 0 and 31 int index2 = index1 + 1; // index between 1 and 32 real_t del = -ttcoe - real_t(index1); @@ -335,25 +335,26 @@ namespace libcloudphxx const quantity &rhod_0 ) { using namespace common; - if (ria == 0) // ice B is formed only if ice A is present + if (ria == 0 || rr == 0) // ice B is formed only if ice A and rain are present { return real_t(0) / si::seconds; } quantity::type, real_t> lambda_r = lambda_rain(rr, rhod_0); - //mean raindrop mass + //mean raindrop mass (eq. A2 in Grabowski 1999) quantity m_r = pi() * moist_air::rho_w() / (real_t(6) * std::pow( lambda_r * si::meters, 3)) * si::cubic_meters; - //mean raindrop velocity + //mean raindrop velocity (eq. A3 in Grabowski 1999) real_t v_r = real_t(251) * std::pow(lambda_r * si::meters * rhod_0 / si::kilograms * si::cubic_meters, real_t(-0.5)); - //mean raindrop radius + //mean raindrop radius (eq. A2 in Grabowski 1999) quantity R_r = real_t(0.5) / lambda_r; //mean ice A particle mass quantity m_a = mass_a(ria, T, rhod_0); //mean ice A velocity real_t v_a = velocity_iceA(ria, rhod_0) * si::seconds / si::meters; //rate of collisions between raindrops and ice A - quantity::type, real_t> N_ra = real_t(1e7) / si::cubic_meters / + real_t N_0r = real_t(1e7); + quantity::type, real_t> N_ra = N_0r / si::cubic_meters / si::meters / lambda_r * std::abs(v_r - v_a) * si::meters / si::seconds * pi() * R_r * R_r * ria / m_a; //nucleation rate: return N_ra * m_r; @@ -369,22 +370,23 @@ namespace libcloudphxx const quantity &rhod_0 ) { using namespace common; - if (ria == 0) // ice B is formed only if ice A is present + if (ria == 0 || rr == 0) // ice B is formed only if ice A and rain are present { return real_t(0) / si::seconds; } quantity::type, real_t> lambda_r = lambda_rain(rr, rhod_0); - //mean raindrop velocity + //mean raindrop velocity (eq. A3 in Grabowski 1999) real_t v_r = real_t(251) * std::pow(lambda_r * si::meters * rhod_0 / si::kilograms * si::cubic_meters, real_t(-0.5)); - //mean raindrop radius + //mean raindrop radius (eq. A2 in Grabowski 1999) quantity R_r = real_t(0.5) / lambda_r; //mean ice A particle mass quantity m_a = mass_a(ria, T, rhod_0); //mean ice A velocity real_t v_a = velocity_iceA(ria, rhod_0) * si::seconds / si::meters; //rate of collisions between raindrops and ice A - quantity::type, real_t> N_ra = real_t(1e7) / si::cubic_meters / + real_t N_0r = real_t(1e7); + quantity::type, real_t> N_ra = N_0r / si::cubic_meters / si::meters / lambda_r * abs(v_r - v_a) * si::meters / si::seconds * pi() * R_r * R_r * ria / m_a; //nucleation rate: return N_ra * m_a; @@ -399,7 +401,7 @@ namespace libcloudphxx const quantity &rhod_0 ) { using namespace common; - if (ria == 0) // calculating melting rate only if ice A is present + if (ria == 0 || T < 273.16 * si::kelvin) // calculating melting rate only if ice A is present and T > 0 C { return real_t(0) / si::seconds; } @@ -426,7 +428,7 @@ namespace libcloudphxx const quantity &rhod_0 ) { using namespace common; - if (rib == 0) // calculating melting rate only if ice B is present + if (rib == 0 || T < 273.16 * si::kelvin) // calculating melting rate only if ice B is present and T > 0 C { return real_t(0) / si::seconds; } @@ -434,7 +436,7 @@ namespace libcloudphxx lambda_ice_b(rib, rhod_0); //mass of average ice B particle: quantity m_b = mass_b(rib, rhod_0); - //diameter of average ice B particle: + //diameter of average ice B particle (eq. A5 in Grabowski 1999): quantity D_b = real_t(1) / lambda_b; //fall velocity of avg. ice B particle: quantity v_b = velocity_iceB(rib, rhod_0); @@ -467,8 +469,8 @@ namespace libcloudphxx real_t alpha = coeff_alpha(T); real_t beta = coeff_beta(T); // growth rate for a single particle: - quantity::type, real_t> dm_dt_AE = (rv - rvsi) / (rvs - rvsi) * alpha * - std::pow(m_a / si::kilograms, beta) * si::kilograms / si::seconds; + quantity::type, real_t> dm_dt_AE = real_t(1e-3) * (rv - rvsi) / (rvs - rvsi) * alpha * + std::pow(m_a * real_t(1e3) / si::kilograms, beta) * si::kilograms / si::seconds; //1e-3 comes from conversion g/sec into kg/sec return ria / m_a * dm_dt_AE; } @@ -494,27 +496,27 @@ namespace libcloudphxx real_t alpha = coeff_alpha(T); real_t beta = coeff_beta(T); // regime AE - quantity::type, real_t> dm_dt_AE = (rv - rvsi) / (rvs - rvsi) * alpha * - std::pow(m_a / si::kilograms, beta) * si::kilograms / si::seconds; + quantity::type, real_t> dm_dt_AE = real_t(1e-3) * (rv - rvsi) / (rvs - rvsi) * alpha * + std::pow(m_a * real_t(1e3) / si::kilograms, beta) * si::kilograms / si::seconds; //1e-3 comes from conversion g/sec into kg/sec // regime BC real_t tan_theta = real_t(1.) + real_t(0.1) * std::log( rhod_0 * rc * real_t(1e3) / si::kilograms * si::cubic_meters); real_t gamma = alpha * std::pow(real_t(5e-8), beta); - quantity::type, real_t> dm_dt_BC = gamma * std::pow( + quantity::type, real_t> dm_dt_BC = real_t(1e-3) * gamma * std::pow( m_a / real_t(5e-11) / si::kilograms, tan_theta) * si::kilograms / si::seconds; //regime CD real_t dzeta = gamma * std::pow(real_t(2e3), tan_theta); real_t xi = std::log(rc * rhod_0 / si::kilograms * si::cubic_meters * real_t(1e9) / dzeta) / std::log( real_t(1e4)); - quantity::type, real_t> dm_dt_CD = dzeta * std::pow( + quantity::type, real_t> dm_dt_CD = real_t(1e-3) * dzeta * std::pow( m_a * real_t(1e7) / si::kilograms, xi) * si::kilograms / si::seconds; //total growth quantity::type, real_t> rima = real_t(0) / si::seconds; - if (m_a > real_t(5e-11) * si::kilograms && m_a < real_t(1e-7) * si::kilograms) + if (m_a > real_t(5e-11) * si::kilograms && m_a <= real_t(1e-7) * si::kilograms) { rima += std::max(real_t(0) * si::kilograms / si::seconds, dm_dt_BC - dm_dt_AE) * ria / m_a; } - if (m_a >= real_t(1e-7) * si::kilograms) + if (m_a > real_t(1e-7) * si::kilograms) { rima += std::max(real_t(0) * si::kilograms / si::seconds, dm_dt_CD - dm_dt_AE) * ria / m_a; } @@ -541,8 +543,8 @@ namespace libcloudphxx real_t alpha = coeff_alpha(T); real_t beta = coeff_beta(T); // growth rate for a single particle: - quantity::type, real_t> dm_dt_AE = (rv - rvsi) / (rvs - rvsi) * alpha * - std::pow(m_b / si::kilograms, beta) * si::kilograms / si::seconds; + quantity::type, real_t> dm_dt_AE = real_t(1e-3) * (rv - rvsi) / (rvs - rvsi) * alpha * + std::pow(m_b * real_t(1e3) / si::kilograms, beta) * si::kilograms / si::seconds; //1e-3 comes from conversion g/sec into kg/sec return rib / m_b * dm_dt_AE; } @@ -567,27 +569,27 @@ namespace libcloudphxx real_t alpha = coeff_alpha(T); real_t beta = coeff_beta(T); // regime AE - quantity::type, real_t> dm_dt_AE = (rv - rvsi) / (rvs - rvsi) * alpha * - std::pow(m_b / si::kilograms, beta) * si::kilograms / si::seconds; + quantity::type, real_t> dm_dt_AE = real_t(1e-3) * (rv - rvsi) / (rvs - rvsi) * alpha * + std::pow(m_b * real_t(1e3) / si::kilograms, beta) * si::kilograms / si::seconds; //1e-3 comes from conversion g/sec into kg/sec // regime BC real_t tan_theta = real_t(1.) + real_t(0.1) * std::log( rhod_0 * rc * real_t(1e3) / si::kilograms * si::cubic_meters); real_t gamma = alpha * std::pow(real_t(5e-8), beta); - quantity::type, real_t> dm_dt_BC = gamma * std::pow( + quantity::type, real_t> dm_dt_BC = real_t(1e-3) * gamma * std::pow( m_b / real_t(5e-11) / si::kilograms, tan_theta) * si::kilograms / si::seconds; //regime CD real_t dzeta = gamma * std::pow(real_t(2e3), tan_theta); real_t xi = std::log(rc * rhod_0 / si::kilograms * si::cubic_meters * real_t(1e9) / dzeta) / std::log( real_t(1e4)); - quantity::type, real_t> dm_dt_CD = dzeta * std::pow( + quantity::type, real_t> dm_dt_CD = real_t(1e-3) * dzeta * std::pow( m_b * real_t(1e7) / si::kilograms, xi) * si::kilograms / si::seconds; //total growth quantity::type, real_t> rimb = real_t(0) / si::seconds; - if (m_b > real_t(5e-11) * si::kilograms && m_b < real_t(1e-7) * si::kilograms) + if (m_b > real_t(5e-11) * si::kilograms && m_b <= real_t(1e-7) * si::kilograms) { rimb += std::max(real_t(0) * si::kilograms / si::seconds, dm_dt_BC - dm_dt_AE) * rib / m_b; } - if (m_b >= real_t(1e-7) * si::kilograms) + if (m_b > real_t(1e-7) * si::kilograms) { rimb += std::max(real_t(0) * si::kilograms / si::seconds, dm_dt_CD - dm_dt_AE) * rib / m_b; } diff --git a/include/libcloudph++/blk_1m/rhs_columnwise.hpp b/include/libcloudph++/blk_1m/rhs_columnwise.hpp index bb0fd956..6ab16ab0 100644 --- a/include/libcloudph++/blk_1m/rhs_columnwise.hpp +++ b/include/libcloudph++/blk_1m/rhs_columnwise.hpp @@ -1,6 +1,6 @@ /** @file * @copyright University of Warsaw - * @brief Rain sedimentation representation for single-moment bulk microphysics + * @brief Rain and ice sedimentation representation for single-moment bulk microphysics * using forcing terms based on the upstrem advection scheme * @section LICENSE * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/) @@ -20,120 +20,168 @@ namespace libcloudphxx template real_t rhs_columnwise( const opts_t &opts, - cont_t &dot_r_cont, + cont_t &dot_rr_cont, const cont_t &rhod_cont, - const cont_t &r_cont, + const cont_t &rr_cont, + const real_t &dz + ) +// + { + using flux_t = quantity::type, real_t>; + + auto dot_rr_unit = si::hertz; + + if (!opts.sedi) return 0; + + // + flux_t flux_in = 0 * si::kilograms / si::cubic_metres / si::seconds; + real_t *dot_rr = NULL; + const real_t zero = 0; + + // this should give zero flux from above the domain top + const real_t *rhod = &*(--(rhod_cont.end())), *rr = &zero; + + auto iter = zip(dot_rr_cont, rhod_cont, rr_cont); + for (auto tup_ptr = iter.end(); tup_ptr != iter.begin();) + { + --tup_ptr; + + const real_t + *rhod_below = &std::get<1>(*tup_ptr), + *rr_below = &std::get<2>(*tup_ptr); + + if (dot_rr != NULL) // i.e. all but first (top) grid cell + { + // terminal momenta at grid-cell edge (to assure precip mass conservation) + flux_t flux_out = -real_t(.5) * ( // averaging + axis orientation + (*rhod_below * si::kilograms / si::cubic_metres) * formulae::v_term( + *rr_below * si::kilograms / si::kilograms, + *rhod_below * si::kilograms / si::cubic_metres, + *rhod_cont.begin() * si::kilograms / si::cubic_metres + ) + + (*rhod * si::kilograms / si::cubic_metres) * formulae::v_term( + *rr * si::kilograms / si::kilograms, + *rhod * si::kilograms / si::cubic_metres, + *rhod_cont.begin() * si::kilograms / si::cubic_metres + ) + ) * (*rr * si::kilograms / si::kilograms) / (dz * si::metres); + + *dot_rr -= (flux_in - flux_out) / (*rhod * si::kilograms / si::cubic_metres) / dot_rr_unit; + flux_in = flux_out; // inflow = outflow from above + } + + dot_rr = &std::get<0>(*tup_ptr); + rhod = rhod_below; + rr = rr_below; + } + + // the bottom grid cell (with mid-cell vterm approximation) + flux_t flux_out = - (*rhod * si::kilograms / si::cubic_metres) * formulae::v_term( + *rr * si::kilograms / si::kilograms, + *rhod * si::kilograms / si::cubic_metres, + *rhod_cont.begin() * si::kilograms / si::cubic_metres + ) * (*rr * si::kilograms / si::kilograms) / (dz * si::metres); + *dot_rr -= (flux_in - flux_out) / (*rhod * si::kilograms / si::cubic_metres) / dot_rr_unit; + + // outflow from the domain + return real_t(flux_out / (si::kilograms / si::cubic_metres / si::seconds)); + } + + +//+ template + real_t rhs_columnwise_ice( + const opts_t &opts, + cont_t &dot_ri_cont, + const cont_t &rhod_cont, + const cont_t &ri_cont, const real_t &dz, - const std::string& precip_type = "rain" + const std::string& ice_type ) // { using flux_t = quantity::type, real_t>; - auto dot_r_unit = si::hertz; + auto dot_ri_unit = si::hertz; if (!opts.sedi) return 0; // flux_t flux_in = 0 * si::kilograms / si::cubic_metres / si::seconds; - real_t *dot_r = NULL; + real_t *dot_ri = NULL; const real_t zero = 0; // this should give zero flux from above the domain top - const real_t *rhod = &*(--(rhod_cont.end())), *r = &zero; + const real_t *rhod = &*(--(rhod_cont.end())), *ri = &zero; - auto iter = zip(dot_r_cont, rhod_cont, r_cont); + auto iter = zip(dot_ri_cont, rhod_cont, ri_cont); for (auto tup_ptr = iter.end(); tup_ptr != iter.begin();) { --tup_ptr; const real_t *rhod_below = &std::get<1>(*tup_ptr), - *r_below = &std::get<2>(*tup_ptr); + *ri_below = &std::get<2>(*tup_ptr); - if (dot_r != NULL) // i.e. all but first (top) grid cell + if (dot_ri != NULL) // i.e. all but first (top) grid cell { flux_t flux_out = real_t(0) * si::kilograms / si::cubic_metres / si::seconds; - if (precip_type == "rain") - { - // terminal momenta at grid-cell edge (to assure precip mass conservation) - flux_out += -real_t(.5) * ( // averaging + axis orientation - (*rhod_below * si::kilograms / si::cubic_metres) * formulae::v_term( - *r_below * si::kilograms / si::kilograms, - *rhod_below * si::kilograms / si::cubic_metres, - *rhod_cont.begin() * si::kilograms / si::cubic_metres - ) + - (*rhod * si::kilograms / si::cubic_metres) * formulae::v_term( - *r * si::kilograms / si::kilograms, - *rhod * si::kilograms / si::cubic_metres, - *rhod_cont.begin() * si::kilograms / si::cubic_metres - ) - ) * (*r * si::kilograms / si::kilograms) / (dz * si::metres); - } - else if (precip_type == "iceA") + if (ice_type == "iceA") { // terminal momenta at grid-cell edge (to assure precip mass conservation) flux_out += -real_t(.5) * ( // averaging + axis orientation (*rhod_below * si::kilograms / si::cubic_metres) * formulae::velocity_iceA( - *r_below * si::kilograms / si::kilograms, + *ri_below * si::kilograms / si::kilograms, *rhod_below * si::kilograms / si::cubic_metres ) + (*rhod * si::kilograms / si::cubic_metres) * formulae::velocity_iceA( - *r * si::kilograms / si::kilograms, + *ri * si::kilograms / si::kilograms, *rhod * si::kilograms / si::cubic_metres ) - ) * (*r * si::kilograms / si::kilograms) / (dz * si::metres); + ) * (*ri * si::kilograms / si::kilograms) / (dz * si::metres); } - else if (precip_type == "iceB") + else if (ice_type == "iceB") { // terminal momenta at grid-cell edge (to assure precip mass conservation) flux_out += -real_t(.5) * ( // averaging + axis orientation (*rhod_below * si::kilograms / si::cubic_metres) * formulae::velocity_iceB( - *r_below * si::kilograms / si::kilograms, + *ri_below * si::kilograms / si::kilograms, *rhod_below * si::kilograms / si::cubic_metres ) + (*rhod * si::kilograms / si::cubic_metres) * formulae::velocity_iceB( - *r * si::kilograms / si::kilograms, + *ri * si::kilograms / si::kilograms, *rhod * si::kilograms / si::cubic_metres ) - ) * (*r * si::kilograms / si::kilograms) / (dz * si::metres); + ) * (*ri * si::kilograms / si::kilograms) / (dz * si::metres); } - *dot_r -= (flux_in - flux_out) / (*rhod * si::kilograms / si::cubic_metres) / dot_r_unit; + *dot_ri -= (flux_in - flux_out) / (*rhod * si::kilograms / si::cubic_metres) / dot_ri_unit; flux_in = flux_out; // inflow = outflow from above } - dot_r = &std::get<0>(*tup_ptr); + dot_ri = &std::get<0>(*tup_ptr); rhod = rhod_below; - r = r_below; + ri = ri_below; } // the bottom grid cell (with mid-cell vterm approximation) flux_t flux_out = real_t(0) * si::kilograms / si::cubic_metres / si::seconds; - if (precip_type == "rain") - { - flux_out += - (*rhod * si::kilograms / si::cubic_metres) * formulae::v_term( - *r * si::kilograms / si::kilograms, - *rhod * si::kilograms / si::cubic_metres, - *rhod_cont.begin() * si::kilograms / si::cubic_metres - ) * (*r * si::kilograms / si::kilograms) / (dz * si::metres); - } - else if (precip_type == "iceA") + if (ice_type == "iceA") { flux_out += - (*rhod * si::kilograms / si::cubic_metres) * formulae::velocity_iceA( - *r * si::kilograms / si::kilograms, + *ri * si::kilograms / si::kilograms, *rhod * si::kilograms / si::cubic_metres - ) * (*r * si::kilograms / si::kilograms) / (dz * si::metres); + ) * (*ri * si::kilograms / si::kilograms) / (dz * si::metres); } - else if (precip_type == "iceB") + else if (ice_type == "iceB") { flux_out += - (*rhod * si::kilograms / si::cubic_metres) * formulae::velocity_iceB( - *r * si::kilograms / si::kilograms, + *ri * si::kilograms / si::kilograms, *rhod * si::kilograms / si::cubic_metres - ) * (*r * si::kilograms / si::kilograms) / (dz * si::metres); + ) * (*ri * si::kilograms / si::kilograms) / (dz * si::metres); } - *dot_r -= (flux_in - flux_out) / (*rhod * si::kilograms / si::cubic_metres) / dot_r_unit; + *dot_ri -= (flux_in - flux_out) / (*rhod * si::kilograms / si::cubic_metres) / dot_ri_unit; // outflow from the domain return real_t(flux_out / (si::kilograms / si::cubic_metres / si::seconds)); } diff --git a/tests/python/unit/api_blk_1m.py b/tests/python/unit/api_blk_1m.py index d36a8a49..4476e062 100644 --- a/tests/python/unit/api_blk_1m.py +++ b/tests/python/unit/api_blk_1m.py @@ -87,7 +87,7 @@ rr = arr_t([0.]) dot_rr_old = dot_rr.copy() -flux = blk_1m.rhs_columnwise(opts, dot_rr, rhod, rr, dz, "rain") +flux = blk_1m.rhs_columnwise(opts, dot_rr, rhod, rr, dz) assert flux == 0 assert dot_rr == dot_rr_old # no rain water -> no precip @@ -104,7 +104,7 @@ assert dot_rib != 0 #testing sedimentation of ice -flux_iceA = blk_1m.rhs_columnwise(opts, dot_ria, rhod, ria, dz, "iceA") -flux_iceB = blk_1m.rhs_columnwise(opts, dot_rib, rhod, rib, dz, "iceB") +flux_iceA = blk_1m.rhs_columnwise_ice(opts, dot_ria, rhod, ria, dz, "iceA") +flux_iceB = blk_1m.rhs_columnwise_ice(opts, dot_rib, rhod, rib, dz, "iceB") assert flux_iceA != 0 assert flux_iceB != 0 \ No newline at end of file From a3f3a6b3d69ddc20c5890088317155fa0831438f Mon Sep 17 00:00:00 2001 From: Agnieszka Makulska Date: Fri, 15 Nov 2024 15:17:08 +0100 Subject: [PATCH 14/17] fixing the limiting of variable changes --- include/libcloudph++/blk_1m/rhs_cellwise.hpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/include/libcloudph++/blk_1m/rhs_cellwise.hpp b/include/libcloudph++/blk_1m/rhs_cellwise.hpp index e0628597..5b5dd172 100644 --- a/include/libcloudph++/blk_1m/rhs_cellwise.hpp +++ b/include/libcloudph++/blk_1m/rhs_cellwise.hpp @@ -119,7 +119,7 @@ namespace libcloudphxx ); // limiting - rr_to_rv = std::min(rr, rr_to_rv) / dt; + rr_to_rv = std::min(rr/dt, rr_to_rv); dot_rv += rr_to_rv; dot_rr -= rr_to_rv; @@ -350,13 +350,13 @@ namespace libcloudphxx } //limiting - rv_to_ria = std::min(rv, rv_to_ria) / dt; - rv_to_rib = std::min(rv, rv_to_rib) / dt; - rc_to_ria = std::min(rc, rc_to_ria) / dt; - rc_to_rib = std::min(rc, rc_to_rib) / dt; - rr_to_rib = std::min(rr, rr_to_rib) / dt; - ria_to_rib = std::min(ria, ria_to_rib) / dt; - ria_to_rr = std::min(ria, ria_to_rr) / dt; + rv_to_ria = std::min(rv / dt, rv_to_ria); + rv_to_rib = std::min(rv / dt, rv_to_rib); + rc_to_ria = std::min(rc / dt, rc_to_ria); + rc_to_rib = std::min(rc / dt, rc_to_rib); + rr_to_rib = std::min(rr / dt, rr_to_rib); + ria_to_rib = std::min(ria / dt, ria_to_rib); + ria_to_rr = std::min(ria / dt, ria_to_rr); dot_rc += - rc_to_ria - rc_to_rib; dot_rv += - rv_to_ria - rv_to_rib; From 240a381608e0e2844bd63326870e4a13abfb792d Mon Sep 17 00:00:00 2001 From: pdziekan Date: Fri, 22 Nov 2024 15:36:08 +0100 Subject: [PATCH 15/17] ice type enum --- include/libcloudph++/blk_1m/rhs_columnwise.hpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/include/libcloudph++/blk_1m/rhs_columnwise.hpp b/include/libcloudph++/blk_1m/rhs_columnwise.hpp index 6ab16ab0..db8db3ee 100644 --- a/include/libcloudph++/blk_1m/rhs_columnwise.hpp +++ b/include/libcloudph++/blk_1m/rhs_columnwise.hpp @@ -14,6 +14,8 @@ namespace libcloudphxx { namespace blk_1m { + enum class ice_t {iceA, iceB}; + // expects the arguments to be columns with begin() pointing to the lowest level // returns rain flux out of the domain //@@ -96,7 +98,7 @@ namespace libcloudphxx const cont_t &rhod_cont, const cont_t &ri_cont, const real_t &dz, - const std::string& ice_type + const ice_t& ice_type ) // { @@ -126,7 +128,7 @@ namespace libcloudphxx if (dot_ri != NULL) // i.e. all but first (top) grid cell { flux_t flux_out = real_t(0) * si::kilograms / si::cubic_metres / si::seconds; - if (ice_type == "iceA") + if (ice_type == ice_t::iceA) { // terminal momenta at grid-cell edge (to assure precip mass conservation) flux_out += -real_t(.5) * ( // averaging + axis orientation @@ -140,7 +142,7 @@ namespace libcloudphxx ) ) * (*ri * si::kilograms / si::kilograms) / (dz * si::metres); } - else if (ice_type == "iceB") + else if (ice_type == ice_t::iceB) { // terminal momenta at grid-cell edge (to assure precip mass conservation) flux_out += -real_t(.5) * ( // averaging + axis orientation @@ -166,14 +168,14 @@ namespace libcloudphxx // the bottom grid cell (with mid-cell vterm approximation) flux_t flux_out = real_t(0) * si::kilograms / si::cubic_metres / si::seconds; - if (ice_type == "iceA") + if (ice_type == ice_t::iceA) { flux_out += - (*rhod * si::kilograms / si::cubic_metres) * formulae::velocity_iceA( *ri * si::kilograms / si::kilograms, *rhod * si::kilograms / si::cubic_metres ) * (*ri * si::kilograms / si::kilograms) / (dz * si::metres); } - else if (ice_type == "iceB") + else if (ice_type == ice_t::iceB) { flux_out += - (*rhod * si::kilograms / si::cubic_metres) * formulae::velocity_iceB( *ri * si::kilograms / si::kilograms, From f95a6f01f69684db2b91e00118ca1366561a63c8 Mon Sep 17 00:00:00 2001 From: Agnieszka Makulska Date: Thu, 28 Nov 2024 16:46:48 +0100 Subject: [PATCH 16/17] fixing problem with enum --- bindings/python/blk_1m.hpp | 2 +- bindings/python/lib.cpp | 4 ++++ tests/python/unit/api_blk_1m.py | 5 +++-- 3 files changed, 8 insertions(+), 3 deletions(-) diff --git a/bindings/python/blk_1m.hpp b/bindings/python/blk_1m.hpp index 3f6c2dab..163fba10 100644 --- a/bindings/python/blk_1m.hpp +++ b/bindings/python/blk_1m.hpp @@ -224,7 +224,7 @@ namespace libcloudphxx { const bp_array &rhod, const bp_array &ri, const typename arr_t::T_numtype &dz, - const std::string ice_type + const b1m::ice_t ice_type ) { arr_t np2bz_dot_ri(np2bz(dot_ri)); diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp index 4bc0e9c0..e2e6b1dd 100644 --- a/bindings/python/lib.cpp +++ b/bindings/python/lib.cpp @@ -177,6 +177,10 @@ BOOST_PYTHON_MODULE(libcloudphxx) bp::def("rhs_cellwise_nwtrph_ice", blk_1m::rhs_cellwise_nwtrph_ice); bp::def("rhs_columnwise", blk_1m::rhs_columnwise); // TODO: handle the returned flux bp::def("rhs_columnwise_ice", blk_1m::rhs_columnwise_ice); + + bp::enum_("ice_t") + .value("iceA", b1m::ice_t::iceA) + .value("iceB", b1m::ice_t::iceB); } // blk_2m stuff diff --git a/tests/python/unit/api_blk_1m.py b/tests/python/unit/api_blk_1m.py index 4476e062..190d1759 100644 --- a/tests/python/unit/api_blk_1m.py +++ b/tests/python/unit/api_blk_1m.py @@ -6,6 +6,7 @@ from libcloudphxx import blk_1m opts = blk_1m.opts_t() +ice_t = blk_1m.ice_t() print("cond =", opts.cond) print("cevp =", opts.cevp) print("revp =", opts.revp) @@ -104,7 +105,7 @@ assert dot_rib != 0 #testing sedimentation of ice -flux_iceA = blk_1m.rhs_columnwise_ice(opts, dot_ria, rhod, ria, dz, "iceA") -flux_iceB = blk_1m.rhs_columnwise_ice(opts, dot_rib, rhod, rib, dz, "iceB") +flux_iceA = blk_1m.rhs_columnwise_ice(opts, dot_ria, rhod, ria, dz, ice_t.iceA) +flux_iceB = blk_1m.rhs_columnwise_ice(opts, dot_rib, rhod, rib, dz, ice_t.iceB) assert flux_iceA != 0 assert flux_iceB != 0 \ No newline at end of file From ccc1216b85098004b493d99d44bdcd723dab5ba0 Mon Sep 17 00:00:00 2001 From: Agnieszka Makulska Date: Fri, 29 Nov 2024 17:25:02 +0100 Subject: [PATCH 17/17] adding missing real_t() --- include/libcloudph++/blk_1m/formulae.hpp | 36 ++++++++++++------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/include/libcloudph++/blk_1m/formulae.hpp b/include/libcloudph++/blk_1m/formulae.hpp index ce12917e..ee6702af 100644 --- a/include/libcloudph++/blk_1m/formulae.hpp +++ b/include/libcloudph++/blk_1m/formulae.hpp @@ -126,17 +126,17 @@ namespace libcloudphxx real_t(-4.99e3) - real_t(4.94e4) * std::log10(real_t(1e3) * IWCS / si::kilograms * si::cubic_meters)) / si::meters; quantity m_as = real_t(2) * pi() * moist_air::rho_i() / ( - std::pow(alpha * si::meters, 3) / si::cubic_meters); + std::pow(alpha * si::meters, real_t(3)) / si::cubic_meters); //calculating the mass of large ice A particle: quantity IWCL = rhod_0 * ria - IWCS; - real_t a_mu = real_t(5.2) + real_t(1.3e-3) * (T / si::kelvin - 273.16); - real_t b_mu = real_t(0.026) - real_t(1.2e-3) * (T / si::kelvin - 273.16); - real_t a_sigma = real_t(0.47) + real_t(2.1e-3) * (T / si::kelvin - 273.16); - real_t b_sigma = real_t(0.018) - real_t(2.1e-4) * (T / si::kelvin - 273.16); + real_t a_mu = real_t(5.2) + real_t(1.3e-3) * (T / si::kelvin - real_t(273.16)); + real_t b_mu = real_t(0.026) - real_t(1.2e-3) * (T / si::kelvin - real_t(273.16)); + real_t a_sigma = real_t(0.47) + real_t(2.1e-3) * (T / si::kelvin - real_t(273.16)); + real_t b_sigma = real_t(0.018) - real_t(2.1e-4) * (T / si::kelvin - real_t(273.16)); real_t mu = a_mu + b_mu * std::log10(real_t(1e3) * IWCL / si::kilograms * si::cubic_meters); real_t sigma = a_sigma + b_sigma * std::log10(real_t(1e3) * IWCL / si::kilograms * si::cubic_meters); quantity m_al = real_t(1.67e17) * pi() * moist_air::rho_i() * std::exp( - 3 * mu + 9 / 2 * std::pow(sigma, 2)) * si::cubic_meters; + real_t(3) * mu + real_t(9) / real_t(2) * std::pow(sigma, real_t(2))) * si::cubic_meters; //mass of the average ice A particle: real_t delta = IWCS / (ria * rhod_0); return delta * m_as + (real_t(1) - delta) * m_al; @@ -278,7 +278,7 @@ namespace libcloudphxx real_t beta = (T > real_t(213.16) * si::kelvins) ? real_t(0.1) + real_t(0.9) * (T - real_t(213.16) * si::kelvins) / (real_t(20) * si::kelvins) : real_t(0.1); - quantity rv_adj = beta * rvs + (1 - beta) * rvsi; + quantity rv_adj = beta * rvs + (real_t(1) - beta) * rvsi; //homogeneous nucleation rate (only if T < -40 C) return (T < real_t(233.16) * si::kelvins) ? (real_t(1) - std::exp(-dt / taunuc)) * std::max(real_t(0) * si::dimensionless(), rv - rv_adj) / dt @@ -335,14 +335,14 @@ namespace libcloudphxx const quantity &rhod_0 ) { using namespace common; - if (ria == 0 || rr == 0) // ice B is formed only if ice A and rain are present + if (ria == real_t(0) || rr == real_t(0)) // ice B is formed only if ice A and rain are present { return real_t(0) / si::seconds; } quantity::type, real_t> lambda_r = lambda_rain(rr, rhod_0); //mean raindrop mass (eq. A2 in Grabowski 1999) quantity m_r = pi() * moist_air::rho_w() / (real_t(6) * std::pow( - lambda_r * si::meters, 3)) * si::cubic_meters; + lambda_r * si::meters, real_t(3))) * si::cubic_meters; //mean raindrop velocity (eq. A3 in Grabowski 1999) real_t v_r = real_t(251) * std::pow(lambda_r * si::meters * rhod_0 / si::kilograms * si::cubic_meters, real_t(-0.5)); @@ -370,7 +370,7 @@ namespace libcloudphxx const quantity &rhod_0 ) { using namespace common; - if (ria == 0 || rr == 0) // ice B is formed only if ice A and rain are present + if (ria == real_t(0) || rr == real_t(0)) // ice B is formed only if ice A and rain are present { return real_t(0) / si::seconds; } @@ -401,7 +401,7 @@ namespace libcloudphxx const quantity &rhod_0 ) { using namespace common; - if (ria == 0 || T < 273.16 * si::kelvin) // calculating melting rate only if ice A is present and T > 0 C + if (ria == real_t(0) || T < real_t(273.16) * si::kelvin) // calculating melting rate only if ice A is present and T > 0 C { return real_t(0) / si::seconds; } @@ -416,7 +416,7 @@ namespace libcloudphxx //ventilation factor: real_t F_a = real_t(0.78) + real_t(0.27) * std::pow(Re, real_t(0.5)); return real_t(4.5e-7) * ria / - m_a * D_a * (273.16 - T / si::kelvins) * F_a * si::kilograms / si::meters / si::seconds; + m_a * D_a * (real_t(273.16) - T / si::kelvins) * F_a * si::kilograms / si::meters / si::seconds; } //melting of ice B (rib -> rr) @@ -428,7 +428,7 @@ namespace libcloudphxx const quantity &rhod_0 ) { using namespace common; - if (rib == 0 || T < 273.16 * si::kelvin) // calculating melting rate only if ice B is present and T > 0 C + if (rib == real_t(0) || T < real_t(273.16) * si::kelvin) // calculating melting rate only if ice B is present and T > 0 C { return real_t(0) / si::seconds; } @@ -445,7 +445,7 @@ namespace libcloudphxx //ventilation factor: real_t F_b = real_t(0.78) + real_t(0.27) * std::pow(Re, real_t(0.5)); return real_t(4.5e-7) * rib / - m_b * D_b * (273.16 - T / si::kelvins) * F_b * si::kilograms / si::meters / si::seconds; + m_b * D_b * (real_t(273.16) - T / si::kelvins) * F_b * si::kilograms / si::meters / si::seconds; } @@ -460,7 +460,7 @@ namespace libcloudphxx const quantity &T, const quantity &rhod_0 ) { - if (ria == 0 || T > 273.16 * si::kelvins) + if (ria == real_t(0) || T > real_t(273.16) * si::kelvins) { return real_t(0) / si::seconds; } @@ -487,7 +487,7 @@ namespace libcloudphxx const quantity &T, const quantity &rhod_0 ) { - if (ria == 0 || T > 273.16 * si::kelvins) + if (ria == real_t(0) || T > real_t(273.16) * si::kelvins) { return real_t(0) / si::seconds; } @@ -534,7 +534,7 @@ namespace libcloudphxx const quantity &T, const quantity &rhod_0 ) { - if (rib == 0 || T > 273.16 * si::kelvins) // calculating growth rate only if ice A is present + if (rib == real_t(0) || T > real_t(273.16) * si::kelvins) // calculating growth rate only if ice A is present { return real_t(0) / si::seconds; } @@ -560,7 +560,7 @@ namespace libcloudphxx const quantity &T, const quantity &rhod_0 ) { - if (rib == 0 || T > 273.16 * si::kelvins) + if (rib == real_t(0) || T > real_t(273.16) * si::kelvins) { return real_t(0) / si::seconds; }