fix some bugs in YSU

erf-model · Jul 17, 2024 · f2d5463 · f2d5463
1 parent 133fe11
commit f2d5463
Show file tree

Hide file tree

Showing 2 changed files with 15 additions and 3 deletions.
diff --git a/Exec/ABL/inputs_GABLS1_ysu b/Exec/ABL/inputs_GABLS1_ysu
@@ -45,7 +45,7 @@ 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
+erf.plot_vars_1     = density x_velocity y_velocity z_velocity pressure theta rhoQKE Kmv Khv Lpbl
 
 
 # SOLVER CHOICE

diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp
@@ -359,7 +359,10 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel,
                 const Real zval_up = use_terrain ? Compute_Zrel_AtCellCenter(i,j,kpbl,z_nd_arr) : gdata.ProbLo(2) + (kpbl + 0.5)*gdata.CellSize(2);
                 const Real zval_dn = use_terrain ? Compute_Zrel_AtCellCenter(i,j,kpbl-1,z_nd_arr) : gdata.ProbLo(2) + (kpbl-1 + 0.5)*gdata.CellSize(2);
                 pblh_arr(i,j,0) = zval_dn + interp_fact*(zval_up-zval_dn);
-                if (pblh_arr(i,j,0) < 0.5*(zval_up+zval_dn) ) {
+
+                const Real zval_0 = use_terrain ? Compute_Zrel_AtCellCenter(i,j,0,z_nd_arr) : gdata.ProbLo(2) + (0.5)*gdata.CellSize(2);
+                const Real zval_1 = use_terrain ? Compute_Zrel_AtCellCenter(i,j,1,z_nd_arr) : gdata.ProbLo(2) + (1.5)*gdata.CellSize(2);
+                if (pblh_arr(i,j,0) < 0.5*(zval_0+zval_1) ) {
                     kpbl = 0;
                 }
                 pbli_arr(i,j,0) = kpbl;
@@ -450,9 +453,18 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel,
                 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::Theta_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);
             });
+
+            // HACK set bottom ghost cell to 1st cell
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                if (k==-1) {
+                    K_turb(i,j,k,EddyDiff::Mom_v) = K_turb(i,j,0,EddyDiff::Mom_v);
+                    K_turb(i,j,k,EddyDiff::Theta_v) = K_turb(i,j,0,EddyDiff::Theta_v);
+                }
+            });
         }
     }
 }