Fractal refinement

CMakeLists.txt

@@ -73,6 +73,12 @@ if(UWLCM_TIMING)
   target_compile_options(uwlcm PRIVATE "-DUWLCM_TIMING")
 endif()
 
+# by default, don't compile 2D and blk, because they dont use refinement and libmpdata++ doesnt compile without refinement
+# also dont compile piggy as it is not compatible with refinement too
+if(NOT UWLCM_DISABLE)
+  set(UWLCM_DISABLE "2D_BLK_1M;3D_BLK_1M;2D_BLK_2M;3D_BLK_2M;2D_LGRNGN;2D_NONE;3D_NONE;PIGGYBACKER")
+endif()
+
 # handle the disable compilation options
 if(UWLCM_DISABLE)
   set(UWLCM_DISABLE_PREFIX "UWLCM_DISABLE_")

src/cases/Anelastic.hpp

@@ -21,7 +21,7 @@ namespace cases
     virtual quantity<si::dimensionless, real_t> r_t(const real_t &z) {throw std::runtime_error("base Anelastic class r_t called");}
 
     // calculate the initial environmental theta and rv profiles
-    // like in Wojtek's BabyEulag
+    // adapted from Wojtek Grabowski's BabyEulag
     void env_prof(detail::profiles_t &profs, int nz)
     {
       using libcloudphxx::common::moist_air::R_d_over_c_pd;
@@ -52,33 +52,57 @@ namespace cases
       real_t d = l_tri<real_t>() / si::joules * si::kilograms / rv;
       real_t f = R_d_over_c_pd<real_t>(); 
 
-      real_t lwp_env = 0;
-
+//      real_t lwp_env = 0;
+//
       for(int k=1; k<nz; ++k)
       {
-        real_t bottom = R_d<real_t>() / si::joules * si::kelvins * si::kilograms * T(k-1) * (1 + 0.61 * profs.rv_e(k-1)); // (p / rho) of moist air at k-1
-        real_t rho1 = profs.p_e(k-1) / bottom; // rho at k-1
-        profs.p_e(k) = profs.p_e(k-1) - rho1 * 9.81 * dz; // estimate of pre at k (dp = -g * rho * dz)
-        real_t thetme = pow(p_1000<real_t>() / si::pascals / profs.p_e(k), f); // 1/Exner
-        real_t thi = 1. / (th_l(k * dz) / si::kelvins); // 1/theta_std
-        real_t y = b * thetme * tt0 * thi; 
-        real_t ees = ee0 * exp(b-y); // saturation vapor pressure (Tetens equation or what?)
-        real_t qvs = a * ees / (profs.p_e(k) - ees);  // saturation vapor mixing ratio = R_d / R_v * ees / p_d
-        // calculate linearized condensation rate
-        real_t cf1 = thetme*thetme*thi*thi;  // T^{-2}
-        cf1 *= c * d * profs.p_e(k) / (profs.p_e(k) - ees); // = l_tri^2 / (C_pd * R_v * T^2) * p/p_d
-        real_t delta = (r_t(k*dz) - qvs) / (1 + qvs * cf1); // how much supersaturated is the air (divided by sth)
-        if(delta < 0.) delta = 0.;
-        profs.rv_e(k) = r_t(k*dz) - delta;
-        profs.rl_e(k) = delta;
-        profs.th_e(k) = th_l(k*dz) / si::kelvins + c * thetme * delta;
-        T(k) = profs.th_e(k) * pow(profs.p_e(k) / (p_1000<real_t>() / si::pascals),  f);
-
-        bottom = R_d<real_t>() / si::joules * si::kelvins * si::kilograms * T(k) * (1 + 0.61 * profs.rv_e(k)); // (p / rho) of moist air at k-1
-        rho1 = profs.p_e(k) / bottom; // rho at k-1
-        lwp_env  += delta * rho1;
+        /// estimate of pressure at k assuming constant pressure (dp = -g * rho * dz) (like in babyEulag)
+        /*
+        real_t Rd_Tv_m1 = R_d<real_t>() / si::joules * si::kelvins * si::kilograms * T(k-1) * (1 + 0.61 * profs.rv(k-1));
+        real_t rho_m1 = profs.p_e(k-1) / Rd_Tv_m1; // rho at k-1
+        profs.p_e(k) = profs.p_e(k-1) - rho1 * 9.81 * dz; 
+        */
+
+        real_t T_midlvl = T(k-1); // pressure calculating assuming constant temperature between levels, this is the first estimate of this constant temp
+        real_t rv_midlvl = profs.rv_e(k-1); // same for water vapor
+        real_t rl_midlvl = 0; // same for liquid water
+
+        for(int i=0; i<100; ++i) // TODO: why 100 iterations?
+        {
+          // pressure profile assuming constant temperature between levels
+          real_t Rd_Tv = R_d<real_t>() / si::joules * si::kelvins * si::kilograms * T_midlvl * (1 + 0.61 * rv_midlvl - rl_midlvl);
+          profs.p_e(k) = profs.p_e(k-1) * exp(-9.81 * dz / Rd_Tv ); // estimate of pre at k assuming constant temperature
+
+          // adjust for supersaturation
+          real_t thetme = pow(p_1000<real_t>() / si::pascals / profs.p_e(k), f); // 1/Exner
+          real_t thi = 1. / (th_l(k * dz) / si::kelvins); // 1/theta_std
+          real_t y = b * thetme * tt0 * thi; 
+          real_t ees = ee0 * exp(b-y); // saturation vapor pressure (Tetens equation or what?)
+          real_t qvs = a * ees / (profs.p_e(k) - ees);  // saturation vapor mixing ratio = R_d / R_v * ees / p_d
+          // calculate linearized condensation rate
+          real_t cf1 = thetme*thetme*thi*thi;  // T^{-2}
+          cf1 *= c * d * profs.p_e(k) / (profs.p_e(k) - ees); // = l_tri^2 / (C_pd * R_v * T^2) * p/p_d
+          real_t delta = (r_t(k*dz) - qvs) / (1 + qvs * cf1); // how much supersaturated is the air (divided by sth)
+          if(delta < 0.) delta = 0.;
+          profs.rv_e(k) = r_t(k*dz) - delta;
+          profs.rl_e(k) = delta;
+          profs.th_e(k) = th_l(k*dz) / si::kelvins + c * thetme * delta;
+          T(k) = profs.th_e(k) * pow(profs.p_e(k) / (p_1000<real_t>() / si::pascals),  f);
+
+          // linear interpolation between levels for the next iteration
+          T_midlvl = (T(k) + T(k-1)) / 2.;
+          rv_midlvl = (profs.rv_e(k) + profs.rv_e(k-1)) / 2.;
+          rl_midlvl = (profs.rl_e(k) + profs.rl_e(k-1)) / 2.;
+        }
+
+        /*
+        Rd_Tv = R_d<real_t>() / si::joules * si::kelvins * si::kilograms * T(k) * (1 + 0.61 * profs.rv_e(k)); // (p / rho) of moist air at k
+        rho_m1 = profs.p_e(k) / Rd_Tv; // rho at k
+        lwp_env  += delta * rho_m1;
+        */
       }
-      lwp_env = lwp_env * 5  * 1e3;
+    //  lwp_env = lwp_env * 5  * 1e3;
+    //
     }
 
     // calculate the initial reference theta and rv profiles
@@ -101,8 +125,8 @@ namespace cases
       st(notopbot) = (profs.th_e(notopbot+1) - profs.th_e(notopbot-1)) / profs.th_e(notopbot);
       real_t st_avg = blitz::sum(st) / (nz-2) / (2.*dz);
       // reference theta
-      profs.th_ref = profs.th_e(0) * (1. + 0.608 * profs.rv_e(0)) * exp(st_avg * k * dz);
-    //  th_ref = th_e(0) * pow(1 + rv_e(0) / a, f) // calc dry theta at z=0 
+      profs.th_reference = profs.th_e(0) * (1. + 0.608 * profs.rv_e(0)) * exp(st_avg * k * dz);
+    //  th_reference = th_e(0) * pow(1 + rv_e(0) / a, f) // calc dry theta at z=0 
     //           * exp(st_avg * k * dz);
       // virtual temp at surface
 

src/cases/CasesCommon.hpp

@@ -38,7 +38,7 @@ namespace cases
       return f_vel_prof(z) * si::seconds / si::meters - mean_vel;
     }
 
-    void init(bool window, quantity<si::length, real_t> Z) 
+    void init(bool window, quantity<si::length, real_t> Z, real_t target_mean_vel = 0) // target_mean_vel is the mean velocity we want to have within the window 
     {
       initialized=true;
       mean_vel = 0;
@@ -48,6 +48,7 @@ namespace cases
         for(int i=0; i < setup::mean_horvel_npts; ++i)
           mean_vel += f_vel_prof(i * (Z / si::meters) / (setup::mean_horvel_npts-1)) * si::seconds / si::meters;
         mean_vel /= setup::mean_horvel_npts;
+        mean_vel -= target_mean_vel;
       }
     }
 
@@ -158,32 +159,33 @@ namespace cases
     }
 
     virtual void update_surf_flux_sens(blitz::Array<real_t, n_dims> surf_flux_sens,
-                                       blitz::Array<real_t, n_dims> th_ground,   
-                                       blitz::Array<real_t, n_dims> U_ground,   
-                                       const real_t &U_ground_z,
-                                       const int &timestep, const real_t &dt, const real_t &dx, const real_t &dy = 0)
+                                       const blitz::Array<real_t, n_dims> &th_ground,   
+                                       const blitz::Array<real_t, n_dims> &U_ground,   
+                                       const int timestep, const real_t dt, 
+                                       const real_t dx, const real_t dy, const real_t U_ground_z)
     {if(timestep==0) surf_flux_sens = 0.;};
 
-    virtual void update_surf_flux_lat(blitz::Array<real_t, n_dims> surf_flux_lat,
-                                       blitz::Array<real_t, n_dims> rt_ground,   
-                                       blitz::Array<real_t, n_dims> U_ground,   
-                                       const real_t &U_ground_z,
-                                       const int &timestep, const real_t &dt, const real_t &dx, const real_t &dy = 0)
+    virtual void update_surf_flux_lat (blitz::Array<real_t, n_dims> surf_flux_lat,
+                                       const blitz::Array<real_t, n_dims> &rt_ground,   
+                                       const blitz::Array<real_t, n_dims> &U_ground,   
+                                       const int timestep, const real_t dt, 
+                                       const real_t dx, const real_t dy, const real_t U_ground_z)
     {if(timestep==0) surf_flux_lat = 0.;};
 
-    virtual void update_surf_flux_uv(blitz::Array<real_t, n_dims> surf_flux_uv,
-                                     blitz::Array<real_t, n_dims> uv_ground,   
-                                     blitz::Array<real_t, n_dims> U_ground,   
-                                     const real_t &U_ground_z,
-                                     const int &timestep, const real_t &dt, const real_t &dx, const real_t &dy = 0, const real_t &uv_mean = 0)
+    virtual void update_surf_flux_uv  (blitz::Array<real_t, n_dims> surf_flux_uv,
+                                       const blitz::Array<real_t, n_dims> &uv_ground,   
+                                       const blitz::Array<real_t, n_dims> &U_ground,   
+                                       const int timestep, const real_t dt, 
+                                       const real_t dx, const real_t dy, const real_t U_ground_z, 
+                                       const real_t uv_mean = 0)
     {if(timestep==0) surf_flux_uv = 0.;};
 
-    virtual void update_rv_LS(blitz::Array<real_t, 1> rv_LS,
-                              int timestep, real_t dt, real_t dz)
+    virtual void update_rv_LS(blitz::Array<real_t, 1> &rv_LS,
+                              const int timestep, const real_t dt, const real_t dz)
     {};
 
-    virtual void update_th_LS(blitz::Array<real_t, 1> th_LS,
-                              int timestep, real_t dt, real_t dz)
+    virtual void update_th_LS(blitz::Array<real_t, 1> &th_LS,
+                              const int timestep, const real_t dt, const real_t dz)
     {};
 
     // ctor

src/cases/CumulusCongestus.hpp

@@ -219,13 +219,12 @@ namespace cases
         profs.th_LS = 0.; // no large-scale horizontal advection
         profs.rv_LS = 0.; 
       }
-
       // functions that set surface fluxes per timestep
       void update_surf_flux_sens(blitz::Array<real_t, n_dims> surf_flux_sens,
-                                       blitz::Array<real_t, n_dims> th_ground,   
-                                       blitz::Array<real_t, n_dims> U_ground,   
-                                       const real_t &U_ground_z,
-                                       const int &timestep, const real_t &dt, const real_t &dx, const real_t &dy) override
+                                 const blitz::Array<real_t, n_dims> &th_ground,   
+                                 const blitz::Array<real_t, n_dims> &U_ground,   
+                                 const int timestep, const real_t dt, 
+                                 const real_t dx, const real_t dy, const real_t U_ground_z) override
       {
         if(timestep == 0) 
           surf_flux_sens = .1 * -1 * (this->rhod_0 / si::kilograms * si::cubic_meters) * theta_std::exner(p_0); // [K kg / (m^2 s)]; -1 because negative gradient of upward flux means inflow
@@ -238,12 +237,11 @@ namespace cases
         }
       }
 
-
       void update_surf_flux_lat(blitz::Array<real_t, n_dims> surf_flux_lat,
-                                       blitz::Array<real_t, n_dims> rt_ground,   
-                                       blitz::Array<real_t, n_dims> U_ground,   
-                                       const real_t &U_ground_z,
-                                       const int &timestep, const real_t &dt, const real_t &dx, const real_t &dy) override
+                                const blitz::Array<real_t, n_dims> &rt_ground,   
+                                const blitz::Array<real_t, n_dims> &U_ground,   
+                                const int timestep, const real_t dt, 
+                                const real_t dx, const real_t dy, const real_t U_ground_z) override
       {
         if(timestep == 0)
           surf_flux_lat = .4e-4 * -1 * (this->rhod_0 / si::kilograms * si::cubic_meters); // [m/s]
@@ -259,18 +257,18 @@ namespace cases
       // one function for updating u or v
       // the n_dims arrays have vertical extent of 1 - ground calculations only in here
       void update_surf_flux_uv(blitz::Array<real_t, n_dims>  surf_flux_uv, // output array
-                               blitz::Array<real_t, n_dims>  uv_ground,    // value of u or v on the ground
-                               blitz::Array<real_t, n_dims>  U_ground,     // magnitude of horizontal ground wind
-                               const real_t &U_ground_z,
-                               const int &timestep, const real_t &dt, const real_t &dx, const real_t &dy, const real_t &uv_mean)
+                               const blitz::Array<real_t, n_dims> &uv_ground,    // value of u or v on the ground
+                               const blitz::Array<real_t, n_dims> &U_ground,     // magnitude of horizontal ground wind
+                               const int timestep, const real_t dt,
+                               const real_t dx, const real_t dy, const real_t U_ground_z,
+                               const real_t uv_mean = 0)
       {
         surf_flux_uv = where(U_ground < 1e-4, 
             - 0.0784 * (uv_ground + uv_mean) / real_t(1e-4) * -1  * (this->rhod_0 / si::kilograms * si::cubic_meters), // 0.0784 m^2 / s^2 is the square of friction velocity = 0.28 m / s
             - 0.0784 * (uv_ground + uv_mean) / U_ground * -1  * (this->rhod_0 / si::kilograms * si::cubic_meters)
           );
       }
 
-
       // ctor
       CumulusCongestusCommon()
       {

src/cases/DYCOMS.hpp

@@ -258,10 +258,10 @@ namespace cases
       }
 
       void update_surf_flux_sens(blitz::Array<real_t, n_dims> surf_flux_sens,
-                                       blitz::Array<real_t, n_dims> th_ground,   
-                                       blitz::Array<real_t, n_dims> U_ground,   
-                                       const real_t &U_ground_z,
-                                       const int &timestep, const real_t &dt, const real_t &dx, const real_t &dy) override
+                                 const blitz::Array<real_t, n_dims> &th_ground,   
+                                 const blitz::Array<real_t, n_dims> &U_ground,   
+                                 const int timestep, const real_t dt, 
+                                 const real_t dx, const real_t dy, const real_t U_ground_z) override
       {
         if(timestep == 0) // TODO: what if this function is not called at t=0? force such call
         {
@@ -272,10 +272,10 @@ namespace cases
       }
 
       void update_surf_flux_lat(blitz::Array<real_t, n_dims> surf_flux_lat,
-                                       blitz::Array<real_t, n_dims> rt_ground,   
-                                       blitz::Array<real_t, n_dims> U_ground,   
-                                       const real_t &U_ground_z,
-                                       const int &timestep, const real_t &dt, const real_t &dx, const real_t &dy) override
+                                const blitz::Array<real_t, n_dims> &rt_ground,   
+                                const blitz::Array<real_t, n_dims> &U_ground,   
+                                const int timestep, const real_t dt, 
+                                const real_t dx, const real_t dy, const real_t U_ground_z) override
       {
         if(timestep == 0) // TODO: what if this function is not called at t=0? force such call
         {
@@ -287,11 +287,12 @@ namespace cases
 
       // one function for updating u or v
       // the n_dims arrays have vertical extent of 1 - ground calculations only in here
-      void update_surf_flux_uv(blitz::Array<real_t, n_dims>  surf_flux_uv, // output array
-                               blitz::Array<real_t, n_dims>  uv_ground,    // value of u or v on the ground
-                               blitz::Array<real_t, n_dims>  U_ground,     // magnitude of horizontal ground wind
-                               const real_t &U_ground_z,
-                               const int &timestep, const real_t &dt, const real_t &dx, const real_t &dy, const real_t &uv_mean)
+      void update_surf_flux_uv  (blitz::Array<real_t, n_dims> surf_flux_uv,
+                                 const blitz::Array<real_t, n_dims> &uv_ground,   
+                                 const blitz::Array<real_t, n_dims> &U_ground,   
+                                 const int timestep, const real_t dt, 
+                                 const real_t dx, const real_t dy, const real_t U_ground_z, 
+                                 const real_t uv_mean = 0) override
       {
         surf_flux_uv = where(U_ground == 0., 0.,
             - 0.0625 * (uv_ground + uv_mean) / U_ground * -1  * (this->rhod_0 / si::kilograms * si::cubic_meters)// 0.0625 m^2 / s^2 is the square of friction velocity = 0.25 m / s; * -1 because negative gradient of upward flux means inflow

src/cases/DryPBL.hpp

@@ -244,31 +244,23 @@ namespace cases
       }
 
       void update_surf_flux_sens(blitz::Array<real_t, n_dims> surf_flux_sens,
-                                 blitz::Array<real_t, n_dims> th_ground,    // value of th on the ground
-                                 blitz::Array<real_t, n_dims> U_ground,     // magnitude of horizontal ground wind
-                                 const real_t &U_ground_z,                   // altituted at which U_ground is diagnosed
-                                 const int &timestep, const real_t &dt, const real_t &dx, const real_t &dy) override
+                                 const blitz::Array<real_t, n_dims> &th_ground,   
+                                 const blitz::Array<real_t, n_dims> &U_ground,   
+                                 const int timestep, const real_t dt, 
+                                 const real_t dx, const real_t dy, const real_t U_ground_z) override
       {
         surf_flux_sens = .01 * -1 * (this->rhod_0 / si::kilograms * si::cubic_meters) * theta_std::exner(p_0); // [K kg / (m^2 s)]; -1 because negative gradient of upward flux means inflow
 
       }
 
-      void update_surf_flux_lat(blitz::Array<real_t, n_dims> surf_flux_lat,
-                                       blitz::Array<real_t, n_dims> rt_ground,   
-                                       blitz::Array<real_t, n_dims> U_ground,   
-                                       const real_t &U_ground_z,
-                                       const int &timestep, const real_t &dt, const real_t &dx, const real_t &dy) override
-      {
-        surf_flux_lat = 0;
-      }
-
       // one function for updating u or v
       // the n_dims arrays have vertical extent of 1 - ground calculations only in here
-      void update_surf_flux_uv(blitz::Array<real_t, n_dims>  surf_flux_uv, // output array
-                               blitz::Array<real_t, n_dims>  uv_ground,    // value of u or v on the ground
-                               blitz::Array<real_t, n_dims>  U_ground,     // magnitude of horizontal ground wind
-                               const real_t &U_ground_z,                   // altituted at which U_ground is diagnosed
-                               const int &timestep, const real_t &dt, const real_t &dx, const real_t &dy, const real_t &uv_mean) override
+      void update_surf_flux_uv  (blitz::Array<real_t, n_dims> surf_flux_uv,
+                                 const blitz::Array<real_t, n_dims> &uv_ground,   
+                                 const blitz::Array<real_t, n_dims> &U_ground,   
+                                 const int timestep, const real_t dt, 
+                                 const real_t dx, const real_t dy, const real_t U_ground_z, 
+                                 const real_t uv_mean = 0) override
       {
         surf_flux_uv = - real_t(0.1) * U_ground * (uv_ground + uv_mean) * -1 * (this->rhod_0 / si::kilograms * si::cubic_meters); // [ kg m/s / (m^2 s) ]
       }