diff --git a/srcF/sfx_sheba_noniterative.f90 b/srcF/sfx_sheba_noniterative.f90 deleted file mode 100644 index c4bbc50c1fdb7a5b5672964eb45eec7462d571da..0000000000000000000000000000000000000000 --- a/srcF/sfx_sheba_noniterative.f90 +++ /dev/null @@ -1,953 +0,0 @@ -#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 - ! -------------------------------------------------------------------------------- -#ifdef SFX_CHECK_NAN - use sfx_common -#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 -#endif - ! -------------------------------------------------------------------------------- - - ! directives list - ! -------------------------------------------------------------------------------- - implicit none - private - ! -------------------------------------------------------------------------------- - - ! public interface - ! -------------------------------------------------------------------------------- - 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 - real(C_FLOAT) :: Pr_t_inf_inv - - real(C_FLOAT) :: alpha_m - real(C_FLOAT) :: alpha_h - real(C_FLOAT) :: alpha_h_fix - real(C_FLOAT) :: Rib_max - real(C_FLOAT) :: gamma - real(C_FLOAT) :: zeta_a - 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 - use, intrinsic :: ISO_C_BINDING, ONLY: C_INT, C_PTR - Import :: sfx_sheba_noit_param_C, sfx_sheba_noit_numericsType_C - implicit none - integer(C_INT) :: grid_size - type(C_PTR), value :: sfx - type(C_PTR), value :: meteo - type(sfx_sheba_noit_param_C) :: model_param - type(sfx_surface_noit_param) :: surface_param - 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) - type (sfx_sheba_noit_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%Pr_t_inf_inv = Pr_t_inf_inv - - sfx_model_param%alpha_m = alpha_m - sfx_model_param%alpha_h = alpha_h - sfx_model_param%alpha_h_fix = alpha_h_fix - sfx_model_param%Rib_max = Rib_max - sfx_model_param%gamma = gamma - sfx_model_param%zeta_a = zeta_a - 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 - ! ---------------------------------------------------------------------------- - type (sfxDataVecType), intent(inout) :: sfx - - type (meteoDataVecType), intent(in) :: meteo - type (numericsType), intent(in) :: numerics - integer, intent(in) :: n - ! ---------------------------------------------------------------------------- - - ! --- local variables - type (meteoDataType) meteo_cell - type (sfxDataType) sfx_cell - integer i - ! ---------------------------------------------------------------------------- -#if defined(INCLUDE_CXX) - 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) - 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_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), & - z0_m = meteo%z0_m(i)) - - 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 - ! -------------------------------------------------------------------------------- - - ! -------------------------------------------------------------------------------- - subroutine get_surface_fluxes(sfx, meteo, numerics) - !< @brief surface flux calculation for single cell - !< @details contains C/C++ interface - ! ---------------------------------------------------------------------------- -#ifdef SFX_CHECK_NAN - use ieee_arithmetic -#endif - - type (sfxDataType), intent(out) :: sfx - - type (meteoDataType), intent(in) :: meteo - type (numericsType), intent(in) :: numerics - ! ---------------------------------------------------------------------------- -<<<<<<< HEAD -======= - ->>>>>>> 9d99a415378a2907d460477f87825d027fae071e - ! --- meteo derived datatype name shadowing - ! ---------------------------------------------------------------------------- - real :: h !< constant flux layer height [m] - real :: U !< abs(wind speed) at 'h' [m/s] - real :: dT !< difference between potential temperature at 'h' and at surface [K] - real :: Tsemi !< semi-sum of potential temperature at 'h' and at surface [K] - real :: dQ !< difference between humidity at 'h' and at surface [g/g] - real :: z0_m !< surface aerodynamic roughness (should be < 0 for water bodies surface) - ! ---------------------------------------------------------------------------- - - ! --- local variables - ! ---------------------------------------------------------------------------- - real z0_t !< thermal roughness [m] - real B !< = ln(z0_m / z0_t) [n/d] - real h0_m, h0_t !< = h / z0_m, h / z0_h [n/d] - - real u_dyn0 !< dynamic velocity in neutral conditions [m/s] - real Re !< roughness Reynolds number = u_dyn0 * z0_m / nu [n/d] - - 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] - - real f_m_conv_lim !< stability function (momentum) value in free convection limit [n/d] - real f_h_conv_lim !< stability function (heat) value in free convection limit [n/d] - - 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 -#endif - ! ---------------------------------------------------------------------------- - -#ifdef SFX_CHECK_NAN - ! --- checking if arguments are finite - if (.not.(is_finite(meteo%U).and.is_finite(meteo%Tsemi).and.is_finite(meteo%dT).and.is_finite(meteo%dQ) & - .and.is_finite(meteo%z0_m).and.is_finite(meteo%h))) then - - NaN = ieee_value(0.0, ieee_quiet_nan) ! setting NaN - sfx = sfxDataType(zeta = NaN, Rib = NaN, & - Re = NaN, B = NaN, z0_m = NaN, z0_t = NaN, & - Rib_conv_lim = NaN, & - Cm = NaN, Ct = NaN, Km = NaN, Pr_t_inv = NaN) - return - end if -#endif - - ! --- shadowing names for clarity - U = meteo%U - Tsemi = meteo%Tsemi - dT = meteo%dT - dQ = meteo%dQ - h = meteo%h - z0_m = meteo%z0_m - - ! --- define surface type - if (z0_m < 0.0) then - surface_type = surface_ocean - else - surface_type = surface_land - end if - - if (surface_type == surface_ocean) then - ! --- define surface roughness [momentum] & dynamic velocity in neutral conditions - call get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_charnock) - ! --- define relative height - h0_m = h / z0_m - endif - if (surface_type == surface_land) then - ! --- define relative height - h0_m = h / z0_m - ! --- define dynamic velocity in neutral conditions - u_dyn0 = U * kappa / log(h0_m) - end if - - ! --- define thermal roughness & B = log(z0_m / z0_h) - Re = u_dyn0 * z0_m / nu_air - call get_thermal_roughness(z0_t, B, z0_m, Re, surface_type) - - ! --- define relative height [thermal] - h0_t = h / z0_t - - ! --- 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) - - ! --- get the fluxes - ! ---------------------------------------------------------------------------- - if (Rib > 0.0) then - ! --- stable stratification block - - ! --- restrict bulk Ri value - ! *: note that this value is written to output - Rib = min(Rib, Rib_max) - - Ribl = (Rib*Pr_t_0_inv) * (1 - z0_t / h) / ((1 - z0_m / h)**2) - - call get_psi_stable(psi_m, psi_h, zeta_a, zeta_a) - call get_psi_stable(psi0_m, psi0_h, zeta_a * z0_m / h, zeta_a * z0_t / h) - - lne = log(h/z0_m) - lnet = log(h/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_stable(psi_m, psi_h, zeta, zeta) - call get_psi_stable(psi0_m, psi0_h, zeta * z0_m / h, zeta * z0_t / h) - - 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) - - Udyn = kappa * U / (log(h / z0_m) - (psi_m - psi0_m)) - Tdyn = kappa * dT * Pr_t_0_inv / (log(h / z0_t) - (psi_h - psi0_h)) - - else if (Rib < Rib_conv_lim) then - ! --- strong instability block - - call get_psi_convection(psi_m, psi_h, zeta, Rib, & - zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, h0_m, h0_t, B, numerics%maxiters_convection) - - fval = (zeta_conv_lim / zeta)**(1.0/3.0) - phi_m = fval / f_m_conv_lim - phi_h = fval / (Pr_t_0_inv * f_h_conv_lim) - - else if (Rib > -0.001) then - ! --- nearly neutral [-0.001, 0] block - - call get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B) - - zeta = 0.0 - phi_m = 1.0 - phi_h = 1.0 / Pr_t_0_inv - else - ! --- weak & semistrong instability block - - call get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, numerics%maxiters_convection) - - phi_m = (1.0 - alpha_m * zeta)**(-0.25) - phi_h = 1.0 / (Pr_t_0_inv * sqrt(1.0 - alpha_h_fix * zeta)) - end if - ! ---------------------------------------------------------------------------- - - ! --- define transfer coeff. (momentum) & (heat) - if(Rib > 0)then - 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 - else - 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) - ! ---------------------------------------------------------------------------- - real :: f_m_conv - real, intent(in) :: zeta - ! ---------------------------------------------------------------------------- - - f_m_conv = (1.0 - alpha_m * zeta)**0.25 - - end function f_m_conv - - function f_h_conv(zeta) - ! ---------------------------------------------------------------------------- - real :: f_h_conv - real, intent(in) :: zeta - ! ---------------------------------------------------------------------------- - - f_h_conv = (1.0 - alpha_h * zeta)**0.5 - - end function f_h_conv - ! -------------------------------------------------------------------------------- - - ! universal functions - ! -------------------------------------------------------------------------------- - subroutine get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B) - !< @brief universal functions (momentum) & (heat): neutral case - ! ---------------------------------------------------------------------------- - real, intent(out) :: psi_m, psi_h !< universal functions - - real, intent(in) :: h0_m, h0_t !< = z/z0_m, z/z0_h - real, intent(in) :: B !< = log(z0_m / z0_h) - ! ---------------------------------------------------------------------------- - - psi_m = log(h0_m) - psi_h = log(h0_t) / Pr_t_0_inv - !*: this looks redundant z0_t = z0_m in case |B| ~ 0 - if (abs(B) < 1.0e-10) psi_h = psi_m / Pr_t_0_inv - - end subroutine - - subroutine get_psi_stable(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 - ! ---------------------------------------------------------------------------- - - - 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)))) - 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))) - - end subroutine - - - subroutine get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, maxiters) - !< @brief universal functions (momentum) & (heat): semi-strong convection case - ! ---------------------------------------------------------------------------- - real, intent(out) :: psi_m, psi_h !< universal functions [n/d] - real, intent(out) :: zeta !< = z/L [n/d] - - real, intent(in) :: Rib !< bulk Richardson number [n/d] - real, intent(in) :: h0_m, h0_t !< = z/z0_m, z/z0_h [n/d] - real, intent(in) :: B !< = log(z0_m / z0_h) [n/d] - integer, intent(in) :: maxiters !< maximum number of iterations - - ! --- local variables - real :: zeta0_m, zeta0_h - real :: f0_m, f0_h - real :: f_m, f_h - - integer :: i - ! ---------------------------------------------------------------------------- - - psi_m = log(h0_m) - psi_h = log(h0_t) - if (abs(B) < 1.0e-10) psi_h = psi_m - - zeta = Rib * Pr_t_0_inv * psi_m**2 / psi_h - - do i = 1, maxiters + 1 - zeta0_m = zeta / h0_m - zeta0_h = zeta / h0_t - if (abs(B) < 1.0e-10) zeta0_h = zeta0_m - - f_m = (1.0 - alpha_m * zeta)**0.25e0 - f_h = sqrt(1.0 - alpha_h_fix * zeta) - - f0_m = (1.0 - alpha_m * zeta0_m)**0.25e0 - f0_h = sqrt(1.0 - alpha_h_fix * zeta0_h) - - f0_m = max(f0_m, 1.000001e0) - f0_h = max(f0_h, 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))) / Pr_t_0_inv - - ! *: don't update zeta = z/L at last iteration - if (i == maxiters + 1) exit - - zeta = Rib * psi_m**2 / psi_h - end do - - end subroutine - - subroutine get_psi_convection(psi_m, psi_h, zeta, Rib, & - zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, & - h0_m, h0_t, B, maxiters) - !< @brief universal functions (momentum) & (heat): fully convective case - ! ---------------------------------------------------------------------------- - real, intent(out) :: psi_m, psi_h !< universal functions [n/d] - real, intent(out) :: zeta !< = z/L [n/d] - - real, intent(in) :: Rib !< bulk Richardson number [n/d] - real, intent(in) :: h0_m, h0_t !< = z/z0_m, z/z0_h [n/d] - real, intent(in) :: B !< = log(z0_m / z0_h) [n/d] - integer, intent(in) :: maxiters !< maximum number of iterations - - real, intent(in) :: zeta_conv_lim !< convective limit zeta - real, intent(in) :: f_m_conv_lim, f_h_conv_lim !< universal function shortcuts in limiting case - ! ---------------------------------------------------------------------------- - - ! --- local variables - real :: zeta0_m, zeta0_h - real :: f0_m, f0_h - real :: p_m, p_h - real :: a_m, a_h - real :: c_lim, f - - integer :: i - ! ---------------------------------------------------------------------------- - - 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)) - - zeta = zeta_conv_lim - do i = 1, maxiters + 1 - zeta0_m = zeta / h0_m - zeta0_h = zeta / h0_t - if (abs(B) < 1.0e-10) zeta0_h = zeta0_m - - f0_m = (1.0 - alpha_m * zeta0_m)**0.25 - f0_h = sqrt(1.0 - 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)) - - c_lim = (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 - psi_h = (f / f_h_conv_lim + p_h + a_h) / Pr_t_0_inv - - ! *: don't update zeta = z/L at last iteration - if (i == maxiters + 1) exit - - zeta = Rib * psi_m**2 / psi_h - end do - - end subroutine - ! -------------------------------------------------------------------------------- - - ! convection limit definition - ! -------------------------------------------------------------------------------- - subroutine get_convection_lim(zeta_lim, Rib_lim, f_m_lim, f_h_lim, & - h0_m, h0_t, B) - ! ---------------------------------------------------------------------------- - real, intent(out) :: zeta_lim !< limiting value of z/L - real, intent(out) :: Rib_lim !< limiting value of Ri-bulk - real, intent(out) :: f_m_lim, f_h_lim !< limiting values of universal functions shortcuts - - real, intent(in) :: h0_m, h0_t !< = z/z0_m, z/z0_h [n/d] - real, intent(in) :: B !< = log(z0_m / z0_h) [n/d] - ! ---------------------------------------------------------------------------- - - ! --- local variables - real :: psi_m, psi_h - real :: f_m, f_h - real :: c - ! ---------------------------------------------------------------------------- - - ! --- define limiting value of zeta = z / L - c = (Pr_t_inf_inv / Pr_t_0_inv)**4 - zeta_lim = (2.0 * alpha_h - c * alpha_m - & - sqrt((c * alpha_m)**2 + 4.0 * c * alpha_h * (alpha_h - alpha_m))) / (2.0 * alpha_h**2) - - f_m_lim = f_m_conv(zeta_lim) - f_h_lim = f_h_conv(zeta_lim) - - ! --- universal functions - f_m = zeta_lim / h0_m - f_h = zeta_lim / h0_t - if (abs(B) < 1.0e-10) f_h = f_m - - f_m = (1.0 - alpha_m * f_m)**0.25 - f_h = sqrt(1.0 - alpha_h_fix * f_h) - - psi_m = 2.0 * (atan(f_m_lim) - atan(f_m)) + alog((f_m_lim - 1.0) * (f_m + 1.0)/((f_m_lim + 1.0) * (f_m - 1.0))) - psi_h = alog((f_h_lim - 1.0) * (f_h + 1.0)/((f_h_lim + 1.0) * (f_h - 1.0))) / Pr_t_0_inv - - ! --- 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 - ! -------------------------------------------------------------------------------- - -end module sfx_sheba_noniterative \ No newline at end of file