Compute Obukhov length based on virtual temperature/flux

ewquon · Aug 15, 2024 · c5cdacd · c5cdacd
1 parent ace1d2a
commit c5cdacd
Show file tree

Hide file tree

Showing 4 changed files with 82 additions and 38 deletions.
diff --git a/Source/BoundaryConditions/ABLMost.H b/Source/BoundaryConditions/ABLMost.H
@@ -217,7 +217,7 @@ public:
             t_star[lev]->setVal(1.E34);
 
             q_star[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
-            q_star[lev]->setVal(1.E34);
+            q_star[lev]->setVal(0.0); // default to dry
 
             olen[lev] = std::make_unique<amrex::MultiFab>(ba2d,dm,ncomp,ng);
             olen[lev]->setVal(1.E34);

diff --git a/Source/BoundaryConditions/ABLMost.cpp b/Source/BoundaryConditions/ABLMost.cpp
@@ -138,19 +138,24 @@ ABLMost::compute_fluxes (const int& lev,
                          bool is_land)
 {
     // Pointers to the computed averages
-    const auto *const tm_ptr  = m_ma.get_average(lev,2);
-    const auto *const umm_ptr = m_ma.get_average(lev,5);
+    const auto *const tm_ptr  = m_ma.get_average(lev,2); // potential temperature
+    const auto *const qvm_ptr = m_ma.get_average(lev,3); // water vapor mixing ratio
+    const auto *const tvm_ptr = m_ma.get_average(lev,4); // virtual potential temperature
+    const auto *const umm_ptr = m_ma.get_average(lev,5); // horizontal velocity magnitude
 
     for (MFIter mfi(*u_star[lev]); mfi.isValid(); ++mfi)
     {
         Box gtbx = mfi.growntilebox();
 
         auto u_star_arr = u_star[lev]->array(mfi);
         auto t_star_arr = t_star[lev]->array(mfi);
+        auto q_star_arr = q_star[lev]->array(mfi);
         auto t_surf_arr = t_surf[lev]->array(mfi);
         auto olen_arr   = olen[lev]->array(mfi);
 
         const auto tm_arr  = tm_ptr->array(mfi);
+        const auto tvm_arr = tvm_ptr->array(mfi);
+        const auto qvm_arr = qvm_ptr->array(mfi);
         const auto umm_arr = umm_ptr->array(mfi);
         const auto z0_arr  = z_0[lev].array();
 
@@ -164,13 +169,13 @@ ABLMost::compute_fluxes (const int& lev,
 
         ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
-            if (is_land && lmask_arr(i,j,k) == 1) {
-                most_flux.iterate_flux(i, j, k, max_iters, z0_arr, umm_arr, tm_arr,
-                                       u_star_arr, t_star_arr, t_surf_arr, olen_arr,
-                                       Hwave_arr, Lwave_arr, eta_arr);
-            } else if (!is_land && lmask_arr(i,j,k) == 0) {
-                most_flux.iterate_flux(i, j, k, max_iters, z0_arr, umm_arr, tm_arr,
-                                       u_star_arr, t_star_arr, t_surf_arr, olen_arr,
+            if (( is_land && lmask_arr(i,j,k) == 1) ||
+                (!is_land && lmask_arr(i,j,k) == 0))
+            {
+                most_flux.iterate_flux(i, j, k, max_iters,
+                                       z0_arr, umm_arr, tm_arr, tvm_arr, qvm_arr,
+                                       u_star_arr, t_star_arr, q_star_arr,  // to be updated
+                                       t_surf_arr, olen_arr,                // to be updated
                                        Hwave_arr, Lwave_arr, eta_arr);
             }
         });

diff --git a/Source/BoundaryConditions/MOSTAverage.cpp b/Source/BoundaryConditions/MOSTAverage.cpp
@@ -124,6 +124,9 @@ MOSTAverage::MOSTAverage (Vector<Geometry>  geom,
             m_averages[lev][iavg]->setVal(1.E34);
         }
 
+        // Default to dry
+        m_averages[lev][3]->setVal(0.0);
+
         if (m_rotate) {
             m_rot_fields[lev][2] = std::make_unique<MultiFab>(ba,dm,ncomp,ng);
             m_rot_fields[lev][3] = std::make_unique<MultiFab>(ba,dm,ncomp,ng);