new profile with density at vector locations

src/cases/Anelastic.hpp

@@ -122,6 +122,14 @@ namespace cases
       profs.rhod = rho_surf * exp(- st_avg * k * dz) * pow(
                1. - cs * (1 - exp(- st_avg * k * dz)), (1. / R_d_over_c_pd<real_t>()) - 1);
 
+      // rhod profile at vector field positions (midway between rhod)
+      // linear interpolation (extrapolation at the edges) from rhod
+      for(int i=1; i<nz; ++i)
+      {
+        profs.rhod_vctr(i) = (profs.rhod(i-1) + profs.rhod(i)) / 2.;
+      }
+      profs.rhod_vctr(0) = profs.rhod(0) - (profs.rhod(1) - profs.rhod(0)) / 2.;
+      profs.rhod_vctr(nz) = profs.rhod(nz-1) + (profs.rhod(nz-1) - profs.rhod(nz-2)) / 2.;
 
       // theta_std env prof to theta_dry_e
 //      for(int k=1; k<nz; ++k)

src/cases/DryThermalGMD2015.hpp

@@ -93,6 +93,7 @@ namespace cases
         parent_t::set_profs(profs, nz, user_params);
 
         profs.rhod = 1;
+        profs.rhod_vctr = 1;
         profs.th_e = 300;
         profs.rv_e = 0; // doesnt matter in dry case, just to have consistent output between runs
 

src/cases/MoistThermalGrabowskiClark99.hpp

@@ -287,6 +287,7 @@ namespace cases
 
         //th_ref = th_std_fctr(th_std_0 / si::kelvins)(k * dz);
         profs.rhod = rho_fctr(rhod_surf)(k * dz); // rhod is dry density profsile?
+        profs.rhod_vctr = rho_fctr(rhod_surf)((k-0.5) * dz); // rhod is dry density profsile?
 
         // turn virtual potential temperature env profsile into env profsile of standard potential temp
         profs.th_e = profs.th_e / (1. + a * profs.rv_e);

src/detail/profiles.hpp

@@ -10,13 +10,12 @@ namespace detail
   // TODO: store profiles in mem (add an alloc 1d function to libmpdata?) and not in params
   struct profiles_t
   {
-    setup::arr_1D_t th_e, p_e, rv_e, rl_e, th_ref, rhod, w_LS, hgt_fctr, th_LS, rv_LS, mix_len, relax_th_rv_coeff;
+    setup::arr_1D_t th_e, p_e, rv_e, rl_e, th_ref, rhod, rhod_vctr, w_LS, hgt_fctr, th_LS, rv_LS, mix_len, relax_th_rv_coeff;
     std::array<setup::arr_1D_t, 2> geostr;
 
     profiles_t(int nz) :
-    // rhod needs to be bigger, cause it divides vertical courant number
-    // TODO: should have a halo both up and down, not only up like now; then it should be interpolated in courant calculation
-      th_e(nz), p_e(nz), rv_e(nz), rl_e(nz), th_ref(nz), rhod(nz+1), w_LS(nz), hgt_fctr(nz), th_LS(nz), rv_LS(nz), mix_len(nz), relax_th_rv_coeff(nz)
+    // rhod_vctr is bigger, because it is density at the position of vector fields (Arakawa-C)
+      th_e(nz), p_e(nz), rv_e(nz), rl_e(nz), th_ref(nz), rhod(nz), rhod_vctr(nz+1), w_LS(nz), hgt_fctr(nz), th_LS(nz), rv_LS(nz), mix_len(nz), relax_th_rv_coeff(nz)
     {
       geostr[0].resize(nz);
       geostr[1].resize(nz);
@@ -32,7 +31,7 @@ namespace detail
 
   struct profile_ptrs_t
   {
-    setup::arr_1D_t *th_e, *p_e, *rv_e, *rl_e, *th_ref, *rhod, *w_LS, *hgt_fctr, *geostr[2], *th_LS, *rv_LS, *mix_len, *relax_th_rv_coeff;
+    setup::arr_1D_t *th_e, *p_e, *rv_e, *rl_e, *th_ref, *rhod, *rhod_vctr, *w_LS, *hgt_fctr, *geostr[2], *th_LS, *rv_LS, *mix_len, *relax_th_rv_coeff;
   };
 
   // copy external profiles into rt_parameters
@@ -41,19 +40,20 @@ namespace detail
   inline void copy_profiles(profiles_t &profs, params_t &p)
   {
     std::vector<std::pair<std::reference_wrapper<setup::arr_1D_t*>, std::reference_wrapper<setup::arr_1D_t>>> tobecopied = {
-      {p.hgt_fctr      , profs.hgt_fctr      },
-      {p.th_e          , profs.th_e          },
-      {p.p_e           , profs.p_e           },
-      {p.rv_e          , profs.rv_e          },
-      {p.rl_e          , profs.rl_e          },
-      {p.th_ref        , profs.th_ref        },
-      {p.rhod          , profs.rhod          },
-      {p.w_LS          , profs.w_LS          },
-      {p.th_LS         , profs.th_LS         },
-      {p.rv_LS         , profs.rv_LS         },
-      {p.geostr[0]     , profs.geostr[0]     },
-      {p.geostr[1]     , profs.geostr[1]     },
-      {p.mix_len       , profs.mix_len       },
+      {p.hgt_fctr          , profs.hgt_fctr          },
+      {p.th_e              , profs.th_e              },
+      {p.p_e               , profs.p_e               },
+      {p.rv_e              , profs.rv_e              },
+      {p.rl_e              , profs.rl_e              },
+      {p.th_ref            , profs.th_ref            },
+      {p.rhod              , profs.rhod              },
+      {p.rhod_vctr         , profs.rhod_vctr         },
+      {p.w_LS              , profs.w_LS              },
+      {p.th_LS             , profs.th_LS             },
+      {p.rv_LS             , profs.rv_LS             },
+      {p.geostr[0]         , profs.geostr[0]         },
+      {p.geostr[1]         , profs.geostr[1]         },
+      {p.mix_len           , profs.mix_len           },
       {p.relax_th_rv_coeff , profs.relax_th_rv_coeff }
     };
 

src/run_hlpr.cpp

@@ -130,7 +130,6 @@ void run(const int (&nps)[n_dims], const user_params_t &user_params)
 
   // reference profiles shared among threads
   detail::profiles_t profs(nz); 
-  // rhod needs to be bigger, cause it divides vertical courant number, TODO: should have a halo both up and down, not only up like now; then it should be interpolated in courant calculation
 
   // assign their values
   case_ptr->set_profs(profs, nz, user_params);

src/solvers/lgrngn/hook_mixed_rhs_ante_step_lgrngn.hpp

@@ -27,11 +27,15 @@ void slvr_lgrngn<ct_params_t>::hook_mixed_rhs_ante_step()
     nancheck(Cy, "Cy after copying from mpdata");
     nancheck(Cz, "Cz after copying from mpdata");
 
-    // ... and now dividing them by this->rhod (TODO: z=0 is located at k=1/2)
+    // ... and now dividing them by density
     {
       Cx.reindex(this->zero) /= (*params.rhod)(this->vert_idx);
       Cy.reindex(this->zero) /= (*params.rhod)(this->vert_idx);
-      Cz.reindex(this->zero) /= (*params.rhod)(this->vert_idx); // TODO: should be interpolated, since theres a shift between positions of rhod and Cz
+      Cz.reindex(this->zero) /= (*params.rhod_vctr)(this->vert_idx); 
+      // NOTE: rhod_vctr is calculated by linear interpolation (extrapolation at the edges) of rhod;
+      //       in libmpdata++, mem->GC (courant * rhod) is linearly interpolated, hence
+      //       courant numbers themselves are not linearly interpolated, because
+      //       mem->GC != rhod_vctr * C_interp, where C_interp is a linearly interpolated courant number
     }
 
     // assuring previous async step finished ...

src/solvers/slvr_common.hpp

@@ -175,6 +175,7 @@ class slvr_common : public slvr_dim<ct_params_t>
       this->record_prof_const("rl_e", params.rl_e->data()); 
       this->record_prof_const("th_ref", params.th_ref->data()); 
       this->record_prof_const("rhod", params.rhod->data()); 
+      this->record_prof_vctr_const("rhod_vctr", params.rhod_vctr->data()); 
       this->record_prof_const("w_LS", params.w_LS->data()); 
       this->record_prof_const("th_LS", params.th_LS->data()); 
       this->record_prof_const("rv_LS", params.rv_LS->data());