QKE src fixes (#1252)

* This runs multiple steps without a negative QKE and the field looks reasonable. * Fix debug modification. * fix sqrt in brunt vaisala freq --------- Co-authored-by: Bruce Perry <[email protected]>
erf-model · Sep 19, 2023 · 114f2cf · 114f2cf
1 parent b16c3e0
commit 114f2cf
Show file tree

Hide file tree

Showing 3 changed files with 15 additions and 15 deletions.
diff --git a/Exec/DevTests/MiguelDev/inputs b/Exec/DevTests/MiguelDev/inputs
@@ -23,7 +23,7 @@ zhi.type = "SlipWall"
 
 # TIME STEP CONTROL
 erf.no_substepping = 0
-erf.fixed_dt       = 1     # fixed time step
+erf.fixed_dt       = 0.5     # fixed time step
 
 # DIAGNOSTICS & VERBOSITY
 erf.sum_interval   = 1       # timesteps between computing mass
@@ -40,8 +40,8 @@ erf.check_int       =-100        # number of timesteps between checkpoints
 # PLOTFILES
 erf.plotfile_type   = netcdf
 erf.plot_file_1     = plt        # prefix of plotfile name
-erf.plot_int_1      = 100         # number of timesteps between plotfiles
-erf.plot_vars_1     = density x_velocity y_velocity z_velocity pressure temp theta pres_hse
+erf.plot_int_1      = 100        # number of timesteps between plotfiles
+erf.plot_vars_1     = density rhotheta rhoQKE x_velocity y_velocity z_velocity pressure temp theta QKE pres_hse
 
 # SOLVER CHOICE
 erf.alpha_T = 0.0
@@ -55,7 +55,7 @@ erf.dynamicViscosity = 0.0 # 5.0
 erf.Cs = 0.16
 
 erf.dycore_horiz_adv_type = "Upwind_3rd"
-erf.dycore_vert_adv_type = "Upwind_5th"
+erf.dycore_vert_adv_type  = "Upwind_5th"
 
 erf.use_coriolis = true
 erf.abl_driver_type = "GeostrophicWind"
@@ -66,11 +66,11 @@ erf.use_rayleigh_damping = true #true
 
 # PROBLEM PARAMETERS (optional)
 prob.rho_0 = 1.0
-prob.T_0 = 300.0
-prob.Xc_0 = 100000
-prob.Yc_0 = 100000
-prob.VMAX = 15
-prob.RMAX = 20000
+prob.T_0   = 300.0
+prob.Xc_0  = 100000
+prob.Yc_0  = 100000
+prob.VMAX  = 15
+prob.RMAX  = 20000
 prob.RZERO = 80000
 prob.ZZERO = 3000
 prob.dampcoef = 0.2

diff --git a/Source/Diffusion/ComputeQKESourceTerm.H b/Source/Diffusion/ComputeQKESourceTerm.H
@@ -48,16 +48,16 @@ ComputeQKESourceTerms (int i, int j, int k,
     dvdz = 0.25*(vvel(i,j,k+1) - vvel(i,j,k-1) + vvel(i,j+1,k+1) - vvel(i,j+1,k-1))*dz_inv;
   }
 
-  // Bouyancy
-  source_term += 2*CONST_GRAV/theta_mean*K_turb(i,j,k,EddyDiff::Theta_v)*dthetadz;
-
   // Production
   source_term += K_turb(i,j,k,EddyDiff::Mom_v) * (dudz*dudz + dvdz*dvdz);
 
+  // Bouyancy
+  source_term -= 2*(CONST_GRAV/theta_mean)*K_turb(i,j,k,EddyDiff::Theta_v)*dthetadz;
+
   // Dissipation
   amrex::Real qke = cell_prim(i,j,k,PrimQKE_comp);
   if (std::abs(qke) > 0.0) {
-    source_term += 2.0 * cell_data(i,j,k,Rho_comp) * std::pow(qke,1.5) /
+    source_term -= 2.0 * cell_data(i,j,k,Rho_comp) * std::pow(qke,1.5) /
       (pbl_B1_l * K_turb(i,j,k,EddyDiff::PBL_lengthscale));
   }
 

diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp
@@ -154,7 +154,7 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
           // Third Length Scale
           amrex::Real l_B;
           if (dthetadz > 0) {
-              amrex::Real N_brunt_vaisala = CONST_GRAV/theta0 * std::sqrt(dthetadz);
+              amrex::Real N_brunt_vaisala = std::sqrt(CONST_GRAV/theta0 * dthetadz);
               if (zeta < 0) {
                   amrex::Real qc = CONST_GRAV/theta0 * surface_heat_flux * l_T;
                   qc = std::pow(qc,1.0/3.0);
@@ -172,7 +172,7 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
           // Compute non-dimensional parameters
           amrex::Real l2_over_q2 = l_comb*l_comb/(qvel(i,j,k)*qvel(i,j,k));
           amrex::Real GM = l2_over_q2 * (dudz*dudz + dvdz*dvdz);
-          amrex::Real GH = -l2_over_q2 / theta0 * dthetadz;
+          amrex::Real GH = -l2_over_q2 * (CONST_GRAV/theta0) * dthetadz;
           amrex::Real E1 = 1.0 + 6.0*A1*A1*GM - 9.0*A1*A2*(1.0-C2)*GH;
           amrex::Real E2 = -3.0*A1*(4.0*A1 + 3.0*A2*(1.0-C5))*(1.0-C2)*GH;
           amrex::Real E3 = 6.0*A2*A1*GM;