Newer
Older

Evgeny Mortikov
committed
module sfx_log

Evgeny Mortikov
committed
! modules used
! --------------------------------------------------------------------------------
#ifdef SFX_CHECK_NAN
use sfx_common
#endif
use sfx_data

Evgeny Mortikov
committed
use sfx_surface

Evgeny Mortikov
committed
use sfx_log_param
! --------------------------------------------------------------------------------
! directives list
! --------------------------------------------------------------------------------
implicit none
private
! --------------------------------------------------------------------------------
! public interface
! --------------------------------------------------------------------------------
public :: get_surface_fluxes
public :: get_surface_fluxes_vec

Виктория Суязова
committed
integer z0m_id
integer z0t_id

Evgeny Mortikov
committed
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
type, public :: numericsType
integer :: maxiters_charnock = 10 !< maximum (actual) number of iterations in charnock roughness

Evgeny Mortikov
committed
end type
! --------------------------------------------------------------------------------
contains
! --------------------------------------------------------------------------------
subroutine get_surface_fluxes_vec(sfx, meteo, numerics, n)
!< @brief surface flux calculation for array data
!< @details contains C/C++ & CUDA interface

Evgeny Mortikov
committed
! ----------------------------------------------------------------------------
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
! ----------------------------------------------------------------------------
do i = 1, n
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), z0_t = meteo%z0_t(i), &
depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))

Evgeny Mortikov
committed
call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
call push_sfx_data(sfx, sfx_cell, i)
end do
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

Evgeny Mortikov
committed
! ----------------------------------------------------------------------------
#ifdef SFX_CHECK_NAN
use ieee_arithmetic
#endif
type (sfxDataType), intent(out) :: sfx
type (meteoDataType), intent(in) :: meteo
type (numericsType), intent(in) :: numerics
! ----------------------------------------------------------------------------
! --- 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)

Виктория Суязова
committed
real :: depth
real :: lai

Evgeny Mortikov
committed
! ----------------------------------------------------------------------------
! --- 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]

Evgeny Mortikov
committed
real u_dyn0 !< dynamic velocity in neutral conditions [m/s]
real Re !< roughness Reynolds number = u_dyn0 * z0_m / nu [n/d]

Evgeny Mortikov
committed

Evgeny Mortikov
committed
real psi_m, psi_h !< universal functions (momentum) & (heat) [n/d]
real phi_m, phi_h !< stability functions (momentum) & (heat) [n/d]

Evgeny Mortikov
committed
real Km !< eddy viscosity coeff. at h [m^2/s]
real Pr_t_inv !< invese Prandt number [n/d]

Evgeny Mortikov
committed
real Cm, Ct !< transfer coeff. for (momentum) & (heat) [n/d]

Evgeny Mortikov
committed
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
#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

Виктория Суязова
committed
depth = meteo%depth
lai = meteo%lai
surface_type=meteo%surface_type
call get_dynamic_roughness_definition(surface_type, ocean_z0m_id, land_z0m_id, lake_z0m_id, snow_z0m_id, &
forest_z0m_id, usersf_z0m_id, ice_z0m_id, mosaic_z0m_id, z0m_id)

Виктория Суязова
committed
call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m1, z0m_id)

Виктория Суязова
committed
call get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t_id, lake_z0t_id, snow_z0t_id, &
forest_z0t_id, usersf_z0t_id, ice_z0t_id, mosaic_z0t_id, z0t_id)

Виктория Суязова
committed

Evgeny Mortikov
committed
Re = u_dyn0 * z0_m / nu_air

Виктория Суязова
committed
call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0_t1, z0t_id)

Виктория Суязова
committed
! --- define relative height
h0_m = h / z0_m

Evgeny Mortikov
committed
! --- define relative height [thermal]
h0_t = h / z0_t
! --- define Ri-bulk
Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / U**2
! --- get the fluxes
! ----------------------------------------------------------------------------
call get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B)
! ----------------------------------------------------------------------------

Evgeny Mortikov
committed
phi_m = 1.0
phi_h = 1.0 / Pr_t_0_inv
! --- define transfer coeff. (momentum) & (heat)
Cm = kappa / psi_m
Ct = kappa / psi_h
! --- 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 = 0.0, Rib = Rib, &

Evgeny Mortikov
committed
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)
end subroutine get_surface_fluxes
! --------------------------------------------------------------------------------
! universal functions
! --------------------------------------------------------------------------------
subroutine get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B)
!< @brief universal functions (momentum) & (heat): neutral case

Evgeny Mortikov
committed
! ----------------------------------------------------------------------------

Evgeny Mortikov
committed
real, intent(in) :: h0_m, h0_t !< = z/z0_m, z/z0_h
real, intent(in) :: B !< = log(z0_m / z0_h)

Evgeny Mortikov
committed
! ----------------------------------------------------------------------------
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

Evgeny Mortikov
committed
if (abs(B) < 1.0e-10) psi_h = psi_m / Pr_t_0_inv
end subroutine
! --------------------------------------------------------------------------------