Merge branch 'Merge2Main' into 'main'

fixes

See merge request inmcm60_pbl/inmcm60_sfx!15

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)

includeCU/FluxComputeFunc.cuh → includeCU/sfx_flux_compute_func.cuh

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

includeCU/MemoryProcessing.cuh → includeCU/sfx_memory_processing.cuh

 #pragma once
-#include "../includeCXX/TemplateParameters.h"
+#include "../includeCXX/sfx_template_parameters.h"
 #include <cstddef>
 
 namespace memproc

includeCXX/SubFunctions.h → 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, 

includeCXX/Flux.h → includeCXX/sfx_flux.h

 #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

includeCXX/FluxComputeFunc.h → includeCXX/sfx_flux_compute_func.h

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

includeCXX/MemoryProcessing.h → includeCXX/sfx_memory_processing.h

 #pragma once
-#include "TemplateParameters.h"
+#include "sfx_template_parameters.h"
 #include <cstddef>
 
 namespace memproc

srcC/SubFunctionsWrapper.c → srcC/sfx_call_cxx.c

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

srcCU/FluxComputeFunc.cu → srcCU/sfx_flux_compute_func.cu

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

srcCU/MemoryProcessing.cu → srcCU/sfx_memory_processing.cu

-#include "../includeCU/MemoryProcessing.cuh"
+#include "../includeCU/sfx_memory_processing.cuh"
 #include <cuda.h>
 #include <cuda_runtime_api.h>
 

srcCXX/SubFunctions.cpp → srcCXX/sfx_call_class_func.cpp

 #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

srcCXX/Flux.cpp → srcCXX/sfx_flux.cpp

 #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

srcCXX/FluxComputeFunc.cpp → srcCXX/sfx_flux_compute_func.cpp

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

srcCXX/MemoryProcessing.cpp → srcCXX/sfx_memory_processing.cpp

-#include "../includeCXX/MemoryProcessing.h"
+#include "../includeCXX/sfx_memory_processing.h"
 #include <cstdlib>
 #include <cstring>
 

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, & 

srcF/FCWrapper.F90 → srcF/sfx_fc_wrapper.F90

 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, &