diff --git a/CMakeLists.txt b/CMakeLists.txt
index df461be3016ff213e48b8ba55601d301bb048d0e..482b97e3f4e3150ff79db67e6495460b530df1c8 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -95,18 +95,16 @@ if(USE_CXX)
srcCXX/model_base.cpp
srcCXX/sfx_esm.cpp
srcCXX/sfx_sheba.cpp
- srcCXX/sfx_compute_sheba.cpp
srcCXX/sfx_call_class_func.cpp
)
set(HEADERS_CXX
includeCU/sfx_surface.cuh
includeCU/sfx_math.cuh
- includeCU/sfx_esm_compute_subfunc.cuh
+ includeCU/sfx_model_compute_subfunc.cuh
includeCXX/model_base.h
includeCXX/sfx_esm.h
includeCXX/sfx_sheba.h
- includeCXX/sfx_compute_sheba.h
includeCXX/sfx_call_class_func.h
)
endif(USE_CXX)
@@ -114,10 +112,10 @@ endif(USE_CXX)
if(INCLUDE_CUDA)
set(SOURCES_CU
srcCU/sfx_esm.cu
- srcCU/sfx_compute_sheba.cu
+ srcCU/sfx_sheba.cu
)
set(HEADERS_CU
- includeCU/sfx_compute_sheba.cuh
+
)
endif(INCLUDE_CUDA)
diff --git a/includeC/sfx_data.h b/includeC/sfx_data.h
index d1bf113c04df7116ac93176c39702f86037fc048..dba42c80adea9bbcfc6fe5292abf65f153396e01 100644
--- a/includeC/sfx_data.h
+++ b/includeC/sfx_data.h
@@ -60,6 +60,7 @@ extern "C" {
int surface_land;
int surface_lake;
+ float kappa;
float gamma_c;
float Re_visc_min;
float h_charnock;
diff --git a/includeCU/sfx_model_compute_subfunc.cuh b/includeCU/sfx_model_compute_subfunc.cuh
new file mode 100644
index 0000000000000000000000000000000000000000..655555fd8f64ada84b398181edb3d6233548dc57
--- /dev/null
+++ b/includeCU/sfx_model_compute_subfunc.cuh
@@ -0,0 +1,272 @@
+#pragma once
+#include "../includeC/sfx_data.h"
+#include "../includeCU/sfx_math.cuh"
+
+template<typename T, class sfx_param>
+FUCNTION_DECLARATION_SPECIFIER void get_convection_lim(T &zeta_lim, T &Rib_lim, T &f_m_lim, T &f_h_lim,
+ const T h0_m, const T h0_t, const T B,
+ const sfx_param& param)
+{
+ T psi_m, psi_h, f_m, f_h, c;
+
+ c = pow(param.Pr_t_inf_inv / param.Pr_t_0_inv, 4);
+ zeta_lim = (2.0 * param.alpha_h - c * param.alpha_m -
+ sqrt( (c * param.alpha_m)*(c * param.alpha_m) + 4.0 * c * param.alpha_h * (param.alpha_h - param.alpha_m))) / (2.0 * param.alpha_h*param.alpha_h);
+
+ f_m_lim = pow(1.0 - param.alpha_m * zeta_lim, 0.25);
+ f_h_lim = sqrt(1.0 - param.alpha_h * zeta_lim);
+
+ f_m = zeta_lim / h0_m;
+ f_h = zeta_lim / h0_t;
+ if (fabs(B) < 1.0e-10) f_h = f_m;
+
+ f_m = pow(1.0 - param.alpha_m * f_m, 0.25);
+ f_h = sqrt(1.0 - param.alpha_h_fix * f_h);
+
+ psi_m = 2.0 * (atan(f_m_lim) - atan(f_m)) + log((f_m_lim - 1.0) * (f_m + 1.0)/((f_m_lim + 1.0) * (f_m - 1.0)));
+ psi_h = log((f_h_lim - 1.0) * (f_h + 1.0)/((f_h_lim + 1.0) * (f_h - 1.0))) / param.Pr_t_0_inv;
+
+ Rib_lim = zeta_lim * psi_h / (psi_m * psi_m);
+}
+
+template<typename T, class sfx_param>
+FUCNTION_DECLARATION_SPECIFIER void get_psi_stable(T &psi_m, T &psi_h, T &zeta,
+ const T Rib, const T h0_m, const T h0_t, const T B,
+ const sfx_param& param)
+{
+ T Rib_coeff, psi0_m, psi0_h, phi, c;
+
+ psi0_m = log(h0_m);
+ psi0_h = B / psi0_m;
+
+ Rib_coeff = param.beta_m * Rib;
+ c = (psi0_h + 1.0) / param.Pr_t_0_inv - 2.0 * Rib_coeff;
+ zeta = psi0_m * (sqrt(c*c + 4.0 * Rib_coeff * (1.0 - Rib_coeff)) - c) / (2.0 * param.beta_m * (1.0 - Rib_coeff));
+
+ phi = param.beta_m * zeta;
+
+ psi_m = psi0_m + phi;
+ psi_h = (psi0_m + B) / param.Pr_t_0_inv + phi;
+}
+
+template<typename T, class sfx_param>
+FUCNTION_DECLARATION_SPECIFIER void get_psi_convection(T &psi_m, T &psi_h, T &zeta,
+ const T Rib, const T h0_m, const T h0_t, const T B,
+ const T zeta_conv_lim, const T f_m_conv_lim, const T f_h_conv_lim,
+ const sfx_param& param,
+ const int maxiters_convection)
+{
+ T zeta0_m, zeta0_h, f0_m, f0_h, p_m, p_h, a_m, a_h, c_lim, f;
+
+ p_m = 2.0 * atan(f_m_conv_lim) + log((f_m_conv_lim - 1.0) / (f_m_conv_lim + 1.0));
+ p_h = log((f_h_conv_lim - 1.0) / (f_h_conv_lim + 1.0));
+
+ zeta = zeta_conv_lim;
+
+ for (int i = 1; i <= maxiters_convection + 1; i++)
+ {
+ zeta0_m = zeta / h0_m;
+ zeta0_h = zeta / h0_t;
+ if (fabs(B) < 1.0e-10)
+ zeta0_h = zeta0_m;
+
+ f0_m = pow(1.0 - param.alpha_m * zeta0_m, 0.25);
+ f0_h = sqrt(1.0 - param.alpha_h_fix * zeta0_h);
+
+ a_m = -2.0*atan(f0_m) + log((f0_m + 1.0)/(f0_m - 1.0));
+ a_h = log((f0_h + 1.0)/(f0_h - 1.0));
+
+ c_lim = pow(zeta_conv_lim / zeta, 1.0 / 3.0);
+ f = 3.0 * (1.0 - c_lim);
+
+ psi_m = f / f_m_conv_lim + p_m + a_m;
+ psi_h = (f / f_h_conv_lim + p_h + a_h) / param.Pr_t_0_inv;
+
+ if (i == maxiters_convection + 1)
+ break;
+
+ zeta = Rib * psi_m * psi_m / psi_h;
+ }
+}
+
+template<typename T, class sfx_param>
+FUCNTION_DECLARATION_SPECIFIER void get_psi_neutral(T &psi_m, T &psi_h, T &zeta,
+ const T h0_m, const T h0_t, const T B,
+ const sfx_param& param)
+{
+ zeta = 0.0;
+ psi_m = log(h0_m);
+ psi_h = log(h0_t) / param.Pr_t_0_inv;
+ if (fabs(B) < 1.0e-10)
+ psi_h = psi_m / param.Pr_t_0_inv;
+}
+
+template<typename T, class sfx_param>
+FUCNTION_DECLARATION_SPECIFIER void get_psi_semi_convection(T &psi_m, T &psi_h, T &zeta,
+ const T Rib, const T h0_m, const T h0_t, const T B,
+ const sfx_param& param,
+ const int maxiters_convection)
+{
+ T zeta0_m, zeta0_h, f0_m, f0_h, f_m, f_h;
+
+ psi_m = log(h0_m);
+ psi_h = log(h0_t);
+
+ if (fabs(B) < 1.0e-10)
+ psi_h = psi_m;
+
+ zeta = Rib * param.Pr_t_0_inv * psi_m * psi_m / psi_h;
+
+ for (int i = 1; i <= maxiters_convection + 1; i++)
+ {
+ zeta0_m = zeta / h0_m;
+ zeta0_h = zeta / h0_t;
+ if (fabs(B) < 1.0e-10)
+ zeta0_h = zeta0_m;
+
+ f_m = pow(1.0 - param.alpha_m * zeta, 0.25e0);
+ f_h = sqrt(1.0 - param.alpha_h_fix * zeta);
+
+ f0_m = pow(1.0 - param.alpha_m * zeta0_m, 0.25e0);
+ f0_h = sqrt(1.0 - param.alpha_h_fix * zeta0_h);
+
+ f0_m = sfx_math::max(f0_m, T(1.000001e0));
+ f0_h = sfx_math::max(f0_h, T(1.000001e0));
+
+ psi_m = log((f_m - 1.0e0)*(f0_m + 1.0e0)/((f_m + 1.0e0)*(f0_m - 1.0e0))) + 2.0e0*(atan(f_m) - atan(f0_m));
+ psi_h = log((f_h - 1.0e0)*(f0_h + 1.0e0)/((f_h + 1.0e0)*(f0_h - 1.0e0))) / param.Pr_t_0_inv;
+
+ if (i == maxiters_convection + 1)
+ break;
+
+ zeta = Rib * psi_m * psi_m / psi_h;
+ }
+}
+
+template<typename T, class sfx_param>
+FUCNTION_DECLARATION_SPECIFIER void get_psi(T &psi_m, T &psi_h,
+ const T zeta,
+ const sfx_param& param)
+{
+ T x_m, x_h;
+ T q_m, q_h;
+
+ if (zeta >= 0.0)
+ {
+ q_m = pow((1.0 - param.b_m) / param.b_m, 1.0 / 3.0);
+ q_h = sqrt(param.c_h * param.c_h - 4.0);
+
+ x_m = pow(1.0 + zeta, 1.0 / 3.0);
+ x_h = zeta;
+
+ psi_m = -3.0 * (param.a_m / param.b_m) * (x_m - 1.0) + 0.5 * (param.a_m / param.b_m) * q_m * (2.0 * log((x_m + q_m) / (1.0 + q_m)) - log((x_m * x_m - x_m * q_m + q_m * q_m) / (1.0 - q_m + q_m * q_m)) + 2.0 * sqrt(3.0) * (atan((2.0 * x_m - q_m) / (sqrt(3.0) * q_m)) - atan((2.0 - q_m) / (sqrt(3.0) * q_m))));
+
+ psi_h = -0.5 * param.b_h * log(1.0 + param.c_h * x_h + x_h * x_h) + ((-param.a_h / q_h) + ((param.b_h * param.c_h) / (2.0 * q_h))) * (log((2.0 * x_h + param.c_h - q_h) / (2.0 * x_h + param.c_h + q_h)) - log((param.c_h - q_h) / (param.c_h + q_h)));
+ }
+ else
+ {
+ x_m = pow(1.0 - param.alpha_m * zeta, 0.25);
+ x_h = pow(1.0 - param.alpha_h * zeta, 0.25);
+
+ psi_m = (4.0 * atan(1.0) / 2.0) + 2.0 * log(0.5 * (1.0 + x_m)) + log(0.5 * (1.0 + x_m * x_m)) - 2.0 * atan(x_m);
+ psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h));
+ }
+}
+
+template<typename T, class sfx_param>
+FUCNTION_DECLARATION_SPECIFIER void get_psi_mh(T &psi_m, T &psi_h,
+ const T zeta_m, const T zeta_h,
+ const sfx_param& param)
+{
+ T x_m, x_h;
+ T q_m, q_h;
+
+ if (zeta_m >= 0.0)
+ {
+ q_m = pow((1.0 - param.b_m) / param.b_m, 1.0 / 3.0);
+ x_m = pow(1.0 + zeta_m, 1.0 / 3.0);
+
+ psi_m = -3.0 * (param.a_m / param.b_m) * (x_m - 1.0) + 0.5 * (param.a_m / param.b_m) * q_m * (2.0 * log((x_m + q_m) / (1.0 + q_m)) - log((x_m * x_m - x_m * q_m + q_m * q_m) / (1.0 - q_m + q_m * q_m)) + 2.0 * sqrt(3.0) * (atan((2.0 * x_m - q_m) / (sqrt(3.0) * q_m)) - atan((2.0 - q_m) / (sqrt(3.0) * q_m))));
+ }
+ else
+ {
+ x_m = pow(1.0 - param.alpha_m * zeta_m, 0.25);
+ psi_m = (4.0 * atan(1.0) / 2.0) + 2.0 * log(0.5 * (1.0 + x_m)) + log(0.5 * (1.0 + x_m * x_m)) - 2.0 * atan(x_m);
+ }
+
+ if (zeta_h >= 0.0)
+ {
+ q_h = sqrt(param.c_h * param.c_h - 4.0);
+ x_h = zeta_h;
+
+ psi_h = -0.5 * param.b_h * log(1.0 + param.c_h * x_h + x_h * x_h) + ((-param.a_h / q_h) + ((param.b_h * param.c_h) / (2.0 * q_h))) * (log((2.0 * x_h + param.c_h - q_h) / (2.0 * x_h + param.c_h + q_h)) - log((param.c_h - q_h) / (param.c_h + q_h)));
+ }
+ else
+ {
+ x_h = pow(1.0 - param.alpha_h * zeta_h, 0.25);
+ psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h));
+ }
+}
+
+
+template<typename T, class sfx_param>
+FUCNTION_DECLARATION_SPECIFIER void get_dynamic_scales(T &Udyn, T &Tdyn, T &Qdyn, T& zeta,
+ const T U, const T Tsemi, const T dT, const T dQ,
+ const T z, const T z0_m, const T z0_t,
+ const T beta,
+ const sfx_param& param,
+ const int maxiters)
+{
+ const T gamma = T(0.61);
+ T psi_m, psi_h, psi0_m, psi0_h, Linv;
+
+ Udyn = param.kappa * U / log(z / z0_m);
+ Tdyn = param.kappa * dT * param.Pr_t_0_inv / log(z / z0_t);
+ Qdyn = param.kappa * dQ * param.Pr_t_0_inv / log(z / z0_t);
+ zeta = 0.0;
+
+ // --- no wind
+ if (Udyn < 1e-5)
+ return;
+
+ Linv = param.kappa * beta * (Tdyn + gamma * Qdyn * Tsemi) / (Udyn * Udyn);
+ zeta = z * Linv;
+
+ // --- near neutral case
+ if (Linv < 1e-5)
+ return;
+
+ for (int i = 0; i < maxiters; i++)
+ {
+ get_psi(psi_m, psi_h, zeta, param);
+ get_psi_mh(psi0_m, psi0_h, z0_m * Linv, z0_t * Linv,
+ param);
+
+ Udyn = param.kappa * U / (log(z / z0_m) - (psi_m - psi0_m));
+ Tdyn = param.kappa * dT * param.Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h));
+ Qdyn = param.kappa * dQ * param.Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h));
+
+ if (Udyn < 1e-5)
+ break;
+
+ Linv = param.kappa * beta * (Tdyn + gamma * Qdyn * Tsemi) / (Udyn * Udyn);
+ zeta = z * Linv;
+ }
+}
+
+template<typename T, class sfx_param>
+FUCNTION_DECLARATION_SPECIFIER void get_phi(T &phi_m, T &phi_h,
+ const T zeta,
+ const sfx_param& param)
+{
+ if (zeta >= 0.0)
+ {
+ phi_m = 1.0 + (param.a_m * zeta * pow(1.0 + zeta, 1.0 / 3.0) ) / (1.0 + param.b_m * zeta);
+ phi_h = 1.0 + (param.a_h * zeta + param.b_h * zeta * zeta) / (1.0 + param.c_h * zeta + zeta * zeta);
+ }
+ else
+ {
+ phi_m = pow(1.0 - param.alpha_m * zeta, -0.25);
+ phi_h = pow(1.0 - param.alpha_h * zeta, -0.5);
+ }
+}
diff --git a/includeCU/sfx_surface.cuh b/includeCU/sfx_surface.cuh
index 2344f41c42c0c6de0b8805863a1972833b59f954..60a0aeba91e6542a9a5766634071f455d58c0513 100644
--- a/includeCU/sfx_surface.cuh
+++ b/includeCU/sfx_surface.cuh
@@ -24,29 +24,28 @@ FUCNTION_DECLARATION_SPECIFIER void get_thermal_roughness(T &z0_t, T &B,
template<typename T>
FUCNTION_DECLARATION_SPECIFIER void get_charnock_roughness(T &z0_m, T &u_dyn0,
const T h, const T U,
- const sfx_esm_param& model_param,
const sfx_surface_param& surface_param,
- const sfx_esm_numericsTypeC& numerics)
+ const int maxiters)
{
T Uc, a, b, c, c_min, f;
Uc = U;
a = 0.0;
b = 25.0;
- c_min = log(surface_param.h_charnock) / model_param.kappa;
+ c_min = log(surface_param.h_charnock) / surface_param.kappa;
- for (int i = 0; i < numerics.maxiters_charnock; i++)
+ for (int i = 0; i < maxiters; i++)
{
f = surface_param.c1_charnock - 2.0 * log(Uc);
- for (int j = 0; j < numerics.maxiters_charnock; j++)
+ for (int j = 0; j < maxiters; j++)
{
- c = (f + 2.0 * log(b)) / model_param.kappa;
+ c = (f + 2.0 * log(b)) / surface_param.kappa;
if (U <= 8.0e0)
- a = log(1.0 + surface_param.c2_charnock * ( pow(b / Uc, 3) ) ) / model_param.kappa;
+ a = log(1.0 + surface_param.c2_charnock * ( pow(b / Uc, 3) ) ) / surface_param.kappa;
c = sfx_math::max(c - a, c_min);
b = c;
}
- z0_m = surface_param.h_charnock * exp(-c * model_param.kappa);
+ z0_m = surface_param.h_charnock * exp(-c * surface_param.kappa);
z0_m = sfx_math::max(z0_m, T(0.000015e0));
Uc = U * log(surface_param.h_charnock / z0_m) / log(h / z0_m);
}
diff --git a/includeCXX/sfx_sheba.h b/includeCXX/sfx_sheba.h
index 676980f1fe273a040481b2a3dc414ab84c7daba8..ae9e7fff1d2503c96aa2885bb01ec2863945af25 100644
--- a/includeCXX/sfx_sheba.h
+++ b/includeCXX/sfx_sheba.h
@@ -1,48 +1,88 @@
#pragma once
#include "sfx_template_parameters.h"
-#include <cstddef>
+#include "../includeCXX/model_base.h"
+#include "../includeC/sfx_data.h"
template<typename T, MemType memIn, MemType memOut, MemType RunMem >
-class FluxSheba
+class FluxShebaBase : public ModelBase<T, memIn, memOut, RunMem>
{
-private:
- T *U, *dT, *Tsemi, *dQ, *h, *in_z0_m;
- T *zeta, *Rib, *Re, *B, *z0_m, *z0_t, *Rib_conv_lim, *Cm, *Ct, *Km, *Pr_t_inv;
+public:
+ using ModelBase<T, memIn, memOut, RunMem>::res_sfx;
+ using ModelBase<T, memIn, memOut, RunMem>::sfx;
+ using ModelBase<T, memIn, memOut, RunMem>::meteo;
+ using ModelBase<T, memIn, memOut, RunMem>::grid_size;
+ using ModelBase<T, memIn, memOut, RunMem>::ifAllocated;
+ using ModelBase<T, memIn, memOut, RunMem>::allocated_size;
- T kappa, Pr_t_0_inv,
- alpha_m, alpha_h,
- a_m, a_h,
- b_m, b_h,
- c_h,
- Re_rough_min,
- B1_rough, B2_rough,
- B_max_land, B_max_ocean, B_max_lake,
- gamma_c,
- Re_visc_min,
- Pr_m, nu_air, g;
+ sfx_surface_param surface_param;
+ sfx_phys_constants phys_constants;
+ sfx_sheba_param model_param;
+ sfx_sheba_numericsTypeC numerics;
- int grid_size;
- bool ifAllocated;
- size_t allocated_size;
-public:
- FluxSheba();
- void set_params(const int grid_size, const T kappa, const T Pr_t_0_inv,
- const T alpha_m, const T alpha_h,
- const float a_m, const float a_h,
- const float b_m, const float b_h,
- const float c_h,
- const T Re_rough_min,
- const T B1_rough, const T B2_rough,
- const T B_max_land, const T B_max_ocean, const T B_max_lake,
- const T gamma_c, const T Re_visc_min,
- const T Pr_m, const T nu_air, const T g);
- ~FluxSheba();
+ FluxShebaBase(sfxDataVecTypeC* sfx,
+ meteoDataVecTypeC* meteo,
+ const sfx_sheba_param model_param,
+ const sfx_surface_param surface_param,
+ const sfx_sheba_numericsTypeC numerics,
+ const sfx_phys_constants phys_constants,
+ const int grid_size);
+ ~FluxShebaBase();
+};
+
+template<typename T, MemType memIn, MemType memOut, MemType RunMem >
+class FluxSheba : public FluxShebaBase<T, memIn, memOut, RunMem>
+{};
- void compute_flux_sheba(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_, T *Rib_conv_lim_, T *Cm_, T *Ct_, T *Km_, T *Pr_t_inv_,
- T *U_, T *dT_, T *Tsemi_, T *dQ_, T *h_, T *in_z0_m_,
- const int maxiters_charnock);
+template<typename T, MemType memIn, MemType memOut >
+class FluxSheba<T, memIn, memOut, MemType::CPU> : public FluxShebaBase<T, memIn, memOut, MemType::CPU>
+{
+ using FluxShebaBase<T, memIn, memOut, MemType::CPU>::res_sfx;
+ using FluxShebaBase<T, memIn, memOut, MemType::CPU>::sfx;
+ using FluxShebaBase<T, memIn, memOut, MemType::CPU>::meteo;
+ using FluxShebaBase<T, memIn, memOut, MemType::CPU>::surface_param;
+ using FluxShebaBase<T, memIn, memOut, MemType::CPU>::phys_constants;
+ using FluxShebaBase<T, memIn, memOut, MemType::CPU>::grid_size;
+ using FluxShebaBase<T, memIn, memOut, MemType::CPU>::ifAllocated;
+ using FluxShebaBase<T, memIn, memOut, MemType::CPU>::allocated_size;
+ using FluxShebaBase<T, memIn, memOut, MemType::CPU>::model_param;
+ using FluxShebaBase<T, memIn, memOut, MemType::CPU>::numerics;
+public:
+ FluxSheba(sfxDataVecTypeC* sfx,
+ meteoDataVecTypeC* meteo,
+ const sfx_sheba_param model_param,
+ const sfx_surface_param surface_param,
+ const sfx_sheba_numericsTypeC numerics,
+ const sfx_phys_constants phys_constants,
+ const int grid_size) : FluxShebaBase<T, memIn, memOut, MemType::CPU>(sfx, meteo, model_param,
+ surface_param, numerics, phys_constants, grid_size) {}
+ ~FluxSheba() = default;
+ void compute_flux();
+};
-private:
- void set_data( T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_, T *Rib_conv_lim_, T *Cm_, T *Ct_, T *Km_, T *Pr_t_inv_,
- T *U_, T *dT_, T *Tsemi_, T *dQ_, T *h_, T *in_z0_m_);
-};
\ No newline at end of file
+#ifdef INCLUDE_CUDA
+template<typename T, MemType memIn, MemType memOut >
+class FluxSheba<T, memIn, memOut, MemType::GPU> : public FluxShebaBase<T, memIn, memOut, MemType::GPU>
+{
+ using FluxShebaBase<T, memIn, memOut, MemType::GPU>::res_sfx;
+ using FluxShebaBase<T, memIn, memOut, MemType::GPU>::sfx;
+ using FluxShebaBase<T, memIn, memOut, MemType::GPU>::meteo;
+ using FluxShebaBase<T, memIn, memOut, MemType::GPU>::surface_param;
+ using FluxShebaBase<T, memIn, memOut, MemType::GPU>::phys_constants;
+ using FluxShebaBase<T, memIn, memOut, MemType::GPU>::grid_size;
+ using FluxShebaBase<T, memIn, memOut, MemType::GPU>::ifAllocated;
+ using FluxShebaBase<T, memIn, memOut, MemType::GPU>::allocated_size;
+ using FluxShebaBase<T, memIn, memOut, MemType::GPU>::model_param;
+ using FluxShebaBase<T, memIn, memOut, MemType::GPU>::numerics;
+public:
+ FluxSheba(sfxDataVecTypeC* sfx,
+ meteoDataVecTypeC* meteo,
+ const sfx_sheba_param model_param,
+ const sfx_surface_param surface_param,
+ const sfx_sheba_numericsTypeC numerics,
+ const sfx_phys_constants phys_constants,
+ const int grid_size) : FluxShebaBase<T, memIn, memOut, MemType::GPU>(sfx, meteo, model_param,
+ surface_param, numerics, phys_constants, grid_size) {}
+ ~FluxSheba() = default;
+ void compute_flux();
+};
+#endif
\ No newline at end of file
diff --git a/srcCU/sfx_esm.cu b/srcCU/sfx_esm.cu
index c2d7aef58647901e46f277a699ce3dea8573a478..49595c249c59ea23e40151477c95c772edacf37c 100644
--- a/srcCU/sfx_esm.cu
+++ b/srcCU/sfx_esm.cu
@@ -2,7 +2,7 @@
#include <cuda_runtime_api.h>
#include "../includeCXX/sfx_esm.h"
-#include "../includeCU/sfx_esm_compute_subfunc.cuh"
+#include "../includeCU/sfx_model_compute_subfunc.cuh"
#include "../includeCU/sfx_surface.cuh"
#include "../includeCU/sfx_memory_processing.cuh"
@@ -46,7 +46,7 @@ __global__ void sfx_kernel::compute_flux(sfxDataVecTypeC sfx,
if (surface_type == surface_param.surface_ocean)
{
- get_charnock_roughness(z0_m, u_dyn0, h, U, model_param, surface_param, numerics);
+ get_charnock_roughness(z0_m, u_dyn0, h, U, surface_param, numerics.maxiters_charnock);
h0_m = h / z0_m;
}
if (surface_type == surface_param.surface_land)
@@ -77,7 +77,7 @@ __global__ void sfx_kernel::compute_flux(sfxDataVecTypeC sfx,
}
else if (Rib < Rib_conv_lim)
{
- get_psi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, model_param, numerics);
+ get_psi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, model_param, numerics.maxiters_convection);
fval = pow(zeta_conv_lim / zeta, 1.0/3.0);
phi_m = fval / f_m_conv_lim;
@@ -92,7 +92,7 @@ __global__ void sfx_kernel::compute_flux(sfxDataVecTypeC sfx,
}
else
{
- get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model_param, numerics);
+ get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model_param, numerics.maxiters_convection);
phi_m = pow(1.0 - model_param.alpha_m * zeta, -0.25);
phi_h = 1.0 / (model_param.Pr_t_0_inv * sqrt(1.0 - model_param.alpha_h_fix * zeta));
diff --git a/srcCU/sfx_sheba.cu b/srcCU/sfx_sheba.cu
new file mode 100644
index 0000000000000000000000000000000000000000..c5c33f806eda335662b1c3d1e399d15bbe7ffd53
--- /dev/null
+++ b/srcCU/sfx_sheba.cu
@@ -0,0 +1,144 @@
+#include <cuda.h>
+#include <cuda_runtime_api.h>
+
+#include "../includeCXX/sfx_sheba.h"
+#include "../includeCU/sfx_model_compute_subfunc.cuh"
+#include "../includeCU/sfx_surface.cuh"
+#include "../includeCU/sfx_memory_processing.cuh"
+
+namespace sfx_kernel
+{
+ template<typename T>
+ __global__ void compute_flux(sfxDataVecTypeC sfx,
+ meteoDataVecTypeC meteo,
+ const sfx_sheba_param model_param,
+ const sfx_surface_param surface_param,
+ const sfx_sheba_numericsTypeC numerics,
+ const sfx_phys_constants phys_constants,
+ const int grid_size);
+}
+
+template<typename T>
+__global__ void sfx_kernel::compute_flux(sfxDataVecTypeC sfx,
+ meteoDataVecTypeC meteo,
+ const sfx_sheba_param model_param,
+ const sfx_surface_param surface_param,
+ const sfx_sheba_numericsTypeC numerics,
+ const sfx_phys_constants phys_constants,
+ const int grid_size)
+{
+ const int index = blockIdx.x * blockDim.x + threadIdx.x;
+ T h, U, dT, Tsemi, dQ, z0_m;
+ T z0_t, B, h0_m, h0_t, u_dyn0, Re,
+ zeta, Rib, Udyn, Tdyn, Qdyn, phi_m, phi_h,
+ Km, Pr_t_inv, Cm, Ct;
+ int surface_type;
+
+ if(index < grid_size)
+ {
+ U = meteo.U[index];
+ Tsemi = meteo.Tsemi[index];
+ dT = meteo.dT[index];
+ dQ = meteo.dQ[index];
+ h = meteo.h[index];
+ z0_m = meteo.z0_m[index];
+
+ surface_type = z0_m < 0.0 ? surface_param.surface_ocean : surface_param.surface_land;
+
+ if (surface_type == surface_param.surface_ocean)
+ {
+ get_charnock_roughness(z0_m, u_dyn0, h, U, surface_param, numerics.maxiters_charnock);
+ h0_m = h / z0_m;
+ }
+ if (surface_type == surface_param.surface_land)
+ {
+ h0_m = h / z0_m;
+ u_dyn0 = U * model_param.kappa / log(h0_m);
+ }
+
+ Re = u_dyn0 * z0_m / phys_constants.nu_air;
+ get_thermal_roughness(z0_t, B, z0_m, Re, surface_param, surface_type);
+
+ // --- define relative height [thermal]
+ h0_t = h / z0_t;
+
+ // --- define Ri-bulk
+ Rib = (phys_constants.g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / (U*U);
+
+ // --- get the fluxes
+ // ----------------------------------------------------------------------------
+ get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta,
+ U, Tsemi, dT, dQ, h, z0_m, z0_t, (phys_constants.g / Tsemi), model_param, 10);
+ // ----------------------------------------------------------------------------
+
+ get_phi(phi_m, phi_h, zeta, model_param);
+ // ----------------------------------------------------------------------------
+
+ // --- define transfer coeff. (momentum) & (heat)
+ Cm = 0.0;
+ if (U > 0.0)
+ Cm = Udyn / U;
+ Ct = 0.0;
+ if (fabs(dT) > 0.0)
+ Ct = Tdyn / dT;
+
+ // --- define eddy viscosity & inverse Prandtl number
+ Km = model_param.kappa * Cm * U * h / phi_m;
+ Pr_t_inv = phi_m / phi_h;
+
+ sfx.zeta[index] = zeta;
+ sfx.Rib[index] = Rib;
+ sfx.Re[index] = Re;
+ sfx.B[index] = B;
+ sfx.z0_m[index] = z0_m;
+ sfx.z0_t[index] = z0_t;
+ sfx.Rib_conv_lim[index] = T(0.0);
+ sfx.Cm[index] = Cm;
+ sfx.Ct[index] = Ct;
+ sfx.Km[index] = Km;
+ sfx.Pr_t_inv[index] = Pr_t_inv;
+ }
+}
+
+template<typename T, MemType memIn, MemType memOut >
+void FluxSheba<T, memIn, memOut, MemType::GPU>::compute_flux()
+{
+ const int BlockCount = int(ceil(float(grid_size) / 1024.0));
+ dim3 cuBlock = dim3(1024, 1, 1);
+ dim3 cuGrid = dim3(BlockCount, 1, 1);
+
+ sfx_kernel::compute_flux<T><<<cuGrid, cuBlock>>>(sfx, meteo, model_param,
+ surface_param, numerics, phys_constants, grid_size);
+
+ if(MemType::GPU != memOut)
+ {
+ const size_t new_size = grid_size * sizeof(T);
+ memproc::memcopy<memOut, MemType::GPU>((void*&)res_sfx->zeta, (void*&)sfx.zeta, new_size);
+ memproc::memcopy<memOut, MemType::GPU>((void*&)res_sfx->Rib, (void*&)sfx.Rib, new_size);
+ memproc::memcopy<memOut, MemType::GPU>((void*&)res_sfx->Re, (void*&)sfx.Re, new_size);
+ memproc::memcopy<memOut, MemType::GPU>((void*&)res_sfx->B, (void*&)sfx.B, new_size);
+ memproc::memcopy<memOut, MemType::GPU>((void*&)res_sfx->z0_m, (void*&)sfx.z0_m, new_size);
+ memproc::memcopy<memOut, MemType::GPU>((void*&)res_sfx->z0_t, (void*&)sfx.z0_t, new_size);
+ memproc::memcopy<memOut, MemType::GPU>((void*&)res_sfx->Rib_conv_lim, (void*&)sfx.Rib_conv_lim, new_size);
+ memproc::memcopy<memOut, MemType::GPU>((void*&)res_sfx->Cm, (void*&)sfx.Cm, new_size);
+ memproc::memcopy<memOut, MemType::GPU>((void*&)res_sfx->Ct, (void*&)sfx.Ct, new_size);
+ memproc::memcopy<memOut, MemType::GPU>((void*&)res_sfx->Km, (void*&)sfx.Km, new_size);
+ memproc::memcopy<memOut, MemType::GPU>((void*&)res_sfx->Pr_t_inv, (void*&)sfx.Pr_t_inv, new_size);
+ }
+}
+
+template class FluxShebaBase<float, MemType::GPU, MemType::GPU, MemType::GPU>;
+template class FluxShebaBase<float, MemType::GPU, MemType::GPU, MemType::CPU>;
+template class FluxShebaBase<float, MemType::GPU, MemType::CPU, MemType::GPU>;
+template class FluxShebaBase<float, MemType::CPU, MemType::GPU, MemType::GPU>;
+template class FluxShebaBase<float, MemType::CPU, MemType::CPU, MemType::GPU>;
+template class FluxShebaBase<float, MemType::CPU, MemType::GPU, MemType::CPU>;
+template class FluxShebaBase<float, MemType::GPU, MemType::CPU, MemType::CPU>;
+
+template class FluxSheba<float, MemType::GPU, MemType::GPU, MemType::GPU>;
+template class FluxSheba<float, MemType::GPU, MemType::GPU, MemType::CPU>;
+template class FluxSheba<float, MemType::GPU, MemType::CPU, MemType::GPU>;
+template class FluxSheba<float, MemType::CPU, MemType::GPU, MemType::GPU>;
+template class FluxSheba<float, MemType::CPU, MemType::CPU, MemType::GPU>;
+template class FluxSheba<float, MemType::CPU, MemType::GPU, MemType::CPU>;
+template class FluxSheba<float, MemType::GPU, MemType::CPU, MemType::CPU>;
\ No newline at end of file
diff --git a/srcCXX/sfx_esm.cpp b/srcCXX/sfx_esm.cpp
index ad77a702428a957d5c18d27f3e227d2f35e0b72b..18e1c7fb2635000ff1f912d8457fb5afbb1c7178 100644
--- a/srcCXX/sfx_esm.cpp
+++ b/srcCXX/sfx_esm.cpp
@@ -8,7 +8,7 @@
#include "../includeCXX/sfx_memory_processing.h"
#include "../includeCU/sfx_surface.cuh"
-#include "../includeCU/sfx_esm_compute_subfunc.cuh"
+#include "../includeCU/sfx_model_compute_subfunc.cuh"
template<typename T, MemType memIn, MemType memOut, MemType RunMem >
FluxEsmBase<T, memIn, memOut, RunMem>::FluxEsmBase(sfxDataVecTypeC* sfx_in,
@@ -49,7 +49,7 @@ void FluxEsm<T, memIn, memOut, MemType::CPU>::compute_flux()
if (surface_type == surface_param.surface_ocean)
{
- get_charnock_roughness(z0_m, u_dyn0, h, U, model_param, surface_param, numerics);
+ get_charnock_roughness(z0_m, u_dyn0, h, U, surface_param, numerics.maxiters_charnock);
h0_m = h / z0_m;
}
if (surface_type == surface_param.surface_land)
@@ -80,7 +80,7 @@ void FluxEsm<T, memIn, memOut, MemType::CPU>::compute_flux()
}
else if (Rib < Rib_conv_lim)
{
- get_psi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, model_param, numerics);
+ get_psi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, model_param, numerics.maxiters_convection);
fval = pow(zeta_conv_lim / zeta, 1.0/3.0);
phi_m = fval / f_m_conv_lim;
@@ -95,7 +95,7 @@ void FluxEsm<T, memIn, memOut, MemType::CPU>::compute_flux()
}
else
{
- get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model_param, numerics);
+ get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model_param, numerics.maxiters_convection);
phi_m = pow(1.0 - model_param.alpha_m * zeta, -0.25);
phi_h = 1.0 / (model_param.Pr_t_0_inv * sqrt(1.0 - model_param.alpha_h_fix * zeta));
diff --git a/srcCXX/sfx_sheba.cpp b/srcCXX/sfx_sheba.cpp
index af52cc2c78fcc2578815b2b9a75e02a96ff04c7e..8417854ff63b0a35d7d2ff7cde63325941fb497e 100644
--- a/srcCXX/sfx_sheba.cpp
+++ b/srcCXX/sfx_sheba.cpp
@@ -1,281 +1,136 @@
#include <iostream>
+#include <cmath>
#include "../includeCXX/sfx_sheba.h"
-#include "../includeCXX/sfx_compute_sheba.h"
#ifdef INCLUDE_CUDA
- #include "../includeCU/sfx_compute_sheba.cuh"
#include "../includeCU/sfx_memory_processing.cuh"
#endif
#include "../includeCXX/sfx_memory_processing.h"
+#include "../includeCU/sfx_surface.cuh"
+#include "../includeCU/sfx_model_compute_subfunc.cuh"
template<typename T, MemType memIn, MemType memOut, MemType RunMem >
-FluxSheba<T, memIn, memOut, RunMem>::FluxSheba()
+FluxShebaBase<T, memIn, memOut, RunMem>::FluxShebaBase(sfxDataVecTypeC* sfx_in,
+ meteoDataVecTypeC* meteo_in,
+ const sfx_sheba_param model_param_in,
+ const sfx_surface_param surface_param_in,
+ const sfx_sheba_numericsTypeC numerics_in,
+ const sfx_phys_constants phys_constants_in,
+ const int grid_size_in) : ModelBase<T, memIn, memOut, RunMem>(sfx_in, meteo_in, grid_size_in)
{
- kappa = 0, Pr_t_0_inv = 0,
- alpha_m = 0, alpha_h = 0,
- a_m = 0, a_h = 0,
- b_m = 0, b_h = 0,
- c_h = 0,
- Re_rough_min = 0,
- B1_rough = 0, B2_rough = 0,
- B_max_land = 0, B_max_ocean = 0, B_max_lake = 0,
- gamma_c = 0,
- Re_visc_min = 0,
- Pr_m = 0, nu_air = 0, g = 0;
-
- grid_size = 0;
-
- ifAllocated = false;
- allocated_size = 0;
+ surface_param = surface_param_in;
+ phys_constants = phys_constants_in;
+ model_param = model_param_in;
+ numerics = numerics_in;
}
template<typename T, MemType memIn, MemType memOut, MemType RunMem >
-void FluxSheba<T, memIn, memOut, RunMem>::set_params(const int grid_size_, const T kappa_, const T Pr_t_0_inv_,
- const T alpha_m_, const T alpha_h_,
- const float a_m_, const float a_h_,
- const float b_m_, const float b_h_,
- const float c_h_,
- const T Re_rough_min_,
- const T B1_rough_, const T B2_rough_,
- const T B_max_land_, const T B_max_ocean_, const T B_max_lake_,
- const T gamma_c_, const T Re_visc_min_,
- const T Pr_m_, const T nu_air_, const T g_)
-{
- kappa = kappa_, Pr_t_0_inv = Pr_t_0_inv_,
- alpha_m = alpha_m_, alpha_h = alpha_h_,
- a_m = a_m_, a_h = a_h_,
- b_m = b_m_, b_h = b_h_,
- c_h = c_h_,
- Re_rough_min = Re_rough_min_, B_max_lake = B_max_lake_,
- B1_rough = B1_rough_, B2_rough = B2_rough_,
- B_max_land = B_max_land_, B_max_ocean = B_max_ocean_, B_max_lake = B_max_lake_,
- gamma_c = gamma_c_,
- Re_visc_min = Re_visc_min_,
- Pr_m = Pr_m_, nu_air = nu_air_, g = g_;
-
- grid_size = grid_size_;
-
- if(RunMem != memIn)
- {
- const size_t new_size = grid_size * sizeof(T);
+FluxShebaBase<T, memIn, memOut, RunMem>::~FluxShebaBase() {}
- allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(U), allocated_size, new_size);
- allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(dT), allocated_size, new_size);
- allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(Tsemi), allocated_size, new_size);
- allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(dQ), allocated_size, new_size);
- allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(h), allocated_size, new_size);
- allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(in_z0_m), allocated_size, new_size);
+template<typename T, MemType memIn, MemType memOut >
+void FluxSheba<T, memIn, memOut, MemType::CPU>::compute_flux()
+{
+ T h, U, dT, Tsemi, dQ, z0_m;
+ T z0_t, B, h0_m, h0_t, u_dyn0, Re,
+ zeta, Rib, Udyn, Tdyn, Qdyn, phi_m, phi_h,
+ Km, Pr_t_inv, Cm, Ct;
+ int surface_type;
- ifAllocated = true;
- }
- if(RunMem != memOut)
+ for (int step = 0; step < grid_size; step++)
{
- const size_t new_size = grid_size * sizeof(T);
- allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(zeta), allocated_size, new_size);
- allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(Rib), allocated_size, new_size);
- allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(Re), allocated_size, new_size);
- allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(Rib_conv_lim), allocated_size, new_size);
- allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(z0_m), allocated_size, new_size);
- allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(z0_t), allocated_size, new_size);
- allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(B), allocated_size, new_size);
- allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(Cm), allocated_size, new_size);
- allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(Ct), allocated_size, new_size);
- allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(Km), allocated_size, new_size);
- allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(Pr_t_inv), allocated_size, new_size);
-
- // T* del = new T[grid_size];
- // for (int i = 0; i < grid_size; i++)
- // del[i] = i;
-
- // memproc::memcopy<RunMem, memOut>(Cm, del, new_size);
-
- // for (int i = 0; i < grid_size; i++)
- // del[i] = 0;
+ U = meteo.U[step];
+ Tsemi = meteo.Tsemi[step];
+ dT = meteo.dT[step];
+ dQ = meteo.dQ[step];
+ h = meteo.h[step];
+ z0_m = meteo.z0_m[step];
- // printf("before\n");
- // for (int i = 0; i < 10; i++)
- // printf("%f\n", del[i]);
+ surface_type = z0_m < 0.0 ? surface_param.surface_ocean : surface_param.surface_land;
- // memproc::memcopy<memOut, RunMem>(del, Cm, new_size);
-
- // for (int i = 0; i < grid_size; i++)
- // printf("%f\n", del[i]);
-
- // delete[] del;
-
- ifAllocated = true;
- }
-}
-
-template<typename T, MemType memIn, MemType memOut, MemType RunMem >
-FluxSheba<T, memIn, memOut, RunMem>::~FluxSheba()
-{
- kappa = 0, Pr_t_0_inv = 0,
- alpha_m = 0, alpha_h = 0,
- a_m = 0, a_h = 0,
- b_m = 0, b_h = 0,
- c_h = 0,
- Re_rough_min = 0,
- B1_rough = 0, B2_rough = 0,
- B_max_land = 0, B_max_ocean = 0, B_max_lake = 0,
- gamma_c = 0,
- Re_visc_min = 0,
- Pr_m = 0, nu_air = 0, g = 0;
-
- grid_size = 0;
-
- if(ifAllocated == true)
- {
- if(RunMem != memIn)
- {
- memproc::dealloc<RunMem>((void*&)U);
- memproc::dealloc<RunMem>((void*&)dT);
- memproc::dealloc<RunMem>((void*&)Tsemi);
- memproc::dealloc<RunMem>((void*&)dQ);
- memproc::dealloc<RunMem>((void*&)h);
+ if (surface_type == surface_param.surface_ocean)
+ {
+ get_charnock_roughness(z0_m, u_dyn0, h, U, surface_param, numerics.maxiters_charnock);
+ h0_m = h / z0_m;
}
- if(RunMem != memOut)
+ if (surface_type == surface_param.surface_land)
{
- memproc::dealloc<RunMem>((void*&)zeta);
- memproc::dealloc<RunMem>((void*&)Rib);
- memproc::dealloc<RunMem>((void*&)Re);
- memproc::dealloc<RunMem>((void*&)Rib_conv_lim);
- memproc::dealloc<RunMem>((void*&)z0_m);
- memproc::dealloc<RunMem>((void*&)z0_t);
- memproc::dealloc<RunMem>((void*&)B);
- memproc::dealloc<RunMem>((void*&)Cm);
- memproc::dealloc<RunMem>((void*&)Ct);
- memproc::dealloc<RunMem>((void*&)Km);
- memproc::dealloc<RunMem>((void*&)Pr_t_inv);
+ h0_m = h / z0_m;
+ u_dyn0 = U * model_param.kappa / log(h0_m);
}
- ifAllocated = false;
- allocated_size = 0;
+ Re = u_dyn0 * z0_m / phys_constants.nu_air;
+ get_thermal_roughness(z0_t, B, z0_m, Re, surface_param, surface_type);
+
+ // --- define relative height [thermal]
+ h0_t = h / z0_t;
+
+ // --- define Ri-bulk
+ Rib = (phys_constants.g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / (U*U);
+
+ // --- get the fluxes
+ // ----------------------------------------------------------------------------
+ get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta,
+ U, Tsemi, dT, dQ, h, z0_m, z0_t, (phys_constants.g / Tsemi), model_param, 10);
+ // ----------------------------------------------------------------------------
+
+ get_phi(phi_m, phi_h, zeta, model_param);
+ // ----------------------------------------------------------------------------
+
+ // --- define transfer coeff. (momentum) & (heat)
+ Cm = 0.0;
+ if (U > 0.0)
+ Cm = Udyn / U;
+ Ct = 0.0;
+ if (fabs(dT) > 0.0)
+ Ct = Tdyn / dT;
+
+ // --- define eddy viscosity & inverse Prandtl number
+ Km = model_param.kappa * Cm * U * h / phi_m;
+ Pr_t_inv = phi_m / phi_h;
+
+ sfx.zeta[step] = zeta;
+ sfx.Rib[step] = Rib;
+ sfx.Re[step] = Re;
+ sfx.B[step] = B;
+ sfx.z0_m[step] = z0_m;
+ sfx.z0_t[step] = z0_t;
+ sfx.Rib_conv_lim[step] = T(0.0);
+ sfx.Cm[step] = Cm;
+ sfx.Ct[step] = Ct;
+ sfx.Km[step] = Km;
+ sfx.Pr_t_inv[step] = Pr_t_inv;
}
-}
-
-template<typename T, MemType memIn, MemType memOut, MemType RunMem >
-void FluxSheba<T, memIn, memOut, RunMem>::set_data(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_, T *Rib_conv_lim_, T *Cm_, T *Ct_, T *Km_, T *Pr_t_inv_,
- T *U_, T *dT_, T *Tsemi_, T *dQ_, T *h_, T *in_z0_m_)
-{
- if(RunMem == memIn)
- {
- U = U_;
- dT = dT_;
- Tsemi = Tsemi_;
- dQ = dQ_;
- h = h_;
- in_z0_m = in_z0_m_;
- }
- else
- {
- const size_t new_size = grid_size * sizeof(T);
-
- memproc::memcopy<RunMem, memIn>(U, U_, new_size);
- memproc::memcopy<RunMem, memIn>(dT, dT_, new_size);
- memproc::memcopy<RunMem, memIn>(Tsemi, Tsemi_, new_size);
- memproc::memcopy<RunMem, memIn>(dQ, dQ_, new_size);
- memproc::memcopy<RunMem, memIn>(h, h_, new_size);
- memproc::memcopy<RunMem, memIn>(in_z0_m, in_z0_m_, new_size);
- }
-
- if(RunMem == memOut)
- {
- zeta = zeta_;
- Rib = Rib_;
- Re = Re_;
- Rib_conv_lim = Rib_conv_lim_;
- z0_m = z0_m_;
- z0_t = z0_t_;
- B = B_;
- Cm = Cm_;
- Ct = Ct_;
- Km = Km_;
- Pr_t_inv = Pr_t_inv_;
- }
-}
-
-template<typename T, MemType memIn, MemType memOut, MemType RunMem >
-void FluxSheba<T, memIn, memOut, RunMem>::compute_flux_sheba(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_, T *Rib_conv_lim_, T *Cm_, T *Ct_, T *Km_, T *Pr_t_inv_,
- T *U_, T *dT_, T *Tsemi_, T *dQ_, T *h_, T *in_z0_m_,
- const int maxiters_charnock)
-{
- set_data(zeta_, Rib_, Re_, B_, z0_m_, z0_t_, Rib_conv_lim_, Cm_, Ct_, Km_, Pr_t_inv_,
- U_, dT_, Tsemi_, dQ_, h_, in_z0_m_);
-
- if(RunMem == MemType::CPU)
- compute_flux_sheba_cpu(zeta, Rib, Re, B, z0_m, z0_t, Rib_conv_lim, Cm, Ct, Km, Pr_t_inv,
- U, dT, Tsemi, dQ, h, in_z0_m,
- kappa, Pr_t_0_inv,
- alpha_m, alpha_h,
- a_m, a_h,
- b_m, b_h,
- c_h,
- Re_rough_min,
- B1_rough, B2_rough,
- B_max_land, B_max_ocean, B_max_lake,
- gamma_c, Re_visc_min,
- Pr_m, nu_air, g,
- maxiters_charnock,
- grid_size);
-#ifdef INCLUDE_CUDA
- else compute_flux_sheba_gpu(zeta, Rib, Re, B, z0_m, z0_t, Rib_conv_lim, Cm, Ct, Km, Pr_t_inv,
- U, dT, Tsemi, dQ, h, in_z0_m,
- kappa, Pr_t_0_inv,
- alpha_m, alpha_h,
- a_m, a_h,
- b_m, b_h,
- c_h,
- Re_rough_min,
- B1_rough, B2_rough,
- B_max_land, B_max_ocean, B_max_lake,
- gamma_c, Re_visc_min,
- Pr_m, nu_air, g,
- maxiters_charnock,
- grid_size);
-#endif
-
- if(RunMem != memOut)
+ if(MemType::CPU != memOut)
{
const size_t new_size = grid_size * sizeof(T);
-
- memproc::memcopy<memOut, RunMem>(zeta_, zeta, new_size);
- memproc::memcopy<memOut, RunMem>(Rib_, Rib, new_size);
- memproc::memcopy<memOut, RunMem>(Re_, Re, new_size);
- memproc::memcopy<memOut, RunMem>(Rib_conv_lim_, Rib_conv_lim, new_size);
- memproc::memcopy<memOut, RunMem>(z0_m_, z0_m, new_size);
- memproc::memcopy<memOut, RunMem>(z0_t_, z0_t, new_size);
- memproc::memcopy<memOut, RunMem>(B_, B, new_size);
- memproc::memcopy<memOut, RunMem>(Cm_, Cm, new_size);
- memproc::memcopy<memOut, RunMem>(Ct_, Ct, new_size);
- memproc::memcopy<memOut, RunMem>(Km_, Km, new_size);
- memproc::memcopy<memOut, RunMem>(Pr_t_inv_, Pr_t_inv, new_size);
+ memproc::memcopy<memOut, MemType::CPU>((void*&)res_sfx->zeta, (void*&)sfx.zeta, new_size);
+ memproc::memcopy<memOut, MemType::CPU>((void*&)res_sfx->Rib, (void*&)sfx.Rib, new_size);
+ memproc::memcopy<memOut, MemType::CPU>((void*&)res_sfx->Re, (void*&)sfx.Re, new_size);
+ memproc::memcopy<memOut, MemType::CPU>((void*&)res_sfx->B, (void*&)sfx.B, new_size);
+ memproc::memcopy<memOut, MemType::CPU>((void*&)res_sfx->z0_m, (void*&)sfx.z0_m, new_size);
+ memproc::memcopy<memOut, MemType::CPU>((void*&)res_sfx->z0_t, (void*&)sfx.z0_t, new_size);
+ memproc::memcopy<memOut, MemType::CPU>((void*&)res_sfx->Rib_conv_lim, (void*&)sfx.Rib_conv_lim, new_size);
+ memproc::memcopy<memOut, MemType::CPU>((void*&)res_sfx->Cm, (void*&)sfx.Cm, new_size);
+ memproc::memcopy<memOut, MemType::CPU>((void*&)res_sfx->Ct, (void*&)sfx.Ct, new_size);
+ memproc::memcopy<memOut, MemType::CPU>((void*&)res_sfx->Km, (void*&)sfx.Km, new_size);
+ memproc::memcopy<memOut, MemType::CPU>((void*&)res_sfx->Pr_t_inv, (void*&)sfx.Pr_t_inv, new_size);
}
}
template class FluxSheba<float, MemType::CPU, MemType::CPU, MemType::CPU>;
-template class FluxSheba<double, MemType::CPU, MemType::CPU, MemType::CPU>;
+template class FluxShebaBase<float, MemType::CPU, MemType::CPU, MemType::CPU>;
#ifdef INCLUDE_CUDA
+ template class FluxShebaBase<float, MemType::GPU, MemType::GPU, MemType::GPU>;
+ template class FluxShebaBase<float, MemType::GPU, MemType::GPU, MemType::CPU>;
+ template class FluxShebaBase<float, MemType::GPU, MemType::CPU, MemType::GPU>;
+ template class FluxShebaBase<float, MemType::CPU, MemType::GPU, MemType::GPU>;
+ template class FluxShebaBase<float, MemType::CPU, MemType::CPU, MemType::GPU>;
+ template class FluxShebaBase<float, MemType::CPU, MemType::GPU, MemType::CPU>;
+ template class FluxShebaBase<float, MemType::GPU, MemType::CPU, MemType::CPU>;
+
template class FluxSheba<float, MemType::GPU, MemType::GPU, MemType::GPU>;
template class FluxSheba<float, MemType::GPU, MemType::GPU, MemType::CPU>;
template class FluxSheba<float, MemType::GPU, MemType::CPU, MemType::GPU>;
@@ -283,12 +138,4 @@ template class FluxSheba<double, MemType::CPU, MemType::CPU, MemType::CPU>;
template class FluxSheba<float, MemType::CPU, MemType::CPU, MemType::GPU>;
template class FluxSheba<float, MemType::CPU, MemType::GPU, MemType::CPU>;
template class FluxSheba<float, MemType::GPU, MemType::CPU, MemType::CPU>;
-
- template class FluxSheba<double, MemType::GPU, MemType::GPU, MemType::GPU>;
- template class FluxSheba<double, MemType::GPU, MemType::GPU, MemType::CPU>;
- template class FluxSheba<double, MemType::GPU, MemType::CPU, MemType::GPU>;
- template class FluxSheba<double, MemType::CPU, MemType::GPU, MemType::GPU>;
- template class FluxSheba<double, MemType::CPU, MemType::CPU, MemType::GPU>;
- template class FluxSheba<double, MemType::CPU, MemType::GPU, MemType::CPU>;
- template class FluxSheba<double, MemType::GPU, MemType::CPU, MemType::CPU>;
#endif
\ No newline at end of file
diff --git a/srcF/sfx_data.f90 b/srcF/sfx_data.f90
index 4d6ceaa686a6ba2f73fa41158c002b6013b5f05a..a3c8a2215aaa489a8189bd00bb3ce9896a1a43cb 100644
--- a/srcF/sfx_data.f90
+++ b/srcF/sfx_data.f90
@@ -137,6 +137,7 @@ module sfx_data
integer(C_INT) :: surface_land
integer(C_INT) :: surface_lake
+ real(C_FLOAT) :: kappa;
real(C_FLOAT) :: gamma_c;
real(C_FLOAT) :: Re_visc_min;
real(C_FLOAT) :: h_charnock;
diff --git a/srcF/sfx_fc_wrapper.F90 b/srcF/sfx_fc_wrapper.F90
index de16ec388125675daba6c2719445c832623174d5..0e84bb3f56d8e6b0134a4e40ea53118d72ff13bc 100644
--- a/srcF/sfx_fc_wrapper.F90
+++ b/srcF/sfx_fc_wrapper.F90
@@ -32,30 +32,6 @@ module C_FUNC
contains
- subroutine set_c_struct_sfx_surface_param_values(surface_param)
- use sfx_data
- use sfx_surface
- implicit none
- type (sfx_surface_param), intent(inout) :: surface_param
- surface_param%surface_ocean = surface_ocean
- surface_param%surface_land = surface_land
- surface_param%surface_lake = surface_lake
-
- surface_param%gamma_c = gamma_c
- surface_param%Re_visc_min = Re_visc_min
- surface_param%h_charnock = h_charnock
- surface_param%c1_charnock = c1_charnock
- surface_param%c2_charnock = c2_charnock
- surface_param%Re_rough_min = Re_rough_min
- surface_param%B1_rough = B1_rough
- surface_param%B2_rough = B2_rough
- surface_param%B3_rough = B3_rough
- surface_param%B4_rough = B4_rough
- surface_param%B_max_lake = B_max_lake
- surface_param%B_max_ocean = B_max_ocean
- surface_param%B_max_land = B_max_land
- end subroutine set_c_struct_sfx_surface_param_values
-
subroutine set_c_struct_sfx_phys_constants_values(constants)
use sfx_data
use sfx_phys_const
diff --git a/srcF/sfx_sheba.f90 b/srcF/sfx_sheba.f90
index e0ae6868c86c8df86afa77b65135040de957bca7..ffa5ae399323d619259a9d845596fad262475432 100644
--- a/srcF/sfx_sheba.f90
+++ b/srcF/sfx_sheba.f90
@@ -90,7 +90,7 @@ contains
call set_meteo_vec_c(meteo, meteo_c)
call set_sfx_vec_c(sfx, sfx_c)
- ! call get_surface_fluxes_sheba(c_loc(sfx_c), c_loc(meteo_c), model_param, surface_param, numerics_c, phys_constants, n)
+ call get_surface_fluxes_sheba(c_loc(sfx_c), c_loc(meteo_c), model_param, surface_param, numerics_c, phys_constants, n)
deallocate(meteo_c)
deallocate(sfx_c)
diff --git a/srcF/sfx_surface.f90 b/srcF/sfx_surface.f90
index d5e369902c2992bfec5620cf08e1e4ea95637fb9..79ebca981a43d1dc6b1bfe91a939f78aacc26b59 100644
--- a/srcF/sfx_surface.f90
+++ b/srcF/sfx_surface.f90
@@ -15,6 +15,9 @@ module sfx_surface
! public interface
! --------------------------------------------------------------------------------
+#if defined(INCLUDE_CXX)
+ public :: set_c_struct_sfx_surface_param_values
+#endif
public :: get_charnock_roughness
public :: get_thermal_roughness
! --------------------------------------------------------------------------------
@@ -104,6 +107,31 @@ contains
end function
+#if defined(INCLUDE_CXX)
+ subroutine set_c_struct_sfx_surface_param_values(surface_param)
+ use sfx_data
+ implicit none
+ type (sfx_surface_param), intent(inout) :: surface_param
+ surface_param%surface_ocean = surface_ocean
+ surface_param%surface_land = surface_land
+ surface_param%surface_lake = surface_lake
+
+ surface_param%kappa = kappa
+ surface_param%gamma_c = gamma_c
+ surface_param%Re_visc_min = Re_visc_min
+ surface_param%h_charnock = h_charnock
+ surface_param%c1_charnock = c1_charnock
+ surface_param%c2_charnock = c2_charnock
+ surface_param%Re_rough_min = Re_rough_min
+ surface_param%B1_rough = B1_rough
+ surface_param%B2_rough = B2_rough
+ surface_param%B3_rough = B3_rough
+ surface_param%B4_rough = B4_rough
+ surface_param%B_max_lake = B_max_lake
+ surface_param%B_max_ocean = B_max_ocean
+ surface_param%B_max_land = B_max_land
+ end subroutine set_c_struct_sfx_surface_param_values
+#endif
! charnock roughness definition
! --------------------------------------------------------------------------------
subroutine get_charnock_roughness(z0_m, u_dyn0, U, h, maxiters)