Skip to content

Commit

Permalink
Merge pull request #1577 from baperry2/ysu-pbl-moreprog
Browse files Browse the repository at this point in the history
More WIP on YSU PBL
  • Loading branch information
baperry2 authored Jul 12, 2024
2 parents c381e25 + b23c0e1 commit 72bfc17
Show file tree
Hide file tree
Showing 7 changed files with 294 additions and 44 deletions.
File renamed without changes.
79 changes: 79 additions & 0 deletions Exec/ABL/inputs_GABLS1_mynn25
Original file line number Diff line number Diff line change
@@ -0,0 +1,79 @@
# ------------------ INPUTS TO MAIN PROGRAM -------------------
stop_time = 32400.0 # 540 min = 9 h (Cuxart et al. 2006)

amrex.fpe_trap_invalid = 0

fabarray.mfiter_tile_size = 1024 1024 1024

# 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 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

zhi.type = "SlipWall"
zhi.theta_grad = 0.01 # [K/m] to match the input sounding

# 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 # largest stable low Mach dt
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

# REFINEMENT / REGRIDDING
amr.max_level = 0 # maximum level number allowed

# CHECKPOINT FILES
erf.check_file = chk # root name of checkpoint file
erf.check_int = -1 # number of timesteps between checkpoints

# PLOTFILES
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_vert_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 (Cuxart et al. 2006)
erf.use_coriolis = true
erf.latitude = 73.0
erf.rotational_time_period = 86455.2516813368

# 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 = "None"

# NOT USED
#erf.rho0_trans = 1.3223 # from Cuxart et al. 2006
#erf.theta_ref = 263.5 # from Cuxart et al. 2006

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
80 changes: 80 additions & 0 deletions Exec/ABL/inputs_GABLS1_ysu
Original file line number Diff line number Diff line change
@@ -0,0 +1,80 @@
# ------------------ INPUTS TO MAIN PROGRAM -------------------
stop_time = 32400.0 # 540 min = 9 h (Cuxart et al. 2006)

amrex.fpe_trap_invalid = 0

fabarray.mfiter_tile_size = 1024 1024 1024

# 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 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

# 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 # largest stable low Mach dt
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

# REFINEMENT / REGRIDDING
amr.max_level = 0 # maximum level number allowed

# CHECKPOINT FILES
erf.check_file = chk # root name of checkpoint file
erf.check_int = -1 # number of timesteps between checkpoints

# PLOTFILES
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_vert_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 (Cuxart et al. 2006)
erf.use_coriolis = true
erf.latitude = 73.0
erf.rotational_time_period = 86455.2516813368

# 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 = "None"

# NOT USED
#erf.rho0_trans = 1.3223 # from Cuxart et al. 2006
#erf.theta_ref = 263.5 # from Cuxart et al. 2006

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
54 changes: 29 additions & 25 deletions Source/DataStructs/TurbStruct.H
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
3 changes: 2 additions & 1 deletion Source/Diffusion/DiffusionSrcForState_N.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,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);
Expand Down
Loading

0 comments on commit 72bfc17

Please sign in to comment.