diff --git a/includeC/sfx_data.h b/includeC/sfx_data.h
index dba42c80adea9bbcfc6fe5292abf65f153396e01..bd14366599263bd14dea299dd94dbf22f7e7ddbb 100644
--- a/includeC/sfx_data.h
+++ b/includeC/sfx_data.h
@@ -85,7 +85,7 @@ extern "C" {
     };
 
     // use sfx_esm_param
-    struct sfx_esm_param
+    struct sfx_esm_param_C
     {
         float kappa;
         float Pr_t_0_inv;
@@ -99,14 +99,14 @@ extern "C" {
         float Rib_max;
     };
 
-    struct sfx_esm_numericsTypeC
+    struct sfx_esm_numericsType_C
     {
         int maxiters_convection;
         int maxiters_charnock;
     };
 
     // use sfx_sheba_param
-    struct sfx_sheba_param
+    struct sfx_sheba_param_C
     {
         float kappa;
         float Pr_t_0_inv;
@@ -120,7 +120,7 @@ extern "C" {
         float c_h;
     };
 
-    struct sfx_sheba_numericsTypeC
+    struct sfx_sheba_numericsType_C
     {
         int maxiters_charnock;
     };
diff --git a/includeCXX/sfx_call_class_func.h b/includeCXX/sfx_call_class_func.h
index 8cdb0a77e3e57a1e9f64fccb2dcb1e67a71cb850..5ba0b86eae57417f520579e00452621d40210107 100644
--- a/includeCXX/sfx_call_class_func.h
+++ b/includeCXX/sfx_call_class_func.h
@@ -7,17 +7,17 @@ extern "C" {
 
     void surf_flux_esm_CXX(struct sfxDataVecTypeC* sfx,
                            struct meteoDataVecTypeC* meteo,
-                           const struct sfx_esm_param* model_param, 
+                           const struct sfx_esm_param_C* model_param, 
                            const struct sfx_surface_param* surface_param,
-                           const struct sfx_esm_numericsTypeC* numerics,
+                           const struct sfx_esm_numericsType_C* numerics,
                            const struct sfx_phys_constants* constants,
                            const int grid_size);
 
     void surf_flux_sheba_CXX(struct sfxDataVecTypeC* sfx,
                            struct meteoDataVecTypeC* meteo,
-                           const struct sfx_sheba_param* model_param, 
+                           const struct sfx_sheba_param_C* model_param, 
                            const struct sfx_surface_param* surface_param,
-                           const struct sfx_sheba_numericsTypeC* numerics,
+                           const struct sfx_sheba_numericsType_C* numerics,
                            const struct sfx_phys_constants* constants,
                            const int grid_size);
 #ifdef __cplusplus
diff --git a/includeCXX/sfx_esm.h b/includeCXX/sfx_esm.h
index d4040389d69569e6ba7799606906abe7f13c8db8..34478320b7c7a5b96955182e239f837c9a4e8c39 100644
--- a/includeCXX/sfx_esm.h
+++ b/includeCXX/sfx_esm.h
@@ -14,17 +14,17 @@ public:
     using ModelBase<T, memIn, memOut, RunMem>::ifAllocated;
     using ModelBase<T, memIn, memOut, RunMem>::allocated_size;
 
-    sfx_surface_param surface_param;
-    sfx_phys_constants phys_constants;
-    sfx_esm_param model_param;
-    sfx_esm_numericsTypeC numerics;
+    sfx_surface_param surface;
+    sfx_phys_constants phys;
+    sfx_esm_param_C model;
+    sfx_esm_numericsType_C numerics;
 
     FluxEsmBase(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 sfx_esm_param_C model, 
+                const sfx_surface_param surface,
+                const sfx_esm_numericsType_C numerics,
+                const sfx_phys_constants phys,
                 const int grid_size);
     ~FluxEsmBase();
 };
@@ -39,22 +39,22 @@ class FluxEsm<T, memIn, memOut, MemType::CPU> : public FluxEsmBase<T, memIn, mem
     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>::surface_param;
-    using FluxEsmBase<T, memIn, memOut, MemType::CPU>::phys_constants;
+    using FluxEsmBase<T, memIn, memOut, MemType::CPU>::surface;
+    using FluxEsmBase<T, memIn, memOut, MemType::CPU>::phys;
     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;
-    using FluxEsmBase<T, memIn, memOut, MemType::CPU>::model_param;
+    using FluxEsmBase<T, memIn, memOut, MemType::CPU>::model;
     using FluxEsmBase<T, memIn, memOut, MemType::CPU>::numerics;
 public:
     FluxEsm(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) : FluxEsmBase<T, memIn, memOut, MemType::CPU>(sfx, meteo, model_param, 
-                                       surface_param, numerics, phys_constants, grid_size) {}
+                const sfx_esm_param_C model, 
+                const sfx_surface_param surface,
+                const sfx_esm_numericsType_C numerics,
+                const sfx_phys_constants phys,
+                const int grid_size) : FluxEsmBase<T, memIn, memOut, MemType::CPU>(sfx, meteo, model, 
+                                       surface, numerics, phys, grid_size) {}
     ~FluxEsm() = default;
     void compute_flux();
 };
@@ -66,22 +66,22 @@ class FluxEsm<T, memIn, memOut, MemType::GPU> : public FluxEsmBase<T, memIn, mem
     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>::surface_param;
-    using FluxEsmBase<T, memIn, memOut, MemType::GPU>::phys_constants;
+    using FluxEsmBase<T, memIn, memOut, MemType::GPU>::surface;
+    using FluxEsmBase<T, memIn, memOut, MemType::GPU>::phys;
     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;
-    using FluxEsmBase<T, memIn, memOut, MemType::GPU>::model_param;
+    using FluxEsmBase<T, memIn, memOut, MemType::GPU>::model;
     using FluxEsmBase<T, memIn, memOut, MemType::GPU>::numerics;
 public:
     FluxEsm(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) : FluxEsmBase<T, memIn, memOut, MemType::GPU>(sfx, meteo, model_param, 
-                                       surface_param, numerics, phys_constants, grid_size) {}
+                const sfx_esm_param_C model, 
+                const sfx_surface_param surface,
+                const sfx_esm_numericsType_C numerics,
+                const sfx_phys_constants phys,
+                const int grid_size) : FluxEsmBase<T, memIn, memOut, MemType::GPU>(sfx, meteo, model, 
+                                       surface, numerics, phys, grid_size) {}
     ~FluxEsm() = default;
     void compute_flux();
 };
diff --git a/includeCXX/sfx_sheba.h b/includeCXX/sfx_sheba.h
index 33c93aa76be928524577774d0831e41b83546c76..832f49a5c0ab515c3d6cb1b54e18b51c7bad5e22 100644
--- a/includeCXX/sfx_sheba.h
+++ b/includeCXX/sfx_sheba.h
@@ -14,17 +14,17 @@ public:
     using ModelBase<T, memIn, memOut, RunMem>::ifAllocated;
     using ModelBase<T, memIn, memOut, RunMem>::allocated_size;
 
-    sfx_surface_param surface_param;
-    sfx_phys_constants phys_constants;
-    sfx_sheba_param model_param;
-    sfx_sheba_numericsTypeC numerics;
+    sfx_surface_param surface;
+    sfx_phys_constants phys;
+    sfx_sheba_param_C model;
+    sfx_sheba_numericsType_C numerics;
 
     FluxShebaBase(sfxDataVecTypeC* sfx,
                 meteoDataVecTypeC* meteo,
-                const sfx_sheba_param model_param, 
-                const sfx_surface_param surface_param,
-                const sfx_sheba_numericsTypeC numerics,
-                const sfx_phys_constants phys_constants,
+                const sfx_sheba_param_C model, 
+                const sfx_surface_param surface,
+                const sfx_sheba_numericsType_C numerics,
+                const sfx_phys_constants phys,
                 const int grid_size);
     ~FluxShebaBase();
 };
@@ -39,22 +39,22 @@ class FluxSheba<T, memIn, memOut, MemType::CPU> : public FluxShebaBase<T, memIn,
     using FluxShebaBase<T, memIn, memOut, MemType::CPU>::res_sfx;
     using FluxShebaBase<T, memIn, memOut, MemType::CPU>::sfx;
     using FluxShebaBase<T, memIn, memOut, MemType::CPU>::meteo;
-    using FluxShebaBase<T, memIn, memOut, MemType::CPU>::surface_param;
-    using FluxShebaBase<T, memIn, memOut, MemType::CPU>::phys_constants;
+    using FluxShebaBase<T, memIn, memOut, MemType::CPU>::surface;
+    using FluxShebaBase<T, memIn, memOut, MemType::CPU>::phys;
     using FluxShebaBase<T, memIn, memOut, MemType::CPU>::grid_size;
     using FluxShebaBase<T, memIn, memOut, MemType::CPU>::ifAllocated;
     using FluxShebaBase<T, memIn, memOut, MemType::CPU>::allocated_size;
-    using FluxShebaBase<T, memIn, memOut, MemType::CPU>::model_param;
+    using FluxShebaBase<T, memIn, memOut, MemType::CPU>::model;
     using FluxShebaBase<T, memIn, memOut, MemType::CPU>::numerics;
 public:
     FluxSheba(sfxDataVecTypeC* sfx,
                 meteoDataVecTypeC* meteo,
-                const sfx_sheba_param model_param, 
-                const sfx_surface_param surface_param,
-                const sfx_sheba_numericsTypeC numerics,
-                const sfx_phys_constants phys_constants,
-                const int grid_size) : FluxShebaBase<T, memIn, memOut, MemType::CPU>(sfx, meteo, model_param, 
-                                       surface_param, numerics, phys_constants, grid_size) {}
+                const sfx_sheba_param_C model, 
+                const sfx_surface_param surface,
+                const sfx_sheba_numericsType_C numerics,
+                const sfx_phys_constants phys,
+                const int grid_size) : FluxShebaBase<T, memIn, memOut, MemType::CPU>(sfx, meteo, model, 
+                                       surface, numerics, phys, grid_size) {}
     ~FluxSheba() = default;
     void compute_flux();
 };
@@ -66,22 +66,22 @@ class FluxSheba<T, memIn, memOut, MemType::GPU> : public FluxShebaBase<T, memIn,
     using FluxShebaBase<T, memIn, memOut, MemType::GPU>::res_sfx;
     using FluxShebaBase<T, memIn, memOut, MemType::GPU>::sfx;
     using FluxShebaBase<T, memIn, memOut, MemType::GPU>::meteo;
-    using FluxShebaBase<T, memIn, memOut, MemType::GPU>::surface_param;
-    using FluxShebaBase<T, memIn, memOut, MemType::GPU>::phys_constants;
+    using FluxShebaBase<T, memIn, memOut, MemType::GPU>::surface;
+    using FluxShebaBase<T, memIn, memOut, MemType::GPU>::phys;
     using FluxShebaBase<T, memIn, memOut, MemType::GPU>::grid_size;
     using FluxShebaBase<T, memIn, memOut, MemType::GPU>::ifAllocated;
     using FluxShebaBase<T, memIn, memOut, MemType::GPU>::allocated_size;
-    using FluxShebaBase<T, memIn, memOut, MemType::GPU>::model_param;
+    using FluxShebaBase<T, memIn, memOut, MemType::GPU>::model;
     using FluxShebaBase<T, memIn, memOut, MemType::GPU>::numerics;
 public:
     FluxSheba(sfxDataVecTypeC* sfx,
                 meteoDataVecTypeC* meteo,
-                const sfx_sheba_param model_param, 
-                const sfx_surface_param surface_param,
-                const sfx_sheba_numericsTypeC numerics,
-                const sfx_phys_constants phys_constants,
-                const int grid_size) : FluxShebaBase<T, memIn, memOut, MemType::GPU>(sfx, meteo, model_param, 
-                                       surface_param, numerics, phys_constants, grid_size) {}
+                const sfx_sheba_param_C model, 
+                const sfx_surface_param surface,
+                const sfx_sheba_numericsType_C numerics,
+                const sfx_phys_constants phys,
+                const int grid_size) : FluxShebaBase<T, memIn, memOut, MemType::GPU>(sfx, meteo, model, 
+                                       surface, numerics, phys, grid_size) {}
     ~FluxSheba() = default;
     void compute_flux();
 };
diff --git a/srcC/sfx_call_cxx.c b/srcC/sfx_call_cxx.c
index 208a8ed67505caacbadf5cc8794d5ebb12dfc65a..7cac2c678fc93526424cac569c5cdf983f262ac2 100644
--- a/srcC/sfx_call_cxx.c
+++ b/srcC/sfx_call_cxx.c
@@ -6,9 +6,9 @@
 
 void get_surface_fluxes_esm (struct sfxDataVecTypeC* sfx,
                                  struct meteoDataVecTypeC* meteo,
-                                 const struct sfx_esm_param* model_param, 
+                                 const struct sfx_esm_param_C* model_param, 
                                  const struct sfx_surface_param* surface_param,
-                                 const struct sfx_esm_numericsTypeC* numerics,
+                                 const struct sfx_esm_numericsType_C* numerics,
                                  const struct sfx_phys_constants* constants,
                                  const int *grid_size)
 {
@@ -17,9 +17,9 @@ void get_surface_fluxes_esm (struct sfxDataVecTypeC* sfx,
 
 void get_surface_fluxes_sheba (struct sfxDataVecTypeC* sfx,
                                  struct meteoDataVecTypeC* meteo,
-                                 const struct sfx_sheba_param* model_param, 
+                                 const struct sfx_sheba_param_C* model_param, 
                                  const struct sfx_surface_param* surface_param,
-                                 const struct sfx_sheba_numericsTypeC* numerics,
+                                 const struct sfx_sheba_numericsType_C* numerics,
                                  const struct sfx_phys_constants* constants,
                                  const int *grid_size)
 {
diff --git a/srcCU/sfx_esm.cu b/srcCU/sfx_esm.cu
index f8028d63a82ee44d920518671ce728799a2030e6..e8aba38f7bf3efc93dec8f530a44603ba0b68ea5 100644
--- a/srcCU/sfx_esm.cu
+++ b/srcCU/sfx_esm.cu
@@ -11,20 +11,20 @@ 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 sfx_esm_param_C model, 
+        const sfx_surface_param surface,
+        const sfx_esm_numericsType_C numerics,
+        const sfx_phys_constants phys,
         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 sfx_esm_param_C model, 
+    const sfx_surface_param surface,
+    const sfx_esm_numericsType_C numerics,
+    const sfx_phys_constants phys,
     const int grid_size)
 {
     const int index = blockIdx.x * blockDim.x + threadIdx.x;
@@ -42,66 +42,66 @@ __global__ void sfx_kernel::compute_flux(sfxDataVecTypeC sfx,
         h = meteo.h[index];
         z0_m = meteo.z0_m[index];
 
-        surface_type = z0_m < 0.0 ? surface_param.surface_ocean : surface_param.surface_land;
+        surface_type = z0_m < 0.0 ? surface.surface_ocean : surface.surface_land;
 
-        if (surface_type == surface_param.surface_ocean)
+        if (surface_type == surface.surface_ocean)
         {
-            get_charnock_roughness(z0_m, u_dyn0, U, h, surface_param, numerics.maxiters_charnock);
+            get_charnock_roughness(z0_m, u_dyn0, U, h, surface, numerics.maxiters_charnock);
             h0_m = h / z0_m;
         }
-        if (surface_type == surface_param.surface_land) 
+        if (surface_type == surface.surface_land) 
         {
             h0_m = h / z0_m;
-            u_dyn0 = U * model_param.kappa / log(h0_m);
+            u_dyn0 = U * model.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);
+        Re = u_dyn0 * z0_m / phys.nu_air;
+        get_thermal_roughness(z0_t, B, z0_m, Re, surface, surface_type);
 
         h0_t = h / z0_t;
-        Rib = (phys_constants.g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / (U*U);
+        Rib = (phys.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);
+                            model);
 
 
         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);
+            Rib = sfx_math::min(Rib, model.Rib_max);
+            get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model);
 
-            fval = model_param.beta_m * zeta;
+            fval = model.beta_m * zeta;
             phi_m = 1.0 + fval;
-            phi_h = 1.0/model_param.Pr_t_0_inv + fval;
+            phi_h = 1.0/model.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.maxiters_convection);
+            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, numerics.maxiters_convection);
 
             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);
+            phi_h = fval / (model.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);
+            get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B, model);
         
             phi_m = 1.0;
-            phi_h = 1.0 / model_param.Pr_t_0_inv;
+            phi_h = 1.0 / model.Pr_t_0_inv;
         }
         else
         {
-            get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model_param, numerics.maxiters_convection);
+            get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model, numerics.maxiters_convection);
             
-            phi_m = pow(1.0 - model_param.alpha_m * zeta, -0.25);
-            phi_h = 1.0 / (model_param.Pr_t_0_inv * sqrt(1.0 - model_param.alpha_h_fix * zeta));
+            phi_m = pow(1.0 - model.alpha_m * zeta, -0.25);
+            phi_h = 1.0 / (model.Pr_t_0_inv * sqrt(1.0 - model.alpha_h_fix * zeta));
         }
 
-        Cm = model_param.kappa / psi_m;
-        Ct = model_param.kappa / psi_h;
+        Cm = model.kappa / psi_m;
+        Ct = model.kappa / psi_h;
 
-        Km = model_param.kappa * Cm * U * h / phi_m;
+        Km = model.kappa * Cm * U * h / phi_m;
         Pr_t_inv = phi_m / phi_h;
 
         sfx.zeta[index]         = zeta;
@@ -125,8 +125,8 @@ void FluxEsm<T, memIn, memOut, MemType::GPU>::compute_flux()
     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);
+    sfx_kernel::compute_flux<T><<<cuGrid, cuBlock>>>(sfx, meteo, model, 
+                                                    surface, numerics, phys, grid_size);
 
     if(MemType::GPU != memOut)
     {
diff --git a/srcCU/sfx_sheba.cu b/srcCU/sfx_sheba.cu
index 2161b18798a531a8c81135c9a7e4c4058ff8bb45..e31f93cd34e2a582eabf301e58ccc7af7c2bae06 100644
--- a/srcCU/sfx_sheba.cu
+++ b/srcCU/sfx_sheba.cu
@@ -11,20 +11,20 @@ namespace sfx_kernel
     template<typename T>
     __global__ void compute_flux(sfxDataVecTypeC sfx,
         meteoDataVecTypeC meteo,
-        const sfx_sheba_param model_param, 
-        const sfx_surface_param surface_param,
-        const sfx_sheba_numericsTypeC numerics,
-        const sfx_phys_constants phys_constants,
+        const sfx_sheba_param_C model, 
+        const sfx_surface_param surface,
+        const sfx_sheba_numericsType_C numerics,
+        const sfx_phys_constants phys,
         const int grid_size);
 }
 
 template<typename T>
 __global__ void sfx_kernel::compute_flux(sfxDataVecTypeC sfx,
     meteoDataVecTypeC meteo,
-    const sfx_sheba_param model_param, 
-    const sfx_surface_param surface_param,
-    const sfx_sheba_numericsTypeC numerics,
-    const sfx_phys_constants phys_constants,
+    const sfx_sheba_param_C model, 
+    const sfx_surface_param surface,
+    const sfx_sheba_numericsType_C numerics,
+    const sfx_phys_constants phys,
     const int grid_size)
 {
     const int index = blockIdx.x * blockDim.x + threadIdx.x;
@@ -43,35 +43,35 @@ __global__ void sfx_kernel::compute_flux(sfxDataVecTypeC sfx,
         h = meteo.h[index];
         z0_m = meteo.z0_m[index];
 
-        surface_type = z0_m < 0.0 ? surface_param.surface_ocean : surface_param.surface_land;
+        surface_type = z0_m < 0.0 ? surface.surface_ocean : surface.surface_land;
 
-        if (surface_type == surface_param.surface_ocean) 
+        if (surface_type == surface.surface_ocean) 
         {
-            get_charnock_roughness(z0_m, u_dyn0, U, h, surface_param, numerics.maxiters_charnock);
+            get_charnock_roughness(z0_m, u_dyn0, U, h, surface, numerics.maxiters_charnock);
             h0_m = h / z0_m;
         }
-        if (surface_type == surface_param.surface_land) 
+        if (surface_type == surface.surface_land) 
         {
             h0_m = h / z0_m;
-            u_dyn0 = U * model_param.kappa / log(h0_m);
+            u_dyn0 = U * model.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);
+        Re = u_dyn0 * z0_m / phys.nu_air;
+        get_thermal_roughness(z0_t, B, z0_m, Re, surface, surface_type);
 
         // --- define relative height [thermal]
         h0_t = h / z0_t;
 
         // --- define Ri-bulk
-        Rib = (phys_constants.g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / (U*U);
+        Rib = (phys.g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / (U*U);
 
         // --- get the fluxes
         // ----------------------------------------------------------------------------
         get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, 
-            U, Tsemi, dT, dQ, h, z0_m, z0_t, (phys_constants.g / Tsemi), model_param, 10);
+            U, Tsemi, dT, dQ, h, z0_m, z0_t, (phys.g / Tsemi), model, 10);
         // ----------------------------------------------------------------------------
 
-        get_phi(phi_m, phi_h, zeta, model_param);
+        get_phi(phi_m, phi_h, zeta, model);
         // ----------------------------------------------------------------------------
 
         // --- define transfer coeff. (momentum) & (heat)
@@ -83,7 +83,7 @@ __global__ void sfx_kernel::compute_flux(sfxDataVecTypeC sfx,
             Ct = Tdyn / dT;
 
         // --- define eddy viscosity & inverse Prandtl number
-        Km = model_param.kappa * Cm * U * h / phi_m;
+        Km = model.kappa * Cm * U * h / phi_m;
         Pr_t_inv = phi_m / phi_h;
 
         sfx.zeta[index]         = zeta;
@@ -107,8 +107,8 @@ void FluxSheba<T, memIn, memOut, MemType::GPU>::compute_flux()
     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);
+    sfx_kernel::compute_flux<T><<<cuGrid, cuBlock>>>(sfx, meteo, model, 
+                                                    surface, numerics, phys, grid_size);
 
     if(MemType::GPU != memOut)
     {
diff --git a/srcCXX/sfx_call_class_func.cpp b/srcCXX/sfx_call_class_func.cpp
index 9659c17b691ecb848a596bbf94de4accd181506f..f021fa0bc001f20184152d108feca9558cdaaa9e 100644
--- a/srcCXX/sfx_call_class_func.cpp
+++ b/srcCXX/sfx_call_class_func.cpp
@@ -8,9 +8,9 @@
 // -------------------------------------------------------------------------- //
 void surf_flux_esm_CXX (sfxDataVecTypeC* sfx,
                            meteoDataVecTypeC* meteo,
-                           const sfx_esm_param* model_param, 
+                           const sfx_esm_param_C* model_param, 
                            const sfx_surface_param* surface_param,
-                           const sfx_esm_numericsTypeC* numerics,
+                           const sfx_esm_numericsType_C* numerics,
                            const sfx_phys_constants* constants,
                            const int grid_size)
 {
@@ -25,9 +25,9 @@ void surf_flux_esm_CXX (sfxDataVecTypeC* sfx,
 
 void surf_flux_sheba_CXX (sfxDataVecTypeC* sfx,
                            meteoDataVecTypeC* meteo,
-                           const sfx_sheba_param* model_param, 
+                           const sfx_sheba_param_C* model_param, 
                            const sfx_surface_param* surface_param,
-                           const sfx_sheba_numericsTypeC* numerics,
+                           const sfx_sheba_numericsType_C* numerics,
                            const sfx_phys_constants* constants,
                            const int grid_size)
 {
diff --git a/srcCXX/sfx_esm.cpp b/srcCXX/sfx_esm.cpp
index c324e80f1fa91fd1fa26e7c5acfaa7dbd1af1713..b704d99aa8114c59fec6e200394a193b98ac5ad3 100644
--- a/srcCXX/sfx_esm.cpp
+++ b/srcCXX/sfx_esm.cpp
@@ -13,15 +13,15 @@
 template<typename T, MemType memIn, MemType memOut, MemType RunMem >
 FluxEsmBase<T, memIn, memOut, RunMem>::FluxEsmBase(sfxDataVecTypeC* sfx_in,
                 meteoDataVecTypeC* meteo_in,
-                const sfx_esm_param model_param_in, 
+                const sfx_esm_param_C model_param_in, 
                 const sfx_surface_param surface_param_in,
-                const sfx_esm_numericsTypeC numerics_in,
+                const sfx_esm_numericsType_C numerics_in,
                 const sfx_phys_constants phys_constants_in,
                 const int grid_size_in) : ModelBase<T, memIn, memOut, RunMem>(sfx_in, meteo_in, grid_size_in)
 {
-    surface_param = surface_param_in;
-    phys_constants = phys_constants_in;
-    model_param = model_param_in;
+    surface = surface_param_in;
+    phys = phys_constants_in;
+    model = model_param_in;
     numerics = numerics_in;
 }
 
@@ -45,66 +45,66 @@ void FluxEsm<T, memIn, memOut, MemType::CPU>::compute_flux()
         h = meteo.h[step];
         z0_m = meteo.z0_m[step];
 
-        surface_type = z0_m < 0.0 ? surface_param.surface_ocean : surface_param.surface_land;
+        surface_type = z0_m < 0.0 ? surface.surface_ocean : surface.surface_land;
 
-        if (surface_type == surface_param.surface_ocean)
+        if (surface_type == surface.surface_ocean)
         {
-            get_charnock_roughness(z0_m, u_dyn0, U, h, surface_param, numerics.maxiters_charnock);
+            get_charnock_roughness(z0_m, u_dyn0, U, h, surface, numerics.maxiters_charnock);
             h0_m = h / z0_m;
         }
-        if (surface_type == surface_param.surface_land) 
+        if (surface_type == surface.surface_land) 
         {
             h0_m = h / z0_m;
-            u_dyn0 = U * model_param.kappa / log(h0_m);
+            u_dyn0 = U * model.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);
+        Re = u_dyn0 * z0_m / phys.nu_air;
+        get_thermal_roughness(z0_t, B, z0_m, Re, surface, surface_type);
 
         h0_t = h / z0_t;
-        Rib = (phys_constants.g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / (U*U);
+        Rib = (phys.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);
+                            model);
 
 
         if (Rib > 0.0) 
         {
-            Rib = std::min(Rib, model_param.Rib_max);
-            get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model_param);
+            Rib = std::min(Rib, model.Rib_max);
+            get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model);
 
-            fval = model_param.beta_m * zeta;
+            fval = model.beta_m * zeta;
             phi_m = 1.0 + fval;
-            phi_h = 1.0/model_param.Pr_t_0_inv + fval;
+            phi_h = 1.0/model.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.maxiters_convection);
+            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, numerics.maxiters_convection);
 
             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);
+            phi_h = fval / (model.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);
+            get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B, model);
         
             phi_m = 1.0;
-            phi_h = 1.0 / model_param.Pr_t_0_inv;
+            phi_h = 1.0 / model.Pr_t_0_inv;
         }
         else
         {
-            get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model_param, numerics.maxiters_convection);
+            get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model, numerics.maxiters_convection);
             
-            phi_m = pow(1.0 - model_param.alpha_m * zeta, -0.25);
-            phi_h = 1.0 / (model_param.Pr_t_0_inv * sqrt(1.0 - model_param.alpha_h_fix * zeta));
+            phi_m = pow(1.0 - model.alpha_m * zeta, -0.25);
+            phi_h = 1.0 / (model.Pr_t_0_inv * sqrt(1.0 - model.alpha_h_fix * zeta));
         }
 
-        Cm = model_param.kappa / psi_m;
-        Ct = model_param.kappa / psi_h;
+        Cm = model.kappa / psi_m;
+        Ct = model.kappa / psi_h;
 
-        Km = model_param.kappa * Cm * U * h / phi_m;
+        Km = model.kappa * Cm * U * h / phi_m;
         Pr_t_inv = phi_m / phi_h;
 
         sfx.zeta[step]         = zeta;
diff --git a/srcCXX/sfx_sheba.cpp b/srcCXX/sfx_sheba.cpp
index abd48c71e3c01b939e5b517525cb4963b8a441c4..a6b74d7854ad72a2d5473097f8c3c5d4e7d3897a 100644
--- a/srcCXX/sfx_sheba.cpp
+++ b/srcCXX/sfx_sheba.cpp
@@ -13,15 +13,15 @@
 template<typename T, MemType memIn, MemType memOut, MemType RunMem >
 FluxShebaBase<T, memIn, memOut, RunMem>::FluxShebaBase(sfxDataVecTypeC* sfx_in,
                 meteoDataVecTypeC* meteo_in,
-                const sfx_sheba_param model_param_in, 
+                const sfx_sheba_param_C model_param_in, 
                 const sfx_surface_param surface_param_in,
-                const sfx_sheba_numericsTypeC numerics_in,
+                const sfx_sheba_numericsType_C numerics_in,
                 const sfx_phys_constants phys_constants_in,
                 const int grid_size_in) : ModelBase<T, memIn, memOut, RunMem>(sfx_in, meteo_in, grid_size_in)
 {
-    surface_param = surface_param_in;
-    phys_constants = phys_constants_in;
-    model_param = model_param_in;
+    surface = surface_param_in;
+    phys = phys_constants_in;
+    model = model_param_in;
     numerics = numerics_in;
 }
 
@@ -46,35 +46,35 @@ void FluxSheba<T, memIn, memOut, MemType::CPU>::compute_flux()
         h = meteo.h[step];
         z0_m = meteo.z0_m[step];
 
-        surface_type = z0_m < 0.0 ? surface_param.surface_ocean : surface_param.surface_land;
+        surface_type = z0_m < 0.0 ? surface.surface_ocean : surface.surface_land;
 
-        if (surface_type == surface_param.surface_ocean) 
+        if (surface_type == surface.surface_ocean) 
         {
-            get_charnock_roughness(z0_m, u_dyn0, U, h, surface_param, numerics.maxiters_charnock);
+            get_charnock_roughness(z0_m, u_dyn0, U, h, surface, numerics.maxiters_charnock);
             h0_m = h / z0_m;
         }
-        if (surface_type == surface_param.surface_land) 
+        if (surface_type == surface.surface_land) 
         {
             h0_m = h / z0_m;
-            u_dyn0 = U * model_param.kappa / log(h0_m);
+            u_dyn0 = U * model.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);
+        Re = u_dyn0 * z0_m / phys.nu_air;
+        get_thermal_roughness(z0_t, B, z0_m, Re, surface, surface_type);
 
         // --- define relative height [thermal]
         h0_t = h / z0_t;
 
         // --- define Ri-bulk
-        Rib = (phys_constants.g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / (U*U);
+        Rib = (phys.g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / (U*U);
 
         // --- get the fluxes
         // ----------------------------------------------------------------------------
         get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, 
-            U, Tsemi, dT, dQ, h, z0_m, z0_t, (phys_constants.g / Tsemi), model_param, 10);
+            U, Tsemi, dT, dQ, h, z0_m, z0_t, (phys.g / Tsemi), model, 10);
         // ----------------------------------------------------------------------------
 
-        get_phi(phi_m, phi_h, zeta, model_param);
+        get_phi(phi_m, phi_h, zeta, model);
         // ----------------------------------------------------------------------------
 
         // --- define transfer coeff. (momentum) & (heat)
@@ -86,7 +86,7 @@ void FluxSheba<T, memIn, memOut, MemType::CPU>::compute_flux()
             Ct = Tdyn / dT;
 
         // --- define eddy viscosity & inverse Prandtl number
-        Km = model_param.kappa * Cm * U * h / phi_m;
+        Km = model.kappa * Cm * U * h / phi_m;
         Pr_t_inv = phi_m / phi_h;
 
         sfx.zeta[step]         = zeta;
diff --git a/srcF/sfx_data.f90 b/srcF/sfx_data.f90
index e4386ecf42ba25d62cd308c08042d6fc34d0089f..9fb30d277d96530758977a44caf2b198b637207a 100644
--- a/srcF/sfx_data.f90
+++ b/srcF/sfx_data.f90
@@ -137,43 +137,6 @@ module sfx_data
         real(C_FLOAT)  :: g;
         real(C_FLOAT)  :: nu_air;
     end type
-
-    ! ESM strcuctures
-    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_esm_numericsTypeC 
-        integer(C_INT) :: maxiters_convection           
-        integer(C_INT) :: maxiters_charnock 
-    end type
-
-    ! SHEBA strcuctures
-    type, BIND(C), public :: sfx_sheba_param 
-        real(C_FLOAT) :: kappa           
-        real(C_FLOAT) :: Pr_t_0_inv
-
-        real(C_FLOAT) :: alpha_m           
-        real(C_FLOAT) :: alpha_h
-        real(C_FLOAT) :: a_m
-        real(C_FLOAT) :: b_m
-        real(C_FLOAT) :: a_h
-        real(C_FLOAT) :: b_h
-        real(C_FLOAT) :: c_h
-    end type
-
-    type, BIND(C), public :: sfx_sheba_numericsTypeC 
-        integer(C_INT) :: maxiters_charnock 
-    end type
 #endif
     ! --------------------------------------------------------------------------------
 
diff --git a/srcF/sfx_esm.f90 b/srcF/sfx_esm.f90
index 9b63b4f4fa63e4422ac3cbf087b504f2f29e67ec..e5312ddee3e4a83c3de898a8762b7eb881d8536f 100644
--- a/srcF/sfx_esm.f90
+++ b/srcF/sfx_esm.f90
@@ -12,7 +12,7 @@ module sfx_esm
     use sfx_surface
     use sfx_esm_param
 #if defined(INCLUDE_CXX)
-    use iso_c_binding, only: C_LOC, C_PTR
+    use iso_c_binding, only: C_LOC, C_PTR, C_INT, C_FLOAT
     use C_FUNC
 #endif
     ! --------------------------------------------------------------------------------
@@ -36,12 +36,48 @@ module sfx_esm
     end type
     ! --------------------------------------------------------------------------------
 
+#if defined(INCLUDE_CXX)
+    type, BIND(C), public :: sfx_esm_param_C 
+        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_esm_numericsType_C 
+        integer(C_INT) :: maxiters_convection           
+        integer(C_INT) :: maxiters_charnock 
+    end type
+
+    INTERFACE
+        SUBROUTINE get_surface_fluxes_esm(sfx, meteo, model_param, surface_param, numerics, constants, grid_size) BIND(C)
+            use sfx_data
+            use, intrinsic :: ISO_C_BINDING, ONLY: C_INT, C_PTR
+            Import :: sfx_esm_param_C, sfx_esm_numericsType_C
+            implicit none
+            integer(C_INT) :: grid_size
+            type(C_PTR), value :: sfx
+            type(C_PTR), value :: meteo
+            type(sfx_esm_param_C) :: model_param
+            type(sfx_surface_param) :: surface_param
+            type(sfx_esm_numericsType_C) :: numerics
+            type(sfx_phys_constants) :: constants
+        END SUBROUTINE get_surface_fluxes_esm
+    END INTERFACE
+#endif 
+
 contains
 
     ! --------------------------------------------------------------------------------
 #if defined(INCLUDE_CXX)
     subroutine set_c_struct_sfx_esm_param_values(sfx_model_param)
-        type (sfx_esm_param), intent(inout) :: sfx_model_param
+        type (sfx_esm_param_C), intent(inout) :: sfx_model_param
         sfx_model_param%kappa = kappa
         sfx_model_param%Pr_t_0_inv = Pr_t_0_inv
         sfx_model_param%Pr_t_inf_inv = Pr_t_inf_inv
@@ -75,9 +111,9 @@ contains
         type (meteoDataVecTypeC), target :: meteo_c         !< meteorological data (input)
         type (sfxDataVecTypeC), target :: sfx_c             !< surface flux data (output)
         type(C_PTR) :: meteo_c_ptr, sfx_c_ptr
-        type (sfx_esm_param) :: model_param
+        type (sfx_esm_param_C) :: model_param
         type (sfx_surface_param) :: surface_param
-        type (sfx_esm_numericsTypeC) :: numerics_c
+        type (sfx_esm_numericsType_C) :: numerics_c
         type (sfx_phys_constants) :: phys_constants
 
         numerics_c%maxiters_convection = numerics%maxiters_convection
diff --git a/srcF/sfx_fc_wrapper.F90 b/srcF/sfx_fc_wrapper.F90
index 0e84bb3f56d8e6b0134a4e40ea53118d72ff13bc..6e2fc0238c09f403eede892ba9d2e517b4347313 100644
--- a/srcF/sfx_fc_wrapper.F90
+++ b/srcF/sfx_fc_wrapper.F90
@@ -1,35 +1,6 @@
 module C_FUNC
 
 #if defined(INCLUDE_CXX)
-
-    INTERFACE
-        SUBROUTINE get_surface_fluxes_esm(sfx, meteo, model_param, surface_param, numerics, constants, grid_size) BIND(C)
-            use sfx_data
-            USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_INT, C_PTR
-            IMPLICIT NONE
-            INTEGER(C_INT) :: grid_size
-            type(C_PTR), value :: sfx
-            type(C_PTR), value :: meteo
-            type(sfx_esm_param) :: model_param
-            type(sfx_surface_param) :: surface_param
-            type(sfx_esm_numericsTypeC) :: numerics
-            type(sfx_phys_constants) :: constants
-        END SUBROUTINE get_surface_fluxes_esm
-
-        SUBROUTINE get_surface_fluxes_sheba(sfx, meteo, model_param, surface_param, numerics, constants, grid_size) BIND(C)
-            use sfx_data
-            USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_INT, C_PTR
-            IMPLICIT NONE
-            INTEGER(C_INT) :: grid_size
-            type(C_PTR), value :: sfx
-            type(C_PTR), value :: meteo
-            type(sfx_sheba_param) :: model_param
-            type(sfx_surface_param) :: surface_param
-            type(sfx_sheba_numericsTypeC) :: numerics
-            type(sfx_phys_constants) :: constants
-        END SUBROUTINE get_surface_fluxes_sheba
-    END INTERFACE
-
     contains
 
     subroutine set_c_struct_sfx_phys_constants_values(constants)
@@ -41,6 +12,7 @@ module C_FUNC
         constants%g = g
         constants%nu_air = nu_air
     end subroutine set_c_struct_sfx_phys_constants_values
+    
 #endif
 end module C_FUNC
 
diff --git a/srcF/sfx_sheba.f90 b/srcF/sfx_sheba.f90
index 7fe86c55da117cd88c6c97e28d8bce38c27cc3f8..7d7cde2bdd513cdc346363de2176874c7266b3d8 100644
--- a/srcF/sfx_sheba.f90
+++ b/srcF/sfx_sheba.f90
@@ -13,7 +13,7 @@ module sfx_sheba
     use sfx_sheba_param
 
 #if defined(INCLUDE_CXX)
-    use iso_c_binding, only: C_LOC, C_PTR
+    use iso_c_binding, only: C_LOC, C_PTR, C_INT, C_FLOAT
     use C_FUNC
 #endif
     ! --------------------------------------------------------------------------------
@@ -36,11 +36,46 @@ module sfx_sheba
     end type
     ! --------------------------------------------------------------------------------
 
+#if defined(INCLUDE_CXX)
+    type, BIND(C), public :: sfx_sheba_param_C
+        real(C_FLOAT) :: kappa           
+        real(C_FLOAT) :: Pr_t_0_inv
+
+        real(C_FLOAT) :: alpha_m           
+        real(C_FLOAT) :: alpha_h
+        real(C_FLOAT) :: a_m
+        real(C_FLOAT) :: b_m
+        real(C_FLOAT) :: a_h
+        real(C_FLOAT) :: b_h
+        real(C_FLOAT) :: c_h
+    end type
+
+    type, BIND(C), public :: sfx_sheba_numericsType_C 
+        integer(C_INT) :: maxiters_charnock 
+    end type
+
+    INTERFACE
+        SUBROUTINE get_surface_fluxes_sheba(sfx, meteo, model_param, surface_param, numerics, constants, grid_size) BIND(C)
+            use sfx_data
+            use, intrinsic :: ISO_C_BINDING, ONLY: C_INT, C_PTR
+            Import :: sfx_sheba_param_C, sfx_sheba_numericsType_C
+            implicit none
+            INTEGER(C_INT) :: grid_size
+            type(C_PTR), value :: sfx
+            type(C_PTR), value :: meteo
+            type(sfx_sheba_param_C) :: model_param
+            type(sfx_surface_param) :: surface_param
+            type(sfx_sheba_numericsType_C) :: numerics
+            type(sfx_phys_constants) :: constants
+        END SUBROUTINE get_surface_fluxes_sheba
+    END INTERFACE
+#endif 
+
 contains
 
 #if defined(INCLUDE_CXX)
     subroutine set_c_struct_sfx_sheba_param_values(sfx_model_param)
-        type (sfx_sheba_param), intent(inout) :: sfx_model_param
+        type (sfx_sheba_param_C), intent(inout) :: sfx_model_param
         sfx_model_param%kappa = kappa
         sfx_model_param%Pr_t_0_inv = Pr_t_0_inv
 
@@ -75,9 +110,9 @@ contains
         type (meteoDataVecTypeC), pointer :: meteo_c         !< meteorological data (input)
         type (sfxDataVecTypeC), pointer :: sfx_c             !< surface flux data (output)
         type(C_PTR) :: meteo_c_ptr, sfx_c_ptr
-        type (sfx_sheba_param) :: model_param
+        type (sfx_sheba_param_C) :: model_param
         type (sfx_surface_param) :: surface_param
-        type (sfx_sheba_numericsTypeC) :: numerics_c
+        type (sfx_sheba_numericsType_C) :: numerics_c
         type (sfx_phys_constants) :: phys_constants
 
         numerics_c%maxiters_charnock = numerics%maxiters_charnock