Skip to content
Snippets Groups Projects
Commit 30a8d300 authored by Anna Shestakova's avatar Anna Shestakova
Browse files

noniterative sheba v.2

parent 214d9533
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
srcF/sfx_sheba_noniterative.f90
+ 4
− 331
View file @ 30a8d300
#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)
......@@ -459,53 +384,19 @@ 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)
......@@ -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
! --------------------------------------------------------------------------------
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment