From ebe20800b5d900dd1803db10293438bc9d6964a9 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Thu, 4 Apr 2024 16:00:00 -0600 Subject: [PATCH 1/9] fix some PBL conditionals - QKE only for MYNN, not YSU --- Source/DataStructs/TurbStruct.H | 54 +++++++++++---------- Source/Diffusion/DiffusionSrcForState_N.cpp | 3 +- Source/Initialization/ERF_init_custom.cpp | 4 +- 3 files changed, 33 insertions(+), 28 deletions(-) diff --git a/Source/DataStructs/TurbStruct.H b/Source/DataStructs/TurbStruct.H index 31a16e391..1b5078b64 100644 --- a/Source/DataStructs/TurbStruct.H +++ b/Source/DataStructs/TurbStruct.H @@ -78,20 +78,22 @@ struct TurbChoice { amrex::Print() << "Selected a PBL model and an LES model: " << "Using PBL for vertical transport, LES for horizontal" << std::endl; } - pp.query("pbl_mynn_A1", pbl_mynn_A1); - pp.query("pbl_mynn_A2", pbl_mynn_A2); - pp.query("pbl_mynn_B1", pbl_mynn_B1); - pp.query("pbl_mynn_B2", pbl_mynn_B2); - pp.query("pbl_mynn_C1", pbl_mynn_C1); - pp.query("pbl_mynn_C2", pbl_mynn_C2); - pp.query("pbl_mynn_C3", pbl_mynn_C3); - pp.query("pbl_mynn_C4", pbl_mynn_C4); - pp.query("pbl_mynn_C5", pbl_mynn_C5); + + if (pbl_type == PBLType::MYNN25) { + pp.query("pbl_mynn_A1", pbl_mynn_A1); + pp.query("pbl_mynn_A2", pbl_mynn_A2); + pp.query("pbl_mynn_B1", pbl_mynn_B1); + pp.query("pbl_mynn_B2", pbl_mynn_B2); + pp.query("pbl_mynn_C1", pbl_mynn_C1); + pp.query("pbl_mynn_C2", pbl_mynn_C2); + pp.query("pbl_mynn_C3", pbl_mynn_C3); + pp.query("pbl_mynn_C4", pbl_mynn_C4); + pp.query("pbl_mynn_C5", pbl_mynn_C5); + } } // Right now, solving the QKE equation is only supported when MYNN PBL is turned on - if (pbl_type == PBLType::MYNN25 || - pbl_type == PBLType::YSU ) { + if (pbl_type == PBLType::MYNN25) { use_QKE = true; pp.query("diffuse_QKE_3D", diffuse_QKE_3D); pp.query("advect_QKE", advect_QKE); @@ -163,20 +165,22 @@ struct TurbChoice { } else if (les_type == LESType::Deardorff) { amrex::Error("It is not recommended to use Deardorff LES and a PBL model"); } - pp.query("pbl_mynn_A1", pbl_mynn_A1, lev); - pp.query("pbl_mynn_A2", pbl_mynn_A2, lev); - pp.query("pbl_mynn_B1", pbl_mynn_B1, lev); - pp.query("pbl_mynn_B2", pbl_mynn_B2, lev); - pp.query("pbl_mynn_C1", pbl_mynn_C1, lev); - pp.query("pbl_mynn_C2", pbl_mynn_C2, lev); - pp.query("pbl_mynn_C3", pbl_mynn_C3, lev); - pp.query("pbl_mynn_C4", pbl_mynn_C4, lev); - pp.query("pbl_mynn_C5", pbl_mynn_C5, lev); + + if (pbl_type == PBLType::MYNN25) { + pp.query("pbl_mynn_A1", pbl_mynn_A1, lev); + pp.query("pbl_mynn_A2", pbl_mynn_A2, lev); + pp.query("pbl_mynn_B1", pbl_mynn_B1, lev); + pp.query("pbl_mynn_B2", pbl_mynn_B2, lev); + pp.query("pbl_mynn_C1", pbl_mynn_C1, lev); + pp.query("pbl_mynn_C2", pbl_mynn_C2, lev); + pp.query("pbl_mynn_C3", pbl_mynn_C3, lev); + pp.query("pbl_mynn_C4", pbl_mynn_C4, lev); + pp.query("pbl_mynn_C5", pbl_mynn_C5, lev); + } } // Right now, solving the QKE equation is only supported when MYNN PBL is turned on - if (pbl_type == PBLType::MYNN25 || - pbl_type == PBLType::YSU ) { + if (pbl_type == PBLType::MYNN25) { use_QKE = true; pp.query("diffuse_QKE_3D", diffuse_QKE_3D, lev); pp.query("advect_QKE" , advect_QKE, lev); @@ -236,7 +240,7 @@ struct TurbChoice { amrex::Print() << "reference theta : " << theta_ref << std::endl; } - if (pbl_type != PBLType::None) { + if (pbl_type == PBLType::MYNN25) { amrex::Print() << "pbl_mynn_A1 : " << pbl_mynn_A1 << std::endl; amrex::Print() << "pbl_mynn_A2 : " << pbl_mynn_A2 << std::endl; amrex::Print() << "pbl_mynn_B1 : " << pbl_mynn_B1 << std::endl; @@ -293,9 +297,9 @@ struct TurbChoice { amrex::Real pbl_ysu_coriolis_freq = 1.0e-4; // TODO: make this consistent with coriolis forcing (but note WRF always uses 1e-4) bool pbl_ysu_over_land = false; // TODO: pull from other inputs and make local amrex::Real pbl_ysu_land_Ribcr = 0.25; // Critical Bulk Richardson number of Land for stable conditions - amrex::Real pbl_ysu_unst_Ribcr = 0.0; + amrex::Real pbl_ysu_unst_Ribcr = 0.0; // Critical Bulk Richardson number for unstable conditions - // QKE stuff - default is not to use it, if MYNN2.5 or YSU PBL is used default is turb transport in Z-direction only + // QKE stuff - default is not to use it, if MYNN2.5 PBL is used default is turb transport in Z-direction only bool use_QKE = false; bool diffuse_QKE_3D = false; bool advect_QKE = true; diff --git a/Source/Diffusion/DiffusionSrcForState_N.cpp b/Source/Diffusion/DiffusionSrcForState_N.cpp index 5e85aa285..33cb3e136 100644 --- a/Source/Diffusion/DiffusionSrcForState_N.cpp +++ b/Source/Diffusion/DiffusionSrcForState_N.cpp @@ -78,7 +78,8 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, bool l_consA = (diffChoice.molec_diff_type == MolecDiffType::ConstantAlpha); bool l_turb = ( (turbChoice.les_type == LESType::Smagorinsky) || (turbChoice.les_type == LESType::Deardorff ) || - (turbChoice.pbl_type == PBLType::MYNN25 ) ); + (turbChoice.pbl_type == PBLType::MYNN25 ) || + (turbChoice.pbl_type == PBLType::YSU ) ); const Box xbx = surroundingNodes(bx,0); const Box ybx = surroundingNodes(bx,1); diff --git a/Source/Initialization/ERF_init_custom.cpp b/Source/Initialization/ERF_init_custom.cpp index 98e935d9f..ca73ac1f8 100644 --- a/Source/Initialization/ERF_init_custom.cpp +++ b/Source/Initialization/ERF_init_custom.cpp @@ -95,8 +95,8 @@ ERF::init_custom (int lev) MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoKE_comp, RhoKE_comp, 1, cons_pert.nGrow()); } - // RhoQKE is only relevant if using MYNN2.5 or YSU - if (solverChoice.turbChoice[lev].pbl_type == PBLType::None) { + // RhoQKE is only relevant if using MYNN2.5 + if (solverChoice.turbChoice[lev].pbl_type != PBLType::MYNN25) { lev_new[Vars::cons].setVal(0.0,RhoQKE_comp,1); } else { MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQKE_comp, RhoQKE_comp, 1, cons_pert.nGrow()); From afabd95c16dbc7a38ba91fa503b82ffd1ca43674 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Thu, 4 Apr 2024 16:53:06 -0600 Subject: [PATCH 2/9] towards GABLS test case for YSU --- ...{inputs_GABLS1 => inputs_GABLS1_deardorff} | 0 Exec/ABL/inputs_GABLS1_mynn25 | 99 +++++++++++++++++++ Exec/ABL/inputs_GABLS1_ysu | 99 +++++++++++++++++++ 3 files changed, 198 insertions(+) rename Exec/ABL/{inputs_GABLS1 => inputs_GABLS1_deardorff} (100%) create mode 100644 Exec/ABL/inputs_GABLS1_mynn25 create mode 100644 Exec/ABL/inputs_GABLS1_ysu diff --git a/Exec/ABL/inputs_GABLS1 b/Exec/ABL/inputs_GABLS1_deardorff similarity index 100% rename from Exec/ABL/inputs_GABLS1 rename to Exec/ABL/inputs_GABLS1_deardorff diff --git a/Exec/ABL/inputs_GABLS1_mynn25 b/Exec/ABL/inputs_GABLS1_mynn25 new file mode 100644 index 000000000..6806e41d7 --- /dev/null +++ b/Exec/ABL/inputs_GABLS1_mynn25 @@ -0,0 +1,99 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +stop_time = 32400.0 # 540 min = 9 h +stop_time = 600.0 # shorter tests + +amrex.fpe_trap_invalid = 0 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY (dx=dy=dz=3.125 m) +geometry.prob_extent = 25 25 400 +amr.n_cell = 4 4 64 +geometry.is_periodic = 1 1 0 + +# MOST BOUNDARY (DEFAULT IS ADIABATIC FOR THETA) +zlo.type = "Most" +erf.most.z0 = 0.1 # from Beare et al. 2006 +# erf.most.zref = 3.125 # == dz/2 +erf.most.surf_temp = 265.0 # initial value, should match input_sounding +erf.most.surf_heating_rate = -0.25 # [K/h] from Beare et al. 2006 + +zhi.type = "SlipWall" +zhi.theta_grad = 0.01 # [K/m] to match the input sounding + +# INITIALIZATION (Beare et al. 2006) +erf.init_type = "input_sounding" +erf.init_sounding_ideal = 1 +erf.input_sounding_file = "input_sounding_GABLS1" + +# TIME STEP CONTROL +erf.fixed_dt = 1.0 +erf.fixed_mri_dt_ratio = 6 + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# planar averaging (with erf.v = 1) +erf.data_log = surf_hist.dat mean_profiles.dat covar_profiles.dat sfs_profiles.dat +erf.profile_int = 20 +erf.interp_profiles_to_cc = false + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = 6000 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 6000 # number of timesteps between plotfiles +erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta +erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta rhoQKE Kmv Khv z_phys + +# SOLVER CHOICE +erf.dycore_horiz_adv_type = "Upwind_3rd" +erf.dycore_vert_adv_type = "Upwind_3rd" +erf.dryscal_horiz_adv_type = "Upwind_3rd" +erf.dryscal_vert_adv_type = "Upwind_3rd" + +erf.molec_diff_type = "None" + +erf.use_gravity = true + +# Coriolis parameter f = 1.39e-4 s^-1, from Beare et al. 2006 +erf.use_coriolis = true +erf.latitude = 73.0 +erf.rotational_time_period = 86455.2516813368 + +# Geostrophic wind (Beare et al. 2006) +erf.abl_driver_type = "GeostrophicWind" +erf.abl_geo_wind = 8.0 0.0 0.0 + +# Turbulence closure +erf.les_type = "Smagorinsky" +erf.Cs = 0.17 # following Gadde et al 2020 BLM +#erf.les_type = "Deardorff" +#erf.Ck = 0.1 # Coefficient in Moeng1984, Eqn. 19 +#erf.Ce = 0.7 # Note: Ce_lcoeff = C_e - 1.9*C_k +#erf.Ce_wall = 3.9 # To account for "wall effects" +#erf.sigma_k = 0.5 # TKE diffusion coeff = 2*Km = Km * inv_sigma_k, Moeng1984, Eqn. 15 + +# NOT USED +#erf.rho0_trans = 1.3223 # from Beare et al. 2006 +#erf.theta_ref = 263.5 # from Beare et al. 2006 + +#erf.use_rayleigh_damping = true +#prob.dampcoef = 0.2 # [1/s] following FastEddy +#prob.zdamp = 100. # from Beare et al. 2006, "most models applied gravity wave damping above 300 m" + +# Initial conditions from Beare et al. 2006 +#prob.KE_0 = 0.4 # [m2/s2] +#prob.KE_decay_height = 250. # [m] +#prob.KE_decay_order = 1 +#prob.T_0_Pert_Mag = 0.1 # [K] +#prob.pert_ref_height = 50. # [m] + +erf.pbl_type = "MYNN2.5" \ No newline at end of file diff --git a/Exec/ABL/inputs_GABLS1_ysu b/Exec/ABL/inputs_GABLS1_ysu new file mode 100644 index 000000000..cdce595f8 --- /dev/null +++ b/Exec/ABL/inputs_GABLS1_ysu @@ -0,0 +1,99 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +stop_time = 32400.0 # 540 min = 9 h +stop_time = 600.0 # shorter tests + +amrex.fpe_trap_invalid = 0 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY (dx=dy=dz=3.125 m) +geometry.prob_extent = 25 25 400 +amr.n_cell = 4 4 64 +geometry.is_periodic = 1 1 0 + +# MOST BOUNDARY (DEFAULT IS ADIABATIC FOR THETA) +zlo.type = "Most" +erf.most.z0 = 0.1 # from Beare et al. 2006 +# erf.most.zref = 3.125 # == dz/2 +erf.most.surf_temp = 265.0 # initial value, should match input_sounding +erf.most.surf_heating_rate = -0.25 # [K/h] from Beare et al. 2006 + +zhi.type = "SlipWall" +zhi.theta_grad = 0.01 # [K/m] to match the input sounding + +# INITIALIZATION (Beare et al. 2006) +erf.init_type = "input_sounding" +erf.init_sounding_ideal = 1 +erf.input_sounding_file = "input_sounding_GABLS1" + +# TIME STEP CONTROL +erf.fixed_dt = 1.0 +erf.fixed_mri_dt_ratio = 6 + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# planar averaging (with erf.v = 1) +erf.data_log = surf_hist.dat mean_profiles.dat covar_profiles.dat sfs_profiles.dat +erf.profile_int = 20 +erf.interp_profiles_to_cc = false + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = 6000 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 6000 # number of timesteps between plotfiles +erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta +erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta rhoQKE Kmv Khv z_phys + +# SOLVER CHOICE +erf.dycore_horiz_adv_type = "Upwind_3rd" +erf.dycore_vert_adv_type = "Upwind_3rd" +erf.dryscal_horiz_adv_type = "Upwind_3rd" +erf.dryscal_vert_adv_type = "Upwind_3rd" + +erf.molec_diff_type = "None" + +erf.use_gravity = true + +# Coriolis parameter f = 1.39e-4 s^-1, from Beare et al. 2006 +erf.use_coriolis = true +erf.latitude = 73.0 +erf.rotational_time_period = 86455.2516813368 + +# Geostrophic wind (Beare et al. 2006) +erf.abl_driver_type = "GeostrophicWind" +erf.abl_geo_wind = 8.0 0.0 0.0 + +# Turbulence closure +erf.les_type = "Smagorinsky" +erf.Cs = 0.17 # following Gadde et al 2020 BLM +#erf.les_type = "Deardorff" +#erf.Ck = 0.1 # Coefficient in Moeng1984, Eqn. 19 +#erf.Ce = 0.7 # Note: Ce_lcoeff = C_e - 1.9*C_k +#erf.Ce_wall = 3.9 # To account for "wall effects" +#erf.sigma_k = 0.5 # TKE diffusion coeff = 2*Km = Km * inv_sigma_k, Moeng1984, Eqn. 15 + +# NOT USED +#erf.rho0_trans = 1.3223 # from Beare et al. 2006 +#erf.theta_ref = 263.5 # from Beare et al. 2006 + +#erf.use_rayleigh_damping = true +#prob.dampcoef = 0.2 # [1/s] following FastEddy +#prob.zdamp = 100. # from Beare et al. 2006, "most models applied gravity wave damping above 300 m" + +# Initial conditions from Beare et al. 2006 +#prob.KE_0 = 0.4 # [m2/s2] +#prob.KE_decay_height = 250. # [m] +#prob.KE_decay_order = 1 +#prob.T_0_Pert_Mag = 0.1 # [K] +#prob.pert_ref_height = 50. # [m] + +erf.pbl_type = "YSU" \ No newline at end of file From fbb32bafd39486f3496239e24677cdadf4e32363 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Thu, 4 Apr 2024 16:53:30 -0600 Subject: [PATCH 3/9] remove a rogue amrex:: --- Source/Diffusion/PBLModels.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp index 75f197ad2..b74d48065 100644 --- a/Source/Diffusion/PBLModels.cpp +++ b/Source/Diffusion/PBLModels.cpp @@ -285,7 +285,7 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel, // Require that MOST zref is 10 m so we get the wind speed at 10 m from most if (most_zref != 10.0) { - amrex::Abort("MOST Zref must be 10m for YSU PBL scheme"); + Abort("MOST Zref must be 10m for YSU PBL scheme"); } // create flattened boxes to store PBL height From 321e46a839cb7899261483d8f3982874e3ef6129 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Thu, 4 Apr 2024 17:31:25 -0600 Subject: [PATCH 4/9] WIP on testing YSU progress + fix for no 'terrain --- Source/Diffusion/PBLModels.cpp | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp index b74d48065..85df90d25 100644 --- a/Source/Diffusion/PBLModels.cpp +++ b/Source/Diffusion/PBLModels.cpp @@ -245,7 +245,8 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel, }); } } else if (turbChoice.pbl_type == PBLType::YSU) { - Error("YSU Model not implemented yet"); + // FIXME + // Error("YSU Model not implemented yet"); /* YSU PBL initially introduced by S.-Y. Hong, Y. Noh, and J. Dudhia, MWR, 2006 [HND06] @@ -280,11 +281,13 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel, //const auto& t_star_arr = most->get_t_star(level)->const_array(mfi); const auto& t_surf_arr = most->get_t_surf(level)->const_array(mfi); //const auto& l_obuk_arr = most->get_olen(level)->const_array(mfi); - const auto& z_nd_arr = z_phys_nd->array(mfi); + const bool use_terrain = (z_phys_nd != nullptr); + const Array4 z_nd_arr = use_terrain ? z_phys_nd->array(mfi) : Array4{}; const Real most_zref = most->get_zref(); // Require that MOST zref is 10 m so we get the wind speed at 10 m from most if (most_zref != 10.0) { + amrex::Print() << "most_zref = " << most_zref << std::endl; Abort("MOST Zref must be 10m for YSU PBL scheme"); } @@ -357,6 +360,9 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel, if (pblh_arr(i,j,0) < 0.5*(zval_up+zval_dn) ) { kpbl = 0; } + + // FIXME: DEBUG + std::cout << i << " " << j << " " << pblh_arr(i,j,0) << std::endl; }); // -- Compute nonlocal/countergradient mixing parameters @@ -367,6 +373,8 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel, // -- Compute coefficients in free stream above PBL + // FIXME: DEBUG + Error("YSU Model not implemented yet"); } } From 4bda1eba2eea3b37bbbf8a9e9aa6e4c241804c8c Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Thu, 4 Apr 2024 18:04:15 -0600 Subject: [PATCH 5/9] fix for most zref check --- Source/Diffusion/PBLModels.cpp | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp index 85df90d25..46b44b223 100644 --- a/Source/Diffusion/PBLModels.cpp +++ b/Source/Diffusion/PBLModels.cpp @@ -286,8 +286,16 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel, const Real most_zref = most->get_zref(); // Require that MOST zref is 10 m so we get the wind speed at 10 m from most - if (most_zref != 10.0) { - amrex::Print() << "most_zref = " << most_zref << std::endl; + bool invalid_zref = false; + if (use_terrain) { + invalid_zref = most_zref != 10.0; + } else { + // zref gets reset to nearest cell center, so assert that zref is in the same cell as the 10m point + Real dz = geom.CellSize(2); + invalid_zref = int((most_zref - 0.5*dz)/dz) != int((10.0 - 0.5*dz)/dz); + } + if (invalid_zref) { + Print() << "most_zref = " << most_zref << std::endl; Abort("MOST Zref must be 10m for YSU PBL scheme"); } From f169c73e6db62aa3ee1d644530bfc175d1e27813 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Wed, 10 Apr 2024 17:56:08 -0600 Subject: [PATCH 6/9] update pbl test files --- Exec/ABL/inputs_GABLS1_mynn25 | 70 +++++++++++++---------------------- Exec/ABL/inputs_GABLS1_ysu | 70 +++++++++++++---------------------- 2 files changed, 50 insertions(+), 90 deletions(-) diff --git a/Exec/ABL/inputs_GABLS1_mynn25 b/Exec/ABL/inputs_GABLS1_mynn25 index 6806e41d7..b6f61e4a8 100644 --- a/Exec/ABL/inputs_GABLS1_mynn25 +++ b/Exec/ABL/inputs_GABLS1_mynn25 @@ -1,33 +1,32 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -stop_time = 32400.0 # 540 min = 9 h -stop_time = 600.0 # shorter tests +stop_time = 32400.0 # 540 min = 9 h (Cuxart et al. 2006) -amrex.fpe_trap_invalid = 0 +amrex.fpe_trap_invalid = 1 fabarray.mfiter_tile_size = 1024 1024 1024 -# PROBLEM SIZE & GEOMETRY (dx=dy=dz=3.125 m) -geometry.prob_extent = 25 25 400 -amr.n_cell = 4 4 64 -geometry.is_periodic = 1 1 0 +# PROBLEM SIZE & GEOMETRY (Cuxart et al. 2006) +geometry.prob_extent = 25 25 400 +amr.n_cell = 4 4 64 + +geometry.is_periodic = 1 1 0 # MOST BOUNDARY (DEFAULT IS ADIABATIC FOR THETA) zlo.type = "Most" -erf.most.z0 = 0.1 # from Beare et al. 2006 -# erf.most.zref = 3.125 # == dz/2 +erf.most.z0 = 0.1 # from Cuxart et al. 2006 erf.most.surf_temp = 265.0 # initial value, should match input_sounding -erf.most.surf_heating_rate = -0.25 # [K/h] from Beare et al. 2006 +erf.most.surf_heating_rate = -0.25 # [K/h] from Cuxart et al. 2006 zhi.type = "SlipWall" zhi.theta_grad = 0.01 # [K/m] to match the input sounding -# INITIALIZATION (Beare et al. 2006) +# INITIALIZATION (Cuxart et al. 2006) erf.init_type = "input_sounding" erf.init_sounding_ideal = 1 erf.input_sounding_file = "input_sounding_GABLS1" # TIME STEP CONTROL -erf.fixed_dt = 1.0 +erf.fixed_dt = 1.0 # largest stable low Mach dt erf.fixed_mri_dt_ratio = 6 # DIAGNOSTICS & VERBOSITY @@ -35,65 +34,46 @@ erf.sum_interval = 1 # timesteps between computing mass erf.v = 1 # verbosity in ERF.cpp amr.v = 1 # verbosity in Amr.cpp -# planar averaging (with erf.v = 1) -erf.data_log = surf_hist.dat mean_profiles.dat covar_profiles.dat sfs_profiles.dat -erf.profile_int = 20 -erf.interp_profiles_to_cc = false - # REFINEMENT / REGRIDDING amr.max_level = 0 # maximum level number allowed # CHECKPOINT FILES erf.check_file = chk # root name of checkpoint file -erf.check_int = 6000 # number of timesteps between checkpoints +erf.check_int = -1 # number of timesteps between checkpoints # PLOTFILES -erf.plot_file_1 = plt # prefix of plotfile name -erf.plot_int_1 = 6000 # number of timesteps between plotfiles -erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta -erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta rhoQKE Kmv Khv z_phys +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 300 # number of timesteps between plotfiles +erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta rhoQKE Kmv Khv + # SOLVER CHOICE -erf.dycore_horiz_adv_type = "Upwind_3rd" erf.dycore_vert_adv_type = "Upwind_3rd" -erf.dryscal_horiz_adv_type = "Upwind_3rd" erf.dryscal_vert_adv_type = "Upwind_3rd" erf.molec_diff_type = "None" erf.use_gravity = true -# Coriolis parameter f = 1.39e-4 s^-1, from Beare et al. 2006 +# Coriolis parameter f = 1.39e-4 s^-1 (Cuxart et al. 2006) erf.use_coriolis = true erf.latitude = 73.0 erf.rotational_time_period = 86455.2516813368 -# Geostrophic wind (Beare et al. 2006) +# Geostrophic wind (Cuxart et al. 2006) erf.abl_driver_type = "GeostrophicWind" erf.abl_geo_wind = 8.0 0.0 0.0 # Turbulence closure -erf.les_type = "Smagorinsky" -erf.Cs = 0.17 # following Gadde et al 2020 BLM -#erf.les_type = "Deardorff" -#erf.Ck = 0.1 # Coefficient in Moeng1984, Eqn. 19 -#erf.Ce = 0.7 # Note: Ce_lcoeff = C_e - 1.9*C_k -#erf.Ce_wall = 3.9 # To account for "wall effects" -#erf.sigma_k = 0.5 # TKE diffusion coeff = 2*Km = Km * inv_sigma_k, Moeng1984, Eqn. 15 +erf.les_type = "None" # NOT USED -#erf.rho0_trans = 1.3223 # from Beare et al. 2006 -#erf.theta_ref = 263.5 # from Beare et al. 2006 +#erf.rho0_trans = 1.3223 # from Cuxart et al. 2006 +#erf.theta_ref = 263.5 # from Cuxart et al. 2006 -#erf.use_rayleigh_damping = true -#prob.dampcoef = 0.2 # [1/s] following FastEddy -#prob.zdamp = 100. # from Beare et al. 2006, "most models applied gravity wave damping above 300 m" +erf.pbl_type = "MYNN2.5" # Initial conditions from Beare et al. 2006 -#prob.KE_0 = 0.4 # [m2/s2] -#prob.KE_decay_height = 250. # [m] -#prob.KE_decay_order = 1 -#prob.T_0_Pert_Mag = 0.1 # [K] -#prob.pert_ref_height = 50. # [m] - -erf.pbl_type = "MYNN2.5" \ No newline at end of file +prob.KE_0 = 0.4 # [m2/s2] +prob.KE_decay_height = 250. # [m] +prob.KE_decay_order = 1 diff --git a/Exec/ABL/inputs_GABLS1_ysu b/Exec/ABL/inputs_GABLS1_ysu index cdce595f8..05dcb0169 100644 --- a/Exec/ABL/inputs_GABLS1_ysu +++ b/Exec/ABL/inputs_GABLS1_ysu @@ -1,33 +1,32 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -stop_time = 32400.0 # 540 min = 9 h -stop_time = 600.0 # shorter tests +stop_time = 32400.0 # 540 min = 9 h (Cuxart et al. 2006) -amrex.fpe_trap_invalid = 0 +amrex.fpe_trap_invalid = 1 fabarray.mfiter_tile_size = 1024 1024 1024 -# PROBLEM SIZE & GEOMETRY (dx=dy=dz=3.125 m) -geometry.prob_extent = 25 25 400 -amr.n_cell = 4 4 64 -geometry.is_periodic = 1 1 0 +# PROBLEM SIZE & GEOMETRY (Cuxart et al. 2006) +geometry.prob_extent = 25 25 400 +amr.n_cell = 4 4 64 + +geometry.is_periodic = 1 1 0 # MOST BOUNDARY (DEFAULT IS ADIABATIC FOR THETA) zlo.type = "Most" -erf.most.z0 = 0.1 # from Beare et al. 2006 -# erf.most.zref = 3.125 # == dz/2 +erf.most.z0 = 0.1 # from Cuxart et al. 2006 erf.most.surf_temp = 265.0 # initial value, should match input_sounding -erf.most.surf_heating_rate = -0.25 # [K/h] from Beare et al. 2006 +erf.most.surf_heating_rate = -0.25 # [K/h] from Cuxart et al. 2006 zhi.type = "SlipWall" zhi.theta_grad = 0.01 # [K/m] to match the input sounding -# INITIALIZATION (Beare et al. 2006) +# INITIALIZATION (Cuxart et al. 2006) erf.init_type = "input_sounding" erf.init_sounding_ideal = 1 erf.input_sounding_file = "input_sounding_GABLS1" # TIME STEP CONTROL -erf.fixed_dt = 1.0 +erf.fixed_dt = 1.0 # largest stable low Mach dt erf.fixed_mri_dt_ratio = 6 # DIAGNOSTICS & VERBOSITY @@ -35,65 +34,46 @@ erf.sum_interval = 1 # timesteps between computing mass erf.v = 1 # verbosity in ERF.cpp amr.v = 1 # verbosity in Amr.cpp -# planar averaging (with erf.v = 1) -erf.data_log = surf_hist.dat mean_profiles.dat covar_profiles.dat sfs_profiles.dat -erf.profile_int = 20 -erf.interp_profiles_to_cc = false - # REFINEMENT / REGRIDDING amr.max_level = 0 # maximum level number allowed # CHECKPOINT FILES erf.check_file = chk # root name of checkpoint file -erf.check_int = 6000 # number of timesteps between checkpoints +erf.check_int = -1 # number of timesteps between checkpoints # PLOTFILES -erf.plot_file_1 = plt # prefix of plotfile name -erf.plot_int_1 = 6000 # number of timesteps between plotfiles -erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta -erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta rhoQKE Kmv Khv z_phys +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 300 # number of timesteps between plotfiles +erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure theta rhoQKE Kmv Khv + # SOLVER CHOICE -erf.dycore_horiz_adv_type = "Upwind_3rd" erf.dycore_vert_adv_type = "Upwind_3rd" -erf.dryscal_horiz_adv_type = "Upwind_3rd" erf.dryscal_vert_adv_type = "Upwind_3rd" erf.molec_diff_type = "None" erf.use_gravity = true -# Coriolis parameter f = 1.39e-4 s^-1, from Beare et al. 2006 +# Coriolis parameter f = 1.39e-4 s^-1 (Cuxart et al. 2006) erf.use_coriolis = true erf.latitude = 73.0 erf.rotational_time_period = 86455.2516813368 -# Geostrophic wind (Beare et al. 2006) +# Geostrophic wind (Cuxart et al. 2006) erf.abl_driver_type = "GeostrophicWind" erf.abl_geo_wind = 8.0 0.0 0.0 # Turbulence closure -erf.les_type = "Smagorinsky" -erf.Cs = 0.17 # following Gadde et al 2020 BLM -#erf.les_type = "Deardorff" -#erf.Ck = 0.1 # Coefficient in Moeng1984, Eqn. 19 -#erf.Ce = 0.7 # Note: Ce_lcoeff = C_e - 1.9*C_k -#erf.Ce_wall = 3.9 # To account for "wall effects" -#erf.sigma_k = 0.5 # TKE diffusion coeff = 2*Km = Km * inv_sigma_k, Moeng1984, Eqn. 15 +erf.les_type = "None" # NOT USED -#erf.rho0_trans = 1.3223 # from Beare et al. 2006 -#erf.theta_ref = 263.5 # from Beare et al. 2006 +#erf.rho0_trans = 1.3223 # from Cuxart et al. 2006 +#erf.theta_ref = 263.5 # from Cuxart et al. 2006 -#erf.use_rayleigh_damping = true -#prob.dampcoef = 0.2 # [1/s] following FastEddy -#prob.zdamp = 100. # from Beare et al. 2006, "most models applied gravity wave damping above 300 m" +erf.pbl_type = "YSU" # Initial conditions from Beare et al. 2006 -#prob.KE_0 = 0.4 # [m2/s2] -#prob.KE_decay_height = 250. # [m] -#prob.KE_decay_order = 1 -#prob.T_0_Pert_Mag = 0.1 # [K] -#prob.pert_ref_height = 50. # [m] - -erf.pbl_type = "YSU" \ No newline at end of file +prob.KE_0 = 0.4 # [m2/s2] +prob.KE_decay_height = 250. # [m] +prob.KE_decay_order = 1 From ef89e82282192c58a9c97de3806d8f9bb7b07acd Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Fri, 12 Jul 2024 11:04:44 -0600 Subject: [PATCH 7/9] working stable ysu --- Exec/ABL/inputs_GABLS1_mynn25 | 2 +- Exec/ABL/inputs_GABLS1_ysu | 3 +- Source/Diffusion/PBLModels.cpp | 98 +++++++++++++++++++++++++++++++--- 3 files changed, 93 insertions(+), 10 deletions(-) diff --git a/Exec/ABL/inputs_GABLS1_mynn25 b/Exec/ABL/inputs_GABLS1_mynn25 index b6f61e4a8..4dd87af26 100644 --- a/Exec/ABL/inputs_GABLS1_mynn25 +++ b/Exec/ABL/inputs_GABLS1_mynn25 @@ -1,7 +1,7 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- stop_time = 32400.0 # 540 min = 9 h (Cuxart et al. 2006) -amrex.fpe_trap_invalid = 1 +amrex.fpe_trap_invalid = 0 fabarray.mfiter_tile_size = 1024 1024 1024 diff --git a/Exec/ABL/inputs_GABLS1_ysu b/Exec/ABL/inputs_GABLS1_ysu index 05dcb0169..00db9c5c5 100644 --- a/Exec/ABL/inputs_GABLS1_ysu +++ b/Exec/ABL/inputs_GABLS1_ysu @@ -1,7 +1,7 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- stop_time = 32400.0 # 540 min = 9 h (Cuxart et al. 2006) -amrex.fpe_trap_invalid = 1 +amrex.fpe_trap_invalid = 0 fabarray.mfiter_tile_size = 1024 1024 1024 @@ -16,6 +16,7 @@ zlo.type = "Most" erf.most.z0 = 0.1 # from Cuxart et al. 2006 erf.most.surf_temp = 265.0 # initial value, should match input_sounding erf.most.surf_heating_rate = -0.25 # [K/h] from Cuxart et al. 2006 +erf.most.zref = 10.0 zhi.type = "SlipWall" zhi.theta_grad = 0.01 # [K/m] to match the input sounding diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp index 46b44b223..69f0dea68 100644 --- a/Source/Diffusion/PBLModels.cpp +++ b/Source/Diffusion/PBLModels.cpp @@ -270,14 +270,12 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel, // Get some data in arrays const auto& cell_data = cons_in.const_array(mfi); - //const auto& K_turb = eddyViscosity.array(mfi); const auto& uvel = xvel.const_array(mfi); const auto& vvel = yvel.const_array(mfi); const auto& z0_arr = most->get_z0(level)->const_array(); const auto& ws10av_arr = most->get_mac_avg(level,4)->const_array(mfi); const auto& t10av_arr = most->get_mac_avg(level,2)->const_array(mfi); - //const auto& u_star_arr = most->get_u_star(level)->const_array(mfi); //const auto& t_star_arr = most->get_t_star(level)->const_array(mfi); const auto& t_surf_arr = most->get_t_surf(level)->const_array(mfi); //const auto& l_obuk_arr = most->get_olen(level)->const_array(mfi); @@ -303,9 +301,11 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel, const GeometryData gdata = geom.data(); const Box xybx = PerpendicularBox(bx, IntVect{0,0,0}); FArrayBox pbl_height(xybx,1); + IArrayBox pbl_index(xybx,1); const auto& pblh_arr = pbl_height.array(); + const auto& pbli_arr = pbl_index.array(); - // -- Diagnose PBL height - starting out assuming non-moist + // -- Diagnose PBL height - starting out assuming non-moist -- // loop is only over i,j in order to find height at each x,y const Real f0 = turbChoice.pbl_ysu_coriolis_freq; const bool over_land = turbChoice.pbl_ysu_over_land; // TODO: make this local and consistent @@ -368,21 +368,103 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel, if (pblh_arr(i,j,0) < 0.5*(zval_up+zval_dn) ) { kpbl = 0; } + pbli_arr(i,j,0) = kpbl; // FIXME: DEBUG std::cout << i << " " << j << " " << pblh_arr(i,j,0) << std::endl; }); - // -- Compute nonlocal/countergradient mixing parameters + // -- Compute nonlocal/countergradient mixing parameters -- + // Not included for stable so nothing to do until unstable treatment is added - // -- Compute entrainment parameters + // -- Compute entrainment parameters -- + // 0 for stable so nothing to do? - // -- Compute diffusion coefficients within PBL + // -- Compute diffusion coefficients -- - // -- Compute coefficients in free stream above PBL + const auto& u_star_arr = most->get_u_star(level)->const_array(mfi); + const auto& l_obuk_arr = most->get_olen(level)->const_array(mfi); + const Array4 &K_turb = eddyViscosity.array(mfi); + + // Dirichlet flags to switch derivative stencil + bool c_ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir) ); + bool c_ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc].lo(5) == ERFBCType::ext_dir) ); + bool u_ext_dir_on_zlo = ( (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir) ); + bool u_ext_dir_on_zhi = ( (bc_ptr[BCVars::xvel_bc].lo(5) == ERFBCType::ext_dir) ); + bool v_ext_dir_on_zlo = ( (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir) ); + bool v_ext_dir_on_zhi = ( (bc_ptr[BCVars::yvel_bc].lo(5) == ERFBCType::ext_dir) ); + + const auto& dxInv = geom.InvCellSizeArray(); + const Real dz_inv = geom.InvCellSize(2); + const int izmin = geom.Domain().smallEnd(2); + const int izmax = geom.Domain().bigEnd(2); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + const Real zval = use_terrain ? Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr) : gdata.ProbLo(2) + (k + 0.5)*gdata.CellSize(2); + const Real rho = cell_data(i,j,k,Rho_comp); + const Real met_h_zeta = use_terrain ? Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd_arr) : 1.0; + const Real dz_terrain = met_h_zeta/dz_inv; + if (k < pbli_arr(i,j,0)) { + // -- Compute diffusion coefficients within PBL + constexpr Real zfacmin = 1e-8; // value from WRF + constexpr Real phifac = 8.0; // value from H10 and WRF + constexpr Real wstar3 = 0.0; // only nonzero for unstable + constexpr Real pfac = 2.0; // profile exponent + const Real zfac = std::min(std::max(1 - zval / pblh_arr(i,j,0), zfacmin ), 1.0); + // Not including YSU top down PBL term (not in H10, added to WRF later) + const Real ust3 = u_star_arr(i,j,0) * u_star_arr(i,j,0) * u_star_arr(i,j,0); + Real wscalek = ust3 + phifac * KAPPA * wstar3 * (1.0 - zfac); + wscalek = std::pow(wscalek, 1.0/3.0); + // stable only + const Real phi_term = 1 + 5 * zval / l_obuk_arr(i,j,0); // phi_term appears in WRF but not papers + wscalek = std::max(u_star_arr(i,j,0) / phi_term, 0.001); // 0.001 limit appears in WRF but not papers + K_turb(i,j,k,EddyDiff::Mom_v) = rho * wscalek * KAPPA * zval * std::pow(zfac, pfac); + K_turb(i,j,k,EddyDiff::Theta_v) = K_turb(i,j,k,EddyDiff::Mom_v); + } else { + // -- Compute coefficients in free stream above PBL + constexpr Real lam0 = 30.0; + constexpr Real min_richardson = -100.0; + constexpr Real prandtl_max = 4.0; + Real dthetadz, dudz, dvdz; + ComputeVerticalDerivativesPBL(i, j, k, + uvel, vvel, cell_data, izmin, izmax, 1.0/dz_terrain, + c_ext_dir_on_zlo, c_ext_dir_on_zhi, + u_ext_dir_on_zlo, u_ext_dir_on_zhi, + v_ext_dir_on_zlo, v_ext_dir_on_zhi, + dthetadz, dudz, dvdz); + const Real shear_squared = dudz*dudz + dvdz*dvdz + 1.0e-9; // 1.0e-9 from WRF to avoid divide by zero + const Real theta = cell_data(i,j,k,RhoTheta_comp) / cell_data(i,j,k,Rho_comp); + Real richardson = CONST_GRAV / theta * dthetadz / shear_squared; + const Real lambdadz = std::min(std::max(0.1*dz_terrain , lam0), 300.0); // in WRF, H10 paper just says use lam0 + const Real lengthscale = lambdadz * KAPPA * zval / (lambdadz + KAPPA * zval); + const Real turbfact = lengthscale * lengthscale * std::sqrt(shear_squared); + + if (richardson < 0) { + richardson = max(richardson, min_richardson); + Real sqrt_richardson = std::sqrt(-richardson); + K_turb(i,j,k,EddyDiff::Mom_v) = rho * turbfact * (1.0 - 8.0 * richardson / (1.0 + 1.746 * sqrt_richardson)); + K_turb(i,j,k,EddyDiff::Theta_v) = rho * turbfact * (1.0 - 8.0 * richardson / (1.0 + 1.286 * sqrt_richardson)); + } else { + const Real oneplus5ri = 1.0 + 5.0 * richardson; + K_turb(i,j,k,EddyDiff::Theta_v) = rho * turbfact / (oneplus5ri * oneplus5ri); + const Real prandtl = std::min(1.0+2.1*richardson, prandtl_max); // limit from WRF + K_turb(i,j,k,EddyDiff::Mom_v) = K_turb(i,j,k,EddyDiff::Theta_v) * prandtl; + } + } + + // limit both diffusion coefficients - from WRF, not documented in papers + constexpr Real ckz = 0.001; + constexpr Real Kmax = 1000.0; + const Real rhoKmin = ckz * dz_terrain * rho; + const Real rhoKmax = rho * Kmax; + K_turb(i,j,k,EddyDiff::Mom_v) = std::max(std::min(K_turb(i,j,k,EddyDiff::Mom_v) ,rhoKmax), rhoKmin); + K_turb(i,j,k,EddyDiff::Mom_v) = std::max(std::min(K_turb(i,j,k,EddyDiff::Theta_v) ,rhoKmax), rhoKmin); + K_turb(i,j,k,EddyDiff::PBL_lengthscale) = pblh_arr(i,j,0); + }); // FIXME: DEBUG - Error("YSU Model not implemented yet"); + // Error("YSU Model not implemented yet"); } } From 49b20ee5b7f8cd2658a6f71b15f3d34dab983a97 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Fri, 12 Jul 2024 11:13:49 -0600 Subject: [PATCH 8/9] remove debug print statement --- Source/Diffusion/PBLModels.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp index 69f0dea68..cffaf33bf 100644 --- a/Source/Diffusion/PBLModels.cpp +++ b/Source/Diffusion/PBLModels.cpp @@ -369,9 +369,6 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel, kpbl = 0; } pbli_arr(i,j,0) = kpbl; - - // FIXME: DEBUG - std::cout << i << " " << j << " " << pblh_arr(i,j,0) << std::endl; }); // -- Compute nonlocal/countergradient mixing parameters -- From b23c0e13e66514f89e8f34fa2d6906834893c125 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Fri, 12 Jul 2024 15:13:25 -0600 Subject: [PATCH 9/9] remove some comments from ysu pbl --- Source/Diffusion/PBLModels.cpp | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp index cffaf33bf..f74f8355f 100644 --- a/Source/Diffusion/PBLModels.cpp +++ b/Source/Diffusion/PBLModels.cpp @@ -245,16 +245,13 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel, }); } } else if (turbChoice.pbl_type == PBLType::YSU) { - // FIXME - // Error("YSU Model not implemented yet"); - /* YSU PBL initially introduced by S.-Y. Hong, Y. Noh, and J. Dudhia, MWR, 2006 [HND06] Further Modifications from S.-Y. Hong, Q. J. R. Meteorol. Soc., 2010 [H10] Implementation follows WRF as of early 2024 with some simplifications - */ + */ #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -276,9 +273,7 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel, const auto& z0_arr = most->get_z0(level)->const_array(); const auto& ws10av_arr = most->get_mac_avg(level,4)->const_array(mfi); const auto& t10av_arr = most->get_mac_avg(level,2)->const_array(mfi); - //const auto& t_star_arr = most->get_t_star(level)->const_array(mfi); const auto& t_surf_arr = most->get_t_surf(level)->const_array(mfi); - //const auto& l_obuk_arr = most->get_olen(level)->const_array(mfi); const bool use_terrain = (z_phys_nd != nullptr); const Array4 z_nd_arr = use_terrain ? z_phys_nd->array(mfi) : Array4{}; const Real most_zref = most->get_zref(); @@ -459,10 +454,6 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel, K_turb(i,j,k,EddyDiff::Mom_v) = std::max(std::min(K_turb(i,j,k,EddyDiff::Theta_v) ,rhoKmax), rhoKmin); K_turb(i,j,k,EddyDiff::PBL_lengthscale) = pblh_arr(i,j,0); }); - - // FIXME: DEBUG - // Error("YSU Model not implemented yet"); } - } }