diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0602036aef85b65ad27a9b5abf6f08477268f42c..4fc203dee267a32d1f4424449f245b2971e2c04a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,7 +1,7 @@
cmake_minimum_required(VERSION 3.19)
option(INCLUDE_CUDA "GPU build in mode" OFF)
-option(INCLUDE_CXX "CXX build in mode" OFF)
+option(USE_CXX "CXX build in mode" OFF)
option(BUILD_DOC "Build documentation" OFF)
option(SFX_CHECK_NAN "Build documentation" OFF)
option(USE_CONFIG_PARSER "Build config parser" OFF)
@@ -31,27 +31,32 @@ set(CMAKE_Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR}/modules)
if(USE_CONFIG_PARSER)
add_subdirectory(parser/)
- list(APPEND RUN_MACRO -DUSE_CONFIG_PARSER)
- set(INCLUDE_CXX ON)
+ add_definitions(-DUSE_CONFIG_PARSER)
+ # list(APPEND RUN_MACRO -DUSE_CONFIG_PARSER)
+ set(USE_CXX ON)
endif(USE_CONFIG_PARSER)
-if(INCLUDE_CXX)
- list(APPEND RUN_MACRO -DINCLUDE_CXX)
-endif(INCLUDE_CXX)
-
if(INCLUDE_CUDA)
+ if(NOT DEFINED CMAKE_CUDA_ARCHITECTURES)
+ message(FATAL_ERROR "
+ CMake will not pass any architecture flags to the compiler
+ because the CUDA architecture is not set. You should specify
+ an architecture: set -DCMAKE_CUDA_ARCHITECTURES=<N>.")
+ endif(NOT DEFINED CMAKE_CUDA_ARCHITECTURES)
+
enable_language(CUDA)
- find_package(CUDA REQUIRED)
- include_directories(${CUDA_INCLUDE_DIRS})
- list(APPEND RUN_MACRO -DINCLUDE_CUDA)
- set(INCLUDE_CXX ON)
+ include_directories(${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES})
+ add_definitions(-DINCLUDE_CUDA)
+ set(USE_CXX ON)
endif(INCLUDE_CUDA)
-if(INCLUDE_CXX)
+if(USE_CXX)
enable_language(C)
enable_language(CXX)
set(CMAKE_CXX_STANDARD 11)
-endif(INCLUDE_CXX)
+ add_definitions(-DINCLUDE_CXX)
+ # list(APPEND RUN_MACRO -DINCLUDE_CXX)
+endif(USE_CXX)
set(SOURCES_F
srcF/sfx_io.f90
@@ -77,43 +82,43 @@ set(HEADERS_F
includeF/sfx_def.fi
)
-if(INCLUDE_CXX)
+if(USE_CXX)
set(SOURCES_C
srcC/sfx_call_cxx.c
)
+ set(HEADERS_C
+ includeC/sfx_data.h
+ )
+
set(SOURCES_CXX
- srcCXX/sfx_surface.cpp
srcCXX/sfx_esm.cpp
- srcCXX/sfx_compute_esm.cpp
srcCXX/sfx_sheba.cpp
srcCXX/sfx_compute_sheba.cpp
srcCXX/sfx_call_class_func.cpp
)
set(HEADERS_CXX
- includeCXX/sfx_surface.h
+ includeCU/sfx_surface.cuh
+ includeCU/sfx_math.cuh
+ includeCU/sfx_esm_compute_subfunc.cuh
includeCXX/sfx_esm.h
- includeCXX/sfx_compute_esm.h
includeCXX/sfx_sheba.h
includeCXX/sfx_compute_sheba.h
includeCXX/sfx_call_class_func.h
)
-endif(INCLUDE_CXX)
+endif(USE_CXX)
if(INCLUDE_CUDA)
set(SOURCES_CU
- srcCU/sfx_surface.cu
- srcCU/sfx_compute_esm.cu
+ srcCU/sfx_esm.cu
srcCU/sfx_compute_sheba.cu
)
set(HEADERS_CU
- includeCU/sfx_surface.cuh
- includeCU/sfx_compute_esm.cuh
includeCU/sfx_compute_sheba.cuh
)
endif(INCLUDE_CUDA)
-if(INCLUDE_CXX OR INCLUDE_CUDA)
+if(USE_CXX OR INCLUDE_CUDA)
set(MEMPROC_SOURCES_CXX
srcCXX/sfx_memory_processing.cpp
)
@@ -130,20 +135,20 @@ if(INCLUDE_CXX OR INCLUDE_CUDA)
includeCU/sfx_memory_processing.cuh
)
endif(INCLUDE_CUDA)
-endif(INCLUDE_CXX OR INCLUDE_CUDA)
+endif(USE_CXX OR INCLUDE_CUDA)
-set(SOURCES ${MEMPROC_HEADERS_CU} ${MEMPROC_SOURCES_CU} ${MEMPROC_HEADERS_CXX} ${MEMPROC_SOURCES_CXX} ${HEADERS_CU} ${SOURCES_CU} ${HEADERS_CXX} ${SOURCES_CXX} ${SOURCES_C} ${HEADERS_F} ${SOURCES_F})
+set(SOURCES ${MEMPROC_HEADERS_CU} ${MEMPROC_SOURCES_CU} ${MEMPROC_HEADERS_CXX} ${MEMPROC_SOURCES_CXX} ${HEADERS_CU} ${SOURCES_CU} ${HEADERS_C} ${HEADERS_CXX} ${SOURCES_CXX} ${SOURCES_C} ${HEADERS_F} ${SOURCES_F})
set(CMAKE_Fortran_FLAGS " -g -fbacktrace -ffpe-trap=zero,overflow,underflow -cpp ")
-if(INCLUDE_CXX OR INCLUDE_CUDA)
+
+if(USE_CXX OR INCLUDE_CUDA)
set(CMAKE_CXX_FLAGS " -g -Wunused-variable ")
set(CMAKE_C_FLAGS " -g ")
set(CMAKE_CUDA_FLAGS " -g ")
-endif(INCLUDE_CXX OR INCLUDE_CUDA)
+endif(USE_CXX OR INCLUDE_CUDA)
add_executable(sfx ${SOURCES})
-add_definitions(${RUN_MACRO})
set_property(TARGET sfx PROPERTY LINKER_LANGUAGE Fortran)
target_include_directories(sfx PUBLIC ${CMAKE_BINARY_DIR}/modules/)
diff --git a/includeC/sfx_data.h b/includeC/sfx_data.h
new file mode 100644
index 0000000000000000000000000000000000000000..2b9e3411ca0b66632779b505fe5653423cf20b62
--- /dev/null
+++ b/includeC/sfx_data.h
@@ -0,0 +1,123 @@
+#pragma once
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+ struct meteoDataType
+ {
+ float h;
+ float U;
+ float dT;
+ float Tsemi;
+ float dQ;
+ float z0_m;
+ };
+
+ struct meteoDataVecTypeC
+ {
+ float *h;
+ float *U;
+ float *dT;
+ float *Tsemi;
+ float *dQ;
+ float *z0_m;
+ };
+
+ struct sfxDataType
+ {
+ float zeta;
+ float Rib;
+ float Re;
+ float B;
+ float z0_m;
+ float z0_t;
+ float Rib_conv_lim;
+ float Cm;
+ float Ct;
+ float Km;
+ float Pr_t_inv;
+ };
+
+ struct sfxDataVecTypeC
+ {
+ float *zeta;
+ float *Rib;
+ float *Re;
+ float *B;
+ float *z0_m;
+ float *z0_t;
+ float *Rib_conv_lim;
+ float *Cm;
+ float *Ct;
+ float *Km;
+ float *Pr_t_inv;
+ };
+
+ // use sfx_esm_param
+ struct sfx_esm_param
+ {
+ float kappa;
+ float Pr_t_0_inv;
+ float Pr_t_inf_inv;
+
+ float alpha_m;
+ float alpha_h;
+ float alpha_h_fix;
+ float beta_m;
+ float beta_h;
+ float Rib_max;
+ };
+
+ struct sfx_esm_numericsTypeC
+ {
+ int maxiters_convection;
+ int maxiters_charnock;
+ };
+
+
+ struct sfx_surface_param
+ {
+ int surface_ocean;
+ int surface_land;
+ int surface_lake;
+
+ float gamma_c;
+ float Re_visc_min;
+ float h_charnock;
+ float c1_charnock;
+ float c2_charnock;
+
+ float Re_rough_min;
+ float B1_rough;
+ float B2_rough;
+ float B3_rough;
+ float B4_rough;
+ float B_max_lake;
+ float B_max_ocean;
+ float B_max_land;
+ };
+
+
+ // TODO: add later after esm implementation
+ // struct sfx_sheba_param
+ // {
+ // float alpha_m;
+ // float alpha_h;
+ // float a_m;
+ // float b_m;
+ // float a_h;
+ // float b_h;
+ // float c_h;
+ // };
+
+ struct sfx_phys_constants
+ {
+ float Pr_m;
+ float g;
+ float nu_air;
+ };
+
+#ifdef __cplusplus
+}
+#endif
\ No newline at end of file
diff --git a/includeCU/sfx_compute_esm.cuh b/includeCU/sfx_compute_esm.cuh
deleted file mode 100644
index df2f94969a65182639b977541485ac4f044c7d2c..0000000000000000000000000000000000000000
--- a/includeCU/sfx_compute_esm.cuh
+++ /dev/null
@@ -1,14 +0,0 @@
-#pragma once
-
-template<typename T>
-void compute_flux_esm_gpu(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_,
- const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_, const T *h_, const T *in_z0_m_,
- const T kappa, const T Pr_t_0_inv, const T Pr_t_inf_inv,
- const T alpha_m, const T alpha_h, const T alpha_h_fix,
- const T beta_m, const T beta_h, const T Rib_max, 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,
- const int maxiters_charnock, const int maxiters_convection,
- const int grid_size);
\ No newline at end of file
diff --git a/includeCU/sfx_esm_compute_subfunc.cuh b/includeCU/sfx_esm_compute_subfunc.cuh
new file mode 100644
index 0000000000000000000000000000000000000000..46812fa9320532355569ec82040a6be9f3a900f5
--- /dev/null
+++ b/includeCU/sfx_esm_compute_subfunc.cuh
@@ -0,0 +1,144 @@
+#pragma once
+#include "../includeC/sfx_data.h"
+#include "../includeCU/sfx_math.cuh"
+
+template<typename T>
+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_esm_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>
+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_esm_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>
+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_esm_param& param,
+ const sfx_esm_numericsTypeC& numerics)
+{
+ 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 <= numerics.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 == numerics.maxiters_convection + 1)
+ break;
+
+ zeta = Rib * psi_m * psi_m / psi_h;
+ }
+}
+
+template<typename T>
+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_esm_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>
+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_esm_param& param,
+ const sfx_esm_numericsTypeC& numerics)
+{
+ 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 <= numerics.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 == numerics.maxiters_convection + 1)
+ break;
+
+ zeta = Rib * psi_m * psi_m / psi_h;
+ }
+}
\ No newline at end of file
diff --git a/includeCU/sfx_math.cuh b/includeCU/sfx_math.cuh
new file mode 100644
index 0000000000000000000000000000000000000000..3f001a76209dbce5de243ffa3598ec397b219cd7
--- /dev/null
+++ b/includeCU/sfx_math.cuh
@@ -0,0 +1,33 @@
+#pragma once
+#ifdef INCLUDE_CUDA
+ #include <cuda.h>
+ #include <cuda_runtime_api.h>
+ #define FUCNTION_DECLARATION_SPECIFIER __host__ __device__
+#else
+ #define FUCNTION_DECLARATION_SPECIFIER
+#endif
+
+namespace sfx_math
+{
+ template<typename T>
+ FUCNTION_DECLARATION_SPECIFIER T min(const T a, const T b);
+
+ template<typename T>
+ FUCNTION_DECLARATION_SPECIFIER T max(const T a, const T b);
+}
+
+template<typename T>
+FUCNTION_DECLARATION_SPECIFIER T sfx_math::min(const T a, const T b)
+{
+ if(a > b)
+ return b;
+ return a;
+}
+
+template<typename T>
+FUCNTION_DECLARATION_SPECIFIER T sfx_math::max(const T a, const T b)
+{
+ if(a > b)
+ return a;
+ return b;
+}
\ No newline at end of file
diff --git a/includeCU/sfx_surface.cuh b/includeCU/sfx_surface.cuh
index 002336c7b2f004c1206b0da8256bf71219d7a2c0..2344f41c42c0c6de0b8805863a1972833b59f954 100644
--- a/includeCU/sfx_surface.cuh
+++ b/includeCU/sfx_surface.cuh
@@ -1,9 +1,55 @@
#pragma once
+#include "../includeC/sfx_data.h"
+#include "../includeCU/sfx_math.cuh"
-// template<typename T>
-// __device__ void get_charnock_roughness(T &z0_m, T &u_dyn0,
-// const T h, const T U,
-// const T kappa,
-// const T h_charnock, const T c1_charnock, const T c2_charnock,
-// const int maxiters);
-
\ No newline at end of file
+template<typename T>
+FUCNTION_DECLARATION_SPECIFIER void get_thermal_roughness(T &z0_t, T &B,
+ const T z0_m, const T Re,
+ const struct sfx_surface_param& param,
+ const int surface_type)
+{
+ // --- define B = log(z0_m / z0_t)
+ B = (Re <= param.Re_rough_min) ? param.B1_rough * log(param.B3_rough * Re) + param.B2_rough :
+ param.B4_rough * (pow(Re, param.B2_rough));
+
+ // --- apply max restriction based on surface type
+ if (surface_type == param.surface_ocean) B = sfx_math::min(B, param.B_max_ocean);
+ if (surface_type == param.surface_lake) B = sfx_math::min(B, param.B_max_land);
+ if (surface_type == param.surface_land) B = sfx_math::min(B, param.B_max_lake);
+
+ // --- define roughness [thermal]
+ z0_t = z0_m / exp(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)
+{
+ 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;
+
+ for (int i = 0; i < numerics.maxiters_charnock; i++)
+ {
+ f = surface_param.c1_charnock - 2.0 * log(Uc);
+ for (int j = 0; j < numerics.maxiters_charnock; j++)
+ {
+ c = (f + 2.0 * log(b)) / model_param.kappa;
+ if (U <= 8.0e0)
+ a = log(1.0 + surface_param.c2_charnock * ( pow(b / Uc, 3) ) ) / model_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 = sfx_math::max(z0_m, T(0.000015e0));
+ Uc = U * log(surface_param.h_charnock / z0_m) / log(h / z0_m);
+ }
+
+ u_dyn0 = Uc / c;
+}
\ No newline at end of file
diff --git a/includeCXX/sfx_call_class_func.h b/includeCXX/sfx_call_class_func.h
index 317f3d92fe009106542225f69b3b2a79a395bc89..973ab7af83735d06e26b164a7eee06be8f2732f3 100644
--- a/includeCXX/sfx_call_class_func.h
+++ b/includeCXX/sfx_call_class_func.h
@@ -1,20 +1,17 @@
#pragma once
+#include "../includeC/sfx_data.h"
#ifdef __cplusplus
extern "C" {
#endif
- void surf_flux_esm_CXX (float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
- float *U_, float *dT_, float *Tsemi_, float *dQ_, float *h_, float *in_z0_m_,
- const float kappa, const float Pr_t_0_inv, const float Pr_t_inf_inv,
- const float alpha_m, const float alpha_h, const float alpha_h_fix,
- const float beta_m, const float beta_h, const float Rib_max, const float Re_rough_min,
- const float B1_rough, const float B2_rough,
- const float B_max_land, const float B_max_ocean, const float B_max_lake,
- const float gamma_c, const float Re_visc_min,
- const float Pr_m, const float nu_air, const float g,
- const int maxiters_charnock, const int maxiters_convection,
- const int grid_size);
+ void surf_flux_esm_CXX(struct sfxDataVecTypeC* sfx,
+ struct meteoDataVecTypeC* meteo,
+ const struct sfx_esm_param* model_param,
+ const struct sfx_surface_param* surface_param,
+ const struct sfx_esm_numericsTypeC* numerics,
+ const struct sfx_phys_constants* constants,
+ const int grid_size);
void surf_flux_sheba_CXX (float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
float *U_, float *dT_, float *Tsemi_, float *dQ_, float *h_, float *in_z0_m_,
diff --git a/includeCXX/sfx_compute_esm.h b/includeCXX/sfx_compute_esm.h
deleted file mode 100644
index 1d7564efa2c9acfe627f8bf75ee044192438b9fe..0000000000000000000000000000000000000000
--- a/includeCXX/sfx_compute_esm.h
+++ /dev/null
@@ -1,14 +0,0 @@
-#pragma once
-
-template<typename T>
-void compute_flux_esm_cpu(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_,
- const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_, const T *h_, const T *in_z0_m_,
- const T kappa, const T Pr_t_0_inv, const T Pr_t_inf_inv,
- const T alpha_m, const T alpha_h, const T alpha_h_fix,
- const T beta_m, const T beta_h, const T Rib_max, 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,
- const int maxiters_charnock, const int maxiters_convection,
- const int grid_size);
\ No newline at end of file
diff --git a/includeCXX/sfx_esm.h b/includeCXX/sfx_esm.h
index 1b7067cee721124781200c70a2838f83d186c6f4..43c54b7c8099331b18fef2fa7395734f19e43d64 100644
--- a/includeCXX/sfx_esm.h
+++ b/includeCXX/sfx_esm.h
@@ -1,43 +1,92 @@
#pragma once
#include "sfx_template_parameters.h"
-#include <cstddef>
+#include "../includeC/sfx_data.h"
+// #include <cstddef>
template<typename T, MemType memIn, MemType memOut, MemType RunMem >
-class FluxEsm
+class FluxEsmBase
{
-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;
+protected:
+ struct sfxDataVecTypeC sfx;
+ struct meteoDataVecTypeC meteo;
+ struct sfx_esm_param model_param;
+ struct sfx_surface_param surface_param;
+ struct sfx_esm_numericsTypeC numerics;
+ struct sfx_phys_constants phys_constants;
- T kappa, Pr_t_0_inv, Pr_t_inf_inv,
- alpha_m, alpha_h, alpha_h_fix,
- beta_m, beta_h,
- Rib_max, 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;
+ struct sfxDataVecTypeC* res_sfx;
int grid_size;
bool ifAllocated;
size_t allocated_size;
public:
- FluxEsm();
- void set_params(const int grid_size, const T kappa, const T Pr_t_0_inv, const T Pr_t_inf_inv,
- const T alpha_m, const T alpha_h, const T alpha_h_fix,
- const T beta_m, const T beta_h, const T Rib_max, 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);
- ~FluxEsm();
+ FluxEsmBase() = default;
+ FluxEsmBase(struct sfxDataVecTypeC* sfx,
+ struct meteoDataVecTypeC* meteo,
+ const struct sfx_esm_param model_param,
+ const struct sfx_surface_param surface_param,
+ const struct sfx_esm_numericsTypeC numerics,
+ const struct sfx_phys_constants phys_constants,
+ const int grid_size);
+ ~FluxEsmBase();
+};
- void compute_flux_esm(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, const int maxiters_convection);
+template<typename T, MemType memIn, MemType memOut, MemType RunMem >
+class FluxEsm : public FluxEsmBase<T, memIn, memOut, RunMem>
+{};
+
+template<typename T, MemType memIn, MemType memOut >
+class FluxEsm<T, memIn, memOut, MemType::CPU> : public FluxEsmBase<T, memIn, memOut, MemType::CPU>
+{
+ using FluxEsmBase<T, memIn, memOut, MemType::CPU>::res_sfx;
+ using FluxEsmBase<T, memIn, memOut, MemType::CPU>::sfx;
+ using FluxEsmBase<T, memIn, memOut, MemType::CPU>::meteo;
+ using FluxEsmBase<T, memIn, memOut, MemType::CPU>::model_param;
+ using FluxEsmBase<T, memIn, memOut, MemType::CPU>::surface_param;
+ using FluxEsmBase<T, memIn, memOut, MemType::CPU>::numerics;
+ using FluxEsmBase<T, memIn, memOut, MemType::CPU>::phys_constants;
+ using FluxEsmBase<T, memIn, memOut, MemType::CPU>::grid_size;
+ using FluxEsmBase<T, memIn, memOut, MemType::CPU>::ifAllocated;
+ using FluxEsmBase<T, memIn, memOut, MemType::CPU>::allocated_size;
+public:
+ FluxEsm() = default;
+ FluxEsm(struct sfxDataVecTypeC* sfx,
+ struct meteoDataVecTypeC* meteo,
+ const struct sfx_esm_param model_param,
+ const struct sfx_surface_param surface_param,
+ const struct sfx_esm_numericsTypeC numerics,
+ const struct sfx_phys_constants phys_constants,
+ const int grid_size) : FluxEsmBase<T, memIn, memOut, MemType::CPU>(sfx, meteo,
+ model_param, surface_param, numerics, phys_constants, grid_size) {};
+ ~FluxEsm() = 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 FluxEsm<T, memIn, memOut, MemType::GPU> : public FluxEsmBase<T, memIn, memOut, MemType::GPU>
+{
+ using FluxEsmBase<T, memIn, memOut, MemType::GPU>::res_sfx;
+ using FluxEsmBase<T, memIn, memOut, MemType::GPU>::sfx;
+ using FluxEsmBase<T, memIn, memOut, MemType::GPU>::meteo;
+ using FluxEsmBase<T, memIn, memOut, MemType::GPU>::model_param;
+ using FluxEsmBase<T, memIn, memOut, MemType::GPU>::surface_param;
+ using FluxEsmBase<T, memIn, memOut, MemType::GPU>::numerics;
+ using FluxEsmBase<T, memIn, memOut, MemType::GPU>::phys_constants;
+ using FluxEsmBase<T, memIn, memOut, MemType::GPU>::grid_size;
+ using FluxEsmBase<T, memIn, memOut, MemType::GPU>::ifAllocated;
+ using FluxEsmBase<T, memIn, memOut, MemType::GPU>::allocated_size;
+public:
+ FluxEsm() = default;
+ FluxEsm(struct sfxDataVecTypeC* sfx,
+ struct meteoDataVecTypeC* meteo,
+ const struct sfx_esm_param model_param,
+ const struct sfx_surface_param surface_param,
+ const struct sfx_esm_numericsTypeC numerics,
+ const struct sfx_phys_constants phys_constants,
+ const int grid_size) : FluxEsmBase<T, memIn, memOut, MemType::GPU>(sfx, meteo,
+ model_param, surface_param, numerics, phys_constants, grid_size) {};
+ ~FluxEsm() = default;
+ void compute_flux();
+};
+#endif
\ No newline at end of file
diff --git a/includeCXX/sfx_surface.h b/includeCXX/sfx_surface.h
deleted file mode 100644
index d82fec0e89e5282782d02690955f5d09b843b1e4..0000000000000000000000000000000000000000
--- a/includeCXX/sfx_surface.h
+++ /dev/null
@@ -1,16 +0,0 @@
-#pragma once
-
-template<typename T>
-void get_charnock_roughness(T &z0_m, T &u_dyn0,
- const T h, const T U,
- const T kappa,
- const T h_charnock, const T c1_charnock, const T c2_charnock,
- const int maxiters);
-
-template<typename T>
-void get_thermal_roughness(T &z0_t, T &B,
- const T z0_m, const T Re,
- const T Re_rough_min,
- const T B1_rough, const T B2_rough, const T B3_rough, const T B4_rough,
- const T B_max_ocean, const T B_max_lake, const T B_max_land,
- const int surface_type);
\ No newline at end of file
diff --git a/srcC/sfx_call_cxx.c b/srcC/sfx_call_cxx.c
index be432977679517756cd92fecc4f560bdca2c2b37..5ec92a4a9eea9636b1ef2ed856ce9700ccc1cb2b 100644
--- a/srcC/sfx_call_cxx.c
+++ b/srcC/sfx_call_cxx.c
@@ -1,31 +1,18 @@
#include <stdlib.h>
#include <stdio.h>
#include "../includeCXX/sfx_call_class_func.h"
+// #include "../includeC/sfx_call_cxx.h"
// -------------------------------------------------------------------------- //
-void get_surface_fluxes_esm (float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
- float *U_, float *dT_, float *Tsemi_, float *dQ_, float *h_, float *in_z0_m_,
- const float *kappa, const float *Pr_t_0_inv, const float *Pr_t_inf_inv,
- const float *alpha_m, const float *alpha_h, const float *alpha_h_fix,
- const float *beta_m, const float *beta_h, const float *Rib_max, const float *Re_rough_min,
- const float *B1_rough, const float *B2_rough,
- const float *B_max_land, const float *B_max_ocean, const float *B_max_lake,
- const float *gamma_c, const float *Re_visc_min,
- const float *Pr_m, const float *nu_air, const float *g,
- const int *maxiters_charnock, const int *maxiters_convection,
- const int *grid_size)
+void get_surface_fluxes_esm (struct sfxDataVecTypeC* sfx,
+ struct meteoDataVecTypeC* meteo,
+ const struct sfx_esm_param* model_param,
+ const struct sfx_surface_param* surface_param,
+ const struct sfx_esm_numericsTypeC* numerics,
+ const struct sfx_phys_constants* constants,
+ const int *grid_size)
{
- surf_flux_esm_CXX ( 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, *Pr_t_inf_inv,
- *alpha_m, *alpha_h, *alpha_h_fix,
- *beta_m, *beta_h, *Rib_max, *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, *maxiters_convection,
- *grid_size);
+ surf_flux_esm_CXX(sfx, meteo, model_param, surface_param, numerics, constants, *grid_size);
}
void get_surface_fluxes_sheba (float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
diff --git a/srcCU/sfx_compute_esm.cu b/srcCU/sfx_compute_esm.cu
deleted file mode 100644
index 354be330976dea575516ee81a927120262738773..0000000000000000000000000000000000000000
--- a/srcCU/sfx_compute_esm.cu
+++ /dev/null
@@ -1,441 +0,0 @@
-#include <iostream>
-#include "../includeCU/sfx_compute_esm.cuh"
-#include "../includeCU/sfx_surface.cuh"
-
-#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
-inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=true)
-{
- if (code != cudaSuccess)
- {
- fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
- if (abort) exit(code);
- }
-}
-
-template<typename T>
-__device__ void get_charnock_roughness(T &z0_m, T &u_dyn0,
- const T h, const T U,
- const T kappa,
- const T h_charnock, const T c1_charnock, const T c2_charnock,
- const int maxiters)
-{
- T Uc, a, b, c, c_min, f;
-
- Uc = U;
- a = 0.0;
- b = 25.0;
- c_min = log(h_charnock) / kappa;
-
- for (int i = 0; i < maxiters; i++)
- {
- f = c1_charnock - 2.0 * log(Uc);
- for (int j = 0; j < maxiters; j++)
- {
- c = (f + 2.0 * log(b)) / kappa;
- if (U <= 8.0e0)
- a = log(1.0 + c2_charnock * ( pow(b / Uc, 3) ) ) / kappa;
- c = max(c - a, c_min);
- b = c;
- }
- z0_m = h_charnock * exp(-c * kappa);
- z0_m = max(z0_m, T(0.000015e0));
- Uc = U * log(h_charnock / z0_m) / log(h / z0_m);
- }
-
- u_dyn0 = Uc / c;
-}
-
-template __device__ void get_charnock_roughness(float &z0_m, float &u_dyn0,
- const float h, const float U,
- const float kappa,
- const float h_charnock, const float c1_charnock, const float c2_charnock,
- const int maxiters);
-template __device__ void get_charnock_roughness(double &z0_m, double &u_dyn0,
- const double h, const double U,
- const double kappa,
- const double h_charnock, const double c1_charnock, const double c2_charnock,
- const int maxiters);
-
-
-template<typename T>
-__device__ 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 T Pr_t_inf_inv, const T Pr_t_0_inv,
- const T alpha_h, const T alpha_m, const T alpha_h_fix)
-{
- T psi_m, psi_h, f_m, f_h, c;
-
- c = pow(Pr_t_inf_inv / Pr_t_0_inv, 4);
- zeta_lim = (2.0 * alpha_h - c * alpha_m - sqrt( (c * alpha_m)*(c * alpha_m) + 4.0 * c * alpha_h * (alpha_h - alpha_m))) / (2.0 * alpha_h*alpha_h);
-
- f_m_lim = pow(1.0 - alpha_m * zeta_lim, 0.25);
- f_h_lim = sqrt(1.0 - 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 - alpha_m * f_m, 0.25);
- f_h = sqrt(1.0 - 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))) / Pr_t_0_inv;
-
- Rib_lim = zeta_lim * psi_h / (psi_m * psi_m);
-}
-
-template __device__ void get_convection_lim(float &zeta_lim, float &Rib_lim, float &f_m_lim, float &f_h_lim,
- const float h0_m, const float h0_t, const float B,
- const float Pr_t_inf_inv, const float Pr_t_0_inv,
- const float alpha_h, const float alpha_m, const float alpha_h_fix);
-template __device__ void get_convection_lim(double &zeta_lim, double &Rib_lim, double &f_m_lim, double &f_h_lim,
- const double h0_m, const double h0_t, const double B,
- const double Pr_t_inf_inv, const double Pr_t_0_inv,
- const double alpha_h, const double alpha_m, const double alpha_h_fix);
-
-template<typename T>
-void __device__ 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 T Pr_t_0_inv, const T beta_m)
-{
- T Rib_coeff, psi0_m, psi0_h, phi, c;
-
- psi0_m = log(h0_m);
- psi0_h = B / psi0_m;
-
- Rib_coeff = beta_m * Rib;
- c = (psi0_h + 1.0) / 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 * beta_m * (1.0 - Rib_coeff));
-
- phi = beta_m * zeta;
-
- psi_m = psi0_m + phi;
- psi_h = (psi0_m + B) / Pr_t_0_inv + phi;
-}
-
-template __device__ void get_psi_stable(float &psi_m, float &psi_h, float &zeta,
- const float Rib, const float h0_m, const float h0_t, const float B,
- const float Pr_t_0_inv, const float beta_m);
-template __device__ void get_psi_stable(double &psi_m, double &psi_h, double &zeta,
- const double Rib, const double h0_m, const double h0_t, const double B,
- const double Pr_t_0_inv, const double beta_m);
-
-template<typename T>
-void __device__ 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 T Pr_t_0_inv,
- const T alpha_h, const T alpha_m, const T alpha_h_fix,
- const int maxiters)
-{
- 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 + 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 - alpha_m * zeta0_m, 0.25);
- f0_h = sqrt(1.0 - 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) / Pr_t_0_inv;
-
- if (i == maxiters + 1)
- break;
-
- zeta = Rib * psi_m * psi_m / psi_h;
- }
-}
-
-template __device__ void get_psi_convection(float &psi_m, float &psi_h, float &zeta,
- const float Rib, const float h0_m, const float h0_t, const float B,
- const float zeta_conv_lim, const float f_m_conv_lim, const float f_h_conv_lim,
- const float Pr_t_0_inv,
- const float alpha_h, const float alpha_m, const float alpha_h_fix,
- const int maxiters);
-template __device__ void get_psi_convection(double &psi_m, double &psi_h, double &zeta,
- const double Rib, const double h0_m, const double h0_t, const double B,
- const double zeta_conv_lim, const double f_m_conv_lim, const double f_h_conv_lim,
- const double Pr_t_0_inv,
- const double alpha_h, const double alpha_m, const double alpha_h_fix,
- const int maxiters);
-
-template<typename T>
-void __device__ get_psi_neutral(T &psi_m, T &psi_h, T &zeta,
- const T h0_m, const T h0_t, const T B,
- const T Pr_t_0_inv)
-{
- zeta = 0.0;
- psi_m = log(h0_m);
- psi_h = log(h0_t) / Pr_t_0_inv;
- if (fabs(B) < 1.0e-10)
- psi_h = psi_m / Pr_t_0_inv;
-}
-
-template __device__ void get_psi_neutral(float &psi_m, float &psi_h, float &zeta,
- const float h0_m, const float h0_t, const float B,
- const float Pr_t_0_inv);
-template __device__ void get_psi_neutral(double &psi_m, double &psi_h, double &zeta,
- const double h0_m, const double h0_t, const double B,
- const double Pr_t_0_inv);
-
-template<typename T>
-void __device__ 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 T Pr_t_0_inv,
- const T alpha_m, const T alpha_h_fix,
- const int maxiters)
-{
- 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 * Pr_t_0_inv * psi_m * psi_m / psi_h;
-
- for (int i = 1; i <= maxiters + 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 - alpha_m * zeta, 0.25e0);
- f_h = sqrt(1.0 - alpha_h_fix * zeta);
-
- f0_m = pow(1.0 - alpha_m * zeta0_m, 0.25e0);
- f0_h = sqrt(1.0 - alpha_h_fix * zeta0_h);
-
- f0_m = max(f0_m, T(1.000001e0));
- f0_h = 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))) / Pr_t_0_inv;
-
- if (i == maxiters + 1)
- break;
-
- zeta = Rib * psi_m * psi_m / psi_h;
- }
-}
-
-template __device__ void get_psi_semi_convection(float &psi_m, float &psi_h, float &zeta,
- const float Rib, const float h0_m, const float h0_t, const float B,
- const float Pr_t_0_inv,
- const float alpha_m, const float alpha_h_fix,
- const int maxiters);
-template __device__ void get_psi_semi_convection(double &psi_m, double &psi_h, double &zeta,
- const double Rib, const double h0_m, const double h0_t, const double B,
- const double Pr_t_0_inv,
- const double alpha_m, const double alpha_h_fix,
- const int maxiters);
-
-template<typename T>
-__global__ void kernel_compute_flux_esm_gpu(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_,
- const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_, const T *h_, const T *in_z0_m_,
- const T kappa, const T Pr_t_0_inv, const T Pr_t_inf_inv,
- const T alpha_m, const T alpha_h, const T alpha_h_fix,
- const T beta_m, const T beta_h, const T Rib_max, 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,
- const int maxiters_charnock, const int maxiters_convection,
- const int grid_size)
-{
- const int index = blockIdx.x * blockDim.x + threadIdx.x;
-
- T h, U, dT, Tsemi, dQ, z0_m;
- T Re, z0_t, B, h0_m, h0_t, u_dyn0, zeta, Rib, zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim, psi_m, psi_h, phi_m, phi_h, Km, Pr_t_inv, Cm, Ct;
- int surface_type;
- T fval;
-
- const T B3_rough = kappa * Pr_m, B4_rough =( 0.14 * ( pow(30.0, B2_rough) ) ) * (pow(Pr_m, 0.8));
- const T h_charnock = 10.0, c1_charnock = log(h_charnock * (g / gamma_c)), c2_charnock = Re_visc_min * nu_air * c1_charnock;
-
- if(index < grid_size)
- {
- U = U_[index];
- Tsemi = Tsemi_[index];
- dT = dT_[index];
- dQ = dQ_[index];
- h = h_[index];
- z0_m = in_z0_m_[index];
-
- if (z0_m < 0.0) surface_type = 0;
- else surface_type = 1;
-
- if (surface_type == 0)
- {
- get_charnock_roughness(z0_m, u_dyn0, h, U, kappa, h_charnock, c1_charnock, c2_charnock, maxiters_charnock);
- h0_m = h / z0_m;
- }
- if (surface_type == 1)
- {
- h0_m = h / z0_m;
- u_dyn0 = U * kappa / log(h0_m);
- }
-
- Re = u_dyn0 * z0_m / nu_air;
-
- if(Re <= Re_rough_min) B = B1_rough * log(B3_rough * Re) + B2_rough;
- else B = B4_rough * (pow(Re, B2_rough));
-
- if (surface_type == 0) B = min(B, B_max_ocean);
- if (surface_type == 1) B = min(B, B_max_land);
- if (surface_type == 2) B = min(B, B_max_lake);
-
- z0_t = z0_m / exp(B);
- h0_t = h / z0_t;
- Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / (U*U);
-
- get_convection_lim(zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim, h0_m, h0_t, B, Pr_t_inf_inv, Pr_t_0_inv, alpha_h, alpha_m, alpha_h_fix);
-
-
- if (Rib > 0.0)
- {
- Rib = min(Rib, Rib_max);
- get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, Pr_t_0_inv, beta_m);
-
- fval = beta_m * zeta;
- phi_m = 1.0 + fval;
- phi_h = 1.0/Pr_t_0_inv + fval;
- }
-
- 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, Pr_t_0_inv, alpha_h, alpha_m, alpha_h_fix, maxiters_convection);
-
- fval = pow(zeta_conv_lim / zeta, 1.0/3.0);
- phi_m = fval / f_m_conv_lim;
- phi_h = fval / (Pr_t_0_inv * f_h_conv_lim);
- }
- else if (Rib > -0.001)
- {
- get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B, Pr_t_0_inv);
-
- phi_m = 1.0;
- phi_h = 1.0 / Pr_t_0_inv;
- }
- else
- {
- get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, Pr_t_0_inv, alpha_m, alpha_h_fix, maxiters_convection);
-
- phi_m = pow(1.0 - alpha_m * zeta, -0.25);
- phi_h = 1.0 / (Pr_t_0_inv * sqrt(1.0 - alpha_h_fix * zeta));
- }
-
- Cm = kappa / psi_m;
- Ct = kappa / psi_h;
-
- Km = kappa * Cm * U * h / phi_m;
- Pr_t_inv = phi_m / phi_h;
-
- zeta_[index] = zeta;
- Rib_[index] = Rib;
- Re_[index] = Re;
- B_[index] = B;
- z0_m_[index] = z0_m;
- z0_t_[index] = z0_t;
- Rib_conv_lim_[index] = Rib_conv_lim;
- Cm_[index] = Cm;
- Ct_[index] = Ct;
- Km_[index] = Km;
- Pr_t_inv_[index] = Pr_t_inv;
- }
-}
-
-template __global__ void kernel_compute_flux_esm_gpu(float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
- const float *U, const float *dt, const float *T_semi, const float *dq, const float *H, const float *in_z0_m,
- const float kappa, const float Pr_t_0_inv, const float Pr_t_inf_inv,
- const float alpha_m, const float alpha_h, const float alpha_h_fix,
- const float beta_m, const float beta_h, const float Rib_max, const float Re_rough_min,
- const float B1_rough, const float B2_rough,
- const float B_max_land, const float B_max_ocean, const float B_max_lake,
- const float gamma_c, const float Re_visc_min,
- const float Pr_m, const float nu_air, const float g,
- const int maxiters_charnock, const int maxiters_convection,
- const int grid_size);
-template __global__ void kernel_compute_flux_esm_gpu(double *zeta_, double *Rib_, double *Re_, double *B_, double *z0_m_, double *z0_t_, double *Rib_conv_lim_, double *Cm_, double *Ct_, double *Km_, double *Pr_t_inv_,
- const double *U, const double *dt, const double *T_semi, const double *dq, const double *H, const double *in_z0_m,
- const double kappa, const double Pr_t_0_inv, const double Pr_t_inf_inv,
- const double alpha_m, const double alpha_h, const double alpha_h_fix,
- const double beta_m, const double beta_h, const double Rib_max, const double Re_rough_min,
- const double B1_rough, const double B2_rough,
- const double B_max_land, const double B_max_ocean, const double B_max_lake,
- const double gamma_c, const double Re_visc_min,
- const double Pr_m, const double nu_air, const double g,
- const int maxiters_charnock, const int maxiters_convection,
- const int grid_size);
-
-template<typename T>
-void compute_flux_esm_gpu(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_,
- const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_, const T *h_, const T *in_z0_m_,
- const T kappa, const T Pr_t_0_inv, const T Pr_t_inf_inv,
- const T alpha_m, const T alpha_h, const T alpha_h_fix,
- const T beta_m, const T beta_h, const T Rib_max, 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,
- const int maxiters_charnock, const int maxiters_convection,
- const int grid_size)
-{
- const int BlockCount = int(ceil(float(grid_size) / 1024.0));
- dim3 cuBlock = dim3(1024, 1, 1);
- dim3 cuGrid = dim3(BlockCount, 1, 1);
-
- kernel_compute_flux_esm_gpu<<<cuGrid, cuBlock>>>(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, Pr_t_inf_inv,
- alpha_m, alpha_h, alpha_h_fix,
- beta_m, beta_h, Rib_max, 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, maxiters_convection,
- grid_size);
-
- gpuErrchk( cudaPeekAtLastError() );
-}
-
-template void compute_flux_esm_gpu(float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
- const float *U, const float *dt, const float *T_semi, const float *dq, const float *H, const float *in_z0_m,
- const float kappa, const float Pr_t_0_inv, const float Pr_t_inf_inv,
- const float alpha_m, const float alpha_h, const float alpha_h_fix,
- const float beta_m, const float beta_h, const float Rib_max, const float Re_rough_min,
- const float B1_rough, const float B2_rough,
- const float B_max_land, const float B_max_ocean, const float B_max_lake,
- const float gamma_c, const float Re_visc_min,
- const float Pr_m, const float nu_air, const float g,
- const int maxiters_charnock, const int maxiters_convection,
- const int grid_size);
-template void compute_flux_esm_gpu(double *zeta_, double *Rib_, double *Re_, double *B_, double *z0_m_, double *z0_t_, double *Rib_conv_lim_, double *Cm_, double *Ct_, double *Km_, double *Pr_t_inv_,
- const double *U, const double *dt, const double *T_semi, const double *dq, const double *H, const double *in_z0_m,
- const double kappa, const double Pr_t_0_inv, const double Pr_t_inf_inv,
- const double alpha_m, const double alpha_h, const double alpha_h_fix,
- const double beta_m, const double beta_h, const double Rib_max, const double Re_rough_min,
- const double B1_rough, const double B2_rough,
- const double B_max_land, const double B_max_ocean, const double B_max_lake,
- const double gamma_c, const double Re_visc_min,
- const double Pr_m, const double nu_air, const double g,
- const int maxiters_charnock, const int maxiters_convection,
- const int grid_size);
\ No newline at end of file
diff --git a/srcCU/sfx_esm.cu b/srcCU/sfx_esm.cu
new file mode 100644
index 0000000000000000000000000000000000000000..c2d7aef58647901e46f277a699ce3dea8573a478
--- /dev/null
+++ b/srcCU/sfx_esm.cu
@@ -0,0 +1,162 @@
+#include <cuda.h>
+#include <cuda_runtime_api.h>
+
+#include "../includeCXX/sfx_esm.h"
+#include "../includeCU/sfx_esm_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_esm_param model_param,
+ const sfx_surface_param surface_param,
+ const sfx_esm_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_esm_param model_param,
+ const sfx_surface_param surface_param,
+ const sfx_esm_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 Re, z0_t, B, h0_m, h0_t, u_dyn0, zeta, Rib, zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim, psi_m, psi_h, phi_m, phi_h, Km, Pr_t_inv, Cm, Ct;
+ int surface_type;
+ T fval;
+
+ 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, model_param, surface_param, numerics);
+ 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);
+
+ h0_t = h / z0_t;
+ Rib = (phys_constants.g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / (U*U);
+
+ get_convection_lim(zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim,
+ h0_m, h0_t, B,
+ model_param);
+
+
+ if (Rib > 0.0)
+ {
+ Rib = sfx_math::min(Rib, model_param.Rib_max);
+ get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model_param);
+
+ fval = model_param.beta_m * zeta;
+ phi_m = 1.0 + fval;
+ phi_h = 1.0/model_param.Pr_t_0_inv + fval;
+ }
+ 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);
+
+ fval = pow(zeta_conv_lim / zeta, 1.0/3.0);
+ phi_m = fval / f_m_conv_lim;
+ phi_h = fval / (model_param.Pr_t_0_inv * f_h_conv_lim);
+ }
+ else if (Rib > -0.001)
+ {
+ get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B, model_param);
+
+ phi_m = 1.0;
+ phi_h = 1.0 / model_param.Pr_t_0_inv;
+ }
+ else
+ {
+ get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model_param, numerics);
+
+ 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));
+ }
+
+ Cm = model_param.kappa / psi_m;
+ Ct = model_param.kappa / psi_h;
+
+ 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] = Rib_conv_lim;
+ 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 FluxEsm<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 FluxEsmBase<float, MemType::GPU, MemType::GPU, MemType::GPU>;
+template class FluxEsmBase<float, MemType::GPU, MemType::GPU, MemType::CPU>;
+template class FluxEsmBase<float, MemType::GPU, MemType::CPU, MemType::GPU>;
+template class FluxEsmBase<float, MemType::CPU, MemType::GPU, MemType::GPU>;
+template class FluxEsmBase<float, MemType::CPU, MemType::CPU, MemType::GPU>;
+template class FluxEsmBase<float, MemType::CPU, MemType::GPU, MemType::CPU>;
+template class FluxEsmBase<float, MemType::GPU, MemType::CPU, MemType::CPU>;
+
+template class FluxEsm<float, MemType::GPU, MemType::GPU, MemType::GPU>;
+template class FluxEsm<float, MemType::GPU, MemType::GPU, MemType::CPU>;
+template class FluxEsm<float, MemType::GPU, MemType::CPU, MemType::GPU>;
+template class FluxEsm<float, MemType::CPU, MemType::GPU, MemType::GPU>;
+template class FluxEsm<float, MemType::CPU, MemType::CPU, MemType::GPU>;
+template class FluxEsm<float, MemType::CPU, MemType::GPU, MemType::CPU>;
+template class FluxEsm<float, MemType::GPU, MemType::CPU, MemType::CPU>;
\ No newline at end of file
diff --git a/srcCU/sfx_surface.cu b/srcCU/sfx_surface.cu
deleted file mode 100644
index d15f64e55dd6ee86701d8b0fa8d2facf84ad72a4..0000000000000000000000000000000000000000
--- a/srcCU/sfx_surface.cu
+++ /dev/null
@@ -1,46 +0,0 @@
-#include "../includeCU/sfx_surface.cuh"
-
-// template<typename T>
-// __device__ void get_charnock_roughness(T &z0_m, T &u_dyn0,
-// const T h, const T U,
-// const T kappa,
-// const T h_charnock, const T c1_charnock, const T c2_charnock,
-// const int maxiters)
-// {
-// T Uc, a, b, c, c_min, f;
-
-// Uc = U;
-// a = 0.0;
-// b = 25.0;
-// c_min = log(h_charnock) / kappa;
-
-// for (int i = 0; i < maxiters; i++)
-// {
-// f = c1_charnock - 2.0 * log(Uc);
-// for (int j = 0; j < maxiters; j++)
-// {
-// c = (f + 2.0 * log(b)) / kappa;
-// if (U <= 8.0e0)
-// a = log(1.0 + c2_charnock * ( pow(b / Uc, 3) ) ) / kappa;
-// c = max(c - a, c_min);
-// b = c;
-// }
-// z0_m = h_charnock * exp(-c * kappa);
-// z0_m = max(z0_m, T(0.000015e0));
-// Uc = U * log(h_charnock / z0_m) / log(h / z0_m);
-// }
-
-// u_dyn0 = Uc / c;
-// }
-
-// template __device__ void get_charnock_roughness(float &z0_m, float &u_dyn0,
-// const float h, const float U,
-// const float kappa,
-// const float h_charnock, const float c1_charnock, const float c2_charnock,
-// const int maxiters);
-// template __device__ void get_charnock_roughness(double &z0_m, double &u_dyn0,
-// const double h, const double U,
-// const double kappa,
-// const double h_charnock, const double c1_charnock, const double c2_charnock,
-// const int maxiters);
-
\ No newline at end of file
diff --git a/srcCXX/sfx_call_class_func.cpp b/srcCXX/sfx_call_class_func.cpp
index 7b307665527d325f3dfab54dedc3d6ed69c1f25e..b61f36b53f852d6f78172726de7505af39df8709 100644
--- a/srcCXX/sfx_call_class_func.cpp
+++ b/srcCXX/sfx_call_class_func.cpp
@@ -6,35 +6,21 @@
#include <vector>
// -------------------------------------------------------------------------- //
-void surf_flux_esm_CXX (float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
- float *U_, float *dT_, float *Tsemi_, float *dQ_, float *h_, float *in_z0_m_,
- const float kappa, const float Pr_t_0_inv, const float Pr_t_inf_inv,
- const float alpha_m, const float alpha_h, const float alpha_h_fix,
- const float beta_m, const float beta_h, const float Rib_max, const float Re_rough_min,
- const float B1_rough, const float B2_rough,
- const float B_max_land, const float B_max_ocean, const float B_max_lake,
- const float gamma_c, const float Re_visc_min,
- const float Pr_m, const float nu_air, const float g,
- const int maxiters_charnock, const int maxiters_convection,
- const int grid_size)
+void surf_flux_esm_CXX (struct sfxDataVecTypeC* sfx,
+ struct meteoDataVecTypeC* meteo,
+ const struct sfx_esm_param* model_param,
+ const struct sfx_surface_param* surface_param,
+ const struct sfx_esm_numericsTypeC* numerics,
+ const struct sfx_phys_constants* constants,
+ const int grid_size)
{
#ifdef INCLUDE_CUDA
- static FluxEsm<float, MemType::CPU, MemType::CPU, MemType::GPU> F;
+ static FluxEsm<float, MemType::CPU, MemType::CPU, MemType::GPU> F(sfx, meteo, *model_param, *surface_param, *numerics, *constants, grid_size);
+ F.compute_flux();
#else
- static FluxEsm<float, MemType::CPU, MemType::CPU, MemType::CPU> F;
+ static FluxEsm<float, MemType::CPU, MemType::CPU, MemType::CPU> F(sfx, meteo, *model_param, *surface_param, *numerics, *constants, grid_size);
+ F.compute_flux();
#endif
-
- F.set_params(grid_size, kappa, Pr_t_0_inv, Pr_t_inf_inv,
- alpha_m, alpha_h, alpha_h_fix,
- beta_m, beta_h, Rib_max, 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);
-
- F.compute_flux_esm(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_,
- maxiters_charnock, maxiters_convection);
}
void surf_flux_sheba_CXX (float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
diff --git a/srcCXX/sfx_compute_esm.cpp b/srcCXX/sfx_compute_esm.cpp
deleted file mode 100644
index 9a30be0e64297d809282d102efba63fdd75e6ff9..0000000000000000000000000000000000000000
--- a/srcCXX/sfx_compute_esm.cpp
+++ /dev/null
@@ -1,330 +0,0 @@
-#include <cmath>
-#include <iostream>
-#include "../includeCXX/sfx_compute_esm.h"
-#include "../includeCXX/sfx_surface.h"
-
-template<typename T>
-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 T Pr_t_inf_inv, const T Pr_t_0_inv,
- const T alpha_h, const T alpha_m, const T alpha_h_fix)
-{
- T psi_m, psi_h, f_m, f_h, c;
-
- c = pow(Pr_t_inf_inv / Pr_t_0_inv, 4);
- zeta_lim = (2.0 * alpha_h - c * alpha_m - sqrt( (c * alpha_m)*(c * alpha_m) + 4.0 * c * alpha_h * (alpha_h - alpha_m))) / (2.0 * alpha_h*alpha_h);
-
- f_m_lim = pow(1.0 - alpha_m * zeta_lim, 0.25);
- f_h_lim = sqrt(1.0 - 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 - alpha_m * f_m, 0.25);
- f_h = sqrt(1.0 - 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))) / Pr_t_0_inv;
-
- Rib_lim = zeta_lim * psi_h / (psi_m * psi_m);
-}
-
-template void get_convection_lim(float &zeta_lim, float &Rib_lim, float &f_m_lim, float &f_h_lim,
- const float h0_m, const float h0_t, const float B,
- const float Pr_t_inf_inv, const float Pr_t_0_inv,
- const float alpha_h, const float alpha_m, const float alpha_h_fix);
-template void get_convection_lim(double &zeta_lim, double &Rib_lim, double &f_m_lim, double &f_h_lim,
- const double h0_m, const double h0_t, const double B,
- const double Pr_t_inf_inv, const double Pr_t_0_inv,
- const double alpha_h, const double alpha_m, const double alpha_h_fix);
-
-template<typename T>
-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 T Pr_t_0_inv, const T beta_m)
-{
- T Rib_coeff, psi0_m, psi0_h, phi, c;
-
- psi0_m = log(h0_m);
- psi0_h = B / psi0_m;
-
- Rib_coeff = beta_m * Rib;
- c = (psi0_h + 1.0) / 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 * beta_m * (1.0 - Rib_coeff));
-
- phi = beta_m * zeta;
-
- psi_m = psi0_m + phi;
- psi_h = (psi0_m + B) / Pr_t_0_inv + phi;
-}
-
-template void get_psi_stable(float &psi_m, float &psi_h, float &zeta,
- const float Rib, const float h0_m, const float h0_t, const float B,
- const float Pr_t_0_inv, const float beta_m);
-template void get_psi_stable(double &psi_m, double &psi_h, double &zeta,
- const double Rib, const double h0_m, const double h0_t, const double B,
- const double Pr_t_0_inv, const double beta_m);
-
-template<typename T>
-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 T Pr_t_0_inv,
- const T alpha_h, const T alpha_m, const T alpha_h_fix,
- const int maxiters)
-{
- 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 + 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 - alpha_m * zeta0_m, 0.25);
- f0_h = sqrt(1.0 - 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) / Pr_t_0_inv;
-
- if (i == maxiters + 1)
- break;
-
- zeta = Rib * psi_m * psi_m / psi_h;
- }
-}
-
-template void get_psi_convection(float &psi_m, float &psi_h, float &zeta,
- const float Rib, const float h0_m, const float h0_t, const float B,
- const float zeta_conv_lim, const float f_m_conv_lim, const float f_h_conv_lim,
- const float Pr_t_0_inv,
- const float alpha_h, const float alpha_m, const float alpha_h_fix,
- const int maxiters);
-template void get_psi_convection(double &psi_m, double &psi_h, double &zeta,
- const double Rib, const double h0_m, const double h0_t, const double B,
- const double zeta_conv_lim, const double f_m_conv_lim, const double f_h_conv_lim,
- const double Pr_t_0_inv,
- const double alpha_h, const double alpha_m, const double alpha_h_fix,
- const int maxiters);
-
-template<typename T>
-void get_psi_neutral(T &psi_m, T &psi_h, T &zeta,
- const T h0_m, const T h0_t, const T B,
- const T Pr_t_0_inv)
-{
- zeta = 0.0;
- psi_m = log(h0_m);
- psi_h = log(h0_t) / Pr_t_0_inv;
- if (fabs(B) < 1.0e-10)
- psi_h = psi_m / Pr_t_0_inv;
-}
-
-template void get_psi_neutral(float &psi_m, float &psi_h, float &zeta,
- const float h0_m, const float h0_t, const float B,
- const float Pr_t_0_inv);
-template void get_psi_neutral(double &psi_m, double &psi_h, double &zeta,
- const double h0_m, const double h0_t, const double B,
- const double Pr_t_0_inv);
-
-template<typename T>
-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 T Pr_t_0_inv,
- const T alpha_m, const T alpha_h_fix,
- const int maxiters)
-{
- 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 * Pr_t_0_inv * psi_m * psi_m / psi_h;
-
- for (int i = 1; i <= maxiters + 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 - alpha_m * zeta, 0.25e0);
- f_h = sqrt(1.0 - alpha_h_fix * zeta);
-
- f0_m = pow(1.0 - alpha_m * zeta0_m, 0.25e0);
- f0_h = sqrt(1.0 - alpha_h_fix * zeta0_h);
-
- f0_m = std::max(f0_m, T(1.000001e0));
- f0_h = std::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))) / Pr_t_0_inv;
-
- if (i == maxiters + 1)
- break;
-
- zeta = Rib * psi_m * psi_m / psi_h;
- }
-}
-
-template void get_psi_semi_convection(float &psi_m, float &psi_h, float &zeta,
- const float Rib, const float h0_m, const float h0_t, const float B,
- const float Pr_t_0_inv,
- const float alpha_m, const float alpha_h_fix,
- const int maxiters);
-template void get_psi_semi_convection(double &psi_m, double &psi_h, double &zeta,
- const double Rib, const double h0_m, const double h0_t, const double B,
- const double Pr_t_0_inv,
- const double alpha_m, const double alpha_h_fix,
- const int maxiters);
-
-template<typename T>
-void compute_flux_esm_cpu(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_,
- const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_, const T *h_, const T *in_z0_m_,
- const T kappa, const T Pr_t_0_inv, const T Pr_t_inf_inv,
- const T alpha_m, const T alpha_h, const T alpha_h_fix,
- const T beta_m, const T beta_h, const T Rib_max, 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,
- const int maxiters_charnock, const int maxiters_convection,
- const int grid_size)
-{
- T h, U, dT, Tsemi, dQ, z0_m;
- T Re, z0_t, B, h0_m, h0_t, u_dyn0, zeta, Rib, zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim, psi_m, psi_h, phi_m, phi_h, Km, Pr_t_inv, Cm, Ct;
- int surface_type;
- T fval;
-
- const T B3_rough = kappa * Pr_m, B4_rough =( 0.14 * ( pow(30.0, B2_rough) ) ) * (pow(Pr_m, 0.8));
- const T h_charnock = 10.0, c1_charnock = log(h_charnock * (g / gamma_c)), c2_charnock = Re_visc_min * nu_air * c1_charnock;
-
- for (int step = 0; step < grid_size; step++)
- {
- U = U_[step];
- Tsemi = Tsemi_[step];
- dT = dT_[step];
- dQ = dQ_[step];
- h = h_[step];
- z0_m = in_z0_m_[step];
-
- if (z0_m < 0.0) surface_type = 0;
- else surface_type = 1;
-
- if (surface_type == 0)
- {
- get_charnock_roughness(z0_m, u_dyn0, h, U, kappa, h_charnock, c1_charnock, c2_charnock, maxiters_charnock);
- h0_m = h / z0_m;
- }
- if (surface_type == 1)
- {
- h0_m = h / z0_m;
- u_dyn0 = U * kappa / log(h0_m);
- }
-
- Re = u_dyn0 * z0_m / nu_air;
-
- if(Re <= Re_rough_min) B = B1_rough * log(B3_rough * Re) + B2_rough;
- else B = B4_rough * (pow(Re, B2_rough));
-
- if (surface_type == 0) B = std::min(B, B_max_ocean);
- if (surface_type == 1) B = std::min(B, B_max_land);
- if (surface_type == 2) B = std::min(B, B_max_lake);
-
- z0_t = z0_m / exp(B);
- h0_t = h / z0_t;
- Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / (U*U);
-
- get_convection_lim(zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim, h0_m, h0_t, B, Pr_t_inf_inv, Pr_t_0_inv, alpha_h, alpha_m, alpha_h_fix);
-
-
- if (Rib > 0.0)
- {
- Rib = std::min(Rib, Rib_max);
- get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, Pr_t_0_inv, beta_m);
-
- fval = beta_m * zeta;
- phi_m = 1.0 + fval;
- phi_h = 1.0/Pr_t_0_inv + fval;
- }
-
- 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, Pr_t_0_inv, alpha_h, alpha_m, alpha_h_fix, maxiters_convection);
-
- fval = pow(zeta_conv_lim / zeta, 1.0/3.0);
- phi_m = fval / f_m_conv_lim;
- phi_h = fval / (Pr_t_0_inv * f_h_conv_lim);
- }
- else if (Rib > -0.001)
- {
- get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B, Pr_t_0_inv);
-
- phi_m = 1.0;
- phi_h = 1.0 / Pr_t_0_inv;
- }
- else
- {
- get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, Pr_t_0_inv, alpha_m, alpha_h_fix, maxiters_convection);
-
- phi_m = pow(1.0 - alpha_m * zeta, -0.25);
- phi_h = 1.0 / (Pr_t_0_inv * sqrt(1.0 - alpha_h_fix * zeta));
- }
-
- Cm = kappa / psi_m;
- Ct = kappa / psi_h;
-
- Km = kappa * Cm * U * h / phi_m;
- Pr_t_inv = phi_m / phi_h;
-
- zeta_[step] = zeta;
- Rib_[step] = Rib;
- Re_[step] = Re;
- B_[step] = B;
- z0_m_[step] = z0_m;
- z0_t_[step] = z0_t;
- Rib_conv_lim_[step] = Rib_conv_lim;
- Cm_[step] = Cm;
- Ct_[step] = Ct;
- Km_[step] = Km;
- Pr_t_inv_[step] = Pr_t_inv;
- }
-}
-
-template void compute_flux_esm_cpu(float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
- const float *U, const float *dt, const float *T_semi, const float *dq, const float *H, const float *in_z0_m,
- const float kappa, const float Pr_t_0_inv, const float Pr_t_inf_inv,
- const float alpha_m, const float alpha_h, const float alpha_h_fix,
- const float beta_m, const float beta_h, const float Rib_max, const float Re_rough_min,
- const float B1_rough, const float B2_rough,
- const float B_max_land, const float B_max_ocean, const float B_max_lake,
- const float gamma_c, const float Re_visc_min,
- const float Pr_m, const float nu_air, const float g,
- const int maxiters_charnock, const int maxiters_convection,
- const int grid_size);
-template void compute_flux_esm_cpu(double *zeta_, double *Rib_, double *Re_, double *B_, double *z0_m_, double *z0_t_, double *Rib_conv_lim_, double *Cm_, double *Ct_, double *Km_, double *Pr_t_inv_,
- const double *U, const double *dt, const double *T_semi, const double *dq, const double *H, const double *in_z0_m,
- const double kappa, const double Pr_t_0_inv, const double Pr_t_inf_inv,
- const double alpha_m, const double alpha_h, const double alpha_h_fix,
- const double beta_m, const double beta_h, const double Rib_max, const double Re_rough_min,
- const double B1_rough, const double B2_rough,
- const double B_max_land, const double B_max_ocean, const double B_max_lake,
- const double gamma_c, const double Re_visc_min,
- const double Pr_m, const double nu_air, const double g,
- const int maxiters_charnock, const int maxiters_convection,
- const int grid_size);
\ No newline at end of file
diff --git a/srcCXX/sfx_compute_sheba.cpp b/srcCXX/sfx_compute_sheba.cpp
index 48463f1a71afdaccde27146f98392f972da6c410..167fcbb65fad8863ca7837721e61ac23ff82c747 100644
--- a/srcCXX/sfx_compute_sheba.cpp
+++ b/srcCXX/sfx_compute_sheba.cpp
@@ -1,7 +1,7 @@
#include <cmath>
#include <iostream>
#include "../includeCXX/sfx_compute_sheba.h"
-#include "../includeCXX/sfx_surface.h"
+// #include "../includeCXX/sfx_surface.h"
template<typename T>
void get_psi_mh(T &psi_m, T &psi_h,
@@ -219,80 +219,80 @@ void compute_flux_sheba_cpu(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_
const int maxiters_charnock,
const int grid_size)
{
- 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;
-
- const T B3_rough = kappa * Pr_m, B4_rough =(0.14 * (pow(30.0, B2_rough))) * (pow(Pr_m, 0.8));
- const T h_charnock = 10.0, c1_charnock = log(h_charnock * (g / gamma_c)), c2_charnock = Re_visc_min * nu_air * c1_charnock;
-
- int surface_type;
-
- for (int step = 0; step < grid_size; step++)
- {
- U = U_[step];
- Tsemi = Tsemi_[step];
- dT = dT_[step];
- dQ = dQ_[step];
- h = h_[step];
- z0_m = in_z0_m_[step];
-
- if (z0_m < 0.0) surface_type = 0;
- else surface_type = 1;
-
- if (surface_type == 0)
- {
- get_charnock_roughness(z0_m, u_dyn0, h, U, kappa, h_charnock, c1_charnock, c2_charnock, maxiters_charnock);
- h0_m = h / z0_m;
- }
- if (surface_type == 1)
- {
- h0_m = h / z0_m;
- u_dyn0 = U * kappa / log(h0_m);
- }
-
- Re = u_dyn0 * z0_m / nu_air;
- get_thermal_roughness(z0_t, B, z0_m, Re, Re_rough_min, B1_rough, B2_rough, B3_rough, B4_rough, B_max_ocean, B_max_lake, B_max_land, surface_type);
-
- // --- define relative height [thermal]
- h0_t = h / z0_t;
-
- // --- define Ri-bulk
- Rib = (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, (g / Tsemi), kappa, Pr_t_0_inv, alpha_m, alpha_h, a_m, a_h, b_m, b_h, c_h, 10);
- // ----------------------------------------------------------------------------
-
- get_phi(phi_m, phi_h, zeta, alpha_m, alpha_h, a_m, a_h, b_m, b_h, c_h);
- // ----------------------------------------------------------------------------
-
- // --- 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 = kappa * Cm * U * h / phi_m;
- Pr_t_inv = phi_m / phi_h;
-
- zeta_[step] = zeta;
- Rib_[step] = Rib;
- Re_[step] = Re;
- B_[step] = B;
- z0_m_[step] = z0_m;
- z0_t_[step] = z0_t;
- Rib_conv_lim_[step] = 0.0;
- Cm_[step] = Cm;
- Ct_[step] = Ct;
- Km_[step] = Km;
- Pr_t_inv_[step] = Pr_t_inv;
- }
+ // 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;
+
+ // const T B3_rough = kappa * Pr_m, B4_rough =(0.14 * (pow(30.0, B2_rough))) * (pow(Pr_m, 0.8));
+ // const T h_charnock = 10.0, c1_charnock = log(h_charnock * (g / gamma_c)), c2_charnock = Re_visc_min * nu_air * c1_charnock;
+
+ // int surface_type;
+
+ // for (int step = 0; step < grid_size; step++)
+ // {
+ // U = U_[step];
+ // Tsemi = Tsemi_[step];
+ // dT = dT_[step];
+ // dQ = dQ_[step];
+ // h = h_[step];
+ // z0_m = in_z0_m_[step];
+
+ // if (z0_m < 0.0) surface_type = 0;
+ // else surface_type = 1;
+
+ // if (surface_type == 0)
+ // {
+ // get_charnock_roughness(z0_m, u_dyn0, h, U, kappa, h_charnock, c1_charnock, c2_charnock, maxiters_charnock);
+ // h0_m = h / z0_m;
+ // }
+ // if (surface_type == 1)
+ // {
+ // h0_m = h / z0_m;
+ // u_dyn0 = U * kappa / log(h0_m);
+ // }
+
+ // Re = u_dyn0 * z0_m / nu_air;
+ // get_thermal_roughness(z0_t, B, z0_m, Re, Re_rough_min, B1_rough, B2_rough, B3_rough, B4_rough, B_max_ocean, B_max_lake, B_max_land, surface_type);
+
+ // // --- define relative height [thermal]
+ // h0_t = h / z0_t;
+
+ // // --- define Ri-bulk
+ // Rib = (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, (g / Tsemi), kappa, Pr_t_0_inv, alpha_m, alpha_h, a_m, a_h, b_m, b_h, c_h, 10);
+ // // ----------------------------------------------------------------------------
+
+ // get_phi(phi_m, phi_h, zeta, alpha_m, alpha_h, a_m, a_h, b_m, b_h, c_h);
+ // // ----------------------------------------------------------------------------
+
+ // // --- 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 = kappa * Cm * U * h / phi_m;
+ // Pr_t_inv = phi_m / phi_h;
+
+ // zeta_[step] = zeta;
+ // Rib_[step] = Rib;
+ // Re_[step] = Re;
+ // B_[step] = B;
+ // z0_m_[step] = z0_m;
+ // z0_t_[step] = z0_t;
+ // Rib_conv_lim_[step] = 0.0;
+ // Cm_[step] = Cm;
+ // Ct_[step] = Ct;
+ // Km_[step] = Km;
+ // Pr_t_inv_[step] = Pr_t_inv;
+ // }
}
template void compute_flux_sheba_cpu(float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
diff --git a/srcCXX/sfx_esm.cpp b/srcCXX/sfx_esm.cpp
index 4e1825a17c1f276aea363c5b8977f6dbd00756da..494cf6c2b602f0f6b570c1e0797fc03a545ea546 100644
--- a/srcCXX/sfx_esm.cpp
+++ b/srcCXX/sfx_esm.cpp
@@ -1,141 +1,128 @@
#include <iostream>
+#include <cmath>
#include "../includeCXX/sfx_esm.h"
-#include "../includeCXX/sfx_compute_esm.h"
#ifdef INCLUDE_CUDA
- #include "../includeCU/sfx_compute_esm.cuh"
#include "../includeCU/sfx_memory_processing.cuh"
#endif
#include "../includeCXX/sfx_memory_processing.h"
+#include "../includeCU/sfx_surface.cuh"
+#include "../includeCU/sfx_esm_compute_subfunc.cuh"
template<typename T, MemType memIn, MemType memOut, MemType RunMem >
-FluxEsm<T, memIn, memOut, RunMem>::FluxEsm()
+FluxEsmBase<T, memIn, memOut, RunMem>::FluxEsmBase(struct sfxDataVecTypeC* sfx_in,
+ struct meteoDataVecTypeC* meteo_in,
+ const struct sfx_esm_param model_param_in,
+ const struct sfx_surface_param surface_param_in,
+ const struct sfx_esm_numericsTypeC numerics_in,
+ const struct sfx_phys_constants phys_constants_in,
+ const int grid_size_in)
{
- kappa = 0, Pr_t_0_inv = 0, Pr_t_inf_inv = 0,
- alpha_m = 0, alpha_h = 0, alpha_h_fix = 0,
- beta_m = 0, beta_h = 0,
- Rib_max = 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;
-}
+ grid_size = grid_size_in;
-template<typename T, MemType memIn, MemType memOut, MemType RunMem >
-void FluxEsm<T, memIn, memOut, RunMem>::set_params(const int grid_size_, const T kappa_, const T Pr_t_0_inv_, const T Pr_t_inf_inv_,
- const T alpha_m_, const T alpha_h_, const T alpha_h_fix_,
- const T beta_m_, const T beta_h_, const T Rib_max_, 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_, Pr_t_inf_inv = Pr_t_inf_inv_,
- alpha_m = alpha_m_, alpha_h = alpha_h_, alpha_h_fix = alpha_h_fix_,
- beta_m = beta_m_, beta_h = beta_h_,
- Rib_max = Rib_max_, 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_,
- gamma_c = gamma_c_,
- Re_visc_min = Re_visc_min_,
- Pr_m = Pr_m_, nu_air = nu_air_, g = g_;
-
- grid_size = grid_size_;
+ model_param = model_param_in;
+ surface_param = surface_param_in;
+ numerics = numerics_in;
+ phys_constants = phys_constants_in;
+ if(RunMem != memOut)
+ res_sfx = sfx_in;
+ else
+ res_sfx = nullptr;
+
if(RunMem != memIn)
{
const size_t new_size = grid_size * sizeof(T);
allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(U), allocated_size, new_size);
+ memproc::realloc<RunMem>((void *&)(meteo.U), allocated_size, new_size);
+ memproc::memcopy<RunMem, memIn>(meteo.U, meteo_in->U, new_size);
allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(dT), allocated_size, new_size);
+ memproc::realloc<RunMem>((void *&)(meteo.dT), allocated_size, new_size);
+ memproc::memcopy<RunMem, memIn>(meteo.dT, meteo_in->dT, new_size);
allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(Tsemi), allocated_size, new_size);
+ memproc::realloc<RunMem>((void *&)(meteo.Tsemi), allocated_size, new_size);
+ memproc::memcopy<RunMem, memIn>(meteo.Tsemi, meteo_in->Tsemi, new_size);
allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(dQ), allocated_size, new_size);
+ memproc::realloc<RunMem>((void *&)(meteo.dQ), allocated_size, new_size);
+ memproc::memcopy<RunMem, memIn>(meteo.dQ, meteo_in->dQ, new_size);
allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(h), allocated_size, new_size);
+ memproc::realloc<RunMem>((void *&)(meteo.h), allocated_size, new_size);
+ memproc::memcopy<RunMem, memIn>(meteo.h, meteo_in->h, new_size);
allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(in_z0_m), allocated_size, new_size);
+ memproc::realloc<RunMem>((void *&)(meteo.z0_m), allocated_size, new_size);
+ memproc::memcopy<RunMem, memIn>(meteo.z0_m, meteo_in->z0_m, new_size);
ifAllocated = true;
}
+ else
+ {
+ meteo = *meteo_in;
+ }
+
if(RunMem != memOut)
{
const size_t new_size = grid_size * sizeof(T);
allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(zeta), allocated_size, new_size);
+ memproc::realloc<RunMem>((void *&)(sfx.zeta), allocated_size, new_size);
allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(Rib), allocated_size, new_size);
+ memproc::realloc<RunMem>((void *&)(sfx.Rib), allocated_size, new_size);
allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(Re), allocated_size, new_size);
+ memproc::realloc<RunMem>((void *&)(sfx.Re), allocated_size, new_size);
allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(Rib_conv_lim), allocated_size, new_size);
+ memproc::realloc<RunMem>((void *&)(sfx.Rib_conv_lim), allocated_size, new_size);
allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(z0_m), allocated_size, new_size);
+ memproc::realloc<RunMem>((void *&)(sfx.z0_m), allocated_size, new_size);
allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(z0_t), allocated_size, new_size);
+ memproc::realloc<RunMem>((void *&)(sfx.z0_t), allocated_size, new_size);
allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(B), allocated_size, new_size);
+ memproc::realloc<RunMem>((void *&)(sfx.B), allocated_size, new_size);
allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(Cm), allocated_size, new_size);
+ memproc::realloc<RunMem>((void *&)(sfx.Cm), allocated_size, new_size);
allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(Ct), allocated_size, new_size);
+ memproc::realloc<RunMem>((void *&)(sfx.Ct), allocated_size, new_size);
allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(Km), allocated_size, new_size);
+ memproc::realloc<RunMem>((void *&)(sfx.Km), allocated_size, new_size);
allocated_size = 0;
- memproc::realloc<RunMem>((void *&)(Pr_t_inv), allocated_size, new_size);
+ memproc::realloc<RunMem>((void *&)(sfx.Pr_t_inv), allocated_size, new_size);
ifAllocated = true;
}
+ else
+ {
+ sfx = *sfx_in;
+ }
}
template<typename T, MemType memIn, MemType memOut, MemType RunMem >
-FluxEsm<T, memIn, memOut, RunMem>::~FluxEsm()
+FluxEsmBase<T, memIn, memOut, RunMem>::~FluxEsmBase()
{
- kappa = 0, Pr_t_0_inv = 0, Pr_t_inf_inv = 0,
- alpha_m = 0, alpha_h = 0, alpha_h_fix = 0,
- beta_m = 0, beta_h = 0,
- Rib_max = 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);
+ memproc::dealloc<RunMem>((void*&)meteo.U);
+ memproc::dealloc<RunMem>((void*&)meteo.dT);
+ memproc::dealloc<RunMem>((void*&)meteo.Tsemi);
+ memproc::dealloc<RunMem>((void*&)meteo.dQ);
+ memproc::dealloc<RunMem>((void*&)meteo.h);
+ memproc::dealloc<RunMem>((void*&)meteo.z0_m);
}
if(RunMem != memOut)
{
- 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);
+ memproc::dealloc<RunMem>((void*&)sfx.zeta);
+ memproc::dealloc<RunMem>((void*&)sfx.Rib);
+ memproc::dealloc<RunMem>((void*&)sfx.Re);
+ memproc::dealloc<RunMem>((void*&)sfx.Rib_conv_lim);
+ memproc::dealloc<RunMem>((void*&)sfx.z0_m);
+ memproc::dealloc<RunMem>((void*&)sfx.z0_t);
+ memproc::dealloc<RunMem>((void*&)sfx.B);
+ memproc::dealloc<RunMem>((void*&)sfx.Cm);
+ memproc::dealloc<RunMem>((void*&)sfx.Ct);
+ memproc::dealloc<RunMem>((void*&)sfx.Km);
+ memproc::dealloc<RunMem>((void*&)sfx.Pr_t_inv);
}
ifAllocated = false;
@@ -143,103 +130,127 @@ FluxEsm<T, memIn, memOut, RunMem>::~FluxEsm()
}
}
-template<typename T, MemType memIn, MemType memOut, MemType RunMem >
-void FluxEsm<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_)
+template<typename T, MemType memIn, MemType memOut >
+void FluxEsm<T, memIn, memOut, MemType::CPU>::compute_flux()
{
- if(RunMem == memIn)
- {
- U = U_;
- dT = dT_;
- Tsemi = Tsemi_;
- dQ = dQ_;
- h = h_;
- in_z0_m = in_z0_m_;
- }
- else
+ T h, U, dT, Tsemi, dQ, z0_m;
+ T Re, z0_t, B, h0_m, h0_t, u_dyn0, zeta, Rib, zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim, psi_m, psi_h, phi_m, phi_h, Km, Pr_t_inv, Cm, Ct;
+ int surface_type;
+ T fval;
+
+ for (int step = 0; step < grid_size; step++)
{
- const size_t new_size = grid_size * sizeof(T);
+ 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];
- 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);
- }
+ surface_type = z0_m < 0.0 ? surface_param.surface_ocean : surface_param.surface_land;
- 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_;
- }
-}
+ if (surface_type == surface_param.surface_ocean)
+ {
+ get_charnock_roughness(z0_m, u_dyn0, h, U, model_param, surface_param, numerics);
+ 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);
+ }
-template<typename T, MemType memIn, MemType memOut, MemType RunMem >
-void FluxEsm<T, memIn, memOut, RunMem>::compute_flux_esm(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, const int maxiters_convection)
-{
- 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_);
+ Re = u_dyn0 * z0_m / phys_constants.nu_air;
+ get_thermal_roughness(z0_t, B, z0_m, Re, surface_param, surface_type);
- if(RunMem == MemType::CPU) compute_flux_esm_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, Pr_t_inf_inv,
- alpha_m, alpha_h, alpha_h_fix,
- beta_m, beta_h, Rib_max, 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, maxiters_convection,
- grid_size);
+ h0_t = h / z0_t;
+ Rib = (phys_constants.g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / (U*U);
-#ifdef INCLUDE_CUDA
- else compute_flux_esm_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, Pr_t_inf_inv,
- alpha_m, alpha_h, alpha_h_fix,
- beta_m, beta_h, Rib_max, 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, maxiters_convection,
- grid_size);
-#endif
+ get_convection_lim(zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim,
+ h0_m, h0_t, B,
+ model_param);
- if(RunMem != memOut)
+
+ if (Rib > 0.0)
+ {
+ Rib = std::min(Rib, model_param.Rib_max);
+ get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model_param);
+
+ fval = model_param.beta_m * zeta;
+ phi_m = 1.0 + fval;
+ phi_h = 1.0/model_param.Pr_t_0_inv + fval;
+ }
+ 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);
+
+ fval = pow(zeta_conv_lim / zeta, 1.0/3.0);
+ phi_m = fval / f_m_conv_lim;
+ phi_h = fval / (model_param.Pr_t_0_inv * f_h_conv_lim);
+ }
+ else if (Rib > -0.001)
+ {
+ get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B, model_param);
+
+ phi_m = 1.0;
+ phi_h = 1.0 / model_param.Pr_t_0_inv;
+ }
+ else
+ {
+ get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model_param, numerics);
+
+ 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));
+ }
+
+ Cm = model_param.kappa / psi_m;
+ Ct = model_param.kappa / psi_h;
+
+ 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] = Rib_conv_lim;
+ sfx.Cm[step] = Cm;
+ sfx.Ct[step] = Ct;
+ sfx.Km[step] = Km;
+ sfx.Pr_t_inv[step] = Pr_t_inv;
+ }
+
+ 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 FluxEsm<float, MemType::CPU, MemType::CPU, MemType::CPU>;
-template class FluxEsm<double, MemType::CPU, MemType::CPU, MemType::CPU>;
+template class FluxEsmBase<float, MemType::CPU, MemType::CPU, MemType::CPU>;
#ifdef INCLUDE_CUDA
+ template class FluxEsmBase<float, MemType::GPU, MemType::GPU, MemType::GPU>;
+ template class FluxEsmBase<float, MemType::GPU, MemType::GPU, MemType::CPU>;
+ template class FluxEsmBase<float, MemType::GPU, MemType::CPU, MemType::GPU>;
+ template class FluxEsmBase<float, MemType::CPU, MemType::GPU, MemType::GPU>;
+ template class FluxEsmBase<float, MemType::CPU, MemType::CPU, MemType::GPU>;
+ template class FluxEsmBase<float, MemType::CPU, MemType::GPU, MemType::CPU>;
+ template class FluxEsmBase<float, MemType::GPU, MemType::CPU, MemType::CPU>;
+
template class FluxEsm<float, MemType::GPU, MemType::GPU, MemType::GPU>;
template class FluxEsm<float, MemType::GPU, MemType::GPU, MemType::CPU>;
template class FluxEsm<float, MemType::GPU, MemType::CPU, MemType::GPU>;
@@ -248,11 +259,11 @@ template class FluxEsm<double, MemType::CPU, MemType::CPU, MemType::CPU>;
template class FluxEsm<float, MemType::CPU, MemType::GPU, MemType::CPU>;
template class FluxEsm<float, MemType::GPU, MemType::CPU, MemType::CPU>;
- template class FluxEsm<double, MemType::GPU, MemType::GPU, MemType::GPU>;
- template class FluxEsm<double, MemType::GPU, MemType::GPU, MemType::CPU>;
- template class FluxEsm<double, MemType::GPU, MemType::CPU, MemType::GPU>;
- template class FluxEsm<double, MemType::CPU, MemType::GPU, MemType::GPU>;
- template class FluxEsm<double, MemType::CPU, MemType::CPU, MemType::GPU>;
- template class FluxEsm<double, MemType::CPU, MemType::GPU, MemType::CPU>;
- template class FluxEsm<double, MemType::GPU, MemType::CPU, MemType::CPU>;
+ // template class FluxEsm<double, MemType::GPU, MemType::GPU, MemType::GPU>;
+ // template class FluxEsm<double, MemType::GPU, MemType::GPU, MemType::CPU>;
+ // template class FluxEsm<double, MemType::GPU, MemType::CPU, MemType::GPU>;
+ // template class FluxEsm<double, MemType::CPU, MemType::GPU, MemType::GPU>;
+ // template class FluxEsm<double, MemType::CPU, MemType::CPU, MemType::GPU>;
+ // template class FluxEsm<double, MemType::CPU, MemType::GPU, MemType::CPU>;
+ // template class FluxEsm<double, MemType::GPU, MemType::CPU, MemType::CPU>;
#endif
\ No newline at end of file
diff --git a/srcCXX/sfx_surface.cpp b/srcCXX/sfx_surface.cpp
deleted file mode 100644
index 676055865f5fa86936b59b45ad588321e865f621..0000000000000000000000000000000000000000
--- a/srcCXX/sfx_surface.cpp
+++ /dev/null
@@ -1,87 +0,0 @@
-#include "../includeCXX/sfx_surface.h"
-#include <cmath>
-#include <iostream>
-
-template<typename T>
-void get_charnock_roughness(T &z0_m, T &u_dyn0,
- const T h, const T U,
- const T kappa,
- const T h_charnock, const T c1_charnock, const T c2_charnock,
- const int maxiters)
-{
- T Uc, a, b, c, c_min, f;
-
- Uc = U;
- a = 0.0;
- b = 25.0;
- c_min = log(h_charnock) / kappa;
-
- for (int i = 0; i < maxiters; i++)
- {
- f = c1_charnock - 2.0 * log(Uc);
- for (int j = 0; j < maxiters; j++)
- {
- c = (f + 2.0 * log(b)) / kappa;
- if (U <= 8.0e0)
- a = log(1.0 + c2_charnock * ( pow(b / Uc, 3) ) ) / kappa;
- c = std::max(c - a, c_min);
- b = c;
- }
- z0_m = h_charnock * exp(-c * kappa);
- z0_m = std::max(z0_m, T(0.000015e0));
- Uc = U * log(h_charnock / z0_m) / log(h / z0_m);
- }
-
- u_dyn0 = Uc / c;
-}
-
-template void get_charnock_roughness(float &z0_m, float &u_dyn0,
- const float h, const float U,
- const float kappa,
- const float h_charnock, const float c1_charnock, const float c2_charnock,
- const int maxiters);
-template void get_charnock_roughness(double &z0_m, double &u_dyn0,
- const double h, const double U,
- const double kappa,
- const double h_charnock, const double c1_charnock, const double c2_charnock,
- const int maxiters);
-
-template<typename T>
-void get_thermal_roughness(T &z0_t, T &B,
- const T z0_m, const T Re,
- const T Re_rough_min,
- const T B1_rough, const T B2_rough, const T B3_rough, const T B4_rough,
- const T B_max_ocean, const T B_max_lake, const T B_max_land,
- const int surface_type)
-{
- // --- define B = log(z0_m / z0_t)
- if (Re <= Re_rough_min)
- B = B1_rough * log(B3_rough * Re) + B2_rough;
- else
- // *: B4 takes into account Re value at z' ~ O(10) z0
- B = B4_rough * (pow(Re, B2_rough));
-
- // --- apply max restriction based on surface type
- if (surface_type == 0)
- B = std::min(B, B_max_ocean);
- else if (surface_type == 2)
- B = std::min(B, B_max_lake);
- else if (surface_type == 1)
- B = std::min(B, B_max_land);
-
- // --- define roughness [thermal]
- z0_t = z0_m / exp(B);
-}
-
-template void get_thermal_roughness(float &z0_t, float &B,
- const float z0_m, const float Re,
- const float Re_rough_min,
- const float B1_rough, const float B2_rough, const float B3_rough, const float B4_rough,
- const float B_max_ocean, const float B_max_lake, const float B_max_land,
- const int surface_type);
-template void get_thermal_roughness(double &z0_t, double &B,
- const double z0_m, const double Re,
- const double Re_rough_min,
- const double B1_rough, const double B2_rough, const double B3_rough, const double B4_rough,
- const double B_max_ocean, const double B_max_lake, const double B_max_land,
- const int surface_type);
\ No newline at end of file
diff --git a/srcF/sfx_data.f90 b/srcF/sfx_data.f90
index c61fd0cf3dff4be2589790dcb730e5cfddc0838c..12dcead2a6bff62e1831c529e210a01f0827acba 100644
--- a/srcF/sfx_data.f90
+++ b/srcF/sfx_data.f90
@@ -4,7 +4,7 @@ module sfx_data
! modules used
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
-
+ use iso_c_binding, only: C_FLOAT, C_INT, C_PTR, C_LOC
! directives list
! --------------------------------------------------------------------------------
implicit none
@@ -14,7 +14,13 @@ module sfx_data
! public interface
! --------------------------------------------------------------------------------
public :: allocate_meteo_vec, deallocate_meteo_vec
+#if defined(INCLUDE_CXX)
+ public :: set_meteo_vec_c
+#endif
public :: allocate_sfx_vec, deallocate_sfx_vec
+#if defined(INCLUDE_CXX)
+ public :: set_sfx_vec_c
+#endif
public :: push_sfx_data
! --------------------------------------------------------------------------------
@@ -23,48 +29,80 @@ module sfx_data
!> @brief meteorological input for surface flux calculation
type, public :: meteoDataType
- real :: h !< constant flux layer height [m]
- real :: U !< abs(wind speed) at 'h' [m/s]
- real :: dT !< difference between potential temperature at 'h' and at surface [K]
- real :: Tsemi !< semi-sum of potential temperature at 'h' and at surface [K]
- real :: dQ !< difference between humidity at 'h' and at surface [g/g]
- real :: z0_m !< surface aerodynamic roughness (should be < 0 for water bodies surface)
+ real(C_FLOAT) :: h !< constant flux layer height [m]
+ real(C_FLOAT) :: U !< abs(wind speed) at 'h' [m/s]
+ real(C_FLOAT) :: dT !< difference between potential temperature at 'h' and at surface [K]
+ real(C_FLOAT) :: Tsemi !< semi-sum of potential temperature at 'h' and at surface [K]
+ real(C_FLOAT) :: dQ !< difference between humidity at 'h' and at surface [g/g]
+ real(C_FLOAT) :: z0_m !< surface aerodynamic roughness (should be < 0 for water bodies surface)
end type
!> @brief meteorological input for surface flux calculation
!> &details using arrays as input
type, public :: meteoDataVecType
-
+#if defined(INCLUDE_CXX)
+ real, pointer :: h(:) !< constant flux layer height [m]
+ real, pointer :: U(:) !< abs(wind speed) at 'h' [m/s]
+ real, pointer :: dT(:) !< difference between potential temperature at 'h' and at surface [K]
+ real, pointer :: Tsemi(:) !< semi-sum of potential temperature at 'h' and at surface [K]
+ real, pointer :: dQ(:) !< difference between humidity at 'h' and at surface [g/g]
+ real, pointer :: z0_m(:) !< surface aerodynamic roughness (should be < 0 for water bodies surface)
+#else
real, allocatable :: h(:) !< constant flux layer height [m]
real, allocatable :: U(:) !< abs(wind speed) at 'h' [m/s]
real, allocatable :: dT(:) !< difference between potential temperature at 'h' and at surface [K]
real, allocatable :: Tsemi(:) !< semi-sum of potential temperature at 'h' and at surface [K]
real, allocatable :: dQ(:) !< difference between humidity at 'h' and at surface [g/g]
real, allocatable :: z0_m(:) !< surface aerodynamic roughness (should be < 0 for water bodies surface)
+#endif
+ end type
+
+#if defined(INCLUDE_CXX)
+ type, public :: meteoDataVecTypeC
+
+ type(C_PTR) :: h !< constant flux layer height [m]
+ type(C_PTR) :: U !< abs(wind speed) at 'h' [m/s]
+ type(C_PTR) :: dT !< difference between potential temperature at 'h' and at surface [K]
+ type(C_PTR) :: Tsemi !< semi-sum of potential temperature at 'h' and at surface [K]
+ type(C_PTR) :: dQ !< difference between humidity at 'h' and at surface [g/g]
+ type(C_PTR) :: z0_m !< surface aerodynamic roughness (should be < 0 for water bodies surface)
end type
+#endif
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
!> @brief surface flux output data
type, public :: sfxDataType
-
- real :: zeta !< = z/L [n/d]
- real :: Rib !< bulk Richardson number [n/d]
- real :: Re !< Reynolds number [n/d]
- real :: B !< = log(z0_m / z0_h) [n/d]
- real :: z0_m !< aerodynamic roughness [m]
- real :: z0_t !< thermal roughness [m]
- real :: Rib_conv_lim !< Ri-bulk convection critical value [n/d]
- real :: Cm !< transfer coefficient for momentum [n/d]
- real :: Ct !< transfer coefficient for heat [n/d]
- real :: Km !< eddy viscosity coeff. at h [m^2/s]
- real :: Pr_t_inv !< inverse turbulent Prandtl number at h [n/d]
+
+ real(C_FLOAT) :: zeta !< = z/L [n/d]
+ real(C_FLOAT) :: Rib !< bulk Richardson number [n/d]
+ real(C_FLOAT) :: Re !< Reynolds number [n/d]
+ real(C_FLOAT) :: B !< = log(z0_m / z0_h) [n/d]
+ real(C_FLOAT) :: z0_m !< aerodynamic roughness [m]
+ real(C_FLOAT) :: z0_t !< thermal roughness [m]
+ real(C_FLOAT) :: Rib_conv_lim !< Ri-bulk convection critical value [n/d]
+ real(C_FLOAT) :: Cm !< transfer coefficient for momentum [n/d]
+ real(C_FLOAT) :: Ct !< transfer coefficient for heat [n/d]
+ real(C_FLOAT) :: Km !< eddy viscosity coeff. at h [m^2/s]
+ real(C_FLOAT) :: Pr_t_inv !< inverse turbulent Prandtl number at h [n/d]
end type
!> @brief surface flux output data
!> &details using arrays as output
type, public :: sfxDataVecType
-
+#if defined(INCLUDE_CXX)
+ real, pointer :: zeta(:) !< = z/L [n/d]
+ real, pointer :: Rib(:) !< bulk Richardson number [n/d]
+ real, pointer :: Re(:) !< Reynolds number [n/d]
+ real, pointer :: B(:) !< = log(z0_m / z0_h) [n/d]
+ real, pointer :: z0_m(:) !< aerodynamic roughness [m]
+ real, pointer :: z0_t(:) !< thermal roughness [m]
+ real, pointer :: Rib_conv_lim(:) !< Ri-bulk convection critical value [n/d]
+ real, pointer :: Cm(:) !< transfer coefficient for momentum [n/d]
+ real, pointer :: Ct(:) !< transfer coefficient for heat [n/d]
+ real, pointer :: Km(:) !< eddy viscosity coeff. at h [m^2/s]
+ real, pointer :: Pr_t_inv(:) !< inverse turbulent Prandtl number at h [n/d]
+#else
real, allocatable :: zeta(:) !< = z/L [n/d]
real, allocatable :: Rib(:) !< bulk Richardson number [n/d]
real, allocatable :: Re(:) !< Reynolds number [n/d]
@@ -76,7 +114,68 @@ module sfx_data
real, allocatable :: Ct(:) !< transfer coefficient for heat [n/d]
real, allocatable :: Km(:) !< eddy viscosity coeff. at h [m^2/s]
real, allocatable :: Pr_t_inv(:) !< inverse turbulent Prandtl number at h [n/d]
+#endif
+ end type
+
+#if defined(INCLUDE_CXX)
+ type, public :: sfxDataVecTypeC
+ type(C_PTR) :: zeta !< = z/L [n/d]
+ type(C_PTR) :: Rib !< bulk Richardson number [n/d]
+ type(C_PTR) :: Re !< Reynolds number [n/d]
+ type(C_PTR) :: B !< = log(z0_m / z0_h) [n/d]
+ type(C_PTR) :: z0_m !< aerodynamic roughness [m]
+ type(C_PTR) :: z0_t !< thermal roughness [m]
+ type(C_PTR) :: Rib_conv_lim !< Ri-bulk convection critical value [n/d]
+ type(C_PTR) :: Cm !< transfer coefficient for momentum [n/d]
+ type(C_PTR) :: Ct !< transfer coefficient for heat [n/d]
+ type(C_PTR) :: Km !< eddy viscosity coeff. at h [m^2/s]
+ type(C_PTR) :: Pr_t_inv !< inverse turbulent Prandtl number at h [n/d]
+ end type
+
+ type, BIND(C), public :: sfx_esm_param
+ real(C_FLOAT) :: kappa
+ real(C_FLOAT) :: Pr_t_0_inv
+ real(C_FLOAT) :: Pr_t_inf_inv
+
+ real(C_FLOAT) :: alpha_m
+ real(C_FLOAT) :: alpha_h
+ real(C_FLOAT) :: alpha_h_fix
+ real(C_FLOAT) :: beta_m
+ real(C_FLOAT) :: beta_h
+ real(C_FLOAT) :: Rib_max
+ end type
+
+ type, BIND(C), public :: sfx_surface_param
+ integer(C_INT) :: surface_ocean
+ integer(C_INT) :: surface_land
+ integer(C_INT) :: surface_lake
+
+ real(C_FLOAT) :: gamma_c;
+ real(C_FLOAT) :: Re_visc_min;
+ real(C_FLOAT) :: h_charnock;
+ real(C_FLOAT) :: c1_charnock;
+ real(C_FLOAT) :: c2_charnock;
+ real(C_FLOAT) :: Re_rough_min;
+ real(C_FLOAT) :: B1_rough;
+ real(C_FLOAT) :: B2_rough;
+ real(C_FLOAT) :: B3_rough;
+ real(C_FLOAT) :: B4_rough;
+ real(C_FLOAT) :: B_max_lake;
+ real(C_FLOAT) :: B_max_ocean;
+ real(C_FLOAT) :: B_max_land;
+ end type
+
+ type, BIND(C), public :: sfx_esm_numericsTypeC
+ integer(C_INT) :: maxiters_convection
+ integer(C_INT) :: maxiters_charnock
end type
+
+ type, BIND(C), public :: sfx_phys_constants
+ real(C_FLOAT) :: Pr_m;
+ real(C_FLOAT) :: g;
+ real(C_FLOAT) :: nu_air;
+ end type
+#endif
! --------------------------------------------------------------------------------
contains
@@ -98,6 +197,24 @@ contains
allocate(meteo%z0_m(n))
end subroutine allocate_meteo_vec
+
+#if defined(INCLUDE_CXX)
+ subroutine set_meteo_vec_c(meteo, meteo_C)
+ !> @brief allocate meteo data vector
+ ! ----------------------------------------------------------------------------
+ type (meteoDataVecType) , intent(in) :: meteo
+ type (meteoDataVecTypeC), pointer, intent(inout) :: meteo_C
+
+ allocate(meteo_C)
+
+ meteo_C%h = c_loc(meteo%h)
+ meteo_C%U = c_loc(meteo%U)
+ meteo_C%dT = c_loc(meteo%dT)
+ meteo_C%Tsemi = c_loc(meteo%Tsemi)
+ meteo_C%dQ = c_loc(meteo%dQ)
+ meteo_C%z0_m = c_loc(meteo%z0_m)
+ end subroutine set_meteo_vec_c
+#endif
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
@@ -139,6 +256,29 @@ contains
allocate(sfx%Pr_t_inv(n))
end subroutine allocate_sfx_vec
+
+#if defined(INCLUDE_CXX)
+ subroutine set_sfx_vec_c(sfx, sfx_C)
+ !> @brief allocate surface fluxes data vector
+ ! ----------------------------------------------------------------------------
+ type (sfxDataVecType), intent(inout) :: sfx
+ type (sfxDataVecTypeC), pointer, intent(inout) :: sfx_C
+
+ allocate(sfx_C)
+
+ sfx_C%zeta = c_loc(sfx%zeta)
+ sfx_C%Rib = c_loc(sfx%Rib)
+ sfx_C%Re = c_loc(sfx%Re)
+ sfx_C%B = c_loc(sfx%B)
+ sfx_C%z0_m = c_loc(sfx%z0_m)
+ sfx_C%z0_t = c_loc(sfx%z0_t)
+ sfx_C%Rib_conv_lim = c_loc(sfx%Rib_conv_lim)
+ sfx_C%Cm = c_loc(sfx%Cm)
+ sfx_C%Ct = c_loc(sfx%Ct)
+ sfx_C%Km = c_loc(sfx%Km)
+ sfx_C%Pr_t_inv = c_loc(sfx%Pr_t_inv)
+ end subroutine set_sfx_vec_c
+#endif
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
diff --git a/srcF/sfx_esm.f90 b/srcF/sfx_esm.f90
index 5c540ccfcf1a6de40a5730bc1a3d787cf879d801..674775849e29a1d03fcc8754f6945730f6f1d2ca 100644
--- a/srcF/sfx_esm.f90
+++ b/srcF/sfx_esm.f90
@@ -11,7 +11,8 @@ module sfx_esm
use sfx_data
use sfx_surface
use sfx_esm_param
-#if defined(INCLUDE_CUDA) || defined(INCLUDE_CXX)
+#if defined(INCLUDE_CXX)
+ use iso_c_binding, only: c_loc
use C_FUNC
#endif
! --------------------------------------------------------------------------------
@@ -38,6 +39,22 @@ module sfx_esm
contains
! --------------------------------------------------------------------------------
+#if defined(INCLUDE_CXX)
+ subroutine set_c_struct_sfx_esm_param_values(sfx_model_param)
+ type (sfx_esm_param), intent(inout) :: sfx_model_param
+ sfx_model_param%kappa = kappa
+ sfx_model_param%Pr_t_0_inv = Pr_t_0_inv
+ sfx_model_param%Pr_t_inf_inv = Pr_t_inf_inv
+
+ sfx_model_param%alpha_m = alpha_m
+ sfx_model_param%alpha_h = alpha_h
+ sfx_model_param%alpha_h_fix = alpha_h_fix
+ sfx_model_param%beta_m = beta_m
+ sfx_model_param%beta_h = beta_h
+ sfx_model_param%Rib_max = Rib_max
+ end subroutine set_c_struct_sfx_esm_param_values
+#endif
+
subroutine get_surface_fluxes_vec(sfx, meteo, numerics, n)
!< @brief surface flux calculation for array data
!< @details contains C/C++ & CUDA interface
@@ -54,19 +71,30 @@ contains
type (sfxDataType) sfx_cell
integer i
! ----------------------------------------------------------------------------
-#if defined(INCLUDE_CUDA) || defined(INCLUDE_CXX)
- call get_surface_fluxes_esm(sfx%zeta, sfx%Rib, sfx%Re, sfx%B, sfx%z0_m, sfx%z0_t, &
- sfx%Rib_conv_lim, sfx%Cm, sfx%Ct, sfx%Km, sfx%Pr_t_inv, &
- meteo%U, meteo%dT, meteo%Tsemi, meteo%dQ, meteo%h, meteo%z0_m, &
- kappa, Pr_t_0_inv, Pr_t_inf_inv, &
- alpha_m, alpha_h, alpha_h_fix, &
- beta_m, beta_h, Rib_max, 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, &
- numerics%maxiters_charnock, numerics%maxiters_convection, &
- n)
+#if defined(INCLUDE_CXX)
+ type (meteoDataVecTypeC), pointer :: meteo_c !< meteorological data (input)
+ type (sfxDataVecTypeC), pointer :: sfx_c !< surface flux data (output)
+ type (sfx_esm_param) :: model_param
+ type (sfx_surface_param) :: surface_param
+ type (sfx_esm_numericsTypeC) :: numerics_c
+ type (sfx_phys_constants) :: phys_constants
+
+ numerics_c%maxiters_convection = numerics%maxiters_convection
+ numerics_c%maxiters_charnock = numerics%maxiters_charnock
+
+ phys_constants%Pr_m = Pr_m;
+ phys_constants%nu_air = nu_air;
+ phys_constants%g = g;
+
+ call set_c_struct_sfx_esm_param_values(model_param)
+ call set_c_struct_sfx_surface_param_values(surface_param)
+ call set_meteo_vec_c(meteo, meteo_c)
+ call set_sfx_vec_c(sfx, sfx_c)
+
+ call get_surface_fluxes_esm(c_loc(sfx_c), c_loc(meteo_c), model_param, surface_param, numerics_c, phys_constants, n)
+
+ deallocate(meteo_c)
+ deallocate(sfx_c)
#else
do i = 1, n
#ifdef SFX_FORCE_DEPRECATED_ESM_CODE
diff --git a/srcF/sfx_fc_wrapper.F90 b/srcF/sfx_fc_wrapper.F90
index 8a9bea82ba137c4a8fb3a6767b7cd157f3ce52af..335056f9ab6e0397a898cfe0b95e6b9727c04ade 100644
--- a/srcF/sfx_fc_wrapper.F90
+++ b/srcF/sfx_fc_wrapper.F90
@@ -1,27 +1,19 @@
module C_FUNC
+
+#if defined(INCLUDE_CXX)
+
INTERFACE
-#if defined(INCLUDE_CUDA) || defined(INCLUDE_CXX)
- SUBROUTINE get_surface_fluxes_esm(zeta, Rib, Re, B, z0_m, z0_t, Rib_conv_lim, &
- Cm, Ct, Km, Prt_inv, &
- U, dT, Tsemi, dQ, h, in_z0_m, &
- kappa, Pr_t_0_inv, Pr_t_inf_inv, &
- alpha_m, alpha_h, alpha_h_fix, &
- beta_m, beta_h, Rib_max, 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, maxiters_convection, &
- grid_size) BIND(C)
- USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_INT
- USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_FLOAT
+ SUBROUTINE get_surface_fluxes_esm(sfx, meteo, model_param, surface_param, numerics, constants, grid_size) BIND(C)
+ use sfx_data
+ USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_INT, C_PTR
IMPLICIT NONE
- INTEGER(C_INT) :: grid_size, maxiters_charnock, maxiters_convection
- REAL(C_FLOAT) :: kappa, Pr_t_0_inv, Pr_t_inf_inv, alpha_m, alpha_h, alpha_h_fix, &
- beta_m, beta_h, Rib_max, 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
- REAL(C_FLOAT), dimension(grid_size) :: U, dT, Tsemi, dQ, h, in_z0_m, zeta, Rib, Re, &
- Rib_conv_lim, z0_m, z0_t, B, Cm, Ct, Km, Prt_inv
+ INTEGER(C_INT) :: grid_size
+ type(C_PTR), value :: sfx
+ type(C_PTR), value :: meteo
+ type(sfx_esm_param) :: model_param
+ type(sfx_surface_param) :: surface_param
+ type(sfx_esm_numericsTypeC) :: numerics
+ type(sfx_phys_constants) :: constants
END SUBROUTINE get_surface_fluxes_esm
SUBROUTINE get_surface_fluxes_sheba(zeta, Rib, Re, B, z0_m, z0_t, Rib_conv_lim, &
@@ -49,7 +41,43 @@ module C_FUNC
REAL(C_FLOAT), dimension(grid_size) :: U, dT, Tsemi, dQ, h, in_z0_m, zeta, Rib, Re, &
Rib_conv_lim, z0_m, z0_t, B, Cm, Ct, Km, Prt_inv
END SUBROUTINE get_surface_fluxes_sheba
-#endif
END INTERFACE
+
+ 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
+ implicit none
+ type (sfx_phys_constants), intent(inout) :: constants
+ constants%Pr_m = Pr_m
+ constants%g = g
+ constants%nu_air = nu_air
+ end subroutine set_c_struct_sfx_phys_constants_values
+#endif
end module C_FUNC