diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index a1ae404800..69640288c3 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -3177,6 +3177,9 @@ List of Parameters +---------------------------------------+----------------------------------------------------------+--------------------+------------------------------------+ | **erf.rad_do_subcol_sampling** | Enable MCICA subcolumn sampling | Boolean | true | +---------------------------------------+----------------------------------------------------------+--------------------+------------------------------------+ +| **erf.rad_use_shoc_cldfrac** | Use Native SHOC diagnosed liquid-cloud fraction in | Boolean | true | +| | RRTMGP when runtime Native SHOC is active | | | ++---------------------------------------+----------------------------------------------------------+--------------------+------------------------------------+ | **erf.rad_orbital_year** | Fixed orbital year for zenith calcs | Integer | < 0 uses timestamp year | +---------------------------------------+----------------------------------------------------------+--------------------+------------------------------------+ | **erf.rad_orbital_eccentricity** | Override orbital eccentricity | Real | < 0 uses computed value | @@ -3226,6 +3229,34 @@ List of Parameters | **erf.four_stream_radiation** | use the four-stream radiation approximation | Boolean | false | +---------------------------------------+----------------------------------------------------------+--------------------+------------------------------------+ +Native SHOC cloud coupling +~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +When ``erf.radiation_model = RRTMGP`` and runtime Native SHOC is selected, +``erf.rad_use_shoc_cldfrac`` defaults to true and uses the Native SHOC +``shoc_cldfrac`` diagnostic. This diagnostic is the liquid-cloud fraction from +the SHOC PDF; it excludes cloud ice. The host ``qc`` field remains the liquid +mixing-ratio input to RRTMGP, and the diagnostic is not a replacement +condensate field. Other radiation models and non-Native PBL paths retain their +existing behavior. + +RRTMGP uses one shared total-cloud mask for the combined liquid and ice optics. +Both liquid and ice paths are therefore normalized with that same total +fraction: the SHOC liquid fraction is used for liquid-only layers, while any +positive ice gives the existing binary ice fraction of one. The existing +in-cloud mixing-ratio floor of ``1e-4`` and cap of ``0.005`` are unchanged, as +is the single conversion from kg/m2 to g/m2. This shared-mask approximation +cannot represent liquid occupying only part of an ice-cloud layer. + +With ``erf.rad_use_shoc_cldfrac = false``, cloud fractions are binary, but +Native SHOC + RRTMGP still runs radiation after the post-SHOC state update; +this is a fraction-only comparison, not a restoration of the pre-coupling +algorithm. With ``erf.rad_do_subcol_sampling = true`` RRTMGP uses the existing +MCICA maximum-random-overlap path. Setting it to false uses a deterministic +clear/cloudy band-to-g-point mapping for binary fractions. The combination of +Native SHOC, ``rad_use_shoc_cldfrac = true``, and sampling disabled is rejected +at startup; use sampling, or disable the SHOC fraction for a binary experiment. + .. _inputs-notes: Notes diff --git a/Docs/sphinx_doc/plotfiles/Plotfile3DReference.rst b/Docs/sphinx_doc/plotfiles/Plotfile3DReference.rst index ccfe5777d6..44b0a235a6 100644 --- a/Docs/sphinx_doc/plotfiles/Plotfile3DReference.rst +++ b/Docs/sphinx_doc/plotfiles/Plotfile3DReference.rst @@ -739,7 +739,9 @@ available. - Native SHOC ``pblh`` is reported in metres above local ground (AGL). * - ``shoc_cldfrac`` - 1 - - Subgrid cloud fraction diagnosed by the native SHOC PDF. + - Liquid-cloud fraction diagnosed by the native SHOC PDF. It excludes + cloud ice and is the optional liquid-fraction input to Native-SHOC + RRTMGP coupling. * - ``shoc_ql`` - kg/kg - Cloud-liquid mixing ratio diagnosed by the native SHOC PDF. diff --git a/Docs/sphinx_doc/theory/PBLschemes.rst b/Docs/sphinx_doc/theory/PBLschemes.rst index 6803857d71..366be11258 100644 --- a/Docs/sphinx_doc/theory/PBLschemes.rst +++ b/Docs/sphinx_doc/theory/PBLschemes.rst @@ -301,6 +301,35 @@ below. Native ``state_update`` currently rejects moisture layouts containing cloud-water or cloud-ice number concentration components because a compatible number closure has not yet been implemented. +Native SHOC and RRTMGP cloud coupling +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +For runtime Native SHOC with RRTMGP, ERF passes the same-step +``shoc_cldfrac`` field to radiation after the Native SHOC state update and its +required state refill. ``shoc_cldfrac`` is the SHOC PDF's liquid-cloud +fraction only; it excludes cloud ice. The host ``qc`` remains the liquid +condensate input to RRTMGP, and ``shoc_ql`` is not substituted for it. + +With the default ``erf.rad_use_shoc_cldfrac = true``, liquid-only layers use +the diagnosed fraction and cloud ice retains its binary fraction. RRTMGP +combines the phases under one total-cloud mask and normalizes both phases with +that same mask. The existing ``1e-4`` fraction floor and ``0.005`` in-cloud +mixing-ratio cap remain active. This bounded shared-mask approximation does +not represent liquid occupying only part of an ice-cloud layer; any positive +ice invokes the binary-ice rule. + +Setting ``erf.rad_use_shoc_cldfrac = false`` selects binary cloud fractions +for the comparison while retaining the post-SHOC radiation ordering. It is a +fraction-only comparison, not a restoration of the pre-PR algorithm. The +``erf.rad_do_subcol_sampling`` combinations are: + +* ``true``: use the existing MCICA maximum-random-overlap sampling path; +* ``false`` with binary fractions: use deterministic clear/cloudy + band-to-g-point mapping; +* ``false`` with Native SHOC fractional coupling: reject the configuration at + startup with guidance to enable sampling or disable + ``rad_use_shoc_cldfrac``. + Native SHOC transport modes ~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/PR3500_FIX_REPORT.md b/PR3500_FIX_REPORT.md new file mode 100644 index 0000000000..9154c0c327 --- /dev/null +++ b/PR3500_FIX_REPORT.md @@ -0,0 +1,136 @@ +# ERF PR3500 Fix Report + +Date: 2026-09-09 +Branch: `Use_SHOC_cldfrac_with_RRTMGP` + +## Merge + +`origin/development` was fetched and merged into this branch. + +- Development tip merged: `489a2aa64b8b49986096a3445285d74bcea73739` +- Pre-merge branch tip: `2515e8db62903e9aa20f28bc389466e9e0883a71` +- Merge base: `0f7b33d46713f464d4eda60453bbca3b60e13d8b` +- Merge commit: `185ae6713` +- Repair commit: the branch `HEAD` containing this report. + +The only textual conflict was `Docs/sphinx_doc/Inputs.rst`. The development +version of the table was retained and the PR3500 `erf.rad_use_shoc_cldfrac` +entry was restored in the merged table. + +## Repair disposition + +### R1: Simple-radiation interface + +`RadiationSimple::Run` now accepts the same optional liquid-cloud-fraction +pointer as the common radiation interface. Simple radiation intentionally +ignores the pointer, preserving its existing physics while restoring the +interface contract. + +### R2: Shared cloud-mass mask + +Both liquid and ice condensate conversions in `Radiation::run_impl` now use +the shared total cloud mask. The existing liquid-only SHOC diagnostic, binary +positive-ice rule, fraction floor, in-cloud mixing-ratio cap, effective radii, +snow handling, and single kg/m2-to-g/m2 conversion are retained. + +The production converter and the regression witness share +`radiation_cloud_mass`. For `q=1e-4`, `rho=1`, and `dz=100`, the witness gives: + +- liquid-only fraction `0.25`: in-cloud mass `0.04`, grid-mean mass `0.01`; +- ice-only or mixed total fraction `1`: in-cloud mass `0.01`, grid-mean mass + `0.01` for each phase; +- fraction floor `1e-4` and `q=1e-8`: in-cloud mass `0.01`, grid-mean mass + `1e-6`; +- cap case `q=0.01`, total fraction `1`: in-cloud and grid-mean mass `0.5`. + +The shared-mask limitation is documented: liquid occupying only part of an +ice-cloud layer cannot be represented by this contract. + +### R3: Sampling flag + +`rad_do_subcol_sampling` now selects the existing MCICA maximum-random-overlap +path when true and a deterministic binary clear/cloudy band-to-g-point map +when false. Both shortwave and longwave paths use the selection. Deterministic +mapping initializes all shortwave optical quantities (`tau`, `ssa`, `g`) and +longwave optical depth for clear and cloudy states. + +Native SHOC plus `rad_use_shoc_cldfrac=true` plus disabled sampling is rejected +at initialization with actionable guidance. Binary cloud fractions remain +available with sampling disabled. + +### R4: Native SHOC ordering and pointer contract + +The merged development ordering is preserved: native SHOC advances first, +updated state and face quantities are refilled, and Native-SHOC RRTMGP receives +the same-step `shoc_cldfrac` diagnostic afterward. Existing optional-pointer +validation continues to check cell centering, component count, BoxArray, and +DistributionMapping compatibility. + +### R5: Documentation + +The Inputs, Plotfile3DReference, PBLschemes, and Native SHOC README documents +now describe the liquid-only diagnostic, host `qc` ownership, shared total +mask, binary ice rule, sampling behavior, and the mixed-phase limitation. + +### R6: Fallback and cadence + +`rad_use_shoc_cldfrac=false` remains a fraction-only binary comparison. It does +not restore pre-coupling radiation ordering. Native SHOC/RRTMGP still runs +radiation after the Native SHOC update and refill. + +### R7: Numeric test portability + +SHOC radiation tests no longer use `EXPECT_DOUBLE_EQ` for `amrex::Real`. +Exact comparisons cast literals to `amrex::Real`; mass-contract checks use a +precision-scaled tolerance. + +### R8: Evidence and review status + +This report records the implementation, build, and test evidence below. The +requested external P3 review was not available in this workspace; no external +review claim is made. Full scientific flux/heating-rate and multi-level +runtime validation remain follow-up work requiring usable RRTMGP coefficient +files and a CUDA device. + +## Verification + +Passed: + +- `bash -n MyBuildcldfrac/cmake.sh` +- `git diff --check` +- Clean requested build using `MyBuildcldfrac/cmake.sh` with its configured + CUDA 12.9, MPI, Kokkos, HDF5/NetCDF, Noah-MP, and RRTMGP options. +- `erf_exec` and `erf_shoc_test` linked successfully. +- `ERF_ShocRadiationCouplingTests.cpp` compiled in the CUDA/unit-test + configuration, including the new shared-mask and sampling tests. +- The focused executable linked successfully using the generated CMake link + script. +- `LD_PRELOAD=/nopt/cuda/12.9/lib64/stubs/libcuda.so ./erf_shoc_unit_tests + --gtest_list_tests` exited 0 and listed the complete SHOC suite, including + `ShocRadiationCloudFraction.*`. + +The focused runtime filter was attempted with: + +```text +LD_PRELOAD=/nopt/cuda/12.9/lib64/stubs/libcuda.so \ + ./erf_shoc_unit_tests --gtest_filter=ShocRadiationCloudFraction.* +``` + +It returned the project skip code 77 because the node reports the CUDA stub +driver and zero devices. Assertions were therefore compiled but not executed +on this host. + +The clean `cmake.sh` build log is: +`MyBuildcldfrac/cmake_build_20260909_203757.log`. + +## Suggested PR text + +Summary: Merge current development into the Native SHOC/RRTMGP branch and +repair PR3500 cloud-fraction coupling. Native SHOC now supplies a same-step +liquid-only PDF diagnostic to RRTMGP after state refill; liquid and ice use one +shared total optical mask; and the sampling flag has explicit MCICA and +deterministic binary behavior with an invalid-combination startup check. + +Tests: clean `MyBuildcldfrac/cmake.sh` build passed; the focused SHOC test +executable compiled and linked; runtime SHOC tests were skipped on this host +because no CUDA device/driver is available. diff --git a/Source/ERF.H b/Source/ERF.H index 3eb0ded76b..00929b0a76 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -595,7 +595,8 @@ public: void advance_radiation (int lev, amrex::MultiFab& cons_in, - const double& dt_advance); + const double& dt_advance, + const amrex::MultiFab* liquid_cloud_fraction = nullptr); #ifdef ERF_USE_EAMXX_SHOC void compute_shoc_tendencies (int lev, diff --git a/Source/PBL/Shoc/README.md b/Source/PBL/Shoc/README.md index 4931d1fc36..166da2340b 100644 --- a/Source/PBL/Shoc/README.md +++ b/Source/PBL/Shoc/README.md @@ -433,6 +433,32 @@ Also run representative non-SHOC regression tests after changing shared coupling microphysics, plotfile, or time-integration paths. Native SHOC must not alter non-SHOC results when it is not selected at runtime. +## RRTMGP radiation coupling + +When Native SHOC is coupled to RRTMGP, ERF passes SHOC's diagnosed liquid-cloud +fraction to radiation when `erf.rad_use_shoc_cldfrac = true` (the default). +Radiation is called after the current Native SHOC update so the diagnostic and +the host thermodynamic state are from the same step. The diagnostic describes +liquid-cloud macrophysics only; cloud ice retains a binary fraction. ERF's host +`qc` field remains the liquid input to RRTMGP, and `shoc_ql` is not substituted +for it by this coupling. + +RRTMGP uses one shared total-cloud mask for the combined liquid and ice optics. +Both phases are normalized with that total fraction, preserving the existing +`1e-4` fraction floor, `0.005` in-cloud mixing-ratio cap, effective radii, +and maximum-random-overlap behavior. Any positive ice invokes the binary-ice +rule. This approximation cannot represent liquid occupying only part of an +ice-cloud layer. + +Setting `erf.rad_use_shoc_cldfrac = false` selects binary fractions while +retaining the post-SHOC radiation ordering; it is a fraction-only comparison, +not a restoration of the pre-PR algorithm. With +`erf.rad_do_subcol_sampling = true`, radiation uses the existing MCICA path. +With sampling disabled, RRTMGP uses deterministic clear/cloudy band-to-g-point +mapping for binary fractions. Native SHOC fractional coupling with sampling +disabled is rejected at startup; enable sampling or disable the SHOC fraction +for a binary experiment. + ## See also User documentation: diff --git a/Source/PhysicsInterfaces/Radiation/ERF_RRTMGP_Interface.H b/Source/PhysicsInterfaces/Radiation/ERF_RRTMGP_Interface.H index d1cd0af421..900875e58f 100644 --- a/Source/PhysicsInterfaces/Radiation/ERF_RRTMGP_Interface.H +++ b/Source/PhysicsInterfaces/Radiation/ERF_RRTMGP_Interface.H @@ -87,6 +87,22 @@ get_subsampled_clouds (const int ncol, real2d_k& cld, real2d_k& p_lay); +optical_props2_t +get_deterministic_clouds (const int ncol, + const int nlay, + const int ngpt, + optical_props2_t& cloud_optics, + gas_optics_t& kdist, + const real2d_k& cld); + +optical_props1_t +get_deterministic_clouds (const int ncol, + const int nlay, + const int ngpt, + optical_props1_t& cloud_optics, + gas_optics_t& kdist, + const real2d_k& cld); + void rrtmgp_initialize (gas_concs_t& gas_concs, const std::string& coefficients_file_sw, const std::string& coefficients_file_lw, @@ -132,7 +148,8 @@ rrtmgp_main (const int ncol, const int nlay, real3d_k& sw_bnd_flux_up , real3d_k& sw_bnd_flux_dn , real3d_k& sw_bnd_flux_dn_dir, real3d_k& lw_bnd_flux_up , real3d_k& lw_bnd_flux_dn, const RealT tsi_scaling, - const bool extra_clnclrsky_diag = false, const bool extra_clnsky_diag = false); + const bool extra_clnclrsky_diag = false, const bool extra_clnsky_diag = false, + const bool do_subcol_sampling = true); void diff --git a/Source/PhysicsInterfaces/Radiation/ERF_RRTMGP_Interface.cpp b/Source/PhysicsInterfaces/Radiation/ERF_RRTMGP_Interface.cpp index ed78143c64..41f00b3bf1 100644 --- a/Source/PhysicsInterfaces/Radiation/ERF_RRTMGP_Interface.cpp +++ b/Source/PhysicsInterfaces/Radiation/ERF_RRTMGP_Interface.cpp @@ -1,3 +1,5 @@ +#include + #include "ERF_RRTMGP_Interface.H" namespace rrtmgp { @@ -220,6 +222,65 @@ get_subsampled_clouds (const int ncol, } +optical_props2_t +get_deterministic_clouds (const int ncol, + const int nlay, + const int ngpt, + optical_props2_t& cloud_optics, + gas_optics_t& kdist, + const real2d_k& cld) +{ + optical_props2_t deterministic_optics; + deterministic_optics.init(kdist.get_band_lims_wavenumber(), + kdist.get_band_lims_gpoint(), + "deterministic_optics"); + deterministic_optics.alloc_2str(ncol, nlay); + + const auto gpoint_bands = kdist.get_gpoint_bands(); + Kokkos::parallel_for(Kokkos::MDRangePolicy>({0, 0, 0}, + {ncol, nlay, ngpt}), + KOKKOS_LAMBDA (int icol, int ilay, int igpt) + { + const auto ibnd = gpoint_bands(igpt); + const bool cloudy = cld(icol,ilay) > zero; + deterministic_optics.tau(icol,ilay,igpt) = + cloudy ? cloud_optics.tau(icol,ilay,ibnd) : zero; + deterministic_optics.ssa(icol,ilay,igpt) = + cloudy ? cloud_optics.ssa(icol,ilay,ibnd) : zero; + deterministic_optics.g(icol,ilay,igpt) = + cloudy ? cloud_optics.g(icol,ilay,ibnd) : zero; + }); + return deterministic_optics; +} + + +optical_props1_t +get_deterministic_clouds (const int ncol, + const int nlay, + const int ngpt, + optical_props1_t& cloud_optics, + gas_optics_t& kdist, + const real2d_k& cld) +{ + optical_props1_t deterministic_optics; + deterministic_optics.init(kdist.get_band_lims_wavenumber(), + kdist.get_band_lims_gpoint(), + "deterministic_optics"); + deterministic_optics.alloc_1scl(ncol, nlay); + + const auto gpoint_bands = kdist.get_gpoint_bands(); + Kokkos::parallel_for(Kokkos::MDRangePolicy>({0, 0, 0}, + {ncol, nlay, ngpt}), + KOKKOS_LAMBDA (int icol, int ilay, int igpt) + { + const auto ibnd = gpoint_bands(igpt); + deterministic_optics.tau(icol,ilay,igpt) = + (cld(icol,ilay) > zero) ? cloud_optics.tau(icol,ilay,ibnd) : zero; + }); + return deterministic_optics; +} + + /* * The following routines provide a simple interface to RRTMGP. These * can be used as-is, but are intended to be wrapped by the SCREAM AD @@ -400,7 +461,8 @@ rrtmgp_main (const int ncol, const int nlay, real3d_k& sw_bnd_flux_up , real3d_k& sw_bnd_flux_dn , real3d_k& sw_bnd_flux_dn_dir, real3d_k& lw_bnd_flux_up , real3d_k& lw_bnd_flux_dn , const RealT tsi_scaling, - const bool extra_clnclrsky_diag, const bool extra_clnsky_diag) + const bool extra_clnclrsky_diag, const bool extra_clnsky_diag, + const bool do_subcol_sampling) { // Setup pointers to RRTMGP SW fluxes fluxes_t fluxes_sw; @@ -463,16 +525,41 @@ rrtmgp_main (const int ncol, const int nlay, optical_props1_t clouds_lw = get_cloud_optics_lw(ncol, nlay, *cloud_optics_lw_k, *k_dist_lw_k, lwp, iwp, rel, rei); - // Do subcolumn sampling to map bands -> gpoints based on cloud fraction and overlap assumption; - // This implements the Monte Carlo Independent Column Approximation by mapping only a single - // subcolumn (cloud state) to each gpoint. + const auto sampling_mode = radiation_cloud_sampling_mode(do_subcol_sampling); + if (sampling_mode == RadiationCloudSamplingMode::DeterministicBinary) { + // The deterministic mode is intentionally binary. Native SHOC's + // fractional mode is rejected during configuration; this guard keeps + // any other fractional supplier from silently becoming a binary run. + auto cldfrac_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), cldfrac); + for (int icol = 0; icol < ncol; ++icol) { + for (int ilay = 0; ilay < nlay; ++ilay) { + const RealT fraction = cldfrac_host(icol,ilay); + if (fraction != zero && fraction != one) { + amrex::Abort("RRTMGP deterministic cloud mapping requires binary " + "total cloud fractions, but a fractional value was supplied."); + } + } + } + } + + // Map bands to gpoints either through the existing MCICA maximum-random + // path or through a deterministic clear/cloudy mapping. Both shortwave + // and longwave use the same dispatch contract. auto nswgpts = k_dist_sw_k->get_ngpt(); - auto clouds_sw_gpt = get_subsampled_clouds(ncol, nlay, nswbands, nswgpts, - clouds_sw, *k_dist_sw_k, cldfrac, p_lay); + auto clouds_sw_gpt = + (sampling_mode == RadiationCloudSamplingMode::MCICA) + ? get_subsampled_clouds(ncol, nlay, nswbands, nswgpts, + clouds_sw, *k_dist_sw_k, cldfrac, p_lay) + : get_deterministic_clouds(ncol, nlay, nswgpts, + clouds_sw, *k_dist_sw_k, cldfrac); // Longwave auto nlwgpts = k_dist_lw_k->get_ngpt(); - auto clouds_lw_gpt = get_subsampled_clouds(ncol, nlay, nlwbands, nlwgpts, - clouds_lw, *k_dist_lw_k, cldfrac, p_lay); + auto clouds_lw_gpt = + (sampling_mode == RadiationCloudSamplingMode::MCICA) + ? get_subsampled_clouds(ncol, nlay, nlwbands, nlwgpts, + clouds_lw, *k_dist_lw_k, cldfrac, p_lay) + : get_deterministic_clouds(ncol, nlay, nlwgpts, + clouds_lw, *k_dist_lw_k, cldfrac); // Do shortwave rrtmgp_sw(ncol, nlay, @@ -1217,4 +1304,3 @@ compute_aerocom_cloudtop (int ncol, int nlay , } } // namespace rrtmgp - diff --git a/Source/PhysicsInterfaces/Radiation/ERF_RRTMGP_Utils.H b/Source/PhysicsInterfaces/Radiation/ERF_RRTMGP_Utils.H index 5a94e6f5f3..0dcba10869 100644 --- a/Source/PhysicsInterfaces/Radiation/ERF_RRTMGP_Utils.H +++ b/Source/PhysicsInterfaces/Radiation/ERF_RRTMGP_Utils.H @@ -4,6 +4,7 @@ #include "rrtmgp_const.h" #include "ERF_Constants.H" +#include "ERF_RadiationCloudFraction.H" namespace rrtmgp { @@ -24,14 +25,9 @@ mixing_ratio_to_cloud_mass (View1 const& mixing_ratio , // NOTE: these thresholds (from E3SM) seem arbitrary, but included here for consistency // This limits in-cloud mixing ratio to amrex::Real(0.005) kg/kg. According to note in cloud_diagnostics // in EAM, this is consistent with limits in MG2. Is this true for P3? - if (cloud_fraction(icol,ilay) > 0) { - // Compute layer-integrated cloud mass (per unit area) - auto incloud_mixing_ratio = std::min(mixing_ratio(icol,ilay) / - std::max(amrex::Real(0.0001), cloud_fraction(icol,ilay)), amrex::Real(0.005)); - cloud_mass(icol,ilay) = incloud_mixing_ratio * rho(icol,ilay) * dz(icol,ilay); - } else { - cloud_mass(icol,ilay) = 0; - } + cloud_mass(icol,ilay) = radiation_cloud_mass(mixing_ratio(icol,ilay), + cloud_fraction(icol,ilay), + rho(icol,ilay), dz(icol,ilay)); }); } diff --git a/Source/PhysicsInterfaces/Radiation/ERF_Radiation.H b/Source/PhysicsInterfaces/Radiation/ERF_Radiation.H index 5a25f241c9..dfa620a499 100644 --- a/Source/PhysicsInterfaces/Radiation/ERF_Radiation.H +++ b/Source/PhysicsInterfaces/Radiation/ERF_Radiation.H @@ -39,6 +39,7 @@ #include #include #include +#include #include @@ -109,6 +110,7 @@ public: const amrex::BoxArray& ba, amrex::Geometry& geom, amrex::MultiFab* cons_in, + const amrex::MultiFab* liquid_cloud_fraction, amrex::iMultiFab* lmask, amrex::MultiFab* t_surf, amrex::Vector& lsm_input_ptrs, @@ -121,7 +123,7 @@ public: const bool updated_lsm) override { set_grids(level, step, time, dt, ba, geom, - cons_in, lmask, t_surf, + cons_in, liquid_cloud_fraction, lmask, t_surf, lsm_input_ptrs, qheating_rates, rad_fluxes, z_phys, lat_ptr, lon_ptr, updated_lsm); @@ -137,6 +139,7 @@ public: const amrex::BoxArray& ba, amrex::Geometry& geom, amrex::MultiFab* cons_in, + const amrex::MultiFab* liquid_cloud_fraction, amrex::iMultiFab* lmask, amrex::MultiFab* t_surf, amrex::Vector& lsm_input_ptrs, @@ -253,6 +256,7 @@ private: bool m_moist = false; bool m_ice = false; int m_qi_comp = -1; + bool m_use_shoc_cldfrac = true; // Do we have a land surface model? bool m_lsm = false; @@ -274,6 +278,9 @@ private: // Pointer to the CC conserved vars amrex::MultiFab* m_cons_in = nullptr; + // Non-owning Native SHOC liquid-cloud fraction input. + const amrex::MultiFab* m_liquid_cloud_fraction = nullptr; + // Pointer to the radiation source terms amrex::MultiFab* m_qheating_rates = nullptr; @@ -426,6 +433,8 @@ private: real2d_k qv_lay; real2d_k qc_lay; real2d_k qi_lay; + real2d_k cldfrac_liq; + real2d_k cldfrac_ice; real2d_k cldfrac_tot; real2d_k eff_radius_qc; real2d_k eff_radius_qi; diff --git a/Source/PhysicsInterfaces/Radiation/ERF_Radiation.cpp b/Source/PhysicsInterfaces/Radiation/ERF_Radiation.cpp index 12dd22a03a..8ace1b0d42 100644 --- a/Source/PhysicsInterfaces/Radiation/ERF_Radiation.cpp +++ b/Source/PhysicsInterfaces/Radiation/ERF_Radiation.cpp @@ -102,6 +102,9 @@ Radiation::Radiation (const int& lev, // Radiation timestep, as a number of atm steps pp.queryAdd("rad_freq_in_steps", m_rad_freq_in_steps); + // Use Native SHOC's diagnosed liquid-cloud fraction when supplied. + pp.query("rad_use_shoc_cldfrac", m_use_shoc_cldfrac); + // Get nvar if specified pp.queryAdd("rad_nvar", m_rad_nvar); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_rad_nvar >= 0, @@ -123,6 +126,15 @@ Radiation::Radiation (const int& lev, // Do MCICA subcolumn sampling pp.queryAdd("rad_do_subcol_sampling", m_do_subcol_sampling); + if (radiation_sampling_config_unsupported( + m_do_subcol_sampling, sc.turbChoice[lev].uses_native_shoc(), + m_use_shoc_cldfrac)) { + amrex::Abort( + "erf.rad_do_subcol_sampling=false does not support Native SHOC " + "fractional-cloud coupling. Set erf.rad_do_subcol_sampling=true, " + "or set erf.rad_use_shoc_cldfrac=false for binary cloud fractions."); + } + // Determine orbital year. If orbital_year is negative, use current year // from timestamp for orbital year; if non-negative, use provided orbital year // for duration of simulation. Note that this is keyed off the value itself, @@ -226,6 +238,9 @@ Radiation::Radiation (const int& lev, << m_nswbands << " " << m_nlwbands << "\n"; Print() << "Number of short/longwave gauss points: " << m_nswgpts << " " << m_nlwgpts << "\n"; + Print() << "Native SHOC cloud fraction for RRTMGP: " + << ((m_use_shoc_cldfrac && sc.turbChoice[lev].uses_native_shoc()) ? "enabled" : "disabled") + << "\n"; Print() << "========================================================\n"; } } @@ -238,6 +253,7 @@ Radiation::set_grids (int& level, const BoxArray& ba, Geometry& geom, MultiFab* cons_in, + const MultiFab* liquid_cloud_fraction, iMultiFab* lmask, MultiFab* t_surf, Vector& lsm_input_ptrs, @@ -256,12 +272,24 @@ Radiation::set_grids (int& level, m_dt = dt; m_geom = geom; m_cons_in = cons_in; + m_liquid_cloud_fraction = liquid_cloud_fraction; m_qheating_rates = qheating_rates; m_rad_fluxes = rad_fluxes; m_z_phys = z_phys; m_lat = lat; m_lon = lon; + if (m_liquid_cloud_fraction != nullptr) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_liquid_cloud_fraction->boxArray().ixType().cellCentered(), + "Radiation liquid cloud fraction must be cell-centered."); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_liquid_cloud_fraction->nComp() > 0, + "Radiation liquid cloud fraction must have at least one component."); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_liquid_cloud_fraction->boxArray() == cons_in->boxArray(), + "Radiation liquid cloud fraction must match the conserved state BoxArray."); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_liquid_cloud_fraction->DistributionMap() == cons_in->DistributionMap(), + "Radiation liquid cloud fraction must match the conserved state DistributionMap."); + } + // Update the day and month time_t timestamp = time_t(time); struct tm *timeinfo = gmtime(×tamp); @@ -359,6 +387,8 @@ Radiation::alloc_buffers () qv_lay = real2d_k("qv" , m_ncol, m_nlay); qc_lay = real2d_k("qc" , m_ncol, m_nlay); qi_lay = real2d_k("qi" , m_ncol, m_nlay); + cldfrac_liq = real2d_k("cldfrac_liq" , m_ncol, m_nlay); + cldfrac_ice = real2d_k("cldfrac_ice" , m_ncol, m_nlay); cldfrac_tot = real2d_k("cldfrac_tot" , m_ncol, m_nlay); eff_radius_qc = real2d_k("eff_radius_qc", m_ncol, m_nlay); eff_radius_qi = real2d_k("eff_radius_qi", m_ncol, m_nlay); @@ -488,6 +518,8 @@ Radiation::dealloc_buffers () qv_lay = real2d_k(); qc_lay = real2d_k(); qi_lay = real2d_k(); + cldfrac_liq = real2d_k(); + cldfrac_ice = real2d_k(); cldfrac_tot = real2d_k(); eff_radius_qc = real2d_k(); eff_radius_qi = real2d_k(); @@ -557,6 +589,8 @@ Radiation::mf_to_kokkos_buffers (iMultiFab* lmask, Table2D qv_lay_tab(qv_lay.data(), {0,0}, {static_cast(qv_lay.extent(0)),static_cast(qv_lay.extent(1))}); Table2D qc_lay_tab(qc_lay.data(), {0,0}, {static_cast(qc_lay.extent(0)),static_cast(qc_lay.extent(1))}); Table2D qi_lay_tab(qi_lay.data(), {0,0}, {static_cast(qi_lay.extent(0)),static_cast(qi_lay.extent(1))}); + Table2D cldfrac_liq_tab(cldfrac_liq.data(), {0,0}, {static_cast(cldfrac_liq.extent(0)),static_cast(cldfrac_liq.extent(1))}); + Table2D cldfrac_ice_tab(cldfrac_ice.data(), {0,0}, {static_cast(cldfrac_ice.extent(0)),static_cast(cldfrac_ice.extent(1))}); Table2D cldfrac_tot_tab(cldfrac_tot.data(), {0,0}, {static_cast(cldfrac_tot.extent(0)),static_cast(cldfrac_tot.extent(1))}); Table2D lwp_tab(lwp.data(), {0,0}, {static_cast(lwp.extent(0)),static_cast(lwp.extent(1))}); @@ -577,6 +611,7 @@ Radiation::mf_to_kokkos_buffers (iMultiFab* lmask, const bool has_lsm = m_lsm; const bool has_lat = m_lat; const bool has_lon = m_lon; + const bool has_shoc_cldfrac = m_use_shoc_cldfrac && (m_liquid_cloud_fraction != nullptr); const bool has_surflayer = (t_surf); int ncol = m_ncol; int nlay = m_nlay; @@ -598,6 +633,8 @@ Radiation::mf_to_kokkos_buffers (iMultiFab* lmask, Array4{}; const Array4& lon_arr = (m_lon) ? m_lon->const_array(mfi) : Array4{}; + const Array4& shoc_cf_arr = has_shoc_cldfrac + ? m_liquid_cloud_fraction->const_array(mfi) : Array4{}; ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // map [i,j,k] 0-based to [icol, ilay] 0-based @@ -610,6 +647,10 @@ Radiation::mf_to_kokkos_buffers (iMultiFab* lmask, Real qv = (moist) ? std::max(cons_arr(i,j,k,RhoQ1_comp)/r,Real(0.)) : Real(0.); Real qc = (moist) ? std::max(cons_arr(i,j,k,RhoQ2_comp)/r,Real(0.)) : Real(0.); Real qi = (ice && qi_comp >= 0) ? std::max(cons_arr(i,j,k,qi_comp)/r,Real(0.)) : Real(0.); + const RadiationCloudFractions cloud_fractions = + radiation_cloud_fractions(qc, qi, + has_shoc_cldfrac ? shoc_cf_arr(i,j,k) : Real(0.), + has_shoc_cldfrac); // EOS avg to z-face Real r_lo = cons_arr(i,j,k-1,Rho_comp); @@ -641,7 +682,9 @@ Radiation::mf_to_kokkos_buffers (iMultiFab* lmask, qv_lay_tab(icol,ilay) = qv; qc_lay_tab(icol,ilay) = qc; qi_lay_tab(icol,ilay) = qi; - cldfrac_tot_tab(icol,ilay) = ((qc+qi)>Real(0.)) ? Real(1.) : Real(0.); + cldfrac_liq_tab(icol,ilay) = cloud_fractions.liquid; + cldfrac_ice_tab(icol,ilay) = cloud_fractions.ice; + cldfrac_tot_tab(icol,ilay) = cloud_fractions.total; // NOTE: These are populated in 'mixing_ratio_to_cloud_mass' lwp_tab(icol,ilay) = Real(0.); @@ -1265,6 +1308,9 @@ Radiation::run_impl () Kokkos::deep_copy(mu0, h_mu0); // Compute layer cloud mass per unit area (populates lwp/iwp) + // Both phases are represented by the same combined optical mask below. Use + // that same fraction for both conversions so the grid-mean condensate mass + // is not inflated when liquid and ice coexist in a layer. rrtmgp::mixing_ratio_to_cloud_mass(qc_lay, cldfrac_tot, r_lay, z_del, lwp); rrtmgp::mixing_ratio_to_cloud_mass(qi_lay, cldfrac_tot, r_lay, z_del, iwp); @@ -1430,7 +1476,8 @@ Radiation::run_impl () lw_clnsky_flux_up_c, lw_clnsky_flux_dn_c, sw_bnd_flux_up_c, sw_bnd_flux_dn_c, sw_bnd_flux_dir_c, lw_bnd_flux_up_c, lw_bnd_flux_dn_c, - eccf, m_extra_clnclrsky_diag, m_extra_clnsky_diag); + eccf, m_extra_clnclrsky_diag, m_extra_clnsky_diag, + m_do_subcol_sampling); // Compute heating rates for this chunk rrtmgp::compute_heating_rate(sw_flux_up_c, sw_flux_dn_c, r_lay_c, z_del_c, sw_heating_c); diff --git a/Source/PhysicsInterfaces/Radiation/ERF_RadiationCloudFraction.H b/Source/PhysicsInterfaces/Radiation/ERF_RadiationCloudFraction.H new file mode 100644 index 0000000000..b384eda6b0 --- /dev/null +++ b/Source/PhysicsInterfaces/Radiation/ERF_RadiationCloudFraction.H @@ -0,0 +1,85 @@ +#ifndef ERF_RADIATION_CLOUD_FRACTION_H +#define ERF_RADIATION_CLOUD_FRACTION_H + +#include +#include + +struct RadiationCloudFractions +{ + amrex::Real liquid; + amrex::Real ice; + amrex::Real total; +}; + +enum class RadiationCloudSamplingMode +{ + MCICA, + DeterministicBinary +}; + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +RadiationCloudSamplingMode +radiation_cloud_sampling_mode (const bool do_subcol_sampling) noexcept +{ + return do_subcol_sampling ? RadiationCloudSamplingMode::MCICA + : RadiationCloudSamplingMode::DeterministicBinary; +} + +inline bool +radiation_sampling_config_unsupported (const bool do_subcol_sampling, + const bool native_shoc, + const bool use_shoc_cldfrac) noexcept +{ + return !do_subcol_sampling && native_shoc && use_shoc_cldfrac; +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real +radiation_cloud_mass (const amrex::Real mixing_ratio, + const amrex::Real cloud_fraction, + const amrex::Real rho, + const amrex::Real dz) noexcept +{ + constexpr amrex::Real min_cloud_fraction = amrex::Real(1.0e-4); + constexpr amrex::Real max_incloud_mixing_ratio = amrex::Real(0.005); + + if (cloud_fraction <= amrex::Real(0.0)) { + return amrex::Real(0.0); + } + + const amrex::Real effective_fraction = + (cloud_fraction < min_cloud_fraction) ? min_cloud_fraction : cloud_fraction; + const amrex::Real in_cloud_mixing_ratio = mixing_ratio / effective_fraction; + const amrex::Real capped_mixing_ratio = + (in_cloud_mixing_ratio < max_incloud_mixing_ratio) + ? in_cloud_mixing_ratio : max_incloud_mixing_ratio; + return capped_mixing_ratio * rho * dz; +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +RadiationCloudFractions +radiation_cloud_fractions (const amrex::Real qc, + const amrex::Real qi, + const amrex::Real supplied_liquid_fraction, + const bool use_supplied_liquid_fraction) noexcept +{ + constexpr amrex::Real min_cloud_fraction = amrex::Real(1.0e-4); + + amrex::Real liquid_cf = amrex::Real(0.0); + if (qc > amrex::Real(0.0)) { + if (use_supplied_liquid_fraction) { + liquid_cf = supplied_liquid_fraction; + liquid_cf = (liquid_cf < amrex::Real(0.0)) ? amrex::Real(0.0) : liquid_cf; + liquid_cf = (liquid_cf > amrex::Real(1.0)) ? amrex::Real(1.0) : liquid_cf; + liquid_cf = (liquid_cf < min_cloud_fraction) ? min_cloud_fraction : liquid_cf; + } else { + liquid_cf = amrex::Real(1.0); + } + } + + const amrex::Real ice_cf = (qi > amrex::Real(0.0)) ? amrex::Real(1.0) : amrex::Real(0.0); + return {liquid_cf, ice_cf, + (liquid_cf > ice_cf) ? liquid_cf : ice_cf}; +} + +#endif diff --git a/Source/PhysicsInterfaces/Radiation/ERF_RadiationInterface.H b/Source/PhysicsInterfaces/Radiation/ERF_RadiationInterface.H index 7fe09c1b49..8d94ee3e31 100644 --- a/Source/PhysicsInterfaces/Radiation/ERF_RadiationInterface.H +++ b/Source/PhysicsInterfaces/Radiation/ERF_RadiationInterface.H @@ -30,6 +30,7 @@ public: const amrex::BoxArray& ba, amrex::Geometry& geom, amrex::MultiFab* cons_in, + const amrex::MultiFab* liquid_cloud_fraction, amrex::iMultiFab* lmask, amrex::MultiFab* t_surf, amrex::Vector& lsm_input_ptrs, diff --git a/Source/PhysicsInterfaces/Radiation/Simple/ERF_RadiationSimple.H b/Source/PhysicsInterfaces/Radiation/Simple/ERF_RadiationSimple.H index 608c6abfa1..adeab61a31 100644 --- a/Source/PhysicsInterfaces/Radiation/Simple/ERF_RadiationSimple.H +++ b/Source/PhysicsInterfaces/Radiation/Simple/ERF_RadiationSimple.H @@ -42,6 +42,7 @@ public: const amrex::BoxArray& ba, amrex::Geometry& geom, amrex::MultiFab* cons_in, + const amrex::MultiFab* liquid_cloud_fraction, amrex::iMultiFab* lmask, amrex::MultiFab* t_surf, amrex::Vector& lsm_input_ptrs, diff --git a/Source/PhysicsInterfaces/Radiation/Simple/ERF_RadiationSimple.cpp b/Source/PhysicsInterfaces/Radiation/Simple/ERF_RadiationSimple.cpp index 43169e7f9e..6958cfa038 100644 --- a/Source/PhysicsInterfaces/Radiation/Simple/ERF_RadiationSimple.cpp +++ b/Source/PhysicsInterfaces/Radiation/Simple/ERF_RadiationSimple.cpp @@ -31,6 +31,7 @@ void RadiationSimple::Run (int& /*level*/, const amrex::BoxArray& /*ba*/, amrex::Geometry& geom, amrex::MultiFab* cons_in, + const amrex::MultiFab* /*liquid_cloud_fraction*/, amrex::iMultiFab* /*lmask*/, amrex::MultiFab* /*t_surf*/, amrex::Vector& /*lsm_input_ptrs*/, diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index 82d2c62920..93d610525d 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -147,7 +147,13 @@ ERF::Advance (int lev, double time, double dt_lev, int iteration, int /*ncycle*/ // ************************************************************************************** // Update the radiation sources with the "old" state // ************************************************************************************** - advance_radiation(lev, S_old, dt_lev); + const bool use_native_shoc_rad_coupling = + solverChoice.rad_type == RadiationType::RRTMGP && + solverChoice.turbChoice[lev].uses_native_shoc(); + + if (!use_native_shoc_rad_coupling) { + advance_radiation(lev, S_old, dt_lev); + } // ************************************************************************************** // Update the "old" state using SHOC @@ -202,6 +208,13 @@ ERF::Advance (int lev, double time, double dt_lev, int iteration, int /*ncycle*/ } } + if (use_native_shoc_rad_coupling) { + AMREX_ALWAYS_ASSERT(native_shoc_driver[lev]); + AMREX_ALWAYS_ASSERT(native_shoc_driver[lev]->has_native_diagnostics()); + advance_radiation(lev, S_old, dt_lev, + &native_shoc_driver[lev]->shoc_cldfrac_diagnostics()); + } + const BoxArray& ba = S_old.boxArray(); const DistributionMapping& dm = S_old.DistributionMap(); diff --git a/Source/TimeIntegration/ERF_AdvanceRadiation.cpp b/Source/TimeIntegration/ERF_AdvanceRadiation.cpp index 46bf4a5bd4..036412a399 100644 --- a/Source/TimeIntegration/ERF_AdvanceRadiation.cpp +++ b/Source/TimeIntegration/ERF_AdvanceRadiation.cpp @@ -4,7 +4,8 @@ using namespace amrex; void ERF::advance_radiation (int lev, MultiFab& cons, - const double& dt_advance) + const double& dt_advance, + const MultiFab* liquid_cloud_fraction) { if (solverChoice.rad_type != RadiationType::None) { #ifdef ERF_USE_NETCDF @@ -40,6 +41,7 @@ void ERF::advance_radiation (int lev, double time_for_rad = t_old[lev] + start_time; rad[lev]->Run(lev, istep[lev], time_for_rad, dt_advance, cons.boxArray(), geom[lev], &(cons), + liquid_cloud_fraction, lmask_lev[lev][0].get(), t_surf, lsm_input_ptrs, lsm_output_ptrs, qheating_rates[lev].get(), rad_fluxes[lev].get(), diff --git a/Tests/Unit/Shoc/CMakeLists.txt b/Tests/Unit/Shoc/CMakeLists.txt index 0a1256b70b..4f6fe1abb8 100644 --- a/Tests/Unit/Shoc/CMakeLists.txt +++ b/Tests/Unit/Shoc/CMakeLists.txt @@ -20,6 +20,7 @@ add_executable(erf_shoc_unit_tests ERF_ShocDriverTests.cpp ERF_ShocPhysicalPropertyTests.cpp ERF_ShocThermoTests.cpp + ERF_ShocRadiationCouplingTests.cpp ) build_erf_exe(erf_shoc_unit_tests NO_ERF_MAIN) @@ -30,6 +31,7 @@ target_include_directories(erf_shoc_unit_tests PRIVATE ${CMAKE_SOURCE_DIR}/Source/PBL ${CMAKE_SOURCE_DIR}/Source/PBL/Shoc ${CMAKE_SOURCE_DIR}/Source/Utils + ${CMAKE_SOURCE_DIR}/Source/PhysicsInterfaces/Radiation ) target_compile_definitions(erf_shoc_unit_tests PRIVATE diff --git a/Tests/Unit/Shoc/ERF_ShocRadiationCouplingTests.cpp b/Tests/Unit/Shoc/ERF_ShocRadiationCouplingTests.cpp new file mode 100644 index 0000000000..59f78ab420 --- /dev/null +++ b/Tests/Unit/Shoc/ERF_ShocRadiationCouplingTests.cpp @@ -0,0 +1,125 @@ +#include "ERF_RadiationCloudFraction.H" + +#include +#include + +namespace { + +using amrex::Real; + +Real +contract_tolerance (const Real expected) +{ + // The cloud-mass helper performs division, a clamp, and two products. + // This is deliberately much smaller than the mixed-phase factor-of-four + // defect that the witness is intended to catch. + const Real scale = (expected > Real(1.0)) ? expected : Real(1.0); + return Real(128.0) * std::numeric_limits::epsilon() * scale; +} + +void +expect_mass_contract (const Real mixing_ratio, + const Real fraction, + const Real expected_grid_mass, + const Real expected_in_cloud_mass, + const Real rho = Real(1.0), + const Real dz = Real(100.0)) +{ + const Real in_cloud_mass = radiation_cloud_mass(mixing_ratio, fraction, rho, dz); + EXPECT_NEAR(in_cloud_mass, expected_in_cloud_mass, + contract_tolerance(expected_in_cloud_mass)); + EXPECT_NEAR(fraction * in_cloud_mass, expected_grid_mass, + contract_tolerance(expected_grid_mass)); +} + +} // namespace + +TEST(ShocRadiationCloudFraction, BinaryClearAndLiquidFallback) +{ + const auto clear = radiation_cloud_fractions(Real(0.0), Real(0.0), Real(0.0), false); + EXPECT_EQ(clear.liquid, Real(0.0)); + EXPECT_EQ(clear.ice, Real(0.0)); + EXPECT_EQ(clear.total, Real(0.0)); + + const auto liquid = radiation_cloud_fractions(Real(1.0e-4), Real(0.0), Real(0.0), false); + EXPECT_EQ(liquid.liquid, Real(1.0)); + EXPECT_EQ(liquid.ice, Real(0.0)); + EXPECT_EQ(liquid.total, Real(1.0)); +} + +TEST(ShocRadiationCloudFraction, SuppliedLiquidFractionIsBounded) +{ + const auto fractional = radiation_cloud_fractions(Real(1.0e-4), Real(0.0), Real(0.25), true); + EXPECT_EQ(fractional.liquid, Real(0.25)); + EXPECT_EQ(fractional.ice, Real(0.0)); + EXPECT_EQ(fractional.total, Real(0.25)); + + const auto below_zero = radiation_cloud_fractions(Real(1.0e-4), Real(0.0), Real(-0.5), true); + EXPECT_EQ(below_zero.liquid, Real(1.0e-4)); + + const auto above_one = radiation_cloud_fractions(Real(1.0e-4), Real(0.0), Real(1.5), true); + EXPECT_EQ(above_one.liquid, Real(1.0)); +} + +TEST(ShocRadiationCloudFraction, IceRemainsBinaryAndPreservesIceOnlyClouds) +{ + const auto ice_only = radiation_cloud_fractions(Real(0.0), Real(1.0e-4), Real(0.0), true); + EXPECT_EQ(ice_only.liquid, Real(0.0)); + EXPECT_EQ(ice_only.ice, Real(1.0)); + EXPECT_EQ(ice_only.total, Real(1.0)); + + const auto mixed = radiation_cloud_fractions(Real(1.0e-4), Real(1.0e-4), Real(0.3), true); + EXPECT_EQ(mixed.liquid, Real(0.3)); + EXPECT_EQ(mixed.ice, Real(1.0)); + EXPECT_EQ(mixed.total, Real(1.0)); + + const auto diagnosed_clear = radiation_cloud_fractions(Real(0.0), Real(0.0), Real(0.7), true); + EXPECT_EQ(diagnosed_clear.liquid, Real(0.0)); + EXPECT_EQ(diagnosed_clear.ice, Real(0.0)); + EXPECT_EQ(diagnosed_clear.total, Real(0.0)); +} + +TEST(ShocRadiationCloudFraction, SharedMaskMassContract) +{ + constexpr Real rho = Real(1.0); + constexpr Real dz = Real(100.0); + + const auto clear = radiation_cloud_fractions(Real(0.0), Real(0.0), Real(0.7), true); + expect_mass_contract(Real(0.0), clear.total, Real(0.0), Real(0.0), rho, dz); + + const auto liquid = radiation_cloud_fractions(Real(1.0e-4), Real(0.0), Real(0.25), true); + expect_mass_contract(Real(1.0e-4), liquid.total, Real(0.01), Real(0.04), rho, dz); + + const auto ice = radiation_cloud_fractions(Real(0.0), Real(1.0e-4), Real(0.0), true); + expect_mass_contract(Real(1.0e-4), ice.total, Real(0.01), Real(0.01), rho, dz); + + const auto mixed = radiation_cloud_fractions(Real(1.0e-4), Real(1.0e-4), Real(0.25), true); + const Real mixed_lwp = radiation_cloud_mass(Real(1.0e-4), mixed.total, rho, dz); + const Real mixed_iwp = radiation_cloud_mass(Real(1.0e-4), mixed.total, rho, dz); + EXPECT_NEAR(mixed.total * mixed_lwp, Real(0.01), contract_tolerance(Real(0.01))); + EXPECT_NEAR(mixed.total * mixed_iwp, Real(0.01), contract_tolerance(Real(0.01))); + + const auto floor_case = radiation_cloud_fractions(Real(1.0e-8), Real(0.0), Real(0.0), true); + EXPECT_EQ(floor_case.total, Real(1.0e-4)); + expect_mass_contract(Real(1.0e-8), floor_case.total, Real(1.0e-6), Real(0.01), rho, dz); + + const auto cap_case = radiation_cloud_fractions(Real(1.0e-2), Real(1.0e-4), Real(0.25), true); + EXPECT_EQ(cap_case.total, Real(1.0)); + expect_mass_contract(Real(1.0e-2), cap_case.total, Real(0.5), Real(0.5), rho, dz); +} + +TEST(ShocRadiationCloudFraction, SamplingFlagSelectsProductionDispatch) +{ + EXPECT_EQ(radiation_cloud_sampling_mode(true), RadiationCloudSamplingMode::MCICA); + EXPECT_EQ(radiation_cloud_sampling_mode(false), RadiationCloudSamplingMode::DeterministicBinary); + + EXPECT_FALSE(radiation_sampling_config_unsupported(true, true, true)); + EXPECT_TRUE(radiation_sampling_config_unsupported(false, true, true)); + EXPECT_FALSE(radiation_sampling_config_unsupported(false, true, false)); + EXPECT_FALSE(radiation_sampling_config_unsupported(false, false, true)); + + // Binary fallback ignores a supplied diagnostic value, as the production + // path does when rad_use_shoc_cldfrac is false. + const auto binary = radiation_cloud_fractions(Real(1.0e-4), Real(0.0), Real(0.25), false); + EXPECT_EQ(binary.total, Real(1.0)); +}