diff --git a/srcF/sfx_sheba_noniterative.f90 b/srcF/sfx_sheba_noniterative.f90
index c4bbc50c1fdb7a5b5672964eb45eec7462d571da..41e0dae3dcef8f568f803b969d5770685ba11d71 100644
--- a/srcF/sfx_sheba_noniterative.f90
+++ b/srcF/sfx_sheba_noniterative.f90
@@ -1,11 +1,6 @@
#include "../includeF/sfx_def.fi"
module sfx_sheba_noniterative
-<<<<<<< HEAD
- !< @brief main Earth System Model surface flux module
-=======
- !< @brief SHEBA surface flux module
->>>>>>> 9d99a415378a2907d460477f87825d027fae071e
! modules used
! --------------------------------------------------------------------------------
@@ -14,12 +9,8 @@ module sfx_sheba_noniterative
#endif
use sfx_data
use sfx_surface
-<<<<<<< HEAD
use sfx_sheba_noit_param
-=======
- use sfx_sheba_param
->>>>>>> 9d99a415378a2907d460477f87825d027fae071e
#if defined(INCLUDE_CXX)
use iso_c_binding, only: C_LOC, C_PTR, C_INT, C_FLOAT
use C_FUNC
@@ -36,24 +27,16 @@ module sfx_sheba_noniterative
! --------------------------------------------------------------------------------
public :: get_surface_fluxes
public :: get_surface_fluxes_vec
-<<<<<<< HEAD
-=======
- public :: get_psi
->>>>>>> 9d99a415378a2907d460477f87825d027fae071e
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
type, public :: numericsType
-<<<<<<< HEAD
integer :: maxiters_convection = 10 !< maximum (actual) number of iterations in convection
-=======
->>>>>>> 9d99a415378a2907d460477f87825d027fae071e
integer :: maxiters_charnock = 10 !< maximum (actual) number of iterations in charnock roughness
end type
! --------------------------------------------------------------------------------
#if defined(INCLUDE_CXX)
-<<<<<<< HEAD
type, BIND(C), public :: sfx_sheba_noit_param_C
real(C_FLOAT) :: kappa
real(C_FLOAT) :: Pr_t_0_inv
@@ -68,28 +51,10 @@ module sfx_sheba_noniterative
end type
type, BIND(C), public :: sfx_sheba_noit_numericsType_C
- integer(C_INT) :: maxiters_convection
-=======
- type, BIND(C), public :: sfx_sheba_param_C
- real(C_FLOAT) :: kappa
- real(C_FLOAT) :: Pr_t_0_inv
-
- real(C_FLOAT) :: alpha_m
- real(C_FLOAT) :: alpha_h
- real(C_FLOAT) :: a_m
- real(C_FLOAT) :: b_m
- real(C_FLOAT) :: a_h
- real(C_FLOAT) :: b_h
- real(C_FLOAT) :: c_h
- end type
-
- type, BIND(C), public :: sfx_sheba_numericsType_C
->>>>>>> 9d99a415378a2907d460477f87825d027fae071e
integer(C_INT) :: maxiters_charnock
end type
INTERFACE
-<<<<<<< HEAD
SUBROUTINE c_sheba_noit_compute_flux(sfx, meteo, model_param, surface_param, numerics, constants, grid_size) BIND(C, &
name="c_sheba_noit_compute_flux")
use sfx_data
@@ -104,7 +69,7 @@ module sfx_sheba_noniterative
type(sfx_sheba_noit_numericsType_C) :: numerics
type(sfx_phys_constants) :: constants
END SUBROUTINE c_sheba_noit_compute_flux
-=======
+
SUBROUTINE c_sheba_compute_flux(sfx, meteo, model_param, surface_param, numerics, constants, grid_size) BIND(C, &
name="c_sheba_compute_flux")
use sfx_data
@@ -119,13 +84,12 @@ module sfx_sheba_noniterative
type(sfx_sheba_numericsType_C) :: numerics
type(sfx_phys_constants) :: constants
END SUBROUTINE c_sheba_compute_flux
->>>>>>> 9d99a415378a2907d460477f87825d027fae071e
+
END INTERFACE
#endif
contains
-<<<<<<< HEAD
! --------------------------------------------------------------------------------
#if defined(INCLUDE_CXX)
subroutine set_c_struct_sfx_sheba_noit_param_values(sfx_model_param)
@@ -143,7 +107,6 @@ contains
end subroutine set_c_struct_sfx_sheba_noit_param_values
#endif
-=======
#if defined(INCLUDE_CXX)
subroutine set_c_struct_sfx_sheba_param_values(sfx_model_param)
type (sfx_sheba_param_C), intent(inout) :: sfx_model_param
@@ -161,7 +124,6 @@ contains
#endif
! --------------------------------------------------------------------------------
->>>>>>> 9d99a415378a2907d460477f87825d027fae071e
subroutine get_surface_fluxes_vec(sfx, meteo, numerics, n)
!< @brief surface flux calculation for array data
!< @details contains C/C++ & CUDA interface
@@ -182,27 +144,24 @@ contains
type (meteoDataVecTypeC), target :: meteo_c !< meteorological data (input)
type (sfxDataVecTypeC), target :: sfx_c !< surface flux data (output)
type(C_PTR) :: meteo_c_ptr, sfx_c_ptr
-<<<<<<< HEAD
type (sfx_sheba_noit_param_C) :: model_param
type (sfx_surface_param) :: surface_param
type (sfx_sheba_noit_numericsType_C) :: numerics_c
type (sfx_phys_constants) :: phys_constants
numerics_c%maxiters_convection = numerics%maxiters_convection
-=======
+
type (sfx_sheba_param_C) :: model_param
type (sfx_surface_param) :: surface_param
type (sfx_sheba_numericsType_C) :: numerics_c
type (sfx_phys_constants) :: phys_constants
->>>>>>> 9d99a415378a2907d460477f87825d027fae071e
numerics_c%maxiters_charnock = numerics%maxiters_charnock
phys_constants%Pr_m = Pr_m;
phys_constants%nu_air = nu_air;
phys_constants%g = g;
-<<<<<<< HEAD
call set_c_struct_sfx_sheba_noit_param_values(model_param)
call set_c_struct_sfx_surface_param_values(surface_param)
call set_meteo_vec_c(meteo, meteo_c)
@@ -213,21 +172,7 @@ contains
call c_sheba_noit_compute_flux(sfx_c_ptr, meteo_c_ptr, model_param, surface_param, numerics_c, phys_constants, n)
#else
do i = 1, n
-#ifdef SFX_FORCE_DEPRECATED_sheba_CODE
-#else
-=======
- call set_c_struct_sfx_sheba_param_values(model_param)
- call set_c_struct_sfx_surface_param_values(surface_param)
- call set_meteo_vec_c(meteo, meteo_c)
- call set_sfx_vec_c(sfx, sfx_c)
-
- meteo_c_ptr = C_LOC(meteo_c)
- sfx_c_ptr = C_LOC(sfx_c)
- call c_sheba_compute_flux(sfx_c_ptr, meteo_c_ptr, model_param, surface_param, numerics_c, phys_constants, n)
-#else
- do i = 1, n
->>>>>>> 9d99a415378a2907d460477f87825d027fae071e
meteo_cell = meteoDataType(&
h = meteo%h(i), &
U = meteo%U(i), dT = meteo%dT(i), Tsemi = meteo%Tsemi(i), dQ = meteo%dQ(i), &
@@ -236,15 +181,9 @@ contains
call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
call push_sfx_data(sfx, sfx_cell, i)
-<<<<<<< HEAD
-#endif
end do
#endif
-=======
- end do
-#endif
->>>>>>> 9d99a415378a2907d460477f87825d027fae071e
end subroutine get_surface_fluxes_vec
! --------------------------------------------------------------------------------
@@ -262,10 +201,6 @@ contains
type (meteoDataType), intent(in) :: meteo
type (numericsType), intent(in) :: numerics
! ----------------------------------------------------------------------------
-<<<<<<< HEAD
-=======
-
->>>>>>> 9d99a415378a2907d460477f87825d027fae071e
! --- meteo derived datatype name shadowing
! ----------------------------------------------------------------------------
real :: h !< constant flux layer height [m]
@@ -288,7 +223,6 @@ contains
real zeta !< = z/L [n/d]
real Rib !< bulk Richardson number
-<<<<<<< HEAD
real zeta_conv_lim !< z/L critical value for matching free convection limit [n/d]
real Rib_conv_lim !< Ri-bulk critical value for matching free convection limit [n/d]
@@ -297,30 +231,22 @@ contains
real psi_m, psi_h !< universal functions (momentum) & (heat) [n/d]
real psi0_m, psi0_h !< universal functions (momentum) & (heat) [n/d]
-=======
+
real Udyn, Tdyn, Qdyn !< dynamic scales
->>>>>>> 9d99a415378a2907d460477f87825d027fae071e
real phi_m, phi_h !< stability functions (momentum) & (heat) [n/d]
real Km !< eddy viscosity coeff. at h [m^2/s]
real Pr_t_inv !< invese Prandt number [n/d]
real Cm, Ct !< transfer coeff. for (momentum) & (heat) [n/d]
-<<<<<<< HEAD
-
- real Udyn, Tdyn
integer surface_type !< surface type = (ocean || land)
real fval !< just a shortcut for partial calculations
real :: C1,A1,A2,lne,lnet,Ribl
-=======
- integer surface_type !< surface type = (ocean || land)
-
->>>>>>> 9d99a415378a2907d460477f87825d027fae071e
#ifdef SFX_CHECK_NAN
real NaN
@@ -379,7 +305,6 @@ contains
! --- define Ri-bulk
Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / U**2
-<<<<<<< HEAD
! --- define free convection transition zeta = z/L value
call get_convection_lim(zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim, &
h0_m, h0_t, B)
@@ -458,54 +383,20 @@ contains
Cm = kappa / psi_m
Ct = kappa / psi_h
end if
-
-=======
- ! --- get the fluxes
- ! ----------------------------------------------------------------------------
- if(Rib > 0)then
- call get_dynamic_scales_noniterative(Udyn, Tdyn, Qdyn, zeta, &
- U, dT, dQ, h, z0_m, z0_t, Rib)
- else
- call get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, &
- U, Tsemi, dT, dQ, h, z0_m, z0_t, (g / Tsemi), 10)
- end if
- ! ----------------------------------------------------------------------------
-
- call get_phi(phi_m, phi_h, zeta)
- ! ----------------------------------------------------------------------------
-
- ! --- define transfer coeff. (momentum) & (heat)
- Cm = 0.0
- if (U > 0.0) then
- Cm = Udyn / U
- end if
- Ct = 0.0
- if (abs(dT) > 0.0) then
- Ct = Tdyn / dT
- end if
-
->>>>>>> 9d99a415378a2907d460477f87825d027fae071e
! --- define eddy viscosity & inverse Prandtl number
Km = kappa * Cm * U * h / phi_m
Pr_t_inv = phi_m / phi_h
! --- setting output
sfx = sfxDataType(zeta = zeta, Rib = Rib, &
-<<<<<<< HEAD
Re = Re, B = B, z0_m = z0_m, z0_t = z0_t, &
Rib_conv_lim = Rib_conv_lim, &
Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv)
-=======
- Re = Re, B = B, z0_m = z0_m, z0_t = z0_t, &
- Rib_conv_lim = 0.0, &
- Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv)
->>>>>>> 9d99a415378a2907d460477f87825d027fae071e
end subroutine get_surface_fluxes
! --------------------------------------------------------------------------------
-<<<<<<< HEAD
! convection universal functions shortcuts
! --------------------------------------------------------------------------------
function f_m_conv(zeta)
@@ -728,224 +619,6 @@ contains
! --- bulk Richardson number
Rib_lim = zeta_lim * psi_h / (psi_m * psi_m)
-=======
- !< @brief get dynamic scales
- ! --------------------------------------------------------------------------------
- subroutine get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, &
- U, Tsemi, dT, dQ, z, z0_m, z0_t, beta, maxiters)
- ! ----------------------------------------------------------------------------
- real, intent(out) :: Udyn, Tdyn, Qdyn !< dynamic scales
- real, intent(out) :: zeta !< = z/L
-
- real, intent(in) :: U !< abs(wind speed) at z
- real, intent(in) :: Tsemi !< semi-sum of temperature at z and at surface
- real, intent(in) :: dT, dQ !< temperature & humidity difference between z and at surface
- real, intent(in) :: z !< constant flux layer height
- real, intent(in) :: z0_m, z0_t !< roughness parameters
- real, intent(in) :: beta !< buoyancy parameter
-
- integer, intent(in) :: maxiters !< maximum number of iterations
- ! ----------------------------------------------------------------------------
-
- ! --- local variables
- real, parameter :: gamma = 0.61
-
- real :: psi_m, psi_h
- real :: psi0_m, psi0_h
- real :: Linv
- integer :: i
- ! ----------------------------------------------------------------------------
-
-
- 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
-
- do i = 1, maxiters
-
- call get_psi(psi_m, psi_h, zeta)
- call get_psi_mh(psi0_m, psi0_h, z0_m * Linv, z0_t * Linv)
-
- 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) exit
-
- Linv = kappa * beta * (Tdyn + gamma * Qdyn * Tsemi) / (Udyn * Udyn)
- zeta = z * Linv
- end do
-
- end subroutine get_dynamic_scales
-
-
- subroutine get_dynamic_scales_noniterative(Udyn, Tdyn, Qdyn, zeta, &
- U, dT, dQ, z, z0_m, z0_t, Rib)
- ! ----------------------------------------------------------------------------
- real, parameter :: gamma = 2.91, zeta_a = 3.6
-
- real, intent(out) :: Udyn, Tdyn, Qdyn !< dynamic scales
- real, intent(out) :: zeta !< = z/L
-
- real, intent(in) :: U !< abs(wind speed) at z
- real, intent(in) :: dT, dQ !< temperature & humidity difference between z and at surface
- real, intent(in) :: z !< constant flux layer height
- real, intent(in) :: z0_m, z0_t !< roughness parameters
- real, intent(in) :: Rib !< bulk Richardson number
-
- ! ----------------------------------------------------------------------------
-
- ! --- local variables
- real :: psi_m, psi_h
- real :: psi0_m, psi0_h
- real :: C1,A1,A2,lne,lnet,Ribl
- ! ----------------------------------------------------------------------------
-
- Ribl = (Rib*Pr_t_0_inv) * (1 - z0_t / z) / ((1 - z0_m / z)**2)
-
- call get_psi(psi_m, psi_h, zeta_a)
- call get_psi_mh(psi0_m, psi0_h, zeta_a * z0_m / z, zeta_a * z0_t / z)
-
- lne = log(z/z0_m)
- lnet = log(z/z0_t)
- C1 = (lne**2)/lnet
- A1 = ((lne - psi_m + psi0_m)**(2*(gamma-1))) &
-& / ((zeta_a**(gamma-1))*((lnet-(psi_h-psi0_h)*Pr_t_0_inv)**(gamma-1)))
- A2 = ((lne - psi_m + psi0_m)**2) / (lnet-(psi_h-psi0_h)*Pr_t_0_inv) - C1
-
- zeta = C1 * Ribl + A1 * A2 * (Ribl**gamma)
-
- call get_psi(psi_m, psi_h, zeta)
- call get_psi_mh(psi0_m, psi0_h, zeta * z0_m / z, zeta * z0_t /z)
-
- 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))
-
-
- end subroutine get_dynamic_scales_noniterative
- ! --------------------------------------------------------------------------------
-
- ! stability functions
- ! --------------------------------------------------------------------------------
- subroutine get_phi(phi_m, phi_h, zeta)
- !< @brief stability functions (momentum) & (heat): neutral case
- ! ----------------------------------------------------------------------------
- real, intent(out) :: phi_m, phi_h !< stability functions
-
- real, intent(in) :: zeta !< = z/L
- ! ----------------------------------------------------------------------------
-
-
- if (zeta >= 0.0) then
- phi_m = 1.0 + (a_m * zeta * (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 = (1.0 - alpha_m * zeta)**(-0.25)
- phi_h = (1.0 - alpha_h * zeta)**(-0.5)
- end if
-
- end subroutine
- ! --------------------------------------------------------------------------------
-
- ! universal functions
- ! --------------------------------------------------------------------------------
- subroutine get_psi(psi_m, psi_h, zeta)
- !< @brief universal functions (momentum) & (heat): neutral case
- ! ----------------------------------------------------------------------------
- real, intent(out) :: psi_m, psi_h !< universal functions
-
- real, intent(in) :: zeta !< = z/L
- ! ----------------------------------------------------------------------------
-
- ! --- local variables
- real :: x_m, x_h
- real :: q_m, q_h
- ! ----------------------------------------------------------------------------
-
-
- if (zeta >= 0.0) then
-
- q_m = ((1.0 - b_m) / b_m)**(1.0 / 3.0)
- q_h = sqrt(c_h * c_h - 4.0)
-
- x_m = (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 = (1.0 - alpha_m * zeta)**(0.25)
- x_h = (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))
- end if
-
- end subroutine
-
-
- subroutine get_psi_mh(psi_m, psi_h, zeta_m, zeta_h)
- !< @brief universal functions (momentum) & (heat): neutral case
- ! ----------------------------------------------------------------------------
- real, intent(out) :: psi_m, psi_h !< universal functions
-
- real, intent(in) :: zeta_m, zeta_h !< = z/L
- ! ----------------------------------------------------------------------------
-
- ! --- local variables
- real :: x_m, x_h
- real :: q_m, q_h
- ! ----------------------------------------------------------------------------
-
-
- if (zeta_m >= 0.0) then
- q_m = ((1.0 - b_m) / b_m)**(1.0 / 3.0)
- x_m = (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 = (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)
- end if
-
- if (zeta_h >= 0.0) then
- 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 = (1.0 - alpha_h * zeta_h)**(0.25)
- psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h))
- end if
->>>>>>> 9d99a415378a2907d460477f87825d027fae071e
end subroutine
! --------------------------------------------------------------------------------