ESM using structs

CMakeLists.txt

 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,31 +135,31 @@ 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/)
+add_executable(drag ${SOURCES})
+set_property(TARGET drag PROPERTY LINKER_LANGUAGE Fortran)
+target_include_directories(drag PUBLIC ${CMAKE_BINARY_DIR}/modules/)
 
 if(USE_CONFIG_PARSER)
-    target_link_libraries(sfx parser_F parser_CXX)
+    target_link_libraries(drag parser_F parser_CXX)
 endif(USE_CONFIG_PARSER)
 
 
 #copy data on post build
 add_custom_command(
-        TARGET sfx POST_BUILD
+        TARGET drag POST_BUILD
         COMMAND ${CMAKE_COMMAND} -E copy_directory
         ${CMAKE_SOURCE_DIR}/data
         ${CMAKE_CURRENT_BINARY_DIR}/data)
\ No newline at end of file

includeC/sfx_data.h

+#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

includeCU/sfx_compute_esm.cuh

-#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

includeCU/sfx_esm_compute_subfunc.cuh

+#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

includeCU/sfx_math.cuh

+#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

includeCU/sfx_surface.cuh

 #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);
+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

includeCXX/sfx_call_class_func.h

 #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, 
+    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_,

includeCXX/sfx_compute_esm.h

-#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

includeCXX/sfx_esm.h

 #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();
+};
+
+template<typename T, MemType memIn, MemType memOut, MemType RunMem >
+class FluxEsm : public FluxEsmBase<T, memIn, memOut, RunMem>
+{};
 
-    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 >
+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_);
+#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

includeCXX/sfx_surface.h

-#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

srcC/sfx_call_cxx.c

 #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, 
+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_,

srcCU/sfx_compute_esm.cu

-#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

srcCU/sfx_esm.cu

+#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

srcCU/sfx_surface.cu

-#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

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, 
+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_,

srcCXX/sfx_compute_esm.cpp

-#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

srcCXX/sfx_compute_sheba.cpp

 #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_,

srcCXX/sfx_esm.cpp

 #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)
+    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++)
     {
-        U = U_;
-        dT = dT_;
-        Tsemi = Tsemi_;
-        dQ = dQ_;
-        h = h_;
-        in_z0_m = in_z0_m_;
+        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];
+
+        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;
         }
-    else
+        if (surface_type == surface_param.surface_land) 
         {
-        const size_t new_size = grid_size * sizeof(T);
-
-        memproc::memcopy<RunMem, memIn>(U, U_, new_size);
-        memproc::memcopy<RunMem, memIn>(dT, dT_, new_size);
-        memproc::memcopy<RunMem, memIn>(Tsemi, Tsemi_, new_size);
-        memproc::memcopy<RunMem, memIn>(dQ, dQ_, new_size);
-        memproc::memcopy<RunMem, memIn>(h, h_, new_size);
-        memproc::memcopy<RunMem, memIn>(in_z0_m, in_z0_m_, new_size);
+            h0_m = h / z0_m;
+            u_dyn0 = U * model_param.kappa / log(h0_m);
         }
 
-    if(RunMem == memOut)
+        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) 
         {
-        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_;
+            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);
         
-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)
+            phi_m = 1.0;
+            phi_h = 1.0 / model_param.Pr_t_0_inv;
+        }
+        else
         {
-    set_data(zeta_, Rib_, Re_, B_, z0_m_, z0_t_, Rib_conv_lim_, Cm_, Ct_, Km_, Pr_t_inv_,
-    U_, dT_, Tsemi_, dQ_, h_, in_z0_m_);
-
-    if(RunMem == MemType::CPU) compute_flux_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);
+            get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model_param, numerics);
             
-#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
+            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));
+        }
 
-    if(RunMem != memOut)
+        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

srcCXX/sfx_surface.cpp

-#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

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,60 +29,132 @@ 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
 
-        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)
+        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)
+    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
 
-        real, allocatable :: zeta(:)            !< = z/L [n/d]
-        real, allocatable :: Rib(:)             !< bulk Richardson number [n/d]
-        real, allocatable :: Re(:)              !< Reynolds number [n/d]
-        real, allocatable :: B(:)               !< = log(z0_m / z0_h) [n/d]
-        real, allocatable :: z0_m(:)            !< aerodynamic roughness [m]
-        real, allocatable :: z0_t(:)            !< thermal roughness [m]
-        real, allocatable :: Rib_conv_lim(:)    !< Ri-bulk convection critical value [n/d]
-        real, allocatable :: Cm(:)              !< transfer coefficient for momentum [n/d]
-        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]
+        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]
+    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 +176,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 +235,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
     ! --------------------------------------------------------------------------------
 
     ! --------------------------------------------------------------------------------