diff --git a/includeCU/sfx-model-compute-subfunc.cuh b/includeCU/sfx-model-compute-subfunc.cuh
index 7acf0340125de4f13acddaf9c105f4e255e83a20..0b45bf965920a4eff027e08eaf158f4ec5c840c1 100644
--- a/includeCU/sfx-model-compute-subfunc.cuh
+++ b/includeCU/sfx-model-compute-subfunc.cuh
@@ -9,22 +9,22 @@ FUCNTION_DECLARATION_SPECIFIER void get_convection_lim(T &zeta_lim, T &Rib_lim,
 {
     T psi_m, psi_h, f_m, f_h, c;
 
-    c = pow(param.Pr_t_inf_inv / param.Pr_t_0_inv, 4);
+    c = powf(param.Pr_t_inf_inv / param.Pr_t_0_inv, 4);
     zeta_lim = (2.0 * param.alpha_h - c * param.alpha_m - 
-                sqrt( (c * param.alpha_m)*(c * param.alpha_m) + 4.0 * c * param.alpha_h * (param.alpha_h - param.alpha_m))) / (2.0 * param.alpha_h*param.alpha_h);
+                sqrtf( (c * param.alpha_m)*(c * param.alpha_m) + 4.0 * c * param.alpha_h * (param.alpha_h - param.alpha_m))) / (2.0 * param.alpha_h*param.alpha_h);
 
-    f_m_lim = pow(1.0 - param.alpha_m * zeta_lim, 0.25);
-    f_h_lim = sqrt(1.0 - param.alpha_h * zeta_lim);
+    f_m_lim = powf(1.0 - param.alpha_m * zeta_lim, 0.25);
+    f_h_lim = sqrtf(1.0 - param.alpha_h * zeta_lim);
 
     f_m = zeta_lim / h0_m;
     f_h = zeta_lim / h0_t;
     if (fabs(B) < 1.0e-10) f_h = f_m;
 
-    f_m = pow(1.0 - param.alpha_m * f_m, 0.25);
-    f_h = sqrt(1.0 - param.alpha_h_fix * f_h);
+    f_m = powf(1.0 - param.alpha_m * f_m, 0.25);
+    f_h = sqrtf(1.0 - param.alpha_h_fix * f_h);
 
-    psi_m = 2.0 * (atan(f_m_lim) - atan(f_m)) + log((f_m_lim - 1.0) * (f_m + 1.0)/((f_m_lim + 1.0) * (f_m - 1.0)));
-    psi_h = log((f_h_lim - 1.0) * (f_h + 1.0)/((f_h_lim + 1.0) * (f_h - 1.0))) / param.Pr_t_0_inv;
+    psi_m = 2.0 * (atanf(f_m_lim) - atanf(f_m)) + logf((f_m_lim - 1.0) * (f_m + 1.0)/((f_m_lim + 1.0) * (f_m - 1.0)));
+    psi_h = logf((f_h_lim - 1.0) * (f_h + 1.0)/((f_h_lim + 1.0) * (f_h - 1.0))) / param.Pr_t_0_inv;
 
     Rib_lim = zeta_lim * psi_h / (psi_m * psi_m);
 }
@@ -36,12 +36,12 @@ FUCNTION_DECLARATION_SPECIFIER void get_psi_stable(T &psi_m, T &psi_h, T &zeta,
 {
     T Rib_coeff, psi0_m, psi0_h, phi, c;
     
-    psi0_m = log(h0_m);
+    psi0_m = logf(h0_m);
     psi0_h = B / psi0_m;
 
     Rib_coeff = param.beta_m * Rib;
     c = (psi0_h + 1.0) / param.Pr_t_0_inv - 2.0 * Rib_coeff;
-    zeta = psi0_m * (sqrt(c*c + 4.0 * Rib_coeff * (1.0 - Rib_coeff)) - c) / (2.0 * param.beta_m * (1.0 - Rib_coeff));
+    zeta = psi0_m * (sqrtf(c*c + 4.0 * Rib_coeff * (1.0 - Rib_coeff)) - c) / (2.0 * param.beta_m * (1.0 - Rib_coeff));
 
     phi = param.beta_m * zeta;
 
@@ -58,8 +58,8 @@ FUCNTION_DECLARATION_SPECIFIER void get_psi_convection(T &psi_m, T &psi_h, T &ze
 {
     T zeta0_m, zeta0_h, f0_m, f0_h, p_m, p_h, a_m, a_h, c_lim, f;
 
-    p_m = 2.0 * atan(f_m_conv_lim) + log((f_m_conv_lim - 1.0) / (f_m_conv_lim + 1.0));
-    p_h = log((f_h_conv_lim - 1.0) / (f_h_conv_lim + 1.0));
+    p_m = 2.0 * atanf(f_m_conv_lim) + logf((f_m_conv_lim - 1.0) / (f_m_conv_lim + 1.0));
+    p_h = logf((f_h_conv_lim - 1.0) / (f_h_conv_lim + 1.0));
 
     zeta = zeta_conv_lim;
 
@@ -70,13 +70,13 @@ FUCNTION_DECLARATION_SPECIFIER void get_psi_convection(T &psi_m, T &psi_h, T &ze
         if (fabs(B) < 1.0e-10) 
             zeta0_h = zeta0_m;
 
-        f0_m = pow(1.0 - param.alpha_m * zeta0_m, 0.25);
-        f0_h = sqrt(1.0 - param.alpha_h_fix * zeta0_h);
+        f0_m = powf(1.0 - param.alpha_m * zeta0_m, 0.25);
+        f0_h = sqrtf(1.0 - param.alpha_h_fix * zeta0_h);
 
-        a_m = -2.0*atan(f0_m) + log((f0_m + 1.0)/(f0_m - 1.0));
-        a_h = log((f0_h + 1.0)/(f0_h - 1.0));
+        a_m = -2.0*atanf(f0_m) + logf((f0_m + 1.0)/(f0_m - 1.0));
+        a_h = logf((f0_h + 1.0)/(f0_h - 1.0));
 
-        c_lim = pow(zeta_conv_lim / zeta, 1.0 / 3.0);
+        c_lim = powf(zeta_conv_lim / zeta, 1.0 / 3.0);
         f = 3.0 * (1.0 - c_lim);
 
         psi_m = f / f_m_conv_lim + p_m + a_m;
@@ -95,8 +95,8 @@ FUCNTION_DECLARATION_SPECIFIER void get_psi_neutral(T &psi_m, T &psi_h, T &zeta,
     const sfx_param& param)
 {
     zeta = 0.0;
-    psi_m = log(h0_m);
-    psi_h = log(h0_t) / param.Pr_t_0_inv;
+    psi_m = logf(h0_m);
+    psi_h = logf(h0_t) / param.Pr_t_0_inv;
     if (fabs(B) < 1.0e-10) 
         psi_h = psi_m / param.Pr_t_0_inv;
 }
@@ -109,8 +109,8 @@ FUCNTION_DECLARATION_SPECIFIER void get_psi_semi_convection(T &psi_m, T &psi_h,
 {
     T zeta0_m, zeta0_h, f0_m, f0_h, f_m, f_h;
 
-    psi_m = log(h0_m);
-    psi_h = log(h0_t);
+    psi_m = logf(h0_m);
+    psi_h = logf(h0_t);
 
     if (fabs(B) < 1.0e-10) 
         psi_h = psi_m;
@@ -124,17 +124,17 @@ FUCNTION_DECLARATION_SPECIFIER void get_psi_semi_convection(T &psi_m, T &psi_h,
         if (fabs(B) < 1.0e-10) 
             zeta0_h = zeta0_m;
 
-        f_m = pow(1.0 - param.alpha_m * zeta, 0.25e0);
-        f_h = sqrt(1.0 - param.alpha_h_fix * zeta);
+        f_m = powf(1.0 - param.alpha_m * zeta, 0.25e0);
+        f_h = sqrtf(1.0 - param.alpha_h_fix * zeta);
 
-        f0_m = pow(1.0 - param.alpha_m * zeta0_m, 0.25e0);
-        f0_h = sqrt(1.0 - param.alpha_h_fix * zeta0_h);
+        f0_m = powf(1.0 - param.alpha_m * zeta0_m, 0.25e0);
+        f0_h = sqrtf(1.0 - param.alpha_h_fix * zeta0_h);
 
         f0_m = sfx_math::max(f0_m, T(1.000001e0));
         f0_h = sfx_math::max(f0_h, T(1.000001e0));
 
-        psi_m = log((f_m - 1.0e0)*(f0_m + 1.0e0)/((f_m + 1.0e0)*(f0_m - 1.0e0))) + 2.0e0*(atan(f_m) - atan(f0_m));
-        psi_h = log((f_h - 1.0e0)*(f0_h + 1.0e0)/((f_h + 1.0e0)*(f0_h - 1.0e0))) / param.Pr_t_0_inv;
+        psi_m = logf((f_m - 1.0e0)*(f0_m + 1.0e0)/((f_m + 1.0e0)*(f0_m - 1.0e0))) + 2.0e0*(atanf(f_m) - atanf(f0_m));
+        psi_h = logf((f_h - 1.0e0)*(f0_h + 1.0e0)/((f_h + 1.0e0)*(f0_h - 1.0e0))) / param.Pr_t_0_inv;
 
         if (i == maxiters_convection + 1) 
             break;
@@ -153,23 +153,23 @@ FUCNTION_DECLARATION_SPECIFIER void get_psi(T &psi_m, T &psi_h,
 
     if (zeta >= 0.0) 
     {
-        q_m = pow((1.0 - param.b_m) / param.b_m, 1.0 / 3.0);
-        q_h = sqrt(param.c_h * param.c_h - 4.0);
+        q_m = powf((1.0 - param.b_m) / param.b_m, 1.0 / 3.0);
+        q_h = sqrtf(param.c_h * param.c_h - 4.0);
 
-        x_m = pow(1.0 + zeta, 1.0 / 3.0);
+        x_m = powf(1.0 + zeta, 1.0 / 3.0);
         x_h = zeta;
 
-        psi_m = -3.0 * (param.a_m / param.b_m) * (x_m - 1.0) + 0.5 * (param.a_m / param.b_m) * q_m * (2.0 * log((x_m + q_m) / (1.0 + q_m)) - log((x_m * x_m - x_m * q_m + q_m * q_m) / (1.0 - q_m + q_m * q_m)) + 2.0 * sqrt(3.0) * (atan((2.0 * x_m - q_m) / (sqrt(3.0) * q_m)) - atan((2.0 - q_m) / (sqrt(3.0) * q_m))));
+        psi_m = -3.0 * (param.a_m / param.b_m) * (x_m - 1.0) + 0.5 * (param.a_m / param.b_m) * q_m * (2.0 * logf((x_m + q_m) / (1.0 + q_m)) - logf((x_m * x_m - x_m * q_m + q_m * q_m) / (1.0 - q_m + q_m * q_m)) + 2.0 * sqrtf(3.0) * (atanf((2.0 * x_m - q_m) / (sqrtf(3.0) * q_m)) - atanf((2.0 - q_m) / (sqrtf(3.0) * q_m))));
 
-        psi_h = -0.5 * param.b_h * log(1.0 + param.c_h * x_h + x_h * x_h) + ((-param.a_h / q_h) + ((param.b_h * param.c_h) / (2.0 * q_h))) * (log((2.0 * x_h + param.c_h - q_h) / (2.0 * x_h + param.c_h + q_h)) - log((param.c_h - q_h) / (param.c_h + q_h)));
+        psi_h = -0.5 * param.b_h * logf(1.0 + param.c_h * x_h + x_h * x_h) + ((-param.a_h / q_h) + ((param.b_h * param.c_h) / (2.0 * q_h))) * (logf((2.0 * x_h + param.c_h - q_h) / (2.0 * x_h + param.c_h + q_h)) - logf((param.c_h - q_h) / (param.c_h + q_h)));
     }
     else
     {
-        x_m = pow(1.0 - param.alpha_m * zeta, 0.25);
-        x_h = pow(1.0 - param.alpha_h * zeta, 0.25); 
+        x_m = powf(1.0 - param.alpha_m * zeta, 0.25);
+        x_h = powf(1.0 - param.alpha_h * zeta, 0.25); 
 
-        psi_m = (4.0 * atan(1.0) / 2.0) + 2.0 * log(0.5 * (1.0 + x_m)) + log(0.5 * (1.0 + x_m * x_m)) - 2.0 * atan(x_m);
-        psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h));
+        psi_m = (4.0 * atanf(1.0) / 2.0) + 2.0 * logf(0.5 * (1.0 + x_m)) + logf(0.5 * (1.0 + x_m * x_m)) - 2.0 * atanf(x_m);
+        psi_h = 2.0 * logf(0.5 * (1.0 + x_h * x_h));
     }
 }
 
@@ -183,28 +183,28 @@ FUCNTION_DECLARATION_SPECIFIER void get_psi_mh(T &psi_m, T &psi_h,
 
     if (zeta_m >= 0.0) 
     {
-        q_m = pow((1.0 - param.b_m) / param.b_m, 1.0 / 3.0);
-        x_m = pow(1.0 + zeta_m, 1.0 / 3.0);
+        q_m = powf((1.0 - param.b_m) / param.b_m, 1.0 / 3.0);
+        x_m = powf(1.0 + zeta_m, 1.0 / 3.0);
 
-        psi_m = -3.0 * (param.a_m / param.b_m) * (x_m - 1.0) + 0.5 * (param.a_m / param.b_m) * q_m * (2.0 * log((x_m + q_m) / (1.0 + q_m)) - log((x_m * x_m - x_m * q_m + q_m * q_m) / (1.0 - q_m + q_m * q_m)) + 2.0 * sqrt(3.0) * (atan((2.0 * x_m - q_m) / (sqrt(3.0) * q_m)) - atan((2.0 - q_m) / (sqrt(3.0) * q_m))));
+        psi_m = -3.0 * (param.a_m / param.b_m) * (x_m - 1.0) + 0.5 * (param.a_m / param.b_m) * q_m * (2.0 * logf((x_m + q_m) / (1.0 + q_m)) - logf((x_m * x_m - x_m * q_m + q_m * q_m) / (1.0 - q_m + q_m * q_m)) + 2.0 * sqrtf(3.0) * (atanf((2.0 * x_m - q_m) / (sqrtf(3.0) * q_m)) - atanf((2.0 - q_m) / (sqrtf(3.0) * q_m))));
     }                                
     else
     {    
-        x_m = pow(1.0 - param.alpha_m * zeta_m, 0.25);
-        psi_m = (4.0 * atan(1.0) / 2.0) + 2.0 * log(0.5 * (1.0 + x_m)) + log(0.5 * (1.0 + x_m * x_m)) - 2.0 * atan(x_m);
+        x_m = powf(1.0 - param.alpha_m * zeta_m, 0.25);
+        psi_m = (4.0 * atanf(1.0) / 2.0) + 2.0 * logf(0.5 * (1.0 + x_m)) + logf(0.5 * (1.0 + x_m * x_m)) - 2.0 * atanf(x_m);
     }
 
     if (zeta_h >= 0.0)
     {    
-        q_h = sqrt(param.c_h * param.c_h - 4.0);
+        q_h = sqrtf(param.c_h * param.c_h - 4.0);
         x_h = zeta_h;
 
-        psi_h = -0.5 * param.b_h * log(1.0 + param.c_h * x_h + x_h * x_h) + ((-param.a_h / q_h) + ((param.b_h * param.c_h) / (2.0 * q_h))) * (log((2.0 * x_h + param.c_h - q_h) / (2.0 * x_h + param.c_h + q_h)) - log((param.c_h - q_h) / (param.c_h + q_h)));
+        psi_h = -0.5 * param.b_h * logf(1.0 + param.c_h * x_h + x_h * x_h) + ((-param.a_h / q_h) + ((param.b_h * param.c_h) / (2.0 * q_h))) * (logf((2.0 * x_h + param.c_h - q_h) / (2.0 * x_h + param.c_h + q_h)) - logf((param.c_h - q_h) / (param.c_h + q_h)));
     }
     else
     {
-        x_h = pow(1.0 - param.alpha_h * zeta_h, 0.25);
-        psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h));
+        x_h = powf(1.0 - param.alpha_h * zeta_h, 0.25);
+        psi_h = 2.0 * logf(0.5 * (1.0 + x_h * x_h));
     }
 }
 
@@ -220,9 +220,9 @@ FUCNTION_DECLARATION_SPECIFIER void get_dynamic_scales(T &Udyn, T &Tdyn, T &Qdyn
     const T gamma = T(0.61);
     T psi_m, psi_h, psi0_m, psi0_h, Linv;
 
-    Udyn = param.kappa * U / log(z / z0_m);
-    Tdyn = param.kappa * dT * param.Pr_t_0_inv / log(z / z0_t);
-    Qdyn = param.kappa * dQ * param.Pr_t_0_inv / log(z / z0_t);
+    Udyn = param.kappa * U / logf(z / z0_m);
+    Tdyn = param.kappa * dT * param.Pr_t_0_inv / logf(z / z0_t);
+    Qdyn = param.kappa * dQ * param.Pr_t_0_inv / logf(z / z0_t);
     zeta = 0.0;
 
     // --- no wind
@@ -242,9 +242,9 @@ FUCNTION_DECLARATION_SPECIFIER void get_dynamic_scales(T &Udyn, T &Tdyn, T &Qdyn
         get_psi_mh(psi0_m, psi0_h, z0_m * Linv, z0_t * Linv, 
             param);
 
-        Udyn = param.kappa * U / (log(z / z0_m) - (psi_m - psi0_m));
-        Tdyn = param.kappa * dT * param.Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h));
-        Qdyn = param.kappa * dQ * param.Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h));
+        Udyn = param.kappa * U / (logf(z / z0_m) - (psi_m - psi0_m));
+        Tdyn = param.kappa * dT * param.Pr_t_0_inv / (logf(z / z0_t) - (psi_h - psi0_h));
+        Qdyn = param.kappa * dQ * param.Pr_t_0_inv / (logf(z / z0_t) - (psi_h - psi0_h));
 
         if (Udyn < 1e-5) 
             break;
@@ -261,12 +261,12 @@ FUCNTION_DECLARATION_SPECIFIER void get_phi(T &phi_m, T &phi_h,
 {
     if (zeta >= 0.0) 
     {
-        phi_m = 1.0 + (param.a_m * zeta * pow(1.0 + zeta, 1.0 / 3.0) ) / (1.0 + param.b_m * zeta);
+        phi_m = 1.0 + (param.a_m * zeta * powf(1.0 + zeta, 1.0 / 3.0) ) / (1.0 + param.b_m * zeta);
         phi_h = 1.0 + (param.a_h * zeta + param.b_h * zeta * zeta) / (1.0 + param.c_h * zeta + zeta * zeta);
     }
     else
     {
-        phi_m = pow(1.0 - param.alpha_m * zeta, -0.25);
-        phi_h = pow(1.0 - param.alpha_h * zeta, -0.5);
+        phi_m = powf(1.0 - param.alpha_m * zeta, -0.25);
+        phi_h = powf(1.0 - param.alpha_h * zeta, -0.5);
     }
 }
diff --git a/includeCU/sfx-surface.cuh b/includeCU/sfx-surface.cuh
index 85908c1abc0a3e9f161ffa1fdbe7ab7dc07a5d29..c4a14c13195558c64cf4a695da731c8f8b0a4c56 100644
--- a/includeCU/sfx-surface.cuh
+++ b/includeCU/sfx-surface.cuh
@@ -8,9 +8,9 @@ FUCNTION_DECLARATION_SPECIFIER void get_thermal_roughness(T &z0_t, T &B,
     const struct sfx_surface_param& param,
     const int surface_type)
 {
-    // --- define B = log(z0_m / z0_t)
-    B = (Re <= param.Re_rough_min) ? param.B1_rough * log(param.B3_rough * Re) + param.B2_rough : 
-                                                 param.B4_rough * (pow(Re, param.B2_rough));
+    // --- define B = logf(z0_m / z0_t)
+    B = (Re <= param.Re_rough_min) ? param.B1_rough * logf(param.B3_rough * Re) + param.B2_rough : 
+                                                 param.B4_rough * (powf(Re, param.B2_rough));
 
     //   --- apply max restriction based on surface type
     if (surface_type == param.surface_ocean)  B = sfx_math::min(B, param.B_max_ocean);
@@ -32,22 +32,22 @@ FUCNTION_DECLARATION_SPECIFIER void get_charnock_roughness(T &z0_m, T &u_dyn0,
     Uc = U;
     a = 0.0;
     b = 25.0;
-    c_min = log(surface_param.h_charnock) / surface_param.kappa;
+    c_min = logf(surface_param.h_charnock) / surface_param.kappa;
 
     for (int i = 0; i < maxiters; i++)
     {
-        f = surface_param.c1_charnock - 2.0 * log(Uc);
+        f = surface_param.c1_charnock - 2.0 * logf(Uc);
         for (int j = 0; j < maxiters; j++)
         {
-            c = (f + 2.0 * log(b)) / surface_param.kappa;
+            c = (f + 2.0 * logf(b)) / surface_param.kappa;
             if (U <= 8.0e0) 
-                a = log(1.0 + surface_param.c2_charnock * ( pow(b / Uc, 3) ) ) / surface_param.kappa;
+                a = logf(1.0 + surface_param.c2_charnock * ( powf(b / Uc, 3) ) ) / surface_param.kappa;
             c = sfx_math::max(c - a, c_min);
             b = c;
         }
         z0_m = surface_param.h_charnock * exp(-c * surface_param.kappa);
         z0_m = sfx_math::max(z0_m, T(0.000015e0));
-        Uc = U * log(surface_param.h_charnock / z0_m) / log(h / z0_m);
+        Uc = U * logf(surface_param.h_charnock / z0_m) / logf(h / z0_m);
     }
     
     u_dyn0 = Uc / c;
diff --git a/srcCU/sfx-esm.cu b/srcCU/sfx-esm.cu
index 993db40535fa7a6da44f3353e8afb921bcd32c4b..67444d72161a309265d4b70062d99f369e301ecd 100644
--- a/srcCU/sfx-esm.cu
+++ b/srcCU/sfx-esm.cu
@@ -52,7 +52,7 @@ __global__ void sfx_kernel::compute_flux(sfxDataVecTypeC sfx,
         if (surface_type == surface.surface_land) 
         {
             h0_m = h / z0_m;
-            u_dyn0 = U * model.kappa / log(h0_m);
+            u_dyn0 = U * model.kappa / logf(h0_m);
         }
 
         Re = u_dyn0 * z0_m / phys.nu_air;
@@ -79,7 +79,7 @@ __global__ void sfx_kernel::compute_flux(sfxDataVecTypeC sfx,
         {
             get_psi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, model, numerics.maxiters_convection);
 
-            fval = pow(zeta_conv_lim / zeta, 1.0/3.0);
+            fval = powf(zeta_conv_lim / zeta, 1.0/3.0);
             phi_m = fval / f_m_conv_lim;
             phi_h = fval / (model.Pr_t_0_inv * f_h_conv_lim);
         }
@@ -94,8 +94,8 @@ __global__ void sfx_kernel::compute_flux(sfxDataVecTypeC sfx,
         {
             get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model, numerics.maxiters_convection);
             
-            phi_m = pow(1.0 - model.alpha_m * zeta, -0.25);
-            phi_h = 1.0 / (model.Pr_t_0_inv * sqrt(1.0 - model.alpha_h_fix * zeta));
+            phi_m = powf(1.0 - model.alpha_m * zeta, -0.25);
+            phi_h = 1.0 / (model.Pr_t_0_inv * sqrtf(1.0 - model.alpha_h_fix * zeta));
         }
 
         Cm = model.kappa / psi_m;
diff --git a/srcCU/sfx-sheba.cu b/srcCU/sfx-sheba.cu
index b4dd14b141801ced4584ab7c503fd3c83baa2e52..cb00722344507909c1699f2e7338186a42d1ede8 100644
--- a/srcCU/sfx-sheba.cu
+++ b/srcCU/sfx-sheba.cu
@@ -53,7 +53,7 @@ __global__ void sfx_kernel::compute_flux(sfxDataVecTypeC sfx,
         if (surface_type == surface.surface_land) 
         {
             h0_m = h / z0_m;
-            u_dyn0 = U * model.kappa / log(h0_m);
+            u_dyn0 = U * model.kappa / logf(h0_m);
         }
 
         Re = u_dyn0 * z0_m / phys.nu_air;
diff --git a/srcCXX/sfx-esm.cpp b/srcCXX/sfx-esm.cpp
index bbbcd2163c9080b18494a6d6cef182aadf3bb83c..5635c05a80a2e6f6bedcabe3114ed64637141904 100644
--- a/srcCXX/sfx-esm.cpp
+++ b/srcCXX/sfx-esm.cpp
@@ -55,7 +55,7 @@ void FluxEsm<T, memIn, memOut, MemType::CPU>::compute_flux()
         if (surface_type == surface.surface_land) 
         {
             h0_m = h / z0_m;
-            u_dyn0 = U * model.kappa / log(h0_m);
+            u_dyn0 = U * model.kappa / logf(h0_m);
         }
 
         Re = u_dyn0 * z0_m / phys.nu_air;
@@ -82,7 +82,7 @@ void FluxEsm<T, memIn, memOut, MemType::CPU>::compute_flux()
         {
             get_psi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, model, numerics.maxiters_convection);
 
-            fval = pow(zeta_conv_lim / zeta, 1.0/3.0);
+            fval = powf(zeta_conv_lim / zeta, 1.0/3.0);
             phi_m = fval / f_m_conv_lim;
             phi_h = fval / (model.Pr_t_0_inv * f_h_conv_lim);
         }
@@ -97,8 +97,8 @@ void FluxEsm<T, memIn, memOut, MemType::CPU>::compute_flux()
         {
             get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, model, numerics.maxiters_convection);
             
-            phi_m = pow(1.0 - model.alpha_m * zeta, -0.25);
-            phi_h = 1.0 / (model.Pr_t_0_inv * sqrt(1.0 - model.alpha_h_fix * zeta));
+            phi_m = powf(1.0 - model.alpha_m * zeta, -0.25);
+            phi_h = 1.0 / (model.Pr_t_0_inv * sqrtf(1.0 - model.alpha_h_fix * zeta));
         }
 
         Cm = model.kappa / psi_m;
diff --git a/srcCXX/sfx-sheba.cpp b/srcCXX/sfx-sheba.cpp
index 7d565650de0bfbba447007d3b14430ef7395dd0c..7901bfd9785dcc4242a4a175ee922bdd19697ea1 100644
--- a/srcCXX/sfx-sheba.cpp
+++ b/srcCXX/sfx-sheba.cpp
@@ -56,7 +56,7 @@ void FluxSheba<T, memIn, memOut, MemType::CPU>::compute_flux()
         if (surface_type == surface.surface_land) 
         {
             h0_m = h / z0_m;
-            u_dyn0 = U * model.kappa / log(h0_m);
+            u_dyn0 = U * model.kappa / logf(h0_m);
         }
 
         Re = u_dyn0 * z0_m / phys.nu_air;