ERF_InitCustomPert_SuperCell.H

Go to the documentation of this file.

 
    // On a mesh other than ConstantDz, heights come from z_cc and the vertical levels,
    //    not from geomdata's prob_lo[2] and CellSize(2)
    const bool use_z_phys = use_physical_z();
 
    const bool use_moisture = (sc.moisture_type != MoistureType::None);
 
    const int khi = geomdata.Domain().bigEnd()[2];
 
    AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1);
 
    ParmParse pp_prob("prob");
    Real q_t = amrex::Real(0.02);
    pp_prob.query("qt_init",q_t);
 
    Real eq_pot_temp = amrex::Real(320.0);
    pp_prob.query("eq_pot_temp",eq_pot_temp);
 
    bool use_empirical = false;
    pp_prob.query("use_empirical_psat",use_empirical);
 
    // This is what we do at k = 0 -- note we assume p = p_0 and T = T_0 at z=0
    const amrex::Real& dz        = geomdata.CellSize()[2];
 
    // Call the routine to calculate the 1d background condition
 
    Vector<Real> h_r(khi+2);
    Vector<Real> h_p(khi+2);
    Vector<Real> h_t(khi+2);
    Vector<Real> h_q_v(khi+2);
 
    amrex::Gpu::DeviceVector<Real> d_r(khi+2);
    amrex::Gpu::DeviceVector<Real> d_p(khi+2);
    amrex::Gpu::DeviceVector<Real> d_t(khi+2);
    amrex::Gpu::DeviceVector<Real> d_q_v(khi+2);
 
    Real height =  amrex::Real(1200.);
    pp_prob.query("height",height);
 
    Real   z_tr = amrex::Real(12000.);
    pp_prob.query("z_tr",z_tr);
 
    Real x_c  = zero  ; pp_prob.query("x_c", x_c);
    Real y_c  = zero  ; pp_prob.query("y_c", y_c);
    Real z_c  = amrex::Real(1.5e3); pp_prob.query("z_c", z_c);
 
    Real r_c = one;
 
    Real x_r  = amrex::Real(10.0e3); pp_prob.query("x_r", x_r);
    Real y_r  = amrex::Real(10.0e3); pp_prob.query("y_r", y_r);
    Real z_r  = amrex::Real(1.5e3); pp_prob.query("z_r", z_r);
 
    Real T_tr     = amrex::Real(213.0); pp_prob.query("T_tr", T_tr);
    Real theta_tr = amrex::Real(343.0); pp_prob.query("theta_tr",theta_tr);
    Real theta_c  =   three; pp_prob.query("theta_c", theta_c);
    Real theta_0  = amrex::Real(300.0); pp_prob.query("theta_0",theta_0);
 
    if (use_z_phys) {
        HSEutils::init_isentropic_hse_no_terrain(h_t.data(), h_r.data(),h_p.data(),h_q_v.data(),m_zlevels_stag[lev].data(),khi,
                                                 q_t,eq_pot_temp,use_empirical,true,height,z_tr,
                                                 theta_0, theta_tr, T_tr);
    } else {
        HSEutils::init_isentropic_hse_no_terrain(h_t.data(), h_r.data(),h_p.data(),h_q_v.data(),dz,geomdata.ProbLo(2),khi,
                                                 q_t,eq_pot_temp,use_empirical,true,height,z_tr,
                                                 theta_0, theta_tr, T_tr);
    }
 
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_r.begin(), h_r.end(), d_r.begin());
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_p.begin(), h_p.end(), d_p.begin());
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_t.begin(), h_t.end(), d_t.begin());
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_q_v.begin(), h_q_v.end(), d_q_v.begin());
 
    Real* t   = d_t.data();
    Real* p   = d_p.data();
 
    ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
    {
        // Geometry (note we must include these here to get the data on device)
        const auto prob_lo  = geomdata.ProbLo();
        const auto dx       = geomdata.CellSize();
        const Real x        = prob_lo[0] + (i + myhalf) * dx[0];
        const Real y        = prob_lo[1] + (j + myhalf) * dx[1];
        const Real z        = use_z_phys ? z_cc(i,j,k) : prob_lo[2] + (k + myhalf) * dx[2];
        // The humidity zones (height, z_tr) are heights above the ground, as in the profile
        //    from init_isentropic_hse_no_terrain; the bubble centre z_c is a coordinate, like x_c
        const Real z_agl    = z - erf_prob_heights::z_ground_cc(i,j,z_sfc_nd);
        Real rad, delta_theta, theta_total, temperature, rho, RH, scalar;
 
        // Introduce the warm bubble. Assume that the bubble is pressure matched with the background
        rad = std::sqrt( amrex::Math::powi<2>((x - x_c)/x_r) +
                         amrex::Math::powi<2>((y - y_c)/y_r) +
                         amrex::Math::powi<2>((z - z_c)/z_r) );
 
        if(rad <= r_c) {
            delta_theta = theta_c * amrex::Math::powi<2>(std::cos(PI*rad/two));
            scalar = amrex::Math::powi<2>(std::cos(PI*rad/two));
        } else {
            delta_theta = zero;
            scalar = zero;
        }
 
        Real scaled_height = z_agl / z_tr;
        int which_zone = -1;
        if (z_agl <= height) {
            which_zone = 1;
        } else if (z_agl <= z_tr) {
            which_zone = 2;
        } else {
            which_zone = 3;
        }
 
        theta_total     = t[k] + delta_theta;
        temperature     = getTgivenPandTh(p[k], theta_total, RdoCp);
        Real T_b        = getTgivenPandTh(p[k], t[k]       , RdoCp);
        RH              = HSEutils::compute_relative_humidity(p[k], T_b, use_empirical, which_zone, scaled_height);
 
        Real q_v_hot;
        if (z_agl < height) {
            q_v_hot = amrex::Real(0.014);
        } else {
            q_v_hot = HSEutils::vapor_mixing_ratio(p[k], T_b, RH, use_empirical, which_zone);
        }
 
        rho = p[k]/(R_d*temperature*(one + RvoRd*q_v_hot));
 
        // Compute background quantities
        Real temperature_back = getTgivenPandTh(p[k], t[k], RdoCp);
        Real T_back           = getTgivenPandTh(p[k], t[k], RdoCp);
        Real RH_back          = HSEutils::compute_relative_humidity(p[k], T_back, use_empirical, which_zone, scaled_height);
 
        Real q_v_back;
        if (z_agl < height) {
            q_v_back = amrex::Real(0.014);
        } else {
            q_v_back = HSEutils::vapor_mixing_ratio(p[k], T_back, RH_back, use_empirical, which_zone);
        }
 
        Real rho_back         = getRhogivenTandPress(temperature_back, p[k], q_v_back);
 
        // This version perturbs rho but not p
        state_pert(i, j, k, RhoTheta_comp) = rho*theta_total - rho_back*t[k]*(one + RvoRd*q_v_back);// rho*d_t[k]*(one + R_v_by_R_d*q_v_hot);
        state_pert(i, j, k, Rho_comp)      = rho - rho_back*(one + q_v_back);
 
        // Set scalar = 0 everywhere
        state_pert(i, j, k, RhoScalar_comp) = rho*scalar;
 
        // mean states
        if (use_moisture) {
            state_pert(i, j, k, RhoQ1_comp) = rho*q_v_hot;
            state_pert(i, j, k, RhoQ2_comp) = zero;
        }
 
    });
 
    // The kernel above reads d_r, d_p, d_t and d_q_v; wait for it before they are freed
    Gpu::streamSynchronize();