diff --git a/src/mam4xx/ndrop.hpp b/src/mam4xx/ndrop.hpp index 97d307b94..60b97e624 100644 --- a/src/mam4xx/ndrop.hpp +++ b/src/mam4xx/ndrop.hpp @@ -1090,69 +1090,64 @@ void update_from_newcld( } // update_from_newcld KOKKOS_INLINE_FUNCTION -void explmix( - const Real qold_km1, // number / mass mixing ratio from previous time step - // at level k-1 [# or kg / kg] - const Real qold_k, // number / mass mixing ratio from previous time step at - // level k [# or kg / kg] - const Real qold_kp1, // number / mass mixing ratio from previous time step - // at level k+1 [# or kg / kg] - Real & - qnew, // OUTPUT, number / mass mixing ratio to be updated [# or kg / kg] - const Real src, // source due to activation/nucleation at level k [# or kg / - // (kg-s)] - const Real - eddy_diff_kp, // zn*zs*density*diffusivity (kg/m3 m2/s) at interface - // [/s]; below layer k (k,k+1 interface) - const Real - eddy_diff_km, // zn*zs*density*diffusivity (kg/m3 m2/s) at interface - // [/s]; above layer k (k,k+1 interface) - const Real overlapp, // cloud overlap below [fraction] - const Real overlapm, // cloud overlap above [fraction] - const Real dtmix // time step [s] +Real explmix(const Real qold_km1, // number / mass mixing ratio from previous + // time step at level k-1 [# or kg / kg] + const Real qold_k, // number / mass mixing ratio from previous time + // step at level k [# or kg / kg] + const Real qold_kp1, // number / mass mixing ratio from previous + // time step at level k+1 [# or kg / kg] + const Real src, // source due to activation/nucleation at level k + // [# or kg / (kg-s)] + const Real eddy_diff_kp, // zn*zs*density*diffusivity (kg/m3 m2/s) + // at interface + // [/s]; below layer k (k,k+1 interface) + const Real eddy_diff_km, // zn*zs*density*diffusivity (kg/m3 m2/s) + // at interface + // [/s]; above layer k (k,k+1 interface) + const Real overlapp, // cloud overlap below [fraction] + const Real overlapm, // cloud overlap above [fraction] + const Real dtmix // time step [s] ) { - - qnew = qold_k + dtmix * (src + eddy_diff_kp * (overlapp * qold_kp1 - qold_k) + - eddy_diff_km * (overlapm * qold_km1 - qold_k)); - + Real qnew = + qold_k + dtmix * (src + eddy_diff_kp * (overlapp * qold_kp1 - qold_k) + + eddy_diff_km * (overlapm * qold_km1 - qold_k)); // force to non-negative qnew = haero::max(qnew, 0); + // OUTPUT, number / mass mixing ratio to be updated [# or kg / kg] + return qnew; } // end explmix KOKKOS_INLINE_FUNCTION -void explmix( - const Real qold_km1, // number / mass mixing ratio from previous time step - // at level k-1 [# or kg / kg] - const Real qold_k, // number / mass mixing ratio from previous time step at - // level k [# or kg / kg] - const Real qold_kp1, // number / mass mixing ratio from previous time step - // at level k+1 [# or kg / kg] - Real & - qnew, // OUTPUT, number / mass mixing ratio to be updated [# or kg / kg] - const Real src, // source due to activation/nucleation at level k [# or kg / - // (kg-s)] - const Real - eddy_diff_kp, // zn*zs*density*diffusivity (kg/m3 m2/s) at interface - // [/s]; below layer k (k,k+1 interface) - const Real - eddy_diff_km, // zn*zs*density*diffusivity (kg/m3 m2/s) at interface - // [/s]; above layer k (k,k+1 interface) - const Real overlapp, // cloud overlap below [fraction] - const Real overlapm, // cloud overlap above [fraction] - const Real dtmix, // time step [s] - const Real qactold_km1, - // optional: number / mass mixing ratio of ACTIVATED species - // from previous step at level k-1 *** this should only be present if - // the current species is unactivated number/sfc/mass - const Real qactold_kp1 - // optional: number / mass mixing ratio of ACTIVATED species - // from previous step at level k+1 *** this should only be present if - // the current species is unactivated number/sfc/mass +Real explmix(const Real qold_km1, // number / mass mixing ratio from previous + // time step at level k-1 [# or kg / kg] + const Real qold_k, // number / mass mixing ratio from previous time + // step at level k [# or kg / kg] + const Real qold_kp1, // number / mass mixing ratio from previous + // time step at level k+1 [# or kg / kg] + const Real src, // source due to activation/nucleation at level k + // [# or kg / (kg-s)] + const Real eddy_diff_kp, // zn*zs*density*diffusivity (kg/m3 m2/s) + // at interface + // [/s]; below layer k (k,k+1 interface) + const Real eddy_diff_km, // zn*zs*density*diffusivity (kg/m3 m2/s) + // at interface + // [/s]; above layer k (k,k+1 interface) + const Real overlapp, // cloud overlap below [fraction] + const Real overlapm, // cloud overlap above [fraction] + const Real dtmix, // time step [s] + const Real qactold_km1, + // optional: number / mass mixing ratio of ACTIVATED species + // from previous step at level k-1 *** this should only be present + // if the current species is unactivated number/sfc/mass + const Real qactold_kp1 + // optional: number / mass mixing ratio of ACTIVATED species + // from previous step at level k+1 *** this should only be present + // if the current species is unactivated number/sfc/mass ) { // the qactold*(1-overlap) terms are resuspension of activated material const Real one = 1.0; - qnew = + Real qnew = qold_k + dtmix * (-src + @@ -1161,6 +1156,8 @@ void explmix( // force to non-negative qnew = haero::max(qnew, 0); + // OUTPUT, number / mass mixing ratio to be updated [# or kg / kg] + return qnew; } // end explmix KOKKOS_INLINE_FUNCTION @@ -1190,10 +1187,8 @@ void update_from_explmix( const ColumnView &overlapm, // cloud overlap involving level kk-1 [fraction] const ColumnView &eddy_diff_kp, // zn*zs*density*diffusivity [/s] const ColumnView &eddy_diff_km, // zn*zs*density*diffusivity [/s] - const ColumnView &qncld, // updated cloud droplet number mixing ratio [#/kg] - const ColumnView &srcn, // droplet source rate [/s] - // source rate for activated number or species mass [/s] - const ColumnView &source) { + const ColumnView &qncld // updated cloud droplet number mixing ratio [#/kg] +) { // threshold cloud fraction to compute overlap [fraction] // BAD CONSTANT @@ -1201,8 +1196,6 @@ void update_from_explmix( const Real zero = 0.0; const Real one = 1.0; - Real tmpa = zero; // temporary aerosol tendency variable [/s] - constexpr int ntot_amode = AeroConfig::num_modes(); // load new droplets in layers above, below clouds Real dtmin = dtmicro; @@ -1219,7 +1212,7 @@ void update_from_explmix( Kokkos::parallel_reduce( Kokkos::TeamVectorRange(team, top_lev, pver_loc), [&](int k, Real &min_val) { - const int kp1 = haero::min(k + 1, pver - 1); + const int kp1 = haero::min(k + 1, pver_loc - 1); const int km1 = haero::max(k - 1, top_lev); // maximum overlap assumption if (cldn(kp1) > overlap_cld_thresh) { @@ -1283,134 +1276,92 @@ void update_from_explmix( // values of nsav and nnew rather than a physical copying. At end of loop // nnew stores index of most recent updated values (either 1 or 2). + team.team_barrier(); for (int isub = 0; isub < nsubmix; isub++) { - Kokkos::parallel_for(Kokkos::TeamVectorRange(team, top_lev, pver_loc), - [&](int kk) { - qncld(kk) = qcld(kk); - srcn(kk) = zero; - }); - // after first pass, switch nsav, nnew so that nsav is the - // recently updated aerosol - team.team_barrier(); - if (isub > 0) { + if (0 < isub) { + // after first pass, switch nsav, nnew so that nsav is the + // recently updated aerosol const int ntemp = nsav; nsav = nnew; nnew = ntemp; - } // end if - - for (int imode = 0; imode < ntot_amode; imode++) { - const int mm = mam_idx[imode][0] - 1; - - // update droplet source - - // rce-comment- activation source in layer k involves particles from k+1 - // srcn(:)=srcn(:)+nact(:,m)*(raercol(:,mm,nsav)) - Kokkos::parallel_for(Kokkos::TeamVectorRange(team, top_lev, pver_loc), - [&](int k) { - const int kp1 = haero::min(k + 1, pver - 1); - srcn(k) += nact(k, imode) * raercol[kp1][nsav](mm); - }); - - // rce-comment- new formulation for k=pver - // srcn( pver )=srcn( pver )+nact( pver ,m)*(raercol(pver,mm,nsav)) - tmpa = raercol[pver - 1][nsav](mm) * nact(pver - 1, imode) + - raercol_cw[pver - 1][nsav](mm) * nact(pver - 1, imode); - srcn(pver - 1) += haero::max(zero, tmpa); - } // end imode - - // qcld == qold - // qncld == qnew - if (enable_aero_vertical_mix) { - team.team_barrier(); - Kokkos::parallel_for(Kokkos::TeamVectorRange(team, top_lev, pver_loc), - [&](int k) { - const int kp1 = haero::min(k + 1, pver - 1); - const int km1 = haero::max(k - 1, top_lev); - explmix(qncld(km1), qncld(k), qncld(kp1), qcld(k), - srcn(k), eddy_diff_kp(k), eddy_diff_km(k), - overlapp(k), overlapm(k), dtmix); - }); } - - team.team_barrier(); - // update aerosol number - // rce-comment - // the interstitial particle mixratio is different in clear/cloudy - // portions of a layer, and generally higher in the clear portion. (we - // have/had a method for diagnosing the the clear/cloudy mixratios.) the - // activation source terms involve clear air (from below) moving into - // cloudy air (above). in theory, the clear-portion mixratio should be - // used when calculating source terms - for (int imode = 0; imode < ntot_amode; imode++) { - const int mm = mam_idx[imode][0] - 1; - // rce-comment - activation source in layer k involves particles from k+1 - // source(:)= nact(:,m)*(raercol(:,mm,nsav)) - Kokkos::parallel_for( - Kokkos::TeamVectorRange(team, top_lev, pver_loc - 1), [&](int k) { - source(k) = nact(k, imode) * raercol[k + 1][nsav](mm); - }); // end k - - tmpa = raercol[pver - 1][nsav](mm) * nact(pver - 1, imode) + - raercol_cw[pver - 1][nsav](mm) * nact(pver - 1, imode); - source(pver - 1) = haero::max(zero, tmpa); - - Kokkos::parallel_for( - Kokkos::TeamVectorRange(team, top_lev, pver_loc), [&](int k) { - const int kp1 = haero::min(k + 1, pver - 1); - const int km1 = haero::max(k - 1, top_lev); - - explmix(raercol_cw[km1][nsav](mm), raercol_cw[k][nsav](mm), - raercol_cw[kp1][nsav](mm), - raercol_cw[k][nnew](mm), // - source(k), eddy_diff_kp(k), eddy_diff_km(k), overlapp(k), - overlapm(k), dtmix); - if (enable_aero_vertical_mix) { - explmix(raercol[km1][nsav](mm), raercol[k][nsav](mm), - raercol[kp1][nsav](mm), - raercol[k][nnew](mm), // output - source(k), eddy_diff_kp(k), eddy_diff_km(k), overlapp(k), - overlapm(k), dtmix, raercol_cw[km1][nsav](mm), - raercol_cw[kp1][nsav](mm)); // optional in - } - }); // end kk - - // update aerosol species mass - for (int lspec = 1; lspec < nspec_amode[imode] + 1; lspec++) { - const int mm = mam_idx[imode][lspec] - 1; - // rce-comment: activation source in layer k involves particles from k+1 - // source(:)= mact(:,m)*(raercol(:,mm,nsav)) - Kokkos::parallel_for( - Kokkos::TeamVectorRange(team, top_lev, pver_loc), [&](int k) { - const int kp1 = haero::min(k + 1, pver - 1); - source(k) = mact(k, imode) * raercol[kp1][nsav](mm); - }); // end k - tmpa = raercol[pver - 1][nsav](mm) * nact(pver - 1, imode) + - raercol_cw[pver - 1][nsav](mm) * nact(pver - 1, imode); - source(pver - 1) = haero::max(zero, tmpa); - Kokkos::parallel_for( - Kokkos::TeamVectorRange(team, top_lev, pver_loc), [&](int k) { - const int kp1 = haero::min(k + 1, pver - 1); - const int km1 = haero::max(k - 1, top_lev); - explmix(raercol_cw[km1][nsav](mm), raercol_cw[k][nsav](mm), - raercol_cw[kp1][nsav](mm), - raercol_cw[k][nnew](mm), // output - source(k), eddy_diff_kp(k), eddy_diff_km(k), overlapp(k), - overlapm(k), dtmix); + Kokkos::parallel_for(Kokkos::TeamVectorRange(team, top_lev, pver_loc), + [&](int k) { qncld(k) = qcld(k); }); + Kokkos::parallel_for( + Kokkos::TeamVectorRange(team, top_lev, pver_loc), [&](int k) { + const int kp1 = haero::min(k + 1, pver_loc - 1); + const int km1 = haero::max(k - 1, top_lev); + const View1D &raercol_km1_nsav = raercol[km1][nsav]; + const View1D &raercol_k_nsav = raercol[k][nsav]; + const View1D &raercol_kp1_nsav = raercol[kp1][nsav]; + const View1D &raercol_k_nnew = raercol[k][nnew]; + const View1D &raercol_cw_km1_nsav = raercol_cw[km1][nsav]; + const View1D &raercol_cw_k_nsav = raercol_cw[k][nsav]; + const View1D &raercol_cw_kp1_nsav = raercol_cw[kp1][nsav]; + const View1D &raercol_cw_k_nnew = raercol_cw[k][nnew]; + // update droplet source + // rce-comment- activation source in layer k involves particles from + // k+1 + // srcn(:)=srcn(:)+nact(:,m)*(raercol(:,mm,nsav)) + Real srcn = zero; + if (k < pver_loc - 1) { + for (int imode = 0; imode < ntot_amode; imode++) { + const int mm = mam_idx[imode][0] - 1; + srcn += nact(k, imode) * raercol_kp1_nsav(mm); + } // end imode + } else { + for (int imode = 0; imode < ntot_amode; imode++) { + const int mm = mam_idx[imode][0] - 1; + const Real tmpa = raercol_k_nsav(mm) * nact(k, imode) + + raercol_cw_k_nsav(mm) * nact(k, imode); + srcn += haero::max(zero, tmpa); + } // end imode + } + // update aerosol number + // rce-comment + // the interstitial particle mixratio is different in clear/cloudy + // portions of a layer, and generally higher in the clear portion. + // (we have/had a method for diagnosing the the clear/cloudy + // mixratios.) the activation source terms involve clear air (from + // below) moving into cloudy air (above). in theory, the + // clear-portion mixratio should be used when calculating source + // terms + // rce-comment: activation source in layer k involves particles from + // k+1 source(:)= mact(:,m)*(raercol(:,mm,nsav)) + if (enable_aero_vertical_mix) { + qcld(k) = + explmix(qncld(km1), qncld(k), qncld(kp1), srcn, eddy_diff_kp(k), + eddy_diff_km(k), overlapp(k), overlapm(k), dtmix); + } + for (int imode = 0; imode < ntot_amode; imode++) { + for (int lspec = 0; lspec < nspec_amode[imode] + 1; lspec++) { + const int mm = mam_idx[imode][lspec] - 1; + Real source = 0; + if (k < pver_loc - 1) { + const Real act = lspec ? mact(k, imode) : nact(k, imode); + source = act * raercol_kp1_nsav(mm); + } else { + const Real tmpa = raercol_k_nsav(mm) * nact(k, imode) + + raercol_cw_k_nsav(mm) * nact(k, imode); + source = haero::max(zero, tmpa); + } + // update aerosol species mass + raercol_cw_k_nnew(mm) = + explmix(raercol_cw_km1_nsav(mm), raercol_cw_k_nsav(mm), + raercol_cw_kp1_nsav(mm), source, eddy_diff_kp(k), + eddy_diff_km(k), overlapp(k), overlapm(k), dtmix); if (enable_aero_vertical_mix) { - explmix(raercol[km1][nsav](mm), raercol[k][nsav](mm), - raercol[kp1][nsav](mm), - raercol[k][nnew](mm), // output - source(k), eddy_diff_kp(k), eddy_diff_km(k), - overlapp(k), overlapm(k), dtmix, - raercol_cw[km1][nsav](mm), - raercol_cw[kp1][nsav](mm)); // optional in + raercol_k_nnew(mm) = + explmix(raercol_km1_nsav(mm), raercol_k_nsav(mm), + raercol_kp1_nsav(mm), source, eddy_diff_kp(k), + eddy_diff_km(k), overlapp(k), overlapm(k), dtmix, + raercol_cw_km1_nsav(mm), raercol_cw_kp1_nsav(mm)); } - }); // end kk - - team.team_barrier(); - } // lspec loop - } // imode loop - } // old_cloud_nsubmix_loop + } // lspec loop + } // imode loop + }); // k loop + team.team_barrier(); + } // old_cloud_nsubmix_loop // evaporate particles again if no cloud Kokkos::parallel_for( @@ -1585,7 +1536,6 @@ void dropmixnuc( csbot_cscen(surface_cell) = one; // Initialize 1D (in space) versions of interstitial and cloud borne aerosol int nsav = 0; - Kokkos::parallel_for( Kokkos::TeamVectorRange(team, top_lev, pver_loc), [&](int k) { for (int imode = 0; imode < ntot_amode; ++imode) { @@ -1676,9 +1626,7 @@ void dropmixnuc( qcld, raercol, raercol_cw, nsav, nnew, nspec_amode, mam_idx, enable_aero_vertical_mix, top_lev, // work vars - overlapp, overlapm, eddy_diff_kp, eddy_diff_km, qncld, - srcn, // droplet source rate [/s] - source); + overlapp, overlapm, eddy_diff_kp, eddy_diff_km, qncld); team.team_barrier(); diff --git a/src/tests/mam4_ndrop_unit_tests.cpp b/src/tests/mam4_ndrop_unit_tests.cpp index 705d86e89..e15aa9665 100644 --- a/src/tests/mam4_ndrop_unit_tests.cpp +++ b/src/tests/mam4_ndrop_unit_tests.cpp @@ -206,14 +206,14 @@ TEST_CASE("test_explmix", "mam4_ndrop") { Real qactold_km1 = 1; Real qactold_kp1 = 1; - ndrop::explmix(qold_km1, qold_k, qold_kp1, q, src, ek_kp1, ek_km1, - overlap_kp1, overlap_km1, dt); + q = ndrop::explmix(qold_km1, qold_k, qold_kp1, src, ek_kp1, ek_km1, + overlap_kp1, overlap_km1, dt); logger.info("q = {}", q); REQUIRE(FloatingPoint::equiv(q, 1.1)); - ndrop::explmix(qold_km1, qold_k, qold_kp1, q, src, ek_kp1, ek_km1, - overlap_kp1, overlap_km1, dt, qactold_km1, qactold_kp1); + q = ndrop::explmix(qold_km1, qold_k, qold_kp1, src, ek_kp1, ek_km1, + overlap_kp1, overlap_km1, dt, qactold_km1, qactold_kp1); logger.info("q = {}", q); REQUIRE(FloatingPoint::equiv(q, 0.9)); diff --git a/src/validation/ndrop/explmix.cpp b/src/validation/ndrop/explmix.cpp index 17ec24172..89c44e333 100644 --- a/src/validation/ndrop/explmix.cpp +++ b/src/validation/ndrop/explmix.cpp @@ -57,11 +57,12 @@ void explmix(Ensemble *ensemble) { if (is_unact) { qactold_km1 = qactold[km1]; qactold_kp1 = qactold[kp1]; - ndrop::explmix(qold_km1, qold_k, qold_kp1, q[k], src, ekkp, ekkm, - overlapp, overlapm, dtmix, qactold_km1, qactold_kp1); + q[k] = + ndrop::explmix(qold_km1, qold_k, qold_kp1, src, ekkp, ekkm, + overlapp, overlapm, dtmix, qactold_km1, qactold_kp1); } else { - ndrop::explmix(qold_km1, qold_k, qold_kp1, q[k], src, ekkp, ekkm, - overlapp, overlapm, dtmix); + q[k] = ndrop::explmix(qold_km1, qold_k, qold_kp1, src, ekkp, ekkm, + overlapp, overlapm, dtmix); } } diff --git a/src/validation/ndrop/update_from_explmix.cpp b/src/validation/ndrop/update_from_explmix.cpp index 1c8da3913..aa8d57e4b 100644 --- a/src/validation/ndrop/update_from_explmix.cpp +++ b/src/validation/ndrop/update_from_explmix.cpp @@ -173,7 +173,7 @@ void update_from_explmix(Ensemble *ensemble) { ndrop::update_from_explmix( team, dtmicro, csbot, cldn, zn, zs, ekd, nact, mact, qcld, raercol, raercol_cw, nsav, nnew, nspec_amode, mam_idx, true, - top_lev, overlapp, overlapm, ekkp, ekkm, qncld, srcn, source); + top_lev, overlapp, overlapm, ekkp, ekkm, qncld); indexes(0) = nnew; indexes(1) = nsav; });