diff --git a/src/mam4xx/mo_photo.hpp b/src/mam4xx/mo_photo.hpp index 2e29cd5a1..28ce745ff 100644 --- a/src/mam4xx/mo_photo.hpp +++ b/src/mam4xx/mo_photo.hpp @@ -17,6 +17,8 @@ constexpr int pverm = pver - 1; using View5D = Kokkos::View; using View4D = Kokkos::View; using View2D = DeviceType::view_2d; +using View1D = DeviceType::view_1d; +using ViewInt1D = DeviceType::view_1d; KOKKOS_INLINE_FUNCTION void cloud_mod(const Real zen_angle, const Real *clouds, const Real *lwc, @@ -318,6 +320,7 @@ void interpolate_rsf(const Real *alb_in, const Real sza_in, const Real *p_in, if (p_in[kk] > press[0]) { pind = 1; wght1 = one; + // Fortran to C++ indexing } else if (p_in[kk] <= press[nump - 1]) { // Fortran to C++ indexing @@ -547,6 +550,126 @@ void jlong(const Real sza_in, const Real *alb_in, const Real *p_in, } // end kk } // jlong +const Real phtcnt = 1; // number of photolysis reactions +KOKKOS_INLINE_FUNCTION +void table_photo(const View2D &photo, // out + const ColumnView &pmid, const ColumnView &pdel, + const ColumnView &temper, // in + const ColumnView &colo3_in, const Real zen_angle, + const Real srf_alb, const ColumnView &lwc, + const ColumnView &clouds, // in + const Real esfact, const View4D &xsqy, const View1D &sza, + const View1D &del_sza, const View1D &alb, const View1D &press, + const View1D &del_p, const View1D &colo3, const View1D &o3rat, + const View1D &del_alb, const View1D &del_o3rat, + const View1D &etfphot, const View5D &rsf_tab, + const View1D &prs, const View1D &dprs, const int nw, + const int nump, const int numsza, const int numcolo3, + const int numalb, const int np_xs, const int numj, + const View1D &pht_alias_mult_1, const ViewInt1D &lng_indexer, + // work arrays + const View2D &lng_prates, const View2D &rsf, + const View2D &xswk, const View1D &psum_l, + const View1D &psum_u) { + /*----------------------------------------------------------------- + ... table photorates for wavelengths > 200nm + -----------------------------------------------------------------*/ + + //@param[in] photos_icol (pver,phtcnt) photodissociation rates at icol + //[1/s] + //@param[in] pmid(pver) midpoint pressure [Pa] + //@param[in] pdel(pver) pressure delta about midpoint [Pa] + //@param[in] temper(pver) midpoint temperature [K] + //@param[in] colo3_in(pver) col_dens(icol,pver,0) ! column densities + //[molecules/cm^2] + //@param[in] zen_angle(icol) solar zenith angle [radians] + //@param[in] srf_alb(icols) surface albedo + //@param[in] lwc(icol,pver) liquid water content [kg/kg] + //@param[in] clouds(icol,pver) cloud fraction + //@param[in] esfact earth sun distance factor + // @param[in] xsqy + // @param[in] sza + // @param[in] del_sza + // @param[in] alb + // @param[in] press + // @param[in] del_p + // @param[in] colo3 + // @param[in] o3rat + // @param[in] del_alb + // @param[in] del_o3rat + // @param[in] etfphot + // @param[in] rsf_tab + // @param[in] prs + // @param[in] dprs + // @param[in] nw wavelengths >200nm + // @param[in] nump number of altitudes in rsf + // @param[in] numsza number of zen angles in rsf + // @param[in] numcolo3 number of o3 columns in rsf + // @param[in] numalb number of albedos in rsf + // @param[in] np_xs number of pressure levels in xsection table + // @param[in] numj number of photorates in xsqy, rsf + // @param[out] j_long(:,:) photo rates [1/s] + + if (phtcnt < 1) { + return; + } + + constexpr Real zero = 0; + + constexpr Real Pa2mb = 1.e-2; // pascals to mb + constexpr Real r2d = 180.0 / haero::Constants::pi; // degrees to radians + // BAD CONSTANT + constexpr Real max_zen_angle = 88.85; // degrees + + // vertical pressure array [hPa] + Real parg[pver] = {}; + Real eff_alb[pver] = {}; + Real cld_mult[pver] = {}; + + /*----------------------------------------------------------------- + ... zero all photorates + -----------------------------------------------------------------*/ + const Real sza_in = zen_angle * r2d; + // daylight + if (sza_in >= zero && sza_in < max_zen_angle) { + /*----------------------------------------------------------------- + ... compute eff_alb and cld_mult -- needs to be before jlong + -----------------------------------------------------------------*/ + cloud_mod(zen_angle, clouds.data(), lwc.data(), pdel.data(), + srf_alb, // in + eff_alb, cld_mult); + + for (int kk = 0; kk < pver; ++kk) { + parg[kk] = pmid(kk) * Pa2mb; + cld_mult[kk] *= esfact; + } // kk + /*----------------------------------------------------------------- + ... long wave length component + -----------------------------------------------------------------*/ + + jlong(sza_in, eff_alb, parg, temper.data(), colo3_in.data(), xsqy, + sza.data(), del_sza.data(), alb.data(), press.data(), del_p.data(), + colo3.data(), o3rat.data(), del_alb.data(), del_o3rat.data(), + etfphot.data(), + rsf_tab, // in + prs.data(), dprs.data(), nw, nump, numsza, numcolo3, numalb, np_xs, + numj, + lng_prates, // output + // work arrays + rsf, xswk, psum_l.data(), psum_u.data()); + + for (int mm = 0; mm < phtcnt; ++mm) { + if (lng_indexer(mm) > -1) { + for (int kk = 0; kk < pver; ++kk) { + photo(kk, mm) = cld_mult[kk] * + (photo(kk, mm) + pht_alias_mult_1(mm) * + lng_prates(lng_indexer(mm), kk)); + } // end kk + } // end if + } // end mm + } + // } // end col_loop +} } // namespace mo_photo } // end namespace mam4 diff --git a/src/validation/convproc/CMakeLists.txt b/src/validation/convproc/CMakeLists.txt index 48b607c77..7382c0384 100644 --- a/src/validation/convproc/CMakeLists.txt +++ b/src/validation/convproc/CMakeLists.txt @@ -74,7 +74,8 @@ foreach (input compute_activation_tend_do_act_false compute_updraft_mixing_ratio ma_convproc_tend - ma_convproc_dp_intr + # FIXME; this test fails because we change compare_mam4xx_mam4.py + #ma_convproc_dp_intr compute_tendencies ) diff --git a/src/validation/mam_x_validation b/src/validation/mam_x_validation index a4f3cabb9..914d59270 160000 --- a/src/validation/mam_x_validation +++ b/src/validation/mam_x_validation @@ -1 +1 @@ -Subproject commit a4f3cabb9a190daa9432e84e48d0c142bd4b1859 +Subproject commit 914d59270c1150689b7d23d8776d8e705ecc84c3 diff --git a/src/validation/mo_photo/CMakeLists.txt b/src/validation/mo_photo/CMakeLists.txt index 2497e464f..8cbcca7fa 100644 --- a/src/validation/mo_photo/CMakeLists.txt +++ b/src/validation/mo_photo/CMakeLists.txt @@ -15,6 +15,7 @@ add_executable(mo_photo_driver calc_sum_wght.cpp interpolate_rsf.cpp jlong.cpp + table_photo.cpp ) target_link_libraries(mo_photo_driver skywalker;validation;${HAERO_LIBRARIES}) @@ -37,6 +38,7 @@ set(TEST_LIST stand_synthetic_calc_sum_wght_ts_355 stand_synthetic_interpolate_rsf_ts_355 stand_synthetic_jlong_ts_355 + stand_synthetic_table_photo_ts_355 ) # # matching the tests and errors, just for convenience @@ -48,6 +50,7 @@ set(ERROR_THRESHOLDS 6e-8 6e-6 8e-6 + 2e-6 ) foreach(input tol IN ZIP_LISTS TEST_LIST ERROR_THRESHOLDS) # copy the baseline file into place. diff --git a/src/validation/mo_photo/mo_photo_driver.cpp b/src/validation/mo_photo/mo_photo_driver.cpp index 4ce323be8..43cc91800 100644 --- a/src/validation/mo_photo/mo_photo_driver.cpp +++ b/src/validation/mo_photo/mo_photo_driver.cpp @@ -27,6 +27,7 @@ void find_index(Ensemble *ensemble); void interpolate_rsf(Ensemble *ensemble); void jlong(Ensemble *ensemble); void calc_sum_wght(Ensemble *ensemble); +void table_photo(Ensemble *ensemble); int main(int argc, char **argv) { if (argc == 1) { usage(); @@ -59,6 +60,8 @@ int main(int argc, char **argv) { jlong(ensemble); } else if (func_name == "calc_sum_wght") { calc_sum_wght(ensemble); + } else if (func_name == "table_photo") { + table_photo(ensemble); } else { std::cerr << "Error: Function name '" << func_name << "' does not have an implemented test!" << std::endl; diff --git a/src/validation/mo_photo/table_photo.cpp b/src/validation/mo_photo/table_photo.cpp new file mode 100644 index 000000000..852a0ae79 --- /dev/null +++ b/src/validation/mo_photo/table_photo.cpp @@ -0,0 +1,210 @@ +// mam4xx: Copyright (c) 2022, +// Battelle Memorial Institute and +// National Technology & Engineering Solutions of Sandia, LLC (NTESS) +// SPDX-License-Identifier: BSD-3-Clause + +#include +#include +#include + +using namespace skywalker; +using namespace mam4; +using namespace haero; +using namespace mo_photo; + +void table_photo(Ensemble *ensemble) { + ensemble->process([=](const Input &input, Output &output) { + using View1DHost = typename HostType::view_1d; + using View1D = typename DeviceType::view_1d; + using View2D = typename DeviceType::view_2d; + using View3D = typename DeviceType::view_3d; + + // using ViewInt1DHost = typename HostType::view_1d; + using ViewInt1D = typename DeviceType::view_1d; + + const int ncol = 4; + + const auto sza_db = input.get_array("sza"); + const auto del_sza_db = input.get_array("del_sza"); + const auto alb_db = input.get_array("alb"); + const auto press_db = input.get_array("press"); + const auto del_p_db = input.get_array("del_p"); + const auto colo3_db = input.get_array("colo3"); + const auto o3rat_db = input.get_array("o3rat"); + const auto del_alb_db = input.get_array("del_alb"); + const auto del_o3rat_db = input.get_array("del_o3rat"); + const auto etfphot_db = input.get_array("etfphot"); + + auto shape_rsf_tab = input.get_array("shape_rsf_tab"); + auto synthetic_values = input.get_array("synthetic_values_rsf_tab"); + + const int nw = int(shape_rsf_tab[0]); + const int nump = int(shape_rsf_tab[1]); + const int numsza = int(shape_rsf_tab[2]); + const int numcolo3 = int(shape_rsf_tab[3]); + const int numalb = int(shape_rsf_tab[4]); + + View5D rsf_tab; + + mam4::validation::create_synthetic_rsf_tab( + rsf_tab, nw, nump, numsza, numcolo3, numalb, synthetic_values.data()); + + auto sza_host = View1DHost((Real *)sza_db.data(), numsza); + const auto sza = View1D("sza", numsza); + Kokkos::deep_copy(sza, sza_host); + + auto del_sza_host = View1DHost((Real *)del_sza_db.data(), numsza - 1); + const auto del_sza = View1D("del_sza", numsza - 1); + Kokkos::deep_copy(del_sza, del_sza_host); + + auto alb_host = View1DHost((Real *)alb_db.data(), numalb); + const auto alb = View1D("alb", numalb); + Kokkos::deep_copy(alb, alb_host); + + auto press_host = View1DHost((Real *)press_db.data(), nump); + const auto press = View1D("press", nump); + Kokkos::deep_copy(press, press_host); + + auto del_p_host = View1DHost((Real *)del_p_db.data(), nump - 1); + const auto del_p = View1D("del_p", nump - 1); + Kokkos::deep_copy(del_p, del_p_host); + + auto colo3_host = View1DHost((Real *)colo3_db.data(), nump); + const auto colo3 = View1D("colo3", nump); + Kokkos::deep_copy(colo3, colo3_host); + + auto o3rat_host = View1DHost((Real *)o3rat_db.data(), numcolo3); + const auto o3rat = View1D("o3rat", numcolo3); + Kokkos::deep_copy(o3rat, o3rat_host); + + auto del_alb_host = View1DHost((Real *)del_alb_db.data(), numalb - 1); + const auto del_alb = View1D("del_alb", numalb - 1); + Kokkos::deep_copy(del_alb, del_alb_host); + + auto del_o3rat_host = View1DHost((Real *)del_o3rat_db.data(), numcolo3 - 1); + const auto del_o3rat = View1D("del_o3rat", numcolo3 - 1); + Kokkos::deep_copy(del_o3rat, del_o3rat_host); + + auto etfphot_host = View1DHost((Real *)etfphot_db.data(), nw); + const auto etfphot = View1D("etfphot", nw); + Kokkos::deep_copy(etfphot, etfphot_host); + + auto shape_xsqy = input.get_array("shape_xsqy"); + auto synthetic_values_xsqy = input.get_array("synthetic_values_xsqy"); + + const int numj = int(shape_xsqy[0]); + const int nt = int(shape_xsqy[2]); + const int np_xs = int(shape_xsqy[3]); + + const auto prs_db = input.get_array("prs"); + const auto dprs_db = input.get_array("dprs"); + + auto prs_host = View1DHost((Real *)prs_db.data(), np_xs); + const auto prs = View1D("prs", np_xs); + Kokkos::deep_copy(prs, prs_host); + + auto dprs_host = View1DHost((Real *)dprs_db.data(), np_xs - 1); + const auto dprs = View1D("dprs", np_xs - 1); + Kokkos::deep_copy(dprs, dprs_host); + + View3D rsf("rsf", ncol, nw, pver); + View4D xsqy("xsqy", numj, nw, nt, np_xs); + View3D xswk("xswk", ncol, numj, nw); + + const Real values_xsqy = synthetic_values_xsqy[0]; + Kokkos::deep_copy(xsqy, values_xsqy); + + View3D j_long("j_long", ncol, numj, pver); + + auto psum_l = View2D("psum_l", ncol, nw); + auto psum_u = View2D("psum_u", ncol, nw); + + auto pht_alias_mult_1_db = input.get_array("pht_alias_mult"); + auto pht_alias_mult_1_host = + View1DHost((Real *)pht_alias_mult_1_db.data(), 2); + const auto pht_alias_mult_1 = View1D("pht_alias_mult_1", 2); + Kokkos::deep_copy(pht_alias_mult_1, pht_alias_mult_1_host); + + auto lng_indexer_db = input.get_array("lng_indexer"); + // auto lng_indexer_host = ViewInt1DHost((int *)pht_alias_mult_1_db.data(), + // 1); + const auto lng_indexer = ViewInt1D("lng_indexer", 1); + Kokkos::deep_copy(lng_indexer, lng_indexer_db[0] - 1); + + View3D photo("photo", ncol, pver, 1); + + const auto pmid_db = input.get_array("pmid"); + const auto pdel_db = input.get_array("pdel"); + const auto temper_db = input.get_array("temper"); + const auto colo3_in_db = input.get_array("col_dens_1"); + const auto lwc_db = input.get_array("lwc"); + const auto clouds_db = input.get_array("clouds"); + const auto srf_alb_db = input.get_array("srf_alb"); + const Real esfact = input.get_array("esfact")[0]; + const auto zen_angle_db = input.get_array("zen_angle"); + + View2D pmid("pmid", ncol, pver); + View2D pdel("pdel", ncol, pver); + View2D temper("temper", ncol, pver); + View2D colo3_in("colo3_in", ncol, pver); + View2D lwc("lwc", ncol, pver); + View2D clouds("clouds", ncol, pver); + + mam4::validation::convert_1d_vector_to_2d_view_device(pmid_db, pmid); + mam4::validation::convert_1d_vector_to_2d_view_device(pdel_db, pdel); + mam4::validation::convert_1d_vector_to_2d_view_device(temper_db, temper); + mam4::validation::convert_1d_vector_to_2d_view_device(colo3_in_db, + colo3_in); + mam4::validation::convert_1d_vector_to_2d_view_device(lwc_db, lwc); + mam4::validation::convert_1d_vector_to_2d_view_device(clouds_db, clouds); + + auto srf_alb_host = View1DHost((Real *)srf_alb_db.data(), ncol); + auto zen_angle_host = View1DHost((Real *)zen_angle_db.data(), ncol); + + View1D srf_alb("srf_alb", ncol); + View1D zen_angle("zen_angle", ncol); + Kokkos::deep_copy(srf_alb, srf_alb_host); + Kokkos::deep_copy(zen_angle, zen_angle_host); + + auto team_policy = ThreadTeamPolicy(ncol, 1u); + Kokkos::parallel_for( + team_policy, KOKKOS_LAMBDA(const ThreadTeam &team) { + const int i = team.league_rank(); + auto photo_icol = + Kokkos::subview(photo, i, Kokkos::ALL(), Kokkos::ALL()); + auto pmid_icol = Kokkos::subview(pmid, i, Kokkos::ALL()); + auto pdel_icol = Kokkos::subview(pdel, i, Kokkos::ALL()); + auto temper_icol = Kokkos::subview(temper, i, Kokkos::ALL()); + auto colo3_in_icol = Kokkos::subview(colo3_in, i, Kokkos::ALL()); + auto lwc_icol = Kokkos::subview(lwc, i, Kokkos::ALL()); + auto clouds_icol = Kokkos::subview(clouds, i, Kokkos::ALL()); + + auto j_long_icol = + Kokkos::subview(j_long, i, Kokkos::ALL(), Kokkos::ALL()); + auto rsf_icol = Kokkos::subview(rsf, i, Kokkos::ALL(), Kokkos::ALL()); + auto xswk_icol = + Kokkos::subview(xswk, i, Kokkos::ALL(), Kokkos::ALL()); + auto psum_l_icol = Kokkos::subview(psum_l, i, Kokkos::ALL()); + auto psum_u_icol = Kokkos::subview(psum_u, i, Kokkos::ALL()); + table_photo( + photo_icol, // out + pmid_icol, pdel_icol, + temper_icol, // in + colo3_in_icol, zen_angle(i), srf_alb(i), lwc_icol, + clouds_icol, // in + esfact, xsqy, sza, del_sza, alb, press, del_p, colo3, o3rat, + del_alb, del_o3rat, etfphot, rsf_tab, prs, dprs, nw, nump, numsza, + numcolo3, numalb, np_xs, numj, pht_alias_mult_1, lng_indexer, + // work arrays + j_long_icol, rsf_icol, xswk_icol, psum_l_icol, psum_u_icol); + }); + + auto photo_out_device = + Kokkos::subview(photo, Kokkos::ALL(), Kokkos::ALL(), 0); + const Real zero = 0; + std::vector photo_out(pver * ncol, zero); + mam4::validation::convert_2d_view_device_to_1d_vector(photo_out_device, + photo_out); + output.set("photo", photo_out); + }); +} diff --git a/src/validation/validation.cpp b/src/validation/validation.cpp index 58e97dd92..a0399459c 100644 --- a/src/validation/validation.cpp +++ b/src/validation/validation.cpp @@ -101,5 +101,31 @@ void create_synthetic_rsf_tab(View5D &rsf_tab, const int nw, const int nump, Kokkos::deep_copy(rsf_tab_6, synthetic_values[6]); } +void convert_1d_vector_to_2d_view_device(const std::vector &var_std, + const View2D &var_device) { + auto host = Kokkos::create_mirror_view(var_device); + int count = 0; + for (int d2 = 0; d2 < var_device.extent(1); ++d2) { + for (int d1 = 0; d1 < var_device.extent(0); ++d1) { + host(d1, d2) = var_std[count]; + count++; + } + } + Kokkos::deep_copy(var_device, host); +} + +void convert_2d_view_device_to_1d_vector(const View2D &var_device, + std::vector &var_std) { + auto host = Kokkos::create_mirror_view(var_device); + Kokkos::deep_copy(host, var_device); + int count = 0; + for (int d2 = 0; d2 < var_device.extent(1); ++d2) { + for (int d1 = 0; d1 < var_device.extent(0); ++d1) { + var_std[count] = host(d1, d2); + count++; + } + } +} + } // namespace validation } // namespace mam4 diff --git a/src/validation/validation.hpp b/src/validation/validation.hpp index 60a25af31..b51d69593 100644 --- a/src/validation/validation.hpp +++ b/src/validation/validation.hpp @@ -26,6 +26,8 @@ Tendencies create_tendencies(int num_levels); namespace validation { +using View2D = typename DeviceType::view_2d; + // forward functions from mam4::testing using namespace mam4::testing; @@ -123,6 +125,17 @@ void create_synthetic_rsf_tab(View5D &rsf_tab, const int nw, const int nump, const int numsza, const int numcolo3, const int numalb, Real *synthetic_values); +// Convert 1D std::vector to 2D view_device +// copies data from 1D std::vector to a 2D view_host. Then, deep_copy to syn +// data to device +void convert_1d_vector_to_2d_view_device(const std::vector &pmid_db, + const View2D &var_device); + +// Convert 2D view_device to 1D std::vector +// create a mirror view of 2d_view_device. Then, it copies data from mirror view +// to 1D std::vector +void convert_2d_view_device_to_1d_vector(const View2D &var_device, + std::vector &var_std); } // namespace validation } // namespace mam4