diff --git a/compile.sh b/compile.sh
index 66a7c5a3b5bab8e4077652612b289600062d5e8c..22b260f68b85411e2cd6580d3157324367508425 100755
--- a/compile.sh
+++ b/compile.sh
@@ -14,5 +14,5 @@ gfortran -c -cpp -Wuninitialized srcF/sfx_esm_param.f90
gfortran -c -cpp -Wuninitialized srcF/sfx_esm.f90
gfortran -c -cpp -Wuninitialized srcF/sfx_main.f90
-gfortran -o sfx.exe sfx_main.o sfx_log.o sfx_log_param.o sfx_esm.o sfx_esm_param.o sfx_surface.o sfx_data.o sfx_io.o sfx_common.o sfx_phys_const.o
+gfortran -o sfx.exe sfx_main.o sfx_log.o sfx_log_param.o sfx_most.o sfx_most_param.o sfx_sheba.o sfx_sheba_param.o sfx_esm.o sfx_esm_param.o sfx_surface.o sfx_data.o sfx_io.o sfx_common.o sfx_phys_const.o
rm *.o *.mod
diff --git a/makefile b/makefile
index 5a0fd984f23733e423903686cbc14074c232f9d2..c14265ba2925727ed10e4ac59049fb4e6db02d46 100644
--- a/makefile
+++ b/makefile
@@ -10,7 +10,7 @@ ifeq ($(COMPILER),gnu)
FC = gfortran
endif
-OBJ_F90 = sfx_phys_const.o sfx_common.o sfx_io.o sfx_data.o sfx_surface.o sfx_log_param.o sfx_log.o sfx_esm_param.o sfx_esm.o sfx_main.o
+OBJ_F90 = sfx_phys_const.o sfx_common.o sfx_io.o sfx_data.o sfx_surface.o sfx_log_param.o sfx_log.o sfx_most_param.o sfx_most.o sfx_sheba_param.o sfx_sheba.o sfx_esm_param.o sfx_esm.o sfx_main.o
OBJ_F =
OBJ = $(OBJ_F90) $(OBJ_F)
diff --git a/srcF/sfx_main.f90 b/srcF/sfx_main.f90
index 7900d484373c3bea0d99d95adb63d86b2a0791a4..40d791265d0e255808181beee797af9dd7739447 100644
--- a/srcF/sfx_main.f90
+++ b/srcF/sfx_main.f90
@@ -13,6 +13,12 @@ program sfx_main
use sfx_log, only: &
get_surface_fluxes_vec_log => get_surface_fluxes_vec, &
numericsType_log => numericsType
+ use sfx_most, only: &
+ get_surface_fluxes_vec_most => get_surface_fluxes_vec, &
+ numericsType_most => numericsType
+ use sfx_sheba, only: &
+ get_surface_fluxes_vec_sheba => get_surface_fluxes_vec, &
+ numericsType_sheba => numericsType
! --------------------------------------------------------------------------------
! directives list
@@ -33,6 +39,8 @@ program sfx_main
character(len = 256) :: model_name
integer, parameter :: model_esm = 0 !< ESM model
integer, parameter :: model_log = 1 !< LOG simplified model
+ integer, parameter :: model_most = 2 !< MOST simplified model
+ integer, parameter :: model_sheba = 3 !< SHEBA simplified model
! input/output data
! --------------------------------------------------------------------------------
@@ -41,8 +49,10 @@ program sfx_main
type(sfxDataVecType) :: sfx !< surface fluxes (output)
- type(numericsType_esm) :: numerics_esm !< surface flux module (ESM) numerics parameters
- type(numericsType_log) :: numerics_log !< surface flux module (LOG) numerics parameters
+ type(numericsType_esm) :: numerics_esm !< surface flux module (ESM) numerics parameters
+ type(numericsType_log) :: numerics_log !< surface flux module (LOG) numerics parameters
+ type(numericsType_most) :: numerics_most !< surface flux module (MOST) numerics parameters
+ type(numericsType_sheba) :: numerics_sheba !< surface flux module (SHEBA) numerics parameters
integer :: num !< number of 'cells' in input
@@ -64,13 +74,15 @@ program sfx_main
character(len = 128), parameter :: arg_key_nmax = '--nmax'
character(len = 128), parameter :: arg_key_help = '--help'
- character(len = 128), parameter :: arg_key_esm = 'esm'
- character(len = 128), parameter :: arg_key_log = 'log'
+ character(len = 128), parameter :: arg_key_model_esm = 'esm'
+ character(len = 128), parameter :: arg_key_model_log = 'log'
+ character(len = 128), parameter :: arg_key_model_most = 'most'
+ character(len = 128), parameter :: arg_key_model_sheba = 'sheba'
- character(len = 128), parameter :: arg_key_mosaic = 'mosaic'
- character(len = 128), parameter :: arg_key_irgason = 'irgason'
- character(len = 128), parameter :: arg_key_sheba = 'sheba'
- character(len = 128), parameter :: arg_key_user = 'user'
+ character(len = 128), parameter :: arg_key_dataset_mosaic = 'mosaic'
+ character(len = 128), parameter :: arg_key_dataset_irgason = 'irgason'
+ character(len = 128), parameter :: arg_key_dataset_sheba = 'sheba'
+ character(len = 128), parameter :: arg_key_dataset_user = 'user'
integer :: is_output_set
integer :: nmax
@@ -97,9 +109,9 @@ program sfx_main
write(*, *) ' --help '
write(*, *) ' print usage options '
write(*, *) ' --model [key]'
- write(*, *) ' key = esm || log'
+ write(*, *) ' key = esm (default) || log || most || sheba'
write(*, *) ' --dataset [key]'
- write(*, *) ' key = mosaic || irgason || sheba || user [files]'
+ write(*, *) ' key = mosaic (default) || irgason || sheba || user [files]'
write(*, *) ' files = in-common-file in-file out-file'
write(*, *) ' --output [file]'
write(*, *) ' set output filename '
@@ -113,10 +125,14 @@ program sfx_main
stop
end if
call get_command_argument(i + 1, arg)
- if (trim(arg) == trim(arg_key_esm)) then
+ if (trim(arg) == trim(arg_key_model_esm)) then
model_id = model_esm
- else if (trim(arg) == trim(arg_key_log)) then
+ else if (trim(arg) == trim(arg_key_model_log)) then
model_id = model_log
+ else if (trim(arg) == trim(arg_key_model_most)) then
+ model_id = model_most
+ else if (trim(arg) == trim(arg_key_model_sheba)) then
+ model_id = model_sheba
else
write(*, *) ' FAILURE! > unknown model [key]: ', trim(arg)
stop
@@ -128,13 +144,13 @@ program sfx_main
stop
end if
call get_command_argument(i + 1, arg)
- if (trim(arg) == trim(arg_key_mosaic)) then
+ if (trim(arg) == trim(arg_key_dataset_mosaic)) then
dataset_id = dataset_MOSAiC
- else if (trim(arg) == trim(arg_key_irgason)) then
+ else if (trim(arg) == trim(arg_key_dataset_irgason)) then
dataset_id = dataset_IRGASON
- else if (trim(arg) == trim(arg_key_sheba)) then
+ else if (trim(arg) == trim(arg_key_dataset_sheba)) then
dataset_id = dataset_SHEBA
- else if (trim(arg) == trim(arg_key_user)) then
+ else if (trim(arg) == trim(arg_key_dataset_user)) then
dataset_id = dataset_USER
if (i + 4 > num_args) then
write(*, *) ' FAILURE! > incorrect arguments for [user] dataset'
@@ -180,6 +196,10 @@ program sfx_main
model_name = "ESM"
else if (model_id == model_log) then
model_name = "LOG"
+ else if (model_id == model_most) then
+ model_name = "MOST"
+ else if (model_id == model_sheba) then
+ model_name = "SHEBA"
else
write(*, *) ' FAILURE! > unknown model id: ', model_id
stop
@@ -270,6 +290,10 @@ program sfx_main
call get_surface_fluxes_vec_esm(sfx, meteo, numerics_esm, num)
else if (model_id == model_log) then
call get_surface_fluxes_vec_log(sfx, meteo, numerics_log, num)
+ else if (model_id == model_most) then
+ call get_surface_fluxes_vec_most(sfx, meteo, numerics_most, num)
+ else if (model_id == model_sheba) then
+ call get_surface_fluxes_vec_sheba(sfx, meteo, numerics_sheba, num)
end if
diff --git a/srcF/sfx_most.f90 b/srcF/sfx_most.f90
new file mode 100644
index 0000000000000000000000000000000000000000..799cf43d671596cd1de134899b53282c23ee2124
--- /dev/null
+++ b/srcF/sfx_most.f90
@@ -0,0 +1,347 @@
+#include "../includeF/sfx_def.fi"
+
+module sfx_most
+ !< @brief MOST surface flux module
+
+ ! modules used
+ ! --------------------------------------------------------------------------------
+#ifdef SFX_CHECK_NAN
+ use sfx_common
+#endif
+ use sfx_data
+ use sfx_surface
+ use sfx_most_param
+ ! --------------------------------------------------------------------------------
+
+ ! directives list
+ ! --------------------------------------------------------------------------------
+ implicit none
+ private
+ ! --------------------------------------------------------------------------------
+
+ ! public interface
+ ! --------------------------------------------------------------------------------
+ public :: get_surface_fluxes
+ public :: get_surface_fluxes_vec
+ ! --------------------------------------------------------------------------------
+
+ ! --------------------------------------------------------------------------------
+ type, public :: numericsType
+ integer :: maxiters_charnock = 10 !< maximum (actual) number of iterations in charnock roughness
+ 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
+ ! ----------------------------------------------------------------------------
+ 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))
+
+ 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
+ ! ----------------------------------------------------------------------------
+#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)
+ ! ----------------------------------------------------------------------------
+
+ ! --- 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
+
+ real Udyn, Tdyn, Qdyn !< dynamic scales
+
+ 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]
+
+ integer surface_type !< surface type = (ocean || land)
+
+
+#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
+
+ ! --- get the fluxes
+ ! ----------------------------------------------------------------------------
+ call get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, &
+ U, Tsemi, dT, dQ, h, z0_m, z0_t, (g / Tsemi), 10)
+ ! ----------------------------------------------------------------------------
+
+ 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
+
+ ! --- 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, &
+ 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
+ ! --------------------------------------------------------------------------------
+
+ !< @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
+ ! --------------------------------------------------------------------------------
+
+ ! 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 + beta_m * zeta
+ phi_h = 1.0 + beta_h * 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
+ ! ----------------------------------------------------------------------------
+
+
+ if (zeta >= 0.0) then
+ psi_m = -beta_m * zeta;
+ psi_h = -beta_h * zeta;
+ 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
+ ! ----------------------------------------------------------------------------
+
+
+ if (zeta_m >= 0.0) then
+ psi_m = -beta_m * zeta_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
+ psi_h = -beta_h * zeta_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
+
+ end subroutine
+ ! --------------------------------------------------------------------------------
+
+end module sfx_most
\ No newline at end of file
diff --git a/srcF/sfx_most_param.f90 b/srcF/sfx_most_param.f90
new file mode 100644
index 0000000000000000000000000000000000000000..51b34fd6610b2ea551017e8ae36e29517ac1ae1d
--- /dev/null
+++ b/srcF/sfx_most_param.f90
@@ -0,0 +1,30 @@
+module sfx_most_param
+ !< @brief MOST BD71 surface flux model parameters
+ !< @details all in SI units
+
+ ! modules used
+ ! --------------------------------------------------------------------------------
+ use sfx_phys_const
+ ! --------------------------------------------------------------------------------
+
+ ! directives list
+ ! --------------------------------------------------------------------------------
+ implicit none
+ ! --------------------------------------------------------------------------------
+
+
+ !< von Karman constant [n/d]
+ real, parameter :: kappa = 0.40
+ !< inverse Prandtl (turbulent) number in neutral conditions [n/d]
+ real, parameter :: Pr_t_0_inv = 1.15
+
+
+ !< stability function coeff. (unstable)
+ real, parameter :: alpha_m = 16.0
+ real, parameter :: alpha_h = 16.0
+
+ !< stability function coeff. (stable)
+ real, parameter :: beta_m = 4.7
+ real, parameter :: beta_h = beta_m
+
+end module sfx_most_param
diff --git a/srcF/sfx_sheba.f90 b/srcF/sfx_sheba.f90
new file mode 100644
index 0000000000000000000000000000000000000000..26231a763c69ec0c2f216c92278ce71268a12c71
--- /dev/null
+++ b/srcF/sfx_sheba.f90
@@ -0,0 +1,379 @@
+#include "../includeF/sfx_def.fi"
+
+module sfx_sheba
+ !< @brief SHEBA surface flux module
+
+ ! modules used
+ ! --------------------------------------------------------------------------------
+#ifdef SFX_CHECK_NAN
+ use sfx_common
+#endif
+ use sfx_data
+ use sfx_surface
+ use sfx_sheba_param
+ ! --------------------------------------------------------------------------------
+
+ ! directives list
+ ! --------------------------------------------------------------------------------
+ implicit none
+ private
+ ! --------------------------------------------------------------------------------
+
+ ! public interface
+ ! --------------------------------------------------------------------------------
+ public :: get_surface_fluxes
+ public :: get_surface_fluxes_vec
+ ! --------------------------------------------------------------------------------
+
+ ! --------------------------------------------------------------------------------
+ type, public :: numericsType
+ integer :: maxiters_charnock = 10 !< maximum (actual) number of iterations in charnock roughness
+ 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
+ ! ----------------------------------------------------------------------------
+ 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))
+
+ 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
+ ! ----------------------------------------------------------------------------
+#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)
+ ! ----------------------------------------------------------------------------
+
+ ! --- 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
+
+ real Udyn, Tdyn, Qdyn !< dynamic scales
+
+ 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]
+
+ integer surface_type !< surface type = (ocean || land)
+
+
+#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
+
+ ! --- get the fluxes
+ ! ----------------------------------------------------------------------------
+ call get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, &
+ U, Tsemi, dT, dQ, h, z0_m, z0_t, (g / Tsemi), 10)
+ ! ----------------------------------------------------------------------------
+
+ 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
+
+ ! --- 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, &
+ 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
+ ! --------------------------------------------------------------------------------
+
+ !< @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
+ ! --------------------------------------------------------------------------------
+
+ ! 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
+
+ end subroutine
+ ! --------------------------------------------------------------------------------
+
+end module sfx_sheba
\ No newline at end of file
diff --git a/srcF/sfx_sheba_param.f90 b/srcF/sfx_sheba_param.f90
new file mode 100644
index 0000000000000000000000000000000000000000..2b0e0411225bd08739fcb44f429f7d5f36e2fa3b
--- /dev/null
+++ b/srcF/sfx_sheba_param.f90
@@ -0,0 +1,32 @@
+module sfx_sheba_param
+ !< @brief SHEBA surface flux model parameters
+ !< @details all in SI units
+
+ ! modules used
+ ! --------------------------------------------------------------------------------
+ use sfx_phys_const
+ ! --------------------------------------------------------------------------------
+
+ ! directives list
+ ! --------------------------------------------------------------------------------
+ implicit none
+ ! --------------------------------------------------------------------------------
+
+
+ !< von Karman constant [n/d]
+ real, parameter :: kappa = 0.40
+ !< inverse Prandtl (turbulent) number in neutral conditions [n/d]
+ real, parameter :: Pr_t_0_inv = 1.15
+
+ !< stability function coeff. (unstable)
+ real, parameter :: alpha_m = 16.0
+ real, parameter :: alpha_h = 16.0
+
+ !< stability function coeff. (stable)
+ real, parameter :: a_m = 5.0
+ real, parameter :: b_m = a_m / 6.5
+ real, parameter :: a_h = 5.0
+ real, parameter :: b_h = 5.0
+ real, parameter :: c_h = 3.0
+
+end module sfx_sheba_param