diff --git a/CMakeLists.txt b/CMakeLists.txt index 9ed5741ad1ed7d36301778922b5eab91f69451c9..954ec6ff6435adb0274c854df1cbea900b4593ae 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -88,6 +88,7 @@ set(SOURCES_F srcF/sfx_most.f90 srcF/sfx_most_param.f90 srcF/sfx_sheba.f90 + srcF/sfx_sheba_noniterative.f90 srcF/sfx_sheba_param.f90 srcF/sfx_sheba_noniterative.f90 srcF/sfx_sheba_noit_param.f90 diff --git a/srcF/sfx_sheba_noniterative.f90 b/srcF/sfx_sheba_noniterative.f90 index 19d01aa62da6109ffdf7831015a05fbd78be5489..c4bbc50c1fdb7a5b5672964eb45eec7462d571da 100644 --- a/srcF/sfx_sheba_noniterative.f90 +++ b/srcF/sfx_sheba_noniterative.f90 @@ -1,7 +1,11 @@ #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 ! -------------------------------------------------------------------------------- @@ -10,7 +14,12 @@ 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 @@ -27,16 +36,24 @@ 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 @@ -52,10 +69,27 @@ module sfx_sheba_noniterative 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 @@ -70,11 +104,28 @@ 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 + use, intrinsic :: ISO_C_BINDING, ONLY: C_INT, C_PTR + Import :: sfx_sheba_param_C, sfx_sheba_numericsType_C + implicit none + INTEGER(C_INT) :: grid_size + type(C_PTR), value :: sfx + type(C_PTR), value :: meteo + type(sfx_sheba_param_C) :: model_param + type(sfx_surface_param) :: surface_param + 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) @@ -92,6 +143,25 @@ 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 + sfx_model_param%kappa = kappa + sfx_model_param%Pr_t_0_inv = Pr_t_0_inv + + sfx_model_param%alpha_m = alpha_m + sfx_model_param%alpha_h = alpha_h + sfx_model_param%a_m = a_m + sfx_model_param%b_m = b_m + sfx_model_param%a_h = a_h + sfx_model_param%b_h = b_h + sfx_model_param%c_h = c_h + end subroutine set_c_struct_sfx_sheba_param_values +#endif + + ! -------------------------------------------------------------------------------- +>>>>>>> 9d99a415378a2907d460477f87825d027fae071e subroutine get_surface_fluxes_vec(sfx, meteo, numerics, n) !< @brief surface flux calculation for array data !< @details contains C/C++ & CUDA interface @@ -112,18 +182,27 @@ 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) @@ -136,6 +215,19 @@ contains 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), & @@ -144,10 +236,15 @@ 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 ! -------------------------------------------------------------------------------- @@ -165,6 +262,10 @@ 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] @@ -187,6 +288,7 @@ 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] @@ -195,12 +297,17 @@ 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 @@ -209,6 +316,11 @@ contains 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 @@ -267,6 +379,7 @@ 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) @@ -346,19 +459,53 @@ contains 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) @@ -581,6 +728,224 @@ 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 ! --------------------------------------------------------------------------------