From 9570b4ec27e945a92169782a91259cd900b690be Mon Sep 17 00:00:00 2001
From: Lizzzka007 <gashchuk2011@mail.ru>
Date: Wed, 20 Dec 2023 00:27:31 +0300
Subject: [PATCH] .
---
includeCU/sfx_compute_sheba.cuh | 17 +++
includeCU/sfx_surface.cuh | 9 ++
includeCXX/sfx_call_class_func.h | 17 ++-
includeCXX/sfx_compute_sheba.h | 11 +-
includeCXX/sfx_esm.h | 2 +-
includeCXX/sfx_sheba.h | 35 +++--
includeCXX/sfx_surface.h | 16 ++
srcC/sfx_call_cxx.c | 36 ++++-
srcCU/sfx_compute_esm.cu | 18 ++-
srcCU/sfx_surface.cu | 46 ++++++
srcCXX/sfx_call_class_func.cpp | 42 +++++-
srcCXX/sfx_compute_esm.cpp | 45 +-----
srcCXX/sfx_compute_sheba.cpp | 250 ++++++++++++++++++++++++++++---
srcCXX/sfx_esm.cpp | 2 +-
srcCXX/sfx_sheba.cpp | 81 ++++++----
srcCXX/sfx_surface.cpp | 87 +++++++++++
srcF/sfx_esm.f90 | 5 +-
srcF/sfx_fc_wrapper.F90 | 30 +++-
18 files changed, 614 insertions(+), 135 deletions(-)
create mode 100644 includeCU/sfx_compute_sheba.cuh
create mode 100644 includeCU/sfx_surface.cuh
create mode 100644 includeCXX/sfx_surface.h
create mode 100644 srcCU/sfx_surface.cu
create mode 100644 srcCXX/sfx_surface.cpp
diff --git a/includeCU/sfx_compute_sheba.cuh b/includeCU/sfx_compute_sheba.cuh
new file mode 100644
index 0000000..acd67d6
--- /dev/null
+++ b/includeCU/sfx_compute_sheba.cuh
@@ -0,0 +1,17 @@
+#pragma once
+
+template<typename T>
+void compute_flux_sheba_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 alpha_m, const T alpha_h,
+ const T a_m, const T a_h,
+ const T b_m, const T b_h,
+ const T c_h,
+ const T Re_rough_min,
+ const T B1_rough, const T B2_rough,
+ const T B_max_land, const T B_max_ocean, const T B_max_lake,
+ const T gamma_c, const T Re_visc_min,
+ const T Pr_m, const T nu_air, const T g,
+ const int maxiters_charnock,
+ const int grid_size);
\ No newline at end of file
diff --git a/includeCU/sfx_surface.cuh b/includeCU/sfx_surface.cuh
new file mode 100644
index 0000000..002336c
--- /dev/null
+++ b/includeCU/sfx_surface.cuh
@@ -0,0 +1,9 @@
+#pragma once
+
+// template<typename T>
+// __device__ void get_charnock_roughness(T &z0_m, T &u_dyn0,
+// const T h, const T U,
+// const T kappa,
+// const T h_charnock, const T c1_charnock, const T c2_charnock,
+// const int maxiters);
+
\ No newline at end of file
diff --git a/includeCXX/sfx_call_class_func.h b/includeCXX/sfx_call_class_func.h
index 08c9164..317f3d9 100644
--- a/includeCXX/sfx_call_class_func.h
+++ b/includeCXX/sfx_call_class_func.h
@@ -4,7 +4,7 @@
extern "C" {
#endif
- 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_,
+ void surf_flux_esm_CXX (float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
float *U_, float *dT_, float *Tsemi_, float *dQ_, float *h_, float *in_z0_m_,
const float kappa, const float Pr_t_0_inv, const float Pr_t_inf_inv,
const float alpha_m, const float alpha_h, const float alpha_h_fix,
@@ -16,6 +16,21 @@ extern "C" {
const int maxiters_charnock, const int maxiters_convection,
const int grid_size);
+ void surf_flux_sheba_CXX (float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
+ 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 alpha_m, const float alpha_h,
+ const float a_m, const float a_h,
+ const float b_m, const float b_h,
+ const float c_h,
+ const float Re_rough_min,
+ const float B1_rough, const float B2_rough,
+ const float B_max_land, const float B_max_ocean, const float B_max_lake,
+ const float gamma_c, const float Re_visc_min,
+ const float Pr_m, const float nu_air, const float g,
+ const int maxiters_charnock,
+ const int grid_size);
+
#ifdef __cplusplus
}
#endif
\ No newline at end of file
diff --git a/includeCXX/sfx_compute_sheba.h b/includeCXX/sfx_compute_sheba.h
index 92516a2..d13bf59 100644
--- a/includeCXX/sfx_compute_sheba.h
+++ b/includeCXX/sfx_compute_sheba.h
@@ -3,12 +3,15 @@
template<typename T>
void compute_flux_sheba_cpu(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_, T *Rib_conv_lim_, T *Cm_, T *Ct_, T *Km_, T *Pr_t_inv_,
const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_, const T *h_, const T *in_z0_m_,
- const T kappa, const T Pr_t_0_inv, const T Pr_t_inf_inv,
- const T alpha_m, const T alpha_h, const T alpha_h_fix,
- const T beta_m, const T beta_h, const T Rib_max, const T Re_rough_min,
+ const T kappa, const T Pr_t_0_inv,
+ const T alpha_m, const T alpha_h,
+ const T a_m, const T a_h,
+ const T b_m, const T b_h,
+ const T c_h,
+ const T Re_rough_min,
const T B1_rough, const T B2_rough,
const T B_max_land, const T B_max_ocean, const T B_max_lake,
const T gamma_c, const T Re_visc_min,
const T Pr_m, const T nu_air, const T g,
- const int maxiters_charnock, const int maxiters_convection,
+ const int maxiters_charnock,
const int grid_size);
\ No newline at end of file
diff --git a/includeCXX/sfx_esm.h b/includeCXX/sfx_esm.h
index 5f47957..1b7067c 100644
--- a/includeCXX/sfx_esm.h
+++ b/includeCXX/sfx_esm.h
@@ -33,7 +33,7 @@ public:
const T Pr_m, const T nu_air, const T g);
~FluxEsm();
- 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_esm(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_, T *Rib_conv_lim_, T *Cm_, T *Ct_, T *Km_, T *Pr_t_inv_,
T *U_, T *dT_, T *Tsemi_, T *dQ_, T *h_, T *in_z0_m_,
const int maxiters_charnock, const int maxiters_convection);
diff --git a/includeCXX/sfx_sheba.h b/includeCXX/sfx_sheba.h
index 3b2cfd3..ce911e0 100644
--- a/includeCXX/sfx_sheba.h
+++ b/includeCXX/sfx_sheba.h
@@ -9,33 +9,38 @@ private:
T *U, *dT, *Tsemi, *dQ, *h, *in_z0_m;
T *zeta, *Rib, *Re, *B, *z0_m, *z0_t, *Rib_conv_lim, *Cm, *Ct, *Km, *Pr_t_inv;
- // T kappa, Pr_t_0_inv, Pr_t_inf_inv,
- // alpha_m, alpha_h, alpha_h_fix,
- // beta_m, beta_h,
- // Rib_max, Re_rough_min,
- // B1_rough, B2_rough,
- // B_max_land, B_max_ocean, B_max_lake,
- // gamma_c,
- // Re_visc_min,
- // Pr_m, nu_air, g;
+ T kappa, Pr_t_0_inv, Pr_t_inf_inv,
+ alpha_m, alpha_h,
+ a_m, a_h,
+ b_m, b_h,
+ c_h,
+ Rib_max, Re_rough_min,
+ B1_rough, B2_rough,
+ B_max_land, B_max_ocean, B_max_lake,
+ gamma_c,
+ Re_visc_min,
+ Pr_m, nu_air, g;
int grid_size;
bool ifAllocated;
size_t allocated_size;
public:
FluxSheba();
- void set_params(const int grid_size, const T kappa, const T Pr_t_0_inv, const T Pr_t_inf_inv,
- const T alpha_m, const T alpha_h, const T alpha_h_fix,
- const T beta_m, const T beta_h, const T Rib_max, const T Re_rough_min,
+ void set_params(const int grid_size, const T kappa, const T Pr_t_0_inv,
+ const T alpha_m, const T alpha_h,
+ const float a_m, const float a_h,
+ const float b_m, const float b_h,
+ const float c_h,
+ const T Re_rough_min,
const T B1_rough, const T B2_rough,
const T B_max_land, const T B_max_ocean, const T B_max_lake,
const T gamma_c, const T Re_visc_min,
const T Pr_m, const T nu_air, const T g);
- ~Flux();
+ ~FluxSheba();
- 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_sheba(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);
+ const int maxiters_charnock);
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_,
diff --git a/includeCXX/sfx_surface.h b/includeCXX/sfx_surface.h
new file mode 100644
index 0000000..d82fec0
--- /dev/null
+++ b/includeCXX/sfx_surface.h
@@ -0,0 +1,16 @@
+#pragma once
+
+template<typename T>
+void get_charnock_roughness(T &z0_m, T &u_dyn0,
+ const T h, const T U,
+ const T kappa,
+ const T h_charnock, const T c1_charnock, const T c2_charnock,
+ const int maxiters);
+
+template<typename T>
+void get_thermal_roughness(T &z0_t, T &B,
+ const T z0_m, const T Re,
+ const T Re_rough_min,
+ const T B1_rough, const T B2_rough, const T B3_rough, const T B4_rough,
+ const T B_max_ocean, const T B_max_lake, const T B_max_land,
+ const int surface_type);
\ No newline at end of file
diff --git a/srcC/sfx_call_cxx.c b/srcC/sfx_call_cxx.c
index 007d00c..be43297 100644
--- a/srcC/sfx_call_cxx.c
+++ b/srcC/sfx_call_cxx.c
@@ -2,7 +2,8 @@
#include <stdio.h>
#include "../includeCXX/sfx_call_class_func.h"
// -------------------------------------------------------------------------- //
-void get_surface_fluxes (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 get_surface_fluxes_esm (float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
float *U_, float *dT_, float *Tsemi_, float *dQ_, float *h_, float *in_z0_m_,
const float *kappa, const float *Pr_t_0_inv, const float *Pr_t_inf_inv,
const float *alpha_m, const float *alpha_h, const float *alpha_h_fix,
@@ -14,7 +15,7 @@ void get_surface_fluxes (float *zeta_, float *Rib_, float *Re_, float *B_, float
const int *maxiters_charnock, const int *maxiters_convection,
const int *grid_size)
{
- surf_flux_CXX ( zeta_, Rib_, Re_, B_, z0_m_, z0_t_, Rib_conv_lim_, Cm_, Ct_, Km_, Pr_t_inv_,
+ surf_flux_esm_CXX ( zeta_, Rib_, Re_, B_, z0_m_, z0_t_, Rib_conv_lim_, Cm_, Ct_, Km_, Pr_t_inv_,
U_, dT_, Tsemi_, dQ_, h_, in_z0_m_,
*kappa, *Pr_t_0_inv, *Pr_t_inf_inv,
*alpha_m, *alpha_h, *alpha_h_fix,
@@ -25,4 +26,35 @@ void get_surface_fluxes (float *zeta_, float *Rib_, float *Re_, float *B_, float
*Pr_m, *nu_air, *g,
*maxiters_charnock, *maxiters_convection,
*grid_size);
+}
+
+void get_surface_fluxes_sheba (float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
+ 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 *alpha_m, const float *alpha_h,
+ const float *a_m, const float *a_h,
+ const float *b_m, const float *b_h,
+ const float *c_h,
+ const float *Re_rough_min,
+ const float *B1_rough, const float *B2_rough,
+ const float *B_max_land, const float *B_max_ocean, const float *B_max_lake,
+ const float *gamma_c, const float *Re_visc_min,
+ const float *Pr_m, const float *nu_air, const float *g,
+ const int *maxiters_charnock,
+ const int *grid_size)
+{
+ surf_flux_sheba_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,
+ *alpha_m, *alpha_h,
+ *a_m, *a_h,
+ *b_m, *b_h,
+ *c_h,
+ *Re_rough_min,
+ *B1_rough, *B2_rough,
+ *B_max_land, *B_max_ocean, *B_max_lake,
+ *gamma_c, *Re_visc_min,
+ *Pr_m, *nu_air, *g,
+ *maxiters_charnock,
+ *grid_size);
}
\ No newline at end of file
diff --git a/srcCU/sfx_compute_esm.cu b/srcCU/sfx_compute_esm.cu
index b9c4e2c..302494a 100644
--- a/srcCU/sfx_compute_esm.cu
+++ b/srcCU/sfx_compute_esm.cu
@@ -1,6 +1,7 @@
#include <cmath>
#include <iostream>
#include "../includeCU/sfx_compute_esm.cuh"
+#include "../includeCU/sfx_surface.cuh"
template<typename T>
__device__ void get_charnock_roughness(T &z0_m, T &u_dyn0,
@@ -45,7 +46,8 @@ template __device__ void get_charnock_roughness(double &z0_m, double &u_dyn0,
const double kappa,
const double h_charnock, const double c1_charnock, const double c2_charnock,
const int maxiters);
-
+
+
template<typename T>
__device__ void get_convection_lim(T &zeta_lim, T &Rib_lim, T &f_m_lim, T &f_h_lim,
const T h0_m, const T h0_t, const T B,
@@ -237,7 +239,7 @@ template __device__ void get_psi_semi_convection(double &psi_m, double &psi_h, d
const int maxiters);
template<typename T>
-__global__ void compute_flux_esm_gpu(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_, T *Rib_conv_lim_, T *Cm_, T *Ct_, T *Km_, T *Pr_t_inv_,
+__global__ void kernel_compute_flux_esm_gpu(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_, T *Rib_conv_lim_, T *Cm_, T *Ct_, T *Km_, T *Pr_t_inv_,
const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_, const T *h_, const T *in_z0_m_,
const T kappa, const T Pr_t_0_inv, const T Pr_t_inf_inv,
const T alpha_m, const T alpha_h, const T alpha_h_fix,
@@ -351,7 +353,7 @@ __global__ void compute_flux_esm_gpu(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_,
}
}
-template __global__ void compute_flux_esm_gpu(float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
+template __global__ void kernel_compute_flux_esm_gpu(float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
const float *U, const float *dt, const float *T_semi, const float *dq, const float *H, const float *in_z0_m,
const float kappa, const float Pr_t_0_inv, const float Pr_t_inf_inv,
const float alpha_m, const float alpha_h, const float alpha_h_fix,
@@ -362,7 +364,7 @@ template __global__ void compute_flux_esm_gpu(float *zeta_, float *Rib_, float *
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_esm_gpu(double *zeta_, double *Rib_, double *Re_, double *B_, double *z0_m_, double *z0_t_, double *Rib_conv_lim_, double *Cm_, double *Ct_, double *Km_, double *Pr_t_inv_,
+template __global__ void kernel_compute_flux_esm_gpu(double *zeta_, double *Rib_, double *Re_, double *B_, double *z0_m_, double *z0_t_, double *Rib_conv_lim_, double *Cm_, double *Ct_, double *Km_, double *Pr_t_inv_,
const double *U, const double *dt, const double *T_semi, const double *dq, const double *H, const double *in_z0_m,
const double kappa, const double Pr_t_0_inv, const double Pr_t_inf_inv,
const double alpha_m, const double alpha_h, const double alpha_h_fix,
@@ -375,7 +377,7 @@ template __global__ void compute_flux_esm_gpu(double *zeta_, double *Rib_, doubl
const int grid_size);
template<typename T>
-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_,
+void compute_flux_esm_gpu(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_, T *Rib_conv_lim_, T *Cm_, T *Ct_, T *Km_, T *Pr_t_inv_,
const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_, const T *h_, const T *in_z0_m_,
const T kappa, const T Pr_t_0_inv, const T Pr_t_inf_inv,
const T alpha_m, const T alpha_h, const T alpha_h_fix,
@@ -391,7 +393,7 @@ void compute_flux_gpu(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_, T *R
dim3 cuBlock = dim3(1024, 1, 1);
dim3 cuGrid = dim3(BlockCount, 1, 1);
- compute_flux<<<cuGrid, cuBlock>>>(zeta_, Rib_, Re_, B_, z0_m_, z0_t_, Rib_conv_lim_, Cm_, Ct_, Km_, Pr_t_inv_,
+ kernel_compute_flux_esm_gpu<<<cuGrid, cuBlock>>>(zeta_, Rib_, Re_, B_, z0_m_, z0_t_, Rib_conv_lim_, Cm_, Ct_, Km_, Pr_t_inv_,
U_, dT_, Tsemi_, dQ_, h_, in_z0_m_,
kappa, Pr_t_0_inv, Pr_t_inf_inv,
alpha_m, alpha_h, alpha_h_fix,
@@ -404,7 +406,7 @@ void compute_flux_gpu(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_, T *R
grid_size);
}
-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_,
+template void compute_flux_esm_gpu(float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
const float *U, const float *dt, const float *T_semi, const float *dq, const float *H, const float *in_z0_m,
const float kappa, const float Pr_t_0_inv, const float Pr_t_inf_inv,
const float alpha_m, const float alpha_h, const float alpha_h_fix,
@@ -415,7 +417,7 @@ template void compute_flux_gpu(float *zeta_, float *Rib_, float *Re_, float *B_,
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(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_esm_gpu(double *zeta_, double *Rib_, double *Re_, double *B_, double *z0_m_, double *z0_t_, double *Rib_conv_lim_, double *Cm_, double *Ct_, double *Km_, double *Pr_t_inv_,
const double *U, const double *dt, const double *T_semi, const double *dq, const double *H, const double *in_z0_m,
const double kappa, const double Pr_t_0_inv, const double Pr_t_inf_inv,
const double alpha_m, const double alpha_h, const double alpha_h_fix,
diff --git a/srcCU/sfx_surface.cu b/srcCU/sfx_surface.cu
new file mode 100644
index 0000000..d15f64e
--- /dev/null
+++ b/srcCU/sfx_surface.cu
@@ -0,0 +1,46 @@
+#include "../includeCU/sfx_surface.cuh"
+
+// template<typename T>
+// __device__ void get_charnock_roughness(T &z0_m, T &u_dyn0,
+// const T h, const T U,
+// const T kappa,
+// const T h_charnock, const T c1_charnock, const T c2_charnock,
+// const int maxiters)
+// {
+// T Uc, a, b, c, c_min, f;
+
+// Uc = U;
+// a = 0.0;
+// b = 25.0;
+// c_min = log(h_charnock) / kappa;
+
+// for (int i = 0; i < maxiters; i++)
+// {
+// f = c1_charnock - 2.0 * log(Uc);
+// for (int j = 0; j < maxiters; j++)
+// {
+// c = (f + 2.0 * log(b)) / kappa;
+// if (U <= 8.0e0)
+// a = log(1.0 + c2_charnock * ( pow(b / Uc, 3) ) ) / kappa;
+// c = max(c - a, c_min);
+// b = c;
+// }
+// z0_m = h_charnock * exp(-c * kappa);
+// z0_m = max(z0_m, T(0.000015e0));
+// Uc = U * log(h_charnock / z0_m) / log(h / z0_m);
+// }
+
+// u_dyn0 = Uc / c;
+// }
+
+// template __device__ void get_charnock_roughness(float &z0_m, float &u_dyn0,
+// const float h, const float U,
+// const float kappa,
+// const float h_charnock, const float c1_charnock, const float c2_charnock,
+// const int maxiters);
+// template __device__ void get_charnock_roughness(double &z0_m, double &u_dyn0,
+// const double h, const double U,
+// const double kappa,
+// const double h_charnock, const double c1_charnock, const double c2_charnock,
+// const int maxiters);
+
\ No newline at end of file
diff --git a/srcCXX/sfx_call_class_func.cpp b/srcCXX/sfx_call_class_func.cpp
index 2090d62..9bc7ac4 100644
--- a/srcCXX/sfx_call_class_func.cpp
+++ b/srcCXX/sfx_call_class_func.cpp
@@ -2,10 +2,11 @@
#include <stdio.h>
#include "../includeCXX/sfx_call_class_func.h"
#include "../includeCXX/sfx_esm.h"
+#include "../includeCXX/sfx_sheba.h"
#include <vector>
// -------------------------------------------------------------------------- //
-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_,
+void surf_flux_esm_CXX (float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
float *U_, float *dT_, float *Tsemi_, float *dQ_, float *h_, float *in_z0_m_,
const float kappa, const float Pr_t_0_inv, const float Pr_t_inf_inv,
const float alpha_m, const float alpha_h, const float alpha_h_fix,
@@ -31,7 +32,44 @@ void surf_flux_CXX (float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_
gamma_c, Re_visc_min,
Pr_m, nu_air, g);
- F.compute_flux(zeta_, Rib_, Re_, B_, z0_m_, z0_t_, Rib_conv_lim_, Cm_, Ct_, Km_, Pr_t_inv_,
+ F.compute_flux_esm(zeta_, Rib_, Re_, B_, z0_m_, z0_t_, Rib_conv_lim_, Cm_, Ct_, Km_, Pr_t_inv_,
U_, dT_, Tsemi_, dQ_, h_, in_z0_m_,
maxiters_charnock, maxiters_convection);
+}
+
+void surf_flux_sheba_CXX (float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
+ 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 alpha_m, const float alpha_h,
+ const float a_m, const float a_h,
+ const float b_m, const float b_h,
+ const float c_h,
+ const float Re_rough_min,
+ const float B1_rough, const float B2_rough,
+ const float B_max_land, const float B_max_ocean, const float B_max_lake,
+ const float gamma_c, const float Re_visc_min,
+ const float Pr_m, const float nu_air, const float g,
+ const int maxiters_charnock,
+ const int grid_size)
+{
+#ifdef INCLUDE_CUDA
+ static FluxSheba<float, MemType::CPU, MemType::CPU, MemType::GPU> F;
+#else
+ static FluxSheba<float, MemType::CPU, MemType::CPU, MemType::CPU> F;
+#endif
+
+ F.set_params(grid_size, kappa, Pr_t_0_inv,
+ alpha_m, alpha_h,
+ a_m, a_h,
+ b_m, b_h,
+ c_h,
+ Re_rough_min,
+ B1_rough, B2_rough,
+ B_max_land, B_max_ocean, B_max_lake,
+ gamma_c, Re_visc_min,
+ Pr_m, nu_air, g);
+
+ F.compute_flux_sheba(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);
}
\ No newline at end of file
diff --git a/srcCXX/sfx_compute_esm.cpp b/srcCXX/sfx_compute_esm.cpp
index b50f49d..9a30be0 100644
--- a/srcCXX/sfx_compute_esm.cpp
+++ b/srcCXX/sfx_compute_esm.cpp
@@ -1,50 +1,7 @@
#include <cmath>
#include <iostream>
#include "../includeCXX/sfx_compute_esm.h"
-
-template<typename T>
-void get_charnock_roughness(T &z0_m, T &u_dyn0,
- const T h, const T U,
- const T kappa,
- const T h_charnock, const T c1_charnock, const T c2_charnock,
- const int maxiters)
-{
- T Uc, a, b, c, c_min, f;
-
- Uc = U;
- a = 0.0;
- b = 25.0;
- c_min = log(h_charnock) / kappa;
-
- for (int i = 0; i < maxiters; i++)
- {
- f = c1_charnock - 2.0 * log(Uc);
- for (int j = 0; j < maxiters; j++)
- {
- c = (f + 2.0 * log(b)) / kappa;
- if (U <= 8.0e0)
- a = log(1.0 + c2_charnock * ( pow(b / Uc, 3) ) ) / kappa;
- c = std::max(c - a, c_min);
- b = c;
- }
- z0_m = h_charnock * exp(-c * kappa);
- z0_m = std::max(z0_m, T(0.000015e0));
- Uc = U * log(h_charnock / z0_m) / log(h / z0_m);
- }
-
- u_dyn0 = Uc / c;
-}
-
-template void get_charnock_roughness(float &z0_m, float &u_dyn0,
- const float h, const float U,
- const float kappa,
- const float h_charnock, const float c1_charnock, const float c2_charnock,
- const int maxiters);
-template void get_charnock_roughness(double &z0_m, double &u_dyn0,
- const double h, const double U,
- const double kappa,
- const double h_charnock, const double c1_charnock, const double c2_charnock,
- const int maxiters);
+#include "../includeCXX/sfx_surface.h"
template<typename T>
void get_convection_lim(T &zeta_lim, T &Rib_lim, T &f_m_lim, T &f_h_lim,
diff --git a/srcCXX/sfx_compute_sheba.cpp b/srcCXX/sfx_compute_sheba.cpp
index ad420c0..1aa7252 100644
--- a/srcCXX/sfx_compute_sheba.cpp
+++ b/srcCXX/sfx_compute_sheba.cpp
@@ -1,18 +1,221 @@
#include <cmath>
#include <iostream>
#include "../includeCXX/sfx_compute_sheba.h"
+#include "../includeCXX/sfx_surface.h"
+
+template<typename T>
+void get_psi_mh(T &psi_m, T &psi_h,
+ const T zeta_m, const T zeta_h,
+ const T alpha_m, const T alpha_h,
+ const T a_m, const T a_h,
+ const T b_m, const T b_h,
+ const T c_h)
+{
+ T x_m, x_h;
+ T q_m, q_h;
+
+ if (zeta_m >= 0.0)
+ {
+ q_m = pow((1.0 - b_m) / b_m, 1.0 / 3.0);
+ x_m = pow(1.0 + zeta_m, 1.0 / 3.0);
+
+ psi_m = -3.0 * (a_m / b_m) * (x_m - 1.0) + 0.5 * (a_m / 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))));
+ }
+ else
+ { x_m = pow(1.0 - 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);
+ }
+
+ if (zeta_h >= 0.0)
+ {
+ q_h = sqrt(c_h * c_h - 4.0);
+ x_h = zeta_h;
+
+ psi_h = -0.5 * b_h * log(1.0 + c_h * x_h + x_h * x_h) + ((-a_h / q_h) + ((b_h * c_h) / (2.0 * q_h))) * (log((2.0 * x_h + c_h - q_h) / (2.0 * x_h + c_h + q_h)) - log((c_h - q_h) / (c_h + q_h)));
+ }
+ else
+ {
+ x_h = pow(1.0 - alpha_h * zeta_h, 0.25);
+ psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h));
+ }
+}
+
+template void get_psi_mh(float &psi_m, float &psi_h,
+ const float zeta_m, const float zeta_h,
+ const float alpha_m, const float alpha_h,
+ const float a_m, const float a_h,
+ const float b_m, const float b_h,
+ const float c_h);
+template void get_psi_mh(double &psi_m, double &psi_h,
+ const double zeta_m, const double zeta_h,
+ const double alpha_m, const double alpha_h,
+ const double a_m, const double a_h,
+ const double b_m, const double b_h,
+ const double c_h);
+
+template<typename T>
+void get_psi(T &psi_m, T &psi_h,
+ const T zeta,
+ const T alpha_m, const T alpha_h,
+ const T a_m, const T a_h,
+ const T b_m, const T b_h,
+ const T c_h)
+{
+ T x_m, x_h;
+ T q_m, q_h;
+
+ if (zeta >= 0.0)
+ {
+ q_m = pow((1.0 - b_m) / b_m, 1.0 / 3.0);
+ q_h = sqrt(c_h * c_h - 4.0);
+
+ x_m = pow(1.0 + zeta, 1.0 / 3.0);
+ x_h = zeta;
+
+ psi_m = -3.0 * (a_m / b_m) * (x_m - 1.0) + 0.5 * (a_m / 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_h = -0.5 * b_h * log(1.0 + c_h * x_h + x_h * x_h) + ((-a_h / q_h) + ((b_h * c_h) / (2.0 * q_h))) * (log((2.0 * x_h + c_h - q_h) / (2.0 * x_h + c_h + q_h)) - log((c_h - q_h) / (c_h + q_h)));
+ }
+ else
+ {
+ x_m = pow(1.0 - alpha_m * zeta, 0.25);
+ x_h = pow(1.0 - 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));
+ }
+}
+
+template void get_psi(float &psi_m, float &psi_h,
+ const float zeta,
+ const float alpha_m, const float alpha_h,
+ const float a_m, const float a_h,
+ const float b_m, const float b_h,
+ const float c_h);
+template void get_psi(double &psi_m, double &psi_h,
+ const double zeta,
+ const double alpha_m, const double alpha_h,
+ const double a_m, const double a_h,
+ const double b_m, const double b_h,
+ const double c_h);
+
+template<typename T>
+void get_dynamic_scales(T &Udyn, T &Tdyn, T &Qdyn, T &zeta,
+ const T U, const T Tsemi, const T dT, const T dQ, const T z, const T z0_m, const T z0_t, const T beta,
+ const T kappa, const T Pr_t_0_inv,
+ const T alpha_m, const T alpha_h,
+ const T a_m, const T a_h,
+ const T b_m, const T b_h,
+ const T c_h,
+ const int maxiters)
+{
+ T psi_m, psi_h, psi0_m, psi0_h, Linv;
+ const T gamma = 0.61;
+
+ Udyn = kappa * U / log(z / z0_m);
+ Tdyn = kappa * dT * Pr_t_0_inv / log(z / z0_t);
+ Qdyn = kappa * dQ * Pr_t_0_inv / log(z / z0_t);
+ zeta = 0.0;
+
+ // --- no wind
+ if (Udyn < 1e-5)
+ return;
+
+ Linv = kappa * beta * (Tdyn + gamma * Qdyn * Tsemi) / (Udyn * Udyn);
+ zeta = z * Linv;
+
+ // --- near neutral case
+ if (Linv < 1e-5)
+ return;
+
+ for (int i = 0; i < maxiters; i++)
+ {
+ get_psi(psi_m, psi_h, zeta, alpha_m, alpha_h,
+ a_m, a_h,
+ b_m, b_h,
+ c_h);
+
+ get_psi_mh(psi0_m, psi0_h, z0_m * Linv, z0_t * Linv,
+ alpha_m, alpha_h,
+ a_m, a_h,
+ b_m, b_h,
+ c_h);
+
+ Udyn = kappa * U / (log(z / z0_m) - (psi_m - psi0_m));
+ Tdyn = kappa * dT * Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h));
+ Qdyn = kappa * dQ * Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h));
+
+ if (Udyn < 1e-5) break;
+
+ Linv = kappa * beta * (Tdyn + gamma * Qdyn * Tsemi) / (Udyn * Udyn);
+ zeta = z * Linv;
+ }
+}
+
+template void get_dynamic_scales(float &Udyn, float &Tdyn, float &Qdyn, float & zeta,
+ const float U, const float Tsemi, const float dT, const float dQ, const float z, const float z0_m, const float z0_t, const float beta,
+ const float kappa, const float Pr_t_0_inv,
+ const float alpha_m, const float alpha_h,
+ const float a_m, const float a_h,
+ const float b_m, const float b_h,
+ const float c_h,
+ const int maxiters);
+template void get_dynamic_scales(double &Udyn, double &Tdyn, double &Qdyn, double & zeta,
+ const double U, const double Tsemi, const double dT, const double dQ, const double z, const double z0_m, const double z0_t, const double beta,
+ const double kappa, const double Pr_t_0_inv,
+ const double alpha_m, const double alpha_h,
+ const double a_m, const double a_h,
+ const double b_m, const double b_h,
+ const double c_h,
+ const int maxiters);
+
+template<typename T>
+void get_phi(T &phi_m, T &phi_h,
+ const T zeta,
+ const T alpha_m, const T alpha_h,
+ const T a_m, const T a_h,
+ const T b_m, const T b_h,
+ const T c_h)
+{
+ if (zeta >= 0.0)
+ {
+ phi_m = 1.0 + (a_m * zeta * pow(1.0 + zeta, 1.0 / 3.0) ) / (1.0 + b_m * zeta);
+ phi_h = 1.0 + (a_h * zeta + b_h * zeta * zeta) / (1.0 + c_h * zeta + zeta * zeta);
+ }
+ else
+ {
+ phi_m = pow(1.0 - alpha_m * zeta, -0.25);
+ phi_h = pow(1.0 - alpha_h * zeta, -0.5);
+ }
+}
+
+template void get_phi(float &phi_m, float &phi_h,
+ const float zeta,
+ const float alpha_m, const float alpha_h,
+ const float a_m, const float a_h,
+ const float b_m, const float b_h,
+ const float c_h);
+template void get_phi(double &phi_m, double &phi_h,
+ const double zeta,
+ const double alpha_m, const double alpha_h,
+ const double a_m, const double a_h,
+ const double b_m, const double b_h,
+ const double c_h);
template<typename T>
void compute_flux_sheba_cpu(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_, T *Rib_conv_lim_, T *Cm_, T *Ct_, T *Km_, T *Pr_t_inv_,
const T *U_, const T *dT_, const T *Tsemi_, const T *dQ_, const T *h_, const T *in_z0_m_,
- const T kappa, const T Pr_t_0_inv, const T Pr_t_inf_inv,
- const T alpha_m, const T alpha_h, const T alpha_h_fix,
- const T beta_m, const T beta_h, const T Rib_max, const T Re_rough_min,
+ const T kappa, const T Pr_t_0_inv,
+ const T alpha_m, const T alpha_h,
+ const T a_m, const T a_h,
+ const T b_m, const T b_h,
+ const T c_h,
+ const T Re_rough_min,
const T B1_rough, const T B2_rough,
const T B_max_land, const T B_max_ocean, const T B_max_lake,
const T gamma_c, const T Re_visc_min,
const T Pr_m, const T nu_air, const T g,
- const int maxiters_charnock, const int maxiters_convection,
+ const int maxiters_charnock,
const int grid_size)
{
T h, U, dT, Tsemi, dQ, z0_m;
@@ -20,6 +223,9 @@ void compute_flux_sheba_cpu(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_
zeta, Rib, Udyn, Tdyn, Qdyn, phi_m, phi_h,
Km, Pr_t_inv, Cm, Ct;
+ const T B3_rough = kappa * Pr_m, B4_rough =(0.14 * (pow(30.0, B2_rough))) * (pow(Pr_m, 0.8));
+ const T h_charnock = 10.0, c1_charnock = log(h_charnock * (g / gamma_c)), c2_charnock = Re_visc_min * nu_air * c1_charnock;
+
int surface_type;
for (int step = 0; step < grid_size; step++)
@@ -34,19 +240,19 @@ void compute_flux_sheba_cpu(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_
if (z0_m < 0.0) surface_type = 0;
else surface_type = 1;
- if (surface_type == surface_ocean)
+ if (surface_type == 0)
{
- get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%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 == surface_land)
+ if (surface_type == 1)
{
h0_m = h / z0_m;
u_dyn0 = U * kappa / log(h0_m);
}
Re = u_dyn0 * z0_m / nu_air;
- get_thermal_roughness(z0_t, B, z0_m, Re, surface_type);
+ get_thermal_roughness(z0_t, B, z0_m, Re, Re_rough_min, B1_rough, B2_rough, B3_rough, B4_rough, B_max_ocean, B_max_lake, B_max_land, surface_type);
// --- define relative height [thermal]
h0_t = h / z0_t;
@@ -56,10 +262,10 @@ void compute_flux_sheba_cpu(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_
// --- get the fluxes
// ----------------------------------------------------------------------------
- get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, U, Tsemi, dT, dQ, h, z0_m, z0_t, (g / Tsemi), 10);
+ get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, U, Tsemi, dT, dQ, h, z0_m, z0_t, (g / Tsemi), kappa, Pr_t_0_inv, alpha_m, alpha_h, a_m, a_h, b_m, b_h, c_h, 10);
// ----------------------------------------------------------------------------
- get_phi(phi_m, phi_h, zeta);
+ get_phi(phi_m, phi_h, zeta, alpha_m, alpha_h, a_m, a_h, b_m, b_h, c_h);
// ----------------------------------------------------------------------------
// --- define transfer coeff. (momentum) & (heat)
@@ -80,7 +286,7 @@ void compute_flux_sheba_cpu(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_
B_[step] = B;
z0_m_[step] = z0_m;
z0_t_[step] = z0_t;
- Rib_conv_lim_[step] = Rib_conv_lim;
+ Rib_conv_lim_[step] = 0.0;
Cm_[step] = Cm;
Ct_[step] = Ct;
Km_[step] = Km;
@@ -90,23 +296,29 @@ void compute_flux_sheba_cpu(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_
template void compute_flux_sheba_cpu(float *zeta_, float *Rib_, float *Re_, float *B_, float *z0_m_, float *z0_t_, float *Rib_conv_lim_, float *Cm_, float *Ct_, float *Km_, float *Pr_t_inv_,
const float *U, const float *dt, const float *T_semi, const float *dq, const float *H, const float *in_z0_m,
- const float kappa, const float Pr_t_0_inv, const float Pr_t_inf_inv,
- const float alpha_m, const float alpha_h, const float alpha_h_fix,
- const float beta_m, const float beta_h, const float Rib_max, const float Re_rough_min,
+ const float kappa, const float Pr_t_0_inv,
+ const float alpha_m, const float alpha_h,
+ const float a_m, const float a_h,
+ const float b_m, const float b_h,
+ const float c_h,
+ const float Re_rough_min,
const float B1_rough, const float B2_rough,
const float B_max_land, const float B_max_ocean, const float B_max_lake,
const float gamma_c, const float Re_visc_min,
const float Pr_m, const float nu_air, const float g,
- const int maxiters_charnock, const int maxiters_convection,
+ const int maxiters_charnock,
const int grid_size);
template void compute_flux_sheba_cpu(double *zeta_, double *Rib_, double *Re_, double *B_, double *z0_m_, double *z0_t_, double *Rib_conv_lim_, double *Cm_, double *Ct_, double *Km_, double *Pr_t_inv_,
const double *U, const double *dt, const double *T_semi, const double *dq, const double *H, const double *in_z0_m,
- const double kappa, const double Pr_t_0_inv, const double Pr_t_inf_inv,
- const double alpha_m, const double alpha_h, const double alpha_h_fix,
- const double beta_m, const double beta_h, const double Rib_max, const double Re_rough_min,
+ const double kappa, const double Pr_t_0_inv,
+ const double alpha_m, const double alpha_h,
+ const double a_m, const double a_h,
+ const double b_m, const double b_h,
+ const double c_h,
+ const double Re_rough_min,
const double B1_rough, const double B2_rough,
const double B_max_land, const double B_max_ocean, const double B_max_lake,
const double gamma_c, const double Re_visc_min,
const double Pr_m, const double nu_air, const double g,
- const int maxiters_charnock, const int maxiters_convection,
+ const int maxiters_charnock,
const int grid_size);
\ No newline at end of file
diff --git a/srcCXX/sfx_esm.cpp b/srcCXX/sfx_esm.cpp
index 5564c28..5d0a9de 100644
--- a/srcCXX/sfx_esm.cpp
+++ b/srcCXX/sfx_esm.cpp
@@ -183,7 +183,7 @@ void FluxEsm<T, memIn, memOut, RunMem>::set_data(T *zeta_, T *Rib_, T *Re_, T *B
}
template<typename T, MemType memIn, MemType memOut, MemType RunMem >
-void FluxEsm<T, memIn, memOut, RunMem>::compute_FluxEsm(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 FluxEsm<T, memIn, memOut, RunMem>::compute_flux_esm(T *zeta_, T *Rib_, T *Re_, T *B_, T *z0_m_, T *z0_t_, T *Rib_conv_lim_, T *Cm_, T *Ct_, T *Km_, T *Pr_t_inv_,
T *U_, T *dT_, T *Tsemi_, T *dQ_, T *h_, T *in_z0_m_,
const int maxiters_charnock, const int maxiters_convection)
{
diff --git a/srcCXX/sfx_sheba.cpp b/srcCXX/sfx_sheba.cpp
index 5edfa73..66d9a81 100644
--- a/srcCXX/sfx_sheba.cpp
+++ b/srcCXX/sfx_sheba.cpp
@@ -12,10 +12,12 @@
template<typename T, MemType memIn, MemType memOut, MemType RunMem >
FluxSheba<T, memIn, memOut, RunMem>::FluxSheba()
{
- kappa = 0, Pr_t_0_inv = 0, Pr_t_inf_inv = 0,
- alpha_m = 0, alpha_h = 0, alpha_h_fix = 0,
- beta_m = 0, beta_h = 0,
- Rib_max = 0, Re_rough_min = 0,
+ kappa = 0, Pr_t_0_inv = 0,
+ alpha_m = 0, alpha_h = 0,
+ a_m = 0, a_h = 0,
+ b_m = 0, b_h = 0,
+ c_h = 0,
+ Re_rough_min = 0,
B1_rough = 0, B2_rough = 0,
B_max_land = 0, B_max_ocean = 0, B_max_lake = 0,
gamma_c = 0,
@@ -29,18 +31,23 @@ FluxSheba<T, memIn, memOut, RunMem>::FluxSheba()
}
template<typename T, MemType memIn, MemType memOut, MemType RunMem >
-void FluxSheba<T, memIn, memOut, RunMem>::set_params(const int grid_size_, const T kappa_, const T Pr_t_0_inv_, const T Pr_t_inf_inv_,
- const T alpha_m_, const T alpha_h_, const T alpha_h_fix_,
- const T beta_m_, const T beta_h_, const T Rib_max_, const T Re_rough_min_,
+void FluxSheba<T, memIn, memOut, RunMem>::set_params(const int grid_size_, const T kappa_, const T Pr_t_0_inv_,
+ const T alpha_m_, const T alpha_h_,
+ const float a_m_, const float a_h_,
+ const float b_m_, const float b_h_,
+ const float c_h_,
+ const T Re_rough_min_,
const T B1_rough_, const T B2_rough_,
const T B_max_land_, const T B_max_ocean_, const T B_max_lake_,
const T gamma_c_, const T Re_visc_min_,
const T Pr_m_, const T nu_air_, const T g_)
{
- kappa = kappa_, Pr_t_0_inv = Pr_t_0_inv_, Pr_t_inf_inv = Pr_t_inf_inv_,
- alpha_m = alpha_m_, alpha_h = alpha_h_, alpha_h_fix = alpha_h_fix_,
- beta_m = beta_m_, beta_h = beta_h_,
- Rib_max = Rib_max_, Re_rough_min = Re_rough_min_, B_max_lake = B_max_lake_,
+ kappa = kappa_, Pr_t_0_inv = Pr_t_0_inv_,
+ alpha_m = alpha_m_, alpha_h = alpha_h_,
+ a_m = a_m_, a_h = a_h_,
+ b_m = b_m_, b_h = b_h_,
+ c_h = c_h_,
+ Re_rough_min = Re_rough_min_, B_max_lake = B_max_lake_,
B1_rough = B1_rough_, B2_rough = B2_rough_,
B_max_land = B_max_land_, B_max_ocean = B_max_ocean_,
gamma_c = gamma_c_,
@@ -99,10 +106,12 @@ void FluxSheba<T, memIn, memOut, RunMem>::set_params(const int grid_size_, const
template<typename T, MemType memIn, MemType memOut, MemType RunMem >
FluxSheba<T, memIn, memOut, RunMem>::~FluxSheba()
{
- kappa = 0, Pr_t_0_inv = 0, Pr_t_inf_inv = 0,
- alpha_m = 0, alpha_h = 0, alpha_h_fix = 0,
- beta_m = 0, beta_h = 0,
- Rib_max = 0, Re_rough_min = 0,
+ kappa = 0, Pr_t_0_inv = 0,
+ alpha_m = 0, alpha_h = 0,
+ a_m = 0, a_h = 0,
+ b_m = 0, b_h = 0,
+ c_h = 0,
+ Re_rough_min = 0,
B1_rough = 0, B2_rough = 0,
B_max_land = 0, B_max_ocean = 0, B_max_lake = 0,
gamma_c = 0,
@@ -183,37 +192,43 @@ void FluxSheba<T, memIn, memOut, RunMem>::set_data(T *zeta_, T *Rib_, T *Re_, T
}
template<typename T, MemType memIn, MemType memOut, MemType RunMem >
-void FluxSheba<T, memIn, memOut, RunMem>::compute_FluxSheba(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 FluxSheba<T, memIn, memOut, RunMem>::compute_flux_sheba(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)
+ const int maxiters_charnock)
{
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_FluxSheba_cpu(zeta, Rib, Re, B, z0_m, z0_t, Rib_conv_lim, Cm, Ct, Km, Pr_t_inv,
+ if(RunMem == MemType::CPU) compute_flux_sheba_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,
+ kappa, Pr_t_0_inv,
+ alpha_m, alpha_h,
+ a_m, a_h,
+ b_m, b_h,
+ c_h,
+ Re_rough_min,
B1_rough, B2_rough,
B_max_land, B_max_ocean, B_max_lake,
gamma_c, Re_visc_min,
Pr_m, nu_air, g,
- maxiters_charnock, maxiters_convection,
+ maxiters_charnock,
grid_size);
#ifdef INCLUDE_CUDA
- else compute_FluxSheba_gpu(zeta, Rib, Re, B, z0_m, z0_t, Rib_conv_lim, Cm, Ct, Km, Pr_t_inv,
- U, dT, Tsemi, dQ, h, in_z0_m,
- kappa, Pr_t_0_inv, Pr_t_inf_inv,
- alpha_m, alpha_h, alpha_h_fix,
- beta_m, beta_h, Rib_max, Re_rough_min,
- B1_rough, B2_rough,
- B_max_land, B_max_ocean, B_max_lake,
- gamma_c, Re_visc_min,
- Pr_m, nu_air, g,
- maxiters_charnock, maxiters_convection,
- grid_size);
+ // else compute_flux_sheba_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,
+ // alpha_m, alpha_h,
+ // a_m, a_h,
+ // b_m, b_h,
+ // c_h,
+ // Re_rough_min,
+ // B1_rough, B2_rough,
+ // B_max_land, B_max_ocean, B_max_lake,
+ // gamma_c, Re_visc_min,
+ // Pr_m, nu_air, g,
+ // maxiters_charnock,
+ // grid_size);
#endif
if(RunMem != memOut)
diff --git a/srcCXX/sfx_surface.cpp b/srcCXX/sfx_surface.cpp
new file mode 100644
index 0000000..2bb012c
--- /dev/null
+++ b/srcCXX/sfx_surface.cpp
@@ -0,0 +1,87 @@
+#include "../includeCXX/sfx_surface.h"
+#include <cmath>
+#include <iostream>
+
+template<typename T>
+void get_charnock_roughness(T &z0_m, T &u_dyn0,
+ const T h, const T U,
+ const T kappa,
+ const T h_charnock, const T c1_charnock, const T c2_charnock,
+ const int maxiters)
+{
+ T Uc, a, b, c, c_min, f;
+
+ Uc = U;
+ a = 0.0;
+ b = 25.0;
+ c_min = log(h_charnock) / kappa;
+
+ for (int i = 0; i < maxiters; i++)
+ {
+ f = c1_charnock - 2.0 * log(Uc);
+ for (int j = 0; j < maxiters; j++)
+ {
+ c = (f + 2.0 * log(b)) / kappa;
+ if (U <= 8.0e0)
+ a = log(1.0 + c2_charnock * ( pow(b / Uc, 3) ) ) / kappa;
+ c = std::max(c - a, c_min);
+ b = c;
+ }
+ z0_m = h_charnock * exp(-c * kappa);
+ z0_m = std::max(z0_m, T(0.000015e0));
+ Uc = U * log(h_charnock / z0_m) / log(h / z0_m);
+ }
+
+ u_dyn0 = Uc / c;
+}
+
+template void get_charnock_roughness(float &z0_m, float &u_dyn0,
+ const float h, const float U,
+ const float kappa,
+ const float h_charnock, const float c1_charnock, const float c2_charnock,
+ const int maxiters);
+template void get_charnock_roughness(double &z0_m, double &u_dyn0,
+ const double h, const double U,
+ const double kappa,
+ const double h_charnock, const double c1_charnock, const double c2_charnock,
+ const int maxiters);
+
+template<typename T>
+void get_thermal_roughness(T &z0_t, T &B,
+ const T z0_m, const T Re,
+ const T Re_rough_min,
+ const T B1_rough, const T B2_rough, const T B3_rough, const T B4_rough,
+ const T B_max_ocean, const T B_max_lake, const T B_max_land,
+ const int surface_type)
+{
+ // --- define B = log(z0_m / z0_t)
+ if (Re <= Re_rough_min)
+ B = B1_rough * log(B3_rough * Re) + B2_rough;
+ else
+ // *: B4 takes into account Re value at z' ~ O(10) z0
+ B = B4_rough * (pow(Re, B2_rough));
+
+ // --- apply max restriction based on surface type
+ if (surface_type == 0)
+ B = std::min(B, B_max_ocean);
+ else if (surface_type == 1)
+ B = std::min(B, B_max_lake);
+ else if (surface_type == 2)
+ B = std::min(B, B_max_land);
+
+ // --- define roughness [thermal]
+ z0_t = z0_m / exp(B);
+}
+
+template void get_thermal_roughness(float &z0_t, float &B,
+ const float z0_m, const float Re,
+ const float Re_rough_min,
+ const float B1_rough, const float B2_rough, const float B3_rough, const float B4_rough,
+ const float B_max_ocean, const float B_max_lake, const float B_max_land,
+ const int surface_type);
+template void get_thermal_roughness(double &z0_t, double &B,
+ const double z0_m, const double Re,
+ const double Re_rough_min,
+ const double B1_rough, const double B2_rough, const double B3_rough, const double B4_rough,
+ const double B_max_ocean, const double B_max_lake, const double B_max_land,
+ const int surface_type);
\ No newline at end of file
diff --git a/srcF/sfx_esm.f90 b/srcF/sfx_esm.f90
index 4b37e91..5c540cc 100644
--- a/srcF/sfx_esm.f90
+++ b/srcF/sfx_esm.f90
@@ -55,7 +55,7 @@ contains
integer i
! ----------------------------------------------------------------------------
#if defined(INCLUDE_CUDA) || defined(INCLUDE_CXX)
- call get_surface_fluxes(sfx%zeta, sfx%Rib, sfx%Re, sfx%B, sfx%z0_m, sfx%z0_t, &
+ call get_surface_fluxes_esm(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, &
@@ -86,7 +86,6 @@ contains
end subroutine get_surface_fluxes_vec
! --------------------------------------------------------------------------------
-#if !defined(INCLUDE_CUDA) && !defined(INCLUDE_CXX)
! --------------------------------------------------------------------------------
subroutine get_surface_fluxes(sfx, meteo, numerics)
!< @brief surface flux calculation for single cell
@@ -261,7 +260,7 @@ contains
end subroutine get_surface_fluxes
! --------------------------------------------------------------------------------
-#endif
+
! convection universal functions shortcuts
! --------------------------------------------------------------------------------
function f_m_conv(zeta)
diff --git a/srcF/sfx_fc_wrapper.F90 b/srcF/sfx_fc_wrapper.F90
index 4e23a17..8a9bea8 100644
--- a/srcF/sfx_fc_wrapper.F90
+++ b/srcF/sfx_fc_wrapper.F90
@@ -1,7 +1,7 @@
module C_FUNC
INTERFACE
#if defined(INCLUDE_CUDA) || defined(INCLUDE_CXX)
- SUBROUTINE get_surface_fluxes(zeta, Rib, Re, B, z0_m, z0_t, Rib_conv_lim, &
+ SUBROUTINE get_surface_fluxes_esm(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, &
@@ -22,7 +22,33 @@ module C_FUNC
B_max_lake, gamma_c, Re_visc_min, Pr_m, nu_air, g
REAL(C_FLOAT), dimension(grid_size) :: U, dT, Tsemi, dQ, h, in_z0_m, zeta, Rib, Re, &
Rib_conv_lim, z0_m, z0_t, B, Cm, Ct, Km, Prt_inv
- END SUBROUTINE get_surface_fluxes
+ END SUBROUTINE get_surface_fluxes_esm
+
+ SUBROUTINE get_surface_fluxes_sheba(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, &
+ alpha_m, alpha_h, &
+ a_m, a_h, &
+ b_m, b_h, &
+ c_h, &
+ Re_rough_min, &
+ B1_rough, B2_rough, &
+ B_max_land, B_max_ocean, B_max_lake, &
+ gamma_c, Re_visc_min, &
+ Pr_m, nu_air, g, &
+ maxiters_charnock, &
+ grid_size) BIND(C)
+ USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_INT
+ USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_FLOAT
+ IMPLICIT NONE
+ INTEGER(C_INT) :: grid_size, maxiters_charnock
+ REAL(C_FLOAT) :: kappa, Pr_t_0_inv, Pr_t_inf_inv, alpha_m, alpha_h, a_m, a_h, b_m, b_h, &
+ c_h, Rib_max, Re_rough_min, B1_rough, B2_rough, B_max_land, B_max_ocean, &
+ B_max_lake, gamma_c, Re_visc_min, Pr_m, nu_air, g
+ REAL(C_FLOAT), dimension(grid_size) :: U, dT, Tsemi, dQ, h, in_z0_m, zeta, Rib, Re, &
+ Rib_conv_lim, z0_m, z0_t, B, Cm, Ct, Km, Prt_inv
+ END SUBROUTINE get_surface_fluxes_sheba
#endif
END INTERFACE
end module C_FUNC
--
GitLab