diff --git a/CMakeLists.txt b/CMakeLists.txt
index 21fe045e2364bc125e8f97805e18a47c76797df3..f38adb2ab6aac4454be4d06ce13764b9f0d132fb 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -45,6 +45,7 @@ if(INCLUDE_CXX)
 endif(INCLUDE_CXX)
 
 set(SOURCES_F 
+    srcF/sfx_io.f90
     srcF/sfx_data.f90
     srcF/sfx_common.f90
     srcF/sfx_esm.f90 
@@ -54,7 +55,7 @@ set(SOURCES_F
     srcF/sfx_main.f90
     srcF/sfx_phys_const.f90
     srcF/sfx_surface.f90
-    srcF/FCWrapper.F90
+    srcF/sfx_fc_wrapper.F90
 )
 
 set(HEADERS_F
@@ -63,45 +64,45 @@ set(HEADERS_F
 
 if(INCLUDE_CXX)
     set(SOURCES_C 
-        srcC/SubFunctionsWrapper.c
+        srcC/sfx_call_cxx.c
     )
 
     set(SOURCES_CXX 
-            srcCXX/Flux.cpp
-            srcCXX/FluxComputeFunc.cpp
-            srcCXX/SubFunctions.cpp
+            srcCXX/sfx_flux.cpp
+            srcCXX/sfx_flux_compute_func.cpp
+            srcCXX/sfx_call_class_func.cpp
     )
     set(HEADERS_CXX 
-            includeCXX/Flux.h
-            includeCXX/FluxComputeFunc.h
-            includeCXX/SubFunctions.h
+            includeCXX/sfx_flux.h
+            includeCXX/sfx_flux_compute_func.h
+            includeCXX/sfx_call_class_func.h
         )
 endif(INCLUDE_CXX)
 
 if(INCLUDE_CUDA)
     set(SOURCES_CU 
-        srcCU/FluxComputeFunc.cu
+        srcCU/sfx_flux_compute_func.cu
     )
     set(HEADERS_CU
-        includeCU/FluxComputeFunc.cuh
+        includeCU/sfx_flux_compute_func.cuh
     )
 endif(INCLUDE_CUDA)
 
 if(INCLUDE_CXX OR INCLUDE_CUDA)
     set(MEMPROC_SOURCES_CXX
-        srcCXX/MemoryProcessing.cpp
+        srcCXX/sfx_memory_processing.cpp
     )
     set(MEMPROC_HEADERS_CXX
-        includeCXX/MemoryProcessing.h
-        includeCXX/TemplateParameters.h
+        includeCXX/sfx_memory_processing.h
+        includeCXX/sfx_template_parameters.h
     )
 
     if(INCLUDE_CUDA)
         set(MEMPROC_SOURCES_CU
-            srcCU/MemoryProcessing.cu
+            srcCU/sfx_memory_processing.cu
         )
         set(MEMPROC_HEADERS_CU
-            includeCU/MemoryProcessing.cuh
+            includeCU/sfx_memory_processing.cuh
         )
     endif(INCLUDE_CUDA)
 endif(INCLUDE_CXX OR INCLUDE_CUDA)
diff --git a/includeCU/FluxComputeFunc.cuh b/includeCU/sfx_flux_compute_func.cuh
similarity index 69%
rename from includeCU/FluxComputeFunc.cuh
rename to includeCU/sfx_flux_compute_func.cuh
index 2c2623e267c6fa9f28678bc1fbfb9517df505688..2628d5d4a44a1e5baedddc7a1913683d9af4c9c6 100644
--- a/includeCU/FluxComputeFunc.cuh
+++ b/includeCU/sfx_flux_compute_func.cuh
@@ -1,8 +1,8 @@
 #pragma once 
 
 template<typename T>
-void compute_flux_gpu(const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_, const T *h_, const T *in_z0_m_,
-    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_,
+void compute_flux_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, 
diff --git a/includeCU/MemoryProcessing.cuh b/includeCU/sfx_memory_processing.cuh
similarity index 91%
rename from includeCU/MemoryProcessing.cuh
rename to includeCU/sfx_memory_processing.cuh
index 49809a2d2502a0f9bccdf14e43e3347b874ce533..c68c858db0764633b37d9cfce0ee0f92bd1d1d69 100644
--- a/includeCU/MemoryProcessing.cuh
+++ b/includeCU/sfx_memory_processing.cuh
@@ -1,5 +1,5 @@
 #pragma once
-#include "../includeCXX/TemplateParameters.h"
+#include "../includeCXX/sfx_template_parameters.h"
 #include <cstddef>
 
 namespace memproc
diff --git a/includeCXX/SubFunctions.h b/includeCXX/sfx_call_class_func.h
similarity index 71%
rename from includeCXX/SubFunctions.h
rename to includeCXX/sfx_call_class_func.h
index e041131b9881dc792faf14210ef4f05f2df1dc44..08c9164ba17a119c28420fa4b35284c07a7cdaff 100644
--- a/includeCXX/SubFunctions.h
+++ b/includeCXX/sfx_call_class_func.h
@@ -4,8 +4,8 @@
 extern "C" {
 #endif
 
-    void surf_flux_CXX (float *U_, float *dT_, float *Tsemi_, float *dQ_, float *h_, float *in_z0_m_,
-    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_,
+    void surf_flux_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, 
diff --git a/includeCXX/Flux.h b/includeCXX/sfx_flux.h
similarity index 78%
rename from includeCXX/Flux.h
rename to includeCXX/sfx_flux.h
index cc7f6c2bb7fd77abcb2e2b21b094f6717a6d755f..8dcb5f09a6e8d239c409701a181bd09e69ec89fe 100644
--- a/includeCXX/Flux.h
+++ b/includeCXX/sfx_flux.h
@@ -1,5 +1,5 @@
 #pragma once
-#include "TemplateParameters.h"
+#include "sfx_template_parameters.h"
 #include <cstddef>
 
 template<typename T, MemType RunMem, MemType memIn>
@@ -31,16 +31,13 @@ public:
     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);
-    void set_data( T *U_, T *dT_, T *Tsemi_, T *dQ_, T *h_, T *in_z0_m_,
-    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_);
     ~Flux();
 
-    void compute_flux(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_, 
+    void compute_flux(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 RunMem, MemType memIn>
-class Compute_Flux
-{
-    static Flux<T, RunMem,memIn> F;
+private:
+    void set_data( T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_, T *Rib_conv_lim_, T *Cm_, T *Ct_, T *Km_, T *Pr_t_inv_,
+    T *U_, T *dT_, T *Tsemi_, T *dQ_, T *h_, T *in_z0_m_);
 };
\ No newline at end of file
diff --git a/includeCXX/FluxComputeFunc.h b/includeCXX/sfx_flux_compute_func.h
similarity index 69%
rename from includeCXX/FluxComputeFunc.h
rename to includeCXX/sfx_flux_compute_func.h
index d74553d8f78b8a4a0d70749abcfc40cdf52ce15a..2cfa0f408af3092c1da29b99e7891f3b8428205c 100644
--- a/includeCXX/FluxComputeFunc.h
+++ b/includeCXX/sfx_flux_compute_func.h
@@ -1,8 +1,8 @@
 #pragma once
 
 template<typename T>
-void compute_flux_cpu(const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_, const T *h_, const T *in_z0_m_,
-    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_,
+void compute_flux_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, 
diff --git a/includeCXX/MemoryProcessing.h b/includeCXX/sfx_memory_processing.h
similarity index 93%
rename from includeCXX/MemoryProcessing.h
rename to includeCXX/sfx_memory_processing.h
index b3db1fece0388f6a471dfbd833fed4e62bd7be8b..2de2fcd729344226dbd8dec326cabe3d1f10caa6 100644
--- a/includeCXX/MemoryProcessing.h
+++ b/includeCXX/sfx_memory_processing.h
@@ -1,5 +1,5 @@
 #pragma once
-#include "TemplateParameters.h"
+#include "sfx_template_parameters.h"
 #include <cstddef>
 
 namespace memproc
diff --git a/includeCXX/TemplateParameters.h b/includeCXX/sfx_template_parameters.h
similarity index 100%
rename from includeCXX/TemplateParameters.h
rename to includeCXX/sfx_template_parameters.h
diff --git a/srcC/SubFunctionsWrapper.c b/srcC/sfx_call_cxx.c
similarity index 69%
rename from srcC/SubFunctionsWrapper.c
rename to srcC/sfx_call_cxx.c
index 91372f764c3fba510487feef8189044181613fed..11ec1f29db53d9705f3e820741aa8c08554f4f65 100644
--- a/srcC/SubFunctionsWrapper.c
+++ b/srcC/sfx_call_cxx.c
@@ -1,9 +1,9 @@
 #include <stdlib.h>
 #include <stdio.h>
-#include "../includeCXX/SubFunctions.h"
+#include "../includeCXX/sfx_call_class_func.h"
 // -------------------------------------------------------------------------- //
-void surf_flux (float *U_, float *dT_, float *Tsemi_, float *dQ_, float *h_, float *in_z0_m_,
-    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_,
+void surf_flux (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, 
@@ -14,8 +14,8 @@ void surf_flux (float *U_, float *dT_, float *Tsemi_, float *dQ_, float *h_, flo
     const int *maxiters_charnock, const int *maxiters_convection, 
     const int *grid_size)
 {
-	surf_flux_CXX ( U_, dT_, Tsemi_, dQ_, h_, in_z0_m_,
-  zeta_, Rib_, Re_, B_, z0_m_, z0_t_, Rib_conv_lim_, Cm_, Ct_, Km_, Pr_t_inv_,
+	surf_flux_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, 
diff --git a/srcCU/FluxComputeFunc.cu b/srcCU/sfx_flux_compute_func.cu
similarity index 68%
rename from srcCU/FluxComputeFunc.cu
rename to srcCU/sfx_flux_compute_func.cu
index b631fec080f5191b2738ba4e5452cd7b82b8af69..ff912b9a4c628ee9e15769e8e6eee1e1b57903c0 100644
--- a/srcCU/FluxComputeFunc.cu
+++ b/srcCU/sfx_flux_compute_func.cu
@@ -1,12 +1,12 @@
 #include <cmath>
 #include <iostream>
-#include "../includeCU/FluxComputeFunc.cuh"
+#include "../includeCU/sfx_flux_compute_func.cuh"
 
 template<typename T>
-__device__ void get_charnock_roughness(const T h, const T U,
+__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,
-    T &z0_m, T &u_dyn0, 
+    const T h_charnock, const T c1_charnock, const T c2_charnock, 
     const int maxiters)
 {
     T Uc, a, b, c, c_min, f;
@@ -35,22 +35,22 @@ __device__ void get_charnock_roughness(const T h, const T U,
     u_dyn0 = Uc / c;
 }
 
-template __device__ void get_charnock_roughness(const float h, const float U,
+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,
-    float &z0_m, float &u_dyn0, 
+    const float h_charnock, const float c1_charnock, const float c2_charnock, 
     const int maxiters);
-template __device__ void get_charnock_roughness(const double h, const double U,
+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,
-    double &z0_m, double &u_dyn0, 
     const int maxiters);
 
 template<typename T>
-__device__ void get_convection_lim(const T h0_m, const T h0_t, const T B,
+__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 &zeta_lim, T &Rib_lim, T &f_m_lim, T &f_h_lim)
+    const T alpha_h, const T alpha_m, const T alpha_h_fix)
 {
     T psi_m, psi_h, f_m, f_h, c;
 
@@ -73,19 +73,19 @@ __device__ void get_convection_lim(const T h0_m, const T h0_t, const T B,
     Rib_lim = zeta_lim * psi_h / (psi_m * psi_m);
 }
 
-template __device__ void get_convection_lim(const float h0_m, const float h0_t, const float B,
+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,
-    float &zeta_lim, float &Rib_lim, float &f_m_lim, float &f_h_lim);
-template __device__ void get_convection_lim(const double h0_m, const double h0_t, const double B,
+    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,
-    double &zeta_lim, double &Rib_lim, double &f_m_lim, double &f_h_lim);
+    const double alpha_h, const double alpha_m, const double alpha_h_fix);
 
 template<typename T>
-void __device__ get_psi_stable(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 &psi_m, T &psi_h, T &zeta)
+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;
     
@@ -102,19 +102,19 @@ void __device__ get_psi_stable(const T Rib, const T h0_m, const T h0_t, const T
     psi_h = (psi0_m + B) / Pr_t_0_inv + phi;
 }
 
-template __device__ void get_psi_stable(const float Rib, const float h0_m, const float h0_t, const float B,
-    const float Pr_t_0_inv, const float beta_m,
-    float &psi_m, float &psi_h, float &zeta);
-template __device__ void get_psi_stable(const double Rib, const double h0_m, const double h0_t, const double B,
-    const double Pr_t_0_inv, const double beta_m,
-    double &psi_m, double &psi_h, double &zeta);
+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(const T Rib, const T h0_m, const T h0_t, const T B, 
+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,
-    T &psi_m, T &psi_h, T &zeta, 
     const int maxiters)
 {
     T zeta0_m, zeta0_h, f0_m, f0_h, p_m, p_h, a_m, a_h, c_lim, f;
@@ -150,23 +150,23 @@ void __device__ get_psi_convection(const T Rib, const T h0_m, const T h0_t, cons
     }
 }
 
-template __device__ void get_psi_convection(const float Rib, const float h0_m, const float h0_t, const float B, 
+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,
-    float &psi_m, float &psi_h, float &zeta, 
     const int maxiters);
-template __device__ void get_psi_convection(const double Rib, const double h0_m, const double h0_t, const double B, 
+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,
-    double &psi_m, double &psi_h, double &zeta, 
     const int maxiters);
 
 template<typename T>
-void __device__ get_psi_neutral(const T h0_m, const T h0_t, const T B,
-    const T Pr_t_0_inv,
-    T &psi_m, T &psi_h, T &zeta)
+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);
@@ -175,18 +175,18 @@ void __device__ get_psi_neutral(const T h0_m, const T h0_t, const T B,
         psi_h = psi_m / Pr_t_0_inv;
 }
 
-template __device__ void get_psi_neutral(const float h0_m, const float h0_t, const float B,
-    const float Pr_t_0_inv,
-    float &psi_m, float &psi_h, float &zeta);
-template __device__ void get_psi_neutral(const double h0_m, const double h0_t, const double B,
-    const double Pr_t_0_inv,
-    double &psi_m, double &psi_h, double &zeta);
+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(const T Rib, const T h0_m, const T h0_t, const T B, 
+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,
-    T &psi_m, T &psi_h, T &zeta,
     const int maxiters)
 {
     T zeta0_m, zeta0_h, f0_m, f0_h, f_m, f_h;
@@ -225,20 +225,20 @@ void __device__ get_psi_semi_convection(const T Rib, const T h0_m, const T h0_t,
     }
 }
 
-template __device__ void get_psi_semi_convection(const float Rib, const float h0_m, const float h0_t, const float B, 
+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,
-    float &psi_m, float &psi_h, float &zeta,
     const int maxiters);
-template __device__ void get_psi_semi_convection(const double Rib, const double h0_m, const double h0_t, const double B, 
+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,
-    double &psi_m, double &psi_h, double &zeta,
     const int maxiters);
 
 template<typename T>
-__global__ void compute_flux(const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_, const T *h_, const T *in_z0_m_,
-    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_,
+__global__ void compute_flux(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, 
@@ -273,7 +273,7 @@ __global__ void compute_flux(const T *U_, const T *dT_, const T *Tsemi_, const T
 
         if (surface_type == 0)
         {
-            get_charnock_roughness(h, U, kappa, h_charnock, c1_charnock, c2_charnock, z0_m, u_dyn0, maxiters_charnock);
+            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) 
@@ -295,13 +295,13 @@ __global__ void compute_flux(const T *U_, const T *dT_, const T *Tsemi_, const T
         h0_t = h / z0_t;
         Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / (U*U);
 
-        get_convection_lim(h0_m, h0_t, B, Pr_t_inf_inv, Pr_t_0_inv, alpha_h, alpha_m, alpha_h_fix, zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim);
+        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(Rib, h0_m, h0_t, B, Pr_t_0_inv, beta_m, psi_m, psi_h, zeta);
+            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;
@@ -310,7 +310,7 @@ __global__ void compute_flux(const T *U_, const T *dT_, const T *Tsemi_, const T
 
         else if (Rib < Rib_conv_lim) 
         {
-            get_psi_convection(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, psi_m, psi_h, zeta, 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, 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;
@@ -318,14 +318,14 @@ __global__ void compute_flux(const T *U_, const T *dT_, const T *Tsemi_, const T
         }
         else if (Rib > -0.001) 
         {
-            get_psi_neutral(h0_m, h0_t, B, Pr_t_0_inv, psi_m, psi_h, zeta);
+            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(Rib, h0_m, h0_t, B, Pr_t_0_inv, alpha_m, alpha_h_fix, psi_m, psi_h, zeta, maxiters_convection);
+            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));
@@ -351,8 +351,8 @@ __global__ void compute_flux(const T *U_, const T *dT_, const T *Tsemi_, const T
     }
 }
 
-template __global__ void compute_flux(const float *U, const float *dt, const float *T_semi, const float *dq, const float *H, const float *in_z0_m,
-    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_,
+template __global__ void compute_flux(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, 
@@ -362,8 +362,8 @@ template __global__ void compute_flux(const float *U, const float *dt, const flo
     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 compute_flux(const double *U, const double *dt, const double *T_semi, const double *dq, const double *H, const double *in_z0_m, 
-    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_,
+template __global__ void compute_flux(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, 
@@ -375,8 +375,8 @@ template __global__ void compute_flux(const double *U, const double *dt, const d
     const int grid_size);
 
 template<typename T>
-void compute_flux_gpu(const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_, const T *h_, const T *in_z0_m_,
-    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_,
+void compute_flux_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, 
@@ -391,8 +391,8 @@ void compute_flux_gpu(const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_,
     dim3 cuBlock = dim3(1024, 1, 1);
 	dim3 cuGrid = dim3(BlockCount, 1, 1);
 
-    compute_flux<<<cuGrid, cuBlock>>>(U_, dT_, Tsemi_, dQ_, h_, in_z0_m_,
-    zeta_, Rib_, Re_, B_, z0_m_, z0_t_, Rib_conv_lim_, Cm_, Ct_, Km_, Pr_t_inv_,
+    compute_flux<<<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, 
@@ -404,8 +404,8 @@ void compute_flux_gpu(const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_,
     grid_size);
 }
 
-template void compute_flux_gpu(const float *U, const float *dt, const float *T_semi, const float *dq, const float *H, const float *in_z0_m,
-    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_,
+template void compute_flux_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, 
@@ -415,8 +415,8 @@ template void compute_flux_gpu(const float *U, const float *dt, const float *T_s
     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_gpu(const double *U, const double *dt, const double *T_semi, const double *dq, const double *H, const double *in_z0_m, 
-    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_,
+template void compute_flux_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, 
diff --git a/srcCU/MemoryProcessing.cu b/srcCU/sfx_memory_processing.cu
similarity index 97%
rename from srcCU/MemoryProcessing.cu
rename to srcCU/sfx_memory_processing.cu
index e00a089256180780d024901afb8325564cff8ec4..7dcc6cb5aedd07261cbe7d38ae1dcb5630859e73 100644
--- a/srcCU/MemoryProcessing.cu
+++ b/srcCU/sfx_memory_processing.cu
@@ -1,4 +1,4 @@
-#include "../includeCU/MemoryProcessing.cuh"
+#include "../includeCU/sfx_memory_processing.cuh"
 #include <cuda.h>
 #include <cuda_runtime_api.h>
 
diff --git a/srcCXX/SubFunctions.cpp b/srcCXX/sfx_call_class_func.cpp
similarity index 68%
rename from srcCXX/SubFunctions.cpp
rename to srcCXX/sfx_call_class_func.cpp
index 010f45af503e0dffa05f15d05f50a09a807441e0..2243e0d8063a473c9044311ccd5ba19c30758516 100644
--- a/srcCXX/SubFunctions.cpp
+++ b/srcCXX/sfx_call_class_func.cpp
@@ -1,7 +1,7 @@
 #include <stdlib.h>
 #include <stdio.h>
-#include "../includeCXX/SubFunctions.h"
-#include "../includeCXX/Flux.h"
+#include "../includeCXX/sfx_call_class_func.h"
+#include "../includeCXX/sfx_flux.h"
 #include <vector>
 
 #ifdef INCLUDE_CUDA
@@ -11,8 +11,8 @@
 #endif
 
 // -------------------------------------------------------------------------- //
-void surf_flux_CXX (float *U_, float *dT_, float *Tsemi_, float *dQ_, float *h_, float *in_z0_m_,
-    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_,
+void surf_flux_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, 
@@ -31,8 +31,7 @@ void surf_flux_CXX (float *U_, float *dT_, float *Tsemi_, float *dQ_, float *h_,
     gamma_c, Re_visc_min,
     Pr_m, nu_air, g);
 
-    F.set_data(U_, dT_, Tsemi_, dQ_, h_, in_z0_m_,
-    zeta_, Rib_, Re_, B_, z0_m_, z0_t_, Rib_conv_lim_, Cm_, Ct_, Km_, Pr_t_inv_);
-
-    F.compute_flux(zeta_, Rib_, Re_, B_, z0_m_, z0_t_, Rib_conv_lim_, Cm_, Ct_, Km_, Pr_t_inv_, maxiters_charnock, maxiters_convection);
+    F.compute_flux(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);
 }
\ No newline at end of file
diff --git a/srcCXX/Flux.cpp b/srcCXX/sfx_flux.cpp
similarity index 85%
rename from srcCXX/Flux.cpp
rename to srcCXX/sfx_flux.cpp
index 08508cc8115da2e0ca3cff54b6e5a0bd2e6c26b0..8f70e1503a7d2e921ccc8f990121120fde053081 100644
--- a/srcCXX/Flux.cpp
+++ b/srcCXX/sfx_flux.cpp
@@ -1,13 +1,13 @@
 #include <iostream>
 
-#include "../includeCXX/Flux.h"
-#include "../includeCXX/FluxComputeFunc.h"
+#include "../includeCXX/sfx_flux.h"
+#include "../includeCXX/sfx_flux_compute_func.h"
 #ifdef INCLUDE_CUDA
-    #include "../includeCU/FluxComputeFunc.cuh"
-    #include "../includeCU/MemoryProcessing.cuh"
+    #include "../includeCU/sfx_flux_compute_func.cuh"
+    #include "../includeCU/sfx_memory_processing.cuh"
 #endif
 
-#include "../includeCXX/MemoryProcessing.h"
+#include "../includeCXX/sfx_memory_processing.h"
 
 template<typename T, MemType RunMem, MemType memIn>
 Flux<T, RunMem, memIn>::Flux()
@@ -134,8 +134,8 @@ Flux<T, RunMem, memIn>::~Flux()
 }
 
 template<typename T, MemType RunMem, MemType memIn>
-void Flux<T, RunMem, memIn>::set_data(T *U_, T *dT_, T *Tsemi_, T *dQ_, T *h_, T *in_z0_m_,
-    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_)
+void Flux<T, RunMem, memIn>::set_data(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_, T *Rib_conv_lim_, T *Cm_, T *Ct_, T *Km_, T *Pr_t_inv_,
+    T *U_, T *dT_, T *Tsemi_, T *dQ_, T *h_, T *in_z0_m_)
 {
     if(RunMem == memIn)
     {
@@ -172,11 +172,15 @@ void Flux<T, RunMem, memIn>::set_data(T *U_, T *dT_, T *Tsemi_, T *dQ_, T *h_, T
 }
 
 template<typename T, MemType RunMem, MemType memIn>
-void Flux<T, RunMem, memIn>::compute_flux(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_, 
+void Flux<T, RunMem, memIn>::compute_flux(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)
 {
-    if(RunMem == MemType::CPU) compute_flux_cpu(U, dT, Tsemi, dQ, h, in_z0_m, 
-    zeta, Rib, Re, B, z0_m, z0_t, Rib_conv_lim, Cm, Ct, Km, Pr_t_inv,
+    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_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, 
@@ -188,8 +192,8 @@ void Flux<T, RunMem, memIn>::compute_flux(T *zeta_, T *Rib_, T *Re_, T *B_, T *z
     grid_size);
 
 #ifdef INCLUDE_CUDA
-    else compute_flux_gpu(U, dT, Tsemi, dQ, h, in_z0_m, 
-    zeta, Rib, Re, B, z0_m, z0_t, Rib_conv_lim, Cm, Ct, Km, Pr_t_inv,
+    else compute_flux_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, 
@@ -230,17 +234,4 @@ template class Flux<double, MemType::CPU, MemType::CPU>;
     template class Flux<double, MemType::GPU, MemType::GPU>;
     template class Flux<double, MemType::GPU, MemType::CPU>;
     template class Flux<double, MemType::CPU, MemType::GPU>;
-#endif 
-
-template class Compute_Flux<float, MemType::CPU, MemType::CPU>;
-template class Compute_Flux<double, MemType::CPU, MemType::CPU>;
-
-#ifdef INCLUDE_CUDA
-    template class Compute_Flux<float, MemType::GPU, MemType::GPU>;
-    template class Compute_Flux<float, MemType::GPU, MemType::CPU>;
-    template class Compute_Flux<float, MemType::CPU, MemType::GPU>;
-
-    template class Compute_Flux<double, MemType::GPU, MemType::GPU>;
-    template class Compute_Flux<double, MemType::GPU, MemType::CPU>;
-    template class Compute_Flux<double, MemType::CPU, MemType::GPU>;
 #endif 
\ No newline at end of file
diff --git a/srcCXX/FluxComputeFunc.cpp b/srcCXX/sfx_flux_compute_func.cpp
similarity index 69%
rename from srcCXX/FluxComputeFunc.cpp
rename to srcCXX/sfx_flux_compute_func.cpp
index 15fbb6c28cffcbd3738db9b2e0f2ad0e0351fee0..1e6bebb16a0add5ff9621a5e52e25fb5b5a1d4c2 100644
--- a/srcCXX/FluxComputeFunc.cpp
+++ b/srcCXX/sfx_flux_compute_func.cpp
@@ -1,12 +1,12 @@
 #include <cmath>
 #include <iostream>
-#include "../includeCXX/FluxComputeFunc.h"
+#include "../includeCXX/sfx_flux_compute_func.h"
 
 template<typename T>
-void get_charnock_roughness(const T h, const T U,
+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,
-    T &z0_m, T &u_dyn0, 
+    const T h_charnock, const T c1_charnock, const T c2_charnock, 
     const int maxiters)
 {
     T Uc, a, b, c, c_min, f;
@@ -35,22 +35,22 @@ void get_charnock_roughness(const T h, const T U,
     u_dyn0 = Uc / c;
 }
 
-template void get_charnock_roughness(const float h, const float U,
+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,
-    float &z0_m, float &u_dyn0, 
     const int maxiters);
-template void get_charnock_roughness(const double h, const double U,
+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,
-    double &z0_m, double &u_dyn0, 
+    const double h_charnock, const double c1_charnock, const double c2_charnock, 
     const int maxiters);
 
 template<typename T>
-void get_convection_lim(const T h0_m, const T h0_t, const T B,
+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 &zeta_lim, T &Rib_lim, T &f_m_lim, T &f_h_lim)
+    const T alpha_h, const T alpha_m, const T alpha_h_fix)
 {
     T psi_m, psi_h, f_m, f_h, c;
 
@@ -73,19 +73,19 @@ void get_convection_lim(const T h0_m, const T h0_t, const T B,
     Rib_lim = zeta_lim * psi_h / (psi_m * psi_m);
 }
 
-template void get_convection_lim(const float h0_m, const float h0_t, const float B,
+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,
-    float &zeta_lim, float &Rib_lim, float &f_m_lim, float &f_h_lim);
-template void get_convection_lim(const double h0_m, const double h0_t, const double B,
+    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,
-    double &zeta_lim, double &Rib_lim, double &f_m_lim, double &f_h_lim);
+    const double alpha_h, const double alpha_m, const double alpha_h_fix);
 
 template<typename T>
-void get_psi_stable(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 &psi_m, T &psi_h, T &zeta)
+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;
     
@@ -102,19 +102,19 @@ void get_psi_stable(const T Rib, const T h0_m, const T h0_t, const T B,
     psi_h = (psi0_m + B) / Pr_t_0_inv + phi;
 }
 
-template void get_psi_stable(const float Rib, const float h0_m, const float h0_t, const float B,
-    const float Pr_t_0_inv, const float beta_m,
-    float &psi_m, float &psi_h, float &zeta);
-template void get_psi_stable(const double Rib, const double h0_m, const double h0_t, const double B,
-    const double Pr_t_0_inv, const double beta_m,
-    double &psi_m, double &psi_h, double &zeta);
+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(const T Rib, const T h0_m, const T h0_t, const T B, 
+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,
-    T &psi_m, T &psi_h, T &zeta, 
     const int maxiters)
 {
     T zeta0_m, zeta0_h, f0_m, f0_h, p_m, p_h, a_m, a_h, c_lim, f;
@@ -150,23 +150,23 @@ void get_psi_convection(const T Rib, const T h0_m, const T h0_t, const T B,
     }
 }
 
-template void get_psi_convection(const float Rib, const float h0_m, const float h0_t, const float B, 
+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,
-    float &psi_m, float &psi_h, float &zeta, 
     const int maxiters);
-template void get_psi_convection(const double Rib, const double h0_m, const double h0_t, const double B, 
+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,
-    double &psi_m, double &psi_h, double &zeta, 
+    const double alpha_h, const double alpha_m, const double alpha_h_fix, 
     const int maxiters);
 
 template<typename T>
-void get_psi_neutral(const T h0_m, const T h0_t, const T B,
-    const T Pr_t_0_inv,
-    T &psi_m, T &psi_h, T &zeta)
+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);
@@ -175,18 +175,18 @@ void get_psi_neutral(const T h0_m, const T h0_t, const T B,
         psi_h = psi_m / Pr_t_0_inv;
 }
 
-template void get_psi_neutral(const float h0_m, const float h0_t, const float B,
-    const float Pr_t_0_inv,
-    float &psi_m, float &psi_h, float &zeta);
-template void get_psi_neutral(const double h0_m, const double h0_t, const double B,
-    const double Pr_t_0_inv,
-    double &psi_m, double &psi_h, double &zeta);
+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(const T Rib, const T h0_m, const T h0_t, const T B, 
+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,
-    T &psi_m, T &psi_h, T &zeta,
     const int maxiters)
 {
     T zeta0_m, zeta0_h, f0_m, f0_h, f_m, f_h;
@@ -225,20 +225,20 @@ void get_psi_semi_convection(const T Rib, const T h0_m, const T h0_t, const T B,
     }
 }
 
-template void get_psi_semi_convection(const float Rib, const float h0_m, const float h0_t, const float B, 
+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,
-    float &psi_m, float &psi_h, float &zeta,
     const int maxiters);
-template void get_psi_semi_convection(const double Rib, const double h0_m, const double h0_t, const double B, 
+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,
-    double &psi_m, double &psi_h, double &zeta,
     const int maxiters);
 
 template<typename T>
-void compute_flux_cpu(const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_, const T *h_, const T *in_z0_m_,
-    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_,
+void compute_flux_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, 
@@ -271,7 +271,7 @@ void compute_flux_cpu(const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_,
 
         if (surface_type == 0)
         {
-            get_charnock_roughness(h, U, kappa, h_charnock, c1_charnock, c2_charnock, z0_m, u_dyn0, maxiters_charnock);
+            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) 
@@ -293,13 +293,13 @@ void compute_flux_cpu(const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_,
         h0_t = h / z0_t;
         Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / (U*U);
 
-        get_convection_lim(h0_m, h0_t, B, Pr_t_inf_inv, Pr_t_0_inv, alpha_h, alpha_m, alpha_h_fix, zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim);
+        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(Rib, h0_m, h0_t, B, Pr_t_0_inv, beta_m, psi_m, psi_h, zeta);
+            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;
@@ -308,7 +308,7 @@ void compute_flux_cpu(const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_,
 
         else if (Rib < Rib_conv_lim) 
         {
-            get_psi_convection(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, psi_m, psi_h, zeta, 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, 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;
@@ -316,14 +316,14 @@ void compute_flux_cpu(const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_,
         }
         else if (Rib > -0.001) 
         {
-            get_psi_neutral(h0_m, h0_t, B, Pr_t_0_inv, psi_m, psi_h, zeta);
+            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(Rib, h0_m, h0_t, B, Pr_t_0_inv, alpha_m, alpha_h_fix, psi_m, psi_h, zeta, maxiters_convection);
+            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));
@@ -349,8 +349,8 @@ void compute_flux_cpu(const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_,
     }
 }
 
-template void compute_flux_cpu(const float *U, const float *dt, const float *T_semi, const float *dq, const float *H, const float *in_z0_m,
-    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_,
+template void compute_flux_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, 
@@ -360,8 +360,8 @@ template void compute_flux_cpu(const float *U, const float *dt, const float *T_s
     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_cpu(const double *U, const double *dt, const double *T_semi, const double *dq, const double *H, const double *in_z0_m, 
-    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_,
+template void compute_flux_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, 
diff --git a/srcCXX/MemoryProcessing.cpp b/srcCXX/sfx_memory_processing.cpp
similarity index 96%
rename from srcCXX/MemoryProcessing.cpp
rename to srcCXX/sfx_memory_processing.cpp
index d3876059de4b7e03dcafef383d200f89c71a8dc8..1d5805a4716eb7c1068a5665c432a8caa95441ba 100644
--- a/srcCXX/MemoryProcessing.cpp
+++ b/srcCXX/sfx_memory_processing.cpp
@@ -1,4 +1,4 @@
-#include "../includeCXX/MemoryProcessing.h"
+#include "../includeCXX/sfx_memory_processing.h"
 #include <cstdlib>
 #include <cstring>
 
diff --git a/srcF/sfx_esm.f90 b/srcF/sfx_esm.f90
index 3484799b7e0a3d10155500c1b8e8bb114788ae97..0eb077c0239b4dc68037a4c7d09f70750f4839e3 100644
--- a/srcF/sfx_esm.f90
+++ b/srcF/sfx_esm.f90
@@ -55,9 +55,9 @@ contains
         integer i
         ! ----------------------------------------------------------------------------
 #if defined(INCLUDE_CUDA) || defined(INCLUDE_CXX)
-        call surf_flux(meteo%U, meteo%dT, meteo%Tsemi, meteo%dQ, meteo%h, meteo%z0_m, &
-        sfx%zeta, sfx%Rib, sfx%Re, sfx%B, sfx%z0_m, sfx%z0_t,  &
+        call surf_flux(sfx%zeta, sfx%Rib, sfx%Re, sfx%B, sfx%z0_m, sfx%z0_t, &
         sfx%Rib_conv_lim, sfx%Cm, sfx%Ct, sfx%Km, sfx%Pr_t_inv, &
+        meteo%U, meteo%dT, meteo%Tsemi, meteo%dQ, meteo%h, meteo%z0_m, &
         kappa, Pr_t_0_inv, Pr_t_inf_inv, & 
         alpha_m, alpha_h, alpha_h_fix, & 
         beta_m, beta_h, Rib_max, Re_rough_min, & 
diff --git a/srcF/FCWrapper.F90 b/srcF/sfx_fc_wrapper.F90
similarity index 91%
rename from srcF/FCWrapper.F90
rename to srcF/sfx_fc_wrapper.F90
index b946e9756a99705fe7fb1d14c5d7ce4b238e83ee..6cd9ea26077d6194f7d3bafb1bc6f285cc9f935a 100644
--- a/srcF/FCWrapper.F90
+++ b/srcF/sfx_fc_wrapper.F90
@@ -1,9 +1,9 @@
 module C_FUNC
     INTERFACE
 #if defined(INCLUDE_CUDA) || defined(INCLUDE_CXX)
-        SUBROUTINE surf_flux(U, dT, Tsemi, dQ, h, in_z0_m, & 
-            zeta, Rib, Re, B, z0_m, z0_t, Rib_conv_lim, &
+        SUBROUTINE surf_flux(zeta, Rib, Re, B, z0_m, z0_t, Rib_conv_lim, &
             Cm, Ct, Km, Prt_inv, & 
+            U, dT, Tsemi, dQ, h, in_z0_m, & 
             kappa, Pr_t_0_inv, Pr_t_inf_inv, & 
             alpha_m, alpha_h, alpha_h_fix, & 
             beta_m, beta_h, Rib_max, Re_rough_min, &